diff --git a/pypsse/cli/explore.py b/pypsse/cli/explore.py index fdb5984..b539a7c 100644 --- a/pypsse/cli/explore.py +++ b/pypsse/cli/explore.py @@ -1,5 +1,5 @@ from pathlib import Path - +import os from loguru import logger import pandas as pd import click @@ -107,6 +107,10 @@ def explore(project_path, simulations_file, export_file_path, load_filter, load, } results = x.sim.read_subsystems(quantities, buses) + # print(results.keys()) + # print(results["LOAD_P"]) + # quit() + had_comp_models = False if "Loads_FmA" in results: had_comp_models = True @@ -129,7 +133,7 @@ def explore(project_path, simulations_file, export_file_path, load_filter, load, is_comp_load[bus] = is_comp load_dict[bus].append(ld_id) - key = f"{ld_id} _{bus}" if len(ld_id) == 1 else f"{ld_id}_{bus}" + key = f"{bus}_{ld_id}" key2 = f"{bus}_{ld_id}".replace(" ", "") load_p = max( results["Loads_MVA"][key].real + results["Loads_IL"][key].real + results["Loads_YL"][key].real, @@ -149,7 +153,7 @@ def explore(project_path, simulations_file, export_file_path, load_filter, load, if bus not in generator_dict: generator_dict[bus] = [] bus_gen[bus] = 0 - key = f"{gen_id} _{bus}" if len(gen_id) == 1 else f"{gen_id}_{bus}" + key = f"{bus}_{gen_id}" generator_dict[bus].append(gen_id) bus_gen[bus] += results["Machines_MVA"][key] @@ -168,7 +172,7 @@ def explore(project_path, simulations_file, export_file_path, load_filter, load, results["is load comp"].append(is_comp_load[bus] if bus in is_comp_load else False) results["total P load [MW]"].append(bus_load_real[bus] if bus in bus_load_real else 0) results["total Q load [MVar]"].append(bus_load_imag[bus] if bus in bus_load_imag else 0) - results["has generation"].append(True if bus in generator_dict else False) + results["has generation"].append(True if (bus in generator_dict and bus_gen[bus] > 0)else False) results["total generation [MVA]"].append(bus_gen[bus] if bus in bus_gen else 0) @@ -200,7 +204,7 @@ def explore(project_path, simulations_file, export_file_path, load_filter, load, results=results[(results["total P load [MW]"] >= load_lower) & (results["total P load [MW]"] <= load_upper)] results=results[(results["total generation [MVA]"] >= gen_lower) & (results["total generation [MVA]"] <= gen_upper)] - print(results) + # print(results) results.to_csv(export_file_path) logger.info(f"Results exported to {export_file_path.absolute()}") diff --git a/pypsse/cli/profiles.py b/pypsse/cli/profiles.py index 826947f..f1be9a8 100644 --- a/pypsse/cli/profiles.py +++ b/pypsse/cli/profiles.py @@ -1,17 +1,17 @@ """ -CLI to run a PyDSS project -""" + CLI to run a PyDSS project + """ from pathlib import Path - + from loguru import logger import click import toml - + from pypsse.models import SimulationSettings from pypsse.common import SIMULATION_SETTINGS_FILENAME from pypsse.profile_manager_interface import ProfileManagerInterface - + @click.argument( "project-path", ) @@ -41,4 +41,3 @@ def get_profiles(project_path, simulations_file=None): profile_interface = ProfileManagerInterface.from_setting_files(file_path) profile_interface.get_profiles() - diff --git a/pypsse/cli/run.py b/pypsse/cli/run.py index 1a67dcf..04df06d 100644 --- a/pypsse/cli/run.py +++ b/pypsse/cli/run.py @@ -7,6 +7,7 @@ from loguru import logger import click import toml +import os from pypsse.models import SimulationSettings from pypsse.common import SIMULATION_SETTINGS_FILENAME @@ -34,11 +35,15 @@ def run(project_path, simulations_file=None): simulation_settiings = toml.load(file_path) simulation_settiings = SimulationSettings(**simulation_settiings) - + if simulation_settiings.log.clear_old_log_file: + log_path = Path(project_path) / "Logs" / "pypsse.log" + if os.path.exists(log_path): + os.remove(log_path) logger.level(simulation_settiings.log.logging_level.value) if simulation_settiings.log.log_to_external_file: log_path = Path(project_path) / "Logs" / "pypsse.log" logger.add(log_path) + x = Simulator.from_setting_files(file_path) diff --git a/pypsse/contingencies.py b/pypsse/contingencies.py index f0faf8f..7b711a0 100644 --- a/pypsse/contingencies.py +++ b/pypsse/contingencies.py @@ -46,7 +46,14 @@ def __init__(self, psse, settings, contingency_type): self.psse = psse self.enabled = False self.tripped = False - + logger.debug( + f"contingency_type : {contingency_type}" + ) + logger.debug( + f"settings : {settings}" + ) + # os.system("PAUSE") + def update(self, t: float): """updates a fault event @@ -55,30 +62,27 @@ def update(self, t: float): """ self.t = t if hasattr(self.settings, "duration"): - if ( - self.settings.time + self.settings.duration - > t - >= self.settings.time - and not self.enabled - ): + if (self.settings.time + self.settings.duration > t >= self.settings.time and not self.enabled): self.enabled = True self.enable_fault() - if ( - t >= self.settings.time + self.settings.duration - and self.enabled - ): + if (t >= self.settings.time + self.settings.duration and self.enabled): self.enabled = False self.disable_fault() - elif ( - not hasattr(self.settings, "duration") - and t >= self.settings.time - and not self.tripped - ): + elif (not hasattr(self.settings, "duration") and t >= self.settings.time and not self.tripped): self.enable_fault() self.tripped = True def enable_fault(self): """enables a fault event""" + data_check = getattr(self.psse, self.fault_method) + logger.debug( + f"data_check : {data_check}" + ) + logger.debug( + f"data_check : {self.fault_method}" + ) + # os.system("PAUSE") + err = getattr(self.psse, self.fault_method)(**self.fault_settings) if err: logger.warning( diff --git a/pypsse/data_writers/hdf5.py b/pypsse/data_writers/hdf5.py index d6e3b3e..5811781 100644 --- a/pypsse/data_writers/hdf5.py +++ b/pypsse/data_writers/hdf5.py @@ -119,7 +119,7 @@ def write( if self.step >= len(self.Timestamp): self.Timestamp.resize((len(self.Timestamp) + 1,)) self.convergence.resize((len(self.convergence) + 1,)) - self.Timestamp[self.step - 1] = np.string_(currenttime.strftime("%Y-%m-%d %H:%M:%S.%f")) + self.Timestamp[self.step - 1] = np.bytes_(currenttime.strftime("%Y-%m-%d %H:%M:%S.%f")) self.convergence[self.step - 1] = convergence # Add object status data to a DataFrame self.store.flush() diff --git a/pypsse/enumerations.py b/pypsse/enumerations.py index 4a6c37c..6b42626 100644 --- a/pypsse/enumerations.py +++ b/pypsse/enumerations.py @@ -28,10 +28,16 @@ class SimulationModes(str, Enum): STATIC = "Steady-state" DYNAMIC = "Dynamic" +class GenerationLevel(str, Enum): + "Valid generation level setting modes" + TRANSMISSION = "transmission" + DISTRIBUTION = "distribution" class HelicsCoreTypes(str, Enum): "HELICS core types" ZMQ = "zmq" + TCP_SS="tcp_ss" + TCP="tcp" class WritableModelTypes(str, Enum): @@ -40,7 +46,9 @@ class WritableModelTypes(str, Enum): PLANT = "Plant" MACHINE = "Machine" GENERATOR = "Induction_machine" - + LOAD_STATUS = "Load_status" + LINE_STATUS = "Line_status" + MACHINE_STATUS = "Machine_status" class ModelTypes(str, Enum): "Supported asset tpyes in PyPSSE" diff --git a/pypsse/helics_interface.py b/pypsse/helics_interface.py index f668bc8..29bd9e8 100644 --- a/pypsse/helics_interface.py +++ b/pypsse/helics_interface.py @@ -1,5 +1,6 @@ import ast - +import os +import time import helics as h import pandas as pd from loguru import logger @@ -52,6 +53,21 @@ def __init__( self.subsystem_info = [] self.publications = {} self.subscriptions = {} + self.load_fault = False + ################################################################# + # FX: add these hardcored load value to better matching + self.load_power = { + tid: { + "transmission_loads_at_fault": at_fault, + "transmission_loads_clear_fault": clear_fault + } + for tid, at_fault, clear_fault in zip( + settings.simulation.transmission_ids, + settings.simulation.transmission_loads_at_fault, + settings.simulation.transmission_loads_clear_fault + ) + } + ################################################################# def enter_execution_mode(self): """Enables federate to enter execution mode""" @@ -66,11 +82,16 @@ def enter_execution_mode(self): def create_federate(self): """Creates a HELICS co-simulation federate""" self.fedinfo = h.helicsCreateFederateInfo() + # logger.debug(f"self.fedinfo: {self.fedinfo}") h.helicsFederateInfoSetCoreName(self.fedinfo, self.settings.helics.federate_name) h.helicsFederateInfoSetCoreTypeFromString(self.fedinfo, self.settings.helics.core_type.value) h.helicsFederateInfoSetCoreInitString(self.fedinfo, "--federates=1") h.helicsFederateInfoSetBroker(self.fedinfo, str(self.settings.helics.broker_ip)) h.helicsFederateInfoSetBrokerPort(self.fedinfo, self.settings.helics.broker_port) + IP = self.settings.helics.broker_ip + Port = self.settings.helics.broker_port + logger.info("Connecting to broker @ {}".format(f"{IP}:{Port}" if Port else IP)) + logger.info(f"Connecting to broker @ {self.settings.helics}") if self.settings.helics.iterative_mode: h.helicsFederateInfoSetTimeProperty( @@ -100,6 +121,7 @@ def register_publications(self, bus_subsystems: dict): self.publications = {} self.pub_struc = [] for publication_dict in self.settings.helics.publications: + # logger.debug(f"publication_dict: {publication_dict}") bus_subsystem_ids = publication_dict.bus_subsystems if not set(bus_subsystem_ids).issubset(self.bus_subsystems): msg = f"One or more invalid bus subsystem ID pass in {bus_subsystem_ids}." @@ -115,6 +137,8 @@ def register_publications(self, bus_subsystems: dict): managed_properties = [ptpy.value for ptpy in getattr(self.export_settings, elm_class.lower())] properties = [p.value for p in publication_dict.asset_properties] + # logger.debug(f"managed_properties: {managed_properties}") + # logger.debug(f"properties: {properties}") if not set(properties).issubset(managed_properties): msg = f"One or more publication property defined for class '{elm_class}' is invalid." @@ -133,6 +157,7 @@ def register_publications(self, bus_subsystems: dict): self.pub_struc.append([{elm_class: properties}, bus_cluster]) temp_res = self.sim.read_subsystems({elm_class: properties}, bus_cluster) temp_res = self.get_restructured_results(temp_res) + for c_name, elm_info in temp_res.items(): for name, v_info in elm_info.items(): for p_name, val in v_info.items(): @@ -235,11 +260,12 @@ def register_subscriptions(self): self.psse_dict[row["bus"]][row["element_type"]][element_id] = {} if isinstance(row["element_property"], str): if row["element_property"] not in self.psse_dict[row["bus"]][row["element_type"]][element_id]: - self.psse_dict[row["bus"]][row["element_type"]][element_id][row["element_property"]] = 0 + # FX: modify this so that multiple feeders can be connected to a Transmission Bus + self.psse_dict[row["bus"]][row["element_type"]][element_id][row["element_property"]] = [] elif isinstance(row["element_property"], list): for r in row["element_property"]: if r not in self.psse_dict[row["bus"]][row["element_type"]][element_id]: - self.psse_dict[row["bus"]][row["element_type"]][element_id][r] = 0 + self.psse_dict[row["bus"]][row["element_type"]][element_id][r] = [] def request_time(self, _) -> (bool, float): """Enables time increment of the federate ina co-simulation." @@ -253,8 +279,10 @@ def request_time(self, _) -> (bool, float): bool: flag for iteration requrirement (rerun same time step) float: current helics time in seconds """ - + r_seconds = self.sim.get_total_seconds() # - self._dss_solver.GetStepResolutionSeconds() + logger.info(f"Time requested: {r_seconds}") + if self.sim.get_time() not in self.all_sub_results: self.all_sub_results[self.sim.get_time()] = {} self.all_pub_results[self.sim.get_time()] = {} @@ -372,20 +400,24 @@ def subscribe(self) -> dict: """ "Subscribes results each iteration and updates PSSE objects accordingly" + + # logger.debug(f"self.psse_dict is {self.psse_dict}") for sub_tag, sub_data in self.subscriptions.items(): if isinstance(sub_data["property"], str): sub_data["value"] = h.helicsInputGetDouble(sub_data["subscription"]) + # logger.debug(f"sub_data is {sub_data}") self.psse_dict[sub_data["bus"]][sub_data["element_type"]][sub_data["element_id"]][ sub_data["property"] - ] = (sub_data["value"], sub_data["scaler"]) + ].append((sub_data["value"], sub_data["scaler"])) elif isinstance(sub_data["property"], list): sub_data["value"] = h.helicsInputGetVector(sub_data["subscription"]) + # logger.debug(f"sub_data is {sub_data}") if isinstance(sub_data["value"], list) and len(sub_data["value"]) == len(sub_data["property"]): for i, p in enumerate(sub_data["property"]): - self.psse_dict[sub_data["bus"]][sub_data["element_type"]][sub_data["element_id"]][p] = ( + self.psse_dict[sub_data["bus"]][sub_data["element_type"]][sub_data["element_id"]][p].append(( sub_data["value"][i], sub_data["scaler"][i], - ) + )) logger.debug("Data received {} for tag {}".format(sub_data["value"], sub_tag)) if self.settings.helics.iterative_mode: @@ -399,31 +431,76 @@ def subscribe(self) -> dict: for i, v_dict in t_info.items(): values = {} j = 0 - for p, v_raw in v_dict.items(): - if isinstance(v_raw, tuple): - v, scale = v_raw - all_values[f"{t}.{b}.{i}.{p}"] = v + for p, v_raws in v_dict.items(): + if isinstance(v_raws, list): if isinstance(p, str): ppty = f"realar{PROFILE_VALIDATION[t].index(p) + 1}" - values[ppty] = v * scale + values[ppty] = 0 + for v_raw in v_raws: + v, scale = v_raw + if isinstance(scale, (float, int)): + values[ppty] += v * scale + else: + values[ppty] += v + all_values[f"{t}.{b}.{i}.{p}"] = values[ppty] elif isinstance(p, list): for _, ppt in enumerate(p): ppty = f"realar{PROFILE_VALIDATION[t].index(ppt) + 1}" - values[ppty] = v * scale + values[ppty] = 0 + for v_raw in v_raws: + v, scale = v_raw + if isinstance(scale, (float, int)): + values[ppty] += v * scale + else: + values[ppty] += v + all_values[f"{t}.{b}.{i}.{p}"] = values[ppty] j += 1 - is_empty = [0 if not vx else 1 for vx in values.values()] + logger.debug(f"{t}.{b}.{i} = {values}") + logger.debug(f"current HELICE time: {self.c_seconds}") + ###################################################### + ## FX: add this for better transit matching + if round(self.c_seconds, 3) == 0.1 and self.settings.simulation.transmission_loads_markup: + logger.debug("the moment of fault") + logger.debug(f"old {t}.{b}.{i} = {values}") + logger.debug(self.load_power) + load_data = self.load_power[b]['transmission_loads_at_fault'] + values['realar1'] = load_data[0] + values['realar2'] = load_data[1] + # os.system("PAUSE") + if round(self.c_seconds, 3) == 0.175 and self.settings.simulation.transmission_loads_markup: + logger.debug("the moment of clearing fault") + logger.debug(f"old {t}.{b}.{i} = {values}") + load_data = self.load_power[b]['transmission_loads_clear_fault'] + values['realar1'] = load_data[0] + values['realar2'] = load_data[1] + # os.system("PAUSE") + ###################################################### if ( sum(is_empty) != 0 and sum(values.values()) < VALUE_UPDATE_BOUND and sum(values.values()) > -VALUE_UPDATE_BOUND + and self.c_seconds > 0.02 ): self.sim.update_object(t, b, i, values) logger.debug(f"{t}.{b}.{i} = {values}") else: logger.debug("write failed") + + ###################################################### + ## FX: clear the result list + for sub_tag, sub_data in self.subscriptions.items(): + if isinstance(sub_data["property"], str): + self.psse_dict[sub_data["bus"]][sub_data["element_type"]][sub_data["element_id"]][sub_data["property"]] = [] + elif isinstance(sub_data["property"], list): + if isinstance(sub_data["value"], list) and len(sub_data["value"]) == len(sub_data["property"]): + for i, p in enumerate(sub_data["property"]): + self.psse_dict[sub_data["bus"]][sub_data["element_type"]][sub_data["element_id"]][p] = [] + logger.debug(f"clear self.psse_dict is {self.psse_dict}") + ###################################################### + # os.system("PAUSE") self.c_seconds_old = self.c_seconds return all_values diff --git a/pypsse/models.py b/pypsse/models.py index 514a236..c125395 100644 --- a/pypsse/models.py +++ b/pypsse/models.py @@ -33,6 +33,7 @@ UseModes, WritableModelTypes, ZoneProperties, + GenerationLevel, ) @@ -55,6 +56,13 @@ class SimSettings(BaseModel): user_models: List[str] = [] setup_files: List[str] = [] simulation_mode: SimulationModes + disable_generation_on_coupled_buses: bool = True + generation_model_level: GenerationLevel = GenerationLevel.TRANSMISSION + generation_std: str = "ieee2800" + transmission_ids: List[int] = [] + transmission_loads_at_fault: List = [] + transmission_loads_clear_fault: List = [] + transmission_loads_markup: bool = False @model_validator(mode="after") def sim_res_smaller_than_sim_time(self): @@ -171,9 +179,9 @@ class HelicsSettings(BaseModel): iterative_mode: bool = False error_tolerance: float = Field(1e-5, g=0) max_coiterations: int = Field(15, ge=1) - broker_ip: IPvAnyAddress = "127.0.0.1" + broker_ip: str="127.0.0.1" #IPvAnyAddress = "127.0.0.1" broker_port: int = 23404 - disable_generation_on_coupled_buses: bool = True + generation_model_level: str = "distribution" publications: List[PublicationDefination] @@ -484,10 +492,92 @@ class ApiWebSocketRequest(BaseModel): class Contingencies(BaseModel): contingencies: List[Union[BusFault, BusTrip, LineFault, LineTrip, MachineTrip]] -class ProfileMap(BaseModel): - id: str - bus: str - multiplier : float = 1 - normalize : bool = False - interpolate : bool = False +class REGCA1_data_model(BaseModel): + iqmax: float=1.5 + iqmin: float=-1.5 + Tg: float=0.02 + xf: float=0.25 + volim0:float=1.0 + volim1:float=1.2 + lvpnt0:float=0.2 + lvpnt1:float=0.8 + khv:float=1.0 + ccflag:bool=False + +class REECA1_data_model(BaseModel): + ppriorityflag:bool=False + Trv:float=0.02 + dbd1:float=-0.02 + dbd2:float=0.02 + Kqv:float=1.0 + iqh1:float=1.2 + iql1:float=-1.2 + Vref:float=1.0 + Tiq:float=0.02 + dpmin:float=-1.2 + dpmax:float=1.2 + Pmax:float=2 + Pmin:float=0 + imax:float=1.2 + Tpord:float=0.02 + ipmax:float=1.2 + ipmin:float=0 + iqmax:float=1.2 + iqmin:float=-1.2 + +class REPCA1_data_model(BaseModel): + refflag:int=0 + fflag:bool=True + Tfilt:float=0.02 + Kp:float=1 + Ki:float=1 + Tft:float=0.02 + Tfv:float=0.05 + emax:float=2 + emin:float=-2 + dbd1:float=-0.02 + dbd2:float=0.02 + qmax:float=2 + qmin:float=-2 + Kpg:float=1 + Kig:float=1 + fdbd1:float=-0.02 + fdbd2:float=0.02 + femax:float=1.2 + femin:float=0 + pmax:float=2 + pmin:float=0 + Tg:float=0.02 + + + +class FastDER_data_model(BaseModel): + Trv:float=0.02 + dbd1:float=-0.025 + dbd2:float=0.025 + Kqv:float=0.02 + Tiq:float=0.02 + Kpg:float=1.0 + Kig:float=1.0 + fdbd1:float=-0.03 + fdbd2:float=0.03 + Trf:float=0.02 + Tpord:float=0.02 + Tg:float=0.02 + iqmax:float=2.0 + iqmin:float=-2.0 + ipmax:float=2.0 + ipmin:float=-2.0 + femax:float=0.2 + femin:float=-0.2 + imax:float=2.0 + ppriorityflag:bool=True + xf: float=0.25 + +class ProfileMap(BaseModel): + id: str + bus: str + multiplier : float = 1 + normalize : bool = False + interpolate : bool = False diff --git a/pypsse/modes/abstract_mode.py b/pypsse/modes/abstract_mode.py index f0bbdac..60fdb0e 100644 --- a/pypsse/modes/abstract_mode.py +++ b/pypsse/modes/abstract_mode.py @@ -2,6 +2,7 @@ import numpy as np from loguru import logger +import os from pypsse.common import CASESTUDY_FOLDER, VALUE_UPDATE_BOUND from pypsse.enumerations import ModelTypes, WritableModelTypes @@ -14,6 +15,7 @@ def __init__( self, psse, dyntools, + der, settings: SimulationSettings, export_settings: ExportFileOptions, subsystem_buses, @@ -33,6 +35,7 @@ def __init__( self.sub_buses = subsystem_buses self.dyntools = dyntools + self.der = der self.settings = settings self.export_settings = export_settings self.func_options = { @@ -264,8 +267,9 @@ def save_model(self): i = 0 file_path = export_path / f"modified_steady_state_{i}.sav" while file_path.exists(): - file_path = export_path / f"modified_steady_state_{i}.sav" i += 1 + file_path = export_path / f"modified_steady_state_{i}.sav" + savfile = export_path / f"modified_steady_state_{i}.sav" rawfile = export_path / f"modified_steady_state_{i}.raw" @@ -680,9 +684,7 @@ def read_subsystems(self, quantities, subsystem_buses, ext_string2_info=None, ma elif func_name in ["loddt2", "lodnofunc", "lodint"]: ierr = self.psse.inilod(int(b)) - ierr, load_id = self.psse.nxtlod(int(b)) - while load_id is not None: if func_name == "lodnofunc": if v in ["BUSNUM", "LOADID"]: @@ -690,18 +692,14 @@ def read_subsystems(self, quantities, subsystem_buses, ext_string2_info=None, ma results = self.add_result(results, q, val, f"{b}_{load_id}") elif v == "BUSNAME": ierr, val = self.psse.notona(int(b)) - results = self.add_result(results, q, val, f"{b}_{load_id}") elif func_name == "loddt2": ierr, val = getattr(self.psse, func_name)(int(b), load_id, v, "ACT") results = self.add_result(results, q, val, f"{b}_{load_id}") elif func_name == "lodint": ierr, val = getattr(self.psse, func_name)(int(b), load_id, v) - results = self.add_result(results, q, val, f"{b}_{load_id}") - ierr, load_id = self.psse.nxtlod(int(b)) - elif func_name in ["macdat", "macdt2", "macnofunc", "macint"]: ierr = self.psse.inimac(int(b)) ierr, mach_id = self.psse.nxtmac(int(b)) @@ -715,33 +713,25 @@ def read_subsystems(self, quantities, subsystem_buses, ext_string2_info=None, ma results = self.add_result(results, q, val, f"{b}_{mach_id}") elif v == "SUBNUMBER": ierr, val = self.psse.busint(int(b), "STATION") - results = self.add_result(results, q, val, f"{b}_{mach_id}") - elif v == "AREANUMBER": ierr, val = self.psse.busint(int(b), "AREA") - results = self.add_result(results, q, val, f"{b}_{mach_id}") - elif v in ["SUBLATITUDE", "SUBLONGITUDE"]: sub_dict = {"SUBLATITUDE": "LATI", "SUBLONGITUDE": "LONG"} ierr, val = self.psse.busint(int(b), "STATION") if val: ierr, val = self.psse.stadat(val, sub_dict[v]) - results = self.add_result(results, q, val, b) else: ierr, val = getattr(self.psse, func_name)(int(b), mach_id, v) - results = self.add_result(results, q, val, f"{b}_{mach_id}") ierr, mach_id = self.psse.nxtmac(int(b)) elif func_name in ["fxsdt2", "fxsnofunc"]: ierr = self.psse.inifxs(int(b)) - ierr, fx_id = self.psse.nxtfxs(int(b)) - while fx_id is not None: if func_name == "fxsnofunc": if v in ["BUSNUM", "FXSHID"]: @@ -749,35 +739,28 @@ def read_subsystems(self, quantities, subsystem_buses, ext_string2_info=None, ma results = self.add_result(results, q, val, f"{b}_{fx_id}") elif v == "BUSNAME": ierr, val = self.psse.notona(int(b)) - results = self.add_result(results, q, val, f"{b}_{fx_id}") else: ierr, val = getattr(self.psse, func_name)(int(b), fx_id, v) - results = self.add_result(results, q, val, f"{b}_{fx_id}") ierr, fx_id = self.psse.nxtfxs(int(b)) elif func_name in ["swsdt1", "swsnofunc"]: if func_name == "swsnofunc": ierr, val = self.psse.swsint(int(b), "STATUS") - else: ierr, val = getattr(self.psse, func_name)(int(b), v) - if ierr == 0: if func_name == "nofunc": if v == "BUSNUM": val = int(b) elif v == "BUSNAME": ierr, val = self.psse.notona(int(b)) - results = self.add_result(results, q, val, b) elif func_name in ["brndat", "brndt2", "brnmsc", "brnint", "brnnofunc"]: ierr = self.psse.inibrx(int(b), 1) - ierr, b1, ickt = self.psse.nxtbrn(int(b)) - ickt_string = str(ickt) while ierr == 0: if func_name == "brnnofunc": @@ -836,7 +819,6 @@ def read_subsystems(self, quantities, subsystem_buses, ext_string2_info=None, ma else: ierr, val = getattr(self.psse, func_name)(int(b), int(b1), str(ickt), v) - if ierr == 0: results = self.add_result( results, q, val, f"{b!s}_{b1!s}_{ickt_string}" @@ -990,21 +972,28 @@ def convert_load(self, bus_subsystem=None): def update_object(self, dtype, bus, element_id, values: dict): val = sum(list(values.values())) if val > -VALUE_UPDATE_BOUND and val < VALUE_UPDATE_BOUND: - if dtype == WritableModelTypes.LOAD.value: + if dtype == WritableModelTypes.LOAD.value or dtype == WritableModelTypes.LOAD_STATUS.value: ierr = self.psse.load_chng_5(ibus=int(bus), id=element_id, **values) elif dtype == WritableModelTypes.GENERATOR.value: ierr = self.psse.induction_machine_data(ibus=int(bus), id=element_id, **values) - elif dtype == WritableModelTypes.MACHINE.value: + elif dtype == WritableModelTypes.MACHINE.value or dtype == WritableModelTypes.MACHINE_STATUS.value: ierr = self.psse.machine_data_2(i=int(bus), id=element_id, **values) elif dtype == WritableModelTypes.PLANT.value: - ierr = self.psse.plant_data_4(ibus=int(bus), inode=element_id, **values) + ierr = self.psse.plant_data_4(ibus=int(bus), inode=0, intgar=[self._i, self._i], **values) + elif dtype == WritableModelTypes.LINE_STATUS.value: + frombus, tobus = bus.split("_") + # old_load = self.psse.brnint(ibus=int(frombus), jbus=int(tobus), ickt=element_id, string='STATUS') + # logger.info(f"old line status {bus}-{element_id}: {old_load}") + ierr = self.psse.branch_data_3(ibus=int(frombus), jbus=int(tobus), ckt=element_id, **values) + # old_load = self.psse.brnint(ibus=int(frombus), jbus=int(tobus), ickt=element_id, string='STATUS') + # logger.info(f"new line status {bus}-{element_id}: {old_load}") else: ierr = 1 if ierr == 0: logger.info(f"Profile Manager: {dtype} '{element_id}' on bus '{bus}' has been updated. {values}") else: - logger.error(f"Profile Manager: Error updating {dtype} '{element_id}' on bus '{bus}'.") - + logger.error(f"Profile Manager: Error updating {dtype} '{element_id}' on bus '{bus}'. Error code is {ierr}") + def has_converged(self): return self.psse.solved() diff --git a/pypsse/modes/constants.py b/pypsse/modes/constants.py index ca91232..a722b68 100644 --- a/pypsse/modes/constants.py +++ b/pypsse/modes/constants.py @@ -205,7 +205,10 @@ def return_data_type(data_list, cls_type): status_translation = {1: "connected", 0: "notconnected"} if param.lower() == "status": for key, value in sub_dict.items(): - sub_dict[key] = status_translation[value] + if value == "connected" or "notconnected": + sub_dict[key] = value + else: + sub_dict[key] = status_translation[value] if class_name in complex_conversion_dict: conv_dict = complex_conversion_dict[class_name] diff --git a/pypsse/modes/dynamic.py b/pypsse/modes/dynamic.py index 80cf661..f8236cb 100644 --- a/pypsse/modes/dynamic.py +++ b/pypsse/modes/dynamic.py @@ -13,12 +13,13 @@ def __init__( self, psse, dyntools, + der, settings: SimulationSettings, export_settings: ExportFileOptions, subsystem_buses, raw_data, ): - super().__init__(psse, dyntools, settings, export_settings, subsystem_buses, raw_data) + super().__init__(psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data) self.time = settings.simulation.start_time self._StartTime = settings.simulation.start_time self.incTime = settings.simulation.simulation_step_resolution @@ -118,7 +119,9 @@ def step(self, t): self.time = self.time + self.incTime self.xTime = 0 - return self.psse.run(0, t, 1, 1, 1) + # self.psse.run(0, t + self.incTime.total_seconds(), 1, 1, 1) + # self.psse.run(0, t + self.incTime.total_seconds(), 1, 1, 1) + return self.psse.run(0, t + self.incTime.total_seconds(), 1, 1, 1) # @kapil do you need this? diff --git a/pypsse/modes/pcm.py b/pypsse/modes/pcm.py index caf5b60..ef01f51 100644 --- a/pypsse/modes/pcm.py +++ b/pypsse/modes/pcm.py @@ -10,8 +10,8 @@ class ProductionCostModel(AbstractMode): - def __init__(self, psse, dyntools, settings, export_settings, subsystem_buses, raw_data): - super().__init__(psse, dyntools, settings, export_settings, subsystem_buses, raw_data) + def __init__(self, psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data): + super().__init__(psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data) self.time = datetime.datetime.strptime(settings["Simulation"]["Start time"], "%m/%d/%Y %H:%M:%S").astimezone( None ) diff --git a/pypsse/modes/snap.py b/pypsse/modes/snap.py index f32d864..15875c9 100644 --- a/pypsse/modes/snap.py +++ b/pypsse/modes/snap.py @@ -1,5 +1,6 @@ import numpy as np from loguru import logger +import os from pypsse.models import SimulationSettings, ExportFileOptions from pypsse.modes.abstract_mode import AbstractMode @@ -8,18 +9,19 @@ class Snap(AbstractMode, DynamicUtils): - "Class defination for snapshat simulation mode (uses snp and sav files)" + "Class defination for snapshot simulation mode (uses snp and sav files)" def __init__( self, psse, dyntools, + der, settings: SimulationSettings, export_settings: ExportFileOptions, subsystem_buses, raw_data, ): - super().__init__(psse, dyntools, settings, export_settings, subsystem_buses, raw_data) + super().__init__(psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data) self.time = settings.simulation.start_time self._StartTime = settings.simulation.start_time self.incTime = settings.simulation.simulation_step_resolution @@ -43,8 +45,7 @@ def init(self, bus_subsystems): # The following logic only runs when the helics interface is enabled self.disable_load_models_for_coupled_buses() - #self.disable_generation_for_coupled_buses() - # self.save_model() + self.disable_generation_for_coupled_buses() ############# ------------------------------------- ############### outx_file = str(self.settings.export.outx_file).split("\\") outx_file[-1] = self.export_settings.filename_prefix + "_" + outx_file[-1] @@ -78,15 +79,19 @@ def init(self, bus_subsystems): self.psse.delete_all_plot_channels() self.setup_all_channels() - + # print("111111111111111111111111111111111111111111111111111111111111111111") logger.debug("pyPSSE initialization complete!") self.initialization_complete = True + # print("222222222222222222222222222222222222222222222222222222222222222222") return self.initialization_complete def step(self, t): "Increments the simulation" self.time = self.time + self.incTime self.xTime = 0 + if self.settings.simulation.disable_generation_on_coupled_buses and self.settings.simulation.generation_model_level.lower() == "transmission" and len(self.settings.simulation.transmission_ibrs) > 0: + # FX: if transmission DER model is disable and ibr is defined + self.der.update_ibr() return self.psse.run(0, t + self.incTime.total_seconds(), 1, 1, 1) def resolve_step(self, t): @@ -168,6 +173,6 @@ def read_subsystems(self, quantities, subsystem_buses, ext_string2_info=None, ma obj_name = f"{bus}_{ld_id}" results[res_base][obj_name] = value else: - logger.warning("Extend function 'read_subsystems' in the Snap class (Snap.py)") + logger.warning(f"{class_name} not in the option") return results diff --git a/pypsse/modes/static.py b/pypsse/modes/static.py index 0368aeb..ecf4c2d 100644 --- a/pypsse/modes/static.py +++ b/pypsse/modes/static.py @@ -8,10 +8,10 @@ class Static(AbstractMode): - def __init__(self, psse, dyntools, settings, export_settings, subsystem_buses, raw_data): + def __init__(self, psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data): "Class defination for steady-state simulation mode" - super().__init__(psse, dyntools, settings, export_settings, subsystem_buses, raw_data) + super().__init__(psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data) self.time = settings.simulation.start_time self._StartTime = settings.simulation.start_time self.incTime = settings.simulation.simulation_step_resolution diff --git a/pypsse/profile_manager/common.py b/pypsse/profile_manager/common.py index 638bc94..2f6433d 100644 --- a/pypsse/profile_manager/common.py +++ b/pypsse/profile_manager/common.py @@ -5,6 +5,9 @@ class ProfileTypes(str, Enum): INDUCTION_MACHINE = 'Induction_machine' MACHINE = 'Machine' PLANT = "Plant" + LOAD_STATUS = "Load_status" + LINE_STATUS = "Line_status" + MACHINE_STATUS = "Machine_status" PROFILE_VALIDATION = { ProfileTypes.LOAD: ["PL", "QL", "IP", "IQ", "YP", "YQ", "PG", "QG"], @@ -52,6 +55,9 @@ class ProfileTypes(str, Enum): "WPF", ], ProfileTypes.PLANT: ["VS", "RMPCT"], + ProfileTypes.LOAD_STATUS: ["STATUS"], + ProfileTypes.LINE_STATUS: ["STATUS"], + ProfileTypes.MACHINE_STATUS: ["STATUS"], } DEFAULT_PROFILE_NAME = "Default" diff --git a/pypsse/profile_manager/profile.py b/pypsse/profile_manager/profile.py index 06d1f9f..a742d15 100644 --- a/pypsse/profile_manager/profile.py +++ b/pypsse/profile_manager/profile.py @@ -8,7 +8,7 @@ class Profile: - "Class defination fora single profile" + "Class defination for single profile" DEFAULT_SETTINGS = {"multiplier": 1, "normalize": False, "interpolate": False} @@ -35,9 +35,12 @@ def __init__(self, profile_obj, solver, mapping_dict, buffer_size=10, neglect_ye def update(self, update_object_properties=True): "Returns value at the current timestep in the given profile" self.time = copy.deepcopy(self.solver.get_time()).astimezone(None) + logger.info(f"self.profile: {self.profile}") if self.time < self.stime or self.time > self.etime: value = np.array([0] * len(self.profile[0])) value1 = np.array([0] * len(self.profile[0])) + valuen1 = np.array([0] * len(self.profile[0])) + dt2 = self.attrs["resTime"] else: dt = (self.time - self.stime).total_seconds() n = int(dt / self.attrs["resTime"]) @@ -46,14 +49,14 @@ def update(self, update_object_properties=True): valuen1 = np.array(list(self.profile[n + 1])) except Exception as _: valuen1 = value - dt2 = ( self.time - (self.stime + datetime.timedelta(seconds=int(n * self.attrs["resTime"]))) ).total_seconds() - value1 = value + (valuen1 - value) * dt2 / self.attrs["resTime"] if update_object_properties: for obj_name in self.value_settings: + if self.value_settings[obj_name]["interpolate"]: + value1 = value + (valuen1 - value) * dt2 / self.attrs["resTime"] bus, object_id = obj_name.split("__") if self.value_settings[obj_name]["interpolate"]: value = value1 @@ -71,6 +74,11 @@ def update(self, update_object_properties=True): def fill_missing_values(self, value): "Fixes issues in profile data" - idx = [f"realar{PROFILE_VALIDATION[self.dtype].index(c) + 1}" for c in self.columns] + # idx = [f"realar{PROFILE_VALIDATION[self.dtype].index(c) + 1}" for c in self.columns] + ## FX: added to handle more data + if self.dtype in ["Load","Induction_machine","Machine","Plant"]: + idx = [f"realar{PROFILE_VALIDATION[self.dtype].index(c) + 1}" for c in self.columns] + else: + idx = [f"intgar{PROFILE_VALIDATION[self.dtype].index(c) + 1}" for c in self.columns] x = dict(zip(idx, list(value))) return x diff --git a/pypsse/profile_manager/profile_store.py b/pypsse/profile_manager/profile_store.py index f132c5d..b59c221 100644 --- a/pypsse/profile_manager/profile_store.py +++ b/pypsse/profile_manager/profile_store.py @@ -40,7 +40,7 @@ def __init__( file_path = settings.simulation.project_path / PROFILES_FOLDER / DEFAULT_PROFILE_STORE_FILENAME if file_path.exists(): - logger.info("Loading existing h5 store") + logger.info(f"Loading existing h5 store: {file_path}") self.store = h5py.File(file_path, mode) else: logger.info("Creating new h5 store") @@ -80,6 +80,7 @@ def setup_profiles(self): self.profiles[f"{group}/{profile_name}"] = Profile( grp[profile_name], self.solver, mapping_dict ) + logger.info(rf"Group {group} \ data set {profile_name} is added") else: logger.warning(rf"Group {group} \ data set {profile_name} not found in the h5 store") else: @@ -266,6 +267,7 @@ def create_metadata( info (str): profile info p_type (ProfileTypes): profile type """ + # logger.debug(start_time.strftime()) metadata = { "sTime": str(start_time), @@ -281,7 +283,7 @@ def create_metadata( } for key, value in metadata.items(): if isinstance(value, str): - value_mod = np.string_(value) + value_mod = np.bytes_(value) else: value_mod = value d_set.attrs[key] = value_mod @@ -294,10 +296,11 @@ def update(self) -> dict: """ results = {} + # logger.info(self.profiles) for profile_name, profile_obj in self.profiles.items(): result = profile_obj.update() results[profile_name] = result - + # logger.info(results) return results # def __del__(self): diff --git a/pypsse/profile_manager_interface.py b/pypsse/profile_manager_interface.py index 43036fa..6543207 100644 --- a/pypsse/profile_manager_interface.py +++ b/pypsse/profile_manager_interface.py @@ -1,5 +1,6 @@ from datetime import timedelta from pathlib import Path +import os from loguru import logger import pandas as pd @@ -37,6 +38,7 @@ def get_profiles(self) -> list: all_datasets = [] for model_type, model_info in self._toml_dict.items(): for profile_id, model_maps in model_info.items(): + logger.info(f"model_type: {model_type}, model_info: {model_info}") for model_map in model_maps: bus_id : str = model_map["bus"] model_id : str = model_map["id"] @@ -44,23 +46,67 @@ def get_profiles(self) -> list: norm: True = model_map.get("normalize") dataset = self._store[f"{model_type}/{profile_id}"] data = np.array(np.array(dataset).tolist()) - - if norm: - data_max = np.array(dataset.attrs["max"]) - data = data / data_max - if mult: - data = data * mult + # logger.info(f"data: {data[:10]}") + # os.system("PAUSE") + if model_type == "Load": + if norm: + data_max = np.array(dataset.attrs["max"]) + data = data / data_max + if mult: + data = data * mult - data = pd.DataFrame(data) - even_sum = data.iloc[:, ::2].sum(axis=1) - odd_sum = data.iloc[:, 1::2].sum(axis=1) - data = [even_sum, odd_sum] + data = pd.DataFrame(data) + P_even_sum = data.iloc[:, ::2].sum(axis=1) + Q_odd_sum = data.iloc[:, 1::2].sum(axis=1) + data = [P_even_sum, Q_odd_sum] + elif model_type == "Machine": + # Added by Fuhong, for a 'good' generation profile. Pgen and Qgen are strictly limited within the geneator maximum and minimum ranges. + data = pd.DataFrame(data) + P_gen = data.iloc[:, 0] + Q_gen = data.iloc[:, 1] + data = [P_gen, Q_gen] + + elif model_type == "Plant": + # Added by Fuhong + data = pd.DataFrame(data) + vs_plant = data.iloc[:, 0] + rmpct_plant = data.iloc[:, 1] + data = [vs_plant, rmpct_plant] + + elif model_type == "Induction_machine": + # TODO + pass + + elif model_type == "Load_status": + data = pd.DataFrame(data) + load_status = data.iloc[:, 0] + data = [load_status] + + elif model_type == "Line_status": + data = pd.DataFrame(data) + line_status = data.iloc[:, 0] + data = [line_status] + + elif model_type == "Machine_status": + data = pd.DataFrame(data) + machine_status = data.iloc[:, 0] + data = [machine_status] + else: + raise TypeError + final_df = pd.concat(data, axis=1) - final_df.columns = [f"{model_type}_{model_id}_{bus_id}_P", f"{model_type}_{model_id}_{bus_id}_Q"] + if model_type == "Load" or model_type == "Machine": + final_df.columns = [f"{model_type}_{model_id}_{bus_id}_P", f"{model_type}_{model_id}_{bus_id}_Q"] + elif model_type == "Plant": + final_df.columns = [f"{model_type}_{model_id}_{bus_id}_V", f"{model_type}_{model_id}_{bus_id}_RMPCT"] + elif model_type == "Load_status" or model_type == "Line_status" or model_type == "Machine_status": + final_df.columns = [f"{model_type}_{model_id}_{bus_id}_Status"] + else: + raise TypeError start_time = str(dataset.attrs["sTime"].decode('utf-8')) end_time = str(dataset.attrs["eTime"].decode('utf-8')) - timestep = timedelta(seconds=dataset.attrs["resTime"]) + timestep = timedelta(seconds=int(dataset.attrs["resTime"])) date_range =pd.date_range( start_time, end_time, diff --git a/pypsse/simulation_controller.py b/pypsse/simulation_controller.py index 84b0063..a27d60d 100644 --- a/pypsse/simulation_controller.py +++ b/pypsse/simulation_controller.py @@ -13,6 +13,7 @@ def sim_controller( psse: object, dyntools: object, + der: object, settings: SimulationSettings, export_settings: ExportFileOptions, subsystem_buses: dict, @@ -33,9 +34,11 @@ def sim_controller( """ sim_modes = {"Dynamic": Dynamic, "Steady-state": Static, "Snap": Snap, "ProductionCostModel": ProductionCostModel} - + # print("0000000000000000000000000000000000000000000000000000") sim = sim_modes[settings.simulation.simulation_mode.value]( - psse, dyntools, settings, export_settings, subsystem_buses, raw_data + psse, dyntools, der, settings, export_settings, subsystem_buses, raw_data ) + # print("33333333333333333333333333333333333333333333333333333") logger.debug(f"Simulator contoller of type {settings.simulation.simulation_mode.value} created") + # print("4444444444444444444444444444444444444444444444444444444") return sim diff --git a/pypsse/simulator.py b/pypsse/simulator.py index e41bc55..8cc508a 100644 --- a/pypsse/simulator.py +++ b/pypsse/simulator.py @@ -17,7 +17,7 @@ import toml from loguru import logger from networkx import Graph - +import pssepath import pypsse.contingencies as c import pypsse.simulation_controller as sc from pypsse.common import ( @@ -48,13 +48,16 @@ from pypsse.result_container import Container from pypsse.models import Contingencies from pypsse.enumerations import PSSE_VERSIONS +from pypsse.utils.dynamic_utils import DynamicUtils + +from aggregatedderapp import App USING_NAERM = 0 N_BUS = 200000 -class Simulator: +class Simulator(DynamicUtils): "Base class for the simulator" _status: SimulationStatus = SimulationStatus.NOT_INITIALIZED @@ -77,6 +80,7 @@ def __init__( self.settings = settings logger.debug(f"Instantiating psse version {psse_version}") + pssepath.add_pssepath(35.4) __import__(psse_version, fromlist=[""]) # noqa: F401 import dyntools @@ -105,6 +109,7 @@ def __init__( self.dyntools = dyntools self.psse = psspy + self.der = App(showProgress=True) # logger.debug('Initializing PSS/E. connecting to license server') self.start_simulation() @@ -188,6 +193,7 @@ def start_simulation(self): self.sim = sc.sim_controller( self.psse, self.dyntools, + self.der, self.settings, self.export_settings, valid_buses, @@ -322,7 +328,7 @@ def run(self): bokeh_server_proc = None logger.debug( - f"Running {self.settings.simulation.simulation_mode.value} simulation for time {self.settings.simulation.simulation_time.total_seconds()} sec" + f"Running {self.settings.simulation.simulation_mode.value} simulation for time {self.settings.simulation.simulation_time.total_seconds()} sec" ) total_simulation_time = ( self.settings.simulation.simulation_time.total_seconds() @@ -376,25 +382,41 @@ def step(self, t: float) -> dict: if self.settings.simulation.use_profile_manager: self.pm.update() ctime = time.time() - self.simStartTime + ierr, rval = self.psse.dsrval('TIME', 0) logger.debug( - f"Simulation time: {t} seconds\nRun time: {ctime}\npsse time: {self.sim.get_time()}" + f"Simulation time: {t} seconds\nRun time: {ctime}\npsse time: {self.sim.get_time()} \nsolver time: {rval}" ) if self.settings.helics and self.settings.helics.cosimulation_mode: - if self.settings.helics.create_subscriptions: + if self.settings.helics.create_subscriptions: + # logger.debug(f"Time requested: {t}") + # self.inc_time, helics_time = self.update_federate_time(t) + # logger.debug(f"Time granted: {helics_time}") self.update_subscriptions() - logger.debug(f"Time requested: {t}") - self.inc_time, helics_time = self.update_federate_time(t) - logger.debug(f"Time granted: {helics_time}") if self.inc_time: + logger.debug(f"run PSSE simulation at time: {t}") self.sim.step(t) + # os.system("PAUSE") else: self.sim.resolve_step() if self.settings.helics and self.settings.helics.cosimulation_mode: self.publish_data() + if self.settings.helics and self.settings.helics.cosimulation_mode: + if self.settings.helics.create_subscriptions: + # self.update_subscriptions() + logger.debug(f"Time requested: {t}") + self.inc_time, helics_time = self.update_federate_time(t) + logger.debug(f"Time granted: {helics_time}") + curr_results = self.update_result_container(t) + # if t >= 0.2: + # The following logic only runs when the helics interface is enabled + # self.disable_load_models_for_coupled_buses() + # self.disable_generation_for_coupled_buses() + # os.system("PAUSE") + return curr_results def update_result_container(self, t: float) -> dict: diff --git a/pypsse/utils/dynamic_utils.py b/pypsse/utils/dynamic_utils.py index 37013a4..a6b4aaa 100644 --- a/pypsse/utils/dynamic_utils.py +++ b/pypsse/utils/dynamic_utils.py @@ -6,7 +6,9 @@ from pypsse.common import MACHINE_CHANNELS from pypsse.modes.constants import dyn_only_options - +import re +import toml +from pypsse.models import REGCA1_data_model, REECA1_data_model, REPCA1_data_model class DynamicUtils: "Utility functions for dynamic simulations" @@ -15,11 +17,8 @@ class DynamicUtils: def disable_generation_for_coupled_buses(self): """Disables generation of coupled buses (co-simulation mode only)""" - if ( - self.settings.helics - and self.settings.helics.cosimulation_mode - and self.settings.helics.disable_generation_on_coupled_buses - ): + if ((self.settings.helics and self.settings.helics.cosimulation_mode and self.settings.simulation.disable_generation_on_coupled_buses) # cosim mode, load in distribution level + or (self.settings.simulation.disable_generation_on_coupled_buses and self.settings.simulation.generation_model_level.lower() == "transmission")): sub_data = pd.read_csv(self.settings.simulation.subscriptions_file) sub_data = sub_data[sub_data["element_type"] == "Load"] generators = {} @@ -32,7 +31,10 @@ def disable_generation_for_coupled_buses(self): for _, row in sub_data.iterrows(): bus = row["bus"] - generators[bus] = generator_list[bus] + if bus in generator_list: + generators[bus] = generator_list[bus] + else: + logger.warning(f"No generators at coupled bus {bus}; skipping generation disable for this bus.") for bus_id, machines in generators.items(): for machine in machines: @@ -58,15 +60,90 @@ def disable_generation_for_coupled_buses(self): ] self.psse.machine_chng_2(bus_id, machine, intgar, realar) logger.info(f"Machine disabled: {bus_id}_{machine}") + + if self.settings.simulation.generation_model_level.lower() == "transmission" and len(self.settings.simulation.transmission_ibrs) > 0: + intgar = [0, self._i, self._i, self._i, self._i, self._i] + realar = [ + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + self._f, + ] + transmission_ibrs_id = self.settings.simulation.transmission_ibrs + transmission_ibrs_std = self.settings.simulation.transmission_ibrs_std + assert len(transmission_ibrs_id)==len(transmission_ibrs_std), f"Input dimension mismatch" + # added by FX for IBR temporarily + ibr_setting_path = r"C:\Users\FXIE\Documents\GitHub\NEARM_TnD\dynamic_cosim_test_cases-model\WECC\Transmission\pypsse_model\Scenarios\IEEE2800\IBR.toml" + ibr_settings = toml.load(ibr_setting_path) + for ibr, std in zip(transmission_ibrs_id, transmission_ibrs_std): + # generation is still in transmission level but use different modeling method + if std.lower() == "ieee2800": + logger.debug(f"Added IBR with IEEE STD 2800 for generator {ibr}") + bus_id, machine = ibr.split("_") + ierr, machine_pg = self.psse.macdat(int(bus_id),machine,'P') + ierr, machine_qg = self.psse.macdat(int(bus_id),machine,'Q') + ierr, machine_base = self.psse.macdat(int(bus_id),machine,'MBASE') + ierr, machine_status = self.psse.macint(int(bus_id),machine,'STATUS') + if ibr in ibr_settings.keys(): + ibr_setting = ibr_settings[ibr] + logger.debug(f"REPCA1: {ibr_setting['REPCA1']}") + param_repca1 = REPCA1_data_model(**ibr_setting["REPCA1"]) + param_reeca1 = REECA1_data_model(**ibr_setting["REECA1"]) + param_regca1 = REGCA1_data_model(**ibr_setting["REGCA1"]) + else: + param_repca1 = None + param_reeca1 = None + param_regca1 = None + if ierr == 0 and machine_status == 1: + self.psse.machine_chng_2(int(bus_id), machine, intgar, realar) + ibr_id = f"ir" + logger.info(f"IBR added: {bus_id}_{machine} (pg={machine_pg},qg={machine_qg},base={machine_base})") + ibr_dt = self.settings.simulation.simulation_step_resolution.total_seconds() + self.der.add_ibr([int(bus_id)],[machine_pg],[machine_qg],[ibr_id],ibr_dt,param_repca1,param_reeca1,param_regca1) + else: + logger.warning(f"Can't add IBR: {bus_id}_{machine})") + else: + # todo + raise Exception("Standard has not been implemented") + logger.debug(f"self.der.ibr['model']: {self.der.ibr['model']}") + + # + # os.system("PAUSE") + + def update_loadchannel_asset_list(self,assset_type, asset): + for n in range(len(self.export_settings.channel_setup)): + channel = self.export_settings.channel_setup[n] + method_type = channel.asset_type + if method_type == assset_type: + logger.info(channel) + asset_list = channel.asset_list + if asset not in asset_list: + asset_list.append(asset) + self.export_settings.channel_setup[n].asset_list = asset_list + def disable_load_models_for_coupled_buses(self): """Disables loads of coupled buses (co-simulation mode only)""" if self.settings.helics and self.settings.helics.cosimulation_mode: sub_data = pd.read_csv(self.settings.simulation.subscriptions_file) sub_data = sub_data[sub_data["element_type"] == "Load"] + ### add by Aadil logger.debug("Implementing load brek logic for dynamic cosimulations") self.break_loads_for_dynamic_cosimulations(loads=None, components_to_replace=["FmD"]) - + self.psse_dict = {} for _, row in sub_data.iterrows(): bus = row["bus"] @@ -77,8 +154,8 @@ def disable_load_models_for_coupled_buses(self): elif ierr == 5: logger.error(f"No dynamic model found for load {load} connected to bus {bus}") else: - raise Exception(f"error={ierr}") - + raise Exception(f"error={ierr}") + def break_loads_for_dynamic_cosimulations(self, loads: list = None, components_to_replace: List[str] = []): """Implements the load split logic @@ -86,8 +163,11 @@ def break_loads_for_dynamic_cosimulations(self, loads: list = None, components_t loads (list, optional): list of coupled loads. Defaults to None. components_to_replace (List[str], optional): components to be simulated on distribution side. Defaults to []. """ - + logger.info("##################### break_loads #########################") components_to_stay = [x for x in self.dynamic_params if x not in components_to_replace] + logger.info(f"components_to_stay: {components_to_stay}") + logger.info(f"components_to_replace: {components_to_replace}") + # os.system("PAUSE") if loads is None: loads = self._get_coupled_loads() logger.debug("Fetching static data for loads") @@ -98,6 +178,7 @@ def break_loads_for_dynamic_cosimulations(self, loads: list = None, components_t loads = self._replicate_coupled_load(loads, components_to_replace) logger.debug("Updating dynamic parameters for load models") self._update_dynamic_parameters(loads, components_to_stay, components_to_replace) + # os.system("PAUSE") def _update_dynamic_parameters(self, loads: dict, components_to_stay: list, components_to_replace: list): """Updates dynamic parameters of composite old / replicated load models @@ -107,19 +188,22 @@ def _update_dynamic_parameters(self, loads: dict, components_to_stay: list, comp components_to_stay (list): components to be simulated on transmission side components_to_replace (list): components to be simulated on distribution side """ - + logger.info(f"_update_dynamic_parameters") new_percentages = {} for load in loads: count = 0 for comp in components_to_stay: count += load[comp] for comp in components_to_stay: - new_percentages[comp] = load[comp] / count + if count == 0: + new_percentages[comp] = 0 + else: + new_percentages[comp] = load[comp] / count for comp in components_to_replace: new_percentages[comp] = 0.0 settings = self._get_load_dynamic_properties(load) - # + for k, v in new_percentages.items(): idx = dyn_only_options["Loads"]["lmodind"][k] settings[idx] = v @@ -127,7 +211,14 @@ def _update_dynamic_parameters(self, loads: dict, components_to_stay: list, comp # self.psse.change_ldmod_con(load['bus'], 'XX' ,r"""CMLDBLU2""" ,idx ,v) values = list(settings.values()) self.psse.add_load_model(load["bus"], "XX", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) - logger.info(f"Dynamic model parameters for load {load['id']} at bus 'XX' changed.") + # self.psse.add_load_model(load["bus"], "A", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) + # self.psse.add_load_model(load["bus"], "B", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) + # self.psse.add_load_model(load["bus"], "C", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) + # self.psse.add_load_model(load["bus"], "D", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) + # self.psse.add_load_model(load["bus"], "F", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) + # self.psse.add_load_model(load["bus"], "S", 0, 1, r"""CMLDBLU2""", 2, [0, 0], ["", ""], 133, values) + logger.info(f"Dynamic model parameters for load {load['id']} at bus 'XX' changed. New settings are {values}.") + def _get_load_dynamic_properties(self, load): "Returns dynamic parameters of composite load models" @@ -140,7 +231,36 @@ def _get_load_dynamic_properties(self, load): assert ierr == 0, f"error={ierr}" settings[i] = value return settings - + + def _get_bus_generation(self, bus): + "Returns the total generation on an certain bus" + generator_list = {} + generators = {} + for gen_bus, gen_id in self.raw_data.generators: + if gen_bus not in generator_list: + generator_list[gen_bus] = [] + generator_list[gen_bus].append(gen_id) + try: + generators[bus] = generator_list[bus] + except: + raise Exception(f"Can't get generator on bus {bus}") + # logger.warning(f"Can't get generator on bus {bus}") + bus_total_p = 0 + bus_total_q = 0 + for bus_id, machines in generators.items(): + for machine in machines: + ierr, ival = self.psse.macint(bus_id, machine, 'STATUS') + if ival == 1: + ierr, cmpval = self.psse.macdt2(bus_id, machine, 'PQ') + assert ierr == 0, f"error={ierr}" + logger.debug(f"{bus_id}.{machine} generation: {cmpval}") + else: + cmpval = complex(0,0) + logger.debug(f"{bus_id}.{machine} offline generation: {cmpval}") + bus_total_p += cmpval.real + bus_total_q += cmpval.imag + return bus_total_p, bus_total_q + def _replicate_coupled_load(self, loads: dict, components_to_replace: list): """create a replica of composite load model @@ -153,12 +273,17 @@ def _replicate_coupled_load(self, loads: dict, components_to_replace: list): """ for load in loads: + logger.info(f"load : {load}") dynamic_percentage = load["FmA"] + load["FmB"] + load["FmC"] + load["FmD"] + load["Fel"] static_percentage = 1.0 - dynamic_percentage + logger.debug(f"Static properties - Load: Static -> value:{static_percentage}") for comp in components_to_replace: static_percentage += load[comp] remaining_load = 1 - static_percentage - total_load = load["MVA"] + # Use TOTAL (MVA + IL + YL) so that constant-current and + # constant-admittance portions are included in the split. + # Using only MVA would drop the I and Y components. + total_load = load["TOTAL"] total_distribution_load = total_load * static_percentage total_transmission_load = total_load * remaining_load # ceate new load @@ -166,26 +291,40 @@ def _replicate_coupled_load(self, loads: dict, components_to_replace: list): load["bus"], "XX", realar=[total_transmission_load.real, total_transmission_load.imag, 0.0, 0.0, 0.0, 0.0], - lodtyp="replica", - ) - # ierr, cmpval = self.psse.loddt2(load["bus"], "XX" ,"MVA" , "ACT") - # modify old load - self.psse.load_data_5( - load["bus"], - str(load["id"]), - realar=[total_distribution_load.real, total_distribution_load.imag, 0.0, 0.0, 0.0, 0.0], - lodtyp="original", + # lodtyp="replica", ) - # ierr, cmpval = self.psse.loddt2(load["bus"], load["id"] ,"MVA" , "ACT") - logger.info(f"Original load {load['id']} @ bus {load['bus']}: {total_load}") - logger.info(f"New load 'XX' @ bus {load['bus']} created successfully: {total_transmission_load}") - logger.info(f"Load {load['id']} @ bus {load['bus']} updated : {total_distribution_load}") - load["distribution"] = total_distribution_load - load["transmission"] = total_transmission_load + if (self.settings.helics and self.settings.helics.cosimulation_mode and self.settings.simulation.disable_generation_on_coupled_buses + and self.settings.helics.generation_model_level == 'distribution'): + total_bus_generation_p, total_bus_generation_q = self._get_bus_generation(load['bus']) + logger.info(f"Generation is modeled in distribution level so transmission load is substituted the generation") + self.psse.load_data_5( + load["bus"], + str(load["id"]), + realar=[total_distribution_load.real-total_bus_generation_p, total_distribution_load.imag-total_bus_generation_q, 0.0, 0.0, 0.0, 0.0], + # lodtyp="original", + ) + logger.info(f"Original load {load['id']} @ bus {load['bus']}: {total_load}") + logger.info(f"New load 'XX' @ bus {load['bus']} created successfully: {total_transmission_load}") + logger.info(f"Load {load['id']} @ bus {load['bus']} updated : ({total_distribution_load.real-total_bus_generation_p},{total_distribution_load.imag-total_bus_generation_q})") + load["distribution"] = complex(total_distribution_load.real-total_bus_generation_p, total_distribution_load.imag-total_bus_generation_q) + load["transmission"] = total_transmission_load + else: + self.psse.load_data_5( + load["bus"], + str(load["id"]), + realar=[total_distribution_load.real, total_distribution_load.imag, 0.0, 0.0, 0.0, 0.0], + # lodtyp="original", + ) + logger.info(f"Original load {load['id']} @ bus {load['bus']}: {total_load}") + logger.info(f"New load 'XX' @ bus {load['bus']} created successfully: {total_transmission_load}") + logger.info(f"Load {load['id']} @ bus {load['bus']} updated : {total_distribution_load}") + load["distribution"] = total_distribution_load + load["transmission"] = total_transmission_load + logger.info(f"{loads}") return loads def _get_coupled_loads(self) -> list: - """Returns a list of all coupled loads ina give simualtion + """Returns a list of all coupled loads in a give simualtion Returns: list: list of coupled loads @@ -193,6 +332,7 @@ def _get_coupled_loads(self) -> list: sub_data = pd.read_csv(self.settings.simulation.subscriptions_file) load = [] + logger.info(f"_get_coupled_loads") for _, row in sub_data.iterrows(): if row["element_type"] == "Load": load.append( @@ -202,6 +342,7 @@ def _get_coupled_loads(self) -> list: "bus": row["bus"], } ) + logger.info(f"load is {load}") return load def _get_load_static_data(self, loads: list) -> dict: @@ -213,15 +354,13 @@ def _get_load_static_data(self, loads: list) -> dict: Returns: dict: mapping load to static values """ - + logger.info(f"_get_load_static_data") values = ["MVA", "IL", "YL", "TOTAL"] for load in loads: for v in values: ierr, cmpval = self.psse.loddt2(load["bus"], str(load["id"]), v, "ACT") load[v] = cmpval - logger.debug(f"Static properties - Load: {v} -> {cmpval}") - return loads def _get_load_dynamic_data(self, loads: list) -> dict: @@ -233,9 +372,11 @@ def _get_load_dynamic_data(self, loads: list) -> dict: Returns: dict: mapping load to dynamic values """ - + logger.info(f"_get_load_dynamic_data") values = dyn_only_options["Loads"]["lmodind"] + loads_with_dynamic_data = [] for load in loads: + has_dynamic = True for v, con_ind in values.items(): ierr = self.psse.inilod(load["bus"]) assert ierr == 0, f"error={ierr}" @@ -243,15 +384,22 @@ def _get_load_dynamic_data(self, loads: list) -> dict: assert ierr == 0, f"error={ierr}" if ld_id is not None: ierr, con_index = self.psse.lmodind(load["bus"], ld_id, "CHARAC", "CON") - assert ierr == 0, f"error={ierr}" + if ierr != 0: + logger.warning( + f"No CHARAC load model at bus {load['bus']} load '{ld_id}' " + f"(lmodind ierr={ierr}). Skipping dynamic data for this load." + ) + has_dynamic = False + break if con_index is not None: act_con_index = con_index + con_ind ierr, value = self.psse.dsrval("CON", act_con_index) assert ierr == 0, f"error={ierr}" load[v] = value logger.debug(f"Dynamic properties - Load: {v} -> index: {act_con_index}, value:{value}") - - return loads + if has_dynamic: + loads_with_dynamic_data.append(load) + return loads_with_dynamic_data def setup_machine_channels(self, machines: dict, properties: list): """sets up machine channels @@ -269,7 +417,7 @@ def setup_machine_channels(self, machines: dict, properties: list): if qty in MACHINE_CHANNELS: self.channel_map[nqty][f"{b}_{mch}"] = [self.chnl_idx] chnl_id = MACHINE_CHANNELS[qty] - logger.info(f"{qty} for machine {b}_{mch} added to channel {self.chnl_idx}") + # logger.info(f"{qty} for machine {b}_{mch} added to channel {self.chnl_idx}") self.psse.machine_array_channel([self.chnl_idx, chnl_id, int(b)], mch, "") self.chnl_idx += 1 @@ -307,12 +455,12 @@ def setup_bus_channels(self, buses: list, properties: list): if qty == "frequency": self.channel_map[qty][b] = [self.chnl_idx] self.psse.bus_frequency_channel([self.chnl_idx, int(b)], "") - logger.info(f"Frequency for bus {b} added to channel { self.chnl_idx}") + # logger.info(f"Frequency for bus {b} added to channel { self.chnl_idx}") self.chnl_idx += 1 elif qty == "voltage_and_angle": self.channel_map[qty][b] = [self.chnl_idx, self.chnl_idx + 1] self.psse.voltage_and_angle_channel([self.chnl_idx, -1, -1, int(b)], "") - logger.info(f"Voltage and angle for bus {b} added to channel {self.chnl_idx} and {self.chnl_idx+1}") + # logger.info(f"Voltage and angle for bus {b} added to channel {self.chnl_idx} and {self.chnl_idx+1}") self.chnl_idx += 2 def poll_channels(self) -> dict: @@ -355,6 +503,8 @@ def setup_all_channels(self): elif method_type == "loads": load_list = [[x, int(y)] for x, y in channel.asset_list] self.setup_load_channels(load_list) + # logger.debug(f"load_list: {load_list}") + # os.system("PAUSE") elif method_type == "machines": machine_list = [[x, int(y)] for x, y in channel.asset_list] self.setup_machine_channels(machine_list, channel.asset_properties) diff --git a/tests/examples/static_example/simulation_settings.toml b/tests/examples/static_example/simulation_settings.toml index 4631004..5ac7f02 100644 --- a/tests/examples/static_example/simulation_settings.toml +++ b/tests/examples/static_example/simulation_settings.toml @@ -5,7 +5,7 @@ psse_solver_timestep = 0.004166667 start_time = "2020-01-01 00:00:00.0" simulation_mode = "Steady-state" use_profile_manager = false -psse_path = "C:/Program Files/PTI/PSSE35/35.4/PSSPY39" +psse_path = "C:/Program Files/PTI/PSSE35/35.6/PSSPY39" project_path = "." case_study = "savnw.sav" raw_file = "savnw.raw"