DeepOPV-V.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "# coding: utf-8\n",
    "# author: Wanjun  HUANG\n",
    "# data: July 4th, 2021 \n",
    "import torch\n",
    "import numpy as np\n",
    "import scipy.io\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import math \n",
    "import torch.utils.data as Data\n",
    "from torch.autograd import Function\n",
    "import gc\n",
    "from scipy import sparse\n",
    "\n",
    "global MAXMIN_V, MAXMIN_PQg,BRANFT,BUS_SLACK\n",
    "global Ybus,baseMVA\n",
    "global Real_Ptrain, Real_Qtrain, Real_Vmtrain, Real_Vatrain, Real_PQdtrain\n",
    "global scale_vm, VmUb, VmLb, bus_slack, bus_PQg\n",
    "global flagVm,flagVa,DELTA,flag_hisv,Nbus,Ntest\n",
    "\n",
    "## whether there is GPU\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "print(device)\n",
    "print(\"Let's use\", torch.cuda.device_count(), \"GPUs!\")\n",
    "\n",
    "## load system name\n",
    "REPEAT = 1 #for speedup test: number of repeated computation \n",
    "s_epoch = 800 # minimum epoch for .pth model saving\n",
    "p_epoch = 10 # print loss for each \"p_epoch\" epoch\n",
    "model_version = 1 # version of model\n",
    "Nbus = 300 # number of buses\n",
    "sys_R = 2 # test case name\n",
    "flag_test = 1 # 0-train model; 1-test well-trained model\n",
    "flag_hisv = 1 # 1-use historical V to calculate dV;0-use predicted V to calculate dV\n",
    "EpochVm = 2 # maximum epoch of Vm\n",
    "EpochVa = 2 # maximum epoch of Va\n",
    "batch_size_training = 100 # mini-batch size for training\n",
    "batch_size_test = 1 # mini-batch size for test\n",
    "\n",
    "## hyperparameters\n",
    "Lrm = 1e-3 # learning rate for Vm\n",
    "Lra = 1e-3 # learning rate for Va\n",
    "k_dV = 1 # coefficient for dVa & dVm in post-processing\n",
    "DELTA = 1e-4 # threshold of violation\n",
    "scale_vm = torch.tensor([10]).float() # scaling of output Vm\n",
    "scale_va = torch.tensor([10]).float() # scaling of output Va\n",
    "\n",
    "# hidden layers for voltage magnitude (Vm) prediction\n",
    "if Nbus == 300:\n",
    "    khidden_Vm = np.array([8,6,4,2], dtype=int)\n",
    "    khidden_Va = np.array([8,6,4,2], dtype=int)\n",
    "    Neach = 12\n",
    "else:\n",
    "    khidden_Vm = np.array([8,4,2], dtype=int)\n",
    "    khidden_Va = np.array([8,4,2], dtype=int) \n",
    "    Neach = 8000\n",
    "    \n",
    "Ntrain = int(4*Neach) \n",
    "Nsample = int(5*Neach)\n",
    "Ntest = int(Neach)\n",
    "\n",
    "#name of hidden layers for result saving\n",
    "nmLm = 'Lm'\n",
    "for i in range(khidden_Vm.shape[0]):\n",
    "    nmLm = nmLm + str(khidden_Vm[i])\n",
    "    \n",
    "Lm = khidden_Vm.shape[0] # number of hidden layers\n",
    "\n",
    "# hidden layers for voltage angles (Va) prediction\n",
    "nmLa = 'La'\n",
    "for i in range(khidden_Va.shape[0]):\n",
    "    nmLa = nmLa + str(khidden_Va[i])    \n",
    "\n",
    "La = khidden_Va.shape[0]  # number of hidden layers  \n",
    "\n",
    "# results name\n",
    "PATHVm = './modelvm'+str(Nbus)+'r'+str(sys_R)+'N'+str(model_version)+nmLm+'E'+str(EpochVm)+'.pth'\n",
    "PATHVa = './modelva'+str(Nbus)+'r'+str(sys_R)+'N'+str(model_version)+nmLa+'E'+str(EpochVa)+'.pth'\n",
    "PATHVms = './modelvm'+str(Nbus)+'r'+str(sys_R)+'N'+str(model_version)+nmLm\n",
    "PATHVas = './modelva'+str(Nbus)+'r'+str(sys_R)+'N'+str(model_version)+nmLa\n",
    "resultnm = './res_'+str(Nbus)+'r'+str(sys_R)+'M'+str(model_version)+'H'+str(flag_hisv)+ 'NT'+str(Ntrain) \\\n",
    "+'B'+str(batch_size_training) +'Em'+str(EpochVm)+'Ea'+str(EpochVa)+nmLm+nmLa+'rp'+str(REPEAT)+'.mat'\n",
    "\n",
    "# load data case 300 large variance\n",
    "mat = scipy.io.loadmat('./data/XY_case300real.mat')\n",
    "load_idx = np.squeeze(mat['load_idx']).astype(int) - 1 \n",
    "matpara = scipy.io.loadmat('./data/pglib_opf_case'+str(Nbus)+'_ieeer'+str(sys_R)+'_para.mat')\n",
    "\n",
    "# power system parameters\n",
    "RPd0 = mat['RPd'] # Pd\n",
    "RQd0 = mat['RQd'] # Qd\n",
    "RPg = mat['RPg'] # Pg\n",
    "RQg = mat['RQg'] # Qg\n",
    "Ybus = matpara['Ybus']\n",
    "Yf = matpara['Yf']\n",
    "Yt = matpara['Yt']\n",
    "bus = matpara['bus']\n",
    "gen = matpara['gen']\n",
    "gencost = matpara['gencost']\n",
    "branch = matpara['branch']\n",
    "baseMVA = matpara['baseMVA']\n",
    "bus_slack = np.where(bus[:, 1] == 3)\n",
    "bus_slack = np.squeeze(bus_slack)\n",
    "BUS_SLACK = torch.from_numpy(bus_slack).long()\n",
    "print('bus_slack', bus_slack)\n",
    "\n",
    "# numpy array to spase matrix\n",
    "Ybus = sparse.csr_matrix(Ybus)\n",
    "Yf = sparse.csr_matrix(Yf) \n",
    "Yt = sparse.csr_matrix(Yt) \n",
    "\n",
    "# output data\n",
    "YVa = mat['RVa']*math.pi/180\n",
    "# do not need to predict voltaeg angles of slack bus \n",
    "YVa = np.delete(YVa, bus_slack, axis=1) \n",
    "YVm = mat['RVm']\n",
    "VmLb = mat['VmLb'] # lower bound of Vm\n",
    "VmUb = mat['VmUb'] # upper bound of Vm\n",
    "kvm = (YVm - VmLb)/(VmUb - VmLb) # scaled Vm\n",
    "\n",
    "# output data\n",
    "yva = YVa\n",
    "yvm = kvm\n",
    "print('yvm',yvm.shape,'yva',yva.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input data: only contain non-zeros loads\n",
    "idx_Pd = np.squeeze(np.where(np.abs(RPd0[0,:])> 0), axis=0)\n",
    "idx_Qd = np.squeeze(np.where(np.abs(RQd0[0,:])> 0), axis=0)\n",
    "bus_Pd = np.squeeze(load_idx[idx_Pd]).astype(int) - 1\n",
    "bus_Qd = np.squeeze(load_idx[idx_Qd]).astype(int) - 1\n",
    "# input data\n",
    "x = np.concatenate((RPd0[:, idx_Pd], RQd0[:, idx_Qd]), axis = 1)/baseMVA \n",
    "print('x', x.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pg Qg data: only contain generators that are turned on\n",
    "idxPg = np.squeeze(np.where(gen[:, 3] > 0), axis=0)\n",
    "idxQg = np.squeeze(np.where(gen[:, 1] > 0), axis=0)\n",
    "\n",
    "# Pd Qd of samples\n",
    "RPd = np.zeros((Nsample,Nbus))\n",
    "RQd = np.zeros((Nsample,Nbus))\n",
    "RPd[:, load_idx] = RPd0[0:Nsample]\n",
    "RQd[:, load_idx] = RQd0[0:Nsample]\n",
    "\n",
    "del RPd0, RQd0\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bus_Pg = gen[idxPg, 0].astype(int) - 1\n",
    "bus_Qg = gen[idxQg, 0].astype(int) - 1\n",
    "bus_PQg = np.concatenate((bus_Pg, bus_Qg + Nbus), axis = 0)\n",
    "MAXMIN_Pg = gen[idxPg, 3:5]/baseMVA\n",
    "MAXMIN_Qg = gen[idxQg, 1:3]/baseMVA\n",
    "MAXMIN_PQg = np.concatenate((MAXMIN_Pg, MAXMIN_Qg), axis = 0) #[Npg+NQg, 2]\n",
    "\n",
    "# numpy to tensor and output data scaling\n",
    "MAXMIN_PQg = torch.from_numpy(MAXMIN_Pg).float()\n",
    "bus = torch.from_numpy(bus).float()\n",
    "bus_PQg = torch.from_numpy(bus_PQg)\n",
    "VmLb = torch.from_numpy(VmLb).float()\n",
    "VmUb = torch.from_numpy(VmUb).float()\n",
    "yvm = torch.from_numpy(yvm).float()\n",
    "yva = torch.from_numpy(yva).float()\n",
    "yvms = yvm*scale_vm\n",
    "yvas = yva*scale_va\n",
    "# constraints violation: Vm \n",
    "MAXMIN_V = (bus[:, 2:4] - VmLb)*scale_vm/(VmUb - VmLb)\n",
    "MAXMIN_Vpu = bus[0,2:4]\n",
    "\n",
    "# branch angle incidence matrix\n",
    "BRANFT = torch.from_numpy(branch[:, 0:2] - 1).long()\n",
    "\n",
    "# loss function\n",
    "criterion = nn.MSELoss()\n",
    "\n",
    "# convert training data to tensor\n",
    "x_tensor = torch.from_numpy(x).float()\n",
    "yvm_tensor = yvms.float()\n",
    "yva_tensor = yvas.float()\n",
    "xtrain = x_tensor[0: Ntrain]\n",
    "yvmtrain = yvm_tensor[0: Ntrain]\n",
    "yvatrain = yva_tensor[0: Ntrain]\n",
    "\n",
    "# batch data\n",
    "training_dataset_vm = Data.TensorDataset(xtrain, yvmtrain)\n",
    "training_loader_vm = Data.DataLoader(\n",
    "        dataset=training_dataset_vm,\n",
    "        batch_size=batch_size_training,\n",
    "        shuffle=False,\n",
    "    )             \n",
    "\n",
    "training_dataset_va = Data.TensorDataset(xtrain, yvatrain)\n",
    "training_loader_va = Data.DataLoader(\n",
    "        dataset=training_dataset_va,\n",
    "        batch_size=batch_size_training,\n",
    "        shuffle=False,\n",
    "    )\n",
    "\n",
    "# test data\n",
    "xtest = x_tensor[Ntrain: Nsample]\n",
    "yvmtest = yvm_tensor[Ntrain: Nsample]\n",
    "yvatest = yva_tensor[Ntrain: Nsample]\n",
    "\n",
    "# batch data\n",
    "batch_size_test = 1\n",
    "test_dataset_vm = Data.TensorDataset(xtest, yvmtest)\n",
    "test_loader_vm = Data.DataLoader(\n",
    "        dataset=test_dataset_vm,\n",
    "        batch_size=batch_size_test,\n",
    "        shuffle=False,\n",
    "    )             \n",
    "\n",
    "test_dataset_va = Data.TensorDataset(xtest, yvatest)\n",
    "test_loader_va = Data.DataLoader(\n",
    "        dataset=test_dataset_va,\n",
    "        batch_size=batch_size_test,\n",
    "        shuffle=False,\n",
    "    )\n",
    "\n",
    "print('yvmtrain',torch.min(yvmtrain), torch.max(yvmtrain))\n",
    "print('yvatrain',torch.min(yvatrain), torch.max(yvatrain))\n",
    "\n",
    "# test load generation\n",
    "Pdtest = RPd[Ntrain: Nsample]/baseMVA\n",
    "Qdtest = RQd[Ntrain: Nsample]/baseMVA\n",
    "Pgtest = RPg[Ntrain: Nsample, idxPg]/baseMVA\n",
    "Qgtest = RQg[Ntrain: Nsample, idxQg]/baseMVA\n",
    "Qgtest = Qgtest.squeeze()\n",
    "\n",
    "# historical voltage\n",
    "Real_Vmtrain = yvmtrain/scale_vm*(VmUb - VmLb) + VmLb\n",
    "Real_Vatrain = yvatrain/scale_va\n",
    "hisVm_max,_ = torch.max(Real_Vmtrain, dim=0)\n",
    "hisVm_min,_ = torch.min(Real_Vmtrain, dim=0)\n",
    "his_Va = np.mean(np.insert(Real_Vatrain.numpy(), bus_slack, values = 0, axis=1),axis=0)\n",
    "his_Vm = np.mean(Real_Vmtrain.numpy(), axis=0)\n",
    "his_V = his_Vm * np.exp(1j * his_Va) # historical V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## other function\n",
    "def get_clamp(Pred, Predmin, Predmax):\n",
    "    # each row is a sample;Predmin and Predmax is the limit for each element of each row\n",
    "    Pred_clip = Pred.clone()\n",
    "    for i in range(Pred.shape[1]):\n",
    "        Pred_clip[:,i] = Pred_clip[:,i].clamp(min = Predmin[i])\n",
    "        Pred_clip[:,i] = Pred_clip[:,i].clamp(max = Predmax[i])\n",
    "    \n",
    "    return Pred_clip\n",
    "\n",
    "# testing Vm\n",
    "def get_mae(real, predict):\n",
    "    '''\n",
    "    mean absolute error\n",
    "    '''\n",
    "    if len(real) == len(predict):\n",
    "        err = torch.mean(torch.abs(real - predict))  \n",
    "        return err\n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "def get_rerr(real, predict):\n",
    "    '''\n",
    "    relative error\n",
    "    '''\n",
    "    if len(real) == len(predict):\n",
    "        err = torch.abs((predict - real) / real)*100\n",
    "        return err\n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "def get_rerr2(real, predict):\n",
    "    '''\n",
    "    relative error\n",
    "    '''\n",
    "    if len(real) == len(predict):\n",
    "        err = (predict - real) / real*100\n",
    "        return err\n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "# power balance\n",
    "def get_PQ(V):\n",
    "    S = np.zeros(V.shape,dtype = 'complex_')\n",
    "    for i in range(V.shape[0]):\n",
    "        I = Ybus.dot(V[i]).conj()\n",
    "        S[i]  = np.multiply(V[i], I)\n",
    "    \n",
    "    P = np.real(S)\n",
    "    Q = np.imag(S) \n",
    "    return P, Q\n",
    "\n",
    "def get_genload(V, Pdtest, Qdtest, bus_Pg, bus_Qg):\n",
    "    S = np.zeros(V.shape,dtype = 'complex_')\n",
    "    for i in range(V.shape[0]):\n",
    "        I = Ybus.dot(V[i]).conj()\n",
    "        S[i]  = np.multiply(V[i], I)\n",
    "    \n",
    "    P = np.real(S)\n",
    "    Q = np.imag(S) \n",
    "    \n",
    "    Pg = P[:, bus_Pg] + Pdtest[:, bus_Pg]\n",
    "    Qg = Q[:, bus_Qg] + Qdtest[:, bus_Qg]   \n",
    "    Pd = -P*1.0\n",
    "    Qd = -Q*1.0\n",
    "    Pd[:, bus_Pg] = Pg - P[:, bus_Pg]\n",
    "    Qd[:, bus_Qg] = Qg - Q[:, bus_Qg]   \n",
    "    return Pg, Qg, Pd, Qd\n",
    "\n",
    "\n",
    "# cost\n",
    "def get_Pgcost(Pg,idxPg,gencost):\n",
    "    cost = np.zeros(Pg.shape[0])\n",
    "    PgMVA = Pg*baseMVA\n",
    "    PgMVA = Pg*baseMVA\n",
    "    for i in range(Pg.shape[0]): \n",
    "        c1 = np.multiply(gencost[idxPg, 0], np.multiply(PgMVA[i,:], PgMVA[i,:]))\n",
    "        c2 = np.multiply(gencost[idxPg, 1], PgMVA[i,:])\n",
    "        cost[i] = np.sum(c1 + c2)\n",
    "            \n",
    "    return cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## NN function\n",
    "class NetVa(nn.Module):\n",
    "    def __init__(self, input_channels, output_channels, hidden_units,khidden):\n",
    "        super(NetVa, self).__init__()\n",
    "        #''' \n",
    "        self.num_layer = khidden.shape[0]\n",
    "        \n",
    "        self.fc1 = nn.Linear(input_channels, khidden[0]*hidden_units)  \n",
    "        if self.num_layer >= 2: \n",
    "            self.fc2 = nn.Linear(khidden[0]*hidden_units, khidden[1]*hidden_units)\n",
    "        \n",
    "        if self.num_layer >= 3:\n",
    "            self.fc3 = nn.Linear(khidden[1]*hidden_units, khidden[2]*hidden_units)\n",
    "        \n",
    "        if self.num_layer >= 4:\n",
    "            self.fc4 = nn.Linear(khidden[2]*hidden_units, khidden[3]*hidden_units)\n",
    "            \n",
    "        if self.num_layer >= 5:\n",
    "            self.fc5 = nn.Linear(khidden[3]*hidden_units, khidden[4]*hidden_units)\n",
    "            \n",
    "        if self.num_layer >= 6:\n",
    "            self.fc6 = nn.Linear(khidden[4]*hidden_units, khidden[5]*hidden_units)\n",
    "            \n",
    "            \n",
    "        self.fcbfend = nn.Linear(khidden[khidden.shape[0]-1]*hidden_units, output_channels)   \n",
    "        self.fcend = nn.Linear(output_channels, output_channels)  \n",
    "\n",
    "    def forward(self, x):\n",
    "        x = F.relu(self.fc1(x))\n",
    "        \n",
    "        if self.num_layer >= 2:\n",
    "            x = F.relu(self.fc2(x))\n",
    "            \n",
    "        if self.num_layer >= 3:\n",
    "            x = F.relu(self.fc3(x))\n",
    "            \n",
    "        if self.num_layer >= 4:\n",
    "            x = F.relu(self.fc4(x))\n",
    "        \n",
    "        if self.num_layer >= 5:\n",
    "            x = F.relu(self.fc5(x))\n",
    "            \n",
    "        if self.num_layer >= 6:\n",
    "            x = F.relu(self.fc6(x))\n",
    "        \n",
    "        # fixed final two layers\n",
    "        x = F.relu(self.fcbfend(x))\n",
    "        x_PredVa = self.fcend(x)\n",
    "                \n",
    "        return x_PredVa\n",
    "        \n",
    "\n",
    "class NetVm(nn.Module):\n",
    "    def __init__(self, input_channels, output_channels, hidden_units,khidden):\n",
    "        super(NetVm, self).__init__()\n",
    "        self.num_layer = khidden.shape[0]\n",
    "        self.fc1 = nn.Linear(input_channels, khidden[0]*hidden_units)\n",
    "        if self.num_layer >= 2: \n",
    "            self.fc2 = nn.Linear(khidden[0]*hidden_units, khidden[1]*hidden_units)\n",
    "        \n",
    "        if self.num_layer >= 3:\n",
    "            self.fc3 = nn.Linear(khidden[1]*hidden_units, khidden[2]*hidden_units)\n",
    "        \n",
    "        if self.num_layer >= 4:\n",
    "            self.fc4 = nn.Linear(khidden[2]*hidden_units, khidden[3]*hidden_units)\n",
    "            \n",
    "        if self.num_layer >= 5:\n",
    "            self.fc5 = nn.Linear(khidden[3]*hidden_units, khidden[4]*hidden_units)\n",
    "            \n",
    "        if self.num_layer >= 6:\n",
    "            self.fc6 = nn.Linear(khidden[4]*hidden_units, khidden[5]*hidden_units)\n",
    "            \n",
    "            \n",
    "        self.fcbfend = nn.Linear(khidden[khidden.shape[0]-1]*hidden_units, output_channels)   \n",
    "        self.fcend = nn.Linear(output_channels, output_channels) \n",
    "        \n",
    "    def forward(self, x):\n",
    "        x = F.relu(self.fc1(x))\n",
    "        \n",
    "        if self.num_layer >= 2:\n",
    "            x = F.relu(self.fc2(x))\n",
    "            \n",
    "        if self.num_layer >= 3:\n",
    "            x = F.relu(self.fc3(x))\n",
    "            \n",
    "        if self.num_layer >= 4:\n",
    "            x = F.relu(self.fc4(x))\n",
    "        \n",
    "        if self.num_layer >= 5:\n",
    "            x = F.relu(self.fc5(x))\n",
    "            \n",
    "        if self.num_layer >= 6:\n",
    "            x = F.relu(self.fc6(x))\n",
    "        \n",
    "        # fixed final two layers\n",
    "        x = F.relu(self.fcbfend(x))\n",
    "        x_PredVm = self.fcend(x)\n",
    "\n",
    "        return x_PredVm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## violation calculation function\n",
    "# Pg Qg violation\n",
    "def get_vioPQg(Pred_Pg, bus_Pg, MAXMIN_Pg, Pred_Qg, bus_Qg, MAXMIN_Qg):\n",
    "    vio_PQgmaxminnum = torch.zeros((Pred_Pg.shape[0],4))\n",
    "    vio_PQgmaxmin = torch.zeros((Pred_Pg.shape[0],4))\n",
    "    vio_PQg = torch.zeros((Pred_Pg.shape[0],2))\n",
    "    lsPg = list()\n",
    "    lsQg = list()\n",
    "    lsidxPg = np.zeros((Pred_Pg.shape[0]),dtype = np.int) # index of sample that has violation\n",
    "    lsidxQg = np.zeros((Pred_Pg.shape[0]),dtype = np.int) # index of sample that has violation\n",
    "    kP = 1 # violated samples index\n",
    "    kQ = 1 # violated samples index\n",
    "    deltaPgL = np.array([[0, 0]])\n",
    "    deltaPgU = np.array([[0, 0]])\n",
    "    deltaQgL = np.array([[0, 0]])\n",
    "    deltaQgU = np.array([[0, 0]])\n",
    "    for i in range(Pred_Pg.shape[0]):\n",
    "        # P\n",
    "        delta = Pred_Pg[i] - MAXMIN_Pg[:, 0]\n",
    "        idxPgUB = np.array(np.where(delta > DELTA))\n",
    "        if np.size(idxPgUB) > 0:\n",
    "            PgUB = np.concatenate((idxPgUB,delta[idxPgUB]),axis=0).T \n",
    "            deltaPgU = np.append(deltaPgU, PgUB, axis = 0)\n",
    "\n",
    "        delta = Pred_Pg[i] - MAXMIN_Pg[:, 1]\n",
    "        idxPgLB = np.array(np.where(delta < -DELTA))\n",
    "        if np.size(idxPgLB) > 0:\n",
    "            PgLB = np.concatenate((idxPgLB,delta[idxPgLB]),axis=0).T\n",
    "            deltaPgL = np.append(deltaPgL, PgLB, axis = 0)\n",
    "        \n",
    "        if np.size(idxPgUB)>0 and np.size(idxPgLB) > 0:\n",
    "            PgLUB = np.concatenate((PgUB,PgLB),axis=0)\n",
    "        elif np.size(idxPgUB) > 0:\n",
    "            PgLUB = PgUB\n",
    "        elif np.size(idxPgLB) > 0:\n",
    "            PgLUB = PgLB\n",
    "        \n",
    "        if (np.size(idxPgUB) + np.size(idxPgLB)) > 0:\n",
    "            PgLUB = PgLUB[PgLUB[:,0].argsort()]\n",
    "            lsPg.append(PgLUB)\n",
    "            lsidxPg[i] = kP # has violation:index of lsPg \n",
    "            kP+= 1#index=lsPg\n",
    "\n",
    "        # Q\n",
    "        delta = Pred_Qg[i] - MAXMIN_Qg[:, 0]\n",
    "        idxQgUB = np.array(np.where(delta > DELTA))\n",
    "        if np.size(idxQgUB) > 0:\n",
    "            QgUB = np.concatenate((idxQgUB,delta[idxQgUB]),axis=0).T\n",
    "            deltaQgU = np.append(deltaQgU, QgUB, axis = 0)\n",
    "\n",
    "        delta = Pred_Qg[i] - MAXMIN_Qg[:, 1]\n",
    "        idxQgLB = np.array(np.where(delta < -DELTA))\n",
    "        if np.size(idxQgLB) > 0:\n",
    "            QgLB = np.concatenate((idxQgLB,delta[idxQgLB]),axis=0).T\n",
    "            deltaQgL = np.append(deltaQgL, QgLB, axis = 0)\n",
    "            \n",
    "        if np.size(idxQgUB) >0 and np.size(idxQgLB) > 0:   \n",
    "            QgLUB = np.concatenate((QgUB,QgLB),axis=0)\n",
    "        elif np.size(idxQgUB) > 0:\n",
    "            QgLUB = QgUB\n",
    "        elif np.size(idxQgLB) > 0:\n",
    "            QgLUB = QgLB\n",
    "         \n",
    "        if (np.size(idxQgUB) + np.size(idxQgLB)) > 0:\n",
    "            QgLUB = QgLUB[QgLUB[:,0].argsort()]\n",
    "            lsQg.append(QgLUB)\n",
    "            lsidxQg[i] = kQ # has violation\n",
    "            kQ+= 1\n",
    "                    \n",
    "        \n",
    "        vio_PQgmaxminnum[i,0] = np.size(idxPgUB)\n",
    "        vio_PQgmaxminnum[i,1] = np.size(idxPgLB)\n",
    "        vio_PQgmaxminnum[i,2] = np.size(idxQgUB)\n",
    "        vio_PQgmaxminnum[i,3] = np.size(idxQgLB)\n",
    "        \n",
    "    # Pg Qg violation ratio\n",
    "    vio_PQgmaxmin[:,0] = (1 - vio_PQgmaxminnum[:, 0]/bus_Pg.shape[0])*100\n",
    "    vio_PQgmaxmin[:,1] = (1 - vio_PQgmaxminnum[:, 1]/bus_Pg.shape[0])*100\n",
    "    vio_PQgmaxmin[:,2] = (1 - vio_PQgmaxminnum[:, 2]/bus_Qg.shape[0])*100\n",
    "    vio_PQgmaxmin[:,3] = (1 - vio_PQgmaxminnum[:, 3]/bus_Qg.shape[0])*100\n",
    "    vio_PQg[:, 0] = (1 - (vio_PQgmaxminnum[:, 0] + vio_PQgmaxminnum[:, 1])/bus_Pg.shape[0])*100\n",
    "    vio_PQg[:, 1] = (1 - (vio_PQgmaxminnum[:, 2] + vio_PQgmaxminnum[:, 3])/bus_Qg.shape[0])*100\n",
    "     \n",
    "    # delete initial \n",
    "    if deltaPgL.shape[0] > 1:\n",
    "        deltaPgL = np.delete(deltaPgL, 0, axis = 0)\n",
    "        \n",
    "    if deltaPgU.shape[0] > 1:    \n",
    "        deltaPgU = np.delete(deltaPgU, 0, axis = 0)\n",
    "    \n",
    "    if deltaQgL.shape[0] > 1:\n",
    "        deltaQgL = np.delete(deltaQgL, 0, axis = 0)\n",
    "    \n",
    "    if deltaQgU.shape[0] > 1:\n",
    "        deltaQgU = np.delete(deltaQgU, 0, axis = 0)\n",
    "    \n",
    "    return lsPg,lsQg,lsidxPg,lsidxQg,vio_PQgmaxmin,vio_PQg, deltaPgL,deltaPgU,deltaQgL,deltaQgU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# modify V theta\n",
    "# dPbus_dV dQbus_dV\n",
    "def dPQbus_dV(his_V,bus_Pg,bus_Qg):\n",
    "    V = his_V.copy()\n",
    "#     Ibus = np.dot(Ybus, his_V)\n",
    "    Ibus = Ybus.dot(his_V).conj()\n",
    "    diagV = np.diag(V)\n",
    "    diagIbus = np.diag(Ibus)\n",
    "    diagVnorm = np.diag(V/np.abs(V))\n",
    "    dSbus_dVm = np.dot(diagV, Ybus.dot(diagVnorm).conj()) + np.dot(diagIbus.conj(), diagVnorm)\n",
    "    dSbus_dVa = 1j*np.dot(diagV, (diagIbus - Ybus.dot(diagV)).conj())\n",
    "    dSbus_dV = np.concatenate((dSbus_dVa, dSbus_dVm), axis=1)\n",
    "    dPbus_dV = np.real(dSbus_dV)\n",
    "    dQbus_dV = np.imag(dSbus_dV)\n",
    "    \n",
    "    return dPbus_dV,dQbus_dV\n",
    "\n",
    "# calculate dV using historical V\n",
    "def get_hisdV(lsPg,lsQg,lsidxPg,lsidxQg,num_viotest,k_dV,bus_Pg,bus_Qg,dPbus_dV,dQbus_dV):\n",
    "    dV = np.zeros((num_viotest,Nbus*2))\n",
    "    j = 0\n",
    "    for i in range(Ntest):\n",
    "        # determin whether there is violation\n",
    "        if (lsidxPg[i] + lsidxQg[i]) > 0:\n",
    "            # calculate dS_dV, dPbus_dV,dQbus_dV\n",
    "            if lsidxPg[i] >0 and lsidxQg[i] >0:\n",
    "                idxPg = lsPg[lsidxPg[i]-1][:,0].astype(np.int32) #note: lsidxPg[i]-1\n",
    "                idxQg = lsQg[lsidxQg[i]-1][:,0].astype(np.int32) #note: lsidxQg[i]-1\n",
    "                busPg = bus_Pg[idxPg]\n",
    "                busQg = bus_Qg[idxQg]\n",
    "                dPQGbus_dV = np.concatenate((dPbus_dV[busPg, :], dQbus_dV[busQg, :]), axis=0) #need bus number pf Pg Qg\n",
    "                dPQg = np.concatenate((lsPg[lsidxPg[i]-1][:,1], lsQg[lsidxQg[i]-1][:,1]), axis=0)\n",
    "            elif lsidxPg[i] >0:  \n",
    "                idxPg = lsPg[lsidxPg[i]-1][:,0].astype(np.int32)\n",
    "                busPg = bus_Pg[idxPg]\n",
    "                dPQGbus_dV = dPbus_dV[busPg, :]\n",
    "                dPQg = lsPg[lsidxPg[i]-1][:,1]\n",
    "            elif lsidxQg[i] >0:     \n",
    "                idxQg = lsQg[lsidxQg[i]-1][:,0].astype(np.int32)\n",
    "                busQg = bus_Qg[idxQg]\n",
    "                dPQGbus_dV = dQbus_dV[busQg, :]\n",
    "                dPQg = lsQg[lsidxQg[i]-1][:,1]\n",
    "\n",
    "            dV[j] = np.dot(np.linalg.pinv(dPQGbus_dV), dPQg*k_dV)           \n",
    "            j+=1   \n",
    "            \n",
    "    return dV\n",
    "\n",
    "# # calculate dV using predicted V\n",
    "# lsidxtest is the test sample that has violation\n",
    "def get_dV(Pred_V,lsPg,lsQg,lsidxPg,lsidxQg,num_viotest,k_dV,bus_Pg,bus_Qg):\n",
    "    dV = np.zeros((num_viotest,Pred_V.shape[1]*2))\n",
    "    j = 0\n",
    "    for i in range(Pred_V.shape[0]):\n",
    "        # determin whether there is violation\n",
    "        if (lsidxPg[i] + lsidxQg[i]) >0:\n",
    "            # dS_dV\n",
    "            V = Pred_V[i].copy()\n",
    "            Ibus = Ybus.dot(his_V).conj()\n",
    "            diagV = np.diag(V)\n",
    "            diagIbus = np.diag(Ibus)\n",
    "            diagVnorm = np.diag(V/np.abs(V))\n",
    "\n",
    "            dSbus_dVm = np.dot(diagV, Ybus.dot(diagVnorm).conj()) + np.dot(diagIbus.conj(), diagVnorm)\n",
    "            dSbus_dVa = 1j*np.dot(diagV, (diagIbus - Ybus.dot(diagV)).conj())\n",
    "            \n",
    "            dSbus_dV = np.concatenate((dSbus_dVa, dSbus_dVm), axis=1)\n",
    "            dPbus_dV = np.real(dSbus_dV)\n",
    "            dQbus_dV = np.imag(dSbus_dV)\n",
    "            # get BusPg\n",
    "            if lsidxPg[i] >0 and lsidxQg[i] >0:\n",
    "                idxPg = lsPg[lsidxPg[i]-1][:,0].astype(np.int32) #note: lsidxPg[i]-1\n",
    "                idxQg = lsQg[lsidxQg[i]-1][:,0].astype(np.int32) #note: lsidxQg[i]-1\n",
    "                busPg = bus_Pg[idxPg]\n",
    "                busQg = bus_Qg[idxQg]\n",
    "                dPQGbus_dV = np.concatenate((dPbus_dV[busPg, :], dQbus_dV[busQg, :]), axis=0) #need bus number pf Pg Qg\n",
    "                dPQg = np.concatenate((lsPg[lsidxPg[i]-1][:,1], lsQg[lsidxQg[i]-1][:,1]), axis=0)\n",
    "            elif lsidxPg[i] >0:  \n",
    "                idxPg = lsPg[lsidxPg[i]-1][:,0].astype(np.int32)\n",
    "                busPg = bus_Pg[idxPg]\n",
    "                dPQGbus_dV = dPbus_dV[busPg, :]\n",
    "                dPQg = lsPg[lsidxPg[i]-1][:,1]\n",
    "            elif lsidxQg[i] >0:     \n",
    "                # get BusQg\n",
    "                idxQg = lsQg[lsidxQg[i]-1][:,0].astype(np.int32)\n",
    "                busQg = bus_Qg[idxQg]\n",
    "                dPQGbus_dV = dQbus_dV[busQg, :]\n",
    "                dPQg = lsQg[lsidxQg[i]-1][:,1]\n",
    "\n",
    "            dV[j] = np.dot(np.linalg.pinv(dPQGbus_dV), dPQg*k_dV)           \n",
    "            j+=1        \n",
    "    return dV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# branch constraints\n",
    "def get_viobran(Pred_V,Pred_Va,branch,Yf,Yt):\n",
    "   #branch format: [branch from bus, branch to bus, branch limit, minang, maxang]\n",
    "    branlp = branch[:, 2]/baseMVA\n",
    "    angminmax = branch[:, 3:5]*math.pi/180\n",
    "    Pred_branang = Pred_Va[:, BRANFT[:, 0]] - Pred_Va[:, BRANFT[:, 1]]\n",
    "    vio_branangnum = torch.zeros(Pred_V.shape[0])\n",
    "    vio_branpfnum = torch.zeros(Pred_V.shape[0])\n",
    "    deltapf = np.array([[0, 0]])\n",
    "    for i in range(Pred_V.shape[0]):\n",
    "    #for i in range(1):\n",
    "        vio_branangnum[i] = np.size(np.where(Pred_branang[i, :] - angminmax[:,0] < -DELTA)) \\\n",
    "        + np.size(np.where(Pred_branang[i, :] - angminmax[:,1] > DELTA))\n",
    "\n",
    "        # branch power flow\n",
    "        fV = Pred_V[i, BRANFT[:, 0]]\n",
    "        tV = Pred_V[i, BRANFT[:, 1]]\n",
    "        fI = Yf.dot(Pred_V[i]).conj()\n",
    "        tI = Yt.dot(Pred_V[i]).conj()       \n",
    "        fS = np.multiply(fV, fI)\n",
    "        tS = np.multiply(tV, tI)\n",
    "        deltafS = np.abs(fS) - branlp\n",
    "        deltatS = np.abs(tS) - branlp\n",
    "        deltafS = np.array(deltafS).ravel()\n",
    "        deltatS = np.array(deltatS).ravel()\n",
    "        idxfs = np.array(np.where(deltafS > DELTA))\n",
    "        idxts = np.array(np.where(deltatS > DELTA))\n",
    "        vio_branpfnum[i] = np.size(idxfs) + np.size(idxfs)\n",
    "        \n",
    "        # save violated samples [index-of-branch violation-degree]\n",
    "        if np.size(idxfs) >= 1:\n",
    "            ii = np.concatenate((idxfs,deltafS[idxfs]),axis=0)\n",
    "            deltapf = np.append(deltapf, ii.T, axis = 0)\n",
    "            \n",
    "        if np.size(idxts) >= 1:\n",
    "            ii = np.concatenate((idxts,deltatS[idxts]),axis=0)\n",
    "            deltapf = np.append(deltapf, ii.T, axis = 0) \n",
    "            \n",
    "    # delete initial\n",
    "    if deltapf.shape[0] > 1:\n",
    "        deltapf = np.delete(deltapf, 0, axis = 0)\n",
    "        deltapfR = deltapf[:, 1]/branlp[0,deltapf[:, 0].astype(int)]*100\n",
    "        deltapf = np.insert(deltapf, 2, values=deltapfR, axis=1)\n",
    " \n",
    "    vio_branang = (1-vio_branangnum/branch.shape[0])*100 \n",
    "    vio_branpf = (1-vio_branpfnum/(branch.shape[0]*2))*100\n",
    "    return vio_branang,vio_branpf,deltapf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# branch constraints with violation details for dV\n",
    "def get_viobran2(Pred_V,Pred_Va,branch,Yf,Yt):\n",
    "   #branch=[branch from bus, branch to bus, branch limit, minang, maxang]\n",
    "    branlp = branch[:, 2]/baseMVA\n",
    "    angminmax = branch[:, 3:5]*math.pi/180\n",
    "    Pred_branang = Pred_Va[:, BRANFT[:, 0]] - Pred_Va[:, BRANFT[:, 1]]\n",
    "    vio_branangnum = torch.zeros(Pred_V.shape[0])\n",
    "    vio_branpfnum = torch.zeros(Pred_V.shape[0])\n",
    "    vio_branpfidx = torch.zeros(Pred_V.shape[0])\n",
    "    lsSf = []\n",
    "    lsSt = []\n",
    "    lsSf_sampidx = []\n",
    "    lsSt_sampidx = []\n",
    "    deltapf = np.array([[0, 0]])\n",
    "    for i in range(Pred_V.shape[0]):\n",
    "    #for i in range(1):\n",
    "        vio_branangnum[i] = np.size(np.where(Pred_branang[i, :] - angminmax[:,0] < -DELTA)) \\\n",
    "        + np.size(np.where(Pred_branang[i, :] - angminmax[:,1] > DELTA))\n",
    "\n",
    "        # branch power flow\n",
    "        fV = Pred_V[i, BRANFT[:, 0]]\n",
    "        tV = Pred_V[i, BRANFT[:, 1]]        \n",
    "        fI = Yf.dot(Pred_V[i]).conj()\n",
    "        tI = Yt.dot(Pred_V[i]).conj()\n",
    "        fS = np.multiply(fV, fI)\n",
    "        tS = np.multiply(tV, tI)\n",
    "        deltafS = np.abs(fS) - branlp\n",
    "        deltatS = np.abs(tS) - branlp\n",
    "        deltafS = np.array(deltafS).ravel()\n",
    "        deltatS = np.array(deltatS).ravel()\n",
    "        idxfs = np.array(np.where(deltafS > DELTA)).reshape(-1, 1)\n",
    "        idxts = np.array(np.where(deltatS > DELTA)).reshape(-1, 1)\n",
    "        vio_branpfnum[i] = np.size(idxfs) + np.size(idxfs)\n",
    "        if np.size(idxfs) >= 1:           \n",
    "            ii = np.concatenate((idxfs,deltafS[idxfs]),axis=1)  \n",
    "            deltapf = np.append(deltapf, ii, axis = 0) \n",
    "            ii = np.concatenate((ii,np.real(fS[idxfs]),np.imag(fS[idxfs])),axis=1)\n",
    "            lsSf.append(ii)\n",
    "            lsSf_sampidx.append(i)           \n",
    "            \n",
    "        if np.size(idxts) >= 1:\n",
    "            ii = np.concatenate((idxts,deltatS[idxts]),axis=1)\n",
    "            deltapf = np.append(deltapf, ii, axis = 0) \n",
    "            ii = np.concatenate((ii,np.real(tS[idxts]),np.imag(tS[idxts])),axis=1)\n",
    "            lsSt.append(ii)\n",
    "            lsSt_sampidx.append(i)\n",
    "\n",
    "        if np.size(idxfs) + np.size(idxts) >= 1:\n",
    "            vio_branpfidx[i] = i+1\n",
    "            \n",
    "    # delete initial\n",
    "    if deltapf.shape[0] > 1:\n",
    "        deltapf = np.delete(deltapf, 0, axis = 0)\n",
    "        deltapfR = deltapf[:, 1]/branlp[deltapf[:, 0].astype(int)]*100\n",
    "        deltapf = np.insert(deltapf, 2, values=deltapfR, axis=1)\n",
    " \n",
    "    vio_branang = (1-vio_branangnum/branch.shape[0])*100 \n",
    "    vio_branpf = (1-vio_branpfnum/(branch.shape[0]*2))*100\n",
    "    return vio_branang,vio_branpf,deltapf,vio_branpfidx,lsSf,lsSt,lsSf_sampidx,lsSt_sampidx\n",
    "\n",
    "def dSlbus_dV(his_V,bus_Va):\n",
    "    V = his_V.copy()\n",
    "    branlp = branch[:, 2]/baseMVA\n",
    "    # branch power flow for FROM branch\n",
    "    fV = V[BRANFT[:, 0]]\n",
    "    fI = Yf.dot(V).conj()\n",
    "    Nsam = 1\n",
    "    dfP_dVm = np.zeros((Nsam, branch.shape[0], Nbus))   \n",
    "    dfQ_dVm = np.zeros((Nsam, branch.shape[0], Nbus))\n",
    "    dfP_dVa = np.zeros((Nsam, branch.shape[0], Nbus-1))    \n",
    "    dfQ_dVa = np.zeros((Nsam, branch.shape[0], Nbus-1))                \n",
    "    diagfI = np.diag(fI)\n",
    "    diagfV = np.diag(fV)\n",
    "    diagVnorm = np.diag(np.true_divide(V, np.abs(V)))\n",
    "#     print('diagfV',diagfV.shape,'diagVnorm',diagVnorm.shape,'Yf',Yf.shape,'diagfI',diagfI.shape,'finc',finc.shape)\n",
    "    dfS_dVm = np.dot(diagfV, Yf.dot(diagVnorm).conj()) + np.dot(diagfI.conj(), np.dot(finc,diagVnorm))\n",
    "    dfP_dVm = np.real(dfS_dVm)\n",
    "    dfQ_dVm = np.imag(dfS_dVm)\n",
    "    diagV = np.diag(V)\n",
    "    dfS_dVa = -1j*np.dot(diagfV, Yf.dot(diagV).conj()) + 1j*np.dot(diagfI.conj(), np.dot(finc, diagV))\n",
    "    dfP_dVa = np.real(dfS_dVa[:, bus_Va])\n",
    "    dfQ_dVa = np.imag(dfS_dVa[:, bus_Va])\n",
    "    dPfbus_dV = np.concatenate((dfP_dVa, dfP_dVm), axis = 1)\n",
    "    dQfbus_dV = np.concatenate((dfQ_dVa, dfQ_dVm), axis = 1)\n",
    "\n",
    "#     for TO branch\n",
    "#     fI = np.zeros((V.shape[0], branch.shape[0]), dtype=complex)\n",
    "#     tI = np.zeros((V.shape[0], branch.shape[0]), dtype=complex)\n",
    "#     tV = V[0, BRANFT[:, 1]]\n",
    "#     tI = Yt.dot(Yt, V).conj()\n",
    "#     tS  = np.multiply(tV, tI)    \n",
    "#     diagtI = np.diag(tI)\n",
    "#     diagtV = np.diag(tV)\n",
    "#     dtP_dVm = np.zeros((Nsam, branch.shape[0], Nbus))\n",
    "#     dtQ_dVm = np.zeros((Nsam, branch.shape[0], Nbus))\n",
    "#     dtP_dVa = np.zeros((Nsam, branch.shape[0], Nbus-1))\n",
    "#     dtQ_dVa = np.zeros((Nsam, branch.shape[0], Nbus-1))\n",
    "#     dtS_dVm = np.dot(diagtV, Yt.dot(diagVnorm).conj()) + np.dot(diagtI.conj(), np.dot(tinc,diagVnorm))\n",
    "#     dtS_dVa = -1j*np.dot(diagtV, Yt.dot(diagV).conj()) + 1j*np.dot(diagtI.conj(), np.dot(tinc, diagV))\n",
    "#     dtP_dVm = np.real(dtS_dVm)\n",
    "#     dtQ_dVm = np.imag(dtS_dVm)\n",
    "#     dtP_dVa = np.real(dtS_dVa[:, bus_Va])\n",
    "#     dtQ_dVa = np.imag(dtS_dVa[:, bus_Va])\n",
    "\n",
    "    return dPfbus_dV, dQfbus_dV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# neural setting\n",
    "input_channels = xtrain.shape[1]\n",
    "output_channels_vm = yvmtrain.shape[1]\n",
    "output_channels_va = yvatrain.shape[1]\n",
    "lossvm = [] # save training losses of Vm\n",
    "lossva = [] # save training losses of Va\n",
    "\n",
    "# determine size of hidden layers\n",
    "if x.shape[1] >= 100:\n",
    "    hidden_units = 128\n",
    "elif x.shape[1] > 30:\n",
    "    hidden_units = 64\n",
    "else:\n",
    "    hidden_units = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train model if it is not test\n",
    "if flag_test == 0:\n",
    "    model_vm = NetVm(input_channels,output_channels_vm,hidden_units,khidden_Vm)\n",
    "    optimizer_vm = torch.optim.Adam(model_vm.parameters(), lr=Lrm)\n",
    "\n",
    "    if torch.cuda.is_available():\n",
    "        model_vm.to(device)\n",
    "        print('model_vm.to(device)')\n",
    "\n",
    "    MAXMIN_V = MAXMIN_V.to(device)\n",
    "    print('*' * 5+'Vm training'+'*' * 5)\n",
    "    #Training process: Voltage magnitude\n",
    "    start_time = time.process_time()\n",
    "    for epoch in range(EpochVm):\n",
    "        running_loss = 0.0\n",
    "        for step, (train_x, train_y) in enumerate(training_loader_vm):\n",
    "            #feedforward\n",
    "            train_x, train_y= train_x.to(device), train_y.to(device)\n",
    "            yvmtrain_hat = model_vm(train_x)    \n",
    "\n",
    "            # if epoch less than specified number/no penalty of V: only MSEloss\n",
    "            loss = criterion(train_y, yvmtrain_hat)\n",
    "            running_loss =  running_loss + loss.item()\n",
    "            \n",
    "            # backproprogate\n",
    "            optimizer_vm.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer_vm.step()      \n",
    "\n",
    "        lossvm.append(running_loss)\n",
    "\n",
    "        if (epoch+1)%p_epoch == 0:          \n",
    "            print('epoch', epoch+1, running_loss,torch.min(yvmtrain_hat).detach(),torch.max(yvmtrain_hat).detach())\n",
    "         \n",
    "        # save trianed model\n",
    "        if (epoch+1)%100 == 0 and (epoch+1) >= s_epoch:            \n",
    "            torch.save(model_vm.state_dict(), PATHVms+'E'+str(epoch+1)+'F'+str(flagVm)+'.pth',_use_new_zipfile_serialization=False)\n",
    "            \n",
    "    time_trianVm = time.process_time()-start_time  \n",
    "    print(\"\\n\")        \n",
    "    print('time_trianVm',round(time_trianVm,5),'seconds')\n",
    "    print(\"\\n\")\n",
    "    \n",
    "    # save trianed model\n",
    "    torch.save(model_vm.state_dict(), PATHVm,_use_new_zipfile_serialization=False)\n",
    "    # plot training loss \n",
    "    fig = plt.figure()\n",
    "    plt.plot(lossvm)\n",
    "    plt.title('lossvm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# testing Vm using training data\n",
    "if flag_test == 0:\n",
    "    print('test training data')\n",
    "    xtrain = xtrain.to(device) \n",
    "    yvmtrain_hat = model_vm(xtrain)\n",
    "    yvmtrain_hat = yvmtrain_hat.cpu()\n",
    "    yvmtrains = yvmtrain/scale_vm*(VmUb - VmLb) + VmLb\n",
    "    yvmtrain_hats = yvmtrain_hat.detach()/scale_vm*(VmUb - VmLb) + VmLb\n",
    "    yvmtrain_hat_clip = get_clamp(yvmtrain_hats, hisVm_min, hisVm_max)\n",
    "\n",
    "    mae_Vmtrain = get_mae(yvmtrains, yvmtrain_hat_clip.detach())\n",
    "    mre_Vmtrain = get_rerr(yvmtrains,yvmtrain_hats.detach())\n",
    "    mre_Vmtrain_clip = get_rerr(yvmtrains,yvmtrain_hat_clip.detach())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if flag_test == 0:\n",
    "    #Training process: voltage angle\n",
    "    model_va = NetVa(input_channels,output_channels_va,hidden_units,khidden_Vm)\n",
    "    optimizer_va = torch.optim.Adam(model_va.parameters(), lr=Lra)\n",
    "    model_va.to(device)\n",
    "    print('model_va.to(device)')\n",
    "    \n",
    "    start_time = time.process_time()\n",
    "    for epoch in range(EpochVa):\n",
    "        running_loss = 0.0\n",
    "        for step, (train_x, train_y) in enumerate(training_loader_va):\n",
    "            #feedforward\n",
    "            train_x, train_y = train_x.to(device), train_y.to(device)\n",
    "            yvatrain_hat = model_va(train_x)\n",
    "            loss = criterion(train_y, yvatrain_hat)\n",
    "            running_loss =  running_loss + loss.item()\n",
    "\n",
    "            # backproprogate\n",
    "            optimizer_va.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer_va.step()\n",
    "\n",
    "        lossva.append(running_loss)\n",
    "\n",
    "        if (epoch+1)%p_epoch == 0:          \n",
    "            print('epoch', epoch+1, running_loss, torch.min(yvatrain_hat).detach(),torch.max(yvatrain_hat).detach())\n",
    "\n",
    "        # save trianed model\n",
    "        if (epoch+1)%100 == 0 and (epoch+1)>= s_epoch:           \n",
    "            torch.save(model_va.state_dict(), PATHVas +'E'+str(epoch+1)+'F'+str(flagVa)+'.pth',_use_new_zipfile_serialization=False)\n",
    "\n",
    "    time_trianVa = time.process_time() - start_time  \n",
    "    \n",
    "    print(\"\\n\")        \n",
    "    print('time_trianVa',round(time_trianVa,5), 'seconds')\n",
    "    print(\"\\n\")  \n",
    "    \n",
    "    # save trianed model\n",
    "    torch.save(model_va.state_dict(), PATHVa, _use_new_zipfile_serialization=False)\n",
    "    fig = plt.figure()\n",
    "    plt.plot(lossva)\n",
    "    plt.title('lossva')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# testing Va using training data\n",
    "if flag_test == 0:\n",
    "    print('test trianing data')\n",
    "    xtrain = xtrain.to(device)\n",
    "    yvatrain_hat = model_va(xtrain)\n",
    "    yvatrain_hat = yvatrain_hat.cpu()\n",
    "    yvatrains = yvatrain/scale_va\n",
    "    yvatrain_hats = yvatrain_hat.detach()/scale_va\n",
    "\n",
    "    mae_Vatrain = get_mae(yvatrains, yvatrain_hats)\n",
    "    mre_Vatrain = get_rerr(yvatrains,yvatrain_hats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## save mat data\n",
    "scale_vm = np.double(scale_vm)\n",
    "scale_va = np.double(scale_va)\n",
    "\n",
    "xtest = xtest.to(device)  \n",
    "\n",
    "# Real_Vatrain = np.insert(Real_Vatrain.numpy(), bus_slack, values = 0, axis=1)\n",
    "# Real_Vtrain = Real_Vmtrain.numpy()* np.exp(1j *Real_Vatrain)\n",
    "# Real_Ptrain, Real_Qtrain = get_PQ(Real_Vtrain)\n",
    "# Real_PQtrain = np.concatenate((Real_Ptrain, Real_Qtrain), axis = 1)\n",
    "# Real_PQtrain = torch.from_numpy(Real_PQtrain)\n",
    "\n",
    "# # train load generation for generation constraints\n",
    "# Real_Pdtrain = RPd[0: Ntrain]/baseMVA\n",
    "# Real_Qdtrain = RQd[0: Ntrain]/baseMVA\n",
    "# Real_PQdtrain = np.concatenate((Real_Pdtrain, Real_Qdtrain), axis = 1)\n",
    "\n",
    "# incidance matricx of branch from\n",
    "finc = np.zeros((branch.shape[0], Nbus), dtype = float)\n",
    "tinc = np.zeros((branch.shape[0], Nbus), dtype = float)\n",
    "for i in range(branch.shape[0]):\n",
    "    finc[i,  branch[i, 0]-1] = 1\n",
    "    tinc[i,  branch[i, 1]-1] = 1\n",
    "# real Vm Va for testing samples\n",
    "yvmtests = yvmtest/scale_vm*(VmUb - VmLb) + VmLb\n",
    "yvatests = yvatest/scale_va\n",
    "    \n",
    "# Va\n",
    "Real_Va = yvatests.clone().numpy()\n",
    "Real_Va = np.insert(Real_Va, bus_slack, values = 0, axis=1)\n",
    "\n",
    "# V\n",
    "Real_V = yvmtests.numpy() * np.exp(1j *Real_Va)\n",
    "# Pg QG\n",
    "Real_Pg, Real_Qg, Real_Pd, Real_Qd = get_genload(Real_V, Pdtest, Qdtest, bus_Pg, bus_Qg) \n",
    "\n",
    "# Jocobian matrix for Pg Qg violation\n",
    "dPbus_dV,dQbus_dV = dPQbus_dV(his_V,bus_Pg,bus_Qg)\n",
    "\n",
    "# Jocobian matrix for branch flow violation\n",
    "bus_Va = np.delete(np.arange(Nbus),bus_slack)\n",
    "dPfbus_dV, dQfbus_dV = dSlbus_dV(his_V,bus_Va)\n",
    "\n",
    "# load trained model if it is testing\n",
    "if flag_test == 1:\n",
    "    # load trained model \n",
    "    print('load trained model: model_vm_load')\n",
    "    model_vm = NetVm(input_channels,output_channels_vm,hidden_units,khidden_Vm)\n",
    "    model_vm.load_state_dict(torch.load(PATHVm, map_location=device))\n",
    "    model_vm.eval()\n",
    "    model_vm.to(device)       \n",
    "    # load trained model \n",
    "    print('load trained model: model_va_load')\n",
    "    model_va = NetVa(input_channels,output_channels_va,hidden_units,khidden_Vm)\n",
    "    model_va.load_state_dict(torch.load(PATHVa, map_location=device))\n",
    "    model_va.eval() \n",
    "    model_va.to(device)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('*'*5+'begin repeated calcualtion for time testing'+'*'*5)  \n",
    "Pred_timeVm_NN_log = np.zeros(REPEAT)\n",
    "Pred_timeVa_NN_log = np.zeros(REPEAT)\n",
    "Pred_PQgtime_log = np.zeros((REPEAT,2))\n",
    "Pred_timeVm_log = np.zeros((REPEAT,2))\n",
    "Pred_timeVa_log = np.zeros((REPEAT,2))\n",
    "for k in range(REPEAT):\n",
    "    if (k+1)%10 == 0:\n",
    "        print('REPEAT',k+1)\n",
    "            \n",
    "    # Predicted data\n",
    "    time_PredVm_NN = 0\n",
    "    yvmtest_hat = torch.zeros((Ntest, Nbus))\n",
    "    for step, (test_x, test_y) in enumerate(test_loader_vm):\n",
    "        test_x = test_x.to(device) \n",
    "        start_time = time.process_time()\n",
    "        yvmtest_hat[step] = model_vm(test_x)\n",
    "        delta_time = time.process_time() - start_time\n",
    "        time_PredVm_NN = time_PredVm_NN + delta_time\n",
    "    \n",
    "    yvmtest_hat = yvmtest_hat.cpu()\n",
    "    start_time = time.process_time()\n",
    "    yvmtest_hats = yvmtest_hat.detach()/scale_vm*(VmUb - VmLb) + VmLb\n",
    "    yvmtest_hat_clip = get_clamp(yvmtest_hats, hisVm_min, hisVm_max)\n",
    "    time_PredVm = time.process_time() - start_time\n",
    "    time_PredVm = time_PredVm + time_PredVm_NN\n",
    "\n",
    "    time_PredVa_NN = 0\n",
    "    yvatest_hat = torch.zeros((Ntest, Nbus-1))\n",
    "    for step, (test_x, test_y) in enumerate(test_loader_va):\n",
    "        test_x = test_x.to(device)\n",
    "        start_time = time.process_time()\n",
    "        yvatest_hat[step] = model_va(test_x)\n",
    "        delta_time = time.process_time() - start_time\n",
    "        time_PredVa_NN = time_PredVa_NN + delta_time\n",
    "\n",
    "    yvatest_hat = yvatest_hat.cpu()\n",
    "    start_time = time.process_time()\n",
    "    yvatest_hats = yvatest_hat.detach()/scale_va\n",
    "    time_PredVa = time.process_time() - start_time\n",
    "    time_PredVa = time_PredVa + time_PredVa_NN\n",
    "\n",
    "    # Va\n",
    "    Pred_Va = yvatest_hats.clone().numpy()\n",
    "    Pred_Va = np.insert(Pred_Va, bus_slack, values = 0, axis=1)\n",
    "    \n",
    "    ## calculate Pg QG\n",
    "    start_time = time.process_time()\n",
    "    Pred_V = yvmtest_hat_clip.clone().numpy()* np.exp(1j * Pred_Va)  # predicted V\n",
    "    Pred_Pg, Pred_Qg, Pred_Pd, Pred_Qd = get_genload(Pred_V, Pdtest, Qdtest, bus_Pg, bus_Qg)\n",
    "    Pred_PQgtime = time.process_time() - start_time    \n",
    "    # Pg Qg violation\n",
    "    lsPg,lsQg,lsidxPg,lsidxQg,vio_PQgmaxmin,vio_PQg,deltaPgL,deltaPgU,deltaQgL,deltaQgU = get_vioPQg(Pred_Pg, bus_Pg, MAXMIN_Pg, Pred_Qg, bus_Qg, MAXMIN_Qg)\n",
    "    lsidxPQg = np.squeeze(np.array(np.where((lsidxPg + lsidxQg) > 0)))\n",
    "    num_viotest = np.size(lsidxPQg) # number of violated samples\n",
    "\n",
    "    # revise power flow (from)\n",
    "    vio_branang,vio_branpf,deltapf,vio_branpfidx,lsSf,_,lsSf_sampidx,_ = get_viobran2(Pred_V,Pred_Va,branch,Yf,Yt)\n",
    "    vio_branpf_num = np.size(np.where(vio_branpfidx>0))\n",
    "    lsSf_sampidx = np.asarray(lsSf_sampidx)\n",
    "\n",
    "    # post-processing of Vm Va \n",
    "    Pred_Va1 = Pred_Va.copy()\n",
    "    Pred_Vm1 = yvmtest_hat_clip.clone().numpy()\n",
    "\n",
    "    start_time = time.process_time()\n",
    "    if flag_hisv:\n",
    "        print('\\nUse historical V to calculate dV')\n",
    "        dV1 = get_hisdV(lsPg,lsQg,lsidxPg,lsidxQg,num_viotest,k_dV,bus_Pg,bus_Qg,dPbus_dV,dQbus_dV)\n",
    "    else:\n",
    "        print('\\nUse predicted V to calculate dV')\n",
    "        dV1 = get_dV(Pred_V,lsPg,lsQg,lsidxPg,lsidxQg,num_viotest,k_dV,bus_Pg,bus_Qg)\n",
    "\n",
    "    if vio_branpf_num > 0:\n",
    "        dV_branch = np.zeros((lsSf_sampidx.shape[0], Nbus*2))\n",
    "        start_time = time.process_time()\n",
    "        for i in range(lsSf_sampidx.shape[0]):\n",
    "            # Pl/deltaSl\n",
    "            mp = np.array(lsSf[i][:, 2]/lsSf[i][:, 1]).reshape(-1, 1)\n",
    "            mq = np.array(lsSf[i][:, 3]/lsSf[i][:, 1]).reshape(-1, 1)\n",
    "            # dPf_dVaVm\n",
    "            dPdV = dPfbus_dV[np.array(lsSf[i][:, 0].astype(int)).squeeze(), :] \n",
    "            # dQf_dVaVm\n",
    "            dQdV = dQfbus_dV[np.array(lsSf[i][:, 0].astype(int)).squeeze(), :]\n",
    "            dmp = mp*dPdV \n",
    "            dmq = mq*dQdV\n",
    "            # M: dS_dVaVm\n",
    "            dmpq_inv = np.linalg.pinv(dmp+dmq)\n",
    "            dV_branch[i] = np.dot(dmpq_inv, np.array(lsSf[i][:, 1])).squeeze()\n",
    "        dV1 = dV1 + dV_branch\n",
    "\n",
    "    # dV\n",
    "    Pred_Va1[lsidxPQg, :] = Pred_Va[lsidxPQg, :] - dV1[:, 0:Nbus]\n",
    "    Pred_Va1[:,bus_slack] = 0\n",
    "    Pred_Vm1[lsidxPQg, :] = yvmtest_hat_clip.numpy()[lsidxPQg, :] - dV1[:, Nbus:2*Nbus]\n",
    "    Pred_Vm1_clip = get_clamp(torch.from_numpy(Pred_Vm1), hisVm_min, hisVm_max)\n",
    "    Pred_V1 = Pred_Vm1_clip.numpy() * np.exp(1j*Pred_Va1) # revised V\n",
    "    Pred_Pg1, Pred_Qg1, Pred_Pd1, Pred_Qd1 = get_genload(Pred_V1, Pdtest, Qdtest, bus_Pg, bus_Qg)\n",
    "    Pred_PQgtime1 = time.process_time() - start_time\n",
    "\n",
    "    # save time\n",
    "    Pred_timeVm_NN_log[k] = time_PredVm_NN/Ntest\n",
    "    Pred_timeVa_NN_log[k] = time_PredVa_NN/Ntest\n",
    "    Pred_PQgtime_log[k,0] = Pred_PQgtime/Ntest\n",
    "    Pred_PQgtime_log[k,1] = Pred_PQgtime1/Ntest\n",
    "    Pred_timeVm_log[k,0] = (Pred_PQgtime + time_PredVm)/Ntest\n",
    "    Pred_timeVa_log[k,0] = (Pred_PQgtime + time_PredVa)/Ntest\n",
    "    Pred_timeVm_log[k,1] = (Pred_PQgtime + Pred_PQgtime1 + time_PredVm)/Ntest\n",
    "    Pred_timeVa_log[k,1] = (Pred_PQgtime + Pred_PQgtime1 + time_PredVa)/Ntest\n",
    "\n",
    "print('*'*5+'end repeated calcualtion for time testing'+'*'*5) \n",
    "\n",
    "# performance evaluation\n",
    "# no revison\n",
    "mae_Vmtest = get_mae(yvmtests, yvmtest_hat_clip.detach())\n",
    "mre_Vmtest_clip = get_rerr(yvmtests,yvmtest_hat_clip.detach())\n",
    "mae_Vatest = get_mae(yvatests, yvatest_hats)\n",
    "mre_Vatest = get_rerr(yvatests, yvatest_hats)\n",
    "# after post-processing\n",
    "mae_Vmtest1 = get_mae(yvmtests, Pred_Vm1_clip)\n",
    "mae_Vatest1 = get_mae(torch.from_numpy(Real_Va).float(), torch.from_numpy(Pred_Va1).float())\n",
    "\n",
    "# load satisfaction\n",
    "mre_Pd = get_rerr(torch.from_numpy(Real_Pd.sum(axis=1)), torch.from_numpy(Pred_Pd.sum(axis=1)))\n",
    "mre_Qd = get_rerr(torch.from_numpy(Real_Qd.sum(axis=1)), torch.from_numpy(Pred_Qd.sum(axis=1)))\n",
    "mre_Pd1 = get_rerr(torch.from_numpy(Real_Pd.sum(axis=1)), torch.from_numpy(Pred_Pd1.sum(axis=1)))\n",
    "mre_Qd1 = get_rerr(torch.from_numpy(Real_Qd.sum(axis=1)), torch.from_numpy(Pred_Qd1.sum(axis=1)))\n",
    "\n",
    "# optimality loss\n",
    "# no revison\n",
    "Pred_cost = get_Pgcost(Pred_Pg,idxPg,gencost)\n",
    "Real_cost = get_Pgcost(Real_Pg,idxPg,gencost)\n",
    "mre_cost = get_rerr2(torch.from_numpy(Real_cost), torch.from_numpy(Pred_cost))\n",
    "# after post-processing\n",
    "Pred_cost1 = get_Pgcost(Pred_Pg1,idxPg,gencost)\n",
    "mre_cost1 = get_rerr2(torch.from_numpy(Real_cost), torch.from_numpy(Pred_cost1))\n",
    "\n",
    "# branch violation\n",
    "vio_branang,vio_branpf,deltapf = get_viobran(Pred_V,Pred_Va,branch,Yf,Yt)\n",
    "if deltapf.shape[1] > 2:\n",
    "    res_branpf = np.array((np.mean(deltapf[:,1]),np.max(deltapf[:,1]),np.mean(deltapf[:,2]),np.max(deltapf[:,2])))\n",
    "else:\n",
    "    res_branpf = np.array((np.mean(deltapf[:,1]),np.max(deltapf[:,1])))\n",
    "\n",
    "# Pg Qg violation degree\n",
    "res_PQUL = np.zeros((2, 8))\n",
    "res_PQUL[0] = np.array((np.mean(deltaPgL[:,1]),np.min(deltaPgL[:,1]),np.mean(deltaPgU[:,1]),np.max(deltaPgU[:,1]), \\\n",
    "                          np.mean(deltaQgL[:,1]),np.min(deltaQgL[:,1]),np.mean(deltaQgU[:,1]),np.max(deltaQgU[:,1])))\n",
    "\n",
    "# Pg Qg violation after post-processing\n",
    "_,_,lsidxPg1,lsidxQg1,vio_PQgmaxmin1,vio_PQg1,deltaPgL1,deltaPgU1,deltaQgL1,deltaQgU1 = get_vioPQg(Pred_Pg1, bus_Pg, MAXMIN_Pg, Pred_Qg1, bus_Qg, MAXMIN_Qg)\n",
    "lsidxPQg1 = np.squeeze(np.array(np.where(lsidxPg1 + lsidxQg1> 0)))\n",
    "num_viotest1 = np.size(lsidxPQg)\n",
    "# Pg Qg violation degree after revision\n",
    "res_PQUL[1] = np.array((np.mean(deltaPgL1[:,1]),np.min(deltaPgL1[:,1]),np.mean(deltaPgU1[:,1]),np.max(deltaPgU1[:,1]), \\\n",
    "                          np.mean(deltaQgL1[:,1]),np.min(deltaQgL1[:,1]),np.mean(deltaQgU1[:,1]),np.max(deltaQgU1[:,1])))\n",
    "\n",
    "# branch violation after post-processing\n",
    "vio_branang1,vio_branpf1,deltapf1 = get_viobran(Pred_V1,Pred_Va1,branch,Yf,Yt)\n",
    "if deltapf1.shape[1] > 2:\n",
    "    res_branpf1 = np.array((np.mean(deltapf1[:,1]),np.max(deltapf1[:,1]),np.mean(deltapf1[:,2]),np.max(deltapf1[:,2])))\n",
    "else:\n",
    "    res_branpf1 = np.array((np.mean(deltapf1[:,1]),np.max(deltapf1[:,1])))\n",
    "\n",
    "# save results\n",
    "resvio = torch.zeros(2, 7)\n",
    "resvio1 = torch.zeros(2, 7)\n",
    "resvio[0] = torch.tensor([torch.mean(mre_cost),torch.mean(mre_Pd),torch.mean(mre_Qd), torch.mean(vio_PQg[:, 0]), torch.mean(vio_PQg[:, 1]), \\\n",
    "                         torch.mean(vio_branang), torch.mean(vio_branpf)]) \n",
    "resvio[1] = torch.tensor([torch.mean(mre_cost1),torch.mean(mre_Pd1),torch.mean(mre_Qd1), torch.mean(vio_PQg1[:, 0]), torch.mean(vio_PQg1[:, 1]), \\\n",
    "                         torch.mean(vio_branang1), torch.mean(vio_branpf1)]) \n",
    "resvio1[0] = torch.tensor([torch.max(mre_cost),torch.max(mre_Pd),torch.max(mre_Qd), torch.min(vio_PQg[:, 0]), torch.min(vio_PQg[:, 1]), \\\n",
    "                         torch.min(vio_branang), torch.min(vio_branpf)])\n",
    "resvio1[1] = torch.tensor([torch.max(mre_cost1),torch.max(mre_Pd1),torch.max(mre_Qd1), torch.min(vio_PQg1[:, 0]), torch.min(vio_PQg1[:, 1]), \\\n",
    "                         torch.min(vio_branang1), torch.min(vio_branpf1)]) \n",
    "\n",
    "resvio = np.around(resvio.numpy(), decimals=2)\n",
    "resvio1 = np.around(resvio1.numpy(), decimals=2)\n",
    "maeV = np.double(torch.tensor([[mae_Vmtest, mae_Vatest],[mae_Vmtest1, mae_Vatest1]]).numpy())\n",
    "lossvm = np.array(lossvm)\n",
    "lossva = np.array(lossva)\n",
    "PredtimeVm = np.mean(Pred_timeVm_log, axis=0)\n",
    "PredtimeVa = np.mean(Pred_timeVa_log, axis=0)\n",
    "mre_cost = np.array(mre_cost)\n",
    "mre_cost1 = np.array(mre_cost1)\n",
    "scipy.io.savemat(resultnm, \\\n",
    "mdict={'resvio': resvio,'resvio1': resvio1,'maeV': maeV,'res_PQUL': res_PQUL,'deltaPgL': deltaPgL,'deltaPgU': deltaPgU,'deltaQgL': deltaQgL,'deltaQgU': deltaQgU, \\\n",
    "    'PredtimeVm': PredtimeVm,'PredtimeVa': PredtimeVa,'Pred_timeVm_log': Pred_timeVm_log,'Pred_timeVa_log': Pred_timeVa_log, \\\n",
    "    'deltaPgL1': deltaPgL1,'deltaPgU1': deltaPgU1,'deltaQgL1': deltaQgL1,'deltaQgU1': deltaQgU1,'deltapf':deltapf,'deltapf1':deltapf1, \\\n",
    "    'res_branpf':res_branpf,'res_branpf1':res_branpf1,'mre_cost':mre_cost,'mre_cost1':mre_cost1,'Pred_timeVm_NN_log':Pred_timeVm_NN_log, \\\n",
    "    'Pred_timeVa_NN_log':Pred_timeVa_NN_log,'Pred_PQgtime_log':Pred_PQgtime_log\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('***************Summary result case',Nbus, 'sys_R', sys_R,'***************')\n",
    "print('training setting: Ntrain',Ntrain,' Ntest',Ntest, 'batch_size',batch_size_training,\\\n",
    "      'Lm', Lm,'La', La,' Lrm', Lrm,' Lra', Lra,' EpochVm',EpochVm,'scale_va',scale_va,'  scale_vm',scale_vm)\n",
    "print('NN layer: hidden_units',hidden_units,' khidden_Vm',khidden_Vm*hidden_units ,'khidden_Va',khidden_Va*hidden_units )\n",
    "print(model_vm)\n",
    "print(model_va)\n",
    "print('\\n','*'*20,'before modification: ','*'*20)\n",
    "print('mae_Vmtest', mae_Vmtest,'mae_Vatest', mae_Vatest)\n",
    "print('num_viotest', num_viotest)\n",
    "print('vio_Pgmax: mean', torch.mean(vio_PQgmaxmin[:,0]).numpy(), '%  max', torch.max(vio_PQgmaxmin[:,0]).numpy(), \\\n",
    "      'vio_Pgmin: mean',torch.mean(vio_PQgmaxmin[:,1]).numpy(), '%  max', torch.max(vio_PQgmaxmin[:,1]).numpy())\n",
    "print('vio_Qgmax: mean', torch.mean(vio_PQgmaxmin[:,2]).numpy(), '%  max', torch.max(vio_PQgmaxmin[:,2]), \\\n",
    "      'vio_Qgmin: mean',torch.mean(vio_PQgmaxmin[:,3]).numpy(), '%  max', torch.max(vio_PQgmaxmin[:,3]).numpy())\n",
    "print('vio_Pgmaxmin: mean', torch.mean(vio_PQg[:, 0]).numpy(), '%  max', torch.max(vio_PQg[:, 0]).numpy(), \\\n",
    "      'vio_Qgmaxmin: mean',torch.mean(vio_PQg[:, 1]).numpy(), '%  max', torch.max(vio_PQg[:, 1]).numpy())\n",
    "print('mean', torch.mean(vio_branang).numpy(), '%  max', torch.max(vio_branang).numpy(),'% ', \\\n",
    "      'vio_branpf: mean',torch.mean(vio_branpf).numpy(), '%  max', torch.max(vio_branpf).numpy(),'% ')\n",
    "print('mean mre_cost', resvio[0,0],'%  ','mre_Pd',resvio[0,1],'%  ','mre_Qd',resvio[0,2],'vio_PQg',resvio[0,3],'%  ', \\\n",
    "      'vio_PQg',resvio[0,4],'%  ','vio_branang',resvio[0,5],'%  ', 'vio_branpf',resvio[0,6],'%  ')\n",
    "\n",
    "print('\\n','*'*20,'after post-processing: ','*'*20)\n",
    "print('mae_Vmtest1', mae_Vmtest1,'mae_Vatest1', mae_Vatest1)\n",
    "print('num_viotest1', num_viotest1) \n",
    "print('vio_Pgmax: mean', torch.mean(vio_PQgmaxmin1[:,0]).numpy(), '%  max', torch.max(vio_PQgmaxmin1[:,0]).numpy(), \\\n",
    "      'vio_Pgmin: mean',torch.mean(vio_PQgmaxmin1[:,1]).numpy(), '%  max', torch.max(vio_PQgmaxmin1[:,1]).numpy())\n",
    "print('vio_Qgmax: mean', torch.mean(vio_PQgmaxmin1[:,2]).numpy(), '%  max', torch.max(vio_PQgmaxmin1[:,2]).numpy(), \\\n",
    "      'vio_Qgmin: mean',torch.mean(vio_PQgmaxmin1[:,3]).numpy(), '%  max', torch.max(vio_PQgmaxmin1[:,3]).numpy())\n",
    "print('vio_Pgmaxmin: mean', torch.mean(vio_PQg1[:, 0]).numpy(), '%  max', torch.max(vio_PQg1[:, 0]).numpy(), \\\n",
    "      'vio_Qgmaxmin: mean',torch.mean(vio_PQg1[:, 1]).numpy(), '%  max', torch.max(vio_PQg1[:, 1]).numpy())\n",
    "print('vio_branang: mean', torch.mean(vio_branang1).numpy(), '%  max', torch.max(vio_branang1).numpy(),'% ', \\\n",
    "      'vio_branpf: mean',torch.mean(vio_branpf1).numpy(), '%  max', torch.max(vio_branpf1).numpy(),'% ')\n",
    "print('mean mre_cost', resvio[1,0],'%  ','mre_Pd',resvio[1,1],'%  ','mre_Qd',resvio[1,2],'%  ', \\\n",
    "      'vio_PQg',resvio[1,3],'%  ','vio_PQg',resvio[1,4],'%  ','vio_branang',resvio1[1,5],'%  ', 'vio_branpf',resvio[1,6],'%  ')\n",
    "\n",
    "if flag_test == 0:\n",
    "    print('\\n', '*'*10,'trian vs test','*'*10)\n",
    "    print('voltage')\n",
    "    print('mae_Vmtrain', mae_Vmtrain)\n",
    "    print('mre_Vmtrain_clip', torch.mean(mre_Vmtrain_clip),'%  ',  torch.max(mre_Vmtrain_clip),'%')\n",
    "    print('mae_Vmtest', mae_Vmtest)\n",
    "    print('mre_Vmtest_clip', torch.mean(mre_Vmtest_clip),'%  ',  torch.max(mre_Vmtest_clip),'%')\n",
    "    print('angle')\n",
    "    print('mae_Vatrain', mae_Vatrain)\n",
    "    print('mre_Vatrain', torch.mean(mre_Vatrain),'%  ',  torch.max(mre_Vatrain),'%')\n",
    "    print('mae_Vatest', mae_Vatest)\n",
    "    print('mre_Vatest', torch.mean(mre_Vatest),'%  ',  torch.max(mre_Vatest),'%')\n",
    "    print('time_trianVm',round(time_trianVm,6),'sec  time_trianVa',round(time_trianVa,6), \\\n",
    "          'sec  \\ntime_trianVm',round(time_trianVm/60,6),'min  time_trianVa',round(time_trianVa/60,6),'min  \\n')\n",
    "    \n",
    "print('\\n','*'*20,'computation time','*'*20)\n",
    "print('Pred_time for each sample:\\nPred_timeVm',np.round(PredtimeVm[0],7),'sec\\nPred_timeVa',np.round(PredtimeVa[0],7),'sec')\n",
    "print('PredtimeVa',PredtimeVa,'PredtimeVm',PredtimeVm)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}