const INCH = 0.0254
struct Node
demand::Float64
elevation::Float64
min_pressure::Float64
max_pressure::Float64
kind::Symbol
end
struct Arc
from_node::Int
to_node::Int
length::Float64
min_diameter::Float64
max_diameter::Float64
max_velocity::Float64
end
struct Network
nodes::Vector{Node}
pipes::Vector{Arc}
diameters::Vector{Float64}
unit_costs::Vector{Float64}
roughness::Float64
positions::Vector{Tuple{Float64,Float64}}
name::String
end
num_nodes(network::Network) = length(network.nodes)
num_pipes(network::Network) = length(network.pipes)
struct Design
cost::Float64
d::Vector{Float64}
q::Vector{Float64}
h::Vector{Float64}
endWater Network Design
Pick commercial pipe diameters for a fixed water-network layout, at the lowest cost that still meets every demand and the nonlinear Hazen-Williams head loss.
Introduction
A water distribution network is a system of hydraulic elements (pipes, pumps, valves, reservoirs) which are connected together to convey given quantities of water, within prescribed pressures, from sources to consumers. Such a system can be represented as a graph in which the nodes correspond to the sources, consumption points and control elements, and the links correspond to the connecting pipes.
The overall planning process of a water distribution network consists of three phases: layout, design, and operation. Although each phase is dependent on the others, they can be formulated and solved as separate problems.
We consider here the design phase. The layout has already been determined during the previous phase, that is, the choice of between which nodes to build a pipe is already fixed. Here, we focus on determining the diameters of the pipes. For simplicity of presentation we consider here simple networks, which do not contain pumps or reservoirs.
Moreover, we deal here with pressurized water networks, where the fluid is transported in pipes with no air contact and thus possibly varying pressure levels. What actually induces a flow between two nodes is explained by a hydraulic head difference in the pipe \(a = (i,j)\) between these nodes, where the flow goes from node \(i\) to node \(j\), which can be modeled using the Hazen-Williams empirical equation:
\[ h_i - h_j = 10.67 \cdot \left(\dfrac{q_a}{K}\right)^{1.852} \cdot L_a \cdot d_a^{-4.87} \]
with
- \(h_i\) and \(h_j\) the hydraulic heads at nodes \(i\) and \(j\), in meters
- \(q_a\) is the discharge, in cubic meters per second
- \(L_a\) is the pipe length, in meters
- \(d_a\) is the pipe diameter, in meters
- \(K\) is the Hazen-Williams coefficient, depending on pipe material
The layout is given, and so is the demand each node draws off or supplies. The ruler under the network is the set of available diameters, and every pipe is fitted with one size off it.
Problem description
Input
- A network represented as a directed graph \(G = (N, A)\), where nodes stand for sources and junctions, and arcs stand for pipes.
- For each node \(i \in N\):
- \(Dem_i\) the demand at node \(i\) (positive if junction, negative if source), in cubic meters per second
- \(E_i\) physical elevation of node \(i\), in meters
- \(P^{\text{min}}_i\) the minimum pressure at node \(i\), in meters
- \(P^{\text{max}}_i\) the maximum pressure at node \(i\), in meters
- For each pipe \(a \in A\):
- \(L_a\) the length of pipe \(a\), in meters
- \(D^{\text{min}}_a\) the minimum diameter of pipe \(a\), in meters
- \(D^{\text{max}}_a\) the maximum diameter of pipe \(a\), in meters
- \(V^{\text{max}}_a\) the flow’s maximum velocity in pipe \(a\)
- A discrete set \(\{D_1, \dots, D_L\}\) of \(L\) commercially-available diameters for the pipes; for each diameter \(D_l\), \(l = 1, \dots, L\), a cost of a linear meter of pipe \(C_l\)
- The Hazen-Williams coefficient \(K\), identical for all pipes in this example.
Problem:
Choose the diameters of the pipes among the set of available diameters and the water flow values such that:
- Each node demand is satisfied
- Water flows satisfy the hydraulic constraints: Hazen-Williams equation, flow conservation and flow bounds
Objective:
Minimize the cost of the pipes (which depends on the selected diameters).
Mixed-integer nonlinear model
Variables
- \(q_a \in \mathbb{R}\), \(a \in A\), the flow in pipe \(a\), in cubic meters per second
- \(d_a\), \(D^{\text{min}}_a \le d_a \le D^{\text{max}}_a\), \(a \in A\), the diameter of pipe \(a\), in meters
- \(c_a \in \mathbb{R}\), \(a \in A\), the cost of a linear meter of pipe \(a\)
- \(h_i \in [P^{\text{min}}_i + E_i; P^{\text{max}}_i + E_i]\), \(i \in N\), the hydraulic head at node \(i\), in meters
- \(x_{a,l} \in \{0,1\}\), \(a \in A\), \(l = 1, \dots, L\): \(x_{a,l} = 1\) iff the diameter of pipe \(a\) is at least diameter \(D_l\)
- \(y_a \in \{0,1\}\), \(a \in A\): \(y_a = 1\) iff the flow in pipe \(a\) runs the way the arc was drawn
Objective: minimize the cost of the pipes
\[ \min \sum_{a \in A} L_a \cdot c_a \]
Constraints
- Diameters are in the discrete set \(\{D_1, \dots, D_L\}\): \[ \forall a \in A, \qquad d_a = D_1 x_{a,1} + \sum_{l=2}^L (D_l - D_{l-1})~x_{a,l} \] \[ \forall a \in A, ~ \forall l = 1, \dots, L-1, \qquad x_{a,l} \ge x_{a,l+1} \]
Note that this incremental modeling allows us to use the branch-and-bound method more effectively, while branching \(d_a \le D_l\) vs \(d_a \ge D_{l+1}\) is achieved through ordinary 0/1 branching on the single binary variable \(x_{a,l}\).
Cost of pipes: \[ \forall a \in A, \qquad c_a = C_1 x_{a,1} + \sum_{l=2}^L (C_l - C_{l-1})~x_{a,l} \]
Flow bounds (dependent on cross-sectional area of pipe): \[ \forall a \in A, \qquad -\frac{\pi}{4} d_a^2 \cdot V^{\text{max}}_a \leq q_a \leq \frac{\pi}{4} d_a^2 \cdot V^{\text{max}}_a \]
Flow conservation: \[ \forall i \in N, \qquad \sum_{a \in \delta^-(i)} q_a - \sum_{a \in \delta^+(i)} q_a = Dem_i \]
Head loss across links (Hazen-Williams equation): \[ \forall a = (i,j) \in A, \qquad h_i - h_j = \text{sign}(q_a) \cdot 10.67 \cdot \left(\dfrac{|q_a|}{K}\right)^{1.852} \cdot L_a \cdot d_a^{-4.87} \]
To implement this constraint, we linearize the sign and the absolute value. To do this, we introduce \(q_a^+ \ge 0\) and \(q_a^- \ge 0\) as two new variables such that \(q_a = q_a^+ - q_a^-\), and we use the binary \(y_a\) to let only one of them be nonzero. Writing \(Q^{\text{max}}_a = \frac{\pi}{4} (D^{\text{max}}_a)^2 V^{\text{max}}_a\) for the most the pipe can carry:
\[ \forall a \in A, \qquad q_a^+ \leq Q^{\text{max}}_a~y_a, \qquad q_a^- \leq Q^{\text{max}}_a~(1 - y_a) \]
Then \(q_a^+ + q_a^-\) is exactly \(|q_a|\), and the Hazen-Williams equation becomes
\[ \forall a = (i,j) \in A, \qquad (h_i - h_j) \cdot d_a^{4.87} = 10.67 \cdot \dfrac{L_a}{K^{1.852}} \cdot q_a \cdot (q_a^+ + q_a^-)^{0.852} \]
Without \(y_a\) both parts can grow together, and the model buys head loss that its flow never produced.
Once the diameters are fixed, flow conservation and the head-loss law determine every flow and every head, so what the model chooses is the diameters, not where the water goes.
The network
A Node is a demand, an elevation and the pressure band it has to stay inside; the source is the one node that supplies rather than draws. An Arc is a fixed pipe run between two nodes, with the diameter range and the velocity limit it has to respect. A Network holds both lists plus the catalogue of commercial diameters, what each costs per metre, the Hazen-Williams roughness, and where to draw each node.
A Design carries its total cost, the diameter chosen for every pipe, the flow through it, and the hydraulic head at every node.
The classic two-loop network of Alperovits and Shamir (1977) is seven nodes, one source at the top and eight pipes. Demands are in cubic metres per second here, since that is what the head loss equation wants, and are shown in cubic metres per hour everywhere a reader sees them.
function two_loop_network()
demand = [-1120, 100, 100, 120, 270, 330, 200] ./ 3600
elevation = [210, 150, 160, 155, 150, 165, 160]
max_pressure = [0, 60, 50, 55, 60, 45, 50]
min_pressure = [0, 30, 30, 30, 30, 30, 30]
kind = [:source; fill(:junction, 6)]
nodes = Node.(demand, elevation, min_pressure, max_pressure, kind)
from_node = [1, 2, 2, 4, 4, 6, 3, 7]
to_node = [2, 3, 4, 5, 6, 7, 5, 5]
min_diameter = [0.3048, 0.1524, 0.254, 0.0762, 0.254, 0.2032, 0.1524, 0.1524]
max_diameter = [0.508, 0.3556, 0.4572, 0.2032, 0.4572, 0.4064, 0.3556, 0.3556]
pipes = Arc.(from_node, to_node, 1000.0, min_diameter, max_diameter, 2.0)
diameters = INCH .* [1, 2, 3, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
costs = [2, 5, 8, 11, 16, 23, 32, 50, 60, 90, 130, 170, 300, 550]
positions = [
(10.0, 10.0),
(5.0, 10.0),
(0.0, 10.0),
(5.0, 5.0),
(0.0, 5.0),
(5.0, 0.0),
(0.0, 0.0),
]
return Network(
nodes, pipes, diameters, costs, 130.0, positions, "Simple water network"
)
end
network = two_loop_network()A pipe is drawn as thick as its diameter, on a scale shared with the catalogue beside the figure, so the choice the model makes is legible as a width rather than as a number to look up. In the solved figure a lighter core runs inside each pipe, as wide as the flow would need at the pipe’s maximum velocity, and the gap between core and wall is the headroom the design bought.
The dashed arcs show the fixed layout and the arrow on each one is the direction it was drawn in, and the flow need not follow it. The arrow at every node carries the demand it draws off, in cubic metres per hour, or supplies at the source. The scale on the right is the catalogue of diameters each pipe has to choose from, drawn at the thickness the figures use.
draw_layout(network)Model implementation
using JuMP
using KNITROminimize_cost(network) builds the model, solves it with Knitro, and returns a Design. The catalogue, the flow balance and the direction binaries are affine, and the same @constraint macro takes the area bound and the Hazen-Williams head loss with their fractional powers written as they read.
The node limit keeps the search short enough to render this page, and max_nodes raises it for the larger network below.
function minimize_cost(network::Network; multistart=true, max_nodes=512)
pipes, nodes = 1:num_pipes(network), 1:num_nodes(network)
diameters, costs = network.diameters, network.unit_costs
roughness = network.roughness
levels = 1:length(diameters)
size_step = [diameters[1]; diff(diameters)]
cost_step = [costs[1]; diff(costs)]
min_diameter = [pipe.min_diameter for pipe in network.pipes]
max_diameter = [pipe.max_diameter for pipe in network.pipes]
min_head = [node.min_pressure + node.elevation for node in network.nodes]
max_head = [node.max_pressure + node.elevation for node in network.nodes]
model = Model(KNITRO.Optimizer)
@variable(model, q[pipes])
@variable(model, c[pipes])
@variable(model, q_plus[pipes] >= 0)
@variable(model, q_minus[pipes] >= 0)
@variable(model, y[pipes], Bin)
@variable(model, x[pipes, levels], Bin)
@variable(model, min_diameter[a] <= d[a in pipes] <= max_diameter[a])
@variable(model, min_head[i] <= h[i in nodes] <= max_head[i])
@objective(model, Min, sum(network.pipes[a].length * c[a] for a in pipes))
for a in pipes
pipe = network.pipes[a]
@constraint(model, d[a] == sum(size_step[k] * x[a, k] for k in levels))
@constraint(model, c[a] == sum(cost_step[k] * x[a, k] for k in levels))
@constraint(model, [k in levels[1:(end - 1)]], x[a, k] >= x[a, k + 1])
section = pi / 4 * pipe.max_velocity * d[a]^2
@constraint(model, q[a] <= section)
@constraint(model, q[a] >= -section)
@constraint(model, q[a] == q_plus[a] - q_minus[a])
max_flow = pi / 4 * pipe.max_diameter^2 * pipe.max_velocity
@constraint(model, q_plus[a] <= max_flow * y[a])
@constraint(model, q_minus[a] <= max_flow * (1 - y[a]))
resistance = 10.67 * pipe.length / roughness^1.852
@constraint(
model,
d[a]^4.87 * (h[pipe.from_node] - h[pipe.to_node]) ==
resistance * q[a] * (q_plus[a] + q_minus[a])^0.852
)
end
@constraint(
model,
[i in nodes],
sum(q[a] for a in pipes if network.pipes[a].to_node == i) -
sum(q[a] for a in pipes if network.pipes[a].from_node == i) ==
network.nodes[i].demand
)
set_attribute(model, "mip_multistart", Int(multistart))
set_attribute(model, "mip_maxnodes", max_nodes)
optimize!(model)
cost = objective_value(model)
return Design(cost, value.(d), value.(q), value.(h))
endOutput visualization
Solving chooses a commercial diameter for each pipe. Each pipe is drawn as thick as the diameter chosen for it, and the light core as thick as the flow it carries would need at its maximum velocity. Beside every pipe sits the pair the design comes down to, the diameter in inches and, under it in blue, the flow through it in cubic metres per hour. The blue arrowhead points the way the water goes, and that is not always the way the arc was drawn.
design = minimize_cost(network);=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 12 variables (7%) and 10 constraints (6%) in 0.00s.
datacheck 0
feastol 1e-06
feastol_abs 1e-06
hessian_no_f 1
mip_maxnodes 512
mip_multistart 1
mip_numthreads 1
ms_numthreads 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
The Knitro mixed integer solver is designed for convex problems.
For non-convex problems it is only a heuristic, and the reported
bounds and optimality claims cannot be verified.
Problem Characteristics | Presolved
-----------------------
Problem type: MINLP
Objective: minimize / linear
Number of variables: 167 | 155
bounds: lower upper range | lower upper range
16 0 134 | 14 0 133
free fixed | free fixed
16 1 | 8 0
cont. binary integer | cont. binary integer
47 120 0 | 36 119 0
Number of constraints: 175 | 165
eq. ineq. range | eq. ineq. range
linear: 31 120 0 | 22 119 0
quadratic: 0 16 0 | 0 16 0
nonlinear: 8 0 0 | 8 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 8 536 | 112 407
quadratic: 0 16 8 | 0 16 8
nonlinear: 0 48 79 | 0 48 72
total: 8 600 79 | 112 468 72
Knitro using Branch and Bound method with 1 thread.
Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.
Coefficient range:
linear objective: [1e+03, 1e+03] | [8e-01, 1e+02]
linear constraints: [3e-02, 2e+02] | [3e-02, 1e+00]
quadratic objective: [0e+00, 0e+00] | [0e+00, 0e+00]
quadratic constraints: [2e+00, 2e+00] | [2e+00, 2e+00]
variable bounds: [8e-02, 2e+02] | [8e-02, 2e+02]
constraint bounds: [3e-02, 4e-01] | [3e-02, 3e-01]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 2.27922e+06 0.657856 39.4898 3.263
1 1.18025e+06 0.330512 9.18097 13.058
2 848732. 0.133435 106.578 13.059
3 652475. 8.42832e-02 10.2552 13.061
4 486740. 4.17998e-02 2.05530 13.062
5 362062. 7.60899e-02 0.889573 13.063
6 366679. 6.64009e-02 38.6113 13.063
7 427385. 0.138880 216.271 13.068
8 568252. 0.332811 145.116 13.068
9 715213. 0.225105 16.0000 13.068
10 846073. 0.199346 50.0604 13.069
11 1.05868e+06 0.229485 21.2122 13.069
12 1.26764e+06 0.307269 51.6973 13.070
13 1.50639e+06 0.351283 51.9160 13.070
14 1.67217e+06 0.361092 51.9826 13.070
15 1.78405e+06 0.361215 100.000 13.071
16 1.57123e+06 0.134125 100.000 13.071
17 1.34223e+06 0.107259 100.000 13.072
18 1.22155e+06 8.18914e-02 100.000 13.072
19 1.33585e+06 2.15891e-02 100.000 13.073
20 1.39031e+06 1.65018e-02 100.000 13.073
30 478977. 3.57268e-04 65.0150 13.081
Root node cutting planes
------------------------
Iter Cuts Best solution Best bound Gap Time
value value (secs)
---- ---- ------------- ---------- --- ------
0 0 -inf 16.130
1 1 -inf 16.177
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 16.200
1 2 464000. FP -inf 16.394
1 2 463000. MPEC -inf 17.195
Knitro deduced that the problem is non-convex.
80 71 463000. MPEC -inf 20.014
100 86 463000. -inf 20.430
114 93 444000. LEAF -inf 20.708
200 109 444000. -inf 21.858
300 109 444000. -inf 22.634
400 108 444000. -inf 23.488
500 108 444000. -inf 24.428
512 108 444000. -inf 24.553
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 4.44000000221460476e+05
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 2
# of restarts = 0
# of nodes processed = 512 (17.179s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 22409 (0.224s)
# of gradient evaluations = 15676 (0.168s)
# of hessian evaluations = 13404 (10.026s)
# of hessian-vector evaluations = 0
# of subproblems processed = 544 (17.997s)
Total program time (secs) = 24.55286 (20.751 CPU time)
Time spent in evaluations (secs) = 10.41726
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 1 / 1
Flow-cover cuts = 0 / 0
Probing cuts = 6 / 18
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 2 / 2 / 0.197s
Rounding heuristic = 3 / 0 / 1.424s
MPEC heuristic = 4 / 2 / 0.591s
Local search heuristic = 9 / 0 / 0.031s
===========================================================================
Show the full outputHide the full output
draw_design(network, design)report_design(network, design)Simple water network, total pipe cost $444,000
pipe diameter flow cost share
(m³/h) ($)
--------------------------------------------------
0 to 1 18" 1,120.0 130,000 29.3%
1 to 2 14" 446.2 60,000 13.5%
1 to 3 16" 573.8 90,000 20.3%
3 to 4 3" 9.6 8,000 1.8%
3 to 5 14" 444.2 60,000 13.5%
5 to 6 8" 114.2 23,000 5.2%
2 to 4 12" 346.2 50,000 11.3%
6 to 4 8" -85.8 23,000 5.2%
node demand elevation head pressure margin
(m³/h) (m) (m) (m) (m)
------------------------------------------------------------
0 -1,120 210 210.00 0.00 0.00
1 100 150 203.25 53.25 23.25
2 100 160 199.08 39.08 9.08
3 120 155 199.78 44.78 14.78
4 270 150 193.55 43.55 13.55
5 330 165 195.64 30.64 0.64
6 200 160 190.54 30.54 0.54
The tightest junction is node 6, 0.54 m above its minimum pressureKnitro puts the money where the flow is. The widest pipes leave the source, and the ones carrying little take the smallest diameter their bounds allow. The margin column is what each node has in hand over its minimum pressure, and the nodes sitting at zero are the ones holding the design where it is. The search stops at the node limit set above, so this is the best design it found rather than one proven optimal.
A larger network
Seven nodes fit on a page. The Hanoi network of Fujiwara and Khang (1990) is the usual next step and a standard benchmark: 32 nodes, 34 pipes, one source held at 100 m, and a catalogue of six diameters from 12 to 40 inches.
Show the Hanoi network
function hanoi_network()
demand =
[
-19940,
890,
850,
130,
725,
1005,
1350,
550,
525,
525,
500,
560,
940,
615,
280,
310,
865,
1345,
60,
1275,
930,
485,
1045,
820,
170,
900,
370,
290,
360,
360,
105,
805,
] ./ 3600
min_pressure = [100; fill(30, 31)]
kind = [:source; fill(:junction, 31)]
nodes = [Node(demand[i], 0.0, min_pressure[i], 100.0, kind[i]) for i in 1:32]
from_node = [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
10,
14,
15,
17,
18,
19,
3,
3,
20,
21,
20,
23,
24,
26,
27,
16,
23,
28,
29,
30,
32,
25,
]
to_node = [
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
16,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
25,
26,
27,
28,
29,
30,
31,
31,
32,
]
length = [
100,
1350,
900,
1150,
1450,
450,
850,
850,
800,
950,
1200,
3500,
800,
500,
550,
2730,
1750,
800,
400,
2200,
1500,
500,
2650,
1230,
1300,
850,
300,
750,
1500,
2000,
1600,
150,
860,
950,
]
max_velocity = [
7,
7,
3,
3,
2.5,
2.5,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
3.5,
3.5,
3,
2,
2,
2,
3,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
]
pipes = Arc.(from_node, to_node, length, 0.3048, 1.016, max_velocity)
diameters = [0.3048, 0.4064, 0.508, 0.6096, 0.762, 1.016]
costs = [45.73, 70.40, 98.39, 129.33, 180.75, 278.28]
positions = [
(12.67525, 0.0),
(12.67525, 3.91355),
(12.5584, 7.41825),
(16.76405, 7.41825),
(19.8014, 7.41825),
(22.72195, 7.41825),
(22.72195, 10.57245),
(22.72195, 14.0771),
(22.72195, 17.52335),
(20.3855, 17.52335),
(20.3855, 20.3271),
(20.3855, 22.3715),
(17.1729, 22.3715),
(17.93225, 17.52335),
(14.83645, 17.52335),
(12.5, 17.52335),
(12.5, 15.24535),
(12.5584, 13.25935),
(12.5584, 10.2804),
(9.17055, 7.41825),
(9.17055, 4.0888),
(9.17055, 0.5257),
(6.3084, 7.41825),
(6.3084, 13.4346),
(6.1916, 17.52335),
(8.46965, 17.52335),
(10.51405, 17.52335),
(3.44625, 7.41825),
(0.0, 7.59345),
(0.0, 12.09115),
(0.0, 17.52335),
(3.44625, 17.52335),
]
return Network(
nodes, pipes, diameters, costs, 130.0, positions, "Hanoi water network"
)
endhanoi_network (generic function with 1 method)hanoi = hanoi_network()draw_layout(hanoi)Thirty-four pipes over six catalogue levels make for a far larger tree, and a deeper search is needed before it runs out of ideas. Too many nodes to carry their demands as well, so the figure below keeps the diameter and the flow on every pipe and leaves the demands to the layout above it.
hanoi_design = minimize_cost(hanoi; max_nodes=2048);=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 56 variables (13%) and 48 constraints (10%) in 0.00s.
datacheck 0
feastol 1e-06
feastol_abs 1e-06
hessian_no_f 1
mip_maxnodes 2048
mip_multistart 1
mip_numthreads 1
ms_numthreads 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
The Knitro mixed integer solver is designed for convex problems.
For non-convex problems it is only a heuristic, and the reported
bounds and optimality claims cannot be verified.
Problem Characteristics | Presolved
-----------------------
Problem type: MINLP
Objective: minimize / linear
Number of variables: 440 | 384
bounds: lower upper range | lower upper range
68 0 303 | 54 0 296
free fixed | free fixed
68 1 | 34 0
cont. binary integer | cont. binary integer
202 238 0 | 153 231 0
Number of constraints: 474 | 426
eq. ineq. range | eq. ineq. range
linear: 134 238 0 | 93 231 0
quadratic: 0 68 0 | 0 68 0
nonlinear: 34 0 0 | 34 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 34 1190 | 204 889
quadratic: 0 68 34 | 0 68 34
nonlinear: 0 204 338 | 0 204 301
total: 34 1462 338 | 204 1146 301
Knitro using Branch and Bound method with 1 thread.
Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.
Coefficient range:
linear objective: [1e+02, 4e+03] | [7e-01, 1e+02]
linear constraints: [1e-01, 1e+02] | [1e-01, 3e+00]
quadratic objective: [0e+00, 0e+00] | [0e+00, 0e+00]
quadratic constraints: [2e+00, 5e+00] | [2e+00, 5e+00]
variable bounds: [3e-01, 1e+02] | [3e-01, 1e+02]
constraint bounds: [2e-02, 6e+00] | [2e-02, 6e+00]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 5.48733e+06 38.3763 33.1366 0.006
1 4.50167e+06 37.8946 8.94413 0.009
2 4.50589e+06 35.2221 24.7210 0.010
3 4.55368e+06 32.4800 16.1406 0.011
4 4.62326e+06 27.6314 15.6910 0.012
5 4.68113e+06 19.6543 15.0828 0.013
6 4.79811e+06 15.3609 10.8317 0.014
7 5.10877e+06 11.0681 13.7016 0.016
8 5.00950e+06 1.88386 1.09439 0.018
9 5.60633e+06 0.547690 12.8639 0.019
10 5.68893e+06 0.358626 10.5301 0.020
11 5.74441e+06 0.293297 26.7861 0.021
12 5.84723e+06 0.205332 26.1380 0.023
13 5.96898e+06 0.263592 40.3661 0.024
14 6.00716e+06 0.194033 34.0655 0.025
15 6.02903e+06 0.163899 58.0006 0.026
16 6.04058e+06 0.150842 45.8800 0.027
17 6.04914e+06 0.100766 0.472275 0.029
18 6.03883e+06 7.90910e-02 3.51037 0.032
19 6.04667e+06 6.94146e-02 3.16603 0.033
20 6.04933e+06 6.51039e-02 3.24381 0.034
30 6.03502e+06 4.67797e-05 1.85419 0.049
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 0.061
1 2 6.37937e+06 FP -inf 0.681
Knitro deduced that the problem is non-convex.
100 84 6.37937e+06 -inf 9.498
200 174 6.37937e+06 -inf 18.853
300 267 6.37937e+06 -inf 27.480
400 354 6.37937e+06 -inf 35.091
500 437 6.37937e+06 -inf 43.371
600 523 6.37937e+06 -inf 52.604
700 608 6.37937e+06 -inf 62.851
750 651 6.27735e+06 MPEC -inf 66.719
800 690 6.27735e+06 -inf 70.375
900 774 6.27735e+06 -inf 81.364
1000 860 6.27735e+06 -inf 92.078
1100 948 6.27735e+06 -inf 102.874
1200 1030 6.27735e+06 -inf 113.214
1300 1121 6.27735e+06 -inf 121.977
1400 1207 6.27735e+06 -inf 131.594
1500 1289 6.27735e+06 -inf 142.319
1600 1370 6.27735e+06 -inf 152.431
1700 1453 6.27735e+06 -inf 162.830
1800 1528 6.27735e+06 -inf 174.088
1900 1608 6.27735e+06 -inf 186.296
2000 1695 6.27735e+06 -inf 198.169
2048 1739 6.27735e+06 -inf 203.227
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 6.27734599959880300e+06
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 2048 (183.628s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 208085 (7.315s)
# of gradient evaluations = 132687 (4.984s)
# of hessian evaluations = 120239 (8.340s)
# of hessian-vector evaluations = 0
# of subproblems processed = 2152 (200.725s)
Total program time (secs) = 203.22716 (203.081 CPU time)
Time spent in evaluations (secs) = 20.64010
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 63 / 63
Flow-cover cuts = 0 / 0
Probing cuts = 0 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 5 / 1 / 1.171s
Rounding heuristic = 1 / 0 / 0.222s
MPEC heuristic = 32 / 1 / 15.961s
Local search heuristic = 7 / 0 / 0.055s
===========================================================================
Show the full outputHide the full output
draw_design(hanoi, hanoi_design)report_design(hanoi, hanoi_design)Hanoi water network, total pipe cost $6,277,346
pipe diameter flow cost share
(m³/h) ($)
--------------------------------------------------
0 to 1 40" 19,940.0 27,828 0.4%
1 to 2 40" 19,049.10 375,678 6.0%
2 to 3 40" 7,960.10 250,452 4.0%
3 to 4 40" 7,830.10 320,022 5.1%
4 to 5 40" 7,105.10 403,506 6.4%
5 to 6 40" 6,100.10 125,226 2.0%
6 to 7 40" 4,750.10 236,538 3.8%
7 to 8 40" 4,200.10 236,538 3.8%
8 to 9 40" 3,675.10 222,624 3.5%
9 to 10 30" 2,000.0 171,712 2.7%
10 to 11 24" 1,500.0 155,196 2.5%
11 to 12 24" 940.0 452,655 7.2%
9 to 13 20" 1,150.10 78,712 1.3%
13 to 14 16" 535.10 35,200 0.6%
14 to 15 16" 255.10 38,720 0.6%
16 to 15 16" 491.3 192,192 3.1%
17 to 16 20" 1,356.3 172,182 2.7%
18 to 17 24" 2,701.3 103,464 1.6%
2 to 18 24" 2,761.3 51,732 0.8%
2 to 19 40" 7,477.7 612,216 9.8%
19 to 20 20" 1,415.0 147,585 2.4%
20 to 21 16" 485.0 35,200 0.6%
19 to 22 40" 4,787.7 737,442 11.7%
22 to 23 30" 3,167.8 222,322 3.5%
23 to 24 30" 2,347.8 234,975 3.7%
25 to 24 16" -832.7 59,840 1.0%
26 to 25 16" 67.3 21,120 0.3%
15 to 26 16" 437.3 52,800 0.8%
22 to 27 16" 574.9 105,600 1.7%
27 to 28 12" 284.9 91,460 1.5%
28 to 29 16" -75.1 112,640 1.8%
29 to 30 16" -435.1 10,560 0.2%
31 to 30 16" 540.1 60,544 1.0%
24 to 31 24" 1,345.1 122,864 2.0%
node demand elevation head pressure margin
(m³/h) (m) (m) (m) (m)
------------------------------------------------------------
0 -19,940 0 100.00 100.00 0.00
1 890 0 97.14 97.14 67.14
2 850 0 61.66 61.66 31.66
3 130 0 56.96 56.96 26.96
4 725 0 51.13 51.13 21.13
5 1,005 0 45.00 45.00 15.00
6 1,350 0 43.56 43.56 13.56
7 550 0 41.85 41.85 11.85
8 525 0 40.50 40.50 10.50
9 525 0 39.50 39.50 9.50
10 500 0 37.94 37.94 7.94
11 560 0 34.51 34.51 4.51
12 940 0 30.30 30.30 0.30
13 615 0 36.10 36.10 6.10
14 280 0 34.57 34.57 4.57
15 310 0 34.14 34.14 4.14
16 865 0 41.25 41.25 11.25
17 1,345 0 51.33 51.33 21.33
18 60 0 58.12 58.12 28.12
19 1,275 0 51.43 51.43 21.43
20 930 0 42.08 42.08 12.08
21 485 0 40.81 40.81 10.81
22 1,045 0 46.03 46.03 16.03
23 820 0 41.30 41.30 11.30
24 170 0 38.43 38.43 8.43
25 900 0 32.55 32.55 2.55
26 370 0 32.57 32.57 2.57
27 290 0 40.81 40.81 10.81
28 360 0 33.10 33.10 3.10
29 360 0 33.23 33.23 3.23
30 105 0 33.54 33.54 3.54
31 805 0 36.21 36.21 6.21
The tightest junction is node 12, 0.30 m above its minimum pressureThe cheapest design published for this network costs about $6.08M, so the node limit leaves a few percent on the table. This is where the nonconvexity shows: Knitro’s bound stays at minus infinity throughout, so there is nothing to measure that gap against from the inside.
References
- “An MINLP Solution Method for a Water Network Problem”, C. Bragalli, C. D’Ambrosio, J. Lee, A. Lodi, P. Toth (2006) DOI
- “Mathematical programming techniques in water network optimization”, C. D’Ambrosio, A. Lodi, S. Wiese, C. Bragalli (2015) DOI
- “Design of Optimal Water Distribution System”, E. Alperovits and U. Shamir (1977) DOI
- “A two-phase decomposition method for optimal design of looped water distribution networks”, O. Fujiwara and D. B. Khang (1990) DOI