**Decision variables**

	- *Vkit* : Number of vehicles of type *k* ∈ K located at node *i* ∈ N at the beginning of period *t* ≤ T .

#### **Attribute and deviation variables**


The four attributes of the model (number of people with high priority evacuated, number of people with normal priority evacuated, total evacuation time and operation cost) are considered jointly within a lexicographical goal programming approach. Each attribute *m* ∈ M is associated with an aspiration level *tgm* that represents a certain value that would be satisfying for the decision maker and we seek to achieve. The unwanted deviations from the aspiration levels *tgm* are measured by the deviation variables *DVm*, in such a way that a zero value in one of them means the associated aspiration level has been achieved. We use a lexicographical approach because some of the considered attributes are incomparably more important than others and thus no trade offs are possible among them. This way, the objective function consists in minimizing lexicographically the four deviation variables following the order of priority of the considered attributes: *DVE* (first level of priority, critical population evacuated), *DVE*' (second level of priority, non critical population evacuated), *DVTM* (third level of priority, total evacuation time), *DVC* (fourth level of priority, operation cost).

#### *3.2. Mathematical Model*

The proposed model is given by Equations (1)–(25): the deviation variables of the objective function (1) are to be minimized lexicographically, subject to constraints (2)–(25), that define feasible evacuation plans.

$$\text{LexMin}\left\{ \left[ DV\_{\mathcal{E}} \right]\_{\prime} \left[ DV\_{\mathcal{E}'} \right]\_{\prime} \left[ DV\_{TM} \right]\_{\prime} \left[ DV\_{\mathcal{C}} \right] \right\} \tag{1}$$

*Mathematics* **2020**, *8*, 648

$$\sum\_{\substack{j|(\vec{\mu})\in\mathcal{E}}} \sum\_{\substack{k}} \sum\_{\substack{t'\le t-\tau^{k}\_{\vec{\mu}}\\t'\le t}} P L^{hk}\_{j\text{it}'} + \sum\_{t'\le t} q p^{h}\_{\vec{\mu}t'} = \sum\_{j|(\vec{\mu})\in\mathcal{E}} \sum\_{k} \sum\_{t'\le t} P L^{hk}\_{\vec{\mu}jt'} + P^{h}\_{\vec{\text{tt}}} \,\forall h, i, t \tag{2}$$

$$wp^hP^h\_{\,it} \subseteq cp^h\_{\,i} \,\forall t, h, i \in \mathcal{N}\mathcal{H} \cup \mathcal{N}\mathcal{S} \tag{3}$$

$$\sum\_{h} wp^{h} PL\_{ijt}^{hk} \le vpc^{k} VL\_{ijt}^{k} \; \forall i, j \in \mathcal{N} \; | \; (ij) \in \mathcal{E} \; \forall k, t \tag{4}$$

$$E = \sum\_{h \mid \forall y p^h = \max\_h \{type^h\}, i \in \mathcal{N}\mathcal{H} \cup \mathcal{N}\mathcal{S}} P\_{i\mathcal{T}}^h \tag{5}$$

$$=\sum\_{h\mid typ^h=\min\_k\{typ^h\}, i\in\mathcal{N}'\mathcal{H}\cup\mathcal{N}'S} P\_{i\mathcal{T}}^h\tag{6}$$

$$\sum\_{\substack{j(\{i\}) \in \mathcal{E} \ t' \le t - \tau\_{ji}^k}} \sum\_{l' \le t - \tau\_{ji}^k} V L\_{j \text{it}'}^{k'} + v a\_i^k = \sum\_{j(\{i\}) \in \mathcal{E}} \sum\_{t' \le t} V L\_{i \text{j} t'}^{k} + V\_{\text{it}}^k \,\forall i, k, t \tag{7}$$

$$\sum\_{i} v a\_{i}^{k} = \sum\_{i} V\_{\text{ir}\mathcal{T}}^{k} \,\forall k \tag{8}$$

$$T = \sum\_{t} t \left( BT\_t - BT\_{t-1} \right) \tag{9}$$

$$BT\_t \ge BT\_{t-1} \,\forall t \tag{10}$$

$$\sum\_{h|\text{type}^h=\max\_h\{type^h\}} \sum\_{i,j|(ii)\in\mathcal{E}} \sum\_k \sum\_{\substack{t'=\lceil t-\tau\_{ij}^k\rceil}} PL\_{ijt'}^{hk} \le (1 - BT\_t) \sum\_{h|\text{type}^h=\max\_h\{type^h\}} \sum\_{i,t'} qp\_{it'}^{h\_i} \forall t \tag{11}$$

*E*

$$T' = \sum\_{t} t (BT\_t' - BT\_{t-1}') \tag{12}$$

$$BT\_t' \ge BT\_{t-1}' \,\forall t \tag{13}$$

$$\sum\_{\substack{h|tpp^h=\min\_h\{typ^h\}\ i,j|(ij)\in\mathcal{E}}} \sum\_{\substack{k}} \sum\_{\substack{t'|t'=\lceil t-\tau\_{ij}^k\rceil}} PL\_{ijt'}^{hk} \le (1-BT\_t') \sum\_{\substack{h|tpp^h=\min\_h\{typ^h\}\ i,t'}} \sum\_{i,t'} qp\_{it'}^{h} \forall t \tag{14}$$

$$TM \ge T \tag{15}$$

$$TM \ge T'\tag{16}$$

$$\text{Cost} = \sum\_{h,i,j,k,t|(ij)\in\mathcal{E}} d\_{ij}^k (f c\_{ij}^k V L\_{ijt}^k + \nu c\_{ij}^k P L\_{ijt}^{hk}) \tag{17}$$

$$Cost \le b$$

$$Cost \ge 0$$

$$E + DV\_E - P\_E = t\mathbb{g}\_E \tag{20}$$

$$E' + DV\_{E'} - P\_{E'} = \text{tg}\_{E'} \tag{21}$$

$$TM - DV\_{TM} + P\_{TM} = t\mathbb{g}\_{TM} \tag{22}$$

$$\text{Cost} - DV\_{\mathbb{C}} + P\_{\mathbb{C}} = \text{tg}\_{\mathbb{C}} \tag{23}$$

$$\mathbb{P}^{h}\_{\text{it}"} \mathbb{P}L^{h}\_{\text{ijt}"} V^{k}\_{\text{it}"} V L^{k}\_{\text{ijt}"} \mathbb{E}, \mathbb{E}', \mathbb{T}, \mathbb{T}', \mathbb{T}' \mathbb{M} \in \mathbb{Z}^{+} \text{ \textquotedbl{}t'}, h, k, t \tag{24}$$

$$BT\_{t\prime}BT\_{t}^{\prime} \in \{0, 1\} \; \forall t \tag{25}$$

Equations (2) take care of the conservation of the flow of people along the network: for each *i*, *h* and *t*, the amount of people of type *h* that have been transported to *i* up to time period *t* plus the amount of people of type *h* that arrived on their own to node *i* up to time period *t* must equal the amount of people of type *h* that leave node *i* up to time *t* plus the ones that stay there. Node and vehicle capacities are imposed by (3) and (4), respectively. As a result, the number of people of each type that can be evacuated is calculated by (5) and (6). Analogously, Equations (7) and (8) are the vehicle flow conservation constraints.

The time period at which every high (respectively normal) priority evacuee has been moved to a safe area is given by Equation (9) (respectively (12)), while (10) and (11) (respectively (13) and (14)) ensure binary variables *BT* (respectively *BT*) are defined properly. The maximum evacuation time is

represented by *TM*, following (15) and (16), and constraints related to the available budget are (17)–(19), defining the total cost as the sum of the fixed cost of moving an empty vehicle and the variable cost for transporting people.

Since we are considering a multi-objective lexicographical goal programming model, the optimal value obtained at each level leads to new constraints that are added to the following levels to ensure that the goals from previous levels are achieved. Objectives with highest priority are those related to the amount of people that is possible to evacuate. Hence, the first and second levels consist in minimizing the deviations from the stated aspiration levels of the total number of people evacuated with high or normal priority, as given in (20) and (21). The third level considers the prompt evacuation of the population, minimizing the deviation with respect to the given aspiration level, as stated in Equation (22). Finally, the minimization of the deviation of the total cost with respect to its aspiration level (see (23)) is taken into account.

#### **4. Case Study: Palu, Indonesia**

This section describes the case study used to evaluate the performance of the proposed model. It is based on the earthquake and tsunami that hit the island of Sulawesi, Indonesia, in September 2018. The Indonesian archipelago is known for its catastrophic earthquakes as well as the occurrence of landslide induced tsunamis. The earthquake and tsunami we are considering were caused by a strike-slip faulting event at shallow depth that occurred within the interior of the Molucca Sea micro-plate, which is part of the Sunda tectonic plate, see [41].

On 28 September 2018, 10:02:44 GMT/UTC hour (local hour, WIB, 18:02:44) a series of strong earthquakes struck mainly Palu, the capital of the Indonesian province of Central Sulawesi, and Donggala. The strongest was a 7.4 M earthquake only 10 km deep with its epicenter close to Palu. The earthquake, according to the Humanitarian Country Team [42], triggered a tsunami, whose waves reached up to three metres, that hit the northern parts of Sulawesi island and resulted in liquefaction and landslides. Furthermore, a total of 170 aftershocks occurred after the main one. As a consequence of the liquefaction and aftershocks, some roads from/to Palu were not accessible for several hours and both the port and airport were servery damaged, highly complicating the Search & Rescue and evacuation management. In particular, the BNPB (National Disaster Management Agency) declared the urgen<sup>t</sup> need for planes with the ability to take off and land in short runway (specifically, Hercules C130) due to the cracks at the airport runway.

Based on the updated information from WHO, see [43], 832 persons died, 580 were injured, more than 16,000 were displaced to temporary shelters and dozens of houses were damaged in the first 24 h. The estimated exposed population was more than 310,00 in Donggala regency and more than 350,000, in Palu.

## *4.1. Data*

The case study has been created by using information obtained from AHA Reports [44–51]. The satellite images and geospatial data of the affected area were obtained from the Copernicus Emergency Management Service website, see [52]. The complete affected area, illustrated in Figure 1, is located in the Sulawesi island and divided into 18 grading maps (named Areas of Interest), covering the cities and villages that suffered the shake of the earthquake and/or the impact along the coast. In particular, we have focused on the two most affected Areas of Interest (AoI), that host the most exposed population (the ones corresponding to the two squares inside the red rectangle of Figure 1), namely AoI 7 (Palu) and AoI 11 (Palu East). Their grading maps are shown in Figure 2, where AoI 7 corresponds to the left-hand-side figure and AoI 11 to the right-hand-side one.

The network of the case study contains 44 nodes representing locations of the city of Palu. According to their characteristics, they are classified as:


**Figure 1.** Areas of Interest of the earthquake and tsunami, Indonesia.

**Figure 2.** Considered areas (Source: [52]). (**a**) AoI 7, Palu and (**b**) AoI 11, East Palu.

The nodes are connected by 76 edges, whose lengths were determined by using GoogleMyMaps technology. According to the information obtained from the maps of Copernicus [52], they represent bridges, level crossings, tunnels or roads which suffered different levels of damage, as blocked or shattered, seriously damaged, partially damaged or passable. Furthermore, they may be highways, single carriageways, local, within towns or unpaved roads. The complete graph is shown in Figure 3.

**Figure 3.** Complete Graph.

The condition of streets and roads, due to the tsunami and the liquefaction of the ground at some areas, made necessary the participation of a diverse number of vehicles for the tasks of people search and evacuation, including helicopters, regular and inflatable boats, army trucks and ambulances. Terrestrial vehicles cannot travel by blocked or shattered arcs and it has been supposed that seriously and partially damaged connections suffer a penalization of 50% and 25% of their maximum velocity, respectively. At the beginning of the operation, ambulances are available at hospitals, while army trucks can be found at several nodes. Moreover, only blocked or shattered roads can be travelled by inflatable boat, because of the floods, water movement and the liquefaction of the terrain. At the beginning of the operation, all inflatable boats are at affected nodes, because they either belong to the neighborhood or the Search & Rescue teams leave them there. Finally, helicopters can travel following

a straight line between any pair of points with heliport available and are located at hospitals and at the airport at the beginning of the operation.

Some characteristics of the vehicles such as capacity (number of not badly injured people that can be transported by each type of vehicle), mean velocity and costs are presented in Table 1, together with the maximum amount of vehicles that can be used in the operation. The fixed cost is independent of the load of the vehicle and is paid only per unit of distance; however, the variable cost depends both on the traversed distance and on the cargo transported, being paid per unit of distance and unit of load (u.l.).


**Table 1.** Characteristics of available vehicles.

The data regarding affected people and their classification into the different categories considered before have been estimated from several public reports. For example, the data about the number of inhabitants and affected in each AoI have been obtained from the maps of Copernicus [52], as displayed in Table 2. The report of BNPB [55], already mentioned, provides the number of evacuees by 1 October 2018: a total of 48,025 persons at hospitals and temporary shelters out of 60,424, the total affected population. Therefore, only around 20% of the total population evacuated to relative homes or hotels.

**Table 2.** Population data obtained from the information from Copernicus.


According to the Indonesia Humanitarian Country Team [56], women, unaccompanied minors, adolescent mothers and other vulnerable groups such as people with disabilities, medically homebound, poor, people living with HIV or AIDS and sexual minorities, are at increased risk after the occurrence of a disaster. In fact, it has been estimated that there were around 350,000 cases of abuse, exploitation, forced or early marriage as a consequence of the particular disaster of this case study. In addition, there were cases of food shortages due to the closure of markets, roads and transports. For this reason, besides the classification according to their capacity, shelters have been classified into secure and non secure for vulnerable groups. These groups, together with severely injured people and pregnan<sup>t</sup> women, have been assigned a high priority. In particular, we have considered the following population groups: severely injured adults (SIA), pregnan<sup>t</sup> women (PW), unaccompanied minors (UM), women susceptible to gender violence (GBV) and the rest of the population (RP). According to Stepanov and Smith [57], most of the remaining population who should be able to evacuate on their own usually move to homes of relatives, hotels or other shelters. The Indonesia Humanitarian Country Team [56] established on 10.05.2018 that during the first periods of the response to this disaster, the risk of violence, exploitation or abuse is heightened and pre-existing gender inequalities may be exacerbated with an environment of impunity, where near 100,000 people were stated as under high risk. It has been assumed that at every pick-up node at the affected area, at the beginning of the operation, there was the same number of persons to be evacuated, following an uniform distribution. There are 14 nodes in AoI 7 and 2 in AoI 11. The details and number of people to be evacuated at the beginning of the operation at each AoI are given in Table 3.


**Table 3.** Population characteristics.

Additionally, people are assumed to be arriving to the pick up nodes dynamically, according to their feeling about the disaster and the Search & Rescue operations. The number of people arriving during each time period is estimated through Equation (26). It is supposed that the feeling regarding the potential danger to the population decreases along time. For this reason, a decreasing exponential function is used. For the first time period, we consider the data of people to be evacuated at the beginning of the operation, displayed in Table 3; in addition, after one sixth of the time span has gone by, we assume that no more people arrive at the pick-up points.

$$q p\_{it}^h = \lfloor q p\_{i0}^h e^{\ln(q p\_{i0}^h \* \frac{1-t}{1/6 \times 7-1})} \rfloor \quad \forall t \le \mathcal{T} / 6, \forall i \in \mathcal{N} \mathcal{A} \tag{26}$$

The capacities of the shelters for different types of people are shown in Table 4. The data was obtained from several reports, such as [53–55]. The authorities prioritized the evacuation of those who were injured or ill. There were flights between Mutiara Sis Al Jufri Airport, in Palu, and nearby airports, mainly at cities like Balikpapan (Sulaiman Seppinggan International Airport), Makassar (Sultan Hasanuddin International Airport), Manado (Sam Ratulangi International Airport) and Jakarta (Halim Perdanakusuma International Airport).


**Table 4.** Shelter characteristics and capacity distributions. PW: pregnan<sup>t</sup> women. UM: unaccompanied minors. GBV: women susceptible to gender violence. RP: rest of the population.

The time horizon considered corresponds to 72 h, discretized into 24 three-hour periods, from 28 September 2018 10 a.m. to 1 October 2018, 10 a.m. According to the Humanitarian Country Team [42], on 9 October 2018, the HCT's Central Sulawesi Earthquake Response Plan requested US\$ 50.5 million for the relief activities after the disaster for the assistance of 191,000 people over three months. As an estimation only for the operations and scope included in our case study, we have considered an available budget of 200,000\$.
