Efficient Simulation Budget Allocation for Ranking the Top 𝑚 Designs

Abstract

We consider the problem of ranking the top 𝑚 designs out of 𝑘 alternatives. Using the optimal computing budget allocation framework, we formulate this problem as that of maximizing the probability of correctly ranking the top 𝑚 designs subject to the constraint of a fixed limited simulation budget. We derive the convergence rate of the false ranking probability based on the large deviation theory. The asymptotically optimal allocation rule is obtained by maximizing this convergence rate function. To implement the simulation budget allocation rule, we suggest a heuristic sequential algorithm. Numerical experiments are conducted to compare the effectiveness of the proposed simulation budget allocation rule. The numerical results indicate that the proposed asymptotically optimal allocation rule performs the best comparing with other allocation rules.

1. Introduction

Discrete event system (DES) simulation has been widely used for analyzing and evaluating complex systems since the assumptions for deriving an analytical solution are rarely satisfied in real situation. While DES simulation has been successful in solving many practical problems in a variety of areas such as supply chain systems, healthcare systems, and manufacturing systems, the concerns on the efficiency have never ended [1]. To obtain a statistically significant value, a large number of simulation replications are needed for each design. The performance of each design is then estimated by its sample mean. The ultimate accuracy of this estimator cannot be improved faster than (1/√𝑁), where 𝑁 is the number of simulation replications.

Ordinal optimization (OO) which aims to obtain a good estimate through ordinal comparisons although the estimated value is still poor emerges as a way to improve the simulation efficiency [2]. If the goal of the simulation experiment is to identify the good designs instead of finding the accurate estimate of the true performance value, which is true in many real applications, OO can reduce the number of simulation replications significantly [3, 4]. Intuitively, a larger portion of the total simulation replications should be allocated to those designs that are critical in identifying the best design in order to achieve a high probability of correct selection. Based on this idea, an optimal computing budget allocation (OCBA) framework has been proposed to enhance the simulation efficiency further [5, 6]. OCBA focuses on the efficiency of simulation by intelligently allocating further replications based on both the mean and the variance. In parallel with OCBA, two other well-known ranking and selection procedures frequently used in simulation are the indifference zone (IZ) procedure and value information procedure (VIP). The IZ procedure focuses on finding a feasible way to guarantee that the prespecified probability of correct selection can be achieved [7]. The VIP uses the Bayesian posterior distribution to describe the evidence of correct selection and allocates further replications based on maximizing the value information [8, 9]. The three popular procedures of finding the best design are compared in [10]. Extensions of the research in OCBA and ranking and selection include subset selection [11–13], selecting the Pareto set for multiple objective functions [14], selecting the best design subject to stochastic constraints [15, 16], and complete ranking [17]. A detailed summary of the existing work in OCBA can be found in [18].

To the best of our knowledge, no previous research has considered the simulation budget allocation for ranking the top 𝑚 designs out of 𝑘 alternatives. The developed top 𝑚 ranking procedure is most useful when the designs for comparison have multiple dimensions of performance measurements with some qualitative criteria such as environmental consideration and political feasibility. Top 𝑚 ranking procedure can help the decision makers to identify the ranking of the top 𝑚 designs based on the quantitative performance measurement. Final decision can be made by incorporating other qualitative performance measurements. Hence, the top 𝑚 ranking gives decision makers a more flexible and people oriented way to support the decision making process. In addition, the top 𝑚 ranking procedure can also be integrated into some evolutionary search algorithms, where the ranking information of the top solutions is needed in order to determine the search direction of next iteration. For example, the ranking of the top candidate solutions may be needed in each iteration of the genetic algorithm since better candidate solutions are given higher probabilities to reproduce in genetic algorithms.

In this paper, we consider the problem of ranking the top 𝑚 designs out of total 𝑘 alternatives, where the performance value of each design can only be estimated with noise via simulation. The objective of this paper is to determine how to allocate the simulation replications among the 𝑘 designs such that the probability of correctly ranking the top 𝑚 designs can be maximized. The organization of this paper is as follows. Section 2 provides the mathematical formulation of the top 𝑚 ranking problem using OCBA framework. In Section 3, we derive the convergence rate function of the false ranking probability. Section 4 derives the asymptotically optimal allocation rule using large deviation theory. A sequential allocation algorithm is proposed in Section 5. Section 6 conducts numerical experiments to demonstrate the effectiveness of the proposed simulation budget allocation rule. Finally, we conclude our paper in Section 7.

2. Problem Formulation

Consider the problem of ranking the top 𝑚 designs out of 𝑘 alternatives. The performance of each design can only be estimated with noise via simulation. The mean performance is used as the ranking criterion. In order to have a steady mean performance value, a large number of simulation replications are needed because of the randomness of individual samples. Given that a total of 𝑛 simulation replications are available, the objective is to find the best allocation of the total 𝑛 replications to the 𝑘 designs in order to maximize the probability of correctly ranking the top 𝑚 designs.

Denote the mean performances of each design by 𝜇₁,…,𝜇_𝑘 such that 𝜇₁<⋯<𝜇_𝑖<⋯<𝜇_𝑘 and 𝜇_𝑖+1−𝜇_𝑖 ≥ 𝛿, 𝑖=1,…,𝑘−1, where 𝛿 is a positive number. Let 𝜶=(𝛼₁,…,𝛼_𝑘) denote the proportion of the total computing budget 𝑛 to be allocated to each design such that ∑^𝑘_𝑖=1𝛼_𝑖 = 1, 𝑖=1,…,𝑘. Let 𝑋_𝑖(𝛼_𝑖𝑛)=(𝛼_𝑖𝑛)⁻¹∑^𝛼𝑖𝑛_𝑗=1𝑋_𝑖𝑗 denote the sample mean of design 𝑖, where (𝑋_𝑖1,…,𝑋_{𝑖,𝛼𝑖𝑛}) are the samples from population 𝑖. We ignore the case when 𝛼_𝑖𝑛 is not an integer because we could let 𝛼_𝑖𝑛 be ⌊𝛼_𝑖𝑛⌋, the greatest integer less than 𝛼_𝑖𝑛. The analysis remains unaffected. The objective of this paper is to find the optimal allocation strategy 𝜶^∗=(𝛼^∗₁,…,𝛼^∗_𝑘) such that the probability of correctly ranking the top 𝑚 designs can be maximized with a fixed limited computing budget 𝑛.

Given that 𝜇₁<⋯<𝜇_𝑖<⋯<𝜇_𝑘, the ranking of the top 𝑚 designs is correctly identified if 𝑋_𝑖(𝛼_𝑖𝑛)≤𝑋_𝑖+1(𝛼_𝑖+1𝑛) for all 𝑖=1,…,𝑚−1 and if 𝑋_𝑚(𝛼_𝑚𝑛)≤𝑋_𝑗(𝛼_𝑗𝑛) for all 𝑗=𝑚+1,…,𝑘. Mathematically, we can write the probability of correctly ranking the top 𝑚 designs as follows: 𝑃CR_𝑚=𝑃^𝑚−1_𝑖=1𝑋_𝑖𝛼_𝑖𝑛≤𝑋_𝑖+1𝛼_𝑖+1𝑛_{𝑘𝑗=𝑚+1}𝑋_𝑚𝛼_𝑚𝑛≤𝑋_𝑗𝛼_𝑗𝑛. (1)

Hence, an optimal computing budget allocation problem can be formulated by maximizing the probability of correctly ranking the top 𝑚 designs: max_{𝛼1,…,𝛼𝑘}𝑃CR_𝑚s.t.𝛼₁+⋯+𝛼_𝑖+⋯+𝛼_𝑘=1,𝛼_𝑖≥0,𝑖=1,…,𝑘. (2)

Maximizing the probability of correctly ranking the top 𝑚 designs is equivalent to minimizing the false ranking probability of the top 𝑚 designs. Using asymptotical analysis, this is also equivalent to maximizing the convergence rate at which the false ranking probability goes to zero. In this paper, we use large deviation theory to derive this convergence rate function and reformulate the optimization model (2) by maximizing the convergence rate function. The assumptions needed in this paper are stated as follows.

Assumption 1. The performance of every design is independently simulated.
The independence of each design ensures that the samples (𝑋_𝑖1,…,𝑋_{𝑖,𝛼𝑖𝑛}) for each 𝑖=1,…,𝑘 generated are independent. Thus, the results we obtained will not be affected by the correlations among different designs.
Define the cumulant generating function of sample mean 𝑋_𝑖(𝑛) to be Λ_𝑖^(𝑛)(𝜃)=ln𝐸(𝑒^{𝜃𝑋𝑖(𝑛)}). The effective domain of any function 𝑓∶𝐷_𝑓→𝑅^∗ is the set {𝑥∈𝐷_𝑓∶𝑓(𝑥)<∞}, while the range is 𝑅^∗=𝑅∪{+∞}. Let 𝐷_Λ𝑖={𝜃∈𝑅∶Λ_𝑖(𝜃)<∞} and 𝐹_𝑖={Λ^′_𝑖(𝜃)∶𝜃∈𝐷⁰_Λ𝑖}. For any set 𝐴, 𝐴^𝑜 denotes its interior and 𝐴 denotes its closure.

Assumption 2. For each 𝑖=1,…,𝑘,(1)the limit Λ_𝑖(𝜃)=lim_𝑛→∞(1/𝑛)Λ_𝑖^(𝑛)(𝑛𝜃) is well defined as an extended real number for all 𝜃;(2)the origin belongs to 𝐷⁰_Λ𝑖;(3)Λ_𝑖(⋅) is strictly convex and steep, that is, lim_𝑛→∞|Λ^′_𝑖(𝜃_𝑛)|=∞, where {𝜃_𝑛} is a sequence converging to the boundary point of 𝐷_Λ𝑖;(4)[𝜇₁,𝜇_𝑘]⊂⋂^𝑘_𝑖=1𝐹⁰_𝑖.

Assumption 2 implies that 𝜇_𝑖=Λ^′_𝑖(0) with 𝑋_𝑖(𝑛)→𝜇_𝑖 almost surely when 𝑛→∞. Furthermore, it indicates that the sample mean 𝑋_𝑖(𝛼_𝑖𝑛) satisfies the large deviation principle. The last condition of Assumption 2 ensures that the sample mean of every design can take any value between 𝜇₁ and 𝜇_𝑘 and 𝑃(𝑋_𝑖(𝛼_𝑖𝑛)≥𝑋_𝑖+1(𝛼_𝑖+1𝑛))>0 [19].

3. Rate Function of the False Ranking Probability

We now derive the probability of false ranking other than the correct ranking probability, followed by the corresponding derivation of its large deviation principle. Recall that the probability of correctly ranking the top 𝑚 designs is defined in (1). The probability of falsely ranking the top 𝑚 designs is simply its complement; that is, 𝑃FR_𝑚=1−𝑃CR_𝑚=𝑃^𝑚−1_𝑖=1𝑋_𝑖𝛼_𝑖𝑛≤𝑋_𝑖+1𝛼_𝑖+1𝑛_{𝑘𝑗=𝑚+1}𝑋_𝑚𝛼_𝑚𝑛≤𝑋_𝑗𝛼_𝑗𝑛^𝐶=𝑃^𝑚−1_𝑖=1𝑋_𝑖𝛼_𝑖𝑛≥𝑋_𝑖+1𝛼_𝑖+1𝑛_{𝑘𝑗=𝑚+1}𝑋_𝑚𝛼_𝑚𝑛≥𝑋_𝑗𝛼_𝑗𝑛, (3) where (⋅)^𝐶 represents the complement of (⋅).

𝑃(FR_𝑚) has a lower bound, lb=maxmax_{1≤𝑖≤𝑚−1}𝑃𝑋_𝑖𝛼_𝑖𝑛≥𝑋_𝑖+1𝛼_𝑖+1𝑛,max_{𝑗=𝑚+1,…,𝑘}𝑃𝑋_𝑚𝛼_𝑚𝑛≥𝑋_𝑗𝛼_𝑗𝑛, (4) and an upper bound ub=(𝑘−1)×lb such that, assuming the limit exists, lim_𝑛→∞1𝑛ln𝑃FR_𝑚=lim_𝑛→∞1𝑛lnlb. (5)

Theorem 3 formally states that the limit exists and the overall convergence rate function is the minimum rate function of each probability. The convergence rate can be understood as the speed at which false ranking probability goes to zero.

Theorem 3. The rate function of 𝑃(𝐹𝑅_𝑚) is given by −lim_𝑛→∞1𝑛ln𝑃𝐹𝑅_𝑚=minmin_{1≤𝑖≤𝑚−1}inf_𝑥𝛼_𝑖𝐼_𝑖(𝑥)+𝛼_𝑖+1𝐼_𝑖+1(𝑥),min_{𝑗=𝑚+1,…,𝑘}inf_𝑥𝛼_𝑗𝐼_𝑗(𝑥)+𝛼_𝑚𝐼_𝑚(𝑥), (6) where 𝐼(𝑥)=sup_𝜃∈𝑅(𝜃𝑥−Λ(𝜃)) is the Fenchel-Legendre transform and Λ(𝜃)=ln𝐸(𝑒^𝜃𝑋).

Proof. If there exists a function 𝐺_𝑖(⋅) such that lim_𝑛→∞1𝑛ln𝑃𝑋_𝑖𝛼_𝑖𝑛≥𝑋_𝑖+1𝛼_𝑖+1𝑛=−𝐺_𝑖𝛼_𝑖,𝛼_𝑖+1, (7) for 𝑖=1,…,𝑚−1, and lim_𝑛→∞1𝑛ln𝑃𝑋_𝑚𝛼_𝑚𝑛≥𝑋_𝑗𝛼_𝑗𝑛=−𝐺_𝑗𝛼_𝑗,𝛼_𝑚, (8) for 𝑗=𝑚+1,…,𝑘.
Then, it can be concluded that lim_𝑛→∞1𝑛ln𝑃FR_𝑚=−minmin_{1≤𝑖≤𝑚−1}𝐺_𝑖𝛼_𝑖,𝛼_𝑖+1,min_{𝑗=𝑚+1,…,𝑘}𝐺_𝑗𝛼_𝑗,𝛼_𝑚. (9)
Now we are in the position to derive the assumed function 𝐺_𝑖(⋅). Define 𝑌_𝑛=(𝑋_𝑖(𝛼_𝑖𝑛),𝑋_𝑖+1(𝛼_𝑖+1𝑛)). The cumulant moment generating function of 𝑌_𝑛 can be written as Λ_𝑛𝜃_𝑖,𝜃_𝑖+1=ln𝐸𝑒^{𝜃𝑖𝑋𝑖(𝛼𝑖𝑛)+𝜃𝑖+1𝑋𝑖+1(𝛼𝑖+1𝑛)}=Λ^(𝛼𝑖_𝑖^𝑛)𝜃_𝑖𝛼_𝑖𝑛+Λ^{(𝛼𝑖+1𝑛)}_𝑖+1𝜃_𝑖+1𝛼_𝑖+1𝑛. (10)
Under Assumption 2, lim_𝑛→∞1𝑛Λ_𝑛𝑛𝜃_𝑖,𝑛𝜃_𝑖+1=𝛼_𝑖Λ_𝑖𝜃_𝑖𝛼_𝑖+𝛼_𝑖+1Λ_𝑖+1𝜃_𝑖+1𝛼_𝑖+1. (11)
By the Gärtner-Ellis Theorem [20], {𝑌_𝑛,𝑛=1,2,…} satisfies large deviation principle with good rate function which can be expressed as follows: 𝐼𝑥_𝑖,𝑥_𝑖+1=sup_{𝜃𝑖,𝜃𝑖+1}𝜃_𝑖𝑥_𝑖+𝜃_𝑖+1𝑥_𝑖+1−𝛼_𝑖Λ_𝑖𝜃_𝑖𝛼_𝑖−𝛼_𝑖+1Λ_𝑖+1𝜃_𝑖+1𝛼_𝑖+1=sup_𝜃𝑖𝜃_𝑖𝑥_𝑖−𝛼_𝑖Λ_𝑖𝜃_𝑖𝛼_𝑖+sup_𝜃𝑖+1𝜃_𝑖+1𝑥_𝑖+1−𝛼_𝑖+1Λ_𝑖+1𝜃_𝑖+1𝛼_𝑖+1=𝛼_𝑖sup_{𝜃𝑖/𝛼𝑖}𝜃_𝑖𝛼_𝑖𝑥_𝑖−Λ_𝑖𝜃_𝑖𝛼_𝑖+𝛼_𝑖+1sup_{𝜃𝑖+1/𝛼𝑖+1}𝜃_𝑖+1𝛼_𝑖+1𝑥_𝑖+1−Λ_𝑖+1𝜃_𝑖+1𝛼_𝑖+1=𝛼_𝑖𝐼_𝑖𝑥_𝑖+𝛼_𝑖+1𝐼_𝑖+1𝑥_𝑖+1. (12)
Hence, from large deviation principle, lim_𝑛→∞1𝑛ln𝑃𝑋_𝑖𝛼_𝑖𝑛≥𝑋_𝑖+1𝛼_𝑖+1𝑛=−inf_𝑥𝛼_𝑖𝐼_𝑖(𝑥)+𝛼_𝑖+1𝐼_𝑖+1(𝑥). (13)
Similarly, lim_𝑛→∞1𝑛ln𝑃𝑋_𝑚𝛼_𝑚𝑛≥𝑋_𝑗𝛼_𝑗𝑛=−inf_𝑥𝛼_𝑗𝐼_𝑗(𝑥)+𝛼_𝑚𝐼_𝑚(𝑥). (14)
Therefore, the convergence rate function of the false ranking probability can be expressed as follows: −lim_𝑛→∞1𝑛ln𝑃FR_𝑚=minmin_{1≤𝑖≤𝑚−1}𝐺_𝑖𝛼_𝑖,𝛼_𝑖+1,min_{𝑗=𝑚+1,…,𝑘}𝐺_𝑗𝛼_𝑗,𝛼_𝑚=minmin_{1≤𝑖≤𝑚−1}inf_𝑥𝛼_𝑖𝐼_𝑖(𝑥)+𝛼_𝑖+1𝐼_𝑖+1(𝑥),min_{𝑗=𝑚+1,…,𝑘}inf_𝑥𝛼_𝑗𝐼_𝑗(𝑥)+𝛼_𝑚𝐼_𝑚(𝑥). (15)

4. Asymptotically Optimal Allocation

The objective is to maximize the probability of correctly ranking the top 𝑚 designs. This can be achieved by minimizing the false ranking probability. It is also equivalent to maximizing the convergence rate of 𝑃(FR_𝑚) subject to ∑^𝑘_𝑖=1𝛼_𝑖=1 and 𝛼_𝑖≥0, for all 𝑖=1,…,𝑘. The optimization model in (1) can be reexpressed as maxminmin_{1≤𝑖≤𝑚−1}inf_𝑥𝛼_𝑖𝐼_𝑖(𝑥)+𝛼_𝑖+1𝐼_𝑖+1(𝑥),min_{𝑗=𝑚+1,…,𝑘}inf_𝑥𝛼_𝑗𝐼_𝑗(𝑥)+𝛼_𝑚𝐼_𝑚(𝑥)s.t.^𝑘_𝑖=1𝛼_𝑖=1,𝛼_𝑖≥0,𝑖∈{1,…,𝑘}. (16)

By [19], 𝛼_𝑖𝐼_𝑖(𝑥)+𝛼_𝑖+1𝐼_𝑖+1(𝑥) is a strictly increasing concave function. The infimum of concave functions is also concave. Likewise, the minimum of concave functions is a concave function too. Define 𝑥(𝛼_𝑖,𝛼_𝑖+1) = arginf_𝑥(𝛼_𝑖𝐼_𝑖(𝑥)+𝛼_𝑖+1𝐼_𝑖+1(𝑥)). As shown in [19], 𝑥(𝛼_𝑖,𝛼_𝑖+1) is the solution to 𝛼_𝑖𝐼^′_𝑖(𝑥)  +  𝛼_𝑖+1𝐼^′_𝑖+1(𝑥)=0 and 𝜕𝛼_𝑖𝐼_𝑖𝑥𝛼_𝑖,𝛼_𝑖+1+𝛼_𝑖+1𝐼_𝑖+1𝑥𝛼_𝑖,𝛼_𝑖+1𝜕𝛼_𝑖=𝐼_𝑖𝑥𝛼_𝑖,𝛼_𝑖+1,𝜕𝛼_𝑖𝐼_𝑖𝑥𝛼_𝑖,𝛼_𝑖+1+𝛼_𝑖+1𝐼_𝑖+1𝑥𝛼_𝑖,𝛼_𝑖+1𝜕𝛼_𝑖+1=𝐼_𝑖+1𝑥𝛼_𝑖,𝛼_𝑖+1. (17)

The result can similarly be applied to 𝛼_𝑗𝐼_𝑗(𝑥)+𝛼_𝑚𝐼_𝑚(𝑥). Therefore, the optimization model (16) is a concave maximization problem and it can be reexpressed as follows: max𝑧s.t.𝛼_𝑖𝐼_𝑖𝑥𝛼_𝑖,𝛼_𝑖+1+𝛼_𝑖+1𝐼_𝑖+1𝑥𝛼_𝑖,𝛼_𝑖+1{}𝛼≥𝑧,𝑖∈1,…,𝑚−1_𝑗𝐼_𝑗𝑥𝛼_𝑗,𝛼_𝑚+𝛼_𝑚𝐼_𝑚𝑥𝛼_𝑗,𝛼_𝑚{}≥𝑧,𝑗∈𝑚+1,…,𝑘_𝑘𝑖=1𝛼_𝑖=1,𝛼_𝑖≥0,𝑖∈{1,…,𝑘}. (18)

Since model (18) is strictly concave and the functions of 𝜶 are continuous, a unique optimal solution must exist and the Karush-Kuhn-Tucker (KKT) conditions are necessary and sufficient for global optimality.

From the KKT conditions on problem (18), we define a new problem (19) by replacing some inequality signs and forcing 𝛼_𝑖 to be strictly positive: max𝑧s.t.𝛼₁𝐼₁𝑥𝛼₁,𝛼₂+𝛼₂𝐼₂𝑥𝛼₁,𝛼₂𝛼=𝑧min_𝑖𝐼_𝑖𝑥𝛼_𝑖,𝛼_𝑖+1+𝛼_𝑖+1𝐼_𝑖+1𝑥𝛼_𝑖,𝛼_𝑖+1,𝛼_𝑖−1𝐼_𝑖−1𝑥𝛼_𝑖−1,𝛼_𝑖+𝛼_𝑖𝐼_𝑖𝑥𝛼_𝑖−1,𝛼_𝑖{}𝛼=𝑧,𝑖∈2,…,𝑚−1_𝑗𝐼_𝑗𝑥𝛼_𝑗,𝛼_𝑚+𝛼_𝑚𝐼_𝑚𝑥𝛼_𝑗,𝛼_𝑚{}=𝑧,𝑗∈𝑚+1,…,𝑘_𝑘𝑖=1𝛼_𝑖=1,𝛼_𝑖>0,𝑖∈{1,…,𝑘}. (19)

Theorem 4. Under Assumptions 1 and 2 in Section 2, problems (18) and (19) are equivalent; that is, a solution 𝜶^∗=(𝛼^∗₁,…,𝛼^∗_𝑘) is the optimal solution to (18) if and only if it is also an optimal solution to (19).

Proof. We assume that a point satisfying the KKT condition of (18) is also feasible to (19). We first prove the forward and backward assertions. We then prove that the assumption that a point satisfying the KKT condition of (18) is also feasible to (19) is indeed correct.
Suppose 𝜶^∗ is the optimal solution to (18). Since the feasible region of (19) is a subset of that of (18), if the optimal solution to (18) is feasible to (19), it must be optimal to (19). Since the KKT conditions are necessary and sufficient for optimality in (18), 𝜶^∗must satisfy the KKT conditions of (18). Hence, 𝜶^∗ is feasible to (19). Therefore, if a point satisfies the KKT condition in (18), it must be optimal to (19).
Suppose the optimal solution to (18) is 𝜶^∗ and the optimal solution to (19) is 𝜶^∗, and 𝜶^∗≠𝜶^∗. Since the KKT conditions are necessary and sufficient condition to (18), thus, 𝜶^∗ must satisfy the KKT conditions. Furthermore, the objective function of (19) is the same as that of (18), and the feasible region of (19) is a subset of that of (18). Therefore, 𝜶^∗ must be infeasible to (19). However, we assumed that a point satisfying the KKT conditions of (18) must be feasible to (19). We have thus reached a contradiction. So we must have 𝜶^∗=𝜶^∗.
We are now in the position to prove that a point satisfying the KKT conditions of (18) must be feasible to (19). If we let 𝛼_𝑖=1/𝑘, we can have 𝑧>0 for problem (18). However, any 𝛼_𝑖=0 for 𝑖=1,…,𝑚−1 will lead to 𝛼_𝑖+1inf_𝑥𝐼_𝑖+1(𝑥)=𝛼_𝑖+1𝐼_𝑖+1(𝜇_𝑖+1)=0. If 𝛼_𝑗=0 for 𝑗=𝑚+1,…,𝑘, we will have 𝛼_𝑚inf_𝑥𝐼_𝑚(𝑥)=𝛼_𝑚𝐼_𝑚(𝜇_𝑚)=0. Therefore, the optimal solution must satisfy 𝛼_𝑖>0 for every 𝑖=1,…,𝑘.
Since the problem (18) is a concave optimization problem, the first order condition is also the optimality condition. According to the KKT conditions, there exist 𝜆_𝑖≥0, 𝜆_𝑗≥0, 𝑖∈{1,…,𝑚−1}, 𝑗∈{𝑚+1,…,𝑘}, and 𝛾>0 such that ^𝑚−1_𝑖=1𝜆_𝑖+^𝑘_𝑗=𝑚+1𝜆_𝑗=1,𝜆₁𝐼₁𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1=𝛾,𝜆_𝑖−1𝐼_𝑖𝑥𝛼^∗_𝑖−1,𝛼^∗_𝑖+𝜆_𝑖𝐼_𝑖𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1{},=𝛾,𝑖∈2,…,𝑚−1𝜆_𝑚−1𝐼_𝑚𝑥𝛼^∗_𝑚−1,𝛼^∗_𝑚+^𝑘_𝑗=𝑚+1𝜆_𝑗𝐼_𝑚𝑥𝛼^∗_𝑗,𝛼^∗_𝑚=𝛾,𝜆_𝑗𝐼_𝑗𝑥𝛼^∗_𝑗,𝛼^∗_𝑚=𝛾,𝑗∈{𝑚+1,…,𝑘},𝜆_𝑖𝛼^∗_𝑖𝐼_𝑖𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1+𝛼^∗_𝑖+1𝐼_𝑖+1𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1{},−𝑧=0,𝑖∈1,…,𝑚−1𝜆_𝑗𝛼^∗_𝑗𝐼_𝑗𝑥𝛼^∗_𝑗,𝛼^∗_𝑚+𝛼^∗_𝑚𝐼_𝑚𝑥𝛼^∗_𝑗,𝛼^∗_𝑚{}.−𝑧=0,𝑗∈𝑚+1,…,𝑘 (20)(21)(22)(23)(24)(25)(26)
If 𝜆₁=0, 𝛾 will be zero. Thus, we could conclude that all 𝜆_𝑖,𝜆_𝑗=0, 𝑗=𝑚+1,…,𝑘, by putting 𝛾=0 into (21) to (24). Substituting 𝜆₁=0 and 𝛾=0 into (22) and (23) will result in 𝜆_𝑖=0, 𝑖=2,…,𝑚−1. However, this contradicts with (20) which requires at least one 𝜆_𝑖,𝜆_𝑗>0. Thus, we conclude that 𝜆₁>0 and 𝛾>0 because 𝐼_𝑖(𝑥(𝛼^∗_𝑖,𝛼^∗_𝑖+1)), 𝑖=1,…,𝑘−1, is strictly positive. Similarly, we could conclude that 𝜆_𝑗>0, 𝑗=𝑚+1,…,𝑘, from (24), and max{𝜆_𝑖,𝜆_𝑖+1}>0, 𝑖=1,…,𝑚−1, from (22).
Based on the results that 𝜆₁>0, 𝜆_𝑗>0, 𝑗∈{𝑚+1,…,𝑘}, max{𝜆_𝑖,𝜆_𝑖+1}>0, 𝑖=2,…,𝑚−1, and constraints of (18), we have the following equality from the complementary slackness condition in (25) and (26). For 𝑖∈{2,…,𝑚−1}, 𝑗∈{𝑚+1,…,𝑘}, 𝑧=𝛼^∗₁𝐼₁𝑥𝛼^∗₁,𝛼^∗₂+𝛼^∗₂𝐼₂𝑥𝛼^∗₁,𝛼^∗₂=𝛼^∗_𝑗𝐼_𝑗𝑥𝛼^∗_𝑗,𝛼^∗_𝑚+𝛼^∗_𝑚𝐼_𝑚𝑥𝛼^∗_𝑗,𝛼^∗_𝑚=min𝛼^∗_𝑖𝐼_𝑖𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1+𝛼^∗_𝑖+1𝐼_𝑖+1𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1,𝛼^∗_𝑖−1𝐼_𝑖−1𝑥𝛼^∗_𝑖−1,𝛼^∗_𝑖+𝛼^∗_𝑖𝐼_𝑖𝑥𝛼^∗_𝑖−1,𝛼^∗_𝑖. (27)
So we have proved the assertion that a point satisfying the KKT conditions of (18) must be feasible to (19). This completes the proof of Theorem 4.

Therefore, we could conclude that the optimal allocation rule to rank the top 𝑚 designs 𝜶^∗=(𝛼^∗₁,…,𝛼^∗_𝑘) solves ^𝑘_𝑖=1𝛼_𝑖=1,𝛼_𝑖𝛼>0,min^∗_𝑖𝐼_𝑖𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1+𝛼^∗_𝑖+1𝐼_𝑖+1𝑥𝛼^∗_𝑖,𝛼^∗_𝑖+1,𝛼^∗_𝑖−1𝐼_𝑖−1𝑥𝛼^∗_𝑖−1,𝛼^∗_𝑖+𝛼^∗_𝑖𝐼_𝑖𝑥𝛼^∗_𝑖−1,𝛼^∗_𝑖=𝛼^∗₁𝐼₁𝑥𝛼^∗₁,𝛼^∗₂+𝛼^∗₂𝐼₂𝑥𝛼^∗₁,𝛼^∗₂=𝛼^∗_𝑗𝐼_𝑗𝑥𝛼^∗_𝑗,𝛼^∗_𝑚+𝛼^∗_𝑚𝐼_𝑚𝑥𝛼^∗_𝑗,𝛼^∗_𝑚,𝑖=2,…,𝑚−1;𝑗=𝑚+1,…,𝑘. (28)

Suppose that the performance of each design follows the normal distribution; that is, 𝑋_𝑖∼𝑁(𝜇_𝑖,𝜎²_𝑖), 𝑖=1,…,𝑘. Equation (28) can be rewritten as follows: ^𝑘_𝑖=1𝛼_𝑖=1,𝛼_𝑖𝜇>0,min^𝑖−𝜇_𝑖+1²2𝜎²_𝑖/𝛼^∗_𝑖+𝜎²_𝑖+1/𝛼^∗_𝑖+1,𝜇_𝑖−𝜇_𝑖−1²2𝜎²_𝑖/𝛼^∗_𝑖+𝜎²_𝑖−1/𝛼^∗_𝑖−1=(𝜇₁−𝜇₂)²2𝜎²₁/𝛼^∗₁+𝜎²₂/𝛼^∗₂=(𝜇_𝑗−𝜇_𝑚)²2𝜎²_𝑗/𝛼^∗_𝑗+𝜎²_𝑚/𝛼^∗_𝑚,𝑖=2,…,𝑚−1;𝑗=𝑚+1,…,𝑘. (29)

5. Sequential Allocation

Based on the results from Theorem 4, we can compute the simulation budget allocation rule using (28) for design performances with any arbitrary distributions or using (29) for normally distributed performances if the parameters of the design performance distribution are given. However, no information on the design performance distribution is known before simulation experiments are conducted. To overcome this dilemma, we suggest a heuristic sequential allocation algorithm in order to implement the allocation rule. The algorithm for sequential allocation is summarized in Algorithm 1.

INPUT: 𝑚,𝑘,𝑛,𝑛0,Δ (𝑛−𝑘𝑛0 is a multiple of Δ and 𝑛0≥5) INITIALIZE: Perform 𝑛0 simulation replications for all designs.      𝑙←0;𝑁𝑙1=⋯=𝑁𝑙𝑖=⋯=𝑁𝑙𝑘 WHILE ∑𝑘𝑖=1𝑁𝑙𝑖<𝑛 DO    (1) Increase the computing budget by Δ    (2) Use (28) to determine the new budget allocation rule, that is, 𝑁𝑖𝑙+1,𝑖=1,…,𝑘.       The population distribution parameters are estimated by the sample statistics.    (3) Simulate additional max(0,𝑁𝑖𝑙+1−𝑁𝑙𝑖) for each design 𝑖=1,…,𝑘. END OF WHILE

Algorithm 1 

Sequential allocation algorithm.

Define 𝑙 to be the iteration number and define 𝑁^𝑙_𝑖, 𝑖=1,…,𝑘 to be the total number of simulation replications that have been allocated to design 𝑖 up to iteration 𝑙. 𝑛 is the total number of simulation replications available. Δ is the number of incremental simulation replications for each iteration. 𝑛₀ is the initial number of simulation runs for each design.

As the simulation continues, design 𝑖 will be ranked number 𝑖 for all 𝑖≤𝑚. The ranking of the top 𝑚 designs may change from iteration to iteration although it will converge to the true ranking when the total computating budget goes to infinity. When the ranking of the top 𝑚 designs changes, the budget allocation in the loop will be applied immediately. Therefore, the actual proportion of budget for every system will converge to the optimal proportion when the number of iterations is sufficiently large.

Furthermore, we need to take note of 𝑛₀, the initial number of replications for every design. 𝑛₀ cannot be too small because the estimation of the rate function can be poor especially when the variance of the performance is large. On the other hand, if 𝑛₀ is too large, some designs will be allocated excessively compared with their optimal allocation numbers. When the total budget is very limited, designs that need more replications may suffer from large 𝑛₀ and this would eventually affect the simulation results. Other than the initial number of replications, the incremental budget Δ is also important in the implementation process. Large Δ results in wasting of budget, while small Δ will lead to expensive computation in the loop.

It is also worthy to note that there will be significant variation in the estimate of the performance value. Denote the empirical cumulant generating function and rate function of system 𝑖 as Λ_𝑖^(𝑛)(𝜃)=ln((1/𝑛)∑^𝑛_𝑗=1𝑒^{𝜃𝑋𝑖𝑗}) and 𝐼_𝑖^(𝑛)(𝑥)=sup_𝜃{𝜃𝑥−Λ_𝑖^(𝑛)(𝜃)}, respectively. In Algorithm 1, we are using 𝐺_𝑖^(𝑛)(𝛼_𝑖,𝛼_𝑖+1) to estimate 𝐺_𝑖(𝛼_𝑖,𝛼_𝑖+1) for the optimal allocation rule. Thus, the true optimal allocation 𝛼^∗_𝑖 is actually estimated by 𝛼^∗_𝑖(𝑛). As argued in [19], the empirical estimation of the optimal allocation is consistent; that is, 𝛼^∗_𝑖(𝑛)→𝛼^∗_𝑖, for all 𝑖=1,…,𝑘 almost surely when 𝑛→+∞.

6. Numerical Experiments

In this section, we test the proposed simulation budget allocation rule for ranking the top 𝑚 designs by comparing it with different allocation rules: equal allocation which simulates each design equally and the OCBA-𝑚 procedure [12]. Although OCBA-𝑚 only considers the selection of the top 𝑚 designs and do not aim to identify the ranking of the top 𝑚 designs, we can use it here for benchmarking purpose. In all the experiments below, the performance of each design is assumed to follow the normal distribution. Therefore, the optimal allocation rule can be obtained by solving (29). The assumption of normal distribution is generally held in simulation experiments since the output is obtained from an average performance or batch means, so that the central limit theorem effect holds.

6.1. Computing Budget Allocation Rules

Equal Allocation. The simulation replications are allocated equally to each design; that is, 𝛼_𝑖=1/𝑘, 𝑖=1,…,𝑘. This is the simplest allocation rule and it can serve as a benchmark for other allocation procedures.

Top 𝑚 Ranking (OCBA-Rm). The simulation budget allocation is derived using (29) and implemented using the sequential allocation algorithm proposed in Section 5. The allocation rule is solved by using the solver “fminimax” in Matlab.

OCBA-m [12]. Sequential allocation algorithm is used to implement this allocation rule. The simulation budget allocation rule is determined by the following equation: 𝛼₁^𝑙+1𝑠²₁/𝛿²₁=𝛼₂^𝑙+1𝑠²₂/𝛿²₂𝛼=⋯=^𝑘𝑙+1𝑠²_𝑘/𝛿²_𝑘, (30) where 𝑠²_𝑖 is the sample variance for design 𝑖, 𝛿²_𝑖=(𝑋_𝑖−𝑐)², and 𝑐=(𝑠_𝑚+1𝑋_𝑚+𝑠_𝑚𝑋_𝑚+1)/(𝑠_𝑚+1+𝑠_𝑚) as defined in [12].

6.2. Numerical Results for Different Allocation Procedures

To compare the performance of the procedures, we carried out numerical experiments for the different allocation procedures discussed above. In comparing the procedures, the effectiveness of the procedures is measured by the probability of correctly ranking the top 𝑚 designs (𝑃(CR_𝑚)) which is estimated by the fraction of the times that the procedure successfully identifies the correct ranking of the top 𝑚 designs out of 10,000 independent simulation runs.

Each of the procedures simulates each of the 𝑘 designs for 𝑛₀=20 replications initially as recommended in [1, 12]. The simulation budget is increased by Δ=40 for each iteration. Table 1 summarizes the mean and variance for the three experiments that will be conducted in this section. In this experiment, the total number of designs is 20, that is, 𝑘=20 and the objective is to identify the ranking of the top 5 designs, that is, 𝑚=5.

Table 1 

Parameters for the numerical experiments.

Figure 1 shows the numerical comparison for the three allocation rules, where (a) is for the equal spacing scenario, (b) is for the equal variance scenario, and (c) is for the increasing spacing but decreasing variance scenario.

(a)

(b)

(c)

Figure 1 

Probability of correctly ranking the top 𝑚 designs under three different allocation rules.

The experiment results show that the proposed simulation budget allocation rule OCBA-Rm performs the best in all three experiments. It is also interesting to note that OCBA-𝑚, which performs significantly better than EA when the objective is selecting the top 𝑚 designs, fares much worse than EA in all three experiments when the objective is to identify the ranking of the top 𝑚 designs. This is because OCBA-𝑚 only focuses on distinguishing the set of the top 𝑚 designs and the set of the nontop 𝑚 designs without considering the ranking within the top 𝑚 designs. Based on the numerical results, we can see that it is important to derive the OCBA-Rm when the ranking of the top 𝑚 alternatives is needed. By using the proposed OCBA-Rm allocation rule, significant number of simulation budgets can be saved compared with EA and OCBA-𝑚.

7. Conclusion

In this paper, we study the problem of simulation budget allocation of ranking the top 𝑚 designs out of 𝑘 alternatives. Based on the large deviation theory, we have derived the asymptotically optimal allocation rule to maximize the probability of correctly ranking the top 𝑚 designs. In addition, a heuristic sequential allocation algorithm is suggested to implement the simulation budget allocation rule. Numerical experiments are conducted to compare the effectiveness of the proposed simulation budget allocation with some existing allocation procedures. The contribution of this paper is twofold. From a ranking and selection perspective, we offer a heuristic for ranking the top 𝑚 designs out of 𝑘 alternatives, where our empirical studies show that it can be more efficient than the existing methods. From the computing budget allocation perspective, our heuristic illustrates how the previous OCBA method for identifying the single best design can be modified to identify the ranking of the top 𝑚 designs. In particular, we derive the asymptotically optimal allocation using large deviation theory, which allows us to remove the assumption that the performance of each design must be normally distributed.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the editor and the anonymous reviewers for the constructive comments which have improved the quality of this paper.