diff --git a/Dockerfile b/Dockerfile index 2666c64b..e81c55bd 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,4 @@ -FROM python:3.7.0-slim +FROM python:3.11.0-slim RUN apt-get update diff --git a/README.md b/README.md index bdbe0fba..a5765236 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,6 @@ # Welcome to the pydss Repository! + **PyDSS** is a high level python interface for **OpenDSS** and provides the following functionalities Documentation on installation, setup and examples can be found here https://natlabrockies.github.io/PyDSS/index.html diff --git a/src/pydss/SolveMode.py b/src/pydss/SolveMode.py index 55c3b069..4ae08157 100644 --- a/src/pydss/SolveMode.py +++ b/src/pydss/SolveMode.py @@ -10,7 +10,8 @@ def GetSolver(settings: SimulationSettingsModel, dssInstance): - logger.info('Setting solver to %s mode.', settings.project.simulation_type.value) + # logger.info('Setting solver to %s mode.', settings.project.simulation_type.value) + logger.info(f'Setting solver to {settings.project.simulation_type.value} mode.') return get_solver_from_simulation_type(settings.project) diff --git a/src/pydss/cli/run.py b/src/pydss/cli/run.py index 4df5b9e5..7096a048 100644 --- a/src/pydss/cli/run.py +++ b/src/pydss/cli/run.py @@ -5,6 +5,7 @@ from pathlib import Path import ast import sys +import os from loguru import logger import click @@ -80,6 +81,11 @@ def run(project_path, options=None, tar_project=False, zip_project=False, verbos logger.error("Logs path %s does not exist", logs_path) sys.exit(1) filename = logs_path / "pydss.log" + if settings.logging.clear_old_log_file: + logs_path = project_path / "Logs" + filename = logs_path / "pydss.log" + if os.path.exists(filename): + os.remove(filename) logger.level(console_level) if filename: @@ -92,8 +98,9 @@ def run(project_path, options=None, tar_project=False, zip_project=False, verbos if not isinstance(options, dict): logger.error("options are invalid: %s", options) sys.exit(1) - + # logger.info(f"indicator 1") project = PyDssProject.load_project(project_path, options=options, simulation_file=simulations_file) + # logger.info(f"indicator 2") project.run(tar_project=tar_project, zip_project=zip_project, dry_run=dry_run) if dry_run: diff --git a/src/pydss/common.py b/src/pydss/common.py index fd1082c6..5d0aecd4 100644 --- a/src/pydss/common.py +++ b/src/pydss/common.py @@ -37,7 +37,6 @@ class ControllerType(enum.Enum): STORAGE_CONTROLLER = "StorageController" THERMOSTATIC_LOAD_CONTROLLER = "ThermostaticLoad" XMFR_CONTROLLER = "xmfrController" - DYNAMIC_VOLTAGE_SUPPORT = "DynamicVoltageSupport" CONTROLLER_TYPES = tuple(x.value for x in ControllerType) CONFIG_EXT = ".toml" diff --git a/src/pydss/dataset_buffer.py b/src/pydss/dataset_buffer.py index 0b60602c..be2f991b 100644 --- a/src/pydss/dataset_buffer.py +++ b/src/pydss/dataset_buffer.py @@ -37,6 +37,7 @@ def __init__( max_chunk_bytes=None, attributes=None, names=None, column_ranges_per_name=None, data=None ): + if max_chunk_bytes is None: max_chunk_bytes = DEFAULT_MAX_CHUNK_BYTES self._buf_index = 0 diff --git a/src/pydss/dssElement.py b/src/pydss/dssElement.py index 377c3fa0..8396e2b6 100644 --- a/src/pydss/dssElement.py +++ b/src/pydss/dssElement.py @@ -1,5 +1,5 @@ import ast - +from loguru import logger from opendssdirect import DSSException from pydss.dssBus import dssBus @@ -60,9 +60,10 @@ def __init__(self, dssInstance): super(dssElement, self).__init__(dssInstance, name, fullName) self._Enabled = dssInstance.CktElement.Enabled() if not self._Enabled: + logger.debug(f"Element isn't defined: {fullName}") return - self._Parameters = {} + logger.debug(fullName) self._NumTerminals = dssInstance.CktElement.NumTerminals() self._NumConductors = dssInstance.CktElement.NumConductors() @@ -135,8 +136,12 @@ def DataLength(self, VarName): return 0, None def GetValue(self, VarName, convert=False): + if self._dssInstance.Element.Name() != self._FullName: self.SetActiveObject() + fullName = self._dssInstance.Element.Name() + # logger.debug("123456789123456789123456789123456789123456789123456789") + # logger.debug(fullName) if VarName in self._Variables: VarValue = self.GetVariable(VarName, convert=convert) elif VarName in self._Parameters: diff --git a/src/pydss/dssInstance.py b/src/pydss/dssInstance.py index 3a08d6c5..86699dab 100644 --- a/src/pydss/dssInstance.py +++ b/src/pydss/dssInstance.py @@ -156,7 +156,11 @@ def _compile_model(self): def _CreateControllers(self, ControllerDict): self._pyControls = {} self._pyControls_types = {} + # logger.info(f'self._dssObjects -> {self._dssObjects}') + # os.system("PAUSE") for ControllerType, ElementsDict in ControllerDict.items(): + # logger.info(f'ControllerType -> {ControllerType}, ElementsDict -> {ElementsDict}') + # os.system("PAUSE") for ElmName, SettingsDict in ElementsDict.items(): Controller = pyControllers.pyController.Create(ElmName, ControllerType, SettingsDict, self._dssObjects, self._dssInstance, self._dssSolver) @@ -205,7 +209,8 @@ def CreateDssObjects(dssBuses): InvalidSelection = ['Settings', 'ActiveClass', 'dss', 'utils', 'PDElements', 'XYCurves', 'Bus', 'Properties'] # TODO: this causes a segmentation fault. Aadil says it may not be needed. #self._dssObjectsByClass={'LoadShape': self._get_relavent_object_dict('LoadShape')} - + # logger.info(f"dss.Circuit.AllElementNames() -> {dss.Circuit.AllElementNames()}") + # os.system("PAUSE") for ElmName in dss.Circuit.AllElementNames(): Class, Name = ElmName.split('.', 1) ClassName = Class + 's' @@ -246,6 +251,10 @@ def RunStep(self, step, updateObjects=None): if self._settings.profiles.use_profile_manager: self.profileStore.update() + if self._settings.helics.co_simulation_mode: + # self._heilcs_interface.updateHelicsPublications() + self._increment_flag, helics_time = self._heilcs_interface.request_time_increment() + if self._settings.helics.co_simulation_mode: self._heilcs_interface.updateHelicsSubscriptions() else: @@ -253,6 +262,17 @@ def RunStep(self, step, updateObjects=None): for object, params in updateObjects.items(): cl, name = object.split('.') self._Modifier.Edit_Element(cl, name, params) + + # pydss_update_agian = "y" + # while pydss_update_agian == "y": + # if self._settings.helics.co_simulation_mode: + # self._heilcs_interface.updateHelicsSubscriptions() + # else: + # if updateObjects: + # for object, params in updateObjects.items(): + # cl, name = object.split('.') + # self._Modifier.Edit_Element(cl, name, params) + # pydss_update_agian = input('enter your pydss update command: ') # run simulation time step and get results time_step_has_converged = True @@ -263,6 +283,7 @@ def RunStep(self, step, updateObjects=None): for i in range(self._settings.project.max_control_iterations): has_converged, error = self._update_controllers(priority, step, i, UpdateResults=False) logger.debug('Control Loop {} convergence error: {}'.format(priority, error)) + logger.debug('Control Loop {} convergence @step {} '.format(priority, step)) if has_converged: priority_has_converged = True break @@ -297,8 +318,9 @@ def RunStep(self, step, updateObjects=None): if self._settings.helics.co_simulation_mode: self._heilcs_interface.updateHelicsPublications() - self._increment_flag, helics_time = self._heilcs_interface.request_time_increment() - + # if step < 1: + # self._increment_flag, helics_time = self._heilcs_interface.request_time_increment_2() + # os.system("PAUSE") return time_step_has_converged def _HandleConvergenceErrorChecks(self, step, error): @@ -316,7 +338,7 @@ def _HandleOpenDSSConvergenceErrorChecks(self, step): self._convergenceErrorsOpenDSS += 1 if self._maxConvergenceErrorCount is not None and self._convergenceErrorsOpenDSS > self._maxConvergenceErrorCount: - logger.error("Exceeded OpenDSS convergence error count threshold at step %s", step) + logger.error(f"Exceeded OpenDSS convergence error count threshold at step {step}") raise OpenDssConvergenceErrorCountExceeded(f"{self._convergenceErrorsOpenDSS} errors occurred") def DryRunSimulation(self, project, scenario): @@ -357,7 +379,7 @@ def RunSimulation(self, project, scenario, MC_scenario_number=None): dss.Solution.Convergence(self._settings.project.error_tolerance) logger.info('Running simulation from {} till {}.'.format(sTime, eTime)) logger.info('Simulation time step {}.'.format(Steps)) - logger.info("Set OpenDSS convergence to %s", dss.Solution.Convergence()) + logger.info(f"Set OpenDSS convergence to {dss.Solution.Convergence()}") logger.info('Max convergence error count {}.'.format(self._maxConvergenceErrorCount)) logger.info("initializing store") self.ResultContainer.InitializeDataStore(project.hdf_store, Steps, MC_scenario_number) @@ -388,12 +410,21 @@ def RunSimulation(self, project, scenario, MC_scenario_number=None): within_range = self._simulation_range.is_within_range(self._dssSolver.GetDateTime()) if within_range: pydss_has_converged = self.RunStep(step) + # logger.info(f'Finish simulation: step {step} 1') + # pydss_has_converged = self.RunStep(step) + # logger.info(f'Finish simulation: step {step} 2') + # pydss_has_converged = self.RunStep(step) + # logger.info(f'Finish simulation: step {step} 3') opendss_has_converged = dss.Solution.Converged() if not opendss_has_converged: - logger.error("OpenDSS did not converge at step=%s pydss_converged=%s", - step, pydss_has_converged) + logger.error(f"OpenDSS did not converge at step={step} pydss_converged={pydss_has_converged}") self._HandleOpenDSSConvergenceErrorChecks(step) - has_converged = pydss_has_converged and opendss_has_converged + + if pydss_has_converged: + has_converged = True + else: + has_converged = pydss_has_converged and opendss_has_converged + if step == 0 and self.ResultContainer is not None: size = make_human_readable_size(self.ResultContainer.max_num_bytes()) logger.info('Storage requirement estimation: %s, estimated based on first time step run.', size) @@ -401,6 +432,10 @@ def RunSimulation(self, project, scenario, MC_scenario_number=None): step, has_converged = self._RunPostProcessors(step, Steps, postprocessors) if self._increment_flag: step += 1 + # if step < 1: + # step += 1 + # else: + # step += 2 # In the case of a frequency sweep, the code updates results at each frequency. # Doing so again would cause a duplicate result. diff --git a/src/pydss/helics_interface.py b/src/pydss/helics_interface.py index 5a5faaa1..404f0077 100644 --- a/src/pydss/helics_interface.py +++ b/src/pydss/helics_interface.py @@ -4,7 +4,7 @@ import helics import os import re - +import pandas as pd from loguru import logger from pydss.simulation_input_models import SimulationSettingsModel @@ -208,6 +208,7 @@ def _create_helics_federate(self): IP = self._settings.helics.broker 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.broker: helics.helicsFederateInfoSetBroker(self.fedinfo, str(self._settings.helics.broker)) if self._settings.helics.broker_port: @@ -261,8 +262,9 @@ def updateHelicsSubscriptions(self): value = helics.helicsInputGetInteger(subscription.sub) if value and value != 0: - if value > 1e6 or value < -1e6: - value = 1.0 + logger.info(f"value is {value}") + if value > 1e6 or value < -1e6 or pd.isna(value): + value = 1.0 value = value * subscription.multiplier subscription.object.SetParameter(subscription.property, value) @@ -337,7 +339,39 @@ def updateHelicsPublications(self): def request_time_increment(self): error = sum([abs(sub.states[0] - sub.states[1]) for sub in self.subscriptions.subscriptions]) - r_seconds = self._dss_solver.GetTotalSeconds() #- self._dss_solver.GetStepResolutionSeconds() + # r_seconds = self._dss_solver.GetTotalSeconds() # parallel + r_seconds = self._dss_solver.GetTotalSeconds() + self._dss_solver.GetStepResolutionSeconds() # transmission priority + # if r_seconds > 0: + # r_seconds = self._dss_solver.GetTotalSeconds() + self._dss_solver.GetStepResolutionSeconds() + # logger.info(f"self._dss_solver.GetTotalSeconds(): {self._dss_solver.GetTotalSeconds()}, self._dss_solver.GetStepResolutionSeconds(): {self._dss_solver.GetStepResolutionSeconds()}") + # os.system('PAUSE') + if not self._settings.helics.iterative_mode: + while self.c_seconds < r_seconds: + self.c_seconds = helics.helicsFederateRequestTime(self._federate, r_seconds) + logger.info('Time requested: {} - time granted: {} '.format(r_seconds, self.c_seconds)) + return True, self.c_seconds + else: + + self.c_seconds, iteration_state = helics.helicsFederateRequestTimeIterative( + self._federate, + r_seconds, + helics.helics_iteration_request_iterate_if_needed + ) + + logger.info('Time requested: {} - time granted: {} error: {} it: {}'.format( + r_seconds, self.c_seconds, error, self.itr)) + if error > -1 and self.itr < self._co_convergance_max_iterations - 1: + self.itr += 1 + return False, self.c_seconds + else: + self.itr = 0 + return True, self.c_seconds + + def request_time_increment_2(self): + error = sum([abs(sub.states[0] - sub.states[1]) for sub in self.subscriptions.subscriptions]) + r_seconds = self._dss_solver.GetTotalSeconds() + self._dss_solver.GetStepResolutionSeconds()/2 + # logger.info('Time requested: {} - self._dss_solver.GetStepResolutionSeconds(): {} '.format(r_seconds, self._dss_solver.GetStepResolutionSeconds())) + # os.system('PAUSE') if not self._settings.helics.iterative_mode: while self.c_seconds < r_seconds: self.c_seconds = helics.helicsFederateRequestTime(self._federate, r_seconds) diff --git a/src/pydss/modes/solver_base.py b/src/pydss/modes/solver_base.py index df8168b1..95fca8ed 100644 --- a/src/pydss/modes/solver_base.py +++ b/src/pydss/modes/solver_base.py @@ -36,7 +36,7 @@ def __init__(self, dssInstance, settings: ProjectModel): #self._dssSolution.DblHour() self.reSolve() - logger.info("%s solver setup complete", settings.simulation_type) + logger.info(f"{settings.simulation_type} solver setup complete") def setFrequency(self, frequency): self._dssSolution.Frequency(frequency) diff --git a/src/pydss/pyControllers/Controllers/MotorStall.py b/src/pydss/pyControllers/Controllers/MotorStall.py index ba8c77db..2b768cec 100644 --- a/src/pydss/pyControllers/Controllers/MotorStall.py +++ b/src/pydss/pyControllers/Controllers/MotorStall.py @@ -3,6 +3,9 @@ import scipy.signal as signal import numpy as np import math +import os +from loguru import logger +import random from pydss.pyControllers.models import MotorStallSettings from pydss.pyControllers.pyControllerAbstract import ControllerAbstract @@ -19,34 +22,61 @@ def __init__(self, motor_obj, settings, dss_instance, elm_object_list, dss_solve self._controlled_element = motor_obj self._settings = MotorStallSettings(**settings) self._dss_solver = dss_solver - self.mode = 3 + self.mode = 1 self.model_mode = self._controlled_element.SetParameter('model', self.mode) # 3 - motor, 1 - Standard constant P+jQ + self.load_mult = dss_instance.Solution.LoadMult() self._controlled_element.SetParameter('vminpu', 0.0) self.kw_rated = self._controlled_element.GetParameter('kw') self.kvar_rated = self._controlled_element.GetParameter('kvar') - self.kva_rated = (self.kw_rated**2 + self.kvar_rated**2)**0.5 + self.rated_pf = self._settings.rated_pf + self.comp_lf = self._settings.comp_lf + # self.kva_rated = self.kw_rated / self.comp_lf + self.kva_rated = math.sqrt(self.kw_rated*self.kw_rated + self.kvar_rated*self.kvar_rated) self.stall_time_start = 0 self.stall = False + self.stall_counting = False + + self.rstrt_time_start = 0 + self.rstrt = True + self.rstrt_counting = False + + self.uv_time_start = 0 + self.uv_trip = False + self.uv_counting = False + self.f_uvr = self._settings.f_uvr + self.uv_tr1 = self._settings.uv_tr1 + self.t_tr1 = self._settings.t_tr1 + + self.vc_1off = self._settings.vc_1off + self.vc_2off = self._settings.vc_2off + self.vc_1on = self._settings.vc_1on + self.vc_2on = self._settings.vc_2on self.r_stall_pu = self._settings.r_stall_pu + self.x_stall_pu = self._settings.x_stall_pu + + self.va_base = 100 * 1e6 self.kvbase = self._controlled_element.sBus[0].GetVariable("kVBase") - self.i_base = 1e3 * self.kw_rated / 1e3 * self.kvbase + self.i_base = self.kva_rated / self.kvbase self.dt = dss_solver.GetStepResolutionSeconds() - self.h = signal.TransferFunction([1], [self._settings.t_th, 1]) - self.r = signal.TransferFunction([1], [self._settings.t_restart, 1]) - - self.u = [0, 0] - self.t_arr = [0, self.dt] - self.x = 0 - - self.i2r = 0 - - self.p_stall = 0 - self.q_stall = 0 + self.t_th = self._settings.t_th + self.i2r_rstr_prev = 1*1*self.r_stall_pu + self.temp_rstr_prev = self.i2r_rstr_prev + self.i2r_nonrstr_prev = 1*1*self.r_stall_pu + self.temp_nonrstr_prev = self.i2r_nonrstr_prev + + self.voltage_prev = 1.0 + self.i2r_rstr = 0 + self.temp_rstr = self.temp_rstr_prev + self.i2r_nonrstr = 0 + self.temp_nonrstr = self.temp_nonrstr_prev + self.trip_rstr = False + self.trip_nonrstr = False + self.thermal_threshold = self._settings.t_th1t + (self._settings.t_th2t-self._settings.t_th1t)*random.random() return @@ -61,74 +91,234 @@ def debugInfo(self): def Update(self, Priority, time, update_results): self.t = self._dss_solver.GetTotalSeconds() - if self.i_base: - self.current_pu = self._controlled_element.GetVariable('CurrentsMagAng')[0] / self.i_base + # logger.info(f"self.t: {self.t}") + # logger.info(f"self._controlled_element: {self._controlled_element}") + # logger.info(f"{self.name} - {self.kw_rated} - {self.kvar_rated} - {self.kvbase} - {self.i_base}") + if self.i_base: + self.p = self._controlled_element.GetVariable('Powers')[0] + self._controlled_element.GetVariable('Powers')[2] + self.q = self._controlled_element.GetVariable('Powers')[1] + self._controlled_element.GetVariable('Powers')[3] self.voltage = self._controlled_element.GetVariable('VoltagesMagAng')[0] - self.p = self._controlled_element.GetVariable('Powers')[0] - self.q = self._controlled_element.GetVariable('Powers')[1] self.voltage_pu = self._controlled_element.sBus[0].GetVariable("puVmagAngle")[0] + # logger.info(f"CurrentsMagAng: {self._controlled_element.GetVariable('CurrentsMagAng')}") + # logger.info(f"voltage_pu: {self.voltage_pu}") + # logger.info(f"Powers: {self._controlled_element.GetVariable('Powers')}") + # logger.info(f"self.kw_rated: {self.kw_rated}") + # logger.info(f"self.kvar_rated: {self.kvar_rated}") - if Priority == 0: - - i2r = self.current_pu ** 2 * self.r_stall_pu - self.i2r = max(self.i2r, i2r) - self.t = np.array([self.t_arr [-1], self.t]) - self.u = np.array([self.u[-1], self.i2r]) - tout, yout, xout = signal.lsim(self.h, self.u, self.t_arr, self.x) - self.x = xout[-1] - i2r_calc = yout[-1] / 50.0 - + comp_lf = self.comp_lf # self.p / self.kw_rated + comp_pf = self.rated_pf # self.p / (self.p**2 + self.q**2)**0.5 - comp_lf = self.p / self.kw_rated - comp_pf = 0.75 #self.p / (self.p**2 + self.q**2)**0.5 - - v_stall_adj = self._settings.v_stall*(1 + self._settings.lf_adj * (comp_lf-1)) - v_break_adj = self._settings.v_break*(1 + self._settings.lf_adj * (comp_lf-1)) + v_stall_adj = self._settings.v_stall*(1 + self._settings.lf_adj * (comp_lf-1)) + v_break_adj = self._settings.v_break*(1 + self._settings.lf_adj * (comp_lf-1)) + # logger.info(f"v_stall_adj: {v_stall_adj}") + # logger.info(f"v_break_adj: {v_break_adj}") + if Priority == 0: + ## stall and restart time clock + if self.voltage_pu < v_stall_adj and not self.stall: + if self.stall_counting: + self.stall_time = self._dss_solver.GetTotalSeconds() - self.stall_time_start + if self.stall_time > self._settings.t_stall and not self.stall: + self.stall = True + self.rstrt = False + else: + self.stall_time_start = self._dss_solver.GetTotalSeconds() + self.stall_counting = True + else: + self.stall_counting = False + if self.voltage_pu > self._settings.v_rstrt and not self.rstrt: + if self.rstrt_counting: + self.rstrt_time = self._dss_solver.GetTotalSeconds() - self.rstrt_time_start + if self.rstrt_time > self._settings.t_restart: + self.rstrt = True + else: + self.rstrt_time_start = self._dss_solver.GetTotalSeconds() + self.rstrt_counting = True + else: + self.rstrt_counting = False + ## uv trip + if self.voltage_pu < self.uv_tr1 and not self.uv_trip: + if self.uv_counting: + self.uv_time = self._dss_solver.GetTotalSeconds() - self.uv_time_start + if self.uv_time > self.t_tr1 and not self.uv_trip: + self.uv_trip = True + self.uv_counting = False + else: + self.uv_time_start = self._dss_solver.GetTotalSeconds() + self.uv_counting = True + if self.uv_trip: + Kthuv = 1.0 - self.f_uvr + else: + Kthuv = 1.0 + # logger.info(f"Kthuv: {Kthuv}") + + ## v contactor trip + if self.voltage_prev <= self.voltage_pu: + ## reconnect + if self.voltage_pu > self.vc_1on: + Kthc = 1.0 + elif self.voltage_pu < self.vc_2on: + Kthc = 0.0 + else: + Kthc = (self.voltage_pu - self.vc_2on)/(self.vc_1on - self.vc_2on) + else: + ## trip + if self.voltage_pu > self.vc_1off: + Kthc = 1.0 + elif self.voltage_pu < self.vc_2off: + Kthc = 0.0 + else: + Kthc = (self.voltage_pu - self.vc_2off)/(self.vc_1off - self.vc_2off) + # logger.info(f"Kthc: {Kthc}") + p0 = 1 - self._settings.k_p1 * (1-v_break_adj)**self._settings.n_p1 q0 = ((1 - comp_pf**2)**0.5 / comp_pf)-self._settings.k_q1*(1-v_break_adj)**self._settings.n_q1 - - p = self.p / self.kw_rated - q = self.q / self.kvar_rated - - if self.voltage_pu > v_break_adj and not self.stall: - p = p0 + self._settings.k_p1*(self.voltage_pu-v_break_adj)**self._settings.n_p1 - q = q0 + self._settings.k_q1*(self.voltage_pu-v_break_adj)**self._settings.n_q1 - self._controlled_element.SetParameter('kw', self.kw_rated * p) - self._controlled_element.SetParameter('kvar', self.kvar_rated * q) - - elif self.voltage_pu <= v_break_adj and not self.stall: - p = p0 + self._settings.k_p2 * (v_break_adj - self.voltage_pu)**self._settings.n_p2 - q = q0 + self._settings.k_q2 * (v_break_adj - self.voltage_pu)**self._settings.n_q2 - - self._controlled_element.SetParameter('kw', self.kw_rated * p ) - self._controlled_element.SetParameter('kvar', self.kvar_rated * q) + # logger.info(f"p0: {p0}") + # logger.info(f"q0: {q0}") + logger.info(f"self.voltage_pu: {self.voltage_pu}") - if self.voltage_pu < v_stall_adj and not self.stall: - self.p_stall = self._controlled_element.GetParameter('kw') - self.q_stall = self._controlled_element.GetParameter('kvar') - self.stall_time_start = self._dss_solver.GetTotalSeconds() - self.stall = True - - if self.voltage_pu > v_stall_adj and self.stall: - self.stall_time = self._dss_solver.GetTotalSeconds() - self.stall_time_start - if self.stall_time < self._settings.t_stall: - self._controlled_element.SetParameter('kw', self.p_stall) - self._controlled_element.SetParameter('kvar', self.q_stall) - else: - if i2r_calc < self._settings.t_th1t: - Kth = 1 - elif i2r_calc > self._settings.t_th2t: - Kth = 0 + # the operation model + if self.stall: + # stage III + p_stall = self.voltage_pu ** 2 * self.r_stall_pu / (self.r_stall_pu ** 2 + self.x_stall_pu ** 2) + q_stall = self.voltage_pu ** 2 * self.x_stall_pu / (self.r_stall_pu ** 2 + self.x_stall_pu ** 2) + # logger.info(f"p_stall: {p_stall}") + # logger.info(f"q_stall: {q_stall}") + if self.rstrt: + logger.info(f"Stage III: Motor stall and {self._settings.f_rst} of load is restarted") + if self.voltage_pu > v_break_adj: + p = p0 + self._settings.k_p1*(self.voltage_pu-v_break_adj)**self._settings.n_p1 + q = q0 + self._settings.k_q1*(self.voltage_pu-v_break_adj)**self._settings.n_q1 else: - m = 1 / (self._settings.t_th1t - self._settings.t_th2t) - c = - m * self._settings.t_th2t - Kth = m * i2r_calc + c + p = p0 + self._settings.k_p2 * (v_break_adj - self.voltage_pu)**self._settings.n_p2 + q = q0 + self._settings.k_q2 * (v_break_adj - self.voltage_pu)**self._settings.n_q2 + # self._controlled_element.SetParameter('kw', Kth * self.kw_rated * p_rstrt * self._settings.f_rst + Kth * p * self.kva_rated * (1-self._settings.f_rst)) + # self._controlled_element.SetParameter('kvar', Kth * self.kvar_rated * q_rstrt * self._settings.f_rst + Kth * q * self.kva_rated * (1-self._settings.f_rst)) + current_pu_nonrstrt = self.voltage_pu / math.sqrt(self.r_stall_pu ** 2 + self.x_stall_pu ** 2) + current_pu_rstrt = p / self.voltage_pu + self.i2r_rstr = current_pu_rstrt * current_pu_rstrt * self.r_stall_pu + self.i2r_nonrstr = current_pu_nonrstrt * current_pu_nonrstrt * self.r_stall_pu + self.temp_rstr = (self.dt*(self.i2r_rstr+self.i2r_rstr_prev)-(self.dt-2*self.t_th)*self.temp_rstr_prev)/(2*self.t_th+self.dt) + self.temp_nonrstr = (self.dt*(self.i2r_nonrstr+self.i2r_nonrstr_prev)-(self.dt-2*self.t_th)*self.temp_nonrstr_prev)/(2*self.t_th+self.dt) + p_rstrt = p * self._settings.f_rst + q_rstrt = q * self._settings.f_rst + p_nonrstrt = p_stall * (1 - self._settings.f_rst) + q_nonrstrt = q_stall * (1 - self._settings.f_rst) + # logger.info(f"p_rstrt: {p_rstrt}") + # logger.info(f"q_rstrt: {q_rstrt}") + # logger.info(f"p_nonrstrt: {p_nonrstrt}") + # logger.info(f"q_nonrstrt: {q_nonrstrt}") + # logger.info(f"current_pu_rstrt: {current_pu_rstrt}") + # logger.info(f"current_pu_nonrstrt: {current_pu_nonrstrt}") + # logger.info(f"self.i2r_nonrstr: {self.i2r_nonrstr}") + # logger.info(f"self.temp_nonrstr: {self.temp_nonrstr}") + else: + logger.info(f"Stage III: Motor stall and not restarted load") + # self._controlled_element.SetParameter('kw', Kth * p_stall * self.kva_rated) + # self._controlled_element.SetParameter('kvar', Kth * q_stall * self.kva_rated) + p_rstrt = p_stall * self._settings.f_rst + q_rstrt = q_stall * self._settings.f_rst + p_nonrstrt = p_stall * (1 - self._settings.f_rst) + q_nonrstrt = q_stall * (1 - self._settings.f_rst) + current_pu = self.voltage_pu / math.sqrt(self.r_stall_pu ** 2 + self.x_stall_pu ** 2) + self.i2r_rstr = current_pu * current_pu * self.r_stall_pu + self.i2r_nonrstr = current_pu * current_pu * self.r_stall_pu + self.temp_rstr = (self.dt*(self.i2r_rstr+self.i2r_rstr_prev)-(self.dt-2*self.t_th)*self.temp_rstr_prev)/(2*self.t_th+self.dt) + self.temp_nonrstr = (self.dt*(self.i2r_nonrstr+self.i2r_nonrstr_prev)-(self.dt-2*self.t_th)*self.temp_nonrstr_prev)/(2*self.t_th+self.dt) + # logger.info(f"p_rstrt: {p_rstrt}") + # logger.info(f"q_rstrt: {q_rstrt}") + # logger.info(f"p_nonrstrt: {p_nonrstrt}") + # logger.info(f"q_nonrstrt: {q_nonrstrt}") + # logger.info(f"current_pu: {current_pu}") + # logger.info(f"self.i2r_rstr: {self.i2r_rstr}") + # logger.info(f"self.temp_rstr: {self.temp_rstr}") + else: + if self.voltage_pu > v_break_adj: + # stage I + logger.info(f"Stage I: normal operation") + p = p0 + self._settings.k_p1*(self.voltage_pu-v_break_adj)**self._settings.n_p1 + q = q0 + self._settings.k_q1*(self.voltage_pu-v_break_adj)**self._settings.n_q1 + # logger.info(f"p1: {p}") + # logger.info(f"q1: {q}") + else: + # stage II or before stall + logger.info(f"Stage II: Motor voltage below the break down voltage") + p = p0 + self._settings.k_p2 * (v_break_adj - self.voltage_pu)**self._settings.n_p2 + q = q0 + self._settings.k_q2 * (v_break_adj - self.voltage_pu)**self._settings.n_q2 + # logger.info(f"p2: {p}") + # logger.info(f"q2: {q}") + current_pu = p / self.voltage_pu + p_rstrt = p * self._settings.f_rst + p_nonrstrt = p * (1 - self._settings.f_rst) + q_rstrt = q * self._settings.f_rst + q_nonrstrt = q * (1 - self._settings.f_rst) + self.i2r_rstr = current_pu * current_pu * self.r_stall_pu + self.i2r_nonrstr = current_pu * current_pu * self.r_stall_pu + self.temp_rstr = (self.dt*(self.i2r_rstr+self.i2r_rstr_prev)-(self.dt-2*self.t_th)*self.temp_rstr_prev)/(2*self.t_th+self.dt) + self.temp_nonrstr = (self.dt*(self.i2r_nonrstr+self.i2r_nonrstr_prev)-(self.dt-2*self.t_th)*self.temp_nonrstr_prev)/(2*self.t_th+self.dt) + # logger.info(f"p_rstrt: {p_rstrt}") + # logger.info(f"q_rstrt: {q_rstrt}") + # logger.info(f"p_nonrstrt: {p_nonrstrt}") + # logger.info(f"q_nonrstrt: {q_nonrstrt}") + # logger.info(f"current_pu: {current_pu}") + # logger.info(f"self.i2r_rstr: {self.i2r_rstr}") + # logger.info(f"self.temp_rstr: {self.temp_rstr}") - self._controlled_element.SetParameter('kw', self.p_stall * Kth ) - self._controlled_element.SetParameter('kvar', self.q_stall * Kth ) - - self.model_mode_old = self.model_mode - + # thermal protection + if self.stall: + if self.trip_rstr: + # logger.info(f"Motor is tripped") + Kth_rstr = 0 + ## for single bus system + elif self.temp_rstr > self._settings.t_th2t: + self.trip_rstr = True + Kth_rstr = 0 + elif self.temp_rstr > self._settings.t_th1t and self.temp_rstr <= self._settings.t_th2t: + Kth_rstr = 1 - (self.temp_rstr - self._settings.t_th1t)/(self._settings.t_th2t - self._settings.t_th1t) + # ## for multiple bus system + # elif self.temp_rstr > self.thermal_threshold: + # self.trip_rstr = True + # Kth_rstr = 0 + else: + Kth_rstr = 1 + + if self.trip_nonrstr: + # logger.info(f"Motor is tripped") + Kth_nonrstr = 0 + ## for single bus system + elif self.temp_nonrstr > self._settings.t_th2t: + self.trip_nonrstr = True + Kth_nonrstr = 0 + elif self.temp_nonrstr > self._settings.t_th1t and self.temp_nonrstr <= self._settings.t_th2t: + Kth_nonrstr = 1 - (self.temp_nonrstr - self._settings.t_th1t)/(self._settings.t_th2t - self._settings.t_th1t) + # # for multiple bus system + # elif self.temp_nonrstr > self.thermal_threshold: + # self.trip_nonrstr = True + # Kth_nonrstr = 0 + else: + Kth_nonrstr = 1 + else: + Kth_rstr = 1 + Kth_nonrstr = 1 + + # logger.info(f"Kth_rstr: {Kth_rstr}") + # logger.info(f"Kth_nonrstr: {Kth_nonrstr}") + + pset = (Kth_rstr*p_rstrt + Kth_nonrstr*p_nonrstrt) * self.kw_rated + qset = (Kth_rstr*q_rstrt + Kth_nonrstr*q_nonrstrt) * self.kw_rated + if self.stall: + qset = (Kth_rstr*q_rstrt + Kth_nonrstr*q_nonrstrt) * self.kw_rated + # logger.info(f"pset: {pset}") + # logger.info(f"qset: {qset}") + + self._controlled_element.SetParameter('kw', Kthc * Kthuv * pset ) + self._controlled_element.SetParameter('kvar', Kthc * Kthuv * qset ) + # os.system("PAUSE") + + self.voltage_prev = self.voltage_pu + self.temp_rstr_prev = self.temp_rstr + self.i2r_rstr_prev = self.i2r_rstr + self.temp_nonrstr_prev = self.temp_nonrstr + self.i2r_nonrstr_prev = self.i2r_nonrstr return 0 diff --git a/src/pydss/pyControllers/Controllers/PvFrequencyRideThru.py b/src/pydss/pyControllers/Controllers/PvFrequencyRideThru.py index 7ad3a866..abf04ebf 100644 --- a/src/pydss/pyControllers/Controllers/PvFrequencyRideThru.py +++ b/src/pydss/pyControllers/Controllers/PvFrequencyRideThru.py @@ -1,7 +1,7 @@ from ipaddress import v4_int_to_packed from pydss.pyControllers.pyControllerAbstract import ControllerAbstract from shapely.geometry import MultiPoint, Polygon, Point, MultiPolygon -from shapely.ops import triangulate, cascaded_union +from shapely.ops import triangulate, unary_union import matplotlib.pyplot as plt import datetime import math @@ -240,11 +240,11 @@ def __CreateOperationRegions(self): MayTripRegion = TotalRegion.difference(intersection) #everything not in the active region is the white "may trip" if self.__Settings['May trip operation'] == 'Trip': - self.CurrLimRegion = cascaded_union([MandatoryRegion1, MandatoryRegion2]) #Naming probably not appropriate with frequency variation. - self.TripRegion = cascaded_union([OFtripRegion, UFtripRegion, MayTripRegion]) + self.CurrLimRegion = unary_union([MandatoryRegion1, MandatoryRegion2]) #Naming probably not appropriate with frequency variation. + self.TripRegion = unary_union([OFtripRegion, UFtripRegion, MayTripRegion]) elif self.__Settings['May trip operation'] == 'Ride-Through': - self.CurrLimRegion = cascaded_union([MandatoryRegion1, MandatoryRegion2, MayTripRegion]) - self.TripRegion = cascaded_union([OFtripRegion, UFtripRegion]) + self.CurrLimRegion = unary_union([MandatoryRegion1, MandatoryRegion2, MayTripRegion]) + self.TripRegion = unary_union([OFtripRegion, UFtripRegion]) else: assert False self.ContinuousRegion = ContinuousRegion diff --git a/src/pydss/pyControllers/Controllers/PvVoltageRideThru.py b/src/pydss/pyControllers/Controllers/PvVoltageRideThru.py index ca285a6d..f58e0d7c 100644 --- a/src/pydss/pyControllers/Controllers/PvVoltageRideThru.py +++ b/src/pydss/pyControllers/Controllers/PvVoltageRideThru.py @@ -2,6 +2,9 @@ from shapely.ops import triangulate, unary_union import datetime import math +import random +from loguru import logger +import os from pydss.pyControllers.pyControllerAbstract import ControllerAbstract from pydss.pyControllers.models import PvVoltageRideThruModel @@ -51,13 +54,39 @@ def __init__(self, pv_object, settings, dss_instance, elm_object_list, dss_solve # Initializing the model #pv_object.SetParameter('kvar', 0) #pv_object.SetParameter('kva', self.model.kva) + self.dt = dss_solver.GetStepResolutionSeconds() self._p_rated = float(pv_object.GetParameter('kW')) + self._pf_rated = float(pv_object.GetParameter('pf')) + self._phase_rated = float(pv_object.GetParameter('Phases')) + + logger.info(f"{self._name} -> _p_rated: {self._p_rated }, _pf_rated: {self._pf_rated}, _phase_rated: {self._phase_rated}, dt: {self.dt}") + # os.system("PAUSE") + + self.t_rv = 0.02 + self.t_v = 0.02 + self.t_g = 0.02 + self.rrpwr = 2 + self.Imax = 1.2 + self.ul0 = 0.44 + self.ul1 = 0.49 + self.uh0 = 1.2 + self.uh1 = 1.15 + self.ul = self.ul0 + (self.ul1-self.ul0)*random.random() + self.uh = self.uh1 + (self.uh0-self.uh1)*random.random() + self.u_in_prev = 1.0 + self.ut_filt_prev = 1.0 + self.Vtrip_ctrl_prev = 1.0 + self.Vmult_prev = 1.0 + self.Ip_out_prev = 1.0 + self.Iq_out_prev = 0.0 + self.Ip_out_filt_prev = 1.0 + self.Iq_out_filt_prev = 0.0 # MISC settings self._trip_deadtime_sec = self.model.reconnect_deadtime_sec self._time_to_p_max_sec = self.model.reconnect_pmax_time_sec #self._p_rated = self.model.max_kw - self._voltage_calc_mode = self.model.voltage_calc_mode + self._voltage_calc_mode = self.model.voltage_calc_mode # max # initialize deadtimes and other variables self._initialize_ride_through_settings() if self.model.follow_standard == PvStandard.IEEE_1547_2003: @@ -177,7 +206,7 @@ def _create_operation_regions(self): self._fault_counter_max = 2 self._fault_counter_clearing_time_sec = 10 - elif self.model.ride_through_category == RideThroughCategory.CATEGORY_I: + elif self.model.ride_through_category == RideThroughCategory.CATEGORY_III: ov2pu_eq = 1.20 ov2sec_eq = 0.16 ov1pu_min = 1.1 @@ -200,41 +229,48 @@ def _create_operation_regions(self): #check overvoltage points if self.model.ov_2_pu != ov2pu_eq: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False if self.model.ov_2_ct_sec != ov2sec_eq: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False if self.model.ov_1_pu < ov1pu_min and self.model.ov_1_pu > ov1pu_max: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False if self.model.ov_1_ct_sec < ov1sec_min and self.model.ov_1_ct_sec > ov1sec_max: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False #check undervoltage points if self.model.uv_2_pu < uv2pu_min and self.model.uv_2_pu > uv2pu_max: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False if self.model.uv_2_ct_sec < uv2sec_min and self.model.uv_2_ct_sec > uv2sec_max: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False if self.model.uv_1_pu < uv1pu_min and self.model.uv_1_pu > uv1pu_max: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False if self.model.uv_1_ct_sec uv1sec_max: - #print("User defined setting outside of IEEE 1547 acceptable range.") + logger.error("User defined setting outside of IEEE 1547 acceptable range.") assert False - self._controlled_element.SetParameter('Model', '7') + # Re-set kV to force RecalcElementData before switching from Model 7 to 1. + # Model 7 internally uses VBase = kV*1000/sqrt(3) even for 1-phase generators, + # but Model 1 expects VBase = kV*1000 for 1-phase. Changing only the 'Model' + # property does not trigger RecalcElementData, leaving a stale VBase that + # causes Model 1 to operate in constant-impedance mode (inflated output). + kv = self._controlled_element.GetParameter('kV') + self._controlled_element.SetParameter('kV', kv) + self._controlled_element.SetParameter('Model', '1') # change to model 1 self._controlled_element.SetParameter('Vmaxpu', V[0]) - self._controlled_element.SetParameter('Vminpu', V[1]) + self._controlled_element.SetParameter('Vminpu', '0.1')# V[1]) contineous_points = [Point(V[0], 0), Point(V[0], t_max), Point(V[1], t_max), Point(V[1], 0)] contineous_region = Polygon([[p.y, p.x] for p in contineous_points]) @@ -260,10 +296,9 @@ def _create_operation_regions(self): may_trip_region = total_region.difference(intersection) if self.model.ride_through_category in [RideThroughCategory.CATEGORY_I, RideThroughCategory.CATEGORY_II]: - if self.model.permissive_operation == PermissiveOperation.CURRENT_LIMITED: + if self.model.permissive_operation == PermissiveOperation.CURRENT_LIMITED: if self.model.may_trip_operation == MayTripOperation.PERMISSIVE_OPERATION: - self.curr_lim_region = unary_union( - [permissive_ov_region, permissive_uv_region, mandatory_region, may_trip_region]) + self.curr_lim_region = unary_union([permissive_ov_region, permissive_uv_region, mandatory_region, may_trip_region]) self.momentary_sucession_region = None self.trip_region = unary_union([ov_trip_region, uv_trip_region]) else: @@ -299,11 +334,15 @@ def Update(self, priority, time, update_results): error = 0 self.time_change = self.time != (priority, time) self.time = time + logger.info(f"self._name: {self._name}") + logger.info(f"priority: {priority}") if priority == 0: self._is_connected = self._connect() + logger.info(f"self._is_connected: {self._is_connected}") if priority == 2: u_in = self._update_violaton_timers() + logger.info(f"Update u_in: {u_in}") if self.model.follow_standard == PvStandard.IEEE_1547_2018: self.voltage_ride_through(u_in) elif self.model.follow_standard == PvStandard.IEEE_1547_2003: @@ -311,7 +350,7 @@ def Update(self, priority, time, update_results): else: raise Exception("Valid standard setting defined. Options are: 1547-2003, 1547-2018") - P = -sum(self._controlled_element.GetVariable('Powers')[::2]) + # P = -sum(self._controlled_element.GetVariable('Powers')[::2]) # self.power_hist.append(P) # self.voltage_hist.append(u_in) # self.timer_hist.append(self._u_violation_time) @@ -330,39 +369,51 @@ def trip(self, u_in): def voltage_ride_through(self, u_in): """ Implementation of the IEEE1587-2018 voltage ride-through requirements for inverter systems """ - self._fault_counter_clearing_time_sec = 1 + # self._fault_counter_clearing_time_sec = 1 Pm = Point(self._u_violation_time, u_in) + logger.info(Pm) if Pm.within(self.curr_lim_region): region = 0 is_in_contioeous_region = False + logger.info(f"curr_lim_region") elif self.momentary_sucession_region and Pm.within(self.momentary_sucession_region): region = 1 is_in_contioeous_region = False - self._trip(self.__dss_solver.GetStepSizeSec(), 0.4, False) + self._trip(self.__dss_solver.GetStepSizeSec(), 0.5, False) # rrpt = 2.0 + logger.info(f"momentary_sucession_region") elif Pm.within(self.trip_region): region = 2 is_in_contioeous_region = False if self.region == [3, 1, 1]: + logger.info(f"self.region 3:1:1 {self.region}") self._trip(self._trip_deadtime_sec, self._time_to_p_max_sec, False, True) - else: - self._trip(self._trip_deadtime_sec, self._time_to_p_max_sec, False) + else: + logger.info(f"self.region else {self.region}") + self._trip(self._trip_deadtime_sec, self._time_to_p_max_sec, True) + logger.info(f"trip_region") else: is_in_contioeous_region = True region = 3 - - self.region = self.region[1:] + self.region[:1] + logger.info(f"is_in_contioeous_region") + + logger.info(f"old self.region: {self.region}") + self.region = self.region[1:] + self.region[:1] # shift [1,2,3]->[2,3,1] self.region[0] = region + logger.info(f"new self.region: {self.region}") if is_in_contioeous_region and not self._is_in_contioeous_region: + logger.info(f"faulr just clear") self._fault_window_clearing_start_time = self.__dss_solver.GetDateTime() clearing_time = (self.__dss_solver.GetDateTime() - self._fault_window_clearing_start_time).total_seconds() + logger.info(f"clearing_time: {clearing_time}") if self._is_in_contioeous_region and not is_in_contioeous_region: if clearing_time <= self._fault_counter_clearing_time_sec: self._fault_counter += 1 if self._fault_counter > self._fault_counter_max: if self.model.multiple_disturdances == MultipleDisturbances.TRIP: + logger.info(f"multiple_disturdances trip") self._trip(self._trip_deadtime_sec, self._time_to_p_max_sec, True) self._fault_counter = 0 else: @@ -371,7 +422,104 @@ def voltage_ride_through(self, u_in): self._fault_counter = 0 self._is_in_contioeous_region = is_in_contioeous_region return + + def DERA_control(self, mode): + u_in = self._controlled_element.GetVariable('VoltagesMagAng')[::2] + u_base = self._controlled_element.sBus[0].GetVariable('kVBase') * 1000 + logger.info(f"{self._name} -> u_in: {u_in}") + Pord = 1.0 + if self._phase_rated == 1: + # single-phase DER + u_in = max(u_in) / u_base + else: + # three-phase DER + u_in = sum(u_in) / u_base / 3.0 + + if self.__dss_solver.GetTotalSeconds() < round(self.dt, 5): + logger.info(f"OpenDSS Initializatin -> time: {self.__dss_solver.GetTotalSeconds()}") + self.u_in_prev = u_in + self.ut_filt_prev = u_in + self.Vtrip_ctrl_prev = 1.0 + self.Vmult_prev = 1.0 + self.Ip_out_prev = 1 / u_in + self.Iq_out_prev = 0.0 + self.Ip_out_filt_prev = 1 / u_in + self.Iq_out_filt_prev = 0.0 + + self.ut_filt = (self.dt*(u_in+self.u_in_prev)-(self.dt-2*self.t_rv)*self.ut_filt_prev)/(2*self.t_rv+self.dt) + Ip = Pord / self.ut_filt + Iq = Pord * math.tan(math.acos(self._pf_rated)) / self.ut_filt + + # Q prority current limit control + if Iq > self.Imax: + Iqcmd = self.Imax + elif Iq < -self.Imax: + Iqcmd = -self.Imax + else: + Iqcmd = Iq + Ipmax = math.sqrt(self.Imax*self.Imax - Iqcmd*Iqcmd) + if Ip > Ipmax: + Ipcmd = Ipmax + elif Ip < 0: + Ipcmd = 0 + else: + Ipcmd = Ip + + # undervoltage multiplier + if self.ut_filt < self.ul0: + v11 = 0 + elif self.ut_filt > self.ul1: + v11 = 1 + else: + v11 = (self.ut_filt - self.ul0)/(self.ul1 - self.ul0) + # overoltage multiplier + if self.ut_filt < self.uh1: + v12 = 1 + elif self.ut_filt > self.uh0: + v12 = 0 + else: + v12 = (self.ut_filt - self.uh1)/(self.uh0 - self.uh1) + Vtrip_ctrl = v11*v12 + self.Vmult = (self.dt*(Vtrip_ctrl+self.Vtrip_ctrl_prev)-(self.dt-2*self.t_v)*self.Vmult_prev)/(2*self.t_v+self.dt) + + Ip_out = self.Vmult*Ipcmd*mode + Iq_out = self.Vmult*Iqcmd*mode + if Ip_out < self.Ip_out_prev: + Ip_out_filt = (self.dt*(Ip_out+self.Ip_out_prev)-(self.dt-2*0.01)*self.Ip_out_filt_prev)/(2*0.01+self.dt) + else: + Ip_out_filt = (self.dt*(Ip_out+self.Ip_out_prev)-(self.dt-2*self.t_g)*self.Ip_out_filt_prev)/(2*self.t_g+self.dt) + Iq_out_filt = (self.dt*(Iq_out+self.Iq_out_prev)-(self.dt-2*self.t_g)*self.Iq_out_filt_prev)/(2*self.t_g+self.dt) + + if self._fault_counter > 0 and (Ip_out_filt > self.Ip_out_filt_prev + self.rrpwr * self.dt): + Ip_out_filt = self.Ip_out_filt_prev + self.rrpwr * self.dt + + logger.info(f"u_in: {u_in}") + logger.info(f"Ipcmd: {Ipcmd}") + logger.info(f"Ip_out: {Ip_out}") + logger.info(f"self.ut_filt: {self.ut_filt}") + logger.info(f"Ip_out_filt: {Ip_out_filt}") + logger.info(f"Ip_out_filt_prev: {self.Ip_out_filt_prev}") + logger.info(f"Iq_out_filt: {Iq_out_filt}") + logger.info(f"self.Vmult: {self.Vmult}") + logger.info(f"mode: {mode}") + + self.DERA_p_limit = math.sqrt(Ip_out_filt*Ip_out_filt + Iq_out_filt*Iq_out_filt) * u_in * self._p_rated + + self.ut_filt_prev = self.ut_filt + self.u_in_prev = u_in + self.Vmult_prev = self.Vmult + self.Vtrip_ctrl_prev = Vtrip_ctrl + self.Ip_out_prev = Ip_out + self.Iq_out_prev = Iq_out + self.Ip_out_filt_prev = Ip_out_filt + self.Iq_out_filt_prev = Iq_out_filt + + + def TODO(self): + #add timer for DERA undervoltage/overvoltage voltage multiplier block + return None + def _connect(self): if not self._is_connected: u_in = self._controlled_element.GetVariable('VoltagesMagAng')[::2] @@ -382,6 +530,12 @@ def _connect(self): self.voltage[0] = u_in u_in = sum(self.voltage) / len(self.voltage) deadtime = (self.__dss_solver.GetDateTime() - self._tripped_start_time).total_seconds() + logger.info(f"_connect u_in: {u_in}") + logger.info(f"deadtime: {deadtime}") + logger.info(f"self._tripped_dead_time: {self._tripped_dead_time}") + # self.DERA_control(0.0) + # logger.info(f"self.DERA_p_limit: {self.DERA_p_limit}") + # os.system("PAUSE") if u_in < self._rvs[0] and u_in > self._rvs[1] and deadtime >= self._tripped_dead_time: self._controlled_element.SetParameter('enabled', True) @@ -392,6 +546,10 @@ def _connect(self): conntime = (self.__dss_solver.GetDateTime() - self._reconnect_start_time).total_seconds() self._p_limit = conntime / self._tripped_p_max_delay * self._p_rated if conntime < self._tripped_p_max_delay \ else self._p_rated + logger.info(f"self._p_limit: {self._p_limit}") + # self.DERA_control(1.0) + # logger.info(f"self.DERA_p_limit: {self.DERA_p_limit}") + # os.system("PAUSE") self._controlled_element.SetParameter('kw', self._p_limit) return self._is_connected @@ -416,12 +574,14 @@ def _trip(self, Deadtime, time2Pmax, forceTrip, permissive_to_trip=False): self._tripped_start_time = self.__dss_solver.GetDateTime() self._tripped_p_max_delay = time2Pmax self._tripped_dead_time = Deadtime + logger.info(f"self._tripped_dead_time: {self._tripped_dead_time}") return def _update_violaton_timers(self): u_in = self._controlled_element.GetVariable('VoltagesMagAng')[::2] u_base = self._controlled_element.sBus[0].GetVariable('kVBase') * 1000 u_in = max(u_in) / u_base if self._voltage_calc_mode == VoltageCalcModes.MAX else sum(u_in) / (u_base * len(u_in)) + logger.info(f"{self._name}: voltage {u_in}") if self.use_avg_voltage: self.voltage = self.voltage[1:] + self.voltage[:1] self.voltage[0] = u_in diff --git a/src/pydss/pyControllers/models.py b/src/pydss/pyControllers/models.py index cae7cb20..d8dc3446 100644 --- a/src/pydss/pyControllers/models.py +++ b/src/pydss/pyControllers/models.py @@ -165,6 +165,38 @@ class MotorStallSimpleSettings(BaseControllerModel): class MotorStallSettings(BaseControllerModel): + t_stall: Annotated[ + float, + Field(0.032, description="Stall time (range) based on laboratory testing"), + ] + t_restart: Annotated[ + float, + Field(0.300, description="Induction motor restart time is relatively short"), + ] + comp_lf: Annotated[ + float, + Field(1.000, description="Motor D real power to motor base ratio"), + ] + rated_pf: Annotated[ + float, + Field(0.980, description="Assumed slightly inductive motors load"), + ] + v_stall: Annotated[ + float, + Field(0.45, ge=0.40, le=0.70, description="Stall voltage (range) based on laboratory testing"), + ] + r_stall_pu: Annotated[ + float, + Field(0.120, description="Based on laboratory testing results of residential air-conditioners."), + ] + x_stall_pu: Annotated[ + float, + Field(0.120, description="Based on laboratory testing results of residential air-conditioners."), + ] + lf_adj: Annotated[ + float, + Field(0.0, description="Load factor adjustment to the stall voltage10"), + ] k_p1: Annotated[ float, Field(0 , description="Real power constant for running state 111"), @@ -197,66 +229,60 @@ class MotorStallSettings(BaseControllerModel): float, Field(2.5, description="Reactive power exponent for running state 2."), ] - t_th: Annotated[ + v_break: Annotated[ float, - Field(4.0, description="Varies based on manufacturer and external factors - sensitivity analysis required"), + Field(0.86, description="Compressor motor 'breakdown' voltage (pu)"), ] f_rst: Annotated[ float, Field(0.2, description="Captures diversity in load; also based on testing (fraction of motors capable of restart)."), ] - lf_adj: Annotated[ + v_rstrt: Annotated[ float, - Field(0.0, description="Load factor adjustment to the stall voltage10"), + Field(0.95, description="Reconnect when acceptable voltage met"), ] - t_th1t: Annotated[ + vc_1off: Annotated[ float, - Field(0.7, description="Assumed tripping starting at 70% temperature"), + Field(0.5, description="Motor D voltage at which contactors start opening, pu"), ] - t_th2t: Annotated[ + vc_2off: Annotated[ float, - Field(1.9, description="Assumed all tripped at 190% temperature"), - ] - p_fault: Annotated[ + Field(0.4, description="Motor D voltage at which all contactors are opened, pu"), + ] + vc_1on: Annotated[ float, - Field(3.5, ge=3.0, le=5.0, description="Active power multiplier post fault."), + Field(0.6, description="Motor D voltage at which all contactors are reclosed, pu"), ] - q_fault: Annotated[ + vc_2on: Annotated[ float, - Field(5.0, ge=3.0, le=7.0, description="Reactive power multiplier post fault."), + Field(0.5, description="Motor D voltage at which contactors start reclosing, pu"), ] - v_stall: Annotated[ + t_th: Annotated[ float, - Field(0.55, ge=0.45, le=0.60, description="Stall voltage (range) based on laboratory testing"), + Field(15.0, description="Varies based on manufacturer and external factors - sensitivity analysis required"), ] - v_break: Annotated[ + t_th1t: Annotated[ float, - Field(0.86, description="Compressor motor 'breakdown' voltage (pu)"), + Field(0.7, description="Assumed tripping starting at 70% temperature"), ] - v_rstrt: Annotated[ + t_th2t: Annotated[ float, - Field(0.95, description="Reconnect when acceptable voltage met"), + Field(1.9, description="Assumed all tripped at 190% temperature"), ] - t_stall: Annotated[ + f_uvr: Annotated[ float, - Field(0.032, description="Stall time (range) based on laboratory testing"), + Field(0.0, ge=0.0, le=0.10, description="Motor D heating time constant, s"), ] - t_restart: Annotated[ + uv_tr1: Annotated[ float, - Field(0.300, description="Induction motor restart time is relatively short"), + Field(0.6, ge=0.45, le=0.70, description="Motor D 1st undervoltage pick-up, pu"), ] - rated_pf: Annotated[ + t_tr1: Annotated[ float, - Field(0.939, description="Assumed slightly inductive motors load"), + Field(0.02, ge=0.01, le=0.05, description="Motor D 1st undervoltage trip delay, s"), ] - r_stall_pu: Annotated[ - float, - Field(0.100, description="Based on laboratory testing results of residential air-conditioners."), - ] - x_stall_pu: Annotated[ - float, - Field(0.100, description="Based on laboratory testing results of residential air-conditioners."), - ] + + class PvControllerModel(BaseControllerModel): diff --git a/src/pydss/pyControllers/pyController.py b/src/pydss/pyControllers/pyController.py index 21d22351..b54b8855 100644 --- a/src/pydss/pyControllers/pyController.py +++ b/src/pydss/pyControllers/pyController.py @@ -19,7 +19,7 @@ def Create(ElmName, ControllerType, Settings, ElmObjectList, dssInstance, dssSol "Please define the controller in ~pydss\pyControllers\Controllers".format( ControllerType ) - + assert (ElmName in ElmObjectList), "'{}' does not exist in the pydss master object dictionary.".format(ElmName) relObject = ElmObjectList[ElmName] diff --git a/src/pydss/pyDSS.py b/src/pydss/pyDSS.py index 9108de24..079f7308 100644 --- a/src/pydss/pyDSS.py +++ b/src/pydss/pyDSS.py @@ -43,7 +43,7 @@ def run_scenario(self, project, scenario, settings: SimulationSettingsModel, dry opendss = OpenDSS(settings) self._dump_scenario_simulation_settings(settings) - logger.info('Running scenario: %s', settings.project.active_scenario) + logger.info(f'Running scenario: {settings.project.active_scenario}') if settings.monte_carlo.num_scenarios > 0: opendss.RunMCsimulation(project, scenario, samples=settings.monte_carlo.num_scenarios) else: diff --git a/src/pydss/pydss_project.py b/src/pydss/pydss_project.py index 26d29dda..4ded8cac 100644 --- a/src/pydss/pydss_project.py +++ b/src/pydss/pydss_project.py @@ -342,6 +342,7 @@ def run(self, logging_configured=True, tar_project=False, zip_project=False, dry else: store_filename = os.path.join(self._project_dir, STORE_FILENAME) self._dump_simulation_settings() + logger.info("indicator 3") driver = None if self._settings.exports.export_data_in_memory: diff --git a/src/pydss/registry.py b/src/pydss/registry.py index de00947c..522d35ad 100644 --- a/src/pydss/registry.py +++ b/src/pydss/registry.py @@ -39,50 +39,13 @@ ), }, ], - ControllerType.PV_VOLTAGE_RIDETHROUGH.value: [ - { - "name": "NO_VRT_DVS_test", - "filename": os.path.join( - os.path.dirname(getattr(pydss, "__path__")[0]), - "pydss/pyControllers/Controllers/Settings/VoltageRideThru.toml" - ), - }, - { - "name": "1547_CAT_III_test", - "filename": os.path.join( - os.path.dirname(getattr(pydss, "__path__")[0]), - "pydss/pyControllers/Controllers/Settings/VoltageRideThru.toml" - ), - }, - { - "name": "1547_CAT_III_DVS_test", - "filename": os.path.join( - os.path.dirname(getattr(pydss, "__path__")[0]), - "pydss/pyControllers/Controllers/Settings/VoltageRideThru.toml" - ), - }, - ], + ControllerType.PV_VOLTAGE_RIDETHROUGH.value: [], ControllerType.SOCKET_CONTROLLER.value: [], ControllerType.STORAGE_CONTROLLER.value: [], ControllerType.XMFR_CONTROLLER.value: [], ControllerType.MOTOR_STALL.value: [], ControllerType.MOTOR_STALL_SIMPLE.value: [], ControllerType.FAULT_CONTROLLER.value: [], - ControllerType.DYNAMIC_VOLTAGE_SUPPORT.value:[{ - "name": "DVS_test", - "filename": os.path.join( - os.path.dirname(getattr(pydss, "__path__")[0]), - "PyDSS/pyControllers/Controllers/Settings/DynamicVoltageSupport.toml" - ), - }, - { - "name": "DVS_VRT_test", - "filename": os.path.join( - os.path.dirname(getattr(pydss, "__path__")[0]), - "PyDSS/pyControllers/Controllers/Settings/DynamicVoltageSupport.toml" - ), - }, - ], }, } diff --git a/src/pydss/simulation_input_models.py b/src/pydss/simulation_input_models.py index 8b2911ce..ca10d6f9 100644 --- a/src/pydss/simulation_input_models.py +++ b/src/pydss/simulation_input_models.py @@ -694,6 +694,31 @@ def check_max_co_iterations(cls, val): raise ValueError(f"max_co_iterations must be between 1 and 1000: {val}") return val +class ThreatModel(InputsBaseModel): + threat_federate: Annotated[ + bool, + Field( + title="threat_federate", + description="Set to true to enable the Threat Federate.", + default=False, + alias="Enable Threat Federate", + )] + federate_name: Annotated[ + str, + Field( + title="federate_name", + description="Name of the federate", + default="threat_test", + alias="Federate Name", + )] + threat_mapping: Annotated[ + Optional[Path], + Field( + None, + title="threat_mapping", + description="Path of the mapping file for threat federate.", + alias="Threat Mapping Path", + )] class LoggingModel(InputsBaseModel): """Defines the user inputs for controlling logging.""" @@ -1131,6 +1156,14 @@ class SimulationSettingsModel(InputsBaseModel): default=HelicsModel(), alias="Helics", )] + threat_settings:Annotated[ + ThreatModel, + Field( + title="threat_settings", + description="Threat settings", + default=ThreatModel(), + alias="Threat", + )] logging: Annotated[ LoggingModel, Field( diff --git a/src/pydss/value_storage.py b/src/pydss/value_storage.py index 1a078ae6..4eb5e498 100644 --- a/src/pydss/value_storage.py +++ b/src/pydss/value_storage.py @@ -355,7 +355,7 @@ def set_nan(self): if np.issubdtype(self._value_type, np.int64): self._value = INTEGER_NAN else: - self._value = np.NaN + self._value = np.nan def set_value(self, value): self._value = value @@ -487,7 +487,7 @@ def set_name(self, name): def set_nan(self): for i in range(len(self._value)): - self._value[i] = np.NaN + self._value[i] = np.nan def set_value(self, value): self._value = value