Efficient Simulation Budget Allocation for Subset Selection Using Regression Metamodels

TL;DR

This paper introduces a regression metamodel-based approach for efficient simulation budget allocation in subset selection problems, improving the probability of correctly identifying top designs under budget constraints.

Contribution

It proposes a novel partitioned quadratic regression metamodel and an approximately optimal budget allocation rule for better simulation efficiency in subset selection.

Findings

Significant improvement in simulation efficiency demonstrated through numerical experiments.

Partitioned quadratic regression models effectively fit the domain within each partition.

The proposed allocation rule outperforms existing methods in accuracy and efficiency.

Abstract

This research considers the ranking and selection (R&S) problem of selecting the optimal subset from a finite set of alternative designs. Given the total simulation budget constraint, we aim to maximize the probability of correctly selecting the top-m designs. In order to improve the selection efficiency, we incorporate the information from across the domain into regression metamodels. In this research, we assume that the mean performance of each design is approximately quadratic. To achieve a better fit of this model, we divide the solution space into adjacent partitions such that the quadratic assumption can be satisfied within each partition. Using the large deviation theory, we propose an approximately optimal simulation budget allocation rule in the presence of partitioned domains. Numerical experiments demonstrate that our approach can enhance the simulation efficiency significantly.

Figures

Equations

$$S_o=\left\{i:\underset{i\in S_o}{\max }f(x_i)<\underset{j\in S_n}{\min }f(x_j),i,j=1,...,t\text{and}|S_o|=m\right\},$$

$$f(x_i)={\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2,i=1,...,t.$$

$$\hat{f}(x_i)={\hat{\beta }}_0+{\hat{\beta }}_1{x}_i+{\hat{\beta }}_2{x}_i^2,i=1,...,t.$$

$$\begin{split}&\max_{n_{i}}\ PCS=P\bigg{\{}\bigcap_{i\in S_{o}}\bigcap_{j\in S_{n}}\hat{f}(x_{i})\leq\hat{f}(x_{j})\bigg{\}}\\ &\text{s.t.}\ \sum_{i=1}^{t}n_{i}=n.\end{split}$$

$$\begin{split}\varsigma_{ij}&=\sigma^{2}(0,x_{i}-x_{j},x_{i}^{2}-x_{j}^{2})(\bm{X}^{T}\bm{X})^{-1}\left(\begin{matrix}0\\ x_{i}-x_{j}\\ x_{i}^{2}-x_{j}^{2}\end{matrix}\right)\\ &={\frac{\sigma^{2}}{n}\Big{(}\frac{\rho_{ij,1}^{2}}{\alpha_{1}}+\frac{\rho_{ij,s}^{2}}{\alpha_{s}}+\frac{\rho_{ij,t}^{2}}{\alpha_{t}}\Big{)}},\end{split}$$

$$\{\begin{array}{l}{\rho }_{ij,1}=\frac{(x_s-x_i)(x_t-x_i)-(x_s-x_j)(x_t-x_j)}{(x_1-x_s)(x_1-x_t)},\\ {\rho }_{ij,s}=\frac{(x_1-x_i)(x_t-x_i)-(x_1-x_j)(x_t-x_j)}{(x_s-x_1)(x_s-x_t)},\\ {\rho }_{ij,t}=\frac{(x_1-x_i)(x_s-x_i)-(x_1-x_j)(x_s-x_j)}{(x_t-x_1)(x_t-x_s)}.\end{array}$$

$$\begin{split}&\left.\begin{aligned} -\lim_{n\rightarrow\infty}\frac{1}{n}\log P\{\hat{f}(x_{i})>\hat{f}(x_{m})\},\ &i\in S_{o},i\neq m\\ -\lim_{n\rightarrow\infty}\frac{1}{n}\log P\{\hat{f}(x_{i})<\hat{f}(x_{m})\},\ &i\in S_{n}\end{aligned}\right\}\\ &=R_{m,i}(\bm{\alpha})=\frac{({f}(x_{m})-{f}(x_{i}))^{2}/2}{{\sigma^{2}}\Big{(}\frac{\rho_{mi,1}^{2}}{\alpha_{1}}+\frac{\rho_{mi,s}^{2}}{\alpha_{s}}+\frac{\rho_{mi,t}^{2}}{\alpha_{t}}\Big{)}},\end{split}$$

$$\begin{array}{r}-\underset{n\to \infty }{\lim }\frac{1}{n}\log PFS=\min \left\{\underset{\genfrac{}{}{0ex}{}{i\in S_o}{i\ne m}}{\min }{R}_{m,i}(\mathbf{\alpha }),\underset{j\in S_n}{\min }{R}_{m,j}(\mathbf{\alpha })\right\}.\end{array}$$

$$\begin{array}{rl}\max \min & \left\{\underset{i\in S_o,i\ne m}{\min }{R}_{m,i}(\mathbf{\alpha }),\underset{j\in S_n}{\min }{R}_{m,j}(\mathbf{\alpha })\right\}\\ & \text{s.t.}{\alpha }_1+{\alpha }_s+{\alpha }_t=1\\ & {\alpha }_1,{\alpha }_s,{\alpha }_t\ge 0.\end{array}$$

$$\displaystyle\begin{split}i^{*}=\arg\min\limits_{i=1,...,t}\bigg{\{}\min\limits_{i\in S_{o},i\neq m}R_{m,i}(\bm{\alpha}),\ \min\limits_{i\in S_{n}}R_{m,i}(\bm{\alpha})\bigg{\}}.\end{split}$$

$${\alpha }_r^{\ast }=\{\begin{array}{ll}\frac{|{\rho }_{m{i}^{\ast },r}|}{|{\rho }_{m{i}^{\ast },1}|+|{\rho }_{m{i}^{\ast },s}|+|{\rho }_{m{i}^{\ast },t}|},& r=1,s,t,\\ 0,& \text{otherwise}.\end{array}$$

$$x_s=\{\begin{array}{ll}{x}_{i^{\ast }}+x_m-x_1,& \frac{3x_1+x_t}{4}\le \frac{x_{i^{\ast }}+x_m}{2}<\frac{x_1+x_t}{2},\\ x_{i^{\ast }}+x_m-x_t,& \frac{x_1+x_t}{2}<\frac{x_{i^{\ast }}+x_m}{2}\le \frac{x_1+3x_t}{4},\\ (x_1+x_t)/2.& \text{otherwise}.\end{array}$$

$$\displaystyle\begin{split}S_{o}=\Bigg{\{}i_{h}:\max_{i_{h}\in S_{o}}f(x_{i_{h}})<\min_{j_{g}\in S_{n}}f(x_{j_{g}}),i_{h}=1_{h},...,k_{h},\\ j_{g}=1_{g},...,k_{g}\ h,g=1,...,l\ \text{and}\ |S_{o}|=m\Bigg{\}}.\end{split}$$

$$\begin{split}PCS=P\bigg{\{}\bigcap_{i_{h}\in S_{o}}\bigcap_{j_{g}\in S_{n}}\hat{f}(x_{i_{h}})\leq\hat{f}(x_{j_{g}})\bigg{\}}.\end{split}$$

$$\begin{array}{rl}& \max PCS\\ & \text{s.t.}\sum _{h=1}^l\sum _{i_h=1_h}^{k_h}{n}_{i_h}=n.\end{array}$$

$$\begin{split}\varsigma_{{ij}_{h}}&=\sigma_{h}^{2}(0,x_{i_{h}}-x_{j_{h}},x_{i_{h}}^{2}-x_{j_{h}}^{2})(\bm{X}_{h}^{T}\bm{X}_{h})^{-1}\left(\begin{aligned} &0\\ x_{i_{h}}&-x_{j_{h}}\\ x_{i_{h}}^{2}&-x_{j_{h}}^{2}\end{aligned}\right)\\ &={\frac{\sigma_{b}^{2}}{\theta_{h}n}\Big{(}\frac{\rho_{{ij}_{h},1}^{2}}{\alpha_{1_{h}}}+\frac{\rho_{{ij}_{h},s}^{2}}{\alpha_{s_{h}}}+\frac{\rho_{{ij}_{h},k}^{2}}{\alpha_{k_{h}}}\Big{)}},\end{split}$$

$$\{\begin{array}{l}{\rho }_{{ij}_h,1}=\frac{(x_{s_h}-x_{i_h})(x_{k_h}-x_{i_h})-(x_{s_h}-x_{j_h})(x_{k_h}-x_{j_h})}{(x_{1_h}-x_{s_h})(x_{1_h}-x_{k_h})},\\ {\rho }_{{ij}_h,s}=\frac{(x_{1_h}-x_{i_h})(x_{k_h}-x_{i_h})-(x_{1_h}-x_{j_h})(x_{k_h}-x_{j_h})}{(x_{s_h}-x_{1_h})(x_{s_h}-x_{k_h})},\\ {\rho }_{{ij}_h,k}=\frac{(x_{1_h}-x_{i_h})(x_{s_h}-x_{i_h})-(x_{1_h}-x_{j_h})(x_{s_h}-x_{j_h})}{(x_{k_h}-x_{1_h})(x_{k_h}-x_{s_h})},\end{array}$$

$$\begin{split}\varsigma_{i_{h},j_{g}}&=\sigma_{h}^{2}(0,x_{i_{h}}-x_{j_{g}},x_{i_{h}}^{2}-x_{j_{g}}^{2})(\bm{X}_{h}^{T}\bm{X}_{h})^{-1}\left(\begin{aligned} &0\\ x_{i_{h}}&-x_{j_{g}}\\ x_{i_{h}}^{2}&-x_{j_{g}}^{2}\end{aligned}\right)\\ &={\frac{\sigma_{h}^{2}}{\theta_{h}n}\!\Big{(}\!\frac{\eta_{i_{h}\!,1}^{2}}{\alpha_{1_{b}}}\!+\!\frac{\eta_{i_{h}\!,s}^{2}}{\alpha_{s_{b}}}\!+\!\frac{\eta_{i_{h}\!,k}^{2}}{\alpha_{k_{b}}}\!\Big{)}\!+\!\frac{\sigma_{g}^{2}}{\theta_{g}n}\!\Big{(}\!\frac{\eta_{j_{g}\!,1}^{2}}{\alpha_{1_{h}}}\!+\!\frac{\eta_{j_{g}\!,s}^{2}}{\alpha_{s_{h}}}\!+\!\frac{\eta_{j_{g}\!,k}^{2}}{\alpha_{k_{h}}}\!\Big{)}},\end{split}$$

$$\{\begin{array}{l}{\eta }_{i_h,1}=\frac{(x_{s_h}-x_{i_h})(x_{k_h}-x_{i_h})}{(x_{1_h}-x_{s_h})(x_{1_h}-x_{k_h})},\\ {\eta }_{i_h,s}=\frac{(x_{1_h}-x_{i_h})(x_{k_h}-x_{i_h})}{(x_{s_h}-x_{1_h})(x_{s_h}-x_{k_h})},\\ {\eta }_{i_h,k}=\frac{(x_{1_h}-x_{i_h})(x_{s_h}-x_{i_h})}{(x_{k_h}-x_{1_h})(x_{k_h}-x_{s_h})}.\end{array}$$

$$\begin{split}&\left.\begin{aligned} \!-\lim_{n\rightarrow\infty}\!\frac{1}{n}\!\log P\{\hat{f}(x_{i_{h}})\!>\!\hat{f}(x_{m_{b}})\},\ i_{h}\!\in\!S_{o}\!\\ \!-\lim_{n\rightarrow\infty}\!\frac{1}{n}\!\log P\{\hat{f}(x_{i_{h}})\!<\!\hat{f}(x_{m_{b}})\},\ i_{h}\!\in\!S_{n}\!\end{aligned}\right\}\!=\!R_{m_{b},i_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}})\\ &=\frac{({f}(x_{m_{b}})\!-\!{f}(x_{i_{h}}))^{2}/2}{\frac{\sigma_{b}^{2}}{\theta_{b}}\!\Big{(}\!\frac{\eta_{m_{b}\!,1}^{2}}{\alpha_{1_{b}}}\!+\!\frac{\eta_{m_{b}\!,s}^{2}}{\alpha_{s_{b}}}\!+\!\frac{\eta_{m_{b}\!,k}^{2}}{\alpha_{k_{b}}}\!\Big{)}\!+\!\frac{\sigma_{h}^{2}}{\theta_{h}}\!\Big{(}\!\frac{\eta_{i_{h}\!,1}^{2}}{\alpha_{1_{h}}}\!+\!\frac{\eta_{i_{h}\!,s}^{2}}{\alpha_{s_{h}}}\!+\!\frac{\eta_{i_{h}\!,k}^{2}}{\alpha_{k_{h}}}\!\Big{)}}\!,h\!\neq\!b,\end{split}$$

$$\begin{split}&\left.\begin{aligned} \!-\!\lim_{n\rightarrow\infty}\frac{1}{n}\log\!P\{\hat{f}(x_{i_{b}})\!>\!\hat{f}(x_{m_{b}})\},\ &i_{b}\!\in\!S_{o},i_{b}\!\neq\!m_{b}\\ \!-\!\lim_{n\rightarrow\infty}\frac{1}{n}\log\!P\{\hat{f}(x_{i_{b}})\!<\!\hat{f}(x_{m_{b}})\},\ &i_{b}\!\in\!S_{n}\end{aligned}\right\}\!=\!R_{m_{b},i_{b}}(\theta_{b},\bm{\alpha_{b}})\\ &=\frac{({f}(x_{m_{b}})-{f}(x_{i_{b}}))^{2}/2}{\frac{\sigma_{b}^{2}}{\theta_{b}}\Big{(}\frac{\rho_{{mi}_{b},1}^{2}}{\alpha_{1_{b}}}+\frac{\rho_{{mi}_{b},s}^{2}}{\alpha_{s_{b}}}+\frac{\rho_{{mi}_{b},k}^{2}}{\alpha_{k_{b}}}\Big{)}}.\end{split}$$

$$\displaystyle\max\Bigg{\{}\!\!\max_{1\leq h\leq l\atop h\neq b}\!\!\bigg{\{}$$

$$\displaystyle\max\limits_{j_{h}\in S_{n}}\!P\{\hat{f}(x_{j_{h}})\!\leq\!\hat{f}(x_{m_{b}})\}\!\bigg{\}},$$

$$\displaystyle\max\bigg{\{}$$

$$\displaystyle\max\limits_{j_{b}\in S_{n}}P\{\hat{f}(x_{j_{b}})\leq\hat{f}(x_{m_{b}})\}\bigg{\}}\Bigg{\}}$$

$$\displaystyle|S_{o}|\!\times\!|S_{n}|\!\times\!\max\Bigg{\{}\max_{1\leq h\leq l\atop h\neq b}\bigg{\{}$$

$$\displaystyle\max\limits_{j_{h}\in S_{n}}P\{\hat{f}(x_{j_{h}})\leq\hat{f}(x_{m_{b}})\}\bigg{\}},$$

$$\displaystyle\max\bigg{\{}$$

$$\displaystyle\max\limits_{j_{b}\in S_{n}}P\{\hat{f}(x_{j_{b}})\leq\hat{f}(x_{m_{b}})\}\bigg{\}}\Bigg{\}}.$$

$$\begin{split}-\lim_{n\rightarrow\infty}\frac{1}{n}\log PFS=\min\!\Bigg{\{}\!\!\!\min_{1\leq h\leq l\atop h\neq b}\!\!\bigg{\{}\!\!\min\limits_{i_{h}\in S_{o}}\!R_{m_{b},i_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}}),\\ \!\min\limits_{j_{h}\in S_{n}}\!\!R_{m_{b},j_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}})\bigg{\}},\\ \min\bigg{\{}\!\!\min\limits_{i_{b}\in S_{o}\atop i_{b}\neq m_{b}}\!R_{m_{b},i_{b}}(\theta_{b},\bm{\alpha_{b}}),\min\limits_{j_{b}\in S_{n}}R_{m_{b},j_{b}}(\theta_{b},\bm{\alpha_{b}})\bigg{\}}\Bigg{\}},\end{split}$$

$$\displaystyle\begin{split}\max&\min\Bigg{\{}\min_{1\leq h\leq l\atop h\neq b}\bigg{\{}\min\limits_{i_{h}\in S_{o}}R_{m_{b},i_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}}),\\ &\min\limits_{j_{h}\in S_{n}}\!\!R_{m_{b},j_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}})\bigg{\}},\\ &\min\bigg{\{}\min\limits_{i_{b}\in S_{o}\atop i_{b}\neq m_{b}}R_{m_{b},i_{b}}(\theta_{b},\bm{\alpha_{b}}),\min\limits_{j_{b}\in S_{n}}R_{m_{b},j_{b}}(\theta_{b},\bm{\alpha_{b}})\bigg{\}}\Bigg{\}}\\ \text{s.t.}\ \ &\alpha_{1_{h}}+\alpha_{s_{h}}+\alpha_{k_{h}}=1,\ h=1,...,l,\\ &\sum_{h=1}^{l}\theta_{h}=1,\\ &\theta_{h},\alpha_{1_{h}},\alpha_{s_{h}},\alpha_{k_{h}}\geq 0,\ h=1,...,l.\end{split}$$

$$\displaystyle\begin{split}\widetilde{R}_{m_{b},i_{h}}(\theta_{h},\bm{\alpha_{h}})=\lim_{l\rightarrow\infty}R_{m_{b},i_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}})\\ =\frac{(f(x_{m_{b}})-f(x_{i_{h}}))^{2}/2}{\frac{\sigma_{h}^{2}}{\theta_{h}}\Big{(}\frac{\eta_{i_{h},1}^{2}}{\alpha_{1_{h}}}+\frac{\eta_{i_{h},s}^{2}}{\alpha_{s_{h}}}+\frac{\eta_{i_{h},k}^{2}}{\alpha_{k_{h}}}\Big{)}}.\end{split}$$

$$\displaystyle\begin{split}\max\min\Bigg{\{}&\!\!\min_{1\leq h\leq l\atop h\neq b}\!\bigg{\{}\!\min\limits_{i_{h}\in S_{o}}\!\!\widetilde{R}_{m_{b},i_{h}}(\theta_{h},\bm{\alpha_{h}}),\min\limits_{j_{h}\in S_{n}}\!\!\widetilde{R}_{m_{b},j_{h}}(\theta_{h},\bm{\alpha_{h}})\bigg{\}},\\ &\min\bigg{\{}\!\!\min\limits_{i_{b}\in S_{o}\atop i_{b}\neq m_{b}}\!\!R_{m_{b},i_{b}}(\theta_{b},\bm{\alpha_{b}}),\min\limits_{j_{b}\in S_{n}}\!\!R_{m_{b},j_{b}}(\theta_{b},\bm{\alpha_{b}})\bigg{\}}\Bigg{\}}\\ &\text{s.t.}\ \alpha_{1_{h}}+\alpha_{s_{h}}+\alpha_{k_{h}}=1,h=1,...,l\\ &\sum_{h=1}^{l}\theta_{h}=1\\ &\theta_{h},\alpha_{1_{h}},\alpha_{s_{h}},\alpha_{k_{h}}\geq 0,h=1,...,l.\end{split}$$

Full text

, , ,

1 Introduction

Discrete-event systems (DES) simulation has played an important role in analyzing modern complex systems and evaluating decision problems, since these systems are usually too difficult to be described using analytical models. DES simulation has been a common analysis method of choice and widely used in many practical applications, such as the queueing systems, electric power grids, air and land traffic control systems, manufacturing plants and supply chains (Xu et al.,, 2015, 2016; Gao and Chen,, 2016). However, running the simulation model is usually time consuming, and a large number of simulation replications are typically required to achieve an accurate estimate of a design decision (Lee et al.,, 2010). In addition, it could be computationally quite expensive to select the best design(s) when the number of alternatives is relatively large.

In this paper, we consider the problem of selecting the optimal subset of the top-$m$ designs out of $t$ alternatives, where the performance of each design is estimated based on simulation. In order to improve the selection efficiency, we aim to intelligently allocate the simulation replications to each design to maximize the probability of correctly selecting all the top-$m$ designs. This problem setting falls in the well-established statistics branch known as ranking and selection (R&S) (Xu et al.,, 2015).

In the literature, several types of efficient R&S procedures have been developed. The indifference-zone (IZ) approach allocates the simulation budget to provide a guaranteed lower bound for the probability of correct selection ($PCS$) (Kim and Nelson,, 2001). Chen et al., (2000) proposed an optimal computing budget allocation (OCBA) approach for R&S problems. The OCBA approach allocates the simulation replications sequentially in order to maximize $PCS$ under a simulation budget constraint. He et al., (2007), Gao and Shi, (2015) and Gao et al., 2017a further developed the OCBA method with the expected opportunity cost ($EOC$) measure, which focuses more on the consequence of a wrong selection compared to $PCS$. Brantley et al., (2013) proposed another approach called optimal simulation design (OSD) to select the best design with regression metamodels. It assumes that all designs fit a single quadratic line and the variance of each design is identically distributed. The OSD approach was further extended in Brantley et al., (2014), Xiao et al., (2015) and Gao et al., (2018) to consider more general problems by dividing the solution space into adjacent partitions. Although these studies are also based on partitioning or metamodeling, they do not aim to select the top $m$ designs, and are therefore different in objectives from this research. Other variants of OCBA include selecting the best design considering resource sharing and allocation (Peng et al.,, 2013) and input uncertainty (Gao et al., 2017b, ).

Most of the existing R&S procedures focus on identifying the best design and return a single choice as the estimated optimum. However, decision makers may prefer to have several good alternatives instead of one and make the final selection by considering some qualitative criteria, such as political feasibility and environmental consideration, which might be neglected by computer models (Gao and Chen,, 2016; Zhang et al.,, 2016). The selection procedure providing the top-$m$ designs can help the decision makers make their final decision in a more flexible way.

The literature on the optimal subset selection is sparse. Chen et al., (2008) and Zhang et al., (2016) considered the optimal subset selection problem using the OCBA framework, which maximizes $PCS$ under a simulation budget constraint. In Gao and Chen, (2015), an optimal subset selection procedure was proposed to minimize the measure of $EOC$. The optimal subset selection problem was further extended in Gao and Chen, (2016) for general underlying distributions using the large deviations theory.

The aforementioned R&S procedures could smartly allocate the computing budget given the simulation results. They estimate the performance of each design only by considering the sample information of the considered design itself. However, the designs nearby could also provide useful information since neighboring designs usually have similar performance. Based on this idea, we aim to improve the selection efficiency by incorporating the information from across the domain into some response surfaces. Unlike traditional R&S methods, the regression based approaches require simulation experiments on only a subset of all the designs under consideration. The performances of the other designs can be inferred based on the sample information of the simulated designs. This provides us an effective way to further improve the efficiency for solving the subset selection problem, which is the motivation of this paper.

In this research, we assume the underlying function is quadratic or approximately quadratic. This assumption could help utilize the structure information of the design space and led to significant improvement of the computational efficiency. It is commonly used in the literature, such as Brantley et al., (2013, 2014); Xiao et al., (2015); McConnell and Servaes, (1990). Based on this assumption we built a quadratic regression metamodel to incorporate the information from across the domain. The first contribution of this work is that we propose an asymptotically optimal allocation rule that determines which designs need to be simulated and the number of simulation budget allocated to them, such that the $PCS$ of the optimal subset could be maximized. We call this procedure the optimal computing budget allocation for selecting the top-$m$ designs with regression (OCBA-mr). In order to further extend the OCBA-mr procedure to more general cases where the underlying function is partially quadratic or non-quadratic, we divide the solution space into adjacent partitions and build a quadratic regression metamodel within each partition. The underlying function in each partition could be well approximated by a quadratic function if the solution space is properly partitioned or each partition is small enough. According to the results in Brantley et al., (2014); Xiao et al., (2015); Gao et al., (2018), the use of partitioned domains along with regression metamodels could significantly improve the simulation efficiency. That means interpolating the solution space can be an effective way for us to have further improvement. For different problems, the solution space could be divided into discrete partitions using different criteria, such as the size of corporations, the type of industries and the temperature of chemical process (Xiao et al.,, 2015). Based on the idea mentioned above, we develop an asymptotically optimal computing budget allocation procedure for selecting the top-$m$ designs with regression in partitioned domains (OCBA-mrp), which is an extension of the OCBA-mr procedure for more general cases. In order to maximize the $PCS$ of the optimal subset, the OCBA-mrp procedure not only determines the optimal simulation budget allocation within each partition but also determines the optimal budget allocation between partitions.

The rest of the paper is organized as follows. In Section 2, we formulate the optimal subset selection problem with regression metamodel and derive an asymptotically optimal simulation budget allocation rule, called OCBA-mr. Section 3 extends the OCBA-mr method for more general cases with partitioned domains and derives another asymptotically optimal simulation budget allocation rule, called OCBA-mrp. The performance of the proposed methods is illustrated with numerical examples in Section 4. Section 5 concludes the paper.

2 Optimal Subset Selection Strategy

In this section, we provide an optimal computing budget allocation rule for the subset selection problem based on the regression metamodel.

2.1 Problem formulation

Without loss of generality, the best design is defined as the design with the smallest mean performance. We introduce the following notations:

• $t$:          total number of designs;
• $x_{i}$:        location of design $i$;
• ${m}$:        design with the $m$-th smallest mean value;
• $f(x_{i})$:    mean performance value for design ${i}$;
• $Y(x_{i})$:   simulation output for design ${i}$;
• $\bm{\beta}$:        vector of $(\beta_{0},\beta_{1},\beta_{2})$ , representing the coefficients of the regression model;
• $\hat{f}(x_{i})$:    estimate of $f(x_{i})$ based on the regression model;
• $S_{o}$:        set of the true top-$m$ designs;
• $S_{n}$:        set of designs not in $S_{o}$ (complement of $S_{o}$);
• $n$:         total number of simulation replications;
• $n_{i}$:        number of simulation replications allocated to design ${i}$;
• $\bm{\alpha}$:         vector of $(\alpha_{1},\alpha_{s},\alpha_{k})$, where $\alpha_{r}$ is the proportion of the total simulation budget $n$ allocated to design ${r}$, i.e., $\alpha_{r}\!\!=\!n_{r}/n$, $r\!=\!1,s,t$.

In the last bulleted item, $r$ takes values only at $1$, $s$ and $t$, where designs $1$ and $t$ are the first and last designs in the solution space, and design $s$ is an intermediate design determined by (5) (the rationale of it will be explained in more detail in Theorem 2). The problem we considered in this paper is to select the top-$m$ designs out of the $t$ alternatives by allocating simulation replications to these three designs. The optimal set is

[TABLE]

where $|S|$ is the size of subset $S$.

We study the problem where the expected performance value across the solution space is quadratic or approximately quadratic in nature. Then, the mean performance of design ${i}$ can be written as

[TABLE]

The coefficients $\bm{\beta}\!=\!({\beta}_{0},{\beta}_{1},{\beta}_{2})$ are unknown beforehand, but can be estimated via simulation samples. Assume that the noises $(\varepsilon)$ of the simulation experiments of design ${i}$ follow normal distribution $N(0,\sigma^{2})$. For each design, the noise is independent from replication to replication. The simulation output is $Y(x_{i})=f(x_{i})+\varepsilon(x_{i})$ with $\varepsilon(x_{i})\sim N(0,\sigma^{2}).$

Given a total of $n$ samples, we define $\bm{Y}$ as the vector of the $n$ simulation samples $Y(x_{i})$ and $\bm{X}$ be an $n\times 3$ matrix with each row $(1,x_{i},x_{i}^{2})$ corresponding to each entry $Y(x_{i})$ of $\bm{Y}$. We estimate the coefficients $\bm{\beta}$ using the ordinary least squares (OLS) method (Hayashi,, 2000) and denote the estimates as $\hat{\bm{\beta}}\!=\!(\hat{\beta}_{0},\hat{\beta}_{1},\hat{\beta}_{2})$. Then, we have $\bm{\hat{\beta}}=(\bm{X}^{T}\bm{X})^{-1}\bm{X}^{T}\bm{Y}$ and $Cov[\bm{\hat{\beta}}]=\sigma^{2}(\bm{X}^{T}\bm{X})^{-1},$ where $T$ denotes transposition and $\bm{X}^{T}\bm{X}$ is known as the information matrix (Kiefer,, 1959). The estimate of ${f}(x_{i})$ can be written as

[TABLE]

We can use the equation (1) which incorporates the sample information of the simulated designs to estimate the expected performance for each design across the solution space. As $\hat{f}(x_{i})$ is a linear combination of $\bm{\hat{\beta}}$, we have $Var[\hat{f}(x_{i})]\!=\!\sigma^{2}\bm{X}_{i}^{T}\!(\bm{X}^{T}\!\bm{X})^{-1}\!\bm{X}_{i},$ where $\bm{X}_{i}^{T}=(1,x_{i},x_{i}^{2})$.

Due to the uncertainty of the estimate of the underlying function, a correct selection of the optimal subset $S_{o}$ may not always occur. Therefore, we introduce a measure, the probability of correct selection ($PCS$), to formulate the R&S problem considered in this paper. The $PCS$ is given by PCS=P\bigg{\{}\bigcap_{i\in S_{o}}\bigcap_{j\in S_{n}}\hat{f}(x_{i})\leq\hat{f}(x_{j})\bigg{\}}.

Given a fixed simulation budget, the optimization problem can be written as follows:

[TABLE]

In this section, we aim to solve problem (2) where the mean performance of each design is estimated using the regression metamodel. Due to the uncertainty in the simulation experiments, multiple simulation replications are needed to generate the estimates of the underlying function $\hat{f}(x_{i})$ accurately. The variance of $\hat{f}(x_{i})$ is a function of the information matrix $(\bm{X}^{T}\bm{X})$ and can be reduced if additional simulation replications are conducted (Xiao et al.,, 2015). We are interested in how to intelligently allocate the simulation budget to proper designs such that $\hat{f}(x_{i})$ can be better estimated using regression metamodel, and the $PCS$ in problem (2) can be maximized.

2.2 Optimization model under large deviations framework

A major difficulty in solving problem (2) is that the objective function $PCS$ does not have a closed-form expression. We seek to solve this optimization problem under an asymptotic framework in which the $PCS$ is maximized or the probability of false selection $(PFS=1-PCS)$ is minimized as $n$ goes to infinity.

The $PCS$ used in this paper is defined based on the quadratic regression model (1), which is constructed using the simulation information of only a fraction of the $t$ designs. We call the designs receiving simulation replications the support designs. In order to construct the quadratic regression model, we need at least three support designs to obtain all of the information in $\bm{X}^{T}\bm{X}$ (Kiefer,, 1959). For simplicity, in this research we let the number of support designs be three, two of which are at the extreme locations, i.e., $x_{1}$ and $x_{t}$ (for this setting, see, e.g., Brantley et al., (2013, 2014); Xiao et al., (2015); Kiefer, (1959)). The case with more than three designs can be similarly analyzed.

Lemma 1: Let $\hat{d}_{i,j}=\hat{f}(x_{i})-\hat{f}(x_{j})$, $i,j=1,2,...,t,i\neq j$. $\hat{d}_{i,j}$ follows normal distribution N\big{(}{f}(x_{i})-{f}(x_{j}),\varsigma_{ij}\big{)}, where

[TABLE]

and

[TABLE]

The proof is similar to that is given in Eq.(22) of Brantley et al., (2013) (becomes the same when $m=1$), and hence is omitted for brevity.

For any $i,j=1,2,...,t,i\neq j$, $\hat{d}_{i,j}$ is a normally distributed random variable. We can use the large deviations theory to derive the convergence rate function of the false selection probability.

Lemma 2: The convergence rate function of incorrect comparison probability for each design is:

[TABLE]

Based on the results in Lemma 2, we can get an explicit expression of the convergence rate function of $PFS$.

Lemma 3: The convergence rate function of $PFS$ is:

[TABLE]

The main assertion of Lemma 3 is that the overall convergence rate of $PFS$ is determined by the minimum convergence rate of the incorrect comparison for each design. Minimizing the $PFS$ is asymptotically equivalent to maximizing the rate at which $PFS$ goes to zero as a function of $\bm{\alpha}$, i.e., maximizing $-\lim\limits_{n\rightarrow\infty}\frac{1}{n}\log PFS$. Based on Lemma 3, the asymptotical version of (2) becomes

[TABLE]

2.3 Asymptotically optimal solution

In this subsection, we seek to derive the optimality conditions for (3). Since the overall convergence rate of $PFS$ is determined by the design with the minimum convergence rate. A false selection is most likely to happen at this key design. Therefore, it is enough for us to investigate the convergence rate of the key design across the solution space. We define $i^{*}$ as the key design. That is

[TABLE]

Theorem 1**.**

The optimization problem (2) can be asymptotically optimized with the following allocation rule:

[TABLE]

The $\rho_{mi^{*},r}$ is also known as the Lagrange interpolating polynomial coefficient (De la Garza et al.,, 1954; Burden and Faires,, 2001). It represents the relative importance of each support design for estimating $\hat{d}_{mi}=\hat{f}(x_{m})-\hat{f}(x_{i})$. Theorem 1 indicates that the support design $r$ will receive more simulation budget if it has larger $\rho_{mi^{*},r}$.

Given the optimal allocation rule (Theorem 1), we next determine the optimal location for the support design $s$.

Theorem 2**.**

The rate function of $PFS$ with allocation rule satisfying (4) can be maximized if the support design $s$ satisfies the following equations.

The support design

[TABLE]

When the $x_{s}$ derived from (5) does not correspond to any design available, we round it to the nearest one. The expression of Theorem 2 is similar to the results in Brantley et al., (2013). The differences are the selection of the key design $x_{i}^{*}$. We define $x_{i}^{*}$ as the design with the minimum convergence rate of PFS, where a false selection of the optimal subset is most likely to happen.

3 Optimal Subset Selection Strategy for Partitioned Domains

The regression based method mentioned above can greatly improve the subset selection efficiency compared to the traditional methods. However, it is constrained with the typical assumptions such as the assumption of quadratic underlying function for the means. It is possible that the underlying function is neither quadratic nor approximately quadratic, so that we will fail to find the top-$m$ designs. In order to extend our method to more general cases, we divide the solution space into adjacent partitions. The quadratic pattern can be expected when the solution space is properly partitioned or each partition is small enough.

3.1 Problem formulation

We first add the following notations for partitioned domains:

• $l$:        number of partitions of the entire domain;
• $k^{h}$:      number of designs in partition $h$, $h=1,2,...,l$;
• ${i_{h}}$:       $i$th design in partition $h$, $i=1,2,...,k^{h}$ (when $i=k^{h}$, we denote design $k^{h}_{h}$​ as design $k_{h}$ for notational simplicity);
• $x_{i_{h}}$:     location of design $i_{h}$;
• $b$:        the partition containing the design with the $m$-th smallest mean value;
• ${m_{b}}$:     design with the $m$-th smallest mean value;
• $\bm{\beta_{h}}$:      vector of $(\beta_{h0},\beta_{h1},\beta_{h2})$ , representing the coefficients of the regression model in partition $h$;
• $n_{h\bullet}$:     number of simulation replications allocated to partition $h$;
• $n_{i_{h}}$:     number of simulation replications allocated to design ${i_{h}}$;
• $\theta_{h}$:      the proportion of simulation budget allocated to partition $h$, i.e., $\theta_{h}=n_{h\bullet}/n$;
• $\bm{\alpha_{h}}$:     vector of $(\alpha_{1_{h}},\alpha_{s_{h}},\alpha_{k_{h}})$, where $\alpha_{r_{h}}$ is the proportion of the simulation budget $n_{h\bullet}$ allocated to design ${r_{h}}$, i.e., $\alpha_{r_{h}}\!\!=\!n_{r_{h}}/n_{h\bullet}$, $r_{h}\!=\!1_{h},s_{h},k_{h}$.

The notations $f(x_{i_{h}})$, $Y(x_{i_{h}})$ and $\hat{f}(x_{i_{h}})$ defined for partitioned domains are similar to those in Section 2, except the design $i$ is replaced by $i_{h}$. The entire domain is divided into $l$ adjacent partitions. Each partition contains $k^{h}$ designs, i.e., there are $\sum_{h=1}^{l}k^{h}=t$ designs in total.

We assume there exists $\epsilon\in\mathbb{R}$ such that there are exactly $m$ designs from the total $t$ designs with mean performances less than $\epsilon$ and the rest $t-m$ designs with mean performances greater than $\epsilon$. It ensures that the optimal subset is well defined and can be distinguished. We define the optimal subset $S_{o}$ as

[TABLE]

In this section, we assume that the expected performance value in each partition is quadratic or approximately quadratic when the solution space is properly partitioned and the mean performance within a partition is continuous and smooth. For this problem setting the $PCS$ is defined as

[TABLE]

Given a fixed simulation budget, the optimization problem can be written as follows:

[TABLE]

In this section, we aim to solve problem (6) in the presence of regression metamodels. We are interested in how to intelligently allocate the simulation budget to proper designs such that $\hat{f}(x_{i_{h}})$ can be better estimated using regression metamodels, and the $PCS$ in problem (6) can be maximized. Note that this model does not require all the designs to be on the same axis, and therefore does not hinder it from being applied for multi-dimensional problems. For multi-dimensional problems, we can treat the range of the underlying function on each dimension as one or more partitions, and then apply this formulation.

3.2 Optimization model under large deviations framework

In order to solve the optimization problem (6), one challenge is how to derive an explicit expression of the $PCS$. We seek to solve (6) under an asymptotic framework in which the probability of false selection ($PFS=1-PCS$) is minimized as $n$ goes to infinity.

Similar to the setting in Section 2, we let the number of support designs in each partition be three, two of which are at the extreme locations, i.e., $x_{1_{h}}$ and $x_{k_{h}}$. We have $\hat{d}_{{ij}_{h}}=\hat{f}(x_{i_{h}})-\hat{f}(x_{j_{h}})$ follow normal distribution N\big{(}{f}(x_{i_{h}})-{f}(x_{m_{h}}),\varsigma_{{ij}_{h}}\big{)} and $\hat{d}_{{i_{h}},j_{g}}=\hat{f}(x_{i_{h}})-\hat{f}(x_{j_{g}})$ follow normal distribution N\big{(}{f}(x_{i_{h}})-{f}(x_{j_{g}}),\varsigma_{i_{h},j_{g}}\big{)} where $h\neq g$.

Similar to the proof of Lemma 1, we have

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

For any $h=1,2,...,l$ and $i_{h}=1_{h},2_{h},...,k_{h}$, $\hat{f}(x_{i_{h}})$ is a normally distributed random variable. We can use the large deviations theory to derive the convergence rate function of the false selection probability.

Lemma 4: The convergence rate functions of incorrect comparison probability for each design are provided as follows:

[TABLE]

and

[TABLE]

According to the Bonferroni inequality, we have $PFS\!\leq\!\!\!\sum_{i_{h}\in S_{o},\atop i_{h}\neq m_{b}}\!\!\!P\{\hat{f}(x_{i_{h}})\!\geq\!\hat{f}(x_{m_{b}})\}+\!\!\!\sum_{j_{g}\in S_{n}}\!\!P\{\hat{f}(x_{j_{g}})\!\leq\!\hat{f}(x_{m_{b}})\}.$ The $PFS$ is bounded below by

[TABLE]

and bounded above by

[TABLE]

Therefore, according to the results in Lemma 4, the convergence rate functions of $PFS$ is given by

[TABLE]

Minimizing the $PFS$ is asymptotically equivalent to maximizing the rate at which $PFS$ goes to zero as a function of $\theta_{h}$ and $\bm{\alpha_{h}}$, $h=1,2,...,l.$ Similarly, the asymptotical version of (6) becomes

[TABLE]

3.3 Asymptotically optimal solution

In this section, we seek to derive the optimality conditions for (9). We want to determine (i) the number of simulation replication allocated to each partition, (ii) the locations of the designs should be simulated in each partition and (iii) the number of simulation replications allocated to those selected designs.

According to Xiao et al., (2015), we have $\theta_{h}/\theta_{b}\rightarrow 0$ as the number of partitions $l$ goes to infinity. That means the fraction of the simulation budget allocated to the partition $b$ containing the $m$-th smallest mean value far exceeds the fraction given to any other partition when $l$ goes to infinity. Given that, the rate function $R_{m_{b},i_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}},\bm{\alpha_{h}})$ converges to $\widetilde{R}_{m_{b},i_{h}}(\theta_{h},\bm{\alpha_{h}})$ where

[TABLE]

The problem (9) can be asymptotically rewritten as

[TABLE]

In order to better analyze problem (10) above, we decompose it as follows:

[TABLE]

for $h=1,2,...,l,h\neq b$, and

[TABLE]

for $h=b$.

Lemma 5: Let $\theta_{h}^{*}$ and $\bm{\alpha_{h}^{*}}$ be the optimal solution to (10). The $\theta_{h}^{*}$ and $\bm{\alpha_{h}^{*}}$ are independent and can be solved separately. In addition, the $l$ $\bm{\alpha_{h}^{*}}$’s $(h=1,2,...,l)$ corresponding to the $l$ partitions are also mutually independent and can be solved separately using (11) and (12).

Similar to Section 2, we define a key design of each partition $h$, denoted as $i_{h}^{*}$. A false selection is most likely to happen at these key designs. That is

[TABLE]

3.3.1 Determine $\bm{\alpha_{h}^{*}}$ and locations of support designs

Given Lemma 5, we can determine $\bm{\alpha_{h}^{*}}$ separately for each partition by solving the optimization problems (11) and (12).

Theorem 3**.**

The optimization problem (6) can be asymptotically optimized with the following allocation rule:

[TABLE]

where

[TABLE]

for $r_{h}=1_{h},s_{h},k_{h}$.

Given the optimal allocation rule (Theorem 3), we next determine the optimal location for the support design $s_{h}$ within each partition.

Theorem 4**.**

The rate function of $PFS$ with allocation rule satisfying (13) can be asymptotically maximized if the support design $s_{h}$ satisfies the following equations.

(a) For all $h\neq b$, support design

[TABLE]

and the optimal allocation rule (13) becomes

[TABLE]

(b) For $h=b$, support design

[TABLE]

and the approximate optimal allocation rule (13) is

[TABLE]

When the $x_{s_{h}}$ derived from (14) and (16) does not correspond to any design available, we round it to the nearest one.

3.3.2 Determine $\theta_{h}^{*}$

In this subsection, we aim to determine the optimal budget allocation rules between the $l$ partitions. Since $\bm{\alpha_{h}^{*}}$ and $\theta_{h}^{*}$, can be solved separately for each partition, $R_{m_{b},i_{h}}(\theta_{b},\theta_{h},\bm{\alpha_{b}^{*}},\bm{\alpha_{h}^{*}})$ is a function of $\theta_{b}$ and $\theta_{h}$ only. Similarly, $R_{m_{b},i_{h}}(\theta_{b},\bm{\alpha_{b}^{*}})$ is a function of $\theta_{b}$ only. Optimization problem (9) can be rewritten as

[TABLE]

By solving the optimization model (18), we can get the optimal $\theta_{h}^{*}$ for all $h=1,...,l$.

Theorem 5**.**

The optimal allocation that asymptotically minimizes the $PFS$ for the problem (6) is that $\theta_{h}^{*}=\frac{\gamma_{h}^{*}}{\sum\limits_{i=1}^{l}\gamma_{i}^{*}},$ where {\gamma_{b}^{*}}=\sigma_{b}\sqrt{\Big{(}\frac{\eta_{m_{b},1}^{2}}{\alpha^{*}_{1_{b}}}\!+\!\frac{\eta_{m_{b},s}^{2}}{\alpha^{*}_{s_{b}}}\!+\!\frac{\eta_{m_{b},k}^{2}}{\alpha^{*}_{k_{b}}}\!\Big{)}\sum\limits_{h\neq b}\frac{{\gamma_{h}^{*}}^{2}}{\sigma_{h}^{2}}} and $\gamma_{h}^{*}=\sigma_{h}^{2}/(f(x_{m_{b}})-f(x_{i_{h}^{*}}))^{2},\ h\neq b.$

The results of Theorem 5 indicate the optimal simulation budget allocation between the $l$ partitions. The results of Theorem 4 determine which designs should be selected for simulation and the number of simulation replications allocated to these selected designs. Since the solution space is divided into adjacent partitions, this optimal budget allocation procedure can be used for more general cases where the underlying function is non-quadratic. By conducting simulation replications on fewer designs, the simulation efficiency is dramatically enhanced compared to the existing methods.

4 Sequential Budget Allocation Algorithm

Based on the discussion above, we propose two sequential budget allocation algorithm to implement the optimality conditions. They are Optimal Optimal Computing Budget Allocation for selecting the top-$m$ designs with regression (OCBA-mr) based on Theorems 1 and 2 and Optimal Optimal Computing Budget Allocation for selecting the top-$m$ designs with regression in partitioned domains (OCBA-mrp) based on Theorems 4 and 5.

(1) OCBA-mr Algorithm

INITIALIZE: $\kappa\leftarrow 0$; Perform $n_{0}$ simulation replications for three design locations in each partition; by convention we use the D-optimal support designs, i.e., $n_{1}^{\kappa}=n_{{((1+t)/2)}}^{\kappa}=n_{t}^{\kappa}=n_{0}$ (round as needed). The incremental simulation budget is $\Delta$.

LOOP: WHILE $\sum_{i=1}^{t}n_{i}^{\kappa}<n$ DO

UPDATE: Estimate a quadratic regression equation $\hat{f}(x_{i})$ using the least squares method based on the information from all prior simulation runs.

Calculate the mean and variance of each design using $\hat{f}(x_{i})=\hat{\beta}_{0}+\hat{\beta}_{1}x_{i}+\hat{\beta}_{2}x_{i}^{2}.$

Determine the observed subset $\hat{S}_{o}$ and $\hat{S}_{n}$; the design with the $m$-th smallest sample mean value; the observed key design $\hat{i}^{*}$.

Calculate $\bm{\hat{\alpha}^{*}}$ and determine the support designs using Theorems 1 and 2.

ALLOCATE: Increase the simulation budget by $\Delta$ and calculate the new budget allocation rule using $\bm{\hat{\alpha}^{*}}$.

SIMULATE: Perform $\max(0,n_{i}^{\kappa}-n_{i}^{\kappa-1})$ simulation replications for design ${i}$ and $i=1,s,t$. $\kappa\leftarrow\kappa+1$.

END OF LOOP.

(2) OCBA-mrp Algorithm

INITIALIZE: $\kappa\leftarrow 0$; Perform $n_{0}$ simulation replications for three design locations in each partition; by convention we use the D-optimal support designs, i.e., $n_{1_{h}}^{\kappa}=n_{{((1+k^{(h)})/2)}_{h}}^{\kappa}=n_{k_{h}}^{\kappa}=n_{0}$ (round as needed). The incremental simulation budget is $\Delta$.

LOOP: WHILE $\sum_{h=1}^{l}\sum_{i_{h}=1_{h}}^{k_{h}}n_{i_{h}}^{\kappa}<n$ DO

UPDATE: Estimate the quadratic regression equations $\hat{f}(x_{i_{h}})$ for each partition using the least squares method based on the information from all prior simulation runs.

Calculate the means and variances of each design using $\hat{f}(x_{i_{h}})=\hat{\beta}_{h0}+\hat{\beta}_{h1}x_{i_{h}}+\hat{\beta}_{h2}x_{i_{h}}^{2},\ i_{h}=1_{h},2_{h},...,k_{h},h=1,2,...,l.$

Determine the observed subset $\hat{S}_{o}$ and $\hat{S}_{n}$; the design with the $m$-th smallest sample mean value; the observed key design $\hat{i_{h}}^{*}$ of each partition.

Calculate $\bm{\hat{\alpha}_{h}^{*}}$ and determine the support designs in each partition using Theorem 4. Calculate ${\hat{\theta}_{h}^{*}}$ using Theorem 5.

ALLOCATE: Increase the simulation budget by $\Delta$ and calculate the new budget allocation rule using $\bm{\hat{\alpha}_{h}^{*}}$ and $\hat{\theta}_{h}^{*}$.

SIMULATE: Perform $\max(0,n_{i_{h}}^{\kappa}-n_{i_{h}}^{\kappa-1})$ simulation replications for design ${i_{h}}$ in each partition $h=1,2,...,l$ and $i_{h}=1_{h},s_{h},k_{h}$. $\kappa\leftarrow\kappa+1$.

END OF LOOP.

5 Numerical Experiments

In this section, we test our proposed simulation budget allocation rules, OCBA-mr and OCBA-mrp, with some existing methods on several typical optimal subset selection problems. We use the following three budget allocation approaches for comparison. The top-$m$ designs are selected based on their mean values.

• •

Equal Allocation (EA): The EA is the most commonly used and simplest method, which allocates the simulation budget equally to each of the design.

• •

OCBAm+ Allocation: The OCBAm+ is an efficient simulation budget allocation procedure for selecting the top $m$ designs. It was developed based on the OCBA framework. For detail see Zhang et al., (2016).

• •

OCBA-ss Allocation: The OCAB-ss procedure was proposed in Gao and Chen, (2016) and seek to solve the optimal subset R&S problems for general underlying distributions using the large deviations theory.

• •

OSD Allocation: The OSD procedure was proposed in Brantley et al., (2013) and seek to solve the single best R&S problems using regression metamodel.

The OCBA-mrp, OCBA-mr and OSD estimate the performance of each design based on some regression metamodels. The EA, OCBAm+ and OCBA-ss use the mean value of each design for comparison, and the mean performance value is computed directly from the simulation output. Unlike our proposed methods and OSD, they do not rely on any response surface to estimate the performance value for any design. In order to compare the performance of these allocation approaches, we test them empirically on the following experiments.

• •

Experiment 1: Consider the optimization problem

[TABLE]

The experiment assumes that the noise for simulation has normal distribution $N(0,2^{2})$ for the objective value. We discretize the domain of the function into 100 evenly spaced points from 0 to 10. Since the underlying function is quadratic, we can use the OCBA-mr method directly for this problem. In order to utilize the OCBA-mrp method, we divide the 100 points into 5 adjacent partitions and there are 20 points in each partition. We want to select the top $m=5$ designs.

• •

Experiment 2: Consider the Griewank function

[TABLE]

We discretize the domain of the function into 100 evenly spaced points from 0 to 20. The nature of this underlying function does not allow us to apply the OCBA-mr method directly. To test the performance of the OCBA-mr method for the non-quadratic problem, we make a slight modification by dividing the solution space into 5 adjacent partitions. Then, we allocate the simulation budget equally to each partition and use the OCBA-mr method within each partition. In this experiment, the OCBA-mrp is also tested under this partition pattern. The experiment assumes that the noise for simulation has normal distribution $N(0,0.2^{2})$ We want to select the top $m=3$ designs.

• •

Experiment 3: Consider the optimization problem

[TABLE]

The entire domain consists of 200 discrete points for $x\in[0,8]$. The 200 points are further divided into 10 adjacent partitions so that there are 20 points in each partition. This experiment assumes that the simulation noise is a standard normal random variable. Similar to the setting in experiment 2, the OCBA-rm method is modified to handle this non-quadratic problem. We want to select the top $m=5$ designs.

• •

Experiment 4: Consider the optimization problem

[TABLE]

The experiment assumes that the noise for simulation has standard normal distribution for the objective measure. We discretize the domain of the function into 200 evenly spaced points from 0 to 2. The 200 points are further divided into 20 adjacent partitions so that there are 10 points in each partition. Similar to the setting in experiment 2, the OCBA-rm method is modified to handle this non-quadratic problem. We want to select the top $m=3$ designs.

• •

Experiment 5: Consider the optimization problem

[TABLE]

The experiment assumes that the noise for simulation has normal distribution $N(0,2^{2})$ for the objective measure. $x_{1}$ and $x_{2}$ are continuous variables with $-5\leq x_{1}\leq 5$ and $-5\leq x_{2}\leq 5$. We discretize the solution space into $11\times 11$ discrete points $x_{1}=\{-5,-4,...,5\}$ and $x_{2}=\{-5,-4,...,5\}$. We divide the $11\times 11$ designs into 11 adjacent partitions with 11 designs in each partition. For each $i=1,2,...,11$, partition $i$ consists of design points $(-5,i-6)$, $(-4,i-6)$, …, $(5,i-6)$. We want to select the top $m=3$ designs, which are $(0,-1)$, $(0,0)$ and $(0,1)$.

• •

Experiment 6: Consider the $(s,S)$ inventory problem in the Simulation Optimization Library (http://simopt.org/wiki/index.php?title=SS_Inventory). The decision variable $(s,S)$ is the inventory strategy. When the inventory of design $(s,S)$ on hand below $s$, an order is incurred to supplement the inventory to $S$. The optimization problem is to minimize the the E[Total cost per period]. The domain is discretized into $20\times 20$ discrete points with $s=\{810,820,...,1000\}$ and $S=\{1510,1520,...,1700\}$. The discrete points are divided into 20 partitions with 20 points in each partition. For each $i=1,2,...,20$, partition $i$ consists of points $(800+10i,1510),(800+10i,1520),...,(800+10i,1700)$. In each partition, we assume the performances of the E[Total cost per period] of each point can be modeled as quadratic functions, where the independent variable is $S$ and $s$ is kept as a constant. We want to select the top $m=3$ designs, which are $(810,1510)$, $(810,1520)$ and $(820,1510)$.

For the budget allocation procedures mentioned above, we let the initial number of replication $n_{0}$ be 10 and the incremental budget $\Delta$ be 100 for all the five experiments. Figures 1 reports the comparison results in terms of $PCS$. The estimate of $PCS$ is based on the average of 2000 independent replications of each procedure for experiments 1-6.

It is observed that all the procedures obtain higher $PCS$ as the available simulation budget increases. Our proposed procedures OCBA-mr and OCBA-mrp perform the best on the tested examples. It shows that our proposed procedures have dramatically improved the selection quality of the top-$m$ designs compared to the existing methods. When the considered problem is quadratic or approximately quadratic, such as example 1, the OCBA-mr method has the best performance. That is because the OCBA-mr method is developed for the quadratic or approximately quadratic problem setting. It can make full use of its advantages by simulating fewer designs once the quadratic assumption is satisfied. When the considered problem is partially quadratic or non-quadratic, such as examples 2-6, the OCBA-mrp method performs better than the modified OCBA-mr method. Because the OCBA-mrp can intelligently allocate the simulation budget both between and within each partition.

Compared to the existing approaches, our proposed approaches incorporate the information from the simulated designs into regression metamodel(s). We only need to simulate a subset of the alternative designs to build up the regression metamodel(s). The performances of the designs that have not been simulated can be inferred based on the metamodel(s). It dramatically improves the selection efficiency.

6 Conclusions

In this study, we further enhanced the simulation efficiency of finding the top-$m$ designs by incorporating quadratic regression equation(s). Using the large deviations theory, we formulate the problem of selecting the top-$m$ designs as that of maximizing the minimum convergence rate of the false selection probability. We derived two asymptotically optimal budget allocation rules, OCBA-mr and OCBA-mrp, for one partition and multi-partition problem setting, respectively. Numerical results suggest that the proposed approaches can dramatically improve the selection efficiency compared to the existing R&S approaches. When the underlying function can be approximated by a quadratic function across the solution space, the OCBA-mr is more efficient. For the non-quadratic problem setting, the OCBA-mrp demonstrates the best performance.

Proof of Lemma 2

Let $\Lambda_{i}^{(d)}(\theta)=\ln E(e^{\theta\hat{d}_{i,m}})$ denote the log-moment generating function of $\hat{d}_{i,m}$ and $I_{i}^{(d)}(\cdot)$ denote the Fenchel-Legendre transform of $\Lambda_{i}^{(d)}$, that is

[TABLE]

For any set $A$, let $A^{\circ}$ denotes its interior and for any function $f(\cdot)$, let $f^{\prime}(x)$ denote the derivative of $f$ at $x$. The effective domain of $\Lambda_{i}$ is $D_{\Lambda_{i}}=\{\theta\in\mathbb{R}:\Lambda_{i}(\theta)<\infty\}$ and $F_{i}=\{\Lambda_{i}^{\prime}(\theta):\theta\in D_{\Lambda_{i}}^{\circ}\}$. Let $f_{min}=\min_{i\in\{1,2,...,t\}}f(x_{i})$ and $f_{max}=\max_{i\in\{1,2,...,t\}}f(x_{i})$. We assume that the interval $[f_{min},f_{max}]\subset\bigcap_{i=1}^{t}F_{i}^{\circ}$. This assumption can be satisfied by some common families of distributions such Normal, Bernoulli, Poisson and Gamma family (Glynn and Juneja,, 2004). It ensures that the range of $\hat{f}(x_{i})$ include $[f_{min},f_{max}]$ and $P(\hat{f}(x_{i})>\hat{f}(x_{m}))$ and $P(\hat{f}(x_{j})<\hat{f}(x_{m}))$ are positive for $i\in S_{o}$ and $j\in S_{n}$.

Then, $\hat{d}_{i,m}=\hat{f}(x_{i})-\hat{f}(x_{m})$ satisfies the large deviation principle with good rate functions

[TABLE]

Since $\hat{d}_{i,m}=\hat{f}(x_{i})-\hat{f}(x_{m})$ are all normally distributed random variables, the rate function can be expressed as follow according to the results in Glynn and Juneja, (2004).

[TABLE]

Proof of Lemma 3

According to the definition of $PCS$, we have

[TABLE]

The last inequality follows from the Bonferroni inequality. Since $PFS=1-PCS$, we can get that

[TABLE]

Therefore, a lower bound on $PFS$ can be derived as

[TABLE]

and a upper bound on $PFS$ can be derived as

[TABLE]

According to Lemma 2, as $n$ goes to infinity we can get the convergence rate function shown in Lemma 3.

Proof of Theorem 1

Given the definition of the key design $i^{*}$, we rewrite (3) as

[TABLE]

Since ${R}_{m,i^{*}}(\bm{\alpha})$ is concave and strictly increasing functions of $\bm{\alpha}$ (Xiao et al.,, 2015; Glynn and Juneja,, 2004), the optimization problem is a convex programming problem. We define a Lagrangian function $L=-{R}_{m,i^{*}}(\bm{\alpha})+\lambda(\alpha_{1}+\alpha_{s}+\alpha_{t}-1)$ and determine the optimal allocation based on its Karush-Kuhn-Tucker (KKT) conditions.

[TABLE]

Plug (20) into $\alpha_{1}^{*}+\alpha_{s}^{*}+\alpha_{t}^{*}=1$, we can establish

[TABLE]

Then, we can get the conclusions in Theorem 1.

Proof of Theorem 2

Plugging (4) into ${R}_{m,i^{*}}(\bm{\alpha})$, we have

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

In order to analyze the optimal location for $x_{s}$, we consider the following five cases.

Case I: $\frac{x_{i^{*}}+x_{m}}{2}<\frac{3x_{1}+x_{t}}{4}$

We have ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}>0$ when $x_{s}<{(x_{1}+x_{t})/2}$; ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}<0$ when $x_{s}>{(x_{1}+x_{t})/2}$. Therefore, when ${(x_{i^{*}}+x_{m})}/{2}<{(3x_{1}+x_{t})}/{4}$, we choose $x_{s}={(x_{1}+x_{t})/2}$.

Case II: $\frac{3x_{1}+x_{t}}{4}\leq\frac{x_{i^{*}}+x_{m}}{2}<\frac{x_{1}+x_{t}}{2}$

We have ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}>0$ when $x_{s}<x_{i^{*}}+x_{m}-x_{1}$; ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}<0$ when $x_{s}>x_{i^{*}}+x_{m}-x_{1}$. Therefore, when ${(3x_{1}+x_{t})}/{4}\leq{(x_{i^{*}}+x_{m})}/{2}<{(x_{1}+x_{t})}/{2}$, we choose $x_{s}=x_{i^{*}}+x_{m}-x_{1}$.

Case III: $\frac{x_{i^{*}}+x_{m}}{2}>\frac{x_{1}+3x_{t}}{4}$

We have ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}>0$ when $x_{s}<{(x_{1}+x_{t})/2}$; ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}<0$ when $x_{s}>{(x_{1}+x_{t})/2}$. Therefore, when ${(x_{i^{*}}+x_{m})}/{2}>{(x_{1}+3x_{t})}/{4}$, we choose $x_{s}={(x_{1}+x_{t})/2}$.

Case IV: $\frac{x_{1}+x_{t}}{2}<\frac{x_{i^{*}}+x_{m}}{2}\leq\frac{x_{1}+3x_{t}}{4}$

We have ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}>0$ when $x_{s}<x_{i^{*}}+x_{m}-x_{t}$; ${\partial{R}_{m,i^{*}}(\bm{\alpha})}/{\partial x_{s}}<0$ when $x_{s}>x_{i^{*}}+x_{m}-x_{t}$. Therefore, when ${(3x_{1}+x_{t})}/{4}\leq{(x_{i^{*}}+x_{m})}/{2}<{(x_{1}+x_{t})}/{2}$, we choose $x_{s}=x_{i^{*}}+x_{m}-\!x_{t}$.

Case V: $\frac{x_{i^{*}}+x_{m}}{2}=\frac{x_{1}+x_{t}}{2}$

The location of $x_{s}$ has no influence on the results. Therefore, we choose $x_{s}={(x_{1}+x_{t})/2}$ for consistency with Case II and Case III.

Analyzing the above the results, we can get the conclusions in Theorem 2.

Proof of Lemma 4

Let $\Lambda_{i_{h}}(\theta)=\ln E(e^{\theta\hat{f}(x_{i_{h}})})$ denote the log-moment generating function of $\hat{f}(x_{i_{h}})$ and $I_{i_{h}}(\cdot)$ denote the Fenchel-Legendre transform of $\Lambda_{i_{h}}$, that is

[TABLE]

Let the scaled cumulant generating function of $(\hat{f}(x_{m_{b}}),\hat{f}(x_{i_{h}}))$ be denoted as follows:

[TABLE]

By the Gärtner-Ellis theorem, $(\hat{f}(x_{m_{b}}),\hat{f}(x_{i_{h}}))$ satisfy the large deviation principle with good rate functions

[TABLE]

Since $(\hat{f}(x_{m_{b}}),\hat{f}(x_{i_{h}}))$ are all normally distributed random variables, the rate function can be expressed as follows according to the results in Glynn and Juneja, (2004). Then, we can get

[TABLE]

The proof of (8) is similar to that is given for Lemma 2, and hence is omitted for brevity.

Proof of Lemma 5

We first rewrite optimization model (10) as

[TABLE]

Similarly, the optimization models (11) and (12) can be rewritten as follows:

[TABLE]

for $h=1,2,...,m,h\neq b$, and

[TABLE]

Let $\bm{\widetilde{\alpha}_{h}^{*}}$ and $\bm{\widetilde{\alpha}_{b}^{*}}$ be the optimal solutions to (23) and (24). The optimization model (22), (23) and (24) have the same objective function, and the domain of (22) is a subset of (23) and (24). Therefore, if $\bm{\widetilde{\alpha}_{h}^{*}}$ and $\bm{\widetilde{\alpha}_{b}^{*}}$ are feasible to (22), they are the optimal solutions of (22). According to the formula of $R_{i_{b}}(\theta_{b},\bm{\alpha_{b}})$ and $\widetilde{R}_{i_{h}}(\theta_{h},\bm{\alpha_{h}})$, we can conclude that $\bm{\widetilde{\alpha}_{h}^{*}}$ and $\bm{\widetilde{\alpha}_{b}^{*}}$ do not depend on the value of $\theta_{h}$. Therefore, we can conclude that $\bm{\widetilde{\alpha}_{h}^{*}}$ and $\bm{\widetilde{\alpha}_{b}^{*}}$ are feasible to (22), and they are optimal to (10). Then we get the conclusions in Lemma 5.

Proof of Theorem 3

In order to determine the optimal $\bm{\alpha_{h}^{*}}$ for each partition, we consider two cases ($h\neq b$ and $h=b$).

For $h\neq b$, we solve the optimization model (11). Given the definition of the key design $i_{h}^{*}$, we rewrite (11) as

[TABLE]

for $h=1,2,...,m$ and $h\neq b$.

Since $\widetilde{R}_{i_{h}^{*}}(\theta_{h},\bm{\alpha_{h}})$ is concave and strictly increasing functions of $\bm{\alpha_{h}}$ (Xiao et al.,, 2015; Glynn and Juneja,, 2004), the optimization problem is a convex programming problem. We define a Lagrangian function $L\!=\!-\widetilde{R}_{m_{b},i_{h}^{*}}(\theta_{h},\bm{\alpha_{h}})+\lambda_{h}(\alpha_{1_{h}}\!+\!\alpha_{s_{h}}\!+\!\alpha_{k_{h}}-1)$ and determine the optimal allocation based on its Karush-Kuhn-Tucker (KKT) conditions.

[TABLE]

[TABLE]

Using the fact that $\alpha_{1_{h}}^{*}+\alpha_{s_{h}}^{*}+\alpha_{k_{h}}^{*}=1$, we can establish

[TABLE]

For $h=b$, we solve the optimization model (12), which could be rewritten as

[TABLE]

The process of solving the (27) is similar to that is given for Theorem 1, and hence is omitted for brevity. Then, we can get the conclusions in Theorem 3.

Proof of Theorem 4

In order to determine the location of the support design $s_{h}$ for each partition, we consider two cases ($h\neq b$ and $h=b$).

For $h\neq b$, we plug (13) into $\widetilde{R}_{m_{b},i_{h}^{*}}(\theta_{h},\bm{\alpha_{h}})$, we have

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

When $x_{s_{h}}>x_{i_{h}^{*}}$,

[TABLE]

When $x_{s_{h}}<x_{i_{h}^{*}}$,

[TABLE]

Therefore, $\widetilde{R}_{m_{b},i_{h}^{*}}(\theta_{h},\bm{\alpha_{h}})$ is maximized when

[TABLE]

According to equation (13), we can get

[TABLE]

When $i_{h}^{*}$ is at the extreme locations ($1_{h}$ and $k_{h}$),

[TABLE]

In this setting, the location of the support design $s_{h}$ has no effect on $\widetilde{R}_{m_{b},i_{h}^{*}}(\theta_{h},\bm{\alpha_{h}})$. For simplicity, we use the D-optimal support design and let $x_{s_{h}}={(x_{1_{h}}+x_{k_{h}})/2}$ (round to the nearest design as needed) (Kiefer,, 1959).

For $h=b$, the proof is similar to that is given for Theorem 2, and hence is omitted for brevity. Then, we can get the conclusions in Theorem 4.

Proof of Theorem 5

Given the definition of key design $i_{h}^{*}$, we rewrite the optimization model (18) as

[TABLE]

The $R_{m_{b},i_{b}^{*}}(\theta_{b},\bm{\alpha_{b}^{*}})$ and $R_{m_{b},i_{h}^{*}}(\theta_{b},\theta_{h},\bm{\alpha_{b}^{*}},\bm{\alpha_{h}^{*}})$ are concave and strictly increasing functions of $\theta_{h}$ and $\theta_{b}$ (Glynn and Juneja,, 2004). Therefore, the optimization problem is a convex programming problem, and the first order condition is also the optimality condition.

Rewrite the optimization model above as

[TABLE]

From the KKT conditions, there exist $\lambda_{h}\geq 0,\ h=1,2,...,l$ and $\mu>0$ such that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From (29), it can be concluded that there must exist some $\lambda_{h}>0,h=1,2,...,l$. If there is one $\lambda_{h}=0,h=1,2,...,l,h\neq b$, it results $\mu=0$ based on (30). That means all the $\lambda_{h}=0$, as ${\partial R_{m_{b},i_{h}^{*}}(\theta_{b},\theta_{h},\bm{\alpha_{b}^{*}},\bm{\alpha_{h}^{*}})}/{\partial\theta_{h}}$ is strictly positive. Therefore, we can conclude that $\lambda_{h}\neq 0,\ h=1,2,...,l,h\neq b.$ According to (32),

[TABLE]

In addition, since $\theta_{h}/\theta_{b}\rightarrow 0$ as $l$ goes to infinity and $R_{m_{b},i_{b}^{*}}(\theta_{b}^{*},\bm{\alpha_{b}^{*}})>R_{m_{b},i_{h}^{*}}(\theta_{b}^{*},\theta_{h}^{*},\bm{\alpha_{b}^{*}},\bm{\alpha_{h}^{*}}),\ h=1,2,...,l,h\neq b$, it can be concluded that $\lambda_{b}=0$ according to (33). Therefore, we have

[TABLE]

According to Theorem 4, we have

[TABLE]

Then, we can get

[TABLE]

and

[TABLE]

As $l$ goes to infinity, plug the results above into (34) and (35), we can get the equations in Theorem 5.