GENETIC ALGORITHMS FOR SUPPLY MANAGEMENT PROBLEM WITH LOWER-BOUNDED DEMANDS

GENETIC ALGORITHMS FOR SUPPLY MANAGEMENT PROBLEM WITH LOWER-BOUNDED DEMANDS Pavel Borisovsky ∗,1 Alexandre Dolgui ∗∗,2 Anton Eremeev ∗∗∗,1 ∗

Omsk State Technical University, 11 pr. Mira, 644050 Omsk, Russia. E-mail: [email protected] ∗∗ Ecole des Mines de Saint Etienne, 158, Cours Fauriel 42023 Saint Etienne cedex 2, France E-mail: [email protected] ∗∗∗ Omsk Branch of Sobolev Institute of Mathematics 13, Pevtsov str. 644099, Omsk, Russia E-mail: [email protected]

Abstract: Two variants of genetic algorithm (GA) for solving the Supply Management Problem with Lower-Bounded Demands (SMPLD) are proposed and experimentally tested. The SMPLD problem consists in planning the shipments from a set of providers to a set of consumers minimizing the total transportation cost, given lower and upper bounds on shipment sizes, lower-bounded consumption and linear delivery costs. The first variant of GA uses the standard binary representation of solutions, the second one is based on the permutation representation and a greedy decoder. The experiments show the behaviour of the algorithms on different problem classes and indicate their competitiveness to the CPLEX solver. c 2006 IFAC Copyright ° Keywords: Genetic algorithms, Mathematical programming, Gradient methods

1. INTRODUCTION This paper presents two variants of genetic algorithm for solving the Supply Management Problem with Lower-Bounded Demands (SMPLD) formulated in (Chauhan et al., 2004). In this problem a set of providers supply some product to a set of consumers, the quantity of product that can be delivered lies between the given minimum and maximum values, and the costs proposed by each provider are linear functions of the quantity being

delivered if the shipment is opened. SMPLD consists in minimizing the total transportation cost: n X m X

zij (aij + cij xij ) → min

(1)

i=1 j=1 n X

xij ≥ Aj , j = 1, ..., m,

(2)

xij ≤ Mi , i = 1, ..., n,

(3)

i=1 m X j=1

1

Partially supported by INTAS grant 03-51-5501 and Russian Foundation for Human Sciences grant 04-0200238A. 2 Partially supported by INTAS grant 03-51-5501.

zij ∈ {0, 1}, i = 1, ..., n; j = 1, ..., m,

(4)

mij zij ≤ xij ≤ Mi zij , i = 1, ..., n; j = 1, ..., m.(5)

Here n is the number of suppliers; m is the number of consumers; variables zij indicate the presence of delivery from a provider i to a consumer j, while xij represent the shipment size. Aj is the minimal amount of product required for the consumer j; mij is the minimum quantity the provider i is prepared to deliver to the consumer j; Mi is the maximum quantity the provider i is able to deliver to any consumer and to all of them in total. Parameters aij , bij , mij , Mi , Aj are supposed to be integer and non-negative. In case all mij = 0, this problem becomes the well-known Fixed Charge Transportation Problem (FCTP) - see e.g. (Barr et al., 1981; Eckert and Gottlieb, 2002). The SMPLD is a modification of the concave cost supply management problem (SMP) considered in (Chauhan and Proth, 2003; Chauhan et al., 2005). Unlike SMPLD, instead of inequality (2), in the SMP there is an equality n X

xij = Aj , j = 1, ..., m,

i=1

and the cost functions are nonnegative and concave. Finding any feasible solution to the SMP is N P -hard even for m = 1, yet there exists a pseudopolynomial-time exact algorithm for it (Chauhan et al., 2005). The SMPLD is also N P -hard (it contains the MINIMUM KNAPSACK problem as a special case) and finding a feasible solution satisfying (2)– (5) is N P -hard even in the case m = 2. However, in the case m = 1 this problem admits a fully polynomial time approximation scheme (FPTAS) and there is a fast greedy 2-approximation algorithm for it (Chauhan et al., 2004). Two variants of genetic algorithm are proposed: the first variant (GAmipx) uses the standard binary representation of solutions and a problemspecific MIP-crossover; the second one (GAgrd) is based on the permutation representation and a greedy decoder. The decoder applies the greedy heuristic separately to a subproblem of each consumer, ordering the consumers according to genotype string. A similar genetic algorithm was developed for the SMP in (Borisovsky, 2004).

operator Cross produces the offspring from two parent individuals by combining and exchanging their elements. The mutation operator Mut adds small random changes to an individual. Let vector Π(t) of N solutions denote the current population on step t in the GA (N is called the population size). The formal scheme of the GA with steady state replacement is as follows: Steady state GA 1. Generate the initial population Π(0) at random. 2. For t := 1 to tmax do: 2.1 Π(t) := Π(t−1) . 2.2 Selection: choose p1 , p2 from Π(t) . 2.4 Apply mutation: p01 = Mut(p1 ), p02 = Mut(p2 ). 2.3 Produce c1 and c2 applying the crossover: (c1 , c2 ) = Cross(p01 , p02 ). 2.5 Choose the two worst individuals q1 , q2 in Π(t) and replace them by c1 , c2 . The choice of each parent on step 2.2 is done by the s-tournament selection: take s individuals by random from Π(t) and select the best one. The encodings of solutions in the GA are usually called genotypes. One of the most essential issues in development of a GA is the choice of representation of solutions in genotypes. Two different representations are considered in this paper.

3. GA WITH BINARY REPRESENTATION A solution may be encoded in the genotype as a Boolean matrix (zij ). It is easy to see that when all values of (zij ) are fixed, the best possible feasible values of (xij ) can be found by solving an LP-problem. However, this problem may have no feasible solutions. In order to deal with such cases we relax the original constraints and add a penalty function to the criterion: n X m X

cij xij + R

i=1 j=1

m X j=1

n X

αj + R

n X

βi → min (6)

i=1

xij ≥ Aj − αj , j = 1, ..., m,

(7)

i=1

2. GENETIC ALGORITHM Genetic algorithm (GA) is a random search algorithm that models a process of evolving a population of individuals (see e.g. (Reeves, 1997)). Each individual corresponds to some solution of the problem (feasible or maybe infeasible) and characterized by its fitness, which reflects the goal function value and the satisfaction of problem constraints. New individuals are built by means of crossover and mutation operators. The crossover

m X

xij ≤ Mi + βi , i = 1, ..., n,

(8)

mij zij ≤ xij ≤ Mi zij , αj ≥ 0,

(9)

i = 1, ..., n, j = 1, ..., m.

(10)

j=1

Let f (z) denote the goal function value of this problem for genotype z. Given a sufficiently large penalty parameter R, the optimal penalty term Pm α in problem (6)–(9) will be equal to 0 iff j j=1

the genotype z can be complemented by a set of values {xij } feasible in problem (1)–(5). If z maps into a solution without a penalty, then its fitness is assumed to be  −1 n X m X φ(z) =  zij (aij + cij xij ) . i=1 j=1

Otherwise assume φ(z) = −f (z). An important property of the binary representation is that for any SMPLD instance there is a genotype z∗ which encodes the optimal solution.

3.1 Standard Crossover and Mutation The first choice for the genetic operators to use with binary representation is usually the onepoint or uniform crossover and standard or onebit-flip mutation (Reeves, 1997). For instance, a modification of the one-point crossover for SMPLD can work as follows. Given two genotypes y1 , y2 , it produces genotypes z1 , z2 , where ½ 2 ½ 1 yij , if j ≤ χ; yij , if j ≤ χ; 2 z = zij = 2 yij , otherwise. yij , otherwise, ij for all i = 1, ..., n, j = 1, ..., m. Here χ is chosen randomly and uniformly between 1 and m − 1. The standard mutation of a genotype y produces random genotype z = Mut(y) so that P {zij = 1 − yij } = pmut , P {zij = yij } = 1 − pmut . for all i = 1, ..., n, j = 1, ..., m. Here parameter pmut defines the mutation intensity. Our preliminary experience with the SMP problem (Borisovsky, 2004) shows that the GAs based on the crossover and mutation of this type are rather inefficient. This is due to relatively high computational cost of each fitness evaluation, which requires to solve problem (6)–(9), and little use of problem specifics. Therefore we have developed a reproduction operator MIP-Crossover, combining both mutation and crossover, specially designed for the SMPLD. 3.2 MIP-Crossover Operator The first goal of the MIP-crossover operator is to find the best possible combination of two given parent genotypes y1 and y2 , if it is possible without extensive computations. This goal is similar to those of the optimized crossover (Aggarwal et al., 1997; Balas and Niehaus, 1998) and the LPcrossover (Eremeev, 1999). The second goal is to add some random changes to the parent solutions and to produce an offspring

different from the parents. This way we aim to introduce the effects of traditional mutation operators, adding diversification to the population. In order to achieve both goals, we formulate a supplementary mixed integer programming problem, obtained from (1)–(5) by deletion of all variables which have zero values in both parent genotypes, keeping a random subset Smut among them: 1 2 zij = 0 for all (i, j) 6∈ Smut : yij = yij = 0. (11)

Here Smut is randomly constructed so that any pair (i, j) is independently included into Smut with probability pmut . Besides that, 10 genotypes, including both parents y1 and y2 and the best individuals of Π(t) are cut off by the Benders cuts, using the approach from (Kolokolov and Levanova, 1996). Let vj and ui be the optimal dual variables linked to constraints (2) and (3) and let sij and tij be the optimal dual variables for the inequalities ”≤” in (5). If f is the cost of the so far best found solution, then the Benders inequality cutting off vector z is defined as follows: n X m X

yij (aij + mij sij − Mi tij )

(12)

i=1 j=1

≤f −1+

m X i=1

Mi ui −

n X

Aj vj .

j=1

If an optimal solution to this problem z∗mipx is found, then it is the output of MIP-crossover. In our implementation the CPLEX 9.0 solver is used to find the optimum of problem (1)–(5), (11), (12). To avoid time-consuming computations we do not start this solver if the number of the variables ν after the deletion exceeds a certain threshold µ (we use µ = 80). If ν > µ, then assume Smut = ∅ and try again. If ν is still greater than µ, then the Boolean variables zij = 1 are fixed for some (i, j) 1 2 such that yij = yij = 1, until ν reaches µ. If ν > µ after such fixing or if problem (1)–(5), (11), (12) is unsolvable then the MIP-crossover just outputs a randomly mutated genotype Mut(y1 ). The GA based on the Boolean representation and MIP-crossover is called GAmipx. The genotypes of initial population in this GA are randomly generated by the same method as in the GAgrd described in Section 4.1. Another option could be to choose the components of the Boolean genotypes independently, assigning them the value 1 with a fixed probability P1 . In the experiments this method turned out to be inefficient due to high probability of infeasible solutions and difficulties with tuning the parameter P1 . In order to ensure the population diversity, we check whether the genotype z∗mipx obtained in the MIP-crossover is already present in the current population Π(t) . If this is not the case, then z∗mipx

replaces the worst individual q in Π(t) . Otherwise, substitute the genotype q by Mut(y1 ). In GAmipx the mutation parameter pmut is adaptively chosen. Every time the MIP-crossover is performed successfully or a duplicate solution is encountered, parameter pmut is multiplied by 1.1. Whenever Mut(y1 ) is inserted into the population instead of a solution to problem (1)–(5), (11), (12) pmut is divided by 2. This adaptive mechanism is implemented to keep the right level of diversity in population so that the number of Boolean variables in the subproblems ν were close to threshold µ.

4. GENETIC ALGORITHM WITH NON-BINARY REPRESENTATION The non-binary representation is based on heuristic decoder which involves consecutive solving the problem separately for each consumer. One way to look at the reconstruction of solution from this encoding is to imagine a queue in which the consumers receive their deliveries, satisfying constraint (2) if it is possible. The genotype in this case is the permutation π = (j1 , ..., jm ) of consumers in the queue. Decoder(π) 1. Assume Mi0 := Mi , i = 1, ..., n. 2. For k := 1 to m do: 2.1 Solve the relaxed SMPLD for one consumer jk . 2.2 Reduce the capacities: Mi0 := Mi0 − xijk , i = 1, ..., n. 2.3 Fix the obtained values for variables xijk , i = 1, ..., n. 3. Find the values {zij } complementary to {xij }. On step 2.1 we solve a one-consumer problem of type (6)–(9) assuming the Boolean values zijk , i = 1, ..., n to be variable as well: let us denote for simplicity A = Ajk , ci = cijk , ai = aijk , mi = mijk , xi = xijk , zi = zijk . Then the relaxed oneconsumer SMPLD is formulated as follows. n X

(ai zi + ci xi ) + Rα → min,

(13)

i=1 n X

xi ≥ A − α,

(14)

i=1

mi zi ≤ xi ≤ Mi0 zi , i = 1, ..., n,

(15)

α ≥ 0, zi ∈ {0, 1}, xi ≥ 0, i = 1, ..., n.

(16)

Here R is supposed to be sufficiently large, so that the goal function value (13) of any feasible solution with α = 0 is less than the value of any

feasible solution with α > 0. This mixed-integer programming problem could be solved exactly but in order to save the time we apply the greedy 2approximation algorithm (Chauhan et al., 2004) to it. This greedy algorithm always returns a feasible solution to (13)–(16) and, it gives a solution with zero penalty term α whenever such solution exists. Here the greedy algorithm is followed by a refining procedure outlined below. 1. Greedy Algorithm Assume A’:=A. While A0 > 0 do: Let Si = min{max{A0 , mi }, Mi0 } ∀i. Choose i such that mi ≤ Mi0 and (ai + ci Si )/ min{A0 , Mi0 } is minimal. Set xi = Si , A0 = A0 − xi , Mi0 = Mi0 − xi . 2. Refining Procedure While A0 < 0 do: Choose i such that xi > mi and ci is maximal. Set xi = xi − min{xi − mi , −A0 }, A0 = A0 + min{xi − mi , −A0 }. If A0 > 0 when the Greedy heuristic terminates, then the solution has a non-zero penalty. It is easy to see that not every feasible solution to the SMPLD problem can be represented by a permutation this way. In fact, an example of solvable SMPLD instance can be constructed showing that even when the decoder solves all subproblems (13)–(16) to optimality, its output may be infeasible for any permutation π. To overcome, at least partially, this deficiency, a simple Local Improvement Heuristic has been proposed. Suppose four elements in matrix (xij ) are chosen at random (at least two of them must be non-zero), then the corresponding subproblem of size 2×2 is easily solved exactly. If a solution is improved, then matrix (xij ) is adjusted accordingly. In the Local Improvement Heuristic this random procedure is applied to each new individual 100 times, so this heuristic may be considered a local search method with random examination of the neighborhood.

4.1 Crossover, Mutation and Initial Population The non-binary representation described above is implemented in the genetic algorithm GAgrd. The initial population Π(0) consists of N randomly chosen permutations. This permutation encoding is similar to the encoding used in many GAs for the traveling salesman and scheduling problems. Thus the well-known partially-mapped crossover (see e.g. (Reeves, 1997)) and exchange mutation are employed. The exchange mutation is performed as follows: given a permutation of consumers π = (j1 , ..., jm ), choose two positions k and l by random end exchange jk and jl . In the

GAgrd the crossover is performed with probability pcross and the probability of mutation is pmut .

5. EXPERIMENTAL RESULTS The GA-Grd was programmed in Visual C++ 6.0 and the GAmipx was coded in OPL Script in ILOG OPL Studio 3.7. The experiments were carried out on Pentium-IV, 3GHz on the instances described below. • Series RAN consists of 8 SMPLD instances obtained from FCTP benchmarks with given n, m and mij = 0, produced by a version of NETGEN generator (Klingman et al., 1974) modified in (Barr et al., 1981) – see also (Gottlieb and Paulmann, 1998). • Series N370e consists of 15 larger FCTP instances with n = 100, m = 50, mij = 0, produced by the same generator in (Sun et al., 1998). Only the first 5 problems of this series are used here. • Series S consists of 8 random SMPLD instances with n = m = 50, Mi ∈ [400, 600], Aj ∈ [350, 550], aij ∈ [500, 4500], cij ∈ [20, 1020], all mij = 10 produced by a generator implemented in (Zaozerskaya, 2003). • Series MIS consists of 5 problems johnson8.4.4, rnd100.1, hamming6.4, c125.9 and mann.a27 obtained from DIMACS collection of MAXIMUM INDEPENDENT SET benchmarks by a reduction similar to the one in Theorem 2 (Chauhan et al., 2004). In these SMPLD instances n = |E| + 1 and m = |V |, where V and E are the sets of vertices and edges in the original DIMACS graph. Thus, n is ranging from 64 to 378 and m is from 507 to 1313. FCTP benchmarks are available at http://www. gamsworld.org/performance/plib/credits.htm, and the graphs from DIMACS set can be found by ftp://dimacs.rutgers.edu/pub/challenge/graph/. In order to put all of the algorithms into equal positions, each of them was given the same amount of computation time. For the smaller instances of series RAN and S the given time was 2 · 103 s. while for the rest of the problems it was 104 s. Within these time limits, the GAs were run independently for equal time slots 10 times on each instance. To compare the GAs to the advanced branch and bound techniques, we ran CPLEX 9.0 solver built into ILOG OPL Studio 3.7 with its default settings for 2 · 103 s. on series RAN and S and for 104 s. on all other problems. In experiments with the GAmipx s = 10, R = 1010 and µ = 80. The average number of iterations on series S and N370e was near 250 and 2600 respectively. The population size was set to 100

Table 1. Results for series RAN f∗ 3252 3664 3714 9711 5247 3823 4270 1373

ran13x13 ran12x21 ran14x18 ran4x64 ran8x32 ran16x16 ran10x26 ran17x17

fcplex 3252 3675 3736 9711 5247 3827 4270 1373

fGAmipx 3252 3664 3715 9711 5247 3823 4270 1373

fGAgrd 3284 3671 3714 9726 5247 3823 4270 1373

Table 2. Results for series S fbest 613951 631219 627294 634595 605145 669346 602737 645843

s1 s2 s3 s4 s5 s6 s7 s8

fcplex 614221 631219 629862 635162 605592 670535 602737 645843

fGAmipx 613951 631219 627294 634609 605471 670225 603114 645843

fGAgrd 614289 631540 627392 634872 605432 670730 603788 647242

Table 3. Results for series N370e e0 e1 e2 e3 e4

fbest 1208672 1204194 1192687 1188904 1209357

fcplex 1308469 1308955 1304851 1293843 1326480

fGAmipx 1289088 1282007 1273117 1276653 1285797

fGAgrd 1284070 1303920 1305210 1305520 1332310

Table 4. Results for series MIS john8.4.4 rnd100.1 hamming6.4 c125.9 mann.a27

f∗ 56 74 60 91 252

fcplex 56 76 60 94 -

fGAmipx 70 100 64 125 378

fGAgrd 56 75 60 93 296

Table 5. Summary of average ”gaps” Series RAN S N370e MIS

rcplex 0.13 % 0.1 % 8.78 % -

rGAmipx 0.003 % 0.03 % 6.71 % 30.83 %

rGAgrd 0.17 % 0.1 % 8.29 % 4.20 %

rGAnlo 1.32 % 4.35 % -

for series RAN, N370e and MIS, while for series S it was set to 500 to make the population converge slower. The inital value of pmut was set to 10−3 . In the GAgrd s = 20, R = 107 , N = 103 , pmut = 0.9, pcross = 0.8 for all of the test problems. The average number of iterations on series S was about 5.5 · 105 and on N370e it was near 1.7 · 106 . The best solutions found over all runs are reported in tables 1 – 4 for each of the algorithms. Here the best-known fbest or optimal f ∗ solution values are given as well. For the instances constructed using reductions the values fbest and f ∗ are obtained from the best-known solutions of the FCTP and the MAXIMUM INDEPENDENT SET problems. For series S the value fbest indicates the best known solution cost. The best values obtained in the described experiments are indicated in bold. As a summary, Table 5 gives the ”gap” values defined as r = (f /fbest − 1) · 100%, where f is the goal function value found by the algorithm being

considered. The ”gap” values are averaged over all problems in a series. Also, in column GAnlo we present the results of the GA NLO-1 from (Eckert and Gottlieb, 2002), especially developed for the FCTP. The run times of GA NLO-1 are unknown to us and its ”gaps” are available for instances ran4x64, ran16x16, ran14x18 and for series N370e. On each problem this GA was run 15 times until 106 non-duplicate solutions were produced. The experiments show that on the problems of small, medium and large size (series RAN, S and n370e) the GAmipx exceeds CPLEX and the GAgrd both in terms of relative error and the number of cases where the best values were reached. The best results for series n370e are obtained by the GA NLO-1, which is specialized for the problems of this type. On series MIS with the largest benchmarks the GAgrd turnes out to be the most robust, outperforming the other algorithms. On mann.a27 CPLEX 9.0 failed to find a feasible solution (most likely due to rounding errors). The GAmipx proved to be very inefficient on MIS, at least with our set of parameters. It is obvious that here the set of solutions was insufficiently covered by the MIP-crossover subproblems. On series MIS most of the subproblems proved to be infeasible, while on the rest of benchmarks this fraction was less than 10%. One of the ways to overcome such problem is to increase the threshold µ. However, this should be combined with usage of CPLEX in sub-optimal mode given a limited crossover time (otherwise the crossover will slow down the GA).

CONCLUSIONS The preliminary results reported here indicate that integration of a modern branch-and-bound MIP-solver into the crossover operator can significantly improve the performance of a GA. This approach may be competitive to a stand-alone branch-and-bound MIP-solver in terms of primary solutions quality but requires an adaptive parameter tuning algorithm. A more simple GA based on a greedy heuristic turns out to be a robust but not so precise method.

REFERENCES Aggarwal, C.C., J.B. Orlin and R.P. Tai (1997). An optimized crossover for maximum independent set. Oper. Res. 45, 225–234. Balas, E. and W. Niehaus (1998). Optimized crossover-based genetic algorithms for the maximum cardinality and maximum weight clique problems. Journ. of Heuristics 4(2), 107–122.

Barr, R.S., R.S. Glover and D. Klingman (1981). A new optimization method for large scale fixed charge transportation problems. Oper. Res. 29(3), 448–463. Borisovsky, P.A. (2004). Genetic algorithms for a supply management problem. In: Proc. Intern. Conf. ”Systems, Mechanisms and Machines Dynamics”. pp. 255–258. Omsk. (In Russian). Chauhan, S.S. and J.-M. Proth (2003). The concave cost supply problem. Europ. Journ. of Oper. Res. 148(2), 374–383. Chauhan, S.S., A.V. Eremeev, A.A. Kolokolov and V.V. Servakh (2005). Concave cost supply management problem for single manufacturing unit. In: Supply Chain Optimisation, Product/Process Design, Facility Location and Flow Control (O. Zaikin A. Dolgui, J. Soldek, Ed.). pp. 167–174. Springer Verlag. Chauhan, S.S., A.V. Eremeev, A.A. Romanova and V.V. Servakh (2004). Approximation of linear cost supply management problem with lower-bounded demands. In: Proc. of Discrete Optimization Methods in Production and Logistics (DOM‘2004). pp. 16–21. Omsk. Eckert, C. and J. Gottlieb (2002). Direct representation and variation operators for the fixed charge transportation problem. In: Proc. of PPSN 2002. pp. 77–87. Eremeev, A.V. (1999). A genetic algorithm with a non-binary representation for the set covering problem. In: Proc. of Operations Research (OR’98). pp. 175–181. Springer Verlag. Gottlieb, J. and L. Paulmann (1998). Genetic algorithms for the fixed charge transportation problem. In: Proc. of 5th IEEE Int. Conf. on Evol. Comp.. pp. 330–335. Klingman, D., A. Napier and J. Stutz (1974). Netgen: A program for generating large scale capacitated assignment, transportation and minimum cost flow networks. Manag. Sci. 20, 814–820. Kolokolov, A.A. and T.V. Levanova (1996). Decomposition and L-class enumeration algorithms for solving some location problems. Vestnik Omskogo Universiteta 1, 21–23. (in Russian). Reeves, C.R. (1997). Genetic algorithms for the operation researcher. INFORMS Journ. on Comput. 9(3), 231–250. Sun, M., J.E. Aronson, P.G. McKeown and D. Drinka (1998). A tabu search heuristic procedure for the fixed charge transportation problem. Europ. Journ. of Oper. Res. 106, 441–456. Zaozerskaya, L.A. (2003). A branch and bound algorithm for solving the concave cost supply management problem. In: Book of Abstr. of GOR Society Int. Conf.. p. 136.