From ee6ebc22e37b4dd0e038c188f45a42733098dcc5 Mon Sep 17 00:00:00 2001 From: Daniel Thom Date: Thu, 15 Jun 2023 11:56:03 -0600 Subject: [PATCH] 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 )