monee.model.phys

monee.model.phys is the pure-math layer of monee: stateless, solver-agnostic functions that take plain variables or expressions and return the relational expressions assembled by the formulation classes in monee.model.formulation. The modules contain no classes, no solver imports and no model state; backend-specific operations are injected as *_impl callables (e.g. cos_impl, log_impl, sqrt_impl, pwl_impl), so the same physics runs on GEKKO, Pyomo or native monee variables.

Constants

Pipe hydraulics

Carrier-independent pipe hydraulics shared by the gas, water and heat models: Reynolds number, friction correlations, piecewise-linear breakpoint generation and sizing helpers.

monee.model.phys.core.hydraulics.calc_max_mass_flow(diameter_m, fluid_density_kg_per_m3, v_max_mps)

Per-pipe mass-flow upper bound [kg/s] from a velocity cap (\(\pi/4 \cdot D^2 \cdot \rho \cdot v_{max}\)). Used to tighten per-pipe big-M relaxations.

monee.model.phys.core.hydraulics.calc_min_diameter_for_mass_flow(mass_flow_kgs, fluid_density_kg_per_m3, v_design_mps)

Smallest pipe diameter [m] carrying mass_flow_kgs at v_design_mps.

Inverse of calc_max_mass_flow(): \(D = \sqrt{4 \cdot m / (\pi \cdot \rho \cdot v)}\). Used to size DHS pipes to their required capacity so the velocity cap does not bind. Returns 0 for degenerate inputs (\(m \le 0\) or non-positive \(\rho/v\)).

monee.model.phys.core.hydraulics.calc_pipe_area(diameter_m)

monee.model.phys.core.hydraulics.filter_near_linear(xs, ys, rtol=1e-06)

monee.model.phys.core.hydraulics.flow_rate_equation(mean_flow_velocity, flow_rate, diameter, fluid_density_kg_per_m3)

monee.model.phys.core.hydraulics.friction_at_high_re(diameter_m: float, roughness_m: float)

Turbulent asymptote \(f = 0.25 / \log_{10}(\varepsilon / (3.7 \cdot D))^2\) (Swamee-Jain at \(Re \to \infty\)). Within a few % of Colebrook for \(Re \gtrsim 10^4\); under-estimates pressure drop in laminar (\(Re < 2300\)). Returns 0 for degenerate inputs (\(D \le 0\), \(\varepsilon \ge 3.7 \cdot D\)).

Return type:

float

monee.model.phys.core.hydraulics.friction_value(Re, D, eps)

monee.model.phys.core.hydraulics.junction_mass_flow_balance(flows)

monee.model.phys.core.hydraulics.logspace(a, b, n)

monee.model.phys.core.hydraulics.phi_pwl_breakpoints(diameter_m: float, roughness_m: float, dynamic_visc_pas: float, pipe_area: float, m_max: float, n_breakpoints: int = 12)

Breakpoints for \(\varphi(m) = friction(Re(m)) \cdot m^2\) on \(m \in [0, m_{max}]\).

Log-spaced \([m_{max} \cdot 10^{-4}, m_{max}]\) plus a 0-anchor, so both the laminar tail (\(\varphi \propto m\)) and turbulent regime (\(\varphi \propto m^2\)) resolve well in a single SOS2.

monee.model.phys.core.hydraulics.piecewise_eq_friction(model, pwl)

monee.model.phys.core.hydraulics.pipe_mass_flow(max_v, min_v, v)

monee.model.phys.core.hydraulics.reynolds_equation(rey_var, mass_flow_kgs, diameter_m, dynamic_visc_pas, pipe_area)

monee.model.phys.core.hydraulics.swamee_jain(reynolds_var, diameter_m, roughness_m, log_func)

AC power flow

monee.model.phys.nonlinear.ac.calc_branch_t(tap, shift)

monee.model.phys.nonlinear.ac.int_flow_from_p(p_from_var, vm_from_pu, vm_to_pu, va_from_rad, va_to_rad, g_branch, b_branch, tap, shift, cos_impl=<built-in function cos>, sin_impl=<built-in function sin>, g_from=0, on_off=1)

monee.model.phys.nonlinear.ac.int_flow_from_q(q_from_var, vm_from_pu, vm_to_pu, va_from_rad, va_to_rad, g_branch, b_branch, tap, shift, cos_impl=<built-in function cos>, sin_impl=<built-in function sin>, b_from=0, on_off=1)

monee.model.phys.nonlinear.ac.int_flow_to_p(p_to_var, vm_from_pu, vm_to_pu, va_from_rad, va_to_rad, g_branch, b_branch, tap, shift, cos_impl=<built-in function cos>, sin_impl=<built-in function sin>, g_to_pu=0, on_off=1)

monee.model.phys.nonlinear.ac.int_flow_to_q(q_to_var, vm_from_pu, vm_to_pu, va_from_rad, va_to_rad, g_branch, b_branch, tap, shift, cos_impl=<built-in function cos>, sin_impl=<built-in function sin>, b_to_pu=0, on_off=1)

monee.model.phys.nonlinear.ac.int_flows(p_from_var, q_from_var, p_to_var, q_to_var, vm_from_pu, vm_to_pu, va_from_rad, va_to_rad, g_branch, b_branch, tap, shift, cos_impl=<built-in function cos>, sin_impl=<built-in function sin>, g_from=0, b_from=0, g_to_pu=0, b_to_pu=0, on_off=1)

All four AC branch flow equations (P/Q at both ends) in one pass.

Mathematically identical to calling int_flow_from_p(), int_flow_from_q(), int_flow_to_p() and int_flow_to_q() separately, but the sub-terms shared across the four (vm_from*vm_to, the sine/cosine of the bus-angle difference, vm_from**2, vm_to**2) are built once and reused. The to-direction reuses the from-direction’s cos/sin via the even/odd identities cos(-x) = cos(x) and sin(-x) = -sin(x). The numeric float coefficients are formed exactly as in the per-equation functions, so the only change is node sharing - the symbolic graph these equations contribute is ~halved, which speeds model assembly on every backend (the equation bodies are backend-agnostic).

Returns [p_from_eq, q_from_eq, p_to_eq, q_to_eq].

monee.model.phys.nonlinear.ac.power_balance_equation(signed_flows)

Gas flow (Weymouth)

monee.model.phys.nonlinear.gf.abs_pressure_pu(pressure_pu, pressure_ambient_pa, pressure_ref_pa)

Absolute pressure (per-unit of pressure_ref_pa) from a gauge pressure_pu.

monee.model.phys.nonlinear.gf.abs_psq_diff_pu(psq_pu_i, psq_pu_j, p_pu_i, p_pu_j, pressure_ambient_pa, pressure_ref_pa)

Absolute squared-pressure difference (per-unit of \(p_{ref}^2\)) when the node pressures are GAUGE.

The Weymouth law needs absolute pressures. With gauge node pressures \(p_g\) and ambient \(p_a\), the absolute is \(p = p_g + p_a\), so

\[p_i^2 - p_j^2 = (p_{g,i}+p_a)^2 - (p_{g,j}+p_a)^2 = (p_{sq,i}-p_{sq,j}) + 2\,a\,(p_{pu,i}-p_{pu,j})\]

with \(a = p_a/p_{ref}\) (the constant \(a^2\) cancels in the difference). With pressure_ambient_pa = 0 this is exactly the historical absolute form \(p_{sq,i} - p_{sq,j}\).

p_pu is whatever pressure the caller supplies: the NLP passes the exact \(\sqrt{p_{sq}}\) intermediate (it tolerates the non-linearity); the MILP and convex MIQCQP pass the LINEARIZED pressure \(p_0 + (p_{sq}-p_0^2)/(2 p_0)\) (affine in the decision variable, the same linearization they use for density), so the offset stays linear / conic.

monee.model.phys.nonlinear.gf.abs_psq_pu(psq_pu, p_pu, pressure_ambient_pa, pressure_ref_pa)

Single-node absolute squared pressure (per-unit), dropping the constant \(a^2\) that cancels in any i-j difference. Short-circuits to the gauge psq_pu when ambient is 0 (no p_pu reference). See abs_psq_diff_pu() for the exact-vs-linearized p_pu convention.

monee.model.phys.nonlinear.gf.calc_C_squared(diameter_m, length_m, t_k, compressibility, r_specific=504.5)

monee.model.phys.nonlinear.gf.pipe_weymouth(p_squared_i, p_squared_j, f_a_pos_sq, f_a_neg_sq, diameter_m, length_m, t_k, compressibility, on_off=1, friction=None, r_specific=504.5, **kwargs)

monee.model.phys.nonlinear.gf.reference_gas_density(grid)

Ideal-gas density [kg/m^3] at the grid’s reference ABSOLUTE pressure.

Absolute = pressure_ref_pa + pressure_ambient_pa (the latter is 0 in the default absolute convention). Used to size per-pipe velocity / mass-flow caps, so the cap reflects the true operating density under a gauge grid.

Water flow (Darcy-Weisbach)

monee.model.phys.nonlinear.wf.darcy_friction(reynolds_var)

monee.model.phys.nonlinear.wf.darcy_weisbach_equation(p_i, p_j, m_pos_sq, m_neg_sq, pipe_length, diameter_m, fluid_density_kg_per_m3, on_off=1, friction=None, **kwargs)

Heat flow

Smooth binary-free hydraulics

Smooth hydraulic primitives for pure-NLP (IPOPT/APOPT) solves: a single signed mass flow with |m| ≈ √(m² + ε²) replaces the direction binary and the positive/negative flow split.

Smooth, binary-free hydraulic primitives for pure-NLP (IPOPT/APOPT) solves.

The MISOCP-shaped models in gf/wf rely on a direction binary, a mass_flow_pos_kgs/mass_flow_neg_kgs split and an epigraph relaxation kept tight by an objective term - none of which survive an NLP relaxation. The helpers here express the same physics with a single signed mass flow and the smooth identity \(f_{pos}^2 - f_{neg}^2 = m \cdot |m|\) (with \(|m| \approx \sqrt{m^2+\varepsilon^2}\)), so the pressure drop is \(C^\infty\), monotonic and odd in the flow.

A pipe’s pressure drop is written as \(\Delta p = -\text{drop\_term}\) (in the carrier’s units), where drop_term is \(friction \cdot m \cdot |m|\) for the analytic friction models and a single odd cubic spline \(\psi(m)\) for the PWL model.

monee.model.phys.nonlinear.smooth.darcy_pressure(p_i, p_j, drop_term, length_m, diameter_m, fluid_density_kg_per_m3)

Signed Darcy-Weisbach: \((p_i-p_j) = -K \cdot \text{drop\_term}\) with \(K = L / (2 \cdot \rho \cdot A^2 \cdot D)\).

monee.model.phys.nonlinear.smooth.drop_term_and_eqs(friction_model, model, dynamic_visc_pas, area, signed, mag, m_max, **kwargs)

Return (drop_term, eqs) for the pressure-drop equation.

constant/nonlinear \(\to\) \(friction \cdot m \cdot |m|\) (plus the Reynolds/blend closure for nonlinear); hybrid \(\to\) a linear laminar term plus the quadratic turbulent term (QCQP, no closure); pwl \(\to\) an odd cubic spline \(\psi(\text{signed})\).

monee.model.phys.nonlinear.smooth.hybrid_drop_term(signed, mag, diameter_m, roughness_m, dynamic_visc_pas, area)

Laminar + fully-rough drop term: the QCQP-friendly “in-between” friction.

Mirrors pandapipes’ default nikuradse (\(f = 64/Re + \lambda_{rough}\)) but folds it into the drop instead of carrying an explicit friction factor. Since the laminar \(64/Re\) cancels one power of \(|m|\) (\((64/Re)\,m|m| = (64\mu A/D)\,m\)), the drop splits into a linear laminar term and the usual quadratic turbulent term:

\[\text{drop} = \underbrace{\tfrac{64\,\mu\,A}{D}}_{K_{lam}}\,m + \lambda_{rough}\,m\,|m|\]

The laminar part is exact Hagen-Poiseuille; the turbulent part is the fully-rough Darcy used by constant. No log/exp, no Reynolds/friction variable, no binaries, no high powers - smooth and representable with the quadratic \(m|m|\) monee already carries, so it ports to the MIQCQP formulations as well as the smooth NLP.

monee.model.phys.nonlinear.smooth.minor_loss_pressure(p_i_pu, p_j_pu, loss_coefficient, signed, mag, diameter_m, fluid_density_kg_per_m3, pressure_ref_pa)

pandapipes heat_exchanger minor loss (no length, no friction factor): \(\Delta p = \zeta \cdot \tfrac{\rho}{2} v^2\). With \(v = \dot m/(\rho A)\) and the signed \(\dot m |\dot m|\) form this is \((p_i-p_j) = -\zeta \cdot \dot m |\dot m| / (2 \rho A^2)\), returned in per-unit. \(\zeta = 0\) pins \(p_i = p_j\) (a lossless injector).

Unlike weymouth_pressure() there is no on_off factor: the only caller (PassiveHeatExchanger) is always on. Thread on_off through here if a switchable minor-loss element is ever added.

monee.model.phys.nonlinear.smooth.signed_psi_breakpoints(diameter_m, roughness_m, dynamic_visc_pas, area, m_max, n_breakpoints=12)

Odd breakpoints for \(\psi(m) = friction(Re(|m|)) \cdot m \cdot |m|\) on \(m \in [-m_{max}, m_{max}]\) - the full Darcy/Weymouth drop term as one spline.

monee.model.phys.nonlinear.smooth.smooth_abs(signed, eps, sqrt_impl=<built-in function sqrt>)

Smooth absolute value \(\sqrt{m^2 + \varepsilon^2}\). \(\varepsilon\) is a small mass-flow scale [kg/s].

monee.model.phys.nonlinear.smooth.smooth_friction_blend(reynolds_scaled, diameter_m, roughness_m, log_impl, exp_impl, re_crit=2300.0, sharpness=400.0)

Darcy friction factor as a single smooth function of Reynolds.

Sigmoid blend of the laminar law 64/Re and the turbulent Swamee-Jain correlation around re_crit - replaces the non-smooth if Re < 2300 switch so IPOPT sees a continuously differentiable friction. reynolds_scaled is Re / REYNOLDS_SCALE (see hyd.REYNOLDS_SCALE).

monee.model.phys.nonlinear.smooth.weymouth_pressure(psq_pu_i, psq_pu_j, p_pu_i, p_pu_j, drop_term, diameter_m, length_m, t_k, compressibility, pressure_ref_pa, on_off, r_specific=504.5, pressure_ambient_pa=0.0)

Signed Weymouth, row-normalised by the per-pipe pressure coefficient.

The raw form \(\Delta(p^2) \cdot C^2 \cdot \text{on\_off} = -\text{drop\_term}\) carries a coefficient \(C^2 \cdot p_{ref}^2 \propto D^5\) on the (dimensionless) squared-pressure difference, which spans ~6 orders across a wide-diameter network - skewing IPOPT’s column scaling and conditioning. Dividing through by that coefficient gives every pipe a unit coefficient on the pressure term (the solution is unchanged):

\[(p_i^2 - p_j^2)_{pu} \cdot \text{on\_off} = -\text{drop\_term} / (C^2 \cdot p_{ref}^2)\]

The \(p_i^2-p_j^2\) uses ABSOLUTE pressures; when the grid carries a nonzero pressure_ambient_pa the node pressures are gauge and the absolute squared difference is formed via gf.abs_psq_diff_pu() (a no-op offset of 0 reproduces the historical absolute form).

MISOCP branch flow

monee.model.phys.misoc.pf.active_power_loss(var_active_power_from, var_active_power_to, var_ell_pu, resistance_r)

monee.model.phys.misoc.pf.gap_expr(var_w_pu_i, var_active_power_ij_pu, var_reactive_power_ij_pu, var_ell_pu, tap=1.0)

monee.model.phys.misoc.pf.reactive_power_loss(var_reactive_power_from, var_reactive_power_to, var_ell_pu, reactance_x)

monee.model.phys.misoc.pf.soc_eq(var_w_pu_i, var_active_power_ij_pu, var_reactive_power_ij_pu, var_ell_pu, tap=1.0)

Exact (non-convex) form of soc_rel(): \(P^2 + Q^2 = (W/tap^2) \cdot ell\). Pins the branch-flow model to the physical surface; requires a global MIQCQP solver (SCIP, Gurobi).

monee.model.phys.misoc.pf.soc_rel(var_w_pu_i, var_active_power_ij_pu, var_reactive_power_ij_pu, var_ell_pu, tap=1.0)

Rotated SOC \(P^2 + Q^2 \le (W/tap^2) \cdot ell\). See soc_rel_lorentz() for the Lorentz reformulation that survives non-convex MIQCP siblings.

monee.model.phys.misoc.pf.soc_rel_lorentz(var_w_pu_i, var_active_power_ij_pu, var_reactive_power_ij_pu, var_ell_pu, var_soc_sum, var_soc_diff, tap=1.0)

Lorentz form of soc_rel(): \(s = (W/tap^2+ell)/2\), \(d = (W/tap^2-ell)/2\), \(P^2 + Q^2 + d^2 \le s^2\) (since \(s^2 - d^2 = (W/tap^2) \cdot ell\)).

monee.model.phys.misoc.pf.voltage_drop(var_w_pu_i, var_w_pu_j, var_active_power_ij_pu, var_reactive_power_ij_pu, var_ell_pu, resistance_r, reactance_x, tap=1.0)

Islanding connectivity

Single-commodity connectivity-flow constraints used for island detection and energization.

monee.model.phys.islanding.angle_lower_bound_energized(theta_i, theta_max, e_i)

theta_i >= -theta_max * e_i (zero angle when de-energised).

monee.model.phys.islanding.angle_upper_bound_energized(theta_i, theta_max, e_i)

theta_i <= theta_max * e_i (zero angle when de-energised).

monee.model.phys.islanding.connectivity_arc_capacity_line(c_ij, y_ij, big_m_conn)

c_ij <= M_conn * y_ij (no flow on open/off branches).

monee.model.phys.islanding.connectivity_arc_capacity_source(c_0s, g_s, big_m_conn)

c_0s <= M_conn * g_s (no flow from disabled grid-forming sources).

monee.model.phys.islanding.connectivity_demand_balance(bus_inflow, bus_outflow, e_i)

Per-bus balance: inflow - outflow == e_i. Each energised bus receives exactly 1 unit of connectivity flow; de-energised buses receive 0.

monee.model.phys.islanding.connectivity_super_source_supply(super_outflow, total_energized_buses)

Super-source supply equals the total number of energised buses: sum_s c_0s == sum_i e_i.

monee.model.phys.islanding.source_reference_angle(theta_s)

Fix the voltage angle at a grid-forming bus: theta_s == 0.