From 31ddfd79924ffde9575ce1063da002991a93147a Mon Sep 17 00:00:00 2001 From: Priti Paudyal Date: Tue, 27 Dec 2022 15:56:56 -0700 Subject: [PATCH 1/5] Copy EV hosting capacity code into disco --- disco/ev/README.md | 1 + disco/ev/__init__.py | 0 disco/ev/feeder_EV_HC.py | 383 ++++++++++++++++++++++++++++++ disco/ev/load_distance_from_SS.py | 37 +++ disco/ev/number_of_ev_chargers.py | 92 +++++++ disco/ev/plot_hosting_capacity.py | 109 +++++++++ 6 files changed, 622 insertions(+) create mode 100644 disco/ev/README.md create mode 100644 disco/ev/__init__.py create mode 100644 disco/ev/feeder_EV_HC.py create mode 100644 disco/ev/load_distance_from_SS.py create mode 100644 disco/ev/number_of_ev_chargers.py create mode 100644 disco/ev/plot_hosting_capacity.py diff --git a/disco/ev/README.md b/disco/ev/README.md new file mode 100644 index 00000000..26674e2e --- /dev/null +++ b/disco/ev/README.md @@ -0,0 +1 @@ +# EV_hosting_capacity \ No newline at end of file diff --git a/disco/ev/__init__.py b/disco/ev/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/disco/ev/feeder_EV_HC.py b/disco/ev/feeder_EV_HC.py new file mode 100644 index 00000000..2ac65608 --- /dev/null +++ b/disco/ev/feeder_EV_HC.py @@ -0,0 +1,383 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Sep 25 07:53:52 2019 + +@author: ppaudyal +""" + +# -*- coding: utf-8 -*- +""" +Created on Wed Sep 18 14:32:29 2019 + +@author: ppaudyal +""" + + +#this code calculates the hosting capcity (checks voltage and thermal violation limits) for a given feeder +# every feeder model should have a 'Master.dss' as main dss file +# + +import os +import pandas as pd +import numpy as np +import math +import opendssdirect as dss + +import load_distance_from_SS +import plot_hosting_capacity +import number_of_ev_chargers + + +# variable quantities +feeder_folder = '477NHS' #'BaseCaseWRR074' #'BaseaseALD095' +low_lmt_0 = 0.950 # lower voltage limit +high_lmt_0 = 1.050 # upper voltage limit +low_lmt_1 = 0.975 # lower voltage limit +high_lmt_1 = 1.05 # upper voltage limit +low_lmt_2 = 0.985 # lower voltage limit +high_lmt_2 = 1.05 # upper voltage limit +v_limit = [[low_lmt_0, high_lmt_0], [low_lmt_1, high_lmt_1], [low_lmt_2, high_lmt_2]] +step_v = 10 +step_thermal = 5 + +extra_percent = [2]#, 25] # in case, this considers extra % for already overloaded elements +extra_limit = [100]#, 120] #120, if we want to consider overload condition when %Normal>= 120 or any other % + +master_dssfile = 'Master.dss' + +# **************** System Information ***************** +MainDir = os.getcwd() +FeederDir = os.path.join(MainDir, feeder_folder) +MasterFile = os.path.join(FeederDir,master_dssfile) +dss.run_command('clear') +dss.run_command('Compile '+ MasterFile) +dss.run_command('solve') +circuit = dss.Circuit +ckt_name = circuit.Name() +print('\n*** Circuit name:',ckt_name,'***\n from ', feeder_folder, ' folder\n') + +AllNodeNames = circuit.YNodeOrder() +pd.DataFrame(AllNodeNames).to_csv('Allnodenames.csv') +#print(type(AllNodeNames)) + +# --------- Voltage Base kV Information ----- +node_number = len(AllNodeNames) +Vbase_allnode = [0] * node_number +ii = 0 +for node in AllNodeNames: + circuit.SetActiveBus(node) + Vbase_allnode[ii] = dss.Bus.kVBase() * 1000 + ii = ii + 1 + +# --------- Capacitor Information ---------- +capNames = dss.Capacitors.AllNames() +hCapNames = [None] * len(capNames) +for i, n in enumerate(capNames): + hCapNames[i] = str(n) +hCapNames = str(hCapNames)[1:-1] + +# --------- Regulator Information ---------- +regNames = dss.RegControls.AllNames() +hRegNames = [None] * len(regNames) +for i, n in enumerate(regNames): + hRegNames[i] = str(n) +hRegNames = str(hRegNames)[1:-1] + +dss.run_command('solve mode=snap') +#dss.run_command('export voltages') +v = np.array(dss.Circuit.AllBusMagPu()) + + +dss.utils.lines_to_dataframe() + +volt_violation = np.any(v>high_lmt_0) or np.any(vhigh_lmt) or np.any(v=new_threshold[case][indx_]: + thermal_violation = True #just exit here (get out of for loop) + #print(i) + break + else: + thermal_violation = False + else: + #check the percentage normal if greater than specified % then only violation + if report[' %Normal'][i]>= th_limit: + thermal_violation = True + break + else: + thermal_violation = False + + else: + #check the percentage normal if greater than specified % then only violation + for i in range(len(report)): + if report[' %Normal'][i]>= th_limit: + thermal_violation = True + break + else: + thermal_violation = False + + + #print(thermal_violation) # for initial testing + if thermal_violation==True: + #print(thermal_violation) # for initial testing + dss.run_command("export Capacity") + dss.run_command("export Currents") + cap_limit = new_kW + dss.run_command('edit Load.' + str(which_load["name"]) + ' kW=' + str(initial_kW)) + dss.run_command('Compile '+ MasterFile) + dss.run_command('solve mode = snap') + else: + tmp_kW = new_kW + + return[cap_limit] + +######################################################################################################################### +# get the load data +[Load,totalkW] = get_loads(dss,circuit,0,'') + +# calculate the hosting capacity based on voltage constraints +v_output_df = [] +for j in range(len(v_limit)): + lmt1 = v_limit[j][0] + lmt2 = v_limit[j][1] + v_lst = [] + v_output_list = [] + v_threshold = [] + v_allow_limit = [] + v_names = [] + v_bus_name = [] + v_default_load = [] + v_maxv = [] + v_minv = [] + for i in range(len(Load)): + v_node_cap = node_V_capacity_check(Load[i], lmt1, lmt2) # v_node_cap is a list + v_allowable_load = v_node_cap[0] - step_v #10 # v_node_cap[0] is a float + v_threshold.append(v_node_cap[0]) + v_allow_limit.append(v_allowable_load) + v_names.append(Load[i]["name"]) + v_bus_name.append(Load[i]["bus1"]) + v_default_load.append(Load[i]["kW"]) + v_maxv.append(v_node_cap[1]) + v_minv.append(v_node_cap[2]) + + v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] + v_output_list = list(map(list, zip(*v_lst))) + #v_output_df[0] = pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage']) + v_output_df.append(pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage'])) + + v_output_df[j].to_csv("Hosting_capacity_test_" + str(feeder_folder) + "_" + str(lmt1) + ".csv") + +# calculate the hosting capacity based on thermal ratings constraint +th_threshold = [[], []] +th_allow_limit = [[], []] +th_names = [[], []] +th_bus_name = [[], []] +th_default_load = [[], []] + +th_output_df = [] +for i in range(len(extra_limit)): + + th_lst = [] + for j in range(len(Load)): + th_node_cap = thermal_overload_check(Load[j], extra_limit[i], i) # th_node_cap is a list + th_allowable_load = th_node_cap[0] - step_thermal #5 # th_node_cap[0] is a float # reduce the value of add_kW + th_threshold[i].append(th_node_cap[0]) + th_allow_limit[i].append(th_allowable_load) + th_names[i].append(Load[j]["name"]) + th_bus_name[i].append(Load[j]["bus1"]) + th_default_load[i].append(Load[j]["kW"]) + + th_lst = [th_names[i], th_bus_name[i], th_default_load[i], th_threshold[i], th_allow_limit[i]] + th_output_list = list(map(list, zip(*th_lst))) + th_output_df.append(pd.DataFrame(th_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Thermal_Violation', 'Maximum_kW'])) + + th_output_df[i].to_csv("Thermal_capacity_test_" + str(feeder_folder) + "_" + str(extra_limit[i]) + ".csv") + + +############################################## for plotting results #################################################### + +load_bus = pd.DataFrame() +load_bus["Load"] = th_output_df[0]["Load"] # +load_bus["Bus"] = th_output_df[0]["Bus"] +node_distance = pd.DataFrame() +node_distance["Node"] = AllNodeNames +node_distance["Distance"] = Bus_Distance + +dist_file = load_distance_from_SS.calc_dist(load_bus, node_distance) + + + +dist_file["Initial_MW"] = th_output_df[0]["Initial_kW"]/1000 +dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"]/1000 +dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"]/1000 +dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"]/1000 + +dist_file["Thermal_Violation_100"] = th_output_df[0]["Thermal_Violation"]/1000 +#dist_file["Thermal_Violation_120"] = th_output_df[1]["Thermal_Violation"]/1000 + +plot_df = dist_file.sort_values(by=['Distance']) + +#plot voltage violation scenarios +plot_hosting_capacity.plot_capacity_V(plot_df, 'Initial_MW', 'Volt_Violation_0.95', 'Volt_Violation_0.975', 'Volt_Violation_0.985', feeder_folder) + +#plot thermal violation +plot_hosting_capacity.plot_capacity_thermal_1(plot_df, 'Initial_MW', 'Thermal_Violation_100', feeder_folder, 100) +#plot_hosting_capacity.plot_capacity_thermal_2(plot_df, 'Initial_MW', 'Thermal_Violation_100', 'Thermal_Violation_120', feeder_folder) + + + +####################### Assuming the hosting capacity is limited by thermal loading ##################################### + +## Difference of initial load and maximum hosting capacity (assuming always thermal limit occurs first) ## + +for i in range(len(extra_limit)): + diff = th_output_df[i]['Thermal_Violation'] - th_output_df[i]['Initial_kW'] + new_df = pd.DataFrame() + + new_df['Load'] = th_output_df[i]['Load'] + new_df['Bus'] = th_output_df[i]['Bus'] + new_df['Initial_kW'] = th_output_df[i]['Initial_kW'] + new_df['Hosting_capacity(kW)'] = diff # additional load it can support + + new_df.to_csv(str(feeder_folder) + "_Additional_HostingCapacity_" + str(extra_limit[i]) + ".csv", index=False) + + + ## find number of ev chargers for each node ## + + chargers_3_2_1 = number_of_ev_chargers.levels_of_charger(th_output_df[i]) + chargers_3_2_1.to_csv(str(feeder_folder) + "_Loadwithlevel3_2_1_" + str(extra_limit[i]) + ".csv") + diff --git a/disco/ev/load_distance_from_SS.py b/disco/ev/load_distance_from_SS.py new file mode 100644 index 00000000..e413714f --- /dev/null +++ b/disco/ev/load_distance_from_SS.py @@ -0,0 +1,37 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Sep 19 09:03:51 2019 + +@author: ppaudyal +""" + +#import pandas as pd + +def calc_dist(loadandbus_, nodeanddistance): + + nodeanddistance["Node_only"] = range(len(nodeanddistance)) + for i in range(len(nodeanddistance)): + tmp = nodeanddistance["Node"][i] + nodeanddistance["Node_only"][i] = tmp[0:-2] + + nodeanddistance_ = nodeanddistance.drop_duplicates(subset = "Node_only", keep = "first") # if subset="Distance" then any two nodes could have same distance + #print(nodeanddistance_.head(), "\n\n") + #print(len(nodeanddistance_)) + + nodeanddistance_.index = range(len(nodeanddistance_)) + # print(nodeanddistance_.head()) + + + loadandbus_["Distance"] = " " + tmp_lst = (nodeanddistance_["Node_only"].tolist()) + # print(tmp_lst) # testing + for i in range(len(loadandbus_)): + if str(loadandbus_["Bus"][i]) in tmp_lst: + tmp_indx = tmp_lst.index(str(loadandbus_["Bus"][i])) #the index of that node in the list + dst = nodeanddistance_["Distance"][tmp_indx] + loadandbus_["Distance"][i] = dst + else: + print("No ", loadandbus_["Bus"][i], "\n") + + loadandbus_.to_csv("Load_distance_from_SS.csv") # save as csv if required + return loadandbus_ \ No newline at end of file diff --git a/disco/ev/number_of_ev_chargers.py b/disco/ev/number_of_ev_chargers.py new file mode 100644 index 00000000..d73d2adc --- /dev/null +++ b/disco/ev/number_of_ev_chargers.py @@ -0,0 +1,92 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Sep 18 16:20:21 2019 + +@author: ppaudyal +""" +# this function counts the number of EVs each node can support; +# first number of xFC (350 kW), then number of level-2 charging (7.2 kW), and then number of level-1 charging (3.3 kW) + +import pandas as pd + +def levels_of_charger(output_file): + load_name = [] + bus_name = [] + xfc = [] + lvl2 = [] + lvl1 = [] + for i in range(len(output_file)): + load_ = output_file['Load'][i] + bus_ = output_file['Bus'][i] + tmp_data = output_file['Maximum_kW'][i]-output_file['Initial_kW'][i] + + if tmp_data >= 350: # + num3 = tmp_data/350 + no_of_xfc = int(num3) + + tmp_diff = tmp_data - no_of_xfc*350 + if tmp_diff >= 7.2: + #then count number of lvl-2 + n2 = tmp_diff/7.2 + no_of_lvl2 = int(n2) + + tmp_d = tmp_diff - no_of_lvl2*7.2 + if tmp_d >= 3.3: + #then number of lvl-1 + n1 = tmp_d/3.3 + no_of_lvl1 = int(n1) + else: + no_of_lvl1 = 0 + + + elif tmp_diff>= 3.3 and tmp_diff < 7.2: #then number of lvl-1 + no_of_lvl2 = 0 + n1 = tmp_d/3.3 + no_of_lvl1 = int(n1) + else: + no_of_lvl2 = 0 + no_of_lvl1 = 0 + + elif tmp_data >= 7.2 and tmp_data < 350 : # if less than 350 more than 7.2 then how many 7.2 + no_of_xfc = 0 + + n2 = tmp_data/7.2 + no_of_lvl2 = int(n2) + + tmp_d = tmp_data - no_of_lvl2*7.2 + if tmp_d >= 3.3: + #then number of lvl-1 + n1 = tmp_d/3.3 + no_of_lvl1 = int(n1) + else: + no_of_lvl1 = 0 + + elif tmp_data >= 3.3 and tmp_data < 7.2 : ### # if less than 7.2 more than 3.3 + no_of_xfc = 0 + no_of_lvl2 = 0 + print(tmp_data) + + n1 = tmp_data/3.3 + no_of_lvl1 = int(n1) + + else: + no_of_xfc = 0 + no_of_lvl2 = 0 + no_of_lvl1 = 0 + + load_name.append(load_) + bus_name.append(bus_) + xfc.append(no_of_xfc) + lvl2.append(no_of_lvl2) + lvl1.append(no_of_lvl1) + + df = pd.DataFrame() + df["load"] = load_name + df["bus"] = bus_name + df["no. of level3"] = xfc + df["no. of level2"] = lvl2 + df["no. of level1"] = lvl1 + + return df + + \ No newline at end of file diff --git a/disco/ev/plot_hosting_capacity.py b/disco/ev/plot_hosting_capacity.py new file mode 100644 index 00000000..010d68f1 --- /dev/null +++ b/disco/ev/plot_hosting_capacity.py @@ -0,0 +1,109 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Sep 18 15:55:53 2019 + +@author: ppaudyal +""" +import matplotlib.pyplot as plt +import os + + +def plot_capacity_V(file, caseA, caseB, caseC, caseD, feeder): + fig, ax1= plt.subplots(nrows=1,figsize=(14, 8)) #figsize=(18,12) + vplot1 = ax1.plot(range(len(file)), file[caseA], c='k', marker = 'o', label = "Initial load") + vplot2 = ax1.plot(range(len(file)), file[caseB], c='r', marker = '+', label = "Volt violation MW [range:(0.95, 1.05)]") + vplot3 = ax1.plot(range(len(file)), file[caseC], c='g', marker = '*', label = "Volt violation MW [range:(0.975, 1.05)]") + vplot4 = ax1.plot(range(len(file)), file[caseD], c='b', marker = '.', label = "Volt violation MW [range:(0.985, 1.05)]") + ax1.margins(x=0) + ax1.tick_params(axis='both', which='major', labelsize=33) + ax1.tick_params(axis='both', which='minor', labelsize=33) + ax1.set_ylabel('Power (MW)', fontsize=33) + ax1.legend(loc='upper right', fontsize=33) + ax1.set_ylim(-2, 65) + ax1.grid() + #ax1.set_xlabel('Load # according to distance from SS', fontsize=40) + ax1.set_xlabel('Feeder Node (where highest is furthest from substation)', fontsize=33) ##'Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)' + # plt.tight_layout() + plt.savefig(str(feeder) + '_Cap_by_V_limit.png', bbox_inches = 'tight', dpi=150) + +def plot_capacity_thermal_2(file, caseA, caseB, caseC, feeder): + fig, ax1= plt.subplots(nrows=1,figsize=(14, 8)) + thplot1 = ax1.plot(range(len(file)), file[caseA], c='k', marker = 'o', label = "Initial load") + thplot2 = ax1.plot(range(len(file)), file[caseB], c='r', marker = '+', label = "Thermal violation MW[100% loading]") + thplot3 = ax1.plot(range(len(file)), file[caseC], c='b', marker = '*', label = "Thermal violation MW[120% loading]") + ax1.margins(x=0) + ax1.tick_params(axis='both', which='major', labelsize=30) + ax1.tick_params(axis='both', which='minor', labelsize=30) + ax1.set_ylabel('Power (MW)', fontsize=30) + ax1.legend(loc='upper right', fontsize=30) + ax1.grid() + #ax1.set_xlabel('Load # according to distance from SS', fontsize=40) + ax1.set_xlabel('Feeder Node (where highest is furthest from substation)', fontsize=30) #ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) + plt.tight_layout() + plt.savefig(str(feeder) + '_Cap_by_thermal_limit.png', dpi=150) + + +def plot_capacity_thermal_1(file, caseA, caseB, feeder, extra): + fig, ax1= plt.subplots(nrows=1,figsize=(14, 8)) + thplot1 = ax1.plot(range(len(file)), file[caseA], c='k', marker = 'o', label = "Initial load") + thplot2 = ax1.plot(range(len(file)), file[caseB], c='r', marker = '+', label = "Thermal violation MW") + ax1.margins(x=0) + ax1.tick_params(axis='both', which='major', labelsize=33) + ax1.tick_params(axis='both', which='minor', labelsize=33) + ax1.set_ylabel('Power (MW)', fontsize=33) + ax1.legend(loc='upper right', fontsize=33) + ax1.grid() + #ax1.set_xlabel('Load # according to distance from SS', fontsize=40) + ax1.set_xlabel('Feeder Node (where highest is furthest from substation)', fontsize=33) #ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) + # fig.tight_layout() + plt.savefig(str(feeder) + '_Cap_by_thermal_limit_' + str(extra) + '.png', bbox_inches = 'tight', dpi=150) + + +if __name__ == "__main__": + + import pandas as pd + import load_distance_from_SS + + limit = 100 + feeder_folder = 'BaseCaseWRR074' + os.chdir(str(os.getcwd()) + '/' + str(feeder_folder)) + #Read Csv + th_output_df=pd.read_csv("Thermal_capacity_test_" + str(feeder_folder) + "_100.csv") + v_output_df = [] + lmt = [0.95, 0.975, 0.985] + for i in range(len(lmt)): + v_output_df.append(pd.read_csv("Hosting_capacity_test_" + str(feeder_folder) + "_" + str(lmt[i]) + ".csv")) + + bus_dist_df = pd.read_csv("node_distance_" + str(feeder_folder) + ".csv", header=None) + bus_dist_df.columns = ['Distance'] + Bus_Distance = bus_dist_df['Distance'].to_list() + node_name_df = pd.read_csv('Allnodenames.csv') + AllNodeNames = node_name_df['Node'].to_list() + + load_bus = pd.DataFrame() + load_bus["Load"] = th_output_df["Load"] + load_bus["Bus"] = th_output_df["Bus"] #load_bus.to_csv(" ") save as csv if required + node_distance = pd.DataFrame() + node_distance["Node"] = AllNodeNames + node_distance["Distance"] = Bus_Distance #node_distance.to_csv(" ") save as csv if required + + dist_file_lst = load_distance_from_SS.calc_dist(load_bus, node_distance) + + + dist_file = dist_file_lst[0] + + dist_file["Initial_MW"] = th_output_df["Initial_kW"]/1000 + dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"]/1000 + dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"]/1000 + dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"]/1000 + + dist_file["Thermal_Violation_100"] = th_output_df["Thermal_Violation"]/1000 + + plot_df = dist_file.sort_values(by=['Distance']) + + os.chdir('../Results_HC/update_for_paper_jan2021') + #plot voltage violation scenarios + plot_capacity_V(plot_df, 'Initial_MW', 'Volt_Violation_0.95', 'Volt_Violation_0.975', 'Volt_Violation_0.985', feeder_folder) + + #plot thermal violation + plot_capacity_thermal_1(plot_df, 'Initial_MW', 'Thermal_Violation_100', feeder_folder, limit) From 18d09f2c235e0a582955d09e83eae3d928e44272 Mon Sep 17 00:00:00 2001 From: Priti Paudyal Date: Tue, 27 Dec 2022 17:07:50 -0700 Subject: [PATCH 2/5] Run black formatter --- disco/ev/feeder_EV_HC.py | 337 +++++++++++++++++------------- disco/ev/load_distance_from_SS.py | 30 +-- disco/ev/number_of_ev_chargers.py | 76 +++---- disco/ev/plot_hosting_capacity.py | 197 ++++++++++------- 4 files changed, 368 insertions(+), 272 deletions(-) diff --git a/disco/ev/feeder_EV_HC.py b/disco/ev/feeder_EV_HC.py index 2ac65608..8b852889 100644 --- a/disco/ev/feeder_EV_HC.py +++ b/disco/ev/feeder_EV_HC.py @@ -13,14 +13,14 @@ """ -#this code calculates the hosting capcity (checks voltage and thermal violation limits) for a given feeder +# this code calculates the hosting capcity (checks voltage and thermal violation limits) for a given feeder # every feeder model should have a 'Master.dss' as main dss file # import os import pandas as pd import numpy as np -import math +import math import opendssdirect as dss import load_distance_from_SS @@ -29,36 +29,38 @@ # variable quantities -feeder_folder = '477NHS' #'BaseCaseWRR074' #'BaseaseALD095' -low_lmt_0 = 0.950 # lower voltage limit -high_lmt_0 = 1.050 # upper voltage limit -low_lmt_1 = 0.975 # lower voltage limit -high_lmt_1 = 1.05 # upper voltage limit -low_lmt_2 = 0.985 # lower voltage limit -high_lmt_2 = 1.05 # upper voltage limit +feeder_folder = "477NHS" #'BaseCaseWRR074' #'BaseaseALD095' +low_lmt_0 = 0.950 # lower voltage limit +high_lmt_0 = 1.050 # upper voltage limit +low_lmt_1 = 0.975 # lower voltage limit +high_lmt_1 = 1.05 # upper voltage limit +low_lmt_2 = 0.985 # lower voltage limit +high_lmt_2 = 1.05 # upper voltage limit v_limit = [[low_lmt_0, high_lmt_0], [low_lmt_1, high_lmt_1], [low_lmt_2, high_lmt_2]] step_v = 10 step_thermal = 5 -extra_percent = [2]#, 25] # in case, this considers extra % for already overloaded elements -extra_limit = [100]#, 120] #120, if we want to consider overload condition when %Normal>= 120 or any other % +extra_percent = [2] # , 25] # in case, this considers extra % for already overloaded elements +extra_limit = [ + 100 +] # , 120] #120, if we want to consider overload condition when %Normal>= 120 or any other % -master_dssfile = 'Master.dss' +master_dssfile = "Master.dss" # **************** System Information ***************** MainDir = os.getcwd() FeederDir = os.path.join(MainDir, feeder_folder) -MasterFile = os.path.join(FeederDir,master_dssfile) -dss.run_command('clear') -dss.run_command('Compile '+ MasterFile) -dss.run_command('solve') +MasterFile = os.path.join(FeederDir, master_dssfile) +dss.run_command("clear") +dss.run_command("Compile " + MasterFile) +dss.run_command("solve") circuit = dss.Circuit ckt_name = circuit.Name() -print('\n*** Circuit name:',ckt_name,'***\n from ', feeder_folder, ' folder\n') +print("\n*** Circuit name:", ckt_name, "***\n from ", feeder_folder, " folder\n") AllNodeNames = circuit.YNodeOrder() -pd.DataFrame(AllNodeNames).to_csv('Allnodenames.csv') -#print(type(AllNodeNames)) +pd.DataFrame(AllNodeNames).to_csv("Allnodenames.csv") +# print(type(AllNodeNames)) # --------- Voltage Base kV Information ----- node_number = len(AllNodeNames) @@ -83,22 +85,22 @@ hRegNames[i] = str(n) hRegNames = str(hRegNames)[1:-1] -dss.run_command('solve mode=snap') -#dss.run_command('export voltages') +dss.run_command("solve mode=snap") +# dss.run_command('export voltages') v = np.array(dss.Circuit.AllBusMagPu()) dss.utils.lines_to_dataframe() -volt_violation = np.any(v>high_lmt_0) or np.any(v high_lmt_0) or np.any(v < low_lmt_0) +print("Initial condition voltage violation: ", volt_violation) Bus_Distance = [] for node in AllNodeNames: circuit.SetActiveBus(node) # print('node', node, 'bus_distance', dss.Bus.Distance()) # testing Bus_Distance.append(dss.Bus.Distance()) -np.savetxt("node_distance_" + str(feeder_folder) + ".csv", Bus_Distance) +np.savetxt("node_distance_" + str(feeder_folder) + ".csv", Bus_Distance) # --------- Load Information ---------- def get_loads(dss, circuit, loadshape_flag, loadshape_folder): @@ -113,18 +115,22 @@ def get_loads(dss, circuit, loadshape_flag, loadshape_folder): "kV": load.kV(), "kW": load.kW(), "PF": load.PF(), - "Delta_conn": load.IsDelta() + "Delta_conn": load.IsDelta(), } cktElement = dss.CktElement bus = cktElement.BusNames()[0].split(".") - datum["kVar"] = float(datum["kW"]) / float(datum["PF"]) * math.sqrt(1 - float(datum["PF"]) * float(datum["PF"])) + datum["kVar"] = ( + float(datum["kW"]) + / float(datum["PF"]) + * math.sqrt(1 - float(datum["PF"]) * float(datum["PF"])) + ) datum["bus1"] = bus[0] datum["numPhases"] = len(bus[1:]) datum["phases"] = bus[1:] if not datum["numPhases"]: datum["numPhases"] = 3 - datum["phases"] = ['1', '2', '3'] + datum["phases"] = ["1", "2", "3"] datum["voltageMag"] = cktElement.VoltagesMagAng()[0] datum["voltageAng"] = cktElement.VoltagesMagAng()[1] datum["power"] = dss.CktElement.Powers()[0:2] @@ -138,136 +144,139 @@ def get_loads(dss, circuit, loadshape_flag, loadshape_folder): ############################################## for voltage violation capacity ###################################### def node_V_capacity_check(which_load, low_lmt, high_lmt): - initial_kW = which_load['kW'] - #print(initial_kW) #just for testing - add_kW = step_v #10 + initial_kW = which_load["kW"] + # print(initial_kW) #just for testing + add_kW = step_v # 10 tmp_kW = initial_kW volt_violation = False - v=None + v = None - while volt_violation==False: - # while the voltages are within limits keep on increasing - new_kW = tmp_kW + add_kW # increase the load by 10 kW each time - dss.run_command('edit Load.' + str(which_load["name"]) + ' kW=' + str(new_kW)) - dss.run_command('solve mode = snap') + while volt_violation == False: + # while the voltages are within limits keep on increasing + new_kW = tmp_kW + add_kW # increase the load by 10 kW each time + dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) + dss.run_command("solve mode = snap") v = np.array(dss.Circuit.AllBusMagPu()) - volt_violation = np.any(v>high_lmt) or np.any(v high_lmt) or np.any(v < low_lmt) + # print(volt_violation) # for initial testing + if volt_violation == True: + # print(volt_violation) # for initial testing vmax = np.max(v) vmin = np.min(v) cap_limit = new_kW - dss.run_command('edit Load.' + str(which_load["name"]) + ' kW=' + str(initial_kW)) - dss.run_command('Compile '+ MasterFile) - dss.run_command('solve mode = snap') + dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) + dss.run_command("Compile " + MasterFile) + dss.run_command("solve mode = snap") else: - tmp_kW = new_kW + tmp_kW = new_kW - return[cap_limit, vmax, vmin] - - -################################# for thermal overload capapcity ####################################################### -dss.run_command('solve mode = snap') + return [cap_limit, vmax, vmin] + + +################################# for thermal overload capapcity ####################################################### +dss.run_command("solve mode = snap") dss.run_command("export Overloads") -#rename this overload csv file +# rename this overload csv file src = circuit.Name() + "_EXP_OVERLOADS.CSV" dst = feeder_folder + "_overloads.csv" if os.path.exists(dst): os.remove(dst) os.rename(src, dst) -#read this overload file and record the ' %Normal' for each line in this file +# read this overload file and record the ' %Normal' for each line in this file overload_df = pd.read_csv(dst) -#len(overload_df) -if(len(overload_df)==0): +# len(overload_df) +if len(overload_df) == 0: print("No thermal violation initially \n") else: print("Thermal violation exists initially \n") print(overload_df.head()) - + elements = [[], []] amps = [[], []] new_threshold = [[], []] for j in range(len(extra_percent)): - + for i in range(len(overload_df)): - overload_df['Element'] = overload_df['Element'].str.strip() # removing the extra spaces at the end (if any) - element = overload_df['Element'][i] - amp = overload_df[' %Normal'][i] + overload_df["Element"] = overload_df[ + "Element" + ].str.strip() # removing the extra spaces at the end (if any) + element = overload_df["Element"][i] + amp = overload_df[" %Normal"][i] element_new_limit = amp + extra_percent[j] elements[j].append(str(element)) amps[j].append(amp) new_threshold[j].append(element_new_limit) - -print('new_threshold', new_threshold) # testing + +print("new_threshold", new_threshold) # testing + def thermal_overload_check(which_load, th_limit, case): - initial_kW = which_load['kW'] - #print(initial_kW) - add_kW = step_thermal #5 + initial_kW = which_load["kW"] + # print(initial_kW) + add_kW = step_thermal # 5 tmp_kW = initial_kW thermal_violation = False - - while thermal_violation==False: - # while the elements loadings are within limits keep on increasing - new_kW = tmp_kW + add_kW # increase the load by 5 kW each time - dss.run_command('edit Load.' + str(which_load["name"]) + ' kW=' + str(new_kW)) - dss.run_command('solve mode = snap') + + while thermal_violation == False: + # while the elements loadings are within limits keep on increasing + new_kW = tmp_kW + add_kW # increase the load by 5 kW each time + dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) + dss.run_command("solve mode = snap") dss.run_command("export Overloads") - report = pd.read_csv(str(ckt_name) + "_EXP_OVERLOADS.CSV") - report['Element'] = report['Element'].str.strip() - - if(len(report)==0): # if no any overload element + report = pd.read_csv(str(ckt_name) + "_EXP_OVERLOADS.CSV") + report["Element"] = report["Element"].str.strip() + + if len(report) == 0: # if no any overload element thermal_violation = False - - elif report['Element'].isin(elements[case]).any(): + + elif report["Element"].isin(elements[case]).any(): for i in range(len(report)): - if report['Element'][i] in elements[case]: - indx_ = elements[case].index(report['Element'][i]) #find the index of - if report[' %Normal'][i]>=new_threshold[case][indx_]: - thermal_violation = True #just exit here (get out of for loop) - #print(i) + if report["Element"][i] in elements[case]: + indx_ = elements[case].index(report["Element"][i]) # find the index of + if report[" %Normal"][i] >= new_threshold[case][indx_]: + thermal_violation = True # just exit here (get out of for loop) + # print(i) break else: thermal_violation = False else: - #check the percentage normal if greater than specified % then only violation - if report[' %Normal'][i]>= th_limit: - thermal_violation = True + # check the percentage normal if greater than specified % then only violation + if report[" %Normal"][i] >= th_limit: + thermal_violation = True break else: - thermal_violation = False - + thermal_violation = False + else: - #check the percentage normal if greater than specified % then only violation + # check the percentage normal if greater than specified % then only violation for i in range(len(report)): - if report[' %Normal'][i]>= th_limit: - thermal_violation = True + if report[" %Normal"][i] >= th_limit: + thermal_violation = True break else: - thermal_violation = False - + thermal_violation = False - #print(thermal_violation) # for initial testing - if thermal_violation==True: - #print(thermal_violation) # for initial testing + # print(thermal_violation) # for initial testing + if thermal_violation == True: + # print(thermal_violation) # for initial testing dss.run_command("export Capacity") dss.run_command("export Currents") cap_limit = new_kW - dss.run_command('edit Load.' + str(which_load["name"]) + ' kW=' + str(initial_kW)) - dss.run_command('Compile '+ MasterFile) - dss.run_command('solve mode = snap') + dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) + dss.run_command("Compile " + MasterFile) + dss.run_command("solve mode = snap") else: - tmp_kW = new_kW + tmp_kW = new_kW + + return [cap_limit] - return[cap_limit] - -######################################################################################################################### -# get the load data -[Load,totalkW] = get_loads(dss,circuit,0,'') + +######################################################################################################################### +# get the load data +[Load, totalkW] = get_loads(dss, circuit, 0, "") # calculate the hosting capacity based on voltage constraints v_output_df = [] @@ -284,21 +293,34 @@ def thermal_overload_check(which_load, th_limit, case): v_maxv = [] v_minv = [] for i in range(len(Load)): - v_node_cap = node_V_capacity_check(Load[i], lmt1, lmt2) # v_node_cap is a list - v_allowable_load = v_node_cap[0] - step_v #10 # v_node_cap[0] is a float + v_node_cap = node_V_capacity_check(Load[i], lmt1, lmt2) # v_node_cap is a list + v_allowable_load = v_node_cap[0] - step_v # 10 # v_node_cap[0] is a float v_threshold.append(v_node_cap[0]) v_allow_limit.append(v_allowable_load) v_names.append(Load[i]["name"]) v_bus_name.append(Load[i]["bus1"]) v_default_load.append(Load[i]["kW"]) - v_maxv.append(v_node_cap[1]) + v_maxv.append(v_node_cap[1]) v_minv.append(v_node_cap[2]) - - v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] - v_output_list = list(map(list, zip(*v_lst))) - #v_output_df[0] = pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage']) - v_output_df.append(pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage'])) - + + v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] + v_output_list = list(map(list, zip(*v_lst))) + # v_output_df[0] = pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage']) + v_output_df.append( + pd.DataFrame( + v_output_list, + columns=[ + "Load", + "Bus", + "Initial_kW", + "Volt_Violation", + "Maximum_kW", + "Max_voltage", + "Min_voltage", + ], + ) + ) + v_output_df[j].to_csv("Hosting_capacity_test_" + str(feeder_folder) + "_" + str(lmt1) + ".csv") # calculate the hosting capacity based on thermal ratings constraint @@ -310,54 +332,70 @@ def thermal_overload_check(which_load, th_limit, case): th_output_df = [] for i in range(len(extra_limit)): - + th_lst = [] for j in range(len(Load)): - th_node_cap = thermal_overload_check(Load[j], extra_limit[i], i) # th_node_cap is a list - th_allowable_load = th_node_cap[0] - step_thermal #5 # th_node_cap[0] is a float # reduce the value of add_kW + th_node_cap = thermal_overload_check(Load[j], extra_limit[i], i) # th_node_cap is a list + th_allowable_load = ( + th_node_cap[0] - step_thermal + ) # 5 # th_node_cap[0] is a float # reduce the value of add_kW th_threshold[i].append(th_node_cap[0]) th_allow_limit[i].append(th_allowable_load) th_names[i].append(Load[j]["name"]) th_bus_name[i].append(Load[j]["bus1"]) th_default_load[i].append(Load[j]["kW"]) - - th_lst = [th_names[i], th_bus_name[i], th_default_load[i], th_threshold[i], th_allow_limit[i]] - th_output_list = list(map(list, zip(*th_lst))) - th_output_df.append(pd.DataFrame(th_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Thermal_Violation', 'Maximum_kW'])) - - th_output_df[i].to_csv("Thermal_capacity_test_" + str(feeder_folder) + "_" + str(extra_limit[i]) + ".csv") - + + th_lst = [th_names[i], th_bus_name[i], th_default_load[i], th_threshold[i], th_allow_limit[i]] + th_output_list = list(map(list, zip(*th_lst))) + th_output_df.append( + pd.DataFrame( + th_output_list, + columns=["Load", "Bus", "Initial_kW", "Thermal_Violation", "Maximum_kW"], + ) + ) + + th_output_df[i].to_csv( + "Thermal_capacity_test_" + str(feeder_folder) + "_" + str(extra_limit[i]) + ".csv" + ) + ############################################## for plotting results #################################################### load_bus = pd.DataFrame() -load_bus["Load"] = th_output_df[0]["Load"] # -load_bus["Bus"] = th_output_df[0]["Bus"] +load_bus["Load"] = th_output_df[0]["Load"] # +load_bus["Bus"] = th_output_df[0]["Bus"] node_distance = pd.DataFrame() node_distance["Node"] = AllNodeNames -node_distance["Distance"] = Bus_Distance +node_distance["Distance"] = Bus_Distance dist_file = load_distance_from_SS.calc_dist(load_bus, node_distance) +dist_file["Initial_MW"] = th_output_df[0]["Initial_kW"] / 1000 +dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"] / 1000 +dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"] / 1000 +dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"] / 1000 -dist_file["Initial_MW"] = th_output_df[0]["Initial_kW"]/1000 -dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"]/1000 -dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"]/1000 -dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"]/1000 - -dist_file["Thermal_Violation_100"] = th_output_df[0]["Thermal_Violation"]/1000 -#dist_file["Thermal_Violation_120"] = th_output_df[1]["Thermal_Violation"]/1000 - -plot_df = dist_file.sort_values(by=['Distance']) +dist_file["Thermal_Violation_100"] = th_output_df[0]["Thermal_Violation"] / 1000 +# dist_file["Thermal_Violation_120"] = th_output_df[1]["Thermal_Violation"]/1000 -#plot voltage violation scenarios -plot_hosting_capacity.plot_capacity_V(plot_df, 'Initial_MW', 'Volt_Violation_0.95', 'Volt_Violation_0.975', 'Volt_Violation_0.985', feeder_folder) +plot_df = dist_file.sort_values(by=["Distance"]) -#plot thermal violation -plot_hosting_capacity.plot_capacity_thermal_1(plot_df, 'Initial_MW', 'Thermal_Violation_100', feeder_folder, 100) -#plot_hosting_capacity.plot_capacity_thermal_2(plot_df, 'Initial_MW', 'Thermal_Violation_100', 'Thermal_Violation_120', feeder_folder) +# plot voltage violation scenarios +plot_hosting_capacity.plot_capacity_V( + plot_df, + "Initial_MW", + "Volt_Violation_0.95", + "Volt_Violation_0.975", + "Volt_Violation_0.985", + feeder_folder, +) +# plot thermal violation +plot_hosting_capacity.plot_capacity_thermal_1( + plot_df, "Initial_MW", "Thermal_Violation_100", feeder_folder, 100 +) +# plot_hosting_capacity.plot_capacity_thermal_2(plot_df, 'Initial_MW', 'Thermal_Violation_100', 'Thermal_Violation_120', feeder_folder) ####################### Assuming the hosting capacity is limited by thermal loading ##################################### @@ -365,19 +403,22 @@ def thermal_overload_check(which_load, th_limit, case): ## Difference of initial load and maximum hosting capacity (assuming always thermal limit occurs first) ## for i in range(len(extra_limit)): - diff = th_output_df[i]['Thermal_Violation'] - th_output_df[i]['Initial_kW'] + diff = th_output_df[i]["Thermal_Violation"] - th_output_df[i]["Initial_kW"] new_df = pd.DataFrame() - - new_df['Load'] = th_output_df[i]['Load'] - new_df['Bus'] = th_output_df[i]['Bus'] - new_df['Initial_kW'] = th_output_df[i]['Initial_kW'] - new_df['Hosting_capacity(kW)'] = diff # additional load it can support - new_df.to_csv(str(feeder_folder) + "_Additional_HostingCapacity_" + str(extra_limit[i]) + ".csv", index=False) + new_df["Load"] = th_output_df[i]["Load"] + new_df["Bus"] = th_output_df[i]["Bus"] + new_df["Initial_kW"] = th_output_df[i]["Initial_kW"] + new_df["Hosting_capacity(kW)"] = diff # additional load it can support + new_df.to_csv( + str(feeder_folder) + "_Additional_HostingCapacity_" + str(extra_limit[i]) + ".csv", + index=False, + ) ## find number of ev chargers for each node ## - - chargers_3_2_1 = number_of_ev_chargers.levels_of_charger(th_output_df[i]) - chargers_3_2_1.to_csv(str(feeder_folder) + "_Loadwithlevel3_2_1_" + str(extra_limit[i]) + ".csv") + chargers_3_2_1 = number_of_ev_chargers.levels_of_charger(th_output_df[i]) + chargers_3_2_1.to_csv( + str(feeder_folder) + "_Loadwithlevel3_2_1_" + str(extra_limit[i]) + ".csv" + ) diff --git a/disco/ev/load_distance_from_SS.py b/disco/ev/load_distance_from_SS.py index e413714f..57e60706 100644 --- a/disco/ev/load_distance_from_SS.py +++ b/disco/ev/load_distance_from_SS.py @@ -5,33 +5,37 @@ @author: ppaudyal """ -#import pandas as pd +# import pandas as pd + def calc_dist(loadandbus_, nodeanddistance): - + nodeanddistance["Node_only"] = range(len(nodeanddistance)) for i in range(len(nodeanddistance)): - tmp = nodeanddistance["Node"][i] + tmp = nodeanddistance["Node"][i] nodeanddistance["Node_only"][i] = tmp[0:-2] - - nodeanddistance_ = nodeanddistance.drop_duplicates(subset = "Node_only", keep = "first") # if subset="Distance" then any two nodes could have same distance - #print(nodeanddistance_.head(), "\n\n") - #print(len(nodeanddistance_)) - + + nodeanddistance_ = nodeanddistance.drop_duplicates( + subset="Node_only", keep="first" + ) # if subset="Distance" then any two nodes could have same distance + # print(nodeanddistance_.head(), "\n\n") + # print(len(nodeanddistance_)) + nodeanddistance_.index = range(len(nodeanddistance_)) # print(nodeanddistance_.head()) - loadandbus_["Distance"] = " " - tmp_lst = (nodeanddistance_["Node_only"].tolist()) + tmp_lst = nodeanddistance_["Node_only"].tolist() # print(tmp_lst) # testing for i in range(len(loadandbus_)): if str(loadandbus_["Bus"][i]) in tmp_lst: - tmp_indx = tmp_lst.index(str(loadandbus_["Bus"][i])) #the index of that node in the list + tmp_indx = tmp_lst.index( + str(loadandbus_["Bus"][i]) + ) # the index of that node in the list dst = nodeanddistance_["Distance"][tmp_indx] loadandbus_["Distance"][i] = dst else: print("No ", loadandbus_["Bus"][i], "\n") - + loadandbus_.to_csv("Load_distance_from_SS.csv") # save as csv if required - return loadandbus_ \ No newline at end of file + return loadandbus_ diff --git a/disco/ev/number_of_ev_chargers.py b/disco/ev/number_of_ev_chargers.py index d73d2adc..6b1a28e4 100644 --- a/disco/ev/number_of_ev_chargers.py +++ b/disco/ev/number_of_ev_chargers.py @@ -4,11 +4,12 @@ @author: ppaudyal """ -# this function counts the number of EVs each node can support; +# this function counts the number of EVs each node can support; # first number of xFC (350 kW), then number of level-2 charging (7.2 kW), and then number of level-1 charging (3.3 kW) import pandas as pd + def levels_of_charger(output_file): load_name = [] bus_name = [] @@ -16,70 +17,71 @@ def levels_of_charger(output_file): lvl2 = [] lvl1 = [] for i in range(len(output_file)): - load_ = output_file['Load'][i] - bus_ = output_file['Bus'][i] - tmp_data = output_file['Maximum_kW'][i]-output_file['Initial_kW'][i] - - if tmp_data >= 350: # - num3 = tmp_data/350 + load_ = output_file["Load"][i] + bus_ = output_file["Bus"][i] + tmp_data = output_file["Maximum_kW"][i] - output_file["Initial_kW"][i] + + if tmp_data >= 350: # + num3 = tmp_data / 350 no_of_xfc = int(num3) - - tmp_diff = tmp_data - no_of_xfc*350 + + tmp_diff = tmp_data - no_of_xfc * 350 if tmp_diff >= 7.2: - #then count number of lvl-2 - n2 = tmp_diff/7.2 + # then count number of lvl-2 + n2 = tmp_diff / 7.2 no_of_lvl2 = int(n2) - - tmp_d = tmp_diff - no_of_lvl2*7.2 - if tmp_d >= 3.3: - #then number of lvl-1 - n1 = tmp_d/3.3 + + tmp_d = tmp_diff - no_of_lvl2 * 7.2 + if tmp_d >= 3.3: + # then number of lvl-1 + n1 = tmp_d / 3.3 no_of_lvl1 = int(n1) else: no_of_lvl1 = 0 - - - elif tmp_diff>= 3.3 and tmp_diff < 7.2: #then number of lvl-1 + + elif tmp_diff >= 3.3 and tmp_diff < 7.2: # then number of lvl-1 no_of_lvl2 = 0 - n1 = tmp_d/3.3 + n1 = tmp_d / 3.3 no_of_lvl1 = int(n1) else: no_of_lvl2 = 0 no_of_lvl1 = 0 - - elif tmp_data >= 7.2 and tmp_data < 350 : # if less than 350 more than 7.2 then how many 7.2 + + elif ( + tmp_data >= 7.2 and tmp_data < 350 + ): # if less than 350 more than 7.2 then how many 7.2 no_of_xfc = 0 - - n2 = tmp_data/7.2 + + n2 = tmp_data / 7.2 no_of_lvl2 = int(n2) - - tmp_d = tmp_data - no_of_lvl2*7.2 + + tmp_d = tmp_data - no_of_lvl2 * 7.2 if tmp_d >= 3.3: - #then number of lvl-1 - n1 = tmp_d/3.3 - no_of_lvl1 = int(n1) + # then number of lvl-1 + n1 = tmp_d / 3.3 + no_of_lvl1 = int(n1) else: no_of_lvl1 = 0 - - elif tmp_data >= 3.3 and tmp_data < 7.2 : ### # if less than 7.2 more than 3.3 + + elif tmp_data >= 3.3 and tmp_data < 7.2: ### # if less than 7.2 more than 3.3 no_of_xfc = 0 no_of_lvl2 = 0 print(tmp_data) - - n1 = tmp_data/3.3 + + n1 = tmp_data / 3.3 no_of_lvl1 = int(n1) - + else: no_of_xfc = 0 no_of_lvl2 = 0 no_of_lvl1 = 0 - + load_name.append(load_) bus_name.append(bus_) xfc.append(no_of_xfc) lvl2.append(no_of_lvl2) lvl1.append(no_of_lvl1) - + df = pd.DataFrame() df["load"] = load_name df["bus"] = bus_name @@ -88,5 +90,3 @@ def levels_of_charger(output_file): df["no. of level1"] = lvl1 return df - - \ No newline at end of file diff --git a/disco/ev/plot_hosting_capacity.py b/disco/ev/plot_hosting_capacity.py index 010d68f1..28f7e1e6 100644 --- a/disco/ev/plot_hosting_capacity.py +++ b/disco/ev/plot_hosting_capacity.py @@ -8,102 +8,153 @@ import os -def plot_capacity_V(file, caseA, caseB, caseC, caseD, feeder): - fig, ax1= plt.subplots(nrows=1,figsize=(14, 8)) #figsize=(18,12) - vplot1 = ax1.plot(range(len(file)), file[caseA], c='k', marker = 'o', label = "Initial load") - vplot2 = ax1.plot(range(len(file)), file[caseB], c='r', marker = '+', label = "Volt violation MW [range:(0.95, 1.05)]") - vplot3 = ax1.plot(range(len(file)), file[caseC], c='g', marker = '*', label = "Volt violation MW [range:(0.975, 1.05)]") - vplot4 = ax1.plot(range(len(file)), file[caseD], c='b', marker = '.', label = "Volt violation MW [range:(0.985, 1.05)]") +def plot_capacity_V(file, caseA, caseB, caseC, caseD, feeder): + fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) # figsize=(18,12) + vplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") + vplot2 = ax1.plot( + range(len(file)), + file[caseB], + c="r", + marker="+", + label="Volt violation MW [range:(0.95, 1.05)]", + ) + vplot3 = ax1.plot( + range(len(file)), + file[caseC], + c="g", + marker="*", + label="Volt violation MW [range:(0.975, 1.05)]", + ) + vplot4 = ax1.plot( + range(len(file)), + file[caseD], + c="b", + marker=".", + label="Volt violation MW [range:(0.985, 1.05)]", + ) ax1.margins(x=0) - ax1.tick_params(axis='both', which='major', labelsize=33) - ax1.tick_params(axis='both', which='minor', labelsize=33) - ax1.set_ylabel('Power (MW)', fontsize=33) - ax1.legend(loc='upper right', fontsize=33) + ax1.tick_params(axis="both", which="major", labelsize=33) + ax1.tick_params(axis="both", which="minor", labelsize=33) + ax1.set_ylabel("Power (MW)", fontsize=33) + ax1.legend(loc="upper right", fontsize=33) ax1.set_ylim(-2, 65) ax1.grid() - #ax1.set_xlabel('Load # according to distance from SS', fontsize=40) - ax1.set_xlabel('Feeder Node (where highest is furthest from substation)', fontsize=33) ##'Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)' + # ax1.set_xlabel('Load # according to distance from SS', fontsize=40) + ax1.set_xlabel( + "Feeder Node (where highest is furthest from substation)", fontsize=33 + ) ##'Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)' # plt.tight_layout() - plt.savefig(str(feeder) + '_Cap_by_V_limit.png', bbox_inches = 'tight', dpi=150) + plt.savefig(str(feeder) + "_Cap_by_V_limit.png", bbox_inches="tight", dpi=150) -def plot_capacity_thermal_2(file, caseA, caseB, caseC, feeder): - fig, ax1= plt.subplots(nrows=1,figsize=(14, 8)) - thplot1 = ax1.plot(range(len(file)), file[caseA], c='k', marker = 'o', label = "Initial load") - thplot2 = ax1.plot(range(len(file)), file[caseB], c='r', marker = '+', label = "Thermal violation MW[100% loading]") - thplot3 = ax1.plot(range(len(file)), file[caseC], c='b', marker = '*', label = "Thermal violation MW[120% loading]") + +def plot_capacity_thermal_2(file, caseA, caseB, caseC, feeder): + fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) + thplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") + thplot2 = ax1.plot( + range(len(file)), + file[caseB], + c="r", + marker="+", + label="Thermal violation MW[100% loading]", + ) + thplot3 = ax1.plot( + range(len(file)), + file[caseC], + c="b", + marker="*", + label="Thermal violation MW[120% loading]", + ) ax1.margins(x=0) - ax1.tick_params(axis='both', which='major', labelsize=30) - ax1.tick_params(axis='both', which='minor', labelsize=30) - ax1.set_ylabel('Power (MW)', fontsize=30) - ax1.legend(loc='upper right', fontsize=30) + ax1.tick_params(axis="both", which="major", labelsize=30) + ax1.tick_params(axis="both", which="minor", labelsize=30) + ax1.set_ylabel("Power (MW)", fontsize=30) + ax1.legend(loc="upper right", fontsize=30) ax1.grid() - #ax1.set_xlabel('Load # according to distance from SS', fontsize=40) - ax1.set_xlabel('Feeder Node (where highest is furthest from substation)', fontsize=30) #ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) + # ax1.set_xlabel('Load # according to distance from SS', fontsize=40) + ax1.set_xlabel( + "Feeder Node (where highest is furthest from substation)", fontsize=30 + ) # ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) plt.tight_layout() - plt.savefig(str(feeder) + '_Cap_by_thermal_limit.png', dpi=150) + plt.savefig(str(feeder) + "_Cap_by_thermal_limit.png", dpi=150) -def plot_capacity_thermal_1(file, caseA, caseB, feeder, extra): - fig, ax1= plt.subplots(nrows=1,figsize=(14, 8)) - thplot1 = ax1.plot(range(len(file)), file[caseA], c='k', marker = 'o', label = "Initial load") - thplot2 = ax1.plot(range(len(file)), file[caseB], c='r', marker = '+', label = "Thermal violation MW") +def plot_capacity_thermal_1(file, caseA, caseB, feeder, extra): + fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) + thplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") + thplot2 = ax1.plot( + range(len(file)), file[caseB], c="r", marker="+", label="Thermal violation MW" + ) ax1.margins(x=0) - ax1.tick_params(axis='both', which='major', labelsize=33) - ax1.tick_params(axis='both', which='minor', labelsize=33) - ax1.set_ylabel('Power (MW)', fontsize=33) - ax1.legend(loc='upper right', fontsize=33) + ax1.tick_params(axis="both", which="major", labelsize=33) + ax1.tick_params(axis="both", which="minor", labelsize=33) + ax1.set_ylabel("Power (MW)", fontsize=33) + ax1.legend(loc="upper right", fontsize=33) ax1.grid() - #ax1.set_xlabel('Load # according to distance from SS', fontsize=40) - ax1.set_xlabel('Feeder Node (where highest is furthest from substation)', fontsize=33) #ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) + # ax1.set_xlabel('Load # according to distance from SS', fontsize=40) + ax1.set_xlabel( + "Feeder Node (where highest is furthest from substation)", fontsize=33 + ) # ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) # fig.tight_layout() - plt.savefig(str(feeder) + '_Cap_by_thermal_limit_' + str(extra) + '.png', bbox_inches = 'tight', dpi=150) + plt.savefig( + str(feeder) + "_Cap_by_thermal_limit_" + str(extra) + ".png", bbox_inches="tight", dpi=150 + ) -if __name__ == "__main__": +if __name__ == "__main__": - import pandas as pd + import pandas as pd import load_distance_from_SS - + limit = 100 - feeder_folder = 'BaseCaseWRR074' - os.chdir(str(os.getcwd()) + '/' + str(feeder_folder)) - #Read Csv - th_output_df=pd.read_csv("Thermal_capacity_test_" + str(feeder_folder) + "_100.csv") + feeder_folder = "BaseCaseWRR074" + os.chdir(str(os.getcwd()) + "/" + str(feeder_folder)) + # Read Csv + th_output_df = pd.read_csv("Thermal_capacity_test_" + str(feeder_folder) + "_100.csv") v_output_df = [] lmt = [0.95, 0.975, 0.985] for i in range(len(lmt)): - v_output_df.append(pd.read_csv("Hosting_capacity_test_" + str(feeder_folder) + "_" + str(lmt[i]) + ".csv")) - - bus_dist_df = pd.read_csv("node_distance_" + str(feeder_folder) + ".csv", header=None) - bus_dist_df.columns = ['Distance'] - Bus_Distance = bus_dist_df['Distance'].to_list() - node_name_df = pd.read_csv('Allnodenames.csv') - AllNodeNames = node_name_df['Node'].to_list() - + v_output_df.append( + pd.read_csv("Hosting_capacity_test_" + str(feeder_folder) + "_" + str(lmt[i]) + ".csv") + ) + + bus_dist_df = pd.read_csv("node_distance_" + str(feeder_folder) + ".csv", header=None) + bus_dist_df.columns = ["Distance"] + Bus_Distance = bus_dist_df["Distance"].to_list() + node_name_df = pd.read_csv("Allnodenames.csv") + AllNodeNames = node_name_df["Node"].to_list() + load_bus = pd.DataFrame() - load_bus["Load"] = th_output_df["Load"] - load_bus["Bus"] = th_output_df["Bus"] #load_bus.to_csv(" ") save as csv if required + load_bus["Load"] = th_output_df["Load"] + load_bus["Bus"] = th_output_df["Bus"] # load_bus.to_csv(" ") save as csv if required node_distance = pd.DataFrame() node_distance["Node"] = AllNodeNames - node_distance["Distance"] = Bus_Distance #node_distance.to_csv(" ") save as csv if required - + node_distance[ + "Distance" + ] = Bus_Distance # node_distance.to_csv(" ") save as csv if required + dist_file_lst = load_distance_from_SS.calc_dist(load_bus, node_distance) - - + dist_file = dist_file_lst[0] - - dist_file["Initial_MW"] = th_output_df["Initial_kW"]/1000 - dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"]/1000 - dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"]/1000 - dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"]/1000 - - dist_file["Thermal_Violation_100"] = th_output_df["Thermal_Violation"]/1000 - - plot_df = dist_file.sort_values(by=['Distance']) - - os.chdir('../Results_HC/update_for_paper_jan2021') - #plot voltage violation scenarios - plot_capacity_V(plot_df, 'Initial_MW', 'Volt_Violation_0.95', 'Volt_Violation_0.975', 'Volt_Violation_0.985', feeder_folder) - - #plot thermal violation - plot_capacity_thermal_1(plot_df, 'Initial_MW', 'Thermal_Violation_100', feeder_folder, limit) + + dist_file["Initial_MW"] = th_output_df["Initial_kW"] / 1000 + dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"] / 1000 + dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"] / 1000 + dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"] / 1000 + + dist_file["Thermal_Violation_100"] = th_output_df["Thermal_Violation"] / 1000 + + plot_df = dist_file.sort_values(by=["Distance"]) + + os.chdir("../Results_HC/update_for_paper_jan2021") + # plot voltage violation scenarios + plot_capacity_V( + plot_df, + "Initial_MW", + "Volt_Violation_0.95", + "Volt_Violation_0.975", + "Volt_Violation_0.985", + feeder_folder, + ) + + # plot thermal violation + plot_capacity_thermal_1(plot_df, "Initial_MW", "Thermal_Violation_100", feeder_folder, limit) From c0b04187ded3a7eb4faa355f0239b48191421713 Mon Sep 17 00:00:00 2001 From: Daniel Thom Date: Tue, 27 Dec 2022 17:51:23 -0700 Subject: [PATCH 3/5] Use dss.Text.Command for error checking --- disco/cli/ev.py | 0 disco/ev/feeder_EV_HC.py | 40 ++++++++++++++++++++-------------------- 2 files changed, 20 insertions(+), 20 deletions(-) create mode 100644 disco/cli/ev.py diff --git a/disco/cli/ev.py b/disco/cli/ev.py new file mode 100644 index 00000000..e69de29b diff --git a/disco/ev/feeder_EV_HC.py b/disco/ev/feeder_EV_HC.py index 8b852889..5e268be1 100644 --- a/disco/ev/feeder_EV_HC.py +++ b/disco/ev/feeder_EV_HC.py @@ -51,9 +51,9 @@ MainDir = os.getcwd() FeederDir = os.path.join(MainDir, feeder_folder) MasterFile = os.path.join(FeederDir, master_dssfile) -dss.run_command("clear") -dss.run_command("Compile " + MasterFile) -dss.run_command("solve") +dss.Text.Command("clear") +dss.Text.Command(f"Compile {MasterFile}") +dss.Text.Command("solve") circuit = dss.Circuit ckt_name = circuit.Name() print("\n*** Circuit name:", ckt_name, "***\n from ", feeder_folder, " folder\n") @@ -85,8 +85,8 @@ hRegNames[i] = str(n) hRegNames = str(hRegNames)[1:-1] -dss.run_command("solve mode=snap") -# dss.run_command('export voltages') +dss.Text.Command("solve mode=snap") +# dss.Text.Command('export voltages') v = np.array(dss.Circuit.AllBusMagPu()) @@ -154,8 +154,8 @@ def node_V_capacity_check(which_load, low_lmt, high_lmt): while volt_violation == False: # while the voltages are within limits keep on increasing new_kW = tmp_kW + add_kW # increase the load by 10 kW each time - dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) - dss.run_command("solve mode = snap") + dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) + dss.Text.Command("solve mode = snap") v = np.array(dss.Circuit.AllBusMagPu()) volt_violation = np.any(v > high_lmt) or np.any(v < low_lmt) # print(volt_violation) # for initial testing @@ -164,9 +164,9 @@ def node_V_capacity_check(which_load, low_lmt, high_lmt): vmax = np.max(v) vmin = np.min(v) cap_limit = new_kW - dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) - dss.run_command("Compile " + MasterFile) - dss.run_command("solve mode = snap") + dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) + dss.Text.Command("Compile " + MasterFile) + dss.Text.Command("solve mode = snap") else: tmp_kW = new_kW @@ -174,8 +174,8 @@ def node_V_capacity_check(which_load, low_lmt, high_lmt): ################################# for thermal overload capapcity ####################################################### -dss.run_command("solve mode = snap") -dss.run_command("export Overloads") +dss.Text.Command("solve mode = snap") +dss.Text.Command("export Overloads") # rename this overload csv file src = circuit.Name() + "_EXP_OVERLOADS.CSV" @@ -223,9 +223,9 @@ def thermal_overload_check(which_load, th_limit, case): while thermal_violation == False: # while the elements loadings are within limits keep on increasing new_kW = tmp_kW + add_kW # increase the load by 5 kW each time - dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) - dss.run_command("solve mode = snap") - dss.run_command("export Overloads") + dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) + dss.Text.Command("solve mode = snap") + dss.Text.Command("export Overloads") report = pd.read_csv(str(ckt_name) + "_EXP_OVERLOADS.CSV") report["Element"] = report["Element"].str.strip() @@ -262,12 +262,12 @@ def thermal_overload_check(which_load, th_limit, case): # print(thermal_violation) # for initial testing if thermal_violation == True: # print(thermal_violation) # for initial testing - dss.run_command("export Capacity") - dss.run_command("export Currents") + dss.Text.Command("export Capacity") + dss.Text.Command("export Currents") cap_limit = new_kW - dss.run_command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) - dss.run_command("Compile " + MasterFile) - dss.run_command("solve mode = snap") + dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) + dss.Text.Command("Compile " + MasterFile) + dss.Text.Command("solve mode = snap") else: tmp_kW = new_kW From 2a15fd544facbaa4a8f066ff38fa91bc10153f16 Mon Sep 17 00:00:00 2001 From: Daniel Thom Date: Wed, 28 Dec 2022 14:50:21 -0700 Subject: [PATCH 4/5] Reorganize code --- disco/cli/disco.py | 2 + disco/cli/ev.py | 160 +++++++ disco/ev/feeder_EV_HC.py | 723 ++++++++++++++++-------------- disco/ev/plot_hosting_capacity.py | 20 +- 4 files changed, 559 insertions(+), 346 deletions(-) diff --git a/disco/cli/disco.py b/disco/cli/disco.py index bbe823b4..26c35630 100644 --- a/disco/cli/disco.py +++ b/disco/cli/disco.py @@ -23,6 +23,7 @@ from disco.cli.config_generic_models import config_generic_models from disco.cli.upgrade_cost_analysis import upgrade_cost_analysis from disco.cli.pydss_hosting_capacity import hosting_capacity_by_timestep +from disco.cli.ev import ev logger = logging.getLogger(__name__) @@ -55,3 +56,4 @@ def cli(): cli.add_command(summarize_hosting_capacity) cli.add_command(upgrade_cost_analysis) cli.add_command(hosting_capacity_by_timestep) +cli.add_command(ev) diff --git a/disco/cli/ev.py b/disco/cli/ev.py index e69de29b..dfdfa3fc 100644 --- a/disco/cli/ev.py +++ b/disco/cli/ev.py @@ -0,0 +1,160 @@ +import fileinput +import logging +import os +import shutil +import sys +from pathlib import Path + +import click + +from jade.loggers import setup_logging +from jade.utils.utils import get_cli_string +from disco.ev.feeder_EV_HC import run +from disco import timer_stats_collector + + +logger = logging.getLogger(__name__) + + +@click.group() +def ev(): + """Run electic vehicle simulations.""" + + +@click.command() +@click.argument("master_file", type=click.Path(exists=True), callback=lambda *x: Path(x[2])) +@click.option( + "-l", "--lower-voltage-limit", default=0.95, type=float, help="Lower voltage limit (P.U.)" +) +@click.option( + "-u", "--upper-voltage-limit", default=1.05, type=float, help="Upper voltage limit (P.U.)" +) +@click.option( + "-v", + "--kw-step-voltage-violation", + default=10.0, + type=float, + show_default=True, + help="kW step value for detecting a voltage violation", +) +@click.option( + "-t", + "--kw-step-thermal-violation", + default=10.0, + type=float, + show_default=True, + help="kW step value for detecting a thermal violation", +) +@click.option( + "-e", + "--extra-percentages-for-existing-overloads", + default=(2.0,), + type=float, + multiple=True, + show_default=True, + help="Considers extra percentages for already overloaded elements", +) +@click.option( + "-T", + "--thermal-loading-limits", + default=(100.0,), + type=float, + multiple=True, + show_default=True, + help="Limits for thermal overloads", +) +@click.option( + "--export-circuit-elements/--no-export-circuit-elements", + default=False, + show_default=True, + help="Export rated values for all circuit elements", +) +# TODO: not implemented +# @click.option( +# "--plot-heatmap/--no-plot-heatmap", +# default=False, +# show_default=True, +# help="Plot heatmap for hosting capacity", +# ) +@click.option( + "-o", + "--output", + default="output", + show_default=True, + callback=lambda *x: Path(x[2]), + help="Output directory", +) +@click.option( + "-f", + "--force", + is_flag=True, + default=False, + show_default=True, + help="Overwrite output directory if it already exists.", +) +@click.option( + "--verbose", + is_flag=True, + default=False, + help="Enable debug logging", +) +def hosting_capacity( + master_file: Path, + lower_voltage_limit: float, + upper_voltage_limit: float, + kw_step_voltage_violation: float, + kw_step_thermal_violation: float, + extra_percentages_for_existing_overloads: tuple[float], + thermal_loading_limits: tuple[float], + export_circuit_elements: bool, + # plot_heatmap: bool, + output: Path, + force: bool, + verbose: bool, +): + """Compute hosting capacity for a feeder.""" + if output.exists(): + if force: + shutil.rmtree(output) + else: + print( + f"output directory {output} already exists. Choose a different path or pass --force.", + file=sys.stderr, + ) + sys.exit(1) + output.mkdir() + + level = logging.DEBUG if verbose else logging.INFO + filename = output / "ev_hosting_capacity.log" + logger = setup_logging( + "disco", filename, console_level=level, file_level=level, packages=["disco"] + ) + logger.info(get_cli_string()) + + backup_file = master_file.with_suffix(".bk") + shutil.copyfile(master_file, backup_file) + with fileinput.input(files=[master_file], inplace=True) as f: + for line in f: + if not line.strip().lower().startswith("solve"): + print(line, end="") + print("Solve mode=snapshot") + + try: + shutil.copyfile(master_file, output / "Master.dss") + run( + master_file, + lower_voltage_limit, + upper_voltage_limit, + kw_step_voltage_violation, + kw_step_thermal_violation, + extra_percentages_for_existing_overloads, + thermal_loading_limits, + export_circuit_elements, + output, + ) + finally: + os.rename(backup_file, master_file) + timer_stats_collector.log_stats(clear=True) + + +ev.add_command(hosting_capacity) diff --git a/disco/ev/feeder_EV_HC.py b/disco/ev/feeder_EV_HC.py index 5e268be1..e46931d9 100644 --- a/disco/ev/feeder_EV_HC.py +++ b/disco/ev/feeder_EV_HC.py @@ -1,8 +1,5 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Sep 25 07:53:52 2019 - -@author: ppaudyal +"""This code calculates the hosting capcity (checks voltage and thermal violation limits) for +a given feeder. """ # -*- coding: utf-8 -*- @@ -13,96 +10,284 @@ """ -# this code calculates the hosting capcity (checks voltage and thermal violation limits) for a given feeder -# every feeder model should have a 'Master.dss' as main dss file -# - +import logging import os +from pathlib import Path + import pandas as pd import numpy as np import math import opendssdirect as dss -import load_distance_from_SS -import plot_hosting_capacity -import number_of_ev_chargers - - -# variable quantities -feeder_folder = "477NHS" #'BaseCaseWRR074' #'BaseaseALD095' -low_lmt_0 = 0.950 # lower voltage limit -high_lmt_0 = 1.050 # upper voltage limit -low_lmt_1 = 0.975 # lower voltage limit -high_lmt_1 = 1.05 # upper voltage limit -low_lmt_2 = 0.985 # lower voltage limit -high_lmt_2 = 1.05 # upper voltage limit -v_limit = [[low_lmt_0, high_lmt_0], [low_lmt_1, high_lmt_1], [low_lmt_2, high_lmt_2]] -step_v = 10 -step_thermal = 5 - -extra_percent = [2] # , 25] # in case, this considers extra % for already overloaded elements -extra_limit = [ - 100 -] # , 120] #120, if we want to consider overload condition when %Normal>= 120 or any other % - -master_dssfile = "Master.dss" - -# **************** System Information ***************** -MainDir = os.getcwd() -FeederDir = os.path.join(MainDir, feeder_folder) -MasterFile = os.path.join(FeederDir, master_dssfile) -dss.Text.Command("clear") -dss.Text.Command(f"Compile {MasterFile}") -dss.Text.Command("solve") -circuit = dss.Circuit -ckt_name = circuit.Name() -print("\n*** Circuit name:", ckt_name, "***\n from ", feeder_folder, " folder\n") - -AllNodeNames = circuit.YNodeOrder() -pd.DataFrame(AllNodeNames).to_csv("Allnodenames.csv") -# print(type(AllNodeNames)) - -# --------- Voltage Base kV Information ----- -node_number = len(AllNodeNames) -Vbase_allnode = [0] * node_number -ii = 0 -for node in AllNodeNames: - circuit.SetActiveBus(node) - Vbase_allnode[ii] = dss.Bus.kVBase() * 1000 - ii = ii + 1 - -# --------- Capacitor Information ---------- -capNames = dss.Capacitors.AllNames() -hCapNames = [None] * len(capNames) -for i, n in enumerate(capNames): - hCapNames[i] = str(n) -hCapNames = str(hCapNames)[1:-1] - -# --------- Regulator Information ---------- -regNames = dss.RegControls.AllNames() -hRegNames = [None] * len(regNames) -for i, n in enumerate(regNames): - hRegNames[i] = str(n) -hRegNames = str(hRegNames)[1:-1] - -dss.Text.Command("solve mode=snap") -# dss.Text.Command('export voltages') -v = np.array(dss.Circuit.AllBusMagPu()) - - -dss.utils.lines_to_dataframe() - -volt_violation = np.any(v > high_lmt_0) or np.any(v < low_lmt_0) -print("Initial condition voltage violation: ", volt_violation) - -Bus_Distance = [] -for node in AllNodeNames: - circuit.SetActiveBus(node) - # print('node', node, 'bus_distance', dss.Bus.Distance()) # testing - Bus_Distance.append(dss.Bus.Distance()) -np.savetxt("node_distance_" + str(feeder_folder) + ".csv", Bus_Distance) +from jade.utils.timing_utils import track_timing, Timer +from disco import timer_stats_collector +from .load_distance_from_SS import calc_dist +from .plot_hosting_capacity import ( + plot_capacity_V, + plot_capacity_thermal_1, + plot_capacity_thermal_2, +) +from .number_of_ev_chargers import levels_of_charger + + +logger = logging.getLogger(__name__) + + +@track_timing(timer_stats_collector) +def run( + master_file: Path, + lower_voltage_limit: float, + upper_voltage_limit: float, + kw_step_voltage_violation: float, + kw_step_thermal_violation: float, + extra_percentages_for_existing_overloads: tuple[float], + thermal_loading_limits: tuple[float], + export_circuit_elements: bool, + # plot_heatmap: bool, + output_dir: Path, +): + dss.Text.Command("clear") + compile_circuit(master_file) + circuit = dss.Circuit + ckt_name = circuit.Name() + logger.info(" Circuit name: %s from %s folder", ckt_name, master_file) + + if export_circuit_elements: + export_circuit_element_properties(output_dir) + + v_limit = [(lower_voltage_limit, upper_voltage_limit)] + # TODO: do we need more limits? + # assert lower_voltage_limit == 0.95 + # v_limit += [(0.975, upper_voltage_limit), (0.985, upper_voltage_limit)] + # TODO: Priti, how to calculate these programmatically? + + AllNodeNames = circuit.YNodeOrder() + + # --------- Voltage Base kV Information ----- + node_number = len(AllNodeNames) + Vbase_allnode = [0] * node_number + ii = 0 + for node in AllNodeNames: + circuit.SetActiveBus(node) + Vbase_allnode[ii] = dss.Bus.kVBase() * 1000 + ii = ii + 1 + + # --------- Capacitor Information ---------- + capNames = dss.Capacitors.AllNames() + hCapNames = [None] * len(capNames) + for i, n in enumerate(capNames): + hCapNames[i] = str(n) + hCapNames = str(hCapNames)[1:-1] + + # --------- Regulator Information ---------- + regNames = dss.RegControls.AllNames() + hRegNames = [None] * len(regNames) + for i, n in enumerate(regNames): + hRegNames[i] = str(n) + hRegNames = str(hRegNames)[1:-1] + + dss.Text.Command("solve mode=snap") + # dss.Text.Command('export voltages') + v = np.array(dss.Circuit.AllBusMagPu()) + + volt_violation = np.any(v > upper_voltage_limit) or np.any(v < lower_voltage_limit) + logger.info("Initial condition voltage violation: %s", volt_violation) + + Bus_Distance = [] + for node in AllNodeNames: + circuit.SetActiveBus(node) + Bus_Distance.append(dss.Bus.Distance()) + np.savetxt(output_dir / "node_distance.csv", Bus_Distance) + + ####### for thermal overload capapcity ####################################################### + dss.Text.Command("solve mode = snap") + overloads_filename = output_dir / "overloads.csv" + dss.Text.Command(f"export Overloads {overloads_filename}") + + # read this overload file and record the ' %Normal' for each line in this file + overload_df = pd.read_csv(overloads_filename) + # len(overload_df) + if len(overload_df) == 0: + logger.info("No thermal violation initially") + else: + logger.info("Thermal violation exists initially") + logger.info("Overloads: %s", overload_df.head()) + + elements = [[], []] + amps = [[], []] + new_threshold = [[], []] + for j in range(len(extra_percentages_for_existing_overloads)): + for i in range(len(overload_df)): + overload_df["Element"] = overload_df["Element"].str.strip() + element = overload_df["Element"][i] + amp = overload_df[" %Normal"][i] + element_new_limit = amp + extra_percentages_for_existing_overloads[j] + elements[j].append(str(element)) + amps[j].append(amp) + new_threshold[j].append(element_new_limit) + + logger.debug("new_threshold=%s", new_threshold) + + ############################################################################################## + # get the load data + [Load, totalkW] = get_loads(dss, circuit, 0, "") + + # calculate the hosting capacity based on voltage constraints + v_output_df = [] + for j in range(len(v_limit)): + lmt1 = v_limit[j][0] + lmt2 = v_limit[j][1] + v_lst = [] + v_output_list = [] + v_threshold = [] + v_allow_limit = [] + v_names = [] + v_bus_name = [] + v_default_load = [] + v_maxv = [] + v_minv = [] + for i in range(len(Load)): + logger.info("Run voltage check on load index=%s", i) + v_node_cap = node_V_capacity_check( + master_file, Load[i], lmt1, lmt2, kw_step_voltage_violation + ) # v_node_cap is a list + v_allowable_load = v_node_cap[0] - kw_step_voltage_violation + v_threshold.append(v_node_cap[0]) + v_allow_limit.append(v_allowable_load) + v_names.append(Load[i]["name"]) + v_bus_name.append(Load[i]["bus1"]) + v_default_load.append(Load[i]["kW"]) + v_maxv.append(v_node_cap[1]) + v_minv.append(v_node_cap[2]) + + v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] + v_output_list = list(map(list, zip(*v_lst))) + # v_output_df[0] = pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage']) + v_output_df.append( + pd.DataFrame( + v_output_list, + columns=[ + "Load", + "Bus", + "Initial_kW", + "Volt_Violation", + "Maximum_kW", + "Max_voltage", + "Min_voltage", + ], + ) + ) + + v_output_df[j].to_csv(output_dir / f"Hosting_capacity_test_{lmt1}.csv") + + # calculate the hosting capacity based on thermal ratings constraint + th_threshold = [[], []] + th_allow_limit = [[], []] + th_names = [[], []] + th_bus_name = [[], []] + th_default_load = [[], []] + + th_output_df = [] + for i in range(len(thermal_loading_limits)): + logger.info("Run thermal check on index=%s limit=%s", i, thermal_loading_limits[i]) + th_lst = [] + for j in range(len(Load)): + if j == 0: + continue + th_node_cap = thermal_overload_check( + Load[j], thermal_loading_limits[i], i, kw_step_thermal_violation, output_dir + ) # th_node_cap is a list + th_allowable_load = ( + th_node_cap[0] - kw_step_thermal_violation + ) # 5 # th_node_cap[0] is a float # reduce the value of add_kW + th_threshold[i].append(th_node_cap[0]) + th_allow_limit[i].append(th_allowable_load) + th_names[i].append(Load[j]["name"]) + th_bus_name[i].append(Load[j]["bus1"]) + th_default_load[i].append(Load[j]["kW"]) + + th_lst = [ + th_names[i], + th_bus_name[i], + th_default_load[i], + th_threshold[i], + th_allow_limit[i], + ] + th_output_list = list(map(list, zip(*th_lst))) + th_output_df.append( + pd.DataFrame( + th_output_list, + columns=["Load", "Bus", "Initial_kW", "Thermal_Violation", "Maximum_kW"], + ) + ) + + th_output_df[i].to_csv( + output_dir / f"Thermal_capacity_test_{thermal_loading_limits[i]}.csv" + ) + + #################### for plotting results #################################################### + + load_bus = pd.DataFrame() + load_bus["Load"] = th_output_df[0]["Load"] # + load_bus["Bus"] = th_output_df[0]["Bus"] + node_distance = pd.DataFrame() + node_distance["Node"] = AllNodeNames + node_distance["Distance"] = Bus_Distance + + dist_file = calc_dist(load_bus, node_distance) + + dist_file["Initial_MW"] = th_output_df[0]["Initial_kW"] / 1000 + dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"] / 1000 + dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"] / 1000 + dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"] / 1000 + + dist_file["Thermal_Violation_100"] = th_output_df[0]["Thermal_Violation"] / 1000 + # dist_file["Thermal_Violation_120"] = th_output_df[1]["Thermal_Violation"]/1000 + + plot_df = dist_file.sort_values(by=["Distance"]) + + # plot voltage violation scenarios + plot_capacity_V( + plot_df, + "Initial_MW", + "Volt_Violation_0.95", + "Volt_Violation_0.975", + "Volt_Violation_0.985", + output_dir, + ) + + # plot thermal violation + plot_capacity_thermal_1(plot_df, "Initial_MW", "Thermal_Violation_100", output_dir, 100) + # plot_capacity_thermal_2(plot_df, 'Initial_MW', 'Thermal_Violation_100', 'Thermal_Violation_120', output_dir) + + ### Assuming the hosting capacity is limited by thermal loading ############################## + + ## Difference of initial load and maximum hosting capacity (assuming always thermal limit occurs first) ## + + for i in range(len(thermal_loading_limits)): + diff = th_output_df[i]["Thermal_Violation"] - th_output_df[i]["Initial_kW"] + new_df = pd.DataFrame() + + new_df["Load"] = th_output_df[i]["Load"] + new_df["Bus"] = th_output_df[i]["Bus"] + new_df["Initial_kW"] = th_output_df[i]["Initial_kW"] + new_df["Hosting_capacity(kW)"] = diff # additional load it can support + + new_df.to_csv( + output_dir / "Additional_HostingCapacity_" + str(thermal_loading_limits[i]) + ".csv", + index=False, + ) + + ## find number of ev chargers for each node ## + + chargers_3_2_1 = levels_of_charger(th_output_df[i]) + chargers_3_2_1.to_csv( + output_dir / "Loadwithlevel3_2_1_" + str(thermal_loading_limits[i]) + ".csv" + ) + # --------- Load Information ---------- +@track_timing(timer_stats_collector) def get_loads(dss, circuit, loadshape_flag, loadshape_folder): data = [] load_flag = dss.Loads.First() @@ -142,283 +327,143 @@ def get_loads(dss, circuit, loadshape_flag, loadshape_folder): return [data, total_load] -############################################## for voltage violation capacity ###################################### -def node_V_capacity_check(which_load, low_lmt, high_lmt): +############################ for voltage violation capacity ###################################### +@track_timing(timer_stats_collector) +def node_V_capacity_check(master_file, which_load, low_lmt, high_lmt, kw_step_voltage_violation): initial_kW = which_load["kW"] - # print(initial_kW) #just for testing - add_kW = step_v # 10 + logger.debug("initial_kW=%s", initial_kW) tmp_kW = initial_kW volt_violation = False v = None - while volt_violation == False: - # while the voltages are within limits keep on increasing - new_kW = tmp_kW + add_kW # increase the load by 10 kW each time - dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) - dss.Text.Command("solve mode = snap") - v = np.array(dss.Circuit.AllBusMagPu()) - volt_violation = np.any(v > high_lmt) or np.any(v < low_lmt) - # print(volt_violation) # for initial testing - if volt_violation == True: - # print(volt_violation) # for initial testing - vmax = np.max(v) - vmin = np.min(v) - cap_limit = new_kW - dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) - dss.Text.Command("Compile " + MasterFile) + while not volt_violation: + with Timer(timer_stats_collector, "check_voltage_violation_after_increasing_kw"): + # while the voltages are within limits keep on increasing + new_kW = tmp_kW + kw_step_voltage_violation + logger.info("Check voltage violation with kW=%s", new_kW) + dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) dss.Text.Command("solve mode = snap") - else: - tmp_kW = new_kW + v = np.array(dss.Circuit.AllBusMagPu()) + volt_violation = np.any(v > high_lmt) or np.any(v < low_lmt) + logger.debug("Voltage violation = %s", volt_violation) + if volt_violation: + vmax = np.max(v) + vmin = np.min(v) + cap_limit = new_kW + # TODO: no point in doing this, right? + # dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) + compile_circuit(master_file) + else: + tmp_kW = new_kW return [cap_limit, vmax, vmin] -################################# for thermal overload capapcity ####################################################### -dss.Text.Command("solve mode = snap") -dss.Text.Command("export Overloads") - -# rename this overload csv file -src = circuit.Name() + "_EXP_OVERLOADS.CSV" -dst = feeder_folder + "_overloads.csv" -if os.path.exists(dst): - os.remove(dst) -os.rename(src, dst) - -# read this overload file and record the ' %Normal' for each line in this file -overload_df = pd.read_csv(dst) -# len(overload_df) -if len(overload_df) == 0: - print("No thermal violation initially \n") -else: - print("Thermal violation exists initially \n") - print(overload_df.head()) - - -elements = [[], []] -amps = [[], []] -new_threshold = [[], []] -for j in range(len(extra_percent)): - - for i in range(len(overload_df)): - overload_df["Element"] = overload_df[ - "Element" - ].str.strip() # removing the extra spaces at the end (if any) - element = overload_df["Element"][i] - amp = overload_df[" %Normal"][i] - element_new_limit = amp + extra_percent[j] - elements[j].append(str(element)) - amps[j].append(amp) - new_threshold[j].append(element_new_limit) - -print("new_threshold", new_threshold) # testing - - -def thermal_overload_check(which_load, th_limit, case): +@track_timing(timer_stats_collector) +def thermal_overload_check( + which_load, th_limit, case, kw_step_thermal_violation, output_dir: Path +): initial_kW = which_load["kW"] - # print(initial_kW) - add_kW = step_thermal # 5 + logger.debug("initial_kW=%s", initial_kW) tmp_kW = initial_kW thermal_violation = False - while thermal_violation == False: - # while the elements loadings are within limits keep on increasing - new_kW = tmp_kW + add_kW # increase the load by 5 kW each time - dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) - dss.Text.Command("solve mode = snap") - dss.Text.Command("export Overloads") - report = pd.read_csv(str(ckt_name) + "_EXP_OVERLOADS.CSV") - report["Element"] = report["Element"].str.strip() - - if len(report) == 0: # if no any overload element - thermal_violation = False - - elif report["Element"].isin(elements[case]).any(): - for i in range(len(report)): - if report["Element"][i] in elements[case]: - indx_ = elements[case].index(report["Element"][i]) # find the index of - if report[" %Normal"][i] >= new_threshold[case][indx_]: - thermal_violation = True # just exit here (get out of for loop) - # print(i) - break + while not thermal_violation: + with Timer(timer_stats_collector, "check_thermal_violation_after_increasing_kw"): + # while the elements loadings are within limits keep on increasing + new_kW = tmp_kW + kw_step_thermal_violation + logger.info("Check thermal violation with kW=%s", new_kW) + dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) + dss.Text.Command("solve mode = snap") + overloads_filename = output_dir / "overloads.csv" + dss.Text.Command("export Overloads {overloads_filename}") + report = pd.read_csv(overloads_filename) + report["Element"] = report["Element"].str.strip() + + if len(report) == 0: # if no any overload element + thermal_violation = False + + elif report["Element"].isin(elements[case]).any(): + for i in range(len(report)): + if report["Element"][i] in elements[case]: + indx_ = elements[case].index(report["Element"][i]) # find the index of + if report[" %Normal"][i] >= new_threshold[case][indx_]: + thermal_violation = True # just exit here (get out of for loop) + break + else: + thermal_violation = False else: - thermal_violation = False - else: - # check the percentage normal if greater than specified % then only violation + # check the percentage normal if greater than specified % then only violation + if report[" %Normal"][i] >= th_limit: + thermal_violation = True + break + else: + thermal_violation = False + + else: + # check the percentage normal if greater than specified % then only violation + for i in range(len(report)): if report[" %Normal"][i] >= th_limit: thermal_violation = True break else: thermal_violation = False - else: - # check the percentage normal if greater than specified % then only violation - for i in range(len(report)): - if report[" %Normal"][i] >= th_limit: - thermal_violation = True - break - else: - thermal_violation = False - - # print(thermal_violation) # for initial testing - if thermal_violation == True: - # print(thermal_violation) # for initial testing - dss.Text.Command("export Capacity") - dss.Text.Command("export Currents") - cap_limit = new_kW - dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) - dss.Text.Command("Compile " + MasterFile) - dss.Text.Command("solve mode = snap") - else: - tmp_kW = new_kW + logger.debug("thermal_violation=%s", thermal_violation) + if thermal_violation: + logger.debug("thermal_violation=%s", thermal_violation) + dss.Text.Command(f"export Capacity {output_dir / 'capacity.csv'}") + dss.Text.Command(f"export Currents {output_dir / 'currents.csv'}") + cap_limit = new_kW + # TODO: no point in doing this, right? + # dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) + compile_circuit(master_file) + else: + tmp_kW = new_kW return [cap_limit] -######################################################################################################################### -# get the load data -[Load, totalkW] = get_loads(dss, circuit, 0, "") - -# calculate the hosting capacity based on voltage constraints -v_output_df = [] -for j in range(len(v_limit)): - lmt1 = v_limit[j][0] - lmt2 = v_limit[j][1] - v_lst = [] - v_output_list = [] - v_threshold = [] - v_allow_limit = [] - v_names = [] - v_bus_name = [] - v_default_load = [] - v_maxv = [] - v_minv = [] - for i in range(len(Load)): - v_node_cap = node_V_capacity_check(Load[i], lmt1, lmt2) # v_node_cap is a list - v_allowable_load = v_node_cap[0] - step_v # 10 # v_node_cap[0] is a float - v_threshold.append(v_node_cap[0]) - v_allow_limit.append(v_allowable_load) - v_names.append(Load[i]["name"]) - v_bus_name.append(Load[i]["bus1"]) - v_default_load.append(Load[i]["kW"]) - v_maxv.append(v_node_cap[1]) - v_minv.append(v_node_cap[2]) - - v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] - v_output_list = list(map(list, zip(*v_lst))) - # v_output_df[0] = pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage']) - v_output_df.append( - pd.DataFrame( - v_output_list, - columns=[ - "Load", - "Bus", - "Initial_kW", - "Volt_Violation", - "Maximum_kW", - "Max_voltage", - "Min_voltage", - ], - ) - ) - - v_output_df[j].to_csv("Hosting_capacity_test_" + str(feeder_folder) + "_" + str(lmt1) + ".csv") - -# calculate the hosting capacity based on thermal ratings constraint -th_threshold = [[], []] -th_allow_limit = [[], []] -th_names = [[], []] -th_bus_name = [[], []] -th_default_load = [[], []] - -th_output_df = [] -for i in range(len(extra_limit)): - - th_lst = [] - for j in range(len(Load)): - th_node_cap = thermal_overload_check(Load[j], extra_limit[i], i) # th_node_cap is a list - th_allowable_load = ( - th_node_cap[0] - step_thermal - ) # 5 # th_node_cap[0] is a float # reduce the value of add_kW - th_threshold[i].append(th_node_cap[0]) - th_allow_limit[i].append(th_allowable_load) - th_names[i].append(Load[j]["name"]) - th_bus_name[i].append(Load[j]["bus1"]) - th_default_load[i].append(Load[j]["kW"]) - - th_lst = [th_names[i], th_bus_name[i], th_default_load[i], th_threshold[i], th_allow_limit[i]] - th_output_list = list(map(list, zip(*th_lst))) - th_output_df.append( - pd.DataFrame( - th_output_list, - columns=["Load", "Bus", "Initial_kW", "Thermal_Violation", "Maximum_kW"], - ) - ) - - th_output_df[i].to_csv( - "Thermal_capacity_test_" + str(feeder_folder) + "_" + str(extra_limit[i]) + ".csv" - ) - - -############################################## for plotting results #################################################### - -load_bus = pd.DataFrame() -load_bus["Load"] = th_output_df[0]["Load"] # -load_bus["Bus"] = th_output_df[0]["Bus"] -node_distance = pd.DataFrame() -node_distance["Node"] = AllNodeNames -node_distance["Distance"] = Bus_Distance - -dist_file = load_distance_from_SS.calc_dist(load_bus, node_distance) - - -dist_file["Initial_MW"] = th_output_df[0]["Initial_kW"] / 1000 -dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"] / 1000 -dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"] / 1000 -dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"] / 1000 - -dist_file["Thermal_Violation_100"] = th_output_df[0]["Thermal_Violation"] / 1000 -# dist_file["Thermal_Violation_120"] = th_output_df[1]["Thermal_Violation"]/1000 - -plot_df = dist_file.sort_values(by=["Distance"]) - -# plot voltage violation scenarios -plot_hosting_capacity.plot_capacity_V( - plot_df, - "Initial_MW", - "Volt_Violation_0.95", - "Volt_Violation_0.975", - "Volt_Violation_0.985", - feeder_folder, -) - -# plot thermal violation -plot_hosting_capacity.plot_capacity_thermal_1( - plot_df, "Initial_MW", "Thermal_Violation_100", feeder_folder, 100 -) -# plot_hosting_capacity.plot_capacity_thermal_2(plot_df, 'Initial_MW', 'Thermal_Violation_100', 'Thermal_Violation_120', feeder_folder) - - -####################### Assuming the hosting capacity is limited by thermal loading ##################################### - -## Difference of initial load and maximum hosting capacity (assuming always thermal limit occurs first) ## - -for i in range(len(extra_limit)): - diff = th_output_df[i]["Thermal_Violation"] - th_output_df[i]["Initial_kW"] - new_df = pd.DataFrame() - - new_df["Load"] = th_output_df[i]["Load"] - new_df["Bus"] = th_output_df[i]["Bus"] - new_df["Initial_kW"] = th_output_df[i]["Initial_kW"] - new_df["Hosting_capacity(kW)"] = diff # additional load it can support - - new_df.to_csv( - str(feeder_folder) + "_Additional_HostingCapacity_" + str(extra_limit[i]) + ".csv", - index=False, +@track_timing(timer_stats_collector) +def compile_circuit(master_file: Path): + logger.info("Compile circuit from %s", master_file) + orig = os.getcwd() + try: + dss.Text.Command(f"Compile {master_file}") + finally: + os.chdir(orig) + + +@track_timing(timer_stats_collector) +def export_circuit_element_properties(output_dir: Path): + export_dir = output_dir / "circuit_elements" + export_dir.mkdir() + exports = ( + ("Capacitor", dss.Capacitors.Count), + ("Fuse", dss.Fuses.Count), + ("Generator", dss.Generators.Count), + ("Isource", dss.Isource.Count), + ("Line", dss.Lines.Count), + ("Load", dss.Loads.Count), + ("Monitor", dss.Monitors.Count), + ("PVSystem", dss.PVsystems.Count), + ("Recloser", dss.Reclosers.Count), + ("RegControl", dss.RegControls.Count), + ("Relay", dss.Relays.Count), + ("Sensor", dss.Sensors.Count), + ("Transformer", dss.Transformers.Count), + ("Vsource", dss.Vsources.Count), + ("XYCurve", dss.XYCurves.Count), ) - ## find number of ev chargers for each node ## + for class_name, count_func in exports: + if count_func() > 0: + filename = export_dir / +f"{class_name}.csv" + df = dss.utils.class_to_dataframe(class_name) + df.to_csv(filename) + logger.info("Exported %s information to %s.", class_name, filename) + else: + logger.info("There are no elements of type=%s to export", class_name) - chargers_3_2_1 = number_of_ev_chargers.levels_of_charger(th_output_df[i]) - chargers_3_2_1.to_csv( - str(feeder_folder) + "_Loadwithlevel3_2_1_" + str(extra_limit[i]) + ".csv" - ) + all_node_names = circuit.YNodeOrder() + pd.DataFrame(all_node_names).to_csv(output_dir / "Allnodenames.csv") diff --git a/disco/ev/plot_hosting_capacity.py b/disco/ev/plot_hosting_capacity.py index 28f7e1e6..78c93d37 100644 --- a/disco/ev/plot_hosting_capacity.py +++ b/disco/ev/plot_hosting_capacity.py @@ -4,11 +4,17 @@ @author: ppaudyal """ -import matplotlib.pyplot as plt +import logging import os +from pathlib import Path + +import matplotlib.pyplot as plt + + +logger = logging.getLogger(__name__) -def plot_capacity_V(file, caseA, caseB, caseC, caseD, feeder): +def plot_capacity_V(file, caseA, caseB, caseC, caseD, output_dir: Path): fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) # figsize=(18,12) vplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") vplot2 = ax1.plot( @@ -44,10 +50,10 @@ def plot_capacity_V(file, caseA, caseB, caseC, caseD, feeder): "Feeder Node (where highest is furthest from substation)", fontsize=33 ) ##'Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)' # plt.tight_layout() - plt.savefig(str(feeder) + "_Cap_by_V_limit.png", bbox_inches="tight", dpi=150) + plt.savefig(output_dir / "Cap_by_V_limit.png", bbox_inches="tight", dpi=150) -def plot_capacity_thermal_2(file, caseA, caseB, caseC, feeder): +def plot_capacity_thermal_2(file, caseA, caseB, caseC, output_dir: Path): fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) thplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") thplot2 = ax1.plot( @@ -75,10 +81,10 @@ def plot_capacity_thermal_2(file, caseA, caseB, caseC, feeder): "Feeder Node (where highest is furthest from substation)", fontsize=30 ) # ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) plt.tight_layout() - plt.savefig(str(feeder) + "_Cap_by_thermal_limit.png", dpi=150) + plt.savefig(output_dir / "Cap_by_thermal_limit.png", dpi=150) -def plot_capacity_thermal_1(file, caseA, caseB, feeder, extra): +def plot_capacity_thermal_1(file, caseA, caseB, output_dir: Path, extra): fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) thplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") thplot2 = ax1.plot( @@ -96,7 +102,7 @@ def plot_capacity_thermal_1(file, caseA, caseB, feeder, extra): ) # ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) # fig.tight_layout() plt.savefig( - str(feeder) + "_Cap_by_thermal_limit_" + str(extra) + ".png", bbox_inches="tight", dpi=150 + output_dir / "Cap_by_thermal_limit_" + str(extra) + ".png", bbox_inches="tight", dpi=150 ) From ee6ebc22e37b4dd0e038c188f45a42733098dcc5 Mon Sep 17 00:00:00 2001 From: Daniel Thom Date: Thu, 15 Jun 2023 11:56:03 -0600 Subject: [PATCH 5/5] Make EV hosting capacity code run faster --- disco/cli/ev.py | 89 ++-- disco/ev/feeder_EV_HC.py | 722 +++++++++++++++++------------- disco/ev/plot_hosting_capacity.py | 50 +-- 3 files changed, 488 insertions(+), 373 deletions(-) diff --git a/disco/cli/ev.py b/disco/cli/ev.py index dfdfa3fc..4533a5be 100644 --- a/disco/cli/ev.py +++ b/disco/cli/ev.py @@ -1,4 +1,3 @@ -import fileinput import logging import os import shutil @@ -10,7 +9,6 @@ from jade.loggers import setup_logging from jade.utils.utils import get_cli_string from disco.ev.feeder_EV_HC import run -from disco import timer_stats_collector logger = logging.getLogger(__name__) @@ -24,10 +22,13 @@ def ev(): @click.command() @click.argument("master_file", type=click.Path(exists=True), callback=lambda *x: Path(x[2])) @click.option( - "-l", "--lower-voltage-limit", default=0.95, type=float, help="Lower voltage limit (P.U.)" + "-c", "--num-cpus", type=int, help="Number of CPUs to use, default is all." ) @click.option( - "-u", "--upper-voltage-limit", default=1.05, type=float, help="Upper voltage limit (P.U.)" + "-l", "--lower-voltage-limit", default=0.95, show_default=True, type=float, help="Lower voltage limit (P.U.)" +) +@click.option( + "-u", "--upper-voltage-limit", default=1.05, show_default=True, type=float, help="Upper voltage limit (P.U.)" ) @click.option( "-v", @@ -47,21 +48,29 @@ def ev(): ) @click.option( "-e", - "--extra-percentages-for-existing-overloads", - default=(2.0,), + "--extra-percentage-for-existing-overloads", + default=2.0, type=float, - multiple=True, show_default=True, - help="Considers extra percentages for already overloaded elements", + help="Considers extra percentage for already overloaded elements", ) @click.option( "-T", - "--thermal-loading-limits", - default=(100.0,), + "--thermal-loading-limit", + default=100.0, type=float, - multiple=True, show_default=True, - help="Limits for thermal overloads", + help="Limit for thermal overloads", +) +@click.option( + "--thermal-tolerance", + type=float, + help="Tolerance to use when finding lowest thermal violations. Uses --kw-step-thermal-violation by default", +) +@click.option( + "--voltage-tolerance", + type=float, + help="Tolerance to use when finding lowest voltage violations. Uses --kw-step-voltage-violation by default", ) @click.option( "--export-circuit-elements/--no-export-circuit-elements", @@ -85,8 +94,7 @@ def ev(): help="Output directory", ) @click.option( - "-f", - "--force", + "--overwrite", is_flag=True, default=False, show_default=True, @@ -100,29 +108,34 @@ def ev(): ) def hosting_capacity( master_file: Path, + num_cpus: int | None, lower_voltage_limit: float, upper_voltage_limit: float, kw_step_voltage_violation: float, kw_step_thermal_violation: float, - extra_percentages_for_existing_overloads: tuple[float], - thermal_loading_limits: tuple[float], + extra_percentage_for_existing_overloads: float, + thermal_loading_limit: float, + voltage_tolerance: float, + thermal_tolerance: float, export_circuit_elements: bool, # plot_heatmap: bool, output: Path, - force: bool, + overwrite: bool, verbose: bool, ): """Compute hosting capacity for a feeder.""" if output.exists(): - if force: + if overwrite: shutil.rmtree(output) else: print( - f"output directory {output} already exists. Choose a different path or pass --force.", + f"{output} already exists. Choose a different path or pass --overwrite.", file=sys.stderr, ) sys.exit(1) output.mkdir() + thermal_tolerance = thermal_tolerance or kw_step_thermal_violation + voltage_tolerance = voltage_tolerance or kw_step_voltage_violation level = logging.DEBUG if verbose else logging.INFO filename = output / "ev_hosting_capacity.log" @@ -131,30 +144,20 @@ def hosting_capacity( ) logger.info(get_cli_string()) - backup_file = master_file.with_suffix(".bk") - shutil.copyfile(master_file, backup_file) - with fileinput.input(files=[master_file], inplace=True) as f: - for line in f: - if not line.strip().lower().startswith("solve"): - print(line, end="") - print("Solve mode=snapshot") - - try: - shutil.copyfile(master_file, output / "Master.dss") - run( - master_file, - lower_voltage_limit, - upper_voltage_limit, - kw_step_voltage_violation, - kw_step_thermal_violation, - extra_percentages_for_existing_overloads, - thermal_loading_limits, - export_circuit_elements, - output, - ) - finally: - os.rename(backup_file, master_file) - timer_stats_collector.log_stats(clear=True) + run( + master_file=master_file, + lower_voltage_limit=lower_voltage_limit, + upper_voltage_limit=upper_voltage_limit, + kw_step_voltage_violation=kw_step_voltage_violation, + voltage_tolerance=voltage_tolerance, + kw_step_thermal_violation=kw_step_thermal_violation, + thermal_tolerance=thermal_tolerance, + extra_percentage_for_existing_overloads=extra_percentage_for_existing_overloads, + thermal_loading_limit=thermal_loading_limit, + export_circuit_elements=export_circuit_elements, + output_dir=output, + num_cpus=num_cpus, + ) ev.add_command(hosting_capacity) diff --git a/disco/ev/feeder_EV_HC.py b/disco/ev/feeder_EV_HC.py index e46931d9..dce6c6d9 100644 --- a/disco/ev/feeder_EV_HC.py +++ b/disco/ev/feeder_EV_HC.py @@ -9,9 +9,13 @@ @author: ppaudyal """ - +import fileinput +import itertools import logging import os +import shutil +import sys +from concurrent.futures import ProcessPoolExecutor from pathlib import Path import pandas as pd @@ -19,8 +23,6 @@ import math import opendssdirect as dss -from jade.utils.timing_utils import track_timing, Timer -from disco import timer_stats_collector from .load_distance_from_SS import calc_dist from .plot_hosting_capacity import ( plot_capacity_V, @@ -33,216 +35,151 @@ logger = logging.getLogger(__name__) -@track_timing(timer_stats_collector) def run( master_file: Path, lower_voltage_limit: float, upper_voltage_limit: float, kw_step_voltage_violation: float, + voltage_tolerance: float, kw_step_thermal_violation: float, - extra_percentages_for_existing_overloads: tuple[float], - thermal_loading_limits: tuple[float], + thermal_tolerance: float, + extra_percentage_for_existing_overloads: float, + thermal_loading_limit: float, export_circuit_elements: bool, # plot_heatmap: bool, output_dir: Path, + num_cpus=None, +): + # This accounts for master files that enable time series mode. + backup_file = master_file.with_suffix(".bk") + shutil.copyfile(master_file, backup_file) + with fileinput.input(files=[master_file], inplace=True) as f: + for line in f: + if not line.strip().lower().startswith("solve"): + print(line, end="") + print("Solve mode=snapshot") + + try: + shutil.copyfile(master_file, output_dir / "Master.dss") + _run( + master_file=master_file, + lower_voltage_limit=lower_voltage_limit, + upper_voltage_limit=upper_voltage_limit, + kw_step_voltage_violation=kw_step_voltage_violation, + voltage_tolerance=voltage_tolerance, + kw_step_thermal_violation=kw_step_thermal_violation, + thermal_tolerance=thermal_tolerance, + extra_percentage_for_existing_overloads=extra_percentage_for_existing_overloads, + thermal_loading_limit=thermal_loading_limit, + export_circuit_elements=export_circuit_elements, + output_dir=output_dir, + num_cpus=num_cpus, + ) + finally: + os.rename(backup_file, master_file) + + +def _run( + master_file: Path, + lower_voltage_limit: float, + upper_voltage_limit: float, + kw_step_voltage_violation: float, + voltage_tolerance: float, + kw_step_thermal_violation: float, + thermal_tolerance: float, + extra_percentage_for_existing_overloads: float, + thermal_loading_limit: float, + export_circuit_elements: bool, + # plot_heatmap: bool, + output_dir: Path, + num_cpus=None, ): dss.Text.Command("clear") compile_circuit(master_file) - circuit = dss.Circuit - ckt_name = circuit.Name() - logger.info(" Circuit name: %s from %s folder", ckt_name, master_file) + ckt_name = dss.Circuit.Name() + logger.info("Circuit name: %s from %s master_file", ckt_name, master_file) if export_circuit_elements: export_circuit_element_properties(output_dir) - v_limit = [(lower_voltage_limit, upper_voltage_limit)] - # TODO: do we need more limits? - # assert lower_voltage_limit == 0.95 - # v_limit += [(0.975, upper_voltage_limit), (0.985, upper_voltage_limit)] - # TODO: Priti, how to calculate these programmatically? - - AllNodeNames = circuit.YNodeOrder() - - # --------- Voltage Base kV Information ----- - node_number = len(AllNodeNames) - Vbase_allnode = [0] * node_number - ii = 0 - for node in AllNodeNames: - circuit.SetActiveBus(node) - Vbase_allnode[ii] = dss.Bus.kVBase() * 1000 - ii = ii + 1 - - # --------- Capacitor Information ---------- - capNames = dss.Capacitors.AllNames() - hCapNames = [None] * len(capNames) - for i, n in enumerate(capNames): - hCapNames[i] = str(n) - hCapNames = str(hCapNames)[1:-1] - - # --------- Regulator Information ---------- - regNames = dss.RegControls.AllNames() - hRegNames = [None] * len(regNames) - for i, n in enumerate(regNames): - hRegNames[i] = str(n) - hRegNames = str(hRegNames)[1:-1] - - dss.Text.Command("solve mode=snap") - # dss.Text.Command('export voltages') - v = np.array(dss.Circuit.AllBusMagPu()) + # TODO: Priti, this code isn't doing anything. Should we delete it? + # node_number = len(AllNodeNames) + # Vbase_allnode = [0] * node_number + # ii = 0 + # for node in AllNodeNames: + # dss.Circuit.SetActiveBus(node) + # Vbase_allnode[ii] = dss.Bus.kVBase() * 1000 + # ii = ii + 1 - volt_violation = np.any(v > upper_voltage_limit) or np.any(v < lower_voltage_limit) - logger.info("Initial condition voltage violation: %s", volt_violation) + # TODO: Priti: these variables are unused. Do we need them? + # hCapNames = [str(x)[1:-1] for x in dss.Capacitors.AllNames()] + # hRegNames = [str(x)[1:-1] for x in dss.RegControls.AllNames()] - Bus_Distance = [] - for node in AllNodeNames: - circuit.SetActiveBus(node) - Bus_Distance.append(dss.Bus.Distance()) - np.savetxt(output_dir / "node_distance.csv", Bus_Distance) + dss.Text.Command("solve mode=snap") + logger.info( + "Initial condition voltage violation: %s", + circuit_has_violations(lower_voltage_limit, upper_voltage_limit), + ) + bus_distances = list_bus_distances() + np.savetxt(output_dir / "node_distance.csv", bus_distances) ####### for thermal overload capapcity ####################################################### - dss.Text.Command("solve mode = snap") - overloads_filename = output_dir / "overloads.csv" - dss.Text.Command(f"export Overloads {overloads_filename}") - - # read this overload file and record the ' %Normal' for each line in this file - overload_df = pd.read_csv(overloads_filename) - # len(overload_df) - if len(overload_df) == 0: - logger.info("No thermal violation initially") - else: + overloads = get_loadings_with_violations(thermal_loading_limit) + if overloads: logger.info("Thermal violation exists initially") - logger.info("Overloads: %s", overload_df.head()) - - elements = [[], []] - amps = [[], []] - new_threshold = [[], []] - for j in range(len(extra_percentages_for_existing_overloads)): - for i in range(len(overload_df)): - overload_df["Element"] = overload_df["Element"].str.strip() - element = overload_df["Element"][i] - amp = overload_df[" %Normal"][i] - element_new_limit = amp + extra_percentages_for_existing_overloads[j] - elements[j].append(str(element)) - amps[j].append(amp) - new_threshold[j].append(element_new_limit) - - logger.debug("new_threshold=%s", new_threshold) - - ############################################################################################## - # get the load data - [Load, totalkW] = get_loads(dss, circuit, 0, "") - - # calculate the hosting capacity based on voltage constraints - v_output_df = [] - for j in range(len(v_limit)): - lmt1 = v_limit[j][0] - lmt2 = v_limit[j][1] - v_lst = [] - v_output_list = [] - v_threshold = [] - v_allow_limit = [] - v_names = [] - v_bus_name = [] - v_default_load = [] - v_maxv = [] - v_minv = [] - for i in range(len(Load)): - logger.info("Run voltage check on load index=%s", i) - v_node_cap = node_V_capacity_check( - master_file, Load[i], lmt1, lmt2, kw_step_voltage_violation - ) # v_node_cap is a list - v_allowable_load = v_node_cap[0] - kw_step_voltage_violation - v_threshold.append(v_node_cap[0]) - v_allow_limit.append(v_allowable_load) - v_names.append(Load[i]["name"]) - v_bus_name.append(Load[i]["bus1"]) - v_default_load.append(Load[i]["kW"]) - v_maxv.append(v_node_cap[1]) - v_minv.append(v_node_cap[2]) - - v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] - v_output_list = list(map(list, zip(*v_lst))) - # v_output_df[0] = pd.DataFrame(v_output_list, columns = ['Load' , 'Bus', 'Initial_kW', 'Volt_Violation', 'Maximum_kW', 'Max_voltage', 'Min_voltage']) - v_output_df.append( - pd.DataFrame( - v_output_list, - columns=[ - "Load", - "Bus", - "Initial_kW", - "Volt_Violation", - "Maximum_kW", - "Max_voltage", - "Min_voltage", - ], - ) - ) + logger.info("Overloads: %s", overloads) + else: + logger.info("No thermal violation initially") - v_output_df[j].to_csv(output_dir / f"Hosting_capacity_test_{lmt1}.csv") + elements_with_extra_threshold = { + k: v + extra_percentage_for_existing_overloads for k, v in overloads.items() + } + + logger.debug("new_threshold=%s", elements_with_extra_threshold) + loads = get_loads() + + v_output_df = calculate_voltage_hosting_capacity( + loads, + master_file, + lower_voltage_limit, + upper_voltage_limit, + kw_step_voltage_violation, + voltage_tolerance, + num_cpus, + ) + v_output_df.to_csv( + output_dir + / f"Hosting_capacity_voltage_test_{lower_voltage_limit}_{upper_voltage_limit}.csv" + ) # calculate the hosting capacity based on thermal ratings constraint - th_threshold = [[], []] - th_allow_limit = [[], []] - th_names = [[], []] - th_bus_name = [[], []] - th_default_load = [[], []] - - th_output_df = [] - for i in range(len(thermal_loading_limits)): - logger.info("Run thermal check on index=%s limit=%s", i, thermal_loading_limits[i]) - th_lst = [] - for j in range(len(Load)): - if j == 0: - continue - th_node_cap = thermal_overload_check( - Load[j], thermal_loading_limits[i], i, kw_step_thermal_violation, output_dir - ) # th_node_cap is a list - th_allowable_load = ( - th_node_cap[0] - kw_step_thermal_violation - ) # 5 # th_node_cap[0] is a float # reduce the value of add_kW - th_threshold[i].append(th_node_cap[0]) - th_allow_limit[i].append(th_allowable_load) - th_names[i].append(Load[j]["name"]) - th_bus_name[i].append(Load[j]["bus1"]) - th_default_load[i].append(Load[j]["kW"]) - - th_lst = [ - th_names[i], - th_bus_name[i], - th_default_load[i], - th_threshold[i], - th_allow_limit[i], - ] - th_output_list = list(map(list, zip(*th_lst))) - th_output_df.append( - pd.DataFrame( - th_output_list, - columns=["Load", "Bus", "Initial_kW", "Thermal_Violation", "Maximum_kW"], - ) - ) - - th_output_df[i].to_csv( - output_dir / f"Thermal_capacity_test_{thermal_loading_limits[i]}.csv" - ) + th_output_df = calculate_thermal_hosting_capacity( + loads, + master_file, + thermal_loading_limit, + kw_step_thermal_violation, + thermal_tolerance, + elements_with_extra_threshold, + num_cpus, + ) + th_output_df.to_csv(output_dir / f"Hosting_capacity_thermal_test_{thermal_loading_limit}.csv") #################### for plotting results #################################################### load_bus = pd.DataFrame() - load_bus["Load"] = th_output_df[0]["Load"] # - load_bus["Bus"] = th_output_df[0]["Bus"] + load_bus["Load"] = th_output_df["Load"] # + load_bus["Bus"] = th_output_df["Bus"] node_distance = pd.DataFrame() - node_distance["Node"] = AllNodeNames - node_distance["Distance"] = Bus_Distance + node_distance["Node"] = dss.Circuit.AllNodeNames() + node_distance["Distance"] = bus_distances dist_file = calc_dist(load_bus, node_distance) - dist_file["Initial_MW"] = th_output_df[0]["Initial_kW"] / 1000 - dist_file["Volt_Violation_0.95"] = v_output_df[0]["Volt_Violation"] / 1000 - dist_file["Volt_Violation_0.975"] = v_output_df[1]["Volt_Violation"] / 1000 - dist_file["Volt_Violation_0.985"] = v_output_df[2]["Volt_Violation"] / 1000 - - dist_file["Thermal_Violation_100"] = th_output_df[0]["Thermal_Violation"] / 1000 - # dist_file["Thermal_Violation_120"] = th_output_df[1]["Thermal_Violation"]/1000 + dist_file["Initial_MW"] = th_output_df["Initial_kW"] / 1000 + dist_file[f"Volt_Violation_{lower_voltage_limit}"] = v_output_df["Volt_Violation"] / 1000 + dist_file[f"Thermal_Violation_{thermal_loading_limit}"] = ( + th_output_df["Thermal_Violation"] / 1000 + ) plot_df = dist_file.sort_values(by=["Distance"]) @@ -250,57 +187,68 @@ def run( plot_capacity_V( plot_df, "Initial_MW", - "Volt_Violation_0.95", - "Volt_Violation_0.975", - "Volt_Violation_0.985", + f"Volt_Violation_{lower_voltage_limit}", output_dir, ) # plot thermal violation - plot_capacity_thermal_1(plot_df, "Initial_MW", "Thermal_Violation_100", output_dir, 100) - # plot_capacity_thermal_2(plot_df, 'Initial_MW', 'Thermal_Violation_100', 'Thermal_Violation_120', output_dir) + # TODO: Priti, is the last parameter correct? + plot_capacity_thermal_1( + plot_df, + "Initial_MW", + f"Thermal_Violation_{thermal_loading_limit}", + output_dir, + thermal_loading_limit, + ) ### Assuming the hosting capacity is limited by thermal loading ############################## - ## Difference of initial load and maximum hosting capacity (assuming always thermal limit occurs first) ## + ## Difference of initial load and maximum hosting capacity (assuming always thermal limit occurs first) + diff = th_output_df["Thermal_Violation"] - th_output_df["Initial_kW"] + new_df = pd.DataFrame() - for i in range(len(thermal_loading_limits)): - diff = th_output_df[i]["Thermal_Violation"] - th_output_df[i]["Initial_kW"] - new_df = pd.DataFrame() + new_df["Load"] = th_output_df["Load"] + new_df["Bus"] = th_output_df["Bus"] + new_df["Initial_kW"] = th_output_df["Initial_kW"] + new_df["Hosting_capacity(kW)"] = diff # additional load it can support - new_df["Load"] = th_output_df[i]["Load"] - new_df["Bus"] = th_output_df[i]["Bus"] - new_df["Initial_kW"] = th_output_df[i]["Initial_kW"] - new_df["Hosting_capacity(kW)"] = diff # additional load it can support + new_df.to_csv( + output_dir / f"Additional_HostingCapacity_{thermal_loading_limit}.csv", + index=False, + ) - new_df.to_csv( - output_dir / "Additional_HostingCapacity_" + str(thermal_loading_limits[i]) + ".csv", - index=False, - ) + # Find number of ev chargers for each node. + chargers_3_2_1 = levels_of_charger(th_output_df) + chargers_3_2_1.to_csv(output_dir / f"Loadwithlevel3_2_1_{thermal_loading_limit}.csv") - ## find number of ev chargers for each node ## - chargers_3_2_1 = levels_of_charger(th_output_df[i]) - chargers_3_2_1.to_csv( - output_dir / "Loadwithlevel3_2_1_" + str(thermal_loading_limits[i]) + ".csv" - ) +def circuit_has_violations(lower_voltage_limit, upper_voltage_limit) -> bool: + """Returns True if the current circuit has voltage violations.""" + v = np.array(dss.Circuit.AllBusMagPu()) + return np.any(v > upper_voltage_limit) or np.any(v < lower_voltage_limit) -# --------- Load Information ---------- -@track_timing(timer_stats_collector) -def get_loads(dss, circuit, loadshape_flag, loadshape_folder): - data = [] +def get_loadings_with_violations(threshold: float) -> dict[str, float]: + """Return a dict of elements with loading violations.""" + return { + x: y + for (x, y) in zip(dss.PDElements.AllNames(), dss.PDElements.AllPctNorm()) + if y >= threshold + } + + +def get_loads() -> list[dict]: + """Return a list of all Loads in the circuit.""" + loads = [] load_flag = dss.Loads.First() - total_load = 0 while load_flag: - load = dss.Loads datum = { - "name": load.Name(), - "kV": load.kV(), - "kW": load.kW(), - "PF": load.PF(), - "Delta_conn": load.IsDelta(), + "name": dss.Loads.Name(), + "kV": dss.Loads.kV(), + "kW": dss.Loads.kW(), + "PF": dss.Loads.PF(), + "Delta_conn": dss.Loads.IsDelta(), } cktElement = dss.CktElement @@ -320,113 +268,278 @@ def get_loads(dss, circuit, loadshape_flag, loadshape_folder): datum["voltageAng"] = cktElement.VoltagesMagAng()[1] datum["power"] = dss.CktElement.Powers()[0:2] - data.append(datum) + loads.append(datum) load_flag = dss.Loads.Next() - total_load += datum["kW"] - return [data, total_load] + return loads + + +def list_bus_distances() -> list[float]: + """Return a list of bus distances.""" + distances = [] + for node in dss.Circuit.AllNodeNames(): + dss.Circuit.SetActiveBus(node) + distances.append(dss.Bus.Distance()) + return distances + + +def calculate_voltage_hosting_capacity( + loads, + master_file, + lower_limit, + upper_limit, + kw_step_voltage_violation, + voltage_tolerance, + num_cpus, +) -> pd.DataFrame: + """Calculate hosting capacity based on voltage and store the result in a DataFrame.""" + v_lst = [] + v_output_list = [] + v_threshold = [] + v_allow_limit = [] + v_names = [] + v_bus_name = [] + v_default_load = [] + v_maxv = [] + v_minv = [] + with ProcessPoolExecutor(max_workers=num_cpus) as executor: + for result in executor.map( + node_V_capacity_check, + loads, + itertools.repeat(master_file), + itertools.repeat(lower_limit), + itertools.repeat(upper_limit), + itertools.repeat(kw_step_voltage_violation), + itertools.repeat(voltage_tolerance), + ): + load, cap_limit, vmax, vmin = result + logger.debug("Ran voltage check on load=%s", load["name"]) + v_allowable_load = cap_limit - kw_step_voltage_violation + v_threshold.append(cap_limit) + v_allow_limit.append(v_allowable_load) + v_names.append(load["name"]) + v_bus_name.append(load["bus1"]) + v_default_load.append(load["kW"]) + v_maxv.append(vmax) + v_minv.append(vmin) + + v_lst = [v_names, v_bus_name, v_default_load, v_threshold, v_allow_limit, v_maxv, v_minv] + v_output_list = list(map(list, zip(*v_lst))) + v_output_df = pd.DataFrame( + v_output_list, + columns=[ + "Load", + "Bus", + "Initial_kW", + "Volt_Violation", + "Maximum_kW", + "Max_voltage", + "Min_voltage", + ], + ) + return v_output_df -############################ for voltage violation capacity ###################################### -@track_timing(timer_stats_collector) -def node_V_capacity_check(master_file, which_load, low_lmt, high_lmt, kw_step_voltage_violation): - initial_kW = which_load["kW"] +def node_V_capacity_check( + load, master_file, lower_limit, upper_limit, kw_step_voltage_violation, tolerance +): + """Returns the lowest capacity at which a violation occurs.""" + compile_circuit(master_file) + initial_kW = load["kW"] + kw_step_voltage_violation + initial_value = initial_kW + new_kW = initial_kW logger.debug("initial_kW=%s", initial_kW) - tmp_kW = initial_kW - volt_violation = False - v = None - - while not volt_violation: - with Timer(timer_stats_collector, "check_voltage_violation_after_increasing_kw"): - # while the voltages are within limits keep on increasing - new_kW = tmp_kW + kw_step_voltage_violation - logger.info("Check voltage violation with kW=%s", new_kW) - dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) - dss.Text.Command("solve mode = snap") - v = np.array(dss.Circuit.AllBusMagPu()) - volt_violation = np.any(v > high_lmt) or np.any(v < low_lmt) - logger.debug("Voltage violation = %s", volt_violation) - if volt_violation: - vmax = np.max(v) - vmin = np.min(v) - cap_limit = new_kW - # TODO: no point in doing this, right? - # dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) - compile_circuit(master_file) - else: - tmp_kW = new_kW + bisector = ViolationBisector( + lower_bound=initial_kW, initial_value=initial_value, tolerance=tolerance + ) + done = False + vmax = 0.0 + vmin = 0.0 + cap_limit = 0.0 + + while not done: + logger.debug("Check voltage violation with kW=%s", new_kW) + voltages = set_load_and_solve(load["name"], new_kW) + last_result_violation = circuit_has_violations(lower_limit, upper_limit) + logger.debug( + "Checked voltage violation load=%s kW=%s last_result_violation=%s", + load["name"], + new_kW, + last_result_violation, + ) + new_kW, done = bisector.get_next_value(last_result_violation=last_result_violation) + + cap_limit = bisector.get_lowest_violation() + set_load_and_solve(load["name"], cap_limit) + voltages = np.array(dss.Circuit.AllBusMagPu()) + vmax = np.max(voltages) + vmin = np.min(voltages) + return load, cap_limit, vmax, vmin + - return [cap_limit, vmax, vmin] +def set_load_and_solve(load_name, kw): + dss.Text.Command(f"edit Load.{load_name} kW={kw}") + dss.Text.Command("solve mode = snap") + + +def calculate_thermal_hosting_capacity( + loads, + master_file, + loading_limit, + kw_step_thermal_violation, + thermal_tolerance, + elements_with_extra_threshold, + num_cpus, +): + """Calculate hosting capacity based on voltage and store the result in a DataFrame.""" + th_threshold = [] + th_allow_limit = [] + th_names = [] + th_bus_name = [] + th_default_load = [] + th_lst = [] + with ProcessPoolExecutor(max_workers=num_cpus) as executor: + for result in executor.map( + thermal_overload_check, + loads, + itertools.repeat(master_file), + itertools.repeat(loading_limit), + itertools.repeat(kw_step_thermal_violation), + itertools.repeat(thermal_tolerance), + itertools.repeat(elements_with_extra_threshold), + ): + load, th_node_cap = result + logger.debug("Ran thermal check on limit=%s load=%s", loading_limit, load["name"]) + # reduce the value of add_kW + th_allowable_load = th_node_cap - kw_step_thermal_violation + th_threshold.append(th_node_cap) + th_allow_limit.append(th_allowable_load) + th_names.append(load["name"]) + th_bus_name.append(load["bus1"]) + th_default_load.append(load["kW"]) + + th_lst = [ + th_names, + th_bus_name, + th_default_load, + th_threshold, + th_allow_limit, + ] + th_output_list = list(map(list, zip(*th_lst))) + th_output_df = pd.DataFrame( + th_output_list, + columns=["Load", "Bus", "Initial_kW", "Thermal_Violation", "Maximum_kW"], + ) + return th_output_df -@track_timing(timer_stats_collector) def thermal_overload_check( - which_load, th_limit, case, kw_step_thermal_violation, output_dir: Path + load, + master_file, + th_limit, + kw_step_thermal_violation, + tolerance, + elements_with_extra_threshold, ): - initial_kW = which_load["kW"] + """Returns the lowest capacity at which a violation occurs.""" + compile_circuit(master_file) + initial_kW = load["kW"] + kw_step_thermal_violation + new_kW = initial_kW logger.debug("initial_kW=%s", initial_kW) - tmp_kW = initial_kW - thermal_violation = False - - while not thermal_violation: - with Timer(timer_stats_collector, "check_thermal_violation_after_increasing_kw"): - # while the elements loadings are within limits keep on increasing - new_kW = tmp_kW + kw_step_thermal_violation - logger.info("Check thermal violation with kW=%s", new_kW) - dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(new_kW)) - dss.Text.Command("solve mode = snap") - overloads_filename = output_dir / "overloads.csv" - dss.Text.Command("export Overloads {overloads_filename}") - report = pd.read_csv(overloads_filename) - report["Element"] = report["Element"].str.strip() - - if len(report) == 0: # if no any overload element - thermal_violation = False - - elif report["Element"].isin(elements[case]).any(): - for i in range(len(report)): - if report["Element"][i] in elements[case]: - indx_ = elements[case].index(report["Element"][i]) # find the index of - if report[" %Normal"][i] >= new_threshold[case][indx_]: - thermal_violation = True # just exit here (get out of for loop) - break - else: - thermal_violation = False - else: - # check the percentage normal if greater than specified % then only violation - if report[" %Normal"][i] >= th_limit: - thermal_violation = True - break - else: - thermal_violation = False + bisector = ViolationBisector( + lower_bound=initial_kW, initial_value=initial_kW, tolerance=tolerance + ) + done = False + load_name = load["name"] + + while not done: + # while the elements loadings are within limits keep on increasing + cmd = f"edit Load.{load_name} kW={new_kW}" + dss.Text.Command(cmd) + dss.Text.Command("solve mode = snap") + loadings = get_loadings_with_violations(th_limit) + + last_result_violation = False + if not loadings: + last_result_violation = False + elif set(loadings).intersection(elements_with_extra_threshold): + for name, pct_loading in loadings.items(): + if name in elements_with_extra_threshold: + if pct_loading >= elements_with_extra_threshold[name]: + last_result_violation = True + break + else: + last_result_violation = True + break + else: + last_result_violation = True + logger.debug( + "Checked thermal violation load=%s kW=%s last_result_violation=%s", + load_name, + new_kW, + last_result_violation, + ) + new_kW, done = bisector.get_next_value(last_result_violation=last_result_violation) + + return load, bisector.get_lowest_violation() + + +class ViolationBisector: + """Helper class to find the lowest value that will cause a violation.""" + + def __init__(self, lower_bound, initial_value=None, tolerance=10): + self._lower_bound = lower_bound + self._initial_value = initial_value or lower_bound + self._cur_value = self._initial_value + self._last_fail_value = None + self._found_first_violation = False + self._lowest_violation = sys.maxsize + self._highest_pass = -1.0 + self._tolerance = tolerance + + def get_next_value(self, last_result_violation: bool): + done = False + if last_result_violation: + if not self._found_first_violation: + self._found_first_violation = True + assert self._cur_value < self._lowest_violation + self._lowest_violation = self._cur_value + if self._lowest_violation <= self._lower_bound: + done = True else: - # check the percentage normal if greater than specified % then only violation - for i in range(len(report)): - if report[" %Normal"][i] >= th_limit: - thermal_violation = True - break - else: - thermal_violation = False - - logger.debug("thermal_violation=%s", thermal_violation) - if thermal_violation: - logger.debug("thermal_violation=%s", thermal_violation) - dss.Text.Command(f"export Capacity {output_dir / 'capacity.csv'}") - dss.Text.Command(f"export Currents {output_dir / 'currents.csv'}") - cap_limit = new_kW - # TODO: no point in doing this, right? - # dss.Text.Command("edit Load." + str(which_load["name"]) + " kW=" + str(initial_kW)) - compile_circuit(master_file) + low = self._highest_pass or self._lower_bound + self._cur_value = low + ((self._lowest_violation - low) // 2) + else: + assert self._highest_pass is None or self._cur_value > self._highest_pass + self._highest_pass = self._cur_value + if self._found_first_violation: + self._cur_value += (self._lowest_violation - self._cur_value) // 2 else: - tmp_kW = new_kW + self._cur_value *= 2 + + if not done and (self._lowest_violation is not None and self._highest_pass is not None): + assert ( + self._lowest_violation > self._highest_pass + ), f"{self._highest_pass=} {self._lowest_violation=}" + done = self._lowest_violation - self._highest_pass <= self._tolerance + + logger.debug( + "end last_result=%s cur_value=%s lowest_fail=%s highest_pass=%s", + last_result_violation, + self._cur_value, + self._lowest_violation, + self._highest_pass, + ) + return self._cur_value, done - return [cap_limit] + def get_lowest_violation(self): + return self._lowest_violation -@track_timing(timer_stats_collector) def compile_circuit(master_file: Path): - logger.info("Compile circuit from %s", master_file) + """Compiles a circuit and ensures that execution remains in the current directory.""" + logger.debug("Compile circuit from %s", master_file) orig = os.getcwd() try: dss.Text.Command(f"Compile {master_file}") @@ -434,7 +547,6 @@ def compile_circuit(master_file: Path): os.chdir(orig) -@track_timing(timer_stats_collector) def export_circuit_element_properties(output_dir: Path): export_dir = output_dir / "circuit_elements" export_dir.mkdir() @@ -465,5 +577,5 @@ def export_circuit_element_properties(output_dir: Path): else: logger.info("There are no elements of type=%s to export", class_name) - all_node_names = circuit.YNodeOrder() + all_node_names = dss.Circuit.AllNodeNames() pd.DataFrame(all_node_names).to_csv(output_dir / "Allnodenames.csv") diff --git a/disco/ev/plot_hosting_capacity.py b/disco/ev/plot_hosting_capacity.py index 78c93d37..7dbd1186 100644 --- a/disco/ev/plot_hosting_capacity.py +++ b/disco/ev/plot_hosting_capacity.py @@ -14,30 +14,32 @@ logger = logging.getLogger(__name__) -def plot_capacity_V(file, caseA, caseB, caseC, caseD, output_dir: Path): - fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) # figsize=(18,12) - vplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") - vplot2 = ax1.plot( +def plot_capacity_V(file, caseA, caseB, output_dir: Path, caseC=None, caseD=None): + _, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) # figsize=(18,12) + ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") + ax1.plot( range(len(file)), file[caseB], c="r", marker="+", - label="Volt violation MW [range:(0.95, 1.05)]", - ) - vplot3 = ax1.plot( - range(len(file)), - file[caseC], - c="g", - marker="*", - label="Volt violation MW [range:(0.975, 1.05)]", - ) - vplot4 = ax1.plot( - range(len(file)), - file[caseD], - c="b", - marker=".", - label="Volt violation MW [range:(0.985, 1.05)]", + label=f"Volt violation MW [range:({caseB})]", # TODO: not a range ) + if caseC is not None: + ax1.plot( + range(len(file)), + file[caseC], + c="g", + marker="*", + label=f"Volt violation MW [range:({caseC})]", # TODO: not a range + ) + if caseD is not None: + ax1.plot( + range(len(file)), + file[caseD], + c="b", + marker=".", + label=f"Volt violation MW [range:({caseD})]", # TODO: not a range + ) ax1.margins(x=0) ax1.tick_params(axis="both", which="major", labelsize=33) ax1.tick_params(axis="both", which="minor", labelsize=33) @@ -85,11 +87,9 @@ def plot_capacity_thermal_2(file, caseA, caseB, caseC, output_dir: Path): def plot_capacity_thermal_1(file, caseA, caseB, output_dir: Path, extra): - fig, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) - thplot1 = ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") - thplot2 = ax1.plot( - range(len(file)), file[caseB], c="r", marker="+", label="Thermal violation MW" - ) + _, ax1 = plt.subplots(nrows=1, figsize=(14, 8)) + ax1.plot(range(len(file)), file[caseA], c="k", marker="o", label="Initial load") + ax1.plot( range(len(file)), file[caseB], c="r", marker="+", label="Thermal violation MW") ax1.margins(x=0) ax1.tick_params(axis="both", which="major", labelsize=33) ax1.tick_params(axis="both", which="minor", labelsize=33) @@ -102,7 +102,7 @@ def plot_capacity_thermal_1(file, caseA, caseB, output_dir: Path, extra): ) # ax1.set_xlabel('Load indices, according to distance from SS \n(load #' + str(len(file)) + ' is the furthest from SS)', fontsize=40) # fig.tight_layout() plt.savefig( - output_dir / "Cap_by_thermal_limit_" + str(extra) + ".png", bbox_inches="tight", dpi=150 + output_dir / f"Cap_by_thermal_limit_{extra}.png", bbox_inches="tight", dpi=150 )