Ghazaleh Ahmadi

Reza Tavakkoli-Moghaddam

Armand Baboli

Mehdi Najafi

A Decision Support Model for Robust Allocation and Routing of Search and Rescue Resources after Earthquake: A Case Study

Keywords: Search and rescue, Earthquakes, Disaster management, Response phase, Robust allocation routing Fatality min. Simulated Annealing No No No No Allocation Deterministic Yes No Travel Rescue costs min. Heuristic

The efficient planning of search and rescue (SAR) operations is highly impactful in the disaster response phase, which offers a limited time window with a declining chance for saving trapped people. The present paper introduces a new robust decision support framework for planning SAR resource deployment in post-disaster districts. A two-stage decomposition approach is applied to formulate the problem as iterative interrelated stages of mixed-integer programming (MIP) models. The first stage presents a robust multi-period allocation model for maximizing fair and effective demand coverage in the affected districts during the entire planning horizon. It takes into account the time-sensitiveness of the operations via a time-dependent demand satisfaction measure and incorporates resource transshipment optimization. The second stage optimizes the routing of the resources allocated in the first stage for each district during the upcoming period. It aims to minimize the weighted sum of SAR demand fulfillment times under consideration of secondary destruction risk, resource collaboration, and rest time requirements. At the end of each period, the proposed framework can be re-executed to capture updated resource, demand, and travel time parameters. To tackle the environment's inherent uncertainty, an interval-based robust optimization approach is adopted. The proposed framework is solved and analyzed for an urban zone in Iran under an earthquake scenario. Results show that the proposed robust models have superior performance compared to a deterministic approach for adaptation to an uncertain disaster environment. More importantly, they prove to be a strong analysis tool for providing helpful managerial insights for the mitigation and preparedness phases.

Introduction

Earthquakes have caused severe damage and numerous fatalities in recent years. They continue to catastrophically threaten millions of people's lives, directly or indirectly, all over the world every year. According to the Significant Earthquake Database 1 (NGDC/WSD), major earthquakes during the last decade caused more than 400,000 deaths and more than 500,000 injuries, and affected 103 countries. The Sichuan earthquake in China in 2008 and earthquakes in Haiti in 2010 and Nepal in 2015 were among the most destructive and affected more than 56 million people (EM-DAT). Immediately after earthquakes strike [START_REF] French | The varied contexts of environmental decision problems and their implications for decision support[END_REF], possible sudden structural collapses trigger search and rescue (SAR) operations in the response phase of the disaster management cycle, which aims at locating, extricating and stabilizing trapped people (OCHA, 2015).

SAR is a fight against time since the survival likelihood of trapped people declines over time [START_REF] Noji | The public health consequences of disasters[END_REF][START_REF] Olson | Urban heavy rescue[END_REF]. Hence, providing rescue operations in a disaster-impacted zone requires the timely availability of sufficient SAR resources based on the damage magnitude. However, given the response phase setting, especially in major disasters that cause damages over a large zone with multiple districts, these resources can face serious shortages for various reasons: failure in the preparedness phase to devise an agile initial response to the disaster; budget limitations; or logistical constraints. Despite the low presence of SAR teams in the vital initial hours, response teams from various local, national, and international governmental or nongovernmental organizations gradually arrive. Therefore, developing well-structured mechanisms to alleviate the intra-and inter-organizational coordination challenges in disaster response is beneficial and important [START_REF] Quarantelli | Disaster crisis management: A summary of research findings[END_REF][START_REF] Altay | OR/MS research in disaster operations management[END_REF][START_REF] Balcik | Coordination in humanitarian relief chains: Practices, challenges and opportunities[END_REF][START_REF] Galindo | Review of recent developments in OR/MS research in disaster operations management[END_REF]. The complexity in the availability of resources, on the one hand, and the diminishing survival chances of trapped people on the other hand, highlight the significance of timely resource allocation in the life-saving endeavors of the response phase [START_REF] Rolland | Decision support for disaster management[END_REF][START_REF] Sawik | An integer programming approach to scheduling in a contaminated area[END_REF][START_REF] Galindo | Review of recent developments in OR/MS research in disaster operations management[END_REF][START_REF] Hoyos | OR models with stochastic components in disaster operations management: A literature survey[END_REF] is a complex and difficult task [START_REF] Clemen | Making hard decisions: an introduction to decision analysis[END_REF]. In such a context, the existence of a decision support framework, which tackles these challenges, can accommodate the decision makers with tools to analyze the situation and to make informed and right choices [START_REF] Thompson | Improving disaster response efforts with decision support systems[END_REF][START_REF] Cioca | Decision support systems used in disaster management[END_REF][START_REF] Othman | An agent-based decision support system for resources' scheduling in emergency supply chains[END_REF].

Inspired by the vulnerability of Tehran in the face of potential earthquakes, and to address the challenging, highly complex and impactful problem of rescue operation management, this study develops a novel robust decision support framework for SAR planning optimization. The proposed framework, which also tackles the inherently uncertain nature of the disaster environment by adopting interval-based robust optimization, integrates two major decisions: (1) allocation of available SAR teams among districts in the affected zone; and (2) routing of allocated resources to serve the demand locations in each district of the zone.

The main contributions of the developed decision support framework can be categorized as follows:

 Considering fairness and operational efficiency in allocation decisions.

 Addressing the risks of disruption in the operations from the occurrence of secondary destruction.  Considering the inherent uncertainty of disasters by adopting an interval robust optimization framework for SAR allocation and routing.  Incorporating two coordination mechanisms to increase management effectiveness in the affected zone. The rest of this paper is structured as follows. Section 2 reviews related studies in the literature. In Section 3, the problem is defined and the proposed framework is presented. Section 4 suggests the robust counterpart models. The application of the proposed framework to an urban district in Tehran is analyzed in Section 5. Conclusions are drawn in Section 6. Tehran and the position of major faults [START_REF] Asadzadeh | An augmented approach for measurement of disaster resilience using connective factor analysis and analytic network process (F'ANP) model[END_REF] 2. Literature review

The present work mainly focuses on search and rescue operations in the earthquake response; but in a broader context, it is related to important decisions about resource allocation, scheduling, and routing in disaster response. The response phase has attracted the attention of the majority of disaster management contributors [START_REF] Goldschmidt | Humanitarian operations and crisis/disaster management: A retrospective review of the literature and framework for development[END_REF], and there is a rich body of research dealing with its different aspects. To not lose focus, the review mainly concerns papers that have studied allocation, scheduling, or routing decisions within the scope of SAR.

Resource allocation

Zhang et al. ( 2012) proposed a deterministic model for the allocation of multiple resources to multiple affected nodes to fully serve demand after primary and possible secondary disasters. They applied a heuristic approach to minimize total costs of resource allocation to primary disasters and opportunity costs of arrival at possible secondary ones. [START_REF] Chu | Post-earthquake allocation approach of medical rescue teams[END_REF] addressed the allocation of medical teams shortly after a severe disaster. They defined saving functions for different degrees of severity and maximized the expected number of survivals via a nonlinear mathematical allocation of available teams at different points. [START_REF] Zhang | A two-stage resource allocation model for lifeline systems quick response with vulnerability analysis[END_REF] developed a deterministic two-stage allocation model for restoring interdependent lifeline infrastructure after disruptions. A cost minimization model selects the destroyed components in the first stage, and then restoration scheduling is determined by a makespan minimization model. A heuristic algorithm was proposed to solve the proposed models. [START_REF] Su | Multiple emergency resource allocation for concurrent incidents in natural disasters[END_REF] addressed the collaboration of various emergency agencies (hospitals, police stations, etc.) to respond to concurrent incidents in disaster zones. They chose the sum of travel times and allocation costs as an objective function, treated the parameters as deterministic, and developed a differential evolution (DE) algorithm as their solution approach. [START_REF] Xiang | A medical resource allocation model for serving emergency victims with deteriorating health conditions[END_REF] developed a queuing network for patients with deteriorating conditions after large-scale disasters and proposed two optimization models for death rate minimization and system time minimization.

Scheduling and Routing

As suggested by [START_REF] Wex | Emergency response in natural disaster management: Allocation and scheduling of rescue units[END_REF], building a service tour for each resource unit in the affected district can be viewed in both the scheduling and routing domains. [START_REF] Fiedrich | Optimized resource allocation for emergency response after earthquake disasters[END_REF] developed a dynamic optimization model to schedule technical equipment in SAR, and also considered construction stabilization and rehabilitation of transportation lifeline tasks to minimize fatalities. [START_REF] Chen | Optimal team deployment in urban search and rescue[END_REF] developed a multi-stage stochastic program to optimize the routing of homogenous SAR teams in a large urban zone to maximize the number of saved people. Taking the arrival time of SAR teams as a priori vs. considering demand, service times, and travel times as random variables, they proposed a sequence of interrelated two-stage stochastic programs as their solution approach and determined the tours of SAR teams. [START_REF] Zheng | Disaster rescue task scheduling: An evolutionary multi-objective optimization approach[END_REF] presented a multi-objective fuzzy rescue task scheduling and risk minimization formulation and solved it with a multi-objective biogeography optimization (BBO) algorithm. They treated task weights, risk factors, and processing and travel time parameters as fuzzy numbers and considered a nonlinear fuzzy risk exposure function for rescue units. [START_REF] Yan | Logistical support scheduling under stochastic travel times given an emergency repair work schedule[END_REF] proposed a cost minimization logistical support scheduling model with stochastic travel time in a given emergency repair network. Adopting average travel times, they simplified the stochastic model to a deterministic model and used a heuristic decomposition approach. [START_REF] Lei | Personnel scheduling and supplies provisioning in emergency relief operations[END_REF] proposed a deterministic service tardiness minimization model to find a timely dispatching schedule for medical supplies originating from different depots by preplanned routes of medical teams in the aftermath of a disaster. [START_REF] Cao | Emergency Rescue Vehicle Dispatch Planning Using A Hybrid Algorithm[END_REF] developed a deterministic multi-objective model for multi-period optimization of emergency rescue vehicles routing in post-disaster operations. In addition to the minimization of the rescue operation time, they considered the accumulated number of an unserved population at each period as the rescue delay cost objective function. Furthermore, by defining a rescue utility score for each demand node, they defined the rescue utility objective function. They used a hybrid NSGA-II and ant colony optimization algorithm to solve the model. [START_REF] Bodaghi | Probabilistic allocation and scheduling of multiple resources for emergency operations; a Victorian bushfire case study[END_REF] presented a deterministic completion time minimization model to determine the scheduling of expandable and nonexpandable resources in emergency operations. To count for the uncertain input parameters, their proposed framework runs the model for a set of plausible scenarios and selects the most common plan. [START_REF] Liu | Decision making on post-disaster rescue routing problems from the rescue efficiency perspective[END_REF] developed a two-phase efficiency-focused routing optimization model. They evaluated the rescue efficiency of each arc based on its travel time and risk as well as the affected population, secondary hazard probability, and damage severity by a DEA model in the first phase. Considering the loss of efficiency for the remaining arcs over the operation progress, the second phase maximizes the total efficiency score of the selected visit sequence. [START_REF] Wex | Emergency response in natural disaster management: Allocation and scheduling of rescue units[END_REF] proposed a deterministic binary quadratic formulation to determine service routes of rescue teams with different capabilities to minimize the sum of the weighted completion time. They addressed different greedy, construction, and improvement heuristics, in addition to a Mont Carlo simulation method, to solve their NP-hard model. Extending their research, [START_REF] Schryen | Resource planning in disaster response decision support models and methodologies[END_REF] presented two deterministic models to coordinate assignment and scheduling of multiskill rescue units, with and without collaboration, for serving multiple incidents with different requirements. They proposed a heuristic solution method, which was later improved by a branchand-price algorithm, developed by [START_REF] Rauchecker | An exact branch-and-price algorithm for scheduling rescue units during disaster response[END_REF]. [START_REF] Shahparvari | Possibilistic scheduling routing for short-notice bush fire emergency evacuation under uncertainties: An Australian case study[END_REF] addressed short-notice evacuation of the population in a bush fire disaster with hard time windows. They presented a possibilistic rescue routing model with the objective function of maximizing the number of transferred evacuees. Although their work is in the evacuation stream, it shares similarities with our focused area because it considers capacitated rescue resources with uncertain demand loads and route traveling times. Table 1 gives an overview of the main characteristics of the above research.

Although there is increasing interest in the planning of resources during SAR operations, the considered studies either overlooked existing uncertainty in the disaster response reality or tackled it mostly with fuzzy logic, stochastic optimization, and scenario-based robust optimization approaches. While these methods add value to deterministic models, which ignore the situation's inherent uncertainty, their results and successful application highly rely on decision-makers' knowledge of uncertain parameters in the form of membership functions, probability distributions, or possible scenarios [START_REF] Bozorgi-Amiri | A multi-objective robust stochastic programming model for disaster relief logistics under uncertainty[END_REF]. However, many underdeveloped or developing countries have difficulties in accessing well-recorded reliable historical data or high-precision simulations. Hence, determining these characteristics for post-disaster timesensitive SAR operations can lead to more inaccuracy or be operationally impractical or inefficient. Furthermore, the use of deterministic models with expected values is accompanied by the risk that solutions will be infeasible and/or non-optimal in cases of uncertain input parameter perturbations [START_REF] Mulvey | Robust optimization of large scale systems[END_REF][START_REF] Ben-Tal | Robust optimization[END_REF].

Such inefficiencies in the management of scarce rescue resources during disaster response with a slipping rate of survival can lead to higher losses caused by an inefficient allocation of resources among demand locations, waste of SAR teams' time in non-optimal routes, and higher safety risks for SAR teams. To close this gap, this paper develops a robust SAR planning model based on the interval robust approach proposed by [START_REF] Bertsimas | The price of robustness[END_REF] to address uncertainty in both the allocation and routing stages. Although a few papers on humanitarian relief [START_REF] Baron | Facility location: A robust Optimization Approach[END_REF][START_REF] Najafi | A multi-objective robust optimization model for logistics planning in the earthquake response phase[END_REF][START_REF] Zhang | A robust counterpart approach to the bi-objective emergency medical service design problem[END_REF][START_REF] Zokaee | A robust optimization model for humanitarian relief chain design under uncertainty[END_REF] have applied interval data robust formulations, to the best of the authors' knowledge, it has not been introduced in the SAR literature. As the primary mission of humanitarian operations is saving lives and alleviating human suffering, the majority of the studies focus on effectiveness rather than cost-related objectives (Falasca et al., 2009). Most of the previous studies have focused on fatality minimization or operation time minimization, and they have mostly neglected the fairness challenge in allocation and scheduling models. The majority of studies have not looked into the collaboration of nonidentical SAR resources in task completion. Furthermore, although task weights can represent task risk rankings, only a few studies have explicitly addressed secondary risks.

Our study is most similar to the work by [START_REF] Schryen | Resource planning in disaster response decision support models and methodologies[END_REF] in the assignment and routing of collaborating teams with different skills. We extend their model and contribute to the SAR literature by:  Incorporating realistic interval uncertainties and developing a novel robust allocation routing framework.  Addressing fairness by maximizing the minimum timely coverage of affected districts by SAR resources in the allocation stage.  Considering operation disruption by including secondary destruction risk of locations with or adjacent to unstable structures.  Coordinating available resources and command centers in different districts of the impacted zone in two ways. First, since urban zones are often districted, the presented decision support framework considers a multi-district multi-period setting. It better corresponds to the resource assignment decision process that arises at a tactical level in the response phase. Second, it incorporates the resource transshipement between affected districts for effective operations management and loss minimization.  Accounting for additional operational constraints in the routing of rescue teams including idling for rest periods and the impact of secondary destructions on the operation duration. As per the definition, the investigated problem belongs to the broader field of resourceconstrained routing and scheduling as well as resource-constrained project scheduling problems (Hartman and Briskorn, 2010;Paraskevopoulus et al. 2013). Addressing both allocation and routing issues in these problems makes them computationally complex problems that belong to NP-hard class [START_REF] Garey | Spatial modelling of urban physical vulnerability to explosion hazards using GIS and fuzzy MCDA[END_REF][START_REF] Miller | A binary-rounding heuristic for multi-period variable-task-duration assignment problems[END_REF][START_REF] Lenstra | Complexity of vehicle and scheduling problems[END_REF][START_REF] Laporte | What you should know about the vehicle routing problem[END_REF]. Furthermore, there is even incremental complexity due to the consideration of real-world constraints reflected from the interviews with the SAR experts, such as the arrival of new resources over the planning horizon, transshipment of resources between districts, resource collaboration in demand locations, and time-sensitive SAR times [START_REF] Schiffer | Vehicle routing and location routing with intermediate stops: a review[END_REF]). In the related literature, the assignment of resources and their service tours in SAR problems and the other variants of resource-constrained project scheduling or routing problems are modeled and solved following simultaneous or two-stages decomposition approaches (Paraskevopoulus et al. 2013;[START_REF] Yalçındağ | Pattern-based decompositions for human resource planning in home health care services[END_REF]. However, the intrinsic complexity of the problem, which contradicts the need for decision-making agility, especially in disaster response, has made the decomposition approaches efficient and interesting methods being applied in numerous studies [START_REF] Schmid | Rich routing problems arising in supply chain management[END_REF][START_REF] Rolland | Decision support for disaster management[END_REF]Nickle et al., 2012;[START_REF] Allaoua | A matheuristic approach for solving a home health care problem[END_REF][START_REF] Yalçındağ | Pattern-based decompositions for human resource planning in home health care services[END_REF][START_REF] Can | Multi-project scheduling with two-stage decomposition[END_REF][START_REF] Misir | An analysis of generalised heuristics for vehicle routing and personnel rostering problems[END_REF][START_REF] Yalçındağ | A decomposition approach for the home health care problem with time windows[END_REF][START_REF] Attia | A decomposition-based heuristic for large employee scheduling problems with inter-department transfers[END_REF][START_REF] Dahmen | A two-stage solution approach for personalized multi-department multi-day shift scheduling[END_REF][START_REF] Memari | Air and ground ambulance locationallocation-routing problem for designing a temporary emergency management system after a disaster[END_REF].

To overcome the drawback of simultaneous modelling of both decisions while benefiting from joint allocation routing decision-making, a two-stage framework is proposed in this study. The first stage takes arriving resource availability times as a priori and builds robust assignment of available and incoming teams to affected districts over a multi-period horizon for SAR response. This stage maximizes allocation fairness by maximizing minimum SAR demand coverage over the whole zone and determines: 1) the allocated number of resources in each district for each period; and 2) resource transfers to other districts or discharging decisions over the planning horizon. The outputs of the first stage for the first period of the planning horizon are fed into the second stage. The second stage models the routing of each district's allocated teams among its demand locations during the several hours' duration of that period. It integrates fatality and time minimization by taking the trapped population in each demand location of the district as the importance weight in an objective function of minimizing the weighted sum of task completion times. The model allows team collaboration and takes into account SAR tasks' non-preemptive characteristics. To have the optimal weighted sum of demand locations' SAR completion times, this stage determine 1) the sequence of demand location visits for each team; and 2) the amount of SAR activity of each team for each demand location on their routes.

After implementing the first period's allocation and routing decisions, the planning framework is run again for the remaining periods of the horizon. In the first stage, the allocation model is solved again to capture the actual status of operation progress and upcoming updates on input parameters. Afterward, with the new allocation scheme, the second stage is run with the updated data to optimize the routing of allocated teams for the next period. This cycle continues until the end of the SAR planning horizon in the response phase.

Proactive approaches in the mitigation and preparedness phases play a key role in reducing disaster loss in the built environment [START_REF] Mojtahedi | Stakeholders' approaches to disaster risk reduction in built environment[END_REF]. Therefore, the proposed decision support framework is designed, not only to help SAR planning during the response phase but more importantly, to also equip relevant disaster stakeholders with a tool for more effective use of resources in investment, training, and capability-building decisions during the risk mitigation and preparedness phases.

Problem description

This study proposes a two-stage approach to optimize the planning of SAR operations in disaster-impacted zones. It assumes that the affected zone is divided into several districts. Each district, with several demand locations, has three types of demand based on the damage severity and requires teams with a compatible set of skills to serve them.

The first stage addresses the coverage of affected districts' SAR demand over the entire SAR timeline in the response phase and formulates it as a multi-period fair resource allocation model (MPFRAM). For each severity type, the model takes the total estimated SAR demand in each district and computes the minimum demand coverage over the entire affected zone. Based on the importance weight of each demand type, the objective function maximizes the allocation fairness by maximizing the weighted sum of all demand types' minimum coverage over the affected zone.

To account for allocation effectiveness about deteriorating chances of saving lives, a decreasing utility parameter is considered for allocation in each period. Incorporation of this parameter, which is driven by the declining survival rate of trapped people after large-scale disasters [START_REF] Fiedrich | Optimized resource allocation for emergency response after earthquake disasters[END_REF], represents the importance of operation effectiveness and favorability of quick demand response [START_REF] Gralla | Assessing trade-offs among multiple objectives for humanitarian aid delivery using expert preferences[END_REF][START_REF] Erbeyoglu | A robust disaster preparedness model for effective and fair disaster response[END_REF]. Having the utility parameters as the weights, each district's demand satisfaction is calculated as the weighted sum of allocated SAR work over the whole horizon. To have optimal utilization of capacities against demand, the model allows for an exchange of SAR units between districts wherever it is beneficial, according to demand coverage status and district distances. As SAR operations evolve and response time shrinks, fewer SAR teams are required and the affected zone prepares for debris management and road restoration. So, the model considers both the deployment and release of SAR units in decision variables. However, due to the usual scarcity of resources vs. overwhelming demand, releasing rescue SAR units are not expected unless in the last periods. The main assumptions and outcomes of this stage can be summarized as follows.

Assumptions:

 SAR demand and travel time parameters are considered uncertain and assumed to be symmetrically distributed at the given intervals.  The number of arriving new SAR teams at each period is assumed to be known.  Each period consists of T hours (e.g., 12 hours), and T can be set based on decision environment circumstances.  The utility of resource allocation in each period represents that period's demand coverage efficiency and declines over time.  SAR teams can travel between districts to improve SAR capacity utilization vs. demand.

However, the travel time between resource-exchanging districts decreases the amount of possible SAR work.

Outputs:

 The allocated number of resources in each district for each period.

 The resource transfer to other districts or discharging decisions over the planning horizon.

The outputs of the first stage for the first period are fed into the second stage. This stage addresses the deployment of allocated teams in the first period of the planning horizon for each district. It formulates the problem as a single period capacitated resource routing model (SPCRRM) and optimizes the weighted sum of SAR task completion times in the district's demand locations. Each demand location needs the SAR operation in one of the three types of demand based on its damage severity. Similarly and with hierarchical levels of skills, each team has a specific capability degree, which qualifies them for serving the respective demand type and less severe situations. Due to the uncertain changing environment in the aftermath of natural disasters, SAR operation times, the number of trapped people, and the travel time between demand locations within each district are assumed to be uncertain, and an interval robust approach is adopted to deal with that.

Some of the demand locations, which include residential buildings, schools, hotels, malls, and urban infrastructure, may be at a higher risk of further destruction and extended operation time due to situations like being located on the coastline with tsunami risk, or being adjunct to partiallycollapsed unstable structures. To adapt the realistic vulnerability of rescue success to operation timing in such situations, an estimated threshold time is considered for high-priority demand locations. These locations, which should be identified based on expert input, are assumed to require longer operation time if they are not served before their threshold time.

The amount of required operation time at each demand location is stated in terms of hours. It may be accomplished by several teams and can be estimated based on building dimensions, approximate trapped population, damage severity, etc. The model adopts the realistic constraint of work-rest periods in difficult working conditions for SAR operations. This assumption becomes vital in team scheduling when the disaster-impacted zone is exposed to infectious disease or dangerous materials [START_REF] Sawik | An integer programming approach to scheduling in a contaminated area[END_REF]. The maximum limit of working hours and the rest duration is assumed a priori, but teams are not allowed to interrupt an unfinished task, even if it exceeds their shift time. If several teams serve at a single demand location, their service times should be overlapping, or follow an end-to-start relationship with at least one other team until the task is completed. Any task that is planned to be started during that period should be completed and removed from the district's total SAR demand. The start time of each team on the last location of their tour, along with their assigned SAR activity hours, determine their availability time for the next period. The main outputs of this stage are as follows:

 The sequence of demand location visits for each team.

 The amount of SAR activity of each team in each demand location on their route.  The expected operation finish time for each demand location based on the start of SAR and the team engagement amounts. This also determines team availability times for the next period, whether in the same district, in the transfer link's destination district, or out of the SAR mission. Due to the very vague status of input parameters over the whole search and rescue phase (120 hours), a dynamically changing environment, and the rapidly incoming updated data, SPCRRM considers a single period of 12 hours duration. After implementing the allocation and routing decisions of the first period, the planning horizon is updated and the proposed framework is rerun to capture the actual operation progress and information updates. By the inclusion of realized remaining demand in each district and other input parameter updated vales, MPFRAM is run for the remaining periods of the horizon. Building on the outputs of the first stage and benefiting from the newly arrived data, the second stage SPCRRM is run for the next period to optimize the routing of allocated teams in each district. This cycle is repeated until the end of the SAR timeline in the response phase. The flow of information between the two stages is outlined in Fig. 2. Feeding updated input data into the models requires a well-designed and stable communication infrastructure, and can highly impact the effectiveness of search and rescue operations [START_REF] Mohamadi | Fuzzy multi-objective stochastic programming model for disaster relief logistics considering telecommunication infrastructures: a case study[END_REF].

Database

The current period's allocation decision

Execution Stage

Implementing the current period's decisions

Fig. 2. SAR two-stage decision support framework

The next two sub-sections present the formulations of the problem as mixed-integer linear programming (MILP) models.

First stage -MPFRAM

The mathematical model in the first stage uses the following notation and decision variables. If there is an increase in the number of allocated teams for demand type s in the affected zone from period t-1 to t Otherwise

The model is formulated as follows:

Max 𝑤 𝑐𝑜𝑣 (1) s.t. 𝑥 = ℎ -𝑓 + 𝑦 - 𝑦 ∀𝑠, 𝑘, 𝑡 = 1 (2) 𝑥 = 𝑥 + ℎ -𝑓 + 𝑦 - 𝑦 ∀𝑠, 𝑘, 𝑡 > 1 (3) ℎ ≤ 𝑐 𝑧 ∀𝑡, 𝑠 (4)
𝑓 ≤ (1 -𝑧 ) 𝑐 ∀𝑠, 𝑡 (5)
𝑐𝑜𝑣 𝑑 - 𝑇𝑥 - 𝑔 𝑦 𝜃 ≤ 0 ∀𝑠, 𝑘 (6)
𝑥 , 𝑦 , ℎ , 𝑓 ≥ 0, 𝑖𝑛𝑡, 𝑐𝑜𝑣 ≥ 0

(7) 𝑧 ∈ {0,1} (8)
To account for fair allocation, the objective function maximizes the demand coverage for the district with the lowest coverage score. This also results in minimizing the gap between the coverage statuses of different districts, and hence, fairness is taken into consideration. Since there are different demand types for severity, the objective function takes the weighted sum of all three types' minimum level of coverage over all the districts in the affected zone.

Constraints (2) and (3) calculate the number of SAR teams in each district. Each district may cover its demand from multiple sources, including teams allocated in the previous period, transferred teams from other districts, and allocation of new teams that reach to the hot zone at each period. Constraint (4) is a capacity constraint for the total number of new teams assigned in each period over the whole affected zone. As operations progress, for each demand type for which the amount of allocated resources exceeds the remaining total demand, the operation center should decide to discharge excess teams (Constraint ( 5)). These teams can either be assigned as fresh resources for satisfying the remaining demand in less severe tasks and be taken into the updated 𝑐 parameter, or dispatched for further operations of the response phase, like relief item distribution.

Constraint (6) computes the minimum coverage level for each demand type over the affected districts. Although the transshipment of resources enables better resource utilization, the time spent on traveling between the districts reduces the remaining time for the operation in the destination district. Hence, the lost time is subtracted from the resources' effective allocated times. To incorporate the reduction in the survival rate of trapped people over response time shrinkage, the effectiveness of allocated resources in each period is considered by the help of utility weights. Sign constraints are given in Constraints ( 7) and ( 8).

After solving the first stage, the number of allocated teams for each demand type and each district in the first period is fed into SPCRRM in the second stage to plan the teams' SAR activities in each district. It is a single-period capacitated routing model that schedule the operations of the allocated teams in each district for the duration (e.g. 12 hours) of the first period of the remaining planning horizon.

Second stage -SPCRRM

The mathematical model in the second stage uses the following notation and decision variables. The objective function (9) minimizes the weighted sum of SAR completion times in affected locations. Since the populations of trapped people in the demand locations are used as importance weights, the objective function can be interpreted as the total man-hours waiting to receive SAR services over all the locations in the district. Constraint ( 10) is an integrality constraint that makes sure that locations with assigned SAR teams are selected in the routing plan. Constraint (11) guarantees that the start time of the SAR mission by each team at each selected demand location cannot exceed the duration of the planning period. Each started SAR activity will continue until completion. The completion time, which might be later than the period's duration, determines the teams' availability times in the next period. Only qualified teams may get involved in serving each demand type, and this case is addressed in the whole model via appropriate condition(s). Teams are ready to start their route once they arrive at the command base, as shown in Constraint (12). Constraint ( 17) ensures that each selected demand location is fully served by contributing teams. Constraint (18) guarantees that if the operation starts at a location, there should be at least one team present to complete the job; in other words, the work cannot be interrupted, and the involved teams may not leave the demand location prior the SAR completion unless another team continues the mission.

𝑏 ≥ 𝐿𝑞 ∀𝑖, 𝑗 ; 𝑗 ≠ 1; 𝑗 ≠ 𝑚 (19) 𝑠′ ≥ 𝑦 𝜖 ∀𝑗; 𝑗 ≠ 1; 𝑗 ≠ 𝑚 (20) 𝑠′ ≤ 𝑇𝑦 + 𝜖 1 -𝑦 ∀𝑗; 𝑗 ≠ 1; 𝑗 ≠ 𝑚 (21)
Constraint ( 19) states that once a team is allocated to a demand location, they should serve in the mission for a minimum amount of time. Constraints (20) and (21) determine whether each location receives its required SAR before the threshold time. Constraints ( 22) to (26) compute the total continuous work of each team, along with their route after visiting each demand location. As shown in Constraints ( 22) to ( 24), the consecutive work of a team after a location is equal to their SAR activity in that location only if this node is visited right after a rest time or at the start of the route. Constraints ( 25) and ( 26) accumulate the SAR activities of the team at the linked locations along the route if the team's working shift is not finished yet. Constraints ( 27) and (28) determine the necessity of rest time for each team after each SAR visit. Constraint (29) ensures the integrality relation of assigning the team for SAR activity on each demand location and their work amount.

Constraint (30) calculates the start time of operation by each team at each demand location along their route. Constraints ( 31) and ( 32) ensure that each selected demand location has a predecessor and successor node on the team's route. The finish time of SAR at each location is calculated in Constraint (33) as the latest time of operation completion by the assigned teams. While the start of SAR at each selected demand location should be within the routing period's duration, the finish time might be later, since no interruption is allowed until the selected location's total demand is satisfied. Sign constraints are given in Constraints (34) and ( 35).

The solution provided by solving this stage flows back to the first stage to complete the planning framework. This feedback, along with the externally updated information, enables the first stage to re-optimize the allocation decisions for the remaining horizon and to continue the operation planning.

In the ever-changing disaster environment, the inherent uncertainty of parameters limits the applicability of deterministic models. Hence, this section introduces the model's interval-based robust counterpart, which keeps the model flexible and robust in real-world applications.

Interval data robust optimization

In practice in an uncertain environment, there are parameters with unknown exact values that may be actualized differently from their estimated nominal values. These gaps are potentially capable of violating some constraints and hurting solution optimality. Therefore, there is a need for implementing solution approaches that remain robust in terms of feasibility and solution quality in the presence of uncertainty [START_REF] Bertsimas | The price of robustness[END_REF]. [START_REF] Soyster | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF] was the first to develop a model to meet this need. However, his model, which constructed a solution feasible for all the data in a convex set, was too conservative. This first attempt was further developed by [START_REF] Ben-Tal | Robust convex optimization[END_REF][START_REF] Ben-Tal | Robust solutions to uncertain programs[END_REF], 2000), El-Ghaoui and Lebret (1997), El-Ghaoui et al. (1998), [START_REF] Bertsimas | The price of robustness[END_REF] and [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF]. [START_REF] Bertsimas | The price of robustness[END_REF] proposed a solution approach for linear mathematical models with an uncertain coefficient matrix. Their approach does not increase the model's complexity and has a flexible in a conservatism level; it is briefly described below. [START_REF] Bertsimas | The price of robustness[END_REF] considered the following model and assumed that the i-th constraint has ⌊𝛤 ⌋ uncertain technological coefficients with values that lie in a stochastic symmetric interval 𝑎 -𝑎 , 𝑎 + 𝑎 , and exactly one parameter with a 𝛤 -⌊𝛤 ⌋ perturbation, where 𝛤 is not necessarily an integer and may take up any value in [0, |𝐽 |] with 𝐽 being the set of uncertain coefficients in constraint i.

Max z=cx 𝑎 𝑥 ≤ 𝑏 l≤ 𝑥 ≤ 𝑢 (36)
To ensure the model stays feasible after the Γi coefficients' perturbations, they rewrote the model with the help of a protection function 𝛽 (𝑥 * . 𝛤 ):

Max z=cx 𝑎 𝑥 + 𝛽 (𝑥 * , 𝛤 ) ≤ 𝑏 𝑙 ≤ 𝑥 ≤ 𝑢 (37)
where

𝛽 (𝑥 * , Γ ) = max 𝑆 ∪ {𝑡 } 𝑆 ⊆ 𝐽 , |𝑆 | = ⌊𝛤 ⌋. 𝑡 ∈ 𝐽 ∖ 𝑆 ∑ a x * + ∈ (𝛤 -⌊𝛤 ⌋)a x * .
𝛽 (𝑥 * , Γ ) protects the feasibility by embedding the largest possible deviation of left-hand side coefficients in i-th constraint. So it needs to find a Γi set of coefficients that can cause the biggest perturbation:

𝛽 (𝑥 * , 𝛤 ) = 𝑀𝑎𝑥 𝑎 𝑥 * 𝑧 ∈ s.t. 𝑧 ∈ ≤ 𝛤 0 ≤ 𝑧 ≤ 1 ∀𝑗 ∈ 𝐽 (38)
In this model, positive values of 𝑧 represent the coefficients in the 𝛤 uncertainty level whose perturbation can cause the biggest risk of violating the right-hand side in the i-th Constraint (36) and ( 37). Since this model is feasible and bounded, according to the strong duality property, its dual problem is also feasible and bounded, and its objective value is equal to the primal. If 𝑧 and 𝑝 are the auxiliary dual variables of Model (38), one can rewrite Model (37) into its ID robust counterpart model (39) by replacing 𝛽 (𝑥 * . 𝛤 ) with the dual model of Model (38). The use of 𝛤 , which is called "the budget of uncertainty" (BoU), facilitates user flexibility in calibrating the robustness of the method based on the preferred conservatism level [START_REF] Bertsimas | Theory and applications of robust optimization[END_REF].

Max 𝑧 = 𝑐𝑥 s.t. 𝑎 𝑥 + 𝑧 Γ + 𝑝 ∈ ≤ 𝑏 𝑧 + 𝑝 ≥ 𝑎 𝑦 -𝑦 ≤ 𝑥 ≤ 𝑦 𝑙 ≤ 𝑥 ≤ 𝑢 𝑧 , 𝑝 , 𝑦 ≥ 0 ∀𝑖, 𝑗 ∈ 𝐽 ∀𝑗 ∀𝑗 (39)
Using 𝛤 as the allowed summation of perturbations in uncertain parameters of row i, this approach enables the user to make the trade-off between the probability of violating the i-th constraint and the effect on the objective function of the nominal problem [START_REF] Bertsimas | Theory and applications of robust optimization[END_REF]). [START_REF] Bertsimas | The price of robustness[END_REF] proved that if more than 𝛤 uncertain parameters change in the i-th constraint, the probability of the constraint's violation is at most 𝐵(𝑛, 𝛤 ), which is approximated by the following bound:

𝐵(𝑛, 𝛤 ) ≈ 1 -𝜙( √ )
where 𝑛 = |𝐽 | and 𝜙( ) is the cumulative distribution function of a standard normal variable.

Hence, one can protect the solution against infeasibility by choosing an appropriate level of 𝛤 , depending on the number of uncertain parameters present in that constraint. Testing different pairs of (𝑛, 𝛤 ) , the proposed approach is capable of guaranteeing a robust feasible solution without sacrificing too much optimality.

In the case of constraints that are similar to Constraint (6), in which there is only one uncertain demand parameter per constraint for each (𝑠, 𝑘), Bertsimas and Sim's approach results in several BoU parameters 𝛤 , each ranging from 0 to 1 for the demand uncertainty in each constraint. This means that even examining only 0 and 1 options for each 𝑑 's BoU leads to 2 * states, which is computationally infeasible to deal with. Therefore, [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF] extended the work of [START_REF] Bertsimas | The price of robustness[END_REF] to account for such situations. Their proposed robust method considered an accumulated conservatism parameter 𝛤 for the relevant set of uncertain parameters, which are dispersed among different rows of constraints.

In the above example, s denotes the demand type and remains in the conservatism level definition. While, the notion of the district is removed, and the method considers a single 𝛤 to accumulate the budget of uncertainty for all districts' demand type s. Therefore, each district's share of demand uncertainty becomes , where R is the number of districts.

In line with the characteristics of both studies, this paper adopts the interval data robust method presented by [START_REF] Bertsimas | The price of robustness[END_REF] and [START_REF] Bertsimas | A robust optimization approach to inventory theory[END_REF] to address uncertain parameters in both stages. SAR demand and the travel time parameters in the first stage, and trapped population, required SAR times, threshold times, and travel time parameters in the second stage are treated as uncertain and are considered in symmetric interval forms. These intervals are fed in by damage estimation applications and systems, expert input, and information updates received during the operations.

MPFRAM robust counterpart:

𝛤 is the number of districts with stochastic travel times from district k. 𝑟 and 𝑣 are dual variables of the protection function of districts' travel times. Since there are R districts considered in the first-stage allocation formulation, is defined as the level of conservatism for demand parameters as stated in Equation ( 40). Furthermore, 𝑞 and 𝑢 are the corresponding regular dual variables for protecting the model's feasibility against demand uncertainty. Equations ( 40) to ( 43) substitute Constraint (6) and build the first stage's robust counterpart.

SPCRRM robust counterpart:

𝛤, which is independently defined for each district in the second stage, is the number of demand locations with uncertain populations. 𝜁 and 𝜁 are defined to address duality in the protective function of population coefficients. Equation ( 9) is replaced with Equations ( 44) and (45). To adopt [START_REF] Bertsimas | The price of robustness[END_REF] approach in the uncertain coefficients in the objective function, the worst case occurs when all the locations have higher populations, worsening the weighted summation of completion times. So, in the presence of conservatism level 𝛤, the objective function can be written as Min = ∑ 𝑝 𝑓 + 𝛽 𝑓 * , Γ , and the rest follows the approach introduced by [START_REF] Bertsimas | The price of robustness[END_REF].

Min
𝛼 + 𝛼 ≥ 𝑏 ′ 𝑞′ ∀𝑗 𝑗 ≠ 1, 𝑚 (47)
𝛽 + 𝛽 ≥ 𝑦 𝑎 ∀𝑗 𝑗 ≠ 1, 𝑚 (48)
Similarly, is considered as the conservatism level for uncertain threshold times in demand locations with the risk of operation expansion. Since uncertain threshold times appear in two different sets of constraints, two protection functions are required. One addresses the feasibility of Constraint (20) with dual variables 𝛾 and 𝛾 by Equations ( 50) and ( 52); and one protects Constraint (21) against the threshold time's uncertainty by 𝛾′ and 𝛾′ in Equations ( 51) and ( 53).

There are (𝑚 -2) possible independent outflow links between the demand locations in the second stage; and therefore, ( ) defines the level of conservatism for travel times for each pair of locations. Since in Equation ( 30), the uncertain travel time parameter is not multiplied by any variable, the constraint's feasibility protection is obtained by the inclusion of travel time perturbation based on its budget of uncertainty. Replacing Equations ( 9), ( 17), ( 18), ( 20), ( 21) and ( 30) with Equations ( 44) to ( 55) builds the second stage's robust counterpart.

Case study

To show the effectiveness of the proposed model, this section describes a robust SAR allocation and routing example based on an earthquake scenario in an urban part of Iran's capital, Tehran. It is the most populated city in the country, with 22 municipal districts, and is located in the southern foothills of the Alborz Mountains, in an earthquake-prone region (Fig. 3). It is surrounded by three major seismically active faults, Mosha to the northeast, North Tehran to the north, and Rei to the south, as well as hidden faults underneath the whole city's sediment layers (JICA, 2000). The city has expanded tremendously over the past 50 years, and now has more than 8,000,000 inhabitants and covers an area of 730 km 2 . It is now considered one of the cities most vulnerable to earthquakes in the world [START_REF] Amini Hosseini | Main challenges on communitybased approaches in earthquake risk reduction: Case study of Tehran, Iran[END_REF].

The scenario considered, called the Tehran floating earthquake scenario (TFES), was developed to assess the city's relative vulnerability and considers the hidden faults underneath the city. In the case of an earthquake caused by these faults, which are simulated to be 13 km in length and 10 km in width, everywhere in the city would have a similar probability of occurrence, with a seismic intensity of 8 to 9 on the MMI scale (JICA, 2000;[START_REF] Ashtari | Microseismicity in the region of Tehran[END_REF]. operations. The number of available teams for the lowest demand type (least severe) is relatively higher than for the other two since volunteer citizens from the neighborhood can cooperate in these cases. For more serious damage, the district's capacity might not be sufficient, and the operations will depend on the arrival of professional resources from other districts, cities or countries, depending on the demand size.

Fig. 4. First-stage results

As can be seen, to restore the maximum minimum weighted coverage among the subdistricts, subdistricts 4 and 6, with relatively more extensive damage, receive a share of available capacity.

Sensitivity analysis

To check the results' sensitivity to different situations, different levels of demand, and capacity are tested. Fig. 5 illustrates the impacts of both SAR resource capacity and improving their availability lead time for subdistrict demand coverage. This analysis, which can be used as a simulation tool, is capable of providing decision-makers with helpful managerial insights into "next steps" for improving the safety of the city's residents. Following are the practical insights that can be deduced from sensitivity analysis of the studied case: Severity degree 3

SD1 SD2 SD3 SD4 SD5 SD6 Weighted coverage
 For all three severity degrees, outcomes are more sensitive to the availability lead time of SAR teams than to increased numbers of resources.  Highly qualified teams that are competent to be involved in severity type 3 missions are more behind on demand compared with the other types; hence, improvements in either capacity or arrival lead time of these teams can significantly improve results.  Although it is very difficult to build capacity or decrease availability lead time of highly qualified teams for all the subdistricts in the district, it is vital for the district's SAR success that the most vulnerable subdistricts (4 and 6) be marked for advanced development and training of potential local and neighbor professionals. As can be seen in the chart, there is a huge 40% opportunity that can be unlocked by a 20% increase in the number of teams and cutting the availability lead time by half.  Lack of visible improvement in the results for the second severity degree group highlights the first priority, which is to fill the major gap between demand and SAR capacity.  Considering that it is easier to build teams for the first-degree severity type, the analysis confirms at least 30% potential improvement if each subdistrict trains eligible and sufficient local teams. Since the main damage to the districts' residents' safety, accessibility of roads, and operability of urban services is caused by medium and severe destruction, the sensitivity of SAR success to the vulnerability of structures to destruction is studied. Fig. 6 presents a simulated glide path for improving district coverage by decreasing damage to buildings (within the same SAR capacity). This glide path, together with district vulnerability and strategic character, can help to leverage relevant managerial decisions about reinforcement/renovation timelines and priorities. Using the proposed allocation model as a simulation tool for sensitivity analysis of SAR demand coverage before an earthquake strikes can be beneficial for disaster management authorities in both the risk mitigation and preparedness phases.

Impact of different levels of uncertainty budgets

Furthermore, different levels of uncertainty budgets are used to show the trade-off between the model's optimality and robustness in the allocation stage (Table 3). This analysis can help decision-makers with setting the model's conservatism to maintain feasibility and closeness-tooptimality when faced with uncertainty. As expected, eliminating deterministic and worst-case scenarios, the objective function follows a mild worsening trend when the model is immunized for higher perturbations and uncertainty budgets.

Performance analysis of proposed robust MPFRAM

To test the performance of the decisions suggested by the proposed robust model vs. the deterministic approach, additional experiments are implemented. To this end, at perturbation levels of 20% and 30%, models with three sets of uncertainty budget s(𝛤 , 𝛤 ) equal to (6, 5), (4, 3) and (3, 2.5), leading respectively to 99%, 90% and 80% feasibility reliability, are implemented. For each perturbation level, the performance of solutions provided by the robust models and the deterministic approach is then evaluated vs. 20 uniformly generated random realizations of the uncertain parameters from their respective intervals. The outputs are compared based on their performances after inserting the proposed allocation decisions in Equations ( 1) to ( 8). For each model, three sub-measures are calculated to define total performance: 1) the new objective function value from the realized 𝑐𝑜𝑣 ; 2) the coverage unfairness at each severity degree based on the standard deviation of realized coverages between the subdistricts; and 3) the violation of Constraint (6). The latter is the difference between a model's output for guaranteed coverage and its realized possible coverage of realized demand by the proposed allocation decisions. This infeasibility, which is of very high importance, means missed coverage promises and represents distorted resource requirement assessment.

To penalize unfavorable sub-measures and to unify the scale of all three parts of the measure, the total of realized demands for each type is used as a multiplier. The measure's mean and its range of changes for different random realizations are calculated in each model.

The results presented in Fig. 7 demonstrate that the robust approach outperforms the deterministic model. In all cases, the robust approaches show a minimum performance higher than the deterministic model's average measure value. Also, this superiority is even reinforced in more volatile environments with higher possible perturbations. At both uncertainty levels, the performance gap between the robust and deterministic approaches grows as the uncertainty budget increases. In other words, higher conservatism and solution immunization lead to lower deviations from promised coverage levels and to fewer resource miscalculations. While the appropriate conservatism level depends on the decision-maker's preference, the assessment shows that the robust allocation model presents a better fit for resource planning in SAR operations.

Second-stage results

To test the model's performance in the second stage, subdistrict 4 (Fig. 3), with an area of 4.9 km 2 and 52,000 habitants, is selected. We consider six demand locations in this subdistrict according to the municipality neighborhood definitions. Demand location 1, with a high-risk CNG station and in proximity to four very tall buildings, is at risk of secondary destruction after the earthquake, and hence is considered for operation time expansion after the threshold time. In the second stage, (6, 4, 1, 1, 10) are chosen as (𝛤, 𝛤 , 𝛤 , 𝛤 , 𝛤 ). This means that all demand locations' population parameters are assumed to be different from the nominal values. 66% of the locations have perturbed SAR processing time, only one of the locations with secondary destruction risk is assumed to have additional processing time and a threshold time different from nominal values, and finally, the travel time on 10 travel links have perturbations. Five teams are considered as the available resources for SAR routing. Teams 1 and 2 are capable only for type 1 severity, teams 3 and 4 are qualified for type 2 severity, and team 5 is capable of the most severe demand type. The other input data 1 used to model the routing can be found in Table 4. The second-stage model is run for a 12-hour period in Subdistrict 4, with 6 locations and 5 available teams, and the routing plan for allocated teams is determined. Fig. 8 shows the activity sequence for each SAR team in Subdistrict 4. Teams 1 and 2 can only satisfy the SAR demand at the first location due to the demand severity compatibility. To shorten the completion time, both teams need to work concurrently. Location 2, with demand type 2, receives its required resources for SAR from both teams 3 and 4. Teams are not allowed to interrupt the operations unless another 1 Data are altered to maintain confidentiality. team gets involved. That is why the model decides to release team 3 to respond to locations 4 and 5, but keeps team 4 involved to fully satisfy location 2's demands. Tasks that take the resource capacities beyond the current planning period get introduced in the next period routing model as new tasks with their remained processing time. Team 3 enters a rest period after it leaves location 2 since it gets to its shift limit. It then continues to visit locations 4 and 5. Team 5 is the only team capable of handling the severe demand type. Hence, this resource capacity is not sufficient for both locations 3 and 6, so the model chooses to serve location 3 due to its higher population. And hence, location 6 remains unserved.

Fig. 8. SAR team' activity sequences in Subdistrict 4

Sensitivity analysis

To estimate the operations' success in minimizing the weighted SAR completion times under different circumstances, sensitivity analysis is completed in different categories. Table 5 shows that the delayed availability of SAR teams can worsen the objective function values. Since location 1, with the first severity level, is at risk for operation time expansion by delayed operation and has a similar population to location 3 in the third severity category, scenarios 3 and 4 lead to similar weighted sums of completion times. Due to different ratios of resource to demand in each severity degree, more insights can be drawn by a detailed investigation of the impact of different availability scenarios on teams' workloads and missed survival opportunities. As can be seen, concurrent delays of all teams without injecting new resources, while keeping the same number of served demand locations, lead 2,00 4,00 6,00 8,00 10,00 12,00 14,00 16,00 18,00 20,00 22,00 Team 1 Team 2 Team 3

Team 4

Team 5

Hour to increased exposure of the task-force due to escalated demand and less remaining time. Moreover, depending on the delayed teams' capability mix, the summation of delayed starts of rescue operations can become higher than arrival delays. It must be addressed during response simulations before disaster strikes. Because there is a higher risk of resource arrival delays in the initial periods, which have the best survival chance.

Although the scenario simulation corresponds to a nighttime earthquake, the study of a daytime occurrence is equally important. Population distribution and density in different places (e.g., residential buildings, educational facilities, offices, or public spaces) differ during the day. Hence, the city's risk exposure and vulnerability can significantly change depending on earthquake onset time [START_REF] Alexander | The health effects of earthquakes in the mid-1990s[END_REF][START_REF] Freire | Integrating population dynamics into mapping human exposure to seismic hazard[END_REF][START_REF] Ara | Impact of Temporal Population Distribution on Earthquake Loss Estimation: A Case Study on Sylhet, Bangladesh[END_REF]. Also, the daytime transportation flow can complicate route accessibility and prolong rescue operations. Given the share of the studied district in the city's daily trips, it is important to evaluate how changes in its connectivity to other parts of the city, which lead to different availability lead times of rescue teams and travel times within rescue tours, can impact operation planning and demand fulfilment.

To this end, greater disruptions in accessibility leading to longer arrival and travel times of teams are tested. Fig. 9a demonstrates that each demand location experiences delayed completion times in the presence of different accessibility prolongation situations. Location 1, with timesensitive demand, experiences the highest impact. Since the disruptions are considered to impact all teams and travel times equally, the delay trends for the rest of the locations follow their priority status, depending on the population and availability of compatible resources. The increasing trends of total delayed completion times are shown in Fig. 9b. This analysis, coupled with survival rate vulnerability, can be used to assess the need for urgent transportation speed-ups.

Impact of different levels of uncertainty budgets

In this section, the impact of different levels of uncertainty budgets on uncertain parameters is tested. For this purpose, the impact of each uncertainty budget on population, travel time and SAR time in the robust model optimal solution is studied for three perturbation levels of 10%, 20% and 30% and at three constraint feasibility reliability levels of 99%, 90% and 80%. The reliability scores are calculated according to [START_REF] Bertsimas | The price of robustness[END_REF] 𝐵(𝑛, 𝛤 ) estimation. Additionally, the impact of concurrent variations of the above parameters in three reliability states is tested. The changes in the optimal solution of the robust model vs. the deterministic model are calculated and presented in Fig. 10. Since the population parameter has a bigger scale, its perturbations result in higher increases in weighted completion times, although the changes remain at almost the same level for different reliabilities. In contrast, while the impacts of SAR and travel times are smaller in magnitude, the solution demonstrates higher sensitivity to their perturbations and immunization efforts. For both SAR and travel time uncertainty, as the reliability level grows, the robust objective function value shows higher favorability vs. the deterministic solution by the increase of perturbations. That is, as the uncertainty budget increases, perturbations of these two uncertain parameters from their values cause greater increases in the objective function for immunization. So any effort to secure better estimates for these two parameters can avoid the destructive effects of their perturbations on the weighted sum of completion times, and allows a lower immunization cost. Dominated by the scale of the population parameter, joint variations in uncertain parameters are similar to the population in magnitude but exhibit higher sensitivity to reliability levels due to the presence of the other two parameters.

Performance analysis of proposed robust SPCRRM

To test the performance of the decisions suggested by the proposed robust routing model vs. the deterministic approach, additional experiments are carried out. First, for concurrent perturbations of uncertain population, SAR time and travel time at both 20% and 30%, three robust models with feasibility reliability levels of 99%, 90% and 80%, and respective uncertainty budgets are established. The outputs of the robust and deterministic models are collected and then tested vs. 20 uniformly generated random realizations of the uncertain parameters from their respective 20% and 30% perturbation intervals. The outputs are compared based on their performances after inserting the proposed routing decisions in Equations ( 9) to (35). For each 2) The infeasibilities in meeting the realized SAR time 𝑏′ of served demand locations vs. the proposed 𝑏 .

3) The infeasibilities in providing uninterrupted operation in a case of collaboration in Constraint (18). 4) The infeasibilities in the proposed operation start time 𝑠′ vs. the realized possible time

𝑠′

in Constraints ( 13) and ( 14). 5) The inefficiencies from both sooner-than-necessary and later-than-possible proposed rest times, based on the assessment of realized continuous workload 𝑤 in a case of following the rest visits 𝑟 proposed by the models. To penalize unfavorable infeasibilities, the exposed population multiplied by a penalty factor is considered and added to the minimization objective function. The measure's mean and its range of changes for different random realizations are calculated in each model and presented in Fig. 11.

The results demonstrate that the robust approach outperforms the deterministic model in all cases. The deterministic model's output varies in a wider range, which decreases its reliability even more. While the robust models maintain consistent performance in varied perturbations, the deterministic approach experiences a larger performance gap by when uncertainty increases. At both perturbation levels, the superiority of the robust model grows as the uncertainty budget and feasibility reliability increase.

The evaluation shows that for optimizing the routing of SAR resources in an uncertain disaster environment, the decisions proposed by the robust approach offer better robustness and optimality. However, the reliability level needs to be tuned based on the decision-maker's preference.

Concluding

This paper presents a new decision support framework with a two-stage approach for robust allocation and routing of rescue teams in an earthquake response environment. Both stages are designed to capture the operations' effectiveness in fighting against time in the battle to extricate trapped people. The first stage seeks a fair allocation scheme for SAR resources over the partitioned affected zone, and the second stage aims at minimizing SAR completion times at demand locations, weighted by the number of trapped people in each. To capture the real-world requirements for providing a robust solution in the presence of uncertain parameters, an interval data robust optimization approach is adopted. Demand, travel time, SAR processing time, and population are considered to be uncertain parameters. The presented decision framework is capable of providing managers with appropriate insights on key points of focus and priorities for improved results before and during operations in the mitigation, preparedness, and response phases. Belonging to the broader family of resource-constrained routing problems, several reallife constraints in the rescue teams' service tours (e.g., idling rest periods, time-dependent service durations, hierarchical skill levels, and collaborative rescue service) are reflected.

The model is applied to a numerical example inspired by an earthquake scenario in an urban district in Iran's capital, Tehran. Comprehensive analysis is done in both stages, and the sensitivity of the solution to different levels of resource capacities, availability lead times, travel times, and budget of uncertainty for uncertain parameters are studied. Additionally, in both stages, the superiority of the proposed robust models' performance vs. the deterministic approach in keeping a near-to-optimal and feasible solution under a perturbed environment is validated.

To solve larger problems, a heuristic or meta-heuristic algorithm needs to be developed and evaluated. Moreover, incorporating path selection for each pair of nodes in SAR routing, and integrating the allocation and routing of SAR teams with medical resource planning, crowd evacuation planning, or lifeline restoration planning in multi-objective setting are other interesting avenues for future research.

Fig. 1 .

 1 Fig. 1. (a) Map of faults and seismicity in Iran (Mojarab et al., 2014); (b) The 22 urban districts ofTehran and the position of major faults[START_REF] Asadzadeh | An augmented approach for measurement of disaster resilience using connective factor analysis and analytic network process (F'ANP) model[END_REF]

  Indices and parameters: 𝑖 SAR teams (𝑖 ∈ {1,2, . . , 𝑁}) 𝑗, 𝑗 Demand locations, where 1 and m are dummy nodes representing the command base (𝑗, 𝑗 ∈ {1,2, . . , 𝑚}) 𝑑 Location j's demand type 𝑑 ∈ {1,2,3} 𝑝 People trapped in location j 𝑐 Team i's capability type 𝑐 ∈ {1,2,3} 𝑔 Travel time between location j and j' 𝜖 Estimated threshold time at which SAR time at location j increases (due to secondary destructions) 𝑏′ Required SAR time in location j before the threshold time 𝑎 Additional SAR time in location j after the threshold time 𝛼 Time at which team i is available in the affected district 𝐿 Minimum acceptable involvement of a team in the operation at a location 𝑇 Duration of the route planning period 𝜏 Shift time 𝜌 Rest time Decision variables: 𝑠 Start time of team i at location j 𝑠′ Start time of SAR at location j 𝑓 Completion time of SAR at location j 𝑏Amount of SAR activity by team i at location j 𝑤 Amount of continuous SAR activity along the route by team i after visiting location jIf team i visits location j' after location j along their route Otherwise

  to (16) calculate the start time of the SAR at each demand location as the earliest time of operation by the allocated qualified teams.

  ; 𝑗 ≠ 1; 𝑗 ≠ 𝑚; 𝑐 ≥ 𝑑 (22) 𝑤 ≤ 𝑏 + 𝑀 2 -𝑣 -𝑟 ∀𝑖, 𝑗, 𝑗 ; 𝑗 ≠ 𝑗′;𝑗, 𝑗 ≠ 1; 𝑗, 𝑗 ≠ 𝑚; 𝑐 ≥ 𝑑 ; 𝑐 ≥ 𝑑 (23) 𝑤 ≤ 𝑏 + 𝑀 1 -𝑣 ∀𝑖, 𝑗′; 𝑗 = 1 ; 𝑗 ≠ 𝑗 ; 𝑐 ≥ 𝑑 (24) 𝑤 ≤ 𝑤 + 𝑏 + 𝑀 1 -𝑣 ∀𝑖, 𝑗, 𝑗 ; 𝑗 ≠ 𝑗 ; 𝑗, 𝑗 ≠ 1; 𝑗, 𝑗 ≠ 𝑚 ; 𝑐 ≥ 𝑑 ; 𝑐 ≥ 𝑑 (25) 𝑤 ≥ 𝑤 + 𝑏 -𝑀 1 + 𝑟 -𝑣 ∀𝑖, 𝑗, 𝑗′ ; 𝑗 ≠ 𝑗′ ; 𝑗, 𝑗 ≠ 1; 𝑗, 𝑗 ≠ 𝑚 ; 𝑐 ≥ 𝑑 ; 𝑐 ≥ 𝑑(26)

  𝑠 + 𝑏 + 𝑔 + 𝑟 𝜌 ≤ 𝑠 + 𝑀 1 -𝑣 ∀𝑖, 𝑗, 𝑗 ; 𝑗 ≠ 𝑗 ; 𝑗 ≠ 𝑚; 𝑐 ≥ 𝑑 ; 𝑐 ≥ 𝑑 + 𝑏 + 𝑀 1 -𝑞′ ∀𝑖, 𝑗; 𝑗 ≠ 1; 𝑗 ≠ 𝑚 (33) 𝑠 , 𝑠′ , 𝑓 , 𝑤 , 𝑏 ≥ 0 (34) 𝑥 , 𝑦 , 𝑞 , 𝑞′ , 𝑣 , 𝑟 ∈ {0,1}

Fig. 3 .

 3 Fig. 3. Tehran urban District 6 with 6 subdistricts (Google map, Ghajari et al., 2017)

Fig. 5 .

 5 Fig. 5. Coverage improvement vs. SAR response

Fig. 6 .

 6 Fig. 6. Coverage improvement glide path

Fig. 7 .

 7 Fig. 7. Performance of robust and deterministic approaches under uncertain parameter realizations at two perturbation levels

Fig. 9 .

 9 Fig. 9. Demand locations' delayed SAR completion with longer arrival and travel times

Fig. 10 .

 10 Fig. 10. Impact of variations in uncertain parameters on the objective function at three reliability levels

Fig. 11 .

 11 Fig. 11. Performance of robust and deterministic approaches under uncertain parameter realizations at two perturbation levels

Table 1 .

 1 Literature review summary

	Publication	Studied decisions	Uncertainty	Secondary risk	Identical resources	Lead time Obj. function	Solution method Real-life case	Fairness	Resource collaboration	Prioritized tasks
	Fiedrich et al., 2000	Scheduling

Table 3 .

 3 First-stage objective function's variations for different levels of BoU

	Uncertainty Budget		Certain	Scenario	Scenario	Base	Scenario	Scenario	Worst-case
	Scenarios		Scenario	1	2	Scenario	3	4	Scenario
	Demand type	1	0	3	4	5	6	6	6
	uncertainty	2	0	2	2	3	4	4	6
	budget	3	0	1	2	2	2	3	6
		1	0	2.4	3.2	4	4.8	5	5
	Travel time with	2	0	1.8	2.4	3	3.6	4.2	5
	other nodes	3	0	1.8	2.4	3	3.6	4.2	5
	uncertainty	4	0	1.8	2.4	3	3.6	4.2	5
	budget	5	0	2.4	3.2	4	4.8	5	5
		6	0	1.8	2.4	3	3.6	4.2	5
	Corresponding	1	43%	35%	33%	31%	30%	30%	30%
	coverage	2	9%	8%	8%	7%	7%	7%	6%
	response for demand type	3	40%	38%	36%	36%	36%	34%	29%
	Run time (Sec.)		99	47	156	101	60	502	200

Table 4 .

 4 Demand locations in the second stage

	Demand location	Demand type	Estimated population nominal value	Threshold time nominal value for SAR demand increase	SAR demand nominal value	Nominal value of SAR time increase after secondary destruction
	1	1	5000	2	6	6
		2	6000	14	23	0
		3	7000	14	28	0
	4	2	2000	14	2	0
	5	2	1800	14	7	0
	6	3	2300	14	9	0

Table 5 .

 5 Sensitivity of objective function to teams' availability lead time

					Arrival time of SAR teams
			Base scenario	Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5
		#1(type 1 capability)	1	1	1	1	3	6
	Teams	#2 (type 1 capability) #3 (type 2 capability) #4 (type 2 capability)	1 1 1	1 1 1	1 3 3	1 3 3	3 3 3	6 6 6
		#5 (type 3 capability)	1	3	3	6	3	6
		Objective function	7,848,718	7,865,518	7,889,038	7,914,238	7,919,038	8,004,128
		Increased workload vs. base scenario	-	-	-	-	6	11.5
		Delayed rescue start vs. base scenario	-	2	8	11	10	26.5
		Run time (Sec.)	6.7	9.4	13.4	7	48.1	52.9

  , five sub-measures aiming at the evaluation of the realized objective function value and imposed infeasibilities are calculated to define total performance:1) The new SAR completion times Min ∑ 𝑝 𝑓 based on the proposed visit sequences.

	Obj. fun. Change vs the	deterministic model	12% 16% 20% 24% 28% 32%	a) Variations in population		Obj. fun. Change vs the	deterministic model	0,4% 0,6% 0,8% 1,0% 1,2% 1,4% 1,6% 1,8%	b) Variations in SAR time
			8%								0,2%
				50%	80%	90%	99%			0,0%
					Robust model reliability level			50%	80% Robust model reliability level 90%	99%
			10% perturbation		20% perturbation		30% perturbation			10% perturbation	20% perturbation	30% perturbation
				c) Variations in travel time				36% d) Joint variations in all uncerain parameters
	Obj. fun. Change vs the	deterministic model	0,000% 0,010% 0,020% 0,030% 0,040%	50%	80%	90%	99%	Obj. func. change vs the	deterministic model	12% 16% 20% 24% 28% 32%
					Robust model reliability level			8%	50%	80%	90%	99%
			10% perturbation		20% perturbation		30% perturbation			Robust model reliability level
											10% perturbation	20% perturbation	30% perturbation

model

Atlas of Tehran metropolis, provided by Tehran municipality,

2 Loss data are altered to maintain confidentiality.

Acknowledgements

The authors would like to thank the Editor-in-Chief and autonomous reviewers for the valuable comments, which greatly improved the quality of this presentation. Also, the support of the Iranian Operations Research Society is highly acknowledged by the second author, as the board member of this society.

Inspired by the results of the TFES simulation for a nighttime earthquake, this section presents the results of MPFAM implementation in one of Tehran's municipal districts (District 6) for a 6period planning horizon and elaborates on the findings from solving SPCRRM in one of its subdistricts for the first period.

District 6 has an area of 21 km 2 and nearly 235,000 inhabitants and is located in the central part of the capital. It is one of the oldest districts in the city and has special characteristics compared to the rest 1 , which highlights the importance of SAR planning for successful loss prevention, not only at the district level but also beyond the district at the city level.

 An average of 9% of its population is above 64 years old, putting this district among the top three for the oldest population and among the most vulnerable districts, compared with an average of 6% in the Tehran megacity.  Up to 26% of its buildings are more than 30 years old.  It hosts the city's second-largest waste collection center and one of its five water treatment plants.  This district is the second frequent traffic destination among Tehran's 22 districts. It is located in a focal part with access to the city's main urban routes and is one of the key passages for delivering aid to the southern districts. Fig. 3 introduce the district with its subdistricts, and Table 2 presents the TFES 2 simulation results for a nighttime earthquake. The nominal values for SAR demand and travel times were extrapolated by SAR experts based on several criteria such as the number and severity of damaged buildings, average number of floors, average number of households in the buildings, and casualty estimates.

First-stage results

To test the results for the first stage, uncertainty budgets (𝛤 , 𝛤 ) are set as (5, 3, 2) and (4, 3, 3, 3, 4, 3) with 20% perturbation from nominal values. So, 5 out of 6 subdistricts have perturbations from their first-type demand nominal values, half have second-type demand different from nominal values, and 33% are not at their nominal third-type demand. Similarly, the travel time from the subdistrict 1 to four other subdistricts and the travel times from subdistrict 4 to three other subdistricts are perturbed.

Fig. 4 depicts the number of SAR units allocated to each subdistrict by the implementation of MPFRAM. The total number of allocated teams in each period and for each demand type represent the considered capacity of SAR resources and their distribution over the planning horizon. In the initial hours after the earthquake, only a small number of teams are available to carry out the SAR