Experimental Design and Analysis

 

Part III Simulation Optimization (NATCORPartIII2015)

10.       Introduction

One of the primary goals of many simulation models is to find the system settings that result in optimal performance. For example

-         Deciding between several different set-ups for a factory or hospital to maximize throughput

-         Choosing the intervention or mix of interventions that leads to the least deaths when treating an infectious disease

-         Optimizing the number of call centre staff on duty to minimize costs subject to constraints on response times

 

Using mathematical notation, we wish to minimize an output f(x), where x is a vector of decision variables and f(x) is the expected value of the random output Y(x),

 

f(x) = Ex[Y(x)].                                                                                                                      (10.1)

 

Nelson and Hong (2009) provide a useful classification for simulation optimization problems into three main groups:

 

1.      The feasible region for x has a small number of solutions; sufficiently small that it is possible to test all of them in the simulation model and select the solution with the best result.

2.      The vector of decision variables, x, is continuous.

3.      The vector of decision variables, x, is discrete and integer ordered.

 

The three examples given above come from each of these classes and are ordered in the same way.

 

In this course, we introduce ranking and selection algorithms that can be used to solve simulation optimization problems of type 1; stochastic approximation algorithms for the solution of problems of type 2 and random search algorithms for the solution of problems of type 3. For more detailed descriptions of all of these methods we would recommend the Handbook of Simulation Optimization edited by Michael Fu.

 

11. Ranking and Selection Algorithms

 

Ranking and selection algorithms are used when the number of options available is small, i.e. x can take only a few values, and it is possible to sample at each of these values. The difficulty lies in the fact that the output of the simulation model is stochastic and that we have only a limited time or computational budget for experiments.

 

We assume that each of the design points x1, x 2, …, x m has a true output associated with it m1, m 2, …, m m. When we run the simulation model, we make ni observations at each of the x i and obtain a set of sample averages Ex[Y(x1)], Ex[Y(x2)], …, Ex[Y(xm)]. These are approximations to the m1, m 2, …, m m.

 

One of the best-known classes of ranking and selection algorithms are Indifference Zone procedures. Indifference zone methods aim to determine the number of observations that need to be made at each design point ni, i = 1, …, m, to ensure that the design point with the smallest output Ex[Y(xk)] has a true output m k, that is within d of the true optimal value m0, such that m0 = min{ mi i = 1, 2, …, m}, with probability greater than 1 – a. The parameter a is the significance level, e.g. 0.05. The variable d is known as the indifference-zone parameter and it is usually set to be the smallest difference that would be practically significant. Writing this more formally, we set the ni such that

 

P{ Ex[Y(xk)] > Ex[Y(xi)] | m(k) –  m(i) ≥ d} ≥ 1 - a.                                                    (11.1)

 

The actual procedures used to determine the ni can be relatively complex and we do not go into details here.

 

 

12. Stochastic Approximation Algorithms

 

Stochastic approximation algorithms are used when the user inputs x, are continuous. They can be thought of as the stochastic equivalent of steepest-descent algorithms, in that they use the gradient of the output f(x) to dictate the direction of movement of the algorithm. 

 

Starting with an initial point x0, the algorithm develops a sequence {xn}, which converges to a unique (local) optimum, i.e. as n tends to infinity, xn will tend to the value of the user input that minimizes f(x).

 

There are two common algorithms on which many of the recent work in this area are based: Robbins-Munro and Kiefer-Wolfowitz. Both have similar iterations and we only include the Kiefer-Wolfowitz equation below,

 

xn+1 = xn – an(Y(xn + cn) – Y(xn – cn))/2cn,

 

where Y(x) is obtained by running the simulation model with input parameters x and consequently is a stochastic variable. The rate of convergence of the algorithm to the optimum is dependent on the choice of the sequences {an} and {cn}, with the optimal choice being somewhat problem-dependent.

 

13. Random Search

 

Random search methods are iterative algorithms that move around the feasible space for the decision variable x in a random fashion. They are useful when the decision variables x are not continuous variables and the number of feasible students is greater than around 20.

 

At each iteration of the algorithm, one or more points will be sampled from the neighbourhood of the current point, and if a sufficient improvement is observed, the algorithm will move to this new point in the next iteration. We introduce a particular example of a random search method termed the Stochastic Ruler method, using the modifications introduced by Alrefaei and Andradσttir (2001).

 

The principle of the Stochastic Ruler method is that in each iteration we compare the value generated by the simulation model with a random value generated from the range of possible values for the objective function. This differentiates it from other random search algorithms in which the comparison is made between realisations of the objective function at different values of the decision variables.

 

Let E[f(x)] be the expected value of the objective function evaluated at x, where we write it as an expectation as it is the average over p iterations. We assume that we know the range of E[f(x)] and that it takes a minimum value of a and a maximum value of b, i.e. the optimal value will lie somewhere on the ruler between a and b. We assume that x can only take values in the space defined by Q and we set up a starting point for the algorithm x0 Ξ Q. 

 

At each iteration of the algorithm, we choose a new point x’ from the neighbourhood of the current point xn. The neighbourhood, N must be defined such that it is equally likely to select x’ when at xn as to select xn when at x’. Examples of neighbourhoods could be any of the other possible values of x or perhaps the values of x that are one different from the current value. Unlike for other random search methods, the choice of neighbourhood structure is not critical to the convergence properties of the algorithm.

 

Algorithm

0.      Define a, b, x0 and the neighbourhood structure N

1.      Generate a candidate solution x’ from the neighbourhood of the current value N(xn, .)

2.      For i = 1 to Mn

Run the simulation model p times to obtain E[f(x)]

Generate a random number U from the uniform distribution U(a, b)

If E[f(x’)] > U then

 xn+1 = xn (i.e. stick with the current solution)

 Exit for

End if

Next i

xn+1 = x’ (i.e. we like this solution, move to this point)

3.      n = n + 1, Go to 1

 

Results suggest that using 1 ≤ M ≤ 5 and p = 1 results in a reasonable performance. The exact choice of M will depend on the neighbourhood structure chosen. If the neighbourhood N is chosen such that few local solutions are generated, a large M is recommended, and vice versa.