Simulation Optimization in Uncertain Environments: An Evolutionary Approach

Proceedings of the 13th IFAC Symposium on Information Control Problems in Manufacturing Moscow, Russia, June 3-5, 2009

Simulation Optimization in Uncertain Environments: An Evolutionary Approach A. Baccouche*. H. Pierreval** A. L. Huyet*** *Research Unit of Technologies Information and Communication, Tunis, Tunisia (e-mail: [email protected]) **LIMOS, UMR CNRS 6158, IFMA, Campus Les Cézeaux, Aubière, France (e-mail: Henri.Pierreval @ifma.fr) ***LIMOS, UMR CNRS 6158, IFMA, Campus Les Cézeaux, Aubière, France (e-mail: [email protected])

Abstract: The design of manufacturing systems often involves simulation optimization approaches to search for the best system performance. Metaheuristic approaches are more and more used to optimize simulation models. In most existing approaches, the optimization is performed with fixed environmental conditions (e. g., part demand, breakdown rates are assumed to be perfectly known). However, in practice the actual data about the system may differ from those used in the simulation model (e. g., modification of the part mix). To cope with this, the candidate solutions, in a simulation optimization process, can be compared on various possible environments, using such principles as those proposed by Taguchi. Unfortunately, when metaheuristics are used, simulation optimization can be very much time consuming, since each solution has to be compared on a number of different environments. In order to provide robust solutions in a more reasonable time, we propose a two stage heuristic search. First, a restricted set of n promising solutions is identified using an evolutionary multimodal simulation optimization process, using the concept of base environmental scenario, recently published. Then robustness evaluation on many environments is performed only on these n promising solutions and the most robust can be chosen. This approach is illustrated on a supply chain problem where several parameters have to be defined. As a result, 13 solutions are found. The most robust solution is not the one that yields the best results in the environment assumed by the decision maker. This result shows how important it is to be able to consider several environments in the simulation optimization process. Keywords: Optimization via simulation, robustness, multimodal optimization, manufacturing systems.

1.

INTRODUCTION

In the design, organization or management of various types of systems, it is very frequent that a number of parameters have to be determined so as to optimize a given performance criterion. In this respect, simulation optimization approaches (Tekin and Sabuncuoglu, 2004), are used to search for the best possible values of a vector of input variables to a simulation model (e. g., size of a buffer, number of ressources and security stock level in a manufacturing system), so as to optimize a function (e. g., expectation) of an output variable (Rosenblatt and al, 1993). This can be formalized as follows: Min/Max f(X), (1) X∈ D where f is the output performance criterion, evaluated through simulation; X = (X1, …, Xn) is a vector that belongs to a domain D = ⊗ Di, where Di represents the respective domain of Xi. In most published works, the simulation model is built from data that characterize given and fixed environmental conditions (e. g., given operating times, transportation delay, etc.). Yet, in practice, disturbances can modify the values, and thus change the environment that was initially expected.

978-3-902661-43-2/09/$20.00 © 2009 IFAC

The expected environment, also called environmental base scenario in Pierreval and Durieux (2007), is the production environment that the domain experts think to be relevant to characterize the future operating conditions of the studied system. In practice, good performance achieved by the system in given conditions can be drastically deteriorated if these conditions change. For example, the production rate of a manufacturing system, evaluated by simulation, can be considered as acceptable, with regard to the customer demand, if breakdowns occur at a given rate, but can become unsatisfactory if the breakdown rate turned out to increase. To cope with this problem, it is essential to take into account uncertainty about data and possible changes in the environment when the performance of a solution is estimated through a simulation model. In this respect, in the optimization simulation area, a few articles suggest approaches based on such approaches as Taguchi’s methodology (Stuart Peace, 1992), response surface methodology (RSM) (Dellino and al, 2008), and robustness measurements based on referee curves (Pierreval and Durieux, 2002). However, when associated with a metaheuristic search, these approaches can require too much computing time. To cope with this problem, we suggest a new robust simulation optimization approach, which aims at taking into account

1504

10.3182/20090603-3-RU-2001.0429

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

several possible environments in the simulation optimization process in a more reasonable computing time. The article is organized as follows. First, related works are introduced and briefly discussed. Then, a novel methodology for robust simulation optimization problems is presented. This methdology uses multimodal optimization, which is adressed in the next section. Our conclusion and proposals for future research terminate the article.

For this they use a referee curve. However, their approach is also very time consuming, since in a first step they use an evolutionary simulation optimization to find the “best solution” and then a second evolutionary search that compares this solution on several environments, which again multiplies the number of evaluation needed for a single solution. 3.

2.

SUGGESTED APPROACH

RELATED WORKS

Robust optimization is a wide research area. The works presented in the literature address a large variety of problems. Beyer and Sendhoff (2007) identify several types of uncertainty involved in robustness optimization and several ways of formulating and addressing the corresponding problems. In this paper, we focus on simulation optimization problems. In this particular area, most works seems to be based on Tagushi’s works. He has addressed the problem of “parameter design” in the presence “noise factors”, which characterize the possible perturbations. He suggested the use of crossed experimental designs (orthogonal arrays) and a signal-to-noise (S/N) ratio or a quadratic loss function (QLF) to measure the robustness of solutions. Taguchian principles have been used for studying robustness of solutions using simulation, for instance in (Lim and al, 1996; Madu and Madu, 1999; Wild and Pignatello, 1991; Benjamin and al, 1995), and to address manufacturing problems, for instance in (Lim and al, 1996; Madu and Madu, 1999; Al Aomar, 2002). In the simulation area, Taguchi’s measures have been exploited by search methods, so as to try to find the most robust solution. Dellino and al (2008) propose to combine Taguchian principles and RSM for the robust simulation optimization of systems. RSM is well known in the simulation optimization area and is not as computationaly expensive as metaheuristics. Unfortunately, it can be quite sensitive to local extrema when the search space is complex (see Tekin and Sabuncuoglu, 2004; Pierreval and Tautou 1997). To avoid this sensitivity to local optima, a genetic search is used in (Pierreval and Durieux, 2002; Al Aomar 2002; Al Aomar, 2006; Yunker and Tew, 1994). Al Aomar (2002) suggests using the S/N ratio and the QLF as fitness functions in a genetic algorithm (GA). However, if this metaheuristic is known to be good at avoiding to be trapped, it is also known as time consuming. This shortcoming occurs in a context where the evaluation of solutions is performed using a simulation run (which often must be long or to be replicated (Kleijnen, 2008)), so that the evaluation of a single solution, say: CostSol, is particularly computationally expensive. Therefore, when cross experimental designs is needed to evaluate a single solution (in the presence of noise), this multiplies the computing cost of evaluating one solution CostSol by NS, the number of scenarios that are defined by the cross experimental design. As a result, the computing time can be prohibitive for certain type of applications. Another evolutionary approach is suggested by Pierreval and Durieux (2002). They do not use Taguchian principles. Their method aims at evaluating what would be lost, in terms of performance, if the operating environment changed and wheither this is acceptable for a decision maker.

The heuristic approach suggested also uses evolutionary computation, to avoid being trapped in local optima, but aims at reducing the computing cost. It relies on the assumption that if we can find a set Q= {S1,.., Sn} of very good solutions to the optimization problem (1), then one can find a robust solution among this set. In other words, in most cases it is reasonable to assume that a robust solution should also perform relatevely well in the expected operational conditions, and thus belong to Q. Let E be the given environmental base scenario. This expected environment is considered as the most relevant to represent the operating conditions in which the system will evolve (see Pierreval and Durieux (2007) for more details about this concept). The proposed approach is based on two stages. In a first stage, we aim at determing the set Q = {S1,.., Sn} of n so called promising solutions. A solution is called promising if it performs well in the conditions expected by the decision makers, i. e., in the given environmental base scenario E. In the second stage, we are interested in selecting the most robust of these n solutions by investigating their performances in various different experimental conditions. The changes can be caused by major modifications or by small disturbances (e.g., modification of the production mix, modification of the raw material delivery conditions, and higher breakdowns rates). Therefore, other possible environments E1, …., Ep have to be defined and implemented as other possible scenarios in the simulation model, so that each of the n solutions can be evaluated on each of the p environments. These environments can be determined using design of experiments approaches, such as Latin hypercube (see Kleijnen (2008) for more details). Certain solutions with excellent performance on the environment E can turn out to be less successful because of the noise that was taken into account in the evaluation of these solutions in the other environments. Several existing methods can be used to compare the n promising solutions in several environments, such as Tagushi’s S/N (Al Aomar, 2002), the referee curve (Pierreval and Durieux, 2002), etc. Once the robustness of every promising solution is evaluated, then we can pick the most robust solution of Q, denoted RS, as the final solution. Let us note that with this approach, the analysis of the solutions robustness is only performed for the n promising points, which are evaluated NS times. This contrasts with the existing evolutionary robust simulation optimization methods that compare every solution of the population NS times. As a consequence, the resulting computing costs are drastically

1505

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

reduced. The decision-makers will also have several possibilities to decide which configuration to adopt. A possible inconvenience of this approach would occur if the problem addressed was such that a given solution S turned out to be very robust but obtained poor results in E (so that it would not be selected in Q). However, since E is the environment that is assumed by the decision maker, this solution S would probably not be acceptable. We have now to focus on how to determine the set Q of promising solutions. Such a problem can be solved using a multimodal optimization search, which allows several local optima to be found during a single search. 4.

MULTIMODAL SIMULATION OPTIMIZATION

Multimodal optimization aims at locating multiple optima in the search space (Fig. 1). In our case solutions are evaluated using simulation on the expected environment and the n located optima will constitute the set Q of the n promising solutions. Several optimization metaheuristics can be relevant for this purpose, such as simulated annealing, tabu search, evolutionary algorithms, etc. The efficiency of multimodal evolutionary optimization has been reported by a number of researchers (e.g., Singh and Deb (2006); Petrowski (1996); Goldberg and Richardson (1987); De Jong (1975)). These approaches generally use niching methods (NMs) as strategies to drive the population of solutions towards several optima.

Singh and Deb, 2006; Qing Ling and al, 2008). The similarity measure is defined using comparisons between individuals in a same niche (i. e., distance metric between two individuals). We use a recently published method suggested by Qing Ling and al (2008), which combines a clustering strategy with a standard crowding method (De Jong, 1975). In fact, to form clusters (niches), the Crowding Clustering GA (CCGA) modifies the process selection proposed in the standard crowding and can eliminate selection errors. CCGA updates populations through competitions firstly between individuals in each cluster (i.e., the fittest solution in each cluster survives to enter the next generation) and secondly between clusters (i.e., aims to prevent multiple clusters from converging to a single optimal). By introducing the concept of peack detection condition proposed by Qing ling and al (2008) and maximising the size of the Crowding factor (CF) of niches, CCGA can avoid genetic drift effectively. The main steps of CCGA are as follows: 1. 2. 3.

4.

5. 6.

Fig. 1. A multi-modal function

7.

NMs aims at maintaining enough diversity in the population and at reducing the effect of genetic drift (e.g. convergence to a single peak) resulting from the selection operator in the standard GAs (Singh and Deb, 2006; Sareni, 1999). Niches divide the population into different subpopulations according to the similarity of the individuals. An efficient NM must be able to identify and maintain multiple diverse final solutions. Various NMs are considered in the literature, which are generally separated in two groups (Sareni, 1999). To form niches, the first group requires a predetermined threshold of similarity fixed for every niche such as the sharing scheme proposed by Goldberg and Richardson (1987) or the clearing procedure suggested by Petrowski (1996). In the second group, the threshold of similarity is not explicitly required. This is the case for such methods as those based on the crowding scheme initially proposed per De Jong (1975) (see

8.

Initialize distributed population. Number of population is PopNum. Recombine and mutate parents to generate PopNum children. Form niches : a clusters. For each parent Pj, j=1,.., PopNum, construct PopNum clusters {Pj, CSj} using a standard crowding model, with CF = PopNum, under a distance metric. For each cluster, select the fittest individual CCj as the center. For the others solutions in the cluster, calculate their distances from CCj, select the largest distance as the radius of cluster CRj. Sort all clusters according to their centers. Define the set of reserved clusters RC. Each element in RC has a center of cluster RCCi and a radius of cluster RCRi. For j=1…PopNum, compare the jth cluster with all clusters RCCi in RC, if for all i, D(CCj, RCCi) > RCRi or Peak (CCj, RCCi) = 1, then place CCj into RC and set the radius of cluster as min (CRj, D(CCj, RCCi)). D is the distance metric and Peak is the peak detection condition. Define the number of elements in RC as NRC, generate (PopNum - NRC). These individuals and the centers of clusters in RC enter the next generation. Repeat step 2 to step 7 until the maximum generation number (MaxGen) reaches.

More information about clustering, exploration and exploitation strategies of CCGA can be found in (Qing Ling and al, 2008). 5.

EXAMPLE

We illustrate our approach through a case inspired from the industry (Huyet and Paris, 2003), related to the relation between a customer (center of distribution) and a supplier, commonly refered to as supply chain loop. The system studied here, is constituted of a buyer group and a supplier which can produce three types of products: Pk for k = 1 to 3. Alorry brings the orders of the buyer group to the supplier

1506

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

and delivers the ordered products to the buyer group. According to the domain experts, the demand of products is decomposed as follows: 65% P1, 25% P2 and 10% P3. The supplier manufacturing system consists of a production line, which carry out three operations. Every operation requires a constant process time. There is a stock that is assumed to be infinite and a fixed number of kanban by product type. The warm-up period is 4000 hours and failures are not considered in the study. This system is presented in Fig. 2. For economical reasons, we have to minimize the costs induce by order-to-delivery lead time and the number of WIPs. The objective function to maximize can be expressed as follows: f(X) = M – Exp[(CWIP* NBE + ¦ Ck *TpsC of product k )] (2) With the given values: • M: large number,

• Ck: cost of one unity of time for orders delay of product k (constant), • TpsCk:time to satisfy a command of product k (random), • k : type of product (constant), In our example, the problem to address is to determine values for the following parameters which characterize X: • KBk =[10…30]: number of kanbans in front of the production line for the product of type k, 3 variables. • TLUk =[1…20]: transport load units associated to product k. • PLUk = [1…10]: production load units associated to product k. Thus, the design problem consist in determining nine parameters, characterized by the vector X = [KBk, TLUk, PLUk] for k = 1 to 3.

• Exp(Z): expectation of Z, • NBE: number of WIP (random), • CWIP : cost of one WIP (constant),

DISTRIBUTION CENTER

SUPPLIER

S T O C K

Process 3 Process 2

Customers orders

Transport Lots

Process1

Delivered Products Kanbans

Production lots Fig. 2. Supply chain loop

Table 1. Promising solutions on E Design Solution X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13

Performance KB1 KB2 KB3 TLU1 TLU2 TLU3 PLU1 PLU2 PLU3 99692 99684 99682 99679 99675 99672 99666 99648 99646 99609 99599 99565 99551

16 14 14 15 20 15 21 10 21 14 26 16 20

12 12 12 13 12 13 12 12 13 11 13 12 12

10 10 10 10 10 12 10 10 12 13 12 10 10

14 14 14 14 17 14 14 14 14 19 14 17 17

1507

10 10 10 10 10 10 10 10 10 5 10 17 17

4 4 4 4 4 4 4 9 4 3 4 14 14

5 5 5 5 5 5 5 5 5 5 5 9 9

4 4 4 4 4 4 4 6 4 4 4 5 5

10 9 10 10 10 9 10 10 9 9 9 8 8

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

The evaluation of the function f is made using simulation. The duration of the simulation is 80000 hours (long run see (Kleijnen, 2008)). The settings of the multimodal GA are: PopNum = 30 and MaxGen = 200. We adopt the linear crossover with random parameter for real coded (Qing Ling and al, 2008). The Euclidean distance as similarity measure and CF = PopNum to form clusters. We obtained 13 promising design solutions (see Table 1). In the expected environment E, the best solution found is number 1. Our system is faced with variation in the market demand. We consider here that the demand ratio of each product is uncertain and subjected to change (uncontrollable environmental factors) and can affect system performance. Table 2. Environments tested

E E1 E2 E3 E4 E5 E6

P1 0,65 0,8 0,65 0,25 0,25 0,65 0,7

P2 0,25 0,1 0,3 0,25 0,45 0,15 0,25

P3 0,1 0,1 0,05 0,5 0,3 0,2 0,05

We are interested in studying other environments with different ratios of orders. Six unexpected environments are considered in the robustness study, as given in Table 2. To evaluate and compare solutions, we run simulations for each promising solution in each unexpected environment and compute the S/N ratio (see Table 3). Using Tagushi’s principles, our problem is a maximisation (2), so we use the so called “larger-the-better S/N”, as follows: p 2 S/N = − 10 ⋅ Log ( 1 / p ¦ 1 / y ) (3) i ij y=1

where p is the number of environments tested for each solution and yij is the performance measure of the solution i under the environment j. The S/N computed in Table 3 shows that the solution X9 is the most robust and is therefore more attractive than X1 for a decision maker who considers that the demand ratio can change. In this example, P2 and P3 have much longuer processing time than P1, so that if their production increase the system performance is modified and the design X9 is more suited than X1. It is mainly in environments E3 and E4, that the deterioration is the wirst. 6.

CONCLUSION

In such problems as manufacturing systems design, simulation optimization is recognized as an important method to determine the right system parameters. In a more and more uncertain world, we have, presented a new simulation optimization approach that aims at finding robust solutions. The presented approach keep the well known benefits of evolutionary methods to avoid being trapped in local extrema, and do not require every solution to be evaluated on each environment. The example given has shown that the best solution found using a “classical” simulation optimization search can not be the most relevant for a decision maker. This two stage method provides several possible solutions, which is useful for the decision maker. Several methods can be used in the second stage to select the final solution, depending on the problem objectives and constraints. Although the proposed approach allows reducing the computing time with regard to other existing robust simulation optimization approaches, the use of evolutionary computation requires many solutions to be evaluated.

Table 3. Performances of promising solutions Design F(Xi ,E) Solution X1 99692 X2 99684 X3 99682 X4 99679 X5 99675 X6 99672 X7 99666 X8 99648 X9 X10 X11 X12 X13

99646 99609 99599 99565 99551

99682,62 99679,49 99678,16 99671,86 99666,63 99662,27 99654,13 99654,76

99695,22 99689,78 99688,91 99682,24 99673,75 99671,02 99664,6 99655,23

84478,91 83732,33 82781,47 83626,42 85104,03 90759,07 85654,4 83345,53

97038,61 95183,24 95594,1 96526,93 97613,32 98270,75 97669,03 93711,13

99673,11 99655,22 99648,24 99653,9 99664,58 99656,37 99655,43 99606,44

99691,35 99689,14 99687,64 99679,81 99670,96 99669,08 99659,89 99656,37

Ratio S/N -0,29 -0,33 -0,34 -0,43 -0,32 -0,17 -0,26 -0,36

99635,23 99627,73 99591,58 99564,19 99550,17

99642,37 99621,8 99594,47 99567,7 99553,99

92322,34 73920,69 92272,44 92070,7 92275,69

99062,79 83610,93 99012,82 98577,13 98665,65

99646,47 99464,77 99601,32 99556,23 99543,49

99639,33 99626,59 99591,96 99570,74 99556,56

-0,13 -0,74 -0,14 -0,15 -0,19

F(Xi ,E1) F(Xi ,E2) F(Xi ,E3) F(Xi ,E4) F(Xi ,E5) F(Xi ,E6)

1508

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

Our future research directions include the distribution or parallelization of this approach (Pierreval and Paris, 2000) so as to reduce the computing costs. REFERENCES Al–Aomar R. (2002). A robust simulation based multicriteria optimization methodology. Proceedings of the 2002 Winter Simualtion Conference, pp. 1931-1939. Al –Aomar R. (2006). Incorporating robustness into Genetic Algorithm search of stochastic simulation outputs. Simulation Modelling Practice and Theory , 14, pp. 201– 223. Benjamin P. C., Erraguntla M., and Mayer R. J. (1995). Using simulation for robust design. Simulation, 65 (2), pp. 116-127. Beyer, H.-G. and Sendhoff, B. (2007). Robust optimization – A comprehensive survey. Computer Methods in Applied Mechanics and Engineering, 196, pp. 3190–3218. De Jong, K.A. (1975). An Analysis of the behavior of a class of genetic adaptative systems. PhD thesis, University of Michigan. Dissertation Abstracts International, 36 (10), 5140B, university Microfilms n° 76-9381. Dellino G., Kleijnen, J.P.C. and Meloni C. (2008). Robust optimization in simulation: Taguchi and Response Surface Methodology. CentER Discussion Paper, 200869, Tilburg University, Netherlands. Goldberg D.E. and Richardson, J. (1987). Genetic algorithms with sharing for multimodal function optimization. Proceedings of the 2nd International Conference on Genetic Algorithms, pp. 41-49. Huyet A.-L. and Paris J.-L. (2003). Setting and analysis of multiproduct supplier links managed by kanban using a learnig+ evolutionnary optimization. Proceedings Of ISC 2003, the Industrial Simulation Conference, pp. 113-117. Kleijnen, J.P.C. (2008). Design and analysis of simulation experiments. Springer-Verlag, Heidelberg. Lim J. M., Kim K. S., Yum B. J. and Hwang H. (1996). Determination of an optimal configuration of operating policies for direct-input-output manufacturing systems using the Taguchi method. Computers Industrial Engineering, 31 (3/4), pp. 555-560. Madu I.E., Madu C. N. (1999). Design optimization using signal-to-noise ratio. Simulation theory and Practice, 7, pp. 349-372. Petrowski A. (1996). A clearing procedure as a niching method for genetic algorithms. Proceedings of the IEEE International Conference on Evolutionary Computation, pp. 798-803. Pierreval H. and Tautou L. (1997). Using evolutionary algorithms and simulation for the optimization of manufacturing systems. IIE Transactions, 29 (3), pp. 181190.

Pierreval H. and Paris J. L. (2000). Distributed evolutionary algorithms for simulation optimization. IEEE transactions on systems man and cybernetics, 30 (1), pp. 15-24. Pierreval H. and Durieux S. (2002). A two-stage simulation optimization heuristic for robust design. Man and Cybernetic, Proceedings of IEEE System, 6 pages. Pierreval H. and Durieux S. (2007). Robust simulation with a base environmental scenario. European Journal of Operational Research, 182, pp.783–793. Qing Ling., Gang, W., Zaiyue, Y. and Qiuping, W. (2008). Crowding clustering genetic algorithm for multimodal function optimization. Applied Soft Computing, 8, pp. 8895. Rosenblatt M.J., Roll Y. and Zyse V. (1993). A combined optimization and simulation approach for designing automated storage/retrieval systems. IIE Transactions, 25 (1), pp. 40-50. Sareni, B. (1999). Méthodes d’optimisation multimodales associées à la modélisation numérique en électromagnétisme. PhD. thesis, Ecole Centrale de Lyon, France. Singh, G. and Deb, K. (2006). Comparison of multi-modal optimization algorithms based on evolutionary algorithms. ACM, GECCO, Seattle, Washington, USA. Stuart Peace G. (1992). Taguchi Methods - A Hands-On Approach. Addison-Wesley Publishing Company, INC, 224 pages. Tekin E. and Sabuncuoglu I. (2004). Simulation optimization : A comprehensive review on theory and applications”. IIE Transactions, 36, pp. 1067-1081. Wild R. H. and Pignatello J. J. (1991). An experimental design strategy for designing robust systems using discrete event simulation. Simulation, 57 (6), pp. 358-368. Yunker, J.M. and Tew, J.D. (1994). Simulation optimization by genetic search. Mathematics and Computers in Simulation, 37 (1), pp. 17 - 28.

1509