Adding calculation for better bigM value.

Author: rcschrg

src/monee/model/formulation/misoc/el.py

@@ -27,16 +27,38 @@ def equations(
         ]
 
 
+def _ell_physics_max(branch, w_max: float) -> float:
+    """Upper bound on per-unit squared current derived from voltage bounds alone.
+
+    From |I_ij| = |V_i - V_j| / |Z_ij| and |V_i|, |V_j| <= sqrt(W_max):
+        ell_ij <= (2*sqrt(W_max))^2 / |Z_ij|^2 = 4*W_max / (r^2 + x^2).
+    """
+    return 4 * w_max / (branch.br_r**2 + branch.br_x**2)
+
+
+def _big_m(w_max: float) -> float:
+    """Compute a tight big-M bound from the voltage bound alone.
+
+    Substituting the physics-based current bound ell_max = 4*W_max/|Z|^2 into
+    the Cauchy-Schwarz result M = (sqrt(W_max) + |Z|*sqrt(ell_max))^2, the
+    impedance cancels and M = 9*W_max, independent of branch impedance.
+    """
+    return 9 * w_max
+
+
 class MISOCPElectricityBranchFormulation(BranchFormulation):
     def ensure_var(self, branch):
-        branch.big_M = 20
         branch.current_pu = Var(1, min=0)
 
     def minimize(self, branch, grid, from_node_model, to_node_model, **kwargs):
         return [branch.current_pu * branch.br_r]
 
     def equations(self, branch, grid, from_node_model, to_node_model, **kwargs):
+        w_max = grid.vm_pu_max**2
+        big_m = _big_m(w_max)
+        ell_phys = _ell_physics_max(branch, w_max)
         return [
+            branch.current_pu <= ell_phys * branch.on_off,
             voltage_drop(
                 from_node_model.vars["vm_pu_squared"],
                 to_node_model.vars["vm_pu_squared"],
@@ -46,7 +68,7 @@ def equations(self, branch, grid, from_node_model, to_node_model, **kwargs):
                 branch.br_r,
                 branch.br_x,
             )
-            <= branch.big_M * (1 - branch.on_off),
+            <= big_m * (1 - branch.on_off),
             voltage_drop(
                 from_node_model.vars["vm_pu_squared"],
                 to_node_model.vars["vm_pu_squared"],
@@ -56,7 +78,7 @@ def equations(self, branch, grid, from_node_model, to_node_model, **kwargs):
                 branch.br_r,
                 branch.br_x,
             )
-            >= -branch.big_M * (1 - branch.on_off),
+            >= -big_m * (1 - branch.on_off),
             soc_rel(
                 from_node_model.vars["vm_pu_squared"],
                 branch.vars["p_from_mw"] / grid.sn_mva,
@@ -75,9 +97,4 @@ def equations(self, branch, grid, from_node_model, to_node_model, **kwargs):
                 branch.current_pu,
                 branch.br_x,
             ),
-            # Calculate the voltage/current error, However using this make the Problem non-convex
-            # branch.gap == gap_expr(from_node_model.vars["vm_pu_squared"],
-            #         branch.vars["p_from_mw"],
-            #         branch.vars["q_from_mvar"],
-            #         branch.current_pu),
         ]

src/monee/model/formulation/qp/gas.py

@@ -44,6 +44,7 @@ class PowerGrid(Grid):
     """
 
     sn_mva: float = 1
+    vm_pu_max: float = 1.5
 
 
 @model

src/monee/model/formulation/qp/water.py

@@ -0,0 +1,91 @@
+"""
+Test that inactive branches (on_off=0) in the MISOCP formulation
+carry exactly zero current and zero power flow.
+
+Topology: ring B0 → B1 → B2 → B0 where the B2→B0 tie-line is a backup
+branch with on_off=0.  All load must flow through B0→B1→B2; the
+backup line should be completely quiescent.
+"""
+
+import math
+
+import monee.model as mm
+from monee.model import Network
+from monee.model.branch import PowerLine
+from monee.model.child import ExtPowerGrid, PowerLoad
+from monee.model.formulation import MISOCP_NETWORK_FORMULATION
+from monee.model.grid import PowerGrid
+from monee.model.node import Bus
+from monee.solver import PyomoSolver
+
+
+def create_ring_with_backup() -> Network:
+    """Three-bus ring; the B2→B0 branch is the inactive backup tie-line."""
+    pn = Network(PowerGrid(name="power", sn_mva=1))
+
+    b0 = pn.node(
+        Bus(base_kv=1),
+        child_ids=[pn.child(ExtPowerGrid(p_mw=0, q_mvar=0, vm_pu=1.0, va_degree=0))],
+        grid=mm.EL,
+    )
+    b1 = pn.node(
+        Bus(base_kv=1),
+        child_ids=[pn.child(PowerLoad(p_mw=0.5, q_mvar=0))],
+        grid=mm.EL,
+    )
+    b2 = pn.node(
+        Bus(base_kv=1),
+        child_ids=[pn.child(PowerLoad(p_mw=0.5, q_mvar=0))],
+        grid=mm.EL,
+    )
+
+    line_kw = dict(r_ohm_per_m=0.00007, x_ohm_per_m=0.00007, parallel=1)
+    pn.branch(PowerLine(length_m=500, **line_kw), b0, b1)
+    pn.branch(PowerLine(length_m=500, **line_kw), b1, b2)
+    pn.branch(PowerLine(length_m=500, backup=True, on_off=0, **line_kw), b2, b0)
+
+    pn.apply_formulation(MISOCP_NETWORK_FORMULATION)
+    return pn
+
+
+def _backup_row(result):
+    """Return the result Series for the backup (on_off=0) line."""
+    lines = result.get(mm.PowerLine)
+    # The backup line has on_off=0; all others have on_off=1
+    return lines[lines["on_off"] == 0].iloc[0]
+
+
+def test_inactive_branch_has_zero_current():
+    result = PyomoSolver().solve(create_ring_with_backup())
+    backup = _backup_row(result)
+    assert math.isclose(backup["current_pu"], 0.0, abs_tol=1e-6), (
+        f"Inactive branch current_pu should be 0, got {backup['current_pu']}"
+    )
+
+
+def test_inactive_branch_has_zero_power_flow():
+    result = PyomoSolver().solve(create_ring_with_backup())
+    backup = _backup_row(result)
+    for col in ("p_from_mw", "p_to_mw", "q_from_mvar", "q_to_mvar"):
+        assert math.isclose(backup[col], 0.0, abs_tol=1e-4), (
+            f"Inactive branch {col} should be 0, got {backup[col]}"
+        )
+
+
+def test_active_lines_carry_full_load():
+    """Both loads must be fully served; only active lines carry current."""
+    result = PyomoSolver().solve(create_ring_with_backup())
+
+    lines = result.get(mm.PowerLine)
+    active = lines[lines["on_off"] == 1]
+
+    # Both active lines must carry non-zero current
+    assert (active["current_pu"] > 1e-4).all(), (
+        f"Active lines should carry current:\n{active['current_pu']}"
+    )
+
+    # All load is served (regulation = 1 for both loads)
+    loads = result.get(mm.PowerLoad)
+    assert (loads["regulation"] > 0.99).all(), (
+        f"All load should be served:\n{loads['regulation']}"
+    )