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
----- ---- ------ ------------ ----------- ----------- ----------- -----------
4 2 0 11481.5 1.00004e-10 1.00004e-10 0.148996 0.225867
1 4 0 11481.5 1.00004e-10 1.00004e-10 0.183721 0.225867
3 1 0 11481.5 8.34888e-10 1.00048e-10 0.158585 0.230613
2 3 0 11481.5 8.34888e-10 1.00048e-10 0.173229 0.230613
11 0 0 11481.5 9.99973e-11 9.99973e-11 3.61979e-03 0.233273
21 1 0 11481.5 1.00303e-10 6.16849e-09 2.26149e-03 0.233273
31 6 0 11481.5 1.00000e-10 1.00000e-10 2.96758e-03 0.233273
8 2 0 11481.5 4.32158e-08 6.09998e-09 2.74566e-02 0.235741
13 4 0 11481.5 9.99854e-11 9.58607e-09 1.57123e-02 0.235741
28 6 0 11481.5 1.00000e-10 1.00000e-10 3.68924e-03 0.235741
10 3 0 11481.5 1.01486e-10 2.03180e-08 1.16563e-02 0.237040
18 7 0 11481.5 5.38730e-08 1.84683e-07 1.44777e-02 0.237040
22 5 0 11481.5 5.83190e-09 5.94482e-08 5.01289e-03 0.237040
24 7 0 11481.5 9.98317e-11 6.27551e-09 7.69471e-03 0.237040
29 5 0 11481.5 1.10084e-06 3.05377e-07 6.69253e-03 0.237040
32 6 0 11481.5 2.07516e-10 1.00014e-10 3.72743e-03 0.237040
33 5 0 11481.5 3.26850e-07 6.07884e-08 1.37972e-02 0.244762
6 7 0 11481.5 9.91847e-11 2.45321e-09 0.189235 0.245087
25 5 0 11481.5 2.83990e-07 2.44632e-08 9.10900e-03 0.245087
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.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 = 782
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 0.26789 ( 4.437 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
----- ---- ------ ------------ ----------- ----------- ----------- -----------
1 4 0 92073.0 3.65253e-09 1.25108e-07 0.132902 0.750686
4 3 0 92073.0 3.65253e-09 1.25108e-07 0.118490 0.750686
9 3 0 85141.0 1.00072e-10 1.00072e-10 9.35853e-02 0.843695
24 1 0 90268.5 3.31440e-09 2.38763e-08 8.68080e-02 0.843695
17 4 0 81567.0 1.00002e-10 1.00002e-10 7.02627e-02 0.863454
18 1 0 90268.5 2.69574e-09 6.98212e-07 0.104894 0.882711
12 2 0 90268.5 1.00017e-10 1.00017e-10 0.122839 0.899489
22 0 0 80778.9 1.27991e-08 4.54936e-07 7.85945e-02 0.899489
11 1 0 80778.9 1.00004e-10 3.67229e-08 0.124490 0.908843
27 6 0 80778.9 8.13725e-09 4.89626e-09 0.163118 0.908843
19 2 0 81567.0 1.19348e-09 4.61667e-09 0.121718 0.918132
0 2 0 85141.0 3.80828e-10 1.08412e-09 0.336746 0.927838
28 1 0 90268.5 8.10142e-08 3.38638e-07 0.155748 0.927838
26 2 0 87912.8 2.12090e-10 3.99564e-09 0.161552 0.938856
30 3 0 92073.0 9.73960e-08 4.08422e-08 0.160745 0.938856
14 0 0 92073.0 1.00005e-10 1.00005e-10 0.144377 0.948611
25 0 0 85141.0 2.71472e-08 1.74499e-08 0.163867 0.948611
6 1 0 80778.9 1.16727e-08 1.60799e-09 0.297335 0.959390
7 5 0 85141.0 1.00000e-10 1.00000e-10 0.317569 0.959390
13 5 0 92073.0 1.00001e-10 1.00001e-10 0.197127 0.959390
29 7 0 85141.0 1.23691e-10 5.83477e-08 0.198921 0.984726
31 4 0 81567.0 1.00049e-10 7.55953e-10 0.242160 0.984726
5 7 0 80778.9 9.99529e-11 6.99401e-10 0.337205 1.07771
10 3 0 85141.0 2.40470e-09 2.77842e-08 0.232763 1.08907
21 4 0 92073.0 5.63887e-10 3.69955e-08 0.243044 1.08907
2 0 0 80778.9 2.57347e-08 7.59701e-09 0.345359 1.10079
3 6 0 80778.9 2.57347e-08 7.59701e-09 0.363153 1.10079
20 3 0 80778.9 1.96451e-10 2.00725e-09 0.232378 1.15101
15 7 0 92073.0 1.69838e-07 1.13552e-08 0.320629 1.22603
16 6 0 90268.5 3.03066e-10 1.33868e-10 0.250007 1.23607
8 4 0 80778.9 9.01195e-08 5.92734e-08 0.282950 1.25033
35 6 0 85141.0 6.93159e-08 6.45247e-09 0.120928 1.26176
23 5 0 85141.0 2.62208e-09 3.37035e-07 0.413426 1.27595
40 5 0 85141.0 1.49062e-07 8.54941e-09 0.109748 1.27595
36 5 0 80778.9 4.41047e-07 1.39258e-07 0.117136 1.36801
39 1 0 81567.0 1.48176e-07 2.76558e-09 0.163682 1.36801
42 1 0 80778.9 1.00000e-10 1.00000e-10 0.109311 1.37957
38 4 0 80778.9 3.49246e-10 4.05584e-08 0.178211 1.39789
32 3 0 92073.0 1.70420e-08 2.81452e-08 0.299416 1.41242
48 1 0 80778.9 3.89264e-10 1.90539e-10 7.49612e-02 1.47193
45 0 0 81567.0 2.98314e-10 1.82385e-10 0.157861 1.49131
44 5 0 80778.9 1.00009e-10 1.00009e-10 0.178687 1.49131
46 2 0 81567.0 7.89574e-10 1.00009e-10 0.166131 1.49131
33 0 0 85141.0 1.00001e-10 1.00001e-10 0.237446 1.55610
49 7 0 90268.5 9.99982e-11 1.19221e-09 0.129124 1.72793
51 0 0 85141.0 9.99973e-11 9.99973e-11 0.139784 1.74011
34 2 0 81567.0 1.18884e-06 7.90391e-09 0.246330 1.75356
37 7 0 80778.9 9.99938e-11 2.32535e-09 0.357936 1.77056
52 5 0 85141.0 5.89443e-10 4.69924e-10 0.177693 1.77056
50 1 0 81567.0 5.63716e-10 3.36818e-08 0.222979 1.78077
47 3 0 80778.9 2.42478e-09 1.09449e-09 0.182538 1.79708
59 1 0 85141.0 6.40284e-10 1.44080e-08 0.141089 1.92244
53 2 0 85141.0 7.54983e-10 2.35559e-10 0.251398 1.93091
57 5 0 85141.0 1.55391e-10 2.03354e-09 0.137919 1.93091
60 4 0 81567.0 6.37854e-09 1.00082e-10 0.175447 1.93091
43 4 0 81567.0 1.78261e-08 1.28455e-09 0.442281 2.06986
62 5 0 80778.9 9.99938e-11 6.01631e-10 0.149706 2.06986
66 7 0 85141.0 7.53062e-10 2.16645e-07 9.70475e-02 2.10687
41 6 0 81567.0 5.96013e-10 5.03365e-07 0.453930 2.16657
70 0 0 92073.0 1.00012e-10 1.00012e-10 8.72472e-02 2.16657
54 3 0 90268.5 9.99965e-11 9.34525e-10 0.328808 2.19215
55 7 0 80778.9 3.97306e-07 2.98519e-07 0.317036 2.21902
63 1 0 80778.9 1.43442e-09 2.14801e-09 0.394437 2.27614
58 6 0 80778.9 1.44678e-09 3.65376e-08 0.313036 2.28762
65 3 0 80778.9 7.79180e-09 2.61610e-10 0.284257 2.28762
68 6 0 80778.9 6.45377e-09 4.56412e-09 0.252404 2.30269
64 4 0 92073.0 9.99973e-11 9.99973e-11 0.313628 2.34992
72 3 0 80778.9 7.05768e-10 1.16163e-07 0.194895 2.34992
75 1 0 85141.0 7.85304e-08 4.29072e-08 0.159337 2.41464
61 2 0 85141.0 3.13776e-10 6.38101e-09 0.397935 2.44349
67 5 0 85141.0 1.00009e-10 1.00009e-10 0.311499 2.44349
71 2 0 85141.0 4.47471e-10 7.26031e-08 0.206724 2.45357
79 1 0 85141.0 2.70679e-07 7.53030e-07 5.92565e-02 2.48157
78 2 0 80778.9 9.99947e-11 9.99947e-11 0.119424 2.48157
69 7 0 92073.0 2.13541e-08 2.13186e-09 0.315113 2.59002
82 1 0 87912.8 1.43700e-09 4.44834e-07 9.33625e-02 2.60033
80 7 0 80778.9 1.00001e-10 1.00001e-10 0.105820 2.60891
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) = 2.82638 ( 24.360 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 13.7419 13.7741
0 6 0 83194.2 6.83105e-05 1.12066e-07 11.6611 13.7745
2 7 0 83091.0 3.04772e-06 4.16383e-07 45.5676 45.5911
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) = 45.59324 ( 76.363 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