Call Admission and Routing in ATM Networks Based on Virtual Path Separation Johan M de Kock Department of Computer Science University of Stellenbosch 7600 Stellenbosch, South Africa
[email protected] Honours Project 1998
Contents List Of Figures
iv
List Of Tables
v
Preface
vi
Acknowledgement
vii
1 Introduction
1.1 Outline of the Project . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Loss Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1 1
2 Computing an Optimal Virtual Path Connection Network by Simulated Annealing 3 2.1 Introduction . . . . . . . . . . . . . . . . . 2.2 Dynamic Recon guration . . . . . . . . . 2.3 Simulated Annealing . . . . . . . . . . . . 2.3.1 Combinatorial Optimization . . . . 2.3.2 Simulated Annealing Theory . . . 2.4 The Simulated Annealing Implementation 2.5 Results . . . . . . . . . . . . . . . . . . . . 2.6 Conclusion . . . . . . . . . . . . . . . . .
3 The Reduced Load Approximation
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. 3 . 4 . 6 . 6 . 6 . 9 . 10 . 11
3.1 A Product-form Loss Network Model . . . . . . . . . . . . . . . . 3.2 The Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 The Blocking Function Eijk . . . . . . . . . . . . . . . . . 3.3.2 Additional Performance Measures . . . . . . . . . . . . . 3.3.3 Solving the Erlang Fixed-point Equations . . . . . . . . . 3.4 Oscillation of the Approximate Link Blocking Probabilities Lm ij;k 3.5 Computer Implementation . . . . . . . . . . . . . . . . . . . . . . 3.6 Numerical Examples and Observations . . . . . . . . . . . . . . . 3.6.1 A Simple Example . . . . . . . . . . . . . . . . . . . . . . 3.6.2 A Larger Network . . . . . . . . . . . . . . . . . . . . . . ii
12
12 13 14 14 15 15 17 22 22 22 23
CONTENTS
iii
4 Route Optimization Using Virtual Subnetworks
29
A B C D
53 56 57 59
4.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Convex Approximation . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Homogeneous Bandwidth Demands . . . . . . . . . . . . . 4.2.2 Heterogeneous Bandwidth Demands . . . . . . . . . . . . 4.3 Re ned Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Algorithmic Implementation . . . . . . . . . . . . . . . . . . . . . 4.4.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Approximation Using Fixed Load Sharing . . . . . . . . . 4.4.3 Re ned Model . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Approximation Using Fixed Load Sharing . . . . . . . . . 4.5.2 Re ned Model . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.1 Farago's Three-Node Networks . . . . . . . . . . . . . . . 4.6.2 An Example Network for which Repeated Substitution Fails 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Convex and Concave Functions The Erlang-B Formula The Basic Stochastic Knapsack Non-linear Optimizers and Non-linear Equation Solvers
29 32 32 33 33 35 35 36 40 41 41 42 46 46 49 52
D.1 Non-linear Optimizers . . . . . . . . . . . . . . . . . . . . . . . . 59 D.2 Non-linear Equation Solvers . . . . . . . . . . . . . . . . . . . . . 59
Bibliography Index
61 63
List of Figures 1.1 A Loss System . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.1 A 49-node Manhattan Grid . . . . . . . . . . . . . . . . . . . . .
4
3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13
A Route Connecting Nodes i and j . . . . . . . Three Routes in a Loss Network . . . . . . . . A Route on which Oscillation Occurs . . . . . . Oscillation of Lm ij;1 . . . . . . . . . . . . . . . . The E ect of Restarting . . . . . . . . . . . . . The E ect of Damping . . . . . . . . . . . . . . A Route on which Oscillation Seems to Occurs Seeming Oscillation of Lm ij;1 . . . . . . . . . . . Convergence of Lm . . ............. ij;1 RLA's route Structure . . . . . . . . . . . . . A Simple 3-node Network . . . . . . . . . . . . The 4 Routes in the 3-node Network . . . . . . The 5-node Network . . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
12 13 17 19 20 21 21 22 23 24 24 24 25
4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9
Relation Between r, r0 and j . . . . . . The Relation Between r, r0 , etc. . . . . . VSNC 's route Structure . . . . . . . . Main Steps of VSNR . . . . . . . . . . VSNR's Principal Data Structures . . . Farago's 3-node Network . . . . . . . . . The Routes in Farago's 3-node Network A Simple 5-node Network . . . . . . . . A Route in the 5-node Network . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
32 37 41 43 44 47 47 49 50
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
A.1 A Convex Function . . . . . . . . . . . . . . . . . . . . . . . . . . 53 A.2 A Concave Function . . . . . . . . . . . . . . . . . . . . . . . . . 54 D.1 The Contours in the xy-plane . . . . . . . . . . . . . . . . . . . . 60
iv
List of Tables 2.1 Results for the Manhattan Grid . . . . . . . . . . . . . . . . . . . 11 2.2 Simulated Annealing Moves for the Manhattan Grid . . . . . . . 11 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11
Lmij;1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 The O ered Loads ijk for the 3-node Network . . . . . . . . . . . 23 The Revenues for the 3-node Network . . . . . . . . . . . . . . . The Link Blocking Probabilities Bijk for the 3-node Network . . . The Route Blocking Probabilities Lkr for the 3-node Network . . The Network Blocking Probabilities Bek for the 3-node Network . The O ered Loads ijk for the 5-node Network . . . . . . . . . . . The Routes in the 5-node Network . . . . . . . . . . . . . . . . . The Link Blocking Probabilities Bijk in the 5-node Network . . . The Route Blocking Probabilities Lkr in the 5-node Network . . The Network Blocking Probabilities Bek and Revenue Rate in the 5-node Network . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 The O ered Loads pk for Farago's 3-node Network . . . . . . . . 4.2 VSNC Revenue Rates for Farago's 3-node Networks . . . . . . . 4.3 Revenue Rates After Applying the Reduced Load Approximation to the VSNC Results for Farago's 3-node Networks . . . . . . . 4.4 VSNR Revenue Rates for Farago's 3-node Networks . . . . . . . 4.5 VSNR Revenue Rates for Farago's 3-node Network (Computed from VSNC Results) . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Revenue Rates for Farago's 3-node Networks . . . . . . . . . . .
v
23 25 25 26 26 26 27 28
28 48 48 48 49 49 50
Preface This document contains the description of an honours project performed in the department of Computer Science at the University of Stellenbosch. The supervisor for the project was prof Anthony E Krzesinski.
Johannes Marthinus de Kock B.Sc. Stellenbosch, South Africa 20 January 1999
vi
Acknowledgement This work was performed within the Siemens-Telkom Centre of Excellence for ATM & Broadband Networks and their Applications and is supported by grants from the South African Foundation for Research Development, Telkom SA Limited and Siemens Telecommunications.
vii
Chapter 1
Introduction This honours project investigates the reduced load approximation and two methods of computing virtual path connection networks (VPCNs) from sparse physical networks.
1.1 Outline of the Project Dynamic Recon guration is a network management control which constructs an optimal virtual path connection network from a sparsely connected physical network. In chapter 2 we de ne the computation of an optimal VPCN in terms of a non-linear optimization problem, which we solve using the combinatorial optimization method of simulated annealing. Chapter 3 is devoted to the reduced load approximation { a method of approximating link blocking probabilities in product-form loss networks. In chapter 4 we discuss the construction of an optimal VPCN by means of virtual subnetworks .
1.2 Loss Networks An important concept in this project is that of a loss network. A loss network is a network that can be modelled as a loss system . Ross [18] de nes a loss system as follows: \A loss system is a collection of resources to which calls, each with an associated holding time and class arrive at random instances. An arriving call either is admitted into the system or is blocked and lost; if the call is admitted, it remains in the system for the duration of its holding time. The admittance decision is based on the call's class and the system's state." Figure 1.1 is a possible way of viewing a loss system (adapted from a diagram in [18]). The vector n = (n1 ; : : : ; nK ) denotes the state of the loss system (where nk denotes the number of class-k calls in the system). In chapter 3 we come across a product-form loss network. The term product-form refers to the equilibrium distribution of such a loss network (see Theorem 5.1 in [18]):
P (X = n) = G1
K n Y k
k
k=1 nk !
1
where n 2 S ;
(1.1)
2
Chapter 1. Introduction
arrival
collection of resources
admitted
departure
lost
Figure 1.1: A Loss System
G=
K n XY k
k
n2S k=1 nk !
where K is the number of call classes, k the o ered load of class-k calls (the arrival rate multiplied by the average holding time) and S the state space of the loss network. Any loss network whose equilibrium distribution is given by (1.1) is referred to as a product-form loss network.
Chapter 2
Computing an Optimal Virtual Path Connection Network by Simulated Annealing Dynamic Recon guration is a network management control which constructs an optimal virtual path connection network (VPCN) from a sparsely connected physical network. We de ne the computation of an optimal VPCN in terms of a non-linear optimization problem, which we solve using the combinatorial optimization method of simulated annealing. We compare simulated annealing with a deterministic algorithm that solves the same model.
2.1 Introduction Virtual paths (VPs) and virtual path connections (VPCs) are important concepts in ATM (Asynchronous Transfer Mode) networks. A VP recognizes the distinct identity of a trac stream between two nodes. Cells belonging to a particular VP are identi ed by a common, xed VP identi er (VPI). A VPC consists of a series of concatenated VPs (possibly with di erent VPIs) and speci es a route to be traversed by a trac stream from an originating node through a number of intermediate nodes to a destination node. Using VPCs allows for faster setup of new connections (along prede ned routes) and rapid movement of trac (ATM cells) with minimal processing at intermediate nodes (which is related to the switching costs at each node). The virtual path connection network (VPCN) design problem can be formulated as follows: given a set of nodes which are interconnected by a physical network, and o ered tracs for each service class, nd the end-to-end routes and bandwidths which maximize network utility. For simplicity, we express network utility as carried trac, but it is emphasized that a large variety of optimization criteria can be deployed. The contents of this chapter has been published in [5].
3
Chapter 2. Computing an Optimal VPCN by Simulated Annealing
4020 140 140 860 4540
3340
2640
140
480
2620 340 140 140 340
1900 280 140
3440
140
2240
260
3140 360 140
140
2320
3660
340
140 3220
860
3860
1520
260
500
1480 260 140 2500
340 140 140 3400
140 140
4700
340
700
140
140
340
3120
440
280
140
140
280
2460
140
140
140
580
3680
140
140
140
260
1480
140
140
340
140
280
2480
800
340
2620
800
1080 4460
4220
4380
4
1060 4800
Figure 2.1: A 49-node Manhattan Grid Dynamic recon guration [3] has been proposed as a simple and robust resource management control to manage ATM networks to respond optimally to slow variations in trac call patterns. The network is logically con gured as a fully-meshed VPCN. In the next section we present a model of dynamic recon guration and a deterministic algorithm (XFG) [2, 3] which can solve this model. Section 2.3 describes the method of simulated annealing, its physical analogy and its modelling by means of a Markov chain. The type of optimization problem this method is intended for is also discussed. In section 2.4 the simulated annealing method is used to solve the dynamic recon guration model. The results obtained with simulated annealing and XFG are discussed and compared in section 2.5.
2.2 Dynamic Recon guration Consider a network that carries S services. Let bs denote the e ective bandwidth of a call of class s (1 s S ). For each origin-destination (OD) pair (i; j ), let Cij denote the capacity (zero if the physical link does not exist) of the unidirectional link i ? j , sij the Poisson arrival rate and 1=sij the mean holding time of class-s calls o ered to OD pair (i; j ). The arrival intensity for calls of class s o ered to OD pair (i; j ) is sij = sij =sij . The rate at which a class-s call o ered to OD pair (i; j ) generates revenue does not depend on how the call is routed (this assumption can be relaxed if cross-connection costs are taken into account) and is given by ijs . The e ective bandwidth of a call of class s is also independent of the call's route. Dynamic recon guration is performed by cross-connecting circuits (using
2.2 Dynamic Recon guration
5
con gurable switches in transit nodes) to form logical direct links (VPCs) between all OD pairs. Let Rij denote the set of admissible routes (including the direct link if it exists) between OD pair (i; j ). A route r = (i; i1; : : : ; ik ; j ) is a sequence of physical links i ? i1 ; i1 ? i2 ; : : : ; ik ? j connecting nodes i and j . A route r between OD pair (i; j ) is admissible if circuits can be reserved on each link along r for the sole use of calls between i and j . A VPC connecting (i; j ) consists of a set Rij of routes where route r 2 Rij has a capacity of xr bandwidth units. A VPC is not reserved for a speci c trac class (service seperation) { all service classes share a VPC (service integration). Let R be the set of all admissible routes in the network and let x = (x1 : : : xR ) where R = jRj. Let Aij denote the set of admissible routes that use the physical link i ? j . The (logical) capacity of the VPC connecting OD pair (i; j ) is given by
Cbij (x) =
X r2R
xr
ij
Let Bijs denote the blocking probability of class-s calls on the VPC connecting OD pair (i; j ). Bijs is a function of the logical capacity Cbij (x), the trac intensities sij and the e ective bandwidths bs . The VPCN design problem can now be speci ed in terms of the following non-linear optimization problem:
Optimization Problem 2.1. Maximize:
F (x) = Subject to:
X r2A
S XX ij s=1
ijs sij (1 ? Bijs )
xr Cij
xr 0
ij
Berezner and Krzesinski [3] developed an ecient deterministic algorithm (XFG) that solves optimization problem 2.1. XFG computes an optimal VPCN by performing a sequence of connections and releases. During a connection circuits are cross-connected along route r 2 Rij and during a release circuits (previously cross-connected along route r) are released and returned to each link i ? j along the route r. XFG selects the route r so that a cross-connection (or release) leads to the largest increase in revenue; XFG terminates when no such route r can be found. Since at each step XFG obtains the largest possible increase in revenue and does not take long-term impacts of the connect/release performed at each step into account, it is a so-called greedy algorithm . Despite being a one-step greedy algorithm, XFG computes a global optimum for all the networks tried thus far. The XFG algorithm is discussed in detail in [2].
6
Chapter 2. Computing an Optimal VPCN by Simulated Annealing
2.3 Simulated Annealing
2.3.1 Combinatorial Optimization
Before introducing the method of simulated annealing, we brie y discuss the type of problem this technique is applicable to. Lawler [11] de nes the related elds of discrete optimization , integer programming and combinatorial optimization as follows:
Discrete optimization encompasses all types of optimization problems for
which discrete solutions are sought. Integer programming and combinatorial optimization are special types of discrete optimization, with contrasting goals and methodologies.
Integer programming deals with mathematical programming problems in
which at least some of the variables are required to be integers. Integer programming algorithms generally exploit only very weak assumptions about the characteristics of problem instances.
Combinatorial optimization deals with problems with signi cant combinatorial structure, often of a type that is awkward or inappropriate to formulate in a standard form such as Ax b. Combinatorial optimization algorithms are designed to exploit the speci c structure of problem instances.
The well-known travelling salesman problem is an example of a combinatorial optimization problem. A salesperson has to visit a speci ed number of cities and return to his or her hometown. Each city has to be visited exactly once and the length of the whole journey has to be minimized. This problem is in general NP-complete but instances of the problem have been solved by simulated annealing. According to Aarts and Korst [1] one of the largest instances of this problem which has been solved is a 442-city instance solved by Grotschel. According to Aarts and Korst [1] simulated annealing (SA) is primarily intended for solving combinatorial optimization problems. However, as shown by Press et al [16], SA can also be applied to continuous optimization problems. The problem addressed in this chapter is combinatorial.
2.3.2 Simulated Annealing Theory
Physicists de ne annealing as a process for obtaining low energy states of a solid in a heat bath by slowly lowering the temperature of the heat bath. If a solid or liquid is slowly cooled and equilibrium is allowed to set in at each temperature, the atoms constituting the object will align themselves in the position of minimum energy. The obverse of annealing is quenching where the temperature is rapidly lowered without reaching thermal equilibrium. Whereas the result of annealing is a perfect crystalline lattice, quenching leads to a polycrystalline or amorphous state with an energy higher than the minimum reached by annealing. The probability Pr (T ) that a system in equilibrium in a heat bath at temperature T will be in a state r with energy Er , is given by the Boltzmann distribution (see for example Mandl [14]):
2.3 Simulated Annealing
7
Pr (T ) = Z (1T ) e?E =kT X ?E =kT Z (T ) = e r
r
r
where k is Boltzmann's constant and Z (T ) is a normalization constant known as the partition function . As reported in [1], the Metropolis algorithm for simulating the evolution of a solid in a heat bath to thermal equilibrium was introduced in 1953. This algorithm generates a sequence of system states by applying a perturbation mechanism and an acceptance criterion (the Metropolis criterion ). If the system is in state i with energy Ei the perturbation mechanism generates a proposed new state k with energy Ek from state i. The Metropolis criterion gives the probability Pik (T ) of accepting the new state. The Metropolis criterion is derived from the Boltzmann distribution and given by
(
? ?E ifif EEk > EEi exp E kT k i Thus if the proposed new state has the same energy or a lower energy than the current state, a state transition always takes places. However, if the proposed state is at a higher energy level, a state transition takes places with probability Pik (T ). In SA the Metropolis algorithm is applied to generate a sequence of solutions of a combinatorial optimization problem. Consider a combinatorial optimization problem where a function f has to be optimized. It can be assumed, without loss of generality, that f has to be minimized (if a function g has to be maximized then set f = ?g). Let S be the set of all possible solutions (all vectors n which satisfy the constraints). In the travelling salesman problem n would correspond to a permutation of the cities and f would be the total length of the salesperson's tour. f is clearly a function of the vector n. Although the cardinality of S (the state space) can be very large, it is nite. The main equivalences [1] between this combinatorial optimization problem and the annealing process are Pik (T ) = 1
i
k
the solutions of the combinatorial optimization problem correspond to the states of the system in the heat bath the cost of the solution (the value of f ) corresponds to the energy of the heat bath.
The analogue of temperature, in the simulated annealing process, is the control parameter c. The Metropolis criterion modi ed for simulated annealing becomes:
(1
if f (m) f (n) Pnm(c) = exp f (n)?f (m) if f (m) > f (n) c
8
Chapter 2. Computing an Optimal VPCN by Simulated Annealing where n is the current solution and m is the proposed new solution. The state transitions from n to m correspond to transitions from one solution to another. A \downhill move" (where f (m) f (n)) which lowers the energy of the system is always accepted and an \uphill move" (where f (m) > f (n)) which increases the energy of the system is accepted with probability Pnm (c). The reason for accepting some \uphill moves" is to avoid inferior local minima. Some optimization algorithms correspond to quenching rather than annealing: they rapidly move \downhill" and frequently terminate in inferior local minima (\poly-crystalline structures" or \amorphous states"). In [1] it is proved that if the Metropolis criterion is used in the above form and a suciently large number of transitions are allowed at a xed value c of the control parameter, the probability that the SA algorithm nds a solution n 2 S is given by the stationary distribution
(2.1) qn (c) = N1(c) exp ? f (cn) where the normalization constant N (c) is analogous to the partition function Z (T ) in the Boltzmann distribution. Aarts and Korst [1] also prove that if the stationary distribution is attained at each value of c, the SA algorithm will
asymptotically converge to the set of globally optimal solutions. The SA algorithm can be modelled by a Markov chain. A homogeneous Markov chain is associated with the transitions at each value of the control parameter. The stationary distribution (2.1) is obtained at a speci c value of c if the homogeneous Markov chain associated with this value satis es the socalled detailed balance equations . This requires in nitely long Markov chains [1]. In practice these chains are approximated by Markov chains of nite length. Thus convergence to a global optimum is no longer guaranteed, although the probability of obtaining the global optimum should be large if a \good" approximation is used. We therefore only obtain quasi-equilibrium at each value of the control parameter. The nite-length Markov chain approximation is implemented by the cooling schedule. The cooling schedule has the following functions [1]: it speci es initial and nal values for the control parameter c it determines how much the control parameter is decreased each time quasi-equilibrium is reached by specifying a decrement function it speci es the number of transitions at each value of the control parameter which determines the length of the Markov chain. The analogy with thermodynamics can be extended further by de ning the entropy of the simulated annealing process. The entropy of an isolated system is de ned [14] as:
S = ?k
X r
Pr (T ) ln Pr (T )
where k is Boltzmann's constant and Pr (T ) is the probability of being in state r when the system is at temperature T . The entropy for the simulated annealing process is now de ned as
2.4 The Simulated Annealing Implementation Sc = ?
X n2S
qn (c) ln qn (c)
As in thermodynamics the entropy can be interpreted as the extent of disorder in the system. The third law of thermodynamics, which states that the entropy tends to zero as the absolute temperature approaches zero, hinges on the conjecture that the energy eigenvalues of the ground state of any system are non-degenerate. In combinatorial optimization problems where f has a unique global minimum (the state n0 2 S ) the third law holds with our de nition of entropy as well:
S0 = ?qn0 (0) ln qn0 (0) = ? ln 1 = 0 Aarts and Korst [1] show that the entropy Sc decreases monotonically during the execution of the SA algorithm { provided that equilibrium is reached at each value of the control parameter.
2.4 The Simulated Annealing Implementation The main reason for using the XFG algorithm is the computational complexity of standard non-linear optimization algorithms when applied to optimally recon gure large networks. The non-linear optimizer CFSQP (Lawrence [12]), for example, requires several hours to solve optimization problem 2.1 for a model of the NSF ATM backbone (Mitra [15]). XFG solves the same model in a few seconds. Both execution times are for a Pentium 100 MHz machine. Although XFG considers all routes in a network, it identi es and works with a small number of active routes . Standard non-linear optimizers require a variable for each possible route in the network. Since a J -node network has in the order of J ! possible routes and few non-linear solvers can handle more than a thousand variables, these optimizers cannot solve optimization problem 2.1 for large networks, whereas XFG can. For veri cation purposes, however, we would like to compare the XFG result with a result computed using another method. We decided on simulated annealing. SA can eciently be build \on top of" XFG: we use XFG's way of identifying active routes, but maximize the revenue by means of SA rather than the XFG algorithm. In contrast with XFG, which always chooses the route r that yields the largest increase in revenue, our SA implementation chooses a route at random. The Metropolis criterion is then applied in order to decide whether or not to perform the proposed action (connection or release) on this route. The SA implementation is such that only revenue-increasing routes are selected for connects but revenue-decreasing routes may be selected for releases. The SA method computes the same optimal revenue as XFG for all the networks tried thus far. However, the same cooling schedule could not be used for all networks.
9
10
Chapter 2. Computing an Optimal VPCN by Simulated Annealing
2.5 Results The 49-node manhattan grid shown in gure 1 presents a network model which is too large to be solved by standard non-linear optimizer. The grid points represent nodes and the lines represent physical links. Each link has a transmission capacity of 10; 000 bandwidth units. Figure 1 shows the transmission capacity reserved on each link for carrying direct trac after recon guration. Clearly, most of the transmission capacity has been allocated to routes forming VPCs. We solve it by SA using the following decrement function (Press et al [16]):
c = c0 (1 ? k=K ) where K is the total number of steps (connections/releases) the algorithm takes to complete. K has to be decided upon before starting the algorithm. Let k denote the number of the current step and c0 denote the initial value of the control parameter. We modify this somewhat and let K be the total number of connections and k the number of connections already performed at this point. The cooling schedule is given by algorithm 2.1. The control parameter c is decreased when L! (where L is the number of physical links in the network) moves have been approved or L moves have been attempted, whichever occurs rst. This determines the number of transitions made at each value of the control parameter, which in turn in uences the lengths of the Markov chains and whether or not a stationary distribution is reached. The 49-node manhattan grid can be recon gured by setting c0 = 1, = 100, ! = 10 and K = 500000. Determining the appropriate value of requires many experiments. [16] provides the following possibilities: = 1; 2; 4. Algorithm 2.1 (Cooling Schedule). Let L be the number of physical links, c0 the initial value of the control parameter, R 2 (0; 1) a random number and K the total number of connects. 1. Set i 0; j 0; k 0; c c0 2. Consider a release or a connect with revenue increment r 3. Metropolis Algorithm Generate R and do one of the following: If r > 0 execute release or connect Set i i + 1; j j + 1 If it is a connect operation, set k k + 1 If R < e =c (\downhill move") execute release or connect Set i i + 1, j j + 1 If it is a connect operation, set k k + 1 If none of the above holds, do not perform the release or connect Set i i + 1 4. If i = L or j = L! Set c c0 (1 ? k=K ) Set i 0; j 0 r
2.6 Conclusion
11
Method XFG SA Revenue 45683:6 45618:9 Blocking Probability 0:0525215 0:0538636 Execution Time 10 min. 5 hours Table 2.1: Results for the Manhattan Grid Operation Connects Releases Revenue-increasing 144890 7817 Revenue-decreasing 0 59801 Total 144890 67618 Attempted \downhill moves" 0 346401 Table 2.2: Simulated Annealing Moves for the Manhattan Grid 5. If the algorithm has converged or k = K , stop. Otherwise go to step 2. During many trials (with di erent initial values of c given by c0 where 1 c0 100) it was observed that = 1 or 2 caused the algorithm to make a large number ( 105) of revenue-decreasing releases which led to inferior local maxima. Although \downhill moves" are intended to avoid inferior local maxima (to \backtrack" from inferior local maxima), too many of them seem to \backtrack" the algorithm from the global maximum. = 4 caused the value of the control parameter to drop much more quickly which in turn led to fewer \downhill moves" being accepted. The results obtained with simulated annealing and XFG for this model are compared in table 2.1. XFG solves this model in about 10 minutes on a Pentuim-Pro 200 MHz machine. The SA method requires more than ve hours on the same machine. The total number of connections and releases made by the SA method are shown in table 2.2. By design, no revenue-decreasing connections are attempted. On the other hand 88% of all releases led to a decrease in revenue. However, only 17% of the attempted \downhill" releases were accepted.
2.6 Conclusion We have used simulated annealing to con rm the results computed by XFG for a large network model. The results were the same but XFG completed the calculation in 10 minutes versus the more than ve hours required by simulated annealing.
Chapter 3
The Reduced Load Approximation Ecient algorithms for calculating exact blocking probabilities (and other performance measures) do not exist for all product-form loss networks (Ross [18]). The reduced load approximation is a method of approximating blocking probabilities in product-form loss networks. In x3.1 and x3.2 we present a model of a product-form loss network and the method of reduced load approximation. Our implementation is discussed in x3.3. Two techniques for dealing with oscillatory behaviour of the repeated substitution algorithm (introduced in x3.3) are surveyed in x3.4. x3.6 is devoted to numerical examples and some observations.
3.1 A Product-form Loss Network Model Consider a product-form loss network (de ned in chapter 1) with L (physical) links which carries K trac classes. Let K = f1; : : :; K g denote the set of trac classes and L the set of links. Class-k calls are o ered to O-D pair (i; j ) according to a Poisson process with rate kij . The average holding time of such a call is 1=kij , yielding an o ered load of ijk = kij =kij . The bandwidth requirement of a class-k call is bk 1 . We de ne a route r = (i; i1; : : : ; ik ; j ), connecting nodes i and j , as a sequence of links i ? i1 ; i1 ? i2 ; : : : ; ik ? j . Figure 3.1 contains an example of a route connecting nodes i and j .
_ i
_ _ _ _ _ _ _ _ _ _ _ _ _
i1
_
_ _ _ _ _ _ _ _ _ _ _ _ _
ik?1
i2
ik
j
Figure 3.1: A Route Connecting Nodes i and j Let Rij denote the set of routes that connect i and j , including the direct link i ? j if it exists (the network is not necessarily fully-meshed). Let R be the 1 Note that the bandwidth requirements are only class-dependent and not link-dependent. In chapter 4 we allow the bandwidth requirements to be link-dependent as well.
12
3.2 The Method
13
set of all routes2 in the network:
R=
[ (i;j)
Rij :
Let kr denote the class-k load o ered to route r. The load o ered to O-D pair (i; j ) is shared equally between all the routes connecting O-D pair (i; j ):
kr = jR1 j ijk for all k 2 K and r 2 Rij . ij
A class-k call connected between O-D pair (i; j ) generates revenue at a rate ijk . For link i ? j , let Aij denote the set of routes that use the link and Cij the capacity (number of bandwidth units) of the link.
3.2 The Method
Our goal is to approximate the link blocking probabilities Bijk for all trac classes k 2 K on all links i ? j 2 L in the network. Consider the routes r1 , r2 and r3 in gure 3.2. r1
J _ _ t _ Jt _ Jt _ _ _ _ jJt _ tJ _ tJ _ _ r2 i r3
Figure 3.2: Three Routes in a Loss Network Let ekij denote the class-k load we attempt to carry on link i ? j :
ekij =
X
r2A
kr :
ij
In calculating the blocking probability for class-k calls on link i ? j we want to take the blocking probabilities on the other links of the routes in Aij (r1 , r2 and r3 in the case of gure 3.2) into account. Therefore we cannot simply use ekij as the class-k load on link i ? j . We make use of the link independence assumption { we assume that the link blocking probabilities on the di erent links are independent. The blocking probability of class-k trac on route r is now given by
( Q 1?
k k i?j 2r (1 ? Bij ) if r > 0; (3.1) 0 if kr = 0: 2 These routes have to be known before applying the reduced load approximation. One option is using least-cost routes supplied by Dijkstra's algorithm (see Manber [13] or Tucker [19] for example).
Lkr =
14
Chapter 3. The Reduced Load Approximation Let bkij be the class-k load o ered to link i ? j when blocking is taken into account:
bkij = =
X r2A
ij
X
r2A
Y
kr
k) (1 ? Bod
o?d2r o?d6=i?j (1 ? Bijk )?1 kr (1 ? Lkr ):
(3.2)
ij
We calculate the link blocking probabilities Bijk using the reduced loads bkij : (3.3) Bijk = Eijk (bij ; b; Cij ) where bij = (bkij j k 2 K), b = (bk j k 2 K) and Eijk is a function returning the blocking probability of class-k calls on link i ? j .
Equations 3.1, 3.2 and 3.3 are clearly interdependent and have to be solved as a system of non-linear equations.
Bijk = Eijk (bij ; b; Cij ) where k 2 K and i ? j 2 L Y (1 ? Bijk ) where k 2 K and r 2 R 1 ? Lkr = bkij =
i?j 2r
X
r2A
(1 ? Bijk )?1 kr (1 ? Lkr ) where k 2 K and i ? j 2 L.
ij
The above equations are known as the Erlang xed-point equations . The approximation (3.2) which led to these equations is known as the reduced load approximation { the load on each physical link is approximated by taking the blocking probabilities on the other links into account.
3.3 Implementation
In this section we discuss our implementation of the method of x3.2. We consider our choice of blocking function, additional performance measures we calculate and methods of solving the Erlang xed-point equations.
3.3.1 The Blocking Function Eijk
There are several options for the blocking function Eijk . We use two, one being the Berezner-Krzesinski stochastic knapsack algorithm (see appendix C) and the other the multi-rate extension of the Erlang-B function given by the following equation:
"
Eijk (bij ; b; Cij ) = 1 ? 1 ? E
K X k=1
bk bkij ; Cij
!#b
k
where E (; C ) is the Erlang-B function , evaluated as described in appendix B.
3.3 Implementation
15
3.3.2 Additional Performance Measures
An additional performance measure we calculate is the network blocking probability per class:
P eBk = Pr2R kr kLkr P (i;j)? kij=jR j P Lk ij r2R r : = (i;j) ijP k ij
(i;j) ij For comparison purposes we introduce the following function for the network's rate of earning revenue:
W=
XX
ijk
X
kr (1 ? Lkr )
r2R XX k k 1 X (1 ? Lkr ): = ij ij jR j ij r2R (i;j) k2K
(i;j) k2K
(3.4)
ij
ij
3.3.3 Solving the Erlang Fixed-point Equations
We employ three di erent methods of solving the Erlang xed-point equations: multi-dimensional Newton's method Broyden's method repeated substitution. Newton's and Broyden's methods are algorithms for solving systems of nonlinear equations (see Press et al [16]). Repeated substitution is a standard method of solving the Erlang xed-point equations (see Ross [18]). We now discuss these two approaches (non-linear equation solvers and repeated substitution).
Non-linear Equation Solvers
The solution of the Erlang xed-point equations is the set of link blocking probabilities Bijk satisfying
Bijk = Eijk (bij ; b; Cij ) where k 2 K and i ? j 2 L Y (1 ? Bijk ) where k 2 K and r 2 R 1 ? Lkr = bkij =
i?j 2r
X
r2A
(1 ? Bijk )?1 kr (1 ? Lkr ) where k 2 K and i ? j 2 L.
ij
These correspond to the following system of non-linear equations:
Lij;k ? Eijk (bij ; b; Cij ) = 0 for all k 2 K and i ? j 2 L:
16
Chapter 3. The Reduced Load Approximation The solution L = (Lij;k j for all i ? j 2 L and k 2 K) of the above equations provides the required link blocking probabilities Bijk . Newton's and Broyden's methods proceed as explained in appendix D; both require an initial point L0 = (L0ij;k j i ? j 2 L and k 2 K) (initial approximations of the link blocking probabilities). We calculate the initial point as follows:
L0ij;k = Eijk (eij ; b; Cij )
(3.5)
where eij = (ekij j k 2 K). The non-linear equation solver generates a sequence of points L1 ; : : : ; Lm of which the last point is an approximate solution to the set of equations. The functions f in (D.1) are de ned as:
fijk = Lij;k ? Eijk (bij ; b; Cij ) for all k 2 K and i ? j 2 L:
For each point Ln we have to calculate the values of all these functions: 1 ? Lkr =
bkij =
Y i?j 2r
X
(1 ? Lnij;k ) where k 2 K and r 2 R (1 ? Lnij;k )?1 kr (1 ? Lkr ) where k 2 K and i ? j 2 L
r2A k fij = Lnij;k ? Eijk (bij ; b; Cij ) for all k 2 K and i ? j 2 L: ij
After the non-linear equation solver has terminated, we set
Bijk = Lmij;k for all k 2 K and i ? j 2 L.
Repeated Substitution
Solving the Erlang xed-point equations, using the method of repeated substitution, involves starting with an initial point L0 = (L0ij;k j i ? j 2 L and k 2 K) (initial approximations of the link blocking probabilities) and generating a sequence of points L0 ; L1 ; : : : Lm ; : : : . We start with the same initial point (3.5) used by the non-linear equation +1 solver. Point Lm+1 = (Lm ij;k j i ? j 2 L and k 2 K) is obtained from point m m L = (Lij;k j i ? j 2 L and k 2 K) in the following way: 1 ? Lkr =
bkij =
Y i?j 2r
X
(1 ? Lm ij;k ) where k 2 K and r 2 R ?1 k k (1 ? Lm ij;k ) r (1 ? Lr ) where k 2 K and i ? j 2 L
r2A m +1 Lij;k = Eijk (bij ; b; Cij ) where k 2 K and i ? j 2 L: ij
(3.6)
3.4 Oscillation of the Approximate Link Blocking Probabilities Lmij;k We stop as soon as
m+1 m Lij;k ? Lij;k < "RLA
for all k 2 K and i ? j 2 L
where "RLA is the convergence criterion and we set
Bijk = Lmij;k+1 for all k 2 K and i ? j 2 L.
3.4 Oscillation of the Approximate Link Blocking Probabilities Lmij;k
Under certain circumstances the approximate link blocking probabilities Lm ij;k in the repeated substitution algorithm oscillate between a high and a low value and never converge. We consider this problem by studying a numerical example. Figure 3.3 shows the only route r between O-D pair (1; 4): R14 = frg. Route r is also the only route using links 1 ? 2, 2 ? 3 and 3 ? 4: A12 = A23 = A34 = frg. The only load o ered to (1; 4) is 141 = 400 which is all o ered to route r: 1r = 400. The bandwidth requirement of class-1 trac is b1 = 1. Since this example only involves single-rate trac, both the multi-rate version of the Erlang-B formula as well as the stochastic knapsack reduce to the single-rate Erlang-B formula. C12 = 100 C23 = 50 C34 = 100
_ _ _ _ _ _ _ _ _ _
1
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2
3
4
Figure 3.3: A Route on which Oscillation Occurs We now apply the reduced load approximation to route r. The initial point
L0 consists of the following components:
1 (400; 1; 100) = 7:508260 10?1 L012;1 = E12 1 (400; 1; 50) = 8:753548 10?1 L023;1 = E23 1 (400; 1; 100) = 7:508260 10?1 L034;1 = E34 which are obtained from (3.5) by noting that eij = (e1ij = 1r ) and b1 = 1.
L1 is now obtained from L0 in the following way: 1 ? L1r =
Y i?j 2r
(1 ? L0ij;1 )
= 7:738931 10?3
b112 = (1 ? L012;1 )?1 1r (1 ? L1r ) = 12:423332
17
18
Chapter 3. The Reduced Load Approximation b123 = (1 ? L023;1 )?1 1r (1 ? L1r ) = 24:835065 b134 = (1 ? L034;1 )?1 1r (1 ? L1r ) = 12:423332 1 (12:423332; 1; 100) = 1:143985 10?54 L112;1 = E12 1 (24:835065; 1; 50) = 3:051124 10?6 L123;1 = E23 1 (12:423332; 1; 100) = 1:143985 10?54 L134;1 = E34
In the following iteration L2 is obtained from L1 : 1 ? L1r =
Y i?j 2r
(1 ? L1ij;1 )
= 9:999969 10?1
b112 = (1 ? L112;1)?1 1r (1 ? L1r ) = 399:998780 b123 = (1 ? L123;1)?1 1r (1 ? L1r ) = 400:000000 b134 = (1 ? L134;1)?1 1r (1 ? L1r ) = 399:998780 1 (399:998780; 1; 100) = 7:508253 10?1 L212;1 = E12 1 (400:000000; 1; 50) = 8:753548 10?1 L223;1 = E23 1 (399:998780; 1; 100) = 7:508253 10?1 L234;1 = E34
Continuing in this fashion, we can draw up table 3.1 for m = 0; : : : ; 5.
m 0 1 2 3 4 5
Lm12;1 = Lm34;1 7:508260 10?1 1:143985 10?54 7:508253 10?1 1:144290 10?54 7:508253 10?1 1:144290 10?54
Lm23;1
8:753548 10?1 3:051124 10?6 8:753548 10?1 3:051590 10?6 8:753548 10?1 3:051590 10?6
Table 3.1: Lm ij;1 It follows that L223;1 = L023;1 and L434;1 = L412;1 = L212;1 = L234;1. Figure 3.4 plots the values of Lm ij;1 for the rst ve iterations. From gure 3.4 and the above calculations it follows that Lm ij;1 oscillates for all links i ? j 2 r. The oscillation of Lm is detected when ij;k
m m?2 Lij;k ? Lij;k < "oscillation
3.4 Oscillation of the Approximate Link Blocking Probabilities Lmij;k where "oscillation is the detection criterion. With "oscillation = 10?10 the oscillation of Lm 23;1 is detected at m = 2 and the oscillation of Lm 12;1 is detected3 at 4 m=3 . 1 i?j =1?2=3?4 i?j =2?3 0:8
Lmij;1
0:6 0:4 0:2 0
0
1
2
m
3
4
5
Figure 3.4: Oscillation of Lm ij;1 We use two ways of dealing with this problem. The rst is to restart the repeated substitution algorithm with new initial values L0ij;k calculated as follows:
L0ij;k =
( jL
M ij;k
?L ?1 j if Lm ij;k is oscillating; 2 M ij;k
LM if Lm ij;k ij;k is not oscillating where M is the value of m when oscillation is detected. In the current example M = 2. Figure 3.5 shows the e ect of this method on the current route. The
repeated substitution algorithm is restarted four times (corresponding to the gaps in the graph). Another way of dealing with oscillation is to apply damping5 to Lm ij;k by setting
and calculating Le
Le0ij;k = L0ij;k for all k 2 K and i ? j 2 L
m+1
from Le in the following way: m
3 Since route r is symmetric the behaviour of Lm and Lm are identical. Therefore all 12;1 34;1 m statements regarding Lm 12;1 apply to L34;1 as well. 4 If our algorithm is instructed to test for oscillation and react to it, it will react as soon as oscillation is detected on a link. In our current example the algorithm will take action at m m = 2 (when Lm 23;1 starts to oscillate) and not wait for L12;1 to start oscillating at m = 3. 5 Dr Peter Taylor from the University of Adelaide in Australia suggested this method.
19
20
Chapter 3. The Reduced Load Approximation
Lmij;1
0:9 0:8 0:7 0:6 0:5 0:4 0:3 0:2 0:1 0
i?j = 1?2 i?j = 2?3
0
20
30
20
m
60
5
Figure 3.5: The E ect of Restarting 1 ? Lkr =
bkij =
Y i?j 2r
X
(1 ? Le m ij;k ) where k 2 K and r 2 R ?1 k k (1 ? Lem ij;k ) r (1 ? Lr ) where k 2 K and i ? j 2 L
r2A m +1 Lij;k = Eijk (bij ; b; Cij ) where k 2 K and i ? j 2 L Le mij;k+1 = (1 ? )Lmij;k+1 + Lemij;k where k 2 K and i ? j 2 L. ij
We stop as soon as
m+1 m Leij;k ? Leij;k < "RLA
for all k 2 K and i ? j 2 L
and again set
Bijk = Le mij;k+1 for all k 2 K and i ? j 2 L. We refer to (0 < 1) as the damping parameter. Figure 3.6 shows the e ect of damping on Lm 23;1. The graph with = 0 corresponds to the undamped
case. Note that restarting is a way of dealing with oscillation when it occurs. Damping is a way of preventing oscillation from occurring. The behaviour of Lm ij;k is sometimes erroneously categorized as oscillatory. We provide a simple example. Consider the links shown in gure 3.7. The only route r connecting O-D pair (1; 3) consists of the links 1 ? 2 and 2 ? 3. Route r is also the only route using the links 1 ? 2 and 2 ? 3.
3.4 Oscillation of the Approximate Link Blocking Probabilities Lmij;k 1
= 0:5 =0
0:8 0:6
Le m
23;1
0:4 0:2 0
01
5
10
m
15
20
23 25
Figure 3.6: The E ect of Damping 1 = 400 and the bandwidth The only load o ered to O-D pair (1; 3) is 13 requirement of class-1 calls is b1 = 1. Since r is the only route in R13 , 1r = 400. We now apply the reduced load approximation to route r. The results are displayed in gure 3.8.
C12 = 100
C23 = 50
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
1
2
3
Figure 3.7: A Route on which Oscillation Seems to Occurs Since L223;1 = L023;1 = 8:753548 10?1 oscillation is detected on link 2 ? 3. However, if we disable oscillation detection (by setting "oscillation = 0) and do not apply damping, the repeated substitution algorithm proceeds beyond m = 2 and Lmij;1 converge for all links (see gure 3.9). If we apply damping the repeated substitution algorithm requires 23 iterations to converge. When damping is not applied but oscillation detection is enabled, oscillation of Lm 23;1 is detected at m = 2. The algorithm is restarted as explained above and requires another 8 iterations to converge. Whichever approach is taken in this example, convergence is attained, however damping and restarting delay convergence. We make one nal observation regarding oscillation. If the convergence criterion has been met for all trac classes on all links, Lm ij;k can still oscillate for some values of i ? j and k but our algorithm will regard it as convergence. However, if the convergence criterion has not been met for some trac classes on some links and there is oscillation between two values Lm+2 and Lm , such ij;k
ij;k
21
22
Chapter 3. The Reduced Load Approximation
Lmij;1
0:9 0:8 0:7 0:6 0:5 0:4 0:3 0:2 0:1 0
i?j = 1?2 i?j = 2?3
0
1
2
m
3
4
Figure 3.8: Seeming Oscillation of Lm ij;1
+2 m+1 that Lm ij;k ? Lij;k < "RLA (for some trac classes k on some links i ? j ) the user can instruct our algorithm to take note of this oscillation or ignore it.
3.5 Computer Implementation
Our computer implementation of the reduced load approximation is called RLA. The model discussed in chapter 4 also makes use of the reduced load approximation. Consequently the data structures and algorithms used in implementing the reduced load approximation are subsets of those used in implementing the model of chapter 4. Therefore we refer the reader to x4.5. The most important data structure used in RLA is the route structure in gure 3.10.
3.6 Numerical Examples and Observations 3.6.1 A Simple Example
Our rst example was inspired by an example due to Farago [7]. The network consists of three nodes and three uni-directional links and is shown in gure 3.11. The network carries two trac classes K = f1; 2g with bandwidth requirements b1 = 1 and b2 = 2. The rate of earning revenue ijk = 1 for all O-D pairs and trac classes. Four routes (shown in gure 3.12) are used to carry the o ered loads which are given by table 3.2. Table 3.3 contains the revenues obtained using the di erent methods of solving the Erlang xed-point equations and the two blocking functions. Tables 3.4 to 3.6 contain various performance measures calculated for the three-node network. In all cases the Erlang xed-point equations were solved
3.6 Numerical Examples and Observations
Lmij;1
0:9 0:8 0:7 0:6 0:5 0:4 0:3 0:2 0:1 0
23
i?j =1?2 i?j =2?3
0
1
2
3
m
4
5
6
7
Figure 3.9: Convergence of Lm ij;1 (i; j ) k = 1 k = 2 (1; 2) 40 80 (1; 3) 80 0 (3; 2) 80 0 Table 3.2: The O ered Loads ijk for the 3-node Network using the method of repeated we did not msubstitution. As explained in x3.4 m ? 2 ?1 m m regard Lij;k as oscillating if Lij;k ? Lij;k < "oscillation but Lij;k ? Lm ij;k < "RLA .
3.6.2 A Larger Network We now consider the ve-node network in gure 3.13 which carries four trac classes K = f1; 2; 3; 4g. Each bi-directional link i ? j has a capacity of Cij = 100 bandwidth units. The revenue factors are ijk = k for all O-D pairs (i; j ) and all trac classes k. The bandwidth requirements are bk = 2k?1 for all trac classes k. The o ered loads are given in table 3.7. There is one route for each Newton Broyden Repeated Substitution Stochastic Knapsack 199:087784 199:087026 199:086971 Erlang-B 188:759184 188:759732 188:759240 Table 3.3: The Revenues for the 3-node Network
24
Chapter 3. The Reduced Load Approximation /* the route struct */ struct route { int O; int D; int length; double *loads; int *links; int *SortedLinks; double *Blocking; };
Figure 3.10: RLA's route Structure ' !& 1"% =#$
C13
C12 = 100
== == == = 100 = = = = == ==
C32 ' !& 3"% #$
' !& 2"% #$
= 100
Figure 3.11: A Simple 3-node Network O-D pair given in table 3.8. We now apply the reduced load approximation to this network, solving the Erlang xed-point equations using the method of repeated substitution. Table 3.9 contains the link blocking probabilities, table 3.10 the route blocking probabilities and table 3.11 the network blocking probability per service as well as the total rate of earning revenue (3.4). Using the stochastic knapsack, the repeated substitution algorithm converged without damping { no oscillation occurred. For the multi-rate Erlang-B formula oscillation was erroneously detected. However, convergence was obtained by disabling oscillation detection. In both cases the convergence criterion was "RLA = 10?6. It is evident from table 3.9 that the link blocking probabilities given by the multi-rate Erlang-B formula are considerably higher than those given by the stochastic knapsack. These in turn in uence the route blocking probabilities (table 3.10) and the network blocking probabilities (table 3.11). The network's route 1 /' !& "% #$ == @2 == == route route 2 = = = = == == ' !& 3%" #$ ' !& 1"% =#$
' !& 1"% =#$
3
== == == route 4 = = = = == ==
' !& "% #$ @2 route ' !& 3"% #$
Figure 3.12: The 4 Routes in the 3-node Network
4
3.6 Numerical Examples and Observations
25
Stochastic Knapsack Erlang-B k=2 k=1 k=2 1 ? 2 5:663644 10?2 1:134130 10?1 7:570045 10?2 1:456703 10?1 1 ? 3 2:459570 10?1 4:356925 10?1 2:808755 10?1 4:828600 10?1 3 ? 2 2:459570 10?1 4:356925 10?1 2:808755 10?1 4:828600 10?1
i?j
k=1
Table 3.4: The Link Blocking Probabilities Bijk for the 3-node Network r 1 2 3 4
Stochastic Knapsack k=2 5:663644 10?2 1:134130 10?1 2:459570 10?1 0 2:459570 10?1 0 4:314191 10?1 6:815571 10?1
k=1
k=1
Erlang-B
k=2 7:570045 10?2 1:456703 10?1 2:808755 10?1 0 2:808755 10?1 0 4:828600 10?1 7:325662 10?1
Table 3.5: The Route Blocking Probabilities Lkr for the 3-node Network rate of earning revenue given by the multi-rate Erlang-B calculation is only 68:3% that of the revenue rate given by the stochastic knapsack calculation. One of these methods is de nitely substantially more accurate than the other when applied to this network.
' !& 1"% ?#$
' !& 2"% #$
?? ?? ?? ?? ?? ?? ? ' !& "% #$ 5
' !& 4"% #$
' !& 3"% #$
Figure 3.13: The 5-node Network
26
Chapter 3. The Reduced Load Approximation
k=1 k=2
Stochastic Knapsack Erlang-B 2:455711 10?1 2:805565 10?1 3:974850 10?1 4:391183 10?1
Table 3.6: The Network Blocking Probabilities Bek for the 3-node Network
(i; j ) k = 1 k = 2 k = 3 k = 4 (1; 2) 20 5 10 5 (1; 3) 10 0 0 5 (1; 4) 5 5 5 5 (1; 5) 20 3 5 3 (2; 3) 10 5 5 5 (2; 4) 20 4 3 10 (2; 5) 20 10 5 5 (3; 4) 80 4 3 5 (3; 5) 10 1 5 1 (4; 5) 10 20 8 1 Table 3.7: The O ered Loads ijk for the 5-node Network
r
1 2 3 4 5 6 7 8 9 10
(i; j ) (1; 2) (1; 3) (1; 4) (1; 5) (2; 3) (2; 4) (2; 5) (3; 4) (3; 5) (4; 5)
links 1?2 1 ? 5; 5 ? 2; 2 ? 3 1?4 1?5 2?3 2?4 2?5 3?4 3 ? 2; 2 ? 5 4 ? 1; 1 ? 5
Table 3.8: The Routes in the 5-node Network
3.6 Numerical Examples and Observations (i ? j ) k Stochastic Knapsack 1?2 1 4:751983 10?2 2 9:433884 10?2 3 1:855257 10?1 4 3:558797 10?1 1?4 1 8:372419 10?2 2 1:623688 10?1 3 3:049231 10?1 4 5:353163 10?1 1?5 1 1:385401 10?1 2 2:603261 10?1 3 4:601690 10?1 4 7:243491 10?1 2?3 1 5:852979 10?2 2 1:153079 10?1 3 2:233421 10?1 4 4:159378 10?1 2?4 1 5:072040 10?2 2 1:001659 10?1 3 1:950413 10?1 4 3:676655 10?1 2?5 1 1:053742 10?1 2 2:017372 10?1 3 3:695390 10?1 4 6:197537 10?1 3?4 1 1:288713 10?1 2 2:441840 10?1 3 4:380988 10?1 4 7:054736 10?1
27 Erlang-B 1:361148 10?1 2:537024 10?1 4:430399 10?1 6:897955 10?1 1:638459 10?1 3:008463 10?1 5:111841 10?1 7:610591 10?1 2:606279 10?1 4:533289 10?1 7:011507 10?1 9:106891 10?1 6:811070 10?2 1:315823 10?1 2:458507 10?1 4:312589 10?1 1:962700 10?1 3:540181 10?1 5:827073 10?1 8:258668 10?1 2:994178 10?1 5:091846 10?1 7:591003 10?1 9:419673 10?1 3:012438 10?1 5:117397 10?1 7:616019 10?1 9:431663 10?1
Table 3.9: The Link Blocking Probabilities Bijk in the 5-node Network
r
k Stochastic Knapsack 1 4:751983 10?2 2 9:433884 10?2 3 1:855257 10?1 4 3:558797 10?1 2 1 2:744238 10?1 1
3
2 3 4 1 2 3 4
Erlang-B 1:361148 10?1 2:537024 10?1 4:430399 10?1 6:897955 10?1 5:172898 10?1 0 0 0 0 9:387814 10?1 9:970522 10?1 8:372419 10?2 1:638459 10?1 ? 1 1:623688 10 3:008463 10?1 3:049231 10?1 5:111841 10?1 5:353163 10?1 7:610591 10?1 continued on next page : : :
Table 3.10: The Route Blocking Probabilities Lkr in the 5-node Network
28
Chapter 3. The Reduced Load Approximation : : : continued from previous page r k Stochastic Knapsack Erlang-B 4 1 1:385401 10?1 2:606279 10?1 2 2:603261 10?1 4:533289 10?1 ? 1 3 4:601690 10 7:011507 10?1 ? 1 4 7:243491 10 9:106891 10?1 ? 2 5 1 5:852979 10 6:811070 10?2 ? 1 2 1:153079 10 1:315823 10?1 ? 1 3 2:233421 10 2:458507 10?1 ? 1 4 4:159378 10 4:312589 10?1 ? 2 6 1 5:072040 10 1:962700 10?1 ? 1 2 1:001659 10 3:540181 10?1 ? 1 3 1:950413 10 5:827073 10?1 ? 1 4 3:676655 10 8:258668 10?1 ? 1 7 1 1:053742 10 2:994178 10?1 ? 1 2 2:017372 10 5:091846 10?1 ? 1 3 3:695390 10 7:591003 10?1 ? 1 4 6:197537 10 9:419673 10?1 ? 1 8 1 1:288713 10 3:012438 10?1 ? 1 2 2:441840 10 5:117397 10?1 ? 1 3 4:380988 10 7:616019 10?1 ? 1 4 7:054736 10 9:431663 10?1 ? 1 9 1 1:577364 10 3:471350 10?1 ? 1 2 2:937832 10 5:737673 10?1 ? 1 3 5:103475 10 8:183256 10?1 ? 1 4 7:779125 10 9:669944 10?1 ? 1 10 1 2:106652 10 3:817710 10?1 ? 1 2 3:804261 10 6:177929 10?1 ? 1 3 6:247759 10 8:539177 10?1 ? 1 4 8:719095 10 9:786600 10?1 Table 3.10: The Route Blocking Probabilities Lkr in the 5-node Network
k
1 2 3 4 revenue rate
Stochastic Knapsack 1:199267 10?1 2:445286 10?1 3:692753 10?1 5:634497 10?1 437:834354
Erlang-B 2:727337 10?1 4:609674 10?1 6:218924 10?1 8:168419 10?1 299:089579
Table 3.11: The Network Blocking Probabilities Bek and Revenue Rate in the 5-node Network
Chapter 4
Route Optimization Using Virtual Subnetworks Farago et al [7] provide a mathematical model of the so-called logical con guration problem of ATM networks. This model is based on the assumption that services with widely di ering bandwidth demands, grades of service (GoS) or congestion control functions should not be completely integrated on the same physical network. Farago et al [7] propose handling these services in di erent logical networks on top of the same physical network and provide a method for solving the logical con guration problem. We have implemented and tested this method. In x4.1 the model from [7] is described. x4.2 and x4.3 describe a solution strategy (proposed by Farago et al [7]) for solving the logical con guration problem. Our implementation of this solution strategy is discussed in x4.4 and x4.5. In x4.6 we present some results obtained with this method.
4.1 The Model Consider a broadband network with N nodes and L physical links which carries K classes of trac. A number of virtual subnetworks will now be constructed on top of the physical network. Let K = f1; : : : ; K g be the set of classes, L = f1; : : :; Lg the set of physical links, V the set of virtual subnetworks and Pv the set of all O-D pairs in virtual subnetwork v 2 V . Each virtual subnetwork consists of virtual links . Let J denote the set of all virtual links over all virtual subnetworks and let J = jJ j. Let Cj denote the logical capacity (number of bandwidth units) of virtual link j . The correspondence between physical and virtual links is expressed by the L J Boolean matrix S , where the `j th element of S is
s`j =
1
0
if physical link ` is used by virtual link j; otherwise:
Let C`phy denote the capacity of physical link `. De ne the following two column vectors: 29
30
Chapter 4. Route Optimization Using Virtual Subnetworks
0C 1 1 B . . C=@ . C A; CJ
0C phy 1 1 B phy C = @ ... C A:
C is known as the logical capacity vector .
CLphy
An important constraint on the logical capacities is that the sum of the logical capacities of the virtual links using a particular physical link cannot exceed the physical capacity of the link:
S C Cphy
(4.1)
where the comparison is component-wise. Central to this model is the notion of a route . A route r is a subset of the set of logical links1 that belong to the same virtual subnetwork. Since each route carries a single trac class, parallel routes are required if two or more classes of trac have to be carried between an O-D pair. The correspondence between routes, virtual links and trac classes is expressed by the indicator variables :
Akjr =
1
0
when route r uses logical link j and carries class-k trac; otherwise:
Let akj denote the number of bandwidth units that a class-k call requires on virtual link j . There are two important implications of this notation:
all routes that carry class-k trac require the same number of bandwidth
units on virtual link j the bandwidth requirement of a call on a given route may vary along the virtual links constituting the route.
Let Rvpk denote the set of routes in virtual subnetwork v carrying class-k trac between O-D pair p. Let R denote the set of routes over all virtual subnetworks and let R = jRj:
R=
[ [ [ v2V p2P k2K
Rvpk :
v
Let Aj = fr 2 R j j 2 rg denote the set of all routes that use virtual link j . The function K : R 7! K maps routes onto classes (K (r) := k if route r carries class-k trac). Assume calls arrive to route r according to a Poisson process with parameter (arrival rate) r . Let 1=r be the average holding time of a call on route r, let r = r =r denote the load o ered to route r and vpk the aggregated o ered load of class-k calls to O-D pair p in virtual subnetwork v. Since (possibly) several routes carry trac of the same class between the same O-D pair, the load sharing or proportion in which the load is shared 1
The words logical link and virtual link are used interchangeably.
4.1 The Model
31
between the di erent routes has to be speci ed. Farago et al [7] de ne the load sharing vector svpk for virtual subnetwork v, O-D pair p and trac class k as follows: vpk svpk = ?svpk r1 : : : sr R
where svpk r = 0 for all r 62 Rvpk (the class-k load o ered to O-D pair p in virtual
subnetwork v can only be shared between class-k routes connecting O-D pair p in v). Since the load sharing factors svpk r represent the proportion in which the o ered load is shared and do not take losses into account, the following constraint holds:
X
r2R
svpk r =
X
r2R
svpk r = 1:
vpk
The logical con guration problem [7] is to dimension the virtual links and set the load sharing factors in such a way that the rate at which the network earns revenue is maximized. In this model the logical con guration problem can be stated as a non-linear maximization problem with the total network revenue rate as the objective function. If a call connected on route r generates revenue at a rate wr , the total network revenue rate W (; C) can be expressed as
W (; C) =
X
r2R
wr r (1 ? Lr (C))
where Lr (C) is the end-to-end blocking probability for trac on route r. This route blocking probability is de ned as the probability of the event that at least one virtual link j along the route r is blocked. The optimization problem is to set the logical link capacities Cj and the load sharing factors svpk r such that the revenue rate W is maximized:
Optimization Problem 4.1.
Maximize:
W=
XXX X v2V p2P k2K r2R v
Subject to:
svpk 0;
X r2R
wr vpk svpk r (1 ? Lr (C))
(4.2)
vpk
S C Cphy ; C 0 svpk r =1
for all v 2 V ; p 2 Pv ; k 2 K:
vpk
It is noted by Farago et al [7] that the objective (4.2) is dicult to optimize due to the dependence of the route blocking probabilities on the logical capacity vector C. Farago et al [7] propose solving optimization problem 4.1 by means of a two-stage method. In stage 1 (described in x4.2) the revenue rate is approximated by a concave function. In stage 2 (described in x4.3) the maximum obtained in stage 1 is used as the starting point for the optimization of a more re ned model.
32
Chapter 4. Route Optimization Using Virtual Subnetworks
4.2 Convex Approximation In the rst stage of the solution method proposed by Farago et al [7], optimization problem 4.1 is modi ed by assuming that the load sharing factors are xed. The revenue rate function is now a function of the logical capacities only:
Optimization Problem 4.2 (Fixed Load Sharing). Maximize:
W (C) =
X r2R
wr r (1 ? Lr (C))
(4.3)
Subject to:
S C Cphy ; C 0: Farago et al [7] derive a convex approximation of the blocking probability
Lr (C) for the cases of homogeneous and heterogeneous bandwidth demands
(single-service and multi-service trac).
4.2.1 Homogeneous Bandwidth Demands In the case of single-service trac, Farago et al [7] introduce a modi ed objective function:
f(C) = W
X r2R
wr r (1 ? Le r (C)):
(4.4)
The only di erence between (4.3) and (4.4) is that Lr (C) is replaced by
1 0 X X r0 ; Cj A Ler (C) = E @ j 2r
r0 2A
j
where E (; C ) is the analytic continuation of Erlang's function as a function of C [7]. The relation between r, r0 and j is clari ed by gure 4.1.
r 0 2 Aj
J _ _ t _ Jt _ Jt _ _ j_ _ Jt _ tJ _ tJ _ _ r r00 2 Aj Figure 4.1: Relation Between r, r0 and j Farago et al [7] prove the following theorem:
4.3 Re ned Model
33
f(C) (from (4.4)) is a concave Theorem 4.1 (Concave Lower Bound). W function of the logical capacity vector C. In the homogeneous case the following inequality holds:
f(C) W (C) W where W (C) is the exact revenue rate (from (4.3)) computed without assuming link independence. Proof. See [7]. From Theorem 4.1 and (4.3) the following inequality is obtained [7]:
X r2R
f(C): wr r W (C) W
According to Farago et al [7], the above inequality implies that the lower f(C) will be close to the exact revenue rate for small blocking probabound W bilities.
4.2.2 Heterogeneous Bandwidth Demands
Farago et al [7] extend the approximation used in the single-service case to multi-service trac. Therefore the modi ed objective function (4.4) is again used. However, Le r (C) is now given by
Ler (C) =
XX j 2r k2K
0 1 X X X vph akj Akjr E @ ahj vph sr0 ; Cj A r0 2A p2P h2K j
(4.5)
v
where r is a route in virtual subnetwork v. Figure 4.1 still applies. In contrast to homogeneous bandwidth demands for which the result of optimization problem 4.2 (using the modi ed objective function (4.4)) is the maximum of the lower bound on the revenue rate, this is not necessarily the case for heterogeneous bandwidth demands [7]. However, Ler (given by (4.5)) f is concave for both single-service and multi-service is convex and therefore W trac.
4.3 Re ned Model In the second stage of the solution method a re ned model is developed employing the reduced load approximation. In virtual subnetwork v, let bjk be the class-k load o ered to virtual link j when blocking is taken into account:
bjk = (1 ? Bjk )?1
X r2A
j
Akjr akj r (1 ? Lr )
(4.6)
34
Chapter 4. Route Optimization Using Virtual Subnetworks where Bjk denotes the blocking probability for class-k trac on virtual link j . In [7] it is assumed that a blocking function Ejk (Ejk (bj1 ; : : : ; bjK ; Cj ) returns the blocking probability of class-k trac on virtual link j ) exists2 . By applying the reduced load approximation and the link independence assumption, Farago et al [7] obtain the following:
Bjk = Ejk (bj1 ; : : : ; bjK ; Cj ) for all k 2 K and j 2 J %r = r (1 ? Lr ) YY 1 ? Lr = (1 ? Bjk )A kjr
bjk =
j 2J k2K
X
r2A
(1 ? Bjk )?1 %r Akjr akj for all k 2 K and j 2 J :
j
(4.7)
A set of auxiliary parameters ckj is now de ned as:
1 0 X X Bw ? X X A a c CC ckj = hkj Ahjr %r B gmr gm gm A @ r m2J g2K m6=j
r2R
h2K
(4.8)
where jh ?1 hkj = (1 ? Bjh )?1 @E @ b (1 ? Bjk ) : jk
The partial derivative of the revenue rate with respect to the logical capacity
Cj is
1 0 XX C; @W = X X A % B w ? Agmr agm cgm C r hj hjr r B A @
@Cj
h2K
r2R
jh hj = ?(1 ? Bjh )?1 @E @C :
m2J g2K m6=j
(4.9)
j
The partial derivative of the revenue rate with respect to the load o ered to route r is
1 0 X X @W = (1 ? L ) @w ? Akjr akj ckj A : r r @r
k2K j 2J
(4.10)
2 Farag o et al [7] mention the Erlang-B formula (for single-service trac), the KaufmanRoberts blocking formula for the basic stochastic knapsack and even a Gaussian approximation as possible candidates.
4.4 Algorithmic Implementation
35
4.4 Algorithmic Implementation In this section we discuss our implementation of the method proposed by Farago et al [7]. The implementation is discussed on the algorithmic level; x4.5 contains the details of the computer implementation.
4.4.1 The Model
The model is essentially the same as in x4.1 except for some restrictions and related changes in notation. The rst restriction is that any given virtual link j can use at most one3 physical link `. The logical capacity vector C is now replaced by the V L logical capacity matrix C . The v`th element of C (cv` ) is the capacity of the virtual link in virtual subnetwork v that uses circuits from physical link ` (this virtual link is called the virtual link corresponding to physical link `). cv` = 0 if O-D pair p (where p is the O-D pair connected by link ` in the physical network) is not connected4 in virtual subnetwork v. In conclusion: the topology of any virtual subnetwork can be obtained by removing zero or more links from the physical network. The constraint on the logical capacities (4.1) is now replaced by
X v2V
cv` C`phy
for all ` 2 L:
Another result of this restriction is that Av` will now denote the set of routes in virtual subnetwork v that use the virtual link corresponding to physical link `. Associated with this restriction we introduce the mapping L : J 7! L, where L (j ) := ` if virtual link j uses physical link `. We also de ne the action of L on a route r = fj1 ; : : : ; jn g as L (r) = f`1; : : : ; `n g if L (ji ) = `i for i = 1; : : : ; n. Let Kv denote the set of trac classes carried by virtual subnetwork v. The second restriction is that the sets of trac classes carried by di erent subnetworks are disjoint:
Kv \ Kv0 = ; for v 6= v0 : This restriction also implies that the variable a introduced in x4.1 can be indexed di erently. ak` now denotes the number of bandwidth units that a class-k call requires on the virtual link corresponding to physical link `5 . The indicator variables are also modi ed: 3 This is not seen as a signi cant restriction since multi-link routes consisting of two or more virtual links are allowed. For example, if trac has to be carried between O-D pair (o; d) via a tandem node c, a route consisting of the virtual links corresponding to the physical links o ? c and c ? d can be used. 4 An O-D pair (o; d) is connected if there is a direct link o ? d in the network. We do not use connected in the graph theoretical sense where connected means that there is a path from o to d in the network (see Chartrand and Oellermann [6] for example). 5 Since there is only one virtual subnetwork v that carries class-k trac, v does not have to be speci ed.
36
Chapter 4. Route Optimization Using Virtual Subnetworks
81 > > < Ak`r = > > :
when route r uses the logical link corresponding to physical link ` and carries class-k trac;
0 otherwise:
Using the same motivation we change the indexing of Rvpk , vpk and svpk r to Rpk , pk and spk r respectively.
4.4.2 Approximation Using Fixed Load Sharing
We now discuss our implementation of stage one (see x4.2) of the method proposed by Farago et al [7]. The load sharing factors are assumed xed. Di erent from x4.2, our route blocking probability function Ler is not convex. With the above restrictions on the virtual links and trac classes, optimization problem 4.2 becomes:
Optimization Problem 4.3 (Fixed Load Sharing).
Maximize:
Subject to:
f(C ) = W X v2V
X r2R
wr r (1 ? Ler (C ))
cv` C`phy
cv` 0
Ler (C ) in (4.11) is given by
(4.11)
for all ` 2 L
for all v 2 V ; ` 2 L:
!#a " eLr (C ) = 1 ? Y 1 ? E X ah`r0 ; cv`
k`
j 2r
r0 2A
(4.12)
v`
where r is a route in virtual subnetwork v, k is the trac class carried on route r, j is the virtual link (in virtual subnetwork v) corresponding to physical link ` and h is the trac class carried by route r0 . The relation between r, r0 and j is illustrated by gure 4.2. Note that using L , (4.12) can also be expressed as:
!#a " eLr (C ) = 1 ? Y 1 ? E X ah`r0 ; cv`
k`
`2L (r)
r0 2A
v`
The partial derivative of the objective (revenue rate) (4.11) with respect to the logical capacity is used in the optimization process. It can easily be veri ed that:
4.4 Algorithmic Implementation
37
r0 2 Av`
J _ _ t _ Jt _ Jt _ _ j_ _ Jt _ tJ _ tJ _ _ r r00 2 Av` Figure 4.2: The Relation Between r, r0 , etc.
f X @ Ler @W @cv` = ? r2R wr r @cv` ; !#a ?1 !" X @ Le r = a @E X a 0 ; c 0 @cv` k` @cv` r0 2A h` r v` 1 ? E r0 2A ah` r ; cv` 2 0 13 k`
v`
a 0
Y 4 X 1?E@ ah`0 r0 ; cv`0 A5 :
v`
k`
`0 20L (r) ` 6=`
r0 2A 0
(4.13)
v`
Theorem 4.2. If r > 0, Ler (C ) given by (4.12) is not a convex function.
Proof. According to Theorem A.1, Ler is convex if and only if the following inequality holds
Ler (C ) Le r (C 0 ) + dLe r (C 0 ) (C ? C 0 ) (4.14) for all logical capacity matrices C and C 0 . We now show that for a certain choice of C and C 0 , (4.14) does not hold. Let cv` denote the v` th component of C ? C 0 , the di erential in (4.14) can then be expressed as
dLe r (C 0 ) (C ? C 0 ) =
X X @ Ler (C 0 )
cv` v2V `2L @cv` X e 0 = @ L@cr (C ) cv` v` `2L X @ Ler (C 0 ) = cv` `2L (r) @cv`
where r is a route in virtual subnetwork v. We have also used the fact that @ Le r (C 0 )=@cv` = P0 for ` 62 L (r). Let ev` = r0 2A ah` r0 (where K (r0 ) = h). ev` > 0 for ` 2 L (r) since r > 0. Using (4.13), the di erential becomes: v`
dLe r (C 0 ) (C ? C 0 ) =
X @ Ler (C 0 ) @cv` cv`
`2L (r)
38
Chapter 4. Route Optimization Using Virtual Subnetworks
X
=
ev` ; c0v` ) [1 ? E (e ; c0 )]a ?1 ak` @E (@c v` v`
`2L (r) Y
`0 20L (r) ` 6=`
k`
v`
[1 ? E (ev`0 ; c0v`0 )]a 0 cv` : k`
From (4.12) Ler (C ) and Le r (C 0 ) can be expressed as
Le r (C ) = 1 ? Ler (C 0 ) = 1 ?
Y `2L (r)
Y
`2L (r)
[1 ? E (ev` ; cv` )]a
k`
[1 ? E (ev` ; c0v` )]a : k`
We now show that if C 0 is the zero logical capacity matrix (all virtual capacities are zero) and C is any non-zero logical capacity matrix, (4.14) does not hold. Using the above expressions for the di erential, Le r (C ) and Le r (C 0 ) as well as E (ev` ; 0) = 1 we obtain
Le r (C ) = 1 ? 0 is a predetermined convergence criterion. As Press et al [16] point out, solving a system of non-linear equations is a dicult problem. Consider a two-dimensional example from [16], the following two equations have to be solved simultaneously:
f (x; y) = 0 g(x; y) = 0:
(D.2) In gure D.1 these functions are plotted in the xy-plane (the gure is similar to one in [16]). The solid lines represent the f = 0 contours and the dashed lines the g = 0 contours. The solutions of (D.2) are the points where the f = 0 and g = 0 contours intersect. A and B are two examples of solutions. However M is not a solution! g>0
g=0
f >0
M
B f =0
f >0
y f 0 g=0 g>0 x
Figure D.1: The Contours in the xy-plane
g
