Two Large-Scale Network Design
Problems
byOlga Raskina

Submitted in partial fulfillment of therequirements for the degreeof Doctor of Philosophy
in the Graduate School of Arts and Sciences

COLUMBIA UNIVERSITY

2003

ABSTRACTTwo Large-Scale Network Design Problems

In this work we consider two large-scale network design problems arising in network
routing and provisioning. The first problem combines long-term business and design decisions on large-scale telecommunication networks. The problem integrates
different elements of network business planning such as capacity expansion, routing
and protection and equipment maintenance as well as a complicated pricing model;
which altogether make the optimization particularly hard. We propose a scalable
optimization algorithm that exploits the structure of the model and produces a fast,high-quality solution. In the second part of this work we consider the maximum
concurrent flow problem and analyze one of the first special-purpose approaches,
proposed in 1971 by Fratta, Gerla and Kleinrock, and show that it yields a fully
polynomial-time approximation scheme.

Contents

1 Introduction

1

2 Business Planning of Optical Networks

8

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2

Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2.1

Price-demand relationship . . . . . . . . . . . . . . . . . . . .

14

2.2.2

Routing and Protection. . . . . . . . . . . . . . . . . . . . .

16

2.2.3

Capacity Planning . . . . . . . . . . . . . . . . . . . . . . . .

21

2.2.4

Cutting plane algorithms for constrained nonlinear optimization 22

2.2.5

Newton’s method for unconstrained nonlinear optimization . .

24

2.2.6

Standard network design problem . . . . . . . . . . . . . . . .

27

The core model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

2.3.1

Model description . . . . . . . . . . . . . . . . . . . . . . . . .

31

Algorithm idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.3

2.4

i

2.5

Difficulty of the core model. . . . . . . . . . . . . . . . . . . . . . .

37

2.6

Protection models . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

2.6.1

Shared protection on single demand (SD) . . . . . . . . . . . .

42

2.6.2

Shared protection across demands (MD) . . . . . . . . . . . .

44

2.6.3

Flows on the paths . . . . . . . . . . . . . . . . . . . . . . . .

47

Algorithm outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

2.7.1

The nonlinear optimization algorithm . . . . . . . . . . . . . .

54

The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

2.7

2.8

2.8.1

Shared protection on single demand with fixed percentages (SDP) 56

2.8.2

Shared protection across demands with fixed percentages (MDP) 65

2.8.3

Shared protection on single demand without percentages (SDN) 66

2.8.4

Shared protection across demands without percentages (MDN)

66

2.8.5

Bounds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

Starting point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

2.9.1

Fixed percentages . . . . . . . . . . . . . . . . . . . . . . . . .

68

2.9.2

No percentages . . . . . . . . . . . . . . . . . . . . . . . . . .

72

2.10 Implementation and numerical results . . . . . . . . . . . . . . . . . .

75

2.10.1 Real-world problems . . . . . . . . . . . . . . . . . . . . . . .

75

2.10.2 Generated problems. . . . . . . . . . . . . . . . . . . . . . .

85

2.11 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

2.9

ii

3 Maximum Concurrent Flow Problem

89

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

3.2

Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

3.3

The central idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

3.4

The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

3.4.1

Analysis of Algorithm FD . . . . . . . . . . . . . . . . . . . . 101

3.4.2

Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . 106

3.4.3

Choosing the stepsize . . . . . . . . . . . . . . . . . . . . . . . 109

3.4.4

Comments on Algorithm FD . . . . . . . . . . . . . . . . . . . 110

3.5Single Frank-Wolfe iterate analysis . . . . . . . . . . . . . . . . . . . 112

3.6

Future research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
3.6.1

Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

3.6.2

Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

3.6.3

Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

iii

List of Figures
2.1

Relationship between demands, prices and capacities . . . . . . . . .

13

2.2

Protection example . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

2.3

An example of the region set . . . . . . . . . . . . . . . . . . . . . . .

59

2.4

Original network for SDP model . . . . . . . . . . . . . . . . . . . . .

69

2.5

Compressed network Ḡ for SDP model . . . . . . . . . . . . . . . . .

69

2.6

Original graph for the MD model . . . . . . . . . . . . . . . . . . . .

71

2.7

Compressed network Ḡ for SDN model . . . . . . . . . . . . . . . . .

72

3.1

Potential function decrease in θ-neighborhood . . . . . . . . . . . . . 121

iv

List of Tables
2.1

Statistics on real-life models . . . . . . . . . . . . . . . . . . . . . . .

38

2.2

Performance comparison (LOQO and SNOPT)- running time (sec) . .

39

2.3

Larger networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

2.4

CPLEX performance - running time (sec) . . . . . . . . . . . . . . . .

41

2.5

Statistics on real-life networks . . . . . . . . . . . . . . . . . . . . . .

77

2.6

Data parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

2.7

Capacity costs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

2.8

Performance comparison - SDP model. . . . . . . . . . . . . . . . .

79

2.9

Average running time (sec) - SDP model . . . . . . . . . . . . . . . .

80

2.10 Heuristic performance - SDP model . . . . . . . . . . . . . . . . . . .

80

2.11 Performance summary - SDN, MDP, MDN models . . . . . . . . . . .

81

2.12 Running time (sec) summary - SDN, MDP, MDN models . . . . . . .

82

2.13 The effect of fixing percentages . . . . . . . . . . . . . . . . . . . . .

82

2.14 Model comparison (%) . . . . . . . . . . . . . . . . . . . . . . . . . .

82v

2.15 Performance comparison (LOQO and SNOPT) - Average optimality gap (%) 84
2.16 Performance comparison (LOQO and SNOPT)- running time (sec) . .

84

2.17 Statistics on the generated networks . . . . . . . . . . . . . . . . . . .

85

2.18 Performance comparison (CPLEX) - running time (sec) . . . . . . . .

86

2.19 Performance summary - Optimality gap (%) . . . . . . . . . . . . . .

87

2.20 Average running time (sec) summary . . . . . . . . . . . . . . . . . .

87vi

ACKNOWLEDGMENTSFirst and foremost I must express my deep an sincere gratitude to my advisor, Professor Daniel Bienstock for his guidance, support and patience, as well as his keensense of humor that led me through the years of my study at Columbia University.
His help and advice, both in technical and personal issues is something I am going toremember and value for the rest of my life.
I would like to thank the members of the committee: Donald Goldfarb, Clifford
Stein, Edward Coffman and Iraj Saniee for valuable remarks that helped to improve
this work.
I am grateful to my mentors at Bell Labs, Iraj Saniee and Qiong Wang for theirassistance and support without which a big part of this dissertation would not havebeen possible. I also thank Ramesh Nagarajan, my mentor at Lucent Technologies,
for my first exposure to the world outside of academia and for his help and patiencethat made my first professional experience valuable and unforgettable. I would liketo thank the wonderful staff of IEOR department office, who had always provided
qualified help with all kinds of problems.
My special thank you to my husband Anton, who has always been there for me
and put up with more than I could ask for during this whole ordeal.
I also thank all my friends for making my student years at Columbia memorablevii

and colorful. I choose not to list all their names, since I do not want to miss any ofthem.
Last, but not least I want to thank my mother who has been the source of mystrength and inspiration through all the years.

viii

CHAPTER 1. INTRODUCTION

1

Chapter 1
Introduction
In this thesis we deal with two network planning and operation problems. The
first problem concerns multiperiod business planning of a telecommunication network
while the second is a version of the maximum throughput routing problem.
Evolving technologies in telecommunication networks offer the capability of verylarge wide-area networks with throughput on the order of terabits (1012 bit) per
second. These technologies are considered very promising for the next generationof transport networks, as they can satisfy the growing demand for bandwidth thatis caused primarily by the explosive growth of web-related services offered over theInternet. The challenge is to react quickly to these increasing bandwidth requirements
while maintaining a reliable service. Networks should be designed and operated so
as to provide adequate capacity in the areas with growing demand without reducing

CHAPTER 1. INTRODUCTION

2

revenue. Models of this kind give rise to very difficult optimization problems of
constantly increasing size and complexity due to both the growth of network size
and traffic throughput, and to increasing concerns about data privacy, reliability,
and security. These decision and optimization problems generally require approaches
that dynamically update static solutions based on the current data and the stateof the network. Despite continuing advances, computer technology is still unable tokeep up with the growth in model and data sizes. At the same time, it is evident
that empirical and intuitive heuristics can no longer fulfill the need for fast andprovably good solutions. Clearly, there is a demand for sophisticated algorithms
employing state-of-the art mathematical and computing techniques to answer thegrowing modeling and optimization needs.
In this work we provide scalable and rigorous algorithms for two difficult largescale network design problems. The first problem deals with long-term planning andoperation of an optical network. The major aspects of such planning involve decisions
spread over several time periods, which concern demand pricing, routing and protection as well as network design and capacity expansion. Each of these components
gives rise to a difficult optimization problem in its own right; thus the integration of
all these components poses a significant challenge to current optimization techniques.
We consider network operators who contemplate capacity planning in their networks
for a period that spans approximately a decade. A key feature in our model is the ex-

CHAPTER 1. INTRODUCTION

3

pectation that, during such a period, several generations of transmission technologies
will emerge, technologies that will outperform prior ones by improved measures ofthroughput such as the maximum number of wavelengths, capacity per wavelength,
or both.
However, there can be incentives to delay the deployment of new technology to
exploit savings from future cost reductions. A dynamic technology environment of
this kind immediately poses two interesting problems. First, there is one of timing:
when should the operator start to deploy new systems and phase out old technologies?
Second, there is one of sizing: what capacity that be deployed on each link during
every time period? These two questions have been referred to as the technology
substitution problem and capacity expansion problem, respectively. In our model,
the two problems are closely interrelated and will be considered together.
The network planning problem above has been discussed extensively in the literature. In many cases, it has been cast as an optimization problem of choosing thecapacity configuration strategy that satisfies forecast demands (which are given as
inputs) at minimum cost (which includes both the cost of deploying new capacities
and the cost of operating old systems); see for example, [5]. While min-cost modelingrationalizes decision-making from an engineering perspective, more comprehensive
approaches have been taken in recent studies to integrate capacity planning into theoverall business optimization strategy [37], [20]. Under the generalized framework,

CHAPTER 1. INTRODUCTION

4

instead of fixing forecast demands, the price-demand relationship is taken into consideration and demands are treated as flexible quantities that can be controlled byprices. As decision variables, prices play a dual role. On the one hand prices determine the revenue for each unit of service rendered. On the other hand demands
cannot deployed capacity. To put it differently, capacity planning is driven not by
the need to meet a fixed demand target, but by the desire of generating more profit.
The objective is to choose the prices of network services and network capacity simultaneously, so as to maximize the overall net profit generated by the network. In a
multiple period model such as ours, the net profit is referred to as Net Present Value(NPV), and is defined to be the total revenue minus the total cost over the planninghorizon, using an appropriate discount factor.
To meet the challenge of solving the generalized multi-period price setting andcapacity planning problem, our work develops a set of efficient, scalable algorithms.
The overall problem is formulated as a nonlinear optimization model with a concave
objective function and linear constraints. While a concave maximization problem is(in theory) solvable using standard techniques, the problem sizes arising in realistic
settings makes such models essentially impossible to solve by conventional methods,
within anything even remotely resembling a reasonable time frame. The algorithm
described in this work provides a considerable improvement in terms of both optimality error and computing time, when compared to state-of-the-art general-purpose

CHAPTER 1. INTRODUCTION

5solvers.
The second problem considered in this work is a version of the maximum throughput routing problem, namely the canonical version, called the maximum concurrent
flow problem. The maximum concurrent flow problem frequently arises in practicalapplications and has received much attention because of its difficulty.
Many telecommunication network design problems contain a multicommodity flow
problem as one of the components. Generally, in a multicommodity flow problem weneed to route a set of demands (or traffic requirements) over a given network withprespecified joint capacities on its edges. The minimization version of the problem
requires a minimum cost routing, while the feasibility version merely tries to finda routing of the demands such that no edge capacity is violated. It is well known
that even the feasibility version is fairly difficult. When the problem turns out to be
infeasible it can be significant to determine the extent of the infeasibility. In thatcase we want to find the maximum throughput of the network, that is the maximum
common fraction of the demands that can be routed simultaneously without exceeding
any of the network capacities. This is called the maximum concurrent flow problem.
This problem arises in telecommunications in many different variations (directedor undirected edges, restrictions on the routings, etc.). Further, many algorithms for
network design rely on solving maximum concurrent flow problems in order to find afeasible routing.

CHAPTER 1. INTRODUCTION

6

What makes this problem especially noteworthy is that it gives rise to extremely
difficult linear programs – instances of a size and type relevant to applications often
prove beyond the reach of state-of-the-art linear programming codes, despite recent
progress in computer technologies. For this reason there has been sustained research
devoted to the development of approximation algorithms. The use of approximation isjustified, as a fast, partially accurate result is much more desirable than an algorithm
that provides an exact solution but requires an excessive amount of time.
In this work we analyze the convergence of one of the first approximation algorithms developed for this type of problem as well as some enhancements of thealgorithm that help to improve performance and simplify implementation. In addition we offer a further refinement of this algorithmic approach as a direction for future
research.
This thesis is organized as follows.
In Chapter 2 we present our work on the first problem – multiperiod business planningof a telecommunication network. Section 2.2 reviews several essential optimization
concepts used in our work, as well as some essential components of telecommunication network design, such as capacity planning and routing and protection. We also
provide a brief review of a recently developed demand-price model used in our design.
Section 2.3 introduces the basic underlying model and discusses the major differences
between this model and classical network design problems. Section 2.4 gives a broad

CHAPTER 1. INTRODUCTION

7

overview of the optimization algorithms and the resulting changes in the model.
Section 2.5 analysis the complexity of the problem and the performance of some wellknown commercial optimizers. In Section 2.6 we expand one of the key concept of
the problem – the survivability of the network, and give precise formulations of thefour different optimization models resulting from applying different protection andintegration techniques. Section 2.7 outlines the stages of the algorithm and presents
the landmarks of each step. Section 2.8 provides the detailed description of the algorithm and the justification of the upper bound technique. Section 2.9 describes theheuristic used to obtain the initial solution. Section 2.10 presents the summary andthe analysis of the extensive experiments conducted with real-life networks and also
large generated networks. Section 2.11 concludes with a discussion of the results.
Chapter 3 proceeds with the second part of this work – the maximum concurrent
flow problem. In Section 3.1 we introduce the problem and give an overview of theprevious work. In Section 3.2 we formulate the problem and introduce some basic
notation used in our work. Section 3.3 outlines the algorithm and explains someessential underlying ideas. Section 3.4 describes the algorithm in detail and thenproceeds with the complexity analysis. In Section 3.5 we describe and analyze aslightly different version of the algorithm. Finally Section 3.6 presents some ideas for
future research.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

8Chapter 2
Business Planning of Optical
Networks

2.1

Introduction

Strategic planning for an optical network involves sophisticated decision making inmany respects. While a major component of this task is to minimize cost by optimizing routing and protection of network traffic and by determining network capacity, thebusiness impact of these schemes should also be put into a broad perspective. Specifically, cost reductions resulting from routing and capacity decisions can be translated
into pricing advantages that enable the carrier to generate more demand and revenue.
Therefore, optimizing pricing and demand should be integrated into the optimization

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

9model. To maximize the total value to the carrier, routing, capacity, and pricing
decisions should not be examined individually, but rather be incorporated into an
integrated business model and optimized systematically.
In this work we develop efficient algorithms for solving capacity planning problems
for long-distance optical transport networks. In our model, a network is formulated
as a set of nodes, representing cities and metropolitan areas, connected by links,
representing physical routes (right-of-way) owned by the network operator. Communications between a node pair are carried by transmission systems deployed on links
that connect the pair. A transmission system is composed of a collection of optical
wavelengths, where each wavelength is a carrier of information with a fixed amountof capacity. The capacity of a system can be set at different levels by configuring different numbers of wavelengths, up to the maximum. Therefore, the cost of deploying
new capacity involves a fixed investment of installing the initial system and a variable
cost of adding wavelengths. Each system or wavelength in use is also associated witha recurring maintenance cost, also defined as the operating expense.
It is natural to plan the capacity expansion over several time periods. During
such a span of time, several generations of transmission technologies will emerge.
A later technology always outperforms previous ones by increasing the maximum
number of wavelengths, capacity per wavelength, or both. Though a newer system
may have a higher fixed cost and/or variable cost per wavelength, the magnitude

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

10of capacity improvement usually will outpace the corresponding cost increase. As
a result, a fully-loaded new system usually has a lower deployment cost per unitof capacity than the older ones, which makes it an attractive candidate for new
deployment. Nevertheless, after a new technology becomes available, the deployment
cost of systems and wavelengths decreases over time due to a learning effect.
Implementing the value-based approach in large-scale optical networks poses significant challenges to optimization techniques. In comparison with the min-cost capacity planning model, the new approach introduces a set of pricing variables, whose
number is on the order of n2 T , where T is the number of planning periods and n is thenumber of nodes in the network (and, consequently, the number of node pairs between
which the price of network service needs to be determined grows proportional to n2 ).
Furthermore, since revenue equals price times demand and demand is a function ofprice, the pricing variables are nonlinear components of the revenue in the objectivefunction. As a result when planning a national or international optical transport
networks (which usually contain 50 − 100 nodes) over a 10-period horizon, applying
the generalized model can easily lead to a nonlinear model with tens of thousands, oreven millions of variables.
The network aspect of the problem brings yet another set of complexities for
optimization. To ensure reliable communications in the presence of possible node andlink failure, optical transport networks are designed to have multiple non-overlapping

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

11

paths between every node pair. The non-uniqueness of routes leads to questions about
which routes should be used to carry demand, and in what quantity. Additional
planning is required to specify the rerouting of traffic when a network failure occurs.
There are several different restoration models that can be employed in this framework.
For example, we can reroute all or some of the commodities to utilize the existing
network capacity. Or, we can only allow to reroute the commodities that used thedisconnected link or node. Or, some resilient capacities are reserved throughout the
network and used only to restore the some or all of the affected traffic. See [7], [31],[1], [2] for the discussions of different models and solution approaches.
These routing and protection decisions need to be modeled in the overall optimization framework. As it will be shown in the following sections, the efficient routing
and protection schemes lead to vastly different formulations of the overall optimization model. Formulations of some complicated schemes require the introduction ofO(n2 T ) additional flow variables. The resulting increase in the problem size makes it
necessary to develop efficient algorithms.
We formulate the general problem as a large-scale non-linear integer model. The
objective is to maximize the carrier’s net present value over an extended period oftime with constraints on network topology (right-of-way), technology roadmap, andprice-demand relationship. In this model, not only demands and capacities but also
end-to-end demand allocations and prices are determined by period, with price being

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

12

a nonlinear function of the demand following a certain model to be discussed later. In
addition, protection channels are to be allocated according to the specified protection
model. Our task is to determine the amount of traffic for each commodity in the
network at each time period and also the amount of fibers and wavelengths required.
Once the equipment is bought we have an option of maintaining it for any numberof periods at some increasing fraction of its initial cost. Furthermore, it can also be
retired at any time at no additional charge and once retired will no longer be available
at any consecutive period.
The modeling of demand as a nonlinear function of price in the context of telecommunication systems was first proposed in [20]. In that work, the solution methodology
could only handle small networks (5 − 10 nodes) which in addition had ring topology.
Further, experimentation showed that the methodology in [20] would not scale. In
this thesis we expand that model to general networks and create a more precise modelof the maintenance cost structure and protection capabilities. In [20] the maintenancecost of any equipment depends only on the current time period. In other words, thecost of maintaining a mile of fiber during any period will be the same for a year old
fiber and for a ten years old fiber. However, this is not the case in real life. Equipment
tends to wear and hence older fibers will usually require higher maintenance expenses.
In our work we view this expenses not as a function of current period, but rather asa function of the current service time. In all fairness, [20] deals mainly with build-

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

13

ing a precise economic model and developing accurate parameters. We instead takethese parameters as input constants and concentrate on the optimization techniques,suitable for large-scale problems.

2.2

Preliminaries

In this section, we discuss several key concepts of the generalized capacity planningmodel. The scope of our study is outlined in Figure 2.1 and is explained as follows.

Figure 2.1: Relationship between demands, prices and capacities

Based on the price-demand relationship, a given choice of prices determines corresponding amounts of demands and revenues. The routing and protection decisions
then map generated demands to link capacity requirements. These requirements are
satisfied by the deployment of transmission systems, whose availability, cost, andcapacity are specified by the technology roadmap. Pricing decisions and capacity
planning are thus interrelated, as demands generated cannot exceed available capac-

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

14

ity.
These three decisions (pricing, routing/protection, and capacity planning) areembedded in the overall optimization framework that seeks to maximize the netpresent value.
In the following, we discuss each of the three decisions as well as some general
optimization techniques and models used in our work.

2.2.1

Price-demand relationshipA modern telecommunication network is determined by continously changing costs
and a rapid reflection of the changes on the amount of traffic through the network. To
reflect these properties in modeling the demands and prices for the optical networks
we first need to introduce the elasticity notion. Let d be the demand and p be theprice per unit of the demand. The demand elasticity  for a given time period isdefined as the negative ratio of the relative change in the demand to the relativechange in the price during that period. That is

=−

∆d/d
∆p/p

(2.2.1)

Here ∆d is the change in demand and ∆p is the relative change in price. Commonly, the demand increases as the price decreases and in that case elasticity  ispositive.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

15

Writing∆d/d = −∆pp(2.2.2)

we can see how different values of  affect the demand-price relationship. If  > 1
then any decrease of price gives a larger increase in demand, whereas if  < 1 thesame reduction in price leads to a smaller increase in demand. For example, if theprice drops by 10% in the first case the demand grows by more than 10% and in thesecond case by less than 10%.
Taking the limit of (2.2.2) we get

d0 /d
=− 0
p /p

(2.2.3)

Let revenue be the product of the demand and price, that is Rev = pd, then from
(2.2.3) we get the expression for the change in revenue, that isRev 0
 − 1 d0p0=(
) = −( − 1)
Rev

d
p

(2.2.4)This relationship demonstrates the effect of the elasticity value on the pricedemand dependence. In particular, for any value  > 1 the decrease in price causes
an increase in the revenue, on the other hand for  < 1 a decrease in price leads to a
reduction of the revenue.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

16

For this reason in our work we assume that the elasticity is a constant greaterthan 1. Then solving (2.2.3) we obtaind = ap−(2.2.5)where a is a positive constant. Using (2.2.5) the revenue can be written as a functionof demand only, giving
Rev = a1/ d1−1/

(2.2.6)

Under the assumption that  > 1, the revenue Rev is a monotonically increasing
and concave function of d. In other words, the network operator derives a higher
revenue by reducing the price and therefore increasing demand, but with diminishingreturn to demand size. This price-demand model was originally developed by Lanning,
et. al. for a simplified network model. More details on the demand-price relationshipand the computation of the elasticity can be found in [20].

2.2.2

Routing and Protection

The reliability or survivability of a network has received a considerable amount of attention in the network design research and is a vital feature of virtually any telecommunication system.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

17

Any fiber cut or other hardware failure even at a single point can result in themalfunction of the whole system, unless a fast and reliable restoration of the traffic isavailable. This restoration usually involves rerouting the demands on an alternative
link or node-disjoint path. The routing and protection, which can be implemented inmany different ways, map end-to-end demands to link capacity requirement and canbe implemented in several different ways.
There are two most commonly used restoration strategies – link restoration andpath restoration.
• Path protection/restoration
In path protection each commodity pair has a statically reserved backup path,that is link- or node-disjoint with all working paths. In path restoration eachcommodity dynamically determines a failed end-to-end path and rerouts the
traffic on a backup path.
– Dedicated-path protection
In this approach a set of disjoint paths is specified for each node pair, andon each of these paths, some amount of bandwidth is allocated exclusively
to the pair. The allocated bandwidth satisfies the requirement that if any
one of these paths is disconnected by a link or node failure, there will still
be enough capacity on remaining paths specified for this pair to carry all
traffic. The scheme provides full protection for every commodity against

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

18

single path failure, is simple for restoring traffic when failure occurs, but
not efficient in bandwidth usage. This model is also referred to as a shared
protection on single demand since the same spare capacity allocated for thecommodity can be used to restore traffic from any of the paths connecting
this pair.
– Shared-path protection
Under this approach, each node pair is assigned a set of working paths
and a set of protection paths. Bandwidth on protection paths can beshared among different node pairs that do not have common links or nodes
on their working paths. With careful planning, shared protection uses
much less redundant bandwidth, as compared to dedicated protection, tofully restore traffic in case of single link failure [14]. Nevertheless, rather
than maintaining the same routing during the failure, as is the case withdedicated protection, shared protection usually requires rerouting of traffic
based on the location of failure. This approach also improves bandwidth
efficiency at the cost of the restoration ability and make the restoration
mechanism more complicated.
– Path restoration:
In this scheme all commodities with traffic on the failed link participate ina distributed algorithm that dynamically discovers a restoration route forCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

19each commodity. If for some commodity no feasible route is found, thisconnection is blocked.
• Link protection/restoration
Link protection/restoration is the scheme by which traffic is rerouted onlyaround the failed link, regardless of the particular demand pairs that contributeto the flow on that link. The corresponding link backup paths are reserved inadvance. The end nodes of the link dynamically determine a failure and reroute
all traffic on the protection path. Similarly, in this approach there are threeprotection schemes.
– Dedicated-link protection
Every working link of the network has a predetermined backup path aroundthis link, on with spare capacity allocated exclusively to the traffic on thelink– Shared-link protection
The capacity on the restoration path around the link can be used by otherbackup path.
– Link restoration
The end nodes of the link discover the restoration path dynamically and
if no feasible path found all connections using the link are blocked.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

20

Generally, path restoration requires more complicated hardware that maintains
the set of restoration paths for each commodity pair and also the information of thesource and destination of each unit of traffic through every link. On the other hand
it has been empirically shown that path restoration requires less total spare capacity [40]. For this reason in this work we consider only path protection/restoration
technique. In addition, in our design we require a one hundred percent path restoration strategy to provide a high degree of network survivability. Thus in this work
we concentrate on only two approaches, dedicated-path protection and shared-path
protection.
We consider two different ways of incorporating routing and protection into ourmodel. The status quo approach is to treat the distribution of traffic on different
paths as a lower-level, separate problem from the business optimization carried out
by our model. Using planning tool specialized for routing and protection design [14],we can obtain the percentages of demand to be carried on each path for every nodepair based on some ”seed demand”. We then apply these fixed percentages to ourmodel as input parameters. We call this method static traffic distribution. The
second method does not use any pre-processor to specify the demand fractions, andviews traffic distribution as decision variables. Consequently, traffic distribution mayvary with time as different amounts of demands are generated in each period. We
define this method as dynamic traffic distribution.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

21

In the next section, we first develop a base formulation that assumes dedicatedprotection and static traffic distribution. We then present formulations of shared
protection and dynamic traffic distribution as extensions to the base model.

2.2.3

Capacity Planning

Capacity planning involves decisions on installation and retirement of transmission
equipment of different technologies over time. The major decisions to be made in thissetting are the amounts of equipment bought, maintained, and retired at every time
period. In our model we assume that any amount of the existing equipment can be
retired at any time at no additional cost, and once retired it is no longer available atany consequent time period. The equipment of different technologies is characterized
by different acquisition costs and amount of wavelengths or capacity.
The cost of deploying a technology depends on technology type, and the timing of
installation. Equipment of a later technology type usually has a better capacity/cost
ratio than older equipment. The cost of given technology decreases over time, and fora given technology, the deployment cost is lower at a later installation date. Furthermore, an optical transmission system includes amplifiers and regenerators that need
to be deployed on the intermediate points on the link between transmission terminals,
and the distance between two adjacent points is approximately the same for all links.
As a result, links of different lengths require different numbers of these components,

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

22

and thus the deployment costs are link-dependent. Likewise, the recurring maintenance costs per period also depends on technology, time, and link. In our models,
we assume that for a given piece of equipment, the deployment cost of a technology
decreases at a fixed rate over time after it becomes available and the maintenance costof a system or wavelength is a fraction of its initial installation cost, and increases
over time at a constant rate.

2.2.4

Cutting plane algorithms for constrained nonlinear optimization

Let f (x) ∈ C1 be a convex function and consider a problem
min f (x)x∈Qfor some convex set Q ∈ Rn . The main idea of the cutting plane method to solve thisproblem is to generate a series of hyperplanes which approximate f (x) from below
and provide a piecewise liner approximation to f (x). This approximation is then
optimized over Q, a new approximation point is generated, and the piecewise linear
approximation is updated using a set of rules specific for each particular method.
Let x1 , . . . xk ∈ Q be a set of iterates with function values f (xi ) and gradientvalues ∇f (xi ). Then the cutting plane model generates a sequencefi (x) = max{f (xj ) + ∇f (xj )T (x − xj )|j = 1 . . . i}

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

23

Clearly fi (x) ≤ f (x) ∀x ∈ Q. A number of ideas for generating the series of iterates x1 , . . . xk ∈ Q have been proposed over the years. Here we review several basic
schemes.
Classical cutting-plane algorithm
This is the simplest algorithm in which the next iterate xi+1 is taken to be argmin{fi (x)|x ∈Q}. It is known to have a very slow convergence [29].
Line search algorithm
Let x∗i = argmin{fi (x)|x ∈ Q} and set
xi+1 = αi xi + (1 − αi )x∗i
where
αi = arg min{f (λxi + (1 − λ)x∗i )|0 ≤ λ < 1}This method is also rarely used because in a general case the computation of αiparameter is very computationally costly. However, for a simple f (x) and a smallenough n this method proves to be quite effective.
Bundle algorithm
Different versions of this algorithm are described in [19]. The general idea is to set
2
xi+1 = argmin{fi (x) + 0.5ui (x − x+i ) }where x+i is chosen among {x1 . . . xi } and ui is the penalty parameter. The idea of
this method is to update x+
i every time there is an adequate decrease in f (xi+1 ) andCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

24

keep x+i unchanged otherwise.
Level algorithm
This model was proposed by Lemarechal et al.[21]At every iteration of the algorithm we define the level parameter li for the current
function value and the next iterate is computed asxi+1 = argmin{0.5(x − xi )2 |x ∈ Q, fi (x) ≤ li }Let f∗ (i) = minx∈Q fi (x) and f ∗ (i) = minj≤ı f (xj ). Define the ith gap as ∆(i) =f ∗ (i) − f∗ (i), then for some 0 < λ < 1 set

li = λf ∗ (i) + (1 − λ)f∗ (i) = f∗ (i) + λ∆(i)

2.2.5

Newton’s method for unconstrained nonlinear optimization

In this section we give a brief description of a second-order algorithm for unconstrained
minimization of a general convex function f (x).
This method is based on the Taylor series approximation of f (x). Let f (x) ∈ C2
for some f : Rn → R and x ∈ Rn . Let p ∈ Rn . Then1f (x + p) = f (x) + g(x)T p + pT G(x)p + o(kpk2 )2

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

25

where
g(x) = ∇f (x)andG(x) = ∇2 f (x)Denote1Q(p) = g(x)T p + pT G(x)p2
Newton’s method proceeds by iteratively optimizing the quadratic approximation
Q(p) to f (x + p) at f (x). We have∇Q(p) = g(x) + G(x)T p
and assuming that f (x) is strictly convex, that is G(x) is nonsingular we get∇Q(p) = 0 iff p = −G(x)−1 g(x)(2.2.7)

Thus the Newton’s method iterates by solving (2.2.7) and updating x ← x + p until
kg(x)k ≤  for some predefined tolerance . Note that in order to solve (2.2.7) wefirst need to invert the Hessian matrix G(x). A commonly used method for implicitly
inverting a Hessian is Cholesky factorization. This method is used to represent a
symmetric positive definite matrix A as
A = RT R

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

26

where R is called Cholessky factor and is a nonsingular upper triangular matrix. The
Cholesky factor can for example be computed using elimination. Given the Cholesky
factor R the solution to a system of linear equations Ax = b can be found by solving
two triangular systems

RT y = b

andRx = yComputation of a Cholesky factor requires approximately half the work needed toperform Gaussian elimination. Moreover, in practical application there are numerous
methods to further speed up the computations. See [3] and [12] for the discussion ofthese methods.
Newton’s method has quadratic convergence when applied to a convex problem.
In fact, it is known that, assuming convexity,nkxn − x∗ k ≤ O(kx0 − x∗ k2 )where x∗ is the optimum of f (x), xn is the current iterate and x0 is the starting
point. See [3] for the details of the proof. Thus the actual rate of convergence heavily
depends on the initial error and whence obtaining a good starting point is essential
for this method. The Newton’s Method can also be used for optimizing a linearly
constrained convex function using a set of sequential unconstrained optimizations.
This approach will be described in more details in Sections 2.7 and 2.8.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.2.6

27

Standard network design problem

The network design or minimum-cost capacity installation problem is the problem of
installing an integer amount of capacities on the links of a given graph such that agiven set of demands can be routed simultaneously.
Consider a graph G = (V, A) where V is the set of nodes and A is the set of arcs
and a let d = {dk } be the given set of demands (commodities) between the nodes of
G. Let K be the number of distinct commodities and for every 1 ≤ k ≤ K let ik andjk be the source and destination of demand dk which in general can be fractional.
The capacity on links can be installed in batches of different size.
Let U be the set of different batch types and M (u) is the size (number of capacity
units) for each type u ∈ U . The cost of installing a batch of capacity of type u onlink l ∈ A is cl (u). The batches are installed in integer amounts and any combinationof batch types is allowed on any link. The goal is to install the batches of capacity
at minimum cost such that all demands are met.
For each commodity k let Pk = {phk , 1 ≤ h ≤ |Pk |} be the set of all paths from ik
to jk , f (phk ) be the amount of flow of commodity dk on path phk and let yl (u) be thenumber of batches of type u installed on link l. The path formulation of the minimum
cost capacity installation problem is thenCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKSmin

X

cl (u)yl (u)Xf (phk ) = dk28

(2.2.8)l,us.t.

∀1 ≤ h ≤ K

hXk,: l∈phkf (phk ) ≤X

M (u)yl (u)uyl (u) ∈ N+∀l ∈ A∀l ∈ A ∀u ∈ U

It is known that the network design problem is NP-hard even for the case when
|U | = 1 and the generic problem is strongly NP-hard [11] as it contains the set coverproblem. This problem has been studied extensively over many years. Magnantiand Mirchandani [25] and Magnanti et al. [26] proposed cut-set inequalities for asingle commodity design with multiple batch types where the batch sizes are integer
multiples of some basic capacity unit. Bienstock and Gunluk [6] generalize theseinequalities to flow-cut-set inequalities for multicommodity design with two batch
types, where the size of one batch is a integer multiple of the size of the other batch.
Gunluk [10] considers multicommodity problem with two batch types and existing
capacity. Chopra et al. [9] proposed inequalities for a single commodity problem
on a directed network with two batch types. Pochet and Wolsey [32] analyzed asingle link network with arbitrary many batch types and Bienstock et al. [5] study amulticommodity design with a single batch type.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.3

29

The core model

Our basic model differs from the standard network design problem in two fundamental
aspects. The first and foremost difference is the way the demand values are defined.
Unlike the basic problem, in our model the demands are not given in advance. Insteadfor every pair of nodes we introduce a variable representing the amount of traffic
between these nodes. Thus we not only minimize the cost of the equipment on thelinks, but at the same we also maximize the revenue from the traffic based on the pricemodel defined in the previous section. The second major difference is that we need tosolve the problem in the span of multiple time periods. For every time period in the
span we again wish to optimize the cost and the revenue of the network. Furthermore
at any time period we are allowed to reuse the equipment bought at earlier periods
provided we pay the price of maintaining it. We can also retire all or any part of theequipment at any time at no additional cost.
Note that even with these changes the complexity of the problem is at least thatof the original network design problem and generally the number of paths in theformulation (2.2.8) can be very large. However, a common approach to survivablenetwork design is to assume a precomputed fixed set of paths for every commodity.
This approach can be justified by several reasons. For instance, a set of paths canbe predefined by a specific protocol or service used to operate the network, or by
an outside provider. For example the OSPF protocol [28] currently employed in

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

30

Internet routing assumes a fixed (shortest for some metric) path for every pair ofnodes. Further, in practice there is an extensive and sometimes idiosyncratic setof requirements on the paths used to carry traffic. Usually those requirements are
non-quantitative and either can not be expressed explicitly as a set of constraints, or
would require a large number of binary and integer variables that would make theproblem practically intractable.
For example, to satisfy certain requirements of quality of service, the carrier needs
to minimize the processing delay incurred at the intermediate nodes. In this case, thecandidate paths may be restricted to those that have no more than a certain number ofhops between the node pair. The length restriction constraint prevents the demands
from taking unacceptably long path in terms of either the number of hops or the actual
path length, thus decreasing the probability of a failure of a path component. In
addition, to guarantee some degree of network survivability an additional set of path
constraints can be imposed. There are several commonly used restrictions utilized in
network modeling. For example, k−connectivity requires the network to contain atleast k node- or link-disjoint path for every commodity pair. Furthermore, concerns
about security and privacy may also limit the carrier’s choice of candidate paths. Seefor example [8].
For these reasons in this work we will only look at a given fixed set of paths forevery commodity. The algorithm in [14] is used to generate the paths with the required

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

31

properties, however any other path generating heuristic can also be employed. For theprotection reasons outlined in Section 2.2.2 and specified below the paths for every
commodity pair are link-disjoint.
We are now in position to formulate our core model.

2.3.1

Model description

Let G = (N, A) be a given network with |N | = n and |A| = m and let T be thenumber of discrete time periods.
For any pair of nodes (ik , jk ) ∈ N and a time period t let dtk be the amount of traffic(demand) between ik and jk at t. We will base our model on the path formulation of
the standard network design problem. Let Pk be the given fixed set of link-disjoint
paths for commodity k, |Pk | = Hk .
Thus we also introduce flow variables f t (phk ) which are the amount of commodity
dtk routed on path phk ∈ Pk , 1 ≤ h ≤ Hk .

For any link l ∈ A and a time period t, let ylt (u) ∈ N+ be the number of batchesof type u bought at t on l, and for any time period q > t, let ylt,q (u) be the amount ofbatches of type u that were bought at t and kept until q. For the consistency of thenotation we let yltt (u) = ylt (u)We are also given the following constant parameters

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

32• Discount factor ht . The discount factor at t is usually computed as the inverseof the predicted interest rate at t. The discount factor is used to compute thepresent value of future income. For a detailed description see [22]. The discountfactor is precomputed and used as an input parameter in our model.
• Buying cost ctl (u). The cost of buying a batch of type u on link l at time periodt. We assume that the buying cost is nonincreasing, that is ctl (u) ≥ ct+1
l (u) andthat the cost per unit of capacity is decreasing as the batch size increases.
• Maintenance fraction µ < 1 and maintenance rate α > 1. The cost of maintain. t
t−q−1
ing a unit of capacity bought at q until t is defined as cq,t
. Thusl = cl ∗ µ ∗ α
the maintenance cost is smaller than the buying cost at the first time period
and increases as t approaches the time horizon T .
• The batch size M (u) is taken to be the same across the network and the time
horizon for every type u.
• Demand-price elasticity tk > 1.
• Demand scaling constant Ak . This constant is usually taken to be equal tothe value of the demand when the price p = 1; thus Ak can be viewed as the. 1/tdemand potential. We define the demand (ik , jk ) coefficient at t as ctk = Ak k
Using the above constants and the price-demand relationship described in the

33

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

previous section we define the revenue at t as
Revt =X

gkt (dtk )(2.3.1)kwhere
gkt (d) = ctk d

1− 1t


k

.

The total equipment cost at t is defined as
Costt =

X

ctl (u)ylt (u) +X q,tcl (u)ylq,t(u)(2.3.2)q<tlwhich is the cost of buying capacity at t plus the cost of maintaining capacities fromthe previous periods.
Our objective is then to maximize the net present value, defined as

NP V =Xt

ht (Revt − Costt ) ≡ F (d, y).

(2.3.3)We introduce two types of constraints. The first set is the capacity constraints
that are inherited directly from the standard network design problem. We haveXk,h:l∈phktf (phk ) ≤

X
u

M (u)



dtk =

ylt (u) +X q,t yl (u) ∀l ∈ A, ∀t
q<t

X
P

f t (phk )CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

34

In addition we restrict the carry over capacity at t not to exceed the capacity at t − 1

ylq,t (u) ≤ ylq,t−1 (u) ∀l, ∀q, ∀t ≥ q + 1
Summarizing all of the above we can now give a complete formulation of the coremodel.

max F (d, y)(2.3.4)s.t.
X(k,h:l∈phkf t (phk ) ≤X

M (u)(ylt (u) +

dtk =

X

f t (phk )

u

q<t

h

2.4

X q,t
yl (u)))∀l ∈ A, ∀t∀k, ∀tylq,t (u) ≤ ylq,t−1 (u)∀ l, ∀q, ∀t ≥ q + 1

ylq,t (u) ∈ N+∀l, ∀q, ∀t ≥ q

Algorithm idea

In this section we will describe the main idea of the algorithm used to solve the above
model. The details of the algorithm will be described in Sections 2.7- 2.9.
The optimization problem described in the previous section is a concave maximization problem over the demand variables d and an integer program over the capacity

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

35

variables y.
A useful heuristic to handle this type of problems is to relax the integrality constraints first and solve the resulting continuous nonlinear program. After the continuous solution is obtained, the continuous variables are fixed at their current values andthe resulting integer problem is solved to optimality using any appropriate method
available.
Generally this approach proves to be inefficient for the small values of the integer
variables since in that case even a fraction of a unit has a significant contribution tothe objective value.
However, in case of telecommunication system the amount of data transmitted
over a network is typically very large due to the recent and continuing advances oftechnology. It is not uncommon to have tens of thousands units of demand between
a single node pair. The results of [20] confirm this assumption for ring networks andwe expect it to be the case for our model as well. As a result, in our case even a
straightforward greedy rounding will produce a solution within a fraction of the lower
bound that was obtained by solving the relaxed problem.
In a view of this we adopt the following approach to solve our model:
1. Relax the integrality constraints on the capacity variables and solve the resulting
continuous nonlinear problem.
2. Fix the values of the demand and flow variables found in Step 1 and find the

36

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

corresponding integer values of the capacity variables.
Note that Step 1 deals with the continuous version of problem (2.3.4). Clearly inthat case the optimum solution will only use the capacity batches with the smallest
per unit cost. Thus we can reduce the core problem to use only the cheapest capacity
on any link at any time period. This is possible since the maintenance fraction andthe maintenance rate are the same across all batch types, links and time periods. We
define the integrated capacity cost as
ctl = minuctl (u),M (u)and the core problem becomes

max

X
t

ht

X
k

gkt (dtk ) −X(ctl ylt +X q,t q,t 
cl y l

(2.4.1)

q<t

l

s.t.
X

k,h : l∈ph
k

X q,t
yl

∀l ∈ A, ∀tf t (phk )∀k, ∀tf t (phk ) ≤ ylt +
dtk =

X

q<t

h

ylq,t ≤ ylq,t−1∀ l, ∀q, ∀t ≥ q + 1(2.4.2)

ylq,t ≥ 0∀l, ∀q, ∀t ≥ q(2.4.3)In what follows we will be looking only at relaxed core problem (2.4.1) and its

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

37

extentions. In the subsequent section we will describe two models based on thisformulation.

2.5

Difficulty of the core model

Having described a core model in the previous section, we will now pause and consider
how one should tackle problems of this sort – a more specific model will be given in
the following section.
In this work we will develop an algorithm approach for optimization problems
derived from the core model. Our algorithms, as we will show, are fast and experimentally scale well. Our algorithms also take advantage of essentially combinatorial
features of the core model. Thus a natural question would be to ask whether standard, e.g. commercial solvers cannot already tackle our problems. In this Section wetake up this issue, to explain why the answer to this question is a clear ’no’.
Problem 2.4.1 is a hard nonlinear optimization problem with a steep objectivefunction and a large number of variables even for a moderate size network. For
example a regular network with 50 nodes and 70 links considered over 14 time periods
can result in a formulation with approximately 60000 variables for the complicatedprotection model, and a larger network with 150 nodes and 800 links over the sametime horizon will already have over one million of variables for the same model. The
large number of variables is due to the presence of a decision variable for every path

38

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Set

Nodes

Links

Vars

Const

Nonz

Paths

Paths

(min)

(max)

snet1

14

22

3584

2310

17472

2

3

snet2

38

48

12782

2940

96222

2

3

snet3

47

55

20909

5775

196623

2

4

snet4

50

64

23870

6720

223272

2

3

snet5

70

94

43680

9870

431382

2

3

Table 2.1: Statistics on real-life modelsof every possible commodity, at every time period.
In order to evaluate the difficulty of the problem we conducted experiments withseveral commercially available solvers on networks of different size. Table 2.1 shows
the statistics on the first set of problems arising from real-life data sets. Here the lasttwo columns show the minimum and maximum number of paths per node pair.
Our testing focused on two well-known nonlinear optimizers: LOQO [39] and
SNOPT [13]. Both are excellent general-purpose systems designed to handle any
nonlinearity in both objective and constraints. Table 2.2 shows the performancecomparison of the solvers on the test networks.
The empty cells indicate that the solver was unable to handle the corresponding
instance, which in most cases means no convergence in a reasonable time frame andCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Solver

snet1

snet2

snet3

snet4

snet5

LOQO

30

600

720

1000

–

SNOPT

90

–

–

–

–

39

Table 2.2: Performance comparison (LOQO and SNOPT)- running time (sec)in some cases insufficient memory for the solver to run, even on a machine with ample
resources. These tests were conducted on 336 Mhz UltraSPARC machine with 1.7Gb
of RAM. What is more, in all of the cases where the solvers were not able to converge
either no feasible solution was found, or the optimality gap was unacceptably large.
This test also indicates that although one of the solvers was able to solve theproblem on the moderate size networks, it still shows a poor performance in terms of
the running time.
Altogether it is evident that the problem cannot be handled by these two commercial optimizers and there is an apparent need for an algorithm that is effective on
our models.
Of course, it is impossible to test every available solver, but the above experiments
lead us to expect that general-purpose solvers will all be ineffective, in particular when
run on much larger problem instances than those in Table 2.1. One of the reasons
underlying our belief concerns the use of Newton’s method (or a similar second-ordermethod) to handle the highly nonlinear objective. This, in turn, requires a good

40

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Set

Nodes

Links

Vars

Const

Nonz

Paths

Paths

(min)

(max)

bnet1

100

500

7455004060005299700

3

6

bnet2

150

900

1659000

88935015340500

3

7

bnet3

200

1500

4336500

1850100

40267500

3

9

bnet4

250

2000

6746250

2852500

69618500

3

10

bnet5

300

4500

18053700 5558700 311100300

412Table 2.3: Larger networks
Cholesky factorization algorithm for large sparse matrices, but may nevertheless exact too high a computational price. To explore this difficulty, we also tested CPLEX[16] which is one of the most efficient linear, integer and quadratic solver available on
the market. In terms of nonlinear programming, CPLEX can only handle (convex)
quadratic programs – but CPLEX also has one of the fastest Cholesky implementations. In order to test CPLEX we replaced the original objective with its second
order Taylor series approximation and generated a large set of data described below.

Table 2.3 shows the statistics on the five large networks generated for this test.
In order to build the second order approximation for these tests we replaced theoriginal objective function 2.4.1 with its piece-wise linear approximation built at t

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Solver

bnet1

bnet2

bnet3

bnet4

CPLEX

2341

7982

14324 18943

41

bnet5
–

Table 2.4: CPLEX performance - running time (sec)uniformly distributed points. The resulting objective function was then used to solve
the optimization problem over all time periods for each commodity separately using
CPLEX linear solver. We used the solution obtained by this method to evaluate theHessian for the second order Taylor series approximation of the objective function.
Note that the purpose of these tests was not to attempt to find a good solution by
optimizing a series of second order approximations, but rather to get an understanding of the nature of problem difficulty and thus we optimized only one quadratic
approximation for each test case.
Table 2.4 shows the performance of CPLEX on the quadratic relaxation of the
problem, when run on the problems in Table 2.3. It is evident that the CPLEX
efficiency deteriorates rapidly as the network size and density increase. Thus clearlya state-of-the-art implementation of Newton’s method, by itself, is unlikely to improve
the performance of a standard nonlinear optimization procedure.
Consequently, our work focuses on an algorithm that exploits the specific structureof the problem constraints and the objective function. Our algorithm will essentiallydecompose the problem into a set of smaller problems and will project out a subset

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

42

of variables; thus drastically decreasing the size of the matrices tackled in Newton’s
method.

2.6

Protection models

The core model described in the previous section allows us to design a network thathas enough capacity to accommodate the required amount of traffic between different
node pairs. However this design does not take into account the reliability of theresulting network. Therefore we need to expand this model to accommodate theprotection/restoration techniques outline in Section 2.2.2.
In other words, we need to allocate additional amount of capacity on some alternative paths for every demand, such that in case of a single or multiple link failure there
is enough slack capacity to restore hundred percent of the affected traffic. In what
follows we provide the detailed description of the two previously described approaches
to reliability requirements.

2.6.1

Shared protection on single demand (SD)In this model we assume that for any commodity not more than one of the commodity
paths can fail at a time. In case of a failure the affected portion of the demand has tobe rerouted through one or more new paths that are link-disjoint with all other paths.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

43

In this case the additional capacity allocated on the new paths can be used to restore
flows from any path of the working set. In principle, this condition can be relaxed
to allow the restoration capacity to be also allocated on the remaining (functioning)paths of the working set. However, in this case this capacity can not be shared acrossall the paths of the set and we do not consider this option in our work.
Let Pk be the given fixed set of link-disjoint paths for demand dk and Q ⊂ Pk be theworking subset used to carry the primary flows of dk . Let f (p) be the flow carried onsome path p ∈ Q and let f be the maximum of these flows, that is f = max p∈Q f (p).
Then the amount of the protection capacity required to restore any primary flow
equals f . This condition can also be viewed as sending an additional f units offlow on some paths of Pk \ Q. In order to formulate this model exactly, one would
need to consider all possible subsets of P which can result in a very large amount ofvariables and constraints even for a moderate size network. Instead, we restrict theset of restoration paths to be just a singleton. With this additional constraint themodel formulation becomes more compact and much less complicated. We replace
the demand constraint dtk =Pt hh f (pk ) with a set of constraints

fkt ≥ f t (phk ) ∀phk ∈ Pk

(2.6.1)CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

dtk =

X
h

2.6.2

f t (phk ) − fkt44(2.6.2)

Shared protection across demands (MD)This model is a generalization of the previous formulation to allow sharing of protection capacities across the multiple demands. Consider a network in Figure 2.2
illustrating shared protection on single demand between nodes i and j. Denote theset of paths Pij for this demand as P .

Figure 2.2: Protection example

Let f (1) = 2, f (2) = 4, f (3) = 4, and d = dij = 6. Therefore 2 units ofcapacity are installed on p1 and 4 units on each p2 and p3 . Either p2 or p3 can beused for restoration, however if either of them is used to restore f (p1 ) the protection
capacity will not be used completely. Therefore even if this path fails we can use theexisting protection capacity to restore the flow from some other failure. Moreover,CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

45

assuming that only one link can fail at a time the entire set of restoration capacities
can be utilized by another demand. Consider again the example in Figure 2.2. Let
p5 = v − i − p3 − j − w and suppose dvw = f (p4 ) = 4 with p4 being a primary path.
Then if p3 is the restoration path for dij the existing protection capacity on p3 canalso be used to restore f (p4 ) in case p4 fails.
The only obstacle here is that it is not known in advance which of the paths will
be used for protection. Consequently, in order to be able to share protection acrossthe demands we need to specify the restoration capacity explicitly.
In order to do that we again consider the given link-disjoint set path P and forevery pair paths p, s ∈ P introduce a variable f (p, s) representing the part of f (p)that will be routed on path s in case p fails. In other words, f (p, s) corresponds tothe protection capacity required on s to restore the traffic from p. Note that unlike
the previous model, here we allow the primary flow to be rerouted through more than
one additional path as well as some of the working paths. Thus the assumption that
only one link can fail at a time is crucial for our formulation. With this informationin hand we can now introduce the set of constraints to formulate the model

1. As in the core model again the demand is equal to the sum of the primary flows
dtk =

X
h

f t (phk )(2.6.3)

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

46

2. Every primary flow has to be restored completelyX

f t (phk ) =0f t (phk , phk )(2.6.4)h6=h0

3. The capacity on any link l has to be sufficient to meet all primary flows plus
the protection flows, taking into account capacity sharing across the demands.
Let zlt be the total protection flow on link l at t.
Then the capacity constraint becomesXk,h:l∈phkf t (phk ) + zlt ≤ ylt +

X q,t
yl
q<t

∀l ∈ A, ∀t(2.6.5)

Consider a failure of some link e. For any demand that has e in one of its
working paths the flow routed on that path has to be restored. Consequently,
the protection capacity on l has to be at least the sum of all restoration flowstfrom e on l. Let zel
represent the protection capacity required on l in case e

fails. Then
t
zel
=

X

0

f t (phk , phk )(2.6.6)h00k, h:e∈phk , h : l∈pk

where the summation is taken over all demands which can be affected by the
failure of e and that can use l in one of the restoration paths. In other words,
the summation is taken over all the demands that have disjoint paths using land e.

47

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Therefore, since we assume that only one link can fail at a time and all theprotection capacity can be shared, the total protection capacity required on l isthen the maximum over all possible failurestzlt = max zele(2.6.7)

Our goal is then to maximize the net present value over (2.6.3)-(2.6.7) and (2.4.2).

2.6.3

Flows on the paths

Generally, telecommunication system design has three major components that contribute to maximization of revenue under given capacity and routing constraints.
1. Strategic planning: decide long-term investment in network systems, manage
capital expenditure on strategical deployment of network capacity.
2. Marketing and sales: set price, acquire customers.
3. Network operation: allocation of existing capacity to customers, traffic engineering and routing. This stage is targeting to make the maximum use of existing
capacity. In other words, at this stage the decisions about network usage are
made. In particular one determines the routing of all the demands, that is thefraction of each demand routed on each of its paths.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS48

The design optimization can be made sequentially over each of the stages. In thatcase each part is optimized separately based on the solution of the previous stages.
This approach allows a relatively easy solution at the expense of the quality of theresult. On the other hand, integrating all three parts together creates a much harder
optimization problem, but can provide a better overall solution. However, the design
stages are usually performed by different teams and the integration is not always easy,
if possible at all.
Therefore we would like to compare different approaches to determine whether it
is worthwhile to go through the costly integration.
In this work we consider two options. The first option is to optimize the low-leveltraffic distribution of stage 3 first and use the solution to solve 1 and 2 together. The
second option is to integrate all three parts of the process.
Note that the models defined in the previous section reflect the second approach.
We now describe the changes to that model that take into account the solution ofstage 2. Recall that stage 3 determines how the demands are routed over the network.
That is, for every demand dk in addition to a fixed set of paths Pk we are given a
fixed fraction or percentage of the demand that has to be routed on this path.
For every phk ∈ Pk we are given a percentage constant rkh representing the fractionof dk that has to be routed on phk . For the SD protection modelCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

X
h

49

rkh = 1 + max{rkh }, ∀khas defined in section (2.6.1). Using the percentage constants we then can represent
the flow variables via the demand variables as

f t (phk ) = rkh dtk
and the optimization problem in this case becomesmaxF (d, y)(2.6.8)s.t.
Xk,h:l∈phkrkh dtk ≤ ylt +

X q,t
yl

∀l ∈ A, ∀t(2.6.9)ylq,t ≤ ylq,t−1∀l, ∀q, ∀t ≥ q + 1(2.6.10)

ylq,t ≥ 0∀l, ∀q, ∀t ≥ q(2.6.11)q<t

For the MD model, in addition to the percentage constant rkh for each path phk ∈ Pk
0

we are also given a set of constants rkh,h for each pair of paths phk , phk ∈ Pk , such that

X
h

X

h6=h0

rkh = 1 ∀k
0rkhh = rkh ∀q ∀k0CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS50

as defined in Section 2.6.2. Using this constants we rewrite problem 2.4.1 as

max F (d, y)s.t.
Xk,h:l∈phkrkh dtk + zlt ≤ ylt +X q,tyl ∀l ∈ A, ∀tq<tylq,t ≤ ylq,t−1
X

t
zel
=

(2.6.12)∀l, ∀q, ∀t ≥ q + 1
0rkhh dtk(2.6.13)(2.6.14)h00k, h:e∈phk , h :l∈pktzlt = max zeleylq,t ≥ 0

(2.6.15)(2.6.16)

Fixing the percentages of the demands over the paths reduces the size of the
problem and makes the optimization process faster and easier. In addition, thismodel allows the network operator to specify his own requirements and constraints
on the distribution of the demands as opposed to choosing the best allocation fromthe profit point of view.
In what follows we propose the algorithms and bounds for each of the four resulting
models and evaluate the implications of the integration.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.7

51

Algorithm outline

Recall from the Section (2.4) the two major solution stages:
1. Relax the integrality constraints on the capacity variables and solve the resulting
continuous nonlinear problem
2. Fix the values of the demand and flow variables found in step 1 and find theappropriate values of the capacity variables.
Before we move to the main algorithm in Step 1 let us first describe the idea
behind step 2. Recall from (2.3.3) that the objective is a concave function of thedemand variables and a linear function of the capacity variables. Therefore after thedemand values are fixed the objective becomes linear. Moreover, since the set of paths
is given and fixed for every demand pair, the values of the flow variables provide thetotal load on every link. The problem in step 2 becomesCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

maxX

ht (−Costt ) ≡ min

t

XXX
t

l

ht (ctl (u)ylt (u) +

u

52

X q,t
cl (u)ylq,t (u))
q<t(2.7.1)s.t.
Xu

M (u)(ylt (u) +X q,tyl (u))) ≥ Ltl∀l ∈ A, ∀t(2.7.2)

ylq,t (u) ≤ ylq,t−1 (u)∀l, ∀q, ∀t ≥ q + 1(2.7.3)

ylq,t (u) ∈ N+∀l, ∀q, ∀t ≥ q(2.7.4)q<tWhere Ltl is the load on link l at t computed as
Ltl =

X

f t (phk )k,h:l∈phk

for the SD protection model and as

Ltl =

X

k,h:l∈phk

for the MD model.

f t (phk ) + zlt

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS53

It is easy to see that problem (2.7.1) is a combination of separate optimization
problems for every link, of a form
min

XX
u

ht (ct (u)y t (u) +tXcqt (u)y qt (u))(2.7.5)

<t

s.t.
X
u

M (u)(y t (u) +Xq<ty qt (u)) ≥ Lt∀t

(2.7.6)y qt (u) ≤ y q,t−1 (u)∀q, ∀t ≥ q + 1(2.7.7)y qt (u) ∈ N+∀q, ∀t ≥ q(2.7.8)Note that for t = 1 the problem 2.7.5 is a regular knapsack problem. For t > 1
it becomes a multiperiod knapsack problem and can be readily solved by a standard
integer programming solver. The size of each problem depends on the number oftime periods different batch types and usually does not exceed 300 − 400 variables
and constraints. Note that in addition the constraint matrix is sparse, thus allowing
to obtain an optimal solution in a matter of seconds.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.7.1

54

The nonlinear optimization algorithm

The algorithm in Step 1 is primarily based on the well-known active set method [3].
Consider a general nonlinear problem
min f (x)s.t. Ax ≤ b
1. The active set algorithm starts with a feasible solution which provides a set A
of tight or active constraints.
2. Ignore the rest of the constraints and optimize the problem only over the current
set of active constraints.
3. If the new solution violates any of the remaining constraints, compute the step
length to the constraint boundary and update the set of active constraints. Otherwise, check the Lagrange multipliers and drop the corresponding constraint
from the set.
4. If no constraints can be dropped, stop. Otherwise, goto Step 2.
Under an appropriate nondegeneracy assumption, this algorithm terminates in a
final number of steps ([3], [23]). In order to carry out Step 2 we usually have to invert
the submatrix corresponding to the active constraints. Let BA be a basis of A and
NA be the remaining columns.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

55

Set
xB = (BA )−1 (b − NA xN ) ≡ BxN + dand optimizef¯ ≡ f (x̄)
where x̄ ≡ BxN + d, xN .
In general, it is not known in advance which variables will be projected out;
each time Steps 2 and 3 are performed the subset of the projected variables andtherefore f¯ may change significantly. This property requires a general optimization
algorithm in Step 2 that can turn out to be both inefficient and hard to implement.
In our algorithm we explore the specific structure of the constraints and the objectivefunction with the result that the same set of variables is projected at each step. Thus,
only minor changes to the structure of f¯ occur at any step of the algorithm. Morespecifically, the values of the demand and flow variables are chosen such that the entire
set of capacity variables can be projected out every time Step 2 of the algorithm iscarried out. Further, since the objective function (2.3.3) is a linear function of thecapacity variables, the nonlinear structure of f¯ remains unaffected throughout thealgorithm. In addition we utilize the properties of f¯ to obtain a tight upper boundat no additional computational expense.

56

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.8

The algorithm

In this section we will describe the optimization algorithm in detail for the SD modelwith fixed percentages and then discuss the adjustments and extensions of the algorithm necessary to apply it to the other three models.

2.8.1

Shared protection on single demand with fixed percentages (SDP)

Let d = {dtk } ∈ Rnd and y = {ylq,t } ∈ Rny , where nd = T n(n−1)
and ny = mT (T2 +1)
2Then the constraints can be rewritten as
Ad + By ≤ 0(2.8.1)

where the entries of A are the corresponding values of the percentages for the link
capacity constraints (2.6.9) and 0 for the maintenance (2.6.10) and nonnegativity(2.6.11) constraints and the entries of B are accordingly 1 in (2.6.9) and 0 − 1 in
(2.6.10) and (2.6.11).
Let d0 be some given values of the demand variables. Then the corresponding
optimal values of the capacity variables y 0 can be easily determined by solving a setof linear programs. Let E 0 be the resulting subset of constraints (2.8.1) that aresatisfied as equalities at (d0 , y 0 ) and I 0 be the remaining constraints. Thus E 0 is theset of all tight or active constraints.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

57

Definition 1 A region is the set of all values of d that induce the same subset of
active constraints.
Let BE 0 be the submatrix of B consisting of the tight constraints. Assume that d0 issuch that it produces a unique optimum set of y 0 , that is BE 0 is a square nonsingularmatrix. This can always be achieved by a small perturbation of d0 .
Any such d0 then belongs to the interior of the unique corresponding region. Letd¯0 be some other set of demand values such that there are exactly two constraints
i and j such that i ∈ E 0 , i 6∈ Ē 0 and j ∈ Ē 0 , j 6∈ E 0 . Then d0 and d¯0 induce
two neighboring regions and their common boundary is defined as E 0 ∪ {j}. The
dimension of this boundary is then nd − 1.
More generally
Definition 2 For any region R defined by a set E of tight constraints, and for any
S ⊆ I, E ∪ S is a boundary of R of dimension nd − |S|
Consider the following example. Let G be a single link network, i.e. n = 2, m = 1,
with t = 2. Let 1 = 2 = 1.5, a = 1, c1 = 10, c2 = 9, µ = 0.3
Then the problem formulation is

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

max(d1 )

1/3+ (d2 )

1/3

− 10y 1 − 9y 2 − 3y 1,2s.t.
d1 ≤ y 1d2 ≤ y 1,2 + y 2
y 1,2 ≤ y 1
y 1 , y 2 , y 1,2 ≥ 0
and nd = 2, ny = 3.
Clearly for any values of d1 , d2 such that d1 < d2 the active constraints are
d1 = y 1
d2 = y 1,2 + y 2(E1 )y 1,2 = y 1

and for any values of d1 , d2 such that d1 > d2 the active constraints are

d1 = y 1
d2 = y 1,2 + y 2
y2 = 0(E2 )58

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

59

These two sets of demand values define two regions E1 and E2 . Let S = {y 2 = 0}.
Then E1 ∪ S is the boundary of E1 and E2 . Thus, the set of demands is split into
two regions by the line d1 = d2 . Note that adding a different constraint to E1 or E2
will define the d1 = 0 or d2 = 0 boundaries. See Figure 2.3

Figure 2.3: An example of the region set

In the interior of the region the values of y can be represented as
y = −BE−10 Ad(2.8.3)Thus in the interior of a region the values of the capacity variables are uniquely
defined by the values of the demand variables, therefore in every region we can project
the space of all variables onto the space of the demand variables. The projectedCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

60

objective function then becomes0

f (d) ≡

X
k

ht



gkt (dtk )− C T BE−10 Ad(2.8.4)

where C = {cq,tl }Region Projection (RP) Algorithm Outline
1. Let d0 be a starting point in the interior of some region R0 and f 0 (d) be theprojected objective function. Let i = 0.

2. Region Optimization. Find maxd∈Ri f i (d)(a) Let A = ∅(b) Find d∗ - the maximum of f i (d) subject to the constraints in A(c) If d∗ ∈ Ri stop, maxd∈Ri f i (d) is obtained.

Otherwise determine thesteplength α to the region boundary. Let C be the constraints defining
the boundary. Reset A ← A ∪ C and goto step (b).
3. Using (2.8.3) find the values of the capacity variables y ∗ corresponding to d∗ .

¯ ȳ) be the optimum of g subject to the original
4. Let g = ∇F (d∗ , y ∗ ) and let (d,constraints (2.6.9)-(2.6.11)

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

61

5. Obtain the optimum steplength β and set (di+1 , y i+1 ) = β(d∗ , y ∗ ) + (1 − β)(d̄, ȳ)

6. If the improvement F (d∗ , y ∗ ) − F (di+1 , y i+1 ) is less than F (d∗ , y ∗ ) for sometolerance  stop. Otherwise, goto step 2.

It is easy to draw a parallel between the above algorithm and the general active set
method described in the previous section. In particular, step 2 of the RP algorithm
corresponds to the optimization of f¯ in the general active set method in case thecurrent solution violates the constraint boundary, whereas step 5 of the RP algorithm
is associated with dropping a constraint from the active set. Note however that in ouralgorithm we allow to remove more than one constraint from A, in fact the entire setof active constraints may change after step 5 of the RP algorithm is performed. Thisin turn may affect the convergence rate of the algorithm if an exact solution is desired.
Nevertheless this algorithm proves to be effective when a reasonable optimality gap
is acceptable.
The constant  in the stopping criterion in step 6 is usually taken to be a fractionof percent. We introduce this criterion in order to prevent tailing off when the solution
is close to the optimum.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

62

Region Optimization
We will now describe our implementation of step 2 of the RP algorithm.

1. Using (2.8.3) we project the original set of constraints (2.8.1) onto the space ofthe demand variables. Thus the constraints become

Ap d ≤ 0
where Ap = A − BBE−10 Ad0
2. Note that when A = ∅, f i (d) is a separable function and hence can be easily
optimized. Let d∗ be the optimum solution. Set dˆ ← d0 .
3. Find the steplength α to the region boundary asT

α = min{1,i−api d∗Tˆ
api (d∗ − d)

ˆ < 0}|api T (d∗ − d)4. If α = 1 exit. The optimum value d∗ does not violate the region boundary andthe region maximum is found.
ˆ
5. Set dˆ ← dˆ + α(d∗ − d)6. Let C be the set of new tight constraints. Set A ← A ∪ C. Note that C might
contain more than one constraint in case of a tie in step 3.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

63

7. Optimize f i (d) subject to the constraints in A. For this step we use an optimization method based on the null space representation of A. Recall that theour constraint set is

Ad = 0Let Z be a basis of the nullspace of A. Recall that A ∈ Rnd x k where k is thenumber of the active constraints. Hence Z ∈ Rnd x(nd −k) . Then any vector din the nullspace of A can be represented as a linear combination of the basis
vectors as d = Zv where v ∈ Rnd −k .
We then proceed by optimizing f i (Zv) using standard Newton’s method.
8. Let v ∗ be the optimum solution. Set d∗ ← Zv ∗ and goto 3.
Bounds
We can now describe the process of obtaining a tight upper bound for the problem(2.6.8)
Recall from the description of the RP that every time step 2 with A = ∅ is carriedout, the global optimum of a region function is found. We show that this optimum,
which may actually be located outside the given region, can be used as an upperbound on F (d, y), using the fact that F (d, y) is a continuous concave function and soCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

64

are all f i (d).
Let F (d) be the projection of F (d, y) onto the space of the demand variables.
Then

F (d) = {f i (d)|d ∈ Ri }is a continuous and concave function of d. Recall from (2.8.4) that all region functions
f i (d) have the same nonlinear part. Thus, F (d) can be defined as the sum of a smoothconcave nonlinear function and a piece-wise linear function. Since F (d) is also smoothand concave, the piece-wise linear function is concave as well. This also follows fromthe fact that the values of the capacity variables in every region are determined as
an optimal solution to a minimum cost LP, which takes the demand values as input
right-hand-side parameters. It is well known that in this case the objective functionof the LP is a convex function of these parameters. Recall that in the original profit
function and therefore in all the projected functions the cost of the capacities is
subtracted from the revenue. Thus, the piece-wise linear function resulting from thecapacity cost is concave.
Therefore, for any region its region function is dominated by the function of any
neighboring region, and since the functions of any two regions intersect only along
the regions’ common boundary we have proved

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

65

Theorem 1 For any two regions Ri and Rj
f j (d) ≤ f i (d) ∀d ∈ Rj
Therefore the global optimum of any region function obtained at step 2 of the RP
algorithm provides a valid upper bound on the optimum value of F (d, y). Note that
this bound is obtained at no additional computational expense and we will show in
Section 2.10 that this is in fact a tight bound.

2.8.2

Shared protection across demands with fixed percentages (MDP)

In this model we have a set of supplementary variables {ztl } representing the maximum
protection load on link l and a set of constraints linking these variables with thedemands. Clearly the values of {ztl } are uniquely defined by the values of {dtk },moreover these auxiliary variables do not affect the objective value and used only todetermine the set of demands that defines the protection capacity on the given link.
Thus we can again project the space of all variables onto the space of the demandvariables and use the algorithm described in the previous section.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.8.3

66

Shared protection on single demand without percentages (SDN)

In this model in addition to the demand and capacity variables we also have flow variables and the constraints reflecting the flow-demand dependence. As in the previous
model the values of the capacity variables are uniquely defined by the flow on thelink. Furthermore, we have an additional set of constraints linking demand and flow
variables. Note that the values of the flow variables uniquely define the values of thedemand variables. Thus we redefine the region concept to capture the flow-demand
relationship. Let f = {f t (phk )} ∈ Rnf where nf = ndP

k |Pk |

Definition 3 A region is the set of all values of f that induce the same subset of
active constraints.
The nonlinear optimization algorithm then proceed in the same manner by projecting out the capacity and demand variables and sequentially optimizing the regions
defined by the flow variables.

2.8.4

Shared protection across demands without percentages(MDN)

For this model we utilize the ideas described in the previous subsection. Using(2.6.3)-(2.6.4) we can immediately project out {dtk } and {f t (phk )} variables. Let

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS0f = {f t (phk , phk )} ∈ Rnfp where nfp = nd

P |Pk |(|Pk |−1)
k

2

67

and redefine regions using

Definition 3 The optimization algorithm then proceeds in the same manner as for theprevious model.

2.8.5

BoundsRecall that the bound for the SDP model is based on the properties of the original
and projected objective functions, specifically the fact that all these functions areconcave and any two projected functions intersect along a single hyperplane. It iseasy to see that these properties are maintained in the other three models as well andhence the same approach to obtaining a valid bound can be utilized. fi

2.9Starting point

Recall from section (2.2.5) that the actual convergence rate of the Newton method
and therefore that of the RP algorithm relies on the quality of the starting point.
On order to obtain a good initial point we start with a first order heuristic. The
central idea of the heuristic is optimize each commodity separately, thus replacing
the large problem with O(n2 ) smaller problems. Clearly, for T = 1, that is fora single time period, optimizing one commodity at a time will provide an overall
optimal solution. On the other hand, in the multi-period problem we have a certain

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

68

amount of capacity maintained from the previous period, that has to be shared bydifferent commodities. One way to resolve the sharing problem would be to split theexisting capacity in the same proportion it was originally used. However, recall thatthe parameters of nonlinear terms of the objective function differ not only between
different commodities, but also between different periods of the same commodity.
Thus in general, this approach will not provide an optimal solution, nevertheless itturns out to be a reasonable starting point in our experiments.

2.9.1

Fixed percentages

We again start with the description of the algorithm for the simplest model, namely
the SD model with fixed percentages. The single commodity problem for this instance
is

max

X

ht

t



gkt (dtk ) −Xl

ctl ylt +X q,t q,t 
cl y lq<ts.t.
Xq:l∈pqr q dtk ≤ ylt +

X q,t
yl
q<t

∀l ∈ A, ∀tylq,t ≤ ylq,t−1∀l, ∀q, ∀t ≥ q + 1

ylq,t ≥ 0∀l, ∀q, ∀t ≥ qConsider the example in Figure 2.4. Let d be the commodity i, j with three given

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

69

Figure 2.4: Original network for SDP modelpaths and the requirement of 50% of d over each path. Again denote the set of paths
Pij as P . Note that the capacity on every path at every time period is uniquely
defined by the value of d at that period. Moreover, the capacity is the same on every
link of a path. Then our problem is equivalent to solving the same problem over thecompressed network Ḡ in Figure 2.5 with the cost of the accumulated linkc̄ t =

X

rh

r̄ =

X

h

X

l : l∈phctl



(2.9.1)and the percentage requirementrhhFigure 2.5: Compressed network Ḡ for SDP model(2.9.2)

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS70

We can now simplify the single commodity instance as

gkt (dtk ) − c̄ t y t +

X

X

∀t

(2.9.3)y q t ≤ y q,t−1∀q, ∀t ≥ q + 1(2.9.4)y qt ≥ 0∀q, ∀t ≥ q(2.9.5)

max

X

ht

t

q<t

qt

c̄ y

qt

s.t.
r̄dtk ≤ M (y t +y q t)

q<tThis is a much smaller problem with only O(T 2 ) variables and can be easily solved
by any of the cutting plane methods described in section (2.2.4).
For the MD model we utilize the same approach of reducing the set of paths to asingle link. However the cost model in this instance is more complicated due to thestructure of the capacity usage.
In order to formulate this model we introduce an additional constant.
Consider the example in Figure 2.6
We are given three commodities
1. d34 with 100% on primary path c and 100% on protection path d − e
2. d12 with 100% on primary path a and 100% on protection path b − e
3. d13 with 100% on primary path a − b and 100% on protection path e

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

71

Figure 2.6: Original graph for the MD model
The protection capacity on link e can be shared by demands d34 and d12 or d34
and d13 , but not d13 and d12 since they share link a on their primary paths. Thus
the protection capacity on link e will be defined by the maximum of d34 and d12 + d13
since in case link a fails both d12 and d13 will be rerouted through link e. In otherwords we can divide the set of the commodities into two groups {d12 , d13 } and {d34 }that share the capacity on link e.
Denote gl the number of groups on link l. In our example gl = 2, gb = gd = 1.
The number of groups can be easily precomputed. Using this constant we then define
the cost c̄t as

t

c̄ =

X X
l

h : l∈phr h ctl +P

th0 : l∈ph0 rhh0 cl 

gl(2.9.6)

Thus the cost of the protection capacity is equally split between the groups that

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

72

share the protection capacity on this link.

2.9.2

No percentages

If the percentages are not given, then the proportion of the demand on each path is
not known, however the load on each link of any path is still the same and the network
in Figure 2.6 becomes a compressed network with 2 nodes (the demand nodes) andHk links.

Figure 2.7: Compressed network Ḡ for SDN model

Then the SD model can be formulated as

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

X

max

ht

t



gkt (dtk ) −HkXc̄ht yht +

X
q<t

h=1




c̄hq t yhq t

73(2.9.7)s.t.
f t (ph ) ≤ M (yht +

X qt
yh ))

1 ≤ h ≤ Hk , ∀tq<tyhqt ≤ yhq,t−1

1 ≤ h ≤ Hk , ∀q, ∀t ≥ q + 1

f t (ph ) ≤ f t
X

dtk =

1 ≤ h ≤ Hk

ph ∈P (d)f t (ph ) − f t ∀tyhqt ≥ 0∀1 ≤ h ≤ Hk , ∀q, ∀t ≥ qwherec̄ht =

Xctl

l : l∈phThis problem is slightly larger than the instance with percentages, neverthelesssince the number of paths is usually small (< 10) it remains an easy problem for the
first-order method.
Note that for the single commodity case there is no essential difference between
the SD and MD models. Thus we can use model (2.9.7) with costsc̄ht =Xl : l∈phr h ctl +P

h0 : l∈ph0 r

gl

hh0 t
cl 

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

74

The solution to model (2.9.7) provides the values of the demand and flows on eachpath for the given commodity. Recall that for the MD model we in addition requireto specify the values of the protection flows explicitly, that is for each path pk we00need the values of the protection flows f t (ph , ph ) for every path ph , h 6= h0 and thevalues of the primary flows f t (ph ) for every path ph such that0f t (ph ) + max
f t (ph , ph ) = f t∗ (ph ) ∀t, 1 ≤ h ≤ Hk0
h 6=h(2.9.8)where f t∗ (ph ) is the optimal solution to (2.9.7)
0Note, however, that any values of f t (ph , ph ) and f t (ph ) that satisfy (2.9.8) will
provide a solution of the same cost. More precisely, let f t∗ = maxh f t∗ (ph ). Clearly
there are always at least two paths that carry the flow of value f t∗ . Let ph(t) be oneof these paths and let

f t (ph ) = f t∗ (ph∀h) 6= h(t)f t (ph(t) ) = 0(20f t (ph , ph ) = 0∀h0 6= h, h(t)f t (ph , ph(t) ) = f t (ph )0Obviously these values of f t (ph , ph ) and f t (ph ) satisfy (2.9.8).

(2(2

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

2.10

75

Implementation and numerical results

In order to evaluate our algorithm and the bound we generated 5 data sets based on
different real-life networks. In addition a set of larger networks was created to test theboundaries of our method. Before we proceed with the descriptions of the data sets
and the computational results let us first give some details on the implementation ofsome parts of the algorithm. The Newton method in step 7 of the RP algorithm was
implemented using the CLAPACK [15] library for sparce matrix operations including
Cholesky factorization. Recall from Section 2.7 that the second part of the algorithm
involves solving a series of a multiperiod knapsack problems. We use CPLEX integer
solver to solve these problems to optimality within a fraction of time required by
the nonlinear optimization algorithm. The resulting integer solution was in all cases
within a fraction of a percent of the nonlinear solution due to the reasons describedin Section 2.4. For this reason we do not provide these results here.
The real-world problems were solved on a 336 Mhz UltraSPARC machine with
1.7Gb of RAM and the test with larger networks were conducted on 1.89Ghz Xeon
with 3Gb of RAM.

2.10.1

Real-world problems

Table 2.5 shows the statistics of the first set of data. Here the last two columns show
the minimum and maximum number of paths per node pair. The small number ofCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

76

paths used is due to the sparsity of the real-life networks. The paths were generatedusing the SPIDER code[14]. The time horizon T was taken to be 14 in all test sets.
Note that this set is identical to the set used in Section 2.5
Recall from Section 2.3 that there are six main data parameters defining the
problem data for every model, whose values we describe next.
The maintenance fraction µ and the maintenance rate α were taken to be equal
to 0.05 and 1.01 correspondingly.
The number of batch types u, the size of each batch M (u) and the initial cost
per unit of capacity are reflected in Table 2.7. The discount factor ht was taken to
be 0.86 for every time period t ≥ 2. The parameters were chosen to approximate thereal data patterns and behavior as closely as possible.
The initial buying costs c1l (u) were also picked to reflect the real cost structure of
the networks based on the installation and mileage costs of the equipment. The buying
costs for all the subsequent time periods t = 2 . . . T were computed as ctl (u) = γct−1l (u)for a given constant γ < 1.
Table 2.6 shows the range of the remaining parameters.
All possible combinations of these parameters were considered (a total of 27) and
50 random sets were generated for each combination. The data corresponding to eachparameter was independently and uniformly distributed with mean shown in table
2.6. Here we take A =P

k Ak .

In all runs the code was terminated when the total

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Set#Nodes#Links#Paths (min)#Paths(max)

snet1

14

22

2

3

snet2

38

48

2

3

snet3

47

55

2

4

snet4

50

64

2

3

snet5

70

94

2

3

77

Table 2.5: Statistics on real-life networks
Ā

γ̄

¯

50000

0.95 1.3

5000000.9

1.45000000 0.85 1.5
Table 2.6: Data parameters
improvement in the last 10-20 steps did not exceed 0.05%.
In what follows we present the detailed numerical analysis of the first algorithm
and the summary of the results for the other three models. In our experiments the
behavior of all models follows the same patterns with only slight variations due tothe different size of the models.
Table 2.8 illustrates the average and boundary values of the performance of theCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

u

1

2

3

4

5

M (u)

198

198

198

762

762

c1 (u)M (u)0.65 0.625 0.42 0.17 0.16

6

78

7

1536 1536
0.1

0.08

Table 2.7: Capacity costs
algorithm on different parameter groups. Table 2.9 shows the comparison of therunning time. In order to single out the pattern for each group of parameters, thedata was averaged on the other two groups.
Table 2.8 clearly shows that the problem becomes easier to solve as the elasticity
 and the yearly reduction rate γ parameters decrease. On the other hand the scaling
constant A does not affect the performance. This observation supports the natural
assumption that the problem becomes harder to solve as the demand-price curve thatis essentially defined by the elasticity becomes steeper. In addition, smaller valuesof γ induce large decrease in cost from one year to another and helps to reduce theproblem degeneracy, thus allowing a faster convergence.
Table 2.10 compares the average optimality gap and the running time for the solution obtained by the starting heuristic with the final solution. The data is averaged
over all parameter groups. This table illustrates the quality of the heuristic which
in most cases provides a good starting solution. The running time of the heuristic isapproximately half of the overall running time. Thus for a case where only a rough79

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

net

snet1

snet2

snet3

snet4

snet5

Gap

Ā

γ̄

¯

(%)

50000

5000005000000

0.9

0.8

0.7

1.3

1.4

1.5

Ave

0.000810.000770.000820.001

0.000780.000620.00098

0.000750.00065Max

0.00850.00820.00880.00980.00810.00680.00960.00820.0071

Min

0

0

0

0

0

0

0

0

0

Ave

0.006

0.00610.0061

0.008

0.006

0.00410.00860.0062

0.004

Max

0.051

0.052

0.051

0.063

0.051

0.039

0.071

0.053

0.04

Min

0

0

0

0

0

0

0

0

0

Ave

0.019

0.018

0.016

0.044

0.017

0.009

0.039

0.018

0.01

Max

0.078

0.076

0.075

0.088

0.076

0.064

0.089

0.075

0.065

Min

0

0

0

0

0

0

0

0

0

Ave

0.019

0.023

0.021

0.031

0.021

0.011

0.033

0.021

0.012

Max

0.11

0.125

0.14

0.25

0.13

0.1

0.21

0.12

0.096

Min

0

0

0

0

0

0

0

0

0

Ave

0.021

0.024

0.025

0.041

0.023

0.012

0.033

0.022

0.013

Max

0.193

0.198

0.21

0.29

0.2

0.176

0.28

0.196

0.18

Min

0.011

0.014

0.009

0.018

0.013

0.007

0.019

0.014

0.007

Table 2.8: Performance comparison - SDP model

80

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

net

Ā

γ̄

¯

50000

5000005000000

0.95

0.9

0.85

1.3

1.4

1.5

snet1

15

15

14

16

15

15

16

16

14

snet2

45

47

45

47

46

44

46

46

44

snet3

100

100

99

102

99

99

100

100

98

snet4

150

152

151

153

151

148

153

152

148

snet5

213

212

214

215

213

213

214

214

212

Table 2.9: Average running time (sec) - SDP model

snet1

snet2

snet3

snet4

snet5

Heur. gap (%)

0.00170.016 0.041

0.049 0.056

Final gap (%)

0.00079 0.006 0.018

0.021 0.023

Heur. time (sec)

7

26

47

77

103

Final time (sec)

15

46

100

152

213

Table 2.10: Heuristic performance - SDP modelCHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Model

SDN

MDP

MDN

Gap(%)

snet1

snet2

snet3

snet4

snet5

Ave

0.32

0.35

0.48

0.47

0.65

Min

0.22

0.17

0.2

0.21

0.4

Max

0.44

0.61

0.88

0.85

1.13

Ave

0.27

0.29

0.33

0.33

0.41

Min

0.16

0.18

0.16

0.18

0.2

Max

0.39

0.43

0.52

0.59

0.64

Ave

1.19

1.77

2.01

2.2

2.3

Min

0.82

1.16

1.29

1.56

1.45

Max

1.56

2.47

2.89

3.08

3.89

81

Table 2.11: Performance summary - SDN, MDP, MDN models
estimate of the solution is required, it can be obtained much faster and improved
later if needed.
The summary performance of the other three models on the same data sets isreflected in tables 2.11 and 2.12. Table 2.11 shows the average and boundary optimality gaps for SD model without percentages (SDN) and MD model with (MDP)
and without (MDN) percentages. Table 2.12 shows the average running time for themodels.
Generally the problems become harder to solve as the problem size grows. Nev-

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Model

snet1

snet2

snet3

snet4

snet5

SDN

10

79

119

108

284

MDP

11

80

126

130

304

MDN

13

93

142

158

332

Table 2.12: Running time (sec) summary - SDN, MDP, MDN models

Model

snet1

snet2

snet3

snet4

snet5

SD

1.5

1.7

1.6

1.8

1.85

MD

1.6

1.9

1.86

1.91

1.9

Table 2.13: The effect of fixing percentages

Model

snet1

snet2

snet3

snet4

snet5

SDP/MDP

1.92

1.79

1.98

2.01

2.02

SDN/MDN

2.4

2.34

2.22

2.3

2.19

Table 2.14: Model comparison (%)

82

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

83

ertheless it is evident that the algorithm can handle large instances of the models
without substantial loss of the solution quality. The additional numerical results
supporting this claim will be shown in the following subsection. However before wemove to the more detailed exploration of the scalability of our algorithm let us firstcompare the performance of the different models from the business solution perspective. Recall that the goal of the optimization algorithm is to create a reliable andsurvivable network. SD and MD models were designed to provide a different degreeof survivability, with higher level available at the expense of receiving a lower overallprofit. In addition a different solution technique, precomputing the demand routing
was proposed, allowing to simplify the model and the algorithm once more at theexpense of the profit. Table 2.13 shows the decrease in the objective value when fixed
percentages are used and Table 2.14 demonstrates the average increase in total profit
when the MD model is used instead of the SD model. The average improvement isabout 2% of the entire business profit for all cases considered. Thus a weaker survivability model offers a considerable increase in company profit consequently providing
the grounds for a serious analysis of the survivability and profit tradeoff. On theother hand the decrease in the profit from using fixed percentages is also about 2%
on average. However it is evident from tables 2.11 and 2.12 that reducing the size ofthe problem by using fixed percentages does not provide a considerable decrease in
the running time or noticeable improvement of the algorithm performance. This in

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Solver

snet1

snet2

snet3

snet4

snet5

LOQO

0

0

0

0

–

SNOPT

0

–

–

–

–

Our Code

0.0008 0.006

0.018 0.021

84

0.023

Table 2.15: Performance comparison (LOQO and SNOPT) - Average optimality gap (%)
Solver

snet1

snet2

snet3

snet4

snet5

LOQO

30

600

720

900

–

SNOPT

90

–

–

–

–

Our Code

15

46

100

152

213

Table 2.16: Performance comparison (LOQO and SNOPT)- running time (sec)turn implies that it is most likely not worthwhile to utilize this technique.
In order to evaluate our algorithm we compared it against LOQO [39] and SNOPT[13] Table 2.15 and 2.16 show the performance comparison of our algorithm andsolvers on the real-life networks. Note that some of these results were also presented
and analysed in Section 2.5.
It is evident that even though the comparison is biased toward our code as it istuned up specifically to solve this problem, the existing solvers are still not able tohandle even a reasonable size network and in cases where LOQO was able to solve

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Set#Nodes#Links#Paths (min)#Paths(max)

bnet1

100

500

3

6

bnet2

150

900

3

7

bnet3

200

1500

3

9

bnet4

250

2000

3

10

bnet5

300

4500

4

12

85

Table 2.17: Statistics on the generated networks
the problem, our algorithm significantly outperforms it in terms of running time with
only a slight increase of the optimality gap.

2.10.2

Generated problems

In this section we analyze the performance of our algorithm on the generated set oflarge dense networks. Table 2.17 shows the statistics on the five networks generated
to test the scalability of the algorithm.
Since these cases were generated to test the quality of the algorithm on the large
data set only one set of parameters was chosen for trials. Namely, Ā = 50000, ¯ = 1.5,γ̄ = 0.95. Recall that these values of the parameters proved to be the hardest forthe optimization algorithm. For the same reason the performance was tested on thelargest model, i.e. MD without percentages.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Solver

bnet1

bnet2

bnet3

bnet4

CPLEX

2341

7982

14324 18943

Our Code

1943

2457

2691

2893

86

bnet5
–
3103

Table 2.18: Performance comparison (CPLEX) - running time (sec)
The remaining parameters were taken the same as in the previous section. Similarly to the previous section, 50 random data sets were generated for each network.
We compared our algorithm with CPLEX [16] which is on of the most efficientlinear, integer and quadratic solver available on the market. CPLEX also has oneof the best Cholesky factorization algorithms built in it. In order to compare ouralgorithm with CPLEX we replace the original objective with its second order Taylor
series approximation. The idea here is to ensure that our implementation of the
Newton method is not the bottleneck of the performance, and that the active set
method considerably reduces the complexity of the problem.
Table 2.18 shows the comparison of our algorithm with the performance of CPLEX
on the quadratic relaxation of the problem. It is evident that the CPLEX efficiency
deteriorates rapidly as the network size and density increase, while our algorithm
shows stable and scalable performance.
Finally, tables 2.19 and 2.20 illustrate the optimality gap and the running timeon all the networks.

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS

Gap (%)

bnet1

bnet2

bnet3

bnet4

bnet5

Ave

2.91

2.99

3.15

3.26

3.51

Max

4.56

4.87

4.91

5.1

5.08

Min

1.71

1.82

1.83

1.93

1.96

87

Table 2.19: Performance summary - Optimality gap (%)
bnet1

bnet2

bnet3

bnet4

bnet5

1943

2457

2691

2893

3103

Table 2.20: Average running time (sec) summary

2.11

Discussion

The previous section demonstrates several important conclusions which we will summarize in this section.
• Algorithmic inferences
The numerical results indicate that the approach proposed in this work leads
to a viable algorithm that provides a high quality solution within a reasonable
timeframe and scales well with the problem size. Note that the specific formof the objective function is not essential for the course of the algorithm. The
idea of projecting out the capacity variables and working only with the demandor flow variables can thus be extended to different pricing and cost models

CHAPTER 2. BUSINESS PLANNING OF OPTICAL NETWORKS88

as long as the resulting objective can be formulated as a separable concave
function. Further, this approach might prove efficient for other network design
problems with liner or concave objective functions. The essential detail of theapproach is the unique dependency of the link capacity and flows (demands),
thus many network design problem that poses this characteristic might benefit
from a similar algorithm.
• Network Planning inferences
From the design point of view we confirmed our initial assumption that a higher
level of protection will result in a lower profit value. The shared protection on
the other hand leads to a more complicated model and as a result requires
longer running time to achieve the similar solution quality. In addition we havedemonstrated that the integration the routing problem into the optimization
model can result in a higher profits as opposed to sequential solution. Thus
even though the integration is not currently employed in existent schemes due
to its complexity and difficulty, it is definitely worthwhile to explore this option
as it is evident that an adequate algorithm can be developed.

Overall we conclude that our results provide justification for further, more extensive study of different protection and integration approaches as well as the expansionof our technique to more general planning, design and pricing problems.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

89

Chapter 3
Maximum Concurrent Flow
Problem

3.1

Introduction

The maximum concurrent flow problem frequently arises in practical applications andhas received much attention due to its difficulty.
The input to the problem is a graph with capacities on its edges and a set of multicommodity demands to be routed. The objective is to find a largest value Γ∗ > 0 suchthat a fraction Γ∗ of every demand can be simultaneously routed without exceeding
the available capacities (possibly Γ∗ > 1). In other words, we want to find a feasible flow with maximum throughput. A precise definition is given in Section 3.2. An

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

90

equivalent problem is the minimum congestion problem: here we must simultaneously
route 100% of every demand so that the maximum load of any edge – the ratio of
total flow to capacity – is minimized. In the literature both problems are sometimes
regarded as one and the same and are jointly referred to as the “maximum concurrent
flow problem”, but in this paper we will work with the maximum throughput version.
Despite large advances in computer technology and optimization software, evenstate-of-the-art linear optimizers are often proven ineffective on very large problems.
Furthermore, the scalability of these codes comes into question. As an example, aconcurrent flow problem on 400 nodes, 1740 arcs and 5000 commodity pairs whichyields a model with 160143 rows, 689056 columns and 2067177 nonzeros is solved
by CPLEX in 9434 sec. on a modern computer, whereas the same problem on 700
nodes, 3120 arcs and 5000 comodity pairs yields a model with 481326 rows, 2131038
columns and 6393176 nonzeros. This larger model already requires 455541 seconds(see [4]). With more extensive experiments it can be shown that the running timeusing CPLEX (dual) grows cubically with the number of columns (on concurrent flow
models). Moreover, during a very large fraction of the running time the optimality
error is rather large (more than one percent, say).
This situation has prompted a long line of research on alternative algorithms
for the maximum concurrent flow problem, in particular, approximation algorithms.
Given the nature of the problem, in a sense an ideal approximation algorithm would

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

91

be one whose precision can be controlled. Given  > 0, a fully polynomial polynomialtime approximation scheme is one that estimates Γ∗ within relative error , and whose
worst-case complexity grows polynomially in −1 and on the size of the graph.
The first such algorithm was given by Shahrokhi and Matula [60] in the late 80’s,
and it applied to the minimum congestion problem. [60] introduced a “potential”
function that exponentially penalizes the load on any edge, so that the approximate
minimization of the potential function yields a flow whose maximum load is provably
near-optimal. To carry out the approximate optimization of the potential function[60] used what amounts to a Frank-Wolfe scheme (see [44], [46], [47]). The overall
algorithm was proved to have a worst-case complexity bound that grows proportionalto −7 (as well as polynomially on the size of the graph). The Sharokhi-Matulaalgorithm can also be viewed as a Lagrangian relaxation approach, where the capacity
constraints are relaxed, and the violation of these constraints is penalized using very
special multipliers. These are exponential multipliers, which, roughly stated, arise in
the Frank-Wolfe scheme as the gradient of the potential function.
The Sharokhi-Matula results catalyzed an intense and still active research effort.
[54] has improved the running time of Sharokhi-Matula algorithm. The work in [54]was later extended ([56]) and improved ([59]). See also [50], [51], [57], [58], [53].
This research primarily showed how to better use an exponential potential function,
still in the Frank-Wolfe Lagrangian relaxation framework, to obtain faster algorithms

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

92

for the minimum congestion problem, and generalizations. A logarithmic “sliding
barrier” function, applied to the minimum congestion problem, is discussed in [61].
The analysis in [62] shows how the exponential potential function naturally arisesin the derandomization of probabilistic algorithms; the algorithms for the maximum
concurrent flow problem described in [52] and [45] are related to the ideas in [62].
The fastest of these algorithms have running time bounded by −2 times a loworder polynomial on the size of the graph and number of commodities. Altogether, thealgorithms in [59], [52] and [45] yield the best bounds for the maximum concurrent
flow problem, but precisely which is “best” depends in a complex manner on thenumber of vertices, edges and commodities. Experimental testing of these algorithms
has yielded implementations that are substantially more effective than commercial
linear programming codes, see [41], [49].
The Sharokhi-Matula algorithm was for a long time considered the first −approximation
algorithm to exploit the idea of a potential function that penalizes the edge load. However, in 1971, Fratta, Gerla and Kleinrock [48] proposed a general algorithmic scheme
towards the solution of various optimization problems arising in telecommunications,among them that of finding a feasible multicommodity flow, which is tantamount tosolving a maximum concurrent flow problem. The approach in [48], which is quite
different from those found in the references cited above, involves a simple idea forincreasing the throughput of a feasible flow while using a rational barrier function

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

93

to prevent flows from exceeding capacities, together with a Frank-Wolfe procedurethat reduces the barrier function. A partial convergence proof was presented in [48].
Although the flow deviation method is well-known and frequently used in the telecommunications community, it is poorly known in the algorithms community.
In this work we show that the flow deviation method, properly implemented,
yields an algorithm that solves the maximum concurrent flow problem to relative
error  by solving O(−2 m3 k 2 + m3 k 2 log m) minimum-cost flow computations, wherem and k are, respectively, the number of edges and commodities. The algorithm we
describe contains the critical algorithmic ideas in [48] together with some refinements
designed to achieve the −2 performance. Without these refinements, the algorithmin [48] can be shown to converge in O(ρ2 m4 −2 k 2 ) shortest path computations, whereρ is the width, which is defined as follows. Given a convex set P ∈ Rn and theset of m inequalities Ax ≤ b, the width of P relative to Ax ≤ b is defined by ρ =maxi maxx∈P ai x/bi [58].

3.2

Definitions

Consider a graph G with n vertices, and a capacity ue > 0 for each edge e. Suppose
we are given a set of k commodities, where for 1 ≤ j ≤ k, commodity j consists ofa pair sj , tj of vertices and a demand amount dj > 0. For 1 ≤ j ≤ k, let Pj denotethe set of paths between sj and tj , and for any edge e let Pe,j denote the subset of

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

94

Pj consisting of those paths that contain e. The maximum concurrent flow problem
is the linear programΓ∗ = maxγ(TF) s.t.
kXX

j=1 p∈Pe,jX

xp ≤ ue ∀ e,
xp = γdj ,p∈Pj

0 ≤ xp(3.2.1)

1 ≤ j ≤ k,

(3.2.2)∀p.

(3.2.3)This definition and formulation apply whether the graph is directed or not. Fromthe point of view of linear programming, this path-based formulation is not the mostefficient – a flow formulation is more compact. Note that for any commodity j alinear program over the constraints (3.2.2), (3.2.3) reduces to a shortest path problem.
Given  > 0, an -approximate solution to TF is a flow x̂ with throughput γ̂, feasible
for TF, such that γ̂ ≥ (1 − )Γ∗ , i.e. γ̂ approximates Γ∗ within relative error at most.
We will say that a flow vector x is f easible if there exists a nonnegative value
γ = γ(x) such that (x, γ) is satisfies (3.2.1 - 3.2.3), in which case we will also say thatx has throughput γ. Given a flow vector x, we will say x has load.
λ(x) = max λe (x)e(3.2.4)

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

95

where for an edge e,.
λe (x) =

Pk

j=1

P

p∈Pe,j xp

ue(3.2.5)is called the load of x on e. Thus x is feasible iff λ(x) ≤ 1, and we will say x is
strictly feasible if λ(x) < 1.
Consider a problem of the formmin {F (x) : x ∈ Q}

(3.2.6)

where Q ⊆ Rn is convex and F () is convex and differentiable over Q. The classicalFrank-Wolfe procedure solves this problem by generating a sequence {xt }, t = 1, 2, · · ·contained in Q. At iteration t, the procedure solves the linear program
min [∇F (xt )]T v : v ∈ Q(3.2.7)

with solution y t ∈ Q, and sets xt+1 = (1 − σ t )xt + σ t y t = xt + σ t (y t − xt ) where0 ≤ σ t ≤ 1 is chosen so as to minimize F ((1 − σ t )xt + σ t y t ). σ t is called the stepsize.
The Frank-Wolfe method is closely related to Danzig-Wolfe decomposition, and to
steepest-descent methods for unconstrained optimization.
For completeness, we state the minimum congestion problem. Given a graph andmulticommodity demands as in the maximum concurrent flow problem, the minimum
congestion problem is the linear programΛ∗ = minλCHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

96

s.t.
k
X
X

j=1 p∈Pe,j

xp − u e λ ≤ 0
X

p∈Pj

∀ e,

xp = dj 1 ≤ j ≤ k,

0 ≤ xp∀p.

(3.2.8)(3.2.9)(3.2.10)

A simple analysis shows that Λ∗ = 1/Γ∗ . In what follows, we will deal exclusively
with the maximum concurrent flow problem.

3.3

The central idea

In this section we describe the main idea in [48].
Consider an instance of the maximum concurrent flow problem, and let x̂ bestrictly feasible with throughput γ̂. To fix ideas, suppose that λ(x̂) = 1/2. Then 2x̂
1x̂ is feasible and has throughput
is feasible and has throughput 2γ̂. In general, y = λ(x̂)
1
γ̂
λ(x̂)

> γ̂ (but also λ(y) = 1). Suppose that by using an appropriate algorithmic

1
construct we obtain, from y, a feasible vector z, still with throughput λ(x̂)γ̂ but withloads that are singnificantly lower than those of y. Then we would have a skeletal
outline for an algorithm:

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

97

OUTLINE
A1. Let x̂ be strictly feasible, with throughput γ̂.
1
B1. Let y = λ(x̂)
x̂.
1C1. Find a feasible z with throughput λ(x̂)γ̂ and λ(z) substantially smaller than λ(y).
1Reset x̂ ← z and γ̂ ← λ(x̂)γ̂, and go to B1.

Thus, each execution of B1 increases throughput while C1 “spreads out” the flow.
Next we describe a concrete implementation of C1.
Let ψ : [0, 1) → [0, +∞) satisfy the properties(i) In the range 0 ≤ t < 1, ψ(t) is increasing, continuous and convex,(ii) ψ(t) → +∞ as t → 1− .
Then one way to accomplish C1 is to define, for any strictly feasible x,
Ψ(x) =X

ψ(λe (x))(3.3.1)e

and to choose z as the solution toγ̂min Ψ(x) : x a feasible flow with γ(x) =
λ(x̂)



(3.3.2)

In fact, by property (ii) of ψ the optimal solution to the optimization problem (3.3.2)“should” have maximum load smaller than λ(x̂). Note, however, that in view of ourCHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

98

construction of y in step B1, Ψ(y) is undefined (because λ(y) = 1), and in general,1
it may be the case that Ψ(x) is undefined for any x ∈ λ(x̂)γ̂P . Thus, we must amend
1
Step B1 slightly: to obtain y, we scale up x̂ by a factor slightly smaller than λ(x̂).

Further, it may not be strictly necessary to minimize Ψ – rather, we should decrease
it from the value Ψ(y). Our algorithmic scheme is now:
Basic Algorithm
A2. Initialization. Let x̂ be strictly feasible with throughput γ̂.
1
is a parameter.
B2. Magnification. Set y = µx̂, where 1 < µ < λ(x̂)

C2. Barrier reduction. Find a feasible z with throughput µγ̂ such that Ψ(z) issufficiently smaller than Ψ(y). If no such z exists, STOP. Otherwise, reset x̂ ← z andγ̂ ← µγ̂, and go to B2.
The algorithm implicit in Steps A2 - C2 is the main contribution in [48]; found in
pp. 112(bott.)–114. It significantly differs from algorithms found in the references
cited above. Other ingredients in [48] are:x
1. The barrier function ψ(x) = 1−x,

2. A Frank-Wolfe procedure to carry out the barrier reduction step C2,
3. A choice for the parameter µ in B2 that provides a sufficiently rapid increasein throughput, andCHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

99

4. A provably good choice for the starting point x̂ in A2.
The particular implementation we describe below differs from that in [48] in two
aspects:• Our Frank-Wolfe iterations are minimum-cost flow problems, as opposed toshortest path computations, and• In C2 we may carry out several Frank-Wolfe iterations (i.e. several gradientsteps) whereas [48] uses a single iteration.
If A2 - C2 is implemented using shortest path computations only and single FrankWolfe iterations in C2, the resulting algorithm can be shown to converge to an approximate solution in O(ρ2 m5 −3 k) shortest paths computations, where ρ is thewidth parameter [58]. Finally, we note that in [48] the term “flow deviation method”
appears to refer to the Frank-Wolfe procedure itself; however it has since become
synonymous with the entire scheme A2 - C2.

3.4

The algorithm

In this section we describe our implementation of A2 - C2. We will use the same
notation as in Section 3.2. We let q be the smallest positive integer such that 2−q ≤ .

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

100

Algorithm FD
Step 0. (Initialization.) To each edge e of the graph assign the length u−1
e . For1 ≤ j ≤ k, let z j denote the flow that carries dj units along the shortest path between
st and tj under this metric, and let z denote the multicommodity flow (z 1 , z 2 , · · · , z k ).
1
Set x0 = 2λ(z)z and t = 0.

Step 1. Set y t =1+ 21
2λ(xt )

 t
x.

Step 2. (Frank-Wolfe procedure.) Write v 0 = y t and γ = γ(y t ).
For h = 0, 1, · · · Do:
A.h For 1 ≤ j ≤ k, let w h,j denote the solution to the minimum-cost flow problem
where we send γdj units of flow from sj to tj and each edge e has cost[∇Ψ(v h )]e and capacity ue . Let w h = (w h,1 , · · · , w h,k ) denote the resulting
multicommodity flow.
B.h Set v h+1 = (1 − σh )v h + σh w h where 0 ≤ σh ≤ 1 is chosen so as to minimizeΨ((1 − σh )v h + σh w h ).
h 2

Ψ(v )
t
h+1
C.h If Ψ(v h ) − Ψ(v h+1 ) < 128(Ψ(vandh )3 +m)k 2 , exit loop: set t ← t + 1, x ← v

go to Step 3.
End
Step 3. If λ(xt ) ≥ 1 − 2m
, or if Ψ(xt ) ≥ m2q+2 , terminate algorithm.

Otherwise to Step 1.
End.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

3.4.1

101

Analysis of Algorithm FD

Given γ > 0 denoteΨ∗ (γ) = min {Ψ(x) : x feasible with γ(x) = γ} .

(3.4.1)

In order to show the correctness of the algorithm, we first state the following
theorem whose proof is deferred:
Theorem 2 Consider an execution of the For loop in Step 2, with input y t . Suppose
the loop exits at iteration h. ThenΨ(v h ) ≤ 2Ψ∗ (γ(y t )).

(3.4.2)

Pending the proof of Theorem 2, the following sequence of lemmas establish thecorrectness and workload of the algorithm.
Lemma 3.4.1 Let r be a nonnegative integer. (i) Suppose γ < (1 − 2−(r+1) )Γ∗ . ThenΨ∗ (γ) < 2r+1 m. (ii) Suppose (1 − 2−r )Γ∗ ≤ γ. Then 2r − 1 ≤ Ψ∗ (γ).
Proof. (i) Consider x∗ feasible with throughput Γ∗ . Then u = Γγ∗ x∗ has throughput
γ and satisfies Ψ(u) ≤ mΨ( Γγ∗ ) < 2r+1 m. (ii) is proved in a similar way.
Lemma 3.4.2 Let 0 < τ < 1, and suppose that an iterate xt in Step 3 satisfiesτλ(xt ) ≥ 1 − 2m
. Thenγ(xt ) ≥ (1 − τ )Γ∗(3.4.3)

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

102

Proof. Assume by contradiction that (3.4.3) does not hold. Let x∗ be feasible witht) ∗throughput Γ∗ . Then u = γ(xx has throughput γ(xt ), and satisfiesΓ∗

λ(u) < 1 − τ.

(3.4.4)

Ψ(u) < m(τ −1 − 1).

(3.4.5)

As a result

Given the value of λ(xt ), we also haveΨ(xt ) ≥

2m
− 1,
τ

(3.4.6)

which together with (3.4.5) contradicts Theorem 2.
Lemmas 3.4.1(i) and 3.4.2 imply:
Corollary 3.4.3 Upon termination of Algorithm FD, we have a feasible flow with
-optimal throughput.
Now we turn to the complexity of Algorithm FD. First we have:∗ΓLemma 3.4.4 The vector x0 has throughput at least 2m.

Proof. This result is implied by weak linear programming duality, but a direct proof
follows. Denote by L∗ the sum, over all commodities, of the shortest path lengths in
Step 0. If x is any multicommodity flow we havemλ(x) ≥ L∗ .

(3.4.7)

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

103

Clearly λ(z) ≤ L∗ . Consequently,
1
1
≥
.
λ(z)
mλ(x)

(3.4.8)∗1
Thus λ(z)≥ Γm , and by construction of x0 , the result follows.

In what follows, for integral 0 ≤ r we will denote by Phase r of the algorithm the setof those iterations t with (1 − 2−r )Γ∗ ≤ γ(xt ) < (1 − 2−(r+1) )Γ∗ . Note that a given
Phase r might be empty, and that the algorithm might perform iterations of Phase r
with r > q.
Lemma 3.4.5 (i) For 0 < r ≤ q, the number of iterations t in Phase r is O(m).
(ii) The number of iterations in Phase 0 is O(m log m). (iii) The total number of
iterations in Phases q + 1, q + 2, · · · is O(m).
−rProof. By Lemma 3.4.2, for r ≥ 0 if γ(xt ) < (1 − 2−(r+1) )Γ∗ we have λ(xt ) < 1 − 24m .
Consequently,
1 − 1/2(1 − λ(xt ))γ(xt )
λ(xt )
1
1=(1 +)γ(xt )2
λ(xt )2−r> (1 +)γ(xt ).
8m

γ(xt+1 ) = γ(y t ) =(3.4.9)(3.4.10)(3.4.11)(i) Suppose 0 < r. By definition we started Phase r with throughput at least(1 − 2−r )Γ∗ . Thus, (3.4.11) implies that this Phase will perform O(m) iterations,as desired.

104

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM(ii)This follows as (i), using Lemma 3.4.4.
(iii) Consider an iteration t during Phase r with r > q. Since the algorithm has not yet−qterminated, we can replace (3.4.11) with the stronger condition γ(y t ) > (1+ 28m )γ(xt ),and again we obtain that there are altogether at most O(m) iterations in Phases
q + 1, q + 2, · · · .
The next lemmas analyze the complexity of each execution of Step 2.
Lemma 3.4.6 Consider an iteration t of Step 1, and let e be any edge. Thent

1 − λe (y t ) ≥ 1−λ2e (x ) .
t)
Proof. We have that λe (y t ) = 1+λ(xλ (xt ). Since by definition λe (xt ) ≤ λ(xt ), the2λ(xt ) e

result follows.
Lemma 3.4.7 Consider an iteration t of Step 1 during Phase r. Then Ψ(y t ) ≤O(m2min{r,q} ).
Proof. Suppose first that r > 0. Then by Theorem 2 and respectively, Lemma 3.4.1(for the case r ≤ q) and the second termination criterion in Step 3 (for the caser > q) we have that Ψ(xt ) ≤ O(m2min{r,q} ). To obtain the desired result, we willshow Ψ(y t ) = O(Ψ(xt )). To see this, consider the contributions of an edge e tott

λe (x )λe (y )
Ψ(xt ) and to Ψ(y t ). These are, respectively, 1−λand 1−λtt . By Lemma 3.4.6,
e (x )e (y )

the denominator in the second expression is at least half of that in the first. To
compare the numerators, note that since this is an iteration during Phase r > 0,

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

105

t

)
λ(xt ) ≥ 1 − 2−r ≥ 1/2 and consequently λe (y t ) = 1+λ(xλ (xt ) ≤ 23 λe (xt ). Hence2λ(xt ) eΨ(y t ) = O(Ψ(xt )), as desired.
Suppose instead that r = 0. We always have (again by Theorem 2 and by the
choice of x0 ) that Ψ(xt ) = O(m). If λ(xt ) ≥ 1/2 the Lemma follows as in the previous
paragraph. If instead λ(xt ) < 1/2, then each edge e will satisfy λe (y t ) ≤ 3/4 and thus
Ψ(y t ) = O(m).
Lemma 3.4.8 Let 0 ≤ r. The number of iterations of A.h-C.h in an execution of
Step 2 during Phase r is O(22 min{r,q} m2 k 2 ).
Proof. We will show that if more than O(22 min{r,q} m2 k 2 ) iterations h achieve
Ψ(v h+1 ) − Ψ(v h ) < −Ψ(v h )2,
128(Ψ(v h )3 + m)k 2(3.4.12)then we will reach a value of Ψ smaller than Ψ∗ , a contradiction. Thus, consider an
iteration h where (3.4.12) holds. Suppose first that
Ψ(v h )3 ≥ m.

(3.4.13)

In this case, the recursion (3.4.12) can be abstracted aszh+1 − zh ≤ −c1,
zh k 2where c > 0 is a constant. This recursion has the property that it reduces zh by afactor of 2 in O(zh2 k 2 ) iterations. Thus, using Lemma 3.4.7, we obtain that there areat most O(22 min{r,q} m2 k 2 ) iterations of A.h-C.h where (3.4.13) holds.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

106

In the remainder of the proof we handle the iterations with Ψ(v h )3 < m. If
r > 0 then using Lemma 3.4.1(ii) we conclude that each iteration satisfying (3.4.12)
decreases Ψ by Ω(1/(mk 2 )), and consequently there are at most O(m5/3 k 2 ) suchiterations (a tighter analysis is possible but not needed). This concludes the proof
if r > 0. Finally, if r = 0 then just as in the previous line we conclude that in∗Γ
at most O(m5/3 k 2 ) iterations we obtain Ψ(v h ) ≤ 1. By Lemma 3.4.4, γ(x0 ) ≥ 2m,

and a variation of the analysis in Lemma 3.4.1 shows that any time during Phase 0,1
1
Ψ∗ ≥ 2m. We conclude that there are at most O( 1/2m
) = O(m) iterations h with

Ψ(v h ) < 1.

Corollary 3.4.9 The total number of Frank-Wolfe iterations A.h over the course of
Algorithm FD is O(k 2 m3 −2 + m3 log m).
In the next section we prove Theorem 2. In Section 3.4.3 we show how to replace
the line-search in step B.h with a provably good stepsize rule that requires O(1) time.

3.4.2

Proof of Theorem 2

Consider an iteration h of the Frank-Wolfe procedure, corresponding to input vector
y t . For simplicity, write γ = γ(y t ). For 0 ≤ σ ≤ 1, write.
g(σ) = Ψ( (1 − σ)v h + σw h ) = Ψ( v h + σ(w h − v h )).

(3.4.14)

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

107

Then g is convex, g(0) = Ψ(v h ),g 0 (0) = [∇Ψ(v h )]T · (w h − v h ), and(3.4.15)

1 00
g (α)σ 2 ,2(3.4.16)g(σ) = g(0) + g 0 (0)σ +where 0 < α < σ (2nd order Taylor expansion). By choice of w h in A.h, the righthand side of (3.4.15) cannot decrease if we replace w h with any other feasible flow
with throughput γ. In particular, if we use a flow w ∗ with throughput γ and suchthat Ψ(w ∗ ) = Ψ∗ (γ), we obtaing(σ) − g(0) ≤ [∇Ψ(v h )]T · (w ∗ − v h )σ + 12 g 00 (α)σ 2 .

(3.4.17)

which, since Ψ is convex, implies:

1g(σ) − g(0) ≤ (Ψ∗ (γ) − Ψ(v h ))σ + g 00 (α)σ 2 .
2(3.4.18)

Next we will bound the quadratic term in (3.4.18). We haveg(α) =Xe

ψ(λe (v h ) + σ[λe (w h ) − λe (v h )]),

(3.4.19)(λe (w h ) − λe (v h ))2(1 − λe (v h ) − α[λe (w h ) − λe (v h )])3

(3.4.20)and consequently
g 00 (α) =

X
e

1
x 3
Since both v h and w h are feasible, (λe (w h )−λe (v h ))2 ≤ k 2 . Further, (1−x)
3 < 8( 1−x ) +8 for all 0 ≤ x < 1. Thus, from (3.4.20) we get thatg 00 (α) ≤ 8(k 2 [Ψ(v h + α(w h − v h )]3 + m).

(3.4.21)

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

108

If σ is chosen so that
Ψ(v h + σ(w h − v h )) ≤ Ψ(v h ),(3.4.22)then (3.4.21) and the convexity of Ψ imply g 00 (α) ≤ 8k 2 (Ψ(v h )3 +m), and substitutingin (3.4.18) we obtain:g(σ) − g(0) =2
Ψ(v h + σ(w h − v h )) − Ψ(v h ) ≤ (Ψ∗ (γ) − Ψ(v h ))σ + 8k 2 (Ψ3 (v h ) + m)σ(3.4.23).

Inequality (3.4.23) holds for any 0 ≤ σ ≤ 1 that satisfies (3.4.22). In Step B.h
we choose σ = σh so that the left-hand side of (3.4.23) is minimized. Instead, let
0 ≤ µ < 1 be the argument that minimizes the quadratic Q on the right-hand sideof (3.4.23) (note: a simple check verifies the stated bounds on µ). We would like toargue that Ψ(v h + σh (w h − v h )) − Ψ(v h ) ≤ Q(µ). This will follow if we can argue
that µ is such that (3.4.23) holds at σ = µ; or, in other words, that σ = µ satisfies(3.4.22).
But note that for 0 ≤ σ small enough (3.4.22) holds by choice of w h (w h − v h isa descent direction). If (3.4.23) holds for all 0 ≤ σ we are done. Otherwise (3.4.22)
holds for some closed interval [0, σ̃] and does not hold outside the interval. Since(3.4.23) holds with equality at 0, and Q0 (0) < 0, clearly µ < σ̃, as desired.

109

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

Thus, from (3.4.23) we get (after a calculation to minimize the quadratic)
Ψ(v h+1 ) − Ψ(v h ) =Ψ(v h + σh (w h − v h )) − Ψ(v h ) ≤ −(Ψ∗ (γ) − Ψ(v h ))2.
32k 2 (Ψ3 (v h ) + m)

(3.4.24)

If Ψ(v h ) > 2Ψ∗ (γ) we therefore obtain

Ψ(v h+1 ) − Ψ(v h ) ≤ −Ψ2 (v h ).
128k 2 (Ψ3 (v h ) + m)

(3.4.25)

The theorem is proved.

3.4.3

Choosing the stepsize

A key step in Algorithm FD is the choice of the stepsize σh in Step B.h. We expectthat from the point of view of computational effectiveness, it will be best to actually
carry out a numerical line-search to estimate σh . In the rest of this section, we show
how to instead choose a value for σh in O(1) time which nevertheless preserves thetheoretical properties of the algorithm.
One way to choose σ so as to minimize the quadratic Q(σ) introduced in the proofof Theorem 2, while assuming Ψ(v h ) > 2Ψ∗ (γ), or in other words, to minimize
Ψ(v h ).
Q̂(σ) = −σ + 8k 2 (Ψ3 (v h ) + m)σ 2
2(3.4.26)

whose decrease was used in the proof of Theorem 2 to bound the decrease in Ψ itself.

110

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

h

Ψ(v )
We have to show that if we alter Algorithm FD to use this stepsize (= 32k2 (Ψ
3 (v h )+m) )we still have a correct algorithm with the same complexity bound.
In fact, we proved that as long as Ψ(v h ) > 2Ψ∗ (γ), choosing σ so as to minimize
2

h

(v )
Q̂(σ) guarantees a decrease in Ψ of at least 128k2Ψ
. Thus the validity of the(Ψ3 (v h )+m)

algorithm is preserved. Further, Lemma 3.4.8 still applies (consider the first line inits proof). This shows that the complexity bound for Algorithm FD still applies.

3.4.4

Comments on Algorithm FD

In [48], the procedures and analysis are not pitched towards provably good approximation – this field existed only in embryonic form in 1971. This is the primary cause
for the differences between Algorithm FD and that in [48]. Nevertheless, the criticalidea in [48] is the key ingredient in Algorithm FD.
(i) The fact that the Frank-Wolfe iterations are minimum-cost flow problems is
only needed in the paragraph immediately following (3.4.20) so that we canguarantee λ(w h ) ≤ 1. If we were to replace the minimum-cost flow calls withshortest path calls, then instead we could only claim λ(w h ) ≤ ρ, where ρ is thewidth ([58] – in the context of this paper, the width of the problem is upperbounded by the ratio of the largest demand dj to the smallest capacity ue ).
This results in an increase in our complexity estimate by a factor of ρ2 . The“trick” of using minimum-cost flow calls to avoid the dependence on ρ is not

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

111

new: it is used in [58] and [50]. Recent algorithms for the maximum concurrent
flow problem avoid the dependence on ρ while only resorting to shortest pathcalls ([52], [45]). We conjecture that a more refined version of our algorithm
will achieve the same behavior.
(ii) The algorithmic outline in [48] uses a single Frank-Wolfe iteration in each Step
2. The algorithm described above can be adapted so as to use single Frank-Wolfe
iterations in each Step 2; it will be shown in Section 3.5 that the dependenceon  of the adapted algorithm also grows as O(−2 ). This adaptation entails afiner tuning of the multiplier in Step 1.
t))(iii) The magnification factor 1−1/2(1−λ(xused in Step 1 is as in [48], except that
λ(xt )

the choice of the parameter 1/2 is ours ([48] only states that “a proper tolerance”
should be used). In Step 0 our choice for flow z is precisely the same as the1
choice in the initialization step in [48]. We start with x0 = 2λ(z)z, rather thanz
with λ(z), so that we can guarantee that Ψ(x0 ) = O(m) (as opposed to O(m2 )).

(iv) Our proof of Theorem 2, up to inequality (3.4.18), follows fairly closely a similar
analysis in [48](pp. 128 - 130). The latter part of our proof follows [48] somewhat
less closely. In particular, [48] does not analyze the quadratic term in (3.4.18),
which is key to the overall complexity estimate. Finally, our analysis prior to
the proof of Theorem 2 is new.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

112(v) Algorithm FD is easily generalized to general packing linear programs [58].
(vi) [48] points out that [43] contains some ideas similar to those in the Flow Deviation method, although in less developed form. Courant [42] has been creditedwith one of the earliest uses of “potential” functions in the solution of systemsof equations.
(vii) Some pointers for the reader who decides to tackle [48]. Here we use the notationfrom [48]. (a) The algorithm in pp. 113-114 is geared towards finding a feasible
flow with throughput 1. This includes the exit condition in Step 1. (b) Step2, and the definition of the F D operator in page 110 inductively show the forn > 0 we always have σn < 1 and consequently REn+1 > REn . (c) There is anerror in the paper for the case n = 0: the exit condition should be σ0 /RE0 ≤ 1.

3.5

Single Frank-Wolfe iterate analysis

In this section we provide the analysis of the original version of the algorithm thatuses a single Frank-Wolfe iteration in each Step 2. First we will show how to tune
the coefficient in Step 1 and the termination criteria in Step 3 and then adjust thecorresponding lemmas to adopt these modifications. The updated FD algorithm now
takes form

113

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

Algorithm FD’
Step 0. (Initialization.) To each edge e of the graph assign the length u−1
e . For1 ≤ j ≤ k, let z j denote the flow that carries dj units along the shortest path between
st and tj under this metric, and let z denote the multicommodity flow (z 1 , z 2 , · · · , z k ).
1
Set x0 = 2λ(z)z and t = 0.
t

2 2

)) t
k
Step 1. Set y t = α̂xt ,where α̂ = 1−α(1−λ(xx and α = ((1−λ(xt ))2am− min{r,q} +m2 ) .
λ(xt )

Step 2. (Frank-Wolfe iteration.) Write v = y t and γ = γ(y t ).
A For 1 ≤ j ≤ k, let w j denote the solution to the minimum-cost flow problem
where we send γdj units of flow from sj to tj and each edge e has cost [∇Ψ(v)]eand capacity ue . Let w = (w 1 , · · · , w k ) denote the resulting multicommodity
flow.
B Set v̄ = (1 − σ)v + σh w where 0 ≤ σ ≤ 1 is chosen so as to minimizeΨ((1 − σ)v + σw).
C Set t ← t + 1, xt ← v̄.
2


Step 3. If t > a m k 2 , where a > 0 is some constant, and either λ(xt ) ≥ 1 − 2m
, or

Ψ(xt ) ≥ m2q+2 , terminate algorithm. Otherwise go to Step 1.
Note that the only change in Step 2 is a single application of Frank-Wolfe iterationand that Step 3 is adjusted to accommodate this modification. The rest of theFD algorithm remains unaffected. We know describe the necessary changes to several
lemmas in our analysis. We should specifically emphasize that the termination criteria

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

114

in Step 3 is chosen such that the result and the proof of Theorem 2 are not affected.
We however will need to modify lemmas 4.6 - 4.9 to analyze the new complexity
results and link the termination criteria of FD and FD’ algorithms.
Lemma 3.5.1 (4.60 ) (i) For 0 < r ≤ q, the number of iterations t in Phase r isO(22 min{r,q} m4 k 2 ). (ii) The number of iterations in Phase 0 is O(m log m). (iii) The
total number of iterations in Phases q + 1, q + 2, · · · is O(22 min{r,q} m4 k 2 ).
Proof. Write r = 2−r . By Lemma 3.4.2, for r ≥ 0 if γ(xt ) < (1 − 2−(r+1) )Γ∗ we haverλ(xt ) < 1 − 4m
. Consequently,

1 − α(1 − λ(xt ))γ(xt )
λ(xt )(1 − α)(1 − λ(xt )))γ(xt )> (1 +
λ(xt )(1 − α)(1 − λ(xt ))> (1 +)γ(xt )
λ(xt )2−3r)γ(xt ).
> (1 +
4m4 k 2γ(xt+1 ) > γ(y t ) >(3.5.1)(3.5.2)(3.5.3)(3.5.4)(i) Suppose 0 < r. By definition we started Phase r with throughput at least (1 −2−r )Γ∗ . Thus, (3.5.4) implies that this Phase will perform O(22r m4 k 2 ) iterations, as
desired. (ii)This follows as (i), using Lemma 3.4.4. (iii) Consider an iteration t during
Phase r with r > q. Since the algorithm has not yet terminated, we can replace (3.5.4)−q2twith the stronger condition γ(y t ) > (1 + 8mk
2 )γ(x ), and again we obtain that there

are altogether at most O(22q m4 k 2 ) iterations in Phases q + 1, q + 2, · · · .
The next lemmas analyze the complexity of each execution of Step 2.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

115− min{r,q}Lemma 3.5.2 (4.70 ) Ψ(y t ) ≤ Ψ(xt ) + 2 amk2

Proof. The increase of the potential on an edge e isλe (α̂xt )λe (xt )λe (xt )(1 − α̂)λ(xt )(1 − α̂)−=≤.
1 − λe (α̂xt ) 1 − λe (xt )(1 − λe (xt ))(1 − α̂λe (xt ))(1 − λ(xt ))(1 − α̂λ(xt ))ThereforeΨ(α̂xt ) − Ψ(xt ) ≤ m

λ(xt )(1 − α̂)
1−α2− min{r,q}=m=(1 − λ(xt ))(1 − α̂λ(xt ))α(1 − λ(xt ))mk 2

Lemma 3.5.3 (4.80 ) Consider an iteration t of Step 1 during Phase r. ThenΨ(y t ) ≤ O(m2min{r,q} ).
Proof. Suppose first that r > 0. Then by Theorem 2 and respectively, Lemma 3.4.1(for the case r ≤ q) and the second termination criterion in Step 3 (for the caser > q) we have that Ψ(xt ) ≤ O(m2min{r,q} ). By Lemma 3.5.2 Ψ(y t ) ≤ Ψ(xt ) + m andthe result follows.
Suppose instead that r = 0. We always have (again by Theorem 2 and by the
choice of x0 ) that Ψ(xt ) = O(m) and the proof follows in the same manner as inLemma 3.4.7
2

Lemma 3.5.4 (4.90 ) For any t > a m k 2 , where a > 0 is some constant, Ψ(v̄) ≤2Ψ∗ (γ(y t )).

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

116

Proof. First we will show that if more than O(22 min{r,q} m2 ) iterations achieve
Ψ(v̄) − Ψ(v) < −Ψ(v)2,
128(Ψ(v)3 + m)k 2(3.5.5)then we will reach a value of Ψ smaller than Ψ∗ , a contradiction. Thus, consider an
iteration where (3.5.5) holds. Suppose first that
Ψ(v)3 ≥ m.

(3.5.6)In this case, using Lemma 3.4.1 and Lemma 3.5.2 and with the appropriate choice ofparameter a the recursion (3.5.5) can be abstracted aszt+1 − zt ≤ −c1,
zt k 2where c is a constant. This recursion has the property that it reduces zt by a factor of
2 in O(k 2 zt2 ) iterations. Thus, using Lemma 3.5.3, we obtain that there are at mostO(22 min{r,q} m2 k 2 ) iterations of Step 1 – Step 3 where (3.5.6) holds.
Next we handle the iterations with Ψ(v)3 < m. If r > 0 then using Lemma 3.4.1(ii)we conclude that each iteration satisfying (3.5.5) decreases Ψ by Ω(1/(mk 2 )), andconsequently there are at most O(m5/3 k 2 ) such iterations. This concludes the proof if
r > 0. If r = 0 then just as in the previous line we conclude that in at most O(m5/3 k 2 )∗Γiterations we obtain Ψ(v) ≤ 1. By Lemma 3.4.4, γ(x0 ) ≥ 2m
, and a variation of the1
analysis in Lemma 3.4.1 shows that any time during Phase 0, Ψ∗ ≥ 2m
. We conclude
1that there are at most O( 1/2m
) = O(m) iterations with Ψ(v) < 1.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

117

2

Thus, after at most a m k 2 iterations of Step 1 – Step 3 we will achieve Ψ(v̄) −2and therefore by Theorem 2 Ψ(v̄) ≤ 2Ψ∗ (γ(y t )). In theΨ(v) < − 128k2Ψ(v)(Ψ(v)3 +m)
2remainder of the proof we will handle iterations with t > a m k 2 . Let t be the first
iteration such that after execution of Step 1 Ψ(y t ) ≥ 2Ψ∗ (γ(y t )). Then
2

Ψ(v)
t
Ψ(v̄) − Ψ(v) < − 128(Ψ(v)
3 +m)k 2 and by Lemma 3.5.2 every iteration increases Ψ(x )− min{r,q}by at most 2 mk2− min{r,q}. Therefore Ψ(y t ) ≤ 2Ψ∗ (γ(y t )) + 2 mk2

and one execution of

Step 2 will achieve Ψ(v̄) ≤ 2Ψ∗ (γ(y t )).
Corollary 3.5.5 (4.100 ) The total number of Frank-Wolfe iterations A over thecourse of Algorithm FD is O(k 2 m4 −2 + m4 k 2 log m).
Thus we have demonstrated that using a single Frank-Wolfe iteration in Step 2
we can achieve the same degree of convergence by tuning the magnification factor in
Step 1

3.6

Future research

In this section we propose some ideas for further enhancement of the algorithm. [55]demonstrated that any Dantzig-Wolf type algorithm for a general fractional packing
problem with high probablity will requier at least Ω(−2 ) iterations. Therefore, inorder to improve the dependency on  we need to move away from the Dantzig-Wolf
framework and develop a different approach to generating a sequience of approximate

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

118

solutions. In what follows we present an approach that can potentially lead to adecrease in the running time dependency on .
Recall that the central operation of the algorithm is to sequentially minimize thefirst-order approximation of the potential function over the set of multicommodity
constraints. The length of the step taken towards the optimum of the gradient essentially reflects the ”quality” of the approximation and produces a new vector within
the neighborhood where the approximation is accurate. The natural modification ofthis scheme would be to replace the first order approximation with a tighter one,
thereby expanding the neighborhood and allowing for longer steps. This in turn will
require a smaller number of iterations and thus will lead to a faster convergence of
the algorithm. The obvious drawback of this approach is the requirement to optimize
a more complex, nonlinear function over the same set of multicommodity constraints.
Further, the nonlinearity of the function will no longer allow for a simple decomposition of the resulting optimization problem. However, if a method to solve thisproblem efficiently is developed, this approach will lead to an overall improvement of
the algorithm performance.
In what follows we consider the second order Taylor approximation of the potentialfunction and offer the analysis of the complexity of the resulting algorithm provided
an efficient method to solve the quadratic multicommodity flow problem is known.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

3.6.1

119

Notations

For some flow f and some constant θ > 0 let fθ be any feasible flow such thatλe (fθ ) = λe (f ) + θte , for some vector t = {te }e with |te | ≤ 1. In other words, fθ
belongs to θ-neighborhood of f with respect to λe (f ).
Let g = f + θt where t is as defined above for some θ-neighborhood of f . Then
θ 3 000
Ψ(g) = Ψ(f ) + θ∇Ψ(f )t + 0.5θ tH(Ψ(f ))t + Ψ (f + αt)t362where α < θ.
Denote.
Q(t) = Ψ(f ) + θ∇Ψ(f )t + 0.5θ 2 tH(Ψ(f ))tand. θ 3 000
E(t) = Ψ (f + αt)t36Then the quadratic multicommodity flow problem in θ-neighborhood of f is defined asmin {Q(t) : g a feasible flow such that g = f + θt}(3.6.1)

Here we define a feasible flow as a flow that satisfies all the constraints except for thejoint capacity constraint. Note that the form of Q(t) does not allow us to decompose(3.6.1) into single-commodity problems.
Let f be a flow such that Ψ(f )  Ψ∗ , then there is a point f¯ in the θ-neighborhoodof f such that Ψ(f¯)  Ψ(f ). Moreover, if the approximation error E(t) in this

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

120neighborhood does not exceed the possible decrease in the potential function, we canconclude that there is a vector t, such that Q(t)  Q(0).

3.6.2

Algorithm

We can now state the new algorithm
Algorithm FDQ
Step 0. (Initialization.)t)) t
Step 1. Set y t = 1−1/2(1−λ(xx . Step 2. Write v 0 = y t and γ = γ(y t ).
λ(xt )For h = 0, 1, · · · Do:
A.h Let w h denote the solution to (3.6.1) with f ≡ v h for some θ to be definedlater.
B.h Set v h+1 = w h .
C.h If Ψ(v h ) − Ψ(v h+1 ) < 7k1.5 (4√2Ψ+√m) , exit loop: set t ← t + 1,
xt ← v h+1 and go to Step 3.
End
Step 3. If λ(xt ) ≥ 1 − 2m
, or if Ψ(xt ) ≥ m2q+2 , terminate algorithm.

Otherwise go to Step 1.
End.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

121

Note that the only difference of FD and FDQ algorithms is the optimization
problem in Step B.h and the termination criteria in C.h. Thus the only additional
analysis required for the FDQ algorithm is the number of iterations in Step 2 andthe validity of its termination criteria.

3.6.3

Analysis

Denote δ = k1.5 (4√1Ψ+√m) and let f ∗ be such that Ψ(f ∗ ) = Ψ∗ (γ) and Ψ(f ) > 2Ψ∗ .
Let Ψ(f ) − Ψ(fθ ) = δ.
Consider the example in Figure 3.1. This example illustrates the main idea of ourapproach, which is the dependence of the decrease in the potential function on thelength of the step θ.

Figure 3.1: Potential function decrease in θ-neighborhood

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

122

More precisely, for every edge eθδ=∗|λe (f ) − λe (f )|Ψ(f ) − Ψ∗Since both f and f ∗ are feasible flows,|λe (f ) − λe (f ∗ )| ≤ k, hence
θ≤

δ
2δ
≤
k 0.5 (Ψ(f ) − Ψ∗ )k 0.5 Ψ(f )

or
θ≤

1
10k 0.5 Ψ(f )1.5

if Ψ(f ) > 2m andθ≤1
10k 0.5 m1.5

otherwise.
Lemma 3.6.1 Let f be a feasible flow and g be a flow such that λ(f ) ≥ 3/5 and1|λe (g) − λe (f )| = θ for every edge e, where θ = 10k0.5 Ψ(fif Ψ(f ) > 2m and)1.5

θ = 10k0.51 m1.5 otherwise. ThenΨ(g) ≤ 3Ψ(f )Proof:

Denote xe ≡ λe (f ), ye ≡ λe (g), x ≡ λ(f ), y ≡ λ(g), then for any edge esuch that λe (g) > λe (f )

ψ(λe (g))
ye (1 − xe )(xe + θ)(1 − xe )===ψ(λe (f ))xe (1 − ye )xe (1 − xe − θ)CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

=1+

123

θ
xe (1 − xe − θ)which is clearly an increasing function of xe for any xe > 3/5. Then for any edge eψ(λe (g))ψ(λ(g))≤.
ψ(λe (e))ψ(λ(f ))Hence
1
1. If Ψ(f ) > 2m, we have θ = 10k0.5 Ψ(f≤ ( 1−x)1.5 . For any x ≥ 3/5
)1.5
10x

(

1 − x 1.5
) ≤ (1 − x)/2.
10x

Hence we get1+θ(1 − x)/2
1≤1+=1+ ≤3x(1 − x − θ)x((1 − x) − (1 − x)/2)xx
1
2. If Ψ ≤ 2m, then 1−x≤ 2m, or x ≤ 1 − 2m
. Then1+1
θ1
5
10m1.5
≤1+
≤1+
≤1+
≤3
1
1
0.5
x(1 − x − θ)
3/5(2m − 1)5m0.5x( 2m − 10m1.5 )ThenΨ(g) =XeX
λe (g)λe (f )≤
3
= 3Ψ(f )1 − λe (g)
1
−
λ
(f
)
e
e


Observation 1 If at some point of the algorithm we achieve λ(f ) ≤ 3/5 then we canterminate Step 2 without disrupting the course of the algorithm. This follows fromthe fact that in this case the termination criteria in Step 3 is not met, and hence theproof of Lemma 3.4.2 is not affected.

124

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

We can now bound the error term E(t) for the quadratic approximation of Ψ(f ).
1
Lemma 3.6.2 For θ = 10k0.5 Ψ(fif Ψ(f ) > 2m and θ = 10k0.51 m1.5 otherwise, we
)1.5

have
5
E(t) ≤ δ
7
where δ =Proof:

k 1.5 (4√1√ )
Ψ(f )+ m

Let v = f + αp. We have000

Ψ (v)t3 =

X
e

X
1
t3e
≤
4
(1 − λe (v))
(1 − λe (v))4e(3.6.2)

1
x 4≤ 8( 1−x) + 8 for all 0 ≤ x < 1. Thus, we get(1−x)4000

Ψ (v) ≤ 8([Ψ(v)]4 + m).

(3.6.3)Using Lemma 3.6.1 we obtain000

1/6Ψ (λ(v))t3 ≤ θ 3 ∗ 8/6([Ψ(v)]4 + m)

(3.6.4)≤ θ 3 ∗ 8/6([3Ψ(f )]4 + m)
= 108θ 3 Ψ(f )4 + 8/6θ 3 m
Then if Ψ(f ) > 2m,
108Ψ(f )48m
1108θ Ψ(f ) + 8/6θ m ≤+≤≤ 1.5≤ 5/7δ
1.5
4.5
1.5
4.5
1000k Ψ(f )1000k Ψ(f )7k Ψ(f )0.5
3

4

3

If Ψ(f ) < 2m then108θ 3 Ψ(f )4 + 8/6θ 3 m ≤

108m4
8m
1
+
≤ 1.5 0.5 ≤ 5/7δ
1.5
4.5
1.5
4.5
1000k m1000k m7k m

125

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM


Thus we have also proved the following result.
Theorem 3 Consider an execution of the For loop in Step 2, with input y t . Suppose
the loop exits at iteration h. ThenΨ(v h ) ≤ 2Ψ∗ (γ(y t )).

(3.6.5)

We now analyze the complexity of Step 2 of the new algorithm.
Lemma 3.6.3 Let 0 ≤ r. The number of iterations of A.h-C.h in an execution of
Step 2 during Phase r of the algorithm FDQ is O(21.5 min{r,q} m1.5 k 1.5 ).
Proof. We will show that if more than O(21.5 min{r,q} m1.5 k 1.5 ) iterations h achieve
Ψ(v h+1 ) − Ψ(v h ) < −2/7δ,(3.6.6)then we will reach a value of Ψ smaller than Ψ∗ , a contradiction.
Thus, consider an iteration h where (3.6.6) holds. Suppose first that2Ψ(v h ) ≥ m.

(3.6.7)In this case, the recursion (3.6.6) can be abstracted as1
zh+1 − zh ≤ −c 0.5 1.5 ,zh kwhere c > 0 is a constant. This recursion has the property that it reduces zh by afactor of 2 in O(zh1.5 k 1.5 ) iterations. Thus, using Lemma 3.4.7, we obtain that there
are at most O(21.5 min{r,q} m1.5 k 1.5 ) iterations of A.h-C.h where (3.6.7) holds.

CHAPTER 3. MAXIMUM CONCURRENT FLOW PROBLEM

126

In the remainder of the proof we handle the iterations with 2Ψ(v h ) < m. If r > 0then using Lemma 3.4.1(ii) we conclude that each iteration satisfying (3.6.6) decreases
Ψ by Ω(1/m0.5 k 1.5 ), and consequently there are at most O(m1.5 k 1.5 ) such iterations.
This concludes the proof if r > 0. Finally, if r = 0 then just as in the previous line
we conclude that in at most O(m1.5 k 1.5 ) iterations we obtain Ψ(v h ) ≤ 1. By Lemma∗Γ
3.4.4, γ(x0 ) ≥ 2m
, and a variation of the analysis in Lemma 3.4.1 shows that any
11. We conclude that there are at most O( 1/2m
) = O(m)time during Phase 0, Ψ∗ ≥ 2m

iterations h with Ψ(v h ) < 1.

Corollary 3.6.4 The total number of iterations A.h over the course of Algorithm
FD is O(m2.5 k 1.5 −1.5 T ) where T is the complexity to optimize the quadratic multicommodity flow problem in θ-neighborhood of f .

Thus it is evident that the algorithm developed in 1971 by Fratta, Gerla andKleinrock not only shows a performance comparable with the most recent algorithms,
but also has a strong potential for further improvement.

Bibliography[1] D. Alevras, M. Grotschel, R. Wessaly, Capacity and survivability models fortelecommunications networks, ZIB preprint SC 97-24, ZIB, Berlin, June 1997.
[2] D. Alevras, M. Grotschel, R. Wessaly, Cost-efficient network synthesis from leasedlines, Annals of Operations research, 76, pp. 1-20, 1998.
[3] N. Bakhvalov, N. Zhidkov, G. Kobelkov, Numerical methods, Physmathlit, 2000.
[4] D. Bienstock, Potential function methods for approximately solving linear programming problems: Theory and Practice, CORE Lecture Series, 2001.
[5] D. Bienstock, O. Gunluk, S. Chopra and C.Y. Tsai, Mininum cost capacity installation for multicommodity flows, Mathematical programming, 81, pp. 177-199,
1998.
[6] D. Bienstock, O. Gunluk, Capacitated network design - polyhedral structure andcomputation, Mathematical programming.

127

BIBLIOGRAPHY

128[7] D. Bienstock, G. Muratore, Strong inequalities for capacitated survivable network
design problems, Mathematical programming, 89, no. 1, Ser. A, pp 127-147, 2000.
[8] D. Bienstock, I. Saniee, ATM network design: traffic models and optimizationbased heuristics, Telecommunication systems, 16:3,4, pp. 399-421, 2001.
[9] S. Chopra, I. Gilboa, S.T. Sastry, Source sink flows with capacity installation inbatches, Discrete Applied Mathematics, 85, pp. 165192, 1998.
[10] O. Gunluk, A branch-and-cut algorithm for capacitated network design problems,Mathematical programming 86, 17-39, 1999.
[11] M. Garey, D. Johnson, Computers and Intractability: A Guide to the Theory of
NP-Completeness, Freeman, 1979.
[12] P. Gill, W. Murray, M.Wright, Numerical linear algebra and optimization,
Addison-Wesley, 1991.
[13] P. Gill, W. Murray, M.Saunders, User’s guide for SNOPT 5.3[14] R. D. Davis, K. Kumaran, G. Liu, I. Saniee, SPIDER: A simple and flexible toolfor design and provisioning of protected lightpaths in optical networks, Bell Labs
Technical Journal, January-June 2001.
[15] CLAPACK library, www.netlib.org/clapack

BIBLIOGRAPHY

129[16] ILOG CPLEX, www.ilog.com/products/cplex[17] L. Fadeev, V. Fadeev, Numerical methods in linear algebra, Physmathgiz, 1960.
[18] J. Kelley, The cutting-plane method for solving convex programs, Journal of
SIAM, 8, pp. 703-712, 1960.
[19] K. Kiwiel, An aggregate subgradient method for nonsmooth convex minimization, Mathematical programming, 27, pp. 320-341, 1983.
[20] S. G. Lanning, D. Mitra, Q. Wang, M. H. Wright, Optimal planning for opticaltransport networks, Philosophical Transactions Royal Society London A, Vol.358,
No.1773, pp.2183-2196, August 2000.
[21] C. Lemarechal, A. Nemirovskii, Y. Nesterov, New variants of bundle methods,
Mathematical programming, 69, pp. 111-147, 1995.
[22] D. Luenberger, Investment science, New York : Oxford University Press, 1998.
[23] D. Luenberger, Introduction to linear and nonlinear programming, AddisonWesley, 1989.
[24] A. Malcev, Linear algebra, Nauka, 1970.
[25] T.L. Magnanti, P. Mirchandani, Shortest paths, single origin-destination network
design, and associated polyhedra, Networks, 23, pp. 103121, 1993

BIBLIOGRAPHY

130[26] T.L. Magnanti, P. Mirchandani, R. Vachani, Modeling and solving the twofacility capacitated network loading problem, Operations Research, 43, pp. 142157,
1995[27] G. Mohan, S. Murthy, A. Somani, Efficient algorithms for routing dependable
connections in WDM optical networks, IEEE/ACM Transactions on Networking,vol. 9, no. 5, pp. 553-566, 2001.
[28] J. T. Moy, OSPF: Anatomy of an Internet Routing Protocol, Addison-Wesley,
1998.
[29] A. Nemirovskii, D. Yudin, Problem complexity and method efficiency in optimization, Wiley-Interscience, 1983.
[30] C. Papadimitriou, K. Steiglitz, Combinatorial optimization, Dover, 1998.
[31] U. Paul, P. Jonas, D. Alevras, M. Grotschel, R. Wessaly, Survivable mobile phone
network architectures: models and solution methods, ZIB preprint SC 96-48, ZIB,
Berlin, December 1996.
[32] Y. Pochet, L.A. Wolsey, Integer knapsack and flow covers with divisible coefficients: Polyhedra, optimization, and separation, Discrete Applied Mathematics,
59, pp. 5774, 1995

BIBLIOGRAPHY

131[33] S. Ramamurthy, B. Makherjee, Survivable WDM mesh networks,Part I - Protection, IEEE INFOCOM, vol. 2, pp. 744-751, 1999.
[34] S. Ramamurthy, B. Makherjee, Survivable WDM mesh networks,Part II Restoration, IEEE ICC, vol. 3, pp. 2023-2030, 1999.
[35] I. Saniee, Optimal routing design in self-healing communications networks, Int.
Trans. Op. Res., vol 3, no. 2, pp. 187-195 1996.
[36] A. Schrijver, Theory of linear and integer programming, Wiley-Interscience,1999.
[37] M. Sridharan, A. Somani, M. Salapaka, Approaches for capacity and revenue
optimization in survivable WDM networks, Journal of High Speed Networks, no. 2,
pp. 109-125 2001.
[38] M. Sridharan, M. Salapaka, A. Somani, Operating mesh-survivable WDM
transport networks, SPIE International Symposium on SPIE Terabit Optical
Netwroking: Terabit Optical Networking, pp. 113-123 2000.
[39] R. Vanderbei, LOQO User’s manual, Statistics and Operations Research, Technical report No SOR-97-08, 1997.
[40] J. Veerasamy, S. Venkatesan, J. C. Shah, Spare capacity assignment in telecom
networks using path restoration, MASCOTS, pp. 370-375 1995.

BIBLIOGRAPHY

132[41] D. Bienstock, Approximately solving large-scale linear programs. I. Strengthening lower bounds and accelerating convergence, CORC Report 1999-1, Columbia
University. (Extended abstract: Proc 11th. Ann. ACM-SIAM Symposium on Discrete Algorithms, January 2000).
[42] R. Courant, Variational methods for the solution of problems of equilibrium andvibration, Bull. Amer. Math. Soc. 49 (1943), 1 – 43 .
[43] S.C. Dafermos and F.T. Sparrow, The traffic assignment problem for a general
network, Journal of Research of the National Bureau of Standards - B, 73B (1969).
[44] G.B. Danzig and P. Wolfe, The decomposition algorithm for linear programming,
Econometrica 29 (1961), 767 – 778.
[45] L. Fleischer, Approximating fractional multicommodity flow independent of thenumber of commodities, SIAM J. Disc. Math. 3 (2000) 505 – 520.
[46] A.V. Fiacco and G.P. McCormick, Nonlinear Programming: Sequential Unconstrained Optimization Techniques, Wiley, New York (1968).
[47] M. Frank and P. Wolfe, An algorithm for quadratic programming, Naval Res.
Logistics Quarterly 3 (1956), 149 – 154.
[48] L. Fratta, M. Gerla and L. Kleinrock, The flow deviation method: an approachto store-and-forward communication network design, Networks 3 (1971), 97 – 133.

BIBLIOGRAPHY

133[49] A.V. Goldberg, J.D. Oldham, S. Plotkin and C. Stein, An implementationof a combinatorial approximation algorithm for minimum-cost multicommodity
flow, in Integer Programming and Combinatorial Optimization (Bixby, Boyd, RiosMercado, eds.) 1412 (1998) 338 - 352.
[50] M.D. Grigoriadis and L.G. Khachiyan Fast approximation schemes for convex
programs with many blocks and couping constraints, SIAM Journal on Optimization 4 (1994) 86 – 107.
[51] M.D. Grigoriadis and L.G. Khachiyan An exponential-function reduction method
for block-angular convex programs, Networks 26 (1995) 59 – 68.
[52] N. Garg and J. Könemann, Faster and simpler algorithms for multicommodity
flow and other fractional packing problems, Proc. 39th Ann. Symp. on Foundationsof Comp. Sci. (1998) 300 – 309.
[53] D. Karger and S. Plotkin, Adding multiple cost constraints to combinatorial
optimization problems, with applications to multicommodity flows, Proc. 27th Ann.
ACM Symp. on Theory of Computing (1995) 18 – 25.
[54] P. Klein, S. Plotkin, C. Stein and E. Tardos, Faster approximation algorithms forthe unit capacity concurrent flow problem with applications to routing and finding
sparse cuts, Proc. 22nd Ann. ACM Symp. on Theory of Computing (1990),310 –321.

BIBLIOGRAPHY

134[55] P. Klein and N. Young, On the number of iterations for Dantzig-Wolfe optimization and packing-covering approximation algorithms.
[56] T. Leighton, F. Makedon, S. Plotkin, C. Stein, E. Tardos and S. Tragoudas, Fast
approximation algorithms for multicommodity flow problems, Proc. 23nd Ann.
ACM Symp. on Theory of Computing (1991), 101-111.
[57] T. Leighton and S. Rao, An approximate max-flow min-cut theorem for uniform
multicommodity flow problems with applications to approximation algorithms,
Proc. FOCS 29 (1988), 422 – 431.
[58] S. Plotkin, D.B. Shmoys and E. Tardos, Fast approximation algorithms for fractional packing and covering problems, Proc. 32nd Annual IEEE Symp. On Foundations of Computer Science, (1991), 495 – 504.
[59] T. Radzik, Fast deterministic approximation for the multicommodity flow problem, Proc. 6th ACM-SIAM Symp. on Discrete Algorithms (1995), 486 – 492.
[60] F. Shahrokhi and D.W. Matula, The maximum concurrent flow problem, Journalof the ACM 37 (1991), 318 – 334.
[61] J. Villavicencio and M.D. Grigoriadis, Approximate Lagrangian decomposition using a modified Kamkarkar logarithmic potential, in Network Optimization

BIBLIOGRAPHY

135(Pardalos and Hager, eds.) Lecture Notes in Economics and Mathematical Systems
450 Springer-Verlag, Berlin (1995), 471 – 485.
[62] N. Young, Randomized rounding without solving the linear program, in The
maximum concurrent flow problem, Proc. 6th ACM-SIAM Symp. on Discrete Algorithms (1995), 170 – 178.