FORWARD: Feasibility Oriented Random-Walk Inspired Algorithm
for Radial Reconfiguration in Distribution Networks

We consider a network-flow distribution problem over a bidirectional distribution network $\mathcal{G}_{D}=\mathcal{G}(\mathcal{V}_{D},\mathcal{E}_{D})$ with $|\mathcal{V}_{D}|=N$ nodes, $|\mathcal{E}_{D}|=m$ edges, a set of $n_{g}$ source nodes $\mathcal{V}_{g}\subset\mathcal{V}_{D}$, and $n_{c}=N-n_{g}$ sink nodes $\mathcal{V}_{c}=\mathcal{V}_{D}\backslash\mathcal{V}_{g}$, each with specified input and output. We let $\mathbf{d}\in_{\geq 0}^{N}$ be the output vector and $\mathbf{g}\in_{\geq 0}^{N}$ be the input vector, where $g_{i}\in_{>0}$ (resp. $d_{j}\in_{>0}$) for $i\in\mathcal{V}_{g}$ (resp. $j\in\mathcal{V}_{c}$) and $g_{i}=0$ (resp. $d_{j}=0$) otherwise. The assumption is that inputs match outputs, i.e., $\sum_{i\in\mathcal{V}_{g}}g_{i}=\sum_{i\in\mathcal{V}_{c}}d_{i}$. The goal is to find an (oriented) radial configuration that delivers the input flow from the source nodes to the sink nodes with minimal overall cost, see Fig. 1. This cost results from the ‘resistance’ or ‘toll’ along the edges, characterized as a quadratic function of the flow across them.

The problem of interest can be formalized as

where $x_{i,j}\in_{>0}$ is the flow across the link $(i,j)$, $A(\mathcal{S})$ is the incidence matrix of the radial configuration (a decision variable of the optimization problem), $C_{i,j}$ is the coefficient of the cost of edge $(i,j)\in\mathcal{S}$. The constraints confine the radial configuration to $\mathcal{F}$ and enforce the flow conservation (Kirchhoff’s law) at the nodes. Throughout this paper, we assume that there are several radial configurations for which the input can meet the specified output. Thus, the feasible solution set of the optimization problem (1) is non-empty.

Problem (1) is highly relevant to various potential flow applications, including natural gas, water, and electricity distribution networks. For instance, in power systems, the efficient and reliable distribution of energy has become increasingly critical with the growing integration of renewable energy and distributed generation sources. Modern power distribution networks, comprising multiple distributed generators, must operate in a radial configuration to adhere to engineering and safety standards. Minimizing energy loss within these networks is vital for both economic viability and environmental sustainability. Consequently, radial configurations cannot be chosen arbitrarily; they must be designed to ensure feasibility and to optimize energy loss (cost of operation).

Optimal radial reconfiguration problems are NP-hard problems, primarily due to the exponential growth in possible configurations as the number of distribution links increases. Specific instance (1) is often cast as Mixed-Integer Non-Linear Programming (MINLP) because of problem-specific side constraints, but for the purposes of this paper we focus on this linear abstraction to address the salient complexities associated with modeling (1b). These problems remain computationally intensive, often requiring extensive computational time and lacking guarantees of optimality [3].

To address these challenges, this paper constructs a polynomial-time solution for (1) by leveraging the graph-like characteristics of distribution systems and how potential flows naturally navigate through a network. Our algorithm uses a greedy radial construction process starting from source nodes, incrementally adding edges while considering the flow across constructed components and the demanded output of the remaining nodes to reach. We draw inspiration from the similarities between electric flow in power networks and random walks [4], developing a novel ‘sampling’ method²²2Although our radial configuration construction is deterministic, we use ‘sampling’ as a conceptual analogy. for constructing radial configurations. In the random walk approach to describing electricity distribution in a given network, a weight proportional to the inverse of the resistance along the corresponding link in the power network is assigned to each link. This creates a notion of the edge weights being the conductance of the edge. Just as electricity flows through paths of least resistance, a random walk probabilistically selects paths based on transition probabilities influenced by edge weights.

Our method assigns weights to edges based on cost and downstream output, ensuring flow distribution aligns with network requirements and effectively delivers a radial configuration that meets demand. This methodology ensures the process aligns with both the physical properties of flow and the optimization constraints of the problem, enabling the efficient identification of feasible and near-optimal radial configurations in polynomial time. Our particular innovation is incorporating a mechanism (Net-Concad function explained in Section II-C) in our incremental ‘sampling’ process that addresses the shortsightedness of greedy selections by informing the process of the comprehensive demand of the remaining nodes. Numerical examples demonstrate the effectiveness of our proposed algorithm.

Related work: Several methods in the literature, especially in power networks applications, leverage topological properties to address the reconfiguration problem [5]. For example, [6] [7] use spectral clustering followed by local greedy search to identify radial configurations. In [8], the reconfiguration problem is linked to the maximum flow problem, based on Ford-Fulkerson’s work [9], which has been extensively studied in both single- and multi-source contexts [10], [11], [12]. These approaches often include a repair procedure for unfeasible solutions, which can be inefficient. While some methods address minimum-cost distribution, radiality remains a critical constraint [13], and cycle-breaking methods used to tackle it offer limited guarantees for large graphs [14]. Other approaches, such as those based on the minimum spanning tree (MST) problem, face challenges due to the constraints of the reconfiguration problem, making the MST approach NP-hard in this context (see Fig. 2). Moreover, although algorithms like Kruskal’s or Prim’s [15] can find MSTs in polynomial time, they lose their optimality when multiple sources are involved, note that the minimum spanning forest (MSF) is not necessarily a subset of the MST (see Fig. 3).

Our proposed algorithm follows a greedy structure, typical of MST algorithms. Starting from an empty edge set $\mathcal{S}$, we iteratively add edges connecting the nodes, guided by a heuristic to ensure feasibility. Inspired by the similarity between flow in distribution networks and random walks [4], we use a dual graph where edge weights represent the inverse of distribution costs. Our random walk starts at source nodes, assigning weights (intuitively interpreted as probabilities) to edges and adding the one with the highest probability of traversal to the growing forest.

We introduce some preliminary steps to sampling to compensate shortsightedness in greedy decision making algorithms by reducing the set of possible edges to select at each iteration. The first step is implementing the Pre-processor function which is followed by the Islander function to simplify the search space. Then, we apply sequentially Net-Concad and Sampler functions for finite number of times to construct the radial configuration. These four functions constitute the elements of our proposed algorithm, which we call FORWARD as defined in the title of the paper. These functions are explained one by one below. Before that we introduce some notations. We define the input vector $\mathbf{p}$, associated with the nodes $\mathcal{V}_{D}$. This vector will be updated at each iteration of the algorithm, thus we denote it by $\mathbf{p}^{t}$ at each iteration $t$. We initialize at $p_{i}^{0}=-d_{i}<0$ for $i\in\mathcal{V}_{c}$ and $p_{i}^{0}=g_{i}>0$ if $i\in\mathcal{V}_{g}$. In what follows, we define the operator $\mathcal{V}(\mathcal{S})$ as the nodes of the sampled edge set $\mathcal{S}$, $\mathcal{E}(\mathcal{V}_{i})$ as the set of edges in $\mathcal{E}_{D}$ interconnecting nodes in $\mathcal{V}_{i}$ and $\mathcal{N}(\mathcal{V}_{i})$ as the set of nodes connected to nodes in $\mathcal{V}_{i}$. Also, $(i\to j)$ denotes a directed edge indicating flow from node $i$ to $j$. Given a directed edge set $\mathcal{S}$, $\mathcal{U}(\mathcal{S})$ returns the undirected edge version of $\mathcal{S}$.

To find a feasible radial distribution configuration from $\mathcal{F}$ that can channel the inputs from the sources to the sinks and meet their demanded output is a non-trivial task. In what follows, we show that the FORWARD algorithm is guaranteed to generate such a feasible radial distribution configuration because of the inherent structure of the condensed dual graph. We start our demonstration with some auxiliary lemmas and remarks.

Lemma III.1 asserts that in the sub-graphs separated from $\mathcal{G}_{D}$ by the Pre-Processor the (oriented) radial configuration to deliver the inputs to the sinks exists and is unique. Therefore, polytree grown by Pre-Processor belongs to all the configurations that are in the feasible set of optimization problem (1), including the optimal configuration.

Note that the Islander function will not affect the sampling procedure as it re-distributes input vector $\mathbf{p}$ without disturbing the flow balance. In a similar way, applying Net-Concad to a graph $\mathcal{G}^{\ell}$ will not modify the sampling space. By construction, edges on the dual condensed graph, created by Net-Concad, $\bar{\mathcal{E}}^{\ell}$ are all the existing edges in $\mathcal{E}_{P}$. Consequently, feasibility is guaranteed if we demonstrate that the FORWARD can always construct a feasible radial configuration in each $\mathcal{G}^{\ell}$ that delivers inputs in the sub-graph to the sinks (i.e., inputs meet the demands at the sinks).

The first result below shows that any digraph created by Sampler function through iterations is a polyforest (a radial configuration without a cycle).

The following lemma shows that any radial digraph created by FORWARD will visit every node in the graph.

In light of the aforementioned, we can conclude the demonstration by noting that the solution provided by FORWARD will be a feasible radial configuration for flow distribution.

We demonstrated the efficacy of the proposed FORWARD algorithm in solving an optimal power distribution problem and discuss how edge weights (used in line 4 of Sampler, Algorithm 4) can be selected for this problem considering the physical nature of the flow in power networks. It is shown in the literature, single-phase direct current (DC) simplifications can be used to approximate the real three-phase alternating current (AC) flow problem, casting the optimal power flow problem as (1). We refer the readers to [19] for specific details of these simplifications.

In a power network, the demand at the sinks (loads) and the input in the sources (generators) are power. The flow in edges is characterized by the demand of the nodes in each distribution line following Ohm’s law. In power network systems, $C_{i,j}$ in the cost function of (1) is the resistance $R_{i,j}$ along the edge $(i,j)$.

Weight design: Given a radial configuration, the flow $x_{i,j}$ at edge $(i\to j)$ is defined as a function of the power demand $d_{i\to j}$ of all the loads in power line extending from the source and containing edge $(i\to j)$ and the voltage in node $j$. Therefore, the cost due to flow through $(i\to j)$ is

FORWARD constructs the radial configuration by iteratively adding the edges. Therefore, the power line connecting the nodes is not known in advance. Consequently, to have a cost-aware edge choice, the greedy selection process must estimate $d_{i\to j}$. For power networks we propose to approximate the energy lost due to flow across an edge $(i\to j)$, when we want to sample this edge and add it to an existing radial configuration (power line) $\mathcal{T}^{\ell}_{k}$, to use

where $p_{j}$ is the power demand at node $j$ and $\widetilde{h}(i)$ is the accumulated cost of power line $\mathcal{T}^{\ell}_{k}$. Thus, a relevant edge weight (probability) for the random walk could be computed as

where $p_{i}$ is the input vector characterized as the power remaining in node $i$ which decreases in each step by subtracting the power demand on $j$. Therefore, the weights on each edge are dynamically changing in function of the steps taken. We point out here that of use of transition probability for (3) is with a degree of abuse of notation as in function of the parameters $w_{i,j}$ can be greater than $1$. Indeed, $w_{i,j}$ is normalized at each iteration considering all the edge candidates to be sampled at time $t$.

Numerical demonstration: In order to show computational benefit of using FORWARD in producing a feasible radial configurations for power network problems and the optimality gap of the algorithm, we consider $6$ different power distribution networks. Specifically, three IEEE graphs have been used: IEEE $13$ with $2$ sources, IEEE $18$ with $2$ sources and IEEE $33$ with $3$ sources. Note that, as mentioned before, these IEEE graphs have been pre-processed to a single-phase DC network structure. Besides, to test scalability, three random distribution systems have been generated following Watts–Strogatz mechanism [20] in order to create small-world graph structures. These graphs are named: WS $120$, WS $240$ and WS $400$, where their number indicates the number of nodes in each network. Networks WS $120$ and WS $240$ where designed with $10$ sources and WS $400$ has $20$ source nodes. We have modeled our numerical problem in PowerDistributionModel (PMD) framework [21] from Los Alamos National Laboratory, which is coded in Julia. This framework uses Knitro, from Artelys, to solve the MINLP. This numerical study was conducted on a MacBook Air with an M3 chip and $24$ GB of RAM. The simulation results are shown in Table I. In this table, power loss is in giga-watts (GW). In time column, $*$ stands for processes which has been manually stopped after certain CPU Time, and $-$ stands for situations where a solution could not be found after $3$h of computation.

According to Table I, Knitro was only able to return a solution for small-sized networks. For these cases, the CPU time of Knitro is significantly more than the CPU time used to implement FORWARD. Knitro (like other commercial solvers) uses heuristics to warm up the process within the feasible domain. Notably, as the network size increases, these heuristics tend to struggle to find a proper initialization within polynomial time due to the large dimension of the combinatorial space. This often results in prolonged CPU times or even the inability to find an initial point. In contrast, FORWARD was able to construct a feasible radial configuration for all cases in noticeably low time, even for a network of 400 nodes.

The results shown in Table I for the cost indicate that FORWARD attains a promising optimality gap. This could be attributed to the physically flow-aware nature of designing the weights and the sampler mechanism of FORWARD. In our future work, we plan to formally characterize the optimality gap of FORWARD.

We presented an algorithm for finding a feasible radial reconfiguration in distribution networks. Promising results from numerical experimentation allows us to believe the proposed method can serve as an effective warm-up strategy for iterative solvers refining solutions toward optimality. Also, it should be emphasized that despite the algorithm is designed for a simplified transportation problem, it can be seen it handles extended problems with complex physics constraints as they are power systems reconfiguration. For future work, it will be shown how, indeed, proposed problem works as an abstraction over other problems with non-linear constraints.