An improved simulated annealing simulation optimization method for discrete parameter stochastic systems

Computers & Operations Research 32 (2005) 343 – 358

www.elsevier.com/locate/dsw

An improved simulated annealing simulation optimization method for discrete parameter stochastic systems Scott L. Rosen∗ , Catherine M. Harmonosky The Harold and Inge Marcus Department of Industrial and Manufacturing Engineering, The Pennsylvania State University, 310 Leonhard Building, University Park, PA 16802, USA

Abstract This paper proposes a new heuristic algorithm for the optimization of a performance measure of a simulation model constrained under a discrete decision space. It is a simulated annealing-based simulation optimization method developed to improve the performance of simulated annealing for discrete variable simulation optimization. This is accomplished by basing portions of the search procedure on inferred statistical knowledge of the system instead of using a strict random search. The proposed method is an asynchronous team-type heuristic that adapts techniques from response surface methodology and simulated annealing. Testing of this method is performed on a detailed simulation model of a semi-conductor manufacturing process consisting of over 40 work-stations with a cost minimization objective. The proposed method is able to obtain superior or equivalent solutions to an established simulated annealing method during each run of the testing experiment. ? 2003 Elsevier Ltd. All rights reserved.

1. Introduction to the simulation optimization problem Analysis of large and complex stochastic systems is a di4cult task because of the complexities that arise when randomness is embedded within a system. Unfortunately, unexplained randomness is a common and unavoidable characteristic among real-world systems. The emergence of simulation modeling as an evaluative tool for stochastic systems has facilitated the ability to obtain performance measure estimates under any given system con6guration. Furthermore, these estimates are much more accurate than when they are estimated by analytical techniques that make unrealistic assumptions about the system. Computer simulations are highly e8ective in answering evaluative questions concerning a stochastic system. However, it is often necessary to locate values of the model decision parameters such that ∗

Corresponding author. Tel.: +1-814-861-4604. E-mail address: [email protected] (S.L. Rosen).

0305-0548/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0305-0548(03)00240-5

344

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

a performance measure is maximized or minimized. This type of problem is referred to as the simulation optimization problem. The simulation optimization problem can be de6ned in the context of a single objective optimization problem in Eq. (1): min:

f()

s:t:

 ∈ ;

(1)

where f() = E[(; !)] is the expected value of the system performance measure estimated by ˆ f(), which is obtained from j sample performances j (; !) of a simulation model observed under an instance of discrete or continuous feasible input parameters  constrained within a set of feasibility  ⊂ Rd . The stochastic e8ects of the system are represented by !. The di4culty of this problem is that f() is an implicit function of the decision parameters where an observation of f() can only be obtained through an execution of the simulation model. Since simulation runs are often expensive, the convergence rate of an algorithm must be considered in addition to the quality of solution generated. In addition, automation should also be considered when developing a simulation optimization technique because it is impractical for a simulation end user to be a participant of a lengthy optimization process. Simulation optimization techniques can be broken down into the following classes: gradient-based search methods (most notably stochastic optimization and response surface methodology), sample path optimization, random search heuristics, simplex search methods, and statistical methods. Gradient-based search methods, simplex search methods and sample path optimization are generally applied to stochastic systems consisting of strictly continuous decision variables. Random search heuristics are commonly used to optimize stochastic systems that contain discrete decision parameters. Statistical methods are also applied when the decision space is discrete and are preferred over the random search methods when the cardinality of the set of alternatives is small. 2. Simulation optimization with simulated annealing This paper considers the discrete parameter simulation optimization problem. The general problem formulation is provided below: min:

f()

s:t:

 ∈ ;

(2)

where  = {1 ; 2 ; : : : ; i ; : : : ; n } and i ∈ Z∀i where Z is the set of integer numbers. The strategy used for optimizing a discrete parameter simulation model is dependent on the size of the solution space. If the solution space is small, for example, 20 or fewer decision solutions, multiple comparison procedures and ranking and selection procedures are often appropriate [1]. However, simulations of systems containing large discrete decision variables are the focus of this research. Discrete parameter systems include manufacturing systems, facility location con6guration, worker allocation and military systems. Simulated annealing is one of the more commonly used heuristics for discrete parameter simulation optimization when the solution space is large. Simulated annealing (SA) is a heuristic search procedure 6rst introduced by Kirkpatrick et al. [2] for complex deterministic optimization problems with a discrete solution space. Bulgak and

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

345

Sanders [3] propose the 6rst heuristic procedure for SA that is used for simulation optimization to optimize bu8er sizes in an automated assembly system. They use con6dence intervals to account for the randomness in the response, which were incorporated into the calculation of acceptance probability of the inferior candidate solutions. Haddock and Mittenhall [4] discuss a slightly di8erent application of SA to simulation optimization. Instead of using con6dence intervals, they use a point estimate of the expected value of the response and employ SA in accordance to Kirkpatrick et al. [4]. The di4culty here is when comparing two solutions based on estimates of expected value one must be careful to apply an appropriate estimate and to run su4cient number of simulations. So and Dowsland [5] consider applying a less accurate estimate of expected value in order to reduce the number of simulation runs. The inaccuracy of the estimates was acknowledged in the calculation of the acceptance probability. More recent SA heuristics for simulation optimization include an adapted SA heuristic to consider constrained discrete parameter simulation optimization problem [6]. Alkhamis et al. [7] present a modi6ed SA algorithm, where the likelihood of an improving move is based on whether the con6dence intervals on the responses under comparison show statistical signi6cance between the expected value of the responses. Alrefaei and Andradottir [8] also consider two SA approaches that use a constant temperature. Most recently, Ahmed and Alkhamis [9] integrate the methods of SA with ranking and selection. 3. Focus of this research SA applies nicely to the discrete parameter simulation optimization problem since it can be performed on any standard “black box” optimization problem and has the capability to locate a satisfying local optimum or a possible global optimum in complex combinatorial optimization problems by avoiding unwanted local optimums. However, since SA is a search based on randomness and no inferred knowledge of the response surface, SA’s convergence to a 6nal solution is often slow. Slow convergence of SA results in a problem when it is applied to the simulation optimization problem because many expensive simulation runs will be required. This paper presents a new SA approach to discrete parameter simulation optimization, which incorporates statistical knowledge of the response to improve its convergence rate. The proposed technique e4ciently performs repeated searches for multiple high quality local optimal solutions by using learned knowledge of the response surface, instead of blindly searching throughout the feasible region. This proposed method is tested against an established SA method for simulation optimization, comparing 6nal solution values and convergence rates. 4. Overview of the proposed method This proposed simulation optimization heuristic technique is referred throughout the paper as the RS team method. The overall methodology behind the RS team method is to search various high quality local optimal solutions of the feasible region with the combined use of 6rst-order linear approximations and SA. The proposed methodology is broken down into two phases. Phase I of this algorithm consists of a search process using linear approximations. The current solution after phase I

346

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

is most likely a solution of high quality. This implies that other high quality solutions could exist in adjacent and nearby neighborhoods. Therefore, these neighborhoods are explored in the next phase of the algorithm. In phase II of the search procedure it is desired to explore a small subset of the feasible region around the phase I 6nal solution in order to locate a solution of possible higher quality that was unreachable by the search based on linear approximations. However, we do not want to move too far away from the phase I solution so that the simulation runs used to locate the phase I 6nal solution are wasted. It should be noted that the phase I solutions are possible local optima solutions. SA works well for phase II because it has the ability to move away from a local optima Phase I solution in order to explore adjacent and nearby neighborhoods. SA also works well for phase II because its parameters can be calibrated to ensure that the search does not venture too far from the phase I solution. Repeated searches for high quality local optimal solutions throughout the feasible region are made by performing phases I and II from di8erent starting solutions. The number of searches and the distance between the starting point of each search can be calibrated in this algorithm to ensure that the majority of the feasible region is explored. The 6nal solution is the best solution obtained from all of the searches. The techniques for constructing linear model approximations and searching along directions of improvement related to the linear models are very similar to techniques used in response surface methodology. However, for the implementation of response surface methodology, it is necessary that each decision variable lies in a continuous space. Therefore, we suggest modi6cations to make it applicable to the problem with a discrete decision space. These modi6cations are discussed later in the paper. The SA phase of the algorithm is similar to Haddock and Mittenhall [4], which modi6es Kirkpatrick et al. [2] by substituting an estimate of the expected value of the response in all places requiring a deterministic objective function value. The techniques utilized to obtain the linear approximations are adaptations of the 6rst phase of Box and Wilson [10]. Hood and Welch [11] also provide detailed background on applying response surface methodology to computer simulation. 4.1. Procedure This section provides formal de6nitions of the parameters used in the RS team method, followed by a list of steps involved in the RS team method. Each step is then discussed in greater detail. The procedure is described in the context of a minimization problem. Parameters m Cs Ns Tl

current step of the algorithm meaning that m − 1 searches have already been attempted the best current solution within step m. The value for Cs is updated after every line search in phase I. The value for Cs is also updated during SA in phase II the set of all neighboring solutions to Cs = {ns :ns is in the neighborhood of Cs } in phase II temperature parameter at iteration l of phase II—e8ects the probability of accepting a non-improving move during SA portion of the algorithm

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

C   

347

cooling parameter—the rate by which the temperature is decremented during the SA portion of the algorithm under a geometric cooling schedule SA stopping parameter—the maximum temperature value where |Ns | consecutive non-improving moves will result in the end of an SA phase distance parameter—the minimum allowable distance between a new starting point and a phase II 6nal solution global optima stopping parameter—the number of searches attempted Initialization

(1) Determine the initial temperature parameter T0 , cooling schedule parameter C, distance parameter , global optima stopping parameter , and the SA stopping parameter . (2) Construct a feasible region for the decision variables by 6rst specifying a minimum (Li ) and a maximum (Ui ) value for each variable i = 1 to n  Li 6 i 6 Ui ∀i. Add any other necessary constraints. (3) Determine an initial starting point at random. Set the simulation response to the initial starting point to Cs . Phase I (4) Fit a linear approximation model around the initial starting point with the use of a 2k factorial design. (5) Advance to a point that is along a direction of improvement in accordance to the linear model and obtain an observation of f() from the integer point closest to the point of advancement. If the response after the 6rst step results in a poorer solution then go to step 7, else set Cs to the best integer solution found along this direction of improvement and go to step 6. (6) Construct a new linear approximation around the stopping point and go to step 5. Phase II (7) Search around a small area of the stopping point using SA until  consecutive non-improving moves occur and all neighbor solutions result in a poorer solution (Cs ¡ ns ∀ns ∈ Ns ). (8) Increase m by one. If the number of searches made is less than  then go to step 9. Otherwise, terminate and record the global optimal solution as the best phase II 6nal solution across all searches. (9) Generate a new random starting point in the feasible region that is distance  away from all phase II 6nal solutions found thus far. Return to step 4. The algorithm begins with the setting of algorithm parameters followed by the user de6ning a feasible region and a random selection of an initial starting point (steps 1–3). The remainder of the algorithm is composed into two phases. Phase I of the algorithm constructs linear approximations and searches through the feasible region using improvement directions generated by local linear models. Phase II uses SA to search the adjacent neighborhoods of the phase I stopping point for solutions of higher quality. Phase I begins by generating a standard experimental design that can provide a least-squares estimation of a linear model (step 4). An orthogonal 6rst-order design centered about the starting

348

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

point is recommended to achieve model coe4cient estimates of minimum variance. The discrete nature of the feasible region must be considered when selecting an appropriate experimental design for generation of the linear models. A su4cient experimental design is the 2k -factorial design and its orthogonality property is not a8ected by the discrete nature of the decision space. Fractions of a 2k -factorial design are encouraged if the user is willing to sacri6ce the probability of obtaining a global optima for a faster convergence rate. After exercising the simulation model at each experimental point of a full 2k design, a vector of responses y = [f()1 ; f()2 ; : : : ; f()2k ] is obtained from the simulation model corresponding to a matrix X representative of the levels of  utilized at each experimental point. The coe4cients to the linear model  are estimated by the least-squares equation in Eq. (3). The resulting linear model g is provided in Eq. (4): ˆ = (X
X )−1 X
y; g = ˆ0 +

n


ˆi i ;

(3) (4)

i=1

where ˆ is a vector of least-squares estimators and the 6rst column of X is a vector of ones. The search for and area of high quality solutions begins by moving in a decent direction of g (step 5). If the feasibility of the direction is not considered and the discrete nature of the feasible region is relaxed, the direction d of maximum decrease of g simpli6es to Eq. (5): −[ˆ1 ; ˆ2 ; : : : ; ˆn ] d=  : ˆ21 + ˆ22 + : : : ˆ2n

(5)

A bene6cial move is to advance from the current solution Cs to the feasible integer solution of minimal distance to the optimal value of f along the line generated by d. This solution is obtained by 6rst locating a value of t that minimizes f(Cs + td) such that Cs + td ∈  and t ¿ 0. Then the new current solution is obtained by taking the closest integer value of each element of Cs + td. Due to the stochastic nature of this implicit function and the computational intensity in obtaining just an observation, there is not a robust procedure for the above operation. Instead, we explore in small steps along in the direction d and take an observation of f at each step through execution of the simulation model. The step sizes of each i obtained by d are kept to a continuous value for each variable, but each variable is rounded o8 to its closest integer value before input into the simulation model. To determine the step sizes for each index of , the user 6rst decides the step size for the i with the largest corresponding value of ˆi . The step-sizes of the other variables are computed by Eq. (6). Oi =

ˆi

(ˆl =Ol )

;

(6)

where Ol is the step size for variable i chosen by the user and Oi is the step size for all other variables i . In order to be economical with the number of simulation runs executed, the search along d simply terminates when a resulting response is worse than the previous step or is infeasible. A new value

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

349

of Cs is then set to the last value of  which resulted in an improvement. A new linear model g is then constructed around Cs by generating a new experiment centered about the point, which resulted in the best solution from the previous search (step 6). Another search is then performed along a new d obtained from terms of the new linear model. These searches are repeated until the response after the 6rst step of a search along d is inferior to the current solution or is infeasible indicating that local linear approximations cannot incorporate enough knowledge of f() to generate a direction of improvement. Phase II begins using SA (step 7) with an initial temperature value T0 that results in a signi6cant probability of accepting a non-improving move in the beginning of SA in order to promote a small exploration of the feasible region from the current solution obtained after phase I. Phase II will also use a value of C that enables quick convergence so that we do not move too far away from the current solution and make waste of the expensive simulation runs used for building the linear approximations. Appropriate values for T0 and C to accomplish this are discussed in more detail later in the paper. The constraints on the feasible region are incorporated during the SA phase by simply transforming f() to a penalty function where an infeasible point will result in a substantially large value of f(). After a phase II 6nal solution is found, a new starting point is generated randomly, which is at least a predetermined distance away from all previously found phase II 6nal solutions (step 9). This is speci6ed by the distance parameter , which is discussed in the next section of the paper. Repeated searches are performed a user-speci6ed number of times and is determined by the global optima stopping parameter , which is also discussed later in the paper. The best solution found from each of the repeated searches is the 6nal solution. 4.2. Generating new starting points The main concept behind this algorithm is to search out di8erent high quality local optimal solutions through the feasible region, with starting points generated in all areas of the feasible region. It is important that a starting point is not generated within a close area of a phase II 6nal solution to avoid convergence at a previously found solution. That would be a waste of simulation runs and time. The value for each decision variable in the starting solution should be generated individually at random, such that every possible value of the decision variable has equal probability of being selected. To ensure that a starting solution s is located an appropriate distance from a phase II 6nal solution a, the sum of the scaled distances Das between the values of the decision variables at the starting solution and the values of the decision variables for a phase II 6nal solution is calculated by Eq. (7): n |ai − si | Das = i=1 ; (7) Ri where n is the number of decision variables, Ri is the range of decision variable i in the solution space, ai is the phase II 6nal solution value for decision variable i, and si is the value of the starting solution for decision variable i. Das should be calculated for each of the phase II 6nal solutions obtained thus far. If each calculated value of Das is greater than the distance parameter value set by the user, then this solution is an

350

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

acceptable starting point. If any calculated value of Das is less than the distance parameter value, then another potential starting point must be generated. 4.3. Termination criterion The RS team method allows the user to control the number of iterative searches that employ both Phases I and II. The user speci6es this when con6guring the global optima stopping parameter  prior to initiation of the algorithm. One iterative search is de6ned by generating a starting solution and executing phases I and II of the algorithm. Each iterative search concludes at the end of phase II when  consecutive non-improving moves occur and all neighbor solutions result in a poorer solution. If a search converges on a phase II solution that has already been found then the search is still counted. The algorithm will terminate after  searches have been performed. 5. Implementation and experimentation To test the performance capabilities of the proposed methodology, an experiment was designed to determine if the RS team method can obtain comparable or superior solutions to SA in fewer simulation runs on various types of systems. SA in accordance to the version of Haddock and Mittenthal [4] was chosen as a comparative algorithm because it is applicable for discrete variable stochastic systems and past research and testing has been previously performed to justify its credibility. 5.1. The testing model The RS team method and SA are tested on a detailed simulation model of an actual semi-conductor manufacturing process containing 399 processing steps, 41 machines, and 50 laborers. More specifically the simulation model pertains to the manufacturing process of a high density interconnect module. The majority of the process is performed across three areas: Carrier Fabrication, Pattern Layering, and Surface Mounting. Within each of the three areas in the modeled system there exists a series of processing modules formed between inspection points, which are comprised of several processing steps. In the pattern layering section of the process, parts pass through these modules repeatedly to create the multiple layers of the end product. Fig. 1 illustrates the system in minor detail. Refer to Harmonosky et al. [13] for a more detailed description of this process model. Due to the size, complexity and stochastic nature of this system, it is best con6gured via simulation. The performance measure of interest for this system is to minimize the cost to manufacture one part. The cost per part objective function is comprised of four major components: work in process cost, labor cost, machine cost and tardiness cost. Based on some preliminary analysis and discussions with managers of the facility, there are 6ve decision variables of the system considered to have a signi6cant e8ect on the cost per part objective: the inter-arrival rate, the number of general workers, the number of class I workers, the number of class II workers, and the number of B machine resources used in the system. Each solution of the system is an allowable parameter setting within the system and the solution space size is de6ned as the number of di8erent possible parameter settings within the system.

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

351

Carrier Preparation

Carrier Fabrication

Pattern Layering

Surface Mounting

Repeat for each layer

Electrical Test

Packaging and Shipping

Fig. 1. Process Pow diagram of testing system.

Although most systems may not have control over their inter-arrival rate of parts into their system, this semi-conductor manufacturing process did have complete control over their inter-arrival rate by using releasing strategies. Since the inter-arrival rate of parts is known to have such a signi6cant e8ect on the total production cost per part it is included as a decision variable. Class I workers, class II workers and the general workers are di8erentiated by the speci6c tasks each group may perform. The B machine resource capacity was included as a decision variable because it is utilized in many of the processing steps and is a potential bottleneck. 5.2. Experimental design A 22 experimental design is used to test the e8ectiveness of the proposed method versus SA. The two factors are the number of decision variables in the system and the solution space size of each of the variables. These two factors are chosen for the experiment because they each have the potential to have a signi6cant impact on convergence time and goodness of solution for each method. This 22 experiment yields four design points, as shown in Fig. 2. 5.3. Setting parameter values (RS team method) The RS team method consists of 6ve parameters that must be con6gured upon initialization: initial temperature (T0 ), cooling rate (C), SA stopping parameter (), global optima stopping parameter () and the distance parameter (). To achieve a good performance from this algorithm, these parameter values must be con6gured properly. This section discusses the con6guration for the values and ranges for each of the 6ve parameters. The parameters T0 ; C, and  a8ect the performance of the SA portion of the algorithm. As discussed earlier, it is important to con6ne the SA search around a small area of the phase I stopping point while ensuring that it leaves the neighborhood of the phase I stopping point. To encourage SA to explore away from the neighborhood of the phase I stopping point, the value of T0 is selected so that the resulting initial probabilities of selecting a non-improving move are relatively high. This research uses Eq. (8) for computing the probability of selecting a non-improving move

352

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358 # OF DECISION VARIABLES 2

5

(Design Point 1)

(Design Point 3)

X1: inter-arrival rate, (1000-1500 min)

X1: inter-arrival rate, (1000-1500 min)

X2: # of general workers, (5-10 workers)

X2: # of general workers, (5-10 workers)

SMALL

X4: # of class II workers (1-5 workers) X5: # of B machine resources (1-5 machines) (Design Point 2)

(Design Point 4)

X1: inter-arrival rate, (0-5000 min)

X1: inter-arrival rate, (0-5000 min)

X2: # of general workers, (3-20 workers)

X2: # of general workers, (3-20 workers) X3: # of class I workers (1-10 workers)

LARGE

SOLUTION SPACE SIZE

X3: # of class I workers (1-5 workers)

X4: # of class II workers (1-5 workers) X5: # of B machine resources (1-5 machines)

Fig. 2. Experimental design points.

at iteration l of the SA phase: e−|((zc −zs )=zs )×100|=Tl ;

(8)

where zc is the objective function value of the candidate solution, zs is the objective function value of the current solution, and Tl is the temperature value at iteration l of SA. This formula [12] uses a percent deviation between the current and candidate solution to calculate the probability of accepting a non-improving move. The advantage of Eq. (8) compared to using only (zc − zs ) in the numerator of the exponent is that the exponent term is a dimensionless value making it independent of problem speci6cations. Therefore, the initial temperature can be chosen so that a desired probability for accepting a speci6c deviation in current and candidate solutions can be obtained for the 6rst iteration. An initial temperature T0 is selected to allow for an inferior solution of 65% to be accepted 90% percent of the time. Rearranging Eq. (8) and substituting T0 for Tl results in Eq. (9): T0 = −&=(ln P);

(9)

where & is the speci6c deviation between current and candidate solutions, such that, & = |((zc − zs )=zs ) × 100| and P is the desired probability of accepting the above deviation. There are di8erent cooling or temperature reduction schedules that can be used for simulated annealing. Each type of cooling schedule has a di8erent e8ect on the behavior of the SA search process. A geometric cooling schedule is chosen for this algorithm. With a geometric cooling schedule, the resulting temperature value Tl from any iteration l is obtained from Eq. (10): Tl = T0 (C)l ;

(10)

where Tl is the temperature at iteration l and C is the cooling parameter chosen by the user. The reason that a geometric cooling schedule is chosen for this algorithm is because it allows for a quick drop in temperature in the beginning SA process followed by increasingly smaller temperature

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

353

drops as the SA process continues. This results in having most of the SA search process at a low temperature, which can keep the exploration from going too far away from the phase I stopping point. A further advantage of the geometric cooling rate is that since the temperature drops largely in the beginning, a high initial temperature can be used. Again, the goal of phase II is to be able to explore around a small area of the Phase I stopping point with a good likelihood of being able to venture away from the neighborhood of phase I. The use of a geometric cooling schedule enables SA to do this. Other common cooling schedules including linear probability cooling and linear temperature cooling do not allow phase II of the algorithm to behave as noted. Logarithmic cooling, however, which has been shown to converge almost surely to a global optima over any subset of  [14], could be calibrated to make SA work e8ectively in this algorithm because it also allows for a high initial temperature and a majority of iterations at high temperature values. The temperature at step l resulting from a logarithmic cooling schedule under a cooling parameter C is given in Eq. (11): Tl =

C : (1 + log(l))

(11)

The logarithmic cooling schedule is a suitable alternative to the geometric cooling schedule for this algorithm. It was not chosen for this algorithm because its a8ects on SA performance are much more sensitive to the value of C chosen by the user than the geometric cooling schedule [15]. Also, though there are mathematical results that show necessary conditions for SA to converge almost surely to a global optima of any subset of  under the use logarithmic cooling, these results provide no extra belief that this optima can be reached any quicker than an exhaustive search. It has been shown through example that convergence times with a logarithmic cooling schedule are exponential for even a simple problem [16]. For both of the two-decision variable experiments, the T0 is set to 5 and C is set to 0.95. In the 6ve-decision variable experiments, the T0 is set to 3 and C is set to 0.95. The SA stopping parameter  is used to determine when to terminate the SA process in phase II and to accept the best solution in the current search as the phase II 6nal solution for the current search. The SA phase ends when every ns ∈ Ns results in a poorer response than Cs and the temperature parameter has dropped below the value of . It is  in conjunction with C and T0 that controls the number of minimal temperature stages and thus the minimal extensiveness of the SA search. Thus, the number of desired minimal temperature stages are 6rst determined and  is set to achieve the number for the previously established values of T0 and C [17]. For each experiment,  is set to a value of 1. This results in an increased allocation of minimal temperature changes for the experiments as the solution space increases. Setting the distance parameter is a subjective task. Intuitively,  should be increased as the overall solution space size increases due to the expanded area of the feasible region that needs to be explored. A minimum value of 0.5 for  should be used [17]. For the two-decision variable, small solution space experiment,  is set to 0.5. For the two-decision variable, large solution space experiment and the 6ve-decision variable, small solution space experiment,  is set to 0.6. For the experiment with 6ve-decision variables and a large solution space,  is set to 0.75. The global optima stopping parameter  is the number of searches attempted by the algorithm. This value is dependent on the solution space of the problem such that a problem with a larger solution space requires a larger value for this parameter. Tremendous Pexibility in deciding the thoroughness

354

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

of the optimization process is achieved through the use of . For example, if the user is interested in a quick solution that requires few simulation runs, the global optima stopping parameter can be set to a small value. For the two-decision variable, small solution space experiment,  is set to 2, since there are only 300 di8erent combinations of solutions. For the two-decision variable, large solution space experiment and the 6ve-decision variable experiment small solution space experiment, which consisted of 10,000 and 37,500 di8erent possible solutions, respectively,  is set to 3. And for the 6ve-decision variable, large solution space experiment, which has over a million di8erent possible solutions,  is set to 5. 5.4. Setting parameters (simulated annealing) This section describes the parameter selection for the simulated annealing algorithm to be used in comparison to the RS team algorithm. The value of T0 is chosen using Eq. (7) as described in Section 5.3. For the experiment containing two-decision variables and a small solution space, T0 is set such that an inferior solution of 5% would be accepted with a probability of 0.9 on the 6rst iteration, resulting in T0 = 50. For the other three SA experiments, which consisted of a larger overall solution space, a T0 =500 was selected so that an inferior solution of 50% would be accepted with a probability of 0.9 on the 6rst iteration, ensuring a longer and more thorough search of the feasible region. These temperature settings were based on recommendations from Parthasarathy and Rajendran [12]. A geometric cooling schedule is used for the simulated annealing algorithm. For the experiment with 6ve-decision variables and a large solution space, C = 99 and 0.95 for the other three experiments. To de6ne the termination criteria for each experiment,  is set in accordance to the values of T0 and C. For the experiment with a small solution space and two-decision variables  = 4. This allows for a minimum of 50 temperature changes, which is a good amount for a system with a small solution space. For the experiment with two-decision variables and a large solution space and the experiment with 6ve decision variables and a small solution space,  = 3. This allows for 100 temperature changes, which is a thorough amount for systems with solution spaces under 100,000 [17]. For the experiment with 6ve-decision variables and a large solution space  = 9, resulting in a minimum of 400 temperature changes, which is adequate for a solution space of that size [17].

6. Results It should 6rst be noted that the values obtained for the cost per part are not the actual true cost per part in the real life system due to proprietary restrictions. However, these cost per part values still allow for relative comparison of the performance of the two tested methods. Results from each of the four experiments are presented in Fig. 3. This 6gure displays the best cost per part solution found by the RS team method and SA, the alternate solutions found by the RS team method and the number of simulation runs that it took for each method to converge. Alternate solutions are possible for the RS team method because it attempts multiple searches. It should also be noted that full factorial designs were implemented in the two-decision variable experiments and 25−2 factorial designs were used for the 6ve-decision variable experiments.

Design Point I

SA

Design Point 2

SA

Design Point 3

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

SA

RS TEAM

RS TEAM

RS TEAM

Design Point 4

SA

RS TEAM

# of simulation runs 71

Best Solution Found $4633.23, X=(1240,6)

# of simulation runs 117

Best Solution Found 1) $4625.03, X=(1240,6)

# of simulation runs 140

Best Solution Found $4633.23, X=(1240,6)

# of simulation runs 150

Best Solution Found 1) $4633.23, X=(1240,6)

# of simulation runs 277

Best Solution Found $6866.32, X=(1240,5,2,2,2)

# of simulation runs 231

Best Solution Found 1) $6679.70, X=(1040,5,2,2,1)

# of simulation runs 532

Best Solution Found $5926.48, X=(3380,3,2,3,2)

# of simulation runs 454

Best Solution Found 1) $5545.55, X=(2230,3,2,3,2)

355

Other Solutions Found 2) $4633.23, X = (1240, 6)

Other Solutions Found 2) $4689.16, X=(1110, 6) 3) $5363.27, X=(3960,4)

Other Solutions Found 2) $6840.36, X=(1430,5,2,3,2) 3) $6845.99, X=(1350,5,2,3,2)

Other Solutions Found 2) $5952.83, X=(140,3,2,3,1) 3) $6427.49, X=(3220,4,2,3,2) 4) $6892.16, X=(4460,3,2,3,2) 5) $6920.58, X=(3470,3,2,2,2)

Fig. 3. Cost per part results for each design point.

The two-decision variable, small solution space experiment contains the smallest overall solution space of all four design points, 300 di8erent possible solutions. The RS team method achieves the best solution of the two methods, with a cost per part of $4625.03 and the decision variables equal to (1010; 6). SA is slightly inferior obtaining a cost per part of $4633.23 with the decision variables equal to (1240; 6). The RS team method, however, requires 117 total simulation runs to obtain the optimal solution. The RS team method’s performance on this experiment can therefore be described as precise, but ine4cient. SA is able to converge in 71 simulation runs. The RS team method attempts two searches within the feasible region and the other search attempt results in the same solution found by SA. However, an advantage of the RS team method is that it provides two solutions that are equivalent or better than the solutions provided by SA. This gives the user Pexibility in his or her decision making process. There are situations where one set of parameter settings that result in a slightly poorer primary objective are preferable over another. For example, the best solution could result in a high variance in the system response or it could provide other inconveniences, such as a drastic reengineering of the current process. Also, the randomly generated initial solution for SA is (1220; 7) and that is located coincidentally close to the 6nal solution of (1240; 6). This contributes to SA’s ability to converge in such few simulation runs, because it needed only three improving moves to locate the optimal solution. If SA started at a location elsewhere in the feasible region, it is likely that convergence would have taken signi6cantly longer.

356

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

The two-decision variable, large solution space experiment consists of a medium overall solution space of 10,000 possible permutations, but with only two-decision variables. This time the RS team method and SA found equivalent solutions with a cost per part of $4633.23 and decision variable values of (1240; 6), shown in Fig. 3. SA converged in 140 simulation runs, while the RS team method converged in 150 runs. The RS team method is able to 6nd an equivalent solution to SA and two other reasonable solutions, $4689.16 (1110; 6) and $5363.27 (3960; 4). The latter solution contains an inter-arrival rate of 3960, which provides the user with a nice alternative of a large inter-arrival rate even though the cost per part is higher. The 6ve-decision variable, small solution space experiment increased the overall solution space of the system to 37,500. The increase in the amount of decision variables has a signi6cant e8ect on the convergence rate of both algorithms, as reported in Fig. 3. The convergence rates for the RS team method and SA are 231 and 271 simulation runs, respectively, which is an approximate 54% and 98% increase, respectively, from the two-decision variable experiments. The RS team method works much more e4ciently for this experimental system, converging 15% faster than SA while obtaining a better solution. The solution for the RS team method is $6679.90 (1040; 5; 2; 2; 1) and the solution for SA is $6866.32 (1240; 5; 2; 3; 2). The two other solutions that the RS team method produces are also better than the SA solution, $6840.36 (1430; 5; 2; 3; 2) and $6845.99 (1350; 5; 2; 3; 2). Thus, the RS team method performs superior to SA in this experiment, obtaining three solutions in di8erent areas of the feasible region that are superior to the SA solution and it did so in 40 fewer simulation runs. The 6ve-decision variable, large solution space experiment increased the number of permutations of the system from 37,500 to 1,125,000. Since there are still 6ve-decision variables in the system, the results from this experiment explicitly show how an increase in solution space a8ects the performance of each algorithm. As shown in Fig. 2, the RS team method converges to a solution of $5545.55 (2330; 3; 2; 3; 2), in 454 runs, which was almost $400 better than SA’s solution, $5926.48 (3380; 3; 2; 3; 2), found with 532 simulation runs. Again, the RS team method provides a superior 6nal solution, while using 78 less simulation runs than SA. Also, the RS team method obtains four additional solutions to the system in fewer combined simulation runs than SA used to obtain one solution. These results from the last two experiments suggest that as the number of permutations of the system increases to a large value (over 37,500) the RS team method’s performance appears superior to SA’s. Therefore, this proposed method is one to consider when implementing a SA algorithm for optimization of a system with a large, discrete solution space. 7. Conclusions A discrete variable simulation optimization method is proposed in this paper, which incorporates statistical knowledge of the response surface into a traditional SA algorithm. This method performs a more thorough and e4cient search of the feasible region than SA by directly searching out multiple high quality local optimal solutions instead of surveying the entire feasible region for the global optimal solution. The method obtains equivalent or superior solutions to SA with improvements ranging from 1% in the two-decision variable, small solution space experiment to 7% in the 6ve-decision

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358

357

variable, large solution space experiment. Additionally, in the 6ve-decision variable experiments it requires 14.6 –17% fewer simulation runs. A potentially signi6cant advantage of the RS team method within the experiment is that it provides the user with alternative solutions that are in di8erent areas of the feasible region. In the 6ve-decision variable, small solution space experiment, both of the alternate solutions obtained by the RS team method are approximately 1% better than the solution obtained by SA. Alternative solutions can be very bene6cial for the user. For example, if it is desirable to minimize the variance involved in a system while optimizing the response of a system, a slightly higher cost solution with a smaller variance located in a di8erent area of the feasible region could be preferable. If the user is provided multiple solutions with the RS team method, he or she has the Pexibility to choose a solution based on its expected value as well as its variance. Other SA simulation optimization techniques do not provide the user with this Pexibility. Another advantage of the RS team method is that the user has the ability to control the convergence rate of the algorithm by varying the number of Phase I/Phase II searches attempted. Considering the 6ve-decision variable, large solution space experiment performed on the RS team method, if only one search is made by the user, convergence would have occurred in roughly 85 runs. SA converged in 532 runs, so that would have been a substantial savings in simulation runs. This provides great Pexibility in the algorithms’ robustness. The more searches that are attempted, the higher the probability that the user has in obtaining the true optimal solution of the system. However, the more searches that are made, the greater number of simulation runs that are needed. This provides the user with a nice alternative on how to use the algorithm. It can either be used as a quick method to 6nd a local optima solution through one search attempt or a thorough method to obtain a high probability at locating the global optimal solution. The results of this experiment show that the di8erence in 6nal solution quality and convergence rates between the RS team method and SA increased as the size of the solution space increased. It is believed that this trend would continue as the solution space increases even more, making it a favorable choice for a discrete parameter simulation optimization problem. Finally, from this experiment it is concluded that this algorithm should be of future consideration when examining SA based simulation optimization techniques for discrete parameter systems. Acknowledgements This research was supported by the Applied Research Laboratory at the Pennsylvania State University. The authors wish to thank Mark Traband at the Applied Research Laboratory for his support and advice. References [1] Goldsman D, Nelson B. Ranking, selection and multiple comparisons in computer simulation. In: Sadowski D, Seila A, Tew J, Manivannan S, editors. Proceedings of the 1994 Winter Simulation Conference. Piscataway, NJ: Institute of Electrical and Electronics Engineers; 1994. p. 192–9. [2] Kirkpatrick S, Gelatt Jr CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220:671–80. [3] Bulgak A, Sanders J. Integrating a modi6ed simulated annealing algorithm with the simulation of a manufacturing system to optimize bu8er sizes in automatic assembly systems. In: Haigh P, Comfort J, Abrams M, editors.

358

[4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

[14] [15] [16] [17]

S.L. Rosen, C.M. Harmonosky / Computers & Operations Research 32 (2005) 343 – 358 Proceedings of the 1988 Winter Simulation Conference. Piscataway, NJ: Institute of Electrical and Electronics Engineers; 1988. p. 684–90. Haddock J, Mittenthal J. Simulation optimization using simulated annealing. Computers and Industrial Engineering 1992;22:387–95. So D, Dowsland K. Simulated annealing: an application to simulation optimization. Proceedings of the OR35 Conference, 1993. Ahmed M, Alkhamis T, Hasan T. Optimizing discrete stochastic systems using simulated annealing and simulation. Computers and Industrial Engineering 1997;32:823–46. Alkhamis T, Ahmed M, Tuan V. Simulated annealing for discrete optimization with estimation. European Journal of Operational Research 1999;116:530–44. Alrefaei M, Andradottir S. A simulated annealing algorithm with constant temperature for discrete stochastic optimization. Management Science 1999;45:748–64. Ahmed M, Alkhamis T. Simulation-based optimization using simulated annealing with ranking and selection. Computers and Operations Research 2002;29:387–402. Box G, Wilson K. On the experimental attainment of optimal conditions. Journal of the Royal Statistical Society 1951;B13:1–45. Hood S, Welch D. Response surface methodology and its application in simulation. In: Russel E, Biles W, Evans G, Mollaghasemi M, editors. Proceedings of the 1993 Winter Simulation Conference. Piscataway, NJ: Institute of Electrical and Electronics Engineers; 1993. p. 115–22. Parthasarathy S, Rajendran C. A simulated annealing heuristic for scheduling to minimize mean weighted tardiness in a Powshop with sequence-dependent setup times of jobs a case study. Production Planning and Control 1997;8: 475–83. Harmonosky C, Miller J, Rosen S, Traband M, Tillotson R, Robbie D. Interfacing simulation with costing software to drive the transformation from prototype manufacturing to high volume manufacturing. In: Sturrock D, Evans G, Farrington P, Nemhard H, editors. Proceedings of the 1999 Winter Simulation Conference. Piscataway, NJ: Institute of Electrical and Electronics Engineers; 1999. p. 1359–64. Hajek B. Cooling schedules for optimal annealing. Mathematics of Operations Research 1988;13:311–29. Sasaki G, Hajek B. The time complexity of maximum matching by simulated annealing. Journal of the Association of Computing Machinery 1988;35:387–403. Johnson D, Aragon C, McGeoch L, Schevon C. Optimization by simulated annealing: an experimental evaluation; Part I, Graph partitioning. Operations Research 1989;6:865–92. Rosen, S. A new simulation optimization approach for discrete parameter stochastic systems. MS Thesis, The Pennsylvania State University, USA, 2000.

Scott L. Rosen is a Ph.D. candidate in the Harold and Inge Marcus Department of Industrial and Manufacturing Engineering at The Pennsylvania State University. His research interests include simulation output analysis/optimization and decision making under uncertainty. Scott Rosen received a B.S. in Industrial Engineering in 1998 from Lehigh University and an M.S. in both Industrial Engineering and Operations Research in 2000 from Penn State University. He has worked on funded projects with Penn State’s Applied Research Lab and is the recipient of Penn State’s Graduate Scholar Fellowship and Graduate Dean’s Fellowship. Catherine M. Harmonosky is an Associate Professor in the Harold and Inge Marcus Department of Industrial and Manufacturing Engineering at The Pennsylvania State University. Dr. Harmonosky’s research interests are manufacturing and production systems analysis, scheduling and simulation. Current and completed funded projects have been in the areas of scheduling/production control, simulation of manufacturing systems and real-time simulation applications. Dr. Harmonosky received her B.S.I.E. from Penn State in 1981 and her M.S.I.E. and Ph.D. from Purdue University in 1982 and 1987, respectively.