Cycle flow formulation of optimal network flow problems for centralized and decentralized solvers

TL;DR

This paper introduces a cycle basis-based formulation for large-scale network flow problems, reducing decision variables and enabling efficient centralized and decentralized solutions, demonstrated on minimum cost flow and power flow problems.

Contribution

It proposes a novel cycle basis approach to reformulate network flow problems, simplifying constraints and facilitating distributed optimization methods.

Findings

Cycle basis reduces decision variables in large network flow problems.

The formulation enables distributed solution via ADMM.

Demonstrated effectiveness on power and minimum cost flow problems.

Abstract

When the underlying physical network layer in optimal network flow problems is a large graph, the associated optimization problem has a large set of decision variables. In this paper, we discuss how the cycle basis from graph theory can be used to reduce the size of this decision variable space. The idea is to eliminate the aggregated flow conservation constraint of these problems by explicitly characterizing its solutions in terms of the span of the columns of the transpose of a fundamental cycle basis matrix of the network plus a particular solution. We show that for any given input/output flow vector, a particular solution can be efficiently constructed from tracing any path that connects a source node to a sink node. We demonstrate our results over a minimum cost flow problem as well as an optimal power flow problem with storage and generation at the nodes. We also show that the new formulation of the minimum cost flow problem based on the cycle basis variables is amenable to a distributed solution. In this regard, we apply our method over a distributed alternating direction method of multipliers (ADMM) solution and demonstrate it over a numerical example.

Figures

Equations

$${\mathbf{x}}^{\star }=\underset{\mathbf{x}{\in }^n}{\mathrm{argmin}}\text{leavevmode}\varphi (\mathbf{x}),\text{leavevmode}\text{leavevmode}\text{s.t.}\mathbf{A}\mathbf{x}=\mathbf{b},\text{leavevmode}\text{leavevmode}\mathbf{g}(\mathbf{x})\le 0,$$

$$\{\mathbf{x}{\in }^n\mid \mathbf{A}\mathbf{x}=\mathbf{b}\}=\{\mathbf{F}\mathbf{z}+{\mathsf{x}}^{\text{p}}\mid \mathbf{z}\in {\mathbb{R}}^{n-\rho }\}.$$

$${\mathbf{z}}^{\star }=$$

$$\bar{\mathbf{g}}(\mathbf{z})=\mathbf{g}(\mathbf{F}\mathbf{z}+{\mathsf{x}}^{\text{p}})\le 0.$$

$${\mathbf{x}}^{\star }$$

$${\sum }_{j=1}^m{\mathfrak{I}}_{ij}^{\text{o}}{x}_j={\mathsf{f}}_i,i\in \{1,\dots ,n\},$$

$${\mathsf{b}}_j\le x_j\le {\mathsf{c}}_j,\text{leavevmode}j\in \{1,\cdots ,m\},$$

$${\mathbf{z}}^{\star }=$$

$${\mathsf{b}}_j\le {\mathbf{z}}^{\top }[{\mathsf{B}}^{\text{of}}{]}_j+{\mathsf{x}}_j^{\text{p}}\le {\mathsf{c}}_j,\text{leavevmode}j\in \{1,\cdots ,m\},$$

$${\mathfrak{I}}^{\text{o}}\mathbf{x}=\mathsf{f}.$$

$$\{{\mathbf{x}}^{\star }(t),{\delta }^{\star }(t),{\mathbf{u}}^{\star }(t),{\mathbf{s}}^{\star }(t),{\theta }^{\star }(t){\}}_{t=1}^T=$$

$$\displaystyle\operatorname{argmin}\leavevmode\nobreak\ \frac{1}{T}\sum\nolimits_{t=1}^{T}\!\Big{(}\sum\nolimits_{j=1}^{n}\!g_{j}(\delta_{j}(t))+\!\sum\nolimits_{i=1}^{m}\!\phi_{i}(x_{i}(t))\Big{)},\leavevmode\nobreak\ \text{s.t.}$$

$$\text{forleavevmodenobreak }t\in \mathcal{T},\text{leavevmode}\text{leavevmode}i\in \{1,\dots ,n\}$$

$$\text{leavevmode}\text{leavevmode}{\sum }_{j=1}^m{\mathfrak{I}}_{ij}^{\text{o}}{x}_i(t)={\delta }_i(t)+u_i(t)+{\mathsf{d}}_i(t),$$

$$\text{leavevmode}\text{leavevmode}{s}_i(t+1)={\lambda }_i{s}_i(t)+u_i(t),$$

$$\text{leavevmode}\text{leavevmode}{\mathbf{B}}_{ij}({\theta }_i(t)-{\theta }_j(t))=x_i(t),j\in {\mathcal{N}}^e(i),$$

$$\text{leavevmode}\text{leavevmode}{\underline{\delta }}_i\le {\delta }_i(t)\le {\bar{\delta }}_i,\text{leavevmode}\text{leavevmode}{\underline{\mathsf{u}}}_i\le u_i(t)\le {\bar{\mathsf{u}}}_i,\text{leavevmode}\text{leavevmode}{\underline{\mathsf{s}}}_i\le s_i(t)\le {\bar{\mathsf{s}}}_i,$$

$$\text{leavevmode}\text{leavevmode}{\mathsf{b}}_j\le x_j(t)\le {\mathsf{c}}_j,\text{leavevmode}\text{leavevmode}j\in \{1,\cdots ,m\},$$

$$\{{\mathbf{z}}^{\star }(t),{\delta }^{\star }(t),{\mathbf{u}}^{\star }(t),{\mathbf{s}}^{\star }(t),{\theta }^{\star }(t){\}}_{t=1}^T=$$

$$\displaystyle\quad\operatorname{argmin}\leavevmode\nobreak\ \frac{1}{T}\sum\nolimits_{t=1}^{T}\!\Big{(}\sum\nolimits_{j=1}^{n}\!g_{j}(\delta_{j}(t))\,+$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \quad\quad\quad\quad\sum\nolimits_{i=1}^{m}\!\phi_{i}(\boldsymbol{\mathbf{z}}(t)^{\top}[\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}]_{i}+\mathsf{x}_{i}^{\text{p}}(\boldsymbol{\mathbf{\delta}}(t),\boldsymbol{\mathbf{u}}(t),\boldsymbol{\mathbf{\mathsf{d}}}(t))\Big{)},\leavevmode\nobreak\ \text{s.t.}$$

$$\text{forleavevmodenobreak }t\in \mathcal{T},\text{leavevmode}\text{leavevmode}i\in \{1,\dots ,n\}$$

$$\text{leavevmode}{\sum }_{j=1}^n({\delta }_j(t)+u_j(t)+{\mathsf{d}}_j(t))=0,$$

$$\text{leavevmode}{s}_i(t+1)={\lambda }_i{s}_i(t)+u_i(t),$$

$$\text{leavevmode}{\mathbf{B}}_{ij}({\theta }_i(t)-{\theta }_j(t))=\mathbf{z}(t{)}^{\top }[{\mathsf{B}}^{\text{of}}{]}_i+{\mathsf{x}}_i^{\text{p}}(\delta (t),\mathbf{u}(t),\mathsf{d}(t)),$$

$$j\in {\mathcal{N}}^e(i),$$

$$\text{leavevmode}{\underline{\delta }}_i\le {\delta }_i(t)\le {\bar{\delta }}_i,\text{leavevmode}\text{leavevmode}{\underline{\mathsf{u}}}_i\le u_i(t)\le {\bar{\mathsf{u}}}_i,\text{leavevmode}\text{leavevmode}{\underline{\mathsf{s}}}_i\le s_i(t)\le {\bar{\mathsf{s}}}_i,$$

$$\text{leavevmode}{\mathsf{b}}_j\le \mathbf{z}(t{)}^{\top }[{\mathsf{B}}^{\text{of}}{]}_j+{\mathsf{x}}_j^{\text{p}}(\delta (t),\mathbf{u}(t),\mathsf{d}(t))\le {\mathsf{c}}_j,$$

$$\text{leavevmode}j\in \{1,\cdots ,m\}.$$

$${\mathfrak{I}}^{\text{o}}\mathbf{x}(t)=\delta (t)+\mathbf{u}(t)+\mathsf{d}(t),t\in \mathcal{T}.$$

$$\displaystyle\big{\{}\boldsymbol{\mathbf{x}}\in^{m}\mid\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\boldsymbol{\mathbf{x}}=\boldsymbol{\mathbf{\delta}}+\boldsymbol{\mathbf{u}}+\boldsymbol{\mathbf{\mathsf{d}}},\boldsymbol{\mathbf{\delta}}\in^{n},\boldsymbol{\mathbf{u}}\in^{n}\}=$$

$$\displaystyle\Big{\{}\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}\boldsymbol{\mathbf{z}}+\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}(\boldsymbol{\mathbf{\delta}},\boldsymbol{\mathbf{u}},\boldsymbol{\mathbf{\mathsf{d}}})\,\Big{|}\,\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}(\boldsymbol{\mathbf{\delta}},\boldsymbol{\mathbf{u}},\boldsymbol{\mathbf{\mathsf{d}}})\!=\!\!\sum\nolimits_{i=1}^{n-1}\!\!(\delta_{i}\!+\!u_{i}\!+\!\mathsf{d}_{i})\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}},$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \sum\nolimits_{i=1}^{n}\!\!(\delta_{i}+u_{i}+\mathsf{d}_{i})=0,\,\boldsymbol{\mathbf{z}}\in\mathbb{R}^{n-m+1},\boldsymbol{\mathbf{\delta}}\in^{n},\boldsymbol{\mathbf{u}}\in^{n}\Big{\}}.$$

$${\theta }_i({\mathbf{y}}_i)=$$

$${\mathbf{y}}^{\star }=\underset{{\mathbf{y}}_1,\cdots ,{\mathbf{y}}_{\mu }}{\mathrm{argmin}}{\sum }_{i=1}^{\mu }{\theta }_i({\mathbf{y}}_i),\text{leavevmodenobreak s.t.}$$

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

$$\{\begin{array}{l}{\mathbf{y}}_i(e_k)=[\{{\bar{y}}_j{\}}_{j\in {\mathcal{C}}^{\text{of}}(e_i)}],\\ {\mathsf{b}}_k\le {\mathbf{y}}_i(e_k{)}^{\top }[\{{\mathsf{B}}_{jk}^{\text{of}}{\}}_{j\in {\mathcal{C}}^{\text{of}}(e_k)}]+{\mathsf{x}}_k^{\text{p}}\le {\mathsf{c}}_k,\end{array}\text{leavevmode}\forall e_k\in {\mathcal{E}}_i^{\mathfrak{C}}.$$

Full text

I Introduction

In a network flow problem, a physical system consisted of several routes between source and sink points transfers input flows from the source points to the sink points. The objective of optimal network flow problems mainly is to minimize the overall cost of transporting flow [1]. Network flow problems appear in many important applications, such as communication networks [2], wireless sensor networks [3], wireless routing and resource allocation [4], transportation systems [5] and power networks [6]. In power network problems, variants of optimal network flow problems also include optimal generation and storage costs in their objectives [7, 8, 9, 10].

With the advent of new technologies, the amount of available data and size of networks have been increasing. Such expansions in the size of physical networks result in increasing the size of optimization problems associated with optimal network flow problems. The number of decision variables has a direct relation with the time and space computation complexity of optimization problems. For large scale optimization problems, there has been efforts to use different variable reduction techniques to reduce problem size. Variable fixing techniques [11], dominance technique [12] and constraint pairing techniques [13] are some general reduction techniques in Integer Quadratic Problems (IQP). Moreover, in [14] a new variable reduction techniques for IQP proposed which fixes some decision variables at zero without loosing optimality. In multi-objective optimization problems also it is shown that using data mining techniques it is possible to reduce less effective variables [15]. For evolutionary optimization problems, [16] 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 to reduce the number of the search variables of the problem (c.f. [17])–

[TABLE]

where $\phi:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $g:\mathbb{R}^{n}\to\mathbb{R}^{m}$ are the cost function and the inequality constraint function, respectively, and $\boldsymbol{\mathbf{A}}\in\mathbb{R}^{p\times n}$ satisfies $\operatorname{rank}(\boldsymbol{\mathbf{A}})=\rho\leq p<n$. The affine feasible set for this optimization problem can be characterized as

[TABLE]

where $\boldsymbol{\mathbf{F}}\in\mathbb{R}^{n\times(n-\rho)}$ is a matrix whose columns expand the null-space of $\boldsymbol{\mathbf{A}}$ and $\boldsymbol{\mathbf{\mathsf{x}}}_{\text{p}}\in\mathbb{R}^{n}$ is a particular solution of $\boldsymbol{\mathbf{A}}\boldsymbol{\mathbf{x}}=\boldsymbol{\mathbf{b}}$. Then, $\boldsymbol{\mathbf{x}}^{\star}$ in (1) satisfies $\boldsymbol{\mathbf{x}}^{\star}=\boldsymbol{\mathbf{F}}\boldsymbol{\mathbf{z}}^{\star}+\boldsymbol{\mathbf{\mathsf{x}}}_{\text{p}}$ where

[TABLE]

Compared to (1), in (3) not only the equality constraint is eliminated but also the number of the search variables are reduced from $n$ to $n\!-\!\rho$. However, the lack of efficient methods to construct matrix $\boldsymbol{\mathbf{F}}$ and particular solution $\boldsymbol{\mathbf{x}}_{\text{p}}$ can be an impediment in use of affine equality constraint elimination method.

Optimal network flow problems are normally cast as a convex optimization problem where the cost is the sum of convex cost of flow through the arcs subject to capacity bounds for each arc and flow conservation equations at each node, resulting in optimization problems of the form (1). 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. The constraints can also be expanded to include other components of the problem. Nevertheless, in all network flow problems, an affine equality constraint that is always present is the flow conservation equation. To reduce the decision variables, one can use the aforementioned affine equality elimination approach, to eliminate the aggregated flow conservation equation from the network flow problems. In this paper, we discuss how the cycle basis structure from graph theory (c.f. [18]) can be used to accomplish this elimination in an efficient manner. Minimum cycle bases have applications in many areas such as electrical circuit theory [19], structural engineering [20], surface reconstruction [21].

In this paper, we show that all the solutions of the flow conservation equation is characterized explicitly in terms of the span of the columns of the transpose of the fundamental cycles basis matrix of the network plus a particular solution. The fundamental cycle basis of a graph can be computed in polynomial time using efficient algorithms such as those in [22, 23, 24] (see Appendix for a breif review). To compute a particular solution, 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 common particular sink node. We demonstrate our flow conservation equation elimination over a minimum cost flow problem as well as an optimal power flow problem with storage and generation at the nodes.

Parallel and distributed solutions are also sought as a method to solve large scale optimal network flow problems in an efficient manner. For example, a minimum cost network flow problem is solved in a distributed manner via dual sub-gradient descent in [1]. For the same problem, a distributed second order method with a better convergence rate is proposed in [25]. In [26], a distributed algorithm based on the local domain ADMM approach is proposed for minimum cost flow problem. These algorithms are all arc-based, i.e., to solve the network flow problem in a distributed manner, each arc or group of arcs are assigned to cyber-layer nodes. Then, the minimum cost network flow optimization problem is cast in a separable manner and solved by cyber-layer nodes in a cooperative way. Although in distributed algorithms the computational cost of the optimal flow problem is distributed among the cyber-layer nodes, the high number of decision variables normally translates to the high number of cyber nodes or large communication overhead between neighboring cyber nodes. Our next contribution in this paper is to show that the new formulation of the minimum cost network flow problem based on the cycle basis variables, which has a reduced set of search variables, is amenable to distributed solutions. Specifically, we demonstrate implantation of a distributed ADMM solution method (c.f. [27] and [28]) over this new formulation. To implement this distributed solution we propose a cyber-layer whose nodes are defined based on the fundamental cycles of a cycle basis of the physical-layer graph. A preliminary version of parts of our results in this paper has appeared in [29].

Notations: ${\mathbb{R}}$, ${\mathbb{R}}_{>0}$, ${\mathbb{R}}_{\geq 0}$, and ≤0 denote the set of real, positive real, non-negative real, and non-positive real numbers, respectively. We let $\boldsymbol{\mathbf{A}}^{\top}$ be the transpose of a matrix $\boldsymbol{\mathbf{A}}$ and $[\boldsymbol{\mathbf{A}}]_{i}$ indicate its $i^{\text{th}}$ column. We let $[\{z_{k}\}_{k=1}^{n}]$ be 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$, defined over $N$ nodes, $N$ arcs, or $N$ cycles, we represent the aggregate vector of these variables by $\boldsymbol{\mathbf{p}}=[\{p_{i}\}_{i=1}^{N}]\in{\mathbb{R}}^{N}$.

II A review of cycle basis in graphs

In this section, following [18], we review our graph related terminology and conventions. We also introduce our graph related notations. We represent a graph of $n$ nodes and $m$ arcs with $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{v_{1},v_{2},\cdots,v_{n}\}$ is the node set and $\mathcal{E}=\{e_{1},\cdots,e_{m}\}\in\mathcal{V}\times\mathcal{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 $\mathcal{G}=(\mathcal{V},\mathcal{E})$, we represent the oriented graph by $\mathcal{G}^{\text{o}}=(\mathcal{V},\mathcal{E}^{\text{o}})$. We write $e_{k}=(v_{i},v_{j})\in\mathcal{E}^{\text{o}}$ if arc $e_{k}$ points from node $v_{i}$ towards node $v_{j}$. If $(v_{i},v_{j})\in\mathcal{E}^{\text{o}}$ then $(v_{j},v_{i})\notin\mathcal{E}^{\text{o}}$, i.e., there is no symmetric arc in the oriented graph. For $\mathcal{G}^{\text{o}}$, the oriented incidence matrix is the matrix $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{E}|}$, where $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}_{ij}=1$ if arc $e_{j}$ leaves node $v_{i}$, $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}_{ij}=-1$ if arc $e_{j}$ enters node $v_{i}$, otherwise $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}_{ij}=0$. For a connected graph of $n$ nodes with a given orientation, the rank of $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}$ is $n-1$.

A cycle of $\mathcal{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 vector $\boldsymbol{\mathbf{\mathsf{c}}}\in^{m}$ is a binary vector with $\boldsymbol{\mathbf{\mathsf{c}}}_{i}=1$ if $e_{i}$ is in the cycle and $\boldsymbol{\mathbf{\mathsf{c}}}_{i}=0$, otherwise. A cycle basis of $\mathcal{G}$ is a set of simple cycles that forms a basis of the cycle space of $\mathcal{G}$. Every cycle in a given cycle basis is called a fundamental cycle. A fundamental cycle basis of a graph is constructed by its spanning tree, in a way that cycles formed by a combination of a path in the tree and a single arc outside of the tree. For every arc outside of the tree, there exist one cycle. Each cycle generated in this way is independent of other cycles, because it has one arc, not exist in other cycles (see Fig. 1). The dimension of cycle basis of a graph is $\mu=m-n+1$. For cycles in an oriented graph $\mathcal{G}^{\text{o}}=(\mathcal{V},\mathcal{E}^{\text{o}})$, we assign the counter clockwise direction as positive cycles orientation and define the oriented cycle vector $\boldsymbol{\mathbf{\mathsf{c}}}^{\text{o}}\in^{m}$ with $\boldsymbol{\mathbf{\mathsf{c}}}^{\text{o}}_{i}=1$ if $e_{i}$ is in the cycle and aligned with its direction, $\boldsymbol{\mathbf{\mathsf{c}}}^{\text{o}}_{i}=-1$ if $e_{i}$ is in the cycle but opposing the direction of the cycle and finally $\boldsymbol{\mathbf{\mathsf{c}}}^{\text{o}}_{i}=0$ if $e_{i}$ is not in the cycle. Given a cycle basis, we define the oriented fundamental cycle basis matrix $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}}\in^{\mu\times m}$ as a matrix whose rows are each the transpose of the oriented cycle vector of the fundamental cycles of this cycle basis. This matrix satisfies $\operatorname{rank}(\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}})=\mu$.

Theorem II.1** (relationship between the oriented incidence matrix and an oriented cycle vector  (c.f. [18])).**

In an oriented graph $\mathcal{G}^{\text{o}}$, every oriented cycle vector $\boldsymbol{\mathbf{\mathsf{c}}}^{o}$ is orthogonal to every row of oriented incident matrix $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}$, i.e., $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\,\boldsymbol{\mathbf{\mathsf{c}}}^{o}=\boldsymbol{\mathbf{0}}_{n}$. $\Box$

III Decision variable reduction in 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 through some form of routes. For example, in a power network the route is a transmission line, while in a transportation network the route is the road connecting two conjunction nodes on the road map. The physical layer topology is described by a connected graph $\mathcal{G}_{\text{physic}}=(\mathcal{V}_{\text{physic}},\mathcal{E}_{\text{physic}})$, where $|\mathcal{V}_{\text{physic}}|=n$ and $|\mathcal{E}_{\text{physic}}|=m$. The flow can travel in both directions in every route, however, we assume a pre-specified positive orientation for each route and based on it we describe the flow network in the physical layer by the oriented version of $\mathcal{G}_{\text{physic}}$, i.e., $\mathcal{G}^{\text{o}}_{\text{physic}}=(\mathcal{V}_{\text{physic}},\mathcal{E}^{\text{o}}_{\text{physic}})$. This physical network transfers flow(s) from a set of source nodes to a set of sink nodes (see physical layers in Fig. 3 and Fig. 4), 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. We let $x_{i}$ be the flow across the arc $e_{i}\in\mathcal{E}^{\text{o}}_{\text{physic}}=\{e_{1},e_{2},\cdots,e_{m}\}$. Every arc $e_{i}\in\mathcal{E}^{\text{o}}_{\text{physic}}$ has a pre-specified capacity, i.e., $\mathsf{b}_{i}\leq x_{i}\leq\mathsf{c}_{i}$, for some known $\mathsf{b}_{i},\mathsf{c}_{i}\in$. For any external flow $\mathsf{f}_{i}$, we use the sign convention of $\mathsf{f}_{i}>0$ for input flow and $\mathsf{f}_{i}<0$ for output flow, and $\mathsf{f}_{i}=0$ otherwise.

III-A Minimum cost flow problem

We consider a minimum cost flow problem over $\mathcal{G}^{\text{o}}_{\text{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}:\to$ associated with flow across each arc $e_{i}\in\mathcal{E}^{o}_{\text{physic}}$, and our objective is to find the network minimizer $\boldsymbol{\mathbf{x}}^{\star}\in^{m}$ in the following optimization problem

[TABLE]

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

Theorem III.1** (Eliminating the flow conservation constraint from (4)).**

Consider the optimal network flow problem (4) over a connected physical network $\mathcal{G}_{\text{\rm{physic}}}^{\text{o}}$. Then, $\boldsymbol{\mathbf{x}}^{\star}$ in (1) satisfies $\boldsymbol{\mathbf{x}}^{\star}=\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}\boldsymbol{\mathbf{z}}^{\star}+\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}$ where

[TABLE]

and $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}$ is an oriented fundamental cycle matrix of $\mathcal{G}^{\text{o}}_{\text{physic}}$ and $\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}$ is a particular solution of $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\,\boldsymbol{\mathbf{x}}=\boldsymbol{\mathbf{\mathsf{f}}}$.

Proof.

The equality constraint (4b) in aggregated form is

[TABLE]

Invoking the same argument that is used to relate solutions of the optimization problem (1) to those of (3), the proof relays on showing that the null-space of $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}$ is spanned by columns of $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}$. By virtue of Theorem II.1, we have $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}=\boldsymbol{\mathbf{0}}$. Recall that for a connected oriented graph $\operatorname{rank}(\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}})=n-1$. Because $\operatorname{rank}(\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top})=m-n+1$, null-space of $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\in^{n\times m}$ is spanned by columns of $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}$. This completes our proof. ∎

The effectiveness of the decision variable reduction method in Theorem III.1 depends on how efficiently one can construct matrix $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}$ and particular solution $\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}$, especially in large scale networks. In regards to matrix $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}$, as reviewed in Appendix, there are efficient algorithms that can construct cycle basis in polynomial time. Next, we propose a simple method to construct a particular solution $\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\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 $\mathsf{f}_{i}=1$ entering the network at node $v_{i}$ and leaving it at node $v_{j}$ can simply be constructed by assuming that $\mathsf{f}_{i}=1$ flows along a path from node $v_{i}$ to node $v_{j}$.

Lemma III.1** (Particular solution of (6)).**

Given an input/output flow vector $\boldsymbol{\mathbf{\mathsf{f}}}$ over $\mathcal{G}_{\text{physic}}^{\text{o}}$ which satisfies $\sum_{i=1}^{n}\mathsf{f}_{i}=0$, a particular solution for (6) is $\boldsymbol{\mathbf{\mathsf{x}}}^{\text{p}}=\sum_{i=1}^{n-1}\mathsf{f}_{i}\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}$. Here, $\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}\in^{m}$, $i\in\{1,\dots,n-1\}$, is constructed from a path that connects node $v_{i}$ to node $v_{n}$ such that $\mathsf{x}^{\text{p},v_{i}}_{j}=1$ (resp. $\mathsf{x}_{j}^{\text{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 $\mathsf{x}_{j}^{\text{p},v_{i}}=0$.

Proof.

For every $i\in\{1,\dots,n-1\}$, consider a virtual scenario where a unit flow $\bar{\mathsf{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{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}\in^{m}$ as described in the statement satisfies $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}=\bar{\boldsymbol{\mathbf{\mathsf{f}}}}^{v_{i}}$, $i\in\{1,\dots,n-1\}$ where $\bar{\mathsf{f}}^{v_{i}}_{i}=1$, $\bar{\mathsf{f}}^{v_{i}}_{n}=-1$ and $\bar{\mathsf{f}}^{v_{i}}_{j}=0$, $j\in\{1,\dots,n\}\backslash\{i,n\}$. For a given network flow vector $\boldsymbol{\mathbf{\mathsf{f}}}$ because $\mathsf{f}_{n}=-\sum_{i=1}^{n-1}\mathsf{f}_{i}$, we can write $\boldsymbol{\mathbf{\mathsf{f}}}=\sum_{i=1}^{n-1}\mathsf{f}_{i}\,\bar{\boldsymbol{\mathbf{\mathsf{f}}}}^{v_{i}}$. Therefore, $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\sum\nolimits_{i=1}^{n-1}\mathsf{f}_{i}\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p}}_{v_{i}}=\sum\nolimits_{i=1}^{n-1}\mathsf{f}_{i}\,\bar{\boldsymbol{\mathbf{\mathsf{f}}}}^{v_{i}}=\boldsymbol{\mathbf{\mathsf{f}}},$ which completes our proof. ∎

A few remarks are in order regarding the particular solution. First, note that construction of the ‘elementary’ particular solution set $\{\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}\}_{i=1}^{n-1}$ is regardless of the value of the network flow vector $\boldsymbol{\mathbf{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 (see Section IV-B for an illustrative numerical example). Finally, to obtain sparse elementary solutions, we can use a shortest path between nodes $v_{i}$ and $v_{n}$.

III-B Optimal power flow with storage and generation at nodes

Next, we consider an optimal power flow problem over a network described by $\mathcal{G}^{\text{o}}_{\text{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 $\mathcal{T}=\{1,\cdots,T\}$. Mathematical modeling of this problem over various scenarios including deterministic and stochastic generators is considered in the literature, e.g., [7, 8, 9, 10]. All these models, at each time $t\in\mathcal{T}$, include a flow conservation equation at each node. As a result for a network with $m$ arcs, the flow conservation equation introduces $m\,|\mathcal{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 [7], which states the problem as a direct current (DC) power flow problem (see (7)). Without loss of generality, we assume that at each time $t\in\mathcal{T}$, each node $v_{i}\in\mathcal{V}_{\text{physic}}$ has a generator which 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 $\mathsf{d}_{i}(t)\in_{\leq 0}$. If a node does not have either of the generation, storage, or load components, we simply remove the respective variables from our formulation below. Then, given a known load profile $\{\boldsymbol{\mathbf{\mathsf{d}}}(t)\}_{t=1}^{T}$, where $\boldsymbol{\mathbf{\mathsf{d}}}(t)=[\{\mathsf{d}_{i}(t)\}_{i=1}^{n}$], the optimization problem of interest is

[TABLE]

where, $g_{j}$, $j\in\{1,\dots,n\}$ and $\phi_{i}$, $i\in\{1,\dots,m\}$ are cost functions for generators and power flows, respectively. Here, $\mathcal{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{\mathsf{s}}_{i},\bar{\mathsf{s}}_{i})\in\times$, $(\underline{\mathsf{u}}_{i},\bar{\mathsf{u}}_{i})\in\times$, $(\underline{\mathsf{\delta}}_{i},\bar{\mathsf{\delta}}_{i})\in\times$, are, respectively, known (lower bound, upper bound) values on storage level, battery charge/discharge and power generation by the generator. Finally $\boldsymbol{\mathbf{B}}_{ij}(\theta_{i}(t)-\theta_{j}(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\mathcal{V}_{\text{physic}}$ at time $t$ and $\boldsymbol{\mathbf{B}}\in{\mathbb{R}}^{n\times n}$ is the imaginary part of the admittance matrix under DC assumption.

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

Proposition III.1** (Eliminating the flow conservation constraint from (7)).**

Consider the optimal power flow problem (7) over a physical network described by $\mathcal{G}_{\text{\rm{physic}}}^{\text{o}}$ with a given set of loads $\{\boldsymbol{\mathbf{\mathsf{d}}}(t)\}_{t=1}^{T}$. Then, $\{\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}\boldsymbol{\mathbf{z}}^{\star}(t)+\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}(\boldsymbol{\mathbf{\delta}}^{\star}(t),\boldsymbol{\mathbf{u}}^{\star}(t),\boldsymbol{\mathbf{\mathsf{d}}}(t)),\boldsymbol{\mathbf{u}}^{\star}(t),\boldsymbol{\mathbf{s}}^{\star}(t),\boldsymbol{\mathbf{\theta}}^{\star}(t)\}_{t=1}^{T}$ is a minimizer of the optimization problem (7) where

[TABLE]

Here, $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}$ is a fundamental cycle basis matrix of $\mathcal{G}^{\text{o}}_{\text{physic}}$ and $\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}(\boldsymbol{\mathbf{\delta}}(t),\boldsymbol{\mathbf{u}}(t),\boldsymbol{\mathbf{\mathsf{d}}}(t))=\sum_{i=1}^{n-1}(\delta_{i}(t)+u_{i}(t)+\mathsf{d}_{i}(t))\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}$, where $\{\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}\}_{i=1}^{n-1}$ is as described in Lemma III.1.

Proof.

The equality constraint (7b) in aggregated form is

[TABLE]

Note that rank of $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\in^{n\times m}$ is $n-1$ and also that $\boldsymbol{\mathbf{1}}_{n}^{\top}\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}=\boldsymbol{\mathbf{0}}$, where $\boldsymbol{\mathbf{1}}_{n}$ is the vector of $n$ ones. Left multiplying (9) by $\boldsymbol{\mathbf{1}}_{n}^{\top}$ results in (8b). Then, for $\boldsymbol{\mathbf{\delta}}(t),\boldsymbol{\mathbf{u}}(t),\boldsymbol{\mathbf{\mathsf{d}}}(t)\in^{n}$ that satisfy  (8b), following the method discussed in Lemma III.1, we can show that $\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}(\boldsymbol{\mathbf{\delta}}(t),\boldsymbol{\mathbf{u}}(t),\boldsymbol{\mathbf{\mathsf{d}}}(t))=\sum_{i=1}^{n-1}(\delta_{i}(t)+u_{i}(t)+\mathsf{d}_{i}(t))\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{i}}$ is a particular solution of (9), i.e., $\boldsymbol{\mathbf{\mathfrak{I}}}^{\text{o}}\,\boldsymbol{\mathbf{\mathsf{x}}}^{\text{\rm{p}}}(\boldsymbol{\mathbf{\delta}}(t),\boldsymbol{\mathbf{u}}(t),\boldsymbol{\mathbf{\mathsf{d}}}(t))=\boldsymbol{\mathbf{\delta}}(t)+\boldsymbol{\mathbf{u}}(t)+\boldsymbol{\mathbf{\mathsf{d}}}(t).$ Therefore, for any given load vector $\boldsymbol{\mathbf{\mathsf{d}}}$ (to simplify the notation we drop argument $t$), we have

[TABLE]

As a result, at each $t\in\mathcal{T}$, we can eliminate the affine equation (7b) and arrive at the equivalent optimization problem (8) whose minimizers are related to minimizers of (7) in a way that is described in the statement. ∎

IV A cycle-basis distributed ADMM algorithm for minimum cost network flow problem

In this section, we consider the minimum cost network flow problem (4) and show that its equivalent form (5) based on the cycle basis variables is amenable to a distributed solution. We start by introducing some notation related to the oriented fundamental cycles of $\mathcal{G}^{\text{o}}$. Let $\mathcal{C}^{\text{of}}=\{\mathfrak{C}^{\text{of}}_{i}\}_{i=1}^{\mu}$, where $\mu=m-n+1$, be the set of fundamental cycles of $\mathcal{G}^{\text{o}}$ whose cycle matrix $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}}$ is used to eliminate the flow conservation equation as explained in Section III. We represent the set of arcs of any $\mathfrak{C}^{\text{of}}_{i}\in\mathcal{C}^{\text{of}}$, $i\in\{1,\dots,\mu\}$, by $\mathcal{E}^{\mathfrak{C}}_{i}=\{e_{j}\in\mathcal{E}^{\text{o}},j\in\{1,\dots,m\}\,|\,\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}}_{ij}\neq 0\}$. For a given cycle basis $\mathcal{C}^{\text{of}}$, we refer to the cycles that share an arc as neighbors and represent the set of (cycle) neighbors of any fundamental cycle $\mathfrak{C}^{\text{of}}_{i}\in\mathcal{C}^{\text{of}}$, $i\in\{1,\dots,\mu\}$, by $\mathcal{N}^{\mathfrak{C}}_{i}=\{j\in\{1,\dots,\mu\}\backslash\{i\}\,|\,\exists\,k\in\{1,\dots,m\}\text{\leavevmode\nobreak\ s.t.\leavevmode\nobreak\ }\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}}_{ik}\neq 0\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}}_{jk}\neq 0\}$. We let $\mathcal{C}^{\text{of}}(e_{i})$ be the set of indexes of the fundamental cycles that arc $e_{i}\in\mathcal{E}^{\text{\text{o}}}$ belongs to them, i.e., $\mathcal{C}^{\text{of}}(e_{i})=\{j\in\{1,\dots,\mu\}\,|\,e_{i}\in\mathcal{E}^{\mathfrak{C}}_{j}\}$.

Recall that to eliminate the flow conservation equation we used $\boldsymbol{\mathbf{x}}=\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}\,{}^{\top}\boldsymbol{\mathbf{z}}+\boldsymbol{\mathbf{\mathsf{x}}}^{\text{p}}$, or equivalently $x_{i}=\boldsymbol{\mathbf{z}}^{\top}[\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}]_{i}+\mathsf{x}_{i}^{\text{p}}$, $i\in\{1,\dots,m\}$. Notice that one can think of every $z_{i}$, $i\in\{1,\dots,\mu\}$ as a cycle flow variable (with positive direction in counterclockwise direction) of the fundamental cycle $\mathfrak{C}^{\text{of}}_{i}$. Recall that, for a given arc $e_{i}\in\mathcal{E}^{\text{o}}_{\text{physic}}$, every element of $\boldsymbol{\mathbf{\mathsf{B}}}^{\text{of}}_{ji}$ is zero except if cycle $\mathfrak{C}^{\text{of}}_{j}$ contains arc $e_{i}$, i.e., $e_{i}\in{\mathcal{E}}^{\mathfrak{C}}_{j}$. As a result, we can deduce that every $x_{i}$, $i\in\{1,\dots,m\}$, is an affine function of $\mathsf{x}_{i}^{\text{p}}$ and $\{z_{k}\}_{k\in\mathcal{C}^{\text{of}}(e_{i})}$, indicating that every arc flow is a function of its particular solution and the cycle flow of fundamental cycles that contain the arc. Given such relationship, then the cost function of every arc as $\phi_{i}(x_{i})=\phi_{i}(\boldsymbol{\mathbf{z}}^{\top}[\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}]_{i}+\mathsf{x}_{i}^{\text{p}})=\psi_{i}(\{z_{k}\}_{k\in\mathcal{C}^{\text{of}}(e_{i})})$.

Cyber layer architecture: based on the observation above, we propose to assign a cyber-layer node to each fundamental cycle (see Fig. 3 as an example). We assume that the cyber-layer nodes of neighboring fundamental cycles can communicate with each other in bi-directional way. For bi-connected physical layer graphs this procedure will result in a connected graph of $\mu$ nodes for cyber layer (see Fig. 3 and Fig. 4 for examples). To obtain fundamental cycles with fewest number of arcs we propose to use minimum weight cycle basis algorithms to generate the fundamental cycle basis for the cyber layer (see Appendix).

IV-A Cycle Basis distributed ADMM solution for minimum cost network flow problem

In this section, we derive an equivalent representation of optimization algorithm (5) which can be solved in a distributed manner using the ADMM algorithm of [26] by our cycle-based cyber layer. To this end, for every cyber-layer node $i\in\{1,\dots,\mu\}$, we define $\boldsymbol{\mathbf{y}}_{i}=(\bar{y}_{i},\tilde{\boldsymbol{\mathbf{y}}}_{i})\in^{|\mathcal{N}^{\mathfrak{C}}_{i}|+1}$, where $\bar{y}_{i}\in$ is the local copy of $z_{i}$ and $\tilde{\boldsymbol{\mathbf{y}}}_{i}$ is the local copy of $\{z_{k}\}_{k\in\mathcal{N}^{\mathfrak{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 its neighbors. Next, we cast the cost function of each cyber node in terms of its decision variable $\boldsymbol{\mathbf{y}}_{i}$. Let $\boldsymbol{\mathbf{y}}_{i}(e_{k})$ be the component(s) of $\boldsymbol{\mathbf{y}}_{i}$ corresponding to $\{z_{j}\}_{\{j\in\mathcal{C}^{\text{of}}(e_{k})\}}$. For every cyber-layer node, we define its cost function as

[TABLE]

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

[TABLE]

In this formulation, every cycle-based cyber node has a copy of the cycle flows that go through its arcs, i.e, $\boldsymbol{\mathbf{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 (11) fits the standard framework developed for distributed ADMM solutions and can be solved for example using the algorithm of [26]. The details are omitted for brevity. In this distributed implementation we assume that elementary particular solution set $\{\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{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 $\boldsymbol{\mathbf{\mathsf{f}}}$ to the cyber-layer agents. For networks with large fundamental cycle sizes, one can split a cycle among several cyber nodes. In this case the length of the $\bar{y}_{i}$ of these agents will be [math] and we can still use the distributed ADMM algorithm to solve the problem. Similarly, two or more cycles can be assigned to one cyber layer.

IV-B Numerical Example

In this section, we demonstrate the use of distributed cycle-based ADMM algorithm for a minimum cost optimal flow over the network shown in Fig. 1. We assign positive flow orientation to the arcs as represented in Fig. 4. We generate the cyber layer based on the minimum weight cycle basis as shown in the bolder network with gray nodes in Fig. 1. In this problem, we set $\mathsf{b}_{i}=-\,\mathsf{c}_{i}$, and $\mathsf{c_{i}}\in{\mathbb{R}}_{>0}$, $i\in\{1,\dots,18\}$. We assume that the cost of the network flow at each arc is given as $\phi_{i}(x_{i})=(\frac{x_{i}}{\mathsf{c}_{i}})^{2}$, where $x_{i}=\boldsymbol{\mathbf{z}}^{\top}[\boldsymbol{\mathbf{\mathsf{B}}}^{\text{\rm{of}}}]_{i}+\mathsf{x}_{i}^{\text{p}}$ is the arc flow and $\boldsymbol{\mathbf{z}}=(z_{1},\cdots,z_{8})^{\top}$ are cycle flows. In the physical layer network in Fig. 4, there are two source nodes $v_{1}$ and $v_{4}$ and two two sink nodes $v_{9}$ and $v_{11}$.

We follow Lemma III.1 to generate a particular solution for given input output flow vector $\boldsymbol{\mathbf{\mathsf{f}}}=(\mathsf{f}_{1},0,0,\mathsf{f}_{4},0,0,0,0,\mathsf{f}_{9},0,-(\mathsf{f}_{1}+\mathsf{f}_{4}+\mathsf{f}_{9}))^{\top}$– recall that $\mathsf{f}_{11}=-(\mathsf{f}_{1}+\mathsf{f}_{4}+\mathsf{f}_{9})$ (recall that input flows have positive and output flows have negative values). We compute the elementary particular solutions for nodes $v_{1}$, $v_{4}$ and $v_{9}$ using shortest path from them to node $v_{11}$: $\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{1}}=(0,0,1,\boldsymbol{\mathbf{0}}_{1\times 4},1,0,1,\boldsymbol{\mathbf{0}}_{1\times 7},1)^{\top}$, $\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{4}}=(\boldsymbol{\mathbf{0}}_{1\times 7},1,0,1,\boldsymbol{\mathbf{0}}_{1\times 7},1)^{\top}$, and $\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{9}}=(\boldsymbol{\mathbf{0}}_{1\times 16},1,0)^{\top}$. Then $\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p}}=\mathsf{f}_{1}\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{1}}+\mathsf{f}_{4}\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p},v_{4}}+\mathsf{f}_{9}\,\bar{\boldsymbol{\mathbf{\mathsf{x}}}}^{\text{p}}$.

In our simulation, the lower and upper capacity bounds are selected uniformly randomly from $[2,50]$, i.e., $\mathsf{c}_{i}\in[2,50]$, $i\in\{1,\dots,m\}$. Specifically, at arcs connected to sink and source points we have $\mathsf{c}_{1}=30$, $\mathsf{c}_{2}=2$, $\mathsf{c}_{3}=27$, $\mathsf{c}_{8}=50$, $\mathsf{c}_{13}=49$, $\mathsf{c}_{14}=23$, $\mathsf{c}_{15}=21$, $\mathsf{c}_{16}=11$, $\mathsf{c}_{17}=16$, and $\mathsf{c}_{18}=37$. For our selected capacity bounds, using Edmonds-Krap algorithm [30] we obtain the maximum input flow $\mathsf{f}_{1}+\mathsf{f}_{4}$ to be $82$. In our simulation, then we set $\mathsf{f}_{1}=32$ and $\mathsf{f}_{4}=50$. The output follows are $\mathsf{f}_{9}=-52$ and $\mathsf{f}_{11}=-30$. After, $50$ iteration the input/output flows change to $\mathsf{f}_{1}=15$, $\mathsf{f}_{4}=30$, $\mathsf{f}_{9}=-15$, and $\mathsf{f}_{11}=-30$.

We use Matlab ‘quadprog’ to solve the problem in a centralized manner to generate reference values to compare the performance of our distributed cycle basis distributed ADMM algorithm as outlined in Section IV. The results are depicted in Fig. 5. In Fig. 5, plots show the distance of arc flow values from their optimum solution during execution of distributed ADMM. During the first 50 iteration the distributed ADMM converges to the optimum solution. Then, for the second external flows, it converges to the optimum solution again.

V Conclusion and future works

We considered optimal network flow problems and investigated how the decision variables 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 flow conservation equation in an efficient manner. In particular, we showed that the computation regarding the proposed variable reduction can be done in a systematic manner, in polynomial time, using existing algorithms. Moreover, we showed that the new formulation of the optimal network flow problems with reduced variables is amenable to distributed solvers. In this regard, we constructed a cyber-layer structure based on cycles in the physical-layer network. We also demonstrated the use of a distributed ADMM solver for minimum cost flow problem.