Allocation based on solving max flow problems

[SouthEndMusic]

The PoC code in https://github.com/Deltares/Ribasim/pull/627 attempts to propagate demands from users to the source. The approach there is flawed, because it does not account for multiple demand sharing paths from users to the source. Therefore here is a new setup for this demand propagation.

For now we assume a single source from which each user gets its demand.

Formulate a max flow problem per priority

A max flow problem consists of a directed graph with capacities associated with the edges, and a single source and target node between which one wants to find the maximum flow as restricted by the capacities of the edges. The Ribasim model graph (which I will call the model graph) can be converted into a graph for the maximum flow problem (which I will call the MFP graph) as follows:

• These are the nodes in the MFP graph:
  • The source node given by the source in the model graph;
  • Nodes in the model graph that form junctions (e.g. nodes 2 and 6 in the model and the end of this comment);
  • A single (sink) node that represents all users of a given priority (see here for multi-sink maximum flow problems);
• Each edge is given by a path that connects 2 nodes mentioned above in the model graph. The capacity of an edge is given by the smallest max_flow_rate of the pumps and outlets on this path. If there are no such nodes on the path, either the capacity is considered infinite, or if this cannot be supported, merge the MFP graph nodes on either side of this edge into one MFP graph node. Other node types that restrict flow (e.g. TabulatedRatingCurve) are not taken into consideration at this stage, as discussed with @gijsber.

There are a few capacity constraints that are not included into the above, which are;

• The sum total demand of the users of a given priority;
• The flow rate of the source.

To include these constraints into the MFP graph, an extra edge (and phantom node) can be added directly after the source node whose capacity is the minimum of the capacity constraints above.

Solve max flow problems

Obtaining the demand propagation from the users to the source then goes as follows:

• Solve the priority 1 max flow problem with an appropriate solver;
• For lower priority users the capacity taken up by higher priority users has to be taken into account. Therefore, there where the current MFP graph and MFP graphs of higher priority have overlap, subtract the capacity used by the solution of the MFP of higher priorities from the capacity of the current MFP graph before solving.

[visr]

Looks good! I was playing around a bit with max flow and plotting, look at this:

This show a max flow solution from node 1 to node 8 on top of the capacities. The flow and capacity is written in the edge labels and indicated in the edge thickness.

#=
https://github.com/MakieOrg/GraphMakie.jl
https://github.com/JuliaGraphs/GraphsFlows.jl
=#

#=
https://github.com/MakieOrg/GraphMakie.jl
https://github.com/JuliaGraphs/GraphsFlows.jl
=#

using GLMakie, GraphMakie, Graphs, GraphsFlows
using SparseArrays

flow_edges = [
    (1, 2, 10),
    (1, 3, 5),
    (1, 4, 15),
    (2, 3, 4),
    (2, 5, 9),
    (2, 6, 15),
    (3, 4, 4),
    (3, 6, 8),
    (4, 7, 16),
    (5, 6, 15),
    (5, 8, 10),
    (6, 7, 15),
    (6, 8, 10),
    (7, 3, 6),
    (7, 8, 10),
]
nv = reduce(max, max(e[1], e[2]) for e in flow_edges)
g = DiGraph(nv)

capacity_matrix = spzeros(Int, nv, nv)  # Create a capacity matrix

for (u, v, f) in flow_edges
    add_edge!(g, u, v)
    capacity_matrix[u, v] = f
end

from = 1
to = 8
f, F = maximum_flow(g, from, to, capacity_matrix)

begin
    ilabels = string.(1:nv)
    node_color = fill(:gray80, nv)
    node_color[from] = :lightgreen
    node_color[to] = :tomato
    elabels =
        [string(F[src(e), dst(e)], " / ", capacity_matrix[src(e), dst(e)]) for e in edges(g)]

    edge_width = 10.0 * [capacity_matrix[src(e), dst(e)] for e in edges(g)] / maximum(capacity_matrix)
    arrow_size = 4 * edge_width

    # draw the edge capacity in grey
    fig = graphplot(
        g;
        edge_width,
        edge_color = :gray80,
        arrow_size,
        arrow_shift=:end,
    )
    # draw the maxflow over the capacity in blue
    edge_width = 10.0 * [F[src(e), dst(e)] for e in edges(g)] / maximum(capacity_matrix)
    graphplot!(
        g;
        ilabels,
        elabels,
        node_marker = Circle,
        node_strokewidth = 1,
        node_color,
        elabels_distance = 15,
        edge_width,
        edge_color = :blue,
        arrow_show=false,
    )
    display(fig)
end

[SouthEndMusic]

Some additional notes on the above proposal:

Combining all user nodes into a single sink in the MFP graph does not directly work, because we care about how the flow is distributed over the different paths that lead to this sink, corresponding to the paths to the different user nodes in the model graph. The way this can be fixed is by keeping the extra edge after the source but only posing the source flow as its capacity, and adding extra edges at the ends of the paths to the sink with a capacity equal to the demand of the corresponding user.

The above can in fact be used to have only a single MFP graph topology covering all MFP problems for each priority. The capacity of these extra user edges can be adapted to the demand of the user of a given priority. This also makes supporting having demands for different priorities from a single user trivial.

I'll add a figure of how I imagine the MFP graph to look for our subnetwork test model.

[SouthEndMusic]

@visr that looks nice. What does the begin ... end scope in your code example do?

[visr]

The additional notes make sense to me. The begin end block doesn't do anything and can be removed. It also doesn't change the scope. I just use it sometimes to group lines that I want to always run together in VSCode. Cell delimiters, special comments like ## also work.

[SouthEndMusic]

[SouthEndMusic]

I thought of another edge case to keep in mind: if the junction itself is a pump/outlet, then a node and edge should be added before it in the MFP graph to account for its capacity restriction.

[SouthEndMusic]

When the flows on the user demand edges in the in the MFP solution are maximal, all demand to the users (of the considered priority) can be provided. If not, the solution already has made a choice about how the available flow is distributed over the users. Not in the ways we have talked about to make it fair, but in a way that takes the capacities in the MFP graph into account. We have to make a decision about how we handle this.

[SouthEndMusic]

This issue starts to look like an epic, the method here does the full allocation instead of only demand propagation from users to the source. Maybe sub-issues are required.

[SouthEndMusic]

I think it is not needed to add the source flow explicitly in the max flow graph. It suffices to add the user demands for a given priority just before the sink. Then iterating over the priorities looks like this:

1. Compute the total source capacity (e.g. flow boundary flow rate at current time in simulation);
2. Set the initial edge capacities of the max flow graph; all are $\infty$ apart from those that have a pump/outlet in the model graph;
3. Assign the priority 1 user demands to the edges in the max flow graph between the users and the sink;
4. Solve the priority 1 max flow problem;
5. Compute the source capacity taken up by priority 1 as the sum of the flows between the users and the sink in the max flow solution;
6. Subtract the source capacity taken up by priority 1 from the total source capacity;
7. Subtract the flows used by the priority 1 max flow solution from the max flow graph edge capacities (edges with infinity capacity will keep their infinity capacity);
8. Repeat steps 3-7 with successive priorities until at step 6 it is detected that the demand that can still reach the users at some priority is larger than the total source capacity that is left.

Note that it can happen that there is not even enough source capacity for the Demands of priority 1. This is not a problem for this approach.

Below is a formal description of the problem of distributing the source capacity that is left over the users.

We are given:

• $P$: the priority at which the sum of flows between the users and the sink in the max flow solution is larger than the source capacity that is left;
• $q_i \geq 0, \quad \text{for }1\le i\le n$: The flows between the users and the sink in the max flow solution of priority $P$, with $n$ the number of users. $q_i$ is bounded above by the demand of user $i$ at priority $P$, but it can be smaller than that because of constraints in the max flow graph;
• $Q_\text{left}$: the source capacity that is left after allocating source capacity to the users for priorities $1,\ldots,P-1$.

We know by definition:

$$ \sum_{i=1}^n qi > Q\text{left}. $$

We want to find:

$0 \le x_i \le q_i, \quad \text{for }1\le i\le n$ (the allocations to the users at priority $P$) such that

$$ \sum_{i=1}^n xi = Q\text{left}, $$

(i.e. all leftover source capacity is used) and the distribution is 'fair'. Maybe a good measure of fairness is

$$ \min_{i: q_i > 0} \frac{x_i}{q_i}, $$

e.g. the solution maximizes the minimum fraction a user gets of its demand of priority $P$ (that is, the part of its demand that can actually reach the user).

[gijsber]

Hi @SouthEndMusic , looks promising. I assume demands are first corrected for local available water, before you enter this process.

[gijsber]

By local available water, I mean volume in the basin above the minimum water level