# Structure of Optimal Transport Networks Subject to a Global Constraint

Marc Durand

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marc Durand. Structure of Optimal Transport Networks Subject to a Global Constraint. Physical Review Letters, 2007, 98 (8), pp.088701. 10.1103/PhysRevLett.98.088701. hal-00089022v3

HAL Id: hal-00089022

https://hal.science/hal-00089022v3

Submitted on 7 Sep 2010

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Structure of optimal pipe networks subject to a global constraint

Marc Durand

Matière et Systèmes Complexes, UMR 7057 CNRS & Université Paris 7 - Denis Diderot, Tour 33/34 - 2ème étage - case 7056, 2 Place Jussieu - 75251 Paris Cedex 05, France (Date textdate; Received textdate; Revised textdate; Accepted textdate; Published textdate)

## Abstract

The structure of pipe networks minimizing the total energy dissipation rate is studied analytically. Among all the possible pipe networks that can be built with a given total pipe volume (or pipe lateral surface area), the network which minimizes the dissipation rate is shown to be loopless. Furthermore, such an optimal network is shown to contain at most $N-2$ nodes in addition to the $N$ sources plus sinks that it connects. These results are valid whether the possible locations for the additional nodes are chosen freely or from a set of nodes (such as points of a grid). Applications of these results to various physical situations and to the efficient computation of optimal pipe networks are also discussed.

PACS numbers: 89.75.Fb, 05.65.+b, 45.70.Vn, 45.70.Qj

Finding the most efficient transport network is an issue arising in a wide variety of contexts [1, 2]. One can cite, among others, the water, natural gas and power supply of a city, telecommunication networks, rail and road traffic, and more recently the design of labs-on-chips or microfluidic devices. Moreover, this problem also appears in theoretical works intending to describe the architecture of the vascular systems of living organisms [3, 4]. Generally speaking, consider a set of sources and sinks embedded in a two- or three-dimensional space, their respective number and locations being fixed. The flow rates into the network from each source, and out of the network through each sink, are also given. The problem consists in interconnecting the sources and sinks via possible intermediate junctions, referred to as additional nodes, in the most efficient way. That is, to minimize a cost function of general form $\sum_k w_k f(i_k)$, where the summation is over all the links that constitute the network. $w_k$ is the “weight” associated with the $k$th link, and $f$ is some function of the flow rate $i_k$ carried by this link. Minimization of the cost function can be done over different optimization parameters and with different constraints.

Here, the structure of pipe networks that minimize the total dissipation rate $U = \sum_k r_k i_k^2$ is studied, where the weight $r_k$ is the “flow resistance” of pipe $k$, defined as:

$$r_k = \frac{\rho l_k}{s_k^m}, \qquad (1)$$

$\rho$ being some positive constant, $l_k$ and $s_k$ the length and cross-sectional area of each pipe respectively, and $m$ a positive constant characterizing the flow profile. For most flows encountered in physics, $m \ge 1$ (some examples of flows are given later). In this letter, two major results are reported. First, among all the possible pipe networks that can be built with a given value of total pipe volume (or total lateral surface area), the network that minimizes $U$ is loopless. This result suggests an explanation for the observed topologies in the vascular systems of various living organisms [3, 5]. Second, the number of

additional nodes in this optimal network cannot exceed $N-2$, where $N$ is the number of initial nodes (sources plus sinks). As a consequence, the number of possible different topologies for the optimal network is finite.

Flow rates in a network are not independent but must satisfy a conservation law at every source, sink, and additional node. That is, the sum of algebraic flow rates at each site must satisfy:

$$\sum_{\text{adjoining pipes}} i_k = \begin{cases} 0 & \text{at every additional node} \\ I_q & \text{at source or sink } q \end{cases} \qquad (2)$$

where $I_q$ is the fixed inflow/outflow at the source/sink $q$ ($I_q > 0$ for a source, $I_q < 0$ for a sink, and $\sum_{\text{sources}} I_q = -\sum_{\text{sinks}} I_q$). The conservation laws 2 alone do not uniquely determine the flow in each pipe of the network. In many situations $i_k$ also derives from a potential function (electrical potential, pressure, concentration, temperature,...) so that the potential difference $v_k$, the flow rate $i_k$, and the resistance $r_k$ of pipe $k$ are related by Ohm's law $v_k = r_k i_k$. In this case, the flow distribution is unique, and each flow rate $i_k$ is an implicit function of the pipe lengths and pipe cross-sections.

Both the network geometry (characterized by the pipe cross-sections and pipe lengths) and topology (the number of pipes and junctions, and their specific arrangement) can be optimized in order to minimize $U$. However, minimization must be done with some constraint on the pipe cross-sections (otherwise, the optimization problem would be trivial: any network connecting the sources to the sinks with infinitely large pipes would be a solution). Here, a global constraint $C_n = \sum_k l_k s_k^n$ on the total volume ($n=1$) or total surface area ($n=1/2$) of the network is considered. Such a global constraint is less restrictive than the local constraint used in other recent studies on optimal networks [5, 6], where every pipe cross-section is fixed.

Let us now prove that, under the assumptions above, the network that minimizes $U$ is loopless. Let us start
with a network of given topology, whose geometrical parameters (pipe cross-sections and pipe lengths) are adjusted to minimize the dissipation rate (while preserving $C_n$). Indeed, optimization of the network geometry has been studied in a previous work [4]: first, when pipe cross-sections are adjusted, the flow rate $i_k$ carried by each pipe in the network scales with its cross-sectional area $s_k$ as:

$$|i_k| = \kappa I s_k^{(m+n)/2}, \quad (3)$$

where $I$ is the total flow rate through the network ($I = \sum_{sources} I_q = -\sum_{sinks} I_q$), and $\kappa$ is a parameter that depends on $m, n$, and the geometry and topology of the network. Thus, the dissipation rate in this network is:

$$U = (\kappa I)^2 C_n. \quad (4)$$

Then, pipe lengths can also be adjusted in order to minimize $U$, while preserving $C_n$ (according to Eq. 4, this is equivalent to minimizing $\kappa$, which still depends on the pipe lengths). Actually, coordinates of the additional nodes are the appropriate independent optimization parameters. When these coordinates can be freely adjusted, it has been shown [4] that the following vector balance is also satisfied at every additional node of the network with optimized cross-sections and node locations:

$$\sum_k s_k^n e_k = 0, \quad (5)$$

where $e_k$ is the outward-pointing unit vector along each adjoining pipe. No such geometrical rule can be established when the locations of the additional nodes must be chosen from a set of nodes (such as points of a grid, or some particular cities of a country). It must also be noted that Eqs. 3 and 5 are necessary conditions for the minimum of $U$ with respect to the geometrical parameters.

Suppose that this network, which satisfies Eq. 3 (and possibly Eq. 5), contains loops. Let us show that, from this original network, a new loopless network with a lower dissipation rate (and with a same value of $C_n$) can be built. Consider an arbitrary loop in this network. To go from a given junction $A$ to another junction $B$ of this loop, there are two different paths, noted $(\alpha)$ and $(\beta)$, as depicted on Fig. 1a. Let us make a shift of material, in such a way that flows in path $(\alpha)$ tend to be strengthened in one direction (say $A$ to $B$) and flows in path $(\beta)$ tend to be strengthened in the opposite direction ($B$ to $A$). That is, the new cross-sectional areas $s'_k$ in the loop are defined as: $s'_k^{(m+n)/2} = s_k^{(m+n)/2} \pm s_0^{(m+n)/2}$ with, for path $(\alpha)$, a plus sign if flow rate in pipe $(i, j)$ is in direction $A \to B$ and a minus sign if the flow rate is in opposite direction, while signs are inverted for path $(\beta)$ (see Fig. 1b). $s_0$ is a positive number smaller than any cross-sectional area $s_k$ of the original loop. Cross-sections outside the loop

remain unaltered ($s'_k = s_k$). Note that flows in a loop cannot turn all clockwise, or counterclockwise (otherwise, the potential difference $V_A - V_B$ between nodes $A$ and $B$, and thus the flow rates in the loop, would trivially be zero). This guarantees that the cross-sectional areas of the loop did not all simultaneously increase (or decrease).

<figure><img src="image_3.png" /><figcaption>FIG. 1: Shift of material in a loop of the network. (a): the original loop, where flow directions in each pipe are indicated with arrows. (b): the same loop, where the cross-sectional area of a pipe is increased when the direction of its carried flow is $A \to B$ along path $(\alpha)$ or $B \to A$ along path $(\beta)$, and decreased otherwise. The other cross-sectional areas in the network remain unaltered.</figcaption></figure>

Such a variation of cross-sectional areas implies a redistribution of flows in the entire network. Let $\{i'_k\}$ be the new distribution of flow rates satisfying Eq. 2 and Ohm's law, $r'_k = \rho l_k / s_k^{im}$ the new resistances, and $U' = \sum_k r'_k i'_k$ the new dissipation rate. Although we do not know the values of the new flow rates, an upper bound on the new dissipation rate $U'$ can be established, using Thomson's principle [7, 8]. Consider a network of given resistances $r'_k$ that connects the sources to the sinks. Thomson's principle states that – among all possible flow distributions $\{j_k\}$ which satisfy the equations of conservation 2 – the actual flow distribution (i.e.: the one deriving from a potential function and satisfying Ohm's law) is the one that makes the function $\sum_k r'_k j_k^2$ an absolute minimum. Let us consider in particular the flow distribution defined as: $j_k = i_k + i_0$ along path $(\alpha)$, $j_k = i_k - i_0$ along path $(\beta)$, and $j_k = i_k$ for any pipe outside the loop. $i_0$ is some positive number, and $\{i_k\}$ is the actual distribution in the original network, the sign of $i_k$ being (re)defined in both paths as positive if directed from $A$ to $B$. The distribution $\{j_k\}$ satisfies conservation equations 2, since the distribution $\{i_k\}$ does. Besides, by choosing $i_0 = \kappa I s_0^{(m+n)/2}$ and using Eq. 3, the flow rate distribution $j_k$ can be rewritten: $j_k = \operatorname{sgn}(i_k) \kappa I s_k^{(m+n)/2}$. Thus, according to Thomson's principle:

$$U' \le (\kappa I)^2 C'_n, \quad (6)$$

with $C'_n = \sum_k l_k s'_k$. Let us now compare the new value of pipe volume/surface area $C'_n$ with the original value $C_n$. This can be done by studying the variation of $C'_n$ with $s_0$. The derivative of this function with respect to
$x = s_0^{(m+n)/2}$ is:

$$ \frac{\partial C'_{n}}{\partial x} = \frac{2n}{m+n} \sum_{path \ \alpha} l_k \left( s_{k}^{(m+n)/2} + x \right)^{(n-m)/(n+m)} - \frac{2n}{m+n} \sum_{path \ \beta} l_k \left( s_{k}^{(m+n)/2} - x \right)^{(n-m)/(n+m)}. \quad (7) $$

Since $m \ge n$, $\frac{\partial C'_{n}}{\partial x}$ is a decreasing function of $x$, and then:

$$ \frac{\partial C'_{n}}{\partial x} \le \left( \frac{\partial C'_{n}}{\partial x} \right)_{x=0}. \quad (8) $$

Therefore, if the bound in inequality 8 is negative, $C'_{n}$ is a decreasing function of $x$. Using Eqs 1 and 3, this bound can be rewritten as:

$$ \left( \frac{\partial C'_{n}}{\partial x} \right)_{x=0} = \frac{2n}{m+n} \left( \sum_{path \ \alpha} r_k |i_k| - \sum_{path \ \beta} r_k |i_k| \right). \quad (9) $$

We could have chosen to reinforce flows in direction $B \to A$ in path $(\alpha)$, and $A \to B$ in path $(\beta)$ instead, which comes to swapping $(\alpha)$ and $(\beta)$ in the calculations above. Inequalities 6 and 8 would still be satisfied for this new shift of material, but this time with an opposite sign for $\left(\frac{\partial C'_{n}}{\partial x}\right)_{x=0}$ (see Eq. 9). So, necessarily $\left(\frac{\partial C'_{n}}{\partial x}\right)_{x=0} \le 0$ for one of the two shifts, and $C'_{n}$ is a decreasing function of $s_0$ for this particular shift, implying that $C'_{n}(s_0) \le C'_{n}(0) = C_n$. From Eqs. 4 and 6, we obtain that the corresponding dissipation rate $U'$ is also lower: $U'(s_0) \le U$ [9].

In a further step, the total volume/surface area can be increased up to its original value $C_n$ by increasing any cross-sectional areas in the network. This will imply a further decrease in $U$ [10]. Thus, we find a small perturbation of the cross-sections such that the dissipation is reduced for a fixed value of $C_n$ [11]. The reasoning above can be applied with increasingly large values of $s_0$, until eventually one of the pipes in the loop has a zero cross-sectional area, and so one of the paths is cut off. Possible dead branches can be removed, the equivalent material being shifted to the rest of the network by increasing any other cross-sectional areas again, so that the constraint stays at its initial value while the dissipation rate is subjected to a further decrease. Finally, the whole procedure can be repeated to eliminate all the duplicate paths until there are no loops in the network. The argument holds even in case of overlapping loops (that is, loops having pipes in common), and more generally for any topology of the original network. Therefore, it comes that the architecture of the network that minimizes $U$ is loopless. Note that condition 5 is not used throughout the reasoning, so the demonstration is valid whether or not the positions of additional nodes can be freely adjusted.

It must be mentioned that the absence of loops in the least dissipative network has already been conjectured

(without proof) in the particular case of a constrained total volume [12]. A similar result has also been obtained in other studies on optimal networks: Banavar et al. [5] analyzed the flow rate distribution minimizing the cost function $\sum_k w_k |i_k|^\gamma$. They showed that the flow pattern formed by this distribution contains no loops if $0 \le \gamma < 1$, and contains loops if $\gamma > 1$. In that study however, both the network topology and the weight $w_k$ of every link are set. The optimization parameters are the flow rates $i_k$, subject to the conservation laws Eq. 2 only (they do not necessarily derive from a potential function and obey Ohm's law). This optimization problem is then very different than the one analyzed in the present letter. In a different study, Xue et al. [6] showed that the network minimizing the cost function $\sum_k l_k |i_k|^\gamma$ is also loopless when $0 \le \gamma < 1$ [13]. In that study, both the network topology and the pipe lengths are free to be adjusted, but the cross-sections are set to 1; instead of a global constraint on the total pipe volume/surface area, Xue et al. consider a more restrictive constraint on every pipe cross-section. Therefore, the absence of loops in the network minimizing $U$ and preserving $C_n$ cannot be deduced from that study. Indeed, using Eqs. 3 and 4, the dissipation rate of a network with optimized cross-sections can be rewritten in terms of the pipe lengths and the flow rates alone: $U = (\kappa I)^\delta \sum_k l_k |i_k|^\gamma$, where $\delta = 2m/(m+n)$ and $\gamma = 2n/(m+n) \le 1$. The prefactor $\kappa$ in this expression is a parameter dependent on the pipe lengths. Thus, $U$ differs clearly from the cost function studied in [6].

<figure><img src="image_4.png" /><figcaption>FIG. 2: (a): the two adjoining pipes of a two-fold junction carry flows in opposite directions in order to satisfy flow rate conservation. (b): the two adjoining pipes can be favorably replaced with a straight one: since the total pipe length is shortened, the dissipation rate will be decreased for a fixed value of $C_n$.</figcaption></figure>

Let us now show that the number of additional nodes is at most $N - 2$ in the optimal network, where $N$ is the total number of sources plus sinks. This result limits the number of possible topologies for the optimal network. Suppose first that the optimal network is a connected loopless network (or tree). According to Euler's formula [14], a tree has one more node than it has links. So, the number of links in a network with $A$ additional nodes is $N + A - 1$. Since each link has two ends, the number of "incident lines", summed over all the nodes, is $2(N + A - 1)$. This number can be evaluated differently: let $N_p$ be the number of sources or sinks with $p$
incident lines. Since each source or sink is linked to the rest of the tree, the smallest value of $p$ for which $N_p$ has a nonzero value is 1, so $\sum_{p \ge 1} N_p = N$. Similarly, let $A_p$ be the number of additional nodes with $p$ incident lines. By definition, two-fold junctions can exist only if its two links are not parallel. Therefore, such junctions cannot exist in a network satisfying Eq. 5. Two-fold junctions could a priori exist if their positions are chosen from a set of nodes. However, they can be favorably (i.e.: with no increase of $U$ and $C_n$) removed and their two adjoining pipes replaced with a straight one, as depicted in Fig. 2. Thus, the smallest value of $p$ for which $A_p$ is not zero is $p = 3$ in both cases, and the total number of incident lines is: $\sum_{p \ge 1} pN_p + \sum_{p \ge 3} pA_p$. Comparing these two expressions for the number of incident lines, and considering that $\sum_{p \ge 1} pN_p \ge N$ and $\sum_{p \ge 3} pA_p \ge 3A$, it appears that:

$$A \le N - 2, \quad (10)$$

as was to be proven. When both the number of sources and the number of sinks are strictly larger than 1, the optimal network might be disconnected. However, using the reasoning above on each of the trees that constitute the optimal network, it comes that the inequality 10 is still satisfied.

Because of the broad definition of the flow resistance (Eq. 1), the results presented in this letter can be applied in various situations. For instance, the $m = 1$ case corresponds to electrical current in wires, liquid flow in porous conducts, mass or heat diffusion in bars (provided that for the latter, the bar lateral surface is insulated). The $m = 2$ case corresponds to the laminar Poiseuille flow in hollow pipes. Minimization can be done for a fixed lateral surface area ($n = 1/2$) if one wants to save the material required to build the hollow pipes, or for a fixed volume ($n = 1$), if one wants to preserve the amount of liquid flowing through the network. Furthermore, these results may also explain the tree-like structure of the circulatory system of various living organisms [3, 5].

Unfortunately, the results presented in this letter do not give insights into the method of building the optimal network practically, or even into the uniqueness of such an optimal network. In fact, as for the Steiner tree problem - which consists in finding the tree of minimal length interconnecting a set of given points - this problem is likely to be NP-hard, meaning that the solution cannot be found without an exhaustive search of all the possible topologies. However, the NP-hardness does not exclude the possibility of establishing basic properties on the geometry and topology of Steiner trees [15]. Similarly, we were able to address features on the structure of pipe networks minimizing the total dissipation rate under a global constraint. Specifically, the upper bound on the number of additional nodes restricts the number of

possible topologies for the optimal network(s). These results make possible the conception of efficient algorithms for computing the optimal pipe network problem [6]. In many situations however, the capacity of the network to resist random injuries may also play a key role in its design. Obviously, a reticulate network containing redundant paths is more adapted than an arborescent one for that purpose. Therefore, it is sometimes essential to look for a compromise between optimization of flow and robustness of the network.

The author thanks B. Abou, S. Durand, T. Forth and A. Rabodzey for careful reading of the manuscript.

[1] P. S. Stevens, *Patterns in Nature* (Little, Brown, Boston, 1974).

[2] P. Ball, *The Self-Made Tapestry: Pattern Formation in Nature* (Oxford University Press, Oxford, 1998).

[3] T. A. McMahon and J. T. Bonner, *On Size and Life* (Scientific American Library, New York, 1983).

[4] M. Durand, *Phys. Rev. E* **73**, 016116 (2006).

[5] J. R. Banavar, F. Colaiori, A. Flammini, A. Maritan, and A. Rinaldo, *Phys. Rev. Lett.* **84**, 4745-4748 (2000).

[6] G. Xue, T. P. Lillys, and D. E. Dougherty, *SIAM J. Optim* **10**, 22-42 (1999).

[7] P.G. Doyle and J.L. Snell, *Random Walks and Electric Networks*, The Carus Mathematical Monograph, Series 22 (The Mathematical Association of America, USA, 1984); Also in arXiv:math.PR/0001057.

[8] J. H. Jeans, *The mathematical theory of electricity and magnetism*, 5th edition, Cambridge University Press (1925).

[9] Note that if all flows in both paths of the loop point towards same direction - say $A$ to $B$ - in the original network, then $|i_k| = i_k$. The sum of $r_k i_k$ along each path is then equal to $V_A - V_B$, so that $(\frac{\partial C_n}{\partial x})_{x=0} = 0$. Thus, either one of the shifts from $(\alpha)$ to $(\beta)$ or from $(\beta)$ to $(\alpha)$ leads to a decrease in the dissipation rate in that specific case.

[10] M. Durand and D. Weaire, *Phys. Rev. E* **70**, 046125 (2004). An increase in any cross-sectional area implies a decrease in the corresponding individual resistance, and thus a reduction in the total dissipation rate.

[11] Incidentally, this proves that Eq. 3, which is properly a necessary condition for extremum of $U$ with respect to the cross-sections, cannot lead to a minimum of $U$ in a network containing loops.

[12] L. Gosselin and A. Bejan, *International Journal of Thermal Sciences* **44**, 53-63 (2005).

[13] This result holds for any flow distribution satisfying the conservation laws Eq. 2 (the flow rates $i_k$ are not optimization parameters, and they do not necessarily obey Ohm's law).

[14] D. Weaire and N. Rivier, *Contemp. Phys.* **25**, 59-99 (1984).

[15] E.N. Gilbert and H.O. Pollak, *SIAM J. Appl. Math.* **16**, 1-29 (1968).