Optimization of Humanitarian Aid Distribution in Case of an Earthquake and Tsunami in the City of Iquique, Chile

: In this paper, we introduce, model, and solve a clustered resource allocation and routing problem for humanitarian aid distribution in the event of an earthquake and subsequent tsunami. First, for the preparedness stage, we build a set of clusters to identify, classify, sort, focus, and prioritize the aid distribution. The clusters are built with k -means method and a modiﬁed version of the capacitated p -median model. Each cluster has a set of beneﬁciaries and candidate delivery aid points. Second, vehicle routes are strategically determined to visit the clusters for the response stage. A mixed integer linear programming model is presented to determine efﬁcient vehicle routes, minimizing the aid distribution times. A vulnerability index is added to our model to prioritize aid distribution. A case study is solved for the city of Iquique, Chile.


Introduction
The disasters have increased 41% between 2005 and 2015 worldwide [1]. Particularly, Chile is exposed to several types of catastrophic events, such as earthquakes, tsunamis, volcanic eruptions, wildfires, and technological events [2,3]. The uncertainty of these events has generated an important interest in humanitarian aid management research to better prepare and respond and, consequently, save lives and reduce the people suffering in case of those catastrophic events.
In this context, humanitarian aid logistics is defined as the efficient and effective logistics assistance in response to humanitarian crises and disasters [4]. The objective is to provide rescue, medical assistance, water, food, and shelter [5]. The three main associated principles are: humanity, neutrality, and impartiality [6,7].
This article studies a clustered resource allocation and routing problem in the humanitarian aid distribution to reduce the distribution times to the final beneficiaries in case of an earthquake and tsunami. We define a Beneficiary Point (BP) as a point (centroid) that concentrates the potential demand of a group of people within a census block. Also, we define an Aid Delivery Point (ADP) as a point that corresponds to a public facility such as: school, stadium, university, health center, municipal organization, etc. The ADPs are in charge of receiving and organizing the aid dispatch. Next, BPs and ADPs are grouped into disjoint clusters using k-means and a capacitated p-median model as clustering methods. Note that the strategic clustering allows estimating, sorting, focusing, and prioritizing the resource allocation [8,9] and identifying some facilities for aid distribution to the potentially affected people, avoiding clutter and excessive delays to reduce people suffering [10]. Moreover, realistic size instances (like the one we address in this study) are best dealt with clustering, downsizing the network [11,12].
In our aid distribution scheme, the humanitarian aid is delivered by vehicles to ADPs. Next, the BPs travel by themselves (using vehicles, walking, bike, etc.) to their allocated nearby ADP. The aid distribution is carried out by vehicles, starting and ending their routes in a distribution center (depot) where the primary aid arrives (e.g., airport, seaport, warehouse). The objective is to minimize the total travel times for humanitarian aid in an earthquake and tsunami. Figure 1 shows a scheme of the problem. The ADPs are the circle nodes, while the BPs are the square nodes. The straight arrows indicate the routes of the vehicles, and dotted arrows indicate the allocations from BP to their closest ADP. The dotted circles are the clusters. The first vehicle distributes aid from a central depot (pentagon node) to ADP 4 in the Cluster 2; hence, this node serves the BP 10. Next, the vehicle serves ADPs 5 and 6 in the Cluster 3 (serving BPs 11 and 12, respectively) and returns to the depot. The second vehicle serves ADPs 3 and 1 in the Cluster 2 and 1, respectively, serving the BPs 7, 8, and 9. Notice that it is not necessary to visit all ADPs (e.g., ADP 2). Moreover, the humanitarian aid could be distributed to more than one ADP within each cluster (e.g., in Figure 1, Cluster 2 is visited by more than one vehicle). A vehicle enters and leaves a cluster just once, and an ADP within a cluster could be visited by more than one vehicle if it is convenient to reduce distribution times or if the capacity of a single vehicle is not enough to distribute all aid of a cluster. It is important to notice that the same network structure could be helpful for other applications such as vaccine distribution, insular waste collection, food, and water distribution, etc. Moreover, other disasters, such as floods, alluviums, wildfires, tornadoes, may use a similar structure for aid distribution.
Our main contribution is to propose a quantitative methodology for organizing and distributing humanitarian aid in case of an earthquake and subsequent tsunami. First, we state strategic clusters to analyze and prepare the area. Second, we determine efficient vehicle routes for humanitarian aid distribution. It should be noticed that the structure of this problem is a generalization of the Generalized Median Tour Problem (GMTP), introduced by [13], since they focus their study on a single tour visiting a set of clusters. Third, we model the problem using a mixed-integer linear programming (MILP) model. Finally, we solve a case study in the city of Iquique, Chile, using several solution strategies.
The rest of the article is organized as follows. Section 2 presents a brief literature review containing up-to-date related articles to our problem. Section 3 shows the methodology for cluster building and the routing problem. Section 4 presents a set of experiments on a case study of Iquique, Chile. Finally, Section 5 shows conclusions and future work.

Literature Review
In the following, we divide the literature review section into two sub-sections. The Section 2.1 is devoted to the related humanitarian aid routing problems. The Section 2.2 is focused on the network design problems that have a similar structure to our problem.

Humanitarian Aid Routing
The humanitarian supply chain aims to supply, quickly as possible, the emergency supplies to the affected people in a disaster [14]. In this context, the last mile distribution in humanitarian aid consists of determining efficient distribution plans from central warehouses to affected people [15]. The main activities are choosing the aid delivery points (shelters, delivery points, etc.), allocating the beneficiaries to aid delivery points, and determining vehicle routes.
Most articles study the aid distribution problem in the context of humanitarian logistics as a Vehicle Routing Problem (VRP) or Traveling Salesman Problem (TSP). The articles [16][17][18][19][20] provide a complete reviews about routing problems in humanitarian logistics. Usually, the authors consider the optimization of one or more objectives such as: minimize the response time, minimize the maximum travel time, maximize the travel reliability, minimize costs, and minimize the unsatisfied demand, among others [6,14,21,22]. Moreover, article [6] indicates that different additional settings are considered for the goods distribution: stochastic demand, multi-commodities, stochastic offer, etc.
Few authors consider the use of clusters in humanitarian logistics aid. Ref. [12] presents a hierarchical clustering and routing procedure for coordinating vehicle routes in post-disaster activities. They consider a multi-level problem, where the first level is to determine the total potential demand in small clusters. Then, the routes are generated iteratively. Ref. [23] prioritizes the most vulnerable population for humanitarian distribution. The demand is grouped into clusters using a fuzzy clustering method. Finally, they develop a multi-objective model to define priority for aid distribution. Ref. [24] provides a k-means clusterization method to locate distribution centers and assign distribution areas in case of a disaster. Ref. [25] states an inter-modal routing model for medical supply distribution. They cluster the potentially affected population using a fuzzy clustering method. Within each cluster, a distribution center is stated for helicopter landing. From this point, a vehicle distributes the supplies to each node in the cluster. Ref. [26] solves a routing problem using the strategy cluster first-route second. The authors use two clustering methods: Clustering Local Observation, and k-means. The vehicle routing model provides routes to serve a set of demand nodes.

Clustered and Generalized Vehicle Routing Problems
The vehicle routing problems that consider the nodes are grouped into clusters are usually known as generalized or clustered vehicle routing problems. The generalized vehicle routing problem (GVRP), proposed first time by [27] consists of determining routes for a set of vehicles that must visit a set of nodes grouped within clusters. Each cluster is visited once by a single vehicle. Thus, a vehicle can visit more than one cluster. A similar problem is the clustered vehicle routing problem (CluVRP), introduced by [8]. In this problem, if a vehicle visits a cluster, it visits all nodes within it. The articles [28,29] present vast reviews about generalized network design problems. In the following, we include the most recent and relevant articles about CluVRP and GVRP.
Ref. [30] provides efficient metaheuristics for solving the CluVRP. Ref. [9] proposes a variant of the CluVRP considering weak cluster constraints. They propose a variable neighborhood search algorithm for solving instances with around 400 nodes. Ref. [31] provides a genetic algorithm for solving the CluVRP. Recently, Ref. [32] implements two heuristics for the CluVRP: particle swarm optimization (PSO) and variable neighborhood search (VNS). They solve instances up to 482 nodes, 97 clusters and 5 vehicles.
Ref. [33] provides a survey of applications for the GVRP. Refs. [28,34] present and compare various mathematical models and exact algorithms for the GVRP. Ref. [35] presents an efficient tabu search heuristic for solving test instances with up to 101 nodes. Ref. [36] addresses the GVRP considering stochastic demands, providing a genetic algorithm for solving instances with 262 nodes. Ref. [37] provides a branch-and-cut-and-price algorithm for the GVRP. Ref. [38] extends the GVRP considering that the clients can belong to more than one cluster. Ref. [39] studies the GVRP with time windows. They propose a column generation-based heuristic for solving instances up to 100 nodes. Recently, Ref. [13] presents the Generalized Median Tour Problem (GMTP). It consists of defining a main route that must visit a set of nodes grouped into disjoint clusters. The route can visit one or more nodes within a cluster. Those nodes not visited by the route are assigned to their closest node in the route (in the same cluster). A hypothetical case study is solved for vaccine distribution in southern Chile.
Unlike the articles mentioned above, we solve a generalized version of the GMTP, considering the use of multiple vehicles and clusters in humanitarian logistics aid. We allow multiple visits from different vehicles to nodes of the clusters. Unlike the CluVRP, in our problem, the vehicles do not necessarily have to visit all stops of a cluster. Like the GMTP, some nodes are assigned to vehicle stops. Moreover, we build clusters using a MILP model and solve a humanitarian distribution aid problem using real data.

Methodology
In the following, we present the methodology used in this study. First, clustering methods based on k-means and p-median are presented. Second, the allocation-routing MILP model is presented. Finally, several solution strategies are presented.

p-Median MILP Model for Cluster Construction
Clustering is a process of grouping data within disjoint areas, usually using a proximity criterion or objects with similar characteristics. A well-known method is k-means, defining clusters using Euclidean distances for partitioning a study area [40].
We propose an ad-hoc capacitated p-median MILP model [41] to build clusters containing ADPs and BPs.This clustering method builds exactly p clusters, minimizing the distances weighted by the demand of BPs plus the distances between ADPs. Figure 2 shows the related scheme to the p-median clustering. Each BP is assigned to a single ADP, and each ADP is assigned to a cluster. Notice that having more than one ADP within a cluster does not imply that all ADPs are selected in the routing stage. In this scheme, ADP 1 is assigned as the cluster's centroid. BPs 4, 5, and 6 are assigned to ADP 1, while ADPs 1 (self-assignment), 2, and 3 are assigned to the same cluster (or ADP 1). Notice that a cluster has limited capacity to accumulate products for serving the potential demand of BPs. Be J the set of ADP, M the set of BP, D i the demand of each BP i ∈ M, Q the maximum demand for each cluster, d ij the shortest distance between nodes i ∈ J ∪ M and j ∈ M ∪ J, and p the number of clusters to construct. The decision variables, objective function, and constraints are defined in the following.

Decision variables
Subject to: Objective (1) minimizes the weighted distance between BPs and ADPs plus the distance among ADPs. Constraints (2) assure that each BP must be allocated with a single ADP. Constraints (3) indicate that each ADP must be allocated to a cluster. Notice that the constraint could include a self-assignment, i.e., a ij = 1 when i = j. Constraint (4) locates exactly p clusters. Constraints (5) indicate that no cluster j is build if there are not allocation variables a. Constraints (6) bound the number of possible allocations to an ADP. Constraints (7) is similar to (6), bounding the maximum number of ADPs on each cluster. Also, those constraints state the relation between allocation and location variables. Constraints (8) state a maximum of allocated demand to a cluster. Finally, constraints (9) indicate the domain of decision variables.

Routing and Allocation MILP Model
The following MILP model allows building aid distribution routes between and within clusters, starting from a central depot (node 0), visiting ADPs, and returning to the depot. Each BP is allocated to an open ADP within a cluster. This model considers that each cluster must be visited by at least one vehicle. A vehicle route within a cluster can visit one or more ADPs. Moreover, given the demand required by the BPs, and the limited capacity of vehicles, an ADP can be visited by more than one vehicle. If a BP is assigned to an ADP within a cluster, a vehicle must serve that ADP.
The notation is the same as the first model plus the following sets and parameters. V the set of vehicles, K the set of clusters, N M k the subset of BPs in the cluster k ∈ K, N J k the subset of ADPs in the cluster k ∈ K, and A is the set of arcs such as A = {(i, j) : i ∈ J ∪ M, j ∈ J ∪ M : i = j} Notice that N J 0 = {0} contains only the depot 0. t ij travel time in the arc (i, j) ∈ A, T is the unloading time on each ADP, D i the demand of each BP i ∈ M, C capacity of each vehicle, and AT is the maximum average travel time allowed for each BP. The decision variables, objective, and constraints are as follows: Decision variables Subject to Objective (10) minimizes the total travel time plus the total unloading time. Constraints (11) and (12) assure that each cluster must be visited and left by at least one vehicle. Constraints (11) force that at least one arc leaves each cluster, while (12) force that at least one arc must enter to each cluster. Constraints (13) assure the connectivity of each vehicle route, stating a flow balance for each ADP and vehicle. In constraints (14) each BP is allocated to a single vehicle and ADP. Constraints (15) and (16)

Vulnerability
We address a vulnerability using the Social Vulnerability Index, SoVi, introduced by [42], to prioritize and focus the humanitarian aid for more vulnerable sectors. It characterizes and measures the economic and social vulnerability conditions before a disaster, interfering with preparation, response, and recovery in a potentially damaged area. In our research, we adopted the SoVi index calculated by [43] for the city of Iquique.
To address the vulnerability in our approach, we use the set of constraints (11)-(24) plus the following parameters and variables: W j is the vulnerability index for the ADP j ∈ J. We assume that all ADPs and BPs within a cluster have the same vulnerability index. The H jv variable indicates the moment where an ADP j ∈ J is served by a vehicle v ∈ V. The µ j variable indicates the maximum service time of ADP j ∈ J. Objective Subject to: (11)-(24) The objective function (25) replaces (10)

Gavish and Graves Constraints
We adapt the well-known Gavish and Graves (GG) flow constraints to our problem [44]. These constraints preclude the subtours apparition, sending one unit of flow to each ADP of the network, and conserving the flow balance for each vehicle and node. Thus, we define the variables f ijv ≥ 0, as the flow from the node i ∈ J to the node j ∈ J of vehicle v ∈ V. The following constraints replace the constraints (21): Constraints (29) state the flow from the central depot. Constraints (30) state indicate the flow balance for all vehicles and nodes. Constraints (31) relate flow and binary variables.
Constraints (32) assure the connectivity for each route, avoiding subtours by tracking the route and time accumulation for each vehicle and node with variables H.
We also try an iterative method (similar to the applied in [13]) to deal with constraints (21). However, the results were not promising, since this method takes a lot of computational time to provide good feasible solutions. So, in the following, we implement and compare both solution strategies (MTZ and GG) in a realistic case study of Iquique, Chile.

Case Study
We address a case study in the city of Iquique, Chile. It has an area of 2242 km 2 , and 191,468 inhabitants [45]. Iquique is a city located in the north of Chile, where there is a seismic gap that could generate an earthquake of M w 8.9 [46] and subsequent tsunami with waves of up to 40 meters [47]. Consequently, intensive preparation in the pre-disaster pre-positioning phase for facility location is required to reduce fatalities and better recovery in case of disasters like earthquakes and tsunamis [48].
The city has 168 educational establishments, which are considered potential ADPs. Notice that 57% of them are under the flooding bound (see the red area of Figure 3, generated in Google Maps [49]). We group the population into 62 groups of 500 meters and preclude all educational establishments in flooding areas. Thus, the instance has 62 BPs and 91 ADPs (see Figure 4). The BPs group 65,264 households demanding around 1135 m 3 of first need kits. Each first need kit has essential items such as food, water, face masks, and blankets.
We chose the International Airport of Iquique Diego Aracena (IQQ), as the central distribution center (depot) where the most humanitarian aid would arrive in a disaster (starting and ending point for vehicles). Note that this airport is located around 40 km from Iquique's downtown.
Twelve vehicles with a capacity of 100 m 3 were considered for distribution. We obtain the travel times for between all pairs of nodes (from Bing Maps) within the city (excluding the airport) using a conservative average speed of 10 km h for travels within the city. Besides, we used an average speed of 50 km h for travels from the airport to the city, because most of these travels are made on the highway connecting the airport with the city. The unloading time for each truck and ADP is estimated in 60 min. Notice that we precluded all floodable segments of routes for each travel time.
For clustering, we used the k-means method and capacitated p-median problem. Notice that a minimum number of 18 clusters is required to assure the feasibility of the clusters methods presented in Section 3.1. Notice the maximum demand accumulated on each cluster is around 100 m 3 . So, we built 18, 20, and 22 clusters to test our methodology under different clustering. We implemented both subtour elimination approaches, GG and MTZ.

Results
All computational experiments were run on a PC with an Intel core i7 processor, and 16GB of RAM. The models were implemented in AMPL and solved with the commercial solver GUROBI 9.1.2., with a time limit of 5 h.
Regarding the clusterization, Figure 5 presents the result of two clusterization (18 clusters). Each cluster is represented with nodes of the same color. The first clusterization (Figure 5a) was built with the capacitated p-median model, while the second was built with the k-means method (Figure 5b). The clustering for both methods is quite similar because both use distance criteria. However, the p-median clustering tends to consider more balanced clusters in terms of accumulated demand of each one. Once determined clusters, the base model is run with GG constraints and MTZ constraints, providing the results of Table 1. #C indicates the number of clusters, Cluster method states the cluster method for each run, STC indicates the constraint set used for subtour elimination, IG is the integrality gap for each run. All the following times are presented in the format hh:mm:ss. ABPT is average travel time for beneficiary points, TVT is the average travel time for the vehicle fleet (discarding the service times), #ADP is the number of visited aid delivery points, AVT is the average travel time per vehicle.
Observe that the travel time is between three to five hours per vehicle in all tests. So, if the 12 vehicles are available and leave the depot simultaneously, the humanitarian aid will be delivered within the critical 72 h after the event. Moreover, it is possible to provide humanitarian aid within 72 h using six vehicles and two trips each. It is important to notice that the most proportion of time is given by unloading times (see the difference between columns TVT and AVT). Regarding the travel times for BPs, they are around 25 min, considering an average speed of 10 km per hour (slow travels in vehicles, motorcycles, bicycles, or walks). The best result for the ABPT is found with 20 clusters, 24 stops, k-means, and GG constraints.
In terms of clustering, the results are similar in magnitude of the average travel time for vehicles and the number of visited ADPs. The better result in 5 h of CPU time is provided using 18 clusters (p-median) and GG constraints.
MTZ and GG constraints show similar behavior concerning the computational performance, presenting an average integrality gap of around 11%. The results indicate that it is less complicated solving the problem with more clusters (observe the IG difference between 18 and 22 clusters). In this case, the clusters are smaller, having fewer sequence combinations within each cluster. However, it does not assure better solutions in terms of the AVT column.
Each vehicle visits at least one cluster in most vehicle routes and three at maximum. Furthermore, we observe that a vehicle visits less than three ADPs within each cluster. Figure 6 shows a vehicle route, visiting the three clusters. It should be noticed that more than one vehicle visits a cluster. It allows to better organize the aid distribution, especially for clusters with many beneficiaries.

Addressing the Vulnerability in Routing
In order to test our methods by considering vulnerability in clusters, we run the same instances of Table 1, using the MTZ constraints. Figure 7a displays the SoVI index in Iquique, where the darker the tonality of the points, the higher the index. Table 2 presents the results of this application, using the same notation of Table 1, and adding the column CPU to present the CPU time required to obtain the optimal solution. This prioritization increases the distribution times because the vehicles need to serve first the more vulnerable clusters (according to the higher SoVi index). Moreover, the number of stops (ADPs) is increased when is compared with the results of Table 1. However, the total distribution time is within the standard time limit of 72 h for first aid kit delivery for all potentially affected populations. In terms of distribution time per vehicle, the best solution is provided by the k-means clusterization with 18 clusters. Figure 7b shows a vehicle route considering vulnerability. Moreover, Table 3 shows the order in which the ADPs are visited on this route.  It can be observed that the model prioritizes these areas with more vulnerability (the bigger the index, the more vulnerable the zone). The ADPs A and B have equal SoVI index and are visited first, while ADPs C and D are visited later.

Scenarios for Road Blocking
Earthquakes can generate damages to the road network. Some roads are blocked, and in others, the speed of traffic is slower [51]. Thus, we generated scenarios where part of the transportation network is damaged. First, we solved the case study considering that the main route from the airport to the city of Iquique is blocked. For this case, we consider a second depot located east of Iquique, corresponding to the Aerodrome of Canchones (see Figure 8). We run an experiment with k-means clustering method, 18 clusters, and the objective (10). The results are presented in Table 4. The results are quite similar (in magnitude) to those obtained in Table 1. The Average time for vehicles (AVT) is around 4 h per vehicle. This result shows that our methodology can be easily adapted to changes in depots.
Moreover, we tested our methodology considering that some road links were damaged. In this case, we consider forbidden links, and some ADPs cannot be used according to the information provided by [48]. Some areas have high potential liquefaction and landslides. Results of this variation are presented in Table 5. Naturally, the times are slightly increased (regarding results of Table 1 because the vehicles have to make detours. Notice that the short distances between ADPs (in the downtown area of Iquique) explain minor variations in distribution times. However, this could be a critical factor in larger and wider areas, or having critical infrastructure like bridges connecting critical areas of a city.

Conclusions
This paper introduced, modeled, and solved a clustered allocation and routing problem for humanitarian aid distribution. We provided a methodology with two stages. First, we clustered the study area in order to organize and focus the humanitarian aid. Second, we modeled the vehicle routing problem using a MILP model. This routing problem considers the design of different routes for vehicles for deliveries within clusters. Each vehicle can visit more than one cluster, entering and leaving a cluster once. Within a cluster, a vehicle can visit more than one delivery point if it is necessary to reduce time. Moreover, two objective functions were tested separately. The first minimized the total distribution time, while the second minimized the maximum distribution time, including a vulnerability index within clusters to prioritize the humanitarian aid.
This methodology was tested for a case study in the city of Iquique, Chile, considering 63 beneficiary points, 92 potential aid delivery points, and 12 vehicles. In this case study, the distribution time is around four hours per vehicle when the objective function is to minimize the distribution times. If the objective function is to minimize the vulnerability to distribute the aid for more vulnerable areas, the distribution time per vehicle is increased to 6-7 h.
We notice that this methodology could be used for aid distribution in different cities where a disaster occurs. Consequently, having a distribution plan before a disaster is crucial for good aid distribution. It allows to estimate the required resources with anticipation and, consequently, opportunely relief the people suffering.
The future work could be extended in several directions. First, heuristics are required to reduce the solution times. Second, more than one depot could be considered. It could support humanitarian aid from the north of Chile or Argentina. Third, in an earthquake, part of the transportation network could be damaged, so a stochastic approach, like the provided by [51], would be helpful to consider this effect. Finally, more than one transportation mode could be used to distribute the aid quickly.