- Email: [email protected]
Beam search heuristic to solve stochastic integer problems under probabilistic constraints
European Journal of Operational Research 167 (2005) 35–47 www.elsevier.com/locate/dsw
Discrete Optimization
Beam search heuristic to solve stochastic integer problems under probabilistic constraints Patrizia Beraldi b
a,*
, Andrzej Ruszczyn´ski
b
a Dipartimento di Elettronica, Informatica e Sistemistica, Universita` della Calabria, 87036 Rende (CS), Italy RUTCOR––Rutgers Center for Operations Research, 640 Bartholomew Road, Piscataway, NJ 08854-8003, USA
Received 27 May 2002; accepted 23 February 2004 Available online 26 June 2004
Abstract This paper proposes a Beam Search heuristic strategy to solve stochastic integer programming problems under probabilistic constraints. Beam Search is an adaptation of the classical Branch and Bound method in which at any level of the search tree only the most promising nodes are kept for further exploration, whereas the remaining are pruned out permanently. The proposed algorithm has been compared with the Branch and Bound method. The numerical results collected on the probabilistic set covering problem show that the Beam Search technique is very efficient and appears to be a promising tool to solve difficult stochastic integer problems under probabilistic constraints. 2004 Elsevier B.V. All rights reserved. Keywords: Stochastic programming; Combinatorial optimization; Beam search heuristic
1. Introduction This paper addresses stochastic programming problems under probabilistic constraints where both the decision variables and the random variables are restricted to be integer. In particular, the class of problems considered throughout can be mathematically formulated as follows:
*
Corresponding author. Tel.: +39 0984 494826; fax: +39 0984 494713. E-mail address: [email protected] (P. Beraldi).
ðSIPCÞ
min
cT x Ax P b
ð1Þ ð2Þ
PfTx P ng P p x P 0 integer:
ð3Þ ð4Þ
Here T is a m · n integer matrix, A is a q · n matrix, c, x 2 Rn, b 2 Rq, n a m-dimensional random vector integer-valued and P denotes probability. Constraints (2) and (3) represent the deterministic and stochastic side of the problem, respectively. In particular, (3) are probabilistic constraints which ensure the satisfaction of the stochastic constraints Tx P n with a prescribed probability level
0377-2217/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2004.02.027
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
36
p 2 (0,1). In the following, we shall use Z and jÆj1 to denote the set of integers and the ‘1 norm, respectively. Furthermore, the inequality ‘‘ P ’’ for vectors will be understood coordinate-wise. SIPC problems arise in a broad spectrum of applications in which integer and combinatorial optimization problems are formulated in the presence of uncertainty and a ‘‘reliable solution’’ is required. Routing and location are two classical examples of applications where probabilistically constrained integer formulations have been applied to define robust strategic plans. The interested readers are referred to the survey papers of [9,13] and to the references therein. As one might expect, the difficulty in addressing SIPC problems is twofold. First, SIPC problems have a combinatorial nature. Second, (3) are joint probabilistic constraints and no assumption on the independence on the components ni, i = 1, . . ., m, of the random vector n, is imposed. To the best of our knowledge, all the probabilistic models proposed in the papers referred above either include individual chance constraints (i.e. probabilistic constraints individually imposed on all the inequalities) or assume independence of the components of the random vector. In both cases, equivalent deterministic formulations of the probabilistic constraints can be easily derived. In particular, the chance constraints PfT i x P ni g P pi ; i ¼ 1; . . . ; m can be replaced by the deterministic constraints: T ix P F i ðp i Þ;
i ¼ 1; . . . ; m;
i
where T denotes the ith row of the matrix T and F i is the pi-quantile of the distribution function of ni (see [17]). In the case of joint probabilistic constraints with independent random components, Eq. (3) can be replaced by m X lnðF i ðT i xÞÞ P ln p; i¼1
and specific reformulation can be then derived for the case of integer random variables as shown in [7]. The previous assumptions are rather restrictive and often not adequate to represent real-life applications. Nevertheless, the derivation of determinis-
tic equivalent formulations for joint probabilistic constraints requires the generation of the set of p-efficient points of the probability distribution function. The cardinality of this set is typically high even for moderate size of the random vector, amplifying the complexity associated with the solution of a single deterministic integer optimization problem. The inherent complexity of the SIPC problems poses the crucial problem to design efficient solution approaches. In this respect the literature is rather scarce. Dentcheva et al. proposed in [7] a cone generation method based on the convexification of a deterministic formulation of the problem. A Branch and Bound method embedding efficient strategies for the determination of lower and upper bound values was proposed in [3]. Both the approaches mentioned above belong to the class of exact solution methods. Because of the complexity of the SIPC problems, the applicability of these methods is typically limited to problems of small and medium sizes. As in deterministic setting, practical-sized instances may be solved rather effectively by resorting to heuristic solution approaches. This paper represents a contribution in this direction. In particular, we propose a heuristic solution approach based on a Beam Search strategy. This is an adaptation of the Branch and Bound method proposed in [3] in which at any level of the search tree only the promising nodes are kept for further exploration, whereas the remaining are pruned off permanently. The rest of the paper is organized as follows. In Section 2 we briefly recall the Branch and Bound method on which the heuristic approach is based on. Section 3 is devoted to the description of the Beam Search heuristic strategy for the solution of SIPC problems. The efficiency of the approach is tested on the probabilistic set covering problem formulated in [2]. The presentation and discussion of the numerical results is reported in Section 4. The paper ends with some concluding remarks.
2. The Branch and Bound method The Branch and Bound method for SIPC problems operates on a deterministic equivalent formu-
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
lation of the original problem which is obtained by means of the so called p-efficient points of the probability distribution function. In particular, given a probability level p 2 (0,1), a point s 2 Zm is called efficient point of the probability distribution function F, if F(s) P p. A point s is said to be p-efficient if it is efficient and satisfies the additional condition that there is no other point y 6 s, y „ s, such that F(y) P p. In [3], the authors propose an alternative characterization of the p-efficient points based on the conditional bounds, which are determined by the conditional marginal distribution functions, i.e. Fi(vijn 6 w) = P{ni 6 vijn 6 w}, for i = 1, . . ., m. The main properties of the conditional bounds, widely used within the Branch and Bound method, are reported in the following lemma (see [3] for the proof). Lemma 2.1. Let p 2 (0,1) and let w 2 Zm be p such that F(w) P p. Define li to be the F ðwÞ -efficient point of the conditional marginal Fi(vijn 6 w), that is li ðwÞ ¼ argminfji jF i ðji jn 6 wÞ P p=F ðwÞg; i ¼ 1; . . . ; m:
ð5Þ
Then (i) for every p-efficient point v 6 w we have v P l(w); (ii) if z 6 w then l(z) P l(w); (iii) w is p-efficient if and only if l(w) = w. The knowledge of the whole set J of the p-efficient points, allows to rewrite the probabilistic constraint Eq. (3) as Tx P sj
for at least one j 2 J :
It is worthwhile noting that the generation of the p-efficient points is computationally intensive and, in addition, the cardinality of the set J can be high even for moderate size of the random vector. This limits the applicability of solution approaches based on complete enumeration strategies only to instances of limited size. The Branch and Bound method relies on a partial enumeration of the p-efficient points which is
37
performed by taking advantage of the lower and upper bounds progressively determined during the computation. In the following, we shall briefly recall the Branch and Bound solution scheme. For a detailed description of the method, the reader is referred to Section 3 of [3]. The Branch and Bound method operates on a special tree where the nodes are candidate points and the arcs represent connections between nodes of consecutive levels. In particular, each node is an integer vector whose entries may take, for assumption, the values 0,1, . . ., hi (see [3]). Thus, each level j of the tree contains the vectors w for which jwj1 = j. For each node v at level j there are arcs (v,w) to all nodes w at level j + 1 which differ from v at exactly one coordinate. The root of the tree is the vector of highest level, whereas the leaves are either not efficient or p-efficient points. The Branch and Bound algorithm operates alternating the phases of enumeration of the efficient points and the solution of the corresponding integer programming problems. In particular, for each level k four main steps are performed: 1. 2. 3. 4.
generation, evaluation, solution, selection.
The first step consists in the generation of the candidate points. The nature of each candidate is established by computing the corresponding lower bound vector (see Lemma 2.1). For all the efficient points, we solve a relaxation of the resulting problem having the lower bound vector as right hand side. Finally, the selection of the promising candidates for further exploration is performed by comparing the value of the relaxation with the best solution found so far (i.e. incumbent value). During the computation, we maintain a list M, storing the p-efficient points. When no candidates for further exploration can be found, we pass to solve the problems in M. The solution process is carried out starting from the point with the lowest bound value and whenever the incumbent value is updated, other problems are removed from M until the list is empty.
38
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
Even with good lower bound and effective selection rules, the Branch and Bound method may require excessively large computation time. This is especially true when the solution time per problem is high. On the other hand, considering simple heuristic procedure to get an approximate solution of the problem (see, for example, those used in [3] to determine the initial incumbent value) may produce poor quality solution. A middle ground between the two extremes of exhaustive search and simple heuristics is to prune the search tree very harshly, allowing only a limited number of candidate nodes to be sprouted. The application of this technique to SIPC problems is described in Section 3.
the problem of finding a good trade off between quick, but poor evaluation and more computationally demanding, but better evaluation. The filtering mechanism represents a good solution to this problem. At first, a simple evaluation function is used to reduce the burden of the beam search. Only nodes not discarded during filtering are subject to a more detailed time consuming evaluation. The number of candidate nodes retained for further evaluation is called the filter width. Fig. 1 illustrates the Beam Search strategy. At the first level of the tree, we generate all the candidate nodes. On the basis of a first screening, we select among these, the most promising nodes to further explore. We observe that the first level of the tree in Fig. 1 actually corresponds to the (highest 1)th level according to the notation used in the generation of the p-efficient points.
3. The beam search strategy Beam search is a heuristic strategy closely related to the Branch and Bound method. It was first used in the field of artificial intelligence where it was applied to solve the speech recognition problem [10]. Successful applications of the Beam Search techniques can be also found in the area of scheduling problems (see for example [19] and the references therein). As the Branch and Bound scheme outlined in Section 2, Beam Search implements a breadth-first search. However, unlike the Branch and Bound, this heuristic technique reduces the width of the search, moving downward only from a limited number of best promising nodes. The number of selected nodes is called beam width. The success of the Beam Search strategy depends on the evaluation function which is used to select nodes that will be further explored. This is a much more crucial issue than in the Branch and Bound method since in an exhaustive search even with poor lower bound the optimal solution is bound to emerge. On the contrary an inadequate evaluation of the candidate nodes in a heuristic strategy may prune potentially promising candidates of the search tree forever, resulting in unsatisfactory solutions. As one might expect, a much careful evaluation requires a higher computational effort. This poses
Fig. 1. The beam search strategy.
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
In Section 3.1 we discuss the main issues in the application of the Beam Search strategy to solve SIPC problems.
3.1. Main issues The main ingredient of the Beam Search strategy is the evaluation function used to select the nodes to explore further at each level of the search tree. In the case of SIPC problems the evaluation function is represented by the lower bound value obtained by solving a relaxation of the problem corresponding to a candidate node. Let us denote by v the candidate node. A lower bound on the optimal solution of all the successors that could be generated starting from v is obtained by solving a relaxation of the problem: min
cT x Ax P b Tx P lðvÞ x P 0 integer;
ð6Þ
where the vector l(v) is computed by using the conditional marginal distribution function (see Lemma 2.1). We observe that the relaxed problems corresponding to the different candidate nodes differ in the right-hand side only, and, thus, in the case of linear relaxation, can be efficiently solved by a dual method exploiting, each time, the advanced basis of the previous optimization step. Even though efficient, the solution phase may require considerable effort because of the large number of candidate nodes and, thus, of problems to solve at each level of the search tree. This suggests the use a simpler evaluation function to use as filtering mechanism. Our algorithm uses a simple lower bound computed by the dual multipliers. In particular, we solve the linear relaxation of the ‘‘worst case’’ instance of the original problem: min
cT x Ax P b Tx P m x P 0 integer;
ð7Þ
where m is such that F(Tx) = F(m) = 1. By using the dual multipliers pi associated to constraints (7), it is possible to derive, for each candidate node v, a lower bound value
39
cðlðvÞÞ ¼ pT lðvÞ: This value is used to perform a first screening of the candidate nodes. More specifically, on the basis of the comparison of this weaker lower bound value with the best solution found so far, we select a first pool of promising candidate (equal to the filter width). Then, the filtered nodes are evaluated again by using the more detailed time consuming cost evaluation based on the solution of the linear relaxation. Among these, the best promising b candidates are selected to be further explored.
4. Computational experiments In this section we report on the computational experiments of the Beam Search algorithm. In order to test the heuristic method we have considered the probabilistic set covering problem proposed in [2]. The motivation in the choice of this problem is twofold. First, set covering is one of the most fundamental models of integer programming whose probabilistic counterpart can be used to model interesting applications in presence of uncertainty. Secondly, it is one of the few problems for which probabilistic instances have been defined. In particular, we have considered three sets of test problems. The first one consists of the large instances defined in [2]. More specifically, we use problems scp41 and scp42 from the Beasleys OR library [1] as deterministic side and Tests 11– 14 of [2] as stochastic side. We shall refer to these problems as Tests 1–8. By varying the deterministic side only, we have obtained other four instances (Tests 9–12). In particular, we have chosen problem scp52 from Beasleys OR library [1]. This is a deterministic set covering problem with m = 200 rows and n = 2000 columns, for which the linear relaxation does not provide an integer solution as for the instances scp41 and scp42. Finally, the third set is defined by instances of larger size. In particular, as deterministic side we have considered the 300 · 3000 set covering problem scpa3 of [1], whereas for the stochastic side we have used the circular and star dependency schemes described in [2]. We have divided the random vector n into l independent groups, each of size r, of
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
40
dependent random variables. By varying l and r, we have obtained four different test problems. Table 1 reports the type of dependency model, the groups size r, and the number of groups l. In the case of the circular dependency, we have assumed that the independent random variables follow a Bernoulli distribution, whereas in the case of star dependency we have considered a Poisson distribution. As far as the choice of the distribution parameters is concerned, for the first 200 components of the random vector we have considered the same parameters reported in Appendix of [1], whereas for the remaining 100 we have considered the parameters reported in Tables 5–8 of Appendix A. By combining the stochastic side and the deterministic one, we have obtained four different instances, refereed as Tests 13–16. All the instances have been solved for two levels of probability p = 0.95 and p = 0.90. All computational experiments have been carried out on a single processor R10K at 195 MHz
Table 1 Problems characteristics
Test Test Test Test
13 14 15 16
Model
Groups size
Number of groups
Star Circular Star Circular
5 5 10 10
60 60 30 30
of the supercomputer SGI Origin 2000, which has four nodes, with 4 Mb cache memory and 128 Mb RAM each. The codes have been implemented in the C language by using the state-ofthe-art LP solver CPLEX 6.5 [6]. The choice of this software, although probably not the most efficient for the solution of deterministic set covering instances, is motivated by its high flexibility. Indeed, it contains callable libraries which allow to perform several operations (such as sensitivity analysis, determination of the dual multipliers, etc.) used in our algorithmic schemes.
Table 2 Numerical results for the first set of test problems p
B&B Sol.
BS b = 5 Time
BS b = 10
BS b = 20
BS b = 30
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Test 1
0.95 0.90
419 413
5.67 60.80
419 419
3.80 22.20
419 419
4.10 27.40
419 417
4.45 38.88
419 414
4.45 45.67
Test 2
0.95 0.90
416 410
7.33 383.02
418 413
3.50 73.00
418 413
3.50 74.00
418 413
3.50 82.00
418 412
3.50 120.00
Test 3
0.95 0.90
416 402
2.99 33.87
416 411
1.80 19.14
416 411
1.92 20.30
416 411
1.92 21.40
416 410
2.10 23.12
Test 4
0.95 0.90
410 404
8.63 380.76
416 411
4.25 48.50
416 411
4.27 48.50
416 411
4.50 50.00
414 410
5.02 70.00
Test 5
0.95 0.90
510 484
7.68 55.61
510 484
4.36 20.82
510 484
4.36 20.82
510 484
4.36 20.82
510 484
4.36 20.82
Test 6
0.95 0.90
503 494
3.12 372.28
503 498
2.74 72.00
503 498
2.74 76.40
503 498
2.74 83.30
503 497
2.74 170.11
Test 7
0.95 0.90
487 447
8.34 62.58
488 479
1.40 28.92
488 463
1.50 33.47
488 460
1.67 43.95
488 460
2.50 59.01
Test 8
0.95 0.90
497 479
11.77 904.28
499 501
6.02 47.30
499 501
6.81 48.00
499 498
7.16 49.60
497 498
9.62 56.30
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
In order to evaluate the performance of our heuristic method, we have executed the algorithm for different beam widths. In particular, we have chosen four different values (5, 10, 20, 30). As far as filtering is concerned, we have adopted the mechanism described in Section 3.2, without fixing any limit as filter width. The numerical results are summarized in Tables 2–4 which report the value of the optimal solution (Sol.) and the computational time in seconds (Time) of both the Branch and Bound method and the Beam Search algorithm for different beam widths. We observe that the results of the Branch and Bound method refer to the hybrid strategy with dual heuristic (see [3] for details). The same strategy has been used in
41
the Beam Search algorithm to determine the initial upper bound value. We note that the running time required to solve the problems of the second set is considerable higher than that required for the problems in the first set. This behaviour is related to the considered problems. As already observed, Tests 9–12 have larger size (in terms of number of variables) than the problems of the first set. Furthermore, the linear relaxation of these problems does not provide an integer solution as for the problems of the first subset. Even higher running times are required to solve the problems in the third set.
4.1. Analysis and discussion of the numerical results
Table 3 Numerical results for the second set of test problems p
B&B
BS b = 5
BS b = 10
BS b = 20
BS b = 30
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Test 9
0.95 0.90
293 283
699.00 2033.00
294 291
300.52 739.01
294 291
353.00 993.12
293 291
474 1052.01
293 287
540.00 1394.00
Test 10
0.95 0.90
292 282
498.00 1842.00
295 291
270.12 412.03
292 290
373.00 522.10
292 290
383.5 814.00
292 288
390.99 598.05
Test 11
0.95 0.90
297 294
817.00 6354.02
299 298
466.17 568.00
297 298
487.70 927.00
297 298
493.90 1258.00
297 295
568.00 2145.12
Test 12
0.95 0.90
301 296
1025.00 1134.12
301 296
650.00 280.00
301 296
710.00 409.00
301 296
750.20 244.5
301 296
770.00 277.10
Table 4 Numerical results for the third set of test problems p
B&B Sol.
BS b = 5 Time
BS b = 10
BS b = 20
BS b = 30
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Test 13
0.95 0.90
224 213
4140.60 7967.60
225 223
3149.50 1367.70
225 223
3225.20 1628.80
225 223
3217.90 1977.30
224 222
3393.70 2771.50
Test 14
0.95 0.90
222 215
6715.10 19519.00
224 221
2888.90 4729.80
224 221
2927.50 4772.40
224 221
2872.90 5102.70
222 220
3354.80 9379.80
Test 15
0.95 0.90
224 214
9143.80 18757.20
226 225
3144.10 3751.20
226 225
3564.1 3917.00
226 225
2731.00 4181.00
226 223
3940.00 4693.30
Test 16
0.95 0.90
220 212
2820.40 7819.20
224 222
1201.40 2224.00
220 222
1694.20 2368.80
220 218
2269.10 3547.50
220 218
2126.00 4443.20
42
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
In evaluating the performance of a heuristic strategy two main issues have to be considered: the reduction of the solution time with the respect to the exact method and the quality of the solution. We analyze the time factor first. It is worthwhile noting that the overall solution time depends on the number of explored points, and, thus, on the time needed for their generation and evaluation. Our approach is heuristic in that it explores a subset of points generated by the Branch and Bound method. For a given test problem, the number of points explored during the computation depends on both the selected probability level p and the beam with b. Figs. 2 and 3
show the time reduction (in percentage) versus the beam width for the test problems belonging to the second set, for p equal to 0.95 and 0.90, respectively. As expected, the smaller is the beam width the higher is the saving in the solution time. We observe, however, that the reduction in the solution time does not decrease linearly with the beam width. For example, for Test 10 (see Fig. 3) the saving in the solution time is about 56% for b = 20 and 67% for b = 30. This phenomena can be explained by observing that different beam widths lead to the selection of different set of candidate nodes. Although the number of nodes to be explored is larger in the the case of larger beam
Fig. 2. Time reduction versus beam width for test problems of the second subset with p = 0.95.
Fig. 3. Time reduction versus beam width for test problems of the second subset with p = 0.90.
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
widths, an early update of the incumbent value may occur, resulting in a reduction of the number of candidate nodes. By the comparison of Figs. 2 and 3, it comes out that higher time reduction can be obtained for the instances solved with a probability level p of 0.90. Obviously, when a lower probability level is imposed a higher number of points have to be processed, and, thus, problems have to be solved. It is worthwhile noting that larger instances are more computationally demanding, since both the generation and the solution phases are much time consuming. Tables 9 and 10 in Appendix B show the time reduction for the three sets of problems, for both probability levels. We note that the problems belonging to the first set are divided into two subsets, namely Sets 1A and 1B. Fig. 4 summarizes the reported results showing the average time reduction as function of the beam width for p = 0.90. We note that typically for a fixed beam width, higher savings can be obtained for larger instances. By analyzing Fig. 4, we observe, however, that for b = 20, time reduction for problems belonging to the first set is higher than that of the third set. This behaviour can be explained by observing that because of the specific problems structure, it may occur that the optimal solution is determined during the computation, thus, reducing the number of explored points. In any case, we expect that better performance can be achieved when larger test problems are considered. Resorting to heuristic approaches is in
43
general inevitable for solving practical-sized instances of most integer programming problems. This requirement is even stronger in our case, since the complexity related to the stochastic nature is added to that related to the solution of integer problems. The other issue to consider in evaluating the efficiency of a heuristic approach is the solutions quality. As one may expect, wider beams lead to better solutions, because the pool of candidate nodes is larger. The numerical results show that the Beam Search strategy exhibits very good performance at least on the test problems considered here. As matter of fact, we observe that even for a beam width of 5 the maximum error (recorded only for one instance) of the heuristic solution with the respect to optimal solution found by the Branch and Bound scheme is around 5%. Furthermore, in some cases the Beam Search algorithm finds the optimal solution. The two main issues discussed above determine the best beam width value. In setting this value, one has to consider the trade-off between the solutions quality and the savings in the solution time. The computational experiments show that the choice of good beam widths depends on the considered problems and requires intensive computational testing and fine tuning.
Fig. 4. Average time reduction versus beam width for the three sets of test problems with p = 0.90.
44
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
The procedure terminates when the current solution may not be improved. Several variants of this basic scheme could be considered (i.e. Tabu´ search, Simulated Annealing). Finally, we observe that in solving stochastic integer problems under probabilistic constraints two levels of approximation may be considered: the first one relies on partially exploring the space of efficient points and second one on heuristically solving the corresponding problems. The design of a two-level heuristic strategy is the subject of ongoing research.
5. Conclusion and future research directions In this paper we have proposed a heuristic strategy to solve stochastic integer problems under probabilistic constraints. The heuristic method implements a beam search strategy which is based on the classical Branch and Bound scheme. At any level of the search tree, only the promising nodes are kept for further exploration, whereas the remaining are pruned out permanently. The computational experiments show that the solution time of the Beam Search is far less than the time required by the Branch and Bound method. Furthermore, the Beam Search algorithm turns out to provide a high quality solution even for low beam width. Hence, the algorithm appears to be efficient and particularly suitable to solve the class of stochastic integer problems under probabilistic constraints. A future research direction is the design and implementation of other heuristic strategies. A first alternative is represented by the classical local search applied in the space of efficient points. In particular, starting from a p-efficient point we generate its neighbourhood, and we move to another point only if an evaluation function is improved.
6. Uncited references [4,5,8,11,12,14–16,18,20]
Appendix A Tables 5–8 report the values of the distribution parameters ai, ki of the circular and star dependency schemes for the additional 100 components of the test problem Tests 13–16. The first 200 components are reported in the tables of Appendix of the paper [2]. The tests have been defined in such a
Table 5 Distribution parameters of Test 13 Group
k1
k2
k3
k4
k5
k6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.790 0.810 0.200 0.810 0.200 0.010 0.809 0.810 0.018 0.601 0.018 0.011 0.011 0.011 0.005 0.090 0.710 0.020 0.011 0.005
0.005 0.010 0.090 0.010 0.890 0.810 0.980 0.810 0.011 0.011 0.811 0.011 0.911 0.011 0.005 0.020 0.140 0.020 0.011 0.005
0.005 0.100 0.001 0.100 0.010 0.900 0.710 0.900 0.011 0.041 0.011 0.089 0.011 0.089 0.805 0.910 0.150 0.810 0.089 0.805
0.990 0.051 0.700 0.951 0.700 0.810 0.020 0.810 0.011 0.611 0.011 0.770 0.811 0.700 0.090 0.050 0.150 0.910 0.700 0.090
0.007 0.010 0.300 0.010 0.300 0.850 0.061 0.005 0.911 0.971 0.011 0.009 0.811 0.009 0.007 0.010 0.010 0.010 0.009 0.007
0.907 0.810 0.300 0.810 0.010 0.101 0.061 0.001 0.615 0.911 0.615 0.011 0.189 0.011 0.900 0.710 0.710 0.010 0.011 0.210
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
45
Table 6 Distribution parameters of Test 14 Group
a1
a2
a3
a4
a5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.010 0.001 0.010 0.500 0.090 0.990 0.010 0.050 0.010 0.010 0.010 0.820 0.030 0.051 0.150 0.500 0.010 0.810 0.030 0.951
0.410 0.100 0.020 0.010 0.900 0.007 0.810 0.050 0.100 0.010 0.410 0.200 0.020 0.410 0.550 0.200 0.410 0.200 0.030 0.010
0.420 0.002 0.030 0.210 0.890 0.907 0.090 0.200 0.002 0.050 0.420 0.090 0.810 0.010 0.090 0.100 0.420 0.090 0.810 0.810
0.051 0.300 0.110 0.010 0.050 0.900 0.090 0.070 0.300 0.010 0.051 0.810 0.010 0.200 0.070 0.110 0.051 0.810 0.010 0.200
0.010 0.500 0.110 0.100 0.030 0.910 0.190 0.900 0.500 0.040 0.010 0.700 0.100 0.090 0.100 0.810 0.010 0.020 0.100 0.890
Table 7 Distribution parameters of Test 15 Group
k1
k2
k3
k4
k5
k6
k7
k8
k9
k10
k11
1 2 3 4 5 6 7 8 9 10
0.890 0.090 0.820 0.710 0.011 0.189 0.911 0.911 0.601 0.011
0.890 0.090 0.061 0.020 0.011 0.011 0.189 0.841 0.011 0.011
0.005 0.790 0.061 0.061 0.089 0.011 0.010 0.611 0.041 0.089
0.005 0.005 0.810 0.961 0.770 0.089 0.710 0.871 0.611 0.770
0.990 0.005 0.810 0.100 0.009 0.770 0.700 0.911 0.971 0.009
0.007 0.990 0.900 0.010 0.011 0.009 0.810 0.018 0.911 0.011
0.907 0.007 0.810 0.100 0.011 0.011 0.050 0.711 0.018 0.011
0.090 0.907 0.050 0.020 0.911 0.011 0.101 0.011 0.811 0.911
0.910 0.900 0.801 0.800 0.011 0.911 0.909 0.011 0.011 0.011
0.010 0.010 0.809 0.500 0.811 0.911 0.890 0.011 0.110 0.811
0.810 0.810 0.890 0.200 0.811 0.911 0.020 0.615 0.615 0.811
Table 8 Distribution parameters of Test 16 Group
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
1 2 3 4 5 6 7 8 9 10
0.100 0.500 0.090 0.910 0.810 0.900 0.810 0.809 0.011 0.099
0.110 0.200 0.890 0.010 0.910 0.300 0.700 0.990 0.911 0.011
0.810 0.100 0.890 0.810 0.100 0.810 0.300 0.710 0.011 0.011
0.810 0.110 0.050 0.900 0.051 0.810 0.910 0.820 0.811 0.911
0.100 0.810 0.005 0.090 0.010 0.100 0.910 0.061 0.811 0.911
0.010 0.810 0.005 0.790 0.810 0.951 0.8100 0.061 0.189 0.911
0.100 0.810 0.990 0.005 0.200 0.010 0.900 0.810 0.011 0.911
0.020 0.910 0.007 0.005 0.090 0.810 0.851 0.810 0.011 0.901
0.300 0.001 0.907 0.990 0.010 0.200 0.850 0.900 0.089 0.189
0.500 0.810 0.900 0.007 0.700 0.890 0.101 0.810 0.770 0.010
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
46
way to reflect the freedom in the users choice; thus, different groups of dependent random variables may have different values of the distribution parameters selected from the range. For each dependency scheme, the test problems are grouped in Tables according to the size r of each group. We note that in the case of star dependency (Tests 13 and 15), we have included an additional random
variable with parameter kr + 1, denoting the center of the star.
Appendix B The following Tables 9 and 10 report the time reduction (in percentage) for the three sets of prob-
Table 9 Time reduction (in percentage) for p = 0.95 Test problem
b=5
b = 10
b = 20
b = 30
Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test
32.98 52.25 39.80 50.75 43.23 12.18 83.21 48.85 57.01 45.76 42.94 36.59 23.94 56.98 65.61 57.40
27.69 52.25 35.79 50.52 43.23 12.18 82.01 42.14 49.50 25.10 40.31 30.73 22.10 56.40 61.02 39.93
21.52 52.25 35.79 97.67 43.23 12.18 79.98 39.17 32.19 22.99 39.55 26.83 22.28 57.22 70.13 19.55
21.52 52.25 29.77 42.06 43.23 12.18 70.02 18.27 22.75 21.49 30.48 24.88 18.04 50.04 56.91 24.62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Table 10 Time reduction (in percentage) for p = 0.90 Test problem
b=5
b = 10
b = 20
b = 30
Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test
63.49 80.94 43.49 87.26 62.56 80.66 53.79 94.77 63.65 77.63 91.06 75.31 82.83 75.77 80.00 71.56
54.93 80.68 40.06 87.26 62.56 79.48 46.52 94.69 51.15 71.66 85.41 63.94 79.56 75.55 79.12 69.71
36.05 78.59 36.82 86.87 92.16 77.62 29.77 94.51 48.25 55.81 80.20 78.44 75.18 73.86 77.71 54.63
24.88 68.67 31.74 81.62 62.56 54.31 5.70 93.77 31.43 67.54 66.24 75.57 65.22 51.95 74.98 43.18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
lems. We observe that the test problems of the first set have been divided into two subsets: Sets 1A and 1B. References [1] J.E. Beasley, OR-library: Distributing test problems by electronic mail, Journal of the Operational Research Society 41 (1990) 1069–1072. [2] P. Beraldi, A. Ruszczyn´ski, The probabilistic set covering problem, Operations Research 50 (2002) 956–967. [3] P. Beraldi, A. Ruszczyn´ski, A Branch and Bound method for stochastic integer problems under probabilistic constraints, Optimization Methods and Software 17 (2002) 359–382. [4] A. Charnes, W.W. Cooper, Chance-constrained programming, Management Science 5 (1958) 73–79. [5] A. Charnes, W.W. Cooper, G.H. Symonds, Cost horizons and certainty equivalents: An approach to stochastic programming of heating oil, Management Science 4 (1958) 235–263. [6] CPLEX, ILOG CPLEX 6.5: Users Manual, CPLEX Optimization, Inc., Incline Village, NV, 1999. [7] D. Dentcheva, A. Pre´kopa, A. Ruszczyn´ski, Bounds for probabilistic integer programming problems, Discrete Applied Mathematics 124 (2002) 55–65. [8] D. Dentcheva, A. Pre´kopa, A. Ruszczyn´ski, Concavity and efficient points of discrete distributions in probabilistic programming, Mathematical Programming 89 (2001) 217–232. [9] M. Gendreau, G. Laporte, R. Se´guin, Stochastic vehicle routing, European Journal of Operational Research 88 (1996) 3–12.
47
[10] B.T. Lowerre, The HARPY speech recognition system, Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA. [11] L.B. Miller, H. Wagner, Change-constrained programming with joint constraints, Operations Research 13 (1965) 930– 945. [12] G.L. Nemhauser, L.A. Wolsey, Integer and Combinatorial Optimization, John Wiley & Sons, New York, 1988. [13] S.H. Owen, M.S. Daskin, Strategic facility location: A review, European Journal of Operational Research 111 (1998) 423–447. [14] A. Pre´kopa, On probabilistic constrained programming, in: Proceedings of the Princeton Symposium on Mathematical Programming, Princeton University Press, 1970, pp. 113– 138. [15] A. Pre´kopa, Contributions to the theory of stochastic programming, Mathematical Programming 4 (1973) 202– 221. [16] A. Pre´kopa, Dual method for solution of one-stage stochastic programming with random RHS obeying a discrete probability distribution, Zeitschrift fu¨r Operations Research 34 (1990) 441–461. [17] A. Pre´kopa, Stochastic Programming, Kluwer Academic Publishers, Boston, MA, 1995. [18] A. Pre´kopa, D. Vizva´ri, T. Badics, Programming under probabilistic constraint with discrete random variable, in: F. Giannesi (Ed.), New Trends in Mathematical Programming, Kluwer Academic Publishers, 1998, pp. 235–255. [19] I. Sabuncuoglu, M. Bayiz, Job shop scheduling with beam search, European Journal of Operational Research 118 (1999) 390–412. [20] S. Sen, Relaxations for the probabilistically constrained programs with discrete random variables, Operations Research Letters 11 (1992) 81–86.
Discrete Optimization
Beam search heuristic to solve stochastic integer problems under probabilistic constraints Patrizia Beraldi b
a,*
, Andrzej Ruszczyn´ski
b
a Dipartimento di Elettronica, Informatica e Sistemistica, Universita` della Calabria, 87036 Rende (CS), Italy RUTCOR––Rutgers Center for Operations Research, 640 Bartholomew Road, Piscataway, NJ 08854-8003, USA
Received 27 May 2002; accepted 23 February 2004 Available online 26 June 2004
Abstract This paper proposes a Beam Search heuristic strategy to solve stochastic integer programming problems under probabilistic constraints. Beam Search is an adaptation of the classical Branch and Bound method in which at any level of the search tree only the most promising nodes are kept for further exploration, whereas the remaining are pruned out permanently. The proposed algorithm has been compared with the Branch and Bound method. The numerical results collected on the probabilistic set covering problem show that the Beam Search technique is very efficient and appears to be a promising tool to solve difficult stochastic integer problems under probabilistic constraints. 2004 Elsevier B.V. All rights reserved. Keywords: Stochastic programming; Combinatorial optimization; Beam search heuristic
1. Introduction This paper addresses stochastic programming problems under probabilistic constraints where both the decision variables and the random variables are restricted to be integer. In particular, the class of problems considered throughout can be mathematically formulated as follows:
*
Corresponding author. Tel.: +39 0984 494826; fax: +39 0984 494713. E-mail address: [email protected] (P. Beraldi).
ðSIPCÞ
min
cT x Ax P b
ð1Þ ð2Þ
PfTx P ng P p x P 0 integer:
ð3Þ ð4Þ
Here T is a m · n integer matrix, A is a q · n matrix, c, x 2 Rn, b 2 Rq, n a m-dimensional random vector integer-valued and P denotes probability. Constraints (2) and (3) represent the deterministic and stochastic side of the problem, respectively. In particular, (3) are probabilistic constraints which ensure the satisfaction of the stochastic constraints Tx P n with a prescribed probability level
0377-2217/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2004.02.027
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
36
p 2 (0,1). In the following, we shall use Z and jÆj1 to denote the set of integers and the ‘1 norm, respectively. Furthermore, the inequality ‘‘ P ’’ for vectors will be understood coordinate-wise. SIPC problems arise in a broad spectrum of applications in which integer and combinatorial optimization problems are formulated in the presence of uncertainty and a ‘‘reliable solution’’ is required. Routing and location are two classical examples of applications where probabilistically constrained integer formulations have been applied to define robust strategic plans. The interested readers are referred to the survey papers of [9,13] and to the references therein. As one might expect, the difficulty in addressing SIPC problems is twofold. First, SIPC problems have a combinatorial nature. Second, (3) are joint probabilistic constraints and no assumption on the independence on the components ni, i = 1, . . ., m, of the random vector n, is imposed. To the best of our knowledge, all the probabilistic models proposed in the papers referred above either include individual chance constraints (i.e. probabilistic constraints individually imposed on all the inequalities) or assume independence of the components of the random vector. In both cases, equivalent deterministic formulations of the probabilistic constraints can be easily derived. In particular, the chance constraints PfT i x P ni g P pi ; i ¼ 1; . . . ; m can be replaced by the deterministic constraints: T ix P F i ðp i Þ;
i ¼ 1; . . . ; m;
i
where T denotes the ith row of the matrix T and F i is the pi-quantile of the distribution function of ni (see [17]). In the case of joint probabilistic constraints with independent random components, Eq. (3) can be replaced by m X lnðF i ðT i xÞÞ P ln p; i¼1
and specific reformulation can be then derived for the case of integer random variables as shown in [7]. The previous assumptions are rather restrictive and often not adequate to represent real-life applications. Nevertheless, the derivation of determinis-
tic equivalent formulations for joint probabilistic constraints requires the generation of the set of p-efficient points of the probability distribution function. The cardinality of this set is typically high even for moderate size of the random vector, amplifying the complexity associated with the solution of a single deterministic integer optimization problem. The inherent complexity of the SIPC problems poses the crucial problem to design efficient solution approaches. In this respect the literature is rather scarce. Dentcheva et al. proposed in [7] a cone generation method based on the convexification of a deterministic formulation of the problem. A Branch and Bound method embedding efficient strategies for the determination of lower and upper bound values was proposed in [3]. Both the approaches mentioned above belong to the class of exact solution methods. Because of the complexity of the SIPC problems, the applicability of these methods is typically limited to problems of small and medium sizes. As in deterministic setting, practical-sized instances may be solved rather effectively by resorting to heuristic solution approaches. This paper represents a contribution in this direction. In particular, we propose a heuristic solution approach based on a Beam Search strategy. This is an adaptation of the Branch and Bound method proposed in [3] in which at any level of the search tree only the promising nodes are kept for further exploration, whereas the remaining are pruned off permanently. The rest of the paper is organized as follows. In Section 2 we briefly recall the Branch and Bound method on which the heuristic approach is based on. Section 3 is devoted to the description of the Beam Search heuristic strategy for the solution of SIPC problems. The efficiency of the approach is tested on the probabilistic set covering problem formulated in [2]. The presentation and discussion of the numerical results is reported in Section 4. The paper ends with some concluding remarks.
2. The Branch and Bound method The Branch and Bound method for SIPC problems operates on a deterministic equivalent formu-
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
lation of the original problem which is obtained by means of the so called p-efficient points of the probability distribution function. In particular, given a probability level p 2 (0,1), a point s 2 Zm is called efficient point of the probability distribution function F, if F(s) P p. A point s is said to be p-efficient if it is efficient and satisfies the additional condition that there is no other point y 6 s, y „ s, such that F(y) P p. In [3], the authors propose an alternative characterization of the p-efficient points based on the conditional bounds, which are determined by the conditional marginal distribution functions, i.e. Fi(vijn 6 w) = P{ni 6 vijn 6 w}, for i = 1, . . ., m. The main properties of the conditional bounds, widely used within the Branch and Bound method, are reported in the following lemma (see [3] for the proof). Lemma 2.1. Let p 2 (0,1) and let w 2 Zm be p such that F(w) P p. Define li to be the F ðwÞ -efficient point of the conditional marginal Fi(vijn 6 w), that is li ðwÞ ¼ argminfji jF i ðji jn 6 wÞ P p=F ðwÞg; i ¼ 1; . . . ; m:
ð5Þ
Then (i) for every p-efficient point v 6 w we have v P l(w); (ii) if z 6 w then l(z) P l(w); (iii) w is p-efficient if and only if l(w) = w. The knowledge of the whole set J of the p-efficient points, allows to rewrite the probabilistic constraint Eq. (3) as Tx P sj
for at least one j 2 J :
It is worthwhile noting that the generation of the p-efficient points is computationally intensive and, in addition, the cardinality of the set J can be high even for moderate size of the random vector. This limits the applicability of solution approaches based on complete enumeration strategies only to instances of limited size. The Branch and Bound method relies on a partial enumeration of the p-efficient points which is
37
performed by taking advantage of the lower and upper bounds progressively determined during the computation. In the following, we shall briefly recall the Branch and Bound solution scheme. For a detailed description of the method, the reader is referred to Section 3 of [3]. The Branch and Bound method operates on a special tree where the nodes are candidate points and the arcs represent connections between nodes of consecutive levels. In particular, each node is an integer vector whose entries may take, for assumption, the values 0,1, . . ., hi (see [3]). Thus, each level j of the tree contains the vectors w for which jwj1 = j. For each node v at level j there are arcs (v,w) to all nodes w at level j + 1 which differ from v at exactly one coordinate. The root of the tree is the vector of highest level, whereas the leaves are either not efficient or p-efficient points. The Branch and Bound algorithm operates alternating the phases of enumeration of the efficient points and the solution of the corresponding integer programming problems. In particular, for each level k four main steps are performed: 1. 2. 3. 4.
generation, evaluation, solution, selection.
The first step consists in the generation of the candidate points. The nature of each candidate is established by computing the corresponding lower bound vector (see Lemma 2.1). For all the efficient points, we solve a relaxation of the resulting problem having the lower bound vector as right hand side. Finally, the selection of the promising candidates for further exploration is performed by comparing the value of the relaxation with the best solution found so far (i.e. incumbent value). During the computation, we maintain a list M, storing the p-efficient points. When no candidates for further exploration can be found, we pass to solve the problems in M. The solution process is carried out starting from the point with the lowest bound value and whenever the incumbent value is updated, other problems are removed from M until the list is empty.
38
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
Even with good lower bound and effective selection rules, the Branch and Bound method may require excessively large computation time. This is especially true when the solution time per problem is high. On the other hand, considering simple heuristic procedure to get an approximate solution of the problem (see, for example, those used in [3] to determine the initial incumbent value) may produce poor quality solution. A middle ground between the two extremes of exhaustive search and simple heuristics is to prune the search tree very harshly, allowing only a limited number of candidate nodes to be sprouted. The application of this technique to SIPC problems is described in Section 3.
the problem of finding a good trade off between quick, but poor evaluation and more computationally demanding, but better evaluation. The filtering mechanism represents a good solution to this problem. At first, a simple evaluation function is used to reduce the burden of the beam search. Only nodes not discarded during filtering are subject to a more detailed time consuming evaluation. The number of candidate nodes retained for further evaluation is called the filter width. Fig. 1 illustrates the Beam Search strategy. At the first level of the tree, we generate all the candidate nodes. On the basis of a first screening, we select among these, the most promising nodes to further explore. We observe that the first level of the tree in Fig. 1 actually corresponds to the (highest 1)th level according to the notation used in the generation of the p-efficient points.
3. The beam search strategy Beam search is a heuristic strategy closely related to the Branch and Bound method. It was first used in the field of artificial intelligence where it was applied to solve the speech recognition problem [10]. Successful applications of the Beam Search techniques can be also found in the area of scheduling problems (see for example [19] and the references therein). As the Branch and Bound scheme outlined in Section 2, Beam Search implements a breadth-first search. However, unlike the Branch and Bound, this heuristic technique reduces the width of the search, moving downward only from a limited number of best promising nodes. The number of selected nodes is called beam width. The success of the Beam Search strategy depends on the evaluation function which is used to select nodes that will be further explored. This is a much more crucial issue than in the Branch and Bound method since in an exhaustive search even with poor lower bound the optimal solution is bound to emerge. On the contrary an inadequate evaluation of the candidate nodes in a heuristic strategy may prune potentially promising candidates of the search tree forever, resulting in unsatisfactory solutions. As one might expect, a much careful evaluation requires a higher computational effort. This poses
Fig. 1. The beam search strategy.
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
In Section 3.1 we discuss the main issues in the application of the Beam Search strategy to solve SIPC problems.
3.1. Main issues The main ingredient of the Beam Search strategy is the evaluation function used to select the nodes to explore further at each level of the search tree. In the case of SIPC problems the evaluation function is represented by the lower bound value obtained by solving a relaxation of the problem corresponding to a candidate node. Let us denote by v the candidate node. A lower bound on the optimal solution of all the successors that could be generated starting from v is obtained by solving a relaxation of the problem: min
cT x Ax P b Tx P lðvÞ x P 0 integer;
ð6Þ
where the vector l(v) is computed by using the conditional marginal distribution function (see Lemma 2.1). We observe that the relaxed problems corresponding to the different candidate nodes differ in the right-hand side only, and, thus, in the case of linear relaxation, can be efficiently solved by a dual method exploiting, each time, the advanced basis of the previous optimization step. Even though efficient, the solution phase may require considerable effort because of the large number of candidate nodes and, thus, of problems to solve at each level of the search tree. This suggests the use a simpler evaluation function to use as filtering mechanism. Our algorithm uses a simple lower bound computed by the dual multipliers. In particular, we solve the linear relaxation of the ‘‘worst case’’ instance of the original problem: min
cT x Ax P b Tx P m x P 0 integer;
ð7Þ
where m is such that F(Tx) = F(m) = 1. By using the dual multipliers pi associated to constraints (7), it is possible to derive, for each candidate node v, a lower bound value
39
cðlðvÞÞ ¼ pT lðvÞ: This value is used to perform a first screening of the candidate nodes. More specifically, on the basis of the comparison of this weaker lower bound value with the best solution found so far, we select a first pool of promising candidate (equal to the filter width). Then, the filtered nodes are evaluated again by using the more detailed time consuming cost evaluation based on the solution of the linear relaxation. Among these, the best promising b candidates are selected to be further explored.
4. Computational experiments In this section we report on the computational experiments of the Beam Search algorithm. In order to test the heuristic method we have considered the probabilistic set covering problem proposed in [2]. The motivation in the choice of this problem is twofold. First, set covering is one of the most fundamental models of integer programming whose probabilistic counterpart can be used to model interesting applications in presence of uncertainty. Secondly, it is one of the few problems for which probabilistic instances have been defined. In particular, we have considered three sets of test problems. The first one consists of the large instances defined in [2]. More specifically, we use problems scp41 and scp42 from the Beasleys OR library [1] as deterministic side and Tests 11– 14 of [2] as stochastic side. We shall refer to these problems as Tests 1–8. By varying the deterministic side only, we have obtained other four instances (Tests 9–12). In particular, we have chosen problem scp52 from Beasleys OR library [1]. This is a deterministic set covering problem with m = 200 rows and n = 2000 columns, for which the linear relaxation does not provide an integer solution as for the instances scp41 and scp42. Finally, the third set is defined by instances of larger size. In particular, as deterministic side we have considered the 300 · 3000 set covering problem scpa3 of [1], whereas for the stochastic side we have used the circular and star dependency schemes described in [2]. We have divided the random vector n into l independent groups, each of size r, of
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
40
dependent random variables. By varying l and r, we have obtained four different test problems. Table 1 reports the type of dependency model, the groups size r, and the number of groups l. In the case of the circular dependency, we have assumed that the independent random variables follow a Bernoulli distribution, whereas in the case of star dependency we have considered a Poisson distribution. As far as the choice of the distribution parameters is concerned, for the first 200 components of the random vector we have considered the same parameters reported in Appendix of [1], whereas for the remaining 100 we have considered the parameters reported in Tables 5–8 of Appendix A. By combining the stochastic side and the deterministic one, we have obtained four different instances, refereed as Tests 13–16. All the instances have been solved for two levels of probability p = 0.95 and p = 0.90. All computational experiments have been carried out on a single processor R10K at 195 MHz
Table 1 Problems characteristics
Test Test Test Test
13 14 15 16
Model
Groups size
Number of groups
Star Circular Star Circular
5 5 10 10
60 60 30 30
of the supercomputer SGI Origin 2000, which has four nodes, with 4 Mb cache memory and 128 Mb RAM each. The codes have been implemented in the C language by using the state-ofthe-art LP solver CPLEX 6.5 [6]. The choice of this software, although probably not the most efficient for the solution of deterministic set covering instances, is motivated by its high flexibility. Indeed, it contains callable libraries which allow to perform several operations (such as sensitivity analysis, determination of the dual multipliers, etc.) used in our algorithmic schemes.
Table 2 Numerical results for the first set of test problems p
B&B Sol.
BS b = 5 Time
BS b = 10
BS b = 20
BS b = 30
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Test 1
0.95 0.90
419 413
5.67 60.80
419 419
3.80 22.20
419 419
4.10 27.40
419 417
4.45 38.88
419 414
4.45 45.67
Test 2
0.95 0.90
416 410
7.33 383.02
418 413
3.50 73.00
418 413
3.50 74.00
418 413
3.50 82.00
418 412
3.50 120.00
Test 3
0.95 0.90
416 402
2.99 33.87
416 411
1.80 19.14
416 411
1.92 20.30
416 411
1.92 21.40
416 410
2.10 23.12
Test 4
0.95 0.90
410 404
8.63 380.76
416 411
4.25 48.50
416 411
4.27 48.50
416 411
4.50 50.00
414 410
5.02 70.00
Test 5
0.95 0.90
510 484
7.68 55.61
510 484
4.36 20.82
510 484
4.36 20.82
510 484
4.36 20.82
510 484
4.36 20.82
Test 6
0.95 0.90
503 494
3.12 372.28
503 498
2.74 72.00
503 498
2.74 76.40
503 498
2.74 83.30
503 497
2.74 170.11
Test 7
0.95 0.90
487 447
8.34 62.58
488 479
1.40 28.92
488 463
1.50 33.47
488 460
1.67 43.95
488 460
2.50 59.01
Test 8
0.95 0.90
497 479
11.77 904.28
499 501
6.02 47.30
499 501
6.81 48.00
499 498
7.16 49.60
497 498
9.62 56.30
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
In order to evaluate the performance of our heuristic method, we have executed the algorithm for different beam widths. In particular, we have chosen four different values (5, 10, 20, 30). As far as filtering is concerned, we have adopted the mechanism described in Section 3.2, without fixing any limit as filter width. The numerical results are summarized in Tables 2–4 which report the value of the optimal solution (Sol.) and the computational time in seconds (Time) of both the Branch and Bound method and the Beam Search algorithm for different beam widths. We observe that the results of the Branch and Bound method refer to the hybrid strategy with dual heuristic (see [3] for details). The same strategy has been used in
41
the Beam Search algorithm to determine the initial upper bound value. We note that the running time required to solve the problems of the second set is considerable higher than that required for the problems in the first set. This behaviour is related to the considered problems. As already observed, Tests 9–12 have larger size (in terms of number of variables) than the problems of the first set. Furthermore, the linear relaxation of these problems does not provide an integer solution as for the problems of the first subset. Even higher running times are required to solve the problems in the third set.
4.1. Analysis and discussion of the numerical results
Table 3 Numerical results for the second set of test problems p
B&B
BS b = 5
BS b = 10
BS b = 20
BS b = 30
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Test 9
0.95 0.90
293 283
699.00 2033.00
294 291
300.52 739.01
294 291
353.00 993.12
293 291
474 1052.01
293 287
540.00 1394.00
Test 10
0.95 0.90
292 282
498.00 1842.00
295 291
270.12 412.03
292 290
373.00 522.10
292 290
383.5 814.00
292 288
390.99 598.05
Test 11
0.95 0.90
297 294
817.00 6354.02
299 298
466.17 568.00
297 298
487.70 927.00
297 298
493.90 1258.00
297 295
568.00 2145.12
Test 12
0.95 0.90
301 296
1025.00 1134.12
301 296
650.00 280.00
301 296
710.00 409.00
301 296
750.20 244.5
301 296
770.00 277.10
Table 4 Numerical results for the third set of test problems p
B&B Sol.
BS b = 5 Time
BS b = 10
BS b = 20
BS b = 30
Sol.
Time
Sol.
Time
Sol.
Time
Sol.
Time
Test 13
0.95 0.90
224 213
4140.60 7967.60
225 223
3149.50 1367.70
225 223
3225.20 1628.80
225 223
3217.90 1977.30
224 222
3393.70 2771.50
Test 14
0.95 0.90
222 215
6715.10 19519.00
224 221
2888.90 4729.80
224 221
2927.50 4772.40
224 221
2872.90 5102.70
222 220
3354.80 9379.80
Test 15
0.95 0.90
224 214
9143.80 18757.20
226 225
3144.10 3751.20
226 225
3564.1 3917.00
226 225
2731.00 4181.00
226 223
3940.00 4693.30
Test 16
0.95 0.90
220 212
2820.40 7819.20
224 222
1201.40 2224.00
220 222
1694.20 2368.80
220 218
2269.10 3547.50
220 218
2126.00 4443.20
42
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
In evaluating the performance of a heuristic strategy two main issues have to be considered: the reduction of the solution time with the respect to the exact method and the quality of the solution. We analyze the time factor first. It is worthwhile noting that the overall solution time depends on the number of explored points, and, thus, on the time needed for their generation and evaluation. Our approach is heuristic in that it explores a subset of points generated by the Branch and Bound method. For a given test problem, the number of points explored during the computation depends on both the selected probability level p and the beam with b. Figs. 2 and 3
show the time reduction (in percentage) versus the beam width for the test problems belonging to the second set, for p equal to 0.95 and 0.90, respectively. As expected, the smaller is the beam width the higher is the saving in the solution time. We observe, however, that the reduction in the solution time does not decrease linearly with the beam width. For example, for Test 10 (see Fig. 3) the saving in the solution time is about 56% for b = 20 and 67% for b = 30. This phenomena can be explained by observing that different beam widths lead to the selection of different set of candidate nodes. Although the number of nodes to be explored is larger in the the case of larger beam
Fig. 2. Time reduction versus beam width for test problems of the second subset with p = 0.95.
Fig. 3. Time reduction versus beam width for test problems of the second subset with p = 0.90.
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
widths, an early update of the incumbent value may occur, resulting in a reduction of the number of candidate nodes. By the comparison of Figs. 2 and 3, it comes out that higher time reduction can be obtained for the instances solved with a probability level p of 0.90. Obviously, when a lower probability level is imposed a higher number of points have to be processed, and, thus, problems have to be solved. It is worthwhile noting that larger instances are more computationally demanding, since both the generation and the solution phases are much time consuming. Tables 9 and 10 in Appendix B show the time reduction for the three sets of problems, for both probability levels. We note that the problems belonging to the first set are divided into two subsets, namely Sets 1A and 1B. Fig. 4 summarizes the reported results showing the average time reduction as function of the beam width for p = 0.90. We note that typically for a fixed beam width, higher savings can be obtained for larger instances. By analyzing Fig. 4, we observe, however, that for b = 20, time reduction for problems belonging to the first set is higher than that of the third set. This behaviour can be explained by observing that because of the specific problems structure, it may occur that the optimal solution is determined during the computation, thus, reducing the number of explored points. In any case, we expect that better performance can be achieved when larger test problems are considered. Resorting to heuristic approaches is in
43
general inevitable for solving practical-sized instances of most integer programming problems. This requirement is even stronger in our case, since the complexity related to the stochastic nature is added to that related to the solution of integer problems. The other issue to consider in evaluating the efficiency of a heuristic approach is the solutions quality. As one may expect, wider beams lead to better solutions, because the pool of candidate nodes is larger. The numerical results show that the Beam Search strategy exhibits very good performance at least on the test problems considered here. As matter of fact, we observe that even for a beam width of 5 the maximum error (recorded only for one instance) of the heuristic solution with the respect to optimal solution found by the Branch and Bound scheme is around 5%. Furthermore, in some cases the Beam Search algorithm finds the optimal solution. The two main issues discussed above determine the best beam width value. In setting this value, one has to consider the trade-off between the solutions quality and the savings in the solution time. The computational experiments show that the choice of good beam widths depends on the considered problems and requires intensive computational testing and fine tuning.
Fig. 4. Average time reduction versus beam width for the three sets of test problems with p = 0.90.
44
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
The procedure terminates when the current solution may not be improved. Several variants of this basic scheme could be considered (i.e. Tabu´ search, Simulated Annealing). Finally, we observe that in solving stochastic integer problems under probabilistic constraints two levels of approximation may be considered: the first one relies on partially exploring the space of efficient points and second one on heuristically solving the corresponding problems. The design of a two-level heuristic strategy is the subject of ongoing research.
5. Conclusion and future research directions In this paper we have proposed a heuristic strategy to solve stochastic integer problems under probabilistic constraints. The heuristic method implements a beam search strategy which is based on the classical Branch and Bound scheme. At any level of the search tree, only the promising nodes are kept for further exploration, whereas the remaining are pruned out permanently. The computational experiments show that the solution time of the Beam Search is far less than the time required by the Branch and Bound method. Furthermore, the Beam Search algorithm turns out to provide a high quality solution even for low beam width. Hence, the algorithm appears to be efficient and particularly suitable to solve the class of stochastic integer problems under probabilistic constraints. A future research direction is the design and implementation of other heuristic strategies. A first alternative is represented by the classical local search applied in the space of efficient points. In particular, starting from a p-efficient point we generate its neighbourhood, and we move to another point only if an evaluation function is improved.
6. Uncited references [4,5,8,11,12,14–16,18,20]
Appendix A Tables 5–8 report the values of the distribution parameters ai, ki of the circular and star dependency schemes for the additional 100 components of the test problem Tests 13–16. The first 200 components are reported in the tables of Appendix of the paper [2]. The tests have been defined in such a
Table 5 Distribution parameters of Test 13 Group
k1
k2
k3
k4
k5
k6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.790 0.810 0.200 0.810 0.200 0.010 0.809 0.810 0.018 0.601 0.018 0.011 0.011 0.011 0.005 0.090 0.710 0.020 0.011 0.005
0.005 0.010 0.090 0.010 0.890 0.810 0.980 0.810 0.011 0.011 0.811 0.011 0.911 0.011 0.005 0.020 0.140 0.020 0.011 0.005
0.005 0.100 0.001 0.100 0.010 0.900 0.710 0.900 0.011 0.041 0.011 0.089 0.011 0.089 0.805 0.910 0.150 0.810 0.089 0.805
0.990 0.051 0.700 0.951 0.700 0.810 0.020 0.810 0.011 0.611 0.011 0.770 0.811 0.700 0.090 0.050 0.150 0.910 0.700 0.090
0.007 0.010 0.300 0.010 0.300 0.850 0.061 0.005 0.911 0.971 0.011 0.009 0.811 0.009 0.007 0.010 0.010 0.010 0.009 0.007
0.907 0.810 0.300 0.810 0.010 0.101 0.061 0.001 0.615 0.911 0.615 0.011 0.189 0.011 0.900 0.710 0.710 0.010 0.011 0.210
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
45
Table 6 Distribution parameters of Test 14 Group
a1
a2
a3
a4
a5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.010 0.001 0.010 0.500 0.090 0.990 0.010 0.050 0.010 0.010 0.010 0.820 0.030 0.051 0.150 0.500 0.010 0.810 0.030 0.951
0.410 0.100 0.020 0.010 0.900 0.007 0.810 0.050 0.100 0.010 0.410 0.200 0.020 0.410 0.550 0.200 0.410 0.200 0.030 0.010
0.420 0.002 0.030 0.210 0.890 0.907 0.090 0.200 0.002 0.050 0.420 0.090 0.810 0.010 0.090 0.100 0.420 0.090 0.810 0.810
0.051 0.300 0.110 0.010 0.050 0.900 0.090 0.070 0.300 0.010 0.051 0.810 0.010 0.200 0.070 0.110 0.051 0.810 0.010 0.200
0.010 0.500 0.110 0.100 0.030 0.910 0.190 0.900 0.500 0.040 0.010 0.700 0.100 0.090 0.100 0.810 0.010 0.020 0.100 0.890
Table 7 Distribution parameters of Test 15 Group
k1
k2
k3
k4
k5
k6
k7
k8
k9
k10
k11
1 2 3 4 5 6 7 8 9 10
0.890 0.090 0.820 0.710 0.011 0.189 0.911 0.911 0.601 0.011
0.890 0.090 0.061 0.020 0.011 0.011 0.189 0.841 0.011 0.011
0.005 0.790 0.061 0.061 0.089 0.011 0.010 0.611 0.041 0.089
0.005 0.005 0.810 0.961 0.770 0.089 0.710 0.871 0.611 0.770
0.990 0.005 0.810 0.100 0.009 0.770 0.700 0.911 0.971 0.009
0.007 0.990 0.900 0.010 0.011 0.009 0.810 0.018 0.911 0.011
0.907 0.007 0.810 0.100 0.011 0.011 0.050 0.711 0.018 0.011
0.090 0.907 0.050 0.020 0.911 0.011 0.101 0.011 0.811 0.911
0.910 0.900 0.801 0.800 0.011 0.911 0.909 0.011 0.011 0.011
0.010 0.010 0.809 0.500 0.811 0.911 0.890 0.011 0.110 0.811
0.810 0.810 0.890 0.200 0.811 0.911 0.020 0.615 0.615 0.811
Table 8 Distribution parameters of Test 16 Group
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
1 2 3 4 5 6 7 8 9 10
0.100 0.500 0.090 0.910 0.810 0.900 0.810 0.809 0.011 0.099
0.110 0.200 0.890 0.010 0.910 0.300 0.700 0.990 0.911 0.011
0.810 0.100 0.890 0.810 0.100 0.810 0.300 0.710 0.011 0.011
0.810 0.110 0.050 0.900 0.051 0.810 0.910 0.820 0.811 0.911
0.100 0.810 0.005 0.090 0.010 0.100 0.910 0.061 0.811 0.911
0.010 0.810 0.005 0.790 0.810 0.951 0.8100 0.061 0.189 0.911
0.100 0.810 0.990 0.005 0.200 0.010 0.900 0.810 0.011 0.911
0.020 0.910 0.007 0.005 0.090 0.810 0.851 0.810 0.011 0.901
0.300 0.001 0.907 0.990 0.010 0.200 0.850 0.900 0.089 0.189
0.500 0.810 0.900 0.007 0.700 0.890 0.101 0.810 0.770 0.010
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
46
way to reflect the freedom in the users choice; thus, different groups of dependent random variables may have different values of the distribution parameters selected from the range. For each dependency scheme, the test problems are grouped in Tables according to the size r of each group. We note that in the case of star dependency (Tests 13 and 15), we have included an additional random
variable with parameter kr + 1, denoting the center of the star.
Appendix B The following Tables 9 and 10 report the time reduction (in percentage) for the three sets of prob-
Table 9 Time reduction (in percentage) for p = 0.95 Test problem
b=5
b = 10
b = 20
b = 30
Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test
32.98 52.25 39.80 50.75 43.23 12.18 83.21 48.85 57.01 45.76 42.94 36.59 23.94 56.98 65.61 57.40
27.69 52.25 35.79 50.52 43.23 12.18 82.01 42.14 49.50 25.10 40.31 30.73 22.10 56.40 61.02 39.93
21.52 52.25 35.79 97.67 43.23 12.18 79.98 39.17 32.19 22.99 39.55 26.83 22.28 57.22 70.13 19.55
21.52 52.25 29.77 42.06 43.23 12.18 70.02 18.27 22.75 21.49 30.48 24.88 18.04 50.04 56.91 24.62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Table 10 Time reduction (in percentage) for p = 0.90 Test problem
b=5
b = 10
b = 20
b = 30
Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test Test
63.49 80.94 43.49 87.26 62.56 80.66 53.79 94.77 63.65 77.63 91.06 75.31 82.83 75.77 80.00 71.56
54.93 80.68 40.06 87.26 62.56 79.48 46.52 94.69 51.15 71.66 85.41 63.94 79.56 75.55 79.12 69.71
36.05 78.59 36.82 86.87 92.16 77.62 29.77 94.51 48.25 55.81 80.20 78.44 75.18 73.86 77.71 54.63
24.88 68.67 31.74 81.62 62.56 54.31 5.70 93.77 31.43 67.54 66.24 75.57 65.22 51.95 74.98 43.18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
P. Beraldi, A. Ruszczyn´ski / European Journal of Operational Research 167 (2005) 35–47
lems. We observe that the test problems of the first set have been divided into two subsets: Sets 1A and 1B. References [1] J.E. Beasley, OR-library: Distributing test problems by electronic mail, Journal of the Operational Research Society 41 (1990) 1069–1072. [2] P. Beraldi, A. Ruszczyn´ski, The probabilistic set covering problem, Operations Research 50 (2002) 956–967. [3] P. Beraldi, A. Ruszczyn´ski, A Branch and Bound method for stochastic integer problems under probabilistic constraints, Optimization Methods and Software 17 (2002) 359–382. [4] A. Charnes, W.W. Cooper, Chance-constrained programming, Management Science 5 (1958) 73–79. [5] A. Charnes, W.W. Cooper, G.H. Symonds, Cost horizons and certainty equivalents: An approach to stochastic programming of heating oil, Management Science 4 (1958) 235–263. [6] CPLEX, ILOG CPLEX 6.5: Users Manual, CPLEX Optimization, Inc., Incline Village, NV, 1999. [7] D. Dentcheva, A. Pre´kopa, A. Ruszczyn´ski, Bounds for probabilistic integer programming problems, Discrete Applied Mathematics 124 (2002) 55–65. [8] D. Dentcheva, A. Pre´kopa, A. Ruszczyn´ski, Concavity and efficient points of discrete distributions in probabilistic programming, Mathematical Programming 89 (2001) 217–232. [9] M. Gendreau, G. Laporte, R. Se´guin, Stochastic vehicle routing, European Journal of Operational Research 88 (1996) 3–12.
47
[10] B.T. Lowerre, The HARPY speech recognition system, Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA. [11] L.B. Miller, H. Wagner, Change-constrained programming with joint constraints, Operations Research 13 (1965) 930– 945. [12] G.L. Nemhauser, L.A. Wolsey, Integer and Combinatorial Optimization, John Wiley & Sons, New York, 1988. [13] S.H. Owen, M.S. Daskin, Strategic facility location: A review, European Journal of Operational Research 111 (1998) 423–447. [14] A. Pre´kopa, On probabilistic constrained programming, in: Proceedings of the Princeton Symposium on Mathematical Programming, Princeton University Press, 1970, pp. 113– 138. [15] A. Pre´kopa, Contributions to the theory of stochastic programming, Mathematical Programming 4 (1973) 202– 221. [16] A. Pre´kopa, Dual method for solution of one-stage stochastic programming with random RHS obeying a discrete probability distribution, Zeitschrift fu¨r Operations Research 34 (1990) 441–461. [17] A. Pre´kopa, Stochastic Programming, Kluwer Academic Publishers, Boston, MA, 1995. [18] A. Pre´kopa, D. Vizva´ri, T. Badics, Programming under probabilistic constraint with discrete random variable, in: F. Giannesi (Ed.), New Trends in Mathematical Programming, Kluwer Academic Publishers, 1998, pp. 235–255. [19] I. Sabuncuoglu, M. Bayiz, Job shop scheduling with beam search, European Journal of Operational Research 118 (1999) 390–412. [20] S. Sen, Relaxations for the probabilistically constrained programs with discrete random variables, Operations Research Letters 11 (1992) 81–86.