[WIP] AC OPF proto with pyoptinterface

cpp/powsybl-cpp/powsybl-cpp.cpp

@@ -587,7 +587,7 @@ std::shared_ptr<sensitivity_analysis_parameters> 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;

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"
 ]
 
 

pypowsybl/opf/__init__.py

@@ -0,0 +1 @@
+from .impl.opf import OptimalPowerFlow, run_ac

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))

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)

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)

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))

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})")

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)

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))

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)
+