from dataclasses import dataclass, field
from plot import draw_network, draw_routing
from report import report_output
@dataclass
class Network:
inlet_contaminants: list
removal: list
triplets: list
max_concentration: list
inlet_flow: list
positions: dict
num_contaminants: int = field(init=False)
num_inlet_splitters: int = field(init=False)
num_units: int = field(init=False)
num_outlet_mixers: int = field(init=False)
def __post_init__(self):
self.num_contaminants = len(self.max_concentration)
self.num_inlet_splitters = len(self.inlet_flow)
self.num_units = len(self.removal)
self.num_outlet_mixers = 1
@property
def num_splitters(self):
return self.num_inlet_splitters + self.num_units
@property
def num_mixers(self):
return self.num_units + self.num_outlet_mixers
@property
def outlet_mixer(self):
return self.num_units
@property
def inlet_splitters(self):
return range(self.num_inlet_splitters)
@property
def outlet_splitters(self):
return range(self.num_inlet_splitters, self.num_splitters)
@property
def splitters(self):
return range(self.num_splitters)
@property
def mixers(self):
return range(self.num_mixers)
@property
def units(self):
return range(self.num_units)
@property
def contaminants(self):
return range(self.num_contaminants)
@dataclass
class Routing:
treated_flow: float
flow: dict
contaminant_flow: dict
unit_flow: dict
unit_inflow: dict
unit_outflow: dictWastewater Treatment Network
Route wastewater through splitters, mixers and treatment units so that the least possible flow goes through costly treatment and the discharge still meets its quality limits.
Introduction
We consider the wastewater treatment network problem. The goal of a wastewater treatment network is to lower the contaminant concentration in wastewater to the required quality standards. Such a network contains splitters, mixers and treatment units. Splitters divide incoming flows, mixers combine them, and treatment units remove a fixed share of each contaminant from the flow passing through them. The problem is to determine how much of the flow from each splitter should go to each mixer, so as to minimize the flow sent to the expensive treatment units while the final output still meets the required quality standards.
Below is the smallest network the page solves, with each arc as wide as the flow it carries and colored by how dirty that flow is. It starts with nothing treated, so both streams run straight to the output and the gauge on the right stands above the dashed rule that marks the discharge limit. Routing part of the flow through the treatment unit brings it back down onto the rule, and no further.
A worked example
First, let’s illustrate the problem on a small instance containing
- 2 inlet splitters,
- 1 triplet inlet mixer / treatment unit / output splitter,
- 1 contaminant.
Let’s visualize the network:
We consider the following input values:
| Input | Value |
|---|---|
| Total flow in Splitter \(i_0\) | \(20{,}000\) kg/h |
| ↳ contaminant flow in Splitter \(i_0\) | \(3\) kg/h |
| Total flow in Splitter \(i_1\) | \(10{,}000\) kg/h |
| ↳ contaminant flow in Splitter \(i_1\) | \(2\) kg/h |
| Quality percentage in the output mixer | \(0.01\%\) of the total input flow, that is \((20{,}000 + 10{,}000) \times 10^{-4} = 3\) kg/h |
| Treatment unit removal ratio | \(0.9\) |
Possible solution. We consider the following solution:
- \(7.41\%\) of the flow from Splitter \(i_0\) goes into Mixer \(0\).
- The rest of the flow from Splitter \(i_0\) goes into Mixer out.
- All the flow from Splitter \(i_1\) goes into Mixer \(0\).
- All the flow from Splitter \(0\) goes to Mixer out.
Let’s visualize this solution:
At the entrance of Mixer \(0\):
- Flow:
- flow from Splitter \(i_1\): \(10{,}000\) kg/h;
- flow from \(7.41\%\) of Splitter \(i_0\): \(1{,}481\) kg/h;
- thus the total flow into Mixer \(0\) is \(10{,}000 + 1{,}481 = 11{,}481\) kg/h.
- Contaminant rate:
- contaminant rate from Splitter \(i_1\): \(2\) kg/h;
- contaminant rate from \(7.41\%\) of Splitter \(i_0\): \(0.22\) kg/h;
- thus the total contaminant rate entering Mixer \(0\) is \(2 + 0.22 = 2.22\) kg/h.
After Treatment Unit \(0\), the contaminant rate is reduced to \(2.22 \times (1 - 0.9) = 0.22\) kg/h.
At the entrance of Mixer out:
- the treated flow, carrying \(0.22\) kg/h of contaminant, is merged with the remaining \(92.59\%\) of the flow from Splitter \(i_0\);
- the remaining contaminant rate from Splitter \(i_0\) is \(2.78\) kg/h;
- thus the total contaminant rate in Mixer out is \(0.22 + 2.78 = 3\) kg/h.
The final contaminant rate meets the required quality of \(3\) kg/h. Finding such a split by hand stops working almost immediately, so let’s now give the formal description of the problem.
Problem description
Input
- A set of contaminants \(C\).
- A set of inlet splitters \(S_I\).
- Total flow \(\beta^\text{IN}_{s_i}\) (kg/h) for each inlet splitter \(s_i \in S_I\).
- Contaminant flow \(\alpha^\text{IN}_{c,s_i}\) (kg/h) for each contaminant \(c \in C\) and each inlet splitter \(s_i \in S_I\).
- An output mixer \(m_O\), with an upper limit of discharge concentration \(\theta^\text{OUT}_c\) for each contaminant \(c \in C\).
- A set of triplets \(M_T \times T \times S_T\), each triplet \((m,t,s)\) corresponding to a treatment unit with its input mixer and output splitter, with a removal ratio \(\gamma_{t,c}\) for each contaminant \(c \in C\) and each treatment unit \(t \in T\).
Problem: choose how to divide the flow at the splitters in order to respect the output quality.
Objective: the goal is to minimize the total flow entering the treatment units, as operational costs are heavily influenced by the flow rates entering these units,
\[ \sum_{t \in T} f_t \]
where \(f_t\) represents the total flow going through treatment unit \(t \in T\).
Mathematical model
Let’s write a mathematical model of the problem.
Variables
The problem is a variant of a network flow problem. There are variables to represent the water flow and the contaminant flow on each arc of the network:
- \(f_{s,m}\), \(s \in S_I \cup S_T\), \(m \in M_T \cup \{m_O\}\): total flow (kg/h) from splitter \(s\) to mixer \(m\).
- \(f_{s,m,c}\), \(s \in S_I \cup S_T\), \(m \in M_T \cup \{m_O\}\), \(c \in C\): flow (kg/h) of contaminant \(c\) from splitter \(s\) to mixer \(m\).
- \(f_t\), \(t \in T\): total flow (kg/h) through treatment unit \(t\).
- \(f^i_{t,c}\), \(t \in T\), \(c \in C\): flow (kg/h) of contaminant \(c\) entering treatment unit \(t\).
- \(f^o_{t,c}\), \(t \in T\), \(c \in C\): flow (kg/h) of contaminant \(c\) leaving treatment unit \(t\).
- \(x_{s,m}\), \(s \in S_I \cup S_T\), \(m \in M_T \cup \{m_O\}\): fraction of the flow at splitter \(s\) that goes to mixer \(m\). This variable is what links the water and contaminant flows.
Objective: minimize the total flow entering the treatment units,
\[ \min \sum_{t \in T} f_t \]
Constraints
Between splitters and mixers.
Mass balance at the inlet splitters: \[ \forall s \in S_I, \quad \forall m \in M_T \cup \{m_O\}, \qquad \beta^\text{IN}_{s} x_{s,m} = f_{s,m} \]
Contaminant mass balance at the inlet splitters: \[ \forall s \in S_I, \quad \forall m \in M_T \cup \{m_O\}, \quad \forall c \in C, \qquad \alpha^\text{IN}_{c,s} x_{s,m} = f_{s,m,c} \]
Mass balance at the outlet splitters: \[ \forall (m,t,s) \in M_T \times T \times S_T, \quad \forall m' \in M_T \cup \{m_O\}, \qquad f_t x_{s,m'} = f_{s,m'} \]
Contaminant mass balance at the outlet splitters: \[ \forall (m,t,s) \in M_T \times T \times S_T, \quad \forall m' \in M_T \cup \{m_O\}, \quad \forall c \in C, \qquad f^o_{t,c} x_{s,m'} = f_{s,m',c} \]
Flow balance: \[ \forall s \in S_I \cup S_T, \qquad \sum_{m \in M_T \cup \{m_O\}} x_{s,m} = 1 \]
Flow conservation at the outlet splitters: \[ \forall (m,t,s) \in M_T \times T \times S_T, \qquad f_t = \sum_{m' \in M_T \cup \{m_O\}} f_{s,m'} \]
Contaminant flow conservation at the outlet splitters: \[ \forall (m,t,s) \in M_T \times T \times S_T, \quad \forall c \in C, \qquad f^o_{t,c} = \sum_{m' \in M_T \cup \{m_O\}} f_{s,m',c} \]
At the inlet mixers.
Flow conservation: \[ \forall (m,t,s) \in M_T \times T \times S_T, \qquad f_t = \sum_{s' \in S_I \cup S_T} f_{s',m} \]
Contaminant flow conservation: \[ \forall (m,t,s) \in M_T \times T \times S_T, \quad \forall c \in C, \qquad f^i_{t,c} = \sum_{s' \in S_I \cup S_T} f_{s',m,c} \]
Removal of contaminants in the treatment units: \[ \forall (m,t,s) \in M_T \times T \times S_T, \quad \forall c \in C, \qquad f^o_{t,c} = (1 - \gamma_{t,c}) f^i_{t,c} \]
Output quality constraint: \[ \forall c \in C, \qquad \sum_{s \in S_I \cup S_T} f_{s,m_o,c} \leq \theta^\text{OUT}_{c} \sum_{s_i \in S_I} \beta^\text{IN}_{s_i} \]
The mass balance constraints are quadratic equalities. Therefore, the problem is a nonconvex quadratically constrained quadratic problem (nonconvex QCQP), and multistart is what gets Knitro past the local optima that come with it.
Input data
A Network names the contaminants, the inlet splitters with their flow and its contaminant content, the treatment triplets with their removal ratios, and the discharge limits. Each of the three solved on this page is written out beside the picture of the network it describes. A Routing holds the flow into the treatment units and every arc’s flow, per contaminant and in total.
An arc bends around any node lying between its two ends, so a flow returning up the network does not read as one crossing a unit. On the larger networks the nodes carry color instead of labels, and one legend names them.
The worked example is the smallest instance the model takes, with two inlet splitters feeding one mixer, one treatment unit, one output splitter and one contaminant. positions is there for the figure alone, since the model never asks where a node sits.
network_0 = Network(
inlet_contaminants=[[3], [2]],
removal=[[0.9]],
triplets=[(0, 0, 2)],
max_concentration=[1e-4],
inlet_flow=[20000, 10000],
positions={
"Splitter i0": [(0, 15), "blue"],
"Splitter i1": [(0, 5), "blue"],
"Mixer 0": [(2, 10), "red"],
"Treatment Unit 0": [(4, 9), "green"],
"Splitter 0": [(6, 10), "orange"],
"Mixer out": [(8, 10), "purple"],
},
)
draw_network(network_0)Model implementation
import pyomo.environ as pyo
SOLVER_NAME = "knitroampl"Everything the inlet splitters do is affine, since their total flow is data. The outlet splitters are where two variables meet, and those products go into the same ConstraintList as written.
def minimize_treatment(network, *, multistart=True, max_solves=None):
splitters, mixers = network.splitters, network.mixers
contaminants, units = network.contaminants, network.units
inlets, triplets = network.inlet_splitters, network.triplets
m_o = network.outlet_mixer
beta, alpha = network.inlet_flow, network.inlet_contaminants
gamma, theta = network.removal, network.max_concentration
model = pyo.ConcreteModel()
model.triplets = pyo.Set(initialize=triplets, dimen=3)
model.fsm = pyo.Var(splitters, mixers, within=pyo.NonNegativeReals)
model.fsmc = pyo.Var(splitters, mixers, contaminants, within=pyo.NonNegativeReals)
model.ft = pyo.Var(units, within=pyo.NonNegativeReals)
model.fitc = pyo.Var(model.triplets, contaminants, within=pyo.NonNegativeReals)
model.fotc = pyo.Var(model.triplets, contaminants, within=pyo.NonNegativeReals)
model.xsm = pyo.Var(splitters, mixers, bounds=(0, 1))
fsm, fsmc, ft, xsm = model.fsm, model.fsmc, model.ft, model.xsm
fitc, fotc = model.fitc, model.fotc
treated = pyo.quicksum(ft[t] for t in units)
model.obj = pyo.Objective(expr=treated, sense=pyo.minimize)
model.con = pyo.ConstraintList()
for s in inlets:
for m in mixers:
model.con.add(beta[s] * xsm[s, m] == fsm[s, m])
for c in contaminants:
model.con.add(alpha[s][c] * xsm[s, m] == fsmc[s, m, c])
for mi, t, s in triplets:
for m1 in mixers:
model.con.add(ft[t] * xsm[s, m1] == fsm[s, m1])
for c in contaminants:
model.con.add(fotc[mi, t, s, c] * xsm[s, m1] == fsmc[s, m1, c])
for s in splitters:
model.con.add(pyo.quicksum(xsm[s, m] for m in mixers) == 1)
for mi, t, s in triplets:
model.con.add(ft[t] == pyo.quicksum(fsm[s, m] for m in mixers))
model.con.add(ft[t] == pyo.quicksum(fsm[sp, mi] for sp in splitters))
for c in contaminants:
f_in, f_out = fitc[mi, t, s, c], fotc[mi, t, s, c]
model.con.add(f_out == pyo.quicksum(fsmc[s, m, c] for m in mixers))
model.con.add(f_in == pyo.quicksum(fsmc[sp, mi, c] for sp in splitters))
model.con.add(f_out == (1 - gamma[t][c]) * f_in)
for c in contaminants:
discharged = pyo.quicksum(fsmc[sp, m_o, c] for sp in splitters)
model.con.add(discharged <= sum(beta[s] * theta[c] for s in inlets))
solver = pyo.SolverFactory(SOLVER_NAME)
solver.options["ms_enable"] = int(multistart)
if max_solves is not None:
solver.options["ms_maxsolves"] = max_solves
solver.solve(model, tee=True)
treated_flow = model.obj()
flows = [{k: v() for k, v in f.items()} for f in (fsm, fsmc, ft, fitc, fotc)]
return Routing(treated_flow, *flows)Solving the worked example
The introduction works that network out by hand. Knitro should reach the same split.
routing_0 = minimize_treatment(network_0)Artelys Knitro 16.0.0: ms_enable=1
=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro using up to 8 threads.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 11 variables (52%) and 11 constraints (52%) in 0.00s.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 0.001
findiff_numthreads 1
hessian_no_f 1
hessopt 1
ms_enable 1
opttol 1e-06
opttol_abs 0.001
Problem Characteristics | Presolved
-----------------------
Problem type: QCQP
Objective: minimize / linear
Number of variables: 21 | 10
bounds: lower upper range | lower upper range
15 0 6 | 6 0 4
free fixed | free fixed
0 0 | 0 0
Number of constraints: 21 | 10
eq. ineq. range | eq. ineq. range
linear: 16 1 0 | 5 1 0
quadratic: 4 0 0 | 4 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 1 45 | 1 23
quadratic: 0 8 4 | 0 8 4
total: 1 53 4 | 1 31 4
Knitro parallel multistart will run with 8 threads.
Return codes description
------------------------
0: The final solution satisfies the termination conditions for verifying optimality.
-100 to -199: A feasible approximate solution was found.
-200 to -299: Knitro terminated at an infeasible point.
-300 to -301: The problem was determined to be unbounded.
-400 to -499: Knitro terminated because it reached a pre-defined limit.
-400 to -409: A feasible point was found.
-410 to -419: No feasible point was found.
-500 to -599: Knitro terminated with an input error or some non-standard error.
A more detailed description of individual return codes and their corresponding
termination messages is provided at
https://www.artelys.com/app/docs/knitro/3_referenceManual/returnCodes.html
Solve Thrd Status Objective FeasError Opt Error Solve Time Real Time
----- ---- ------ ------------ ----------- ----------- ----------- -----------
1 3 0 11481.5 1.00004e-10 1.00004e-10 3.24011e-02 0.152909
4 1 0 11481.5 1.00004e-10 1.00004e-10 3.63932e-02 0.152909
2 7 0 11481.5 8.34888e-10 1.00048e-10 4.11059e-02 0.166615
3 5 0 11481.5 8.34888e-10 1.00048e-10 4.22931e-02 0.166615
11 5 0 11481.5 9.99973e-11 9.99973e-11 3.58299e-02 0.167661
21 7 0 11481.5 1.00303e-10 6.16849e-09 1.80995e-03 0.167661
31 3 0 11481.5 1.00000e-10 1.00000e-10 1.48299e-02 0.167661
13 4 0 11481.5 9.99854e-11 9.58607e-09 9.95272e-03 0.168013
8 3 0 11481.5 4.32158e-08 6.09998e-09 6.00920e-02 0.168013
28 3 0 11481.5 1.00000e-10 1.00000e-10 3.04912e-03 0.168013
10 7 0 11481.5 1.01486e-10 2.03180e-08 1.86347e-02 0.168365
18 2 0 11481.5 5.38730e-08 1.84683e-07 5.85879e-03 0.168365
22 7 0 11481.5 5.83190e-09 5.94482e-08 3.77070e-03 0.168365
24 6 0 11481.5 9.98317e-11 6.27551e-09 4.25190e-03 0.168365
29 6 0 11481.5 1.10084e-06 3.05377e-07 4.84280e-03 0.168365
32 6 0 11481.5 2.07516e-10 1.00014e-10 2.71818e-03 0.168365
6 6 0 11481.5 9.91847e-11 2.45321e-09 4.95243e-02 0.169011
25 6 0 11481.5 2.83990e-07 2.44632e-08 6.78933e-03 0.169011
15 2 0 11481.5 2.93572e-10 2.26112e-08 1.06101e-02 0.172904
9 1 0 11481.5 1.00004e-10 1.00004e-10 5.84700e-02 0.172904
MULTISTART: Best locally optimal solution is returned.
EXIT: Multi-start stopped because of a low estimated probability of finding
an unobserved solution. Set ms_terminate=0 to disable multi-start rule-based
termination procedure.
20 solve(s) returned satisfactory solutions.
Final Statistics
----------------
Final objective value = 1.14814814687429e+04
Final feasibility error (abs / rel) = 1.10e-06 / 6.82e-11
Final optimality error (abs / rel) = 3.05e-07 / 3.05e-07
# of iterations = 878
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 0.19335 ( 8.607 CPU time)
================================================================================
Show the full outputHide the full output
draw_routing(network_0, routing_0, contaminants=[0])We found the same result as previously.
Larger example
The first part features a highly minimalist model. The next one introduces more realistic elements, with three inlet splitters, three treatment units and three contaminants. Every splitter may feed every mixer, so the network below is dense on purpose, since it is the set of routes on offer and not the ones that get used. The dataset comes from Cheng and Li, 2020 DOI.
network_1 = Network(
inlet_contaminants=[[12, 10, 10], [6, 3, 1.5], [1, 5, 1]],
removal=[[0.9, 0, 0], [0, 0.99, 0], [0, 0, 0.8]],
triplets=[(0, 0, 3), (1, 1, 4), (2, 2, 5)],
max_concentration=[1e-4, 1e-4, 1e-4],
inlet_flow=[20000, 15000, 5000],
positions={
"Splitter i0": [(0, 10), "blue"],
"Splitter i1": [(0, 8), "blue"],
"Splitter i2": [(0, 6), "blue"],
"Mixer 0": [(2, 18), "red"],
"Mixer 1": [(2, 14), "red"],
"Mixer 2": [(2, 1), "red"],
"Treatment Unit 0": [(4, 17), "green"],
"Treatment Unit 1": [(4, 13), "green"],
"Treatment Unit 2": [(4, 0), "green"],
"Splitter 0": [(6, 18), "orange"],
"Splitter 1": [(6, 14), "orange"],
"Splitter 2": [(6, 1), "orange"],
"Mixer out": [(8, 8), "purple"],
},
)
draw_network(network_1)routing_1 = minimize_treatment(network_1)Artelys Knitro 16.0.0: ms_enable=1
=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro using up to 8 threads.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 57 variables (40%) and 57 constraints (41%) in 0.00s.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 0.001
findiff_numthreads 1
hessian_no_f 1
hessopt 1
ms_enable 1
opttol 1e-06
opttol_abs 0.001
Problem Characteristics | Presolved
-----------------------
Problem type: QCQP
Objective: minimize / linear
Number of variables: 141 | 84
bounds: lower upper range | lower upper range
117 0 24 | 63 0 21
free fixed | free fixed
0 0 | 0 0
Number of constraints: 138 | 81
eq. ineq. range | eq. ineq. range
linear: 87 3 0 | 30 3 0
quadratic: 48 0 0 | 48 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 3 348 | 3 234
quadratic: 0 96 48 | 0 96 48
total: 3 444 48 | 3 330 48
Knitro parallel multistart will run with 8 threads.
Return codes description
------------------------
0: The final solution satisfies the termination conditions for verifying optimality.
-100 to -199: A feasible approximate solution was found.
-200 to -299: Knitro terminated at an infeasible point.
-300 to -301: The problem was determined to be unbounded.
-400 to -499: Knitro terminated because it reached a pre-defined limit.
-400 to -409: A feasible point was found.
-410 to -419: No feasible point was found.
-500 to -599: Knitro terminated with an input error or some non-standard error.
A more detailed description of individual return codes and their corresponding
termination messages is provided at
https://www.artelys.com/app/docs/knitro/3_referenceManual/returnCodes.html
Solve Thrd Status Objective FeasError Opt Error Solve Time Real Time
----- ---- ------ ------------ ----------- ----------- ----------- -----------
4 3 0 92073.0 3.65253e-09 1.25108e-07 4.45657e-02 0.491302
1 0 0 92073.0 3.65253e-09 1.25108e-07 4.80090e-02 0.491302
9 0 0 85141.0 1.00072e-10 1.00072e-10 8.63341e-02 0.544927
24 2 0 90268.5 3.31440e-09 2.38763e-08 5.54337e-02 0.544927
17 4 0 81567.0 1.00002e-10 1.00002e-10 5.57973e-02 0.557122
18 2 0 90268.5 2.69574e-09 6.98212e-07 6.97831e-02 0.569622
12 2 0 90268.5 1.00017e-10 1.00017e-10 7.93532e-02 0.579028
22 6 0 80778.9 1.27991e-08 4.54936e-07 7.22430e-02 0.579028
11 4 0 80778.9 1.00004e-10 3.67229e-08 8.37134e-02 0.585587
27 6 0 80778.9 8.13725e-09 4.89626e-09 0.112363 0.585587
19 3 0 81567.0 1.19348e-09 4.61667e-09 9.63295e-02 0.592150
0 2 0 85141.0 3.80828e-10 1.08412e-09 0.192531 0.599020
28 7 0 90268.5 8.10142e-08 3.38638e-07 0.115459 0.599020
26 2 0 87912.8 2.12090e-10 3.99564e-09 0.115139 0.606933
30 4 0 92073.0 9.73960e-08 4.08422e-08 0.100862 0.606933
14 1 0 92073.0 1.00005e-10 1.00005e-10 0.128034 0.614096
25 3 0 85141.0 2.71472e-08 1.74499e-08 0.120544 0.614096
6 4 0 80778.9 1.16727e-08 1.60799e-09 0.165301 0.621959
7 6 0 85141.0 1.00000e-10 1.00000e-10 0.197927 0.621959
13 6 0 92073.0 1.00001e-10 1.00001e-10 0.132059 0.621959
29 5 0 85141.0 1.23691e-10 5.83477e-08 0.146435 0.646748
31 0 0 81567.0 1.00049e-10 7.55953e-10 0.149354 0.646748
5 1 0 80778.9 9.99529e-11 6.99401e-10 0.204945 0.731050
10 0 0 85141.0 2.40470e-09 2.77842e-08 0.176206 0.738220
21 0 0 92073.0 5.63887e-10 3.69955e-08 0.169203 0.738220
2 7 0 80778.9 2.57347e-08 7.59701e-09 0.214300 0.745968
3 5 0 80778.9 2.57347e-08 7.59701e-09 0.218876 0.745968
20 4 0 80778.9 1.96451e-10 2.00725e-09 0.167116 0.754555
15 7 0 92073.0 1.69838e-07 1.13552e-08 0.208454 0.815835
16 5 0 90268.5 3.03066e-10 1.33868e-10 0.221211 0.825798
8 3 0 80778.9 9.01195e-08 5.92734e-08 0.246860 0.832430
35 2 0 85141.0 6.93159e-08 6.45247e-09 0.112644 0.839613
23 1 0 85141.0 2.62208e-09 3.37035e-07 0.307745 0.847015
40 7 0 85141.0 1.49062e-07 8.54941e-09 7.62113e-02 0.847015
36 7 0 80778.9 4.41047e-07 1.39258e-07 0.108071 0.923423
39 1 0 81567.0 1.48176e-07 2.76558e-09 0.107867 0.923423
42 0 0 80778.9 1.00000e-10 1.00000e-10 7.48216e-02 0.935952
38 0 0 80778.9 3.49246e-10 4.05584e-08 0.113848 0.950568
32 5 0 92073.0 1.70420e-08 2.81452e-08 0.174267 0.961352
48 4 0 80778.9 3.89264e-10 1.90539e-10 4.99990e-02 0.999490
45 5 0 81567.0 2.98314e-10 1.82385e-10 0.108388 1.01540
44 3 0 80778.9 1.00009e-10 1.00009e-10 0.133783 1.01540
46 7 0 81567.0 7.89574e-10 1.00009e-10 0.128065 1.01540
33 3 0 85141.0 1.00001e-10 1.00001e-10 0.164625 1.05312
49 6 0 90268.5 9.99982e-11 1.19221e-09 8.86041e-02 1.18242
51 4 0 85141.0 9.99973e-11 9.99973e-11 9.45903e-02 1.19054
34 6 0 81567.0 1.18884e-06 7.90391e-09 0.211987 1.19934
37 4 0 80778.9 9.99938e-11 2.32535e-09 0.212877 1.21101
52 3 0 85141.0 5.89443e-10 4.69924e-10 0.139560 1.21101
50 5 0 81567.0 5.63716e-10 3.36818e-08 0.174258 1.21921
47 0 0 80778.9 2.42478e-09 1.09449e-09 0.123709 1.23090
59 3 0 85141.0 6.40284e-10 1.44080e-08 9.51969e-02 1.31607
57 1 0 85141.0 1.55391e-10 2.03354e-09 9.88014e-02 1.32242
53 7 0 85141.0 7.54983e-10 2.35559e-10 0.238719 1.32242
60 5 0 81567.0 6.37854e-09 1.00082e-10 0.106267 1.32242
43 1 0 81567.0 1.78261e-08 1.28455e-09 0.282413 1.43704
62 3 0 80778.9 9.99938e-11 6.01631e-10 0.113838 1.43704
66 0 0 85141.0 7.53062e-10 2.16645e-07 6.40093e-02 1.44400
41 2 0 81567.0 5.96013e-10 5.03365e-07 0.304608 1.51012
70 4 0 92073.0 1.00012e-10 1.00012e-10 5.67633e-02 1.51012
54 6 0 90268.5 9.99965e-11 9.34525e-10 0.256591 1.52503
55 0 0 80778.9 3.97306e-07 2.98519e-07 0.246437 1.53870
63 7 0 80778.9 1.43442e-09 2.14801e-09 0.266141 1.57911
58 2 0 80778.9 1.44678e-09 3.65376e-08 0.228208 1.58790
65 6 0 80778.9 7.79180e-09 2.61610e-10 0.200385 1.58790
68 2 0 80778.9 6.45377e-09 4.56412e-09 0.195163 1.59711
64 5 0 92073.0 9.99973e-11 9.99973e-11 0.224808 1.62604
72 5 0 80778.9 7.05768e-10 1.16163e-07 0.146624 1.62604
75 2 0 85141.0 7.85304e-08 4.29072e-08 0.119658 1.68792
61 1 0 85141.0 3.13776e-10 6.38101e-09 0.347872 1.69524
67 3 0 85141.0 1.00009e-10 1.00009e-10 0.217187 1.69524
71 6 0 85141.0 4.47471e-10 7.26031e-08 0.147092 1.70285
78 0 0 80778.9 9.99947e-11 9.99947e-11 8.25595e-02 1.72501
79 5 0 85141.0 2.70679e-07 7.53030e-07 3.74026e-02 1.72501
69 0 0 92073.0 2.13541e-08 2.13186e-09 0.219363 1.81873
82 0 0 87912.8 1.43700e-09 4.44834e-07 6.22764e-02 1.82824
80 6 0 80778.9 1.00001e-10 1.00001e-10 8.27178e-02 1.83646
MULTISTART: Best locally optimal solution is returned.
EXIT: Multi-start stopped because of a low estimated probability of finding
an unobserved solution. Set ms_terminate=0 to disable multi-start rule-based
termination procedure.
77 solve(s) returned satisfactory solutions.
Final Statistics
----------------
Final objective value = 8.07789201450768e+04
Final feasibility error (abs / rel) = 4.41e-07 / 2.88e-11
Final optimality error (abs / rel) = 1.39e-07 / 1.39e-07
# of iterations = 16620
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 1.97062 ( 21.829 CPU time)
================================================================================
Show the full outputHide the full output
draw_routing(network_1, routing_1)report_output(network_1, routing_1)contaminant discharged limit of limit
(kg/h) (kg/h)
--------------------------------------------------
1 4.00 4.00 100.0%
2 4.00 4.00 100.0%
3 4.00 4.00 100.0%The flow entering Splitter \(i_0\) is highly concentrated in contaminants 1, 2 and 3, so it is essential that it first passes through the three treatment units to reduce the concentration enough to meet the quality criteria. The flow entering Splitter \(i_1\) is very concentrated in contaminant 1, which is why a portion of it goes to Treatment Unit 0; a small portion can go straight to the output mixer since it is not very concentrated overall, so the criteria are met without sending the whole flow through the units. Splitter \(i_2\) is highly concentrated, particularly in contaminant 2, so it should mainly pass through Treatment Unit 1, though it must still pass through the other units to help stay below the threshold. The threshold is reached exactly for all three contaminants.
Even larger example
Let’s try another example with a higher number of units, ten inlet splitters, eight treatment units and five contaminants. The contaminant flows, the removal ratios and the inlet flows were drawn at random, and the nodes sit on a regular grid. Every splitter may feed every mixer here too, for 178 arcs on offer, and the figure below draws every one of them. It is dense on purpose, since what it shows is the choice the model faces rather than a route through it. The solution that follows keeps flow on about a fifth of the splitter-to-mixer connections and leaves the rest empty. Left to itself, Knitro’s multistart tries a few hundred starting points on a network this size, taking minutes and buying half a percent, so the number of solves is capped here.
positions_2 = {f"Splitter i{s}": [(0, 18 - 2 * s), "blue"] for s in range(10)}
positions_2.update({f"Mixer {m}": [(2, 16 - 2 * m), "red"] for m in range(8)})
positions_2.update(
{f"Treatment Unit {t}": [(4, 15 - 2 * t), "green"] for t in range(8)}
)
positions_2.update({f"Splitter {s}": [(6, 16 - 2 * s), "orange"] for s in range(8)})
positions_2["Mixer out"] = [(8, 8), "purple"]
network_2 = Network(
inlet_contaminants=[
[2, 9, 1, 4, 1],
[7, 7, 7, 6, 3],
[1, 7, 0, 6, 6],
[9, 0, 7, 4, 3],
[9, 1, 5, 0, 0],
[0, 8, 0, 6, 3],
[6, 0, 8, 3, 7],
[7, 8, 3, 5, 3],
[3, 7, 4, 0, 6],
[8, 1, 2, 4, 1],
],
removal=[
[0.91, 0.69, 0.7, 0.64, 0.84],
[0.52, 0.92, 0.67, 0.78, 0.73],
[0.9, 0.54, 0.72, 0.71, 0.86],
[0.64, 0.73, 0.66, 0.75, 0.81],
[0.64, 0.86, 0.84, 0.8, 0.84],
[0.68, 0.74, 0.66, 0.72, 0.8],
[0.68, 0.71, 0.68, 0.67, 0.68],
[0.75, 0.64, 0.77, 0.92, 0.82],
],
triplets=[(t, t, t + 10) for t in range(8)],
max_concentration=[1e-4] * 5,
inlet_flow=[8087, 4725, 9336, 7910, 9972, 7826, 6102, 9671, 9923, 5457],
positions=positions_2,
)
draw_network(network_2)routing_2 = minimize_treatment(network_2, max_solves=3)Artelys Knitro 16.0.0: ms_enable=1
ms_maxsolves=3
=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro using up to 8 threads.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 540 variables (44%) and 540 constraints (48%) in 0.00s.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 0.001
findiff_numthreads 1
hessian_no_f 1
hessopt 1
ms_enable 1
ms_maxsolves 3
opttol 1e-06
opttol_abs 0.001
Problem Characteristics | Presolved
-----------------------
Problem type: QCQP
Objective: minimize / linear
Number of variables: 1222 | 682
bounds: lower upper range | lower upper range
1060 0 162 | 527 0 155
free fixed | free fixed
0 0 | 0 0
Number of constraints: 1131 | 591
eq. ineq. range | eq. ineq. range
linear: 694 5 0 | 154 5 0
quadratic: 432 0 0 | 432 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 8 3164 | 8 2084
quadratic: 0 864 432 | 0 864 432
total: 8 4028 432 | 8 2948 432
Knitro parallel multistart will run with 8 threads.
Return codes description
------------------------
0: The final solution satisfies the termination conditions for verifying optimality.
-100 to -199: A feasible approximate solution was found.
-200 to -299: Knitro terminated at an infeasible point.
-300 to -301: The problem was determined to be unbounded.
-400 to -499: Knitro terminated because it reached a pre-defined limit.
-400 to -409: A feasible point was found.
-410 to -419: No feasible point was found.
-500 to -599: Knitro terminated with an input error or some non-standard error.
A more detailed description of individual return codes and their corresponding
termination messages is provided at
https://www.artelys.com/app/docs/knitro/3_referenceManual/returnCodes.html
Solve Thrd Status Objective FeasError Opt Error Solve Time Real Time
----- ---- ------ ------------ ----------- ----------- ----------- -----------
1 5 0 82667.7 3.96885e-05 2.12792e-07 11.9542 11.9697
0 4 0 83194.2 6.83105e-05 1.12066e-07 9.73859 11.9701
2 0 0 83091.0 3.04772e-06 4.16383e-07 39.3045 39.3191
MULTISTART: Best locally optimal solution is returned.
EXIT: All multi-start solves have terminated.
3 solve(s) returned satisfactory solutions.
Final Statistics
----------------
Final objective value = 8.26676604389921e+04
Final feasibility error (abs / rel) = 3.97e-05 / 4.65e-09
Final optimality error (abs / rel) = 2.13e-07 / 2.13e-07
# of iterations = 8502
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 39.32070 ( 66.760 CPU time)
================================================================================
Show the full outputHide the full output
draw_routing(network_2, routing_2)report_output(network_2, routing_2)contaminant discharged limit of limit
(kg/h) (kg/h)
--------------------------------------------------
1 7.90 7.90 100.0%
2 7.90 7.90 100.0%
3 7.90 7.90 100.0%
4 7.90 7.90 100.0%
5 6.19 7.90 78.3%References
- “A multi-commodity flow formulation for the optimal design of wastewater treatment networks”, X. Cheng and X. Li (2020) DOI