Merge pull request #57 from NatLabRockies/develop_fx_bugfix

enables cosimulation with hpc settings
allows non-composite loads
starts on steady state enabling

Author: nadiavp

pypsse/cli/explore.py

@@ -1,5 +1,5 @@
 from pathlib import Path
-
+import os
 from loguru import logger
 import pandas as pd
 import click
@@ -129,7 +129,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 +149,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 +168,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 +200,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()}")
 

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()
-    

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,13 @@ def run(project_path, simulations_file=None):
 
     simulation_settiings = toml.load(file_path)
     simulation_settiings = SimulationSettings(**simulation_settiings)
-    
     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)
+        if simulation_settiings.log.clear_old_log_file:
+            logger.add(log_path, mode="w")
+        else:
+            logger.add(log_path)
 
     x = Simulator.from_setting_files(file_path)
 

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(

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()

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"

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
+        #################################################################
+        # add these hardcored load values for 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}."
@@ -133,6 +155,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 +258,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 +277,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 +398,23 @@ def subscribe(self) -> dict:
         """
 
         "Subscribes results each iteration and updates PSSE objects accordingly"
+
         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 +428,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}")
+                    ######################################################
+                    ## 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")
+        
+        ######################################################
+        ## 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