diff --git a/cpp/powsybl-cpp/powsybl-cpp.cpp b/cpp/powsybl-cpp/powsybl-cpp.cpp index 40f2d4a65c..f6e29b0f99 100644 --- a/cpp/powsybl-cpp/powsybl-cpp.cpp +++ b/cpp/powsybl-cpp/powsybl-cpp.cpp @@ -587,7 +587,7 @@ std::shared_ptr SensitivityAnalysisParameters:: delete ptr; }); } - + void SensitivityAnalysisParameters::load_to_c_struct(sensitivity_analysis_parameters& params) const { params.flow_flow_sensitivity_value_threshold = flow_flow_sensitivity_value_threshold; params.voltage_voltage_sensitivity_value_threshold = voltage_voltage_sensitivity_value_threshold; diff --git a/pypowsybl/__init__.py b/pypowsybl/__init__.py index 40f64ebdd5..1b953797e8 100644 --- a/pypowsybl/__init__.py +++ b/pypowsybl/__init__.py @@ -23,7 +23,8 @@ grid2op, rao, report, - voltage_initializer + voltage_initializer, + opf ) from pypowsybl.network import per_unit_view @@ -48,7 +49,8 @@ "grid2op", "dynamic", "rao", - "report" + "report", + "opf" ] diff --git a/pypowsybl/opf/__init__.py b/pypowsybl/opf/__init__.py new file mode 100644 index 0000000000..e498e9b9f0 --- /dev/null +++ b/pypowsybl/opf/__init__.py @@ -0,0 +1 @@ +from .impl.opf import OptimalPowerFlow, run_ac \ No newline at end of file diff --git a/pypowsybl/opf/impl/__init__.py b/pypowsybl/opf/impl/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pypowsybl/opf/impl/bounds/__init__.py b/pypowsybl/opf/impl/bounds/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pypowsybl/opf/impl/bounds/battery_power_bounds.py b/pypowsybl/opf/impl/bounds/battery_power_bounds.py new file mode 100644 index 0000000000..8b16ba69ef --- /dev/null +++ b/pypowsybl/opf/impl/bounds/battery_power_bounds.py @@ -0,0 +1,29 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + +class BatteryPowerBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + # battery active and reactive power bounds + for bat_num, row in enumerate(network_cache.batteries.itertuples()): + if row.bus_id: + p_bounds = Bounds(row.min_p, row.max_p).mirror() + logger.log(TRACE_LEVEL, f"Add active power bounds {p_bounds} to battery '{row.Index}' (num={bat_num})") + bat_p_index = variable_context.bat_p_num_2_index[bat_num] + model.set_variable_bounds(variable_context.bat_p_vars[bat_p_index], *Bounds.fix(row.Index, p_bounds.min_value, p_bounds.max_value)) + + bat_q_index = variable_context.bat_q_num_2_index[bat_num] + if bat_q_index != -1: # valid + q_bounds = Bounds.get_reactive_power_bounds(row).reduce(parameters.reactive_bounds_reduction).mirror() + logger.log(TRACE_LEVEL, f"Add reactive power bounds {q_bounds} to battery '{row.Index}' (num={bat_num})") + model.set_variable_bounds(variable_context.bat_q_vars[bat_q_index], *Bounds.fix(row.Index, q_bounds.min_value, q_bounds.max_value)) diff --git a/pypowsybl/opf/impl/bounds/bus_voltage_bounds.py b/pypowsybl/opf/impl/bounds/bus_voltage_bounds.py new file mode 100644 index 0000000000..1d357a8518 --- /dev/null +++ b/pypowsybl/opf/impl/bounds/bus_voltage_bounds.py @@ -0,0 +1,23 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + + +class BusVoltageBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + for bus_num, row in enumerate(network_cache.buses.itertuples()): + v_bounds = Bounds.get_voltage_bounds(row.low_voltage_limit, row.high_voltage_limit, parameters.default_voltage_bounds) + logger.log(TRACE_LEVEL, f"Add voltage magnitude bounds {v_bounds} to bus '{row.Index}' (num={bus_num})'") + model.set_variable_bounds(variable_context.v_vars[bus_num], + *Bounds.fix(row.Index, v_bounds.min_value, v_bounds.max_value)) + model.set_variable_start(variable_context.v_vars[bus_num], 1.0) diff --git a/pypowsybl/opf/impl/bounds/dangling_line_voltage_bounds.py b/pypowsybl/opf/impl/bounds/dangling_line_voltage_bounds.py new file mode 100644 index 0000000000..f70c21113d --- /dev/null +++ b/pypowsybl/opf/impl/bounds/dangling_line_voltage_bounds.py @@ -0,0 +1,25 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + + +class DanglingLineVoltageBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + for dl_num, row in enumerate(network_cache.dangling_lines.itertuples()): + if row.bus_id: + v_bounds = Bounds.get_voltage_bounds(None, None, parameters.default_voltage_bounds) + logger.log(TRACE_LEVEL, f"Add voltage magnitude bounds {v_bounds} to dangling line bus '{row.Index}' (num={dl_num})'") + dl_index = variable_context.dl_num_2_index[dl_num] + model.set_variable_bounds(variable_context.dl_v_vars[dl_index], + *Bounds.fix(row.Index, v_bounds.min_value, v_bounds.max_value)) + model.set_variable_start(variable_context.dl_v_vars[dl_index], 1.0) diff --git a/pypowsybl/opf/impl/bounds/generator_power_bounds.py b/pypowsybl/opf/impl/bounds/generator_power_bounds.py new file mode 100644 index 0000000000..7c672dc512 --- /dev/null +++ b/pypowsybl/opf/impl/bounds/generator_power_bounds.py @@ -0,0 +1,29 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + +class GeneratorPowerBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + # generator active and reactive power bounds + for gen_num, row in enumerate(network_cache.generators.itertuples()): + if row.bus_id: + p_bounds = Bounds(row.min_p, row.max_p).mirror() + logger.log(TRACE_LEVEL, f"Add active power bounds {p_bounds} to generator '{row.Index}' (num={gen_num})") + gen_p_index = variable_context.gen_p_num_2_index[gen_num] + model.set_variable_bounds(variable_context.gen_p_vars[gen_p_index], *Bounds.fix(row.Index, p_bounds.min_value, p_bounds.max_value)) + + gen_q_index = variable_context.gen_q_num_2_index[gen_num] + if gen_q_index != -1: # valid + q_bounds = Bounds.get_reactive_power_bounds(row).reduce(parameters.reactive_bounds_reduction).mirror() + logger.log(TRACE_LEVEL, f"Add reactive power bounds {q_bounds} to generator '{row.Index}' (num={gen_num})") + model.set_variable_bounds(variable_context.gen_q_vars[gen_q_index], *Bounds.fix(row.Index, q_bounds.min_value, q_bounds.max_value)) diff --git a/pypowsybl/opf/impl/bounds/slack_bus_angle_bounds.py b/pypowsybl/opf/impl/bounds/slack_bus_angle_bounds.py new file mode 100644 index 0000000000..d7317dc7c5 --- /dev/null +++ b/pypowsybl/opf/impl/bounds/slack_bus_angle_bounds.py @@ -0,0 +1,20 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + +class SlackBusAngleBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + # slack bus angle forced to 0 + slack_bus_id = network_cache.slack_terminal.iloc[0].bus_id if len(network_cache.slack_terminal) > 0 else network_cache.buses.iloc[0].name + slack_bus_num = network_cache.buses.index.get_loc(slack_bus_id) + model.set_variable_bounds(variable_context.ph_vars[slack_bus_num], 0.0, 0.0) + logger.info(f"Angle reference is at bus '{slack_bus_id}' (num={slack_bus_num})") diff --git a/pypowsybl/opf/impl/bounds/transformer_3w_middle_voltage_bounds.py b/pypowsybl/opf/impl/bounds/transformer_3w_middle_voltage_bounds.py new file mode 100644 index 0000000000..24ec8db7ee --- /dev/null +++ b/pypowsybl/opf/impl/bounds/transformer_3w_middle_voltage_bounds.py @@ -0,0 +1,25 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + + +class Transformer3wMiddleVoltageBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + for t3_num, t3_row in enumerate(network_cache.transformers_3w.itertuples()): + if t3_row.bus1_id or t3_row.bus2_id or t3_row.bus3_id: + v_bounds = Bounds.get_voltage_bounds(None, None, parameters.default_voltage_bounds) + logger.log(TRACE_LEVEL, f"Add voltage magnitude bounds {v_bounds} to 3 windings transformer middle '{t3_row.Index}' (num={t3_num})'") + t3_index = variable_context.t3_num_2_index[t3_num] + model.set_variable_bounds(variable_context.t3_middle_v_vars[t3_index], + *Bounds.fix(t3_row.Index, v_bounds.min_value, v_bounds.max_value)) + model.set_variable_start(variable_context.t3_middle_ph_vars[t3_index], 1.0) diff --git a/pypowsybl/opf/impl/bounds/vsc_cs_power_bounds.py b/pypowsybl/opf/impl/bounds/vsc_cs_power_bounds.py new file mode 100644 index 0000000000..7b81572969 --- /dev/null +++ b/pypowsybl/opf/impl/bounds/vsc_cs_power_bounds.py @@ -0,0 +1,34 @@ +import logging + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + +class VscCsPowerBounds(VariableBounds): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + # VSC converter station active and reactive power bounds + for vsc_cs_num, row in enumerate(network_cache.vsc_converter_stations.itertuples()): + if row.bus_id: + p_bounds = Bounds(-row.max_p, row.max_p).mirror() + logger.log(TRACE_LEVEL, + f"Add active power bounds {p_bounds} to VSC converter station '{row.Index}' (num={vsc_cs_num})") + model.set_variable_bounds(variable_context.vsc_cs_p_vars[vsc_cs_num], + *Bounds.fix(row.Index, p_bounds.min_value, p_bounds.max_value)) + + q_bounds = Bounds.get_reactive_power_bounds(row).reduce( + parameters.reactive_bounds_reduction).mirror() + logger.log(TRACE_LEVEL, + f"Add reactive power bounds {q_bounds} to VSC converter station '{row.Index}' (num={vsc_cs_num})") + if abs(q_bounds.max_value - q_bounds.min_value) < 1.0 / network_cache.network.nominal_apparent_power: + logger.error( + f"Too small reactive power bounds {q_bounds} for VSC converter station '{row.Index}' (num={vsc_cs_num})") + model.set_variable_bounds(variable_context.vsc_cs_q_vars[vsc_cs_num], + *Bounds.fix(row.Index, q_bounds.min_value, q_bounds.max_value)) diff --git a/pypowsybl/opf/impl/constraints/__init__.py b/pypowsybl/opf/impl/constraints/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pypowsybl/opf/impl/constraints/branch_flow_constraints.py b/pypowsybl/opf/impl/constraints/branch_flow_constraints.py new file mode 100644 index 0000000000..75dc51ae99 --- /dev/null +++ b/pypowsybl/opf/impl/constraints/branch_flow_constraints.py @@ -0,0 +1,137 @@ +from math import hypot, atan2 + +from pyoptinterface import ipopt, nl + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + +R2 = 1.0 +A2 = 0.0 + +class BranchFlowConstraints(Constraints): + @staticmethod + def add_closed_branch_constraint(a1: float, b1: float, b2: float, g1: float, g2: float, ksi: float, model, p1_var, p2_var, ph1_var, + ph2_var, q1_var, q2_var, r1: float, v1_var, v2_var, y: float): + with nl.graph(): + sin_ksi = nl.sin(ksi) + cos_ksi = nl.cos(ksi) + theta1 = ksi - a1 + A2 - ph1_var + ph2_var + theta2 = ksi + a1 - A2 + ph1_var - ph2_var + sin_theta1 = nl.sin(theta1) + cos_theta1 = nl.cos(theta1) + sin_theta2 = nl.sin(theta2) + cos_theta2 = nl.cos(theta2) + + p1_eq = r1 * v1_var * (g1 * r1 * v1_var + y * r1 * v1_var * sin_ksi - y * R2 * v2_var * sin_theta1) - p1_var + q1_eq = r1 * v1_var * (-b1 * r1 * v1_var + y * r1 * v1_var * cos_ksi - y * R2 * v2_var * cos_theta1) - q1_var + p2_eq = R2 * v2_var * (g2 * R2 * v2_var - y * r1 * v1_var * sin_theta2 + y * R2 * v2_var * sin_ksi) - p2_var + q2_eq = R2 * v2_var * (-b2 * R2 * v2_var - y * r1 * v1_var * cos_theta2 + y * R2 * v2_var * cos_ksi) - q2_var + + model.add_nl_constraint(p1_eq == 0.0) + model.add_nl_constraint(q1_eq == 0.0) + model.add_nl_constraint(p2_eq == 0.0) + model.add_nl_constraint(q2_eq == 0.0) + + @staticmethod + def add_open_side1_branch_constraint(b1: float, b2: float, g1: float, g2: float, ksi: float, model, p2_var, q2_var, v2_var, + y: float): + with nl.graph(): + sin_ksi = nl.sin(ksi) + cos_ksi = nl.cos(ksi) + + shunt = (g1 + y * sin_ksi) * (g1 + y * sin_ksi) + (-b1 + y * cos_ksi) * (-b1 + y * cos_ksi) + p2_eq = R2 * R2 * v2_var * v2_var * ( + g2 + y * y * g1 / shunt + (b1 * b1 + g1 * g1) * y * sin_ksi / shunt) - p2_var + q2_eq = -R2 * R2 * v2_var * v2_var * ( + b2 + y * y * b1 / shunt - (b1 * b1 + g1 * g1) * y * cos_ksi / shunt) - q2_var + + model.add_nl_constraint(p2_eq == 0.0) + model.add_nl_constraint(q2_eq == 0.0) + + @staticmethod + def add_open_side2_branch_constraint(b1: float, b2: float, g1: float, g2: float, ksi: float, model, p1_var, q1_var, + r1: float, v1_var, y: float): + with nl.graph(): + sin_ksi = nl.sin(ksi) + cos_ksi = nl.cos(ksi) + + shunt = (g2 + y * sin_ksi) * (g2 + y * sin_ksi) + (-b2 + y * cos_ksi) * (-b2 + y * cos_ksi) + p1_eq = r1 * r1 * v1_var * v1_var * ( + g1 + y * y * g2 / shunt + (b2 * b2 + g2 * g2) * y * sin_ksi / shunt) - p1_var + q1_eq = -r1 * r1 * v1_var * v1_var * ( + b1 + y * y * b2 / shunt - (b2 * b2 + g2 * g2) * y * cos_ksi / shunt) - q1_var + + model.add_nl_constraint(p1_eq == 0.0) + model.add_nl_constraint(q1_eq == 0.0) + + @staticmethod + def add_branch_constraint(branch_index: int, bus1_id: str, bus2_id: str, network_cache: NetworkCache, model, + r: float, x: float, g1: float, b1: float, g2: float, b2: float, r1: float, a1: float, + variable_context: VariableContext) -> None: + z = hypot(r, x) + y = 1.0 / z + ksi = atan2(r, x) + + if bus1_id and bus2_id: + bus1_num = network_cache.buses.index.get_loc(bus1_id) + bus2_num = network_cache.buses.index.get_loc(bus2_id) + v1_var = variable_context.v_vars[bus1_num] + v2_var = variable_context.v_vars[bus2_num] + ph1_var = variable_context.ph_vars[bus1_num] + ph2_var = variable_context.ph_vars[bus2_num] + p1_var = variable_context.closed_branch_p1_vars[branch_index] + q1_var = variable_context.closed_branch_q1_vars[branch_index] + p2_var = variable_context.closed_branch_p2_vars[branch_index] + q2_var = variable_context.closed_branch_q2_vars[branch_index] + + BranchFlowConstraints.add_closed_branch_constraint(a1, b1, b2, g1, g2, ksi, model, p1_var, p2_var, + ph1_var, ph2_var, q1_var, q2_var, r1, v1_var, v2_var, y) + elif bus2_id: + bus2_num = network_cache.buses.index.get_loc(bus2_id) + v2_var = variable_context.v_vars[bus2_num] + p2_var = variable_context.open_side1_branch_p2_vars[branch_index] + q2_var = variable_context.open_side1_branch_q2_vars[branch_index] + + BranchFlowConstraints.add_open_side1_branch_constraint(b1, b2, g1, g2, ksi, model, p2_var, q2_var, v2_var, y) + elif bus1_id: + bus1_num = network_cache.buses.index.get_loc(bus1_id) + v1_var = variable_context.v_vars[bus1_num] + p1_var = variable_context.open_side2_branch_p1_vars[branch_index] + q1_var = variable_context.open_side2_branch_q1_vars[branch_index] + + BranchFlowConstraints.add_open_side2_branch_constraint(b1, b2, g1, g2, ksi, model, p1_var, q1_var, r1, + v1_var, y) + + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for branch_num, row in enumerate(network_cache.lines.itertuples(index=False)): + r, x, g1, b1, g2, b2 = row.r, row.x, row.g1, row.b1, row.g2, row.b2 + r1 = 1.0 + a1 = 0.0 + branch_index = variable_context.branch_num_2_index[branch_num] + self.add_branch_constraint(branch_index, row.bus1_id, row.bus2_id, network_cache, model, + r, x, g1, b1, g2, b2, r1, a1, + variable_context) + + for transfo_num, row in enumerate(network_cache.transformers_2w.itertuples(index=True)): + r, x, g, b, rho, alpha = row.r_at_current_tap, row.x_at_current_tap, row.g_at_current_tap, row.b_at_current_tap, row.rho, row.alpha + if parameters.twt_split_shunt_admittance: + g1 = g / 2 + g2 = g / 2 + b1 = b / 2 + b2 = b / 2 + else: + g1 = g + g2 = 0 + b1 = b + b2 = 0 + r1 = rho + a1 = alpha + branch_num = len(network_cache.lines) + transfo_num + branch_index = variable_context.branch_num_2_index[branch_num] + self.add_branch_constraint(branch_index, row.bus1_id, row.bus2_id, network_cache, model, + r, x, g1, b1, g2, b2, r1, a1, + variable_context) + diff --git a/pypowsybl/opf/impl/constraints/current_limit_constraints.py b/pypowsybl/opf/impl/constraints/current_limit_constraints.py new file mode 100644 index 0000000000..25bbd2b485 --- /dev/null +++ b/pypowsybl/opf/impl/constraints/current_limit_constraints.py @@ -0,0 +1,40 @@ +import logging +from typing import Any + +import pyoptinterface as poi +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + +class CurrentLimitConstraints(Constraints): + + @staticmethod + def _create_limit_constraints_for_side(side: str, branch_row, current_limits, model, p: Any, q: Any): + if branch_row.Index in current_limits.index: + limit_row = current_limits.loc[branch_row.Index] + if branch_row.bus1_id and branch_row.bus2_id: + i = poi.ExprBuilder() + i += p * p + q * q + model.add_quadratic_constraint(i, poi.Leq, limit_row.value * limit_row.value) + logger.log(TRACE_LEVEL, + f"Add side {side} current limit constraint (value={limit_row.value}) to branch '{branch_row.Index}'") + + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + current_limits1 = network_cache.current_limits1 + current_limits2 = network_cache.current_limits2 + for branch_num, branch_row in enumerate(network_cache.lines.itertuples()): + branch_index = variable_context.branch_num_2_index[branch_num] + + self._create_limit_constraints_for_side('1', branch_row, current_limits1, model, + variable_context.closed_branch_p1_vars[branch_index], + variable_context.closed_branch_q1_vars[branch_index]) + self._create_limit_constraints_for_side('2', branch_row, current_limits2, model, + variable_context.closed_branch_p2_vars[branch_index], + variable_context.closed_branch_q2_vars[branch_index]) diff --git a/pypowsybl/opf/impl/constraints/dangling_line_flow_constraints.py b/pypowsybl/opf/impl/constraints/dangling_line_flow_constraints.py new file mode 100644 index 0000000000..aeb425eccb --- /dev/null +++ b/pypowsybl/opf/impl/constraints/dangling_line_flow_constraints.py @@ -0,0 +1,38 @@ +from math import hypot, atan2 + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.constraints.branch_flow_constraints import R2, A2, BranchFlowConstraints +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class DanglingLineFlowConstraints(Constraints): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for dl_num, row in enumerate(network_cache.dangling_lines.itertuples(index=False)): + r, x, g, b, bus_id = row.r, row.x, row.g, row.b, row.bus_id + if bus_id: + g1 = g + b1 = b + g2 = 0 + b2 = 0 + r1 = 1.0 + a1 = 0.0 + dl_index = variable_context.dl_num_2_index[dl_num] + z = hypot(r, x) + y = 1.0 / z + ksi = atan2(r, x) + bus_num = network_cache.buses.index.get_loc(bus_id) + v1_var = variable_context.v_vars[bus_num] + ph1_var = variable_context.ph_vars[bus_num] + v2_var = variable_context.dl_v_vars[dl_index] + ph2_var = variable_context.dl_ph_vars[dl_index] + p1_var = variable_context.dl_branch_p1_vars[dl_index] + q1_var = variable_context.dl_branch_q1_vars[dl_index] + p2_var = variable_context.dl_branch_p2_vars[dl_index] + q2_var = variable_context.dl_branch_q2_vars[dl_index] + + BranchFlowConstraints.add_closed_branch_constraint(a1, b1, b2, g1, g2, ksi, model, p1_var, p2_var, ph1_var, ph2_var, q1_var, q2_var, r1, v1_var, v2_var, y) diff --git a/pypowsybl/opf/impl/constraints/hvdc_line_constraints.py b/pypowsybl/opf/impl/constraints/hvdc_line_constraints.py new file mode 100644 index 0000000000..80d09013fa --- /dev/null +++ b/pypowsybl/opf/impl/constraints/hvdc_line_constraints.py @@ -0,0 +1,43 @@ +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + + +def hvdc_line_losses(p, r, sb): + return r * p * p / sb + + +def add_converter_losses(p, loss_factor): + return p * (1.0 - loss_factor / 100.0) + + +class HvdcLineConstraints(Constraints): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for hvdc_line_row in network_cache.hvdc_lines.itertuples(index=False): + cs1_id, cs2_id, r, nominal_v = hvdc_line_row.converter_station1_id, hvdc_line_row.converter_station2_id, hvdc_line_row.r, hvdc_line_row.nominal_v + cs1_num = network_cache.vsc_converter_stations.index.get_loc(cs1_id) + cs2_num = network_cache.vsc_converter_stations.index.get_loc(cs2_id) + cs1_row = network_cache.vsc_converter_stations.loc[cs1_id] + cs2_row = network_cache.vsc_converter_stations.loc[cs2_id] + cs1_index = variable_context.vsc_cs_num_2_index[cs1_num] + cs2_index = variable_context.vsc_cs_num_2_index[cs2_num] + p1_var = variable_context.vsc_cs_p_vars[cs1_index] + p2_var = variable_context.vsc_cs_p_vars[cs2_index] + + loss_factor1 = cs1_row.loss_factor + loss_factor2 = cs2_row.loss_factor + sb = network_cache.network.nominal_apparent_power + if NetworkCache.is_rectifier(cs1_id, hvdc_line_row): + p_rectifier = add_converter_losses(p1_var, loss_factor1) + p_eq = add_converter_losses(p_rectifier - hvdc_line_losses(p_rectifier, r, sb), + loss_factor2) + p2_var + else: + p_rectifier = add_converter_losses(p2_var, loss_factor2) + p_eq = add_converter_losses(p_rectifier - hvdc_line_losses(p_rectifier, r, sb), + loss_factor1) + p1_var + + model.add_quadratic_constraint(p_eq == 0.0) diff --git a/pypowsybl/opf/impl/constraints/power_balance_constraints.py b/pypowsybl/opf/impl/constraints/power_balance_constraints.py new file mode 100644 index 0000000000..9e9b573449 --- /dev/null +++ b/pypowsybl/opf/impl/constraints/power_balance_constraints.py @@ -0,0 +1,171 @@ +import pyoptinterface as poi +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + + +class PowerBalanceConstraints(Constraints): + + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for bus_expr in self.create_bus_expr_list(network_cache, variable_context): + model.add_linear_constraint(bus_expr, poi.Eq, 0.0) + + + class BusesBalance: + def __init__(self, bus_count: int): + self.p_gen = [[] for _ in range(bus_count)] + self.q_gen = [[] for _ in range(bus_count)] + self.p_load = [0.0 for _ in range(bus_count)] + self.q_load = [0.0 for _ in range(bus_count)] + + @staticmethod + def _to_expr(bus_gen: list, bus_load: list) -> list[poi.ExprBuilder]: + return [poi.ExprBuilder() + poi.quicksum(gen) - load + for gen, load in zip(bus_gen, bus_load)] + + def to_expr(self) -> list[poi.ExprBuilder]: + return self._to_expr(self.p_gen, self.p_load) + self._to_expr(self.q_gen, self.q_load) + + + @classmethod + def create_bus_expr_list(cls, network_cache, variable_context): + buses_balance = cls.BusesBalance(len(network_cache.buses)) + + # branches + for branch_num, row in enumerate(network_cache.branches.itertuples(index=False)): + branch_index = variable_context.branch_num_2_index[branch_num] + if row.bus1_id and row.bus2_id: + bus1_num = network_cache.buses.index.get_loc(row.bus1_id) + bus2_num = network_cache.buses.index.get_loc(row.bus2_id) + buses_balance.p_gen[bus1_num].append(variable_context.closed_branch_p1_vars[branch_index]) + buses_balance.q_gen[bus1_num].append(variable_context.closed_branch_q1_vars[branch_index]) + buses_balance.p_gen[bus2_num].append(variable_context.closed_branch_p2_vars[branch_index]) + buses_balance.q_gen[bus2_num].append(variable_context.closed_branch_q2_vars[branch_index]) + elif row.bus2_id: + bus2_num = network_cache.buses.index.get_loc(row.bus2_id) + buses_balance.p_gen[bus2_num].append(variable_context.open_side1_branch_p2_vars[branch_index]) + buses_balance.q_gen[bus2_num].append(variable_context.open_side1_branch_q2_vars[branch_index]) + elif row.bus1_id: + bus1_num = network_cache.buses.index.get_loc(row.bus1_id) + buses_balance.p_gen[bus1_num].append(variable_context.open_side2_branch_p1_vars[branch_index]) + buses_balance.q_gen[bus1_num].append(variable_context.open_side2_branch_q1_vars[branch_index]) + + # generators + for gen_num, gen_row in enumerate(network_cache.generators.itertuples(index=False)): + bus_id = gen_row.bus_id + if bus_id: + bus_num = network_cache.buses.index.get_loc(bus_id) + gen_p_index = variable_context.gen_p_num_2_index[gen_num] + gen_q_index = variable_context.gen_q_num_2_index[gen_num] + buses_balance.p_gen[bus_num].append(variable_context.gen_p_vars[gen_p_index]) + if gen_q_index == -1: # invalid + buses_balance.q_load[bus_num] += gen_row.target_q + else: + buses_balance.q_gen[bus_num].append(variable_context.gen_q_vars[gen_q_index]) + + # batteries + for bat_num, bat_row in enumerate(network_cache.batteries.itertuples(index=False)): + bus_id = bat_row.bus_id + if bus_id: + bus_num = network_cache.buses.index.get_loc(bus_id) + bat_p_index = variable_context.bat_p_num_2_index[bat_num] + bat_q_index = variable_context.bat_q_num_2_index[bat_num] + buses_balance.p_gen[bus_num].append(variable_context.bat_p_vars[bat_p_index]) + if bat_q_index == -1: # invalid + buses_balance.q_load[bus_num] += bat_row.target_q + else: + buses_balance.q_gen[bus_num].append(variable_context.bat_q_vars[bat_q_index]) + + # static var compensators + for svc_num, row in enumerate(network_cache.static_var_compensators.itertuples(index=False)): + bus_id = row.bus_id + if bus_id: + svc_index = variable_context.svc_num_2_index[svc_num] + bus_num = network_cache.buses.index.get_loc(bus_id) + buses_balance.q_gen[bus_num].append(variable_context.svc_q_vars[svc_index]) + + # aggregated loads + loads_sum = network_cache.loads.groupby("bus_id", as_index=False).agg({"p0": "sum", "q0": "sum"}) + for row in loads_sum.itertuples(index=False): + bus_id = row.bus_id + if bus_id: + bus_num = network_cache.buses.index.get_loc(bus_id) + buses_balance.p_load[bus_num] -= row.p0 + buses_balance.q_load[bus_num] -= row.q0 + + # shunts + for shunt_num, row in enumerate(network_cache.shunts.itertuples(index=False)): + bus_id = row.bus_id + if bus_id: + shunt_index = variable_context.shunt_num_2_index[shunt_num] + bus_num = network_cache.buses.index.get_loc(bus_id) + buses_balance.p_gen[bus_num].append(variable_context.shunt_p_vars[shunt_index]) + buses_balance.q_gen[bus_num].append(variable_context.shunt_q_vars[shunt_index]) + + # VSC converter stations + for vsc_cs_num, row in enumerate(network_cache.vsc_converter_stations.itertuples(index=False)): + bus_id = row.bus_id + if bus_id: + vsc_cs_index = variable_context.vsc_cs_num_2_index[vsc_cs_num] + bus_num = network_cache.buses.index.get_loc(bus_id) + buses_balance.p_gen[bus_num].append(variable_context.vsc_cs_p_vars[vsc_cs_index]) + buses_balance.q_gen[bus_num].append(variable_context.vsc_cs_q_vars[vsc_cs_index]) + + # dangling lines + dl_buses_balance = cls.BusesBalance(len(variable_context.dl_v_vars)) + for dl_num, row in enumerate(network_cache.dangling_lines.itertuples(index=False)): + bus_id = row.bus_id + if bus_id: + dl_index = variable_context.dl_num_2_index[dl_num] + bus_num = network_cache.buses.index.get_loc(bus_id) + buses_balance.p_gen[bus_num].append(variable_context.dl_branch_p1_vars[dl_index]) + buses_balance.q_gen[bus_num].append(variable_context.dl_branch_q1_vars[dl_index]) + dl_buses_balance.p_gen[dl_index].append(variable_context.dl_branch_p2_vars[dl_index]) + dl_buses_balance.q_gen[dl_index].append(variable_context.dl_branch_q2_vars[dl_index]) + dl_buses_balance.p_load[dl_index] -= row.p0 + dl_buses_balance.q_load[dl_index] -= row.q0 + + # 3 windings transformers + t3_buses_balance = cls.BusesBalance(len(variable_context.t3_middle_v_vars)) + for t3_num, (t3_id, t3_row) in enumerate(network_cache.transformers_3w.iterrows()): + t3_index = variable_context.t3_num_2_index[t3_num] + if t3_row.bus1_id or t3_row.bus2_id or t3_row.bus3_id: + leg1_index = variable_context.t3_leg1_num_2_index[t3_num] + leg2_index = variable_context.t3_leg2_num_2_index[t3_num] + leg3_index = variable_context.t3_leg3_num_2_index[t3_num] + + if t3_row.bus1_id: + bus1_num = network_cache.buses.index.get_loc(t3_row.bus1_id) + buses_balance.p_gen[bus1_num].append(variable_context.t3_closed_branch_p1_vars[leg1_index]) + buses_balance.q_gen[bus1_num].append(variable_context.t3_closed_branch_q1_vars[leg1_index]) + t3_buses_balance.p_gen[t3_index].append(variable_context.t3_closed_branch_p2_vars[leg1_index]) + t3_buses_balance.q_gen[t3_index].append(variable_context.t3_closed_branch_q2_vars[leg1_index]) + else: + t3_buses_balance.p_gen[t3_index].append(variable_context.t3_open_side1_branch_p2_vars[leg1_index]) + t3_buses_balance.q_gen[t3_index].append(variable_context.t3_open_side1_branch_q2_vars[leg1_index]) + + if t3_row.bus2_id: + bus2_num = network_cache.buses.index.get_loc(t3_row.bus2_id) + buses_balance.p_gen[bus2_num].append(variable_context.t3_closed_branch_p1_vars[leg2_index]) + buses_balance.q_gen[bus2_num].append(variable_context.t3_closed_branch_q1_vars[leg2_index]) + t3_buses_balance.p_gen[t3_index].append(variable_context.t3_closed_branch_p2_vars[leg2_index]) + t3_buses_balance.q_gen[t3_index].append(variable_context.t3_closed_branch_q2_vars[leg2_index]) + else: + t3_buses_balance.p_gen[t3_index].append(variable_context.t3_open_side1_branch_p2_vars[leg2_index]) + t3_buses_balance.q_gen[t3_index].append(variable_context.t3_open_side1_branch_q2_vars[leg2_index]) + + if t3_row.bus3_id: + bus3_num = network_cache.buses.index.get_loc(t3_row.bus3_id) + buses_balance.p_gen[bus3_num].append(variable_context.t3_closed_branch_p1_vars[leg3_index]) + buses_balance.q_gen[bus3_num].append(variable_context.t3_closed_branch_q1_vars[leg3_index]) + t3_buses_balance.p_gen[t3_index].append(variable_context.t3_closed_branch_p2_vars[leg3_index]) + t3_buses_balance.q_gen[t3_index].append(variable_context.t3_closed_branch_q2_vars[leg3_index]) + else: + t3_buses_balance.p_gen[t3_index].append(variable_context.t3_open_side1_branch_p2_vars[leg3_index]) + t3_buses_balance.q_gen[t3_index].append(variable_context.t3_open_side1_branch_q2_vars[leg3_index]) + + return buses_balance.to_expr() + dl_buses_balance.to_expr() + t3_buses_balance.to_expr() diff --git a/pypowsybl/opf/impl/constraints/shunt_flow_constraints.py b/pypowsybl/opf/impl/constraints/shunt_flow_constraints.py new file mode 100644 index 0000000000..d16b1936ba --- /dev/null +++ b/pypowsybl/opf/impl/constraints/shunt_flow_constraints.py @@ -0,0 +1,25 @@ +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + + +class ShuntFlowConstraints(Constraints): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for shunt_num, row in enumerate(network_cache.shunts.itertuples(index=False)): + g, b, bus_id = row.g, row.b, row.bus_id + if bus_id: + shunt_index = variable_context.shunt_num_2_index[shunt_num] + p_var = variable_context.shunt_p_vars[shunt_index] + q_var = variable_context.shunt_q_vars[shunt_index] + bus_num = network_cache.buses.index.get_loc(bus_id) + v_var = variable_context.v_vars[bus_num] + + p_eq = g * v_var * v_var - p_var + q_eq = -b * v_var * v_var - q_var + + model.add_quadratic_constraint(p_eq == 0.0) + model.add_quadratic_constraint(q_eq == 0.0) diff --git a/pypowsybl/opf/impl/constraints/static_var_compensator_reactive_limits_constraints.py b/pypowsybl/opf/impl/constraints/static_var_compensator_reactive_limits_constraints.py new file mode 100644 index 0000000000..f739d40a5f --- /dev/null +++ b/pypowsybl/opf/impl/constraints/static_var_compensator_reactive_limits_constraints.py @@ -0,0 +1,27 @@ +import pyoptinterface as poi +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + + +class StaticVarCompensatorReactiveLimitsConstraints(Constraints): + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for svc_num, row in enumerate(network_cache.static_var_compensators.itertuples(index=False)): + b_min, b_max, bus_id = row.b_min, row.b_max, row.bus_id + if bus_id: + svc_index = variable_context.svc_num_2_index[svc_num] + q_var = variable_context.svc_q_vars[svc_index] + bus_num = network_cache.buses.index.get_loc(bus_id) + v_var = variable_context.v_vars[bus_num] + q_min_expr = poi.ExprBuilder() + q_min_expr += b_min * v_var * v_var + q_min_expr -= q_var + model.add_quadratic_constraint(q_min_expr, poi.Leq, 0.0) + q_max_expr = poi.ExprBuilder() + q_max_expr += b_max * v_var * v_var + q_max_expr -= q_var + model.add_quadratic_constraint(q_max_expr, poi.Geq, 0.0) diff --git a/pypowsybl/opf/impl/constraints/transformer_3w_flow_constraints.py b/pypowsybl/opf/impl/constraints/transformer_3w_flow_constraints.py new file mode 100644 index 0000000000..240b329f9f --- /dev/null +++ b/pypowsybl/opf/impl/constraints/transformer_3w_flow_constraints.py @@ -0,0 +1,70 @@ +from math import hypot, atan2 + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.constraints.branch_flow_constraints import BranchFlowConstraints +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class Transformer3wFlowConstraints(Constraints): + @staticmethod + def _create_leg_constraint(leg_r: float, leg_x: float, leg_g: float, leg_b: float, rho: float, alpha: float, + bus_id: str, t3_index: int, leg_index: int, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model): + r = leg_r + x = leg_x + if parameters.twt_split_shunt_admittance: + g1 = leg_g / 2 + g2 = leg_g / 2 + b1 = leg_b / 2 + b2 = leg_b / 2 + else: + g1 = leg_g + g2 = 0.0 + b1 = leg_b + b2 = 0.0 + r1 = rho + a1 = alpha + z = hypot(r, x) + y = 1.0 / z + ksi = atan2(r, x) + if bus_id: + bus_num = network_cache.buses.index.get_loc(bus_id) + v1_var = variable_context.v_vars[bus_num] + ph1_var = variable_context.ph_vars[bus_num] + v2_var = variable_context.t3_middle_v_vars[t3_index] + ph2_var = variable_context.t3_middle_ph_vars[t3_index] + p1_var = variable_context.t3_closed_branch_p1_vars[leg_index] + q1_var = variable_context.t3_closed_branch_q1_vars[leg_index] + p2_var = variable_context.t3_closed_branch_p2_vars[leg_index] + q2_var = variable_context.t3_closed_branch_q2_vars[leg_index] + + BranchFlowConstraints.add_closed_branch_constraint(a1, b1, b2, g1, g2, ksi, model, p1_var, p2_var, ph1_var, ph2_var, q1_var, q2_var, r1, v1_var, v2_var, y) + else: + v2_var = variable_context.t3_middle_v_vars[t3_index] + ph2_var = variable_context.t3_middle_ph_vars[t3_index] + p2_var = variable_context.t3_open_side1_branch_p2_vars[leg_index] + q2_var = variable_context.t3_open_side1_branch_q2_vars[leg_index] + BranchFlowConstraints.add_open_side1_branch_constraint(b1, b2, g1, g2, ksi, model, p2_var, q2_var, v2_var, y) + + def add(self, parameters: ModelParameters, network_cache: NetworkCache, + variable_context: VariableContext, model: ipopt.Model) -> None: + for t3_num, row in enumerate(network_cache.transformers_3w.itertuples(index=False)): + t3_index = variable_context.t3_num_2_index[t3_num] + leg1_r, leg2_r, leg3_r = row.r1_at_current_tap, row.r2_at_current_tap, row.r3_at_current_tap + leg1_x, leg2_x, leg3_x = row.x1_at_current_tap, row.x2_at_current_tap, row.x3_at_current_tap + leg1_g, leg2_g, leg3_g = row.g1_at_current_tap, row.g2_at_current_tap, row.g3_at_current_tap + leg1_b, leg2_b, leg3_b = row.b1_at_current_tap, row.b2_at_current_tap, row.b3_at_current_tap + rho1, rho2, rho3 = row.rho1, row.rho2, row.rho3 + alpha1, alpha2, alpha3 = row.alpha1, row.alpha2, row.alpha3 + bus1_id, bus2_id, bus3_id = row.bus1_id, row.bus2_id, row.bus3_id + if bus1_id or bus2_id or bus3_id: + leg1_index = variable_context.t3_leg1_num_2_index[t3_num] + leg2_index = variable_context.t3_leg2_num_2_index[t3_num] + leg3_index = variable_context.t3_leg3_num_2_index[t3_num] + self._create_leg_constraint(leg1_r, leg1_x, leg1_g, leg1_b, rho1, alpha1, bus1_id, t3_index, leg1_index, parameters, network_cache, variable_context, model) + self._create_leg_constraint(leg2_r, leg2_x, leg2_g, leg2_b, rho2, alpha2, bus2_id, t3_index, leg2_index, parameters, network_cache, variable_context, model) + self._create_leg_constraint(leg3_r, leg3_x, leg3_g, leg3_b, rho3, alpha3, bus3_id, t3_index, leg3_index, parameters, network_cache, variable_context, model) diff --git a/pypowsybl/opf/impl/costs/__init__.py b/pypowsybl/opf/impl/costs/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pypowsybl/opf/impl/costs/minimal_active_power_cost_function.py b/pypowsybl/opf/impl/costs/minimal_active_power_cost_function.py new file mode 100644 index 0000000000..466b16b835 --- /dev/null +++ b/pypowsybl/opf/impl/costs/minimal_active_power_cost_function.py @@ -0,0 +1,18 @@ +import pyoptinterface as poi +from pyoptinterface import ExprBuilder + +from pypowsybl.opf.impl.model.cost_function import CostFunction +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class MinimalActivePowerCostFunction(CostFunction): + def __init__(self): + super().__init__('Minimal active power') + + def create(self, network_cache: NetworkCache, variable_context: VariableContext) -> ExprBuilder: + cost = poi.ExprBuilder() + for gen_num in range(len(variable_context.gen_p_vars)): + a, b, c = 0, 1.0, 0 # TODO + cost += a * variable_context.gen_p_vars[gen_num] * variable_context.gen_p_vars[gen_num] + b * variable_context.gen_p_vars[gen_num] + c + return cost diff --git a/pypowsybl/opf/impl/costs/minimal_losses_cost_function.py b/pypowsybl/opf/impl/costs/minimal_losses_cost_function.py new file mode 100644 index 0000000000..d088466a12 --- /dev/null +++ b/pypowsybl/opf/impl/costs/minimal_losses_cost_function.py @@ -0,0 +1,18 @@ +import pyoptinterface as poi +from pyoptinterface import ExprBuilder + +from pypowsybl.opf.impl.model.cost_function import CostFunction +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class MinimalLossesCostFunction(CostFunction): + def __init__(self): + super().__init__('Minimal losses power') + + def create(self, network_cache: NetworkCache, variable_context: VariableContext) -> ExprBuilder: + cost = poi.ExprBuilder() + for branch_index in range(len(variable_context.closed_branch_p1_vars)): + cost += variable_context.closed_branch_p1_vars[branch_index] - variable_context.closed_branch_p2_vars[branch_index] + return cost + diff --git a/pypowsybl/opf/impl/costs/minimize_against_reference_cost_function.py b/pypowsybl/opf/impl/costs/minimize_against_reference_cost_function.py new file mode 100644 index 0000000000..9f2f0d380c --- /dev/null +++ b/pypowsybl/opf/impl/costs/minimize_against_reference_cost_function.py @@ -0,0 +1,51 @@ +import pyoptinterface as poi +from pandas import DataFrame +from pyoptinterface import ExprBuilder + +from pypowsybl.opf.impl.model.cost_function import CostFunction +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class MinimizeAgainstReferenceCostFunction(CostFunction): + def __init__(self): + super().__init__('Minimize against reference') + + def create(self, network_cache: NetworkCache, variable_context: VariableContext) -> ExprBuilder: + cost = poi.ExprBuilder() + + for gen_num, gen_row in enumerate(network_cache.generators.itertuples(index=False)): + if gen_row.bus_id: + gen_p_expr = poi.ExprBuilder() + gen_p_expr += variable_context.gen_p_vars[gen_num] + gen_p_expr += gen_row.target_p + cost += gen_p_expr * gen_p_expr + if gen_row.voltage_regulator_on: + bus_num = network_cache.buses.index.get_loc(gen_row.bus_id) + v_var = variable_context.v_vars[bus_num] + cost += (v_var - gen_row.target_v) * (v_var - gen_row.target_v) + + for bat_num, bat_row in enumerate(network_cache.batteries.itertuples(index=False)): + if bat_row.bus_id: + bat_p_expr = poi.ExprBuilder() + bat_p_expr += variable_context.bat_p_vars[bat_num] + bat_p_expr += bat_row.target_p + cost += bat_p_expr * bat_p_expr + if bat_row.voltage_regulator_on: + bus_num = network_cache.buses.index.get_loc(bat_row.bus_id) + v_var = variable_context.v_vars[bus_num] + cost += (v_var - bat_row.target_v) * (v_var - bat_row.target_v) + + for vsc_cs_num, vsc_cs_row in enumerate(network_cache.vsc_converter_stations.itertuples()): + if vsc_cs_row.bus_id: + if NetworkCache.is_rectifier(vsc_cs_row.Index, vsc_cs_row): + vsc_cs_p_expr = poi.ExprBuilder() + vsc_cs_p_expr += variable_context.vsc_cs_p_vars[vsc_cs_num] + vsc_cs_p_expr -= vsc_cs_row.target_p + cost += vsc_cs_p_expr * vsc_cs_p_expr + if vsc_cs_row.voltage_regulator_on: + bus_num = network_cache.buses.index.get_loc(vsc_cs_row.bus_id) + v_var = variable_context.v_vars[bus_num] + cost += (v_var - vsc_cs_row.target_v) * (v_var - vsc_cs_row.target_v) + + return cost diff --git a/pypowsybl/opf/impl/costs/redispatching_cost_function.py b/pypowsybl/opf/impl/costs/redispatching_cost_function.py new file mode 100644 index 0000000000..faffea5a83 --- /dev/null +++ b/pypowsybl/opf/impl/costs/redispatching_cost_function.py @@ -0,0 +1,43 @@ +import pyoptinterface as poi +from pyoptinterface import ExprBuilder + +from pypowsybl.opf.impl.model.cost_function import CostFunction +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class RedispatchingCostFunction(CostFunction): + def __init__(self, target_p_l1_weight=0.1, target_p_l2_weight=1.0, target_v_l2_weight=0.5): + super().__init__('Redispatching') + self._target_p_l1_weight = target_p_l1_weight + self._target_p_l2_weight = target_p_l2_weight + self._target_v_l2_weight = target_v_l2_weight + + def create(self, network_cache: NetworkCache, variable_context: VariableContext) -> ExprBuilder: + cost = poi.ExprBuilder() + for gen_num, gen_row in enumerate(network_cache.generators.itertuples(index=False)): + if gen_row.bus_id: + gen_p_expr = poi.ExprBuilder() + gen_p_expr += variable_context.gen_p_vars[gen_num] + gen_p_expr += gen_row.target_p + cost += self._target_p_l2_weight * gen_p_expr * gen_p_expr + cost += self._target_p_l1_weight * gen_p_expr + if gen_row.voltage_regulator_on: + bus_num = network_cache.buses.index.get_loc(gen_row.bus_id) + v_var = variable_context.v_vars[bus_num] + cost += self._target_v_l2_weight * (v_var - gen_row.target_v) * (v_var - gen_row.target_v) + + for vsc_cs_num, vsc_cs_row in enumerate(network_cache.vsc_converter_stations.itertuples()): + if vsc_cs_row.bus_id: + if NetworkCache.is_rectifier(vsc_cs_row.Index, vsc_cs_row): + vsc_cs_p_expr = poi.ExprBuilder() + vsc_cs_p_expr += variable_context.vsc_cs_p_vars[vsc_cs_num] + vsc_cs_p_expr -= vsc_cs_row.target_p + cost += self._target_p_l2_weight * vsc_cs_p_expr * vsc_cs_p_expr + cost += self._target_p_l1_weight * vsc_cs_p_expr + if vsc_cs_row.voltage_regulator_on: + bus_num = network_cache.buses.index.get_loc(vsc_cs_row.bus_id) + v_var = variable_context.v_vars[bus_num] + cost += self._target_v_l2_weight * (v_var - vsc_cs_row.target_v) * (v_var - vsc_cs_row.target_v) + + return cost diff --git a/pypowsybl/opf/impl/model/__init__.py b/pypowsybl/opf/impl/model/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pypowsybl/opf/impl/model/bounds.py b/pypowsybl/opf/impl/model/bounds.py new file mode 100644 index 0000000000..259bbc5d1e --- /dev/null +++ b/pypowsybl/opf/impl/model/bounds.py @@ -0,0 +1,56 @@ +import logging +from typing import Any + +logger = logging.getLogger(__name__) + + +class Bounds: + def __init__(self, min_value: float, max_value: float, margin: float = 1e-6): + self._min_value = min_value + self._max_value = max_value + self._margin = margin + + @property + def min_value(self) -> float: + return self._min_value + + @property + def min_value_with_margin(self) -> float: + return self._min_value - self._margin + + @property + def max_value_with_margin(self) -> float: + return self._max_value + self._margin + + @property + def max_value(self) -> float: + return self._max_value + + def contains(self, value: float) -> bool: + return self.min_value_with_margin <= value <= self.max_value_with_margin + + def reduce(self, reduction: float) -> 'Bounds': + min_sign = 1 if self._min_value >= 0 else -1 + max_sign = 1 if self._max_value >= 0 else -1 + return Bounds(self._min_value * (1 + min_sign * reduction), self._max_value * (1 - max_sign * reduction)) + + def mirror(self) -> 'Bounds': + return Bounds(-self._max_value, -self._min_value) + + def __repr__(self) -> str: + return f"[{self._min_value}, {self._max_value}]" + + @staticmethod + def get_voltage_bounds(_low_voltage_limit: float | None, _high_voltage_limit: float | None, default_voltage_bounds: 'Bounds'): + return default_voltage_bounds # FIXME get from voltage level dataframe + + @staticmethod + def get_reactive_power_bounds(row: Any) -> 'Bounds': + return Bounds(row.min_q_at_target_p, row.max_q_at_target_p) + + @staticmethod + def fix(id:str, lb: float, ub: float) -> tuple[float, float]: + if lb > ub: + logger.warning(f"{id}, lower bound {lb} is greater than upper bound {ub}") + return ub, lb + return lb, ub diff --git a/pypowsybl/opf/impl/model/constraints.py b/pypowsybl/opf/impl/model/constraints.py new file mode 100644 index 0000000000..d22158ed0d --- /dev/null +++ b/pypowsybl/opf/impl/model/constraints.py @@ -0,0 +1,17 @@ +from abc import ABC, abstractmethod + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + + +class Constraints(ABC): + @abstractmethod + def add(self, + parameters: ModelParameters, + network_cache: NetworkCache, + variable_context: VariableContext, + model: ipopt.Model) -> None: + pass diff --git a/pypowsybl/opf/impl/model/cost_function.py b/pypowsybl/opf/impl/model/cost_function.py new file mode 100644 index 0000000000..b64447322b --- /dev/null +++ b/pypowsybl/opf/impl/model/cost_function.py @@ -0,0 +1,19 @@ +from abc import ABC, abstractmethod + +from pyoptinterface import ExprBuilder + +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.variable_context import VariableContext + + +class CostFunction(ABC): + def __init__(self, name): + self._name = name + + @property + def name(self): + return self._name + + @abstractmethod + def create(self, network_cache: NetworkCache, variable_context: VariableContext) -> ExprBuilder: + pass diff --git a/pypowsybl/opf/impl/model/model_parameters.py b/pypowsybl/opf/impl/model/model_parameters.py new file mode 100644 index 0000000000..ae4545dc65 --- /dev/null +++ b/pypowsybl/opf/impl/model/model_parameters.py @@ -0,0 +1,23 @@ +from pypowsybl.opf.impl.model.bounds import Bounds + + +class ModelParameters: + def __init__(self, + reactive_bounds_reduction: float, + twt_split_shunt_admittance: bool, + default_voltage_bounds: Bounds) -> None: + self._reactive_bounds_reduction = reactive_bounds_reduction + self._twt_split_shunt_admittance = twt_split_shunt_admittance + self._default_voltage_bounds = default_voltage_bounds + + @property + def reactive_bounds_reduction(self) -> float: + return self._reactive_bounds_reduction + + @property + def twt_split_shunt_admittance(self) -> bool: + return self._twt_split_shunt_admittance + + @property + def default_voltage_bounds(self) -> Bounds: + return self._default_voltage_bounds \ No newline at end of file diff --git a/pypowsybl/opf/impl/model/network_cache.py b/pypowsybl/opf/impl/model/network_cache.py new file mode 100644 index 0000000000..f7cabff7fe --- /dev/null +++ b/pypowsybl/opf/impl/model/network_cache.py @@ -0,0 +1,462 @@ +import math + +import pandas as pd +from pandas import DataFrame + +from pypowsybl.network import Network +from pypowsybl.opf.impl.model.bounds import Bounds + + +class NetworkCache: + def __init__(self, network: Network) -> None: + self._network = network + self._network.per_unit = True + self._voltage_levels = self._build_voltage_levels(network) + self._buses = self._build_buses(network, self._voltage_levels) + self._generators = self._build_generators(network, self.buses) + self._loads = self._build_loads(network, self.buses) + self._shunts = self._build_shunt_compensators(network, self.buses) + self._static_var_compensators = self._build_static_var_compensators(network, self.buses) + self._vsc_converter_stations = self._build_vsc_converter_stations(network, self.buses) + self._lcc_converter_stations = self._build_lcc_converter_stations(network, self.buses) + self._hvdc_lines = self._build_hvdc_lines(network, self._vsc_converter_stations, self._lcc_converter_stations, self.buses) + self._lines = self._build_lines(network, self.buses) + self._transformers_2w = self._build_2w_transformers(network, self.buses) + self._transformers_3w = self._build_3w_transformers(network, self.buses) + self._branches = self._build_branches(network, self.buses) + self._dangling_lines = self._build_dangling_lines(network, self.buses) + self._batteries = self._build_batteries(network, self.buses) + self._current_limits1, self._current_limits2 = self._build_current_limits(network) + self._slack_terminal = self._build_stack_terminal(network, self.buses) + + @staticmethod + def _filter_injections(injections: DataFrame, buses: DataFrame) -> DataFrame: + if len(injections) == 0: + return injections + injections_and_buses = pd.merge(injections, buses, right_index=True, left_on='bus_id', how='left') + return injections_and_buses[ + (injections_and_buses['connected_component'] == 0) & (injections_and_buses['synchronous_component'] == 0)] + + @staticmethod + def _filter_branches(branches: DataFrame, buses: DataFrame) -> DataFrame: + if len(branches) == 0: + return branches + branches_and_buses = pd.merge(branches, buses, left_on='bus1_id', right_index=True, how='left') + branches_and_buses = pd.merge(branches_and_buses, buses, left_on='bus2_id', right_index=True, + suffixes=('', '_2'), how='left') + return branches_and_buses[ + (branches_and_buses['connected_component'] == 0) & (branches_and_buses['synchronous_component'] == 0) | ( + branches_and_buses['connected_component_2'] == 0) & ( + branches_and_buses['synchronous_component_2'] == 0)] + + @staticmethod + def _filter_3w_transformers(transformers_3w: DataFrame, buses: DataFrame) -> DataFrame: + if len(transformers_3w) == 0: + return transformers_3w + transfos_and_buses = pd.merge(transformers_3w, buses, left_on='bus1_id', right_index=True, how='left') + transfos_and_buses = pd.merge(transfos_and_buses, buses, left_on='bus2_id', right_index=True, suffixes=('', '_2'), how='left') + transfos_and_buses = pd.merge(transfos_and_buses, buses, left_on='bus3_id', right_index=True, suffixes=('', '_3'), how='left') + return transfos_and_buses[ + (transfos_and_buses['connected_component'] == 0) & (transfos_and_buses['synchronous_component'] == 0) + | (transfos_and_buses['connected_component_2'] == 0) & (transfos_and_buses['synchronous_component_2'] == 0) + | (transfos_and_buses['connected_component_3'] == 0) & (transfos_and_buses['synchronous_component_3'] == 0)] + + @staticmethod + def _build_branches(network: Network, buses: DataFrame): + branches = network.get_branches(attributes=['bus1_id', 'bus2_id']) + return NetworkCache._filter_branches(branches, buses) + + @staticmethod + def _build_2w_transformers(network: Network, buses: DataFrame): + transfos = network.get_2_windings_transformers( + attributes=['bus1_id', 'bus2_id', 'rho', 'alpha', 'r_at_current_tap', 'x_at_current_tap', 'g_at_current_tap', 'b_at_current_tap']) + return NetworkCache._filter_branches(transfos, buses) + + @staticmethod + def _build_3w_transformers(network: Network, buses: DataFrame): + transfos = network.get_3_windings_transformers( + attributes=['bus1_id', 'bus2_id', 'bus3_id', + 'rated_u0', + 'rho1', 'rho2', 'rho3', + 'alpha1', 'alpha2', 'alpha3', + 'r1_at_current_tap', 'r2_at_current_tap', 'r3_at_current_tap', 'x1_at_current_tap', 'x2_at_current_tap', 'x3_at_current_tap', + 'g1_at_current_tap', 'g2_at_current_tap', 'g3_at_current_tap', 'b1_at_current_tap', 'b2_at_current_tap', 'b3_at_current_tap']) + return NetworkCache._filter_3w_transformers(transfos, buses) + + @staticmethod + def _build_lines(network: Network, buses: DataFrame): + lines = network.get_lines(attributes=['bus1_id', 'bus2_id', 'r', 'x', 'g1', 'g2', 'b1', 'b2']) + return NetworkCache._filter_branches(lines, buses) + + @staticmethod + def _build_shunt_compensators(network: Network, buses: DataFrame): + shunts = network.get_shunt_compensators(attributes=['bus_id', 'g', 'b']) + return NetworkCache._filter_injections(shunts, buses) + + @staticmethod + def _build_static_var_compensators(network: Network, buses: DataFrame): + svcs = network.get_static_var_compensators(attributes=['bus_id', 'b_min', 'b_max', 'regulated_bus_id']) + return NetworkCache._filter_injections(svcs, buses) + + @staticmethod + def _filter_hvdc_lines(hvdc_lines: DataFrame, vsc_converter_stations: DataFrame, + lcc_converter_stations: DataFrame, buses: DataFrame) -> DataFrame: + if len(hvdc_lines) == 0: + return hvdc_lines + converter_stations = pd.concat([vsc_converter_stations[['bus_id']], lcc_converter_stations[['bus_id']]]) + hvdc_lines_and_buses = pd.merge(hvdc_lines, converter_stations, left_on='converter_station1_id', + right_index=True, how='left') + hvdc_lines_and_buses = pd.merge(hvdc_lines_and_buses, converter_stations, left_on='converter_station2_id', + right_index=True, suffixes=('', '_2'), how='left') + hvdc_lines_and_buses = pd.merge(hvdc_lines_and_buses, buses, left_on='bus_id', right_index=True, how='left') + hvdc_lines_and_buses = pd.merge(hvdc_lines_and_buses, buses, left_on='bus_id_2', right_index=True, + suffixes=('', '_2'), how='left') + return hvdc_lines_and_buses[ + (hvdc_lines_and_buses['connected_component'] == 0) & (hvdc_lines_and_buses['synchronous_component'] == 0) & ( + hvdc_lines_and_buses['connected_component_2'] == 0) & ( + hvdc_lines_and_buses['synchronous_component_2'] == 0)] + + @staticmethod + def _build_hvdc_lines(network: Network, vsc_converter_stations: DataFrame, lcc_converter_stations: DataFrame, + buses: DataFrame): + hvdc_lines = network.get_hvdc_lines(attributes=['converter_station1_id', 'converter_station2_id', 'converters_mode', 'nominal_v', 'r']) + return NetworkCache._filter_hvdc_lines(hvdc_lines, vsc_converter_stations, lcc_converter_stations, buses) + + @staticmethod + def _build_vsc_converter_stations(network: Network, buses: DataFrame): + stations = network.get_vsc_converter_stations(attributes=['bus_id', 'loss_factor', 'min_q_at_target_p', 'max_q_at_target_p', + 'target_v', 'target_q', 'voltage_regulator_on', 'hvdc_line_id', + 'regulated_bus_id']) + if len(stations) > 0: + hvdc_lines = network.get_hvdc_lines(attributes=['converter_station1_id', 'converter_station2_id', 'converters_mode', 'target_p', 'max_p']) + stations = pd.merge(stations, hvdc_lines, left_on='hvdc_line_id', right_index=True, how='left') + return NetworkCache._filter_injections(stations, buses) + + @staticmethod + def _build_lcc_converter_stations(network: Network, buses: DataFrame): + stations = network.get_lcc_converter_stations(attributes=['bus_id', 'loss_factor', 'power_factor', 'hvdc_line_id']) + if len(stations) > 0: + hvdc_lines = network.get_hvdc_lines(attributes=['converters_mode', 'target_p', 'max_p']) + stations = pd.merge(stations, hvdc_lines, left_on='hvdc_line_id', right_index=True, how='left') + return NetworkCache._filter_injections(stations, buses) + + @staticmethod + def _build_loads(network: Network, buses: DataFrame): + loads = network.get_loads(attributes=['bus_id', 'p0', 'q0']) + return NetworkCache._filter_injections(loads, buses) + + @staticmethod + def _build_generators(network: Network, buses: DataFrame): + generators = network.get_generators( + attributes=['bus_id', 'min_p', 'max_p', 'min_q_at_target_p', 'max_q_at_target_p', 'target_p', 'target_q', + 'target_v', 'voltage_regulator_on', 'regulated_bus_id']) + return NetworkCache._filter_injections(generators, buses) + + @staticmethod + def _build_buses(network: Network, voltage_levels: DataFrame): + buses = network.get_buses(attributes=['voltage_level_id', 'connected_component', 'synchronous_component']) + buses = buses[(buses['synchronous_component'] == 0) & (buses['connected_component'] == 0)] + return pd.merge(buses, voltage_levels, left_on='voltage_level_id', right_index=True, how='left') + + @staticmethod + def _build_voltage_levels(network: Network) -> DataFrame: + return network.get_voltage_levels(attributes=['low_voltage_limit', 'high_voltage_limit', 'nominal_v']) + + @staticmethod + def _build_dangling_lines(network: Network, buses: DataFrame): + dangling_lines = network.get_dangling_lines(attributes=['bus_id', 'r', 'x', 'g', 'b', 'p0', 'q0', 'paired']) + return NetworkCache._filter_injections(dangling_lines, buses) + + @staticmethod + def _build_batteries(network: Network, buses: DataFrame): + batteries = network.get_batteries(attributes=['bus_id', 'min_p', 'max_p', 'min_q_at_target_p', 'max_q_at_target_p', 'target_p', 'target_q']) + if len(batteries): + voltage_regulation = network.get_extensions('voltageRegulation') + batteries = pd.merge(batteries, voltage_regulation, left_index=True, right_on='id', how='left') + batteries['voltage_regulator_on'] = batteries['voltage_regulator_on'].fillna(False) + batteries = NetworkCache._filter_injections(batteries, buses) + # FIXME to remove when extensions will be per united + batteries['target_v'] /= batteries['nominal_v'] + return batteries + + @staticmethod + def _build_current_limits(network: Network) -> tuple[DataFrame, DataFrame]: + limits = network.get_operational_limits(attributes=['name', 'value']).reset_index( + level=['side', 'type', 'acceptable_duration', 'group_name']) + limits = limits[(limits['type'] == 'CURRENT') & (limits['name'] == 'permanent_limit')] + return limits[limits['side'] == 'ONE'][['value']], limits[limits['side'] == 'TWO'][['value']] + + @staticmethod + def _build_stack_terminal(network, buses: DataFrame): + slack_terminal = network.get_extensions('slackTerminal') + return NetworkCache._filter_injections(slack_terminal, buses) + + @property + def network(self) -> Network: + return self._network + + @property + def buses(self) -> DataFrame: + return self._buses + + @property + def generators(self) -> DataFrame: + return self._generators + + @property + def loads(self) -> DataFrame: + return self._loads + + @property + def shunts(self) -> DataFrame: + return self._shunts + + @property + def static_var_compensators(self) -> DataFrame: + return self._static_var_compensators + + @property + def vsc_converter_stations(self) -> DataFrame: + return self._vsc_converter_stations + + @property + def lcc_converter_stations(self) -> DataFrame: + return self._lcc_converter_stations + + @property + def lines(self) -> DataFrame: + return self._lines + + @property + def transformers_2w(self) -> DataFrame: + return self._transformers_2w + + @property + def transformers_3w(self) -> DataFrame: + return self._transformers_3w + + @property + def branches(self) -> DataFrame: + return self._branches + + @property + def hvdc_lines(self) -> DataFrame: + return self._hvdc_lines + + @property + def dangling_lines(self) -> DataFrame: + return self._dangling_lines + + @property + def batteries(self) -> DataFrame: + return self._batteries + + @property + def current_limits1(self) -> DataFrame: + return self._current_limits1 + + @property + def current_limits2(self) -> DataFrame: + return self._current_limits2 + + @property + def slack_terminal(self) -> DataFrame: + return self._slack_terminal + + @staticmethod + def is_rectifier(vsc_cs_id: str, hvdc_line_row) -> bool: + return ((vsc_cs_id == hvdc_line_row.converter_station1_id and hvdc_line_row.converters_mode == 'SIDE_1_RECTIFIER_SIDE_2_INVERTER') + or (vsc_cs_id == hvdc_line_row.converter_station2_id and hvdc_line_row.converters_mode == 'SIDE_1_INVERTER_SIDE_2_RECTIFIER')) + + def is_too_small_reactive_power_bounds(self, q_bounds: Bounds): + return abs(q_bounds.max_value - q_bounds.min_value) < 1.0 / self._network.nominal_apparent_power + + def update_generators(self, + connected_gen_ids: list[str], + connected_gen_target_p: list[float], + connected_gen_target_q: list[float], + connected_gen_target_v: list[float], + connected_gen_voltage_regulator_on: list[bool], + connected_gen_p: list[float], + connected_gen_q: list[float], + disconnected_gen_ids: list[str] = None): + + if connected_gen_ids: + self._network.update_generators(id=connected_gen_ids, + target_p=connected_gen_target_p, + target_q=connected_gen_target_q, + target_v=connected_gen_target_v, + voltage_regulator_on=connected_gen_voltage_regulator_on, + p=connected_gen_p, + q=connected_gen_q) + + if disconnected_gen_ids: + self._network.update_generators(id=disconnected_gen_ids, + p=[0.0] * len(disconnected_gen_ids), + q=[0.0] * len(disconnected_gen_ids)) + + self._generators = self._build_generators(self._network, self._buses) + + def update_batteries(self, + connected_bat_ids: list[str], + connected_bat_target_p: list[float], + connected_bat_target_q: list[float], + connected_bat_target_v: list[float], + connected_bat_voltage_regulator_on: list[bool], + connected_bat_p: list[float], + connected_bat_q: list[float], + disconnected_bat_ids: list[str] = None): + if connected_bat_ids: + self._network.update_batteries(id=connected_bat_ids, + target_p=connected_bat_target_p, + target_q=connected_bat_target_q, + p=connected_bat_p, + q=connected_bat_q) + self._network.update_extensions("voltageRegulation", + id=connected_bat_ids, + voltage_regulator_on=connected_bat_voltage_regulator_on, + target_v=connected_bat_target_v) + + if disconnected_bat_ids: + self._network.update_batteries(id=disconnected_bat_ids, + p=[0.0] * len(disconnected_bat_ids), + q=[0.0] * len(disconnected_bat_ids)) + + self._batteries = self._build_batteries(self._network, self._buses) + + def update_vsc_converter_stations(self, + connected_vsc_cs_ids: list[str], + connected_vsc_cs_target_q: list[float], + connected_vsc_cs_target_v: list[float], + connected_vsc_cs_voltage_regulator_on: list[bool], + connected_vsc_cs_p: list[float], + connected_vsc_cs_q: list[float], + disconnected_vsc_cs_ids: list[str]): + + if connected_vsc_cs_ids: + self._network.update_vsc_converter_stations(id=connected_vsc_cs_ids, + target_q=connected_vsc_cs_target_q, + target_v=connected_vsc_cs_target_v, + voltage_regulator_on=connected_vsc_cs_voltage_regulator_on, + p=connected_vsc_cs_p, + q=connected_vsc_cs_q) + + if disconnected_vsc_cs_ids: + self._network.update_vsc_converter_stations(id=disconnected_vsc_cs_ids, + p=[0.0] * len(disconnected_vsc_cs_ids), + q=[0.0] * len(disconnected_vsc_cs_ids)) + + self._vsc_converter_stations = self._build_vsc_converter_stations(self._network, self._buses) + + def update_hvdc_lines(self, hvdc_line_ids: list[str], hvdc_line_target_p: list[float]): + self._network.update_hvdc_lines(id=hvdc_line_ids, target_p=hvdc_line_target_p) + self._hvdc_lines = self._build_hvdc_lines(self._network, self._vsc_converter_stations, self._lcc_converter_stations, self._buses) + + def update_static_var_compensators(self, + connected_svc_ids: list[str], + connected_svc_target_q: list[float], + connected_svc_target_v: list[float], + connected_svc_regulation_mode: list[str], + connected_svc_p: list[float], + connected_svc_q: list[float], + disconnected_svc_ids: list[str]): + if connected_svc_ids: + self._network.update_static_var_compensators(id=connected_svc_ids, + target_q=connected_svc_target_q, + target_v=connected_svc_target_v, + regulation_mode=connected_svc_regulation_mode, + p=connected_svc_p, + q=connected_svc_q) + + if disconnected_svc_ids: + self._network.update_static_var_compensators(id=disconnected_svc_ids, + p=[0.0] * len(disconnected_svc_ids), + q=[0.0] * len(disconnected_svc_ids)) + + self._static_var_compensators = self._build_static_var_compensators(self._network, self._buses) + + def update_shunt_compensators(self, + connected_shunt_ids: list[str], + connected_shunt_p: list[float], + connected_shunt_q: list[float], + disconnected_shunt_ids: list[str]): + if connected_shunt_ids: + self._network.update_shunt_compensators(id=connected_shunt_ids, + p=connected_shunt_p, + q=connected_shunt_q) + + if disconnected_shunt_ids: + self._network.update_shunt_compensators(id=disconnected_shunt_ids, + p=[0.0] * len(disconnected_shunt_ids), + q=[0.0] * len(disconnected_shunt_ids)) + + self._shunts = self._build_shunt_compensators(self._network, self._buses) + + def update_branches(self, + branch_ids: list[str], + branch_p1: list[float], + branch_p2: list[float], + branch_q1: list[float], + branch_q2: list[float]): + self._network.update_branches(id=branch_ids, + p1=branch_p1, + p2=branch_p2, + q1=branch_q1, + q2=branch_q2) + + self._branches = self._build_branches(self._network, self._buses) + + def update_transformers_3w(self, + t3_ids: list[str], + t3_p1: list[float], + t3_p2: list[float], + t3_p3: list[float], + t3_q1: list[float], + t3_q2: list[float], + t3_q3: list[float], + t3_v: list[float], + t3_angle: list[float]): + self._network.update_3_windings_transformers(id=t3_ids, + p1=t3_p1, + p2=t3_p2, + p3=t3_p3, + q1=t3_q1, + q2=t3_q2, + q3=t3_q3) + self._network.add_elements_properties(id=t3_ids, + v=t3_v, + angle=t3_angle) + + self._transformers_3w = self._build_3w_transformers(self._network, self._buses) + + def update_buses(self, + bus_ids: list[str], + bus_v_mag: list[float], + bus_v_angle: list[float]): + self._network.update_buses(id=bus_ids, v_mag=bus_v_mag, v_angle=bus_v_angle) + self._buses = self._build_buses(self._network, self._voltage_levels) + + def update_dangling_lines(self, + connected_dl_ids: list[str], + connected_dl_v: list[float], + connected_dl_angle: list[float], + connected_dl_p: list[float], + connected_dl_q: list[float], + disconnected_dl_ids: list[str]): + if connected_dl_ids: + self._network.update_dangling_lines(id=connected_dl_ids, + p=connected_dl_p, + q=connected_dl_q) + self._network.add_elements_properties(id=connected_dl_ids, + v=connected_dl_v, + angle=connected_dl_angle) + + if disconnected_dl_ids: + self._network.update_dangling_lines(id=disconnected_dl_ids, + p=[0.0] * len(disconnected_dl_ids), + q=[0.0] * len(disconnected_dl_ids)) + self._network.add_elements_properties(id=disconnected_dl_ids, + v=[math.nan] * len(disconnected_dl_ids), + angle=[math.nan] * len(disconnected_dl_ids)) + + self._dangling_lines = self._build_dangling_lines(self._network, self._buses) diff --git a/pypowsybl/opf/impl/model/opf_model.py b/pypowsybl/opf/impl/model/opf_model.py new file mode 100644 index 0000000000..ba551f1719 --- /dev/null +++ b/pypowsybl/opf/impl/model/opf_model.py @@ -0,0 +1,91 @@ +import logging +import time + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.cost_function import CostFunction +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_bounds import VariableBounds +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache + +logger = logging.getLogger(__name__) + +class OpfModel: + def __init__(self, network_cache: NetworkCache, model: ipopt.Model, variable_context: VariableContext): + self._network_cache = network_cache + self._model = model + self._variable_context = variable_context + + @property + def network_cache(self) -> NetworkCache: + return self._network_cache + + @property + def model(self) -> ipopt.Model: + return self._model + + @property + def variable_context(self) -> VariableContext: + return self._variable_context + + @classmethod + def build(cls, network_cache: NetworkCache, parameters: ModelParameters, + variable_bounds: list[VariableBounds], constraints: list[Constraints], + cost_function: CostFunction) -> 'OpfModel': + logger.info("Building model...") + start = time.perf_counter() + + model = ipopt.Model() + model.set_raw_parameter("print_user_options", "yes") + model.set_raw_parameter("print_advanced_options", "yes") + + # create variables + variable_context = VariableContext.build(network_cache, model) + + # variable bounds + for variable_bounds in variable_bounds: + variable_bounds.add(parameters, network_cache, variable_context, model) + + # constraints + for constraint in constraints: + constraint.add(parameters, network_cache, variable_context, model) + + # cost function + logger.debug(f"Using cost function: '{cost_function.name}'") + cost = cost_function.create(network_cache, variable_context) + model.set_objective(cost) + + logger.info(f"Model built in {time.perf_counter() - start:.3f} seconds") + + return OpfModel(network_cache, model, variable_context) + + def update_network(self): + self.variable_context.update_network(self.network_cache, self.model) + + def analyze_violations(self, parameters: ModelParameters) -> None: + # check voltage bounds + for bus_num, (bus_id, row) in enumerate(self.network_cache.buses.iterrows()): + v = self.model.get_value(self.variable_context.v_vars[bus_num]) + v_bounds = Bounds.get_voltage_bounds(row.low_voltage_limit, row.high_voltage_limit, parameters.default_voltage_bounds) + if not v_bounds.contains(v): + logger.error(f"Voltage magnitude violation: bus '{bus_id}' (num={bus_num}) {v} not in {v_bounds}") + + # check generator limits + for gen_num, (gen_id, row) in enumerate(self.network_cache.generators.iterrows()): + if row.bus_id: + gen_p_index = self.variable_context.gen_p_num_2_index[gen_num] + p = self.model.get_value(self.variable_context.gen_p_vars[gen_p_index]) + + p_bounds = Bounds(row.min_p, row.max_p).mirror() + if not p_bounds.contains(p): + logger.error(f"Generator active power violation: generator '{gen_id}' (num={gen_num}) {p} not in [{-row.max_p}, {-row.min_p}]") + + gen_q_index = self.variable_context.gen_q_num_2_index[gen_num] + if gen_q_index != -1: # valid + q = self.model.get_value(self.variable_context.gen_q_vars[gen_q_index]) + q_bounds = Bounds.get_reactive_power_bounds(row).mirror() + if not q_bounds.contains(q): + logger.error(f"Generator reactive power violation: generator '{gen_id}' (num={gen_num}) {q} not in {q_bounds}") diff --git a/pypowsybl/opf/impl/model/variable_bounds.py b/pypowsybl/opf/impl/model/variable_bounds.py new file mode 100644 index 0000000000..775e8ab38e --- /dev/null +++ b/pypowsybl/opf/impl/model/variable_bounds.py @@ -0,0 +1,17 @@ +from abc import ABC, abstractmethod + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.variable_context import VariableContext +from pypowsybl.opf.impl.model.network_cache import NetworkCache + + +class VariableBounds(ABC): + @abstractmethod + def add(self, + parameters: ModelParameters, + network_cache: NetworkCache, + variable_context: VariableContext, + model: ipopt.Model): + pass diff --git a/pypowsybl/opf/impl/model/variable_context.py b/pypowsybl/opf/impl/model/variable_context.py new file mode 100644 index 0000000000..6a3087cd9a --- /dev/null +++ b/pypowsybl/opf/impl/model/variable_context.py @@ -0,0 +1,652 @@ +import logging +import math +from dataclasses import dataclass +from typing import Any + +from pyoptinterface import ipopt + +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.util import TRACE_LEVEL + +logger = logging.getLogger(__name__) + + +@dataclass +class VariableContext: + v_vars: Any + ph_vars: Any + gen_p_vars: Any + gen_q_vars: Any + bat_p_vars: Any + bat_q_vars: Any + shunt_p_vars: Any + shunt_q_vars: Any + svc_q_vars: Any + vsc_cs_p_vars: Any + vsc_cs_q_vars: Any + closed_branch_p1_vars: Any + closed_branch_q1_vars: Any + closed_branch_p2_vars: Any + closed_branch_q2_vars: Any + open_side1_branch_p2_vars: Any + open_side1_branch_q2_vars: Any + open_side2_branch_p1_vars: Any + open_side2_branch_q1_vars: Any + dl_v_vars: Any + dl_ph_vars: Any + dl_branch_p1_vars: Any + dl_branch_p2_vars: Any + dl_branch_q1_vars: Any + dl_branch_q2_vars: Any + t3_middle_v_vars: Any + t3_middle_ph_vars: Any + t3_closed_branch_p1_vars: Any + t3_closed_branch_p2_vars: Any + t3_closed_branch_q1_vars: Any + t3_closed_branch_q2_vars: Any + t3_open_side1_branch_p2_vars: Any + t3_open_side1_branch_q2_vars: Any + branch_num_2_index: list[int] + gen_p_num_2_index: list[int] + gen_q_num_2_index: list[int] + bat_p_num_2_index: list[int] + bat_q_num_2_index: list[int] + shunt_num_2_index: list[int] + svc_num_2_index: list[int] + vsc_cs_num_2_index: list[int] + dl_num_2_index: list[int] + t3_num_2_index: list[int] + t3_leg1_num_2_index: list[int] + t3_leg2_num_2_index: list[int] + t3_leg3_num_2_index: list[int] + + @staticmethod + def build(network_cache: NetworkCache, model: ipopt.Model) -> 'VariableContext': + bus_count = len(network_cache.buses) + v_vars = model.add_m_variables(bus_count, name="v") + ph_vars = model.add_m_variables(bus_count, name="ph") + + gen_count = len(network_cache.generators) + gen_p_nums = [] + gen_q_nums = [] + gen_p_num_2_index = [-1] * gen_count + gen_q_num_2_index = [-1] * gen_count + too_small_q_bounds_generator_ids = [] + for gen_num, row in enumerate(network_cache.generators.itertuples()): + if row.bus_id: + gen_p_num_2_index[gen_num] = len(gen_p_nums) + gen_p_nums.append(gen_num) + + q_bounds = Bounds.get_reactive_power_bounds(row).mirror() + if network_cache.is_too_small_reactive_power_bounds(q_bounds): + too_small_q_bounds_generator_ids.append(row.Index) + logger.log(TRACE_LEVEL, f"Too small reactive power bounds {q_bounds} for generator '{row.Index}' (num={gen_num})") + else: + gen_q_num_2_index[gen_num] = len(gen_q_nums) + gen_q_nums.append(gen_num) + + gen_p_vars = model.add_m_variables(len(gen_p_nums), name="gen_p") + gen_q_vars = model.add_m_variables(len(gen_q_nums), name="gen_q") + + bat_count = len(network_cache.batteries) + bat_p_nums = [] + bat_q_nums = [] + bat_p_num_2_index = [-1] * bat_count + bat_q_num_2_index = [-1] * bat_count + for bat_num, row in enumerate(network_cache.batteries.itertuples()): + if row.bus_id: + bat_p_num_2_index[bat_num] = len(bat_p_nums) + bat_p_nums.append(bat_num) + + q_bounds = Bounds.get_reactive_power_bounds(row).mirror() + if network_cache.is_too_small_reactive_power_bounds(q_bounds): + too_small_q_bounds_generator_ids.append(row.Index) + logger.log(TRACE_LEVEL, f"Too small reactive power bounds {q_bounds} for battery '{row.Index}' (num={bat_num})") + else: + bat_q_num_2_index[bat_num] = len(bat_q_nums) + bat_q_nums.append(bat_num) + + bat_p_vars = model.add_m_variables(len(bat_p_nums), name="bat_p") + bat_q_vars = model.add_m_variables(len(bat_q_nums), name="bat_q") + + if too_small_q_bounds_generator_ids: + logger.warning(f"{len(too_small_q_bounds_generator_ids)} generators|batteries have too small reactive power bounds") + + shunt_nums = [] + shunt_count = len(network_cache.shunts) + shunt_num_2_index = [-1] * shunt_count + for shunt_num, row in enumerate(network_cache.shunts.itertuples()): + if row.bus_id: + shunt_num_2_index[shunt_num] = len(shunt_nums) + shunt_nums.append(shunt_num) + shunt_p_vars = model.add_m_variables(len(shunt_nums), name="shunt_p") + shunt_q_vars = model.add_m_variables(len(shunt_nums), name="shunt_q") + + svc_nums = [] + svc_count = len(network_cache.static_var_compensators) + svc_num_2_index = [-1] * svc_count + for svc_num, row in enumerate(network_cache.static_var_compensators.itertuples()): + if row.bus_id: + svc_num_2_index[svc_num] = len(svc_nums) + svc_nums.append(svc_num) + svc_q_vars = model.add_m_variables(len(svc_nums), name="svc_q") + + vsc_cs_nums = [] + vsc_cs_count = len(network_cache.vsc_converter_stations) + vsc_cs_num_2_index = [-1] * vsc_cs_count + for vsc_cs_num, row in enumerate(network_cache.vsc_converter_stations.itertuples()): + if row.bus_id: + vsc_cs_num_2_index[vsc_cs_num] = len(vsc_cs_nums) + vsc_cs_nums.append(vsc_cs_num) + vsc_cs_p_vars = model.add_m_variables(len(vsc_cs_nums), name="vsc_cs_p_vars") + vsc_cs_q_vars = model.add_m_variables(len(vsc_cs_nums), name="vsc_cs_q_vars") + + closed_branch_nums = [] + open_side1_branch_nums = [] + open_side2_branch_nums = [] + branch_count = len(network_cache.lines) + len(network_cache.transformers_2w) + branch_num_2_index = [-1] * branch_count + for branch_num, row in enumerate(network_cache.branches.itertuples(index=False)): + if row.bus1_id and row.bus2_id: + branch_num_2_index[branch_num] = len(closed_branch_nums) + closed_branch_nums.append(branch_num) + elif row.bus2_id: + branch_num_2_index[branch_num] = len(open_side1_branch_nums) + open_side1_branch_nums.append(branch_num) + elif row.bus1_id: + branch_num_2_index[branch_num] = len(open_side2_branch_nums) + open_side2_branch_nums.append(branch_num) + closed_branch_p1_vars = model.add_m_variables(len(closed_branch_nums), name='closed_branch_p1') + closed_branch_q1_vars = model.add_m_variables(len(closed_branch_nums), name='closed_branch_q1') + closed_branch_p2_vars = model.add_m_variables(len(closed_branch_nums), name='closed_branch_p2') + closed_branch_q2_vars = model.add_m_variables(len(closed_branch_nums), name='closed_branch_q2') + open_side1_branch_p2_vars = model.add_m_variables(len(open_side1_branch_nums), name='open_side1_branch_p2') + open_side1_branch_q2_vars = model.add_m_variables(len(open_side1_branch_nums), name='open_side1_branch_q2') + open_side2_branch_p1_vars = model.add_m_variables(len(open_side2_branch_nums), name='open_side2_branch_p1') + open_side2_branch_q1_vars = model.add_m_variables(len(open_side2_branch_nums), name='open_side2_branch_q1') + + dl_count = len(network_cache.dangling_lines) + dl_nums = [] + dl_num_2_index = [-1] * dl_count + for dl_num, row in enumerate(network_cache.dangling_lines.itertuples()): + if row.bus_id: + dl_num_2_index[dl_num] = len(dl_nums) + dl_nums.append(dl_num) + + dl_v_vars = model.add_m_variables(len(dl_nums), name="dl_v") + dl_ph_vars = model.add_m_variables(len(dl_nums), name="dl_ph") + dl_branch_p1_vars = model.add_m_variables(len(dl_nums), name="dl_branch_p1") + dl_branch_p2_vars = model.add_m_variables(len(dl_nums), name="dl_branch_p2") + dl_branch_q1_vars = model.add_m_variables(len(dl_nums), name="dl_branch_q1") + dl_branch_q2_vars = model.add_m_variables(len(dl_nums), name="dl_branch_q2") + + t3_count = len(network_cache.transformers_3w) + t3_nums = [] + t3_num_2_index = [-1] * t3_count + t3_closed_branch_leg_nums = [] + t3_open_side2_leg_nums = [] + t3_leg1_num_2_index = [-1] * t3_count + t3_leg2_num_2_index = [-1] * t3_count + t3_leg3_num_2_index = [-1] * t3_count + for t3_num, t3_row in enumerate(network_cache.transformers_3w.itertuples()): + if t3_row.bus1_id or t3_row.bus2_id or t3_row.bus3_id: + t3_num_2_index[t3_num] = len(t3_nums) + t3_nums.append(t3_num) + + if t3_row.bus1_id: + t3_leg1_num_2_index[t3_num] = len(t3_closed_branch_leg_nums) + t3_closed_branch_leg_nums.append(t3_num) + else: + t3_leg1_num_2_index[t3_num] = len(t3_open_side2_leg_nums) + t3_open_side2_leg_nums.append(t3_num) + + if t3_row.bus2_id: + t3_leg2_num_2_index[t3_num] = len(t3_closed_branch_leg_nums) + t3_closed_branch_leg_nums.append(t3_num) + else: + t3_leg2_num_2_index[t3_num] = len(t3_open_side2_leg_nums) + t3_open_side2_leg_nums.append(t3_num) + + if t3_row.bus3_id: + t3_leg3_num_2_index[t3_num] = len(t3_closed_branch_leg_nums) + t3_closed_branch_leg_nums.append(t3_num) + else: + t3_leg3_num_2_index[t3_num] = len(t3_open_side2_leg_nums) + t3_open_side2_leg_nums.append(t3_num) + + t3_middle_v_vars = model.add_m_variables(len(t3_nums), name="t3_v") + t3_middle_ph_vars = model.add_m_variables(len(t3_nums), name="t3_ph") + t3_closed_branch_p1_vars = model.add_m_variables(len(t3_closed_branch_leg_nums), name="t3_closed_branch_p1") + t3_closed_branch_p2_vars = model.add_m_variables(len(t3_closed_branch_leg_nums), name="t3_closed_branch_p2") + t3_closed_branch_q1_vars = model.add_m_variables(len(t3_closed_branch_leg_nums), name="t3_closed_branch_q1") + t3_closed_branch_q2_vars = model.add_m_variables(len(t3_closed_branch_leg_nums), name="t3_closed_branch_q2") + t3_open_side1_p2_vars = model.add_m_variables(len(t3_open_side2_leg_nums), name="t3_open_side1_branch_p2") + t3_open_side1_q2_vars = model.add_m_variables(len(t3_open_side2_leg_nums), name="t3_open_side1_branch_q2") + + return VariableContext(v_vars, ph_vars, + gen_p_vars, gen_q_vars, + bat_p_vars, bat_q_vars, + shunt_p_vars, shunt_q_vars, + svc_q_vars, + vsc_cs_p_vars, vsc_cs_q_vars, + closed_branch_p1_vars, closed_branch_q1_vars, + closed_branch_p2_vars, closed_branch_q2_vars, + open_side1_branch_p2_vars, open_side1_branch_q2_vars, + open_side2_branch_p1_vars, open_side2_branch_q1_vars, + dl_v_vars, dl_ph_vars, + dl_branch_p1_vars, dl_branch_p2_vars, + dl_branch_q1_vars, dl_branch_q2_vars, + t3_middle_v_vars, t3_middle_ph_vars, + t3_closed_branch_p1_vars, t3_closed_branch_p2_vars, + t3_closed_branch_q1_vars, t3_closed_branch_q2_vars, + t3_open_side1_p2_vars, t3_open_side1_q2_vars, + branch_num_2_index, + gen_p_num_2_index, gen_q_num_2_index, + bat_p_num_2_index, bat_q_num_2_index, + shunt_num_2_index, + svc_num_2_index, + vsc_cs_num_2_index, + dl_num_2_index, + t3_num_2_index, t3_leg1_num_2_index, t3_leg2_num_2_index, t3_leg3_num_2_index) + + def _update_generators(self, network_cache: NetworkCache, model: ipopt.Model): + connected_gen_ids = [] + connected_gen_target_p = [] + connected_gen_target_q = [] + connected_gen_target_v = [] + connected_gen_voltage_regulator_on = [] + connected_gen_p = [] + connected_gen_q = [] + disconnected_gen_ids = [] + for gen_num, (gen_id, row) in enumerate(network_cache.generators.iterrows()): + bus_id = row.bus_id + if bus_id: + connected_gen_ids.append(gen_id) + + gen_p_index = self.gen_p_num_2_index[gen_num] + gen_q_index = self.gen_q_num_2_index[gen_num] + + p = model.get_value(self.gen_p_vars[gen_p_index]) + target_p = -p + connected_gen_target_p.append(target_p) + connected_gen_p.append(p) + + if gen_q_index == -1: + target_q = row.target_q + q = -row.target_q + target_v = row.target_v + voltage_regulator_on = False + else: + q = model.get_value(self.gen_q_vars[gen_q_index]) + target_q = -q + + regulated_bus_num = network_cache.buses.index.get_loc(row.regulated_bus_id) + v = model.get_value(self.v_vars[regulated_bus_num]) + target_v = v + + q_bounds = Bounds.get_reactive_power_bounds(row).mirror() + voltage_regulator_on = q_bounds.contains(q) + + connected_gen_target_q.append(target_q) + connected_gen_target_v.append(target_v) + connected_gen_voltage_regulator_on.append(voltage_regulator_on) + connected_gen_q.append(q) + + logger.log(TRACE_LEVEL, f"Update generator '{gen_id}' (num={gen_num}): target_p={target_p}, target_q={target_q}, target_v={target_v}, voltage_regulator_on={voltage_regulator_on}, p={p}, q={q}") + else: + disconnected_gen_ids.append(gen_id) + logger.log(TRACE_LEVEL, f"Update disconnected generator '{gen_id}' (num={gen_num})") + + network_cache.update_generators(connected_gen_ids, + connected_gen_target_p, + connected_gen_target_q, + connected_gen_target_v, + connected_gen_voltage_regulator_on, + connected_gen_p, + connected_gen_q, + disconnected_gen_ids) + + def _update_batteries(self, network_cache: NetworkCache, model: ipopt.Model): + connected_bat_ids = [] + connected_bat_target_p = [] + connected_bat_target_q = [] + connected_bat_target_v = [] + connected_bat_voltage_regulator_on = [] + connected_bat_p = [] + connected_bat_q = [] + disconnected_bat_ids = [] + for bat_num, (bat_id, row) in enumerate(network_cache.batteries.iterrows()): + bus_id = row.bus_id + if bus_id: + connected_bat_ids.append(bat_id) + + bat_p_index = self.bat_p_num_2_index[bat_num] + bat_q_index = self.bat_q_num_2_index[bat_num] + + p = model.get_value(self.bat_p_vars[bat_p_index]) + target_p = -p + connected_bat_target_p.append(target_p) + connected_bat_p.append(p) + + if bat_q_index == -1: + target_q = row.target_q + q = -row.target_q + target_v = row.target_v + voltage_regulator_on = False + else: + q = model.get_value(self.bat_q_vars[bat_q_index]) + target_q = -q + + bus_num = network_cache.buses.index.get_loc(row.bus_id) + v = model.get_value(self.v_vars[bus_num]) * row.nominal_v # FIXME + target_v = v + + q_bounds = Bounds.get_reactive_power_bounds(row).mirror() + voltage_regulator_on = q_bounds.contains(q) + + connected_bat_target_q.append(target_q) + connected_bat_target_v.append(target_v) + connected_bat_voltage_regulator_on.append(voltage_regulator_on) + connected_bat_q.append(q) + + logger.log(TRACE_LEVEL, f"Update battery '{bat_id}' (num={bat_num}): target_p={target_p}, target_q={target_q}, target_v={target_v}, voltage_regulator_on={voltage_regulator_on}, p={p}, q={q}") + else: + disconnected_bat_ids.append(bat_id) + logger.log(TRACE_LEVEL, f"Update disconnected battery '{bat_id}' (num={bat_num})") + + network_cache.update_batteries(connected_bat_ids, + connected_bat_target_p, + connected_bat_target_q, + connected_bat_target_v, + connected_bat_voltage_regulator_on, + connected_bat_p, + connected_bat_q, + disconnected_bat_ids) + + def _update_vsc_converter_stations(self, network_cache: NetworkCache, model: ipopt.Model): + connected_vsc_cs_ids = [] + connected_vsc_cs_target_q = [] + connected_vsc_cs_target_v = [] + connected_vsc_cs_voltage_regulator_on = [] + connected_vsc_cs_p = [] + connected_vsc_cs_q = [] + disconnected_vsc_cs_ids = [] + for vsc_cs_num, (vsc_cs_id, row) in enumerate(network_cache.vsc_converter_stations.iterrows()): + bus_id = row.bus_id + if bus_id: + connected_vsc_cs_ids.append(vsc_cs_id) + vsc_cs_index = self.vsc_cs_num_2_index[vsc_cs_num] + + p = model.get_value(self.vsc_cs_p_vars[vsc_cs_index]) + connected_vsc_cs_p.append(p) + + q = model.get_value(self.vsc_cs_q_vars[vsc_cs_index]) + target_q = -q + connected_vsc_cs_target_q.append(target_q) + connected_vsc_cs_q.append(q) + + regulated_bus_num = network_cache.buses.index.get_loc(row.regulated_bus_id) + v = model.get_value(self.v_vars[regulated_bus_num]) + target_v = v + connected_vsc_cs_target_v.append(target_v) + + q_bounds = Bounds.get_reactive_power_bounds(row).mirror() + voltage_regulator_on = q_bounds.contains(q) + connected_vsc_cs_voltage_regulator_on.append(voltage_regulator_on) + + logger.log(TRACE_LEVEL, f"Update VSC converter station '{vsc_cs_id}' (num={vsc_cs_num}): target_q={target_q}, target_v={target_v}, voltage_regulator_on={voltage_regulator_on}, p={p}, q={q}") + else: + disconnected_vsc_cs_ids.append(vsc_cs_id) + logger.log(TRACE_LEVEL, f"Update disconnected VSC converter station '{vsc_cs_id}' (num={vsc_cs_num})") + + network_cache.update_vsc_converter_stations(connected_vsc_cs_ids, + connected_vsc_cs_target_q, + connected_vsc_cs_target_v, + connected_vsc_cs_voltage_regulator_on, + connected_vsc_cs_p, + connected_vsc_cs_q, + disconnected_vsc_cs_ids) + + def _update_hvdc_lines(self, network_cache: NetworkCache, model: ipopt.Model): + hvdc_line_ids = [] + hvdc_line_target_p = [] + for hvdc_line_num, (hvdc_line_id, hvdc_line_row) in enumerate(network_cache.hvdc_lines.iterrows()): + hvdc_line_ids.append(hvdc_line_id) + vsc_cs1_num = network_cache.vsc_converter_stations.index.get_loc(hvdc_line_row.converter_station1_id) + vsc_cs2_num = network_cache.vsc_converter_stations.index.get_loc(hvdc_line_row.converter_station2_id) + p1 = model.get_value(self.vsc_cs_p_vars[vsc_cs1_num]) + p2 = model.get_value(self.vsc_cs_p_vars[vsc_cs2_num]) + target_p = abs(p1) if NetworkCache.is_rectifier(hvdc_line_row.converter_station1_id, hvdc_line_row) else abs(p2) + hvdc_line_target_p.append(target_p) + + logger.log(TRACE_LEVEL, f"Update HVDC line '{hvdc_line_id}': target_p={target_p}") + + network_cache.update_hvdc_lines(hvdc_line_ids, hvdc_line_target_p) + + def _update_static_var_compensators(self, network_cache: NetworkCache, model: ipopt.Model): + connected_svc_ids = [] + connected_svc_target_q = [] + connected_svc_target_v = [] + connected_svc_regulation_mode = [] + connected_svc_p = [] + connected_svc_q = [] + disconnected_svc_ids = [] + for svc_num, (svc_id, row) in enumerate(network_cache.static_var_compensators.iterrows()): + bus_id = row.bus_id + if bus_id: + connected_svc_ids.append(svc_id) + + connected_svc_p.append(0.0) + + svc_index = self.svc_num_2_index[svc_num] + q = model.get_value(self.svc_q_vars[svc_index]) + target_q = -q + connected_svc_target_q.append(target_q) + connected_svc_q.append(q) + + regulated_bus_num = network_cache.buses.index.get_loc(row.regulated_bus_id) + v = model.get_value(self.v_vars[regulated_bus_num]) + target_v = v + connected_svc_target_v.append(target_v) + + q_bounds = Bounds(row.b_min * v * v, row.b_max * v * v) + regulation_mode = 'VOLTAGE' if q_bounds.contains(q) else 'REACTIVE_POWER' + connected_svc_regulation_mode.append(regulation_mode) + + logger.log(TRACE_LEVEL, f"Update SVC '{svc_id}' (num={svc_num}): target_q={target_q}, target_v={target_v}, regulation_mode={regulation_mode}") + else: + disconnected_svc_ids.append(svc_id) + logger.log(TRACE_LEVEL, f"Update disconnected SVC '{svc_id}' (num={svc_num})") + + network_cache.update_static_var_compensators(connected_svc_ids, + connected_svc_target_q, + connected_svc_target_v, + connected_svc_regulation_mode, + connected_svc_p, + connected_svc_q, + disconnected_svc_ids) + + def _update_shunt_compensators(self, network_cache: NetworkCache, model: ipopt.Model): + connected_shunt_ids = [] + connected_shunt_p = [] + connected_shunt_q = [] + disconnected_shunt_ids = [] + for shunt_num, (shunt_id, row) in enumerate(network_cache.shunts.iterrows()): + bus_id = row.bus_id + if bus_id: + shunt_index = self.shunt_num_2_index[shunt_num] + p = model.get_value(self.shunt_p_vars[shunt_index]) + q = model.get_value(self.shunt_q_vars[shunt_index]) + connected_shunt_ids.append(shunt_id) + connected_shunt_p.append(p) + connected_shunt_q.append(q) + + logger.log(TRACE_LEVEL, f"Update shunt '{shunt_id}' (num={shunt_num}): p={p} q={q}") + else: + disconnected_shunt_ids.append(shunt_id) + logger.log(TRACE_LEVEL, f"Update disconnected shunt '{shunt_id}' (num={shunt_num})") + + network_cache.update_shunt_compensators(connected_shunt_ids, connected_shunt_p, connected_shunt_q, disconnected_shunt_ids) + + def _update_branches(self, network_cache: NetworkCache, model: ipopt.Model): + branch_ids = [] + branch_p1 = [] + branch_p2 = [] + branch_q1 = [] + branch_q2 = [] + for branch_num, (branch_id, row) in enumerate(network_cache.branches.iterrows()): + branch_index = self.branch_num_2_index[branch_num] + branch_ids.append(branch_id) + if row.bus1_id and row.bus2_id: + p1 = model.get_value(self.closed_branch_p1_vars[branch_index]) + p2 = model.get_value(self.closed_branch_p2_vars[branch_index]) + q1 = model.get_value(self.closed_branch_q1_vars[branch_index]) + q2 = model.get_value(self.closed_branch_q2_vars[branch_index]) + elif row.bus2_id: + p1 = 0.0 + p2 = model.get_value(self.open_side1_branch_p2_vars[branch_index]) + q1 = 0.0 + q2 = model.get_value(self.open_side1_branch_q2_vars[branch_index]) + elif row.bus1_id: + p1 = model.get_value(self.open_side2_branch_p1_vars[branch_index]) + p2 = 0.0 + q1 = model.get_value(self.open_side2_branch_q1_vars[branch_index]) + q2 = 0.0 + else: + p1 = 0.0 + p2 = 0.0 + q1 = 0.0 + q2 = 0.0 + + branch_p1.append(p1) + branch_p2.append(p2) + branch_q1.append(q1) + branch_q2.append(q2) + + logger.log(TRACE_LEVEL, f"Update branch '{branch_id}': p1={p1} p2={p2} q1={q1} q2={q2}") + + network_cache.update_branches(branch_ids, branch_p1, branch_p2, branch_q1, branch_q2) + + def _update_transformers_3w(self, network_cache: NetworkCache, model: ipopt.Model): + t3_ids = [] + t3_p1 = [] + t3_p2 = [] + t3_p3 = [] + t3_q1 = [] + t3_q2 = [] + t3_q3 = [] + t3_v = [] + t3_angle = [] + for t3_num, (t3_id, t3_row) in enumerate(network_cache.transformers_3w.iterrows()): + t3_ids.append(t3_id) + t3_index = self.t3_num_2_index[t3_num] + if t3_row.bus1_id or t3_row.bus2_id or t3_row.bus3_id: + leg1_index = self.t3_leg1_num_2_index[t3_num] + leg2_index = self.t3_leg2_num_2_index[t3_num] + leg3_index = self.t3_leg3_num_2_index[t3_num] + + if t3_row.bus1_id: + p1 = model.get_value(self.t3_closed_branch_p1_vars[leg1_index]) + q1 = model.get_value(self.t3_closed_branch_q1_vars[leg1_index]) + else: + p1 = 0 + q1 = 0 + + if t3_row.bus2_id: + p2 = model.get_value(self.t3_closed_branch_p1_vars[leg2_index]) + q2 = model.get_value(self.t3_closed_branch_q1_vars[leg2_index]) + else: + p2 = 0 + q2 = 0 + + if t3_row.bus3_id: + p3 = model.get_value(self.t3_closed_branch_p1_vars[leg3_index]) + q3 = model.get_value(self.t3_closed_branch_q1_vars[leg3_index]) + else: + q3 = 0 + p3 = 0 + + v = model.get_value(self.t3_middle_v_vars[t3_index]) + angle = model.get_value(self.t3_middle_ph_vars[t3_index]) + else: + p1 = 0.0 + p2 = 0.0 + p3 = 0.0 + q1 = 0.0 + q2 = 0.0 + q3 = 0.0 + v = math.nan + angle = math.nan + + t3_p1.append(p1) + t3_p2.append(p2) + t3_p3.append(p3) + t3_q1.append(q1) + t3_q2.append(q2) + t3_q3.append(q3) + t3_v.append(v * t3_row.rated_u0) + t3_angle.append(math.degrees(angle)) + + logger.log(TRACE_LEVEL, f"Update 3 windings transformer '{t3_id}': p1={p1} p2={p2} p3={p3} q1={q1} q2={q2} q3={q3} v={v} angle={angle}") + + network_cache.update_transformers_3w(t3_ids, t3_p1, t3_p2, t3_p3, t3_q1, t3_q2, t3_q3, t3_v, t3_angle) + + def _update_buses(self, network_cache: NetworkCache, model: ipopt.Model): + bus_ids = [] + bus_v_mag = [] + bus_v_angle = [] + for bus_num, (bus_id, row) in enumerate(network_cache.buses.iterrows()): + bus_ids.append(bus_id) + v = model.get_value(self.v_vars[bus_num]) + bus_v_mag.append(v) + angle = model.get_value(self.ph_vars[bus_num]) + bus_v_angle.append(angle) + + logger.log(TRACE_LEVEL, f"Update bus '{bus_id}' (num={bus_num}): v={v}, angle={angle}") + + network_cache.update_buses(bus_ids, bus_v_mag, bus_v_angle) + + def _update_dangling_lines(self, network_cache: NetworkCache, model: ipopt.Model): + connected_dl_ids = [] + connected_dl_v = [] + connected_dl_angle = [] + connected_dl_p = [] + connected_dl_q = [] + disconnected_dl_ids = [] + for dl_num, (dl_id, row) in enumerate(network_cache.dangling_lines.iterrows()): + dl_index = self.dl_num_2_index[dl_num] + + if row.bus_id: + v = model.get_value(self.dl_v_vars[dl_index]) + angle = model.get_value(self.dl_ph_vars[dl_index]) + connected_dl_ids.append(dl_id) + connected_dl_v.append(v * row.nominal_v) + connected_dl_angle.append(math.degrees(angle)) + p = model.get_value(self.dl_branch_p1_vars[dl_index]) + q = model.get_value(self.dl_branch_q1_vars[dl_index]) + connected_dl_p.append(p) + connected_dl_q.append(q) + logger.log(TRACE_LEVEL, f"Update dangling line '{dl_id}' (num={dl_num}): v={v}, angle={angle}, p={p}, q={q}") + else: + disconnected_dl_ids.append(dl_id) + logger.log(TRACE_LEVEL, f"Update disconnected dangling line '{dl_id}' (num={dl_num})") + + network_cache.update_dangling_lines(connected_dl_ids, + connected_dl_v, + connected_dl_angle, + connected_dl_p, + connected_dl_q, + disconnected_dl_ids) + + def update_network(self, network_cache: NetworkCache, model: ipopt.Model) -> None: + self._update_generators(network_cache, model) + self._update_batteries(network_cache, model) + self._update_vsc_converter_stations(network_cache, model) + self._update_hvdc_lines(network_cache, model) + self._update_static_var_compensators(network_cache, model) + self._update_shunt_compensators(network_cache, model) + self._update_branches(network_cache, model) + self._update_transformers_3w(network_cache, model) + self._update_buses(network_cache, model) + self._update_dangling_lines(network_cache, model) diff --git a/pypowsybl/opf/impl/network_statistics.py b/pypowsybl/opf/impl/network_statistics.py new file mode 100644 index 0000000000..f6a2426532 --- /dev/null +++ b/pypowsybl/opf/impl/network_statistics.py @@ -0,0 +1,48 @@ +import logging +from typing import Optional + +import numpy as np +from pandas import Series +from tabulate import tabulate + +from pypowsybl._pypowsybl import ElementType +from pypowsybl.opf.impl.model.network_cache import NetworkCache + +logger = logging.getLogger(__name__) + + +class NetworkStatistics: + def __init__(self, network_cache: NetworkCache) -> None: + self._network_cache = network_cache + self._initial_values = {} + + def _get_column(self, element_type: ElementType, attribute_id: str) -> Optional[Series]: + if element_type == ElementType.GENERATOR: + return self._network_cache.generators[attribute_id] if len(self._network_cache.generators) > 0 else None + elif element_type == ElementType.BATTERY: + return self._network_cache.batteries[attribute_id] if len(self._network_cache.batteries) > 0 else None + elif element_type == ElementType.VSC_CONVERTER_STATION: + return self._network_cache.vsc_converter_stations[attribute_id] if len(self._network_cache.vsc_converter_stations) > 0 else None + else: + raise ValueError(f"Unknown element type: {element_type}") + + def add(self, element_type: ElementType, attribute_id: str): + column = self._get_column(element_type, attribute_id) + if column is not None: + self._initial_values[(element_type, attribute_id)] = column.copy() + + def print(self): + for (element_type, attribute_id), initial_values in self._initial_values.items(): + final_values = self._get_column(element_type, attribute_id) + + table_data = [ + ['Initial', initial_values.min(), initial_values.max(), initial_values.mean(), + initial_values.std(), np.median(initial_values)], + ['Final', final_values.min(), final_values.max(), final_values.mean(), + final_values.std(), np.median(final_values)] + ] + + headers = ['', 'Min', 'Max', 'Mean', 'Std', 'Median'] + table = tabulate(table_data, headers=headers, floatfmt='.3f', tablefmt='simple') + + logger.info(f"\nStatistics: {element_type} - {attribute_id}\n{table}") diff --git a/pypowsybl/opf/impl/opf.py b/pypowsybl/opf/impl/opf.py new file mode 100644 index 0000000000..9c75dfb277 --- /dev/null +++ b/pypowsybl/opf/impl/opf.py @@ -0,0 +1,117 @@ +import logging +import time + +import pyoptinterface as poi +from pypowsybl._pypowsybl import ElementType + +from pypowsybl.network import Network +from pypowsybl.opf.impl.bounds.battery_power_bounds import BatteryPowerBounds +from pypowsybl.opf.impl.bounds.bus_voltage_bounds import BusVoltageBounds +from pypowsybl.opf.impl.bounds.dangling_line_voltage_bounds import DanglingLineVoltageBounds +from pypowsybl.opf.impl.bounds.generator_power_bounds import GeneratorPowerBounds +from pypowsybl.opf.impl.bounds.slack_bus_angle_bounds import SlackBusAngleBounds +from pypowsybl.opf.impl.bounds.transformer_3w_middle_voltage_bounds import Transformer3wMiddleVoltageBounds +from pypowsybl.opf.impl.bounds.vsc_cs_power_bounds import VscCsPowerBounds +from pypowsybl.opf.impl.constraints.branch_flow_constraints import BranchFlowConstraints +from pypowsybl.opf.impl.constraints.current_limit_constraints import CurrentLimitConstraints +from pypowsybl.opf.impl.constraints.dangling_line_flow_constraints import DanglingLineFlowConstraints +from pypowsybl.opf.impl.constraints.hvdc_line_constraints import HvdcLineConstraints +from pypowsybl.opf.impl.constraints.power_balance_constraints import PowerBalanceConstraints +from pypowsybl.opf.impl.constraints.shunt_flow_constraints import ShuntFlowConstraints +from pypowsybl.opf.impl.constraints.static_var_compensator_reactive_limits_constraints import \ + StaticVarCompensatorReactiveLimitsConstraints +from pypowsybl.opf.impl.constraints.transformer_3w_flow_constraints import Transformer3wFlowConstraints +from pypowsybl.opf.impl.costs.minimize_against_reference_cost_function import MinimizeAgainstReferenceCostFunction +from pypowsybl.opf.impl.costs.redispatching_cost_function import RedispatchingCostFunction +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.model.constraints import Constraints +from pypowsybl.opf.impl.model.model_parameters import ModelParameters +from pypowsybl.opf.impl.model.network_cache import NetworkCache +from pypowsybl.opf.impl.model.opf_model import OpfModel +from pypowsybl.opf.impl.network_statistics import NetworkStatistics +from pypowsybl.opf.impl.parameters import OptimalPowerFlowParameters, OptimalPowerFlowMode + +logger = logging.getLogger(__name__) + + +# pip install pyoptinterface llvmlite tccbox +# +# git clone https://github.com/coin-or-tools/ThirdParty-Mumps.git +# cd ThirdParty-Mumps +# ./get.Mumps +# ./configure --prefix $HOME/mumps +# make -j 8 +# make install +# +# git clone https://github.com/coin-or/Ipopt +# cd Ipopt/ +# ./configure --prefix $HOME/ipopt --with-mumps-cflags="-I$HOME/mumps/include/coin-or/mumps/" --with-mumps-lflags="-L$HOME/mumps/lib -lcoinmumps" +# make -j 8 +# make install +# +# export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/mumps/lib:$HOME/ipopt/lib +# +class OptimalPowerFlow: + def __init__(self, network: Network) -> None: + self._network = network + + def run(self, parameters: OptimalPowerFlowParameters) -> bool: + network_cache = NetworkCache(self._network) + + variable_bounds = [BusVoltageBounds(), + SlackBusAngleBounds(), + GeneratorPowerBounds(), + BatteryPowerBounds(), + VscCsPowerBounds(), + DanglingLineVoltageBounds(), + Transformer3wMiddleVoltageBounds()] + constraints: list[Constraints] = [BranchFlowConstraints(), + ShuntFlowConstraints(), + StaticVarCompensatorReactiveLimitsConstraints(), + HvdcLineConstraints(), + PowerBalanceConstraints(), + DanglingLineFlowConstraints(), + Transformer3wFlowConstraints()] + if parameters.mode == OptimalPowerFlowMode.REDISPATCHING: + constraints.append(CurrentLimitConstraints()) + cost_function = RedispatchingCostFunction(1.0, 1.0, 1.0) + else: + cost_function = MinimizeAgainstReferenceCostFunction() + model_parameters = ModelParameters(parameters.reactive_bounds_reduction, + parameters.twt_split_shunt_admittance, + Bounds(parameters.default_voltage_bounds[0], parameters.default_voltage_bounds[1])) + opf_model = OpfModel.build(network_cache, model_parameters, variable_bounds, constraints, cost_function) + + network_stats = NetworkStatistics(network_cache) + network_stats.add(ElementType.GENERATOR, 'target_v') + network_stats.add(ElementType.BATTERY, 'target_v') + network_stats.add(ElementType.VSC_CONVERTER_STATION, 'target_v') + network_stats.add(ElementType.GENERATOR, 'target_p') + network_stats.add(ElementType.BATTERY, 'target_p') + network_stats.add(ElementType.VSC_CONVERTER_STATION, 'target_p') + + logger.info("Starting optimization...") + start = time.perf_counter() + + opf_model.model.set_model_attribute(poi.ModelAttribute.Silent, False) + opf_model.model.optimize() + status = opf_model.model.get_model_attribute(poi.ModelAttribute.TerminationStatus) + + logger.info(f"Optimization ends with status {status} in {time.perf_counter() - start:.3f} seconds.") + + # for debugging + opf_model.analyze_violations(model_parameters) + + # update network + opf_model.update_network() + + network_stats.print() + + network_cache.network.per_unit = False # FIXME design to improve + + return status == poi.TerminationStatusCode.LOCALLY_SOLVED + + +def run_ac(network: Network, parameters: OptimalPowerFlowParameters = OptimalPowerFlowParameters()) -> bool: + opf = OptimalPowerFlow(network) + return opf.run(parameters) diff --git a/pypowsybl/opf/impl/parameters.py b/pypowsybl/opf/impl/parameters.py new file mode 100644 index 0000000000..0219b574d2 --- /dev/null +++ b/pypowsybl/opf/impl/parameters.py @@ -0,0 +1,34 @@ +from enum import Enum + + +class OptimalPowerFlowMode(Enum): + LOADFLOW = "LOADFLOW" + REDISPATCHING = "REDISPATCHING" + + +class OptimalPowerFlowParameters: + def __init__(self, + reactive_bounds_reduction: float = 0.1, + twt_split_shunt_admittance = False, + mode: OptimalPowerFlowMode = OptimalPowerFlowMode.LOADFLOW, + default_voltage_bounds: tuple[float, float] = (0.8, 1.1)) -> None: + self._reactive_bounds_reduction = reactive_bounds_reduction + self._twt_split_shunt_admittance = twt_split_shunt_admittance + self._mode = mode + self._default_voltage_bounds = default_voltage_bounds + + @property + def reactive_bounds_reduction(self) -> float: + return self._reactive_bounds_reduction + + @property + def twt_split_shunt_admittance(self) -> bool: + return self._twt_split_shunt_admittance + + @property + def mode(self) -> OptimalPowerFlowMode: + return self._mode + + @property + def default_voltage_bounds(self) -> tuple[float, float]: + return self._default_voltage_bounds diff --git a/pypowsybl/opf/impl/util.py b/pypowsybl/opf/impl/util.py new file mode 100644 index 0000000000..08fa274bc1 --- /dev/null +++ b/pypowsybl/opf/impl/util.py @@ -0,0 +1,4 @@ +import logging + +TRACE_LEVEL = 5 +logging.addLevelName(TRACE_LEVEL, "TRACE") diff --git a/pytest.ini b/pytest.ini index 31ddc98a40..ed2ff3b990 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,4 +1,4 @@ [pytest] #Change those options to get logs -log_cli=false -log_level=ERROR +log_cli=true +log_level=DEBUG diff --git a/requirements.txt b/requirements.txt index 9e208d9b1d..cce4ee1ea7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,6 +5,8 @@ networkx matplotlib==3.9.2 pybind11[global]==2.13.6 cmake>=3.20 +pyoptinterface==0.5.1 +tabulate==0.9.0 # optional dependencies pandapower==3.1.2; python_version<'3.13' # TO REMOVE WHEN COMPATIBLE VERSION AVAILABLE diff --git a/tests/test_network_extensions.py b/tests/test_network_extensions.py index f163238d74..1be8f0ce2d 100644 --- a/tests/test_network_extensions.py +++ b/tests/test_network_extensions.py @@ -589,6 +589,30 @@ def test_batteries_voltage_regulation(): assert network.get_extensions('voltageRegulation').empty +def test_batteries_voltage_regulation(): + network = pn.load(str(TEST_DIR.joinpath('battery.xiidm'))) + assert network.get_extensions('voltageRegulation').empty + + network.create_extensions('voltageRegulation', id='BAT', voltage_regulator_on=True, target_v=400.0) + e = network.get_extensions('voltageRegulation') + expected = pd.DataFrame( + index=pd.Series(name='id', data=['BAT']), + columns=['voltage_regulator_on', 'target_v'], + data=[[True, 400.0]]) + pd.testing.assert_frame_equal(expected, e, check_dtype=False) + + network.update_extensions('voltageRegulation', id=['BAT'], voltage_regulator_on=False, target_v=399.0) + e = network.get_extensions('voltageRegulation') + expected = pd.DataFrame( + index=pd.Series(name='id', data=['BAT']), + columns=['voltage_regulator_on', 'target_v'], + data=[[False, 399.0]]) + pd.testing.assert_frame_equal(expected, e, check_dtype=False) + + network.remove_extensions('voltageRegulation', ['BAT']) + assert network.get_extensions('voltageRegulation').empty + + def test_get_extensions_information(): extensions_information = pypowsybl.network.get_extensions_information() assert extensions_information.loc['cgmesMetadataModels']['detail'] == 'Provides information about CGMES metadata models' diff --git a/tests/test_opf.py b/tests/test_opf.py new file mode 100644 index 0000000000..958c565440 --- /dev/null +++ b/tests/test_opf.py @@ -0,0 +1,204 @@ +# +# Copyright (c) 2025, RTE (http://www.rte-france.com) +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +import logging + +import numpy as np +import pandas as pd +import pytest +from pandas import DataFrame + +import pypowsybl as pp +from pypowsybl.network import Network +from pypowsybl.opf.impl.model.bounds import Bounds +from pypowsybl.opf.impl.opf import OptimalPowerFlowParameters +from pypowsybl.opf.impl.parameters import OptimalPowerFlowMode + + +@pytest.fixture(autouse=True) +def set_up(): + pp.set_config_read(False) + logging.basicConfig(level=logging.DEBUG) + logging.getLogger('powsybl').setLevel(1) + pd.options.display.max_columns = None + pd.options.display.expand_frame_repr = False + + +def test_reactive_range(): + b = Bounds(10, 20) + assert b.min_value == 10 + assert b.max_value == 20 + assert str(b) == "[10, 20]" + assert b.mirror().min_value == -20 + assert b.mirror().max_value == -10 + assert b.mirror().min_value_with_margin == -20.000001 + assert b.mirror().max_value_with_margin == -9.999999 + assert b.contains(10) + assert b.contains(20) + assert not b.contains(5) + assert not b.contains(21) + assert b.contains(9.99999999) + assert not b.contains(9.999) + assert b.contains(20.00000001) + assert not b.contains(20.001) + assert b.reduce(0.1).min_value == 11.0 + assert b.reduce(0.1).max_value == 18.0 + assert b.mirror().reduce(0.1).min_value == -18.0 + assert b.mirror().reduce(0.1).max_value == -11.0 + assert b.mirror().reduce(0.1).min_value_with_margin == -18.000001 + assert b.mirror().reduce(0.1).max_value_with_margin == -10.999999 + + +def create_loadflow_parameters(): + return pp.loadflow.Parameters(voltage_init_mode=pp.loadflow.VoltageInitMode.DC_VALUES, + hvdc_ac_emulation=False, + provider_parameters={'plausibleActivePowerLimit': '10000.0', + 'svcVoltageMonitoring': 'false'}) + + +def run_opf_then_lf(network: pp.network.Network, + opf_parameters: OptimalPowerFlowParameters = OptimalPowerFlowParameters(), + iteration_count: int = 1): + lf_parameters = create_loadflow_parameters() + lf_result = pp.loadflow.run_ac(network, lf_parameters) + assert lf_result[0].status == pp.loadflow.ComponentStatus.CONVERGED + + assert pp.opf.run_ac(network, opf_parameters) + + validate(network) + + lf_parameters.voltage_init_mode = pp.loadflow.VoltageInitMode.PREVIOUS_VALUES + lf_result = pp.loadflow.run_ac(network, lf_parameters) + assert lf_result[0].status == pp.loadflow.ComponentStatus.CONVERGED + assert lf_result[0].iteration_count == iteration_count + + +def validate(network: Network): + validation_parameters = pp.loadflow.ValidationParameters(threshold=1, check_main_component_only=True) + validation_types = [ + pp.loadflow.ValidationType.BUSES, + pp.loadflow.ValidationType.FLOWS, + # pp.loadflow.ValidationType.SHUNTS, to fix because active power is expected nan + # pp.loadflow.ValidationType.TWTS, to fix + pp.loadflow.ValidationType.TWTS3W, + # pp.loadflow.ValidationType.GENERATORS # to fix because remote voltage not taken into account + ] + result = pp.loadflow.run_validation(network, validation_types, validation_parameters) + # print(result.buses[result.buses['validated'] == False]) + assert result.valid + + +def calculate_overloading(network: Network) -> DataFrame: + sa = pp.security.create_analysis() + lf_parameters = create_loadflow_parameters() + sa_parameters = pp.security.Parameters(lf_parameters) + sa_results = sa.run_ac(network, sa_parameters) + print(sa_results.pre_contingency_result.limit_violations) + limit_violations = sa_results.limit_violations + limit_violations = limit_violations[ + (limit_violations['limit_type'] == 'CURRENT') & (limit_violations['limit_name'] == 'permanent')] + limit_violations['loading'] = limit_violations['value'] / limit_violations['limit'] + limit_violations = limit_violations[['side', 'loading']] + return limit_violations + + +def increase_load(network: Network, value: float): + loads = network.get_loads() + load_sum = loads['p0'].sum() + loads['p0'] = loads['p0'] * (1 + value / load_sum) + network.update_loads(id=loads.index, p0=loads['p0']) + + +def test_esg_tuto1(): + run_opf_then_lf(pp.network.create_eurostag_tutorial_example1_network()) + + +def test_ieee9(): + run_opf_then_lf(pp.network.create_ieee9()) + + +def test_ieee14(): + run_opf_then_lf(pp.network.create_ieee14()) + + +def test_ieee14_redispatching(): + network = pp.network.create_ieee14() + + # add a current limit on L7-9-1 + network.create_operational_limits(pd.DataFrame.from_records(index='element_id', data=[ + {'element_id': 'L7-9-1', 'name': 'permanent', 'side': 'ONE', + 'type': 'CURRENT', 'value': 1000, 'acceptable_duration': np.inf, + 'fictitious': False, 'group_name': 'GROUP1'}, + ])) + network.update_lines(id=['L7-9-1'], selected_limits_group_1=['GROUP1']) + + lf_parameters = create_loadflow_parameters() + lf_result = pp.loadflow.run_ac(network, lf_parameters) + assert lf_result[0].status == pp.loadflow.ComponentStatus.CONVERGED + + assert network.get_lines(attributes=['i1']).loc['L7-9-1'].i1 == pytest.approx(1113.514, rel=1e-3, abs=1e-3) + assert len(calculate_overloading(network)) == 1 + + opf_parameters = OptimalPowerFlowParameters(mode=OptimalPowerFlowMode.REDISPATCHING) + assert pp.opf.run_ac(network, opf_parameters) + + lf_parameters.voltage_init_mode = pp.loadflow.VoltageInitMode.PREVIOUS_VALUES + lf_result = pp.loadflow.run_ac(network, lf_parameters) + assert lf_result[0].status == pp.loadflow.ComponentStatus.CONVERGED + assert lf_result[0].iteration_count == 1 + + assert network.get_lines(attributes=['i1']).loc['L7-9-1'].i1 == pytest.approx(980.865, rel=1e-3, abs=1e-3) + assert len(calculate_overloading(network)) == 0 + + +def test_ieee14_open_side_1_branch(): + network = pp.network.create_ieee14() + network.update_lines(id=['L3-4-1'], connected1=[False]) + run_opf_then_lf(network) + row = network.get_branches(id=['L3-4-1'], attributes=['p1', 'q1', 'p2', 'q2' ]).head(1) + assert row.p1.values[0] == 0 + assert row.q1.values[0] == 0 + assert row.p2.values[0] != 0 + assert row.q2.values[0] != 0 + + +def test_ieee14_open_side_2_branch(): + network = pp.network.create_ieee14() + network.update_lines(id=['L3-4-1'], connected2=[False]) + run_opf_then_lf(network) + row = network.get_branches(id=['L3-4-1'], attributes=['p1', 'q1', 'p2', 'q2' ]).head(1) + assert row.p1.values[0] != 0 + assert row.q1.values[0] != 0 + assert row.p2.values[0] == 0 + assert row.q2.values[0] == 0 + + +def test_ieee30(): + run_opf_then_lf(pp.network.create_ieee30()) + + +def test_ieee57(): + run_opf_then_lf(pp.network.create_ieee57()) + + +def test_ieee118(): + run_opf_then_lf(pp.network.create_ieee118()) + + +def test_ieee300(): + run_opf_then_lf(pp.network.create_ieee300()) + + +def test_metrix_tutorial_six_buses(): + run_opf_then_lf(pp.network.create_metrix_tutorial_six_buses_network()) + + +def test_micro_grid_be(): + run_opf_then_lf(pp.network.create_micro_grid_be_network()) + + +def test_micro_grid_nl(): + run_opf_then_lf(pp.network.create_micro_grid_nl_network())