Allocation

1 Introduction

This document provides a comprehensive technical reference for the allocation implementation in Ribasim. It bridges the gap between the mathematical formulation described in the concept documentation and the actual code implementation using JuMP.jl.

The allocation algorithm solves linear optimization problems to distribute water among competing demands. This document explains:

1. How mathematical formulations translate to JuMP code
2. Problem Building (placeholder initialization → real value updates)
3. Integration between the physical layer and optimization layer
4. Data structures and indexing patterns
5. Implementation details for each node type

2 Architecture Overview

2.1 Problem Building

The allocation optimization problems are build using a two-phase approach:

2.1.1 Phase 1: Problem Structure Initialization (allocation_init.jl)

During initialization, JuMP variables and constraints are created with placeholder values. These placeholders define the structure and relationships but not the actual values:

# Example from add_basin! function
current_storage = 1000.0 # Placeholder value
max_storage = 5000.0 # Placeholder value
problem[:basin_storage_change] = JuMP.@variable(
    problem,
    -current_storage / scaling.storage ≤
    basin_storage_change[node_id = basin_ids_subnetwork] ≤
    (max_storage - current_storage) / scaling.storage
)

Why placeholders? JuMP requires the problem structure (variables, constraints, objective) to be defined upfront. The actual physical values (water levels, flows, demands) are not yet known during initialization and will change at each allocation timestep.

2.1.2 Phase 2: Real Value Updates (allocation_optim.jl)

Before each optimization, the set_simulation_data! functions update constraints with real values from the physical layer:

# From set_simulation_data! for Basin
for basin_id in basin_ids_subnetwork
    basin_idx = basin_id.idx
    current_storage_basin = current_storage[basin_idx]

    # Update the variable bounds with actual storage values
    storage_change_variable = storage_change[basin_id]
    JuMP.set_lower_bound(
        storage_change_variable,
        -current_storage_basin / scaling.storage
    )
    # ... more updates
end

This pattern allows efficient repeated optimization: the problem structure remains fixed while only the constraint coefficients and bounds are updated.

2.2 Physical Layer Integration

The allocation layer operates on a linearized version of the physical layer. Here’s how physical quantities map to optimization variables:

2.2.1 Basin Levels and Storage

In the physical layer, basin level \(h\) is a nonlinear function of storage \(S\) determined by the basin profile (area-storage relationship). For optimization, we linearize around the current state:

Mathematical formulation: \[h^{n+1} \approx h^n + \frac{1}{A^n}(S^{n+1} - S^n)\]

where \(A^n\) is the basin area at the current timestep.

Code implementation:

# In set_simulation_data! for Basin
# Get linearization point
h_n = basin.current_level[basin_idx]
A_n = current_area[basin_idx]  # Area at current level
S_n = current_storage[basin_idx]

# Define the decision variable for storage change
# ΔS = S^{n+1} - S^n (scaled)
storage_change = problem[:basin_storage_change]

# The level at end of timestep is:
# h^{n+1} = h^n + (1/A^n) * ΔS
# This relationship is embedded in constraint coefficients (see below)

2.2.2 Connector Nodes (Flow as Function of Levels)

Nodes like TabulatedRatingCurve, LinearResistance, and ManningResistance have flow \(Q\) that depends on upstream level \(h_a\) and downstream level \(h_b\):

Mathematical formulation: \[Q^{n+1} \approx Q^n + \frac{\partial Q}{\partial h_a}(h_a - h_a^n) + \frac{\partial Q}{\partial h_b}(h_b - h_b^n)\]

For a basin downstream: \(h_b^{n+1} - h_b^n \approx \frac{1}{A_b}(S_b^{n+1} - S_b^n)\)

Code implementation:

# From linearize_connector_node!
t_after = t + Δt_allocation

# Evaluate flow and derivatives at the linearization point
h_a = get_level(p, inflow_id, t_after)
h_b = get_level(p, outflow_id, t_after)
q0 = flow_function(connector_node, node_id, h_a, h_b, p, t_after)

# Compute partial derivatives with automatic differentiation
∂q∂h_a = forward_diff(
    level_a -> flow_function(connector_node, node_id, level_a, h_b, p, t_after), h_a,
)
∂q∂h_b = forward_diff(
    level_b -> flow_function(connector_node, node_id, h_a, level_b, p, t_after), h_b,
)

# Update constraint: Q = q0 + ∂q/∂h_a * Δh_a + ∂q/∂h_b * Δh_b
# When h_b is a Basin level, substitute the linearized profile:
# Q = q0 + (∂q/∂h_b / A_b) * ΔS_b + ...

The function set_partial_derivative_wrt_level! handles converting level derivatives to storage derivatives:

function set_partial_derivative_wrt_level!(
    allocation_model::AllocationModel,
    node_id::NodeID,
    ∂q∂h::Float64,  # Partial derivative of flow w.r.t. level
    p::Parameters,
    constraint::JuMP.ConstraintRef,
)::Nothing
    (; problem, scaling) = allocation_model
    (; current_area) = p.state_and_time_dependent_cache

    # Convert ∂Q/∂h to ∂Q/∂S using ∂h/∂S = 1/A, and account for the fact that the
    # flow and storage variables are scaled differently
    storage_change = problem[:basin_storage_change][node_id]
    JuMP.set_normalized_coefficient(
        constraint,
        storage_change,
        -∂q∂h * scaling.storage / (scaling.flow * current_area[node_id.idx]),
    )
    return nothing
end

3 Mathematical Formulation to Code Mapping

3.1 Decision Variables

3.1.1 Flow Variables

Mathematical formulation: For each link \((i,j)\) in the network, a flow variable \(Q_{ij}\) represents the volumetric flow rate.

Code implementation:

# From add_flow!
problem[:flow] = JuMP.@variable(problem, flow[link = flow_links_subnetwork])

Key points: - Indexed by link::Tuple{NodeID, NodeID} representing (source, destination) - Declared without bounds; the capacities depend on the allocation timestep and the current physical state, and are therefore set before every optimization by update_flow_variable_bounds! (see Flow capacity bounds) - Scaled by scaling.flow for numerical stability - Stored in sparse array for efficient access: problem[:flow][link]

3.1.2 Basin Storage Change Variables

Mathematical formulation: For each basin \(b\), a storage change variable \(\Delta S_b = S_b^{n+1} - S_b^n\) represents the change in storage over the allocation timestep.

Code implementation:

# From add_basin!
problem[:basin_storage_change] = JuMP.@variable(
    problem,
    -current_storage / scaling.storage ≤
    basin_storage_change[node_id = basin_ids_subnetwork] ≤
    (max_storage - current_storage) / scaling.storage
)

Key points: - Lower bound prevents storage from going negative - Upper bound prevents exceeding maximum basin capacity - Placeholder values replaced before each optimization - Scaled by scaling.storage for numerical conditioning

3.1.3 Allocated Flow Variables (Demand Nodes)

Mathematical formulation: For each demand node \(v\) and priority \(p \in P_v\): \[0 \le F^p_v \le d_v^p\]

Code implementation:

# From add_user_demand!
user_demand_allocated = JuMP.@variable(
    problem,
    0 ≤
        user_demand_allocated[
            node_id = user_demand_ids_subnetwork,
            DemandPriorityIterator(node_id, p_independent),
        ] ≤
        d
)

Key points: - Doubly-indexed by (node_id, demand_priority) - DemandPriorityIterator(node_id, p_independent) yields only the priorities for which this node actually has a demand, so no variables are created for unused combinations - Bounded by the demand value (a placeholder d, updated before optimization by set_demands!) - Similar structures for FlowDemand and for secondary network connections

3.1.4 Error Variables

Mathematical formulation: For UserDemand/FlowDemand: \[E^p_v \ge 0, \quad \overline{E}^p_v \ge 0\]

Code implementation:

# From add_user_demand!
user_demand_error = JuMP.@variable(
    problem,
    0 ≤
        user_demand_error[
            node_id = user_demand_ids_subnetwork,
            DemandPriorityIterator(node_id, p_independent),
            [:first, :second],
        ] ≤
        1
)

Key points: - Triple-indexed: (node_id, demand_priority, objective_ord) - objective_ord = :first: relative error \(E^p_v\) used by the first (demand) objective - objective_ord = :second: fairness error \(\overline{E}^p_v\) used by the second objective - Bounded between 0 and 1 because these are relative errors

3.2 Constraints

3.2.1 Flow Conservation

Mathematical formulation: For conservative nodes (pumps, outlets, resistances), inflow equals outflow: \[Q_{\text{in}} = Q_{\text{out}}\]

Code implementation:

# From add_flow_conservation!
problem[Symbol(constraint_name)] = JuMP.@constraint(
    problem,
    [node_id = node_ids],
    flow[inflow_link[node_id.idx].link] == flow[outflow_link[node_id.idx].link],
    base_name = "flow_conservation_$node_name"
)

Key points: - Simple equality constraint between two flow variables - Applied to Pump, Outlet, LinearResistance, ManningResistance, TabulatedRatingCurve - Ensures mass conservation through the node

3.2.2 Volume Conservation (Basin Water Balance)

Mathematical formulation: \[\frac{dS}{dt} = \sum_{k=1}^{N_l} Q_k + f_{\text{pos}} - f_{\text{neg}}\]

Discretized with backward Euler: \[\frac{S^{n+1} - S^n}{\Delta t} = \sum_{k} Q_k^{n+1} + f^{n+1}_{\text{pos}} - f^{n+1}_{\text{neg}}\]

Or in terms of storage change \(\Delta S = S^{n+1} - S^n\): \[\Delta S = \Delta t \left(\sum_{k} Q_k^{n+1} + f^{n+1}_{\text{pos}} - f^{n+1}_{\text{neg}}\right)\]

Code implementation:

# From add_conservation!
problem[:volume_conservation] = JuMP.@constraint(
    problem,
    [node_id = basin_ids_subnetwork],
    storage_change[node_id] ==
        scaling.flow / scaling.storage *
        (
        f_pos - f_neg * low_storage_factor[node_id] + inflow_sum[node_id] -
            outflow_sum[node_id]
    ),
    base_name = "volume_conservation"
)

Key points: - The low_storage_factor multiplies the implicit negative forcing (evaporation, infiltration) to prevent negative storage - Forcing terms (f_pos, f_neg) represent precipitation, drainage, surface runoff, evaporation and infiltration - The coefficient scaling.flow / scaling.storage converts the scaled flow rates to a scaled storage change over the allocation timestep - Updated each optimization with current forcing values

3.2.3 UserDemand Allocated Sum

Mathematical formulation: The sum of allocated flows over all priorities equals the total inflow, summed over all inflow links of the UserDemand node: \[\sum_{b \in \text{in}(v)} F_{(b, v)} = \sum_{p \in P_v} F^p_v\]

Code implementation:

# From add_user_demand!
problem[:user_demand_allocated_sum_constraint] = JuMP.@constraint(
    problem,
    [node_id = user_demand_ids_subnetwork],
    sum(flow[link_metadata.link] for link_metadata in inflow_links[node_id.idx]) == sum(
        user_demand_allocated[node_id, demand_priority] for
            demand_priority in DemandPriorityIterator(node_id, p_independent)
    );
    base_name = "user_demand_allocated_sum_constraint"
)

Key points: - Ensures only demanded water enters UserDemand nodes - A UserDemand node can have multiple inflow links, so the left-hand side is a sum - Sum on the right is over priorities where the node has a demand - Links the flow variables to the allocated flow variables

3.2.4 Error Constraints

Mathematical formulation: For UserDemand/FlowDemand: \[d_v^p \cdot E_v^p \ge d_v^p - F_v^p\]

Code implementation:

# From add_user_demand!
problem[:user_demand_relative_error_constraint] = JuMP.@constraint(
    problem,
    [
        node_id = user_demand_ids_subnetwork,
        demand_priority = DemandPriorityIterator(node_id, p_independent),
    ],
    d * user_demand_error[node_id, demand_priority, :first] ≥
        d - user_demand_allocated[node_id, demand_priority],
    base_name = "user_demand_relative_error_constraint"
)

Key points: - Multiplication by demand d makes this effectively an absolute error - The error variable is forced to be at least the unmet demand fraction - Updated before each optimization with current demands

For LevelDemand:

Mathematical formulation: \[E^p_{b, \text{lower}} \ge s(h^p_{b, \min}) - (s(h_b^\text{init}) + \Delta S_b)\]

Code implementation:

# From add_level_demand!
problem[:storage_constraint_lower] = JuMP.@constraint(
    problem,
    [
        node_id = basin_ids_subnetwork_with_level_demand,
        demand_priority = DemandPriorityIterator(node_id, p_independent),
    ],
    level_demand_error[node_id, demand_priority, :lower, :first] ≥
        minimum_storage - (starting_storage + storage_change[node_id]),
    base_name = "storage_constraint_lower"
)

Key points: - The error represents the storage deficit below the minimum level - level_demand_error is indexed by (node_id, demand_priority, side, objective_ord) with side ∈ (:lower, :upper) - Separate constraints exist for the upper side, and the demands themselves only appear as constraint constants (they are not decision variables)

3.2.5 Return Flow Constraints

Mathematical formulation: For UserDemand with return factor \(r_i(t)\): \[Q_{\text{out}} = r_i \sum_{b \in \text{in}(v)} Q_{(b, v)}\]

Code implementation:

# From add_user_demand!
problem[:user_demand_return_flow] = JuMP.@constraint(
    problem,
    [node_id = user_demand_ids_subnetwork],
    flow[outflow_link[node_id.idx].link] ==
        return_factor * sum(flow[lm.link] for lm in inflow_links[node_id.idx]),
    base_name = "user_demand_return_flow"
)

Key points: - return_factor is a placeholder, updated before optimization with the value at the end of the allocation timestep - Links the total inflow and the outflow of a UserDemand node - Models consumptive use (when the return factor < 1)

3.2.6 Linearized Connector Node Constraints

Mathematical formulation: For nodes like TabulatedRatingCurve: \[Q^{n+1} \approx Q^n + \frac{\partial Q}{\partial h_a}(h_a^{n+1} - h_a^n) + \frac{\partial Q}{\partial h_b}(h_b^{n+1} - h_b^n)\]

Code implementation:

# From add_linearized_connector_node!
problem[Symbol(constraint_name)] = JuMP.@constraint(
    problem,
    [node_id = node_ids_subnetwork],
    flow[inflow_link[node_id.idx].link] ==
        q0 +  # Flow at linearization point
        (upstream_is_basin ?
            (∂q∂h_upstream / A_upstream) * storage_change[upstream_id] :
            0.0) +
        (downstream_is_basin ?
            (∂q∂h_downstream / A_downstream) * storage_change[downstream_id] :
            0.0),
    base_name = constraint_name
)

Key points: - q0, ∂q∂h_upstream, ∂q∂h_downstream are placeholders - Derivatives divided by area when basin is involved (converts \(\partial Q/\partial h\) to \(\partial Q/\partial S\)) - Updated in set_simulation_data! before each optimization - Handles various combinations (basin-basin, basin-boundary, etc.)

3.3 Objectives

The allocation optimization solves multiple objectives in sequence using lexicographic goal programming. Each objective is optimized; before the next non-empty objective is solved, the preceding result is constrained with an inequality in the direction of improvement.

3.3.1 Primary Objective: Minimize Demand-Weighted Error

Mathematical formulation: For flow demands (UserDemand, FlowDemand) at priority \(p\): \[ \min \sum_{v \in V_p} w_v^p E_v^p, \qquad w_v^p = \frac{d_v^p}{\sum_{u \in V_p} d_u^p}. \]

For level demands at priority \(p\): \[\min \sum_{b; p\in P^\min_b} E^p_{b, \text{lower}} + \sum_{b; p\in P^\max_b} E^p_{b, \text{upper}}\]

Code implementation:

# From add_demand_objectives!
# For flow demands
first_objective_expression = JuMP.AffExpr(0.0)
for node_id in user_demand_ids_subnetwork
    if demand_priority in DemandPriorityIterator(node_id, p_independent)
        demand = get_demand(user_demand, node_id, demand_priority, t)
        JuMP.add_to_expression!(
            first_objective_expression,
            demand * user_demand_error[node_id, demand_priority, :first]
        )
    end
end

# Store both demand goals for later optimization
second_objective_expression = JuMP.AffExpr(0.0)  # Fairness-error expression
push!(
    objectives,
    AllocationObjective(
        AllocationObjectiveType.demand_flow,
        demand_priority,
        demand_priority_idx,
        [first_objective_expression, second_objective_expression],
    )
)

Key points: - Error variables are weighted by their relative demand. This is equivalent to minimizing the total absolute unmet demand because the objective is the original demand-weighted error multiplied by the positive constant \(1 / \sum_{v \in V_p} d_v^p\). - The normalization is applied after all demand coefficients for a priority have been updated, including demands collected from secondary networks. It keeps each first objective and its retained lexicographic constraint numerically well scaled when the flow scale is much larger than the active demands. - Separate objectives for each priority - Objectives are stored and optimized in order by optimize_multi_objective!

3.3.2 Secondary Objective: Minimize Fairness Error

Mathematical formulation: After minimizing the total error, minimize deviations from the average allocation rate.

For flow demands: \[\min \sum_{v,p} \overline{E}_v^p\]

where \(\overline{E}_v^p \ge E_v^p - G^p\) and \(G^p\) is the demand-weighted mean relative error: \[ G^p = \sum_{v \in V_p} w_v^p E_v^p. \]

Code implementation:

# From add_demand_objectives!
# Compute average error (constraint on average error variable)
average_flow_unit_error_constraint =
    @constraint(
        problem,
        [demand_priority = demand_priorities_all; ...],
        average_flow_unit_error[demand_priority] ==
            sum(weight * error[node_id, demand_priority, :first]
                for node_id, weight in demand_weights),
        base_name = "average_flow_unit_error"
    )

# Fairness error: overline_E >= E - G
# (Implicitly defined through constraints)

# Objective sums fairness errors
objective_expression_fairness = sum(
    user_demand_error[node_id, demand_priority, :second]
    for node_id in user_demand_ids_subnetwork if ...
)

Key points: - average_flow_unit_error represents \(G^p\) - Fairness errors (the :second error variables) penalize being worse than average - Optimized after the primary objective is satisfied

3.3.3 LevelDemand Fairness Normalization

LevelDemand first-objective errors represent storage volumes. For a priority \(p\) and side \(s\), the implementation derives the mean current area over participating Basins,

\[ \bar A_{p,s} = \frac{1}{n_{p,s}}\sum_{b \in B_{p,s}} A_b, \]

then defines the mean storage error with \(n_{p,s}\bar E_{p,s} = \sum_b E^p_{b,s}\). The second objective penalizes storage errors after converting them by the relative area scale:

\[ \overline E^p_{b,s} \ge \frac{\bar A_{p,s}}{A_b}E^p_{b,s} - \bar E_{p,s}. \]

The ratios are recomputed for every allocation solve from current Basin areas. This gives the same relative level-error comparison as dividing each storage error by its area, but avoids mixing very large raw-area and very small reciprocal-area coefficients in the LP.

3.3.4 Route Priority Objective

Mathematical formulation: \[\min \sum w_i \cdot Q_i\]

where \(w_i\) is the route priority weight (cost) for node \(i\).

Code implementation:

# From add_route_priority_objective!
objective_expression = JuMP.AffExpr(0.0)
for link in flow_links
    # Get the route priority weight for the source node
    weight = route_priority_weight[link[1]]
    if weight > 0
        JuMP.add_to_expression!(
            objective_expression,
            weight * flow[link]
        )
    end
end

Key points: - Optimized after all demand objectives are satisfied - Only affects routing, not how much is allocated - Higher weight = less preferred route (higher cost)

3.3.5 Low Storage Factor Objective

Mathematical formulation: \[\max \sum_{b \in \text{basins}} \alpha_b\]

where \(\alpha_b\) is the low storage factor for basin \(b\).

Code implementation:

# From add_low_storage_factor_objective!
# Note: Maximizing low_storage_factor is equivalent to minimizing its negative
expression = -sum(
    low_storage_factor[node_id]
    for node_id in basin_ids_subnetwork
)

push!(
    objectives,
    AllocationObjective(
        type = AllocationObjectiveType.low_storage_factor,
        expressions = [expression],
    ),
)

Key points: - Maximizes the low storage factor to allow more outflow from basins - Optimized last to avoid infeasibility from emptying basins - Special priority (-1) ensures it runs after demand objectives

4 Updating Constraints Before Optimization

Before each optimization run, constraints must be updated with current values from the physical layer. This is done through a series of set_simulation_data! functions.

4.1 Basin Updates

function set_simulation_data!(
    allocation_model::AllocationModel,
    basin::Basin,
    p::Parameters,
    t::Float64,
    Δt_allocation::Float64,
    du::CVector,
)::Bool
    # ... setup code ...

    for basin_id in basin_ids_subnetwork
        idx = basin_id.idx
        storage_now = current_storage[idx]
        storage_max = storage_to_level[idx].t[end]

        # Storage may never become negative
        Δstorage = storage_change[basin_id]
        JuMP.set_lower_bound(Δstorage, -storage_now / scaling.storage)

        # The upper bound leaves room for the storage change predicted by the physical
        # layer, so a Basin that is (nearly) full can still receive water
        Δstorage_predicted =
            formulate_dstorage_wrt_time(du, p.p_independent, t, basin_id) * Δt_allocation
        Δstorage_upper = if storage_now > storage_max
            max(2 * Δstorage_predicted, 0.0)
        else
            max(2 * Δstorage_predicted, storage_max - storage_now)
        end
        JuMP.set_upper_bound(Δstorage, Δstorage_upper / scaling.storage)
        # ... rest of constraint updates ...
    end
    return errors
end

Key concepts: - Variable bounds are updated each timestep based on the current Basin storage - The bounds are conservative: they describe the storage range the Basin can reach within this allocation timestep. The same interval is reused to derive the flow capacity bounds of connector nodes - Forcing terms (precipitation, evaporation) are computed and added to constraints - Errors are detected (e.g. storage outside the valid range) and reported

4.2 Connector Node Updates (Linearization)

function linearize_connector_node!(
    allocation_model::AllocationModel,
    connector_node::AbstractParameterNode,
    flow_constraint,
    flow_function::Function,
    p::Parameters,
    t::Float64,
    Δt_allocation,
)
    (; scaling) = allocation_model
    (; inflow_link, outflow_link) = connector_node

    # Evaluate at the end of the allocation timestep (backward Euler)
    t_after = t + Δt_allocation

    for node_id in only(flow_constraint.axes)
        inflow_id = inflow_link[node_id.idx].link[1]
        outflow_id = outflow_link[node_id.idx].link[2]

        # Levels at the linearization point
        h_a = get_level(p, inflow_id, t_after)
        h_b = get_level(p, outflow_id, t_after)

        # Flow at the linearization point becomes the constraint right-hand side
        constraint = flow_constraint[node_id]
        q0 = flow_function(connector_node, node_id, h_a, h_b, p, t_after)
        JuMP.set_normalized_rhs(constraint, q0 / scaling.flow)

        # Only linearize if the level comes from a Basin, since only then it is a
        # decision variable of the problem
        if inflow_id.type == NodeType.Basin
            ∂q∂h_a = forward_diff(
                level_a -> flow_function(connector_node, node_id, level_a, h_b, p, t_after),
                h_a,
            )
            set_partial_derivative_wrt_level!(
                allocation_model, inflow_id, ∂q∂h_a, p, constraint,
            )
        end
        # ... same for outflow_id ...
    end
end

Key concepts: - Derivatives are computed with automatic differentiation (forward_diff) of the same flow functions the physical layer uses, so the allocation layer cannot drift from the physics - Only Basin levels are decision variables; LevelBoundary and Terminal levels enter through the constant term q0 - Backward Euler: everything is evaluated at the end of the timestep (\(t + \Delta t\)). For Basins get_level returns the level at the start of the timestep, which is exactly the point we want to linearize around

4.3 Demand Updates

function set_demands!(
    allocation_model::AllocationModel,
    node::Union{UserDemand, FlowDemand},
    # ... other parameters
)::Nothing
    for node_id in demand_node_ids_subnetwork
        for demand_priority in DemandPriorityIterator(node_id, p_independent)
            # Get current demand from the physical layer, in scaled units
            demand = get_demand(node, node_id, demand_priority, t) / scaling.flow

            # Update error constraints with the current demand
            for objective_ord in (:first, :second)
                error_constraint =
                    node_relative_error_constraint[node_id, demand_priority, objective_ord]
                JuMP.set_normalized_coefficient(
                    error_constraint,
                    node_error[node_id, demand_priority, objective_ord],
                    demand,
                )
                JuMP.set_normalized_rhs(error_constraint, demand)
            end

            # Update upper bound on allocated variable
            JuMP.set_upper_bound(node_allocated[node_id, demand_priority], demand)
        end
    end
    return nothing
end

Key concepts: - Demands are time-varying and interpolated from input tables - Constraint coefficients are updated to reflect current demands - The demands are stored in the problem already divided by scaling.flow, so the bounds of the allocated variables are directly usable as scaled flow bounds elsewhere (see Flow capacity bounds) - Upper bounds on allocation variables ensure allocated ≤ demanded

4.4 Flow capacity bounds

Flow variables are declared without bounds, because their capacities depend on the current state of the physical layer and on the allocation timestep. They are set before every optimization, in this order:

1. update_flow_variable_bounds! sets the physical capacity of every non-fixed flow variable;
2. update_user_demand_flow_bounds! tightens the bounds of the links around UserDemand nodes to their total demand;
3. update_flow_demand_variable_bounds! derives the bounds of the FlowDemand helper variables from the now finalized flow bounds.

All of these use the helper set_variable_bounds!, where an infinite bound means that no JuMP bound is set at all:

function set_variable_bounds!(
    variable::JuMP.VariableRef,
    lower::Float64,
    upper::Float64,
)::Nothing
    JuMP.has_lower_bound(variable) && JuMP.delete_lower_bound(variable)
    JuMP.has_upper_bound(variable) && JuMP.delete_upper_bound(variable)
    isfinite(lower) && JuMP.set_lower_bound(variable, lower)
    isfinite(upper) && JuMP.set_upper_bound(variable, upper)
    return nothing
end

4.4.1 Physical capacities

Each flow link is bounded by both of the nodes it connects, so the bounds are the intersection of the intervals returned by connector_flow_capacity_bounds for either node:

function update_flow_variable_bounds!(
    allocation_model::AllocationModel,
    integrator::DEIntegrator,
    Δt_allocation::Float64,
)::Nothing
    (; p, t) = integrator
    (; problem, scaling, flow_links_subnetwork) = allocation_model
    flow = problem[:flow]

    for flow_link in flow_links_subnetwork
        flow_var = flow[flow_link]
        # Flows out of FlowBoundary nodes are fixed to the boundary flow
        JuMP.is_fixed(flow_var) && continue

        lower, upper = -Inf, Inf
        for node_id in flow_link
            node_lower, node_upper = connector_flow_capacity_bounds(
                allocation_model, node_id, p, t, Δt_allocation,
            )
            lower = max(lower, node_lower)
            upper = min(upper, node_upper)
        end
        lower ≤ upper ||
            error("Empty flow capacity interval [$lower, $upper] for $flow_link.")
        set_variable_bounds!(flow_var, lower / scaling.flow, upper / scaling.flow)
    end
    return nothing
end

connector_flow_capacity_bounds dispatches on the node type:

• Pump and Outlet: the tightest min_flow_rate and max_flow_rate over \([t, t + \Delta t]\), since allocation assigns a single flow rate for the whole timestep;
• LinearResistance, ManningResistance and non-allocation-controlled TabulatedRatingCurve: linearized_flow_bounds, see below;
• allocation controlled TabulatedRatingCurve: \([0, Q(h_{a,\max})]\) where \(h_{a,\max}\) is the highest reachable upstream level;
• UserDemand: \([0, \infty)\), tightened afterwards by update_user_demand_flow_bounds!;
• anything else: \((-\infty, \infty)\), i.e. no bound from this node.

4.4.2 Bounds from the linearized flow relation

For nodes whose flow follows from a linearized flow-level relation, the tightest possible bounds follow from evaluating that same linearization over the reachable level changes:

function linearized_flow_bounds(
    allocation_model, connector_node, flow_function, node_id, p, t, Δt_allocation,
)::Tuple{Float64, Float64}
    t_after = t + Δt_allocation
    inflow_id = connector_node.inflow_link[node_id.idx].link[1]
    outflow_id = connector_node.outflow_link[node_id.idx].link[2]

    h_a = get_level(p, inflow_id, t_after)
    h_b = get_level(p, outflow_id, t_after)
    q0 = flow_function(connector_node, node_id, h_a, h_b, p, t_after)

    lower, upper = q0, q0
    ∂q∂h_a = forward_diff(
        level_a -> flow_function(connector_node, node_id, level_a, h_b, p, t_after), h_a,
    )
    ∂q∂h_b = forward_diff(
        level_b -> flow_function(connector_node, node_id, h_a, level_b, p, t_after), h_b,
    )
    for (∂q∂h, level_id) in ((∂q∂h_a, inflow_id), (∂q∂h_b, outflow_id))
        iszero(∂q∂h) && continue
        Δh_min, Δh_max = allocation_level_change_bounds(allocation_model, p, level_id)
        Δq_min, Δq_max = minmax(∂q∂h * Δh_min, ∂q∂h * Δh_max)
        lower += Δq_min
        upper += Δq_max
    end
    return (lower, upper)
end

The reachable level change follows from the bounds on the storage change decision variable, converted with the same linearized profile \(\Delta h = \Delta S / A^n\) the constraints use:

function allocation_level_change_bounds(
    allocation_model::AllocationModel, p::Parameters, node_id::NodeID,
)::Tuple{Float64, Float64}
    (; problem, scaling, node_ids_in_subnetwork) = allocation_model

    # Levels of all other node types are data rather than decision variables,
    # so they cannot change within one optimization
    node_id.type == NodeType.Basin || return (0.0, 0.0)

    # `flow_links_subnetwork` can contain links to Basins outside this subnetwork
    (node_id ∈ node_ids_in_subnetwork.basin_ids_subnetwork) || return (-Inf, Inf)

    storage_change = problem[:basin_storage_change][node_id]
    level_per_storage =
        scaling.storage / p.state_and_time_dependent_cache.current_area[node_id.idx]
    return (
        JuMP.lower_bound(storage_change) * level_per_storage,
        JuMP.upper_bound(storage_change) * level_per_storage,
    )
end

Because these bounds are derived from the same linearization that appears in the flow constraint, they can never conflict with the constraint and therefore never make the problem infeasible.

Important

For allocation controlled Pump and Outlet nodes there is no fallback bound, so a finite max_flow_rate is required. This is checked by valid_allocation_flow_capacity in core/src/validation.jl, both for the current parameters and for all DiscreteControl states, so a control state change cannot make the problem unbounded later in the simulation.

4.4.3 Demand-based bounds

update_user_demand_flow_bounds! bounds all inflow links and the outflow link of a UserDemand node between zero and its total demand over all demand priorities. For the outflow the return factor is conservatively assumed to be 1. Since the demand values stored in the problem are already divided by scaling.flow, the upper bounds of the allocated variables can be used directly.

update_flow_demand_variable_bounds! then takes the (now final) bounds of the flow variable of the node with a flow demand and transfers them to the helper variables: the allocated variable of the earliest demand priority may go down to \(\min(0, Q^\min)\) to represent negative flow, and the ‘demand priority 0’ surplus variable is bounded above by \(\max(0, Q^\max)\). If the flow variable has no bound on a side, the corresponding helper variable bound is deleted as well.

5 The Optimization Loop

5.1 Orchestration: update_allocation!

update_allocation! is the entry point that is called every allocation timestep. It first refreshes the physical state (water_balance!) and the scaling factors, and then per subnetwork performs the update sequence:

update_control_states!(allocation_model, p_independent)
set_simulation_data!(allocation_model, integrator, Δt)
reset_demand_coefficients(allocation_model)
set_demands!(allocation_model, integrator, Δt)
update_flow_variable_bounds!(allocation_model, integrator, Δt)
update_user_demand_flow_bounds!(allocation_model, p_independent)
update_flow_demand_variable_bounds!(allocation_model, p_independent)
normalize_flow_demand_objectives!(allocation_model)
warm_start!(allocation_model, integrator, Δt)

The order matters: set_demands! must run before the flow bounds are set, because update_user_demand_flow_bounds! reads the demands from the bounds of the allocated variables, and update_flow_demand_variable_bounds! in turn reads the finalized flow variable bounds.

Secondary networks are updated first. If a primary network is present, its update sequence is interrupted after reset_demand_coefficients to run the demand collection for every secondary network.

Afterwards, for every subnetwork in turn the problem is optimized (optimize_multi_objective!), the results are parsed and saved, the primary network allocation is distributed over the secondary networks (allocate_flows_to_subnetwork), and the results are applied to the physical layer with apply_control_from_allocation!.

5.2 Lexicographic optimization

The main optimization loop in optimize_multi_objective! solves objectives in sequence:

function optimize_multi_objective!(
    model::AllocationModel,
    config::Config,
    t::Number,
    # ...
)::Nothing
    (; problem, objectives, temporary_constraints) = model

    latest_expression = JuMP.AffExpr()
    latest_bound = 0.0
    for objective in objectives
        for (expression_idx, expression) in enumerate(objective.expressions)
            iszero(expression) && continue
            if !iszero(latest_expression) && !latest_expression_is_constrained
                push!(
                    temporary_constraints,
                    JuMP.@constraint(problem, latest_expression <= latest_bound),
                )
            end
            JuMP.@objective(problem, Min, expression)
            JuMP.optimize!(problem)
            latest_constraint =
                isempty(temporary_constraints) ? nothing : last(temporary_constraints)
            parse_termination_status(
                model, objective, expression, latest_constraint, config, t,
            )
            if objective.retain_expressions[expression_idx]
                latest_bound = JuMP.objective_value(problem)
                latest_expression = expression
                latest_expression_is_constrained = false
            end
        end
    end

    # Clean up temporary constraints after optimization
    delete_temporary_constraints!(allocation_model)

    return nothing
end

Key concepts: - Lexicographic optimization: Each objective is optimized in order - Before optimizing objective \(i + 1\), add constraint: \(\text{obj}_i \le \text{optimal}_i\) - This ensures later objectives don’t degrade earlier ones - An expression can be non-retained. It is optimized as a deterministic tie-break, but is not constrained for the next expression. - Temporary constraints are removed after all objectives are solved - The termination status is checked after every objective optimization; an infeasible or other non-optimal solve stops allocation immediately

6 Scaling for Numerical Stability

The optimization uses scaling factors to improve numerical conditioning:

@kwdef mutable struct ScalingFactors
    flow::Float64 = 1.0e3            # Typical flow rate (m³/s)
    storage::Float64 = 1.0e6         # Typical storage value (m³)
    mean_half_storage::Float64 = 0.0 # Initialization-time reference storage (m³)
end

function update_scaling!(p::Parameters, Δt::Float64)
    (; p_independent, state_and_time_dependent_cache) = p
    (; allocation_models) = p_independent.allocation
    (; current_storage) = state_and_time_dependent_cache
    for allocation_model in allocation_models
        (; node_ids_in_subnetwork, scaling) = allocation_model
        (; basin_ids_subnetwork) = node_ids_in_subnetwork

        storage_sum = sum(current_storage[id.idx] for id in basin_ids_subnetwork)
        scaling.storage =
            (scaling.mean_half_storage + storage_sum / length(basin_ids_subnetwork)) / 2
        scaling.flow = scaling.storage / Δt
    end
    return
end

mean_half_storage is computed once when the allocation problem is built, as the mean over the Basins in the subnetwork of half their maximum storage.

Key concepts: - Storage variables are divided by scaling.storage - Flow variables are divided by scaling.flow - The storage scale is updated for each allocation solve from the current mean storage and the initialization-time mean half-storage - The flow scale equals scaling.storage / Δt, so that a typical flow over an allocation timestep amounts to a typical storage change - Per-priority flow-demand objectives are normalized independently, because a global flow scale can still make a small active demand numerically insignificant - Improves numerical stability and solver performance - Must scale/unscale when communicating with physical layer

6.1 Example: Scaled Volume Conservation

Unscaled formulation: \[\Delta S \text{ [m³]} = \Delta t \text{ [s]} \cdot \left( Q_{\text{in}} - Q_{\text{out}} \text{ [m³/s]} \right)\]

Scaled formulation: \[\tilde{\Delta S} \cdot S_{\text{scale}} = \Delta t \cdot \left( \tilde{Q}_{\text{in}} \cdot Q_{\text{scale}} - \tilde{Q}_{\text{out}} \cdot Q_{\text{scale}} \right)\]

where \(\tilde{\Delta S} = \Delta S / S_{\text{scale}}\) and \(\tilde{Q} = Q / Q_{\text{scale}}\).

Dividing through by \(S_{\text{scale}}\): \[\tilde{\Delta S} = \frac{\Delta t \cdot Q_{\text{scale}}}{S_{\text{scale}}} \left( \tilde{Q}_{\text{in}} - \tilde{Q}_{\text{out}} \right)\]

In the code:

storage_change[node_id] ==
    scaling.flow / scaling.storage * (
        f_pos - f_neg * low_storage_factor[node_id] +
        inflow_sum[node_id] - outflow_sum[node_id]
    )

Since scaling.flow = scaling.storage / Δt, this single coefficient carries both the unit conversion from flow rate to volume over the allocation timestep and the ratio of the two scaling factors. It is updated together with the rest of the constraint in set_simulation_data!.

7 Data Structures and Indexing

7.1 NodeID

Nodes are identified by the NodeID struct:

struct NodeID
    type::NodeType.T  # e.g., NodeType.Basin, NodeType.UserDemand
    idx::Int32        # Index into the type-specific data arrays
    value::Int32      # User-facing ID from input file
end

Key concepts: - idx is used to index into type-specific arrays (e.g., basin.storage[basin_id.idx]) - value is the ID from the input file (for error messages, output) - type distinguishes different node types

7.2 Links

Links are represented as tuples of NodeIDs:

link::Tuple{NodeID, NodeID}  # (source, destination)

Key concepts: - Used to index flow variables: problem[:flow][link] - Direction matters: (A, B) ≠ (B, A) - Stored in graph edge metadata

7.3 JuMP Sparse Arrays

JuMP uses SparseAxisArray for multi-dimensional variables:

# Example: doubly-indexed variable
user_demand_allocated[node_id, demand_priority]

# Example: triply-indexed variable
user_demand_error[node_id, demand_priority, objective_ord]  # objective_ord ∈ (:first, :second)

Key concepts: - Only creates variables for valid index combinations - Filtered by iterating over DemandPriorityIterator(node_id, p_independent), which yields only the demand priorities the node actually has - Efficient for sparse problems (not all nodes have all priorities)

7.4 AllocationModel Struct

The main data structure:

@kwdef mutable struct AllocationModel
    subnetwork_id::Int32
    node_ids_in_subnetwork::NodeIDsInSubnetwork
    problem::JuMP.Model  # The JuMP optimization problem
    # Time since the last time results were written
    Δt_since_last_record::Float64 = 0.0
    has_demand_priority::Vector{Bool}
    # Optimization objectives
    objectives::Vector{AllocationObjective} = AllocationObjective[]
    # Forcing volumes (precipitation, drainage, evaporation, infiltration)
    explicit_positive_forcing_volume::OrderedDict{NodeID, Float64} = OrderedDict()
    implicit_negative_forcing_volume::OrderedDict{NodeID, Float64} = OrderedDict()
    # Cumulative tracking for output
    cumulative_supplied_volume::OrderedDict{Tuple{NodeID, NodeID}, Float64} = OrderedDict()
    sources::OrderedDict{Int32, NodeID} = OrderedDict()
    secondary_network_demand::OrderedDict{Tuple{NodeID, NodeID}, Vector{Float64}} =
        OrderedDict()
    # All links with at least one node in this subnetwork
    flow_links_subnetwork::Vector{Tuple{NodeID, NodeID}} = Tuple{NodeID, NodeID}[]
    scaling::ScalingFactors = ScalingFactors()
    # Current LevelDemand fairness normalization by demand priority
    level_demand_area_sum::Dict{Int32, Float64} = Dict()
    level_demand_count::Dict{Int32, Int} = Dict()
    level_demand_area_scale::Dict{Int32, Float64} = Dict()
    # Temporary constraints for lexicographic optimization
    temporary_constraints::Vector{JuMP.ConstraintRef} = JuMP.ConstraintRef[]
end

Note

The allocation timestep Δt_allocation is not a field of AllocationModel. Because the timestep is adaptive, it is passed as a function argument to all functions that need it. Δt_since_last_record is unrelated: it accumulates the time since results were last written.

8 Secondary Networks and Primary-Secondary Connections

8.1 Primary Network

The primary network (subnetwork_id = 1) represents the main water system. It:

• Contains its own demand nodes
• Connects to secondary networks via Pump or Outlet nodes
• Treats each secondary network as a single demand node

8.2 Secondary Networks

Secondary networks (subnetwork_id > 1):

• Can only connect to the primary network
• Cannot connect to each other directly
• Have their own internal allocation optimization

8.3 The Two-Stage Process

8.3.1 Stage 1: Demand Collection

For each secondary network:

1. Allow inflow from the primary network up to the physical capacity of the connecting links
2. Solve the demand and fairness objectives for each priority
3. Minimize total primary-to-secondary inlet flow as a non-retained tie-break
4. Use that minimum inlet flow as the secondary network’s request for the priority

The tie-break removes solver-dependent excess flow from demand collection. It is non-retained because later priorities may legitimately need more inlet flow.

function preprocess_demand_collection!(
    allocation_model::AllocationModel,
    p_independent::ParametersIndependent,
)::Nothing
    (; problem, subnetwork_id) = allocation_model
    @assert !is_primary_network(subnetwork_id)
    flow = problem[:flow]

    # The capacity bounds of the primary network connections were set from the current
    # allocation data, so demands are collected up to the physical capacity. The primary
    # network allocation will subsequently restrict these links to its result.
    for link in p_independent.allocation.primary_network_connections[subnetwork_id]
        flow_variable = flow[link]
        lower_bound =
            JuMP.has_lower_bound(flow_variable) ? JuMP.lower_bound(flow_variable) : 0.0
        JuMP.set_lower_bound(flow_variable, max(0.0, lower_bound))
    end

    return nothing
end

Note that only the lower bound needs adjusting here: the upper bound already is the physical capacity of the link, set by update_flow_variable_bounds!. Since the connecting node is always an allocation controlled Pump or Outlet, that capacity is finite.

8.3.2 Stage 2: Allocation

In primary network:

1. Solve allocation considering:
   • Primary network’s own demands
   • Aggregated demands from secondary networks
2. Allocated amounts to secondary networks determine their inflow limits

In each secondary network:

1. Set inflow constraint to the allocated amount from primary network
2. Solve allocation to distribute water to internal demands

function allocate_flows_to_subnetwork(
    allocation_models::Vector{AllocationModel},
    primary_network_connections,
)::Nothing
    primary_network = get_primary_network(allocation_models)
    primary_problem = primary_network.problem

    for secondary_network in get_secondary_networks(allocation_models)
        # Get allocated flow from primary network
        for link in primary_network_connections[secondary_network.subnetwork_id]
            allocated_flow = JuMP.value(primary_problem[:flow][link])

            # Set as capacity for secondary network
            secondary_problem = secondary_network.problem
            JuMP.set_upper_bound(
                secondary_problem[:flow][link],
                allocated_flow
            )
        end
    end
end

9 Warm Start

To improve solver performance, the optimization can be warm-started with flow rates from the physical layer:

function warm_start!(
    allocation_model::AllocationModel,
    integrator::DEIntegrator,
    Δt_allocation::Float64,
)::Nothing
    (; p, t) = integrator
    (; problem, scaling, node_ids_in_subnetwork) = allocation_model
    (; basin_ids_subnetwork) = node_ids_in_subnetwork
    flow = problem[:flow]
    storage_change = problem[:basin_storage_change]
    du = get_du(integrator)
    (; link_to_state_idx) = p.p_independent

    # Extrapolate the current instantaneous flow rates from the physical layer
    for link in only(flow.axes)
        state_index = get_state_index(getaxes(du), link_to_state_idx, link)
        if !isnothing(state_index)
            JuMP.set_start_value(flow[link], du[state_index] / scaling.flow)
        end
    end

    # Extrapolate the current instantaneous storage rates from the physical layer
    for node_id in basin_ids_subnetwork
        JuMP.set_start_value(
            storage_change[node_id],
            formulate_dstorage_wrt_time(du, p.p_independent, t, node_id) * Δt_allocation /
                scaling.storage,
        )
    end

    return nothing
end

Key concepts: - Starting values guide the solver to a good initial solution - Can significantly reduce solver iterations - Based on extrapolating current physical state forward - Not every flow link corresponds to a state of the physical layer, hence the isnothing check

10 Output and Communication with Physical Layer

10.1 Parsing Results

After optimization, results must be extracted and stored. parse_allocations! has a dispatcher method plus one method for UserDemand and FlowDemand nodes and one for LevelDemand nodes:

function parse_allocations!(
    integrator::DEIntegrator,
    allocation_model::AllocationModel,
    Δt_allocation,
)::Nothing
    # Dispatches to the UserDemand, FlowDemand and LevelDemand methods
end

# For UserDemand and FlowDemand
for node_id in demand_node_ids_subnetwork
    for demand_priority in DemandPriorityIterator(node_id, p_independent)
        # Get allocated amount from the optimization result and unscale it
        allocated = JuMP.value(node_allocated[node_id, demand_priority]) * scaling.flow

        # Store on the node for use by the physical layer
        node.allocated[node_id.idx, demand_priority_idx] = allocated

        # Store for output in the demand results table
        push!(record_demand, DemandRecordDatum(...))
    end
end

For LevelDemand nodes there is no allocated flow variable; instead the volume supplied to the Basin over the allocation timestep is derived from the optimized storage change and recorded.

10.2 Applying Allocation Results

Allocated amounts must be communicated to the physical layer.

10.2.1 UserDemand Nodes

The allocated flow rates per demand priority are written to user_demand.allocated. In the physical layer, the abstraction of a UserDemand node is limited by the sum of these allocated flow rates rather than by its demand, so that the node never takes more than it was allocated.

10.2.2 Pump, Outlet and TabulatedRatingCurve Control

When a Pump, Outlet or TabulatedRatingCurve has allocation_controlled set to true:

function apply_control_from_allocation!(
    node::Union{Pump, Outlet, TabulatedRatingCurve},
    allocation_model::AllocationModel,
    integrator::DEIntegrator,
)::Nothing
    (; problem, scaling) = allocation_model
    flow = problem[:flow]

    for node_id in controlled_node_ids
        # Get optimized flow from allocation
        link = (inflow_id(node_id), outflow_id(node_id))
        optimized_flow = JuMP.value(flow[link]) * scaling.flow

        # Set as target flow rate in physical layer
        node.flow_rate[node_id.idx] = optimized_flow
    end
end

The physical layer still applies the usual reduction factors (low storage, level difference, maximum downstream level) to this flow rate, so the realized flow can be lower than the allocated flow rate.

11 Handling Infeasibility

When the optimization problem is infeasible, diagnostic tools help identify the cause:

11.1 Infeasibility Analysis

function analyze_infeasibility(
    allocation_model::AllocationModel,
    t::Float64,
    config::Config,
)::JuMP.TerminationStatusCode
    (; problem) = allocation_model

    # Find Irreducible Inconsistent Subsystem (IIS)
    data_infeasibility = MathOptAnalyzer.analyze(
        MathOptAnalyzer.Infeasibility.Analyzer(),
        problem;
        optimizer = get_optimizer(),
    )

    # Extract conflicting constraints
    violated_constraints = [...]

    # Try relaxing constraints to identify issues
    constraint_to_penalty = Dict(
        constraint => (isempty(JuMP.name(constraint)) ? 1.0 : 0.5)
        for constraint in violated_constraints
    )
    constraint_to_slack = JuMP.relax_with_penalty!(problem, constraint_to_penalty)
    JuMP.optimize!(problem)

    # Report which constraints are in conflict
    for irreducible_infeasible_subset in data_infeasibility.iis
        constraint_violations = [...]
        @error "Set of incompatible constraints found" constraint_violations
    end

    return JuMP.INFEASIBLE
end

Key concepts: - IIS = minimal set of constraints that cannot be satisfied simultaneously - Constraint relaxation helps identify which constraints are problematic - Named constraints are more informative for debugging

11.2 Common Infeasibility Causes

1. Basin overflows: In allocation all variables must have upper bounds. The upper bound of the basin level is defined by the user’s input. If more water flows into the basin than is possible based on the user defined basin dimensions, the model is infeasible.

12 Performance Considerations

12.1 Problem Size

• Variables: O(links + basins + demand_nodes × priorities)
• Constraints: O(links + basins + demand_nodes × priorities)
• Objectives: O(priorities)

For a network with 100 nodes, 200 links, 20 demand nodes, and 4 priorities: - ~300 variables - ~400 constraints - 4-8 objectives (depending on objectives configured)

12.2 Solver Performance

The HiGHS solver is configured for allocation problems:

function get_optimizer()
    JuMP.optimizer_with_attributes(
        HiGHS.Optimizer,
        "log_to_console" => false,
        "time_limit" => 60.0,
        "random_seed" => 0,
        "small_matrix_value" => 1e-12,  # Numerical threshold
    )
end

Tips for performance: - Use scaling to keep values O(1) - Minimize number of priorities (each adds objectives) - Warm start with physical layer flows - Simplify network topology where possible

12.3 Allocation Timestep Selection

By default the allocation timestep \(\Delta t_{\text{allocation}}\) is determined adaptively from bounds on the linearization errors, see Adaptive timestepping and compute_adaptive_allocation_Δt. It can be overridden with a fixed value using the allocation.dt configuration setting. The trade-off is:

• Too small: Frequent optimization overhead, linearization always valid
• Too large: Linearization may be inaccurate, less responsive to changes

13 Debugging Tips

13.1 Inspecting the Problem

# Write problem to file for manual inspection
write_problem_to_file(problem, config)

# Check variable values
println(JuMP.value(flow[link]))
println(JuMP.value(storage_change[basin_id]))

# Check constraint values
println(JuMP.normalized_rhs(constraint))
println(JuMP.normalized_coefficient(constraint, variable))

13.2 Common Issues

13.2.1 Variables at Bounds

function get_bounds_hit(variable::JuMP.VariableRef)::Tuple{Bool, Bool}
    hit_lower_bound = if JuMP.has_lower_bound(variable)
        JuMP.value(variable) ≤ JuMP.lower_bound(variable)
    else
        false
    end

    hit_upper_bound = if JuMP.has_upper_bound(variable)
        JuMP.value(variable) ≥ JuMP.upper_bound(variable)
    else
        false
    end

    return hit_lower_bound, hit_upper_bound
end

If many variables are at bounds, the problem may be under-constrained or infeasible.

13.2.2 Numerical Scaling Issues

function analyze_scaling(
    allocation_model::AllocationModel,
    t::Float64,
    config::Config,
)::Nothing
    data_numerical = MathOptAnalyzer.analyze(
        MathOptAnalyzer.Numerical.Analyzer(),
        problem;
        threshold_small = 1e-12,
        threshold_large = 1e6,
    )

    # Check for poorly scaled coefficients
    for data in data_numerical.matrix_small
        @error "Too small coefficient" data.coefficient
    end

    for data in data_numerical.matrix_large
        @error "Too large coefficient" data.coefficient
    end
end

14 Summary

This technical reference provides the link between mathematical formulation and code implementation:

1. Two-phase approach: Structure defined at initialization with placeholders, values updated before each optimization
2. Physical layer integration: Linearization of basin profiles and connector nodes around current state
3. Decision variables: Flows, storage changes, allocated amounts, errors
4. Flow capacity bounds: Derived before each optimization from the physical node properties and the reachable Basin levels
5. Constraints: Flow conservation, volume conservation, demand relationships, linearized physics
6. Objectives: Lexicographic goal programming for demands, fairness, route priorities
7. Scaling: Improves numerical stability by keeping values O(1)
8. Data structures: NodeID-based indexing, sparse JuMP arrays, link tuples
9. Primary-secondary networks: Two-stage allocation process

For user-facing documentation, see Allocation Concept. For implementation code, see core/src/allocation_init.jl, core/src/allocation_optim.jl, core/src/allocation_util.jl, and core/src/validation.jl.