A sequential-design metamodeling strategy for simulation optimization

Abstract

This paper examines the feasibility and worthiness of a sequential simulation optimization strategy using nonparametric metamodeling. Incorporating notions advanced in the nonparametric statistics literature, the procedure starts with a uniform grid of points, and then adds points based on the solution of a mathematical programming problem involving quantiles of the squared second derivative of a thin-plate spline metamodel. Termination is reached based on two user-specified criteria.

A feasibility study is conducted, generating 21,000 nonparametric metamodels to fit seven different, basic simulation surfaces. The appropriateness of metamodel fits is judged using recently published criteria. It is concluded that the nonparametric thin-plate spline sequential procedure faithfully reproduces the test case response surfaces and terminates reasonably. However, it is also seen that misleading results may be obtained in systems heavily constrained by budget, and that splines may do a poor job fitting plateaus due to their inherent predisposition to “create ripples.”

Introduction

Consider a continuous process industrial plant for which management develops a software simulation model whose purpose is to determine operating conditions that will maximize plant profits. Certain parameters of the plant's operation are not variable—they may be built into the structure of the plant (for example, the number of various types of machines). These parameters represent uncontrollable inputs in the simulation model. Other parameters may be manipulated, such as valve and pressure settings, and represent controllable inputs in the simulation model. The simulation model may be run for a specified time period at chosen values of the controllable inputs, and the resulting profit calculated. This profit would be the response variable or output from the simulation model, which in a realistic simulation model is usually a random variable—random due to stochastic processing times, queues due to machine breakdowns, etc. Due to its stochastic nature, repeated runs of the model even at the same values of controllable inputs lead to different outputs. Often the average value of the output is calculated, and deviations from the average are analyzed (modeled) subsequently as “error.”

Simulation optimization is the process of determining those values of the controllable input variables that optimize the values of the stochastic output variables. In terms of the industrial scenario mentioned, the simulation optimization problem might be to find the pressure valve settings in the plant that lead to maximal expected profit. The simulation model itself can be thought of as a complex function mapping controllable input values to response values. The optimization is not only stochastic, but is constrained, as the controllable inputs usually have practical limits on their ranges that define a feasible region for exploration. The simulation optimization problem is often characterized as a stochastic search over this feasible exploration region.

This problem is difficult to solve because the function represented by the simulation model is almost always unknown, as is the magnitude and distribution of the error in the response. Further, the size and shape of the error distribution may vary over the ranges of the controllable input variables. Often the problem is made more amenable to solution through the assumption of one or more of the following:

• Domain (“expert”) knowledge can be used to restrict the area of search to a smaller region where the global optimum is known to exist.

• The response function is assumed to be linear or quadratic over portions of the region.

• The error is assumed to be small, homogeneous over the region, and normally distributed.

If all of the above hold, then a well-known technique that uses linear and quadratic regression models such as response surface methodology (RSM) may be successfully employed (see Box and Draper [1], for instance). Other techniques have also been suggested. For a survey of techniques see Meketon [2], Jacobsen and Schruben [3], Safizadeh [4], Azadivar [5] and Fu [6].

In one approach to solving this search problem, a design strategy is employed to determine outputs of the simulation model at carefully chosen input locations within the feasible region. Then a metamodel is fit to these data points and assumed to be a good representation of the simulation model output (which is often termed the true response surface). As such, the metamodel is a surrogate for the simulation model. Decisions such as where to search next or whether the search is over are then often based on the metamodel, rather than the simulation model. The advantage of a metamodel is that it is a simpler representation than the simulation model, and hence is more amenable to ready solution or calculation. The disadvantage is the resultant loss of accuracy; e.g., a metamodel derived from a simulation model of the industrial process described above would be expected to be a less faithful rendition than the simulation model with respect to the plant's “true” profit behavior.

Metamodels can be both parametric, as in the case of RSM, and nonparametric, which are models that do not a priori postulate a functional form. Comparatively little research has been done with nonparametric metamodels for simulation optimization.

Our interest is in approaches that are robust to violations of the above assumptions, that is, it is in those lines of attack that can model response surfaces even with multiple modes (“peaks” or “valleys”) and/or heterogeneous error distributions and variance. Barton [7], [8] provides reviews of such approaches. He notes that nonparametric regression techniques such as kernel and spline smoothing appear to be particularly promising on these tougher problems.

Another factor of interest to us is budget-constrained simulation optimization. Under this scenario, a modeler is limited in the number of data points that can be amassed in search of the optimum, for any of a variety of reasons, including the intrusive nature of the collection process, or its expense.

The optimization of simulation models is an important subject, as several have indicated (see, e.g., Law and Kelton [9], Sargent [10]). For example, as early as 1971 Myers notes the “successful application of known RSM techniques in such areas as chemistry, engineering, biology, agronomy, textiles, the food industry, education, psychology, and others” [11]. Monetary ramifications can be significant; industry generally will not invest in simulation optimization unless the cost implications are sizable. Moreover, the solution of the tougher surfaces can be even more important, as by definition, the peaks and valleys of such systems can be more numerous and difficult to distinguish, and a satisficing or optimal operating region thus more difficult to ascertain. Thus a greedy simulation optimization procedure robust enough to solve tough surfaces is worthy from both a practical as well as a theoretical viewpoint.

Several factors can make metamodeling inappropriate or cumbersome. For example, surface modality (i.e., the number of peaks and valleys in the range of interest), the magnitude of the variance and whether it is homogeneous or heterogeneous over the region, and the distributional form of the error, can render parametric metamodeling ineffective. Moreover, the presence of singularities, discontinuities, and discrete regions in the problem space is problematic.

For the case of fitting nonparametric metamodels, Keys et al. [12] report that error distribution and magnitude are less significant than modality and the basic form of the response function. As a first step, therefore, this research focuses on this prime determinant of metamodeling difficulty for nonparametric approaches, surface modality. We establish a small set of “tough surfaces” for scrutiny, and include surfaces with two, three, and four peaks. Plateaus, ridges and saddle points are included as they are common surface characteristics. Finally, to make sure that this approach is not a step backwards, but that “simple” cases can still be handled, a monomodal (one-peaked) surface is inserted into the study set. Not included in this paper, but left for future work, are surfaces with discontinuities, singularities, and surfaces that are discrete. We note in passing that the user often knows a priori if these are present. Moreover, unlike some metamodeling approaches that assume strict continuity, there is no reason to suspect that nonparametric models will not be able to manage surfaces with these characteristics.

As mentioned, nonparametric metamodels can sometimes function when assumptions prohibit the sanctioned use of parametric approaches. Because no functional form is assumed in advance, and because of the manner in which some nonparametric regressors are calculated, multiple modes, nonnormal error distributions, etc., are often amenable to solution.

Two common nonparametric metamodels are the smoothing spline (see Appendix A) and the kernel smoother (see, e.g., Härdle [13]). Barton [7] mentions these as particularly promising for simulation optimization. We use the thin-plate smoothing spline (the higher dimensional analog of the smoothing spline) in our work. This choice is guided by results obtained by Keys [14], which indicate the properties of the thin plate spline that make it preferable for our work. One property of the spline metamodel surface is that it is smooth and devoid of small features resulting from error terms in the responses. Another property is that the metamodels shift their representation of the response surface when a design point is added; that is, the addition of a data point does not result in a local change to an existing surface, but the generation of a new global solution for the surface. Finally, the thin-plate spline has a well-developed set of theoretical results [13], [15], [16]. Thin plate splines converge asymptotically to the true function at a rate of n^−3/(6+d) where d is the dimension of input space using uniform designs [13], which is a little slower than the parametric regression rate of n^−1/2. However, these results are based on the convergence of the mean squared error (MSE) to zero, and so do not give any guidance as to when the thin-plate spline will represent the true form of the response surface, only that it will converge to it asymptotically. When modeling multimodal surfaces, one cannot assume that the metamodel of the surface faithfully renders the features of the true response surface. A metamodel that does not have the same number of optima in approximately the same locations as the true response surface is termed an alias.

The objective of this research is to develop a sequential, nonparametric algorithm to be used in budget-constrained “tough” cases, and then to test the algorithm as to its ability to determine optima in relatively few computer runs. The study is important in a practical sense in that the benefits of simulation optimization can be extended to scenarios presently not accessible. Of academic import is the realization that the whole field of nonparametric simulation optimization is worthy of increased study.

The rest of paper is organized as follows. The next section contains the sequential algorithm developed for simulation optimization with nonparametric, thin-plate splines. This is followed by an evaluation of the algorithm on seven “tough surfaces,” alluded to above. The paper then sums up, draws conclusions, and discusses future work.