A robust bi-objective location-routing model for providing emergency medical services

Disastrous situations are always full of ambiguities in terms of their frequency and severity. Therefore, relieving the affected people needs effective planning in terms of the service resources to decrease injuries and losses. In a crisis, wasted time may change life to death (Altay, 2012).

The process of planning, management and the control of the flow of resources to provide services for people after the occurrence of disasters is known as emergency or relief logistics (Sheu, 2007; Kim and Jayakrishnan, 2010; Santos et al., 2010). As the severity of disasters has been increased exponentially around the world since 1950 (Özdamar and Ertem, 2015), the importance of accurate planning with the consideration of its complexity and uncertain nature has become more than ever. In such a situation, there are challenges which make the problem highly complex. Challenges can relate to the planning to provide services for injured people, limited time and resources and the inability to estimate an accurate number of injured people.

In general, the humanitarian relief chains (HRCs) seek high speed when providing emergency requirements for injured people in order to minimize the number of injuries and mortality. This can be done through the efficient and effective allocation of limited resources (Tofighi et al., 2016). Relief logistics include all of the processes of estimating, supplying, transporting, maintaining and distributing the commodities, equipment, services and all of the necessities to both injured people and relief groups. All of these necessities must be provided in the minimum possible time (proper time) to the designated locations (proper location) and in the required quantities for the intended people and groups (certain individuals) via the scientific and accurate method (proper method) (Özdamar et al., 2004).

Three steps are considered for the response process, including (1) demand point collection, (2) demand classification based on severity or urgency and (3) solving the routing problem and sending the vehicles (Luca et al., 2015). In this regard, the location and routing decisions related to relief distribution planning in a post-disaster situation play a key role in providing timely emergency services.

Barbarosoğlu et al. (2002) focused on modeling the helicopter relief operation and flight planning in the Turkish Army during a disaster and they indicated its role in planning within the operational areas and in limiting the service-providing time. Yi and Kumar (2007) solved the logistics problem in disaster relief operations. They proposed ant colony optimization (ACO) to provide an optimal dispatching planning of commodities to the distribution centers in the affected areas. Two main research studies on disaster relief and evacuation logistics were suggested by Yi and Ozdamar (2007) and Özdamar and Yi (2008). Their studies aimed to minimize the total service delay. An integer programming model was then proposed by Altay (2012) for the capability-based location-allocation problem considering multiple resources. Salman and Gul (2014) presented a novel multi-period mixed-integer linear programming (MILP) model for locating additional emergency units and transporting the casualties to hospitals by ambulance. They investigated a case study of the Istanbul earthquake.

Leiras et al. (2014) provided a comprehensive literature review on humanitarian logistics research including trends and challenges. Rath and Gutjahr (2014) proposed a routing-location model for relief in disastrous conditions. They expressed that intermediate warehouses should be constructed as temporary and outpatient hospitals in order to perform the relief steps and to meet the injured people's needs at a higher speed.

Halskau (2014) investigated marine helicopter routing through a hub location problem. The objective was to minimize the average number of casualties in emergency conditions. Ahmadi et al. (2015) presented a new structured method for solving the routing problem in disastrous conditions. According to their proposed approach, information was imported gradually over the time and it attempted to design the routes in such a way that minimal changes occurred in routes at the planning horizon time. This method was investigated in several numerical examples. Knyazkov et al. (2015) evaluated dynamic ambulance routing for transferring patients in St. Petersburg. In fact, the objective of the problem was to reduce the transferring times of injured parties to hospitals to save the patients' lives. This was accompanied by selecting an optimal route and a target hospital for an ambulance in a large city, resulting in many complications. They considered dynamic features for an urban environment, including traffic flow, population density and hospital capabilities. This led to the determination of emergency medical services (EMS) quality. They investigated the problem through a case study in St. Petersburg and concluded that the traffic conditions had a significant effect on decision-making and that it might change the entire route as well as the target hospital. Chen and Yu (2016) investigated the problem of a temporary facility location for EMS during a disaster concerning the transportation and demand rate in this situation. Maghfiroh et al. (2018) investigated a dynamic vehicle routing problem (VRP) for the last mile distribution during a disaster. They developed a modified simulated annealing (SA) and variable neighborhood search (VNS) algorithm to solve the problem with the aim of minimizing the total traveling time. An improved shuffled frog leaping algorithm (SFLA) was proposed by Duan et al. (2018) to study the dynamic emergency vehicle dispatching problem. They considered response time, accident severity and accident time windows as the main factors to develop an emergency vehicle dispatching model.

Zhang et al. (2018) presented a sustainable multi-depot location-routing problem (LRP) considering uncertain information. They concurrently took into account the travel time, emergency relief costs and CO₂ emissions in their model via uncertainty theory. They sought to solve the problem using a hybrid genetic algorithm (GA) and uncertain simulations. An integrated location-inventory-routing problem (LIRP) was proposed for pre- and post-disaster management by Tavana et al. (2018). They considered a multi-echelon humanitarian logistics network to provide an appropriate flow of relief products. They solved the problem using an improved non-dominated sorting genetic algorithm II (NSGA-II) algorithm. According to the discussed literature up to this stage, our model considers refueling the vehicles and the use of helicopters together with ambulances.

In the field of special studies on robust optimization (RO) in relief logistics, Bozorgi-Amiri et al. (2011), for the first time, introduced a multi-objective robust stochastic programming approach for relief logistics under uncertain conditions. A dynamic multi-objective robust LRP for relief distribution was suggested by Bozorgi-Amiri and Khorsi (2016) considering uncertain demands, travel time and the cost parameters. They implemented the ε-constraint method to cope with the multi-objectiveness of the problem. They tested the applicability of their proposed approach on a real case study. A robust stochastic vehicle routing and scheduling method for bushfire emergency evacuation were investigated by Shahparvari and Abbasi (2017) related to an Australian case study. They took into account the road availability and disruptions in their proposed model and they developed a greedy solution method to solve the problem. Veysmoradi et al. (2018) introduced a multi-objective open LRP for relief distribution networks considering split delivery and multi-mode transportation under uncertain conditions. They formulated the problem as a mixed-integer non-linear programming (MINLP) model, and they provided a hybrid solution using RO and fuzzy multi-objective programming. Liu et al. (2018) developed a RO model for relief logistics planning considering uncertain demands and transportation time. They implemented their proposed model in a case study problem in a city of China that had recently suffered from an earthquake. They proposed optimal management policies using sensitivity analyses.

After reviewing the literature, it can be found that considering the concurrent transportation planning of ambulances and helicopters has not been investigated adequately. On the other hand, there has been a lack of attention to response time as the main objective against the total cost minimization. Time plays a critical role in the relief logistics system. More importantly, during a disaster, vehicles should not be confronted with fuel supply problems since this disrupts the relief operations and it increases the response time. This subject was not considered in most studies as well. Moreover, uncertainty is regarded as one of the most basic features of the relief logistics system which is studied in our problem using the RO technique.

Since the investigated problem is an extension of LRP and it is classified as an NP-hard problem (Ghaffari-Nasab et al., 2013), it cannot be solved using exact methods within a reasonable computational time. For this purpose, many researchers have worked on different multi-objective algorithms such as greedy-search-based multi-objective GA (Chang et al., 2014), multiobjective particle swarm optimization (Ghasemi et al., 2019) and an NSGA-II (Vahdani et al., 2018) in order to solve the multi-objective problems that are particular to relief logistics. In this research, a multi-objective algorithm for the SFLA is developed to solve the problem approximately. SFLA is known as a high-quality algorithm that is used to solve optimization problems (Rahimi-Vahed and Mirzaei, 2007; Niknam et al., 2011; Arshi et al., 2014). Furthermore, the efficiency of the algorithm is evaluated in comparison with the ε-constraint method and the NSGA-II for small- and large-sized problems respectively.

In other words, a robust multi-objective optimization model is essential to provide effective, efficient and equitable relief logistics. Efficient algorithms are also needed to find good quality solutions in the aftermath of the disaster. In this study, we focus on these requirements as the main contributions.

The remaining sections of this paper are organized as follows. Section 2 describes the problem and it presents the proposed robust mathematical model. The suggested solution approaches include the ε-constraint method, SFLA and NSGA-II as introduced in Section 3. Section 4 represents the computational results including the algorithm comparison, robust and deterministic problem comparison and sensitivity analysis. Finally, the concluding remarks and outlook of research are explained in Section 5.

In this paper, a relief network is considered which includes fixed hospitals, potential temporary hospitals, transfer points, ambulance stations, helicopter stations, demand points and refueling points, and the number and location of each one are predetermined. This network is defined within the operational area where many injured people are waiting for relief. To simplify the problem, only two types of injured people are defined, namely, those who are treated at the site of the incident (green-type) and those who should be transferred to the hospital, due to the severity of their injury (red-type) (Tofighi et al., 2016). Moreover, the potential points for facility establishment are predetermined in the relief network.

The overall structure of the proposed supply chain is as shown in Figure 1. In this figure, the red star points represent the red-type patients (acute care needed), green star points represent the green-type patients, blue circle points represent ambulance stations, blue triangle points represent transfer stations, blue hexagonal points represent potential temporary hospitals, blue square points represent fixed hospitals, blue rectangular points represent helicopter stations and blue dodecagon points represent gas stations. Since the landing and taking-off of the helicopters occur at pre-defined locations, we need to consider transfer stations in the network (Bozorgi-Amiri et al., 2017). In fact, transferring injured people to the helicopter depot is done by an ambulance and then, transferring them to the hospital is performed by a helicopter.

In conclusion, our proposed model attempts to concurrently improve the efficiency, effectiveness and equity of the allocation of critical resources during the response phase. The measure of efficiency is denoted by the maximum time required by ambulances to leave disaster-affected points (Huang et al., 2012). The measure of effectiveness can be represented by the speed and sufficiency of the deliveries (Huang et al., 2012), which is studied by considering the assignment of helicopters to provide instantaneous services. Finally, the equity measure is addressed by the time window of each point to provide fair distribution among injured people (Gutjahr and Nolz, 2016). The arrival time of the ambulances at each point should occur within a given interval.

The main assumptions of the problem are as follows:

1. The distance between different points of the network is deterministic.

2. Ambulances and helicopters have limited capacity for providing services for patients.

3. Ambulances have a limitation of fuel consumption.

4. There are several potential points for the construction of temporary hospitals and transfer points.

5. Patients are divided into two groups of green and red.

6. Each ambulance can carry one or more red patients up to its capacity to the hospitals or transfer points.

7. Each ambulance can initially handle a green-type patient and then carry the red-type patient to the hospital or transfer point.

8. Hospitals have enough capacity to service all the red-type patients.

9. The number of patients after the incident is uncertain.

10. A hard time window is defined for each affected area (demand point).

11. Each helicopter should land at a transfer station and then transfer the patients to the hospitals.

12. The time spent by helicopters is instantaneous.

The following symbols are used in this paper. Note that min, lit and man are the abbreviations of minute, liter and the number of people, respectively.

$I$⁠: Set of all points (⁠$i,i′,j,j′∈I$⁠)

$H$⁠: Set of fixed hospitals (⁠$H⊂I$⁠)

$TH$⁠: Set of potential temporary hospitals (⁠$TH⊂I$⁠)

$TP$⁠: Set of transfer stations (⁠$TP⊂I$⁠)

$WS$⁠: Set of ambulance stations (⁠$WS⊂I$⁠)

$HS$⁠: Set of helicopter stations (⁠$HS⊂I$⁠)

$AA$⁠: Set of demand point (⁠$AA⊂I$⁠)

$FP$⁠: Set of gas stations (⁠$FP⊂I$⁠)

$K$⁠: Set of ambulances (⁠$k∈K$⁠)

$HC$⁠: Set of helicopters $h∈HC$

$tij$ (min): Travel time between the two points i and j by each of the ambulances

$QAk$ (man): Capacity of the kth ambulance

$QHh$ (man): Capacity of the hth helicopter

$DRi$ (man): Number of red-type patients at the ith disaster-affected point, which is defined in an uncertainty range

DRi (man): Fluctuation rate of the number of red-type patients at the ith point

DRi (man): Estimated number of red-type patients at the ith point

$DGi$ (man): Number of green-type patients at the ith point

$STi$ (min): Duration of providing services for green-type patients at the ith point

$SGi$ (min): Duration of preparing the red-type patients for being transferred from the ith point

$FAk$ (lit): Primary fuel of the kth ambulance

$FBk$ (lit): Fuel of each ambulance after refueling at the gas station

$HPhi$⁠: Binary parameter will be equal to 1 if the helicopter h is at the ith helicopter station at the starting time

$α$ (lit/min): Conversion factor of the elapsed time to the consumed fuel in ambulances

$cxij$ ($): Cost of covering the i-to-j route by each ambulance

$cyij$ ($): Cost of covering the i-to-j route by each helicopter

$fci$ ($): Cost of constructing the transfer station i

$fdi$ ($): Cost of constructing the temporary hospital i

$Ei$ (min): Lower bound defined for the time window of point i

$Li$ (min): Upper bound defined for the time window of point i

$PkiS$⁠: Binary variable is equal to 1 if the ambulance k is at the ith station at the starting time

$PkiH$⁠: Binary variable is equal to 1 if the ambulance k is at the ith hospital at the starting time

$yi$⁠: Binary variable is equal to 1 if the transfer station i is established

$y'i$⁠: Binary variable is equal to 1 if the temporary hospital i is established

$xijk$⁠: Binary variable is equal to 1 if the ambulance k covers the route between points i and j

$zijh$⁠: Binary variable is equal to 1 if the helicopter h covers the route between points i and j

$δki$ (lit): Residual fuel of the kth ambulance at the time of arriving at the ith point

$tkiarrive$(min): Arrival time of the ambulance k to the ith point

$tkileft$ (min): Time of leaving the disaster-affected point i by ambulance k

$nxki$ (man): Number of red-type patients transferred from the ith disaster-affected point by the kth ambulance

$rxki$ (man): Number of red-type patients delivered by ambulance k to the transfer station i

$arki$ (man): Number of red-type patients delivered by ambulance k to the hospital i

$agki$ (man): Number of green-type patients at the area i treated by ambulance k

$x'ijh$ (man): Number of red-type patients transferred by helicopter h from the transfer point i to the hospital j

Accordingly, the deterministic mathematical model is presented as follows:

The first objective function minimizes relief time, which includes minimization of the maximum time required by the ambulance to leave the disaster-affected points. The second objective function is related to the minimization of the costs, which includes minimization of the routing costs of vehicles and locating temporary hospitals and transfer stations.

Constraint (1) states that an ambulance is at one of the stations, or it is being used in one of the hospitals at the beginning of the time horizon. Constraint (2) states that each ambulance in any station must start its trip from the same point. Constraint (3) states that each ambulance in any hospital must start its trip from the same point. Constraint (4) states that a loop cannot be generated at any point. Constraint (5) states that an ambulance should leave the point that it previously entered. Constraint (6) states that the destination of each ambulance is the transfer stations or the hospitals. Constraints (7) and (8) determine the primary required fuel for ambulances. Constraint (9) calculates the amount of fuel consumption in the routes which are covered by each ambulance. Constraint (10) specifies the refueling of ambulances at the gas station. Constraint (11) determines the arrival time to and departure time from each of the affected points for the ambulances. Constraint (12) states that all the red-type patients must be transferred from the affected areas by ambulances. Constraint (13) states that if a transfer station is established, then the red-type patients can be transferred to it. Constraint (14) states that the ambulances can only travel to the established stations. Constraint (15) specifies the number of red-type patients at each transfer station. Constraint (16) determines the number of red-type patients that are directly transferred to the hospitals by ambulances. Constraint (17) considers the departure time of each ambulance from the affected points equal to its arrival time plus the time of providing reliefs for the green-type patients. Constraint (18) states that all of the green-type patients should be treated. Constraint (19) represents the capacity limitation of ambulances. Constraint (20) states that helicopters are allowed to enter the transfer stations that were previously established. Constraint (21) states that if a helicopter enters a certain transfer station, it must leave there toward another transfer station or hospital. Constraint (22) states that only those helicopters are allowed to go to the hospital that previously stopped at one of the transfer stations. Constraint (23) represents the capacity limitation of the helicopters. Constraint (24) indicates that at each transfer station, all of the red-type patients brought by the ambulances should be transferred to the hospitals by helicopters. Constraint (25) indicates that a helicopter should land at a transfer station and then transfer the patients to the hospitals. Constraints (26) and (27) imply that the patients are only sent to the temporary hospitals which are already established. Constraint (28) states that the helicopters can only go to those temporary hospitals that are already established. Constraints (29) and (30) represent the type of decision variables.

In the proposed mathematical model, Constraints (15) and (16) cause the non-linearity of the proposed mathematical model. Using the definition of the variable, these constraints can be linearized and as a result, a complex integer programming model can be presented. Since the value of $∑i∈Ixijk$ is equal to 0 or 1, the following conversions can be made (Tirkolaee et al., 2017).

For Constraint (16), we have:

RO is one of the most efficient approaches for investigating the effects of real-world uncertainty in a problem (Tirkolaee et al., 2018, 2020; Golpîra and Tirkolaee, 2019). In this section, the robust mathematical model of the problem is presented based on the Bertsimas and Sim's (2004) model as follows:

The number of red-type patients at the ith point is indicated by $DRi$⁠. Suppose that the number of red-type at the disaster affected points is not clear and constant and might be afflicted by uncertainty based on the real-world conditions. To consider such conditions, an uncertainty interval is considered for this parameter; i.e. [$DR˜i−DRˆi,DR˜i+DRˆi$] based on the robust approach of Bertsimas and Sim (2004). Following the interval uncertainty space, each of the uncertain $DRˆi$ is in the form of a symmetric and finite distance with DRi as its center. Moreover, $DRˆi$ is the fluctuation rate of this parameter, which is obtained from the equation $DRˆi$ = ρ$DRˆi$⁠, and ρ > 0 is the level of uncertainty.

The objectives of the problem do not change and the added constraints are presented as follows:

Constraints (1) to (14)

Constraints (17) to (30)

Constraints (31) to (35)

Constraints (36) to (40).

Constraint (12) states that all of the red-type patients should be transferred from the affected areas by ambulances. After developing a robust model, this constraint transforms into Equation (41) (Bertsimas and Sim, 2004), and then Constraints (41)–(45) enter into the problem instead of Equation (12).

It should be noted that $Γi$ in Equation (41), takes a value between 0 and the number of red-type patients at the ith point, which is equal to $DRˆi$⁠. If $Γi$ is equal to $DRˆi$⁠, it indicates that an improper limit was selected for the parameter. This also which implies the highest conservatism, while a value of zero indicates the lack of control over uncertainty. The decision-maker can consider a value for it regarding the conservatism level. Furthermore, since the parameter $DRˆi$ is the only right-side vector of Constraint (12) without being multiplied by any variable, variable $ψi$ and constraint $ψi=1,∀i∈AA$ are considered conventionally to represent the structure of the robust model proposed by Bertsimas and Sim (2004).

In this section, the developed solution methods for this problem are described. Furthermore, in order to validate the methods, the validation of the solution method is done through randomly generated instances. Our proposed solution methods include the ε-constraint method and the SFLA. To solve the problem, ten breakpoints are regarded for each objective function, and accordingly, ten Pareto points are generated for each instance (Tirkolaee et al., 2019a).

The ε-constraint technique is one of the well-known methods for dealing with multi-objective problems. This solves this type of problem by transferring all of the objective functions, except for one of them, into the constraints in each step (Ehrgott and Gandibleux, 2003; Bérubé et al., 2009). One of the advantages of the ε-constraint method, compared to other exact methods such as the objectives' weighted sum method, is the ability to control the number of generated efficient solutions in the ε-constraint method, which is impossible in the weighted sum method (Mavrotas, 2009). By considering the values selected for $εi$⁠, the problem is solved.

In the proposed ε-constraint method, the first objective function is considered to be the main objective and the second objective function is considered to be the second objective. The results obtained using the ε-constraint method are presented in the “Numerical Results” section.

SFLA is the developed version of the shuffled complex evolution algorithm (SCE or SCE-UA). This is an algorithm with a relatively long background in smart optimization (Eusuff et al., 2006). SCE-UA is a combination of the evolutionary capabilities of GA with the random searching capability of the controlled random searching (CRS). Therefore, it can be classified among memetic algorithms.

By adding elitism and swarm intelligence capabilities to the SCE algorithm, the SFLA with several structural similarities is constructed. The similarity and relationship are so much that during the implementation of each of these algorithms, by applying a few changes, the computer program and implementation of each algorithm can be converted to the other.

In memetic algorithms, in contrast with GA, the attributes and capabilities are inherited to the children by parents. In memetic algorithms, each person (with respect to Lamarck's evolutionary theory) obtains useful attributes and characteristics by searching in their neighborhood (local search). This means that, in addition to the evolution of the population, the evolution also proceeds individually. Memetic algorithms are henceforth sometimes known as hybrid algorithms (HA) and local genetic algorithms (LGA).

The main steps of the algorithm in problem-solving are as follows:

• Step 0: Selecting m and n, where m and n represent the number of groups (memeplex) and the number of frogs in each group, respectively; therefore, the total population is calculated through F = mn.

• Step 1: Generating primary random possible solutions in the number of F (primary population) and calculating the objective function of each solution.

Moreover, each solution consists of the following information:

1. Number of demand points allocated to each ambulance

2. Route of each ambulance

3. Location of each ambulance at the starting moment

4. Destination of each ambulance

5. Ambulances' arrival time to demand points

6. Ambulances' departure time from demand points

7. Gas stations that each ambulance should visit

8. Established temporary hospitals

9. Established transfer points

10. Demand collected at each transfer point

11. Number of transfer points allocated to each helicopter

12. Route of each helicopter

13. Location of each helicopter at any starting moment (this is additional information and is given for ease of calculation)

14. Destination of each helicopter

15. Travel costs of each ambulance

16. Establishment costs of temporary hospitals

17. Establishment costs of transfer points

18. Travel costs of each helicopter

• Step 2: Sort the frogs based on their fitness. Store the best solution (F_g) in the total population.

• Step 3: Put it = 1 (it = 1, 2, …, N), where N represents the maximum search iterations. Divide the frogs into the groups (Mk) as follows:

For example, in case of dividing 30 frogs into six groups with five frogs, the indices 1, 7, 13, 19 and 25 are located in the first group, and the indices 4, 10, 16, 22 and 28 are located in the fourth group.

• Step 4: Frog locating algorithm (FLA):

  • Substep 4. 1: Put im = 1 (im = 1, 2, ..., m).

  • Substep 4. 2: Put in = 1 (in = 1, 2, …, b) (counter of the number of search iterations in each group).

  • Substep 4. 3: Select q frogs from among the frogs of the im group randomly with the following probability (sub-grouping); pj = 2(n + 1 - j)/n(n + 1), (j = 1, …, n).

  • Substep 4. 4: Put k = 1 (k = 1, 2, ..., q).

  • Substep 4. 5: Select an improvement algorithm (among existing algorithms) randomly with probabilities of (0.2, 0.2, 0.4 and 0.2).

  • Substep 4. 6: Implement the selected improvement algorithm on the solution of the kth frog and obtain a new solution. If the new possible solution is better than the solution of the kth frog, then replace the new solution for the worst solution in the sub-group and then go to Substep 4. 8. Otherwise (if the new solution is impossible, or if it is not better than the solution of the kth frog), select another improvement algorithm and re-implement this step.

If no new possible solution is obtained with all the 4 improvement algorithms, go to Substep 4. 7.

• Substep 4. 7: Replace a new random solution for the worst solution in the subgroup and then go to Substep 4. 8.

• Substep 4. 8: Put k = k + 1 and if k ≤ q, go back to Substep 4.4; otherwise, go to Substep 4. 9.

• Substep 4. 9: Put in = in + 1 and if in ≤ b, go back to Substep 4. 2; otherwise, go to Substep 4. 10.

• Substep 4. 10: Put im = im + 1 and if im ≤ m, go back to Substep 4.1; otherwise, go to Step 5.

• Step 5: Re-sort the existing solutions based on their fitness. Update the best solution in the total population (Fg). Set it = it + 1 and if it ≤ N, go to Step 3; otherwise, stop.

The random solutions are generated through a heuristic method as follows:

• Step 1: Generating a random formation

In this step, all the demand points are divided among the ambulances. In other words, when we have D demand points and |K| ambulances, then the total number of formations will be equal to D|K| and thus in this step, we select a formation randomly.

For example, to divide 30 demand points among 4 ambulances, the random formations of 8, 2, 15 and 5 is obtained in this step.

• Step 2: Allocating the demand points to ambulances randomly based on the formations of Step 1 and determining the route of each ambulance

In this step, first, the demand points are sorted randomly. Then, they are allocated to the ambulances based on the formations of step 1. For the previous example, after random sorting of the demand points, the first eight points are allocated to ambulance 1, 2, next points to ambulance 2, 15, next points to ambulance three and, the last 5 demand points to ambulance 4.

• Step 3: Checking the capacity of the ambulances

In the previous step, the route of each ambulance was obtained. Consequently, the capacity of the ambulances should be investigated. If the capacity of an ambulance allows going to the allocated demand points, we go to the next step. Otherwise, we return to the previous step.

• Step 4: Checking the time window constraint of demand points and calculating the time and cost

In this step, by calculating the arrival time of each ambulance to each demand point as well as its departure time from that point, the time window constraint is investigated. Also, residual fuel of the ambulances is taken into account and if necessary, the time of going to the gas station for refueling will be taken into consideration as well.

Furthermore, in this step, the primary location of the ambulances (an ambulance station or hospital) is also calculated and determined; moreover, the destination of the ambulances is determined randomly in this step.

In this step, if the time window constraint is not met for a demand point, we go back to Step 1. Otherwise, we will go to Step 5.

• Step 5: Generating the route of helicopters

In the previous step, the destination of the ambulances was determined where some of the ambulances delivered their patients to the transfer points; therefore, the route of helicopters should be determined to take these patients to the (fixed or temporary) hospitals. For this purpose, firstly, the demand collected at the transfer points is randomly allocated to the helicopters regarding their capacity. This procedure is performed almost similar to the solution of the transportation problem, but randomly, not based on the transportation cost. After allocating the demand collected at the transfer points among the helicopters, the cost of traveling through the routes (from the primary location of the helicopter to the transfer point, between transfer points, from transfer points to fix or temporary hospitals) is also calculated.

Two ambulances are selected randomly. In the primary random formation, a demand point is taken from an ambulance and then it is allocated to another ambulance. Then, this change is applied on the route of those two ambulances as well.

Such changes might result in altering the capacity constraint of these two ambulances; thus, the capacity constraint is re-checked for these two ambulances.

Moreover, the arrival and departure time from the demand points, as well as the traveling costs, are re-calculated for these two ambulances; moreover, the time window constraint is re-investigated.

Finally, helicopters are re-routed as well.

An ambulance is selected randomly and then two or more demand points are displaced in the route. For example, if the route of the selected ambulance is:

8, 14, 5 and 15

Random displacement of 3 points: The points are displaced with random order of 1, 4, 3.

1. Point 3 is replaced by point 1.

2. Point 4 is replaced by point 3.

3. Point 1 is replaced by point 4.

5, 14, 15 and 8

After applying this change, the arrival and departure time are recalculated for this ambulance as well as the costs; furthermore, since no change is made in the destination of the ambulances, there is no need for re-routing the helicopters.

In this algorithm, the destination of all ambulances is determined randomly again and the previous destinations are eliminated. Then, due to the probable changes in the amount of the demand collected at the transfer points, the helicopters are re-routed as well.

In this algorithm, the demand collected at the transfer points is randomly divided among the helicopters and the helicopters are re-routed.

Flowchart of the proposed SFLA algorithm is shown in Figure 2.

Srinivas and Deb (1994) introduced the NSGA. They divided the evolutionary group into several levels according to a dominance relation for selection and solution. An operational NSGA scale was then optimized by Deb et al. (2002) to apply the elite mechanism instead of sharing coefficient of the density function. This algorithm is called NSGA-II. The suggested pseudo-code of NSGA-II in this research is given in Figure 3.

Solution representation and the main mechanism of the proposed NSGA-II are discussed in the following.

In the proposed algorithms, to demonstrate an initial solution, a six-row matrix is used, such that all of its entries are 0 and 1. Figure 4 shows the defined string for the GA.

In Figure 4, in the first row of the matrix, there are |K|*|I | columns and if the ambulance k is at the ith station at the starting moment, its relevant entry in this row takes the value of 1. In the second row of the matrix, there are |K|*|I| columns and if the ambulance k is at the ith hospital at the starting moment, then its relevant entry in this row takes the value of 1.

In the third row of the matrix, there are |I| columns and if the transfer station i is established, then the ith entry of this row takes the value of 1. The fourth row of the matrix has I columns and if the temporary hospital j is established, then the ith entry of this row takes the value of 1. The fifth row of the matrix has |I|*|J|*|K| columns and the ijkth entry of this row is equal to 1 when the ambulance k travels the route between i and j. Similarly, the fifth row of the matrix has |I|*|J|*|H| columns and the ijhth entry of this row is equal to 1 when the helicopter h travels the route between i and j.

In addition, for crossover operator in the NSGA-II, a two-point crossover operator is used, while for mutation operator, a single-point mutation operator is used.

In this section, in order to analyze the validation of the presented mathematical model, several instances of the problem are generated randomly in different sizes. The information of each instance is as shown in Table 1. In order to generate instances in a small size (1st instance), medium size (2nd and 3rd instances) and large size (4th instance) 2D space, n nodes are considered such as the nodes of the fixed hospitals, nodes of potential temporary hospitals, nodes of transfer stations, nodes of ambulance stations, nodes of helicopter stations, demand nodes and nodes of gas station points. By considering the different values for them, the instances are generated. In fact, different sized instances can be generated considering different sizes for the relief network (see Table 1). In each random instance, a unique uniform distribution is used which is shown in Table 2. Most of these values are adapted from previous research such as that by Najafi et al. (2013). Different values can yield different results and the impact of the parameters' change can be studied using sensitivity analysis. In each random instance, the value of $Γi$ is equal to half of the total number of red-type injured people rounded upwardly and calculations are performed for the uncertainty levels of 0.2, 0.4 and 0.5. These values are selected because the maximum uncertainty level of 50% is probable in relief logistics (Zokaee et al., 2016).

To solve the ε-constraint solution method, the GAMS (V-24.1) software was used and the algorithm was encoded using the MATLAB (V-16.0) software. The problem is solved separately in order to determine the Pareto points based on each objective function. It should be noted that all runs of the problem in the GAMS software are performed considering a time limitation of 3,600 s (Bozorgi-Amiri et al., 2011). This time limitation is only set to test the performance of the proposed algorithms in small-sized instances. In fact, after executing the 3600-s time limitation, the best-obtained solution is reported which might be the overall optimal or local optimal solution. It should be noted that 4th instance in the ε-constraint method does not yield a solution within this time limitation. This instance can be solved only by the SFLA algorithm.

In Table 1, from left to right, the columns represent the rank of the generated random instances, the number of nodes in the network, the number of available ambulances, the number of available helicopters, the number of disaster-affected points within the network, the number of fixed hospitals, the number of potential temporary hospitals, the number of transfer stations, the number of ambulance stations, the number of helicopter stations and the number of gas stations respectively.

The results are reported for various levels of uncertainty in Tables 3–5 respectively. In these tables, the value reported for each objective function is obtained by solving the instances with that specific single objective. The gap reported in Tables 3–5 is calculated as follows:

Then, with regard to the ε-constraint method for different levels of uncertainty, the values of the epsilons are determined for the second objective, which are presented in Table 6 for the different breakpoints. In Table 6, each instance is solved for the ten breakpoints and three levels of uncertainty so then the values of $ε2$ are obtained.

As seen in Table 6, the uncertainty level has a direct impact on the value of the objectives. It demonstrates that higher uncertainty levels lead to higher (worse) values for the objectives by keeping in mind the feasibility of the solutions. The main reason is the presence of uncertainty in the number of injured people which requires more resources and facilities for transportation planning. Consequently, the service time increases as well.

The Pareto fronts of the proposed 1st-4th instances are presented in Figures 5–8 for the uncertainty level of 0.5. The 1st and 2nd instances are solved using the exact method, while the 3rd and 4th instances are solved using the SFLA algorithm.

According to Figures 5–8, different instances have different Pareto fronts with different dispersion rates. The obtained trade-offs indicate that the objectives interact reversely. In other words, by increasing one of them, the other decreases. Moreover, the dispersion of the Pareto fronts in each instance is different. To study the trade-offs and to choose the best values of these points, the reference point approach is implemented.

In this subsection, an efficient approach is introduced to choose the best Pareto solutions. The selected Pareto solution can be investigated and implemented as part of the optimal planning of the relief logistics system. Deb and Sundar (2006) presented a reference point approach in the context of multi-objective problems in order to determine the best point or Pareto point. This method can be applied by either weighing the objectives or not. The main idea of this method is to highlight the solutions that are close to the reference points.

The normalized Euclidean distance (dev_j) between the solutions of the best non-dominated level and the reference point is calculated using Equation (48). Each solution with a lower value of dev_j has a higher priority for the decision-maker.

In Equation (47), $wi$ is the weight of objective i in each solution, $fi$ is the value of objective i, $z¯i$ is the value of reference objective, and $fimax$ and $fimin$ are the maximum and minimum values for objective i respectively. The values of the parameters related to the determination of the best Pareto solution are presented in Table 7.

The results obtained from solving three instances are presented as the best Pareto solution at the break points and these results are compared with those of the meta-heuristic algorithm. Tables 8 and 9 represent the results obtained from both the ε-constraint method and the SFLA algorithm. In these tables, the values of the first and second objectives as well as the average solution time are presented for the different levels of uncertainty.

The results obtained from solving the problem through the exact method are presented in Figure 9 in different sizes for the different levels of uncertainty.

As shown in Figure 9, the value of the second objective increases from the small to large instances, meaning that increasing the size of the problem leads to the increase of the second objective. Since the first objective represents the minimization of the ambulances' maximum departure time from the affected points, it may be reduced by using a higher number of ambulances. However, there is a trade-off between these total cost increases and the number of injured people. The results of the SFLA algorithm are presented as follows.

The results obtained from solving the problem using the SFLA algorithm in different size for different levels of uncertainty are presented in Figure 10.

After reviewing the results in Tables 8 and 9 as well as in Figures 9 and 10, the results of the exact method and SFLA algorithm are close to each other and the proposed algorithm has a high performance. Furthermore, it yields better results for the uncertainty levels of 0.2 and 0.4 in the 3rd instance than in the applied exact method. The results of the NSGA-II algorithm are presented as follows.

According to Figure 11 and Table 10, the results of the NSGA-II algorithm are highly similar to the exact method and SFLA algorithm. Nevertheless, in order to specify the proximity of our meta-heuristic algorithms to the exact method and to evaluate their performance, several well-known indicators of the multi-objective algorithms are calculated in the following.

In this section, the results are compared with those of the NSGA-II and SFLA algorithms and then the algorithms were validated. In fact, results of the NSGA-II and SFLA algorithms in instances with small and medium sizes are compared to the exact method by the CPLEX solver of GAMS software (Tirkolaee et al., 2020). Since it is not logical to solve the large-sized instances with exact methods, they were solved using meta-heuristic algorithms.

For a more accurate validation of the proposed algorithms and also to compare their performance, the indicators of the multi-objective algorithms were used. For this purpose, the four indicators of dispersion measurement (DM), mean ideal distance (MID), space measure (SM) and the number of Pareto solutions (NPS) were calculated and then the performance of the proposed algorithms was assessed. Higher values of DM and NPS show the higher ability of the algorithm to find the Pareto front. However, this is true for MID and SM with their lower values (Sadeghi and Niaki, 2015).

As shown in Tables A1–A3 ( Appendix), the meta-heuristic algorithms used in this research perform very closely to the exact algorithm and their performance charts are very similar. Thus they are appropriately efficient for finding near-optimal solutions. Accordingly, due to such an appropriate performance, they can be used to solve the large-sized instances that cannot be solved by the ε-constraint method.

More specifically, in terms of the diversity criterion, the SFLA algorithm has a much better performance than the NSGA-II algorithm. This means that the SFLA algorithm is capable of generating Pareto fronts that include a wider range of possible solutions. However, the SFLA algorithm has functional fluctuations in terms of the diversity criterion. Thus, it does not have a better performance than the NSGA-II algorithm in the 1st instance with an uncertainty level of 0.5. It works so much better in the other instances. On this basis, disregarding the exceptions, it can be concluded that the SFLA algorithm exhibits a higher performance in this index.

In the criterion of distance from an ideal point, the SFLA algorithm has a better performance. This means that in this algorithm, the distance from the ideal solution is less than that of the other algorithms. In contrast to the diversity index, there is no exception in the index of distance from the ideal point and the SFLA algorithm has a better performance.

Furthermore, based on the distance criterion, the SFLA algorithm represents a better performance and it generates more approximate solutions. Thus it can be said that it has a more uniform front compared to the NSGA-II algorithm. Besides, the SFLA algorithm always succeeds to find more Pareto solutions. In fact, the main idea of implementing the NSGA-II algorithm is to test the performance of the SFLA algorithm in large-sized instances.

For a general comparison of the algorithms, the simple additive weighting (SAW) method is used (Zitzler et al., 2000). Accordingly, all of the values obtained for the indicators are normalized first. Then for each algorithm in each problem, the average of the indicators is calculated (with the same weight) and then considered as the performance of that algorithm in the given problem. These comparisons are shown in Figures 12–14 for the uncertainty levels of 0.2, 0.4 and 0.5. Regarding Figures 12–14, the NSGA-II and SFLA algorithms perform similarly to the optimal method to a large extent.

However, to solve the problem in a large size (4th instance), the proposed meta-heuristic algorithms were used because the ε-constraint method cannot provide accurate solutions for such problems. To determine which algorithm has a better performance in large-sized instances, the values of the indicators are investigated in Figures 12–14. These figures reveal that the SFLA algorithm has a much better performance in general.

Figure 15 also shows the solution time for different instances. According to Figure 14, increasing the dimensions of the problem leads to a significant increase in the solution time of the exact method. As a result, in the 4th instance, the ε-constraint method could not solve the problem in the given time limitation. On the other hand, the meta-heuristic algorithms were capable of solving the problem in a much shorter time. The NSGA-II algorithm requires less time to find the Pareto solutions compared to the SFLA algorithm.

This subsection provides an analytical comparison between the results obtained by the deterministic and robust problems. In fact, as a great managerial implication, it would be useful to study the effects of the uncertainty on the problem which may even make the problem infeasible due to resource shortages. To achieve a better perception of the effect of uncertain parameters and concerning the use of the RO approach to solve the problem, the results of the problems in deterministic conditions are compared the robust problem at an uncertainty level of 0.5 using the SFLA algorithm in Table 11.

According to Table 11, the first objective of the robust model is lower than the first objective in all of the deterministic instances, which is one of the benefits of RO. This indicates that in robust instances, the focus is on the reduction of departure maximum time from the given points in total, which can only be achieved by spending higher costs. In contrast to the first objective, the second objective in all of the robust instances is much more than the second objective of the deterministic instances, implying that having a robust problem and solution requires considering higher resources and subsequently, a higher total cost. Such an increase in cost is observable in the last four columns of Table 11.

One of the major managerial achievements is the special attention paid to the uncertain conditions in the real world, which would lead to the changes in the solution space as well as the results. Robustness also has a direct impact on the service-providing facilities and equipment in the case of disasters and eventually the total costs. Most importantly, this impacts on the response time. This can be associated with increasing the costs in order to impose a subsequent reduction in response time.

This subsection provides a quick reference guide to the different statuses. To this end, Table 12 is presented as a summary table showing the superior performance of each solution method.

According to Table 12, it is concluded that the $ε$-constraint method is regarded as the best solver for small- and medium-sized instances and that SFLA is the best for large-sized instances.

In this subsection, the impacts of changing the transferring times; i.e. $tij$ on the objective functions is studied using sensitivity analysis. In other words, $tij$ can be regarded as a determinant parameter from the perspective of the operations managers. Hence, management wants to know how much this change may have a positive/negative effect on the objective functions. To this end, the fluctuation range of −20%, −10%, 0%, +10 and 20% are taken into account for $tij$⁠. Then the obtained values of both objective functions were evaluated for the robust problem using the SFLA algorithm. The 4th instance was considered to conduct this sensitivity analysis. The results are represented in Table 13 and Figures 16–17.