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 knitroEverything the inlet splitters do is affine, since their total flow is data. The outlet splitters are where two variables meet, and add_constraint takes those products 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
arcs = [(s, m) for s in splitters for m in mixers]
arc_loads = [(s, m, c) for s, m in arcs for c in contaminants]
unit_loads = [(*k, c) for k in triplets for c in contaminants]
prob = knitro.Problem()
fsm = {a: prob.add_variable(lb=0) for a in arcs}
fsmc = {k: prob.add_variable(lb=0) for k in arc_loads}
ft = {t: prob.add_variable(lb=0) for t in units}
fitc = {k: prob.add_variable(lb=0) for k in unit_loads}
fotc = {k: prob.add_variable(lb=0) for k in unit_loads}
xsm = {a: prob.add_variable(lb=0, ub=1) for a in arcs}
prob.add_objective(prob.nsum(ft[t] for t in units))
for s in inlets:
for m in mixers:
prob.add_constraint(beta[s] * xsm[s, m] == fsm[s, m])
for c in contaminants:
prob.add_constraint(alpha[s][c] * xsm[s, m] == fsmc[s, m, c])
for mi, t, s in triplets:
for m1 in mixers:
prob.add_constraint(ft[t] * xsm[s, m1] == fsm[s, m1])
for c in contaminants:
prob.add_constraint(fotc[mi, t, s, c] * xsm[s, m1] == fsmc[s, m1, c])
for s in splitters:
prob.add_constraint(prob.nsum(xsm[s, m] for m in mixers) == 1)
for mi, t, s in triplets:
prob.add_constraint(ft[t] == prob.nsum(fsm[s, m] for m in mixers))
prob.add_constraint(ft[t] == prob.nsum(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]
prob.add_constraint(f_out == prob.nsum(fsmc[s, m, c] for m in mixers))
prob.add_constraint(f_in == prob.nsum(fsmc[sp, mi, c] for sp in splitters))
prob.add_constraint(f_out == (1 - gamma[t][c]) * f_in)
for c in contaminants:
discharged = prob.nsum(fsmc[sp, m_o, c] for sp in splitters)
prob.add_constraint(discharged <= sum(beta[s] * theta[c] for s in inlets))
prob.set_param(knitro.KN_PARAM_MS_ENABLE, int(multistart))
if max_solves is not None:
prob.set_param(knitro.KN_PARAM_MS_MAXSOLVES, max_solves)
prob.solve()
treated_flow = prob.get_attr(knitro.KN_ATTR_OBJ_VALUE)
flows = [{k: v.value 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)=======================================
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
feastol 1e-06
feastol_abs 0.001
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 expressions: 64 | 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
----- ---- ------ ------------ ----------- ----------- ----------- -----------
4 7 0 11481.5 1.97906e-09 3.43297e-08 0.129678 0.180313
1 3 0 11481.5 1.97906e-09 3.43297e-08 0.130028 0.180313
3 1 0 11481.5 9.99969e-11 9.99967e-11 0.132175 0.187040
2 2 0 11481.5 9.99969e-11 9.99967e-11 0.132522 0.187040
14 3 0 11481.5 9.10831e-11 5.49294e-08 2.83785e-03 0.187040
13 5 0 11481.5 1.00014e-10 1.00014e-10 4.43251e-03 0.190530
7 4 0 11481.5 1.82482e-10 1.00010e-10 0.137887 0.192095
9 3 0 11481.5 9.99973e-11 9.99974e-11 1.12145e-02 0.192095
18 2 0 11481.5 6.21018e-08 9.08553e-07 4.18704e-03 0.192095
25 3 0 11481.5 9.99765e-11 4.45056e-10 5.41833e-03 0.192095
33 3 0 11481.5 6.36739e-09 9.62100e-08 3.75472e-03 0.192095
23 0 0 11481.5 1.00019e-10 1.00019e-10 5.75907e-03 0.192442
22 4 0 11481.5 5.15865e-09 1.05347e-08 6.83681e-03 0.192442
32 0 0 11481.5 1.43555e-07 5.12142e-07 5.89243e-03 0.192442
6 5 0 11481.5 6.22277e-10 4.11470e-07 0.140828 0.192544
10 1 0 11481.5 6.49711e-09 1.06138e-09 1.10832e-02 0.192544
11 2 0 11481.5 1.05138e-09 4.75666e-10 1.41867e-02 0.192671
19 7 0 11481.5 2.63753e-09 7.55205e-07 8.83170e-03 0.192671
30 3 0 11481.5 5.04510e-08 3.00784e-08 1.15812e-02 0.192671
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.
19 solve(s) returned satisfactory solutions.
Final Statistics
----------------
Final objective value = 1.14814814659550e+04
Final feasibility error (abs / rel) = 5.16e-09 / 6.36e-13
Final optimality error (abs / rel) = 1.05e-08 / 1.05e-08
# of iterations = 800
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 0.21373 ( 1.483 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)=======================================
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
feastol 1e-06
feastol_abs 0.001
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 expressions: 404 | 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 87912.8 1.41960e-06 1.01201e-07 2.21374e-02 0.705323
1 2 0 87912.8 1.41960e-06 1.01201e-07 2.20884e-02 0.705323
27 0 0 90268.5 1.00001e-10 1.00000e-10 6.44518e-02 0.774819
8 3 0 85141.0 1.00004e-10 1.00004e-10 7.16249e-02 0.797981
6 6 0 80778.9 1.08634e-07 2.09926e-07 0.135253 0.821372
28 1 0 80778.9 1.00022e-10 1.00022e-10 0.130093 0.832700
10 3 0 81567.0 5.67525e-10 1.71998e-08 0.120261 0.846738
11 6 0 85141.0 7.13044e-10 9.27011e-08 0.153385 0.854518
23 1 0 80778.9 4.14730e-10 6.09038e-08 0.128373 0.854518
24 2 0 90268.5 1.08230e-10 6.30852e-09 0.136911 0.854518
21 7 0 87912.8 1.00081e-10 2.70703e-09 0.141494 0.863833
7 4 0 90268.5 3.18107e-08 4.44850e-07 0.193178 0.872264
17 2 0 80778.9 7.96974e-07 4.77736e-07 0.147476 0.872264
26 7 0 81567.0 9.03316e-09 4.21987e-10 0.139394 0.872264
29 2 0 85141.0 5.70290e-09 1.68179e-08 0.140855 0.872264
22 4 0 87912.8 6.75943e-10 1.00155e-10 0.188208 0.882425
0 5 0 80778.9 4.84565e-08 5.44635e-08 0.175512 0.893583
15 7 0 80778.9 1.00012e-10 1.00011e-10 0.251494 0.893583
31 3 0 85141.0 2.80852e-09 2.49094e-10 0.167362 0.893583
16 6 0 80778.9 2.19141e-09 3.37137e-09 0.221505 0.901707
20 3 0 80778.9 1.00011e-10 1.00010e-10 0.214607 0.994241
9 2 0 85141.0 1.52795e-10 5.23607e-09 0.330787 1.00110
14 3 0 80778.9 4.34165e-09 1.37050e-09 0.250840 1.00110
13 4 0 81567.0 1.31242e-10 3.18161e-10 0.291471 1.00110
5 7 0 87912.8 7.65431e-09 3.59132e-07 0.214527 1.01123
18 0 0 85141.0 5.07536e-09 1.65140e-07 0.246679 1.01123
25 6 0 87912.8 5.02041e-10 9.03387e-08 0.275960 1.13884
12 5 0 90268.5 1.00000e-10 1.44236e-10 0.279081 1.14954
36 5 0 80778.9 2.01287e-08 3.51702e-09 8.93135e-02 1.15820
38 6 0 87912.8 1.00008e-10 1.00009e-10 9.94493e-02 1.17625
19 5 0 85141.0 1.78110e-07 1.15965e-08 0.271406 1.21506
33 7 0 80778.9 9.99485e-11 4.54702e-09 0.129077 1.21506
2 0 0 81567.0 1.88191e-08 2.65984e-10 0.374525 1.37058
3 1 0 81567.0 1.88191e-08 2.65984e-10 0.500166 1.37058
32 2 0 80778.9 5.83809e-09 9.72856e-07 0.165634 1.37058
39 5 0 90268.5 2.06972e-10 6.94377e-10 0.156106 1.40691
50 7 0 81567.0 8.39320e-08 1.16412e-08 5.97782e-02 1.40691
37 1 0 90268.5 1.00001e-10 1.00001e-10 0.248550 1.40718
46 6 0 81567.0 9.99734e-11 2.77607e-10 8.22596e-02 1.40718
48 6 0 90268.5 2.03727e-10 8.52015e-09 6.52714e-02 1.40718
34 0 0 81567.0 1.00000e-10 1.00000e-10 0.228830 1.64302
40 6 0 90268.5 1.92676e-08 4.69240e-07 0.201725 1.64302
30 4 -202 95219.0 0.850000 1.00000 0.564262 1.66153
42 2 0 80778.9 9.99991e-11 1.35692e-09 0.267536 1.67931
43 0 0 87912.8 5.12955e-10 6.39197e-08 0.244106 1.70266
45 1 0 80778.9 7.57508e-10 5.27721e-08 0.204724 1.72094
51 7 0 90268.5 4.09760e-09 8.26611e-09 0.124157 1.72094
41 7 0 85141.0 9.99991e-11 9.99982e-11 0.323641 1.73992
60 5 0 81567.0 1.91411e-08 3.46307e-07 8.74286e-02 1.76135
53 6 0 87912.8 2.04153e-08 5.31911e-09 0.141196 1.79682
59 3 0 90268.5 9.97922e-11 2.65179e-08 0.172015 1.80920
56 6 0 85141.0 8.69797e-08 7.89940e-07 0.119324 1.81668
57 0 0 80778.9 1.20977e-07 9.94228e-08 0.143392 1.81668
52 0 0 87912.8 8.98980e-09 3.52932e-08 0.155419 1.85155
54 1 0 80778.9 1.00022e-09 6.89295e-08 0.190248 2.00010
55 7 0 80778.9 1.00038e-10 1.54621e-10 0.166336 2.00921
58 1 0 85141.0 6.23030e-09 4.54816e-07 0.253473 2.05722
49 2 0 80778.9 1.00002e-10 1.00001e-10 0.298056 2.06756
47 4 0 80778.9 7.12134e-10 5.71748e-08 0.410043 2.07337
68 1 0 80778.9 1.00001e-10 1.00001e-10 0.129791 2.17088
44 5 0 90268.5 8.39998e-09 3.47814e-07 0.464181 2.18986
61 2 0 80778.9 1.00000e-10 1.00000e-10 0.234511 2.18986
69 5 0 80778.9 1.24335e-07 1.24932e-07 0.139720 2.18986
63 5 0 85141.0 4.00178e-10 9.49102e-09 0.168755 2.19560
62 4 0 80778.9 5.21424e-08 5.73732e-10 0.211760 2.19560
35 3 -200 5619.37 1.00000 1.00000 0.714198 2.28684
71 2 0 87912.8 1.00030e-10 7.74027e-10 0.218114 2.31843
70 4 0 92073.0 2.14161e-07 9.20440e-10 0.229333 2.34079
64 7 0 85141.0 1.13505e-09 3.90294e-10 0.270782 2.36763
66 0 0 85141.0 1.00000e-10 1.00000e-10 0.266173 2.36763
77 2 0 80778.9 4.47414e-07 1.95391e-08 9.00776e-02 2.37883
73 1 0 85141.0 1.01923e-07 1.02514e-07 0.157972 2.44127
72 5 0 80778.9 1.00004e-10 1.00005e-10 0.147815 2.44892
76 4 0 85141.0 9.99494e-11 1.04492e-07 0.159564 2.46252
75 7 0 92073.0 3.81813e-10 9.88263e-10 0.186273 2.47523
84 0 0 92073.0 4.71810e-09 2.03227e-09 0.153876 2.48560
79 5 0 80778.9 1.30849e-08 3.67806e-07 0.175931 2.49501
81 7 0 80778.9 4.55293e-09 2.83804e-07 0.185632 2.49501
85 4 0 92073.0 1.28283e-08 6.02099e-08 7.87998e-02 2.49501
86 5 0 81567.0 9.99982e-11 9.99982e-11 8.17760e-02 2.57589
74 0 0 81567.0 9.98925e-11 5.88317e-10 0.217758 2.58796
82 6 0 80778.9 1.07102e-08 6.85128e-08 0.171354 2.58796
67 3 0 85141.0 3.80591e-09 8.27698e-07 0.380896 2.60255
87 3 0 80778.9 2.95108e-09 8.16885e-08 8.65301e-02 2.62537
83 2 0 87912.8 9.99965e-11 2.13277e-10 0.178562 2.63606
65 6 0 85141.0 2.42858e-09 4.79786e-07 0.535562 2.64553
91 5 0 92073.0 1.01944e-08 7.62865e-09 7.34007e-02 2.64553
78 3 0 85141.0 1.00019e-10 1.00019e-10 0.201644 2.65636
94 7 0 80778.9 5.68114e-10 1.60420e-08 8.04506e-02 2.65636
80 1 0 92073.0 8.09558e-08 9.16096e-08 0.220517 2.67121
90 2 0 85141.0 2.43708e-08 5.26357e-07 9.99208e-02 2.68077
95 1 0 80778.9 1.00017e-10 1.00017e-10 8.63292e-02 2.69301
92 3 0 80778.9 9.99245e-11 8.43465e-09 0.115067 2.69301
99 1 0 80778.9 2.94676e-10 3.06161e-08 0.127851 2.76304
97 2 0 92073.0 2.02999e-09 2.80550e-08 0.187287 2.76304
98 7 0 80778.9 9.99894e-11 7.58079e-10 0.196273 2.79922
96 5 0 85141.0 1.00000e-10 1.00000e-10 0.233732 2.79936
89 6 0 85141.0 5.99670e-10 1.36885e-08 0.427916 2.89081
93 4 0 85141.0 1.00012e-10 1.00012e-10 0.200016 2.89094
88 0 0 85141.0 1.00002e-10 1.98997e-09 0.306019 2.89135
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.
98 solve(s) returned satisfactory solutions.
1 solve(s) converged to infeasible points.
1 solve(s) reached no improvement limit.
Final Statistics
----------------
Final objective value = 8.07789201479966e+04
Final feasibility error (abs / rel) = 7.97e-07 / 4.57e-11
Final optimality error (abs / rel) = 4.78e-07 / 4.78e-07
# of iterations = 25879
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 2.89208 ( 21.053 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)=======================================
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
feastol 1e-06
feastol_abs 0.001
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 expressions: 3343 | 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 7 0 83194.2 3.40085e-05 2.02630e-07 6.61499 6.62665
2 5 0 83091.0 5.05818e-05 1.11795e-08 9.69882 9.71052
0 0 0 83194.2 2.20664e-06 1.84475e-08 9.84660 9.85797
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.30910025896985e+04
Final feasibility error (abs / rel) = 5.06e-05 / 6.04e-09
Final optimality error (abs / rel) = 1.12e-08 / 1.12e-08
# of iterations = 3195
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 9.85828 ( 26.194 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 5.35 7.90 67.8%References
- “A multi-commodity flow formulation for the optimal design of wastewater treatment networks”, X. Cheng and X. Li (2020) DOI