A Distributed Method for Optimal Capacity Reservation
Nicholas Moehle∗

Xinyue Shen†

Zhi-Quan Luo‡

Stephen Boyd§

arXiv:1705.00677v1 [cs.DC] 1 May 2017

May 3, 2017

Abstract
We consider the problem of reserving link capacity in a network in such a way
that any of a given set of flow scenarios can be supported. In the optimal capacity
reservation problem, we choose the reserved link capacities to minimize the reservation
cost. This problem reduces to a large linear program, with the number of variables and
constraints on the order of the number of links times the number of scenarios. Small
and medium size problems are within the capabilities of generic linear program solvers.
We develop a more scalable, distributed algorithm for the problem that alternates
between solving (in parallel) one flow problem per scenario, and coordination steps,
which connect the individual flows and the reservation capacities.

1

Introduction

We address the capacity reservation problem, the problem of reserving link capacity in a
network in order to support multiple possible flow patterns. We are given a description of
the network, and a collection of traffic demand scenarios, which specify the amount of flow
to be routed from some sources to some sinks. We must reserve enough link capacity so that
for any of these scenarios, there exist flows that route the traffic demand while using no more
than the reserved capacity on each link. Subject to this requirement, we seek to minimize a
linear reservation cost. This problem is also referred to as the network design problem.
The capacity reservation problem can be expressed as a linear program (LP), a fact first
recognized by Gomory and Hu in 1962 [GH62]. For moderate problem sizes it can be solved
using generic LP solvers. However, problems with a large number of scenarios (and especially
those that do not fit in a single computer’s memory) are beyond the reach of such solvers.
We present an iterative method for solving the capacity reservation problem based on the
alternating direction method of multipliers (ADMM). Each iteration of the algorithm involves
solving, in parallel, a single minimum-cost flow problem for each scenario. This means that
∗

Mechanical Engineering Department, Stanford University. moehle@stanford.edu
Electronic Engineering Department, Tsinghua University.
‡
Electrical and Computer Engineering, University of Minnesota.
§
Electrical Engineering Department, Stanford University.
†

1

our method can easily exploit parallel computing architectures, and can scale to enormous
problem instances that do not fit in the memory of a single computer. Unlike general ADMM,
our method maintains a feasible point at each iteration (and therefore an upper bound on
the problem value). We can also compute a lower bound at modest cost, which can be used
to bound the suboptimality of the current iterate, allowing for non-heuristic termination of
the algorhtm.

1.1

Previous work

Capacity reservation problem with finite source set. The optimal reservation problem has a long history, and was first proposed by Gomory and Hu in 1962 [GH62], who
consider the case of a finite source set, and a single commodity. In their extension to the
multicommodity case, they propose a decomposition method based on duality [GH64]. They
also note that the structure of the capacity reservation problem makes decomposition techniques attractive. Labbé et al. [LSSW99] propose two decomposition techinques, one of
which is based on Lagrangian duality and foreshadows the method we propose here. Other
decomposition approaches are based on cutting plane methods, which can be interpreted as
using only a small subset of the source set, and adding additional source vectors as necessary. Examples of this approach include the work of Minoux et al. [Min81] and Petrou et al.
[PLO07].
Capacity reservation problem with infinite source set. Several other source set
descriptions are given in the literature, but are typically intractable. For example, when the
source set is given by polyhedron (described by a set of inequalities), the capacity reservation
problem is known to be NP-hard [CSOS07]. However, for this and other cases, such as
ellipsoidal and conic source sets, several effective heuristics exist. One tractable approach is
to restrict the flows to be an affine function of the demands. This approach was first proposed
by Ben-Ameur and Kerivin [BAK03, BAK05], and is called static routing, oblivious routing,
dynamic routing, or affine routing, depending on the form of the affine function used. For a
summary, see [PR11]. Such a heuristic can be viewed as an application of affinely adjustable
robust optimization, has been studied extensively in robust optimization; see [BTGN09].
Convex optimization. In the case of a finite source set, the capacity reservation problem can be reduced to a convex optimization problem (indeed, a linear program). These
problems are tractable and mature software exists that can be used to solve them [BV04].
Our approach to the capacity reservation problem is based on a particular algortithm for
solving convex optimization problems, called the alternating direction method of multipliers (ADMM) [PB14, BPC+ 11], which is well-suited for solving large convex problems on
distributed computers.

2

1.2

Outline

We start by describing the capacity reservation problem in §2, and a simple heuristic for
solving it. In §3, we provide a duality-based method for obtaining lower bounds on the problem value. In §4, we illustrate the heuristic, as well as the lower bounds, on two illustrative
examples. In §5, we give an ADMM-based distributed algorithm for solving the capacity
reservation problem. We conclude with a numerical example in §6.

2

Capacity reservation problem

We consider a single-commodity flow on a network. The network is modeled as a directed
connected graph with n nodes and m edges, described by its incidence matrix A ∈ Rn×m ,

 1 edge j points to node i
−1 edge j points from node i
Aij =

0 otherwise.
We let f ∈ Rm denote the vector of edge flows, which we assume are nonnegative and limited
by a given edge capacity, i.e., 0 ≤ f ≤ c, where the inequalities are interpreted elementwise,
and c ∈ Rm is the vector of edge capacities. We let s ∈ Rn denote the vector of sources,
i.e., flows that enter (or leave, when negative) the node. Flow conservation at each node is
expressed as Af + s = 0. (This implies that 1T s = 0.) We say that f is a feasible flow for
the source s if there exists f that satisfies Af + s = 0 and 0 ≤ f ≤ c. For a given source
vector s, finding a feasible flow (or determining that none exists) reduces to solving a linear
program (LP).
We consider a set of source vectors S ⊂ Rn , which we call the source set. A flow policy
is a function F : S → Rm . We say that a flow policy is feasible if for each s ∈ S, F(s) is a
feasible flow for s. Roughly speaking, a flow policy gives us a feasible flow for each possible
source. Our focus is on the selection of a flow policy, given the network (i.e., A and c) and
a description of the source set S.
A reservation (or reserved capacity) is a vector r ∈ Rm that satisfies 0 ≤ r ≤ c. We
say that the reservation r supports a flow policy F if for every source vector s ∈ S, we have
F(s) ≤ r. Roughly speaking, we reserve a capacity rj on edge j; then, for any possible source
s, we can find a feasible flow f that satisfies fj ≤ rj , i.e., it does not use more than the
reserved capacity on each edge. The cost of the reservation r is pT r, where p ∈ Rm , p ≥ 0 is
the vector of edge reservation prices. The optimal capacity reservation (CR) problem is to
find a flow policy F and a reservation r that supports it, while minimizing the reservation
cost.
The CR problem can be written as
minimize pT r
subject to AF(s) + s = 0,
r ≤ c.
3

0 ≤ F(s) ≤ r,

∀s ∈ S

(1)

The variables are the reservation r ∈ Rm and the flow policy F : S → Rm . The data are A,
S, c, and p. This problem is convex, but when S is infinite, it is infinite dimensional (since
the variables include a function on an infinite set), and contains a semi-infinite constraint,
i.e., linear inequalities indexed by an infinite set. We let J ? denote the optimal value of the
CR problem.
We will focus on the special case when S is finite, S = {s(1) , . . . , s(K) }. We can think of
s(k) as the source vector in the kth scenario. The CR problem (1) is tractable in this case,
indeed, an LP:
minimize pT r
subject to Af (k) = s(k) ,
r ≤ c.

0 ≤ f (k) ≤ r,

k = 1, . . . , K,

(2)

The variables are the reservations r ∈ Rm and the scenario flow vectors f (k) ∈ Rm for
each scenario k. This LP has m(K + 1) scalar variables, nK linear equality constraints,
and m(2K + 1) linear inequality constraints. In the remainder of the paper, we assume the
problem data are such that this problem is feasible, which occurs if and only if for each
scenario there is a feasible flow. (This can be determined by solving K independent LPs.)
We note that a solution to (2) is also a solution for the problem (1) with the (infinite)
source set S = Co{s(1) , . . . , s(K) }, where Co denotes the convex hull. In other words, a
reservation that is optimal for (2) is also optimal for the CR problem in which the finite
source set S is replaced with its convex hull. To see this, we note that the optimal value of
(2) is a lower bound on the optimal value of (1) with S the convex hull, since the feasible
set of the former includes the feasible set of the latter. But we can extend a solution of the
LP (2) to a feasible point for the CR problem (1) with the same objective, which therefore
is optimal. To create such as extension, we define F(s) for any s ∈ S. We define θ(s) as
the unique least Euclidean norm vector
with θ(s) ≥ 0, 1T θ(s) = 1 (with 1 the vector with
P
all entries one) that satisfies s = k θk s(k) . (These are barycentric coordinates.) We then
P
(k)
.
define F(s) = K
k=1 θk (s)f
We note for future use another form of the CR problem (2). We eliminate r using
r = maxk f (k) , where the max is understood to apply elementwise. Then we obtain the
problem
minimize pT maxk f (k)
(3)
subject to Af (k) = s(k) , 0 ≤ f (k) ≤ c, k = 1, . . . , K,
where the variables are f (k) , k = 1, . . . , K. This is a convex optimization problem.
We will also express the CR problem using matrices, as
minimize pT max(F )
subject to AF + S = 0,
0≤F ≤C

(4)

with variable F ∈ Rm×K . The columns of F are the flow vectors f (1) , . . . , f (K) . The columns
of S are the source vectors s(1) , . . . , s(K) , and the columns of C are all the capacity vector
4

c, i.e., C = c1T , where 1 is the vector with all entries one. The inequalities above are
interpreted elementwise, and the max of a matrix is taken over the rows, i.e., it is the vector
containing the largest element in each row.
Heuristic flow policy. We mention a simple heuristic for the reservation problem that
foreshadows our method. We choose each scenario flow vector f (k) as a solution to the
capacitated minimum-cost flow problem
minimize pT f
subject to Af = s(k) ,

(5)

0 ≤ f ≤ c,

with variable f ∈ Rm . These flow vectors are evidently feasible for (3), and therefore the
associated objective value J heur is an upper bound on the optimal value of the problem, i.e.,
J ? ≤ J heur = pT max f (k) .
k

Note that finding the flow vectors using this heuristic involves solving K smaller LPs independently, instead of the one large LP (2). This heuristic flow policy greedily minimizes the
cost for each source separately; but it does not coordinate the flows for the different sources
to reduce the reservation cost.
We will now show that we also have
J heur /K ≤ J ? ,
i.e., J heur /K is a lower bound on J ? . This shows that the heuristic policy is no more than
K-suboptimal. We first observe that for any nonnegative vectors f (k) , we have
!
K
K
X
X
f (k) ≤ pT max f (k) ,
(6)
pT f (k) = pT (1/K)
(1/K)
k=1

k=1

k

which follows from the fact that the maximum of a set of K numbers is at least as big as
their mean. Now minimize the lefthand and righthand sides over the feasible set of (2). The
the righthand side becomes J ? , and the lefthand side becomes J heur .
In §4.2 we give an example showing that this suboptimality bound is tight, i.e., the
heuristic flow policy can indeed be up to K suboptimal.

3

Lower bounds and optimality conditions

The heuristic flow policy described above provides an upper and lower bound on J ? . These
bounds are computationally appealing because they involve solving K small LPs independently. Here we extend the basic heuristic method to one in which the different flows are
computed using different prices. This gives a parametrized family of lower and upper bounds
on J ? . Using ideas based on linear programming duality, we can show that this family is
tight, i.e., there is a choice of prices for each scenario that yields lower and upper bound J ? .
This duality idea will form the basis for our distributed solution to (2), described in §5.
5

Scenario prices. P
We consider nonnegative scenario pricing vectors π (k) ∈ Rm , for k =
(k)
1, . . . , K, satisfying K
= p.
k=1 π
Then we have
K
X
pT max f (k) ≥
π (k)T f (k) .
(7)
k

k=1

It follows that the optimal value of the LP
PK (k)T (k)
minimize
f
k=1 π
subject to Af (k) = s(k) , 0 ≤ f (k) ≤ c,

k = 1, . . . , K

(8)

is a lower bound on J ? . This problem is separable; it can be solved by solving K small
LPs (capacited flow problems) independently. In other words, to obtain the bound, we
first decompose the reservation prices on each edge into scenario prices. Then, for each
scenario k, the flow f (k) is chosen (independently of the other scenarios) to be a minimumcost flow according to price vector π (k) . The scenario price vectors can be interpreted as a
decomposition of the reservation price vector, i.e., the price of reservation along each edge
can be decomposed into a collection of flow prices for that edge, with one for each scenario.
Note that the inequality (6) is a special case of (7), obtained with scenario prices π (k) =
(1/K)p. In fact, with these scenario prices, the heuristic subproblem (5) is also just a special
case of (8), i.e., the family of lower bounds obtained by solving (8) with different scenario
pricing vectors includes the single lower bound obtained using the heuristic method.
Given a set of scenario prices, the flows found as above are feasible for the CR problem.
Therefore pT maxk f (k) is an upper bound on J ? .
Tightness of the bound and optimality conditions. There exist some scenario prices
for which the optimal value J ? of the CR problem (2) and its relaxation (8) are the same.
This is a consequence of linear programming duality; a (partial) dual of (2) is the problem of
determining the scenario pricing vectors that maximize the optimal value of (8). In addition,
given such scenario prices, any optimal flow policy for (2) is also optimal for (8) with these
scenario prices.
We can formalize this notion as follows. A flow policy f (k) is optimal if and only if there
exists a collection of scenario prices π (k) such that the following conditions hold.
1. The vectors π (k) must be valid scenario prices, i.e.,
K
X

(k)

πj

= πj ,

π (k) ≥ 0,

k = 1, . . . , K.

k=1

2. The inequality (7) must be tight:
pT max f (k) =
k

K
X
k=1

6

π (k)T f (k) .

(9)

3. The flow vectors f (k) , for k = 1, . . . , K, must be optimal for (5), with scenario prices
π (k) .
Condition 2 has the interpretation that the optimal reservation cost can be written as a
sum of K cost terms corresponding to the scenarios, i.e., the reservation cost can be allocated
to the K scenarios. This can be made the basis for a payment scheme, in which scenario k is
charged with a fraction π (k)T f (k) of the total reservation cost. Note that because the scenario
prices that show optimality of a collection of scenario flow vectors are not necessarily unique,
these cost allocations may not be unique.
Condition 2 also implies the complimentarity condition
(k)T

πj

(f (k) − r)j = 0

for all k and j. In other words, a positive scenario price on an edge implies that the
corresponding scenario edge flow is equal to the optimal reservation, i.e., this scenario edge
flow must contribute to the reservation cost. Similarly, if a scenario edge flow is not equal
to that edge’s reservation, then the scenario price must be zero for that edge.
Condition 3 means that, given the scenario prices, the flow vector for each scenario is
a minimum-cost flow using the scenario price vector associated with that scenario. We can
interpret the price vectors as weighting the contribution of a scenario edge flow to the total
(k)
reservation cost. In the case that πj = 0, then the flow over edge j under scenario k does
(k)
not contribute to the total reservation cost. If πj > 0, then the flow over edge j under
(k)
scenario k contributes to the total reservation cost, and from before, we have fj = rj .
Additionally, if π (k) = 0 for some k, then scenario k is irrelevant i.e., removing it has no
effect on the reservation cost. If instead π (k) = p, we conclude that scenario k is the only
relevant scenario.

4
4.1

Examples
Illustrative example

Here we consider a simple example with n = 5 nodes, m = 10 edges, and K = 8 scenarios.
The graph was generated randomly, from a uniform distribution on all connected graphs
with 5 nodes and 10 edges, and is shown in figure 1. We use price vector p = 1 and capacity
vector c = 1. The 8 scenario source vectors were chosen according to s(k) = −Az (k) , where
each element of z (k) was randomly chosen from {0, 1/3, 2/3, 1}. This guarantees that there
is a feasible flow for each scenario.
The optimal reservation cost is 6.0, and the objective of the heuristic policy is 7.6. (The
lower bound from the heuristic policy is 2.3.) The optimal and heuristic flow policies are
shown in figure (2). The upper plot shows the optimal policy, and the lower plot shows the
heuristic policy. For each plot, the bars show the flow policy; the 10 groups are the edges,
and the 8 bars are the edge flows under each scenario. The line above each group of bars is
the reservation for that edge.
7

1

8
2

10
2

5
1

6

5

7

9

3
3

4

4

Figure 1: The graph of the illustrative example.
Optimal scenario prices are given in table 1. First, we confirm that the row sums are
indeed p. The fact that only some scenario prices are positive for each edge indicates that no
scenario is relevant for every edge, i.e., that there is no single scenario that is being planned
for. However, for every scenario, the scenario prices are positive for at least one edge, which
means that every scenario is potentially relevant to the optimal reservation.

4.2

K-suboptimal example

In §2, we showed that the heuristic method is not more than K suboptimal, and the previous
example showed that the heuristic approach is not, in fact, optimal. Here we give an example
that shows that the K-suboptimality result is in fact tight, i.e., for every positive integer K,
there is an instance of (2) for which the heuristic is (nearly) K suboptimal.
We consider a graph with n = 2a + 1 nodes, for some positive integer a, with the nodes
arranged in three layers, with a, a, and 1 nodes, respectively. (We show the graph for a = 3
in figure 3.) There is an edge connecting every node in the first layer with every node in
the second layer, i.e., for every i = 1, . . . , a and every j = a + 1, . . . , 2a, the edge (i, j) is in
the graph. The price in this edge is  if i = j, and 2 if i 6= j. Similarly, every node in the
second layer connects to the last node, i.e., for every i = a + 1, . . . , 2a, the edge (i, 2a + 1)
is in the graph. The price for this edge is 1.
We consider K = a scenarios. In each scenario, one unit of flow in injected into a single
node in the first layer, and must be routed to the last node. In other words, the source
vector s(k) is all zeros, except for the kth element, which is 1, and the last element, which is

8

flows (optimal)
flows (naive)

1.2
1.0
0.8
0.6
0.4
0.2
0.0
1.2
1.0
0.8
0.6
0.4
0.2
0.0

2

4

2

4

edge

6

8

10

6

8

10

Figure 2: The optimal and naive flow policies.

Edge

1
1
2
3
4
5
6
7
8
9
10

2

3
1.0

0.33

Scenario
4
5
0.33

0.38
1.0
1.0
1.0
0.38

6

7

8

0.33
0.28

0.33

0.62
0.33

0.33

0.33

0.33

0.33
1.0
0.33

Table 1: The scenario prices π (1) , . . . , π (8) . Column k in the table is the scenario price vector
π (k) . Blank spaces correspond to 0 entries.
9



1
2
2
2

2

4
1



5

2
2
3

1

7

1
2


6

Figure 3: Graph example for a = 3. The labels give the reservation price for that edge.
−1.
A feasible (but not necessarily optimal) point for the CR problem (2) is to route all flow
equally through node a + 1, the first node in the second layer. We must then reserve one
unit of capacity for edge (a + 1, 2a + 1), and one unit of capacity from each node in the first
layer to node a + 1, i.e., for edges (i, a + 1), for i = 1, . . . , a. The total reservation cost is
(2a − 1) + 1.
We now consider the (unique) heuristic solution. The solution to (5) under scenario k is
to route all flow from node k to node a + k, and finally to node 2a + 1, Because no edges
support flow under any two different scenarios, the heuristic therefore achieves a reservation
cost of K( + 1).
As  → 0, the cost difference between the feasible reservation for (2) and the heuristic
reservation approaches K. Because the feasible reservation provides an upper bound on the
optimal reservation cost, the difference between the optimal reservation and the heuristic
reservation also approaches K.

5

Distributed solution method

In this section we provide a parallelizable algorithm, based on the alternating direction
method of multipliers (ADMM), to solve (2). (For details on ADMM, see [BPC+ 11].) We
first express the problem (4) in the consensus form
minimize pT max(F̃ ) + g(F )
subject to F = F̃

(10)

with variables F ∈ Rm×K and F̃ ∈ Rm×K . The function g is the indicator function for
feasible flows,

0 AF + S = 0, 0 ≤ F ≤ C
g(F ) =
∞ otherwise.
In the consensus problem (10) we have replicated the flow matrix variable, and added a
consensus constraint, i.e., the requirement that the two variables representing flow policy
10

must agree.
The augmented Lagrangian for problem (10) is
L(F, F̃ , Π) = pT max(F̃ ) + g(F ) + Tr ΠT (F − F̃ ) + (ρ/2)kF − F̃ k2F ,
where k · kF is the Frobenius norm, and ρ > 0 is a parameter. Here Π ∈ Rm×K is the dual
variable.
ADMM. The ADMM algorithm (with over-relaxation) for (10) is given below. First we
initialize the iterates F̃ (0) and Π(0). Then we carry out the following steps:

F (l + 1) = argmin L F, F̃ (l), Π(l)
F

F̃ (l + 1) = argmin L αF (l + 1) + (1 − α)F̃ (l), F̃ , Π(l)
F̃


Π(l + 1) = Π(l) + ρ αF (l + 1) + (1 − α)F̃ (l) − F̃ (l + 1) ,
where α is an algorithm parameter in (0, 2); the argument l is the iteration number. For
reasons that will become clear, we call these three steps the flow policy update, the reservation
update, and the price update, respectively.
Convergence. Given our assumption that (2) is feasible, we have F (l) → F ? , F̃ (l) → F ? ,
and Π(l) → Π? , where (F ? , Π? ) is a primal-dual solution to (4). This follows from standard
results on ADMM [BPC+ 11].
Flow policy update. Minimizing the Lagrangian over F requires minimizing a quadratic
function over the feasible set of (4). This can be interpreted as the solution of a minimumcost flow problem for each of the K scenarios. More specifically, column k of the matrix
F (l + 1) is the unique solution to the quadratic minimum-cost flow problem:
minimize π (k)T f + (ρ/2)kf − f˜(k) k2
subject to Af + s(k) = 0
0 ≤ f ≤ c.

(11)

Here f˜(k) and π (k) are the kth columns of the matrices F̃ (l), and Π(l), respectively. These
K flow problems can be solved in parallel to update the K columns of the iterate F (l + 1).
Note that the objective can be interpreted as a sum of the edge costs, using the estimate of
the scenario prices, plus a quadratic regularization term that penalizes deviation of the flow
vector from the previous flow values f˜(k) .
Reservation update. Minimizing the Lagrangian over F̃ requires minimizing the sum of
a quadratic function of F̃ plus a positive weighted sum of its row-wise maxima. This step

11

decouples into m parallel steps. More specifically, The jth row of F̃ (l + 1) is the unique
solution to
 2

(12)
,
minimize pj max f˜ − πjT f˜ + (ρ/2) f˜ − αfj + (1 − α)f˜j
where f˜ is the variable, and fj , f˜j , and πj are the jth rows of F (l + 1), F̃ (l), and Π(l),
respectively. (The scalar pj is the jth element of the vector p.) We interpret the max of a
vector to be the largest element of that vector. This step can be interpreted as an implicit
update of the reservation vector. The first two terms are the difference between the two sides
of (7). The problem can therefore be interpreted as finding flows for each edge that minimize
the looseness of this bound while not deviating much from a combination of the previous
flow values fj and f˜j . (Recall that at a primal-dual solution to (2), this inequality becomes
tight.) There is a simple solution for each of these subproblems, given in appendix A; the
computational complexity for each subproblem scales like K log K.
Price update. In the third step, the dual variable Π(l + 1) is updated.
Iterate properties. Note that the iterates F (l), for l = 1, 2 . . . , are feasible for (4), i.e.,
the columns of F (l) are always feasible flow vectors. This is a result of the simple fact that
each column is the solution of the quadratic minimum-cost flow problem (11). This means
that U (l) = pT max F (l) is an upper bound on the optimal problem value. Furthermore,
because F (l) converges to a solution of (4), this upper bound U (l) converges to J ? (though
not necessarily monotonically).
It can be shown that Π(l) ≥ 0 and Π(l)1 = p for l = 1, 2, . . . , i.e., the columns of Π(l)
are feasible scenario price vectors satisfying optimality condition 1 of §3. This is proven
in Appendix B. We can use this fact to obtain a lower bound on the problem value, by
computing the optimal value of (8), where the scenario pricing vectors π (k) are given by the
columns of Π(l), which is a lower bound on the optimal value of (4). We call this bound L(l).
Because Π(l) converges to an optimal dual variable of (10), L(l) converges to J ? (though it
need not converge monotonically).
Additionally, the optimality condition 2 is satisfied by iterates F̃ (l + 1) and Π(l + 1), for
l = 1, 2, . . . , if we take f (k) and π (k) to be the columns F̃ (l + 1) and Π(l + 1), respectively.
This is shown in Appendix B.
Stopping criterion. Several reasonable stopping criteria are known for ADMM. One standard stopping criterion is that that the norms of the primal residual kF (l) − F̃ (l)k and the
dual residual kΠ(l + 1) − Π̃(l)k are small. (For a discussion, see [BPC+ 11, §3.3.1].)
In our case, we can use the upper and lower bounds U (l) and L(l) to bound the suboptimality of the current feasible point F (l). More specifically, we stop when
U (l) − L(l) ≤ rel L(l)

(13)

which guarantees a relative error not exceeding rel . Note that computing L(l) requires
solving K small LPs, which has roughly the same computational cost as one ADMM iteration;
12

it may therefore be beneficial to check this stopping criterion only occasionally, instead of
after every iteration.
Parameter choice. Although our algorithm converges for any positive ρ and any α ∈
(0, 2), the values of these parameters strongly impact the number of iterations required. In
many numerical experiments, we have found that choosing α = 1.8 works well. The choice
of the parameter ρ is more critical. Our numerical experiments suggest choosing ρ as
ρ=µ

1T p
,
max(1T F heur )

(14)

where F heur is a solution to the heuristic problem described in §2, and µ is a positive constant.
With this choice of ρ, the primal iterates F (l) and F̃ (l) are invariant under positive scaling
of p, and scale with S and c. (Similarly, the dual iterates Π(l) scale with p, but are invariant
to scaling of S and c.) Thus the choice (14) renders the algorithm invariant to scaling of p,
and also to scaling of S and c. Numerical experiments suggest that values of µ in the range
between 0.01 and 0.1 work well in practice.
Initialization. Convergence of F (l), F̃ (l), and Π(l) to optimal values is guaranteed for any
initial choice of F̃ (0) and Π(0). However we propose initializing F̃ (0) as a solution to the
heuristic method described in §2, and Π(0) as (1/K)p1T . (In this case, we have F (1) = F̃ (0),
which means the initial feasible point is no more than K suboptimal.)

6

Numerical example

Problem instance. We consider a numerical example with n = 2000 nodes, m = 5000
edges, and K = 1000 scenarios, for which the total number of (flow) variables is 5 million.
The graph is randomly chosen from a uniform distribution on all graphs with 2000 nodes
and 5000 edges. The elements of the price vector p are uniformly randomly distributed over
the interval [0, 1], and c = 1. The elements of the source vectors are generated according to
s(k) = −Az (k) , where the elements of z (k) are uniformly randomly distributed over the interval
[0, 1]. The optimal value is 18053, and the objective obtained by the heuristic flow policy is
24042. The initial lower bound is 5635, so the initial relative gap is (U (0)−L(0))/L(0) = 3.26.
Convergence. We use α = 1.8, with ρ chosen according to (14) with µ = 0.05. We
terminate when (13) is satisfied with rel = 0.01, i.e., when we can guarantee that the flow
policy is no more than 1% suboptimal. We check the stopping criterion at every iteration;
however, we only update the lower bound (i.e., evaluate L(l)) when l is a multiple of 10.
The algorithm terminates in 95 iterations. The convergence is shown in figure 4. The
upper plot shows the convergence of the upper and lower bounds to the optimal problem
value, and the lower plot shows the relative gap (U (l) − L(l))/L(l) and the realtive suboptimality (U (l) − J ? )/J ? (computed after the algorithm has converged to very high accuracy).
13

25000
20000
15000

Upper bound
Lower bound

10000
5000

0

20

40

60

80

100

101
10

Rel. gap
Rel. subopt.

0

10−1
10−2
10−3

0

20

40

60

80

100

Figure 4: Top: Upper and lower bounds U (l) and L(l) on the optimal value. Bottom:
Estimated relative gap (U (l) − L(l))/L(l) and relative suboptimality (U (l) − J ? )/J ? .
(Note that both the upper and the lower bounds were updated by the best values to date.
i.e., they were non-increasing and non-decreasing, respectively. The lower bound was only
updated once every 10 iterations.) Our algorithm guarantees that we terminate with a policy
no more than 1% suboptimal; in fact, the final policy is around 0.4% suboptimal.
Timing results. We implemented our method in Python, using the multiprocessing package. The flow update subproblems were solved using the commercial interior-point solver
Mosek [ApS15], interfaced using CVXPY, a Python-based modeling language for convex
optimization [DB16]. The reservation update subproblems were solved using the simple
method given in Appendix A. (The simulation was also run using the open source solver
ECOS [DCB13], which however is single-threaded, and unable to solve the full problem directly. Mosek is able to solve the full problem directly, which allows us to compare the solve
time with our distributed method.)
We used a Linux computer with 4 Intel Xeon E5-4620 CPUs, with 8 cores per CPU and
2 virtual cores per physical core. The flow update and reservation update subproblems were
solved in parallel with 64 hyperthreads.

14

method
flow update
distributed method
full CR problem

solve time (s)
1.2
6268
45953

Table 2: The timing results for the numerical example.
The computing times are listed in table 2. The first line gives the average time to solve a
single flow update subproblem (11). The second line gives the total time taken to solve the
problem using our distributed method (which took 95 iterations). The third line gives the
time taken to solve the full CR problem directly. For this simulation, we terminated Mosek
when the relative surrogate duality gap in the interior-point method reached 0.01. (Although
this is a similar stopping criterion to the one used for our method, it is not equivalent to,
nor a sufficient condition for, having a relative gap less than 0.02.)
For this example, our distributed method is about 7× faster than Mosek. This timing
advantage will of course scale with the number of cores or threads, since each iteration of
our method is trivially parallelizable, whereas the sparse matrix factorization that dominates
each iteration of an interior-point method is not. By scaling the problem size just a bit
larger, we reach a point where direct solution of the CR is problem is no longer possible; the
distributed method, however, will easily handle (much) larger problems.

7

Extensions

Capacitated nodes. In our formulation, only the links have capacity limits. However, in
many practical cases, nodes may also have capacity limits, with a corresponding reservation
cost.
Given a problem with a capacitated node, we show how to formulate an equivalent
problem, with only capacitated edges. Consider a capacitated node i with capacity a. We
can define new graph, with node i replaced by two nodes, i1 and i2 ; all edges previously
entering i now enter i1 , and all edges previously leaving i now leave i2 . We then connect i1
to i2 with a single new edge, with capacity a and edge reservation price 0.
General reservation costs. The cost of reserving capacity on edge j can be any convex
function φj of the reservation rj .

15

Multicommodity network. We can extend (2) to handle T different types of commodities. The multicommodity capacity reservation problem is
minimize pT r
subject to Af (k,t) = s(k,t) , 0 ≤ f (k,t)
PT
(k,t)
≤ r,
t=1 f
r ≤ c.

k = 1, . . . , K,
k = 1, . . . , K,

t = 1, . . . , T

(15)

The variables are the reservations r ∈ Rm and f (k,t) ∈ Rm for k = 1, . . . , K and t = 1, . . . , T .
The vector is s(k,t) is the source vector associated with type t under scenario k. We note
that the scenario prices of §3 extend immediately to the multicommodity case.

16

References
[ApS15]

MOSEK ApS. The MOSEK optimization toolbox for MATLAB manual. Version
7.1 (Revision 28), 2015.

[BAK03]

W. Ben-Ameur and H. Kerivin. New economical virtual private networks. Communications of the ACM, 46(6):69–73, 2003.

[BAK05]

W. Ben-Ameur and H. Kerivin. Routing of uncertain traffic demands. Optimization and Engineering, 6(3):283–313, 2005.

[BPC+ 11] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers.
Foundations and Trends in Machine Learning, 3(1):1–122, 2011.
[BTGN09] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski. Robust optimization. Princeton
University Press, 2009.
[BV04]

S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University
Press, 2004.

[CSOS07] C. Chekuri, F. B. Shepherd, G. Oriolo, and M. G. Scutellá. Hardness of robust
network design. Networks, 50(1):50–54, 2007.
[DB16]

S. Diamond and S. Boyd. CVXPY: A Python-embedded modeling language for
convex optimization. Journal of Machine Learning Research, 17(83):1–5, 2016.

[DCB13]

A. Domahidi, E. Chu, and S. Boyd. ECOS: An SOCP solver for embedded
systems. In Proceedings of the 12th European Control Conference, pages 3071–
3076. IEEE, 2013.

[GH62]

R. E. Gomory and T. C. Hu. An application of generalized linear programming
to network flows. Journal of the Society for Industrial and Applied Mathematics,
10(2):260–283, 1962.

[GH64]

R. E. Gomory and T. C. Hu. Synthesis of a communication network. Journal of
the Society for Industrial and Applied Mathematics, 12(2):348–369, 1964.

[LSSW99] M. Labbé, R. Séguin, P. Soriano, and C. Wynants. Network synthesis with nonsimultaneous multicommodity flow requirements: bounds and heuristics. Technical report, Institut de Statistique et de Recherche Opérationnelle, Université
Libre de Bruxelles, 1999.
[Min81]

M. Minoux. Optimum synthesis of a network with non-simultaneous multicommodity flow requirements. North-Holland Mathematics Studies, 59:269–277,
1981.

17

[PB14]

N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):123–231, 2014.

[PLO07]

G. Petrou, C. Lemaréchal, and A. Ouorou. An approach to robust network design
in telecommunications. RAIRO-Operations Research, 41(4):411–426, 2007.

[PR11]

M. Poss and C. Raack. Affine recourse for the robust network design problem:
between static and dynamic routing. In Network Optimization, pages 150–155.
Springer, 2011.

18

A

Solution of the reservation update subproblem

Here we give a simple method to carry out the reservation update problem (12). This section
is mostly taken from [PB14, §6.4.1].
We can express (12) as
minimize βt + (1/2)kx − uk2
subject to xi ≤ t, i = 1, . . . K,

(16)

Here there variables are x ∈ RK and t ∈ R. (The variable x corresponds to fj in (12), and
the parameters u and β correspond to αfj + (1 − αf˜j ) + πj /ρ and pj /ρ, respectively.)
The optimality conditions are
x?i ≤ t? ,

µ?i ≥ 0,

µ?i (x?i − t? ) = 0,

(x?i − ui ) + µ?i = 0,

1T µ? = β,

where x? and t? are the optimal primal variables, and µ? is the optimal dual variable. If
xi < t, the third condition implies that µ?i = 0, and if x?i = t? , the fourth condition implies
that µ?i = ui − t? . Because µ?i ≥ 0, this implies µ?i = (ui − t? )+ . Substituting for µ?i in the
fifth condition gives
n
X
(ui − t? )+ = β,
i=1
?

We can solve this equation for t by first finding an index k such that
u[k+1] − u[k] ≤ (u[k] − β)/k ≤ u[k] − u[k−1]
where u[i] is the sum of the i largest elements of u. This can be done efficiently by first
sorting the elements of u, then computing the cumulative sums in the descending order until
the left inequality in the above holds. Note that in this way there is no need to check the
second inequality, as it will always hold. (The number of computational operations required
to sort u is K log K, making this the most computationally expensive step of the solution.)
With this index k computed, we have
t? = (u[k] − β)/k.
We then recover x? as x?i = min{t? , ui }.

B

Iterate properties

Here we prove that, for any l = 1, 2, . . . , if π (k) are the columns of Π(l), optimality condition
1 is satisfied, i.e., Π(l)1 = p and Π(l) ≥ 0. Additionally, if f (k) are taken to be the columns
of F (l), then condition 2 is also satisfied, i.e.,
m X
K
X

Π(l)jk F̃ (l)jk = pT max F̃ (l).

j=1 k=1

19

To see this, we first note that, by defining h : Rm×K → R such that h(F ) = pT max F ,
then the subdifferential of h at F is
∂h(F ) = {Π ∈ Rm×K | Π1 = p, Π ≥ 0, Πjk > 0 =⇒ (max F )j = Fjk },
i.e., it is the set of matrices whose columns are scenario prices and whose nonzero elements
coincide with the elements of F that are equal to the row-wise maximum value.
Now, from the reservation update equation, we have
F̃ (l + 1) = argmin pT max(F̃ ) + (ρ/2)kαF (l + 1) + (1 − α)F̃ (l) + (1/ρ)Π(l) − F̃ k2F .
F̃

The optimality conditions are


Π(l) + ρ − F̃ (l + 1) + αF (l + 1) + (1 − α)F̃ (l) ∈ ∂h F̃ (l + 1) .
Using the price update equation, this is

Π(l + 1) ∈ ∂h F̃ (l + 1) .
This shows that the columns of Π(l + 1) are indeed valid scenario prices. Additionally,
we have (dropping the indices (l + 1) for Π(l + 1) and F̃ (l + 1))
m X
K
X
j=1 k=1

Πjk F̃jk =

m
X
j=1

(max F̃ )j

K
X

Πjk =

k=1

m
X

(max F̃ )j pj = pT max F̃

j=1

The first step follows from the fact that the only terms for which Π(l + 1) are positive, and
thus the only terms that contribute to the sum are those for which the corresponding element
in F̃ (l + 1) is equal to the maximum in that row. The second step follows from the fact that
Π(l + 1)1 = p.

20

