Shuva GHOSH,Abhijit GOSAVI

Department of Engineering Management and Systems Engineering,Missouri University of Science and Technology,Rolla,MO 65409,U.S.A.

A semi-Markov model for post-earthquake emergency response in a smart city

Shuva GHOSH,Abhijit GOSAVI?

Department of Engineering Management and Systems Engineering,Missouri University of Science and Technology,Rolla,MO 65409,U.S.A.

An earthquake significant on the Richter scale occurring in an area with a high population density requires an effective and equitable emergency response plan.Emergency resources are usually located in so-called responding centers.One of the first problems faced by disaster-response management personnel in the rapidly degrading post-earthquake conditions is to gage the hazard rate to which the disaster-affected area is subjected,estimate the time taken to bring the situation under control,also called restoration time,and select the appropriate responding center for relief-and-rescue activities.In this paper,we propose an elaborate semi-Markov model to capture the stochastic dynamics of the events that follow an earthquake,which will be used to quantify the hazard rate to which people are exposed and estimate the restoration time.The model will be further used,via dynamic programming,to determine the appropriate responding center.Our proposed model can be employed in conjunction with a variety of hazard scales and by collecting data on a few parameters related to emergency management.The model will be particularly useful in a smart city,where historic data on events following an earthquake would be systematically and accurately recorded.

Emergency management,earthquake,hazard,Markov chains,smart city,degradation

1 Introduction

Earthquakes,hurricanes,tornadoes,volcanic eruptions,fires,floods,snow blizzards,and droughts are examples of natural disasters that can cause deaths and property damage.While each natural disaster has its own characteristic features,earthquakes all over the world share some common characteristics,such as gas leakage,fire,building collapse,and in the worst case floods.When an earthquake considered significant on the Richter scale strikes,disaster-response personnel send resources,e.g.,rescue and medical teams,ambulances,and medicines,to the disaster affected zone.These actions are aimed at evacuating people from the affected area,preventing further damage,and ensuring that people are not hurt in anysecondary disasters that may follow the primary one.Examples of recent natural disasters include the 1999 earthquake in Izmit,Turkey,the 2008 earthquake in Wenchuan,China,the earthquake in Haitiin 2010,the earthquake-tsunami in Japan in 2011,and the recent earthquake in Nepal in 2015.The magnitude of devastation caused by recent earthquakes has intensified the need forsystematicdisaster response planning.

Since resources are limited and must be deployed urgently,disaster response managers have to formulate sensible strategies/policies for allocating resources to affected areas.Response planners,hence,need tools that can help them understand the different eventualities that they might encounter.Operations research/optimization has been playing a key role in developing tools to help understand disaster management(see e.g.,[1]and[2]for a comprehensive review of operations research for disaster management)and optimize the decision-processes associated to responding to disasters(see e.g.,[3]).

In recent times,urban and rural living has been becoming increasingly “smart”.Smartcities are already a reality[4],and recent literature strongly suggests that smart cities are making intelligent use of data gathered from sensors to make them more liveable and sustainable[5–7]and more risk-free[8].Cyber-physical systems are increasingly being viewed as examples of smart cities where the cyber-components can enable the connections[9].In a world of big data,Internet of Things(IoT),and smart living,the need for automated response to emergencies is being felt strongly:see in particular[10]and[11].In these smart/connected cities,numerous data streams are collected[4],which can be used to model the stochastic dynamics of the events following an earthquake,as well as the performance characteristics of responding centers,e.g.,how long it takes to arrive at the scene and how long it takes to perform a restoration activity like dousing a fire etc.

An important question studied here is:Where should relief and rescue teams come from?i.e.,alocalresponding center or to thecentral responding center.The local responding center is likely to arrive on the scene quickly,but is typically slower,while the opposite would be true of the central responding center located further away but one that has additional resources.There are situations where both the responding centers are from a neighboring area with differing capabilities,and the decision has to be made regarding which center to choose.To be more specific,the model we propose should be able to i)select the appropriate responding center,depending on the situation at hand,and ii)answer what-if type of questions that can be used to compare and contrast different response strategies.

The events that follow an earthquake can have varying levels of the hazard/risk that they pose to the individuals in the disaster zone.For instance,the four major post-earthquake events that we consider in this paper,namely,gas/chemical materialleakage,fire,building/structure collapse,and flooding,do not pose the same level of danger.Thus,a gas leakage is usually not as dangerous as a fire;a building collapse is likely to be more hazardous than a fire or a gas leakage,while catastrophic floods are possibly the worst in this hazard spectrum,where helicopters may eventually be needed for rescue.Indeed,the actual conditions that follow an earthquake are often a combination of many of these events.Therefore,the decision-making processes involved must distinguish between these events in a manner such that the resulting rescue strategy is the most effective possible and the response is as well-coordinated as it can be.

A well-coordinated response would also require for instance the automatic shutting down of gas lines and nuclear reactors and rasing readiness level of fire and ambulance services,apart from announcements on radio/TV,stopping trains,and cars pulling over to make way for rescue teams.While most of these measures can be implemented as an immediate response to an earthquake,clearly,the model we propose would require assessment of the current state of the system and then data processing,i.e.,solving a semi-Markov model(which can be done instantaneously in modern computers).Fig.1 explains how the semi-Markov model would be implemented in a smart environment to generate an automatic response that would notify the appropriate response center.

Fig.1 Automated response mechanism in a smart city.

The three key modeling issues related to comparing different strategies of resource allocation after an earthquake are:i)How does one quantify the level of danger(developing a hazard scale)posed by the events(e.g.,building collapse,fires etc)?ii)How does one estimate performance characteristics(variables)related to the smart city?iii)How does one model the dynamics of the events that follow the earthquake?Fig.2 depicts the modeling process that would be needed.A goal of this paper is to address these questions by mathematically capturing the stochastic dynamics of the events via the semi-Markov model,defining the input data needed for the model,and using a flexible scale for measuring the hazards/dangers associated to each event(gas leakage/fire/building collapse/flooding).The model will specify the optimal responding center and will also be able to answer what-if type questions related to estimating the time it takes to bring the situation under control and the mean hazard exposure under any given strategy.

Fig.2 Modeling response to the earthquake in a smart city.

Contributions of our paperThe overarching objective of this research is to develop a stochastic,control-theoretic(semi-Markov),discrete-event model that studies the degradation of a system subjected to an earthquake and how the model can be used to help the system survive and be restored to its original functional state in an effective manner.The specific contributions of our paper can be summarized as follows:i)To the best of our knowledge,our work is the first to develop a comprehensive solution methodology that quantifies the level of hazard(danger)via a hazard scale in an earthquake-response model.ii)Our model also captures the stochastic dynamics of four major earthquake hazards,namely,fire,gas leakage,structural collapse,and floods,that should help in decision-making of smart city disaster managers—in the three or four days immediately following a disaster.iii)Finally,our model develops equations that estimate restoration times and the average hazard exposure of a given strategy,thereby allowing managers to conduct what-if analyses of different strategies;this is critical for disaster managers in better understanding different strategies and their impact and thus plan for response resources and their locations.

The rest of this paper is organized as follows.Section 2 provides a literature review.In Section 3,we present details of the problem formulation.The algorithm and the resulting what-if analysis are detailed in Section 4.Numerical results are discussed in Section 5.Section 6 concludes this paper with some remarks on this research and suggestions for future research.

2 Literature review

While the literature on disaster management is extensive,we focus on the literature on conditions encountered immediately after earthquakes and the response strategies needed.Swersey[12]was one of the first to develop a Markov chain model for emergency resources deployment needed in responding to fire signals received by the local fire fighting services.Fiedrich et al.[13]is one of the few works that directly employ probabilities of occurrence of incidents,such as dam failures and fire,in their model.They develop a comprehensive model for resource allocation to reduce the number of casualties following an earthquake.They also state that the three days after an earthquake are crucial for performing relief activities,which necessitates a model that captures the dynamics of the events that unfold in the short time following an earthquake.Also of interest to us here is the domino effect that occurs when one event leads to another.Khan and Abbasi[14]have identified fire,explosions,toxic release,etc.,as events that can initiate the domino-effect and their probabilities.The book of Sekijawa et al.[15]shows how to estimate the probability of fire after an earthquake.Other related works include Chang[16],who presents an urban transportation planning model that computes a per-formance measure for accessibility of an area after an earthquake,which can be used to study loss of transportation services,as well as post-disaster detouring and congestion.The book by Kreimer at al.[8]studies how cities can be made safer by studying degradation and survivability.Zhu et al.[17]develop a resource allocation model for local resource reserve depots after a disaster.

The work closest to ours is that of Wei et al.[18],which considers a two-state Markov chain model that accounts for two types of incidents following the earthquake in Wenchuan,China:a fire and a gas leakage.Thus,with the exception of[18]that accounts for fire and gas leakage,what the literature lacks is a stochastic model that also accounts for theothertwo major disaster-inducing events that typically occur after the earthquake,i.e.,collapsing of buildings and floods,as well as the combinations of the occurrence of the four disaster-inducing events(fire,gas leakage,building collapse,and floods).Our paper seeks to fill this gap in the literature.

3 Problem formulation

This section formulates the problem under study and discusses the development of the hazard scale needed for the solution.Section 3.1 provides details of the different incidents that can follow an earthquake,while the hazard scale is described in Section 3.2.The semi-Markov model is presented in Section 3.3.

3.1 Earthquake dynamics and control

In this section,we focus on the events that follow an earthquake and how its effects can be controlled in order to set the stage for the subsequent model.In order to develop a comprehensive model for an earthquake,it is necessary to take into account most,if not all,of the incidents that are likely to occur following the earthquake[13,19].The federal emergency management agency(FEMA)of the U.S.in its earthquake safety guide for homeowners[20]lists gas and chemical material leakage,fire,buildings/bridges/structures collapsing,floods(dam breakage),etc.,as some of these incidents.In particular,gas leakages are known to lead to fires and structural collapses,e.g.,dam breakages,to floods.The report mentions that 25%of the fires after an earthquake are due to gas leakage(see Page 25 of[20]).Our model will account for four key incidents:i)gas and chemical material leakage,ii)fire,iii)building(including structures such as dams)collapse,and iv)floods.Typically,the situation following an earthquake can worsen quickly–with one incident leading to another with similar or worse consequences.This phenomenon of the primary incident,leading to another incident,namely the secondary incident,with a“magnitude equal or higher than that of the primary one”is often called the domino-effect[21,22].The probability of one incident leading to another differs from area to area,depending on the area’s population,the height of buildings,and the number of earthquake-prone structures,amongst other factors.These probabilities can be obtained from subject matter experts who are usually civil engineers working for disaster management authorities.

Immediately after a significant earthquake strikes,the authorities in charge respond by dispatching teams of ambulances,medical personnel,fire fighters,rescue personnel,forklifts,excavators,towing trucks,concrete cutters,etc.,to the affected location.When the earthquake strikes at an area,the decision-maker has two choices:call the nearest responding agency(local agency)or the central agency.Consider Fig.3.If the emergency response is needed at CityXor CityY,one could reach out to either the local agency,which is nearer,or the central agency,which is further away.Typically,emergency personnel from the nearest agency can arrive sooner,but they may have fewer resources.On the other hand,the central agency is likely to take longer to arrive on the scene,but once it does arrive,it can bring the situation under control in less time,due to the fact that it has access to a greater volume of resources.

Fig.3 A map of two cities,X and Y,served by a local agency and a central agency.

3.2 The hazard scale

A scale that quantifies the level of risk/hazard posed by a given incident will be necessary for our resource allocation model and what-if analyses of different strategies.Hazard scales have been used in other contexts,e.g.,to measure risks associated to terrorism,where there are four levels:green,yellow,orange and red[23].In the context of incidents that occur after a disaster,Lin and Jarrett[24]have analyzed quantitative risk assessment methods.See also Antonioni et al.[25]for a methodology to perform quantitative risk assessment of major accidents occurring in an industrial complex due to earthquake impacts.In general,the hazard exposure in an area associated to an event is dependent on the nature of the incident as well as the population density of the area.We now describe a flexible numerical scale for the hazard-inducing incidents in an earthquake.

Naturally,the level of danger posed by each incident tends to be different.Thus,for instance,gas and chemical material leakage alone do not pose the kind of risk created by floods.A quantitative value is hence needed to be assigned to each incident–in order to differentiate between the incidents.Clearly,the assignment of values to each level of hazard will be up to the analyst under the premise that higher the value,the greater the level of danger.Furthermore,it will also depend on the population density of the area under consideration and the susceptibility of the population to the hazard under consideration.Thus,for instance,the St.Louis region in the state of Missouri,U.S.A.is in a region where a significant earthquake is expected to occur.It is an urban region that is densely populated,and is on the banks of a major river,Mississippi.Thus,the St.Louis region is likely to have a higher hazard susceptibility than a rural region in Arizona,U.S.A.with a significantly lower population density.

We provide two examples for hazard scales in Table 1,and now explain this concept.The value of the hazard,assigned in the hazard scale,will be a metric for the level of danger posed to an area in unit time.In the first example,scale I,the incident of gas and chemical material leakage is assigned the minimum level of hazard of 2,while the incident of floods is assigned the maximum value of 24;the incidents of fire and building collapse will have intermediate values of 22and 23respectively.Thus,for example,in this scale,spending unit time in an area under floods would be 8 times more dangerous than spending the same amount of time in an area where the only incident is gas and chemical material leakage.It is also to be understood that if an area has multiple incidents occurring simultaneously,the total rate of hazard will be additive in the incidents.Thus,scale I would assign a hazard value of(4+8)=12 to an area where fires and floods are raging.As we will see later,multiple combinations are possible,and our overall goal is to minimize the time rate of hazard posed to people affected.Anotherscale,named scale II,is shown in Table 1.In this scale,the values of hazard also increase with the level of danger posed,but are not powers of 2;this is provided as another example to show that the scale is flexible and the values can be assigned by the analyst.We also note that our model is general enough to handleanygiven scale provided for the incidents,i.e.,any reasonable scale that accounts for the nature of the incidents and the susceptibility and density of the population exposed to the incident.In general,our model will designed to minimize the rate of hazard,and in particular,its objective function can be expressed as

where E is the expectation operator.

Table 1 Hazard scales.

3.3 Semi-Markov decision process

In this section,we describe the semi-Markov decision process(SMDP)model that will be used in this paper.We provide details of the state and action space,as well as the transition probability,cost,and time models.However first,we formally present SMDP-related notation that will be needed in the re-mainder of the paper.

3.3.1 Notation

.S:Set of decision-making states.

.S′:Set of natural states.

.A:Set of actions in any decision-making state.

.μ(i):The action chosen in stateiwhen policy μ is followed.

.TCM(·,·,·):S × A × S → R:One-step immediate cost incurred in going from a decision-making state to another under a given action.

.TTM(·,·,·):S × A × S → R:The time spent in one transition from one decision-making state to another under a given action.

.TPM(·,·,·):S × A × S → [0,1]:The associated transition probability of going from one decision making state to another under a given action.

.ρμ:The average cost per unit time of the policy μ,which is also our objective function in the model.

3.3.2 State and action space

As stated above,the state will be defined by the combination of incidents,defined in Table 1,that can occur simultaneously.The sixteen resulting states in our model are explicitly defined in Table 2.Underlying an SMDP,one has the so-called decision-making process and the natural process[26].The natural process accounts for all the states in the system,including those in which no decision-making is performed,where S ? S′.The states numbered 1 through 8 in Table 2 are the decision-making states,while the remaining states in the table are of the non-decision-making type(also called natural states).

As stated earlier,the selection of a particular responding center denotes an action in our SMDP model.There are two possible actions:Action 1 will select the local responding center,and action 2 will select the central responding center.Thus,A={1,2}is the action space.

Table 2 States:The states numbered 1 through 8 are decision-making states and the remaining states are of the non-decision making type(natural states).

3.3.3 Transition matrices

We now develop the transition probability,cost,and time matrices needed for the SMDP model.Three of the key elements of an SMDP model are the three types of matrices,one for each action:the transition probability matrix(TPM),the transition time matrix(TTM),and the transition cost matrix(TCM).These matrices essentially define the dynamics of the system as the states evolve.The evolution of the states in our model can be explained as follows.See Fig.4.The system initiates in the fully functional state,also called stable state,which is numbered 1.The earthquake occurs,and the system transitions immediately to any one of the states numbered 2 through 8,denoted byiin our notation and in Fig.4.At this point,a decision is made,i.e.,the responding center is selected.The system then transitions to a non-decision-making(or natural)state,which denoted byi′in our notation and Fig.4.This transition occurs during the time interval in which the responders are traveling from the response center to the affected location.Now,after the responders arrive and the response is complete,the system transitions from the non-decision-making state to the fully functional(or stable)state with probability 1.

Fig.4 A schematic showing the evolution of the decision–making and natural states.

LetTPM(i,a,j)denote the transition probability that the system goes from a decision-making stateito a decision-making statejin one transition under actiona.Then,for both actions,the TPM has the following general structure for our model:

where a?indicates a non-zero value.Note that,as discussed above,the system will transition from every state,except for state 1,to state 1 with probability1.It also needs to be noted that our TPM will only account for the eight decision-making states,which are the states numbered from 1 through 8.

The probabilities of going from a decision-making state to a natural state will be referred to as domino probabilities.The domino probability of transiting from a decision-making stateito a natural statei′when actionawas pursued in stateiwill be denoted byD(i,a,i′).The probabilities in the TPM and theDmatrices can be obtained from data collection that can be done in a smart city.

The transition time can be defined as the time taken by the systemingoing from one decision-making state to another in one transition.The transition time will depend on the travel times and the response times.TTM(i,a,j)will define the expected transition time taken by the system to go from stateito statejunder actionain one transition.In order to define the transition time matrix,we will need additional notation:We denote the mean travel time for an actiona,i.e.,the time taken by the responders to arrive at the scene after the responding agency is notified,by TT(a).Also,the response time will be defined as the mean time taken by a responding center to bring the situation under control after the responders from that center arrive on the scene.Clearly,the response time will depend on the state of the system in which the center finds the system when it arrives and the responding agency.We will use RT(i′,a)to represent the mean response time under actionawhen the system is in statei′,one of the natural states.The transition times in the transition time matrix can then be defined as

When statei∈{2,3,...,8},

For all other cases,the transition times will equal zero,i.e.,TTM(·,·,·)=0.

Equation(2)can be explained as follows for our model.Note that the stable state(1)is assumed to instantaneously transfer to any one of the states numbered 2 through 8 and hence the transition time from 1 is 0.The transition time from any one of the states numbered 2 through 8 is a sum of the travel time(under the action chosen) and the expected response time from the natural state,to which it transitions(i′),to the stable state.Clearly,since the system can transition from stateito any one of the states in the natural process,one must compute an expectation over all natural states,i.e.,over S′.

To define the transition cost matrix,we also need some additional notation.LetH(φ)denote the hazard level of incident φ,where φ ∈ I.Hence,using our hazard scale I in Table1:

When statei∈{2,3,...,8}:

For all other cases,TCM(·,·,·)=0.

Equation(3)follows in a manner similar to that explained for equation(2),and hence we relegate the explanation of equation(3)to Appendix A.Also,how the transition and domino probabilities are estimated is also explained in Appendix B for the reader’s reference.

4 Solution technique

In this section,we first describe the policy iteration algorithm used to optimally solve an average cost SMDP via dynamic programming[26,27],and then discuss key what-if analyses that can be performed in this context by the disaster manager to evaluate the efficacy of different strategies in terms of metrics such as restoration time and the average hazard exposure.

4.1 Algorithm steps

The SMDP is aimed at finding the best policy μ?,i.e.,one that minimizes the average cost,ρμ.The optimal average cost will be denoted by ρ?in what follows.

Step 1Letkdenote the iteration number,which is initialized to 1. Select any policy in an arbitrary manner and denote it by μk.Let μ?denote the optimal policy.

Step 2(Policy evaluation) Solve the following linear system of equations.For alli∈S:In the above linear system of equations,the unknowns are thehkterms and the term ρk.Any one of thehkterms should be set to 0 in the above to obtain a unique solution.

Step3(Policy improvement)Choose a new policy μk+1so that for alli∈ S

If possible,one should set μk+1= μk.

Step 4If the new policy is identical to the old one,i.e.,if μk+1(i)= μk(i)for eachi,then stop,and set μ?(i)= μk(i)for everyi.Otherwise,incrementkby 1,and return to Step 2.

4.2 What-if analysis

Finally,in this section,we will present formulations that can be used for performance evaluation and for answering what-if questions, which typically do not require optimization. Two important issues that disaster response managers would be interested in here are i)the time it takes to bring the situation under control and ii)the expected hazard rate to which the area would be subjected.It would be of interest here to compute these metrics for a given center,either the local or the central one.As such,we now present the formulas needed for evaluating both of these metrics.

Restoration timeWe define restoration time as the amount of time needed to bring the situation under control after an earthquake,i.e.,all the affected people are taken to hospitals,in addition to water,electricity,and road connections getting restored.The recovery time is a different concept,and is usually defined as the time it takes to rebuild all the structures that have collapsed,which can take many months or even years.Studying recovery time is beyond the scope of this paper.For the model that we have developed,the restoration time is the time elapsed since the earthquake strikes until the situation is brought under control.Since,the system goes from the stable state to a starting state(depicted byiin Fig.4),the restoration time will depend on the starting state as well as the action chosen in the starting state.Let τa(i)denote the restoration time from starting stateiwhen actionais chosen.Then,

Note that in the above,a=1 when the local center is used,whilea=2 when the central center is used.

Average hazard exposureFor a given policy μ,the average hazard exposure is defined by ρμ.To compute the value of ρμ,we first compute the steady-state probabilities of the Markov chain underlying policy μ by solving the following system of linear equations:

where Πμ(i)equals the steady-state probability of stateiassociated with policy μ and denotes the unknown in the equation above.Then,we have that

We will illustrate the use of the formulas in equations(5)and(6)in the following section on numerical results.

5 Numerical results

In this section,we illustrate our model via numerical examples.The first example(test case)we study is named Case 1 and will serve as our baseline case.This will be used as a vehicle to provide all the details for calculations,so readers can duplicate them for their own research.Thereafter,other test cases will be described.We first describe the inputs needed in our experiments.

5.1 Inputs

For using our model,we need the following inputs:theTPMand theDmatrices,as well the travel times and response times that we describe later;the values of all of these can be obtained as discussed above.We will use the following matrix as the TPM for Case 1;note that this matrix is the same for both actions.

The other key input for the model are the values of theDmatrices.Dawill denote the matrix for actiona,where the term in theith row andi′th column ofDawill essentially be equivalent toD(i,a,i′).Note thatiwill take values from SDwhilei′will take values from S.Thus,for our case,we will have an 8×16 matrix.For Case 1,we will use the following two matrices:

For Case 1,the travel times used in our calculations will be as follows:TT(1)=1/2 hour and TT(2)=3/2 hours,where TT(a)denotes the travel time for actiona.We now present the data for the response times.Let Ta(φ)denote the response time when incident φ is involved and actionais chosen.The response times for action 1 in hours that we used for Case 1 are as follows:

The response times for action 2 in Case 1 in hours are as follows:

Finally,the hazard levels,TH(·),under the two scales,were defined in Table2.

For test Cases numbered 2 through 8,we experimented with a different set of response times to see how the algorithm behaved;all other data remain identical to those of Case 1.Cases 2 and 6:T2(FL)=150;Cases 3 and 7:T2(FL)=155;Cases 4 and 8:T2(FL)=200.Also,Cases 1 through 4 use hazard scale I,while the remaining cases use scale II.

Cases 9 and 10 model a scenario based from the state of Missouri.The two cities for which this model could be applied are Rolla and Columbia within the state of Missouri.See Fig.5,where the local agency for Rolla(or Columbia)is in Jefferson City and the central agency is in St.Louis.We use data for the city of Rolla in our experiments.The distance between Jefferson City and Rolla can be traversed in approximately 1 hour(i.e.,TT(1)=1hr),while the same between

Rolla and St.Louis can be traversed in about 2.5 hours(i.e.,TT(2)=2.5hrs).Also for Cases 9 and 10,we use the following data for the response times(in hours)based on local estimates:

The response times for action 2 in Cases 9 and 10 are as follows:

All other inputs for Cases 9 and 10 are as defined in Case 1.The value of α in equation(3)was set to 0.5 in all our experiments(Cases 1 through 10);this value can be changed by the analyst depending on historical data.

Fig.5 Map of a part of the state of Missouri,U.S.A.,where for the city of Rolla(and Columbia),Jefferson City serves as the local agency and the city of St.Louis serves as the central agency.

5.2 Outputs

The inputs have to be now pre-processed for some other data that will be needed by our algorithm.We first describe how to obtain the response times for each state-action pair,RT(·,·).The response time for a statei′under actionawill be the sum of the response times under actionafor all the incidents involved in statei′.Mathematically,

where I(i′)denotes the set of incidents associated to statei′,the natural state.Thus,for example in Case 1,the response time for state 8 under action 1,RT(8,1),which involves,gas leakage,fire,and building collapse,will be 12+36+60=108 hours while the same under action 2 will be 12+24+48=84 hours.Similarly,the response time for state 16,which involves all four incidents(gas leakage,fire,building collapse and floods),will be 12+36+60+200=308 under action 1 and 12+24+48+100=184 under action 2.

Using the inputs described above and the values of the matrix RT(·,·),we computed the TCM(·,·)and the TTM(·,·)matrix via equations(3)and(2)respectively,in a computer program and ran our algorithm in the same program.The policy iteration algorithm was also run on each case to determine the optimal solution;the algorithm converged to the optimal solution in each case in less than 2 seconds on an Intel Pentium Processor with a speed of 2.66GHz on a 64-bit operating system.The software used for the computer programs was MATLAB.Results of running our algorithm on all the cases are provided in Table3;note that in the table,a policy μ is denoted by〈μ(1),μ(2),...,μ(|S|)〉,in order to specify the action chosen by the policy in each decision-making state.Thus,for example,for Case 1,for states 1 through 3,the optimal action is 1(i.e.,choose the local agency),while for states 4 through 8,the optimal action is 2(i.e.,choose the central agency).

Finally,in Tables 4 through 6,we demonstrate the answers of some what-if questions.Tables 4 and 5 show the values of the restoration time(defined in equation(5))for sample cases.In particular,the interpretation of the results in these tables is that given the initial unstable state,the manager can estimate the amount of time(in hours)needed to bring the situation under control(restoration time)starting from when the agency concerned is notified.Table 6 shows the values of the average hazard rate for sample cases,should some policy other than the optimal be chosen.The value of the average hazard rate can be computed via equation(6).Analysis of this nature can be of great help to disaster managers for resource capacity planning,as well as location planning of resources.Further,analysis of this nature also helps in superior preparation for a disaster and helps in better coordination of restoration activities when a disaster strikes.

Table 3 Results on 10 test cases where ρ?denotes the average cost of the optimal policy.

Table 4 Restoration times(in hours)on sample cases.

Table 5 Restoration times(in hours)on cases from Missouri.

Table 6 What-if analysis of sample policies.

6 Conclusions and future research

The problem of coordination of activities that are required in the immediate aftermath of an earthquake has been a topic of interest for many years now.However,there is no mathematical model in the literature that can capture thestochastic dynamics of the events that occur immediately after an earthquake,e.g.,gas and chemical material leakages,fire,structural collapse,and floods.More generally,it is necessary to develop a model that studies how a complex system degrades after an earthquake and how disaster response can be carried out to increase the system’s survivability.This paper seeks to fill this gap in the literature and presents a comprehensive model that not only captures the dynamics of a post-earthquake scenario but also delivers a strategy for selecting the responding center in an optimal manner.In particular,the model is geared towards development of an automatic response system that will be needed in smart and connected cities.Further,the model will also be effective as part of a pre-disaster response planning exercise in any city.The modeling approach was geared towards minimizing the mean hazard rate to which a disaster-affected zone is subjected in a post-earthquake scenario.A generalized framework was developed, which was flexible enough for the user to change inputs depending on ground realities.

A number of directions for future research can be explored.First,large-scale simulation-based models for an entire state(a large geographical area)in a country,which could exploit the domino-effect information with less analytical effort than is needed here,should be very attractive in practice.Second,developing a transportation plan based on the optimal results from this model should be of interest to practitioners.

[1]H.O.Mete,Z.B.Zabinsky.Stochastic optimization of medical supply location and distribution in disaster management.Internation Journal of Production Economics,2010,126(1):76–84.

[2]N.Altay,W.G.Green III.OR/MS research in disaster operations management.European Journal of Operational Research,2006,175(1):475–493.

[3]O.L.Huibregtse,S.P.Hoogendoorn,A.Hegyi,et al.A method to optimize evacuation instructions.OR Spectrum,2011,33(3):595–627.

[4]B.Marr.How big data and the internet of things create smarter cities.Forbes,2015:http://www.forbes.com/sites/bernardmarr/2015/05/19/how-big-data-and-the-internet-ofthings-create-smarter-cities/#1c98a7e663d8.

[5]M.Batty.Big data,smart cities and city planning.Dialogues in Human Geography,2013,3(3):274–279.

[6]R.Kitchin.The real-time city?Big data and smart urbanism.Geo Journal,2014,79(1):1–14.

[7]A.M.Townsend.Smart Cities:Big Data,Civic Hackers,the Quest for a New Utopia.New York:W.W.Norton&Company,2013.

[8]A.Kreimer,M.Arnold,A.Carlin.Building Safer Cities:The Future of Disaster Risk.Washington,D.C.:World Bank Publications,2003.

[9]L.Gurgen,O.Gunalp,Y.Benazzouz,et al.Self-aware cyberphysical systems and applications in smart buildings and cities.Design,Automation&Test in Europe Conference&Exhibition,Grenoble,France:EDA Consortium,2013:1149–1154.

[10]S.Menosono.In the Wake of Japan’s Earthquake and Tsunami.2015:http://www.hitachiconsulting.co.jp/hc_images/news/pdf/news_120628_en.pdf.

[11]E.Asimakopoulou,N.Bessis.Buildings and crowds:Forming smartcities for more effective disaster management.Proceedings of the 50th International Conference on Innovative Mobile and Internet Services in Ubiquitous Computing,Seoul:IEEE Computer Society,2011:229–234.

[12]A.J.Swersey.A Markovian decision model for deciding how many fire companies to dispatch.Management Science,1982,28(4):352–365.

[13]F.Fiedrich,F.Gehbauer,U.Rickers.Optimized resource allocation for emergency response after earthquake disasters.Safety Science,35:41–57,2000.

[14]F.I.Khan,S.A.Abbasi.Models for domino effect analysis in chemical process industries.Process Safety Progress,1998,17(2):107–123.

[15]A.Sekijawa,M.Ebihara,H.Notake.Development of seismic-induced fire risk assessment method for a building.Proceedings of the 5th International Conference on Probabilistic SafetyAssessment andManagement,Tokyo:Universal Academy Press,2003:1927–1932.

[16]S.E.Chang.Transportation planning for disasters:An accessibility approach.Environment and Planning A,2003,35(6):1051–1072.

[17]J.Zhu,J.Huang,D.Liu,et al.Resources allocation problem for local reserve depots in disaster management based on scenario analysis.Proceedings of the 7th International Symposium on Operations Research and Its Applications,Lijiang:World Publishing Corporation,2008:395–407.

[18]W.Wei,L.Mao,W.Li.Dynamic optimization methods of emergency resources deployment based on Markov decision process for Wenchuan earthquake.Proceedings of the 2nd International Workshop on Database Technology and Applications,Wuhan:IEEE,2010.

[19]C.Scawthorn,J.M.Eidinger,A.J.Schiff.Fire Following Earthquake.Reston,U.S.A.:American Society of Civil Engineers(ASCE),2005.

[20]FEMA.Earthquake Safety Guide for Homeowners.2005:http://www.sjsu.edu/ihouse/docs/fema-530.pdf.

[21]C.Delvosalle.Domino effect phenomena:Definition,overview,classification.Proceedings of the 1st European Seminar on Domino Effects,Leuven,Belgium:Federal Ministry of Employment,Safety Administration,Direction Chemical Risks,1996:5–15.

[22]F.Clini,R.M.Darbra,J.Casal.Historical analysis of accidents involving domino effect.Proceedings of the 4th International Conference on Safety&Environment in Process Industry,Florence,Italy:Aidic Servizi Srl,2010:335–340.

[23]L.A.Cox.What’s wrong with risk matrices?Risk Analysis,2008,28(2):497–512.

[24]X.Lin,R.Jarrett.A practical approach to quantitative risk assessment.Proceedings of the Risk Conference.Wellington,Newzealand:Society for Risk Analysis,2009.

[25]G.Antonioni,G.Spadoni,V.Cozzani.A methodology for the quantitative risk assessment of major accidents triggered by seismicevents.Journal of Hazardous Materials,2007,147(1/2):48–59.

[26]D.P.Bertsekas.Dynamic Programming and Optimal Control.2nd ed.Belmont:Athena,2000.

[27]A. Gosavi.Simulation-basedOptimization:Parametric Optimization Techniques and Reinforcement Learning.New York:Springer,2003.

Appendix

A1 Details on equation(3)

Equation(3)can be explained as follows.The time it takes to transition fromito state 1 is essentially comprised of two epochs:the travel time and the response time(see Fig.4).Hence the total hazard encountered when going from stateito state 1 is the sum of the hazard encountered during the travel time and the response time.The first term in equation(3)equals the hazard during the travel time,while the second term equals the total hazard encountered in the statei′during the response time.

A2 Estimating transition and domino probabilities

Various elements of the model presented above can be estimated from the smart data in an intelligent city and from data obtained from subject matter experts.The TPMs can be estimated as follows:For every(i,a,j)combination:

whereN(i,a,j)is the number of recorded transitions from decision-making stateito decision-making statejunder actiona.The domino transition probabilities can be estimated as follows:

whereM(i,a,i′)denotes the number of recorded transitions from decision-making stateito natural statei′under actiona.Data of this nature,i.e.,N(·,·,·)andM(·,·,·),can be extracted from past earthquakes.Indeed,it is expected that in a smart city,data of this kind would be recorded and stored after each disaster.The values of the vector TT(·)can be obtained from geographical distances and those of the vector RT(·)can be obtained from subject matter experts.

15 April 2016;revised 5 December 2016;accepted 5 December 2016

DOI 10.1007/s11768-017-6060-y

?Corresponding author.

E-mail:gosavia@mst.edu.Tel.:+573-341-4624.

This work was partially supported by the University of Missouri Research Board(No.09-01-2010).

?2017 South China University of Technology,Academy of Mathematics and Systems Science,CAS,and Springer-Verlag Berlin Heidelberg

Shuva GHOSHis an Assistant Professor at Bangladesh University of Engineering and Technology.He obtained a Ph.D.in Engineering Management from the Missouri University of Science and Technology in 2013.His research interests are in Markov processes,simulation,and optimization.He is a member of INFORMS.E-mail:shuva99@gmail.com.

Abhijit GOSAVIis an Associate Professor at Missouri University of Science and Technology in the Department of Engineering Management and Systems Engineering.He obtained a Ph.D.in Industrial Engineering from the University of South Florida in 1999.His research interests are in simulation-based optimization and Markov decision processes.He is the author of the Springer text,Simulation-based Optimization,currently in its second edition.He is a member of INFORMS,IIE,and ASEM.E-mail:gosavia@mst.edu.