WIXLIB

IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 26, NO. 4, PP. 1605-1618, JUNE 2018contrast, the first order Lagrangian dual type methods do notneed global network topology information. Rather, each nodeonly needs the queue length information of its neighbors.This paper proposes a new first order Lagrangian dualtype backpressure algorithm that is as simple as the existingalgorithms in [5], [6], [8] but has a better utility-delay tradeoff.The new backpressue algorithm achieves a vanishing utilityoptimality gap that decays like O(1/t), where t is the numberof iterations. It also guarantees that the queue length at eachnode is always bounded by a fixed constant of the same orderas the optimal Lagrange multiplier of the network optimizationproblem. This improves on the utility-delay tradeoffs of priorwork. In particular, it improves the steady-state [O( ), O(1/ 2 )]utility-delay tradeoff in [6] and the [O( ), O(1/ )] utility-delaytradeoff of the drift-plus-penalty algorithm in [5], both ofwhich yield an unbounded queue length to have a vanishingutility optimality gap. Indeed, the steady-state utility-delaytradeoff of our algorithm is [0, O(1)]. The convergence timeto reach this limiting performance is also faster than priorwork. Our algorithm achieves a zero utility gap and O(1) queuelengths with a fast O(1/t) convergence rate.The new backpressure algorithm differs from existing firstorder backpressure algorithms in the following aspects:1) The “backpressure” quantities in this paper are withrespect to newly introduced weights. These are differentfrom queues used in other backpressure algorithms, butcan still be locally tracked and updated.2) The rate control and routing decision rule involves aquadratic term that is similar to a term used in proximalalgorithms [20].Note that the benefit of introducing a quadratic term innetwork optimization has been observed in [21]. Work [21]developed a distributive rate control algorithm for networkutility maximization (NUM) problems with given routing pathsthat can be reformulated as a special case of the problemtreated in this paper. The algorithm of [21] considers a fixedset of predetermined paths for each session and does notscale well when treating all (typically exponentially many)possible paths of a general network. The algorithm proposedin [21] is not a backpressure type and and uses virtual queuesrather than actual queues. It is interesting to note that thevirtual queues in [21] remain O(1).1 However, O(1) virtualqueue size does not imply the achieved utility has a fastconvergence (the convergence rate of [21] remains an openquestion). In contrast, the algorithm in the current paper usesboth a quadratic term and a modified weight to get O(1) actualqueue size with a fast O(1/t) convergence rate.In our conference version [1], the proposed backpressure algorithm has a global algorithm parameter α, which is requiredto be chosen based on the total number of links and sessionsin the network. In this paper, the proposed backpressure1 A technique in the proof of Prop 6.3.2 in [22] suggests that in the case offixed path routing of [21] and under several other assumptions that includethe implementation of a closest-to-source priority rule, there exists a constantthat bounds the deviation between virtual and actual queues. The technique in[22] relies on a fixed-path assumption and has bounds that can grow at eachstage k of the path.algorithm allows each node n to locally determine its algorithmparameter αn based on the number of local link connections.The network algorithm developed in this paper is based onour recent work [23] that develops a new convex programmingmethod with an O(1/t) error decay for general constrainedconvex programs (including problems that are nonsmooth, notstrongly convex and have nonlinear constraints). After our conference version [1] and the arXiv preprint [24] of this paper,Wang and Shroff in [25] studied the same joint rate control androuting problem considered in this paper and propose anotherbackpressure algorithm by using the generalized alternatingdirection method of multipliers (ADMM) from [26], [27], [28].While that ADMM approach is very different from ours, theobtained backpressure algorithm in [25] is remarkably similarto our algorithm in this paper: It uses identical virtual queues(up to a constant scaling) and solves source rate and link ratesubproblems with identical structures. One difference is thetwo algorithms use different weights for these subproblems.The other main difference is that at each iteration the algorithmin [25] requires to update link rates first and then updatesource rates based on the new link rate values, while ouralgorithm can update link rates and source rates in parallelsince these subproblems are fully decoupled. This is becausethe algorithm in [25] is restricted to the sequential updateprocedure of ADMM, while our algorithm uses the parallelmethods developed by us in [23] together with the naturalseparability of the linear node flow balance constraints.The algorithm in [25] achieves source rates that satisfyx[t] x , although a convergence rate for general utilityfunctions is not given. 2 In contrast, our analysis provesan explicit O(1/t) optimality deviation for the time averagenetwork utility. This time average starts at time 0 and providesa performance guarantee for online implementation duringthe entire network operation. However, a remarkable propertyof the algorithm in [25] is that when the utility is smooth(i.e., differentiable with Lipschitz continuous gradients) andstrongly convex, then the algorithm in [25] can ensure thefinal value of the iterate x[t] converges exponentially fast toan optimal value x . It remains a promising research directionto investigate whether the fast final iteration convergence alsoholds for our algorithm, and to investigate the effect of theparallel source and link rate update property and differentupdate rules for weight parameters in the subproblems.II. S YSTEM M ODEL AND P ROBLEM F ORMULATIONConsider a slotted data network with normalized time slotst {0, 1, 2, . . .}. This network is represented by a graphG (N, L), where N is the set of nodes and L N Nis the set of directed links. Let N N and L L. Thisnetwork is shared by F end-to-end sessions denoted by a setF . For each end-to-end session f F , the source node Src( f )2 According to [26], the authors in [25] conclude their algorithm has2“O(1/t) non-ergodic convergence” defined as kx[t] x[t 1] k O(1/t),or equivalently, kx[t] x[t 1] k O(1/ t), in [26]. However, this unusualdefinition of non-ergodic convergence rate does not imply any convergencerate of kx[t] x k or utility optimality gap U(x[t]) U(x ) and hence does notprovide insight into convergence behavior that network optimization concerns.Thus, it has nothing to do with the O(1/t) convergence to optimal utility weshow in the current paper.

IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 26, NO. 4, PP. 1605-1618, JUNE 2018and destination node Dst( f ) are given but the routes are notspecified. Each session f has a continuous and concave utilityfunction U f (x f ) that represents the “satisfaction” received byaccepting x f amount of data for session f into the networkat each slot. Unlike [6], [11] where U f (·) is assumed tobe differentiable and strongly concave, this paper considersgeneral concave utility functions U f (·), including those thatare neither differentiable nor strongly concave. Formally, eachutility function U f is defined over an interval dom(U f ), calledthe domain of the function.Denote the capacity of link l as Cl and assume it is a fixed(f )and positive constant.3 Define µl as the amount of session f ’sdata routed at link l that is to be determined by our algorithm.Note that in general, the network may be configured such thatsome session f is forbidden to use link l. For each link l,define Sl F as the set of sessions that are allowed to uselink l. The case of unrestricted routing is treated by definingSl F for all links l.(f )Note that if l (n, m) with n, m N , then µl and Cl(f )can also be respectively written as µ(n,m) and C(n,m) . For eachnode n N , denote the sets of its incoming links and outgoinglinks as I(n) and O(n), respectively. Note that x f , f F and(f )µl , l L, f F are the decision variables of a joint ratecontrol and routing algorithm. If the global network topologyinformation is available, the optimal joint rate control androuting can be formulated as the following multi-commoditynetwork flow problem:ÕmaxU f (x f )(1)(f )x f ,µlf Fs.t.x f 1 {n Src( f )} Õ(f )µl l I(n)Õ(f )µll O(n) f F , n N \ {Dst( f )}Õf F(f )µl(f )µl(f )µl(2) Cl, l L,(3) 0, l L, f Sl,(4) 0, l L, f F \ Sl,(5)x f dom(U f ), f F(6)where 1 {·} is an indicator function; (2) represents the nodeflow conservation constraints relaxed by replacing the equalitywith an inequality, meaning that the total rate of flow f intonode n is less than or equal to the total rate of flow f outof the node (since, in principle, we can always send fakedata for departure links when the inequality is loose); and (3)represents link capacity constraints. Note that for each flow f ,there is no constraint (2) at its destination node Dst( f ) sinceall incoming data are consumed by this node.The above formulation includes network utility maximization with fixed paths as special cases. In the case when eachsession only has one single given path, e.g., the network utilitymaximization problem considered in [30], we could modify thesets Sl used in constraints (4) and (5) to reflect this fact. Forexample, if link l1 is only used for sessions f1 and f2 , then3 As stated in [11], this is a suitable model for wireline networks andwireless networks with fixed transmission power and orthogonal channels.Sl1 { f1, f2 }. Similarly, the case [21] where each session isrestricted to using links from a set of predefined paths can betreated by modifying the sets Sl accordingly. See AppendixA for more discussions.The solution to problem (1)-(6) corresponds to the optimaljoint rate control and routing. However, to solve this convexprogram at a single computer, we need to know the globalnetwork topology and the solution is a centralized one, whichis not practical for large data networks. As observed in [10],[6], [8], [11], various versions of backpressure algorithms canbe interpreted as distributed solutions to problem (1)-(6) fromfirst order Lagrangian dual type methods.Assumption 1: (Feasibility) Problem (1)-(6) has at least one( f ), optimal solution vector [x f ; µl ] f F,l L .Assumption 2: (Existence of Lagrange multipliers) Assume the convex program (1)-(6) has Lagrange multipliersattaining the strong duality. Specifically, define convex set(f )C {[x f ; µl ] f F,l L : (3)-(6) hold}. Assume there exists( f ), a Lagrange multiplier vector λ [λn ] f F,n N\{Dst( f )} 0such thatq(λ ) max{(1) : (2)-(6)} Íwhereq(λ) max[x ;µ( f ) ] C f F U f (x f )f lÍÍÍ(f )(f ) l I(n) µlf )} λn x f 1 {n Src( f )}Í f F n N\{Dst((f ) is the Lagrangian dual function of probleml O(n) µl(1)-(6) by treating (3)-(6) as a convex set constraint.Assumptions 1 and 2 hold in most cases of interest. Forexample, Slater’s condition guarantees Assumption 2. Sincethe constraints (2)-(6) are linear, Proposition 6.4.2 in [31]ensures that Lagrange multipliers exist whenever constraints(2)-(6) are feasible and when the utility functions U f areeither defined over open sets (such as U f (x) log(x) withdom(U f ) (0, )) or can be concavely extended to opensets, meaning that there is an 0 and a concave functionef : ( , ) R such that Uef (x) U f (x) whenever x 0.4UFact 1: (Replacing inequality with equality) If Assumption 1 holds, problem (1)-(6) has an optimal solution vector( f ), [x f ; µl ] f F,l L such that all constraints (2) take equalities.(f )Proof: Note that each µl can appear on the left side inat most one constraint (2) and appear on the right side in at( f ), most one constraint (2). Let [x f , µl ] f F,l L be an optimalsolution vector such that at least one inequality constraint (2)( f ), is loose. Note that we can reduce the value of µlon theright side of a loose (2) until either that constraint holds with( f ), equality, or until µlreduces to 0. The objective functionvalue does not change, and no constraints are violated. Wecan repeat the process until all inequality constraints (2) aretight.4 If dom(U ) [0, ), such concave extension is possible if the rightfderivative of U f at x 0 is finite (such as for U f (x) log(1 x) or U f (x) min[x, 3]). Such an extension is impossible for the example U f (x) xbecause the slope is infinite at x 0. Nevertheless, Lagrange multipliersoften exist even for these utility functions, such as when Slater’s conditionholds [31].

IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 26, NO. 4, PP. 1605-1618, JUNE 2018III. T HE N EW BACKPRESURE A LGORITHMA. Discussion of Various Queueing ModelsAt each node, an independent queue backlog is maintainedfor each session. At each slot t, let x f [t] be the source session(f )rates; and let µl [t] be the link session rates. Some prior work(f )enforces the constraint (2) via virtual queues Yn [t] of thefollowing form:n(f )(f )Yn [t 1] max Yn [t] x f [t]1 {n Src( f )}oÕÕ(f )(f ) µl [t] µl [t], 0 .(7)l I(n)l O(n)While this virtual equation is a meaningful approximation, itdiffers from reality in that new injected data are allowed tobe transmitted immediately, or equivalently, a single packetis allowed to enter and leave many nodes within the sameslot. Further, there is no clear connection between the virtual(f )queues Yn [t] in (7) and the actual queues in the network.Indeed, it is easy to construct examples that show there can be(f )an arbitrarily large difference between the Yn [t] value in (7)and the physical queue size in actual networks (see AppendixB).(f )An actual queueing network has queues Zn [t] with thefollowing dynamics:noÕ(f )(f )(f )Zn [t 1] max Zn [t] µl [t], 0In short, our algorithm can be faithfully implemented withrespect to actual queueing networks, and converges to exactoptimality on those networks.The next lemma shows that if an algorithm can guarantee(f )virtual queues Q n [t] defined in (9) are bounded, then actualphysical queues satisfying (8) are also bounded.Lemma 1: Consider a network flow problem describedby problem (1)-(6). For all l L and f F , let(f )µl [t], x f [t] be decisions yielded by a dynamic algorithm.(f )(f )(f )Suppose Yn [t], Zn [t], Q n [t] evolve by (7)-(9) with initial(f )(f )(f )conditions Yn [0] Zn [0] Q n [0] 0. If there exists a(f )constant B 0 such that Q n [t] B, t, then we also haveÍ(f )1) Zn [t] 2B l O(n) Cl for all t {0, 1, 2, . . .}.Í(f )2) Yn [t] 2B l O(n) Cl for all t {0, 1, 2, . . .}.Proof:1) Fix f F , n N \ {Dst( f )}. Define an auxiliaryb( f )b( f )virtualÍ queue Q n [t] that is initialized by Q n [0] B l O(n) Cl and evolves according to (9). It followsb(nf ) [t] Q(nf ) [t] B Íl O(n) Cl, t. Since Q(nf ) [t] that Qb(nf ) [t] Íl O(n) Cl B, t by assumption, we have QÍ(f )b( f )l O(n) µl [t], t. This implies that Q n [t] also satisfies:onÕ(f )b(nf ) [t] b(nf ) [t 1] max Qµl [t], 0Ql O(n) x f [t]1 {n Src( f )} l O(n) x f [t]1 {n Src( f )} (f )µl [t].(8)This is faithful to actual queue dynamics and does not allowdata to be retransmitted over multiple hops in one slot. Notethat (8) is an inequality because the new arrivals from otherÍ(f )nodes may be strictly less than l I(n) µl [t] because thoseother nodes may not have enough backlog to send. The model(8) allows for any decisions to be made to fill the transmissionÍ(f )(f )(f )values µl [t] in the case that Zn [t] l O(n) µl [t],provided that (8) holds.This paper develops an algorithm that converges to theoptimal utility defined by problem (1)-(6), and that producesworst-case bounded queues on the actual queueing network,that is, with actual queues that evolve as given in (8). Tobegin, it is convenient to introduce the following virtual queueequation(f )(f )Õl O(n)(f )µl [t] x f [t]1 {n Src( f )} Õ(f )µl [t]l I(n)(9)(f )(f )µl [t], t (10)l I(n)Õl I(n)Q n [t 1] Q n [t] Õwhere Q n [t] represents a virtual queue value associatedwith session f at node n. At first glance, this model (9)appears to be only an approximation, perhaps even a worse(f )approximation than (7), because it allows the Q n [t] values(f )to be negative. Indeed, we use Q n [t] only as virtual queuesto inform the algorithm and do not treat them as actualqueues. However, this paper shows that using these virtualqueues to choose the µ[t] decisions ensures not only thatthe desired constraints (2) are satisfied, but that the resulting(f )µ[t] decisions create bounded queues Zn [t] in the actualnetwork, where the actual queues evolve according to (8).which is identical to (8) except the inequality is replaced(f )b(nf ) [0]; and Qb(nf ) [t]by an equality. Since Zn [0] 0 Q(f )(f )bsatisfies (10), by inductions, Zn [t] Q n [t], t.b(nf ) [t] Q(nf ) [t] B Íl O(n) Cl, t and Q(nf ) [t] Since Qb(nf ) [t] 2B Íl O(n) Cl, t. It followsB, t, we have QÍ(f )that Zn [t] 2B l O(n) Cl, t.2) The proof of part (2) is similar and is in Appendix C.B. The New Backpressure AlgorithmIn this subsection, we propose a new backpressure algorithmthat yields source session rates x f [t] and link session rates(f )µl [t] at each slot such that the physical queues for eachsession at each node are bounded by a constant and the timeaverage utility satisfiest 1Õ1ÕÕU f (x f [t]) U f (x f ) O(1/t), t 1t τ 0f Ff Fwhere x f are from the optimal solution to (1)-(6). Note thatJensen’s inequality further implies thatÕf FUft 1 Õ1Õx f [τ] U f (x f ) O(1/t), t 1t τ 0f FThe new backpressure algorithm is described in Algorithm1. Similar to existing backpressure algorithms, the updates inAlgorithm 1 at each node n are fully distributed and onlydepend on weights at itself and its neighbor nodes. Unlikeexisting backpressure algorithms, the weights used to update

IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 26, NO. 4, PP. 1605-1618, JUNE 2018(f )decision variables x f [t] and µl [t] are not the virtual queues(f )(f )Q n [t] themselves, rather, they are augmented values Wn [t]equal to the sum of the virtual queues and the amount of netinjected data in the previous slot t 1. In addition, the updatesinvolve an additional quadratic term, which is similar to a termused in proximal algorithms [20]. Note that the source rateupdate (in the second bullet) and the link rate update (in thethird bullet) are fully decoupled and hence can be performedin parallel or in a swapped order.the equation ψ(x f ) 0 and can be found by a bisectionsearch.2) Suppose dom(U f ) (0, ) and U f (x f ) w f log(x f ) forsome weight w f 0. Then:Algorithm 1 The New Backpressure AlgorithmLet αn 0, n N be constant parameters. Initialize(f )(f )x f [ 1] 0, µl [ 1] 0, f F , l L and Q n [0] 0, n N, f F . At each time t {0, 1, 2, . . .}, each noden does the following: For each f F , if node n is not the destination node of(f )session f , i.e., n , Dst( f ), then define weight Wn [t]:Proof: Omitted for brevity.The problem (14)-(17) can be represented as follows by(f )eliminating µ(n,m), f S(n,m) , completing the square andreplacing maximization with minimization. (Note that K S(n,m) F .)(f )l I(n)(11)xf (13)s.t.1Õ(zk ak )22 k 1KÕzk b(18)(19)zk 0, k {1, 2, . . . , K }(20)Lemma 3: The solution to problem (18)-(20) is given byzk max{0, ak θ }, k {1, 2, . . . , K } where θ 0 can befound either by a bisection search (See Appendix D) or byAlgorithm 2 with complexity O(K log K).Proof: A similar problem where (19) is replaced with anequality constraint is considered in [32]. The optimal solutionto this quadratic program is characterized by its KKT conditionand a corresponding algorithm can be developed to obtain itsKKT point. A complete proof is presented in Appendix D.(f )For all (n, m) O(n), choose {µ(n,m) [t], f F } as thesolution to the following convex program:Õ(f )(f ) (f )maxWn [t] Wm [t] µ(n,m)(f )µ(n, m) f F αn αm Õ(f )(f )µ(n,m) µ(n,m) [t 1] 2(14)f Fs.t.Õ(f )µ(n,m) C(n,m)f F(f )µ(n,m)(f )µ(n,m) Kminl O(n)x f dom(U f )4αnk 1If node n is the destination node, i.e., n Dst( f ), then(f )define Wn [t] 0. Notify neighbor nodes (nodes k thatcan send session f data to node n, i.e., k such that(f )f S(k,n) ) about this new Wn [t] value.For each f F , if node n is the source node of sessionf , i.e., n Src( f ), choose x f [t] as the solution to 2(f )max U f (x f ) Wn [t]x f αn x f x f [t 1](12)s.t. 2αn x f [t 1] Wn [t]4αnq(f )(Wn [t] 2αn x f [t 1])2 8αn w f(f )Wn [t] Q n [t] x f [t 1]1 {n Src( f )}ÕÕ(f )(f ) µl [t 1] µl [t 1]. (f )x̂ f (15) 0, f S(n,m)(16) 0, f S(n,m)(17)For each f F , if node n is not the destination of f ,(f )i.e., n , Dst( f ), update virtual queue Q n [t 1] by (9).Algorithm 2 Algorithm to solve problem (18)-(20)ÍK1) Check if k 1max{0, ak } b holds. If yes, let θ 0 and zk max{0, ak }, k {1, 2, . . . , K } and terminatethe algorithm; else, continue to the next step.2) Sort all ak , {1, 2, . . . , K } in a decreasing order π suchthat aπ(1) aπ(2) · · · aπ(K) . Define S0 0.3) For k 1 to KS b Let Sk Sk 1 ak . Let θ kk . If θ 0, aπ(k) θ 0 and aπ(k 1) θ 0,then terminate the loop; else, continue to the nextiteration in the loop.4) Let zk max{0, ak θ }, k {1, 2, . . . , K } andterminate the algorithm.Note that step (3) in Algorithm 2 has complexity O(K) andhence the overall complexity of Algorithm 2 is dominated bythe sorting step (2) with complexity O(K log(K)).C. Almost Closed-Form Updates in Algorithm 1(f )This subsection shows the decisions x f [t] and µl [t] inAlgorithm 1 have either closed-form solutions or “almost”closed-form solutions at each iteration t.Lemma 2: Let x̂ f x f [t] denote the solution to (12)-(13).1) Suppose dom(U f ) [0, ) and U f (x f ) is differentiable.(f )Let ψ(x f ) U f0 (x f ) 2αn x f 2αn x f [t 1] Wn [t].If ψ(0) 0, then x̂ f 0; otherwise x̂ f is the root toIV. P ERFORMANCE A NALYSIS OF A LGORITHM 1A. PreliminariesThis subsection introduces facts from convex analysis andsome useful additional notation.Definition 1 (Strongly Concave Functions): Let Z Rn bea convex set. Function f is said to be strongly concave on

IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 26, NO. 4, PP. 1605-1618, JUNE 2018Z with modulus α if there exists a constant α 0 such thatf (z) 12 αkzk 2 is concave on Z.By the definition of strongly concave functions, it is easy toshow that if f (z) is concave and α 0, then f (z) αkz z0 k 2is strongly concave with modulus 2α for any constant z0 . Thenext lemma shows that a strongly concave function deviatesat least quadratically from its maximum value.Lemma 4 (Corollary 1 in [23]): Let Z Rn be a convexset. Let function f be strongly concave on Z with modulus αand zopt be a global maximum of f on Z. Then, f (zopt ) f (z) α2 kzopt zk 2 for all z Z.(f )Define column vector y [x f ; µl ] f F,l L , which aggregates all optimization variables in problem (1)-(6). For eachf F , n N \ {Dst( f )}, define column vector((f )[x f ; µl ]l I(n) O(n) if n Src( f ),(f )yn (21)(f )[µl ]l I(n) O(n)else,which is composed of the control actions appearing in each(f )constraint (2); and introduce a function with respect to yn asÕÕ(f )(f ) (f )(f )(22)µlgn (yn ) x f 1 {n Src( f )} µl l O(n)l I(n)Thus, constraint (2) can be rewritten as(f )Proof: Fix f F and n N \ Dst( f ), we have 2 1 (f ) 21 (f )Q n [t 1] Q n [t]22 2 1 (f ) 2(a) 1(f )(f ) (f ) Q n [t] gn (yn [t]) Q n [t]22 21 (f ) f(f )(f ) f Q n [t]gn (yn [t]) gn (yn [t])2where (a) follows from (23).By the definition of [t], we have Õ 21(f )(f ) 2 [t] Q n [t 1] Q n [t]2f F,n N\Dst( f )(a) Õ(f )(f )Note

O„ "utility optimality gap with an O„1š 2"queue length. That is, a small utility optimality gap (corresponding to a small ) is available only at the cost of a large queue length. The drift-plus-penalty backpressure algorithm [5], which has the best utility-delay tradeoff among all existing first order