Manpower allocation using genetic annealing

European Journal of Operational Research 111 (1998) 405±420

Theory and Methodology

Manpower allocation using genetic annealing Nicolas Abboud a

a,1

, Masahiro Inuiguchi a, Masatoshi Sakawa

a,*

, Yoshio Uemura

b,2

Department of Industrial and Systems Engineering, Faculty of Engineering, Hiroshima University, 4-1 Kagamiyama 1 chome, Higashi Hiroshima-shi, Hiroshima 739, Japan b Juken Sangyo Co., Ltd., 1-1 Mokuzaikou-Minami, Hatsukaichi-shi, Hiroshima 738, Japan Received 17 March 1997; accepted 8 July 1997

Abstract In this paper, we focus on a real size manpower allocation problem. It was modeled after a real world problem of distributing the salesmen force over the branches of a company. The problem includes multiple objectives and the number of salesmen at each branch is unspeci®ed. Conventional integer programming approach and conventional metaheuristics seem to have problems with solving the large size version of this problem. The versatility of our proposed heuristics based on a modi®cation of genetic annealing is exempli®ed through solving the real size manpower allocation problem. For comparison sake, several small sized versions were solved using our method, conventional integer programming approach, and some well known metaheuristics. Ó 1998 Elsevier Science B.V. All rights reserved. Keywords: Assignment; Genetic annealing; Genetic algorithm; Simulated annealing; Fuzzy programming

1. Introduction One of the managerial problems faced by most companies is manpower allocation. Depending on the company and the purposes of the allocation, the problem de®nition and, in particular, the constraints would vary. However, it all boils down to assigning a number of workers to a number of jobs. An instance of this, known as the assignment problem, is a special case of the transportation

*

Corresponding author. Fax: +81 824 24 7694. Fax: +81-824-22-7195; e-mail: [email protected] 2 Fax: +81-829-32-6237. 1

0377-2217/98/$19.00 Ó 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 7 7 - 2 2 1 7 ( 9 7 ) 0 0 3 2 6 - 3

problem. The problem is to assign n persons to n jobs. This problem is a very simple model of manpower allocation, so that it can be solved by Kuhn's Hungarian Algorithm (Bazaraa et al., 1990). However, in a real world manpower allocation problem, several factors, such as individual ability, individual preference and satisfaction, etc., are usually simultaneously present. Moreover, the required number of workers for each job is not always speci®ed. Generally speaking, such a real world problem is not easy to solve. To tackle such a problem, one of two avenues may be chosen. The ®rst is to develop a mathematical analysis to obtain the optimal solution to an even more complex manpower allocation problem. Konno and Ishii's research (Konno and Ishii, 1995a, b) can be

406

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

included in this category. The second is a heuristic approach for obtaining an approximately optimal solution. In this paper, we opt for the latter. The particular kind of the manpower allocation problem that we are treating in this paper was modeled after a real world problem. The problem is to distribute the salesmen force over the branches of a company satisfying several goals simultaneously and considering the salesmen's abilities, satisfactions and preferences. In this problem, the required number of salesmen for each branch is not speci®ed. Instead, for each branch, a fuzzy target gross sales is given. A total target gross sales of all branches is also fuzzily speci®ed. To tackle this problem, ®rst, representing all fuzzy target gross sales as fuzzy sets (fuzzy goals), we formulate the problem as a fuzzy combinatorial programming problem. Applying the Bellman and Zadeh's fuzzy decision method (Bellman and Zadeh, 1970), the problem is transformed into a mixed integer programming problem (MIP). Such an MIP includes a lot of 0±1 decision variables and has a special structure, called the Special Ordered Sets (SOS) structure (Ogryczak, 1996). This results in a problem that is not easily solved. The branch and bound method, usually used in solving MIP's, becomes practically ine€ective in solving such kind of problems when considering the time it takes to reach optimality. Judging these facts to be a serious performance hindrance, we developed a heuristic approach delivering a nearly optimal solution in a practical amount of time. This was in accordance with what takes place in reality, since, most of the time, an exact optimal solution is not what the management is looking for but a guiding solution for their decisions. Our heuristic approach was based on the genetic annealing (GAn) metaheuristic (Price, 1994). It was adopted due to its versatility as an extension of genetic algorithm (GA) and simulated annealing (SA). However, we did some modi®cations to GAn allowing us to obtain an earlier, better solution. Our paper is organized as follows. Section 2 would deal with the problem statement and formulation. Section 3 would introduce the used metaheuristic, namely, genetic annealing. Section 4 would discuss the coding method, the parameters settings, and the constraints handling. Section 5

would compare the approaches taken to solve the small size problems. Section 6 would explain the modi®cation done to the genetic annealing metaheuristic and to genetic algorithm metaheuristic to improve their respective performance when applied to the real size problem. In Section 7, we would present some concluding remarks. 2. Problem statement and formulation In our problem, the aims of the headquarters are to: · maximize the total gross sales goal of all branches, · maximize the gross sales of each branch, · maximize the satisfaction of each salesman. A salesman satisfaction is set according to his preference list specifying his degree of satisfaction when assigned to any of the branches. Furthermore, an ability coecient is set for each worker. Thus, the gross sales for every salesman when assigned to a branch was calculated using the following equation: gross sales of salesmani in branchj ˆ satisfaction factor  total gross sales aspiration level  ability coefficient = total number of salesmen: On the other hand, the two constraints imposed by the headquarters are that: · each salesman should be assigned to only one branch, · at least two salesmen should be assigned per branch. Our problem can be formulated as a multiobjective 0±1 integer problem as follows: n X m X gij xij …1† maximize z1 …x† ˆ iˆ1 jˆ1

maximize

n X zj2 …x† ˆ gij xij ;

j ˆ 1; . . . ; m

…2†

i ˆ 1; . . . ; n

…3†

iˆ1

maximize subject to

zi3 …x† ˆ

m X sij xij ; jˆ1

m X xij ˆ 1; jˆ1

i ˆ 1; . . . ; n;

…4†

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420 n X

xij P 2;

j ˆ 1; . . . ; m;

…5†

iˆ1

xij 2 f0; 1g;

i ˆ 1; . . . ; nand j ˆ 1; . . . ; m;

…6†

where n is the total number of salesmen, m the total number of branches, gij the gross sales of the ith salesman when assigned to the jth branch, sij the satisfaction of the ith salesman when assigned to the jth branch, and xij the binary decision variable whose values are set according to the following:

xij ˆ

8 0 > > > > <

if ith salesman is not assigned to

> 1 > > > :

if ith salesman is assigned to

j th branch

the constraints (4)±(6) can be reduced to the following single objective MIP: maximize k subject to lGtg lGgj lGsi

…7†

j th branch:

As for the ®rst objective function, z1 …x†, it is to maximize the total gross sales of all branches. The second group of objective functions, the zj2 …x†'s, are to maximize the gross sales of each branch, and the third group of objective functions, the zi3 …x†'s, are to maximize the satisfaction of each worker. To deal with the multiple objectives in our problem (1)±(6), we introduce a goal to each objective function. However, the target values cannot be speci®ed precisely. Consequently, we represent them as the fuzzy goals Gtg , Ggj where j ˆ 1; . . . ; m, and Gsi where i ˆ 1; . . . ; n, and apply the fuzzy programming approach (Sakawa, 1993; Zimmermann, 1996). Based on Bellman and Zadeh's fuzzy decision method (Bellman and Zadeh, 1970), the objective functions can be aggregated using the minimum operator as  lD …x† ˆ min lGtg …z1 …x††; min lGgj …z2j …x††; jˆ1;...;m  …8† min lGsi …z3i …x†† ; iˆ1;...;n

where lGtg , lGgj and lGsi are the membership functions of Gtg , Ggj , and Gsi , respectively. Our problem becomes that of maximizing the aggregated objective function (8) under the constraints (4)±(6). Introducing an auxiliary variable k, the problem of maximizing the objective function (8) under

407

n X m X

…9† !

gij xij

iˆ1 jˆ1 n X

P k;

…10†

!

gij xij

P k;

j ˆ 1; . . . ; m;

…11†

i ˆ 1; . . . ; n;

…12†

iˆ1

! m X sij xij P k; jˆ1

m X xij ˆ 1;

i ˆ 1; . . . ; n;

…13†

j ˆ 1; . . . ; m;

…14†

jˆ1 n X xij P 2; iˆ1

xij 2 f0; 1g;

i ˆ 1; . . . ; n and j ˆ 1; . . . ; m;

0 6 k 6 1;

…15† …16†

where k is the only continuous variable. In this problem, k shows the membership degree of the aggregated objective function, lD , implying the addition of constraint (16). However, if the membership functions of the fuzzy goals are tightly de®ned, this problem may have no feasible solutions. On the other hand, if the membership functions are too relaxed, this problem may have a lot of nearly optimal solutions. Therefore, we consider the problem without constraint (16), so that we can obtain the nearest solution possible to the reservation level in the tight membership function case, and the farthest solution possible from the aspiration level in the relaxed membership function case. By doing this, we can supply relevant information to the decision maker in which direction to modify his membership functions. Due to their tractability and their ease of use, linear fuzzy membership functions were used to represent the fuzzy goals. They could be written as:

l…r† ˆ

8 0 > < 0 rÿy

y 1 ÿy 0 > : 1

if r < y 0 ; if y 0 6 r 6 y 1 ; 1

if r > y :

…17†

408

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

Eq. (17) corresponds to the linguistic expression ``roughly larger than or equal to y 1 '' (see Fig. 1). From this equation, we can see that a linear membership function is de®ned by two parameters y 1 , the aspiration level, and y 0 , the reservation level, which we must specify. Those parameters are determined by the corresponding goal given by the headquarters. Thus, Eqs. (10)±(12) would have a similar membership function to Eq. (17). Adopting such membership functions and dropping constraint (16), problem (9)±(16) becomes: maximize k subject to P P … i j gij xij † ÿ gt0 …

gt1 ÿ gt0

P i

gij xij † ÿ gj0 gj1 ÿ gj0

…

Pm

m X jˆ1 n X

jˆ1 sij xij † s1i ÿ s0i

…18†

P k;

P k;

ÿ s0i

P k;

…19†

j ˆ 1; . . . ; m; i ˆ 1; . . . ; n;

…20†

xij ˆ 0 or xij ˆ 1:

…21†

Applied on an SOS variable leads to the dichotomy

xij ˆ 1;

i ˆ 1; . . . ; n;

…22†

xij P 2;

j ˆ 1; . . . ; m;

…23†

iˆ1

xij 2 f0; 1g;

i ˆ 1; . . . ; n and j ˆ 1; . . . ; m:

…24†

When compared to the problems in Bazaraa et al. (1990) Konno and Ishii (1995a, b), (18)±(24) is quite complex and presents more diculties to be

Fig. 1. Roughly larger than or equal to y 1 .

dealt with. To get the exact optimal solution of such a problem using the branch and bound method would not be practical. This is due to the fact that a problem with a couple of hundred salesmen and a couple of dozens branches would result in a couple of thousands 0±1 variables that would require a solution time beyond any practicality. A second diculty comes from the n constraints (22) which are creating n Special Ordered Sets (SOS) (Ogryczak, 1996). SOS: In fact, the n constraints (22) are nothing but an algebraic representation of our n logical multiple choice requirements. Put in other words, for an i, a salesman, only one xij is set to 1, corresponding to only one j of the m branches. SOS creates the following problem with the standard branch and bound algorithm for the MIP. The standard branching rule on 0±1 variables is …25†

xi1 ‡


‡ xi…jÿ1† ‡ xi…j‡1† ‡


‡ xim ˆ 1 or xij ˆ 1: This creates an extremely unbalanced branching on the set of alternatives. It causes a low e€ectiveness of the branch and bound algorithm which may lead to an enormously long computation process, making impossible to solve even quite small problems in a reasonable amount of time. In Ogryczak (1996), a modeling technique, called Special Ordered Inequalities (SOI), was introduced to solve the SOS problem. The way SOI works is by introducing a new 0±1 variable, yij , de®ned as: for i ˆ 1; . . . ; n; ( j ˆ 1; xi1 ; …26† yij ˆ xi…jÿ1† ‡ xij ; j ˆ 2; . . . ; m: Replacing xij by yij would transform the SOS problem into an SOI problem where branching on the yij 's follows the standard branching rule (25). However, since an optimal solution is not sought and is practically too much to ask for due to the size of the problem, metaheuristics, such as genetic annealing, can be a very practical alternative to provide a near optimal solution.

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

409

3. Genetic annealing In Price (1994), an approach called genetic annealing (GAn), combining both the GA (Genetic Algorithm) and the SA (Simulated Annealing) techniques, was claimed to take the best from both genetic and annealing worlds. Such a claim made us opt for this metaheuristic. Later on, we would see that it was a favorable choice. A closer look at GAn reveals that it is an extension of the Creutz ``demon'' algorithm (Creutz, 1983) using threshold acceptance similar to Dueck and Scheuer (1990). However, the use of thermodynamic annealing criteria is not unique to GAn (Lin et al., 1993). Like GA (Goldberg, 1989; Michalewicz, 1996), GAn is a population based metaheuristic using operators to cross-breed or to mutate the population members. Three operations are used in GAn: swapping, re¯ection, and splicing. Every time we try to mutate a string, one of the three operations is chosen at random according to their respective probabilities. Swapping: When using the ®rst operation, we can swap a pair of genes, or we can explore remote regions of solution space more e€ectively if we occasionally swap more than one pair of genes at a time. To decide on the number of taken swaps, ®rst, we choose two random and di€erent genes, make a swap, then generate a random number, 0 6 rx 6 1, and compare rx to a preset constant number EXSW 2 ‰0; 1†. Only if rx 6 EXSW, another swap is made. An example of a mutation consisting of two swaps is illustrated in Fig. 2. Re¯ection: The second operation can generate large scale variations in a chromosome by re¯ecting two groups of genes. To determine the size of the re¯ected groups of genes, ®rst, we choose a random re¯ection axis, then re¯ect the ®rst genes on either sides of the axis. After that, a random number, 0 6 rx 6 1, is generated and compared to a preset constant number EXR 2 ‰0; 1†. Only if rx 6 EXR, the second genes on either sides of the axis are re¯ected; and so on. An example of a mutation consisting of a three genes re¯ection is illustrated in Fig. 2. Splicing: The third operation parallels crossover in GA. It proceeds as follows. First, we choose a

Fig. 2. GAn operators.

random donor chromosome, then choose a random starting point and a random length for the splice. After that, the splice taken from the donor replaces the corresponding portion in the target chromosome. However, the donor chromosome remains unaltered. An example of a splice of length four is illustrated in Fig. 2. In GAn, we assign an energy threshold, Thi , to each chromosome. Initially, each threshold equals the energy (the ®tness function) of the chromosome to which it is assigned. In our problem, the energy of a chromosome is calculated using (8). Acceptance criteria: A chromosome threshold determines which trial mutations constitute acceptable improvements. If the energy of a mutant is less than the threshold of the chromosome that created it, we reject the mutant and move on to the next chromosome. However, if its energy is greater than or equal to the threshold, we accept the mutant as a replacement for its creator. GAn uses

410

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

an energy bank, DE, to keep track of the energy liberated by successful mutants. Whenever a mutant passes the threshold test, we add the di€erence between the threshold and the mutant energy to DE for temporary storage. Then, we make the threshold equal to the accepted mutant energy and move on to the next chromosome. Cooling: After each chromosome has been subjected to a random mutation, we cool down the population by lowering each threshold by a quantity de ˆ DE  C=N , where C is the cooling constant taken between 0 and 1, and N is the population size. Annealing results from repeated cycles of spontaneous heating of successful mutants and then a uniform cooling by lowering the threshold energy of each population member equally. After they have been cooled down, thresholds are lower than the energies of the chromosomes to which they have been assigned. This means that, sometimes, we are forced to accept a mutant with a lower energy than the energy of the chromosome it is replacing. In essence, the entire population acts like a giant heat reservoir that exchanges energy among its members. Thus, contrary to GA, GAn uses some measure, based on the comparison of the mutant energy (®tness) with its parents threshold, to decide whether to accept or to reject the mutant. This way, a more strict survival of the ®ttest, relative to that of GA, is imposed. This measure is similar to SA's. However, traditional SA (Aarts and Korst, 1989) is based on the Metropolis algorithm, which is computationally expensive, since the mutant acceptance criterion is a comparison between a randomly generated number and an exponential term. Con®gurations that are not much worse than their creators may be rejected depending on the generated random number. Whereas, genetic annealing accepts every con®guration that is not much worse than its creator based on the result of the energy-threshold comparison that requires neither a random number generator nor an exponential term. Traditional SA also relies on an empirically derived ``annealing schedule'' to control the rate at which the Metropolis temperature should be lowered. The advantage of using thresholds to track the time-averaged gain of spontaneous energy from individual con®gura-

tions is that we can maintain equilibrium at any temperature just by depleting energy at each threshold at the same average rate that the ensemble gains it. 3.1. The algorithm In the following, we present the algorithm of GAn. PS, PR, and PX are, respectively, the probabilities that a mutation will swap random genes, re¯ect a group of genes, or splice a group of genes. Naturally, PS ‡ PR ‡ PX ˆ 1. Step 1. Randomly create an initial population of N con®gurations. Step 2. For (i ˆ 1 to N) initialize the ith threshold, Th‰iŠ, with the energy of the ith con®guration. Step 3. DE ˆ 0 /*empty the energy bank*/. Step 4. /*begin heating loop*/ For (i ˆ 1 to N ) Splice, re¯ect, or swap the ith con®guration Compute the energy, E, of the resulting mutant If E < Th‰iŠ then restore the old con®guration If E P Th‰iŠ then: a. DE ˆ DE ‡ Th‰iŠ ÿ E /*add energy di€erence to DE*/ b. Th‰iŠ ˆ E /*reset threshold*/ c. Replace old con®guration with successful mutant Mutate next con®guration or end heating loop. Step 5. de ˆ DE  C=N /*compute cooling decrement*/. Step 6. /*begin cooling loop*/ For (i ˆ 1 to N ) Th‰iŠ ˆ Th‰iŠ ‡ de /*add de to each threshold*/ Return to Step 3 once all thresholds have been cooled down or end if the maximum number of iterations is reached. 3.2. Numerical example To make everything very clear, let us consider the following simple application of the algorithm as illustrated in Fig. 3. Step 1: Create an initial population consisting of three chromosomes, ch1, ch2, and ch3. Let us

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

411

Fig. 3. GAn example.

assume a very simple ®tness function whose value is the sum of the genes. Thus, for example, the energy(®tness) of ch1 is E ˆ 4 ‡ 1 ‡ 2 ‡ 3 ˆ 10. Step 2: Initially, the thresholds are equal to the energies: Th‰1Š ˆ 10, Th‰2Š ˆ 9, and Th‰3Š ˆ 9. Step 3: DE ˆ 0. For the sake of this example, let us apply only the splicing operator. Step 4: ch1 gets a splice of ch2 (the lighter shaded second and third genes). E ˆ 4 ‡ 2 ‡ 2 ‡3 ˆ 11 which is greater than Th‰1Š. Then, DE ˆ DE ‡ Th‰1Š ÿ E ˆ 0 ‡ 10 ÿ 11 ˆ ÿ1, Th‰1Š ˆ E ˆ 11, and ch1 is replaced by its mutant (appearing in Population 2 of Fig. 3). ch2 gets a splice of ch3 (the darker shaded third, fourth and ®rst genes). E ˆ 3 ‡ 2 ‡ 1 ‡ 1 ˆ 7 which is less than Th‰2Š. Then, restore old con®guration of ch2 (as appearing in Population 2 of Fig. 3). ch3 gets a splice of ch1 (the underlined third and fourth genes). E ˆ 3 ‡ 4 ‡ 2 ‡ 3 ˆ 12 which is greater than Th‰3Š. Then, DE ˆ DE ‡ Th‰3Š ÿ E ˆ ÿ1 ‡ 9 ÿ 12 ˆ ÿ4, Th‰3Š ˆ E ˆ 12, and ch3 is replaced by its mutant (as appearing in Population 2 of Fig. 3). Step 5: Assuming that C ˆ 0.96, de ˆ DE C=N ˆ ÿ4  0:96=3 ˆ ÿ1:28.

Step 6: Th‰1Š ˆ Th‰1Š ‡ de ˆ 11 ÿ 1:28 ˆ 9:72, Th‰2Š ˆ Th‰2Š ‡ de ˆ 9 ÿ 1:28 ˆ 7:72, and Th‰3Š ˆ Th‰3Š ‡ de ˆ 12 ÿ 1:28 ˆ 10:72. Then we will return to Step 3 and so on. . . 4. Application of the methods 4.1. Test problems settings We used the following membership functions settings for all our problems. Reservation and aspiration levels: The reservation and aspiration level for each branch was set as 0% and +10%, respectively, from its estimated gross sales (EGS) for big cities …EGS P 200†, and as )10% and +10% respectively from its EGS for other cities …EGS < 200†. The total gross sales reservation and aspiration levels are calculated as )5% and 10% from the sum of all EGS's. As we have said before, a salesman satisfaction is set according to his preference list specifying his degree of satisfaction when assigned to any of the branches. Hidden constraints: It was deemed practical and in accordance with what takes place in the real

412

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

world that the salesmen would not be asked to make a lengthy preference list of all branches. Since human beings usually make preferences based on a subgroup of the available alternatives, a rather shorter version of the preference list was constructed. An example is found in Table 1, where, for each salesman, the ®rst column contains the chosen branch index, and the second column contains the corresponding gross sales if assigned to that branch. As we can see from Table 2, p1 has ®ve branches (m ˆ 5) and a preference list length of 3 (r ˆ 3). Reducing the length of the list from 5 to 3 could seem a little bit trivial, but if we consider the problems of Table 3, it would become obvious that it is only natural to make a preference list of length 10 (r ˆ 10) out of 33 possible branches (m ˆ 33). However, this added a hidden constraint to the problem, speci®cally, only a subgroup of the branches is feasible for a salesman. In order to

check for an even worse case that is possible in the real world, we added one more hidden constraint where all salesmen would include all big cities branches in their preference list. For example, in Table 1, the only big city, 3, is included in all the preference lists. The preference lists were ordered according to the decreasing order of preference, with the highest preference being at the top of the list. A satisfaction coecient was set accordingly. The satisfaction factor for rank 0 in the list is 1.0, for rank 1 is 0.95, for rank 2 is 0.9, for ranks 3±f is 0:9 ÿ …ranknumber ÿ 2†  0:2, where f is the length of the preference list. After this, we make sure that every branch was chosen by at least two salesmen; otherwise, any such branch is appended to every salesman list with 0.7 satisfaction factor. Furthermore, an ability coecient was set for each salesman as follows: 1.5 to the ®rst 10% of the

Table 1 Preference lists for problem p1 Salesman-1

Salesman-7

Salesman-13

Salesman-19

Salesman-25

Salesman-31

3 35.2917 2 33.5271 1 31.7625

3 29.4097 2 27.9392 1 26.4688

3 23.5278 4 22.3514 1 21.175

3 23.5278 0 22.3514 2 21.175

3 23.5278 4 22.3514 2 21.175

3 17.6458 2 16.7635 0 15.8812

Salesman-2

Salesman-8

Salesman-14

Salesman-20

Salesman-26

Salesman-32

3 35.2917 1 33.5271 2 31.7625

3 29.4097 2 27.9392 4 26.4688

3 23.5278 4 22.3514 0 21.175

3 23.5278 4 22.3514 0 21.175

3 23.5278 0 22.3514 1 21.175

3 17.6458 4 16.7635 1 15.8812

Salesman-3

Salesman-9

Salesman-15

Salesman-21

Salesman-27

Salesman-33

3 35.2917 0 33.5271 4 31.7625

3 29.4097 2 27.9392 4 26.4688

3 23.5278 1 22.3514 4 21.175

3 23.5278 1 22.3514 0 21.175

3 23.5278 2 22.3514 0 21.175

3 17.6458 1 16.7635 2 15.8812

Salesman-4

Salesman-10

Salesman-16

Salesman-22

Salesman-28

Salesman-34

3 35.2917 2 33.5271 0 31.7625

3 29.4097 1 27.9392 4 26.4688

3 23.5278 1 22.3514 4 21.175

3 23.5278 1 22.3514 4 21.175

3 23.5278 2 22.3514 1 21.175

3 11.7639 1 11.1757 0 10.5875

Salesman-5

Salesman-11

Salesman-17

Salesman-23

Salesman-29

Salesman-35

3 29.4097 2 27.9392 0 26.4688

3 23.5278 1 22.3514 0 21.175

3 23.5278 4 22.3514 1 21.175

3 23.5278 2 22.3514 0 21.175

3 17.6458 1 16.7635 2 15.8812

3 11.7639 2 11.1757 1 10.5875

Salesman-6

Salesman-12

Salesman-18

Salesman-24

Salesman-30

Salesman-36

3 29.4097 0 27.9392 2 26.4688

3 23.5278 2 22.3514 1 21.175

3 23.5278 4 22.3514 2 21.175

3 23.5278 0 22.3514 4 21.175

3 17.6458 2 16.7635 1 15.8812

3 11.7639 0 11.1757 2 10.5875

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

413

Table 2 Small scale problems comparison MOMIP

SA

GA

GAn

p1 (n ˆ 36) (m ˆ 5) (p ˆ 3)

Problem LB UB ) Time

0.683455 0.890655 ) 7 h 28 min

Max Min Avg Var

0.749635 0.713543 0.7184901 1.27466 ´ 10ÿ4

0.808455 0.738606 0.8008153 4.8199 ´ 10ÿ4

0.808455 0.72022 0.7786768 1.506072 ´ 10ÿ3

p2 (n ˆ 26) (m ˆ 4) (p ˆ 3)

LB UB ) Time

0.735593 0.903415 ) 3 h 23 min

Max Min Avg Var

0.735593 0.731281 0.7350543 1.87 ´ 10ÿ6

0.735593 0.731276 0.7317121 1.86 ´ 10ÿ6

0.735593 0.731281 0.7351618 1.86 ´ 10ÿ6

p3 (n ˆ 20) (m ˆ 3) (p ˆ 3)

LB UB ) Time

0.808458 0.824844 ) 2 h 51 min

Max Min Avg Var

0.808458 0.808456 0.8084578 4.0 ´ 10ÿ13

0.808458 0.808458 0.808458 0

0.808458 0.793753 0.8069867 2.16 ´ 10ÿ5

p4 (n ˆ 19) (m ˆ 3) (p ˆ 3)

LB UB ) Time

0.614004 0.7285 ) 2 h 39 min

Max Min Avg Var

0.614033 0.614033 0.614033 4.93 ´ 10ÿ17

0.614033 0.614031 0.6140328 4.0 ´ 10ÿ13

0.614033 0.614033 0.614033 4.93 ´ 10ÿ17

p5 (n ˆ 18) (m ˆ 3) (p ˆ 3)

Opt. ) ) Time

0.833856 ) ) 2 h 06 min

Max Min Avg Var

0.833856 0.770833 0.8008165 8.17005 ´ 10ÿ4

0.833856 0.833856 0.833856 0

0.833856 0.770833 0.8212514 7.06115 ´ 10ÿ4

p5-b (n ˆ 18) (m ˆ 3) (p ˆ 3)

LB UB ) Time

0.777414 0.926533 ) 4 h 20 min

Max Min Avg Var

0.777414 0.777414 0.777414 0

0.777414 0.759153 0.7737618 5.93 ´ 10ÿ5

0.777414 0.777414 0.777414 0

p6 (n ˆ 17) (m ˆ 3) (p ˆ 3)

Opt. ) ) Time

0.878681 ) ) 1 h 40 min

Max Min Avg Var

0.878681 0.836217 0.8447095 3.20545 ´ 10ÿ4

0.878681 0.793753 0.8574478 9.01539 ´ 10ÿ4

0.878681 0.836217 0.8701864 3.20533 ´ 10ÿ4

p7 (n ˆ 13) (m ˆ 2) (p ˆ 2)

Opt. ) ) Time

0.879114 ) ) 15 s

Max Min Avg Var

0.879114 0.879114 0.879114 5.92 ´ 10ÿ16

0.879114 0.879114 0.879114 5.92 ´ 10ÿ16

0.879114 0.879114 0.879114 5.92 ´ 10ÿ16

p7-c (n ˆ 13) (m ˆ 2) (p ˆ 2)

Opt. ) ) Time

0.917148 ) ) 2s

Max Min Avg Var

0.917148 0.917148 0.917148 1.97 ´ 10ÿ16

0.917148 0.917148 0.917148 1.97 ´ 10ÿ16

0.917148 0.904809 0.9134463 3.55 ´ 10ÿ5

salesman, 1.25 to the next 15%, 1.0 to the next 50%, 0.75 to the next 15%, and 0.5 to the last 10%. Thus, the gross sales for every salesman when assigned to a branch was calculated using the following equation: gij ˆ

satisfaction factor  ability coefficient  gt1 : n …1†

4.2. Parameters of the methods To check the performance of GAn, we compared it to its two predecessors, GA and SA, on one hand, and to the solution given by momip, 3 an 3 Available for academic use from http://www.iiasa.ac.at/ marek/.

414

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

Table 3 Real size problem: Standard vs. modi®ed Problem

GAn

GAnm

GA

SA

GAn

GAna

p8 (n ˆ 215) (m ˆ 33) (r ˆ 10)

Max Min Avg Var

0.22544 0.14024 0.18453 0.00079

0.36030 0.18322 0.28635 0.00272

)0.47312 )0.64930 )0.54345 0.00308

)1.25828 )1.48449 )1.38441 0.00563

0.10516 0.01320 0.05193 0.00123

0.14458 0.01499 0.08584 0.00225

p8-r (n ˆ 215) (m ˆ 33) (r ˆ 10)

Max Min Avg Var

0.29526 0.21172 0.24995 0.00060

0.38570 0.27877 0.34092 0.001411

)0.19689 )0.32529 )0.27463 0.00203

)0.68434 )0.89478 )0.81521 0.00545

0.20210 0.0455 0.13961 0.00237

0.33790 0.18085 0.24850 0.00182

MIP software that is based on branch and bound and that takes care of the SOS structure. Let's ®rst start by considering the setting of the parameters of the above methods. Stopping criteria: The stopping criteria of SA, GA, and GAn was reaching the maximum number of iterations, while the stopping criteria of momip was ®nding an optimal solution, reaching a preset limit on the number of nodes, or reaching a preset limit on the number of visited nodes without improving the current solution. Coding and constraints handling: Along with operators de®nition, coding is one of the most important factor in¯uencing the performance of a metaheuristic. The ®rst kind of coding, that directly came to our mind, was to make a chromosome where each gene would represent a salesman; whereas, the value stored in each gene would be the index of the branch to which the salesman is assigned. This way, we will satisfy automatically constraint (22). However, the preference lists of the salesmen contain only a subgroup of all branches, the probability of creating an infeasible mutant after applying any of the operators of any of the metaheuristics is very high. This would require a mutant repairing procedure that would be de®nitely computationally expensive. A better coding would be to store in each gene the index of the preference list of each salesman. This way, we are still satisfying constraint (22), but we do not have the infeasibility problem of the previous coding. Nevertheless, we needed a procedure that would ensure the satisfaction of constraint (23). Such a procedure makes a list, vl , of all branches violating Eq. (23), then ®xes the violation of each of these branches one at a time as follows. As an

example, let us consider branchv from the vl list. We look for salesmen who have branchv on their preference list, and pick up the one among them such that, if his assignment is changed from his current branchc into branchv , branchc would not violate constraint (23). Initial population: It is randomly created by assigning random preference list indexes to the genes. However, the satisfaction of constraint (23) was ensured for each initial chromosome using the procedure described above. SA: A classical version of SA was used. The ®rst implementation, that came to our mind, was to use an operator that would perform a single swap of a pair of genes. However, such an implementation would limit the ability of improving on the initial solution, due to a simple fact that we are ®xing the number of each preference index included in the chromosome. For example, if we had three salesmen with a 0-index assignment in our initial solution, we will have only three salesmen (not necessarily the same ones) with a 0index assignment in all of our generated solutions. Thus, the alternative was to have two operators: the ®rst would still perform a single swap; whereas, the second would mutate one randomly chosen gene. Choosing between the two was done randomly with the ®rst having a probability of 0.6 and the second 0.4. As for the rest of the SA parameters, the cooling rate (a) was set to 0.96 and the temperature at each cycle was updated using Tnew ˆ a  Told . The temperature cycle length was determined using two parameters, the maximum number of swaps and of successful swaps per cycle. The transition probability was calculated as

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

 pt ˆ exp

 lmutant ÿ lparent D D : Tcurrent

GA: We used a simple GA with roulette selection and an elitist model, where the best two individuals survive unconditionally. The population size was set to 100. The probability of crossover and mutation was set to their respective classical standard, i.e., 0.7 and 1/n. We used a two points crossover. The ®tness function of an individual, fitind , was calculated using a scaling of Eq. (8) to render it positive by adding its population minimum if negative: if b ˆ

min …lind D † < 0; then fitind

population

ˆ fitind ÿ b: GAn: For GAn, after some experimentation with its parametric settings, the probability of splicing, re¯ecting, and swapping were set to 0.8, 0.1, and 0.1, respectively; whereas, EXSW and EXR were both set to 0.35. This is very much in correspondence with the standard GA crossover and mutation probability settings, where in GAn, splicing corresponds to crossover and re¯ection and swapping correspond to mutation. The population size was set to 50; whereas, the cooling constant, C, was set to 0.96, same setting as SA. We do not claim to tune SA and GA parameters to their best performance, but after some numerical experimentation, we opted to set them to their respective classical standards. That is why, we have let SA and GA run for twice the running time of GAn, and set the population size to 50 and 100 for GAn and GA, respectively. In the following, we will present some numerical simulations. All our programs were run on a 133 MHz Pentium PC using Borland C++ version 4.5 under Windows 95.

5. Comparisons of small size problems For the purpose of checking the proximity of GAn solution to the optimal solution, we used a set of small size problems for our simulations (Table 2). This scaling down of the real problem size is done due to the fact that the branch and

415

bound solution time increases exponentially with the size of the problem. Even for p1, which has only 180 0±1 variables, momip took 7 h and 28 min without reaching optimality. On the other hand, GAn took 2 min and SA and GA took 4 min for each random seed run. In the second column of Table 2, Opt. stands for optimal, LB for lower bound, and UB for upper bound values. LB and UB provided limit guidance when optimality was not reached. Some statistical results of runs with the same 10 random seeds for all three metaheuristics are shown in Table 2. For the problems in Table 2, the used metaheuristics took a couple of seconds to a couple of minutes; whereas, momip took hours for all but the smallest problems, p7 and p7-c. The number of solutions visited by SA, GA, and GAn were about 1.75, 1.5, and 1.25 million respectively. As can be seen from Table 2, except for p1, all metaheuristics gave comparable results for all problems. In the case of p1, as the number of 0±1 variables increased, SA started falling behind both GA and GAn. The di€erence in the solution quality would become more prominent as we will address the real world problem (Table 3). To have a clearer idea of the near optimality expectancy of the three methods at any of the random seeds, we plotted all 10 runs' results for all problems (Fig. 4). In Fig. 4, the broken line squares are the optimality area determined using momip's LB and UB, when momip did not reach optimality. However, in case momip reached optimality, a broken line ran horizontally to group all optimal answers of the di€erent methods. As can be seen from Fig. 4, all methods had at least one solution in the optimal area for all problems. That would allow us to draw the conclusion that SA, GA, and GAn can provide nearly optimal, if not optimal, solution in the case of small size problems. 6. Modi®cation of GAn and GA After the above successful simulations results, we felt more con®dent and moved on to consider problems with a more realistic size. To do so, we considered a real case with 215 salesmen and 33

416

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

Fig. 4. 10 runs comparison for problems of Table 2.

branches. We created two problems, p8 and p8-r. To test the e€ect of our second hidden constraint, i.e., all salesmen will include all big cities in their preference list, we created the preference lists of p8 accordingly; whereas, the preference lists of p8-r were created using a selection probability equal to the ratio of the branch estimated gross sales over the sum of estimated gross sales. The results of runs with 10 di€erent random seeds for GAn, GA, and their modi®cations, along with SA, are shown in Table 3. Each random seed running time for all methods was set to 37 min. This allowed SA, GA, GAn, GAnm, GAn , and GAna to visit about 1.2, 1, 0.85, 0.5, 1, and 0.5 million solutions, respectively. From the results in Table 3, we can see that the fact of including all big cities in the salesmen preference list has a

constraining e€ect on the problem. This shows in a lower maximum and average for all tested methods. Let us consider those modi®ed methods and the reasons why we made them. Our observation that GAn is spending a long time in local minima prompted us to consider modifying GAn into GAnm. This can be seen in Fig. 5 as longer plateaus for GAn compared to GAnm's. We deduced that the introduction of a diversi®cation technique to GAn is in order. Therefore, we modi®ed GAn into GAnm as follows. Since splicing is the operator that would spread uniformity in the population, we decided to couple it with a mutation. Thus, we check if, after splicing, the mutant is the same as the parent, we apply either re¯ection or swapping according to

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

417

Fig. 5. Solution quality vs. Time for the methods of Table 3.

their respective probabilities. This way, we try to escape local minima by pairing a mutation with the splicing. As we can see from Table 3, this modi®cation improved the solution of both p8 and p8-r, where the maximum, the minimum, and the average of GAnm are higher than GAn's. To further check our modi®cation, we applied GAnm to the problems in Table 2. GAnm proved not only to be robust by giving the same result for all the 10 runs of each problem, but also its results were all equal to the maximums reached by GAn. Since SA's and GA's performance in the case of small size problems were comparable to GAn's, we applied them to p8 and p8-r. However, all their solutions were infeasible. This shows as the nega-

tive k values in Table 3. This prompted us to consider why would GAn succeed where SA and GA had failed. In the case of SA, even the use of two operators, namely, swapping and mutation, didn't help much. SA lagged behind both GA and GAn. The cause can be analyzed as a lack of diversity creation power on the side of SA, i.e., one swap or one mutation performed on a single individual compared to a population undergoing crossovers, mutations, splicing, re¯ection, and multiple swapping in the case of GA and GAn. Put in other words, SA, being a local neighborhood search method, was not able to explore much of the huge solution space (more than 7000 0±1 variables as compared to less or equal to 180 for

418

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

the small size problems) in the given limited time. The fact that we were starting from an initial random solution resulted in a start in the infeasible area of our problem solution space. Thus, if the used metaheuristic, like SA, didn't have enough exploration power, it was not able to escape the infeasible area. However, the puzzle was the failure of GA and the success of GAn, eventhough both of them are population based metaheuristics. The main two di€erences between the two methods can be summarized as follows: The operators: In GA, we have two operators: crossover and mutation. Crossover would modify both parents into two new children; whereas, mutation would introduce some new foreign random genes to a child chromosome, i.e., not belonging to the genes that already constituted the chromosome in question. While in GAn, we have three operators: splicing, re¯ection, and swapping. Splicing would modify only one of the parents, the target, leaving the second parent, the donor, untouched. Both re¯ection and swapping would not introduce any foreign genes to a chromosome, but rather change the order of the already existing genes in the chromosome. In GA, both operations are performed on every mutant. While in GAn, only one of the three operations is performed on each mutant. The selection and acceptance method: In GA, individuals are allowed to produce o€springs according to a probability based on their ®tness function. This is known as roulette selection. Thus, not all individuals have the chance to reproduce. Moreover, all mutants are accepted unconditionally. To safeguard the best individuals, the elitist model insure the unconditional survival of the best two individuals. While, in GAn, every individual is allowed to reproduce, children are not accepted unconditionally. On the contrary, only if a child passes the threshold test, making him not a lot worse than his parent, he is allowed to replace his parent. Not adopting the elitist model nor the roulette selection, allowed all individuals to strive to be the ®ttest. Therefore, the survival of the ®ttest is not anymore restricted to the top two individuals but open to everybody. This would give more possibilities and more power to a population to improve.

To check the e€ect of each of those two main di€erences, we created GAn , which is a GA that uses GAn's operators, and GAna , which is a GAn that uses GA's operators. GAn would show us the e€ect of GAn's operators; whereas, GAna would show us the e€ect of GAn's selection and acceptance method. As can be seen from Table 3, both GA's modi®cations, GAn and GAna , gave feasible solutions with GAna giving a better solution than GAn . This would let us infer that the operators themselves have quite an in¯uence on the solution quality, but the higher e€ect is due to the selection and acceptance method. To check further each method's solution quality, we plotted its best solution as a function of time. As shown in Fig. 5, all methods gave a similar behavior to a GA with an elite model. From both Table 3 and Fig. 5, we can see that GAnm gave the best results, followed by GAna . We can attribute their common success to both their common selection and acceptance method and their common pairing of the breeding operation with a mutation operation. On the other hand, the lead of GAnm over GAna can be solely attributed to its operators. Therefore, to have an idea of how our best method's (GAnm's) computation time would change with the size of the problem, we plotted in Fig. 6 the time taken by 100 generations against the problem size. From Fig. 6, it is obvious that the computation time, as a function of the problem size, increases in a linear fashion. From this data, we can extrapolate that an increase of 1000 in the problem size would incur an increase of 3 s in the computation time of 100 generations. 7. Conclusion In this paper, we have examined the way of solving a real world manpower allocation problem that has more diculties than the problems solved in the previous literature. In the process, we have used several metaheuristics and created some modi®ed versions of most of them. By the numerical simulation of several small size problems, we have checked how much the metaheuristics

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

419

Fig. 6. GAnm computation time vs. Problem size.

solutions are close to the optimum. After our data analysis, we found out that our modi®ed version of the adopted genetic annealing metaheuristic gave the best results. We analyzed the reasons behind the success of our method where the simple genetic algorithm failed. Our work put once more under the light the importance of the encoding, the chosen operators, and the selection and acceptance method in a metaheuristic. References Aarts, E., Korst, J., 1989. Simulated Annealing and Boltzmann Machine. Wiley, New York. Bazaraa, M., Jarvis, J., Sherali, H., 1990. Linear Programming and Network Flows. Wiley, New York.

Bellman, R.E., Zadeh, L.A., 1970. Decision making in a fuzzy environment. Management Science 17, 141±164. Creutz, M., 1983. Micro canonical Monte Carlo simulation. Physics Review Letters 50 (19), 1411±1414. Dueck, G., Scheuer, T., 1990. Threshold accepting: A general purpose optimization algorithm appearing superior to simulated annealing. Journal of Computational Physics 90, 161±175. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison±Wesley, Reading, MA. Konno, T., Ishii, H., 1995a. Fuzzy distribution problem of workers (in Japanese). Journal of Japan Society for Fuzzy Theory and Systems 7 (3), 162±167. Konno, T., Ishii, H., 1995b. Fuzzy distribution problem with two types of workers (in Japanese). Journal of Japan Society for Fuzzy Theory and Systems 7 (3), 168±174. Lin, F., Kao, C., Hsu, C., 1993. Applying the genetic approach to simulated annealing in solving some NP-hard problems.

420

N. Abboud et al. / European Journal of Operational Research 111 (1998) 405±420

IEEE Transactions on Systems, Man, and Cybernetics 23 (6), 1752±1767. Michalewicz, Z., 1996. Genetic Algorithms + Data Structures ˆ Evolution Programs. Springer, Berlin. Ogryczak, W., 1996. A note on modeling multiple choice requirements for simple mixed integer programming solvers. Computers & Operations Research 23 (2), 199±205.

Price, K., 1994. Genetic annealing. Dr. Dobb's Journal 19 (10), 127±132. Sakawa, M., 1993. Fuzzy Sets and Interactive Multiobjective Optimization. Plenum Press, New York. Zimmermann, H.-J., 1996. Fuzzy Set Theory and Its Applications, 3rd edition. Kluwer Academic Publishers, Boston, MA.