struct Network
inlet_contaminants::Matrix{Float64}
removal::Matrix{Float64}
triplets::Vector{Tuple{Int,Int,Int}}
max_concentration::Vector{Float64}
inlet_flow::Vector{Float64}
positions::Vector{Tuple{String,Tuple{Float64,Float64},String}}
end
num_contaminants(network::Network) = length(network.max_concentration)
num_inlet_splitters(network::Network) = length(network.inlet_flow)
num_units(network::Network) = size(network.removal, 1)
function num_splitters(network::Network)
return num_inlet_splitters(network) + num_units(network)
end
num_mixers(network::Network) = num_units(network) + 1
outlet_mixer(network::Network) = num_mixers(network)
inlet_splitters(network::Network) = 1:num_inlet_splitters(network)
function outlet_splitters(network::Network)
return (num_inlet_splitters(network) + 1):num_splitters(network)
end
splitters(network::Network) = 1:num_splitters(network)
mixers(network::Network) = 1:num_mixers(network)
units(network::Network) = 1:num_units(network)
contaminants(network::Network) = 1:num_contaminants(network)
struct Routing
treated_flow::Float64
flow::Dict{Tuple{Int,Int},Float64}
contaminant_flow::Dict{Tuple{Int,Int,Int},Float64}
unit_flow::Dict{Int,Float64}
unit_inflow::Dict{Tuple{Int,Int,Int,Int},Float64}
unit_outflow::Dict{Tuple{Int,Int,Int,Int},Float64}
endWastewater 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(
[3.0; 2.0;;],
[0.9;;],
[(1, 1, 3)],
[1e-4],
[20000.0, 10000.0],
[
("Splitter i0", (0.0, 15.0), "blue"),
("Splitter i1", (0.0, 5.0), "blue"),
("Mixer 0", (2.0, 10.0), "red"),
("Treatment Unit 0", (4.0, 9.0), "green"),
("Splitter 0", (6.0, 10.0), "orange"),
("Mixer out", (8.0, 10.0), "purple"),
],
)
draw_network(network_0)Model implementation
using JuMP
using KNITROEverything the inlet splitters do is affine, since their total flow is data. The outlet splitters are where two variables meet, and the same @constraint macro takes those products as written.
function create_model(network::Network)
splitters, mixers = 1:num_splitters(network), 1:num_mixers(network)
contaminants, units = 1:num_contaminants(network), 1:num_units(network)
triplets = network.triplets
beta, alpha = network.inlet_flow, network.inlet_contaminants
gamma, theta = network.removal, network.max_concentration
m_o = outlet_mixer(network)
model = Model(KNITRO.Optimizer)
@variable(model, fsm[splitters, mixers] >= 0)
@variable(model, fsmc[splitters, mixers, contaminants] >= 0)
@variable(model, ft[units] >= 0)
@variable(model, fitc[triplets, contaminants] >= 0)
@variable(model, fotc[triplets, contaminants] >= 0)
@variable(model, 0 <= xsm[splitters, mixers] <= 1)
@objective(model, Min, sum(ft[t] for t in units))
for s in inlet_splitters(network), m in mixers
x = xsm[s, m]
@constraint(model, beta[s] * x == fsm[s, m])
@constraint(model, [c in contaminants], alpha[s, c] * x == fsmc[s, m, c])
end
for mts in triplets, m in mixers
_, t, s = mts
x = xsm[s, m]
@constraint(model, ft[t] * x == fsm[s, m])
@constraint(model, [c in contaminants], fotc[mts, c] * x == fsmc[s, m, c])
end
@constraint(model, [s in splitters], sum(xsm[s, m] for m in mixers) == 1)
for mts in triplets
mi, t, s = mts
@constraint(model, ft[t] == sum(fsm[s, m] for m in mixers))
@constraint(model, ft[t] == sum(fsm[sp, mi] for sp in splitters))
for c in contaminants
@constraint(model, fotc[mts, c] == sum(fsmc[s, m, c] for m in mixers))
end
for c in contaminants
@constraint(model, fitc[mts, c] == sum(fsmc[sp, mi, c] for sp in splitters))
end
for c in contaminants
@constraint(model, fotc[mts, c] == (1 - gamma[t, c]) * fitc[mts, c])
end
end
for c in contaminants
discharged = sum(fsmc[sp, m_o, c] for sp in splitters)
@constraint(model, discharged <= sum(beta) * theta[c])
end
return model, (fsm=fsm, fsmc=fsmc, ft=ft, fitc=fitc, fotc=fotc)
end
function get_routing(network::Network, model, variables)
splitters, mixers = 1:num_splitters(network), 1:num_mixers(network)
contaminants, units = 1:num_contaminants(network), 1:num_units(network)
triplets = network.triplets
arcs = [(s, m) for s in splitters for m in mixers]
fsm, fsmc, ft, fitc, fotc = variables
flow = Dict(a => value(fsm[a...]) for a in arcs)
load = Dict((s, m, c) => value(fsmc[s, m, c]) for (s, m) in arcs, c in contaminants)
unit_flow = Dict(t => value(ft[t]) for t in units)
inflow = Dict((k..., c) => value(fitc[k, c]) for k in triplets, c in contaminants)
outflow = Dict((k..., c) => value(fotc[k, c]) for k in triplets, c in contaminants)
return Routing(objective_value(model), flow, load, unit_flow, inflow, outflow)
end
function minimize_treatment(network::Network; multistart=true, max_solves=nothing)
model, variables = create_model(network)
set_attribute(model, "ms_enable", Int(multistart))
max_solves === nothing || set_attribute(model, "ms_maxsolves", max_solves)
optimize!(model)
return get_routing(network, model, variables)
endSolving 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.
datacheck 0
feastol 1e-06
feastol_abs 0.001
hessian_no_f 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 5 0 11481.5 7.63100e-08 1.85986e-07 3.36250e-02 6.28828e-02
4 6 0 11481.5 7.63100e-08 1.85986e-07 3.36042e-02 6.28828e-02
3 7 0 11481.5 2.44249e-10 1.00023e-10 3.58307e-02 6.48537e-02
2 2 0 11481.5 2.44249e-10 1.00023e-10 3.58473e-02 6.48537e-02
11 7 0 11481.5 7.37588e-08 8.70487e-10 3.24078e-03 6.48537e-02
14 1 0 11481.5 3.76531e-10 5.12737e-08 8.52801e-03 6.48537e-02
0 0 0 11481.5 9.99800e-11 4.29948e-10 3.51257e-02 6.75725e-02
5 1 0 11481.5 1.00000e-10 1.00000e-10 3.96432e-02 6.75725e-02
15 6 0 11481.5 2.46470e-10 1.00010e-10 3.31004e-03 6.75725e-02
28 2 0 11481.5 1.58333e-10 3.00552e-07 2.72269e-03 6.75725e-02
24 1 0 11481.5 4.23825e-10 1.92595e-10 6.54458e-03 6.75725e-02
17 4 0 11481.5 4.77303e-09 4.35968e-08 4.09121e-03 6.96935e-02
18 5 0 11481.5 2.52728e-10 1.00063e-10 3.19308e-03 6.96935e-02
22 5 0 11481.5 8.63133e-10 2.31317e-07 3.02788e-03 6.96935e-02
23 2 0 11481.5 2.09664e-08 7.63667e-10 3.89265e-03 6.96935e-02
26 0 0 11481.5 1.33636e-06 1.32552e-07 3.40574e-03 6.96935e-02
9 6 0 11481.5 6.51548e-10 4.73840e-09 6.71640e-03 7.15451e-02
10 0 0 11481.5 2.66919e-08 4.52090e-07 6.08771e-03 7.15451e-02
32 1 0 11481.5 1.00000e-10 1.00000e-10 7.73364e-03 7.15451e-02
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.14814808186549e+04
Final feasibility error (abs / rel) = 1.34e-06 / 1.57e-10
Final optimality error (abs / rel) = 1.33e-07 / 1.33e-07
# of iterations = 618
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 0.11575 ( 0.837 CPU time)
================================================================================
Show the full outputHide the full output
draw_routing(network_0, routing_0; shown=[1])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(
[12.0 10.0 10.0; 6.0 3.0 1.5; 1.0 5.0 1.0],
[0.9 0.0 0.0; 0.0 0.99 0.0; 0.0 0.0 0.8],
[(1, 1, 4), (2, 2, 5), (3, 3, 6)],
[1e-4, 1e-4, 1e-4],
[20000.0, 15000.0, 5000.0],
[
("Splitter i0", (0.0, 10.0), "blue"),
("Splitter i1", (0.0, 8.0), "blue"),
("Splitter i2", (0.0, 6.0), "blue"),
("Mixer 0", (2.0, 18.0), "red"),
("Mixer 1", (2.0, 14.0), "red"),
("Mixer 2", (2.0, 1.0), "red"),
("Treatment Unit 0", (4.0, 17.0), "green"),
("Treatment Unit 1", (4.0, 13.0), "green"),
("Treatment Unit 2", (4.0, 0.0), "green"),
("Splitter 0", (6.0, 18.0), "orange"),
("Splitter 1", (6.0, 14.0), "orange"),
("Splitter 2", (6.0, 1.0), "orange"),
("Mixer out", (8.0, 8.0), "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.
datacheck 0
feastol 1e-06
feastol_abs 0.001
hessian_no_f 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 1 0 80778.9 1.60975e-09 2.50203e-09 3.96377e-02 0.450845
1 6 0 80778.9 1.60975e-09 2.50203e-09 3.98747e-02 0.450845
5 3 0 85141.0 1.79063e-07 8.52051e-08 9.24546e-02 0.491857
10 3 0 80778.9 4.10587e-09 9.80479e-09 6.16518e-02 0.491857
11 0 0 92073.0 1.00004e-10 1.00004e-10 5.22828e-02 0.491857
21 7 0 85141.0 2.08820e-09 8.53150e-07 4.24945e-02 0.491857
13 5 0 90268.5 4.30392e-09 8.15743e-10 8.69525e-02 0.499586
6 2 0 87912.8 9.99920e-11 9.99911e-11 0.163100 0.533003
14 7 0 80778.9 9.27685e-10 2.64195e-08 0.107534 0.533003
12 3 0 85141.0 9.12119e-08 1.27145e-07 0.114329 0.542142
16 0 0 80778.9 1.49157e-10 7.77327e-09 8.98407e-02 0.542142
15 2 0 87912.8 1.00004e-10 1.18698e-09 0.126075 0.542142
27 0 0 92073.0 5.63179e-10 2.98969e-09 9.94350e-02 0.542142
23 0 0 90268.5 8.02065e-08 2.39920e-07 0.118812 0.557410
9 6 0 80778.9 7.06481e-07 7.80372e-08 0.151816 0.575558
25 7 0 87912.8 1.00001e-10 1.00001e-10 0.116491 0.575558
0 0 0 85141.0 1.03026e-07 4.06162e-07 0.126530 0.586403
3 5 0 87912.8 8.03845e-09 3.06149e-10 0.160818 0.596831
2 7 0 87912.8 8.03845e-09 3.06149e-10 0.160939 0.596831
29 5 0 81567.0 2.19945e-08 1.79154e-09 0.154033 0.596831
24 2 0 85141.0 5.62487e-08 6.15085e-09 0.144673 0.604662
28 4 0 90268.5 1.31422e-09 1.61549e-09 0.164194 0.604662
31 2 0 80778.9 7.31234e-10 2.74033e-07 0.138097 0.604662
7 4 0 85141.0 1.25146e-09 2.44982e-08 0.216406 0.693744
18 1 0 90268.5 9.99982e-11 4.34284e-09 0.157136 0.693744
26 1 0 90268.5 7.37985e-09 5.40682e-07 0.155626 0.693744
8 1 0 85141.0 5.02041e-10 1.19508e-08 0.166929 0.780710
19 4 0 85141.0 1.73713e-08 4.77724e-07 0.176214 0.786333
20 5 0 90268.5 7.41188e-09 3.28938e-08 0.160099 0.797532
17 6 0 85141.0 1.67347e-10 1.76369e-08 0.248890 0.797532
32 4 0 80778.9 3.20142e-10 4.54055e-09 6.04675e-02 0.820954
34 6 0 92073.0 9.99769e-11 1.42115e-09 9.51900e-02 0.820954
30 7 0 85141.0 5.16957e-09 3.14773e-08 0.234395 0.827852
37 5 0 80778.9 2.29193e-10 1.79430e-08 6.26336e-02 0.841985
36 1 0 92073.0 9.99991e-11 1.00000e-10 7.59329e-02 0.892031
39 7 0 85141.0 1.94996e-09 7.44572e-08 7.63220e-02 0.892031
33 2 0 80778.9 1.00000e-10 1.00000e-10 0.123457 0.926430
41 1 0 81567.0 1.00848e-08 1.59868e-10 9.96395e-02 0.933821
46 4 0 90268.5 1.45064e-07 2.14513e-07 0.115843 0.949090
38 4 0 92073.0 9.99982e-11 9.99973e-11 0.129140 0.964547
45 1 0 92073.0 9.99982e-11 9.99982e-11 0.118703 0.964547
35 0 0 80778.9 1.25086e-08 3.90018e-08 0.181069 0.973296
52 4 0 81567.0 2.21968e-10 5.92615e-09 6.67524e-02 0.994120
48 3 0 81567.0 4.26908e-08 1.12027e-09 0.110687 1.05158
42 6 0 85141.0 1.24055e-09 2.55978e-07 0.167367 1.06944
40 5 0 80778.9 1.00004e-10 1.00004e-10 0.200802 1.07671
43 2 0 80778.9 7.91260e-10 1.00917e-07 0.181032 1.07671
44 7 0 80778.9 1.34311e-08 1.56790e-08 0.204041 1.09609
47 0 0 85141.0 9.99361e-11 1.91049e-08 0.189322 1.12319
49 6 0 80778.9 7.69794e-09 3.02619e-08 0.133983 1.13084
51 1 0 85141.0 1.00055e-10 2.34661e-10 0.158775 1.15202
56 4 0 81567.0 8.90168e-10 5.02200e-08 8.60897e-02 1.16404
53 2 0 90268.5 2.83762e-10 4.59038e-08 0.145299 1.17496
55 3 0 85141.0 5.52973e-10 3.75883e-08 0.117981 1.20727
22 3 -200 375.740 1.00000 1.00000 0.569008 1.23039
57 0 0 87912.8 9.99574e-11 1.83009e-09 0.149079 1.23039
54 7 0 92073.0 5.68980e-09 7.68037e-07 0.178786 1.26574
60 4 0 80778.9 9.99920e-11 4.77686e-09 9.80796e-02 1.26574
59 1 0 80778.9 9.99822e-11 3.55445e-09 0.108011 1.38710
50 5 0 85141.0 5.02562e-07 2.86174e-07 0.227071 1.40854
58 6 0 90268.5 6.00267e-09 9.96602e-09 0.158350 1.41391
67 6 0 85141.0 3.67436e-09 1.25079e-07 7.31021e-02 1.43220
63 5 0 85141.0 1.33407e-10 2.46397e-09 0.111041 1.44004
65 7 0 92073.0 1.50247e-07 5.27726e-08 0.130085 1.45000
64 0 0 81567.0 2.09621e-10 9.99805e-11 0.146912 1.47075
75 6 0 87912.8 1.00004e-10 6.15398e-09 6.48225e-02 1.47075
62 2 0 85141.0 1.28796e-08 2.79554e-08 0.185346 1.55106
77 7 0 80778.9 4.72983e-09 3.23307e-07 0.109896 1.57280
66 4 0 85141.0 1.00013e-10 1.00014e-10 0.179654 1.58305
70 6 0 92073.0 1.00001e-10 1.00002e-10 0.160671 1.58305
76 3 0 80778.9 9.99956e-11 9.99965e-11 0.128172 1.58305
71 7 0 80778.9 9.99139e-11 3.68967e-08 0.151600 1.64841
68 1 0 92073.0 4.55209e-10 1.73932e-10 0.254527 1.64841
73 0 0 80778.9 1.00002e-10 2.65572e-10 0.162557 1.64841
81 5 0 80778.9 1.00012e-10 1.00012e-10 8.93056e-02 1.64841
82 7 0 80778.9 3.75317e-09 5.18224e-09 7.56755e-02 1.64841
78 1 0 80778.9 1.00001e-10 1.00000e-10 0.143692 1.77887
91 6 0 80778.9 7.65340e-10 3.24936e-08 3.42676e-02 1.80879
80 6 0 80778.9 8.73415e-10 9.99982e-11 0.138429 1.81649
69 5 0 85141.0 2.11003e-10 4.89891e-08 0.263125 1.84405
79 0 0 80778.9 9.95417e-11 7.69796e-09 0.195780 1.84405
86 6 0 90268.5 1.00000e-10 1.00000e-10 9.31472e-02 1.84405
93 6 0 85141.0 7.83194e-10 5.38615e-09 6.27388e-02 1.84405
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.
82 solve(s) returned satisfactory solutions.
1 solve(s) converged to infeasible points.
Final Statistics
----------------
Final objective value = 8.07789201484892e+04
Final feasibility error (abs / rel) = 7.06e-07 / 4.54e-11
Final optimality error (abs / rel) = 7.80e-08 / 7.80e-08
# of iterations = 17636
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 1.96639 ( 14.968 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 = [
[("Splitter i$(s)", (0.0, 18.0 - 2s), "blue") for s in 0:9]
[("Mixer $(m)", (2.0, 16.0 - 2m), "red") for m in 0:7]
[("Treatment Unit $(t)", (4.0, 15.0 - 2t), "green") for t in 0:7]
[("Splitter $(s)", (6.0, 16.0 - 2s), "orange") for s in 0:7]
[("Mixer out", (8.0, 8.0), "purple")]
]
network_2 = Network(
[
2.0 9.0 1.0 4.0 1.0;
7.0 7.0 7.0 6.0 3.0;
1.0 7.0 0.0 6.0 6.0;
9.0 0.0 7.0 4.0 3.0;
9.0 1.0 5.0 0.0 0.0;
0.0 8.0 0.0 6.0 3.0;
6.0 0.0 8.0 3.0 7.0;
7.0 8.0 3.0 5.0 3.0;
3.0 7.0 4.0 0.0 6.0;
8.0 1.0 2.0 4.0 1.0
],
[
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
],
[(t, t, t + 10) for t in 1:8],
fill(1e-4, 5),
[8087.0, 4725.0, 9336.0, 7910.0, 9972.0, 7826.0, 6102.0, 9671.0, 9923.0, 5457.0],
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.
datacheck 0
feastol 1e-06
feastol_abs 0.001
hessian_no_f 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
----- ---- ------ ------------ ----------- ----------- ----------- -----------
2 1 0 83194.2 3.45214e-05 7.03156e-07 11.0075 14.2909
1 2 0 83194.2 1.19045e-05 9.43883e-09 16.0658 16.0769
0 0 0 83194.2 8.92644e-05 6.37347e-08 29.6417 29.6521
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.31942264087918e+04
Final feasibility error (abs / rel) = 8.93e-05 / 1.03e-08
Final optimality error (abs / rel) = 6.37e-08 / 6.37e-08
# of iterations = 7158
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 29.65248 ( 56.754 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.60 7.90 70.9%References
- “A multi-commodity flow formulation for the optimal design of wastewater treatment networks”, X. Cheng and X. Li (2020) DOI