diff --git a/.github/workflows/add-to-project.yml b/.github/workflows/add-to-project.yml index 500f80d..998f17f 100644 --- a/.github/workflows/add-to-project.yml +++ b/.github/workflows/add-to-project.yml @@ -1,31 +1,31 @@ -name: Auto Assign to Project(s) +# name: Auto Assign to Project(s) -on: - issues: - types: [opened, labeled] - pull_request: - types: [opened, labeled] - issue_comment: - types: [created] -env: - MY_GITHUB_TOKEN: ${{ secrets.GRIDAPPSD_PAT }} +# on: +# issues: +# types: [opened, labeled] +# pull_request: +# types: [opened, labeled] +# issue_comment: +# types: [created] +# env: +# MY_GITHUB_TOKEN: ${{ secrets.GRIDAPPSD_PAT }} -jobs: - assign_one_project: - runs-on: ubuntu-latest - name: Assign to One Project - steps: - - name: Assign NEW issues and NEW pull requests to project 2 - uses: srggrs/assign-one-project-github-action@1.3.1 - if: github.event.action == 'opened' - with: - project: 'https://github.com/orgs/GRIDAPPSD/projects/7' +# jobs: +# assign_one_project: +# runs-on: ubuntu-latest +# name: Assign to One Project +# steps: +# - name: Assign NEW issues and NEW pull requests to project 2 +# uses: srggrs/assign-one-project-github-action@1.3.1 +# if: github.event.action == 'opened' +# with: +# project: 'https://github.com/orgs/GRIDAPPSD/projects/7' - - name: Assign issues and pull requests with `bug` label to project 3 - uses: srggrs/assign-one-project-github-action@1.3.1 - if: | - contains(github.event.issue.labels.*.name, 'bug') || - contains(github.event.pull_request.labels.*.name, 'bug') - with: - project: 'https://github.com/orgs/GRIDAPPSD/projects/7' - column_name: 'Labeled' +# - name: Assign issues and pull requests with `bug` label to project 3 +# uses: srggrs/assign-one-project-github-action@1.3.1 +# if: | +# contains(github.event.issue.labels.*.name, 'bug') || +# contains(github.event.pull_request.labels.*.name, 'bug') +# with: +# project: 'https://github.com/orgs/GRIDAPPSD/projects/7' +# column_name: 'Labeled' diff --git a/.github/workflows/auto-assign-project.yml b/.github/workflows/auto-assign-project.yml new file mode 100644 index 0000000..9cb9424 --- /dev/null +++ b/.github/workflows/auto-assign-project.yml @@ -0,0 +1,22 @@ +name: Add bugs to bugs project + +on: + issues: + types: + - opened + +jobs: + add-to-project: + name: Add issue to project + runs-on: ubuntu-latest + steps: + - uses: actions/add-to-project@v0.3.0 + with: + # You can target a repository in a different organization + # to the issue + project-url: https://github.com/orgs/GRIDAPPSD/projects/7 + # project-url: https://github.com/orgs//projects/ + # github-token: ${{ secrets.ADD_TO_PROJECT_PAT }} + github-token: ${{ secrets.AUTO_ASSIGN_PAT }} + # labeled: bug, needs-triage + # label-operator: OR diff --git a/.gitignore b/.gitignore index 7f4c4a4..f4c7230 100644 --- a/.gitignore +++ b/.gitignore @@ -31,6 +31,9 @@ wheels/ *.egg-info/ .installed.cfg *.egg +*output/** +*data/** +*plots/** # PyInstaller # Usually these files are written by a python script from a template diff --git a/micro-apps/carbon-management-app/carbon_management_app.py b/micro-apps/carbon-management-app/carbon_management_app.py new file mode 100644 index 0000000..f106cb0 --- /dev/null +++ b/micro-apps/carbon-management-app/carbon_management_app.py @@ -0,0 +1,779 @@ +import time +import os +import math +from cmath import exp +from pprint import pprint +# import networkx as nx +# import numpy as np +import cvxpy as cp +# from pandas import DataFrame +from typing import Any, Dict +import importlib +from datetime import datetime, timedelta, timezone +from dataclasses import dataclass, asdict +from argparse import ArgumentParser +import json +from typing import Dict +from cimgraph import utils +from tabulate import tabulate +from cimgraph.databases import ConnectionParameters, BlazegraphConnection +from cimgraph.models import FeederModel +import cimgraph.data_profile.rc4_2021 as cim +from cimgraph.data_profile.rc4_2021 import EnergyConsumer +from cimgraph.data_profile.rc4_2021 import BatteryUnit +from cimgraph.data_profile.rc4_2021 import PowerElectronicsConnection +from gridappsd import GridAPPSD, topics, DifferenceBuilder +from gridappsd.simulation import Simulation +from gridappsd.simulation import PowerSystemConfig +from gridappsd.simulation import SimulationConfig +from gridappsd.simulation import SimulationArgs +from gridappsd.simulation import ModelCreationConfig +from gridappsd import topics as t +import csv +from pathlib import Path + +# IEEE123_APPS = "_E3D03A27-B988-4D79-BFAB-F9D37FB289F7" +IEEE13 = "_49AD8E07-3BF9-A4E2-CB8F-C3722F837B62" +IEEE123 = "_C1C3E687-6FFD-C753-582B-632A27E28507" +IEEE123PV = "_A3BC35AA-01F6-478E-A7B1-8EA4598A685C" +IEEE123_APPS = "_F49D1288-9EC6-47DB-8769-57E2B6EDB124" +ieee123_apps_feeder_head_measurement_mrids = [ + "_903a0e85-6a11-4b8e-82cf-163376df7764", + "_6d82c9d9-5f99-4356-8a8d-acda9cd6f17b", + "_971765b5-da69-4ea8-a0cb-d2b9613b8ee3" +] +SOURCE_BUS = '150r' +OUT_DIR = "outputs" +ROOT = os.getcwd() + + +class CarbonManagementApp(object): + """Carbon Management Class + This class implements a centralized algorithm for carbon management in a grid by controlling batteries in the distribution model. + + Attributes: + sim_id: The simulation id that the instance will be perfoming control on. If the sim_id is None then the + application will assume it is performing control on the actual field devices. + Type: str. + Default: None. + gad_obj: An instatiated object that is connected to the gridappsd message bus usually this should be the same object which subscribes, but that + isn't required. + Type: GridAPPSD. + Default: None. + model_id: The mrid of the cim model to perform carbon management control on. + Type: str. + Default: None. + simulation: Simulation object of the identified cim model. + Type: Simulation. + Default: None. + network: An instatiated object of knowledge graph class for distribution feeder objects. + Type: FeederModel. + Default: None. + Methods: + initialize(headers:Dict[Any], message:Dict[Any]): This callback function will trigger the application to + read the model database to find all the measurements in the model as well as batteries in the sytem. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + on_measurement(headers:Dict[Any], message:Dict[Any]): This callback function will be used to capture the + measurements dictionary needed for control. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + optimize_battery(): This is the main function for controlling the battery. + """ + + def __init__(self, + gad_obj: GridAPPSD, + model_id: str, + network: FeederModel, + sim_id: str = None, + simulation: Simulation = None): + if not isinstance(gad_obj, GridAPPSD): + raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + if not isinstance(model_id, str): + raise TypeError(f'model_id must be a str type.') + if model_id is None or model_id == '': + raise ValueError(f'model_id must be a valid uuid.') + if not isinstance(sim_id, str) and sim_id is not None: + raise TypeError(f'sim_id must be a string type or {None}!') + if not isinstance(simulation, Simulation) and simulation is not None: + raise TypeError(f'The simulation arg must be a Simulation type or {None}!') + self.simulation = simulation + self.gad_obj = gad_obj + self.init_batt_dis = True + self._count = 0 + self._publish_to_topic = topics.simulation_input_topic(simulation.simulation_id) + self._init_batt_diff = DifferenceBuilder(simulation.simulation_id) + + self.Battery = {} + self.Solar = {} + self.EnergyConsumer = {} + self.peak_va_measurements_A = {} + self.peak_va_measurements_B = {} + self.peak_va_measurements_C = {} + self.has_batteries = True + self.has_power_electronics = True + self.has_energy_consumers = True + if PowerElectronicsConnection not in network.graph: + self.has_power_electronics = False + self.has_batteries = False + # raise ValueError("No power electronic devices in network.") + elif len(network.graph[PowerElectronicsConnection].keys()) == 0: + self.has_power_electronics = False + self.has_batteries = False + if EnergyConsumer not in network.graph: + self.has_energy_consumers = False + + if self.has_power_electronics: + self._collect_power_electronic_devices(network) + # raise ValueError("No power electronic devices in network.") + + if len(self.Battery) == 0: + self.has_batteries = False + print("No batteries in system.") + # raise ValueError("No batteries in network.") + + self.network = network + self.findFeederHeadLoadMeasurements() + + if self.has_energy_consumers: + self._collect_energy_consumers(network) + simulation.add_onmeasurement_callback(self.on_measurement) + # data output + out_dir = Path(__file__).parent/"output" + feeder = network.graph.get(cim.Feeder, {}).get(model_id) + self.model_results_path = out_dir/feeder.name + self.model_results_path.mkdir(parents=True, exist_ok=True) + # self.findFeederHeadLoadMeasurements() + self.init_output_file() + + def init_output_file(self): + for child in self.model_results_path.iterdir(): + if not child.is_dir(): + child.unlink() + with open(self.model_results_path / "voltages.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "voltage"] + writer.writerow(header) + with open(self.model_results_path / "total_load.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "phase", "p", "q"] + writer.writerow(header) + with open(self.model_results_path / "optimization_summary.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "load_pv_a", "load_pv_b", "load_pv_c", "load_pv_batt_a", "load_pv_batt_b", "load_pv_batt_c", "status"] + writer.writerow(header) + with open(self.model_results_path / "optimization_result.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "battery", "phases", "p_batt"] + writer.writerow(header) + with open(self.model_results_path / "simulation_table.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header=["time","mrid", "battery","phases","p_a","p_b","p_c","soc"] + writer.writerow(header) + with open(self.model_results_path / "dispatches.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "value"] + writer.writerow(header) + with open(self.model_results_path / "feeder_head.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "phase", "p", "q"] + writer.writerow(header) + + def get_feeder_head_measurements(self, measurements): + total_loads = {} + if self.network.container.mRID == IEEE123_APPS: + for mrid in ieee123_apps_feeder_head_measurement_mrids: + phase = self.network.graph.get(cim.Analog, {}).get(mrid).phases.value + s_magnitude = measurements.get(mrid).get("magnitude") + s_angle_rad = measurements.get(mrid).get("angle") + total_loads[phase] = s_magnitude*exp(1j*s_angle_rad) + return total_loads + + def save_total_load_data(self, time, total_loads): + file_path = self.model_results_path/"total_load.csv" + header = ["time", "phase", "p", "q"] + for phase, load in total_loads.items(): + data = [time, phase, load.real, load.imag] + add_data_to_csv(file_path, data, header=header) + + def save_voltage_data(self, time): + file_path = self.model_results_path/"voltages.csv" + header = ["time", "mrid", "node", "phase", "voltage"] + for mrid in self.pnv_measurements.keys(): + v = self.pnv_measurements_pu[mrid] + phase = self.pnv_measurements[mrid]["measurement_object"].phases.value + node = self.pnv_measurements[mrid]["measurement_object"].Terminal.ConnectivityNode.name + if v is not None: + data = [time, mrid, node, phase, v] + add_data_to_csv(file_path, data, header=header) + + def _collect_power_electronic_devices(self, network): + for pec in network.graph.get(PowerElectronicsConnection, {}).values(): + # inv_mrid = pec.mRID + for unit in pec.PowerElectronicsUnit: + unit_mrid = unit.mRID + if isinstance(unit, BatteryUnit): + self.Battery[unit_mrid] = {'phases': [], 'measurementType': [], 'measurementmRID': [], + 'measurementPhases': []} + self.Battery[unit_mrid]['name'] = unit.name + self.Battery[unit_mrid]['ratedS'] = float(pec.ratedS) / 1000 + self.Battery[unit_mrid]['ratedE'] = float(unit.ratedE) / 1000 + else: + self.Solar[unit_mrid] = {'phases': [], 'measurementType': [], 'measurementmRID': [], + 'measurementPhases': []} + self.Solar[unit_mrid]['name'] = pec.name + self.Solar[unit_mrid]['ratedS'] = float(pec.ratedS) / 1000 + + if not pec.PowerElectronicsConnectionPhases: + if unit_mrid in self.Battery: + self.Battery[unit_mrid]['phases'] = 'ABC' + if unit_mrid in self.Solar: + self.Solar[unit_mrid]['phases'] = 'ABC' + else: + phases = [] + for phase in pec.PowerElectronicsConnectionPhases: + phases.append(phase.phase.value) + if unit_mrid in self.Battery: + self.Battery[unit_mrid]['phases'] = phases + if unit_mrid in self.Solar: + self.Solar[unit_mrid]['phases'] = phases + + for measurement in pec.Measurements: + if unit_mrid in self.Battery: + self.Battery[unit_mrid]['measurementType'].append(measurement.measurementType) + self.Battery[unit_mrid]['measurementmRID'].append(measurement.mRID) + if measurement.phases.value is not None: + self.Battery[unit_mrid]['measurementPhases'].append(measurement.phases.value) + if unit_mrid in self.Solar: + self.Solar[unit_mrid]['measurementType'].append(measurement.measurementType) + self.Solar[unit_mrid]['measurementmRID'].append(measurement.mRID) + if measurement.phases.value is not None: + self.Solar[unit_mrid]['measurementPhases'].append(measurement.phases.value) + + def _collect_energy_consumers(self, network): + for ld in network.graph.get(EnergyConsumer, {}).values(): + ld_mrid = ld.mRID + self.EnergyConsumer[ld_mrid] = {'phases': [], 'measurementType': [], 'measurementmRID': [], + 'measurementPhases': []} + self.EnergyConsumer[ld_mrid]['name'] = ld.name + if not ld.EnergyConsumerPhase: + self.EnergyConsumer[ld_mrid]['phases'] = 'ABC' + else: + phases = [] + for phase in ld.EnergyConsumerPhase: + phases.append(phase.phase.value) + self.EnergyConsumer[ld_mrid]['phases'] = phases + for measurement in ld.Measurements: + self.EnergyConsumer[ld_mrid]['measurementType'].append(measurement.measurementType) + self.EnergyConsumer[ld_mrid]['measurementmRID'].append(measurement.mRID) + if measurement.phases.value is not None: + self.EnergyConsumer[ld_mrid]['measurementPhases'].append(measurement.phases.value) + + def find_phase(self, degree): + ref_points = [0, 120, -120] + closest = min(ref_points, key=lambda x: abs(degree - x)) + phases = {0: 'A', -120: 'B', 120: 'C'} + return phases[closest] + + def pol2cart(self, mag, angle_deg): + # Convert degrees to radians. GridAPPS-D spits angle in degrees + angle_rad = math.radians(angle_deg) + p = mag * math.cos(angle_rad) + q = mag * math.sin(angle_rad) + return p, q + + def find_injection(self, object, measurements): + + for item in object: + object[item]['P_inj'] = [0, 0, 0] + object[item]['Q_inj'] = [0, 0, 0] + meas_type = object[item]['measurementType'] + va_idx = [i for i in range(len(meas_type)) if meas_type[i] == 'VA'] + pnv_idx = [i for i in range(len(meas_type)) if meas_type[i] == 'PNV'] + soc_idx = [i for i in range(len(meas_type)) if meas_type[i] == 'SoC'] + if soc_idx: + soc = measurements[object[item]['measurementmRID'][soc_idx[0]]]['value'] + self.Battery[item]['soc'] = soc + if 's' in object[item]['phases'][0]: + angle = measurements[object[item]['measurementmRID'][pnv_idx[0]]]['angle'] + rho = measurements[object[item]['measurementmRID'][va_idx[0]]]['magnitude'] + phi = measurements[object[item]['measurementmRID'][va_idx[0]]]['angle'] + p, q = self.pol2cart(rho, phi) + if self.find_phase(angle) == 'A': + object[item]['P_inj'][0] = 2 * p / 1000 + object[item]['Q_inj'][0] = 2 * q / 1000 + object[item]['phases'] = 'A' + elif self.find_phase(angle) == 'B': + object[item]['P_inj'][1] = 2 * p / 1000 + object[item]['Q_inj'][1] = 2 * q / 1000 + object[item]['phases'] = 'B' + else: + object[item]['P_inj'][2] = 2 * p / 1000 + object[item]['Q_inj'][2] = 2 * q / 1000 + object[item]['phases'] = 'C' + elif object[item]['phases'] == 'ABC': + for k in range(3): + rho = measurements[object[item]['measurementmRID'][va_idx[k]]]['magnitude'] + phi = measurements[object[item]['measurementmRID'][va_idx[k]]]['angle'] + p, q = self.pol2cart(rho, phi) + object[item]['P_inj'][k] = p / 1000 + object[item]['Q_inj'][k] = q / 1000 + else: + rho = measurements[object[item]['measurementmRID'][va_idx[0]]]['magnitude'] + phi = measurements[object[item]['measurementmRID'][va_idx[0]]]['angle'] + p, q = self.pol2cart(rho, phi) + if object[item]['phases'][0] == 'A': + object[item]['P_inj'][0] = p / 1000 + object[item]['Q_inj'][0] = q / 1000 + elif object[item]['phases'][0] == 'B': + object[item]['P_inj'][1] = p / 1000 + object[item]['Q_inj'][1] = q / 1000 + else: + object[item]['P_inj'][2] = p / 1000 + object[item]['Q_inj'][2] = q / 1000 + + def optimize_battery(self, timestamp): + # Define optimization variables + n_batt = len(self.Battery) + p_flow_A = cp.Variable(integer=False, name='p_flow_A') + p_flow_B = cp.Variable(integer=False, name='p_flow_B') + p_flow_C = cp.Variable(integer=False, name='p_flow_C') + p_flow_mod_A = cp.Variable(integer=False, name='p_flow_mod_A') + p_flow_mod_B = cp.Variable(integer=False, name='p_flow_mod_B') + p_flow_mod_C = cp.Variable(integer=False, name='p_flow_mod_C') + p_batt = cp.Variable(n_batt, integer=False, name='p_batt') + soc = cp.Variable(n_batt, integer=False, name='soc') + lambda_c = cp.Variable(n_batt, boolean=True, name='lambda_c') + lambda_d = cp.Variable(n_batt, boolean=True, name='lambda_d') + P_batt_A = cp.Variable(n_batt, integer=False, name='P_batt_A') + P_batt_B = cp.Variable(n_batt, integer=False, name='P_batt_B') + P_batt_C = cp.Variable(n_batt, integer=False, name='P_batt_C') + deltaT = 0.25 + n_batt_ABC = 0 + for batt in self.Battery: + if 'ABC' in self.Battery[batt]['phases']: + n_batt_ABC += 1 + # Defining variable to constraint same sign for multiple 3-ph batteries + if n_batt_ABC > 0: + b = cp.Variable(n_batt_ABC, boolean=True, name='b') + + constraints = [] + sum_flow_A, sum_flow_B, sum_flow_C = 0.0, 0.0, 0.0 + for pv in self.Solar: + sum_flow_A += self.Solar[pv]['P_inj'][0] + sum_flow_B += self.Solar[pv]['P_inj'][1] + sum_flow_C += self.Solar[pv]['P_inj'][2] + for load in self.EnergyConsumer: + sum_flow_A += self.EnergyConsumer[load]['P_inj'][0] + sum_flow_B += self.EnergyConsumer[load]['P_inj'][1] + sum_flow_C += self.EnergyConsumer[load]['P_inj'][2] + + idx = 0 + idx_b = 0 + # For now, we assume battery to be 100% efficient. Need to rewrite the soc constraints if using different efficiency + for batt in self.Battery: + constraints.append( + soc[idx] == self.Battery[batt]['soc'] / 100 + p_batt[idx] * deltaT / self.Battery[batt]['ratedE']) + constraints.append(p_batt[idx] <= lambda_c[idx] * self.Battery[batt]['ratedS']) + constraints.append(p_batt[idx] >= - lambda_d[idx] * self.Battery[batt]['ratedS']) + constraints.append(lambda_c[idx] + lambda_d[idx] <= 1) + constraints.append(soc[idx] <= 0.9) + constraints.append(soc[idx] >= 0.2) + # Tracking dispatch in A, B, and C phase by different batteries + if 'ABC' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == p_batt[idx] / 3) + constraints.append(P_batt_B[idx] == p_batt[idx] / 3) + constraints.append(P_batt_C[idx] == p_batt[idx] / 3) + constraints.append(p_batt[idx] >= - b[idx_b] * self.Battery[batt]['ratedS']) + constraints.append(p_batt[idx] <= self.Battery[batt]['ratedS'] * (1 - b[idx_b])) + idx_b += 1 + else: + if 'A' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == p_batt[idx]) + constraints.append(P_batt_B[idx] == 0.0) + constraints.append(P_batt_C[idx] == 0.0) + if 'B' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == 0.0) + constraints.append(P_batt_B[idx] == p_batt[idx]) + constraints.append(P_batt_C[idx] == 0.0) + if 'C' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == 0.0) + constraints.append(P_batt_B[idx] == 0.0) + constraints.append(P_batt_C[idx] == p_batt[idx]) + idx += 1 + # Ensuring three phase batteries have the same sign. We don't want one charging and another discharging + # Although, mathematically it might sound correct, it is not worth to do such dispatch. + for k in range(n_batt_ABC - 1): + constraints.append(b[k] == b[k + 1]) + + # Constraints for flow in phase ABC + constraints.append(p_flow_A == sum_flow_A + sum(P_batt_A[k] for k in range(n_batt))) + constraints.append(p_flow_B == sum_flow_B + sum(P_batt_B[k] for k in range(n_batt))) + constraints.append(p_flow_C == sum_flow_C + sum(P_batt_C[k] for k in range(n_batt))) + + # Individual modulus is needed to make sure one phase doesn't compensate for other + constraints.append(p_flow_mod_A >= p_flow_A) + constraints.append(p_flow_mod_A >= - p_flow_A) + + constraints.append(p_flow_mod_B >= p_flow_B) + constraints.append(p_flow_mod_B >= - p_flow_B) + + constraints.append(p_flow_mod_C >= p_flow_C) + constraints.append(p_flow_mod_C >= - p_flow_C) + + # Objective function and invoke a solver + objective = (p_flow_mod_A + p_flow_mod_B + p_flow_mod_C) + problem = cp.Problem(cp.Minimize(objective), constraints) + problem.solve() + + # Extract optimization solution + idx = 0 + dispatch_batteries = {} + optimization_solution_table = [] + for batt in self.Battery: + name = self.Battery[batt]['name'] + dispatch_batteries[batt] = {} + if p_batt[idx].value is None: + continue + dispatch_batteries[batt]['p_batt'] = p_batt[idx].value * 1000 + optimization_solution_table.append([name, self.Battery[batt]['phases'], p_batt[idx].value]) + data = [timestamp, batt, name, self.Battery[batt]['phases'], p_batt[idx].value] + header = ["time", "battery", "phases", "p_batt"] + add_data_to_csv(self.model_results_path/"optimization_result.csv", data, header=header) + idx += 1 + print('Optimization Solution') + print(tabulate(optimization_solution_table, headers=['Battery', 'phases', 'P_batt (kW)'], tablefmt='psql')) + if problem.status == "optimal": + load_pv = ['{:.3f}'.format(sum_flow_A), '{:.3f}'.format(sum_flow_B), '{:.3f}'.format(sum_flow_C)] + load_pv_batt = ['{:.3f}'.format(p_flow_A.value), '{:.3f}'.format(p_flow_B.value), + '{:.3f}'.format(p_flow_C.value)] + optimization_summary = [] + optimization_summary.append([load_pv, load_pv_batt, problem.status]) + data = [timestamp] + data.extend(load_pv) + data.extend(load_pv_batt) + data.append(problem.status) + header = ["time", "load_pv_a", "load_pv_b", "load_pv_c", "load_pv_batt_a", "load_pv_batt_b", "load_pv_batt_c", "status"] + add_data_to_csv(self.model_results_path/"optimization_summary.csv", data, header=header) + print(tabulate(optimization_summary, headers=['Load+PV (kW)', 'Load+PV+Batt (kW)', 'Status'], tablefmt='psql')) + return dispatch_batteries + + def findFeederHeadLoadMeasurements(self): + energySources = self.network.graph.get(cim.EnergySource, {}) + feederLoadObject = None + for eSource in energySources.values(): + feederLoadObject = findFeederHeadPowerMeasurmentsObject(eSource) + feederLoadMeasurements = feederLoadObject.Measurements + for measurement in feederLoadMeasurements: + if measurement.measurementType == 'VA': + if measurement.phases in [cim.PhaseCode.A, cim.PhaseCode.AN]: + self.peak_va_measurements_A[measurement.mRID] = {'object': measurement, 'value': None} + elif measurement.phases in [cim.PhaseCode.B, cim.PhaseCode.BN]: + self.peak_va_measurements_B[measurement.mRID] = {'object': measurement, 'value': None} + elif measurement.phases in [cim.PhaseCode.C, cim.PhaseCode.CN]: + self.peak_va_measurements_C[measurement.mRID] = {'object': measurement, 'value': None} + if not self.peak_va_measurements_A or not self.peak_va_measurements_B or not self.peak_va_measurements_C: + raise RuntimeError(f'feeder {self.graph_model.container.mRID}, has no measurements associated with the ' + 'feeder head transformer!') + + def on_measurement(self, sim: Simulation, timestamp: dict, measurements: dict) -> None: + if self.simulation is not None: + self.simulation.pause() + self.timestamp = timestamp + total_loads = self.get_feeder_head_measurements(measurements) + self.save_total_load_data(timestamp, total_loads) + for mrid in self.peak_va_measurements_A.keys(): + measurement = measurements.get(self.peak_va_measurements_A[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_A[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "A", load.real, load.imag]) + for mrid in self.peak_va_measurements_B.keys(): + measurement = measurements.get(self.peak_va_measurements_B[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_B[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "B", load.real, load.imag]) + for mrid in self.peak_va_measurements_C.keys(): + measurement = measurements.get(self.peak_va_measurements_C[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_C[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "C", load.real, load.imag]) + if not self.has_batteries: + return + if self._count % 10 == 0: + # Call function to extract simulation measurements from injection sources + self.find_injection(self.Battery, measurements) + self.find_injection(self.EnergyConsumer, measurements) + self.find_injection(self.Solar, measurements) + + # Display real-time battery injection and current SOC from measurements + simulation_table_batteries = [] + for batt in self.Battery: + name = self.Battery[batt]['name'] + phases = self.Battery[batt]['phases'] + data = [timestamp, batt, name, phases] + simulation_table_batteries.append( + [name, phases, self.Battery[batt]['P_inj'], self.Battery[batt]['soc']]) + data.extend(self.Battery[batt]['P_inj']) + data.extend([self.Battery[batt]['soc']]) + header=["time","mrid", "battery","phases","p_a","p_b","p_c","soc"] + add_data_to_csv(self.model_results_path/"simulation_table.csv", data, header=header) + print(f'\n.......Curren timestamp: {timestamp}.......\n') + print('Simulation Table') + print(tabulate(simulation_table_batteries, headers=['Battery', 'phases', 'P_batt (kW)', 'SOC'], + tablefmt='psql')) + + # Invoke optimization for given grid condition + dispatch_values = self.optimize_battery(timestamp) + + # Dispatch battery. Note that -ve values in forward difference means charging batteries + # Make necessary changes in sign convention from optimization values + for unit in dispatch_values: + self._init_batt_diff.add_difference(unit, 'PowerElectronicsConnection.p', + - dispatch_values[unit]['p_batt'], 0.0) + msg = self._init_batt_diff.get_message() + with open(self.model_results_path/f"dispatch_{timestamp}.json", "w") as f: + json.dump(msg, f, indent=4) + + header = ["time", "mrid", "value"] + data = [timestamp, unit, -dispatch_values[unit]['p_batt']] + add_data_to_csv(self.model_results_path/"dispatches.csv", data) + self.gad_obj.send(self._publish_to_topic, json.dumps(msg)) + + self._count += 1 + if self.simulation is not None: + self.simulation.resume() + + +def findFeederHeadPowerMeasurmentsObject(cimObj: object): + ''' + Helper function to find feeder transformer + ''' + if not isinstance(cimObj, cim.EnergySource): + raise TypeError('findFeederHeadPowerMeasurmentsObject(): cimObj must be an instance of cim.EnergySource!') + equipmentToCheck = [cimObj] + feederHeadLoadMeasurementsObject = None + i = 0 + while not feederHeadLoadMeasurementsObject and i < len(equipmentToCheck): + equipmentToAdd = [] + for eq in equipmentToCheck[i:]: + powerMeasurementsCount = 0 + for meas in eq.Measurements: + if meas.measurementType == 'VA': + powerMeasurementsCount += 1 + if powerMeasurementsCount == 3: + feederHeadLoadMeasurementsObject = eq + break + else: + terminals = eq.Terminals + connectivityNodes = [] + for t in terminals: + if t.ConnectivityNode not in connectivityNodes: + connectivityNodes.append(t.ConnectivityNode) + for cn in connectivityNodes: + for t in cn.Terminals: + if t.ConductingEquipment not in equipmentToCheck and t.ConductingEquipment not in equipmentToAdd: + equipmentToAdd.append(t.ConductingEquipment) + i = len(equipmentToCheck) + equipmentToCheck.extend(equipmentToAdd) + if not feederHeadLoadMeasurementsObject: + raise RuntimeError('findFeederPowerRating(): No feeder head object with power measurements could be found for ' + f'EnergySource, {cimObj.name}!') + return feederHeadLoadMeasurementsObject + + +def createSimulation(gad_obj: GridAPPSD, model_info: Dict[str, Any]) -> Simulation: + # if not isinstance(gad_obj, GridAPPSD): + # raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + # if not isinstance(model_info, dict): + # raise TypeError(f'The model_id must be a dictionary type.') + line_name = model_info.get('modelId') + subregion_name = model_info.get('subRegionId') + region_name = model_info.get('regionId') + sim_name = model_info.get('modelName') + # if line_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelId.') + # if subregion_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key subRegionId.') + # if region_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key regionId.') + # if sim_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelName.') + # psc = PowerSystemConfig(Line_name=line_name, + # SubGeographicalRegion_name=subregion_name, + # GeographicalRegion_name=region_name) + # start_time = int(datetime.utcnow().replace(microsecond=0).timestamp()) + # sim_args = SimulationArgs( + # start_time=f'{start_time}', + # # duration=f'{24*3600}', + # duration=120, + # simulation_name=sim_name, + # run_realtime=False + # # pause_after_measurements=True + # ) + # sim_config = SimulationConfig(power_system_config=psc, simulation_config=sim_args) + # sim_obj = Simulation(gapps=gad_obj, run_config=sim_config) + system_config = PowerSystemConfig( + GeographicalRegion_name=region_name, + SubGeographicalRegion_name=subregion_name, + Line_name=line_name + ) + + model_config = ModelCreationConfig( + load_scaling_factor=1, + schedule_name="ieeezipload", + z_fraction=0, + i_fraction=1, + p_fraction=0, + randomize_zipload_fractions=False, + use_houses=False + ) + + start = datetime(2023, 1, 1, 0, tzinfo=timezone.utc) + epoch = datetime.timestamp(start) + duration = timedelta(hours=24).total_seconds() + # using hadrcode epoch to ensure afternoon time + # epoch = 1704207000 + + sim_args = SimulationArgs( + start_time=epoch, + duration=duration, + simulator="GridLAB-D", + timestep_frequency=1000, + timestep_increment=1000, + run_realtime=False, + simulation_name=sim_name, + power_flow_solver_method="NR", + model_creation_config=asdict(model_config), + pause_after_measurements=True, + ) + + sim_config = SimulationConfig( + power_system_config=asdict(system_config), + simulation_config=asdict(sim_args) + ) + + sim_obj = Simulation(gad_obj, asdict(sim_config)) + return sim_obj + + +def main(control_enabled: bool, start_simulations: bool, model_id: str = None): + if not isinstance(control_enabled, bool): + raise TypeError(f'Argument, control_enabled, must be a boolean.') + if not isinstance(start_simulations, bool): + raise TypeError(f'Argument, start_simulations, must be a boolean.') + if not isinstance(model_id, str) and model_id is not None: + raise TypeError( + f'The model id passed to the convervation voltage reduction application must be a string type or {None}.') + + cim_profile = 'rc4_2021' + cim = importlib.import_module('cimgraph.data_profile.' + cim_profile) + feeder = cim.Feeder(mRID=model_id) + params = ConnectionParameters( + url="http://localhost:8889/bigdata/namespace/kb/sparql", + cim_profile=cim_profile) + + gapps = GridAPPSD(username='system', password='manager') + bg = BlazegraphConnection(params) + network = FeederModel( + connection=bg, + container=feeder, + distributed=False) + utils.get_all_data(network) + + local_simulations = {} + app_instances = {'field_instances': {}, 'external_simulation_instances': {}, 'local_simulation_instances': {}} + response = gapps.query_model_info() + models = response.get('data', {}).get('models', []) + model_is_valid = False + if model_id is not None: + for m in models: + m_id = m.get('modelId') + if model_id == m_id: + model_is_valid = True + break + if not model_is_valid: + raise ValueError(f'The model id provided does not exist in the GridAPPS-D plaform.') + else: + models = [m] + if start_simulations: + for m in models: + local_simulations[m.get('modelId', '')] = createSimulation(gapps, m) + else: + # TODO: query platform for running simulations which is currently not implemented in the GridAPPS-D Api + pass + # Create an cvr controller instance for all the real systems in the database + # for m in models: + # m_id = m.get('modelId') + # app_instances['field_instances'][m_id] = ConservationVoltageReductionController(gad_object, m_id) + + for m_id, simulation in local_simulations.items(): + measurements_topic = topics.simulation_output_topic(simulation) + simulation.start_simulation() + c_mapp = CarbonManagementApp( + gapps, m_id, network, simulation=simulation) + # gapps.subscribe(measurements_topic, c_mapp) + try: + while True: + time.sleep(0.1) + except KeyboardInterrupt: + print("Simulation manually stopped.") + finally: + print(" -- Stopping simulation -- ") + simulation.stop() + + +def add_data_to_csv(file_path, data, header=None): + # Check if the file exists + file_exists = os.path.isfile(file_path) + + # Open the CSV file in append mode (or create it if it doesn't exist) + with open(file_path, mode='a', newline='') as file: + writer = csv.writer(file) + # If the file doesn't exist, write the header first + if not file_exists and header: + writer.writerow(header) + # Write the data (a list or tuple representing a row) + writer.writerow(data) + + +if __name__ == "__main__": + parser = ArgumentParser() + parser.add_argument('model_id', nargs='?', default=None, help='The unit_mrid of the cim model to perform cvr on.') + parser.add_argument('-s', + '--start_simulations', + action='store_true', + help='Flag to have application start simulations') + parser.add_argument('-d', + '--disable_control', + action='store_true', + help='Flag to disable control on startup by default.') + args = parser.parse_args() + # args.start_simulations = True + # args.model_id = '_EE71F6C9-56F0-4167-A14E-7F4C71F10EAA' + main(args.disable_control, args.start_simulations, args.model_id) diff --git a/micro-apps/conservation-voltage-reduction-app/cvr_main.py b/micro-apps/conservation-voltage-reduction-app/cvr_main.py new file mode 100644 index 0000000..f8a42df --- /dev/null +++ b/micro-apps/conservation-voltage-reduction-app/cvr_main.py @@ -0,0 +1,1352 @@ +import importlib +import logging +import math +import os +import csv +from pathlib import Path +from cmath import exp +from argparse import ArgumentParser +from datetime import datetime, timezone, timedelta +from typing import Any, Dict +from uuid import uuid4 + +import cimgraph.utils as cimUtils +from cimgraph.databases import ConnectionParameters +from cimgraph.databases.blazegraph.blazegraph import BlazegraphConnection +from cimgraph.databases.gridappsd.gridappsd import GridappsdConnection +from cimgraph.models import FeederModel +from gridappsd import DifferenceBuilder, GridAPPSD, topics +from gridappsd.simulation import * +from gridappsd.utils import ProcessStatusEnum +from opendssdirect import dss + +DB_CONNECTION = None +CIM_GRAPH_MODELS = {} +ieee123_apps_mrid = "_F49D1288-9EC6-47DB-8769-57E2B6EDB124" +ieee123_apps_feeder_head_measurement_mrids = [ + "_903a0e85-6a11-4b8e-82cf-163376df7764", + "_6d82c9d9-5f99-4356-8a8d-acda9cd6f17b", + "_971765b5-da69-4ea8-a0cb-d2b9613b8ee3" +] + +logging.basicConfig(format='%(asctime)s::%(levelname)s::%(name)s::%(filename)s::%(lineno)d::%(message)s', + filename='cvr.log', + filemode='w', + level=logging.INFO, + encoding='utf-8') +logger = logging.getLogger(__name__) +cim = None +DB_CONNECTION = None +CIM_GRAPH_MODELS = {} + + +class ConservationVoltageReductionController(object): + """CVR Control Class + This class implements a centralized algorithm for conservation voltage reduction by controlling regulators and + capacitors in the distribution model. + + Attributes: + sim_id: The simulation id that the instance will be perfoming cvr on. If the sim_id is None then the + application will assume it is performing cvr on the actual field devices. + Type: str. + Default: None. + period: The amount of time to wait before calculating new cvr setpoints. + Type: int. + Default: None. + Methods: + initialize(headers:Dict[Any], message:Dict[Any]): This callback function will trigger the application to + read the model database to find all the measurements in the model as well as the controllable + regulators and capacitors in the sytem. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + enableControl(headers:Dict[Any], message:Dict[Any]): This callback function will tell the application that + it is allowed to perform cvr on the systems or simulations its tied to. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + disableControl(headers:Dict[Any], message:Dict[Any]): This callback function will prevent the application + from performing cvr on the systems or simulations its tied to. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + on_measurement(headers:Dict[Any], message:Dict[Any]): This callback function will be used to update the + applications measurements dictionary needed for cvr control. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + cvr_control(): This is the main function for performing the cvr control. + """ + period = 3600 + lower_voltage_limit_pu = 0.9 + max_violation_time = 300 + + def __init__(self, + gad_obj: GridAPPSD, + model_id: str, + period: int = None, + low_volt_lim: float = None, + sim_id: str = None, + simulation: Simulation = None): + if not isinstance(gad_obj, GridAPPSD): + raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + if not isinstance(model_id, str): + raise TypeError(f'model_id must be a str type.') + if model_id is None or model_id == '': + raise ValueError(f'model_id must be a valid uuid.') + if not isinstance(period, int) and period is not None: + raise TypeError(f'period must be an int type or {None}!') + if not isinstance(low_volt_lim, float) and period is not None: + raise TypeError(f'low_volt_lim must be an float type or {None}!') + if not isinstance(sim_id, str) and sim_id is not None: + raise TypeError(f'sim_id must be a string type or {None}!') + if not isinstance(simulation, Simulation) and simulation is not None: + raise TypeError(f'The simulation arg must be a Simulation type or {None}!') + self.model_results_path = None + self.timestamp = 0 + self.model_id = model_id + self.id = uuid4() + self.platform_measurements = {} + self.desired_setpoints = {} + self.controllable_regulators = {} + self.controllable_capacitors = {} + self.pnv_measurements = {} + self.pnv_measurements_pu = {} + self.va_measurements = {} + self.pos_measurements = {} + self.link_measurements = {} + self.source_measurements = {} + self.peak_va_measurements_A = {} + self.peak_va_measurements_B = {} + self.peak_va_measurements_C = {} + if period is None: + self.period = ConservationVoltageReductionController.period + else: + self.period = period + if low_volt_lim is None: + self.low_volt_lim = ConservationVoltageReductionController.lower_voltage_limit_pu + else: + self.low_volt_lim = low_volt_lim + self.measurements_topic = None + self.setpoints_topic = None + self.simulation = None + self.simulation_id = None #TODO: figure out what simulation_id should be when deployed in the field. + if simulation is not None: + self.simulation = simulation + self.simulation_id = simulation.simulation_id + self.simulation.add_onmeasurement_callback(self.on_measurement) + if sim_id is None: + self.measurements_topic = topics.field_output_topic(None, None) + else: + self.simulation_id = sim_id + self.measurements_topic = topics.simulation_output_topic(sim_id) + self.gad_obj = gad_obj + self.log = self.gad_obj.get_logger() + if self.measurements_topic: + self.gad_obj.subscribe(self.measurements_topic, self.on_measurement_callback) + if self.simulation_id: + self.setpoints_topic = topics.simulation_input_topic(self.simulation_id) + else: + self.setpoints_topic = topics.field_input_topic() + # Read model_id from cimgraph to get all the controllable regulators and capacitors, and measurements. + self.graph_model = buildGraphModel(model_id) + self.findFeederHeadLoadMeasurements() + # Store controllable_capacitors + self.controllable_capacitors = self.graph_model.graph.get(cim.LinearShuntCompensator, {}) + # print(self.controllable_capacitors.keys()) + # Store controllable_regulators + powerTransformers = self.graph_model.graph.get(cim.PowerTransformer, {}) + ratioTapChangers = self.graph_model.graph.get(cim.RatioTapChanger, {}) + # print(powerTransformers.keys()) + # print(ratioTapChangers.keys()) + for mRID, powerTransformer in powerTransformers.items(): + for powerTransformerEnd in powerTransformer.PowerTransformerEnd: + if powerTransformerEnd.RatioTapChanger is not None: + if mRID not in self.controllable_regulators.keys(): + self.controllable_regulators[mRID] = {} + self.controllable_regulators[mRID]['RatioTapChangers'] = [] + self.controllable_regulators[mRID]['PhasesToName'] = {} + self.controllable_regulators[mRID]['regulator_object'] = powerTransformer + self.controllable_regulators[mRID]['RatioTapChangers'].append(powerTransformerEnd.RatioTapChanger) + self.controllable_regulators[mRID]['PhasesToName'] = { + 'A': powerTransformer.name, + 'B': powerTransformer.name, + 'C': powerTransformer.name + } + for transformerTank in powerTransformer.TransformerTanks: + for transformerEnd in transformerTank.TransformerTankEnds: + if transformerEnd.RatioTapChanger is not None: + if mRID not in self.controllable_regulators.keys(): + self.controllable_regulators[mRID] = {} + self.controllable_regulators[mRID]['RatioTapChangers'] = [] + self.controllable_regulators[mRID]['PhasesToName'] = {} + self.controllable_regulators[mRID]['regulator_object'] = powerTransformer + self.controllable_regulators[mRID]['RatioTapChangers'].append(transformerEnd.RatioTapChanger) + self.controllable_regulators[mRID]['PhasesToName'][transformerEnd.phases.value] = \ + transformerTank.name + + # Store measurements of voltages, loads, pv, battery, capacitor status, regulator taps, switch states. + # print(self.controllable_regulators.keys()) + measurements = {} + measurements.update(self.graph_model.graph.get(cim.Analog, {})) + measurements.update(self.graph_model.graph.get(cim.Discrete, {})) + + # print(measurements.keys()) + for meas in measurements.values(): + if meas.measurementType == 'PNV': # it's a voltage measurement. store it. + mrid = meas.mRID + self.pnv_measurements[mrid] = {'measurement_object': meas, 'measurement_value': None} + elif meas.measurementType == 'VA': # it's a power measurement. + if isinstance(meas.PowerSystemResource, (cim.EnergyConsumer, cim.PowerElectronicsConnection)): + mrid = meas.PowerSystemResource.mRID + if mrid not in self.va_measurements.keys(): + self.va_measurements[mrid] = {'measurement_objects': {}, 'measurement_values': {}} + self.va_measurements[mrid]['measurement_objects'][meas.mRID] = meas + self.va_measurements[mrid]['measurement_values'][meas.mRID] = None + if isinstance(meas.PowerSystemResource, (cim.ACLineSegment, cim.PowerTransformer)): + mrid = meas.PowerSystemResource.mRID + if mrid not in self.va_measurements.keys(): + self.link_measurements[mrid] = {'measurement_objects': {}, 'measurement_values': {}} + self.link_measurements[mrid]['measurement_objects'][meas.mRID] = meas + self.link_measurements[mrid]['measurement_values'][meas.mRID] = None + elif meas.measurementType == 'Pos': # it's a positional measurement. + if isinstance(meas.PowerSystemResource, (cim.Switch, cim.PowerTransformer, cim.LinearShuntCompensator)): + mrid = meas.PowerSystemResource.mRID + if mrid not in self.pos_measurements.keys(): + self.pos_measurements[mrid] = {'measurement_objects': {}, 'measurement_values': {}} + self.pos_measurements[mrid]['measurement_objects'][meas.mRID] = meas + self.pos_measurements[mrid]['measurement_values'][meas.mRID] = None + # for r in self.controllable_regulators.values(): + # print(json.dumps(r['PhasesToName'], indent=4, sort_keys=True)) + # data output + out_dir = Path(__file__).parent/"output" + feeder = self.graph_model.graph.get(cim.Feeder, {}).get(self.model_id) + self.model_results_path = out_dir/feeder.name + logger.info(self.model_results_path) + self.model_results_path.mkdir(parents=True, exist_ok=True) + # self.findFeederHeadLoadMeasurements() + self.init_output_file() + if self.simulation is not None: + self.simulation.start_simulation() + self.simulation_id = self.simulation.simulation_id + self.differenceBuilder = DifferenceBuilder(self.simulation_id) + self.dssContext = dss.NewContext() + self.create_opendss_context() + self.next_control_time = 0 + self.voltage_violation_time = -1 + self.isValid = True + self.isFirstRun = True + + def init_output_file(self): + for child in self.model_results_path.iterdir(): + child.unlink() + with open(self.model_results_path / "voltages.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "voltage"] + writer.writerow(header) + with open(self.model_results_path / "total_load.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "phase", "p", "q"] + writer.writerow(header) + with open(self.model_results_path / "caps.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "status"] + writer.writerow(header) + with open(self.model_results_path / "regs.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "tap"] + writer.writerow(header) + with open(self.model_results_path / "feeder_head.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "phase", "p", "q"] + writer.writerow(header) + with open(self.model_results_path / "caps_meas.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "status"] + writer.writerow(header) + with open(self.model_results_path / "regs_meas.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "tap"] + writer.writerow(header) + + def get_feeder_head_measurements(self, measurements): + total_loads = {} + if self.graph_model.container.mRID == ieee123_apps_mrid: + for mrid in ieee123_apps_feeder_head_measurement_mrids: + phase = self.graph_model.graph.get(cim.Analog, {}).get(mrid).phases.value + s_magnitude = measurements.get(mrid).get("magnitude") + s_angle_rad = measurements.get(mrid).get("angle") + total_loads[phase] = s_magnitude*exp(1j*math.radians(s_angle_rad)) + return total_loads + + def on_measurement(self, sim: Simulation, timestamp: str, measurements: Dict[str, Dict]): + self.timestamp = timestamp + self.desired_setpoints.clear() + if isinstance(sim, Simulation): + self.setpoints_topic = topics.simulation_input_topic(sim.simulation_id) + self.simulation.pause() + elif self.simulation_id: + self.setpoints_topic = topics.simulation_input_topic(self.simulation_id) + else: + self.setpoints_topic = topics.field_input_topic() + + + for mrid in self.peak_va_measurements_A.keys(): + measurement = measurements.get(self.peak_va_measurements_A[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_A[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "A", load.real, load.imag]) + for mrid in self.peak_va_measurements_B.keys(): + measurement = measurements.get(self.peak_va_measurements_B[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_B[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "B", load.real, load.imag]) + for mrid in self.peak_va_measurements_C.keys(): + measurement = measurements.get(self.peak_va_measurements_C[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_C[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "C", load.real, load.imag]) + + for mrid in self.pnv_measurements.keys(): + meas = measurements.get(mrid) + if meas is not None: + self.pnv_measurements[mrid]['measurement_value'] = meas + for psr_mrid in self.va_measurements.keys(): + for mrid in self.va_measurements[psr_mrid]['measurement_values'].keys(): + meas = measurements.get(mrid) + if meas is not None: + self.va_measurements[psr_mrid]['measurement_values'][mrid] = meas + for psr_mrid in self.pos_measurements.keys(): + for mrid in self.pos_measurements[psr_mrid]['measurement_values'].keys(): + meas = measurements.get(mrid) + if meas is not None: + self.pos_measurements[psr_mrid]['measurement_values'][mrid] = meas + if isinstance(self.pos_measurements[psr_mrid]['measurement_objects'][mrid].PowerSystemResource, (cim.LinearShuntCompensator, cim.PowerTransformer)): + pos_measurement_debug_data = { + "object_type": f"{type(self.pos_measurements[psr_mrid]['measurement_objects'][mrid].PowerSystemResource).__name__}", + "name": self.pos_measurements[psr_mrid]['measurement_objects'][mrid].PowerSystemResource.name, + "phase": self.pos_measurements[psr_mrid]['measurement_objects'][mrid].phases.value, + "value": int(meas.get("value")), + "timestamp": timestamp + } + logger.info(f"Pos Measurement: {json.dumps(pos_measurement_debug_data)}") + self.calculate_per_unit_voltage() + self.save_voltage_data(timestamp) + self.update_opendss_with_measurements() + if int(timestamp) >= self.next_control_time: + self.cvr_control() + self.next_control_time = int(timestamp) + self.period + self.voltage_violation_time = -1 + else: + total_violation_time = 0 + voltage_violation_phases = '' + for mrid, val in self.pnv_measurements_pu.items(): + if val is not None: + if val < self.low_volt_lim: + if self.voltage_violation_time < 0: + self.voltage_violation_time = int(timestamp) + meas = self.pnv_measurements.get(mrid, {}).get('measurement_object', None) + if isinstance(meas, cim.Measurement): + measPhase = meas.phases + if measPhase in [ + cim.PhaseCode.s1, cim.PhaseCode.s12, cim.PhaseCode.s1N, cim.PhaseCode.s12N, + cim.PhaseCode.s2, cim.PhaseCode.s2N + ]: + measPhase = findPrimaryPhase(meas.PowerSystemResource) + if measPhase.value not in voltage_violation_phases: + voltage_violation_phases += measPhase.value + else: + total_violation_time = int(timestamp) - self.voltage_violation_time + break + if total_violation_time > ConservationVoltageReductionController.max_violation_time: + if voltage_violation_phases is None: + logger.info(f"voltage_violation_phases is None: phase of voltage violation could not be found.") + self.cvr_control(voltage_violation_phases) + self.voltage_violation_time = int(timestamp) + if isinstance(sim, Simulation): + self.simulation.resume() + + def findFeederHeadLoadMeasurements(self): + energySources = self.graph_model.graph.get(cim.EnergySource, {}) + feederLoadObject = None + for eSource in energySources.values(): + feederLoadObject = findFeederHeadPowerMeasurmentsObject(eSource) + feederLoadMeasurements = feederLoadObject.Measurements + for measurement in feederLoadMeasurements: + if measurement.measurementType == 'VA': + if measurement.phases in [cim.PhaseCode.A, cim.PhaseCode.AN]: + self.peak_va_measurements_A[measurement.mRID] = {'object': measurement, 'value': None} + elif measurement.phases in [cim.PhaseCode.B, cim.PhaseCode.BN]: + self.peak_va_measurements_B[measurement.mRID] = {'object': measurement, 'value': None} + elif measurement.phases in [cim.PhaseCode.C, cim.PhaseCode.CN]: + self.peak_va_measurements_C[measurement.mRID] = {'object': measurement, 'value': None} + if not self.peak_va_measurements_A or not self.peak_va_measurements_B or not self.peak_va_measurements_C: + raise RuntimeError(f'feeder {self.graph_model.container.mRID}, has no measurements associated with the ' + 'feeder head transformer!') + + def calculate_per_unit_voltage(self): + for mrid in self.pnv_measurements.keys(): + self.pnv_measurements_pu[mrid] = None + meas_base = None + meas = self.pnv_measurements.get(mrid) + if (meas is None) or (meas.get('measurement_value') is None): + self.log.warning(f'Measurement {mrid} is missing. It will be ignored in cvr control.') + continue + meas_value = meas.get('measurement_value').get('magnitude') + if meas_value is None: + self.log.error(f'An unexpected value of None was recieved for PNV measurement {mrid}. It will be ' + 'ignored in cvr control.') + continue + meas_obj = meas.get('measurement_object') + if (meas_obj is None) or (not isinstance(meas_obj, cim.Measurement)): + self.log.error(f'The cim.Measurement object {mrid} associated with this measurement value is missing! ' + 'It will be ignored in cvr control.') + logger.error(f'The measurement dictionary for mrid {mrid} is missing from the CIM database.') + continue + + meas_term = meas_obj.Terminal + if isinstance(meas_obj.PowerSystemResource, cim.PowerElectronicsConnection): + meas_base = float(meas_obj.PowerSystemResource.ratedU) + elif isinstance(meas_term, cim.Terminal): + if meas_term.TransformerEnd: + if isinstance(meas_term.TransformerEnd[0], cim.PowerTransformerEnd): + meas_base = float(meas_term.TransformerEnd[0].ratedU) + elif isinstance(meas_term.TransformerEnd[0], cim.TransformerTankEnd): + if isinstance(meas_term.TranformerEnd[0].BaseVoltage, cim.BaseVoltage): + meas_base = float(meas_term.TransformerEnd[0].BaseVoltage.nominalVoltage) + else: + self.log.error('meas_term.TransformerEnd[0].BaseVoltage is None') + logger.error('meas_term.TransformerEnd[0].BaseVoltage is None') + raise RuntimeError('PNV measurement BaseVoltage for ConnectivityNode ' + f'{meas_term.ConnectivityNode.name} is None') + elif isinstance(meas_term.ConductingEquipment, cim.ConductingEquipment): + if isinstance(meas_term.ConductingEquipment.BaseVoltage, cim.BaseVoltage): + meas_base = float(meas_term.ConductingEquipment.BaseVoltage.nominalVoltage) + else: + logger.error(f'meas_term.ConductingEquipment.BaseVoltage is None') + else: + logger.error(f'meas_term.ConductingEquipment is None') + if (meas_base is None) or (meas_base < 1e-10): + self.log.error(f'Unable to get the nominal voltage for measurement with mrid {mrid}.') + logger.error('Voltage Measurement has no associated nominal Voltage.\nMeasurement: ' + f'{meas_obj.name}\nTerminal: {meas_obj.Terminal.name}\n') + continue + pnv_meta_data = {'measured PNV': meas_value, 'measured LLV': meas_value * math.sqrt(3.0), 'rated V': meas_base} + pu_vals = [meas_value / meas_base, meas_value * math.sqrt(3.0) / meas_base] + diff = max(pu_vals) + for i in range(len(pu_vals)): + if abs(1.0-pu_vals[i]) < diff: + diff = abs(1.0-pu_vals[i]) + self.pnv_measurements_pu[mrid] = pu_vals[i] + pnv_meta_data['pu V'] = self.pnv_measurements_pu[mrid] + logger.info(f'Voltage meta data for {meas_obj.Terminal.ConnectivityNode.name}:{json.dumps(pnv_meta_data,indent=4)}') + + def save_total_load_data(self, time, total_loads): + file_path = self.model_results_path/"total_load.csv" + header = ["time", "phase", "p", "q"] + for phase, load in total_loads.items(): + data = [time, phase, load.real, load.imag] + add_data_to_csv(file_path, data, header=header) + + def save_voltage_data(self, time): + file_path = self.model_results_path/"voltages.csv" + header = ["time", "mrid", "node", "phase", "voltage"] + for mrid in self.pnv_measurements.keys(): + v = self.pnv_measurements_pu[mrid] + phase = self.pnv_measurements[mrid]["measurement_object"].phases.value + node = self.pnv_measurements[mrid]["measurement_object"].Terminal.ConnectivityNode.name + if v is not None: + data = [time, mrid, node, phase, v] + add_data_to_csv(file_path, data, header=header) + + def on_measurement_callback(self, header: Dict[str, Any], message: Dict[str, Any]): + timestamp = message.get('message', {}).get('timestamp', '') + measurements = message.get('message', {}).get('measurements', {}) + self.on_measurement(None, timestamp, measurements) + + def create_opendss_context(self): + message = { + 'configurationType': 'DSS Base', + 'parameters': { + 'i_fraction': '1.0', + 'z_fraction': '0.0', + 'model_id': self.graph_model.container.mRID, + 'load_scaling_factor': '1.0', + 'schedule_name': 'ieeezipload', + 'p_fraction': '0.0' + }, + 'resultFormat': 'JSON' + } + base_dss_response = self.gad_obj.get_response(topics.CONFIG, message, timeout=60) + base_dss_dict = json.loads(base_dss_response.get('message', ''), strict=False) + base_dss_str = base_dss_dict.get('data', '') + fileDir = Path(__file__).parent / 'cvr_app_instances' / f'{self.id}' / 'master.dss' + fileDir.parent.mkdir(parents=True, exist_ok=True) + endOfBase = base_dss_str.find('calcv\n') + len('calcv\n') + with fileDir.open(mode='w') as f_master: + f_master.write(base_dss_str[:endOfBase]) + self.dssContext.Command(f'Redirect {fileDir}') + self.dssContext.Command(f'Compile {fileDir}') + self.dssContext.Solution.SolveNoControl() + if not self.dssContext.Solution.Converged(): + raise RuntimeError("Internal OpenDSS model fails to converge!") + + def update_opendss_with_measurements(self): + for psr_mrid in self.va_measurements.keys(): + va_val = complex(0.0, 0.0) + val_is_valid = True + for meas_mrid in self.va_measurements[psr_mrid]['measurement_values'].keys(): + meas_val = self.va_measurements[psr_mrid]['measurement_values'][meas_mrid] + if meas_val is not None: + va_val += complex( + meas_val.get('magnitude') * math.cos(math.radians(meas_val.get('angle'))), + meas_val.get('magnitude') * math.sin(math.radians(meas_val.get('angle')))) / 1000.0 + else: + val_is_valid = False + break + if val_is_valid: + meas_obj = self.va_measurements[psr_mrid]['measurement_objects'][meas_mrid].PowerSystemResource + name = meas_obj.name + if isinstance(meas_obj, cim.EnergyConsumer): + load = self.dssContext.Loads.First() + while load: + if self.dssContext.Loads.Name() == name: + break + else: + load = self.dssContext.Loads.Next() + old_val = complex(self.dssContext.Loads.kW(), self.dssContext.Loads.kvar()) + self.dssContext.Command(f'Load.{name}.kw={va_val.real}') + self.dssContext.Command(f'Load.{name}.kvar={va_val.imag}') + logger.info(f'Changing load, {name}.\nOld: {old_val}\nNew: ' + f'{complex(self.dssContext.Loads.kW(), self.dssContext.Loads.kvar())}') + elif isinstance(meas_obj, cim.PowerElectronicsConnection): + if isinstance(meas_obj.PowerElectronicsUnit[0], cim.PhotovoltaicUnit): + pv = self.dssContext.PVsystems.First() + while pv: + if self.dssContext.PVsystems.Name() == name: + break + else: + pv = self.dssContext.PVsystems.Next() + old_val = complex(self.dssContext.PVsystems.kW(), self.dssContext.PVsystems.kvar()) + self.dssContext.Command(f'PVSystem.{name}.pmpp={va_val.real}') + self.dssContext.Command(f'PVSystem.{name}.kvar={va_val.imag}') + logger.info(f'Changing PV, {name}.\nOld: {old_val}\nNew: ' + f'{complex(self.dssContext.PVsystems.kW(), self.dssContext.PVsystems.kvar())}') + elif isinstance(meas_obj.PowerElectronicsUnit[0], cim.BatteryUnit): + batt = self.dssContext.Storages.First() + while batt: + if self.dssContext.Storages.Name() == name: + break + else: + batt = self.dssContext.Storages.Next() + self.dssContext.Command(f'Storage.{name}.kw={va_val.real}') + self.dssContext.Command(f'Storage.{name}.kvar={va_val.imag}') + #TODO: update wind turbines in opendss + #TODO: update other generator/motor type object in opendss + for psr_mrid in self.pos_measurements.keys(): + num_of_meas = len(self.pos_measurements[psr_mrid]['measurement_objects'].keys()) + pos_val = [0] * num_of_meas + val_is_valid = True + key_list = list(self.pos_measurements[psr_mrid]['measurement_values'].keys()) + for meas_mrid in key_list: + meas_val = self.pos_measurements[psr_mrid]['measurement_values'][meas_mrid] + if meas_val is not None: + pos_val[key_list.index(meas_mrid)] = int(meas_val.get('value')) + else: + val_is_valid = False + break + if val_is_valid: + meas_obj = self.pos_measurements[psr_mrid]['measurement_objects'][meas_mrid].PowerSystemResource + meas_phase = self.pos_measurements[psr_mrid]['measurement_objects'][meas_mrid].phases.value + name = meas_obj.name + if isinstance(meas_obj, cim.Switch): + switch_state = sum(pos_val) + if switch_state == 0: + logger.info(f'Opening switch {name}') + self.dssContext.Command(f'open Line.{name} 1') + else: + logger.info(f'Closing switch {name}') + self.dssContext.Command(f'close Line.{name} 1') + elif isinstance(meas_obj, cim.LinearShuntCompensator): + cap = self.dssContext.Capacitors.First() + while cap: + if self.dssContext.Capacitors.Name() == name: + break + else: + cap = self.dssContext.Capacitors.Next() + if not cap: + logger.error(f'There is no capacitor with the name {name} in the OpenDSS model!') + self.log.error(f'There is no capacitor with the name {name} in the OpenDSS model!') + else: + capacitor_state = sum(pos_val) + if capacitor_state == 0: + logger.info(f'Disengaging shunt capacitor.') + self.dssContext.Capacitors.Open() + else: + logger.info(f'Engaging shunt capacitor.') + self.dssContext.Capacitors.Close() + # self.dssContext.Command(f'Capacitor.{name}.states={pos_val}') + if isinstance(meas_obj, cim.PowerTransformer): + name = self.controllable_regulators.get(psr_mrid, {}).get('PhasesToName', {}).get(meas_phase) + if name is not None: + reg = self.dssContext.Transformers.First() + while reg: + if self.dssContext.Transformers.Name() == name: + break + else: + reg = self.dssContext.Transformers.Next() + if not reg: + self.log.error(f'There is not regulator with the name {name} in the OpenDSS model!') + else: + regulationPerStep = (self.dssContext.Transformers.MaxTap() + - self.dssContext.Transformers.MinTap()) \ + / self.dssContext.Transformers.NumTaps() + tapStep = 1.0 + (pos_val[0] * regulationPerStep) + self.dssContext.Command(f'Transformer.{name}.Taps=[1.0, {tapStep}]') + self.dssContext.Solution.SolveNoControl() + if not self.dssContext.Solution.Converged(): + raise RuntimeError("Internal OpenDSS model fails to converge after updating with measurements!") + + def cvr_control(self, phases: str = None): + capacitor_list = [] + for cap_mrid, cap in self.controllable_capacitors.items(): + cap_meas_dict = self.pos_measurements.get(cap_mrid) + pos = None + if cap_meas_dict is not None: + cap_meas_val_dict = cap_meas_dict.get('measurement_values') + pos_list = list(cap_meas_val_dict.values()) + if pos_list: + pos = pos_list[0] + if pos is not None: + # For a capacitor, pos = 1 means the capacitor is connected to the grid. Otherwise, it's 0. + capacitor_list.append((cap_mrid, cap, cap.bPerSection, pos['value'], True)) + regulator_list = [] + for psr_mrid, psr_dict in self.controllable_regulators.items(): + tapChangersList = self.controllable_regulators.get(psr_mrid, {}).get('RatioTapChangers', []) + if tapChangersList: + regulator_list.append((psr_mrid, tapChangersList, phases)) + if not (capacitor_list or regulator_list): + return + meas_list = [] + # FIXMEOr: Is this going to pnv_measurements_pu too many times? + for meas_mrid in self.pnv_measurements_pu.keys(): + if self.pnv_measurements_pu[meas_mrid] is not None: + meas_list.append(self.pnv_measurements_pu[meas_mrid]) + if meas_list: + if min(meas_list) > self.low_volt_lim: + if capacitor_list: + sorted(capacitor_list, key=lambda x: x[2], reverse=True) + local_capacitor_list = [] + for element_tuple in capacitor_list: + if element_tuple[3] == 1: + local_capacitor_list.append(element_tuple) + self.desired_setpoints.update(self.decrease_voltage_capacitor(local_capacitor_list)) + if regulator_list: + self.desired_setpoints.update(self.decrease_voltage_regulator(regulator_list)) + else: + if capacitor_list: + sorted(capacitor_list, key=lambda x: x[2]) + local_capacitor_list = [] + for element_tuple in capacitor_list: + if element_tuple[3] == 0: + local_capacitor_list.append(element_tuple) + self.desired_setpoints.update(self.increase_voltage_capacitor(local_capacitor_list)) + if regulator_list: + self.desired_setpoints.update(self.increase_voltage_regulator(regulator_list)) + + # FIXMEOr: Check the update_opendss_with_measurements function. + + if self.desired_setpoints: + self.send_setpoints() + + return + + def increase_voltage_capacitor(self, cap_list: list) -> dict: + return_dict = {} + success = False + while cap_list and not success: + element_tuple = cap_list.pop(0) + if not element_tuple[4]: + continue + cap_mrid = element_tuple[0] + cap_obj = element_tuple[1] + pos_val = [1] + #saved_states = self.dssContext.Command(f'Capacitor.{cap_obj.name}.states') + cap = self.dssContext.Capacitors.First() + while cap: + if cap_obj.name == self.dssContext.Capacitors.Name(): + break + else: + cap = self.dssContext.Capacitors.Next() + if cap: + saved_states = self.dssContext.Capacitors.States() + else: + saved_states = [1] + if saved_states[0] != pos_val[0]: + self.dssContext.Capacitors.Close() + self.dssContext.Solution.SolveNoControl() + converged = self.dssContext.Solution.Converged() + if converged: + logger.info(f'Powerflow converged when closing cap {cap_mrid}.') + return_dict[cap_mrid] = {'setpoint': 1, 'object': cap_obj} + if min(self.dssContext.Circuit.AllBusMagPu()) > self.low_volt_lim: + logger.info(f'Voltage violations were eliminated when closing cap {cap_mrid}.') + success = True + else: + self.dssContext.Command(f'Capacitor.{cap_obj.name}.states={saved_states}') + return return_dict + + def decrease_voltage_capacitor(self, cap_list: list) -> dict: + return_dict = {} + while cap_list: + element_tuple = cap_list.pop(0) + if not element_tuple[4]: + continue + cap_mrid = element_tuple[0] + cap_obj = element_tuple[1] + pos_val = [0] + # saved_states = self.dssContext.Command(f'Capacitor.{cap_obj.name}.states') + cap = self.dssContext.Capacitors.First() + while cap: + if cap_obj.name == self.dssContext.Capacitors.Name(): + break + else: + cap = self.dssContext.Capacitors.Next() + if cap: + saved_states = self.dssContext.Capacitors.States() + logger.info(f'saved states:{saved_states}') + else: + saved_states = [0] + if saved_states[0] != pos_val[0]: + self.dssContext.Command(f'Capacitor.{cap_obj.name}.states={pos_val}') + self.dssContext.Solution.SolveNoControl() + converged = self.dssContext.Solution.Converged() + success = False + if converged: + logger.info(f'Powerflow converged when opening cap {cap_mrid}.') + if min(self.dssContext.Circuit.AllBusMagPu()) > self.low_volt_lim: + logger.info(f'Opening cap {cap_mrid} caused no voltage violations.') + return_dict[cap_mrid] = {'setpoint': 0, 'object': cap_obj} + success = True + if not success: + logger.info(f'Opening cap {cap_mrid} caused voltage violations or powerflow did not converge.') + self.dssContext.Command(f'Capacitor.{cap_obj.name}.states={saved_states}') + logger.info(f'return_dict length: {len(return_dict)}') + return return_dict + + def decrease_voltage_regulator(self, reg_list: list) -> dict: + tap_budget = 5 + return_dict = {} + while reg_list: + element_tuple = reg_list.pop(0) + psr_mrid = element_tuple[0] + oldSetpointRatio = [None] * len(element_tuple[1]) + newSetpointRatio = [None] * len(element_tuple[1]) + stopDecrease = [False] * len(element_tuple[1]) + for i in range(tap_budget): + if False not in stopDecrease: + break + for j in range(len(element_tuple[1])): + if stopDecrease[j]: + continue + rtp = element_tuple[1][j] + xfmrEnd = rtp.TransformerEnd + rtpName = '' + if isinstance(xfmrEnd, cim.PowerTransformerEnd): + rtpName = xfmrEnd.PowerTransformer.name + elif isinstance(xfmrEnd, cim.TransformerTankEnd): + rtpName = xfmrEnd.TransformerTank.name + reg = self.dssContext.Transformers.First() + while reg: + if self.dssContext.Transformers.Name() == rtpName: + break + else: + reg = self.dssContext.Transformers.Next() + if reg: + self.dssContext.Transformers.Wdg(2) + if i == 0: + oldSetpointRatio[j] = self.dssContext.Transformers.Tap() + newSetpointRatio[j] = oldSetpointRatio[j] + minTapRatio = self.dssContext.Transformers.MinTap() + rtpCurrentTap = newSetpointRatio[j] + tapStep = float(rtp.stepVoltageIncrement) / 100.0 + if rtpCurrentTap <= minTapRatio: + break + self.dssContext.Transformers.Tap(rtpCurrentTap - tapStep) + self.dssContext.Solution.SolveNoControl() + converged = self.dssContext.Solution.Converged() + if converged: + logger.info(f'Powerflow converged when modifying regulator {psr_mrid}.') + if min(self.dssContext.Circuit.AllBusMagPu()) > self.low_volt_lim: + logger.info(f'modifying regulator {psr_mrid} did not cause a voltage violation.') + if rtp.mRID not in return_dict.keys(): + return_dict[rtp.mRID] = { + 'setpoint': tapRatioToTapPosition(rtpCurrentTap + tapStep, rtp), + 'old_setpoint': tapRatioToTapPosition(oldSetpointRatio[j], rtp), + 'object': rtp + } + else: + return_dict[rtp.mRID]['setpoint'] = tapRatioToTapPosition( + rtpCurrentTap - tapStep, rtp) + rtpCurrentTap = rtpCurrentTap - tapStep + newSetpointRatio[j] = rtpCurrentTap + else: + stopDecrease[j] = True + self.dssContext.Transformers.Tap(newSetpointRatio[j]) + else: + self.dssContext.Transformers.Tap(newSetpointRatio[j]) + break + return return_dict + + def increase_voltage_regulator(self, reg_list: list) -> dict: + tap_budget = 5 + return_dict = {} + success = False + while reg_list and not success: + element_tuple = reg_list.pop(0) + psr_mrid = element_tuple[0] + phases = element_tuple[2] + oldSetpointRatio = [None] * len(element_tuple[1]) + newSetpointRatio = [None] * len(element_tuple[1]) + for i in range(tap_budget): + if success: + break + for j in range(len(element_tuple[1])): + rtp = element_tuple[1][j] + rtpPhases = getRatioTapChangerPhases(rtp) + if phases is None: + phases = cim.PhaseCode.ABC.value + if not (rtpPhases in phases or phases in rtpPhases): + continue + if success: + break + xfmrEnd = rtp.TransformerEnd + rtpName = '' + if isinstance(xfmrEnd, cim.PowerTransformerEnd): + rtpName = xfmrEnd.PowerTransformer.name + elif isinstance(xfmrEnd, cim.TransformerTankEnd): + rtpName = xfmrEnd.TransformerTank.name + reg = self.dssContext.Transformers.First() + while reg: + if self.dssContext.Transformers.Name() == rtpName: + break + else: + reg = self.dssContext.Transformers.Next() + if reg: + self.dssContext.Transformers.Wdg(2) + if i == 0: + oldSetpointRatio[j] = self.dssContext.Transformers.Tap() + newSetpointRatio[j] = oldSetpointRatio[j] + maxTapRatio = self.dssContext.Transformers.MaxTap() + rtpCurrentTap = newSetpointRatio[j] + tapStep = float(rtp.stepVoltageIncrement) / 100.0 + if rtpCurrentTap >= maxTapRatio: + break + self.dssContext.Transformers.Tap(rtpCurrentTap + tapStep) + self.dssContext.Solution.SolveNoControl() + converged = self.dssContext.Solution.Converged() + if converged: + logger.info(f'Powerflow converged when modifying regulator {psr_mrid}.') + if rtp.mRID not in return_dict.keys(): + return_dict[rtp.mRID] = { + 'setpoint': tapRatioToTapPosition(rtpCurrentTap + tapStep, rtp), + 'old_setpoint': tapRatioToTapPosition(oldSetpointRatio[j], rtp), + 'object': rtp + } + else: + return_dict[rtp.mRID]['setpoint'] = tapRatioToTapPosition(rtpCurrentTap + tapStep, rtp) + rtpCurrentTap = rtpCurrentTap + tapStep + newSetpointRatio[j] = rtpCurrentTap + if min(self.dssContext.Circuit.AllBusMagPu()) > self.low_volt_lim: + logger.info(f'Voltage violations were eliminated when modifying regulator {psr_mrid}.') + success = True + break + else: + self.dssContext.Transformers.Tap(rtpCurrentTap) + break + return return_dict + + def reg_openDSS_to_CIMHub(self, turnsRatio): + return turnsRatio + + def send_setpoints(self): + self.differenceBuilder.clear() + for mrid, dictVal in self.desired_setpoints.items(): + cimObj = dictVal.get('object') + newSetpoint = dictVal.get('setpoint') + oldSetpoint = dictVal.get('old_setpoint') + if isinstance(cimObj, cim.LinearShuntCompensator): + phases = "" + for shuntCompensatorPhase in cimObj.ShuntCompensatorPhase: + phases += shuntCompensatorPhase.phase.value + if not phases: + phases = "ABC" + pos_setpoint_command_debug_data = { + "object_type": f"{type(cimObj).__name__}", + "name": cimObj.name, + "phase": phases, + "value": newSetpoint + } + logger.info(f"Pos Setpoint Change: {json.dumps(pos_setpoint_command_debug_data)}") + self.differenceBuilder.add_difference(mrid, 'ShuntCompensator.sections', newSetpoint, + int(not newSetpoint)) + data = [self.timestamp, mrid, "", "", newSetpoint] + add_data_to_csv(self.model_results_path/"caps.csv", data=data) + elif isinstance(cimObj, cim.RatioTapChanger): + phases = "" + if isinstance(cimObj.TransformerEnd, cim.PowerTransformerEnd): + phases = "ABC" + else: + phases = cimObj.TransformerEnd.phases.value + pos_setpoint_command_debug_data = { + "object_type": f"{type(cimObj).__name__}", + "name": cimObj.name, + "phase": phases, + "value": newSetpoint + } + logger.info(f"Pos Setpoint Change: {json.dumps(pos_setpoint_command_debug_data)}") + self.differenceBuilder.add_difference(cimObj.mRID, 'TapChanger.step', newSetpoint, oldSetpoint) + data = [self.timestamp, mrid, "", "", newSetpoint] + add_data_to_csv(self.model_results_path/"regs.csv", data=data) + else: + self.log.warning(f'The CIM object with mRID, {mrid}, is not a cim.LinearShuntCompensator or a ' + f'cim.PowerTransformer. The object is a {type(cimObj)}. This application will ignore ' + 'sending a setpoint to this object.') + setpointMessage = json.dumps(self.differenceBuilder.get_message(), indent=4, sort_keys=True) + logger.info(f"send_setpoint():\ntopic: {self.setpoints_topic}\nmessage: {type(setpointMessage)}") + self.gad_obj.send(self.setpoints_topic, setpointMessage) + + def simulation_completed(self, sim: Simulation): + self.log.info(f'Simulation for ConservationVoltageReductionController:{self.id} has finished. This application ' + 'instance can be deleted.') + self.isValid = False + + def __del__(self): + directoryToDelete = Path(__file__).parent / 'cvr_app_instances' / f'{self.id}' + removeDirectory(directoryToDelete) + + +def findFeederHeadPowerMeasurmentsObject(cimObj: object): + ''' + Helper function to find feeder transformer + ''' + if not isinstance(cimObj, cim.EnergySource): + raise TypeError('findFeederHeadPowerMeasurmentsObject(): cimObj must be an instance of cim.EnergySource!') + equipmentToCheck = [cimObj] + feederHeadLoadMeasurementsObject = None + i = 0 + while not feederHeadLoadMeasurementsObject and i < len(equipmentToCheck): + equipmentToAdd = [] + for eq in equipmentToCheck[i:]: + powerMeasurementsCount = 0 + for meas in eq.Measurements: + if meas.measurementType == 'VA': + powerMeasurementsCount += 1 + if powerMeasurementsCount == 3: + feederHeadLoadMeasurementsObject = eq + break + else: + terminals = eq.Terminals + connectivityNodes = [] + for t in terminals: + if t.ConnectivityNode not in connectivityNodes: + connectivityNodes.append(t.ConnectivityNode) + for cn in connectivityNodes: + for t in cn.Terminals: + if t.ConductingEquipment not in equipmentToCheck and t.ConductingEquipment not in equipmentToAdd: + equipmentToAdd.append(t.ConductingEquipment) + i = len(equipmentToCheck) + equipmentToCheck.extend(equipmentToAdd) + if not feederHeadLoadMeasurementsObject: + raise RuntimeError('findFeederPowerRating(): No feeder head object with power measurements could be found for ' + f'EnergySource, {cimObj.name}!') + return feederHeadLoadMeasurementsObject + + +def findPrimaryPhase(cimObj): + ''' + Helper function for finding the primary phase an instance of cim.ConductingEquipment on the secondary + system is connected to. + ''' + if not isinstance(cimObj, cim.ConductingEquipment): + raise TypeError('findPrimaryPhase(): cimObj must be an instance of cim.ConductingEquipment!') + equipmentToCheck = [cimObj] + phaseCode = None + xfmr = None + i = 0 + while not xfmr and i < len(equipmentToCheck): + equipmentToAdd = [] + for eq in equipmentToCheck[i:]: + if isinstance(eq, cim.PowerTransformer): + xfmr = eq + break + else: + terminals = eq.Terminals + connectivityNodes = [] + for t in terminals: + if t.ConnectivityNode not in connectivityNodes: + connectivityNodes.append(t.ConnectivityNode) + for cn in connectivityNodes: + for t in cn.Terminals: + if t.ConductingEquipment not in equipmentToCheck and t.ConductingEquipment not in equipmentToAdd: + equipmentToAdd.append(t.ConductingEquipment) + i = len(equipmentToCheck) + equipmentToCheck.extend(equipmentToAdd) + if not xfmr: + raise RuntimeError('findPrimaryPhase(): no upstream centertapped transformer could be found for secondary ' + f'system object {cimObj.name}!') + for tank in xfmr.TransformerTanks: + if phaseCode is not None: + break + for tankEnd in tank.TransformerTankEnds: + if tankEnd.phases not in [ + cim.PhaseCode.none, cim.PhaseCode.s1, cim.PhaseCode.s12, cim.PhaseCode.s12N, cim.PhaseCode.s1N, + cim.PhaseCode.s2, cim.PhaseCode.s2N + ]: + phaseCode = tankEnd.phases + break + if not phaseCode: + raise RuntimeError('findPrimaryPhase(): the upstream centertapped transformer has no primary phase defined!?') + return phaseCode + + +def getRatioTapChangerPhases(ratioTapChanger) -> str: + if not isinstance(ratioTapChanger, cim.RatioTapChanger): + raise TypeError('getRatioTapChanger(): Argument ratioTapChanger must be an instance of cim.RatioTapChanger!') + phases = None + transformerEnd = ratioTapChanger.TransformerEnd + if isinstance(transformerEnd, cim.PowerTransformerEnd): + phases = cim.PhaseCode.ABC.value + elif isinstance(transformerEnd, cim.TransformerTankEnd): + phases = transformerEnd.phases.value + if phases is None: + raise ValueError('getRatioTapChanger(): The TransformerEnd associated with the RatioTapChanger is expected to ' + 'be an instance of PowerTransformerEnd or TransformerTankEnd!') + return phases + + +def tapPositionToTapRatio(tapPosition: int, ratioTapChanger) -> float: + '''Converts integer tap position to a pu turns ratio for a given cim.RatioTapChanger instance''' + if not isinstance(tapPosition, int): + raise TypeError('tapPositionToTapRatio(): tapPosition must be an int.') + if not isinstance(ratioTapChanger, cim.RatioTapChanger): + raise TypeError('tapPositionToTapRatio(): ratioTapChanger must be a cim.RatioTapChanger.') + maxTapPosition = int(ratioTapChanger.highStep) + minTapPosition = int(ratioTapChanger.lowStep) + neutralTapPosition = int(ratioTapChanger.neutralStep) + stepRatio = float(ratioTapChanger.stepVoltageIncrement) / 100.0 + if tapPosition > maxTapPosition or tapPosition < minTapPosition: + raise ValueError('tapPostionToTapRatio(): tapPosition must be within the maximum and minimum tap positions of ' + 'the ratioTapChanger!') + if neutralTapPosition is None: + neutralTapPosition = 0 + tapPostionRatio = 1.0 + (float(tapPosition - neutralTapPosition) * stepRatio) + return tapPostionRatio + + +def tapRatioToTapPosition(tapPositionRatio: float, ratioTapChanger) -> float: + '''Converts pu turn ratio tap setting to an integer tap position for a given cim.RatioTapChanger instance''' + if not isinstance(tapPositionRatio, float): + raise TypeError('tapRatioToTapPosition(): tapPositionRatio must be a float.') + if not isinstance(ratioTapChanger, cim.RatioTapChanger): + raise TypeError('tapRatioToTapPosition(): ratioTapChanger must be a cim.RatioTapChanger.') + maxTapPosition = int(ratioTapChanger.highStep) + minTapPosition = int(ratioTapChanger.lowStep) + neutralTapPosition = int(ratioTapChanger.neutralStep) + stepRatio = float(ratioTapChanger.stepVoltageIncrement) / 100.0 + if tapPositionRatio > tapPositionToTapRatio(maxTapPosition, + ratioTapChanger) or tapPositionRatio < tapPositionToTapRatio( + minTapPosition, ratioTapChanger): + raise ValueError('tapRatioToTapPosition(): tapPositionRatio must be within the maximum and minimum tap ' + 'positions of the ratioTapChanger!') + if neutralTapPosition is None: + neutralTapPosition = 0 + if tapPositionRatio < 1.0: + tapPosition = neutralTapPosition - math.floor(((1.0 - tapPositionRatio) / stepRatio) + 0.5) + else: + tapPosition = neutralTapPosition + math.floor(((tapPositionRatio - 1.0) / stepRatio) + 0.5) + if tapPosition > maxTapPosition: + tapPosition = maxTapPosition + elif tapPosition < minTapPosition: + tapPosition = minTapPosition + return tapPosition + + +def buildGraphModel(mrid: str) -> FeederModel: + global CIM_GRAPH_MODELS + if not isinstance(mrid, str): + raise TypeError(f'The mrid passed to the convervation voltage reduction application must be a string.') + if mrid not in CIM_GRAPH_MODELS.keys(): + feeder_container = cim.Feeder(mRID=mrid) + graph_model = FeederModel(connection=DB_CONNECTION, container=feeder_container, distributed=False) + # graph_model.get_all_edges(cim.PowerTransformer) + # graph_model.get_all_edges(cim.TransformerTank) + # graph_model.get_all_edges(cim.Asset) + # graph_model.get_all_edges(cim.LinearShuntCompensator) + # graph_model.get_all_edges(cim.PowerTransformerEnd) + # graph_model.get_all_edges(cim.TransformerEnd) + # graph_model.get_all_edges(cim.TransformerMeshImpedance) + # graph_model.get_all_edges(cim.TransformerTankEnd) + # graph_model.get_all_edges(cim.TransformerTankInfo) + # graph_model.get_all_edges(cim.RatioTapChanger) + # graph_model.get_all_edges(cim.TapChanger) + # graph_model.get_all_edges(cim.TapChangerControl) + # graph_model.get_all_edges(cim.TapChangerInfo) + # graph_model.get_all_edges(cim.LinearShuntCompensatorPhase) + # graph_model.get_all_edges(cim.Terminal) + # graph_model.get_all_edges(cim.ConnectivityNode) + # graph_model.get_all_edges(cim.BaseVoltage) + # graph_model.get_all_edges(cim.EnergySource) + # graph_model.get_all_edges(cim.EnergyConsumer) + # graph_model.get_all_edges(cim.ConformLoad) + # graph_model.get_all_edges(cim.NonConformLoad) + # graph_model.get_all_edges(cim.EnergyConsumerPhase) + # graph_model.get_all_edges(cim.LoadResponseCharacteristic) + # graph_model.get_all_edges(cim.PowerCutZone) + # graph_model.get_all_edges(cim.Analog) + # graph_model.get_all_edges(cim.Discrete) + cimUtils.get_all_data(graph_model) + CIM_GRAPH_MODELS[mrid] = graph_model + return CIM_GRAPH_MODELS[mrid] + + +def createSimulation(gad_obj: GridAPPSD, model_info: Dict[str, Any]) -> Simulation: + if not isinstance(gad_obj, GridAPPSD): + raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + if not isinstance(model_info, dict): + raise TypeError(f'The model_id must be a dictionary type.') + line_name = model_info.get('modelId') + subregion_name = model_info.get('subRegionId') + region_name = model_info.get('regionId') + sim_name = model_info.get('modelName') + if line_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelId.') + if subregion_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key subRegionId.') + if region_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key regionId.') + if sim_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelName.') + psc = PowerSystemConfig(Line_name=line_name, + SubGeographicalRegion_name=subregion_name, + GeographicalRegion_name=region_name) + # start_time = int(datetime.utcnow().replace(microsecond=0).timestamp()) + start = datetime(2023, 1, 1, 0, tzinfo=timezone.utc) + epoch = datetime.timestamp(start) + duration = timedelta(hours=24).total_seconds() + sim_args = SimulationArgs( + start_time=epoch, + duration=duration, + simulation_name=sim_name, + run_realtime=False, + pause_after_measurements=False) + sim_config = SimulationConfig(power_system_config=psc, simulation_config=sim_args) + sim_obj = Simulation(gapps=gad_obj, run_config=sim_config) + return sim_obj + + +def createGadObject() -> GridAPPSD: + gad_user = os.environ.get('GRIDAPPSD_USER') + if gad_user is None: + os.environ['GRIDAPPSD_USER'] = 'system' + gad_password = os.environ.get('GRIDAPPSD_PASSWORD') + if gad_password is None: + os.environ['GRIDAPPSD_PASSWORD'] = 'manager' + gad_app_id = os.environ.get('GRIDAPPSD_APPLICATION_ID') + if gad_app_id is None: + os.environ['GRIDAPPSD_APPLICATION_ID'] = 'ConservationVoltageReductionApplication' + return GridAPPSD() + + +def removeDirectory(directory: Path | str): + if isinstance(directory, str): + directory = Path(directory) + for item in directory.iterdir(): + if item.is_dir(): + removeDirectory(item) + else: + item.unlink() + directory.rmdir() + + +def add_data_to_csv(file_path, data, header=None): + # Check if the file exists + file_exists = os.path.isfile(file_path) + # Open the CSV file in append mode (or create it if it doesn't exist) + with open(file_path, mode='a', newline='') as file: + writer = csv.writer(file) + # If the file doesn't exist, write the header first + if not file_exists and header: + writer.writerow(header) + # Write the data (a list or tuple representing a row) + writer.writerow(data) + + +def main(control_enabled: bool, start_simulations: bool, model_id: str = None): + if not isinstance(control_enabled, bool): + raise TypeError(f'Argument, control_enabled, must be a boolean.') + if not isinstance(start_simulations, bool): + raise TypeError(f'Argument, start_simulations, must be a boolean.') + if not isinstance(model_id, str) and model_id is not None: + raise TypeError( + f'The model id passed to the convervation voltage reduction application must be a string type or {None}.') + global cim, DB_CONNECTION + cim_profile = 'rc4_2021' + iec61970_301 = 7 + cim = importlib.import_module(f'cimgraph.data_profile.{cim_profile}') + # params = ConnectionParameters(cim_profile=cim_profile, iec61970_301=iec61970_301) + # DB_CONNECTION = GridappsdConnection(params) + params = ConnectionParameters(url='http://localhost:8889/bigdata/namespace/kb/sparql', + cim_profile=cim_profile, + iec61970_301=iec61970_301) + DB_CONNECTION = BlazegraphConnection(params) + gad_object = createGadObject() + gad_log = gad_object.get_logger() + platform_simulations = {} + local_simulations = {} + app_instances = {'field_instances': {}, 'external_simulation_instances': {}, 'local_simulation_instances': {}} + response = gad_object.query_model_info() + models = response.get('data', {}).get('models', []) + model_is_valid = False + if model_id is not None: + for m in models: + m_id = m.get('modelId') + if model_id == m_id: + model_is_valid = True + break + if not model_is_valid: + raise ValueError(f'The model id provided does not exist in the GridAPPS-D plaform.') + else: + models = [m] + if start_simulations: + for m in models: + local_simulations[m.get('modelId', '')] = createSimulation(gad_object, m) + else: + #TODO: query platform for running simulations which is currently not implemented in the GridAPPS-D Api + pass + # Create an cvr controller instance for all the real systems in the database + # for m in models: + # m_id = m.get('modelId') + # app_instances['field_instances'][m_id] = ConservationVoltageReductionController(gad_object, m_id) + for sim_id, m_id in platform_simulations.items(): + app_instances['external_simulation_instances'][sim_id] = ConservationVoltageReductionController(gad_object, + m_id, + sim_id=sim_id) + for m_id, simulation in local_simulations.items(): + app_instances['local_simulation_instances'][m_id] = ConservationVoltageReductionController( + gad_object, m_id, simulation=simulation) + app_instances_exist = False + if len(app_instances['field_instances']) > 0: + app_instances_exist = True + elif len(app_instances['external_simulation_instances']) > 0: + app_instances_exist = True + elif len(app_instances['local_simulation_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if app_instances_exist: + gad_log.info('ConservationVoltageReductionApplication successfully started.') + gad_object.set_application_status(ProcessStatusEnum.RUNNING) + + while gad_object.connected: + try: + app_instances_exist = False + #TODO: check platform for any running external simulations to control. + invalidInstances = [] + for m_id, app in app_instances.get('field_instances', {}).items(): + if not app.isValid: + invalidInstances.append(m_id) + for m_id in invalidInstances: + del app_instances['field_instances'][m_id] + invalidInstances = [] + for sim_id, app in app_instances.get('external_simulation_instances', {}).items(): + #TODO: check if external simulation has finished and shutdown the corresponding app instance. + if not app.isValid: + invalidInstances.append(sim_id) + for sim_id in invalidInstances: + del app_instances['external_simulation_instances'][sim_id] + invalidInstances = [] + for m_id, app in app_instances.get('local_simulation_instances', {}).items(): + if not app.isValid: + invalidInstances.append(m_id) + for m_id in invalidInstances: + del app_instances['local_simulation_instances'][m_id] + if len(app_instances['field_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if len(app_instances['external_simulation_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if len(app_instances['local_simulation_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if not app_instances_exist: + application_downtime = int(time.time()) - application_uptime + if application_downtime > 3600: + gad_log.info('There have been no running instances ConservationVoltageReductionController ' + 'instances for an hour. Shutting down the application.') + gad_object.set_application_status(ProcessStatusEnum.STOPPING) + gad_object.disconnect() + except KeyboardInterrupt: + gad_log.info('Manually exiting ConservationVoltageReductionApplication') + gad_object.set_application_status(ProcessStatusEnum.STOPPING) + appList = list(app_instances.get('field_instances', {})) + for app_mrid in appList: + del app_instances['field_instances'][app_mrid] + appList = list(app_instances.get('external_simulation_instances', {})) + for app_mrid in appList: + del app_instances['external_simulation_instances'][app_mrid] + appList = list(app_instances.get('local_simulation_instances', {})) + for app_mrid in appList: + del app_instances['local_simulation_instances'][app_mrid] + gad_object.disconnect() + + +if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('model_id', nargs='?', default=None, help='The mrid of the cim model to perform cvr on.') + parser.add_argument('-s', + '--start_simulations', + action='store_true', + help='Flag to have application start simulations') + parser.add_argument('-d', + '--disable_control', + action='store_true', + help='Flag to disable control on startup by default.') + args = parser.parse_args() + main(args.disable_control, args.start_simulations, args.model_id) diff --git a/micro-apps/conservation-voltage-reduction-app/pyproject.toml b/micro-apps/conservation-voltage-reduction-app/pyproject.toml new file mode 100644 index 0000000..ef760bc --- /dev/null +++ b/micro-apps/conservation-voltage-reduction-app/pyproject.toml @@ -0,0 +1,40 @@ +[tool.poetry] +name = "conservation-voltage-reduction-app" +version = "0.1.0" +description = "" +authors = ["afisher1 <4552674+afisher1@users.noreply.github.com>"] +readme = "README.md" + +[tool.poetry.dependencies] +python = "^3.10" +gridappsd-python = {version = "^2024", allow-prereleases = true} +cim-graph = "^0.1.2a0" +cvxpy = "^1.5.2" +tabulate = "^0.9.0" +networkx = "^3.3" +plotly = "^5.24.1" +pandas = "^2.2.3" +matplotlib = "^3.9.2" + +[tool.poetry.group.dev.dependencies] +pre-commit = "^3.7.0" +yapf = "^0.40.2" +mypy = "^1.9.0" + +[build-system] +requires = ["poetry-core"] +build-backend = "poetry.core.masonry.api" + +[tool.yapfignore] +ignore_patterns = [ + ".venv/**", + ".pytest_cache/**", + "dist/**", + "docs/**" +] + +[tool.yapf] +based_on_style = "pep8" +spaces_before_comment = 4 +column_limit = 120 +split_before_logical_operator = true diff --git a/micro-apps/no-app/no_app.py b/micro-apps/no-app/no_app.py new file mode 100644 index 0000000..eb7f9a7 --- /dev/null +++ b/micro-apps/no-app/no_app.py @@ -0,0 +1,599 @@ +import time +import os +import math +from pprint import pprint +# import networkx as nx +# import numpy as np +import cvxpy as cp +# from pandas import DataFrame +from typing import Any, Dict +import importlib +from datetime import datetime, timedelta +from dataclasses import dataclass, asdict +from argparse import ArgumentParser +import json +from typing import Dict +from cimgraph import utils +from tabulate import tabulate +from cimgraph.databases import ConnectionParameters, BlazegraphConnection +from cimgraph.models import FeederModel +from cimgraph.data_profile.rc4_2021 import EnergyConsumer +from cimgraph.data_profile.rc4_2021 import BatteryUnit +from cimgraph.data_profile.rc4_2021 import PowerElectronicsConnection +from gridappsd import GridAPPSD, topics, DifferenceBuilder +from gridappsd.simulation import Simulation +from gridappsd.simulation import PowerSystemConfig +from gridappsd.simulation import SimulationConfig +from gridappsd.simulation import SimulationArgs +from gridappsd.simulation import ModelCreationConfig +from gridappsd import topics as t +import csv + +IEEE123_APPS = "_E3D03A27-B988-4D79-BFAB-F9D37FB289F7" +IEEE13 = "_49AD8E07-3BF9-A4E2-CB8F-C3722F837B62" +IEEE123 = "_C1C3E687-6FFD-C753-582B-632A27E28507" +IEEE123PV = "_A3BC35AA-01F6-478E-A7B1-8EA4598A685C" +SOURCE_BUS = '150r' +OUT_DIR = "outputs" +ROOT = os.getcwd() + + +class CarbonManagementApp(object): + """Carbon Management Class + This class implements a centralized algorithm for carbon management in a grid by controlling batteries in the distribution model. + + Attributes: + sim_id: The simulation id that the instance will be perfoming control on. If the sim_id is None then the + application will assume it is performing control on the actual field devices. + Type: str. + Default: None. + gad_obj: An instatiated object that is connected to the gridappsd message bus usually this should be the same object which subscribes, but that + isn't required. + Type: GridAPPSD. + Default: None. + model_id: The mrid of the cim model to perform carbon management control on. + Type: str. + Default: None. + simulation: Simulation object of the identified cim model. + Type: Simulation. + Default: None. + network: An instatiated object of knowledge graph class for distribution feeder objects. + Type: FeederModel. + Default: None. + Methods: + initialize(headers:Dict[Any], message:Dict[Any]): This callback function will trigger the application to + read the model database to find all the measurements in the model as well as batteries in the sytem. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + on_measurement(headers:Dict[Any], message:Dict[Any]): This callback function will be used to capture the + measurements dictionary needed for control. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + optimize_battery(): This is the main function for controlling the battery. + """ + + def __init__(self, + gad_obj: GridAPPSD, + model_id: str, + network: FeederModel, + sim_id: str = None, + simulation: Simulation = None): + if not isinstance(gad_obj, GridAPPSD): + raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + if not isinstance(model_id, str): + raise TypeError(f'model_id must be a str type.') + if model_id is None or model_id == '': + raise ValueError(f'model_id must be a valid uuid.') + if not isinstance(sim_id, str) and sim_id is not None: + raise TypeError(f'sim_id must be a string type or {None}!') + if not isinstance(simulation, Simulation) and simulation is not None: + raise TypeError(f'The simulation arg must be a Simulation type or {None}!') + self.gad_obj = gad_obj + self.init_batt_dis = True + self._count = 0 + self._publish_to_topic = topics.simulation_input_topic(simulation.simulation_id) + self._init_batt_diff = DifferenceBuilder(simulation.simulation_id) + + self.Battery = {} + self.Solar = {} + self.EnergyConsumer = {} + self.has_batteries = True + self.has_power_electronics = True + self.has_energy_consumers = True + if PowerElectronicsConnection not in network.graph: + self.has_power_electronics = False + self.has_batteries = False + # raise ValueError("No power electronic devices in network.") + elif len(network.graph[PowerElectronicsConnection].keys()) == 0: + self.has_power_electronics = False + self.has_batteries = False + if EnergyConsumer not in network.graph: + self.has_energy_consumers = False + + if self.has_power_electronics: + self._collect_power_electronic_devices(network) + # raise ValueError("No power electronic devices in network.") + + if len(self.Battery) == 0: + self.has_batteries = False + print("No batteries in system.") + # raise ValueError("No batteries in network.") + + if self.has_energy_consumers: + self._collect_energy_consumers(network) + + simulation.add_onmeasurement_callback(self.on_measurement) + # simulation.start_simulation() + + def _collect_power_electronic_devices(self, network): + for pec in network.graph[PowerElectronicsConnection].values(): + # inv_mrid = pec.mRID + for unit in pec.PowerElectronicsUnit: + unit_mrid = unit.mRID + print(type(unit)) + if isinstance(unit, BatteryUnit): + self.Battery[unit_mrid] = {'phases': [], 'measurementType': [], 'measurementmRID': [], + 'measurementPhases': []} + self.Battery[unit_mrid]['name'] = unit.name + self.Battery[unit_mrid]['ratedS'] = float(pec.ratedS) / 1000 + self.Battery[unit_mrid]['ratedE'] = float(unit.ratedE) / 1000 + else: + self.Solar[unit_mrid] = {'phases': [], 'measurementType': [], 'measurementmRID': [], + 'measurementPhases': []} + self.Solar[unit_mrid]['name'] = pec.name + self.Solar[unit_mrid]['ratedS'] = float(pec.ratedS) / 1000 + + if not pec.PowerElectronicsConnectionPhases: + if unit_mrid in self.Battery: + self.Battery[unit_mrid]['phases'] = 'ABC' + if unit_mrid in self.Solar: + self.Solar[unit_mrid]['phases'] = 'ABC' + else: + phases = [] + for phase in pec.PowerElectronicsConnectionPhases: + phases.append(phase.phase.value) + if unit_mrid in self.Battery: + self.Battery[unit_mrid]['phases'] = phases + if unit_mrid in self.Solar: + self.Solar[unit_mrid]['phases'] = phases + + for measurement in pec.Measurements: + if unit_mrid in self.Battery: + self.Battery[unit_mrid]['measurementType'].append(measurement.measurementType) + self.Battery[unit_mrid]['measurementmRID'].append(measurement.mRID) + if measurement.phases.value is not None: + self.Battery[unit_mrid]['measurementPhases'].append(measurement.phases.value) + if unit_mrid in self.Solar: + self.Solar[unit_mrid]['measurementType'].append(measurement.measurementType) + self.Solar[unit_mrid]['measurementmRID'].append(measurement.mRID) + if measurement.phases.value is not None: + self.Solar[unit_mrid]['measurementPhases'].append(measurement.phases.value) + + def _collect_energy_consumers(self, network): + for ld in network.graph[EnergyConsumer].values(): + ld_mrid = ld.mRID + self.EnergyConsumer[ld_mrid] = {'phases': [], 'measurementType': [], 'measurementmRID': [], + 'measurementPhases': []} + self.EnergyConsumer[ld_mrid]['name'] = ld.name + if not ld.EnergyConsumerPhase: + self.EnergyConsumer[ld_mrid]['phases'] = 'ABC' + else: + phases = [] + for phase in ld.EnergyConsumerPhase: + phases.append(phase.phase.value) + self.EnergyConsumer[ld_mrid]['phases'] = phases + for measurement in ld.Measurements: + self.EnergyConsumer[ld_mrid]['measurementType'].append(measurement.measurementType) + self.EnergyConsumer[ld_mrid]['measurementmRID'].append(measurement.mRID) + if measurement.phases.value is not None: + self.EnergyConsumer[ld_mrid]['measurementPhases'].append(measurement.phases.value) + + def find_phase(self, degree): + ref_points = [0, 120, -120] + closest = min(ref_points, key=lambda x: abs(degree - x)) + phases = {0: 'A', -120: 'B', 120: 'C'} + return phases[closest] + + def pol2cart(self, mag, angle_deg): + # Convert degrees to radians. GridAPPS-D spits angle in degrees + angle_rad = math.radians(angle_deg) + p = mag * math.cos(angle_rad) + q = mag * math.sin(angle_rad) + return p, q + + def find_injection(self, object, measurements): + + for item in object: + object[item]['P_inj'] = [0, 0, 0] + object[item]['Q_inj'] = [0, 0, 0] + meas_type = object[item]['measurementType'] + va_idx = [i for i in range(len(meas_type)) if meas_type[i] == 'VA'] + pnv_idx = [i for i in range(len(meas_type)) if meas_type[i] == 'PNV'] + soc_idx = [i for i in range(len(meas_type)) if meas_type[i] == 'SoC'] + if soc_idx: + soc = measurements[object[item]['measurementmRID'][soc_idx[0]]]['value'] + self.Battery[item]['soc'] = soc + if 's' in object[item]['phases'][0]: + angle = measurements[object[item]['measurementmRID'][pnv_idx[0]]]['angle'] + rho = measurements[object[item]['measurementmRID'][va_idx[0]]]['magnitude'] + phi = measurements[object[item]['measurementmRID'][va_idx[0]]]['angle'] + p, q = self.pol2cart(rho, phi) + if self.find_phase(angle) == 'A': + object[item]['P_inj'][0] = 2 * p / 1000 + object[item]['Q_inj'][0] = 2 * q / 1000 + object[item]['phases'] = 'A' + elif self.find_phase(angle) == 'B': + object[item]['P_inj'][1] = 2 * p / 1000 + object[item]['Q_inj'][1] = 2 * q / 1000 + object[item]['phases'] = 'B' + else: + object[item]['P_inj'][2] = 2 * p / 1000 + object[item]['Q_inj'][2] = 2 * q / 1000 + object[item]['phases'] = 'C' + elif object[item]['phases'] == 'ABC': + for k in range(3): + rho = measurements[object[item]['measurementmRID'][va_idx[k]]]['magnitude'] + phi = measurements[object[item]['measurementmRID'][va_idx[k]]]['angle'] + p, q = self.pol2cart(rho, phi) + object[item]['P_inj'][k] = p / 1000 + object[item]['Q_inj'][k] = q / 1000 + else: + rho = measurements[object[item]['measurementmRID'][va_idx[0]]]['magnitude'] + phi = measurements[object[item]['measurementmRID'][va_idx[0]]]['angle'] + p, q = self.pol2cart(rho, phi) + if object[item]['phases'][0] == 'A': + object[item]['P_inj'][0] = p / 1000 + object[item]['Q_inj'][0] = q / 1000 + elif object[item]['phases'][0] == 'B': + object[item]['P_inj'][1] = p / 1000 + object[item]['Q_inj'][1] = q / 1000 + else: + object[item]['P_inj'][2] = p / 1000 + object[item]['Q_inj'][2] = q / 1000 + + def optimize_battery(self, timestamp): + # Define optimization variables + n_batt = len(self.Battery) + p_flow_A = cp.Variable(integer=False, name='p_flow_A') + p_flow_B = cp.Variable(integer=False, name='p_flow_B') + p_flow_C = cp.Variable(integer=False, name='p_flow_C') + p_flow_mod_A = cp.Variable(integer=False, name='p_flow_mod_A') + p_flow_mod_B = cp.Variable(integer=False, name='p_flow_mod_B') + p_flow_mod_C = cp.Variable(integer=False, name='p_flow_mod_C') + p_batt = cp.Variable(n_batt, integer=False, name='p_batt') + soc = cp.Variable(n_batt, integer=False, name='soc') + lambda_c = cp.Variable(n_batt, boolean=True, name='lambda_c') + lambda_d = cp.Variable(n_batt, boolean=True, name='lambda_d') + P_batt_A = cp.Variable(n_batt, integer=False, name='P_batt_A') + P_batt_B = cp.Variable(n_batt, integer=False, name='P_batt_B') + P_batt_C = cp.Variable(n_batt, integer=False, name='P_batt_C') + deltaT = 0.25 + n_batt_ABC = 0 + for batt in self.Battery: + if 'ABC' in self.Battery[batt]['phases']: + n_batt_ABC += 1 + # Defining variable to constraint same sign for multiple 3-ph batteries + if n_batt_ABC > 0: + b = cp.Variable(n_batt_ABC, boolean=True, name='b') + + constraints = [] + sum_flow_A, sum_flow_B, sum_flow_C = 0.0, 0.0, 0.0 + for pv in self.Solar: + sum_flow_A += self.Solar[pv]['P_inj'][0] + sum_flow_B += self.Solar[pv]['P_inj'][1] + sum_flow_C += self.Solar[pv]['P_inj'][2] + for load in self.EnergyConsumer: + sum_flow_A += self.EnergyConsumer[load]['P_inj'][0] + sum_flow_B += self.EnergyConsumer[load]['P_inj'][1] + sum_flow_C += self.EnergyConsumer[load]['P_inj'][2] + + idx = 0 + idx_b = 0 + # For now, we assume battery to be 100% efficient. Need to rewrite the soc constraints if using different efficiency + for batt in self.Battery: + constraints.append( + soc[idx] == self.Battery[batt]['soc'] / 100 + p_batt[idx] * deltaT / self.Battery[batt]['ratedE']) + constraints.append(p_batt[idx] <= lambda_c[idx] * self.Battery[batt]['ratedS']) + constraints.append(p_batt[idx] >= - lambda_d[idx] * self.Battery[batt]['ratedS']) + constraints.append(lambda_c[idx] + lambda_d[idx] <= 1) + constraints.append(soc[idx] <= 0.9) + constraints.append(soc[idx] >= 0.2) + # Tracking dispatch in A, B, and C phase by different batteries + if 'ABC' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == p_batt[idx] / 3) + constraints.append(P_batt_B[idx] == p_batt[idx] / 3) + constraints.append(P_batt_C[idx] == p_batt[idx] / 3) + constraints.append(p_batt[idx] >= - b[idx_b] * self.Battery[batt]['ratedS']) + constraints.append(p_batt[idx] <= self.Battery[batt]['ratedS'] * (1 - b[idx_b])) + idx_b += 1 + else: + if 'A' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == p_batt[idx]) + constraints.append(P_batt_B[idx] == 0.0) + constraints.append(P_batt_C[idx] == 0.0) + if 'B' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == 0.0) + constraints.append(P_batt_B[idx] == p_batt[idx]) + constraints.append(P_batt_C[idx] == 0.0) + if 'C' in self.Battery[batt]['phases']: + constraints.append(P_batt_A[idx] == 0.0) + constraints.append(P_batt_B[idx] == 0.0) + constraints.append(P_batt_C[idx] == p_batt[idx]) + idx += 1 + # Ensuring three phase batteries have the same sign. We don't want one charging and another discharging + # Although, mathematically it might sound correct, it is not worth to do such dispatch. + for k in range(n_batt_ABC - 1): + constraints.append(b[k] == b[k + 1]) + + # Constraints for flow in phase ABC + constraints.append(p_flow_A == sum_flow_A + sum(P_batt_A[k] for k in range(n_batt))) + constraints.append(p_flow_B == sum_flow_B + sum(P_batt_B[k] for k in range(n_batt))) + constraints.append(p_flow_C == sum_flow_C + sum(P_batt_C[k] for k in range(n_batt))) + + # Individual modulus is needed to make sure one phase doesn't compensate for other + constraints.append(p_flow_mod_A >= p_flow_A) + constraints.append(p_flow_mod_A >= - p_flow_A) + + constraints.append(p_flow_mod_B >= p_flow_B) + constraints.append(p_flow_mod_B >= - p_flow_B) + + constraints.append(p_flow_mod_C >= p_flow_C) + constraints.append(p_flow_mod_C >= - p_flow_C) + + # Objective function and invoke a solver + objective = (p_flow_mod_A + p_flow_mod_B + p_flow_mod_C) + problem = cp.Problem(cp.Minimize(objective), constraints) + problem.solve() + + # Extract optimization solution + idx = 0 + dispatch_batteries = {} + optimization_solution_table = [] + for batt in self.Battery: + name = self.Battery[batt]['name'] + dispatch_batteries[batt] = {} + dispatch_batteries[batt]['p_batt'] = p_batt[idx].value * 1000 + optimization_solution_table.append([name, self.Battery[batt]['phases'], p_batt[idx].value]) + data = [timestamp, name, self.Battery[batt]['phases'], p_batt[idx].value] + header = ["time", "battery", "phases", "p_batt"] + add_data_to_csv("optimization_result.csv", data, header=header) + idx += 1 + print('Optimization Solution') + print(tabulate(optimization_solution_table, headers=['Battery', 'phases', 'P_batt (kW)'], tablefmt='psql')) + + load_pv = ['{:.3f}'.format(sum_flow_A), '{:.3f}'.format(sum_flow_B), '{:.3f}'.format(sum_flow_C)] + load_pv_batt = ['{:.3f}'.format(p_flow_A.value), '{:.3f}'.format(p_flow_B.value), + '{:.3f}'.format(p_flow_C.value)] + optimization_summary = [] + optimization_summary.append([load_pv, load_pv_batt, problem.status]) + data = [timestamp] + data.extend(load_pv) + data.extend(load_pv_batt) + data.append(problem.status) + header = ["time", "load_pv_a", "load_pv_b", "load_pv_c", "load_pv_batt_a", "load_pv_batt_b", "load_pv_batt_c", "status"] + add_data_to_csv("optimization_summary.csv", data, header=header) + print(tabulate(optimization_summary, headers=['Load+PV (kW)', 'Load+PV+Batt (kW)', 'Status'], tablefmt='psql')) + return dispatch_batteries + + def on_measurement(self, sim: Simulation, timestamp: dict, measurements: dict) -> None: + if not self.has_batteries: + return + if self._count % 10 == 0: + # Call function to extract simulation measurements from injection sources + self.find_injection(self.Battery, measurements) + self.find_injection(self.EnergyConsumer, measurements) + self.find_injection(self.Solar, measurements) + + # Display real-time battery injection and current SOC from measurements + simulation_table_batteries = [] + for batt in self.Battery: + name = self.Battery[batt]['name'] + phases = self.Battery[batt]['phases'] + data = [timestamp, name, phases] + simulation_table_batteries.append( + [name, phases, self.Battery[batt]['P_inj'], self.Battery[batt]['soc']]) + data.extend(self.Battery[batt]['P_inj']) + data.extend([self.Battery[batt]['soc']]) + header=["time","battery","phases","p_batt_a","p_batt_b","p_batt_c","soc"] + add_data_to_csv("simulation_table.csv", data, header=header) + print(f'\n.......Curren timestamp: {timestamp}.......\n') + print('Simulation Table') + print(tabulate(simulation_table_batteries, headers=['Battery', 'phases', 'P_batt (kW)', 'SOC'], + tablefmt='psql')) + + # Invoke optimization for given grid condition + dispatch_values = self.optimize_battery(timestamp) + + # Dispatch battery. Note that -ve values in forward difference means charging batteries + # Make necessary changes in sign convention from optimization values + for unit in dispatch_values: + self._init_batt_diff.add_difference(unit, 'PowerElectronicsConnection.p', + - dispatch_values[unit]['p_batt'], 0.0) + msg = self._init_batt_diff.get_message() + self.gad_obj.send(self._publish_to_topic, json.dumps(msg)) + self._count += 1 + + +def createSimulation(gad_obj: GridAPPSD, model_info: Dict[str, Any]) -> Simulation: + # if not isinstance(gad_obj, GridAPPSD): + # raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + # if not isinstance(model_info, dict): + # raise TypeError(f'The model_id must be a dictionary type.') + line_name = model_info.get('modelId') + subregion_name = model_info.get('subRegionId') + region_name = model_info.get('regionId') + sim_name = model_info.get('modelName') + # if line_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelId.') + # if subregion_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key subRegionId.') + # if region_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key regionId.') + # if sim_name is None: + # raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelName.') + # psc = PowerSystemConfig(Line_name=line_name, + # SubGeographicalRegion_name=subregion_name, + # GeographicalRegion_name=region_name) + # start_time = int(datetime.utcnow().replace(microsecond=0).timestamp()) + # sim_args = SimulationArgs( + # start_time=f'{start_time}', + # # duration=f'{24*3600}', + # duration=120, + # simulation_name=sim_name, + # run_realtime=False + # # pause_after_measurements=True + # ) + # sim_config = SimulationConfig(power_system_config=psc, simulation_config=sim_args) + # sim_obj = Simulation(gapps=gad_obj, run_config=sim_config) + system_config = PowerSystemConfig( + GeographicalRegion_name=region_name, + SubGeographicalRegion_name=subregion_name, + Line_name=line_name + ) + + model_config = ModelCreationConfig( + load_scaling_factor=1, + schedule_name="ieeezipload", + z_fraction=0, + i_fraction=1, + p_fraction=0, + randomize_zipload_fractions=False, + use_houses=False + ) + + start = datetime(2024, 1, 1) + epoch = datetime.timestamp(start) + duration = timedelta(minutes=7).total_seconds() + # using hadrcode epoch to ensure afternoon time + epoch = 1704207000 + + sim_args = SimulationArgs( + start_time=epoch, + duration=duration, + simulator="GridLAB-D", + timestep_frequency=1000, + timestep_increment=1000, + run_realtime=True, + simulation_name=sim_name, + power_flow_solver_method="NR", + model_creation_config=asdict(model_config) + ) + + sim_config = SimulationConfig( + power_system_config=asdict(system_config), + simulation_config=asdict(sim_args) + ) + + sim_obj = Simulation(gad_obj, asdict(sim_config)) + return sim_obj + + +def main(control_enabled: bool, start_simulations: bool, model_id: str = None): + if not isinstance(control_enabled, bool): + raise TypeError(f'Argument, control_enabled, must be a boolean.') + if not isinstance(start_simulations, bool): + raise TypeError(f'Argument, start_simulations, must be a boolean.') + if not isinstance(model_id, str) and model_id is not None: + raise TypeError( + f'The model id passed to the convervation voltage reduction application must be a string type or {None}.') + + cim_profile = 'rc4_2021' + cim = importlib.import_module('cimgraph.data_profile.' + cim_profile) + feeder = cim.Feeder(mRID=model_id) + params = ConnectionParameters( + url="http://localhost:8889/bigdata/namespace/kb/sparql", + cim_profile=cim_profile) + + gapps = GridAPPSD(username='system', password='manager') + bg = BlazegraphConnection(params) + network = FeederModel( + connection=bg, + container=feeder, + distributed=False) + utils.get_all_data(network) + + local_simulations = {} + app_instances = {'field_instances': {}, 'external_simulation_instances': {}, 'local_simulation_instances': {}} + response = gapps.query_model_info() + models = response.get('data', {}).get('models', []) + model_is_valid = False + if model_id is not None: + for m in models: + m_id = m.get('modelId') + if model_id == m_id: + model_is_valid = True + break + if not model_is_valid: + raise ValueError(f'The model id provided does not exist in the GridAPPS-D plaform.') + else: + models = [m] + if start_simulations: + for m in models: + local_simulations[m.get('modelId', '')] = createSimulation(gapps, m) + else: + # TODO: query platform for running simulations which is currently not implemented in the GridAPPS-D Api + pass + # Create an cvr controller instance for all the real systems in the database + # for m in models: + # m_id = m.get('modelId') + # app_instances['field_instances'][m_id] = ConservationVoltageReductionController(gad_object, m_id) + + for m_id, simulation in local_simulations.items(): + measurements_topic = topics.simulation_output_topic(simulation) + simulation.start_simulation() + c_mapp = CarbonManagementApp( + gapps, m_id, network, simulation=simulation) + # gapps.subscribe(measurements_topic, c_mapp) + try: + while True: + time.sleep(0.1) + except KeyboardInterrupt: + print("Simulation manually stopped.") + finally: + print(" -- Stopping simulation -- ") + simulation.stop() + + +def add_data_to_csv(file_path, data, header=None): + # Check if the file exists + file_exists = os.path.isfile(file_path) + + # Open the CSV file in append mode (or create it if it doesn't exist) + with open(file_path, mode='a', newline='') as file: + writer = csv.writer(file) + # If the file doesn't exist, write the header first + if not file_exists and header: + writer.writerow(header) + # Write the data (a list or tuple representing a row) + writer.writerow(data) + + +if __name__ == "__main__": + parser = ArgumentParser() + parser.add_argument('model_id', nargs='?', default=None, help='The unit_mrid of the cim model to perform cvr on.') + parser.add_argument('-s', + '--start_simulations', + action='store_true', + help='Flag to have application start simulations') + parser.add_argument('-d', + '--disable_control', + action='store_true', + help='Flag to disable control on startup by default.') + args = parser.parse_args() + # args.start_simulations = True + # args.model_id = '_EE71F6C9-56F0-4167-A14E-7F4C71F10EAA' + main(args.disable_control, args.start_simulations, args.model_id) diff --git a/micro-apps/__init__.py b/micro-apps/peak-shaving-app/__init__.py similarity index 100% rename from micro-apps/__init__.py rename to micro-apps/peak-shaving-app/__init__.py diff --git a/micro-apps/peak-shaving-app/peak_shaving_main.py b/micro-apps/peak-shaving-app/peak_shaving_main.py new file mode 100644 index 0000000..93f1a50 --- /dev/null +++ b/micro-apps/peak-shaving-app/peak_shaving_main.py @@ -0,0 +1,1679 @@ +import importlib +import logging +import math +import os +import csv +from cmath import exp +from argparse import ArgumentParser +from datetime import datetime, timedelta, timezone +from typing import Any, Dict +from uuid import uuid4 + +import cimgraph.utils as cimUtils +from cimgraph.databases import ConnectionParameters +from cimgraph.databases.blazegraph.blazegraph import BlazegraphConnection +from cimgraph.databases.gridappsd.gridappsd import GridappsdConnection +from cimgraph.models import FeederModel +from gridappsd import DifferenceBuilder, GridAPPSD, topics +from gridappsd.simulation import * +from gridappsd.utils import ProcessStatusEnum + +logging.basicConfig(format='%(asctime)s::%(levelname)s::%(name)s::%(filename)s::%(lineno)d::%(message)s', + filename='peak_shaving.log', + filemode='w', + level=logging.INFO, + encoding='utf-8') +logger = logging.getLogger(__name__) +cim = None +DB_CONNECTION = None +CIM_GRAPH_MODELS = {} + + +IEEE123_APPS = "_F49D1288-9EC6-47DB-8769-57E2B6EDB124" +ieee123_apps_feeder_head_measurement_mrids = [ + "_903a0e85-6a11-4b8e-82cf-163376df7764", + "_6d82c9d9-5f99-4356-8a8d-acda9cd6f17b", + "_971765b5-da69-4ea8-a0cb-d2b9613b8ee3" +] + +def add_data_to_csv(file_path, data, header=None): + # Check if the file exists + file_exists = os.path.isfile(file_path) + + # Open the CSV file in append mode (or create it if it doesn't exist) + with open(file_path, mode='a', newline='') as file: + writer = csv.writer(file) + # If the file doesn't exist, write the header first + if not file_exists and header: + writer.writerow(header) + # Write the data (a list or tuple representing a row) + writer.writerow(data) + + +def findMeasurement(measurementsList, type: str, phases): + rv = None + for measurement in measurementsList: + if measurement.measurementType == type and measurement.phases in phases: + rv = measurement + break + if rv is None: + raise RuntimeError(f'findMeasurement(): No measurement of type {type} exists for any of the given phases!.') + return rv + + +def findMeasurements(measurementsList, type: str): + rv = [] + for measurement in measurementsList: + if measurement.measurementType == type: + rv.append({'object': measurement, 'value': None}) + if not rv: + raise RuntimeError(f'findMeasurements(): No measurement of type {type} exists!.') + return rv + + +def findPrimaryPhase(cimObj): + ''' + Helper function for finding the primary phase an instance of cim.ConductingEquipment on the secondary + system is connected to. + ''' + if not isinstance(cimObj, cim.ConductingEquipment): + raise TypeError('findPrimaryPhase(): cimObj must be an instance of cim.ConductingEquipment!') + equipmentToCheck = [cimObj] + phaseCode = None + xfmr = None + i = 0 + while not xfmr and i < len(equipmentToCheck): + equipmentToAdd = [] + for eq in equipmentToCheck[i:]: + if isinstance(eq, cim.PowerTransformer): + xfmr = eq + break + else: + terminals = eq.Terminals + connectivityNodes = [] + for t in terminals: + if t.ConnectivityNode not in connectivityNodes: + connectivityNodes.append(t.ConnectivityNode) + for cn in connectivityNodes: + for t in cn.Terminals: + if t.ConductingEquipment not in equipmentToCheck and t.ConductingEquipment not in equipmentToAdd: + equipmentToAdd.append(t.ConductingEquipment) + i = len(equipmentToCheck) + equipmentToCheck.extend(equipmentToAdd) + if not xfmr: + raise RuntimeError('findPrimaryPhase(): no upstream centertapped transformer could be found for secondary ' + f'system object {cimObj.name}!') + for tank in xfmr.TransformerTanks: + if phaseCode is not None: + break + for tankEnd in tank.TransformerTankEnds: + if tankEnd.phases not in [ + cim.PhaseCode.none, cim.PhaseCode.s1, cim.PhaseCode.s12, cim.PhaseCode.s12N, cim.PhaseCode.s1N, + cim.PhaseCode.s2, cim.PhaseCode.s2N + ]: + phaseCode = tankEnd.phases + break + if not phaseCode: + raise RuntimeError('findPrimaryPhase(): the upstream centertapped transformer has no primary phase defined!?') + return phaseCode + + +def findFeederPowerRating(cimObj): + ''' + Helper function for finding the feeder's rated power from and instance of cim.EnergySource. + ''' + if not isinstance(cimObj, cim.EnergySource): + raise TypeError('findPrimaryPhase(): cimObj must be an instance of cim.EnergySource!') + equipmentToCheck = [cimObj] + xfmr = None + feederPowerRating = None + i = 0 + while not xfmr and i < len(equipmentToCheck): + equipmentToAdd = [] + for eq in equipmentToCheck[i:]: + if isinstance(eq, cim.PowerTransformer): + xfmr = eq + break + else: + terminals = eq.Terminals + connectivityNodes = [] + for t in terminals: + if t.ConnectivityNode not in connectivityNodes: + connectivityNodes.append(t.ConnectivityNode) + for cn in connectivityNodes: + for t in cn.Terminals: + if t.ConductingEquipment not in equipmentToCheck and t.ConductingEquipment not in equipmentToAdd: + equipmentToAdd.append(t.ConductingEquipment) + i = len(equipmentToCheck) + equipmentToCheck.extend(equipmentToAdd) + if not xfmr: + raise RuntimeError('findFeederPowerRating(): No feeder head transformer could be found for EnergySource, ' + f'{cimObj.name}!') + powerTransformerEnds = xfmr.PowerTransformerEnd + if powerTransformerEnds: + feederPowerRating = float(powerTransformerEnds[0].ratedS) + else: + raise RuntimeError('findFeederPowerRating(): The found at the feeder head is not a three phase transformer!') + logger.info(f'Feeder Power Rating: {feederPowerRating} W') + return feederPowerRating + + +def findFeederHeadPowerMeasurmentsObject(cimObj: object): + ''' + Helper function to find feeder transformer + ''' + if not isinstance(cimObj, cim.EnergySource): + raise TypeError('findFeederHeadPowerMeasurmentsObject(): cimObj must be an instance of cim.EnergySource!') + equipmentToCheck = [cimObj] + feederHeadLoadMeasurementsObject = None + i = 0 + while not feederHeadLoadMeasurementsObject and i < len(equipmentToCheck): + equipmentToAdd = [] + for eq in equipmentToCheck[i:]: + powerMeasurementsCount = 0 + for meas in eq.Measurements: + if meas.measurementType == 'VA': + powerMeasurementsCount += 1 + if powerMeasurementsCount == 3: + feederHeadLoadMeasurementsObject = eq + break + else: + terminals = eq.Terminals + connectivityNodes = [] + for t in terminals: + if t.ConnectivityNode not in connectivityNodes: + connectivityNodes.append(t.ConnectivityNode) + for cn in connectivityNodes: + for t in cn.Terminals: + if t.ConductingEquipment not in equipmentToCheck and t.ConductingEquipment not in equipmentToAdd: + equipmentToAdd.append(t.ConductingEquipment) + i = len(equipmentToCheck) + equipmentToCheck.extend(equipmentToAdd) + if not feederHeadLoadMeasurementsObject: + raise RuntimeError('findFeederPowerRating(): No feeder head object with power measurements could be found for ' + f'EnergySource, {cimObj.name}!') + return feederHeadLoadMeasurementsObject + + +def buildGraphModel(mrid: str) -> FeederModel: + global CIM_GRAPH_MODELS + if not isinstance(mrid, str): + raise TypeError(f'The mrid passed to the convervation voltage reduction application must be a string.') + if mrid not in CIM_GRAPH_MODELS.keys(): + feeder_container = cim.Feeder(mRID=mrid) + graph_model = FeederModel(connection=DB_CONNECTION, container=feeder_container, distributed=False) + cimUtils.get_all_data(graph_model) + CIM_GRAPH_MODELS[mrid] = graph_model + return CIM_GRAPH_MODELS[mrid] + + +def createSimulation(gad_obj: GridAPPSD, model_info: Dict[str, Any]) -> Simulation: + if not isinstance(gad_obj, GridAPPSD): + raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + if not isinstance(model_info, dict): + raise TypeError(f'The model_id must be a dictionary type.') + line_name = model_info.get('modelId') + subregion_name = model_info.get('subRegionId') + region_name = model_info.get('regionId') + sim_name = model_info.get('modelName') + if line_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelId.') + if subregion_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key subRegionId.') + if region_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key regionId.') + if sim_name is None: + raise ValueError(f'Bad model info dictionary. The dictionary is missing key modelName.') + psc = PowerSystemConfig(Line_name=line_name, + SubGeographicalRegion_name=subregion_name, + GeographicalRegion_name=region_name) + start_time = int(datetime.utcnow().replace(microsecond=0).timestamp()) + start = datetime(2023, 1, 1, 0, tzinfo=timezone.utc) + epoch = datetime.timestamp(start) + duration = timedelta(minutes=1).total_seconds() + sim_args = SimulationArgs(start_time=f'{epoch}', + duration=f'{duration}', + simulation_name=sim_name, + run_realtime=True, + pause_after_measurements=False) + sim_config = SimulationConfig(power_system_config=psc, simulation_config=sim_args) + sim_obj = Simulation(gapps=gad_obj, run_config=sim_config) + return sim_obj + + +def createGadObject() -> GridAPPSD: + gad_user = os.environ.get('GRIDAPPSD_USER') + if gad_user is None: + os.environ['GRIDAPPSD_USER'] = 'system' + gad_password = os.environ.get('GRIDAPPSD_PASSWORD') + if gad_password is None: + os.environ['GRIDAPPSD_PASSWORD'] = 'manager' + gad_app_id = os.environ.get('GRIDAPPSD_APPLICATION_ID') + if gad_app_id is None: + os.environ['GRIDAPPSD_APPLICATION_ID'] = 'PeakShavingApplication' + return GridAPPSD() + + +def removeDirectory(directory: Path | str): + if isinstance(directory, str): + directory = Path(directory) + for item in directory.iterdir(): + if item.is_dir(): + removeDirectory(item) + else: + item.unlink() + directory.rmdir() + + +def pol2cart(mag, angle_deg): + # Convert degrees to radians. GridAPPS-D spits angle in degrees + angle_rad = math.radians(angle_deg) + p = mag * math.cos(angle_rad) + q = mag * math.sin(angle_rad) + return p, q + + +class PeakShavingController(object): + """Peak Shaving Control Class + This class implements a centralized algorithm for feeder peak shaving by controlling batteries in the + distribution model. + + Attributes: + sim_id: The simulation id that the instance will be perfoming peak shaving on. If the sim_id is None then the + application will assume it is performing peak shaving on the actual field devices. + Type: str. + Default: None. + Methods: + on_measurement(headers:Dict[Any], message:Dict[Any]): This callback function will be used to update the + applications measurements dictionary needed for peak shaving control. + Arguments: + headers: A dictionary containing header information on the message received from the GridAPPS-D + platform. + Type: dictionary. + Default: NA. + message: A dictionary containing the message received from the GridAPPS-D platform. + Type: dictionary. + Default: NA. + Returns: NA. + peak_shaving_control(): This is the main function for performing the peak shaving control. + """ + + def __init__(self, + gad_obj: GridAPPSD, + model_id: str, + peak_setpoint: float = None, + base_setpoint: float = None, + sim_id: str = None, + simulation: Simulation = None): + if not isinstance(gad_obj, GridAPPSD): + raise TypeError(f'gad_obj must be an instance of GridAPPSD!') + if not isinstance(model_id, str): + raise TypeError(f'model_id must be a str type.') + if model_id is None or model_id == '': + raise ValueError(f'model_id must be a valid uuid.') + if not isinstance(peak_setpoint, float) and peak_setpoint is not None: + raise TypeError(f'peak_setpoint must be of float type or {None}!') + if not isinstance(base_setpoint, float) and base_setpoint is not None: + raise TypeError(f'base_setpoint must be of float type or {None}!') + if not isinstance(sim_id, str) and sim_id is not None: + raise TypeError(f'sim_id must be a string type or {None}!') + if not isinstance(simulation, Simulation) and simulation is not None: + raise TypeError(f'The simulation arg must be a Simulation type or {None}!') + self.model_results_path = None + self.timestamp = 0 + self.model_id = model_id + self.id = uuid4() + self.desired_setpoints = {} + self.pnv_measurements = {} + self.pnv_measurements_pu = {} + self.va_measurements = {} + self.pos_measurements = {} + self.link_measurements = {} + self.controllable_batteries_ABC = {} + self.controllable_batteries_A = {} + self.controllable_batteries_B = {} + self.controllable_batteries_C = {} + self.peak_va_measurements_A = {} + self.peak_va_measurements_B = {} + self.peak_va_measurements_C = {} + self.peak_setpoint_A = None + self.peak_setpoint_B = None + self.peak_setpoint_C = None + self.base_setpoint_A = None + self.base_setpoint_B = None + self.base_setpoint_C = None + self.measurements_topic = None + self.setpoints_topic = None + self.simulation = None + self.simulation_id = None # TODO: figure out what simulation_id should be when deployed in the field. + if simulation is not None: + self.simulation = simulation + self.simulation_id = simulation.simulation_id + self.simulation.add_onmeasurement_callback(self.on_measurement) + if sim_id is None: + self.measurements_topic = topics.field_output_topic(None, sim_id) + else: + self.simulation_id = sim_id + self.measurements_topic = topics.simulation_output_topic(sim_id) + + self.gad_obj = gad_obj + self.log = self.gad_obj.get_logger() + if self.measurements_topic: + self.gad_obj.subscribe(self.measurements_topic, self.on_measurement_callback) + if self.simulation_id: + self.setpoints_topic = topics.simulation_input_topic(self.simulation_id) + else: + self.setpoints_topic = topics.field_input_topic() + # Read model_id from cimgraph to get all the controllable batteries, and measurements. + self.graph_model = buildGraphModel(model_id) + self.installed_battery_capacity_ABC = 0.0 + self.installed_battery_capacity_A = 0.0 + self.installed_battery_capacity_B = 0.0 + self.installed_battery_capacity_C = 0.0 + self.installed_battery_power_ABC = 0.0 + self.installed_battery_power_A = 0.0 + self.installed_battery_power_B = 0.0 + self.installed_battery_power_C = 0.0 + self.configureBatteryProperties() + # feederLoadMeasurements(measurements.keys()) + if peak_setpoint is not None: + self.peak_setpoint_A = peak_setpoint / 3.0 + self.peak_setpoint_B = peak_setpoint / 3.0 + self.peak_setpoint_C = peak_setpoint / 3.0 + else: + self.configurePeakShavingSetpoint() + if base_setpoint is not None: + self.base_setpoint_A = base_setpoint / 3.0 + self.base_setpoint_B = base_setpoint / 3.0 + self.base_setpoint_C = base_setpoint / 3.0 + else: + self.base_setpoint_A = 0.8 * self.peak_setpoint_A + self.base_setpoint_B = 0.8 * self.peak_setpoint_B + self.base_setpoint_C = 0.8 * self.peak_setpoint_C + self.findFeederHeadLoadMeasurements() + + + measurements = {} + measurements.update(self.graph_model.graph.get(cim.Analog, {})) + measurements.update(self.graph_model.graph.get(cim.Discrete, {})) + + # print(measurements.keys()) + for meas in measurements.values(): + if meas.measurementType == 'PNV': # it's a voltage measurement. store it. + mrid = meas.mRID + self.pnv_measurements[mrid] = {'measurement_object': meas, 'measurement_value': None} + elif meas.measurementType == 'VA': # it's a power measurement. + if isinstance(meas.PowerSystemResource, (cim.EnergyConsumer, cim.PowerElectronicsConnection)): + mrid = meas.PowerSystemResource.mRID + if mrid not in self.va_measurements.keys(): + self.va_measurements[mrid] = {'measurement_objects': {}, 'measurement_values': {}} + self.va_measurements[mrid]['measurement_objects'][meas.mRID] = meas + self.va_measurements[mrid]['measurement_values'][meas.mRID] = None + if isinstance(meas.PowerSystemResource, (cim.ACLineSegment, cim.PowerTransformer)): + mrid = meas.PowerSystemResource.mRID + if mrid not in self.va_measurements.keys(): + self.link_measurements[mrid] = {'measurement_objects': {}, 'measurement_values': {}} + self.link_measurements[mrid]['measurement_objects'][meas.mRID] = meas + self.link_measurements[mrid]['measurement_values'][meas.mRID] = None + elif meas.measurementType == 'Pos': # it's a positional measurement. + if isinstance(meas.PowerSystemResource, (cim.Switch, cim.PowerTransformer, cim.LinearShuntCompensator)): + mrid = meas.PowerSystemResource.mRID + if mrid not in self.pos_measurements.keys(): + self.pos_measurements[mrid] = {'measurement_objects': {}, 'measurement_values': {}} + self.pos_measurements[mrid]['measurement_objects'][meas.mRID] = meas + self.pos_measurements[mrid]['measurement_values'][meas.mRID] = None + + # data output + out_dir = Path(__file__).parent/"output" + feeder = self.graph_model.graph.get(cim.Feeder, {}).get(self.model_id) + self.model_results_path = out_dir/feeder.name + self.model_results_path.mkdir(parents=True, exist_ok=True) + # self.findFeederHeadLoadMeasurements() + self.init_output_file() + if self.simulation is not None: + self.simulation.start_simulation() + self.simulation_id = self.simulation.simulation_id + self.differenceBuilder = DifferenceBuilder(self.simulation_id) + self.isValid = True + self.first_message = True + + def init_output_file(self): + for child in self.model_results_path.iterdir(): + if not child.is_dir(): + child.unlink() + with open(self.model_results_path / "voltages.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "node", "phase", "voltage"] + writer.writerow(header) + with open(self.model_results_path / "total_load.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "phase", "p", "q"] + writer.writerow(header) + with open(self.model_results_path / "feeder_head.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "phase", "p", "q"] + writer.writerow(header) + with open(self.model_results_path / "battery_power.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "battery", "p_a", "p_b", "p_c", "soc"] + writer.writerow(header) + with open(self.model_results_path / "dispatches.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "mrid", "name", "phase", "value"] + writer.writerow(header) + with open(self.model_results_path / "loads_no_batt.csv", mode='a', newline='') as file: + writer = csv.writer(file) + header = ["time", "a", "b", "c"] + writer.writerow(header) + + def get_feeder_head_measurements(self, measurements): + total_loads = {} + if self.graph_model.container.mRID == IEEE123_APPS: + for mrid in ieee123_apps_feeder_head_measurement_mrids: + phase = self.graph_model.graph.get(cim.Analog, {}).get(mrid).phases.value + s_magnitude = measurements.get(mrid).get("magnitude") + s_angle_rad = measurements.get(mrid).get("angle") + total_loads[phase] = s_magnitude*exp(1j*math.radians(s_angle_rad)) + return total_loads + + def on_measurement(self, sim: Simulation, timestamp: str, measurements: Dict[str, Dict]): + self.timestamp = timestamp + self.desired_setpoints.clear() + if isinstance(sim, Simulation): + self.setpoints_topic = topics.simulation_input_topic(sim.simulation_id) + sim.pause() + total_loads = self.get_feeder_head_measurements(measurements) + self.save_total_load_data(timestamp, total_loads) + for mrid in self.peak_va_measurements_A.keys(): + measurement = measurements.get(self.peak_va_measurements_A[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_A[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "A", load.real, load.imag]) + for mrid in self.peak_va_measurements_B.keys(): + measurement = measurements.get(self.peak_va_measurements_B[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_B[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "B", load.real, load.imag]) + for mrid in self.peak_va_measurements_C.keys(): + measurement = measurements.get(self.peak_va_measurements_C[mrid]['object'].mRID) + if measurement is not None: + self.peak_va_measurements_C[mrid]['value'] = measurement + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load = mag * exp(1j*math.radians(ang_in_deg)) + add_data_to_csv(self.model_results_path/"feeder_head.csv", [self.timestamp, mrid, "C", load.real, load.imag]) + for mrid in self.controllable_batteries_ABC.keys(): + power_data_to_save = [0, 0, 0] + for measDict in self.controllable_batteries_ABC[mrid]['power_measurements']: + measurement = measurements.get(measDict['object'].mRID) + if measurement is not None: + measDict['value'] = measurement + phase = measDict["object"].phases.value + p, q = pol2cart(measurement.get("magnitude"), measurement.get("angle")) + power_data_to_save[("ABC").index(phase)] = p + measurement = measurements.get(self.controllable_batteries_ABC[mrid]['soc_measurement']['object'].mRID) + if measurement is not None: + self.controllable_batteries_ABC[mrid]['soc_measurement']['value'] = measurement + header = ["time", "mrid", "battery", "p_a", "p_b", "p_c", "soc"] + data = [self.timestamp, mrid, measDict['object'].name] + data.extend(power_data_to_save) + data.append(self.controllable_batteries_ABC[mrid]['soc_measurement']['value'].get("value")) + add_data_to_csv(self.model_results_path/"battery_power.csv", data) + for mrid in self.controllable_batteries_A.keys(): + for measDict in self.controllable_batteries_A[mrid]['power_measurements']: + measurement = measurements.get(measDict['object'].mRID) + if measurement is not None: + measDict['value'] = measurement + measurement = measurements.get(self.controllable_batteries_A[mrid]['soc_measurement']['object'].mRID) + if measurement is not None: + self.controllable_batteries_A[mrid]['soc_measurement']['value'] = measurement + for mrid in self.controllable_batteries_B.keys(): + for measDict in self.controllable_batteries_B[mrid]['power_measurements']: + measurement = measurements.get(measDict['object'].mRID) + if measurement is not None: + measDict['value'] = measurement + measurement = measurements.get(self.controllable_batteries_B[mrid]['soc_measurement']['object'].mRID) + if measurement is not None: + self.controllable_batteries_B[mrid]['soc_measurement']['value'] = measurement + for mrid in self.controllable_batteries_C.keys(): + for measDict in self.controllable_batteries_C[mrid]['power_measurements']: + measurement = measurements.get(measDict['object'].mRID) + if measurement is not None: + measDict['value'] = measurement + measurement = measurements.get(self.controllable_batteries_C[mrid]['soc_measurement']['object'].mRID) + if measurement is not None: + self.controllable_batteries_C[mrid]['soc_measurement']['value'] = measurement + + for mrid in self.pnv_measurements.keys(): + meas = measurements.get(mrid) + if meas is not None: + self.pnv_measurements[mrid]['measurement_value'] = meas + for psr_mrid in self.va_measurements.keys(): + for mrid in self.va_measurements[psr_mrid]['measurement_values'].keys(): + meas = measurements.get(mrid) + if meas is not None: + self.va_measurements[psr_mrid]['measurement_values'][mrid] = meas + for psr_mrid in self.pos_measurements.keys(): + for mrid in self.pos_measurements[psr_mrid]['measurement_values'].keys(): + meas = measurements.get(mrid) + if meas is not None: + self.pos_measurements[psr_mrid]['measurement_values'][mrid] = meas + if isinstance(self.pos_measurements[psr_mrid]['measurement_objects'][mrid].PowerSystemResource, (cim.LinearShuntCompensator, cim.PowerTransformer)): + pos_measurement_debug_data = { + "object_type": f"{type(self.pos_measurements[psr_mrid]['measurement_objects'][mrid].PowerSystemResource).__name__}", + "name": self.pos_measurements[psr_mrid]['measurement_objects'][mrid].PowerSystemResource.name, + "phase": self.pos_measurements[psr_mrid]['measurement_objects'][mrid].phases.value, + "value": int(meas.get("value")), + "timestamp": timestamp + } + logger.info(f"Pos Measurement: {json.dumps(pos_measurement_debug_data)}") + + self.calculate_per_unit_voltage() + self.save_voltage_data(timestamp) + self.peak_shaving_control() + if isinstance(sim, Simulation): + sim.resume() + + def findFeederHeadLoadMeasurements(self): + energySources = self.graph_model.graph.get(cim.EnergySource, {}) + feederLoadObject = None + for eSource in energySources.values(): + feederLoadObject = findFeederHeadPowerMeasurmentsObject(eSource) + feederLoadMeasurements = feederLoadObject.Measurements + for measurement in feederLoadMeasurements: + if measurement.measurementType == 'VA': + if measurement.phases in [cim.PhaseCode.A, cim.PhaseCode.AN]: + self.peak_va_measurements_A[measurement.mRID] = {'object': measurement, 'value': None} + elif measurement.phases in [cim.PhaseCode.B, cim.PhaseCode.BN]: + self.peak_va_measurements_B[measurement.mRID] = {'object': measurement, 'value': None} + elif measurement.phases in [cim.PhaseCode.C, cim.PhaseCode.CN]: + self.peak_va_measurements_C[measurement.mRID] = {'object': measurement, 'value': None} + if not self.peak_va_measurements_A or not self.peak_va_measurements_B or not self.peak_va_measurements_C: + raise RuntimeError(f'feeder {self.graph_model.container.mRID}, has no measurements associated with the ' + 'feeder head transformer!') + + def calculate_per_unit_voltage(self): + for mrid in self.pnv_measurements.keys(): + self.pnv_measurements_pu[mrid] = None + meas_base = None + meas = self.pnv_measurements.get(mrid) + if (meas is None) or (meas.get('measurement_value') is None): + self.log.warning(f'Measurement {mrid} is missing. It will be ignored in cvr control.') + continue + meas_value = meas.get('measurement_value').get('magnitude') + if meas_value is None: + self.log.error(f'An unexpected value of None was recieved for PNV measurement {mrid}. It will be ' + 'ignored in cvr control.') + continue + meas_obj = meas.get('measurement_object') + if (meas_obj is None) or (not isinstance(meas_obj, cim.Measurement)): + self.log.error(f'The cim.Measurement object {mrid} associated with this measurement value is missing! ' + 'It will be ignored in cvr control.') + logger.error(f'The measurement dictionary for mrid {mrid} is missing from the CIM database.') + continue + + meas_term = meas_obj.Terminal + if isinstance(meas_obj.PowerSystemResource, cim.PowerElectronicsConnection): + meas_base = float(meas_obj.PowerSystemResource.ratedU) + elif isinstance(meas_term, cim.Terminal): + if meas_term.TransformerEnd: + if isinstance(meas_term.TransformerEnd[0], cim.PowerTransformerEnd): + meas_base = float(meas_term.TransformerEnd[0].ratedU) + elif isinstance(meas_term.TransformerEnd[0], cim.TransformerTankEnd): + if isinstance(meas_term.TranformerEnd[0].BaseVoltage, cim.BaseVoltage): + meas_base = float(meas_term.TransformerEnd[0].BaseVoltage.nominalVoltage) + else: + self.log.error('meas_term.TransformerEnd[0].BaseVoltage is None') + logger.error('meas_term.TransformerEnd[0].BaseVoltage is None') + raise RuntimeError('PNV measurement BaseVoltage for ConnectivityNode ' + f'{meas_term.ConnectivityNode.name} is None') + elif isinstance(meas_term.ConductingEquipment, cim.ConductingEquipment): + if isinstance(meas_term.ConductingEquipment.BaseVoltage, cim.BaseVoltage): + meas_base = float(meas_term.ConductingEquipment.BaseVoltage.nominalVoltage) + else: + logger.error(f'meas_term.ConductingEquipment.BaseVoltage is None') + else: + logger.error(f'meas_term.ConductingEquipment is None') + if (meas_base is None) or (meas_base < 1e-10): + self.log.error(f'Unable to get the nominal voltage for measurement with mrid {mrid}.') + logger.error('Voltage Measurement has no associated nominal Voltage.\nMeasurement: ' + f'{meas_obj.name}\nTerminal: {meas_obj.Terminal.name}\n') + continue + pnv_meta_data = {'measured PNV': meas_value, 'measured LLV': meas_value * math.sqrt(3.0), 'rated V': meas_base} + pu_vals = [meas_value / meas_base, meas_value * math.sqrt(3.0) / meas_base] + diff = max(pu_vals) + for i in range(len(pu_vals)): + if abs(1.0-pu_vals[i]) < diff: + diff = abs(1.0-pu_vals[i]) + self.pnv_measurements_pu[mrid] = pu_vals[i] + pnv_meta_data['pu V'] = self.pnv_measurements_pu[mrid] + logger.info(f'Voltage meta data for {meas_obj.Terminal.ConnectivityNode.name}:{json.dumps(pnv_meta_data,indent=4)}') + + def save_total_load_data(self, time, total_loads): + file_path = self.model_results_path/"total_load.csv" + header = ["time", "phase", "p", "q"] + for phase, load in total_loads.items(): + data = [time, phase, load.real, load.imag] + add_data_to_csv(file_path, data, header=header) + + def save_voltage_data(self, time): + file_path = self.model_results_path/"voltages.csv" + header = ["time", "mrid", "node", "phase", "voltage"] + for mrid in self.pnv_measurements.keys(): + v = self.pnv_measurements_pu[mrid] + phase = self.pnv_measurements[mrid]["measurement_object"].phases.value + node = self.pnv_measurements[mrid]["measurement_object"].Terminal.ConnectivityNode.name + if v is not None: + data = [time, mrid, node, phase, v] + add_data_to_csv(file_path, data, header=header) + + def configureBatteryProperties(self): + '''This function uses cimgraph to populate the following properties: + + self.controllable_batteries_A + self.controllable_batteries_B + self.controllable_batteries_C + ''' + powerElectronicsConnections = self.graph_model.graph.get(cim.PowerElectronicsConnection, {}) + for pec in powerElectronicsConnections.values(): + isBatteryInverter = False + batteryCapacity = 0.0 + for powerElectronicsUnit in pec.PowerElectronicsUnit: + if isinstance(powerElectronicsUnit, cim.BatteryUnit): + isBatteryInverter = True + batteryCapacity += float(powerElectronicsUnit.ratedE) + if isBatteryInverter: + inverterPhases = self.getInverterPhases(pec) + if inverterPhases in [cim.PhaseCode.A, cim.PhaseCode.AN]: + self.installed_battery_capacity_A += batteryCapacity + self.installed_battery_power_A += float(pec.ratedS) + self.controllable_batteries_A[pec.mRID] = { + 'object': pec, + 'maximum_power': float(pec.ratedS), + 'power_measurements': findMeasurements(pec.Measurements, 'VA'), + 'soc_measurement': { + 'object': findMeasurement(pec.Measurements, 'SoC', [cim.PhaseCode.none]), + 'value': None + } + } + elif inverterPhases in [cim.PhaseCode.B, cim.PhaseCode.BN]: + self.installed_battery_capacity_B += batteryCapacity + self.installed_battery_power_B += float(pec.ratedS) + self.controllable_batteries_B[pec.mRID] = { + 'object': pec, + 'maximum_power': float(pec.ratedS), + 'power_measurements': findMeasurements(pec.Measurements, 'VA'), + 'soc_measurement': { + 'object': findMeasurement(pec.Measurements, 'SoC', [cim.PhaseCode.B, cim.PhaseCode.BN]), + 'value': None + } + } + elif inverterPhases in [cim.PhaseCode.C, cim.PhaseCode.CN]: + self.installed_battery_capacity_C += batteryCapacity + self.installed_battery_power_C += float(pec.ratedS) + self.controllable_batteries_C[pec.mRID] = { + 'object': pec, + 'maximum_power': float(pec.ratedS), + 'power_measurements': findMeasurements(pec.Measurements, 'VA'), + 'soc_measurement': { + 'object': findMeasurement(pec.Measurements, 'SoC', [cim.PhaseCode.C, cim.PhaseCode.CN]), + 'value': None + } + } + elif inverterPhases in [cim.PhaseCode.ABC, cim.PhaseCode.ABCN]: + self.installed_battery_capacity_ABC += batteryCapacity + self.installed_battery_power_ABC += float(pec.ratedS) + self.controllable_batteries_ABC[pec.mRID] = { + 'object': pec, + 'maximum_power': float(pec.ratedS), + 'power_measurements': findMeasurements(pec.Measurements, 'VA'), + 'soc_measurement': { + 'object': findMeasurement(pec.Measurements, 'SoC', [cim.PhaseCode.none]), + 'value': None + } + } + + def getInverterPhases(self, cimObj): + # algorithm that attempts to find the upstream centertapped transformer feeding secondary system inverters. + phaseCode = None + if not isinstance(cimObj, cim.PowerElectronicsConnection): + raise TypeError('PeakShavingController.getInverterPhases(): cimObj must be an instance of ' + 'cim.PowerElectronicsConnection!') + strPhases = '' + if len(cimObj.PowerElectronicsConnectionPhases) > 0: + for pecp in cimObj.PowerElectronicsConnectionPhases: + if pecp.phase != cim.SinglePhaseKind.s1 and pecp.phase != cim.SinglePhaseKind.s2: + strPhases += pecp.phase.value + else: + strPhases = pecp.phase.value + else: + strPhases = 'ABC' + if 'A' in strPhases or 'B' in strPhases or 'C' in strPhases: + phaseCodeStr = '' + if 'A' in strPhases: + phaseCodeStr += 'A' + if 'B' in strPhases: + phaseCodeStr += 'B' + if 'C' in strPhases: + phaseCodeStr += 'C' + phaseCode = cim.PhaseCode(phaseCodeStr) + else: # inverter is on the secondary system need to trace up to the centertapped tansformer. + phaseCode = self.findPrimaryPhase(cimObj) + return phaseCode + + def configurePeakShavingSetpoint(self): + energySources = self.graph_model.graph.get(cim.EnergySource, {}) + feederPowerRating = 0.0 + for source in energySources.values(): + feederPowerRating = findFeederPowerRating(source) + self.peak_setpoint_A = max(0.5 * (feederPowerRating / 3.0), + (feederPowerRating / 3.0) - (self.installed_battery_capacity_A * 0.95)) + self.peak_setpoint_B = max(0.5 * (feederPowerRating / 3.0), + (feederPowerRating / 3.0) - (self.installed_battery_capacity_B * 0.95)) + self.peak_setpoint_C = max(0.5 * (feederPowerRating / 3.0), + (feederPowerRating / 3.0) - (self.installed_battery_capacity_C * 0.95)) + + def on_measurement_callback(self, header: Dict[str, Any], message: Dict[str, Any]): + timestamp = message.get('message', {}).get('timestamp', '') + measurements = message.get('message', {}).get('measurements', {}) + self.on_measurement(None, timestamp, measurements) + + def peak_shaving_control(self): + lower_limit = 20.0 + upper_limit = 80.0 + real_load_A, real_load_B, real_load_C = self.get_load_minus_batteries() + add_data_to_csv(self.model_results_path/"loads_no_batt.csv", [self.timestamp, real_load_A, real_load_B, real_load_C]) + power_diff_A = 0.0 + power_diff_B = 0.0 + power_diff_C = 0.0 + deadband_ABC = True + if real_load_A > self.peak_setpoint_A: + power_diff_A = real_load_A - self.peak_setpoint_A + elif real_load_A < self.base_setpoint_A: + power_diff_A = real_load_A - self.base_setpoint_A + if real_load_B > self.peak_setpoint_B: + power_diff_B = real_load_B - self.peak_setpoint_B + elif real_load_B < self.base_setpoint_B: + power_diff_B = real_load_B - self.base_setpoint_B + if real_load_C > self.peak_setpoint_C: + power_diff_C = real_load_C - self.peak_setpoint_C + elif real_load_C < self.base_setpoint_C: + power_diff_C = real_load_C - self.base_setpoint_C + min_power_diff = min(power_diff_A, power_diff_B, power_diff_C) + max_power_diff = max(power_diff_A, power_diff_B, power_diff_C) + if min_power_diff > 1e-6: + deadband_ABC = False + control_dict, actual_power = self.calc_batt_discharge_ABC(3.0 * min_power_diff, lower_limit) + self.desired_setpoints.update(control_dict) + power_diff_A -= actual_power / 3.0 + power_diff_B -= actual_power / 3.0 + power_diff_C -= actual_power / 3.0 + elif max_power_diff < -1e-6: + deadband_ABC = False + control_dict, actual_power = self.calc_batt_charge_ABC(3.0 * abs(max_power_diff), upper_limit) + self.desired_setpoints.update(control_dict) + power_diff_A += abs(actual_power) / 3.0 + power_diff_B += abs(actual_power) / 3.0 + power_diff_C += abs(actual_power) / 3.0 + if power_diff_A > 1e-6: + control_dict = self.calc_batt_discharge_A(power_diff_A, lower_limit) + self.desired_setpoints.update(control_dict) + elif power_diff_A < -1e-6: + control_dict = self.calc_batt_charge_A(abs(power_diff_A), upper_limit) + self.desired_setpoints.update(control_dict) + if power_diff_B > 1e-6: + control_dict = self.calc_batt_discharge_B(power_diff_B, lower_limit) + self.desired_setpoints.update(control_dict) + elif power_diff_B < -1e-6: + control_dict = self.calc_batt_charge_B(abs(power_diff_B), upper_limit) + self.desired_setpoints.update(control_dict) + if power_diff_C > 1e-6: + control_dict = self.calc_batt_discharge_C(power_diff_C, lower_limit) + self.desired_setpoints.update(control_dict) + elif power_diff_C < -1e-6: + control_dict = self.calc_batt_charge_C(abs(power_diff_C), upper_limit) + self.desired_setpoints.update(control_dict) + if deadband_ABC: + for batt_id in self.controllable_batteries_ABC.keys(): + measurements = self.controllable_batteries_ABC[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + break + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + if (abs(current_power) > 1e-6): + self.desired_setpoints[batt_id] = { + 'object': self.controllable_batteries_ABC[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + if self.desired_setpoints: + self.send_setpoints() + + def get_load_minus_batteries(self): + load_A = 0.0 + load_B = 0.0 + load_C = 0.0 + for xfmr_id in self.peak_va_measurements_A.keys(): + measurement = self.peak_va_measurements_A[xfmr_id].get('value') + if measurement is not None: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_A += mag * math.cos(math.radians(ang_in_deg)) + for xfmr_id in self.peak_va_measurements_B.keys(): + measurement = self.peak_va_measurements_B[xfmr_id].get('value') + if measurement is not None: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_B += mag * math.cos(math.radians(ang_in_deg)) + for xfmr_id in self.peak_va_measurements_C.keys(): + measurement = self.peak_va_measurements_C[xfmr_id].get('value') + if measurement is not None: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_C += mag * math.cos(math.radians(ang_in_deg)) + for batt_id_3ph in self.controllable_batteries_ABC.keys(): + for measDict in self.controllable_batteries_ABC[batt_id_3ph]['power_measurements']: + measurement = measDict.get('value') + measurementObj = measDict.get('object') + if (measurement is None) or (measurementObj is None): + continue + measurementPhase = measurementObj.phases + if measurementPhase in [cim.PhaseCode.A, cim.PhaseCode.AN]: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_A -= mag * math.cos(math.radians(ang_in_deg)) + elif measurementPhase in [cim.PhaseCode.B, cim.PhaseCode.BN]: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_B -= mag * math.cos(math.radians(ang_in_deg)) + elif measurementPhase in [cim.PhaseCode.C, cim.PhaseCode.CN]: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_C -= mag * math.cos(math.radians(ang_in_deg)) + for batt_id_1ph in self.controllable_batteries_A.keys(): + for measDict in self.controllable_batteries_A[batt_id_1ph]['power_measurements']: + measurement = measDict.get('value') + if measurement is not None: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_A -= mag * math.cos(math.radians(ang_in_deg)) + for batt_id_1ph in self.controllable_batteries_B.keys(): + for measDict in self.controllable_batteries_B[batt_id_1ph]['power_measurements']: + measurement = measDict.get('value') + if measurement is not None: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_B -= mag * math.cos(math.radians(ang_in_deg)) + for batt_id_1ph in self.controllable_batteries_C.keys(): + for measDict in self.controllable_batteries_C[batt_id_1ph]['power_measurements']: + measurement = measDict.get('value') + if measurement is not None: + mag = measurement.get('magnitude') + ang_in_deg = measurement.get('angle') + if ((mag is not None) and (ang_in_deg is not None)): + load_C -= mag * math.cos(math.radians(ang_in_deg)) + return (load_A, load_B, load_C) + + def calc_batt_discharge_ABC(self, power_to_discharge_ABC, lower_limit): + if not isinstance(power_to_discharge_ABC, float): + raise TypeError('calc_batt_discharge_ABC(): power_to_discharge_ABC must be an instance of float!') + if power_to_discharge_ABC < 0.0: + raise ValueError(f'calc_batt_discharge_ABC(): power_to_discharge_ABC must be nonnegative!') + if not isinstance(lower_limit, float): + raise TypeError('calc_batt_discharge_ABC(): lower_limit must be an instance of float!') + if (lower_limit > 100.0) or (lower_limit < 0.0): + raise ValueError(f'calc_batt_discharge_ABC(): lower_limit must belong to the [0, 1] interval!') + return_dict = {} + power_acc = 0.0 + available_capacity_ABC = self.installed_battery_power_ABC + for batt_id in self.controllable_batteries_ABC.keys(): + measurements = self.controllable_batteries_ABC[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_ABC[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc < lower_limit: + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + else: + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + if available_capacity_ABC <= 0.0: + return return_dict, power_acc + for batt_id in self.controllable_batteries_ABC.keys(): + measurements = self.controllable_batteries_ABC[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + break + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_ABC[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc < lower_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_ABC[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_discharge_ABC / + available_capacity_ABC) * self.controllable_batteries_ABC[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_ABC[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1.0e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_ABC[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': new_power + } + power_acc += new_power + else: + power_acc += current_power + return return_dict, power_acc + + def calc_batt_charge_ABC(self, power_to_charge_ABC, upper_limit): + if not isinstance(power_to_charge_ABC, float): + raise TypeError('calc_batt_charge_ABC(): power_to_charge_ABC must be an instance of float!') + if power_to_charge_ABC < 0.0: + raise ValueError(f'calc_batt_charge_ABC(): power_to_charge_ABC must be nonnegative!') + if not isinstance(upper_limit, float): + raise TypeError('calc_batt_charge_ABC(): upper_limit must be an instance of float!') + if (upper_limit > 100.0) or (upper_limit < 0.0): + raise ValueError(f'calc_batt_charge_ABC(): upper_limit must belong to the [0, 1] interval!') + return_dict = {} + power_acc = 0.0 + available_capacity_ABC = self.installed_battery_power_ABC + for batt_id in self.controllable_batteries_ABC.keys(): + measurements = self.controllable_batteries_ABC[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_ABC[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc > upper_limit: + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + else: + available_capacity_ABC -= self.controllable_batteries_ABC[batt_id]['maximum_power'] + if available_capacity_ABC <= 0.0: + return return_dict, power_acc + for batt_id in self.controllable_batteries_ABC.keys(): + measurements = self.controllable_batteries_ABC[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + break + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_ABC[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc > upper_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_ABC[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_charge_ABC / + available_capacity_ABC) * self.controllable_batteries_ABC[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_ABC[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_ABC[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': -new_power + } + power_acc += new_power + else: + power_acc += current_power + return return_dict, power_acc + + def calc_batt_discharge_A(self, power_to_discharge_A, lower_limit): + if not isinstance(power_to_discharge_A, float): + raise TypeError('calc_batt_discharge_A(): power_to_discharge_A must be an instance of float!') + if power_to_discharge_A < 0.0: + raise ValueError(f'calc_batt_discharge_A(): power_to_discharge_A must be nonnegative!') + if not isinstance(lower_limit, float): + raise TypeError('calc_batt_discharge_A(): lower_limit must be an instance of float!') + if (lower_limit > 100.0) or (lower_limit < 0.0): + raise ValueError(f'calc_batt_discharge_A(): lower_limit must belong to the [0, 1] interval!') + return_dict = {} + available_capacity_A = self.installed_battery_power_A + for batt_id in self.controllable_batteries_A.keys(): + measurements = self.controllable_batteries_A[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_A[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc < lower_limit: + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + else: + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + if available_capacity_A <= 0.0: + return return_dict + for batt_id in self.controllable_batteries_A.keys(): + measurements = self.controllable_batteries_A[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + continue + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_A[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc < lower_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_A[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_discharge_A / + available_capacity_A) * self.controllable_batteries_A[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_A[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_A[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': new_power + } + return return_dict + + def calc_batt_charge_A(self, power_to_charge_A, upper_limit): + if not isinstance(power_to_charge_A, float): + raise TypeError('calc_batt_charge_A(): power_to_charge_A must be an instance of float!') + if power_to_charge_A < 0.0: + raise ValueError(f'calc_batt_charge_A(): power_to_charge_A must be nonnegative!') + if not isinstance(upper_limit, float): + raise TypeError('calc_batt_charge_A(): upper_limit must be an instance of float!') + if (upper_limit > 100.0) or (upper_limit < 0.0): + raise ValueError(f'calc_batt_charge_A(): upper_limit must belong to the [0, 1] interval!') + return_dict = {} + available_capacity_A = self.installed_battery_power_A + for batt_id in self.controllable_batteries_A.keys(): + measurements = self.controllable_batteries_A[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_A[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc > upper_limit: + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + else: + available_capacity_A -= self.controllable_batteries_A[batt_id]['maximum_power'] + if available_capacity_A <= 0.0: + return return_dict + for batt_id in self.controllable_batteries_A.keys(): + measurements = self.controllable_batteries_A[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + continue + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_A[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc > upper_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_A[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_charge_A / + available_capacity_A) * self.controllable_batteries_A[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_A[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_A[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': -new_power + } + return return_dict + + def calc_batt_discharge_B(self, power_to_discharge_B, lower_limit): + if not isinstance(power_to_discharge_B, float): + raise TypeError('calc_batt_discharge_B(): power_to_discharge_B must be an instance of float!') + if power_to_discharge_B < 0.0: + raise ValueError(f'calc_batt_discharge_B(): power_to_discharge_B must be nonnegative!') + if not isinstance(lower_limit, float): + raise TypeError('calc_batt_discharge_B(): lower_limit must be an instance of float!') + if (lower_limit > 100.0) or (lower_limit < 0.0): + raise ValueError(f'calc_batt_discharge_B(): lower_limit must belong to the [0, 1] interval!') + return_dict = {} + available_capacity_B = self.installed_battery_power_B + for batt_id in self.controllable_batteries_B.keys(): + measurements = self.controllable_batteries_B[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_B[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc < lower_limit: + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + else: + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + if available_capacity_B <= 0.0: + return return_dict + for batt_id in self.controllable_batteries_B.keys(): + measurements = self.controllable_batteries_B[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + continue + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_B[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc < lower_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_B[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_discharge_B / + available_capacity_B) * self.controllable_batteries_B[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_B[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_B[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': new_power + } + return return_dict + + def calc_batt_charge_B(self, power_to_charge_B, upper_limit): + if not isinstance(power_to_charge_B, float): + raise TypeError('calc_batt_charge_B(): power_to_charge_B must be an instance of float!') + if power_to_charge_B < 0.0: + raise ValueError(f'calc_batt_charge_B(): power_to_charge_B must be nonnegative!') + if not isinstance(upper_limit, float): + raise TypeError('calc_batt_charge_B(): upper_limit must be an instance of float!') + if (upper_limit > 100.0) or (upper_limit < 0.0): + raise ValueError(f'calc_batt_charge_B(): upper_limit must belong to the [0, 1] interval!') + return_dict = {} + available_capacity_B = self.installed_battery_power_B + for batt_id in self.controllable_batteries_B.keys(): + measurements = self.controllable_batteries_B[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_B[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc > upper_limit: + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + else: + available_capacity_B -= self.controllable_batteries_B[batt_id]['maximum_power'] + if available_capacity_B <= 0.0: + return return_dict + for batt_id in self.controllable_batteries_B.keys(): + measurements = self.controllable_batteries_B[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + continue + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_B[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc > upper_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_B[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_charge_B / + available_capacity_B) * self.controllable_batteries_B[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_B[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_B[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': -new_power + } + return return_dict + + def calc_batt_discharge_C(self, power_to_discharge_C, lower_limit): + if not isinstance(power_to_discharge_C, float): + raise TypeError('calc_batt_discharge_C(): power_to_discharge_C must be an instance of float!') + if power_to_discharge_C < 0.0: + raise ValueError(f'calc_batt_discharge_C(): power_to_discharge_C must be nonnegative!') + if not isinstance(lower_limit, float): + raise TypeError('calc_batt_discharge_C(): lower_limit must be an instance of float!') + if (lower_limit > 100.0) or (lower_limit < 0.0): + raise ValueError(f'calc_batt_discharge_C(): lower_limit must belong to the [0, 1] interval!') + return_dict = {} + available_capacity_C = self.installed_battery_power_C + for batt_id in self.controllable_batteries_C.keys(): + measurements = self.controllable_batteries_C[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_C[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc < lower_limit: + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + else: + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + if available_capacity_C <= 0.0: + return return_dict + for batt_id in self.controllable_batteries_C.keys(): + measurements = self.controllable_batteries_C[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + continue + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_C[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc < lower_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_C[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_discharge_C / + available_capacity_C) * self.controllable_batteries_C[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_C[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_C[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': new_power + } + return return_dict + + def calc_batt_charge_C(self, power_to_charge_C, upper_limit): + if not isinstance(power_to_charge_C, float): + raise TypeError('calc_batt_charge_C(): power_to_charge_C must be an instance of float!') + if power_to_charge_C < 0.0: + raise ValueError(f'calc_batt_charge_C(): power_to_charge_C must be nonnegative!') + if not isinstance(upper_limit, float): + raise TypeError('calc_batt_charge_C(): upper_limit must be an instance of float!') + if (upper_limit > 100.0) or (upper_limit < 0.0): + raise ValueError(f'calc_batt_charge_C(): upper_limit must belong to the [0, 1] interval!') + return_dict = {} + available_capacity_C = self.installed_battery_power_C + for batt_id in self.controllable_batteries_C.keys(): + measurements = self.controllable_batteries_C[batt_id]['power_measurements'] + mag = None + ang_in_deg = None + for measurement in measurements: + if measurement.get('value') is None: + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + break + mag = measurement.get('value', {}).get('magnitude') + ang_in_deg = measurement.get('value', {}).get('angle') + if (mag is None) or (ang_in_deg is None): + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + break + if (mag is None) or (ang_in_deg is None): + continue + batt_soc = self.controllable_batteries_C[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if batt_soc > upper_limit: + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + else: + available_capacity_C -= self.controllable_batteries_C[batt_id]['maximum_power'] + if available_capacity_C <= 0.0: + return return_dict + for batt_id in self.controllable_batteries_C.keys(): + measurements = self.controllable_batteries_C[batt_id]['power_measurements'] + mag = [] + ang_in_deg = [] + for measurement in measurements: + if measurement.get('value') is None: + continue + mag.append(measurement.get('value', {}).get('magnitude')) + ang_in_deg.append(measurement.get('value', {}).get('angle')) + if (not mag) or (not ang_in_deg) or (None in mag) or (None in ang_in_deg): + continue + current_power = 0.0 + for i in range(len(mag)): + current_power += mag[i] * math.cos(math.radians(ang_in_deg[i])) + batt_soc = self.controllable_batteries_C[batt_id]['soc_measurement'].get('value', {}).get('value') + if batt_soc is not None: + if (batt_soc > upper_limit) and (abs(current_power) > 1e-6): + return_dict[batt_id] = { + 'object': self.controllable_batteries_C[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': 0.0 + } + continue + else: + continue + new_power = (power_to_charge_C / + available_capacity_C) * self.controllable_batteries_C[batt_id]['maximum_power'] + new_power = min(new_power, self.controllable_batteries_C[batt_id]['maximum_power']) + if abs(new_power - current_power) > 1e-6: + return_dict[batt_id] = { + 'object': self.controllable_batteries_C[batt_id]['object'], + 'old_setpoint': current_power, + 'setpoint': -new_power + } + return return_dict + + def send_setpoints(self): + self.differenceBuilder.clear() + for dictVal in self.desired_setpoints.values(): + cimObj = dictVal.get('object') + newSetpoint = dictVal.get('setpoint') + oldSetpoint = dictVal.get('old_setpoint') + if isinstance(cimObj, cim.PowerElectronicsConnection): + self.differenceBuilder.add_difference(cimObj.PowerElectronicsUnit[0].mRID, + 'PowerElectronicsConnection.p', newSetpoint, oldSetpoint) + name = cimObj.Terminals[0].ConnectivityNode.name + data = [self.timestamp, cimObj.PowerElectronicsUnit[0].mRID, name, "ABC", newSetpoint] + add_data_to_csv(self.model_results_path/"dispatches.csv", data) + else: + self.log.warning( + f'The CIM object with mRID, {cimObj.mRID}, is not a cim.PowerElectronicsConnection. The ' + f'object is a {type(cimObj)}. This application will ignore sending a setpoint to this ' + 'object.') + setpointMessage = json.dumps(self.differenceBuilder.get_message(), indent=4, sort_keys=True) + self.gad_obj.send(self.setpoints_topic, setpointMessage) + + def simulation_completed(self, sim: Simulation): + self.log.info(f'Simulation for PeakShavingController:{self.id} has finished. This application ' + 'instance can be deleted.') + self.isValid = False + + +def main(control_enabled: bool, start_simulations: bool, model_id: str = None): + if not isinstance(control_enabled, bool): + raise TypeError(f'Argument, control_enabled, must be a boolean.') + if not isinstance(start_simulations, bool): + raise TypeError(f'Argument, start_simulations, must be a boolean.') + if not isinstance(model_id, str) and model_id is not None: + raise TypeError( + f'The model id passed to the convervation voltage reduction application must be a string type or {None}.') + global cim, DB_CONNECTION + cim_profile = 'rc4_2021' + iec61970_301 = 7 + cim = importlib.import_module(f'cimgraph.data_profile.{cim_profile}') + # params = ConnectionParameters(cim_profile=cim_profile, iec61970_301=iec61970_301) + # DB_CONNECTION = GridappsdConnection(params) + params = ConnectionParameters(url='http://localhost:8889/bigdata/namespace/kb/sparql', + cim_profile=cim_profile, + iec61970_301=iec61970_301) + DB_CONNECTION = BlazegraphConnection(params) + gad_object = createGadObject() + gad_log = gad_object.get_logger() + platform_simulations = {} + local_simulations = {} + app_instances = {'field_instances': {}, 'external_simulation_instances': {}, 'local_simulation_instances': {}} + response = gad_object.query_model_info() + models = response.get('data', {}).get('models', []) + model_is_valid = False + if model_id is not None: + for m in models: + m_id = m.get('modelId') + if model_id == m_id: + model_is_valid = True + break + if not model_is_valid: + raise ValueError(f'The model id provided does not exist in the GridAPPS-D plaform.') + else: + models = [m] + if start_simulations: + for m in models: + local_simulations[m.get('modelId', '')] = createSimulation(gad_object, m) + else: + # TODO: query platform for running simulations which is currently not implemented in the GridAPPS-D Api + pass + # Create an peak shaving controller instance for all the real systems in the database + # for m in models: + # m_id = m.get('modelId') + # app_instances['field_instances'][m_id] = PeakShavingController(gad_object, m_id) + for sim_id, m_id in platform_simulations.items(): + app_instances['external_simulation_instances'][sim_id] = PeakShavingController( + gad_object, + m_id, + sim_id=sim_id, + peak_setpoint = 0.35e6*3, + base_setpoint = 0.2e6*3, + ) + for m_id, simulation in local_simulations.items(): + app_instances['local_simulation_instances'][m_id] = PeakShavingController( + gad_object, + m_id, + peak_setpoint = 0.35e6*3, + base_setpoint = 0.2e6*3, + simulation=simulation + ) + app_instances_exist = False + if len(app_instances['field_instances']) > 0: + app_instances_exist = True + elif len(app_instances['external_simulation_instances']) > 0: + app_instances_exist = True + elif len(app_instances['local_simulation_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if app_instances_exist: + gad_log.info('PeakShavingApplication successfully started.') + gad_object.set_application_status(ProcessStatusEnum.RUNNING) + + while gad_object.connected: + try: + app_instances_exist = False + # TODO: check platform for any running external simulations to control. + invalidInstances = [] + for m_id, app in app_instances.get('field_instances', {}).items(): + if not app.isValid: + invalidInstances.append(m_id) + for m_id in invalidInstances: + del app_instances['field_instances'][m_id] + invalidInstances = [] + for sim_id, app in app_instances.get('external_simulation_instances', {}).items(): + # TODO: check if external simulation has finished and shutdown the corresponding app instance. + if not app.isValid: + invalidInstances.append(sim_id) + for sim_id in invalidInstances: + del app_instances['external_simulation_instances'][sim_id] + invalidInstances = [] + for m_id, app in app_instances.get('local_simulation_instances', {}).items(): + if not app.isValid: + invalidInstances.append(m_id) + for m_id in invalidInstances: + del app_instances['local_simulation_instances'][m_id] + if len(app_instances['field_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if len(app_instances['external_simulation_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if len(app_instances['local_simulation_instances']) > 0: + app_instances_exist = True + application_uptime = int(time.time()) + if not app_instances_exist: + application_downtime = int(time.time()) - application_uptime + if application_downtime > 3600: + gad_log.info('There have been no running instances PeakShavingController ' + 'instances for an hour. Shutting down the application.') + gad_object.set_application_status(ProcessStatusEnum.STOPPING) + gad_object.disconnect() + except KeyboardInterrupt: + gad_log.info('Manually exiting PeakShavingApplication') + gad_object.set_application_status(ProcessStatusEnum.STOPPING) + appList = list(app_instances.get('field_instances', {})) + for app_mrid in appList: + del app_instances['field_instances'][app_mrid] + appList = list(app_instances.get('external_simulation_instances', {})) + for app_mrid in appList: + del app_instances['external_simulation_instances'][app_mrid] + appList = list(app_instances.get('local_simulation_instances', {})) + for app_mrid in appList: + del app_instances['local_simulation_instances'][app_mrid] + gad_object.disconnect() + + +if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('model_id', + nargs='?', + default=None, + help='The mrid of the cim model to perform peak shaving on.') + parser.add_argument('-s', + '--start_simulations', + action='store_true', + help='Flag to have application start simulations') + parser.add_argument('-d', + '--disable_control', + action='store_true', + help='Flag to disable control on startup by default.') + args = parser.parse_args() + main(args.disable_control, args.start_simulations, args.model_id) diff --git a/pyproject.toml b/micro-apps/peak-shaving-app/pyproject.toml similarity index 64% rename from pyproject.toml rename to micro-apps/peak-shaving-app/pyproject.toml index 71acec6..40e3f17 100644 --- a/pyproject.toml +++ b/micro-apps/peak-shaving-app/pyproject.toml @@ -1,19 +1,23 @@ [tool.poetry] -name = "micro-apps" +name = "peak-shaving-app" version = "0.1.0" description = "" -authors = ["Your Name "] +authors = ["afisher1 <4552674+afisher1@users.noreply.github.com>"] readme = "README.md" [tool.poetry.dependencies] python = "^3.10" -gridappsd-python = "^2023.12.1" -cim-graph = "^0.1.2a0" +gridappsd-python = {version = "^2024", allow-prereleases = true} +cim-graph = "^0.1.5a0" [tool.poetry.group.dev.dependencies] -pre-commit = "^3.6.0" +pre-commit = "^3.7.0" yapf = "^0.40.2" -mypy = "^1.8.0" +mypy = "^1.9.0" + +[build-system] +requires = ["poetry-core"] +build-backend = "poetry.core.masonry.api" [tool.yapfignore] ignore_patterns = [ @@ -26,14 +30,5 @@ ignore_patterns = [ [tool.yapf] based_on_style = "pep8" spaces_before_comment = 4 -column_limit = 99 +column_limit = 120 split_before_logical_operator = true - -[tool.mypy] -show_error_context = true -pretty = true -show_column_numbers = true - -[build-system] -requires = ["poetry-core"] -build-backend = "poetry.core.masonry.api" diff --git a/micro-apps/power_factor/__init__.py b/micro-apps/power_factor/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/micro-apps/power_factor/main.py b/micro-apps/power_factor/main.py new file mode 100644 index 0000000..7f00176 --- /dev/null +++ b/micro-apps/power_factor/main.py @@ -0,0 +1,500 @@ +import os +import json +import traceback +import logging +from pprint import pprint +import networkx as nx +import numpy as np +import cvxpy as cp +import importlib +from datetime import datetime, timedelta, timezone +from dataclasses import dataclass, asdict +import cimgraph.utils +from cimgraph.databases import ConnectionParameters +from cimgraph.databases import BlazegraphConnection +from cimgraph.databases import GridappsdConnection +from cimgraph.models import FeederModel +import cimgraph.data_profile.rc4_2021 as cim_dp +from gridappsd import GridAPPSD, DifferenceBuilder +from gridappsd.simulation import Simulation +from gridappsd.simulation import PowerSystemConfig +from gridappsd.simulation import ApplicationConfig +from gridappsd.simulation import SimulationConfig +from gridappsd.simulation import SimulationArgs +from gridappsd.simulation import ModelCreationConfig +from gridappsd.simulation import ServiceConfig +from gridappsd import topics as t +import models +import query + +IEEE123_APPS = "_E3D03A27-B988-4D79-BFAB-F9D37FB289F7" +IEEE13 = "_49AD8E07-3BF9-A4E2-CB8F-C3722F837B62" +IEEE123 = "_C1C3E687-6FFD-C753-582B-632A27E28507" +IEEE123PV = "_E407CBB6-8C8D-9BC9-589C-AB83FBF0826D" +SOURCE_BUS = "150r" +OUT_DIR = "outputs" +ROOT = os.getcwd() + + +logging.basicConfig(level=logging.DEBUG, + filename=f"{ROOT}/{OUT_DIR}/log.txt", filemode='w') +log = logging.getLogger(__name__) + + +def remap_phases( + info: models.MeasurementInfo, + map: models.PhaseMap) -> models.MeasurementInfo: + if info.phase == cim_dp.PhaseCode.s1: + info.phase = map.s1 + return info + + if info.phase == cim_dp.PhaseCode.s2: + info.phase = map.s2 + return info + + return info + + +def dist_matrix(ders: dict, dists: dict) -> np.ndarray: + phases = 3 + map = [1/(1+dists[bus]) for der, bus in ders.items()] + map = map / np.linalg.norm(map) + mat = np.array([map,]*phases).transpose() + return mat + + +def dict_to_array(d: dict) -> np.ndarray: + results = d.values() + data = list(results) + return np.array(data) + + +def save_data(results: models.DataInfo, filename: str) -> None: + with open(f"{OUT_DIR}/{filename}.json", "w") as f: + f.write(json.dumps(asdict(results))) + + +class PowerFactor(object): + gapps: GridAPPSD + sim: Simulation + network: cimgraph.GraphModel + compensators: models.Compensators + consumers: models.Consumers + electronics: models.PowerElectronics + source_mrid: str + data_info: models.DataInfo + + def __init__(self, gapps: GridAPPSD, sim: Simulation, network: cimgraph.GraphModel): + self.simulation = sim + self.gapps = gapps + self.switches = query.get_switches(network) + self.compensators = query.get_compensators(network) + self.consumers = query.get_consumers(network) + self.electronics = query.get_power_electronics(network) + self.tranformers = query.get_Transformers(network) + self.data_info = models.DataInfo() + + self.data_info.compensators_ratings = self.compensators.ratings + self.data_info.pecs_ratings = self.electronics.ratings + self.temp_source_bus = ["_9d721fcc-7229-42f5-b2e2-d9b0a8adecfd", + "_81930576-6f0d-4c4f-a857-8e359b5210cd", "_425640dd-811d-483f-9675-03d9b516b7f8"] + + (graph, source_bus, mrid) = query.generate_graph(network) + self.source_mrid = mrid + paths = nx.shortest_path_length(graph, source=source_bus) + pec_map = query.map_power_electronics(network) + self.dist = dist_matrix(pec_map, paths) + self.data_info.pecs_distance = {k: (1/v[0]) - 1 + for k, v in zip(pec_map.keys(), self.dist)} + + sim.add_onmeasurement_callback(self.on_measurement) + self.diff_builder = DifferenceBuilder(sim.simulation_id) + self.topic = t.simulation_input_topic(sim.simulation_id) + self.last_dispatch = 0 + + # self.init_electronics() + + def on_measurement(self, sim: Simulation, timestamp: dict, measurements: dict) -> None: + print(f"{timestamp}: on_measurement") + if timestamp - self.last_dispatch >= 6: + + pecs = {key: models.PhasePower() + for key in self.electronics.ratings.keys()} + + self.data_info.data.append(models.Data( + timestamp=timestamp, + total_load=models.PhasePower(), + net_load=models.PhasePower(), + pec_dispatch=models.PhasePower(), + pecs=pecs + )) + self.last_dispatch = timestamp + control = self.dispatch() + print("Total Reactive PECS: ", np.sum(control, axis=0)) + self.set_electronics(control) + + for k, v in measurements.items(): + if k in self.temp_source_bus: + pass + # print("SOURCE: ", k, v) + + if k in self.switches.measurement_map: + val = models.SimulationMeasurement(**v) + if val.value is not None: + info = self.switches.measurement_map[k] + + if k in self.consumers.measurement_map: + val = models.SimulationMeasurement(**v) + if val.value is not None: + continue + + info = self.consumers.measurement_map[k] + if info.value_type == "PNV": + old = self.consumers.measurements_pnv[info.mrid] + new = models.update_pnv(info, val, old) + self.consumers.measurements_pnv[info.mrid] = new + + if info.value_type == "VA": + map = self.consumers.measurements_pnv[info.mrid] + info = remap_phases(info, map) + self.consumers.measurement_map[k] = info + + old = self.consumers.measurements_va[info.mrid] + new = models.update_va(info, val, old) + self.consumers.measurements_va[info.mrid] = new + + if k in self.compensators.measurement_map: + val = models.SimulationMeasurement(**v) + if val.value is not None: + continue + + info = self.compensators.measurement_map[k] + if info.value_type == "PNV": + old = self.compensators.measurements_pnv[info.mrid] + new = models.update_pnv(info, val, old) + self.compensators.measurements_pnv[info.mrid] = new + + if info.value_type == "VA": + map = self.compensators.measurements_pnv[info.mrid] + info = remap_phases(info, map) + self.compensators.measurement_map[k] = info + + old = self.compensators.measurements_va[info.mrid] + new = models.update_va(info, val, old) + self.compensators.measurements_va[info.mrid] = new + + if k in self.tranformers.measurement_map: + val = models.SimulationMeasurement(**v) + if val.value is not None: + continue + + info = self.tranformers.measurement_map[k] + if info.value_type == "PNV": + old = self.tranformers.measurements_pnv[info.mrid] + new = models.update_pnv(info, val, old) + self.tranformers.measurements_pnv[info.mrid] = new + + if info.value_type == "VA": + map = self.tranformers.measurements_pnv[info.mrid] + info = remap_phases(info, map) + self.tranformers.measurement_map[k] = info + + old = self.tranformers.measurements_va[info.mrid] + new = models.update_va(info, val, old) + self.tranformers.measurements_va[info.mrid] = new + + if k in self.electronics.measurement_map: + val = models.SimulationMeasurement(**v) + if val.value is not None: + continue + + info = self.electronics.measurement_map[k] + if info.value_type == "PNV": + old = self.electronics.measurements_pnv[info.mrid] + new = models.update_pnv(info, val, old) + self.electronics.measurements_pnv[info.mrid] = new + + if info.value_type == "VA": + map = self.electronics.measurements_pnv[info.mrid] + info = remap_phases(info, map) + self.electronics.measurement_map[k] = info + + old = self.electronics.measurements_va[info.mrid] + new = models.update_va(info, val, old) + self.electronics.measurements_va[info.mrid] = new + if self.simulation is not None: + self.simulation.resume() + + def toggle_switches(self) -> None: + v: models.MeasurementInfo + for k, v in self.switches.measurement_map.items(): + if v.phase == "A": + self.diff_builder.add_difference( + k, models.Switches.open_attribute, True, False + ) + + topic = t.simulation_input_topic(sim.simulation_id) + message = self.diff_builder.get_message() + # self.gapps.send(topic, message) + + def set_compensators(self) -> np.ndarray: + real_loads = dict_to_array(self.consumers.measurements_va)[:, :, 0] + total_real = np.sum(real_loads, axis=0).tolist() + + imag_loads = dict_to_array(self.consumers.measurements_va)[:, :, 1] + total_imag = np.sum(imag_loads, axis=0).tolist() + + starting_load = np.sum(imag_loads, axis=0) + total_load = np.sum(imag_loads, axis=0) + + self.data_info.data[-1].total_load.set_phases(total_real, total_imag) + + comps_abc = {k: v for k, v in self.compensators.ratings.items() + if v.a.imag != 0 and v.b.imag != 0 and v.c.imag != 0} + comps_a = {k: v for k, v in self.compensators.ratings.items() + if v.b.imag == 0 and v.c.imag == 0} + comps_b = {k: v for k, v in self.compensators.ratings.items() + if v.a.imag == 0 and v.c.imag == 0} + comps_c = {k: v for k, v in self.compensators.ratings.items() + if v.a.imag == 0 and v.b.imag == 0} + + for k, v in comps_abc.items(): + vars = np.sum(np.array(v), axis=1) + net = total_load + vars + if net.all() < 0: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 0, 1 + ) + else: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 1, 0 + ) + total_load = net + + for k, v in comps_a.items(): + vars = np.sum(np.array(v), axis=1) + net = total_load + vars + if net[0] < 0: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 0, 1 + ) + else: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 1, 0 + ) + total_load = net + + for k, v in comps_b.items(): + vars = np.sum(np.array(v), axis=1) + net = total_load + vars + if net[1] < 0: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 0, 1 + ) + else: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 1, 0 + ) + total_load = net + + for k, v in comps_c.items(): + vars = np.sum(np.array(v), axis=1) + net = total_load + vars + if net[2] < 0: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 0, 1 + ) + else: + self.diff_builder.add_difference( + k, models.Compensators.sections_attribute, 1, 0 + ) + total_load = net + + topic = t.simulation_input_topic(sim.simulation_id) + message = self.diff_builder.get_message() + # self.gapps.send(topic, message) + net = starting_load - total_load + self.data_info.data[-1].net_load.set_phases( + total_real, net.tolist()) + return total_load + + def init_electronics(self) -> None: + for idx, key in enumerate(self.electronics.ratings.keys()): + mrid = self.electronics.units[key] + self.diff_builder.add_difference( + mrid, models.PowerElectronics.p_attribute, 0.0, 0.0 + ) + topic = t.simulation_input_topic(sim.simulation_id) + message = self.diff_builder.get_message() + # self.gapps.send(topic, message) + + def set_electronics(self, dispatch: np.ndarray) -> None: + real = dict_to_array(self.electronics.measurements_va)[:, :, 0] + total_real = np.sum(real, axis=0) + total_imag = np.sum(dispatch, axis=0) + self.data_info.data[-1].pec_dispatch.set_phases(total_real, total_imag) + + for idx, key in enumerate(self.electronics.ratings.keys()): + mrid = self.electronics.units[key] + self.data_info.data[-1].pecs[key].set_phases( + real[idx], dispatch[idx].tolist()) + value = np.sum(dispatch[idx]).item() + self.diff_builder.add_difference( + mrid, models.PowerElectronics.q_attribute, value, 0.0 + ) + topic = t.simulation_input_topic(sim.simulation_id) + message = self.diff_builder.get_message() + # self.gapps.send(topic, message) + + def get_phase_vars(self) -> np.ndarray: + load = dict_to_array(self.consumers.measurements_va)[:, :, 1] + total_load = np.sum(load, axis=0) + + comp = dict_to_array(self.compensators.measurements_va)[:, :, 1] + total_comp = np.sum(comp, axis=0) + + print("Total consumers: ", total_load) + print("Total compensators: ", total_comp) + + return total_load + total_comp + + def get_bounds(self) -> np.ndarray: + apparent = dict_to_array(self.electronics.ratings)[:, :, 0] + real = dict_to_array(self.electronics.measurements_va)[:, :, 0] + total_real = np.sum(real, axis=0).tolist() + print("Total Apparent PECS: ", np.sum(apparent, axis=0)) + print("Total Real PECS: ", np.sum(real, axis=0)) + reactive = np.sqrt(abs(real**2 - apparent**2)) + return reactive + + def dispatch(self) -> np.ndarray: + bounds = self.get_bounds() + A = np.clip(bounds, a_min=0, a_max=1) + b = self.set_compensators() + + # constrain values to lesser of ders and network + total_der = np.sum(bounds, axis=0) + direction = np.clip(b, a_max=1, a_min=-1) + b = np.array([min(vals) for vals in zip(abs(total_der), abs(b))]) + b = b*direction + + # Construct the problem. + m, n = np.shape(A) + x = cp.Variable((m, n)) + + cost = cp.sum(self.dist.T@cp.abs(x)) + objective = cp.Minimize(cost) + constraints = [ + cp.sum(A.T@x, axis=0) == b, + -bounds <= x, + x <= bounds] + + prob = cp.Problem(objective, constraints) + + # The optimal objective is returned by prob.solve(). + prob.solve(solver=cp.CLARABEL, verbose=False) + + if prob.status == 'optimal': + return np.round(x.value, 0) + + return np.zeros_like(x) + + +@dataclass +class ModelInfo: + modelName: str + modelId: str + stationName: str + stationId: str + subRegionName: str + subRegionId: str + regionName: str + regionId: str + + +if __name__ == "__main__": + try: + cim_profile = 'rc4_2021' + cim = importlib.import_module('cimgraph.data_profile.' + cim_profile) + mrid = IEEE123PV + + feeder = cim.Feeder(mRID=mrid) + + params = ConnectionParameters( + url="http://localhost:8889/bigdata/namespace/kb/sparql", + cim_profile=cim_profile) + bgc = BlazegraphConnection(params) + + # params = ConnectionParameters( + # cim_profile=cim_profile, + # iec61970_301=7 + # ) + # gac = GridappsdConnection(params) + gapps = GridAPPSD(username='system', password='manager') + network = FeederModel( + connection=bgc, + container=feeder, + distributed=False) + + cimgraph.utils.get_all_data(network) + cimgraph.utils.get_all_measurement_data(network) + + model_info = gapps.query_model_info()['data']['models'] + for model in model_info: + if model['modelId'] == mrid: + system = ModelInfo(**model) + + system_config = PowerSystemConfig( + GeographicalRegion_name=system.regionId, + SubGeographicalRegion_name=system.subRegionId, + Line_name=system.modelId + ) + + model_config = ModelCreationConfig( + load_scaling_factor=1, + schedule_name="ieeezipload", + z_fraction=0, + i_fraction=1, + p_fraction=0, + randomize_zipload_fractions=False, + use_houses=False + ) + + start = datetime(2023, 1, 1, 4) + epoch = datetime.timestamp(start) + duration = timedelta(hours=24).total_seconds() + + start = datetime(2023, 1, 1, 12, tzinfo=timezone.utc) + epoch = datetime.timestamp(start) + duration = timedelta(hours=1).total_seconds() + sim_args = SimulationArgs( + start_time=epoch, + duration=duration, + simulator="GridLAB-D", + timestep_frequency=1000, + timestep_increment=1000, + run_realtime=False, + simulation_name=system.modelName, + power_flow_solver_method="NR", + model_creation_config=asdict(model_config), + pause_after_measurements=True, + ) + + sim_config = SimulationConfig( + power_system_config=asdict(system_config), + simulation_config=asdict(sim_args) + ) + sim = Simulation(gapps, asdict(sim_config)) + + app = PowerFactor(gapps, sim, network) + + sim.run_loop() + + save_data(app.data_info, "summary") + + except KeyboardInterrupt: + sim.stop() + + except Exception as e: + log.debug(e) + log.debug(traceback.format_exc()) diff --git a/micro-apps/power_factor/models.py b/micro-apps/power_factor/models.py new file mode 100644 index 0000000..63c4535 --- /dev/null +++ b/micro-apps/power_factor/models.py @@ -0,0 +1,234 @@ +from dataclasses import dataclass, field, astuple +import cimgraph.data_profile.rc4_2021 as cim_dp +import numpy as np + + +@dataclass +class PhaseSwitch: + a: bool = False + b: bool = False + c: bool = False + + def __array__(self): + return np.array(astuple(self)) + + def __len__(self): + return astuple(self).__len__() + + def __getitem__(self, item): + return astuple(self).__getitem__(item) + + +@dataclass +class ComplexPower: + real: float = 0.0 + imag: float = 0.0 + + def __array__(self): + return np.array(astuple(self)) + + def __len__(self): + return astuple(self).__len__() + + def __getitem__(self, item): + return astuple(self).__getitem__(item) + + +@dataclass +class PhasePower: + a: ComplexPower = field(default=ComplexPower()) + b: ComplexPower = field(default=ComplexPower()) + c: ComplexPower = field(default=ComplexPower()) + + def __array__(self): + return np.array(astuple(self)) + + def __len__(self): + return astuple(self).__len__() + + def __getitem__(self, item): + return astuple(self).__getitem__(item) + + def set_phases(self, real: list[float], imag: list[float]) -> None: + if len(real) == 0: + real = [self.a.real, self.b.real, self.c.real] + + if len(imag) == 0: + imag = [self.a.imag, self.b.imag, self.c.imag] + + assert (3 == len(real)) + assert (3 == len(imag)) + print(real, imag) + self.a = ComplexPower( + real=real[0], imag=imag[0]) + self.b = ComplexPower( + real=real[1], imag=imag[1]) + self.c = ComplexPower( + real=real[2], imag=imag[2]) + + +@dataclass +class PhaseMap: + s1: cim_dp.PhaseCode | None = None + s2: cim_dp.PhaseCode | None = None + + +@dataclass +class SimulationMeasurement: + measurement_mrid: str + value: int | None = None + angle: float | None = None + magnitude: float | None = None + + +@dataclass +class MeasurementInfo: + mrid: str + value_type: str + phase: cim_dp.PhaseCode + + +def update_va( + info: MeasurementInfo, + val: SimulationMeasurement, + power: PhasePower) -> PhasePower: + if val.angle is None or val.magnitude is None: + return power + + mag = float(val.magnitude) + ang = float(val.angle) + val = ComplexPower( + real=mag * np.cos(np.deg2rad(ang)), + imag=mag * np.sin(np.deg2rad(ang)) + ) + + if info.phase == cim_dp.PhaseCode.A: + power.a = val + return power + + if info.phase == cim_dp.PhaseCode.B: + power.b = val + return power + + if info.phase == cim_dp.PhaseCode.C: + power.c = val + return power + + return power + + +def find_nearest(array: np.array, value: float) -> int: + return (np.abs(array - value)).argmin() + + +def update_pnv( + info: MeasurementInfo, + val: SimulationMeasurement, + phase_map: PhaseMap) -> PhaseMap: + if val.angle is None or val.magnitude is None: + return phase_map + + ang = float(val.angle) + phases = np.array([0.0, -120.0, 120.0]) + phase_id = find_nearest(phases, ang) + + if info.phase == cim_dp.PhaseCode.s1: + if phase_id == 0: + phase_map.s1 = cim_dp.PhaseCode.A + return phase_map + + if phase_id == 1: + phase_map.s1 = cim_dp.PhaseCode.B + return phase_map + + if phase_id == 2: + phase_map.s1 = cim_dp.PhaseCode.C + return phase_map + + if info.phase == cim_dp.PhaseCode.s2: + if phase_id == 0: + phase_map.s2 = cim_dp.PhaseCode.A + return phase_map + + if phase_id == 1: + phase_map.s2 = cim_dp.PhaseCode.B + return phase_map + + if phase_id == 2: + phase_map.s2 = cim_dp.PhaseCode.C + return phase_map + + return phase_map + + +@dataclass +class Compensators: + ratings: dict[PhasePower] = field(default_factory=dict) + measurements_va: dict[PhasePower] = field(default_factory=dict) + measurements_pnv: dict[PhaseMap] = field(default_factory=dict) + measurement_map: dict[MeasurementInfo] = field(default_factory=dict) + sections_attribute: str = "ShuntCompensator.sections" + + +@dataclass +class PowerElectronics: + units: dict[str] = field(default_factory=dict) + ratings: dict[PhasePower] = field(default_factory=dict) + measurements_va: dict[PhasePower] = field(default_factory=dict) + measurements_pnv: dict[PhaseMap] = field(default_factory=dict) + measurement_map: dict[MeasurementInfo] = field(default_factory=dict) + p_attribute: str = "PowerElectronicsConnection.p" + q_attribute: str = "PowerElectronicsConnection.q" + + +@dataclass +class Generators: + units: dict[str] = field(default_factory=dict) + ratings: dict[PhasePower] = field(default_factory=dict) + measurements_va: dict[PhasePower] = field(default_factory=dict) + measurements_pnv: dict[PhaseMap] = field(default_factory=dict) + measurement_map: dict[MeasurementInfo] = field(default_factory=dict) + p_attribute: str = "RotatingMachine.p" + q_attribute: str = "RotatingMachine.q" + + +@dataclass +class Switches: + normal: dict[PhaseSwitch] = field(default_factory=dict) + measurements_va: dict[PhaseSwitch] = field(default_factory=dict) + measurements_pnv: dict[PhaseMap] = field(default_factory=dict) + measurement_map: dict[MeasurementInfo] = field(default_factory=dict) + open_attribute: str = "Switch.open" + + +@dataclass +class Transformers: + measurements_va: dict[PhaseSwitch] = field(default_factory=dict) + measurements_pnv: dict[PhaseMap] = field(default_factory=dict) + measurement_map: dict[MeasurementInfo] = field(default_factory=dict) + open_attribute: str = "TapChanger.step" + + +@dataclass +class Consumers: + ratings: dict[PhasePower] = field(default_factory=dict) + measurements_va: dict[PhasePower] = field(default_factory=dict) + measurements_pnv: dict[PhaseMap] = field(default_factory=dict) + measurement_map: dict[MeasurementInfo] = field(default_factory=dict) + + +@dataclass +class Data: + timestamp: int + total_load: PhasePower + net_load: PhasePower + pec_dispatch: PhasePower + pecs: dict[PhasePower] = field(default_factory=dict) + + +@dataclass +class DataInfo: + pecs_distance: dict[float] = field(default_factory=dict) + pecs_ratings: dict[PhasePower] = field(default_factory=dict) + compensators_ratings: dict[PhasePower] = field(default_factory=dict) + data: list[Data] = field(default_factory=list) diff --git a/micro-apps/power_factor/plot.py b/micro-apps/power_factor/plot.py new file mode 100644 index 0000000..29ae565 --- /dev/null +++ b/micro-apps/power_factor/plot.py @@ -0,0 +1,269 @@ +import os +import json +import math +import traceback +import networkx as nx +import numpy as np +# from pandas import DataFrame +import matplotlib.pyplot as plt +from matplotlib.collections import PolyCollection +from datetime import timedelta +from models import DataInfo, PhasePower, Data + +ROOT = os.getcwd() +IN_DIR = f"{ROOT}/inputs" +OUT_DIR = f"{ROOT}/outputs" + + +def load_summary(filepath: str) -> DataInfo: + with open(filepath) as f: + data = json.load(f) + return DataInfo(**data) + + +def extract_real(data: list[PhasePower]) -> (list[float], list[float], list[float]): + a = [var['a']['real']/1000 for var in data] + b = [var['b']['real']/1000 for var in data] + c = [var['c']['real']/1000 for var in data] + return (a, b, c) + + +def extract_imag(data: list[PhasePower]) -> (list[float], list[float], list[float]): + a = [var['a']['imag']/1000 for var in data] + b = [var['b']['imag']/1000 for var in data] + c = [var['c']['imag']/1000 for var in data] + return (a, b, c) + + +def extract_comps_real(data: dict[PhasePower]) -> (dict[float], dict[float], dict[float]): + a = {} + b = {} + c = {} + for k, v in data.items(): + if v['a']['real'] != 0.0: + a[k] = v['a']['real']/1000 + if v['b']['real'] != 0.0: + b[k] = v['b']['real']/1000 + if v['c']['real'] != 0.0: + c[k] = v['c']['real']/1000 + return (a, b, c) + + +def extract_comps_imag(data: dict[PhasePower]) -> (dict[float], dict[float], dict[float]): + a = {} + b = {} + c = {} + for k, v in data.items(): + if v['a']['imag'] != 0.0: + a[k] = v['a']['imag']/1000 + if v['b']['imag'] != 0.0: + b[k] = v['b']['imag']/1000 + if v['c']['imag'] != 0.0: + c[k] = v['c']['imag']/1000 + return (a, b, c) + + +def extract_pecs_real(data: dict[PhasePower]) -> (dict[float], dict[float], dict[float]): + a = {} + b = {} + c = {} + for k, v in data.items(): + if v['a']['real'] != 0.0: + a[k] = v['a']['real']/1000 + if v['b']['real'] != 0.0: + b[k] = v['b']['real']/1000 + if v['c']['real'] != 0.0: + c[k] = v['c']['real']/1000 + return (a, b, c) + + +def extract_pecs_imag(data: dict[PhasePower]) -> (dict[float], dict[float], dict[float]): + a = {} + b = {} + c = {} + for k, v in data.items(): + if v['a']['imag'] != 0.0: + a[k] = v['a']['imag']/1000 + if v['b']['imag'] != 0.0: + b[k] = v['b']['imag']/1000 + if v['c']['imag'] != 0.0: + c[k] = v['c']['imag']/1000 + return (a, b, c) + + +def percent_avail(s: float, p: float, q: float) -> float: + if q == 0.0: + return 0 + + num = abs(q) + den = math.sqrt(s**2 - p**2) + return num/den + + +def gather_pecs( + dist: dict[float], + nameplate: dict[float], + real: dict[float], + imag: dict[float]) -> (list[float], list[float]): + x = [] + y = [] + for k, d in dist.items(): + if k in nameplate: + va = nameplate[k] + + w = 0.0 + if k in real: + w = real[k] + + var = 0.0 + if k in imag: + var = imag[k] + + y.append(percent_avail(va, w, var)) + x.append(d) + return (x, y) + + +def plot_pecs(summary: DataInfo, idx: int) -> None: + + dist = data.pecs_distance + (np_a, np_b, np_c) = extract_pecs_real(data.pecs_ratings) + pecs = data.data[idx]['pecs'] + (imag_a, imag_b, imag_c) = extract_pecs_imag(pecs) + (real_a, real_b, real_c) = extract_pecs_real(pecs) + + ax, ay = gather_pecs(dist, np_a, real_a, imag_a) + bx, by = gather_pecs(dist, np_b, real_b, imag_b) + cx, cy = gather_pecs(dist, np_c, real_c, imag_c) + + x = ax + bx + cx + y = ay + by + cy + + fig, axis = plt.subplots() + axis.scatter(x, y) + + axis.set(xlabel='Distance (km)', ylabel='Reactive Power (% avail)') + plt.savefig(f'{OUT_DIR}/dist_var_{idx}.png', dpi=400) + + +def plot(summary: DataInfo) -> None: + + total_load = [data['total_load'] for data in summary.data[2:]] + (total_a, total_b, total_c) = extract_imag(total_load) + + net_load = [data['net_load'] for data in summary.data[2:]] + (net_a, net_b, net_c) = extract_imag(net_load) + + pecs = [data['pec_dispatch'] for data in summary.data[2:]] + (pecs_a, pecs_b, pecs_c) = extract_imag(pecs) + + fig, ax = plt.subplots(3, sharex=True) + ax[0].plot(range(len(total_a)), total_a, label='total') + ax[0].plot(range(len(net_a)), net_a, label='net') + ax[0].plot(range(len(pecs_a)), pecs_a, label='pecs') + ax[1].plot(range(len(total_b)), total_b, label='total') + ax[1].plot(range(len(net_b)), net_b, label='net') + ax[1].plot(range(len(pecs_b)), pecs_b, label='pecs') + ax[2].plot(range(len(total_c)), total_c, label='total') + ax[2].plot(range(len(net_c)), net_c, label='net') + ax[2].plot(range(len(pecs_c)), pecs_c, label='pecs') + + phases = ["A", "B", "C"] + for idx, axis in enumerate(ax.flat): + axis.set(xlabel='Time (H)', ylabel=f"{phases[idx]} (kVAR)") + axis.label_outer() + + fig.legend(["Total", "Compensators", "PECs"], loc='upper center', ncols=3) + plt.savefig(f'{OUT_DIR}/total_var.png', dpi=400) + + +def plot_graph(G: nx.Graph, dist: dict, pos: dict) -> None: + + # nodes + nx.draw_networkx_nodes( + G, + pos, + node_size=20, + nodelist=list(dist.keys()), + node_color=list(dist.values()), + cmap=plt.cm.plasma + ) + + # node labels + nx.draw_networkx_labels(G, pos, font_size=2, font_family="sans-serif") + + # edges + nx.draw_networkx_edges(G, pos, width=2, alpha=0.4) + + # edge weight labels + edge_labels = nx.get_edge_attributes(G, "weight") + nx.draw_networkx_edge_labels(G, pos, edge_labels, font_size=2) + + ax = plt.gca() + ax.margins(0.08) + plt.axis("off") + plt.tight_layout() + plt.savefig("outputs/graph.png", dpi=400) + + +def polygon_under_graph(x: list[float], y: list[float]): + """ + Construct the vertex list which defines the polygon filling the space under + the (x, y) line graph. This assumes x is in ascending order. + """ + return [(x[0], 0.), *zip(x, y), (x[-1], 0.)] + + +def plot_3d(summary: DataInfo) -> None: + dist = data.pecs_distance + (np_a, np_b, np_c) = extract_pecs_real(data.pecs_ratings) + + axis = plt.figure().add_subplot(projection='3d') + + x = [] + y = [] + z = [] + for step in data.data: + ts = step['timestamp'] + if ts % 3600 != 0: + continue + + pecs = step['pecs'] + (imag_a, imag_b, imag_c) = extract_pecs_imag(pecs) + (real_a, real_b, real_c) = extract_pecs_real(pecs) + ax, ay = gather_pecs(dist, np_a, real_a, imag_a) + bx, by = gather_pecs(dist, np_b, real_b, imag_b) + cx, cy = gather_pecs(dist, np_c, real_c, imag_c) + combined = list(zip(ax+bx+cx, ay+by+cy)) + sorted_combined = sorted(combined) + d, q = zip(*sorted_combined) + x.append(d) + y.append(q) + z.append(timedelta(seconds=ts).seconds//3600) + + verts = [polygon_under_graph(a, b) for a, b in zip(x, y)] + facecolors = plt.colormaps['viridis_r'](np.linspace(0, 1, len(verts))) + poly = PolyCollection(verts, facecolors=facecolors, alpha=.25) + axis.add_collection3d(poly, zs=z, zdir='y') + for idx, hour in enumerate(z): + z = [hour]*len(x[idx]) + axis.plot(x[idx], hour, y[idx], c=facecolors[idx], alpha=1) + + axis.set( + xlabel='Dist (km)', + zlabel='VAR (%)', + ylabel='Time (HH)') + axis.invert_xaxis() + + plt.savefig(f'{OUT_DIR}/pecs_3d.png', dpi=400) + + +if __name__ == "__main__": + try: + data = load_summary(f"{OUT_DIR}/summary.json") + plot(data) + plot_3d(data) + + except Exception as e: + print(e) + print(traceback.format_exc()) diff --git a/micro-apps/power_factor/query.py b/micro-apps/power_factor/query.py new file mode 100644 index 0000000..212b016 --- /dev/null +++ b/micro-apps/power_factor/query.py @@ -0,0 +1,471 @@ +import networkx as nx +import cimgraph.utils +import cimgraph.models +import cimgraph.data_profile.rc4_2021 as cim_dp +import models +import logging + +log = logging.getLogger(__name__) + + +def generate_graph(network: cimgraph.models.GraphModel) -> (nx.Graph, str, str): + graph = nx.Graph() + + network.get_all_edges(cim_dp.ACLineSegment) + network.get_all_edges(cim_dp.EnergySource) + network.get_all_edges(cim_dp.ConnectivityNode) + network.get_all_edges(cim_dp.Terminal) + + source_bus = "" + next_bus = "" + mrid = "" + source: cim_dp.EnergySource + for source in network.graph[cim_dp.EnergySource].values(): + terminals = source.Terminals + bus_1 = terminals[0].ConnectivityNode.name + source_bus = bus_1 + + xfmr: cim_dp.PowerTransformer + for xfmr in network.graph[cim_dp.PowerTransformer].values(): + ends = xfmr.Terminals + buses = [] + for end in ends: + buses.append(end.ConnectivityNode.name) + + buses = list(set(buses)) + bus_1 = buses[0] + bus_2 = buses[-1] + if source_bus == bus_1: + mrid = xfmr.mRID + next_bus = bus_2 + if source_bus == bus_2: + mrid = xfmr.mRID + next_bus = bus_1 + graph.add_edge(bus_1, bus_2, weight=0.0) + + line: cim_dp.ACLineSegment + for line in network.graph[cim_dp.ACLineSegment].values(): + terminals = line.Terminals + bus_1 = terminals[0].ConnectivityNode.name + bus_2 = terminals[1].ConnectivityNode.name + if next_bus == bus_1: + mrid = line.mRID + if next_bus == bus_2: + mrid = line.mRID + graph.add_edge(bus_1, bus_2, weight=float(line.length)) + + switch: cim_dp.LoadBreakSwitch + for switch in network.graph[cim_dp.LoadBreakSwitch].values(): + # if not switch.normalOpen: + terminals = switch.Terminals + bus_1 = terminals[0].ConnectivityNode.name + bus_2 = terminals[1].ConnectivityNode.name + if source_bus == bus_1: + mrid = switch.mRID + if source_bus == bus_2: + mrid = switch.mRID + graph.add_edge(bus_1, bus_2, weight=0.0) + + return (graph, source_bus, mrid) + + +def get_graph_positions(network: cimgraph.models.GraphModel) -> dict: + loc = cimgraph.utils.get_all_bus_locations(network) + return {d['name']: [float(d['x']), float(d['y'])] + for d in loc.values()} + + +def get_compensators(network: cimgraph.models.GraphModel) -> models.Compensators: + compensators = models.Compensators() + if cim_dp.LinearShuntCompensator not in network.graph: + return compensators + + shunt: cim_dp.LinearShuntCompensator + for shunt in network.graph[cim_dp.LinearShuntCompensator].values(): + + mrid = shunt.mRID + nom_v = float(shunt.nomU) + p_imag = -float(shunt.bPerSection) * nom_v**2 + + if not shunt.ShuntCompensatorPhase: + measurement = cim_dp.Measurement + for measurement in shunt.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + compensators.measurement_map[measurement.mRID] = info + + compensators.measurements_va[mrid] = models.PhasePower() + compensators.measurements_pnv[mrid] = models.PhaseMap() + + power = models.ComplexPower(real=0.0, imag=p_imag/3) + phases = models.PhasePower(a=power, b=power, c=power) + compensators.ratings[mrid] = phases + + else: + ratings = models.PhasePower() + phase: models.LinearShuntCompensatorPhase + for phase in shunt.ShuntCompensatorPhase: + power = models.ComplexPower(real=0.0, imag=p_imag) + if phase.phase == cim_dp.SinglePhaseKind.A: + ratings.a = power + if phase.phase == cim_dp.SinglePhaseKind.B: + ratings.b = power + if phase.phase == cim_dp.SinglePhaseKind.C: + ratings.c = power + + compensators.ratings[mrid] = ratings + + measurement = cim_dp.Measurement + for measurement in shunt.Measurements: + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + if measurement.measurementType == "PNV": + s1 = measurement.phases == cim_dp.PhaseCode.s1 + s2 = measurement.phases == cim_dp.PhaseCode.s2 + if not (s1 or s2): + print("IGNORING: ", mrid, measurement.name, + measurement.measurementType, measurement.phases) + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + compensators.measurement_map[measurement.mRID] = info + + compensators.measurements_va[mrid] = models.PhasePower() + compensators.measurements_pnv[mrid] = models.PhaseMap() + + return compensators + + +def get_consumers(network: cimgraph.GraphModel) -> models.Consumers: + consumers = models.Consumers() + if cim_dp.EnergyConsumer not in network.graph: + return consumers + + load: cim_dp.EnergyConsumer + for load in network.graph[cim_dp.EnergyConsumer].values(): + mrid = load.mRID + + if not load.EnergyConsumerPhase: + measurement = cim_dp.Measurement + for measurement in load.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + consumers.measurement_map[measurement.mRID] = info + + consumers.measurements_va[mrid] = models.PhasePower() + consumers.measurements_pnv[mrid] = models.PhaseMap() + + p = float(load.p) + q = float(load.q) + power = models.ComplexPower(real=p/3, imag=q/3) + ratings = models.PhasePower(a=power, b=power, c=power) + consumers.ratings[mrid] = ratings + + else: + ratings = models.PhasePower() + phase: cim_dp.EnergyConsumerPhase + + for phase in load.EnergyConsumerPhase: + power = models.ComplexPower(real=phase.p, imag=phase.q) + + if phase.phase == cim_dp.SinglePhaseKind.A: + ratings.a = power + if phase.phase == cim_dp.SinglePhaseKind.B: + ratings.b = power + if phase.phase == cim_dp.SinglePhaseKind.C: + ratings.c = power + + consumers.ratings[mrid] = ratings + + measurement = cim_dp.Measurement + for measurement in load.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + if measurement.measurementType == "PNV": + s1 = measurement.phases == cim_dp.PhaseCode.s1 + s2 = measurement.phases == cim_dp.PhaseCode.s2 + if not (s1 or s2): + print("IGNORING: ", mrid, measurement.name, + measurement.measurementType, measurement.phases) + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + consumers.measurement_map[measurement.mRID] = info + + consumers.measurements_va[mrid] = models.PhasePower() + consumers.measurements_pnv[mrid] = models.PhaseMap() + + return consumers + + +def get_switches(network: cimgraph.GraphModel) -> models.Switches: + switches = models.Switches() + if cim_dp.LoadBreakSwitch not in network.graph: + return switches + + switch: cim_dp.LoadBreakSwitch + for switch in network.graph[cim_dp.LoadBreakSwitch].values(): + mrid = switch.mRID + + if not switch.SwitchPhase: + measurement = cim_dp.Measurement + for measurement in switch.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + switches.measurement_map[measurement.mRID] = info + + switches.measurements_va[mrid] = models.PhaseSwitch() + switches.measurements_pnv[mrid] = models.PhaseMap() + + value = str(switch.normalOpen).capitalize() == "True" + switches.normal[mrid] = models.PhaseSwitch( + a=value, b=value, c=value) + else: + normal = models.PhaseSwitch() + phase: cim_dp.SwitchPhase + for phase in switch.SwitchPhase: + value = str(switch.normalOpen).capitalize() == "True" + if phase.phaseSide1 == cim_dp.SinglePhaseKind.A: + normal.a = value + if phase.phaseSide1 == cim_dp.SinglePhaseKind.B: + normal.b = value + if phase.phaseSide1 == cim_dp.SinglePhaseKind.C: + normal.c = value + switches.normal[mrid] = normal + + measurement = cim_dp.Measurement + for measurement in switch.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + if measurement.measurementType == "PNV": + s1 = measurement.phases == cim_dp.PhaseCode.s1 + s2 = measurement.phases == cim_dp.PhaseCode.s2 + if not (s1 or s2): + print("IGNORING: ", mrid, measurement.name, + measurement.measurementType, measurement.phases) + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + switches.measurement_map[measurement.mRID] = info + + switches.measurements_va[mrid] = models.PhaseSwitch() + switches.measurements_pnv[mrid] = models.PhaseMap() + + return switches + + +def get_power_electronics(network: cimgraph.GraphModel) -> models.PowerElectronics: + electronics = models.PowerElectronics() + if cim_dp.PowerElectronicsConnection not in network.graph: + return electronics + + pec: cim_dp.PowerElectronicsConnection + for pec in network.graph[cim_dp.PowerElectronicsConnection].values(): + mrid = pec.mRID + + for unit in pec.PowerElectronicsUnit: + electronics.units[mrid] = unit.mRID + + if not pec.PowerElectronicsConnectionPhases: + measurement = cim_dp.Measurement + for measurement in pec.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + electronics.measurement_map[measurement.mRID] = info + + electronics.measurements_va[mrid] = models.PhasePower() + electronics.measurements_pnv[mrid] = models.PhaseMap() + + rated_s = float(pec.ratedS) + power = models.ComplexPower(real=rated_s/3, imag=rated_s/3) + ratings = models.PhasePower(a=power, b=power, c=power) + electronics.ratings[mrid] = ratings + + else: + ratings = models.PhasePower() + phase: cim_dp.PowerElectronicsConnectionPhase + for phase in pec.PowerElectronicsConnectionPhases: + power = models.ComplexPower( + real=float(phase.p), imag=float(phase.q)) + if phase.phase == cim_dp.SinglePhaseKind.A: + ratings.a = power + if phase.phase == cim_dp.SinglePhaseKind.B: + ratings.b = power + if phase.phase == cim_dp.SinglePhaseKind.C: + ratings.c = power + electronics.ratings[mrid] = ratings + + measurement = cim_dp.Measurement + for measurement in pec.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + if measurement.measurementType == "PNV": + s1 = measurement.phases == cim_dp.PhaseCode.s1 + s2 = measurement.phases == cim_dp.PhaseCode.s2 + if not (s1 or s2): + print("IGNORING: ", mrid, measurement.name, + measurement.measurementType, measurement.phases) + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + electronics.measurement_map[measurement.mRID] = info + + electronics.measurements_va[mrid] = models.PhasePower() + electronics.measurements_pnv[mrid] = models.PhaseMap() + + return electronics + + +def get_Transformers(network: cimgraph.GraphModel) -> models.Transformers: + transformers = models.Transformers() + if cim_dp.PowerTransformer not in network.graph: + return transformers + + xfmr: cim_dp.PowerTransformer + for xfmr in network.graph[cim_dp.PowerTransformer].values(): + mrid = xfmr.mRID + + measurement = cim_dp.Measurement + for measurement in xfmr.Measurements: + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + transformers.measurement_map[measurement.mRID] = info + + transformers.measurements_va[mrid] = models.PhasePower() + transformers.measurements_pnv[mrid] = models.PhaseMap() + + return transformers + + +def get_power_electronics(network: cimgraph.GraphModel) -> models.PowerElectronics: + electronics = models.PowerElectronics() + if cim_dp.PowerElectronicsConnection not in network.graph: + return electronics + + pec: cim_dp.PowerElectronicsConnection + for pec in network.graph[cim_dp.PowerElectronicsConnection].values(): + mrid = pec.mRID + + for unit in pec.PowerElectronicsUnit: + electronics.units[mrid] = unit.mRID + + if not pec.PowerElectronicsConnectionPhases: + measurement = cim_dp.Measurement + for measurement in pec.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + electronics.measurement_map[measurement.mRID] = info + + electronics.measurements_va[mrid] = models.PhasePower() + electronics.measurements_pnv[mrid] = models.PhaseMap() + + rated_s = float(pec.ratedS) + power = models.ComplexPower(real=rated_s/3, imag=rated_s/3) + ratings = models.PhasePower(a=power, b=power, c=power) + electronics.ratings[mrid] = ratings + + else: + ratings = models.PhasePower() + phase: cim_dp.PowerElectronicsConnectionPhase + for phase in pec.PowerElectronicsConnectionPhases: + power = models.ComplexPower( + real=float(phase.p), imag=float(phase.q)) + if phase.phase == cim_dp.SinglePhaseKind.A: + ratings.a = power + if phase.phase == cim_dp.SinglePhaseKind.B: + ratings.b = power + if phase.phase == cim_dp.SinglePhaseKind.C: + ratings.c = power + electronics.ratings[mrid] = ratings + + measurement = cim_dp.Measurement + for measurement in pec.Measurements: + + pnv = measurement.measurementType == "PNV" + va = measurement.measurementType == "VA" + if not (pnv or va): + continue + + if measurement.measurementType == "PNV": + s1 = measurement.phases == cim_dp.PhaseCode.s1 + s2 = measurement.phases == cim_dp.PhaseCode.s2 + if not (s1 or s2): + print("IGNORING: ", mrid, measurement.name, + measurement.measurementType, measurement.phases) + continue + + info = models.MeasurementInfo( + mrid=mrid, value_type=measurement.measurementType, phase=measurement.phases) + electronics.measurement_map[measurement.mRID] = info + + electronics.measurements_va[mrid] = models.PhasePower() + electronics.measurements_pnv[mrid] = models.PhaseMap() + + return electronics + + +def map_power_electronics(network: cimgraph.GraphModel) -> dict: + map = {} + if cim_dp.PowerElectronicsConnection not in network.graph: + return map + + pec: cim_dp.PowerElectronicsConnection + for pec in network.graph[cim_dp.PowerElectronicsConnection].values(): + name = pec.mRID + bus = pec.Terminals[0].ConnectivityNode.name + map[name] = bus + return map diff --git a/micro-apps/power_factor/test.py b/micro-apps/power_factor/test.py new file mode 100644 index 0000000..ab9d5f0 --- /dev/null +++ b/micro-apps/power_factor/test.py @@ -0,0 +1,73 @@ +import numpy as np +import cvxpy as cp + + +if __name__ == "__main__": + # Problem data. + cost = [[0.14898845, 0.14898845, 0.14898845], + [0.19368499, 0.19368499, 0.19368499], + [0.1638873, 0.1638873, 0.1638873], + [0.19368499, 0.19368499, 0.19368499], + [0.22348268, 0.22348268, 0.22348268], + [0.26817921, 0.26817921, 0.26817921], + [0.25328037, 0.25328037, 0.25328037], + [0.22348268, 0.22348268, 0.22348268], + [0.28307806, 0.28307806, 0.28307806], + [0.34267344, 0.34267344, 0.34267344], + [0.1638873, 0.1638873, 0.1638873], + [0.2979769, 0.2979769, 0.2979769], + [0.28307806, 0.28307806, 0.28307806], + [0.04469654, 0.04469654, 0.04469654], + [0.34267344, 0.34267344, 0.34267344], + [0.26817921, 0.26817921, 0.26817921], + [0.11919076, 0.11919076, 0.11919076], + [0.1638873, 0.1638873, 0.1638873], + [0.11919076, 0.11919076, 0.11919076]] + cost = np.array(cost) + bounds = [[250000., 0., 0.], + [0., 150000., 0.], + [250000., 0., 0.], + [0., 0., 300000.], + [0., 260000., 0.], + [120000., 0., 0.], + [350000., 0., 0.], + [150000., 0., 0.], + [0., 400000., 0.], + [0., 200000., 0.], + [250000., 0., 0.], + [130000., 0., 0.], + [0., 0., 280000.], + [0., 150000., 0.], + [0., 0., 100000.], + [0., 0., 120000.], + [0., 0., 260000.], + [125000., 0., 0.], + [0., 300000., 0.]] + bounds = np.array(bounds) +A = np.clip(bounds, a_min=0, a_max=1) +b = [274999.99995318, 14999.99995318, 129999.99995318] +b = np.array(b) + +# Construct the problem. +m, n = np.shape(A) +x = cp.Variable((m, n)) + +# *, +, -, / are overloaded to construct CVXPY objects. +objective = cp.Minimize(cp.sum(cost.T@cp.abs(x))) +# <=, >=, == are overloaded to construct CVXPY constraints. +constraints = [cp.sum(A.T@x, axis=0) == b, -bounds <= x, x <= bounds] +prob = cp.Problem(objective, constraints) + +# The optimal objective is returned by prob.solve(). +result = prob.solve(solver=cp.CLARABEL, verbose=True) +print(result) + +# The optimal value for x is stored in x.value. +print(b) +print(cost) +print(np.round(x.value, 0)) +print(np.sum(A.T@x.value, axis=0)) + +# The optimal Lagrange multiplier for a constraint +# is stored in constraint.dual_value. +# print(constraints[0].dual_value)