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 pyomo.environ as pyo
SOLVER_NAME = "knitroampl"minimize_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 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]
model = pyo.ConcreteModel()
model.pipes = pyo.RangeSet(0, network.num_pipes - 1)
model.nodes = pyo.RangeSet(0, network.num_nodes - 1)
model.levels = pyo.RangeSet(0, len(diameters) - 1)
model.q = pyo.Var(model.pipes)
model.q_plus = pyo.Var(model.pipes, within=pyo.NonNegativeReals)
model.q_minus = pyo.Var(model.pipes, within=pyo.NonNegativeReals)
model.d = pyo.Var(
model.pipes, bounds=lambda m, a: (min_diameter[a], max_diameter[a])
)
model.h = pyo.Var(model.nodes, bounds=lambda m, i: (min_head[i], max_head[i]))
model.x = pyo.Var(model.pipes, model.levels, within=pyo.Binary)
model.y = pyo.Var(model.pipes, within=pyo.Binary)
model.c = pyo.Var(model.pipes)
total_cost = pyo.quicksum(model.c[a] * network.pipes[a].length for a in model.pipes)
model.obj = pyo.Objective(expr=total_cost, sense=pyo.minimize)
model.con = pyo.ConstraintList()
for a in pipes:
pipe = network.pipes[a]
model.con.add(model.d[a] == sum(size_step[k] * model.x[a, k] for k in levels))
model.con.add(model.c[a] == sum(cost_step[k] * model.x[a, k] for k in levels))
for k in levels[:-1]:
model.con.add(model.x[a, k] >= model.x[a, k + 1])
model.con.add(-math.pi / 4 * model.d[a] ** 2 * pipe.max_velocity <= model.q[a])
model.con.add(model.q[a] <= math.pi / 4 * model.d[a] ** 2 * pipe.max_velocity)
model.con.add(model.q[a] == model.q_plus[a] - model.q_minus[a])
max_flow = math.pi / 4 * pipe.max_diameter**2 * pipe.max_velocity
model.con.add(model.q_plus[a] <= max_flow * model.y[a])
model.con.add(model.q_minus[a] <= max_flow * (1 - model.y[a]))
i, j = pipe.from_node, pipe.to_node
magnitude = model.q_plus[a] + model.q_minus[a]
gradient = 10.67 * model.q[a] * magnitude**0.852 / roughness**1.852
head_drop = model.h[i] - model.h[j]
model.con.add(model.d[a] ** 4.87 * head_drop == gradient * pipe.length)
for i in nodes:
inflow = sum(model.q[a] for a in pipes if network.pipes[a].to_node == i)
outflow = sum(model.q[a] for a in pipes if network.pipes[a].from_node == i)
model.con.add(inflow - outflow == network.nodes[i].demand)
solver = pyo.SolverFactory(SOLVER_NAME)
solver.options["mip_multistart"] = int(multistart)
solver.options["mip_maxnodes"] = max_nodes
solver.solve(model, tee=True)
cost = model.obj()
return Design(
cost,
[model.d[a]() for a in pipes],
[model.q[a]() for a in pipes],
[model.h[i]() 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)Artelys Knitro 16.0.0: mip_multistart=1
mip_maxnodes=512
=======================================
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
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 1
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 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, 4e-01] | [3e-02, 3e-01]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 2.42348e+06 0.674180 42.3928 0.077
1 1.22629e+06 0.281618 11.9243 0.090
2 827503. 0.135787 59.4059 0.090
3 593071. 6.82983e-02 0.594011 0.092
4 549645. 3.83092e-02 13.1237 0.093
5 590136. 2.03345e-02 17.1714 0.093
6 492211. 9.84040e-04 9.04826 0.094
7 451733. 2.38631e-04 9.58173 0.094
8 427860. 8.25818e-05 3.37539 0.095
9 420826. 1.86729e-05 5.07084 0.095
10 420233. 6.45367e-06 0.474508 0.095
11 420226. 7.45846e-06 0.346456 0.096
12 420226. 1.01450e-06 3.91670e-03 0.096
13 420226. 1.54462e-07 5.10997e-03 0.097
14 420226. 1.05348e-09 1.90367e-06 0.097
15 420226. 1.05348e-09 1.90367e-06 0.098
16 420226. 1.05348e-09 1.90367e-06 1.119
Root node cutting planes
------------------------
Iter Cuts Best solution Best bound Gap Time
value value (secs)
---- ---- ------------- ---------- --- ------
0 0 -inf 1.712
1 1 -inf 1.779
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 1.806
2 3 473000. DURD -inf 1.823
Knitro deduced that the problem is non-convex.
10 11 464000. MPEC -inf 1.926
15 16 455000. MPEC -inf 1.971
125 95 447000. FCRD -inf 2.322
148 104 444000. LEAF -inf 2.384
429 110 444000. -inf 3.056
519 110 444000. -inf 3.169
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 4.44000000231566373e+05
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 3
# of restarts = 0
# of nodes processed = 519 (8.564s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 24210 (0.166s)
# of gradient evaluations = 18099 (0.070s)
# of hessian evaluations = 15238 (0.186s)
# of hessian-vector evaluations = 0
# of subproblems processed = 639 (10.015s)
Total program time (secs) = 3.17261 (9.369 CPU time)
Time spent in evaluations (secs) = 0.42216
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 2 / 2
Flow-cover cuts = 0 / 0
Probing cuts = 6 / 3
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 2 / 2 / 0.418s
Rounding heuristic = 4 / 1 / 0.068s
MPEC heuristic = 5 / 2 / 0.875s
Local search heuristic = 10 / 0 / 0.041s
===========================================================================
WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
- termination condition: maxIterations
- message from solver: Knitro 16.0.0\x3a MIP\x3a Node limit reached.
Integer feasible point found.; objective 444000.0002315664; optimality
gap Infinity; 519 nodes; 639 subproblem solvesShow 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)Artelys Knitro 16.0.0: mip_multistart=1
mip_maxnodes=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
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 1
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 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.46850e+06 38.3533 12.0343 0.079
1 4.76699e+06 74.0594 10.9364 0.083
2 4.62278e+06 59.5365 10.9364 0.084
3 4.58408e+06 58.6353 10.9364 0.085
4 3.80691e+06 31.7799 1.89770 0.088
5 3.50661e+06 25.0087 2.45898 0.091
6 3.31767e+06 22.2640 2.90373 0.094
7 3.22788e+06 3.12823 1.96132 0.097
8 3.43911e+06 2.81615 20.9679 0.098
9 3.44730e+06 2.66484 76.7046 0.100
10 3.43226e+06 2.53350 0.654835 0.103
11 3.58460e+06 2.32841 1.66944 0.104
12 3.61442e+06 2.24506 0.772870 0.107
13 3.66139e+06 2.19626 13.4477 0.108
14 3.74666e+06 2.05567 12.5245 0.109
15 3.77851e+06 2.00297 11.3482 0.110
16 3.81253e+06 1.94451 12.6505 0.111
17 3.82386e+06 1.92413 10.2036 0.112
18 3.85778e+06 1.86801 9.90568 0.113
19 3.89481e+06 1.79206 5.39902 0.117
20 3.92479e+06 1.74255 6.03900 0.118
30 4.73907e+06 0.781960 13.2103 0.135
40 5.40067e+06 0.469952 10.9902 0.145
50 5.80045e+06 6.36481e-02 17.9570 0.160
60 6.02460e+06 4.74049e-03 0.412528 0.176
70 6.02384e+06 9.77666e-11 1.13560e-05 0.189
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 -inf 0.194
Knitro deduced that the problem is non-convex.
78 72 7.26392e+06 FP -inf 2.793
600 508 7.26392e+06 -inf 8.259
608 516 6.38036e+06 MPEC -inf 8.346
906 761 6.37946e+06 MPEC -inf 11.628
957 804 6.22575e+06 MPEC -inf 12.062
1303 1083 6.22575e+06 -inf 15.089
2004 1664 6.22575e+06 -inf 21.782
2052 1705 6.22575e+06 -inf 22.259
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 6.22574779999295250e+06
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 2052 (111.537s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 132351 (3.569s)
# of gradient evaluations = 90666 (1.303s)
# of hessian evaluations = 79562 (3.398s)
# of hessian-vector evaluations = 0
# of subproblems processed = 2181 (128.255s)
Total program time (secs) = 22.26912 (131.613 CPU time)
Time spent in evaluations (secs) = 8.27100
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 117 / 108
Flow-cover cuts = 0 / 0
Probing cuts = 0 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 4 / 1 / 2.888s
Rounding heuristic = 1 / 0 / 0.243s
MPEC heuristic = 35 / 3 / 13.864s
Local search heuristic = 9 / 0 / 0.106s
===========================================================================
WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
- termination condition: maxIterations
- message from solver: Knitro 16.0.0\x3a MIP\x3a Node limit reached.
Integer feasible point found.; objective 6225747.7999929525; optimality
gap Infinity; 2052 nodes; 2181 subproblem solvesShow the full outputHide the full output
draw_design(hanoi, hanoi_design)report_design(hanoi, hanoi_design)Hanoi water network, total pipe cost $6,225,748
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,897.1 250,452 4.0%
3 to 4 40" 7,767.1 320,022 5.1%
4 to 5 40" 7,042.1 403,506 6.5%
5 to 6 40" 6,037.1 125,226 2.0%
6 to 7 40" 4,687.1 236,538 3.8%
7 to 8 40" 4,137.1 236,538 3.8%
8 to 9 40" 3,612.1 222,624 3.6%
9 to 10 30" 2,000.0 171,712 2.8%
10 to 11 24" 1,500.0 155,196 2.5%
11 to 12 24" 940.0 452,655 7.3%
9 to 13 20" 1,087.1 78,712 1.3%
13 to 14 12" 472.1 22,865 0.4%
14 to 15 16" 192.1 38,720 0.6%
16 to 15 12" 270.6 124,843 2.0%
17 to 16 20" 1,135.6 172,182 2.8%
18 to 17 24" 2,480.6 103,464 1.7%
2 to 18 20" 2,540.6 39,356 0.6%
2 to 19 40" 7,762.3 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" 5,072.3 737,442 11.8%
22 to 23 30" 3,455.7 222,322 3.6%
23 to 24 30" 2,635.7 234,975 3.8%
25 to 24 20" -1,117.3 83,632 1.3%
26 to 25 12" -217.3 13,719 0.2%
15 to 26 16" 152.7 52,800 0.8%
22 to 27 16" 571.6 105,600 1.7%
27 to 28 12" 281.6 91,460 1.5%
28 to 29 16" -78.4 112,640 1.8%
29 to 30 16" -438.4 10,560 0.2%
31 to 30 20" 543.4 84,615 1.4%
24 to 31 24" 1,348.4 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.03 57.03 27.03
4 725 0 51.29 51.29 21.29
5 1,005 0 45.26 45.26 15.26
6 1,350 0 43.85 43.85 13.85
7 550 0 42.18 42.18 12.18
8 525 0 40.86 40.86 10.86
9 525 0 39.90 39.90 9.90
10 500 0 38.34 38.34 8.34
11 560 0 34.91 34.91 4.91
12 940 0 30.70 30.70 0.70
13 615 0 36.84 36.84 6.84
14 280 0 31.93 31.93 1.93
15 310 0 31.68 31.68 1.68
16 865 0 41.24 41.24 11.24
17 1,345 0 48.49 48.49 18.49
18 60 0 54.29 54.29 24.29
19 1,275 0 50.69 50.69 20.69
20 930 0 41.35 41.35 11.35
21 485 0 40.08 40.08 10.08
22 1,045 0 44.69 44.69 14.69
23 820 0 39.13 39.13 9.13
24 170 0 35.57 35.57 5.57
25 900 0 32.15 32.15 2.15
26 370 0 31.45 31.45 1.45
27 290 0 39.52 39.52 9.52
28 360 0 31.98 31.98 1.98
29 360 0 32.12 32.12 2.12
30 105 0 32.43 32.43 2.43
31 805 0 33.34 33.34 3.34
The tightest junction is node 12, 0.70 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