Driver-123NF_processes.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The jupyternotify extension is already loaded. To reload it, use:\n",
      "  %reload_ext jupyternotify\n"
     ]
    }
   ],
   "source": [
    "%load_ext jupyternotify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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\n",
    "import my_process_funcs as prc # new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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_accur.xlsx\"\n",
    "#modelpath = filepath + \"004_GB_IEEE123_OPAL_zeroCrossZ.xlsx\"\n",
    "\n",
    "    # All GridBright load files should be in the following folder\n",
    "    #loadfolder = \"/Users/jasperpakshong/Documents/Berkeley/ENERGISE/IEEE13/\"\n",
    "    #loadpath = loadfolder + \"IEEE13testload_w_extreme_act.xlsx\"\n",
    "\n",
    "    # filepath = \"IEEE13/\"\n",
    "    # modelpath = filepath + \"001 phasor08_IEEE13_OPAL.xls\"\n",
    "    # loadfolder = \"IEEE13/\"\n",
    "    # loadpath = loadfolder + \"001_phasor08_IEEE13_norm03_HIL_7_1.xlsx\"\n",
    "\n",
    "    #filepath = \"AL0001/\"\n",
    "    #modelpath = filepath + \"AL0001_OPAL_working.xls\"\n",
    "    #loadfolder = \"AL0001/\"\n",
    "    #loadpath = loadfolder + \"AL0001_tvload_afternoon1h.csv\"\n",
    "      \n",
    "    #filepath = \"13NF_balanced/\"\n",
    "    #modelpath = filepath + \"016 GB_IEEE13_balance_all_ver2.xls\"\n",
    "    #loadfolder = \"13NF_balanced/\"\n",
    "    #loadpath = loadfolder + \"016 GB_IEEE13_balance all ver2_time_sigBuilder_secondWise_norm03.csv\"\n",
    "\n",
    "#==========================================================================================================\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [ESSENTIAL] create feeder object\n",
    "\n",
    "fin_feeder = ff.feeder_init(modelpath,'','',timesteps,Vbase_ll,Sbase,depths,leaves)\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",
    "#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",
    "print(\"First ten nodes are:\")\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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%notify\n",
    "\n",
    "# TEMP\n",
    "# --------- run eval_config on CPP config --------------\n",
    "# need check that end-config of CPP process is feasible for both control frameworks\n",
    "parmObj=hm.configParms()\n",
    "# parmObj.set_ctrlTypes(['VVC', 'VVC', 'VWC', 'VWC', 'VWC', 'VVC', 'VVC', 'VVC'])\n",
    "# parmObj.set_version(2) # 1 for PBC\n",
    "\n",
    "#parmObj.set_ctrlTypes(['PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC'])\n",
    "parmObj.set_ctrlTypes(['PBC', 'PBC'])\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "\n",
    "#all_act_locs = ['bus_152', 'bus_21', 'bus_58', 'bus_70', 'bus_4', 'bus_82', 'bus_109', 'bus_47']\n",
    "# all_act_locs = ['bus_47','bus_58','bus_4','bus_109','bus_21','bus_70','bus_152','bus_82']\n",
    "# perf_nodes=['bus_47','bus_58','bus_4','bus_109','bus_21','bus_70','bus_152','bus_82']\n",
    "\n",
    "all_act_locs = ['bus_81','bus_44']\n",
    "perf_nodes=['bus_81','bus_44']\n",
    "# print(len(all_act_locs))\n",
    "# print(len(foo))\n",
    "feas, maxError,numfeas,bestF,indicMat =prc.eval_config(parmObj,fin_feeder, all_act_locs,perf_nodes, node_index_map,substation_name,depths,file_name,Vbase_ll, Sbase)\n",
    "\n",
    "\n",
    "# we have F as 762x762, and number of 1-eigs is 750, so 12 less"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------- TEMP ----------\n",
    "## ------------ run placeMaxColocActs_stopAtInfeas, aka overloadCPP process -----------\n",
    "# place colocated actuators until an infeasible loc is tested, then call find_good_colocated and return \n",
    "# produces colormap associated with \"motex\" motivating example\n",
    "\n",
    "t = time.time()\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "max_act_config_stopInfeas,parmObj=prc.placeMaxColocActs_stopAtInfeas(parmObj,fin_feeder, file_name, node_index_map, depths, substation_name,Vbase_ll, Sbase)\n",
    "elapsed = time.time() - t\n",
    "print('Time elapsed=',elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------- Optional -----------------\n",
    "# Run if want to diaply generated .png files in the python notebook\n",
    "from IPython.display import display\n",
    "from PIL import Image\n",
    "\n",
    "# shorten png printing commands into one function call\n",
    "def printpng(filename):\n",
    "    im = Image.open(filename) \n",
    "    #print(im.size) # (pixel width,pixel height)\n",
    "    width=im.size[0]\n",
    "    height=im.size[1]\n",
    "    im.draft('RGB',(width/1,height/1)) \n",
    "    px = im.load()\n",
    "    display(im)\n",
    "    # reference on img.draft func: https://stackoverflow.com/questions/57663734/how-to-speed-up-image-loading-in-pillow-python\n",
    "    \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------ run find_good_colocated -----------------\n",
    "# assessed feas of placing 1 co-located actuator at each loc and makes colormap\n",
    "\n",
    "# all_act_locs = []\n",
    "# perf_nodes = []\n",
    "# here we set act_locs=[] to create heatmap for empty network\n",
    "# feas_configs, heatMapNames=prc.find_good_colocated(fin_feeder, [], node_index_map, substation_name, depths,file_name, Vbase_ll, Sbase)\n",
    "\n",
    "# print('----------- Feas configs are: ---------')\n",
    "# print(feas_configs)\n",
    "# Source.from_file(heatMapNames[0]) # display graph in notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create figure 2, colormap of phase coupling \n",
    "ratios = vis.phaseCouplingPerNode(fin_feeder,depths,file_name)\n",
    "for key, val in ratios.items():\n",
    "    ratios[key] = np.mean(val) # take the mean of phase coupling across all three phases\n",
    "    \n",
    "vis.createColorMap(fin_feeder, ratios, 'fig2_123NF')\n",
    "ff.clear_graph(fin_feeder)\n",
    "Source.from_file('colorMap_fig2_123NFwide') # display graph in notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- Run computeFParamSpace, for debugging ---------\n",
    "# all_act_locs = ['bus_25','bus_39','bus_56','bus_108'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_25','bus_39','bus_56','bus_108']\n",
    "\n",
    "# Fq_ub,Fp_ub=ctrl.computeFParamSpace_v2(fin_feeder, all_act_locs, perf_nodes,R,X,depths,node_index_map)\n",
    "# # this function prints Zgood (from R and X matrices) and Z_toSubst, where Z_toSubst is correct but doesnt match Zgood. Investigate\n",
    "# print('(Fp_ub,Fq_ub)=(',Fp_ub,',',Fq_ub,')')\n",
    "\n",
    "#vis.markActuatorConfig(all_act_locs, fin_feeder, file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------- Config Chunk (Assign act and perf locs here) ----------\n",
    "\n",
    "# binary heat map that shows good locations to place a new actuator when\n",
    "# a) there are 0 existing actuators, b) there is 1 existing actuator, and c) there are two existing actuators.\n",
    "# actuator list = complex config 1, i.e. at each step we choose one of the good places and buildup to complex config 1\n",
    "\n",
    "# #complex config 1 (feas)\n",
    "#all_act_locs = ['bus_49','bus_76','bus_65'] # 49, then 76, then 65\n",
    "#perf_nodes = ['bus_49','bus_76','bus_65']\n",
    "# 16336 seconds \n",
    "\n",
    "#complex config 2 (feas)\n",
    "# all_act_locs = ['bus_25','bus_39','bus_56','bus_108'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_25','bus_39','bus_56','bus_108']\n",
    "# 30206 seconds\n",
    "\n",
    "# #complex config 3 (feas)\n",
    "# all_act_locs = ['bus_105','bus_112','bus_100','bus_91'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_105','bus_112','bus_100','bus_91']\n",
    "# # time taken = 21674 seconds\n",
    "\n",
    "# #complex config 4, many act 1 perf\n",
    "#all_act_locs = ['bus_42','bus_42','bus_42','bus_42'] # 49, then 76, then 65\n",
    "#perf_nodes = ['bus_151','bus_46','bus_251','bus_8']\n",
    "# # time taken = XXX seconds\n",
    "\n",
    "# # complex config 3\n",
    "# all_act_locs = ['bus_105','bus_112','bus_100','bus_91'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_105','bus_112','bus_100','bus_91']\n",
    "\n",
    "# evaluate this config\n",
    "# all_act_locs = ['bus_60','bus_58','bus_20','bus_72'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_300','bus_58','bus_20','bus_72']\n",
    "\n",
    "# Adam sim\n",
    "# all_act_locs = ['bus_49','bus_76','bus_152'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_49','bus_76','bus_65']\n",
    "\n",
    "# # eval neighborhood config\n",
    "# all_act_locs = ['bus_82','bus_87','bus_76','bus_49','bus_41','bus_46','bus_60'] # 49, then 76, then 65\n",
    "# perf_nodes = ['bus_77','bus_77','bus_77','bus_44','bus_44','bus_44','bus_66']\n",
    "\n",
    "# run RHP on neighborhood config for 7th step only:\n",
    "# set_act_locs = ['bus_82','bus_87','bus_76','bus_49','bus_41','bus_46']\n",
    "# set_perf = ['bus_77','bus_77','bus_77','bus_44','bus_44','bus_44']\n",
    "# addon_act_nodes = ['bus_60']\n",
    "# addon_perf_nodes = ['bus_66']\n",
    "\n",
    "# heatmap/collect stable configs for steps 1-7\n",
    "# set_act_locs = []\n",
    "# set_perf = []\n",
    "# addon_act_nodes = ['bus_82','bus_87','bus_76','bus_49','bus_41','bus_46','bus_60']\n",
    "# addon_perf_nodes = ['bus_77','bus_77','bus_77','bus_44','bus_44','bus_44','bus_66']\n",
    "\n",
    "# # run RHP for step 3 of 10-step placeMaxColoc, not going in paper\n",
    "# set_act_locs = ['bus_58','bus_20','bus_72']\n",
    "# set_perf = ['bus_58','bus_20','bus_72']\n",
    "# addon_act_nodes = ['bus_300']\n",
    "# addon_perf_nodes = ['bus_300']\n",
    "\n",
    "# motivating example (motex), list is equal to output 'max_act_config_stopInfeas'\n",
    "# set_act_locs = ['bus_55','bus_33','bus_5','bus_96','bus_113','bus_197']\n",
    "# set_perf = ['bus_55','bus_33','bus_5','bus_96','bus_113','bus_197']\n",
    "# addon_act_nodes = ['bus_60']\n",
    "# addon_perf_nodes = ['bus_60']\n",
    "\n",
    "\n",
    "# Auto CPP configs:\n",
    "# seed=1\n",
    "# set_act_locs = []\n",
    "# set_perf = []\n",
    "# addon_act_nodes = ['bus_114','bus_52','bus_85','bus_8','bus_75','bus_106','bus_22','bus_113','bus_4','bus_7','bus_4','bus_3','bus_2','bus_5']\n",
    "# addon_perf_nodes = ['bus_114','bus_52','bus_85','bus_8','bus_75','bus_106','bus_22','bus_113','bus_4','bus_7','bus_4','bus_3','bus_2','bus_5']\n",
    "\n",
    "# # seed=2\n",
    "# set_act_locs = []\n",
    "# set_perf = []\n",
    "# addon_act_nodes = ['bus_96','bus_85','bus_97','bus_83','bus_105','bus_11','bus_152','bus_197','bus_6','bus_102','bus_58','bus_4','bus_2']\n",
    "# addon_perf_nodes = ['bus_96','bus_85','bus_97','bus_83','bus_105','bus_11','bus_152','bus_197','bus_6','bus_102','bus_58','bus_4','bus_2']\n",
    "\n",
    "# seed=3\n",
    "# set_act_locs = []\n",
    "# set_perf = []\n",
    "# addon_act_nodes = ['bus_197','bus_55','bus_5','bus_113','bus_33','bus_96','bus_45','bus_63','bus_59','bus_106','bus_10','bus_91','bus_46','bus_25','bus_6','bus_17','bus_80','bus_49','bus_4','bus_73',]\n",
    "# addon_perf_nodes = ['bus_96','bus_85','bus_97','bus_83','bus_105','bus_11','bus_152','bus_197','bus_6','bus_102','bus_58','bus_4','bus_2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------- mark actuator config -------------\n",
    "# assign 'all_act_locs' in config chunk\n",
    "vis.markActuatorConfig(all_act_locs,fin_feeder, file_name)\n",
    "\n",
    "# convert the generated .dot to.png:\n",
    "# (graph,) = pydot.graph_from_dot_file('generated_figs/'+'actConfig_'+ file_name) # requires 'import pydot'\n",
    "# graph.write_png('generated_figs/'+'config_60.png') # can change the name\n",
    "# printpng('generated_figs/config_60.png') # can change the name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------- run eval_config on CPP config --------------\n",
    "# need check that end-config of CPP process is feasible for both control frameworks\n",
    "\n",
    "parmObj=hm.configParms()\n",
    "# parmObj.set_ctrlTypes(['VVC', 'VVC', 'VWC', 'VWC', 'VWC', 'VVC', 'VVC', 'VVC'])\n",
    "# parmObj.set_version(2) # 1 for PBC\n",
    "\n",
    "parmObj.set_ctrlTypes(['PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC'])\n",
    "parmObj.set_ctrlTypes(['PBC', 'PBC'])\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "\n",
    "#all_act_locs = ['bus_152', 'bus_21', 'bus_58', 'bus_70', 'bus_4', 'bus_82', 'bus_109', 'bus_47']\n",
    "all_act_locs = ['bus_47','bus_58','bus_4','bus_109','bus_21','bus_70','bus_152','bus_82']\n",
    "perf_nodes=['bus_47','bus_58','bus_4','bus_109','bus_21','bus_70','bus_152','bus_82']\n",
    "\n",
    "all_act_locs = ['bus_81','bus_44']\n",
    "perf_nodes=['bus_81','bus_44']\n",
    "# print(len(all_act_locs))\n",
    "# print(len(foo))\n",
    "feas, maxError,numfeas,bestF,indicMat=prc.eval_config(parmObj,fin_feeder, all_act_locs,perf_nodes, node_index_map,substation_name,depths,file_name,Vbase_ll, Sbase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -------- run eval_config for a config with actuators at ALL nodes of the network  -----------\n",
    "# from running this, hypothesis is: multiple actautors on every node tracking single perf is feas for ANY choice of perf\n",
    "\n",
    "# assign 'all_act_locs' and 'perf_nodes' in config chunk\n",
    "graphNodes_nosub = hm.remove_subst_nodes(fin_feeder, file_name) # dont consider co-located at substation nodes, node 650 and 651\n",
    "all_act_locs=graphNodes_nosub[:] # place actuators at ALL non-substation nodes\n",
    "n=len(all_act_locs)\n",
    "print('acts=',all_act_locs)\n",
    "\n",
    "perf='bus_84' # you choose\n",
    "perf_nodes=[perf] * n\n",
    "print('perfs=',perf)\n",
    "feas, maxError,numfeas,indicMat=prc.eval_config(parmObj,fin_feeder, all_act_locs,perf_nodes, node_index_map,substation_name,depths,file_name,Vbase_ll, Sbase)\n",
    "\n",
    "#printpng('generated_figs/heatmap_colocated_123NF.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- good vs. bad branch scenario, 1-step heatmap\n",
    "set_acts = ['bus_8','bus_53','bus_57','bus_66'] # has bad branch\n",
    "addon_acts= ['bus_152'] # arbitrary\n",
    "\n",
    "parmObj=hm.configParms()\n",
    "foo=['VVC','VWC','VVC','VWC']\n",
    "parmObj.set_ctrlTypes(foo)\n",
    "parmObj.set_version(2) # 1 for PBC\n",
    "\n",
    "# here we set act_locs as the existing actuators\n",
    "feas_configs, heatMapNames=prc.find_good_colocated(parmObj,fin_feeder, set_acts, addon_acts, node_index_map, substation_name, depths,file_name, Vbase_ll, Sbase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- good vs. bad branch scenario, 1-step heatmap\n",
    "set_acts = ['bus_8','bus_53','bus_57','bus_74'] # has good branch\n",
    "addon_acts= ['bus_152'] # arbitrary\n",
    "\n",
    "parmObj=hm.configParms()\n",
    "foo=['VVC','VWC','VVC','VWC']\n",
    "parmObj.set_ctrlTypes(foo)\n",
    "parmObj.set_version(2) # 1 for PBC\n",
    "\n",
    "# here we set act_locs as the existing actuators\n",
    "feas_configs, heatMapNames=prc.find_good_colocated(parmObj,fin_feeder, set_acts, addon_acts,node_index_map, substation_name, depths,file_name, Vbase_ll, Sbase)\n",
    "\n",
    "# result looks like: 'heatmap_colocated' + '_' + file_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# ---------------- run CPP (formely called find_good_coloated) ---------------------------\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_ctrlTypes(['VVC', 'VVC', 'VWC', 'VWC', 'VWC', 'VVC', 'VVC', 'VVC']) # format is ctrl types of [set_acts addon_acts]\n",
    "parmObj.set_version(2) # 1 for PBC\n",
    "\n",
    "# parmObj.set_ctrlTypes(['PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC', 'PBC'])\n",
    "# parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "# set_acts = []\n",
    "# addon_acts=['bus_152', 'bus_21', 'bus_58', 'bus_70', 'bus_4', 'bus_82', 'bus_109', 'bus_47']\n",
    "set_acts = []\n",
    "addon_acts = ['bus_47','bus_58','bus_4','bus_109','bus_21','bus_70','bus_152','bus_82']\n",
    "\n",
    "\n",
    "t = time.time()\n",
    "lst_feas_configs, lzn_error_run_sum, heatMapNames = prc.find_good_colocated(parmObj,fin_feeder, set_acts, addon_acts, node_index_map,substation_name, depths, file_name, Vbase_ll, Sbase)\n",
    "# lst_feas_configs is list of dictionaries\n",
    "elapsed = time.time() - t\n",
    "print('Time elapsed=',elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## ------------ run placeMaxColocActs_stopAtInfeas, aka overloadCPP process -----------\n",
    "# place colocated actuators until an infeasible loc is tested, then call find_good_colocated and return \n",
    "# produces colormap associated with \"motex\" motivating example\n",
    "\n",
    "t = time.time()\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "max_act_config_stopInfeas,parmObj=prc.placeMaxColocActs_stopAtInfeas(parmObj,fin_feeder, file_name, node_index_map, depths, substation_name,Vbase_ll, Sbase)\n",
    "elapsed = time.time() - t\n",
    "print('Time elapsed=',elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------ run place_max_colocated_acts, aka auto_overloadCPP process -------------\n",
    "t = time.time()\n",
    "seed=2 # choose\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_version(2) # 1 for PBC\n",
    "print('----------- Seed='+str(seed)+' ----------------------------')\n",
    "max_act_config,parmObj = prc.place_max_coloc_acts(parmObj,seed,fin_feeder, file_name, node_index_map, depths, substation_name,Vbase_ll, Sbase)\n",
    "elapsed = time.time() - t\n",
    "print('Time elapsed=',elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%notify\n",
    "# ------- run runHeatMapProcess (RHP) ---------------\n",
    "# 6.18.22\n",
    "# run RHP on neighborhood config for 1st to 7th step:\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_ctrlTypes(['PBC','PBC','PBC','PBC','PBC','PBC','PBC'])\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "# set_act_locs = []\n",
    "# set_perf = []\n",
    "# addon_act_nodes = ['bus_82','bus_87','bus_76','bus_49','bus_41','bus_46','bus_60']\n",
    "# addon_perf_nodes = ['bus_77','bus_77','bus_77','bus_44','bus_44','bus_44','bus_66']\n",
    "\n",
    "# run RHP on neighborhood config for 7th step only:\n",
    "set_act_locs = ['bus_82','bus_87','bus_76','bus_49','bus_41','bus_46']\n",
    "set_perf = ['bus_77','bus_77','bus_77','bus_44','bus_44','bus_44']\n",
    "addon_act_nodes = ['bus_60']\n",
    "addon_perf_nodes = ['bus_66']\n",
    "\n",
    "t = time.time()\n",
    "lst_feas_configs, lzn_error_run_sum, heatMapNames = prc.runHeatMapProcess(parmObj,fin_feeder, set_act_locs, set_perf, addon_act_nodes, addon_perf_nodes, node_index_map, substation_name, depths, file_name,Vbase_ll, Sbase)\n",
    "# lst_feas_configs is list of dictionaries\n",
    "elapsed = time.time() - t\n",
    "print('Time elapsed=',elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Save RHP variable 'lst_feas_configs' to .pkl file, so dont have to run again in future\n",
    "import pickle\n",
    "filename='CPPseed2_lstfeasconfigs_droop_mod1.6.pkl'\n",
    "with open(filename, 'wb') as f:  # Python 3: open(..., 'wb')\n",
    "    pickle.dump(lst_feas_configs, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------- Setup branch analysis on 'lst_feas_configs' populated by RHP -------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the data 'lst_feas_configs' from .pkl file\n",
    "# import pickle\n",
    "# filename='nbhs_RHP_lstfeasconfigs.pkl'\n",
    "# with open(filename,'rb') as f:  # Python 3: open(..., 'rb')\n",
    "#     lst_feas_configs = pickle.load(f)\n",
    "     \n",
    "    \n",
    "import os\n",
    "import pickle\n",
    "scores = {} # scores is an empty dict already\n",
    "filename='CPPseed2_lstfeasconfigs_PBC_mod1.1.pkl'\n",
    "print(os.path.getsize(filename))\n",
    "if os.path.getsize(filename) > 0:      \n",
    "    with open(filename, \"rb\") as f:\n",
    "        unpickler = pickle.Unpickler(f)\n",
    "        # if file is not empty scores will be equal\n",
    "        # to the value unpickled\n",
    "        lst_feas_configs = unpickler.load()\n",
    "        \n",
    "print(len(lst_feas_configs)) # just to see if it loaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to save and restore variables..\n",
    "# https://stackoverflow.com/questions/2960864/how-can-i-save-all-the-variables-in-the-current-python-session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# ------------ create histogram of numfeas --------------\n",
    "numFeasList=[]\n",
    "#print([d['numfeas'] for d in lst_feas_configs])\n",
    "for cfg in lst_feas_configs:\n",
    "    a=cfg['numfeas']\n",
    "    b=a[0].item() # remove layers of brackets\n",
    "    numFeasList.append(b)\n",
    "\n",
    "#print(numFeasList)\n",
    "plt.title(\"numFeas Histogram\")\n",
    "plt.xlabel(\"number of feasible Fs\")\n",
    "plt.ylabel(\"number of configs\")\n",
    "plt.hist(numFeasList, bins=40) # -10: is last 10 ele\n",
    "#plt.hist(numFeasList[-10:], bins=40) # -10: is last 10 ele\n",
    "plt.savefig('numFeasHist.png') # modify for each feeder\n",
    "plt.show() # need to savefig before plt.show\n",
    "\n",
    "print(len(numFeasList))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert lst_feas_configs (list of dictionaries) to newlst_feas_configs (a list of lists)\n",
    "# input a list of feasible actuator configurations where each configuration is its own list (ie a list of lists)\"\n",
    "length=len(lst_feas_configs)\n",
    "newlst_feas_configs = []\n",
    "numact_lst=[] # store number of actuators in separate list\n",
    "i=0\n",
    "for dic in lst_feas_configs:\n",
    "    #print('dic=',dic)\n",
    "    #print(type(dic))\n",
    "    i+=1\n",
    "    if i<length-1:\n",
    "        #print('dic[act]=',str(dic['act']))\n",
    "        newlst_feas_configs += [dic['act']]\n",
    "        numact_lst+=[dic['numfeas']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create 'placed' variable\n",
    "# placed is an ordered list of lists, where the k'th ele of the list is the actuators placed at that step\n",
    "set_acts=[]\n",
    "addon_acts=['bus_152', 'bus_21', 'bus_58', 'bus_70', 'bus_4', 'bus_82', 'bus_109', 'bus_47']\n",
    "acts=set_acts + addon_acts\n",
    "\n",
    "placed=[]\n",
    "k=0\n",
    "while k<=len(acts):\n",
    "    placed+=[acts[0:k]]\n",
    "    k+=1\n",
    "print(placed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- run plot_actuator_num_histogram ---------\"\n",
    "# input list of lists of feasible actuator configurations\"\n",
    "# output histogram reflecting number of actuators in each feas configuration\"\n",
    "vis.plot_actuator_num_histogram(newlst_feas_configs, file_name)  \n",
    "# note the max number in the histogram\n",
    "\n",
    "# ----- run assign_network_branches ---------\"\n",
    "branches = vis.assign_network_branches(fin_feeder, substation_name)\n",
    "vis.mark_network_branches(fin_feeder, branches, file_name, substation_name,depths)\n",
    "#Source.from_file('branch_key:' + file_name)\n",
    "print('branches table:') # format into a table with index and branch nodes\n",
    "i=0\n",
    "for br in branches:\n",
    "    print('idx',i,': ',br)\n",
    "    #print('phase',vis.getPhases(br,fin_feeder))\n",
    "    i+=1\n",
    "\n",
    "\n",
    "# ----- run find_good_branches ---------\"\n",
    "# the 'best' branch is the branch with the most actuators from feas configs on it across all configurations\n",
    "num_good_bad_branches = 4 # return the n best and worst branches\n",
    "vis.find_good_branches(newlst_feas_configs, branches, num_good_bad_branches, placed)\n",
    "\n",
    "# ----- run determine_good_or_bad_branch ---------\"\n",
    "#print('\\nBranches are ', branches)\n",
    "branch = branches[-2]  # you choose as input\n",
    "print('\\nBranch chosen=',branch)\n",
    "good_branch_threshold = 3 # threshhold # configs that must use the branch (i.e. an act on that branch) for the branch to be \"good\"\n",
    "configs_with_branch = vis.determine_good_or_bad_branch(branch, newlst_feas_configs, good_branch_threshold)\n",
    "# output: configs that use the branch, and if that # exceeds threshhold, call the branch \"good\"\n",
    "print('Num feas configs =', len(newlst_feas_configs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -------- Example experiment for Table 1; phase ratio mod2x, PBC, m=5 -----------\n",
    "# Desc: find an $m$-node configuration under PBC or droop control that is \\emph{good} on the orig 123NF \n",
    "# and record the best associated $F$. Then we apply this best $F$ to the same configuration \n",
    "# on a modified 123NF, where every line has a higher R/X or phase ratio with the same line length\n",
    "#If the configuration on the higher-ratio 123NF is no longer \\emph{good}, the result supports the remarks. \n",
    "# Then we do the reverse.\n",
    "\n",
    "i=0\n",
    "support=0 # count how many times inc feeder ratio loses stability for same config/kgain set\n",
    "not_support=0\n",
    "while i<15: # try design config for a number of trials\n",
    "    print('-----------------------------')\n",
    "    \n",
    "    # --------- Initialize orig feeder -------------\n",
    "    modelpath = filepath + \"004_GB_IEEE123_OPAL.xls\"\n",
    "    file_name = '123NForig'\n",
    "    fin_feeder = ff.feeder_init(modelpath,loadfolder,loadpath,timesteps,Vbase_ll,Sbase,depths,leaves)\n",
    "    node_index_map = hm.createNodeIndexMap(fin_feeder) #node indices for indicMat and F matrix\n",
    "\n",
    "    parmObj=hm.configParms()\n",
    "    parmObj.set_version(1) # 1 for PBC\n",
    "    numAct=5\n",
    "    seedkey=i # different random number seed for each trial, so get different sets of random configs\n",
    "    parmObj,act_locs,bestFparm,indicMat=prc.design_config(parmObj,seedkey,numAct,fin_feeder, file_name, node_index_map, depths, substation_name,Vbase_ll, Sbase)\n",
    "    \n",
    "    if bestFparm is not -1: # proceed if designing config found a config\n",
    "        # --------- Initialize mod feeder -------------\n",
    "        modelpath = filepath + \"123NF_phrat_mod2.xls\"\n",
    "    #     modelpath = filepath + \"123NF_phrat_mod1.5.xls\"\n",
    "    #     modelpath = filepath + \"123NF_rxrat_mod2.xls\"\n",
    "    #     modelpath = filepath + \"123NF_rxrat_mod1.5.xls\"\n",
    "        file_name = '123NFmod2'\n",
    "        fin_feeder = ff.feeder_init(modelpath,loadfolder,loadpath,timesteps,Vbase_ll,Sbase,depths,leaves)\n",
    "        node_index_map = hm.createNodeIndexMap(fin_feeder) #node indices for indicMat and F matrix\n",
    "\n",
    "        feas=ctrl.eval_config_F(parmObj,bestFparm,indicMat,fin_feeder,depths,node_index_map)\n",
    "        if feas==False: # if bestF on mod feeder is infeas, break\n",
    "            print('-----> Support! F on modded feeder is infeas')\n",
    "            support+=1\n",
    "        else:\n",
    "            print('-----> inconclusive, F on modded feeder is still feas')\n",
    "\n",
    "        \n",
    "    # now do the reverse...\n",
    "    parmObj,act_locs,bestFparm,indicMat=prc.design_config(parmObj,seedkey,numAct,fin_feeder, file_name, node_index_map, depths, substation_name,Vbase_ll, Sbase)\n",
    "    \n",
    "    if bestFparm is not -1: # proceed if designing config found a config\n",
    "        print(act_locs)  # print config found\n",
    "        print(parmObj.get_ctrlTypes())\n",
    "\n",
    "        # --------- Initialize orig feeder -------------\n",
    "        modelpath = filepath + \"004_GB_IEEE123_OPAL.xls\"\n",
    "        file_name = '123NForig'\n",
    "        fin_feeder = ff.feeder_init(modelpath,loadfolder,loadpath,timesteps,Vbase_ll,Sbase,depths,leaves)\n",
    "        node_index_map = hm.createNodeIndexMap(fin_feeder) #node indices for indicMat and F matrix\n",
    "\n",
    "        feas=ctrl.eval_config_F(parmObj,bestFparm,indicMat,fin_feeder,depths,node_index_map)\n",
    "        if feas==False: # if bestF on mod feeder is infeas, break\n",
    "            print('-----> Not support - F on orig feeder is infeas')\n",
    "            not_support+=1\n",
    "        else:\n",
    "            print('-----> inconclusive, F on orig feeder is still feas')\n",
    "    \n",
    "    if support+not_support>=10:\n",
    "        break\n",
    "        \n",
    "    i+=1\n",
    "#----------------------------------------------------------------------------\n",
    "# Print results\n",
    "print('total number evaluated:',i+1)\n",
    "print('support=',support)\n",
    "print('not support=',not_support)\n",
    "try: \n",
    "    percent_support=support/(support+not_support)\n",
    "except: \n",
    "    print('error: support+not_support = 0')\n",
    "\n",
    "print('percent_support=',percent_support)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------- TEMP, test out new Gdisk computF_parm and assignF --------------\n",
    "# Gdisk sandbox code is in dir: C:\\Users\\jaimi\\Documents\\MATLAB\\Research\\topology_work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # # Save RHP variable 'indicMat' to .pkl file, so dont have to run again in future\n",
    "# import pickle\n",
    "# filename='indicMat_8.25.pkl'\n",
    "# with open(filename, 'wb') as f:  # Python 3: open(..., 'wb')\n",
    "#     pickle.dump(indicMat, f)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load indicMat\n",
    "import os\n",
    "import pickle\n",
    "scores = {} # scores is an empty dict already\n",
    "filename='indicMat_8.25.pkl'\n",
    "print(os.path.getsize(filename))\n",
    "if os.path.getsize(filename) > 0:      \n",
    "    with open(filename, \"rb\") as f:\n",
    "        unpickler = pickle.Unpickler(f)\n",
    "        # if file is not empty scores will be equal\n",
    "        # to the value unpickled\n",
    "        indicMat = unpickler.load()\n",
    "        \n",
    "print(len(indicMat)) # just to see if it loaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%notify\n",
    "# ------ run detControlMatExistence ---------\n",
    "t = time.time()\n",
    "\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_ctrlTypes(['PBC'])\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "n=int(len(indicMat)/6)\n",
    "A,B=hm.setupStateSpace(parmObj,fin_feeder, n,node_index_map, depths)\n",
    "t = time.time()\n",
    "#f_ub_list=ctrl.computeFParamSpace_v3(B,indicMat)\n",
    "feas,feasFs,numfeas,numTried,num_act,bestF,indicMat=ctrl.detControlMatExistence(parmObj, A, B, indicMat)\n",
    "elapsed = time.time() - t\n",
    "print('Time elapsed=',elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%notify\n",
    "# ------ run eval_config ---------\n",
    "parmObj=hm.configParms()\n",
    "# parmObj.set_ctrlTypes(['VVC', 'VVC', 'VWC', 'VWC', 'VWC', 'VVC', 'VVC', 'VVC'])\n",
    "# parmObj.set_version(2) # 1 for PBC\n",
    "\n",
    "parmObj.set_ctrlTypes(['PBC', 'PBC', 'PBC'])\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "\n",
    "#all_act_locs = ['bus_152', 'bus_21', 'bus_58', 'bus_70', 'bus_4', 'bus_82', 'bus_109', 'bus_47']\n",
    "all_act_locs = ['bus_47','bus_58','bus_4']\n",
    "perf_nodes=['bus_47','bus_58','bus_4']\n",
    "\n",
    "# print(len(all_act_locs))\n",
    "# print(len(foo))\n",
    "feas, maxError,numfeas,bestF,indicMat=prc.eval_config(parmObj,fin_feeder, all_act_locs,perf_nodes, node_index_map,substation_name,depths,file_name,Vbase_ll, Sbase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next steps:\n",
    "# use N2 to try different perturbations in same percent level\n",
    "# decide what it means for config to be 'feasible' now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply a terrible config and see how many of the 15 work\n",
    "parmObj=hm.configParms()\n",
    "parmObj.set_ctrlTypes(['PBC', 'PBC', 'PBC'])\n",
    "parmObj.set_version(1) # 1 for PBC\n",
    "all_act_locs = ['bus_59','bus_13','bus_53']\n",
    "perf_nodes=['bus_52','bus_60','bus_40']\n",
    "feas, maxError,numfeas,bestF,indicMat=prc.eval_config(parmObj,fin_feeder, all_act_locs,perf_nodes, node_index_map,substation_name,depths,file_name,Vbase_ll, Sbase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}