Internal functions

EMB.variables_node(m, 𝒩ᶜ::Vector{<:Compressor}, 𝒯, modeltype::EMB.EnergyModel)

When the node vector is a Vector{<:Compressor} the potential increase (potential_Δ) variables are created for each compressor and timestep.

EMB.variables_flow_resource(m, 𝒩::Vector{<:EMB.Node}, 𝒫::Vector{<:ResourcePressure}, 𝒯, modeltype::EMB.EnergyModel)

Define the potential variables potential_in and potential_out of the nodes if the system has ResourcePressures. The blend composition variables proportion_out are created together with the link variables, as they depend on the network topology.

EMB.variables_flow_resource(m, ℒ::Vector{<:EMB.Link}, 𝒫::Vector{<:ResourcePressure}, 𝒯, modeltype::EnergyModel)
EMB.variables_flow_resource(m, ℒ::Vector{<:EMB.Link}, 𝒫::Vector{<:ResourcePooling}, 𝒯, modeltype::EnergyModel)

Define additional pressure-related variables for links if there are ResourcePressure in the system. For ResourcePooling blends, define flow_component variables on each blend-carrying link.

For the ResourcePooling dispatch, proportion_out[n, t, p] is also created for the nodes that track the blend composition explicitly, see blend_share_plan.

For the ResourcePressure dispatch, the following variables are created:

• link_potential_in[l ∈ ℒ, 𝒯, inputs(l.to)] — inlet potential for each link (all links, all resources of the destination node).
• link_potential_out[l ∈ ℒ, 𝒯, outputs(l.from)] — outlet potential for each link (all links, all resources of the source node).
• has_flow[l ∈ ℒ_pressure, 𝒯] — binary; equals 1 when a link is active (carrying flow). Used as a big-M switch in pressure-balance and capacity constraints to deactivate idle links. Note: has_flow = 1 is not sufficient to guarantee actual flow — it is only a necessary condition.

The binary variable is declared only for ℒ_pressure — the subset of links carrying at least one CompoundResource (ResourcePressure or ResourcePooling{<:ResourcePressure}). Links that carry only non-pressure resources such as ResourceCarrier (e.g., Power) do not require pressure binaries and are excluded to keep the MINLP compact.

EMB.variables_ext_data(m, _::Type{<:BlendStorageData}, 𝒩ˢᵗᵒʳ::Vector{<:EMB.Node}, 𝒯, 𝒫, modeltype::EnergyModel)

Create the component inventory variables of the storages 𝒩ˢᵗᵒʳ with BlendStorageData, indexed by the sub-resources p of their blend:

• stor_level_component[n, t, p] ≥ 0 is the stored amount of p at the end of t.
• stor_level_component_Δ_op[n, t, p] is the rate of change of the stored amount of p.
• stor_level_component_Δ_rp[n, t_rp, p] is the change of the stored amount of p in the representative period t_rp. It is only created if 𝒯 has representative periods.

EMB.constraints_resource(m, n::EMB.Node, 𝒯, 𝒫::Vector{:ResourcePressure}, modeltype::EMB.EnergyModel)
EMB.constraints_resource(m, n::EMB.Node, 𝒯, 𝒫::Vector{:ResourcePooling{Any}}, modeltype::EMB.EnergyModel)  
EMB.constraints_resource(m, n::EMB.Node, 𝒯, 𝒫::Vector{:ResourcePooling{ResourcePressure}}, modeltype::EMB.EnergyModel)

Add blending and/or pressure related constraints to node n based on specific resource types.

• If 𝒫::Vector{<:ResourcePressure} then it adds only pressure related constraints
• If 𝒫::Vector{<:ResourcePooling{Any}} then it adds only blending related constraints
• If 𝒫::Vector{<:ResourcePooling{ResourcePressure}} then it adds both pressure and blending related constraints

Note! The blending constraints for nodes require ℒ:Vector{<:EMB.Link} to be passed as argument. Thus, all of them are defined in constraints_couple_resource() functions.

EMB.constraints_resource(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:ResourcePressure}, modeltype::EMB.EnergyModel)
EMB.constraints_resource(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:ResourcePooling{Any}}, modeltype::EMB.EnergyModel)  
EMB.constraints_resource(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:ResourcePooling{ResourcePressure}}, modeltype::EMB.EnergyModel)

Add blending and/or pressure related constraints to node l based on specific resource types transported through l`.

• If 𝒫::Vector{<:ResourcePressure} then it adds only pressure related constraints
• If 𝒫::Vector{<:ResourcePooling{Any}} then it adds only blending related constraints
• If 𝒫::Vector{<:ResourcePooling{ResourcePressure}} then it adds both pressure and blending related constraints

Note! The blending constraints for nodes require ℒ:Vector{<:EMB.Link} to be passed as argument. Thus, all of them are defined in constraints_couple_resource() functions.

EMB.constraints_couple_resource(m, 𝒩::Vector{<:EMB.Node}, ℒ::Vector{<:EMB.Link}, 𝒫::Vector{<:ResourcePressure}, 𝒯, modeltype::EMB.EnergyModel)
EMB.constraints_couple_resource(m, 𝒩::Vector{<:EMB.Node}, ℒ::Vector{<:EMB.Link}, 𝒫::Vector{<:ResourcePooling{Any}}, 𝒯, modeltype::EMB.EnergyModel)
EMB.constraints_couple_resource(m, 𝒩::Vector{<:EMB.Node}, ℒ::Vector{<:EMB.Link}, 𝒫::Vector{<:ResourcePooling{ResourcePressure}}, 𝒯, modeltype::EMB.EnergyModel)

Add blending and/or pressure related coupling constraints between nodes and links based on specific resource types.

EMB.create_link(m, l::ZeroDrop, 𝒯, 𝒫, modeltype::EMB.EnergyModel)

Dispatched function for setting the flow-balance constraints for a ZeroDrop link.

ZeroDrop is a zero-pressure-drop pipeline: outlet flow equals inlet flow for every resource. It is declared as <: EMB.Link (not <: EMB.Direct) so that the bidirectional mutual-exclusion constraint in constraints_bidirectional_pressure still applies when two anti-parallel ZeroDrop links model a reversible pipeline.

EMB.create_link(m, l::CapDirect, 𝒯, 𝒫::Vector{<:CompoundResource}, modeltype::EMB.EnergyModel)

Dispatched function for setting the constraints for a link of type CapDirect.

EMB.constraints_capacity(m, l::CapDirect, 𝒯, modeltype::EMB.EnergyModel)

Function for creating the constraints on the maximum capacity of a link CapDirect.

EMB.constraints_capacity(m, n::UnitConversion, 𝒯, modeltype::EMB.EnergyModel)

Remove the capacity constraints for UnitConversion nodes. These nodes do not have :cap_use.

EMB.constraints_opex_var(m, n::UnitConversion, 𝒯, modeltype::EMB.EnergyModel)

Remove the opex variable constraints for UnitConversion nodes.

EMB.constraints_opex_fixed(m, n::UnitConversion, 𝒯, modeltype::EMB.EnergyModel)

Remove the opex fixed constraints for UnitConversion nodes.

EMB.constraints_flow_in(m, n::PoolingNode, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Function for creating the constraint on the inlet flow of a PoolingNode. The sum of the flows from all input links must be equal to the used capacity of the node. It differs from generic nodes in that it is not defined as a proportion of the :input and :output fields.

EMB.constraints_flow_in(m, n::DiscretizedPoolingNode, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Same as PoolingNode: the sum of all input flows equals the used capacity of the node. DiscretizedPoolingNode accepts mixed resource inputs (e.g. CH4 and Blend) and the total throughput is the sum over all inlet resources, not a per-resource equality. Without this dispatch the default NetworkNode implementation would generate one constraint per input resource (flow_in[p] == cap_use), incorrectly forcing all resource flows equal.

EMB.constraints_flow_in(m, n::TransitNode, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Same as PoolingNode: the sum of all input flows equals the used capacity of the node. TransitNode accepts mixed resource inputs (H2, CH4, Blend) and treats their total as the throughput, identical to PoolingNode.

EMB.constraints_flow_in(m, n::BlendingNode, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Same as PoolingNode and TransitNode: the sum of all input flows equals the used capacity. BlendingNode accepts mixed resource inputs (H2, CH4, Blend) in variable proportions.

EMB.constraints_flow_in(m, n::DeblendNode, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Volumetric mass balance for DeblendNode: all incoming blend equals the outgoing blend plus the H₂ extracted from the blend. Extraction > delivery when psa_efficiency < 1:

flow_in[Blend, t] == flow_out[Blend, t] + Σ_{p ∈ extracted} flow_out[p, t] / η

where η = psa_efficiency(n). The PSA consumes (1−η)/η × flow_out[H₂] internally.

EMB.constraints_flow_in(m, n::UnitConversion, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Remove the flowin constraints for UnitConversion nodes. These nodes do not use :capuse, instead flows are calculated based on their data extension in EMB.constraintsextdata.

EMB.constraints_flow_out(m, n::DeblendNode, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Set cap_use equal to the total H₂ extracted from the blend (i.e. delivered ÷ η). This links the investable capacity (cap_inst) to the PSA extraction rate via the standard EMB constraint cap_use ≤ cap_inst.

cap_use[n, t] == Σ_{p ∈ extracted} flow_out[p, t] / η

where η = psa_efficiency(n). The PSA is sized by its extraction throughput, not by the net delivery, so capacity accounts for the internally consumed fraction.

The outgoing blend flow is NOT constrained here — it is determined by the volumetric mass balance in constraints_flow_in.

EMB.constraints_flow_out(m, n::UnitConversion, 𝒯::TimeStructure, modeltype::EMB.EnergyModel)

Remove the flowout constraints for UnitConversion nodes. These nodes do not use :capuse, instead flows are calculated based on their data extension in EMB.constraintsextdata.

constraints_balance_pressure(m, n::EMB.Node, 𝒯, 𝒫::Vector)
constraints_balance_pressure(m, n::EMB.NetworkNode, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_balance_pressure(m, n::Compressor, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_balance_pressure(m, n::Compressor, 𝒯, 𝒫::Vector{<:ResourcePressure})
constraints_balance_pressure(m, n::DeblendNode, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_balance_pressure(m, n::PoolingNode, 𝒯, 𝒫::Vector{<:ResourcePressure})
constraints_balance_pressure(m, n::PoolingNode, 𝒯, 𝒫::Vector{<:ResourcePooling{<:ResourcePressure}})
constraints_balance_pressure(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:CompoundResource})

Set internal balance pressures between potential_in and potential_out in Nodes n and Links l`. The balance will depend on type of nodes:

• For NetworkNode nodes, the inlet and outlet potentials are equal.
• For Compressor nodes, the inlet potential is lower than or equal to the outlet potential.
• For DeblendNode nodes, the Blend output equals the Blend input (no pressure drop through the separator); the extracted H₂ output pressure is unconstrained here.
• For PoolingNode nodes, the inlet potential of all input resources is equal to the outlet potential of all output resources.
• For Links, the inlet potential is higher than or equal to the outlet potential. if there is no flow through l, both potentials are set to zero.

Note: Sinks and Source nodes do not have internal pressure balances, as their potentials are only defined at inlet or outlet respectively.

constraints_balance_pressure(m, n::Compressor, 𝒯, 𝒫::Vector{<:ResourcePressure})

For a Compressor carrying a pure ResourcePressure gas (e.g. H₂ without blending), the outlet potential must be at least as large as the inlet potential — the compressor only raises pressure, never lowers it.

This mirrors constraints_balance_pressure(m, n::SimpleCompressor, ..., CompoundResource). The equality potential_out = potential_in + potential_Δ is enforced in constraints_pressure_couple.

constraints_pressure_bounds_element(m, x, 𝒯, 𝒫::Vector{<:CompoundResource})

This function calls subfunctions to set the pressure bounds for each AbstractPressureData assigned to element x.

constraints_pressure_bounds(m, n::Node, data::T, 𝒯, 𝒫::Vector{<:CompoundResource}) where {T<:PressureData}
constraints_pressure_bounds(m, n::Sink, data::T, 𝒯, 𝒫::Vector{<:CompoundResource}) where {T<:PressureData}
constraints_pressure_bounds(m, l::Link, data::T, 𝒯, 𝒫::Vector{<:CompoundResource}) where {T<:PressureData}

Set the pressure limits according to the PressureData data assigned to Node n and Link l. If n is a Sink, the limits are applied to the inlet potential. Otherwise, the limits are applied to the outlet potential. For Links, the limits are applied to the outlet potential only if the link has flow.

The data can be of type:

• MaxPressureData, which will set a maximum pressure bound.
• MinPressureData, which will set a minimum pressure bound.
• FixPressureData, which will set an equality pressure bound.

constraints_pressure_couple(m, n::EMB.Source, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_pressure_couple(m, n::EMB.NetworkNode, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_pressure_couple(m, n::Compressor, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_pressure_couple(m, n::Compressor, ℒ, 𝒯, 𝒫::Vector{<:ResourcePressure})
constraints_pressure_couple(m, n::PoolingNode, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_pressure_couple(m, n::DiscretizedPoolingNode, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_pressure_couple(m, n::TransitNode, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})
constraints_pressure_couple(m, n::EMB.Sink, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})

Constraints coupling node potentials to the potentials of their connected links.

Each dispatch sets the node's potential_in and potential_out equal (within a big-M relaxation) to the corresponding link_potential_out and link_potential_in of connected pressure-carrying links. The coupling logic varies by node type:

• EMB.Source: outlet potential equals each outgoing link's inlet potential (equality, no big-M needed since sources always inject).
• EMB.NetworkNode, PoolingNode, TransitNode, EMB.Sink: uses has_flow big-M for both the upper and lower bound on potential_in[n]. When has_flow=1, both bounds bind, forcing potential_in[n] = link_potential_out[l] for each active link. When multiple links are active, their outlet potentials are all forced equal to the node inlet potential — correct junction physics assuming pressure reduction valves reduce any higher-pressure incoming flows to the node pressure.
• Compressor: same as NetworkNode, but incoming Power links are filtered out before applying any pressure coupling (pressure coupling applies only to the gas resource, not the power resource).

Link filtering

In all dispatches, ℒᵗᵒ and ℒᶠʳᵒᵐ are filtered to include only links that carry at least one CompoundResource. This excludes Power-only links (e.g., compressor power inputs, output links from unit-conversion nodes) which have no link_potential_* variables and whose has_flow binary is not declared.

Pressure coupling

For each active incoming link (has_flow=1), the node inlet potential is bounded above by the link outlet potential via a has_flow-gated big-M constraint:

potential_in[n] ≤ link_potential_out[l] + 1e4*(1 - has_flow[l])

When has_flow=1 this forces potential_in[n] ≤ link_potential_out[l]. When has_flow=0 the bound is inactive. With multiple active incoming links the node inlet potential is bounded above by the minimum of all active link outlet potentials. The optimizer naturally drives the node pressure up to that minimum (to maximise downstream flow capacity), giving "lowest incoming pressure" semantics without a separate lower bound constraint.

No explicit lower bound is added for potential_in[n] from incoming links; the lower bound comes from downstream constraints (Weymouth flow requirements, MinPressureData).

constraints_pressure_couple(m, n::Compressor, ℒ, 𝒯, 𝒫::Vector{<:ResourcePressure})

Pressure coupling for a Compressor carrying a pure ResourcePressure gas (e.g. pure H₂, no blending). Mirrors the CompoundResource dispatch but filters for ResourcePressure links, excluding Power-only links that carry no pressure potential.

Enforces:

• potential_in[n] ≤ link_potential_out[l_to] + BigM×(1−has_flow[l_to]) (inlet coupling)
• potential_out[n] = potential_in[n] + potential_Δ[n] (pressure rise)
• potential_Δ[n] ≤ max_incr_potential[n] (compressor limit)
• link_potential_in[l_from] = potential_out[n] (bigM equality) (outlet coupling)

constraints_pressure_couple(m, n::TransitNode, ℒ, 𝒯, 𝒫::Vector{<:CompoundResource})

Couple inlet/outlet potentials of links to the potential_in/potential_out of a TransitNode. Uses has_flow big-M for both upper and lower bounds, coupling each active link's outlet potential to the node inlet potential.

constraints_flow_capacity(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:CompoundResource})

Limit flow to the installed capacity and deactivate it when has_flow[l, t] is zero.

The installed-capacity constraint remains time-dependent. The has_flow constraint uses a finite upper bound on the maximum installed capacity as its big-M. If no finite bound exists, the constraint can only be omitted when PressureLinkData is present. Taylor cuts then force zero flow at zero potentials, and constraints_flow_pressure adds a PWA-derived big-M for blend links.

_maximum_link_capacity(m, l, 𝒯)

Return a finite upper bound on the installed capacity of link l, or nothing.

Investment links use the upper bounds on link_cap_current. Fixed-capacity links use fixed values or upper bounds on link_cap_inst. The physical capacity is the final fallback.

constraints_flow_pressure(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:ResourcePressure})
constraints_flow_pressure(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:ResourcePooling{<:ResourcePressure}})
constraints_flow_pressure(m, l::EMB.Link, 𝒯, 𝒫::Vector{<:Resource})
constraints_flow_pressure(m, l::EMB.Direct, 𝒯, 𝒫::Vector)

Set Weymouth constraints in link l, calculating the flow based on the pressure difference between inlet and outlet.

The Weymouth equation is approximated using the first-order Taylor expansion when the resource is a ResourcePressure. For a ResourcePooling of ResourcePressure sub-resources (blend links), a Piecewise Affine Approximation (PWA) in (p_in, p_out, proportion_H2) is used when a BlendData object is attached to the link. If no BlendData is present — for example when blend: nothing is set in YAML for a link that still carries a blend resource due to its node types — the function falls back to the 2D Taylor approximation at the reference (pure-component) composition encoded in PressureLinkData. This fallback is conservative: it uses the capacity at zero H₂ fraction, which is the minimum possible capacity (since adding H₂ lightens the gas and increases throughput).

No-op methods are defined for EMB.Direct links and plain Resource vectors where no Weymouth constraint is needed.

constraints_pwa(m, l::Link, p_blend::ResourcePooling, p_track::ResourcePressure, 𝒯, plane, pwa::PWAFunc)

constraints_energy_potential(m, n::SimpleCompressor, 𝒯, 𝒫, modeltype::EMB.EnergyModel)
constraints_energy_potential(m, n::EMB.Node, 𝒯, 𝒫, modeltype::EMB.EnergyModel)

Sets the relationship between energy needs and pressure increase, or any other parameter in Compressor n that determines its energy consumption. Skip if the node is not a type Compressor.

If n is a SimpleCompressor, the energy consumption will be calculated as the product of the flow through the compressor, and the energy input required per unit of flow. This is defined with the default constraint constraints_flow_in().

Note! If new Compressor types are created with different relationships between energy flow and pressure increase, this function should be updated to include the new type and relationship.

constraints_energy_potential(m, n::FlowPressureCompressor, 𝒯, 𝒫, modeltype::EMB.EnergyModel)

For a FlowPressureCompressor, the energy input must cover both the flow-dependent and the pressure-dependent components of compression work:

flow_in[n, t, PowerRes] ≥ flow_weight(n, t) × cap_use[n, t]
                          + pressure_weight(n, t) × potential_Δ[n, t]

The energy resource PowerRes is identified as the input resource that is NOT in the gas outputs (i.e., setdiff(inputs(n), outputs(n))), exactly as in SimpleCompressor.

" constraintsbidirectionalpressure(m, l::EMB.Link, ℒ, 𝒯, 𝒫) constraintsbidirectionalpressure(m, l::EMB.Direct, ℒ, 𝒯, 𝒫)

Ensure that parallel links defining bidirectionality cannot have flow at the same time. The link without flow will automatically have zero linkinpotential and linkoutpotential due to the constraints_balance_pressure constraints.

blend_share_plan(𝒩, ℒ, 𝒫)

Decide which nodes track the blend composition explicitly with proportion_out.

The component flows on the links give the composition wherever blend mixes, so a share is only needed where a composition must be shared or read:

• Composition sources (_is_composition_source): pooling nodes and blend storages in mode :prescribed or :fully_mixed, whose outflows all carry the composition of the node.
• Weymouth pipes with BlendLinkData, whose PWA depends on the share of the tracked resource. The share is read from
  1. the inflow node, if it is a composition source, or else
  2. the first pooling node downstream, if the pipe is its only blend inflow, possibly through passthrough nodes with a single blend inflow and outflow. The composition is then the same along the whole path. Or else
  3. the inflow node, if all gas reaching it comes from composition sources through passthrough nodes, which then copy the share (see _constraints_blend_passthrough). The copy requires simultaneously active inflows of these nodes to have the same composition.

All other nodes, e.g. sinks, merges and passthrough nodes off such paths, only balance the component flows. A pipe without a share raises an ArgumentError.

Returns a NamedTuple with

• provenance: for each blend, the classification of each node as proportion-capable (_PROPORTION_CAPABLE), component-flow-only (_COMPONENT_FLOW_ONLY) or both,
• share_nodes: the nodes with proportion_out,
• copy_nodes: for each blend, the passthrough nodes copying the share from upstream,
• pwa_reader: for each (link, blend) pair with BlendLinkData, the node whose share the PWA uses.

validate_blend_topology(𝒩, ℒ, 𝒫; switched_links = nothing)

Reject blend topologies that the component-flow formulation can not represent, and return the provenance of blend_share_plan:

• A Storage of a blend requires BlendStorageData.
• A node with several outgoing blend links must give them the same composition. Pooling nodes and blend storages with a shared composition do so. A passthrough node (see _is_blend_passthrough) is only accepted if all its outgoing blend links are in switched_links, the links with a has_flow binary; _constraints_passthrough_single_outflow then allows at most one of them to carry flow in each period. Other nodes, e.g. BlendingNode, are rejected.
• Each link with BlendLinkData needs a share, see blend_share_plan.

constraints_blend_proportion(m, n, ℒ, 𝒯, 𝒫)

Set the proportion_out variable and any associated constraints for node n.

• PoolingNode (bilinear/MINLP): for each blend output link, flow_component[l_out, t, p] = proportion_out[n, t, p] × link_in[l_out, t, blend]

• DiscretizedPoolingNode (piecewise linear/MILP): same semantics as PoolingNode but uses PiecewiseLinearOpt.piecewiselinear to approximate the bilinear product with a piecewise linear function. Controlled by DiscretizationData.

• TransitNode and generic blend-passthrough NetworkNode (linear): proportion_out[n, t, p] = proportion_out[upstream_node, t, p] where upstream_node is the node at the other end of the primary incoming blend link.

• All other node types: no-op.

Normalisation (Σ proportion_out = 1) is applied globally by constraints_component_blend and is not repeated here.

constraints_blend_proportion(m, n::BlendingNode, ℒ, 𝒯, 𝒫)

No-op for BlendingNode. The blend composition on the single outgoing link is fully determined by the linear component mass balance in constraints_node_component_balance. No bilinear proportion_out × link_in constraint is needed or generated.

Quality bounds are handled separately by constraints_quality_blend using direct linear flow-component bounds.

constraints_blend_proportion(m, n::DiscretizedPoolingNode, ℒ, 𝒯, 𝒫)

MILP-compatible mixing constraint for DiscretizedPoolingNode.

Dispatches to one of three sub-implementations based on the method field of DiscretizationData:

• :pwl — 2D piecewise-linear approximation via PiecewiseLinearOpt.piecewiselinear (the original formulation). Bilinear flow_component = proportion_out × link_in is approximated on a K×K grid. Uses ⌈log₂(K−1)⌉ binary variables per product with the :Logarithmic method.

• :bigm — standard BigM MILP formulation. ONE set of K binary variables δ[t, k] per time step (shared across all outgoing links). The composition proportion_out is set directly to a discrete level α_k = (k-1)/(K-1), which is a linear JuMP expression. This also makes downstream pressure constraints MILP-linear. N-1 + 1 = N outgoing link flows are fully determined; the last component flow is inferred from the node mass balance.

• :bigm_log — logarithmic BigM MILP formulation. Uses n_bits = ⌈log₂(K)⌉ binary variables per time step via binary encoding. The actual number of levels is K_actual = 2^n_bits ≥ n_points. For n_points = 10: 4 bits → 16 levels.

See _constraints_dpn_pwl, _constraints_dpn_bigm for details.

constraints_blend_proportion(m, n::EMB.Storage, ℒ, 𝒯, 𝒫)

Set the composition of the gas discharged from a blend storage n according to the mode of its BlendStorageData:

• :relaxed: no constraint, the component inventories bound the discharge.
• :prescribed: flow_component[l, t, p] = q_p(t) × link_in[l, t, blend] on each discharge link l, and proportion_out[n, t, p] = q_p(t).
• :fully_mixed: flow_component[l, t, p] = proportion_out[n, t, p] × link_in[l, t, blend] on each discharge link l and stor_level_component[n, t, p] = proportion_out[n, t, p] × stor_level[n, t] (bilinear). The composition is that of the stored gas at the end of period t.

No-op for storages without BlendStorageData.

constraints_quality_blend(m, n::BlendingNode, ℒ, 𝒯, 𝒫)

Enforce blend composition bounds for BlendingNode using linear flow-component bounds on the single outgoing blend link.

Since BlendingNode has no proportion_out variable, bounds from RefBlendData are expressed directly in terms of flow_component and link_in:

flow_component[l_out, t, H₂] ≤ max_h2 × link_in[l_out, t, Blend]   (LINEAR)
flow_component[l_out, t, H₂] ≥ min_h2 × link_in[l_out, t, Blend]   (LINEAR)

With big-M relaxation for has_flow = 0 (same pattern as the generic proportion_out quality dispatch):

flow_component[l_out, t, H₂] ≤ max_h2 × link_in[l_out, t, Blend] + M × (1 - has_flow)

where M = 1 (the maximum possible flowcomponent / linkin ratio).

No-op when no RefBlendData is attached.

constraints_quality_blend(m, n, ℒ, 𝒯, 𝒫)

Enforce the blend composition bounds from RefBlendData attached to node n.

Since proportion_out[n] is now defined for all nodes with blend in inputs or outputs (including Sink and UnitConversion), the bounds are applied directly to proportion_out[n]. This is always linear.

No-op when no RefBlendData is attached.

constraints_node_component_balance(m, n, ℒ, 𝒯, 𝒫)

Linear per-component mass balance at a node n.

For each subresource p of every blend: Σ{lout blend} flowcomponent[lout, t, p] == flowin[n, t, p] (if p ∈ inputs(n), else 0) + Σ{lin blend} flowcomponent[l_in, t, p]

This is a linear identity that connects the raw inflow of p (e.g., pure H2 injected at a PoolingNode) with the component-tracked flows on the blend links.

No-op for Source nodes (no blend output) and Sink nodes (no blend output to balance).

constraints_node_component_balance(m, n::DeblendNode, ℒ, 𝒯, 𝒫)

Component balance for DeblendNode:

• Extracted sub-resource (e.g. H₂ that also appears in outputs(n)): flow_component[l_out, t, H₂] = flow_component[l_in, t, H₂] - flow_out[n, t, H₂]

• Passthrough sub-resources (e.g. CH₄): flow_component[l_out, t, CH₄] = flow_component[l_in, t, CH₄]

constraints_node_component_balance(m, n::EMB.Storage, ℒ, 𝒯, 𝒫)

Set the rate of change of the component inventories of a blend storage n from the component flows on its charge links ℒⁱⁿ and discharge links ℒᵒᵘᵗ:

stor_level_component_Δ_op[n, t, p] =
    Σ_{l ∈ ℒⁱⁿ} flow_component[l, t, p] - Σ_{l ∈ ℒᵒᵘᵗ} flow_component[l, t, p]

Together with the level balance in EMB.constraints_ext_data and the bound stor_level_component ≥ 0, the storage can not discharge more of a component than was charged and stored.

No-op for storages without BlendStorageData.

constraints_link_component_balance(m, ℒ, 𝒯, 𝒫)

For every link carrying a ResourcePooling blend, enforce:

link_in[l, t, blend] == Σ_p flow_component[l, t, p]

so that the sum of all tracked component flows equals the total blend flow.

constraints_component_blend(m, 𝒩, ℒ, 𝒯, 𝒫)

Entry point for the component-flow pooling formulation.

Adds, in order:

1. The checks of validate_blend_topology and, at passthrough nodes with several outgoing blend links, at most one active outflow.
2. Linear link component balances (link_in[blend] == Σ flow_component[p]).
3. Linear per-component node mass balances.
4. The composition constraints of the nodes that track it (see blend_share_plan): bilinear or discretized mixing at pooling nodes, blend storage, and the share copied by passthrough nodes from upstream.
5. Linear quality bounds from RefBlendData.
6. Deferred unit conversion and deblending constraints on component flows.
7. Normalisation Σ proportion_out[n,t,p] = 1 at the passthrough nodes copying a share.

EMB.constraints_ext_data(m, n::EMB.Storage, 𝒯, 𝒫, modeltype::EnergyModel, data::BlendStorageData)

Add the level balance of each component inventory of the blend storage n:

stor_level_component[n, t, p] = previous level + stor_level_component_Δ_op[n, t, p] × duration(t)

The previous level follows the same rules as the total level of the storage in EnergyModelsBase, including the cyclic constraints of the storage behaviour. The component levels add up to the total level:

Σ_p stor_level_component[n, t, p] = stor_level[n, t]

The rate of change stor_level_component_Δ_op is set from the component flows on the links in constraints_node_component_balance, as the link variables do not exist yet when this function is called.

EMB.constraints_ext_data(m, n::Node, 𝒯, 𝒫, modeltype::EnergyModel, data::FlowToEnergyData)

Function to convert flow rates to energy rates (e.g., Sm3/d to MWh/d). The time basis of the flow rates is preserved in the conversion.

When the node receives a ResourcePooling blend input with a per-sub-resource LHV dict, the constraint would be quadratic (LHV(proportion_out) × flow_in). In that case this function defers constraint creation — the linearised form is added later in constraints_component_blend once flow_component link variables exist:

flow_out = Σ_sub  lhv_sub × flow_component[l_in, t, sub]

For scalar LHV (single-resource input), the constraint is created here directly:

flow_out == LHV × flow_in

The scalar constraint is not applied to the CO₂ output of an EmissionsConversionData.

EMB.check_node_data(n::EMB.Node, data::BlendStorageData, 𝒯, modeltype::EMB.EnergyModel, check_timeprofiles::Bool)

Check that a BlendStorageData is valid:

• the node is a Storage with a Cyclic behaviour that stores the blend of the data,
• the node has only one BlendStorageData,
• the sub-resources of the blend are not inputs or outputs of the node, and
• in mode :prescribed, the composition is given for all sub-resources, lies in [0, 1] and sums to 1 in every operational period.

nodes_upstream_of(n::Node, ℒ::Vector{<:Link})

Tracks all nodes associated with a given node n through the links in ℒ. We refer to associated to all the nodes that lead to n following the direction of the links.

Links with a negative integer id are treated as reverse CapDirect links (bidirectional models use id = -(forward_id) for reverse directions) and are skipped during traversal so that upstream source sets remain directionally correct.

sources_upstream_of(n::Node, ℒ::Vector{<:Link})
sources_upstream_of(n::EMB.Node, ℒ::Vector{<:EMB.Link}, resources::Vector{EMB.Resources})

Tracks all nodes associated to n and filter by source nodes. If resources is provided, only sources which output any of the resources in sub_res are returned.

get_links_to_node_blend(n::Node, ℒ::Vector{<:EMB.Link}, sub_res, blend)

Gets the links into node n which transport any of the resources in sub_res or the blend resource.

define_points_curve(x1, x2, x3)

Defines the points (inlet and outlet pressures and proportion) for the surface for the PWA.

res_types(𝒫::Array{<:Resource})

Return the unique resource types in an Array of resources 𝒫.

res_types_seg(𝒫::Array{<:Resource})

Return a Vector-of-Vectors of resources segmented by the sub-types.

Get optimizer; error if not set.

link_data(l::CapDirect)

Returns the [ExtensionData] array of link l.

It overwrites the EMB.link_data(l::Link) method, which returns an empty ExtensionData vector.

normalised_weymouth(weymouth, molar_fraction_hydrogen)

Calculate the normalised flow constant with respect to the specific gravity using specific operating points.

calculate_flow_to_approximate(constant, x1, x2, x3)

Calculates the flow of gas with the Weymouth equation using the normalised weymouth constant.

Typically, the constant in the Weymouth equation depends on the specific gravity of the gas. Here, we need instead a normalised it with respect to the specific gravity. This allows to calculate the flows considering the different proportions of the components.

Variables

• weymouth_ct::Float64 -> Normalised Weymouth constant
• x1::Float64 -> Inlet pressures
• x2::Float64 -> Outlet pressures
• x3::Float64 -> Proportion of tracking component
• molmass_other::Float64 -> Molar mass of the other component in the blend (e.g., for methane is 16.042 g/mol)
• molmass_track::Float64 -> Molar mass of the tracking component (e.g., for hydrogen is 2.016 g/mol)

get_step_pressure()

Return the current pressure step (in the same units as the inputs, e.g. bars) used when building piecewise-affine approximations of the Weymouth relation.

get_pwa(data_pressure::PressureLinkData, data_blend::BlendData, optimizer; resolution_prop=0.01)
get_pwa(l::EMB.Link, optimizer; resolution_prop=0.01)

Generate or retrieve the cached PWA (piecewise-affine approximation) of the Weymouth flow equation for a link with pressure and blending data.

The approximation is a Concave piecewise-affine function of three inputs: (p_in, p_out, proportion_tracked). It is generated by fitting planes to exact Weymouth evaluations on a pressure grid, with zero-flow anchor points at (0, 0, α) added to suppress infeasible positive intercepts.

Results are cached to a JSON file via the Scratch.jl scratch space (keyed on the fit parameters). See get_input_fn, read_from_json, and write_to_json.

The fit is seeded from the cache key, so the same data gives the same PWA in every process. Fitted and cached PWAs are checked with _pwa_fit_ok. A fit that underestimates the exact flow by more than 10 % of the maximum flow is refitted (up to _PWA_FIT_ATTEMPTS times, with a different seed each time) and never cached; if no attempt passes, an error is thrown. An ArgumentError is thrown if the pressure step (see set_step_pressure!) leaves fewer than three pressure levels in the range of the link.

If get_clamp_pwa_intercepts returns true (the default), plane intercepts are clamped via _clamp_pwa_intercepts after fitting to ensure feasibility on inactive reverse links.

get_LHV(data::FlowToEnergyData)
get_LHV(data::FlowToEnergyData, p::EMB.Resource)

Collects the resources if Dict{<:Resource,<:Real} (used for ResourcePooling) or the LHV if is a Real (for other types of resources) If the resource p is specified, retrieves the LHV for that specific resource from the dictionary.

resource_lhv(n::Node, 𝒫ˡʰᵛ::Vector{ResourcePooling}, data::FlowToEnergyData)
resource_lhv(n::EMB.Node, 𝒫ˡʰᵛ::Vector{Resource}, data::FlowToEnergyData)

Function to calculate the LHV of the input resource for a given node with FlowToEnergyData. If the input resource is of type ResourcePooling, the LHV is calculated as the weighted average of the LHV of the individual resources in the blend. If the input resource is any other type of Resource, the LHV is simply retrieved from the data.

get_input_fn(x, y)

Generate a unique filename in the cache directory for the given input parameters x and y.

read_from_json(fn)

Reads the PiecewiseAffineApprox.PWAFunc objects from a JSON file fn.

write_to_json(fn, pwa)

Writes the PiecewiseAffineApprox.PWAFunc object pwa to a JSON file fn.

_constraints_flow_pressure_taylor(m, l, 𝒯, 𝒫, pressure_data)

Add first-order Taylor linearisation cuts of the Weymouth equation for link l.

Generates 149 linearisation points along the pressure envelope defined by pressure_data (from PIn down to POut) and adds one constraint per point, per time period, and per resource in 𝒫:

link_in[l,t,p] ≤ √W · [∂f/∂p_in · link_potential_in[l,t,p] + ∂f/∂p_out · link_potential_out[l,t,p]]

where W is the Weymouth constant and the partial-derivative coefficients are evaluated at each linearisation point. Called by constraints_flow_pressure for ResourcePressure links and as a fallback for blend links without BlendData.

_clamp_pwa_intercepts(pwa::PWAFunc{Concave, D}) where {D}

Return a copy of pwa with all plane intercepts clamped to β ≤ 0.

For a concave PWA, the per-plane constraint added by PiecewiseAffineApprox.constr is:

z ≤ -dot(plane.α, x) - plane.β

At x = (0, 0, 0) (link potentials forced to zero by constraints_balance_pressure when a reverse link is inactive): z ≤ -plane.β. With z = 0 (no flow on the inactive link), this becomes 0 ≤ -plane.β. If any plane has plane.β > 0 due to extrapolation outside the fitting domain, this constraint is infeasible.

Clamping sets β_k ← min(0, β_k) for each plane. This:

• Makes 0 ≤ -β_k always satisfiable (since -β_k ≥ 0 after clamping).
• Preserves the outer approximation: decreasing β_k increases -β_k, loosening the upper bound and keeping approx(x) ≥ f(x) on the fitting domain.

Called by get_pwa when get_clamp_pwa_intercepts() returns true (the default). Use set_clamp_pwa_intercepts! to control this behaviour.

_constraints_blend_passthrough(m, n, ℒ, 𝒯, 𝒫)

Linear proportion passthrough for nodes that carry blend unchanged (no injection). Only the passthrough nodes that blend_share_plan selects to copy the share, i.e. those between composition sources and a link with BlendLinkData, get these constraints.

Single-link case (no circularity possible): When a node has exactly one incoming blend link, a hard equality is used: proportionout[n, t, p] == proportionout[upstream, t, p]

Multi-link case (bidirectional / cyclic networks): When a node has two or more incoming blend links, a big-M formulation conditioned on has_flow is used instead of a hard equality. For each incoming link l_up:

proportion_out[n, t, p] - proportion_out[upstream, t, p] ∈ [-(1 - has_flow), (1 - has_flow)]

When has_flow[l_up, t] = 1 (link is active) this collapses to a hard equality, propagating the upstream composition. When has_flow[l_up, t] = 0 (link is inactive) the constraint is trivially satisfied and proportion_out[n] is not pinned to that (inactive) upstream.

This matters in bidirectional networks: a pair {A→B, B→A} carries has_flow[A→B] + has_flow[B→A] ≤ 1, guaranteeing at most one direction is active. Without conditioning, iterating both links would add proportion[n] = proportion[A] and proportion[n] = proportion[B] — contradictory if A and B carry different compositions, and circular (free) if they share the same cycle. Conditioning on has_flow breaks the contradiction: only the active-direction link constrains the proportion.

proportion_out ∈ [0, 1], so M = 1 is the tightest valid big-M coefficient.

Fallback: if has_flow is not in the model (networks without ResourcePressure), hard equalities are used in all cases — safe for the acyclic / unidirectional topologies those models represent.

A copy node may have several simultaneously active inflows (a merge). The constraints then require their compositions to be equal, which restricts the solution; place a pooling node at such merges instead.

_blend_provenance(𝒩, blend, incoming)

Classify each node with blend as output as proportion-capable (all its gas comes from composition sources), component-flow-only (its composition is only known through component flows), or both. The fixed-point propagation through passthrough nodes is cycle-safe, including networks with anti-parallel links.

_downstream_share_anchor(l, sub, incoming, outgoing)

Return the pooling node downstream of link l whose share equals the composition on l, or nothing. Starting at l.to, the path may pass through passthrough nodes with exactly one blend inflow and outflow. Each node on the path, including the pooling node, must have the previous link as its only blend inflow and no inputs of the sub-resources sub, so that no other gas mixes in.

_constraints_dpn_bigm(m, n, ℒ, 𝒯, 𝒫, K, max_prop; log_encoding=false)

BigM MILP formulation for DiscretizedPoolingNode using the composition approach.

Composition Approach

Rather than per-link split fractions, a SINGLE set of binary variables selects the composition (proportion) level for the pooling node at each time step. This is shared across ALL outgoing links — the same composition applies to all downstream pipelines, which is the physical interpretation of a pooling node.

max_prop limits the proportion axis to [0, max_prop]. All K levels span this range, so every level is within the feasible H₂ fraction bound. With max_prop = 1.0 (default), levels span [0%, 100%]; with max_prop = 0.20, they span [0%, 20%].

Standard BigM (log_encoding = false):

• K binary variables δ[t, k] per time step, exactly one active: Σ_k δ[t,k] = 1
• K discrete levels α_k = (k-1)/(K-1) × max_prop for k=1..K
• proportion_out[n, t, sub[1]] = Σ_k α_k · δ[t,k] ← linear!
• proportion_out[n, t, sub[2]] = 1 - proportion_out[n, t, sub[1]]
• For each outgoing link l and level k:
  • w[l,t,k] approximates δ[t,k] · link_in[l,t,blend] via BigM:
    • w ≤ F_max · δ[t,k]
    • w ≥ link_in - F_max · (1 - δ[t,k])
    • 0 ≤ w ≤ link_in
  • flow_component[l,t,sub[1]] = Σ_k α_k · w[l,t,k]
  • flow_component[l,t,sub[2]] = link_in[l,t,blend] - flow_component[l,t,sub[1]]

Logarithmic BigM (log_encoding = true):

• n_bits = ⌈log₂(K)⌉ binary variables b[t, j] per time step
• Actual levels: K_actual = 2^n_bits (rounded up from K)
• Level index: k = Σ_j 2^(j-1) · b[t,j] (binary encoding)
• proportion_out[n, t, sub[1]] = Σ_j 2^(j-1)/(K_actual-1) × max_prop · b[t,j] ← linear!
• For each outgoing link l and bit j:
  • w_j[l,t] approximates b[t,j] · link_in[l,t,blend] via BigM
  • flow_component[l,t,sub[1]] = Σ_j 2^(j-1)/(K_actual-1) × max_prop · w_j[l,t]

Why This Is Equivalent to the User's N-1 Split Fraction Description

For a 2-component blend and N outgoing links: the composition approach uses ONE binary set (not N-1) because all N links share the same pooled composition. The N-th link's component flows are then inferred by the node mass balance (already a linear constraint).

Restriction

Currently supports only 2-component blends (M = 2). Asserts length(sub) == 2.

_constraints_dpn_pwl(m, n, ℒ, 𝒯, 𝒫, K, max_prop)

Original PWL formulation for DiscretizedPoolingNode using PiecewiseLinearOpt.

For each outgoing blend link l, subresource p, and time t, approximates the bilinear product proportion_out[n,t,p] × link_in[l,t,blend] on a K×K grid using the :Logarithmic method (⌈log₂(K−1)⌉ binary variables per product).

max_prop limits the proportion axis to [0, max_prop] instead of [0, 1], placing all K breakpoints within the physically feasible H₂ fraction range.

_constraints_deblend_extraction(m, 𝒩, ℒ, 𝒯, 𝒫)

For each DeblendNode in 𝒩, add the fraction-extraction upper bound on the H₂ extracted from the blend (delivered ÷ η):

flow_out[n, t, p] / η ≤ max_extract_fraction × Σ_{l_in} flow_component[l_in, t, p]

where η = psa_efficiency(n). This must be called after flow_component variables have been created.

_constraints_unit_conversion_blend_lhv(m, 𝒩, ℒ, 𝒯, 𝒫)

Add the linearised FlowToEnergyData conversion constraint for every UnitConversion node whose input is a ResourcePooling blend and whose FlowToEnergyData specifies per-sub-resource LHV values.

This must run after link variables (flow_component) have been declared, which is why it cannot be placed in constraints_ext_data (the node constraint phase).

For each qualifying node the constraint added is:

flow_out[n, t, p_out] == Σ_{l_in} Σ_sub  lhv_sub × flow_component[l_in, t, sub]

which is linear because flow_component[l_in, t, sub] is an explicit variable (it already encodes proportion_out × flow_in).

Outputs that appear as keys in an EmissionsConversionData on the same node are skipped here — they are handled separately by _constraints_unit_conversion_co2_output.

_constraints_unit_conversion_co2_output(m, 𝒩, ℒ, 𝒯, 𝒫)

Add the linearised CO₂ emission output constraint for every UnitConversion node that has an EmissionsConversionData attached and a blend (ResourcePooling) input.

Must run after flow_component link variables exist (i.e., in the deferred blend constraint phase, called from constraints_component_blend).

For each qualifying node and each resource p whose value appears as an output in the node's output dict AND whose key is in EmissionsConversionData.co2_intensity, the constraint is:

flow_out[n, t, p] == Σ_{l_in} Σ_sub  co2_intensity[sub] × flow_component[l_in, t, sub]

Units: flow_component [MSm³/d], co2_intensity [tCO₂/MSm³] → flow_out [tCO₂/d].

_blend_links_in(n, ℒ, links_with_fc, blend)

Return all links in ℒ delivering blend to node n. Throws an error if there are none, or if any of them lacks flow_component variables, since the energy and CO₂ outputs of n would then not be fully constrained.

_is_blend_passthrough(n)

Return true if the blend composition at node n is set by _constraints_blend_passthrough, i.e. copied from upstream without a per-link mixing constraint.

_is_composition_source(n)

Return true if node n sets the composition of its blend outputs itself, so that the composition can be shared by several outgoing links and by composition-dependent BlendLinkData links. These are the pooling nodes, and blend storages in mode :prescribed or :fully_mixed.

_blend_storage_mode(n)

Return the mode of the BlendStorageData of node n, or nothing if n is not a Storage or has no such data.

_constraints_proportion_normalisation(m, n, sub, ℒ_out, 𝒯)

Require the proportions proportion_out[n, t, p] of the sub-resources sub to sum to 1.

When has_flow variables are present (pressure networks with ZeroDrop / Pipeline links), a conditional big-M formulation is used:

Σ_p proportion_out[n, t, p] ≤ 1
Σ_p proportion_out[n, t, p] ≥ has_flow[l, t]   for each outgoing blend link l ∈ ℒ_out

This relaxes the equality when no gas flows (has_flow = 0), allowing the NLP relaxation to set proportion_out = 0. Without it, a node that only receives a single sub-resource (e.g. a pure-H₂ electrolyser) would propagate proportion_out[H2] = 1.0 even at zero flow, violating quality bounds at downstream consumers and causing the NLP relaxation used by MINLP solvers to report spurious local infeasibility. In the MILP (has_flow ∈ {0,1}), the two constraints collapse to the equality whenever an outgoing blend link is active.

_switched_links(m, ℒ)

Return the links whose flow forces has_flow = 1, or nothing if the model has no has_flow variables. These are the links carrying a ResourcePressure or a ResourcePooling{<:ResourcePressure}, for which constraints_flow_capacity adds link_in ≤ M·has_flow. Links with carrier-only blends also get a has_flow binary, but nothing ties it to their flow.

_constraints_passthrough_single_outflow(m, 𝒩, ℒ, 𝒯, 𝒫)

For every passthrough node with several outgoing links carrying a blend, allow at most one of them to be active in each period:

Σ_{l_out} has_flow[l_out, t] ≤ 1

Flow on a pressure link forces has_flow = 1 (constraints_flow_capacity), so the blend leaves the node through a single link, whose component flows equal those entering the node. This keeps the passthrough linear while allowing bidirectional pipelines, where a junction has an outgoing link towards each neighbour but only one carries flow.

_dpn_flow_bound(m, n::DiscretizedPoolingNode, l, 𝒯)

Return a finite upper bound on the blend flow on the outgoing link l of n: the smaller of the installed-capacity bounds of n and l. Used for the flow breakpoints (:pwl) and big-M coefficients (:bigm, :bigm_log).

_maximum_node_capacity(m, n, 𝒯)

Return a finite upper bound on the installed capacity of node n, or nothing.

Same logic as _maximum_link_capacity, using cap_current (investments) and cap_inst.

_pwa_flow_bound(pwa::PWAFunc)

Return a finite upper bound on the flow allowed by the concave PWA pwa over the box p_in, p_out ∈ [0, 1e4] (the big-M bound on link potentials in constraints_balance_pressure) and proportion ∈ [0, 1].

Each plane k gives z ≤ -α_k⋅x - β_k. The maximum of an affine function over a box lies at a corner, and the PWA lies below every plane, so min_k max_corners (-α_k⋅x - β_k) is a valid bound.

_pwa_fit_ok(pwa::PWAFunc, points, z; rtol = 0.1)

Return true if the PWA underestimates the flow z at the fitting points by at most rtol × maximum(z).

_pwa_underestimate(pwa::PWAFunc, points, z)

Return the largest underestimate z_k - pwa(x_k) over the fitting points, relative to maximum(z). A concave PWA gives pwa(x) = min_k (-α_k⋅x - β_k).

_with_seed(f, seed)

Call f() with the global RNG seeded with seed, and restore the caller's RNG state afterwards. Tasks spawned inside f (e.g. parallel Cluster trials) derive their RNG from the seeded state, so the result is reproducible.