diff --git a/LocalFeeder/FeederSimulator.py b/LocalFeeder/FeederSimulator.py index f98ecd1..99caa63 100644 --- a/LocalFeeder/FeederSimulator.py +++ b/LocalFeeder/FeederSimulator.py @@ -1,5 +1,6 @@ """Core class to abstract OpenDSS into Feeder class.""" +import csv import json import logging import math @@ -8,7 +9,7 @@ import time from enum import Enum from time import strptime -from typing import Dict, List, Optional, Set +from typing import Dict, List, Optional, Set, Tuple import boto3 import numpy as np @@ -70,6 +71,15 @@ class FeederConfig(BaseModel): topology_output: str = "topology.json" use_sparse_admittance: bool = False tap_setting: Optional[int] = None + open_lines: Optional[List[str]] = None + + +# Open Lines: + +# "Line.padswitch(r:p9udt496-p9udt527)p9u_166790", +# "Line.padswitch(r:p9udt527-p9udt528)p9u_166794" is better for large SMART-DS. + +# "Line.goab_disswitch(r:p1udt1425-p1udt881)p1u_12301", # SMALL class FeederMapping(BaseModel): @@ -143,6 +153,7 @@ def __init__(self, config: FeederConfig): else: self._feeder_file = config.existing_feeder_file + self.open_lines = config.open_lines self.load_feeder() if self._sensor_location is None: @@ -153,16 +164,16 @@ def __init__(self, config: FeederConfig): def forcast_pv(self, steps: int) -> list: """ - Forecasts day ahead PV generation for the OpenDSS feeder. The OpenDSS file is run and the + Forecasts day ahead PV generation for the OpenDSS feeder. The OpenDSS file is run and the average irradiance is computed over all PV systems for each time step. This average irradiance is used to compute the individual PV system power output """ - cmd = f'Set stepsize={self._simulation_time_step} Number=1' + cmd = f"Set stepsize={self._simulation_time_step} Number=1" dss.Text.Command(cmd) forecast = [] for k in range(steps): dss.Solution.Solve() - + # names of PV systems and forecasted power output pv_names = [] powers = [] @@ -171,7 +182,7 @@ def forcast_pv(self, steps: int) -> list: flag = dss.PVsystems.First() avg_irradiance = dss.PVsystems.IrradianceNow() while flag: - avg_irradiance = (avg_irradiance + dss.PVsystems.IrradianceNow())/2 + avg_irradiance = (avg_irradiance + dss.PVsystems.IrradianceNow()) / 2 flag = dss.PVsystems.Next() # now compute the power output from the evaluated average irradiance @@ -180,7 +191,7 @@ def forcast_pv(self, steps: int) -> list: pv_names.append(f"PVSystem.{dss.PVsystems.Name()}") powers.append(dss.PVsystems.Pmpp() * avg_irradiance) flag = dss.PVsystems.Next() - + forecast.append(xr.DataArray(powers, coords={"ids": pv_names})) return forecast @@ -192,7 +203,7 @@ def snapshot_run(self): assert self._state != OpenDSSState.UNLOADED, f"{self._state}" self.reenable() dss.Text.Command("CalcVoltageBases") - dss.Text.Command("solve mode=snapshot") + dss.Text.Command("solve mode=snapshot number=1") self._state = OpenDSSState.SNAPSHOT_RUN def reenable(self): @@ -302,6 +313,24 @@ def get_node_names(self): """Get node names in order.""" return self._AllNodeNames + def get_bus_coords(self) -> Dict[str, Tuple[float, float]] | None: + """Load bus coordinates from OpenDSS.""" + bus_path = os.path.join(os.path.dirname(self._feeder_file), "Buscoords.dss") + if not os.path.exists(bus_path): + self.bus_coords = None + return self.bus_coords + with open(bus_path, "r") as f: + bus_coord_csv = csv.reader(f, delimiter=" ") + bus_coords = {} + for row in bus_coord_csv: + try: + identifier, x, y = row + bus_coords[identifier] = (float(x), float(y)) + except ValueError as e: + logging.warning(f"Unable to parse row in bus coords: {row}, {e}") + return None + return bus_coords + def load_feeder(self): """Load feeder once downloaded. Relies on legacy mode.""" # Real solution is kvarlimit with kvarmax @@ -337,6 +366,10 @@ def load_feeder(self): if self.tap_setting is not None: # Doesn't work with AutoTrans or 3-winding transformers. dss.Text.Command(f"batchedit transformer..* wdg=2 tap={self.tap_setting}") + + if self.open_lines is not None: + for l in self.open_lines: + self.open_line(l) self._state = OpenDSSState.LOADED def disable_elements(self): @@ -359,7 +392,7 @@ def disabled_run(self): dss.Text.Command("CalcVoltageBases") dss.Text.Command("set maxiterations=20") # solve - dss.Text.Command("solve") + dss.Text.Command("solve number=1") self._state = OpenDSSState.DISABLED_RUN def get_y_matrix(self): @@ -384,7 +417,7 @@ def get_load_y_matrix(self): dss.Text.Command("CalcVoltageBases") dss.Text.Command("set maxiterations=20") # solve - dss.Text.Command("solve") + dss.Text.Command("solve number=1") Ysparse = csc_matrix(dss.YMatrix.getYsparse()) Ymatrix = Ysparse.tocoo() @@ -397,7 +430,7 @@ def get_load_y_matrix(self): dss.Text.Command("CalcVoltageBases") dss.Text.Command("set maxiterations=20") - dss.Text.Command("solve") + dss.Text.Command("solve number=1") self._state = OpenDSSState.SOLVE_AT_TIME return coo_matrix( @@ -432,7 +465,7 @@ def just_solve(self): self._state != OpenDSSState.UNLOADED and self._state != OpenDSSState.DISABLED_RUN ), f"{self._state}" - dss.Text.Command("solve") + dss.Text.Command("solve number=1") def solve(self, hour, second): """Solve at specified time. Must not be unloaded or disabled.""" @@ -445,7 +478,7 @@ def solve(self, hour, second): f"set mode=yearly loadmult=1 number=1 hour={hour} sec={second} " f"stepsize=0" ) - dss.Text.Command("solve") + dss.Text.Command("solve number=1") self._state = OpenDSSState.SOLVE_AT_TIME def _ready_to_load_power(self, static): @@ -568,6 +601,11 @@ def get_PQs_gen(self, static=False): ) return pq_xr.sortby(pq_xr.ids) + def open_line(self, line_name: str): + """Open a line in the circuit.""" + dss.Circuit.SetActiveElement(line_name) + dss.CktElement.Open(2, 0) + def get_PQs_cap(self, static=False): """Get active and reactive power of Capacitors as xarray.""" self._ready_to_load_power(static) diff --git a/LocalFeeder/sender_cosim.py b/LocalFeeder/sender_cosim.py index 4fa5aaf..c2112fc 100644 --- a/LocalFeeder/sender_cosim.py +++ b/LocalFeeder/sender_cosim.py @@ -352,10 +352,13 @@ def go_cosim( h.helicsFederateEnterExecutingMode(vfed) initial_data = get_initial_data(sim, config) + topology_dict = initial_data.topology.dict() + topology_dict["bus_coords"] = sim.get_bus_coords() + topology_json = json.dumps(topology_dict) logger.info("Sending topology and saving to topology.json") with open(config.topology_output, "w") as f: - f.write(initial_data.topology.json()) - pub_topology.publish(initial_data.topology.json()) + f.write(topology_json) + pub_topology.publish(topology_json) # Publish the forecasted PV outputs as a list of MeasurementArray logger.info("Evaluating the forecasted PV") @@ -366,6 +369,9 @@ def go_cosim( granted_time = -1 request_time = 0 + initial_timestamp = datetime.strptime( + config.start_date, "%Y-%m-%d %H:%M:%S" + ) while request_time < int(config.number_of_timesteps): granted_time = h.helicsFederateRequestTime(vfed, request_time) @@ -394,14 +400,15 @@ def go_cosim( for pv_set in pv_sets: sim.set_pv_output(pv_set[0].split(".")[1], pv_set[1], pv_set[2]) + current_hour = 24*(floored_timestamp.date() - initial_timestamp.date()).days + floored_timestamp.hour logger.info( - f"Solve at hour {floored_timestamp.hour} second " + f"Solve at hour {current_hour} second " f"{60*floored_timestamp.minute + floored_timestamp.second}" ) sim.snapshot_run() sim.solve( - floored_timestamp.hour, + current_hour, 60 * floored_timestamp.minute + floored_timestamp.second, ) diff --git a/recorder/record_subscription.py b/recorder/record_subscription.py index 517b484..abd5e38 100644 --- a/recorder/record_subscription.py +++ b/recorder/record_subscription.py @@ -2,6 +2,7 @@ import json import logging from datetime import datetime +import time import helics as h import numpy as np