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 8.65072e-02 0.117444
4 4 0 11481.5 7.63100e-08 1.85986e-07 8.69289e-02 0.117444
2 1 0 11481.5 2.44249e-10 1.00023e-10 8.86184e-02 0.119498
3 3 0 11481.5 2.44249e-10 1.00023e-10 8.87454e-02 0.119498
11 1 0 11481.5 7.37588e-08 8.70487e-10 3.22739e-03 0.119498
14 1 0 11481.5 3.76531e-10 5.12737e-08 3.04266e-03 0.119498
0 0 0 11481.5 9.99800e-11 4.29948e-10 8.82218e-02 0.121470
5 7 0 11481.5 1.00000e-10 1.00000e-10 8.85829e-02 0.121470
15 4 0 11481.5 2.46470e-10 1.00010e-10 3.31832e-03 0.121470
24 1 0 11481.5 4.23825e-10 1.92595e-10 3.50946e-03 0.121470
28 0 0 11481.5 1.58333e-10 3.00552e-07 2.53171e-03 0.121470
18 6 0 11481.5 2.52728e-10 1.00063e-10 3.34189e-03 0.122843
17 1 0 11481.5 4.77303e-09 4.35968e-08 5.18020e-03 0.122843
22 6 0 11481.5 8.63133e-10 2.31317e-07 3.13684e-03 0.122843
23 3 0 11481.5 2.09664e-08 7.63667e-10 3.75414e-03 0.122843
26 7 0 11481.5 1.33636e-06 1.32552e-07 3.38040e-03 0.122843
9 4 0 11481.5 6.51548e-10 4.73840e-09 7.71682e-03 0.172917
10 0 0 11481.5 2.66919e-08 4.52090e-07 6.66217e-03 0.172917
32 5 0 11481.5 1.00000e-10 1.00000e-10 3.88057e-03 0.172917
31 0 0 11481.5 6.62881e-09 3.29025e-10 1.99900e-02 0.172917
30 7 0 11481.5 1.00035e-10 1.04823e-10 2.30386e-02 0.172917
33 2 0 11481.5 1.06821e-10 2.30778e-08 5.38477e-02 0.172917
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.
22 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 = 1310
# of CG iterations = 0
# of function evaluations = 0
# of gradient evaluations = 0
# of Hessian evaluations = 0
Total program time (secs) = 0.17327 ( 1.226 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 4 0 80778.9 1.60975e-09 2.50203e-09 2.86115e-02 0.447264
1 6 0 80778.9 1.60975e-09 2.50203e-09 2.97026e-02 0.447264
5 3 0 85141.0 1.79063e-07 8.52051e-08 8.85839e-02 0.485505
10 3 0 80778.9 4.10587e-09 9.80479e-09 6.02945e-02 0.485505
11 0 0 92073.0 1.00004e-10 1.00004e-10 5.26793e-02 0.485505
21 2 0 85141.0 2.08820e-09 8.53150e-07 4.29649e-02 0.485505
13 3 0 90268.5 4.30392e-09 8.15743e-10 8.78469e-02 0.491488
6 2 0 87912.8 9.99920e-11 9.99911e-11 0.130140 0.525177
14 7 0 80778.9 9.27685e-10 2.64195e-08 0.109345 0.525177
12 2 0 85141.0 9.12119e-08 1.27145e-07 0.115145 0.534560
16 5 0 80778.9 1.49157e-10 7.77327e-09 9.12983e-02 0.534560
15 1 0 87912.8 1.00004e-10 1.18698e-09 0.128431 0.534560
27 7 0 92073.0 5.63179e-10 2.98969e-09 0.100811 0.534560
23 7 0 90268.5 8.02065e-08 2.39920e-07 0.120903 0.550156
9 6 0 80778.9 7.06481e-07 7.80372e-08 0.160830 0.568852
25 1 0 87912.8 1.00001e-10 1.00001e-10 0.120179 0.568852
0 0 0 85141.0 1.03026e-07 4.06162e-07 0.127535 0.579870
2 1 0 87912.8 8.03845e-09 3.06149e-10 0.167731 0.591025
3 5 0 87912.8 8.03845e-09 3.06149e-10 0.168242 0.591025
29 4 0 81567.0 2.19945e-08 1.79154e-09 0.175595 0.591025
24 2 0 85141.0 5.62487e-08 6.15085e-09 0.149631 0.598945
31 0 0 80778.9 7.31234e-10 2.74033e-07 0.138975 0.598945
28 3 0 90268.5 1.31422e-09 1.61549e-09 0.168084 0.598945
7 7 0 85141.0 1.25146e-09 2.44982e-08 0.167232 0.688833
18 6 0 90268.5 9.99982e-11 4.34284e-09 0.157833 0.688833
26 6 0 90268.5 7.37985e-09 5.40682e-07 0.159399 0.688833
8 4 0 85141.0 5.02041e-10 1.19508e-08 0.175304 0.771663
19 4 0 85141.0 1.73713e-08 4.77724e-07 0.198982 0.779116
20 3 0 90268.5 7.41188e-09 3.28938e-08 0.164459 0.790243
17 0 0 85141.0 1.67347e-10 1.76369e-08 0.249529 0.790243
32 6 0 80778.9 3.20142e-10 4.54055e-09 6.10759e-02 0.813894
34 2 0 92073.0 9.99769e-11 1.42115e-09 9.63600e-02 0.813894
30 1 0 85141.0 5.16957e-09 3.14773e-08 0.241114 0.820779
37 0 0 80778.9 2.29193e-10 1.79430e-08 6.34717e-02 0.849921
36 4 0 92073.0 9.99991e-11 1.00000e-10 7.79767e-02 0.878101
39 1 0 85141.0 1.94996e-09 7.44572e-08 7.82925e-02 0.878101
33 3 0 80778.9 1.00000e-10 1.00000e-10 0.125053 0.913160
41 4 0 81567.0 1.00848e-08 1.59868e-10 0.101643 0.920525
46 6 0 90268.5 1.45064e-07 2.14513e-07 0.115728 0.953023
38 6 0 92073.0 9.99982e-11 9.99973e-11 0.128016 0.966833
45 7 0 92073.0 9.99982e-11 9.99982e-11 0.119658 0.966833
35 7 0 80778.9 1.25086e-08 3.90018e-08 0.166787 0.975349
52 6 0 81567.0 2.21968e-10 5.92615e-09 6.65710e-02 0.993261
48 5 0 81567.0 4.26908e-08 1.12027e-09 0.111572 1.05322
42 2 0 85141.0 1.24055e-09 2.55978e-07 0.170410 1.06785
40 0 0 80778.9 1.00004e-10 1.00004e-10 0.201671 1.07454
43 3 0 80778.9 7.91260e-10 1.00917e-07 0.184182 1.07454
44 1 0 80778.9 1.34311e-08 1.56790e-08 0.217510 1.09198
47 4 0 85141.0 9.99361e-11 1.91049e-08 0.189167 1.12390
49 2 0 80778.9 7.69794e-09 3.02619e-08 0.134611 1.13150
51 7 0 85141.0 1.00055e-10 2.34661e-10 0.158958 1.14777
56 5 0 81567.0 8.90168e-10 5.02200e-08 8.56202e-02 1.16174
53 3 0 90268.5 2.83762e-10 4.59038e-08 0.146120 1.17245
55 1 0 85141.0 5.52973e-10 3.75883e-08 0.120202 1.20259
22 5 -200 375.740 1.00000 1.00000 0.586701 1.22717
22 5 * 32170.4 0.998949
57 4 0 87912.8 9.99574e-11 1.83009e-09 0.147948 1.22717
54 6 0 92073.0 5.68980e-09 7.68037e-07 0.178394 1.26378
60 5 0 80778.9 9.99920e-11 4.77686e-09 9.74064e-02 1.26378
59 7 0 80778.9 9.99822e-11 3.55445e-09 0.109656 1.38717
50 0 0 85141.0 5.02562e-07 2.86174e-07 0.226235 1.41358
58 2 0 90268.5 6.00267e-09 9.96602e-09 0.159256 1.42095
67 5 0 85141.0 3.67436e-09 1.25079e-07 7.38068e-02 1.43961
63 0 0 85141.0 1.33407e-10 2.46397e-09 0.110616 1.44744
65 6 0 92073.0 1.50247e-07 5.27726e-08 0.129178 1.45742
64 4 0 81567.0 2.09621e-10 9.99805e-11 0.145822 1.47174
75 5 0 87912.8 1.00004e-10 6.15398e-09 6.57326e-02 1.47174
62 1 0 85141.0 1.28796e-08 2.79554e-08 0.185815 1.55938
77 6 0 80778.9 4.72983e-09 3.23307e-07 0.111093 1.58148
66 7 0 85141.0 1.00013e-10 1.00014e-10 0.198778 1.59178
70 5 0 92073.0 1.00001e-10 1.00002e-10 0.161319 1.59178
76 3 0 80778.9 9.99956e-11 9.99965e-11 0.130850 1.59178
71 6 0 80778.9 9.99139e-11 3.68967e-08 0.152655 1.65333
68 2 0 92073.0 4.55209e-10 1.73932e-10 0.257446 1.65333
73 4 0 80778.9 1.00002e-10 2.65572e-10 0.162768 1.65333
81 0 0 80778.9 1.00012e-10 1.00012e-10 8.87893e-02 1.65333
82 6 0 80778.9 3.75317e-09 5.18224e-09 7.66362e-02 1.65333
78 2 0 80778.9 1.00001e-10 1.00000e-10 0.144902 1.79214
91 7 0 80778.9 7.65340e-10 3.24936e-08 3.54388e-02 1.82188
80 5 0 80778.9 8.73415e-10 9.99982e-11 0.139012 1.82961
69 0 0 85141.0 2.11003e-10 4.89891e-08 0.262943 1.85743
79 4 0 80778.9 9.95417e-11 7.69796e-09 0.194742 1.85743
86 5 0 90268.5 1.00000e-10 1.00000e-10 9.39708e-02 1.85743
93 7 0 85141.0 7.83194e-10 5.38615e-09 6.30501e-02 1.85743
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.98166 ( 15.193 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 7 0 83194.2 3.45214e-05 7.03156e-07 11.3509 15.1969
1 5 0 83194.2 1.19045e-05 9.43883e-09 17.0528 17.0630
0 0 0 83194.2 8.92644e-05 6.37347e-08 30.5921 30.6023
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) = 30.60266 ( 59.046 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