diff --git a/.github/workflows/docker-test.yml b/.github/workflows/docker-test.yml index 8544016..d834fa5 100644 --- a/.github/workflows/docker-test.yml +++ b/.github/workflows/docker-test.yml @@ -14,8 +14,7 @@ jobs: - name: build docker container shell: bash run: | - echo "$SSH_KEY" > gadal_docker_key - docker build --secret id=gadal_github_key,src=gadal_docker_key --progress=plain -t gadal-example:0.0.0 . + docker build --progress=plain -t oedisi-example:0.0.0 . env: SSH_KEY: ${{secrets.SGIDAL_CLONE_KEY}} DOCKER_BUILDKIT: '1' @@ -23,8 +22,8 @@ jobs: shell: bash run: | mkdir outputs_build - docker volume create --name gadal_output --opt type=none --opt device=$(pwd)/outputs_build --opt o=bind - docker run --rm --mount source=gadal_output,target=/simulation/outputs gadal-example:0.0.0 + docker volume create --name oedisi_output --opt type=none --opt device=$(pwd)/outputs_build --opt o=bind + docker run --rm --mount source=oedisi_output,target=/simulation/outputs oedisi-example:0.0.0 - name: Set up Python ${{ matrix.python-version }} uses: conda-incubator/setup-miniconda@v2 with: @@ -34,9 +33,7 @@ jobs: shell: bash -l {0} run: | pip install matplotlib pyarrow numpy matplotlib pandas - eval `ssh-agent -s` - ssh-add - <<< '${{ secrets.SGIDAL_CLONE_KEY }}' - pip install git+ssh://git@github.com/openEDI/GADAL@v0.2.2 + pip install oedisi==1.0.0 python post_analysis.py outputs_build - name: Archive logs diff --git a/.github/workflows/test-api.yml b/.github/workflows/test-api.yml index afefbad..c3641a8 100644 --- a/.github/workflows/test-api.yml +++ b/.github/workflows/test-api.yml @@ -24,14 +24,11 @@ jobs: shell: bash -l {0} run: | pip install -r requirements.txt - eval `ssh-agent -s` - ssh-add - <<< '${{ secrets.SGIDAL_CLONE_KEY }}' - pip install git+ssh://git@github.com/openEDI/GADAL@v0.2.2 - name: Run example shell: bash -l {0} run: | - python test_full_systems.py - helics run --path=build/test_system_runner.json + oedisi build --system scenarios/docker_system.json + oedisi run python post_analysis.py - name: Archive logs uses: actions/upload-artifact@v2 diff --git a/.github/workflows/test-local-feeder.yml b/.github/workflows/test-local-feeder.yml new file mode 100644 index 0000000..a6366b9 --- /dev/null +++ b/.github/workflows/test-local-feeder.yml @@ -0,0 +1,35 @@ +name: TestLocalFeeder + +on: [push] +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest] + python-version: ['3.10'] + #include: + #- os: ubuntu-latest + #python-version: 3.10 + + steps: + - uses: actions/checkout@v2 + - name: Set up Python ${{ matrix.python-version }} + uses: conda-incubator/setup-miniconda@v2 + with: + auto-update-conda: true + python-version: ${{ matrix.python-version }} + - name: Install python dependencies + shell: bash -l {0} + run: | + pip install -r requirements.txt + pip install plotille pytest + - name: Run pytest on LocalFeeder + shell: bash -l {0} + run: | + cd LocalFeeder + pytest -s . + env: + PYTEST_ADDOPTS: "--color=yes" + PYTHONIOENCODING: UTF-8 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..d06cec1 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,31 @@ +# See https://pre-commit.com for more information +# See https://pre-commit.com/hooks.html for more hooks +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v3.2.0 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-yaml + - id: check-added-large-files +- repo: https://github.com/pycqa/flake8 + rev: 6.0.0 + hooks: + - id: flake8 +- repo: https://github.com/psf/black + rev: 19.10b0 + hooks: + - id: black +- repo: https://github.com/PyCQA/isort + rev: 5.9.3 + hooks: + - id: isort +- repo: https://github.com/pre-commit/mirrors-mypy + rev: v0.981 + hooks: + - id: mypy + additional_dependencies: [pydantic] # add if use pydantic +- repo: https://github.com/pycqa/pydocstyle + rev: 6.3.0 + hooks: + - id: pydocstyle diff --git a/Dockerfile b/Dockerfile index a8ecf1b..fbfc6a2 100644 --- a/Dockerfile +++ b/Dockerfile @@ -3,14 +3,11 @@ FROM python:3.10.6-slim-bullseye RUN apt-get update && apt-get install -y git ssh RUN mkdir -p /root/.ssh -ENV GIT_SSH_COMMAND="ssh -i /run/secrets/gadal_github_key" -RUN ssh-keyscan -t rsa github.com >> ~/.ssh/known_hosts -RUN --mount=type=secret,id=gadal_github_key pip install git+ssh://git@github.com/openEDI/GADAL@v0.2.2 WORKDIR /simulation -COPY test_full_systems.py . COPY scenarios/docker_system.json docker_system.json +COPY components.json . COPY LocalFeeder LocalFeeder COPY README.md . COPY measuring_federate measuring_federate @@ -18,8 +15,9 @@ COPY wls_federate wls_federate COPY recorder recorder RUN mkdir -p outputs build -RUN python test_full_systems.py --system docker_system.json COPY requirements.txt . RUN pip install -r requirements.txt -ENTRYPOINT ["helics", "run", "--path=build/test_system_runner.json"] + +RUN oedisi build --system docker_system.json +ENTRYPOINT ["oedisi", "run"] diff --git a/LocalFeeder/.gitignore b/LocalFeeder/.gitignore index 595781d..f594f53 100644 --- a/LocalFeeder/.gitignore +++ b/LocalFeeder/.gitignore @@ -1,4 +1,7 @@ base_*.csv +base_*.npz +opendss +profiles # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/LocalFeeder/FeederSimulator.py b/LocalFeeder/FeederSimulator.py index 9824d00..fac257e 100644 --- a/LocalFeeder/FeederSimulator.py +++ b/LocalFeeder/FeederSimulator.py @@ -1,37 +1,39 @@ - -# -*- coding: utf-8 -*- -# import helics as h -import opendssdirect as dss -import numpy as np -import csv -import time -from time import strptime -from scipy.sparse import coo_matrix +"""Core class to abstract OpenDSS into Feeder class.""" +import json +import logging +import math import os import random -import math -import logging -import json -import sys +import time +from enum import Enum +from time import strptime +from typing import Dict, List, Set, Optional + import boto3 +import numpy as np +import opendssdirect as dss +import xarray as xr from botocore import UNSIGNED from botocore.config import Config - +from oedisi.types.data_types import Command, InverterControl, InverterControlMode from pydantic import BaseModel +from scipy.sparse import coo_matrix, csc_matrix from dss_functions import ( - snapshot_run, parse_Ymatrix, get_loads, get_pvSystems, get_Generator, - get_capacitors, get_voltages, get_y_matrix_file, get_vnom, get_vnom2 + get_capacitors, + get_generators, + get_loads, + get_pvsystems, + get_voltages, ) -import dss_functions logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) logger.setLevel(logging.INFO) + def permutation(from_list, to_list): - """ - Create permutation representing change in from_list to to_list + """Create permutation representing change in from_list to to_list. Specifically, if `permute = permutation(from_list, to_list)`, then `permute[i] = j` means that `from_list[i] = to_list[j]`. @@ -41,187 +43,316 @@ def permutation(from_list, to_list): You may view the permutation as a function from `from_list` to `to_list`. """ - #return [to_list.find(v) for v in enumerate(from_list)] index_map = {v: i for i, v in enumerate(to_list)} return [index_map[v] for v in from_list] - -def check_node_order(l1, l2): - logger.debug('check order ' + str(l1 == l2)) - - class FeederConfig(BaseModel): + """JSON configuration. Special cases S3 sources right now.""" + name: str use_smartds: bool = False - profile_location: str - opendss_location: str - sensor_location: str = '' + profile_location: str + opendss_location: str + existing_feeder_file: Optional[str] = None + sensor_location: Optional[str] = None start_date: str number_of_timesteps: float - run_freq_sec: float = 15*60 + run_freq_sec: float = 15 * 60 start_time_index: int = 0 topology_output: str = "topology.json" - use_sparse_admittance = False + use_sparse_admittance: bool = False -class FeederSimulator(object): - """ A simple class that handles publishing the solar forecast - """ +class OpenDSSState(Enum): + """Enum of all OpenDSSStates traversed in a simulation.""" - def __init__(self, config): - """ Create a ``FeederSimulator`` object + UNLOADED = 1 + LOADED = 2 + SNAPSHOT_RUN = 3 + SOLVE_AT_TIME = 4 + DISABLED_RUN = 5 + DISABLED_SOLVE = 6 + DISABLED = 7 - """ - self._feeder_file = None - self._simulation_time_step = None + +class FeederSimulator(object): + """A simple class that handles publishing the solar forecast.""" + + # Private variables initialized later + _feeder_file: str + _AllNodeNames: List[str] + _source_indexes: List[int] + _nodes_index: List[int] + _name_index_dict: Dict[str, int] + _inverter_to_pvsystems: Dict[str, Set[str]] + _pvsystem_to_inverter: Dict[str, str] + _pvsystems: Set[str] + _inverters: Set[str] + _inverter_counter: int + _xycurve_counter: int + + def __init__(self, config: FeederConfig): + """Create a ``FeederSimulator`` object.""" + self._state = OpenDSSState.UNLOADED self._opendss_location = config.opendss_location self._profile_location = config.profile_location self._sensor_location = config.sensor_location self._use_smartds = config.use_smartds - self._circuit=None - self._AllNodeNames=None - self._source_indexes=None - # self._capacitors=[] - self._capNames = [] - self._regNames = [] + self._inverter_to_pvsystems = {} + self._pvsystem_to_inverter = {} + self._inverters = set() + self._inverter_counter = 0 + self._xycurve_counter = 0 - # timegm(strptime('2019-07-23 14:50:00 GMT', '%Y-%m-%d %H:%M:%S %Z')) - self._start_time = int(time.mktime(strptime(config.start_date, '%Y-%m-%d %H:%M:%S'))) + self._start_time = int( + time.mktime(strptime(config.start_date, "%Y-%m-%d %H:%M:%S")) + ) self._run_freq_sec = config.run_freq_sec self._simulation_step = config.start_time_index self._number_of_timesteps = config.number_of_timesteps self._vmult = 0.001 - self._nodes_index = [] - self._name_index_dict = {} + self._simulation_time_step = "15m" + if config.existing_feeder_file is None: + if self._use_smartds: + self._feeder_file = os.path.join("opendss", "Master.dss") + self.download_data("oedi-data-lake", update_loadshape_location=True) + else: + self._feeder_file = os.path.join("opendss", "master.dss") + self.download_data("gadal") + else: + self._feeder_file = config.existing_feeder_file + + self.load_feeder() - if self._use_smartds: - self.download_data('oedi-data-lake','Master.dss',True) - self.load_feeder() + if self._sensor_location is None: self.create_measurement_lists() - else: - self.download_data('gadal','master.dss') - self.load_feeder() - def download_data(self, bucket_name, master_name, update_loadshape_location=False): + self.snapshot_run() + assert self._state == OpenDSSState.SNAPSHOT_RUN, f"{self._state}" + + def snapshot_run(self): + """Run snapshot of simuation without specifying a time. - #Equivalent to --no-sign-request - s3_resource = boto3.resource('s3',config=Config(signature_version=UNSIGNED)) + Used for initialization. + """ + assert self._state != OpenDSSState.UNLOADED, f"{self._state}" + self.reenable() + dss.Text.Command("CalcVoltageBases") + dss.Text.Command("solve mode=snapshot") + self._state = OpenDSSState.SNAPSHOT_RUN + + def reenable(self): + dss.Text.Command("Batchedit Load..* enabled=yes") + dss.Text.Command("Batchedit Vsource..* enabled=yes") + dss.Text.Command("Batchedit Isource..* enabled=yes") + dss.Text.Command("Batchedit Generator..* enabled=yes") + dss.Text.Command("Batchedit PVsystem..* enabled=yes") + dss.Text.Command("Batchedit Capacitor..* enabled=yes") + dss.Text.Command("Batchedit Storage..* enabled=no") + + def download_data(self, bucket_name, update_loadshape_location=False): + """Download data from bucket path.""" + logging.info(f"Downloading from bucket {bucket_name}") + # Equivalent to --no-sign-request + s3_resource = boto3.resource("s3", config=Config(signature_version=UNSIGNED)) bucket = s3_resource.Bucket(bucket_name) opendss_location = self._opendss_location profile_location = self._profile_location sensor_location = self._sensor_location - self._feeder_file = os.path.join('opendss',master_name) - self._simulation_time_step = '15m' for obj in bucket.objects.filter(Prefix=opendss_location): - output_location = os.path.join('opendss',obj.key.replace(opendss_location,'').strip('/')) + output_location = os.path.join( + "opendss", obj.key.replace(opendss_location, "").strip("/") + ) os.makedirs(os.path.dirname(output_location), exist_ok=True) - bucket.download_file(obj.key,output_location) + bucket.download_file(obj.key, output_location) - modified_loadshapes = '' - os.makedirs(os.path.join('profiles'), exist_ok=True) + modified_loadshapes = "" + os.makedirs(os.path.join("profiles"), exist_ok=True) if update_loadshape_location: all_profiles = set() - with open(os.path.join('opendss','LoadShapes.dss'),'r') as fp_loadshapes: + with open(os.path.join("opendss", "LoadShapes.dss"), "r") as fp_loadshapes: for row in fp_loadshapes.readlines(): - new_row = row.replace('../','') - new_row = new_row.replace('file=','file=../') - for token in new_row.split(' '): - if token.startswith('(file='): - location = token.split('=../profiles/')[1].strip().strip(')') + new_row = row.replace("../", "") + new_row = new_row.replace("file=", "file=../") + for token in new_row.split(" "): + if token.startswith("(file="): + location = ( + token.split("=../profiles/")[1].strip().strip(")") + ) all_profiles.add(location) - modified_loadshapes=modified_loadshapes+new_row - with open(os.path.join('opendss','LoadShapes.dss'),'w') as fp_loadshapes: + modified_loadshapes = modified_loadshapes + new_row + with open(os.path.join("opendss", "LoadShapes.dss"), "w") as fp_loadshapes: fp_loadshapes.write(modified_loadshapes) for profile in all_profiles: - s3_location = f'{profile_location}/{profile}' - bucket.download_file(s3_location,os.path.join('profiles',profile)) + s3_location = f"{profile_location}/{profile}" + bucket.download_file(s3_location, os.path.join("profiles", profile)) else: for obj in bucket.objects.filter(Prefix=profile_location): - output_location = os.path.join('profiles',obj.key.replace(profile_location,'').strip('/')) + output_location = os.path.join( + "profiles", obj.key.replace(profile_location, "").strip("/") + ) os.makedirs(os.path.dirname(output_location), exist_ok=True) - bucket.download_file(obj.key,output_location) + bucket.download_file(obj.key, output_location) - if sensor_location != '': - output_location = os.path.join('sensors',os.path.basename(sensor_location)) + if sensor_location is not None: + output_location = os.path.join("sensors", os.path.basename(sensor_location)) if not os.path.exists(os.path.dirname(output_location)): os.makedirs(os.path.dirname(output_location)) - bucket.download_file(sensor_location,output_location) - - def create_measurement_lists(self, - percent_voltage=75, - percent_real=75, - percent_reactive=75, - voltage_seed=1, - real_seed=2, - reactive_seed=3 - ): - + bucket.download_file(sensor_location, output_location) + + def create_measurement_lists( + self, + percent_voltage=75, + percent_real=75, + voltage_seed=1, + real_seed=2, + reactive_seed=3, + ): + """Initialize list of sensor locations for the measurement federate.""" random.seed(voltage_seed) - os.makedirs('sensors', exist_ok=True) - voltage_subset = random.sample(self._AllNodeNames,math.floor(len(self._AllNodeNames)*float(percent_voltage)/100)) - with open(os.path.join('sensors','voltage_ids.json'),'w') as fp: - json.dump(voltage_subset,fp,indent=4) + os.makedirs("sensors", exist_ok=True) + voltage_subset = random.sample( + self._AllNodeNames, + math.floor(len(self._AllNodeNames) * float(percent_voltage) / 100), + ) + with open(os.path.join("sensors", "voltage_ids.json"), "w") as fp: + json.dump(voltage_subset, fp, indent=4) random.seed(real_seed) - real_subset = random.sample(self._AllNodeNames,math.floor(len(self._AllNodeNames)*float(percent_real)/100)) - with open(os.path.join('sensors','real_ids.json'),'w') as fp: - json.dump(real_subset,fp,indent=4) + real_subset = random.sample( + self._AllNodeNames, + math.floor(len(self._AllNodeNames) * float(percent_real) / 100), + ) + with open(os.path.join("sensors", "real_ids.json"), "w") as fp: + json.dump(real_subset, fp, indent=4) random.seed(reactive_seed) - reactive_subset = random.sample(self._AllNodeNames,math.floor(len(self._AllNodeNames)*float(percent_voltage)/100)) - with open(os.path.join('sensors','reactive_ids.json'),'w') as fp: - json.dump(reactive_subset,fp,indent=4) - - def snapshot_run(self): - snapshot_run(dss) + reactive_subset = random.sample( + self._AllNodeNames, + math.floor(len(self._AllNodeNames) * float(percent_voltage) / 100), + ) + with open(os.path.join("sensors", "reactive_ids.json"), "w") as fp: + json.dump(reactive_subset, fp, indent=4) def get_circuit_name(self): + """Get name of current opendss circuit.""" return self._circuit.Name() def get_source_indices(self): + """Get indcies of slack buses.""" return self._source_indexes def get_node_names(self): + """Get node names in order.""" return self._AllNodeNames - def load_feeder(self): - dss.Basic.LegacyModels(True) - dss.run_command("clear") - result = dss.run_command("redirect " + self._feeder_file) - if not result == '': - raise ValueError("Feeder not loaded: "+result) + """Load feeder once downloaded. Relies on legacy mode.""" + # Real solution is kvarlimit with kvarmax + dss.Basic.LegacyModels(True) + dss.Text.Command("clear") + dss.Text.Command("redirect " + self._feeder_file) + result = dss.Text.Result() + if not result == "": + raise ValueError("Feeder not loaded: " + result) self._circuit = dss.Circuit self._AllNodeNames = self._circuit.YNodeOrder() self._node_number = len(self._AllNodeNames) self._nodes_index = [self._AllNodeNames.index(ii) for ii in self._AllNodeNames] - self._name_index_dict = {ii: self._AllNodeNames.index(ii) for ii in self._AllNodeNames} + self._name_index_dict = { + ii: self._AllNodeNames.index(ii) for ii in self._AllNodeNames + } self._source_indexes = [] for Source in dss.Vsources.AllNames(): - self._circuit.SetActiveElement('Vsource.' + Source) + self._circuit.SetActiveElement("Vsource." + Source) Bus = dss.CktElement.BusNames()[0].upper() for phase in range(1, dss.CktElement.NumPhases() + 1): - self._source_indexes.append(self._AllNodeNames.index(Bus.upper() + '.' + str(phase))) - + self._source_indexes.append( + self._AllNodeNames.index(Bus.upper() + "." + str(phase)) + ) self.setup_vbase() + self._pvsystems = set() + for PV in get_pvsystems(dss): + self._pvsystems.add("PVSystem." + PV["name"]) + self._state = OpenDSSState.LOADED + + def disable_elements(self): + """Disable most elements. Used in disabled_run.""" + assert self._state != OpenDSSState.UNLOADED, f"{self._state}" + # dss.Text.Command("batchedit transformer..* wdg=2 tap=1") + dss.Text.Command("batchedit regcontrol..* enabled=false") + dss.Text.Command("batchedit vsource..* enabled=false") + dss.Text.Command("batchedit isource..* enabled=false") + dss.Text.Command("batchedit load..* enabled=false") + dss.Text.Command("batchedit generator..* enabled=false") + dss.Text.Command("batchedit pvsystem..* enabled=false") + dss.Text.Command("Batchedit Capacitor..* enabled=false") + dss.Text.Command("batchedit storage..* enabled=false") + self._state = OpenDSSState.DISABLED + + def disabled_run(self): + """Disable most elements and solve. Used for most Y-matrix needs.""" + self.disable_elements() + assert self._state == OpenDSSState.DISABLED, f"{self._state}" + dss.Text.Command("CalcVoltageBases") + dss.Text.Command("set maxiterations=20") + # solve + dss.Text.Command("solve") + self._state = OpenDSSState.DISABLED_RUN + def get_y_matrix(self): - get_y_matrix_file(dss) - Ymatrix = parse_Ymatrix('base_ysparse.csv', self._node_number) + """Calculate Y-matrix as a coo-matrix. Disables some elements.""" + self.disabled_run() + self._state = OpenDSSState.DISABLED_RUN + + Ysparse = csc_matrix(dss.YMatrix.getYsparse()) + Ymatrix = Ysparse.tocoo() new_order = self._circuit.YNodeOrder() permute = np.array(permutation(new_order, self._AllNodeNames)) - #inv_permute = np.array(permutation(self._AllNodeNames, new_order)) - return coo_matrix((Ymatrix.data, (permute[Ymatrix.row], permute[Ymatrix.col])), shape=Ymatrix.shape) + return coo_matrix( + (Ymatrix.data, (permute[Ymatrix.row], permute[Ymatrix.col])), + shape=Ymatrix.shape, + ) + + def get_load_y_matrix(self): + """Calculate Y-matrix as a coo-matrix. Disables most except load.""" + assert self._state == OpenDSSState.SOLVE_AT_TIME, f"{self._state}" + self.disable_elements() + dss.Text.Command("batchedit Load..* enabled=true") + dss.Text.Command("CalcVoltageBases") + dss.Text.Command("set maxiterations=20") + # solve + dss.Text.Command("solve") + + Ysparse = csc_matrix(dss.YMatrix.getYsparse()) + Ymatrix = Ysparse.tocoo() + new_order = self._circuit.YNodeOrder() + permute = np.array(permutation(new_order, self._AllNodeNames)) + + dss.Text.Command("batchedit Load..* enabled=false") + self._state = OpenDSSState.DISABLED_RUN + self.reenable() + + dss.Text.Command("CalcVoltageBases") + dss.Text.Command("set maxiterations=20") + dss.Text.Command("solve") + self._state = OpenDSSState.SOLVE_AT_TIME + return coo_matrix( + (Ymatrix.data, (permute[Ymatrix.row], permute[Ymatrix.col])), + shape=Ymatrix.shape, + ) def setup_vbase(self): + """Load base voltages into feeder.""" self._Vbase_allnode = np.zeros((self._node_number), dtype=np.complex_) self._Vbase_allnode_dict = {} for ii, node in enumerate(self._AllNodeNames): @@ -229,199 +360,395 @@ def setup_vbase(self): self._Vbase_allnode[ii] = dss.Bus.kVBase() * 1000 self._Vbase_allnode_dict[node] = self._Vbase_allnode[ii] - def get_G_H(self, Y11_inv): - Vnom = self.get_vnom() - # ys=Y11 - # R = np.linalg.inv(ys).real - # X = np.linalg.inv(ys).imag - R=Y11_inv.real - X=Y11_inv.imag - G = (R * np.diag(np.cos(np.angle(Vnom)) / abs(Vnom)) - X * np.diag(np.sin(np.angle(Vnom)) / Vnom)).real - H = (X * np.diag(np.cos(np.angle(Vnom)) / abs(Vnom)) - R * np.diag(np.sin(np.angle(Vnom)) / Vnom)).real - return Vnom, G, H - - def get_vnom2(self): - _Vnom, Vnom_dict = get_vnom2(dss) - Vnom = np.zeros((len(self._AllNodeNames)), dtype=np.complex_) - for voltage_name in Vnom_dict.keys(): - Vnom[self._name_index_dict[voltage_name]] = Vnom_dict[voltage_name] - # Vnom(1: 3) = []; - Vnom = np.concatenate((Vnom[:self._source_indexes[0]], Vnom[self._source_indexes[-1] + 1:])) - return Vnom - - def get_vnom(self): - _Vnom, Vnom_dict = get_vnom(dss) - Vnom = np.zeros((len(self._AllNodeNames)), dtype=np.complex_) - # print([name_voltage_dict.keys()][:5]) - for voltage_name in Vnom_dict.keys(): - Vnom[self._name_index_dict[voltage_name]] = Vnom_dict[voltage_name] - # Vnom(1: 3) = []; - logger.debug(Vnom[self._source_indexes[0]:self._source_indexes[-1]]) - Vnom = np.concatenate((Vnom[:self._source_indexes[0]], Vnom[self._source_indexes[-1] + 1:])) - Vnom = np.abs(Vnom) - return Vnom - - def get_PQs_load(self,static=False): - num_nodes = len(self._name_index_dict.keys()) - - PQ_names= ['' for i in range(num_nodes)] - PQ_types= ['' for i in range(num_nodes)] - PQ_load = np.zeros((num_nodes), dtype=np.complex_) - for ld in get_loads(dss,self._circuit): - self._circuit.SetActiveElement('Load.' + ld["name"]) - for ii in range(len(ld['phases'])): - name = ld['bus1'].upper() + '.' + ld['phases'][ii] - index = self._name_index_dict[name] + def initial_disabled_solve(self): + """If run is disabled, then we can still solve at 0.0.""" + assert self._state == OpenDSSState.DISABLED_RUN, f"{self._state}" + hour = 0 + second = 0 + dss.Text.Command( + f"set mode=yearly loadmult=1 number=1 hour={hour} sec={second} " + f"stepsize=0" + ) + dss.Text.Command("solve") + self._state = OpenDSSState.DISABLED_SOLVE + + def just_solve(self): + """Solve without setting time or anything. Useful for commands.""" + assert ( + self._state != OpenDSSState.UNLOADED + and self._state != OpenDSSState.DISABLED_RUN + ), f"{self._state}" + dss.Text.Command("solve") + + def solve(self, hour, second): + """Solve at specified time. Must not be unloaded or disabled.""" + assert ( + self._state != OpenDSSState.UNLOADED + and self._state != OpenDSSState.DISABLED_RUN + ), f"{self._state}" + + dss.Text.Command( + f"set mode=yearly loadmult=1 number=1 hour={hour} sec={second} " + f"stepsize=0" + ) + dss.Text.Command("solve") + self._state = OpenDSSState.SOLVE_AT_TIME + + def _ready_to_load_power(self, static): + """Check if opendss state can actually calculate power.""" + if static: + assert self._state != OpenDSSState.UNLOADED, f"{self._state}" + else: + assert self._state == OpenDSSState.SOLVE_AT_TIME, f"{self._state}" + + def get_PQs_load(self, static=False): + """Get active and reactive power of loads as xarray.""" + self._ready_to_load_power(static) + + all_node_names = set(self._AllNodeNames) + PQs: List[complex] = [] + node_names: List[str] = [] + pq_names: List[str] = [] + # PQ_load = np.zeros((num_nodes), dtype=np.complex_) + for ld in get_loads(dss, self._circuit): + self._circuit.SetActiveElement("Load." + ld["name"]) + current_pq_name = dss.CktElement.Name() + for ii in range(len(ld["phases"])): + node_name = ld["bus1"].upper() + "." + ld["phases"][ii] + assert ( + node_name in all_node_names + ), f"{node_name} for {current_pq_name} not found" if static: - power = complex(ld['kW'],ld['kVar']) - PQ_load[index] += power/len(ld['phases']) + power = complex(ld["kW"], ld["kVar"]) + PQs.append(power / len(ld["phases"])) else: power = dss.CktElement.Powers() - PQ_load[index] += complex(power[2 * ii], power[2 * ii + 1]) - PQ_names[index] = name - PQ_types[index] = 'Load' - - return PQ_load,PQ_names,PQ_types - - def get_PQs_pv(self,static=False): - num_nodes = len(self._name_index_dict.keys()) - - PQ_names= ['' for i in range(num_nodes)] - PQ_types= ['' for i in range(num_nodes)] - PQ_PV = np.zeros((num_nodes), dtype=np.complex_) - for PV in get_pvSystems(dss): - bus = PV["bus"].split('.') + PQs.append(complex(power[2 * ii], power[2 * ii + 1])) + pq_names.append(current_pq_name) + node_names.append(node_name) + pq_xr = xr.DataArray( + PQs, + dims=("eqnode",), + coords={ + "equipment_ids": ("eqnode", pq_names), + "ids": ("eqnode", node_names), + }, + ) + return pq_xr.sortby(pq_xr.ids) + + def get_PQs_pv(self, static=False): + """Get active and reactive power of PVSystems as xarray.""" + self._ready_to_load_power(static) + + all_node_names = set(self._AllNodeNames) + PQs: List[complex] = [] + node_names: List[str] = [] + pq_names: List[str] = [] + for PV in get_pvsystems(dss): + bus = PV["bus"].split(".") if len(bus) == 1: - bus = bus + ['1', '2', '3'] - self._circuit.SetActiveElement('PVSystem.'+PV["name"]) + bus = bus + ["1", "2", "3"] + self._circuit.SetActiveElement("PVSystem." + PV["name"]) + current_pq_name = dss.CktElement.Name() for ii in range(len(bus) - 1): - name = bus[0].upper() + '.' + bus[ii + 1] - index = self._name_index_dict[name] + node_name = bus[0].upper() + "." + bus[ii + 1] + assert ( + node_name in all_node_names + ), f"{node_name} for {current_pq_name} not found" if static: - power = complex(-1*PV['kW'],-1*PV['kVar']) #-1 because injecting - PQ_PV[index] += power/(len(bus)-1) + power = complex( + -1 * PV["kW"], -1 * PV["kVar"] + ) # -1 because injecting + PQs.append(power / (len(bus) - 1)) else: power = dss.CktElement.Powers() - PQ_PV[index] += complex(power[2 * ii], power[2 * ii + 1]) - PQ_names[index] = name - PQ_types[index] = 'PVSystem' - return PQ_PV,PQ_names,PQ_types - - def get_PQs_gen(self,static=False): - num_nodes = len(self._name_index_dict.keys()) - - PQ_names= ['' for i in range(num_nodes)] - PQ_types= ['' for i in range(num_nodes)] - PQ_gen = np.zeros((num_nodes), dtype=np.complex_) - for PV in get_Generator(dss): - bus = PV["bus"] - self._circuit.SetActiveElement('Generator.'+PV["name"]) + PQs.append(complex(power[2 * ii], power[2 * ii + 1])) + pq_names.append(current_pq_name) + node_names.append(node_name) + pq_xr = xr.DataArray( + PQs, + dims=("eqnode",), + coords={ + "equipment_ids": ("eqnode", pq_names), + "ids": ("eqnode", node_names), + }, + ) + return pq_xr.sortby(pq_xr.ids) + + def get_PQs_gen(self, static=False): + """Get active and reactive power of Generators as xarray.""" + self._ready_to_load_power(static) + + all_node_names = set(self._AllNodeNames) + PQs: List[complex] = [] + node_names: List[str] = [] + pq_names: List[str] = [] + for gen in get_generators(dss): + bus = gen["bus"].split(".") + if len(bus) == 1: + bus = bus + ["1", "2", "3"] + self._circuit.SetActiveElement("Generator." + gen["name"]) + current_pq_name = dss.CktElement.Name() for ii in range(len(bus) - 1): - name = bus[0].upper() + '.' + bus[ii + 1] - index = self._name_index_dict[name] + node_name = bus[0].upper() + "." + bus[ii + 1] + assert ( + node_name in all_node_names + ), f"{node_name} for {current_pq_name} not found" if static: - power = complex(-1*PV['kW'],-1*PV['kVar']) #-1 because injecting - PQ_gen[index] += power/(len(bus)-1) + power = complex( + -1 * gen["kW"], -1 * gen["kVar"] + ) # -1 because injecting + PQs.append(power / (len(bus) - 1)) else: power = dss.CktElement.Powers() - PQ_gen[index] += complex(power[2 * ii], power[2 * ii + 1]) - PQ_names[index] = name - PQ_types[index] = 'Generator' - return PQ_gen,PQ_names,PQ_types - - def get_PQs_cap(self,static=False): - num_nodes = len(self._name_index_dict.keys()) - - PQ_names= ['' for i in range(num_nodes)] - PQ_types= ['' for i in range(num_nodes)] - PQ_cap = np.zeros((num_nodes), dtype=np.complex_) + PQs.append(complex(power[2 * ii], power[2 * ii + 1])) + pq_names.append(current_pq_name) + node_names.append(node_name) + pq_xr = xr.DataArray( + PQs, + dims=("eqnode",), + coords={ + "equipment_ids": ("eqnode", pq_names), + "ids": ("eqnode", node_names), + }, + ) + return pq_xr.sortby(pq_xr.ids) + + def get_PQs_cap(self, static=False): + """Get active and reactive power of Capacitors as xarray.""" + self._ready_to_load_power(static) + + all_node_names = set(self._AllNodeNames) + PQs: List[complex] = [] + node_names: List[str] = [] + pq_names: List[str] = [] for cap in get_capacitors(dss): + current_pq_name = cap["name"] for ii in range(cap["numPhases"]): - name = cap["busname"].upper() + '.' + cap["busphase"][ii] - index = self._name_index_dict[name] + node_name = cap["busname"].upper() + "." + cap["busphase"][ii] + assert ( + node_name in all_node_names + ), f"{node_name} for {current_pq_name} not found" if static: - power = complex(0,-1*cap['kVar']) #-1 because it's injected into the grid - PQ_cap[index] += power/cap["numPhases"] + power = complex( + 0, -1 * cap["kVar"] + ) # -1 because it's injected into the grid + PQs.append(power / cap["numPhases"]) else: - PQ_cap[index] = complex(0,cap["power"][2 * ii + 1]) - PQ_names[index] = name - PQ_types[index] = 'Capacitor' - - return PQ_cap,PQ_names,PQ_types - - - - - def get_loads(self): - loads = get_loads(dss, self._circuit) - self._load_power = np.zeros((len(self._AllNodeNames)), dtype=np.complex_) - load_names = [] - load_powers = [] - load = loads[0] - for load in loads: - for phase in load['phases']: - self._load_power[ - self._name_index_dict[load['bus1'].upper() + '.' + phase] - ] = complex(load['power'][0], load['power'][1]) - load_names.append(load['bus1'].upper() + '.' + phase) - load_powers.append(complex(load['power'][0], load['power'][1])) - return self._load_power, load_names, load_powers - - def get_load_sizes(self): - return dss_functions.get_load_sizes(dss,self._loads) + PQs.append(complex(0, cap["power"][2 * ii + 1])) + pq_names.append(current_pq_name) + node_names.append(node_name) + pq_xr = xr.DataArray( + PQs, + dims=("eqnode",), + coords={ + "equipment_ids": ("eqnode", pq_names), + "ids": ("eqnode", node_names), + }, + ) + return pq_xr.sortby(pq_xr.ids) + + def get_base_voltages(self): + """Get base voltages xarray. Can be uesd anytime.""" + return xr.DataArray(self._Vbase_allnode, {"ids": self._AllNodeNames}) + + def get_disabled_solve_voltages(self): + """Get voltage xarray when elements are disabled.""" + assert self._state == OpenDSSState.DISABLED_SOLVE, f"{self._state}" + return self._get_voltages() + + def get_voltages_snapshot(self): + """Get voltage xarray in snapshot run.""" + assert self._state == OpenDSSState.SNAPSHOT_RUN, f"{self._state}" + return self._get_voltages() def get_voltages_actual(self): - ''' - - :return voltages in actual values: - ''' - _, name_voltage_dict = get_voltages(self._circuit) + """Get voltages xarray at current time.""" + assert self._state == OpenDSSState.SOLVE_AT_TIME, f"{self._state}" + return self._get_voltages() + + def _get_voltages(self): + """Get powers and put in order of self._AllNodeNames.""" + assert ( + self._state != OpenDSSState.DISABLED_RUN + and self._state != OpenDSSState.UNLOADED + ), f"{self._state}" + name_voltage_dict = get_voltages(self._circuit) res_feeder_voltages = np.zeros((len(self._AllNodeNames)), dtype=np.complex_) for voltage_name in name_voltage_dict.keys(): - res_feeder_voltages[self._name_index_dict[voltage_name]] = name_voltage_dict[voltage_name] + res_feeder_voltages[ + self._name_index_dict[voltage_name] + ] = name_voltage_dict[voltage_name] - return res_feeder_voltages + return xr.DataArray( + res_feeder_voltages, {"ids": list(name_voltage_dict.keys())} + ) - def get_voltages(self): - ''' + def change_obj(self, change_commands: List[Command]): + """set/get an object property. + + Parameters + ---------- + change_commands: List[Command] + + + Examples + -------- + ``change_obj([Command('PVsystem.pv1','kVAr',25)])`` + """ + assert self._state != OpenDSSState.UNLOADED, f"{self._state}" + for entry in change_commands: + dss.Circuit.SetActiveElement( + entry.obj_name + ) # make the required element as active element + # dss.CktElement.Properties(entry.obj_property).Val = entry.val + # dss.Properties.Value(entry.obj_property, str(entry.val)) + properties = dss.CktElement.AllPropertyNames() + element_name = dss.CktElement.Name() + assert entry.obj_property.lower() in map( + lambda x: x.lower(), properties + ), f"{entry.obj_property} not in {properties} for {element_name}" + dss.Text.Command(f"{entry.obj_name}.{entry.obj_property}={entry.val}") + + def create_inverter(self, pvsystem_set: Set[str]): + """Create new inverter from set of pvsystem. + + Parameters + ---------- + pvsystem_set: Set[str] + Set of pvsystems to assign to inverter. Currently only 1 or all + can be added. This cannot be changed or deleted at present. + + Returns + ------- + name of inverter with ``InvControl.`` prefix: str + + Warnings and Caveats + -------------------- + Using multiple pvsystems will issue a warning, since + more OpenDSS nodes can be created, which will break + more getter functions. + + New inverters cannot be assigned new lists. + + This only works with the legacy InvControl settings + """ + assert all( + pvsystem not in self._pvsystem_to_inverter and pvsystem in self._pvsystems + for pvsystem in pvsystem_set + ), f"PVsystem(s) {pvsystem_set} is already assigned inverter or may not exist" + name = f"InvControl.invgenerated{self._inverter_counter}" + # May need PVsystemlist + assert all( + pv_name.split(".")[0].lower() == "pvsystem" for pv_name in pvsystem_set + ) + if pvsystem_set == self._pvsystems: + # This provides more stable behavior in the "default" case. + pvlist = "" + else: + if len(pvsystem_set) != 1: + logging.error( + """Controlling mulitple pvsystems manually results in unstable + behavior when the number of phases differ""" + ) + pvlist = ", ".join(pv_name.split(".")[1] for pv_name in pvsystem_set) + dss.Text.Command(f"New {name} PVsystemList=[{pvlist}]") + self._inverter_counter += 1 + self._inverters.add(name) + for pvsystem in pvsystem_set: + self._pvsystem_to_inverter[pvsystem] = name + self._inverter_to_pvsystems[name] = pvsystem_set + return name + + def create_xy_curve(self, x, y): + """Create xy curve given two lists or arrays of floats. + + Returns + ------- + name : str + """ + name = f"XYcurve.xygenerated{self._xycurve_counter}" + npts = len(x) + assert len(x) == len(y), "Length of curves do not match" + x_str = ",".join(str(i) for i in x) + y_str = ",".join(str(i) for i in y) + dss.Text.Command(f"New {name} npts={npts} Yarray=({y_str}) Xarray=({x_str})") + self._xycurve_counter += 1 + return name + + def set_properties_to_inverter(self, inverter: str, inv_control: InverterControl): + """Modify a legacy InvControl object.""" + if inv_control.vvcontrol is not None: + vvc_curve = self.create_xy_curve( + inv_control.vvcontrol.voltage, inv_control.vvcontrol.reactive_response + ) + dss.Text.Command(f"{inverter}.vvc_curve1={vvc_curve.split('.')[1]}") + dss.Text.Command( + f"{inverter}.deltaQ_factor={inv_control.vvcontrol.deltaq_factor}" + ) + dss.Text.Command( + f"{inverter}.VarChangeTolerance={inv_control.vvcontrol.varchangetolerance}" + ) + dss.Text.Command( + f"{inverter}.VoltageChangeTolerance={inv_control.vvcontrol.voltagechangetolerance}" + ) + dss.Text.Command( + f"{inverter}.VV_RefReactivePower={inv_control.vvcontrol.vv_refreactivepower}" + ) + if inv_control.vwcontrol is not None: + vw_curve = self.create_xy_curve( + inv_control.vwcontrol.voltage, inv_control.vwcontrol.power_response, + ) + dss.Text.Command(f"{inverter}.voltwatt_curve={vw_curve.split('.')[1]}") + dss.Text.Command( + f"{inverter}.deltaP_factor={inv_control.vwcontrol.deltap_factor}" + ) + if inv_control.mode == InverterControlMode.voltvar_voltwatt: + dss.Text.Command(f"{inverter}.CombiMode = VV_VW") + else: + dss.Text.Command(f"{inverter}.Mode = {inv_control.mode.value}") - :return per unit voltages: - ''' - res_feeder_voltages = np.abs(self.get_voltages_complex()) - return res_feeder_voltages + def apply_inverter_control(self, inv_control: InverterControl): + """Apply inverter control to OpenDSS. - def get_voltages_complex(self): - ''' + Create InvControl if necessary and modifies it. - :return per unit voltages: - ''' - #voltages = np.array(self._circuit.AllBusVolts()) - #voltages = voltages[::2] + 1j * voltages[1::2] - #assert len(voltages) == len(self._AllNodeNames) - #return voltages + Parameters + ---------- + inv_control: InverterControl - _, name_voltage_dict = get_voltages(self._circuit) - # print('Publish load ' + str(temp_feeder_voltages.real[0])) - # set_load(dss, data_df['load'], current_time, loads) - # feeder_voltages = [0j] * len(self._AllNodeNames) - res_feeder_voltages = np.zeros((len(self._AllNodeNames)), dtype=np.complex_) - # print([name_voltage_dict.keys()][:5]) - temp_array_pu = np.zeros((len(self._AllNodeNames)), dtype=np.complex_) - for voltage_name in name_voltage_dict.keys(): + Returns + ------- + name of inverter with ``InvControl.`` prefix: str - self._circuit.SetActiveBus(voltage_name) - # temp1 = dss.Bus.puVmagAngle() - temp = dss.Bus.PuVoltage() - temp = complex(temp[0], temp[1]) - temp_array_pu[self._name_index_dict[voltage_name]] = temp - res_feeder_voltages[self._name_index_dict[voltage_name]] = name_voltage_dict[voltage_name] / (self._Vbase_allnode_dict[voltage_name]) + Warnings and Caveats + -------------------- + Using multiple pvsystems will issue a warning, since + more OpenDSS nodes can be created, which will break + more getter functions. - return res_feeder_voltages + New inverters cannot be assigned new lists. - def run_command(self, cmd): - dss.run_command(cmd) + This only works with the legacy InvControl settings, + volt var, volt watt, and volt-var volt-watt combined mode. + """ + if inv_control.pvsystem_list is None: + pvsystem_set = self._pvsystems + else: + pvsystem_set = set(inv_control.pvsystem_list) + inverter_set = set( + self._pvsystem_to_inverter[pvsystem] + for pvsystem in pvsystem_set + if pvsystem in self._pvsystem_to_inverter + ) + if len(inverter_set) == 1: + (inverter,) = inverter_set + else: + inverter = self.create_inverter(pvsystem_set) - def solve(self,hour,second): - dss.run_command(f'set mode=yearly loadmult=1 number=1 hour={hour} sec={second} stepsize={self._simulation_time_step} ') - dss.run_command('solve') + assert ( + self._inverter_to_pvsystems[inverter] == pvsystem_set + ), f"{self._inverter_to_pvsystems[inverter]} does not match {pvsystem_set} for {inverter}" + self.set_properties_to_inverter(inverter, inv_control) + return inverter diff --git a/LocalFeeder/README.md b/LocalFeeder/README.md index 2089f34..50c9c65 100644 --- a/LocalFeeder/README.md +++ b/LocalFeeder/README.md @@ -9,7 +9,5 @@ https://data.openei.org/s3_viewer?bucket=oedi-data-lake&prefix=SMART-DS%2Fv1.0%2 Or under the s3 link, ``` -oedi-data-lake/SMART-DS/v1.0/2017/SFO/P1U/scenarios/base_timeseries/opendss/p1uhs0_1247/p1uhs0_1247--p1udt942 +oedi-data-lake/SMART-DS/v1.0/2017/SFO/P1U/scenarios/base_timeseries/opendss/p1uhs0_1247/p1uhs0_1247--p1udt942 ``` - - diff --git a/LocalFeeder/component_definition.json b/LocalFeeder/component_definition.json index b9445cc..2c88254 100644 --- a/LocalFeeder/component_definition.json +++ b/LocalFeeder/component_definition.json @@ -8,13 +8,18 @@ {"type": "", "port_id": "number_of_timesteps"}, {"type": "", "port_id": "start_time_index"} ], - "dynamic_inputs": [], + "dynamic_inputs": [ + {"type": "CommandList", "port_id": "change_commands", "optional": true}, + {"type": "InverterControlList", "port_id": "inv_control", "optional": true} + ], "dynamic_outputs": [ {"type": "VoltagesMagnitude", "port_id": "voltages_magnitude"}, {"type": "VoltagesReal", "port_id": "voltages_real"}, {"type": "VoltagesImaginary", "port_id": "voltages_imag"}, {"type": "PowersReal", "port_id": "powers_real"}, {"type": "PowersImaginary", "port_id": "powers_imag"}, - {"type": "Topology", "port_id": "topology"} + {"type": "Topology", "port_id": "topology"}, + {"type": "Injection", "port_id": "injections"}, + {"type": "", "port_id": "load_y_matrix"} ] } diff --git a/LocalFeeder/dss_functions.py b/LocalFeeder/dss_functions.py index aaab04f..3b04219 100644 --- a/LocalFeeder/dss_functions.py +++ b/LocalFeeder/dss_functions.py @@ -1,160 +1,9 @@ -import csv +"""OpenDSS functions. Mutates global state, originally from GO-Solar project.""" import math -import numpy as np -import pandas as pd -from scipy import sparse as sparse -import cmath -from scipy.sparse import csc_matrix, save_npz - -# from simulator.simulator import node_number - -vmult = 0.001 -pos120 = complex (-0.5, 0.5 * cmath.sqrt(3.0)); # For Phase C -neg120 = complex (-0.5, -0.5 * cmath.sqrt(3.0)); # For Phase B -phase_shift = {'A':1,'B':neg120,'C':pos120,'s1':1,'s2':1} -lookup = {'A': '1', 'B': '2', 'C': '3', 'N':'4','S1':'1','S2':'2', 's1':'1','s2':'2', 's1\ns2': ['1','2'],'s2\ns1': ['2','1'],'': '1.2.3'} - - -def get_y_matrix_file(dss): - # dss.Circuit.SetActiveClass("Vsource") - # - # dss.Circuit.SetActiveElement("source") - # batchedit transformer..* wdg=2 tap=1 - # batchedit regcontrol..* enabled=false - # batchedit vsource..* enabled=false - # batchedit isource..* enabled=false - # batchedit load..* enabled=false - # batchedit generator..* enabled=false - # batchedit pvsystem..* enabled=false - # batchedit storage..* enabled=false - # solve - # export y triplet base_ysparse.csv - # export ynodelist base_nodelist.csv - # export summary base_summary.csv - dss.run_command('batchedit transformer..* wdg=2 tap=1') - dss.run_command('batchedit regcontrol..* enabled=false') - dss.run_command('batchedit vsource..* enabled=false') - dss.run_command('batchedit isource..* enabled=false') - dss.run_command('batchedit load..* enabled=false') - dss.run_command('batchedit generator..* enabled=false') - dss.run_command('batchedit pvsystem..* enabled=false') - dss.run_command('Batchedit Capacitor..* enabled=false') - dss.run_command('batchedit storage..* enabled=false') - dss.run_command('CalcVoltageBases') - dss.run_command('set maxiterations=20') - # solve - dss.run_command('solve') - dss.run_command('export y triplet base_ysparse.csv') - dss.run_command('export ynodelist base_nodelist.csv') - dss.run_command('export summary base_summary.csv') - Ysparse = csc_matrix(dss.YMatrix.getYsparse()) - save_npz('base_ysparse.npz', Ysparse) - - # dss.run_command('show Y') - # dss.run_command('solve mode=snap') - -def get_vnom2(dss): - circuit = dss.Circuit - AllNodeNames = circuit.AllNodeNames() - # This or DSSCircuit.AllBusVmagPu - Vnom = circuit.AllBusMagPu() # circuit.AllBusVmagPu() - - vmags = circuit.AllBusVMag() - Vnom2 = circuit.AllBusVolts() - test_Vnom2 = np.array([complex(Vnom2[i], Vnom2[i + 1]) for i in range(0, len(Vnom2), 2)]) - test_Vnom2_dict = {AllNodeNames[ii].upper(): test_Vnom2[ii] for ii in range(len(test_Vnom2))} - - test_vmag_volts_result = np.allclose(vmags, np.abs(test_Vnom2)) - logger.debug('test_vmag_volts_result', test_vmag_volts_result) - - AllNodeNamesY = circuit.YNodeOrder() - yv = circuit.YNodeVArray() - test_yv = np.array([complex(yv[i], yv[i + 1]) for i in range(0, len(yv), 2)]) - test_yv_dict = {AllNodeNamesY[ii]: test_yv[ii] for ii in range(len(test_yv))} - test_yv_result = np.allclose(vmags, np.abs(test_yv)) - logger.debug('test_yv_result', test_yv_result) - - logger.debug('Test dictionary') - for i in test_yv_dict.keys(): - if abs(abs(test_Vnom2_dict[i]) - abs(test_yv_dict[i])) > .0001: - logger.debug(i, abs(test_Vnom2_dict[i]), abs(test_yv_dict[i])) - # for t1, t2 in zip(np.abs(test_Vnom2), np.abs(test_yv)): - # np.testing.assert_array_almost_equal(np.abs(test_Vnom2), np.abs(test_yv),decimal=5) - - V = np.ones(len(Vnom) // 2, dtype=np.complex_) - for i in range(len(V)): - V[i] = complex(Vnom[2 * i], Vnom[2 * i + 1]) - vnom_dict = {AllNodeNames[ii].upper(): V[ii] for ii in range(len(V))} - return V, vnom_dict - -def get_vnom(dss): - dss.run_command('BatchEdit Load..* enabled=no') - dss.run_command('BatchEdit Generator..* enabled=no') - dss.run_command('solve mode=snap') - - circuit = dss.Circuit - AllNodeNames = circuit.AllNodeNames() - Vnom = circuit.AllBusVolts() - bases = dss.Settings.VoltageBases() - base = bases[-1] / math.sqrt(3) - V = np.ones(len(Vnom) // 2, dtype=np.complex_) - for i in range(len(V)): - V[i] = complex(Vnom[2 * i], Vnom[2 * i + 1]) / (base * 1000) - - vnom_dict = {AllNodeNames[ii].upper(): V[ii] for ii in range(len(V))} - dss.run_command('BatchEdit Load..* enabled=yes') - dss.run_command('BatchEdit Generator..* enabled=yes') - dss.run_command('solve mode=snap') - return V, vnom_dict - - -def snapshot_run(dss): - """ Function used for configuring the daily mode simulation run - """ - - dss.run_command('Batchedit Load..* enabled=yes') - dss.run_command('Batchedit Vsource..* enabled=yes') - dss.run_command('Batchedit Isource..* enabled=yes') - dss.run_command('Batchedit Generator..* enabled=yes') - dss.run_command('Batchedit PVsystem..* enabled=yes') - dss.run_command('Batchedit Capacitor..* enabled=yes') - dss.run_command('Batchedit Storage..* enabled=no') - dss.run_command('CalcVoltageBases') - dss.run_command('solve mode=snapshot') - - -def parse_Ymatrix(Ysparse, totalnode_number): - ''' - Read platform y maxtrix file. It has no '[', ']' or '=' - Handels slack nodes positions not at the top-left of the matrix - :param Ysparse: - :param slack_no: - :param totalnode_number: - :return: Y00,Y01,Y10,Y11,Y11_sparse,Y11_inv - ''' - - Ymatrix = sparse.lil_matrix((totalnode_number, totalnode_number), dtype=np.complex_) - - with open(Ysparse, 'r') as csvfile: - reader = csv.reader(csvfile, delimiter=',', ) - next(reader, None) # skip the headers - for row in reader: - row_value = int(row[0]) - column_value = int(row[1]) - r = row_value - 1 - c = column_value - 1 - g = float(row[2]) - b = float(row[3]) - Ymatrix[r, c] = complex(g, b) - Ymatrix[c, r] = Ymatrix[r, c] - - from scipy.sparse import load_npz - Ymatrix_new = load_npz('base_ysparse.npz') - return Ymatrix_new.tocoo() - def get_loads(dss, circuit): + """Get list of load dicts from OpenDSS circuit.""" data = [] load_flag = dss.Loads.First() dss.Circuit.SetActiveClass("Load") @@ -165,34 +14,34 @@ def get_loads(dss, circuit): "kV": load.kV(), "kW": load.kW(), "PF": load.PF(), - "Delta_conn": load.IsDelta() + "Delta_conn": load.IsDelta(), } - indexCktElement = circuit.SetActiveElement("Load.%s" % datum["name"]) + _ = circuit.SetActiveElement("Load.%s" % datum["name"]) cktElement = dss.CktElement bus = cktElement.BusNames()[0].split(".") - datum["kVar"] = float(datum["kW"]) / float(datum["PF"]) * math.sqrt(1 - float(datum["PF"]) * float(datum["PF"])) + datum["kVar"] = ( + float(datum["kW"]) + / float(datum["PF"]) + * math.sqrt(1 - float(datum["PF"]) * float(datum["PF"])) + ) datum["bus1"] = bus[0] datum["numPhases"] = len(bus[1:]) datum["phases"] = bus[1:] if not datum["numPhases"]: datum["numPhases"] = 3 - datum["phases"] = ['1', '2', '3'] - # if not datum["numPhases"] == [u'']: - # datum["numPhases"] = 3 - # datum['busphase'] = ['1', '2', '3'] - # else: - # datum['busphase'] = [lookup[phs] for phs in datum["phases"]] - # datum["numPhases"] = len(datum['busphase']) + datum["phases"] = ["1", "2", "3"] datum["voltageMag"] = cktElement.VoltagesMagAng()[0] datum["voltageAng"] = cktElement.VoltagesMagAng()[1] - datum["power"] = dss.CktElement.Powers()[0:2] + datum["power"] = dss.CktElement.Powers()[:2] data.append(datum) load_flag = dss.Loads.Next() return data -def get_pvSystems(dss): + +def get_pvsystems(dss): + """Get list of PVSystem dicts from OpenDSS circuit.""" data = [] PV_flag = dss.PVsystems.First() while PV_flag: @@ -207,67 +56,60 @@ def get_pvSystems(dss): NumPhase = dss.CktElement.NumPhases() bus = dss.CktElement.BusNames()[0] - #PVkV = dss.run_command('? ' + PVname + '.kV') #Not included in PVsystems commands for some reason - + # PVkV = dss.run_command('? ' + PVname + '.kV') + # Not included in PVsystems commands for some reason datum["name"] = PVname datum["bus"] = bus datum["phases"] = bus[1:] datum["Pmpp"] = PVpmpp datum["pf"] = PVpf - #datum["kV"] = PVkV datum["kW"] = PVkW datum["kVar"] = PVkvar datum["kVARated"] = PVkVARated datum["numPhase"] = NumPhase datum["numPhases"] = NumPhase - datum["power"] = dss.CktElement.Powers()[0:2*NumPhase] + datum["power"] = dss.CktElement.Powers()[: 2 * NumPhase] data.append(datum) PV_flag = dss.PVsystems.Next() return data -def get_Generator(dss): +def get_generators(dss): + """Get list of Generator dicts from OpenDSS circuit.""" data = [] gen_flag = dss.Generators.First() dss.Circuit.SetActiveClass("Generator") while gen_flag: - # gen = dss.Generators.Name() - # print(gen) - datum = {} - # GENname = dss.CktElement.Name() GENname = dss.Generators.Name() NumPhase = dss.CktElement.NumPhases() bus = dss.CktElement.BusNames()[0] - GENkVar = dss.Generators.kvar() GENkW = dss.Generators.kW() GENpf = dss.Generators.PF() GENkV = dss.Generators.kV() - datum["name"] = GENname - # datum["bus"] = bus - bus = bus.split('.') + bus = bus.split(".") if len(bus) == 1: - bus = bus + ['1', '2', '3'] - # else: - # bus = [bus[1]] - datum["bus"] = bus - datum["phases"] = bus[1:] - datum["name_bus"] = datum["name"] + '.' + bus[0] ## TOOO multiple phases - datum["kW"] = GENkW - datum["kVar"] = dss.Generators.kvar() - datum["pf"] = GENpf - datum["kV"] = GENkV - #datum["kVA"] = GENkVA - datum["numPhase"] = NumPhase - datum["numPhases"] = NumPhase - #datum["power"] = dss.CktElement.Powers()[0:2*NumPhase] + bus = bus + ["1", "2", "3"] + datum = { + "name": GENname, + "bus": bus, + "phases": bus[1:], + "name_bus": GENname + "." + bus[0], + "kW": GENkW, + "kVar": dss.Generators.kvar(), + "pf": GENpf, + "kV": GENkV, + "numPhase": NumPhase, + "numPhases": NumPhase, + } data.append(datum) gen_flag = dss.Generators.Next() return data def get_capacitors(dss): + """Get list of Capacitor dicts from OpenDSS circuit.""" data = [] cap_flag = dss.Capacitors.First() dss.Circuit.SetActiveClass("Capacitor") @@ -278,40 +120,26 @@ def get_capacitors(dss): bus = dss.CktElement.BusNames()[0] kvar = dss.Capacitors.kvar() datum["name"] = capname - temp = bus.split('.') + temp = bus.split(".") datum["busname"] = temp[0] datum["busphase"] = temp[1:] if not datum["busphase"]: - datum["busphase"] = ['1','2','3'] + datum["busphase"] = ["1", "2", "3"] datum["kVar"] = kvar datum["numPhases"] = NumPhase - datum["power"] = dss.CktElement.Powers()[0:2 * NumPhase] + datum["power"] = dss.CktElement.Powers()[: 2 * NumPhase] data.append(datum) cap_flag = dss.Capacitors.Next() return data def get_voltages(circuit): + """Get dict of names to voltages from OpenDSS circuit.""" temp_Vbus = circuit.YNodeVArray() AllNodeNames = circuit.YNodeOrder() node_number = len(AllNodeNames) - name_voltage_dict = {AllNodeNames[ii]:complex(temp_Vbus[ii * 2], temp_Vbus[ii * 2 + 1]) for ii in range(node_number)} - feeder_voltages = np.array([complex(temp_Vbus[ii * 2], temp_Vbus[ii * 2 + 1]) for ii in range(node_number)]) - # voltage_pu = list(map(lambda x: abs(x[0]) / x[1], zip(feeder_voltages, BASEKV))) - # print(feeder_voltages) - # print(BASEKV) - return feeder_voltages, name_voltage_dict - - -def get_load_sizes(dss, loads): - load_size_dict = {} - for load in loads: - bus_phases = load['bus1'] +'.' + '.'.join(load['phases']) - load_size_dict[load['name']] = {'numPhases':load['numPhases'], - 'bus':load['bus1'], - 'bus_phases': bus_phases, - 'phases':load['phases'], - 'kVar': load['kVar'], - 'kV': load['kV'], - 'kW': load['kW']} - return load_size_dict + name_voltage_dict = { + AllNodeNames[ii]: complex(temp_Vbus[ii * 2], temp_Vbus[ii * 2 + 1]) + for ii in range(node_number) + } + return name_voltage_dict diff --git a/LocalFeeder/mypy.ini b/LocalFeeder/mypy.ini new file mode 100644 index 0000000..e743215 --- /dev/null +++ b/LocalFeeder/mypy.ini @@ -0,0 +1,19 @@ +[mypy] +plugins = pydantic.mypy + +ignore_missing_imports = True +follow_imports = skip +warn_redundant_casts = True +warn_unused_ignores = True +disallow_any_generics = True +#check_untyped_defs = True # hard with opendss +no_implicit_reexport = True + +# for strict mypy: (this is the tricky one :-)) +#disallow_untyped_defs = True + +[pydantic-mypy] +init_forbid_extra = True +init_typed = True +warn_required_dynamic_aliases = True +warn_untyped_fields = True diff --git a/LocalFeeder/pytest.ini b/LocalFeeder/pytest.ini new file mode 100644 index 0000000..e610b3c --- /dev/null +++ b/LocalFeeder/pytest.ini @@ -0,0 +1,3 @@ +[pytest] +log_cli = True +log_cli_level = INFO diff --git a/LocalFeeder/sender_cosim.py b/LocalFeeder/sender_cosim.py index f720377..380104e 100644 --- a/LocalFeeder/sender_cosim.py +++ b/LocalFeeder/sender_cosim.py @@ -1,260 +1,448 @@ +"""HELICS wrapper for OpenDSS feeder simulation.""" +import json import logging +from dataclasses import dataclass +from datetime import datetime, timedelta +from typing import Any, Dict, List + import helics as h -import opendssdirect as dss -import pandas as pd -import json -from dss_functions import snapshot_run -from FeederSimulator import FeederSimulator, FeederConfig -from pydantic import BaseModel -from typing import List, Tuple import numpy as np -from datetime import datetime, timedelta -from gadal.gadal_types.data_types import ( - Complex,Topology,VoltagesReal,VoltagesImaginary,PowersReal, - PowersImaginary, AdmittanceMatrix, VoltagesMagnitude, - VoltagesAngle, Injection, AdmittanceSparse +import numpy.typing as npt +import xarray as xr +from oedisi.types.data_types import ( + AdmittanceMatrix, + AdmittanceSparse, + CommandList, + EquipmentNodeArray, + Injection, + InverterControlList, + MeasurementArray, + PowersImaginary, + PowersReal, + Topology, + VoltagesAngle, + VoltagesImaginary, + VoltagesMagnitude, + VoltagesReal, ) +from scipy.sparse import coo_matrix + +from FeederSimulator import FeederConfig, FeederSimulator logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) logger.setLevel(logging.DEBUG) -test_se = False -def numpy_to_y_matrix(array): - return [ - [(element.real, element.imag) for element in row] - for row in array - ] +def numpy_to_y_matrix(array: npt.NDArray[np.complex64]): + """Convert 2d numpy array to list of lists.""" + return [[(element.real, element.imag) for element in row] for row in array] -def sparse_to_admittance_sparse(array, unique_ids): + +def sparse_to_admittance_sparse(array: coo_matrix, unique_ids: List[str]): + """Convert coo sparse array to AdmittanceSparse type.""" return AdmittanceSparse( - from_equipment = [unique_ids[i] for i in array.row], - to_equipment = [unique_ids[i] for i in array.col], - admittance_list = [(data.real, data.imag) for data in array.data] + from_equipment=[unique_ids[i] for i in array.row], + to_equipment=[unique_ids[i] for i in array.col], + admittance_list=[(data.real, data.imag) for data in array.data], ) -def setup_sim(config: FeederConfig): - sim = FeederSimulator(config) - - snapshot_run(dss) - return sim - def get_true_phases(angle): - if np.abs(angle-0)<0.2: + """Round complex angles to predefined set of phases.""" + if np.abs(angle - 0) < 0.2: return 0 - elif np.abs(angle-np.pi/3)<0.2: - return np.pi/3 - elif np.abs(angle-2*np.pi/3)<0.2: - return 2*np.pi/3 - elif np.abs(angle-3*np.pi/3)<0.2: - return 3*np.pi/3 - elif np.abs(angle-(-np.pi/3))<0.2: - return -np.pi/3 - elif np.abs(angle-(-2*np.pi/3))<0.2: - return -2*np.pi/3 - elif np.abs(angle-(-3*np.pi/3))<0.2: - return -3*np.pi/3 + elif np.abs(angle - np.pi / 3) < 0.2: + return np.pi / 3 + elif np.abs(angle - 2 * np.pi / 3) < 0.2: + return 2 * np.pi / 3 + elif np.abs(angle - 3 * np.pi / 3) < 0.2: + return 3 * np.pi / 3 + elif np.abs(angle - (-np.pi / 3)) < 0.2: + return -np.pi / 3 + elif np.abs(angle - (-2 * np.pi / 3)) < 0.2: + return -2 * np.pi / 3 + elif np.abs(angle - (-3 * np.pi / 3)) < 0.2: + return -3 * np.pi / 3 else: logger.debug("error") -def go_cosim(sim, config: FeederConfig): +def xarray_to_dict(data): + """Convert xarray to dict with values and ids for JSON serialization.""" + coords = {key: list(data.coords[key].data) for key in data.coords.keys()} + return {"values": list(data.data), **coords} - deltat = 0.01 - fedinitstring = "--federates=1" - logger.info("Creating Federate Info") - fedinfo = h.helicsCreateFederateInfo() - h.helicsFederateInfoSetCoreName(fedinfo, config.name) - h.helicsFederateInfoSetCoreTypeFromString(fedinfo, "zmq") - h.helicsFederateInfoSetCoreInitString(fedinfo, fedinitstring) - h.helicsFederateInfoSetTimeProperty(fedinfo, h.helics_property_time_delta, deltat) - vfed = h.helicsCreateValueFederate(config.name, fedinfo) +def xarray_to_powers(data, **kwargs): + """Conveniently turn xarray into PowersReal and PowersImaginary.""" + powersreal = PowersReal(**xarray_to_dict(data.real), **kwargs) + powersimag = PowersImaginary(**xarray_to_dict(data.imag), **kwargs) + return powersreal, powersimag - pub_voltages_real = h.helicsFederateRegisterPublication(vfed, "voltages_real", h.HELICS_DATA_TYPE_STRING, "") - pub_voltages_imag = h.helicsFederateRegisterPublication(vfed, "voltages_imag", h.HELICS_DATA_TYPE_STRING, "") - pub_voltages_magnitude = h.helicsFederateRegisterPublication(vfed, "voltages_magnitude", h.HELICS_DATA_TYPE_STRING, "") - pub_powers_real = h.helicsFederateRegisterPublication(vfed, "powers_real", h.HELICS_DATA_TYPE_STRING, "") - pub_powers_imag = h.helicsFederateRegisterPublication(vfed, "powers_imag", h.HELICS_DATA_TYPE_STRING, "") - pub_topology = h.helicsFederateRegisterPublication(vfed, "topology", h.HELICS_DATA_TYPE_STRING, "") - h.helicsFederateEnterExecutingMode(vfed) +def concat_measurement_arrays(*ps: MeasurementArray): + """Concatenate list of measurements into one.""" + accuracy = None + if all((p.accuracy is not None for p in ps)): + accuracy = [e for p in ps for e in p.accuracy] - Y = sim.get_y_matrix() - #logger.debug("Eigenvalues and vectors") - #logger.debug(np.linalg.eig(Y.toarray())) - #y_matrix = numpy_to_y_matrix(Y.toarray()) - unique_ids = sim._AllNodeNames + assert all(ps[0].units == p.units for p in ps) - if config.use_sparse_admittance: - admittancematrix = sparse_to_admittance_sparse( - Y, - unique_ids + bad_data_threshold = None + if all((p.bad_data_threshold is not None for p in ps)): + bad_data_threshold = [e for p in ps for e in p.bad_data_threshold] + + assert all(ps[0].time == p.time for p in ps) + + if all((isinstance(p, EquipmentNodeArray) for p in ps)): + equipment_ids = [e for p in ps for e in p.equipment_ids] + + return ps[0].__class__( + values=[v for p in ps for v in p.values], + ids=[id for p in ps for id in p.ids], + equipment_ids=equipment_ids, + units=ps[0].units, + accuracy=accuracy, + bad_data_threshold=bad_data_threshold, + time=ps[0].time, ) else: - admittancematrix = AdmittanceMatrix( - admittance_matrix= numpy_to_y_matrix( - Y.toarray() - ), - ids=unique_ids + return ps[0].__class__( + values=[v for p in ps for v in p.values], + ids=[id for p in ps for id in p.ids], + units=ps[0].units, + accuracy=accuracy, + bad_data_threshold=bad_data_threshold, + time=ps[0].time, ) - logger.debug("_Vbase_allnode") - logger.debug(sim._Vbase_allnode) +def get_powers(PQ_load, PQ_PV, PQ_gen, PQ_cap): + """Turn xararys into PowersReal and PowersImaginary.""" + PQ_load_real, PQ_load_imag = xarray_to_powers(PQ_load) + PQ_PV_real, PQ_PV_imag = xarray_to_powers(PQ_PV) + PQ_gen_real, PQ_gen_imag = xarray_to_powers(PQ_gen) + PQ_cap_real, PQ_cap_imag = xarray_to_powers(PQ_cap) + + power_real = concat_measurement_arrays( + PQ_load_real, PQ_PV_real, PQ_gen_real, PQ_cap_real + ) + power_imag = concat_measurement_arrays( + PQ_load_imag, PQ_PV_imag, PQ_gen_imag, PQ_cap_imag + ) + return power_real, power_imag + + +@dataclass +class InitialData: + """Initial data from start of simulation.""" + + Y: Any + topology: Topology - slack_bus = [ - sim._AllNodeNames[i] for i in range( - sim._source_indexes[0], sim._source_indexes[-1] + 1 + +def get_initial_data(sim: FeederSimulator, config: FeederConfig): + """Get and calculate InitialData from simulation.""" + Y = sim.get_y_matrix() + unique_ids = sim._AllNodeNames + + if config.use_sparse_admittance: + admittancematrix = sparse_to_admittance_sparse(Y, unique_ids) + else: + admittancematrix = AdmittanceMatrix( + admittance_matrix=numpy_to_y_matrix(Y.toarray()), ids=unique_ids ) + + slack_ids = [ + sim._AllNodeNames[i] + for i in range(sim._source_indexes[0], sim._source_indexes[-1] + 1) ] - unique_ids = sim._AllNodeNames - snapshot_run(sim) - - logger.debug("slack_bus") - logger.debug(slack_bus) - logger.debug("unique_ids") - logger.debug(unique_ids) - - - all_PQs = {} - all_PQs['load'] = sim.get_PQs_load(static=True) # Return type is PQ_values, PQ_names, PQ_types all with same size - all_PQs['pv'] = sim.get_PQs_pv(static=True) - all_PQs['gen'] = sim.get_PQs_gen(static=True) - all_PQs['cap'] = sim.get_PQs_cap(static=True) - - PQ_real = [] - PQ_imaginary = [] - PQ_names = [] - PQ_types = [] - for i in range(len(all_PQs['load'][0])): - for key in all_PQs: - if all_PQs[key][1][i] !='': - PQ_real.append(-1*all_PQs[key][0][i].real) # injections are negative singe PQ values are positive - PQ_imaginary.append(-1*all_PQs[key][0][i].imag) # injections are negative singe PQ values are positive - PQ_names.append(all_PQs[key][1][i]) - PQ_types.append(all_PQs[key][2][i]) - power_real = PowersReal(ids = PQ_names, values = PQ_real, equipment_type = PQ_types) - power_imaginary = PowersImaginary(ids = PQ_names, values = PQ_imaginary, equipment_type = PQ_types) - injections = Injection(power_real=power_real, power_imaginary = power_imaginary) - - sim.solve(0,0) - feeder_voltages = sim.get_voltages_actual() - phases = list(map(get_true_phases, np.angle(feeder_voltages))) - base_voltages = list(sim._Vbase_allnode) + base_voltages = sim.get_base_voltages() base_voltagemagnitude = VoltagesMagnitude( - values = [abs(i) for i in base_voltages], - ids = unique_ids + values=list(np.abs(base_voltages).data), ids=list(base_voltages.ids.data) ) - base_voltageangle = VoltagesAngle( - values = phases, - ids = unique_ids - ) + # We have to do snapshot run so we can re-enable things properly. + # Technically we don't have to solve. + sim.snapshot_run() + PQ_load = sim.get_PQs_load(static=True) + PQ_PV = sim.get_PQs_pv(static=True) + PQ_gen = sim.get_PQs_gen(static=True) + PQ_cap = sim.get_PQs_cap(static=True) + + sim.solve(0, 0) + power_real, power_imaginary = get_powers(-PQ_load, -PQ_PV, -PQ_gen, -PQ_cap) + injections = Injection(power_real=power_real, power_imaginary=power_imaginary) + + feeder_voltages = sim.get_voltages_actual() + feeder_angles: npt.NDArray[np.float64] = np.angle(feeder_voltages.data) + phases = list(map(get_true_phases, feeder_angles)) + base_voltageangle = VoltagesAngle(values=phases, ids=list(feeder_voltages.ids.data)) - topology = Topology( admittance=admittancematrix, base_voltage_angles=base_voltageangle, injections=injections, base_voltage_magnitudes=base_voltagemagnitude, - slack_bus=slack_bus, + slack_bus=slack_ids, ) + return InitialData(Y=Y, topology=topology) + + +def agg_to_ids(x: xr.core.dataarray.DataArray, ids): + """Aggregate xarray to ids. Specialized to equipment node arrays.""" + target = xr.zeros_like(ids, dtype=np.float64) + if x.shape == (0,): + return target + + _, x_grouped = xr.align(ids, x.groupby("ids").sum(), join="left", fill_value=0.0) + return x_grouped + + +@dataclass +class CurrentData: + """Current data at time t. ``arr.ids`` gives bus ids.""" + + feeder_voltages: xr.core.dataarray.DataArray + PQ_injections_all: xr.core.dataarray.DataArray + calculated_power: xr.core.dataarray.DataArray + injections: Injection + load_y_matrix: Any + + +def get_current_data(sim: FeederSimulator, Y): + """Construct current data from simulation after having solved.""" + feeder_voltages = sim.get_voltages_actual() + PQ_load = sim.get_PQs_load(static=False) + PQ_PV = sim.get_PQs_pv(static=False) + PQ_gen = sim.get_PQs_gen(static=False) + PQ_cap = sim.get_PQs_cap(static=False) + + # Assumes everything is controllable! + power_real, power_imaginary = get_powers(-PQ_load, -PQ_PV, -PQ_gen, -PQ_cap) + injections = Injection(power_real=power_real, power_imaginary=power_imaginary) + + ids = xr.DataArray(sim._AllNodeNames, coords={ + "ids": sim._AllNodeNames, + }) + PQ_injections_all = ( + agg_to_ids(PQ_load, ids) + + agg_to_ids(PQ_PV, ids) + + agg_to_ids(PQ_gen, ids) + + agg_to_ids(PQ_cap, ids) + ) + + PQ_injections_all = PQ_injections_all.assign_coords(equipment_ids=('ids', list(map(lambda x: x.split(".")[0], sim._AllNodeNames)))) + calculated_power = ( + feeder_voltages * (Y.conjugate() @ feeder_voltages.conjugate()) / 1000 + ) + + PQ_injections_all[sim._source_indexes] = -calculated_power[sim._source_indexes] + + Y_load = sim.get_load_y_matrix() + return CurrentData( + feeder_voltages=feeder_voltages, + PQ_injections_all=PQ_injections_all, + calculated_power=calculated_power, + injections=injections, + load_y_matrix=Y_load, + ) + + +def where_power_unbalanced(PQ_injections_all, calculated_power, tol=1): + """Find errors where PQ_injectinos does not match calculated power.""" + errors = PQ_injections_all + calculated_power + (indices,) = np.where(np.abs(errors) > tol) + return errors.ids[indices] + + +def go_cosim(sim: FeederSimulator, config: FeederConfig, input_mapping: Dict[str, str]): + """Run HELICS federate with FeederSimulation. + + TODO: Maybe this should be a class or a coroutine or something cleaner. + There are many options. + """ + deltat = 0.01 + fedinitstring = "--federates=1" + + logger.info("Creating Federate Info") + fedinfo = h.helicsCreateFederateInfo() + h.helicsFederateInfoSetCoreName(fedinfo, config.name) + h.helicsFederateInfoSetCoreTypeFromString(fedinfo, "zmq") + h.helicsFederateInfoSetCoreInitString(fedinfo, fedinitstring) + h.helicsFederateInfoSetTimeProperty(fedinfo, h.helics_property_time_delta, deltat) + vfed = h.helicsCreateValueFederate(config.name, fedinfo) + + pub_voltages_real = h.helicsFederateRegisterPublication( + vfed, "voltages_real", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_voltages_imag = h.helicsFederateRegisterPublication( + vfed, "voltages_imag", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_voltages_magnitude = h.helicsFederateRegisterPublication( + vfed, "voltages_magnitude", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_powers_real = h.helicsFederateRegisterPublication( + vfed, "powers_real", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_powers_imag = h.helicsFederateRegisterPublication( + vfed, "powers_imag", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_topology = h.helicsFederateRegisterPublication( + vfed, "topology", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_injections = h.helicsFederateRegisterPublication( + vfed, "injections", h.HELICS_DATA_TYPE_STRING, "" + ) + pub_load_y_matrix = h.helicsFederateRegisterPublication( + vfed, "load_y_matrix", h.HELICS_DATA_TYPE_STRING, "" + ) + + command_set_key = ( + "unused/change_commands" + if "change_commands" not in input_mapping + else input_mapping["change_commands"] + ) + sub_command_set = vfed.register_subscription(command_set_key, "") + sub_command_set.set_default("[]") + sub_command_set.option["CONNECTION_OPTIONAL"] = 1 + + inv_control_key = ( + "unused/inv_control" + if "" not in input_mapping + else input_mapping["inv_control"] + ) + sub_invcontrol = vfed.register_subscription(inv_control_key, "") + sub_invcontrol.set_default("[]") + sub_invcontrol.option["CONNECTION_OPTIONAL"] = 1 + + h.helicsFederateEnterExecutingMode(vfed) + initial_data = get_initial_data(sim, config) logger.info("Sending topology and saving to topology.json") with open(config.topology_output, "w") as f: - f.write(topology.json()) - pub_topology.publish(topology.json()) - + f.write(initial_data.topology.json()) + pub_topology.publish(initial_data.topology.json()) granted_time = -1 - current_hour = 0 - current_second = 0 - current_index = config.start_time_index for request_time in range(0, int(config.number_of_timesteps)): - while granted_time < request_time: - granted_time = h.helicsFederateRequestTime(vfed, request_time) - - logger.info('start time: '+str(datetime.now())) - current_index+=1 - current_timestamp = datetime.strptime(config.start_date, '%Y-%m-%d %H:%M:%S') + timedelta(seconds = current_index*config.run_freq_sec) - current_second+=config.run_freq_sec - if current_second >=60*60: - current_second = 0 - current_hour+=1 - logger.info(f'Get Voltages and PQs at {current_index} {granted_time} {request_time}') - - sim.solve(current_hour,current_second) - - feeder_voltages = sim.get_voltages_actual() - all_PQs = {} - all_PQs['load'] = sim.get_PQs_load(static=False) # Return type is PQ_values, PQ_names, PQ_types all with same size - all_PQs['pv'] = sim.get_PQs_pv(static=False) - all_PQs['gen'] = sim.get_PQs_gen(static=False) - all_PQs['cap'] = sim.get_PQs_cap(static=False) - - PQ_injections_all = all_PQs['load'][0]+ all_PQs['pv'][0] + all_PQs['gen'][0]+ all_PQs['cap'][0] - - logger.debug("Feeder Voltages") - logger.debug(feeder_voltages) - logger.debug("PQ") - logger.debug(PQ_injections_all) - logger.debug("Calculated Power") - Cal_power = feeder_voltages * (Y.conjugate() @ feeder_voltages.conjugate()) / 1000 - errors = PQ_injections_all + Cal_power - sort_errors = np.sort(np.abs(errors)) - logger.debug("errors") - logger.debug(errors) - if np.any(sort_errors[:-3] > 1): - raise ValueError('Power balance does not hold') - PQ_injections_all[sim._source_indexes] = -Cal_power[sim._source_indexes] - power_balance = (feeder_voltages * (Y.conjugate() @ feeder_voltages.conjugate()) / 1000) - logger.debug(power_balance) - indices, = np.nonzero(np.abs(errors) > 1) - logger.debug("Indices with error > 1") - logger.debug(indices) - logger.debug([sim._AllNodeNames[i] for i in indices]) - logger.debug("Power, Voltages, and Calculated Power at Indices") - logger.debug(PQ_injections_all[indices]) - logger.debug(feeder_voltages[indices]) - logger.debug(power_balance[indices]) - - - phases = list(map(get_true_phases, np.angle(feeder_voltages))) - - base_voltageangle = VoltagesAngle( - values = phases, - ids = unique_ids + granted_time = h.helicsFederateRequestTime(vfed, request_time) + + current_index = int(granted_time) # floors + current_timestamp = datetime.strptime( + config.start_date, "%Y-%m-%d %H:%M:%S" + ) + timedelta(seconds=granted_time * config.run_freq_sec) + floored_timestamp = datetime.strptime( + config.start_date, "%Y-%m-%d %H:%M:%S" + ) + timedelta(seconds=current_index * config.run_freq_sec) + + change_obj_cmds = CommandList.parse_obj(sub_command_set.json) + sim.change_obj(change_obj_cmds.__root__) + + inverter_controls = InverterControlList.parse_obj(sub_invcontrol.json) + for inv_control in inverter_controls.__root__: + sim.apply_inverter_control(inv_control) + + logger.info( + f"Solve at hour {floored_timestamp.hour} second " + f"{60*floored_timestamp.minute + floored_timestamp.second}" ) - logger.debug('Publish load ' + str(feeder_voltages.real[0])) - voltage_magnitudes = np.abs(feeder_voltages.real + 1j* feeder_voltages.imag) - pub_voltages_magnitude.publish(VoltagesMagnitude(values=list(voltage_magnitudes), ids=sim._AllNodeNames, time = current_timestamp).json()) - pub_voltages_real.publish(VoltagesReal(values=list(feeder_voltages.real), ids=sim._AllNodeNames, time = current_timestamp).json()) - pub_voltages_imag.publish(VoltagesImaginary(values=list(feeder_voltages.imag), ids=sim._AllNodeNames, time = current_timestamp).json()) - pub_powers_real.publish(PowersReal(values=list(PQ_injections_all.real), ids=sim._AllNodeNames, time = current_timestamp).json()) - pub_powers_imag.publish(PowersImaginary(values=list(PQ_injections_all.imag), ids=sim._AllNodeNames, time = current_timestamp).json()) + sim.snapshot_run() + sim.solve( + floored_timestamp.hour, + 60 * floored_timestamp.minute + floored_timestamp.second, + ) - logger.info('end time: '+str(datetime.now())) + current_data = get_current_data(sim, initial_data.Y) + + bad_bus_names = where_power_unbalanced( + current_data.PQ_injections_all, current_data.calculated_power + ) + if len(bad_bus_names) > 0: + raise ValueError( + f""" + Bad buses at {bad_bus_names.data} + + OpenDSS PQ + {current_data.PQ_injections_all.loc[bad_bus_names]} + + PowerBalance PQ + {current_data.calculated_power.loc[bad_bus_names]} + """ + ) + + logger.debug( + f"Publish load {current_data.feeder_voltages.ids.data[0]} " + f"{current_data.feeder_voltages.data[0]}" + ) + voltage_magnitudes = np.abs(current_data.feeder_voltages) + pub_voltages_magnitude.publish( + VoltagesMagnitude( + **xarray_to_dict(voltage_magnitudes), time=current_timestamp, + ).json() + ) + pub_voltages_real.publish( + VoltagesReal( + **xarray_to_dict(current_data.feeder_voltages.real), + time=current_timestamp, + ).json() + ) + pub_voltages_imag.publish( + VoltagesImaginary( + **xarray_to_dict(current_data.feeder_voltages.imag), + time=current_timestamp, + ).json() + ) + pub_powers_real.publish( + PowersReal( + **xarray_to_dict(current_data.PQ_injections_all.real), + time=current_timestamp, + ).json() + ) + pub_powers_imag.publish( + PowersImaginary( + **xarray_to_dict(current_data.PQ_injections_all.imag), + time=current_timestamp, + ).json() + ) + pub_injections.publish(current_data.injections.json()) + + if config.use_sparse_admittance: + pub_load_y_matrix.publish( + sparse_to_admittance_sparse( + current_data.load_y_matrix, sim._AllNodeNames + ).json() + ) + else: + pub_load_y_matrix.publish( + AdmittanceMatrix( + admittance_matrix=numpy_to_y_matrix( + current_data.load_y_matrix.toarray() + ), + ids=sim._AllNodeNames + ).json() + ) + + logger.info("end time: " + str(datetime.now())) h.helicsFederateDisconnect(vfed) h.helicsFederateFree(vfed) h.helicsCloseLibrary() -class FeederCosimConfig(BaseModel): - feeder_config: FeederConfig - - def run(): - with open('static_inputs.json') as f: + """Load static_inputs and input_mapping and run JSON.""" + with open("static_inputs.json") as f: parameters = json.load(f) + with open("input_mapping.json") as f: + input_mapping = json.load(f) config = FeederConfig(**parameters) - sim = setup_sim(config) - go_cosim(sim, config) + sim = FeederSimulator(config) + go_cosim(sim, config, input_mapping) + -if __name__ == '__main__': +if __name__ == "__main__": run() diff --git a/LocalFeeder/tests/test_feeder.py b/LocalFeeder/tests/test_feeder.py new file mode 100644 index 0000000..f26e8e2 --- /dev/null +++ b/LocalFeeder/tests/test_feeder.py @@ -0,0 +1,594 @@ +import logging +import os +import sys + +import numpy as np +import pandas as pd +import plotille +import pytest +import xarray as xr +from oedisi.types.data_types import ( + EquipmentNodeArray, + InverterControl, + InverterControlMode, + VVControl, + VWControl, +) + +sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +import FeederSimulator +import sender_cosim +from sender_cosim import agg_to_ids + + +@pytest.fixture(scope="session", autouse=True) +def init_federate_simulation(): + federate_config = FeederSimulator.FeederConfig( + **{ + "use_smartds": False, + "profile_location": "gadal_ieee123/profiles", + "opendss_location": "gadal_ieee123/qsts", + "sensor_location": "gadal_ieee123/sensors.json", + "start_date": "2017-01-01 00:00:00", + "number_of_timesteps": 96, + "run_freq_sec": 900, + "topology_output": "../../outputs/topology.json", + "name": "feeder", + } + ) + logging.info("Initializating simulation data") + _ = FeederSimulator.FeederSimulator(federate_config) + return + + +@pytest.fixture() +def federate_config(): + return FeederSimulator.FeederConfig( + **{ + "use_smartds": False, + "profile_location": "gadal_ieee123/profiles", + "opendss_location": "gadal_ieee123/qsts", + "sensor_location": "gadal_ieee123/sensors.json", + "existing_feeder_file": "opendss/master.dss", + "start_date": "2017-01-01 00:00:00", + "number_of_timesteps": 96, + "run_freq_sec": 900, + "topology_output": "../../outputs/topology.json", + "name": "feeder", + } + ) + + +def plot_y_matrix(Y): + Y_max = np.max(np.abs(Y)) + + width = Y.shape[0] // 2 + height = Y.shape[1] // 4 + cvs = plotille.Canvas(width, height, mode="rgb") + arrayY = np.log(np.abs(Y.toarray())[: (4 * height), : (2 * width)] + 1) + flatY = arrayY.reshape(arrayY.size) + max_flatY = flatY.max() + cvs.braille_image([int(255 * np.abs(y) / max_flatY) for y in flatY]) + print(f"Y max: {Y_max}") + print(cvs.plot()) + + +def test_ordering(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + Y = sim.get_y_matrix() + + with pytest.raises(AssertionError, match=".*DISABLED_RUN.*"): + sim.solve(0, 0) + with pytest.raises(AssertionError, match=".*DISABLED_RUN.*"): + _ = sim.get_voltages_actual() + + sim.snapshot_run() + base_voltages = sim.get_voltages_snapshot() + assert np.max(base_voltages) > 100 + Y2 = sim.get_y_matrix() + assert np.max(np.abs(Y - Y2)) < 0.001 + + sim.snapshot_run() + sim.solve(0, 0) + feeder_voltages = sim.get_voltages_actual() + assert np.max(feeder_voltages) > 100 + with pytest.raises(AssertionError, match=".*SOLVE_AT_TIME.*"): + _ = sim.get_voltages_snapshot() + + sim.solve(10, 0) + Y_load = sim.get_load_y_matrix() + plot_y_matrix(Y) + plot_y_matrix(Y_load - Y) + + +def rtheta_to_xy(r, theta): + x = np.array(r * np.cos(theta)) + y = np.array(r * np.sin(theta)) + return x, y + + +def test_voltages(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + startup(sim) + + # Get Voltages + base = sim.get_base_voltages() + sim.disabled_run() + sim.initial_disabled_solve() + disabled_solve = sim.get_disabled_solve_voltages() + sim.snapshot_run() + snapshot = sim.get_voltages_snapshot() + sim.solve(0, 0) + actuals = sim.get_voltages_actual() + + # Plot magnitudes + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + # fig.set_x_limits(min_=-3, max_=3) + fig.set_y_limits(min_=0, max_=2) + fig.color_mode = "byte" + fig.plot( + range(len(base)), np.abs(disabled_solve / base).data, lc=50, label="Disabled" + ) + fig.plot(range(len(base)), np.abs(snapshot / base).data, lc=75, label="Snapshot") + fig.plot(range(len(base)), np.abs(actuals / base).data, lc=100, label="Actuals") + print("\n" + fig.show(legend=True)) + + # Plot magnitudes + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + # fig.set_x_limits(min_=-3, max_=3) + fig.set_y_limits(min_=0.9, max_=1.1) + fig.color_mode = "byte" + fig.plot( + range(len(base)), np.abs(disabled_solve / base).data, lc=50, label="Disabled" + ) + fig.plot(range(len(base)), np.abs(snapshot / base).data, lc=75, label="Snapshot") + fig.plot(range(len(base)), np.abs(actuals / base).data, lc=100, label="Actuals") + print("\n" + fig.show(legend=True)) + + # Plot angles better + r = xr.DataArray(range(len(base)), {"ids": base.ids.data}) + 100 + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + fig.set_x_limits(min_=-400, max_=400) + fig.set_y_limits(min_=-400, max_=400) + fig.color_mode = "byte" + fig.scatter(*rtheta_to_xy(r, np.angle(base)), lc=50, label="Disabled") + fig.scatter(*rtheta_to_xy(r, np.angle(snapshot)), lc=75, label="Snapshot") + fig.scatter(*rtheta_to_xy(r, np.angle(actuals)), lc=100, label="Actuals") + print("\n" + fig.show(legend=True)) + + # Plot angles better + r = xr.DataArray(range(len(base)), {"ids": base.ids.data}) + 100 + fig = plotille.Figure() + fig.width = 90 + fig.height = 30 + fig.set_x_limits(min_=0, max_=400) + fig.set_y_limits(min_=-50, max_=50) + fig.color_mode = "byte" + fig.scatter(*rtheta_to_xy(r, np.angle(base)), lc=50, label="Disabled") + fig.scatter(*rtheta_to_xy(r, np.angle(snapshot)), lc=75, label="Snapshot") + fig.scatter(*rtheta_to_xy(r, np.angle(actuals)), lc=100, label="Actuals") + print("\n" + fig.show(legend=True)) + + +def test_xarray_translation(): + x = xr.DataArray( + [0 + 1j, 1 + 2j, 2 + 3j, 3 + 4j], + dims=("ids",), + coords={ + "ids": ("ids", ["1", "1", "2", "3"]), + "equipment_ids": ("ids", ["0", "1", "2", "3"]), + }, + ) + powerreal, powerimag = sender_cosim.xarray_to_powers(x) + assert powerreal.ids == ["1", "1", "2", "3"] + assert powerreal.values == [0.0, 1.0, 2.0, 3.0] + assert powerimag.values == [1.0, 2.0, 3.0, 4.0] + + +def test_simulation(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + startup(sim) + Y = sim.get_y_matrix() + plot_y_matrix(Y) + sim.snapshot_run() # You need to re-enable! + getting_and_concatentating_data(sim) + initial_data(sim, federate_config) + sim.snapshot_run() + simulation_middle(sim, Y) + + +def startup(sim): + assert sim._feeder_file is not None + assert sim._AllNodeNames is not None + assert sim._circuit is not None + assert sim._AllNodeNames is not None + + +def getting_and_concatentating_data(sim): + PQ_load = sim.get_PQs_load(static=True) + PQ_PV = sim.get_PQs_pv(static=True) + PQ_gen = sim.get_PQs_gen(static=True) + PQ_cap = sim.get_PQs_cap(static=True) + + pv_real, pv_imag = sender_cosim.xarray_to_powers(PQ_PV) + gen_real, gen_imag = sender_cosim.xarray_to_powers(PQ_gen) + test_real = sender_cosim.concat_measurement_arrays(pv_real, gen_real) + assert test_real.values[5] == PQ_PV.data[5].real + assert test_real.ids[5] == PQ_PV.ids.data[5] + + ids = xr.DataArray(sim._AllNodeNames, coords={"ids": sim._AllNodeNames}) + PQ_injections_all = ( + agg_to_ids(PQ_load, ids) + + agg_to_ids(PQ_PV, ids) + + agg_to_ids(PQ_gen, ids) + + agg_to_ids(PQ_cap, ids) + ) + assert sorted(list(PQ_injections_all.ids.data)) == sorted(sim._AllNodeNames) + + +def initial_data(sim, federate_config): + initial_data = sender_cosim.get_initial_data(sim, federate_config) + # Plot magnitudes + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + # fig.set_x_limits(min_=-3, max_=3) + fig.set_y_limits(min_=0, max_=3000) + fig.color_mode = "byte" + fig.plot( + range(len(initial_data.topology.base_voltage_magnitudes.ids)), + np.sort(initial_data.topology.base_voltage_magnitudes.values), + lc=100, + label="Base", + ) + print("\n" + fig.show(legend=True)) + + # Plot angles better + r = np.array(range(len(initial_data.topology.base_voltage_magnitudes.ids))) + 100 + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + fig.set_x_limits(min_=-400, max_=400) + fig.set_y_limits(min_=-400, max_=400) + fig.color_mode = "byte" + fig.scatter( + *rtheta_to_xy(r, initial_data.topology.base_voltage_angles.values), + lc=50, + label="Base", + ) + print("\n" + fig.show(legend=True)) + + df = pd.DataFrame( + { + "voltages": initial_data.topology.base_voltage_magnitudes.values, + "angles": initial_data.topology.base_voltage_angles.values, + } + ) + logging.info(df.describe()) + df = pd.DataFrame( + {"injections": initial_data.topology.injections.power_real.values} + ) + logging.info(df.describe()) + assert initial_data is not None + + +def plot_complex_array(data, label="Voltages"): + print("Index vs Complex Magnitude") + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + # fig.set_x_limits(min_=-3, max_=3) + fig.set_y_limits(min_=0, max_=np.max(np.abs(data)) * 1.1) + fig.color_mode = "byte" + fig.plot(range(len(data)), np.abs(data), lc=100, label=label) + print("\n" + fig.show(legend=True)) + + print("Index vs Complex Angle (radial plot)") + r = np.array(range(len(data))) + 100 + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + fig.set_x_limits(min_=-400, max_=400) + fig.set_y_limits(min_=-400, max_=400) + fig.color_mode = "byte" + fig.scatter(*rtheta_to_xy(r, np.angle(data)), lc=50, label=label) + print("\n" + fig.show(legend=True)) + + print("Log magnitude with angle (radial plot)") + r = np.log(np.array(np.abs(data)) + 1.0) + fig = plotille.Figure() + fig.width = 60 + fig.height = 30 + fig.set_x_limits(min_=-10, max_=10) + fig.set_y_limits(min_=-10, max_=10) + fig.color_mode = "byte" + fig.scatter(*rtheta_to_xy(r, np.angle(data)), lc=50, label=label) + print("\n" + fig.show(legend=True)) + + +def simulation_middle(sim, Y): + # this one may need to be properly ordered + logging.info(f"Current directory : {os.getcwd()}") + sim.solve(0, 0) + + current_data = sender_cosim.get_current_data(sim, Y) + assert len(current_data.injections.power_real.values) == len( + current_data.injections.power_real.ids + ) + + current_data_again = sender_cosim.get_current_data(sim, Y) + assert np.allclose( + current_data_again.feeder_voltages, + current_data.feeder_voltages + ) + + assert '113' in current_data.PQ_injections_all.equipment_ids.data + df = pd.DataFrame( + { + "p": current_data.PQ_injections_all.real, + "q": current_data.PQ_injections_all.imag, + "voltages": np.abs(current_data.feeder_voltages), + "phases": np.angle(current_data.feeder_voltages), + } + ) + logging.info(df.describe()) + print("Feeder Voltages") + plot_complex_array(current_data.feeder_voltages.data, label="Feeder Voltages") + assert np.max(current_data.feeder_voltages) > 50 + + plot_complex_array(current_data.PQ_injections_all.data, label="PQ injection") + + diff = np.abs( + current_data.PQ_injections_all - (-current_data.calculated_power) + ) * np.exp( + 1j + * ( + np.angle(current_data.PQ_injections_all) + - np.angle(-current_data.calculated_power) + ) + ) + + plot_complex_array(diff.data, label="Calculated - Injected") + + bad_bus_names = sender_cosim.where_power_unbalanced( + current_data.PQ_injections_all, current_data.calculated_power + ) + assert len(bad_bus_names) == 0 + + +def equipment_indices_on_equipment_node_array( + equipment_node_array: EquipmentNodeArray, equipment_type: str +): + return map( + lambda iv: iv[0], + filter( + lambda iv: iv[1].startswith(equipment_type), + enumerate(equipment_node_array.equipment_ids), + ), + ) + + +def test_controls(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + Y = sim.get_y_matrix() + sim.snapshot_run() # Needed to bring out of disabled state + sim.solve(8, 0) + sim.just_solve() + current_data = sender_cosim.get_current_data(sim, Y) + assert len(current_data.injections.power_real.values) == len( + current_data.injections.power_real.ids + ) + + bad_bus_names = sender_cosim.where_power_unbalanced( + current_data.PQ_injections_all, current_data.calculated_power + ) + assert len(bad_bus_names) == 0 + + # Find first with equipment type = PVSystem + power_real = current_data.injections.power_real + pv_system_indices = list( + equipment_indices_on_equipment_node_array(power_real, "PVSystem") + ) + max_index = np.argmax(np.abs([power_real.values[i] for i in pv_system_indices])) + pv_system_index = pv_system_indices[max_index] + pv_system_index = None + for i in range(len(power_real.ids)): + if power_real.ids[i] == "113.1" and power_real.equipment_ids[i].startswith( + "PVSystem" + ): + pv_system_index = i + break + assert pv_system_index is not None + + print( + f"{power_real.ids[pv_system_index]} {power_real.equipment_ids[pv_system_index]}" + ) + # Try setting current power to half of that. + assert abs(power_real.values[pv_system_index]) > 0.01 + sim.change_obj( + [ + FeederSimulator.Command( + obj_name="PVSystem.113", + obj_property="%Pmpp", + val=5, # power_real.values[pv_system_index] / 2, + ) + ] + ) + # dss.Text.Command("PVsystem.113.PF=0.01") + # Check properties in AllPropertyNames in CktElement or just Element + # Solve and observe current power. It should change. + sim.just_solve() + new_data = sender_cosim.get_current_data(sim, Y) + new_power_real = new_data.injections.power_real + print(power_real.values[pv_system_index]) + print(power_real.ids[pv_system_index]) + print(new_power_real.values[pv_system_index]) + print(new_power_real.ids[pv_system_index]) + assert ( + np.abs( + new_power_real.values[pv_system_index] - power_real.values[pv_system_index] + ) + > 1 + ) + (bad_indices,) = np.where( + np.abs(np.array(new_power_real.values) - np.array(power_real.values)) > 1 + ) + + for i in bad_indices: + print( + f"Old: {power_real.equipment_ids[i]} {power_real.ids[i]} " + f"{power_real.values[i]}" + ) + print( + f"New: {new_power_real.equipment_ids[i]} {new_power_real.ids[i]} " + f"{new_power_real.values[i]}" + ) + + assert bad_indices == [pv_system_index] + # Run another time step. What happens? + sim.solve(9, 0) + next_data = sender_cosim.get_current_data(sim, Y) + next_power_real = next_data.injections.power_real + print(f"8,0: {power_real.values[pv_system_index]}") + print(f"8,0: {new_power_real.values[pv_system_index]}") + print(f"9,0: {next_power_real.values[pv_system_index]}") + + +def test_inv_control(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + test_inverter_data = InverterControl( + pvsystem_list=["PVSystem.113"], + vvcontrol=VVControl( + deltaq_factor=0.7, + voltage=[0.5, 0.945, 0.97, 0.99, 1.045, 1.50], + reactive_response=[1, 1, 1, 1, 1, 1], + ), + ) + + _ = sim.get_y_matrix() + sim.snapshot_run() # Needed to bring out of disabled state + sim.solve(8, 0) + + old_voltages = sim.get_voltages_actual() + + sim.apply_inverter_control(test_inverter_data) + sim.just_solve() + + new_voltages = sim.get_voltages_actual() + + assert float( + np.sum(np.abs(new_voltages.loc["113.1"] - old_voltages.loc["113.1"])) + ) > 0.01 * float(np.abs(old_voltages.loc["113.1"])) + + sim.apply_inverter_control(test_inverter_data) + sim.just_solve() + + even_newer_voltages = sim.get_voltages_actual() + + assert ( + float( + np.sum(np.abs(new_voltages.loc["113.1"] - even_newer_voltages.loc["113.1"])) + ) + < 1 + ) + + test_inverter_data = InverterControl( + pvsystem_list=["PVSystem.113"], + vvcontrol=VVControl( + deltaq_factor=0.7, + voltage=[0.5, 0.945, 0.97, 0.99, 1.045, 1.50], + reactive_response=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ), + ) + + sim.apply_inverter_control(test_inverter_data) + sim.just_solve() + + nocontrol_voltages = sim.get_voltages_actual() + + assert ( + float( + np.sum(np.abs(old_voltages.loc["113.1"] - nocontrol_voltages.loc["113.1"])) + ) + < 1 + ) + + +def test_inv_control_full(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + test_inverter_data = InverterControl( + vvcontrol=VVControl( + voltage=[0.5, 0.945, 0.97, 0.99, 1.045, 1.50], + reactive_response=[1, 1, 1, 1, 1, 1], + ), + ) + + _ = sim.get_y_matrix() + sim.snapshot_run() # Needed to bring out of disabled state + sim.solve(8, 0) + + old_voltages = sim.get_voltages_actual() + + sim.apply_inverter_control(test_inverter_data) + sim.just_solve() + + new_voltages = sim.get_voltages_actual() + + assert float( + np.sum(np.abs(new_voltages.loc["113.1"] - old_voltages.loc["113.1"])) + ) > 0.01 * float(np.abs(old_voltages.loc["113.1"])) + + assert float( + np.sum(np.abs(new_voltages.loc["87.3"] - old_voltages.loc["87.3"])) + ) > 0.01 * float(np.abs(old_voltages.loc["87.3"])) + + +def test_inv_combined_control(federate_config): + logging.info("Loading sim") + sim = FeederSimulator.FeederSimulator(federate_config) + + _ = sim.get_y_matrix() + sim.snapshot_run() # Needed to bring out of disabled state + sim.solve(8, 0) + + old_voltages = sim.get_voltages_actual() + + test_inverter_data = InverterControl( + vvcontrol=VVControl( + deltaq_factor=0.7, + voltage=[0.5, 0.945, 0.97, 0.99, 1.045, 1.50], + reactive_response=[1, 1, 1, 1, 1, 1], + ), + vwcontrol=VWControl( + deltap_factor=1.0, + voltage=[0.5, 0.945, 0.97, 0.99, 1.045, 1.50], + power_response=[0.1, 0.1, 0.1, 0.1, 0.1, 0.1], + ), + mode=InverterControlMode.voltvar_voltwatt, + ) + + sim.apply_inverter_control(test_inverter_data) + sim.just_solve() + + new_voltages = sim.get_voltages_actual() + + assert float( + np.sum(np.abs(new_voltages.loc["113.1"] - old_voltages.loc["113.1"])) + ) > 0.01 * float(np.abs(old_voltages.loc["113.1"])) + + assert float( + np.sum(np.abs(new_voltages.loc["87.3"] - old_voltages.loc["87.3"])) + ) > 0.01 * float(np.abs(old_voltages.loc["87.3"])) diff --git a/README.md b/README.md index 143f779..3bcbc77 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# GADAL-example +# oedisi-example This example shows how to use the GADAL api to manage simulations. We also use it as a testing ground for the testing the combination of feeders, @@ -6,24 +6,17 @@ state estimation, and distributed OPF. # Install and Running Locally -1. Install the GADAL `componentframework` using `python setup.py develop` or using -`pip install git+ssh://git@github.com/openEDI/GADAL@v0.2.2`. - -To run the simulation, you'll need several libraries such as OpenDSSDirect.py and pyarrow. +1. To run the simulation, you'll need several libraries such as OpenDSSDirect.py and pyarrow. ``` pip install -r requirements.txt ``` -For analysis, you'll also need matplotlib. -``` -pip install matplotlib` -``` -2. Run `python test_full_systems.py` to initialize the system -defined in `test_system.json` in a `build` directory. +2. Run `oedisi build --system scenarios/docker_system.json` to initialize the system +defined in `scenarios/test_system.json` in a `build` directory. You can specify your own directory with `--system` and your own system json with `--system`. -3. Run `helics run --path=build/test_system_runner.json` +3. Run `oedisi run` 4. Analyze the results using `python post_analysis.py` This computes some percentage relative errors in magnitude (MAPE) and angle (MAE), @@ -38,10 +31,15 @@ If the simulation fails, you may **need** to kill the `helics_broker` manually b When debugging, you should check the `.log` files for errors. Error code `-9` usually occurs when it is killed by the broker as opposed to failing directly. -# Components +You can use the `oedisi` CLI tools to help debug specific components or timing. + +- `oedisi run-with-pause` +- `oedisi debug-component --foreground feeder` + +# Components All the required components are defined in folders within this repo. Each component -pulls types from `gadal.data_types`. +pulls types from `oedisi.types.data_types`. ![Block diagram of simulation](sgidal-example.png) @@ -87,7 +85,7 @@ information about the inputs and outputs of each component. We created component python scripts that matched these component descriptions and followed the GADAL API for configuration. -In order to use the data types from other federates, the `gadal.gadal_types` +In order to use the data types from other federates, the `oedisi.types` module is critical. If additional data is needed, then we recommend subclassing the pydantic models and adding the data in the required federates as needed. This ensures that others should still be able to parse your types if @@ -114,19 +112,14 @@ for each component with the right configuration. # Docker Container -One of the downsides of having `gadal` as a private library at present is that it complicates automated -installation. We could do this by copying over a tar file, or we can use it by pip installing -it with an SSH key. Currently, we use an SSH key. - -Assuming the github SSH key is in the current directory `gadal_docker_key`, we can build the docker image with ``` -docker build --secret id=gadal_github_key,src=gadal_docker_key -t gadal-example:0.0.0 . +docker build -t oedisi-example:0.0.0 . ``` To get a docker volume pointed at the right place locally, we have to run more commands ``` mkdir outputs_build -docker volume create --name gadal_output --opt type=none --opt device=$(PWD)/outputs_build --opt o=bind +docker volume create --name oedisi_output --opt type=none --opt device=$(PWD)/outputs_build --opt o=bind ``` If `pwd` is unavailable on your system, then you must specify the exact path. On windows, this will end up @@ -134,7 +127,7 @@ being `/c/Users/.../outputs_builds/`. You must use forward slashes. Then we can run the docker image: ``` -docker run --rm --mount source=gadal_output,target=/simulation/outputs gadal-example:0.0.0 +docker run --rm --mount source=oedisi_output,target=/simulation/outputs oedisi-example:0.0.0 ``` You can omit the docker volume parts as well as `--mount` if you do not care about the exact outputs. diff --git a/components.json b/components.json new file mode 100644 index 0000000..599dfdf --- /dev/null +++ b/components.json @@ -0,0 +1,6 @@ +{ + "Recorder": "recorder/component_definition.json", + "LocalFeeder": "LocalFeeder/component_definition.json", + "MeasurementComponent": "measuring_federate/component_definition.json", + "StateEstimatorComponent": "wls_federate/component_definition.json" +} diff --git a/errors.png b/errors.png index 39275ca..6623ab1 100644 Binary files a/errors.png and b/errors.png differ diff --git a/measuring_federate/generate_test_config.py b/measuring_federate/generate_test_config.py index 2d95e70..50eae92 100644 --- a/measuring_federate/generate_test_config.py +++ b/measuring_federate/generate_test_config.py @@ -3,9 +3,209 @@ from typing import List -BUSES = ['P1UDT882LV.1', 'P1UDT882LV.2', 'P1ULV4109.1', 'P1ULV4109.2', 'P1UDM3853.1', 'P1UDM3853.2', 'P1ULV4110.1', 'P1ULV4110.2', 'P1UDT883LV.1', 'P1UDT883LV.2', 'P1UDM3854.1', 'P1UDM3854.2', 'P1ULV4111.1', 'P1ULV4111.2', 'P1UDT884LV.1', 'P1UDT884LV.2', 'P1ULV4112.1', 'P1ULV4112.2', 'P1UDM3856.1', 'P1UDM3856.2', 'P1ULV4113.1', 'P1ULV4113.2', 'P1UDT893LV.1', 'P1UDT893LV.2', 'P1ULV4136.1', 'P1ULV4136.2', 'P1UDM3883.1', 'P1UDM3883.2', 'P1ULV4137.1', 'P1ULV4137.2', 'P1UDT894LV.1', 'P1UDT894LV.2', 'P1UDM3884.1', 'P1UDM3884.2', 'P1ULV4138.1', 'P1ULV4138.2', 'P1UDT895LV.1', 'P1UDT895LV.2', 'P1ULV7543.1', 'P1ULV7543.2', 'P1UDM3885.1', 'P1UDM3885.2', 'P1ULV4139.1', 'P1ULV4139.2', 'P1UDM3887.1', 'P1UDM3887.2', 'P1ULV4140.1', 'P1ULV4140.2', 'P1UDM3889.1', 'P1UDM3889.2', 'P1ULV7633.1', 'P1ULV7633.2', 'P1UDT896LV.1', 'P1UDT896LV.2', 'P1UDM3888.1', 'P1UDM3888.2', 'P1ULV4197.1', 'P1ULV4197.2', 'P1UDT909LV.1', 'P1UDT909LV.2', 'P1ULV4155.1', 'P1ULV4155.2', 'P1UDT942LV.1', 'P1UDT942LV.2', 'P1ULV4193.1', 'P1ULV4193.2', 'P1UDT943LV.1', 'P1UDT943LV.2', 'P1ULV4196.1', 'P1ULV4196.2', 'P1UDT1431LV.1', 'P1UDT1431LV.2', 'P1ULV7544.1', 'P1ULV7544.2', 'P1UDT1432LV.1', 'P1UDT1432LV.2', 'P1ULV7546.1', 'P1ULV7546.2', 'P1UDT1433LV.1', 'P1UDT1433LV.2', 'P1ULV7547.1', 'P1ULV7547.2', 'P1UDT908-P1UDT909X.2', 'P1UDT909.2', 'P1UDT1433-P1UDT882X.3', 'P1UDT882.3', 'P1UDT1432-P1UDT881X.2', 'P1UDT1432.2', 'P1UDT896-P1UDT938X.3', 'P1UDT896.3', 'P1UDT942-P1UHS0_1247X.1', 'P1UDT942-P1UHS0_1247X.2', 'P1UDT942-P1UHS0_1247X.3', 'P1UDT942.1', 'P1UDT942.2', 'P1UDT942.3', 'P1UDT907LV.1', 'P1UDT907LV.2', 'P1UDT907LV.3', 'P1ULV4153.1', 'P1ULV4153.2', 'P1ULV4153.3', 'P1UDT908LV.1', 'P1UDT908LV.2', 'P1UDT908LV.3', 'P1ULV4154.1', 'P1ULV4154.2', 'P1ULV4154.3', 'P1UDT907.1', 'P1UDT907.2', 'P1UDT907.3', 'P1UDT908.1', 'P1UDT908.2', 'P1UDT908.3', 'P1UDT893.1', 'P1UDT893.2', 'P1UDT893.3', 'P1UDT894.1', 'P1UDT894.2', 'P1UDT894.3', 'P1UDT1425LV.1', 'P1UDT1425LV.2', 'P1UDT1425LV.3', 'P1ULV7531.1', 'P1ULV7531.2', 'P1ULV7531.3', 'P1UDT881LV.1', 'P1UDT881LV.2', 'P1UDT881LV.3', 'P1ULV4108.1', 'P1ULV4108.2', 'P1ULV4108.3', 'P1UDT1425-P1UDT881X.1', 'P1UDT1425-P1UDT881X.2', 'P1UDT1425-P1UDT881X.3', 'P1UDT881.1', 'P1UDT881.2', 'P1UDT881.3', 'P1UDT1433.1', 'P1UDT1433.2', 'P1UDT1433.3', 'P1UDT1425.1', 'P1UDT1425.2', 'P1UDT1425.3', 'P1UDT883.1', 'P1UDT883.2', 'P1UDT883.3', 'P1UDT884.1', 'P1UDT884.2', 'P1UDT884.3', 'P1UDT1431.1', 'P1UDT1431.2', 'P1UDT1431.3', 'P1UDT1431-P1UDT895X.1', 'P1UDT1431-P1UDT895X.2', 'P1UDT1431-P1UDT895X.3', 'P1UDT894-P1UDT895X.1', 'P1UDT894-P1UDT895X.2', 'P1UDT894-P1UDT895X.3', 'P1UDT895.1', 'P1UDT895.2', 'P1UDT895.3', 'P1UDT894-P1UDT942X.1', 'P1UDT894-P1UDT942X.2', 'P1UDT894-P1UDT942X.3', 'P1UDT939LV.1', 'P1UDT939LV.2', 'P1UDT939LV.3', 'P1ULV4189.1', 'P1ULV4189.2', 'P1ULV4189.3', 'P1UDT938LV.1', 'P1UDT938LV.2', 'P1UDT938LV.3', 'P1ULV4188.1', 'P1ULV4188.2', 'P1ULV4188.3', 'P1UDT892LV.1', 'P1UDT892LV.2', 'P1UDT892LV.3', 'P1ULV4132.1', 'P1ULV4132.2', 'P1ULV4132.3', 'P1UDT938.1', 'P1UDT938.2', 'P1UDT938.3', 'P1UDT892.1', 'P1UDT892.2', 'P1UDT892.3', 'P1UDT943.1', 'P1UDT943.2', 'P1UDT943.3', 'P1UDT939.1', 'P1UDT939.2', 'P1UDT939.3'] +BUSES = [ + "P1UDT882LV.1", + "P1UDT882LV.2", + "P1ULV4109.1", + "P1ULV4109.2", + "P1UDM3853.1", + "P1UDM3853.2", + "P1ULV4110.1", + "P1ULV4110.2", + "P1UDT883LV.1", + "P1UDT883LV.2", + "P1UDM3854.1", + "P1UDM3854.2", + "P1ULV4111.1", + "P1ULV4111.2", + "P1UDT884LV.1", + "P1UDT884LV.2", + "P1ULV4112.1", + "P1ULV4112.2", + "P1UDM3856.1", + "P1UDM3856.2", + "P1ULV4113.1", + "P1ULV4113.2", + "P1UDT893LV.1", + "P1UDT893LV.2", + "P1ULV4136.1", + "P1ULV4136.2", + "P1UDM3883.1", + "P1UDM3883.2", + "P1ULV4137.1", + "P1ULV4137.2", + "P1UDT894LV.1", + "P1UDT894LV.2", + "P1UDM3884.1", + "P1UDM3884.2", + "P1ULV4138.1", + "P1ULV4138.2", + "P1UDT895LV.1", + "P1UDT895LV.2", + "P1ULV7543.1", + "P1ULV7543.2", + "P1UDM3885.1", + "P1UDM3885.2", + "P1ULV4139.1", + "P1ULV4139.2", + "P1UDM3887.1", + "P1UDM3887.2", + "P1ULV4140.1", + "P1ULV4140.2", + "P1UDM3889.1", + "P1UDM3889.2", + "P1ULV7633.1", + "P1ULV7633.2", + "P1UDT896LV.1", + "P1UDT896LV.2", + "P1UDM3888.1", + "P1UDM3888.2", + "P1ULV4197.1", + "P1ULV4197.2", + "P1UDT909LV.1", + "P1UDT909LV.2", + "P1ULV4155.1", + "P1ULV4155.2", + "P1UDT942LV.1", + "P1UDT942LV.2", + "P1ULV4193.1", + "P1ULV4193.2", + "P1UDT943LV.1", + "P1UDT943LV.2", + "P1ULV4196.1", + "P1ULV4196.2", + "P1UDT1431LV.1", + "P1UDT1431LV.2", + "P1ULV7544.1", + "P1ULV7544.2", + "P1UDT1432LV.1", + "P1UDT1432LV.2", + "P1ULV7546.1", + "P1ULV7546.2", + "P1UDT1433LV.1", + "P1UDT1433LV.2", + "P1ULV7547.1", + "P1ULV7547.2", + "P1UDT908-P1UDT909X.2", + "P1UDT909.2", + "P1UDT1433-P1UDT882X.3", + "P1UDT882.3", + "P1UDT1432-P1UDT881X.2", + "P1UDT1432.2", + "P1UDT896-P1UDT938X.3", + "P1UDT896.3", + "P1UDT942-P1UHS0_1247X.1", + "P1UDT942-P1UHS0_1247X.2", + "P1UDT942-P1UHS0_1247X.3", + "P1UDT942.1", + "P1UDT942.2", + "P1UDT942.3", + "P1UDT907LV.1", + "P1UDT907LV.2", + "P1UDT907LV.3", + "P1ULV4153.1", + "P1ULV4153.2", + "P1ULV4153.3", + "P1UDT908LV.1", + "P1UDT908LV.2", + "P1UDT908LV.3", + "P1ULV4154.1", + "P1ULV4154.2", + "P1ULV4154.3", + "P1UDT907.1", + "P1UDT907.2", + "P1UDT907.3", + "P1UDT908.1", + "P1UDT908.2", + "P1UDT908.3", + "P1UDT893.1", + "P1UDT893.2", + "P1UDT893.3", + "P1UDT894.1", + "P1UDT894.2", + "P1UDT894.3", + "P1UDT1425LV.1", + "P1UDT1425LV.2", + "P1UDT1425LV.3", + "P1ULV7531.1", + "P1ULV7531.2", + "P1ULV7531.3", + "P1UDT881LV.1", + "P1UDT881LV.2", + "P1UDT881LV.3", + "P1ULV4108.1", + "P1ULV4108.2", + "P1ULV4108.3", + "P1UDT1425-P1UDT881X.1", + "P1UDT1425-P1UDT881X.2", + "P1UDT1425-P1UDT881X.3", + "P1UDT881.1", + "P1UDT881.2", + "P1UDT881.3", + "P1UDT1433.1", + "P1UDT1433.2", + "P1UDT1433.3", + "P1UDT1425.1", + "P1UDT1425.2", + "P1UDT1425.3", + "P1UDT883.1", + "P1UDT883.2", + "P1UDT883.3", + "P1UDT884.1", + "P1UDT884.2", + "P1UDT884.3", + "P1UDT1431.1", + "P1UDT1431.2", + "P1UDT1431.3", + "P1UDT1431-P1UDT895X.1", + "P1UDT1431-P1UDT895X.2", + "P1UDT1431-P1UDT895X.3", + "P1UDT894-P1UDT895X.1", + "P1UDT894-P1UDT895X.2", + "P1UDT894-P1UDT895X.3", + "P1UDT895.1", + "P1UDT895.2", + "P1UDT895.3", + "P1UDT894-P1UDT942X.1", + "P1UDT894-P1UDT942X.2", + "P1UDT894-P1UDT942X.3", + "P1UDT939LV.1", + "P1UDT939LV.2", + "P1UDT939LV.3", + "P1ULV4189.1", + "P1ULV4189.2", + "P1ULV4189.3", + "P1UDT938LV.1", + "P1UDT938LV.2", + "P1UDT938LV.3", + "P1ULV4188.1", + "P1ULV4188.2", + "P1ULV4188.3", + "P1UDT892LV.1", + "P1UDT892LV.2", + "P1UDT892LV.3", + "P1ULV4132.1", + "P1ULV4132.2", + "P1ULV4132.3", + "P1UDT938.1", + "P1UDT938.2", + "P1UDT938.3", + "P1UDT892.1", + "P1UDT892.2", + "P1UDT892.3", + "P1UDT943.1", + "P1UDT943.2", + "P1UDT943.3", + "P1UDT939.1", + "P1UDT939.2", + "P1UDT939.3", +] -SLACK_BUS = ['P1UDT942-P1UHS0_1247X.1', 'P1UDT942-P1UHS0_1247X.2', 'P1UDT942-P1UHS0_1247X.3'] +SLACK_BUS = [ + "P1UDT942-P1UHS0_1247X.1", + "P1UDT942-P1UHS0_1247X.2", + "P1UDT942-P1UHS0_1247X.3", +] class MeasurementConfig(BaseModel): @@ -20,10 +220,12 @@ class MeasurementConfig(BaseModel): seed(12345) for i in SLACK_BUS: BUSES.remove(i) - print(MeasurementConfig( - name = "MeasurementFederate", - gaussian_variance = 0.1, - voltage_ids = sample(BUSES, 3*len(BUSES) // 4), - real_power_ids = sample(BUSES, 3*len(BUSES) // 4), - reactive_power_ids = sample(BUSES, 3*len(BUSES) // 4) - ).json()) + print( + MeasurementConfig( + name="MeasurementFederate", + gaussian_variance=0.1, + voltage_ids=sample(BUSES, 3 * len(BUSES) // 4), + real_power_ids=sample(BUSES, 3 * len(BUSES) // 4), + reactive_power_ids=sample(BUSES, 3 * len(BUSES) // 4), + ).json() + ) diff --git a/measuring_federate/measuring_federate.py b/measuring_federate/measuring_federate.py index 20b2324..5d5acae 100644 --- a/measuring_federate/measuring_federate.py +++ b/measuring_federate/measuring_federate.py @@ -6,17 +6,19 @@ import scipy.io import json from datetime import datetime -from gadal.gadal_types.data_types import MeasurementArray +from oedisi.types.data_types import MeasurementArray, EquipmentNodeArray logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) logger.setLevel(logging.INFO) + class MeasurementConfig(BaseModel): name: str gaussian_variance: float measurement_file: str random_percent: float + run_freq_time_step: float = 1.0 def get_indices(labelled_array, indices): @@ -25,29 +27,46 @@ def get_indices(labelled_array, indices): return [inv_map[v] for v in labelled_array.ids] -def reindex(measurement_array, indices): +def reindex(measurement_array: MeasurementArray, indices): inv_map = {v: i for i, v in enumerate(measurement_array.ids)} - for i in inv_map: - logger.debug(i) - return MeasurementArray(values=[ - measurement_array.values[inv_map[i]] for i in indices - ], ids=indices, units = measurement_array.units, equipment_type = measurement_array.equipment_type, time = measurement_array.time) + if isinstance(measurement_array, EquipmentNodeArray): + return measurement_array.__class__( + values=[measurement_array.values[inv_map[i]] for i in indices], + ids=indices, + units=measurement_array.units, + equipment_ids=[measurement_array.equipment_ids[inv_map[i]] for i in indices], + time=measurement_array.time, + ) + else: + return measurement_array.__class__( + values=[measurement_array.values[inv_map[i]] for i in indices], + ids=indices, + units=measurement_array.units, + time=measurement_array.time, + ) -def apply(f, measurement_array): - return MeasurementArray( - values=list(map(f, measurement_array.values)), - ids=measurement_array.ids, - units = measurement_array.units, - equipment_type = measurement_array.equipment_type, - time = measurement_array.time - ) +def apply(f, measurement_array: MeasurementArray): + if isinstance(measurement_array, EquipmentNodeArray): + return measurement_array.__class__( + values=list(map(f, measurement_array.values)), + ids=measurement_array.ids, + units=measurement_array.units, + equipment_ids=measurement_array.equipment_ids, + time=measurement_array.time, + ) + else: + return measurement_array.__class__( + values=list(map(f, measurement_array.values)), + ids=measurement_array.ids, + units=measurement_array.units, + time=measurement_array.time, + ) class MeasurementRelay: def __init__(self, config: MeasurementConfig, input_mapping): self.rng = np.random.default_rng(12345) - deltat = 0.01 # deltat = 60. # Create Federate Info object that describes the federate properties # @@ -58,7 +77,7 @@ def __init__(self, config: MeasurementConfig, input_mapping): logger.debug(config.name) h.helicsFederateInfoSetTimeProperty( - fedinfo, h.helics_property_time_delta, deltat + fedinfo, h.helics_property_time_delta, config.run_freq_time_step ) self.vfed = h.helicsCreateValueFederate(config.name, fedinfo) @@ -69,7 +88,7 @@ def __init__(self, config: MeasurementConfig, input_mapping): input_mapping["subscription"], "" ) - #TODO: find better way to determine what the name of this federate instance is than looking at the subscription + # TODO: find better way to determine what the name of this federate instance is than looking at the subscription self.pub_measurement = self.vfed.register_publication( "publication", h.HELICS_DATA_TYPE_STRING, "" ) @@ -82,7 +101,7 @@ def transform(self, measurement_array: MeasurementArray, unique_ids): new_array = reindex(measurement_array, unique_ids) return apply( lambda x: x + self.rng.normal(scale=np.sqrt(self.gaussian_variance)), - new_array + new_array, ) def run(self): @@ -92,11 +111,14 @@ def run(self): granted_time = h.helicsFederateRequestTime(self.vfed, h.HELICS_TIME_MAXTIME) while granted_time < h.HELICS_TIME_MAXTIME: - logger.info('start time: '+str(datetime.now())) + logger.info("start time: " + str(datetime.now())) json_data = self.sub_measurement.json - measurement = MeasurementArray(**json_data) + if "equipment_ids" in json_data: + measurement = EquipmentNodeArray.parse_obj(json_data) + else: + measurement = MeasurementArray.parse_obj(json_data) - with open(self.measurement_file,'r') as fp: + with open(self.measurement_file, "r") as fp: self.measurement = json.load(fp) measurement_transformed = self.transform(measurement, self.measurement) logger.debug("measured transformed") @@ -105,7 +127,7 @@ def run(self): self.pub_measurement.publish(measurement_transformed.json()) granted_time = h.helicsFederateRequestTime(self.vfed, h.HELICS_TIME_MAXTIME) - logger.info('end time: '+str(datetime.now())) + logger.info("end time: " + str(datetime.now())) self.destroy() diff --git a/post_analysis.py b/post_analysis.py index c81bc15..b8b6837 100644 --- a/post_analysis.py +++ b/post_analysis.py @@ -5,52 +5,86 @@ import pandas as pd import json import os -from gadal.gadal_types.data_types import MeasurementArray, AdmittanceMatrix, Topology +from oedisi.types.data_types import MeasurementArray, AdmittanceMatrix, Topology import argparse -parser = argparse.ArgumentParser(description='Create plots') -parser.add_argument("directory", nargs="?", default="outputs", - help="directory to look for voltage_measurements") +parser = argparse.ArgumentParser(description="Create plots") +parser.add_argument( + "directory", + nargs="?", + default="outputs", + help="directory to look for voltage_measurements", +) args = parser.parse_args() -voltage_real = feather.read_feather(os.path.join(args.directory, "voltage_real.feather")) -voltage_imag = feather.read_feather(os.path.join(args.directory, "voltage_imag.feather")) +voltage_real = feather.read_feather( + os.path.join(args.directory, "voltage_real.feather") +) +voltage_imag = feather.read_feather( + os.path.join(args.directory, "voltage_imag.feather") +) with open(os.path.join(args.directory, "topology.json")) as f: topology = Topology.parse_obj(json.load(f)) - base_voltages = np.array(topology.base_voltage_magnitudes.values) + base_voltage_df = pd.DataFrame( + { + "id": topology.base_voltage_magnitudes.ids, + "value": topology.base_voltage_magnitudes.values, + } + ) + base_voltage_df.set_index("id", inplace=True) + base_voltages = base_voltage_df["value"] + +true_voltages = voltage_real.drop("time", axis=1) + 1j * voltage_imag.drop( + "time", axis=1 +) +true_voltages["time"] = voltage_real["time"] +true_voltages.set_index("time", inplace=True) + +voltage_mag = feather.read_feather(os.path.join(args.directory, "voltage_mag.feather")) +estimated_time = voltage_mag["time"] +voltage_mag.drop("time", axis=1) +voltage_angle = feather.read_feather( + os.path.join(args.directory, "voltage_angle.feather") +).drop("time", axis=1) -true_voltages = voltage_real.drop('time', axis=1) + 1j * voltage_imag.drop('time', axis=1) -time = voltage_real["time"] +estimated_voltages = voltage_mag * np.exp(1j * voltage_angle) +estimated_voltages["time"] = estimated_time +estimated_voltages.set_index("time", inplace=True) -voltage_mag = feather.read_feather(os.path.join(args.directory, "voltage_mag.feather")).drop('time', axis=1) -voltage_angle = feather.read_feather(os.path.join(args.directory, "voltage_angle.feather")).drop('time', axis=1) +time_intersection = pd.merge( + true_voltages, estimated_voltages, left_index=True, right_index=True +).index.to_numpy() + +estimated_voltages = estimated_voltages.loc[time_intersection, :] +true_voltages = true_voltages.loc[time_intersection, :] + +estimated_voltages = estimated_voltages.reindex(true_voltages.columns, axis=1) -estimated_voltages = voltage_mag * np.exp(1j * voltage_angle) def plots(true_voltages, estimated_voltages, time=0, unit="kV"): n_nodes = true_voltages.shape[0] x_axis = np.arange(n_nodes) - fig1, ax = plt.subplots(figsize=(10,10)) + fig1, ax = plt.subplots(figsize=(10, 10)) ax.bar(x_axis, np.angle(estimated_voltages)) ax.bar(x_axis, np.angle(true_voltages), width=0.5) - #ax.set_xticks(x_axis, true_voltages.index, rotation=-90, fontsize=5) - #ax.set_tick_params(axis='x', labelsize=5, rotation=-90) - ax.set_xlabel('Node number') - ax.set_ylabel('Voltage Angles') - ax.legend(['Estimated', 'True']) + # ax.set_xticks(x_axis, true_voltages.index, rotation=-90, fontsize=5) + # ax.set_tick_params(axis='x', labelsize=5, rotation=-90) + ax.set_xlabel("Node number") + ax.set_ylabel("Voltage Angles") + ax.legend(["Estimated", "True"]) ax.set_title(f"Voltage Angles at t={time}") - fig2, ax = plt.subplots(figsize=(10,10)) - ax.plot(x_axis, np.abs(estimated_voltages), '-o') - ax.plot(x_axis, np.abs(true_voltages), '-o') - #ax.set_xticks(x_axis, true_voltages.index, rotation=-90, fontsize=5) - ax.set_xlabel('Node number') - ax.set_ylabel(f'Voltage Magnitudes ({unit})') - ax.legend(['Estimated', 'True']) + fig2, ax = plt.subplots(figsize=(10, 10)) + ax.plot(x_axis, np.abs(estimated_voltages), "-o") + ax.plot(x_axis, np.abs(true_voltages), "-o") + # ax.set_xticks(x_axis, true_voltages.index, rotation=-90, fontsize=5) + ax.set_xlabel("Node number") + ax.set_ylabel(f"Voltage Magnitudes ({unit})") + ax.legend(["Estimated", "True"]) ax.set_title(f"Voltage Magnitudes at t={time}") return fig1, fig2 @@ -59,20 +93,23 @@ def errors(true_voltages, estimated_voltages): true_mag = np.abs(true_voltages) nonzero_parts = true_mag != 0.0 MAPE = np.mean( - np.array(np.abs(true_mag - np.abs(estimated_voltages)) - / true_mag)[nonzero_parts] + np.array(np.abs(true_mag - np.abs(estimated_voltages)) / true_mag)[ + nonzero_parts + ] * 100 ) angle_difference = np.abs(np.angle(true_voltages) - np.angle(estimated_voltages)) - angle_difference[angle_difference >= np.pi] = 2*np.pi - angle_difference[angle_difference >= np.pi] + angle_difference[angle_difference >= np.pi] = ( + 2 * np.pi - angle_difference[angle_difference >= np.pi] + ) MAE = np.mean(np.array(angle_difference)[nonzero_parts] * 180 / np.pi) return MAPE, MAE def error_table(true_voltages, estimated_voltages): error_table = [] - for i, t in enumerate(time): - MAPE, MAE = errors(true_voltages.iloc[i,:], estimated_voltages.iloc[i,:]) + for i, t in enumerate(true_voltages.index): + MAPE, MAE = errors(true_voltages.iloc[i, :], estimated_voltages.iloc[i, :]) error_table.append({"t": t, "MAPE": MAPE, "MAE": MAE}) return pd.DataFrame(error_table) @@ -88,15 +125,24 @@ def plot_errors(err_table): ax.set_xticks(err_table["t"][::5], err_table["t"][::5], rotation=-25, fontsize=5) return fig + err_table = error_table(true_voltages, estimated_voltages) plot_errors(err_table).savefig("errors.png") MAPE, MAE = errors(true_voltages, estimated_voltages) print(f"MAPE = {MAPE}, MAE={MAE}") -fig1, fig2 = plots(true_voltages.iloc[0,:] / base_voltages, estimated_voltages.iloc[0,:] / base_voltages, unit="p.u.") +fig1, fig2 = plots( + true_voltages.iloc[0, :] / base_voltages, + estimated_voltages.iloc[0, :] / base_voltages, + unit="p.u.", +) fig1.savefig("voltage_angles_0.png") fig2.savefig("voltage_magnitudes_0.png") -if len(true_voltages) >=96: - fig1, fig2 = plots(true_voltages.iloc[95,:] / base_voltages, estimated_voltages.iloc[95,:] / base_voltages, time=95, unit="p.u.") +if len(true_voltages) >= 94: + fig1, fig2 = plots( + true_voltages.iloc[93, :] / base_voltages, + estimated_voltages.iloc[93, :] / base_voltages, + time=93, + unit="p.u.", + ) fig1.savefig("voltage_angles_95.png") fig2.savefig("voltage_magnitudes_95.png") - diff --git a/recorder/record_subscription.py b/recorder/record_subscription.py index f4cbfed..7038491 100644 --- a/recorder/record_subscription.py +++ b/recorder/record_subscription.py @@ -8,12 +8,13 @@ import csv import pyarrow as pa from datetime import datetime -from gadal.gadal_types.data_types import MeasurementArray +from oedisi.types.data_types import MeasurementArray logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) logger.setLevel(logging.INFO) + class Recorder: def __init__(self, name, feather_filename, csv_filename, input_mapping): self.rng = np.random.default_rng(12345) @@ -35,9 +36,7 @@ def __init__(self, name, feather_filename, csv_filename, input_mapping): logger.info("Value federate created") # Register the publication # - self.sub = self.vfed.register_subscription( - input_mapping["subscription"], "" - ) + self.sub = self.vfed.register_subscription(input_mapping["subscription"], "") self.feather_filename = feather_filename self.csv_filename = csv_filename @@ -50,34 +49,39 @@ def run(self): start = True granted_time = h.helicsFederateRequestTime(self.vfed, h.HELICS_TIME_MAXTIME) - with pa.OSFile(self.feather_filename, 'wb') as sink: + with pa.OSFile(self.feather_filename, "wb") as sink: writer = None while granted_time < h.HELICS_TIME_MAXTIME: - logger.info('start time: '+str(datetime.now())) + logger.info("start time: " + str(datetime.now())) logger.debug(granted_time) # Check that the data is a MeasurementArray type json_data = self.sub.json - json_data['time'] = granted_time + json_data["time"] = granted_time measurement = MeasurementArray(**self.sub.json) - measurement_dict = {key: value for key, value in zip(measurement.ids,measurement.values)} - measurement_dict['time'] = measurement.time.strftime("%Y-%m-%d %H:%M:%S") + measurement_dict = { + key: value + for key, value in zip(measurement.ids, measurement.values) + } + measurement_dict["time"] = measurement.time.strftime( + "%Y-%m-%d %H:%M:%S" + ) logger.debug(measurement.time) if start: schema_elements = [(key, pa.float64()) for key in measurement.ids] - schema_elements.append(('time',pa.string())) + schema_elements.append(("time", pa.string())) schema = pa.schema(schema_elements) writer = pa.ipc.new_file(sink, schema) start = False cnt = 0 - writer.write_batch(pa.RecordBatch.from_pylist([ - measurement_dict - ])) + writer.write_batch(pa.RecordBatch.from_pylist([measurement_dict])) - granted_time = h.helicsFederateRequestTime(self.vfed, h.HELICS_TIME_MAXTIME) - logger.info('end time: '+str(datetime.now())) + granted_time = h.helicsFederateRequestTime( + self.vfed, h.HELICS_TIME_MAXTIME + ) + logger.info("end time: " + str(datetime.now())) if writer is not None: writer.close() diff --git a/requirements.txt b/requirements.txt index 53ef9ca..41cfdcd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,12 @@ -helics==3.2.1.post15 -helics-apps==3.2.1 +helics==3.4.0 +helics-apps==3.4.0 pydantic pyarrow scipy matplotlib numpy pandas -OpenDSSDirect.py[extras] +OpenDSSDirect.py[extras]==0.7.0 boto3 +xarray +oedisi==1.0.0 diff --git a/test_full_systems.py b/test_full_systems.py deleted file mode 100644 index 7f169ad..0000000 --- a/test_full_systems.py +++ /dev/null @@ -1,61 +0,0 @@ -#!/usr/bin/env python3 -# Load component framework -from gadal.componentframework.basic_component import component_from_json -from gadal.componentframework.system_configuration import ( - generate_runner_config, - WiringDiagram, -) -import argparse - -parser = argparse.ArgumentParser(description="Build system") -parser.add_argument( - "--target-directory", - type=str, - default="build", - help="Target directory to put the system in", - metavar="PARAM", -) - -parser.add_argument( - "--system", - type=str, - default="scenarios/docker_system.json", - help="Wiring diagram json to build", - metavar="PARAM", -) - -args = parser.parse_args() - -def bad_type_checker(type, x): - "Doesn't do any type checking on the exchange types" - return True - -# We make classes for each component using a type checker -LocalFeeder = component_from_json( - "LocalFeeder/component_definition.json", bad_type_checker -) -MeasurementComponent = component_from_json( - "measuring_federate/component_definition.json", bad_type_checker -) -StateEstimatorComponent = component_from_json( - "wls_federate/component_definition.json", bad_type_checker -) -Recorder = component_from_json( - "recorder/component_definition.json", bad_type_checker -) - -# Dictionary used to interpret test_system.json -component_types = { - "LocalFeeder": LocalFeeder, - "MeasurementComponent": MeasurementComponent, - "StateEstimatorComponent": StateEstimatorComponent, - "Recorder": Recorder -} - -# Read wiring diagram (lists components, links, and parameters) -wiring_diagram = WiringDiagram.parse_file(args.system) -#wiring_diagram.clean_model() -# Generate runner config using wiring diagram and component types -runner_config = generate_runner_config(wiring_diagram, component_types, target_directory=args.target_directory) -with open(f"{args.target_directory}/test_system_runner.json", "w") as f: - f.write(runner_config.json(indent=2)) diff --git a/voltage_angles_0.png b/voltage_angles_0.png index 455e6ab..ec48f80 100644 Binary files a/voltage_angles_0.png and b/voltage_angles_0.png differ diff --git a/voltage_angles_95.png b/voltage_angles_95.png index 4869c31..7334b6c 100644 Binary files a/voltage_angles_95.png and b/voltage_angles_95.png differ diff --git a/voltage_magnitudes_0.png b/voltage_magnitudes_0.png index 2fb9d3a..9d1cc91 100644 Binary files a/voltage_magnitudes_0.png and b/voltage_magnitudes_0.png differ diff --git a/voltage_magnitudes_95.png b/voltage_magnitudes_95.png index 24ecd14..cc13493 100644 Binary files a/voltage_magnitudes_95.png and b/voltage_magnitudes_95.png differ diff --git a/wls_federate/state_estimator_federate.py b/wls_federate/state_estimator_federate.py index 4ea0cbe..77b6e3f 100644 --- a/wls_federate/state_estimator_federate.py +++ b/wls_federate/state_estimator_federate.py @@ -16,47 +16,65 @@ from typing import List, Optional, Union from scipy.optimize import least_squares from datetime import datetime -from gadal.gadal_types.data_types import AdmittanceSparse, MeasurementArray, AdmittanceMatrix, Topology, Complex, VoltagesMagnitude, VoltagesAngle, VoltagesReal, VoltagesImaginary, PowersReal, PowersImaginary +from oedisi.types.data_types import ( + AdmittanceSparse, + MeasurementArray, + AdmittanceMatrix, + Topology, + Complex, + VoltagesMagnitude, + VoltagesAngle, + VoltagesReal, + VoltagesImaginary, + PowersReal, + PowersImaginary, +) import scipy.sparse logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) logger.setLevel(logging.INFO) + def cal_h(knownP, knownQ, knownV, Y, deltaK, VabsK, num_node): - h1 = (VabsK[knownV]).reshape(-1,1) + h1 = (VabsK[knownV]).reshape(-1, 1) Vp = VabsK * np.exp(1j * deltaK) S = Vp * (Y.conjugate() @ Vp.conjugate()) P, Q = S.real, S.imag - h2, h3 = P[knownP].reshape(-1,1), Q[knownQ].reshape(-1,1) + h2, h3 = P[knownP].reshape(-1, 1), Q[knownQ].reshape(-1, 1) h = np.concatenate((h1, h2, h3), axis=0) return h.reshape(-1) + def cal_H(X0, z, num_node, knownP, knownQ, knownV, Y): deltaK, VabsK = X0[:num_node], X0[num_node:] num_knownV = len(knownV) - #Calculate original H1 + # Calculate original H1 H11, H12 = np.zeros((num_knownV, num_node)), np.zeros(num_knownV * num_node) - H12[np.arange(num_knownV)*num_node + knownV] = 1 + H12[np.arange(num_knownV) * num_node + knownV] = 1 H1 = np.concatenate((H11, H12.reshape(num_knownV, num_node)), axis=1) Vp = VabsK * np.exp(1j * deltaK) -##### S = np.diag(Vp) @ Y.conjugate() @ Vp.conjugate() -###### Take gradient with respect to V - H_pow2 = ( - Vp.reshape(-1, 1) * Y.conjugate() * np.exp(-1j * deltaK).reshape(1, -1) + - np.exp(1j * deltaK) * np.diag(Y.conjugate() @ Vp.conjugate()) - ) + ##### S = np.diag(Vp) @ Y.conjugate() @ Vp.conjugate() + ###### Take gradient with respect to V + H_pow2 = Vp.reshape(-1, 1) * Y.conjugate() * np.exp(-1j * deltaK).reshape( + 1, -1 + ) + np.exp(1j * deltaK) * np.diag(Y.conjugate() @ Vp.conjugate()) # Take gradient with respect to delta H_pow1 = ( - 1j * Vp.reshape(-1, 1) * (np.diag(Y.conjugate() @ Vp.conjugate()) - - Y.conjugate() * Vp.conjugate().reshape(1, -1)) - ) - + 1j + * Vp.reshape(-1, 1) + * ( + np.diag(Y.conjugate() @ Vp.conjugate()) + - Y.conjugate() * Vp.conjugate().reshape(1, -1) + ) + ) + H2 = np.concatenate((H_pow1.real, H_pow2.real), axis=1)[knownP, :] H3 = np.concatenate((H_pow1.imag, H_pow2.imag), axis=1)[knownQ, :] - H = np.concatenate((H1, H2, H3), axis=0) + H = np.concatenate((H1, H2, H3), axis=0) return -H + def residual(X0, z, num_node, knownP, knownQ, knownV, Y): delta, Vabs = X0[:num_node], X0[num_node:] h = cal_h(knownP, knownQ, knownV, Y, delta, Vabs, num_node) @@ -66,30 +84,26 @@ def residual(X0, z, num_node, knownP, knownQ, knownV, Y): logger.debug(z) logger.debug("h") logger.debug(h) - return z-h + return z - h def get_y(admittance: Union[AdmittanceMatrix, AdmittanceSparse], ids: List[str]): if type(admittance) == AdmittanceMatrix: assert ids == admittance.ids - return matrix_to_numpy( - admittance.admittance_matrix - ) + return matrix_to_numpy(admittance.admittance_matrix) elif type(admittance) == AdmittanceSparse: node_map = {name: i for (i, name) in enumerate(ids)} return scipy.sparse.coo_matrix( ( - [ - v[0] + 1j*v[1] - for v in admittance.admittance_list - ], + [v[0] + 1j * v[1] for v in admittance.admittance_list], ( [node_map[r] for r in admittance.from_equipment], - [node_map[c] for c in admittance.to_equipment] - ) + [node_map[c] for c in admittance.to_equipment], + ), ) ).toarray() + def matrix_to_numpy(admittance: List[List[Complex]]): "Convert list of list of our Complex type into a numpy matrix" return np.array([[x[0] + 1j * x[1] for x in row] for row in admittance]) @@ -100,9 +114,11 @@ def get_indices(topology, measurement): inv_map = {v: i for i, v in enumerate(topology.base_voltage_magnitudes.ids)} return [inv_map[v] for v in measurement.ids] + class UnitSystem(str, Enum): - SI = 'SI' - PER_UNIT = 'PER_UNIT' + SI = "SI" + PER_UNIT = "PER_UNIT" + class AlgorithmParameters(BaseModel): tol: float = 5e-7 @@ -112,7 +128,17 @@ class AlgorithmParameters(BaseModel): class Config: use_enum_values = True -def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_ang=0, initial_V=1, slack_index=0): + +def state_estimator( + parameters: AlgorithmParameters, + topology, + P, + Q, + V, + initial_ang=0, + initial_V=1, + slack_index=0, +): """Estimates voltage magnitude and angle from topology, partial power injections P + Q i, and lossy partial voltage magnitude. @@ -139,27 +165,34 @@ def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_ knownV = get_indices(topology, V) if parameters.units == UnitSystem.SI: - z = np.concatenate(( - V.array, -1000*np.array(P.array), -1000*np.array(Q.array) - ), axis=0) + z = np.concatenate( + (V.array, -1000 * np.array(P.array), -1000 * np.array(Q.array)), axis=0 + ) Y = get_y(topology.admittance, ids) elif parameters.units == UnitSystem.PER_UNIT: base_power = 100 if parameters.base_power != None: base_power = parameters.base_power - z = np.concatenate(( - V.values / base_voltages[knownV], - -np.array(P.values) / base_power, - -np.array(Q.values) / base_power - ), axis=0) + z = np.concatenate( + ( + V.values / base_voltages[knownV], + -np.array(P.values) / base_power, + -np.array(Q.values) / base_power, + ), + axis=0, + ) Y = get_y(topology.admittance, ids) # Hand-crafted unit conversion (check it, it works) - Y = base_voltages.reshape(1, -1) * Y * \ - base_voltages.reshape(-1, 1) / (base_power * 1000) + Y = ( + base_voltages.reshape(1, -1) + * Y + * base_voltages.reshape(-1, 1) + / (base_power * 1000) + ) else: raise Exception(f"Unit system {parameters.units} not supported") tol = parameters.tol - + if type(initial_ang) != np.ndarray: delta = np.full(num_node, initial_ang) else: @@ -177,22 +210,22 @@ def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_ logging.debug(delta) X0 = np.concatenate((delta, Vabs)) logging.debug(X0) - ang_low = np.concatenate(([-1e-5], np.ones(num_node -1)* (- np.inf))) - ang_up = np.concatenate(([1e-5], np.ones(num_node -1)* ( np.inf))) - mag_low = np.ones(num_node)* (- np.inf) - mag_up = np.ones(num_node)* (np.inf) + ang_low = np.concatenate(([-1e-5], np.ones(num_node - 1) * (-np.inf))) + ang_up = np.concatenate(([1e-5], np.ones(num_node - 1) * (np.inf))) + mag_low = np.ones(num_node) * (-np.inf) + mag_up = np.ones(num_node) * (np.inf) low_limit = np.concatenate((ang_low, mag_low)) up_limit = np.concatenate((ang_up, mag_up)) # Weights are ignored since errors are sampled from Gaussian # Real dimension of solutions is # 2 * num_node - len(knownP) - len(knownV) - len(knownQ) if len(knownP) + len(knownV) + len(knownQ) < num_node * 2: - #If not observable + # If not observable res_1 = least_squares( residual, X0, jac=cal_H, - #bounds = (low_limit, up_limit), + # bounds = (low_limit, up_limit), # method = 'lm', verbose=2, ftol=tol, @@ -205,8 +238,8 @@ def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_ residual, X0, jac=cal_H, - bounds = (low_limit, up_limit), - #method = 'lm', + bounds=(low_limit, up_limit), + # method = 'lm', verbose=2, ftol=tol, xtol=tol, @@ -223,13 +256,15 @@ def state_estimator(parameters: AlgorithmParameters, topology, P, Q, V, initial_ if parameters.units == UnitSystem.SI: return vmagestDecen, vangestDecen elif parameters.units == UnitSystem.PER_UNIT: - return vmagestDecen*(base_voltages), vangestDecen - + return vmagestDecen * (base_voltages), vangestDecen class StateEstimatorFederate: "State estimator federate. Wraps state_estimation with pubs and subs" - def __init__(self, federate_name, algorithm_parameters: AlgorithmParameters, input_mapping): + + def __init__( + self, federate_name, algorithm_parameters: AlgorithmParameters, input_mapping + ): "Initializes federate with name and remaps input into subscriptions" deltat = 0.1 @@ -281,9 +316,10 @@ def run(self): topology = Topology.parse_obj(self.sub_topology.json) ids = topology.base_voltage_magnitudes.ids logger.info("Topology has been read") - slack_index = None - if (not isinstance(topology.admittance, AdmittanceMatrix) and - not isinstance(topology.admittance, AdmittanceSparse)): + slack_index = None + if not isinstance(topology.admittance, AdmittanceMatrix) and not isinstance( + topology.admittance, AdmittanceSparse + ): raise "Weighted Least Squares algorithm expects AdmittanceMatrix/Sparse as input" for i in range(len(ids)): @@ -293,46 +329,55 @@ def run(self): while granted_time < h.HELICS_TIME_MAXTIME: if not self.sub_voltages_magnitude.is_updated(): - granted_time = h.helicsFederateRequestTime(self.vfed, h.HELICS_TIME_MAXTIME) + granted_time = h.helicsFederateRequestTime( + self.vfed, h.HELICS_TIME_MAXTIME + ) continue - logger.info('start time: '+str(datetime.now())) + logger.info("start time: " + str(datetime.now())) - voltages = VoltagesMagnitude.parse_obj(self.sub_voltages_magnitude.json) - power_Q = PowersImaginary.parse_obj(self.sub_power_Q.json) power_P = PowersReal.parse_obj(self.sub_power_P.json) + power_Q = PowersImaginary.parse_obj(self.sub_power_Q.json) knownP = get_indices(topology, power_P) knownQ = get_indices(topology, power_Q) knownV = get_indices(topology, voltages) - + if self.initial_V is None: - #Flat start or using average measurements + # Flat start or using average measurements if len(knownP) + len(knownV) + len(knownQ) > len(ids) * 2: self.initial_V = 1.0 else: - self.initial_V = np.mean(np.array(voltages.values) / np.array(topology.base_voltage_magnitudes.values)[knownV]) + self.initial_V = np.mean( + np.array(voltages.values) + / np.array(topology.base_voltage_magnitudes.values)[knownV] + ) if self.initial_ang is None: self.initial_ang = np.array(topology.base_voltage_angles.values) voltage_magnitudes, voltage_angles = state_estimator( self.algorithm_parameters, - topology, power_P, power_Q, voltages, initial_V=self.initial_V, - initial_ang=self.initial_ang, slack_index=slack_index + topology, + power_P, + power_Q, + voltages, + initial_V=self.initial_V, + initial_ang=self.initial_ang, + slack_index=slack_index, + ) + # self.initial_V = voltage_magnitudes + # self.initial_ang = voltage_angles + self.pub_voltage_mag.publish( + VoltagesMagnitude( + values=list(voltage_magnitudes), ids=ids, time=voltages.time + ).json() ) - #self.initial_V = voltage_magnitudes - #self.initial_ang = voltage_angles - self.pub_voltage_mag.publish(VoltagesMagnitude( - values=list(voltage_magnitudes), - ids=ids, - time = voltages.time - ).json()) - self.pub_voltage_angle.publish(VoltagesAngle( - values=list(voltage_angles), - ids=ids, - time = voltages.time - ).json()) - logger.info('end time: '+str(datetime.now())) + self.pub_voltage_angle.publish( + VoltagesAngle( + values=list(voltage_angles), ids=ids, time=voltages.time + ).json() + ) + logger.info("end time: " + str(datetime.now())) self.destroy() @@ -354,13 +399,8 @@ def destroy(self): else: parameters = AlgorithmParameters.parse_obj({}) - with open("input_mapping.json") as f: input_mapping = json.load(f) - sfed = StateEstimatorFederate( - federate_name, - parameters, - input_mapping - ) + sfed = StateEstimatorFederate(federate_name, parameters, input_mapping) sfed.run()