Non-Concave Network Utility Maximization in Connectionless Networks: A Fully Distributed Traffic Allocation Algorithm
Jingyao Wang Mahmoud Ashour Constantino Lagoa Necdet Aybat Hao Che Zhisheng Duan


theoremTheorem
lemmaLemma
propositionProposition
corollaryCorollary
definitionDefinition
remarkRemark

#1⟨#1⟩












 


 [address_wang]Jingyao Wang and Zhisheng Duan are with the Department of Mechanics and Engineering Science, Peking University, Beijing 100871, China, (yayale.8@163.com, duanzs@pku.edu.cn).
 

[address_mahmound]Mahmoud Ashour and Constantino Lagoa are with the Department of Electrical Engineering and Computer Science, Pennsylvania State University, University Park, PA 16802, USA, (mma240@psu.edu, lagoa@engr.psu.edu).


[address_necdet]Necdet Aybat is with the Department of Industrial Engineering, Pennsylvania State University,
      University Park, PA 16802, USA, (nsa10@psu.edu).
        
[address_hao]Hao Che is with the Department of Computer Science and Engineering, University of Texas at Arlington, Arlington, TX 769019, USA, (hche@cse.uta.edu).     
 
 
[acknowledgement] This work was partially supported by NSF grants CNS-1329422, CMMI-1635106,
FCC-1629625, NNSF of China grants 61673026 and the China Scholarship Council. 
 





This paper considers the optimization-based traffic allocation problem among multiple end points in connectionless networks. The network utility function is modeled as a non-concave function, since it is the best description of the quality of service perceived by users with inelastic applications, such as video and audio streaming. However, the resulting non-convex optimization problem, is challenging and requires new analysis and solution techniques. To overcome these challenges, we first propose a hierarchy of problems whose optimal value converges to the optimal value of the non-convex optimization problem as the number of moments tends to infinity.
From this hierarchy of problems, we obtain a convex relaxation of the original non-convex optimization problem by considering truncated moment sequences. For solving the convex relaxation, we propose a fully distributed iterative algorithm, which enables each node to adjust its date allocation/ rate adaption among any given set of next hops solely based on information from the neighboring nodes. Moreover, the proposed traffic allocation algorithm converges to the optimal value of the convex relaxation at a O(1/K) rate, where K is the iteration counter, with a bounded optimality. At the end of this paper, we perform numerical simulations to demonstrate the soundness of the developed algorithm.












§ INTRODUCTION

Applications and services supported by modern communication networks have diverse requirements, e.g., high throughput and low latency. Traffic engineering (TE) has long been used to optimize the utilization of the limited network resources so that such requirements are fulfilled. This entails developing data rate allocation algorithms and congestion control protocols capable of maximizing a given network utility subject to network resource constraints <cit.>. Many problems of recent interest arising in diverse fields can be cast as an optimization problem, and network utility maximization (NUM) is no different.

In large-scale networks, the size of the optimization problems rapidly increases as the number of nodes and links increase. This stimulates the necessity of developing decentralized control algorithms capable of decomposing the high-dimensional problem into separate moderate-size subproblems that can be solved independently and locally at various network nodes. The main idea behind such decentralized control algorithms is to distribute the computations required for the solution of the optimization problem among various nodes <cit.>-<cit.>. This approach exploits local information available at each node. Nevertheless, information exchange among different nodes is inevitable since distinct data flows share the same network resources. Therefore, distributed optimization approaches not only aim at decomposing the problem, but also minimizing the communication overhead.

In the benchmark work by Kelly et. al. <cit.>, the optimization of the utility of a large-scale broadband network with limited bandwidth resources is considered. The authors propose two classes of rate control algorithms by casting the NUM problem in both primal and dual forms. In <cit.>, a family of decentralized sending rate control laws are proposed to steer the traffic allocation to an optimal operating point while avoiding congestion. A non-linear control theoretic approach is employed in <cit.> to derive adaptation laws that enable each node to independently distribute its traffic optimally among any given set of next hops. More recently, reference <cit.> considers the NUM, derives its dual problem, and uses a distributed gradient-based approach for its solution. A similar approach appears in <cit.>. In spite of the existence of a relatively dense literature on NUM, most available results consider only the optimization of concave utility functions. However, it has been shown that the reward experienced by the users of real-time applications, such as video and audio streaming, cannot be accurately modeled using concave functions. Reference <cit.> shows that the video quality perceived by users on a mobile device is a non-decreasing and step-like function with respect to the data rate, because users have almost similar quality of experience on 3 Mbps and 1 Mbps <cit.>. This observation motivates considering the optimization of non-concave network utility functions, which constitutes a main focus of this paper.

Non-concave NUM is a non-convex optimization problem; hence, it is difficult to solve. Nevertheless, there exist some attempts in the literature for deriving algorithms that provide near-optimal solutions. Reference <cit.> develops a centralized algorithm that solves the NUM problem with polynomial utilities. Reference <cit.> determines the conditions under which the standard distributed dual-based algorithm can still converge to the global optimal solution with non-concave utilities.

This paper develops a distributed iterative algorithm for the optimization of a generalized class of non-concave network utility functions that capture a wide variety of real-world applications. In particular, we focus on connectionless networks, where each node is required to distribute its traffic among a set of next hops without prior arrangement so that the network utility is maximized. We handle the challenge posed by the non-convexity of the optimization problem by developing a sequence of convex relaxations whose solution converges to that of the original problem. We use results on polynomial optimization and moment sequences to derive the convex relaxations <cit.>. Furthermore, we propose an iterative primal-dual algorithm <cit.> that enables each node to distribute its traffic among the set of next hops. We emphasize on the distributed nature of the algorithm, where each node uses its local information and need not communicate with other nodes except its direct neighbors.





§ NOTATION


Throughout this paper, the traffic flows are assumed to be described by a fluid flow model, and the only resource constraint taken into account is link bandwidth. In the remainder of this paper, call and flow will be used interchangeably.


Let 𝒩 denote the set of nodes in the network,
and ℒ⊂𝒩×𝒩 denote the set of links connecting particular pairs of nodes. We assume that each link l ∈ℒ has a finite capacity c_l>0. Moreover, let 𝒮≜{ s_1,s_2,…,s_n} and 𝒟≜{ d_1,d_2,…,d_n} denote respectively the set of source nodes and the set of destination nodes contained in 𝒩 such that 𝒮∩𝒟= ∅. The intended destination for each source node s_i is d_i for i∈ℐ≜{1,…,n}, i.e., without loss of generality, we assume that there is a one-to-one correspondence between  𝒮 and 𝒟, and ℐ denotes the set of different flow (call) types in the network. Given source node s∈𝒮, let ℒ_s denote the set of links connected to it. Let the sending data rate through link l∈ℒ_s be x_s,l^out, and all such sending data rates be 𝐱_s^out≜ [x_s,l^out]_l∈ℒ_s. We define the aggregate sending data rate of s∈𝒮 be denoted by r_s≜∑_l∈ℒ_sx_s,l^out. Also, let ℬ≜𝒩∖  (𝒮∪𝒟)= {b_1,b_2,…,b_m} denote the set of forwarding nodes contained in 𝒩. Given b∈ℬ, let ℐ_b be the set of flows visiting node b, and ℒ_b⊆ℒ denote the set of links connected to it. Suppose ℒ_b,i^out⊆ℒ_b denote the set of outgoing links from b associated with calls (flows) of type i∈ℐ_b. Similarly, let ℒ_b,i^in⊂ℒ_b denote the set of incoming links to b associated with calls (flows) of type i ∈ℐ_b. Furthermore, given b∈ℬ, for each i∈ℐ_b and l∈ℒ_b,i^out, let x_i,b,l^out denote the data rate of call type i ∈ℐ_b, associated with  s_i and d_i, forwarded from node b through link l∈ℒ_b,i^out.
The above notation is exemplified in Fig. 1 for the case of allocating flows associated with two source nodes, s_1 and s_2, and two destination nodes, d_1 and d_2.



Given b∈ℬ and l∈ℒ_b, 
 let ℐ_b,l^in⊂ℐ be the set of call types forwarded to node b through link l, and ℐ_b,l^out⊆ℐ_b be the set of call types forwarded from node b through link l. Moreover, given node b∈ℬ and link l∈ℒ_b, let e_l(b) denote the adjacent node to b through link l.
We summarize all the notation for the communication network in Table I for the convenience of the reader.


Now, given node b∈ℬ, let the vector containing all flow rates departing from node b through link l∈ℒ_b be denoted by 𝐱_b,l^out≜ [x_i,b,l^out]_i∈ℐ_b,l^out∈R_+^|ℐ_b,l^out|, where |.| denotes the cardinality of a set.


Given node b∈ℬ and l∈ℒ_b, let 1_b,l∈R^1× |ℐ_b,l^out| be the row vector with all elements equal to 1. In a similar way, let δ_b,l∈R^1× |ℐ_b,l^in| be the row vector with all elements equal to 1 if link l is bidirectional, and 0 otherwise.




Also, let . denote the Euclidean norm. Given a convex set 𝒜, let 𝐼_𝒜(.) denote the indicator function of 𝒜, i.e., 𝐼_𝒜(ω)=0 for ω∈𝒜 and equal to + ∞ otherwise, and let 𝑃_𝒜(ω)≜argmin{υ-ω: υ∈𝒜} denote the projection onto 𝒜. Given a closed convex set 𝒜, we define the distance function as d_𝒜(ω)≜𝑃_𝒜(ω)-ω. Also, 𝐈_n is the n× n identity matrix.



















§ PROBLEM FORMULATION


Consider a communication network consisting of a set of source nodes 𝒮. Each source node s∈𝒮 has a local utility function U_s(r_s):R_+→R_+ of its sending data rate r_s. For a fixed order ℓ >0, the utility function is defined as a general non-concave polynomial-like function in the form

    U_s(r_s)≜∑_j=0^ℓ p_s,j (r_s)^j/ℓ.

This particular form of objective functions is so flexible that it can be used to approximate a wide variety of functions arising in practical applications such as step functions for the video streaming case <cit.>.

The objective of this paper is to design a data rate allocation algorithm for the communication network such that the utilization of resources is maximized, while satisfying the network resource constraints. The network resource constraints considered in this paper include link capacity constraints, Minimum Rate Guaranteed and Upper Bounded Rate Service (MRGUBRS) requirements, and flow conservation constraints through nodes.

More precisely, for any link l∈ℒ, the aggregated flows going through this link should not exceed the link capacity. For example, in Fig. 1, the bidirectional link l_3 is shared by flows belonging to two source nodes. The data rates x_1,b_2,l_3^out and x_2,b_3,l_3^out going through this link should satisfy that

    x_1,b_2,l_3^out + x_2,b_3,l_3^out≤ c_l_3.

For the unidirectional link l_2,
node b_2 forwards data rate x_2,b_2,l_2^out through this link. Then, x_2,b_2,l_2^out is upper bounded by c_l_2.



Given flows of type i∈ℐ, recall that flows of type i∈ℐ is associated with source/destination pairs s_i/d_i.   For fixed link l∈ℒ_s_i, the corresponding data rate x_i,l^out is determined at source node s_i ∈𝒮 and multiple paths are available for transporting these flows.
More precisely, each node on these paths divide incoming traffic into available links by striving to conserve the flows belonging to each source node (i.e., aims at no losses) and to avoid link congestion. In Fig. 1, node b_3 tries to satisfy

    x_1,b_2,l_3^out=x_1,b_3,l_4^out+x_1,b_3,l_7^out.


Finally, flows belonging to each source node s∈𝒮 is assumed to be of the MRGUBS category, i.e., for some 0< ξ_s<ζ_s and s∈𝒮,

    ξ_s≤ r_s≤ζ_s.



Now, considering the above constrains and assumptions, we can formulate the problem of optimal
traffic allocation as follows:

    maximize∑_s∈𝒮 U_s(r_s),

subject to the network capacity constraints [Note that the formulation in this paper allows for the existence of bidirectional links.]

    ∑_i ∈ℐ_b,l^out x_i,b,l^out+∑_i ∈ℐ_b,l^inx_i,e_l(b),l^out≤ c_l, l ∈ℒ_b, b ∈ℬ,

the flow conservation constraints at each node

    ∑_l ∈ℒ_b,i^in x_i,e_l(b),l^out- ∑_l̃∈ℒ_b,i^outx_i,b,l̃^out=0, i ∈ℐ_b,  b ∈ℬ,

the non-negativity of forwarded data rates constraints

    x_i,b,l^out≥ 0, i∈ℐ_b,l^out, l∈ℒ_b,  b∈ℬ,

and the MRGUBS requirements

    ( 𝐱_s^out, r_s)∈𝒳_s, s∈𝒮,

where the set 𝒳_s is defined as

    𝒳_s≜{  ( 𝐱_s^out, r_s)∈R_+^|ℒ_s|×R_+:
    ξ_s ≤ r_s≤ζ_s, r_s= ∑_l∈ℒ_s x_i,l^out}.







Most literature in the context of NUM considers maximizing concave diminishing functions.
However,
modern communication networks are dominated by various inelastic applications, such as internet video and audio streaming. Users' satisfaction for these applications cannot be modeled with concave functions. It is better to be described as non-concave functions. For instance, the utility for voice applications is a sigmoidal function <cit.>. Thus, we consider users' perceived qualification of Cost of Service (CoS) and model the utility function as a general class of non-concave polynomial functions.
Moreover, the challenges of attempting to solve the resulting traffic allocation problem (<ref>) are two-fold.
First, the optimization problem obviously constitutes a non-convex problem since its objective function is non-concave. Second, global information on fast timescale events, as required in the above formulation, is not generally available. The latter fact stimulates the necessity of developing a distributed algorithm that converges to the optimal data rate allocation of the non-convex NUM problem.




§ MAIN RESULTS

In this section, we present our approach used to overcome the challenges opposed by the non-convexity of the optimization problem. In particular, we first present a convex relaxation to the non-convex NUM problem (<ref>). This convex relaxation is chosen from a hierarchy of optimization problems whose optimal value converges to the optimal value of problem (<ref>) as the number of moments tends to infinity. For solving the convex relaxation problem, we propose a distributed primal-dual algorithm (DPDA), which enables all nodes to update their data rate allocation solely using immediate local information.
A salient feature of the proposed algorithm is that the iterate sequence converges to the optimal solution at a O(1/K) rate, where K is the iteration counter, with a bounded optimality. 




 §.§ NUM convex relaxation

The non-convexity of the optimization problem (<ref>) opposes challenges for us to analysis and solve the traffic allocation problem. However, the following proposition provides a hierarchy of optimization problems whose optimal value converges to the optimal value of the non-convex optimal problem (<ref>). For solving the traffic allocation problem, we choose a convex one from this hierarchy of problem by truncating the number of moments to the finite case. This proposition is one of the main results of this paper.






The solution of the following optimization problem converges to the solution
of the non-convex NUM problem (<ref>) with non-concave user utility functions of the form (<ref>) as the positive parameter α→∞. Moreover, problem (<ref>) is convex if α≤ℓ.

    𝐱maximize      ∑_s ∈𝒮𝐩_s^T𝐦_s
       subject to      m_s,0=1,  s∈𝒮,
             𝐌(0, α,𝐦_s)≽ 0,  s∈𝒮, 
             β_s𝐌(0, α-2,𝐦_s)- 𝐌(2, α,𝐦_s) ≽ 0, 
              m_s,j≤ (r_s)^j/ℓ, j∈{1,…,α}, s ∈𝒮,
              x_s,l^out≤ c_l, l ∈ℒ_s,  s∈𝒮,
             1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out≤ c_l,  l ∈ℒ_b, b ∈ℬ,
              B 𝐱=0,
              ( 𝐱_s^out, r_s)∈𝒳_s, s∈𝒮,
             𝐱_b,l^out≽ 0, l∈ℒ_b,  b∈ℬ.

The objective function is a linear function of variables 𝐦_s=[m_s,j]_j ∈{0,,α} with parameters 𝐩_s=[p_s,j]_j ∈{0,,α}.
The decision variable 𝐱 of problem (<ref>) is a vector consisting of the data rate x_s,l^out, r_s and 𝐦_s for each s∈𝒮, and the sending data rate
𝐱_b,l^out, b∈ℬ for each l∈ℒ. More precisely, the dimension of vector 𝐱 is ∑_s∈𝒮 (|ℒ_s |+α+2)+∑_b∈ℬ∑_i∈ℐ_b|ℒ_b,i^out|.
In the constraints,
B∈R^(∑_b∈ℬ |ℐ_b|)×(∑_s∈𝒮 (|ℒ_s |+α+2)+∑_b∈ℬ∑_i∈ℐ_b|ℒ_b,i^out| ) denotes the edge-node-like incidence matrix, i.e., the entry B_(s,b,l),ω, corresponding to flow-node-link triplet (s,b,l)∈𝒮×ℬ×ℒ and ω∈𝐱, equal to 1 if the data rate ω of flows belonging to source node s is forwarded from node b through link l, -1 if the the data rate ω is received at node b, and 0 otherwise.
 β_s is a known upper bound on the aggregate data rate of source s∈𝒮,
and the moment matrices 𝐌∈R^h+1×R^h+1 are of the form

    𝐌(k,k+2h,𝐦_s) =
    [[     m_s,k   m_s,k+1         …   m_s,k+h;   m_s,k+1         ⋱         ⋱ m_s,k+h+1;         ⋮         ⋱         ⋱         ⋮;   m_s,k+h         …         …  m_s,k+2h ]].





The proof is shown in Appendix A.


Hereafter, we use α = ℓ.
It is worth mentioning that the result of Proposition 1 holds for the even order ℓ. Nonetheless, similar results can be derived for the odd ℓ, which is omitted for brevity. The proposed problem (<ref>) constitutes a convex optimization problem, because it maximizes the sum of linear functions subject to convex constraints. Therefore, it can be easily solved if global information is available. Nevertheless, the objective of this paper is to solve this problem in a distributed fashion that leverages per hop information available at each node.

Before moving on, we
introduce some notation that renders the formulation of (<ref>) conveniently compact.
For every s ∈𝒮, let the set 𝒜_s be defined as

    𝒜_s={  (𝐱_s^out, 𝐦_s,  r_s)∈R_+^|ℒ_s|×R^ℓ+1×R_+: m_s,0=1, 
       𝐌(0, ℓ, 𝐦_s)
    ≽ 0,  β_s𝐌(0, ℓ-2, 𝐦_s)-  𝐌(2,ℓ, 𝐦_s) ≽ 0,
        x_s,l^out≤ c_l, l∈ℒ_s, m_s, j≤ (r_s)^j/ℓ, j∈{1,…,ℓ}, ( 𝐱_s^out, r_s)∈𝒳_s}.





 §.§ Algorithm DPDA

The constrains set of convex relaxation (<ref>) consists of local constraints, e.g., capacity constraints and global constraints, e.g., flow conservation constraints through nodes. The existence of global constraints renders difficulty for us to solve problem (<ref>) in a distributed fashion. However, 
the primal-dual method, proposed by Chambolle and Pock in <cit.> for solving convex-concave saddle point problems makes it possible. This algorithm can be adapted to solve the multi-agent consensus optimization problem as discussed in <cit.>. We also use the distributed primal-dual algorithm in <cit.> to solve our traffic allocation problem (<ref>). We present the resulting iterative algorithm, i.e., DPDA, of which iterate sequence converges to the solution of (<ref>). The details of developing DPDA can be found in Appendix B.





The suboptimality and feasibility of the DPDA
iterate sequence can be bounded as in the following theorem.




Given the communication network and the convex optimization problem (<ref>). Let d_s>0, s∈𝒮 and d_i,b,l>0, i∈ℐ_b,l^out, l∈ℒ_b, b∈ℬ be given (sufficiently large) constants. Recall that the decision variable 𝐱 of problem (<ref>) is a vector consisting of the data rate 𝐱_s,l^out, r_s and 𝐦_s for each s∈𝒮, and the sending data rate
𝐱_b,l^out, b∈ℬ for each l∈ℒ. Also recall that vector variables λ,θ are the dual variables associated with the capacity constraints and the flow conservation constraints at nodes, respectively. Let (𝐱^⋆,λ^⋆,θ^⋆) be an arbitrary saddle-point for the Lagrange function of problem (<ref>), and
{𝐱^k}_k≥ 0 be the iterate sequence generated using Algorithm DPDA, initialized from an arbitrary 𝐱^0 and [λ_b,l^0]_l∈ℒ_b,b∈ℬ=0. 
 Let the primal-dual step sizes [τ_s]_s∈𝒮,[τ_i,b,l]_i∈ℐ_b,l^out,l∈ℒ_b,b∈ℬ and γ be positive constants satisfying the following inequalities
 
    1/τ_s-γ(4+d_s ) ≥ 0,

for all s∈𝒮, and
 
    1/κ_b, l(1/τ_i,b,l-γ(4+d_i,b,l ) )≥ m_l+1,

for all i∈ℐ_b,l^out,l∈ℒ_b,b∈𝒩, where m_l is the total number of sources using link l to transport flows. Denote the average of sending data rates
by 𝐱̅^K≜1/K∑_k=1^K 𝐱^k, where K≥ 1. Then, {𝐱̅^K} converges to the maximum of the utility function of the problem (<ref>) subject to the resource allocation constraints. In particular, the average of the iterative sequence asymptotically converges to the feasible solution, i.e.,

    θ^⋆B𝐱̅^K+  ∑_b∈ℬ∑_l∈ℒ_bλ^⋆_b,l h(𝐱̅_b,l^out,𝐱̅_e_l(b),l^out)
    ≤Θ_1/K, ∀ K≥ 1.

It also asymptotically maximizes the utility function of the problem (<ref>), i.e.,

    | ∑_s∈𝒮𝐩_s^T(𝐦̅_s-𝐦_s^⋆)|
    ≤Θ_1/K, ∀ K≥ 1,

where the notation h(𝐱̅_b,l^out,𝐱̅_e_l(b),l^out) and Θ_1 is defined in Appendix C.






The proof is presented in Appendix C.



Algorithm DPDA is a fully distributed traffic allocation algorithm. 
This point can be verified by looking through the implementation procedure. 
The step-size parameters are decided before implementing the algorithm. It is given in Theorem <ref> that those parameters  satisfy conditions (<ref>) and (<ref>), both of which are local conditions. Thus, choosing the parameters requires no global information. In the first step, the variables z_s,l, l∈ℒ_s, s∈𝒮 and z_i,b,l, i ∈ℐ_b,l^out, l ∈ℒ_b, b∈ℬ are local variables respectively introduced for each source node and each forwarding node. It is worth noting that giving the initial state value of x_s,l^out, l∈ℒ_s, s∈𝒮 and x_i,b,l^out, i ∈ℐ_b,l^out, l ∈ℒ_b, b∈ℬ to those introduced variables is also a local operation. For the first iteration, i.e., K=1, in steps 3 and 4, DPDA enables all nodes to update their sending data rates in parallel. Each node solely uses immediate information from its neighboring nodes to perform all computations. In step 5, the link price λ_b, l^k+1, l ∈ℒ_b,  b∈ℬ is updated with new local data rate allocation solution. This step can be performed at both end points that each link connects, which just uses their local information. Step 6
updates the introduced local variables with the new local data rate allocation solution. The iterative procedure continues until the iterate sequence converges to the optimal solution.




It follows from inequalities (<ref>) and (<ref>) that DPDA converges at the rate of 𝑂(1/K), where K is the number of iterations.





If the problem (<ref>) has a unique solution, then the sequence of sample averages converges to that solution.








§ SIMULATION RESULTS

In this section, we present some simulation results which exemplify the behavior of the proposed algorithm, i.e., Algorithm DPDA. The simulations show that the final data rate allocation results in a value of the utility function barely distinguishable from the optimal one.

We consider the network model shown in Fig. 2, where we also show all the links' bandwidths, and source-destination pairs. The network model allows for multiple paths available for flows belonging to each source node. We consider a total of 8 different combinations of source/destination
nodes. Moreover, we list the prescribed next hops for all forwarding nodes b_i, i=1,…,8, in Table II. For example, the upper left cell means that node b_1 forwards the data of source s_1 to nodes b_2 and b_7.

The objective throughout the simulation is to maximize the sum utility of source nodes, where source s_i, i=1,…,8, has the utility function given by

    U_s_i(r_s_i)
    =     1.763(r_s_i)^1/6
    -20.718(r_s_i)^2/6
    +88.568(r_s_i)^3/6
       -169.102(r_s_i)^4/6
    +145.167(r_s_i)^5/6-44.677(r_s_i)^6/6.

U_s_i(r_s_i) is
a step-like non-concave polynomial-like function. We consider to optimize a step-like non-concave function, because it is more likely to describe the video quality perceived by a user in a video streaming application <cit.>.
Moreover, we obtain the resource constraints information from Fig. 2 and Table. II, and impose the lower and upper bounds on the aggregate data rate of each user as ξ_s_i=0 and ζ_s_i=10, i=1,…,8, respectively.


Given the network topology shown in Fig. 2, we choose the step-size parameters to satisfy the convergence condition set forth by Theorem 2. All step-size parameters are chosen locally using local information.
Fig. 3 shows the performance of Algorithm DPDA for these step-size parameters. It can be seen that the utility function converges to the optimal one, which is obtained by using Genetic Algorithm while assuming the availability of global information. Although all the computations of DPDA are performed locally at each node, it attains almost the same network utility obtained by a centralized optimization algorithm. This implies that the iterate sequence of Algorithm DPDA can indeed converge to the optimal traffic allocation.

Fig. 4 shows the representative data rate trajectories for MRGUBS flows belonging to source nodes s_3
and s_4. Both data rate sequences are generated by DPDA. It can be seen from Fig. 4 that the MRGUBS requirements are satisfied.












§ CONCLUSIONS AND DIRECTIONS FOR FUTURE RESEARCH

In this paper, we proposed a distributed traffic allocation algorithm, i.e., DPDA, to allow distributed optimal traffic engineering in a connectionless autonomous network. DPDA is distributed and converges at a 𝑂(1/K) rate, where K is the number of iterations. Moreover, numerical simulation results showed that the behavior of DPDA mimics the optimal traffic distribution.

The results presented in this paper are just the first step towards the implementation of an optimal fast distributed algorithm for traffic engineering. There are many issues that need further consideration. In particular, efforts should be put on testing the implementation in large-scale network settings.




























§ APPENDIX A. PRELIMINARY RESULTS AND PROOF OF PROPOSITION 1

In this Appendix, we include some results from real analysis theory and the main steps of proving Proposition 1.







 §.§ Preliminary results

In this subsection, we first recall some results from real analysis theory which are fundamental for the traffic allocation in connectionless networks.

Let f be an arbitrary real-valued function, ℱ be a compact set, not necessarily convex, and μ be a probability measure. Then,

    xinf{ f(x) : x ∈ℱ} =
    μinf{∫ f d μ : supp(μ) ⊂ℱ},

where supp(μ) denotes the support of the measure μ.


For the sake of completeness, we briefly mention the main steps of this well-known fact. Let x^⋆∈ℱ be a minimizer of f such that f(x)≥ f(x^⋆) for every x ∈ℱ. Then, we have ∫ f d μ≥ f(x^⋆) hold for every probability measure μ with supp(μ) ⊂ℱ. That is to say, we have the following inequality hold

    xinf{ f(x) : x ∈ℱ}≤μinf{∫ f d μ : supp(μ) ⊂ℱ}.


On the other hand, we have ∫ f d δ_x^⋆= f(x^⋆), where δ_x^⋆ is the Dirac measure of x^⋆ on the set ℱ. Since δ_x^⋆ is a particular probability measure with supp(δ_x^⋆) ⊂ℱ and ∫ f d δ_x^⋆= f(x^⋆), we have 

    xinf{ f(x) : x ∈ℱ}≥μinf{∫ f d μ : supp(μ) ⊂ℱ}.


In conclusion, the result of Lemma<ref> is established by (<ref>) and (<ref>).


We proceed with the following theorem <cit.> that provides necessary and sufficient conditions for the existence of Borel measures whose support is included in bounded symmetric intervals of the real line.



Given a sequence 𝐭≜{t_j}_j=1^ℓ and a scalar ϵ>0, there exists a Borel measure μ(.) with support contained in 𝒴≐ [-ϵ,ϵ] such that μ(𝒴)=1 and
t_j=E_μ[y^j]=∫_𝒴 y^jμ(dy) is true if and only if


  * when ℓ=2k+1 (odd case), the following holds

    ϵ𝐌(0,2k,𝐭) ≽𝐌(1,2k+1,𝐭)  
       𝐌(1,2k+1,𝐭) ≽ -ϵ𝐌(0,2k,𝐭),



  * when ℓ=2k (even case), the following holds

    𝐌(0,2k,𝐭) ≽ 0  
       ϵ^2𝐌(0,2k-2,𝐭) ≽𝐌(2,2k,𝐭),


where 𝐌(k,k+2h,𝐭) ∈ℝ^(h+1) × (h+1) is a Hankel matrix of the form

    𝐌(k,k+2h,𝐭) =
    [[     t_k   t_k+1       …   t_k+h;   t_k+1       ⋱       ⋱ t_k+h+1;       ⋮       ⋱       ⋱       ⋮;   t_k+h       …       …  t_k+2h ]].

and t_0=1.



Direct application of Theorem III.2.3 and Theorem III.2.4 in <cit.>.








 §.§ Proof of Proposition 1





We note that the problem can be converted into a  polynomial optimization form with a change of variables y_s=(r_s)^(1/ℓ).
The equivalent problem is stated as follows.

    𝐱maximize      ∑_s∈𝒮∑_j=0^ℓ p_s,j(y_s)^j
       subject to       y_s≤ (r_s)^(1/ℓ),  s ∈𝒮,
             ∑_i ∈ℐ_b,l^out x_i,b,l^out+∑_i ∈ℐ_b,l^inx_i,e_l(b),l^out≤ c_l, l ∈ℒ_b, b ∈ℬ,
             ∑_l ∈ℒ_b,i^out x_i,b,l^out- ∑_l̃∈ℒ_bx_i,e_l̃(b),l̃^out=0, i ∈ℐ_b,  b ∈ℬ,
              ( 𝐱_s^out, r_s)∈𝒳_s,  s ∈𝒮,
              x_i,b,l^out≥ 0, i∈ℐ_b,l^out, l∈ℒ_b,  b∈ℬ.

Note that the feasible set in (<ref>) is convex.
However, the equivalent problem is still a non-convex problem, because of the non-concavity of the utility function. Then, instead of working with y_s, we optimize over moments of probability distributions in the space of y_s. More precisely, suppose y_s is a random variable and we denote by m_s,j the j-th moment of y_s for some probability measure μ, i.e., m_s,j=E_μ[y_s^j].

Now, we consider transforming problem (<ref>) into an optimization problem over the space of probability measures of y_s with a support contained in the feasible set of
(<ref>). 

  
  * Based on Lemma 1, the objective function becomes

    ∫∑_s∈𝒮∑_j=0^ℓ p_s,j(y_s)^j d μ_i
    = ∑_s ∈𝒮𝐩_s^T𝐦_s.

  
  * The first three constraints in (<ref>) are justified by Theorem 1.
  
  * We use the set of constraints

    m_s,j≤ r_s^j / ℓ,   j ∈{1,,α}

to approximate the constraint y_s≤ (r_s)^(1/ℓ).
  
  * The left hand of each constraint ∑_i ∈ℐ_b,l^out x_i,b,l^out+∑_i ∈ℐ_b,l^inx_i,e_l(b),l^out≤ c_l for l∈ℒ_b, b∈ℬ,
is written as 1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out. In a similar way,
we rewrite constrains ∑_l ∈ℒ_b,i^out x_i,b,l^out- ∑_l̃∈ℒ_b,i^outx_i,e_l̃(b),l̃^out=0, i ∈ℐ_b,  b ∈ℬ in a matrix form, i.e., B𝐱=0.




In conclusion, Lemma 1, Theorem 1 and (<ref>) establish the result of Proposition 1.




















§ APPENDIX B. DERIVATION OF DPDA

The constrains set of convex relaxation (<ref>) consists of local constraints, e.g., capacity constraints and global constraints, e.g., flow conservation constraints through nodes. The existence of global constraints renders difficulty for us to solve problem (<ref>) in a distributed fashion. However, 
the primal-dual method, proposed by Chambolle and Pock in <cit.> for solving convex-concave saddle point problems makes it possible. This algorithm can be adapted to solve the multi-agent consensus optimization problem as discussed in <cit.>. We also use the distributed primal-dual algorithm in <cit.> to solve our traffic allocation problem (<ref>). This Appendix aims at developing the distributed algorithm that converges to the solution of (<ref>).


The optimization problem (<ref>) can be compactly stated as

    𝐱maximize      ∑_s ∈𝒮𝐩_s^T𝐦_s
       subject to      1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out-c_l≤ 0, l∈ℒ_b, b ∈ℬ
             B𝐱=0,
             (𝐱_s^out, 𝐦_s, r_s)∈𝒜_s,  s∈𝒮,
             𝐱_b,l^out≽ 0, l∈ℒ_b,  b∈ℬ,

where 𝒜_s are the set of local constraints for each source node s ∈𝒮, as defined in (<ref>).


We introduce the convex-concave saddle-point form of the primal problem (<ref>),

    𝐱minλ, θmaxL(𝐱, λ, θ),

 where L(𝐱, λ, θ) is the Lagrangian function given by

    L(𝐱, λ, θ)
    =   -∑_s ∈𝒮 (𝐩_s^T𝐦_s
    -I_𝒜_s(𝐱_s^out, 𝐦_s, r_s)  )
       +∑_b∈ℬ∑_l∈ℒ_b I_R_+^|ℐ_b,l^out|(𝐱_b,l^out)-∑_b∈ℬ∑_l∈ℒ_b I_R_+(λ_b,l)
    
       +∑_b∈ℬ∑_l∈ℒ_b1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out - c_l, λ_b,l
       +B𝐱, θ.

θ∈R^∑_b∈ℬ |I_b| is the vector of dual variables associated with the flow conservation constraint at nodes B𝐱=0. Given l∈ℒ_b and b∈ℬ, the dual variable λ_b,l is introduced for the capacity inequality constrains 1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out≤ c_l. Moreover, λ=[λ_b,l]_l∈ℒ_b,b∈ℬ.

Now, given the initial iterates 𝐱^0, λ^0, θ^0 and parameters
γ>0, τ_s>0 for all s∈𝒮, τ_i,b,l>0, κ_b,l>0 for all i∈ℐ_b,l^out, l∈ℒ_b and b∈ℬ,
we present the following primal-dual iterations to solve (<ref>):

    𝐱^k+1←   𝐱argmin-∑_s ∈𝒮 (𝐩_s^T𝐦_s-I_𝒜_s(𝐱_s^out, 𝐦_s, r_s)  )
       +∑_b∈ℬ∑_l∈ℒ_b I_R_+^|ℐ_b,l^out|(𝐱_b,l^out)
    
       +∑_b∈ℬ∑_l∈ℒ_b1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out - c_l, λ_b,l^k
       +B𝐱, θ^k+∑_b∈ℬ∑_l∈ℒ_b∑_i∈ℐ_b,l^out1/2τ_i,b,l(x_i,b,l^out-x_i,b,l^out,k)^2
        + ∑_s ∈𝒮1/2τ_s( (x_s,l^out- x_s,l^out,k)^2 + 𝐦_s-𝐦_s^k_2^2 +  (r_s-r_s^k)_2^2   )
    ;
    λ_b, l^k+1←   λ_b, largmax -  I_R_+(λ_b,l)+⟨1_b,l(2 𝐱_b,l^out, k+1-𝐱_b,l^out, k)
       +
    δ_b,l(2 𝐱_e_l(b),l^out, k+1-𝐱_e_l(b),l^out, k) - c_l, λ_b, l⟩
       -1/2κ_b, l(λ_b, l-λ_b, l^k)^2,  l∈ℒ_b,   b∈ℬ;
    θ^k+1←   θargmaxB(2𝐱^k+1-𝐱^k), θ-1/2γθ-θ^k_2^2
        =θ^k+γ B (2𝐱^k+1-𝐱^k).



Although the convergence to the optimal traffic allocation is guaranteed under the primal-dual method, it is still not a distributed algorithm. In fact, solving the optimization problem involved in the primal variables 𝐱^k+1 update rule requires global information about the network due to the presence of the term  B𝐱, θ^k, which is associated with the flow conservation constraints at nodes. Moreover, computing the term 1_b,l𝐱_b,l^out+δ_b,l𝐱_e_l(b),l^out - c_l, λ_b,l^k,  l∈ℒ_b, b∈ℬ forces neighboring nodes to exchange information, because bidirectional links are allowed to exist in the model. This fact hinder us from directly implementing the primal-dual iterations. Nevertheless, we exploit the structure of the inner product B𝐱, θ^k and note that this term is a summation of local linear functions of the local variables. In addition, the sending data rates of neighboring nodes is local information. These observations indicates that it is possible to develop an optimal decentralized traffic allocation algorithm.






Using recursion in θ update rule in (<ref>), we can write θ^k+1 as a partial summation of previous primal variable 𝐱^k iterations, i.e., θ^k=θ^0+γ∑_n=0^k-1 B (2𝐱^n+1-𝐱^n).
Let θ^0 be
γ B 𝐱^out,0, 𝐳^0 be 𝐱^0 and 𝐳^k ≜𝐱^k+∑_n=1^k𝐱^n for k≥ 1.
Then we get

    B𝐱, θ^k
    =    γ𝐱^out,B^T B𝐳^k
    
    =   γ∑_b∈𝒩∑_l ∈ℒ_b∑_i ∈ I_b,l^out x_i, b, l^out( ∑_l̃∈ℒ_b,i^out z_i, b, l̃^k-  ∑_l̅∈ℒ_b z_i, e_l̅(b), l̅^k).


The quadratic operation for updating λ_b,l^k+1 in (<ref>) entails solving the following projection
problem:

    λ_b,l^k+1←𝑃_R_+(λ_b,l^k+ 1_b,l(2 𝐱_b,l^out, k+1-𝐱_b,l^out, k)
    +δ_b,l(2 𝐱_e_l(b),l^out, k+1-𝐱_e_l(b),l^out, k) - c_l).




Substituting (<ref>) and (<ref>) into (<ref>) yields a distributed traffic allocation algorithm shown in Algorithm 1.




§ APPENDIX C. PROOF OF THEOREM 2

In this section, we present the Proof of Theorem 2.



Due to space limitations, we only prove that if conditions (<ref>) and (<ref>) hold, the following inequality is true:

    Q(A,B)≜[ [  D_τ -A^T -B^T;   -A  D_κ    0;   -B    0  D_γ ]] ≽ 0

where
D_κ≜diag([1/κ_b, l]_ l∈ℒ_b, b∈ℬ), D_γ≜1/γ I_∑_b∈𝒩 |I_b|, and
D_τ≜diag([v_sτ^T, v_bτ^T]^T ) where v_sτ≜ [1/τ_s1_(|ℒ_s| +ℓ+2)× 1]_s∈𝒮 and v_bτ=[1/τ_i,b,l]_ i∈ℐ_b,l^out, l∈ℒ_b, b∈ℬ. Moreover,
A≜diag([A_l]_l∈ℒ_b, b∈ℬ), where A_ l is a row vector with the same dimension as vector variable 𝐱, and the i-th entry of vector A_ l, equals to 1 if the data rate denoted by the i-th element is transported through link l, 0 otherwise.




Based on “Schur complement Lemma", we have Q(A,B) ≽ 0 holds if and only if

    [[  D_τ -A^T;   -A  D_κ ]]
    -γ[[ B^TB    0;    0    0 ]]≽ 0.

Moreover, since D_κ≽ 0, again using “Schur complement Lemma", one can conclude that (<ref>) holds if and only if

    D_τ- γ B^TB- A^TD_κ^-1A ≽ 0.

 Denote matrix B^TB by Ω, and we can write
Ω into the sum of two matrices, i.e.,

    Ω=diag([ω_i,b,l]_i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩)+E,

where ω_i,b,l=1 if node b∈𝒮⋃ℬ_e where ℬ_e is the set of nodes that forward traffic to destination nodes, ω_i,b,l=2 if b∈ℬ⋂ℬ_e^c where ℬ_e^c is the complement of ℬ_e, otherwise 0. Also, all the diagonal elements of matrix E are equal to 0, and the  non-diagonal element E_(i,b_1,l_1),(j,b_2,l_2), corresponding to data rates x_i,b_1,l_1^out∈𝐱 and x_i,b_2,l_2^out∈𝐱, equals to 1 if both data rates belong to the same source node and they are forwarded from the same node, i.e., i=j∈ℐ and b_1=b_2∈𝒩, -1 if both data rates belong to the same source node and nodes b_1 and b_2 are neighboring, and 0 otherwise. Based on “Gershgorin Circle Theorem" <cit.>, we have

    diag([ω_i,b,l]_ i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩)+
    diag([d_i,b,l]_ i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩)
     -E ≽ 0,
 since d_i,b,l is chosen to be large enough. Therefore,

    Ω≼  2diag([ω_i,b,l]_ i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩)+
    diag([d_i,b,l]_ i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩)

Moreover,

    Ω≼diag([4+d_i,b,l]_ i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩).

Hence, it is sufficient to have

    τ- γdiag([4+d_i,b,l]_ i ∈ℐ_b,l^out, l∈ℒ_b, b∈𝒩) - A^TD_κ^-1A ≽ 0 ,

 and this condition holds if the inequalities (<ref>) and (<ref>) in the statement of Theorem 1 are true.


Let (𝐱^⋆,λ^⋆,θ^⋆) be an arbitrary saddle-point for the Lagrange function of problem (<ref>), and
{𝐱^k}_k≥ 0 be the iterate sequence generated using Algorithm DPDA, initialized from an arbitrary 𝐱^0 and [λ_b,l^0]_l∈ℒ_b,b∈ℬ=0.
Denote the average of sending data rates
by 𝐱̅^K≜1/K∑_k=1^K 𝐱^k, where K≥ 1.
Then, following the proof in <cit.>, we have that {𝐱̅^K} converges to the maximum of the utility function of the problem (<ref>) subject to the resource allocation constraints. In particular, the following error bounds hold for all K≥ 1:

    θ^⋆B𝐱̅^K+    ∑_b∈ℬ∑_l∈ℒ_bλ^⋆_b,l h(𝐱̅_b,l^out,𝐱̅_e_l(b),l^out)
    ≤Θ_1/K,
       |  ∑_s∈𝒮𝐩_s^T(𝐦̅_s-𝐦_s^⋆)|
    ≤Θ_1/K,

where h(𝐱̅_b,l^out,𝐱̅_e_l(b),l^out) denotes the distance function d_R_-
(1_b,l𝐱̅_b,l^out + δ_b,l𝐱̅_e_l(b),l^out-c_l), and Θ_1≜2/γθ^⋆ ^2-γ/2B𝐱̅^0^2
+ ∑_b∈ℬ∑_l∈ℒ_b (∑_i∈ℐ_b,l^out1/2τ_i,b,l(x_i,b,l^out,⋆-x_i,b,l^out,0)^2+ 1/2κ_b,l (λ^⋆_b,l)^2 ) + ∑_s∈𝒮1/2τ_s(𝐦_s^⋆-𝐦_s^0 ^2+(r_s^⋆-r_s^0  )^2+ ∑_l∈ℒ_s(x_s,l^out,⋆-x_s,l^out,0)^2 ).





§ REFERENCES

99


Kelly_JORS_1998 F. P. Frank, A. K. Maulloo and D. K. H. Tan,  Rate control for communication networks: shadow prices, proportional fairness and stability, Journal of the Operational Research society, Springer, vol. 49, no. 3, pp. 237–252, 1998.


lagoa_2004_N C. M. Lagoa, H. Che and B. A Movsichoff, Adaptive control algorithm for dencentralized optimal traffic engineering in the Internet, IEEE/ACM Transactions on  Networking, vol. 12, no. 3, pp. 415–428, June 2004.



lagoa_hop B. A. Movsichoff, A. Bernardo, C. M. Lagoa and H. Che, Decentralized optimal traffic engineering in connectionless networks, IEEE Journal on Selected Areas in Communications, vol. 23, no. 2, pp. 293–303, 2005.




beck2013optimalA. Beck, A.  Nedic, A. Ozdaglar and M. Teboulle, Optimal distributed gradient methods for network resource allocation problems, Submitted for publication, 2013.


nekouei E. Nekouei, G. Nair and T.  Alpcan, Convergence analysis of quantized primal-dual algorithm in quadratic network utility maximization problems, IEEE Conference on Decision and Control CDC, pp. 2655–2660, 2015.



Yin_2015_SC X. Q. Yin, A. Jindal, V. Sekar and B. Sinopoli, A control-theoretic approach for dynamic adaptive video streaming over HTTP, SIGCOMM, London, United Kingdom, pp. 325–338, August 2015.





Fazel_CDC_2005 M. Fazel and M. Chiang, Network utility maximization with nonconcave utilities using sum-of-squares method, Proceedings of the 44th IEEE Conference on Decision and Control, 2005, pp. 1867–1874.






Hande_Net_2007 P. Hande, S. Y. Zhang and M. Chiang, Distributed rate allocation for inelastic flows, IEEE/ACM Transactions on Networking, vol. 15, no. 6, pp. 1240–1253, 2007.







lasserre J. B. Lasserre, Global optimization with polynomials and the problem of moments, SIAM Journal on Optimization, vol. 11, no. 3, pp. 796–817, 2001.






laurent M. Laurent, Sums of squares, moment matrices and optimization over polynomials, Emerging applications of algebraic geometry, Springer New York, 2009: 157–270.







Aybat_arxiv_2016 N. S. Aybat and E. F. Hamedani, A primal-dual method for conic constrained distributed optimization problems, Advances in Neural Information Processing Systems 29, edited by D. D. Lee and M. Sugiyama and U. V. Luxburg and I. Guyon and R. Garnett, Curran Associates, Inc., pp. 5049–5057, 2016, http://papers.nips.cc/paper/6242-a-primal-dual-method-for-conic-constrained-distributed-optimization-problems.pd.


ozay N. Ozay, C. M. Lagoa and M. Sznaier, Set membership identification of switched linear systems with known number of subsystems, Automatica, vol. 51, pp. 180–191, 2015.


krein M. G. krein and A. A. Nudelman, The markov moment problem and extremal problems, volume 50 of translations of mathematical monographs, American Mathematical Society, Providence, Rhode Island, 1977.




Chambolle_mathe_2015 A. Chambolle and T. Pock, On the ergodic convergence rates of a first-order primal-dual algorithm, Mathematical Programming, Springer, 2015, pp. 1–35.



Golub:1996:MC:248979 G. H. Golub and C. F. Van Loan, Matrix Computations (3rd Ed.), Johns Hopkins University Press, Baltimore, MD, USA, 1996.






