Computers & Industrial Engineering 36 (1999) 409±426 www.elsevier.com/locate/compindeng
A new hybrid optimization algorithm L. Shia,*, S. OÂ lafssonb, Q. Chena
aDepartment of Industrial Engineering, University of Wisconsin-Madison, Madison, WI 53706, USA bDepartment of Industrial and Manufacturing Systems Engineering, Iowa State University, Ames, IA 50011, USA
Abstract
We develop a new optimization algorithm that combines the genetic algorithm and a recently proposed global optimization algorithm called the nested partitions method. The resulting hybrid algorithm retains the global perspective of the nested partitions method and the local search capabilities of the genetic algorithm. We also present a detailed application of the new algorithm to a NP-hard product design problem and it is found empirically to outperform a pure genetic algorithm implementation, particularly for large problems. 5 1999 Elsevier Science Ltd. All rights reserved.
Keywords: Combinatorial optimization; Genetic algorithms; Nested partitions method; Product design
1. Introduction
Many industrial engineering applications involve combinatorial optimization, that is, obtaining the optimal solution among a ®nite set of alternatives. Such optimization problems are notoriously dicult to solve. One of the primary reasons is that in most applications the number of alternatives is extremely large and only a fraction of them can be considered within a reasonable amount of time. As a result, heuristic algorithms are often applied in combinatorial optimization. Among the most successful such heuristics are evolutionary algorithms, genetic algorithms, tabu search, and neural networks. All of these algorithms are sequential in the sense that they move iteratively between single solutions or sets of solutions. However, in some applications it may be desirable to maintain a more global perspective, that is, to consider the entire solution space in each iteration. In this paper we propose a new
* Corresponding author.
E-mail address: leyuan@ie.engr.wisc.edu (L. Shi)
0360-8352/99/$ – see front matter 5 1999 Elsevier Science Ltd. All rights reserved. PII: S0360-8352(99)00140-0
410 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
optimization algorithm to address this dicult class of problems. The new method converges to a global optimum of combinatorial optimization problems in ®nite time and is very robust. We also present a detailed application of the method to a product design problem. Numerical results demonstrate the eectiveness of our proposed method.
The new method combines two optimization algorithms: the nested partitions (NP) method and the genetic algorithm (GA). The NP method is a randomized optimization method that has recently been developed for global optimization [13]. This method has been found to be promising for dicult combinatorial optimization problems such as the traveling salesman problem [14], and production scheduling problems [12]. The NP method may be described as an adaptive sampling method that uses partitioning to concentrate the sampling eort in those subsets of the feasible region that are considered the most promising. It combines global search through global sampling of the feasible region, and local search that is used to guide where the search should be concentrated. The method is ̄exible in that it can be applied to a variety of optimization problems, and it can also be combined with other methods.
The genetic algorithm (GA) is also a random search method, which is based on the concept of natural selection. It starts from an initial population and then uses a mixture of reproduction, crossover, and mutation to create new, and hopefully better, populations. The GA usually works with an encoded feasible region, where each feasible point is represented as a string, which is commonly referred to as a chromosome. Each chromosome consists of a number of smaller strings called genes. The reproduction involves selecting the best, or ®ttest, chromosomes from the current population, the crossover involves generating new chromosomes from this selection, and ®nally mutation is used to increase variety by making small changes to the genes of a given chromosome. The GA was originally proposed in [8] and has in recent years been found to be very eective for a variety of industrial engineering problems [3]. Furthermore, various hybrid genetic algorithms have sometimes been found to outperform pure genetic algorithms [2,5,11]. Although it works with a population of strings rather than a single string, the GA is essentially a sequential heuristic and is not guaranteed to converge to a global optimum. The NP method, on the other hand, maintains a global perspective in every iteration and converges in ®nite time to a global optimum [13]. The hybrid algorithm proposed in this paper retains the bene®ts of both methods. For more detailed descriptions of the NP and GA methods we refer the reader to [13] and [4], respectively.
The remainder of the paper is organized as follows. In Section 2 we develop the new hybrid algorithm and state the basic convergence results. In Section 3 we provide a detailed application of the algorithm to a product design problem, and Section 4 contains some concluding remarks and future research directions.
2. Algorithm development
The problem that we want to solve is to maximize an objective function f: R R, over a
feasible region , that is, to solve:
maxfy: 1 y2Y
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 411
We assume that is ®nite and has been encoded in a form suitable for genetic algorithms, such as the binary string form. In theory, this problem can be solved by enumeration of the feasible region, but for most realistic applications this is a notoriously dicult problem, characterized by the combinatorial explosion of the size of feasible region, and lack of structure that can be used to guide the search. It is for such unstructured large problems that both the NP method and the GA are particularly well suited, and the same holds for the hybrid algorithm introduced here.
2.1. The nested partitions method
As stated in the introduction, the nested partitions (NP) method is a new optimization approach that has been proposed for solving global optimization problems [13]. This method can be brie ̄y described as follows. In each iteration we assume that we have a region s
that is considered the most promising. We partition this most promising region into M subregions and aggregate the entire surrounding region 7s into one region. At each iteration, we therefor look at M + 1 disjoint subsets that cover the feasible region. Each of these M+1 regions is sampled using some random sampling scheme, and the estimated performance function values at randomly selected points are used to estimate the promising index for each region. This index determines which region becomes the most promising region in the next iteration. If one of the subregions is found to be best this region becomes the most promising region. If the surrounding region is found to be best the method backtracks to a larger region. To choose this larger region we use a ®xed backtracking rule. The new most promising region is then partitioned and sampled in a similar fashion.
The methodology described above may be divided into four main steps that constitute the NP method. Each of these steps can be implemented in a generic fashion, but can also be combined with other optimization methods and adapted to take advantage of any special structure of a given problem.
1. Partitioning. The ®rst step of each iteration is to partition the current most promising region into several subregions and aggregate the surrounding region into one region. The partitioning strategy imposes a structure on the feasible region and is therefore very important for the speed of convergence of the algorithm. If the partitioning is such that most of the good solutions tend to be clustered together in the same subregions, it is likely that the algorithm quickly concentrates the search in these subsets of the feasible region. It should be noted that since the feasible region is ®nite the partitioning can be done by grouping arbitrary points together in each subregion. Therefore, a good partitioning strategy always exists, although it may not be easy to identify.
2. Random sampling. The next step of the algorithm is to randomly sample from each of the subregions and from the aggregated surrounding region. This can be done in almost any fashion. The only condition is that each solution in a given sampling region should be selected with a positive probability. Clearly uniform sampling can always be used. However, it may often be worthwhile to incorporate special structures into the sampling procedure. The aim of such a sampling method should be to select good solutions with a higher probability than poor solutions.
412 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
3. Calculation of promising index. Once each region has been sampled the next step is to use the sample points to calculate the promising index of each region. Again the NP methodology oers a great deal of ̄exibility. The only requirement imposed on a promising index is that it should agree with the original performance function on singleton regions.
4. Backtracking. If one of the subregions has the best promising index, the algorithm moves to this region and considers it to be the most promising region in the next iteration. If the surrounding region has the best promising index the algorithm backtracks to a larger region.
The NP method described here can be applied, in its generic form, to a wide range of problems [12±14]. It is also capable of taking advantage of special structures by incorporating them into the partitioning, sampling, and promising index steps. The partitioning approach also makes the algorithm very compatible with emerging parallel computing capabilities. Each region can be treated independently and in parallel, with only a minor coordination overhead. The method also uniquely combines global search through partitioning and sampling, and local search through estimation of the promising index.
2.1.1. Example (partitioning)
As an illustration of the NP method, consider the following simple example. Assume that ={1, 2, 3, 4} contains only four points. Then the feasible region can for example be partitioned into M = 2 subregions {1, 2} and {3, 4} as shown in Fig. 1. These two regions can then each be partitioned further into two subregions, which then are singleton regions and cannot be partitioned further. This creates a partitioning tree as Fig. 1 illustrates. The NP method moves from one region to another in this tree, considering in each iteration one of these regions to be most promising. q
We will combine the NP method with the GA as follows. First recall that when applying GA each solution is assumed to have a given ®tness value, and the GA search starts with an initial population and imitates natural selection to improve the ®tness of the population with the overall goal of ®nding the ®ttest solution. Therefore, we let the promising index of a region be the ®tness value of the ®ttest solution in the region and we use GA search to estimate the
Fig. 1. Partitioning tree for the NP method.
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 413
promising index, that is, to seek the ®ttest solution in each region. The random sampling step then becomes equivalent to obtaining an initial population for the GA search and the ®nal population is used for the promising index value. The partitioning and backtracking steps remain unchanged. This is a natural combination of the two methods and the resulting hybrid algorithm retains the global perspective of the NP method and the powerful local search capabilities of the GA.
2.2. Hybrid NP/GA algorithm
We now describe the hybrid NP/GA algorithm in detail. In the kth iteration we assume that there is a subregion s(k)
of the feasible region that is considered to be the ®ttest. Initially we assume that nothing is known about the ®ttest region so s(0)=. The ®rst step is to partition the ®ttest region into M disjoint subregions and aggregate the surrounding region 7s(k) into one region. Hence, in each iteration we consider a partitioning of the entire feasible region. The second step is to use some randomized method to obtain an initial population of solutions or chromosomes from each region. This should be done in such a way that each chromosome has a positive probability of being selected in the initial population. The third step is to apply the GA search to each initial population individually. This search should be constrained to stay within the region where the initial population was obtained. The forth step is to calculate a summary statistic of the ®nal population in each region and use that as a measure of the overall ®tness of the region. This summary statistic is usually the performance of the ®ttest chromosome in the ®nal population. The ®fth and ®nal step is to determine the ®ttest region for the next iteration. If one of the subregions is estimated to have the best overall ®tness, this region becomes the ®ttest region in the next iteration. The new ®ttest region is thus nested within the last. On the other hand, if the surrounding region is found to have the best overall ®tness, the algorithm backtracks to a larger region that contains the ®ttest chromosome in the surrounding region. The partitioning continues until singleton regions are obtained and no further partitioning is possible.
A potential drawback to the pure NP method is that it may get stuck for considerable time at good but suboptimal solutions at maximum depth. In other words if s(k) is a singleton region that corresponds to a solution that is better than most but not all of the other solutions, it may take many iterations before the NP algorithm backtracks from this region. In the hybrid NP/GA algorithm we avoid this problem as follows. We note that once maximum depth is reached the surrounding region contains almost all of the feasible points. Furthermore, in the ®rst iteration the entire feasible region is sampled with equal intensity. In this step, after applying GA search to each subregion we can therefore take the best chromosome from each region and form a set of M high quality chromosomes that, since they all come from dierent regions, are assured to have diverse genetic material. These chromosomes can then be reproduced into the initial population whenever the surrounding region of a maximum depth region is sampled. This helps the algorithm to backtrack from suboptimal maximum depth regions.
To state the hybrid NP/GA algorithm more precisely we need a few de®nitions. First we need to de®ne the space of subregions on which the algorithm moves.
414 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
De®nition 1. A region constructed using the partitioning scheme described above is called a valid region given the ®xed partition. The collection of all valid regions forms a space of subsets, which we denote by . Singleton regions are of special interest and we let 0 denote the collection of all such valid regions.
Thus the hybrid NP/GA algorithm moves from one element on to another, eventually aiming at locating a singleton element, that is, an element in 0 that contains a global optimum. To keep track of the status of the algorithm it is convenient to de®ne the regions by their depth in the partitioning tree and by how they are obtained by partitioning regions of one less depth.
De®nition 2. We call the singleton regions in 0 regions of maximum depth, and more generally, talk about the depth of any region. The depth function, d: R R, is de®ned iteratively in the obvious manner, with having depth zero and so forth. We assume that the maximum depth is uniform, that is, d =d(s ) for all s 6 0.
De®nition 3. If a valid region s 6 is formed by partitioning a valid region Z 6 , then s is called a subregion of region Z, and region Z is called a superregion of region s 6 .
Using the notation and de®nitions above the hybrid NP/GA algorithm is given below.
2.3. Algorithm NP/GA
Step 0. Initialization. Set k = 0 and s(k )=.
Step 1. Partitioning. If d(s(k))$d, that is, s(k)@0, partition the ®ttest region, s(k), into Ms(k ) subregions s1(k),FFF,sMs(k ) (k). If d(s(k))=d then let Ms(k )=1 and s1(k)=s(k).
If d(s(k))$0, that is, s(k)$, aggregate the surrounding region 7s(k) into one region
sMk1 k.
Step 2. Initial population. If d(s(k )) $ d use a randomized method to obtain an initial
population of Nj strings from each of the regions sj(k ), j = 1, 2, F F F , Ms(k )+1, POPj hyj1, yj2, FFF, yjNji, j1, 2, FFF, Msk 1:
2 If d(s(k ))=d then obtain a population POP j of Nj M strings as above and let the initial
IIII
j j I0
population be POP I POP I0 [ POP0, where POP0 is the diverse high quality population found in the ®rst iteration (see Step 6). Note that we implicitly assume that Nj>M.
Step 3. GA search. Apply the GA to each initial population POP Ij individually, obtaining a ®nal population for each region sIj(k ), j = 1, 2, F F F , Ms(k )+1,
POPFj hyj1, yj2, FFF, yjNji, j1, 2, FFF, Msk 1: 3 FFF
If k=0 then go to Step 6, otherwise go to Step 4.
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 415
Step 4. Overall ®tness. Estimate the overall ®tness of the region by the performance of the ®ttest chromosome in the ®nal population. That is, the overall ®tness of each region is estimated by
Fsj max fyji, j1, 2, FFF, Msk 1: 4
F i2f1, 2, FFF, Nj g
Step 5. Update ®ttest region. Calculate the index of the region with the best overall ®tness,
j 2 arg max Fs : 5 kj
j2f1, 2, FFF, Msk1g
If more than one region is equally ®t, the tie can be broken arbitrarily, except at maximum depth where ties are broken by staying in the current ®ttest region. If this index corresponds to a region that is a subregion of s(k ), then let this be the ®ttest region in the next iteration. That is s(k + 1)=sjÃk(k ) if jÃk E Ms(k ).
If the index corresponds to the surrounding region, backtrack to a region Z 6 that is de®ned by containing the ®ttest chromosome in the surrounding region and being of depth h less than the current most promising region. Note that this ®ttest chromosome y®t is known as the argument where the minimum (4) is realized for j=Ms(k )+1. In other words, s(k + 1)=Z where
dZ dsk D 6 and
yfit 2 Z: 7
This uniquely de®nes the next most promising region. Let k=k+1. Go back to Step 1.
Step 6. Initial diverse population. Let y j {j
be the ®ttest chromosome from the jth region
8
{j arg max y ji,
F i2f1, 2, FFF, Nj g
for j = 1, 2, F F F , M. Let POP0 be the set of the ®ttest chromosome from each region
POP0 hy1{1 , F F F , yM{M i: 9
Go back to Step 4.
By Step 1 of the algorithm partitioning is performed if the current most promising region is
not of maximum depth. In practice we usually stop partitioning D>0 depth levels above the maximum depth, that is, partition if and only if d(s(k )) < d D. The reason for this is that close to the maximum depth the number of solutions in each subregion is small and can be explored eciently. We also note that due to its iterative nature, some chromosomes may be generated more than once in Step 2 of the algorithm.
As the algorithm evolves a sequence of ®ttest regions {s(k )}1k = 0 is generated. This sequence
FF
416 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
is a stochastic process and may be used to analyze the convergence of the algorithm. This topic
is brie ̄y discussed in the next section.
2.4. Convergence
Due to the systematic partitioning of the feasible region the hybrid NP/GA algorithm can be seen to converge to a global optimum. The main convergence results is stated in the following theorem.
Theorem 1. The hybrid NP/GA algorithm converges to a global optimum in ®nite time. That is, let denote the set of strings that are globally optimal. Then with probability one there exists y 6 and k < 1 such that s(k)={y}, for all kek.
Proof. Due to the random method of obtaining initial populations, and the randomness involved in the GA search, the sequence of the ®ttest region in each iteration {s(k)}1k = 0 is a stochastic process. Furthermore, since the ®ttest region in any iteration only depends on which region was the ®ttest in the last iteration, {s(k)}1k = 0 is a Markov chain. It is also clear that for any y 6 , if s(k1)={y} 6 0 for some k1>0, then since s(k1)={y} } 6 S0 6 is a global optimum and times at the maximum depth are broken by favoring the current ®ttest region, s(k )={y } for all k e k1. Hence {y } 6 0 is an absorbing state. Furthermore, if y 6 7 and s(k)={y}, then in the next iteration there is a positive probability that the ®nal population in the surrounding region contains a string from and the algorithm backtracks. Therefore, y is not absorbing. Hence, the strings in correspond to the only absorbing states of the Markov chain and the chain will eventually be absorbed in one of these states with probability one.
Finally, since all the non-absorbing states are transient each can only be visited a ®nite number of times; and since there is a ®nite number of transient states the absorbion time must be ®nite. Therefore, with probability one there exists y 6 and k < 1 such that s(k )={y }, for all kek.q
In addition to proving global convergence, the Markov structure of the hybrid NP/GA algorithm can also be used to derive results about the expected number of iterations until absorbion, that is, the expected time until the global optimum is found [13]. This in turn can be used as a basis for stopping criteria and to provide guidelines for how long the GA search should be applied in each step of the algorithm. For related results we refer the interested reader to [13].
3. Application to a product design problem
In this section we apply the hybrid NP/GA algorithm to a product design problem where the objective is to design a single product that will maximize the market share of the producer.
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 417
The problem is important to market research and needs to be addressed both when entering a new market, and for existing markets since to stay competitive a producer must constantly reevaluate the product line and ad or delete products if warranted. This is generally referred to as the product design problem [1,6,7,9,15].
The product design problem involves determining the levels of the attributes of a new or redesigned product in such a way that it maximizes a given objective function. The attributes can for example be color and shape, that may take levels such as blue, green, or red for the color, and oval, round, or rectangular for the shape. We assume that the preferences or utilities of individual customers have been elicited for each level of every attribute, for example using conjoint or hybrid conjoint analysis [7], and are stored in a consumers utility matrix. We also assume that a customer will choose a particular product if its utility is higher than that of a competing status quo product, which may be dierent for each customer. This problem has been shown to be NP-hard [10].
Mathematically, the product design problem can be stated as follows. A product has K attributes of interest, each of which can be set to Lk levels, k=1, 2,FFF, K. Without loss of generality, in all the numerical results below we assume that Lk=L is the same for each attribute. There are N customers, or market segments, that each has given utilities for every level of each attribute, and the problem is to select levels for each attribute to satisfy the producers objective. Each customer i6{1, 2,FFF, N} is assumed to have a known utility uikl associated with level l of the kth attribute, 1EkEK, 1ElELk. Given these utilities, the objective may be to maximize market share, pro®t, or customer utility; but here we focus exclusively on the market share objective.
To formulate this problem mathematically we de®ne the decision variables as x={xkl}, where
xkl 1, if attribute k is set to level l, 10 0, otherwise:
The objective is to maximize the market share and a customer i will purchase the oered product if and only if this customer's total utility for this product is greater than the total utility qi of this customer's status quo product, 1 E i E N. The total utility the ith customer obtains from this product is
XK XLk
uikl xkl, k1 l1
and this customer therefore purchases the oered product if and only if
XK XLk
uikl xkl > qi, k1 l1
so the performance function for the ith customer can be written as
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
8> ( XK XLk
>max 0, uiklxklqi
, if
fix ><>: XK XLk k1 l1 : 11
> uikl xkl qi k1 l1
0, otherwise
This performance function assigns one to all the customers that purchase the oered product and zero to all other customers. Now the product design problem may be stated as the following mathematical programming problem,
XN i1
XLk l1
xkl 0, 1, 8k, l: 14
Since this problem is NP-hard it is typically solved using some heuristic method [1,6,9], and to demonstrate the necessity of using heuristics for this problem we solve it using the CPLEX mixed integer programming (MIP) solver for small instances of N, K, and L. We compare the CPU time until the MIP solver, the pure GA, and the hybrid NP/GA ®nd the known optimum. The implementation of the hybrid NP/GA algorithm is as described in Section 3.1 below, and the parameter settings for the pure GA and the hybrid algorithm can be found in Section 3.2 below. All of the algorithms are terminated once the known optimum is found. The results can be found in Table 1 and it is evident that the heuristics outperform the MIP solver for this problem. We can clearly only expect to solve relatively small problems to optimality and for realistically sized applications the use of heuristics is imperative. Also note that the
418
maxfx x
subject to
k1l1
Kk
)
XXL
fix, 12 xkl1, k1,2,FFF,K 13
Table 1
Comparison of CPU times
Problem size Solution time (s) K L N vv CPLEX MIP
3 3 20 27 0.14
4 4 50 256 9.62
5 5 100 3125 1750.44
6 6 200 46,656 109,55.75
a Less than 0.01 s.
Hybrid NP/GA
Pure GA
a a
uiklxkl6qi
0.17
0.25
0.37 0.06 1.46 3.10
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 419
pure GA performs the best for the three smaller problems but the hybrid method performs the best for the largest problem. This is consistent with the more extensive numerical results reported in Section 3.2 below. Table 1 also shows the size vv of the feasible region, and it is clear that as K and L grow the feasible region quickly becomes very large.
To apply the GA we encode the solution space by using a binary string representation. Each string or chromosome y consists of K substrings or genes, corresponding to the K attributes. The length of the kth substring is equal to Lk, the number of levels for the kth attribute. Each of these attribute substrings consists of Lk 1 zeros and a single one, indicating the level of the attribute. In terms of the decision variables (10) above the genes are represented by strings
Ak xk1xk2 FFFxkLk,
and the chromosomes are represented by strings
y A1A2 FFFAK
x11x12 FFFx1L1x21x22 FFFx2L2 FFFxK1xK2 FFFxKLK:
15
The total length of the string is denoted by P = Kk = 1Lk, so each chromosome y is represented by a P 1 matrix. The feasible region is the collection of all chromosomes y that satisfy (15) and the constraints (13) and (14) above. It is clear from equation (15) that the size of this feasible region is vv =L1L2FFFLK, or if we assume that Lk=L is constant, vv=LK. It is evident that vv grows very rapidly in the input parameters, especially as the umber of attributes K increases.
3.1. Hybrid NP/GA algorithm for the product design problem
We now describe how the hybrid NP/GA algorithm applies to the product design problem. In particular, we need to determine how to implement the partitioning step, how to obtain an initial population from each region, and how to implement the GA search. The remaining steps, namely calculation of the overall ®tness and selection of a new ®ttest region, follow the general formulas given by Eqs. (4)±(7) in the algorithm description.
3.1.1. Partitioning
At each iteration the current ®ttest region is partitioned into Lk subregions by ®xing the
kth attribute to one of its Lk possible levels. This implies that the maximum depth is d=K. The sequence in which the attributes are selected to be ®xed is arbitrary, but may aect the performance of the algorithm. In general, one would attempt to ®x the `most critical’ attribute ®rst but to simplify the notation we assume that the attributes have been ordered in the sequence in which they are to be ®xed. This imposes no loss of generality. We can write a string y 6 in terms of its K attributes as y=[Ay1Ay2FFFAyK ], and a region s 6 of depth d(s) is therefore determined by
snyAy1Ay2FFFAyK2Y:Ayi asi,i1,2, FFF,dso 16 where
420 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 asi xi1xi2 FFFxiLi
is a ®xed value for the ith attribute, that is xijs=1 for a ®xed js, and xij=0 for all j$js. Example (partitioning). Assume that there are K = 3 attributes with L1=3, L2=2, and L3=3.
The feasible region is therefore given by
Y x11x12x13x21x22x31x32x33:x11 x12 x13 1, x21 x22 1, x31 x32 x33 1:
If the current most promising region is the entire feasible region s (k)=, then there are L1=3 subregions,
s1k ny AyAyAy 2 Y:Ay as1ko 12311
ny Ay1Ay2Ay3 2 Y:Ay1 1 0 0o
y1 0 0 x21x22x31x32x33:x21 x22 1, x31 x32 x33 1
s2k ny AyAyAy 2 Y:Ay as2ko 12311
ny Ay1Ay2Ay3 2 Y:Ay1 0 1 0o
y0 1 0 x21x22x31x32x33:x21 x22 1, x31 x32 x33 1
s3k ny AyAyAy 2 Y:Ay as3ko 12311
ny Ay1Ay2Ay3 2 Y:Ay1 0 0 1o
y0 0 1 x21x22x31x32x33:x21 x22 1, x31 x32 x33 1:
Now if in the next iteration s(k + 1)=s2(k ) then there are L2=2 subregions s1k 1 ny AyAyAy 2 Y:Ay as1k1, Ay as1k1o
nyAy1Ay2Ay32Y:Ay1 0 1 0, Ay2 1 0o y0 1 0 1 0 x31x32x33:x31 x32 x33 1
123 11 22
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 421 s2k 1 ny AyAyAy 2 Y:Ay as2k1, Ay as2k1o
nyAy1Ay2Ay32Y:Ay1 0 1 0, Ay2 1 0o y0 1 0 0 1 x31x32x33:x31 x32 x33 1:
Each of these regions can then be partitioned into L3=3 subregions, which are also maximum depth regions and cannot be partitioned further.q
j1 asj as,
123 11 22
In general, the partitioning step involves partitioning a region s 6 into Ld(s )+1 subregions, si nyAyAyFFFAy2Y:Ay asi,i1,2, FFF,ds1o, i1,2, FFF,Lds1:
The general partitioning procedure is given by the following algorithm.
Algorithm. Partitioning
Step 1. Attributes that are ®xed in s are ®xed to the same levels for each of the subregions fsj gLds1 . That is, let
12Kii
ii
for i = 1, 2,FFF, d(s) and j = 1, 2,FFF, Ld(s )+1.
Step 2. Fix the next attribute of each subregion to the corresponding level asj bj, j1, 2, FFF, Lds1
ds1
where b( j ) is a binary coding of the jth level.
We note that here the binary coding b(j) of the jth level is a string of length Ld(s)+1 that has one in the jth place and zeros everywhere else,
3.1.2. Initial population
The initial population POPj from a subregion sj 6 can be selected by uniform sampling of I
the region. We note that sj is de®ned by d(sj) of the K attributes being ®xed to levels asj asj FFFasj ,
1 2 ndsj o sj yAyAyFFFAy2Y:Ay asj, i1,2, FFF,dsj ,
12Kii
so this sampling procedure should select levels for the remaining Kd(sj) attributes. By
422 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
uniform sampling we mean that each of the Lk possible levels of the kth attribute has equal
probability of being selected. A procedure for generating a uniform sample chromosome or
string ys Ays Ays F F F Ays from a region sj 6 70 is therefore given by the algorithm below. 12K
The output of this algorithm is a single chromosome ys, so it should be repeated until the predetermined size of the initial population is reached.
Algorithm. Uniform sampling
Step 0. Let Ays asj for all iEd(s ). iij
Set k=d(sj )+1.
Step 1. Generating a uniform random variable u from
(0, 1).
Step 2. If (i 1)/LkEu < i/Lk then set the kth attribute to level i, i = 1, 2,FFF, Lk.
Ays bi, k
where as before b(i ) is the binary coding of the ith level.
Step 3. If k=K stop, otherwise let k=k + 1 and go back to Step 1. Here
(0, 1) denotes the uniform distribution on the unit interval.
Example (sampling) Assume, as in the example in Section 3.1.1, that there are K=3 attributes with L1=3, L2=2, and L3=3. We saw previously that if s(k)=, then there are L1=3 subregions
s1k ny Ay1Ay2Ay3 2 Y:Ay1 1 0 0o
s2k ny Ay1Ay2Ay3 2 Y:Ay1 0 1 0o
s3k ny Ay1Ay2Ay3 2 Y:Ay1 0 0 1o
Each of these subregions is sampled by determining the value of Ay2 and Ay3, which involves
two iterations of the above algorithm. In the ®rst iteration we generate a uniform random
variable u 6 (0, 1). If u < 12then Ay2=[1 0], and if 12 Eu < 1 then Ay2=[0 1]. Note that the
cuto point is determined by L2=2. In the second iteration we generate a new uniform
random variable u6(0, 1). If u<13 then Ay3=[1 0 0], if 13Eu<23then Ay3=[0 1 0], and if
23 Eu < 1 then Ay3=[0 0 1]. As before, the cuto points are determined by L3=3. This
procedure is identical for each of the three regions. For example, when we sample s1(k), if
u < 1 in the ®rst iteration and 1 Eu < 2 in the second iteration, then the sample chromosome 2 ysysys33
generatedisys=[A1 A2 A3 ]=[10010010].q
Alternatively, we can assign unequal weight to each level based on some heuristic reasoning. This leads to a weighted sampling strategy for obtaining the initial population. The algorithm
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 423
for weighted sampling is identical to the uniform sampling algorithm, except instead of the uniform probabilities used in Step 2, each attribute may have a dierent probability of being selected.
3.1.3. GA search
Once an initial population POPj has been obtained from a subregion sj, then the GA search I
proceeds as follows. First the Nj/2 ®ttest chromosomes are selected for reproduction from the population. These chromosomes will also survive intact in the next population. Secondly, parts of chromosomes are selected randomly from the reproduced chromosomes. Each attribute, that is not ®xed in sj, of the selected pair has a given probability of participating in a crossover, that is, being swapped with the corresponding attribute of the other chromosome. Finally, each chromosome may be selected with a given probability as a candidate for mutation. If a chromosome is mutated, an attribute from that chromosome is selected at random and assigned a random level. As before, only attributes that are not ®xed in sj may be selected to be mutated. The following algorithm gives the exact procedure for the GA search in a region sj 6 .
Algorithm. GA search
Step 0. Initialization. Let POP0=POPjI. Let k=0.
Step 1. Reproduction. Sort the
according to their ®tness
fy1 efy2 e F F F efyNj :
kkk
Add the ®ttest strings to the new population, POPk1 hy1y2FFFyNj=2i:
in the
current
population
POPk=[y1 y2 F F F yNj ] kkk
strings
kkk
Step 2. Crossover. Select two random indices 1Ei, jENj/2, i$j. Select a random attribute
k where d(sj ) < k E K, and create two new chromosomes by replacing the kth attribute of yi [j] k
with the kth attribute of yk , and vice versa. Add these new chromosomes to the new population POPk + 1. Repeat until Nj/2 new chromosomes have been created.
Step 3. Mutation. Selct a random index 1EiENj. Select a random attribute k where d(sj ) < k E K, and mutate yi by replacing its kth attribute with a random level. Repeat.
k
Step 4. Let k=k + 1. If k>kmax, the total number iterations, go to Step 5. Otherwise, go
back to Step 1.
Step 5. Let POPFj POPk.
Since the GA search algorithm is called M = 1 times for each iteration, where M is the number of subregions, the hybrid NP/GA algorithm provides an automatic parallelization of the GA search eort. Each region can be treated independently and in parallel, and the hybrid
424 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
Fig. 2. Comparison of pure GA and hybrid NP/GA for 10 and 20 attribute problems.
algorithm is therefore very compatible with parallel computation. This may be useful for very large problems.
3.2. Numerical results
In this section we present numerical results that compare the performance of the hybrid NP/ GA algorithm to a pure GA implementation. The implementation of the hybrid NP/GA algorithm is as described in the last section, with the following parameters. The initial population size for each region is taken to be 20 and it is obtained using uniform sampling. The mutation rate is 0.2, that is, a chromosome has a 20% chance of being mutated. The GA search is applied for a ®xed number of 20 iterations for each region in each step except at maximum depth where it is increased to 1000 to enable to algorithm to backtrack more easily. If backtracking is chosen in the kth iteration then the algorithm backtracks h=d(s(k ))/2 steps, or to a depth that is half of the current depth. The pure GA algorithm is implemented with the same parameters as the GA search in the hybrid algorithm but is not terminated after a ®xed number of iterations. Both the hybrid NP/GA and the pure GA are terminated after a ®xed time interval.
The results are for four dierent problems. Due to the size of these problems it is not possible to solve them exactly within reasonable time and the true optimum is therefore unknown. All the problems have N = 200 customers and for simplicity we let all the attributes have the same number of levels. The ®rst problem has K = 10 attributes and L = 10 levels for each attribute. The second problem has K = 20 attributes and L = 20 levels for each attribute. The third problem has K = 30 attributes and L = 30 levels for each attribute. The fourth problem has K = 50 attributes and L = 20 levels for each attribute. The utilities of each level of each attributes are generated uniformly for each customers. The status quo prices are also generated uniformly for each customer.
The performance of both algorithms may be found in Fig. 2 for the two smaller problems,
L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426 425
and Fig. 3 for the two larger problems. These graphs show the best objective function value found as a function of the CPU time in seconds. Due to the size of these problems it is not known if the best solutions reported here are the true optima. As is to be expected, due to increased overhead cost of the hybrid algorithm, the pure GA search obtains good solutions sooner, and for the smaller problems performs better for the ®rst 30±60 s. However, the hybrid NP/GA algorithm catches up to it very rapidly and then consistently outperforms the pure GA. The bene®ts of the hybrid algorithm increase as the problem size increases, and for the two larger problems the solution found by the hybrid algorithm outperforms the ones obtained by the pure GA immediately from the time when the ®rst solution is reported by the hybrid method. These results give a strong indication that the hybrid NP/GA algorithm may be very useful in applications were the problem size is very large.
4. Conclusions and future research
We have presented a new optimization algorithm that combines a recently proposed global optimization method called the nested partitions (NP) method and genetic algorithms (GA) in a novel way. The resulting algorithm retains the bene®ts of both methods, that is, the global perspective and convergence of the NP method and the powerful local search capabilities of the GA. This new optimization algorithm was shown empirically to be capable of outperforming pure GA search for a dicult product design problem. This is especially true for very large problems where the global perspective of the NP method is of particular importance.
However, further theoretical and empirical development is needed for the algorithm. Product design problems provide a rich class of NP-hard problems that will be used for extensive empirical testing of the algorithm. In particular, we will consider the more dicult problem
Fig. 3. Comparison of pure GA and hybrid NP/GA for 30 and 50 attribute problems.
426 L. Shi et al. / Computers & Industrial Engineering 36 (1999) 409±426
where a product line, rather than a single product, must be designed. We will also consider dierent objective functions, such as total customer utility and pro®t.
The theoretical development of the hybrid NP/GA algorithm will focus on using the Markov structure to obtain results about the time until convergence and to develop stopping rules. This is also expected to be of practical importance since good stopping rules are imperative in applications.
References
[1] Balakrishnan PV, Jacob VS. Genetic algorithms for product design. Management Science 1996;42:1105±17.
[2] Chen R, Gen M. Parallel machine scheduling problems using memetic algorithms. Computers & Industrial
Engineering 1997;33:761±4.
[3] Gen M, Wasserman GS, Smith AE. Special issue on genetic algorithms and industrial engineering. Computers
& Industrial Engineering 1996;30:000.
[4] Goldberg DE. Genetic algorithms in search, optimization, and machine learning. Reading, MA: Addison-
Wesley, 1989.
[5] Gong D, Gen M, Yamazaki G, Xu W. Hybrid evolutionary method for capacitated location-allocation
problem. Computers & Industrial Engineering 1997;33:577±80.
[6] Green PE, Krieger AM. Recent contribution to optimal product positioning and buyer segmentation. European
Journal of Operations Research 1989;41:127±41.
[7] Green PE, Srinivasan V. Conjoint analysis in consumer research: new developments and directions. Journal of
Marketing 1990;54:3±19.
[8] Holland JH. Adaptation in natural and arti®cial systems. Ann Arbor, MI: The University of Michigan Press,
1975.
[9] Kohli R, Khrisnamurti R. A heuristic approach to product design. Management Science 1987;33:1523±33.
[10] Kohli R, Khrisnamurti R. Optimal product design using conjoint analysis: computational complexity and algor- ithms. European Journal of Operations Research 1989;40:186±95.
[11] Murata T, Ishibuchi H, Tanaka H. Genetic algorithms for ̄owshop scheduling problems. Computers & Industrial Engineering 1996;30:1061±72.
[12] OÂ lafsson S, Shi L. A method for scheduling in parallel manufacturing systems with ̄exible resources, IIE Transactions 1998 (in press).
[13] Shi L, OÂ lafsson S. Nested partitions method for global optimization. Operations Research 1998 (in press).
[14] Shi L, OÂ lafsson S, Sun N. New parallel randomized algorithms for the traveling salesman problem. Computers
& Operations Research 1999;26:371±9.
[15] Zufryden F. A conjoint-measurement-based approach to optimal new product design and market segmentation.
In: Stocker AD, editor. Analytical approaches to product and market planning. Cambridge, MA: Marketing Science Institute, 1977.