A simulation-based optimization approach for the calibration of a discrete event simulation model of an emergency department
A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Roma
AA simulation-based optimization approach for thecalibration of a discrete event simulation model of anemergency department
Alberto De Santis · Tommaso Giovannelli · Stefano Lucidi · Mauro Messedaglia · Massimo RomaAbstract
Accurate modeling of the patient flow within an Emergency Depart-ment (ED) is required by all studies dealing with the increasing and well-knownproblem of overcrowding. Since Discrete Event Simulation (DES) models are oftenadopted with the aim of assessing solutions for reducing the impact of this world-wide phenomenon, an accurate estimation of the service time of the ED processesis necessary to guarantee the reliability of the results. However, simulation modelsconcerning EDs are frequently affected by data quality problems, thus requiring aproper estimation of the missing parameters.In this paper, a simulation-based optimization approach is used to estimate theincomplete data in the patient flow within an ED by adopting a model calibrationprocedure. The objective function of the resulting minimization problem repre-sents the deviation between simulation output and real data, while the constraintsensure that the response of the simulation is sufficiently accurate according to theprecision required. Data from a real case study related to a big ED in Italy is usedto test the effectiveness of the proposed approach. The experimental results showthat the model calibration allows recovering the missing parameters, thus leadingto an accurate DES model.
Keywords
Simulation-based Optimization · Emergency Department · DiscreteEvent Simulation · Model Calibration · Derivative-Free Optimization
Over the last few years, Emergency Departments (EDs) have been raising anincreasing attention in the Operations Research and Management Science com-
Alberto De Santis, Stefano Lucidi, Tommaso Giovannelli, Massimo RomaDipartimento di Ingegneria Informatica Automatica e Gestionale “A. Ruberti”, SAPIENZAUniversit`a di Roma, via Ariosto, 25 – 00185 Roma, Italy.E-mail: [email protected], [email protected], [email protected],[email protected] MessedagliaACTOR Start up of SAPIENZA Universit`a di Roma, via Nizza 45, 00198 Roma, Italy. E-mail:[email protected] a r X i v : . [ m a t h . O C ] F e b A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Roma munities due to the international phenomenon of the overcrowding [19,30,36,11,28], which leads to longer waiting times, higher mortality rates, and lower pa-tient satisfaction [31]. One of the most popular tools adopted to study this crucialproblem is Discrete Event Simulation (DES) [39,20,37,7,23,13], which is used torepresent the complex and stochastic patient flows in the ED in place of analyticalmodels (for a recent literature review of simulation modeling in EDs, see [32]).Frequently, DES models are combined with an optimization algorithm to makeoptimal decisions in relation to some Key Performance Indicators (KPIs). The re-sulting Simulation-Based Optimization (SBO) approach [14,15,3], whose effectiveapplication relies on the accuracy of the simulation model, has been frequentlyapplied to hospital EDs (for a recent review, see, e.g., [38]).Many papers in the literature deal with the quality of the input data used inED simulation models. This issue, which strongly affects the reliability of the re-sults, has been carefully analyzed. Indeed, the cost, time, and challenges requiredby collecting ED empirical data represents a serious limitation for every simulationmodel [30]. In order to accurately replicate and predict the patient flows within theED, [12] develops a process mining approach which handles the noise factors in thedataset after introducing assumptions on how to interpret the data. In particular,these factors include the following cases: starting and ending time of each activityperformed in the ED may be unknown; information about urgent patients may beregistered after the activity is completed and, in general, the timestamps of eachactivity may not be promptly recorded at the right time; for each patient, thesame activity may be recorded multiple times for technical reasons, giving rise tomisleading information. A framework to categorize all the ED data quality issuesis proposed in [34], which also provides assessment techniques for each data qualityproblem category. Moreover, this paper highlights that most of the works in theliterature focus on the problem of missing data, which is the problem addressed inthis paper. It consists in the lack of information on some or all the key timestampsthat define the activities performed in the ED, i.e., starting and ending time. Thiswell-known issue prevents gaining knowledge about the duration of each activity,which is required to estimate the corresponding probability distributions to usein the simulation model. Several SBO approaches are proposed in the literatureto tackle this crucial problem. The idea behind each approach is to leverage theknown information in order to estimate the parameters of the probability distri-butions underlying the missing data. This is accomplished by comparing the KeyPerformance Indicators (KPIs) computed through the simulation model with thecorresponding values derived by the data collected in the real system. The resultingprocedure is known as model calibration .The mathematical formulation of the optimization problem and the specificalgorithm used for determining the optimal solution are the two features that dis-tinguish the papers dealing with missing data. For instance, both [24] and [18]adopt similar approaches that leverage the time differences between the knowntimestamps. The former proposes an unconstrained optimization problem wherethe objective function is a consistency measure that compares the average, the standard deviation, and the proportions of the time differences for each triage tag.The latter uses a constrained optimization problem where the objective functionis based on a modified chi–square goodness of fit and the constraints are intro-duced to guarantee each time difference the same level of accuracy. As regards theapproaches used to solve the problem, both papers use metaheuristic procedures. imulation-based optimization for ED DES model calibration 3
In particular, [24] considers both a descent method and a simulated annealingalgorithm, while [18] adopts an approach that combines a genetic algorithm withsimulated annealing and optimal computing budget allocation. Moreover, in [18]the authors point out that the approach in [24] could be improved in several ways:including in the objective function all the time differences defined by the availabletimestamps, since delays in one activity may impact on downstream activities inthe patient flow; modeling the possibility for each patient to have more than onemedical visit, which may affect the time differences; using a more formal objectivefunction and solving the resulting optimization problem through a more efficientalgorithm.Another paper in the literature on missing data that is worth being mentionedis [26], which assumes an agent–based simulation model, as opposed to the DESmodel considered in [24] and [18]. Although the framework underlying this workis different, a simulation-based optimization problem is used for the same goalof minimizing the deviation between real data and simulation output in order toestimate the missing parameters of the simulation model. The Length Of Stay(LOS), i.e., the difference between the discharge time and the arrival time to theED, is the time difference considered in the objective function, which adopts theJensen–Shannon divergence. A systematic method based on the pattern searchmethod APPSPACK [16] is used to find the optimal configuration of parameters.The approach proposed in this paper aims to handle the problem of missingtimestamps, which affects many simulation models dealing with EDs, as evidencedby the papers in the specific literature. Among the sources of noise affecting thequality of the input data necessary to build a reliable simulation model, the prob-lem of missing timestamps has a strong impact on the overall accuracy of the sim-ulation, since it prevents the knowledge of the values used to derive appropriateprobability distributions. Therefore, while other issues may be resolved by eithercarefully cleaning the dataset or introducing assumptions on how to interpret thedata, the unavailability of timestamps requires more sophisticated procedures.Compared to the other papers dealing with missing data in ED datasets, theapproach discussed in this paper aims both to propose a new formulation of theresulting optimization problem and to improve on the optimization strategies typ-ically used in the literature. Since building a simulation model is a process thatrequires a considerable amount of effort and thus is not expected to be completedin a short time, an exact algorithm providing optimal solutions may be preferableto metaheuristic procedures, whose final solutions are returned faster but withoutoptimality guarantees. However, although global convergent algorithms appear tobe a reasonable choice, metaheuristics are the methods mainly adopted in theliterature dedicated to missing data (see, e.g., [24,18]). To fill this gap, in thispaper a SBO approach is developed both to propose an alternative version of theoptimization problems used in [24,18] and to use a Derivative-Free Optimization(DFO) strategy based on global convergence, which allows the algorithm to findan optimal solution with optimality guarantees. Similarly to [24] and differentlyfrom [18], Weibull distributions are adopted to generate the values of the activity service times associated with missing data. Indeed, this probability distribution isconsidered suitable when data is unavailable (see, e.g., [25]). Moreover, the mostcritical patients, who are excluded from the two approaches developed by [24]and [18], are included to prevent the simulation model from returning inaccurateresults for this important class of patients. Finally, in the proposed constrained
A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Roma optimization problem, the objective function is based on the comparison betweenreal and simulated Empirical Cumulative Distribution Functions (ECDFs), whichdo not rely on intervals, unlike the functions adopted in [24,18,26].This paper is organized as follows. Section 2 reports the typical structure ofthe datasets related to the patient flow within an ED, highlighting the availableinformation and the KPIs of interest. Section 3 defines the SBO problem used tocalibrate a general DES model concerning an ED. Section 4 describes the real EDused as a case study and the DES model implemented. Section 5 reports the resultsobtained by applying the SBO approach to the case study. Finally, in Section 6we draw some concluding remarks and we discuss directions for future work.
The timestamps defining the activities in the patient flow are represented in Figure2.1 and described in Table 2.1. Although these timestamps are commonly recorded
Fig. 2.1
Timestamps collected throughout the patient flow.
Table 2.1
Description of the timestamps collected throughout the patient flow. t starting time of triage. t ending time of triage. t starting time of medical visit. t ending time of medical visit and starting time of examinations. t ending time of the last exam or reassessment. t patient receives the last medical report. t patient is discharged and leaves the ED. by the electronic systems used to collect data, the actual patient flow may includeextra times omitted from this scheme. For example, since the arrival time to theED is usually not registered, before triage there may be an additional waiting time not included in any record. Moreover, extra waiting times and observationperiods may be present both in the examination part and in the discharge phase.Since usually there is no available data for these times, in Figure 2.1 both waitingand service times for exams are included in a single box between t and t , whichare associated with the end of the medical visit and the end of the last exam or imulation-based optimization for ED DES model calibration 5 reassessment after a visit, respectively. After receiving the medical report, patientsmay be either subject to additional medical visits or discharged after a final waitingtime, which is due to a further observation or to the time required by the ED staffto prepare the discharge. It is important to point out that the disposition, whichrefers to the decision to admit a patient to a hospital ward, is associated with t ,like the discharge. This implies that the boarding time of the admitted patients(i.e., the waiting time before being transferred to the hospital ward) is assumed tostart at t , thus being excluded from the ED scope. This choice is in accordancewith the case study considered in Section 4. Contrarily, in some EDs the dischargeof patients admitted may be considered as the time of admission to the hospitalward (see, e.g., [6]).In most cases, the available timestamps are t , t , and t , which are markedin the scheme as known . Since this is a typical setting in practice, throughout thispaper they are supposed to be the only timestamps recorded by the ED along with t (which may not be always available). As a consequence, for each patient it ispossible to compute – the Door-to-Doctor Time ( DOT ), which is the time difference between thestart of the triage and the start of the medical visit, namely t − t ; – the Doctor-to-Discharge Time ( DIT ), which is the time difference between thestart of the medical visit and the discharge, namely t − t if t is available, t − t otherwise.It is important to remark that when patients require additional medical visits afterthe exams, DIT can be interpreted as the time difference between the timestampof the last medical report (or the discharge if t is not available) and the startingtime of the first medical visit. Note that t and t are equal if a patient does notneed additional medical visits.Since in many cases only a subset of the timestamps related to the patientflow are known [34], the service time cannot be computed for all the activities.Apart from the final waiting time, which can be recovered for each patient throughthe difference t − t , note that the durations of all the other activities are notcompletely defined. In fact, in case of triage and medical visit, the ending time ismissing, while for exams both timestamps are unknown. The goal of the approach proposed is to recover the information needed to buildan accurate simulation model by leveraging the known information through theminimization of the deviation between real data and simulation output. In order todefine the simulation-based optimization problem, let us introduce the sets below. – Let C be the set of the triage tags. – Let U ( c ) be the set of units where patients with tag c ∈ C can be visited andtreated. – Let T = { DOT, DIT } be the set of the time differences considered.In the absence of data, suitable probability distributions are Weibull and lognormal(see [25]). Let us arbitrarily consider a Weibull distribution, whose probability A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Roma density function is f ( w ) = (cid:40) αβ − α w α − e − ( w/β ) α if w > , α > β > c and on the unit u . This choiceleads to a number of different pairs of parameters equal to (cid:80) c ∈ C | U ( c ) | , where | U ( c ) | refers to the cardinality of the set U ( c ). Let us denote as x ∈ R n , y ∈ R n ,and z ∈ R n the corresponding probability distribution parameters for the servicetimes considered, where n v ≤ (cid:80) c ∈ C | U ( c ) | with v ∈ { , , } . In particular, forall c ∈ C and for all u ∈ U ( c ), the shape and scale parameters of the Weibulldistributions are – x c u and x c u for the triage probability distribution, – y c u and y c u for the medical visit probability distribution, – z c u and z c u for the exams probability distribution.Let F simc u i and F realc u i be the ECDFs of the values of the simulated and real timedifference i ∈ T for patients with tag c visited in unit u . Moreover, let k simc u i and k realc u i be the number of such patients from the simulation and from the real dataset.Hence, for all j ∈ { sim, real } , we have that F jc u i ( t ) = 1 k jc u i k jc u i (cid:88) h =1 { τ c u i h ≤ t } with t ≥ , where τ c u i h is the value of the time difference i recorded for the h –patient, with h ∈ { , . . . , k jc u i } , considered from j ∈ { sim, real } . It is important to point outthat the values τ c u i h of the time differences computed from the simulation dependon the service times of triage, medical visit, and exams drawn from the Weibulldistributions described above. Therefore, this dependence can be written explicitlyas F simc u i ( t ; x, y, z ), where x , y , and z are the vectors containing all the associatedshape and scale parameters.The mathematical problem formulation is reported as follows min x, y, z (cid:88) c ∈ C (cid:88) u ∈ U ( c ) (cid:88) i ∈T (cid:16) (cid:90) ∞ ( F simc u i ( t ; x, y, z ) − F realc u i ( t )) dt (cid:17) s.t. x ∈ P , (3.1) imulation-based optimization for ED DES model calibration 7 where the feasible set P = (cid:110) ( x, y, z ) ∈ R n × R n × R n | g c u i ( x, y, z ) ≤ ,h c u i ( x, y, z ) ≤ ,l x ≤ x ≤ u x l y ≤ y ≤ u y l z ≤ z ≤ u z for all c ∈ C, u ∈ U ( c ) , i ∈ T (cid:111) is defined by the following functions – g c u i ( x, y, z ) = (cid:12)(cid:12) µ simc u i ( x,y,z ) − µ realc u i µ realc u i (cid:12)(cid:12) − tol c u iµ , which compares the sample means µ simc u i ( x, y, z ) and µ realc u i of the time difference i computed through the simulated(by averaging over the independent replications) and real data, respectively, – h c u i ( x, y, z ) = (cid:12)(cid:12) σ simc u i ( x,y,z ) − σ realc u i σ realc u i (cid:12)(cid:12) − tol c u iσ , which compares the sample stan- dard deviations σ simc u i ( x, y, z ) and σ realc u i of the time difference i computed through the simulated (by averaging over the independent replications) and real data,respectively,and l x , l y , l z , u x , u y , and u z are vectors defining bound constraints. Positive thresh-olds, namely tol c u iµ and tol c u iσ , are used to state the degree of accuracy requiredfor the simulation model. The objective function is the sum of the integrals ofthe squared difference between F simc u i and F realc u i over the sets C , U ( c ), and T .The choice of using ECDFs instead of histograms, which are commonly adoptedin the literature [24,18], is due to the different reliability of the information pro-vided. In particular, compared to ECDFs, the description of the data provided byhistograms is strongly affected by the choice of the width of each bin, which isnot straightforward and may lead to distributions with different shapes based onhow data is grouped. The decision variables are the Weibull distribution param-eters contained in the vectors x, y , and z . For the sake of simplicity, the previousformulation does not include variables that are not parameters of probability dis-tributions, although they may be present. It is important to remark that in thegeneral framework described, the dependence of each pair of parameters on both c and u implies that the triage service time is affected by the triage tag and the EDunit. This assumption, which turns out to be reasonable when applied to the du-ration of medical visit and exams, may lead to an excessive number of variables forthe triage service time, if it does not hold. Indeed, while the color tag significantlyaffects the triage duration since urgent patients undergo a faster triage than less critical patients, the impact of the unit may be negligible. However, as concernsthe case study described in Section 4, interviews with the ED staff have shown thatslightly different procedures are adopted by the nurses in charge of triage basedon the unit where a patient is assigned. This means that final conclusions shouldbe drawn after analyzing the specific system considered. Since the probability dis-tribution parameters affect the dimension of the optimization problem, avoiding A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Roma unnecessary variables allows the algorithm to benefit from a lower computationalcost.
The case study concerns the ED of
Policlinico Umberto I , which is a very largehospital in Rome, Italy. It is the biggest ED in the region of Lazio in terms ofnumber of patient arrivals per year (about 140,000 on average). By using the datacollected from the patient flow through the ED for the whole year 2018, this casestudy is adopted to test the effectiveness of the approach proposed in this paper.4.1 Description of the EDThe ED of Policlinico Umberto I is divided into several areas, each one associatedwith a medical specialty. The backbone of the ED is the central area , which isdevoted to treating diseases and disorders related to internal medicine and generalsurgery, which affect the majority of patients. Separated from this main area,there are other parts of the ED that deal with the following medical specialties:ophthalmology, obstetrics, pediatrics, hematology, and dentistry. The focus of thispaper is on the central area, whose number of arrivals in 2018 amounted to morethan 50,000.In order to gain a complete understanding of the ED processes, the descriptionof the ED units and staff and the summary of the patient flow are reported in thesequel. Other than a triage area , where each incoming patient is assigned a colortag by a nurse in charge of this task, the ED is composed by • a Medical Unit (MU) , devoted to patients needing specialized medications andtreatments, with areas dedicated to the most critical patients; • a Surgical Unit (SU) , devoted to patients needing either to receive a surgi-cal operation or to recover from it, with areas dedicated to the most criticalpatients; • a Resuscitation Area (RA) , for the most acutely ill and injured patients, whoneed timely treatments; • a Minor Injuries Unit (MIU) , for the least urgent patients, whose treatmentcan be delayed or deferred; • an Orthopedic Unit (OU) , for patients suffering from orthopedic disorders.Moreover, all of these units have rooms where patients can either wait for examsor stay for observation. Red-tagged patients can be visited and treated in RA orin dedicated areas within MU and SU, which are open 24 hours a day and areprovided with equipment and staff specialized for dealing with life-threateningillnesses and injuries. In particular, 1 and 2 seats are available in the dedicatedareas of MU and SU, respectively, and further 2 seats are available in RA. As concerns the medical treatment of the other patients, MU and SU can host upto three and two patients during the day (8.00 a.m.–8.00 p.m.), respectively. Atnight, MU can host two patients, while in SU one seat is available. Moreover, MIUhas two seats, which are available from 8.00 a.m. to 8.00 p.m., Monday throughSaturday. When patients experience excessive waiting times, two additional seats imulation-based optimization for ED DES model calibration 9 may be used to visit up to four patients simultaneously. All this information issummarized in Tables 4.1–4.2.
Table 4.1
Number of seats available for medical visit and treatment in MU, SU, and MIU,the latter being open from Monday to Saturday.
MU SU MIUDay (8.00 a.m.–8.00 p.m.) 3 2 2Night (8.00 p.m.–8.00 a.m.) 2 1 0
Table 4.2
Number of seats available for medical visit and treatment of red-tagged patientsin RA and in the dedicated areas of MU and SU.
MU SU RADay (8.00 a.m.–8.00 p.m.) 1 2 2Night (8.00 p.m.–8.00 a.m.) 1 2 2
Table 4.3
Feasible assignments of patients to the ED units according to the color tag. A crossat the entry ( i, j ) indicates that a patient with color tag i can be assigned to the unit j . MU SU RA MIU OU
White - - - X X
Green
X X - X X
Yellow
X X - - X
Red
X X X - -
As regards the patient flow, after arriving autonomously or by emergency med-ical vehicles, all the incoming patients are admitted to the triage area, where anurse assigns the color tag. After the triage, the patients are visited and treatedin one of the units previously described, according to the color tag assigned andthe severity of the illness/injury. Table 4.3 represents a scheme showing the unitswhere patients may be assigned based on the color tags. In case of red-tagged pa-tients, the medical visit is timely performed in RA or in the dedicated areas of MUand SU. As concerns the other color tags, the yellow and green tagged patientsshare the same resources in MU and SU. However, while the former patients canbe assigned only to MU and SU, the latter may be sent to MIU during its open-ing hours if their health conditions are deemed as not likely to worsen. If MIU isclosed, all the green-tagged patients are visited and treated in MU and SU. Thisdiversion allows the ED to reduce the occurrence of work overload in SU and MU,which may give rise to overcrowded units. Moreover, it is important to point outthat the white triage tag is assigned only if MIU is open, otherwise the green tag is used.In many cases, after the medical visit, additional exams may be required. Otherthan performing reassessments of the patients and requiring additional exams,physicians may also require further observation periods. Finally, at the end of thepathway, patients are discharged from the ED. This last stage includes a final waiting time whose duration depends on the type of outcome. Indeed, a longerwait is expected for patients that need to either be hospitalized at a hospital wardor transferred to another hospital, while patients discharged home can usuallyleave the ED in shorter time.4.2 The Discrete Event Simulation modelThis section describes the DES model of the ED of Policlinico Umberto I, which hasbeen implemented by using Simmer [33], a process-oriented and trajectory-basedDES package for R . In the DES model, patients are represented by the model entities , which are created according to the statistical model used for modelingthe arrival process. After arriving at the ED, the simulated patient flow followedby each entity is based on trajectories determined by the logical rules used tobuild the model. Each trajectory is associated with the pathway followed by apatient according both to the color tag received at triage and to the ED unitassigned. Figures 4.1–4.5 show the patient flow from the arrival to the dischargeaccording to the color tag assigned at triage. In particular, Figure 4.1 focuses onthe logic underlying the assignment of color tag at triage. The other figures dealwith the segments of the patient flow following the triage, providing in bracketsthe information about the resources seized in case the corresponding process is nota simple delay. Such resources represent the seats at each ED unit, whose capacitycan be either fixed, as is the case with RA and the areas of MU and SU dedicatedto red patients, or based on schedule, as for MU, SU, and MIU.The patient flow through the model can be described as follows. After beingcreated at the beginning of the simulation, the entities corresponding to deceasedpatients leave the model, while the remaining entities enter the triage area for theassignment of both color tag and ED unit, which are stored as entity attributes.The discrete uniform probability distributions used for assigning the color tag andthe ED unit vary according to the time at which the entity starts the triage. Ifthis event happens in the daytime, the corresponding probability distributionsinclude also the white color tag and MIU among the alternatives to sample. Thetriage phase is represented by as many delay processes as the number of colortags and units. Each of these delays is associated with a probability distributionthat returns the value of the triage duration for each patient. In particular, the8 different processes considered correspond to the possible pairs of color tags andunits reported in Table 4.3, with the exception of the pair given by the green colortag and MIU. Indeed, the triage duration for a green-tagged patient eligible forMIU depends on the specialty required, whether medical or surgical. Since thisinformation is unknown, in the DES model half of the MIU patients are subjectto the triage duration of MU patients, while the other half are subject to thetriage duration of SU patients. It is important to remark that, since the processesare simple delay, no queue is generated before the triage. The reason underlyingthis choice is twofold: on the one hand, this allows for alignment between data e simulation model, since the timestamps denoted as t in Figure 2.1 are used toderive the arrival probability distribution, which consequently returns the startingtime of triage (and not the starting time of the waiting time for triage); on theother hand, the waiting time before triage can be considered negligible, accordingto the ED staff. imulation-based optimization for ED DES model calibration 11 Fig. 4.1
Plot of the patient flow in the ED from the arrival to the assignment of the colortag.
After the triage process, some entities leave the model, thus representing thepatients who leave without being seen, while the other ones start waiting for themedical visit. The entity selected for the visit depends both on the priority class,i.e., the color tag, and on the FIFO criterion. Then, when the visit begins, oneseat in the corresponding unit is seized and the duration of the visit is returnedby a suitable probability distribution based on the entity color tag.After the medical visit, the following phase includes exams and reassessments ,whose duration is generated by means of a suitable probability distribution. Likethe triage, a single delay process is used for this phase as well, due to the lack ofknowledge of the resources required. It is important to point out that all possiblefurther treatments are considered included in the service time of this process.After the exams, the entities corresponding to patients refusing hospitalization , leaving during exams , and being transferred to another hospital are removed fromthe model, while the other entities proceed to the next stage, i.e., the final waitingtime before discharge.Since in the dataset some timestamps defining the starting and ending time oftriage, medical visit, and exams are not available, the service times correspondingto these activities cannot be recovered directly. Moreover, additional delay pro-cesses are considered to reproduce all the service times that cannot be directlycomputed through the data, such as the setup times that are sometimes neededfor sanitizing the ED areas and the idle times caused by unexpected requests forpersonnel from other ED units, which give rise to sudden activity disruption. Inorder to effectively model these times, which impact on all the patients withoutpredictable patterns, uniform probability distributions are considered before thequeue for the medical visit, according to the color tag.As regards the KPIs of interest, the focus is on the time differences DOT and
DIT described in Figure 2.1, which can be computed through the data of thereal system and then compared with the corresponding values returned by thesimulation. Therefore, the KPIs of interest are – DOT , which is the time difference between the starting time of the triage andthe starting time of the visit, namely t − t in Figure 2.1; – DIT , which is the time difference between the starting time of the visit andthe time of the discharge, namely t − t in Figure 2.1. Fig. 4.2
Plot of the red-tagged patient flow.
Fig. 4.3
Plot of the yellow taggedpatient flow.
Fig. 4.4
Plot of the white tagged pa-tient flow.imulation-based optimization for ED DES model calibration 13
Fig. 4.5
Plot of the green-tagged patient flow.
Fig. 4.6
Plot of the weekly average hourly arrival rate for the first 13 weeks of the year. accordance with the literature (see, i.e., [22]), the average arrival rate among thedays of the week is significantly different. Therefore, since averaging over thesedays would lead to inaccurate results, the different days of the week must beconsidered separately.
The stochastic process defined by the arrivals to the ED is assumed to be wellmodeled by a Nonhomogeneous Poisson Process (NHPP), which is a standard as-sumption in the literature (see, e.g.,[35,24,39,2,1,17]). The validity of this hypoth-esis for the case study considered is supported by the extensive analysis performedin [9,10], where statistical hypothesis testing is adopted. To achieve an accuraterepresentation of the arrival rate, 24 time slots of unitary length are considered foreach day of the week, thus obtaining a piecewise constant approximation. Whilethis approach allows taking into account the within-day variation, the day-to-dayvariation is considered by estimating the hourly arrival rate function separatelyfor each day of the week.Although the collected data concerns the whole year 2018, the DES modelof the ED has been built by using the data related to January. The choice offocusing on this month stems from the will of the ED management to reduce theovercrowding level observed in the winter season, which exhibits longer waitingtimes, as emerged by interviewing the ED staff. Among the winter months, Januaryis observed to suffer from the heaviest workload, which puts a strain on the EDprocesses, thus requiring a careful analysis. However, the simulation model couldbe also easily adapted to include input parameters estimated from data related todifferent months.From 00:00 of January 1 to 23:59 of January 31, 2018, the total number ofpatients arrived to the ED is 4192. The timestamps recorded are the ones markedas known in Figure 2.1. In Table 4.4, the number and the percentage of colortags assigned at triage are reported along with the number of patients who leavewithout being seen (LWBS). Since from Monday to Saturday MIU is open from8.00 a.m. to 8.00 p.m., the daytime and the night are analyzed separately to reflectthe change in the ED setting. Although, in some cases, re-evaluations can lead todifferent patient tags at discharge with respect to those assigned at triage, in thedataset considered this information is unavailable. Finally, Tables 4.5–4.6 show
Table 4.4
Number and percentage of color tags assigned at triage in the daytime (columns2-3) and at night (column 4-5).
TAGS ASSIGNED AT TRIAGEDay (8.00 a.m.–8.00 p.m.) Night (8.00 p.m.–8.00 a.m.) LWBS
White
65 2.17% - - 1
Green
Yellow
Red
157 5.33% 95 7.71 % -2948 100% 1231 100% 17 the number and proportion of the color tags among the ED units. Although OUis out of the scope of this analysis, Table 4.5 includes this unit since OU patientsshare the triage station with the other patients, thus affecting the counts reportedin Table 4.4.
The percentage of patients requiring more than one medical visit is reportedin Table 4.7 for each color tag and ED unit. Since the number of patients needingmore than visit is small, in the DES model all the entities are assumed to undergoone medical visit and, accordingly, each patient flow is represented sequentially.The probability distributions used to generate values for the final waiting time imulation-based optimization for ED DES model calibration 15
Table 4.5
Number (and percentage) of patients assigned to the ED units for each color tag.
MU SU RA MIU OU
White - - - 47 (73.44 %) 17 (26.56 %)
Green
248 (13.75 %) 628 (34.81 %) - 260 (14.41 %) 668 (37.03 %)
Yellow
Red
191 (75.79 %) 45 (17.86 %) 16 (6.35 %) - -1755 1366 16 307 734
Table 4.6
Number (and percentage) of green-tagged patients assigned to MU, SU, and MIUin the daytime (8.00 a.m. – 8.00 p.m.) and at night (8.00 p.m. – 8.00 a.m.).
MU SU MIU
Daytime
132 (16.60 %) 403 (50.69 %) 260 (32.70 %)
Night
116 (34.02 %) 225 (65.98 %) -248 628 260
Table 4.7
Percentage of patients requiring more than one visit and probability distributionof the final waiting time before discharge for each color tag and ED unit.
MORE THAN ONE VISIT FINAL WAITING TIME
White
MIU 0% Weib(1.12, 0.41)MIU 1.61% Weib(0.83, 0.57)
Green
SU 2.44% Weib(0.61, 1.07)MU 4.8% Weib(0.57, 4.91)
Yellow
SU 2.04% Weib(0.63, 2.33)MU 3.86% Weib(0.62, 7.22)RA 6.25% Weib(0.67, 9.26)
Red
SU 6.52% Weib(0.71, 4.47)MU 7.9% Weib(0.69, 6.65) before the discharge are reported in Table 4.7 as well (note that the base unit fortime adopted in the DES model is the hour).4.4 Model verification/validation and design of experimentsThe aim of the DES model of the ED of Policlinico Umberto I is to provide areliable tool for assessing the impact of changes to the current settings on theovercrowding level. To this end, verification and validation are important steps tobe carefully considered. While standard techniques have been adopted to verifythe model, such as debugging and model trace, validation has required a moresophisticated analysis. Indeed, a model calibration approach has been used toachieve an accurate simulation output that well reproduces the real system data.Important KPIs to describe the ED status are the two time differences con-tained in T and the number of entities generated by the simulation and associated with each color tag. While the former KPIs are part of the objective function usedin the calibration procedure and they are calculated across the iterations of theoptimization algorithm, the latter do not require a continuous monitoring, sincetheir values are not affected by the calibration. Indeed, the only variations are dueto the different number of entities discharged before the simulation ends, which leads to different results in terms of KPI computation. Table 4.8 reports the val-ues of the patient counts obtained from the simulation of the starting point of thecalibration procedure along with their confidence intervals. By comparing thesevalues with the counts in Table 4.5, it can be observed that the simulation modelprovides an accurate output in terms of number of patients. Table 4.8
Output values (with confidence interval) for the number of patients returned bythe simulation of the starting point of the calibration procedure.
MU SU MIU
White - - ± . Green ± .
33 612 ± .
73 239 ± . Yellow ± .
55 716 ± . - Red ± .
39 57 ± . - Since the objective and constraint functions compare the real system data andthe simulation output, achieving adequate accuracy is imperative for a fair com-parison. To this end, the simulation output is estimated through 30 independentsimulation replications, each of them 38 days long, with a warm up period of 7days, thus matching the 31 days of January.
The sets defined in Section 3 can be adapted to the case study as follows. – Let C = { W,G,Y,R } be the set of the triage tags. – Let U ( c ) ⊆ { MU , SU , MIU } be the set of the ED units where patients with tag c ∈ C can be visited and treated. In particular, – U ( W ) = { MIU } , – U ( G ) = { MU , SU , MIU } , – U ( Y ) = { MU , SU } , – U ( R ) = { MU , SU } .For the specific instance represented by this case study, 8 different pairs of param-eters are considered for the Weibull distributions representing the service timesof medical visit and exams. Each pair is associated with an element of the sets U ( c ). Instead, 7 pairs are used for the probability distributions of triage for thereasons stated in Section 4.2, thus resulting in 23 overall pairs of parameters, i.e.,46 unknown parameters to be determined through the calibration procedure.A crucial and preliminary step is the choice of the starting point of the opti-mization, given by x , y , and z (where x , y , and z are the vectors containingall the corresponding shape and scale parameters, as defined in Section 3), whichappears not to be straightforward since the missing data prevents the computa-tion of good initial values for the parameters. Instead of starting from randomly generated values, a better strategy is to leverage the known information in a rea-sonable manner. While the parameters x c u of the triage probability distributionsare arbitrarily fixed so that the generated service times are in accordance withthe ED staff suggestions, for the parameters of the medical visit and exams dis-tributions a different approach is followed. Indeed, in addition to the timestamps imulation-based optimization for ED DES model calibration 17 shown in Figure 2.1, which are the most commonly known timestamps in everyED, for this specific case study further available information is represented bythe time at which the physician requires exams. In general, this happens duringthe medical visit, although sometimes exams are required throughout the patientflow when periodic check-ups are performed. To provide the parameters y c u of themedical visit probability distributions with initial values, a good starting pointcan be obtained by computing for each patient the duration between the startof the visit and the latest time at which an exam is required. The latter time isconsidered within a reasonable period from the start of the visit that should reflectthe medical visit’s maximum service time. Once this timestamp is identified andthe durations are consequently available, their values can be used to initialize theparameters by fitting the corresponding Weibull probability distributions. More-over, the time between the presumptive end of the visit and the discharge can beused in the same manner to obtain good initial values for the parameters z c u ofthe exams distributions.Although Problem 3.1 is a continuous optimization problem, its variables areconsidered as discrete to efficiently solve the problem. In particular, at the time ofsolving the problem, x c up , y c up , and z c up , with p ∈ { , } , are treated as granularvariables, i.e., variables with a controlled number of decimals (see, e.g., [5]). Thischoice is motivated by the fact that the output of the simulation model is insensi-tive to small changes in the values of the decision variables, thus making discretevariables preferable. To this end, a specific granularity is considered based on therole of the parameter associated with the variable in the corresponding probabilitydistribution. All the variables are pairs of shape and scale parameters, which havea different impact on the simulation output. Indeed, while the former parametersdetermine the shape of the probability distribution, the latter affect the scale ofthe values generated. With respect to the Weibull distribution adopted, in prac-tice the behavior observed is that the larger the value of the shape parameter, thelarger the dispersion of the values of the random variables associated. Moreover,such random variables take on values with a larger scale as the value of the scaleparameter increases. Since some preliminary analyses have shown that the impacton the simulation output, measured in terms of the two KPIs DOT and
DIT , ismore noticeable for variations in the scale parameter, for each pair of variablesthe granularity δ minp of the values is fixed to 10 − if p = 1 (i.e., shape parameter)and to 10 − if p = 2 (i.e., scale parameter). Note that treating the variables asgranular requires the new variables to be equal to x c up /δ minp ∈ Z , thus meaningthat x c up have to take on real values that are multiple of δ minp . The same type ofconstraint applies to the variables of the visit and exams probability distributions.To solve the SBO problem for calibrating the simulation model of the EDin hand, the approach of the Sample Average Approximation (SAA) [29,21] isadopted. As a consequence, the resulting optimization problem is deterministicand it can be solved by applying an algorithm from the class of DFO [8,4]. Inparticular, the optimization algorithm proposed in [27] is used for solving theproblem by adopting its default parameters. This algorithm has been successfully applied to the same case study in [10] and it is suited for efficiently solving integerblack-box constrained problems, like Problem 3.1, thanks to unconventional searchdirections, an effective penalty approach, and strong global convergence proper-ties. The maximum number of function evaluations, which represents the stoppingcondition, is set to 3000. Note that by using the SAA approach, the empirical cumulative distribution functions used in the optimization problem are estimatedthrough the corresponding sample means over the 30 independent simulation repli-cations.As concerns the optimization problem to solve, the values of the lower andupper bounds introduced in Problem (3.1) are reported in Table 5.1. Moreover,the tolerance τ c u iµ and τ c u iσ are both fixed to 0 .
35 for yellow and green-taggedpatients visited and treated in MU and SU, 0 . Table 5.1
Values of the lower and upper bounds for each pair of variables ( x c u , x c u ),( y c u , y c u ), and ( z c u , z c u ), for all c ∈ C and u ∈ U ( c ). l x c u l x c u l y c u l y c u l z c u l z c u u x c u u x c u u y c u u y c u u z c u u z c u Given these settings, the experimental results are obtained by using a PC withIntel Core i7-4790K quad-core 4.00 GHz Processor and 32 GB RAM. Table 5.2reports the service time distributions of triage, medical visit, and exams and re-assessments at the final solution determined by the optimization algorithm (werecall that the base unit for time adopted in the DES model is the hour). Twotypes of plots are shown in Appendices A and B to assess the accuracy of thecalibrated simulation model. The former appendix focuses on the KPI given bythe average hourly number of patients inside the ED after the start of the medi-cal visit, namely, patients satisfying two conditions: ( i ) they have already startedthe medical visit; ( ii ) they have not been discharged yet. In particular, the plotsreport the comparison between the value of the KPI computed through the realdata and the corresponding value derived from the simulation output along withits 95% confidence interval. The latter appendix reports the comparison betweenthe Empirical Cumulative Distribution Functions (ECDFs) of the values of thetime differences DOT and
DIT computed for each patient through the real dataand the simulation output resulting from the model calibration procedure. In bothappendices, the comparisons are performed for all the color tags c ∈ C and forall the ED units u ∈ U ( c ) where a patient with tag c can be assigned. Moreover,both the average hourly number of patients inside the ED after the start of themedical visit and the ECDFs are obtained as an average over the independentsimulation replications. Both types of plots show that, in general, the values ofthe KPIs corresponding to the system data are well reproduced by the simulationoutput. However, some dissimilarities may be still present in spite of the calibra-tion, especially for green and yellow-tagged patients. This is due to the difficultyin reproducing the patient flow when sharing of the resources is involved betweenpatients with different triage tag, as is the case for the medical visit of green andyellow-tagged patients at both M U and SU . Moreover, the experimental results show that, in general, the comparison between real data and simulation output ismore accurate for the time differences DIT. Indeed, larger errors are observed inboth types of plots when time differences DOT are considered due to the difficultyin reproducing the actual waiting times of patients, which depend both on theservice time of the medical visit and on the number of patients that are in queue. imulation-based optimization for ED DES model calibration 19 Table 5.2
Service time distributions of triage, medical visit, and exams and reassessments atthe final solution determined by the optimization algorithm. Note that the triage service timeof green-tagged patients in MIU is missing since in the DES model these patients are subjectto the Weibull distributions used for SU and MU.
TRIAGE MEDICAL VISIT EXAMS
White
MIU Weib(0.61, 0.42) Weib(0.99, 0.57) Weib(0.62, 0.23)MIU - Weib(1.06, 0.41) Weib(1.29, 1.61)
Green
SU Weib(1000, 0.50) Weib(0.64, 0.31) Weib(0.68, 2.11)MU Weib(0.60, 0.49) Weib(1000, 1.63) Weib(0.69, 12.26)
Yellow
SU Weib(1000, 0.28) Weib(0.51, 0.22) Weib(0.72, 6.32)MU Weib(0.55, 0.50) Weib(0.62, 7.22) Weib(0.78, 25.46)
Red
SU Weib(0.92, 0.10) Weib(0.88, 0.75) Weib(1.41, 19.34)MU Weib(0.80, 0.19) Weib(2.24, 0.36) Weib(0.78, 30.92)
Hence, it is the result of a seize-delay-release process used within the simulation,while DIT is given by the mere sum of service times.By observing the ECDF plots in Appendix B, it is possible to gain useful insightthat histograms, which are commonly adopted in the literature, fail to provide. Forinstance, the ECDF related to yellow-tagged patients in SU shows that the timedifferences DOT collected from the simulation model are systematically larger thanthe corresponding real values, meaning that either the estimated parameters of theprobability distribution of the medical visit service time or the modeling of the in-teraction with green-tagged patients in SU may be improved with further research.To avoid ambiguities, it is important to remark that since the simulated ECDFsare obtained as an average over independent simulation replications, the jumps(or steps) in these simulated functions correspond to patients observed across allthe replications. This is a consequence of averaging different step functions, suchas the ECDFs, which have jumps at different points in their domain.The same conclusions on the accuracy of the simulation model can be drawnalso by observing Figures 5.1–5.2, which report the comparison between the aver-age current and simulated values of DOT and DIT for all the triage tags and allthe ED units. As already mentioned, the simulation output related to the time dif-ferences DIT is a good approximation of the real system values since for each colortag the current values are either within the corresponding confidence interval orclose to it. Instead, even though for the time differences DOT the accuracy of thesimulation responses as estimates of the current values is lower, the relative erroris within the tolerances τ c u iµ chosen. Indeed, all the constraints of Problem (3.1)are (of course) satisfied at the optimal solution obtained by applying the modelcalibration approach. Note that reducing the values of the tolerances potentiallyallows achieving more accurate solutions, however, the optimization algorithm isnot guaranteed to determine a feasible solution.The experimental results discussed above show that the model calibration pro-cedure enables estimating the missing parameters in order for the simulation model to satisfactorily reproduce the ED operations despite the unavailable timestamps,which represent the main hurdle for achieving a high level of accuracy. The con-straints on the sample means used in Problem (3.1) ensure that the responsesof the simulation model are on average close to the corresponding real values.Moreover, the constraints on the sample variances of the KPIs help the optimiza- Fig. 5.1
Plot of current values and simulation output of DOT with the confidence interval.
Fig. 5.2
Plot of current values and simulation output of DIT with the confidence interval. tion algorithm to avoid points associated with inaccurate results. However, whenlooking more closely into the details of the numerical results, some dissimilaritiesbetween real data and simulation output may be observed as a result of the impactof the problem of missing data, which necessarily undermines the overall accuracy.Notwithstanding, the simulation model may be deemed sufficiently reliable withrespect to the specific objectives of the analysis.
The SBO approach proposed in this paper addresses one of the data quality prob-lems that frequently affect the DES models concerning EDs. In particular, thefocus is on the problem of missing data, which consists in the unavailability of data related to some of the starting and ending times of the activities performedin the ED. This well-known issue is responsible for the lack of knowledge of theservice time of some processes, which is required to estimate the correspondingprobability distributions to use in the simulation model. For this purpose, afterassuming suitable probability distributions for the processes associated with the imulation-based optimization for ED DES model calibration 21 unknown service times, a model calibration procedure is proposed to estimate themissing parameters of such distributions.The preliminary experimental results on the ED of a big hospital in Italy showthat the model calibration procedure enables estimating the missing parametersin order for the simulation model to satisfactorily reproduce the ED operationsdespite the unavailable timestamps, which represent the main hurdle for achievinga high level of accuracy. Although proper constraints ensure that the responses ofthe simulation model are on average close to the corresponding real values, whenlooking more closely into the details of the experimental results, some dissimilar-ities between real data and simulation output may be observed as a result of theimpact of the problem of missing data, which necessarily undermines the over-all accuracy. Notwithstanding, the simulation model may be deemed sufficientlyreliable with respect to the specific objectives of the analysis. Since there is stilla margin for improvement, several ideas are considered for further research. Forinstance, as the worst results are observed for the time differences DOT, a signifi-cant improvement may be obtained by using different weights for the terms of theobjective function of Problem (3.1) and by assigning larger values to the termsassociated with DOT. Moreover, different objective functions and starting pointsmay be taken into account also to assess the robustness of the approach. Finally,increasing the number of function evaluations used as stopping condition may leadto better solutions allowing a more thorough exploration of the feasible region.
Acknowledgment
The authors are grateful to Prof. F. Romano (Chief Medical Officer) and Dr. L.De Vito (Medical Director of ED) of
Policlinico Umberto I of Rome for theiravailability to carry out this study.
A Model calibration - Plots I
This appendix focuses on the KPI given by the average hourly number of patients satisfying twoconditions: ( i ) they have already started the medical visit; ( ii ) they have not been dischargedyet. In particular, the plots report the comparison between the value of the KPI computedthrough the real data and the corresponding value derived from the simulation output alongwith its 95% confidence interval. The comparisons are performed for all color tags c ∈ C andfor all the ED units u ∈ U ( c ) where a patient with color tag c can be assigned.2 A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Romaimulation-based optimization for ED DES model calibration 23 B Model calibration - Plots II
This appendix reports the comparison between the Empirical Cumulative Distribution Func-tions (ECDFs) of the time differences
DOT and
DIT computed through the real data and thesimulation output resulting from the calibration procedure. The comparisons are performedfor all color tags c ∈ C and for all ED units u ∈ U ( c ). The colored curves correspond to thesimulation output.4 A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Romaimulation-based optimization for ED DES model calibration 256 A. De Santis, T. Giovannelli, S. Lucidi, M. Messedaglia, M. Roma Conflict of interest
The authors declare that they have no conflict of interest.
References
1. Ahalt, V., Argon, N., Strickler, J., Mehrotra, A.: Comparison of emergency departmentcrowding scores: a discrete-event simulation approach. Health Care Management Science , 144–155 (2018)2. Ahmed, M.A., Alkhamis, T.M.: Simulation optimization for an emergency departmenthealthcare unit in Kuwait. European Journal of Operational Research (3), 936 – 942(2009)3. Amaran, S., Sahinidis, N., Sharda, B., Bury, S.: Simulation optimization: a review ofalgorithms and applications. 4OR , 301–333 (2014)4. Audet, C., Hare, W.: Derivative-Free and Blackbox Optimization. Springer InternationalPublishing (2017). DOI 10.1007/978-3-319-68913-55. Audet, C., Le Digabel, S., Tribes, C.: The mesh adaptive direct search algorithm forgranular and discrete variables. SIAM Journal on Optimization , 1164–1189 (2019).DOI 10.1137/18M1175872imulation-based optimization for ED DES model calibration 276. Baumlin, K., Shapiro, J., Weiner, C., Gottlieb, B., Chawla, N., Richardson, L.: Clinicalinformation system and process redesign improves emergency department efficiency. TheJoint Commission Journal on Quality and Patient Safety (4), 179–185 (2010)7. Bedoya-Valencia, L., Kirac, E.: Evaluating alternative resource allocation in an emergencydepartment using discrete event simulation. Simulation , 1041–1051 (2016)8. Conn, A., Scheinber, K., Vicente, L.: Introduction to Derivative-Free Optimization. Societyfor Industrial and Applied Mathematics, Philadelphia (2009)9. De Santis, A., Giovannelli, T., Lucidi, S., Messedaglia, M., Roma, M.: An optimalnon–uniform piecewise constant approximation for the patient arrival rate for a moreefficient representation of the emergency departments arrival process. Technical Report1–2020, Dipartimento di Ingegneria Informatica Automatica e Gestionale “A. Ruberti”,SAPIENZA Universit`a di Roma (2020)10. De Santis, A., Giovannelli, T., Lucidi, S., Messedaglia, M., Roma, M.: Determining theoptimal piecewise constant approximation for the Nonhomogeneous Poisson Process rateof Emergency Department patient arrivals. arXiv e-prints arXiv:2101.11138 (2021)11. Di Somma, S., Paladino, L., Vaughan, L., Lalle, I., Magrini, L., Magnanti, M.: Overcrowd-ing in emergency department: an international issue. Internal and Emergency Medicine (2014). DOI 10.1007/s11739-014-1154-812. Duma, D., Aringhieri, R.: An ad hoc process mining approach to discover patient pathsof an emergency department. Flexible Services and Manufacturing Journal (2018).DOI 10.1007/s10696-018-9330-113. Fava, G., Giovannelli, T., Messedaglia, M., Roma, M.: Effect of different patient peakarrivals on an Emergency Department via discrete event simulation. arXiv e-printsarXiv:2101.12432 (2021)14. Fu, M.C. (ed.): Handbook of Simulation Optimization. Springer, New York (2015)15. Gosavi, A.: Simulation-Based Optimization: Parametric Optimization Techniques and Re-inforcement Learning, 2nd edn. Springer Publishing Company, Incorporated (2014)16. Gray, G.A., Kolda, T.G.: Algorithm 856: Appspack 4.0: Asynchronous parallel patternsearch for derivative-free optimization. ACM Trans. Math. Softw. (3), 485–507 (2006).DOI 10.1145/1163641.1163647. URL https://doi.org/10.1145/1163641.1163647
17. Guo, H., Gao, S., Tsui, K., Niu, T.: Simulation optimization for medical staff configurationat emergency department in Hong Kong. IEEE Transactions on Automation Science andEngineering (4), 1655–1665 (2017)18. Guo, H., Goldsman, D., Tsui, K.L., Zhou, Y., Wong, S.Y.: Using simulation and optimi-sation to characterise durations of emergency department service times with incompletedata. International Journal of Production Research (21), 6494–6511 (2016)19. Hoot, N., Aronsky, D.: Systematic review of emergency department crowding: causes,effects, and solutions. Annals of Emergency Medicine (2), 126–136 (2008)20. Joshi, V., Lim, C., Teng, S.G.: Simulation study: Improvement for non-urgent patientprocesses in the emergency department. Engineering Management Journal , 145–157(2016)21. Kim, S., Pasupathy, R., Henderson, S.G.: A Guide to Sample Average Approximation, pp.207–243. Springer New York, New York, NY (2015). DOI 10.1007/978-1-4939-1384-8 8.URL https://doi.org/10.1007/978-1-4939-1384-8_8
22. Kim, S.H., Whitt, W.: Are call center and hospital arrivals well modeled by nonhomo-geneous Poisson process ? Manufactory & Service Operations Management , 464–480(2014)23. Kramer, A., Dosi, C., Iori, M., Vignoli, M.: Successful implementation of discrete eventsimulation: the case of an italian emergency department (2020)24. Kuo, Y.H., Rado, O., Lupia, B., Leung, J.M.Y., Graham, C.A.: Improving the efficiency ofa hospital emergency department: a simulation study with indirectly imputed service-timedistributions. Flexible Services and Manufacturing Journal (1), 120–147 (2016)25. Law, A.: Simulation Modeling and Analysis, fifth edn. McGraw-Hill (2015)26. Liu, Z., Rexachs, D., Epelde, F., Luque, E.: A simulation and optimization based methodfor calibrating agent-based emergency department models under data scarcity. Comput.Ind. Eng. (C), 300–309 (2017). DOI 10.1016/j.cie.2016.11.036. URL https://doi.org/10.1016/j.cie.2016.11.036
27. Liuzzi, G., Lucidi, S., Rinaldi, F.: An algorithmic framework based on primitive di-rections and nonmonotone line searches for black-box optimization problems with inte-ger variables. Mathematical Programming Computation (4), 673–702 (2020). DOI10.1007/s12532-020-00182-7. URL https://doi.org/10.1007/s12532-020-00182-7 , 368–376 (2017)29. Pagnoncelli, B.K., Ahmed, S., Shapiro, A.: Sample Average Approximation Method forChance Constrained Programming: Theory and Applications. Journal of OptimizationTheory and Applications (2), 399–416 (2009). DOI 10.1007/s10957-009-9523-630. Paul, S.A., Reddy, M.C., De Flitch, C.J.: A systematic review of simulation studies inves-tigating emergency department overcrowding. Simulation (8-9), 559–571 (2010)31. Pines, J.M., Iyer, S., Disbot, M., Hollander, J.E., Shofer, F.S., Datner, E.M.: The effect ofemergency department crowding on patient satisfaction for admitted patients. AcademicEmergency Medicine (9), 825–831 (2008). DOI https://doi.org/10.1111/j.1553-2712.2008.00200.x32. Salmon, A., Rachuba, S., Briscoe, S., Pitt, M.: A structured literature review of simula-tion modelling applied to emergency departments: Current patterns and emerging trends.Operations Research for Health Care , 1–13 (2018)33. Ucar, I., Smeets, B., Azcorra, A.: simmer: Discrete-event simulation for R. Journal ofStatistical Software (2), 1–30 (2019). DOI 10.18637/jss.v090.i0234. Vanbrabant, L., Martin, N., Ramaekers, K., Braekers, K.: Quality of input data in emer-gency department simulations: Framework and assessment techniques. Simulation Mod-elling Practice and Theory , 83 – 101 (2019)35. Whitt, W., Zhang, X.: A data-driven model of an emergency department. OperationsResearch for Health Care , 1 – 15 (2017)36. Wiler, J., Griffey, R., Olsen, T.: Review of modeling approaches for emergency departmentpatient flow and crowding research. Academic Emergency Medicine , 1371–1379 (2011)37. Wong, Z.S.Y., Lit, A.C.H., Leung, S.Y., Tsui, K.L., Chin, K.S.: A discrete-event simulationstudy for emergency room capacity management in a hong kong hospital. 2016 WinterSimulation Conference (WSC) pp. 1970–1981 (2016)38. Yousefi, M., Yousefi, M., Fogliatto, F.S.: Simulation-based optimization methods appliedin hospital emergency departments: A systematic review. SIMULATION (10), 791–806(2020). DOI 10.1177/003754972094448339. Zeinali, F., Mahootchi, M., Sepehri, M.: Resource planning in the emergency departments:a simulation-base metamodeling approach. Simulation Modelling Practice and Theory53