diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 9eba54c..38869fc 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -23,6 +23,30 @@ TODO HVDC in Jacobian (see pandapower) TODO: in ContingencyAnalysisCpp: add back the `if(!ac_solver_used)` inside the `remove_from_Ybus` in order to perform the "invertibility" check TODO: in `main.cpp` check the returned policy of pybind11 and also the `py::call_guard()` stuff +TODO: a cpp class that is able to compute (DC powerflow) ContingencyAnalysis and TimeSeries using PTDF and LODF +TODO: integration test with pandapower (see `pandapower/contingency/contingency.py` and import `lightsim2grid_installed` and check it's True) + +[0.9.1] 2024-xx-yy +-------------------------- +- [FIXED] a bug due to wrong type (in a numpy array) for the element name which lead in turn + to a fail assertion (equality between two numpy arrays returning a bool and not an array) +- [FIXED] a bug when init a grid from pypowsybl: the wrong value was used for trafos `h` (double) +- [FIXED] a bug when init a grid from pypowsybl: wrong values for `_ls_to_orig` and `_orig_to_ls` + was set (and later used) +- [FIXED] yet another bug when init a grid from pypowsybl: the voltage in kV (not in pu) + could be set due to "wrong" labelling of the bus ids +- [FIXED] yet another bug when init a grid from pypowsybl: the ratio of the transformers + sent in lightsim2grid did not take into account the "`rated_u1` `rated_u2`" on both side + (only used on one side) +- [FIXED] yet another bug when init a grid from pypowsybl: the ratio of the transformers + sent in lightsim2grid did not take into account the ratio in the `pypow_net.get_ratio_tap_changers()` +- [ADDED] a method for the `ContingencyAnalysisCPP` class that returns, for all contingencies + in the contingency list, which will be simulated and which causes the grid to be disconnected. +- [ADDED] it is now possible to use "one substation" (voltage level) pypowsybl side is + one substation in lightsim2grid. +- [IMPROVED] removing a weird `1j * h_` when initializing powerlines and transformers. This was + part of a pandapower "hack" which is not present anymore (see + https://github.com/BDonnot/lightsim2grid/issues/88#issue-2443299039) [0.9.0] 2024-07-29 -------------------------- diff --git a/docs/conf.py b/docs/conf.py index 9dc3eb5..cabfe19 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -22,7 +22,7 @@ author = 'Benjamin DONNOT' # The full version, including alpha/beta/rc tags -release = "0.9.0" +release = "0.9.1" version = '0.9' # -- General configuration --------------------------------------------------- diff --git a/docs/lightsimbackend.rst b/docs/lightsimbackend.rst index 62f333c..3f0ac1c 100644 --- a/docs/lightsimbackend.rst +++ b/docs/lightsimbackend.rst @@ -43,7 +43,7 @@ For standard grid2op environment, you can use it like: # the `LightSimBackend` will be used to carry out the powerflow computation instead # of the default grid2op `PandaPowerBackend` -Customization of the backend +Customization of the solver ------------------------------- .. warning:: Use grid2op > 1.7.1 for this feature to work properly. Otherwise some bugs (hard to detect) will occur. @@ -63,6 +63,9 @@ You can customize the way the backend behaves in different ways: - `dist_slack_non_renew`: by default in most grid2op environment, the slack bus is "centralize" / "single slack". This parameters allows to bypass this restriction and use all non renewable generators (and turned on and with > 0.) in a distributed slack bus setting. It might change the default `solver_type` used. +- \* `use_static_gen`: bool=False, DO NOT USE AT THE MOMENT. When it will be available, you will be able to loader_kwargs + both "static" generators (pq generators) and "regular" (pv generators) as generators in lightsim2grid. It does + not work at the moment and has no effect. - \* `detailed_infos_for_cascading_failures`: for exhaustivity, do not modify. - \* `can_be_copied`: for exhaustivity, do not modify. @@ -84,12 +87,56 @@ The easiest way to customize your backend is when you create the grid2op environ ) ) +Customization of the input format +---------------------------------- + +For a few versions now, we try to extend the capability of lightsim2grid and make it work +with other data "reader". We started this process by allowing to initialize a +lightsim2grid `GridModel` from a pypowsybl network. + +For example, if you environment contains a grid in the iidm format (native format of pypowsybl networks), +you can load it with: + +.. code-block:: python + + import grid2op + from lightsim2grid.LightSimBackend import LightSimBackend + from grid2op.Agent import RandomAgent + + # create an environment + env_with_iidm_as_the_grid_description = ... + env = grid2op.make(env_name, + backend=LightSimBackend(loader_method="pypowsybl") + ) + +You can also customize the way lightsim2grid works with some extra options: + + +- `loader_method`: Literal["pandapower", "pypowsybl"]: from which grid "file description" + the grid will be loaded. If you use `pandapower` then pandapower needs to be installed. + If you specified `pypowsybl` then pypowsybl needs to be installed on your machine. +- `loader_kwargs` : ``dict``: some customization to use when loading the grid. It is not + not used when loading the grid from `pandapower`. Please refer to the documentation of + :attr:`LightSimBackend._loader_kwargs` for more information. + +Other Customization +-------------------- + +Here are some other extra features you can use in lightsim2grid (but that are not yet supported by grid2op +so not really usable...) : + +- stop_if_load_disco : Optional[bool] = True: whether to stop the computation (returning a `DivergingPowerflow` exception) + if a load is disconnected. +- stop_if_gen_disco : Optional[bool] = True: whether to stop the computation (returning a `DivergingPowerflow` exception) + if a generator is disconnected. + Detailed documentation -------------------------- .. automodule:: lightsim2grid.lightSimBackend :members: :autosummary: + :private-members: * :ref:`genindex` diff --git a/lightsim2grid/__init__.py b/lightsim2grid/__init__.py index c109ead..8922c7c 100644 --- a/lightsim2grid/__init__.py +++ b/lightsim2grid/__init__.py @@ -6,7 +6,7 @@ # SPDX-License-Identifier: MPL-2.0 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. -__version__ = "0.9.0" +__version__ = "0.9.1" __all__ = ["newtonpf", "SolverType", "ErrorType", "solver", "compilation_options"] diff --git a/lightsim2grid/elements/__init__.py b/lightsim2grid/elements/__init__.py index 7e795ee..d13b34c 100644 --- a/lightsim2grid/elements/__init__.py +++ b/lightsim2grid/elements/__init__.py @@ -6,29 +6,33 @@ # SPDX-License-Identifier: MPL-2.0 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. -__all__ = ["DataGen", +__all__ = ["GeneratorContainer", "GenInfo", - "DataSGen", + "SGenContainer", "SGenInfo", - "DataLoad", + "LoadContainer", "LoadInfo", - "DataShunt", + "ShuntContainer", "ShuntInfo", - "DataTrafo", + "TrafoContainer", "TrafoInfo", - "DataLine", + "LineContainer", "LineInfo", + "DCLineContainer", + "DCLineInfo", ] -from lightsim2grid_cpp import DataGen +from lightsim2grid_cpp import GeneratorContainer from lightsim2grid_cpp import GenInfo -from lightsim2grid_cpp import DataSGen +from lightsim2grid_cpp import SGenContainer from lightsim2grid_cpp import SGenInfo -from lightsim2grid_cpp import DataLoad +from lightsim2grid_cpp import LoadContainer from lightsim2grid_cpp import LoadInfo -from lightsim2grid_cpp import DataShunt +from lightsim2grid_cpp import ShuntContainer from lightsim2grid_cpp import ShuntInfo -from lightsim2grid_cpp import DataTrafo +from lightsim2grid_cpp import TrafoContainer from lightsim2grid_cpp import TrafoInfo -from lightsim2grid_cpp import DataLine +from lightsim2grid_cpp import LineContainer from lightsim2grid_cpp import LineInfo +from lightsim2grid_cpp import DCLineContainer +from lightsim2grid_cpp import DCLineInfo diff --git a/lightsim2grid/gridmodel/from_pandapower/_aux_add_line.py b/lightsim2grid/gridmodel/from_pandapower/_aux_add_line.py index 98accdf..061d9d5 100644 --- a/lightsim2grid/gridmodel/from_pandapower/_aux_add_line.py +++ b/lightsim2grid/gridmodel/from_pandapower/_aux_add_line.py @@ -7,8 +7,11 @@ # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. import numpy as np +from packaging import version +import pandapower as pp from ._pp_bus_to_ls_bus import pp_bus_to_ls +from ._my_const import _MIN_PP_VERSION_ADV_GRID_MODEL def _aux_add_line(converter, model, pp_net, pp_to_ls=None): """ @@ -29,21 +32,41 @@ def _aux_add_line(converter, model, pp_net, pp_to_ls=None): "Some pp_net.line[\"parallel\"] != 1 it is not handled by lightsim yet.") #### find the right powerline parameters - line_r, line_x, line_h = \ - converter.get_line_param( - pp_net.line["r_ohm_per_km"].values * pp_net.line["length_km"].values, - pp_net.line["x_ohm_per_km"].values * pp_net.line["length_km"].values, - pp_net.line["c_nf_per_km"].values * pp_net.line["length_km"].values, - pp_net.line["g_us_per_km"].values * pp_net.line["length_km"].values, - pp_net.bus.loc[pp_net.line["from_bus"]]["vn_kv"], - pp_net.bus.loc[pp_net.line["to_bus"]]["vn_kv"], - ) - - ### add them to the grid - model.init_powerlines(line_r, line_x, line_h, - pp_bus_to_ls(pp_net.line["from_bus"].values, pp_to_ls), - pp_bus_to_ls(pp_net.line["to_bus"].values, pp_to_ls) - ) + if version.parse(pp.__version__) >= _MIN_PP_VERSION_ADV_GRID_MODEL: + # new pandapower with support for different h at both side + line_r, line_x, line_h_or, line_h_ex = \ + converter.get_line_param( + pp_net.line["r_ohm_per_km"].values * pp_net.line["length_km"].values, + pp_net.line["x_ohm_per_km"].values * pp_net.line["length_km"].values, + pp_net.line["g_us_per_km"].values * pp_net.line["length_km"].values, + pp_net.line["c_nf_per_km"].values * pp_net.line["length_km"].values, + pp_net.bus.loc[pp_net.line["from_bus"]]["vn_kv"], + pp_net.bus.loc[pp_net.line["to_bus"]]["vn_kv"], + ) + + ### add them to the grid + model.init_powerlines_full(line_r, line_x, line_h_or, line_h_ex, + pp_bus_to_ls(pp_net.line["from_bus"].values, pp_to_ls), + pp_bus_to_ls(pp_net.line["to_bus"].values, pp_to_ls) + ) + + else: + # legacy pandapower, when they did not support lines with different h both side + line_r, line_x, line_h = \ + converter.get_line_param_legacy( + pp_net.line["r_ohm_per_km"].values * pp_net.line["length_km"].values, + pp_net.line["x_ohm_per_km"].values * pp_net.line["length_km"].values, + pp_net.line["g_us_per_km"].values * pp_net.line["length_km"].values, + pp_net.line["c_nf_per_km"].values * pp_net.line["length_km"].values, + pp_net.bus.loc[pp_net.line["from_bus"]]["vn_kv"], + pp_net.bus.loc[pp_net.line["to_bus"]]["vn_kv"], + ) + + ### add them to the grid + model.init_powerlines(line_r, line_x, line_h, + pp_bus_to_ls(pp_net.line["from_bus"].values, pp_to_ls), + pp_bus_to_ls(pp_net.line["to_bus"].values, pp_to_ls) + ) for line_id, is_connected in enumerate(pp_net.line["in_service"].values): if not is_connected: # powerline is deactivated diff --git a/lightsim2grid/gridmodel/from_pandapower/_aux_add_trafo.py b/lightsim2grid/gridmodel/from_pandapower/_aux_add_trafo.py index f6aa299..1ed2d46 100644 --- a/lightsim2grid/gridmodel/from_pandapower/_aux_add_trafo.py +++ b/lightsim2grid/gridmodel/from_pandapower/_aux_add_trafo.py @@ -8,8 +8,11 @@ import warnings import numpy as np -from ._pp_bus_to_ls_bus import pp_bus_to_ls +import pandapower as pp +from packaging import version +from ._pp_bus_to_ls_bus import pp_bus_to_ls +from ._my_const import _MIN_PP_VERSION_ADV_GRID_MODEL def _aux_add_trafo(converter, model, pp_net, pp_to_ls): """ @@ -69,30 +72,57 @@ def _aux_add_trafo(converter, model, pp_net, pp_to_ls): tap_angles_ = np.deg2rad(tap_angles_) # compute physical parameters - trafo_r, trafo_x, trafo_b = \ - converter.get_trafo_param(tap_step_pct, - tap_pos, - tap_angles_, # in radian ! - is_tap_hv_side, - pp_net.bus.loc[pp_net.trafo["hv_bus"]]["vn_kv"], - pp_net.bus.loc[pp_net.trafo["lv_bus"]]["vn_kv"], - pp_net.trafo["vk_percent"].values, - pp_net.trafo["vkr_percent"].values, - pp_net.trafo["sn_mva"].values, - pp_net.trafo["pfe_kw"].values, - pp_net.trafo["i0_percent"].values, - ) - - # initialize the grid - model.init_trafo(trafo_r, - trafo_x, - trafo_b, - tap_step_pct, - tap_pos, - shift_, - is_tap_hv_side, - pp_bus_to_ls(pp_net.trafo["hv_bus"].values, pp_to_ls), - pp_bus_to_ls(pp_net.trafo["lv_bus"].values, pp_to_ls)) + if version.parse(pp.__version__) >= _MIN_PP_VERSION_ADV_GRID_MODEL: + # TODO + trafo_r, trafo_x, trafo_b = \ + converter.get_trafo_param(tap_step_pct, + tap_pos, + tap_angles_, # in radian ! + is_tap_hv_side, + pp_net.bus.loc[pp_net.trafo["hv_bus"]]["vn_kv"], + pp_net.bus.loc[pp_net.trafo["lv_bus"]]["vn_kv"], + pp_net.trafo["vk_percent"].values, + pp_net.trafo["vkr_percent"].values, + pp_net.trafo["sn_mva"].values, + pp_net.trafo["pfe_kw"].values, + pp_net.trafo["i0_percent"].values, + ) + + # initialize the grid + model.init_trafo(trafo_r, + trafo_x, + trafo_b, + tap_step_pct, + tap_pos, + shift_, + is_tap_hv_side, + pp_bus_to_ls(pp_net.trafo["hv_bus"].values, pp_to_ls), + pp_bus_to_ls(pp_net.trafo["lv_bus"].values, pp_to_ls)) + else: + trafo_r, trafo_x, trafo_b = \ + converter.get_trafo_param_legacy(tap_step_pct, + tap_pos, + tap_angles_, # in radian ! + is_tap_hv_side, + pp_net.bus.loc[pp_net.trafo["hv_bus"]]["vn_kv"], + pp_net.bus.loc[pp_net.trafo["lv_bus"]]["vn_kv"], + pp_net.trafo["vk_percent"].values, + pp_net.trafo["vkr_percent"].values, + pp_net.trafo["sn_mva"].values, + pp_net.trafo["pfe_kw"].values, + pp_net.trafo["i0_percent"].values, + ) + + # initialize the grid + model.init_trafo(trafo_r, + trafo_x, + trafo_b, + tap_step_pct, + tap_pos, + shift_, + is_tap_hv_side, + pp_bus_to_ls(pp_net.trafo["hv_bus"].values, pp_to_ls), + pp_bus_to_ls(pp_net.trafo["lv_bus"].values, pp_to_ls)) for tr_id, is_connected in enumerate(pp_net.trafo["in_service"].values): if not is_connected: diff --git a/lightsim2grid/gridmodel/from_pandapower/_my_const.py b/lightsim2grid/gridmodel/from_pandapower/_my_const.py new file mode 100644 index 0000000..403044e --- /dev/null +++ b/lightsim2grid/gridmodel/from_pandapower/_my_const.py @@ -0,0 +1,11 @@ +# Copyright (c) 2024, RTE (https://www.rte-france.com) +# See AUTHORS.txt +# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0. +# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file, +# you can obtain one at http://mozilla.org/MPL/2.0/. +# SPDX-License-Identifier: MPL-2.0 +# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. + +from packaging import version +# pandapower version with more advanced grid modelling +_MIN_PP_VERSION_ADV_GRID_MODEL = version.parse("2.15") diff --git a/lightsim2grid/gridmodel/from_pypowsybl/_from_pypowsybl.py b/lightsim2grid/gridmodel/from_pypowsybl/_from_pypowsybl.py index 7159d1a..dd06184 100644 --- a/lightsim2grid/gridmodel/from_pypowsybl/_from_pypowsybl.py +++ b/lightsim2grid/gridmodel/from_pypowsybl/_from_pypowsybl.py @@ -7,11 +7,15 @@ # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. import warnings +import copy import numpy as np import pandas as pd import pypowsybl as pypo -from typing import Union +from typing import Optional, Union +from packaging import version +PP_BUG_RATIO_TAP_CHANGER = version.parse("1.7.0") +PYPOWSYBL_VER = version.parse(pypo.__version__) from lightsim2grid_cpp import GridModel @@ -36,22 +40,21 @@ def _aux_get_bus(bus_df, df, conn_key="connected", bus_key="bus_id"): return bus_id, mask_disco.values -def init(net : pypo.network, +def init(net : pypo.network.Network, gen_slack_id: Union[int, str] = None, slack_bus_id: int = None, sn_mva = 100., sort_index=True, f_hz = 50., # unused + net_pu : Optional[pypo.network.Network] = None, only_main_component=True, - return_sub_id=False): + return_sub_id=False, + n_busbar_per_sub=None, # new in 0.9.1 + buses_for_sub=None, # new in 0.9.1 + ): model = GridModel() # model.set_f_hz(f_hz) - # for substation - # network.get_voltage_levels()["substation_id"] - # network.get_substations() - # network.get_busbar_sections() - if gen_slack_id is not None and slack_bus_id is not None: raise RuntimeError("Impossible to intialize a grid with both gen_slack_id and slack_bus_id") @@ -61,29 +64,108 @@ def init(net : pypo.network, bus_df = bus_df_orig.sort_index() else: bus_df = bus_df_orig - bus_df["bus_id"] = np.arange(bus_df.shape[0]) - bus_df_orig["bus_id"] = bus_df.loc[bus_df_orig.index]["bus_id"] - voltage_levels = net.get_voltage_levels() + + if sort_index: + voltage_levels = net.get_voltage_levels().sort_index() + else: + voltage_levels = net.get_voltage_levels() + + all_buses_vn_kv = voltage_levels.loc[bus_df["voltage_level_id"].values]["nominal_v"].values + sub_unique = None + sub_unique_id = None + if n_busbar_per_sub is not None and buses_for_sub is not None: + # I am in a compatibility mode, + # I need to use the same convention as grid2op + # for the buses labelling + bus_df = bus_df.sort_values("voltage_level_id", kind="stable") + sub_unique = bus_df["voltage_level_id"].unique() + nb_sub_unique = sub_unique.shape[0] + sub_unique_id = np.arange(nb_sub_unique) + bus_per_sub = bus_df[["voltage_level_id", "name"]].groupby("voltage_level_id").count() + if (bus_per_sub["name"] > n_busbar_per_sub).any(): + max_bb = bus_per_sub["name"].max() + raise RuntimeError(f"Impossible configuration: we found a substation with {max_bb} " + f"while asking for {n_busbar_per_sub}. We cannot load a grid with these " + f"kwargs. If you use LightSimBackend, you need to change `loader_kwargs` " + f"and especially the `n_busbar_per_sub` to be >= {max_bb}") + bus_local_id = np.concatenate([np.arange(el) for el in bus_per_sub.values]) + sub_id_duplicate = np.repeat(sub_unique_id, bus_per_sub.values.ravel()) + bus_global_id = bus_local_id * nb_sub_unique + sub_id_duplicate + bus_df["bus_id"] = bus_global_id + all_buses_vn_kv = 1. * voltage_levels["nominal_v"].values + ls_to_orig = np.zeros(n_busbar_per_sub * nb_sub_unique, dtype=int) - 1 + ls_to_orig[bus_df["bus_id"].values] = np.arange(bus_df.shape[0]) + n_sub = nb_sub_unique + n_bb_per_sub = n_busbar_per_sub + bus_df = bus_df.sort_index() + else: + if buses_for_sub is not None: + raise NotImplementedError("This is not implemented at the moment") + # retrieve the labels from the buses in the original grid + bus_df["bus_id"] = -1 + bus_df.loc[bus_df_orig.index, "bus_id"] = np.arange(bus_df.shape[0]) + bus_df = bus_df.sort_values("bus_id") + ls_to_orig = 1 * bus_df["bus_id"].values + + n_sub = bus_df.shape[0] + n_bb_per_sub = None + if n_busbar_per_sub is None: + warnings.warn("You should avoid using this function without at least `buses_for_sub` or `n_busbar_per_sub`. " + "Setting automatically n_busbar_per_sub=1") + n_busbar_per_sub = 1 + + # used to be done in the Backend previously, now we do it here instead + bus_init = 1. * all_buses_vn_kv + bus_doubled = np.concatenate([bus_init for _ in range(n_busbar_per_sub)]) + ls_to_orig = np.concatenate((ls_to_orig, np.full((n_busbar_per_sub - 1) * n_sub, fill_value=-1, dtype=ls_to_orig.dtype))) + + all_buses_vn_kv = np.concatenate([all_buses_vn_kv for _ in range(n_busbar_per_sub)]) model.set_sn_mva(sn_mva) model.set_init_vm_pu(1.06) - model.init_bus(voltage_levels.loc[bus_df["voltage_level_id"].values]["nominal_v"].values, + model.init_bus(all_buses_vn_kv, 0, 0 # unused ) - model._orig_to_ls = 1 * bus_df_orig["bus_id"].values + model._ls_to_orig = ls_to_orig + model.set_n_sub(n_sub) + if n_bb_per_sub is not None: + model._max_nb_bus_per_sub = n_busbar_per_sub # do the generators if sort_index: df_gen = net.get_generators().sort_index() else: df_gen = net.get_generators() + # to handle encoding in 32 bits and overflow when "splitting" the Q values among - min_q = df_gen["min_q"].values.astype(np.float32) - max_q = df_gen["max_q"].values.astype(np.float32) - min_q[~np.isfinite(min_q)] = np.finfo(np.float32).min * 1e-4 + 1. - max_q[~np.isfinite(max_q)] = np.finfo(np.float32).max * 1e-4 - 1. + min_float_value = np.finfo(np.float32).min * 1e-4 + 1. + max_float_value = np.finfo(np.float32).max * 1e-4 + 1. + min_q_aux = 1. * df_gen["min_q"].values + too_small = min_q_aux < min_float_value + min_q_aux[too_small] = min_float_value + min_q = min_q_aux.astype(np.float32) + + max_q_aux = 1. * df_gen["max_q"].values + too_big = np.abs(max_q_aux) > max_float_value + max_q_aux[too_big] = np.sign(max_q_aux[too_big]) * max_float_value + max_q = max_q_aux.astype(np.float32) + min_q[~np.isfinite(min_q)] = min_float_value + max_q[~np.isfinite(max_q)] = max_float_value gen_bus, gen_disco = _aux_get_bus(bus_df, df_gen) + + # dirty fix for when regulating elements are not the same + bus_reg = copy.deepcopy(df_gen["regulated_element_id"].values) + vl_reg = copy.deepcopy(df_gen["voltage_level_id"].values) + mask_ref_bbs = bus_reg != df_gen.index + + bbs_df = net.get_busbar_sections().copy() + if not (np.isin(bus_reg[mask_ref_bbs], bbs_df.index)).all(): + raise RuntimeError("At least some generator are in 'remote control' mode " + "and does not control a busbar section, this is not supported " + "at the moment.") + vl_reg[mask_ref_bbs] = bbs_df.loc[bus_reg[mask_ref_bbs], "voltage_level_id"].values model.init_generators_full(df_gen["target_p"].values, - df_gen["target_v"].values / voltage_levels.loc[df_gen["voltage_level_id"].values]["nominal_v"].values, + # df_gen["target_v"].values / voltage_levels.loc[df_gen["voltage_level_id"].values]["nominal_v"].values, + df_gen["target_v"].values / voltage_levels.loc[vl_reg]["nominal_v"].values, df_gen["target_q"].values, df_gen["voltage_regulator_on"].values, min_q, @@ -116,30 +198,44 @@ def init(net : pypo.network, else: df_line = net.get_lines() # per unit - branch_from_kv = voltage_levels.loc[df_line["voltage_level1_id"].values]["nominal_v"].values - branch_to_kv = voltage_levels.loc[df_line["voltage_level2_id"].values]["nominal_v"].values + if net_pu is None: + try: + net_pu = copy.deepcopy(net) + net_pu.per_unit = True + except Exception as exc_: + from pypowsybl.network import PerUnitView + net_pu = PerUnitView(net) + warnings.warn("The `PerUnitView` (python side) is less efficient and less " + "tested that the equivalent java class. Please upgrade pypowsybl version") + df_line_pu = net_pu.get_lines().loc[df_line.index] + # branch_from_kv = voltage_levels.loc[df_line["voltage_level1_id"].values]["nominal_v"].values + # branch_to_kv = voltage_levels.loc[df_line["voltage_level2_id"].values]["nominal_v"].values # only valid for lines with same voltages at both side... - # branch_from_pu = branch_from_kv * branch_from_kv / sn_mva - # line_r = df_line["r"].values / branch_from_pu - # line_x = df_line["x"].values / branch_from_pu - # line_h_or = (1j*df_line["g1"].values + df_line["b1"].values) * branch_from_pu - # line_h_ex = (1j*df_line["g2"].values + df_line["b2"].values) * branch_from_pu + # # branch_from_pu = branch_from_kv * branch_from_kv / sn_mva + # # line_r = df_line["r"].values / branch_from_pu + # # line_x = df_line["x"].values / branch_from_pu + # # line_h_or = (1j*df_line["g1"].values + df_line["b1"].values) * branch_from_pu + # # line_h_ex = (1j*df_line["g2"].values + df_line["b2"].values) * branch_from_pu # real per unit conversion # see https://github.com/powsybl/pypowsybl/issues/642 # see https://github.com/powsybl/powsybl-core/blob/266442cbbd84f630acf786018618eaa3d496c6ba/ieee-cdf/ieee-cdf-converter/src/main/java/com/powsybl/ieeecdf/converter/IeeeCdfImporter.java#L347 # for right formula - v1 = branch_from_kv - v2 = branch_to_kv - line_r = sn_mva * df_line["r"].values / v1 / v2 - line_x = sn_mva * df_line["x"].values / v1 / v2 - tmp_ = np.reciprocal(df_line["r"].values + 1j*df_line["x"].values) - b1 = df_line["b1"].values * v1*v1/sn_mva + (v1-v2)*tmp_.imag*v1/sn_mva - b2 = df_line["b2"].values * v2*v2/sn_mva + (v2-v1)*tmp_.imag*v2/sn_mva - g1 = df_line["g1"].values * v1*v1/sn_mva + (v1-v2)*tmp_.real*v1/sn_mva - g2 = df_line["g2"].values * v2*v2/sn_mva + (v2-v1)*tmp_.real*v2/sn_mva - line_h_or = (b1 + 1j * g1) - line_h_ex = (b2 + 1j * g2) + # v1 = branch_from_kv + # v2 = branch_to_kv + # line_r = sn_mva * df_line["r"].values / v1 / v2 + # line_x = sn_mva * df_line["x"].values / v1 / v2 + # tmp_ = np.reciprocal(df_line["r"].values + 1j*df_line["x"].values) + # b1 = df_line["b1"].values * v1*v1/sn_mva + (v1-v2)*tmp_.imag*v1/sn_mva + # b2 = df_line["b2"].values * v2*v2/sn_mva + (v2-v1)*tmp_.imag*v2/sn_mva + # g1 = df_line["g1"].values * v1*v1/sn_mva + (v1-v2)*tmp_.real*v1/sn_mva + # g2 = df_line["g2"].values * v2*v2/sn_mva + (v2-v1)*tmp_.real*v2/sn_mva + # line_h_or = (g1 + 1j * b1) + # line_h_ex = (g2 + 1j * b2) + line_r = df_line_pu["r"].values + line_x = df_line_pu["x"].values + line_h_or = (df_line_pu["g1"].values + 1j * df_line_pu["b1"].values) + line_h_ex = (df_line_pu["g2"].values + 1j * df_line_pu["b2"].values) lor_bus, lor_disco = _aux_get_bus(bus_df, df_line, conn_key="connected1", bus_key="bus1_id") lex_bus, lex_disco = _aux_get_bus(bus_df, df_line, conn_key="connected2", bus_key="bus2_id") model.init_powerlines_full(line_r, @@ -159,28 +255,66 @@ def init(net : pypo.network, df_trafo = net.get_2_windings_transformers().sort_index() else: df_trafo = net.get_2_windings_transformers() - # TODO net.get_ratio_tap_changers() - # TODO net.get_phase_tap_changers() + + df_trafo_pu = net_pu.get_2_windings_transformers().loc[df_trafo.index] + ratio_tap_changer = net_pu.get_ratio_tap_changers() + if net.get_phase_tap_changers().shape[0] > 0: + raise RuntimeError("Phase tap changer are not handled by the pypowsybl converter (but they are by lightsim2grid)") + # phase_tap_changer = net_pu.get_phase_tap_changers() + shift_ = np.zeros(df_trafo.shape[0]) - tap_pos = 1.0 * shift_ - is_tap_hv_side = np.ones(df_trafo.shape[0], dtype=bool) # TODO + tap_position = 1.0 * shift_ + is_tap_hv_side = np.zeros(df_trafo.shape[0], dtype=bool) # TODO # per unit - trafo_from_kv = voltage_levels.loc[df_trafo["voltage_level1_id"].values]["nominal_v"].values - trafo_to_kv = voltage_levels.loc[df_trafo["voltage_level2_id"].values]["nominal_v"].values - trafo_to_pu = trafo_to_kv * trafo_to_kv / sn_mva + # trafo_from_kv = voltage_levels.loc[df_trafo["voltage_level1_id"].values]["nominal_v"].values + # trafo_to_kv = voltage_levels.loc[df_trafo["voltage_level2_id"].values]["nominal_v"].values + # trafo_to_pu = trafo_to_kv * trafo_to_kv / sn_mva # tap - tap_step_pct = (df_trafo["rated_u1"] / trafo_from_kv - 1.) * 100. - has_tap = tap_step_pct != 0. - tap_pos[has_tap] += 1 - tap_step_pct[~has_tap] = 1.0 # or any other values... + # tap_step_pct = (df_trafo["rated_u1"] / trafo_from_kv - 1.) * 100. + # has_tap = tap_step_pct != 0. + # tap_pos[has_tap] += 1 + # tap_step_pct[~has_tap] = 1.0 # or any other values... + # trafo_r = df_trafo["r"].values / trafo_to_pu + # trafo_x = df_trafo["x"].values / trafo_to_pu + # trafo_h = (df_trafo["g"].values + 1j * df_trafo["b"].values) * trafo_to_pu + trafo_r = df_trafo_pu["r"].values + trafo_x = df_trafo_pu["x"].values + trafo_h = (df_trafo_pu["g"].values + 1j * df_trafo_pu["b"].values) + + # now get the ratio + # in lightsim2grid (cpp) + # RealVect ratio = my_one_ + 0.01 * trafo_tap_step_pct.array() * trafo_tap_pos.array(); + # in powsybl (https://javadoc.io/doc/com.powsybl/powsybl-core/latest/com/powsybl/iidm/network/TwoWindingsTransformer.html) + # rho = transfo.getRatedU2() / transfo.getRatedU1() + # * (transfo.getRatioTapChanger() != null ? transfo.getRatioTapChanger().getCurrentStep().getRho() : 1); + # * (transfo.getPhaseTapChanger() != null ? transfo.getPhaseTapChanger().getCurrentStep().getRho() : 1); + ratio = 1. * (df_trafo_pu["rated_u2"].values / df_trafo_pu["rated_u1"].values) + has_r_tap_changer = np.isin(df_trafo_pu.index, ratio_tap_changer.index) + + if PYPOWSYBL_VER <= PP_BUG_RATIO_TAP_CHANGER: + # bug in per unit view in both python and java + ratio[has_r_tap_changer] = 1. * ratio_tap_changer.loc[df_trafo_pu.loc[has_r_tap_changer].index, "rho"].values + else: + ratio[has_r_tap_changer] *= 1. * ratio_tap_changer.loc[df_trafo_pu.loc[has_r_tap_changer].index, "rho"].values + no_tap = ratio == 1. + tap_neg = ratio < 1. + tap_positive = ratio > 1. + tap_step_pct = 1. * ratio + tap_step_pct[tap_positive] -= 1. + tap_step_pct[tap_positive] *= 100. + tap_step_pct[tap_neg] = (ratio[tap_neg] - 1.)*100. + tap_step_pct[no_tap] = 100. + tap_position[tap_positive] += 1 + tap_position[tap_neg] += 1 + tor_bus, tor_disco = _aux_get_bus(bus_df, df_trafo, conn_key="connected1", bus_key="bus1_id") tex_bus, tex_disco = _aux_get_bus(bus_df, df_trafo, conn_key="connected2", bus_key="bus2_id") - model.init_trafo(df_trafo["r"].values / trafo_to_pu, - df_trafo["x"].values / trafo_to_pu, - 2.*(1j*df_trafo["g"].values + df_trafo["b"].values) * trafo_to_pu, + model.init_trafo(trafo_r, + trafo_x, + trafo_h, tap_step_pct, - tap_pos, + tap_position, shift_, is_tap_hv_side, tor_bus, # TODO do I need to change hv / lv @@ -267,8 +401,8 @@ def init(net : pypo.network, else: try: gen_slack_id_int = int(gen_slack_id) - except Exception: - raise RuntimeError("'slack_bus_id' should be either an int or a generator names") + except Exception as exc_: + raise RuntimeError("'slack_bus_id' should be either an int or a generator names") from exc_ if gen_slack_id_int != gen_slack_id: raise RuntimeError("'slack_bus_id' should be either an int or a generator names") model.add_gen_slackbus(gen_slack_id_int, 1.) @@ -281,6 +415,8 @@ def init(net : pypo.network, for gen_id, is_slack in enumerate(gen_is_conn_slack): if is_slack: model.add_gen_slackbus(gen_id, 1. / nb_conn) + else: + raise RuntimeError("You need to provide at least one slack with `gen_slack_id` or `slack_bus_id`") # TODO # sgen => regular gen (from net.get_generators()) with voltage_regulator off TODO @@ -296,13 +432,19 @@ def init(net : pypo.network, # and now deactivate all elements and nodes not in the main component if only_main_component: model.consider_only_main_component() + model.init_bus_status() # automatically disconnect non connected buses if not return_sub_id: # for backward compatibility return model else: # voltage_level_id is kind of what I call "substation" in grid2op - vl_unique = bus_df["voltage_level_id"].unique() - sub_df = pd.DataFrame(index=np.sort(vl_unique), data={"sub_id": np.arange(vl_unique.size)}) + if sub_unique is None: + vl_unique = bus_df["voltage_level_id"].unique() + sub_df = pd.DataFrame(index=np.sort(vl_unique), data={"sub_id": np.arange(vl_unique.size)}) + else: + # already computed before when innitializing the buses and substations + # I don't do it again to ensure consistency + sub_df = pd.DataFrame(index=sub_unique, data={"sub_id": sub_unique_id}) buses_sub_id = pd.merge(left=bus_df, right=sub_df, how="left", left_on="voltage_level_id", right_index=True)[["bus_id", "sub_id"]] gen_sub = pd.merge(left=df_gen, right=sub_df, how="left", left_on="voltage_level_id", right_index=True)[["sub_id"]] load_sub = pd.merge(left=df_load, right=sub_df, how="left", left_on="voltage_level_id", right_index=True)[["sub_id"]] diff --git a/lightsim2grid/lightSimBackend.py b/lightsim2grid/lightSimBackend.py index 5738f60..35e91e7 100644 --- a/lightsim2grid/lightSimBackend.py +++ b/lightsim2grid/lightSimBackend.py @@ -24,7 +24,7 @@ import grid2op from grid2op.Action import CompleteAction from grid2op.Backend import Backend -from grid2op.Exceptions import BackendError +from grid2op.Exceptions import BackendError, Grid2OpException from grid2op.dtypes import dt_float, dt_int, dt_bool try: from grid2op.Action._backendAction import _BackendAction @@ -72,32 +72,66 @@ def __init__(self, stop_if_load_disco : Optional[bool] = True, stop_if_gen_disco : Optional[bool] = True, ): + #: ``int`` maximum number of iteration allowed for the solver + #: if the solver has not converge after this, it will + #: send a "divergence error" self.max_it = max_iter + + #: ``float`` tolerance of the solver self.tol = tol # tolerance for the solver + self._check_suitable_solver_type(solver_type, check_in_avail_solver=False) self.__current_solver_type = solver_type - # does the "turned off" generators (including when p=0) - # are pv buses + #: does the "turned off" generators (including when p=0) + #: are pv buses self._turned_off_pv = turned_off_pv - # distributed slack, on non renewable gen with P > 0 + #: distributed slack, on non renewable gen with P > 0 self._dist_slack_non_renew = dist_slack_non_renew # add the static gen to the list of controlable gen in grid2Op self._use_static_gen = use_static_gen # TODO implement it + #: For now, you can initialize a "lightsim2grid" LightsimBackend + #: either from pypowsybl or from pandapower. + #: Use with `LightsimBackend(..., loader_method="pypowsybl")` + #: or `LightsimBackend(..., loader_method="pandapower")` (default) self._loader_method = loader_method + #: Which key-word arguments will be used to initialize the Gridmodel + #: either from pandapower or lightsim2grid. + #: It is not currently used for pandapower. + #: + #: For pypowsybl it can contain the following keys: + #: + #: - `n_busbar_per_sub` (``int``): number of independant buses for + #: each substation in the GridModel. + #: - `use_buses_for_sub` (``bool``): whether to use the buses (in the + #: pypowsybl Network) to initialize the "substation" in the lightsim2grid + #: Gridmodel (if ``True``). If ``False`` it will use the `voltage_levels` + #: of the pypowsybl Network. + #: - `gen_slack_id` (``int``): which generator will be used for the slack + #: - `use_grid2op_default_names` (``bool``): whether to use the "default names" + #: assigne by grid2op or to read them from the the iidm grid. + #: - `reconnect_disco_gen` (``bool``): whether to reconnect the disconnected + #: generators in the iidm grid. If set to ``True`` then the generators will be + #: reconnected with a p setpoint of 0. MW and a voltage setpoint of 1 pu + #: - `reconnect_disco_load` (``bool``): whether to reconnec the disconnected + #: load from in the iidm grid. If reconnected, load will have a target + #: p set to 0. and a target q set to 0. + #: self._loader_kwargs = loader_kwargs #: .. versionadded:: 0.8.0 - #: if set to `True` (default) then the backend will raise a + #: + #: if set to ``True`` (default) then the backend will raise a #: BackendError in case of disconnected load self._stop_if_load_disco = stop_if_load_disco #: .. versionadded:: 0.8.0 - #: if set to `True` (default) then the backend will raise a + #: + #: if set to ``True`` (default) then the backend will raise a #: BackendError in case of disconnected generator self._stop_if_gen_disco = stop_if_gen_disco @@ -119,10 +153,12 @@ def __init__(self, self._can_be_copied = can_be_copied #: .. versionadded:: 0.8.0 + #: #: Which type of grid format can be read by your backend. #: It is "json" if loaded from pandapower or #: "xiidm" if loaded from pypowsybl. self.supported_grid_format = None + if loader_method == "pandapower": self.supported_grid_format = ("json", ) # new in 0.8.0 elif loader_method == "pypowsybl": @@ -135,7 +171,7 @@ def __init__(self, from grid2op.Space import GridObjects # lazy import self.__has_storage = hasattr(GridObjects, "n_storage") if not self.__has_storage: - warnings.warn("Please upgrade your grid2Op to >= 1.5.0. You are using a backward compatibility " + warnings.warn("Please upgrade your grid2op to >= 1.5.0. You are using a backward compatibility " "feature that will be removed in further lightsim2grid version.") if version.parse(grid2op.__version__) < grid2op_min_shunt_cls_properly_handled: warnings.warn(f"You are using a legacy grid2op version. It is not possible to deactivate the shunts in lightsim2grid. " @@ -266,6 +302,7 @@ def __init__(self, self._storage_res = None self._reset_res_pointers() self._debug_Vdc = None # use only for debug ! + self._orig_grid_pypowsybl = None def _aux_init_super(self, detailed_infos_for_cascading_failures, @@ -507,8 +544,6 @@ def _should_not_have_to_do_this(self, path=None, filename=None): # included in grid2op now ! # but before `make_complete_path` was introduced we need to still # be able to use lightsim2grid - import os - from grid2op.Exceptions import Grid2OpException if path is None and filename is None: raise Grid2OpException( "You must provide at least one of path or file to load a powergrid." @@ -523,25 +558,41 @@ def _should_not_have_to_do_this(self, path=None, filename=None): raise Grid2OpException('There is no powergrid at "{}"'.format(full_path)) return full_path - def _aux_pypowsybl_init_substations(self, loader_kwargs): - # now handle the legacy "make as if there are 2 busbars per substation" - # as it is done with grid2Op simulated environment - if (("double_bus_per_sub" in loader_kwargs and loader_kwargs["double_bus_per_sub"]) or - ("n_busbar_per_sub" in loader_kwargs and loader_kwargs["n_busbar_per_sub"])): - bus_init = self._grid.get_bus_vn_kv() - orig_to_ls = np.array(self._grid._orig_to_ls) - bus_doubled = np.concatenate([bus_init for _ in range(self.n_busbar_per_sub)]) - self._grid.init_bus(bus_doubled, 0, 0) - if hasattr(type(self), "can_handle_more_than_2_busbar"): - for i in range(self.__nb_bus_before * (self.n_busbar_per_sub - 1)): - self._grid.deactivate_bus(i + self.__nb_bus_before) - else: - for i in range(self.__nb_bus_before): - self._grid.deactivate_bus(i + self.__nb_bus_before) - new_orig_to_ls = np.concatenate([orig_to_ls + i * self.__nb_bus_before - for i in range(self.n_busbar_per_sub)] - ) - self._grid._orig_to_ls = new_orig_to_ls + def _aux_get_substation_handling_from_loader_kwargs(self, loader_kwargs) -> Union[int, None]: + res = None + if "double_bus_per_sub" in loader_kwargs and loader_kwargs["double_bus_per_sub"]: + res = DEFAULT_N_BUSBAR_PER_SUB + if "n_busbar_per_sub" in loader_kwargs and loader_kwargs["n_busbar_per_sub"]: + if res is not None: + raise BackendError("When intializing a grid from pypowsybl, you cannot " + "set both `double_bus_per_sub` and `n_busbar_per_sub` " + "in the `loader_kwargs`. " + "You can only set `n_busbar_per_sub` in this case.") + res = int(loader_kwargs["n_busbar_per_sub"]) + if loader_kwargs["n_busbar_per_sub"] != res: + raise BackendError("When initializing a grid from pypowsybl, the `n_busbar_per_sub` " + "loader kwargs should be properly convertible to an int " + "giving the default number of busbar sections per substation.") + if res is not None: + if self.n_busbar_per_sub != res: + if self.n_busbar_per_sub != DEFAULT_N_BUSBAR_PER_SUB: + warnings.warn(f"You specified `n_busbar={self.n_busbar_per_sub}` when calling `grid2op.make(...)` " + f"but also specified `n_busbar_per_sub` (or `double_bus_per_sub`) in the " + f"`loader_kwargs` when building the `LightSimBackend` " + f"(**eg** `LighSimBackend(..., loader_kwargs={{'n_busbar_per_sub': {res} }})`). " + f"IN THIS CASE {self.n_busbar_per_sub} WILL BE USED (the things specified in the " + f"`grid2op.make(...)`)") + res = self.n_busbar_per_sub + elif res != DEFAULT_N_BUSBAR_PER_SUB: + warnings.warn(f"You specified `n_busbar_per_sub={res}` (or `double_bus_per_sub`) " + f"in the `load_kwargs` when building the LightSimBackend " + f"(**eg** `LighSimBackend(..., loader_kwargs={{'n_busbar_per_sub': {res} }})`) " + f"without modifying the `n_busbar=...` when calling `grid2op.make(...)`. " + f"LightSimBackend WILL OVERRIDE THE MAXIMUM NUMBER OF BUSES PER SUBSTATION " + f"and set it to {res}. You can silence this warning by creating the env with " + f"`grid2op.make(..., n_busbar={res})`") + self.n_busbar_per_sub = res + return res def _get_subid_from_buses_legacy(self, buses_sub_id, el_sub_df): # buses_sub_id is the first element as returned by from_pypowsybl / init function @@ -564,14 +615,35 @@ def _load_grid_pypowsybl(self, path=None, filename=None): full_path = self._should_not_have_to_do_this(path, filename) grid_tmp = pypow_net.load(full_path) + self._orig_grid_pypowsybl = grid_tmp gen_slack_id = None if "gen_slack_id" in loader_kwargs: gen_slack_id = loader_kwargs["gen_slack_id"] - self._grid, subs_id = init_pypow(grid_tmp, gen_slack_id=gen_slack_id, sort_index=True, return_sub_id=True) + + df = grid_tmp.get_substations() + buses_for_sub = True + if "use_buses_for_sub" in loader_kwargs and loader_kwargs["use_buses_for_sub"]: + df = grid_tmp.get_buses() + buses_for_sub = False + self.n_sub = df.shape[0] + self.name_sub = ["sub_{}".format(i) for i, _ in enumerate(df.iterrows())] + + n_busbar_per_sub = self._aux_get_substation_handling_from_loader_kwargs(loader_kwargs) + if n_busbar_per_sub is None: + n_busbar_per_sub = self.n_busbar_per_sub + self._grid, subs_id = init_pypow(grid_tmp, + gen_slack_id=gen_slack_id, + sort_index=True, + return_sub_id=True, + n_busbar_per_sub=n_busbar_per_sub, + buses_for_sub=buses_for_sub, + ) (buses_sub_id, gen_sub, load_sub, (lor_sub, tor_sub), (lex_sub, tex_sub), batt_sub, sh_sub, hvdc_sub_from_id, hvdc_sub_to_id) = subs_id - self.__nb_bus_before = len(self._grid.get_bus_vn_kv()) - self._aux_pypowsybl_init_substations(loader_kwargs) + if not buses_for_sub: + self.__nb_bus_before = self._grid.get_n_sub() + else: + self.__nb_bus_before = grid_tmp.get_buses().shape[0] self._aux_setup_right_after_grid_init() # mandatory for the backend @@ -584,16 +656,7 @@ def _load_grid_pypowsybl(self, path=None, filename=None): else: self.n_shunt = None - df = grid_tmp.get_substations() - from_sub = True - - if "use_buses_for_sub" in loader_kwargs and loader_kwargs["use_buses_for_sub"]: - df = grid_tmp.get_buses() - from_sub = False - self.n_sub = df.shape[0] - self.name_sub = ["sub_{}".format(i) for i, _ in enumerate(df.iterrows())] - - if not from_sub: + if not buses_for_sub: # consider that each "bus" in the powsybl grid is a substation # this is the "standard" behaviour for IEEE grid in grid2op # but can be considered "legacy" behaviour for more realistic grid @@ -613,14 +676,19 @@ def _load_grid_pypowsybl(self, path=None, filename=None): self.storage_to_subid = np.array(this_batt_sub, dtype=dt_int) self.shunt_to_subid = np.array(this_sh_sub, dtype=dt_int) else: - # TODO get back the sub id from the grid_tmp.get_substations() + # consider effectively that each "voltage_levels" in powsybl grid + # is a substation (in the underlying gridmodel) + + # TODO get back the sub id from the grid_tmp.get_voltage_levels() # need to work on that grid2op side: different make sure the labelling of the buses are correct ! - (buses_sub_id, gen_sub, load_sub, (lor_sub, tor_sub), (lex_sub, tex_sub), batt_sub, sh_sub, hvdc_sub_from_id, hvdc_sub_to_id) = subs_id - raise NotImplementedError("Today the only supported behaviour is to consider the 'buses' of the powsybl grid " - "are the 'substation' of this backend. " - "This will change in the future, but in the meantime please add " - "'use_buses_for_sub' = True in the `loader_kwargs` when loading " - "a lightsim2grid grid.") + + self.load_to_subid = np.array(load_sub.values.ravel(), dtype=dt_int) + self.gen_to_subid = np.array(gen_sub.values.ravel(), dtype=dt_int) + self.line_or_to_subid = np.concatenate((lor_sub.values.ravel(), tor_sub.values.ravel())).astype(dt_int) + self.line_ex_to_subid = np.concatenate((lex_sub.values.ravel(), tex_sub.values.ravel())).astype(dt_int) + self.storage_to_subid = np.array(batt_sub.values.ravel(), dtype=dt_int) + self.shunt_to_subid = np.array(sh_sub.values.ravel(), dtype=dt_int) + self.n_sub = grid_tmp.get_voltage_levels().shape[0] # the names use_grid2op_default_names = True @@ -628,31 +696,42 @@ def _load_grid_pypowsybl(self, path=None, filename=None): use_grid2op_default_names = False if use_grid2op_default_names: - self.name_load = np.array([f"load_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_loads())]) - self.name_gen = np.array([f"gen_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_generators())]) + self.name_load = np.array([f"load_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_loads())]).astype(str) + self.name_gen = np.array([f"gen_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_generators())]).astype(str) self.name_line = np.array([f"{el.bus_or_id}_{el.bus_ex_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_lines())] + - [f"{el.bus_hv_id}_{el.bus_lv_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_trafos())]) - self.name_storage = np.array([f"storage_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_storages())]) - self.name_shunt = np.array([f"shunt_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_shunts())]) + [f"{el.bus_hv_id}_{el.bus_lv_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_trafos())]).astype(str) + self.name_storage = np.array([f"storage_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_storages())]).astype(str) + self.name_shunt = np.array([f"shunt_{el.bus_id}_{id_obj}" for id_obj, el in enumerate(self._grid.get_shunts())]).astype(str) else: - self.name_load = np.array(load_sub.index) - self.name_gen = np.array(gen_sub.index) - self.name_line = np.concatenate((lor_sub.index, tor_sub.index)) - self.name_storage = np.array(batt_sub.index) - self.name_shunt = np.array(sh_sub.index) - self.name_sub = np.array(buses_sub_id.index) + self.name_load = np.array(load_sub.index.astype(str)) + self.name_gen = np.array(gen_sub.index.astype(str)) + self.name_line = np.concatenate((lor_sub.index.astype(str), tor_sub.index.astype(str))) + self.name_storage = np.array(batt_sub.index.astype(str)) + self.name_shunt = np.array(sh_sub.index.astype(str)) + if not buses_for_sub: + self.name_sub = np.array(buses_sub_id.index.astype(str)) + else: + self.name_sub = np.array(grid_tmp.get_buses().loc[buses_sub_id.drop_duplicates("sub_id").index, "voltage_level_id"].values) + + # and now things needed by the backend (legacy) + self.prod_pu_to_kv = 1.0 * self._grid.get_bus_vn_kv()[[el.bus_id for el in self._grid.get_generators()]] + self.prod_pu_to_kv = self.prod_pu_to_kv.astype(dt_float) if "reconnect_disco_gen" in loader_kwargs and loader_kwargs["reconnect_disco_gen"]: for el in self._grid.get_generators(): if not el.connected: self._grid.reactivate_gen(el.id) self._grid.change_bus_gen(el.id, self.gen_to_subid[el.id]) + self._grid.change_p_gen(el.id, 0.) + self._grid.change_v_gen(el.id, self.prod_pu_to_kv[el.id]) if "reconnect_disco_load" in loader_kwargs and loader_kwargs["reconnect_disco_load"]: for el in self._grid.get_loads(): if not el.connected: self._grid.reactivate_load(el.id) self._grid.change_bus_load(el.id, self.load_to_subid[el.id]) + self._grid.change_p_load(el.id, 0.) + self._grid.change_q_load(el.id, 0.) # complete the other vectors self._compute_pos_big_topo() @@ -678,7 +757,8 @@ def _load_grid_pypowsybl(self, path=None, filename=None): self._aux_finish_setup_after_reading() def _aux_setup_right_after_grid_init(self): - self._grid.set_n_sub(self.__nb_bus_before) + if self._orig_grid_pypowsybl is None: + self._grid.set_n_sub(self.__nb_bus_before) self._handle_turnedoff_pv() self.available_solvers = self._grid.available_solvers() @@ -699,7 +779,7 @@ def _aux_setup_right_after_grid_init(self): self._grid.change_solver(self.__current_solver_type) # handle multiple busbar per substations - if hasattr(type(self), "can_handle_more_than_2_busbar"): + if hasattr(type(self), "can_handle_more_than_2_busbar") and self._orig_grid_pypowsybl is None: # grid2op version >= 1.10.0 then we use this self._grid._max_nb_bus_per_sub = self.n_busbar_per_sub @@ -1056,7 +1136,7 @@ def apply_action(self, backendAction: Union["grid2op.Action._backendAction._Back # print(f" load p {backendAction.load_p.values[backendAction.load_p.changed]}") # TODO DEBUG WINDOWS # print(f" prod_p p {backendAction.prod_p.values[backendAction.prod_p.changed]}") # TODO DEBUG WINDOWS - # update the injections + # update the injections try: self._grid.update_gens_p(backendAction.prod_p.changed, backendAction.prod_p.values) @@ -1434,6 +1514,7 @@ def copy(self) -> Self: res._backend_action_class = self._backend_action_class # this is const res.__init_topo_vect = self.__init_topo_vect res.available_solvers = self.available_solvers + res._orig_grid_pypowsybl = self._orig_grid_pypowsybl res._reset_res_pointers() res._fetch_grid_data() # assign back "self" attributes diff --git a/lightsim2grid/tests/pandapower_contingency.py b/lightsim2grid/tests/pandapower_contingency.py new file mode 100644 index 0000000..49a5891 --- /dev/null +++ b/lightsim2grid/tests/pandapower_contingency.py @@ -0,0 +1,526 @@ +# -*- coding: utf-8 -*- + +# Copyright (c) 2016-2023 by University of Kassel and Fraunhofer Institute for Energy Economics +# and Energy System Technology (IEE), Kassel. All rights reserved. +# copied on August, 27th, 2024 from + +import numpy as np +import pandas as pd +import warnings + +import pandapower as pp + +############################# +# way to import lightsim2grid modified to work correctly independantly +# of lightsim2grid version +# see https://github.com/e2nIEE/pandapower/issues/2361#issuecomment-2263091364 +# for more information + +try: + # lightsim2grid version >= 0.9.0 + from lightsim2grid.gridmodel import init_from_pandapower as init_ls2g + lightsim2grid_installed = True +except ImportError: + try: + # lightsim2grid <= 0.8.2 + from lightsim2grid.gridmodel import init as int_ls2g + lightsim2grid_installed = True + except ImportError: + # lightsim2grid is most likely not installed + lightsim2grid_installed = False + +if lightsim2grid_installed: + try: + # ligthsim2grid >= 0.8.0 + from lightsim2grid.contingencyAnalysis import ContingencyAnalysisCPP + except ImportError: + try: + # ligthsim2grid < 0.8.0 + from lightsim2igrd.securityAnalysis import SecurityAnalysisCPP as ContingencyAnalysisCPP + except ImportError: + # issue a warning here maybe ? Most likely lightsim2grid version is too old + lightsim2grid_installed = False + try: + from lightsim2grid.solver import SolverType + except ImportError: + # issue a warning here maybe ? Most likely lightsim2grid version is too old + lightsim2grid_installed = False +# end of the modifications +########################### + +try: + from lightsim2grid_cpp import KLUSolver, KLUSolverSingleSlack + + KLU_solver_available = True +except ImportError: + KLU_solver_available = False + +try: + import pandaplan.core.pplog as logging +except ImportError: + import logging + +logger = logging.getLogger(__name__) + + +def run_contingency(net, nminus1_cases, pf_options=None, pf_options_nminus1=None, write_to_net=True, + contingency_evaluation_function=pp.runpp, **kwargs): + """ + Obtain either loading (N-0) or max. loading (N-0 and all N-1 cases), and min/max bus voltage magnitude. + The variable "temperature_degree_celsius" can be used in addition to "loading_percent" to obtain max. temperature. + In the returned dictionary, the variable "loading_percent" represents the loading in N-0 case, + "max_loading_percent" and "min_loading_percent" represent highest and lowest observed loading_percent among all + calculated N-1 cases. The same convention applies to "temperature_degree_celsius" when applicable. + This function can be passed through to pandapower.timeseries.run_timeseries as the run_control_fct argument. + + INPUT + ---------- + **net** - pandapowerNet + **nminus1_cases** - dict + describes all N-1 cases, e.g. {"line": {"index": [1, 2, 3]}, "trafo": {"index": [0]}, "trafo3w": {"index": [1]}} + **pf_options** - dict + options for power flow calculation in N-0 case + **pf_options_nminus1** - dict + options for power flow calculation in N-1 cases + **write_to_net** - bool + whether to write the results of contingency analysis to net (in "res_" tables). The results will be written for + the following additional variables: table res_bus with columns "max_vm_pu", "min_vm_pu", + tables res_line, res_trafo, res_trafo3w with columns "max_loading_percent", "min_loading_percent", + "causes_overloading", "cause_element", "cause_index", table res_line with columns + "max_temperature_degree_celsius", "min_temperature_degree_celsius" (if "tdpf" set to True) + "causes_overloading": does this element, when defining the N-1 case, cause overloading of other elements? the + overloading is defined by net.line["max_loading_percent_nminus1"] (if set) or net.line["max_loading_percent"] + "cause_element": element ("line", "trafo", "trafo3w") that causes max. loading of this element + "cause_index": index of the element ("line", "trafo", "trafo3w") that causes max. loading of this element + **contingency_evaluation_function** - func + function to use for power flow calculation, default pp.runpp + + OUTPUT + ------- + **contingency_results** - dict + dict of arrays per element for index, min/max result + """ + # set up the dict for results and relevant variables + # ".get" in case the options have been set in pp.set_user_pf_options: + raise_errors = kwargs.get("raise_errors", False) + if "recycle" in kwargs: kwargs["recycle"] = False # so that we can be sure it doesn't happen + if pf_options is None: pf_options = net.user_pf_options.get("pf_options", net.user_pf_options) + if pf_options_nminus1 is None: pf_options_nminus1 = net.user_pf_options.get("pf_options_nminus1", + net.user_pf_options) + + contingency_results = {element: {"index": net[element].index.values} + for element in ("bus", "line", "trafo", "trafo3w") if len(net[element]) > 0} + for element in contingency_results.keys(): + if element == "bus": + continue + contingency_results[element].update( + {"causes_overloading": np.zeros_like(net[element].index.values, dtype=bool), + "cause_element": np.empty_like(net[element].index.values, dtype=object), + "cause_index": np.empty_like(net[element].index.values, dtype=np.int64)}) + result_variables = {**{"bus": ["vm_pu"]}, + **{key: ["loading_percent"] for key in ("line", "trafo", "trafo3w") if len(net[key]) > 0}} + if len(net.line) > 0 and (net.get("_options", {}).get("tdpf", False) or + pf_options.get("tdpf", False) or pf_options_nminus1.get("tdpf", False)): + result_variables["line"].append("temperature_degree_celsius") + + # for n-1 + for element, val in nminus1_cases.items(): + for i in val["index"]: + if ~net[element].at[i, "in_service"]: + continue + net[element].at[i, 'in_service'] = False + try: + contingency_evaluation_function(net, **pf_options_nminus1, **kwargs) + _update_contingency_results(net, contingency_results, result_variables, nminus1=True, + cause_element=element, cause_index=i) + except Exception as err: + logger.error(f"{element} {i} causes {err}") + if raise_errors: + raise err + finally: + net[element].at[i, 'in_service'] = True + + # for n-0 + contingency_evaluation_function(net, **pf_options, **kwargs) + _update_contingency_results(net, contingency_results, result_variables, nminus1=False) + + if write_to_net: + for element, element_results in contingency_results.items(): + index = element_results["index"] + for var, val in element_results.items(): + if var == "index" or var in net[f"res_{element}"].columns.values: + continue + net[f"res_{element}"].loc[index, var] = val + + return contingency_results + + +def run_contingency_ls2g(net, nminus1_cases, contingency_evaluation_function=pp.runpp, **kwargs): + """ + Execute contingency analysis using the lightsim2grid library. This works much faster than using pandapower. + Limitation: the results for branch flows are valid only for the "from_bus" of lines and "hv_bus" of transformers. + This can lead to a small difference to the results using pandapower. + The results are written in pandapower results tables. + Make sure that the N-1 cases do not lead to isolated grid, otherwise results with pandapower and this function will + be different. Reason: pandapower selects a different gen as slack if the grid becomes isolated, but + lightsim2grid would simply return nan as results for such a contingency situation. + WARNING: continuous bus indices, 0-start, are required! + This function can be passed through to pandapower.timeseries.run_timeseries as the run_control_fct argument. + + The results will written for the + following additional variables: table res_bus with columns "max_vm_pu", "min_vm_pu", + tables res_line, res_trafo, res_trafo3w with columns "max_loading_percent", "min_loading_percent", + "causes_overloading", "cause_element", "cause_index", table res_line with columns "max_temperature_degree_celsius", + "min_temperature_degree_celsius" (if "tdpf" set to True) + "causes_overloading": does this element, when defining the N-1 case, cause overloading of other elements? the + overloading is defined by net.line["max_loading_percent_nminus1"] (if set) or net.line["max_loading_percent"] + "cause_element": element ("line", "trafo", "trafo3w") that causes max. loading of this element + "cause_index": index of the element ("line", "trafo", "trafo3w") that causes max. loading of this element + "congestion_caused_mva": overall congestion in the grid in MVA during the N-1 case due to the failure of the element + + INPUT + ---------- + **net** - pandapowerNet + **nminus1_cases** - dict + describes all N-1 cases, e.g. {"line": {"index": [1, 2, 3]}, "trafo": {"index": [0]}} + Note: trafo3w is not supported + **contingency_evaluation_function** - func + function to use for power flow calculation, default pp.runpp (but only relevant for N-0 case) + """ + if not lightsim2grid_installed: + raise UserWarning("lightsim2grid package not installed. " + "Install lightsim2grid e.g. by running 'pip install lightsim2grid' in command prompt.") + # check for continuous bus index starting with 0: + n_bus = len(net.bus) + last_bus = net.bus.index[-1] + if net.bus.index[0] != 0 or last_bus != n_bus - 1 or sum(net.bus.index) != last_bus * n_bus / 2: + raise UserWarning("bus index must be continuous and start with 0 (use pandapower.create_continuous_bus_index)") + contingency_evaluation_function(net, **kwargs) + + trafo_flag = False + if np.any(net.trafo.tap_phase_shifter): + trafo_flag = True + tap_phase_shifter, tap_pos, shift_degree = _convert_trafo_phase_shifter(net) + + # setting "slack" back-and-forth is due to the difference in interpretation of generators as "distributed slack" + if net._options.get("distributed_slack", False): + slack_backup = net.gen.slack.copy() + net.gen.loc[net.gen.slack_weight != 0, 'slack'] = True + msg = "LightSim cannot handle multiple slack bus at the moment. Only the first " \ + "slack bus of pandapower will be used." + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", msg) + lightsim_grid_model = init_ls2g(net) + net.gen['slack'] = slack_backup + solver_type = SolverType.KLU if KLU_solver_available else SolverType.SparseLU + else: + lightsim_grid_model = init_ls2g(net) + solver_type = SolverType.KLUSingleSlack if KLU_solver_available else SolverType.SparseLUSingleSlack + + if trafo_flag: + net.trafo.tap_phase_shifter = tap_phase_shifter + net.trafo.tap_pos = tap_pos + net.trafo.shift_degree = shift_degree + + n_lines = len(net.line) + n_lines_cases = len(nminus1_cases.get("line", {}).get("index", [])) + n_trafos = len(net.trafo) + n_trafos_cases = len(nminus1_cases.get("trafo", {}).get("index", [])) + + # todo: add option for DC power flow + sa = ContingencyAnalysisCPP(lightsim_grid_model) + sa.change_solver(solver_type) + + map_index = {} + for element, values in nminus1_cases.items(): + index = np.array(_get_iloc_index(net[element].index.values, values["index"]), dtype=np.int64) + map_index[element] = index.copy() # copy() because += n_lines happens later + if element == "trafo3w": + raise NotImplementedError("trafo3w not implemented for lightsim2grid contingency analysis") + elif element == "trafo": + index += n_lines + sa.add_multiple_n1(index) + + # s.add_multiple_n1(net.line.index.values.astype(int)) + v_init = net._ppc["internal"]["V"] + sa.compute(v_init, net._options["max_iteration"], net._options["tolerance_mva"]) + v_res = sa.get_voltages() + sa.compute_flows() + kamps_all = sa.get_flows() + vm_pu = np.abs(v_res) + vm_pu[np.abs(vm_pu) <= 1e-6] = np.nan # voltages are 0. in case of divergence + net.res_bus["max_vm_pu"] = np.nanmax(vm_pu, axis=0) + net.res_bus["min_vm_pu"] = np.nanmin(vm_pu, axis=0) + + max_i_ka_limit_all = np.r_[ + net.line.max_i_ka.values, net.trafo.sn_mva.values / (net.trafo.vn_hv_kv.values * np.sqrt(3))] + max_loading_limit_all = np.r_[net.line["max_loading_percent_nminus1"] + if "max_loading_percent_nminus1" in net.line.columns + else net.line["max_loading_percent"] if n_lines > 0 else [], + net.trafo["max_loading_percent_nminus1"] + if "max_loading_percent_nminus1" in net.trafo.columns + else net.trafo["max_loading_percent"] if n_trafos > 0 else []] + voltage_all = np.r_[net.bus.loc[net.line.from_bus.values, "vn_kv"].values if n_lines > 0 else [], + net.trafo.vn_hv_kv if n_trafos > 0 else []] + flows_all_mva = np.nan_to_num(kamps_all * voltage_all * np.sqrt(3)) + flows_limit_all = np.nan_to_num(max_loading_limit_all / 100 * max_i_ka_limit_all * voltage_all * np.sqrt(3)) + + big_number = 1e6 + for element in ("line", "trafo"): + if len(net[element]) == 0: + continue + if element == "line": + kamps_element = kamps_all[:, 0:n_lines] + kamps_element_cause = kamps_all[0:n_lines_cases, :] + flows_element_mva = flows_all_mva[0:n_lines_cases, :] + max_i_ka_limit = max_i_ka_limit_all[0:n_lines] + else: + kamps_element = kamps_all[:, n_lines:n_lines + n_trafos] + kamps_element_cause = kamps_all[n_lines_cases:n_lines_cases + n_trafos_cases, :] + flows_element_mva = flows_all_mva[n_lines_cases:n_lines_cases + n_trafos_cases, :] + max_i_ka_limit = max_i_ka_limit_all[n_lines:n_lines + n_trafos] + # max_i_ka = np.nanmax(kamps_line, where=~np.isnan(kamps_line), axis=0, initial=0) + cause_index = np.nanargmax(kamps_element, axis=0) + cause_element = np.where(cause_index < n_lines_cases, "line", "trafo") + max_i_ka = np.nanmax(kamps_element, axis=0) + net[f"res_{element}"]["max_loading_percent"] = max_i_ka / max_i_ka_limit * 100 + min_i_ka = np.nanmin(kamps_element, axis=0, where=kamps_element != 0, + initial=big_number) # this is quite daring tbh + min_i_ka[min_i_ka == big_number] = 0 + # min_i_ka = np.nanmin(kamps_line, axis=0) + net[f"res_{element}"]["min_loading_percent"] = min_i_ka / max_i_ka_limit * 100 + causes_overloading = np.any(kamps_element_cause > max_loading_limit_all * max_i_ka_limit_all / 100, axis=1) + net[f"res_{element}"]["causes_overloading"] = False + if element in nminus1_cases: + # order of n-1 cases is always sorted, so "vertical" sorting is different than "horizontal" + net[f"res_{element}"].loc[net[element].index.values[ + np.sort(map_index[element])], "causes_overloading"] = causes_overloading + cause_mask = cause_element == "line" + if "line" in map_index: + cause_index[cause_mask] = net.line.index.values[np.sort(map_index["line"])[cause_index[cause_mask]]] + if "trafo" in map_index: + cause_index[~cause_mask] = net.trafo.index.values[ + np.sort(map_index["trafo"])[cause_index[~cause_mask] - n_lines_cases]] + net[f"res_{element}"]["cause_index"] = cause_index + net[f"res_{element}"]["cause_element"] = cause_element + + congestion_mva = flows_element_mva - flows_limit_all + congestion_mva = np.where(congestion_mva < 0, 0, congestion_mva) + congestion_caused = congestion_mva.sum(axis=1) + if element in nminus1_cases: + # order of n-1 cases is always sorted, so "vertical" sorting is different than "horizontal" + net[f"res_{element}"].loc[net[element].index.values[ + np.sort(map_index[element])], "congestion_caused_mva"] = congestion_caused + + +def _convert_trafo_phase_shifter(net): + tap_phase_shifter = net.trafo.tap_phase_shifter.values.copy() + # vn_hv_kv = net.trafo.vn_hv_kv.values.copy() + shift_degree = net.trafo.shift_degree.values.copy() + + tap_pos = net.trafo.tap_pos.values + tap_neutral = net.trafo.tap_neutral.values + tap_diff = tap_pos - tap_neutral + tap_step_degree = net.trafo.tap_step_degree.values.copy() + + net.trafo.loc[tap_phase_shifter, 'shift_degree'] += tap_diff[tap_phase_shifter] * tap_step_degree[tap_phase_shifter] + net.trafo["tap_pos"] = 0 + net.trafo["tap_phase_shifter"] = False + + return tap_phase_shifter, tap_pos, shift_degree + + +def _update_contingency_results(net, contingency_results, result_variables, nminus1, cause_element=None, + cause_index=None): + for element, vars in result_variables.items(): + for var in vars: + val = net[f"res_{element}"][var].values + if nminus1: + if var == "loading_percent": + s = 'max_loading_percent_nminus1' \ + if 'max_loading_percent_nminus1' in net[element].columns \ + else 'max_loading_percent' + # this part with cause_mask and max_mask is not very efficient nor orderly + loading_limit = net[element].loc[contingency_results[element]["index"], s].values + cause_mask = val > loading_limit + if np.any(cause_mask): + contingency_results[cause_element]["causes_overloading"][ + contingency_results[cause_element]["index"] == cause_index] = True + max_mask = val > contingency_results[element].get("max_loading_percent", np.full_like(val, -1)) + if np.any(max_mask): + contingency_results[element]["cause_index"][max_mask] = cause_index + contingency_results[element]["cause_element"][max_mask] = cause_element + for func, min_max in ((np.fmax, "max"), (np.fmin, "min")): + key = f"{min_max}_{var}" + func(val, + contingency_results[element].setdefault(key, np.full_like(val, np.nan, dtype=np.float64)), + out=contingency_results[element][key], + where=net[element]["in_service"].values & ~np.isnan(val)) + else: + contingency_results[element][var] = val + + +def get_element_limits(net): + """ + Construct the dictionary of element limits + + INPUT + ---------- + **net** - pandapowerNet + + OUTPUT + ------- + **element_limits** - dict + """ + element_limits = {} + if "max_vm_pu" in net.bus and "min_vm_pu" in net.bus: + fill_val = np.full_like(net.bus.index, np.nan, np.float64) + bus_index = net.bus.loc[~pd.isnull(net.bus.get("max_vm_pu", fill_val)) & + ~pd.isnull(net.bus.get("min_vm_pu", fill_val))].index.values + if len(bus_index) != 0: + element_limits.update({"bus": { + "index": bus_index, + "max_limit": net.bus.loc[bus_index, "max_vm_pu"].values, + "min_limit": net.bus.loc[bus_index, "min_vm_pu"].values, + "max_limit_nminus1": + net.line.loc[bus_index, "max_vm_nminus1_pu"].values + if "max_vm_nminus1_pu" in net.bus.columns + else net.bus.loc[bus_index, "max_vm_pu"].values, + "min_limit_nminus1": + net.line.loc[bus_index, "min_vm_nminus1_pu"].values + if "min_vm_nminus1_pu" in net.bus.columns + else net.bus.loc[bus_index, "min_vm_pu"].values}}) + + for element in ("line", "trafo", "trafo3w"): + if net[element].empty or ("max_loading_percent" not in net[element] and + "max_temperature_degree_celsius" not in net[element]): + continue + + fill_val = np.full_like(net[element].index, np.nan, np.float64) + element_index = net[element].loc[ + ~pd.isnull(net[element].get("max_loading_percent", fill_val)) | + ~pd.isnull(net[element].get("max_temperature_degree_celsius", fill_val))].index.values + if len(element_index) == 0: + continue + + d = {element: {"index": element_index}} + + if "max_loading_percent" in net[element].columns: + d[element].update({ + "max_limit": net[element].loc[element_index, "max_loading_percent"].values, + "min_limit": -net[element].loc[element_index, "max_loading_percent"].values, + "max_limit_nminus1": + net[element].loc[element_index, "max_loading_nminus1_percent"].values + if "max_loading_nminus1_percent" in net[element].columns + else net[element].loc[element_index, "max_loading_percent"].values, + "min_limit_nminus1": + -net[element].loc[element_index, "max_loading_nminus1_percent"].values + if "max_loading_nminus1_percent" in net[element].columns + else - net[element].loc[element_index, "max_loading_percent"].values}) + + if element == "line" and "max_temperature_degree_celsius" in net[element].columns: + col = "max_temperature_degree_celsius" + d[element].update({"max_temperature_degree_celsius": net.line.loc[element_index, col].values, + "min_temperature_degree_celsius": -net.line.loc[element_index, col].values}) + + element_limits.update(d) + return element_limits + + +def check_elements_within_limits(element_limits, contingency_results, nminus1=False, branch_tol=1e-3, bus_tol=1e-6): + """ + Check if elements are within limits + + INPUT + ---------- + **element_limits** - dict + **contingency_results** - dict + **nminus1** - bool + **branch_tol** - float + tolerance of the limit violation check for branch limits + **bus_tol** - float + tolerance of the limit violation check for bus limits + + OUTPUT + ------- + True if all within limits (no violations), False if any limits violated + """ + for element, values in contingency_results.items(): + limit = element_limits[element] + if 'max_temperature_degree_celsius' in limit and 'temperature_degree_celsius' in values: + if np.any(values['temperature_degree_celsius'] > limit['max_temperature_degree_celsius'] + branch_tol): + return False + if nminus1 and np.any(values['max_temperature_degree_celsius'] > + limit['max_temperature_degree_celsius'] + branch_tol): + return False + if "max_limit" not in limit: + continue + if element == "bus": + if np.any(values["vm_pu"] > limit["max_limit"] + bus_tol) or \ + np.any(values["vm_pu"] < limit["min_limit"] - bus_tol): + return False + if nminus1 and (np.any(values["max_vm_pu"] > limit["max_limit_nminus1"] + bus_tol) or + np.any(values["min_vm_pu"] < limit["min_limit_nminus1"] - bus_tol)): + return False + continue + if np.any(values['loading_percent'] > limit["max_limit"] + branch_tol): + return False + if nminus1 and np.any(values['max_loading_percent'] > limit["max_limit_nminus1"] + branch_tol): + return False + return True + + +def _get_sub_index(index, sub_index): + lookup_index = np.zeros(index.max() + 1, dtype=np.int64) + lookup_index[index] = np.arange(len(index), dtype=np.int64) + return lookup_index[sub_index] + + +def _get_iloc_index(index, sub_index): + # index_dict = {k: i for i, k in enumerate(index)} + # return np.array([index_dict[i] for i in sub_index], dtype=np.int64) + # should be the same: + continuous_index = np.arange(len(index), dtype=np.int64) + return continuous_index[_get_sub_index(index, sub_index)] + + +def _log_violation(element, var, val, limit_index, mask): + if np.any(mask): + s = ' (N-1)' if 'max' in var else '' + with np.printoptions(precision=3, suppress=True): + logger.info(f"{element}: {var}{s} violation at index {limit_index[mask]} ({val[mask]})") + + +def report_contingency_results(element_limits, contingency_results, branch_tol=1e-3, bus_tol=1e-6): + """ + Print log messages for elements with violations of limits + + INPUT + ---------- + **element_limits** - dict + **contingency_results** - dict + **branch_tol** - float + tolerance for branch results + **bus_tol** - float + tolerance for bus results + """ + for element, results in contingency_results.items(): + limit = element_limits[element] + index = _get_sub_index(results["index"], limit["index"]) + for var, val in results.items(): + tol = bus_tol if element == "bus" else branch_tol + if var == "index": + continue + if "min" in var: + mask = val[index] < limit['min_limit_nminus1'] - tol + _log_violation(element, var, val[index], limit["index"], mask) + elif "max" in var: + mask = val[index] > limit['max_limit_nminus1'] + tol + _log_violation(element, var, val[index], limit["index"], mask) + elif "cause" in var: + continue + else: + mask_max = val[index] > limit['max_limit'] + tol + _log_violation(element, var, val[index], limit["index"], mask_max) + mask_min = val[index] < limit['min_limit'] - tol + _log_violation(element, var, val[index], limit["index"], mask_min) diff --git a/lightsim2grid/tests/test_DataConverter.py b/lightsim2grid/tests/test_DataConverter.py index 0e6c5de..f477a72 100644 --- a/lightsim2grid/tests/test_DataConverter.py +++ b/lightsim2grid/tests/test_DataConverter.py @@ -8,13 +8,12 @@ import unittest import numpy as np -import pdb -from numpy.lib.arraysetops import isin import pandapower.networks as pn from lightsim2grid_cpp import PandaPowerConverter -class MakeTests(unittest.TestCase): +class TestDataConverterLegacy(unittest.TestCase): + """test the converter used for legacy pandapower version (before 2.14.something)""" def setUp(self): self.converter = PandaPowerConverter() self.tol = 1e-8 @@ -27,17 +26,17 @@ def test_case6_data(self): net = pn.case6ww() self.converter.set_sn_mva(net.sn_mva) # TODO raise an error if not set ! self.converter.set_f_hz(net.f_hz) - line_r, line_x, line_h = self.converter.get_line_param( + line_r, line_x, line_h = self.converter.get_line_param_legacy( net.line["r_ohm_per_km"].values * net.line["length_km"].values, net.line["x_ohm_per_km"].values * net.line["length_km"].values, - net.line["c_nf_per_km"].values * net.line["length_km"].values, net.line["g_us_per_km"].values * net.line["length_km"].values, + net.line["c_nf_per_km"].values * net.line["length_km"].values, net.bus.loc[net.line["from_bus"]]["vn_kv"], net.bus.loc[net.line["to_bus"]]["vn_kv"] ) res_r = np.array([0.001, 0.0005, 0.001, 0.0008, 0.0005, 0.0005, 0.001, 0.0007, 0.0012, 0.0002, 0.002]) res_x = np.array([0.002, 0.002, 0.003, 0.003, 0.0025, 0.001, 0.003, 0.002, 0.0026, 0.001, 0.004]) - res_h = np.array([4.+0.j, 4.+0.j, 6.+0.j, 6.+0.j, 6.+0.j, 2.+0.j, 4.+0.j, 5.+0.j, 5.+0.j, 2.+0.j, 8.+0.j]) + res_h = 1j * np.array([4.+0.j, 4.+0.j, 6.+0.j, 6.+0.j, 6.+0.j, 2.+0.j, 4.+0.j, 5.+0.j, 5.+0.j, 2.+0.j, 8.+0.j]) # new in pandapower 2.7.0 : order changed, sn_mva changed! order_270 = [0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 2] # for el_r, el_x in zip(line_r, line_x): print(np.where((np.abs(res_r*100.-el_r) <= 1e-6) & (np.abs(res_x*100.-el_x) <= 1e-6))[0]) @@ -49,11 +48,11 @@ def test_case30_data(self): net = pn.case30() self.converter.set_sn_mva(net.sn_mva) # TODO raise an error if not set ! self.converter.set_f_hz(net.f_hz) - line_r, line_x, line_h = self.converter.get_line_param( + line_r, line_x, line_h = self.converter.get_line_param_legacy( net.line["r_ohm_per_km"].values * net.line["length_km"].values, net.line["x_ohm_per_km"].values * net.line["length_km"].values, - net.line["c_nf_per_km"].values * net.line["length_km"].values, net.line["g_us_per_km"].values * net.line["length_km"].values, + net.line["c_nf_per_km"].values * net.line["length_km"].values, net.bus.loc[net.line["from_bus"]]["vn_kv"], net.bus.loc[net.line["to_bus"]]["vn_kv"] ) @@ -69,12 +68,12 @@ def test_case30_data(self): 0.0018, 0.0027, 0.0033, 0.0038, 0.0021, 0.004 , 0.0042, 0.006 , 0.0045, 0.002 , 0.002 , 0.0006, 0.0018, 0.0004, 0.0012, 0.0008, 0.0004]) - res_h = np.array([3.+0.j, 2.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, - 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 2.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, - 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, - 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, - 0.+0.j, 2.+0.j, 2.+0.j, 1.+0.j, 2.+0.j, 0.+0.j, 1.+0.j, 1.+0.j, - 0.+0.j]) + res_h = 1j * np.array([3.+0.j, 2.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, + 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 2.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, + 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, + 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, + 0.+0.j, 2.+0.j, 2.+0.j, 1.+0.j, 2.+0.j, 0.+0.j, 1.+0.j, 1.+0.j, + 0.+0.j]) # new in pandapower 2.7.0 : order changed, sn_mva changed! # for el_r, el_x, el_h in zip(line_r, line_x, line_h): # print(np.where((np.abs(res_r*100.-el_r) <= 1e-6) & @@ -90,11 +89,11 @@ def test_case118_data(self): net = pn.case118() self.converter.set_sn_mva(net.sn_mva) # TODO raise an error if not set ! self.converter.set_f_hz(net.f_hz) - line_r, line_x, line_h = self.converter.get_line_param( + line_r, line_x, line_h = self.converter.get_line_param_legacy( net.line["r_ohm_per_km"].values * net.line["length_km"].values, net.line["x_ohm_per_km"].values * net.line["length_km"].values, - net.line["c_nf_per_km"].values * net.line["length_km"].values, net.line["g_us_per_km"].values * net.line["length_km"].values, + net.line["c_nf_per_km"].values * net.line["length_km"].values, net.bus.loc[net.line["from_bus"]]["vn_kv"], net.bus.loc[net.line["to_bus"]]["vn_kv"] ) @@ -156,41 +155,41 @@ def test_case118_data(self): 1.500e-03, 1.350e-04, 5.610e-04, 3.760e-04, 2.000e-04, 9.860e-04, 6.820e-04, 3.020e-04, 9.190e-04, 9.190e-04, 2.180e-03, 1.170e-03, 1.015e-03, 2.778e-03, 3.240e-03, 1.270e-03, 4.115e-03]) - res_h = np.array([ 2.54 +0.j, 1.082+0.j, 0.502+0.j, 0.878+0.j, 4.88 +0.j, - 4.444+0.j, 1.178+0.j, 3.368+0.j, 3.6 +0.j, 12.4 +0.j, - 1.034+0.j, 3.68 +0.j, 10.38 +0.j, 1.572+0.j, 4.978+0.j, - 1.264+0.j, 0.648+0.j, 4.72 +0.j, 2.28 +0.j, 1.87 +0.j, - 8.174+0.j, 3.796+0.j, 2.58 +0.j, 3.48 +0.j, 4.06 +0.j, - 1.234+0.j, 2.76 +0.j, 2.76 +0.j, 4.7 +0.j, 1.934+0.j, - 5.28 +0.j, 10.6 +0.j, 2.14 +0.j, 5.48 +0.j, 4.14 +0.j, - 0.874+0.j, 3.268+0.j, 2.18 +0.j, 4.06 +0.j, 1.876+0.j, - 1.11 +0.j, 4.94 +0.j, 5.44 +0.j, 2.3 +0.j, 2.54 +0.j, - 2.86 +0.j, 1.876+0.j, 5.46 +0.j, 4.72 +0.j, 6.04 +0.j, - 1.474+0.j, 2.4 +0.j, 4.76 +0.j, 2.16 +0.j, 3.28 +0.j, - 1.464+0.j, 2.94 +0.j, 1.816+0.j, 5.36 +0.j, 5.41 +0.j, - 4.07 +0.j, 4.08 +0.j, 6.2 +0.j, 0.986+0.j, 1.434+0.j, - 4.72 +0.j, 1.844+0.j, 4.72 +0.j, 6.268+0.j, 0.76 +0.j, - 4.61 +0.j, 2.02 +0.j, 2. +0.j, 6.2 +0.j, 0.768+0.j, - 5.18 +0.j, 1.628+0.j, 1.972+0.j, 0.276+0.j, 5.02 +0.j, - 3.58 +0.j, 1.198+0.j, 1.356+0.j, 2.14 +0.j, 4.44 +0.j, - 0.21 +0.j, 4.66 +0.j, 1.298+0.j, 1.142+0.j, 2.98 +0.j, - 1.01 +0.j, 2.16 +0.j, 2.46 +0.j, 4.04 +0.j, 4.98 +0.j, - 8.64 +0.j, 2.84 +0.j, 17.64 +0.j, 2.16 +0.j, 2.38 +0.j, - 51.4 +0.j, 90.8 +0.j, 3.99 +0.j, 0.83 +0.j, 11.73 +0.j, - 2.51 +0.j, 1.926+0.j, 1.426+0.j, 3.194+0.j, 6.32 +0.j, - 0.268+0.j, 1.318+0.j, 3.66 +0.j, 0.568+0.j, 0.984+0.j, - 2.7 +0.j, 4.2 +0.j, 42.2 +0.j, 0.55 +0.j, 1.552+0.j, - 1.222+0.j, 4.66 +0.j, 3.44 +0.j, 6.068+0.j, 4.226+0.j, - 2.24 +0.j, 3.32 +0.j, 3.16 +0.j, 4.72 +0.j, 116.2 +0.j, - 1.604+0.j, 8.6 +0.j, 8.6 +0.j, 4.44 +0.j, 1.258+0.j, - 1.874+0.j, 3.42 +0.j, 1.396+0.j, 4.058+0.j, 3.1 +0.j, - 123. +0.j, 7.38 +0.j, 7.3 +0.j, 2.02 +0.j, 0.732+0.j, - 0.374+0.j, 2.42 +0.j, 3.32 +0.j, 2.42 +0.j, 1.788+0.j, - 5.98 +0.j, 1.748+0.j, 5.69 +0.j, 5.36 +0.j, 5.646+0.j, - 3.76 +0.j, 3.88 +0.j, 1.456+0.j, 1.468+0.j, 0.98 +0.j, - 21.6 +0.j, 104.6 +0.j, 1.738+0.j, 38. +0.j, 2.48 +0.j, - 2.48 +0.j, 5.78 +0.j, 3.1 +0.j, 2.682+0.j, 7.092+0.j, - 8.28 +0.j, 12.2 +0.j, 10.198+0.j]) + res_h = 1j * np.array([ 2.54 +0.j, 1.082+0.j, 0.502+0.j, 0.878+0.j, 4.88 +0.j, + 4.444+0.j, 1.178+0.j, 3.368+0.j, 3.6 +0.j, 12.4 +0.j, + 1.034+0.j, 3.68 +0.j, 10.38 +0.j, 1.572+0.j, 4.978+0.j, + 1.264+0.j, 0.648+0.j, 4.72 +0.j, 2.28 +0.j, 1.87 +0.j, + 8.174+0.j, 3.796+0.j, 2.58 +0.j, 3.48 +0.j, 4.06 +0.j, + 1.234+0.j, 2.76 +0.j, 2.76 +0.j, 4.7 +0.j, 1.934+0.j, + 5.28 +0.j, 10.6 +0.j, 2.14 +0.j, 5.48 +0.j, 4.14 +0.j, + 0.874+0.j, 3.268+0.j, 2.18 +0.j, 4.06 +0.j, 1.876+0.j, + 1.11 +0.j, 4.94 +0.j, 5.44 +0.j, 2.3 +0.j, 2.54 +0.j, + 2.86 +0.j, 1.876+0.j, 5.46 +0.j, 4.72 +0.j, 6.04 +0.j, + 1.474+0.j, 2.4 +0.j, 4.76 +0.j, 2.16 +0.j, 3.28 +0.j, + 1.464+0.j, 2.94 +0.j, 1.816+0.j, 5.36 +0.j, 5.41 +0.j, + 4.07 +0.j, 4.08 +0.j, 6.2 +0.j, 0.986+0.j, 1.434+0.j, + 4.72 +0.j, 1.844+0.j, 4.72 +0.j, 6.268+0.j, 0.76 +0.j, + 4.61 +0.j, 2.02 +0.j, 2. +0.j, 6.2 +0.j, 0.768+0.j, + 5.18 +0.j, 1.628+0.j, 1.972+0.j, 0.276+0.j, 5.02 +0.j, + 3.58 +0.j, 1.198+0.j, 1.356+0.j, 2.14 +0.j, 4.44 +0.j, + 0.21 +0.j, 4.66 +0.j, 1.298+0.j, 1.142+0.j, 2.98 +0.j, + 1.01 +0.j, 2.16 +0.j, 2.46 +0.j, 4.04 +0.j, 4.98 +0.j, + 8.64 +0.j, 2.84 +0.j, 17.64 +0.j, 2.16 +0.j, 2.38 +0.j, + 51.4 +0.j, 90.8 +0.j, 3.99 +0.j, 0.83 +0.j, 11.73 +0.j, + 2.51 +0.j, 1.926+0.j, 1.426+0.j, 3.194+0.j, 6.32 +0.j, + 0.268+0.j, 1.318+0.j, 3.66 +0.j, 0.568+0.j, 0.984+0.j, + 2.7 +0.j, 4.2 +0.j, 42.2 +0.j, 0.55 +0.j, 1.552+0.j, + 1.222+0.j, 4.66 +0.j, 3.44 +0.j, 6.068+0.j, 4.226+0.j, + 2.24 +0.j, 3.32 +0.j, 3.16 +0.j, 4.72 +0.j, 116.2 +0.j, + 1.604+0.j, 8.6 +0.j, 8.6 +0.j, 4.44 +0.j, 1.258+0.j, + 1.874+0.j, 3.42 +0.j, 1.396+0.j, 4.058+0.j, 3.1 +0.j, + 123. +0.j, 7.38 +0.j, 7.3 +0.j, 2.02 +0.j, 0.732+0.j, + 0.374+0.j, 2.42 +0.j, 3.32 +0.j, 2.42 +0.j, 1.788+0.j, + 5.98 +0.j, 1.748+0.j, 5.69 +0.j, 5.36 +0.j, 5.646+0.j, + 3.76 +0.j, 3.88 +0.j, 1.456+0.j, 1.468+0.j, 0.98 +0.j, + 21.6 +0.j, 104.6 +0.j, 1.738+0.j, 38. +0.j, 2.48 +0.j, + 2.48 +0.j, 5.78 +0.j, 3.1 +0.j, 2.682+0.j, 7.092+0.j, + 8.28 +0.j, 12.2 +0.j, 10.198+0.j]) # new in pandapower 2.7.0 : order changed, sn_mva changed! # for el_r, el_x in zip(line_r, line_x): print(np.where((np.abs(res_r*100.-el_r) <= 1e-6) & (np.abs(res_x*100.-el_x) <= 1e-6))[0]) # for el in [el for el in np.arange(line_r.shape[0]) if not np.isin(el, order_270)]: print(el) @@ -243,13 +242,13 @@ def test_case118_data(self): 4.04933224e-05, 3.88000000e-04, 3.75000000e-04, 3.86000000e-04, 2.68000000e-04, 3.70000000e-04, 1.59594718e-04, 3.70000000e-04, 2.01181945e-04]) - trafo_h_res = np.array([ 0. -0.j , 0. -0.j , - 0. -0.j , 4.4602909 -0.00140652j, - 16.40272367-0.00022869j, 0. -0.j , - 0. -0.j , 0. -0.j , - 0. -0.j , 0. -0.j , - 63.96323106-0.01411497j, 0. -0.j , - 81.1310369 -0.02879733j]) + trafo_h_res = 1j * np.array([ 0. -0.j , 0. -0.j , + 0. -0.j , 4.4602909 -0.00140652j, + 16.40272367-0.00022869j, 0. -0.j , + 0. -0.j , 0. -0.j , + 0. -0.j , 0. -0.j , + 63.96323106-0.01411497j, 0. -0.j , + 81.1310369 -0.02879733j]) # new in pandapower 2.7.0 : order changed, sn_mva changed! # for el_r, el_x in zip(trafo_r, trafo_x): print(np.where((np.abs(trafo_r_res*100.-el_r) <= 1e-6) & (np.abs(trafo_x_res*100.-el_x) <= 1e-6))[0]) # for el in [el for el in np.arange(line_r.shape[0]) if not np.isin(el, order_270)]: print(el) diff --git a/lightsim2grid/tests/test_SameResPP.py b/lightsim2grid/tests/test_SameResPP.py index 4ecad36..dd995c2 100644 --- a/lightsim2grid/tests/test_SameResPP.py +++ b/lightsim2grid/tests/test_SameResPP.py @@ -385,15 +385,15 @@ def _aux_test(self, pn_net): vn_trafo_hv, vn_trafo_lv, shift_pp = _calc_tap_from_dataframe(pp_net, trafo_df) ratio = _calc_nominal_ratio_from_dataframe(ppc, trafo_df, vn_trafo_hv, vn_trafo_lv, bus_lookup) r_t, x_t, b_t = _calc_r_x_y_from_dataframe(pp_net, trafo_df, vn_trafo_lv, vn_lv, pp_net.sn_mva) - + b_t *= 1j # to fix https://github.com/BDonnot/lightsim2grid/issues/88 + # check where there are mismatch if any val_r_pp = r_t val_r_me = trafo_r all_equals_r = np.abs(val_r_pp - val_r_me) <= self.tol if not np.all(all_equals_r): - test_ok = False - print(f"\t Error: some trafo resistance are not equal, max error: {np.max(np.abs(val_r_pp - val_r_me)):.5f}") - + raise AssertionError(f"\t Error: some trafo resistance are not equal, max error: {np.max(np.abs(val_r_pp - val_r_me)):.5f}") + val_x_pp = x_t val_x_me = trafo_x all_equals_x = np.abs(val_x_pp - val_x_me) <= self.tol diff --git a/lightsim2grid/tests/test_init_from_pypowsybl.py b/lightsim2grid/tests/test_init_from_pypowsybl.py index a966719..691b377 100644 --- a/lightsim2grid/tests/test_init_from_pypowsybl.py +++ b/lightsim2grid/tests/test_init_from_pypowsybl.py @@ -29,7 +29,8 @@ def get_pypo_grid(self): return pp.network.create_ieee14() def get_slackbus_id(self): - # id in pandapower + # id in pandapower which is the same than the ID in the pypowsybl network when loaded from disk + # but might not be the same as the lightsim2grid gridmodel (if sort_index is True, which is the default) return 0 def get_equiv_pdp_grid(self): @@ -64,7 +65,7 @@ def setUp(self) -> None: self.ref_samecase = None # init lightsim2grid model - self.gridmodel = init_from_pypowsybl(self.network_ref, slack_bus_id=self.get_slackbus_id()) + self.gridmodel, self.el_ids = init_from_pypowsybl(self.network_ref, slack_bus_id=self.get_slackbus_id(), return_sub_id=True) # use some data self.nb_bus_total = self.network_ref.get_buses().shape[0] @@ -128,6 +129,8 @@ def test_compare_pp(self): # # self.pp_samecase["_ppc"]["internal"]["Bbus"] # self.pp_samecase["_ppc"]["internal"]["Ybus"][64,67] # self.pp_samecase["_ppc"]["internal"]["bus"] + # bus_id = self.el_ids[0] + # reorder = np.argsort([int(el.lstrip("VL").rstrip("0").rstrip("_")) for el in bus_id.index]).reshape(1, -1) if v_ls_ref is not None: max_ = np.abs(v_ls[reorder] - v_ls_ref).max() assert max_ <= self.tol_eq, f"error for vresults for dc: {max_:.2e}" diff --git a/lightsim2grid/tests/test_integration_pandapower.py b/lightsim2grid/tests/test_integration_pandapower.py new file mode 100644 index 0000000..d2c1a9c --- /dev/null +++ b/lightsim2grid/tests/test_integration_pandapower.py @@ -0,0 +1,149 @@ +# Copyright (c) 2024, RTE (https://www.rte-france.com) +# See AUTHORS.txt +# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0. +# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file, +# you can obtain one at http://mozilla.org/MPL/2.0/. +# SPDX-License-Identifier: MPL-2.0 +# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform + +import copy +import unittest +import numpy as np +from packaging import version + +import pandapower as pp +from pandapower_contingency import run_contingency, run_contingency_ls2g + +from lightsim2grid.gridmodel import init_from_pandapower +from lightsim2grid.gridmodel.from_pandapower._my_const import _MIN_PP_VERSION_ADV_GRID_MODEL + + +def run_for_from_bus_loading(net, **kwargs): + # copied from https://github.com/pawellytaev/pandapower/blob/5abbea2573a85e7fd3e2d21f3ea6c1aee5ffa3cf/pandapower/test/contingency/test_contingency.py#L331 + pp.runpp(net, **kwargs) + net.res_line["loading_percent"] = net.res_line.i_from_ka / net.line.max_i_ka * 100. + # 8, 37, 67, 70, 85 + # print(net.res_bus.iloc[37]["vm_pu"]) + if len(net.trafo) > 0: + max_i_ka_limit = net.trafo.sn_mva.values / (net.trafo.vn_hv_kv.values * np.sqrt(3)) + net.res_trafo["loading_percent"] = net.res_trafo.i_hv_ka / max_i_ka_limit * 100. + + +class GridModelInit(unittest.TestCase): + def setUp(self) -> None: + if version.parse(pp.__version__) < _MIN_PP_VERSION_ADV_GRID_MODEL: + self.skipTest(f"Pandapower too old for this test: requires >={_MIN_PP_VERSION_ADV_GRID_MODEL}, found {pp.__version__}") + return super().setUp() + + def _aux_test_init(self, case=None): + if case is None: + case = pp.networks.case118() + + gridmodel = init_from_pandapower(case) + + # run powerflow + pp.runpp(case) + resV = gridmodel.ac_pf(np.ones(case.bus.shape[0], dtype=np.complex128), 10, 1e-8) + + # check Ybus + Ybus_pp = case["_ppc"]["internal"]["Ybus"] + Ybus_ls = gridmodel.get_Ybus() + assert Ybus_pp.shape == Ybus_ls.shape + assert Ybus_ls.size == Ybus_pp.size + assert (Ybus_ls.nonzero()[0] == Ybus_pp.nonzero()[0]).all() + assert (Ybus_ls.nonzero()[1] == Ybus_pp.nonzero()[1]).all() + if np.abs(Ybus_ls - Ybus_pp).max() >= 1e-8: + tol = 1e-1 + while True: + row_, col_ = (np.abs(Ybus_ls - Ybus_pp).toarray() >= tol).nonzero() + if len(row_) >= 1: + raise AssertionError("Error for some Ybus coeffs. " + f"There are {(np.abs(Ybus_ls - Ybus_pp) >= 1e-8).sum()} difference (>= 1e-8). " + f"At tolerance {tol:.1e}, issues are with {row_}, {col_}") + tol /= 10 + # check Sbus + # TODO + + # check resulting voltages + # TODO + + # check flows + # TODO + + def test_case118(self): + net = pp.networks.case118() + self._aux_test_init(net) + + +class PandapowerContingencyIntegration(unittest.TestCase): + def setUp(self) -> None: + if version.parse(pp.__version__) < _MIN_PP_VERSION_ADV_GRID_MODEL: + self.skipTest(f"Pandapower too old for this test: requires >={_MIN_PP_VERSION_ADV_GRID_MODEL}, found {pp.__version__}") + return super().setUp() + + def _aux_test_case(self, case=None, nminus1_cases=None): + """test inspired from the test provided in the github issue https://github.com/BDonnot/lightsim2grid/issues/88#issuecomment-2265150641 + linked https://github.com/pawellytaev/pandapower/blob/da4b5ae6acd42e75dd37ef20d4dcbd823fef48d3/pandapower/test/contingency/test_contingency.py#L156""" + if case is None: + case = pp.networks.case118() + net = copy.deepcopy(case) + net2 = copy.deepcopy(net) + if nminus1_cases is None: + nminus1_cases = {"line": {"index": net.line.index.values}, + "trafo": {"index": net.trafo.index.values}} + + res = run_contingency(net2, nminus1_cases, contingency_evaluation_function=run_for_from_bus_loading) + + run_contingency_ls2g(net, nminus1_cases) + for s in ("min", "max"): + # import pdb + # pdb.set_trace() + # (np.abs(res["bus"][f"{s}_vm_pu"] - net.res_bus[f"{s}_vm_pu"].values) >= 1e-5).nonzero() + assert np.allclose(res["bus"][f"{s}_vm_pu"], net.res_bus[f"{s}_vm_pu"].values, atol=1e-8, rtol=0), f'for bus voltages for {s}, max error {np.abs(res["bus"][f"{s}_vm_pu"] - net.res_bus[f"{s}_vm_pu"].values).max():.2e}' + assert np.allclose(np.nan_to_num(res["line"][f"{s}_loading_percent"]), + net.res_line[f"{s}_loading_percent"].values, atol=1e-6, rtol=0), s + if len(net.trafo) > 0: + assert np.allclose(np.nan_to_num(res["trafo"][f"{s}_loading_percent"]), + net.res_trafo[f"{s}_loading_percent"].values, atol=1e-6, rtol=0), s + + def test_case118(self): + net = pp.networks.case118() + cont_lines = net.line.index.values + cont_trafo = net.trafo.index.values + + # not the same behaviour in pp and ls for contingency of these lines / trafo + lines_diff_behaviour = [6, 7, 103, 121, 163, 164, 170] + trafos_diff_behaviour = [11, 12] + ################# + # check that indeed the behaviour is different (ie that at least one bus is disconnected in pandapower) + for el in lines_diff_behaviour: + test_net = copy.deepcopy(net) + test_net.line.loc[el, "in_service"] = False + pp.runpp(test_net) + assert (~np.isfinite(test_net.res_bus["vm_pu"])).any() + for el in trafos_diff_behaviour: + test_net = copy.deepcopy(net) + test_net.trafo.loc[el, "in_service"] = False + pp.runpp(test_net) + assert (~np.isfinite(test_net.res_bus["vm_pu"])).any() + # end of the "test of the integration test" + ############# + # now the contingency analysis + cont_lines = np.delete(cont_lines, lines_diff_behaviour) + cont_trafo = np.delete(cont_trafo, trafos_diff_behaviour) + nminus1_cases = {"line": {"index": cont_lines}, + "trafo": {"index": cont_trafo} + } + self._aux_test_case(net, nminus1_cases) + + def test_case14(self): + net = pp.networks.case14() + nminus1_cases = {"line": {"index": net.line.index.values}, + "trafo": {"index": np.delete(net.trafo.index.values, -2) # not the same behaviour in pp and ls for contingency of penultimate trafo + }} + self._aux_test_case(net, nminus1_cases) + + +if __name__ == "__main__": + unittest.main() + \ No newline at end of file diff --git a/lightsim2grid/tests/test_issue_30.py b/lightsim2grid/tests/test_issue_30.py index 9c6d178..7ad3116 100644 --- a/lightsim2grid/tests/test_issue_30.py +++ b/lightsim2grid/tests/test_issue_30.py @@ -59,3 +59,8 @@ def test_diverge(self): # and it used to make this diverge sim_o, sim_r, sim_d, sim_i = obs.simulate(self.env.action_space()) assert not sim_d, "should have converged" + + +if __name__ == "__main__": + unittest.main() + \ No newline at end of file diff --git a/lightsim2grid/tests/test_issue_56.py b/lightsim2grid/tests/test_issue_56.py index 412615b..d3589b7 100644 --- a/lightsim2grid/tests/test_issue_56.py +++ b/lightsim2grid/tests/test_issue_56.py @@ -90,5 +90,4 @@ def setUp(self) -> None: if __name__ == "__main__": unittest.main() - \ No newline at end of file diff --git a/setup.py b/setup.py index a5f6182..d80be10 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ from pybind11.setup_helpers import Pybind11Extension, build_ext -__version__ = "0.9.0" +__version__ = "0.9.1" KLU_SOLVER_AVAILABLE = False # Try to link against SuiteSparse (if available) diff --git a/src/DataConverter.cpp b/src/DataConverter.cpp index 49c5f52..9ef5309 100644 --- a/src/DataConverter.cpp +++ b/src/DataConverter.cpp @@ -7,6 +7,8 @@ // This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. #include "DataConverter.h" +#include + void PandaPowerConverter::_check_init(){ if(sn_mva_ <= 0.){ @@ -20,17 +22,17 @@ void PandaPowerConverter::_check_init(){ std::tuple - PandaPowerConverter::get_trafo_param(const RealVect & tap_step_pct, - const RealVect & tap_pos, - const RealVect & tap_angles, - const std::vector & is_tap_hv_side, - const RealVect & vn_hv, // nominal voltage of hv bus - const RealVect & vn_lv, // nominal voltage of lv bus - const RealVect & trafo_vk_percent, - const RealVect & trafo_vkr_percent, - const RealVect & trafo_sn_trafo_mva, - const RealVect & trafo_pfe_kw, - const RealVect & trafo_i0_pct) + PandaPowerConverter::get_trafo_param_legacy(const RealVect & tap_step_pct, + const RealVect & tap_pos, + const RealVect & tap_angles, + const std::vector & is_tap_hv_side, + const RealVect & vn_hv, // nominal voltage of hv bus + const RealVect & vn_lv, // nominal voltage of lv bus + const RealVect & trafo_vk_percent, + const RealVect & trafo_vkr_percent, + const RealVect & trafo_sn_trafo_mva, + const RealVect & trafo_pfe_kw, + const RealVect & trafo_i0_pct) { //TODO consistency: move this class outside of here _check_init(); @@ -81,8 +83,6 @@ std::tuple res = @@ -126,12 +126,12 @@ std::tuple - PandaPowerConverter::get_line_param(const RealVect & branch_r, - const RealVect & branch_x, - const RealVect & branch_c, - const RealVect & branch_g, // TODO - const RealVect & branch_from_kv, - const RealVect & branch_to_kv) + PandaPowerConverter::get_line_param_legacy(const RealVect & branch_r, + const RealVect & branch_x, + const RealVect & branch_g, + const RealVect & branch_c, + const RealVect & branch_from_kv, + const RealVect & branch_to_kv) { //TODO does not use c at the moment! _check_init(); @@ -141,11 +141,50 @@ std::tuple(); + // b = 2 * net.f_hz * math.pi * line["c_nf_per_km"].values * 1e-9 * baseR * length_km * parallel + // g = line["g_us_per_km"].values * 1e-6 * baseR * length_km * parallel + // g + 1j . b + CplxVect powerlines_h = CplxVect::Constant(nb_line, 0.); + powerlines_h.array() += (1e-6 * branch_g.array().cast() + + my_i * 2.0 * f_hz_ * M_PI * 1e-9 * branch_c.array().cast()); powerlines_h.array() *= branch_from_pu.array().cast(); std::tuple res = std::tuple (std::move(powerlines_r), std::move(powerlines_x), std::move(powerlines_h)); return res; } + +std::tuple + PandaPowerConverter::get_line_param(const RealVect & branch_r, + const RealVect & branch_x, + const RealVect & branch_g, + const RealVect & branch_c, + const RealVect & branch_from_kv, + const RealVect & branch_to_kv) +{ + _check_init(); + const int nb_line = static_cast(branch_r.size()); + RealVect branch_from_pu = branch_from_kv.array() * branch_from_kv.array() / sn_mva_; + + RealVect powerlines_r = branch_r.array() / branch_from_pu.array(); + RealVect powerlines_x = branch_x.array() / branch_from_pu.array(); + + // TODO + // b = 2 * net.f_hz * math.pi * line["c_nf_per_km"].values * 1e-9 * baseR * length_km * parallel + // g = line["g_us_per_km"].values * 1e-6 * baseR * length_km * parallel + // g + 1j . b + CplxVect powerlines_h_or = CplxVect::Constant(nb_line, 0.); + powerlines_h_or.array() += (1e-6 * branch_g.array().cast() + + my_i * 2.0 * f_hz_ * M_PI * 1e-9 * branch_c.array().cast()); + powerlines_h_or.array() *= branch_from_pu.array().cast() * 0.5; + CplxVect powerlines_h_ex = powerlines_h_or; + // END TODO + + std::tuple res = std::tuple( + std::move(powerlines_r), std::move(powerlines_x), std::move(powerlines_h_or), std::move(powerlines_h_ex)); + return res; +} + diff --git a/src/DataConverter.h b/src/DataConverter.h index d8106da..b315fdc 100644 --- a/src/DataConverter.h +++ b/src/DataConverter.h @@ -33,7 +33,7 @@ class PandaPowerConverter : public BaseConstants // data converters /** - This converts the trafo from pandapower to r, x and h (pair unit) + This converts the trafo from pandapower to r, x and h (pair unit) (for legacy (<= 2.14.somthing) pandapower) **/ std::tuple + get_trafo_param_legacy(const RealVect & tap_step_pct, + const RealVect & tap_pos, + const RealVect & tap_angles, + const std::vector & is_tap_hv_side, + const RealVect & vn_hv, // nominal voltage of hv bus + const RealVect & vn_lv, // nominal voltage of lv bus + const RealVect & trafo_vk_percent, + const RealVect & trafo_vkr_percent, + const RealVect & trafo_sn_trafo_mva, + const RealVect & trafo_pfe_kw, + const RealVect & trafo_i0_pct); + + + /** + pair unit properly the powerlines (for legacy (<= 2.14.somthing) pandapower) + **/ + std::tuple + get_line_param_legacy(const RealVect & branch_r, + const RealVect & branch_x, + const RealVect & branch_g, + const RealVect & branch_c, + const RealVect & branch_from_kv, + const RealVect & branch_to_kv); + /** + pair unit properly the powerlines (for most recent pandapower) + **/ + std::tuple get_line_param(const RealVect & branch_r, const RealVect & branch_x, + const RealVect & branch_g, const RealVect & branch_c, - const RealVect & branch_g, //TODO g is not supported atm! const RealVect & branch_from_kv, const RealVect & branch_to_kv); diff --git a/src/GridModel.cpp b/src/GridModel.cpp index 558ee79..b4c826e 100644 --- a/src/GridModel.cpp +++ b/src/GridModel.cpp @@ -267,10 +267,11 @@ void GridModel::set_ls_to_orig(const IntVect & ls_to_orig){ throw std::runtime_error("Impossible to set the converter ls_to_orig: the provided vector has not the same size as the number of bus on the grid."); _ls_to_orig = ls_to_orig; const auto size = ls_to_orig.lpNorm(); - _orig_to_ls = IntVect::Zero(size); - _orig_to_ls.array() -= 1; - for(auto i = 0; i < size; ++i){ - _orig_to_ls[i] = _ls_to_orig[i]; + _orig_to_ls = IntVect::Constant(size + 1, -1); + int i = 0; + for(auto el : _ls_to_orig){ + if(el != -1) _orig_to_ls[el] = i; + ++i; } } @@ -287,7 +288,7 @@ void GridModel::set_orig_to_ls(const IntVect & orig_to_ls){ } if(nb_bus_ls != bus_vn_kv_.size()) throw std::runtime_error("Impossible to set the converter orig_to_ls: the number of 'non -1' component in the provided vector does not match the number of buses on the grid."); - _ls_to_orig = IntVect::Zero(nb_bus_ls); + _ls_to_orig = IntVect::Constant(nb_bus_ls, -1); Eigen::Index ls2or_ind = 0; for(auto or2ls_ind = 0; or2ls_ind < nb_bus_ls; ++or2ls_ind){ const auto my_ind = _orig_to_ls[or2ls_ind]; diff --git a/src/GridModel.h b/src/GridModel.h index 8bc311c..0c50403 100644 --- a/src/GridModel.h +++ b/src/GridModel.h @@ -972,6 +972,10 @@ class GridModel : public GenericContainer { n_sub_ = n_sub; } + int get_n_sub() const + { + return n_sub_; + } void set_max_nb_bus_per_sub(int max_nb_bus_per_sub) { if(bus_vn_kv_.size() != n_sub_ * max_nb_bus_per_sub){ diff --git a/src/batch_algorithm/ContingencyAnalysis.cpp b/src/batch_algorithm/ContingencyAnalysis.cpp index 6a2581a..4191019 100644 --- a/src/batch_algorithm/ContingencyAnalysis.cpp +++ b/src/batch_algorithm/ContingencyAnalysis.cpp @@ -115,6 +115,20 @@ bool ContingencyAnalysis::remove_from_Ybus(Eigen::SparseMatrix & Ybus return check_invertible(Ybus); } +IntVect ContingencyAnalysis::is_grid_connected_after_contingency(){ + const bool ac_solver_used = _solver.ac_solver_used(); + Eigen::SparseMatrix Ybus = ac_solver_used ? _grid_model.get_Ybus_solver() : _grid_model.get_dcYbus_solver(); + IntVect res = IntVect::Constant(_li_coeffs.size(), 0); + int cont_id = 0; + for(const auto & coeffs_modif: _li_coeffs){ + if(remove_from_Ybus(Ybus, coeffs_modif)) res(cont_id) = 1; + else res(cont_id) = 0; + readd_to_Ybus(Ybus, coeffs_modif); + ++cont_id; + } + return res; +} + void ContingencyAnalysis::readd_to_Ybus(Eigen::SparseMatrix & Ybus, const std::vector & coeffs) const { diff --git a/src/batch_algorithm/ContingencyAnalysis.h b/src/batch_algorithm/ContingencyAnalysis.h index e23b423..a67bf45 100644 --- a/src/batch_algorithm/ContingencyAnalysis.h +++ b/src/batch_algorithm/ContingencyAnalysis.h @@ -103,7 +103,8 @@ class ContingencyAnalysis: public BaseBatchSolverSynch // make the computation void compute(const CplxVect & Vinit, int max_iter, real_type tol); - + IntVect is_grid_connected_after_contingency(); + Eigen::Ref compute_flows() { compute_flows_from_Vs(); clean_flows(); diff --git a/src/element_container/LineContainer.cpp b/src/element_container/LineContainer.cpp index fdcf6bc..5dd6758 100644 --- a/src/element_container/LineContainer.cpp +++ b/src/element_container/LineContainer.cpp @@ -143,8 +143,8 @@ void LineContainer::_update_model_coeffs() // for AC // see https://matpower.org/docs/MATPOWER-manual.pdf eq. 3.2 const cplx_type ys = 1. / (powerlines_r_(i) + my_i * powerlines_x_(i)); - const cplx_type h_or = my_i * powerlines_h_or_(i); - const cplx_type h_ex = my_i * powerlines_h_ex_(i); + const cplx_type h_or = powerlines_h_or_(i); + const cplx_type h_ex = powerlines_h_ex_(i); yac_ff_(i) = (ys + h_or); yac_tt_(i) = (ys + h_ex); yac_tf_(i) = -ys; @@ -287,8 +287,8 @@ void LineContainer::fillBp_Bpp(std::vector > & Bp, yft_bp = -ys_bp_r; ytf_bp = -ys_bp_r; const real_type ys_bpp_r = std::imag(ys_bpp); - yff_bpp = ys_bpp_r + std::imag(my_i * powerlines_h_or_(line_id)); - ytt_bpp = ys_bpp_r + std::imag(my_i * powerlines_h_ex_(line_id)); + yff_bpp = ys_bpp_r + std::imag(powerlines_h_or_(line_id)); + ytt_bpp = ys_bpp_r + std::imag(powerlines_h_ex_(line_id)); yft_bpp = -ys_bpp_r; ytf_bpp = -ys_bpp_r; diff --git a/src/element_container/TrafoContainer.cpp b/src/element_container/TrafoContainer.cpp index 488c2e7..acb850d 100644 --- a/src/element_container/TrafoContainer.cpp +++ b/src/element_container/TrafoContainer.cpp @@ -133,7 +133,7 @@ void TrafoContainer::_update_model_coeffs() // for AC // see https://matpower.org/docs/MATPOWER-manual.pdf eq. 3.2 const cplx_type ys = 1. / (r_(i) + my_i * x_(i)); - const cplx_type h = my_i * h_(i) * 0.5; + const cplx_type h = h_(i) * 0.5; double tau = ratio_(i); if(!is_tap_hv_side_[i]) tau = my_one_ / tau; real_type theta_shift = shift_(i); @@ -452,8 +452,8 @@ void TrafoContainer::fillBp_Bpp(std::vector > & Bp, ytf_bp = -std::imag(ys_bp * emitheta_shift_bp); yft_bp = -std::imag(ys_bp * eitheta_shift_bp); const real_type ys_bpp_r = std::imag(ys_bpp); - yff_bpp = (ys_bpp_r + std::imag(0.5 * my_i * h_(tr_id))) / (tau_bpp * tau_bpp); - ytt_bpp = ys_bpp_r + std::imag(0.5 * my_i * h_(tr_id)); + yff_bpp = (ys_bpp_r + std::imag(0.5 * h_(tr_id))) / (tau_bpp * tau_bpp); + ytt_bpp = ys_bpp_r + std::imag(0.5 * h_(tr_id)); ytf_bpp = -ys_bpp_r / tau_bpp; yft_bpp = -ys_bpp_r / tau_bpp; diff --git a/src/help_fun_msg.cpp b/src/help_fun_msg.cpp index 2434b38..4107a58 100644 --- a/src/help_fun_msg.cpp +++ b/src/help_fun_msg.cpp @@ -16,7 +16,7 @@ const std::string DocSolver::get_J_python = R"mydelimiter( The "jacobian" matrix is only available for some powerflow (the one based on the Newton Raphson algorithm) and we provide it only for the last computed iteration. - .. info:: + .. note:: It is using the "solver" labelling, as this is accessed from the solvers. .. seealso:: @@ -31,7 +31,7 @@ const std::string DocSolver::get_J_python = R"mydelimiter( const std::string DocSolver::get_Va = R"mydelimiter( Returns the voltage angles for each buses as a numpy vector of real number. - .. info:: + .. note:: It is using the "solver" labelling, as this is accessed from the solvers. .. seealso:: @@ -46,7 +46,7 @@ const std::string DocSolver::get_Va = R"mydelimiter( const std::string DocSolver::get_Vm = R"mydelimiter( Returns the voltage magnitude for each buses as a numpy vector of real number. - .. info:: + .. note:: It is using the "solver" labelling, as this is accessed from the solvers. .. seealso:: @@ -60,7 +60,7 @@ const std::string DocSolver::get_Vm = R"mydelimiter( const std::string DocSolver::get_V = R"mydelimiter( Returns the complex voltage for each buses as a numpy vector of complex number. - .. info:: + .. note:: It is using the "solver" labelling, as this is accessed from the solvers. .. seealso:: @@ -1624,7 +1624,8 @@ const std::string DocIterator::DCLineInfo = R"mydelimiter( inverted, then flows should also be inverted (`xx` becomes `yy` and reciprocally). For the sake of simplicity, you can only control the active value at the `or` side of the dc powerline. The - active value at the `ex` side + active value at the `ex` side is then computed by the GridModel. + )mydelimiter"; const std::string DocIterator::target_p_or_mw = R"mydelimiter( @@ -1638,13 +1639,16 @@ const std::string DocIterator::target_p_or_mw = R"mydelimiter( const std::string DocIterator::target_vm_or_pu = R"mydelimiter( The target active voltage setpoint at the `or` side of the powerline (in pu NOT in kV). + )mydelimiter"; const std::string DocIterator::target_vm_ex_pu = R"mydelimiter( The target active voltage setpoint at the `ex` side of the powerline (in pu NOT in kV). + )mydelimiter"; const std::string DocIterator::dc_line_formula = R"mydelimiter( + .. note:: A DC line is modeled by two connected generators and some losses to convert the power from one to the other. @@ -1658,7 +1662,7 @@ const std::string DocIterator::dc_line_formula = R"mydelimiter( - if `or_mw` is positive, then `ex_mw = -1.0 * (or_mw - loss_mw) * (1.0 - 0.01 * loss_percent)` - if `or_mw` is negative, then `ex_mw = -1.0 * or_mw / (1.0 - 0.01 * loss_percent) + loss_mw - Where `or_mw` denotes the power injected at the origin side and `ex_mw` the power injected at the `extremity` + Where `or_mw` denotes the power injected at the `origin` side and `ex_mw` the power injected at the `extremity` side. .. note:: @@ -1733,10 +1737,12 @@ const std::string DocIterator::res_theta_ex_deg_dcline = R"mydelimiter( const std::string DocIterator::gen_or = R"mydelimiter( Direct access to the "or" generators, directly returns a :class:`lightsim2grid.elements.GenInfo` + )mydelimiter" + DocIterator::dc_line_formula; const std::string DocIterator::gen_ex = R"mydelimiter( - Direct access to the "ex" generators, directly returns a :class:`lightsim2grid.elements.GenInfo` + Direct access to the "ex" generators, directly returns a :class:`lightsim2grid.elements.GenInfo` + )mydelimiter" + DocIterator::dc_line_formula; const std::string DocGridModel::GridModel = R"mydelimiter( @@ -2065,7 +2071,7 @@ const std::string DocGridModel::get_Va = R"mydelimiter( You can :func:`lightsim2grid.gridmodel.GridModel.get_Va_solver` to get the previous (before 0.9.0) behaviour. - .. info:: + .. note:: You can use the :attr:`lightsim2grid.gridmodel.GridModel.id_ac_solver_to_me` (or :attr:`lightsim2grid.gridmodel.GridModel.id_dc_solver_to_me`) to know at which bus (on the grid) they corresponds. @@ -2087,7 +2093,7 @@ const std::string DocGridModel::get_Vm = R"mydelimiter( You can :func:`lightsim2grid.gridmodel.GridModel.get_Vm_solver` to get the previous (before 0.9.0) behaviour. - .. info:: + .. note:: You can use the :attr:`lightsim2grid.gridmodel.GridModel.id_ac_solver_to_me` (or :attr:`lightsim2grid.gridmodel.GridModel.id_dc_solver_to_me`) to know at which bus (on the grid) they corresponds. @@ -2148,7 +2154,7 @@ const std::string DocGridModel::get_Va_solver = R"mydelimiter( lightsim2grid version. The new version of :func:`lightsim2grid.gridmodel.GridModel.get_Va` now returns the id labelled with the gridmodel convention (for consistency). - .. info:: + .. note:: You can use the :attr:`lightsim2grid.gridmodel.GridModel.id_ac_solver_to_me` (or :attr:`lightsim2grid.gridmodel.GridModel.id_dc_solver_to_me`) to know at which bus (on the grid) they corresponds. @@ -2169,7 +2175,7 @@ const std::string DocGridModel::get_Vm_solver = R"mydelimiter( lightsim2grid version. The new version of :func:`lightsim2grid.gridmodel.GridModel.get_Vm` now returns the id labelled with the gridmodel convention (for consistency). - .. info:: + .. note:: You can use the :attr:`lightsim2grid.gridmodel.GridModel.id_ac_solver_to_me` (or :attr:`lightsim2grid.gridmodel.GridModel.id_dc_solver_to_me`) to know at which bus (on the grid) they corresponds. @@ -2189,7 +2195,7 @@ const std::string DocGridModel::get_V_solver = R"mydelimiter( lightsim2grid version. The new version of :func:`lightsim2grid.gridmodel.GridModel.get_V` now returns the id labelled with the gridmodel convention (for consistency). - .. info:: + .. note:: You can use the :attr:`lightsim2grid.gridmodel.GridModel.id_ac_solver_to_me` (or :attr:`lightsim2grid.gridmodel.GridModel.id_dc_solver_to_me`) to know at which bus (on the grid) they corresponds. diff --git a/src/main.cpp b/src/main.cpp index ac31441..7aa5326 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -659,8 +659,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def(py::init<>()) .def("set_f_hz", &PandaPowerConverter::set_f_hz) .def("set_sn_mva", &PandaPowerConverter::set_sn_mva) + .def("get_line_param_legacy", &PandaPowerConverter::get_line_param_legacy) .def("get_line_param", &PandaPowerConverter::get_line_param) - .def("get_trafo_param", &PandaPowerConverter::get_trafo_param); + .def("get_trafo_param", &PandaPowerConverter::get_trafo_param) + .def("get_trafo_param_legacy", &PandaPowerConverter::get_trafo_param_legacy); py::class_(m, "SolverControl", "TODO") .def(py::init<>()) @@ -945,6 +947,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) // auxiliary functions .def("set_n_sub", &GridModel::set_n_sub, DocGridModel::_internal_do_not_use.c_str()) + .def("get_n_sub", &GridModel::get_n_sub, DocGridModel::_internal_do_not_use.c_str()) .def("set_max_nb_bus_per_sub", &GridModel::set_max_nb_bus_per_sub, DocGridModel::_internal_do_not_use.c_str()) .def("set_load_pos_topo_vect", &GridModel::set_load_pos_topo_vect, DocGridModel::_internal_do_not_use.c_str()) .def("set_gen_pos_topo_vect", &GridModel::set_gen_pos_topo_vect, DocGridModel::_internal_do_not_use.c_str()) @@ -1023,6 +1026,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) // inspect the class .def("my_defaults", &ContingencyAnalysis::my_defaults_vect, DocSecurityAnalysis::my_defaults_vect.c_str()) + .def("is_grid_connected_after_contingency", &ContingencyAnalysis::is_grid_connected_after_contingency, DocGridModel::_internal_do_not_use.c_str()) // TODO // perform the computation .def("compute", &ContingencyAnalysis::compute, py::call_guard(), DocSecurityAnalysis::compute.c_str())