Distributed Combined Space Partitioning and Network Flow Optimization: an Optimal Transport Approach (Extended Version)

TL;DR

This paper introduces a novel approach combining space partitioning and network flow optimization using semi-discrete optimal transport, providing a distributed algorithm with proven duality properties, applicable to large-scale systems like power and communication networks.

Contribution

It formulates a coupled SDOT and flow optimization problem, proves strong duality, and develops a distributed dual gradient algorithm for large-scale network applications.

Findings

Algorithm performs well in simulations.

Applicable to power distribution network reconfiguration.

Provides theoretical guarantees via duality analysis.

Abstract

This paper studies a combined space partitioning and network flow optimization problem, with applications to large-scale power, transportation, or communication systems. In dense wireless networks, one may want to simultaneously optimize the assignment of many spatially distributed users to base stations and route the resulting communication traffic through the backbone network. We formulate the overall problem by coupling a semi-discrete optimal transport (SDOT) problem, capturing the space partitioning component, with a minimum-cost flow problem on a discrete network. This formulation jointly optimizes the assignment of a continuous demand distribution to certain endpoint network nodes and the routing of flows over the network to serve the demand, under capacity constraints. As for SDOT problems, we show that the formulation of our problem admits a tight relaxation taking the form of an infinite-dimensional linear program, derive a finite-dimensional dual problem, and show that strong duality holds. We leverage these results to design a distributed dual gradient ascent algorithm to solve the problem, where nodes in the graph perform computations based solely on locally available information. Simulation results illustrate the algorithm performance and its applicability to an electric power distribution network reconfiguration problem. This version extends the CDC 2025 conference paper with additional proof sketches.

Figures

Equations

$$\underset{\begin{array}{c}T:X\to \mathcal{S}\\ p\in {\mathbb{R}}^{|\mathcal{A}|}\end{array}}{\inf }$$

$$\sum _{(i,j)\in \mathcal{A}}{p}_{ij}-\sum _{(j,i)\in \mathcal{A}}{p}_{ji}=s_i-(T_{\#}\mu {)}_i,\forall i\in \mathcal{S},$$

$$\sum _{(i,j)\in \mathcal{A}}{p}_{ij}-\sum _{(j,i)\in \mathcal{A}}{p}_{ji}=s_i,\forall i\in \mathcal{N}\setminus \mathcal{S},$$

$$a_{ij}\le p_{ij}\le b_{ij},\forall (i,j)\in \mathcal{A},$$

$$\underset{\begin{array}{c}\pi \in {\mathcal{M}}_{+}(X\times \mathcal{S})\\ p\in {\mathbb{R}}^{|\mathcal{A}|}\end{array}}{\inf }$$

$$\sum _{(i,j)\in \mathcal{A}}{p}_{ij}-\sum _{(j,i)\in \mathcal{A}}{p}_{ji}=s_i-\pi (X,i),\forall i\in \mathcal{S},$$

$$\sum _{(i,j)\in \mathcal{A}}{p}_{ij}-\sum _{(j,i)\in \mathcal{A}}{p}_{ji}=s_i,\forall i\in \mathcal{N}\setminus \mathcal{S},$$

$$\sum _{i\in \mathcal{S}}\pi (A,i)=\mu (A),\forall A\subset \mathcal{B}(X),$$

$$a_{ij}\le p_{ij}\le b_{ij},\forall (i,j)\in \mathcal{A},$$

$$\underset{\psi \in {\mathbb{R}}^{|\mathcal{N}|}}{\sup }q(\psi ),$$

$$q(\psi )$$

$$+\sum _{(i,j)\in \mathcal{A}}\underset{p\in [a_{ij},b_{ij}]}{\min }\{c_{ij}(p)-({\psi }_i-{\psi }_j)p\}.$$

$${\mathrm{Lag}}_i(\psi )=\left\{x\in X:c(x,i)-{\psi }_i\le c(x,j)-{\psi }_j,\forall j\in \mathcal{S}\right\}.$$

$$p_{ij}(\psi ):=\arg \underset{p\in [a_{ij},b_{ij}]}{\min }\{c_{ij}(p)-({\psi }_i-{\psi }_j)p\},$$

$$g(\psi {)}_i=s_i-\mu ({\mathrm{Lag}}_i(\psi ))-\sum _{(i,j)\in \mathcal{A}}{p}_{ij}(\psi )+\sum _{(j,i)\in \mathcal{A}}{p}_{ji}(\psi ),$$

$$g(\psi {)}_i=s_i-\sum _{(i,j)\in \mathcal{A}}{p}_{ij}(\psi )+\sum _{(j,i)\in \mathcal{A}}{p}_{ji}(\psi ).$$

$$q_{\text{SDOT}}(\psi )$$

$$q_{\text{MCF}}(\psi )=\sum _{i\in \mathcal{N}}{\psi }_i{s}_i+\sum _{(i,j)\in \mathcal{A}}\underset{p\in [a_{ij},b_{ij}]}{\min }\{c_{ij}(p)-({\psi }_i-{\psi }_j)p\},$$

$$q_{\text{SDOT}}(\psi )=\sum _{i\in \mathcal{S}}{\int }_{{\mathrm{Lag}}_i(\psi )}(c(x,i)-{\psi }_i)\mathrm{d}\mu (x).$$

$${\psi }^{k+1}={\psi }^k+{\gamma }_{k}g({\psi }^k),$$

$$\sum _{k=0}^{\infty }{\gamma }_k=+\infty \text{and}\sum _{k=0}^{\infty }{\gamma }_k^2<+\infty .$$

$$p_{ij}\leftarrow \arg \underset{p\in [a_{ij},b_{ij}]}{\min }\{c_{ij}(p)-({\psi }_i-{\psi }_j)p\}$$

$$c(x,y)+c(y,i)-{\psi }_i$$

$$\le c(x,j)-{\psi }_j$$

$$\le c(x,y)+c(y,j)-{\psi }_j,$$

$$\inf_{z\in E}\Bigl{\{}\Theta(z)+\Xi(z)\Bigr{\}}=-\max_{z^{*}\in E^{*}}\Bigl{\{}-\Theta^{*}(-z^{*})-\Xi^{*}(z^{*})\Bigr{\}}.$$

$$⟨(u,v),(\pi ,p)⟩={\int }_{X\times \mathcal{S}}u(x,i)\mathrm{d}\pi (x,i)+\sum _{(i,j)\in \mathcal{A}}{v}_{ij}{p}_{ij}.$$

$${\Theta }_{\text{OT}}(u)=\{\begin{array}{ll}0,& \text{if }u(x,i)\ge -c(x,i)\forall (x,i)\in X\times \mathcal{S},\\ +\infty ,& \text{otherwise,}\end{array}$$

$$\Xi (u,v)={\int }_X\phi (x)\mathrm{d}\mu (x)+\sum _{i\in \mathcal{N}}{\psi }_i{g}_i,$$

$$\underset{(u,v)\in E}{\inf }\{\Theta (u,v)+\Xi (u,v)\}=-\underset{\psi \in {\mathbb{R}}^{|N|}}{\sup }q(\psi ).$$

$$\Theta^{*}\bigl{(}-(\pi,p)\bigr{)}=\int_{X\times\mathcal{S}}c(x,i)\,\mathrm{d}\pi(x,i)+\sum_{(i,j)\in\mathcal{A}}c_{ij}(p_{ij}),$$

$$q({\psi }^{\ast })=$$

$$+\sum _{i\in \mathcal{S}}{\psi }_i^{\ast }\left[-\sum _{(i,j)\in \mathcal{A}}{p}_{ij}^{\ast }+\sum _{(j,i)\in \mathcal{A}}{p}_{ji}^{\ast }+s_i-{\pi }^{\ast }(X,i)\right]$$

$$+\sum _{i\in \mathcal{N}\setminus \mathcal{S}}{\psi }_i^{\ast }\left[-\sum _{(i,j)\in \mathcal{A}}{p}_{ij}^{\ast }+\sum _{(j,i)\in \mathcal{A}}{p}_{ji}^{\ast }+s_i\right].$$

$$B_{ij}:=$$

$$-\underset{p\in [a_{ij},b_{ij}]}{\min }\{c_{ij}(p)-({\psi }_i^{\ast }-{\psi }_j^{\ast })p\},$$

$$\left\{(x,i)\in X\times \mathcal{S}:c(x,i)-{\psi }_i^{\ast }=\underset{k\in \mathcal{S}}{\min }\{c(x,k)-{\psi }_k^{\ast }\}\right\}.$$

Full text

I Introduction

Optimizing large-scale networks requires solving two fundamental problems: allocating end resources to serve the demand for a product or service, and efficiently routing flows through the network to ultimately meet that demand. For instance, in electric power systems, consumers are assigned to a substation of the distribution system and the necessary power flows are routed through the transmission network to feed the substations. Similarly, in cellular communication networks, users must be matched to base stations and the resulting traffic must be routed through the backbone network. Although the assignment and routing problems are often treated separately for computational efficiency, considering both jointly can improve performance [1, 2].

This work addresses such a combined assignment and flow optimization problem, from an optimal transport (OT) point of view [3, 4]. Specifically, we couple a semi-discrete optimal transport (SDOT) problem  [3, Chapter 5] with a minimum-cost flow (MCF) problem [5], and develop a distributed algorithm to solve the combined problem. In the standard SDOT problem, one seeks an optimal transport map from an arbitrary source measure, representing for instance a demand distribution over a geographic area, to a discrete target measure, e.g., a set of endpoint nodes capable of serving this demand. This map in effect partitions the space supporting the demand distribution into different cells, such that the total demand in each cell equals the capacity of the associated endpoint node. Meanwhile, the MCF problem allows us to determine how to route flows efficiently through a network connecting the endpoints to supply nodes, subject to edge capacity constraints.

Problems related to the one formulated in this paper have been considered in several application domains. For example, in electrical distribution networks comprising substations and consumers connected via lines equipped with tie and sectionalizing switches, the radial reconfiguration problem [6], further discussed in Section V-B, aims to selectively open and close switches to form a radial topology, with substations acting both as root nodes of the distribution network and as interfaces to the transmission network where generation is available. This reconfiguration process results in a spatial partitioning that guarantees demand satisfaction while minimizing power losses. Conventional approaches to address it rely on solving computationally difficult mixed-integer programs or use faster heuristics that yield suboptimal solutions [6]. Recent studies have explored clustering [7, 8] and partitioning techniques [9] to mitigate computational challenges, but retain a purely discrete network modelling approach. In contrast, the SDOT framework adopted here can provide good heuristic solutions for asymptotically large distribution networks by considering the limit of a continuous distribution of customers.

In communication networks, related problems combining user association and resource allocation have been formulated as difficult combinatorial problems [2, 10]. Alternatively, SDOT has been used for space partitioning, e.g., for C-RAN device association [11] and communication with UAVs[12]. However, these papers did not integrate network flow optimization. Another related problem is branched optimal transport (BOT) [13], which uses subadditive costs in an OT formulation to encourage mass transport along paths, effectively yielding branched networks. However, in BOT, the backbone network needs to be designed, whereas it is fixed and given here. Moreover, exact methods are limited to small instances and heuristics are generally required [13].

To ensure scalability to large-scale networks, we propose a distributed optimization method to solve our problem, where an agent placed at each node of the backbone network updates its local variables by communicating only with neighbouring agents. Distributed methods for various OT problems have been proposed in recent years. Decentralized alternating direction methods of multipliers (ADMM)-based methods can be used in discrete settings [14], while [15] propose a distributed online optimization and control strategy for a continuous problem. For SDOT, dual deterministic and stochastic gradient ascent methods can be implemented in a distributed manner by computing assignment cell areas exactly or via sampling methods [16, 17]. MCF problems have similarly been solved using distributed auction-based algorithms [18] or distributed dual gradient ascent when cost functions are strictly convex [5]. Under additional regularity conditions, [19] develops a distributed Newton method. The approach proposed here relies on a dual gradient ascent to accommodate both the SDOT and MCF components.

The main contributions of this work are threefold. First, we propose a new framework combining SDOT and MCF, enabling joint space partitioning and flow optimization for the backbone network. Specifically, our model extends SDOT by allowing the target distribution to vary according to network constraints, a significant departure from standard formulations where this distribution is fixed. Second, leveraging optimal transport theory, we derive the dual problem and show that strong duality holds. Third, we develop a distributed dual ascent algorithm to solve the combined SDOT–MCF problem, suitable for large-scale networks.

The rest of this paper is organized as follows. Section II formulates the problem. Section III-A introduces a Kantorovich-type relaxation [4] of the problem, derives its dual, shows that strong duality holds, and that the relaxation is tight. Section IV develops the distributed dual gradient ascent algorithm. Simulation results are presented in Section V, and concluding remarks are provided in Section VI.

II Problem formulation

We now formally define our optimization problem. We consider a space $X$ equipped with a non-negative measure $\mu$, which models, for instance, a demand distribution for a certain product, e.g., electric power, by a population of consumers spatially distributed over $X$. In addition, we consider a network modelled by a directed graph $G=(\mathcal{N},\mathcal{A})$, where $\mathcal{N}\subset\mathbb{N}$ is a finite set of nodes and $\mathcal{A}\subset\mathcal{N}\times\mathcal{N}$ is a set of arcs. Each node $i\in\mathcal{N}$ has an associated net supply value $s_{i}\in\mathbb{R}$ (with $s_{i}<0$ allowed, indicating a demand node). A subset of the nodes, denoted by $\mathcal{S}\subset\mathcal{N}$, corresponds to endpoint stations effectively serving the consumers, e.g., distribution substations in power systems. Serving a unit of demand at $x\in X$ by endpoint $i\in\mathcal{S}$ incurs a cost $c(x,i)$, where $c:X\times\mathcal{N}\to\mathbb{R}\cup\{+\infty\}$ is a given function, taking possibly infinite values. Moreover, the product can travel through the network, e.g., the electric transmission system, with the variables $p=(p_{ij})_{(i,j)\in\mathcal{A}}\in\mathbb{R}^{|\mathcal{A}|}$ representing the flows along the arcs, where $|\cdot|$ denotes the cardinality of a set. Sending a flow $p_{ij}$ along arc $(i,j)\in\mathcal{A}$ incurs a cost $c_{ij}(p_{ij})$, e.g., corresponding to active power losses, for given functions $c_{ij}:\mathbb{R}\to\mathbb{R}\cup\{+\infty\}$. Illustrative examples are shown on Figures 1(a) and 3 in Section V.

The overall objective is to simultaneously determine how to assign each consumer of $X$ to an endpoint in $\mathcal{S}$ and route the product through the network, so as to satisfy demand and minimize the combined costs of assignment and routing. A deterministic assignment of consumers in $X$ to endpoints in $\mathcal{S}$ can be described by a map $T:X\to\mathcal{S}$. The total demand assigned to endpoint $i\in\mathcal{S}$ is then given by $(T_{\#}\mu)_{i}=\mu(\{x\in X:T(x)=i\})=\mu(T^{-1}(i))$, where $T_{\#}\mu$ denotes the pushforward measure of $\mu$ by $T$. The set $T^{-1}(i)\subset X$ defines a cell of points assigned to endpoint $i$, and these cells for $i\in\mathcal{S}$ form a partition of $X$. The combined spatial assignment and network flow optimization problem is then formulated as follows

[TABLE]

where in (1d) the scalars $a_{ij}$ and $b_{ij}$, for $(i,j)\in\mathcal{A}$, are lower and upper limits on arc flows, and (1b) and (1c) represent the conservation of flow at each node. In particular, we distinguish between the flow constraints (1b) at the endpoints $i\in\mathcal{S}$, which must serve the demand $(T_{\#}\mu)_{i}$, and the flow constraints (1c) at the other nodes that only route traffic through the network.

In this paper, we make the following assumptions.

Assumption 1**.**

We assume that:

1. (a)

The assignment cost function $c:X\times\mathcal{S}\to\mathbb{R}\cup\{+\infty\}$ is lower semi-continuous and bounded from below. 2. (b)

The domain $X$ is compact. 3. (c)

The arc cost functions $c_{ij}:\mathbb{R}\to\mathbb{R}\cup\{+\infty\}$ are convex and lower semi-continuous for all $(i,j)\in\mathcal{A}$. 4. (d)

Problem (1) admits a feasible solution.

In Assumption 1, (a) and (b) are standard in the OT literature, with (b) being satisfied in practical applications and useful to simplify technical arguments [4]. Property (c) is necessary to leverage duality results and also standard for MCF problems. Assumption 1-(d) is non-trivial and implies in particular, by summing constraints (1b) and (1c), that we must have $\mu(X)=\sum_{i\in\mathcal{N}}s_{i}$, i.e., an equilibrium between demand and supply, in order for the flow conservation constraints to admit a feasible solution.

The objective of this work is to develop an algorithm to solve the optimization problem (1), which moreover admits a distributed implementation by the nodes of the network, i.e., each node should execute operations that only require information exchanges with their neighbours in the graph $G$. In addition, the method of assigning consumers to endpoints should also scale to large-scale problems.

III Characterization of an Optimal Solution

III-A Kantorovich Relaxation and its Dual

From an OT point of view, (1) is a type of Monge problem [4], because the assignment of consumers to endpoints takes the form of a map. For such problems, it is generally useful to introduce the corresponding Kantorovich relaxation, by replacing the transport map $T$ with a transport plan $\pi\in\mathcal{M}_{+}(X\times\mathcal{S})$, i.e., a non-negative measure on $X\times\mathcal{S}$. Intuitively, this change allows for randomized assignments of consumers to endpoints. The relaxed problem reads

[TABLE]

where $\mathcal{B}(X)$ denotes the $\mu$-measurable subsets of $X$. Constraint (2c) ensures that $\mu$ is the first marginal of $\pi$, and in (2b) the quantity $\pi(X,i)$ represents again the total demand assigned to $i$. Problem (2) is a relaxation of (1), because (1) is obtained by restricting $\pi$ to measures induced by deterministic maps, i.e., of the form $\pi=(\text{id}_{X},T)_{\#}\mu$, where $\text{id}_{X}$ is the identity map of $X$ and $T:X\to\mathcal{S}$.

The benefit of this relaxation is that (2) is now a linear program, although still infinite-dimensional in general, which satisfies useful duality properties [4]. In particular, we establish the following key duality result, whose proof is sketched in Appendix -A.

Theorem 1**.**

Under Assumption 1, the optimal value of (2) is equal to

[TABLE]

where the dual function $q:\mathbb{R}^{|\mathcal{N}|}\to\mathbb{R}$ is given by

[TABLE]

Moreover, the infimum in (2) is attained at an optimal solution $(\pi^{*},p^{*})\in\mathcal{M}_{+}(X\times\mathcal{S})\times\mathbb{R}^{|\mathcal{A}|}$.

Note that (3) is a finite-dimensional optimization problem, in contrast to the primal problem (2), a feature that we exploit to design our algorithm in Section IV.

III-B Reconstructing the Primal Optimal Solution

To recover an optimal solution for the original problem (1), we introduce a few additional assumptions. The first is commonly assumed to guarantee the tightness of the Kantorovich relaxation for SDOT problems [16].

Assumption 2**.**

For every pair of distinct nodes $i,j\in\mathcal{S}$ and every real number $r\in\mathbb{R}$, the set $\{x\in X:c(x,i)-c(x,j)=r\}$ has $\mu$-measure zero.

Define, for each $i\in\mathcal{S}$ and $\psi\in\mathbb{R}^{|\mathcal{N}|}$, the generalized Laguerre cell [20, 3]

[TABLE]

These cells are used below to define the optimal assignment map. Assumption 2 ensures that the intersection of two distinct Laguerre cells, where assignment randomization could be beneficial, has $\mu$-measure [math] and hence does not contribute to the overall cost. Note that the Laguerre cells are polytopes when $X=\mathbb{R}^{d}$ and $x\mapsto c(x,i)$ is the squared Euclidean distance between $x$ and the position of node $i\in\mathcal{S}$ [20]. Next, we introduce the following technical assumption.

Assumption 3**.**

A dual optimal solution $\psi^{*}\in\mathbb{R}^{|\mathcal{N}|}$ maximizing (3) exists.

Explicit conditions under which Assumption 3 is satisfied will be developed in future work. It is empirically satisfied in our numerical simulations in Section V, and it is known to be satisfied for the standard SDOT problem under Assumption 1-(b), see [21, Lemma 9]. Finally, the next assumption, which strengthens Assumption 1, is often made to simplify the analysis of the MCF component of the problem.

Assumption 4**.**

For each arc $(i,j)\in\mathcal{A}$, the cost function $c_{ij}:\mathbb{R}\to\mathbb{R}$ is strictly convex.

Assumption 4 ensures that for all $(i,j)\in\mathcal{A}$, the minimizer

[TABLE]

exists and is unique. The following proposition provides a means to compute an optimal solution for the original problem (1), assuming an optimal dual solution is known. It also shows that the relaxation (2) is tight.

Proposition 1**.**

Under Assumptions 1, 2, 3, and 4, let $\psi^{*}$ be a maximizer of (3), let $p^{*}_{ij}\coloneqq p_{ij}(\psi^{*})$ for each arc $(i,j)\in\mathcal{A}$, computed from (5), and let $p^{*}\in\mathbb{R}^{|\mathcal{A}|}$ be the vector with components $p^{*}_{ij}$. Define the assignment map $T^{*}(x):=\arg\min_{i\in\mathcal{S}}\{c(x,i)-\psi^{*}_{i}\}$, for every $x\in X$, with ties between endpoints in the minimization broken arbitrarily, if any. Then, the pair $(T^{*},p^{*})$ is an optimal solution for (1).

Note that the map $T^{*}$ in Proposition 1 specifies that all points in $\mathrm{Lag}_{i}(\psi^{*})$ should be assigned to endpoint $i\in\mathcal{S}$. The proof of this proposition is sketched in Appendix -B.

IV Distributed Algorithm

In this section we propose a distributed algorithm to solve problem (1). Based on the previous analysis, one can compute an optimal dual solution for (3) and then recover an optimal primal solution using Proposition 1.

IV-A Concavity and Supergradient of the Dual Function

In this section we show formally that the dual function (4) is concave, which is expected from duality theory. We also provide an explicit expression for its supergradient, which we exploit later to design a gradient ascent algorithm.

Proposition 2**.**

The dual function $q$ defined in (4) is concave. Moreover, under Assumptions 1, 2, and 4, a supergradient $g(\psi)$ at $\psi\in\mathbb{R}^{|\mathcal{N}|}$ has components, for each $i\in\mathcal{S}$,

[TABLE]

and, for each $i\in\mathcal{N}\setminus\mathcal{S}$,

[TABLE]

Proof.

We decompose the dual function $q$ in (4) as $q(\psi)=q_{\text{SDOT}}(\psi)+q_{\text{MCF}}(\psi)$, with

[TABLE]

and,

[TABLE]

by analogy with the dual functions for the SDOT problem [16] (without the linear term $\sum_{i\in\mathcal{S}}\psi_{i}s_{i}$) and the MCF problem [5]. Then, $q_{\text{MCF}}$ is concave as a sum of linear terms and the minimum of affine functions. Under Assumption 4, its supergradient at $\psi$ is given by the expression (7) for each $i\in\mathcal{N}$, see [18, Chapter 5.5.]. Moreover, under Assumption 2, we can write

[TABLE]

As shown in [16, Theorem 4], this function is concave and its supergradient at $\psi$ has as component $-\mu(\mathrm{Lag}_{i}(\psi))$ for each $i\in\mathcal{S}$. By addition, $q$ is concave and its supergradient is given by (6) for each $i\in\mathcal{S}$ and  (7) for each $i\in\mathcal{N}\setminus S$. ∎

IV-B Supergradient Algorithm and Distributed Implementation

Starting from an arbitrary initial value $\psi^{0}\in\mathbb{R}^{|\mathcal{N}|}$, an iterative supergradient ascent algorithm to maximize (3) takes the form

[TABLE]

with $g(\psi)$ given in Proposition 2 and the positive scalar stepsizes $\{\gamma_{k}\}_{k\geq 0}$ satisfying the standard conditions

[TABLE]

Under these conditions, and noting that the gradients in Proposition 2 are uniformly bounded, the sequence $q(\psi^{k})$ converges to an optimal value $q^{*}$ of (3) as $k\to+\infty$. Moreover, under Assumption 3, the iterates $\psi_{k}$ also converge to an optimal solution $\psi^{*}$ [22, Proposition 8.2.6], which can then be used to compute the assignments and flows using Proposition 1. We can now describe a distributed method to solve problem (1), outlined in Algorithm 1.

Assume that a computing agent at each node $i\in\mathcal{N}$ within the network updates the dual variable $\psi_{i}$ according to (8), for which it needs to compute the supergradient component $g_{i}(\psi)$ given in Proposition 2. At each iteration the agents compute the current flows $\{p_{ij}\}_{(i,j)\in\mathcal{A}}$ on Lines 5-7 by exchanging information only with their outgoing and incoming neighbours. For endpoint nodes $i\in\mathcal{S}$, an additional step on Line 9 is required to determine the measure $m_{i}=\mu(\mathrm{Lag}_{i}(\psi))$ of its current Laguerre cell in order to compute the supergradient component (6). This step depends a priori on the variables $\psi_{j}$ of possibly all the other endpoints $j\in\mathcal{S}$. Most practical instances feature a finite set of customers and a discrete measure $\mu$ encoding their individual demands. At iteration $k$, each endpoint $i\in\mathcal{S}$ needs to compare the adjusted costs $c(x,i)-\psi_{i}^{k}$ with those of the other endpoints. A customer at $x$ is matched to the endpoint that offers the lowest value. Then, we can compute $m_{i}=\sum_{x\in\mathrm{Lag}_{i}(\psi)}\mu(x)$ the total demand of the customers lying in the Laguerre cell of $i$. For some cost functions $c(x,i)$, e.g., the squared distance between $x$ and the position of node $i$, the geometric properties of the Laguerre cells can be exploited so that an endpoint only needs to compare the adjusted costs with other endpoints that share a cell boundary with it, see, e.g., [23].

Remark 1**.**

The discretization of the space $X$ to compute the masses $m_{i}$ deterministically can be replaced by a stochastic integration method, generating samples according to $\mu$ for which the endpoints then compare their adjusted costs, as discussed in [16] for SDOT. The iterates (8) then become a stochastic gradient algorithm, again with convergence guarantees.

Finally, the agents need to detect convergence in a distributed manner. A possible test used on Line 15 is that each node stops updating its dual variable $\psi_{i}$ when its variation stays below a certain threshold $\varepsilon$. Each agent may also wait until the variables of its neighbours have converged as well before recomputing the final flows (as on Line 6) and Laguerre cells. Regarding the assignment computation, each endpoint can determine the points $x\in X$ belonging to its Laguerre cell. Alternatively in some applications, the final weights $\psi^{k}$ can be broadcast to the customers, who can then determine their assigned endpoint automatically.

V Numerical Simulations

In this section, we illustrate our method through two numerical examples. First, we test Algorithm 1 on a synthetic example. Second, we briefly explore the applicability of the method to a more complex scenario relevant to electric power distribution networks.

V-A Synthetic Example

We first consider a simple scenario shown on Figure 1(a), defined on a square domain $X\subset\mathbb{R}^{2}$ of side length $L=100$. A continuous consumer demand distribution is discretized on a $200\times 200$ grid, for a total of 40,000 points. Each consumer is assigned a mass from a truncated Gaussian density with mean $(50,75)$ and standard deviation $25$. The backbone network consists of two nodes $\{S_{1},S_{2}\}$ with a supply $s_{i}=0.5$, called source nodes, and four nodes $\{I_{1},I_{2},I_{3},I_{4}\}$ with zero supply, called interconnection nodes. Node positions and arcs are also shown on Figure 1(a).

The assignment cost $c(x,i)$ between consumer demand at $x$ and node $i$ is given by the Euclidean distance between $x$ and the position of node $i$. The arc costs are assumed to be quadratic functions of the flows, i.e., $c_{ij}(p)\coloneqq d_{ij}\,p^{2}$, where $d_{ij}$ is the Euclidean distance between nodes $i$ and $j$. We assume $a_{ij}=-1$ and $b_{ij}=1$ for all $(i,j)\in\mathcal{A}$, and obtain the closed-form local flow update at Line 6 of Algorithm 1 given by $p_{ij}=\operatorname{proj}_{[-1,1]}\left(\frac{\psi_{i}-\psi_{j}}{2\,d_{ij}}\right)$. The algorithm is run for $300$ iterations with diminishing step-size $\gamma_{k}=1/(1+0.01k)$. Figure 2 shows the convergence of the components of the supergradient $g$ toward zero and of the dual variables $\psi$. Figure 1(b) illustrates the resulting routing and partitioning solutions. Each cell is coloured according to its assigned endpoint node, and the network is represented with annotated flows on its arcs. Source nodes are labelled with their fixed supply $0.5$, and endpoints display their assigned demand

V-B Electric Power Network Example

We now discuss an application to an electrical power network consisting of a transmission system linking generators to substations and a distribution network connecting substations to consumers. Consumer demand, initially discrete, is approximated by a continuous distribution. Substations are endpoint nodes $\mathcal{S}$ in our formulation.

An important problem in electrical networks is the radial reconfiguration of the distribution network: activating or deactivating switches on lines to ensure that each consumer is served via exactly one distribution path to a unique substation, forming tree-like (radial) subnetworks rooted at substations, see Figure 3. Solving large-scale reconfiguration problems exactly is computationally challenging [6], so that the coupling with the transmission network optimization problem is usually neglected. Here we use Algorithm 1 to develop a heuristic producing an initial assignment of customers to substations while minimizing transmission costs. This enables subsequent parallelized radial reconfiguration per substation, greatly reducing the computational burden.

We assume that the arc costs on the transmission network are $c_{ij}(p)=r_{ij}\,p^{2}$, representing power losses, where $r_{ij}$ is the line resistance and $p$ is the power flow. We take $r_{ij}$ as the distance $d_{ij}$ in this example. Next, to approximately capture the power losses in the distribution network through the simple assignment cost $\int c(x,T(x))\,\mathrm{d}\mu(x)$, we choose for $c(x,i)$ the geodesic (shortest‑path) resistance through this network between a customer $x$ and the substation $i$. The quality of this heuristic choice will be explored in future work. This choice also ensures an essential connectivity property for the resulting Laguerre cells: all consumers assigned to a substation can be physically connected to it via paths lying entirely within their assigned cell. Briefly, if a consumer $x$ belongs to the Laguerre cell $\mathrm{Lag}_{i}(\psi)$ associated with substation $i$ and vector $\psi\in\mathbb{R}^{|\mathcal{S}|}$, all intermediate nodes $y$ along the shortest path from $x$ to $i$ also belongs to $\mathrm{Lag}_{i}(\psi)$. This follows directly from the definitions of Laguerre cells and shortest paths, because for all $j\in\mathcal{S}$:

[TABLE]

which implies $c(y,i)-\psi_{i}\leq c(y,j)-\psi_{j}$, hence $y\in\mathrm{Lag}_{i}(\psi)$.

Figure 3 illustrates the partitioning and routing solutions obtained for this power network example. The transmission network topology mirrors the synthetic example described in Section V-A. The distribution network was generated by randomly placing 1000 consumer nodes within the domain, each assigned a uniformly random demand value. To ensure a realistic network structure, each consumer node was connected to its two or three geographically closest neighbours. Figure 3 displays the resulting Laguerre cells, with consumer nodes coloured according to their assigned substation, alongside optimized flows within the transmission network. The algorithm parameters used match those of the synthetic example.

VI Conclusion

We have studied a combined space partitioning and network flow optimization problem motivated by the operation of large-scale networks serving spatially distributed consumers. By adopting an optimal transport perspective, we derived a dual formulation and designed a distributed algorithm relying only on local information exchanges. Numerical experiments confirmed the effectiveness of the method and illustrated its relevance to power distribution network reconfiguration. Future work will investigate extensions of the framework, in particular incorporating flexible generation through additional decision variables at network nodes.

-A Proof Sketch of Theorem 1

Theorem 1 can be proved by following a general methodology for OT problems, as in [4, Theorem 1.3]. Recall that for any function $f:E\to\mathbb{R}\cup\{+\infty\}$ on a normed vector space $E$ with dual $E^{*}$, its conjugate $f^{*}:E^{*}\to\mathbb{R}\cup\{+\infty\}$ is defined by $f^{*}(y)=\sup_{x\in E}\left\{\langle x,y\rangle-f(x)\right\}.$

The methodology relies on the following version of the Fenchel–Rockafellar duality theorem [4, Theorem 1.9].

Theorem 2** (Fenchel–Rockafellar).**

Let $E$ be a normed vector space and $E^{*}$ its dual. If $\Theta,\Xi:E\to\mathbb{R}\cup\{+\infty\}$ are convex functions and if there exists $z_{0}\in E$ such that $\Theta(z_{0})$ and $\Xi(z_{0})$ are finite and $\Theta$ is continuous at $z_{0}$, then

[TABLE]

To prove Theorem 1, we define suitable convex functions $\Theta$ and $\Xi$ on a suitable space $E$ and apply Theorem 2. In our case, we take $E=\mathcal{C}_{b}(X\times\mathcal{S})\times\mathbb{R}^{|\mathcal{A}|}$, where $\mathcal{C}_{b}(X\times\mathcal{S})$ denotes the space of continuous bounded functions on $X\times\mathcal{S}$ equipped with the supremum norm. By Riesz’s theorem, its topological dual space can be identified with the space of finite signed Radon measures on $X\times\mathcal{S}$, denoted by $\mathcal{M}(X\times\mathcal{S})$. Thus, the dual space is $E^{*}=\mathcal{M}(X\times\mathcal{S})\times\mathbb{R}^{|\mathcal{A}|}$, with the duality pairing between $E$ and $E^{*}$ given by

[TABLE]

We then define $\Theta(u,v)=\Theta_{\text{OT}}(u)+\Theta_{\text{MCF}}(v)$, where

[TABLE]

as in [4], and $\Theta_{\text{MCF}}(v)=\sum_{(i,j)\in\mathcal{A}}\tilde{c}_{ij}^{*}(-v_{ij})$, with $\tilde{c}_{ij}(p)=c_{ij}(p)+\delta_{[a_{ij},b_{ij}]}(p).$

As for the function $\Xi:E\to\mathbb{R}\cup\{+\infty\}$, extending the idea in [4, Theorem 1.3], if there exist functions $\varphi\in C_{b}(X)$ and $\psi\in\mathbb{R}^{|\mathcal{N}|}$ such that $u(x,i)=\varphi(x)+\psi_{i}$ for all $(x,i)\in X\times\mathcal{S}$ $v_{ij}=\psi_{j}-\psi_{i}$ for all $(i,j)\in\mathcal{A}$, then

[TABLE]

otherwise, we set $\Xi(u,v)=+\infty$. The function $\Xi$ is well-defined [4, p. 27]. With some sign changes in the variable definitions, we have that

[TABLE]

Next, after some calculations, leveraging $\tilde{c}_{ij}^{**}=\tilde{c}_{ij}$ under Assumption 1-(c) (convexity), we can show that

[TABLE]

provided that $\pi\in\mathcal{M}(X\times\mathcal{S})$, $\pi\geq 0$ and $p_{ij}\in[a_{ij},b_{ij}]$, i.e., matching (2a) under (2d), and $+\infty$ otherwise. For $\Xi^{*}$, by reparameterizing any $(u,v)$ with finite $\Xi(u,v)$ in terms of $\varphi$ and $\psi$, one obtains that $\Xi^{*}(\pi,p)$ enforces the marginal condition $\pi(A\times\mathcal{S})=\mu(A)$ for all $A\subset\mathcal{B}(X)$ (see (2c)) and the flow balance constraints (see (2b)); that is, $\Xi^{*}(\pi,p)=0$ if these conditions hold and $+\infty$ otherwise.

It follows that $\sup_{(\pi,p)\in E^{*}}\{-\Theta^{*}(-(\pi,p))-\Xi^{*}(\pi,p)\}$ is exactly the negative of the primal problem (2). Thus, by applying Theorem 2 we obtain the desired duality result. Assumption 1, with the key hypotheses that $c$ and $c_{ij}$ are lower semicontinuous and convex, $X$ is compact, and the feasibility condition $\mu(X)=\sum_{i\in N}g_{i}$, ensures that $\Theta$ and $\Xi$ satisfy the convexity, continuity and qualification conditions required to apply the theorem.

-B Proof Sketch of Proposition 1

Let $(\pi^{*},p^{*})$ be an optimal solution for (2) and $\psi^{*}$ be an optimal dual solution for (3). Since $(\pi^{*},p^{*})$ is feasible, the flow constraints hold, so by multiplying each constraint in the primal problem by $\psi^{*}_{i}$ and using strong duality, we have

[TABLE]

Substituting $q(\psi^{*})$ by its definition (4) and bringing all the terms to the right-hand side, the $\psi_{i}^{*}s_{i}$ terms vanish and we obtain $\sum_{(i,j)\in\mathcal{A}}B_{ij}+\int_{X\times\mathcal{S}}A(x,i)\,\mathrm{d}\pi^{*}(x,i)=0$, with

[TABLE]

and, $A(x,i)=[c(x,i)-\psi^{*}_{i}]-\min_{k\in\mathcal{S}}\{c(x,k)-\psi^{*}_{k}\},$ using the marginal constraint (2c). Now, because $A(x,i)\geq 0$ and $B_{ij}\geq 0$, each term must be zero.

For the flow terms, $B_{ij}(p^{*}_{ij})=0$ for every arc $(i,j)\in\mathcal{A}$ implies $p^{*}_{ij}=p_{ij}(\psi^{*})$, i.e., the minimizer in (5) for $\psi^{*}$. For the integral term, we must have $\pi^{*}$-almost everywhere that $A(x,i)=0$, i.e., $\pi^{*}$ is concentrated on the set

[TABLE]

Hence, for $\mu$-almost all $x$, when the minimizing index $i$ is unique $\pi^{*}$ must assign $i$ to $x$, which corresponds to the map $T^{*}$ in Proposition 1. When the minimizer is not unique, $\pi^{*}$ may randomize between them, however under Assumption 2, choosing the deterministic assignment $T^{*}$ also for such $x$ has no impact on the cost. Overall, this shows that $(\operatorname{id}_{X},T^{*})_{\#}\mu$ and $p^{*}$ is optimal for (2), so $(T^{*},p^{*})$ is optimal for (1) and the relaxation is tight.