Wasserversorgung

Problemstellung

Wir betrachten ein System von Aquädukten, die Wasser von drei Flüssen (Knoten F1, F2 und F3) in eine Großstadt (Knoten T) transportieren. Die anderen Knoten sind Verbindungspunkte im Wasserversorgungssystem der Stadt. Jeder Aquädukt hat eine gewisse Transportkapazität in Tausend m\(^3\) Wasser pro Tag. Die Stadt möchte einen Durchflussplan bestimmen, der den Wasserdurchfluss zur Stadt maximiert, siehe Abbildung Abbildung 1.

Abbildung 1: Wasserversorgung

Quelle: Hillier, Lieberman: Introduction to Operations Research. 10th edition, 2015.Hillier, Lieberman: Introduction to Operations Research. 10th edition, 2015. Problem 10.5-3, Seite 428 f.

Modellierung

Um dieses Problem als ein Maximum Flow Problem zu formulieren, fügen wir dem gerichteten Graphen noch eine Quelle S hinzu, die mit den Knoten F1, F2 und F3 verbunden ist. Als Kantenkapazitäten verwenden wir die gesamten Wassermengen der Flüsse, siehe Abbildung Abbildung 2.

Abbildung 2: Wasserversorgung mit Quelle

Wir formulieren das Maximum Flow Problem auf dem Graphen mit Kantenmenge \(E\) und Knotenmenge \(N\):

  • Entscheidungsvariablen: Fluss \(x_e\) durch Kante \(e \in E\)
  • Zielfunktion: Maximiere \(x_{D,T} + x_{E,T} + x_{F,T}\) oder äquivalent \(x_{S,F1} + x_{S,F2} + x_{S,F3}\)
  • Nebenbedingungen:
    • Flussbilanz für jeden Knoten außer Quelle und Senke
    • Kapazitätsbeschränkungen für jede Kante

Daten

Code 
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pyomo.environ as pyo
Code 
capacity = {
    ("S", "F1"): 140,
    ("S", "F2"): 150,
    ("S", "F3"): 150,
    ("F1", "A"): 75,
    ("F1", "B"): 65,
    ("F2", "A"): 40,
    ("F2", "B"): 50,
    ("F2", "C"): 60,
    ("F3", "B"): 80,
    ("F3", "C"): 70,
    ("A",  "D"): 60,
    ("A",  "E"): 45,
    ("B",  "D"): 70,
    ("B",  "E"): 55,
    ("B",  "F"): 45,
    ("C",  "E"): 70,
    ("C",  "F"): 90,
    ("D",  "T"): 120,
    ("E",  "T"): 190,
    ("F",  "T"): 130}

Implementierung

Code 
nodes = ["S", "F1", "F2", "F3", "A", "B", "C", "D", "E", "F", "T"]
junctions  = ["F1", "F2", "F3", "A", "B", "C", "D", "E", "F"]
edges = list(capacity.keys())
edges
[('S', 'F1'),
 ('S', 'F2'),
 ('S', 'F3'),
 ('F1', 'A'),
 ('F1', 'B'),
 ('F2', 'A'),
 ('F2', 'B'),
 ('F2', 'C'),
 ('F3', 'B'),
 ('F3', 'C'),
 ('A', 'D'),
 ('A', 'E'),
 ('B', 'D'),
 ('B', 'E'),
 ('B', 'F'),
 ('C', 'E'),
 ('C', 'F'),
 ('D', 'T'),
 ('E', 'T'),
 ('F', 'T')]
Code 
model = pyo.ConcreteModel()

model.N = pyo.Set(initialize=nodes)      # set of nodes
model.J = pyo.Set(initialize=junctions)  # set of junctions
model.E = pyo.Set(initialize=edges)      # set of edges

model.x = pyo.Var(model.E, domain=pyo.NonNegativeReals)

@model.Constraint(model.E)
def flow_upper_bound(model, *e):
    # Note: The asterisk in *e is needed to
    # pass the input argument e as a tuple
    return model.x[e] <= capacity[e]

model.flow = pyo.Objective(expr=
    sum(model.x["S", n] for n in ["F1", "F2", "F3"]),
    sense=pyo.maximize)

@model.Constraint(model.J)
def node_constraint(model, j):
    return sum(model.x[j, n] for n in model.N if (j, n) in model.E) - \
           sum(model.x[n, j] for n in model.N if (n, j) in model.E) == 0

# model.pprint()
Code 
solver = pyo.SolverFactory('cbc')
# solver = pyo.SolverFactory('glpk')
# solver = pyo.SolverFactory('appsi_highs')
# solver = pyo.SolverFactory('gurobi')

results = solver.solve(model, tee=False)
print(f"status = {results.solver.status}")

print(f"maximum flow = {pyo.value(model.flow):.2f} qm/d")
status = ok
maximum flow = 395.00 qm/d

Ergebnisse

Code 
x_sol_dict = model.x.extract_values()
# x_sol_dict

# formatted print out of solution:
for edge in edges:
    print(f"flow on edge {edge} = {x_sol_dict[edge]:.2f} qm/d")
flow on edge ('S', 'F1') = 140.00 qm/d
flow on edge ('S', 'F2') = 130.00 qm/d
flow on edge ('S', 'F3') = 125.00 qm/d
flow on edge ('F1', 'A') = 75.00 qm/d
flow on edge ('F1', 'B') = 65.00 qm/d
flow on edge ('F2', 'A') = 20.00 qm/d
flow on edge ('F2', 'B') = 50.00 qm/d
flow on edge ('F2', 'C') = 60.00 qm/d
flow on edge ('F3', 'B') = 55.00 qm/d
flow on edge ('F3', 'C') = 70.00 qm/d
flow on edge ('A', 'D') = 50.00 qm/d
flow on edge ('A', 'E') = 45.00 qm/d
flow on edge ('B', 'D') = 70.00 qm/d
flow on edge ('B', 'E') = 55.00 qm/d
flow on edge ('B', 'F') = 45.00 qm/d
flow on edge ('C', 'E') = 70.00 qm/d
flow on edge ('C', 'F') = 60.00 qm/d
flow on edge ('D', 'T') = 120.00 qm/d
flow on edge ('E', 'T') = 170.00 qm/d
flow on edge ('F', 'T') = 105.00 qm/d