Distributed Algorithm for Real-Time Energy Optimal Routing Based on Dual Decomposition of Linear Programming

Abstract

This work proposes a novel in-network distributed algorithm for real-time energy optimal routing in ad hoc and sensor networks for systems with linear cost functions and constant communication delays. The routing problem is described as a minimum-cost multicommodity network flow problem by linear programming and modified by network replication to a real-time aware form. Based on the convex programming theory we use dual decomposition to derive the distributed algorithm. Thanks to the exact mathematical derivation, the algorithm computes the energy optimal real-time routing. It uses only peer-to-peer communication between neighboring nodes and does not need any central node or knowledge about the whole network structure. Each node knows only the produced and collected data flow and the costs of its outgoing communication links. According to our knowledge, this work is the first, which solves the real-time routing problem with linear cost functions and constant communication delays, using the dual decomposition.

1. Introduction

Our work is focused on a distributed algorithm for data flow routing through multihop ad hoc and sensor networks, where all data has to be delivered to the destinations in time. An example of a target application is a network periodically sensing some consumption variables (like gas consumption, water consumption, etc.) in large objects. Each sensing device produces a data flow of a particular volume, which is supposed to be routed through the network. The objective is to optimize the energy consumption for the data transfer (minimal possible consumption in the whole network), while constrained by the communication capacities (maximum data volume which can be transferred by a device per time unit) and by maximum communication delay.

We assume a time division multiple access (TDMA) protocol (e.g., GTS allocation in IEEE 802.15.4 [1]) which ensures collision-free communication and causes constant communication delay. This approach is used quite often in industry and time-critical systems. Due to the TDMA mechanism assumed, the worst-case delay from the source node to the sink node is the sum of the particular delays for each of the hops, assumed to be an integer (derived from the parameters like TDMA period, worst-case execution time of the communication stack…). In a particular setting, we may assume a unit hop delay (e.g., the same TDMA period). In this paper we assume the unit hop delay (the deadlines are expressed as the number of communication hops between devices) which is very transparent for the reader.

There are many communication protocols designed to find the exact energy optimal data routing in wireless sensor networks. However, to comply with the communication capacities, they usually need a central computational point with the knowledge of the actual network structure and parameters (e.g., [2]). The existence of such a computational point decreases the robustness of the system against the network damage. Furthermore, the routing of information about the actual network structure and parameters has to be solved in the case of the centralized algorithm.

In this paper, we propose a distributed algorithm, which computes the energy optimal real-time data routing without the need of any central computational or data point. The algorithm supposes that each node knows only the cost (energy consumption per transmitted data unit) of its outgoing communication links and the data which is supposed to send and receive. The whole algorithm is mathematically derived using the convex programming theory. The main purpose of this paper is to present the principle of new distributed routing algorithms rather than to present an application-ready algorithm. We believe that the presented approach can lead to new, efficient, and highly adaptive routing algorithms for ad hoc and sensor networks. Moreover, the approach used in this work can be adapted for general distribution of convex optimization problems into the network.

2. Related Works

Traditionally, the routing problems for data networks are formulated as linear or convex multicommodity network flow routing problems, for example [2–4]. One of the advantages of this method is that several cost functions and constraints can be put together (e.g., different types of capacity and energy consumption constraints and real-time constraints). Using the same underlying model, we can easily combine the solution of different works focused on partial problems.

Several papers have been performed in the area of real-time routing in multihop wireless sensor networks. In [5], a well-known soft real-time communication protocol SPEED is presented. Several works use the relation between the message propagation speed and transmitting energy to balance the trade-off between the energy consumption and communication delay [6–8]. In [9] the authors deal with real-time communications over cluster-tree sensor networks, where they evaluate the end-to-end communication delay. However, none of these algorithms can ensure real-time and energy optimal data routing especially in high-loaded networks.

There are several works, which focus on the decomposition of network problems described by strictly convex optimization. A systematic presentation and classification of the decomposition techniques for network utility maximization (NUM) is presented in [10–12].

In [13–15], the authors use dual decomposition to decompose cross-layer optimization problems into the optimization of separated layers. The presented approaches lead to structural decomposition (e.g., to the routing layer, the capacity layer…) which is not suitable for the derivation of the in-network distributed algorithm. In [16], a general distributed algorithm for a strictly convex optimization problem with a common parameter for all nodes is presented. The decomposition of an optimal routing problem is presented, for example, in [17, 18], where the authors focus on the node-path formulation and based on dual decomposition derive an algorithm suitable for networks with a small number of communication paths.

In [19, 20], the authors use the duality theorem and node-link problem formulation to analytically derive the maximum network lifetime for linear networks. In [21], the authors use the node-link problem formulation and derive the distributed routing algorithm with an extension for the queuing delay. However, all these algorithms are limited to strictly convex cost functions and fail in the case of linear cost functions.

This paper also continues our previous work. It joins and extends the results from [22, 23]. In [22], we have derived a distributed in-network routing algorithm based on dual decomposition. The algorithm consists of one iteration loop and is even able to solve problems with a linear objective function. In [23], we have presented a centralized algorithm for energy optimal real-time routing, which is represented as a min-cost multicommodity network flow routing problem with side constraints.

2.1. Contribution and Outline

The main contributions of this paper are

• (1)

introduction of a new mathematically derived, distributed algorithm for energy optimal real-time routing based on network replication and dual decomposition;

(2)

introduction of a new approach for an in-network distribution of real-time routing problems with constant communication delay and linear objective function;

(3)

application of a distribution procedure for routing problems presented in [22] on a real-time routing problem [23].

The paper is organized as follows. Section 3 briefly describes the multicommodity network flow model. In Section 4 we define the real-time routing model. In Section 5, the distributed algorithm and its derivation are presented. An example and computational complexity experiments are given in Section 6. Section 7 concludes the paper and mentions the future work. The proof of the algorithm convergence is presented in the Appendix.

3. Multicommodity Network Flow Model

In this section, we briefly summarize the basic terminology and specify the multicommodity network flow model. For more details see, for example, [3, 4].

The network is represented by an oriented graph, where for each device able to send or receive data, a node of the graph exists. The nodes are labeled as $n=1,\dots ,N$ . Directed communication links are represented as ordered pairs $(n_1,n_2)$ of distinct nodes. The links are labeled as $l=1,\dots ,L$ . We define the set of the links l leaving the node n as $𝒪(n)$ and the set of the links l incoming to node n as $ℐ(n)$ . The network structure is described with two incidence matrices in node-link form. The matrix of the incoming links is denoted as $A^{+}$ , and the matrix of the outgoing links is denoted as $A^{-}$ :

$\begin{array}{cc}\hfill A_{n,l}^{+}& =\{\begin{array}{cc}1,\hfill & \hspace{0.17em}\hspace{0.17em}l\in \mathcal{I}\left(n\right)\hspace{0.17em}\hspace{0.17em}\left(\text{link}\hspace{0.17em}\hspace{0.17em}l\hspace{0.17em}\hspace{0.17em}\text{enters}\hspace{0.17em}\hspace{0.17em}\text{node}\hspace{0.17em}\hspace{0.17em}n\right),\hfill \\ 0,\hfill & \hspace{0.17em}\hspace{0.17em}\text{otherwise},\hfill \end{array}\hfill \\ \hfill A_{n,l}^{-}& =\{\begin{array}{cc}1,\hfill & \hspace{0.17em}\hspace{0.17em}l\in \mathcal{O}\left(n\right)\hspace{0.17em}\hspace{0.17em}\left(\text{link}\hspace{0.17em}\hspace{0.17em}l\hspace{0.17em}\hspace{0.17em}\text{leaves}\hspace{0.17em}\hspace{0.17em}\text{node}\hspace{0.17em}\hspace{0.17em}n\right),\hfill \\ 0,\hfill & \hspace{0.17em}\hspace{0.17em}\text{otherwise}.\hfill \end{array}\hfill \end{array}$

(1)

By m we denote an index of the communication demand and by ℳ we denote a set of all communication demands. The communication demands can be seen as the flow of various commodities incoming/leaving the network in some nodes. The flow of each communication demand has to satisfy the flow conservation law at each node (for a given commodity, the sum of the flow incoming to the node is equal to the sum of the flow leaving the node)

$\begin{array}{c}{A}^{-}{\vec{x}}^{\hspace{0.17em}\left(m\right)}+{\vec{s}}_{\text{out}}^{\hspace{0.17em}\left(m\right)}=A^{+}{\vec{x}}^{\hspace{0.17em}\left(m\right)}+{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m\right)}\hspace{1em}\forall m\in \mathcal{M},\end{array}$

(2)

where the column vector ${\vec{s}}_{\text{in}}^{\hspace{0.17em}(m)}\ge \vec{0}$ denotes the flow coming into the network, the ${\vec{s}}_{\text{out}}^{\hspace{0.17em}(m)}\ge \vec{0}$ denotes the flow leaving the network, and the ${\vec{x}}^{\hspace{0.17em}(m)}\ge \vec{0}$ denotes the flow routed through the network for demand m. Notice that a multisource multisink problem can be described in this way (e.g., the data gathering problem).

In this work, the node capacity constraints are used. They are the most common capacity constraint used for wireless networks. The total volume of the flow leaving one node has to satisfy the capacity constraint $D{\sum }_{m\in ℳ}{\vec{x}}^{(m)}\le \vec{\mu }$ , where $\vec{\mu }\ge \vec{0}$ is the column vector of the node capacities and matrix D describes the constraint structure (i.e., $D=A^{-}$ ).

In summary, the network flow model imposes the following constraints on the network flow variables ${\vec{x}}^{(m)}$ :

$\begin{array}{c}{A}^{-}{\vec{x}}^{\hspace{0.17em}\left(m\right)}+{\vec{s}}_{\text{out}}^{\hspace{0.17em}\left(m\right)}=A^{+}{\vec{x}}^{\hspace{0.17em}\left(m\right)}+{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m\right)}\hspace{1em}\forall m\in \mathcal{M},\\ D{\displaystyle \sum }_{m\in \mathcal{M}}{\vec{x}}^{\left(m\right)}\le \vec{\mu },\\ {\vec{x}}^{\left(m\right)}\ge \vec{0}\hspace{1em}\forall m\in \mathcal{M}.\end{array}$

(3)

The task of the total energy minimization is to minimize the cost function $f_{\text{cost}}={\vec{c}}^{\hspace{0.17em}T}{\sum }_{m\in ℳ}{\vec{x}}^{(m)}$ by setting the flow vector ${\vec{x}}^{(m)}$ for all $m\in ℳ$ subject to the system of inequalities (3). Vector $\vec{c}>0$ is a column vector of the energy consumption per data unit transmitted. The components of vector $\vec{c}$ correspond with the energy consumption for the individual links in the network. Their values are usually determined directly from the transmission energy level in the nodes, which is needed for a reliable connection. This is done by the MAC layer of the communication protocol.

Another target application of this approach is a network, where the energy supplies of the nodes can be recharged with different maintenance costs (e.g., depending on the node reachability). In such an application, the recharge price is included in vector $\vec{c}$ , and we minimize the network long-time maintenance cost.

4. Real-Time Multicommodity Network Flow Model

In this section, we extend the multicommodity flow model by real-time constraints, which guarantee a controllable routing delay through the network. Each communication demand has its own deadline, and the communication delay of each demand has to be shorter than its deadline. We model the hop delay as an integer value associated with each communication link. For transparent model derivation, we assume the same communication delay over the entire network (i.e., each communication hop causes a delay equal to one). However, the model can be extended to a more general form, where the communication delays are integers.

4.1. Mathematical Model of Real-Time Routing

Let vector ${\vec{x}}^{(m,w)}\in R_{+}^L$ denote the flow of the communication demand m with an integer communication delay w in the network. Then the flow vector ${\vec{x}}^{(m)}$ , independent of the flow delay of demand m, is equal to the sum of the flow vectors over all acceptable delays: ${\vec{x}}^{(m)}={\sum }_{w=1}^{d^{(m)}}{\vec{x}}^{(m,w)}$ , where $d^{(m)}$ denotes the deadline of the communication demand m. Using this equation, we can rewrite the capacity constraint from (3) into a new form:

$\begin{array}{c}D{\displaystyle \sum }_{m\in \mathcal{M}}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\vec{x}}^{\left(m,w\right)}\le \vec{\mu }.\end{array}$

(4)

Vector ${\vec{s}}^{\hspace{0.17em}(m,w)}\in R_{+}^N$ stands for the flow of the demand m leaving the network with the communication delay w, and vector ${\vec{s}}_{\text{in}}^{\hspace{0.17em}(m,w)}\in R_{+}^N$ denotes the flow of demand m coming into the network with the initial delay w. The flow of each demand may come into the network and leave it in more nodes. If all flow fragments of one demand coming into the network in different nodes have the same initial delay, then ${\vec{s}}_{\text{in}}^{\hspace{0.17em}(m,0)}={\vec{s}}_{\text{in}}^{\hspace{0.17em}(m)}$ , and ${\vec{s}}_{\text{in}}^{\hspace{0.17em}(m,w)}=0$ for $w>0$ . The flow of demand m leaving the network prior to the deadline is

$\begin{array}{c}{\vec{s}}_{\text{out}}^{\hspace{0.17em}\left(m\right)}={\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}\hspace{1em}\forall m\in \mathcal{M}.\end{array}$

(5)

Through (4) and (5), we have converted the real-time constraint (i.e., the delay has to be shorter than the deadline) to the structural constraint. Only the flow, whose delay is shorter than the deadline, is represented. The flow, which does not meet the deadline, violates the flow conservation law, and then the network flow constraints are not satisfied; that is, this solution is not feasible.

If the flow is sent through the network, the flow delay is increased by each communication hop. The flow of demand m coming into node n with communication delay w has to either leave the network in node n with the same delay w or reach the neighbor node with delay $w+1$ . The flow conservation law from (2) can be rewritten in the delay awareness form as

$\begin{array}{c}{A}^{-}{\vec{x}}^{\left(m,w+1\right)}+{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}=A^{+}{\vec{x}}^{\left(m,w\right)}+{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m,w\right)},\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\forall m\in \mathcal{M},\hspace{1em}0\le w\le d^{\left(m\right)}.\end{array}$

(6)

In summary, the real-time energy optimal multicommodity flow routing problem can be written as

$\begin{array}{cc}\hfill & \underset{\vec{x},\vec{s}}{\min }{\vec{c}}^{\hspace{0.17em}T}{\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{x}}^{\left(m,w\right)}\hfill \\ \hfill \text{subject to:}& \hspace{1em}{A}^{-}{\vec{x}}^{\left(m,w+1\right)}+{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}=A^{+}{\vec{x}}^{\left(m,w\right)}+{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m,w\right)},\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{0.17em}\forall m\in \mathcal{M},\hspace{1em}0\le w\le d^{\left(m\right)},\hfill \\ \hfill & \hspace{1em}\hspace{1em}D{\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\vec{x}}^{\left(m,w\right)}\le \vec{\mu },\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{s}}^{\left(m,w\right)}={\vec{s}}_{\text{out}}^{\left(m\right)}\hspace{1em}\forall m\in \mathcal{M},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{x}}^{\left(m,d^{\left(m\right)}+1\right)}=\vec{0},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{x}}^{\left(m,0\right)}=\vec{0},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{x}}^{\left(m,w\right)}\ge \vec{0},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}\ge \vec{0}.\hfill \end{array}$

(7)

Problem (7) describes the problem in a form, which can be solved in centralized way by any linear programming solver.

4.2. Intuitive Presentation of Extended Graph

In this section, we illustrate, in an intuitive way, the graph transformation, which has been discussed in Section 4.1 by mathematical equations. New variables have appeared for each communication link as well as new constraining equations for each node. These variables and constraints can be seen as virtual layers of the network where each layer represents a different communication delay w. The number of the network layers is equal to the integer deadline of the demand m plus one (the number of allowed communication hops plus a zero layer). As consistent with the structure of (6), all communication links are redirected to the nodes in the higher layer, which means that the flow is routed not only in node space but also in delay space. Because the number of layers is limited by the deadline and the flow can leave the network only in virtual nodes of the sink nodes, all possible routings through this transformed network hold the deadlines. An easy example of the graph replication is shown in Figure 1.

Figure 1

Example for the intuitive presentation.

5. Decomposition of the Routing Problem

To decompose the routing algorithm, we use the gradient optimization method to solve its dual problem. However, the linearity of the cost function of the problem (7) would cause oscillations in the gradient algorithm and prevents us finding the optimal solution. Therefore, we use the proximal-point method (for details see [4]) to modify the problem into a strictly convex form, which allows the usage of the gradient method. Moreover, we rewrite the problem into an equality form for a more transparent presentation:

$\begin{array}{cc}\hfill & \underset{{\vec{x}}^{\mathrm{\prime \prime }},{\vec{z}}^{\mathrm{\prime \prime }},{\vec{s}}^{\mathrm{\prime \prime }}}{\min }\hspace{0.17em}\underset{\vec{x},\vec{z},\vec{s}}{\min }\hspace{0.17em}g\left(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\mathrm{\prime \prime }}\right)\hfill \\ \hfill \text{subject to:}& \hspace{1em}{A}^{-}{\vec{x}}^{\left(m,w+1\right)}+{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}=A^{+}{\vec{x}}^{\left(m,w\right)}+{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m,w\right)},\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{0.17em}\forall m\in \mathcal{M},\hspace{1em}0\le w\le d^{\left(m\right)},\hfill \\ \hfill & \hspace{1em}\hspace{1em}D{\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\vec{x}}^{\left(m,w\right)}=\vec{\mu }-\vec{z},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}={\vec{s}}_{\text{out}}^{\hspace{0.17em}\left(m\right)}\hspace{1em}\forall m\in \mathcal{M},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{x}}^{\left(m,d^{\left(m\right)}+1\right)}=\vec{0},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{x}}^{\left(m,0\right)}=\vec{0},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{x}}^{\left(m,w\right)}\ge \vec{0},\hfill \\ \hfill & \hspace{1em}\hspace{1em}{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}\ge \vec{0},\hfill \end{array}$

(8)

where the objective function $\text{g}(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\hspace{0.17em}\mathrm{\prime \prime }})$ is

$\begin{array}{cc}\hfill & g\left(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\hspace{0.17em}\mathrm{\prime \prime }}\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}={\vec{c}}^T{\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\vec{x}}^{\left(m,w\right)}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}+\hspace{0.17em}\hspace{0.17em}\varepsilon {\left(\vec{z}-{\vec{z}}^{\hspace{0.17em}\mathrm{\prime \prime }}\right)}^T\left(\vec{z}-{\vec{z}}^{\hspace{0.17em}\mathrm{\prime \prime }}\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}+\varepsilon {\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\left({\vec{x}}^{\left(m,w\right)}-{\vec{x}}^{\mathrm{\prime \prime }\left(m,w\right)}\right)}^T\left({\vec{x}}^{\left(m,w\right)}-{\vec{x}}^{\mathrm{\prime \prime }\left(m,w\right)}\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}+\varepsilon {\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\left({\vec{s}}^{\hspace{0.17em}\left(m,w\right)}-{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }\left(m,w\right)}\right)}^T\left({\vec{s}}^{\hspace{0.17em}\left(m,w\right)}-{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }\left(m,w\right)}\right).\hfill \end{array}$

(9)

We have added slack variables $\vec{z}\ge \vec{0}$ into problem (8) to convert the problem into the equality form. The variables ${\vec{x}}^{\mathrm{\prime \prime }}$ , ${\vec{s}}^{\mathrm{\prime \prime }}$ , and ${\vec{z}}^{\mathrm{\prime \prime }}$ in (9) are the proximal-point variables, which have been added to convert the objective function from a linear form to a strictly convex one. Please notice that the set of optimal solutions for problem (8) is the same as for the problem (7). For the optimal solution of problem (8) it holds $\vec{x}={\vec{x}}^{\mathrm{\prime \prime }}$ , $\vec{s}={\vec{s}}^{\mathrm{\prime \prime }}$ , $\vec{z}={\vec{z}}^{\mathrm{\prime \prime }}$ . In this way the routing problem has been separated into two nested subproblems. The internal subproblem is the minimization over the variables $\vec{x}$ , $\vec{s}$ , $\vec{z}$ , and it is strictly convex. The outer subproblem minimizes the internal one by the proximal-point variables ${\vec{x}}^{\mathrm{\prime \prime }}$ , ${\vec{s}}^{\mathrm{\prime \prime }}$ , ${\vec{z}}^{\mathrm{\prime \prime }}$ .

5.1. Dual Problem

To solve the internal subproblem of (8) (minimization over the variables $\vec{x}$ , $\vec{z}$ , $\vec{s}$ ) we present its dual problem, which allows us to derive the distributable gradient algorithm. According to Slater's conditions (see e.g., [24]) the optimal solutions of the dual and primal problems have the same optimal values in this case.

The Lagrangian function of the problem (8) is

$\begin{array}{cc}\hfill & L\left(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\mathrm{\prime \prime }},\vec{\theta },\vec{\lambda },\vec{\gamma }\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}=g\left(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\mathrm{\prime \prime }}\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{\theta }}^{{\left(m,w\right)}^T}(A^{-}{\vec{x}}^{\left(m,w+1\right)}-A^{+}{\vec{x}}^{\left(m,w\right)}+{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}-{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m,w\right)})\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}+\hspace{0.17em}\hspace{0.17em}{\vec{\lambda }}^T\left(D{\displaystyle \sum }_{m\in \mathcal{M}}\hspace{0.17em}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\vec{x}}^{\left(m,w\right)}+\vec{z}-\vec{\mu }\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}+\hspace{0.17em}\hspace{0.17em}{\vec{\gamma }}^{{\left(m\right)}^T}\left({\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{s}}^{\hspace{0.17em}\left(m,w\right)}-{\vec{s}}_{\text{out}}^{\hspace{0.17em}\left(m\right)}\right),\hfill \end{array}$

(10)

where $\vec{x}\ge \vec{0}$ , $\vec{s}\ge \vec{0}$ , and $\vec{z}\ge \vec{0}$ are primal variables and $\vec{\theta }$ , $\vec{\lambda }$ and $\vec{\gamma }$ are dual variables. The dual function W is

$\begin{array}{cc}\hfill & W\left({\vec{x}}^{\mathrm{\prime \prime }},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},{\vec{z}}^{\mathrm{\prime \prime }},\vec{\theta },\vec{\lambda },\vec{\gamma }\right)=\underset{\vec{x},\vec{s},\vec{z}\ge \vec{0}}{\min }\hspace{0.17em}L\left(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\mathrm{\prime \prime }},\vec{\theta },\vec{\lambda },\vec{\gamma }\right).\hfill \end{array}$

(11)

The minimizers of the dual function (11) are

$\begin{array}{cc}\hfill {\vec{x}}_{\min }^{\left(m,w\right)}& =[{\vec{x}}^{\mathrm{\prime \prime }\left(m,w\right)}\hfill \\ \hfill & \hspace{0.17em}\hspace{1em}{\hspace{0.17em}\hspace{0.17em}-\frac{1}{2\varepsilon }\left(A^{{-}^T}{\vec{\theta }}^{\left(m,w-1\right)}-A^{{+}^T}{\vec{\theta }}^{\left(m,w\right)}+D^T\vec{\lambda }+\vec{c}\right)]}^{+},\hfill \\ \hfill {\vec{z}}_{\min }& ={\left[{\vec{z}}^{\mathrm{\prime \prime }}-\frac{1}{2\varepsilon }\vec{\lambda }\right]}^{+},\hfill \\ \hfill {\vec{s}}_{\min }^{\hspace{0.17em}\left(m,w\right)}& ={\left[{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }\left(m,w\right)}-\frac{1}{2\varepsilon }{\vec{\gamma }}^{\left(m\right)}\right]}^{+},\hfill \end{array}$

(12)

where symbol $[\cdots {]}^{+}$ denotes a positive or zero value in each component of vector $[\cdots {]}^{+}=\max \hspace{0.17em}(\vec{0},..)$ and ${\vec{x}}^{(m,0)}={\vec{x}}^{(m,d^{(m)}+1)}=\vec{0}$ .

The dual problem of the internal subproblem of (8) is

$\begin{array}{c}U\left({\vec{x}}^{\mathrm{\prime \prime }},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},{\vec{z}}^{\mathrm{\prime \prime }}\right)=\underset{\vec{\theta },\vec{\lambda },\vec{\gamma }}{\max }W\left({\vec{x}}^{\mathrm{\prime \prime }},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},{\vec{z}}^{\mathrm{\prime \prime }},\vec{\theta },\vec{\lambda },\vec{\gamma }\right).\end{array}$

(13)

The gradients of the dual function (11) are

$\begin{array}{cc}\hfill {\nabla }_{\theta }{W}^{\left(m,w\right)}& =A^{-}{\vec{x}}_{\min }^{\left(m,w+1\right)}-A^{+}{\vec{x}}_{\min }^{\left(m,w\right)}+{\vec{s}}_{\min }^{\hspace{0.17em}\left(m,w\right)}-{\vec{s}}_{\text{in}}^{\hspace{0.17em}\left(m,w\right)},\hfill \\ \hfill {\nabla }_{\lambda }{W}^{\left(m\right)}& =D{\displaystyle \sum }_{m\in \mathcal{M}}{\displaystyle \sum }_{w=1}^{d^{\left(m\right)}}{\vec{x}}_{\min }^{\left(m,w\right)}+{\vec{z}}_{\min }-\vec{\mu },\hfill \\ \hfill {\nabla }_{\gamma }W& ={\displaystyle \sum }_{w=0}^{d^{\left(m\right)}}{\vec{s}}_{\min }^{\hspace{0.17em}\left(m,w\right)}-{\vec{s}}_{\text{out}}^{\hspace{0.17em}\left(m\right)}.\hfill \end{array}$

(14)

The gradients of the dual problem (13) are

$\begin{array}{cc}\hfill {\nabla }_{x^{\prime \prime }}{U}^{\left(m,w\right)}& =-\hspace{0.17em}\hspace{0.17em}2\varepsilon \left({\vec{x}}_{\min }^{\left(m,w\right)}-{\vec{x}}^{\mathrm{\prime \prime }\left(m,w\right)}\right),\hfill \\ \hfill {\nabla }_{{\vec{z}}^{\mathrm{\prime \prime }}}U& =-\hspace{0.17em}\hspace{0.17em}2\varepsilon \left({\vec{z}}_{\min }-{\vec{z}}^{\mathrm{\prime \prime }}\right),\hfill \\ \hfill {\nabla }_{{\vec{s}}^{\mathrm{\prime \prime }}}{U}^{\left(m,w\right)}& =-\hspace{0.17em}\hspace{0.17em}2\varepsilon \left({\vec{s}}_{\min }^{\hspace{0.17em}\left(m,w\right)}-{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }\left(m,w\right)}\right).\hfill \end{array}$

(15)

5.2. Dual Gradient Algorithm

Using the dual problem (13) and the dual function (11) we rewrite the routing problem (8) into the form:

$\begin{array}{c}\underset{{\vec{x}}^{\mathrm{\prime \prime }},{\vec{z}}^{\mathrm{\prime \prime }},{\vec{s}}^{\mathrm{\prime \prime }}}{\min }\hspace{0.17em}\underset{\vec{\theta },\vec{\lambda },\vec{\gamma }}{\max }\hspace{0.17em}\underset{\vec{x},\vec{s},\vec{z}\ge \vec{0}}{\min }\hspace{0.17em}L\left(\vec{x},{\vec{x}}^{\mathrm{\prime \prime }},\vec{s},{\vec{s}}^{\hspace{0.17em}\mathrm{\prime \prime }},\vec{z},{\vec{z}}^{\mathrm{\prime \prime }},\vec{\theta },\vec{\lambda },\vec{\gamma }\right).\end{array}$

(16)

A gradient algorithm created from problem (16) consists of 2 nested loops, and it is described as (17). The internal loop solves the dual problem (13) using the gradients (14). The outer loop minimizes problem (16) over the proximal-point variables using the gradients (15):

$\begin{array}{cc}\hfill & \text{LOOP}\hspace{0.17em}\hspace{0.17em}1\hfill \\ \hfill & \hspace{1em}\text{LOOP}\hspace{0.17em}\hspace{0.17em}2\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\text{Compute}\hspace{0.17em}\hspace{0.17em}\text{the}\hspace{0.17em}\hspace{0.17em}{\vec{x}}_{\min },\hspace{0.17em}\hspace{0.17em}{\vec{z}}_{\min },\hspace{0.17em}\hspace{0.17em}{\vec{s}}_{\min }^{\hspace{0.17em}\left(m,w\right)}\hspace{0.17em}\hspace{0.17em}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\text{according}\hspace{0.17em}\hspace{0.17em}\text{to}\hspace{0.17em}\hspace{0.17em}\text{Equations}\hspace{0.17em}\hspace{0.17em}\left(13\right),\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\begin{array}{c}{\vec{\theta }}^{\left(m,w\right)}={\vec{\theta }}^{\left(m,w\right)}+\alpha {\nabla }_{\theta }{W}^{\left(m,w\right)}\\ {\vec{\gamma }}^{\left(m\right)}={\vec{\gamma }}^{\left(m\right)}+\alpha {\nabla }_{\gamma }{W}^{\left(m\right)}\\ \vec{\lambda }=\vec{\lambda }+\alpha {\nabla }_{\lambda }W\end{array}\hfill \\ \hfill & \hspace{1em}\text{END}\hspace{0.17em}\hspace{0.17em}2\hspace{0.17em}\hspace{0.17em}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\begin{array}{ccc}\hspace{1em}\hspace{1em}\hspace{1em}{\vec{x}}^{\mathrm{\prime \prime }\left(m,w\right)}={\vec{x}}^{\mathrm{\prime \prime }\left(m,w\right)}+\alpha {\nabla }_{x^{\prime \prime }}{U}^{\left(m,w\right)}& & \\ \hspace{1em}\hspace{1em}\hspace{1em}{\vec{z}}^{\mathrm{\prime \prime }}={\vec{z}}^{\mathrm{\prime \prime }}+\alpha {\nabla }_{{\vec{z}}^{\mathrm{\prime \prime }}}U& & \\ \hspace{1em}\hspace{1em}\hspace{1em}{\vec{s}}^{\prime \prime \left(m,w\right)}={\vec{s}}^{\prime \prime \left(m,w\right)}+\alpha {\nabla }_{{\vec{s}}^{\prime \prime }}{U}^{\left(m,w\right)}& & \end{array}\hfill \\ \hfill & \text{END}\hspace{0.17em}\hspace{0.17em}1\hspace{0.17em}\hspace{0.17em}\hfill \end{array}$

(17)

The $\alpha >0$ is a constant step size of the algorithm.

The important property of algorithm (17) is that it can be distributed as an in-network algorithm. Each node is responsible for the computation of the flow volume of the links leaving the node and for computation of the other corresponding variables. All the components of the variable vectors are a function of the node local variables and of the variables of the neighboring nodes that are within one hop communication distance. Only peer-to-peer communication during the algorithm is needed.

However, the distributed version of such an algorithm would have problems with the termination of the internal loop and with the synchronization of the loops between the nodes. Moreover the nested loops would cause an increase in the iterations.

To overcome these problems, we join both loops of the gradient algorithm into one loop, where we update all the variables simultaneously. Using (12), (14), and (15) we derive the iterative algorithm which is presented in Table 1, where k denotes the iteration number.

Table 1

Distributed routing algorithm.

(1) Initialize variables ${\vec{x}}_0^{\mathrm{\prime \prime }(m,w)},{\vec{s}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)},{\vec{z}}_0^{\mathrm{\prime \prime }(m,w)},{\vec{\theta }}_0^{(m,w)},{\vec{\lambda }}_0,{\vec{\gamma }}_0^{(m)}$ .

(2) Compute primal variables ${\vec{x}}_k$ , ${\vec{z}}_k$ , ${\vec{s}}_k$ according to:

${\vec{x}}_k^{(m,w)}=[{\vec{x}}_k^{\mathrm{\prime \prime }(m,w)}-(1/2\varepsilon )(A^{{-}^T}{\vec{\theta }}_k^{(m,w-1)}-A^{{+}^T}{\vec{\theta }}_k^{(m,w)}+D^T{\vec{\lambda }}_k+\vec{c}){]}^{+}$

${\vec{z}}_k=[{\vec{z}}_k^{\mathrm{\prime \prime }}-(1/2\varepsilon ){\vec{\lambda }}_k{]}^{+}$

${\vec{s}}_k^{\hspace{0.17em}(m,w)}=[{\vec{s}}_k^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}-(1/2\varepsilon ){\vec{\gamma }}_k^{(m)}{]}^{+}$

${\vec{x}}_k^{(m,0)}={\vec{x}}_k^{(m,d^{(m)}+1)}\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\vec{0}$ .

(3) Send/Receive the primal variables ${\vec{x}}_k,{\vec{z}}_k,{\vec{s}}_k$ to/from neighboring nodes.

(4) Compute dual variables ${\vec{\theta }}_{k+1},{\vec{\lambda }}_{k+1},{\vec{\gamma }}_{k+1}$ :

${\vec{\theta }}_{k+1}^{(m,w)}={\vec{\theta }}_k^{(m,w)}+\alpha (A^{-}{\vec{x}}_k^{(m,w+1)}-A^{+}{\vec{x}}_k^{(m,w)}+{\vec{s}}_k^{\hspace{0.17em}(m,w)}-{\vec{s}}_{\text{in}}^{\hspace{0.17em}(m,w)})$

${\vec{\lambda }}_{k+1}={\vec{\lambda }}_k+\alpha (D{\sum }_{m\in ℳ}{\sum }_{w=1}^{d^{(m)}}{\vec{x}}_k^{(m,w)}+{\vec{z}}_k-\vec{\mu })$

${\vec{\gamma }}_{k+1}^{(m)}={\vec{\gamma }}_k^{(m)}+\alpha ({\sum }_{w=0}^{d^{(m)}}{\vec{s}}_k^{\hspace{0.17em}(m,w)}-{\vec{s}}_{\text{out}}^{\hspace{0.17em}(m)})$ .

(5) Compute proximal-point variables ${\vec{x}}_{k+1}^{\mathrm{\prime \prime }},{\vec{s}}_{k+1}^{\hspace{0.17em}\mathrm{\prime \prime }},{\vec{z}}_{k+1}^{\mathrm{\prime \prime }}$

${\vec{x}}_{k+1}^{\mathrm{\prime \prime }(m,w)}={\vec{x}}_k^{\mathrm{\prime \prime }(m,w)}+\alpha (-2\varepsilon ({\vec{x}}_k^{(m,w)}-{\vec{x}}_k^{\mathrm{\prime \prime }(m,w)}))$

${\vec{z}}_{k+1}^{\mathrm{\prime \prime }}={\vec{z}}_k^{\mathrm{\prime \prime }}+\alpha (-2\varepsilon ({\vec{z}}_k-{\vec{z}}_k^{\mathrm{\prime \prime }}))$

${\vec{s}}_{k+1}^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}={\vec{s}}_k^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}+\alpha (-2\varepsilon ({\vec{s}}_k^{\hspace{0.17em}(m,w)}-{\vec{s}}_k^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}))$

(6) Send/Receive the dual variables ${\vec{\theta }}_{k+1},{\vec{\lambda }}_{k+1},{\vec{\gamma }}_{k+1}$ to/from the neighboring nodes.

(7) Set $k=k+1$ and start new iteration in step 2.

The correctness of such an algorithm is not seen directly from its derivation and has to be proven. The proof of the algorithm convergence is not a trivial problem, and its simplified version is presented in the Appendix. A necessary condition for the algorithm convergence assumed in the proof is $\alpha <1/2\varepsilon$ . Moreover, we have performed several simulation experiments to test the algorithm convergence in Section 6.

The initial variables ${\vec{x}}_0^{\mathrm{\prime \prime }(m,w)}$ , ${\vec{s}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}$ , ${\vec{z}}_0^{\mathrm{\prime \prime }(m,w)}$ , ${\vec{\theta }}_0^{(m,w)}$ , ${\vec{\lambda }}_0$ , ${\vec{\gamma }}_0^{(m)}$ are set to arbitrary initial values. The closer the values are to the final solution, the faster the algorithm converges. This property can be used in the case of minor changes of the network structure during its operation or in case of a precomputed routing, for example, based on Dijkstra's algorithm.

5.3. Variables Initialization

As we have mentioned in Section 5.2, the number of iterations depends on the initial variables ${\vec{x}}_0^{\mathrm{\prime \prime }(m,w)}$ , ${\vec{s}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}$ , ${\vec{z}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}$ , ${\vec{\theta }}_0^{(m,w)}$ , ${\vec{\lambda }}_0$ , ${\vec{\gamma }}_0^{(m)}$ . The closer the initial variables are to the final solution the less number of iterations are needed. In this work, we focus only on the dual variables. The proximal-point variables ${\vec{x}}_0^{\mathrm{\prime \prime }(m,w)}$ , ${\vec{s}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}$ , ${\vec{z}}_0^{\mathrm{\prime \prime }(m,w)}$ are set to zero. The variables ${\vec{\lambda }}_0$ , ${\vec{\gamma }}_0^{(m)}$ cannot be estimated without the knowledge of the final routing. We set them to zero.

Let us suppose some initial variable ${\vec{\theta }}_0$ which satisfies a condition:

$\begin{array}{c}{A}^{{-}^T}{\vec{\theta }}_0^{(m,w-1)}-A^{{+}^T}{\vec{\theta }}_0^{\left(m,w\right)}+D^T{\vec{\lambda }}_0+\vec{c}\ge \vec{0}.\end{array}$

(18)

According to (12), such initial variables do not cause any changes of the algorithm variables on its own (i.e., no variable changes if ${\vec{s}}_{\text{out}}^{\hspace{0.17em}(m)}={\vec{s}}_{\text{in}}^{\hspace{0.17em}(m)}=0$ ).

According to the complementary slackness condition (for details see, e.g., [4]) for the optimal solution holds for two neighboring nodes: if $({\theta }_{k,l^{+}}^{(m,w-1)}-{\theta }_{k,l^{-}}^{(m,w)})<c_l$ then $x_{k,l}^{(m,w)}=0$ , if $({\theta }_{k,l^{+}}^{(m,w-1)}-{\theta }_{k,l^{-}}^{(m,w)})=c_l$ then $x_{k,l}^{(m,w)}\ge 0$ , and the link l is a part of the shortest path for the demand $m\in ℳ$ . (We use the index l to denote the vector component corresponding to the link l, $l^{+}$ to denote end node of the link, and $l^{-}$ to denote the start node of the link.) This fact leads us directly to Dijkstra's algorithm which can be used in a distributed way. If for all sink nodes $n_{\text{out}}$ we set ${\theta }_{0,n_{\text{out}}}^{(m,w)}=0$ for all $m\in ℳ,\hspace{0.17em}\hspace{0.17em}0\le w\le d^{(m)}$ and for the other nodes n we set ${\theta }_{0,n}^{(m,w)}=the\hspace{0.17em}\hspace{0.17em}shortest\hspace{0.17em}\hspace{0.17em}distance\hspace{0.17em}\hspace{0.17em}to\hspace{0.17em}\hspace{0.17em}the\hspace{0.17em}\hspace{0.17em}sink\hspace{0.17em}\hspace{0.17em}node$ , we get an initial setting, which satisfies condition (18). Moreover, this initial setting is much closer to the optimal solution than the setting ${\theta }_{0,n}^{(m,w)}=0$ for all $n,m,w$ . In Section 6.3, we present several experiments to evaluate the heuristic behavior in comparison with the zero initial setting. Please notice that due to the capacity constraints, the heuristic solution does not need to be equal to the final optimal solution. According to our experiments it rapidly decreases the number of iterations.

6. Experiments

To demonstrate the behavior and the correctness of the distributed routing algorithm, we have performed several experiments in Matlab. We have focused on a basic problem, where for each communication demand one node sends data flow to one sink node (i.e., a multicommodity mono-source, mono-sink problem). Therefore, $s_{\text{in},n_1}^{(m)}=1$ for the source node $n_1$ and $s_{\text{out},n_2}^{(m)}=1$ for the sink node $n_2$ of the communication demand m.

The random networks for the experiments have been constructed as follows. We consider a square field of size [size × size], where the size is changing during the experiments. The field is divided into subsquares of size $[1\times 1]$ . One node is randomly placed into each subsquare, and the communication distance is set to 1.7 (i.e., node $n_1$ can communicate with node $n_2$ , if and only if their Euclidean distance is less than 1.7). Please notice that our network is close to the “unit-disk network” [2]. The communication costs $\vec{c}$ per transmitted data flow unit have been set as the square power of the distance between the nodes. The link capacities have been set to two ${\mu }_n=2$ , and the maximum number of communication hops to $d^{(m)}=6$ for all $m\in ℳ$ . The constants of the algorithm have been set as $\alpha =0.03$ and $\epsilon =0.3$ . The initial values ${\vec{x}}_0^{\mathrm{\prime \prime }(m,w)}$ , ${\vec{s}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)}$ , ${\vec{\theta }}_0^{(m,w)}$ , ${\vec{\lambda }}_0$ , ${\vec{\gamma }}_0^{(m)}$ have been set to 0 and $z_{0,n}^{\mathrm{\prime \prime }}={\mu }_n$ for all experiments except for those in Sections 6.2 and 6.4. Only feasible problems are used.

During the experiments we evaluate k as a number of iterations needed to achieve less than 1% deviation of the cost function from the optimal value, less than 1% capacity violation, and less than 1% flow conservation law violation during the last 500 iterations. The 500 iterations are included in the statistics. (the optimal value was computed separately by a centralized algorithm for evaluation purposes only).

Unfortunately, as the presented approach of dual decomposition of the routing problems in node-link form is a new technique, there are only few works in this area at this moment. The most similar work can be found in [21]. However, it is focused on slightly different problems. It uses different energy consumption and communication delay models, and it is limited only on problems with strictly convex objective functions. Our work is focused on the problems with linear objective functions and presents different qualities than [21]. Due to the differences in the used models and algorithms derivations the experimental comparison would have a disputable contribution and would be strongly dependent on the chosen problems.

6.1. Example

To present the resulting optimal data flow routing in the network we have performed an experiment based on the network described above. The field size has been set to 9 (i.e., 81 nodes in the network), and the number of communication demands has been set to 8. The link capacities have been set to ${\mu }_n=4$ for this problem. The initial variables ${\vec{\theta }}_0$ have been set according to Section 5.3.

The optimal data flow routing in the network is shown in Figure 2. The 3D routing model in the node-delay space is presented in Figure 3 where the vertical axis represents the communication delay. No communication demand routing has more than 6 communication hops.

Figure 2

Optimal data flow routing (multicommodity, mono-source, mono-sink problem).

Figure 3

Optimal data flow routing in node-delay space.

The algorithm behavior can be seen in [25] on video. The flow routing in the network and in the node-delay space and the progress of the flow conservation law are presented.

6.2. Algorithm Convergence

We have performed a set of experiments with random initial variables on random networks, and we have evaluated the algorithm convergence. The field size has been set to 7. There are 8 communication demands in the network. The initial values have been set randomly from intervals: ${\vec{x}}_0^{\mathrm{\prime \prime }(m,w)},{\vec{s}}_0^{\hspace{0.17em}\mathrm{\prime \prime }(m,w)},{\vec{z}}_0^{\mathrm{\prime \prime }(m,w)}\in 〈\mathrm{0,2}〉$ for proximal-point variables, and ${\vec{\theta }}_0^{(m,w)},{\vec{\lambda }}_0,{\vec{\gamma }}_0^{(m)}\in 〈\mathrm{0,20}〉$ for dual variables. The intervals have been chosen as the double value of typical optimal values.

The algorithm has been run 2000 times on random networks. The results are presented in Figure 4. There, the number of iterations is placed on the horizontal axis, and the number of experiments, which has been finished before the number of iterations, is on the vertical axis.

Figure 4

Algorithm convergence with random initial variables.

This experiment provides an important practical verification of the theoretical proof of the algorithm convergence. 95% of the experiments have been finished before 14300 iterations.

6.3. Number of Iterations

To demonstrate the statistical behavior of the algorithm, we have gradually increased the number of nodes from 16 to 100 (i.e., the field size from 4 to 10) for 8 communication demands. The computation has been repeated, on random networks 1000 times for each number of nodes.

The results have been evaluated as a maximum, average, and minimum number of iterations needed to achieve the optimal value, and it is presented in Figure 5. There, the experiment progress is presented for the basic algorithm without the initial heuristic in black (the initial variables have been set to zero) and for the algorithm with the initial heuristic according to Section 5.3 in gray.

Figure 5

The number of iterations in relation to the number of nodes. In black, for zero initial variables. In gray, for initialization according to Section 5.3.

Similar experiments have been performed for the iterations dependence on the number of communication demands. We have gradually increased the number of communication demands from 1 to 16 in networks with 49 nodes. The computation has been repeated, on random networks 1000 times for each number of communication demands. The results are presented in Figure 6 for both algorithms with and without the initial heuristic.

Figure 6

The number of iterations in relation to the number of communication demands. In black, for zero initial variables. In gray, for initialization according to Section 5.3.

In Table 2, we present the percentages of infeasible problems for each number of communication demands. The values say how much percent of the generated problems have not been feasible due to the capacity and real-time constraints. From Table 2, it follows that the network has been close to saturation and that the constraints affect the final routings.

Table 2

Percentage of infeasible problems in the experiment for the number of communication demands.

Number of comm. demands	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
Infeasible problems [%]	0.5	2.0	3.3	4.0	3.9	5.8	7.0	8.6	12.0	13.8	17.5	20.7	27.1	33.1	39.3	48.1

The important outcome of these experiments is the observation that the number of iterations is approximately linear. It follows that the algorithm is easily applied to networks with many nodes. We can see a significant improvement of the initial heuristic in the Figures.

6.4. Network Change

The advantage of the one-loop algorithm presented in this work is that it can automatically adjust the data routing in case of network changes. In order to evaluate the algorithm behavior in this case, we have simulated a dying node in the network as follow.

(1) We generated a random network with communication demands and found the optimal solution. (2) In the original problem from step 1, we removed one node and measured the number of iterations to find the new optimal solution. (3) We repeated step 2 for 15 different nodes with the largest data flow.

The simulation has been performed for 942 random original networks with field size set to 7 and with 8 communication demands (i.e., $942\hspace{0.17em}\hspace{0.17em}\times \hspace{0.17em}\hspace{0.17em}15=14130\hspace{0.17em}\hspace{0.17em}\text{experiments}$ ).

The results are presented in Figure 7. There, the number of iterations is placed on the horizontal axis, and the number of experiments, which has been finished before the number of iterations, is on the vertical axis.

Figure 7

Algorithm convergence for network change simulation.

$95\text{\%}$ of the experiments have been finished before 4466 iterations. Moreover, less than $15.0\text{\%}$ of the dual variables and $2.0\text{\%}$ of the primal variables have been changed during the experiments, which significantly decrease the amount of transmitted data.

Video presentation can be seen in [25].

7. Conclusion

We have presented a distributed algorithm for the real-time energy optimal routing in ad hoc and sensor networks for systems with linear cost functions and constant communication delays. Such networks are used very often in the industrial environment, where TDMA-derived mechanisms or mechanisms with a significant stack computation delay are used. We have described the routing problem as a multicommodity network flow optimization problem, modified it into a time aware form, and used the dual decomposition method to derive the distributed algorithm. The algorithm does not need any central computational node with knowledge about the whole network structure, and it only uses local communication between the neighboring nodes. This rapidly increases the robustness of the algorithm in the case of partial network damage. We have performed several simulations to evaluate the algorithm behavior and to test its convergence. The mathematical proof of the algorithm convergence is available in the Appendix.

The main purpose of this paper was to present the basic concept of a new in-network distributed routing algorithms for real-time data and its derivation. From the experimental section it is seen that the algorithm is not application-ready because of the high number of iterations, which leads to high number of communications. However, the main strength of the algorithm is not to find the whole optimal routing in an unknown network, but to adapt an existing routing in case of local network changes (dead/new node, loss of connection, communication costs, or capacities change, etc.) where the number of data communications could be significantly decreased. Moreover, the algorithm can easily adapt to slow network changes. The algorithm is suitable for static networks with sporadic changes or for networks with slow continuous changes, so the routing can be adjusted.

According to our preliminary experiments the number of iterations can be significantly decreased in future work. The algorithm can be extended by heuristics based on the partial knowledge about the network structure (e.g., node geographical position) and heuristics based on Newton's method. The results of our preliminary testing indicate that Newton's method-based heuristics can decrease the number of iterations more than 3 times.

Considering the fact that the algorithm is based on Linear programming formulation, we believe that the principle of the algorithm and the approach used to its derivation can be used to solve many different problems in the sensor networks area, like resource sharing, network localization, object tracking, and so forth.