Cycle Flow Formulation of Optimal Network Flow Problems and Respective Distributed Solutions

I. Introduction

IN a network flow problem, a physical system consisting of several routes between source and sink points transfers input flows from the source points to the sink points. The primary objective of optimal network flow problems is to minimize the overall cost of transporting flow [1]. Network flow problems appear in many important applications, such as software-defined networking [2], wireless sensor networks [3], transportation systems [4], [5], and power networks [6]–[8]. In power network problems, variants of optimal network flow problems also aim to include optimal generation and storage costs [9]–[12].

With the advent of new technologies, the amount of available data and network size has been increasing, which necessitate various performance improvement techniques in cyber-physical systems [13], [14], and increase the size of optimization problems. However, the number of decision variables is directly related with the time and space (resources^①) computation complexity of optimization solvers. Decision variable reduction techniques, as well as parallel/distributed optimization solutions, are therefore investigated to manage the complexity of large-scale optimization problems. Distributed solutions are of particular interest in decentralized cyber-physical operations, which aim to solve network flow problems in a scalable, fast, and decentralized manner. Although the computational cost of distributed algorithms in an optimal network flow problem is distributed among the cyber-layer nodes, the high number of decision variables normally translates to a large number of cyber nodes or large in-network communication overhead. Additionally limitations of these solutions include network congestion and the direct effect of communication cost on life-expectancy of battery-operated cyber nodes. To address these limitations, our objective is to construct a variable reduction technique that reduces the size of the decision variables in optimal network flow problems, while maintaining a sparse structure amenable to distributed solutions at a lower communication cost.

Variable fixing techniques [15], dominance techniques [16], and constraint pairing techniques [17] are general variable reduction techniques in integer quadratic problems. Also, in multi-objective optimization problems data mining techniques are used to reduce less effective variables [18]. For evolutionary optimization problems, [19] presents how variable reduction techniques can be applied to obtain the variable relations from the partial derivatives of an optimization function. For optimization problems of the form (1), eliminating affine equality constraint, as discussed below, is also a method for reducing the number of the search variables of the problem (see [20])

where $ \phi: {\mathbb{R}}^{m}\rightarrow {\mathbb{R}} $ and $ g:{\mathbb{R}}^m\to{\mathbb{R}}^r $ are the cost function and the inequality constraint function, respectively, and $ {\bf{A}} \in {\mathbb{R}}^{n \times m} $ satisfies $ {\rm rank}({\bf{A}}) = \rho\leq n < m $.

The affine feasible set for this optimization problem can be characterized as

where $ {\bf{F}} \in {\mathbb{R}}^{m \times (m-\rho)} $ is a matrix whose columns span the null-space of $ {\bf{A}} $ and $ {\bf{x}}^{\rm{p}} \in {\mathbb{R}}^{m} $ is a particular solution of $ {\bf{A}}{\bf{x}} = {\bf{b}} $. Then, $ {\bf{x}}^\star $ in (1) satisfies $ {\bf{x}}^{\star} = {\bf{F}}{\bf{z}}^\star+{\bf{x}}^{\rm{p}} $ where

Compared to (1), in (3), the equality constraint is eliminated and the number of the decision variables are reduced from $ m $ to $ m\!-\!\rho $.

Optimal network flow problems are normally in the form of the optimization problem (1), where the cost is the sum of the convex cost of the flow through the arcs subject to capacity bounds for each arc and flow conservation equations at each node. In variations of the optimal network flow problem, the cost can be augmented to include the cost of, e.g., generation and storage at nodes, and the constraints can be expanded to include other components. Nevertheless, an affine equality constraint that always is present in network flow problems is the flow conservation equation. The coefficient matrix of the linear equation describing the flow conservation constraints has a nullity of $ m-n+1 $. Therefore, one can employ the affine equality elimination method (2) to reduce the decision variables of the optimal flow problems by $ \rho = m-n+1 $ counts. This variable reduction can be substantial, since in network flow problems, the number of arcs $ m $ is usually much larger than the number of the nodes $ n $. In time-varying problems where a network flow problem is solved over a time horizon of length $ N $, the variable reduction is of $ N(m-n+1) $ counts. However, the lack of efficient methods to construct matrix $ {\bf{F}} $ and particular solution $ {\bf{x}}^{\rm{p}} $ in (2) can be an impediment when using the affine equality constraint elimination method.

Various matrix factorization techniques, including LU, QR, LQ, SVD, and Gauss-Jordan elimination are proposed to compute the span of the null-space of a matrix $ {\bf{A}}\in {{\mathbb{R}}}^{n\times m} $ [21]–[23]. However, such techniques result in either a full matrix $ {\bf{F}} $, which destroys the sparsity of the constraints and separability of the cost function of the optimal network flow problems, or a sparse $ {\bf{F}} $, whose sparsity is not aligned with the physical layer of the network, and as a result produces an equivalent optimization problem that is hard to decentralize. Distributed solutions for network flow problems, e.g., in [1], [13], [24], [25] all take advantage of the sparsity of the flow conservation equation and separability of the cost function equations of the optimal network flow formulations in their developments. Moreover, arc-based DC and AC power flow with ADMM is studied in [26], [27].

Statement of Contribution: In this paper, we use a graph theoretic approach to reduce the decision variables of optimal network flow problems by eliminating the flow conservation constraint. In particular, we show that all solutions of the flow conservation equation are characterized explicitly in terms of the span of the columns of the transpose of a cycle basis matrix of the oriented network plus a particular solution. Cycle basis of a graph can be computed in polynomial time using efficient algorithms such as those in [28], [29]. To compute a particular solution, we then propose a graph theoretic approach. We show that for any given input/output flow vector, a particular solution can be efficiently constructed from a set of elementary solutions, each obtained from tracing a flow of value $ 1 $ over the network from each node to a single common sink node. An advantage of our proposed variable reduction method is that, for fixed physical layer networks, we only need to compute the cycle basis matrix and the fundamental particular solutions once. That is, the constituting components of our method are independent of the numerical values of the cost, the arc capacities or the input/output flow values. We demonstrate the application of our flow conservation equation elimination over a minimum cost flow problem (i.e., a static optimization problem), as well as an optimal power flow problem with storage and generation at the nodes (i.e., a time-varying optimization problem). Our next contribution is to show that our proposed cycle basis reduced-decision-variable formulation is amenable to distributed solutions–this is due to the sparse and particular composition of the cycle basis matrices that result in a tangible and topologically traceable cycle flow definition for the variables in the reduced space. In this regard, we demonstrate an implantation of a distributed ADMM solution method (see [30] and [31]) for our reduced-variable reformulation of the minimum cost network flow problem. To implement this distributed solution, we propose a cyber-layer whose nodes are defined based on the cycles of the cycle basis of the physical-layer graph, used in the variable reduction stage. This architecture is of particular interest for physical networks with articulation points^②. For such networks, we show that the minimum cost optimal network flow problem decomposes into two or more decoupled smaller optimization problems, and each can be solved in a distributed manner independent from the others. Our final contribution is to show that distributed solutions for reduced-decision-variable reformulations can also be developed for conventional node-based cyber-layers. We use an example to show that the communication cost of the distributed ADMM algorithm for our proposed cycle-based formulation of the minimum cost flow problem is lower than that of a distributed ADMM algorithm for the original arc-based formulation. A preliminary version of parts of our results in this paper has appeared in [33].

Notations: $ {{\mathbb{R}}} $, $ {{\mathbb{R}}}_{>0} $, $ {{\mathbb{R}}}_{\ge 0} $, and $ {{\mathbb{R}}}_{\leq 0} $ denote the set of real, positive real, non-negative real, and non-positive real numbers, respectively. $ {\bf{1}}_n $ is the vector of $ n $ ones. $ {\bf{A}}^\top $ is the transpose of a matrix $ {\bf{A}} $ and $ [{\bf{A}}]_i $ indicate its $ ith$ column. $ [\{z_k\}_{k = 1}^n] $ is the column vector obtained from stacking the elements of an ordered set $ \{z_k\}_{k = 1}^n $. For network variables $ \{p_i\}_{i = 1}^N\!\subset\! {{\mathbb{R}}} $, defined over $ N $ nodes, $ N $ arcs, or $ N $ cycles, the aggregate vector of these variables is $ {\bf{p}}\! = \![\{p_i\}_{i = 1}^N]\!\in\! {{\mathbb{R}}}^N $. The relative complement of a set $ {\cal{B}} $ in a set $ {\cal{A}} $ is $ {\cal{A}}\backslash{\cal{B}}\! = \!\{x\!\in\!{\cal{A}}\,|\,x\not\in{\cal{B}}\} $.

Organization of the Paper: Section II defines graph-related terminologies and some of the basic properties of the cycle basis method. Section III formally presents our decision variable reduction method and its application in a minimum cost flow and an optimal power flow problems. Section IV presents a distributed ADMM solution for the cycle-flow-formulated minimum cost network flow problem and discusses the communication reduction on some example networks. Section V demonstrates the performance of the distributed cycle basis ADMM on a numerical example in terms of convergence to optimum solution and the number of communications. Section VI provides the conclusions.

II. A Review of Cycle Basis in Graphs

In this section, following [34], we review our graph-related terminology and conventions, and introduce our graph-related notations. We represent a graph of $ n $ nodes and $ m $ arcs with $ {\cal{G}} = ({\cal{V}},{\cal{E}}) $, where $ {\cal{V}} = \{ v_1,v_2, \dots , v_n\} $ is the node set and $ {\cal{E}} = \{e_1,\dots,e_m\}\in {\cal{V}}\times {\cal{V}} $ is the arc set. The graph is assumed to be undirected and with no self-loop. A walk is an alternating sequence of nodes and connecting arcs. A path is a walk that does not include any node twice, except for its first and last nodes which can be the same. A graph is connected if there is a path from its every node to every other node. The degree of a node in a graph is the total number of arcs connected to that node.

When there is an orientation assigned to the arcs of a graph $ {\cal{G}} = ( {\cal{V}},{\cal{E}}) $, we represent the oriented graph by $ {\cal{G}}^{\rm o} = ({\cal{V}},{\cal{E}}^{\rm{o}}) $. We write $ e_k = (v_i,v_j)\in{\cal{E}}^{\rm{o}} $ if arc $ e_k $ points from node $ v_i $ towards node $ v_j $. If $ (v_i,v_j)\in{\cal{E}}^{\rm{o}} $ then $ (v_j,v_i)\notin{\cal{E}}^{\rm{o}} $, i.e., there is no symmetric arc in the oriented graph. For $ {\cal{G}}^{\rm{o}} $, the oriented incidence matrix is the matrix $ {\cal{I}}^{\rm o} \in {\mathbb{R}}^{|{\cal{V}}| \times |{\cal{E}}|} $, where $ {\cal{I}}^{\rm o}_{ij} = 1 $ if arc $ e_j $ leaves node $ v_i $, $ {\cal{I}}^{\rm o}_{ij} = -1 $ if arc $ e_j $ enters node $ v_i $, otherwise $ {\cal{I}}^{\rm o}_{ij} = 0 $. For a connected graph of $ n $ nodes with a given orientation, the rank of $ {\cal{I}}^{\rm o} $ is $ n-1 $.

A cycle of $ {\cal{G}} $ is any sub-graph in which each node has even degree. A simple cycle is a path that begins and ends on the same node with no other repetitions of nodes. A cycle in the oriented graph $ {\cal{G}}^{\rm{o}} $ is a cycle in the underlying undirected graph $ {\cal{G}} $. We assign the counter clockwise direction as positive cycles orientation. Associated with each cycle is a oriented cycle vector $ {\bf{c}}^{\rm o}\in {{\mathbb{R}}}^m $ with $ {\bf{c}}^{\rm o}_{i} = 1 $ if $ e_i $ is in the cycle and aligned with its direction, $ {\bf{c}}^{\rm o}_{i} = -1 $ if $ e_i $ is in the cycle but opposing the direction of the cycle and finally $ {\bf{c}}^{\rm o}_{i} = 0 $ if $ e_i $ is not in the cycle. The following relationship exists between every oriented cycle vector and the oriented incidence matrix.

Lemma 1 (Relationship between the oriented incidence matrix and an oriented cycle vector (see [34] for proof)): In an oriented graph $ {\cal{G}}^ o $, every oriented cycle vector $ {\bf{c}}^o $ is orthogonal to every row of oriented incident matrix $ {\cal{I}}^{ o } $, i.e., $ {\cal{I}}^{ o }\,{\bf{c}}^o = {\bf{0}}_n $.■

Next, note that the vector space over Q, the Galois field of 2, generated by the oriented cycle vectors is the cycle space of $ {\cal{G}}^{\rm{o}} $. A set of cycles is called a cycle basis if it forms a basis for this vector space. When $ {\cal{G}} $ is connected, the cycle space has dimension $ \mu = m-n+1 $. Given a cycle basis, we define an oriented cycle basis matrix $ {\bf{B}}^{\rm o}\in {{\mathbb{R}}}^{\mu \times m} $ as a matrix whose rows are each the transpose of the oriented cycle vector of the cycles of this cycle basis. This matrix satisfies $ {\rm rank}({\bf{B}}^{\rm o}) = m-n+1 $. We associate a weight $ 1 $ to each arc of $ {\cal{G}}^{\rm o} $. The weight of a cycle basis is the sum of the weights of its cycles. A minimum cycle basis of $ {\cal{G}}^{\rm o} $ is a cycle basis of minimum weight.

III. Decision Variable Reduction in Optimal Network Flow Problems

In this section, we show how two well-known network flow problems can benefit from affine equality elimination method to reduce their search variables. We study our optimal network flow problems of interest over a network of $ n $ nodes, where each node is connected to a subset of other nodes with an arc. For example, in a power network the arc is a transmission line, while in a transportation network the arc is the road connecting two conjunction nodes on the road map. The physical layer topology is described by a connected graph $ {\cal{G}}_{\rm{physic}} = ( {\cal{V}}_{\rm{physic}}, {\cal{E}}_{\rm{physic}}) $, where $ | {\cal{V}}_{\rm{physic}}| = n $ and $ | {\cal{E}}_{\rm{physic}}| = m $. The flow can travel in both directions in every arc, however, we assume a pre-specified positive orientation for each arc and based on it we describe the flow network in the physical layer by the oriented version of $ {\cal{G}}_{\rm{physic}} $, i.e., $ {\cal{G}}^{\rm o}_{\rm{physic}} = ( {\cal{V}}_{\rm{physic}}, {\cal{E}}^{\rm o}_{\rm{physic}}) $. This physical network transfer flow(s) from a set of source nodes to a set of sink nodes (see physical layer in Fig. 1), while respecting the conservation of the flow constraints, i.e., the total inflow into each node must be equal to the total outflow from that node. Given $ x_i $ be the flow across the arc $ e_i\in {\cal{E}}^{\rm o}_{\rm{physic}} = \{e_1,e_2,\dots,e_m\} $, every arc $ e_i\in {\cal{E}}^{\rm o}_{\rm{physic}} $ has a pre-specified capacity, i.e., $ {\rm{b}}_{i}\leq x_{i}\leq {\rm{c}}_{i} $, for some known $ {\rm{b}}_{i},{\rm{c}}_{i}\in {{\mathbb{R}}} $. For any external flow $ {\rm{f}}_i $, we use the sign convention of $ {\rm{f}}_i>0 $ for input flow and $ {\rm{f}}_i<0 $ for output flow, and $ {\rm{f}}_i = 0 $ otherwise. Using the oriented incidence matrix, the conservation of the flow constraint at each node can be represented as

Figure  1.  A physical network with two cycle-based cyber-layer network overlaid atop. The physical layer network has two source nodes $v_1$ and $v_4$ and one sink node $v_{13}$. The physical layer has $19$ arcs and $7$ cycles. The lower cyber-layer belongs to a case that cyber-layer nodes are solving minimum network flow problem (8) using our proposed cycle-flow based distributed ADMM algorithm. The top cyber-layer belongs to a case that cyber-layer nodes are solving problem (7) using an arc-flow based distributed ADMM algorithm proposed in [24]. The connection lines in each cyber-layer are based on the variable dependency between the cyber nodes.

DownLoad: Full-Size Img PowerPoint

which in an aggregated form is represented with

where $ {\bf{f}} = ({\rm{f}}_1,\cdots,{\rm{f}}_n)^\top $. Note that since $ {\rm rank}({\cal{I}}^{\rm{o}}) = n-1 $ and $ {\bf{1}}_n^\top{\cal{I}}^{\rm{o}} = {\bf{0}} $, every set of valid exogenous input/output flow should satisfy $ \sum_{i = 1}^n {\rm{f}}_i = 0 $. The following result uses Lemma 1, to characterize the solutions of (5).

Theorem 1 (Solution set of (5)): Let $ \sum_{i = 1}^n {\rm{f}}_i = 0 $. Then, the solution set of the aggregated conservation of flow (5) is given by

where $ {\bf{x}}^{ {\rm{p}} } $ is a particular solution of (5).

Proof: The proof relies on showing that the null-space of $ {\cal{I}}^{\rm{o}} $ is spanned by columns of $ {\bf{B}}^{\rm o}\,^\top $. By virtue of Lemma 1, we have $ {\cal{I}}^{\rm{o}}{\bf{B}}^{\rm o}\,^\top = {\bf{0}} $. Recall that for a connected oriented graph $ {\rm rank}({\cal{I}}^{\rm{o}}) = n-1 $. Because $ {\rm rank}({\bf{B}}^{\rm o}\,^\top) = m-n+1 $, null-space of $ {\cal{I}}^{\rm{o}}\in {{\mathbb{R}}}^{n\times m} $ is spanned by columns of $ {\bf{B}}^{\rm o}\,^\top $.■

In what follows, we use the result of Theorem 1 to develop an affine equality elimination method to reduce the decision variables of two optimal network flow problems. However, the effectiveness of this approach depends on how efficiently one can construct matrix $ {\bf{B}}^{\rm o} $ and particular solution $ {\bf{x}}^{\rm{p}} $, especially in large scale networks. In regards to $ {\bf{B}}^{\rm o} $, there are several efficient polynomial time algorithms in the literature (see e.g., [28], [29]); see Section V for the cycle matrix corresponding to the minimum weight cycle basis of the network of Fig. 1.

Next, we propose a simple method to construct a particular solution $ {\bf{x}}^{\rm{p}} $ using graph topology. Our method relies on the superposition property of linear algebra equations, and the fact that a particular solution for a unit flow $ {\rm{f}}_i = 1 $ entering the network at node $ v_i $ and leaving it at node $ v_j $ can simply be constructed by assuming that $ {\rm{f}}_i = 1 $ flows along a path from node $ v_i $ to node $ v_j $.

Lemma 2: (Particular solution of (5)). Given an input/output flow vector $ {\bf{f}} $ over $ {\cal{G}}_{ physic }^{ o } $ which satisfies $ \sum_{i = 1}^n {\rm{f}}_i = 0 $, a particular solution for (5) is $ {\bf{x}}^{ p } = \sum_{i = 1}^{n-1}{\rm{f}}_i\, {\bar{\bf{x}}}^{ p,v_i} $. Here, $ {\bar{\bf{x}}}^{ p ,v_i}\in {{\mathbb{R}}}^m $, $ i\in\{1,\dots,{n-1}\} $, is constructed from a path that connects node $ v_i $ to node $ v_n $ such that $ {\rm{x}}^{ p ,v_i}_j = 1 $ (resp. $ {\rm{x}}_{j}^{ p ,v_i} = -1 $) if $ e_j $, $ j\in\{1,\dots,{m}\} $, is on this path and is along (resp. opposing) the direction of the path, otherwise $ {\rm{x}}_{j}^{p ,v_i} = 0 $.

Proof: For every $ i\in\{1,\dots,{n-1}\} $, consider a virtual scenario where a unit flow $ \bar{{\rm{f}}}_i = 1 $ enters the network at node $ v_i $ and leaves it at node $ v_n $. Using a simple flow tracing over the network we can see that $ {\bar{\bf{x}}}^{{\rm{p}},v_i}\in {{\mathbb{R}}}^m $ as described in the statement satisfies $ {\cal{I}}^{\rm{o}}\,{\bar{\bf{x}}}^{{\rm{p}},v_i} = {\bar{\bf{f}}}^{v_i} $, $ i\in\{1,\dots,{n-1}\} $ where $ \bar{{\rm{f}}}^{v_i}_i = 1 $, $ \bar{{\rm{f}}}^{v_i}_n = -1 $ and $ \bar{{\rm{f}}}^{v_i}_j = 0 $, $ j\in\{1,\dots,{n}\}\backslash\{i,n\} $. For a given network flow vector $ {\bf{f}} $ because $ {\rm{f}}_n = -\sum_{i = 1}^{n-1} {\rm{f}}_i $, we can write $ {\bf{f}} = \sum_{i = 1}^{n-1} {\rm{f}}_i \,{\bar{\bf{f}}}^{v_i} $. Therefore, $ {\cal{I}}^{\rm{o}} \sum\nolimits_{i = 1}^{n-1}{\rm{f}}_i\, {\bar{\bf{x}}}^{\rm{p}}_{v_i} = \sum\nolimits_{i = 1}^{n-1} {\rm{f}}_i \,{\bar{\bf{f}}}^{v_i} = {\bf{f}}. $■

A few remarks are in order regarding the particular solution. First, note that construction of the `fundamental' particular solution set $ \{{\bar{\bf{x}}}^{{\rm{p}},v_i}\}_{i = 1}^{n-1} $ is independent of the value of the network flow vector $ {\bf{f}} $. Second, for problems with single source and single sink nodes, we can label the source node $ v_1 $ and the sink node $ v_n $ and compute the elementary solution set only for node $ v_1 $. More particularly, if in a given network flow problem the sink and the source nodes are fixed we only need to compute the elementary particulars solution set for the collection of sink and source nodes. Finally, to obtain sparse elementary solutions, we can use a shortest path between nodes $ v_i $ and $ v_n $. A numerical example discussing the particular solution of the flow conservation equation for the network of Fig. 1 is given in Section V.

A Minimum Cost Flow Problem

We consider a minimum cost flow problem over $ {\cal{G}}^{\rm o}_{\rm{physic}} $ with a given set of input and output flows at specific source and sink points. In this problem, there is a convex cost $ \phi_{i}: {{\mathbb{R}}}\to {{\mathbb{R}}} $ associated with flow across each arc $ e_i\in{\cal{E}}^{o}_{\rm{physic}} $, and our objective is to find the network minimizer $ {\bf{x}}^\star\in {{\mathbb{R}}}^m $ in the following optimization problem

where $ {\bf{f}} = ({\rm{f}}_1,\cdots,{\rm{f}}_n)^\top $ is the given input/output flow vector which satisfies $ \sum_{i = 1}^n {\rm{f}}_i = 0 $. Here, (7b) captures the flow conservation at nodes across the network and (7c) describes the arc capacity constraints. The number of search variables in the optimization problem (7) is equal to the number of the arcs of the network, i.e., $ |{\cal{E}}^{o}_{\rm{physic}}| = m $. Our result below uses Theorem 1 to eliminate the affine equality constrains (7b) and reduce the search variables to $ m-n+1 $.

Theorem 2: (Eliminating the flow conservation constraint from (7)). Consider the optimal network flow problem (7) over a connected physical network $ {\cal{G}}_{ {\rm{physic}} }^ o $. Then, $ {\bf{x}}^\star $ in (1) satisfies $ {\bf{x}}^{\star} = {\bf{B}}^ {\rm{o}} \,^\top{\bf{z}}^\star+{\bf{x}}^{ {\rm{p}} } $ where

and $ {\bf{B}}^ {\rm{o}} $ is an oriented cycle basis matrix of $ {\cal{G}}^ o _{ physic } $ and $ {\bf{x}}^{ {\rm{p}} } $ is a particular solution of $ {\cal{I}}^{ o }\,{\bf{x}} = {\bf{f}} $.

Proof: The proof follows from invoking the same argument that is used to relate solutions of the optimization problem (1) to those of (3), and the result of Theorem 1.■

B Optimal Power Flow With Storage and Generation at Nodes

Next, we consider an optimal power flow problem over a network described by $ {\cal{G}}^{\rm o}_{\rm{physic}} $ with storage, generation and load at its nodes (see Fig. 2). The objective in this problem is to minimize the cost of power generation along with energy loss at the transmission lines over some finite time interval $ {\cal{T}} = \{1,\dots,T\} $. Mathematical modeling of this problem over various scenarios, including deterministic and stochastic generators, is considered in the literature [9]–[12].

Figure  2.  A schematic representation of a network flow problem with generators and energy storage and loads in a network

DownLoad: Full-Size Img PowerPoint

All of these models, at each time $ t\in{\cal{T}} $, include a flow conservation equation at each node. As a result for a network with $ m $ arcs, the flow conservation equation introduces $ m\,|{\cal{T}}| $ decision variables into the optimization problem.

In our study below, without loss of generality, we use the deterministic form (no renewable generation source) of the optimal network flow problem studied in [9], which states the problem as a direct current (DC) power flow problem (see (9)). Without loss of generality, let at each time $ t\in{\cal{T}} $, each node $ v_i\in {\cal{V}}_{\rm{physic}} $ have a generator that supplies a bounded $ \delta_{i}(t) $ power, a battery with a bounded storage level $ s_i(t) $ and a charge/discharge variable $ u_i(t) $, and a known demand $ {\rm{d}}_i(t)\in {{\mathbb{R}}}_{\leq0} $. If a node does not have any of the generation, storage, or load components, the respective variables should simply be removed from the formulation below. Then, given a known load profile $ \{{\bf{d}}(t)\}_{t = 1}^T $, where $ {\bf{d}}(t) = [\{{\rm{d}}_i(t)\}_{i = 1}^n $], the optimization problem of interest in [9] is

where the cost function includes generator costs $ g_j $, $ j\in\{1,\dots,{n}\} $, and transmission cost $ \phi_i $, $ i\in\{1,\dots,{m}\} $, following DC power grid model in [35]. Here, $ {\cal{N}}^{e}(i) $ is the set of the nodes that are connected to node $ v_i $ through an arc, and $ \lambda_i \in (0,1] $ is the storage energy dissipation factor. Moreover, $ (\underline{{\rm{s}}}_i,\bar{{\rm{s}}}_i)\in {{\mathbb{R}}}\times {{\mathbb{R}}} $, $ (\underline{{\rm{u}}}_i,\bar{{\rm{u}}}_i)\in {{\mathbb{R}}}\times {{\mathbb{R}}} $, $ (\underline{\delta}_i,\bar{\delta}_i)\in {{\mathbb{R}}}\times {{\mathbb{R}}} $, are known as lower bound and upper bound values on storage level, battery charge/discharge, and power generation by the generator, respectively. Finally $ {\bf{B}}_{ik}(\theta_i(t) - \theta_k(t)) $ is the DC approximation for alternating current power flow. Here, $ \theta_i(t) $ is the voltage phase angle of node (bus) $ v_i\in {\cal{V}}_{\rm{physic}} $ at time $ t $ and $ {\bf{B}} \in {{\mathbb{R}}}^{n \times n} $ is the imaginary part of the admittance matrix under DC assumption (for more details see [9]).

The following result shows that the number of search variables related to the flow in the optimization problem (9) can be reduced from $ m|{\cal{T}}| $ to $ (m-n+1)|{\cal{T}}| $ via eliminating the flow conservation constraint at each $ t\in{\cal{T}} $. An interesting observation in the result below is that in order to eliminate the flow conservation equations (9c), we need to introduce a new set of affine equality constraint (10b) which ensure balance between the external input and output flows.

Proposition 1: (Eliminating the flow conservation constraint from (9)). Consider the optimal power flow problem (9) over a physical network described by $ {\cal{G}}_{ {\rm{physic}} }^ o $ with a given set of loads $ \{{\bf{d}}(t)\}_{t = 1}^T $. Then, $ \{{\bf{B}}^ {\rm{o}} \,^\top{\bf{z}}^\star(t)+{\bf{x}}^{ {\rm{p}} }({{\delta}}^{\star}(t),{\bf{u}}^{\star}(t),$${\bf{d}}(t)),{\bf{u}}^{\star}(t),{\bf{s}}^\star(t),{{\theta}}^{\star}(t)\}_{t = 1}^T $ is a minimizer of the optimization problem (9) where

Here, $ {\bf{B}}^ {\rm{o}} $ is a cycle basis matrix of $ {\cal{G}}^ o _{ physic } $ and ${\bf{x}}^{ {\rm{p}} }({{\delta}}(t),{\bf{u}}(t),$$ {\bf{d}}(t)) = \displaystyle\sum_{i = 1}^{n-1}(\delta_i(t) + u_i(t)+{\rm{d}}_i(t))\, {\bar{\bf{x}}}^{ p,v_i} $, where $ \{{\bar{\bf{x}}}^{ p ,v_i}\}_{i = 1}^{n-1} $ is as described in Lemma 1.

Proof: The equality constraint (9c) in aggregated form is

Note that rank of $ {\cal{I}}^{\rm{o}}\in {{\mathbb{R}}}^{n\times m} $ is $ n-1 $ and $ {\bf{1}}_n^\top{\cal{I}}^{\rm{o}} = {\bf{0}} $, where $ {\bf{1}}_n $ is the vector of $ n $ ones. Left multiplying (11) by $ {\bf{1}}_n^\top $ results in (10b). Then, for a given $ {{\delta}}(t),{\bf{u}}(t),{\bf{d}}(t)\in {{\mathbb{R}}}^n $ that satisfy (10b), following the method discussed in Lemma 1, we can show that $ {\bf{x}}^{\rm{p}}({{\delta}}(t),{\bf{u}}(t),{\bf{d}}(t)) = \sum_{i = 1}^{n-1}(\delta_i(t) + u_i(t)+{\rm{d}}_i(t))\, {\bar{\bf{x}}}^{{\rm{p}},v_i} $ is a particular solution of (11), i.e., $ {\cal{I}}^{\rm{o}}\,{\bf{x}}^{\rm{p}}({{\delta}}(t),{\bf{u}}(t),{\bf{d}}(t)) = {{\delta}}(t)+{\bf{u}}(t)+$${\bf{d}}(t). $ Therefore, for any given load vector $ {\bf{d}} $ (to simplify the notation we drop argument $ t $), we have

As a result, at each $ t\in{\cal{T}} $, we can eliminate the affine equation (9c) and obtain the equivalent optimization problem (10). As it is described in the statement, the minimizers of problem (10) and (9) are equal.■

IV. A Cycle-Basis Distributed ADMM Algorithm for Minimum Cost Network Flow Problem

In this section, our objective is to demonstrate the effectiveness of the cycle basis decision variable reduction technique in decreasing the communication cost of distributed solutions of optimal network flow problems. We use the minimum network flow problem (7) as our demonstration case. For this problem, we first develop a distributed ADMM algorithm to solve (8). Then, we compare the communication cost of this algorithm to that of a distributed ADMM solution for (7).

To develop our results, we introduce first some notations related to the oriented cycles of $ {\cal{G}}^{\rm o}_{\rm{physic}} $. Suppose $ \mu = m-n+1 $. Let $ {\cal{C}} = \{{\cal{C}}_i\}_{i = 1}^{\mu} $, be the set of cycles of corresponding to $ {\bf{B}}^{\rm o} $ that is used to eliminate the flow conservation equation as explained in Section III. We represent the set of arcs of any $ {\cal{C}}_i\in{\cal{C}} $, $ i\in\{1,\dots,{\mu}\} $, by $ {\cal{E}}^{\cal{C}}_i = \{e_j\in {\cal{E}}^{\rm o},j\in\{1,\dots,{m}\}\,|\,{\bf{B}}^{\rm{o}}_{ij}\neq 0\} $. For a given cycle basis $ {\cal{C}} $, we refer to the cycles that share an arc as neighbors and represent the set of neighboring cycles of any cycle $ {\cal{C}}_i\in{\cal{C}} $, $ i\in\{1,\dots,{\mu}\} $, by $ {\cal{N}}^{{\cal{C}}}_i = \{j\in\{1,\dots,{\mu}\}\backslash\{i\}\,|\, \exists\, k\in$$\{1,\dots,{m}\}\;{\rm s.t.}\; {\bf{B}}^{\rm o}_{ik}\neq 0 \;{\rm and}\; {\bf{B}}^{\rm o}_{jk}\neq 0 \} $.

We let $ {\cal{I}}^{\rm{c}}(e_i) $ be the set of indexes of the cycles that arc $ e_i\in {\cal{E}}_{\rm{physic}} $ belongs to them, i.e., $ {\cal{I}}^{\rm{c}}(e_i) = \{j\in\{1,\dots,{\mu}\}\,|\, e_i\in{\cal{E}}^{\cal{C}}_j\} $. We let $ {\cal{J}}^{\cal{E}}(v_i) $ be the set of all the arcs that go through node $ v_i\in {\cal{V}}_{\rm{physic}} $, i.e., $ {\cal{J}}^{\cal{E}}(v_i) = \{v_j\in {\cal{V}}_{\rm{physic}}\backslash\{v_i\}\,|\,(v_i,v_j)\in{\cal{E}}_{\rm{physic}}\} $.

Recall that to eliminate the flow conservation equation we used $ {\bf{x}} \!=\! {\bf{B}}^{\rm o}\,^\top{\bf{z}}\!\!+{\bf{x}}^{\rm{p}} $, or equivalently $ x_i \!=\! {\bf{z}}^\top[{\bf{B}}^{\rm o}]_i\!+\!{\rm{x}}_i^{\rm{p}} $, $ i\in\{1,\dots,{m}\} $. Since, for a given arc $ e_i\in {\cal{E}}^{\rm o}_{\rm{physic}} $, every element of $ {\bf{B}}^{\rm o}_{ji} $ is zero except if cycle $ {\cal{C}}^{\rm of}_j $ contains arc $ e_i $, every $ x_i $, $ i\in\{1,\dots,{m}\} $, is an affine function of $ {\rm{x}}_i^{\rm{p}} $ and $ \{z_k\}_{k\in{\cal{I}}^{\rm{c}}(e_i)} $. If we think of every $ z_i $, $ i\in\{1,\dots,{\mu}\} $ as a cycle flow variable (with positive direction in counterclockwise direction) of the cycle $ {\cal{C}}_i $, then every arc flow $ x_j $, $ j\in\{1,\dots,{m}\} $ is a function of its corresponding component of the particular solution and the cycle flows of the cycles that contain the arc. Given such relationship, the cost function of every arc is

Next, we derive equivalent representations of optimization problem (8) that can be solved in a distributed manner using the ADMM algorithm of [13] via two different cyber-layer architectures.

Cycle-Based Cyber-Layer: to develop our first distributed solution to solve (8), we assign a cyber-layer node to each cycle of the cycle basis that we used in our decision variable reduction stage (see Fig. 1 as an example). We refer to this architecture as cycle-based cyber-layer. We assume that the cyber-layer nodes of the neighboring cycles can communicate with each other in a bi-directional way. This procedure results in a cyber-layer with $ \mu $ nodes. To obtain cycle basis with fewest number of arcs in each cycle we can use a minimum weight cycle basis algorithms (see e.g. [29]).

Now, for every cyber-layer node $ i\in\{1,\dots,{\mu}\} $, we define $ {\bf{y}}_i = (\bar{y}_i,\tilde{{\bf{y}}}_i)\in {{\mathbb{R}}}^{|{\cal{N}}^{{\cal{C}}}_i|+1} $, where $ \bar{y}_i\in {{\mathbb{R}}} $ is the local copy of $ z_i $ and $ \tilde{{\bf{y}}}_i $ is the local copy of $ \{z_k\}_{k\in{\cal{N}}^{{\cal{C}}}_i} $ at cyber node $ i $. With this definition, we assume that every cyber node besides its own corresponding cycle flow has also a copy of cycle flow variable of each of its neighbors. Next, we cast the cost function of each cyber node in terms of its decision variable $ {\bf{y}}_i $. Let $ {\bf{y}}_i(e_k) $ be the component(s) of $ {\bf{y}}_i $ corresponding to $ \{z_j\}_{j\in{\cal{I}}^{\rm{c}}(e_k)} $. For every cyber-layer node, we define its cost function as

Then, we can cast the minimum cost network flow problem (8) in the following equivalent form

Set of constraints at each cyber agent $ i\in\{1,\dots,{\mu}\}:$

In this formulation, every cycle-based cyber node has a copy of the cycle flows that go through its arcs, i.e, $ {\bf{y}}_i $. The equality constraint at each node ensures that local copies of the cycle flows of the neighboring agents are the same, while the inequality constraint ensures that the flow through the arcs of each cycle respect the capacity bounds. The new formulation (13) fits the standard framework developed for distributed ADMM solutions and can be solved for example using the algorithm of [13]. The details are omitted for brevity. Once every cyber node $ i\!\in\{1,\dots,{\mu}\} $ computes $ {\bf{y}}_i^\star $, then it can use $ x_k^\star = {\bf{y}}^\star_i(e_k)^\top[\{{\bf{B}}^{\rm o}_{jk}\}_{j\in{\cal{I}}^{\rm{c}}(e_k)}]+{\rm{x}}^{\rm p}_k $ to obtain the optimal cycle flows through its arcs $ e_k\in{\cal{E}}^{\cal{C}}_i $. In this distributed implementation, we assume that elementary particular solution set $ \{{\bar{\bf{x}}}^{{\rm{p}},v_i}\}_{i = 1}^{n-1} $ are computed off-line and are available at cyber nodes. At operation times, we only need to broadcast the input/output flow vector $ {\bf{f}} $ to the cyber-layer agents. For large cycles, one can split the cycle among several cyber nodes. A numerical example demonstrating our proposed distributed ADMM algorithm over a cycle-based cyber-layer can be found in authors' preliminary work [33].

One can solve also the original minimum cost network flow problem (7) in a distributed manner by a cycle-based cyber-layer using, for example, the distributed ADMM algorithm of [24]. To do so, every cyber-layer node needs a local copy of the arc flows across all the nodes in its corresponding cycle. The local copies are required to create a local copy of the flow conservation equation of the nodes. Then, to coordinate these local copies, the cyber-layer agents that share a node are required to communicate with each other. Since there are more cycles that are connected to each other through the nodes than those connected by arcs, as shown in Fig. 1, this distributed solution requires more connections in the cyber-layer than the distributed solution for (8). Moreover, because there are more arc flows than cycle flows, the message sizes exchanged among neighboring cyber agents solving (7) is larger. Therefore, the distributed solution for (7) will be less favorable from network congestion and communication energy consumption perspective.

Remark 1: (Example networks). In a square mesh physical layer network with $ n = N^2 $ nodes, the total number of arcs is $ m = 2(N^2 - N) $, and the total number of cycles is $ c = N^2 - 2N + 1 $, see the top plots in Fig. 3. In the mesh networks of Fig. 3, cyber-layer constructed based on the cycle-basis partitioning of the physical layer is shown via the thick blue network. For a cyber-layer which solves the alternative optimal network flow problem (8) via an ADMM distributed solution, a cyber agent only needs to communicate with neighboring agents which have a common arc. The number of communication links in cyber-layer is $ m_1 = $ $2((N-1)^2 - (N-1)) =2(N^2 - 3N + 2) $. For a cyber-layer which solves the original optimal network flow problem (7) via an ADMM distributed solution, any cyber agent needs to communicate with neighboring agents which has a common node in the graph. Therefore, the number of communication links in cyber-layer is $ m_2 \!=\! ((N-1)^2 \!-\! 2(N-1) + 1)2 \!+\! 2(N^2 - $$3N + 2) =4N^2 - 14N + 12 $. Consequently, we can obtain that $ \lim_{n \rightarrow \infty} \dfrac{m_2}{m_1} = 2 $. An interesting case involving our proposed cycle-based distributed ADMM solution of (13) is when the physical layer network has articulation points (e.g., node $ 27 $ in IEEE bus system 30 in Fig. 4 and node $ 9 $ of the physical layer in Fig. 1 are articulated points). The physical layer graph with an articulation point is consisted of subgraphs that are connected to each other through simple nodes with no common arcs between them. Therefore, when our proposed aforementioned ADMM algorithm is implemented, there is no variable dependency between the cyber nodes of these subgraphs. This means that the original fully coupled optimization problem is now divided into smaller optimization problems that each can be solved independently. For example in the network of Fig. 1, cyber-layer node 1 can obtain the optimal solution across its arcs in one ADMM iteration without a need to coordinate with other cyber agents. An extreme example case is also shown in the bottom plots of Fig. 3. In this example there is no common arc among cycles, but they all have one common node. The number of communication links for cycle-based cyber-layer graph is zero. In this case, a cycle-based ADMM algorithm finds the optimum solution in one step, without any communication. However, in arc-based network flow, the cyber-layer is a complete graph as all of the agents are connected to each other.

Figure  4.  The graph related to IEEE bus system 30 with 41 arcs, 30 nodes and 12 cycles.

DownLoad: Full-Size Img PowerPoint

Figure  3.  Examples of physical layer networks (black network) with the corresponding cycle-based cyber-layer atop (blue thicker network).

DownLoad: Full-Size Img PowerPoint

Node-Based Cyber-Layer: We can solve problem (8) in a distributed manner with the conventional node-based cyber-layer, as well. In a node-based cyber-layer architecture, a cyber agent is assigned to each node of the physical layer to compute the flow across all the incident arcs of the corresponding physical layer node. For example, in traffic networks, a cyber agent is assigned to each intersection, which is a node in the physical layer. In this case, the topology of the cyber-layer is exactly the same as the physical layer. To solve problem (8) using a distributed ADMM algorithm over such cyber-layer, we assume that each cyber agent has a local variable $ {\bf{y}}_i $ consisted of local copies of the cycle flows that goes across its arcs. We also split the flow cost across each arc between each end nodes. Then, the total local cost at each cyber node is

Then, we can cast the minimum cost network flow problem (8) in the following equivalent form

Set of constraints at each cyber agent $ i\in\{1,\dots,{n}\}:$

Given this equivalent representation, similar to the method described for the cycle-based cyber-layer, we can now solve the optimization problem (8) using a distributed ADMM algorithm. In a cycle flow representation, it is more likely to have fewer cycles that go through a node than the arcs that are incident at that node. Therefore, in a distributed ADMM solution, the size of the broadcast messages of a cyber agent solving (8) is more likely to be less than the size of the broadcast messages when we solve (7).

V. Demonstrative Examples

We demonstrate the use of distributed cycle-flow-based ADMM algorithm for a minimum cost optimal flow over the network shown in Fig. 1, and compare it with the arc-flow-based ADMM distributed algorithm [24]. The arrows in the physical layer show the positive flow orientation assigned to the arcs. The cyber-layer is generated based in the minimum weight cycle basis partitioning of the physical layer as shown in the bolder network with blue nodes in Fig. 1. In this problem, we set capacity bounds $ {\rm{b}}_i = -\,{\rm{c}}_i $, and $ {\rm{c}}_i\in [1,20],\; i\in\{1,\dots,{19}\} $. We assume that the cost of the network flow at each arc is given as $ \phi_i(x_i) = \frac{1}{2}{\rm{a}}_ix_i^2 +{\rm{b}}_ix_i $, where $ x_{i} = {\bf{z}}^\top[{\bf{B}}^{\rm o}]_i+{\rm{x}}_i^{\rm{p}} $ is the arc flow and $ {\bf{z}} = (z_1,\cdots, z_7)^\top $ are the cycle flows. The parameters of the cost functions are chosen randomly from $ {\rm{a}}_i\in[1, 20],\; i\in\{1,\dots,{19}\} $ and $ {\rm{b}}_i\in[1, 20],\; i\in\{1,\dots,{19}\} $. In the physical layer network in Fig. 1, there are two source nodes $ v_1 $ and $ v_4 $ and one sink node $ v_{13} $. For our selected capacity bounds, using Edmonds-Krap algorithm [36] we obtain the maximum input flow $ {\rm{f}}_1+{\rm{f}}_4 $ to be $ 30 $. In our simulation, then we set $ {\rm{f}}_1 = 10 $ and $ {\rm{f}}_4 = 15 $. The output follows are $ {\rm{f}}_{13} = -25 $. For the network of Fig. 1, there are 7 cycles and 19 edges in the network and the orientated cycle basis matrix is, (recall that we have assumed the positive cycle flow direction to be counterclockwise). We follow Lemma III.1 to generate a particular solution for given input output flow vector $ {\bf{f}} = ({\rm{f}}_1,0,0,{\rm{f}}_4,0,0,0,0,0,0,0,0,{\rm{f}}_{13})^\top $, where f₁₃ = –(f₁ + f₄) (recall that input flows have positive and output flows have negative values). We compute the elementary particular solutions for nodes $ v_1 $ and $ v_4 $ using shortest path from them to node $ v_{13} $: $ {\bar{\bf{x}}}^{{\rm{p}},v_1} = (1,{\bf{0}}_{1\times 4},1,{\bf{0}}_{1\times 6},1, 0, 0, 1, 0, 1, 0)^\top $, $ {\bar{\bf{x}}}^{{\rm{p}},v_4} = ({\bf{0}}_{1\times 7}, $$1,0,1,{\bf{0}}_{1\times 4},1, 1, 0, 1, 0)^\top$. Then $ {\bar{\bf{x}}}^{\rm{p}} = {\rm{f}}_1\,{\bar{\bf{x}}}^{{\rm{p}},v_1}+{\rm{f}}_4\,{\bar{\bf{x}}}^{{\rm{p}},v_4} $.

The central optimum solution of our example is ${\bf{x}}^\star \! = \{ 9.64, $2.35, –2, 0.03, 3.70, 5.90, 2.39, 13, 3, 10, 3.14, 5.94, 9.05, 5.94, 10, 14, 11, 14, 11} with objective function value of $ {f}^\star = 9316.8 $ (also, $ {\bf{z}}^\star\! = \! \{0.35, -2, 0.39, -5, 4.1, 0.95, 11\} $). Fig. 5 shows the absolute percentage error, defined as $ \epsilon_f = ({f}^k - {f}^\star) / {f}^\star$ at iteration $ k $ of the cycle-flow-based and arc-flow-based distributed ADMM algorithms. The result shows that the convergence of the cycle-flow-based ADMM is faster than arc-flow-based ADMM when the problem is solved via the cycle-based cyber-layers shown in Fig. 1 for each of these approaches. To compare the number of communications between cycle-flow-based and arc-flow-based ADMM algorithms, we set a threshold of $ \epsilon_f = 0.01 $ for the absolute percentage error. Distributed cycle-flow-based ADMM converges after $ 5 $ itereation to $ \epsilon_f = 0.01 $ absolute percentage error while the distributed arc-flow-based ADMM converges after $ 13 $ iterations. Each iteration of the ADMM algorithm needs to send primal variables twice for every communication link $ (i,j) $ in cyber network from node $ i $ to node $ j $ and vice versa. The number of communication links in cycle-flow-based and arc-flow-based cyber-layer are $ 7 $ and $ 14 $ respectively. As a result, the number of communications for cycle-flow-based ADMM is $ 35 $ and for arc-flow-based ADMM is $ 112 $. In Fig. 6, we compare the number of communications between cycle-flow-based and arc-flow-based distributed ADMM for the cyber nodes shown in Fig. 1. The number of communications is obtained by adding up the communication incidences at each iteration, until we have $ \| {\bf{x}}_i^k-{\bf{x}}_i^\star\| < \epsilon = 0.1 $, where $ {\bf{x}}_i^\star $ is the optimum solution of the central problem. Fig. 6 shows that the cycle-flow based algorithm requires fewer communications.

Figure  6.  Comparison of the number of communications between the cyber-layer nodes when they solve the optimal network flow problem using a distributed cycle-flow-based and arc-flow-based ADMM algorithms until we have $\| {\bf{x}}_i^k-{\bf{x}}_i^\star\| < \epsilon=0.1$, $i \in \left\{ {1, \ldots ,19} \right\}$. As seen, in the arc-flow-based algorithms the cyber nodes communicate more.

DownLoad: Full-Size Img PowerPoint

Figure  5.  $ \epsilon_f = ({f}^k - {f}^\star )/ {f}^\star$ versus iteration number $k$ for the arc-flow-based and cycle-flow-based distributed ADMM algorithms over the cyber-layer architectures shown in Fig. 1.

DownLoad: Full-Size Img PowerPoint

VI. Conclusion

We have considered optimal network flow problems and investigated how the decision variables size of these problems can be reduced by eliminating the affine flow conservation equations. Our study was based on exploiting cycle basis concept from graph theory to eliminate the flow conservation equation in an efficient manner. In particular, we have shown that the components of our variable reduction method can be obtained in a systematic manner using graph theoretic approaches. Moreover, we have shown that the new formulation of the optimal network flow problems with reduced variables has a sparse structure and can be solved via distributed optimization solvers. In this regard, we have demonstrated the use of a distributed ADMM solver for the cycle-flow-based minimum cost flow formulation, and showed that this distributed operation leads to a reduced communication cost among the cyber-layer nodes.