Optimizing resource allocation in service systems via simulation: A Bayesian formulation

Abstract

The service sector has become increasingly important in today's economy. To meet the rising expectation of high‐quality services, efficiently allocating resources is vital for service systems to balance service qualities with costs. In particular, this paper focuses on a class of resource allocation problems where the service‐level objective and constraints are in the form of probabilistic measures. Further, process complexity and system dynamics in service systems often render their performance evaluation and optimization challenging and relying on simulation models. To this end, we propose a generalized resource allocation model with probabilistic measures, and subsequently, develop an optimal computing budget allocation (OCBA) formulation to select the optimal solution subject to random noises in simulation. The OCBA formulation minimizes the expected opportunity cost that penalizes based on the quality of the selected solution. Further, the formulation takes a Bayesian approach to consider the prior knowledge and potential performance correlations on candidate solutions. Then, the asymptotic optimality conditions of the formulation are derived, and an iterative algorithm is developed accordingly. Numerical experiments and a case study inspired by a real‐world problem in a hospital emergency department demonstrate the effectiveness of the proposed algorithm for solving the resource allocation problem via simulation.

INTRODUCTION

The past few decades have witnessed a rapid growth of the service economy globally (Buckley & Majumdar, 2018). Nowadays, the service sector plays a vital role in a wide range of industries, including retail, healthcare, hospitality, finance, and IT. Traditional industries such as manufacturing have also seen the trend of servitization (Örsdemir et al., 2019) by introducing more service components in product offerings, requiring manufacturing processes to be more dynamic and flexible (Dmitry et al., 2021). Such service‐oriented systems are often tailored to meet different customer requirements and expectations, and are complicated in process dependencies and dynamics. Further, many activities in service systems deal directly with humans rather than machines, thereby introducing various forms of variability and challenges in measuring system performance. This paper focuses on service systems without tractable mathematical structures in view of the operational complexity, process nonlinearity, and performance variability, thus relying on simulation models in system evaluation and optimization.

A notable challenge of the service systems under study is to meet rising expectations on services without overspending on relevant resources. On the one hand, service providers are required or incentivized to meet and improve service quality as measured by certain metrics, such as waiting time (WT) in call centers and hospitals, on‐time ratio in online food ordering and delivery platforms, and availability percentage in IT services. These service goals are expressly provisioned in a service contract or perceived as crucial to customer satisfaction and market share. On the other hand, service providers work with a given set of resources, such as personnel with specific skills or training, space and equipment availability, and an operating budget. With the limitation on resources, they are motivated to improve service by intelligently allocating resources in place (e.g., functional areas) and time (e.g., shifts). For example, a hospital can improve care delivery by managing bed assignments and nurse schedules (Best et al., 2015; Kim & Mehrotra, 2015); a retail store can optimize product variety and quantity by properly utilizing storage and shelf spaces (Ton & Raman, 2010). This paper studies a class of service problems where one service level is optimized, and several others are imposed as constraints, subject to given resource limitations.

A case in point is a real‐world resource allocation problem in a hospital emergency department (ED) (see Section 3.1 for the formulation and Section 5.2 for the case study). Common resources in an ED include specialized medical staff, beds, medical equipment, and salary budgets. The ED has established a set of achievable service goals as measured by WTs for different categories of patients. While the target service levels on urgent patients are high and satisfiable, the ED has recognized the need to improve the service quality on less‐urgent patients. Therefore, it is desirable to seek resource allocation scenarios that improve the “nice‐to‐have” service goal while continuing to meet all other “must‐have” requirements. Furthermore, the imposed service levels on WTs are expressed in the form of probabilistic measures; that is, a service level is defined as the percentage of customers not waiting longer than a given threshold.

Such probabilistic constraints (or chance constraints) are often used in measuring quality of service (Hong et al., 2015), and thus problems with similar structures appear in other service systems. Call centers may operate under a delay‐percentile contract, where delays for a certain percentage of contract customers should be within a fixed delay bound, and noncontract customers have no guarantees on delay (Milner & Olsen, 2008). A call center may seek an operator schedule to minimize the delay for noncontract customers, while satisfying the requirement for contract customers. Such problems are also prevailing in public service offices, such as governmental departments that issue passports, driver licenses, and ID cards. They need to determine the staffing level for a specific time slot, such that at least a certain percentage of customers would wait less than a given time limit (Taigel et al., 2018). In this paper, we will use the term “resource allocation” for the aforesaid staffing schedules, equipment assignments, and other allocation of resources in place and time.

To study these service systems with probabilistic measures, simulation models can be deployed to evaluate the performance of different resource allocation scenarios so as to identify the best one. However, one notable drawback of performing analysis using simulation is the existence of random noises and thus the need to run a large number of simulation replications in order to evaluate each scenario accurately. With sizable alternatives, it will be very time consuming to run sufficient replications for every scenario, and thus an efficient strategy is necessary to search for the optimal one subject to the computing budget (i.e., the total number of simulation replications). This line of research falls under the domain of ranking‐and‐selection (R&S). To this end, we will develop an optimal computing budget allocation (OCBA) formulation using the expected opportunity cost (EOC). Further, the proposed formulation will follow the Bayesian framework and impose a correlated prior belief, as treated in Frazier et al. (2009), on the performance of resource allocation scenarios. As such, it utilizes prior knowledge accumulated by the service provider on these scenarios and captures the performance correlation on “neighboring” scenarios.

The remainder of this paper is organized as follows. Section 2 summarizes related literature on applications of stochastic service systems and methodologies of simulation optimization. Section 3 introduces the resource allocation model and the OCBA formulation, with the solution methodology developed in Section 4. Section 5 presents the numerical results, and Section 6 concludes this paper. Proofs of lemmas and theorems, as well as additional analysis and numerical results, are enclosed in the Supporting Information.

LITERATURE REVIEW

In this section, we will review relevant literature from both application and methodology viewpoints. As the problem under study considers probabilistic measures, we limit the review to those dealing with stochastic systems.

From the application point of view, most modeling papers studying stochastic resource allocation in service systems focus on problems with tractable mathematical structures. These well‐structured problems can be modeled and solved as stochastic optimization problems, such as stochastic programming (Bodur & Luedtke, 2017), robust optimization (Mattia et al., 2017), and Markov decision process (Huh et al., 2013). These problems typically have one or more stochastic constraints and/or objective with others being deterministic. While the stochastic constraints/objective are often in the form of expectations (see, e.g., Bodur & Luedtke, 2017; Gans et al., 2015), some problems study probabilistic constraints as those formulated in this paper. For example, Beraldi et al. (2004) developed a stochastic programming framework to determine the emergency medical service site locations and the number of emergency vehicles, where the service‐level constraints are probabilistic constraints; Robbins and Harrison (2010) studied a call center scheduling problem subject to a service constraint that a proportion of calls must be answered within a fixed time.

For stochastic service systems that do not possess tractable structures, simulation has been the tool of choice for analysis. Simulation can be used to evaluate the stochastic feasibility requirements, such as passenger processing targets at an airport (Mason et al., 1998), and average WT requirements for patients in EDs (Izady & Worthington, 2012). Simulation can also be combined into optimization procedures to search for the optimal solution, known as simulation(‐based) optimization, or optimization via simulation. For example, Atlason et al. (2008) and Cezik and L'Ecuyer (2008) solved call center staffing problems with non–closed‐form service‐level functions evaluated via simulation; Tsai and Zheng (2013) addressed a two‐echelon inventory problem, where one constraint requires the expected response time at each depot to be within a threshold; Guo et al. (2017) and Chen et al. (2020) studied staffing problems in hospital EDs while satisfying stochastic service quality requirements. This paper contributes to the above literature on studying a generalized resource allocation problem in service systems, where two unique structures exist: (1) The performance of neighboring resource allocation scenarios may be correlated and the modeler can obtain prior knowledge on their joint distribution; and (2) the objective and multiple constraints are stochastic and in the form of probabilistic measures.

To tackle these special structures, this paper contributes to the R&S literature from the methodology point of view, particularly to those constrained R&S methods for solving problems with stochastic constraints and/or objective. One stream of research on constrained R&S is based on the sequential indifference‐zone (IZ) framework and aims to identify the optimal solution with a probability exceeding a stipulated threshold (Andradóttir & Kim, 2010; Batur & Kim, 2010; Healey et al., 2014; Hong et al., 2015). Another type of formulation is based on the OCBA framework and aims to maximize the probability of correct selection (PCS) (Lee et al., 2012). A heuristic procedure was developed in Choi et al. (2021) by considering the precision of the sample means when maximizing PCS, so as to improve the algorithm performance when large simulation noise exists. Finally, large‐deviation approaches have also been investigated, where Hunter and Pasupathy (2013) and Pasupathy et al. (2015) seek to maximize the rate of decay of the probability of false selection (PFS). Gao et al. (2019) further incorporated quadratic regression metamodels in the large‐deviation framework to improve the search efficiency.

While these existing methods optimize the PCS or PFS measure, this paper develops a new Bayesian OCBA formulation using the EOC measure, motivated by the aforesaid structures of the resource allocation problem under study. First, the OCBA formulations with the PCS or PFS measure aim to identify the best solution with the maximum probability, treating all nonbest solutions equally. However, for the service systems under study, the industry is more interested in the economic value of the selection than its statistical significance. In other words, the decision makers are more concerned with the loss incurred in selecting a solution, rather than the statistical significance of selecting exactly the best one. As such, the OCBA formulation proposed in this paper uses the EOC measure as the objective, which penalizes for the quality of the selection. While EOC has been used as an objective in unconstrained R&S problems (Chick & Wu, 2005; Gao et al., 2017), it has rarely been studied in the constrained R&S counterparts. Second, in order to capture the prior knowledge on the solution performances and the correlation of neighboring solutions, the proposed formulation uses the Bayesian framework. While Pujowidianto et al. (2013) proposed a frequentist version of EOC for a constrained R&S problem, to the best of our knowledge, this is the first Bayesian EOC formulation developed in the constrained R&S literature that accounts for the prior belief and correlation of solutions. Third, the probabilistic measures in the objective and constraints of the problem provide convenience in the formulation and its theoretical properties, including a strong finite‐time convergence property of the proposed algorithm. We mention that although the probabilistic measures can be written as expectations of an indicator function following a Bernoulli distribution, the information on the Bernoulli distribution makes our method more efficient and robust, as pointed out in Hong et al. (2015).

Finally, this paper is more tangentially related to the literature on stochastically constrained optimization via simulation (COvS) and stochastically constrained Bayesian optimization (CBO). Both COvS and CBO are constrained optimization problems where samples from the objective and constraint measures can be treated as being generated from a black box. This black box is the simulation model for COvS and is the real system for CBO. In COvS and CBO, the set of candidate solutions is much larger, and search‐based heuristic methods need to be developed to find the optimal or near‐optimal solutions. Such heuristics can be based on penalty functions (Park & Kim, 2015), gradients (Luo & Lim, 2013; Nagaraj & Pasupathy, 2013), random search strategies (Chen et al., 2020; Gao & Chen, 2016), or acquisition functions (Gardner et al., 2014; Letham et al., 2019; Ungredda & Branke, 2021). A valuable future research is to extend the methodology developed in this paper to solve the COvS and CBO problems.

PROBLEM DESCRIPTION

In this section, we first introduce a resource allocation problem in a hospital, and then generalize the model to other service systems. Then, a Bayesian OCBA model is formulated to efficiently identify the optimal solution of the model via simulation.

A resource allocation model with probabilistic measures

A motivating example of our model is a resource allocation problem in the ED of a large public hospital in Hong Kong. This hospital had large patient flows with an average of 350 to 430 patients per day to the ED. Sixty‐five percent of the patients were walk‐in patients, while the others were brought in by ambulance. All patients were labeled by triage categories based on their clinical conditions, including categories I (critical), II (emergency), III (urgent), IV (semiurgent), and V (nonurgent), of which categories III and IV patients made up the major proportions (about 90%). The hospital committed to provide high‐quality services to ED patients, where WT serves as a key metric. An example of such service‐level requirements measured by the patient WT and length of stay (LOS) is shown in Table 1.

TABLE 1

Time limits and service levels for an ED in Hong Kong

Patient category	Expected time limit	Desired service level
I	Immediately	Whenever possible
II	WT ⩽ 15 min	≥ 95%
III	WT ⩽ 30 min	≥ 90%
IV	WT ⩽ 2 h	≥ 75%
I–V	LOS ⩽ 4 h	≥ 98%

As seen from Table 1, patients in category I should be treated immediately with any available resources. The WT and service‐level expectations for patients in categories II, III, and IV lower as the urgency drops, while no requirement is imposed on category V patients alone. The ED also maintains a goal of treating at least 98% of patients in all categories within 4 h, measured by LOS. While these requirements were established based on benchmarks and attainability, the focus was mainly on urgent patients. The hospital leadership believed that, with an efficient allocation of resources, it is possible to improve the service level for category IV patients (semiurgent and 48.76% of total) while meeting all other service requirements. Indeed, different categories of patients undergo different treatment procedures, and occupy different medical resources such as manpower, machines, and spaces. By properly assigning medical staff to each service area and each shift, the ED can find a scenario to improve the WT of category IV patients while maintaining the service for urgent patients. More specifically, define

$W{T}_2$

,

$W{T}_3$

,

$W{T}_4$

, and

$LOS$

as the daily average patient WT for categories II, III, IV, and the daily average LOS for all patients, respectively. Then, the aforementioned problem can be modeled as follows.

1

2

3

4

5

where ω represents random noises caused by uncertainties in patient arrivals and service durations, and

$\mathbf{x}$

is a vector of decision variables defining a staff assignment scenario that satisfies all the deterministic (nonprobabilistic) constraints defined by

$\mathcal{X}$

. For the above example, such deterministic constraints include the total budget constraint of personnel, and the minimum and maximum number of specialists required at each service area during each shift (e.g., no more than two nurses but no less than one nurse at the admission desk). Note that although the service level of category I patients is not explicitly modeled in the above formulation, any incoming category I patient should be treated immediately with available resources, thereby reducing the resource availability to other patients and their attendant service levels.

The aforementioned decision problem can be generalized to other service systems, such as call centers and public service offices. To this end, we define a stochastic performance measure

$F_h(\mathbf{x},\omega )$

for

$h=0,1,\text{…},H$

, where H is the total number of constraints and 0 indicates the objective function. Specifically, each

$F_h(\mathbf{x},\omega )$

is a stochastic metric given a resource allocation scenario

$\mathbf{x}$

. For example,

$F_h(\mathbf{x},\omega )$

may capture the average daily WT or peak daily LOS in the ED as required, where the value of this metric observed varies due to random noise ω. Correspondingly, each performance goal is defined as

$d_h$

for

$h=0,1,\text{…},H$

, and the desired service level as

$1-{\alpha }_h$

for

$h=1,\text{…},H$

. Without loss of generality, we assume that (1) the lower

$d_h$

is, the better performance it has; and (2) the higher

$1-{\alpha }_h$

is, the better service level the system provides. Hence, the generalized resource allocation problem for service systems can be formulated as follows:

$$(\mathcal{P})\underset{\mathbf{x}}{\min }P\{F_0(\mathbf{x},\omega )>d_0\}$$

6

7

8

where

$\mathcal{X}$

is a finite discrete solution space of

$\mathbf{x}$

, determined by the deterministic constraints of the service system, such as budget and capacity constraints. In the remainder of this paper, we will focus on the problem

$\mathcal{P}$

defined by Equations (6)–(8) where

$\mathcal{X}$

is finite, and will refer to a resource allocation scenario

$\mathbf{x}$

as a solution of

$\mathcal{P}$

.

A simulation model will be used to evaluate the above probabilistic measures,

$P\{F_h(\mathbf{x},\omega )>d_h\}$

, for every solution

$\mathbf{x}$

. Since most simulation models are expensive and time consuming to run and the decision may need to be made within a given time period, the challenge of solving problem

$\mathcal{P}$

using a simulation model in this context is the one of allocating the limited computing budget to select the optimal solution under observation noises, thereby falling under the OCBA regime. As mentioned in Section 2, the characteristics of problem

$\mathcal{P}$

warrant a new Bayesian OCBA formulation using the EOC measure, to be introduced in the sequel.

A Bayesian OCBA formulation

Suppose there is a total of K alternative solutions (resource allocation scenarios) in

$\mathcal{X}$

, and

${\mathbf{x}}_i$

represents one such solution for

$i\in \mathrm{\Theta }=\{1,\text{…},K\}$

. For a given

${\mathbf{x}}_i$

, each simulation replication j provides one estimate,

$F_h({\mathbf{x}}_i,{\omega }_{ij})$

, for each performance measure h where

${\omega }_{ij}$

is the random noise observed in the j‐th replication. Here, we assume that for any solution

${\mathbf{x}}_i$

and probabilistic measure h,

$F_h({\mathbf{x}}_i,{\omega }_{ij})$

's are identically distributed for

$j=1,2,\text{…}$

Further define an indicator function

9

It is seen that

$X_{hij}$

follows a Bernoulli distribution with a success probability

$p_{hi}\doteq P\{F_h({\mathbf{x}}_i,{\omega }_{ij})>d_h\}=\mathbb{E}[X_{hij}]$

. Apparently,

$p_{hi}$

is unknown in practice, and can be best estimated by its sample average

${\overline{p}}_{hi}\doteq \frac{1}{n}{\sum }_{j=1}^n{X}_{hij}$

. To solve problem

$\mathcal{P}$

via simulation given the computing budget T, one needs to determine the number of replications allocated to evaluate each solution

${\mathbf{x}}_i$

(denoted by

$N_i$

), so as to best identify the true optimal solution

${\mathbf{x}}_b$

where

$b=\underset{i\in \mathrm{\Theta }:p_{hi}\le {\alpha }_h,h=1,\text{…},H}{arg\; min}{p}_{0i}$

.

To capture the prior knowledge on solutions and correlated samples of neighboring solutions, we proceed to develop an OCBA formulation from a Bayesian perspective. Specifically, since

$p_{hi}$

is practically unknown, the algorithm starts with a prior belief over each

$p_{hi}$

and updates its posterior distribution as simulation samples are gathered. To this end, let

${\mathbf{\Psi }}_h=({\mathrm{\Psi }}_{h1},\cdots ,{\mathrm{\Psi }}_{hK})$

be a random variable following a multivariate distribution of

${\pi }_h^t({\mathit{\psi }}_h)$

, that is,

${\mathbf{\Psi }}_h\sim {\pi }_h^t({\mathit{\psi }}_h)$

, where

${\mathit{\psi }}_h=({\psi }_{h1},\cdots ,{\psi }_{hK})\in {[0,1]}^K$

is a realization of

${\mathbf{\Psi }}_h$

and

${\pi }_h^t({\mathit{\psi }}_h)$

is the Bayesian belief over

$(p_{h1},\cdots ,p_{hK})$

for service level h after t simulation replications are gathered. Thus,

${\pi }_h^0({\mathit{\psi }}_h)$

and

${\pi }_h^T({\mathit{\psi }}_h)$

represent the prior distribution before simulation and the posterior distribution after the computing budget T exhausts, respectively. Note that

${\pi }_h^0({\mathit{\psi }}_h)$

allows

${\mathrm{\Psi }}_{hi}$

and

${\mathrm{\Psi }}_{h{i}^{\prime }}$

to be correlated based on prior knowledge on

${\mathbf{x}}_i$

and

${\mathbf{x}}_{i^{\prime }}$

. For each realization

${\mathit{\psi }}_h$

of

${\mathbf{\Psi }}_h$

for

$h=0,1,\cdots ,H$

, define

$\stackrel{\sim }{b}=\underset{i\in \mathrm{\Theta }:{\psi }_{hi}\le {\alpha }_h,h=1,\text{…},H}{arg\; min}{\psi }_{0i}$

as the best solution under realization

${\mathit{\psi }}_h$

. Noting that

$(p_{h1},\text{…},p_{hK})$

is one possible realization of

${\mathbf{\Psi }}_h$

, if

${\psi }_{hi}=p_{hi}$

for

$h=0,1,\text{…},H$

and

$i=1,\cdots ,K$

, then

$\stackrel{\sim }{b}=b$

.

Here, we first make a mild assumption on the prior distribution (Russo, 2020) in Assumption 1. Assumption 1

The prior joint distributions

${\pi }_h^0({\mathit{\psi }}_h)$

(

$h=0,1,\cdots ,H$

) are continuous and uniformly bounded away from 0 and ∞, that is, there exist

$0<\underline{c}<\overline{c}<\infty$

such that

$\underline{c}\le {\pi }_h^0({\mathit{\psi }}_h)\le \overline{c}$

for all

${\mathit{\psi }}_h\in {[0,1]}^K$

and all

$h=0,1,\cdots ,H$

.

Remark 1

Assumption 1 is a mild assumption that covers a large class of distributions, such as the multivariate normal distribution (Howard, 1998) that we will use in the numerical experiments in Section 5. In particular, for any continuous function

$r(\mathit{\psi })$

satisfying

$0<\underline{r}\le r(\mathit{\psi })\le \overline{r}<\infty$

for

$\psi \in {[0,1]}^K$

, we can convert

$r(\mathit{\psi })$

into a prior

${\pi }_h^0({\mathit{\psi }}_h)=\frac{r({\mathit{\psi }}_h)}{{\int }_{{[0,1]}^K}r({\mathit{\psi }}_h^{\prime })d{\mathit{\psi }}_h^{\prime }}$

that satisfies Assumption 1. When no prior knowledge of the system is available, the modeler can use the uncorrelated uniform prior instead of the correlated prior, or estimate the prior distribution using the maximum likelihood estimation through initial samples (Rasmussen & Williams, 2006). In addition, it should be noted that the lower bound

$\underline{c}$

in Assumption 1 does not have to hold for the entire region of [0, 1] ^K for

${\mathit{\psi }}_h$

. In fact, we can relax this condition and only restrict

$\underline{c}\le {\pi }_h^0({\mathit{\psi }}_h)$

to hold for certain subregions of [0, 1] ^K depending on h and solution

${\mathbf{x}}_i$

. This will become more apparent in the proof of Theorem 1, and we only mention this fact here without stating the detailed conditions. However, since such tighter conditions would overcomplicate the elaboration of Lemma 2, we prefer the succinct statement in Assumption 1 and will stick to it throughout the paper.

Let

$N_i^{(t)}$

be the number of simulation replications allocated to solution

${\mathbf{x}}_i$

when t replications are gathered, and

$L_t({\psi }_{hi})$

be the corresponding likelihood function (of observations). Note that

$N_i^{(0)}=0$

,

$N_i^{(T)}=N_i$

, and

$\sum _{i\in \mathrm{\Theta }}{N}_i^{(t)}=t$

. Further,

$L_t({\psi }_{hi})={({\psi }_{hi})}^{{\sum }_{j=1}^{N_i^{(t)}}{X}_{hij}}{(1-{\psi }_{hi})}^{{\sum }_{j=1}^{N_i^{(t)}}(1-X_{hij})}$

. The posterior distribution can be updated according to the Bayes' rule as follows:

10

The posterior distribution in Equation (10) may not have a closed form if the prior distribution is correlated. In this case, we need numerical integration methods to solve integrals in Equation (10). In the case of high‐dimensional integrals if the total number of solutions, K, is large, Monte Carlo estimation and the importance sampling method can provide robust estimates with high accuracy.

We are now ready to introduce the OCBA formulation that minimizes the Bayesian EOC. More specifically, based on T simulation samples, the selected solution

${\mathbf{x}}_{\widehat{b}}$

is the one estimated with the minimum posterior mean in the objective function while satisfying all constraints; that is,

$\widehat{b}=\underset{i\in \mathrm{\Theta }:{\mathbb{E}}_T[{\mathrm{\Psi }}_{hi}]\le {\alpha }_h,h=1,\text{…},H}{arg\; min}{\mathbb{E}}_T[{\mathrm{\Psi }}_{0i}]$

, where

11

is the posterior mean. The Bayesian EOC is then defined as

12

where

$O{C}_{h\widehat{b}}={\varpi }_h\max ({\psi }_{h\widehat{b}}-{\alpha }_h,0)$

,

$h=1,\cdots ,H$

, and

13

Here, the weight

$0<{\varpi }_h<\infty$

for

$h=1,\text{…},H$

reflects the unit cost of violating constraint h relative to the unit cost of the objective. Specifically, the opportunity cost

$O{C}_{0\widehat{b}}$

is the penalty for the objective of the selected solution

${\mathbf{x}}_{\widehat{b}}$

, and the cost

$O{C}_{h\widehat{b}}$

penalizes the loss incurred if

${\mathbf{x}}_{\widehat{b}}$

violates the constraint on service level h. Note that when

$\underset{i\in \mathrm{\Theta }:{\psi }_{hi}\le {\alpha }_h,h=1,\text{…},H}{arg\; min}{\psi }_{0i}$

does not exist, we still recommend the solution

${\mathbf{x}}_{\widehat{b}}$

when the constrained optimization problem is infeasible. Hence, we set the attendant opportunity cost to be 1 to penalize such selection, noting that the opportunity cost of

$\max \{{\psi }_{0\widehat{b}}-{\psi }_{0\stackrel{\sim }{b}},0\}\le 1$

when

$\stackrel{\sim }{b}$

exists. A large weight

${\varpi }_h$

can be imposed to penalize against selecting an infeasible solution without affecting the asymptotic results to be derived later.

Subsequently, the Bayesian OCBA formulation of problem

$\mathcal{P}$

is defined as follows:

14

We assume that in problem

$\mathcal{P}$

, the constraint measures are not exactly equal to the constraint limit on the right‐hand side; that is,

$p_{hi}\ne {\alpha }_h$

for all constraints in any solution. This is a common assumption made in the OCBA literature, which ensures that

$EO{C}_{Bayes}$

approaches 0 as the computing budget increases.

The

$\mathcal{OCBA}\mathcal{-}\mathcal{P}$

formulation distinguishes from other OCBA models, such as the one in Lee et al. (2012) for constrained R&S, not only in the use of the EOC objective, but also in capturing the performance correlation of solutions rather than assuming independence in solution measures. Intuitively speaking, capturing such correlation enables more effective learning of solution performances from samples, since the performance of a solution is not only learned using its own samples, but also using samples from its neighboring solutions. On the other hand, it brings challenges in developing theoretical results, which will be addressed next.

SOLUTION METHODOLOGY

Although

$EO{C}_{Bayes}$

in Equation (12) is very complex to analyze (Chick et al., 2010), we will show that as

$T\to \infty$

,

$\widehat{b}$

exists almost surely, and

$EO{C}_{Bayes}$

is well defined and decays in a deterministic rate. Therefore, we will next solve

$\mathcal{OCBA}\mathcal{-}\mathcal{P}$

by optimizing the approximated rate of decay of

$EO{C}_{Bayes}$

.

Rate of decay of Bayesian EOC

Considering the definition of the selected best

$\widehat{b}$

,

$\widehat{b}$

does not always exist due to simulation noises, even if the optimal solution exists for problem

$\mathcal{P}$

. Particularly, if

$N_i$

remains fixed (

$N_i<\infty$

) for some solution

${\mathbf{x}}_i$

as

$T\to \infty$

,

$\widehat{b}$

may not exist since

$\{i\in \mathrm{\Theta }:{\mathbb{E}}_T[{\mathrm{\Psi }}_{hi}]\le {\alpha }_h,h=1,\text{…},H\}$

could be an empty set. To avoid this unfavorable situation, we make the following assumption: Assumption 2

As

$T\to \infty$

,

$N_i\to \infty$

and

${\rho }_i\doteq \underset{T\to \infty }{\lim }\frac{N_i}{T}$

exists for all

$i=1,\cdots ,K$

.

Remark 2

In fact, under Assumption 2, we will show in the sequel that

$\widehat{b}$

exists given that

$\mathcal{P}$

is feasible and

$\underset{T\to \infty }{\lim }\frac{1}{T}\log EO{C}_{Bayes}\le 0$

as

$T\to \infty$

. Meanwhile, it can be shown that if

$N_i<\infty$

for some

${\mathbf{x}}_i$

as

$T\to \infty$

, either

$\widehat{b}$

does not exist or

$\underset{T\to \infty }{\lim }\frac{1}{T}\log EO{C}_{Bayes}=0$

, which is an inferior allocation compared to the allocation under Assumption 2.

To facilitate the analysis, we first divide the set of indices for nonbest solutions into two exhaustive and mutually exclusive sets,

${\mathrm{\Theta }}_O$

and

${\mathrm{\Theta }}_F$

, such that

$\mathrm{\Theta }={\mathrm{\Theta }}_O\cup {\mathrm{\Theta }}_F\cup \{b\}$

.

${\mathrm{\Theta }}_O$

is the set of indices for nonbest solutions that are feasible, whereas

${\mathrm{\Theta }}_F$

is the set of indices for nonbest solutions with at least one constraint in Equation (7) violated. Further defining the set of constraints violated for any solution

${\mathbf{x}}_i$

as

${\mathrm{\Pi }}_i=\{h:p_{hi}>{\alpha }_h,h\in \{1,\text{…},H\}\}$

, we have

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_O\doteq \{i\in \mathrm{\Theta }:i\ne b,{\mathrm{\Pi }}_i=\varnothing \},{\mathrm{\Theta }}_F\doteq \{i\in \mathrm{\Theta }:{\mathrm{\Pi }}_i\ne \varnothing \}.\end{array}$$

15

Furthermore, for any

$i\in {\mathrm{\Theta }}_F$

, we define

$h_i^{\ast }=\underset{h\in {\mathrm{\Pi }}_i}{arg\; max}(p_{hi}-{\alpha }_h)$

as the “most violated” constraint. Lemma 1 then provides an approximate and upper bound of

$EO{C}_{Bayes}$

(called

$AEO{C}_{Bayes}$

). Lemma 1

$EO{C}_{Bayes}$

in Equation (12) is bounded from above by the following

$AEO{C}_{Bayes}$

almost surely:

16

where

$h_i^{\ast }=\underset{h\in {\mathrm{\Pi }}_i}{arg\; max}(p_{hi}-{\alpha }_h)$

for any

$i\in {\mathrm{\Theta }}_F$

.

Proof

The proof is given in the Supporting Information EC.1.1 to this paper.

$\square$

Remark 3

Although

$AEO{C}_{Bayes}$

is an upper bound of

$EO{C}_{Bayes}$

, it is constructed to be as tight as possible such that it closely approximates

$EO{C}_{Bayes}$

. The proof of Lemma 1 in the Supporting Information EC.1.1 shows why this bound is tight, where the Bonferroni inequality is applied to reach Equation (EC.4), and subspace is constructed to provide a tight bound on the EOC in Equation (EC.5) for both cases of

$i\in {\mathrm{\Theta }}_O$

and

$i\in {\mathrm{\Theta }}_F$

. Further,

$AEO{C}_{Bayes}$

becomes tighter as the total simulation budget becomes larger, and

$AEO{C}_{Bayes}\to 0$

as

$T\to \infty$

(Russo, 2020).

We now proceed to analyze the rate of decay of

$AEO{C}_{Bayes}$

. We first present Lemma 2 and Lemma 3, and then develop Theorem 1 on the rate of decay using the lemmas. Lemma 2

For any

$0\le a_1<a_2\le 1$

,

$$\begin{array}{c}\hfill \underset{T\to \infty }{\lim }\frac{1}{T}\log {\int }_{\{{\mathit{\psi }}_h:a_1\le {\psi }_{hi}\le a_2\}}{\pi }_h^T({\mathit{\psi }}_h)d{\mathit{\psi }}_h=-{\rho }_i\underset{{\psi }_{hi}\in [a_1,a_2]}{inf}{I}_{hi}({\psi }_{hi}),\end{array}$$

17

where

$$\begin{array}{c}\hfill I_{hi}({\psi }_{hi})=p_{hi}\log \frac{p_{hi}}{{\psi }_{hi}}+(1-p_{hi})\log \frac{1-p_{hi}}{1-{\psi }_{hi}}.\end{array}$$

18

Proof

The proof is given in the Supporting Information EC.1.2 to this paper.

$\square$

Lemma 3 Ganesh et al. (2004)

Consider positive sequences

$a_j(n)$

,

$j=1,\cdots ,m$

. If

$\underset{n\to \infty }{\lim }\frac{1}{n}\log a_j(n)$

exists for all j, then

$\underset{n\to \infty }{\lim }\frac{1}{n}\log ({\sum }_{j=1}^m{a}_j(n))=\underset{j\in \{1,\text{…},m\}}{\max }\{\underset{n\to \infty }{\lim }\frac{1}{n}\log a_j(n)\}$

.

Theorem 1

The rate of decay of

$AEO{C}_{Bayes}$

in Equation (16) is bounded by the following

$AEOC\mathrm{-}B$

:

19

where

$p_{0b}<c_i<p_{0i}$

for any

$i\in {\mathrm{\Theta }}_O$

, and

$I_{hi}(x)$

is given by Equation (18) for any

$h=0,1,\cdots ,H$

and

$i=1,\cdots ,K$

.

Proof

By Equation (16) and Lemma 3, we have

20

21

22

Consider the term in Equation (20) for any

$h=1,\cdots ,H$

. Since

$1\le 1+{\varpi }_h({\psi }_{hb}-{\alpha }_h)\le 1+{\varpi }_h$

for the domain

$\{{\mathit{\psi }}_h:{\psi }_{hb}>{\alpha }_h\}$

, we have

23

and

24

resulting

25

Next, consider the term in Equation (21) for any

$i\in {\mathrm{\Theta }}_O$

. For the domain

$\{{\mathit{\psi }}_0:{\psi }_{0i}<{\psi }_{0b}\}$

, there exists a

$\delta >0$

such that

$\delta \le {\psi }_{0b}-{\psi }_{0i}\le 1$

. Following the similar proof as above, we have

26

Considering the integrand on the right‐hand side of Equation (26), define a constant

$c_i$

where

$p_{0b}<c_i<p_{0i}$

. Since

$\{{\psi }_{0i}<{\psi }_{0b}\}\subset \{{\psi }_{0i}\le c_i\}\cup \{c_i\le {\psi }_{0b}\}$

, we have

27

From Lemma 2,

$\underset{T\to \infty }{\lim }\frac{1}{T}\log {\int }_{\{{\mathit{\psi }}_h:a_1\le {\psi }_{hi}\le a_2\}}{\pi }_h^T({\mathit{\psi }}_h)d{\mathit{\psi }}_h$

exists for any

$0\le a_1<a_2\le 1$

. Thus, by Equation (26), Equation (27), and Lemma 3,

28

Substituting Equations (25) and (28) into Equations (20)–(22), we have

29

Applying Lemma 2 to the right‐hand side of Equation (29), Equation (19) readily follows, since

$I_{hi}(x)$

is a convex function with the minimum achieved at

$x=p_{hi}$

.

$\square$

We mention that the only approximation used in Theorem 1 to arrive at

$AEOC\mathrm{-}B$

is Equation (28). The value of

$c_i$

should depend only on

$p_{0b}$

and

$p_{0i}$

and affects the tightness of

$AEOC\mathrm{-}B$

. We will discuss the choice of

$c_i$

in the Supporting Information EC.2 to this paper.

Optimality conditions

Since neither minimizing

$EO{C}_{Bayes}$

in problem

$\mathcal{OCBA}\mathcal{-}\mathcal{P}$

nor maximizing its rate of decay is mathematically tractable, we instead maximize the bound of the rate of decay of

$AEO{C}_{Bayes}$

, that is,

$AEOC\mathrm{-}B$

defined in the right‐hand side of Equation (19). Letting

$z=\min \{\underset{h=1,\cdots ,H}{\min }{\rho }_b{I}_{hb}({\alpha }_h),\underset{i\in {\mathrm{\Theta }}_O}{\min }{\rho }_b{I}_{0b}(c_i),\underset{i\in {\mathrm{\Theta }}_O}{\min }{\rho }_i{I}_{0i}(c_i),\underset{i\in {\mathrm{\Theta }}_F}{\min }{\rho }_i{I}_{h_i^{\ast }i}({\alpha }_{h_i^{\ast }})\}$

, the resulting problem can be expressed as

$$(\mathcal{OCBA}\mathcal{-}\mathcal{P}\mathcal{-}\mathcal{B})\underset{z,{\rho }_1,\text{…},{\rho }_K}{\max }z$$

30

$$\text{s.t.}z\le {\rho }_b\min \left\{\underset{h=1,\cdots ,H}{\min }{I}_{hb}({\alpha }_h),\underset{k\in {\mathrm{\Theta }}_O}{\min }{I}_{0b}(c_k)\right\},$$

31

$$z\le {\rho }_i{I}_{0i}(c_i),\forall i\in {\mathrm{\Theta }}_O,$$

32

$$z\le {\rho }_i{I}_{h_i^{\ast }i}({\alpha }_{h_i^{\ast }}),\forall i\in {\mathrm{\Theta }}_F,$$

33

$$\sum _{i=1}^K{\rho }_i=1,$$

34

$${\rho }_i\ge 0,\forall i\in \mathrm{\Theta }.$$

35

Problem

$\mathcal{OCBA}\mathcal{-}\mathcal{P}\mathcal{-}\mathcal{B}$

in Equations (30)–(35) is a linear optimization problem, whose optimal solution stated in Theorem 2 provides an asymptotic optimality condition for problem

$\mathcal{OCBA}\mathcal{-}\mathcal{P}$

. Theorem 2

Let

${\eta }_i$

be as follows:

36

where

$I_{hi}(·)$

is given by Equation (18). As

$T\to \infty$

,

$AEOC\mathrm{-}B$

is minimized if

$$\begin{array}{c}\hfill \frac{{\rho }_i^{\ast }}{{\rho }_j^{\ast }}=\frac{{\eta }_j}{{\eta }_i},\forall i,j\in \mathrm{\Theta },\end{array}$$

37

and thus, the asymptotic optimal sample allocation strategy is given by

$$\begin{array}{c}\hfill N_i^{\ast }=\frac{1/{\eta }_i}{{\sum }_{k=1}^{K}1/{\eta }_k}T,\forall i\in \mathrm{\Theta }.\end{array}$$

38

Proof

Since

${\rho }_i\doteq {\lim }_{T\to \infty }\frac{N_i}{T}$

exists from Assumption 2 and

${\sum }_{i=1}^K{\rho }_i=1$

with

${\rho }_i\ge 0$

, we can treat

${\rho }_i$

as continuous variables between 0 and 1. Therefore, problem

$\mathcal{OCBA}\mathcal{-}\mathcal{P}\mathcal{-}\mathcal{B}$

becomes a linear program with continuous variables. As such, the optimal solution is obtained when the right‐hand sides of all constraints (31)–(33) are equal. To see why, z can only be decreased further from this solution if all

${\rho }_i$

for

$i\in \mathrm{\Theta }$

are decreased simultaneously so as to decrease all right‐hand sides in Equations (31)–(33). However, decreasing one

${\rho }_i$

will increase another since

${\rho }_b+{\sum }_{i\in {\mathrm{\Theta }}_O}{\rho }_i+{\sum }_{i\in {\mathrm{\Theta }}_F}{\rho }_i=1$

from Equation (34). Hence, by the definition in Equation (36), we have

$$\begin{array}{c}\hfill {\rho }_i^{\ast }{\eta }_i={\rho }_j^{\ast }{\eta }_j,\forall i,j\in \mathrm{\Theta }.\end{array}$$

39

Note that

$I_{hi}({\psi }_{hi})$

in Equation (18) is a convex function with its minimum value of 0 achieved at

${\psi }_{hi}=p_{hi}$

. By the definition of problem

$\mathcal{OCBA}\mathcal{-}\mathcal{P}$

,

$p_{hi}\ne {\alpha }_h$

for all h and i, and

$p_{0b}<c_i<p_{0i}$

by choice, we have

${\eta }_i>0$

. Thus, the optimality condition in Equation (37) readily follows, so as the asymptotically optimal sample allocation strategy in Equation (38) considering the constraint (34).

$\square$

Intuitively speaking, the simulation replications allocated to solution

${\mathbf{x}}_i$

should be proportional to the inverse of its rate function

${\eta }_i$

given by Equation (38). Observe that the rate function in Equation (18) is a convex function with its minimum achieved at

${\psi }_{hi}=p_{hi}$

. The closer

${\psi }_{hi}$

is to

$p_{hi}$

, the smaller

$I_{hi}({\psi }_{hi})$

is. Therefore, the intuition is that we should increase the sample size to a solution if either optimality or feasibility is harder to detect. More specifically, the sample size allocated to solution

${\mathbf{x}}_i$

should become larger as

$c_i$

gets closer to

$p_{0i}$

for any

$i\in {\mathrm{\Theta }}_O$

, and as the constraint violation becomes less substantial for any

$i\in {\mathrm{\Theta }}_F$

; for the best solution

${\mathbf{x}}_b$

, its sample size should be increased if one of the attendant service levels

$p_{hb}$

is close to the required service level

${\alpha }_h$

, or if one of the

$c_i$

gets close to

$p_{0b}$

. This intuition is consistent with the standard OCBA allocation strategy for normally distributed random variables in Chen et al. (2000), where more samples are allocated to solutions whose optimality is harder to decide, or whose variance is larger.

Theorem 2 provides a simple closed‐form formula for allocating the total computing budget to all solutions, such that the

$AEOC\mathrm{-}B$

is asymptotically optimized. That is, such an allocation strategy optimizes the Bayesian EOC of selecting the estimated best solution

$\widehat{b}$

given a sufficiently large computing budget. The allocation strategy in Theorem 2 may not be optimal when T is relatively small. However, the performance shall get better as T becomes larger. Remark 4

It may be of interest to consider two degenerate cases of problem

$\mathcal{P}$

. The first case is an unconstrained version of problem

$\mathcal{P}$

, where stochastic constraints in Equation (7) no longer exist. In this case, Theorem 2 still holds, except that

$H=0$

,

${\mathrm{\Theta }}_F=\varnothing$

(no infeasible solution for the unconstrained problem), and

${\mathrm{\Theta }}_O=\mathrm{\Theta }\setminus \{b\}$

, resulting in a degenerate case of Equation (36) given by

40

The second case is a feasibility detection variant of problem

$\mathcal{P}$

, where the objective function in Equation (6) no longer exists. The goal is to find all feasible solutions meeting the service targets, that is,

41

Let

${\mathcal{S}}_f=\{i\in \mathrm{\Theta }:{\mathbf{x}}_i\in \mathcal{X}\text{and}P\{F_h({\mathbf{x}}_i,\omega )>d_h\}\le {\alpha }_h\text{for}h=1,\text{…},H\}$

represent the set of true feasible solutions and let

${\widehat{\mathcal{S}}}_f=\{i\in \mathrm{\Theta }:{\mathbf{x}}_i\in \mathcal{X}\text{and}{\mathbb{E}}_T[{\mathrm{\Psi }}_{hi}]\le {\alpha }_h\text{for}h=1,\text{…},H\}$

represent the set of feasible solutions estimated from samples. The Bayesian EOC for problem

${\mathcal{P}}_f$

can be defined as

42

where

$O{C}_{hi}^{f_1}={\varpi }_h\max \{{\psi }_{hi}-{\alpha }_h,0\}$

and

$O{C}_{hi}^{f_2}={\varpi }_h\max \{{\alpha }_h-{\psi }_{hi},0\}$

,

$h=1,\cdots ,H$

. Using analysis similar to Theorem 1, we can show that the rate of decay of

$EO{C}_{f:Bayes}$

is bounded by

$$\begin{array}{c}\hfill \min \left\{\underset{i\in {\mathcal{S}}_f}{\min }{\rho }_i\underset{h=1,\cdots ,H}{\min }{I}_{hi}({\alpha }_h),\underset{i\notin {\mathcal{S}}_f}{\min }{\rho }_i{I}_{h_i^{\ast }i}({\alpha }_{h_i^{\ast }})\right\}.\end{array}$$

43

Maximizing Equation (43), the asymptotic optimal sample allocation strategy has the same form as Equation (38), where

${\eta }_i$

in Equation (36) becomes

44

Iterative algorithm and convergence analysis

Recall that

${\eta }_i$

in Theorem 2 is computed using the true success probability of Bernoulli variables,

$p_{hi}$

. However,

$p_{hi}$

is practically unknown, but can be estimated by the posterior mean. Therefore, we propose an iterative algorithm in Algorithm 1 to sequentially allocate simulation replications to the solutions in consideration, by implementing Theorem 2.

ALGORITHM 1

Iterative Sample Allocation Algorithm for Solving

$\mathcal{OCBA}\mathcal{-}\mathcal{P}\mathcal{-}\mathcal{B}$

0.	Initialize: Set iteration counter $r=0$ and the incremental budget Δ. Set initial sample size $N_i^r=n_0$ and incremental sample size ${\delta }_i^r=N_i^r$ for all $i\in \mathrm{\Theta }$ . Specify the prior ${\pi }_h^0({\mathit{\psi }}_h)$ for $h=0,1,\cdots ,H$ .
1.	Simulate: Run simulation for each solution ${\mathbf{x}}_i$ , $i\in \mathrm{\Theta }$ , for ${\delta }_i^r$ replications. For the j‐th simulation replication of ${\mathbf{x}}_i$ , record the simulation output $F_h({\mathbf{x}}_i,{\omega }_{ij})$ , and compute $X_{hij}$ using Equation (9) for all h.
2.	Estimate: Update the posterior distributions ${\pi }_h^t({\mathit{\psi }}_h)$ using Equation (10) for all h, where $t={\sum }_{i=1}^K{N}_i^r$ . Estimate the success probability of service level h for solution ${\mathbf{x}}_i$ as ${\widehat{p}}_{hi}={\mathbb{E}}_t[{\mathrm{\Psi }}_{hi}]$ using Equation (11) for all h and $i\in \mathrm{\Theta }$ .
3.	Divide: Compute the set of violated constraints for any $i\in \mathrm{\Theta }$ as ${\widehat{\mathrm{\Pi }}}_i=\{h\in \{1,\text{…},H\}:{\widehat{p}}_{hi}>{\alpha }_h\}$ . Then, find the sample best solution: $$\begin{array}{c}\hfill \widehat{b}=\underset{i\in \mathrm{\Theta }:{\widehat{\mathrm{\Pi }}}_i=\varnothing }{arg\; min}{\widehat{p}}_{0i},\end{array}$$ 45 and divide the non‐best solutions into two sets: $$\begin{array}{c}\hfill {\widehat{\mathrm{\Theta }}}_O=\{i|i\ne \widehat{b},{\widehat{\mathrm{\Pi }}}_i=\varnothing \},{\widehat{\mathrm{\Theta }}}_F=\{i|{\widehat{\mathrm{\Pi }}}_i\ne \varnothing \}.\end{array}$$ 46
4.	Stop: If $t\ge T$ , terminate and output $\widehat{b}$ as the best solution; otherwise, continue.
5.	Update: For each $i\in {\widehat{\mathrm{\Theta }}}_O$ , choose ${\widehat{c}}_i=w_i{\widehat{p}}_{0\widehat{b}}+(1-w_i){\widehat{p}}_{0i}$ where $0<w_i<1$ with $w_i=\frac{1}{2}$ recommended; and for each $i\in {\widehat{\mathrm{\Theta }}}_F$ , choose ${\widehat{h}}_i^{\ast }=\underset{h\in {\mathrm{\Pi }}_i}{arg\; max}({\widehat{p}}_{hi}-{\alpha }_h)$ . Then, for each $i\in \mathrm{\Theta }$ , compute ${\widehat{\eta }}_i^r$ using Equation (36), where ${\mathrm{\Theta }}_F$ and ${\mathrm{\Theta }}_O$ are respectively approximated by ${\widehat{\mathrm{\Theta }}}_F$ and ${\widehat{\mathrm{\Theta }}}_O$ , and $I_{hi}(·)$ is computed using Equation (18) with $p_{hi}$ approximated by ${\widehat{p}}_{hi}$ . If ${\widehat{\eta }}_i^r=0$ , bound ${\widehat{\eta }}_i^r$ by a sufficiently small positive number κ: ${\widehat{\eta }}_i^r\leftarrow \kappa$ .
6.	Allocate: Increase the computing budget by Δ, and determine the new allocation using $$\begin{array}{c}\hfill {\widehat{N}}_i^{(\ast ,r)}=\frac{1/{\widehat{\eta }}_i^r}{{\sum }_{k=1}^{K}1/{\widehat{\eta }}_k^r}(n_{0}K+(r+1)\mathrm{\Delta }),\forall i\in \mathrm{\Theta }.\end{array}$$ 47 Allocate additional ${\delta }_i^{r+1}=\max \{{\widehat{N}}_i^{(\ast ,r)}-N_i^r,0\}$ simulation replications to solution ${\mathbf{x}}_i$ . Update $N_i^{r+1}\leftarrow N_i^r+{\delta }_i^{r+1}$ , and $r\leftarrow r+1$ . Go to Step 1.

In Step 5 of Algorithm 1, for

$i\in {\widehat{\mathrm{\Theta }}}_O$

,

${\widehat{c}}_i$

can be chosen between

${\widehat{p}}_{0\widehat{b}}$

and

${\widehat{p}}_{0i}$

, and we will show in the Supporting Information EC.2 that

${\widehat{c}}_i=\frac{{\widehat{p}}_{0\widehat{b}}+{\widehat{p}}_{0i}}{2}$

is a recommended choice. Further, parameter

${\eta }_i$

is strictly positive by definition, and thus

${\widehat{\eta }}_i^r$

is bounded by a sufficiently small positive number in Step 5. Finally, in Step 6, it may be necessary to round the values computed by Equation (47) to integers, for example, by preserving the total sum while minimizing the maximum relative deviation of such rounding.

While the asymptotic performance of Algorithm 1 is guaranteed by Theorem 2, a more important practical concern is its finite‐time performance. Specifically, two metrics are critical: (1) How close is the posterior mean

${\widehat{p}}_{hi}$

to the true success probability

$p_{hi}$

? and (2) How close is the empirical sample allocation ratio in Equation (47) to the asymptotically optimal one in Equation (38)? For the latter, we define

${\rho }_i^{\ast }=\frac{1/{\eta }_i}{{\sum }_{k=1}^{K}1/{\eta }_k}$

as the asymptotic optimal ratio of samples allocated to solution

${\mathbf{x}}_i$

, and

${\widehat{\rho }}_i^r=\frac{1/{\widehat{\eta }}_i^r}{{\sum }_{k=1}^{K}1/{\widehat{\eta }}_k^r}$

as the corresponding empirical ratio at iteration r. To this end, Theorem 3 characterizes the finite‐time performance of Algorithm 1, with the proof supported by Lemma 4. Lemma 4

Suppose the prior joint distributions

${\pi }_h^0({\mathit{\psi }}_h)$

(

$h=0,1,\cdots ,H$

) satisfy Assumption 1. Further, suppose at iteration r of Algorithm 1,

$N_i^r$

samples have been allocated to

${\mathbf{x}}_i$

out of a total of t samples, and define the estimated success probability

${\widehat{p}}_{hi}\doteq {\mathbb{E}}_t[{\mathrm{\Psi }}_{hi}]$

from Equation (11). Then, for any

$\gamma >0$

,

$\mathbb{E}[e^{\gamma W}]<\infty$

, where

$$\begin{array}{c}\hfill W=\underset{r\in \mathbb{N}\cup \{0\}}{\max }\underset{i\in \mathrm{\Theta }}{\max }\underset{h=0,1,\cdots ,H}{\max }\sqrt{\frac{N_i^r+1}{\log (e+N_i^r)}}|{\widehat{p}}_{hi}-p_{hi}|.\end{array}$$

48

Proof

The proof is given in the Supporting Information EC.1.3 to this paper.

$\square$

Theorem 3

Suppose Assumption 1 holds. For any sufficiently small

$\epsilon >0$

, let

49

be the earliest iteration such that for all iterations

$r\ge R^{\epsilon }$

in Algorithm 1, the posterior mean

${\widehat{p}}_{hi}$

and the empirical allocation ratio

${\widehat{\rho }}_i^r$

are sufficiently close to their asymptotic counterparts,

$p_{hi}$

and

${\rho }_i^{\ast }$

, respectively. Further, as the computing budget

$T\to \infty$

,

$\mathbb{E}[R^{\epsilon }]<\infty$

and

$Var[R^{\epsilon }]<\infty$

.

Proof

The proof is given in the Supporting Information EC.1.4 to this paper.

$\square$

Remark 5

By the Chebyshev inequality, the finite expectation and variance indicate that by a very high probability,

$EO{C}_{Bayes}$

will converge at the almost optimal convergence rate after a finite number of iterations. For example, let

$r^{\ast }=\mathbb{E}[R^{\epsilon }]+\sqrt{10Var[R^{\epsilon }]}$

, the Chebyshev inequality tells us that

$|{\widehat{p}}_{hi}-p_{hi}|\le \epsilon$

and

$|{\rho }_i^{\ast }-{\widehat{\rho }}_i^r|\le \epsilon$

for all

$i\in \mathrm{\Theta },h=0,1,\cdots ,H$

, and

$r>r^{\ast }$

with a probability greater than 90%.

Note that by Equation (19), the empirical value of

$AEOC\mathrm{-}B$

provided by Algorithm 1 converges to its theoretically optimal counterpart as the sample allocation converges. That is,

$|{\rho }_i^{\ast }-{\widehat{\rho }}_i^r|\le \epsilon$

means

50

where ε decreases to zeros as

$\epsilon \to 0$

. In other words, Algorithm 1 achieves the almost optimal value of

$AEOC\mathrm{-}B$

when

$r\ge r^{\ast }$

.

NUMERICAL EXPERIMENTS

In this section, we provide numerical results that show the effectiveness of the proposed algorithm. To this end, we test randomly generated test cases and a case study inspired by the real‐world staff allocation problem introduced in Section 3.1.

Random test cases

We first test the performance of Algorithm 1 (short for PB in the following context) using randomly generated test cases. Three benchmark algorithms are selected in comparison, namely, the OCBA for constrained optimization in Lee et al. (2012) (short for CO), the stochastically constrained R&S via SCORE in Pasupathy et al. (2015) (short for SC), and the heuristic procedure based on statistical hypothesis tests in Choi et al. (2021) (short for HE). Although the CO, SC, and HE algorithms do not consider the performance correlations in solutions as PB does, they are the most relevant algorithms developed to solve the constrained R&S problems under study in this paper. As such, we test the PB algorithm with an uncorrelated uniform prior (short for PB _u ) as a direct comparison to the three benchmarks. Further, we test the PB algorithm with a correlated prior (short for PB _c ) in order to quantify the benefit of using a correlated prior in the PB implementation. To this end, we generate random test cases with different parameter settings, and compute the optimal solutions using the aforesaid five algorithms by varying computing budgets. Note that the true optimum for each random test case can be analytically computed (to become apparent once we introduce the generation procedure for random cases below) and will be used as ground truths for performance comparison. Note further that our implementations of all algorithms strictly use sample approximations without leveraging the knowledge of the true optimum.

We compare the algorithm performance under different values of

$|\mathcal{X}|$

, H, and T. Recalling that the CO and HE algorithms were developed based on the assumption that samples follow a normal underlying distribution, we will further compare the performance with different distributions, including normal, uniform, and exponential. Given each parameter setting (

$|\mathcal{X}|$

, H, T, and an underlying distribution), 1,000 random test cases were generated based on the following procedure:

• 1.

  Generate the value of

  $d_h$
  for
  $h=0,1,\text{…},H$
  , each following a discrete uniform distribution
  $\mathcal{DU}(1,200)$
  .

• 2.

  Generate the value of

  ${\alpha }_h$
  for
  $h=1,\text{…},H$
  , each following a uniform distribution
  $\mathcal{U}(0.01,0.1)$
  .

• 3.

  For each solution

  ${\mathbf{x}}_i\in \mathcal{X}$
  for
  $i=1,\text{…},K$
  , generate the value of
  $F_h({\mathbf{x}}_i,\omega )$
  for
  $h=0,1,\text{…},H$
  , each following the chosen underlying distribution.

For each random test case, the value of

$\mathbb{E}[F_h({\mathbf{x}}_i,\omega )]$

depends on the randomly generated distribution parameters, while

${\mathbb{E}}_s[\mathbb{E}[F_h({\mathbf{x}}_i,\omega )]]$

increases with i where

${\mathbb{E}}_s[·]$

denotes the expectation over the randomness caused by the distribution parameters. As seen, the test cases were generated to mimic the correlation in solution performances in a way that the performance of a solution may be more correlated with its neighboring ones.

In the PB _c implementation, we set the prior distribution of

$\log ({\mathrm{\Psi }}_{hi}/(1-{\mathrm{\Psi }}_{hi}))$

,

$i\in \mathrm{\Theta }$

as a multivariate normal (Howard, 1998). The covariance function is squared exponential (Rasmussen & Williams, 2006); that is,

$\mathrm{\Sigma }({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }})={\tau }^2{e}^{-\varphi \parallel {\mathbf{x}}_i-{\mathbf{x}}_{i^{\prime }}{\parallel }_2^2}$

with

$\tau >0$

and

$\varphi >0$

. Let

$Cov$

be the covariance matrix whose

$(i,i^{\prime })$

entry is

${(Cov)}_{i{i}^{\prime }}=\mathrm{\Sigma }({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }})$

. Then, the prior belief has a density

${\pi }_h^0({\mathit{\psi }}_h)$

proportional to

51

where y is the expectation of the random vector

$\log ({\mathbf{\Psi }}_h/(1-{\mathbf{\Psi }}_h))=(\log ({\mathrm{\Psi }}_{h1}/(1-{\mathrm{\Psi }}_{h1})),\cdots$

,

$\log ({\mathrm{\Psi }}_{hK}/(1-{\mathrm{\Psi }}_{hK})))$

. In the experiments, we set

$\tau =1$

,

$\varphi =1$

, and

$\mathit{y}={(y_1,y_2,\text{…},y_K)}^{\top }$

to be

$y_i=\log (\frac{{\mathbb{E}}_s[p_{hi}]}{1-{\mathbb{E}}_s[p_{hi}]})$

, where

$p_{hi}=P\{F_h({\mathbf{x}}_i,\omega )>d_h\}$

. Here, the choice of y mimics practical situations when decision makers may have some prior knowledge on

$p_{hi}$

, but the prior is inaccurate if

$p_{hi}$

is predicted by

${\mathbb{E}}_s[p_{hi}]$

. We mention that the posterior mean does not have a closed form and is approximated using the Markov chain Monte Carlo approach.

For each random case, its true (theoretical) optimal solution for the attendant problem

$\mathcal{P}$

can be analytically computed based on the distribution information. Then, the same problem is solved using the five algorithms given a computing budget T. The solution obtained is compared to the true optimal solution. Three metrics are computed based on the 1,000 cases:

• Probability of correct selection (PCS): the percentage of cases where a true optimal solution is identified using the chosen algorithm;

• Probability of selecting an infeasible solution (PIF): the percentage of cases where an infeasible solution is identified as an optimal solution using the chosen algorithm;

• Expect opportunity cost (EOC): the expected performance gap between the selected solution and the true optimum. Specifically, the frequentist opportunity cost of selecting

  ${\mathbf{x}}_i\ne {\mathbf{x}}_b$
  is defined as
  $O{C}_i^{\ast }={\sum }_{h=0}^{H}O{C}_{hi}^{\ast }$
  , where
  $O{C}_{hi}^{\ast }=\max \{p_{hi}-{\alpha }_h,0\}$
  ,
  $h=1,2,\cdots ,H$
  , and
  $O{C}_{0i}^{\ast }=\max \{p_{0i}-p_{0b},0\}$
  . Let
  $O{C}_{no}^{\ast }=0.1$
  when
  $\widehat{b}$
  does not exist. Then,
  $EOC=O{C}_{no}^{\ast }P(\widehat{b}\text{does}\text{not}\text{exist})+{\sum }_{i\ne b}O{C}_i^{\ast }P(\widehat{b}=i)$
  , where the simulation results are used to estimate
  $P(\widehat{b}\text{does}\text{not}\text{exist})$
  and
  $P(\widehat{b}=i)$
  .

Here, PCS measures an algorithm's expected performance (i.e., how good it is in selecting the optimal solution); PIF measures its worst‐case performance (i.e., how bad it is in selecting an infeasible solution); and EOC measures the expected performance gap between the selected solution and the true optimum (i.e., the overall performance on the solution quality).

Table 2 exhibits the performance comparison results for the normal underlying distribution. Specifically, the table shows the PCS, PIF, and EOC metrics for three groups of settings,

$(K=200,H=1)$

,

$(K=500,H=3)$

, and

$(K=500,H=1)$

, where the difficulties of problems in each group increase in the aforesaid order. It is straightforward to see that group

$(K=500,H=1)$

is much more challenging than group

$(K=200,H=1)$

, while it may be counterintuitive to see that group

$(K=500,H=3)$

is also easier than

$(K=500,H=1)$

. The latter is due to the fact that each problem in group

$(K=500,H=1)$

is expected to have more feasible solutions than that in group

$(K=500,H=3)$

. Since identifying an infeasible solution is typically less computationally expensive than verifying a feasible solution, group

$(K=500,H=1)$

is much harder than group

$(K=500,H=3)$

. In summary, we characterize three groups of problems as follows:

$(K=200,H=1)$

has a small size of candidate solutions;

$(K=500,H=3)$

has a small size of feasible solutions and a large size of infeasible ones; and

$(K=500,H=1)$

has large sizes of feasible and infeasible solutions. Note that, for all algorithms, we set

$n_0=20$

and

$\mathrm{\Delta }=4,000$

for

$(K=200,H=1)$

, and

$\mathrm{\Delta }=10,000$

for both

$(K=500,H=3)$

and

$(K=500,H=1)$

. For each group, we further vary the values of computing budget T to show the convergence trend.

TABLE 2

Performance comparison with normal underlying distributions (

${\alpha }_h\sim \mathcal{U}(0.01,0.1)$

)

Setting		PCS (%)					PIF (%)					EOC (× 1,000)
$(K,H)$	T	CO	SC	HE	PB u	PB c	CO	SC	HE	PB u	PB c	CO	SC	HE	PB u	PB c
(200,1)	8,000	14.9	15.5	15.0	28.6	37.2	38.3	37.3	36.7	12.6	17.0	27.1	26.1	26.3	16.9	15.1
	16,000	35.2	29.1	29.7	54.2	62.7	22.6	24.0	28.2	4.2	5.8	13.7	13.8	14.9	7.8	6.2
	24,000	55.1	52.0	48.5	64.7	73.3	13.7	11.8	15.1	3.5	3.7	8.5	6.2	7.5	4.8	3.4
	32,000	64.4	68.0	66.7	74.6	77.7	11.2	7.1	6.3	2.2	2.8	6.7	3.3	4.2	3.6	2.9
	40,000	69.5	77.0	75.2	78.1	82.3	10.0	4.3	4.7	2.0	2.3	6.3	2.0	3.4	2.7	2.2
(500,3)	20,000	27.9	19.9	20.5	38.2	42.7	46.2	56.2	57.7	14.2	20.6	30.6	36.4	34.3	28.8	26.1
	40,000	59.1	28.1	28.1	69.6	74.3	17.8	44.5	45.2	5.6	7.3	15.1	26.7	25.8	10.8	8.5
	60,000	73.0	32.2	35.0	81.3	85.0	10.6	38.8	36.5	3.3	4.1	11.8	22.2	20.9	6.4	4.0
	80,000	77.6	37.4	38.8	85.1	88.7	8.2	34.4	33.6	2.5	2.6	10.2	20.7	19.2	5.3	2.9
	100,000	80.3	42.5	42.5	86.8	90.3	7.8	29.0	29.3	2.6	2.2	9.9	18.0	17.8	4.2	2.3
(500,1)	20,000	14.4	15.4	17.2	34.9	44.1	39.3	38.1	34.7	10.5	13.5	27.2	27.2	25.1	14.5	12.3
	40,000	40.2	23.2	24.8	63.3	74.1	20.1	30.2	27.1	5.6	3.6	11.9	16.5	15.3	5.5	3.3
	60,000	61.4	29.7	30.4	77.5	83.3	11.9	24.2	24.5	2.6	1.5	7.2	11.9	11.8	2.6	1.8
	80,000	71.9	34.5	33.3	82.9	86.7	9.9	20.2	21.6	1.9	1.4	5.9	10.0	9.9	2.1	1.4
	100,000	74.9	37.6	38.0	86.9	90.6	9.3	19.4	19.9	1.2	1.0	5.7	8.8	8.0	1.2	0.8

From Table 2, it is seen that, for the PCS metric, PB _u clearly outperforms three benchmarks for all three groups of problems; particularly, for the hardest (third) group, PB _u reaches a PCS of 77.5% when

$T=60,000$

, a value that the benchmark algorithms fail to reach even when

$T=100,000$

. For the PIF metric, PB _u shows good controls of not choosing an infeasible solution especially as T becomes sufficiently large, while the performance of the benchmarks are not as satisfactory. For the EOC metric, PB _u again outperforms benchmarks by a substantial margin for the harder problems (second and third groups), indicating a strong overall performance by taking into account the quality of the selected solution. Further, PB _c outperforms PB _u for all test cases in terms of PCS and EOC, showing the benefit of considering a correlated prior when performance correlations indeed exist. In terms of PIF, PB _c performs better in harder problems, while PB _u is superior in easy problems. Note that the advantage of PB _c over PB _u tends to diminish as T increases.

To check the robustness of PB across different underlying distributions, Table 3 and Table 4 exhibit performance comparisons with uniform and exponential underlying distributions, respectively. It is seen that PB is superior to all three benchmarks robustly across various underlying distributions. We also point out that the advantage of PB is more significant for the hardest group of cases, indicating that PB is a better choice when the size of problem becomes large.

TABLE 3

Performance comparison with uniform underlying distributions (

${\alpha }_h\sim \mathcal{U}(0.01,0.1)$

)

Setting		PCS (%)					PIF (%)					EOC (× 1,000)
$(K,H)$	T	CO	SC	HE	PB u	PB c	CO	SC	HE	PB u	PB c	CO	SC	HE	PB u	PB c
(200,1)	8,000	21.6	23.1	23.9	39.9	45.7	39.8	39.0	38.8	11.0	14.9	29.8	29.5	28.7	18.3	16.4
	16,000	44.9	42.3	36.7	61.7	73.3	23.5	24.3	30.2	5.7	6.0	14.6	14.0	15.8	8.0	4.6
	24,000	60.4	61.8	55.3	73.5	82.2	13.5	13.0	17.2	2.8	2.3	10.3	5.6	7.8	5.1	2.5
	32,000	70.0	74.2	68.2	79.9	84.9	10.8	8.3	11.2	2.6	2.2	8.8	3.2	4.2	3.7	1.8
	40,000	73.2	81.0	78.0	83.5	87.5	10.3	5.5	6.4	1.6	2.2	8.4	1.8	2.7	2.8	1.5
(500,3)	20,000	31.7	24.5	25.9	45.8	50.6	45.4	53.4	52.5	12.5	17.0	36.6	38.9	36.7	31.9	26.0
	40,000	59.0	57.2	33.7	76.5	78.4	18.0	21.9	42.6	4.8	5.0	20.8	17.8	26.8	9.8	8.8
	60,000	69.9	70.6	43.5	85.3	87.6	11.1	10.4	33.1	2.7	3.3	16.7	14.1	22.8	5.6	4.8
	80,000	75.3	77.1	55.3	88.3	92.0	8.6	6.2	20.8	1.7	2.6	14.6	12.7	16.9	4.6	2.8
	100,000	77.7	80.8	64.6	90.6	93.3	7.9	5.1	14.8	1.3	1.7	13.2	11.8	14.4	3.5	2.6
(500,1)	20,000	21.5	23.9	24.2	47.8	56.6	39.3	39.5	36.8	7.2	10.7	31.1	28.8	29.8	16.1	11.8
	40,000	47.7	43.9	38.2	77.7	82.1	21.4	23.0	27.4	3.1	2.6	13.3	13.4	15.4	4.1	2.6
	60,000	68.0	65.1	60.4	86.1	90.4	12.2	9.9	15.5	1.7	0.9	7.4	5.4	7.1	2.2	0.9
	80,000	77.2	79.9	74.6	91.0	92.7	9.8	4.8	8.4	1.1	0.9	6.3	2.6	4.0	1.5	0.7
	100,000	80.0	87.5	84.8	92.6	94.7	9.8	2.7	4.0	0.4	0.3	6.3	1.8	2.5	1.3	0.6

TABLE 4

Performance comparison with exponential underlying distributions (

${\alpha }_h\sim \mathcal{U}(0.01,0.1)$

)

Setting		PCS (%)					PIF (%)					EOC (× 1,000)
$(K,H)$	T	CO	SC	HE	PB u	PB c	CO	SC	HE	PB u	PB c	CO	SC	HE	PB u	PB c
(200,1)	8,000	22.2	22.6	20.9	42.3	44.8	39.3	41.2	40.7	10.6	13.0	31.7	30.6	30.3	17.4	17.6
	16,000	44.6	43.2	37.7	66.5	70.8	21.9	23.3	27.9	5.2	4.0	15.5	14.0	15.1	7.1	6.1
	24,000	60.3	62.3	56.1	76.4	82.3	13.5	11.7	16.1	2.1	2.2	9.9	6.7	8.5	4.4	3.4
	32,000	69.1	74.5	66.8	82.0	85.7	10.7	6.6	9.8	1.4	2.7	8.5	3.6	6.3	3.0	2.7
	40,000	73.1	82.7	76.2	85.6	88.4	10.2	4.7	6.3	1.2	1.4	7.9	2.3	4.3	2.5	2.3
(500,3)	20,000	31.8	25.0	22.2	43.9	50.2	48.9	54.6	54.2	15.0	17.8	34.0	36.7	38.9	31.2	24.0
	40,000	62.0	55.1	34.1	76.3	79.7	17.5	22.2	40.8	4.3	5.9	16.9	17.7	26.5	8.6	7.2
	60,000	74.4	72.1	44.0	84.2	88.1	11.4	10.5	31.9	2.9	3.3	12.8	11.8	20.2	5.6	4.1
	80,000	80.4	79.0	57.2	88.0	91.7	6.8	6.3	20.4	2.3	1.9	11.3	10.5	16.0	4.2	2.8
	100,000	81.5	82.8	66.9	90.5	93.4	7.3	4.5	14.5	1.9	1.9	10.9	9.5	13.3	2.9	2.3
(500,1)	20,000	21.7	21.6	20.4	45.0	55.0	39.6	42.6	40.7	8.6	11.7	29.8	29.8	30.0	18.6	12.7
	40,000	48.3	45.1	37.3	75.1	80.6	21.9	22.4	26.9	2.8	3.6	14.6	13.1	14.5	4.9	3.1
	60,000	66.9	66.3	55.1	85.0	90.3	14.7	10.0	18.3	2.2	1.3	8.9	5.6	8.4	2.2	1.7
	80,000	73.2	80.5	72.7	88.8	91.6	12.2	4.7	8.2	1.5	1.3	8.0	3.0	5.1	1.7	1.4
	100,000	75.7	87.1	81.7	91.4	93.4	11.4	3.3	5.5	1.1	0.7	7.6	1.9	3.4	1.2	1.1

In the Supporting Information EC.3, we modified the procedure of generating random test cases by setting the values of

${\alpha }_h(h=1,\cdots ,H)$

to follow a uniform distribution

$\mathcal{U}(0.01,0.3)$

instead of

$\mathcal{U}(0.01,0.1)$

, capturing a wider range of service levels seen in applications. With this modification, the random test cases are expected to have more feasible solutions in general. The performance comparisons are shown in Tables EC.1, EC.2, and EC.3 in the Supporting Information, and the observations are quite similar to those in Tables 2, 3, and 4, showing the robustness of PB across different settings of problem

$\mathcal{P}$

by varying

${\alpha }_h$

. Further, we varied the incremental budget Δ in the PB algorithm implementation. The impact of Δ on performance and suggestions on the choice of Δ are also discussed in the Supporting Information EC.3.

In summary, these numerical results show that, compared to CO, SC, and HE, PB is more likely to correctly select the optimal solution, and less likely to select an infeasible one; even if an optimal solution is not selected, PB is much more likely to select a near‐optimal one. All these are desired characteristics for an R&S algorithm. Further, the performance improvement is quite significant even when the computing budget is relatively small. This will benefit many applications as each simulation replication can be expensive and time consuming to run and the decision‐making time frame may be short (e.g., within hours or half a day). Finally, when prior knowledge is available for a problem, it would be beneficial to incorporate a correlated prior.

An ED case study

In this section, we will demonstrate how the proposed model and algorithm can be used to improve services in the ED example introduced in Section 3.1. A discrete‐event simulation model is used to evaluate the performance of the highly complex and dynamic processes in ED. The data used in this case study were collected from July 1, 2009, to June 30, 2010, and the service process and simulation model used have been verified in Guo et al. (2017) and Chen et al. (2020).

From the historical data obtained, the ED received various patients with the volume varying from 350 to 430 per day, where the patients fall into five categories with the following percentages on average (from categories I to V): 0.87%, 1.45%, 42.09%, 48.76%, and 6.83%. In the simulation model, the patient arrivals follow a nonhomogeneous Poisson process, whose hourly arrival rates were calculated empirically based on the historical data and are shown in the Supporting Information EC.4.

The discrete event simulation model follows the actual service processes (Chen et al., 2020), where the detailed process map and the distributions used to generate the service times at each process are displayed in the Supporting Information EC.4. Specifically, an incoming patient, depending on the situation, may go through different service areas, including registration and triage, resuscitation (for critical patients), consultation, laboratory test and imaging investigation, and observation. Main medical resources required include different types of nurses, doctors with different seniority levels and expertise, different medical devices (e.g., electrocardiography [ECG] and ultrasonography [USG]), as well as different spaces (offices, rooms, and labs). Medical staff worked three shifts around the clock: morning shift (from 8 a.m. to 4 p.m.), evening shift (from 4 p.m. to 12 p.m.), and night shift (from 0 a.m. to 8 a.m.).

The simulation model was implemented in Arena, and was validated using the 1‐year data collected. For the model validation, we set the length of the simulation to be 1 month with a 7‐day warm‐up period and ran the simulation model for 500 replications. We used three indicators, waiting time at the triage (WTT), waiting time for physicians (WTP), and length of stay (LOS), to assess the validity of outputs, as shown in Table 5. Note that the data of WTT and WTP for category I patients are not available from the cooperative hospital. It is seen that the simulation can closely mimic the real data.

TABLE 5

Simulation model validation results

		CAT I	CAT II	CAT III	CAT IV	CAT V
WTT (h)	Real data	–	0.04	0.09	0.11	0.11
	Simulated result	–	0.03	0.11	0.11	0.11
WTP (h)	Real data	–	0.11	0.31	1.65	2.07
	Simulated result	–	0.07	0.27	1.64	2.17
LOS (h)	Real data	1.30	1.23	3.53	4.15	3.13
	Simulated result	1.56	1.51	3.49	4.37	3.20

Next, we experiment with the simulation model to identify the optimal solution

${\mathbf{x}}^{\ast }$

that optimizes the resource allocation problem defined in Equations (1)–(5). Recall that each candidate solution

$\mathbf{x}$

consists of the medical staff assignment for three 8‐h shifts of a day. The medical staff is composed of the admission nurse (AN), triage nurse (TN), junior nurse (JN), senior nurse (SN), junior doctor (JD), and senior doctor (SD). To protect the data privacy, we set the labor cost for each type of staff using the cost unit (CU) as follows: 1 CU per AN, 1 CU per TN, 1.5 CUs per JN, 2 CUs per SN, 3 CUs per JD, and 4 CUs per SD. For the actual solution implemented in the ED, the total staffing cost is 85.5 CUs. By applying the following filters to all possible solutions, we identified 73 promising candidate solutions for the R&S problem (

$|\mathcal{X}|=73$

).

• Solutions were filtered by the deterministic constraints specified in Equation (5), including the personnel and budget constraints. Specifically, the personnel constraint imposes practical bounds on the number of medical staff of each type and shift. Any solution resulting in a staffing cost of more than 85.5 CUs was also filtered, since the ED preferred maintaining the current budget.

• Solutions dominated by others were filtered. For example, consider solution A with 2 SDs for the night shift and solution B with 1 SD for the night shift, while all other staff assignments are equal. We expect solution A to dominate solution B without evaluating them since the objective is to find a solution that maximizes the service level of type IV patients, thereby excluding solution B from the solution space.

• Solutions were further filtered by consulting domain experts, who would remove nonpromising solutions based on their prior experience on staff schedules or preliminary simulation runs.

The full set of the 73 candidate solutions is included in the Supporting Information EC.4. We mention that the above filtering process substantially reduces the solution space, making it suitable to construct an R&S problem studied in this paper, while the original problem without filtering is more suitable to be solved using a search‐based algorithm, such as the one developed in Chen et al. (2020).

We compare the performance of PB to two benchmarks: One is the CO in Lee et al. (2012) and the other is the equal allocation (EA). The EA approach is a simple strategy of allocating the computing budget equally among all candidate solutions, and is often used in the literature as a baseline. The computations were conducted in Arena 14.0 via VBA codes. The prior belief was set to have a density

${\pi }_h^0({\mathit{\psi }}_h)$

proportional to Equation (51), where

$\tau =100$

,

$\varphi =1$

, and y is a vector of zeros. We set the total computing budget

$T=2,117$

, and therefore each solution was simulated for 29 replications using EA. For PB and CO, we ran

$n_0=10$

initial replications for each solution, and set the incremental budget

$\mathrm{\Delta }=73$

. Note that among the set of solution candidates, only one of them has been implemented in practice and thus no ground‐truth values are known for these solutions. To obtain approximated ground‐truth values, we simulated each solution for 500 replications (sufficiently large) and selected the best solution as the estimated optimal solution. Then, we ran each algorithm to select the best solution and compared it to this estimated optimal solution. For each algorithm, we ran the same experiment for 60 times, and Figure 1 plots the average values of PCS and EOC evolving as simulation replications increase.

FIGURE 1

Performance comparison on the ED example

It is seen from Figure 1 that the PCS converges to 1 and EOC converges to 0 faster when using PB than using CO, while the PCS and EOC do not converge within the tested range of the computing budget when using EA. EA is apparently not a good strategy, and needs a much larger computing budget to conclude a reliable solution. To further illustrate that the optimal solution obtained by PB is indeed better than the actual staff allocation schedule, we simulated both solutions for 500 replications. Table 6 shows the performance comparison between the actual solution implemented by the hospital and the optimal solution obtained by PB, where

$P_{W{T}_2}$

,

$P_{W{T}_3}$

,

$P_{LOS}$

, and

$P_{W{T}_4}$

are the probabilistic measures defined on the left‐hand side of Equations (2), (3), (4), and (1), respectively. It is seen that the optimized solution meets all probabilistic service levels specified in Table 1, while the actual solution fails to do so (with a small margin). Further, the optimal solution performs much better in terms of the objective (1) (i.e.,

$P_{W{T}_4}$

) compared to the actual solution, with a 31.80% of improvement on this measure. The total staff costs for both solutions are very close and within the specified budget. The comparison shows the superiority of the optimized solution, and confirms the necessity of optimizing the resource allocation for this ED.

TABLE 6

Comparison between the actual and optimal solutions

		Staff allocation						Service level
	Shift	AN	TN	JN	SN	JD	SD	$P_{W{T}_2}$	$P_{W{T}_3}$	$P_{LOS}$	$P_{W{T}_4}$	Cost
Actual solution	Morning	1	2	2	4	2	3	93.79%	89.10%	88.31%	67.59%	85.5
	Evening	1	2	2	4	2	3
	Night	1	1	1	2	2	2
Optimal solution	Morning	2	2	3	2	3	3	97.08%	91.41%	98.46%	99.39%	84.5
	Evening	1	2	2	2	2	3
	Night	1	1	2	2	2	2

CONCLUSION

The purpose of this research is to formulate a model to study resource allocation problems in service systems, where the system performance is evaluated via simulation. Particularly, the objective and constraints in the class of problems under study are in the form of probabilistic measures. Simulation models are known to be capable of capturing complex processes and system dynamics, which are hard to model in mathematical programs but are often present in service systems. However, a major challenge in optimizing systems using simulation models, known as the research field of simulation optimization or optimization via simulation, is the trade‐off between simulation accuracy and the affordable computing budget. To this end, we formulate an OCBA model to efficiently balance the effort of exploring the solution space and exploiting given solutions with more simulation replications. From the methodological point of view, the OCBA problem formulated in this paper possesses a couple of unique features compared to the literature: (1) Instead of maximizing the PCS, the problem under study aims to minimize the EOC, considering that the service provider is typically more concerned with the quality of the selected solution than the statistical significance of selecting the true optimal one; and (2) the proposed model considers the prior knowledge accumulated on solutions, as well as the performance correlations on solutions, and hence uses a Bayesian modeling approach. To the best of our knowledge, such a Bayesian OCBA model has not been studied in the literature. Next, based on the property of the probabilistic measures, we derive the asymptotic optimality conditions of the OCBA formulation. Then, an efficient iterative algorithm is developed to sequentially allocate the simulation budget to the most necessary solutions, guided by the optimality conditions. The efficiency of the algorithm is shown via the theoretical finite‐time convergence analysis, as well as the numerical experiments.

From a practical vintage point, the generalized resource allocation problem is motivated by a real‐life resource allocation problem stemming from a hospital ED, and is widely applicable to service systems, such as hospitals, call centers, and public service offices. The proposed algorithm can expedite the identification of an optimal resource allocation scenario using a simulation model, with the presence of probabilistic constraints and objective. Extensive numerical examples with known true optimum conducted in this paper have demonstrated the superiority of the proposed algorithm compared to benchmarks, evaluated by the PCS, PIF, and EOC metrics. Further, we apply the algorithms to the aforesaid staff allocation problem in the ED in Hong Kong, whose simulation model is driven by real data. This case study has further validated the practicality of the proposed model and algorithm.

Several future research directions are available to explore. Sometimes, the service system may want to guarantee the PCS (e.g., 90%) and minimize the required computing budget (Chen et al., 2014). Thus, a future research direction is to explore alternative simulation optimization methods, such as the IZ approaches. Further, this research considers problems with a manageable size of solution space. In the proposed algorithm, every solution is simulated for a minimum number of replications and evaluated to obtain an initial assessment. For a large‐size problem, it may not be computationally feasible to provide an initial evaluation for every solution. As such, a metamodel‐based or random search–based method needs to be developed to balance the exploration of unknown solutions and the exploitation of known solutions with more simulation replications. Finally, when the resource constraints in Equations (1)–(5) or Equations (6)–(8) can be expressed as mathematical constraints, it would be interesting to explore a combined simulation optimization and mathematical programming approach.

Footnotes

ACKNOWLEDGMENTS

We gratefully thank the department editor, the senior editor, and the two reviewers for their valuable comments that have greatly improved the paper. Siyang Gao thanks the support by the City University of Hong Kong (Grants 7005269 and 7005568). Jianzhong Du thanks the support by the National Natural Science Foundation of China (Grant 72091211).

ORCID iD

Weiwei Chen

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.49 MB