import math
from dataclasses import dataclass
from plot import draw_design, draw_layout
from report import report_design
INCH = 0.0254
@dataclass
class Node:
demand: float
elevation: float
min_pressure: float
max_pressure: float
kind: str
@dataclass
class Arc:
from_node: int
to_node: int
length: float
min_diameter: float
max_diameter: float
max_velocity: float
@dataclass
class Network:
nodes: list
pipes: list
diameters: list
unit_costs: list
roughness: float
positions: list
name: str
@property
def num_nodes(self):
return len(self.nodes)
@property
def num_pipes(self):
return len(self.pipes)
@dataclass
class Design:
cost: float
d: list
q: list
h: listWater 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.
def two_loop_network():
demand = [d / 3600 for d in [-1120, 100, 100, 120, 270, 330, 200]]
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"] + ["junction"] * 6
nodes = list(map(Node, demand, elevation, min_pressure, max_pressure, kind))
from_node = [0, 1, 1, 3, 3, 5, 2, 6]
to_node = [1, 2, 3, 4, 5, 6, 4, 4]
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[a], to_node[a], 1000, min_diameter[a], max_diameter[a], 2)
for a in range(8)
]
diameters = [INCH * n for n in [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, 10), (5, 10), (0, 10), (5, 5), (0, 5), (5, 0), (0, 0)]
return Network(
nodes, pipes, diameters, costs, 130, positions, "Simple water network"
)
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
from itertools import pairwise
import 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 add_constraint 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.
def minimize_cost(network, *, multistart=True, max_nodes=512):
pipes = range(network.num_pipes)
nodes = range(network.num_nodes)
diameters, costs = network.diameters, network.unit_costs
roughness = network.roughness
levels = range(len(diameters))
size_step = [diameters[0], *(b - a for a, b in pairwise(diameters))]
cost_step = [costs[0], *(b - a for a, b in pairwise(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]
binary = knitro.KN_VARTYPE_BINARY
prob = knitro.Problem()
q = [prob.add_variable() for _ in pipes]
q_plus = [prob.add_variable(lb=0) for _ in pipes]
q_minus = [prob.add_variable(lb=0) for _ in pipes]
d = [prob.add_variable(lb=min_diameter[a], ub=max_diameter[a]) for a in pipes]
h = [prob.add_variable(lb=min_head[i], ub=max_head[i]) for i in nodes]
x = {(a, k): prob.add_variable(vtype=binary) for a in pipes for k in levels}
y = [prob.add_variable(vtype=binary) for _ in pipes]
c = [prob.add_variable() for _ in pipes]
prob.add_objective(prob.nsum(c[a] * network.pipes[a].length for a in pipes))
for a in pipes:
pipe = network.pipes[a]
prob.add_constraint(d[a] == prob.nsum(size_step[k] * x[a, k] for k in levels))
prob.add_constraint(c[a] == prob.nsum(cost_step[k] * x[a, k] for k in levels))
for k in levels[:-1]:
prob.add_constraint(x[a, k] >= x[a, k + 1])
prob.add_constraint(-math.pi / 4 * d[a] ** 2 * pipe.max_velocity <= q[a])
prob.add_constraint(q[a] <= math.pi / 4 * d[a] ** 2 * pipe.max_velocity)
prob.add_constraint(q[a] == q_plus[a] - q_minus[a])
max_flow = math.pi / 4 * pipe.max_diameter**2 * pipe.max_velocity
prob.add_constraint(q_plus[a] <= max_flow * y[a])
prob.add_constraint(q_minus[a] <= max_flow * (1 - y[a]))
i, j = pipe.from_node, pipe.to_node
gradient = 10.67 * q[a] * (q_plus[a] + q_minus[a]) ** 0.852 / roughness**1.852
prob.add_constraint(d[a] ** 4.87 * (h[i] - h[j]) == gradient * pipe.length)
for i in nodes:
inflow = prob.nsum(q[a] for a in pipes if network.pipes[a].to_node == i)
outflow = prob.nsum(q[a] for a in pipes if network.pipes[a].from_node == i)
prob.add_constraint(inflow - outflow == network.nodes[i].demand)
prob.set_param(knitro.KN_PARAM_MIP_MULTISTART, int(multistart))
prob.set_param(knitro.KN_PARAM_MIP_MAXNODES, max_nodes)
prob.solve()
cost = prob.get_attr(knitro.KN_ATTR_OBJ_VALUE)
return Design(
cost,
[d[a].value for a in pipes],
[q[a].value for a in pipes],
[h[i].value for i in nodes],
)Output 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.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
mip_maxnodes 512
mip_multistart 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 expressions: 768 | 107
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 64 | 0 48 59
total: 8 600 64 | 112 468 59
Knitro using Branch and Bound method with 8 threads.
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, 3e-01] | [3e-02, 3e-01]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 2.03491e+06 0.674180 33.8057 0.111
1 1.38285e+06 0.291103 14.9026 0.112
2 1.19891e+06 0.240401 29.6168 0.114
3 879763. 0.121619 137.138 0.115
4 759249. 5.74066e-02 40.6096 0.115
5 692192. 3.20136e-02 200.102 0.116
6 615203. 1.97273e-02 1.13754 0.117
7 576681. 1.04844e-02 0.661926 0.118
8 483792. 7.43952e-04 11.6312 0.119
9 447390. 2.56561e-04 43.9700 0.119
10 434238. 5.50950e-05 4.50925 0.120
11 432904. 2.91056e-05 8.91488 0.120
12 430553. 2.87074e-05 0.988761 0.121
13 430037. 2.78354e-05 1.37292 0.121
14 429986. 2.58505e-05 2.59755 0.122
15 429907. 2.13305e-05 3.67867 0.122
16 429784. 1.20048e-05 3.25976 0.123
17 429677. 2.61106e-06 0.681554 0.123
18 429647. 4.05387e-07 1.87600e-02 0.123
19 429642. 1.22467e-07 5.75130e-04 0.124
20 429641. 4.00418e-08 5.35282e-05 0.124
Root node cutting planes
------------------------
Iter Cuts Best solution Best bound Gap Time
value value (secs)
---- ---- ------------- ---------- --- ------
0 0 -inf 0.821
1 1 -inf 0.872
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 0.889
2 3 465000. MPEC -inf 0.942
Knitro deduced that the problem is non-convex.
123 98 447000. FCRD -inf 1.585
139 111 444000. LEAF -inf 1.618
514 124 444000. -inf 2.256
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 4.44000000036096200e+05
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 2
# of restarts = 0
# of nodes processed = 514 (7.690s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 21680 (0.168s)
# of gradient evaluations = 16844 (0.159s)
# of hessian evaluations = 14030 (0.291s)
# of hessian-vector evaluations = 0
# of subproblems processed = 615 (8.762s)
Total program time (secs) = 2.25905 (8.860 CPU time)
Time spent in evaluations (secs) = 0.61852
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 9 / 9
Flow-cover cuts = 0 / 0
Probing cuts = 6 / 12
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 2 / 2 / 0.329s
Rounding heuristic = 4 / 1 / 0.021s
MPEC heuristic = 3 / 1 / 0.740s
Local search heuristic = 8 / 0 / 0.034s
===========================================================================
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
def hanoi_network():
demand = [
d / 3600
for d in [
-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,
]
]
min_pressure = [100] + [30] * 31
kind = ["source"] + ["junction"] * 31
nodes = [Node(demand[i], 0, min_pressure[i], 100, kind[i]) for i in range(32)]
from_node = [
0,
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
9,
13,
14,
16,
17,
18,
2,
2,
19,
20,
19,
22,
23,
25,
26,
15,
22,
27,
28,
29,
31,
24,
]
to_node = [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
24,
25,
26,
27,
28,
29,
30,
30,
31,
]
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[a], to_node[a], length[a], 0.3048, 1.016, max_velocity[a])
for a in range(34)
]
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, positions, "Hanoi water network"
)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.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
mip_maxnodes 2048
mip_multistart 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 expressions: 2097 | 444
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 272 | 0 204 243
total: 34 1462 272 | 204 1146 243
Knitro using Branch and Bound method with 8 threads.
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.36381e+06 38.3745 35.8178 0.005
1 4.31976e+06 37.8556 7.35799 0.008
2 4.36618e+06 33.9685 84.9390 0.009
3 4.39301e+06 27.5006 40.9296 0.010
4 4.49105e+06 24.4181 34.7411 0.011
5 4.67517e+06 16.8567 13.9834 0.012
6 4.91990e+06 11.8304 14.7923 0.013
7 5.06154e+06 9.74299 16.3380 0.014
8 5.49396e+06 5.41242 13.6300 0.016
9 5.48754e+06 5.25743 12.2948 0.017
10 5.67397e+06 3.39670 11.9458 0.018
11 5.78344e+06 2.30572 16.2924 0.019
12 5.86814e+06 1.81686 12.5267 0.020
13 5.93546e+06 1.18391 8.71263 0.022
14 5.97788e+06 1.01103 17.1525 0.023
15 6.00562e+06 0.781613 11.6426 0.024
16 6.01228e+06 0.744578 16.6390 0.025
17 6.02935e+06 0.666555 15.4640 0.026
18 6.07832e+06 0.403155 33.0912 0.027
19 6.08704e+06 0.372200 13.9177 0.028
20 6.09849e+06 7.54804e-02 0.304853 0.033
30 6.03499e+06 6.82347e-08 0.592700 0.053
40 6.03499e+06 1.05950e-11 0.223740 0.079
WARNING: Evaluations seem expensive and concurrent evaluations are disabled,
please consider reducing the number of threads for the branch-and-bound
using option 'mip_numthreads'.
Recommended value: 'mip_numthreads=7'
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 0.091
Knitro deduced that the problem is non-convex.
15 14 6.37937e+06 FP -inf 0.630
293 247 6.30282e+06 MPEC -inf 3.871
487 410 6.30282e+06 -inf 5.900
961 810 6.30282e+06 -inf 13.203
1448 1226 6.30282e+06 -inf 19.967
1862 1580 6.27134e+06 MPEC -inf 25.407
1935 1633 6.27134e+06 -inf 26.357
2049 1728 6.27134e+06 -inf 27.922
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 6.27134449998810515e+06
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 2049 (177.279s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 222203 (8.207s)
# of gradient evaluations = 138311 (6.208s)
# of hessian evaluations = 124736 (11.470s)
# of hessian-vector evaluations = 0
# of subproblems processed = 2152 (194.161s)
Total program time (secs) = 27.93290 (197.955 CPU time)
Time spent in evaluations (secs) = 25.88571
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 46 / 46
Flow-cover cuts = 0 / 0
Probing cuts = 0 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 5 / 1 / 1.321s
Rounding heuristic = 1 / 0 / 0.336s
MPEC heuristic = 37 / 2 / 15.592s
Local search heuristic = 8 / 0 / 0.088s
===========================================================================
Show the full outputHide the full output
draw_design(hanoi, hanoi_design)report_design(hanoi, hanoi_design)Hanoi water network, total pipe cost $6,271,344
pipe diameter flow cost share
(m³/h) ($)
--------------------------------------------------
0 to 1 40" 19,940.0 27,828 0.4%
1 to 2 40" 19,050.0 375,678 6.0%
2 to 3 40" 7,737.8 250,452 4.0%
3 to 4 40" 7,607.8 320,022 5.1%
4 to 5 40" 6,882.8 403,506 6.4%
5 to 6 40" 5,877.8 125,226 2.0%
6 to 7 40" 4,527.8 236,538 3.8%
7 to 8 40" 3,977.8 236,538 3.8%
8 to 9 40" 3,452.8 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 16" 927.8 56,320 0.9%
13 to 14 16" 312.8 35,200 0.6%
14 to 15 16" 32.8 38,720 0.6%
16 to 15 16" 503.8 192,192 3.1%
17 to 16 20" 1,368.8 172,182 2.7%
18 to 17 24" 2,713.8 103,464 1.6%
2 to 18 24" 2,773.8 51,732 0.8%
2 to 19 40" 7,688.4 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,998.4 737,442 11.8%
22 to 23 30" 3,366.3 222,322 3.5%
23 to 24 30" 2,546.3 234,975 3.7%
25 to 24 20" -1,043.4 83,632 1.3%
26 to 25 12" -143.4 13,719 0.2%
15 to 26 16" 226.6 52,800 0.8%
22 to 27 16" 587.1 105,600 1.7%
27 to 28 12" 297.1 91,460 1.5%
28 to 29 16" -62.9 112,640 1.8%
29 to 30 16" -422.9 10,560 0.2%
31 to 30 16" 527.9 60,544 1.0%
24 to 31 24" 1,332.9 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 57.20 57.20 27.20
4 725 0 51.68 51.68 21.68
5 1,005 0 45.89 45.89 15.89
6 1,350 0 44.55 44.55 14.55
7 550 0 42.99 42.99 12.99
8 525 0 41.76 41.76 11.76
9 525 0 40.88 40.88 10.88
10 500 0 39.32 39.32 9.32
11 560 0 35.89 35.89 5.89
12 940 0 31.68 31.68 1.68
13 615 0 34.11 34.11 4.11
14 280 0 33.55 33.55 3.55
15 310 0 33.54 33.54 3.54
16 865 0 40.99 40.99 10.99
17 1,345 0 51.24 51.24 21.24
18 60 0 58.09 58.09 28.09
19 1,275 0 50.89 50.89 20.89
20 930 0 41.54 41.54 11.54
21 485 0 40.27 40.27 10.27
22 1,045 0 45.04 45.04 15.04
23 820 0 39.75 39.75 9.75
24 170 0 36.41 36.41 6.41
25 900 0 33.40 33.40 3.40
26 370 0 33.07 33.07 3.07
27 290 0 39.61 39.61 9.61
28 360 0 31.28 31.28 1.28
29 360 0 31.37 31.37 1.37
30 105 0 31.67 31.67 1.67
31 805 0 34.23 34.23 4.23
The tightest junction is node 28, 1.28 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