*3.1. Background*

Mozambique is a coastal country in Southern Africa, which borders on the south with South Africa, on the southwest with Eswatini, to the west with Zimbabwe, Zambia, and Malawi, and to the north with Tanzania. On the east side stands the Mozambique Canal, with Madagascar and the Comoros Islands as overseas neighbors. According to the 2019 Human Development Index (HDI) [19], Mozambique is among the 10 least developed countries, ranking 180 out of 189. The population of the country is of 28.9 million inhabitants, approximately (Mozambique National Institute of Statistics) and the most populous province is Nampula, representing almost 20% of the total. Figure 1 shows the location of the Nampula Province withing the Mozambican territory. The Nampula Province is prone to frequent floods due to tropical storms. The frequency and intensity of these phenomena has been growing in the last years. In fact, since 2015, the Nampula Province has been hit by the Tropical Storm Chedza (January 2015), the Tropical Cyclone Dineo (February 2017), a tropical depression (January 2018), the Tropical Cyclone Idai (March 2019), and the Tropical Cyclone Kenneth (April 2019).

This case study focuses on the 2018 low-pressure system that formed in the Mozambique Channel on 13 January 2018 and evolved into the tropical depression stage on 16 January. The tropical depression penetrated the territory of Mozambique from the Nampula Province, more specifically, from the Mossuril district. The storm system consisted of heavy rain, winds of 85 km/h, and affected the provinces of Nampula, Niassa and Cabo Delgado, accumulating 400 mm of rain in less than four days. At the same time, this disaster was exacerbated by the Congo air masses and Tropical Cyclone Berguitta. Figure 2 shows the evolution of the phenomenon through time.

**Figure 1.** Territory of Mozambique and of the Nampula Province, in red. (Source: Wikimedia Commons, license CC BY-SA 3.0)

**Figure 2.** Evolution of the 2018 low-pressure system that formed in the Mozambique Channel on 13 January 2018. (Source: Meteo France, license [20])

In terms of damage, according to the INGC [21], the floods affected more than 80,000 people and killed 34, mainly in the Nampula Province, which accounted for 73,240 of the victims. Additionally, many roads were shut, which hampered the response and rescue operations.

## *3.2. Dataset*

In the dataset, the set *V* is comprised of 38 vertices representing the 18 districts and the five municipalities in the Nampula Province, and 15 intersections. The set *A* consists of 106 arcs, corresponding to the 53 road sections that connect the vertices. The graph is illustrated in Figure 3. At the time of the emergency, the inventories were located in the cities of Nampula and Nacala, corresponding to vertices 1 and 19 in the figure. The traversal time of the arcs, *lengthij*, have been calculated using Google Maps and are expressed in minutes.

 **Figure 3.** Graph representing the Nampula Province.

The report by the INGC [21] provides information relative to the impact of the tropical storm in the north of Mozambique, which includes the Nampula Province. The report illustrates the affected population in each municipality and the roads that have been shut and, therefore, it is possible to devise a single scenario that represents the real impact of the tropical storm in the area. From this point onward, we will refer to this scenario as the *deterministic scenario*.

The *stochastic scenarios* have been generated under the assumption that the tropical storm penetrated the territory of Mozambique from each of the coastal towns in the Nampula Province (i.e., Angoche, Ilha de Mozambique, Larde, Memba, Mogincual, Moma, Mossuril, Nacala, Nacala-a-Velha, and Lunga), including Mossuril that is the original entry point. Thus, 10 scenarios are included in the dataset. Each scenario *ω* ∈ Ω is characterized by a probability, *p<sup>ω</sup>*, the vertices demands, *demand<sup>ω</sup>i* , and the availability of the road sections, *saf eωij* . Due to the lack of information, the scenarios are considered equiprobable, i.e., *pω* = 110 . The demand at each vertex and the availability of the arcs have been estimated by means of regression models that made use of the information presented in the report by the INGC regarding the effects of the tropical storm [21], the 2017 Census [22], and data provided by the National Road Administration of Mozambique. More in detail, the demand of the vertices corresponding to intersections is set to zero. On the other hand, for all the vertices corresponding to population centers and for each scenario, the demand has been obtained from a linear regression model that expresses the logarithm of the ratio of the affected population in each community as a function of the following variables:


All the independent variables are statistically significant and the coefficient of determination is *R*<sup>2</sup> = 0.9221, meaning that the model explains most of the variability in the data.

To determine the availability of the arcs, an exploratory analysis of logistic regression models has been done. Different types of independent variables have been considered: related to the road (e.g., length, materials, type, conditions), related to bridges in the road (e.g., number of bridges, total length of bridges, materials), and geographical (e.g., distance to the entry point, height). Unfortunately, no variable was statistically significant. This is probably due to the large number of missing values in the data, despite having it obtained from official sources. Therefore, the availability of all the pairs of arcs {(*<sup>i</sup>*, *j*),(*j*, *i*)} in the graph is determined using Bernoulli trials with a fixed probability of 0.16981. This probability corresponds to the ratio of road sections interdicted by the storm system over the total number in the Nampula Province, i.e., 953= 0.16981.

 The dataset is publicly available and can be downloaded from [23].

#### **4. Computational Experiments**

In this section, the experiments and their results are presented and discussed. Two groups of experiments are considered. The first group concerns the deterministic scenario. This means that, in this setting, the set of scenarios Ω is comprised of a single scenario that represents the real impact that the tropical storm had on the Nampula Province in terms of affected population and road interdicted. The second group is run on the stochastic scenarios. In this context, the model optimizes the decision based on the average performance over the randomly generated scenarios.

Regarding the objective functions, the model is first solved with respect to the unsatisfied demand (Equation (20)). Then, to understand the relationship between the two time measures considered (Equations (21) and (22)), the payoff matrix is calculated. Calculating the payoff matrix requires solving the model four times. Therefore, a single run involves solving the problem five times. This is explained in Algorithm 1. In the first step, the model is solved while optimizing the unsatisfied demand (Equation (20)). The value obtained is fixed. Therefore, from this point on, the model will produce only solutions with that specific value of unsatisfied demand. Then, objective function (21) is optimized, giving the ideal value of the maximum arrival time (max time). The maximum arrival time is set to max time and objective function (22) is optimized, obtaining the anti-ideal value of the total distribution time (total time). Next, the maximum arrival time is unfixed and the total distribution time is optimized (Equation (22)) to calculate its ideal value, total time. Finally, the total distribution time is fixed to total time and the maximum arrival time is optimized (Equation (21)), obtaining the anti-ideal value max time.



The optimization model has been programmed in GAMS (ver. 29.1.0) and solved using CPLEX (ver. 12.9.0.0) on a Dell Precision 5540 (sourced from Cardiff, UK), equipped with Intel-R CoreTM i9-9880H CPU @ 2.30GHz × 16 and 16 GB RAM. The multithreading option of CPLEX has been used. ACPUtimelimitof1800shasbeensetonalltheoptimizationprocesses.

 The following experiments are carried out:


#### *4.1. Deterministic Scenario: Fixed Inventories*

In this experiment, only the deterministic scenario is considered and the inventories are located in Nampula and Nacala as in the real disaster. The model can determine the road sections to fortify. Table 1 presents the results.


**Table 1.** Results for the experiment on the deterministic scenario and with fixed inventories.

The first column shows the number of fortified road sections. The second column presents the value of the unsatisfied demand (Equation (20)). The third and the fourth columns show the ideal (overlined) and anti-ideal (underlined) value for the maximum arrival time objective function (Equation (21)). The last two columns illustrate the ideal (overlined) and anti-ideal (underlined) value for the total time objective function (Equation (22)).

The first line in Table 1 can be used as a benchmark as it corresponds to the real setting: no road section had been fortified prior to the disaster and there are two emergency inventories located in Nampula and Nacala. The unsatisfied demand is equal to zero, so all the demand center could have been reached from at least one inventory. By increasing the number of fortified road sections (i.e., *Q* ≥ 0), it can be observed that both the maximum time and the total time decrease. However, the maximum improvement is achieved for *Q* = 2, so it would not have been beneficial to fortify more than two road sections. Finally, we can observe that the ideal and the anti-ideal values of both time objectives are the same. Thus, no trade-off between the two time objectives is detected.

#### *4.2. Deterministic Scenario: Model-Defined Inventory Locations*

In this experiment, only the deterministic scenario is considered. However, differently from the previous one, the model can decide both the inventory locations and the road sections to be fortified. Table 2 illustrates the results of the experiment.


**Table 2.** Results for the experiment on the deterministic scenario with model-defined inventory locations.

The first two columns shows the number of inventories and fortified road sections, respectively. The third column presents the value of the unsatisfied demand (Equation (20)). The fourth and the fifth columns show the ideal (overlined) and anti-ideal (underlined) value for the maximum arrival time objective function (Equation (21)). The last two columns illustrate the ideal (overlined) and anti-ideal (underlined) value for the total time objective function (Equation (22)). The first line shows the benchmark from the previous experiment, that is, the solution obtained without fortifying any road segmen<sup>t</sup> and fixing the inventories to the real locations.

As expected, the unsatisfied demand is equal to zero as all the demand centers are connected to at least one inventory. Again it is possible to identify a maximum value for *Q*. According to the results, there is no point in fortifying more than three road sections, as the solution values do not improve. Differently from the previous case, a trade-off between the time measures can be detected. Also, it can be observed that the solution with one inventory and no fortified roads has a lower total distribution time than the benchmark. This result emphasizes the importance of an adequate pre-location of the emergency inventories, as it can lead to faster distribution with fewer resources. Finally, the decrease in the solution times becomes less prominent for the solutions with more than two inventories and one fortified road.

#### *4.3. Stochastic Scenarios*

In the following experiments the set Ω is comprised of the 10 scenarios generated as explained in Section 3.2. Table 3 presents the solution values for different configurations of the parameters *P* and *Q*.


**Table 3.** Solution values for the experiments on the stochastic scenarios.

The first two columns shows the number of inventories and fortified road sections, respectively. The third column presents the expected unsatisfied demand (Equation (20)). The fourth and the fifth columns show the ideal (overlined) and anti-ideal (underlined) value for the maximum arrival time objective function (Equation (21)). The last two columns illustrate the ideal (overlined) and anti-ideal (underlined) value for the total time objective function (Equation (22)). Solution values with an asterisk (\*) are sub-optimal as the execution was halted due to the time limit. The first line presents the benchmark, obtained without fortifying any road segmen<sup>t</sup> and fixing the inventories to the real locations.

Considering the humanitarian context of the problem, it is imperative to satisfy all the demand. Therefore, according to the results, the configurations that should be considered for implementation are:


The choice between them should be driven by the costs of fortifying a road and opening an emergency inventory. In terms of the time objectives, it is possible to detect a trade-off. However, all the instances that present different ideal and anti-ideal values could not identify at least one of the optimal solutions within the time limit. Therefore, these results are not conclusive. Finally, the solution of the instance (*<sup>P</sup>*, *Q*)=(2, 0) can be compared with the benchmark. When using the proposed model to define the location of the inventories, the expected unsatisfied demand improves by 69.94%, corresponding to 6,125.8 people rescued on average. The distribution times are not comparable since the solution of the model reaches more demand centers than the benchmark, which implies a larger distribution operation in terms of demand centers relieved and emergency goods delivered. However, despite that, the expected total delivery time still improves by 43.69%.

#### *4.4. Comparing Deterministic and Stochastic Solutions*

This subsection compares the solutions obtained by the stochastic model with those of a deterministic model that considers one scenario at a time, to understand the usefulness of the presented approach. Table 4 presents the solutions obtained by the two models. The parameters (*<sup>P</sup>*, *Q*)=(4, 2) and only the solutions corresponding to the ideal total time and the anti-ideal max time have been chosen for illustrative purposes. However, similar results have been obtained with all the configurations.


**Table 4.** Solution comparison between the stochastic and the deterministic model, for (*<sup>P</sup>*, *Q*)=(4, <sup>2</sup>).

The table illustrates in the first row the solution of the stochastic model and in the following 10 rows the individual solutions of the deterministic model, one for each scenario. The first column identifies the scenario considered. The second and the third column present the set of inventory locations and the set of fortified edges, respectively.

From the table, it is possible to observe that addressing each scenario separately would not allow the identification of the best global solution, as each deterministic solution is *ad hoc*. In fact, there is little overlap between the deterministic inventory locations and the stochastic one. This is clearly illustrated in Figure 4, which shows the bar plot of the frequencies of the inventory locations in the deterministic solution.

**Figure 4.** Bar plot of the frequencies of the inventory locations in the deterministic solutions. The x-axis represent the locations and the y-axis their frequencies in the deterministic solutions.

According to the plot, if a decision-maker were to use the most frequent locations to select the best configuration of inventory locations, the set {3, 9, 15, 20} would be chosen. However, the performance of this solution is strongly sub-optimal, as shown in Table 5.


**Table 5.** Objective function values comparison between the stochastic and the deterministic solutions.

The table presents the objective function values of the stochastic solution, {3, 14, 16, 19}, and of the deterministic solution, {3, 9, 15, 20}. The first column shows the model considered. The second column illustrates the inventory locations identified by the model. The remaining columns present the corresponding objective function values. The solutions are compared on the stochastic model, that is, considering all the scenarios. Generally speaking, the deterministic solution is expected to perform worse than the stochastic, as the former is a heuristic solution obtained by solving each scenario independently and then choosing the most frequent locations across all the solutions (Figure 4), while the latter is the optimal solution to the stochastic model. From the table it can be seen that the deterministic solution is strongly sub-optimal. In fact, in terms of unsatisfied demand, the stochastic solution allows to relieve 1,175.3 more affected persons, on average, than the deterministic, corresponding to an improvement of 93.99%. Despite attending more demand, the stochastic solution is also more efficient in terms of distribution time, improving by 2.26% the max time and by 52.90% the total time of the deterministic solution.

Although this analysis only considered the inventory locations, similar conclusions can be drawn on the road fortifications.

## **5. Conclusions**

This paper presents the first pre-disaster stochastic model that jointly optimizes the location of emergency inventories and the fortification of transportation network elements in Humanitarian Logistics. Another important feature of the model is that it is parsimonious in terms of data requirements. This is fundamental to be able to realistically apply the model to underdeveloped and developing countries, that usually have limited information on resources available and infrastructure status. The problem has been modeled as a two-level lexicographic problem, where the first level minimizes the expected unsatisfied demand, while the second level considers two time measures: expected maximum arrival time at the demand vertices, and expected total distribution time. As a weak trade-off is detected between the time measures, their relationship is assessed using the pay-off matrix, rather than relying on more complex multicriteria methods. The model is tested on a novel dataset based on a case study on the storm system that hit the Nampula Province (Mozambique) in January 2018. The model is used to evaluate the real response to the emergency and identify what the best course of action should have been. The experiments show that the solution obtained by the model is better than the current policy and improves on the expected unsatisfied demand by 69.94% (corresponding to 6,125.8 more people relieved, on average) and on the expected total delivery time by 43.69%. Also, the model provides two solutions that allow to service all the demand under all the scenarios. The choice between these two solutions depends on the costs of opening an emergency inventory and fortifying a road section, and it is left to the decision-maker. Finally, we illustrated empirically how useful it is to integrate uncertainty in the optimization.

Differently from other investigators, our future research plans do not involve making the model more complicated by factoring in additional operations that, in turn, require more information, as this would be against the rationale behind this contribution. Our objective is to devise models that help decision-makers in emergency managemen<sup>t</sup> as much as possible while keeping the data required to a realistic level. For this same reason, we are currently considering to represent the uncertainty in the disaster outcome by applying robust optimization rather than stochastic optimization. Moving from an expected cost objective to a minimax formulation has two major advantages. First, it does not need the probability distributions of the uncertain parameters. Second, it results in more conservative

solutions that cater for worst-case outcomes [24], which is a desirable feature in the context of a disaster. A different approach which would be worth exploring is formulating the problem to consider Recoverable Robustness [25]. A solution is recovery robust if it can be recovered by limited means in all likely scenarios. This provides a middle ground between classical Robust Optimization and Stochastic Programming. Another interesting extension to the model would be introducing a temporal dimension to consider multiple emergencies over time. This is motivated by the application context, as many disasters (e.g., floods and hurricanes) are seasonal. Therefore, the decision-maker could plan the fortification of the road network over multiple periods, while prepositioning the inventories before every emergency. Finally, an alternative line of research concerns the generation of stochastic scenarios. Namely, we are interested in looking for data-efficient ways of building realistic scenarios which take into account correlations between vertices and arcs disruptions.

**Author Contributions:** Conceptualization, F.L., J.M. and B.V.; methodology, F.L. and J.M.; software, J.M., F.L. and B.V.; validation, B.V., J.M. and F.L.; formal analysis, F.L., B.V. and J.M.; investigation, B.V. and F.L.; resources, B.V.; data curation, J.M.; writing—original draft preparation, F.L.; writing—review and editing, F.L. and B.V.; visualization, J.M. and B.V.; supervision, B.V.; project administration, F.L. and B.V.; funding acquisition, B.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was partially funded by the European Commission's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie, gran<sup>t</sup> number MSCA-RISE 691161 (GEO-SAFE); the Government of Spain, gran<sup>t</sup> MTM2015-65803-R; the Complutense University of Madrid, scholarship *Ayudas para la realización de estancias en proyectos de cooperación para el desarrollo sostenible*. All the financial support is gratefully acknowledged.

**Acknowledgments:** The authors would like to thank the following entities, agencies and organizations for all the help and support provided: National Institute of Statistics (Mozambique), National Road Administration (Mozambique), National Institute of Disaster Management (Mozambique), Red Cross (Spain), Maputo International Airport, and the NGO ONGAWA.

**Conflicts of Interest:** The authors declare no conflict of interest.
