key: cord-0046632-f20cavib
authors: Ballarini, Paolo; Duma, Davide; Horváth, Andras; Aringhieri, Roberto
title: Petri Nets Validation of Markovian Models of Emergency Department Arrivals
date: 2020-06-02
journal: Application and Theory of Petri Nets and Concurrency
DOI: 10.1007/978-3-030-51831-8_11
sha: 9e67b608ac2296e0cc6898f274e967f9e64a7814
doc_id: 46632
cord_uid: f20cavib

Modeling of hospital’s Emergency Departments (ED) is vital for optimisation of health services offered to patients that shows up at an ED requiring treatments with different level of emergency. In this paper we present a modeling study whose contribution is twofold: first, based on a dataset relative to the ED of an Italian hospital, we derive different kinds of Markovian models capable to reproduce, at different extents, the statistical character of dataset arrivals; second, we validate the derived arrivals model by interfacing it with a Petri net model of the services an ED patient undergoes. The empirical assessment of a few key performance indicators allowed us to validate some of the derived arrival process model, thus confirming that they can be used for predicting the performance of an ED.

Optimisation of hospital's Emergency Departments (ED) is concerned with maximising the efficiency of health services offered to patients that shows up at an ED requiring treatments with different level of emergency. In this respect the ability to build formal models that faithfully reproduce the behaviour of an ED and to analyse their performances is vital. A faithful model of an ED must consist of at least two components: a model of the patients arrival and a model of the ED services (the various activities that an ED patient may undergo throughout his permanence at the ED). The composition of these two models must be such that the model of patients arrival is used to feed in the services model and the analysis of the performances of the composed model, through meaningful key performance indicators (KPIs), may yield essential indications as to where intervene so to improve the overall efficiency of the ED.

In this paper we focus on the problem of deriving a stochastic model of the patient arrivals that is capable of correctly reproducing the statistical character of a given arrivals dataset. To this aim we considered a real dataset relative to the activity of the ED of an hospital of the city of Cantù in northern Italy and which contains different data including the complete set of timestamps of the arrival of each patients over one year (i.e. the year 2015). Based on the statistical analysis of the dataset we derived different instances of Markovian models belonging to different classes (i.e. Markov renewal processes, Markov arrival processes, hidden Markov models) and we showed that only one model instance, amongst those derived, is capable of fully reproducing the statistical character of the dataset, including the periodicity of patients arrivals exhibited by the dataset. To further validate the derived models we developed a formal encoding in terms of stochastic Petri nets. This allowed us to compare the performances of each arrival models, coupled with a model of ED services, through assessment of a KPI formally encoded and accurately assessed through a statistical model checking tool. The obtained results evidenced that one amongst the derived arrival models faithfully match the dataset behaviour hence it can reliably be used for analysing the performances of different ED services.

Paper Organisation. The paper is organised as follows: In Sect. 2 we present the statistical analysis of the dataset and discuss the main statistical characteristics of the patient's arrival process by considering different perspectives. In Sect. 3 we describe the derivation of different kinds of Markovian models of the patient arrival's process, specifically we discuss the class of renewal Markov processes given in terms of continuous phase type (CPH) distributions, in Sect. 3.1, the class of Markov arrival process (MAP), in Sect. 3.2 and finally the class of hidden Markov models (HMM), in Sect. 3.3. In Sect. 4 we give the formal encoding of each of the derived models in terms of stochastic Petri nets. Finally in Sect. 5 we present the results of some experiments aimed at assessing a KPI against the model given by the composition of patient arrival models and a simple model of the ED services. Conclusive remarks wrap up the paper in Sect. 6. Related Work. Literature on modeling of hospital EDs is very wide, faces different kinds of problems and cannot be reviewed in brief. We mention only some works to position ourselves in the literature and point out differences. Many works, such as [2, 3, 12, 14, 21, 24] , deal with forecasting the number of patients arrivals but do not experiment with Markovian models. The models proposed in this paper can also be applied to forecasting but we cannot elaborate on this issue due to space limits. Several works, e.g., [20, 22] , build a queuing model of ED but do not deal with data driven modeling of the arrival process. Few works, e.g., [31] , are concerned with the data-driven development of stochastic models capturing both the patient arrival as well as the services received at the ED. Our paper falls in this category but with respect to other we experiment with a wide range of Markovian processes to model patient arrivals and provide also their Petri net counterpart which can be derived in an automatic manner.

One way of representing the patient arrival flow is to register the time elapsed between consecutive show ups. To this end let X i , i = 1, 2, 3, ... denote the ith inter-arrival time. Another possibility is to count the number of patients that arrive in an interval of a certain length. We will use intervals whose length is one hour, one day or one week. The associated series are Y H,i , Y D,i and Y W,i with i = 1, 2, 3, ... which denote the number of patient arrivals in the ith hour, day and week, respectively. Clearly, the second approach keeps less information with respect to the first but it has the advantage that it relates easily to the periodic nature of the patient arrival process. For example, the series Y H,3+24j with j = 0, 1, 2, ... provides the number of arrivals between 2am and 3am during the 1st, the 2nd, the 3rd day, and so on.

We start off by analyzing the patient arrival pattern inside a day by considering the average and the variance of the number of patient arrivals during the ith hour of the day. I.e., Fig. 1 , E H,i and V ar H,i are depicted as a function of i. The highest number of patient arrivals is registered between 9am and 10am with about 5.9 patients on average. Between 5am and 6am the average number of patients is instead less than 1. As expected, the time of the day has a strong impact on the number of patient arrivals. The variance, V ar H,i , has very similar values to those of the mean E H,i , which is compatible with a well-known characteristic of the Poisson process, hence suggesting that the patient arrival process may be adequately modelled through a Poisson process with time inhomogeneous intensity. To further investigate this issue, we calculate the distribution of the number of arrivals during a given hour of the day and compare it the Poisson distribution with the same mean. The resulting probability mass function (pmf) is depicted in Fig. 2 for the interval 1pm-2pm and 5am-6am. The interval 1pm-2pm is the hour of the day where the Poisson distribution differs most from the actual distribution of the number of arrivals (the difference is calculated as the sum of absolute differences in the pmf) while the interval 5am-6am is the hour where the Poisson distribution is most similar to the experimental distribution. The hourly autocorrelation function (acf) is defined as

and it is depicted in Fig. 3 . The oscillation present in the autocorrelation is due to the changing patient arrival intensity during the day (Fig. 1) and it remains unaltered even for very large values of n.

We study now Y D,i in order to determine if there is a pattern in the patient arrivals during a week. In Fig. 4 we depict E D,i = E [{Y D,i+k |k = 0, 7, 14, ...}] and V ar D,i = V ar [{Y D,i+k |k = 0, 7, 14, ...}] with i = 1, 2, ..., 7 and observe that there are not large differences between days of the week. The variance differs to a large extent from the average on most days of the week. We calculated also the acf of Y D,i and saw that there is only a very mild correlation between number of arrivals of consecutive days. For this reason we do not aim to model daily correlation with our models. Finally, we look at the series containing the inter-arrival times themselves, X i with i = 1, 2, 3.... The average inter-arrival time is 19.43 and its variance is 483.18. The largest value in the series is 337 min. In Fig. 5 we provide the front and the tail of the probability density function (pdf) of the inter-arrival times, respectively. The autocorrelation present in X i is depicted in Fig. 6 where there is oscillation due to the fact that a period with high arrival intensity (8am-22pm) are followed by a period with low arrival intensity (22pm-8am). This oscillation however is less easy to interpret than that in Fig. 3 and fades away as n grows because for large values of n there is no deterministic connection between the time of the day of the ith and the (i + n)th arrival.

Application to the Cantù Dataset

The simplest arrival models are renewal processes, i.e., processes in which consecutive inter-arrival times are independent and share the same distribution. To have a Markovian renewal process, we can approximate the observed inter-arrival time distribution by phase type (PH) distributions. A degree-n PH distributions is given by the distribution of time to absorption in a Markov chain (MC) with n transient states (the phases) and one absorbing state [25] . A PH distribution is continuous if a continuous time MC (CTMC) is used and it is discrete if a discrete time MC (DTMC) is applied. In this work we apply continuous PH (CPH) distributions.

In order to understand basic concepts about CPH models, consider the degree-3 CPH distribution depicted in Fig. 7 where numbers inside a state indicate the initial probability of the state and numbers on the arcs are transition intensities. The gray state is the absorbing one. A CPH distribution is conveniently described by the vector of initial probabilities of the transient states denoted by α (we assume in this work that the initial probability of the absorbing state is 0) and the matrix, denoted by Q, containing the transition intensities among the transient states and the opposite of the sum of the intensities of the outgoing transitions in the diagonal (this allows to determine the transition intensities toward the absorbing state). For the CPH in Fig. 7 we have α = (0.2, 0.8, 0),

The pdf of a CPH distribution can be obtained based on basic properties of CTMCs as f (x) = αe Qx (−Q)1I where 1I is the column vector of 1s. From the previous the moments can be derived and have the form m i = i!α(−Q) −i 1I. We depicted the pdf of the above example in Fig. 7 . Several well-known distributions are contained in the family of CPH distributions. The exponential distribution is a degree-1 CPH distribution. Also hyper-exponential and Erlang distributions are simple to represent as CPH distributions.

Our aim is to find α and Q such that the pdf of the associated CPH distribution is a good approximation of the pdf of the patient inter-arrival times. There are two families of approaches to face this problem. The first is based on the maximum likelihood principle, i.e., having set the degree n to a given value, α and Q are searched for such that it is most likely to reproduce the empirical value at hand. Techniques based on this approach suffers from two drawbacks: first, the number of parameters in α and Q grows quadratically with n; second, the vector-matrix representation provided by α and Q is redundant and this gives rise to difficult optimization problems. One of the first methods in this line is described in [9] and an associated tool is presented in [17] . The second family of approaches is based on matching a set of statistical parameters of the available data. The set usually includes only moments [5, 18] but can refer also to characteristics of the pdf [4] or the cumulative distribution function at some points [15] . We apply here moment matching techniques which require much less computational effort, result in good fit of the inter-arrival distribution under study and have the advantage of considering a non-redundant representation (based on the moments every entry of α and Q is determined directly). It is known that a degree-n CPH distribution is determined by 2n−1 moments (see, e.g., [11] ). Not any valid 2n−1 moments can be realized however by a degreen CPH distribution. The approaches proposed in [5, 18] construct a degree-n CPH distribution matching 2n − 1 moments if the moments can be realized by a degree-n CPH distribution. The methods are implemented in the BuTools package (http://webspn.hit.bme.hu/ ∼ butools/) and we applied them with n = 1, 2 and 3. In the following we report the initial probabilities and the transition intensities of the resulting CPH distributions: The pdf of the above CPH approximations is depicted in Fig. 8 . With n = 1, the resulting exponential distribution, which captures only the first moment, is a bad approximation both of the front and of the tail. With n = 2 the tail is approximated well but the shape of the front is rather different. With n = 3 both the front and the tail are captured to a satisfactory extent. It is evident that, even if the inter-arrival time distribution is well approximated, a renewal process is significantly different from the data set under study because it cannot exhibit correlation and periodicity. Nevertheless, as we will show in Sect. 5, the renewal process with a proper PH distribution can give good approximation of some performance indices. In the following we propose models in which there is correlation in the inter-arrival time sequence.

Continuous time Markovian arrival processes (MAP), introduced in [26] , are a wide class of point processes in which the events (arrivals) are governed by a background CTMC and they can be seen as the generalization of the Poisson process allowing for non-exponential and correlated inter-arrival times. We provide here a brief introduction to MAPs.

The infinitesimal generator matrix of the underlying CTMC of an n-state MAP will be denoted by the n × n matrix D. Arrivals can be generated in two ways. First, in each state i of the CTMC a Poisson process with intensity λ i is active and can give rise to arrivals during a sojourn in the state. Second, when in the CTMC a transition from state i to j occurs an arrival is generated with probability p i,j . A convenient and compact notation to describe a MAP is obtained by collecting all intensities in two matrices, D (0) and D (1) , in such a way that D (0) contains the transition intensities that do not generate an event and D (1) those that give rise to one, including the intensities of the Poisson processes in the diagonal of D (1) . Moreover, the diagonal entries of D (0) are set in such a way that we have D (0) + D (1) = D. Accordingly, the entries of D (0) and D (1) are

In most of the literature, the vector of steady state probabilities of the background chain, denoted by γ = (γ 1 , ..., γ n ), is used as initial probability vector of the model 1 . In this paper we will do the same. As an example, consider the 2-state MAP described by

which implies that during a sojourn in state 1 no arrivals are generated while during a sojourn in state 2 a Poisson process with intensity equal to 5 is active. Moreover, a transition from state 1 to state 2 generates an arrival with probability 2/3 while transitions from state 2 to 1 are always associated with an arrival. The steady state probabilities of the background chain are γ = (0.4, 0.6).

The well-known Poisson process is a MAP with a single state (n = 1). The renewal process used in Sect. 3.1 can be expressed in terms of a MAP by setting

The Markov modulated Poisson process, in which a Poisson process is active in every state of a background CTMC, is a MAP whose D (1) matrix contains non-zero entries only in its diagonal.

As for CPH distributions, two families of parameter estimation methods have been developed in the literature: the first based on the maximum-likelihood principle ( [13, 27, 29] ) and the second based on matching a few statistical parameters of the arrival process. The representation given by D (0) and D (1) contains 2n 2 −n parameters (in every row of D = D (0) +D (1) the sum of the entries must be zero) and it is redundant as an n-state MAP is determined by n 2 + 2n − 1 parameters (2n − 1 moments of the inter-arrival times and n 2 joint moments of consecutive inter-arrival times; see [30] for details). Maximum-likelihood based methods suffers from the same drawbacks described before in case of CPH distributions. In this paper we experiment with methods that belong to the second family of approaches.

A 2-state MAP is determined by 3 moments of the inter-arrival times and the lag-1 auto-correlation of the inter-arrival time sequence [10, 30] . Our sequence has however such a lag-1 auto-correlation that cannot be realized with only 2 states. In [19] a method was proposed that creates a MAP with any 3 inter-arrival time moments and lag-1 auto-correlation. This method, implemented in the BuTools package, provides the following 6-state MAP: As opposed to the approach used in Sect. 3.1, in the arrival process generated by a MAP there can be correlation between subsequent inter-arrival times. In Fig. 9 we depict the acf in the inter-arrival time sequence generated by the above 6-state MAP and that computed on the available dataset. As guaranteed by the applied method for n = 1 the auto-correlation is matched exactly but then it fails to follow the auto-correlation of the data. The acf of the sequence counting the number of arrival per hour generated by the same MAP is given instead in Fig. 10 which is very different from the that of the dataset. In the number of patients per day sequence the auto-correlation of the 6-state MAP with n = 1 is 0.00469433, i.e., it is negligible, as opposed to that in the data where it is 0.22. As seen above, 6 states are necessary to match three moments and just the lag-1 autocorrelation of the inter-arrival times. This indicates that a MAP with large number of states is necessary to capture the peculiar statistical features of the patient arrival process. General purpose MAP fitting techniques are not applicable however with such large number of states.

A Hidden Markov model (HMM) [28] can be thought of as a generalisation of a DTMC for which an external observer cannot directly see the states but only observe some output whose probability to be emitted depends on the state. In practice an HMM is characterised by: i) a set of states S = {s 1 , . . . , s n }; ii) a set of possible observations O = {o 1 , . . . , o m }; iii) a n × n state-transition probability matrix A = {a ij } with a ij being the probability of transitioning from state s i to s j ; iv) a n × m state-observation probability matrix B = {b ij } where b ij is the probability of observing o j when the chain enters state s i ; v) an initial distribution over the set of states S denoted by α. Accordingly, an HMM is completely determined by the triple (α, A, B) .

There exist three classical problems for HMM, all of which requiring a sequence of observations O = (o 1 , . . . o k ) . The first one is the evaluation problem: given an HMM by (α, A, B) and a sequence O, calculate the probability that the HMM produces O. The second, called decoding problem, consists of determining the most probable state sequence given a HMM and a sequence O. The third one, the learning problem, has only O as input and is about finding such triple (α, A, B) with which observing O is most probable.

In this paper we focus on the third type of problem and in particular we develop and study two HMMs aimed at reproducing the statistical characteristics of the patients' arrival dataset (described in Sect. 2). We start off with a rather coarse-grained (only 3-state) but general HMM, which reveals to be capable of reproducing only marginally the auto-correlation characteristics of the patients arrivals data. Then we propose a HMM with particular underlying DTMC that shows good agreement with the dataset from a statistical point of view. interpreted as the number of patient arrivals per hour 2 . This means that the sequence generated by the HMM has to be post-processed if we need to specify the exact arrival instance for each patient. This post-processing will consist of distributing the arrivals in uniform manner inside the hour.

There are 2 free parameters in π, 3 × 2 = 6 in A and 3 × 14 = 42 in B (because π must be normalized and also the rows both in A and B must be normalized). We applied the Baum-Welch algorithm [8] to determine the optimal parameters starting from several initial parameter sets chosen randomly. With this relatively small model, the final optimal parameters obtained by the Baum-Welch algorithm are independent of the initial values (apart from permutations of the states). The obtained HMM is with π = (0, 1, 0) and 

The mean number of arrivals per hour with the above HMM is 3.0862 and its variance is 5.67669 while in the data trace the mean and the variance are 3.08676 and 5.69061, respectively. In Fig. 11 (left) we depict the autocorrelation of the number of arrivals per hour. As one can expect, the 3-state HMM is not able to reproduce the sustained oscillation present in the data shown in Fig. 3 (autocorrelation is negligible for n ≥ 10) but does much better than the 6-state MAP given in (3-4) (see Fig. 10 ).

A 24 States HMM. In order to have a model that is able to exhibit oscillation in the autocorrelation function of the number of arrivals per hour, we define a 24-state HMM in which each state corresponds to an hour of the day and the transition probabilities are such that the process deterministically cycles through the 24 states. Accordingly, A is 24 × 24 and its entry in position (i, j) is 1 in j = (i + 1) (mod 24) and it is 0 otherwise. The initial probability vector is set to π = (1, 0, ..., 0). As before, the possible observations are O = {0, 1, ..., 14}.

The Baum-Welch algorithm is such that parameters set to 0 initially remain 0. Consequently, the algorithm does not change the matrix A and the vector π and has an effect only on the entries of B. The number of parameters is larger, it is 24 × 14 = 336, but thanks to the deterministic behavior of the underlying DTMC the Baum-Welch algorithm determines the parameters in a single iteration. With the resulting HMM the mean number of arrivals per hour is 3.08619 while the variance is 5.67465. The hourly autocorrelation is shown in Fig. 11 (center). This model provides very similar hourly autocorrelation to that of the data set.

We experimented also with a 24-state HMM letting the Baum-Welch algorithm to change any entry of the matrix A (i.e., the initial entries of A are random strictly greater than 0 and strictly smaller than 1). This way the number of free parameters is 24 × 23 + 24 × 14 = 888. With this large number of parameters the Baum-Welch algorithm performs 4798 iterations and requires about 3 min of computation time on a standard portable computer. The resulting HMM has mean and variance equal to 3.08539 and 5.67468, respectively. The matrix A has a similar structure to the one before but the probabilities of going to the next state is not 1 but a value between 0.6 and 0.99. The resulting autocorrela- tion structure is depicted in 11 (right). With the non-deterministic underlying Markov chain the autocorrelation cannot exhibit sustained oscillation and the autocorrelation vanishes after about 150 hours (i.e., 6-7 days).

In order to assess the validity of the patients arrivals models presented in Sects. 3.1-3.3, we define a formal encoding of, respectively, the renewal Markov arrival model, the MAP and the HMM models, in terms of a superset of the Generalised Stochastic Petri Net (GSPN) [23] formalism, that we refer to as extended GSPN (eGSPN). Specifically, eGSPN is a class of stochastic Petri nets that, like GSPN, allows for combining immediate and stochastic timed transitions, but that, differently from the original GSPN, is not constrained to exponentially distributed timed transitions. We will use the eGSPN models described in this section to run a validation through the assessment of a number of key performance indicators (KPIs, see Sect. 5). For the sake of space we only give the syntactic elements necessary for characterizing an eGSPN, while we omit the formal semantics. We remark that if it is well known that the semantics of a GSPN can be reduced (by elimination of vanishing markings) to a CTMC process, equivalently it can be shown that the semantics of an eGSPN model corresponds to a generalised semi-Markov process (GSMP) [16] , i.e. a larger class of stochastic processes that subsumes CTMCs. For what concerns pol, single server and infinite server semantics are as usual [23] .

In the following we define a Petri net that represents the patient arrival process in which inter-arrival times follow a degree-n PH distribution given by the vectormatrix pair, (α, Q) (as described in Sect. 3.1). The set of places is P = {P 0 , P 1 , ..., P n , P n+1 , Arrs} where P 0 will be used to start the process according to α after each arrivals, P 1 , ..., P n correspond to the transient states of the PH distribution, P n+1 corresponds to the absorbing state and Arrs is the place where patient arrivals are accumulated. The set of immediate transitions is T I = {t arr ∪ {t 0,i : 1 ≤ i ≤ n, α i > 0}} where the weight of transition t arr is 1 and the weight associated with transitions t i , 1 ≤ i ≤ n, α i > 0 is α i . Firing of t arr will generate an arrival while firing of a transition t i will imply that phase i is chosen as initial state in the MC associated with the PH distribution. Exponential transitions are used to model the transitions among the phases. Their set is

where the first set in the union contains transitions among the transient states while the second those to the absorbing state. Accordingly, the rate associated with these transitions is Q i,j for transition t i,j if 1 ≤ i, j ≤ n and it is − n j=1 Q i,j for transition t i,n+1 with 1 ≤ i ≤ n. The set of transitions is T = T I ∪ T E . Transition t i,j , 0 ≤ i, j ≤ n + 1, if present, is with input bag {P i } and output bag {P j }. Input bag of transition t arr is {P n+1 } and its output bag is {P 0 , Arrs}. In Fig. 12 we depict the PN representation of the arrival process in which inter-arrival times are PH according to the example proposed in Sect. 3.1. For immediate transitions (black bars) we specify the weight while for exponential ones (white rectangles) the rate.

Now we turn our attention to MAPs defined by the pair of matrices D (0) and D (1) as described in Sect. 3.2.

In case of an n-state MAP, the set of places is P = {P 0 , P 1 , ..., P n , Arrs} where P 0 is used to start to model, places P 0 , ..., P n correspond to the states of the MAP and place Arrs is where the patients are gathered. A set of immediate transitions is used to start the process according to the steady state probabilities, γ. This set is T I = {t 0,i : 1 ≤ i ≤ n, γ i > 0} in which each transition t 0,i is associated with weight γ i . In the set of exponential transitions, a transition is associated with each positive entry of the matrices D (0) and D (1) . Accordingly, the set is

where the second set in the union is the set of transitions that generate arrival. The rates of transitions t i,j and t * i,j are D (0) i,j and D (1) i,j , respectively. The set Figure 13 shows the eGSPN encoding of the arrival process using the MAP specified in (2). 

In the following we provide the PN encoding of an n-state HMM given by the triple (π, A, B) , assuming that the observations are interpreted as number of arrivals in an hour and distributing the arrivals within the hour in a uniform manner. The set of places is P = {P 0 , P 1 , ..., P n , O 1 , ..., O n , W, Arrs} where P 0 is used to start the model according to the initial probabilities given in π, places P 1 , ..., P n correspond to the states of the background chain, places O i , 1 ≤ i ≤ n are used to emit the observation when the background chain enters state i, place W collect arrivals that has to be distributed in a given hour, and place Arrs gathers all arrivals. A set of immediate transitions, T I,1 = {t 0,i : 1 ≤ i ≤ n, π i > 0}, is used to start the process according to π. A set of deterministic transitions models the background Markov chain:

in which all transitions are associated with fixed delay of one time unit and t i,j , if present, is with weight A i,j . Another set of immediate transitions,

i,j is with weight B i,j , is used to emit observations. In order to distribute the arrivals inside an hour a uniform transition is used t dist whose minimal firing time is 0 and maximal firing time In Fig. 14 we show a part of the PN that represents the 3-state HMM given in (5) . Deterministic transitions are drawn as gray rectangles and are provided with their weight. The single uniform transition is drawn as a black rectangle and labeled with its firing interval. Number of occurrences of a place in a bag is written as a label to the corresponding arc if it is different from 1. In the figure we omitted transitions t (o) i,• with i = 1, 2 and connecting arcs.

In order to validate the different models of patients arrival process we have coupled each of them with an elementary model of the ED services a patient may go through during his permanence at the ED and assessed meaningful KPI on the coupled models. For the sake of simplicity we assumed the ED patient flow consisting of a simple pipeline of 3 services: triage, visit and discharge for each of which we considered exponentially distributed service time with the following settings 3 : triage∼Exp (12) , the visit∼Exp (6) and discharge∼Exp(60) (i.e. visit is assumed to be the slowest while discharge the fastest service). Figure 15 shows the eGSPN encoding of the ED service model: notice that place Arrs is shared with the eGSPN models of the patient arrivals (representing the composition of the two models). To validate each arrival model we compared their performance with that resulting by the simulation of the ED service model (Fig. 15) with the actual arrivals of the dataset (for this we generated an eGSPN encoding of the dataset arrivals). One relevant KPI that we considered for comparing the performances of the different patient arrival models coupled with the ED services model is:

φ 1 ≡ pmf of the number of patients in the ED during a periodic time − window.

We formally specified φ 1 by means of the Hybrid Automata Specification Language (HASL) [7] and assessed it through the statistical model checking platform Cosmos [1, 6] . HASL model checking is a procedure that takes a stochastic model (in terms of an eGSPN) as input as well as a linear hybrid automaton (LHA) together with a list of quantities Z i to be estimated. The procedure uses the LHA as a monitor (i.e., a filter) to select trajectories which are automatically sampled from the eGSPN. The statistics, which are stored in the LHA variables, are collected on the trajectories and used to obtain a confidence interval for each quantity Z i associated with the LHA monitor (we refer the reader to [7] for more details). Figure 16 depicts the HASL specification for φ 1 . The LHA consists of 4 states (l 0 ,l on ,l off ,l end ), 1 clock variable (t), 2 stopwatches (t a , t p ), M + 1 real valued variables x i (0 ≤ i ≤ M , whose final value is the probability that i patients have been observed in the ED during the periodic time windows), plus a number of auxiliary variables and parameters. The LHA is designed so that measuring, along a trajectory, is periodically switched ON/OFF (parameter TP being the ON period duration) and stops as soon as the NP -th period has occurred. In the initial state l 0 the occur- and in so doing the timer t a (which stores the duration of the occupation at i patients since the last arrival/departure) is reset so that it can correctly be used in the freshly started ON-period. Finally, the LHA stops monitoring as soon as NP ON-periods have been observed along the monitored trajectory: at that moment each x i is normalised w.r.t. the sum of all x i , hence on ending the monitoring of trajectory x i is assigned with the probability that the i patients have been observed in the ED during NP observed periods

. Observe that such an LHA can also be used for "non-periodic" measures: it suffices to set TP = 0. The HASL specification for φ 1 is completed by the list of HASL expressions AVG(last(x i )) which indicate that a confidence interval for average of the last value that x i has at the end of an accepted trajectory is computed by COSMOS. Figure 17 depicts plots computed with the COSMOS tool and resulting from assessing specification φ 1 (i.e., the pmf of the number of patients in the ED) against the eGSPN models of the patient arrivals coupled with the 3-services model of the ED (Fig. 15 ). Plotted results refer to a 1-year (365 days) observation window and have been computed as 99.99% confidence interval of 10 −3 width. The plot on the left refers to the pmf measured without taking into account any specific time-window over the day (continuous measuring over 24 h for 365 days), while the plots in the center and on the right refer to periodic measuring over a 1-hour period, at low arrival intensity (from 2am to 3am, center picture) Fig. 17 . Pmf of the num. of patients in the ED computed through 3-phases PH renewal model, 6-state MAP, 3-state HMM and 24-state HMM versus dataset model over the whole day (left), during a low-arrival hour (center) and high-arrivals hour (right). and at high arrival intensity (from 11am to 12pm, right picture), respectively. The obtained results witness the clear advantage of the HMM24 model over the rest: if when no specific hour of the day is considered the pmf of the 3-phases renewal model and of the 6-states MAP provide an acceptable approximation of the pmf computed w.r.t. to dataset (red dashed plot) and the pmf of the HMM3 and HMM24 are essentially indistinguishable, this is no longer the case when the pmf is computed on a specific time window as shown in the center and right plots. Only the pmf of the HMM24 matches that of the dataset during a low-arrivals [2am,3am] window and a high-arrival [11am,12pm] window: the pmf of the renewal and map models instead are essentially identical to those measured continuously while that for the HMM3 exhibits a slight tendency towards matching the pmf of the dataset at low-arrivals but completely fails to do so when arrivals become intense.

In this paper we have considered the problem of deriving stochastic models that are capable to accurately reproduce the statistical nature of patients arrivals process observed on a real dataset of the ED of an Italian hospital. Starting from the statistical analysis of the dataset we figured out relevant characteristics of the patients arrival process, e.g., the hour of the day with highest/lowest number of arrivals, the moments, the variance and, most importantly, the periodic nature of number of arrivals-per-hour, the latter shown by auto-correlation function estimated on the dataset. In quest for a model capable of capturing all of these aspects we have considered different classes of Markovian models namely, renewal Markov processes, Markov arrival processes and hidden Markov models. Based on the considered dataset, we derived for each class a few model instances (1, 2 and 3-phase PH renewal process, a 6-state MAP, a 3-state and a 24-state HMM). We assessed the quality of each model by comparing their basic statistical characteristics (i.e., moments, auto-correlation) with those assessed initially on the dataset. To complete the validation process we then considered the coupling of patient's arrival model with a simple Petri net model of the ED services. To this aim we gave a formal encoding of the different kind of patients arrival Markov models in form of stochastic Petri nets. We then evaluated each model by estimation of the pmf of the number of patients in the ED during a given period of the day which we formally encoded and accurately assessed through HASL statistical model checking. The obtained results showed that the 24-states HMM is the only model, amongst those considered, capable of reproducing the statistical characteristics of the dataset, including the hourly periodicity of the arrivals. Therefore such a model could be validly used for predicting the performances of a realistic ED design once coupled with a more realistic model of the internal ED services.

Short-term forecasting of emergency inpatient flow

Forecasting the emergency department patients flow

Constructing matrix exponential distributions by moments and behavior around zero

On moments based Padé approximations of ruin probabilities

COSMOS: a statistical model checker for the hybrid automata stochastic logic

HASL: a new approach for performance evaluation and model checking from concepts to experimentation

A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains

A benchmark for PH estimation algorithms: results for Acyclic-PH

A Markovian canonical form of second-order matrix-exponential processes

Moment characterization of matrix exponential and Markovian arrival processes

Regression forecasting of patient admission data

An EM-algorithm for MAP fitting from real traffic data

Fitting mixtures of exponentials to long-tail distributions to analyze network performance models

Stochastic Petri Nets -Modelling, Stability, Simulation

PhFit: a general phase-type fitting tool

Matching more than three moments with acyclic phase type distributions

Matching marginal moments and lag autocorrelations with maps

Applying queueing theory to the study of emergency department operations: a survey and a discussion of comparable simulation studies

A multivariate time series approach to modeling and forecasting demand in the emergency department

Estimating the waiting time of multi-priority emergency patients with downstream blocking

Modelling with Generalized Stochastic Petri Nets, 1st edn

Review of predicting number of patients in the queue in the hospital using Monte Carlo simulation

Probability distributions of phase type

A versatile Markovian point process

Faster maximum likelihood estimation algorithms for Markovian arrival processes

A tutorial on Hidden Markov models and selected applications in speech recognition

An EM algorithm for estimation in Markov-modulated poisson processes

A minimal representation of Markov arrival processes and a moments matching method

A data-driven model of an emergency department