Designing optimal networks for multicommodity transport problem

Designing and optimizing different flows in networks is a relevant problem in many contexts. While a number of methods have been proposed in the physics and optimal transport literature for the one-commodity case, we lack similar results for the multicommodity scenario. In this paper we present a model based on optimal transport theory for finding optimal multicommodity flow configurations on networks. This model introduces a dynamics that regulates the edge conductivities to achieve, at infinite times, a minimum of a Lyapunov functional given by the sum of a convex transport cost and a concave infrastructure cost. We show that the long-time asymptotics of this dynamics are the solutions of a standard constrained optimization problem that generalizes the one-commodity framework. Our results provide insights into the nature and properties of optimal network topologies. In particular, they show that loops can arise as a consequence of distinguishing different flow types, complementing previous results where loops, in the one-commodity case, were obtained as a consequence of imposing dynamical rules on the sources and sinks or when enforcing robustness to damage. Finally, we provide an efficient implementation of our model which converges faster than standard optimization methods based on gradient descent.

Figure 1

Multicommodity problem illustration for $M=3$. Left: topology of the graph; numbers inside nodes correspond to their indices $v$. Right: an admissible configuration of fluxes. Here, numbers inside nodes correspond to their mass inflow, and we assume each commodity to have its mass concentrated in a single vertex; in the gray node, no mass is entering or exiting, i.e., $S_3^i=0$ for every $i$. Widths of edges are drawn proportional to $F_e^i$, and in the case where all $F_e^i=0$, links are not drawn; arrows denote the direction of the colored fluxes $F_e^i$. Notice how blue mass and orange mass share the same edges, thus creating possible traffic congestion. The inflow-outflow rates of each color are $S^1=(1,0,0,-1)$ (blue), $S^2=(0,3,0,-3)$ (green), and $S^3=(-2,0,0,2)$ (orange).

Figure 2

Toy model where loops are optimal. Here, $M=2$; hence only two colors move though the network. The triangle network has source vectors $S^1=(+1,-1,0),S^2=(-1,+2,-1)$, ${\ell }_a={\ell }_b=1.5,{\ell }_c=1$. The gray patches denote the net loads of each node when ignoring the colors. At the bottom we show one loopy solution on the left and three trees on the right. The green and blue arrows denote the orientation of the loop defined in Appendix pp4 and of the edges, respectively. In detail, the loop has fluxes $F_a=(0,-1),F_b=(0,-1),F_c=(-1,0)$, the leftmost tree has $F_a=(0,0),F_b=(0,-2),F_c=(-1,+1)$ (similarly for the other two). The green star refers to the topology of the toy model used in Fig. 3.

Figure 3

Phase diagram in $\ell$. ${\mathrm{\Gamma }}_{\mathrm{crit}}$ denotes the minimum value of $\mathrm{\Gamma }$ above which loops are optimal. The setup is the same as for the toy model in Fig. 2. Values of ${\mathrm{\Gamma }}_{\mathrm{crit}}$ are found by solving Eq. (13); ${\ell }_c=1$. The area under the triangular surface with white background is not allowed as the triangular geometry is not defined there. One can notice that there is an entire region where ${\mathrm{\Gamma }}_{\mathrm{crit}}\le 1$; inside it, loops are optimal. On the right-hand side of the figure we show three different topologies, i.e., choices of ${\ell }_a,{\ell }_b$ for which optimal solutions can be loopy; these are associated with the markers drawn on the heat map, and the green star is the configuration of the toy model in Fig. 2. In particular, running our dynamics on the pink-diamond graph (respectively, the yellow-plus graph) leads to a loopy configuration for $\mathrm{\Gamma }>{\mathrm{\Gamma }}_{\mathrm{crit}}$ or to a tree if $\mathrm{\Gamma }<{\mathrm{\Gamma }}_{\mathrm{crit}}$. Running our dynamics on the blue-cross graph returns always since its ${\mathrm{\Gamma }}_{\mathrm{crit}}$ is equal to 1. Widths of the edges are proportional to the final $\left|\right|F_e{\left|\right|}_2$, and edges that are not visible have negligible fluxes; nodes' dimensions are proportional to their inflowing mass. In the bottom right portion of the panel we report the values of ${\mathrm{\Gamma }}_{\mathrm{crit}}$ obtained solving Eq. (13) fixing ${\ell }_a,{\ell }_b$ as given by the markers.

Figure 4

Interaction between commodities. The figure shows how changing the load of one commodity influences the path taken by others. (a) We plot the heat map of ${\mathrm{\Gamma }}_{\mathrm{crit}}$ obtained using Eq. (13), with $S_1^1=-S_3^1=+1$ and for different configurations of the purple commodity ($i=2$): (i) $S_2^2=+2$, $S_3^2=S_1^2=-1$; (ii) $S_2^2=+4$, $S_3^2=S_1^2=-2$. The areas under the white surface correspond to regions where node 1 is a sink, respectively, nodes 2 and 3 are sources, for the purple mass. The green and the red triangular markers denote the configurations we discuss in Sec. 3c and in the rest of the panel. (b) Optimal graphs for the green and the red triangle, fixing $\mathrm{\Gamma }=0.6$ and $\mathrm{\Gamma }=0.8$. The width of each edge $e$ is proportional to $|F_e^i|$ for each color. (c) One-commodity solutions obtained injecting only the yellow or the purple mass in the triangle network.

Figure 5

Phase diagram in $S$. This can be numerically found by fixing $\left\{{\ell }_e\right\}$ as in Fig. 2, $S_2^1=-1-S_3^1,S_1^2=-2-S_3^2$, while varying $S_3^1$ and $S_3^2$. The color bar denotes the value of ${\mathrm{\Gamma }}_{\mathrm{crit}}$ above which optimal solutions can be loopy. The green star denotes the configuration of $S$ used in Fig. 2.

Figure 6

Sketch of the trimming procedure described in Appendix pp4, on the toy model of Fig. 2. The tree obtained in the colored case is not optimal, while the one obtained in the colorless case has lower cost but violates Kirchhoff's law.

Figure 7

Computational complexity. (a) Running time (in seconds) as a function of system size $N$. (b) Running time as a function of the ratio $M/N$ between the number of commodities and system size. GD denotes gradient descent, implemented with Eqs. (E3) and (E4); we only show it for $\beta =0.5$ as for the other values it fails to converge within a reasonable time. Similarly, for $\beta =1$ and $M/N<1$, Optimization fails to converge, and hence we only report Dynamics.