diff --git a/src/GridCal/Gui/Main/SubClasses/Model/diagrams.py b/src/GridCal/Gui/Main/SubClasses/Model/diagrams.py index 0a5c18592..e5c0f6950 100644 --- a/src/GridCal/Gui/Main/SubClasses/Model/diagrams.py +++ b/src/GridCal/Gui/Main/SubClasses/Model/diagrams.py @@ -52,7 +52,6 @@ from GridCal.Gui.Diagrams.MapWidget.Tiles.TileProviders.blue_marble import BlueMarbleTiles from GridCal.Gui.Diagrams.MapWidget.Tiles.TileProviders.cartodb import CartoDbTiles from GridCal.Gui.Diagrams.MapWidget.Tiles.TileProviders.open_street_map import OsmTiles -from trunk.qt_related.slippy_map import MapWidget ALL_EDITORS = Union[SchematicWidget, GridMapWidget] ALL_EDITORS_NONE = Union[None, SchematicWidget, GridMapWidget] diff --git a/src/GridCal/Gui/Main/icons_rc.py b/src/GridCal/Gui/Main/icons_rc.py index 2514ba890..09dc88c67 100644 --- a/src/GridCal/Gui/Main/icons_rc.py +++ b/src/GridCal/Gui/Main/icons_rc.py @@ -38287,7 +38287,7 @@ \x00\x00\x12\xb6\x00\x01\x00\x00\x00\x01\x00\x06\xc5\xc7\ \x00\x00\x01\x88\xae\xf9[%\ \x00\x00\x07\xa8\x00\x00\x00\x00\x00\x01\x00\x02\xda$\ -\x00\x00\x01\x92\xb3\xa2\x89y\ +\x00\x00\x01\x92\xb8\x97\x16\xad\ \x00\x00\x010\x00\x00\x00\x00\x00\x01\x00\x00N\xe0\ \x00\x00\x01\x90N\xda\x14\x95\ \x00\x00\x12\x84\x00\x01\x00\x00\x00\x01\x00\x06\xb0\xef\ diff --git a/src/GridCal/__version__.py b/src/GridCal/__version__.py index cf9fa03bf..f692c7465 100644 --- a/src/GridCal/__version__.py +++ b/src/GridCal/__version__.py @@ -16,7 +16,7 @@ _current_year_ = datetime.datetime.now().year # do not forget to keep a three-number version!!! -__GridCal_VERSION__ = "5.1.53" +__GridCal_VERSION__ = "5.1.55" url = 'https://github.com/SanPen/GridCal' diff --git a/src/GridCalEngine/Devices/Branches/transformer3w.py b/src/GridCalEngine/Devices/Branches/transformer3w.py index b0710b01b..381f18ced 100644 --- a/src/GridCalEngine/Devices/Branches/transformer3w.py +++ b/src/GridCalEngine/Devices/Branches/transformer3w.py @@ -29,6 +29,7 @@ def delta_to_star(z12: float, z23: float, z31: float) -> Tuple[float, float, float]: """ Perform the delta->star transformation + See: https://www.electronics-tutorials.ws/dccircuits/dcp_10.html :param z12: 1 to 2 delta value :param z23: 2 to 3 delta value :param z31: 3 to 1 delta value @@ -44,17 +45,19 @@ def delta_to_star(z12: float, z23: float, z31: float) -> Tuple[float, float, flo return 1e-20, 1e-20, 1e-20 -def start_to_delta(z1: float, z2: float, z3: float): +def star_to_delta(z1: float, z2: float, z3: float) -> Tuple[float, float, float]: """ Perform the star->delta transformation - :param z1: - :param z2: - :param z3: - :return: + See: https://www.electronics-tutorials.ws/dccircuits/dcp_10.html + :param z1: 0->1 impedance + :param z2: 0->2 impedance + :param z3: 0->3 impedance + :return: z12, z23, z31 impedances """ - z12 = z1 + z2 - z23 = z2 + z3 - z31 = z3 + z1 + zt = z1 * z2 + z2 * z3 + z3 * z1 + z12 = zt / z3 + z23 = zt / z2 + z31 = zt / z1 return z12, z23, z31 @@ -427,8 +430,8 @@ def fill_from_star(self, r1: float, r2: float, r3: float, x1: float, x2: float, :param x2: reactance of the branch 2 (p.u.) :param x3: reactance of the branch 3 (p.u.) """ - self._r12, self._r23, self._r31 = start_to_delta(z1=r1, z2=r2, z3=r3) - self._x12, self._x23, self._x31 = start_to_delta(z1=x1, z2=x2, z3=x3) + self._r12, self._r23, self._r31 = star_to_delta(z1=r1, z2=r2, z3=r3) + self._x12, self._x23, self._x31 = star_to_delta(z1=x1, z2=x2, z3=x3) self.winding1.R = r1 self.winding1.X = x1 diff --git a/src/GridCalEngine/IO/cim/cgmes/cgmes_v3_0_0/devices/sv_power_flow.py b/src/GridCalEngine/IO/cim/cgmes/cgmes_v3_0_0/devices/sv_power_flow.py index e27eb7a40..c957ce1ff 100644 --- a/src/GridCalEngine/IO/cim/cgmes/cgmes_v3_0_0/devices/sv_power_flow.py +++ b/src/GridCalEngine/IO/cim/cgmes/cgmes_v3_0_0/devices/sv_power_flow.py @@ -14,17 +14,19 @@ # You should have received a copy of the GNU Lesser General Public License # along with this program; if not, write to the Free Software Foundation, # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +from __future__ import annotations + from GridCalEngine.IO.base.units import UnitMultiplier, UnitSymbol from GridCalEngine.IO.cim.cgmes.base import Base -from GridCalEngine.IO.cim.cgmes.cgmes_enums import cgmesProfile, UnitSymbol +from GridCalEngine.IO.cim.cgmes.cgmes_enums import UnitSymbol class SvPowerFlow(Base): - def __init__(self, rdfid, tpe, resources=list(), class_replacements=dict()): + def __init__(self, rdfid, tpe="SvPowerFlow", resources=list(), class_replacements=dict()): Base.__init__(self, rdfid=rdfid, tpe=tpe, resources=resources, class_replacements=class_replacements) - self.p: float = None - self.q: float = None + self.p: float = 0.0 + self.q: float = 0.0 from GridCalEngine.IO.cim.cgmes.cgmes_v3_0_0.devices.terminal import Terminal self.Terminal: Terminal | None = None diff --git a/src/GridCalEngine/IO/cim/cgmes/gridcal_to_cgmes.py b/src/GridCalEngine/IO/cim/cgmes/gridcal_to_cgmes.py index 16f48ce30..d4db71d3c 100644 --- a/src/GridCalEngine/IO/cim/cgmes/gridcal_to_cgmes.py +++ b/src/GridCalEngine/IO/cim/cgmes/gridcal_to_cgmes.py @@ -1,35 +1,46 @@ +# GridCal +# Copyright (C) 2015 - 2024 Santiago Peñate Vera +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 3 of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. import numpy as np -from GridCalEngine.DataStructures.numerical_circuit import \ - compile_numerical_circuit_at, NumericalCircuit +from datetime import datetime +from GridCalEngine.DataStructures.numerical_circuit import (compile_numerical_circuit_at, NumericalCircuit) from GridCalEngine.Devices import MultiCircuit from GridCalEngine.Devices.Substation.bus import Bus -from GridCalEngine.IO.cim.cgmes.base import get_new_rdfid, form_rdfid, \ - rfid2uuid +from GridCalEngine.IO.cim.cgmes.base import get_new_rdfid, form_rdfid from GridCalEngine.IO.cim.cgmes.cgmes_circuit import CgmesCircuit from GridCalEngine.IO.cim.cgmes.cgmes_enums import (cgmesProfile, WindGenUnitKind, RegulatingControlModeKind, UnitMultiplier, TransformerControlMode) -from GridCalEngine.IO.cim.cgmes.cgmes_utils import find_object_by_uuid, \ - find_object_by_vnom, find_object_by_cond_eq_uuid, \ - get_ohm_values_power_transformer, find_tn_by_name, find_object_by_attribute -from GridCalEngine.IO.cim.cgmes.cgmes_v2_4_15.devices.full_model import \ - FullModel +from GridCalEngine.IO.cim.cgmes.cgmes_utils import (find_object_by_uuid, + find_object_by_vnom, find_object_by_cond_eq_uuid, + get_ohm_values_power_transformer, find_tn_by_name, + find_object_by_attribute) +from GridCalEngine.IO.cim.cgmes.cgmes_v2_4_15.devices.full_model import FullModel from GridCalEngine.IO.cim.cgmes.base import Base import GridCalEngine.Devices as gcdev -from GridCalEngine.Simulations.PowerFlow.NumericalMethods.common_functions import \ - compute_zip_power -from GridCalEngine.Simulations.PowerFlow.power_flow_results import \ - PowerFlowResults +from GridCalEngine.Simulations.PowerFlow.NumericalMethods.common_functions import compute_zip_power +from GridCalEngine.Simulations.PowerFlow.power_flow_results import PowerFlowResults from GridCalEngine.enumerations import CGMESVersions, TapChangerTypes -from GridCalEngine.IO.cim.cgmes.cgmes_enums import ( - SynchronousMachineOperatingMode, - SynchronousMachineKind) +from GridCalEngine.IO.cim.cgmes.cgmes_enums import (SynchronousMachineOperatingMode, SynchronousMachineKind) from GridCalEngine.data_logger import DataLogger -from typing import Dict, List, Tuple, Union +from typing import List, Union # region create new classes for CC @@ -39,8 +50,16 @@ def create_cgmes_headers(cgmes_model: CgmesCircuit, desc: str = "", scenariotime: str = "", modelingauthorityset: str = "", version: str = ""): - from datetime import datetime + """ + :param cgmes_model: + :param profiles_to_export: + :param desc: + :param scenariotime: + :param modelingauthorityset: + :param version: + :return: + """ if cgmes_model.cgmes_version == CGMESVersions.v2_4_15: fm_list = [FullModel(rdfid=get_new_rdfid(), tpe="FullModel"), FullModel(rdfid=get_new_rdfid(), tpe="FullModel"), @@ -63,7 +82,7 @@ def create_cgmes_headers(cgmes_model: CgmesCircuit, fm.scenarioTime = scenariotime if modelingauthorityset != "": fm.modelingAuthoritySet = modelingauthorityset - current_time = datetime.utcnow() + current_time = datetime.now() formatted_time = current_time.strftime("%Y-%m-%dT%H:%M:%S.%fZ") fm.created = formatted_time fm.version = version @@ -154,33 +173,26 @@ def create_cgmes_headers(cgmes_model: CgmesCircuit, return cgmes_model -# def create_multic_cn_nodes_from_buses(multi_circuit_model: MultiCircuit, -# logger: DataLogger) -> None: -# -# for cgmes_elm in multi_circuit_model.buses: -# gcdev_elm = gcdev.ConnectivityNode( -# idtag=cgmes_elm.uuid, -# code=cgmes_elm.description, -# name=cgmes_elm.name, -# dc=False, -# default_bus=bus -# ) -# -# multi_circuit_model.connectivity_nodes.append(gcdev_elm) -# -# def create_cgmes_terminal(mc_bus: Bus, seq_num: Union[int, None], cond_eq: Union[None, Base], cgmes_model: CgmesCircuit, logger: DataLogger): - """ Creates a new Terminal in CGMES model, - and connects it the relating Topologinal Node """ + """ + Creates a new Terminal in CGMES model, + and connects it the relating Topologinal Node + :param mc_bus: + :param seq_num: + :param cond_eq: + :param cgmes_model: + :param logger: + :return: + """ new_rdf_id = get_new_rdfid() terminal_template = cgmes_model.get_class_type("Terminal") term = terminal_template(new_rdf_id) - term.name = f'{cond_eq.name} - T{seq_num}' + term.name = f'{cond_eq.name} - T{seq_num}' if cond_eq is not None else "" # term.shortName = if seq_num is not None: term.sequenceNumber = seq_num @@ -215,6 +227,13 @@ def create_cgmes_terminal(mc_bus: Bus, def create_cgmes_load_response_char(load: gcdev.Load, cgmes_model: CgmesCircuit, logger: DataLogger): + """ + + :param load: + :param cgmes_model: + :param logger: + :return: + """ new_rdf_id = get_new_rdfid() lrc_template = cgmes_model.get_class_type("LoadResponseCharacteristic") lrc = lrc_template(rdfid=new_rdf_id) @@ -434,6 +453,13 @@ def create_cgmes_current_limit(terminal, def create_operational_limit_set(terminal, cgmes_model: CgmesCircuit, logger: DataLogger): + """ + + :param terminal: + :param cgmes_model: + :param logger: + :return: + """ new_rdf_id = get_new_rdfid() object_template = cgmes_model.get_class_type("OperationalLimitSet") op_lim_set = object_template(rdfid=new_rdf_id) @@ -452,6 +478,13 @@ def create_operational_limit_set(terminal, def create_operational_limit_type(mc_elm: gcdev.Line, cgmes_model: CgmesCircuit, logger: DataLogger): + """ + + :param mc_elm: + :param cgmes_model: + :param logger: + :return: + """ new_rdf_id = get_new_rdfid() object_template = cgmes_model.get_class_type("OperationalLimitSet") op_lim_type = object_template(rdfid=new_rdf_id) @@ -465,6 +498,15 @@ def add_location(cgmes_model: CgmesCircuit, longitude: float, latitude: float, logger: DataLogger): + """ + + :param cgmes_model: + :param device: + :param longitude: + :param latitude: + :param logger: + :return: + """ object_template = cgmes_model.get_class_type("Location") location = object_template(rdfid=get_new_rdfid(), tpe="Location") @@ -495,6 +537,13 @@ def add_location(cgmes_model: CgmesCircuit, def get_cgmes_geograpical_regions(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger): + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ for mc_class in [multi_circuit_model.countries, multi_circuit_model.areas]: for mc_elm in mc_class: object_template = cgmes_model.get_class_type("GeographicalRegion") @@ -513,6 +562,13 @@ def get_cgmes_geograpical_regions(multi_circuit_model: MultiCircuit, def get_cgmes_subgeograpical_regions(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger): + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ for mc_class in [multi_circuit_model.communities, multi_circuit_model.zones]: for mc_elm in mc_class: @@ -540,7 +596,7 @@ def get_cgmes_subgeograpical_regions(multi_circuit_model: MultiCircuit, else: try: sub_geo_region.Region = \ - cgmes_model.cgmes_assets.GeographicalRegion_list[0] + cgmes_model.cgmes_assets.GeographicalRegion_list[0] except: sub_geo_region.Region = None logger.add_warning( @@ -556,6 +612,12 @@ def get_cgmes_subgeograpical_regions(multi_circuit_model: MultiCircuit, def get_base_voltage_from_boundary(cgmes_model: CgmesCircuit, vnom: float): + """ + + :param cgmes_model: + :param vnom: + :return: + """ bv_list = cgmes_model.elements_by_type_boundary.get("BaseVoltage") if bv_list is not None: for bv in bv_list: @@ -567,6 +629,13 @@ def get_base_voltage_from_boundary(cgmes_model: CgmesCircuit, vnom: float): def get_cgmes_base_voltages(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger) -> None: + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ base_volt_set = set() for bus in multi_circuit_model.buses: @@ -587,6 +656,13 @@ def get_cgmes_base_voltages(multi_circuit_model: MultiCircuit, def get_cgmes_substations(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger) -> None: + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ for mc_elm in multi_circuit_model.substations: object_template = cgmes_model.get_class_type("Substation") substation = object_template(rdfid=form_rdfid(mc_elm.idtag)) @@ -602,8 +678,7 @@ def get_cgmes_substations(multi_circuit_model: MultiCircuit, substation.Region = region else: try: - substation.Region = \ - cgmes_model.cgmes_assets.SubGeographicalRegion_list[0] + substation.Region = cgmes_model.cgmes_assets.SubGeographicalRegion_list[0] except: substation.Region = None logger.add_warning(msg='Region not found for Substation', @@ -618,6 +693,13 @@ def get_cgmes_substations(multi_circuit_model: MultiCircuit, def get_cgmes_voltage_levels(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger) -> None: + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ for mc_elm in multi_circuit_model.voltage_levels: object_template = cgmes_model.get_class_type("VoltageLevel") vl = object_template(rdfid=form_rdfid(mc_elm.idtag)) @@ -655,6 +737,13 @@ def get_cgmes_voltage_levels(multi_circuit_model: MultiCircuit, def get_cgmes_tn_nodes(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger) -> None: + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ for bus in multi_circuit_model.buses: if not bus.is_internal: @@ -701,6 +790,13 @@ def get_cgmes_tn_nodes(multi_circuit_model: MultiCircuit, def get_cgmes_cn_nodes_from_tn_nodes(multi_circuit_model: MultiCircuit, cgmes_model: CgmesCircuit, logger: DataLogger) -> None: + """ + + :param multi_circuit_model: + :param cgmes_model: + :param logger: + :return: + """ for tn in cgmes_model.cgmes_assets.TopologicalNode_list: new_rdf_id = get_new_rdfid() object_template = cgmes_model.get_class_type("ConnectivityNode") @@ -862,7 +958,7 @@ def get_cgmes_ac_line_segments(multicircuit_model: MultiCircuit, line.length = mc_elm.length current_rate = mc_elm.rate * 1e3 / ( - mc_elm.get_max_bus_nominal_voltage() * 1.73205080756888) + mc_elm.get_max_bus_nominal_voltage() * 1.73205080756888) current_rate = np.round(current_rate, 4) create_cgmes_current_limit(line.Terminals[0], current_rate, cgmes_model, logger) @@ -1058,7 +1154,7 @@ def get_cgmes_power_transformers(multicircuit_model: MultiCircuit, ) current_rate = mc_elm.rate * 1e3 / ( - mc_elm.get_max_bus_nominal_voltage() * 1.73205080756888) + mc_elm.get_max_bus_nominal_voltage() * 1.73205080756888) current_rate = np.round(current_rate, 4) create_cgmes_current_limit(cm_transformer.Terminals[0], current_rate, cgmes_model, logger) @@ -1417,9 +1513,9 @@ def get_cgmes_sv_voltages(cgmes_model: CgmesCircuit, sv_voltage = object_template(rdfid=new_rdf_id, tpe='SvVoltage') tp_list_with_boundary = ( - cgmes_model.cgmes_assets.TopologicalNode_list + - cgmes_model.elements_by_type_boundary.get( - 'TopologicalNode', None)) + cgmes_model.cgmes_assets.TopologicalNode_list + + cgmes_model.elements_by_type_boundary.get( + 'TopologicalNode', None)) sv_voltage.TopologicalNode = tp_list_with_boundary[i] # as the ORDER of the results is the same as the order of buses (=tn) @@ -1492,12 +1588,12 @@ def get_cgmes_sv_power_flow(multi_circuit: MultiCircuit, Y0=nc.load_data.Y, Vm=np.abs(pf_results.voltage[nc.load_data.get_bus_indices()]) ) - + for (mc_load_like, load_power) in zip(load_objects, load_power): term = find_object_by_cond_eq_uuid( object_list=cgmes_model.cgmes_assets.Terminal_list, - cond_eq_target_uuid=mc_load_like.idtag # missing Load uuid + cond_eq_target_uuid=mc_load_like.idtag # missing Load uuid ) if isinstance(term, cgmes_model.get_class_type("Terminal")): @@ -1620,7 +1716,6 @@ def get_cgmes_sv_tap_step(multi_circuit: MultiCircuit, cgmes_model: CgmesCircuit, pf_results: PowerFlowResults, logger: DataLogger) -> None: - branch_objects = multi_circuit.get_branches_wo_hvdc() tap_modules = pf_results.tap_module diff --git a/src/GridCalEngine/Utils/MIP/pulp_interface.py b/src/GridCalEngine/Utils/MIP/pulp_interface.py index 28289d014..ce2d32e49 100644 --- a/src/GridCalEngine/Utils/MIP/pulp_interface.py +++ b/src/GridCalEngine/Utils/MIP/pulp_interface.py @@ -209,7 +209,12 @@ def solve(self, robust: bool = False, show_logs: bool = False, progress_text(f"Solving model with {self.solver_type.value}...") # solve the model - status = self.model.solve(solver=self.get_solver(show_logs=show_logs)) + try: + status = self.model.solve(solver=self.get_solver(show_logs=show_logs)) + except pulp.PulpSolverError as e: + self.logger.add_error(msg=str(e), ) + # Retry with Highs + status = self.model.solve(solver=HiGHS(mip=self.model.isMIP(), msg=show_logs)) if status != self.OPTIMAL: self.originally_infeasible = True diff --git a/src/GridCalEngine/__version__.py b/src/GridCalEngine/__version__.py index 68d871283..401f2b3e5 100644 --- a/src/GridCalEngine/__version__.py +++ b/src/GridCalEngine/__version__.py @@ -16,7 +16,7 @@ _current_year_ = datetime.datetime.now().year # do not forget to keep a three-number version!!! -__GridCalEngine_VERSION__ = "5.1.53" +__GridCalEngine_VERSION__ = "5.1.55" url = 'https://github.com/SanPen/GridCal' diff --git a/src/GridCalServer/__version__.py b/src/GridCalServer/__version__.py index 800d1754e..c7b6eba04 100644 --- a/src/GridCalServer/__version__.py +++ b/src/GridCalServer/__version__.py @@ -16,7 +16,7 @@ _current_year_ = datetime.datetime.now().year # do not forget to keep a three-number version!!! -__GridCalServer_VERSION__ = "5.1.53" +__GridCalServer_VERSION__ = "5.1.55" url = 'https://github.com/SanPen/GridCal' diff --git a/src/tests/test_admittance_tap_derivatives.py b/src/tests/test_admittance_tap_derivatives.py index ab2c90175..11c2c6583 100644 --- a/src/tests/test_admittance_tap_derivatives.py +++ b/src/tests/test_admittance_tap_derivatives.py @@ -18,12 +18,289 @@ import os import numpy as np import GridCalEngine.api as gce +from scipy import sparse as sp from GridCalEngine.enumerations import TapPhaseControl, TapModuleControl from GridCalEngine.DataStructures.numerical_circuit import compile_numerical_circuit_at, NumericalCircuit -from trunk.acopf.acopf_admittance_tap_derivation import (compute_finitediff_admittances, - compute_finitediff_admittances_2dev, - compute_analytic_admittances, - compute_analytic_admittances_2dev) +from GridCalEngine.Simulations.OPF.NumericalMethods.ac_opf import run_nonlinear_opf +from GridCalEngine.Simulations.OPF.opf_options import OptimalPowerFlowOptions + + +def example_3bus_acopf(): + """ + + :return: + """ + + grid = gce.MultiCircuit() + + b1 = gce.Bus(is_slack=True) + b2 = gce.Bus() + b3 = gce.Bus() + + grid.add_bus(b1) + grid.add_bus(b2) + grid.add_bus(b3) + + # grid.add_line(gce.Line(bus_from=b1, bus_to=b2, name='line 1-2', r=0.001, x=0.05, rate=100)) + grid.add_line(gce.Line(bus_from=b2, bus_to=b3, name='line 2-3', r=0.001, x=0.05, rate=100)) + # grid.add_line(gce.Line(bus_from=b3, bus_to=b1, name='line 3-1_1', r=0.001, x=0.05, rate=100)) + # grid.add_line(Line(bus_from=b3, bus_to=b1, name='line 3-1_2', r=0.001, x=0.05, rate=100)) + + grid.add_load(b3, gce.Load(name='L3', P=50, Q=20)) + grid.add_generator(b1, gce.Generator('G1', vset=1.00, Cost=1.0, Cost2=2.0)) + grid.add_generator(b2, gce.Generator('G2', P=10, vset=0.995, Cost=1.0, Cost2=3.0)) + + tr1 = gce.Transformer2W(b1, b2, 'Trafo1', tap_phase_control_mode=TapPhaseControl.Pf, + tap_module=1.1, tap_phase=0.02, r=0.001, x=0.05) + grid.add_transformer2w(tr1) + + tr2 = gce.Transformer2W(b3, b1, 'Trafo1', tap_phase_control_mode=TapPhaseControl.Pf, + tap_module=1.05, tap_phase=-0.02, r=0.001, x=0.05) + grid.add_transformer2w(tr2) + + nc = compile_numerical_circuit_at(circuit=grid) + + A, B, C, D, E, F = compute_analytic_admittances(nc) + + A_, B_, C_, D_, E_, F_ = compute_finitediff_admittances(nc) + + G, H, I, J, K, L, M, N, O, P, Q, R = compute_analytic_admittances_2dev(nc) + + G_, H_, I_, J_, K_, L_, M_, N_, O_, P_, Q_, R_ = compute_finitediff_admittances_2dev(nc) + + options = gce.PowerFlowOptions(gce.SolverType.NR, verbose=False) + power_flow = gce.PowerFlowDriver(grid, options) + power_flow.run() + # print('\n\n', grid.name) + # print('\tConv:\n', power_flow.results.get_bus_df()) + # print('\tConv:\n', power_flow.results.get_branch_df()) + opf_options = OptimalPowerFlowOptions() + + pf_options = gce.PowerFlowOptions(solver_type=gce.SolverType.NR, verbose=3) + run_nonlinear_opf(grid=grid, + opf_options=opf_options, + pf_options=pf_options, + plot_error=True) + + +def compute_analytic_admittances(nc: NumericalCircuit): + """ + + :param nc: + :return: + """ + k_m = nc.k_m + k_tau = nc.k_tau + k_mtau = nc.k_mtau + + tapm = nc.branch_data.tap_module + tapt = nc.branch_data.tap_angle + + Cf = nc.Cf + Ct = nc.Ct + ys = 1.0 / (nc.branch_data.R + 1.0j * nc.branch_data.X + 1e-20) + + # First partial derivative with respect to tap module + mp = tapm[k_m] + tau = tapt[k_m] + ylin = ys[k_m] + + dYffdm = np.zeros(len(tapm), dtype=complex) + dYftdm = np.zeros(len(tapm), dtype=complex) + dYtfdm = np.zeros(len(tapm), dtype=complex) + dYttdm = np.zeros(len(tapm), dtype=complex) + + dYffdm[k_m] = -2 * ylin / (mp * mp * mp) + dYftdm[k_m] = ylin / (mp * mp * np.exp(-1.0j * tau)) + dYtfdm[k_m] = ylin / (mp * mp * np.exp(1.0j * tau)) + + dYfdm = sp.diags(dYffdm) * Cf + sp.diags(dYftdm) * Ct + dYtdm = sp.diags(dYtfdm) * Cf + sp.diags(dYttdm) * Ct + + dYbusdm = Cf.T * dYfdm + Ct.T * dYtdm # Cf_m.T and Ct_m.T included earlier + + # First partial derivative with respect to tap angle + mp = tapm[k_tau] + tau = tapt[k_tau] + ylin = ys[k_tau] + + dYffdt = np.zeros(len(tapm), dtype=complex) + dYftdt = np.zeros(len(tapm), dtype=complex) + dYtfdt = np.zeros(len(tapm), dtype=complex) + dYttdt = np.zeros(len(tapm), dtype=complex) + + dYftdt[k_tau] = -1j * ylin / (mp * np.exp(-1.0j * tau)) + dYtfdt[k_tau] = 1j * ylin / (mp * np.exp(1.0j * tau)) + + dYfdt = sp.diags(dYffdt) * Cf + sp.diags(dYftdt) * Ct + dYtdt = sp.diags(dYtfdt) * Cf + sp.diags(dYttdt) * Ct + + dYbusdt = Cf.T * dYfdt + Ct.T * dYtdt + + return dYbusdm, dYfdm, dYtdm, dYbusdt, dYfdt, dYtdt + + +def compute_finitediff_admittances(nc: NumericalCircuit, tol=1e-5): + """ + + :param nc: NumericalCircuit + :param tol: tolerance + :return: + """ + indices = nc.get_simulation_indices() + k_m = indices.k_m + k_tau = indices.k_tau + + # base values + adm0 = nc.get_admittance_matrices() + + # Modify tap modules + nc.branch_data.tap_module[k_m] += tol + adm1 = nc.get_admittance_matrices() + nc.branch_data.tap_module[k_m] -= tol + + dYf_dm = (adm1.Yf - adm0.Yf) / tol + dYt_dm = (adm1.Yt - adm0.Yt) / tol + dY_dm = (adm1.Ybus - adm0.Ybus) / tol + + # modify tap angles + nc.branch_data.tap_angle[k_tau] += tol + adm2 = nc.get_admittance_matrices() + nc.branch_data.tap_angle[k_tau] -= tol + + dYf_dt = (adm2.Yf - adm0.Yf) / tol + dYt_dt = (adm2.Yt - adm0.Yt) / tol + dY_dt = (adm2.Ybus - adm0.Ybus) / tol + + return dY_dm, dYf_dm, dYt_dm, dY_dt, dYf_dt, dYt_dt + + +def compute_analytic_admittances_2dev(nc: NumericalCircuit): + """ + + :param nc: + :return: + """ + indices = nc.get_simulation_indices() + k_m = indices.k_m + k_tau = indices.k_tau + k_mtau = indices.k_mtau + + tapm = nc.branch_data.tap_module + tapt = nc.branch_data.tap_angle + + Cf = nc.Cf + Ct = nc.Ct + ys = 1.0 / (nc.branch_data.R + 1.0j * nc.branch_data.X + 1e-20) + + # Second partial derivative with respect to tap module + mp = tapm[k_m] + tau = tapt[k_m] + ylin = ys[k_m] + + # primitived + dyff_dmdm = np.zeros(nc.nbr, dtype=complex) + dyft_dmdm = np.zeros(nc.nbr, dtype=complex) + dytf_dmdm = np.zeros(nc.nbr, dtype=complex) + dytt_dmdm = np.zeros(nc.nbr, dtype=complex) + + dyff_dmdm[k_m] = 6 * ylin / (mp * mp * mp * mp) + dyft_dmdm[k_m] = -2 * ylin / (mp * mp * mp * np.exp(-1.0j * tau)) + dytf_dmdm[k_m] = -2 * ylin / (mp * mp * mp * np.exp(1.0j * tau)) + + dYf_dmdm = (sp.diags(dyff_dmdm) * Cf + sp.diags(dyft_dmdm) * Ct) + dYt_dmdm = (sp.diags(dytf_dmdm) * Cf + sp.diags(dytt_dmdm) * Ct) + + dY_dmdm = (Cf.T * dYf_dmdm + Ct.T * dYt_dmdm) + + # Second partial derivative with respect to tap angle + mp = tapm[k_tau] + tau = tapt[k_tau] + ylin = ys[k_tau] + + # primitives + dyff_dtdt = np.zeros(nc.nbr, dtype=complex) + dyft_dtdt = np.zeros(nc.nbr, dtype=complex) + dytf_dtdt = np.zeros(nc.nbr, dtype=complex) + dytt_dtdt = np.zeros(nc.nbr, dtype=complex) + + dyft_dtdt[k_tau] = ylin / (mp * np.exp(-1.0j * tau)) + dytf_dtdt[k_tau] = ylin / (mp * np.exp(1.0j * tau)) + + dYf_dtdt = sp.diags(dyff_dtdt) * Cf + sp.diags(dyft_dtdt) * Ct + dYt_dtdt = sp.diags(dytf_dtdt) * Cf + sp.diags(dytt_dtdt) * Ct + + dY_dtdt = Cf.T * dYf_dtdt + Ct.T * dYt_dtdt + + # Second partial derivative with respect to both tap module and angle + mp = tapm[k_mtau] + tau = tapt[k_mtau] + ylin = ys[k_mtau] + + # primitives + dyffdmdt = np.zeros(nc.nbr, dtype=complex) + dyft_dmdt = np.zeros(nc.nbr, dtype=complex) + dytf_dmdt = np.zeros(nc.nbr, dtype=complex) + dyttdmdt = np.zeros(nc.nbr, dtype=complex) + + dyft_dmdt[k_mtau] = 1j * ylin / (mp * mp * np.exp(-1.0j * tau)) + dytf_dmdt[k_mtau] = -1j * ylin / (mp * mp * np.exp(1.0j * tau)) + + dYf_dmdt = sp.diags(dyffdmdt) * Cf + sp.diags(dyft_dmdt) * Ct + dYt_dmdt = sp.diags(dytf_dmdt) * Cf + sp.diags(dyttdmdt) * Ct + + dY_dmdt = Cf.T * dYf_dmdt + Ct.T * dYt_dmdt + + dYf_dtdm = dYf_dmdt.copy() + dYt_dtdm = dYt_dmdt.copy() + dY_dtdm = dY_dmdt.copy() + + return (dY_dmdm, dYf_dmdm, dYt_dmdm, dY_dmdt, dYf_dmdt, dYt_dmdt, + dY_dtdm, dYf_dtdm, dYt_dtdm, dY_dtdt, dYf_dtdt, dYt_dtdt) + + +def compute_finitediff_admittances_2dev(nc: NumericalCircuit, tol=1e-5): + """ + + :param nc: + :param tol: + :return: + """ + indices = nc.get_simulation_indices() + k_m = indices.k_m + k_tau = indices.k_tau + + # Refference + dY0_dm, dYf0_dm, dYt0_dm, dY0_dt, dYf0_dt, dYt0_dt = compute_finitediff_admittances(nc) + + # Modify the tap module + nc.branch_data.tap_module[k_m] += tol + dY_dm, dYf_dm, dYt_dm, dY_dt, dYf_dt, dYt_dt = compute_finitediff_admittances(nc) + nc.branch_data.tap_module[k_m] -= tol + + dYf_dmdm = (dYf_dm - dYf0_dm) / tol + dYt_dmdm = (dYt_dm - dYt0_dm) / tol + dY_dmdm = (dY_dm - dY0_dm) / tol + + dYf_dtdm = (dYf_dt - dYf0_dt) / tol + dYt_dtdm = (dYt_dt - dYt0_dt) / tol + dY_dtdm = (dY_dt - dY0_dt) / tol + + # Modify the tap angle + nc.branch_data.tap_angle[k_tau] += tol + dY_dm, dYf_dm, dYt_dm, dY_dt, dYf_dt, dYt_dt = compute_finitediff_admittances(nc) + nc.branch_data.tap_angle[k_tau] -= tol + + dYf_dmdt = (dYf_dm - dYf0_dm) / tol + dYt_dmdt = (dYt_dm - dYt0_dm) / tol + dY_dmdt = (dY_dm - dY0_dm) / tol + + dYf_dtdt = (dYf_dt - dYf0_dt) / tol + dYt_dtdt = (dYt_dt - dYt0_dt) / tol + dY_dtdt = (dY_dt - dY0_dt) / tol + + return (dY_dmdm, dYf_dmdm, dYt_dmdm, dY_dmdt, dYf_dmdt, dYt_dmdt, + dY_dtdm, dYf_dtdm, dYt_dtdm, dY_dtdt, dYf_dtdt, dYt_dtdt) def case_3bus() -> NumericalCircuit: