1

A New Backpressure Algorithm for Joint Rate
Control and Routing with Vanishing Utility
Optimality Gaps and Finite Queue Lengths
Hao Yu and Michael J. Neely
Department of Electrical Engineering
University of Southern California

arXiv:1701.04519v1 [cs.NI] 17 Jan 2017

Abstract
The backpressure algorithm has been widely used as a distributed solution to the problem of joint rate control and routing
in multi-hop data networks. By controlling a parameter V in the algorithm, the backpressure algorithm can achieve an arbitrarily
small utility optimality gap. However, this in turn brings in a large queue length at each node and hence causes large network delay.
This phenomenon is known as the fundamental utility-delay tradeoff. The best known utility-delay tradeoff for general networks
is [O(1/V ), O(V )] and is attained by a backpressure algorithm based on a drift-plus-penalty technique. This may suggest that
to achieve an arbitrarily small utility optimality gap, the existing backpressure algorithms necessarily yield an arbitrarily large
queue length. However, this paper proposes a new backpressure algorithm that has a vanishing utility optimality gap, so utility
converges to exact optimality as the algorithm keeps running, while queue lengths are bounded throughout by a finite constant.
The technique uses backpressure and drift concepts with a new method for convex programming.

I. I NTRODUCTION
In multi-hop data networks, the problem of joint rate control and routing is to accept data into the network to maximize
certain utilities and to make routing decisions at each node such that all accepted data are delivered to intended destinations
without overflowing any queue in intermediate nodes. The original backpressure algorithm proposed in the seminal work [1]
by Tassiulas and Ephremides addresses this problem by assuming that incoming data are given and are inside the network
stability region and develops a routing strategy to deliver all incoming data without overflowing any queue. In the context of
[1], there is essentially no utility maximization consideration in the network. The backpressure algorithm is further extended
by a drift-plus-penalty technique to deal with data network with both utility maximization and queue stability considerations
[2], [3], [4]. Alternative extensions for both utility maximization and queue stabilization are developed in [5], [6], [7], [8].
The above extended backpressure algorithms have different dynamics and/or may yield different utility-delay tradeoff results.
However, all of them rely on “backpressure” quantities, which are the differential backlogs between neighboring nodes.
It has been observed in [9], [5], [7], [10] that the drift-plus-penalty and other alternative algorithms can be interpreted as
first order Lagrangian dual type methods for constrained optimization. In addition, these backpressure algorithms follow certain
fundamental utility-delay tradeoffs. For instance, the primal-dual type backpressure algorithm in [5] achieves an O(1/V ) utility
optimality gap with an O(V 2 ) queue length, where V is an algorithm parameter. By controlling parameter V , a small utility
optimality gap is available only at the cost of a large queue length. The drift-plus-penalty backpressure algorithm [4], which
has the best utility-delay tradeoff among all existing first order Lagrangian dual type methods for general networks, can only
achieve an O(1/V ) utility optimality gap with an O(V ) queue length. Under certain restrictive assumptions over the network,
a better [O(1/V ), O(log(V ))] tradeoff is achieved via an exponential Lyapunov function in [11], and an [O(1/V ), O(log2 (V ))]
tradeoff is achieved via a LIFO-backpressure algorithm in [12]. The existing utility-delay tradeoff results seem to suggest that
a large queueing delay is unavoidable if a small utility optimality gap is demanded.
Recently, there have been many attempts in obtaining new variations of backpressure algorithms by applying Newton’s
method to the Lagrangian dual function. In the recent work [10], the authors develop a Newton’s method for joint rate control
and routing. However, the utility-delay tradeoff in [10] is still [O(1/V ), O(V 2 )]; and the algorithm requires a centralized
projection step (although Newton directions can be approximated in a distributed manner). Work [13] considers a network
flow control problem where the path of each flow is given (and hence there is no routing part in the problem), and proposes
a decentralized Newton based algorithm for rate control. Work [14] considers network routing without an end-to-end utility
and only shows the stability of the proposed Newton based backpressure algorithm. All of the above Netwon’s method based
algorithms rely on distributed approximations for the inverse of Hessians, whose computations still require certain coordinations
for the local information updates and propagations and do not scale well with the network size. In contrast, the first order
Lagrangian dual type methods do not need global network topology information. Rather, each node only needs the queue
length information of its neighbors.
This paper proposes a new first order Lagrangian dual type backpressure algorithm that is as simple as the existing algorithms
in [4], [5], [7] but has a better utility-delay tradeoff. The new backpressue algorithm achieves a vanishing utility optimality gap
that decays like O(1/t), where t is the number of iterations. It also guarantees that the queue length at each node is always

2

bounded by a fixed constant of the same order as the optimal Lagrange multiplier of the network optimization problem. This
improves on the utility-delay tradeoffs of prior work. In particular, it improves the [O(1/V ), O(V 2 )] utility-delay tradeoff in
[5] and the [O(1/V ), O(V )] utility-delay tradeoff of the drift-plus-penalty algorithm in [4], both of which yield an unbounded
queue length to have a vanishing utility optimality gap. The new backpressure algorithm differs from existing first order
backpressure algorithms in the following aspects:
1) The “backpressure” quantities in this paper are with respect to newly introduced weights. These are different from queues
used in other backpressure algorithms, but can still be locally tracked and updated.
2) The rate control and routing decision rule involves a quadratic term that is similar to a term used in proximal algorithms
[15].
Note that the benefit of introducing a quadratic term in network optimization has been observed in [16]. Work [16] considers
a network utility maximization problem with given routing paths that is a special case of the problem treated in this paper.
The algorithm of [16] considers a fixed set of predetermined paths for each session and does not scale well when treating all
(typically exponentially many) possible paths of a general network. The algorithm proposed in [16] is not a backpressure type
and hence is fundamentally different from ours. For example, the algorithm in [16] needs to update the primal variables (source
session rates for each path) at least twice per iteration, while our algorithm only updates the primal variables (source session
rates and link session rates) once per iteration. The prior work [16] shows that the utility optimality gap is asymptotically zero
without analyzing the decay rate, while this paper shows the utility optimality gap decays like O(1/t).
II. S YSTEM M ODEL AND P ROBLEM F ORMULATION
Consider a slotted data network with normalized time slots t ∈ {0, 1, 2, . . .}. This network is represented by a graph
G = (N , L), where N is the set of nodes and L ⊆ N × N is the set of directed links. Let |N | = N and |L| = L. This network
is shared by F end-to-end sessions denoted by a set F. For each end-to-end session f ∈ F, the source node Src(f ) and
destination node Dst(f ) are given but the routes are not specified. Each session f has a continuous and concave utility function
Uf (xf ) that represents the “satisfaction” received by accepting xf amount of data for session f into the network at each
slot. Unlike [5], [10] where Uf (·) is assumed to be differentiable and strongly concave, this paper considers general concave
utility functions Uf (·), including those that are neither differentiable nor strongly concave. Formally, each utility function Uf
is defined over an interval dom(Uf ), called the domain of the function. It is assumed throughout that either dom(Uf ) = [0, ∞)
or dom(Uf ) = (0, ∞), the latter being important for proportionally fair utilities [17] Uf (x) = log(x) that have singularities
at x = 0 .
(f )
Denote the capacity of link l as Cl and assume it is a fixed and positive constant.1 Define µl as the amount of session
f ’s data routed at link l that is to be determined by our algorithm. Note that in general, the network may be configured such
that some session f is forbidden to use link l. For each link l, define Sl ⊆ F as the set of sessions that are allowed to use
link l. The case of unrestricted routing is treated by defining Sl = F for all links l.
(f )
(f )
Note that if l = (n, m) with n, m ∈ N , then µl and Cl can also be respectively written as µ(n,m) and C(n,m) . For each
node n ∈ N , denote the sets of its incoming links and outgoing links as I(n) and O(n), respectively. Note that xf , ∀f ∈ F
(f )
and µl , ∀l ∈ L, ∀f ∈ F are the decision variables of a joint rate control and routing algorithm. If the global network topology
information is available, the optimal joint rate control and routing can be formulated as the following multi-commodity network
flow problem:
X
max
Uf (xf )
(1)
(f )

xf ,µl

s.t.

f ∈F

xf 1{n=Src(f )} +

X

(f )

µl

l∈I(n)

X

(f )

µl

X

≤

(f )

µl , ∀f ∈ F, ∀n ∈ N \ {Dst(f )}

(2)

l∈O(n)

≤ Cl , ∀l ∈ L,

(3)

f ∈F
(f )

µl

≥ 0, ∀l ∈ L, ∀f ∈ Sl ,

(4)

(f )
µl

= 0, ∀l ∈ L, ∀f ∈ F \ Sl ,

(5)

xf ∈ dom(Uf ), ∀f ∈ F

(6)

where 1{·} is an indicator function; (2) represents the node flow conservation constraints relaxed by replacing the equality
with an inequality, meaning that the total rate of flow f into node n is less than or equal to the total rate of flow f out of the
node (since, in principle, we can always send fake data 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 ) since all incoming
data are consumed by this node.
1 As

stated in [10], this is a suitable model for wireline networks and wireless networks with fixed transmission power and orthogonal channels.

3

The above formulation includes network utility maximization with fixed paths as special cases. In the case when each
session only has one single given path, e.g., the network utility maximization problem considered in [18], we could modify
the sets Sl used in constraints (4) and (5) to reflect this fact. For example, if link l1 is only used for sessions f1 and f2 , then
Sl1 = {f1 , f2 }. Similarly, the case [16] where each flow is restricted to using links from a set of predefined paths can be
treated by modifying the sets Sl accordingly. See Appendix A for more discussions.
The solution to problem (1)-(6) corresponds to the optimal joint rate control and routing. However, to solve this convex
program at a single computer, we need to know the global network topology and the solution is a centralized one, which is
not practical for large data networks. As observed in [9], [5], [7], [10], various versions of backpressure algorithms can be
interpreted as distributed solutions to problem (1)-(6) from first order Lagrangian dual type methods.
(f ),∗
Assumption 1: (Feasibility) Problem (1)-(6) has at least one optimal solution vector [x∗f ; µl
]f ∈F ,l∈L .
Assumption 2: (Existence of Lagrange multipliers) Assume the convex program (1)-(6) has Lagrange multipliers attaining
(f )
the strong duality. Specifically, define convex set C = {[xf ; µl ]f ∈F ,l∈L : (3)-(6) hold}. Assume there exists a Lagrange
(f ),∗
∗
multiplier vector λ = [λn ]f ∈F ,n∈N \{Dst(f )} ≥ 0 such that
q(λ∗ ) = sup{(1) : (2)-(6)}
P
P
P
P
P
(f ) 
(f )
(f ) 
xf 1{n=Src(f )} + l∈I(n) µl − l∈O(n) µl
where q(λ) = sup[x ;µ(f ) ]∈C
f ∈F Uf (xf ) −
f ∈F
n∈N \{Dst(f )} λn
f
l
is the Lagrangian dual function of problem (1)-(6) by treating (3)-(6) as a convex set constraint.
Assumptions 1 and 2 hold in most cases of interest. For example, Slater’s condition guarantees Assumption 2. Since the
constraints (2)-(6) are linear, Proposition 6.4.2 in [19] ensures that Lagrange multipliers exist whenever constraints (2)-(6) are
feasible and when the utility functions Uf are either defined over open sets (such as Uf (x) = log(x) with dom(Uf ) = (0, ∞))
ef : (−, ∞) → R such
or can be concavely extended to open sets, meaning that there is an  > 0 and a concave function U
ef (x) = Uf (x) whenever x ≥ 0.2
that U
(f ),∗
Fact 1: (Replacing inequality with equality) If Assumption 1 holds, problem (1)-(6) has an optimal solution vector [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 in at 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 optimal solution vector such that at least one inequality constraint (2)
(f ),∗
is loose. Note that we can reduce the value of µl
on the right side of a loose (2) until either that constraint holds with
(f ),∗
equality, or until µl
reduces to 0. The objective function value does not change, and no constraints are violated. We can
repeat the process until all inequality constraints (2) are tight.
III. T HE N EW BACKPRESURE A LGORITHM
A. Discussion of Various Queueing Models
At each node, an independent queue backlog is maintained for each session. At each slot t, let xf [t] be the source session
(f )
(f )
rates; and let µl [t] be the link session rates. Some prior work enforces the constraint (2) via virtual queues Yn [t] of the
following form:
n
o
X (f )
X (f )
Yn(f ) [t + 1] = max Yn(f ) [t] + xf [t]1{n=Src(f )} +
µl [t] −
µl [t], 0 .
(7)
l∈I(n)

l∈O(n)

While this virtual equation is a meaningful approximation, it differs from reality in that new injected data are allowed to
be transmitted immediately, or equivalently, a single packet is allowed to enter and leave many nodes within the same slot.
(f )
Further, there is no clear connection between the virtual queues Yn [t] in (7) and the actual queues in the network. Indeed,
(f )
it is easy to construct examples that show there can be an arbitrarily large difference between the Yn [t] value in (7) and the
physical queue size in actual networks (see Appendix B).
(f )
An actual queueing network has queues Zn [t] with the following dynamics:
n
o
X (f )
X (f )
Zn(f ) [t + 1] ≤ max Zn(f ) [t] −
µl [t], 0 + xf [t]1{n=Src(f )} +
µl [t].
(8)
l∈O(n)

l∈I(n)

This is faithful to actual queue dynamics and does not allow data to be retransmitted over multiple hops in one slot. Note that
P
(f )
(8) is an inequality because the new arrivals from other nodes may be strictly less than l∈I(n) µl [t] because those other
nodes may not have enough backlog to send. The model (8) allows for any decisions to be made to fill the transmission values
P
(f )
(f )
(f )
µl [t] in the case that Zn [t] ≤ l∈O(n) µl [t], provided that (8) holds.
2 If dom(U ) = [0, ∞), such concave extension is possible if the right-derivative of U at x = 0 is finite (such as for U (x) = log(1 + x) or
f
f
f
√
Uf (x) = min[x, 3]). Such an extension is impossible for the example Uf (x) = x because the slope is infinite at x = 0. Nevertheless, Lagrange multipliers
often exist even for these utility functions, such as when Slater’s condition holds [19].

4

This paper develops an algorithm that converges to the optimal utility defined by problem (1)-(6), and that produces worstcase bounded queues on the actual queueing network, that is, with actual queues that evolve as given in (8). To begin, it is
convenient to introduce the following virtual queue equation
X (f )
X (f )
)
(f )
Q(f
µl [t] + xf [t]1{n=Src(f )} +
µl [t],
(9)
n [t + 1] =Qn [t] −
l∈O(n)

l∈I(n)

(f )

where Qn [t] represents a virtual queue value associated with session f at node n. At first glance, this model (9) appears to
(f )
be only an approximation, perhaps even a worse approximation than (7), because it allows the Qn [t] values to be negative.
(f )
Indeed, we use Qn [t] only as virtual queues to inform the algorithm and do not treat them as actual queues. However, this
paper shows that using these virtual queues to choose the µ[t] decisions ensures not only that the desired constraints (2) are
(f )
satisfied, but that the resulting µ[t] decisions create bounded queues Zn [t] in the actual network, where the actual queues
evolve according to (8). In short, our algorithm can be faithfully implemented with respect to actual queueing networks, and
converges to exact optimality on those networks.
(f )
The next lemma shows that if an algorithm can guarantee virtual queues Qn [t] defined in (9) are bounded, then actual
physical queues satisfying (8) are also bounded.
(f )
Lemma 1: Consider a network flow problem described by problem (1)-(6). For all l ∈ L and f ∈ F, let µl [t], xf [t]
(f )
(f )
(f )
be decisions yielded by a dynamic algorithm. Suppose Yn [t], Zn [t], Qn [t] evolve by (7)-(9) with initial conditions
(f )
(f )
(f )
(f )
Vn [0] = Zn [0] = Qn [0] = 0. If there exists a constant B > 0 such that |Qn [t]| ≤ B, ∀t, then
P
(f )
1) Zn [t] ≤ 2B + l∈O(n) Cl for all t ∈ {0, 1, 2, . . .}.
P
(f )
2) Yn [t] ≤ 2B + l∈O(n) Cl for all t ∈ {0, 1, 2, . . .}.
Proof:
P
)
b (f )
b (f
1) Fix f ∈ F, n ∈ N \ {Dst(f )}. Define an auxiliary virtual queue Q
n [t] that is initialized by Qn [0] = B +
l∈O(n) Cl
P
(f
)
(f
)
(f
)
b n [t] = Qn [t] + B +
and evolves according to (9). It follows that Q
l∈O(n) Cl , ∀t. Since Qn [t] ≥ −B, ∀t by
P
P
(f )
(f )
)
b
b (f
assumption, we have Qn [t] ≥
Cl ≥
µ [t], ∀t. This implies that Q
n [t] also satisfies:
l∈O(n)

l∈O(n)

l

n
o
X (f )
X (f )
)
(f )
b (f
b
Q
[t
+
1]
=
max
Q
[t]
−
µ
[t],
0
+ xf [t]1{n=Src(f )} +
µl [t], ∀t
n
n
l
l∈O(n)

(10)

l∈I(n)

)
(f )
b (f )
b (f
which is identical to (8) except the inequality is replaced by an equality. Since Zn [0] = 0 < Q
n [0]; and Qn [t]
)
(f )
b (f
satisfies (10), by inductions, Zn [t] ≤ Q
n [t], ∀t.
P
P
)
(f )
(f )
(f )
b (f
b
Since Qn [t] = Qn [t] + B + l∈O(n) Cl , ∀t and Qn [t] ≤ B, ∀t, we have Q
n [t] ≤ 2B +
l∈O(n) Cl , ∀t. It follows
P
(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 Algorithm
In this subsection, we propose a new backpressure algorithm that yields source session rates xf [t] and link session rates
at each slot such that the physical queues for each session at each node are bounded by a constant and the time average
utility satisfies
t−1
X
1XX
Uf (xf [t]) ≥
Uf (x∗f ) − O(1/t), ∀t
t τ =0
(f )
µl [t]

f ∈F

f ∈F

where x∗f are from the optimal solution to (1)-(6). Note that Jensen’s inequality further implies that
X
f ∈F

Uf

t−1
 X
1X
xf [τ ] ≥
Uf (x∗f ) − O(1/t), ∀t
t τ =0
f ∈F

The new backpressure algorithm is described in Algorithm 1. Similar to existing backpressure algorithms, the updates in
Algorithm 1 at each node n are fully distributed and only depend on weights at itself and its neighbor nodes. Unlike existing
(f )
(f )
backpressure algorithms, the weights used to update decision variables xf [t] and µl [t] are not the virtual queues Qn [t]
(f )
themselves, rather, they are augmented values Wn [t] equal to the sum of the virtual queues and the amount of net injected
data in the previous slot t − 1. In addition, the updates involve an additional quadratic term, which is similar to a term used
in proximal algorithms [15].

5

Algorithm 1 The New Backpressure Algorithm
(f )

(f )

Let αn > 0, ∀n ∈ N be constant parameters. Initialize xf [−1] = 0, µl [−1] = 0, ∀f ∈ F, ∀l ∈ L and Qn [0] = 0, ∀n ∈
N , ∀f ∈ F. At each time t ∈ {0, 1, 2, . . .}, each node n does the following:
(f )
• For each f ∈ F, if node n is not the destination node of session f , i.e., n 6= Dst(f ), then define weight Wn [t]:
X (f )
X (f )
)
Wn(f ) [t] =Q(f
µl [t − 1] −
µl [t − 1],
(11)
n [t] + xf [t − 1]1{n=Src(f )} +
l∈I(n)

l∈O(n)

(f )

•

If node n is the destination node, i.e., n = Dst(f ), then define Wn [t] = 0. Notify neighbor nodes (nodes k that can
(f )
send session f to node n, i.e., ∀k such that f ∈ S(k,n) ) about this new Wn [t] value.
For each f ∈ F, if node n is the source node of session f , i.e., n = Src(f ), choose xf [t] as the solution to
2
max Uf (xf ) − Wn(f ) [t]xf − αn xf − xf [t − 1]
(12)
xf

s.t.
•

xf ∈ dom(Uf )

(13)

(f )

For all (n, m) ∈ O(n), choose {µ(n,m) [t], ∀f ∈ F} as the solution to the following convex program:
X
 (f )
 X (f )
2
(f )
(f )
max
Wn(f ) [t] − Wm
[t] µ(n,m) − αn + αm
µ(n,m) − µ(n,m) [t − 1]

(14)

(f )

µ(n,m) f ∈F

s.t.

X

f ∈F
(f )

µ(n,m) ≤ C(n,m)

(15)

f ∈F
(f )

•

µ(n,m) ≥ 0, ∀f ∈ S(n,m)

(16)

(f )
µ(n,m)

(17)

= 0, ∀f 6∈ S(n,m)
(f )

For each f ∈ F, if node n is not the destination of f , i.e., n 6= Dst(f ), update virtual queue Qn [t + 1] by (9).

C. Almost Closed-Form Updates in Algorithm 1
(f )

This subsection shows the decisions xf [t] and µl [t] in Algorithm 1 have either closed-form solutions or “almost” closedform solutions at each iteration t.
Lemma 2: Let x̂f ≡ xf [t] denote the solution to (12)-(13).
(f )

1) Suppose dom(Uf ) = [0, ∞) and Uf (xf ) is differentiable. Let h(xf ) = Uf0 (xf ) − 2αn xf + 2αn xf [t − 1] − Wn [t]. If
h(0) < 0, then x̂f = 0; otherwise x̂f is the root to the equation h(xf ) = 0 and can be found by a bisection search.
2) Suppose dom(Uf ) = (0, ∞) and Uf (xf ) = wf log(xf ) for some weight wf > 0. Then:
q
(f )
(f )
(Wn [t] − 2αn xf [t − 1])2 + 8αn wf
2αn xf [t − 1] − Wn [t]
x̂f =
+
4αn
4αn
Proof: Omitted for brevity.
(f )
The problem (14)-(17) can be represented as follows by eliminating µ(n,m) , f 6∈ S(n,m) , completing the square and replacing
maximization with minimization. (Note that K = |S(n,m) | ≤ |F|.)
K

min

s.t.

1X
(zk − ak )2
2

(18)

zk ≤ b

(19)

k=1
K
X

k=1

zk ≥ 0, ∀k ∈ {1, 2, . . . , K}

(20)

Lemma 3: The solution to problem (18)-(20) is given by zk∗ = max{0, ak − θ∗ }, ∀k ∈ {1, 2, . . . , K} where θ∗ ≥ 0 can be
found either by a bisection search (See Appendix D) or by Algorithm 2 with complexity O(K log K).
Proof: A similar problem where (19) is replaced with an equality constraint in considered in [20]. The optimal solution
to this quadratic program is characterized by its KKT condition and a corresponding algorithm can be developed to obtain its
KKT point. A complete proof is presented in Appendix D.
Note that step (3) in Algorithm 2 has complexity O(K) and hence the overall complexity of Algorithm 2 is dominated by
the sorting step (2) with complexity O(K log(K)).

6

Algorithm 2 Algorithm to solve problem (18)-(20)
PK
1) Check if k=1 max{0, ak } ≤ b holds. If yes, let θ∗ = 0 and zk∗ = max{0, ak }, ∀k ∈ {1, 2, . . . , K} and terminate the
algorithm; else, continue to the next step.
2) Sort all ak , ∈ {1, 2, . . . , K} in a decreasing order π such that aπ(1) ≥ aπ(2) ≥ · · · ≥ aπ(K) . Define S0 = 0.
3) For k = 1 to K
S −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 next iteration in the
loop.
4) Let zk∗ = max{0, ak − θ∗ }, ∀k ∈ {1, 2, . . . , K} and terminate the algorithm.

IV. P ERFORMANCE A NALYSIS OF A LGORITHM 1
A. Basic Facts from Convex Analysis
Definition 1 (Lipschitz Continuity): Let Z ⊆ Rn be a convex set. Function h : Z → Rm is said to be Lipschitz continuous
on Z with modulus β if there exists β > 0 such that kh(z1 ) − h(z2 )k ≤ βkz1 − z2 k for all z1 , z2 ∈ Z.
Definition 2 (Strongly Concave Functions): Let Z ⊆ Rn be a convex set. Function h is said to be strongly concave on Z
with modulus α if there exists a constant α > 0 such that h(z) + 21 αkzk2 is concave on Z.
By the definition of strongly concave functions, it is easy to show that if h(z) is concave and α > 0, then h(z) − αkz − z0 k2
is strongly concave with modulus 2α for any constant z0 .
Lemma 4: Let Z ⊆ Rn be a convex set. Let function h be strongly concave on Z with modulus α and zopt be a global
maximum of h on Z. Then, h(zopt ) ≥ h(z) + α2 kzopt − zk2 for all z ∈ Z.
B. Preliminaries
(f )

Define column vector y = [xf ; µl ]f ∈F ,l∈L . For each f ∈ F, n ∈ N \ {Dst(f )}, define column vector
(
(f )
[xf ; µl ]l∈I(n)∪O(n) if n = Src(f ),
(f )
yn =
(f )
[µl ]l∈I(n)∪O(n)
else,

(21)
(f )

which is composed by the control actions appearing in each constraint (2); and introduce a function with respect to yn as
X (f )
X (f )
gn(f ) (yn(f ) ) = xf 1{n=Src(f )} +
µl −
µl
(22)
l∈I(n)

l∈O(n)

Thus, constraint (2) can be rewritten as
gn(f ) (yn(f ) ) ≤ 0, ∀f ∈ F, ∀n ∈ N \ {Dst(f )}.
(f )

Note that each vector yn is a subvector of y and has length dn + 1 where dn is the degree of node n (the total number
of outgoing links and incoming links) if node n is the source of session f ; and has length dn if node n is not the source of
session f .
(f )
(f )
Fact 2: Each function gn (·) defined in (22) is Lipschitz continuous with respect to vector yn with modulus
p
βn ≤ dn + 1.
where dn is the degree of node n.
(f )
(f )
(f )
Proof: This fact can be easily shown by noting that each gn (yn ) is a linear function with respect to vector yn and
has at most dn + 1 non-zero coefficients that are equal to ±1.
Note that virtual queue update equation (9) can be rewritten as:
)
(f )
(f )
(f )
Q(f
n [t + 1] = Qn [t] + gn (yn [t]),

(23)

and weight update equation (11) can be rewritten as:
)
(f )
(f )
Wn(f ) [t] = Q(f
n [t] + gn (yn [t − 1]).

(24)

Define
L(t) =

1X
2

X

f ∈F n∈N \Dst(f )

2
)
Q(f
n [t]

(25)

7

and call it a Lyapunov function. In the remainder of this paper, double summations are often written compactly as a single
summation, e.g.,
X
X
X
∆

· =
· .
f ∈F n∈N \Dst(f )

f ∈F ,
n∈N \Dst(f )

Define the Lyapunov drift as
∆[t] = L(t + 1) − L(t).
The following lemma follows directly from equation (23).
Lemma 5: At each iteration t ∈ {0, 1, . . .} in Algorithm 1, the Lyapunov drift is given by

X
1 (f ) f 2 
)
(f )
f
Q(f
g (yn [t])
∆[t] =
.
n [t]gn (yn [t]) +
2 n

(26)

f ∈F ,
n∈N \Dst(f )

Proof: Fix f ∈ F and n ∈ N \ Dst(f ), we have
2 1
2
1 (f )
)
Qn [t + 1] − Q(f
n [t]
2
2
2 1
2
(a) 1
(f )
(f )
(f )
)
= Qn [t] + gn (yn [t]) − Q(f
n [t]
2
2
1 (f ) f 2
(f )
(f )
f
=Qn [t]gn (yn [t]) + gn (yn [t])
2

(27)

where (a) follows from (23).
By the definition of ∆[t], we have
∆[t] =
(a)

=

1
2



X

)
Q(f
n [t + 1]

2

2 
)
− Q(f
n [t]

f ∈F ,
n∈N \Dst(f )

X
f ∈F ,
n∈N \Dst(f )



)
(f )
f
Q(f
n [t]gn (yn [t]) +

1 (f ) f 2 
g (yn [t])
2 n

where (a) follows from
P (27).
Define f (y) = f ∈F Uf (xf ). At each time t, consider choosing a decision vector y[t] that includes elements in each
(f )
subvector yn [t] to solve the following problem:
X
X
X (f )

(f )
max f (y) −
Wn(f ) [t]gn(f ) (yn(f ) ) + αn kyn(f ) − yn(f ) [t − 1]k2 −
αn
(µl − µl [t − 1])2
(28)
y

f ∈F ,
n∈N \Dst(f )

f ∈F ,
n=Dst(f )

l∈I(n)

s.t. (3)-(6)

(29)

The expression (28) is a modified drift-plus-penalty expression. Unlike the standard drift-plus-penalty expressions from [4],
(f )
(f )
(f )
(f )
the above expression uses weights Wn [t], which arguments each Qn [t] by gn (yn [t − 1]), rather than virtual queues
(f )
Qn [t]. It also includes a “prox”-like term that penalizes deviation from the previous y[t − 1] vector. This results in the novel
backpressure-type algorithm of Algorithm 1. Indeed, the decisions in Algorithm 1 were derived as the solution to the above
problem (28)-(29). This is formalized in the next lemma.
Lemma 6: At each iteration t ∈ {0, 1, . . .}, the action y[t] jointly chosen in Algorithm 1 is the solution to problem (28)-(29).
(f )
Proof: The proof involves collecting terms associated with the xf [t] and µl [t] decisions. See Appendix E for details.
Furthermore, the next lemma summarizes that the action y[t] jointly chosen in Algorithm 1 provides a lower bound for the
drift-plus-penalty expression at each iteration t ∈ {0, 1, . . .}.
(f ),∗
(f )
(f ),∗
Lemma 7: Let y∗ = [x∗f ; µl
]f ∈F ,l∈L be an optimal solution to problem (1)-(6) given in Fact 1, i.e., gn (yn ) = 0, ∀f ∈
(f )
F, ∀n ∈ N \Dst(f ). If αn ≥ 12 (dn +1), ∀n ∈ N , where dn is the degree of node n, then the action y[t] = [xf [t]; µl [t]]f ∈F ,l∈L
jointly chosen in Algorithm 1 at each iteration t ∈ {0, 1, . . .} satisfies
−∆[t] + f (y[t]) ≥ f (y∗ ) + Φ[t] − Φ[t − 1]

P
P
(f ),∗
(f )
(f ),∗
(f )
where Φ[t] = f ∈F ,n∈N αn 1{n6=Dst(f )} kyn
− yn [t]k2 + αn 1{n=Dst(f )} l∈I(n) (µl
− µl [t])2 .
Proof: See Appendix F.
It remains to show that this modified backpressure algorithm leads to fundamentally improved performance.

8

C. Utility Optimality Gap Analysis
 (f ) 
(f )
Define column vector Q[t] = Qn [t] f ∈F ,n∈N \{Dst(f )} as the stacked vector of all virtual queues Qn [t] defined in (9).
(f )

(f )

Note that (25) can be rewritten as L(t) = 21 kQ[t]k2 . Define vectorized constraints (2) as g(y) = [gn (yn )]f ∈F ,n∈N \Dst(f ) .
(f ),∗
(f )
(f ),∗
Lemma 8: Let y∗ = [x∗f ; µl
]f ∈F ,l∈L be an optimal solution to problem (1)-(6) given in Fact 1, i.e., gn (yn ) =
1
0, ∀f ∈ F, ∀n ∈ N \ Dst(f ). If αn ≥ 2 (dn + 1), ∀n ∈ N in Algorithm 1, where dn is the degree of node n, then for all t ≥ 1,
t−1
X

1
f (y[τ ]) ≥ tf (y∗ ) − ζ + kQ[t]k2 .
2
τ =0
P
P
(f ),∗
(f ),∗ 2 
where ζ = Φ[−1] = f ∈F ,n∈N αn 1{n6=Dst(f )} kyn k2 + αn 1{n=Dst(f )} l∈I(n) (µl
) is a constant.
Proof: By Lemma 7, we have −∆[τ ] + f (y[τ ]) ≥ f (y∗ ) + Φ[t] − Φ[t − 1], ∀τ ∈ {0, 1, . . . , t − 1}. Summing over
τ ∈ {0, 1, . . . , t − 1} yields
t−1
X

f (y[τ ]) −

τ =0

t−1
X

∆[τ ]

τ =0

≥tf (y∗ ) +

t−1
X

Φ[τ ] − Φ[τ − 1]



τ =0

=tf (y∗ ) + Φ[t] − Φ[−1]
(a)

≥ tf (y∗ ) − Φ[−1]

where (a) follows from the fact that Φ[t] ≥ 0, ∀t.
Recall ∆[τ ] = L[τ + 1] − L[τ ], simplifying summations and rearranging terms yields
t−1
X

f (y[τ ]) ≥tf (y∗ ) − Φ[−1] + L[t] − L[0]

τ =0

1
= tf (y∗ ) − Φ[−1] + kQ[t]k2
2

(a)

where (a) follows from the fact that L[0] = 0 and L[t] = 12 kQ[t]k2 .
The next theorem summarizes that Algorithm 1 yields a vanishing utility optimality gap that approaches zero like O(1/t).
(f ),∗
(f ),∗
(f )
Theorem 1: Let y∗ = [x∗f ; µl
]f ∈F ,l∈L be an optimal solution to problem (1)-(6) given in Fact 1, i.e., gn (yn ) =
1
0, ∀f ∈ F, ∀n ∈ N \ Dst(f ). If αn ≥ 2 (dn + 1), ∀n ∈ N in Algorithm 1, where dn is the degree of node n, then for all
t ≥ 1, we have
t−1
X
1XX
1
Uf (xf [τ ]) ≥
Uf (x∗f ) − ζ,
t τ =0
t
f ∈F

f ∈F

where ζ is a constant defined in Lemma 8. Moreover, if we define xf [t] =
X

Uf (xf [t]) ≥

f ∈F

Proof: Recall that f (y) =

P

f ∈F

X
f ∈F

1
t

Pt−1

τ =0

xf [τ ], ∀f ∈ F, then

1
Uf (x∗f ) − ζ.
t

Uf (xf ). By Lemma 8, we have

t−1 X
X

Uf (xf [τ ]) ≥t

τ =0 f ∈F

X
f ∈F

(a)

≥t

X

1
Uf (x∗f ) − ζ + kQ[t]k2
2
Uf (x∗f ) − ζ.

f ∈F

where (a) follows from the trivial fact that kQ[t]k2 ≥ 0.
Dividing both sides by a factor t yields the first inequality in this theorem. The second inequality follows from the concavity
of Uf (·) and Jensen’s inequality.

9

D. Queue Stability Analysis
Lemma 9: Let Q[t], t ∈ {0, 1, . . .} be the virtual queues in Algorithm 1. For any t ≥ 1,
Q[t] =

t−1
X

g(y[τ ])

τ =0

Proof: This lemma follows directly from the fact that Q[0] = 0 and queue update equation (9) can be written as
Q[t + 1] = Q[t] + g(y[t − 1]).
(f )
The next theorem shows the boundedness of all virtual queues Qn [t] in Algorithm 1.
(f ),∗
(f )
(f ),∗
∗
∗
Theorem 2: Let y = [xf ; µl
]f ∈F ,l∈L be an optimal solution to problem (1)-(6) given in Fact 1, i.e., gn (yn ) =
∗
0, ∀f ∈ F, ∀n ∈ N \ Dst(f ), and λ be a Lagrange multiplier vector given in Assumption 2. If αn ≥ 21 (dn + 1)2 , ∀n ∈ N in
Algorithm 1, where dn is the degree of node n, then for all t ≥ 1,
p
∗
)
|Q(f
2ζ, ∀f ∈ F, ∀n ∈ N \ {Dst(f )}.
n [t]| ≤ 2kλ k +
where ζ is a constant defined in Lemma
8.

Proof: Let q(λ) = supy∈C f (y) − λT g(y)
{0, 1, . . . , }, by Assumption 2, we have

be the Lagrangian dual function defined in Assumption 2. For all τ ∈
(a)

f (y∗ ) = q(λ∗ ) ≥ f (y[τ ]) − λ∗,T g(y[τ ])
where (a) follows from the definition of q(λ∗ ). Rearranging terms yields
f (y[τ ]) ≤ f (y∗ ) + λ∗,T g(y[τ ]), ∀τ ∈ {0, 1, . . .}.
Fix t > 0. Summing over τ ∈ {0, 1, . . . , t − 1} yields
t−1
X

f (y[τ ]) ≤tf (y∗ ) +

τ =0

t−1
X

λ∗,T g(y[τ ])

τ =0

=tf (y∗ ) + λ∗,T

t−1
X

g(y[τ ])

τ =0
(a)

= tf (y∗ ) + λ∗,T Q[t]

(b)

≤tf (y∗ ) + kλ∗ kkQ[t]k

where (a) follows form Lemma 9 and (b) follows from Cauchy-Schwarz inequality.
On the other hand, by Lemma 8, we have
t−1
X

1
f (y[τ ]) ≥ tf (y∗ ) − ζ + kQ[t]k2 .
2
τ =0
Combining the last two inequalities and cancelling the common terms yields
1
kQ[t]k2 − ζ ≤ kλ∗ kkQ[t]k
2
2
⇒ kQ[t]k − kλ∗ k ≤ kλ∗ k2 + 2ζ
q
⇒kQ[t]k ≤ kλ∗ k + kλ∗ k2 + 2ζ
p
(a)
⇒kQ[t]k ≤ 2kλ∗ k + 2ζ
√
√
√
where (a) follows from the basic inequality a + b ≤ a + b for any a, b ≥ 0.
Thus, for any f ∈ F and n ∈ N \ {Dst(f )}, we have
p
∗
)
|Q(f
2ζ.
n [t]| ≤ kQ[t]k ≤ 2kλ k +
√
(f )
This theorem shows that the absolute values of all virtual queues Qn [t] are bounded by a constant B = 2kλ∗ k + 2ζ from
(f )
(f )
above.P
By Lemma 1 and discussions in Section III-A, the actual physical queues Zn [t] evolving via (8) satisfy Zn [t] ≤
2B + l∈O(n) Cl , ∀t. This is summarized in the next corollary.
(f ),∗

(f )

(f ),∗

Corollary 1: Let y∗ = [x∗f ; µl
]f ∈F ,l∈L be an optimal solution to problem (1)-(6) given in Fact 1, i.e., gn (yn ) =
0, ∀f ∈ F, ∀n ∈ N \ Dst(f ), and λ∗ be a Lagrange multiplier vector given in Assumption 2. If αn ≥ 21 (dn + 1)2 , ∀n ∈ N in

10

(f )

Algorithm 1, where dn is the degree of node n, then all actual physical queues Zn [t], ∀f ∈ F, ∀n ∈ N \ {Dst(f )} in the
network evolving via (8) satisfy
X
p
Cl , ∀t.
Zn(f ) [t] ≤4kλ∗ k + 2 2ζ +
l∈O(n)

where ζ is a constant defined in Lemma 8.
E. Performance of Algorithm 1
Theorems 1 and 2 together imply that Algorithm 1 with αn ≥ 12 (dn + 1), ∀n ∈ N can achieve a vanishing utility optimality
gap that decays like O(1/t), where t is number of iterations, and guarantees the physical queues at each node for each session
are always bounded by a constant that is independent of the utility optimality gap.
This is superior to existing backpressure algorithms from [5], [4], [10] that can achieve an O(1/V ) utility gap only at
the cost of an O(V 2 ) or O(V ) queue length, where V is an algorithm parameter. To obtain a vanishing utility gap, existing
backpressure algorithms in [5], [4], [10] necessarily yield unbounded queues. To obtain a vanishing utility gap, existing
backpressure algorithms in [5], [4] yield unbounded queues. We also comment that O(V 2 ) queue bound in the primal-dual
type backpressure algorithm [5] is actually of the order V 2 kλ∗ k + B1 where λ∗ is the Lagrangian multiplier vector attaining
strong duality and B1 is a constant determined by the problem parameters. A recent work [21] also shows that the O(V ) queue
bound in the backpressure algorithm from drift-plus-penalty is of the order V kλ∗ k+B2 where B2 is also a constant determined
by the problem parameters. Since λ∗ is a constant vector independent of V , both algorithms are claimed to have O(V 2 ) or
O(V ) queue bounds. By Corollary 1, Algorithm 1 guarantees physical queues at each node are bounded by 4kλ∗ k + B3 ,
where B3 is constant given a problem. Thus, the constant queue bound guaranteed by Algorithm 1 is typically smaller than
the O(V 2 ) or O(V ) queue bounds from [5] and [21] even for a small V . (A small V can yield a poor utility performance in
the backpressure algorithms in [5], [4].)
V. N UMERICAL E XPERIMENT
In this section, we consider a simple network with 6 nodes and 8 links and 2 sessions as described in Figure 1. This network
has two sessions: session 1 from node 1 to node 6 has utility function log(x1 ) and session 2 from node 3 to node 4 has utility
function 1.5 log(x2 ). (The log utilities are widely used as metrics of proportional fairness in the network [17].) The routing
path of each session is arbitrary as long as data can be delivered from the source node to the destination node. For simplicity,
assume that each link has capacity 1. The optimal source session rate to problem (1)-(6) is x∗1 = 1.2 and x∗2 = 1.8 and link
session rates, i.e., static routing for each session, is drawn in Figure 2.

2

l3

4
l7

l1

11

l4

l2
Session 1: 1->6

3

l5

6

l6

5

l8

Session 2: 3->4
Fig. 1.

A simple network with 6 nodes, 8 links and 2 sessions.

To compare the convergence performance of Algorithm 1 and the backpressure algorithm in [4] (with the best utility-delay
tradeoff among all existing backpressure algorithms), we run both Algorithm 1 with αn = 12 dn + 1), ∀n ∈ N and the
backpressure algorithm in [4] with V = 500 to plot Figure 3. It can be observed from Figure 3 that Algorithm 1 converges to
the optimal source session rates faster than the backpressure algorithm in [4]. The backpressure algorithm in [4] with V = 400
takes around 2500 iterations to converges to source rates close to (1.2, 1.8) while Algorithm 1 only takes around 800 iterations
to converges to (1.2, 1.8) (as shown in the zoom-in subfigure at the top right corner.) In fact, the backpressure algorithm in [4]
with V = 500 can not converge to the exact optimal source session rate (1.2, 1.8) but can only converge to its neighborhood
with a distance gap determined by the value of V . This is an effect from the fundamental [O(1/V ), O(V )] utility-delay tradeoff

11

2

1

4

11

6

0.2

Session 1: 1->6

4

1

1

0.8

1

0.2

3

0.2

5

0.8

3

5

Session 2: 3->4
Fig. 2.

The optimal routing for the network in Figure 1.

of the the backpressure algorithm in [4]. In contrast, Algorithm 1 can eventually converge to the the exact optimal source
session rate (1.2, 1.8). A zoom-in subfigure at the bottom right corner in Figure 1 verifies this and shows that the source
rate for Session 1 in Algorithm 1 converges to 1.2 while the source rate in the backpressure algorithm in [4] with V = 500
oscillates around a point slightly larger than 1.2.
5
4.5

Backpressure Alg in [4]

New Backpressure Alg: Algorithm 1

2
4

1.5
1

3.5

Source Rates

0.5
3

0

Session 2: 3->4

0

200

400

600

800

2.5
2
1.5
1
1.205
1.2
1.195
3000

0.5
Session 1: 1->6
0

0

500

1000

1500

2000

2500

3500

3000

4000
3500

4000

Iterations
Fig. 3.

Convergence performance comparison between Algorithm 1 and the backpressure algorithm in [4].

√
∗
∗
PCorollary 1 shows that Algorithm 1 guarantees each actual queue in the network is bounded by constant 4kλ k+2 2ξky k+
l∈O(n) Cl . Recall that the backpressure algorithm in [4] can guarantee the actual queues in the network are bounded by a
constant of order V kλ∗ k. Figure 4 plots the sum of actual queue length at each node for Algorithm 1 and the backpressure
algorithm in [4] with V = 10, 100 and 500. (Recall a larger V in the backpressure algorithm in [4] yields a smaller utility
gap but a larger queue length.) It can be observed that Algorithm 1 has the smallest actual queue length (see the zoom-in
subfigure) and the actual queue length of the backpressure algorithm in [4] scales linearly with respect to V .
VI. C ONCLUSION
This paper develops a new first-order Lagrangian dual type backpressure algorithm for joint rate control and routing in multihop data networks. The new backpressure algorithm can achieve vanishing utility optimality gaps and finite queue lengths.
This improves the state-of-art [O(1/V ), O(V 2 )] or [O(1/V ), O(V )] utility-delay tradeoff attained by existing backpressure
algorithms [5], [9], [7], [10].

12

1400

Network Sum Physical Queue Length

1200

Backpressure Alg in [4] (V=500)

1000

800
40
30

600

20
10

400

0

0

50

100

150

200

250

300

Backpressure Alg in [4] (V=100)

200

Backpressure Alg in [4] (V=10)
0

New Backpressure Alg: Algorithm 1
0

500

1000

1500

2000

2500

3000

Iterations
Fig. 4.

Actual Queue length comparison between Algorithm 1 and the backpressure algorithm in [4].

A PPENDIX A
N ETWORK U TILITY M AXIMIZATION WITH P REDETERMINED M ULTI -PATH
Consider multi-path network utility maximization in [16] where each session has multiple given paths, then the source
session rate xf in problem (1)-(6) becomes a vector xf = [xf,j ]j∈Pf where Pf is the set of paths for session f and the
(f )
(f )
(f,j)
(f )
link session rate µl becomes a vector µl = [µl ]j∈Pl . Define Sl as the set of paths for session f that are allowed
(f )
to use link l. Note that if all paths of session f are forbidden to use link l, then Sl = ∅. The multi-path network utility
maximization problem can be formulated as follows:
X
X
max
Uf (
xf,j )
f ∈F

s.t.

X

j∈Pf

xf,j 1{n=Src(f )} +

j∈Pl

X X

(f )

µl

j∈Pl l∈I(n)

X X

(f,j)
µl

≤

X X

(f )

µl , ∀f ∈ F, ∀n ∈ N \ {Dst(f )}

j∈Pl l∈O(n)

≤ Cl , ∀l ∈ L,

f ∈F j∈Pl
(f,j)

µl

(f )

≥ 0, ∀l ∈ L, ∀f ∈ F, ∀j ∈ Sl ,

(f,j)

(f )

µl
= 0, ∀l ∈ L, ∀f ∈ F, ∀j ∈ Pf \ Sl ,
X
xf,j ∈ dom(Uf ), ∀f ∈ F,
j∈Pl

xf,j ≥ 0, ∀f ∈ F, j ∈ Pf
The above formulation is in the form of problem (1)-(6) except that the variable dimension is extended.
A PPENDIX B
A N E XAMPLE I LLUSTRATING THE P OSSIBLY L ARGE G AP BETWEEN M ODEL (7) AND M ODEL (8)
Consider a network example shown in Figure 5. The network has 3k + 1 nodes where only node 0 is a destination; and
ai , i ∈ {1, 2, . . . , k} and bi , i ∈ {1, 2, . . . , k} can have exogenous arrivals. Assume all link capacities are equal to 1; and the
exogenous arrivals are periodic with period 2k, as follows:
• Time slot 1: One packet arrives at node a1 .
• Time slot 2: One packet arrives at node a2 .

13

···
Time slot k: One packet arrives at node ak .
• Time slot k + 1: One packet arrives at node b1 .
• Time slot k + 2: One packet arrives at node b2 .
• ···
• Time slot 2k: One packet arrives at node bk .
Under dynamics (7), each packet arrives on its own slot and traverses all links of its path to exit on the same slot it arrived.
The queue backlog in each node is 0 for all time.
Under dynamics (8), the first packet arrives at time slot 1 to node a1 . This packet visits node a2 at time slot 2, when the
second packet also arrives at a2 . One of these packets is delivered to node a3 at time slot 3, and another packet also arrives to
node 3. The nodes {1, . . . , k} do not have any exogenous arrivals and act only to delay the delivery of all packets from the ai
nodes. It follows that the link from node k to node 0 will send exactly one packet over each slots t ∈ {2k+1, 2k+2, . . . , 2k+k}.
Similarly, the link from bk to 0 sends exactly one packet to node 0 over each of these same slots. Thus, node 0 receives 2
packets on each slot t ∈ {2k + 1, 2k + 2, . . . , 2k + k}, but can only output 1 packet per slot. The queue backlog in this node
grows linearly and reaches k + 1 at time 2k + k. Thus, the backlog in node 0 can be arbitrarily large when k is large. This
example demonstrates that, even when there is only one destination, the deviation between virtual queues under dynamics (7)
and actual queues under dynamics (8) can be arbitrarily large, even when the in-degree and out-degree of 1 and an in-degree
of at most 2.
•
•

a1

ak-1

a2

ak
1
2
k

b1
Fig. 5.

b2

bk-1

bk

0

An example illustrating the possibly large gap between queue model (7) and queue model (8)

A PPENDIX C
P ROOF OF PART (2) IN L EMMA 1
Fix f ∈ F, n ∈ N \ {Dst(f )}. By (10),
)
b (f
Q
n [t + 1]
X (f )
X (f )
 (f )
b n [t] −
= max Q
µl [t], 0 + xf [t]1{n=Src(f )} +
µl [t]
l∈O(n)

l∈I(n)

n
X (f )
X (f )
X (f ) o
)
b (f
= max Q
µl [t] −
µl [t], xf [t]1{n=Src(f )} +
µl [t]
n [t] + xf [t]1{n=Src(f )} +
l∈I(n)

l∈O(n)

l∈I(n)

(a)

n
o
X (f )
X (f )
)
b (f
≥ max Q
µl [t] −
µl [t], 0
n [t] + xf [t]1{n=Src(f )} +
l∈I(n)

l∈O(n)

(f )
where (a) follows from the fact that µl [t], xf [t], ∀f, l, t are non-negative.
(f )
)
b (f
is identical to the right side of (7) and recall that Yn [0] = 0 < Q
n [0].
P
(f
)
(f
)
(f
)
b n [t] = Qn [t] + B +
Since Q
l∈O(n) Cl , ∀t and Qn [t] ≤ B, ∀t, we have
P
(f )
Yn [t] ≤ 2B + l∈O(n) Cl , ∀t.

Note that the right side of the above equation
(f )
)
b (f
By inductions, we have Yn [t] ≤ Q
n [t], ∀t.
P
(f
)
b n [t] ≤ 2B +
Q
l∈O(n) Cl , ∀t. It follows that

A PPENDIX D
P ROOF OF L EMMA 3
Note that problem (18)-(20) satisfies Slater’s condition. So the optimal solution toP
problem (18)-(20) is characterized by KKT
K
conditions [22]. Introducing Lagrange multipliers θ ∈ R+ for inequality constraint k=1 zk ≤ b and ν = [ν1 , . . . , νK ]T ∈ RK
+
∗
∗
∗ T
for inequality constraints zk ≥ 0, k ∈ {1, 2, . . . , K}. Let z = [z1 , . . . , zK ] and (θ∗ , ν ∗ ) be any primal and dual pair
PK
with the zero duality
gap. By KKT conditions, we have zk∗ − ak + θ∗ − νk∗ = 0, ∀k ∈ {1, 2, . . . , K}; k=1 zk∗ ≤ b; θ∗ ≥

P
K
∗
∗
∗
∗ ∗
0; θ∗
k=1 zk − b = 0; zk ≥ 0, ∀k ∈ {1, 2, . . . , K}; νk ≥ 0, ∀k ∈ {1, 2, . . . , K}; νk zk = 0, ∀k ∈ {1, 2, . . . , K}.

14

PK
PK
∗
νk∗ , ∀k ∈ {1, 2, . . . , K} in all equations yields θ∗ ≥ ak −zk∗ , k ∈ {1, 2, . . . , K}; k=1 zk∗ ≤ b; θ∗ ≥ 0; θ∗
k=1 zk −
 Eliminating
∗
∗
∗ ∗
b = 0; zk ≥ 0, ∀k ∈ {1, 2, . . . , K}; (zk − ak + θ )zk = 0, ∀k ∈ {1, 2, . . . , K}.
For all k ∈ {1, 2, . . . , K}, we consider θ∗ < ak and θ∗ ≥ ak separately:
1) If θ∗ < ak , then θ∗ ≥ ak − zk∗ holds only when zk∗ > 0, which by (zk∗ − ak + θ∗ )zk∗ = 0 implies that zk∗ = ak − θ∗ .
2) If θ∗ ≥ ak , then zk∗ > 0 is impossible, because zk∗ > 0 implies that zk∗ − ak + θ∗ > 0, which together with zk∗ > 0
contradicts the slackness condition (zk∗ − ak + θ∗ )zk∗ = 0. Thus, if θ∗ ≥ ak , we must have zk∗ = 0.
PK
Summarizing both cases, wehave zk∗ = max{0, ak − θ∗ }, ∀k ∈ {1, 2, . . . , K}, where θ∗ is chosen such that k=1 zk∗ ≤ b,
P
K
∗
θ∗ ≥ 0 and θ∗
k=1 zk − b = 0.

PK
∗
zk∗ − b is guaranteed to hold
To find such θ , we first check if θ∗ = 0. If θ∗ = 0 is true, the slackness condition θ∗
k=1P
PK
P
K
K
and we need to further require k=1 zk∗ = k=1 max{0, ak } ≤ b. Thus θ∗ = 0 if and only if k=1 max{0, ak } ≤ b. Thus,
PK
Algorithm 2 check if k=1 max{0, ak } ≤ b holds at the first step and if this is true, then we conclude θ∗ = 0 and we are
done!

PK
PK
PK
Otherwise, we know θ∗ > 0. By the slackness condition θ∗
zk∗ −b = 0, we must have k=1 zk∗ = k=1 max{0, ak −
k=1
PK
θ∗ } = b. To find θ∗ > 0 such that k=1 max{0, ak − θ∗ } = b, we could apply a bisection search by noting that all zk∗ are
decreasing with respect to θ∗ .
Another algorithm of finding θ∗ is inspired by the observation that if aj ≥ ai , ∀i, j ∈ {1, 2, . . . , K}, then zj∗ ≥ zi∗ . Thus, we
first sort all ak in a decreasing order, say π is the permutation such that aπ(1) ≥ aπ(2) ≥ · · · ≥ aπ(K) ; and then sequentially
check if k ∈ {1, 2, . . . , K} is the index such that aπ(k) −θ∗ ≥ 0 and aπ(k+1) −θ∗ < 0. To check this, we first assume k is indeed
Pk
such an index and solve the equation j=1 (aπ(j) − θ∗ ) = b to obtain θ∗ ; (Note that in Algorithm 2, to avoid recalculating the
Pk
Pk
partial sum j=1 aπ(j) for each k, we introduce the parameter Sk = j=1 aπ(j) and update Sk incrementally. By doing this,
the complexity of each iteration in the loop is only O(1).) then verify the assumption by checking if θ∗ ≥ 0, aπ(k) − θ∗ ≥ 0
and aπ(k+1) − θ∗ ≤ 0. The algorithm is described in Algorithm 2 and has complexity O(K log(K)). The overall complexity
is dominated by the step of sorting all ak .
A PPENDIX E
P ROOF OF L EMMA 6
The objective function (28) can be rewritten as
X

f (y) −
Wn(f ) [t]gn(f ) (yn(f ) ) + αn kyn(f ) − yn(f ) [t − 1]k2 −
f ∈F ,
n∈N \Dst(f )
(a)

=

X
f ∈F

−

l∈I(n)



f ∈F ,
n∈N \Dst(f )

−

f ∈F ,n=Dst(f )
(b)

=

X

X

αn (xf − xf [t − 1])2 1{n=Src(f )} +

X

X

(f )

(µl

l∈I(n)

(f )

(µl

(f )

− µl [t − 1])2 +

X

(f )

(µl

(f )

− µl [t − 1])2



l∈O(n)

(f )

− µl [t − 1])2

l∈I(n)


(f )
Uf (xf ) − WSrc(f ) [t]xf − αSrc(f ) (xf − xf [t − 1])2 +

f ∈F

−

(f )

− µl [t − 1])2

l∈O(n)

l∈I(n)

αn

(f )

(µl


X (f )
X (f ) 
Wn(f ) [t] xf 1{n=Src(f )} +
µl −
µl

f ∈F ,
n∈N \Dst(f )

X

X

αn

f ∈F ,n=Dst(f )

X

Uf (xf ) −

X

X

X

 (f )
(f )
Wn(f ) [t] − Wm
[t] µ(n,m)

(n,m)∈L f ∈F

X

X

(αn + αm )

(n,m)∈L

(f )
(µ(n,m)

−

(f )
µ(n,m) [t

2

− 1])

(30)

f ∈F

P
P
(f )
(f )
(f )
(f )
(f )
(f )
where (a) follows from the fact that gn (yn ) = xf 1{n=Src(f )} + l∈I(n) µl − l∈O(n) µl and kyn − yn [t − 1]k2 =
P
P
(f )
(f )
(f )
(f )
(xf − xf [t − 1])2 1{n=Src(f )} + l∈I(n) (µl − µl [t − 1])2 + l∈O(n) (µl − µl [t − 1])2 ; and (b) follows by collecting
(f )

(f )

(f )

(f )

each linear term µl and each quadratic term (µl − µl [t − 1])2 . Note that each link session rate µl appears twice with
P
P
P
(f )
(f )
(f ) 
opposite signs in the summation term f ∈F ,n∈N \{Dst(f )} Wn [t] xf 1{n=Src(f )} + l∈I(n) µl − l∈O(n) µl
unless link
(f )

l flows into Dst(f ) and recall that WDst(f ) = 0, ∀f ∈ F. The quadratic terms are collected in a similar way. Note that the
P
P
(f )
(f )
term f ∈F ,n=Dst(f ) αn l∈I(n) (µl − µl [t − 1])2 introduced to the objective function (28) is necessary to guarantee each
(f )

(f )

quadratic term (µ(m,n) − µ(m,n) [t − 1])2 with the same link index (n, m) but different flow indices f ∈ F have the same
coefficient αn + αm in the last line of (30).
(f )
Note that equation (30) is now separable for each scalar xf and vector [µ(n,m) ]f ∈F . Thus, problem (28)-(29) can be
decomposed into independent smaller optimization problems in the form of problem (12)-(13) with respect to each scalar xf ,
(f )
and in the form of problem (14)-(17) with respect to each vector [µ(n,m) ]f ∈F .

15

A PPENDIX F
P ROOF OF L EMMA 7
(f )

(f )

Note that Wn [t] appears as a known constant in (12). Since Uf (xf ) is concave and Wn [t]xf is linear, it follows that
(12) is strongly concave with respect to xf with modulus 2αn . Since xf [t] is chosen to solve (12)-(13), by Lemma 4, ∀f ∈ F,
we have
(f )

(f )

Uf (xf [t]) − WSrc(f ) [t]xf [t] − αn (xf [t] − xf [t − 1])2 ≥ Uf (x∗f ) − WSrc(f ) [t]x∗f − αn (x∗f − xf [t − 1])2 + αn (x∗f − xf [t])2 .
{z
} |
{z
}
|
(31)-I

(31)-II

(31)
Similarly, we know (14) is strongly concave with respect to vector [µf(n,m) ]f ∈F with modulus 2(αn + αm ). By Lemma 4,
∀(n, m) ∈ O(n), we have
X
 (f )
 X (f )
2
(f )
(f )
Wn(f ) [t] − Wm
[t] µ(n,m) [t] − αn + αm
µ(n,m) [t] − µ(n,m) [t − 1]
f ∈F

f ∈F

|
≥

{z

}

(32)-I

X

 (f ),∗
 X (f ),∗
2
 X (f ),∗
2
(f )
(f )
(f )
Wn(f ) [t] − Wm
[t] µ(n,m) − αn + αm
µ(n,m) − µ(n,m) [t − 1] + αn + αm
µ(n,m) − µ(n,m) [t] .

f ∈F

f ∈F

|

f ∈F

{z

}

(32)-II

(32)
(f )

Recall that each column vector yn defined in (21) is composed by control actions that appear in each constraint (2); column
P
(f )
vector y = [xf ; µl ]f ∈F ,l∈L is the collection of all control actions; and f (y) = f ∈F Uf (xf ). Summing term (31)-I over
all f ∈ F and term (32)-I over all (n, m) ∈ L and using an argument similar to the proof of Lemma 6 (Recall that y[t] is
jointly chosen to minimize (28) by Lemma 6.) yields
X
X
(31)-I +
(32)-I
f ∈F

(n,m)∈N

X

=f (y[t]) −


Wn(f ) [t]gn(f ) (yn(f ) [t]) + αn kyn(f ) [t] − yn(f ) [t − 1]k2 −

X
f ∈F ,
n=Dst(f )

f ∈F ,
n∈N \Dst(f )

(f )

X

αn

(f )

(µl [t] − µl [t − 1])2 . (33)

l∈I(n)


P
P
(f ),∗
(f )
(f )
(f ),∗
− µl [t])2 . Summing
− yn [t]k2 + αn 1{n=Dst(f )} l∈I(n) (µl
Recall that Φ[t] = f ∈F ,n∈N αn 1{n6=Dst(f )} kyn
term (31)-II over all f ∈ F and term (32)-II over all (n, m) ∈ L yields
X
X
X
(31)-II +
(32)-II = f (y∗ ) + Φ[t] − Φ[t − 1] −
Wn(f ) [t]gn(f ) (yn(f ),∗ ),
(34)
f ∈F

f ∈F ,
n∈N \Dst(f )

(n,m)∈N

Combining (31)-(34) and rearranging terms yields
f (y[t])
X

≥f (y∗ ) + Φ[t] − Φ[t − 1] −

f ∈F ,
n∈N \Dst(f )

+

X
f ∈F ,
n=Dst(f )

αn

X

(f )

X

Wn(f ) [t]gn(f ) (yn(f ),∗ ) +

Wn(f ) [t]gn(f ) (yn(f ) [t]) + αn kyn(f ) [t] − yn(f ) [t − 1]k2



f ∈F ,
n∈N \Dst(f )

(f )

(µl [t] − µl [t − 1])2

l∈I(n)

(a)

X

≥ f (y∗ ) + Φ[t] − Φ[t − 1] +

Wn(f ) [t]gn(f ) (yn(f ) [t]) + αn kyn(f ) [t] − yn(f ) [t − 1]k2



f ∈F ,
n∈N \Dst(f )
(b)

X

=f (y∗ ) + Φ[t] − Φ[t − 1] +

)
(f )
(f )
(f )
(f )
(f )
(f )
(f )
(f )
2
Q(f
n [t]gn (yn [t]) + gn (yn [t − 1])gn (yn [t]) + αn kyn [t] − yn [t − 1]k



f ∈F ,
n∈N \Dst(f )

(35)
(f )

(f ),∗

where (a) follows because gn (yn
0; (b) follows from the fact that

) = 0, ∀f ∈ F, ∀n ∈ N \Dst(f ), and

(f )
Wn [t]

=

(f )
Qn [t]

+

(f )
(f )
gn (yn [t

− 1]).

P

f ∈F ,n=Dst(f )

αn

(f )
(f )
2
l∈I(n) (µl [t]−µl [t−1])

P

≥

16

Recall that uT1 u2 = 12 u21 + 21 u22 − 12 (u1 − u2 )2 for any u1 , u2 ∈ R. Thus, for all f ∈ F, n ∈ N \ Dst(f ), we have
gn(f ) (yn(f ) [t − 1])gn(f ) (yn(f ) [t]) =

2 1
2 1
2
1 (f ) (f )
g (yn [t − 1]) + gn(f ) (yn(f ) [t]) − gn(f ) (yn(f ) [t − 1]) − gn(f ) (yn(f ) [t]) .
2 n
2
2

(36)

Substituting (36) into (35) yields
f (y[t]) ≥f (y∗ ) + Φ[t] − Φ[t − 1] +

X



)
(f )
(f )
Q(f
n [t]gn (yn [t]) +

f ∈F ,
n∈N \Dst(f )

2 1
2
1 (f ) (f )
gn (yn [t − 1]) + gn(f ) (yn(f ) [t])
2
2


2
1 (f ) (f )
gn (yn [t − 1]) − gn(f ) (yn(f ) [t]) + αn kyn(f ) [t] − yn(f ) [t − 1]k2
2

X
(a)
2 1
2
1 (f ) (f )
)
(f )
(f )
gn (yn [t − 1]) + gn(f ) (yn(f ) [t])
Q(f
≥ f (y∗ ) + Φ[t] − Φ[t − 1] +
n [t]gn (yn [t]) +
2
2
−

f ∈F ,
n∈N \Dst(f )


1 
+ αn − βn2 kyn(f ) [t] − yn(f ) [t − 1]k2
2

X
(b)
1 (f ) (f ) 2 
)
(f )
(f )
Q(f
≥f (y∗ ) + Φ[t] − Φ[t − 1] +
[t]g
(y
[t])
+
g (yn [t])
n
n
n
2 n

(37)

f ∈F ,
n∈N \Dst(f )
(f )

where (a) follows from the Fact 2, i.e., each gn (·) is Lipschitz with modulus βn and (b) follows because αn ≥ 12 (dn + 1),
2
√
(f )
(f )
βn ≤ dn + 1 and 12 gn (yn [t − 1]) ≥ 0.
Subtracting (26) from (37) and cancelling the common terms on both sides yields
−∆[t] + f (y[t]) ≥ f (y∗ ) + Φ[t] − Φ[t − 1].
R EFERENCES
[1] L. Tassiulas and A. Ephremides, “Stability properties of constrained queueing systems and scheduling policies for maximum throughput in multihop
radio networks,” IEEE Transactions on Automatic Control, vol. 37, no. 12, pp. 1936–1948, 1992.
[2] M. J. Neely, “Dynamic power allocation and routing for satellite and wireless networks with time varying channels,” Ph.D. dissertation, Massachusetts
Institute of Technology, 2003.
[3] L. Georgiadis, M. J. Neely, and L. Tassiulas, “Resource allocation and cross-layer control in wireless networks,” Foundations and Trends in Networking,
2006.
[4] M. J. Neely, Stochastic network optimization with application to communication and queueing systems. Morgan & Claypool Publishers, 2010.
[5] A. Eryilmaz and R. Srikant, “Joint congestion control, routing, and mac for stability and fairness in wireless networks,” IEEE Journal on Selected Areas
in Communications, vol. 24, no. 8, pp. 1514–1524, 2006.
[6] A. L. Stolyar, “Maximizing queueing network utility subject to stability: Greedy primal-dual algorithm,” Queueing Systems, vol. 50, no. 4, pp. 401–457,
2005.
[7] X. Lin and N. B. Shroff, “Joint rate control and scheduling in multihop wireless networks,” in Proceedings of IEEE Conference on Decision and Control
(CDC), 2004.
[8] J.-W. Lee, R. R. Mazumdar, and N. B. Shroff, “Opportunistic power scheduling for dynamic multi-server wireless systems,” IEEE Transactions on
Wireless Communications, vol. 5, no. 6, pp. 1506–1515, 2006.
[9] M. J. Neely, E. Modiano, and C. E. Rohrs, “Dynamic power allocation and routing for time-varying wireless networks,” IEEE Journal on Selected Areas
in Communications, vol. 23, no. 1, pp. 89–103, 2005.
[10] J. Liu, N. B. Shroff, C. H. Xia, and H. D. Sherali, “Joint congestion control and routing optimization: An efficient second-order distributed approach,”
IEEE/ACM Transactions on Networking, vol. 24, no. 3, pp. 1404–1420, 2015.
[11] M. J. Neely, “Super-fast delay tradeoffs for utility optimal fair scheduling in wireless networks,” IEEE Journal on Selected Areas in Communications,
vol. 24, no. 8, pp. 1489–1501, 2006.
[12] L. Huang, S. Moeller, M. J. Neely, and B. Krishnamachari, “LIFO-backpressure achieves near-optimal utility-delay tradeoff,” IEEE/ACM Transactions
on Networking, vol. 21, no. 3, pp. 831–844, 2013.
[13] E. Wei, A. Ozdaglar, and A. Jadbabaie, “A distributed Newton method for network utility maximization–I: algorithm,” IEEE Transactions on Automatic
Control, vol. 58, no. 9, pp. 2162–2175, 2013.
[14] M. Zargham, A. Ribeiro, and A. Jadbabaie, “Accelerated backpressure algorithm,” in Proceedings of IEEE Global Communications Conference
(GLOBECOM), 2013.
[15] N. Parikh and S. Boyd, “Proximal algorithms,” Foundations and Trends in Optimization, vol. 1, no. 3, pp. 123–231, 2013.
[16] X. Lin and N. B. Shroff, “Utility maximization for communication networks with multipath routing,” IEEE Transactions on Automatic Control, vol. 51,
no. 5, pp. 766–781, 2006.
[17] F. P. Kelly, A. K. Maulloo, and D. K. Tan, “Rate control for communication networks: Shadow prices, proportional fairness and stability,” Journal of
the Operational Research Society, pp. 237–252, 1998.
[18] S. H. Low and D. E. Lapsley, “Optimization flow control—I: basic algorithm and convergence,” IEEE/ACM Transactions on Networking, vol. 7, no. 6,
pp. 861–874, 1999.
[19] D. P. Bertsekas, A. Nedić, and A. E. Ozdaglar, Convex Analysis and Optimization. Athena Scientific, 2003.
[20] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra, “Efficient projections onto the l1 -ball for learning in high dimensions,” in Proceedings of
International Conference on Machine learning (ICML), 2008.
[21] M. J. Neely, “A simple convergence time analysis of drift-plus-penalty for stochastic optimization and convex programs,” arXiv preprint arXiv:1412.079,
2014.
[22] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, 2004.

