Driver-123NF_network_data_processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2.27.21\n",
    "# This 123NF driver has all the funcs, whereas driver_code_forPaper-blank.ipynb has only the 123NF for heatmap paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import importlib\n",
    "import numpy as np\n",
    "import math as m\n",
    "import statistics as st\n",
    "import cmath\n",
    "import matplotlib.pyplot as plt \n",
    "import itertools\n",
    "import random\n",
    "from operator import add\n",
    "from graphviz import Source, render\n",
    "import pydot\n",
    "\n",
    "import datetime\n",
    "import time\n",
    "\n",
    "import sys\n",
    "sys.path.append('py_modules') # below modules are in this folder\n",
    "#print(sys.path)\n",
    "import setup_nx # your own module, setup.nx.py\n",
    "importlib.reload(setup_nx)\n",
    "import my_feeder_funcs as ff\n",
    "import my_impedance_funcs as imp\n",
    "import my_configVis_funcs as vis\n",
    "import my_detControlMatExistence_funcs as ctrl\n",
    "import my_detLznRange_funcs as lzn\n",
    "import my_heatmapSetup_funcs as hm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-06-25 16:51:07.151079\n"
     ]
    }
   ],
   "source": [
    "# [Essential] specify input feeder data for IEEE 123-node test feeder\n",
    "\n",
    "# Enter the path/name of the impedance model data (excel file)\n",
    "filepath = \"feeder_impedance_models/\"\n",
    "#modelpath = filepath + \"004_GB_IEEE123_OPAL.xls\"\n",
    "#modelpath = filepath + \"004_GB_IEEE123_OPAL_stripped.xls\"\n",
    "modelpath = filepath + \"004_GB_IEEE123_OPAL_accur.xlsx\" # accurate wrt public 123NF data\n",
    "#modelpath = filepath + \"004_GB_IEEE123_OPAL_zeroCrossZ.xlsx\" # accurate wrt public 123NF data\n",
    "\n",
    "#file_name = string specifying name of dot file created when make_graph() is called\n",
    "file_name = '123NF'\n",
    "#file_name = 'demo123NF.dot'\n",
    "\n",
    "# Specify substation kV, kVA bases, name, and the number of timesteps in the load data'\n",
    "Vbase_ll = 4160\n",
    "Vbase = Vbase_ll / np.sqrt(3)\n",
    "Sbase = 5000/3\n",
    "substation_name = 'bus_150'\n",
    "timesteps = 1\n",
    "\n",
    "# initialize some variables\n",
    "ts = time.time()\n",
    "print(datetime.datetime.fromtimestamp(ts))\n",
    "plot = 0 #turn plot on/off\n",
    "depths = {}\n",
    "leaves = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished initializing feeder\n",
      "n= 127\n",
      "n= 127\n"
     ]
    }
   ],
   "source": [
    "# [ESSENTIAL] create feeder object\n",
    "\n",
    "#fin_feeder = ff.feeder_init(modelpath,loadfolder,loadpath,timesteps,Vbase_ll,Sbase,depths,leaves)\n",
    "fin_feeder = ff.feeder_init(modelpath,'','',timesteps,Vbase_ll,Sbase,depths,leaves)\n",
    "\n",
    "print(\"Finished initializing feeder\")\n",
    "ff.make_graph(fin_feeder, file_name)\n",
    "node_index_map = hm.createNodeIndexMap(fin_feeder) #node indices for indicMat and F matrix\n",
    "R,X=hm.createRXmatrices_3ph(fin_feeder, depths,file_name)\n",
    "n=round(len(R)/3)\n",
    "print('n=',n)\n",
    "\n",
    "# #print('depths=',depths) # check this is populated, lists how far each node is from substation\n",
    "# #print('depths length=',len(depths))\n",
    "\n",
    "# # print list of first 10 buses in network\n",
    "# count = 0 \n",
    "# for i in fin_feeder.network:    \n",
    "#     print(i) \n",
    "#     count += 1\n",
    "#     if count >= 10:\n",
    "#         break\n",
    "    \n",
    "# import csv  # jaimie do you need this?\n",
    "# graph = fin_feeder.network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for converting a .dot to .png\n",
    "(graph,) = pydot.graph_from_dot_file('generated_figs/good_bad_branch/badbranch_PBC')\n",
    "graph.write_png('generated_figs/good_bad_branch/badbranch_PBC_cb.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write busnames into a csv\n",
    "import csv\n",
    "graphNodes_noSub=hm.remove_subst_nodes(fin_feeder, file_name) # remove substation nodes, busList will have as many entries as R and X matrix length\n",
    "assert(len(graphNodes_noSub)*3==len(R)) # *3 is because R is 3ph\n",
    "with open(\"123NF_busList.csv\", 'w', newline='') as csvfile:\n",
    "    spamwriter = csv.writer(csvfile, delimiter='-',\n",
    "                            quotechar='|', quoting=csv.QUOTE_MINIMAL)\n",
    "    spamwriter.writerows(graphNodes_noSub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(node_index_map)\n",
    "R.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print 3-phase node_index table\n",
    "my_buses=list(node_index_map.keys())\n",
    "my_list=list(node_index_map.values())\n",
    "list1 = [(i+1) * 3-2 for i in my_list]\n",
    "list2 = [(i+1) * 3-1 for i in my_list]\n",
    "list3 = [(i+1) * 3-0 for i in my_list]\n",
    "for my_buses,list1, list2,list3 in zip(my_buses,list1,list2,list3):\n",
    "    print(my_buses,list1, list2,list3)\n",
    "#Table format: [bus_name, phaseA idx, phaseB idx, phaseC idx]\n",
    "# indices are shifted by 1 to be for MATLAB (starts at 1 not 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# R is of size 387, which is 129*3, so set n=129\n",
    "n=129\n",
    "parmObj=hm.configParms()\n",
    "foo=['PBC','PBC','PBC','PBC']\n",
    "parmObj.set_ctrlTypes(foo)\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "A, B = hm.setupStateSpace(parmObj,fin_feeder, n,node_index_map,depths,file_name)\n",
    "print('R = '+str(R))\n",
    "print('X = '+str(X))\n",
    "print('A = '+str(A))\n",
    "print('B = '+str(B))\n",
    "print(A.shape)\n",
    "print(B.shape)\n",
    "\n",
    "# write A and B matrices to csv\n",
    "#np.savetxt(\"123NF_Amat.csv\", A, delimiter=\",\")\n",
    "#np.savetxt(\"123NF_Bmat.csv\", B, delimiter=\",\")\n",
    "\n",
    "# Save R and X matrices to csv to import into matlab\n",
    "# np.savetxt reference: https://thispointer.com/how-to-save-numpy-array-to-a-csv-file-using-numpy-savetxt-in-python/\n",
    "np.savetxt('Rmat_123NF_accur.csv', R, delimiter=',')\n",
    "np.savetxt('Xmat_123NF_accur.csv', X, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [Optional] run impedance-related functions\n",
    "slack_bus = None\n",
    "for bus_name, depth in depths.items():\n",
    "    if depth == 0:\n",
    "        slack_bus = bus_name\n",
    "        break\n",
    "        \n",
    "# -------------------- now we call functions: ---------------------------------------------\n",
    "#print(depths)\n",
    "# modify node names when change feeders\n",
    "#plot_histogram_RX_ratios(fin_feeder, leaves_only = True)\n",
    "print('Z between buses:')\n",
    "print(np.around(imp.get_total_impedance_between_two_buses(fin_feeder, 'bus_37', 'bus_15',depths),2))\n",
    "print('\\nBus 49 Z to substation:')\n",
    "print(np.around(imp.get_total_impedance_from_substation(fin_feeder, 'bus_49',depths),2))\n",
    "print('\\nR/X ratio of bus 49 to substation:')\n",
    "print(imp.get_RX_ratio_tosubst(fin_feeder,'bus_49',depths))\n",
    "print('R/X Ratios:')\n",
    "print('\\n67 to 79: ')\n",
    "print(imp.get_RX_ratio_between_two_buses(fin_feeder, 'bus_67', 'bus_79',depths))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print z_to_sub for all buses to a csv\n",
    "graphNodes_nosub = hm.remove_subst_nodes(fin_feeder, file_name) # dont consider substation nodes, node 650 and 651 for 13NF\n",
    "\n",
    "myabsZs=np.array([], dtype=np.int64).reshape(0,3)\n",
    "\n",
    "for busName in graphNodes_nosub: # for all non-slack buses\n",
    "    z2sub=np.around(imp.get_total_impedance_from_substation(fin_feeder, busName,depths),2)\n",
    "    zself=[[z2sub[0][0], z2sub[1][1], z2sub[2][2]]]\n",
    "    print('zzelf=',zself)\n",
    "    abs_z=np.absolute(zself)\n",
    "    print('abs_z=',abs_z)\n",
    "    myabsZs=np.append(myabsZs,np.array(abs_z),axis=0) \n",
    "\n",
    "print(myabsZs)\n",
    "#print(np.array(myabsZs))\n",
    "np.savetxt('z2sub2.csv', np.array(myabsZs), delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[4.4  4.48 4.79]]\n",
      "[[4.6  4.36 4.41]]\n"
     ]
    }
   ],
   "source": [
    "# compute sum of z2subs for a given config\n",
    "cfg1=['bus_300','bus_197','bus_104','bus_108','bus_106'] # 2 redorange, 3 orange\n",
    "z2sub_cfg1=imp.get_z2sub_for_config(fin_feeder,cfg1,depths)\n",
    "print(np.around(z2sub_cfg1,2))\n",
    "\n",
    "cfg2=['bus_76','bus_49','bus_109','bus_251','bus_66'] # \n",
    "z2sub_cfg2=imp.get_z2sub_for_config(fin_feeder,cfg2,depths)\n",
    "print(np.around(z2sub_cfg2,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Color Bin Key for z2sub_heatmap\n",
      "darkgreen: (0.08845934148212913 --> 0.23706737443536646)\n",
      "limegreen: (0.23706737443536646 --> 0.38567540738860384)\n",
      "yellowgreen: (0.38567540738860384 --> 0.5342834403418412)\n",
      "yellow: (0.5342834403418412 --> 0.6828914732950785)\n",
      "gold: (0.6828914732950785 --> 0.8314995062483158)\n",
      "orange: (0.8314995062483158 --> 0.9801075392015532)\n",
      "tomato: (0.9801075392015532 --> 1.1287155721547906)\n",
      "firebrick: (1.1287155721547906 --> 1.277323605108028)\n"
     ]
    }
   ],
   "source": [
    "# create colormap of z2sub\n",
    "z2sub_dict={} # will populate\n",
    "graphNodes_nosub = hm.remove_subst_nodes(fin_feeder, file_name) # dont consider substation nodes, node 650 and 651 for 13NF\n",
    "for node_name in graphNodes_nosub:\n",
    "    z2sub=imp.get_total_impedance_from_substation(fin_feeder, node_name,depths)\n",
    "    foo=list((imp.get_RX_ratio_tosubst(fin_feeder, node_name, depths)).values())\n",
    "    zself=[[z2sub[0][0], z2sub[1][1], z2sub[2][2]]] #take the self impedances\n",
    "    zmag_per_node=max(np.absolute(zself)) # max |z| across phases, max to handle single-phase nodes\n",
    "    z2sub_dict[node_name]=zmag_per_node\n",
    "    \n",
    "vis.createColorMap(fin_feeder, z2sub_dict, 'z2sub_heatmap') # input values are scalar reals, one for each node"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------- THIS ONE ------------------\n",
    "# Rx ratio and phase coupling ratio informs bad branches\n",
    "branches = vis.assign_network_branches(fin_feeder, substation_name)\n",
    "for branch in branches:\n",
    "    if 'bus_150'in branch:\n",
    "        branch.remove('bus_150')\n",
    "        branch.remove('bus_149')\n",
    "rx_dic = {}\n",
    "phase_ratios = vis.phaseCouplingPerNode(fin_feeder, depths, file_name)\n",
    "graphNodes_nosub = hm.remove_subst_nodes(fin_feeder, file_name) # dont consider substation nodes, node 650 and 651 for 13NF\n",
    "for node_name in graphNodes_nosub:\n",
    "    foo=list((imp.get_RX_ratio_tosubst(fin_feeder, node_name, depths)).values())\n",
    "    rx_3ph=np.array(foo).astype('float')\n",
    "    rx_3ph[rx_3ph==0]=np.nan # assign NaN to all zeros\n",
    "    rx_dic[node_name]=np.nanmean(rx_3ph) # nanmean isgnores nans, now scalar metric for each node\n",
    "    \n",
    "    ph_3ph=phase_ratios[node_name].astype('float')\n",
    "    ph_3ph[ph_3ph==0]=np.nan # assign NaN to all zeros\n",
    "    # take mean across phases, then take mag of the complex scalar\n",
    "    phase_ratios[node_name]=np.absolute(np.nanmean(ph_3ph)) # nanmean isgnores nans\n",
    "    \n",
    "vis.createColorMap(fin_feeder, phase_ratios, 'phase_coupling') # input values are scalar reals, one for each node\n",
    "vis.createColorMap(fin_feeder, rx_dic, 'PQ_coupling')\n",
    "\n",
    "#------------ make mixed colormap -----------\n",
    "normd_phase_ratios = {}\n",
    "normd_rx = {}\n",
    "avg_phase = np.nanmean(list(phase_ratios.values()))\n",
    "avg_rx = np.nanmean(list(rx_dic.values()))\n",
    "\n",
    "for key, val in phase_ratios.items():\n",
    "    if not np.isnan(val):\n",
    "        normd_phase_ratios[key] = val/avg_phase\n",
    "for key, val in rx_dic.items():\n",
    "    if not np.isnan(val):\n",
    "        normd_rx[key] = val/avg_rx\n",
    "\n",
    "mixed_vals = {}\n",
    "for node in graphNodes_nosub:\n",
    "    mixed_vals[node] = 1.7*normd_phase_ratios[node] + 0.3*normd_rx[node]\n",
    "vis.createColorMap(fin_feeder, mixed_vals, 'mixed')\n",
    "# colormaps show up in generated_figs folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collect cr,cx,and d for each branch --> put into coupling paper table\n",
    "# this code doenst account for single and double phase, averages across the zeros in those cases\n",
    "\n",
    "node1='bus_68'\n",
    "node2='bus_71'\n",
    "\n",
    "ph_3ph=imp.get_phase_ratio_between_two_buses(fin_feeder, node1,node2, depths)\n",
    "print('complex 3ph=',ph_3ph)\n",
    "cr=np.real(ph_3ph)\n",
    "cx=np.imag(ph_3ph)\n",
    "phrat_rscalar=np.absolute(np.nanmean(ph_3ph)) # real scalar metric\n",
    "print('real scalar=',phrat_rscalar)\n",
    "\n",
    "rx_3ph=imp.get_RX_ratio_between_two_buses(fin_feeder, node1,node2, depths)\n",
    "print(rx_3ph.values())\n",
    "\n",
    "z_brhead_subst=imp.get_total_impedance_from_substation(fin_feeder, node1, depths)\n",
    "import scipy.linalg as LA\n",
    "u, s, vh=LA.svd(z_brhead_subst)\n",
    "len2sub=s[0] # largest singular value is length metric for matrix\n",
    "\n",
    "print('------------------')\n",
    "print('d=',np.around(sum(rx_3ph.values())/len(rx_3ph),3))\n",
    "print('cr=',np.around(np.nanmean(cr),3))\n",
    "print('cx=',np.around(np.nanmean(cx),3))\n",
    "print('length=',np.around(len2sub,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------- make table of branches, ratios, and num used ------------\n",
    "#3/17: this code chunk not used yet, made table in coupling figs spreadsheet instead\n",
    "\n",
    "# fix: make list of branches only branches with 3 phase lines, since single and double phase laterals wont be used much \n",
    "\n",
    "# print('branches table:') # format into a table with index and branch nodes\n",
    "# i=0\n",
    "# for br in branches:\n",
    "#     print(i,',',br)\n",
    "#     i+=1\n",
    "    \n",
    "num_used=np.full((1, len(graphNodes_nosub)),-1) # placeholder\n",
    "\n",
    "print('bus | rx_rat | ph_rat | num_used')\n",
    "print('large and small rx_rat:')\n",
    "foo1=sorted(normd_rx, key=normd_rx.get, reverse=True)[:3] # get keys of 3 largest values of dict\n",
    "for node in foo1:\n",
    "    print(node,np.around(normd_rx[node],3),np.around(normd_phase_ratios[node],3),-1)\n",
    "\n",
    "foo2=sorted(normd_rx, key=normd_rx.get, reverse=False)[:3] # get keys of 3 smallest values of dict\n",
    "for node in foo2:\n",
    "    print(node,np.around(normd_rx[node],3),np.around(normd_phase_ratios[node],3),-1)\n",
    "\n",
    "print('large and small ph_rat:')\n",
    "foo3=sorted(normd_phase_ratios, key=normd_phase_ratios.get, reverse=True)[:3] # get keys of 5 largest values of dict\n",
    "for node in foo3:\n",
    "    print(node,np.around(normd_rx[node],3),np.around(normd_phase_ratios[node],3),-1)\n",
    "\n",
    "foo4=sorted(normd_phase_ratios, key=normd_phase_ratios.get, reverse=False)[:3] # get keys of 5 smallest values of dict\n",
    "for node in foo4:\n",
    "    print(node,np.around(normd_rx[node],3),np.around(normd_phase_ratios[node],3),-1)\n",
    "\n",
    "print('large and small sum_rat:')\n",
    "foo3=sorted(mixed_vals, key=mixed_vals.get, reverse=True)[:3] # get keys of 5 largest values of dict\n",
    "for node in foo3:\n",
    "    print(node,np.around(normd_rx[node],3),np.around(normd_phase_ratios[node],3),-1)\n",
    "\n",
    "foo4=sorted(mixed_vals, key=mixed_vals.get, reverse=False)[:3] # get keys of 5 smallest values of dict\n",
    "for node in foo4:\n",
    "    print(node,np.around(normd_rx[node],3),np.around(normd_phase_ratios[node],3),-1)\n",
    "\n",
    "    \n",
    "# for node in graphNodes_nosub:\n",
    "#     print(node,np.around(normd_phase_ratios[node],3),np.around(normd_rx[node],3),-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Source.from_file('colorMap_rx') # display colormap to see good/bad branches \n",
    "Source.from_file('colorMap_phase_coupl') # display colormap to see good/bad branches\n",
    "Source.from_file('colorMap_mixed') # display colormap to see good/bad branches  "
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}