- Email: [email protected]
The quadratic three-dimensional assignment problem: Exact and approximate solution methods
European Journal of Operational Research 184 (2008) 416–428 www.elsevier.com/locate/ejor
Discrete Optimization
The quadratic three-dimensional assignment problem: Exact and approximate solution methods Peter M. Hahn a,*, Bum-Jin Kim a, Thomas Stu¨tzle b, Sebastian Kanthak c, William L. Hightower d, Harvind Samra e, Zhi Ding e, Monique Guignard f a
f
Electrical and Systems Engineering, University of Pennsylvania, Philadelphia, PA 19104-6315, USA b CoDE, IRIDIA, CP 194/6, Universite´ Libre de Bruxelles, 1050 Brussels, Belgium c Computer Science, Darmstadt University of Technology, 64283 Darmstadt, Germany d Mathematics and Computer Science, High Point University, High Point, NC 27262, USA e Electrical and Computer Engineering, University of California, Davis, Davis, CA 95616, USA Operations and Information Management, The Wharton School, University of Pennsylvania, Philadelphia, PA 19104-6340, USA Received 8 November 2005; accepted 14 November 2006 Available online 11 January 2007
Abstract This paper reports on algorithm development for solving the quadratic three-dimensional assignment problem (Q3AP). The Q3AP arises, for example, in the implementation of a hybrid ARQ (automatic repeat request) scheme for enriching diversity among multiple packet re-transmissions, by optimizing the mapping of data bits to modulation symbols. Typical practical problem sizes would be 8, 16, 32 and 64. We present an exact solution method based upon a reformulation linearization technique that is one of the best available for solving the quadratic assignment problem (QAP). Our current exact algorithm is useful for Q3AP instances of size 13 or smaller. We also investigate four stochastic local search algorithms that provide optimum or near optimum solutions for large and difficult QAP instances and adapt them for solving the Q3AP. The results of our experiments make it possible to get good solutions to signal mapping problems of size 8 and 16. 2006 Elsevier B.V. All rights reserved. Keywords: Quadratic assignment; Branch-and-bound; Heuristics; Three-index assignment; OR in telecommunications
1. Introduction Pierskalla (1967) introduced the quadratic three-dimensional assignment problem (Q3AP) in a technical memorandum. The work was never published in the open literature. Since then, nothing on the subject has * Corresponding author. Present address: 2127 Tryon Street, Philadelphia, PA 19146-1228, USA. Tel.: +1 215 546 3413; fax: +1 215 546 4043. E-mail addresses: [email protected] (P.M. Hahn), [email protected] (B.-J. Kim), [email protected] (T. Stu¨tzle), [email protected] (S. Kanthak), [email protected] (W.L. Hightower), [email protected] (H. Samra), [email protected] (Z. Ding), [email protected] (M. Guignard).
0377-2217/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2006.11.014
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
417
appeared in the literature. The authors re-discovered the Q3AP while working together on a problem arising in data transmission system design. We will explain this problem, after giving the formulation of the quadratic assignment problem (QAP) and the Q3AP. The Lawler (1963) formulation of the QAP is given by ( ) N X N X N X N X C ijkn xij xkn : x 2 X; x binary ; ð1Þ min i¼1
j¼1
k¼1 n¼1
PN PN where x 2 X fx P 0 : i¼1 xij ¼ 1 for j = 1, . . . , N; j¼1 xij ¼ 1 for i = 1, . . . , N}. The Koopmans and Beckmann (1957) formulation of the QAP assumes that the four-dimensional cost matrix C is generated from a pair of N · N matrices, F and D, with Cijkn = fik Æ djn (i, j, k, n = 1, . . . , N). The Q3AP is an extension of the QAP that can be represented as a three-dimensional assignment problem and it is given by the following formulation: 9 8 > > > > = > > > i¼1 j¼1 p¼1 k¼1 n¼1 q¼1 ; : i¼1 j¼1 p¼1 k6¼i
n6¼j
q6¼p
where I, J and P are disjoint sets of constraints: ( ) N X N X I¼ xP0: xijp ¼ 1 for i ¼ 1; . . . ; N ; j¼1
( J¼
xP0:
P¼
N X N X i¼1
( xP0:
p¼1
p¼1
N X N X i¼1
ð3Þ
) xijp ¼ 1 for j ¼ 1; . . . ; N ;
ð4Þ
) xijp ¼ 1 for p ¼ 1; . . . ; N :
ð5Þ
j¼1
The problem is so named because the object is to optimize a quadratic function over the three-dimensional assignment polytope x 2 I \ J \ P. The quadratic expression in the objective function need not include any xijpxknq terms for which i = k or j = n or p = q unless all three equalities hold. If i = k and j = n and p = q, then xijpxknq = xijp; otherwise xijpxknq = 0. Being an extension of the NP-hard problems QAP and the threedimensional assignment problem (3AP) (Garey and Johnson, 1979), it is easy to see that the Q3AP is also NP-hard. The Q3AP was motivated, in part, by a problem arising in wireless data transmission systems, where information data bits for transmission are organized into fixed length packets. Just before transmission, several contiguous bits are mapped into distinct discrete symbols during modulation. At the receiver, these symbols are demodulated and mapped back into the original bit strings, which are then put back into data packets (for delivery by service providers). Due to interference, receiver noise, channel fading, and channel distortions, erroneous data packets may be received in these systems, which can be detected by the receiver through a cyclical redundancy checking (CRC) code inserted at the transmitter. To ensure data service quality, re-transmission of the erroneous packet may be requested upon a CRC failure by sending, over a feedback channel, an automatic repeat request (ARQ). Recent research works (Samra and Ding, 2003a,b,c) have shown that the mapping of binary data in the retransmitted packet into symbols should be changed for each re-transmission in order to significantly improve the joint detection of the original erroneous bit strings from multiple observations. Loosely speaking, for this ARQ mechanism to be as effective as possible, the symbol mappings for the repeated transmission should provide as much diversified representation (or statistical independence) as possible among the resulting received transmission symbols. In the parlance of communication systems, providing statistical independence among signals representing the same data is called transmission diversity. Providing such diversity between re-transmissions of the same intended data can be done on a pair-wise basis and would constitute solving a QAP. In this case, the cost coefficients {Cijkn} indicate the cost of simultaneously assigning bit strings i and k to the transmitted symbols j and n. However, better performance is achieved by making the
418
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
original transmission and the two following repeats as statistically independent as possible. This latter optimization requires the Q3AP.1 In that case, the cost coefficients {Cijpknq} refer to the simultaneous assignment of the symbols j and n to the bit strings i and k in the first re-transmission and the assignment of symbols p and q in the second re-transmission. In both cases, the cost coefficients for the QAP and the Q3AP can be computed, based on some usual assumptions on the error model of the transmission errors (corruption by additive white Gaussian noise with fixed variance). For details of these computations, we refer to (Samra et al., 2005). Note, that in the QAP case, there are N! possible mappings, while in the Q3AP case, we have N! · N! possible mappings. Fortunately, the optimization of the diversity method needs to be performed only once in the design phase of any actual implementation. There is no need for real time optimization. Once the optimal mapping is acquired for a given type of wireless fading channel, the transmitter/receiver can simply use a lookup table to determine and agree with each other on the next mapping. Results in Samra et al. (2003, 2005) demonstrate tremendous gains, even with suboptimal QAP mapping. For four common wireless fading scenarios from flat fading to frequency selective channels, their results show more than 3–4 dB of signal-to-noise ratio (SNR) gain, which implies more than 50% power savings or 100% increase in signal coverage area. Moreover, it is shown in Samra et al. (2003, 2005) that mapping diversity with only one re-transmission delivers better performance than simple ARQ with three re-transmissions. Again, there is a 50% reduction in bandwidth usage or 100% improvement in capacity. In this paper, we study algorithms for tackling the Q3AP, including exact algorithms and heuristics. Our prior computational experience in solving the QAP forms the basis for this paper. We enjoyed marked success by solving difficult QAPLIB (Burkard et al., 1991) instances that previously were unsolved, including the Krarup30a (Hahn and Krarup, 2001) and the Nugent et al. (1968) instance of size N = 25. To solve both, we implemented a bound computation (Hahn et al., 1998) based on a dual ascent procedure that efficiently solves the level-1 RLT formulation (Adams and Sherali, 1986, 1990). In Guignard et al. (2005) and Kim (2006), we extended these bounds to the Q3AP and we conducted experimental tests to determine the utility of level-1 RLT bounds for solving the Q3AP. We designed and implemented a basic branch-and-bound algorithm, and applied it to a set of benchmark instances created from difficult benchmark QAP instances found on QAPLIB. These chosen academic QAP problem instances are noted for their difficulty. Thus, the Q3AP instances derived from them are also difficult. We also applied the branch-and-bound algorithm to a size 8 Q3AP corresponding to the dual mappings required for two 8-PSK re-transmissions. This latter instance is hereafter referred to as REAL-8. Another realistic problem of size 16 was generated and attempts were made to solve it exactly. This problem (denoted REAL-16) represented the dual mappings required for two 16-QAM (quadrature amplitude modulation) re-transmissions. So far, we have experimental exact solution results to the four Q3AP instances derived from QAP, of sizes N = 8, 9, 12 and 13 and to the REAL-8 instance. Due to the eight symmetries, the REAL-8 instance is solved in less time than required to solve the Nug8-derived instance, which has only four symmetries. However, the computation times for all these instances are substantial and therefore we also studied the application of stochastic local search (SLS) methods (Hoos and Stu¨tzle, 2004). In fact, computational experience with SLS methods on the QAP has shown that they typically reach high quality or optimal solutions much quicker than exact algorithms. Here, we adapted previously proposed SLS algorithms for the QAP to solve the Q3AP, including a simulated annealing (SA) algorithm of Connolly (1990), the Robust Tabu Search algorithm of Taillard (1991), the FANT algorithm by Taillard (1998), and an iterated local search algorithm (ILS) by Stu¨tzle (2006). As shown in an experimental study of these algorithms, they can reach high quality solutions in reasonable computation times. In fact, the best performing of these methods improves the best-known objective function value for the REAL-16 instance by about 19%, when compared with the best results of other heuristic methods. In Section 2, we present a level-1 reformulation linearization technique (RLT) relaxation on the Q3AP and in Section 3 we demonstrate a lower bounding method based on level-1 RLT. RLT bounding methods have proved successful in branch-and-bound algorithms for the QAP and appear to work well on the 1
In fact, performance is improved by making the following 3, 4, . . . , X re-transmissions as statistically independent as possible, which would require Q4AP, Q5AP, . . . , QXAP optimization. However, doing this would lead to diminishing returns, at the cost of increased system complexity. Our results show that there is still a significant benefit when moving from the optimization of one re-transmission to the optimization of the two re-transmissions, that is, when moving from the QAP to the Q3AP.
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
419
Q3AP. Experimental results of our exact solution method are presented in Section 4. In Section 5, we introduce four SLS algorithms for this very difficult problem. Results of our experiments with SLS algorithms are presented in Section 6. Section 7 presents our conclusions and our plans for future work. 2. Reformulation linearization technique (RLT) on the Q3AP The RLT was first introduced in Adams and Sherali (1986) and Adams and Sherali (1990) and best described in Sherali and Adams (1999). It is a technique for generating a hierarchy of polyhedral representations for linear and polynomial 0–1 programming problems. These new representations are extended linear programs that permit the construction of Lagrangean duals that step-by-step come closer to solving exactly the original problem. This is demonstrated in Adams et al. (2007) for RLT level-2. In this paper we take advantage of only the first step (level-1) of this hierarchy. As suggested by its name, the RLT methodology consists of two steps: reformulation and linearization. We describe in this section how this methodology applies to the Q3AP, and refer to Sherali and Adams (1990, 1994, 1999) for more general and detailed explanations. Consider the Q3AP as presented in Eqs. (2)–(5). The level-1 RLT formulation is generated via the following two steps: Reformulation Multiply each of the 3N equations and each of the N3 non-negative restrictions x P 0 defining I, J, P in (3)–(5) by each binary variable xknq, and append these new restrictions to (3)–(5), respectively. Substitute xknq ¼ x2knq throughout the objective function and constraints. When a variable xijp in a given constraint is multiplied by xknq, express the resulting product as xijpxknq in that order. (Observe that for binary variables, such multiplications, together with the assignment constraints, imply that xijpxknq = 0 if i = k or j = n or p = q unless all three equalities hold.) Linearization Linearize the resulting problem by substituting, for every occurrence of each product xijpxknq with i 5 k, j 5 n, and p 5 q, the continuous variable vijpknq. Enforce the trivial but extremely important restrictions that vijpknq = vknqijp for all (i, j, p, k, n, q) with i < k, j 5 n and p 5 q. These restrictions define what we call complementary pairs. When applied to the bound computation, these have a dramatic effect on increasing the bound. The level-1 RLT formulation for the Q3AP is generated as instructed, which results in: RLT: XXX XXXXXX Minimize : bijp xijp þ C ijpknq vijpknq i
Subject to :
j
XX k k6¼i
n n6¼j
k k6¼i
q q6¼p
n n6¼j
q q6¼p
XX XX
p
i
j
p
k k6¼i
n n6¼j
ð6Þ
q q6¼p
vijpknq ¼ xijp
8ði; j; p; qÞ; q 6¼ p;
ð7Þ
vijpknq ¼ xijp
8ði; j; p; nÞ; n 6¼ j;
ð8Þ
vijpknq ¼ xijp
8ði; j; p; kÞ; k 6¼ i;
ð9Þ
vijpknq ¼ vknqijp 8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q; vijpknq P 0 8ði; j; p; k; n; qÞ; i 6¼ k; j 6¼ n; p 6¼ q; x 2 I \ J \ P; x binary:
ð10Þ ð11Þ ð12Þ
Problem Q3AP and RLT are equivalent in the following sense. Given any feasible solution x to the Q3AP, there exists a v such that (x, v) is feasible to the RLT with the same objective value. Conversely, for any feasible solution (x, v) to the RLT, the corresponding x is feasible to the Q3AP with the same objective value. The proof is straightforward and follows the standard RLT proof arguments.
420
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
One can show that the problem RLT possesses a block-diagonal structure that can be exploited. In particular, when constraints (10) are not present, RLT reduces to N3 + 1 three-dimensional assignment problems, or 3AP’s, one over each (i, j, p) triplet in (7)–(9), and one over x 2 I \ J \ P. This suggests a generalization of the Gilmore–Lawler bound (GLB) for the QAP (Gilmore, 1962). However, the dual procedure bound proposed in Hahn and Grant (1998), of which the GLB is a special case, performs significantly better. 3. A level-1 RLT lower bound for the Q3AP In order to develop a Lagrangean dual procedure to solve the RLT, we first obtain a smaller reformulation via the substitution of variables suggested by constraint (10), thereby reducing the number of variables and constraints, in addition to halving the number of non-negativity restrictions in (11). Now, the formulation becomes: RLT1: XXX XXXXXX e ijpknq vijpknq Minimize : bijp xijp þ ð13Þ C i
Subject to :
j
p
i
j
p
k k>i
e ijpknq ¼ C ijpknq þ C knqijp ; where C XX XX vijpknq þ vknqijp ¼ xijp k k>i
n n6¼j
XX k k>i
q q6¼p
n n6¼j
q q6¼p
n n6¼j
q q6¼p
XX XX
q q6¼p
8ði; j; p; qÞ; q 6¼ p;
ð14Þ
8ði; j; p; nÞ; n 6¼ j;
ð15Þ
n n6¼j
k k
vijpknq þ
n n6¼j
XX k k
vknqijp ¼ xijp
q q6¼p
vijpknq ¼ xijp
8ði; j; p; kÞ; k > i;
ð16Þ
vknqijp ¼ xijp
8ði; j; p; kÞ; k < i;
ð16aÞ
8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q;
ð17Þ
vijpknq P 0
x 2 I \ J \ P; x binary:
ð18Þ
Now we consider the following Lagrangean dual problem LD1, whereby constraints (14), (15), (16), (16a) are placed into the objective function using multipliers q, k, g and u, respectively. LD1: maximize hðq; k; g; uÞ where 9 8 > > = > n q ; : i j p i j p k k>i
n6¼j
q6¼p
where ~ bijp ¼ bijp þ
X q q6¼p
qijpq þ
X n n6¼j
kijpn þ
X k k>i
gijpk þ
X
uijpk
8ði; j; pÞ;
ð20Þ
k k
e ijpknq qijpq qknqp kijpn kknqj gijpk uknqi C ijpknq ¼ C
8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q;
ð21Þ
subject to constraints (17) and (18). Next, we consider a further Lagrangean dual problem LD2, whereby equality constraints in (18) on x variables are also placed into the objective function using multipliers a, b and c, respectively.
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
421
LD2: maximize hða; b; c; q; k; g; uÞ where 9 8 > > = > n q ; : i j p i j p i j p k k>i n6¼j q6¼p
ð22Þ
where bijp ¼ ~ bijp ai bj cp
8ði; j; pÞ;
ð23Þ
subject to: vijpknq P 0 0 6 xijp 6 1
8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q;
ð24Þ
8ði; j; pÞ:
ð25Þ
Notice that we now relax the equivalent formulation by allowing xijp to be fractional between 0 and 1. 4. Level-1 RLT branch-and-bound results To demonstrate the potential of the approach, we implemented a branch-and-bound algorithm that uses the level-1 RLT bounds. Our branch-and-bound algorithm follows the same strategies as those utilized for the level-1 RLT branch-and-bound algorithm for solving large QAP instances (Hahn et al., 2001). Table 1 summarizes our experimental results on five Q3AP instances of size 8, 9, 12, and 13. One of these instances (designated REAL-8) is a size 8 symbol-mapping problem (or the 8 phase shift keying modulation known as 8-PSK), while the other four instances were derived from QAP instances. Three of these instances were created from the now famous Nugent problem instances (Nugent et al., 1968) found in QAPLIB (Bur´ ric Taillard. This latter kard et al., 1991). One instance, ‘‘Tai9’’ was derived from an instance posed by E instance involves the placement of 10 incrementally heavier blades on the periphery of a turbine shaft. Since the placement of the first blade is arbitrary, this instance reduces to a QAP of size 9. All these QAP instances are Koopmans–Beckmann. (For a definition of Koopmans–Beckmann, see the first paragraph of Section 1.) To generate a six-dimensional cost matrix for one of our Q3AP instances, we used the following relationship Cijpknq = fik Æ djn Æ fik Æ dpq; (i, j, p, k, n, q = 1, . . . , N). One might wonder about the validity of testing a general Q3AP algorithm on problem instances that are Koopmans–Beckmann. In fact, our experimental results show that the Koopmans–Beckmann generated problems are harder to solve than non-Koopmans–Beckmann problems of the same size arising from wireless communication examples. The size 8, 9 and 12 problem instances were run on a Sun Ultra 10 workstation with a single 360 MHz CPU. The size 13 problem instance, derived from the Nug15 instance, was run on a single 733 MHz CPU of a Dell 7150 server. The runtimes in Table 1 are normalized for the Dell server. The Nugent derived Q3AP test problems required significantly more time to solve than the QAP instances from which they were created, with the size 8 instance taking longer than it would to exhaustively enumerate Table 1 Experimental results on Level-1 RLT branch-and-bound algorithm
Root lower bound Optimum value No. of optima Starting upper bound Dell 7150 runtime Number of nodes a
8-PSK
Nug8 derived
Tai9 derived
Nug12 derived
Nug13a derived
5,909,755,156 15,328,457,328 6144 9E + 18 155.8 seconds 267,264
91.5 134.0 2 9E + 18 48.5 seconds 621
55,250,527,200,558 74,755,545,715,200 32 9E + 18 10.3 minutes 12,666
229.8 580.0 3 9E + 18 14.8 hours 65,597
899.7 1912.0 2 9E + 18 23.1 days 1,790,840
Cutting down the flow and distance matrices of the Nug15 problem created the Nug13 instance.
422
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
the 8! · 8! feasible solutions (at 10 nanoseconds per solution). However, the 14.8 hours it takes to solve the size 12 problem instance is more reasonable, considering that exhaustive enumeration would take over 70 years. 5. Approximate solution methods Stochastic local search methods belong to the most successful techniques for the approximate solution of hard optimization problems (Aarts and Lenstra, 1997; Hoos and Stu¨tzle, 2004). These methods are also the best way to provide feasible solutions in large industrial problems with tight constraints. In this work, we adapted four successful SLS algorithms for the QAP to the Q3AP: the simulated annealing algorithm by Connolly (1990), the robust tabu search (RoTS) algorithm and fast ant system (FANT) by Taillard (1991, 1998), and an iterated local search (ILS) algorithm by Stu¨tzle (2006). Computational experience has shown that these algorithms can find optimal solutions to a wide variety of QAP instances much faster than the best performing exact algorithms. For example, the ILS algorithm finds the best-known solution to the above mentioned Nugent 30 instance in about 1.3 CPU seconds, on average (measured on a 1.2 GHz MP Athlon CPU). These results on the QAP suggest that if one wants to solve the Q3AP with sizes as large as N = 64, these methods are probably the only means to get good quality solutions. Although, these methods cannot prove optimality of the solutions, they can be useful in providing good initial upper bounds for exact methods. 5.1. Iterative improvement for the Q3AP All four SLS algorithms we implemented rely on an underlying iterative improvement algorithm. Local search algorithms for the QAP usually exploit a solution representation based on permutations. In fact, an alternative representation of the general QAP formulation in Eq. (1) is given by min
N X N X i¼1
ð26Þ
C i/ðiÞk/ðkÞ ;
k¼1
where / is a permutation of the set of integers {1, . . . , N}. ForP QAP instances having the Koopmans–Beckman N PN property, this would result in the minimization of the term i¼1 k¼1 F ik D/ðiÞ/ðkÞ , where Cijkl = FikDjl. Similarly, a solution to the Q3AP, using the formulation given by Eqs. (2)–(5), can be represented using two permutations / and u of the set of integers {1, . . . , N} as ( ) N N X N X X min f ð/; uÞ ¼ bi/ðiÞuðiÞ þ C i/ðiÞ/ðjÞjuðiÞuðjÞ : /; u are permutations over the set f1; . . . ; N g : i¼1
i¼1
j¼1
ð27Þ Given this formulation, it is also easy to see that the search space size of the Q3AP is N! · N! Iterative improvement algorithms for the QAP typically make use of the 2-exchange neighborhood relation, in which two solutions are neighbored if their permutations differ in exactly two positions. That is, the neighborhood N(/) of permutation / is given by 0
0
0
N ð/Þ ¼ f/0 j/ðrÞ ¼ /ðsÞ; /ðsÞ ¼ /ðrÞ; r 6¼ s and /ðiÞ ¼ /ðiÞ 8i 6¼ r; sg:
ð28Þ
For the Q3AP, one can apply this 2-exchange neighborhood to each of / and u. Hence, two solutions for the Q3AP are neighbored, if either / or u differ in exactly two positions and we have that N ð/; uÞ ¼ fð/0 ; u0 Þj/ðrÞ0 ¼ /ðsÞ; /ðsÞ0 ¼ /ðrÞ; r 6¼ s and /ðiÞ0 ¼ /ðiÞ8i 6¼ r; s and uðiÞ0 ¼ uðiÞ8i or uðrÞ0 ¼ uðsÞ; uðsÞ0 ¼ uðrÞ; r 6¼ s and uðiÞ0 ¼ uðiÞ8i 6¼ r; s and /ðiÞ0 ¼ /ðiÞ8ig:
ð29Þ
Essential to the performance of iterative improvement algorithms can also be the strategy that is used to explore the neighborhood. The most important of these are the first-improvement strategy that accepts an improving neighboring solution as soon as it is found, and the best-improvement strategy that first examines
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
423
the whole neighborhood and then chooses a neighboring solution that gives the largest improvement in the objective function value. For the Q3AP, one may distinguish additional strategies for exploring the neighborhood, depending on whether / and u are explored in sequence or in parallel when searching for an improved neighboring solution. Finally, it has to be mentioned that the evaluation of neighboring solutions also needs to be adapted from the QAP case, which can be done in a rather straightforward way. The necessary adaptations of the four SLS algorithms we employ here concern mainly the details of the local search, as outlined above. 5.2. Simulated annealing (SA) for the Q3AP In SA, worsening moves are accepted probabilistically by a criterion that takes into account the deterioration incurred in the objective function value and a time varying parameter T called temperature. We adapted the simulated annealing of Connolly (1990), which is known to be one of the best performing SA algorithms for the QAP, to the Q3AP. In this algorithm, the neighborhood is scanned in a fixed order. A neighboring solution s 0 is accepted deterministically, if it is better than or equal to the current one, s, or, if s 0 is worse, with a probability given by exp[(f(s 0 ) f(s))/T], where f is the objective function (Metropolis et al., 1953). In our SA algorithm (SA-Q3AP), we use the same parameter settings as that of the best performing SA variant identified by Connolly for the QAP. At each iteration, SA-Q3AP generates two neighboring solutions: One, which differs in the position of two elements in permutation / and one, which differs in the position of two elements in permutation u. Then, the better of the two neighbors is used for the standard acceptance test of the SA algorithm. Preliminary experiments showed that this strategy for examining the neighborhood gives best performance compared to generating only one neighbor at each iteration or when evaluating all neighbors and then using the best for the acceptance test. In preliminary experiments, SA-Q3AP was shown to outperform an adaptation of the Simulated Annealing algorithm of Burkard and Rendl (1984) to the Q3AP. 5.3. Tabu search for the Q3AP In tabu search, escapes from local optima are possible due to the prohibition of moves to recently visited solutions. This is done by keeping track of parts of the search history, mainly in the form of attributes of recently performed moves between neighboring solutions. The tabu search algorithm we use (RoTS-Q3AP) is a straightforward extension of Taillard’s robust tabu search algorithm (Taillard, 1991). RoTS-Q3AP first evaluates all possible neighboring solutions and then moves to the best ‘‘non-tabu’’ neighbor. A neighboring solution is considered ‘‘tabu’’, if during the last tt iterations we had that p(k) = i and p(l) = j and that in the current solution p(k) 5 i and p(l) 5 j, where p = /, u and tt is a parameter called tabu tenure. In other words, it is forbidden to assign two positions in either of the two permutations the same value they had in the previous tt iterations. Additionally, if for r iterations we have p(k) 5 i, i.e., the value i was not assigned to position k in the previous r iterations, all solutions that do not assign value i to position k are forbidden; this actually forces an assignment of the value i to p(k) (forced moves). The parameter tt was set as originally proposed in the robust tabu search algorithm, while, based on some preliminary experiments, r was set to 10n3. Additionally, an aspiration criterion is used that overrides the tabu status of a solution, if it is better than the previously best-found solution. 5.4. Ant colony optimization for the Q3AP Ant algorithms and, in particular, Ant Colony Optimization (ACO) (Dorigo and Stu¨tzle, 2004) are algorithmic techniques that take inspiration from real ants behavior. In ACO, artificial ants stochastically construct solutions biased by (artificial) pheromone trails, which serve as a numerical information that is modified at the algorithm’s execution time to reflect the ants’ search experience. We applied fast ant system (FANT), which was previously studied for the QAP (Taillard, 1998). For an overview of ant algorithm applications to the QAP we refer to Stu¨tzle and Dorigo (1999). As in all other applications of ant algorithms to the QAP, the pheromone value sij refers to setting p(i) = j. In our adaptation of FANT, FANT-Q3AP, we use two pheromone matrices, one for generating each of the two permutations / and u. For the solution construction, the same parameters as proposed in Taillard (1998) are used. The only new feature we use is analogous to the
424
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
‘‘forced moves’’ of robust tabu search: If for r iterations p(k) 5 i, this assignment is deterministically adopted in the solution construction. As in RoTS-Q3AP, a good parameter setting was found to be r = 10n3. Once a solution is constructed, it is improved using a straightforward adaptation to the Q3AP of the local search procedure originally proposed for FANT. 5.5. Iterated local search for the Q3AP Iterated local search (ILS) is a conceptually simple SLS method that builds a chain of local optima by cycling through phases of solution perturbations, which are applied to local optima, local search, and acceptance test for defining the incumbent solution (Lourenco et al., 2002). For the Q3AP, we implemented a variant of the ILS algorithms studied by Stu¨tzle (2006). The resulting ILS-Q3AP algorithm works as follows. It first generates two random permutations for / and u and applies a first-improvement local search that systematically examines the 2-exchanges in / and u. The main loop of the ILS algorithm then first applies a perturbation, consisting of a random exchange of a number of positions in each of the two permutations, to a local optimum (/*, u*). This results in a perturbed solution (/ 0 , u 0 ) that is the starting solution for the next local search. Once a new locally optimal solution (/ 0 *, u 0 *) is reached, an acceptance criterion decides whether the next perturbation is applied to (/ 0 *, u 0 *) or to the previous local optimum (/*, u*). In ILS-Q3AP, worse solutions are accepted with a probability given by exp{(f((/ 0 *, u 0 *)) f((/*, u*)))/T}, where T is a parameter that is set such that a decrease of the objective function of the current solution (/*, u*) by 1% is accepted with a probability of 0.5. If for a large number of iterations (typically a few thousand) of the main ILS loop, no improved solution was found; then the search is restarted from a new, randomly generated initial solution. 6. Approximate solution method results We tested the four SLS methods on all Q3AP instances of Section 4 plus an additional one of size 15, which was derived from the Nug15 instance from QAPLIB. For each of the four algorithms, 8 runs were performed on a single 440 MHz CPU of an HPJ5000 machine, which has almost exactly the same speed for our codes as the Dell 7150 machine on which the branch-and-bound algorithm was run. On the QAP derived instance of size 8 and 9, all the algorithms found the optimal solution within a few seconds of computation time, and the results are not shown explicitly here. The experimental results on the instances of size 12–15 are summarized in Tables 2–4. In each of these tables is given the average (over 8 runs) objective function value found by the four SLS algorithms after 1 minute, 10 minutes, 1 hour and 10 hours of computation time. The best average objective value for each time limit is in bold face. For the derived Nug12 and Nug13 instances, several of the SLS algorithms can reach the optimal solution. Overall, ILS-Q3AP appears to be the best-performing algorithm. It is the only heuristic to find, in every single run, the optimum solution for the size 8, 9, 12 and 13 instances and the best-known solution for the size 15 instance. It is also interesting to compare the performance of the four SLS algorithms to the branch-and-bound algorithm. Although SLS algorithms cannot prove optimality of the solutions they find, one can compare the computation times the different algorithms need to reach the first optimal solution. To do so, we ran ILS-Q3AP and RoTS-Q3AP 25 times and measured the average time, as well as the standard deviation (SD), each
Table 2 Average objective function values for SA-Q3AP, FANT-Q3AP, TS-Q3AP, and ILS-Q3AP for eight runs on Nug12 using different time limits from 1 minute to 10 hours
FANT-Q3AP RoTS-Q3AP SA-Q3AP ILS-Q3AP a
1 minute run
10 minute run
1 hour run
10 hour run
878 779 771 796
857 748 705 634
808 607 626 589
750 580a 583 580a
Optimum objective function value.
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
425
Table 3 Average objective function values for SA-Q3AP, FANT-Q3AP, TS-Q3AP, and ILS-Q3AP for eight runs on Nug13 using different time limits from 1 minute to 10 hour
FANT-Q3AP RoTS-Q3AP SA-Q3AP ILS-Q3AP a
1 minute run
10 minute run
1 hour run
10 hour run
2322 2118 2102 2120
2130 1995 2024 2014
2107 1918 1967 1962
2054 1912a 1912a 1912a
Optimum objective function value.
Table 4 Average objective function values for SA-Q3AP, FANT-Q3AP, TS-Q3AP, and ILS-Q3AP for 8 runs on Nug15 using different time limits from 1 minute to 10 hours
FANT-Q3AP RoTS-Q3AP SA-Q3AP ILS-Q3AP
1 minute run
10 minute run
1 hour run
10 hour run
3172 2843 2752 2751
3019 2642 2631 2568
2943 2547 2504 2399
2703 2324 2400 2230
algorithm took to reach the optimal solution. We ran one run of the branch-and-bound algorithm, with the upper bound set to an extremely large value and reported the time it took to reach the first optimum. The results are presented in Table 5. From the computational results in Table 5, one can observe that the SLS algorithms reach optimal solutions about one to two orders of magnitude faster than the branch-and-bound algorithm for the instance sizes tested. We conjecture that a simple combination of the two, where an SLS algorithm is used to determine good upper bounds for the branch-and-bound, is highly desirable. In this context, it is interesting to note that the time-to-completion of the branch-and-bound algorithm was between a factor of 2 (Nug12 derived instance), to a factor of 12 (Tail9 derived instance) larger than the time to find the first occurrence of an optimal solution. Hence, if this factor is small, then such a combination can also speed up significantly the time to prove optimality. The excellent performance of SLS algorithms is confirmed on a Q3AP symbol-mapping instance (REAL16) of size 16. For REAL16, the matrix entries were originally given as double precision floating-point numbers but then scaled to truncated integers for running the SLS algorithms. The results for REAL16 are given in Table 6, again measured on a single 440 MHz CPU of an HPJ5000 machine. Best solutions are shown in bold face entries. Second best solutions are shown in italics. As can be seen, the ILS algorithm gives the best performance, if 10 minutes or more are allowed. In fact, for the largest allowed computation time of 10 hours, the worst solution identified by ILS-Q3AP across eight runs was even better than the best solution found by either Table 5 Runtime to reach optimum for ILS-Q3AP, RoTS-Q3AP and b-and-b algorithms for four Q3AP instances (run-times of ILS-Q3AP and RoTS-Q3AP are across 25 independent runs) Real size 8
Tai9 derived
NUG12 derived
NUG13 derived
ILS-Q3AP
Average time SD time Best time Worst time
0.043 seconds 0.023 seconds 0.01 seconds 0.09 seconds
1.65 seconds 1.31 seconds 0.2 seconds 5.13 seconds
2162.2 seconds (36.03 minutes) 2138.07 seconds 53 seconds 8404.57 seconds
7519.2 seconds (2.09 hours) 6801.27 seconds 69.44 seconds 27194.4 seconds
RoTS-Q3AP
Average time SD time Best time Worst time
0.44 seconds 0.58 seconds 0.01 seconds 2.57 seconds
4.06 seconds 2.86 seconds 0.04 seconds 14.64 seconds
9104.3 seconds (2.53 hours) 8107.64 seconds 169.58 seconds 35733.93 seconds
3002.0 seconds (50.03 minutes) 2752.14 seconds 128.24 seconds 9719.47 seconds
B&B
Time to 1st Opt
19.2 seconds
50.6 seconds
10.6 hours
165.35 hours
426
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
Table 6 Average objective function values and standard deviation for SA-Q3AP, TS-Q3AP and ILS-Q3AP for eight runs on REAL16 for different time limits from 1 minute to 10 hours
ILS-Q3AP RoTS-Q3AP SA-Q3AP
AVE SD AVE SD AVE SD
1 minute run
10 minute run
30 minute run
1 hour run
10 hours run
4.13730E+15 4.65885E + 14 3.61907E + 15 5.13961E + 14 3.40559E + 15 1.08742E + 14
2.88551E + 15 6.07584E + 14 3.17127E + 15 2.84639E + 14 3.04805E + 15 1.44465E + 14
2.66533E + 15 4.16076E + 14 2.92033E + 15 6.75849E + 13 2.82990E + 15 3.49652E + 14
2.46528E + 15 4.08002E + 14 2.86702E + 15 6.01477E + 13 2.73253E + 15 3.04737E + 14
2.20841E + 15 6.81110E + 13 2.73298E + 15 8.41498E + 13 2.58397E + 15 1.47923E + 14
The Quadratic Three-dimensional Assignment Problem: 4.5E+15 4E+15 3.5E+15 3E+15 2.5E+15 2E+15 1.5E+15 1E+15 5E+14 0
ILS I_SA TABU
1 min
10 min 30 min
1 hour 10 hour
Time
Fig. 1. Development of the average objective function value for ILS-Q3AP, SA-Q3AP, and RoTS-Q3AP applied to the hard, symbolmapping Q3AP instance of size N = 16.
SA-Q3AP or RoTS-Q3AP. The computational results are further illustrated in Fig. 1. It is noteworthy to mention that, after only 1 hour, the best result found by ILS-Q3AP is about 19% better than the previously bestknown solution. 7. Conclusions and future directions We present a new development of exact and local search algorithms for solving the Q3AP. These algorithms were adapted from their QAP counterparts. While for the local search algorithms these adaptations could be done in a reasonably straightforward way, for the exact algorithms this required the development and implementation of effective lower bounds for the 3AP. Although our computational results are encouraging, they also illustrate the level of difficulty associated with the Q3AP. Presently, our exact solution method, based upon the best techniques available for solving the QAP, is computationally feasible only for Q3AP instances of size 13 or smaller. Stochastic local search techniques, which provide optimum or near optimum solutions for large and difficult QAP instances, are essential for reaching high quality solutions to reasonable size Q3AP instances. Among the various SLS algorithms tested, the overall best performance was given by an iterated local search algorithm. However, thus far, even the best performing SLS algorithm incurs for small Q3AP instances computation times for reaching optimal or best-known solutions that are orders of magnitude larger than for similar size QAP instances. Clearly, much more work is needed on this challenging and yet important new combinatorial optimization problem. The current application for this optimization technique is the practical implementation of wireless communication systems that exploit adaptive symbol mapping to increase statistical independence over multiple packet re-transmissions. Commercial as well as military wireless systems could immediately benefit from optimum symbol mappings, since there would be fewer bit errors and hence fewer packet errors when communication channels are subjected to noise, interference or jamming. Lower packet errors lead to fewer (less vulnerable) re-transmissions and higher data throughput given fixed bandwidth. The recent emergence of
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
427
wireless sensor networks puts some more urgency into the adoption of this optimization technique. Moreover, the development of efficient Q3AP algorithms could enable the extension of mapping diversity to other communication scenarios. A prime example of such a scenario is the design of space-time codes for multi-antenna transceivers. Another example is the design of low-complexity error correction codes. Again, each optimal mapping can be pre-computed against several common classes of fading channels, so that the transmitters and receivers need only a simple lookup table to take advantage of the optimization. There are a number of directions for future work. Firstly, we require a better understanding of why the Q3AP is actually so much more difficult to solve than the QAP. Search space size is here certainly not the most important issue, since QAP instances with similar search space sizes are much more efficiently solved by exact and local search algorithms. One difficulty arising in the Q3AP for exact algorithms is that in the lower bound computations for the Q3AP, the NP-hard 3AP would need to be solved, while in the QAP case the lower bound computations involve solving linear assignment problems, which can be solved in O(n3). On the local search side, more research efforts are required to understand in more detail what features of the Q3AP cause the difficulties encountered by the SLS algorithms. Further work should also address the refinement of the algorithms. For example, a direct next target is to solve the size N = 16 instance REAL16 to determine whether the solution found by the SLS algorithms is actually optimal. For tackling this instance and even larger ones with exact algorithms, many considerations must be addressed, including computational strategies, choice of branching rules, frequency of bound computations, and the methods for data storage and retrieval. We expect that considerable improvements on the lower bounds currently in use are possible and that it is also worthwhile to explore other techniques for lower bound computations. Additionally, further improvements may result by taking into account the fact that wireless communication problems have a high degree of symmetry. For example, in phase shift keying (PSK) instances, the placement of the first signal vector is entirely arbitrary, so that a problem of size N automatically reduces to a problem of size N 1. Concerning local search, other schemes for the exploration of the neighborhoods may be worth being studied. Another direction for further research is the study of hybrid SLS algorithms that combine, for example, ILS with a tabu search local search. Computational results on some classes of QAP instances suggest that these may be a promising also for the Q3AP. We expect that this paper and our results will engender interest in the Q3AP by the research community. This in turn should help to find new and better ways to solve combinatorial optimization problems in such areas as facilities planning, layout of communication and transportation networks and efficient layout of microcircuits. The same benefits should accrue to various 0–1 mixed integer problems as those found in airline scheduling, organization of manufacturing operations and supply chain management. Acknowledgements This material is based upon work supported by the National Science Foundation under Grant Nos. 9900376 and 0400155. University of Pennsylvania Ph.D. candidate Yi-rong Zhu provided the derivation for the Lagrangean relaxation, from which the Q3AP lower bound calculations were made. Thomas Stu¨tzle acknowledges support of the Belgian FNRS, of which he is a research associate. We especially thank the anonymous referee who gave valuable suggestions to improve the presentation of the paper. References Aarts, E.H.L., Lenstra, J.K., 1997. Editors, Local Search in Combinatorial Optimization. John Wiley & Sons, Chichester, UK. Adams, W.P., Sherali, H.D., 1986. A tight linearization and an algorithm for zero–one quadratic programming problems. Management Science 32, 1274–1290. Adams, W.P., Sherali, H.D., 1990. Linearization strategies for a class of zero–one mixed integer programming problems. Operations Research 38, 217–226. Adams, W.P., Guignard, M., Hahn, P.M., Hightower, W.L., 2007. A level-2 reformulation–linearization technique bound for the quadratic assignment problem. European Journal of Operational Research 180, 13–26. doi:10.1016/j.ejor.2006.03.051. Burkard, R.E., Rendl, F., 1984. A thermodynamically motivated simulation procedure for combinatorial optimization problems. European Journal of Operational Research 17, 169–174.
428
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
Burkard, R.E., Karisch, S.E., Rendl, F., 1991. QAPLIB – A quadratic assignment problem library. European Journal of Operational Research 55, 115–119, QAPLIB is found on the web at. Samra, H., Ding, Z., 2003a. Symbol mapping diversity design for packet retransmissions through fading channels. In: Proceedings of the Globecom 2003 Conference, 1989–1993. San Francisco, CA. Samra, H., Ding, Z., 2003b. Effective ARQ protocols using adaptive modulation and symbol mapping diversity. In: Proceedings of the 37th Asilomar Conference on Signals, Systems and Computers, Monterey, CA, pp. 74–78. Samra, H., Ding, Z., 2003c. Symbol mapping diversity in iterative decoding/demodulation of ARQ systems. In: Proceedings of the 2003 International Conference on Communications, Anchorage, AL, pp. 3585–3589. Samra, H., Ding, Z., Hahn, P.M., 2003. Optimal symbol mapping diversity for multiple packet transmissions. In: Proceedings of the 2003 International Conference on Acoustics, Speech, and Signal Processing, Hong Kong, IV, pp. 181–184. Samra, H., Ding, Z., Hahn, P.M., 2005. Symbol mapping diversity design for multiple packet transmissions. IEEE Transactions on Communications 53 (5), 810–817. Sherali, H.D., Adams, W.P., 1990. A hierarchy of relaxations between the continuous and convex hull representations for zero–one programming problems. SIAM Journal on Discrete Mathematics 3, 411–430. Sherali, H.D., Adams, W.P., 1994. A hierarchy of relaxations and convex hull characterizations for mixed-integer zero–one programming problems. Discrete Applied Mathematics 52, 83–106. Sherali, H.D., Adams, W.P., 1999. A Reformulation–Linearization Technique for Solving Discrete and Continuous Non-convex Problems. first ed. Kluwer Academic Publishers, Norwell, MA. Stu¨tzle, T., 2006. Iterated local search for the quadratic assignment problem. European Journal of Operational Research 174 (3), 1519– 1539. An earlier version is available as Technical Report AIDA-99–03. FG Intellektik, FB Informatik, TU, Darmstadt, Germany 1999. Stu¨tzle, T., Dorigo, M., 1999. Ant algorithms for the quadratic assignment problem. In: Corne, D., Dorigo, M., Glover, F. (Eds.), New Ideas in Optimization. McGraw-Hill, London, UK. Taillard, E., 1991. Robust taboo search for the quadratic assignment problem. Parallel Computing 17, 443–455. Taillard, E.D., 1998. FANT: Fast ant system. Technical Report IDSIA-46–98, IDSIA, Lugano, Switzerland.
Discrete Optimization
The quadratic three-dimensional assignment problem: Exact and approximate solution methods Peter M. Hahn a,*, Bum-Jin Kim a, Thomas Stu¨tzle b, Sebastian Kanthak c, William L. Hightower d, Harvind Samra e, Zhi Ding e, Monique Guignard f a
f
Electrical and Systems Engineering, University of Pennsylvania, Philadelphia, PA 19104-6315, USA b CoDE, IRIDIA, CP 194/6, Universite´ Libre de Bruxelles, 1050 Brussels, Belgium c Computer Science, Darmstadt University of Technology, 64283 Darmstadt, Germany d Mathematics and Computer Science, High Point University, High Point, NC 27262, USA e Electrical and Computer Engineering, University of California, Davis, Davis, CA 95616, USA Operations and Information Management, The Wharton School, University of Pennsylvania, Philadelphia, PA 19104-6340, USA Received 8 November 2005; accepted 14 November 2006 Available online 11 January 2007
Abstract This paper reports on algorithm development for solving the quadratic three-dimensional assignment problem (Q3AP). The Q3AP arises, for example, in the implementation of a hybrid ARQ (automatic repeat request) scheme for enriching diversity among multiple packet re-transmissions, by optimizing the mapping of data bits to modulation symbols. Typical practical problem sizes would be 8, 16, 32 and 64. We present an exact solution method based upon a reformulation linearization technique that is one of the best available for solving the quadratic assignment problem (QAP). Our current exact algorithm is useful for Q3AP instances of size 13 or smaller. We also investigate four stochastic local search algorithms that provide optimum or near optimum solutions for large and difficult QAP instances and adapt them for solving the Q3AP. The results of our experiments make it possible to get good solutions to signal mapping problems of size 8 and 16. 2006 Elsevier B.V. All rights reserved. Keywords: Quadratic assignment; Branch-and-bound; Heuristics; Three-index assignment; OR in telecommunications
1. Introduction Pierskalla (1967) introduced the quadratic three-dimensional assignment problem (Q3AP) in a technical memorandum. The work was never published in the open literature. Since then, nothing on the subject has * Corresponding author. Present address: 2127 Tryon Street, Philadelphia, PA 19146-1228, USA. Tel.: +1 215 546 3413; fax: +1 215 546 4043. E-mail addresses: [email protected] (P.M. Hahn), [email protected] (B.-J. Kim), [email protected] (T. Stu¨tzle), [email protected] (S. Kanthak), [email protected] (W.L. Hightower), [email protected] (H. Samra), [email protected] (Z. Ding), [email protected] (M. Guignard).
0377-2217/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2006.11.014
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
417
appeared in the literature. The authors re-discovered the Q3AP while working together on a problem arising in data transmission system design. We will explain this problem, after giving the formulation of the quadratic assignment problem (QAP) and the Q3AP. The Lawler (1963) formulation of the QAP is given by ( ) N X N X N X N X C ijkn xij xkn : x 2 X; x binary ; ð1Þ min i¼1
j¼1
k¼1 n¼1
PN PN where x 2 X fx P 0 : i¼1 xij ¼ 1 for j = 1, . . . , N; j¼1 xij ¼ 1 for i = 1, . . . , N}. The Koopmans and Beckmann (1957) formulation of the QAP assumes that the four-dimensional cost matrix C is generated from a pair of N · N matrices, F and D, with Cijkn = fik Æ djn (i, j, k, n = 1, . . . , N). The Q3AP is an extension of the QAP that can be represented as a three-dimensional assignment problem and it is given by the following formulation: 9 8 > > > > =
n6¼j
q6¼p
where I, J and P are disjoint sets of constraints: ( ) N X N X I¼ xP0: xijp ¼ 1 for i ¼ 1; . . . ; N ; j¼1
( J¼
xP0:
P¼
N X N X i¼1
( xP0:
p¼1
p¼1
N X N X i¼1
ð3Þ
) xijp ¼ 1 for j ¼ 1; . . . ; N ;
ð4Þ
) xijp ¼ 1 for p ¼ 1; . . . ; N :
ð5Þ
j¼1
The problem is so named because the object is to optimize a quadratic function over the three-dimensional assignment polytope x 2 I \ J \ P. The quadratic expression in the objective function need not include any xijpxknq terms for which i = k or j = n or p = q unless all three equalities hold. If i = k and j = n and p = q, then xijpxknq = xijp; otherwise xijpxknq = 0. Being an extension of the NP-hard problems QAP and the threedimensional assignment problem (3AP) (Garey and Johnson, 1979), it is easy to see that the Q3AP is also NP-hard. The Q3AP was motivated, in part, by a problem arising in wireless data transmission systems, where information data bits for transmission are organized into fixed length packets. Just before transmission, several contiguous bits are mapped into distinct discrete symbols during modulation. At the receiver, these symbols are demodulated and mapped back into the original bit strings, which are then put back into data packets (for delivery by service providers). Due to interference, receiver noise, channel fading, and channel distortions, erroneous data packets may be received in these systems, which can be detected by the receiver through a cyclical redundancy checking (CRC) code inserted at the transmitter. To ensure data service quality, re-transmission of the erroneous packet may be requested upon a CRC failure by sending, over a feedback channel, an automatic repeat request (ARQ). Recent research works (Samra and Ding, 2003a,b,c) have shown that the mapping of binary data in the retransmitted packet into symbols should be changed for each re-transmission in order to significantly improve the joint detection of the original erroneous bit strings from multiple observations. Loosely speaking, for this ARQ mechanism to be as effective as possible, the symbol mappings for the repeated transmission should provide as much diversified representation (or statistical independence) as possible among the resulting received transmission symbols. In the parlance of communication systems, providing statistical independence among signals representing the same data is called transmission diversity. Providing such diversity between re-transmissions of the same intended data can be done on a pair-wise basis and would constitute solving a QAP. In this case, the cost coefficients {Cijkn} indicate the cost of simultaneously assigning bit strings i and k to the transmitted symbols j and n. However, better performance is achieved by making the
418
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
original transmission and the two following repeats as statistically independent as possible. This latter optimization requires the Q3AP.1 In that case, the cost coefficients {Cijpknq} refer to the simultaneous assignment of the symbols j and n to the bit strings i and k in the first re-transmission and the assignment of symbols p and q in the second re-transmission. In both cases, the cost coefficients for the QAP and the Q3AP can be computed, based on some usual assumptions on the error model of the transmission errors (corruption by additive white Gaussian noise with fixed variance). For details of these computations, we refer to (Samra et al., 2005). Note, that in the QAP case, there are N! possible mappings, while in the Q3AP case, we have N! · N! possible mappings. Fortunately, the optimization of the diversity method needs to be performed only once in the design phase of any actual implementation. There is no need for real time optimization. Once the optimal mapping is acquired for a given type of wireless fading channel, the transmitter/receiver can simply use a lookup table to determine and agree with each other on the next mapping. Results in Samra et al. (2003, 2005) demonstrate tremendous gains, even with suboptimal QAP mapping. For four common wireless fading scenarios from flat fading to frequency selective channels, their results show more than 3–4 dB of signal-to-noise ratio (SNR) gain, which implies more than 50% power savings or 100% increase in signal coverage area. Moreover, it is shown in Samra et al. (2003, 2005) that mapping diversity with only one re-transmission delivers better performance than simple ARQ with three re-transmissions. Again, there is a 50% reduction in bandwidth usage or 100% improvement in capacity. In this paper, we study algorithms for tackling the Q3AP, including exact algorithms and heuristics. Our prior computational experience in solving the QAP forms the basis for this paper. We enjoyed marked success by solving difficult QAPLIB (Burkard et al., 1991) instances that previously were unsolved, including the Krarup30a (Hahn and Krarup, 2001) and the Nugent et al. (1968) instance of size N = 25. To solve both, we implemented a bound computation (Hahn et al., 1998) based on a dual ascent procedure that efficiently solves the level-1 RLT formulation (Adams and Sherali, 1986, 1990). In Guignard et al. (2005) and Kim (2006), we extended these bounds to the Q3AP and we conducted experimental tests to determine the utility of level-1 RLT bounds for solving the Q3AP. We designed and implemented a basic branch-and-bound algorithm, and applied it to a set of benchmark instances created from difficult benchmark QAP instances found on QAPLIB. These chosen academic QAP problem instances are noted for their difficulty. Thus, the Q3AP instances derived from them are also difficult. We also applied the branch-and-bound algorithm to a size 8 Q3AP corresponding to the dual mappings required for two 8-PSK re-transmissions. This latter instance is hereafter referred to as REAL-8. Another realistic problem of size 16 was generated and attempts were made to solve it exactly. This problem (denoted REAL-16) represented the dual mappings required for two 16-QAM (quadrature amplitude modulation) re-transmissions. So far, we have experimental exact solution results to the four Q3AP instances derived from QAP, of sizes N = 8, 9, 12 and 13 and to the REAL-8 instance. Due to the eight symmetries, the REAL-8 instance is solved in less time than required to solve the Nug8-derived instance, which has only four symmetries. However, the computation times for all these instances are substantial and therefore we also studied the application of stochastic local search (SLS) methods (Hoos and Stu¨tzle, 2004). In fact, computational experience with SLS methods on the QAP has shown that they typically reach high quality or optimal solutions much quicker than exact algorithms. Here, we adapted previously proposed SLS algorithms for the QAP to solve the Q3AP, including a simulated annealing (SA) algorithm of Connolly (1990), the Robust Tabu Search algorithm of Taillard (1991), the FANT algorithm by Taillard (1998), and an iterated local search algorithm (ILS) by Stu¨tzle (2006). As shown in an experimental study of these algorithms, they can reach high quality solutions in reasonable computation times. In fact, the best performing of these methods improves the best-known objective function value for the REAL-16 instance by about 19%, when compared with the best results of other heuristic methods. In Section 2, we present a level-1 reformulation linearization technique (RLT) relaxation on the Q3AP and in Section 3 we demonstrate a lower bounding method based on level-1 RLT. RLT bounding methods have proved successful in branch-and-bound algorithms for the QAP and appear to work well on the 1
In fact, performance is improved by making the following 3, 4, . . . , X re-transmissions as statistically independent as possible, which would require Q4AP, Q5AP, . . . , QXAP optimization. However, doing this would lead to diminishing returns, at the cost of increased system complexity. Our results show that there is still a significant benefit when moving from the optimization of one re-transmission to the optimization of the two re-transmissions, that is, when moving from the QAP to the Q3AP.
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
419
Q3AP. Experimental results of our exact solution method are presented in Section 4. In Section 5, we introduce four SLS algorithms for this very difficult problem. Results of our experiments with SLS algorithms are presented in Section 6. Section 7 presents our conclusions and our plans for future work. 2. Reformulation linearization technique (RLT) on the Q3AP The RLT was first introduced in Adams and Sherali (1986) and Adams and Sherali (1990) and best described in Sherali and Adams (1999). It is a technique for generating a hierarchy of polyhedral representations for linear and polynomial 0–1 programming problems. These new representations are extended linear programs that permit the construction of Lagrangean duals that step-by-step come closer to solving exactly the original problem. This is demonstrated in Adams et al. (2007) for RLT level-2. In this paper we take advantage of only the first step (level-1) of this hierarchy. As suggested by its name, the RLT methodology consists of two steps: reformulation and linearization. We describe in this section how this methodology applies to the Q3AP, and refer to Sherali and Adams (1990, 1994, 1999) for more general and detailed explanations. Consider the Q3AP as presented in Eqs. (2)–(5). The level-1 RLT formulation is generated via the following two steps: Reformulation Multiply each of the 3N equations and each of the N3 non-negative restrictions x P 0 defining I, J, P in (3)–(5) by each binary variable xknq, and append these new restrictions to (3)–(5), respectively. Substitute xknq ¼ x2knq throughout the objective function and constraints. When a variable xijp in a given constraint is multiplied by xknq, express the resulting product as xijpxknq in that order. (Observe that for binary variables, such multiplications, together with the assignment constraints, imply that xijpxknq = 0 if i = k or j = n or p = q unless all three equalities hold.) Linearization Linearize the resulting problem by substituting, for every occurrence of each product xijpxknq with i 5 k, j 5 n, and p 5 q, the continuous variable vijpknq. Enforce the trivial but extremely important restrictions that vijpknq = vknqijp for all (i, j, p, k, n, q) with i < k, j 5 n and p 5 q. These restrictions define what we call complementary pairs. When applied to the bound computation, these have a dramatic effect on increasing the bound. The level-1 RLT formulation for the Q3AP is generated as instructed, which results in: RLT: XXX XXXXXX Minimize : bijp xijp þ C ijpknq vijpknq i
Subject to :
j
XX k k6¼i
n n6¼j
k k6¼i
q q6¼p
n n6¼j
q q6¼p
XX XX
p
i
j
p
k k6¼i
n n6¼j
ð6Þ
q q6¼p
vijpknq ¼ xijp
8ði; j; p; qÞ; q 6¼ p;
ð7Þ
vijpknq ¼ xijp
8ði; j; p; nÞ; n 6¼ j;
ð8Þ
vijpknq ¼ xijp
8ði; j; p; kÞ; k 6¼ i;
ð9Þ
vijpknq ¼ vknqijp 8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q; vijpknq P 0 8ði; j; p; k; n; qÞ; i 6¼ k; j 6¼ n; p 6¼ q; x 2 I \ J \ P; x binary:
ð10Þ ð11Þ ð12Þ
Problem Q3AP and RLT are equivalent in the following sense. Given any feasible solution x to the Q3AP, there exists a v such that (x, v) is feasible to the RLT with the same objective value. Conversely, for any feasible solution (x, v) to the RLT, the corresponding x is feasible to the Q3AP with the same objective value. The proof is straightforward and follows the standard RLT proof arguments.
420
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
One can show that the problem RLT possesses a block-diagonal structure that can be exploited. In particular, when constraints (10) are not present, RLT reduces to N3 + 1 three-dimensional assignment problems, or 3AP’s, one over each (i, j, p) triplet in (7)–(9), and one over x 2 I \ J \ P. This suggests a generalization of the Gilmore–Lawler bound (GLB) for the QAP (Gilmore, 1962). However, the dual procedure bound proposed in Hahn and Grant (1998), of which the GLB is a special case, performs significantly better. 3. A level-1 RLT lower bound for the Q3AP In order to develop a Lagrangean dual procedure to solve the RLT, we first obtain a smaller reformulation via the substitution of variables suggested by constraint (10), thereby reducing the number of variables and constraints, in addition to halving the number of non-negativity restrictions in (11). Now, the formulation becomes: RLT1: XXX XXXXXX e ijpknq vijpknq Minimize : bijp xijp þ ð13Þ C i
Subject to :
j
p
i
j
p
k k>i
e ijpknq ¼ C ijpknq þ C knqijp ; where C XX XX vijpknq þ vknqijp ¼ xijp k k>i
n n6¼j
XX k k>i
q q6¼p
n n6¼j
q q6¼p
n n6¼j
q q6¼p
XX XX
q q6¼p
8ði; j; p; qÞ; q 6¼ p;
ð14Þ
8ði; j; p; nÞ; n 6¼ j;
ð15Þ
n n6¼j
k k
vijpknq þ
n n6¼j
XX k k
vknqijp ¼ xijp
q q6¼p
vijpknq ¼ xijp
8ði; j; p; kÞ; k > i;
ð16Þ
vknqijp ¼ xijp
8ði; j; p; kÞ; k < i;
ð16aÞ
8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q;
ð17Þ
vijpknq P 0
x 2 I \ J \ P; x binary:
ð18Þ
Now we consider the following Lagrangean dual problem LD1, whereby constraints (14), (15), (16), (16a) are placed into the objective function using multipliers q, k, g and u, respectively. LD1: maximize hðq; k; g; uÞ where 9 8 > > =
n6¼j
q6¼p
where ~ bijp ¼ bijp þ
X q q6¼p
qijpq þ
X n n6¼j
kijpn þ
X k k>i
gijpk þ
X
uijpk
8ði; j; pÞ;
ð20Þ
k k
e ijpknq qijpq qknqp kijpn kknqj gijpk uknqi C ijpknq ¼ C
8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q;
ð21Þ
subject to constraints (17) and (18). Next, we consider a further Lagrangean dual problem LD2, whereby equality constraints in (18) on x variables are also placed into the objective function using multipliers a, b and c, respectively.
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
421
LD2: maximize hða; b; c; q; k; g; uÞ where 9 8 > > =
ð22Þ
where bijp ¼ ~ bijp ai bj cp
8ði; j; pÞ;
ð23Þ
subject to: vijpknq P 0 0 6 xijp 6 1
8ði; j; p; k; n; qÞ; i < k; j 6¼ n; p 6¼ q;
ð24Þ
8ði; j; pÞ:
ð25Þ
Notice that we now relax the equivalent formulation by allowing xijp to be fractional between 0 and 1. 4. Level-1 RLT branch-and-bound results To demonstrate the potential of the approach, we implemented a branch-and-bound algorithm that uses the level-1 RLT bounds. Our branch-and-bound algorithm follows the same strategies as those utilized for the level-1 RLT branch-and-bound algorithm for solving large QAP instances (Hahn et al., 2001). Table 1 summarizes our experimental results on five Q3AP instances of size 8, 9, 12, and 13. One of these instances (designated REAL-8) is a size 8 symbol-mapping problem (or the 8 phase shift keying modulation known as 8-PSK), while the other four instances were derived from QAP instances. Three of these instances were created from the now famous Nugent problem instances (Nugent et al., 1968) found in QAPLIB (Bur´ ric Taillard. This latter kard et al., 1991). One instance, ‘‘Tai9’’ was derived from an instance posed by E instance involves the placement of 10 incrementally heavier blades on the periphery of a turbine shaft. Since the placement of the first blade is arbitrary, this instance reduces to a QAP of size 9. All these QAP instances are Koopmans–Beckmann. (For a definition of Koopmans–Beckmann, see the first paragraph of Section 1.) To generate a six-dimensional cost matrix for one of our Q3AP instances, we used the following relationship Cijpknq = fik Æ djn Æ fik Æ dpq; (i, j, p, k, n, q = 1, . . . , N). One might wonder about the validity of testing a general Q3AP algorithm on problem instances that are Koopmans–Beckmann. In fact, our experimental results show that the Koopmans–Beckmann generated problems are harder to solve than non-Koopmans–Beckmann problems of the same size arising from wireless communication examples. The size 8, 9 and 12 problem instances were run on a Sun Ultra 10 workstation with a single 360 MHz CPU. The size 13 problem instance, derived from the Nug15 instance, was run on a single 733 MHz CPU of a Dell 7150 server. The runtimes in Table 1 are normalized for the Dell server. The Nugent derived Q3AP test problems required significantly more time to solve than the QAP instances from which they were created, with the size 8 instance taking longer than it would to exhaustively enumerate Table 1 Experimental results on Level-1 RLT branch-and-bound algorithm
Root lower bound Optimum value No. of optima Starting upper bound Dell 7150 runtime Number of nodes a
8-PSK
Nug8 derived
Tai9 derived
Nug12 derived
Nug13a derived
5,909,755,156 15,328,457,328 6144 9E + 18 155.8 seconds 267,264
91.5 134.0 2 9E + 18 48.5 seconds 621
55,250,527,200,558 74,755,545,715,200 32 9E + 18 10.3 minutes 12,666
229.8 580.0 3 9E + 18 14.8 hours 65,597
899.7 1912.0 2 9E + 18 23.1 days 1,790,840
Cutting down the flow and distance matrices of the Nug15 problem created the Nug13 instance.
422
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
the 8! · 8! feasible solutions (at 10 nanoseconds per solution). However, the 14.8 hours it takes to solve the size 12 problem instance is more reasonable, considering that exhaustive enumeration would take over 70 years. 5. Approximate solution methods Stochastic local search methods belong to the most successful techniques for the approximate solution of hard optimization problems (Aarts and Lenstra, 1997; Hoos and Stu¨tzle, 2004). These methods are also the best way to provide feasible solutions in large industrial problems with tight constraints. In this work, we adapted four successful SLS algorithms for the QAP to the Q3AP: the simulated annealing algorithm by Connolly (1990), the robust tabu search (RoTS) algorithm and fast ant system (FANT) by Taillard (1991, 1998), and an iterated local search (ILS) algorithm by Stu¨tzle (2006). Computational experience has shown that these algorithms can find optimal solutions to a wide variety of QAP instances much faster than the best performing exact algorithms. For example, the ILS algorithm finds the best-known solution to the above mentioned Nugent 30 instance in about 1.3 CPU seconds, on average (measured on a 1.2 GHz MP Athlon CPU). These results on the QAP suggest that if one wants to solve the Q3AP with sizes as large as N = 64, these methods are probably the only means to get good quality solutions. Although, these methods cannot prove optimality of the solutions, they can be useful in providing good initial upper bounds for exact methods. 5.1. Iterative improvement for the Q3AP All four SLS algorithms we implemented rely on an underlying iterative improvement algorithm. Local search algorithms for the QAP usually exploit a solution representation based on permutations. In fact, an alternative representation of the general QAP formulation in Eq. (1) is given by min
N X N X i¼1
ð26Þ
C i/ðiÞk/ðkÞ ;
k¼1
where / is a permutation of the set of integers {1, . . . , N}. ForP QAP instances having the Koopmans–Beckman N PN property, this would result in the minimization of the term i¼1 k¼1 F ik D/ðiÞ/ðkÞ , where Cijkl = FikDjl. Similarly, a solution to the Q3AP, using the formulation given by Eqs. (2)–(5), can be represented using two permutations / and u of the set of integers {1, . . . , N} as ( ) N N X N X X min f ð/; uÞ ¼ bi/ðiÞuðiÞ þ C i/ðiÞ/ðjÞjuðiÞuðjÞ : /; u are permutations over the set f1; . . . ; N g : i¼1
i¼1
j¼1
ð27Þ Given this formulation, it is also easy to see that the search space size of the Q3AP is N! · N! Iterative improvement algorithms for the QAP typically make use of the 2-exchange neighborhood relation, in which two solutions are neighbored if their permutations differ in exactly two positions. That is, the neighborhood N(/) of permutation / is given by 0
0
0
N ð/Þ ¼ f/0 j/ðrÞ ¼ /ðsÞ; /ðsÞ ¼ /ðrÞ; r 6¼ s and /ðiÞ ¼ /ðiÞ 8i 6¼ r; sg:
ð28Þ
For the Q3AP, one can apply this 2-exchange neighborhood to each of / and u. Hence, two solutions for the Q3AP are neighbored, if either / or u differ in exactly two positions and we have that N ð/; uÞ ¼ fð/0 ; u0 Þj/ðrÞ0 ¼ /ðsÞ; /ðsÞ0 ¼ /ðrÞ; r 6¼ s and /ðiÞ0 ¼ /ðiÞ8i 6¼ r; s and uðiÞ0 ¼ uðiÞ8i or uðrÞ0 ¼ uðsÞ; uðsÞ0 ¼ uðrÞ; r 6¼ s and uðiÞ0 ¼ uðiÞ8i 6¼ r; s and /ðiÞ0 ¼ /ðiÞ8ig:
ð29Þ
Essential to the performance of iterative improvement algorithms can also be the strategy that is used to explore the neighborhood. The most important of these are the first-improvement strategy that accepts an improving neighboring solution as soon as it is found, and the best-improvement strategy that first examines
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
423
the whole neighborhood and then chooses a neighboring solution that gives the largest improvement in the objective function value. For the Q3AP, one may distinguish additional strategies for exploring the neighborhood, depending on whether / and u are explored in sequence or in parallel when searching for an improved neighboring solution. Finally, it has to be mentioned that the evaluation of neighboring solutions also needs to be adapted from the QAP case, which can be done in a rather straightforward way. The necessary adaptations of the four SLS algorithms we employ here concern mainly the details of the local search, as outlined above. 5.2. Simulated annealing (SA) for the Q3AP In SA, worsening moves are accepted probabilistically by a criterion that takes into account the deterioration incurred in the objective function value and a time varying parameter T called temperature. We adapted the simulated annealing of Connolly (1990), which is known to be one of the best performing SA algorithms for the QAP, to the Q3AP. In this algorithm, the neighborhood is scanned in a fixed order. A neighboring solution s 0 is accepted deterministically, if it is better than or equal to the current one, s, or, if s 0 is worse, with a probability given by exp[(f(s 0 ) f(s))/T], where f is the objective function (Metropolis et al., 1953). In our SA algorithm (SA-Q3AP), we use the same parameter settings as that of the best performing SA variant identified by Connolly for the QAP. At each iteration, SA-Q3AP generates two neighboring solutions: One, which differs in the position of two elements in permutation / and one, which differs in the position of two elements in permutation u. Then, the better of the two neighbors is used for the standard acceptance test of the SA algorithm. Preliminary experiments showed that this strategy for examining the neighborhood gives best performance compared to generating only one neighbor at each iteration or when evaluating all neighbors and then using the best for the acceptance test. In preliminary experiments, SA-Q3AP was shown to outperform an adaptation of the Simulated Annealing algorithm of Burkard and Rendl (1984) to the Q3AP. 5.3. Tabu search for the Q3AP In tabu search, escapes from local optima are possible due to the prohibition of moves to recently visited solutions. This is done by keeping track of parts of the search history, mainly in the form of attributes of recently performed moves between neighboring solutions. The tabu search algorithm we use (RoTS-Q3AP) is a straightforward extension of Taillard’s robust tabu search algorithm (Taillard, 1991). RoTS-Q3AP first evaluates all possible neighboring solutions and then moves to the best ‘‘non-tabu’’ neighbor. A neighboring solution is considered ‘‘tabu’’, if during the last tt iterations we had that p(k) = i and p(l) = j and that in the current solution p(k) 5 i and p(l) 5 j, where p = /, u and tt is a parameter called tabu tenure. In other words, it is forbidden to assign two positions in either of the two permutations the same value they had in the previous tt iterations. Additionally, if for r iterations we have p(k) 5 i, i.e., the value i was not assigned to position k in the previous r iterations, all solutions that do not assign value i to position k are forbidden; this actually forces an assignment of the value i to p(k) (forced moves). The parameter tt was set as originally proposed in the robust tabu search algorithm, while, based on some preliminary experiments, r was set to 10n3. Additionally, an aspiration criterion is used that overrides the tabu status of a solution, if it is better than the previously best-found solution. 5.4. Ant colony optimization for the Q3AP Ant algorithms and, in particular, Ant Colony Optimization (ACO) (Dorigo and Stu¨tzle, 2004) are algorithmic techniques that take inspiration from real ants behavior. In ACO, artificial ants stochastically construct solutions biased by (artificial) pheromone trails, which serve as a numerical information that is modified at the algorithm’s execution time to reflect the ants’ search experience. We applied fast ant system (FANT), which was previously studied for the QAP (Taillard, 1998). For an overview of ant algorithm applications to the QAP we refer to Stu¨tzle and Dorigo (1999). As in all other applications of ant algorithms to the QAP, the pheromone value sij refers to setting p(i) = j. In our adaptation of FANT, FANT-Q3AP, we use two pheromone matrices, one for generating each of the two permutations / and u. For the solution construction, the same parameters as proposed in Taillard (1998) are used. The only new feature we use is analogous to the
424
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
‘‘forced moves’’ of robust tabu search: If for r iterations p(k) 5 i, this assignment is deterministically adopted in the solution construction. As in RoTS-Q3AP, a good parameter setting was found to be r = 10n3. Once a solution is constructed, it is improved using a straightforward adaptation to the Q3AP of the local search procedure originally proposed for FANT. 5.5. Iterated local search for the Q3AP Iterated local search (ILS) is a conceptually simple SLS method that builds a chain of local optima by cycling through phases of solution perturbations, which are applied to local optima, local search, and acceptance test for defining the incumbent solution (Lourenco et al., 2002). For the Q3AP, we implemented a variant of the ILS algorithms studied by Stu¨tzle (2006). The resulting ILS-Q3AP algorithm works as follows. It first generates two random permutations for / and u and applies a first-improvement local search that systematically examines the 2-exchanges in / and u. The main loop of the ILS algorithm then first applies a perturbation, consisting of a random exchange of a number of positions in each of the two permutations, to a local optimum (/*, u*). This results in a perturbed solution (/ 0 , u 0 ) that is the starting solution for the next local search. Once a new locally optimal solution (/ 0 *, u 0 *) is reached, an acceptance criterion decides whether the next perturbation is applied to (/ 0 *, u 0 *) or to the previous local optimum (/*, u*). In ILS-Q3AP, worse solutions are accepted with a probability given by exp{(f((/ 0 *, u 0 *)) f((/*, u*)))/T}, where T is a parameter that is set such that a decrease of the objective function of the current solution (/*, u*) by 1% is accepted with a probability of 0.5. If for a large number of iterations (typically a few thousand) of the main ILS loop, no improved solution was found; then the search is restarted from a new, randomly generated initial solution. 6. Approximate solution method results We tested the four SLS methods on all Q3AP instances of Section 4 plus an additional one of size 15, which was derived from the Nug15 instance from QAPLIB. For each of the four algorithms, 8 runs were performed on a single 440 MHz CPU of an HPJ5000 machine, which has almost exactly the same speed for our codes as the Dell 7150 machine on which the branch-and-bound algorithm was run. On the QAP derived instance of size 8 and 9, all the algorithms found the optimal solution within a few seconds of computation time, and the results are not shown explicitly here. The experimental results on the instances of size 12–15 are summarized in Tables 2–4. In each of these tables is given the average (over 8 runs) objective function value found by the four SLS algorithms after 1 minute, 10 minutes, 1 hour and 10 hours of computation time. The best average objective value for each time limit is in bold face. For the derived Nug12 and Nug13 instances, several of the SLS algorithms can reach the optimal solution. Overall, ILS-Q3AP appears to be the best-performing algorithm. It is the only heuristic to find, in every single run, the optimum solution for the size 8, 9, 12 and 13 instances and the best-known solution for the size 15 instance. It is also interesting to compare the performance of the four SLS algorithms to the branch-and-bound algorithm. Although SLS algorithms cannot prove optimality of the solutions they find, one can compare the computation times the different algorithms need to reach the first optimal solution. To do so, we ran ILS-Q3AP and RoTS-Q3AP 25 times and measured the average time, as well as the standard deviation (SD), each
Table 2 Average objective function values for SA-Q3AP, FANT-Q3AP, TS-Q3AP, and ILS-Q3AP for eight runs on Nug12 using different time limits from 1 minute to 10 hours
FANT-Q3AP RoTS-Q3AP SA-Q3AP ILS-Q3AP a
1 minute run
10 minute run
1 hour run
10 hour run
878 779 771 796
857 748 705 634
808 607 626 589
750 580a 583 580a
Optimum objective function value.
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
425
Table 3 Average objective function values for SA-Q3AP, FANT-Q3AP, TS-Q3AP, and ILS-Q3AP for eight runs on Nug13 using different time limits from 1 minute to 10 hour
FANT-Q3AP RoTS-Q3AP SA-Q3AP ILS-Q3AP a
1 minute run
10 minute run
1 hour run
10 hour run
2322 2118 2102 2120
2130 1995 2024 2014
2107 1918 1967 1962
2054 1912a 1912a 1912a
Optimum objective function value.
Table 4 Average objective function values for SA-Q3AP, FANT-Q3AP, TS-Q3AP, and ILS-Q3AP for 8 runs on Nug15 using different time limits from 1 minute to 10 hours
FANT-Q3AP RoTS-Q3AP SA-Q3AP ILS-Q3AP
1 minute run
10 minute run
1 hour run
10 hour run
3172 2843 2752 2751
3019 2642 2631 2568
2943 2547 2504 2399
2703 2324 2400 2230
algorithm took to reach the optimal solution. We ran one run of the branch-and-bound algorithm, with the upper bound set to an extremely large value and reported the time it took to reach the first optimum. The results are presented in Table 5. From the computational results in Table 5, one can observe that the SLS algorithms reach optimal solutions about one to two orders of magnitude faster than the branch-and-bound algorithm for the instance sizes tested. We conjecture that a simple combination of the two, where an SLS algorithm is used to determine good upper bounds for the branch-and-bound, is highly desirable. In this context, it is interesting to note that the time-to-completion of the branch-and-bound algorithm was between a factor of 2 (Nug12 derived instance), to a factor of 12 (Tail9 derived instance) larger than the time to find the first occurrence of an optimal solution. Hence, if this factor is small, then such a combination can also speed up significantly the time to prove optimality. The excellent performance of SLS algorithms is confirmed on a Q3AP symbol-mapping instance (REAL16) of size 16. For REAL16, the matrix entries were originally given as double precision floating-point numbers but then scaled to truncated integers for running the SLS algorithms. The results for REAL16 are given in Table 6, again measured on a single 440 MHz CPU of an HPJ5000 machine. Best solutions are shown in bold face entries. Second best solutions are shown in italics. As can be seen, the ILS algorithm gives the best performance, if 10 minutes or more are allowed. In fact, for the largest allowed computation time of 10 hours, the worst solution identified by ILS-Q3AP across eight runs was even better than the best solution found by either Table 5 Runtime to reach optimum for ILS-Q3AP, RoTS-Q3AP and b-and-b algorithms for four Q3AP instances (run-times of ILS-Q3AP and RoTS-Q3AP are across 25 independent runs) Real size 8
Tai9 derived
NUG12 derived
NUG13 derived
ILS-Q3AP
Average time SD time Best time Worst time
0.043 seconds 0.023 seconds 0.01 seconds 0.09 seconds
1.65 seconds 1.31 seconds 0.2 seconds 5.13 seconds
2162.2 seconds (36.03 minutes) 2138.07 seconds 53 seconds 8404.57 seconds
7519.2 seconds (2.09 hours) 6801.27 seconds 69.44 seconds 27194.4 seconds
RoTS-Q3AP
Average time SD time Best time Worst time
0.44 seconds 0.58 seconds 0.01 seconds 2.57 seconds
4.06 seconds 2.86 seconds 0.04 seconds 14.64 seconds
9104.3 seconds (2.53 hours) 8107.64 seconds 169.58 seconds 35733.93 seconds
3002.0 seconds (50.03 minutes) 2752.14 seconds 128.24 seconds 9719.47 seconds
B&B
Time to 1st Opt
19.2 seconds
50.6 seconds
10.6 hours
165.35 hours
426
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
Table 6 Average objective function values and standard deviation for SA-Q3AP, TS-Q3AP and ILS-Q3AP for eight runs on REAL16 for different time limits from 1 minute to 10 hours
ILS-Q3AP RoTS-Q3AP SA-Q3AP
AVE SD AVE SD AVE SD
1 minute run
10 minute run
30 minute run
1 hour run
10 hours run
4.13730E+15 4.65885E + 14 3.61907E + 15 5.13961E + 14 3.40559E + 15 1.08742E + 14
2.88551E + 15 6.07584E + 14 3.17127E + 15 2.84639E + 14 3.04805E + 15 1.44465E + 14
2.66533E + 15 4.16076E + 14 2.92033E + 15 6.75849E + 13 2.82990E + 15 3.49652E + 14
2.46528E + 15 4.08002E + 14 2.86702E + 15 6.01477E + 13 2.73253E + 15 3.04737E + 14
2.20841E + 15 6.81110E + 13 2.73298E + 15 8.41498E + 13 2.58397E + 15 1.47923E + 14
The Quadratic Three-dimensional Assignment Problem: 4.5E+15 4E+15 3.5E+15 3E+15 2.5E+15 2E+15 1.5E+15 1E+15 5E+14 0
ILS I_SA TABU
1 min
10 min 30 min
1 hour 10 hour
Time
Fig. 1. Development of the average objective function value for ILS-Q3AP, SA-Q3AP, and RoTS-Q3AP applied to the hard, symbolmapping Q3AP instance of size N = 16.
SA-Q3AP or RoTS-Q3AP. The computational results are further illustrated in Fig. 1. It is noteworthy to mention that, after only 1 hour, the best result found by ILS-Q3AP is about 19% better than the previously bestknown solution. 7. Conclusions and future directions We present a new development of exact and local search algorithms for solving the Q3AP. These algorithms were adapted from their QAP counterparts. While for the local search algorithms these adaptations could be done in a reasonably straightforward way, for the exact algorithms this required the development and implementation of effective lower bounds for the 3AP. Although our computational results are encouraging, they also illustrate the level of difficulty associated with the Q3AP. Presently, our exact solution method, based upon the best techniques available for solving the QAP, is computationally feasible only for Q3AP instances of size 13 or smaller. Stochastic local search techniques, which provide optimum or near optimum solutions for large and difficult QAP instances, are essential for reaching high quality solutions to reasonable size Q3AP instances. Among the various SLS algorithms tested, the overall best performance was given by an iterated local search algorithm. However, thus far, even the best performing SLS algorithm incurs for small Q3AP instances computation times for reaching optimal or best-known solutions that are orders of magnitude larger than for similar size QAP instances. Clearly, much more work is needed on this challenging and yet important new combinatorial optimization problem. The current application for this optimization technique is the practical implementation of wireless communication systems that exploit adaptive symbol mapping to increase statistical independence over multiple packet re-transmissions. Commercial as well as military wireless systems could immediately benefit from optimum symbol mappings, since there would be fewer bit errors and hence fewer packet errors when communication channels are subjected to noise, interference or jamming. Lower packet errors lead to fewer (less vulnerable) re-transmissions and higher data throughput given fixed bandwidth. The recent emergence of
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
427
wireless sensor networks puts some more urgency into the adoption of this optimization technique. Moreover, the development of efficient Q3AP algorithms could enable the extension of mapping diversity to other communication scenarios. A prime example of such a scenario is the design of space-time codes for multi-antenna transceivers. Another example is the design of low-complexity error correction codes. Again, each optimal mapping can be pre-computed against several common classes of fading channels, so that the transmitters and receivers need only a simple lookup table to take advantage of the optimization. There are a number of directions for future work. Firstly, we require a better understanding of why the Q3AP is actually so much more difficult to solve than the QAP. Search space size is here certainly not the most important issue, since QAP instances with similar search space sizes are much more efficiently solved by exact and local search algorithms. One difficulty arising in the Q3AP for exact algorithms is that in the lower bound computations for the Q3AP, the NP-hard 3AP would need to be solved, while in the QAP case the lower bound computations involve solving linear assignment problems, which can be solved in O(n3). On the local search side, more research efforts are required to understand in more detail what features of the Q3AP cause the difficulties encountered by the SLS algorithms. Further work should also address the refinement of the algorithms. For example, a direct next target is to solve the size N = 16 instance REAL16 to determine whether the solution found by the SLS algorithms is actually optimal. For tackling this instance and even larger ones with exact algorithms, many considerations must be addressed, including computational strategies, choice of branching rules, frequency of bound computations, and the methods for data storage and retrieval. We expect that considerable improvements on the lower bounds currently in use are possible and that it is also worthwhile to explore other techniques for lower bound computations. Additionally, further improvements may result by taking into account the fact that wireless communication problems have a high degree of symmetry. For example, in phase shift keying (PSK) instances, the placement of the first signal vector is entirely arbitrary, so that a problem of size N automatically reduces to a problem of size N 1. Concerning local search, other schemes for the exploration of the neighborhoods may be worth being studied. Another direction for further research is the study of hybrid SLS algorithms that combine, for example, ILS with a tabu search local search. Computational results on some classes of QAP instances suggest that these may be a promising also for the Q3AP. We expect that this paper and our results will engender interest in the Q3AP by the research community. This in turn should help to find new and better ways to solve combinatorial optimization problems in such areas as facilities planning, layout of communication and transportation networks and efficient layout of microcircuits. The same benefits should accrue to various 0–1 mixed integer problems as those found in airline scheduling, organization of manufacturing operations and supply chain management. Acknowledgements This material is based upon work supported by the National Science Foundation under Grant Nos. 9900376 and 0400155. University of Pennsylvania Ph.D. candidate Yi-rong Zhu provided the derivation for the Lagrangean relaxation, from which the Q3AP lower bound calculations were made. Thomas Stu¨tzle acknowledges support of the Belgian FNRS, of which he is a research associate. We especially thank the anonymous referee who gave valuable suggestions to improve the presentation of the paper. References Aarts, E.H.L., Lenstra, J.K., 1997. Editors, Local Search in Combinatorial Optimization. John Wiley & Sons, Chichester, UK. Adams, W.P., Sherali, H.D., 1986. A tight linearization and an algorithm for zero–one quadratic programming problems. Management Science 32, 1274–1290. Adams, W.P., Sherali, H.D., 1990. Linearization strategies for a class of zero–one mixed integer programming problems. Operations Research 38, 217–226. Adams, W.P., Guignard, M., Hahn, P.M., Hightower, W.L., 2007. A level-2 reformulation–linearization technique bound for the quadratic assignment problem. European Journal of Operational Research 180, 13–26. doi:10.1016/j.ejor.2006.03.051. Burkard, R.E., Rendl, F., 1984. A thermodynamically motivated simulation procedure for combinatorial optimization problems. European Journal of Operational Research 17, 169–174.
428
P.M. Hahn et al. / European Journal of Operational Research 184 (2008) 416–428
Burkard, R.E., Karisch, S.E., Rendl, F., 1991. QAPLIB – A quadratic assignment problem library. European Journal of Operational Research 55, 115–119, QAPLIB is found on the web at