Evaluating relief center designs for disaster relief distribution

Merve Ozen (Department of Industrial and Systems Engineering, Madison Graduate School, University of Wisconsin, Madison, Wisconsin, USA)
Ananth Krishnamurthy (Department of Industrial and Systems Engineering, Madison Graduate School, University of Wisconsin, Madison, Wisconsin, USA)

Journal of Humanitarian Logistics and Supply Chain Management

ISSN: 2042-6747

Publication date: 3 April 2018



Relief item distribution to victims is a key activity during disaster response. Currently many humanitarian organizations follow simple guidelines based on experience to assess need and distribute relief supplies. However, the interviews with practitioners suggest a problem in efficiency in relief distribution efforts. The purpose of this paper is to develop a model and solution methodology that can estimate relief center (RC) performance, measured by waiting time for victims and throughput, for any RC design and analyze the impact of key design decisions on these performance measures.


Interviews with practitioners and current practice guidelines are used to understand relief distribution and a queuing network model is used to represent the relief distribution. Finally, the model is applied to data from the 2015 Nepal earthquake.


The findings identify that dissipating congestion created by crowds, varying item assignment decisions to points of distribution, limiting the physical RC capacity to control congestion and using triage queue to balance distribution times, are effective strategies that can improve RC performance.

Research limitations/implications

This research bases the RC designs on Federal Emergency Management Agency guidelines and assumes a certain area and volunteer availability.


This paper contributes to humanitarian logistics by discussing useful insights that can impact how relief agencies set up and operate RCs. It also contributes to the queuing literature by deriving analytic solutions for the steady state probabilities of finite capacity, state dependent queues with blocking.



Ozen, M. and Krishnamurthy, A. (2018), "Evaluating relief center designs for disaster relief distribution", Journal of Humanitarian Logistics and Supply Chain Management, Vol. 8 No. 1, pp. 22-48. https://doi.org/10.1108/JHLSCM-03-2017-0012

Download as .RIS



Emerald Publishing Limited

Copyright © 2018, Emerald Publishing Limited

1. Introduction

Each year disasters affect millions of people worldwide, and effective disaster response is essential to minimize the suffering of communities. Disaster response activities include, search and rescue operations, first aid medical attention, distribution of critical relief supplies and restoration of critical infrastructure. This paper focuses on the distribution of critical relief supplies to victims. Prior to being distributed to disaster victims, relief supplies go through the relief distribution supply chain. The relief distribution supply chain has three tiers: procurement, transportation, and distribution. This paper focuses on the distribution of relief items to victims at the disaster site.

Relief items are distributed to victims via temporary structures called relief centers (RCs), where victims arrive and wait in line to receive the supplies they need. The RCs are usually setup in large, open areas to allow for crowd access, while the particular setup of an RC can vary. Crowding experienced during distribution is a significant problem that can hinder the efficiency of relief distribution. In the case of Salvation Army, 2010 Haiti earthquake response, relief distribution efforts in Port-au-Prince came to halt due to the crowds converging around the RC. One of the Salvation Army responders remembers it as: “Too many people gathered too fast. People were crushed and we had to shut down the distribution (Ozen and Krishnamurthy, 2017).” After several days a new RC design was implemented where the crowd was more strictly controlled, and victims were accepted one by one to the RC. As in the example of Haiti relief distribution, RC design and staffing changes can help manage the crowds better, and decrease victim waiting times. In this paper we focus on the problem of how to design and setup RCs to increase distribution efficiency.

Previous research on relief distribution focused on procuring and transporting the relief supplies to the disaster affected region, and lacked focus on the distribution of the supplies at the disaster site. The current relief distribution practice is based on prior experience of relief agencies. To understand the current practice and the challenges faced, we interviewed practitioners and organized working group sessions with representatives from the Salvation Army, Red Cross, South East Wisconsin Citizens and Organizations Active in Disasters, and Federal Emergency Management Agency (FEMA). The interviews revealed that relief agencies including FEMA, Red Cross and the Salvation Army have guidelines for setting up RCs and distributing relief (Federal Emergency Management Agency, 2008), US-AID (United States Agency for International Development, 2005) and International Federation of Red Cross (IFRC) (International Federation of Red Cross and Red Crescent Societies, 2008). However, they emphasized that the guidelines merely provide a starting point and there is currently no tool that can be used to evaluate alternative RC designs based on the unique needs of a disaster. In this paper, we aim to address this need to quantify the impact of various RC design strategies. In particular, we aim to answer the following research questions:


What is an effective way to evaluate performance of different RC layout designs?


What is the impact of various RC design decisions (RC layout, RC capacity, item assignment, use of triage) on relief distribution performance?

To answer these research questions, we model a RC as a queuing network, where each point of distribution is represented as a finite capacity queue, where victims arrive and wait to receive supplies. In this queuing network, we explicitly model the level of crowding and we propose a solution approach to effectively compute performance (throughput and waiting time) of any RC layout design using a decomposition approach. Using our solution approach, we conduct numerical experiments where we compare different layouts and investigate the components of a layout that can reduce effects of crowding. We use data from the 2015 Nepal earthquake to estimate victim arrival rates and relief item distribution rates to test our model. Our main findings are: given the same resources, how a RC is setup has a significant impact on relief distribution performance; dissipating the crowd by leveraging multiple distribution points within the RC decreases waiting times considerably; triage queues can be leveraged to balance the workload in the RC to alleviate bottlenecks; and limiting number of victims allowed within the RC at once can decrease victim waiting times, however if decreased beyond a certain level, it can impair the throughput of the RC. All of our findings point to the importance of decision support models that can help practitioners decide on RC design strategies to serve more victims in a shorter amount of time.

The rest of this paper is organized as follows. Section 2 provides a review of the literature and identifies the gaps addressed by this paper. Section 3 describes two RC layouts (adapted from FEMA and an alternative layout) and models the associated queuing networks. Section 4 describes the solution methodology including the decomposition algorithm and the analytic formulas used for solving the decomposed queues. Section 5 applies the solution methodology to numerically investigate the effects of crowding in both layouts using Nepal earthquake data and Section 6 concludes the paper.

2. Literature review

Simpson and Hancock (2009), Altay and Green (2006), Galindo and Batta (2013) and Larson et al. (2006) provide extensive surveys of contributions of the operations research literature to disaster response. These studies underline the need for more research related to the last mile operations during disaster response. One of the key activities at the last mile of disaster response operations is distributing relief items to affected victims. The relief items are first procured from strategic suppliers or pre-positioned stock. Next, they are transported through a network and finally received at the disaster region. Once the items arrive, they are distributed to the victims through RCs. Figure 1 shows this relief distribution supply chain. The prior literature has mainly focused on the procurement and transportation-related problems, while the challenges of distributing the items once they arrive at the disaster site, has received little attention. Below we provide a survey of the literature in all three areas of the relief distribution supply chain to highlight this gap.

2.1 Literature on procurement

Procuring the relief supplies is the first step in the relief distribution supply chain. Many papers in the literature focus on procurement-related problems, such as selecting suppliers, pre-positioning stock and managing inventory levels. Balcik and Ak (2014) discuss supplier selection and supplier agreements for relief supplies, Wang et al. (2015) look at pre-purchasing decisions using an options contract framework for humanitarian supply chains, and Duran et al. (2011) discuss pre-positioning of relief items for CARE International. Akkihal (2006) investigates pre-positioning of non-consumable items for humanitarian operations. Adida et al. (2011) focus specifically on stockpiling at hospitals for disaster planning and Acimovic and Goentzel (2016) investigate stockpiling decisions (for non-food items including blankets, buckets, mosquito nets) using stochastic optimization models. Sodhi and Tang (2014) focus on integrating local micro-retailers while solving the pre-positioning problem to support local economies affected by disasters.

Other papers study the inventory management problem related to disasters. Whybark (2007) discusses the challenges related to inventory management following a disaster. Das and Hanaoka (2014) and Ozbay and Ozguven (2008) present stochastic inventory replenishment models for inventory control. Natarajan and Swaminathan (2014) construct a multi-period stochastic inventory model with financial constraints to investigate the impact funding uncertainty on relief item inventories. Wu et al. (2010) combine inventory management with a demand forecasting approach, to improve the forecasts along with the emergency supplies management system. Xu et al. (2010), Deqiang et al. (2011), Zhao and Cao (2015), Sheu (2010) and Li (2010) discuss forecasting approaches applicable in a disaster setting. Taskin and Lodree (2011) use forecasting data for hurricanes to manage emergency supply.

2.2 Literature on transportation

Once the relief items are procured, the next step is to transport the supplies from the stock points to the last mile distribution points, which usually involves a multi-tier supply chain. The problems faced at this stage include, deciding on facility (warehouse, staging area) locations and capacity and transportation routes between these facilities and the final distribution points. Jia et al. (2007) provide a survey for facility location problems and proposes a model specifically for large-scale emergencies. The proposed model can be cast as a generalization of the covering, P-median, and P-center models taking into account characteristics of large-scale emergencies. Balcik and Beamon (2008) analyze a distribution network design problem, which combines facility location and stock pre-positioning decisions. Hong et al. (2015) use a stochastic model to determine the size and the location of the response facilities and the inventory levels of relief supplies at each facility. Iakovou et al. (1996) focus on finding the optimal location for cleanup equipment for responding to possible oil spills using a linear integer programming model, while Batta and Mannur (1990) study a covering-location problem for ambulances and fire trucks.

Barbarosoglu and Arda (2004) model the transportation of first-aid commodities during a disaster as a multi-commodity, multi-modal two-stage stochastic program. Campbell et al. (2008) analyze routing problems relevant to disaster relief operations using the objective of minimizing arrival times. Manoj et al. (2016) combine decision for locating staging areas, assigning inventory to staging locations and routing trucks from staging areas to distribution sites. Han et al. (2011) use an optimization method to decide on warehouse locations, fleet routing and scheduling.

2.3 Literature on distribution

The last step in the relief distribution supply chain is to distribute the relief items to victims at the disaster regions. Relief items are distributed via temporary structures setup by relief agencies called RCs. The location of the RCs, the design of the RCs, its staffing and operational policies are all problems that require attention to make this process efficient. To the best of our knowledge, Roy et al. (2011) and Ozen and Krishnamurthy (2013) are the only prior work in the literature that focus on relief supply distribution at the disaster site. These papers both highlight the importance of RC design and emphasize the impact large crowds can have on relief distribution performance. Roy et al. (2011) compare different layouts using a network of M/G/C/C queues. Our work differs from Roy et al. (2011) by using a different queuing approach, where the blocking relationship between queues in the network is captured explicitly. Therefore, the model proposed in this paper can provide more realistic performance measure estimates. This work also differs from Ozen and Krishnamurthy (2013) by modeling effects of crowding using state dependent queues. Due to these enhancements, new analytic formulas are derived in this paper to solve the resulting queuing network.

In addition to the theoretical contributions, this paper also contributes to the practice by evaluating layouts that are currently used in the field based on FEMA guidelines and by proposing a new layout design to improve the current practice. This paper also uses Nepal earthquake data, to numerically evaluate the performance of both layouts and identifies trade-offs in design that can help practitioners increase relief distribution efficiency. Hence, this work aims to bridge the gap between theory and practice of relief distribution.

2.4 Queuing network literature

We model victims arriving and waiting at an RC to receive relief items as a network of finite capacity queues with state dependent service rates. Prior studies on finite capacity queuing networks with blocking have used decomposition methods. These decomposition methods are based on the idea of decomposing the system down to two-queue pairs, determining sub-system parameters and ensuring consistency and conservation of flow. Dallery and Frein (1993), Gershwin (1987), Brandwajn and Jow (1988), Altiok (1982) provide three popular decomposition methods based on three approaches: downstream blocking, upstream starvation and both blocking and starvation (symmetric view). However, once the system is decomposed into individual queues, the resulting finite capacity state dependent queues need to be analyzed. We contribute to queuing literature by providing analytic solutions for the state dependent finite capacity queues. This contribution is to the subset of the queuing literature that use decomposition methods to analyze more complex systems such as those with unreliable servers (Bihan and Dallery, 2000), closed-loop networks (Maggio et al., 2009), rework loops (Li, 2004) and assembly-type systems (Gershwin and Burman, 2000).

3. The RC model

RCs have points of distribution (PODs), triage queues, volunteers and inventory and replenishment areas. The efficiency of RC operations is impacted by the following criteria:

  1. Number of PODs: the number of PODs to open in a RC can depend on number available volunteers, volume of relief supplies available at the RC, the area of the RC and the number of victims the RC is serving.

  2. Item assignment to PODs: the decision on how to assign the relief supplies to the PODs, will depend on the types of items and the total number of items being distributed, as well as the packaging of the delivered items. As an example, tarpaulins are large and bulky, hence it might be preferable to distribute them separately. Some items such as cleaning supplies come as a set and are distributed together. The assignment decision of items to PODs, can also impact the time it takes to distribute the items. If an RC decides to distribute tarpaulins with cleaning sets at a POD, the time it takes to distribute items from that POD will increase (to the summation of both service times), compared to distributing these items separately. We will refer to this phenomenon as the additive service time assumption.

  3. Physical capacity of the RC: the term physical RC capacity, refers to the number of victims that can be physically present in the RC at a time. Note that this is different than the RC throughput. This physical capacity will depend on the RC area, the number of victims the RC is serving, the number of available volunteers and the level of available supply. The decision parameter here is the number of victims allowed to queue within the RC at once. Note that, increasing the allowed RC capacity will increase the congestion within the RC, and can increase waiting times of victims and decrease throughput. We investigate this trade-off using our model.

  4. Triage queue: the triage queue is used to check victim information prior to item distribution. Note that triage, as it is used in relief distribution, is different from how it is used in medical emergency response (Iserson and Moskop, 2007; Ergun et al., 2014). In medical emergency response the aim is to prioritize victims according to the urgency of the problem. In contrast, at an RC, the triage is used to check information and route the victims to the correct locations to minimize congestion. How much service should be provided at the triage queue is an important decision parameter for RC design.

  5. Walkways in the RC: the victims access the RC using walkways, move between PODs via walkways and exits the RC via walkways. When congestion in a walkway increases, the rate at which victims can move decreases significantly. This phenomenon is known as effects of crowding and has been quantified via empirical studies (Tregenza, 1976). The number, size and length of the walkways will all impact how fast victims can move within the RC. Moreover, expanding the allowed physical capacity of the RC by accepting more victims to access the RC at once, will increase crowding and slow down movement of victims. In fact, this was the reason why the Salvation Army Port-au-Prince distribution came to a halt.

In order to capture the effects of these design criteria on the performance of the RC, we model the RCs as a network of finite capacity queues where walkways, PODs and the triage are all represented as separate queues in the network. The number of PODs and the number of walkways will impact the number of queues in the network. The assignment of items to PODs and the tasks associated with triage will impact the distribution speed at these queues. The walkway queue service rates are dependent on the congestion in the RC and are modeled using state dependent finite capacity queues. Next we discuss two layouts (one adapted from FEMA and one alternative layout we designed) using these decision criteria.

The Guide to PODs by FEMA and US Army Corps of Engineers (Federal Emergency Management Agency, 2008) provide one of the most detailed relief distribution guidelines for the current practice. Figure 2 shows the design and summarizes the parameters provided by FEMA (Federal Emergency Management Agency, 2008).

In this RC layout, there are four PODs staffed with a volunteer each, all PODs distribute all four of the available items and all victims access the RC using the same walkway. In the aftermath of a disaster, the number of victims that need supplies is often higher than the available distribution capacity. In a design where all victims converge to queue along the same entry walkway, congestion will increase and slow down the movement of the victims. Moreover, due to all PODs distributing all items, the service time at each POD is equal, and higher compared to only distributing a subset of items at a POD.

Based on these observations and the decision criteria discussed previously, we design an alternative RC layout shown in Figure 3. In this alternative layout, we keep the overall RC area, the total number of PODs, the total number of volunteers and the walkway widths equal to the FEMA guidelines. This allows us to separate design decisions from effects of resource availability. In the alternative design, we change the item assignment to PODs and increase the number of walkways in the RC. Our aim here is to: decrease the service time at each POD by distributing a subset of items (due to the additive service time assumption); and create alternative routes for victims which can “thin” the crowd within the RC, and as a result can increase the service rate of the walkways (note the initial division of crowd to two). Note that, this design may effectively spread the crowd, however, it also results in victims going through a larger network of queues. We aim to investigate the trade-off between victims going through more queues which are less congested, and victims going through less queues which are more congested. Moreover, in our numerical studies we vary the physical RC capacity and work assigned to triage queue and investigate the impact these decisions have on both RC layouts. In particular we aim to investigate the impact of the following design decisions: what is the impact of crowd dissipating layout designs on victim waiting times and RC throughput? What is the impact of item assignment to PODs and number of walkways in an RC design on performance? What is the impact of the physical RC capacity on performance? How can triaging be utilized to improve RC performance? Our methodology can be used to evaluate a broad class of layouts in addition to the ones shown in Figures 2 and 3. In Section 5 we show how our method can be used to evaluate other layout designs.

3.1 Queuing model representation of the RC designs

We model each RC design, as a queuing network of finite capacity queues. We explicitly model both the walkways that provide access to the PODs and the PODs themselves as queues in the network. We model the walkways as finite capacity queues with state dependent service rates, where the rate of service is dependent on the speed at which pedestrians are able to move on the walkway. As the number of victims on the walkway (congestion) increases, the movement speed will slow down. We use empirical studies by Tregenza (1976) to model this relationship. We model the PODs as finite capacity queues with steady service rates where the service time only depends on the items being distributed.

Let the set i=1, 2, i, ..., N represent the queues in the RC and Ki represent the finite capacity of each queue. Also let µi(ki) and µi represent the exponential service rate parameters for walkways and PODs, respectively, where ki=1, 2, , ki, ..., Ki.

For analytic tractability, we assume that the external arrival rate of victims to the RC follows a Poisson process with parameter λ. Although, the arrival rate to the RC will be the same regardless of the design, the arrival rate to a particular POD depends on the set of items available at that POD and the demand for that subset of items. For the layout in Figure 3, a victim needing all items will go through both PODs 1 and 2, while a victim needing only item 1 will only visit POD 1. The queuing model therefore allows for converging and diverging flow of victims while keeping the overall flow of victims consistent.

Figure 4 depicts the resulting queuing network for the RC layout shown in Figure 2 which represents the current practice. We only model the first POD based on the assumption that victims going through the other three PODs will experience similar waiting times. Modeling the first POD alone will provide an optimistic view of current practice, since we omit the additional congestion experienced by victims receiving service from the PODs at the end of the linear placement (due to going through a higher number of queues). Figure 5 depicts the resulting queuing network from Figure 3 for a victim requiring all items. By leveraging the symmetry of the RC design we only model PODs 1 and 2.

Note that in the layout design given by Figure 2 (queuing network representation shown in Figure 4) all victims, irrespective of the set of items they require, use the same walkway to access the RC and the same long walkway is used to access each POD. In this design, even though victims go through less number of queues in total, the queues they go through are more congested. Moreover, a victim needing a single item or all items will go through the same network of queues. In the proposed design (see Figure 3 for layout and Figure 5 for queuing network), there are two entry points to the RC following the triage queue aiming to divide the crowd from the beginning. In addition, item assignment to each POD is varied so that a victim needing a single item need not go through all queues in the RC, and a victim needing all items will go through less congested queues. Moreover, item assignment to PODs, not only creates alternative routes for victims but also can decrease the service time at the PODs. By distributing fewer number of items, distribution time can be shortened, resulting in faster queues. In Section 5 we consider various item assignments to analyze the impact of both varied routing and service rate balancing can have on relief distribution performance.

It is possible to analyze the layout in Figure 3, as a multi-class queuing network. However, in our solution approach, we use input aggregation and calculate the arrival rate and service rate at each queue for a typical aggregate victim. See Whitt (1983) for the details and accuracy of the aggregation approach.

Next we discuss the methodology proposed to solve the resulting queuing networks. Using this methodology, we estimate the performance measures (average time a victim going through the RC experiences and throughput of the RC) and we determine them using a decomposition methodology since exact analytic solutions are hard to obtain. The overall performance measure we use (victim suffering objective) combines these two measures.

4. Performance evaluation

The RC designs result in open queuing networks with finite capacity queues and blocking after service. Let ki denote the number of victims at queue i and let P=(k1, k2, , ki, …, kN) represent the steady state probability for the network. To estimate the network’s performance measures, we first need to solve this queuing network for its steady state probabilities, which we do so following a three step process. We first decompose the network of queues into two-queue subsystems, each characterized with an effective arrival rate and a two-phase Coxian effective service rate. We then analyze each sub-system in isolation using the analytic formulas derived in Section 4.1. Next, we use a recursive algorithm (see Section 4.2), to link back the network and capture the inter-dependencies via blocking probability estimations to calculate accurate performance measures for the network.

Figure 6 shows the network after decomposition, where each sub-system is composed of a queue and its downstream, to take into account the probability of the upstream queue being blocked. The two main assumptions made during decomposition are: a queue may only be blocked by its immediate successor queue; and the arrival rate to each queue is Poisson. Each individual queue obtained from the decomposition is defined by its effective arrival rate and the effective service rate. The effective arrival rate must ensure that the system throughput stays constant, while the effective service rate must account for the probability of blocking and the associated blocking service time (see Figure 6).

Blocking in the system occurs if, following a service completion at queue i, the downstream queue i+1 is at capacity. Note that the downstream queue being at capacity does not always imply blocking. To distinguish the two states, we increase the queue capacities by 1 to represent the blocked victim in queue. Then the probability of blocking for queue i can be written as ai=P(ki=Ki+1) for i=2, 3, , N. Next we describe the effective arrival and service rate calculations.

The external arrival rate to the network is known and is assumed to be Poisson with parameter λ. Then the throughput of the network can be written as TH1=λ×(1−a1), where a1 is the probability that queue 1 is at capacity. This implies that the RC is unable to provide service to λ×a1 victims. Then, the effective arrival rate to each subsequent queue is defined as λi=TH1/(1−ai), i=2, , N, where ai, i=2, , N is the probability queue i blocks its upstream.

The service time at each queue is assumed to be exponential with a fixed rate for the distribution PODs, and state dependent exponential for queues representing RC walkways. Due to the memory-less property of the exponential distribution, the time a station is blocked can also be characterized by the exponential distribution of its downstream. This phenomenon is shown in Figure 7, where queue i is blocked with probability ai+1, and the victim blocked waits until a service completion happens at the downstream queue, which is exponential with parameter µi+1. This implies that effective service time at queue i has a two-phase Coxian distribution, where the first phase represents the actual service time and the second phase represents the time spent blocked, and ai+1 represents the probability of blocking. Note that all queues in the system besides the final queue can be characterized as an M/C2(ki)/1/Ki queue. Also note that, a queue with service time that is independent of the state ki, is a special case of a queue with state dependent service times. Since the last queue cannot be blocked it is characterized as an M/M(ki)/1/Ki queue. Next we derive the formulas required to solve M/C2(ki)/1/Ki. The derivation of formulas to solve M/M(ki)/1/Ki queues is given in the Appendix, Section A.1.

4.1 Analysis of queues in isolation

We defined each queue in the network as a standalone queue with Poisson arrivals with rate λi=TH1/(1−ai) and a two-phase Coxian distributed service time with parameters µi and (µi+1 × ai+1) characterizing the first and second phases of the service time, respectively. To calculate the steady state probabilities of the M/C2(ki)/1/K queue in isolation, we analyze the Markov Chain embedded at departure epochs. Let π(ki),ki=0, 1, , Ki−1 denote the steady state probabilities of the embedded Markov chain and let P(ki), ki=0, 1, 2, , Ki represent the steady state probabilities of the M/C2(ki)/1/Ki queue. Note that for simplicity the subscript i representing a specific queue in the network will be suppressed throughout the remainder of this section.

Let Xk denote the state of the Markov chain representing the number of victims present in the queue right after the kth departure and let qk,n represent the probability that, following a departure leaving the queue in state k, exactly n victims arrive during the service time and let rk,n represent, r k , n = j = n q k , j . Then, the transition probability matrix for the Markov chain embedded at departure instances can be written as:

Q = [ q 0,0 q 0,1 q 0,2 q 0, K 2 r 0, K 1 q 1,0 q 1,1 q 1,2 q 1, K 2 r 1, K 1 0 q 2,0 q 2,1 q 2, K 3 r 2, K 2 0 0 0 q K 1,0 r K 1,1 ]

To estimate qk,n, we repeatedly use the memory-less property of the exponential distribution along with the probability representing the minimum of two exponential random variables, where µd(Kd+1) represents the service time of the downstream queue when it blocks its upstream queue. Equation (1) is based on estimating the probability that the service time ends, before an arrival can happen, whether the queue is blocked or not. Equation (2) is based on computing the probabilities that, if the queue is not blocked, an arrival happens before the service time ends, and if the queue is blocked, an arrival happens either during the first service time or the blocked service time. Finally, Equation (3) generalizes this idea to the case where, multiple arrivals happen prior to a service completion. If the queue is not-blocked, n arrivals should happen prior to the service completion. And the service rate needs to be updated after each arrival due to having a state dependent service rate. If the queue is blocked, then a total of n arrivals in any combination can happen during the first and the blocked service times. Again we track the combination, as well as the sequence of arrivals, to update the service rate parameter accordingly.

For k⩾0 and n=0:

(1) q n , k = ( 1 a ) ( μ ( k ) λ + μ ( k ) ) + a ( μ ( k ) λ + μ ( k ) ) ( μ d ( K d + 1 ) λ + μ d ( K d + 1 ) )

For k⩾0 and n=1:

(2) q n , k = ( 1 a ) ( λ λ + μ ( k ) ) ( μ ( k + 1 ) λ + μ ( k + 1 ) ) + a [ ( λ λ + μ ( k ) μ ( k + 1 ) λ + μ ( k + 1 ) μ d ( K d + 1 ) λ + μ d ( K d + 1 ) ) + ( μ ( k ) λ + μ ( k ) λ λ + μ d ( K d + 1 ) μ d ( K d + 1 ) λ + μ d ( K d + 1 ) ) ]

For k⩾0 and n>1:

(3) q k , n = ( 1 a ) j = k k + n 1 λ λ + μ ( j ) μ ( n + k ) λ + μ ( n + k ) + a b=1 n [ ( j = k k + b 1 λ λ + μ ( j ) ) μ ( k + b ) λ + μ ( k + b ) ( j = k + b k + n 1 λ λ + μ d ( K d + 1 ) ) μ d ( K d + 1 ) λ + μ d ( K d + 1 ) ] + a [ ( j = k k + n 1 λ λ + μ ( j ) ) μ ( n + k ) λ + μ ( n + k ) μ d ( K d + 1 ) λ + μ d ( K d + 1 ) ] + a [ μ ( k ) λ + μ ( k ) ( j = k k + n 1 λ λ + μ d ( K d + 1 ) ) μ d ( K d + 1 ) λ + μ d ( K d + 1 ) ]

Now that the matrix Q is completely defined, we solve the equation π=πQ to determine the steady state probabilities of the Markov chain, π:

(4) π ( k ) = P ( k ) 1 P ( K ) for k=0,1,2, , K 1 λ ( 1 P ( K ) ) = i=1 K μ e f f ( i ) P ( i )

Then using π(k), k=1, 2, …, K−1 we get the steady state probabilities P(k) using Equation (4), where µeff(k) is defined as ((1)/(μ(k)) + (ad)/(μd(Kd + 1))) − 1. Then we can use Equation (5) to get the steady state probabilities from departure epoch probabilities:

(5) P ( k ) = π ( k ) π(0) + ρ , k=0,1,2, , K 1 P ( K ) = ρ 1 + P ( 0 ) ρ

Finally, we calculate the throughput, average number of victims in each queue and the average time spent at each queue using Equation (6) for each queue i. THi in Equation (6), represents the throughput of queue i and is equal to the number of victims that enter the RC from queue 1. Li represents the average number of victims in queue i and is calculated by taking an expectation over the state space for queue i. Finally W ¯ i represents the average time a victim spends in queue i and is computed using Little’s Law:

(6) T H i = λ i × ( 1 a i ) L i = j=0 K i P i ( j ) j , i = 1 , 2 , ... , N j = 0 W ¯ i = L i / T H , i=1,2 , ..., N

4.2 Iterative algorithm to link queues

The steady state probabilities of any queue i in the network is dependent on the probability that their downstream queue (i+1) is full, and the probability that the first queue is full (since it governs the overall system throughput). Formally, the probability ai is a function of a1 and ai+1, ∀i=1, 2, 3, …, N−1 and aN is solely dependent on a1. Therefore, to accurately estimate the steady state probabilities of each queue, we solve each queue in isolation assuming starting blocking probabilities (using equations derived in Section 4) and iterate until all blocking probability estimates converge. The iterative algorithm is summarized below:

Inputs: λ, Ki, µi(ki), for all i=1, 2, …, N, initial estimates a i ( 0 ) , i=1,2,.. N and Step1: For m=1, 2,, M execute steps 2-4.

Step2: For queue i=N:

2.1 Analyze queue N as an M/M(kN)/1/KN queue with arrival rate λN, service rate µN(kN) and capacity KN and obtain steady state probabilities PN(kN=j), j=1, 2, , KN.

2.2 Update arrival rate λN with revised estimate of aN=PN(kN=KN).

2.3 Iterate until estimate of PN(kN=KN) converges.

2.4 Let a N ( m ) correspond to this estimate of PN(kN=KN).

Step3: For queues i=N−1, N−2,, 2:

3.1 Analyze queue i as an M/C2(k)/1/Ki queue with arrival rate λi, service rate µi(ki) and capacity Ki+1 and obtain steady state probabilities Pi(ki=j), j=1, 2, …, Ki+1.

3.2 Update arrival rate λi with revised estimate of ai=Pi(ki=Ki+1).

3.3 Iterate until estimate of Pi(ki=Ki+1) converges.

3.4 Let a i ( m ) correspond to this estimate of Pi(ki=Ki+1).

Step4: For queue i=1:

4.1 Analyze queue 1 as an M/C2(k)/1/K1 queue with arrival rate λ, service rate µ1(k1) and capacity K1+1 and obtain steady state probabilities:

P 1 ( k 1 = j ) , j=1,2, , K 1 + 1.

4.2 Compute the absolute difference δ = | a 1 ( m ) P 1 ( k 1 = K 1 + 1 ) | .

4.3 δ⩽ϵ, go to step 5.

4.4 If δ>ϵ, set m=m+1 and repeat steps 2-4.

Step 5: Record steady state probabilities Pi(ki), ki=1, 2,, Ki+1, ∀i.

Step 6: Compute performance measures using Equation (6).

For a given value of a1, we start with the last queue N. We guess an initial value for aN and solve this queue for its steady state probabilities and obtain a revised estimate of aN=PN(nN=KN+1). We iterate over values of aN until the estimate of aN is consistent with λN. Using the solution for queue N, we analyze queue N−1. This process continues until we get a revised estimate of a1. If the error between the current and previous estimate of a1 is greater than, we repeat the solution process for queue N through queue 1, else we exit the loops. Note that for the search for a1 we use a grid search, while the search for all other ai follows a fixed point iteration algorithm (Conte and de Boor, 1980). In the global loop we apply a grid search over values of a1 as it allows us to solve for cases where the external arrival rate greatly exceeds the achievable throughput of the network.

4.3 Objective function for measuring performance of an RC

The two main RC performance measures are; average time to service a victim and the hourly throughput of the RC, since minimizing suffering of victims requires both reaching to all victims and reaching them as quickly as possible. Therefore, maximizing the RC throughput and minimizing the waiting time of victims is the most preferred outcome, however there are decisions that require trade-offs between these two measures. Moreover, the evaluation of the throughput measure is dependent on the total population the RC needs to serve. The throughput measure may be critical if the target population is large and many victims face the risk of being denied service. Therefore a single objective function capturing these trade-offs is needed. There has been many disaster response research on appropriate humanitarian objective functions such as minimizing waiting time, minimizing unsatisfied demand and ensuring equity. Among all proposed objectives, the human suffering objective proposed by Holguin-Veras et al. (2013) is a good candidate for evaluating RCs because it quantifies victim suffering using willingness to pay (WTP) functions for a given relief item, which increases with the time of deprivation. The objective function takes into account: the waiting time experienced by victims (W); and the hourly rate at which the RC serves victims (TH); and if supply (or capacity) is less than demand, the number of victims the RC was unable to serve. Equation (7) gives the objective function. The parameters used are: P the affected population in need of items, ω per hour cost of a volunteer, n the number of volunteers and τ the time period in hours.

The objective function is written as:

(7) γ ( W ¯ , T H , ω , τ ) = ϕ ( W ¯ , T H , τ ) + φ ( W ¯ , T H , τ ) + κ ( ω , τ , n )

The components of the objective function ϕ ( W ¯ , T H ) , φ ( W ¯ , T H ) and κ(ω, τ, n) represent the suffering function for victims (measured by waiting time (W)) that receive service at the RC (characterized by throughput (TH)), the suffering function for victims that could not receive service at the RC and the labor cost function, respectively. The logistics cost for the RC is assumed to be the staffing cost alone, since we assume that the items to be distributed are already in the affected area at the time of RC setup. The functions ϕ and φ have the exponential structure following the WTP function reported in Holguin-Veras et al. (2013). Equations (8) and (9) show the exponential structure of the WTP function and the parameters a and b are obtained from Holguin-Veras et al. (2013). In Section 7.3 we perform sensitivity analysis to understand the effects of this particular parameter choice on RC decisions. Equation (10) computes the logistics cost of an RC by multiplying number of volunteers with the hourly cost of a volunteer and the total hours worked. Note that the objective function uses the waiting time and throughput of the RC as inputs, which are the outputs of our model that evaluates a particular RC design, obtained using the algorithm described in Section 4.

Suffering function for victims that receive service at the RC:

(8) ϕ ( W ¯ , T H , τ ) = [ exp ( a + ( b × W ¯ ) ) exp ( a ) ] × T H × τ

Suffering function for victims that do not receive service at the RC:

(9) φ ( W ¯ , T H , τ ) = [ exp ( a + ( b × τ ) ) exp ( a ) ] × ( ( P ( T H × τ ) ) )

Labor cost function at the RC:

(10) κ ( ω , τ , n ) = ω × n × τ

5. Case study: Nepal Earthquake (2015)

In this section, we discuss the application of our model using the 2015 Nepal earthquake as a case study. The Nepal earthquake is a good case for analyzing relief distribution, because it resulted in the displacement of three million people, making relief item distribution a top priority during the immediate response. Moreover, due to inadequate infrastructure and limited number of volunteers and agencies, the efficient distribution of relief was a significant challenge. Using Nepal data from Humanitarian Data Exchange Database (2016), we estimate arrival and service rates for RCs and use our model to evaluate the RC designs given in Figures 2 and 3. We use numerical experiments to: evaluate the impact of item assignment decisions to PODs; quantify the trade-off between going through more queues vs going through less queues with more congestion; evaluate the impact of the decision to limit the physical RC capacity on RC performance; and investigate alternative uses for the triage queue. We use the multi-criteria suffering objective function to evaluate the RC designs and operational strategies. Sections 5.2, 5.3 and 5.4 present numerical analysis that answer these research questions.

5.1 Parameter estimation

To solve the queuing model and estimate RC performance measures, arrival rate to the RC and service rate at the PODs are required input parameters. We use data from the Humanitarian Data Exchange Database, which is an open platform for humanitarian data sharing (https://data.humdata.org/) to estimate these key parameters. We use two data sets: the distributed items data set; and the demographics and casualty information data set. The distributed items data set provides information on types of items distributed, the distributing agency, number of households served and the dates the distribution took place. The data set is prepared by the Global Shelter Cluster (GSC) which is a coordinating agency that brings together 35 global partners to coordinate and respond to disasters and conflict situations. The GSC is co-chaired by the IFRC and the United Nations High Commissioner for Refugees. The demographics and casualty data set provides population, number of households, total number of deaths and injured victims per district. The data set is provided by the Nepal United Nations Office for the Coordination of Humanitarian Affairs with data contributions from the Nepal Ministry of Home Affairs and the Nepali Police.

The data set on distributed supplies contains RCs that were active for seven months, distributing different types of items. Our study focuses on the items that support the immediate needs of victims following an earthquake such as tarpaulins, blankets, kitchen supplies (for water and food) and hygiene kits. Therefore we limit the analysis time frame to the first three weeks. Note that the distribution of perishable and repeat items such as water and food are not part of this analysis or data set. Next we describe the process used to estimate the arrival and service rates from the raw data. Note that as described in Section 4, these form crucial inputs to the queuing model.

5.1.1 Arrival rate estimation

The arrival rate per RC per day is estimated for each district, as the total survived population divided by the total RC days of active item distribution. To calculate the total survived population we use the demographics information of each district, and assume all earthquake survivors needed critical relief supplies. This gives us the total demand rate. Estimating the arrival rate per hour per RC was challenging. To do so, we used the total RC days parameter for each district, which refers to the total number of days an RC was open, summed over all RCs. This parameter allows us to convert total demand to an hourly arrival rate for an RC, by assuming that the distribution was carried out for 16 hours each day.

While the demand per district is closely related to its population, the total RC days varies between districts due to the accessibility of the district, and the differences in the number of agencies operating in that district. Hence, the arrival rate to an RC varies between different districts. Table I provides the resulting arrival rate per RC, per hour, for the top ten most affected districts. In our analysis we first discuss insights with respect to Bhaktapur, since its arrival rate represents the average. Subsequently in Section 5.5 we discuss insights for Nuwakot and Kathmandu to capture the effects of increasing arrival rates.

5.1.2 Service time estimation

There are no standard times provided for disaster relief item distribution. Therefore, we estimate the time it takes to distribute an item of a particular type, using the data set of distributed supplies. The data set contains the number of households served in a given day for an agency, district, item triplet. In our analysis, we distinguish item types but we do not distinguish districts or agencies since the time to distribute a tarpaulin will be the same regardless of location or agency. Finally, we estimate number of households served per day for a given item type. The number of items required of a given Table I: Nepal Data: Arrival Rate Estimation type, per household, is different for different items. For example every household (made up of five people on average) needs one tarpaulin while needing two blankets. To calculate the distribution time per item, we need number of items (for a given item type) a household needs. We obtain this information from Red Cross Nepal Emergency Appeal Operation Update Reports as shown in Table II. Then we calculate the time it takes to distribute a specific item. Table III provides the average service rates per item calculated over the ten most affected districts.

5.2 Comparing RC layout designs

In the first set of experiments we compare layouts in Figure 2 (representing current practice) and Figure 3 (representing alternative design) and two alternatives to evaluate the separate impact of item assignment decisions and crowd dissipation strategies. Layout alternative 1 is the queuing representation of current practice adapted from FEMA guidelines, as previously given by Figures 2 and 4. Layout alternative 4 leverages the multi-access points to the RC idea, however distributes all items from all PODs as described in FEMA guidelines. Layout alternatives 1 and 4 are shown side by side in Figure 8. Both these layouts result in a 4 queue network composed of the entry walkway to triage, the triage queue, entry walkway to POD 1 and service queue for POD 1.

Layout alternatives 2 and 3 represent the queuing network representation of the alternative design, where in alternative 2, POD 1 distributes tarpaulins and kitchen sets and POD 2 distributes blankets and sleeping mats. Whereas in design 3, POD 1 distributes tarpaulins and blankets while POD 2 distributes kitchen sets and sleeping mats. Layout alternatives 2 and 3 are shown side by side in Figure 9. Both these layouts result in a seven queue network composed of the entry walkway to triage, the triage queue, entry walkway to POD 1, service queue for POD 1, connecting walkway to POD 2, entry walkway to POD 2 and service queue for POD 2. Note that in all four designs, the total area of the RC and the walkway widths for pedestrians are the same, as specified by FEMA documentation.

Comparing layouts 1 and 4, we aim to investigate the impact of dissipating the crowd, on RC performance. In layout design 1, the entry walkway will be congested due to the single access point. In layout 4, we remedy this by creating multi-access points in the layout. Next, we compare layouts 1 and 2. In layout 2, we assign a subset of items to each POD, and dissipate the crowd entering the RC. By comparing these two layouts, we aim to compare the combined impact of item assignment and crowd dissipation. Assigning a subset of items to each POD, is expected to decrease the distribution time. Note that in layout 2, victims will be going through less congested queues and receive faster service at the PODs. However, they will have to go through a higher number of queues (7 queues as opposed to 4). Through numerical experiments, we aim to compare the impact between going through a shorter network of queues, and going through a larger network with less congested, faster queues. Finally, we compare layouts 2 and 3 to investigate the impact of item assignment alone on RC performance. In layout 3, distributing tarpaulins and blankets together will decrease the service time at POD 1. This results in lower blocking in the system, and therefore a higher throughput. However, a higher throughput, increases the number of victims present in the RC at once (hence increases the congestion), which in turn, creates higher victim waiting times.

Table IV provides the parameters for the numerical experiment. We use the arrival rate of Bhaktapur (58/hour) as a representation of the average arrival rate. Note that the overall arrival to the RC is the same for all designs. However, the arrival rate to each queue may differ depending on the demand for items being distributed at a particular POD. To calculate arrival rate per POD, we use the data provided in Table II. We use the distribution rate per item as given in Table III and assume the rates are additive. Additivity is a reasonable assumption (if not necessary) for items like tarpaulins and kitchen sets where handling takes considerable time. Note that for PODs at which not all arriving victims need all items, we use a weighted average to calculate the service rate. For the walkway service rates we use one of the empirical curves that characterize the exponential relationship between crowd density and walking speed given by Tregenza (1976). The exact exponential decay function used in these experiments is given in Table IV. To compute the suffering objective, we assume 16 hours of operation per day and for victims unable to receive service at the RC, we assume a 24 hour daily deprivation. For the labor cost we assume $10/hour per volunteer. Note that since the number of volunteers in all layout designs are the same, the layout design evaluations are insensitive to this assumption.

Table V provides the results for the four different RC designs. In the table we report the average time experienced by a victim ( W ¯ ) , the total hourly throughput (TH) of the RC, and the suffering objective. Note that in the Appendix in Section A.2 we validate our solution approach by comparing the results obtained via the decomposition algorithm with simulation results.

Comparing the results for layouts 1 and 4, we observe that dissipating the crowd by creating multiple access points, resulted in an improvement in victim waiting times. Here it is important to note that, the waiting time estimate for layout 1 is for victims going through POD 1 alone. Victims going through PODs 2, 3 and 4 in layout 1 will experience higher waiting times, especially under congestion. Therefore, we conclude that when congestion is high, dissipating the crowd can have a considerable impact on victim waiting times. Since the item assignment to PODs is the same for layouts 1 and 4, the bottleneck service rate remains the same. And that is why, the throughput for both layouts is equal. Next, we compare results for layouts 1 and 2, and observe a significant improvement for both victim waiting times and the RC throughput. This shows that changing item assignment to PODs combined with crowd dissipation strategies (under the additive service time assumption) can drastically improve RC performance. This is because in layout 2, congestion in walkways is lower (crowd dissipation) and service time at the PODs is lower (item assignment). Hence, all queues in the network can flow faster, resulting in increased throughput and decreased waiting times. Finally, we compare the results for layouts 2 and 3. We observe that, item assignment can also impact both the throughput and service time measures.

The results obtained in this experiment set lead us to several insights. First, going through a larger number of queues can decrease the time victims spend in the RC. This is a counter intuitive result, since waiting in more lines should result in an increase in time spent waiting. However, our numerical experiments show that, if creating more queues in the RC layout can be leveraged to decrease the congestion and balance service times, larger networks can be preferred. Second, the decision of how to assign items to PODs can have a significant impact on performance. Most often, having all items available at each POD is easier to setup. However, our results show that it is important to think about whether this decision is creating significantly higher service times at the PODs. This may not be the case for smaller items or bundled items where the additive service time assumption does not hold. However, quantifying the impact of assignment decision for cases where it will matter, can increase the awareness of practitioners making these decisions. Lastly, some decisions can have an opposite impact on throughput and waiting time measures. As an example, the throughput increase in layout 3, comes at the expense of an increased waiting time. In such cases, having a performance measure such as the suffering objective, that encompasses multiple criteria will be highly useful. Comparing the suffering function objective for layouts 2 and 3, we observe the cost of layout 3 is much lower, even though it creates a higher waiting time for victims at the RC. This is because, victims receiving their supplies waiting an additional hour has less weight compared to 13 victims (every hour) not receiving items at all. This showcases the need for humanitarian multi-criteria objective functions.

We would also like to note that in these experiments, we did not consider the impact of number of available volunteers or available RC area on RC design decisions. If there are limited number of volunteers or if the area is limited which limits the number of PODs the RC can operate, distributing all items from a single point may be the only option. However, using models like the one proposed in this paper, can estimate the impact these limitations will have on expected waiting times, and when options present themselves, these models can help make more efficient decisions.

5.3 Effects of allowed queue capacity

The physical size of the RC and the walkways are given. However, the number of victims allowed to queue inside the RC is a decision parameter. We investigate the effectiveness of this decision called allowed physical queue capacity Ki, i=1, 2, , N, as an alternative congestion control strategy. Since the actual queue capacity (for walkways) is dependent on the size of the walkway, which we do not change for these experiments, by varying Ki between 25 and 75 for all four layout designs, we impact the congestion of queues in the RC. We keep the rest of the parameters constant as given in Table IV. Figure 10 presents the resulting changes in RC throughput and average service time while Table VI reports the suffering objective results.

We observe that, by limiting the number of victims allowed to wait at each queue, we can control congestion and decrease waiting times for victims. The throughput is unaffected due to being bounded by the service rate at the distribution PODs. (Service rate for Tarpaulins is 10.7 per hour while the arrival rate to the RC is 58 victims per hour). Note that, the suffering function captures the cost of victims not served by the RC for all cases. For a setting where the bottleneck is the walkways instead of the distribution PODs, we would expect the throughput to improve as well. An example of this strategy from practice is, after the Nepal Earthquake, Red Cross relief distribution volunteers gated some of their RCs and allowed victims to enter the RC one by one. On the other hand there were many other RCs which were flooded with victims. Through decision support models, this impact could be quantified and can help humanitarian organizations make informed decisions. By quantifying the effects of allowed capacity with these experiments, we conclude that the best practice is to keep the walkway density between 0.5-0.7 victims/sqm (N=25 for walkway sizes in this experiment).

The key insight to take away from this experiment is the fact that congestion has a high impact on victim waiting times. And when possible, controlling the victim flow within the RC can help manage high waiting times. However, it is important to note that, controlling the victim flow in practice will depend on how well the RC can be monitored, the available staff to gate the RC and the particular RC layout. Moreover, the extent at which crowding happens given a victim arrival rate, will depend on the overall available RC area and walkway widths. However, given these limitations, models that can quantify impact of decisions under many different scenarios will be helpful. Next, we investigate how best to leverage the triage queue to improve performance.

5.4 Effect of triage

The triage queue is modeled in all designs and its primary function is to check identifications and route victims to correct PODs. In our next experiment, we investigate whether the triage queue can also be used to balance the queuing network and alleviate bottlenecks. To balance the workload in the RC network, we look at transferring a portion of the work performed at the PODs to the triage queue. In this scenario, the triage queue will keep the victim routing function but might also hand out a common item. We use RC layout design 2, and assume work transfer from the bottleneck queue (K, S) to the triage queue, in increasing percentages.

We observe that work transfer whenever possible can be a very effective strategy to both increase the throughput and decrease the average system time. However, it is crucial to be aware of the risk of making triage the bottleneck in this process. Deciding the amount of work to transfer (or initially assign to triage) is not trivial, mainly because it is not only dependent on the rate of service, but also dependent on the rates of arrival to each queue, and the particular RC design (area and length of queues). Hence, models and fast solution algorithms such as those proposed in Section 4, are valuable tools to relief agencies. Figure 11 presents the changes in throughput and average system time as a result of percentage work transfer from service queue distributing kitchen supplies and sleeping mats to the triage queue. The results of the experiments reveal an interesting queuing behavior. We observe that the effects of moving work to triage queue can be characterized in three zones. Zone 1 (in this experiment until 25 percent) increases throughput while decreasing system time. This is due to the common impact of bottleneck service rate increase and decreased blocking. As the transferred work amount increases in Zone 2, the throughput decreases. This decrease in throughput, decreases system congestion, resulting in lower waiting times. Finally in Zone 3, the triage becomes a major bottleneck, hurting both throughput and system time measures. Such threshold behavior is useful to know when making decisions on site to improve RC performance.

The main limitation for moving work to the triage queue is that, the optimal work distribution might not always be feasible. As an example, moving 20 percent of tarpaulin distribution time, to the triage queue is not an option. However, distributing water in addition to checking IDs can be practical, can eliminate the need for an additional queue and can speed distribution.

Figure 12 provides the suffering objective results for this experiment. Due to the limited available capacity to serve the affected crowd, the current objective puts more emphasis on higher throughput. Therefore, the decrease in throughput in zone 2 is penalized, even with the decrease in waiting time. Appendix in Section A.3 provides a sensitivity analysis on WTP function parameters used in the suffering level calculations, and concludes that the results obtained are not very sensitive to these parameters.

5.5 Performance under increasing arrival rates

In previous experiments we used the arrival rate observed at Bhaktapur. However there were districts observing much higher arrival rates due to higher earthquake damage and less available resources. In this section we run experiments for Nuwakot and Kathmandu to understand how the RC performance results would change under heavier arrival rates. Figure 13 shows the system time results for Nuwakot and Kathmandu for all four layout designs shown in Figure 8, and Table VII shows the suffering function results.

We observe that the throughput remains the same with increased arrivals due to the limit on capacity (note that the service PODs are the bottleneck). Due to finite capacity queues, arrivals beyond capacity are lost, keeping the throughput measure steady. Figure 13 shows the average system time change for increasing arrival rates. Increasing rates of arrival, increases the congestion causing the system to be slower, resulting in a much lower performance and a much higher suffering function. Strategies involving limiting the allowed queue capacity, and using layouts that leverage crowd dissipation, becomes increasingly important for cases such as Kathmandu and Nuwakot.

6. Conclusions

In this paper, we present a modeling and solution approach for RCs that distribute aid to victims at the disaster site. We model the relief distribution operations through a RC using a finite capacity queuing network, which explicitly models effects of crowding via state dependent service rates. To solve this queuing network, we derive new analytic formulas for steady state probabilities of state dependent Coxian queues, and estimate the throughput of the RC and the waiting time of victims. We use this solution methodology and Nepal earthquake data, to numerically investigate four layout designs. Using these layout designs we analyze critical decision variables that can impact RC performance including; item assignment decisions to PODs in the RC, the effects of crowding on RC performance, the impact of physical RC capacity, and the impact of alternative uses for the triage queue. We emphasize that, the methodology of this paper is not specific to the layouts discussed in this paper, and it can be used to evaluate the performance of any RC design.

Our findings from the analysis conducted have significant practical impacts. First, we find that, given a certain level of resource availability (available area, number of volunteers, number of PODs), how the RC is setup can have a significant impact on performance. Our conversations with practitioners regarding the current practice suggests that, there is a lack of awareness with respect to the impact certain decisions can have of relief distribution performance. We believe studies like this one, that quantify trade-offs and effects of design decisions can impact the way practitioners think about RC setup and operations.

We suggest that practitioners think about ways to dissipate the crowd within the RC area. We find that, even when crowd dissipating designs create more queues for victims to go through, the decrease in congestion can actually shorten the victim waiting times despite the longer route. We also suggest to be mindful of the types of items an RC will be distributing and adjusting the design of the RC accordingly. When an RC is distributing large items such as tarpaulins, one-stop-shop approaches can be more harmful than helpful for performance. We also suggest to limit physical RC capacity. This can be counter intuitive when lots of victims are waiting to access the RC. However, due to increasing congestion, allowing as many victims as possible to enter the RC at once will increase the waiting time of victims, not decrease it. Finally, we suggest practitioners to look for ways to balance the workload between the PODs in the RC. We performed numerical experiments using the triage queue to balance the workload, however this is not the only way. The findings apply to balancing the service time between PODs as well.

We believe that the findings of this paper will: create awareness for the impact design decisions have on relief distribution performance among practitioners; highlight the need for models that can quantify impact of RC design decisions; and highlight some of the essential trade-offs (such as impact of crowd dissipation, impact of item assignment, impact of allowed RC capacity) in RC design. This study also contributes to the body of work that can increase the efficiency of relief distribution during disaster response and relieve some of the suffering caused by disasters.


The relief distribution supply chain

Figure 1

The relief distribution supply chain

RC layout: current practice

Figure 2

RC layout: current practice

RC layout: 1alternative layout

Figure 3

RC layout: 1alternative layout

Current practice: queuing network representation

Figure 4

Current practice: queuing network representation

Proposed design: queuing network representation

Figure 5

Proposed design: queuing network representation

Queuing network after decomposition

Figure 6

Queuing network after decomposition

Model of a queue in isolation

Figure 7

Model of a queue in isolation

Queuing network representation of layout alternatives

Figure 8

Queuing network representation of layout alternatives

Queuing network representation of layout alternatives

Figure 9

Queuing network representation of layout alternatives

Increasing queue capacity effects of the RC

Figure 10

Increasing queue capacity effects of the RC

Work transfer to triage queue: RC performance results

Figure 11

Work transfer to triage queue: RC performance results

Numerical experiment 3: suffering objective results

Figure 12

Numerical experiment 3: suffering objective results

Average system time: increasing rates of arrival

Figure 13

Average system time: increasing rates of arrival

Willingness to pay functions

Figure A1

Willingness to pay functions

Nepal data: arrival rate estimation

District Arrival/RC day Arrival/RC hour
Bhaktapur 925 58
Dhading 465 29
Dolakha 461 29
Gorkha 168 10
Kathmandu 3,158 197
Lalitpur 643 40
Nuwakot 1,535 96
Ramechhap 368 23
Rasuwa 245 15
Sindhupalchok 138 9

Items required per household

Item Conversion parameter
Tarpaulin 1/household
Blanket 2/household
Kitchen sets 1/household
Sleeping mat 5/household

Service rates for items distributed

Item Average distribution rate
Tarpaulin 10.7/hr
Blanket 31.6/hr
Kitchen sets 10.2/hr
Sleeping mat 20.2/hr

Numerical experiment 1 parameters

Parameter Value
Arrival rate λ=58/hour (for Bhaktapur)
Service rates μT=10.7/hr, μB=31.6/hr, μK=10.2/hr, μS=20.2/hr
Queue capacity Ni=50, ∀i
Walkway service rate μw(ki)=60×exp((−c/d)×(ki)1/4), b=40/60, c=−log(b), d=(Ni)1/4

Numerical experiment 1 results

Layout design W ¯ TH ϕ ( W ¯ , T H ) φ ( W ¯ , T H ) κ(ω, τ, n) Total
Layout 1 6.63 hrs 27.14/hr $2,294 $25,744 $800 $28,838
Layout 2 2.79 hrs 33.64/hr $936 $18,424 $800 $20,160
Layout 3 3.59 hrs 46.86/hr $1,763 $3,536 $800 $6,099
Layout 4 6.39 hrs 27.14/hr $2,176 $25,744 $800 $28,720

Numerical experiment 2: suffering objective results

Layout design N = 25 N = 50 N = 75
Layout 1 $27,734 $28,838 $30,777
Layout 2 $19,713 $20,160 $20,676
Layout 3 $5,473 $6,099 $6,711
Layout 4 $27,436 $28,720 $30,670

Impact of increasing arrival rates on the suffering objective

Layout design λ=58 λ=96 λ=197
Layout 1 $28,838 $80,835 $195,932
Layout 2 $20,160 $71,984 $186,355
Layout 3 $6,099 $57,907 $172,323
Layout 4 $28,720 $80,615 $195,403

Accuracy of algorithm

Model results Simulation results
Layout design W ¯ TH W ¯ TH
Layout 1 6.63 hrs 27.14/hr 6.90 hrs 26.74/hr
Layout 2 2.79 hrs 33.64/hr 3.70 hrs 46.54/hr
Layout 3 3.59 hrs 46.86/hr 3.07 hrs 33.29/hr
Layout 4 6.39 hrs 27.14/hr 6.46 hrs 26.8/hr

Sensitivity analysis: suffering objective

Function g1 Function g2
Work transfer to triage (%) Objective value Ranking Objective value Ranking
0 $20,164 8 $1,664 8
5 $18,187 6 $1,586 7
10 $15,986 5 $1,502 5
15 $13,572 3 $1,407 3
20 $11,783 2 $1,311 2
25 $9,765 1 $1,262 1
30 $14,715 4 $1,422 4
35 $19,034 7 $1,579 6
40 $22,602 9 $1,710 9
45 $25,796 10 $1,872 10
50 $28,550 11 $2,002 11
55 $30,712 12 $2,085 12
60 $32,575 13 $2,159 13


A.1 Steady State probabilities of the M/M(k)/1/K queue

The steady state probabilities of the M/M(k)/1/K queue can be obtained by constructing and solving the balance equations given by the following equation:

(A1) P ( 0 ) λ = P ( 1 ) μ ( 1 ) P ( 1 ) λ + P ( 1 ) μ ( 1 ) = P ( 0 ) λ + P ( 2 ) μ ( 2 ) P ( 2 ) λ + P ( 1 ) μ ( 1 ) = P ( 1 ) λ + P ( 3 ) μ ( 3 ) P ( K 1 ) λ = P ( K ) μ ( K )

Solving the above equations, we can obtain:

(A2) P ( 0 ) = 1 1 + n=1 K λ n μ ( 1 ) μ ( n )
(A3) P ( k ) = P ( 0 ) λ n μ ( 1 ) μ ( k )

We then use Equation (6) to obtain the throughput (TH) and total time (W) performance measures.

A.2 Validation of the model

We provide simulation results for the experiments to validate our mathematical model and solution approach in Table AI. We use Arena simulation (version 15.0) of the queuing model for various cases where each run was for 100, 16 hour days and the results show the average of 10 replications. The error for throughput on average was between 0-2 percent and the errors for system time was on average 1-10 percent.

A.3 Impact of suffering function parameters

In the suffering objective function, we use a WTP function to express waiting (deprivation) time (both for served and not served victims) in monetary terms so that system time, throughput and logistics cost can be expressed with a single measure (dollar value). However, it is important to note that, the resulting measure is dependent on the form of the particular function used for WTP. The function we used for our numerical experiments are based on the WTP function for water defined in Holguin-Veras et al. (2013). Given that the relief items in our data set are tarpaulins and kitchen sets, which are less sensitive to deprivation we ran experiments to investigate how the objective function results change for a WTP function that is less steep. For the new WTP function we decreased parameter (b) to obtain a less steep exponential curve. Figure A1 graphs both functions g1 and g2, where g1 represents the function from Holguin-Veras et al. (2013) and function g2 represents a more gradual incline.

We re-calculated all objective functions for all experiments in Sections 5.2, 5.3 and 5.4 using the WTP function g2. The only differences in decisions we observed happened for the triage experiments. Table AII compares the objective function values for triage experiments for both WTP functions. The steepness of the function used effects the trade-off between throughput and average system time. This phenomenon changes the ranking of 5 and 35 percent work removal to triage scenarios. Under the second exponential function (g2) the decrease in system time outweighs the loss in throughput.

We conclude that the decisions are not very sensitive to the parameters of the WTP function. What is more important is to use an objective function that takes into account both victim suffering and labor costs. Moreover, the multi-criteria used by the suffering objective allows for easier decision making for relief agencies. If future empirical research provides item specific WTP parameters, we expect that these improved estimates will only provide more robust decisions.


Acimovic, J. and Goentzel, J. (2016), “Models and metrics to assess humanitarian response activity”, Journal of Operations Management, Vol. 45, pp. 11-29.

Adida, E., DeLaurentis, P.-C.C. and Lawley, M.A. (2011), “Hospital stockpiling for disaster planning”, IIE Transactions, Vol. 43 No. 5, pp. 348-362.

Akkihal, A.R. (2006), “Inventory pre-positioning for humanitarian operations”, Master’s thesis, MIT, Boston, MA.

Altay, N. and Green, W.G. (2006), “Or/ms research in disaster operations management”, European Journal of Operational Research, Vol. 175 No. 1, pp. 475-493.

Altiok, T. (1982), “Approximate analysis of exponential tandem queues with blocking”, European Journal of Operational Research, Vol. 11 No. 4, pp. 390-398.

Balcik, B. and Beamon, B.M. (2008), “Facility location in humanitarian relief”, International Journal of Logistics: Research and Applications, Vol. 11 No. 2, pp. 101-121.

Balcik, B. and Ak, D. (2014), “Supplier selection for framework agreements in humanitarian relief”, Production and Operations Management, Vol. 23 No. 6, pp. 1028-1041.

Barbarosoglu, G. and Arda, Y. (2004), “A two-stage stochastic programming framework for transportation planning in disaster response”, Journal of the Operational Research Society, Vol. 55 No. 1, pp. 43-53.

Batta, R. and Mannur, N. (1990), “Covering location models for emergency situations that require multiple response units”, Management Science, Vol. 36 No. 1, pp. 16-23.

Bihan, H.L. and Dallery, Y. (2000), “A robust decomposition method for the analysis of production lines with unreliable machines and finite buffers”, Annals of Operations Research, Vol. 93 No. 1, pp. 265-297.

Brandwajn, A. and Jow, Y.L. (1988), “An approximation method for tandem queues with blocking”, Operations Research, Vol. 36 No. 1, pp. 73-83.

Campbell, A.M., Vandenbussche, D. and Hermann, W. (2008), “Routing for relief efforts”, Transportation Science, Vol. 42 No. 2, pp. 127-145.

Conte, S.D. and de Boor, C. (1980), Elementary Numerical Analysis, McGraw-Hill, New York, NY.

Dallery, Y. and Frein, Y. (1993), “On decomposition methods for tandem queuing networks with blocking”, Operations Research, Vol. 41 No. 2, pp. 386-399.

Das, R. and Hanaoka, S. (2014), “Relief inventory modelling with stochastic lead-time and demand”, European Journal of Operational Research, Vol. 235, pp. 616-623.

Deqiang, F., Yun, L. and Changbing, L. (2011), “Forecasting the demand of emergency supplies: based on the CBR theory and BP neural network”, in Kaminishi, K., Duysters, G. and Hoyos, A.D. (Eds), Proceedings of the 8th International Conference on Innovation and Management, Kitakyushu, pp. 700-704.

Duran, S., Gutierrez, M.A. and Keskinocak, P. (2011), “Pre-positioning of emergency items worldwide for care international”, Interfaces.

Ergun, O., Gui, L., Stamm, J.L.H., Keskinocak, P. and Swann, J. (2014), “Improving humanitarian operations through technology-enabled collaboration”, Production and Operations Management, Vol. 23 No. 6, pp. 1002-1014.

Federal Emergency Management Agency (2008), IS-26 Guide to Points of Distribution, FEMA.

Galindo, G. and Batta, R. (2013), “Review of recent developments in or/ms research in disaster operations management”, European Journal of Operational Research, Vol. 230 No. 1, pp. 201-211.

Gershwin, S.B. (1987), “An efficient decomposition method for the approximate evaluation of tandem queues with finite storage space and blocking”, Operations Research, Vol. 35 No. 2, pp. 291-305.

Gershwin, S.B. and Burman, M.H. (2000), “A decomposition method for analyzing inhomogeneous assembly/disassembly systems”, Annals of Operations Research, Vol. 93 No. 1, pp. 91-115.

Han, Y., Guan, X. and Shi, L. (2011), “Optimization based method for supply location selection and routing in large scale emergency material delivery”, IEEE Transactions on Automation Science and Engineering, Vol. 8 No. 4, pp. 683-693.

Holguin-Veras, J., Perez, N., Jaller, M., Luk, N., Wassenhove, V. and Aros-Vera, F. (2013), “On the appropriate objective function for post-disaster humanitarian logistics models”, Journal of Operations Management, Vol. 31 No. 5, pp. 262-280.

Hong, X., Lejeune, M. and Noyan, N. (2015), “Stochastic network design for disaster preparedness”, IIE Transactions, Vol. 47 No. 4, pp. 329-357.

Humanitarian Data Exchange Database (2016), available at: https://data.humdata.org/ (accessed March 2016).

Iakovou, E., Chi, M.I., Douligeris, C. and Korde, A. (1996), “Optimal location and capacity of emergency cleanup equipment for oil spill response”, European Journal of Operational Research, Vol. 96 No. 1, pp. 72-80.

International Federation of Red Cross and Red Crescent Societies (2008), “Guidelines for assessment in emergencies”, technical report.

Iserson, K.V. and Moskop, J.C. (2007), “Triage in medicine, part i: concept, history and types”, Annals of Emergency Medicine, Vol. 49 No. 3, pp. 275-281.

Jia, H., Ordonez, F. and Dessouky, M. (2007), “A modeling framework for facility location of medical services for large-scale emergencies”, IIE Transactions, Vol. 39 No. 1, pp. 41-55.

Larson, R.C., Metzger, M.D. and Cahn, M.F. (2006), “Responding to emergencies: lessons learned and the need for analysis”, Interfaces, Vol. 36 No. 6, pp. 486-501.

Li, J. (2004), “Performance analysis of production systems with rework loops”, IIE Transactions, Vol. 36 No. 8, pp. 755-765.

Li, R. (2010), “Research on emergency supplies demand forecasting model”, 2nd International Conference on Industrial Mechatronics and Automation, Vol. 1, pp. 244-247.

Maggio, N., Matta, A., Gershwin, S.B. and Tolio, T. (2009), “A decomposition approximation for three-machine closed-loop production systems with unreliable machines, finite buffers and a fixed population”, IIE Transactions, Vol. 41 No. 6, pp. 562-574.

Manoj, V., Subodha, K. and Sushil, G. (2016), “An integrated logistic model for predictable disasters”, Production and Operations Management, Vol. 25 No. 5, pp. 791-811.

Natarajan, K.V. and Swaminathan, J.M. (2014), “Inventory management in humanitarian operations: impact of amount, schedule and uncertainty in funding”, Manufacturing and Service Operations Management, Vol. 16 No. 4, pp. 595-603.

Ozbay, K. and Ozguven, E.E. (2008), “Stochastic humanitarian inventory control model for disaster planning”, Transportation Research Record, Vol. 1517 No. 2022, pp. 63-75.

Ozen, M. and Krishnamurthy, A. (2013), “Models to evaluate layout and staffing policies at relief centers”, in Krishnamurthy, A. and Chan, W.K.V. (Eds), Proceedings of the 2013 Industrial and Systems Engineering Research Conference, San Juan, pp. 1235-1244.

Ozen, M. and Krishnamurthy, A. (2017), “In-kind donations for disaster response: Interviews with practitioners”, technical report, University of Wisconsin Madison, Madison, WI.

Roy, D., Krishnamurthy, A. and Bhat, S. (2011), Humanitarian and Relief Logistics: Research Issues, Case Studies and Future Trends, Springer Series Operations Research/Computer Science (ORCS), New York, NY.

Sheu, J.B. (2010), “Dynamic relief-demand management for emergency logistics operations under large-scale disasters”, Transportation Research Part E: Logistics and Transportation Review, Vol. 46 No. 1, pp. 1-17.

Simpson, N.C. and Hancock, P.G. (2009), “Responding to emergencies: Lessons learned and the need for analysis”, Journal of the Operational Research Society, Vol. 60 No. 1, pp. 126-139.

Sodhi, M.S. and Tang, C.S. (2014), “Buttressing supply chains against floods in Asia for humanitarian relief and economic recovery”, Production and Operations Management, Vol. 23 No. 6, pp. 938-950.

Taskin, S. and Lodree, E.J. (2011), “A Bayesian decision model with hurricane forecast updates for emergency supplies inventory management”, Journal of the Operational Research Society, Vol. 62 No. 6, pp. 1098-1108.

Tregenza, P. (1976), The Design of Interior Circulation, Van Nostrand Reinhold, New York, NY.

United States Agency for International Development (2005), “Field operations guide for disaster assessment and response”, technical report, United States Agency for International Development, Washington, DC.

Wang, X., Li, F., Liang, L., Huang, Z. and Ashley, A. (2015), “Pre-purchasing with option contract and coordination in a relief supply chain”, International Journal of Production Economics, Vol. 167, pp. 170-176.

Whitt, W. (1983), “The queuing network analyzer”, The Bell System Technical Journal, Vol. 62 No. 9, pp. 2779-2815.

Whybark, C.D. (2007), “Issues in managing disaster relief inventories”, International Journal of Production Economics, Vol. 108 No. 1, pp. 228-235.

Wu, S., Ru, Y. and Li, H. (2010), “A study on inventory management method in emergency logistics based on natural disasters”, 2010 International Conference on EProduct E-Service and E-Entertainment, ICEEE2010, November 7-November 9.

Xu, X., Qi, Y. and Hua, Z. (2010), “Forecasting demand of commodities after natural disasters”, Expert Systems with Applications, Vol. 37 No. 6, pp. 4313-4317.

Zhao, J. and Cao, C. (2015), “Review of relief demand forecasting problem in emergency logistic system”, Journal of Service Science and Management, Vol. 8 No. 1, pp. 92-98.

Further reading

Cho, S.-H., Jang, H., Lee, T. and Turner, J. (2014), “Simultaneous location of trauma centers and helicopters for emergency medical service planning”, Operations Research, Vol. 62 No. 4, pp. 751-771.

Corresponding author

Merve Ozen can be contacted at: mozen@wisc.edu