Determination of server location in emergency care systems: an index proposal using data envelopment analysis and the hypercube queuing model

ECS location problem and Hypercube queuing model

A common goal in managing ECS location problems is optimizing the coverage of a given location, such as a city. An area is considered covered when it is located less than a critical distance from at least one of the existing servers. Location problems based on cluster coverage are prescriptive and their objective is to find the optimal location for server units, in order to guarantee total coverage and the maximization of performance measures (Brotcorne, Laporte & Semet, 2003).

An important advance in the theory related to ECS location problems was the proposal of probabilistic coverage models, which consider randomness in server availability as a significant factor. The advantage of these models is that they consider emergency call uncertainties, such as when a call enters the system and its location and duration (Swersey, 1994).

The first probabilistic ECS location problem solutions focused on maximizing system coverage, considering the probability of a given server unit being busy. Daskin (1983) presented a pioneering study in which location models were extended considering this possibility. ReVelle & Hogan (1989) presented a model whose objective was to maximize population coverage within a response time established by a level of reliability. Recent research has sought to relax the three simplifying assumptions presented in previous studies, namely: server independence, identical server occupancy rates, and server occupancy rate invariant to location.

As a way of relaxing these assumptions, some studies have used the hypercube queuing model to calculate the probabilities of server occupation. The hypercube utilizes spatially distributed queues in probabilistic location scenarios. The queuing system performance indicators offered by the hypercube model can be used in prescriptive models to assess the optimal location of a server (Brotcorne, Laporte & Semet, 2003).

From the metrics produced by the hypercube model, it is feasible to aggregate inputs and outputs into a composite index. This process can be accomplished through the implementation of Data Envelopment Analysis (DEA).

Data envelopment analysis (DEA) and benefit of the doubt (BoD) approach

DEA is a technique developed by Charnes, Cooper & Rhodes (1978) to determine the efficiency of each unit under evaluation (the so-called decision-making unit - DMU) by aggregating the multiple inputs (raw materials, energy, capital, etc.) and outputs (products, income, customers served, etc.) of a system. Since its inception, the field of application for DEA has expanded, with new models and applications of the technique appearing every year (Liu et al., 2013).

For example, recent efforts have been made to evaluate the efficiency of transport systems using DEA in combination with other techniques. Taletović & Sremac (2023) utilized Principal Components Analysis (PCA), while Blagojević et al. (2020) employed fuzzy Analytic Hierarchy Process (AHP). However, these methodologies are primarily employed to streamline the quantity of input and output factors or to determine the optimal prioritization of metrics.

One of the most relevant DEA applications is in the construction of composite indexes (CIs) to measure phenomena considering multiple dimensions (Cherchye et al., 2006; Mariano, Ferraz & Gobbo, 2021). A composite index (CI) is an indicator that aggregates several independent variables into a single measure. According to Booysen (2002), the construction of a CI involves five steps, which are: (1) indicator selection (redundancy must be avoided and data availability ensured); (2) standardization (a set of procedures to group indicators under the same unit of measurement); (3) weighting (how to determine the degree of importance of each indicator); (4) aggregation (methods for combining indicators into a single measure); and (5) validation, both internal (index consistency) and external (comparison with other indicators). DEA provides solutions for the normalization, weighting, and aggregation steps of the indicators.

The application of DEA for the construction of composite indices is called benefit of the Doubt (BoD), as systematized by Cherchye et al. (2006). There are several different models of BoD and Mariano, Ferraz & Gobbo (2021) have systematized and compared most of these models, using them to recalculate the Human Development Index (HDI). Mariano, Sobreiro & Rebelatto (2015) states that DEA/BoD can be used to build CIs in two ways: (a) in CIs composed only of desirable attributes (in this case, a constant input equal to 1 is adopted for all DMUs); and (b) in CIs expressed as a ratio of desirable and undesirable indicators, being the approach used in this work.

The measures within this approach are treated as outputs when they contribute positively to ECS performance and as inputs when they are undesirable for ECS performance.

To measure the performance of an ECS, two metrics are considered:

1. Response time, considering that a quick response can save many lives; and

2. Workload, which may compromise service quality due to employee stress or burnout or equipment overload etc.

Response time is considered a desirable indicator because it is a core characteristic of an ECS; workload, on other hand, is considered an undesirable indicator because it is an unavoidable negative consequence of an ECS.

Proposed composite index

The index proposed in this study integrates the hypercube model and DEA/BoD to evaluate the location of servers in emergency care systems. Four vital steps are required to build the index: (a) enumerating all possible configurations for the ECS studied (DMU determination); (b) applying the hypercube queuing model to determine performance indicators (inputs and outputs) for each configuration; (c) using DEA/BoD to build a composite index in order to rank different ECS configurations (CI construction); and (d) evaluating the CI of each configuration and the ranking thus obtained to determine the best solution for the real problem (validation). Figure 1 summarizes the four stages of CI construction.

At the conclusion of the four steps, the configurations can be ranked according to their performance. From this classification, an analysis of whether there are clear patterns between the best and worst configurations presented can be made, which will guide the future actions of the emergency system decision makers. It is important to note that the implementation of any solution suggested by the model must be analyzed in terms of its technical and economic feasibility. For example, in the case of a reallocation of servers, the potential costs of adapting infrastructure and expanding staff should be considered.

Enumeration of configurations (DMU determination)

All location patterns must be compared to define the best configurations. These configurations are defined from a list of all the distribution possibilities of the ‘n’ servers in the ‘A’ atoms, which are spatial subdivisions of the region studied. Using the concept of mathematical partition, all permutations for the system studied can be defined.

Firstly, all partitions of the number of servers ‘n’ must be obtained. For example, if there are n = 4 servers available, five partitions can be obtained: “4”, “3 + 1”, “2 + 2”, “2 + 1 + 1” and “1 + 1 + 1 + 1”. These partitions have a number of terms equal to 1, 2, 3 and 4, respectively. It should be noted that the number of servers used in this calculation includes only those that can be allocated to different atoms. Thus, considering the previous example, if the number of atoms were A = 3, the partition “1 + 1 + 1 + 1” could not be applied to the problem in question.

The remaining ‘P’ feasible partitions must be filled with zeros, so all partitions have a number ‘A’ of terms. So that, in the provided example, the following four adjusted partitions can be observed: P1 = “4 + 0 + 0”, P2 = “3 + 1 + 0”, P3 = “2 + 2 + 0” and P4 = “2 + 1 + 1”.

The total number of configurations can be obtained by permuting the terms of the filled partitions. For instance, in the case of partition P1, three potential configurations are present: “4 + 0 + 0”, “0 + 4 + 0” and “0 + 0 + 4”. Each permutation, that is, each distribution of ambulances to an atom, will be a different DMU, to be evaluated by DEA in later stages.

The number of configurations (permutations) for each adjusted partition ‘ P_k’ can be obtained according to Eq. (1). The total number of configurations of a problem is the sum of the configurations generated by all partitions P_k, as shown in Eq. (2). (1)${\displaystyle C_k=\frac{A!}{\prod _{\forall j}{R}_j!.}}$

where:

C_k: number of configurations of the partition “P_k”;

C: total of configurations;

P: quantity of feasible partitions;

A: number of atoms;

R_j: number of repetitions of the term j in the partition.

For the case with n = 4 basic servers and A = 3 atoms, the adjusted partition P₄ = “2 + 1 + 1”, generates a number of permutations equal to $\frac{3!}{2!}=3$; in this case, R₁ = 2 because the term “1” is repeated twice. The adjusted partition P₂ = “3 + 1 + 0”, on the other hand, generates a number of permutations equal to 3! = 6 because it does not include repetitions in the terms of the partition. All the 15 configurations generated in this hypothetical case are expressed in Table 1.

If the system consists of distinct types of servers, e.g., Advanced Server Units (ASU) and Basic Server Units (BSU), which differ with respect to the equipment and staff sent to the call, the enumeration must be performed separately for each type of server and, at the end, the results must be multiplied. If, for example, in addition to the four BSUs there was one ASU, three configurations would be generated (see Table 1). The total number of configurations to be evaluated, therefore, would be 45, equivalent to 15 multiplied by three.

As the number of configurations grows exponentially with the number of servers and atoms, the problem may become computationally unviable. Table 2 shows the configurations resulting from different numbers of atoms and servers. For example, in the situation with 10 atoms and 10 servers the number of configurations is 92,423. This computational constraint might not make applications in real problems unfeasible, since, even for large cities, the problem could be divided into service regions, with fewer atoms. This subdivision of large cities is to be expected, as an ambulance would hardly answer a call in a place very far from its atom.

Based on the definition of possible permutations, performance measures are then calculated. Each permutation represents one possible configuration of the distribution of servers in an atom. The hypercube model will be applied for each configuration, determining its performance measures. A description of how the hypercube model should be built can be seen in the subsequent step.

Inputs and outputs determination

The second stage involves the application of the hypercube model for each of the configurations built in the previous stage to determine the response times (outputs) and the workloads (inputs).

The original hypercube model is a descriptive model used as a tool for the analysis and planning of urban emergency systems (Larson, 1974). In addition to considering uncertainties regarding the origin of calls, service times and the availability of servers, the model addresses the geographic and temporal complexities of a region, based on spatially distributed queues. Basically, the idea is to expand the state space of an M/M/m queue system to represent each server individually, which may include more complicated dispatch policies.

The solution of the model is initiated by the construction of a set of equilibrium equations for the system. These equations describe, in terms of probabilities, the transition between states that have reached a steady state (Larson, 1974; Larson & Odoni, 2007); that is, the probabilities are independent of time (Shortle et al., 2018). For the hypercube model, a state describes the status of given servers. For example, the state {110} of a system with three servers represents the situation where two servers are busy, and one is free. An equilibrium equation is established for each state transition (for example, {001} to {011}). It is important to mention that the number of states of the hypercube model grow exponentially according to the number of servers, in the order of 2n.

The results of the equilibrium equations are based on the values of the steady-state probabilities of the system, allowing for the calculation of performance measures, such as: server workload, average response time of the system or of each server, and dispatch frequencies, among others (Larson & Odoni, 2007).

To apply the classic hypercube model, nine critical assumptions must be verified, as described by Larson & Odoni (2007). Extensions for relaxing these hypotheses, providing a more accurate representation of reality, are found in the literature. Most of these extensions are related to the dispatch policy and operational characteristics of care systems, as can be seen in the work of Takeda, Widmer & Morabito (2007), Souza et al. (2015) and Chiyoshi, Iannoni & Morabito (2011).

In our proposal, focused on problems of the localization of emergency systems, the extension presented by Beojone, Máximo de Souza & Iannoni (2020) was used (see Appendix). The proposed extension considers the aggregation of servers to reduce the number of equilibrium equations necessary to represent the queuing system; it is especially important in large systems with multiple servers of the same type in one place and allows increasing computational efficiency. The use of servers fully dedicated to answering calls with high priority (where there is a risk of death for the patient) is also another important aspect to be considered in this model.

To use the server aggregation approach, Eq. (3) must be applied to calculate the number of states |M| for a system with T groups of servers. These groups are the partition terms that were defined in the previous section. Servers belonging to the same group must be homogeneous and indistinguishable; that is, they serve the same type of occurrence and are positioned within the same atom. (3)${\displaystyle \left|M\right|=\prod _{t=1}^T\left(n_t+1\right)+\left|Q\right|}$

Each group t has n_t servers; |Q| represents the number of queue states; that is, all possible combinations of server occupancy and queue size. It is noted that the number of equilibrium equations in the server aggregation model is significantly less, especially in large systems, than the number of equilibrium equations in the exact hypercube model, as presented by Beojone, Máximo de Souza & Iannoni (2020).

As in this model, if it is considered that there are servers dedicated to calls with high priority in emergency systems, queues can form, even with servers available (Iannoni, Chiyoshi & Morabito, 2015; Rodrigues et al., 2017). This change in the configuration of the exact hypercube model requires important adjustments to the equilibrium equations and to the queue. Thus, Eq. (4) calculates the number of queue states: (4)${\displaystyle \left|Q\right|=\left|Q_{ds}\right|+\left|Q_{pr}\right|=n_{ds}{Q}_l+\left(\begin{array}{c}\hfill 2+Q_l\hfill \\ \hfill Q_l\hfill \end{array}\right)-1}$

where:

• |Q_ds| is the number of queue states that arise from the model due to fully dedicated servers;

• |Q_pr| is the number of queue states in saturated circumstances; that is, when all servers are busy and there are priority levels in the queue;

• n_ds is the size of the group of fully dedicated servers;

• Q_l is the queue capacity.

One way to differentiate users and represent priorities in the hypercube model is to subdivide atoms into sub atoms (“layering”). This approach was used by Takeda (2000) and Iannoni (2005). Following the strategy presented in Larson & Odoni (2007), sub atoms are created by combining each geographic atom with different types of calls in terms of priority. The integrated model presented here adopts two levels of priority for regular and serious calls. Thus, the number of sub atoms will be double the number of atoms. The Appendix presents the equations referring to the generalization of the state-to-state transitions of the hypercube model (Beojone, Máximo de Souza & Iannoni, 2020) proposed in this integration.

Various internal and external performance measures can be calculated from the solution obtained, usually by some numerical method for the solution of linear equation systems (e.g., Gauss-Siedel). These measures include workload and response times, for each atom, server and region, and dispatch frequencies. The workload is calculated according to Eq. (5) (Souza et al., 2015). (5)${\displaystyle {\rho }_{it}=\sum _{r\in M}\frac{n_{tr}\cdot P_r}{n_t}}$

where:

ρ is the workload calculated for server i in group t;

r is a particular state in the set of states M;

n_tr is the number of servers in group k occupied in state r;

P_r is the probability of the state r.

Workload and average response time are the performance measures used in the DEA model proposed. The average response times for each atom and server were obtained by calibrating the original response time for each configuration. The calibrating procedure can be viewed in Chiyoshi, Galvão & Morabito (2001) and Galvão & Morabito (2008).

Composite index determination

In order not to make the review too extensive, this section will focus on describing the specific elements of the model adopted in this research. For basic definitions of DEA, it is recommended an introductory bibliography on the subject (Cooper, Seiford & Tone, 2007; Cooper, Seiford & Zhu, 2011).

The inputs for the model will be the workloads of each server. It is important for servers to be able to perform their functions effectively, without excessive workload. The inputs will be represented as ρ_ik, where the index k will represent each configuration (DMU) and i each server in this configuration; The number of inputs “I” is equal to the number of servers “n”.

The outputs of the model, on the other hand, are the reciprocal response times for each server to each sub atom, as one of the main objectives of an emergency system is to promote rapid response (Halkos & Petrou, 2019). The use of the reciprocal, as suggested by Golany & Roll (1989), was possible because the variables used did not present null values and because there was no need to maintain the proportionality of the scales for the purpose of this study (You & Yan, 2011). As it is an ECS location problem that deals with lives, the adoption of the reciprocal of the average response time of the servers can be justified, because it penalizes more longer response times.

In the present approach, for each configuration “k”, three groups of outputs related to the inverse of the average response times were adopted: (a) the inverse of the average response time for an atom h related to regular calls, which is represented by y_hk; (b) the inverse of the average response time for an atom h related to serious calls, which is represented by $y_{hk}^S$; and (c) the inverse of the average response time of a server I, which is represented by z_ik. The model considers A atoms and n servers. Therefore, the number of outputs “O” in the model is equal to two times the number of atoms (as they are divided into two sub atoms - regular and serious calls) plus the number of servers: O = 2A + n.

There are several models related to Data Envelopment Analysis (DEA), which, according to Mariano, Sobreiro & Rebelatto (2015), vary according to the establishment of different assumptions, such as the type of scale return, the type of efficiency measure (radial or non-radial), and the type of orientation (for input, output or both), among others. As the configuration for this integration proposal is to determine the efficiency of emergency services systems, the orientation to the output, that is, to increase the reciprocal of the response time to the user, is preferable. This is due to the social character of this type of system, in which the priority is to provide a better service level even at the expense of a higher input utilization rate. The primary objective of a health system is to guarantee quality of life through more treatment and assistance. Therefore, from a moral point of view, for this type of system the orientation to the output seems more adequate (Stefko, Gavurova & Kocisova, 2018).

In addition, the BCC model (Banker, Charnes & Cooper, 1984) seems to be the most appropriate choice, since it is expected that it will become increasingly costly, in terms of additional workload, to decrease the response time by keeping the other characteristics of the system constant. This premise is supported by the healthcare queuing literature (Palvannan & Teow, 2012).

To improve the model, the addition of weight constraints, which establish limits on the weighting of inputs and outputs, was necessary (Cherchye et al., 2006). The addition of weight constraints introduces subjectivity to the model and should be justified. In the case of emergency systems, it is usual to prioritize calls with respect to their severity and urgency. It makes sense, then, to include in the model the constraint that, for each atom, the DEA must give greater weight to the response times for the most urgent calls. This type of constraint is called the assurance region and was initially proposed by Thompson et al. (1986). This constraint is expressed in Eq. (12).

To prevent the model from assigning excessive weights to the outputs of a certain group and ignoring another, additional constraints can be inserted by establishing a minimum contribution limit for the efficiency of a given set of weights. This additional constraint considers that the various outputs can be divided into two distinct groups: (a) those referring to the inverse of the response time to atoms (both normal and serious calls) and (b) that of the inverse of the average response time of servers. If the DEA model assigned 100% weight to one group or another, efficient scenarios, in which response times for servers or atoms would be critically high, could be suggested as efficient configurations (Pedraja-Chaparro, Salinas-Jimenez & Smith, 1997).

Thus, in order to avoid the model attributing an excessive weight to the outputs of a certain group and ignoring another, it was established that the set of variables linked to the inverse of the response times to atoms should contribute with at least P% of the efficiency, while the variables related to the inverse of the response times of each server should also contribute with at least P% of efficiency, leaving the model free to allocate the missing (100-2P)%. This approach of attributing constraint to the relative contribution of groups of variables was used by Morais & Camanho (2011), inspired by the work of Wong & Beasley (1990).

An advantage of including constraints is to reduce the number of efficient tied DMUs, improving model discrimination (Wang, Luo & Liang, 2009). These constraints are expressed in Eqs. (9), (10) and (11). If a user chooses to reduce the subjectivity of the model, it is possible to omit restrictions Eqs. (9), (10) and (11), considering that this may result in the omission of some variables and an increase in the number of ties. To make the determination of P more objective, multicriteria decision-making methods can be applied (Kuo, Wang & Tien, 2010).

For each configuration “k” generated in the first step (C configurations in total), the linear programming problem related to the output-oriented BCC model with weight constraints must be solved. Equations (6)–(17) present a possible model in the multiplier form, taking two levels of urgency (normal and serious) and limiting the contribution to the efficiency of each group of outputs by at least P%. (6)${\displaystyle \mathrm{Min}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\gamma }_{\mathrm{i}}{\rho }_{\mathrm{i}0}-\mathrm{w}}$

Subject to

where:

• γ_i is the weight for the workload of server i;

• α_h is the weight for the inverse of the response time for regular calls in an atom h;

• ${\alpha }_h^S$ is the weight for the inverse of the response time for serious calls in an atom h;

• β_i is the weight for the inverse of the response time for a server i;

• ρ_i0 isthe workload of server i for the DMU under analysis;

• ρ_ik isthe workload of server i for a DMU k;

• y_h0 is the inverse of the response time for regular calls in an h atom, considering the DMU under analysis;

• $y_{h0}^S$ is the inverse of the response time for serious calls in an h atom, considering the DMU under analysis;

• y_hk is the inverse of the response time for regular calls in an atom h, considering a DMU k;

• $y_{hk}^S$ is the inverse of the response time for serious calls in an atom h, considering a DMU k;

• z_i0 is the inverse of the response time for a server i for the DMU under analysis;

• z_ik is the inverse of the response time for a server i, considering a DMU k;

• w is the scale factor;

• P is the minimum contribution to efficiency from each group of response times;

• ɛ is a non-Archimedean constant.

The values of the objective function are the efficiencies of each configuration. Since many configurations can result in the same efficiency values, it is necessary to conduct a last discrimination step to obtain a final efficiency ranking.

To deal with cases with many ties between configurations, there are several objective methods to differentiate them, such as the inverted frontier, cross-evaluation, or the triple index (Mariano, Sobreiro & Rebelatto, 2015), depending exclusively on the configuration studied. Of all these methods, the inverted frontier stands out for its simplicity, since it consists of just conducting the analysis using outputs as inputs and vice versa. The first to identify the properties of the inverted frontier was Yamada, Matui & Sugiyama (1994). Years later, Leta et al. (2005), Zhou, Ang & Poh (2007) and Mariano & Rebelatto (2014) proposed three different forms to use the inverted frontier as a discrimination method (Mariano, Ferraz & Gobbo, 2021).

The Leta et al. (2005) approach is the simplest of these methods; it is based on the construction of a composite index that is the normalized mean between the result of the standard frontier and one minus the result of the inverted frontier. It is important to mention that the inverted frontier proposed by Leta et al. (2005) was adopted in this research, but other approaches can also be applied in this step.