- Email: [email protected]
Exact solution of the SONET Ring Loading Problem
Operations Research Letters 25 (1999) 119–129
www.elsevier.com/locate/orms
Exact solution of the SONET Ring Loading Problem Mauro Dell’Amico a; ∗ , Martine Labbe b , Francesco Maoli c a Dipartimento b
di Scienze e Metodi Dell’Ingegneria, UniversitÂa di Modena e Reggio Emilia, Viale Allegri, 15, 42100 Reggio Emilia, Italy Institut de Statistique et Recherche OpÃerationnelle and Service de MathÃematiques de la Gestion, UniversitÃe Libre de Bruxelles, CP 210/01, Boulevard de Triomphe, B-1050 Bruxelles, Belgium c Dipartimento di Elettronica e Informazione, Politecnico di Milano, P.zza L. da Vinci 32, 20133 Milano, Italy Received 1 June 1996; received in revised form 1 February 1999
Abstract In this paper we address the problem of planning the capacity of the local rings in synchronous optical networks (SONET). We present ecient lower and upper bound procedures and a branch and bound algorithm which is able to nd the exact c 1999 Elsevier Science B.V. All rights reserved. solution of large instances, employing short computing times. Keywords: Communication; Programming; Linear algorithm; Programming; Integer; Branch and bound
1. Introduction One of the most successful ways to exploit the great potential of ber optics technology is o ered by the recent introduction of the synchronous digital hierarchy (SDH) in Europe and synchronous optical network (SONET) in the US. This has suggested a number of new optimization problems, most of them in the eld of combinatorial optimization. We address here the problem of planning the capacity of a single bi-directional ring of n ber links connecting n add-and-drop multiplexers (ADM) corresponding to n customer (or set of customers) locations. This problem was rst studied in [1] where it is called the SONET Ring Loading Problem. To the author’s knowledge, no exact approach has yet been proposed for this problem, which is known to be ∗ Corresponding author. Fax: +39-0522 496 466. E-mail address: [email protected] (M. Dell’Amico)
NP-hard [1] but is not yet satisfactorily classi ed with respect to its approximability status: in particular, it is still unknown if a polynomial-time approximation scheme may exist. The problem is known to be polynomially solvable only in the case of unit [2]. In [4] a very ecient heuristic has been proposed for the SONET Ring Loading Problem: it is a simple Greedy rounding algorithm which achieves a load that is guaranteed to exceed the optimum by at most 32 times the maximum demand. Moreover, this algorithm nds the maximal solution whenever all demands are equal to one. In Section 2 we formulate the problem and present a simple heuristic having an approximation guarantee factor of 2. In Section 3 we give a strongly polynomial algorithm for the continuous relaxation of the integer formulation. This is exploited in Section 4 within an exact branch-and-bound approach, whose performance is extensively tested as shown by the computational results in Section 5.
c 1999 Elsevier Science B.V. All rights reserved. 0167-6377/99/$ - see front matter PII: S 0 1 6 7 - 6 3 7 7 ( 9 9 ) 0 0 0 3 1 - 0
120
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
2. Formulation and heuristic solution We are given a ring graph G = (V; E) that is a simple circuit with node set V = {1; : : : ; n} and edge set E = {1; : : : ; n}. We also consider a set of m demands D = {1; : : : ; m}. Each demand i is de ned by an origin node oi , a destination node di and a weight wi which represent the amount of trac which must be routed from oi to di along the ring. The entire traf c of a demand must be routed either clockwise or counterclockwise and it can not be split to be transmitted in both directions simultaneously. Each edge ‘ ∈ E may be initially loaded with a xed trac t‘ ¿0, which is independent of the demands in D (this trac is due to pre-routed demands). The problem we consider is to determine the routing of each demand in D, in order to minimize the total load of an edge. Without loss of generality, we can assume that all the input data are integers. In order to present a mathematical model we need to introduce some notation. For each demand i let i+ denote the path from node oi to node di , clockwise, and i− the path from oi to di , counterclockwise. For each demand i ∈ D let xi be a boolean variable which takes value one if the trac wi is routed clockwise, value zero otherwise. We also use the continuous variable z to identify the maximum load of an edge. The model is as follows: (P)
min z
(1)
s:t: z¿
X
wk xk
k:‘∈k+
+
X
wk (1 − xk ) + t‘ ;
∀‘ ∈ E;
(2)
k:‘∈k−
xk ∈ {0; 1}; ∀k ∈ D:
(3)
Objective function (1) minimizes the maximum load of an edge. Constraints (2) ensure that the value of variable z is not less than the load of each edge ‘ ∈ E. In particular, the rst sum in (2) gives the total load of edge ‘ due to the demands routed clockwise; the second sum gives the load due to the demands routed counterclockwise. Constraints (3) impose that the variables be binary.
The continuous relaxation C(P) of problem P is obtained by allowing the variables xi to assume any real value in [0; 1], i.e. substituting (3) with 06xk 61;
∀k ∈ D:
(30 )
A demand can then be split by routing part of the trac wi clockwise and the remaining trac counterclockwise. A feasible solution for problem P can be obtained from the solution of the continuous relaxation by rounding each variable with fractional value to the closest integer. This corresponds to routing each fractional demand k (i.e a demand associated with a value xk with 0 ¡ xk ¡ 1) clockwise (counterclockwise) if xk ¿0:5 (xk ¡ 0:5). We call HR the heuristic algorithm which determines an approximate solution by applying this simple rounding technique. Theorem 2.1. The value of the solution obtained by algorithm HR is at most twice the value of the optimal solution. Proof. Algorithm HR routes the entire trac of a demand k along the path k+ if xk ¿0:5 and along k− , otherwise. It follows that the trac on an edge in the approximate solution is at most twice the trac routed on that edge in the continuous relaxation. Thus the maximum load of an edge in the heuristic solution is at most 2L, where L is the value of the continuous relaxation. Since L is a lower bound on the optimum solution value, the result follows. Algorithm HR can be improved by introducing a post-processing of the rst K demands with largest weight, having fractional values. In particular, we start by rounding the values of the variables of all demands, except that of the K selected fractional ones. Then we enumerate the 2K solutions obtained with all the possible routing combinations of the K fractional demands and we choose the best integer solution. If K is a constant independent of the input data, the procedure can be implemented to run in polynomial time. We call this improved algorithm KHR. In the next section we show how to solve problem C(P) eciently. The e ectiveness of procedure KHR is evaluated in Section 5.
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
3. A strongly polynomial algorithm for the continuous relaxation Before giving the details of our ecient algorithm to solve problem C(P), we need to introduce some further notation. For each pair of distinct edges r ∈ E; s ∈ E let Cr; s denote the set of crossing demands, i.e. the set of demands which must be routed exactly on one of the two edges r; s. Formally, Cr; s = {k ∈ D : r ∈ k+ and s ∈ k− or r ∈ k− and s ∈ k+ }:
Any pair of distinct edges (r; s) partitions the node set into two subsets, and is thus called a straight cut of G. In the following we will speak indi erently of the pair of distinct edges (r; s) or of the cut (r; s). Let denote the set of all possible straight cuts of G. Given a demand index i ∈ D and a partial solution xk ; k = 1; : : : ; i, let X X wk xk + wk (1 − xk ); ‘ ∈ E F‘i = k6i:‘∈k+
k6i:‘∈k−
denote the total trac on edge ‘, due to the rst i demands. In a solution to C(P), the load of an edge ‘ is therefore F‘n + t‘ . Let zC∗ be the optimal solution to C(P). We start the description of our algorithm by introducing a simple lower bound on zC∗ . Given a cut (r; s), it can be seen that the trac of each crossing demand must be routed through edge r or edge s, so X 1 wk + tr + ts (4) L1 = max 2 (r; s)∈ k∈Cr; s
is a valid lower bound. Since the optimal load of an edge ‘ ∈ E is at least t‘ , an improved bound is L = max{L1; L2}
(5)
with L2 = max‘∈E t‘ . The bound L is the optimal solution value for problem C(P). To prove this, we give a strongly polynomial time algorithm, called CSONET, which determines a feasible solution to C(P). Our algorithm considers one demand at a time and routes it optimally by determining the maximum trac which can be routed clockwise (or counterclockwise) without exceeding the bound L. More precisely, we have m iterations, one for each demand. At iteration
121
i, the routing of the rst i − 1 demands is already xed so that the edges of the ring have preassigned loads and we consider demand i and the path i+ . For all edges and straight cuts in i+ we compute two values C1 and C2, similar to L1 and L2, which give the total amount of trac which must certainly be routed on some edge of i+ , disregarding demand i. We do the same for path i− and we determine the gaps + and − between the lower bound L and the above quantities. If + ¿− then we route the trac equal to the minimum of wi and + along i+ , and the possibly remaining trac of demand i along i− . A similar routing is made if − ¿+ , but with paths i+ and i− exchanged. The details of procedure CSONET are given in the pseudocode below. Algorithm CSONET begin P compute L1 = 12 max k∈Cr; s wk + tr + ts ; (r; s)∈
compute L = max(L1; max‘∈E t‘ ); for each ‘ ∈ E do F‘0 :=0; for i:=1 to m do P 1 wk + Fri−1 C1:= 2 max(r; s)∈i+ k¿i k∈Cr; s i−1 +Fs + tr + ts ; C2:=max‘∈i+ (t‘ + F‘i−1 )
+ :=L − max(C1; C2); P wk + Fri−1 C1:= 12 max(r; s)∈− k¿i i k∈Cr; s +Fsi−1 + tr + ts ; C2:=max‘∈− (t‘ + F‘i−1 ) i
− :=L − max(C1; C2); if + ¿− then + xi :=min(1; wi ); else − xi :=max(0; 1 − wi ); endif endfor end. The following theorem holds.
122
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
Theorem 3.1. Algorithm CSONET determines in strongly polynomial time an optimal solution for problem C(P) and its value; zC∗ ; is equal to L. The proof of the theorem is given in the appendix. Alternative proofs were independently developed, during the writing of this paper (see [4,5]). 4. An exact algorithm The results of the previous sections have been used to develop a depth- rst branch and bound algorithm for the exact solution of problem P. Our basic algorithm is as follows. We use a binary decision tree. At each node of the tree we solve problem C(P) modi ed by xing the routing of a set of demands. Among the demands routed with a fractional trac in the solution to C(P), we determine the one having maximum trac. Two child nodes are generated by assigning this demand either clockwise or counterclockwise. At the root node we set t‘ = 0 for each ‘ ∈ E, then we apply algorithm CSONET of the previous section to solve problem C(P) and to obtain the lower bound L. The continuous solution is used by algorithm KHR of Section 2 to obtain a heuristic solution of value UB. If L = UB we have nished, otherwise we start the exploration of the branch-decision-tree. When a demand i is routed along a path (i.e. a forward step is performed) we increase by wi the values t‘ for each ‘ ∈ . When a backtracking occurs the value wi is subtracted from t‘ , for all ‘ ∈ . Observe that the loads t‘ completely describe the trac due to a partial solution, therefore no modi cation in the algorithms of Sections 2 and 3 is required in order to compute lower and upper bounds for the intermediate nodes of the tree. Preliminary computational experiments performed with instances randomly generated as in [1] show that the time required by our algorithm to nd an optimal solution is one order of magnitude less than the time required by CPLEX to nd an integer solution which is within 0.5% of the optimum (for more details on the experiments with CPLEX see [1]). Stimulated by these results, we tried to improve the performance of our basic algorithm with several techniques.
Presorting of the demands: We observed experimentally that the value of the heuristic solution obtained with KHR is very close to the optimal value. Therefore we decided to change the branching rule. We select not the fractional variable with maximum trac, but the rst fractional variable in a presorted list. This list is de ned in such a way that the rst descent in the branch-decision-tree, from the root node to a leaf, builds the heuristic solution. Our guess is that a limited number of backtracking (i.e. routing changes) are required to obtain the optimal solution. We also tried a presorting by nonincreasing weights. Improvement of the lower bound: Given the solution to C(P), we consider the K largest demands with fractional value (where K is the same parameter introduced in Section 2 to obtain the improved heuristic algorithm). Similar to procedure KHR, we enumerate the 2K partial solutions obtained with all the possible routing combinations of these K demands and we apply the lower bound procedure to the remaining m−K demands of each partial solution. The minimum of these bounds is a valid lower bound which improves bound L. We apply this technique to all the nodes of the tree. Second lower bound: Given a straight cut (r; s) and the crossing demands Cr; s , we know that each of these demands must be routed through edge r or s. Thus, minimizing the maximum load on r and s, with an integer assignment of the demands in Cr; s , gives a valid lower bound on the optimum solution value. This problem can be transformed into a subset-sum problem in which a set of q items is given, each of which has an associated weight, and a subset S of the q items which maximizes the sum of the weights of the item in S without exceeding a given capacity c is sought. Given a straight cut (r; s), we de ne P an instance of subset-sum with capacity c = 12 i∈Cr; s wi and with an item i, having weight wi , associated with each demand i ∈ Cr; s . The optimal solution to the subset-sum problem minimizes the di erence between c and the sum of the weights of the items selected, so it also minimizes the difference between the items selected and discarded. In terms of the SONET problem, this solution minimizes the di erence of the trac loads on r and s and thus the maximum of the two. For each of the n(n − 1)=2 straight cuts we solve the corresponding subset-sum problem with a C translation of the
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
FORTRAN routine given in the book by Martello and Toth [3]. We apply this bound to all the nodes of the tree. All the above improvements of the basic branch and bound have been tested with several computational experiments. The impact of each technique on the computational times is as follows: • Each kind of presorting leads to smaller computational times for some instances and larger for other instances, so, on average, no improvement can be obtained with this technique. • The improvement of the lower bound halves the computational times. • The second lower bound has a larger value than that of the continuous relaxation only for a few instances. In these cases, it drastically reduces the size of the branch-decision-tree, but on average the computational time required to compute the bound is not compensated by the reduction of the tree. From the above preliminary computational experiments we have seen that signi cantly improvements are obtained only with a strengthening of the lower bound coming from continuous relaxation. In order to obtain the maximum information from this bound we developed the following ÿxing procedure which tries to x the values of the fractional variables to either zero or one. The procedure performs a series of iterations. At each iteration we consider, in turn, each demand, say j, with fractional value in the optimal solution to C(P). We route the entire trac of demand j clockwise updating the values t‘ accordingly and we again compute the lower bound L. If the lower bound is larger than or equal to the current best solution value then the demand must be routed counterclockwise to have a chance of obtaining a better solution value. If the value of the fractional variable has not been xed by this attempt, we route the demand j counterclockwise and we repeat the same steps to try to x it. When all the fractional demands have been considered we compute the new solution of the continuous relaxation of the problem obtained from P by xing the routing of some demands, and a new iteration is started. The algorithm stops when, at one iteration, no demand has been xed. The pseudocode follows.
123
Procedure FIXING input: upper bound value U ; repeat compute lower bound L; nd the optimal solution x∗ of C(P); NOFIX := TRUE; for each demand j with xj∗ 6= bxj∗ c do route demand j clockwise; compute lower bound L; if L¿U then x xj∗ :=0; let NOFIX := FALSE; else route demand j counterclockwise; compute lower bound L; if L¿U then x xj∗ :=1; let NOFIX:=FALSE; endif endfor until (NOFIX = TRUE) The time complexity of procedure FIXING is O(I (mn2 + Fmn)) where I is the number of iterations performed and F the number of fractional variables. In the worst case it is O(m3 n2 ), but in practice, since the number of fractional demands is small, it is quite fast. The use of this procedure at each node of the branch-decision-tree reduced the computational times by one order of magnitude. We also noted that with this procedure the improved lower bound no longer has e ect. So our nal code is the basic branch and bound with procedure FIXING applied at each node. 5. Computational results We coded in C language the heuristic procedure KHR of Section 2, the lower bound procedure CSONET of Section 3 and the branch and bound of the previous Section. We also implemented the two greedy heuristics described in [1]. The rst heuristic, that we call 1PG, is a One Phase Greedy which considers one demand at a time and tentatively routes it clockwise and counterclockwise. The demand is routed on the path which determines the minimum increase in the current solution value. The second heuristic, say 2PG, is a Two-Phase Greedy which tentatively routes the current demand clockwise and counterclockwise, then completes the partial solution with 1PG and routes the demand along the path
124
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
associated with the minimum solution value. This algorithm is slightly di erent from that in [1], since we use 1PG to determine the routing of each demand, whereas Cosares and Saniee apply 1PG only if there is a tie in the partial load determined by the routing of the current demand clockwise or counterclockwise. We generated random instances as in [1]: the weights of the demands are randomly selected according to the uniform distribution in [5,100] and origins and destinations are di erent integers also uniformly randomly selected. The results for this class of problems, with values of n=(20; 40; 60; 80; 100) and values of m ranging from 50 to 1000, are reported in Table 1. Let z H be the solution value of the generic heuristic H for a speci c instance, and let z ∗ be its optimal solution value (computed using our branch and bound algorithm). For each instance, we compute the percentage error with respect to the optimal solution, i.e. err = 100(z H − z ∗ )=z ∗ . The columns labeled avg and max , in the table, report, respectively, the average and the maximum percentage error. The columns labeled time give the average computing time in seconds on a Pentium PC with a clock at 100 MHz. For each entry in the tables 20 instances were generated and solved. The labels Favg and Fmax indicate, respectively, the average and maximum number of fractional variables obtained in the solution of the continuous relaxation. For these experiments we xed the parameter K of heuristic KHR at 5. Procedure KHR dominates the other heuristics and it is fast. Its average errors are very small, exceeding 1% only for small instances (n = 20; m6100, and n = 40; m = 100). The maximum error never exceeds a few percent. Heuristics 1PG and 2PG have errors one order of magnitude larger than that of KHR. 1PG has negligible computing time, but produces the worst solutions. The average error of 2PG is one-third to one-half less that of 1PG, but its computing time is two orders of magnitude larger. KHR has computing times larger than that of 1PG, but still very small: less than 5 s for the largest instances. It is worth noting that currently the real instances have n620 and m6200 and the largest weight of a demand is bounded by 100, so the rst three rows of our table give a sample of real-life problems. The same instances of Table 1 have been solved with the branch-and-bound algorithm of the previous
section. The results are given in Table 2. For each pair n; m we give the average and maximum computing time (over 20 instances), the average and maximum number of explored nodes, the average and maximum depth of the branch-decision-tree and the ratio between the value of the lower bound at the root node (L1) and the optimal solution value (z ∗ ). All the instances with size corresponding to that of the largest real-life problems (n620 and m6200) are solved by the exact approach with small computing times. Larger instances are also solved in reasonable computing times, but the maximum times increases up to 1 h. The average number of explored nodes is, in general, small. This is mainly due to our xing procedure (see Section 4). For large instances, however, the lower bound at the root node, in many cases, is equal to the optimum solution value, but the algorithm determines a feasible solution with this value only after visiting a large number of nodes. To better understand the behaviour of the heuristic KHR and of the branch-and-bound procedure we generated instances from three further classes: (i) wj ∈ [5; 1000]; (ii) wj ∈ [50; 100] and (iii) wj ∈ [500; 1000]. The performances of the three heuristics are very similar to those obtained with weights wj ∈ [5; 100], but the new instances are more dicult for the branch and bound. The results are reported in Tables 3–5: the numbers in brackets give the number of instances which were not solved within 1 h of CPU time (in these cases the averages are computed on the number of solved instances). Acknowledgements This work has been partially sponsored by Human Capital and Mobility Contract number ERBCHRCXT930087. The rst and the third authors have been partially supported by MURST and CNR of Italy. The third author also thanks the NSFR of Belgium for nancial support. Appendix Proof of Theorem 3.1. We rst show that CSONET runs in strongly polynomial time, than we prove that L = zC∗ .
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
125
Table 1 Uniformly random instances with wj ∈ [5; 100]: heuristic solutions KHR n 20 40 60 80 100
1PG
2PG
m
avg
max
Favg
Fmax
Time
avg
max
Time
avg
max
Time
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
2.3 1.2 0.5 1.1 0.6 0.3 0.7 0.5 0.2 0.5 0.3 0.1 0.6 0.2 0.1
6.3 2.7 1.0 2.9 1.7 0.6 1.8 1.0 0.5 1.4 0.7 0.3 1.2 0.6 0.2
1.55 1.60 1.20 1.75 1.40 1.30 1.50 1.45 1.00 1.65 0.95 1.15 1.65 1.11 1.12
3 3 3 3 2 3 5 3 1 3 2 3 3 2 2
0.01 0.04 0.04 0.08 0.20 0.46 0.36 0.78 1.60 0.90 1.81 3.58 1.72 2.29 4.81
16.8 20.2 21.2 19.2 24.0 25.8 18.5 23.6 25.0 21.1 24.1 26.7 23.2 24.0 27.3
40.0 50.1 32.2 29.3 31.6 31.2 31.6 32.6 32.6 33.9 37.6 31.0 28.8 27.2 35.9
0.01 0.01 0.01 0.01 0.01 0.04 0.01 0.04 0.12 0.01 0.08 0.28 0.04 0.09 0.22
9.7 12.0 12.0 12.3 14.3 17.9 9.6 14.6 18.0 13.4 16.0 20.4 15.6 15.9 20.8
20.1 21.6 18.0 22.7 19.6 24.1 17.5 19.6 25.0 21.3 21.1 23.7 21.1 20.2 23.0
0.04 0.12 0.69 0.33 1.44 5.80 1.19 4.86 19.44 2.18 8.73 34.85 4.23 16.71 67.01
Note: Pentium/100 s; averages over 20 instances. Table 2 Uniformly random instances with wj ∈ [5; 100]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
0.26 1.02 2.68 2.49 7.71 36.40 9.83 24.17 233.94 25.30 135.31 251.03 41.97 162.16 304.84
0.82 2.55 7.86 7.70 17.71 81.89 22.16 71.92 2541.64 109.95 682.79 949.85 129.26 752.41 1117.81
35.8 44.6 23.0 46.0 38.0 87.0 42.3 43.1 674.4 56.8 185.3 44.8 52.0 253.8 72.0
215 127 77 209 133 323 103 263 12761 507 1637 219 208 2191 338
6.3 10.2 9.3 10.1 12.2 21.4 14.1 13.5 19.1 13.7 19.2 18.5 15.4 22.7 24.0
12 22 33 25 28 57 33 43 72 50 48 80 32 63 114
0.994 0.997 0.993 0.998 0.996 0.996 0.999 0.999 0.999 0.999 0.999 0.999 0.992 0.999 0.999
Note: Pentium/100 s; averages over 20 instances.
It is not dicult to be convinced that values L1; L; C1 and C2 can be computed in polynomial time, with respect to n and m. Moreover, since the main loop for executes m iterations, the overall algorithm terminates within a time bounded by a polynomial in n and m. Hence the algorithm is strongly polynomial. We now show that the value of the solution determined by CSONET is equal to L.
At the end of each iteration i of the main loop for in CSONET, the following two conditions hold: X 1 wk + Fri + Fsi + tr + ts ∀(r; s) ∈ ; L¿ 2 k¿i k∈Cr; s (6) L¿t‘ + F‘i ;
∀‘ ∈ E:
(7)
126
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129 Table 3 Uniformly random instances with wj ∈ [5; 1000]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
0.16 0.47 3.93 2.99 20.27 39.60 15.09 24.04 55.00 20.43 73.86 260.32 50.29 161.07 374.79
1.57 1.57 15.71 6.29 103.71 80.14 33.00 73.86 128.86 48.71 198.00 (1)831.29 136.71 781.00 1145.57
399.9 376.3 636.7 516.3 1396.8 557.8 903.6 310.6 160.3 333.6 308.1 4813.8 541.6 782.0 197.1
1569 1195 1632 1123 9094 1989 2490 1414 666 894 943 33127 1608 6512 666
199.7 187.9 318.1 257.9 698.0 278.6 451.7 155.0 79.8 166.5 153.8 2406.6 270.6 390.9 98.3
784 597 816 561 4547 994 1245 707 333 447 471 16563 804 3256 333
0.998 0.999 1.000 0.997 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.998 1.000 1.000 1.000
Note: Pentium/100 s; averages over 20 instances. Table 4 Uniformly random instances with wj ∈ [50; 100]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
0.31 2.20 1.26 4.09 16.81 33.52 8.96 42.74 106.86 60.34 46.04 143.47 348.07 79.99 229.74
3.14 7.86 4.71 17.29 106.86 (1)138.29 26.71 262.43 660.00 180.71 161.86 424.29 2669.86 220.00 534.29
330.3 916.2 91.8 498.3 1068.5 746.5 479.6 735.7 195.7 1181.6 366.8 145.9 42582.5 85.2 135.7
2382 3244 441 1735 7902 4394 2020 6272 1516 3952 2605 609 394236 229 717
164.9 458.0 45.7 249.0 534.1 373.3 239.6 367.6 97.8 590.7 183.2 72.8 21291.1 42.3 67.8
1191 1622 220 867 3951 2197 1010 3136 758 1976 1302 304 197118 114 358
0.994 1.000 0.999 1.000 0.999 0.999 0.998 0.999 1.000 0.999 0.999 1.000 0.999 1.000 1.000
Note: Pentium/100 s; averages over 20 instances.
If condition (7) is satis ed at the end of the nth iteration, then no edge ‘ has load greater than L. Since value L is a lower bound on the optimal solution value to C(P), then the solution provided by algorithm CSONET is optimal and the result follows. We now prove, by induction, that conditions (6) – (7) hold at each iteration. We will give the details of
the proof for the case + ¿− , the proof for the other case (+ ¡ − ) can be obtained by using similar arguments. The two conditions (6) and (7) are certainly satis ed at the beginning of the procedure (iteration 0), by de nition of L. Now suppose they hold at the end of iteration i − 1 (0 ¡ i6n): we will show that they also hold at the end of iteration i.
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
127
Table 5 Uniformly random instances with wj ∈ [500; 1000]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
4.40 4.09 29.86 15.09 47.14 348.07 107.49 247.19 226.60 234.77 454.61 1331.85 480.80 256.14 1002.67
15.71 14.14 103.71 48.71 201.14 1796.14 271.86 595.57 1293.29 895.71 3034.43 (2)2742.14 1535.28 3177.45 (2)2819.56
3442.7 1838.3 6109.6 2843.7 4274.4 11493.3 7317.2 24586.2 9066.6 6001.7 24827.5 8428.1 7704.1 11912.3 32456.2
11066 6947 20117 10667 13821 57066 22844 193262 32528 12460 77485 22283 36110 52994 113245
1721.1 919.0 3054.7 1421.6 2137.0 5746.5 3658.4 12293.0 4533.1 3000.7 12413.5 4214.3 3851.8 5956.0 11087.4
5533 3474 10058 5333 6910 28533 11422 96631 16264 6230 38742 11141 18055 26497 112341
0.995 0.997 0.999 0.998 0.997 1.000 0.999 1.000 1.000 0.998 1.000 0.998 0.998 0.999 0.99
Note: Pentium/100 s; averages over 20 instances.
From the de nition of C1 and C2, and the fact that conditions (6) and (7) are satis ed at the end of iteration i−1, it follows that C16L; C26L, so + ¿0 and the algorithm routes wi xi (=min(wi ; + )) trac units of demand i, on each edge ‘ ∈ i+ . We rst prove that the thesis holds for the edges and the straight cuts in i+ . From the above considerations we have wi xi 6+ 6L − C2 and C2¿t‘ + F‘i−1 for each ‘ ∈ i+ , so that L¿t‘ + F‘i−1 + wi xi . Since F‘i = F‘i−1 + wi xi , condition (7) immediately holds. Similarly, for condition (6) we have wi xi 6 + 6L − C1 1X wk + Fri−1 + Fsi−1 + tr + ts 6L− 2 k¿i k∈Cr; s
for each r ∈ i+ , s ∈ i+ . Since demand i is not included in Cr; s for all r ∈ i+ and s ∈ i+ we have X 1 wk + Fri−1 + Fsi−1 + 2wi xi + tr + ts L¿ 2 k¿i =
k∈Cr; s
1X wk + Fri + Fsi + tr + ts 2 k¿i k∈Cr; s
and condition (6) holds for any straight cut in i+ .
We now prove that the two conditions (6), (7) hold for the edges and the straight cuts in i− . If xi = 1 then no trac from demand i is routed on path i− , so the only remaining case is 0 ¡ + ¡ wi (i.e. the algorithm routes + trac units along path i+ and wi − + along i− ). By de nition of + we know that after iteration i at least one of the following two cases arises: case (i) an edge ‘ in i+ is saturated i.e. t‘ +F‘i−1 + wi xi = t‘ + F‘i = L; case (ii) a straight cut (r; s) with r ∈ i+ , s ∈ i+ is saturated, i.e. X 1 wk + Fri−1 + Fsi−1 + tr + ts + 2+ 2 k¿i k∈Cr; s
X 1 wk + Fri + Fsi + tr + ts = L: = 2 k¿i k∈Cr; s
In the remainder of the proof we refer to four particular edges r; s; u; v (see Fig. 1), with r ∈ i+ , s ∈ i+ , u ∈ i− , v ∈ i− and such that edge r precedes edge s, if we move clockwise from oi to di , and u precedes v, if we move counter-clockwise from oi to di .
128
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
for j = v and j = u, respectively, and subtracting (8) we obtain X 1X wk + wk + 2wi + Fui−1 + Fvi−1 2 k¿i k¿i k∈Cu; r
k∈Cv; r
+tu + tv − 2+ 6L: Recalling again that Fui = Fui−1 + wi − + and that Fvi = Fvi−1 + wi − + , we obtain X 1X wk + wk + Fui + Fvi + tu + tv 6L: 2 k¿i k¿i
Fig. 1. The four particular edges r; s; u and v.
k∈Cu; r
Case (i): Let r denote the edge in the routing of demand i, i.e
i+
saturated by
tr + Fri−1 + + = L:
(8)
Since condition (6) holds, at the end of iteration i − 1, for any cut (r; j) with r ∈ i+ and j ∈ i− , we have X 1 wk + Fji−1 + Fri−1 + tj + tr 2 k¿i k∈Cj; r
X 1 wk + wi + Fji−1 + Fri−1 + tj + tr = 2 k¿i
It is now sucient to observe that Cu; v ⊆ Cu; r ∪ Cv; r in order to derive X 1 wk + Fui + Fvi + tu + tv 6L; 2 k¿i k∈Cu; v
which shows that condition (6) holds for any straight cut (u; v) in i− . Case(ii): Let (r; s) with r ∈ i+ and s ∈ i+ denote the saturated cut, i.e. X 1 wk + Fri−1 + Fsi−1 + tr + ts + + = L 2 k¿i k∈Cr; s
(10)
k∈Cj; r
6L:
k∈Cv; r
(9)
Multiplying inequality (9) by 2 and subtracting Eq. (8) we obtain X wk + wi + Fji−1 − + + tj 6L: k¿i k∈Cj; r
Observing that Fji = Fji−1 + wi − + , for each edge j ∈ i− , we immediately conclude that tj + Fji 6L, and condition (7) holds. Condition (6) must be considered for straight cuts in i− only if i− contains at least two distinct edges, say u and v. Therefore we can refer to the two cuts (u; r) and (v; r). Adding two inequalities of type (9)
i−
and the two cuts (j; r) and consider an edge j ∈ and (j; s). Note that inequality (9) can be rewritten with s instead of r, so that by adding the two inequalities of type (9) associated with cuts (j; r) and (j; s), and observing that: Cr; s ⊆ Cj; r ∪ Cj; s we obtain tj + wi + Fji−1 X 1 + wk + Fri−1 + Fsi−1 + tr + ts 62L: 2 k¿i k∈Cr; s
Subtracting Eq. (10) and recalling that: (a) demand i 6∈ Cr; s ; (b) Fji =Fji−1 +wi −+ , we obtain tj +Fji 6L and condition (7) holds. We conclude the proof by showing that condition (6) also holds for any cut (u; v) with u ∈ i− and
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
v ∈ i− . Consider the two cuts (v; r) and (u; s) and observe that Cu; v ∪ Cr; s ⊆ Cv; r ∪ Cu; s . Adding the two inequalities of type (9) associated to the cuts (u; s) and (v; r), and subtracting Eq. (10) we obtain 1X wk + Fui−1 + Fvi−1 +tu +tv + 2(wi − + ) 2 k¿i k∈Cu; v
6L: Since F‘i = F‘i−1 + w‘ − + , ∀‘ ∈ i− condition (6) immediately follows. A straightforward implementation of algorithm CSONET performs m iterations (one for each demand) and at each iteration the quantities C1 and C2 are computed twice. Each evaluation of C2 requires O(n) time, whereas we can certainly determine C1 in O(mn2 ). Hence the naive implementation of CSONET runs in O(m2 n2 ), but this time complexity can be reduced to O(n3 ) as follows. Observe that the correctness of CSONET does not depend on the order in which the demands are examined (for loop), therefore we can use a demand ordering enables more than one demand at a time to be routed. We need to consider two possibilities. First, let us suppose that L = L2 and let ‘ˆ be the edge which determines the value of L2. Since we know that the nal solution value is L, then no demand can be routed through this edge and a complete solution can be immediately obtained by routing each demand along the ˆ If otherwise L = L1, let path which does not contain ‘. (r; s) be the straight cut which determines value L1. Given any demand k such that r ∈ k+ and s ∈ k+ (resp. r ∈ k− and s ∈ k− ), CSONET de nes + = 0
129
(resp. − = 0) and it routes the entire demand k along k− (resp. k+ ). Therefore no demand k ∈ D\Cr; s can be routed through edge r or s. If D\Cr; s =∅, a complete feasible solution can be immediately obtained by setting xi = 0:5 for each demand i ∈ Cr; s . Otherwise, we order the demands so that those belonging to D \ Cr; s are considered rst, hence routed at the same time. If we have not xed the routing of all demands we iteratively continue to apply the above reasoning and order the remaining demands accordingly. The maximum number of iterations performed by the above implementation is n. Indeed one can prove that, at a certain iteration, at least one of the two edges of the cut which determines the value L1(=L) does not belong to a cut considered in a previous iteration. Since it is possible to compute the values of all the straight cuts in O(n2 ) (see e.g. [5]), CSONET can be implemented to run in O(n3 ) time. References [1] S. Cosares, I. Saniee, An optimization problem related to balancing loads on SONET rings, Telecommunication Systems 3 (1994) 165–181. [2] A. Frank, T. Nishiseki, N. Saito, H. Suzuki, E. Tardos, Algorithms for routing around a rectangle, Discrete Appl. Math. 40 (1992) 363–378. [3] S. Martello, P. Toth, Knapsack Problems, Algorithms and Computer Implementations, Wiley, Chichester, 1992. [4] A. Schrijver, P. Seymour, P. Winkler, The ring loading problem, SIAM J. Discrete Math. 11 (1999) 1–14. [5] R. Vachani, A. Shulman, P. Kubat, J. Ward, Multicommodity
ows in ring networks, INFORMS J. Comput. 8 (1996) 235–242.
www.elsevier.com/locate/orms
Exact solution of the SONET Ring Loading Problem Mauro Dell’Amico a; ∗ , Martine Labbe b , Francesco Maoli c a Dipartimento b
di Scienze e Metodi Dell’Ingegneria, UniversitÂa di Modena e Reggio Emilia, Viale Allegri, 15, 42100 Reggio Emilia, Italy Institut de Statistique et Recherche OpÃerationnelle and Service de MathÃematiques de la Gestion, UniversitÃe Libre de Bruxelles, CP 210/01, Boulevard de Triomphe, B-1050 Bruxelles, Belgium c Dipartimento di Elettronica e Informazione, Politecnico di Milano, P.zza L. da Vinci 32, 20133 Milano, Italy Received 1 June 1996; received in revised form 1 February 1999
Abstract In this paper we address the problem of planning the capacity of the local rings in synchronous optical networks (SONET). We present ecient lower and upper bound procedures and a branch and bound algorithm which is able to nd the exact c 1999 Elsevier Science B.V. All rights reserved. solution of large instances, employing short computing times. Keywords: Communication; Programming; Linear algorithm; Programming; Integer; Branch and bound
1. Introduction One of the most successful ways to exploit the great potential of ber optics technology is o ered by the recent introduction of the synchronous digital hierarchy (SDH) in Europe and synchronous optical network (SONET) in the US. This has suggested a number of new optimization problems, most of them in the eld of combinatorial optimization. We address here the problem of planning the capacity of a single bi-directional ring of n ber links connecting n add-and-drop multiplexers (ADM) corresponding to n customer (or set of customers) locations. This problem was rst studied in [1] where it is called the SONET Ring Loading Problem. To the author’s knowledge, no exact approach has yet been proposed for this problem, which is known to be ∗ Corresponding author. Fax: +39-0522 496 466. E-mail address: [email protected] (M. Dell’Amico)
NP-hard [1] but is not yet satisfactorily classi ed with respect to its approximability status: in particular, it is still unknown if a polynomial-time approximation scheme may exist. The problem is known to be polynomially solvable only in the case of unit [2]. In [4] a very ecient heuristic has been proposed for the SONET Ring Loading Problem: it is a simple Greedy rounding algorithm which achieves a load that is guaranteed to exceed the optimum by at most 32 times the maximum demand. Moreover, this algorithm nds the maximal solution whenever all demands are equal to one. In Section 2 we formulate the problem and present a simple heuristic having an approximation guarantee factor of 2. In Section 3 we give a strongly polynomial algorithm for the continuous relaxation of the integer formulation. This is exploited in Section 4 within an exact branch-and-bound approach, whose performance is extensively tested as shown by the computational results in Section 5.
c 1999 Elsevier Science B.V. All rights reserved. 0167-6377/99/$ - see front matter PII: S 0 1 6 7 - 6 3 7 7 ( 9 9 ) 0 0 0 3 1 - 0
120
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
2. Formulation and heuristic solution We are given a ring graph G = (V; E) that is a simple circuit with node set V = {1; : : : ; n} and edge set E = {1; : : : ; n}. We also consider a set of m demands D = {1; : : : ; m}. Each demand i is de ned by an origin node oi , a destination node di and a weight wi which represent the amount of trac which must be routed from oi to di along the ring. The entire traf c of a demand must be routed either clockwise or counterclockwise and it can not be split to be transmitted in both directions simultaneously. Each edge ‘ ∈ E may be initially loaded with a xed trac t‘ ¿0, which is independent of the demands in D (this trac is due to pre-routed demands). The problem we consider is to determine the routing of each demand in D, in order to minimize the total load of an edge. Without loss of generality, we can assume that all the input data are integers. In order to present a mathematical model we need to introduce some notation. For each demand i let i+ denote the path from node oi to node di , clockwise, and i− the path from oi to di , counterclockwise. For each demand i ∈ D let xi be a boolean variable which takes value one if the trac wi is routed clockwise, value zero otherwise. We also use the continuous variable z to identify the maximum load of an edge. The model is as follows: (P)
min z
(1)
s:t: z¿
X
wk xk
k:‘∈k+
+
X
wk (1 − xk ) + t‘ ;
∀‘ ∈ E;
(2)
k:‘∈k−
xk ∈ {0; 1}; ∀k ∈ D:
(3)
Objective function (1) minimizes the maximum load of an edge. Constraints (2) ensure that the value of variable z is not less than the load of each edge ‘ ∈ E. In particular, the rst sum in (2) gives the total load of edge ‘ due to the demands routed clockwise; the second sum gives the load due to the demands routed counterclockwise. Constraints (3) impose that the variables be binary.
The continuous relaxation C(P) of problem P is obtained by allowing the variables xi to assume any real value in [0; 1], i.e. substituting (3) with 06xk 61;
∀k ∈ D:
(30 )
A demand can then be split by routing part of the trac wi clockwise and the remaining trac counterclockwise. A feasible solution for problem P can be obtained from the solution of the continuous relaxation by rounding each variable with fractional value to the closest integer. This corresponds to routing each fractional demand k (i.e a demand associated with a value xk with 0 ¡ xk ¡ 1) clockwise (counterclockwise) if xk ¿0:5 (xk ¡ 0:5). We call HR the heuristic algorithm which determines an approximate solution by applying this simple rounding technique. Theorem 2.1. The value of the solution obtained by algorithm HR is at most twice the value of the optimal solution. Proof. Algorithm HR routes the entire trac of a demand k along the path k+ if xk ¿0:5 and along k− , otherwise. It follows that the trac on an edge in the approximate solution is at most twice the trac routed on that edge in the continuous relaxation. Thus the maximum load of an edge in the heuristic solution is at most 2L, where L is the value of the continuous relaxation. Since L is a lower bound on the optimum solution value, the result follows. Algorithm HR can be improved by introducing a post-processing of the rst K demands with largest weight, having fractional values. In particular, we start by rounding the values of the variables of all demands, except that of the K selected fractional ones. Then we enumerate the 2K solutions obtained with all the possible routing combinations of the K fractional demands and we choose the best integer solution. If K is a constant independent of the input data, the procedure can be implemented to run in polynomial time. We call this improved algorithm KHR. In the next section we show how to solve problem C(P) eciently. The e ectiveness of procedure KHR is evaluated in Section 5.
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
3. A strongly polynomial algorithm for the continuous relaxation Before giving the details of our ecient algorithm to solve problem C(P), we need to introduce some further notation. For each pair of distinct edges r ∈ E; s ∈ E let Cr; s denote the set of crossing demands, i.e. the set of demands which must be routed exactly on one of the two edges r; s. Formally, Cr; s = {k ∈ D : r ∈ k+ and s ∈ k− or r ∈ k− and s ∈ k+ }:
Any pair of distinct edges (r; s) partitions the node set into two subsets, and is thus called a straight cut of G. In the following we will speak indi erently of the pair of distinct edges (r; s) or of the cut (r; s). Let denote the set of all possible straight cuts of G. Given a demand index i ∈ D and a partial solution xk ; k = 1; : : : ; i, let X X wk xk + wk (1 − xk ); ‘ ∈ E F‘i = k6i:‘∈k+
k6i:‘∈k−
denote the total trac on edge ‘, due to the rst i demands. In a solution to C(P), the load of an edge ‘ is therefore F‘n + t‘ . Let zC∗ be the optimal solution to C(P). We start the description of our algorithm by introducing a simple lower bound on zC∗ . Given a cut (r; s), it can be seen that the trac of each crossing demand must be routed through edge r or edge s, so X 1 wk + tr + ts (4) L1 = max 2 (r; s)∈ k∈Cr; s
is a valid lower bound. Since the optimal load of an edge ‘ ∈ E is at least t‘ , an improved bound is L = max{L1; L2}
(5)
with L2 = max‘∈E t‘ . The bound L is the optimal solution value for problem C(P). To prove this, we give a strongly polynomial time algorithm, called CSONET, which determines a feasible solution to C(P). Our algorithm considers one demand at a time and routes it optimally by determining the maximum trac which can be routed clockwise (or counterclockwise) without exceeding the bound L. More precisely, we have m iterations, one for each demand. At iteration
121
i, the routing of the rst i − 1 demands is already xed so that the edges of the ring have preassigned loads and we consider demand i and the path i+ . For all edges and straight cuts in i+ we compute two values C1 and C2, similar to L1 and L2, which give the total amount of trac which must certainly be routed on some edge of i+ , disregarding demand i. We do the same for path i− and we determine the gaps + and − between the lower bound L and the above quantities. If + ¿− then we route the trac equal to the minimum of wi and + along i+ , and the possibly remaining trac of demand i along i− . A similar routing is made if − ¿+ , but with paths i+ and i− exchanged. The details of procedure CSONET are given in the pseudocode below. Algorithm CSONET begin P compute L1 = 12 max k∈Cr; s wk + tr + ts ; (r; s)∈
compute L = max(L1; max‘∈E t‘ ); for each ‘ ∈ E do F‘0 :=0; for i:=1 to m do P 1 wk + Fri−1 C1:= 2 max(r; s)∈i+ k¿i k∈Cr; s i−1 +Fs + tr + ts ; C2:=max‘∈i+ (t‘ + F‘i−1 )
+ :=L − max(C1; C2); P wk + Fri−1 C1:= 12 max(r; s)∈− k¿i i k∈Cr; s +Fsi−1 + tr + ts ; C2:=max‘∈− (t‘ + F‘i−1 ) i
− :=L − max(C1; C2); if + ¿− then + xi :=min(1; wi ); else − xi :=max(0; 1 − wi ); endif endfor end. The following theorem holds.
122
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
Theorem 3.1. Algorithm CSONET determines in strongly polynomial time an optimal solution for problem C(P) and its value; zC∗ ; is equal to L. The proof of the theorem is given in the appendix. Alternative proofs were independently developed, during the writing of this paper (see [4,5]). 4. An exact algorithm The results of the previous sections have been used to develop a depth- rst branch and bound algorithm for the exact solution of problem P. Our basic algorithm is as follows. We use a binary decision tree. At each node of the tree we solve problem C(P) modi ed by xing the routing of a set of demands. Among the demands routed with a fractional trac in the solution to C(P), we determine the one having maximum trac. Two child nodes are generated by assigning this demand either clockwise or counterclockwise. At the root node we set t‘ = 0 for each ‘ ∈ E, then we apply algorithm CSONET of the previous section to solve problem C(P) and to obtain the lower bound L. The continuous solution is used by algorithm KHR of Section 2 to obtain a heuristic solution of value UB. If L = UB we have nished, otherwise we start the exploration of the branch-decision-tree. When a demand i is routed along a path (i.e. a forward step is performed) we increase by wi the values t‘ for each ‘ ∈ . When a backtracking occurs the value wi is subtracted from t‘ , for all ‘ ∈ . Observe that the loads t‘ completely describe the trac due to a partial solution, therefore no modi cation in the algorithms of Sections 2 and 3 is required in order to compute lower and upper bounds for the intermediate nodes of the tree. Preliminary computational experiments performed with instances randomly generated as in [1] show that the time required by our algorithm to nd an optimal solution is one order of magnitude less than the time required by CPLEX to nd an integer solution which is within 0.5% of the optimum (for more details on the experiments with CPLEX see [1]). Stimulated by these results, we tried to improve the performance of our basic algorithm with several techniques.
Presorting of the demands: We observed experimentally that the value of the heuristic solution obtained with KHR is very close to the optimal value. Therefore we decided to change the branching rule. We select not the fractional variable with maximum trac, but the rst fractional variable in a presorted list. This list is de ned in such a way that the rst descent in the branch-decision-tree, from the root node to a leaf, builds the heuristic solution. Our guess is that a limited number of backtracking (i.e. routing changes) are required to obtain the optimal solution. We also tried a presorting by nonincreasing weights. Improvement of the lower bound: Given the solution to C(P), we consider the K largest demands with fractional value (where K is the same parameter introduced in Section 2 to obtain the improved heuristic algorithm). Similar to procedure KHR, we enumerate the 2K partial solutions obtained with all the possible routing combinations of these K demands and we apply the lower bound procedure to the remaining m−K demands of each partial solution. The minimum of these bounds is a valid lower bound which improves bound L. We apply this technique to all the nodes of the tree. Second lower bound: Given a straight cut (r; s) and the crossing demands Cr; s , we know that each of these demands must be routed through edge r or s. Thus, minimizing the maximum load on r and s, with an integer assignment of the demands in Cr; s , gives a valid lower bound on the optimum solution value. This problem can be transformed into a subset-sum problem in which a set of q items is given, each of which has an associated weight, and a subset S of the q items which maximizes the sum of the weights of the item in S without exceeding a given capacity c is sought. Given a straight cut (r; s), we de ne P an instance of subset-sum with capacity c = 12 i∈Cr; s wi and with an item i, having weight wi , associated with each demand i ∈ Cr; s . The optimal solution to the subset-sum problem minimizes the di erence between c and the sum of the weights of the items selected, so it also minimizes the difference between the items selected and discarded. In terms of the SONET problem, this solution minimizes the di erence of the trac loads on r and s and thus the maximum of the two. For each of the n(n − 1)=2 straight cuts we solve the corresponding subset-sum problem with a C translation of the
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
FORTRAN routine given in the book by Martello and Toth [3]. We apply this bound to all the nodes of the tree. All the above improvements of the basic branch and bound have been tested with several computational experiments. The impact of each technique on the computational times is as follows: • Each kind of presorting leads to smaller computational times for some instances and larger for other instances, so, on average, no improvement can be obtained with this technique. • The improvement of the lower bound halves the computational times. • The second lower bound has a larger value than that of the continuous relaxation only for a few instances. In these cases, it drastically reduces the size of the branch-decision-tree, but on average the computational time required to compute the bound is not compensated by the reduction of the tree. From the above preliminary computational experiments we have seen that signi cantly improvements are obtained only with a strengthening of the lower bound coming from continuous relaxation. In order to obtain the maximum information from this bound we developed the following ÿxing procedure which tries to x the values of the fractional variables to either zero or one. The procedure performs a series of iterations. At each iteration we consider, in turn, each demand, say j, with fractional value in the optimal solution to C(P). We route the entire trac of demand j clockwise updating the values t‘ accordingly and we again compute the lower bound L. If the lower bound is larger than or equal to the current best solution value then the demand must be routed counterclockwise to have a chance of obtaining a better solution value. If the value of the fractional variable has not been xed by this attempt, we route the demand j counterclockwise and we repeat the same steps to try to x it. When all the fractional demands have been considered we compute the new solution of the continuous relaxation of the problem obtained from P by xing the routing of some demands, and a new iteration is started. The algorithm stops when, at one iteration, no demand has been xed. The pseudocode follows.
123
Procedure FIXING input: upper bound value U ; repeat compute lower bound L; nd the optimal solution x∗ of C(P); NOFIX := TRUE; for each demand j with xj∗ 6= bxj∗ c do route demand j clockwise; compute lower bound L; if L¿U then x xj∗ :=0; let NOFIX := FALSE; else route demand j counterclockwise; compute lower bound L; if L¿U then x xj∗ :=1; let NOFIX:=FALSE; endif endfor until (NOFIX = TRUE) The time complexity of procedure FIXING is O(I (mn2 + Fmn)) where I is the number of iterations performed and F the number of fractional variables. In the worst case it is O(m3 n2 ), but in practice, since the number of fractional demands is small, it is quite fast. The use of this procedure at each node of the branch-decision-tree reduced the computational times by one order of magnitude. We also noted that with this procedure the improved lower bound no longer has e ect. So our nal code is the basic branch and bound with procedure FIXING applied at each node. 5. Computational results We coded in C language the heuristic procedure KHR of Section 2, the lower bound procedure CSONET of Section 3 and the branch and bound of the previous Section. We also implemented the two greedy heuristics described in [1]. The rst heuristic, that we call 1PG, is a One Phase Greedy which considers one demand at a time and tentatively routes it clockwise and counterclockwise. The demand is routed on the path which determines the minimum increase in the current solution value. The second heuristic, say 2PG, is a Two-Phase Greedy which tentatively routes the current demand clockwise and counterclockwise, then completes the partial solution with 1PG and routes the demand along the path
124
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
associated with the minimum solution value. This algorithm is slightly di erent from that in [1], since we use 1PG to determine the routing of each demand, whereas Cosares and Saniee apply 1PG only if there is a tie in the partial load determined by the routing of the current demand clockwise or counterclockwise. We generated random instances as in [1]: the weights of the demands are randomly selected according to the uniform distribution in [5,100] and origins and destinations are di erent integers also uniformly randomly selected. The results for this class of problems, with values of n=(20; 40; 60; 80; 100) and values of m ranging from 50 to 1000, are reported in Table 1. Let z H be the solution value of the generic heuristic H for a speci c instance, and let z ∗ be its optimal solution value (computed using our branch and bound algorithm). For each instance, we compute the percentage error with respect to the optimal solution, i.e. err = 100(z H − z ∗ )=z ∗ . The columns labeled avg and max , in the table, report, respectively, the average and the maximum percentage error. The columns labeled time give the average computing time in seconds on a Pentium PC with a clock at 100 MHz. For each entry in the tables 20 instances were generated and solved. The labels Favg and Fmax indicate, respectively, the average and maximum number of fractional variables obtained in the solution of the continuous relaxation. For these experiments we xed the parameter K of heuristic KHR at 5. Procedure KHR dominates the other heuristics and it is fast. Its average errors are very small, exceeding 1% only for small instances (n = 20; m6100, and n = 40; m = 100). The maximum error never exceeds a few percent. Heuristics 1PG and 2PG have errors one order of magnitude larger than that of KHR. 1PG has negligible computing time, but produces the worst solutions. The average error of 2PG is one-third to one-half less that of 1PG, but its computing time is two orders of magnitude larger. KHR has computing times larger than that of 1PG, but still very small: less than 5 s for the largest instances. It is worth noting that currently the real instances have n620 and m6200 and the largest weight of a demand is bounded by 100, so the rst three rows of our table give a sample of real-life problems. The same instances of Table 1 have been solved with the branch-and-bound algorithm of the previous
section. The results are given in Table 2. For each pair n; m we give the average and maximum computing time (over 20 instances), the average and maximum number of explored nodes, the average and maximum depth of the branch-decision-tree and the ratio between the value of the lower bound at the root node (L1) and the optimal solution value (z ∗ ). All the instances with size corresponding to that of the largest real-life problems (n620 and m6200) are solved by the exact approach with small computing times. Larger instances are also solved in reasonable computing times, but the maximum times increases up to 1 h. The average number of explored nodes is, in general, small. This is mainly due to our xing procedure (see Section 4). For large instances, however, the lower bound at the root node, in many cases, is equal to the optimum solution value, but the algorithm determines a feasible solution with this value only after visiting a large number of nodes. To better understand the behaviour of the heuristic KHR and of the branch-and-bound procedure we generated instances from three further classes: (i) wj ∈ [5; 1000]; (ii) wj ∈ [50; 100] and (iii) wj ∈ [500; 1000]. The performances of the three heuristics are very similar to those obtained with weights wj ∈ [5; 100], but the new instances are more dicult for the branch and bound. The results are reported in Tables 3–5: the numbers in brackets give the number of instances which were not solved within 1 h of CPU time (in these cases the averages are computed on the number of solved instances). Acknowledgements This work has been partially sponsored by Human Capital and Mobility Contract number ERBCHRCXT930087. The rst and the third authors have been partially supported by MURST and CNR of Italy. The third author also thanks the NSFR of Belgium for nancial support. Appendix Proof of Theorem 3.1. We rst show that CSONET runs in strongly polynomial time, than we prove that L = zC∗ .
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
125
Table 1 Uniformly random instances with wj ∈ [5; 100]: heuristic solutions KHR n 20 40 60 80 100
1PG
2PG
m
avg
max
Favg
Fmax
Time
avg
max
Time
avg
max
Time
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
2.3 1.2 0.5 1.1 0.6 0.3 0.7 0.5 0.2 0.5 0.3 0.1 0.6 0.2 0.1
6.3 2.7 1.0 2.9 1.7 0.6 1.8 1.0 0.5 1.4 0.7 0.3 1.2 0.6 0.2
1.55 1.60 1.20 1.75 1.40 1.30 1.50 1.45 1.00 1.65 0.95 1.15 1.65 1.11 1.12
3 3 3 3 2 3 5 3 1 3 2 3 3 2 2
0.01 0.04 0.04 0.08 0.20 0.46 0.36 0.78 1.60 0.90 1.81 3.58 1.72 2.29 4.81
16.8 20.2 21.2 19.2 24.0 25.8 18.5 23.6 25.0 21.1 24.1 26.7 23.2 24.0 27.3
40.0 50.1 32.2 29.3 31.6 31.2 31.6 32.6 32.6 33.9 37.6 31.0 28.8 27.2 35.9
0.01 0.01 0.01 0.01 0.01 0.04 0.01 0.04 0.12 0.01 0.08 0.28 0.04 0.09 0.22
9.7 12.0 12.0 12.3 14.3 17.9 9.6 14.6 18.0 13.4 16.0 20.4 15.6 15.9 20.8
20.1 21.6 18.0 22.7 19.6 24.1 17.5 19.6 25.0 21.3 21.1 23.7 21.1 20.2 23.0
0.04 0.12 0.69 0.33 1.44 5.80 1.19 4.86 19.44 2.18 8.73 34.85 4.23 16.71 67.01
Note: Pentium/100 s; averages over 20 instances. Table 2 Uniformly random instances with wj ∈ [5; 100]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
0.26 1.02 2.68 2.49 7.71 36.40 9.83 24.17 233.94 25.30 135.31 251.03 41.97 162.16 304.84
0.82 2.55 7.86 7.70 17.71 81.89 22.16 71.92 2541.64 109.95 682.79 949.85 129.26 752.41 1117.81
35.8 44.6 23.0 46.0 38.0 87.0 42.3 43.1 674.4 56.8 185.3 44.8 52.0 253.8 72.0
215 127 77 209 133 323 103 263 12761 507 1637 219 208 2191 338
6.3 10.2 9.3 10.1 12.2 21.4 14.1 13.5 19.1 13.7 19.2 18.5 15.4 22.7 24.0
12 22 33 25 28 57 33 43 72 50 48 80 32 63 114
0.994 0.997 0.993 0.998 0.996 0.996 0.999 0.999 0.999 0.999 0.999 0.999 0.992 0.999 0.999
Note: Pentium/100 s; averages over 20 instances.
It is not dicult to be convinced that values L1; L; C1 and C2 can be computed in polynomial time, with respect to n and m. Moreover, since the main loop for executes m iterations, the overall algorithm terminates within a time bounded by a polynomial in n and m. Hence the algorithm is strongly polynomial. We now show that the value of the solution determined by CSONET is equal to L.
At the end of each iteration i of the main loop for in CSONET, the following two conditions hold: X 1 wk + Fri + Fsi + tr + ts ∀(r; s) ∈ ; L¿ 2 k¿i k∈Cr; s (6) L¿t‘ + F‘i ;
∀‘ ∈ E:
(7)
126
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129 Table 3 Uniformly random instances with wj ∈ [5; 1000]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
0.16 0.47 3.93 2.99 20.27 39.60 15.09 24.04 55.00 20.43 73.86 260.32 50.29 161.07 374.79
1.57 1.57 15.71 6.29 103.71 80.14 33.00 73.86 128.86 48.71 198.00 (1)831.29 136.71 781.00 1145.57
399.9 376.3 636.7 516.3 1396.8 557.8 903.6 310.6 160.3 333.6 308.1 4813.8 541.6 782.0 197.1
1569 1195 1632 1123 9094 1989 2490 1414 666 894 943 33127 1608 6512 666
199.7 187.9 318.1 257.9 698.0 278.6 451.7 155.0 79.8 166.5 153.8 2406.6 270.6 390.9 98.3
784 597 816 561 4547 994 1245 707 333 447 471 16563 804 3256 333
0.998 0.999 1.000 0.997 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.998 1.000 1.000 1.000
Note: Pentium/100 s; averages over 20 instances. Table 4 Uniformly random instances with wj ∈ [50; 100]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
0.31 2.20 1.26 4.09 16.81 33.52 8.96 42.74 106.86 60.34 46.04 143.47 348.07 79.99 229.74
3.14 7.86 4.71 17.29 106.86 (1)138.29 26.71 262.43 660.00 180.71 161.86 424.29 2669.86 220.00 534.29
330.3 916.2 91.8 498.3 1068.5 746.5 479.6 735.7 195.7 1181.6 366.8 145.9 42582.5 85.2 135.7
2382 3244 441 1735 7902 4394 2020 6272 1516 3952 2605 609 394236 229 717
164.9 458.0 45.7 249.0 534.1 373.3 239.6 367.6 97.8 590.7 183.2 72.8 21291.1 42.3 67.8
1191 1622 220 867 3951 2197 1010 3136 758 1976 1302 304 197118 114 358
0.994 1.000 0.999 1.000 0.999 0.999 0.998 0.999 1.000 0.999 0.999 1.000 0.999 1.000 1.000
Note: Pentium/100 s; averages over 20 instances.
If condition (7) is satis ed at the end of the nth iteration, then no edge ‘ has load greater than L. Since value L is a lower bound on the optimal solution value to C(P), then the solution provided by algorithm CSONET is optimal and the result follows. We now prove, by induction, that conditions (6) – (7) hold at each iteration. We will give the details of
the proof for the case + ¿− , the proof for the other case (+ ¡ − ) can be obtained by using similar arguments. The two conditions (6) and (7) are certainly satis ed at the beginning of the procedure (iteration 0), by de nition of L. Now suppose they hold at the end of iteration i − 1 (0 ¡ i6n): we will show that they also hold at the end of iteration i.
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
127
Table 5 Uniformly random instances with wj ∈ [500; 1000]: exact solution Time n 20 40 60 80 100
Exp. Nodes
Depth
m
Avg
Max
Avg
Max
Avg
Max
L1=z ∗
50 100 200 100 200 400 150 300 600 200 400 800 250 500 1000
4.40 4.09 29.86 15.09 47.14 348.07 107.49 247.19 226.60 234.77 454.61 1331.85 480.80 256.14 1002.67
15.71 14.14 103.71 48.71 201.14 1796.14 271.86 595.57 1293.29 895.71 3034.43 (2)2742.14 1535.28 3177.45 (2)2819.56
3442.7 1838.3 6109.6 2843.7 4274.4 11493.3 7317.2 24586.2 9066.6 6001.7 24827.5 8428.1 7704.1 11912.3 32456.2
11066 6947 20117 10667 13821 57066 22844 193262 32528 12460 77485 22283 36110 52994 113245
1721.1 919.0 3054.7 1421.6 2137.0 5746.5 3658.4 12293.0 4533.1 3000.7 12413.5 4214.3 3851.8 5956.0 11087.4
5533 3474 10058 5333 6910 28533 11422 96631 16264 6230 38742 11141 18055 26497 112341
0.995 0.997 0.999 0.998 0.997 1.000 0.999 1.000 1.000 0.998 1.000 0.998 0.998 0.999 0.99
Note: Pentium/100 s; averages over 20 instances.
From the de nition of C1 and C2, and the fact that conditions (6) and (7) are satis ed at the end of iteration i−1, it follows that C16L; C26L, so + ¿0 and the algorithm routes wi xi (=min(wi ; + )) trac units of demand i, on each edge ‘ ∈ i+ . We rst prove that the thesis holds for the edges and the straight cuts in i+ . From the above considerations we have wi xi 6+ 6L − C2 and C2¿t‘ + F‘i−1 for each ‘ ∈ i+ , so that L¿t‘ + F‘i−1 + wi xi . Since F‘i = F‘i−1 + wi xi , condition (7) immediately holds. Similarly, for condition (6) we have wi xi 6 + 6L − C1 1X wk + Fri−1 + Fsi−1 + tr + ts 6L− 2 k¿i k∈Cr; s
for each r ∈ i+ , s ∈ i+ . Since demand i is not included in Cr; s for all r ∈ i+ and s ∈ i+ we have X 1 wk + Fri−1 + Fsi−1 + 2wi xi + tr + ts L¿ 2 k¿i =
k∈Cr; s
1X wk + Fri + Fsi + tr + ts 2 k¿i k∈Cr; s
and condition (6) holds for any straight cut in i+ .
We now prove that the two conditions (6), (7) hold for the edges and the straight cuts in i− . If xi = 1 then no trac from demand i is routed on path i− , so the only remaining case is 0 ¡ + ¡ wi (i.e. the algorithm routes + trac units along path i+ and wi − + along i− ). By de nition of + we know that after iteration i at least one of the following two cases arises: case (i) an edge ‘ in i+ is saturated i.e. t‘ +F‘i−1 + wi xi = t‘ + F‘i = L; case (ii) a straight cut (r; s) with r ∈ i+ , s ∈ i+ is saturated, i.e. X 1 wk + Fri−1 + Fsi−1 + tr + ts + 2+ 2 k¿i k∈Cr; s
X 1 wk + Fri + Fsi + tr + ts = L: = 2 k¿i k∈Cr; s
In the remainder of the proof we refer to four particular edges r; s; u; v (see Fig. 1), with r ∈ i+ , s ∈ i+ , u ∈ i− , v ∈ i− and such that edge r precedes edge s, if we move clockwise from oi to di , and u precedes v, if we move counter-clockwise from oi to di .
128
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
for j = v and j = u, respectively, and subtracting (8) we obtain X 1X wk + wk + 2wi + Fui−1 + Fvi−1 2 k¿i k¿i k∈Cu; r
k∈Cv; r
+tu + tv − 2+ 6L: Recalling again that Fui = Fui−1 + wi − + and that Fvi = Fvi−1 + wi − + , we obtain X 1X wk + wk + Fui + Fvi + tu + tv 6L: 2 k¿i k¿i
Fig. 1. The four particular edges r; s; u and v.
k∈Cu; r
Case (i): Let r denote the edge in the routing of demand i, i.e
i+
saturated by
tr + Fri−1 + + = L:
(8)
Since condition (6) holds, at the end of iteration i − 1, for any cut (r; j) with r ∈ i+ and j ∈ i− , we have X 1 wk + Fji−1 + Fri−1 + tj + tr 2 k¿i k∈Cj; r
X 1 wk + wi + Fji−1 + Fri−1 + tj + tr = 2 k¿i
It is now sucient to observe that Cu; v ⊆ Cu; r ∪ Cv; r in order to derive X 1 wk + Fui + Fvi + tu + tv 6L; 2 k¿i k∈Cu; v
which shows that condition (6) holds for any straight cut (u; v) in i− . Case(ii): Let (r; s) with r ∈ i+ and s ∈ i+ denote the saturated cut, i.e. X 1 wk + Fri−1 + Fsi−1 + tr + ts + + = L 2 k¿i k∈Cr; s
(10)
k∈Cj; r
6L:
k∈Cv; r
(9)
Multiplying inequality (9) by 2 and subtracting Eq. (8) we obtain X wk + wi + Fji−1 − + + tj 6L: k¿i k∈Cj; r
Observing that Fji = Fji−1 + wi − + , for each edge j ∈ i− , we immediately conclude that tj + Fji 6L, and condition (7) holds. Condition (6) must be considered for straight cuts in i− only if i− contains at least two distinct edges, say u and v. Therefore we can refer to the two cuts (u; r) and (v; r). Adding two inequalities of type (9)
i−
and the two cuts (j; r) and consider an edge j ∈ and (j; s). Note that inequality (9) can be rewritten with s instead of r, so that by adding the two inequalities of type (9) associated with cuts (j; r) and (j; s), and observing that: Cr; s ⊆ Cj; r ∪ Cj; s we obtain tj + wi + Fji−1 X 1 + wk + Fri−1 + Fsi−1 + tr + ts 62L: 2 k¿i k∈Cr; s
Subtracting Eq. (10) and recalling that: (a) demand i 6∈ Cr; s ; (b) Fji =Fji−1 +wi −+ , we obtain tj +Fji 6L and condition (7) holds. We conclude the proof by showing that condition (6) also holds for any cut (u; v) with u ∈ i− and
M. Dell’Amico et al. / Operations Research Letters 25 (1999) 119–129
v ∈ i− . Consider the two cuts (v; r) and (u; s) and observe that Cu; v ∪ Cr; s ⊆ Cv; r ∪ Cu; s . Adding the two inequalities of type (9) associated to the cuts (u; s) and (v; r), and subtracting Eq. (10) we obtain 1X wk + Fui−1 + Fvi−1 +tu +tv + 2(wi − + ) 2 k¿i k∈Cu; v
6L: Since F‘i = F‘i−1 + w‘ − + , ∀‘ ∈ i− condition (6) immediately follows. A straightforward implementation of algorithm CSONET performs m iterations (one for each demand) and at each iteration the quantities C1 and C2 are computed twice. Each evaluation of C2 requires O(n) time, whereas we can certainly determine C1 in O(mn2 ). Hence the naive implementation of CSONET runs in O(m2 n2 ), but this time complexity can be reduced to O(n3 ) as follows. Observe that the correctness of CSONET does not depend on the order in which the demands are examined (for loop), therefore we can use a demand ordering enables more than one demand at a time to be routed. We need to consider two possibilities. First, let us suppose that L = L2 and let ‘ˆ be the edge which determines the value of L2. Since we know that the nal solution value is L, then no demand can be routed through this edge and a complete solution can be immediately obtained by routing each demand along the ˆ If otherwise L = L1, let path which does not contain ‘. (r; s) be the straight cut which determines value L1. Given any demand k such that r ∈ k+ and s ∈ k+ (resp. r ∈ k− and s ∈ k− ), CSONET de nes + = 0
129
(resp. − = 0) and it routes the entire demand k along k− (resp. k+ ). Therefore no demand k ∈ D\Cr; s can be routed through edge r or s. If D\Cr; s =∅, a complete feasible solution can be immediately obtained by setting xi = 0:5 for each demand i ∈ Cr; s . Otherwise, we order the demands so that those belonging to D \ Cr; s are considered rst, hence routed at the same time. If we have not xed the routing of all demands we iteratively continue to apply the above reasoning and order the remaining demands accordingly. The maximum number of iterations performed by the above implementation is n. Indeed one can prove that, at a certain iteration, at least one of the two edges of the cut which determines the value L1(=L) does not belong to a cut considered in a previous iteration. Since it is possible to compute the values of all the straight cuts in O(n2 ) (see e.g. [5]), CSONET can be implemented to run in O(n3 ) time. References [1] S. Cosares, I. Saniee, An optimization problem related to balancing loads on SONET rings, Telecommunication Systems 3 (1994) 165–181. [2] A. Frank, T. Nishiseki, N. Saito, H. Suzuki, E. Tardos, Algorithms for routing around a rectangle, Discrete Appl. Math. 40 (1992) 363–378. [3] S. Martello, P. Toth, Knapsack Problems, Algorithms and Computer Implementations, Wiley, Chichester, 1992. [4] A. Schrijver, P. Seymour, P. Winkler, The ring loading problem, SIAM J. Discrete Math. 11 (1999) 1–14. [5] R. Vachani, A. Shulman, P. Kubat, J. Ward, Multicommodity
ows in ring networks, INFORMS J. Comput. 8 (1996) 235–242.