Integrated System Design for Post-Disaster Management: Multi-Facility, Multi-Period, and Bi-Objective Optimization Approach

: Disaster management requires e ﬃ cient allocation of essential facilities with consideration of various objectives. During the response and recovery phase of disaster management (RRDM), various types of missions occur in multiple periods, and each of them needs di ﬀ erent support from facilities. In this study, a bi-objective mathematical model was derived to support multi-period RRDM by optimal allocation of required facilities such as drone stations, shelters, emergency medical facilities, and warehouses according to the mission life cycle. Connectivity between facilities was considered to ensure inter-facility complementarity. For e ﬃ cient derivation of Pareto solutions, a modi ﬁ ed epsilon-constraint algorithm for bi-objective optimization was developed. The algorithm was tested with a realistic disaster simulation scenario using HAZUS 4.0 as a demonstration of the bene ﬁ ts of the proposed approach. With the simulation experiments, the proposed approach was expected to provide e ﬃ cient operational plans and guidelines to decision makers for the bi-objective optimization problem in RRDM systems. In addition, the consideration of inter-facility connectivity can play an important role in the RRDM, especially for robustness and readiness.


Introduction
According to the Centre for Research on the Epidemiology of Disasters (CRED), disasters inflict high losses of both life and property.Disasters occurring between 2005 and 2014 caused damage to property estimated at 1.4 trillion USD and affected 1.7 billion people [1].More seriously, large-scale disasters will occur more frequently in the future due to unplanned urbanization, environmental degradation, pandemics, and climate change.Thus, disaster management systems that can perform their required missions more efficiently are absolutely necessary in order to recover more agilely from damage.
Comprehensive disaster management is commonly described in terms of four programmatic phases: preparation, response, recovery, and mitigation [2,3].The first phase of the disaster management cycle is preparation, in which, prior to a disaster's occurrence, plans are prepared for various types of disasters that could strike within the area of responsibility.The next phase is response, which takes place immediately after a disaster and involves a wide range of activities such as warning, evacuation, life saving, detection of further damage, distribution of emergency relief, etc. Next is the recovery phase, which focuses on longer-term responses to the disaster such as cleanup and rebuilding.The final phase is mitigation, which is activated to prevent similar disasters and their damages from reoccurring.
Over the course of all four phases of disaster management, different types of missions are accomplished.Figure 1 shows the various missions conducted sequentially as part of Over the course of all four phases of disaster management, different types of missions are accomplished.Figure 1 shows the various missions conducted sequentially as part of the overall disaster management process.Although every phase is vital in disaster management, the combined response and recovery phase of disaster management (RRDM) focuses crucially on post-disaster activities that can minimize the economic and human impacts of disasters and accelerate, as well as facilitate, a return to ordinary life.Especially in the RRDM, warnings to and evacuation of threatened populations are conducted immediately after a disaster occurs.Shelters are opened, and emergency rescue, medical care, and relief distributions are conducted simultaneously.Also, infrastructure restoration and rebuilding is performed in sequence.Due to the dual, simultaneous, and timesequential nature of the RRDM, it is essential to understand required missions that sequentially occur and to construct facility networks in a timely manner in order to efficiently handle the disaster and minimize the impact on society.The majority of missions in the RRDM require high flexibility for facility allocation and accessibility to disaster areas as access thereto generally is limited due to destruction of infrastructure such as road and telecommunication systems.Also, in disaster areas, there are many uncertainties, such as the risk of further disaster, the occurrence of unexpected events generally, and emergencies.One of the significant advances that has improved flexibility and accessibility is unmanned aerial vehicles (UAVs), also known as drones.UAVs are widely used for both commercial and civil missions; for example, they provide border security, coastguard, firefighting, emergency rescue, distribution infrastructure, monitoring, aerial photography, and communication relay services [5].Such advantages recommend UAVs for use in disaster management, especially for integration in RRDM scenarios to prevent the spread of damage and to facilitate quick recovery.For example, Vincenzi et al. [6] applied UAVs to the RRDM in search, rescue, and disaster relief roles and verified the benefits of UAVs in the management of historical disaster events (such as California wildfires, hurricanes in Mississippi, and Washington state landslides) generally and in real-time imagery, firefighting, and search and rescue functions specifically.
However, drones and small-sized UAVs have a fundamental dependence on limited energy sources (e.g., batteries or fuels), which shortens their endurance and, so too, their capacity to continuously execute missions [7].These limitations can be alleviated by locating battery recharge/replacement stations across the disaster area.In such cases, discharged drones will visit one of these stations, recharge, and return to their missions.The goal of this study was to optimize facility allocation for application of UAVs in the RRDM.The facility allocation problem includes various facility types (such as drone stations, The majority of missions in the RRDM require high flexibility for facility allocation and accessibility to disaster areas as access thereto generally is limited due to destruction of infrastructure such as road and telecommunication systems.Also, in disaster areas, there are many uncertainties, such as the risk of further disaster, the occurrence of unexpected events generally, and emergencies.One of the significant advances that has improved flexibility and accessibility is unmanned aerial vehicles (UAVs), also known as drones.UAVs are widely used for both commercial and civil missions; for example, they provide border security, coastguard, firefighting, emergency rescue, distribution infrastructure, monitoring, aerial photography, and communication relay services [5].Such advantages recommend UAVs for use in disaster management, especially for integration in RRDM scenarios to prevent the spread of damage and to facilitate quick recovery.For example, Vincenzi et al. [6] applied UAVs to the RRDM in search, rescue, and disaster relief roles and verified the benefits of UAVs in the management of historical disaster events (such as California wildfires, hurricanes in Mississippi, and Washington state landslides) generally and in real-time imagery, firefighting, and search and rescue functions specifically.
However, drones and small-sized UAVs have a fundamental dependence on limited energy sources (e.g., batteries or fuels), which shortens their endurance and, so too, their capacity to continuously execute missions [7].These limitations can be alleviated by locating battery recharge/replacement stations across the disaster area.In such cases, discharged drones will visit one of these stations, recharge, and return to their missions.The goal of this study was to optimize facility allocation for application of UAVs in the RRDM.The facility allocation problem includes various facility types (such as drone stations, shelters, emergency medical facilities, and warehouses) to support search, rescue, evacuation, and medical services in the RRDM.As candidate locations for facilities, existing social infrastructures including auditoriums, schools, hospitals, and warehouses will be Systems 2024, 12, 69 3 of 20 considered.In order to find optimal solutions, a bi-objective mathematical model along with a modified epsilon-constraint algorithm is proposed.

Literature Review
With the increasing importance of disaster management, various studies have been conducted to address the key issue of important facility allocation.In this context, the determination of shelter, emergency medical facilities, and warehouse locations has been intensively studied to minimize loss in terms of both human lives and economic damage and to facilitate recovery from disasters.The goal of the facility location problem is to determine the location of shelter facilities to which threatened people in dangerous places can be evacuated [8].Sherali et al. [9] developed nonlinear mixed-integer programming to minimize total congestion-related evacuation time under hurricane and flood conditions.For optimality, an exact algorithm was developed based on the Benders' decomposition method, and the proposed model and algorithm were verified with a Virginia Beach test problem.Kongsomsaksakul et al. [10] developed a bi-level mathematical model to design a shelter location and allocation network.In the upper-level problem, the authority determines the location of shelters with the objective of minimizing total evacuation time.In the lower level, evacuees choose the shelter and evacuation route.As the solution approached, a genetic algorithm was developed.Kar and Hodgson [11] suggested the use of a geographic information system (GIS) to score the suitability of emergency evacuation shelters in Southern Florida.Saadatseresht et al. [8] developed a multi-objective evacuation planning method for the determination of shelter utilization and travelling distance, along with an evolutionary algorithm to find Pareto-optimal solutions.
Locating emergency medical resources is another significant issue that has a direct association with saving human lives, and thus, many studies have been conducted to locate emergency medical resources, such as ambulances and medical facilities.Alsalloum and Rand [12] suggested goal programming as a means of locating emergency vehicles and thereby maximizing coverage of expected demand while minimizing vehicle number.Their developed model has been used to evaluate locations in Saudi Arabia.Rajagopalan et al. [13] considered dynamic redeployment of ambulances with fluctuating demand for time.Their proposed model pursued the minimum number of ambulances and their locations for each time cluster and validated their initial findings by simulation experiments.From the perspective of emergency medical facilities, Paluzzi [14] addressed a P-median-based heuristic location model for locating emergency service facilities in Carbondale.The goal of this model was to determine the optimal locations of emergency service facilities while minimizing the total aggregated distance between demand points and facilities.The effectiveness of the proposed model was validated by comparison with different types of location models.Jia et al. [15] proposed a general medical-service facility location model for large-scale emergencies such as earthquakes and terrorist attacks.The model can be cast as a covering model, a P-median model or a P-center model.The model was tested with illustrative examples in the Los Angeles area.Ko et al. [16] designed a system for determining the location, capacity, and capability of emergency medical facilities, considering multiple types of emergency diseases with corresponding survival rate functions.Fiedrich et al. [17] addressed the issue of emergency resource allocation after an earthquake disaster to minimize the total number of fatalities during the initial search and rescue period.A dynamic optimization model was developed to calculate the resource performance and efficiency for tasks related to the response.In addition, Fiedrich and Burghardt [18] introduced the use of an agent-based system for disaster management for more timely and enhanced data acquisition, information production, decision support, and action coordination.More recently, Sharma et al. [19] addressed location and allocation issues regarding temporary blood facilities for post-disaster periods with the objective of minimizing maximum distance between new facilities, permanent blood centers, and temporary blood centers.
The main goal of response missions in the RRDM is to provide in-time relief goods such as food, medicine, and hygiene products to affected areas [20,21].Therefore, warehouse location design must be optimized to enable efficient transport of relief goods from suppliers to disaster victims.Barzinpour and Esmaeili [22] developed a location model for relief chains in disaster management and addressed a multi-objective problem in order to minimize setup costs for facilities and transportation costs while maximizing cumulative coverage for a population.Mete and Zabinsky [23] considered the issue of medical supply location and distribution in disaster management, following a two-stage stochastic programming approach.In the first stage, the locations of warehouses and their inventory decisions were determined with the objective of minimizing the total cost of warehouse operation.In the second stage, transportation plans and demand satisfaction were addressed by minimizing the total transportation cost and unsatisfied demand.Moshref-Javadi and Lee [24] proposed a location-routing model for the minimization of total latency as a tool for the distribution of supplies to affected areas in post-disaster relief contexts.As a solution approach, they proposed a memetic algorithm and a recursive granular algorithm to find the most feasible solutions within a reasonable time.Rabiei et al. [25] focused on the assignment of volunteers in the post-disaster phase.Unmet triage needs and unsatisfied preferences of volunteers are formulated and considered as objective functions.Recently, Wang et al. [26] developed a two-stage distributionally robust optimization model for disaster relief logistics.The fixed costs of making contracts with suppliers, relief resource reservation, and pre-positioning relief resources are considered as systemic costs of relief logistics.To handle probability distribution of demand, a two-stage robust optimization model is developed.
In order to perform RRDM missions in a timely manner, UAV applications are strongly emphasized, due specifically to their high flexibility and accessibility.Vincenzi et al. [6] introduced the role of UAVs in RRDM efforts such as search, rescue, and disaster relief.Various UAV applications have been investigated: rescue missions in disasters [27], aerial imaging for post-disaster assessment [28], routing for continuous monitoring [29], communication systems for disaster recovery [30], and others.However, the issue of locating UAV battery recharge/replacement stations in disaster management has rarely been addressed.Kim and Morrison [31] tackled the UAV system design issue for search and patrol missions, deriving station location, UAV number per station, and UAV schedules simultaneously.Bor-Yaliniz et al. [32] suggested an approach for determining the efficient placement of aerial base stations.They considered the use of low-altitude UAVs to construct cellular networks for cases of unexpected events and developed and tested a mathematical formulation for the location of aerial base stations to maximize network revenue.However, all of the above studies took a short-term perspective unsuitable for long-duration RRDM.Also, there is uncertainty with respect to UAV flight distance.Kim et al. [33] addressed the stochastic facility location issue for UAVs with consideration of uncertain flight distances.Masroor et al. [34] suggested a multi-criterion optimization model for efficient deployment of UAVs for disaster management.The model pursues the maximization of users' connectivity while considering the minimum number of UAVs deployed.And recently, the effectiveness of drones in emergency situations was investigated by Zhang et al. [35].The use and effectiveness of drones for firefighting in high-rise buildings and logistics supports was validated based on the hybrid multi-criterion model of a fuzzy analytical hierarchy process, analytical network process, and decision-making trial and evaluation laboratory approaches.
Clearly, as seen in the review of the above-noted studies, location allocation issues regarding various types of essential facilities (e.g., shelters, emergency medical facilities, and warehouses) have been comprehensively addressed.However, most of the relevant studies are limited to the individual facility perspective and lack a holistic view of an integrated system with its inter-facility connectivity.For example, Yahyaei and Bozorgi-Amiri [36] emphasized the importance of integrated allocation of shelters and supply facilities.Furthermore, various RRDM missions (e.g., evacuation, search, transportation of injured people and relief goods) occur throughout the entire period, and thus it is necessary to allocate facilities that are close enough for inter-connectivity.Additionally to the issues of connectivity between facilities and multi-period missions, the location of UAV stations for replenishment and recharging must be considered simultaneously so as to facilitate the efficient use of UAVs in disaster management.
This paper addresses the optimization problem of facility allocation for multi-period missions in the RRDM in consideration of two objectives: minimization of total distance and minimization of the sum of setup and operation costs.First, in Section 3, we present the characteristics, constraints, and objectives of RRDM missions.In Section 4, we present a corresponding bi-objective mathematical model, and in Section 5, we propose a suitable solution approach.In Section 6, we verify the applicability of the proposed approach with simulation experiments and consider inter-facility connectivity as well.

Multi-Period Missions
Throughout the course of the RRDM, different types of missions must be pursued over multiple periods.Wallace and Baligh [37] classified post-disaster management missions into three temporal categories, within 0.1 year, from 0.1 year to 1 year, and after 1 year.For example, within 0.1 years of the occurrence of a disaster, the duties of damage assessment, epidemiological surveillance, and evacuation are performed, and search and rescue, treatment of injured, and relief logistics are performed sequentially thereafter.Similarly, Moe and Pathranarakul [38] suggested that disaster management phases (prediction warning, emergency relief, rehabilitation and reconstruction missions, in order) should occur during the RRDM.The study matches each phase with appropriate timing and activities for effective post-disaster management.
Balcik and Beamon [39] proposed a relief mission life cycle with different amounts of relative mission frequencies according to the phases of disaster management.As shown in Figure 2, it is clear that each phase has a different life cycle and required resources during the RRDM.The missions covered in the present study are limited to search, evacuation, emergency patient transfer, and relief distribution.It is assumed that the RRDM proceeds over T periods.Each mission has a different occurrence period and life cycle.A mission consists of a set of tasks that can be defined by location, period, and demand.A task in a specific location can be defined over several time periods and can have different amounts of demand for each time period.
Amiri [36] emphasized the importance of integrated allocation of shelters and supply facilities.Furthermore, various RRDM missions (e.g., evacuation, search, transportation of injured people and relief goods) occur throughout the entire period, and thus it is necessary to allocate facilities that are close enough for inter-connectivity.Additionally to the issues of connectivity between facilities and multi-period missions, the location of UAV stations for replenishment and recharging must be considered simultaneously so as to facilitate the efficient use of UAVs in disaster management.
This paper addresses the optimization problem of facility allocation for multi-period missions in the RRDM in consideration of two objectives: minimization of total distance and minimization of the sum of setup and operation costs.First, in Section 3, we present the characteristics, constraints, and objectives of RRDM missions.In Section 4, we present a corresponding bi-objective mathematical model, and in Section 5, we propose a suitable solution approach.In Section 6, we verify the applicability of the proposed approach with simulation experiments and consider inter-facility connectivity as well.

Multi-Period Missions
Throughout the course of the RRDM, different types of missions must be pursued over multiple periods.Wallace and Baligh [37] classified post-disaster management missions into three temporal categories, within 0.1 year, from 0.1 year to 1 year, and after 1 year.For example, within 0.1 years of the occurrence of a disaster, the duties of damage assessment, epidemiological surveillance, and evacuation are performed, and search and rescue, treatment of injured, and relief logistics are performed sequentially thereafter.Similarly, Moe and Pathranarakul [38] suggested that disaster management phases (prediction warning, emergency relief, rehabilitation and reconstruction missions, in order) should occur during the RRDM.The study matches each phase with appropriate timing and activities for effective post-disaster management.
Balcik and Beamon [39] proposed a relief mission life cycle with different amounts of relative mission frequencies according to the phases of disaster management.As shown in Figure 2, it is clear that each phase has a different life cycle and required resources during the RRDM.The missions covered in the present study are limited to search, evacuation, emergency patient transfer, and relief distribution.It is assumed that the RRDM proceeds over T periods.Each mission has a different occurrence period and life cycle.A mission consists of a set of tasks that can be defined by location, period, and demand.A task in a specific location can be defined over several time periods and can have different amounts of demand for each time period.

Facility Candidates
In a disaster situation, undamaged or damaged but usable civil facilities are preferentially considered as facility candidates.Coppola [40] suggested a number of candidate  1, both schools and warehouses are suitable candidates for evacuation and relief distribution.In addition, in order to serve search operations such as surveillance, facilities with large spaces such as schools and parks can be considered as drone stations.In this study, existing and usable civil facilities in the disaster area were considered as priority candidate sites.One candidate can be used for various purposes, as noted above, and the usability of a candidate site is known in advance.The available capacity and operation costs of each candidate site also are given in advance.

Objectives
The goal of this study was to achieve two objectives: minimization of total transportation distance (for serving of all tasks) and minimization of the sum of the setup and operating costs of facilities.Minimizing the total distance and cost is critical and standard in disaster management (see [14,19,22,23,26], etc.).Furthermore, the two objective functions are in inverse relation, providing a suitable basis for analyzing trade-offs among Pareto solutions.The choice of the objective function can vary depending on the nature of the problem and the decision maker's subjectivity.Additionally, it is worth considering that, if computational power allows, the development of a multi-objective mathematical model with more than two objective functions is also feasible.Minimization of evacuation time, resources, and maximizing coverage can be an appropriate candidate for more objective functions, as seen in the literature [9,10,12,17,22].
The best allocation of facilities may change as the RRDM proceeds.In this light, Figure 3 provides an illustrative example of a feasible solution to a problem instance with 30 tasks, 5 candidate sites, and 10 time periods.As can be seen, the entire RRDM entails 30 tasks that occur within the 10 time periods of the mission life cycle.There are 5 facility candidates, namely a school, hospital, two warehouses, and a drone station.The different types of tasks (i.e., search, evacuation, medical treatment, and relief distribution) are distinguished by different colors.

Facility Candidates
In a disaster situation, undamaged or damaged but usable civil facilities are preferentially considered as facility candidates.Coppola [40] suggested a number of candidate shelter and emergency medical facilities.Nurre et al. [41] recommended well-known, central entities in the city as candidates for relief goods warehouses.Table 1 summarizes possible facility candidates.A single candidate can be used for a variety of purposes.For example, in Table 1, both schools and warehouses are suitable candidates for evacuation and relief distribution.In addition, in order to serve search operations such as surveillance, facilities with large spaces such as schools and parks can be considered as drone stations.In this study, existing and usable civil facilities in the disaster area were considered as priority candidate sites.One candidate can be used for various purposes, as noted above, and the usability of a candidate site is known in advance.The available capacity and operation costs of each candidate site also are given in advance.

Objectives
The goal of this study was to achieve two objectives: minimization of total transportation distance (for serving of all tasks) and minimization of the sum of the setup and operating costs of facilities.Minimizing the total distance and cost is critical and standard in disaster management (see [14,19,22,23,26], etc.).Furthermore, the two objective functions are in inverse relation, providing a suitable basis for analyzing trade-offs among Pareto solutions.The choice of the objective function can vary depending on the nature of the problem and the decision maker's subjectivity.Additionally, it is worth considering that, if computational power allows, the development of a multi-objective mathematical model with more than two objective functions is also feasible.Minimization of evacuation time, resources, and maximizing coverage can be an appropriate candidate for more objective functions, as seen in the literature [9,10,12,17,22].
The best allocation of facilities may change as the RRDM proceeds.In this light, Figure 3 provides an illustrative example of a feasible solution to a problem instance with 30 tasks, 5 candidate sites, and 10 time periods.As can be seen, the entire RRDM entails 30 tasks that occur within the 10 time periods of the mission life cycle.There are 5 facility candidates, namely a school, hospital, two warehouses, and a drone station.The different types of tasks (i.e., search, evacuation, medical treatment, and relief distribution) are distinguished by different colors.

Mathematical Formulation
We present a mathematical model of bi-objective optimization of facility allocation for post-disaster missions utilizing drones.The indices, sets, parameters, and decision variables are defined in the following section.

Formulation
Minimize Systems 2024, 12, 69 The proposed mathematical model serves the two objectives noted above: minimization of total distance and minimization of the sum of setup and operation costs.Equations ( 1) and ( 2) show the bi-objective functions of the proposed study, F 1 and F 2 .This study deals with four types of mission, each of which includes a transportation context.Equation ( 1) is used to minimize the total transportation distance (F 1 ) while serving all tasks.Equation ( 2) is formulated to minimize the total setup cost and operation cost of facilities (F 2 ).Equation ( 3) represents the single-assignment constraint.Via Equation ( 3), each task should be served by a single facility.The role of Equation ( 4) is to match the task with the facility type.The variable m i signifies the task type, thereby having one of the values among {S, E, M, R}, and n j is defined as a subset of {S, E, M, R} (refer to Appendices A and B).If candidate facility j ∈ J is not able to serve task i with m i , which means that the value of m i is not an element of set n j , Equation (4) makes the assignment link impossible.Using Equation (5), the assignment link between task i ∈ I and candidate j ∈ J becomes valid only for the time period t ∈ u i .Equation (6) shows the capacity restriction of each facility.Equation ( 7) is formulated to reduce the search effort of locating and operating facility candidate j ∈ J during time period t ∈ T. If there is no demand for task type n j during t, the candidate facility j ∈ J will not be opened during t ∈ T for type n j .This is a useful constraint because each mission in this study has a different life cycle.During t ∈ T, the facility candidate j ∈ J can be opened for a single purpose out of n j via Equation (8).The occurrence of facility setup is confirmed via Equation (9).In this study, it is assumed that setup occurs only once for each facility j ∈ J and type k ∈ n j .Equations ( 10) and (11) are formulated to limit the service ranges of drone stations and emergency medical facilities.These ranges are limited by the drone station service radius (R) and the maximum allowed travel time of emergency patients (S max ).Decision variables are denoted in Equations ( 12)-( 14).

Estimation of Facility Service Range
The service ranges of emergency medical facilities are limited by the degree of injury of victims, and UAVs' service time is limited by their energy source.To consider such limitations, the service ranges of drone stations and emergency medical facilities are limited by drone station service radius (R) and maximum allowed travel time of emergency patients (S max ), as just noted.R and S max will guarantee feasible and realistic facility allocation that can serve every search and patient transfer task in the system.For estimation of R, let V be the set of heterogeneous drones that are used in RRDM.FQ v and TS v stand for the flight capacity (minutes) and travel speed (meters per minute) of drone v ∈ V. ∼ P i means the service time (minutes) of a search task, and it takes values between [P i − Pi and P i + Pi ].The notations P i and Pi denote the nominal value and the maximum deviation from the nominal value, respectively.Equation ( 15) describes the service range of a drone station, and minimum estimation of R can be obtained when Now, using Equation (15), Equation ( 10) can be reformulated to Equation ( 16) as follows: For estimation of S max , there are many guidelines regarding emergency patient transfer.Typically, many patients in a disaster area suffer traumatic injury.The terminology "golden hour" represents the first 60 min after traumatic injury has been incurred.The idea is that trauma patients have better outcomes if they are provided definitive care within 60 min of their injuries [42].In this context, Song et al. [43] theoretically estimated the service ranges of UAV service stations.Also, legislation can be used to facilitate the estimation of S max .For example, in June 2005, New Hampshire law was amended to enable alternative health care to participate in inter-facility transfer if the availability of conventional providers exceeds 30 min according to the National Highway Traffic Safety Administration [44].

Pareto Optimality
This section reviews the main concepts and definitions pertaining to the bi-objective optimization problem (BOP).In general, BOP minimization is formulated as Subject to x ∈ D; where x = (x 1 , x 2 , . .., x n ) is the vector representing the decision variables, and D is the set of feasible solutions.A BOP usually does not have a unique optimal solution but rather a set of solutions known as the Pareto-optimal set.Each Pareto-optimal solution represents a compromise between different objectives, and the components of the corresponding vector of objectives cannot all be simultaneously improved.In the minimization setting, Pareto dominance and Pareto optimality are defined as follows.
Definition 1. (Pareto dominance).A given vector u ∈ D dominates a vector v ∈ D in the Pareto sense, if and only if u is partially less than v (u < v), i.e., Definition 2. (Pareto-optimal solution).A solution u ∈ D is a Pareto-optimal solution, if and only if there is no v ∈ D such that v dominates u.Pareto-optimal solutions are also called efficient non-dominated solutions.
Definition 3. (Pareto-optimal set).The Pareto-optimal set or the efficient set is defined as P = {x ∈ D: x is a Pareto-optimal solution in D}.

Modified Epsilon-Constraint Algorithm
Chankong and Haimes [45] proved that for the general multi-objective optimization problem, exact Pareto-optimal solutions can be found by the epsilon-constraint method.This methodology is still widely used and has been improved for various multi-objective optimization problems [46][47][48].In the present study, a modified epsilon-constraint method was introduced to find every Pareto-optimal solution of the proposed bi-objective mathematical problem.
The epsilon-constraint method solves a constrained single-objective problem set that is obtained by choosing one objective as the only objective to optimize, incorporating epsilon constraints for the remaining objectives, as shown in Equation ( 17), and assigning different values to the components of the epsilon vectors.In this study, an exact Pareto set was obtained by solving the following epsilon-constraint model (M EC ) problem iteratively while changing the value of ε.
F 1 and F 2 are depicted in Equations ( 1) and ( 2), respectively.D is the feasible region defined by constraints (3) to (14).The value of F 2 has a strict upper bound.The upper bound of F 2 , UB(F 2 ), is obtained by solving max x∈D F 2 (x).
Let t and ε t be the iteration index and the value of ε at the tth iteration of M EC , respectively.The proposed epsilon-constraint method starts with ε 0 = UB(F 2 ) and iteratively decreases ε using the F 2 value of the optimum of the previous single-objective run.The solution procedure continues as long as the resulting problem M EC is feasible.Proposition 1 demonstrates the Pareto relationship between F 1 and F 2 .Proposition 2 guarantees that the proposed iterative procedure will search all Pareto solutions.The optimal solution for (M EC ,ε t ) is denoted by opt(M EC ,ε t ).The values of F 1 and F 2 in opt(M EC ,ε t ) are denoted by F 1 *(M EC ,ε t ) and F 2 (M EC ,ε t ), respectively.Note that F 2 (M EC ,ε t ) = ε t+1 .Let D(M EC , ε t ) be the feasible region defined by (M EC , ε t ).
Proposition 1. F 1 *(M EC ,ε t ) monotonically increases as t increases.And there exists a Pareto relationship between F 1 and F 2 .
Proof.The monotonic increase in F 1 *(M EC ,ε t ) with respect to t comes from the differences in the feasible region.As t increases, ε t (=F 2 *(M EC ,ε t−1 )) decreases, thus limiting the feasible region.Because As a result, Pareto optimality exists between F 1 and F 2 .
Proof.Let vector b be the opt(M EC ,ε t−1 ).Then, dominate the optimal solution vector b (M EC ,ε t−1 ) in the Pareto sense.
Algorithm 1 summarizes the overall solution procedure of the proposed epsilonconstraint method.In this algorithm, set S stands for the set of Pareto solutions, and s ∈ S is a single Pareto solution.Steps 6 and 7 in Algorithm 1 are required to handle multioptimality on F 1 .In the multiple optimal solutions on F 1 , the same F 1 value can be achieved with different F 2 values.However, only one, that with the lowest F 2 value, is selected for the Pareto solution among the multiple optimal solutions on F 1 .As Algorithm 1 proceeds, the value of F 2 (M EC ,ε t ) gradually diminishes, ultimately causing the right-hand side of Equation ( 17) to decrease, resulting in the emergence of an infeasible case.

Illustrative Case Study
The simulation data on the earthquake in San Diego were generated using Hazus 4.0 [49], which is a software package with a geographic information system (GIS) called ArcGIS.The earthquake module in Hazus can evaluate a wide range of losses resulting from earthquake scenarios and can be used to transform the information of ArcGIS into quantitative estimates of damage and loss.The damaged area and number of casualties and shelters are calculated by two modules: the direct damage module and the direct loss module.
To generate the damage data on the earthquake, we made a hazard scenario of magnitude 9.0 for San Diego by setting an epicenter.The detailed parameters are provided in Table 2. Based on the demographic and building information, four types of demand (Search, Evacuation, Medical treatment, and Relief Logistics) were calculated with the Hazus software.Also, five facility types (Auditorium, Hospital, Park, School, and Warehouse) were selected from the actual data on the facilities.Detailed facility and task information can be found in Appendices A and B, respectively.All of the demand and facility locations for the earthquake scenario are indicated in Figure 4 below.
The results of the experiment using the proposed algorithm are summarized in Table 3.The algorithm is coded in C++ with IBM ILOG CPLEX 12.4 and runs on an Intel Xeon E5-1620 3.6 GHz processor with 16 GB RAM.The service radius of a drone station (R) and the maximum allowed transportation time for emergency patients (S max ) are specified as 7 km and 20 min, respectively.Based on these Pareto solutions, a Pareto curve is drawn in Figure 5.This result shows that the proposed algorithm can find Pareto solutions effectively within a reasonable time.
Based on the demographic and building information, four types of demand (Search, Evacuation, Medical treatment, and Relief Logistics) were calculated with the Hazus software.Also, five facility types (Auditorium, Hospital, Park, School, and Warehouse) were selected from the actual data on the facilities.Detailed facility and task information can be found in Appendices A and B, respectively.All of the demand and facility locations for the earthquake scenario are indicated in Figure 4 below.The results of the experiment using the proposed algorithm are summarized in Table 3.The algorithm is coded in C++ with IBM ILOG CPLEX 12.4 and runs on an Intel Xeon E5-1620 3.6 GHz processor with 16 GB RAM.The service radius of a drone station (R) and the maximum allowed transportation time for emergency patients (Smax) are specified as 7 km and 20 min, respectively.Based on these Pareto solutions, a Pareto curve is drawn in Figure 5.This result shows that the proposed algorithm can find Pareto solutions effectively within a reasonable time.For example, Figure 6 illustrates the results of facility allocations of a Pareto solution wherein F1 and F2 are 1030.67 and 78,684, respectively.According to the solution, Christian Unified Schools supports evacuation and relief distribution by changing setups three times, while Liberty Charter School supports the same missions by changing setups twice.In addition, John F Kennedy Park is allocated as a drone station only for the three beginning periods even though it is available for all time periods, due to the life cycle of the For example, Figure 6 illustrates the results of facility allocations of a Pareto solution wherein F 1 and F 2 are 1030.67 and 78,684, respectively.According to the solution, Christian Unified Schools supports evacuation and relief distribution by changing setups three times, while Liberty Charter School supports the same missions by changing setups twice.In addition, John F Kennedy Park is allocated as a drone station only for the three beginning periods even though it is available for all time periods, due to the life cycle of the search mission and the relatively shorter distances from searching locations.Throughout the case study, it was demonstrated that the proposed mathematical model can utilize existing facilities for various purposes and change their roles during a long-duration RRDM in order to efficiently achieve missions.Also, the Pareto set was derived in a reasonable time to support decision making, particularly considering the trade-offs between objectives.

Results of Large-Sized Experiments
To evaluate the applicability for more realistic problems of greater complexity, we designed large-sized experiments by increasing the numbers of tasks and facilities.With the various sizes of task sets and candidate site sets, Pareto solutions were generated, and the ranges of F1 and F2 values, the number of Pareto solutions and the computation times were confirmed for each example.T was set to 10 for every example, and the numbers of the four tasks were set as 25, 50, 75, and 100, respectively.The number of facilities was determined to be 10, 20, 30, 40, and 50 for each experiment and quartered into four types.The demands of tasks and the capacities of facilities were generated from the range of Hazus simulation results.The density of tasks in the experiment was defined as ∑  ∈ .The environment and computer configurations for the experiments were the same as described in the precious section.
Table 4 summarizes the results of the large-sized implementations.Because an increase in |J| causes the expansion of the solution space, the minimum values of F1 and F2 were consistently decreased as |J| increased for each case of |I| = (25,50,75,100).In addition, the larger solution space generated a greater number of Pareto solutions for the computation cost.The total computation time and computation time per Pareto solution increased continuously as the problem size became larger.In the largest case, with (|I|, |J|) = (100, 50), the solution procedure required 6125 s (>1.7 h), which is reasonable in terms of RRDM system design.

Results of Large-Sized Experiments
To evaluate the applicability for more realistic problems of greater complexity, we designed large-sized experiments by increasing the numbers of tasks and facilities.With the various sizes of task sets and candidate site sets, Pareto solutions were generated, and the ranges of F 1 and F 2 values, the number of Pareto solutions and the computation times were confirmed for each example.T was set to 10 for every example, and the numbers of the four tasks were set as 25, 50, 75, and 100, respectively.The number of facilities was determined to be 10, 20, 30, 40, and 50 for each experiment and quartered into four types.The demands of tasks and the capacities of facilities were generated from the range of Hazus simulation results.The density of tasks in the experiment was defined as ∑ i∈I u i .The environment and computer configurations for the experiments were the same as described in the precious section.
Table 4 summarizes the results of the large-sized implementations.Because an increase in |J| causes the expansion of the solution space, the minimum values of F 1 and F 2 were consistently decreased as |J| increased for each case of |I| = (25, 50, 75, 100).In addition, the larger solution space generated a greater number of Pareto solutions for the computation cost.The total computation time and computation time per Pareto solution increased continuously as the problem size became larger.In the largest case, with (|I|, |J|) = (100, 50), the solution procedure required 6125 s (>1.7 h), which is reasonable in terms of RRDM system design.

Connectivity for Complementary Response
In a previous study [50], two facilities were defined as 'connected' if the distance between them was less than a given value.The concept of connectivity is critical to the provision of seamless RRDM operations when certain important facilities are not available.In addition, nearby facilities of the same type can share their emergency supplies, enabling more flexible responses.In order to improve preparedness by considering connectivity, we propose the following additional constraint where the new parameter D jj represents the distance between facilities j and j in J, V j is defined as a set of facilities within a given range MR for facility j ( V j = j D jj ≤ MR ), and the new variable n determines the number of facilities within the range.The new Equation ( 18) was applied to the case study discussed in Section 6.1.By adding the constraint with a range (MR = 5 km) and n = 2, a new Pareto solution was derived with the same epsilon (ε = 80,000), as shown in Figure 7.By adding this additional constraint, the computation time was slightly increased (by 0.23 s in CPU time) per solution under the same epsilon.To graphically investigate the connectivity, two minimum F 1 solutions of time period 1 (with/without connectivity constraint) were plotted, as depicted in Figures 8 and 9.As shown in Figure 8, some facilities in the solution were distant and without connectivity; in fact, due to their long distances, they would not be able to provide complementary responses in continuous operations when uncertainties (e.g., beyondprediction demand, additional damage to facilities, etc.) arose.But with constraint (18), more than one facility can easily collaborate in disaster scenarios, as shown in Figure 9, with the benefit of connectivity.However, in the solution shown in Figure 9, the facilities conceded the distance to the task nodes in order to ensure connectivity.As a consequence, the total distance to the task nodes increased slightly (to 114.37 km) compared with the original solution.This result exemplifies the trade-offs between two operational issues, connectivity and service quality.Consideration of connectivity can improve the preparedness and robustness of solutions by connecting facilities of the same type within a given distance.However, it may also cause inefficiency, in terms of service and response quality, by locating facilities in long-distance areas.Thus, the results showed that decision makers can use the proposed model to conduct a quantitative analysis of potential increases in objective value to improve connectivity.

Concluding Remarks
Disasters impart great impacts on human life and civil infrastructure.To minimize such harms and promote a return to ordinary life, various missions are conducted during the RRDM.The key aspects of the RRDM are concurrency and a long-term perspective.Various missions occur simultaneously, and each of them has a different mission life cycle.Furthermore, some missions need high flexibility and accessibility.As a result, RRDM system design requires a flexible, integrated, and long-term approach.In this study, an integrated system design was pursued for the entire RRDM.Locations and mission allocations of drone stations, shelters, emergency medical facilities, and warehouses were determined via a mathematical optimization approach.By considering the life cycle of each mission and the usability of existing social infrastructure facilities, the proposed system provides real-world applicability as well as academic value.
To derive a complete Pareto set with the proposed bi-objective mathematical model, a modified epsilon-constraint algorithm was developed.The mathematical model and algorithm proposed in Sections 4 and 5, respectively, were tested using simulation experiments based on HAZUS 4.0, and the applicability for large-sized problems was demonstrated via large-size implementation.With the simulation experiments, the proposed approach is expected to provide efficient operational plans and guidelines to decision makers for the bi-objective optimization problem in RRDM systems.Additionally, in terms of the inter-facility connectivity as a constraint, the trade-offs between connectivity and service quality issue were examined to provide complementary and flexible responses.The consideration of inter-facility connectivity can play an important role in the RRDM, especially for robustness and preparedness.
Possible future research may include operational levels such as detailed routing of vehicles (e.g., trucks).In addition, various scenarios in earthquake-prone countries such as RISK-UE can be considered for additional verifications.Lastly, for optimally efficient RRDM, a set of heterogeneous vehicles (e.g., trucks, drones, and mobile robots) would be worthy of further investigation.

Figure 3 .
Figure 3. Illustrative example of changes in tasks for facilities (30 tasks and 5 candidate sites over 10 time periods).

Figure 4 .
Figure 4. Illustration of facilities and tasks for earthquake scenario in San Diego.In the right-side figure, the different circle sizes represent the demand levels of tasks.

Figure 4 .
Figure 4. Illustration of facilities and tasks for earthquake scenario in San Diego.In the right-side figure, the different circle sizes represent the demand levels of tasks.

Figure 5 .
Figure 5. Pareto curve for earthquake scenario in San Diego.

Figure 5 .
Figure 5. Pareto curve for earthquake scenario in San Diego.

Figure 8 .
Figure 8. Illustration of facility types at time 1 without consideration of connectivity.

Figure 8 .
Figure 8. Illustration of facility types at time 1 without consideration of connectivity.Figure 8. Illustration of facility types at time 1 without consideration of connectivity.

Figure 8 .
Figure 8. Illustration of facility types at time 1 without consideration of connectivity.Figure 8. Illustration of facility types at time 1 without consideration of connectivity.

Systems 2024 , 21 Figure 9 .
Figure 9. Illustration of facility types at time 1 with consideration of connectivity.
[41]emergency medical facilities.Nurre et al.[41]recommended well-known, central entities in the city as candidates for relief goods warehouses.Table1summarizes possible facility candidates.A single candidate can be used for a variety of purposes.For example, in Table shelter

Table 1 .
Examples of facility candidates for each type.

Table 1 .
Examples of facility candidates for each type.
∈ u i .It is equal to 0 for t / ∈ u i ; Q j,k: Capacity of facility j when it is used as type k ∈ n j .It is equal to 0 for k / ∈ n j ; C j,k ij : Travelling time between task i and candidate facility j.T ij is not proportional to D ij ; A i t : Demand quantity of task i, t t : Operation cost of facility j of type k ∈ n j for period t ∈ T. It is equal to 0 for k / ∈ n j ; S j,k : Setup cost of candidate j for type k ∈ n j facility; R : Service radius of drone station; S max : Maximum allowed transportation time of emergency patients; M : A large positive real number; O j,k : Binary setup decision variable.It is equal to 1 if candidate facility j is used for type k ∈ n j facility; Y j,k t : Binary location decision variable.It is equal to 1 if candidate facility j of type k ∈ n j is opened for period t; X i,j t : Binary assignment decision variable.It is equal to 1 if task i of k = m i , t ∈ U i is served by facility j;

Table 2 .
Parameters for case study.

Table 3 .
Result summary of earthquake scenario generated from Hazus.

Table 3 .
Result summary of earthquake scenario generated from Hazus.

i∈I u i ) |J| Value of F 1 (km) Value of F 2 (1000 USD) Number of PS
Systems 2024, 12, x FOR PEER REVIEW 13 of 21

Table 4 .
Results of large-sized experiments.

Table 4 .
Results of large-sized experiments.