Flow communication between primary and secondary networks

[SouthEndMusic]

Fixes #2327
Fixes #2230
Fixes #2383

[SouthEndMusic]

I stumbled upon a problem when working on this. We determine the demand(s) of a secondary network to the primary network by optimizing within the secondary network, minimizing how much flow is taken from the main network and interpreting the flow into the secondary network from the primary network as the demand of the secondary network.

The problem is that with the current formulation it is not possible to split up that demand into the different priorities, the optimization only yields the total flow into the secondary network from the primary network. We can't use a global (multi commodity) water balance type computation for this either, because we allow multiple inlets from the primary network into secondary networks.

We can add decision variables per demand priority at each secondary network inlet, but then we need to have a way to constrain those. I cannot quickly come up with a way to do that without adding a new objective per demand priority.

@gijsber @jarsarasty @simulutions

[SouthEndMusic]

I'm having trouble coming up with something that works within the lexicographic multi-objective optimization. Say we add a flow variable per demand priority for each primary network to secondary network inlet. A naive solution would be to, after optimizing for a certain demand priority when collecting demands in a secondary network, minimize the sum of the inflows variables from the primary network corresponding to this demand priority. Then however the solver can 'abuse' the inlet variables for later priorities to allocate for this priority. I don't see how this can be avoided other than by manually adjusting constraints/bounds in between optimizations. @jarsarasty do you have a clever trick up your sleeve for this?

[jarsarasty]

Hi @SouthEndMusic

Looking at the challenge of splitting secondary network demands by priority and making this solution information available to the primary network, I see three possible approaches:

Alternative 1: Post-Processing

Core Idea:

Extract priority-based demands from secondary network solutions after optimization but before primary network optimization.

Implementation Outline:

1. Optimize secondary networks with current variable structure (no changes)
2. Extract priority demands from error variable solutions:
   • priority_1_allocated = total_priority_1_demand - JuMP.value(error_priority_1_first)
   • priority_2_allocated = total_priority_2_demand - JuMP.value(error_priority_2_first)
3. Communicate to primary network: Pass [priority_1: X, priority_2: Y] demands
4. Primary network optimizes with priority-specific demand information

Advantages:

No additional optimization variables, minimal code changes, maintains solver performance

Limitations or disasvantages:

The demand priorities may not be acuratelly reconstructed in cases of multiple inlets per secondary network. However, multiple-inlet splitting not be immediately relevant as we are currently not considering this case, as far as I remember.

Alternative 2: Priority-indexed flow variables

Core Idea:

Add priority-indexed flow variables for primary→secondary connections, letting MOA.jl naturally discover optimal priority distribution.

Implementation Outline:

1. Add priority-indexed inlet variables (the missing link):
   secondary_inlet_flow_by_priority[inlet, priority] = @variable(problem, ...)
2. Keep existing demand variables unchanged (they're already priority-indexed):
3. Add flow conservation constraint:

  @constraint(problem,
      [inlet = connecting_links],
      flow[inlet] == sum(secondary_inlet_flow_by_priority[inlet, p] for p in priorities)
  )

4. Link priority inlet flows to existing demand variables:

@constraint(problem,
    [priority = priorities],
    sum(secondary_inlet_flow_by_priority[inlet, priority] for inlet in
secondary_inlets) ==
    sum(user_demand_allocated[node, priority] for node in secondary_demand_nodes)
)

5. Objectives remain completely unchanged (they already work with existing priority-indexed variables!)

Key Insight: Most of the priority structure already exists in Ribasim - we're just
"filling in the missing link" between total inlet flows and priority-indexed demands.

Advantages:

Mathematically optimal, handles complex topologies naturally, and works with multiple inlets

Limitations or disasvantages:

The number of variables would increase in around 20% to 60% in some cases, and this could impact the performance. However, this impact would not be very drastic as the current models are linear.

Alternative 3: Dual Decomposition

Core Idea:

Use Lagrangian relaxation to coordinate between primary and secondary networks through price signals.

Implementation Outline:

1. Initialize dual multipliers (λ) representing "water prices" for each priority
2. Iterate between:
   - Subproblems: Optimize each secondary network with priority-specific prices
   - Master problem: Optimize primary network allocation based on secondary demands
   - Dual updates: Adjust λ based on capacity constraint violations
3. Convergence: Stop when dual multipliers stabilize

Advantages:

Theoretically elegant, proven convergence under convexity assumption (current situation), handles complex coupling naturally

Disadvantages:

Most complex to implement, requires dual update tuning, may need many iterations

Recommendation

Given that current test models use single inlet per secondary network, I'd suggest starting with Alternative 1 (Post-Processing) for immediate implementation, since:

• The multiple-inlet complexity doesn't currently exist in practice
• Simpler to implement and test
• Can be extended to Alternative 2 later if multiple inlets become needed

Alternative 2 would be the best long-term solution for mathematical rigor and extensibility. However, we should evaluate its impact on performance

Alternative 3 is more of academic interest at this moment, unless we encounter some limitations with the other approaches.

Shall we plan a meeting to discuss these approaches?

[SouthEndMusic]

Hi @jarsarasty, thank you very much for your elaborate response. Here's my response per alternative:

Alternative 1: Post-Processing

We do want to support multiple inlets from the primary network to the secondary network (I'm pretty sure this is a requirement for the Dutch schematization). Another disadvantage of this approach is that it does not take the topology and constraints between the demands and secondary network inlets into account.

Alternative 2:

Here I don't fully understand 4. It looks like generally you want to equate the total allocated flow for a certain priority in a secondary network to the inflow over all inlets from the primary network indexed for that priority. There are some issues with this:

• It assumes that the inlets from the primary network are the only source, but there can be sources within the secondary network
• This doesn't work (directly) for LevelDemand within the secondary network

Alternative 3:

I don't understand this well enough at the moment to comment

[jarsarasty]

Hi @SouthEndMusic ,

Regarding your questions about the multi-commotidy formulation (priority-indexed flow variables approach ):

• It assumes that the inlets from the primary network are the only source, but there can be sources within the secondary network

I had not considered that case, but the formulation can easily be extended to include sources within the same network.

• This doesn't work (directly) for LevelDemand within the secondary network

Right. Some more elaboration is necessary to make it work with LevelDemands, but the basic concept remains the same.

I think the main challenge is the added size of the problem.

[SouthEndMusic]

The above problem was discussed with @jarsarasty, @simulutions and @gijsber. We came to the conclusion that we can do an approximation of the distribution of the flow from the primary network to the secondary network over the priorities, by:

1. Compute the total inflow over all inlets
2. Distribute the leftover demands over that total flow
3. Distribute those demands proportionally over the inlets

[gijsber]

Just to add:
the approximation will not work for very well for multiple inlets. Therefore we envision the end situation: #2563

[SouthEndMusic]

More detailed steps on how to proceed @simulutions:

• Remove the sources field of the AllocationModel struct (somewhat unrelated, but I'm pretty sure we don't need it anymore)
• In preprocess_demand_collection, set the capacity of the primary network to secondary network inlets to MAX_ABS_FLOW / scaling.flow (effectively unbounded). These are the flow variables with indices given by keys(allocation_model.subnetwork_demand)
• In postprocess_demand_collection:
  • compute the the demands as stated above, and store them in allocation_model.subnetwork_demand
  • set the demands in the problem for the primary network (which can be obtained as first(p_independent.allocation.allocation_problems), would be good to @assert that it is the primary network with is_primary_network). See for inspiration set_demands! and add_subnetwork_demand! where the relevant parameters are defined

I think that should be it for the demand communication. I also added #2230 and #2383 to the PR description because I think this PR is the last piece needed for those.

[SouthEndMusic]

Plans have changed again. It was decided that callbacks are not within the scope of MultiObjectiveAlgorithms.jl. We can just make our own loop over the objectives. So here's the approach:

• We need new objectives for minimizing the flow from the inlets into the secondary network from the primary network (only if a primary network is present?). To this end, in add_demand_objectives!, after pushing the *_objective_expression to objective_expressions_all, add a new expression with the sum of all the flow variables for the primary network to secondary network inlets (these are not divided into the various priorities themselves). Store somewhere in the allocation_model which indices in the objective_expressions_all correspond to the new objectives (and possibly for which demand priority). Beware that this can give problems with the fact that empty objective expressions are filtered out at the end of the AllocationModel constructor.
• Optimization in the secondary network comes in the collect_demands and allocate types (see ``@enumx AllocationOptimizationType). When using the allocate` type (in the last loop in `update_allocation!`), the new objectives are not needed, so they need to be filtered out before passing the objectives to the model in `optimize!` (the function in Ribasim, not the one from JuMP).
• When optimizing with the collect_demand types, we need to write our own Lexicographic algorithm. I don't think we need preprocess_demand_allocation and postprocess_demand_collection! anymore, maybe it's best to write a dedicated optimize! method for this. In this method:
  • Loop over the objectives
  • Apply the Lexicographic algorithm by adding a constraint objective = objective_value after optimizing for the objective (apart from for the new objectives)
  • When the objective is of the newly added type, the relevant flow values are the sum of the demands from the primary network. Compute the demand of the current priority and set it in allocation_model.subnetwork_demand and in the problem for the primary network (see my previous post for details on how to do this)
  • Remove the constraints that you added for the Lexicographic algorithm (if that is not possible, you can also add all constraints at initialization and only modify them in and after the Lexicographic algorithm)
• Work out the function communicate_secondary_network_allocations!; here the amounts allocated to the secondary networks in the primary networks are set as bounds on the the inflows from the primary network to the secondary network in the secondary network problem

@simulutions I'm probably missing some details here, feel free to reach out.

[visr]

🎉

[jarsarasty]

🎉