# Multigraph Modeling for Urban Distribution of Emergency Commodities with Semisoft Time Windows under Uncertainty

Authors: Hamid Tikani; Mostafa Setak; Darya Abbasi
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

In this paper, we studied a stochastic bi-objective mathematical model for effective and reliable rescue operations in multigraph network. The problem is addressed by a two-stage stochastic nonlinear mixed-integer program where the reliability of routes is explicitly traded-off with total weighted completion time. The underlying transportation network is able to keep a group of multiattribute parallel arcs between every pair of nodes. By this, the proposed model should consider the routing decision in logistic planning along with the path selection in an uncertain condition. The first stage of the model concerns with the vehicle routing decisions which is not involved with random parameters; besides, the second stage of the model involves with the departure time at each demand node and path finding decisions after observation of random vectors in the first stage considering a finite number of scenarios. To efficiently solve the presented model, an enhanced nondominated sorting genetic algorithm II (NSGA-II) is proposed. The effectiveness of the introduced method is then evaluated by conducting several numerical examples. The results implied the high performance of our method in comparison to the standard NSGA-II. In further analyses, we investigated the beneficiary of using multigraph setting and showed the applicability of the proposed model using a real transportation case.


---

## Body

## 1. Introduction

Natural and human-caused disasters are unforeseen and unexpected events, which would cause catastrophic loss of life, physical destruction, and massive economic disruption. After any emergency events, the demands for medical equipment, sanitation medicine, water, food, tents, blankets, and other relief commodities are expected to highly increase. Emergency Events Database (EM-DAT) was launched by the Center for Research in 1988 on the Epidemiology of Disasters. EM-DAT is a global database on disasters containing essential core data on the occurrence and impacts of mass technological and natural disasters from 1900 to 2015 [1]. Figure 1 shows that in the year 2015, natural disasters had a devastating impact on human societies. In the year 2015, 376 reported natural disasters caused the fatality of 22,765 people, made 110.3 million victims, and caused US $ 70.3 billion damages [1]. According to the statistics, designing an effective logistic network for rescue operations in the early post-disaster situations can significantly reduce the number of victims.Figure 1
Trends in occurrence and victims of disasters [1].In the traditional vehicle routing problem, the shortest route between nodes is taken into consideration as the distance between two different nodes. Considering only one arc between each pair of nodes may result in overlooking some solutions from the solution region. In the real world, many transportation networks are multigraph which means it is possible to reach the destination from an origin by different paths with multiple features (cost, reliability, time, and so on) and conditions (such as traffic jams, weather conditions, or accidents). Besides, disasters can seriously damage infrastructures leading to destroy the distribution network. Therefore, considering the availability of alternative paths in crisis management is a challenging but important issue. The reason is that considering parallel edges in the transportation network can substantially enhance the achieved solutions’ quality and can accelerate the emergency response in disaster situations. The current study extends a bi-objective optimization model based on Tikani and Setak [2] for the vehicle routing and distribution problem for relief operations with the presence of multiple links between vertexes where the links are differentiated by two independent attributes including different travel time and reliability.Crises are often unpredictable; as a result, one of the most important issues which relief organizations are confronting is the high degree of uncertainty intrinsic in disasters. The mathematical model proposed by Tikani and Setak [2] does not reflect uncertainty in the transportation of post-disaster. Thus, in order to allow the model approache the reality, we developed two-stage mathematical models in an uncertain condition. The uncertainty in the reliability of links, demands of affected areas, and congestion of roads are captured by a finite number of possible scenarios.Some disasters impact the transportation infrastructures and links. Reliable logistic planning in post-disaster holds a pivotal role in ensuring the timely transferring of relief goods to disaster areas and protecting the rescue team [3]. In this study, the reliability of the route is taken into account as an objective that scarcely has been addressed in the literature. We studied the problem in a multigraph network under uncertain conditions, and we have more than one connection between nodes with various traffic patterns and reliabilities. Considering reliability in objective along with other objectives will lead to more nondominated solutions which give managers more confidence to make decisions.The remainder of this study is organized as follows. The related literature and innovations of this paper are reviewed in Section2. In Section 3, the problem description and mathematical model are presented. In Section 4, the solution method is introduced. In Section 5, computational results and case study are presented. Finally, the conclusion is stated in Section 6.


## 2. Literature Review

Recently, emergency logistics management has received considerable attention from researchers. In this section, the literature review is split into four subsections. In the first part, the literature relevant to the post-disaster logistics planning and reliability of distribution is discussed. In the second subsection, we review the works relevant to the vehicle routing problem with semisoft time windows. In the third part of this section, the studies on routing problems in multigraph networks are reviewed. The fourth subsection focuses on works considering uncertainty in post-disaster. Then, the contributions of the paper are presented.


### 2.1. Related Works on Reliable Distribution in Post-Disaster

In post-disaster conditions, critical infrastructures might fail. Therefore, in emergency logistics management, the expected reliability of selected roads is one of the vital items which should be considered as it holds an important role in the timely and safe delivery of required facilities. Although designing a reliable and safe dispatching process is important, only a few works in the literature considered the concept of reliability in the emergency dispatching. In this regard, Vitoriano et al. [4] developed the multiobjective mathematical model for relief operation problem that takes equity, time, priority, cost, security, and reliability into consideration. A goal programming method was employed to cope with different objectives. In another study, Hamedi et al. [5] proposed a model considering scheduling of supplies and reliable routing together over a time-dependent network. A three-objective integer nonlinear programming model was presented by Wang et al. [3] for the relief distribution problem. The objective functions were total cost, links reliability, and travel time. Liberatore et al. [6] proposed a distribution model considering practical factors such as security and reliability. A multiobjective model was developed by Vahdani et al. [7]. They considered vehicle routing, facility location, and emergency roadway repair operations in their model. Zhou et al. [8] investigated resource distribution in post-disaster by presenting a multiobjective mathematical model. They considered roadway availability in a multiperiod dynamic emergency resource-scheduling problem. Tikani and Setak [2] formulated a reliable time-dependent vehicle routing problem with hard time windows in a multigraph and developed metaheuristics to minimize delays in transferring essential goods in a disaster.


### 2.2. Related Works on Uncertainty in Post-Disaster

Disasters are generally characterized by a high level of uncertainty. Various parameters such as time and demands are hard to predict; hence, decisions should be taken in an uncertain environment. Two-stage stochastic programming is one of the most well-known modeling methods which supports decision-making under uncertainty. Bozorgi-Amiri et al. [9] proposed a multiobjective robust stochastic model for relief logistic under uncertainty. They took demand, purchase, supply, and transportation costs into account as uncertain parameters. They developed their model based on two stages. In the first stage, the distribution center locations and required inventory of any relief goods under storage were determined. In the second stage, the level of items delivered from relief centers to disaster regions was specified. The proposed model was in accord with the assumption that disaster information is not dependent on routing the vehicles and time. Bozorgi-Amiri and Khorsi [10] developed a multiobjective dynamic location-routing model for midplanning and short-term regarding relief under uncertainty. Travel, cost, time, and demand were considered as uncertain parameters. The aims of their model were minimization of the maximum shortage among the disaster regions in all periods, total cost, and travel time. The ɛ-constraint method was employed to solve the proposed model. Tofighi et al. [11] formulated a two-echelon relief logistics network design problem with determining the number of local distribution centers and central warehouses. Also, a relief distribution program for various crisis scenarios was proposed. They developed a new two-stage scenario-based possibilistic-stochastic model. A metaheuristic algorithm was employed to solve the model. A robust model for a three-level relief chain was proposed by Zokaee et al. [12]. They considered relief distribution centers, suppliers, and the disaster regions in an uncertain environment. The model was employed to increase the satisfaction of people by minimizing the shortage of relief goods and minimizing the relief operations’ costs. The capability of their method was verified by a case study. Vahdani et al. [13] presented a two-stage, mixed-integer, multiobjective, multicommodity, and multiperiod model in the three-level relief chain. In the first phase, locations of the distribution warehouses and centers with different capacity levels were taken into consideration. In the second phase, they took the restricted hard time windows into account. They also considered the reliability of the route in their study and employed two metaheuristic algorithms to solve the proposed problem. Li et al. [14] proposed a bi-objective rescue routing model for urban emergency logistics based on travel time reliability in a stochastic network. In this study, road travel time was assumed to be strongly influenced by urban storm flooding disasters. To solve the model, a hybrid metaheuristic that integrated Tabu Search and ant colony optimization was designed. It is worth mentioning that their model aims to find the best alternative path with time-independent travel times from a warehouse to only one disaster point considering that a vehicle can provide rescue services to only one affected area.


### 2.3. Vehicle Routing and Scheduling Problem with Semisoft Time Windows

In resource distribution in post-disaster, it is momentous to optimize the delivery plans to improve the time of emergency packages, reducing distribution costs and increasing the packages safety of delivering process. Vehicle routing problem with semisoft time windows is an extension of vehicle routing problem, including the vehicle routing problem with hard and soft time windows. In vehicle routing problem with soft time windows, maximum and minimum limits of time window are possible to be violated at a penalty. A problem that only takes penalties on late arrival into account is named the vehicle routing and scheduling problem with semisoft time windows [15].Qureshi et al. [15, 16] proposed a novel column generation according to the exact approach for the scheduling and vehicle routing problem with semisoft time windows. In their study, an upper bound on the lateness time window deviation was considered. The proposed solution algorithm is proved to be effective in solving medium-sized test problems. Tas et al. [17] formulated a vehicle routing problem with stochastic travel times and soft time windows. They considered service and transportation costs in their model. A Tabu Search algorithm was used as a solution method. Beheshti and Hejazi [18] presented a mathematical model for the vehicle routing problem with general soft time window. To solve the model, a hybrid column generation-metaheuristic algorithm was proposed. Bhusiri et al. [19] studied the vehicle routing problem with soft time windows. They applied the branch-and-price method which led to a set partitioning master problem and its new subproblem. They presented the new methods to solve the subproblem.


### 2.4. Vehicle Routing Problem in Multigraph Networks

Prior studies in the vehicle routing problem assume that only one edge exists among two sequential nodes. Owing to the intricacy of urban environments, it is likely to exist more than one route among two nodes. Neglecting the available parallel arcs in the transportation network in such problems may adversely affect the quality of obtained solutions. Hence, this concept is important to be taken into account by decision makers. A few works in the literature considered multigraph network in vehicle routing problem. Garaix et al. [20] focused on the routing problem with the existence of alternative paths through a multigraph setting. Lai et al. [21] developed the heterogeneous vehicle routing model in a multigraph network. Alternative arcs connecting customers were characterized by various features (cost and time). Ticha et al. [22] for solving the vehicle routing problem with time windows in single and multigraph networks presented two branch-and-price algorithms. However, only two studies focused on emergency logistics planning in a multigraph network including Tikani and Setak [2] and Setak et al. [23]. Setak et al. [23] modeled a time-dependent vehicle routing problem in a multigraph-based network for emergency logistics operations in the natural disaster occurrences. Also, they considered FIFO property in their model and concentrated on the minimization of total number of vehicles and waiting time. Tikani and Setak [2] investigated the availability of parallel links between nodes considering multiple attributes including different reliabilities and travel times. They proposed a nonlinear mathematical model which minimized the total weighted completion time where the reliability of routes was controlled to be more than a predetermined threshold value. Two other works which focused on routing application in multigraphs are Tikani and Setak [24] and Tikani and Setak [25].In the current study, we extend the model of Tikani and Setak [2] from several aspects and study a stochastic reliable time-dependent relief routing problem in multigraphs with semisoft time windows (abbreviated by SR-TDRRPM-SSTW) for early post-disaster operations. The main contributions of the current study are as follows:(i)
This is the first study that investigates the availability of alternative paths in urban emergency logistics under uncertain conditions. To do this, we applied two-stage stochastic programming where the demands of affected areas, reliability of links, traffic congestions, and service time duration are considered as a finite number of possible realizations.(ii)
Due to the uncertainty of vehicles’ speeds aftermath of a disaster, it is possible that the emergency packages do not deliver to the victims and wounded people after a desirable time. To handle this issue, we partially relaxed the time window and penalized a time-dependent late arrival cost. In detail, the semisoft time window is incorporated to allow and also control the delays in servicing demand nodes.(iii)
The problem is studied under multiple objectives. The first objective minimizes the total transportation costs including operating costs and late arrival penalties, while the second one maximizes the reliability of constructed routes.(iv)
Due to the nonlinearity and NP-hardness of the problem, it is computationally difficult to solve the model. For this purpose, we present a modified multiobjective metaheuristic based on NSGA-II to tackle the problem. We compare the performance of modified NSGA-II with the standard NSGA-II to prove the efficiency of the developed approach.


## 2.1. Related Works on Reliable Distribution in Post-Disaster

In post-disaster conditions, critical infrastructures might fail. Therefore, in emergency logistics management, the expected reliability of selected roads is one of the vital items which should be considered as it holds an important role in the timely and safe delivery of required facilities. Although designing a reliable and safe dispatching process is important, only a few works in the literature considered the concept of reliability in the emergency dispatching. In this regard, Vitoriano et al. [4] developed the multiobjective mathematical model for relief operation problem that takes equity, time, priority, cost, security, and reliability into consideration. A goal programming method was employed to cope with different objectives. In another study, Hamedi et al. [5] proposed a model considering scheduling of supplies and reliable routing together over a time-dependent network. A three-objective integer nonlinear programming model was presented by Wang et al. [3] for the relief distribution problem. The objective functions were total cost, links reliability, and travel time. Liberatore et al. [6] proposed a distribution model considering practical factors such as security and reliability. A multiobjective model was developed by Vahdani et al. [7]. They considered vehicle routing, facility location, and emergency roadway repair operations in their model. Zhou et al. [8] investigated resource distribution in post-disaster by presenting a multiobjective mathematical model. They considered roadway availability in a multiperiod dynamic emergency resource-scheduling problem. Tikani and Setak [2] formulated a reliable time-dependent vehicle routing problem with hard time windows in a multigraph and developed metaheuristics to minimize delays in transferring essential goods in a disaster.


## 2.2. Related Works on Uncertainty in Post-Disaster

Disasters are generally characterized by a high level of uncertainty. Various parameters such as time and demands are hard to predict; hence, decisions should be taken in an uncertain environment. Two-stage stochastic programming is one of the most well-known modeling methods which supports decision-making under uncertainty. Bozorgi-Amiri et al. [9] proposed a multiobjective robust stochastic model for relief logistic under uncertainty. They took demand, purchase, supply, and transportation costs into account as uncertain parameters. They developed their model based on two stages. In the first stage, the distribution center locations and required inventory of any relief goods under storage were determined. In the second stage, the level of items delivered from relief centers to disaster regions was specified. The proposed model was in accord with the assumption that disaster information is not dependent on routing the vehicles and time. Bozorgi-Amiri and Khorsi [10] developed a multiobjective dynamic location-routing model for midplanning and short-term regarding relief under uncertainty. Travel, cost, time, and demand were considered as uncertain parameters. The aims of their model were minimization of the maximum shortage among the disaster regions in all periods, total cost, and travel time. The ɛ-constraint method was employed to solve the proposed model. Tofighi et al. [11] formulated a two-echelon relief logistics network design problem with determining the number of local distribution centers and central warehouses. Also, a relief distribution program for various crisis scenarios was proposed. They developed a new two-stage scenario-based possibilistic-stochastic model. A metaheuristic algorithm was employed to solve the model. A robust model for a three-level relief chain was proposed by Zokaee et al. [12]. They considered relief distribution centers, suppliers, and the disaster regions in an uncertain environment. The model was employed to increase the satisfaction of people by minimizing the shortage of relief goods and minimizing the relief operations’ costs. The capability of their method was verified by a case study. Vahdani et al. [13] presented a two-stage, mixed-integer, multiobjective, multicommodity, and multiperiod model in the three-level relief chain. In the first phase, locations of the distribution warehouses and centers with different capacity levels were taken into consideration. In the second phase, they took the restricted hard time windows into account. They also considered the reliability of the route in their study and employed two metaheuristic algorithms to solve the proposed problem. Li et al. [14] proposed a bi-objective rescue routing model for urban emergency logistics based on travel time reliability in a stochastic network. In this study, road travel time was assumed to be strongly influenced by urban storm flooding disasters. To solve the model, a hybrid metaheuristic that integrated Tabu Search and ant colony optimization was designed. It is worth mentioning that their model aims to find the best alternative path with time-independent travel times from a warehouse to only one disaster point considering that a vehicle can provide rescue services to only one affected area.


## 2.3. Vehicle Routing and Scheduling Problem with Semisoft Time Windows

In resource distribution in post-disaster, it is momentous to optimize the delivery plans to improve the time of emergency packages, reducing distribution costs and increasing the packages safety of delivering process. Vehicle routing problem with semisoft time windows is an extension of vehicle routing problem, including the vehicle routing problem with hard and soft time windows. In vehicle routing problem with soft time windows, maximum and minimum limits of time window are possible to be violated at a penalty. A problem that only takes penalties on late arrival into account is named the vehicle routing and scheduling problem with semisoft time windows [15].Qureshi et al. [15, 16] proposed a novel column generation according to the exact approach for the scheduling and vehicle routing problem with semisoft time windows. In their study, an upper bound on the lateness time window deviation was considered. The proposed solution algorithm is proved to be effective in solving medium-sized test problems. Tas et al. [17] formulated a vehicle routing problem with stochastic travel times and soft time windows. They considered service and transportation costs in their model. A Tabu Search algorithm was used as a solution method. Beheshti and Hejazi [18] presented a mathematical model for the vehicle routing problem with general soft time window. To solve the model, a hybrid column generation-metaheuristic algorithm was proposed. Bhusiri et al. [19] studied the vehicle routing problem with soft time windows. They applied the branch-and-price method which led to a set partitioning master problem and its new subproblem. They presented the new methods to solve the subproblem.


## 2.4. Vehicle Routing Problem in Multigraph Networks

Prior studies in the vehicle routing problem assume that only one edge exists among two sequential nodes. Owing to the intricacy of urban environments, it is likely to exist more than one route among two nodes. Neglecting the available parallel arcs in the transportation network in such problems may adversely affect the quality of obtained solutions. Hence, this concept is important to be taken into account by decision makers. A few works in the literature considered multigraph network in vehicle routing problem. Garaix et al. [20] focused on the routing problem with the existence of alternative paths through a multigraph setting. Lai et al. [21] developed the heterogeneous vehicle routing model in a multigraph network. Alternative arcs connecting customers were characterized by various features (cost and time). Ticha et al. [22] for solving the vehicle routing problem with time windows in single and multigraph networks presented two branch-and-price algorithms. However, only two studies focused on emergency logistics planning in a multigraph network including Tikani and Setak [2] and Setak et al. [23]. Setak et al. [23] modeled a time-dependent vehicle routing problem in a multigraph-based network for emergency logistics operations in the natural disaster occurrences. Also, they considered FIFO property in their model and concentrated on the minimization of total number of vehicles and waiting time. Tikani and Setak [2] investigated the availability of parallel links between nodes considering multiple attributes including different reliabilities and travel times. They proposed a nonlinear mathematical model which minimized the total weighted completion time where the reliability of routes was controlled to be more than a predetermined threshold value. Two other works which focused on routing application in multigraphs are Tikani and Setak [24] and Tikani and Setak [25].In the current study, we extend the model of Tikani and Setak [2] from several aspects and study a stochastic reliable time-dependent relief routing problem in multigraphs with semisoft time windows (abbreviated by SR-TDRRPM-SSTW) for early post-disaster operations. The main contributions of the current study are as follows:(i)
This is the first study that investigates the availability of alternative paths in urban emergency logistics under uncertain conditions. To do this, we applied two-stage stochastic programming where the demands of affected areas, reliability of links, traffic congestions, and service time duration are considered as a finite number of possible realizations.(ii)
Due to the uncertainty of vehicles’ speeds aftermath of a disaster, it is possible that the emergency packages do not deliver to the victims and wounded people after a desirable time. To handle this issue, we partially relaxed the time window and penalized a time-dependent late arrival cost. In detail, the semisoft time window is incorporated to allow and also control the delays in servicing demand nodes.(iii)
The problem is studied under multiple objectives. The first objective minimizes the total transportation costs including operating costs and late arrival penalties, while the second one maximizes the reliability of constructed routes.(iv)
Due to the nonlinearity and NP-hardness of the problem, it is computationally difficult to solve the model. For this purpose, we present a modified multiobjective metaheuristic based on NSGA-II to tackle the problem. We compare the performance of modified NSGA-II with the standard NSGA-II to prove the efficiency of the developed approach.


## 3. Problem Description and Modeling

Distribution of emergency relief after the occurrence of a disaster has a significant impact on reducing casualties. Due to the uncertain situations that occur in such a complex situation, we focus on the stochastic version of the problem to achieve reliable and prompt distribution of urgent commodities in inner-city areas. The proposed SR-TDRRPM-STW model is developed based on a time-dependent vehicle routing problem and presented in a special transportation network named multigraph networkMG=V,E where the set of vertexes V includes one central depot (DC) along with a set of demand nodes N=1,…,n. Moreover, the available connections between the nodes are presented by E=i,j,m:i,j∈ν,i≠j which shows the link between nodes i and j. Each link has two attributes including different traffic patterns and reliability. In the following, we discuss the properties of SR-TDRR PM-STW.


### 3.1. FIFO (First-In, First-Out) Property Satisfaction

The nonpassing or the first-in, first-out (FIFO) property guarantees that between two identical vehiclesA and B, if vehicle A begins to traverse the arc i,j earlier, then it arrives at its destination j before vehicle B. Primitive works on time-dependent vehicle routing problem did not satisfy FIFO property [26]. In fact, they utilized a discrete function of time to model the travel times. An example of their approach is illustrated in Figure 2 [26]. As it can be seen from the figure if a vehicle departs the origin at 4 : 30, it will arrive at its destination in 150 miles away at 9 : 30. Nonetheless, if a vehicle departs at 6 : 00, it will arrive to its destination at 9 : 00. It is clear that this situation is not feasible in the real word and the result does not conform the FIFO property.Figure 2
An example of ‘‘passing” with discrete travel time function.Recently, scholars employ continuous travel time functions over the time horizon. In this paper, we employed the approach presented by Setak et al. [27] to convert a travel speed function of Ichoua et al. [26] to a continuous travel time function. An example of this approach is provided in Figure 3. In this method, a vector U is defined as initial distinct intervals for the mth arc between two nodes i and j. Therefore, at first, T can be represented by T=T1,T2,…,TUmij. Aftermath, new time intervals must be calculated for all links with various traffic congestions. According to Setak et al. [27], equation (1) is employed to achieve new time intervals. In this equation, vijmiju represents the travel speed at the mth link i,j in uth time interval. In the following, the number of time intervals in travel time function (H) is calculated by H=2U−1+1 where U shows the number of time intervals in the function of time speed. Finally, two coefficients a and b are calculated according to equations (5) and (6).(1)Tnewijm=Tnewijm1,Tnewijm2,Tnewijm3,…,TnewijmH−1,TnewijmH=T1,T2−T1−dijmvijm1,…,TUm−TUm−1−,TUm.Figure 3
Travel time function corresponds to travel speed function in Figure2.IfTnewijmh<Tnewijmh−1, then Tnewijmh=Tnewijmh−1.In addition,(2)xijm=dijmvijm1,dijmvijm2,…,dijmvijmU,(3)yijm=xijm1,xijm1,xijm2,xijm2…,xijmU,(4)bijm=yijm2−yijm1Tnewijm2−Tnewijm1,yijm3−yijm2Tnewijm3−Tnewijm2,…,yijmHm−yijmHm−1TnewijmHm−TnewijmHm–1,(5)aijm=yijm1−bijm×Tnewijm1,…,yijmhm−bijm×TnewijmHm.Henceforth, we can calculate the travel time at themth edge among nodes i and j by the following formula:(6)t=aijmh+bijmh×t0,where t0 is the vehicle departure time from node i to node j throughout time interval [Tnewijmh, Tnewijmh+1].


### 3.2. Reliability of Distribution

In the proposed problem, a discrete function with a finite set of independent scenarios is considered to formulate the reliability of each route. In this regard, the reliability ofmth arc between vertexes i and j in scenario ω is considered as successful transportation through mth arc between two successive nodes i and j presented by rijmω. In other words, the reliability is addressed as the probability for vehicles to traverse the connections between two demand nodes in time aftermath of a disaster. Evidently, after the occurrence of disasters, evaluating the status of links (including reliability and traffic congestions) is more complicated due to the dynamic situation. However, according to literatures [2, 4], the values of reliability and related traffic congestions on links can be estimated based on subjective perceptions instead of observed data. Thereupon, the reliability of kth route performed by vehicle k in scenario ω can be computed by knowing the reliability of each employed link as follows:(7)reliabilityroutek=∏i,j,m∈E|∑h∈Hmxijmhkωrijmω.


### 3.3. Semisoft Time Window

Using a semisoft time window, a demand node is permitted to be served after the defined time window (TW) by considering a penalty for the delayed arrival time [28]. These penalties can also be assumed to vary from one affected area to another based on their priorities. Figure 4 shows the nondecreasing time-dependent arrival penalty in SR-TDRRPM-SSTW.Figure 4
Cost function for delayed services.Consider that the service completion time at the demandi is tikω. Moreover, C represents the unit penalty cost for the late arrival. The cost inducing from delayed service at node i in scenario ω is formulated as a function of tikω using the following:(8)penalty cost=0,iftikω≤TWiα,C.tikω−TWiα,ifTWiα≤tikω≤TWiβ,∞,ifTWiβ≤tikω.


## 3.1. FIFO (First-In, First-Out) Property Satisfaction

The nonpassing or the first-in, first-out (FIFO) property guarantees that between two identical vehiclesA and B, if vehicle A begins to traverse the arc i,j earlier, then it arrives at its destination j before vehicle B. Primitive works on time-dependent vehicle routing problem did not satisfy FIFO property [26]. In fact, they utilized a discrete function of time to model the travel times. An example of their approach is illustrated in Figure 2 [26]. As it can be seen from the figure if a vehicle departs the origin at 4 : 30, it will arrive at its destination in 150 miles away at 9 : 30. Nonetheless, if a vehicle departs at 6 : 00, it will arrive to its destination at 9 : 00. It is clear that this situation is not feasible in the real word and the result does not conform the FIFO property.Figure 2
An example of ‘‘passing” with discrete travel time function.Recently, scholars employ continuous travel time functions over the time horizon. In this paper, we employed the approach presented by Setak et al. [27] to convert a travel speed function of Ichoua et al. [26] to a continuous travel time function. An example of this approach is provided in Figure 3. In this method, a vector U is defined as initial distinct intervals for the mth arc between two nodes i and j. Therefore, at first, T can be represented by T=T1,T2,…,TUmij. Aftermath, new time intervals must be calculated for all links with various traffic congestions. According to Setak et al. [27], equation (1) is employed to achieve new time intervals. In this equation, vijmiju represents the travel speed at the mth link i,j in uth time interval. In the following, the number of time intervals in travel time function (H) is calculated by H=2U−1+1 where U shows the number of time intervals in the function of time speed. Finally, two coefficients a and b are calculated according to equations (5) and (6).(1)Tnewijm=Tnewijm1,Tnewijm2,Tnewijm3,…,TnewijmH−1,TnewijmH=T1,T2−T1−dijmvijm1,…,TUm−TUm−1−,TUm.Figure 3
Travel time function corresponds to travel speed function in Figure2.IfTnewijmh<Tnewijmh−1, then Tnewijmh=Tnewijmh−1.In addition,(2)xijm=dijmvijm1,dijmvijm2,…,dijmvijmU,(3)yijm=xijm1,xijm1,xijm2,xijm2…,xijmU,(4)bijm=yijm2−yijm1Tnewijm2−Tnewijm1,yijm3−yijm2Tnewijm3−Tnewijm2,…,yijmHm−yijmHm−1TnewijmHm−TnewijmHm–1,(5)aijm=yijm1−bijm×Tnewijm1,…,yijmhm−bijm×TnewijmHm.Henceforth, we can calculate the travel time at themth edge among nodes i and j by the following formula:(6)t=aijmh+bijmh×t0,where t0 is the vehicle departure time from node i to node j throughout time interval [Tnewijmh, Tnewijmh+1].


## 3.2. Reliability of Distribution

In the proposed problem, a discrete function with a finite set of independent scenarios is considered to formulate the reliability of each route. In this regard, the reliability ofmth arc between vertexes i and j in scenario ω is considered as successful transportation through mth arc between two successive nodes i and j presented by rijmω. In other words, the reliability is addressed as the probability for vehicles to traverse the connections between two demand nodes in time aftermath of a disaster. Evidently, after the occurrence of disasters, evaluating the status of links (including reliability and traffic congestions) is more complicated due to the dynamic situation. However, according to literatures [2, 4], the values of reliability and related traffic congestions on links can be estimated based on subjective perceptions instead of observed data. Thereupon, the reliability of kth route performed by vehicle k in scenario ω can be computed by knowing the reliability of each employed link as follows:(7)reliabilityroutek=∏i,j,m∈E|∑h∈Hmxijmhkωrijmω.


## 3.3. Semisoft Time Window

Using a semisoft time window, a demand node is permitted to be served after the defined time window (TW) by considering a penalty for the delayed arrival time [28]. These penalties can also be assumed to vary from one affected area to another based on their priorities. Figure 4 shows the nondecreasing time-dependent arrival penalty in SR-TDRRPM-SSTW.Figure 4
Cost function for delayed services.Consider that the service completion time at the demandi is tikω. Moreover, C represents the unit penalty cost for the late arrival. The cost inducing from delayed service at node i in scenario ω is formulated as a function of tikω using the following:(8)penalty cost=0,iftikω≤TWiα,C.tikω−TWiα,ifTWiα≤tikω≤TWiβ,∞,ifTWiβ≤tikω.


## 4. Mathematical Formulation

The purpose of the proposed SR-TDRRPM-SSTW is to maximize the reliability of distribution together with minimizing the sum of service completion times in delivering prioritized items considering the late arrival penalties in the post-disaster situation under uncertainty. In the proposed problem, the following assumptions are considered:(i)
Vehicles are assumed to be homogenous, depart from the distribution center in the first horizon, and finally go back to it once the distribution is finished(ii)
The locations of affected nodes and candidate central depot are known prior(iii)
The probability of successful transportations between vertexes is also estimated for each of the scenarios in advance(iv)
The problem is studied in multigraph setting where the alternative links are differentiated by two attributes(v)
The reliability of links, distance, and time-dependent traffic congestions can be estimated by transportation research community’s expertise or using advanced disaster detection technology in real-time(vi)
The quantities of delivery items demanded by shelters for each scenario are known(vii)
The SR-TDRRPM-SSTW only considers the affected areas that are accessible by vehicles in inner-city areasTo deal with the uncertainty of the parameters in SR-TDRRPM-SSTW, it is modeled by two-stage stochastic formulation. To do this, a scenario tree including a set of independent scenarios is considered to estimate the stochastic demands, reliability of links, service duration time, and travel speeds. The decision variables can be segmented into two groups consisting scenario-independent variablesxijk which are known as first-stage decisions and scenario-dependent variables xijmhkω,tikω which are known as second-stage decisions. The second-stage decisions are influenced by the uncertain parameters and made after realization of the first-stage actions [29]. The first-stage problem is a capacitated vehicle routing problem where the objectives attempt to minimize the expected operational costs and late arrival penalties by equation (9) as well as maximize the expected reliability of constructed links related to the second-stage subproblems P2 with equation (10).Using the notations presented in Notation and Definitions, the mathematic formulation of stage one (P1) and stage two (P2) of SR-TDRRPM-SSTW can be modeled as follows.


### 4.1. First Stage (P1)

(9)Objective1:Eξhx,ξ,(10)Objective2:Eξgx,ξ.Constrains are as follows:(11)∑i∈0∪​N,i≠j∑k∈Kxijk=1,∀j∈N,(12)∑j∈n+1∪​N,i≠j∑k∈Kxijk=1,∀i∈N,(13)∑i∈0∪​Nxijk=∑i∈n+1∪​Nxjik,∀j∈N,∀k∈K,(14)∑j∈0∪​N∑k∈Kx0jk≤K,(15)xiik=0,∀i∈N∪​P,∀k∈K,(16)xi0k=0,∀i∈N∪​P,∀k∈K,(17)xn+1jk=0,∀j∈N,∀k∈K,(18)xijk=0,1,∀i,j∈N∪​P,∀k∈K.Constraints (11)–(13) define the flow on the path of vehicle k. Constraint (14) limits the total number of available vehicles. Constraints (15)–(18) guarantee the accuracy of achieved solutions. Next, for each scenario ω∈Ω, P2 can be written as follows.


### 4.2. Second Stage (P2)

(19)hx,ξω=MinC1∑k∈Ktn+1kω+C2∑i∈0∪​N∑k∈Kθikω,(20)gx,ξω=Max∏i,j,m∈E|∑h∈Hmxijmhkωrijmω,∀k∈K,(21)∑i∈0∪​N,i≠j∑m∈Mij∑h∈Hmxijmhkω=xijk,∀i,j∈N,∀k∈K,(22)tjkω−tikω≥aijmhω+bijmhωtikω+sjω+xijmhkω−1E,(23)tikω≤E.∑j∈N∪​n+1xijk,∀i∈N,∀k∈K,(24)tikω−siω≤TWiβ,∀i∈N,∀k∈K,(25)tikω−siω−TWiα≤θikω,∀i∈N,∀k∈K,(26)∑i∈NDiω∑j∈N∪​n+1∑m∈Mij∑h∈Hmxijmhkω≤Q,∀k∈K,(27)xijmhkω=0,1,∀i,j∈N∪​P,∀m∈Mij,∀h∈Hm,∀k∈K,(28)tikω≥0,∀i∈N∪​P,∀k∈K,(29)θikω≥0,∀i∈N∪​P,∀k∈K.Constraint (21) is employed to restrict the vehicle k to select only one path between two nodes i and j. Constraints (22)–(23) compute the service completion time at each demand node in each scenario. Constraint (24) imposes that the maximum service completion time should be less than TWiβ. Moreover, constraint (25) computes the allowable late arrival time at each node. Constraint (26) controls the capacity of vehicles along each constructed route. Constraints (27)–(29) define the type of variables.


## 4.1. First Stage (P1)

(9)Objective1:Eξhx,ξ,(10)Objective2:Eξgx,ξ.Constrains are as follows:(11)∑i∈0∪​N,i≠j∑k∈Kxijk=1,∀j∈N,(12)∑j∈n+1∪​N,i≠j∑k∈Kxijk=1,∀i∈N,(13)∑i∈0∪​Nxijk=∑i∈n+1∪​Nxjik,∀j∈N,∀k∈K,(14)∑j∈0∪​N∑k∈Kx0jk≤K,(15)xiik=0,∀i∈N∪​P,∀k∈K,(16)xi0k=0,∀i∈N∪​P,∀k∈K,(17)xn+1jk=0,∀j∈N,∀k∈K,(18)xijk=0,1,∀i,j∈N∪​P,∀k∈K.Constraints (11)–(13) define the flow on the path of vehicle k. Constraint (14) limits the total number of available vehicles. Constraints (15)–(18) guarantee the accuracy of achieved solutions. Next, for each scenario ω∈Ω, P2 can be written as follows.


## 4.2. Second Stage (P2)

(19)hx,ξω=MinC1∑k∈Ktn+1kω+C2∑i∈0∪​N∑k∈Kθikω,(20)gx,ξω=Max∏i,j,m∈E|∑h∈Hmxijmhkωrijmω,∀k∈K,(21)∑i∈0∪​N,i≠j∑m∈Mij∑h∈Hmxijmhkω=xijk,∀i,j∈N,∀k∈K,(22)tjkω−tikω≥aijmhω+bijmhωtikω+sjω+xijmhkω−1E,(23)tikω≤E.∑j∈N∪​n+1xijk,∀i∈N,∀k∈K,(24)tikω−siω≤TWiβ,∀i∈N,∀k∈K,(25)tikω−siω−TWiα≤θikω,∀i∈N,∀k∈K,(26)∑i∈NDiω∑j∈N∪​n+1∑m∈Mij∑h∈Hmxijmhkω≤Q,∀k∈K,(27)xijmhkω=0,1,∀i,j∈N∪​P,∀m∈Mij,∀h∈Hm,∀k∈K,(28)tikω≥0,∀i∈N∪​P,∀k∈K,(29)θikω≥0,∀i∈N∪​P,∀k∈K.Constraint (21) is employed to restrict the vehicle k to select only one path between two nodes i and j. Constraints (22)–(23) compute the service completion time at each demand node in each scenario. Constraint (24) imposes that the maximum service completion time should be less than TWiβ. Moreover, constraint (25) computes the allowable late arrival time at each node. Constraint (26) controls the capacity of vehicles along each constructed route. Constraints (27)–(29) define the type of variables.


## 5. Solution Methodology: Modified NSGA-II

The NSGA-II is a popular multiobjective evolutionary algorithm which was first offered by Deb et al. [30]. It provides an approximation of the nondominated Pareto front by using the successive computation of a series of generations of solutions. In each iteration, this algorithm generates new populations by employing crossover and mutation operators. Afterward, members of the current and new generations are mixed together. Finally, population members are ranked according to the crowding distance and nondominance concepts to select best solutions. Solutions which lie above the Pareto front are nonoptimal solutions; however, solutions that lie under the Pareto front are infeasible. All solutions on the Pareto front are optimal by considering a specified objective weight preference. Herein, first, the solution representation is presented and then the procedures of the operators are briefly described.


### 5.1. Solution Representation

Solution scheme in the modified NSGA-II method is designated to specify the assignment of demand points to various routes and determine the order of vertexes for each vehicle. Herein, we applied a continuous representation form (CRF) as aS×N+K−1 dimension matrix which includes numbers in the interval [0, 1]. In order to obtain the sequence of vertexes, the numbers in the first row of CRF are sorted decently; then, the nodes’ sequences can be achieved by comparing the location of each number in the initial CRF and in the sorted version. Next, the m link between two sequential nodes is selected if the value of gene be in the interval of m−1/Mij≤value<m/Mij. Figure 5 illustrates a sample of the decoding method for SR-TDRRPM-SSTW with 6 vertexes and 2 vehicles in two scenarios with multiple arcs between points.Figure 5
Solution representation and decoding process.


### 5.2. Crossover Operator

Crossover is one of the fundamental operators in the genetic algorithm which is also referred to as recombination. It is used to explore new solution space. This operator employs two distinct chromosomes as parents to generate new offspring. The obtained offspring inherit the characteristics of their parents with a certain method. There exist various methods to combine two given parents. In this study, we applied two-point crossover to generate new individuals. In this approach, we select two points on the parent chromosome strings. Then, genes between these two points are swapped between the parent structures. An example process of crossover operator is provided in Figure6.Figure 6
Example of crossover operator.


### 5.3. Mutation Operators

Mutation operator is mainly responsible for maintaining diversity in the population. Mutation is conducted via altering a randomly selected gene. Here, we used four mutation operators which are defined as follows:Two-point swap: this operator selects two alleles of the chromosome randomly. Then, for creating one new chromosome, two selected alleles will be swapped.Insertion: it moves a randomly selected gene (allele) to another randomly selected place.Reversion: this operator reverses two randomly selected-point cuts.Arc exchange: in this operator, the links without changing the sequence of nodes will be randomly exchanged.


### 5.4. Migration Operator

We also applied another operator called immigration. As it occurs in many societies, some foreign individuals with various genetic characteristics are added to the current population periodically. Immigration operator is performed to escape local optima by exploring new solution space. To do this, a part of poor individuals at each generation is replaced with some new randomly created immigrants [31].


### 5.5. Caching Mechanism

Kratica et al. [32] employed the caching method in the genetic algorithm (GA) for the first time. This technique improves the algorithm’s running time through storing a certain number of already measured fitness functions in a cache table. In this study, for caching, a practical strategy is used called Least Recently Used (LRU) employing a hash-queue data structure [33]. Specifically, the fitness function of all chromosomes is gathered in a cache table, then, in case of subsequent iteration of a solution, the corresponding objective value is rapidly retrieved.


### 5.6. Parameters Tuning

The Taguchi method is an efficient method for determining parameters of a metaheuristic algorithm. This method uses the signal to noise (S/N) ratio to find the best level of each parameter. Table1 reports the input parameters of the proposed algorithm and their levels. Figure 7 demonstrates the obtained S/N ratio corresponding to each defined parameter. In order to implement the Taguchi method, the Spacing Metric (SM) is considered as a measure for evaluating the optimal results [34]. The SM factor is achieved by equation SM=∑i=1n−1d¯−di/n−1d¯ where n is total number of solutions in the final Pareto front, di represents the Euclidian distance between every successive Pareto solutions, and the average of all di is denoted by d¯. Higher values in Figure 7 represent the more appropriate levels for intended parameters.Table 1
Parameters and their level in GA.
ParameterA: crossover percentageB: mutation percentageC: immigration percentageD: NPOPE: MaxItLevel #10.600.050.056075Level #20.650.100.108090Level #30.700.150.15100100NPOP: number of population; MaxIt: maximum number of iterations.Figure 7
S/N ratio plots for the modified NSGA-II.


## 5.1. Solution Representation

Solution scheme in the modified NSGA-II method is designated to specify the assignment of demand points to various routes and determine the order of vertexes for each vehicle. Herein, we applied a continuous representation form (CRF) as aS×N+K−1 dimension matrix which includes numbers in the interval [0, 1]. In order to obtain the sequence of vertexes, the numbers in the first row of CRF are sorted decently; then, the nodes’ sequences can be achieved by comparing the location of each number in the initial CRF and in the sorted version. Next, the m link between two sequential nodes is selected if the value of gene be in the interval of m−1/Mij≤value<m/Mij. Figure 5 illustrates a sample of the decoding method for SR-TDRRPM-SSTW with 6 vertexes and 2 vehicles in two scenarios with multiple arcs between points.Figure 5
Solution representation and decoding process.


## 5.2. Crossover Operator

Crossover is one of the fundamental operators in the genetic algorithm which is also referred to as recombination. It is used to explore new solution space. This operator employs two distinct chromosomes as parents to generate new offspring. The obtained offspring inherit the characteristics of their parents with a certain method. There exist various methods to combine two given parents. In this study, we applied two-point crossover to generate new individuals. In this approach, we select two points on the parent chromosome strings. Then, genes between these two points are swapped between the parent structures. An example process of crossover operator is provided in Figure6.Figure 6
Example of crossover operator.


## 5.3. Mutation Operators

Mutation operator is mainly responsible for maintaining diversity in the population. Mutation is conducted via altering a randomly selected gene. Here, we used four mutation operators which are defined as follows:Two-point swap: this operator selects two alleles of the chromosome randomly. Then, for creating one new chromosome, two selected alleles will be swapped.Insertion: it moves a randomly selected gene (allele) to another randomly selected place.Reversion: this operator reverses two randomly selected-point cuts.Arc exchange: in this operator, the links without changing the sequence of nodes will be randomly exchanged.


## 5.4. Migration Operator

We also applied another operator called immigration. As it occurs in many societies, some foreign individuals with various genetic characteristics are added to the current population periodically. Immigration operator is performed to escape local optima by exploring new solution space. To do this, a part of poor individuals at each generation is replaced with some new randomly created immigrants [31].


## 5.5. Caching Mechanism

Kratica et al. [32] employed the caching method in the genetic algorithm (GA) for the first time. This technique improves the algorithm’s running time through storing a certain number of already measured fitness functions in a cache table. In this study, for caching, a practical strategy is used called Least Recently Used (LRU) employing a hash-queue data structure [33]. Specifically, the fitness function of all chromosomes is gathered in a cache table, then, in case of subsequent iteration of a solution, the corresponding objective value is rapidly retrieved.


## 5.6. Parameters Tuning

The Taguchi method is an efficient method for determining parameters of a metaheuristic algorithm. This method uses the signal to noise (S/N) ratio to find the best level of each parameter. Table1 reports the input parameters of the proposed algorithm and their levels. Figure 7 demonstrates the obtained S/N ratio corresponding to each defined parameter. In order to implement the Taguchi method, the Spacing Metric (SM) is considered as a measure for evaluating the optimal results [34]. The SM factor is achieved by equation SM=∑i=1n−1d¯−di/n−1d¯ where n is total number of solutions in the final Pareto front, di represents the Euclidian distance between every successive Pareto solutions, and the average of all di is denoted by d¯. Higher values in Figure 7 represent the more appropriate levels for intended parameters.Table 1
Parameters and their level in GA.
ParameterA: crossover percentageB: mutation percentageC: immigration percentageD: NPOPE: MaxItLevel #10.600.050.056075Level #20.650.100.108090Level #30.700.150.15100100NPOP: number of population; MaxIt: maximum number of iterations.Figure 7
S/N ratio plots for the modified NSGA-II.


## 6. Computational Results

### 6.1. An Illustrative Example

In this section, we provide a small-sized sample in a deterministic situation to illustrate that how the model proposes a distribution plan considering different traffic conditions and reliabilities of links. This instance contains one distribution center and seven beneficiary regions. Table2 gives the gives the demands of each region. Moreover, the related weights are calculated by wi=Dij/Dj×10. For all of nodes, the time window and service time are set to 15 and 1, respectively. One vehicle with capacity 400 is utilized to transship emergency supplies from the depot to various destinations. As depicted in Figure 8(a), the locations of DC and demand nodes are dispersed in a two-dimensional space and the distances are corresponding to the Euclidean distance between locations. The reliability of the connections between different points is categorized to three cases. Furthermore, within the planning horizon, four numbers of time intervals are taken into consideration for every passable links. The time intervals are [0, 3], [3, 6], [6, 9], and more than 9. The upper limit of the last period (the fourth time interval) is assumed to be a large number to guaranty that the vehicle’s travel time is finished in the last interval. We assume three different scenarios for traffic congestions. In the first scenario, the traffic is high, and most of the times, the speeds are low. While in the second and third scenarios, we encounter medium- and low-traffic congestions with higher traveling speeds. The speed patterns under various scenarios are shown in Figure 9 by different colors where dotted arcs depict less reliable connections. Figure 8(b) shows one of the routing schemes with reliability=0.9025.Table 2
The data of the small-scale instance.
Number (i)1234567Di50806040704555wi1.2521.511.751.1251.375Figure 8
(a) Sample multigraph network; (b) optimum solution.Figure 9
(a) Speed changes for the first scenario-high traffic; (b) speed changes for the second scenario-medium traffic; (c) speed changes for the third scenario-low traffic.
(a)(b)(c)


### 6.2. Computational Experiments

Here, a numerical analysis is performed to investigate the performance of the modified metaheuristic. In this regard, we utilized the test instances prepared by Tikani and Setak [2] with some adaptations. Same as Tikani and Setak [2], two alternative links between all vertexes are considered in a multigraph setting. In this study, we consider two scenarios with the respective probabilities of 60% and 40%. The reliability of links, service time duration, and demands of the affected areas in the first scenario are in accordance with Tikani and Setak [2], and the parameters in the second scenario are also generated by the same distribution. We assume that the second scenario has a lower traffic condition in comparison to the first one. To this end, the speeds are generated in the interval [20, 35]. Moreover, the threshold in the semisoft TW is defined as same as the upper bound of Tikani and Setak [2] and TWiα=0.8×TWiβ. Moreover, C1=10andC2=2.To measure the effectiveness of the proposed NSGA-II, the performance of this algorithm is evaluated with standard NSGA-II through five metrics. Table3 summarizes the performance measurements and their definitions. As mentioned in Table 3, the higher values of QM indicate the supremacy of the proposed algorithm as well as the lower values of SM, MID, and computation time [35]. Table 4 shows the comparison results of the described metrics. Moreover, Figure 10 visualizes the box plots of the performance measures. According to the obtained results, modified NSGA-II shows the better performance compared with NSGA-II in all defined metrics.Table 3
Definition of performance metrics.
MeasureDefinitionCriterion∗FormulaNOSIndicates the number of achieved solution in the Pareto frontL—Spacing metric (SM)Shows that whether the achieved solutions in the Pareto front are distributed evenlySSM=∑i=1n−1d¯−di/n−1d¯Mean ideal distance (MID)Presents the mean distance among the solutions of Pareto set from the ideal solutionSMID=∑i=1nF1i−F1best/F1Totalmax−F1Totalmin2+F2i−F2best/F2Totalmax−F2Totalmin2/mQuality metric (QM)Presents the contribution of the algorithm regarding to all nondominated solutionsLQM=number of non_dominated solutions of the method/total number of non_dominated solutionsSolving time (CPU time)Running time of the algorithmS—S: smaller value is better; L: larger value is better.Table 4
Comparison metrics for different size instances.
Problem no.Modified NSGA-IIStandard NSGA-IINOSSMMIDQMCPU time (s)NOSSMMIDQMCPU time (s)140.1280.506154.2020.1020.5450.5067.23250.3600.4490.7156.3240.3310.6510.5069.17350.3900.445160.4130.9580.4660.6073.52440.0960.368164.1730.1720.4690.7577.16560.2940.4350.8566.5830.5180.6870.1179.54650.0590.3200.6272.1850.2370.6330.3780.69760.0480.290175.3940.1980.5910.6686.02860.2700.125176.0930.6100.3530.5089.16970.1360.2260.8777.3940.7480.5090.2590.031060.0680.4320.7580.1150.4050.6310.5087.511180.4730.2650.8084.0650.4620.4460.6097.621280.4660.326185.8460.7110.5080.63100.41390.1140.1950.8187.1150.4590.3590.45103.91480.3100.3550.8889.7750.5360.2050.11106.715100.1510.2090.7692.0170.5040.2440.38111.61690.2560.2490.8196.2160.4250.3020.54113.417110.3370.2740.9196.9980.2970.3970.41115.018100.1640.3210.8398.1970.4140.3130.58117.819110.2200.2980.78100.180.6020.4260.50119.320120.1530.3601104.990.4420.3490.58122.721110.3270.3010.84106.3110.5260.3970.53124.622130.5560.3250.76108.1100.4420.4010.47125.023160.2880.1940.94110.0110.8990.2730.52127.624150.3800.2200.83114.3120.5640.2580.61129.025160.2880.3621117.9130.4420.3590.69130.8Average8.840.2530.3140.8786.986.3600.4800.4300.49101.81Figure 10
Comparison metrics for the standard NSGA-II (red boxes) and the modified NSGA-II (green boxes).
(a)(b)(c)(d)


### 6.3. Statistical Analysis

In order to present a more precise comparison, we performed a statistical analysis using the pairedt-test for each of the metrics. The null hypothesis of the test assumes that the difference between the proposed NSGA-II and standard NSGA-II for each specific metric is insignificant, while the alternative hypothesis states that the proposed NSGA-II yields better overall outcomes than the standard NSGA-II regarding the studied performance metrics.Herein, at first, the Anderson–Darling test is utilized to check the normality of all metrics, which shows the normal distribution of applied metrics’ values. The results of the applied pairedt-test are reported in Table 5, and it shows that in all tests, null hypothesis is not accepted at confidence level 0.95. The statistical analysis confirms the superiority of the modified NSGA-II compared with its standard version.Table 5
TheP values of the paired t-test for different metrics.
Metric’s nameT valueP valueTest resultsNOS10.70×SM−5.820×MID−5.150×QM12.080×CPU time−22.340××: null hypothesis is not accepted.


### 6.4. Significance of considering Nondominated Alternative Links

In order to evaluate the importance of using multigraph setting in the SR-TDRRPM-SSTW, in this section, we discard the availability of multiple paths between nodes. To this end, we randomly remove one of the parallel links among the vertexes. Then, the problems are resolved in the achieved single graph. The best obtained value regarding each objective function is described in Table6. As we expected, considering nondominated extra edges between nodes can substantially improve the quality of solutions. The mean of average achievements for four instances is relatively 10.4% in which improvements are calculated by the following:(30)improvement%=100×OBJM−OBJSOBJS,where OBJM is the objective of multigraph and OBJS is the objective of simple graph network.Table 6
The best achieved values of each objective in multigraph network versus simple graph.
Problem no.MultigraphSimple graphAverage of improvementsOF 1OF 2OF 1OF 24167.10.92184.20.8012.1410190.30.84206.80.769.2516218.00.75234.60.688.6822287.90.74324.30.6611.66


### 6.5. A Real Transportation Network

Iran experiences a variety of different natural disasters on a frequent basis. This country is located in a very active seismic zone, called Alpide belt. During the past decades, several major and numerous minor earthquakes have occurred in this country. Generally, in Iran, earthquakes with more than two Richter can result in huge social and economic losses as well as casualties [12]. In this research, to demonstrate the feasibility of the logistics model, we studied a possible earthquake in the real-world network of Isfahan. The province of Isfahan is one of the most historical and industrial regions of Iran. Location of Isfahan on the map and the major active faults of Iran are shown in Figure 11(a). According to geographical research in a 100-kilometer radius of Isfahan city, 63 faults were found nearby this city; some of these faults were hundreds of kilometers in length [36]. The fault density map and earthquakes detected in the mentioned region in 2000–2003 are displayed in Figure 11(b).Figure 11
(a) Major active faults and location of Isfahan in the map of Iran; (b) satellite image for fault density in 100 km radius of Isfahan and detected earthquakes in 2000–2003.
(a)(b)In this study, as illustrated in Figure12, the underlying transportation network comprises 6 shelters and one distribution center. The distribution center is situated in Bakhtiar-Dasht, and also the shelters are in Isfahan’s disparate areas. To deliver emergency packages from the depot to the different shelters, three similar vehicles with a capacity of 2000 exist. Required information regarding the problem is provided in Table 7. The shelters’ demands are assessed according to the polls provided by area’s crisis management experts. Times of service and priorities are set based on the nodes’ demands in accordance with two formulations Si1=1.5h+0.05×∑ωDiω/2 and Si2=2h+0.05×∑ωDiω/2 and wi=10×∑ωDiω/2×∑ω∑jDjω, respectively. In order to study the case study in a larger scale, we added fourteen potential shelter locations by consultation with experts (including some existing stadiums, schools, and college locations in Isfahan), and their related demands are estimated on the interval of [500, 900]. Other required parameters are also computed based on the abovementioned formulations.Figure 12
Traffic zones with various traffic patterns and representation of the logistics network (Google Maps).Table 7
Information regarding the main shelters in the network.
SheltersLocationsDemandsPrioritiesService times (hour)Scenario 1Scenario 2Scenario 1Scenario 21Bagh-Fadak105012502.03542.272.572Pirouzi Stadium70011001.59292.152.453Enghelab Stadium6608401.32742.072.374Ghadir Park100011001.85842.222.525Vaziche Stadium7508501.41592.102.406Abrisham85011501.76992.202.50The reliability of links between nodes and the speed pattern in the disaster aftermath are set in accord with deliberation with specialists in crisis management and civil architecture fields. We classified the region based on three traffic zones of heavy, medium, and light to determine the traffic pattern, which are shown in Figure12. This classification is performed based on existing transportation statistics of Isfahan. For each category of traffic, Figure 13 depicts the speed patterns for four time intervals [0, 3], [3, 6], [6, 9], and more than 9. It is notable that in some cases a route from one node to another one may include various traffic regions. In this situation, it is required that some breakpoints be specified in the indicated route where the vehicles enter to another traffic area. Hence, by considering these breakpoints, travel times are calculated using the method of Section 3.1.Figure 13
(a) Speed changes of the light traffic zone; (b) speed changes of the medium traffic zone; (c) speed changes of the heavy traffic zone (continues line: first scenario; dash line: second scenario).
(a)(b)(c)Owing to the urban networks’ structure, it is evident that the distances among nodes are not symmetric. In this regard, two alternative links employed between each pair of nodes are obtained from Google Maps. We report the details of links’ distances and their reliabilities between the main shelter locations in Table8. For the rest of the demand points, detailed information is not presented for the sake of brevity.Table 8
Details of distances (km) for two parallel links between main shelter locations in multigraph.
Depoti1i2i3i4i5i6Depot×20.1-122.8-123.7-125.7-116.2-129.6-0.9524.9-125.5-123.9-131.4-116.4-133.1-0.95i123.3-1×11.4-120.7-0.9514.3-120.1-130.1-0.9526.9-115.9-121.1-0.9515.9-121.2-137.3-0.95i226.0-112.4-1×7.90-0.953.2-119.2-129.2-0.9531.4-115.2-113.9-0.953.9-124.0-132.3-0.90i324.8-0.9519.5-0.959.0-0.95×6.7-0.9514.0-0.9515.2-137.9-0.9522.2-0.9514.9-0.957.5-0.9020.7-0.8525.8-1i425.7-116.3-12.6-17.8-0.95×25.1-124.4-0.9531.9-120.3-13.0-113.7-0.9529.9-0.932.1-0.95i517.0-118.6-120.0-113.5-0.9025.8-1×17.2-0.9017.6-119.9-120.5-114.6-0.9526.3-123.6-0.95i617.3-0.9030.1-0.9525.9-0.9516.3-123.6-0.9518.8-0.90×35.5-0.9546.9-0.9530.7-0.9524.6-131.3-0.9524.8-0.95We solved the proposed model once for a logistics network which includes only the main shelter locations (Case I) and once for a logistics network considering both main and potential shelter locations (Case II). By solving the proposed model, a series of solutions as a Pareto front are achieved. In Figures14 and 15, we depicted the obtained results for Case I and Case II, respectively. As can be seen from the figures, increasing the reliability of routes yields to increase the first objective function. It accentuates the conflict of objectives. Herein, a decision maker should select one of the solutions in the Pareto front according to his preference. By comparing Figures 14 and 15, it becomes clear that by increasing the number of demand points in Case II, the overall reliability (objective 2) of routes decreases and cost of traveling (objective 1) increases since a vehicle should serve more demand nodes.Figure 14
Approximate Pareto optimal for the studied Case I.Figure 15
Approximate Pareto optimal for the studied Case II.


## 6.1. An Illustrative Example

In this section, we provide a small-sized sample in a deterministic situation to illustrate that how the model proposes a distribution plan considering different traffic conditions and reliabilities of links. This instance contains one distribution center and seven beneficiary regions. Table2 gives the gives the demands of each region. Moreover, the related weights are calculated by wi=Dij/Dj×10. For all of nodes, the time window and service time are set to 15 and 1, respectively. One vehicle with capacity 400 is utilized to transship emergency supplies from the depot to various destinations. As depicted in Figure 8(a), the locations of DC and demand nodes are dispersed in a two-dimensional space and the distances are corresponding to the Euclidean distance between locations. The reliability of the connections between different points is categorized to three cases. Furthermore, within the planning horizon, four numbers of time intervals are taken into consideration for every passable links. The time intervals are [0, 3], [3, 6], [6, 9], and more than 9. The upper limit of the last period (the fourth time interval) is assumed to be a large number to guaranty that the vehicle’s travel time is finished in the last interval. We assume three different scenarios for traffic congestions. In the first scenario, the traffic is high, and most of the times, the speeds are low. While in the second and third scenarios, we encounter medium- and low-traffic congestions with higher traveling speeds. The speed patterns under various scenarios are shown in Figure 9 by different colors where dotted arcs depict less reliable connections. Figure 8(b) shows one of the routing schemes with reliability=0.9025.Table 2
The data of the small-scale instance.
Number (i)1234567Di50806040704555wi1.2521.511.751.1251.375Figure 8
(a) Sample multigraph network; (b) optimum solution.Figure 9
(a) Speed changes for the first scenario-high traffic; (b) speed changes for the second scenario-medium traffic; (c) speed changes for the third scenario-low traffic.
(a)(b)(c)


## 6.2. Computational Experiments

Here, a numerical analysis is performed to investigate the performance of the modified metaheuristic. In this regard, we utilized the test instances prepared by Tikani and Setak [2] with some adaptations. Same as Tikani and Setak [2], two alternative links between all vertexes are considered in a multigraph setting. In this study, we consider two scenarios with the respective probabilities of 60% and 40%. The reliability of links, service time duration, and demands of the affected areas in the first scenario are in accordance with Tikani and Setak [2], and the parameters in the second scenario are also generated by the same distribution. We assume that the second scenario has a lower traffic condition in comparison to the first one. To this end, the speeds are generated in the interval [20, 35]. Moreover, the threshold in the semisoft TW is defined as same as the upper bound of Tikani and Setak [2] and TWiα=0.8×TWiβ. Moreover, C1=10andC2=2.To measure the effectiveness of the proposed NSGA-II, the performance of this algorithm is evaluated with standard NSGA-II through five metrics. Table3 summarizes the performance measurements and their definitions. As mentioned in Table 3, the higher values of QM indicate the supremacy of the proposed algorithm as well as the lower values of SM, MID, and computation time [35]. Table 4 shows the comparison results of the described metrics. Moreover, Figure 10 visualizes the box plots of the performance measures. According to the obtained results, modified NSGA-II shows the better performance compared with NSGA-II in all defined metrics.Table 3
Definition of performance metrics.
MeasureDefinitionCriterion∗FormulaNOSIndicates the number of achieved solution in the Pareto frontL—Spacing metric (SM)Shows that whether the achieved solutions in the Pareto front are distributed evenlySSM=∑i=1n−1d¯−di/n−1d¯Mean ideal distance (MID)Presents the mean distance among the solutions of Pareto set from the ideal solutionSMID=∑i=1nF1i−F1best/F1Totalmax−F1Totalmin2+F2i−F2best/F2Totalmax−F2Totalmin2/mQuality metric (QM)Presents the contribution of the algorithm regarding to all nondominated solutionsLQM=number of non_dominated solutions of the method/total number of non_dominated solutionsSolving time (CPU time)Running time of the algorithmS—S: smaller value is better; L: larger value is better.Table 4
Comparison metrics for different size instances.
Problem no.Modified NSGA-IIStandard NSGA-IINOSSMMIDQMCPU time (s)NOSSMMIDQMCPU time (s)140.1280.506154.2020.1020.5450.5067.23250.3600.4490.7156.3240.3310.6510.5069.17350.3900.445160.4130.9580.4660.6073.52440.0960.368164.1730.1720.4690.7577.16560.2940.4350.8566.5830.5180.6870.1179.54650.0590.3200.6272.1850.2370.6330.3780.69760.0480.290175.3940.1980.5910.6686.02860.2700.125176.0930.6100.3530.5089.16970.1360.2260.8777.3940.7480.5090.2590.031060.0680.4320.7580.1150.4050.6310.5087.511180.4730.2650.8084.0650.4620.4460.6097.621280.4660.326185.8460.7110.5080.63100.41390.1140.1950.8187.1150.4590.3590.45103.91480.3100.3550.8889.7750.5360.2050.11106.715100.1510.2090.7692.0170.5040.2440.38111.61690.2560.2490.8196.2160.4250.3020.54113.417110.3370.2740.9196.9980.2970.3970.41115.018100.1640.3210.8398.1970.4140.3130.58117.819110.2200.2980.78100.180.6020.4260.50119.320120.1530.3601104.990.4420.3490.58122.721110.3270.3010.84106.3110.5260.3970.53124.622130.5560.3250.76108.1100.4420.4010.47125.023160.2880.1940.94110.0110.8990.2730.52127.624150.3800.2200.83114.3120.5640.2580.61129.025160.2880.3621117.9130.4420.3590.69130.8Average8.840.2530.3140.8786.986.3600.4800.4300.49101.81Figure 10
Comparison metrics for the standard NSGA-II (red boxes) and the modified NSGA-II (green boxes).
(a)(b)(c)(d)


## 6.3. Statistical Analysis

In order to present a more precise comparison, we performed a statistical analysis using the pairedt-test for each of the metrics. The null hypothesis of the test assumes that the difference between the proposed NSGA-II and standard NSGA-II for each specific metric is insignificant, while the alternative hypothesis states that the proposed NSGA-II yields better overall outcomes than the standard NSGA-II regarding the studied performance metrics.Herein, at first, the Anderson–Darling test is utilized to check the normality of all metrics, which shows the normal distribution of applied metrics’ values. The results of the applied pairedt-test are reported in Table 5, and it shows that in all tests, null hypothesis is not accepted at confidence level 0.95. The statistical analysis confirms the superiority of the modified NSGA-II compared with its standard version.Table 5
TheP values of the paired t-test for different metrics.
Metric’s nameT valueP valueTest resultsNOS10.70×SM−5.820×MID−5.150×QM12.080×CPU time−22.340××: null hypothesis is not accepted.


## 6.4. Significance of considering Nondominated Alternative Links

In order to evaluate the importance of using multigraph setting in the SR-TDRRPM-SSTW, in this section, we discard the availability of multiple paths between nodes. To this end, we randomly remove one of the parallel links among the vertexes. Then, the problems are resolved in the achieved single graph. The best obtained value regarding each objective function is described in Table6. As we expected, considering nondominated extra edges between nodes can substantially improve the quality of solutions. The mean of average achievements for four instances is relatively 10.4% in which improvements are calculated by the following:(30)improvement%=100×OBJM−OBJSOBJS,where OBJM is the objective of multigraph and OBJS is the objective of simple graph network.Table 6
The best achieved values of each objective in multigraph network versus simple graph.
Problem no.MultigraphSimple graphAverage of improvementsOF 1OF 2OF 1OF 24167.10.92184.20.8012.1410190.30.84206.80.769.2516218.00.75234.60.688.6822287.90.74324.30.6611.66


## 6.5. A Real Transportation Network

Iran experiences a variety of different natural disasters on a frequent basis. This country is located in a very active seismic zone, called Alpide belt. During the past decades, several major and numerous minor earthquakes have occurred in this country. Generally, in Iran, earthquakes with more than two Richter can result in huge social and economic losses as well as casualties [12]. In this research, to demonstrate the feasibility of the logistics model, we studied a possible earthquake in the real-world network of Isfahan. The province of Isfahan is one of the most historical and industrial regions of Iran. Location of Isfahan on the map and the major active faults of Iran are shown in Figure 11(a). According to geographical research in a 100-kilometer radius of Isfahan city, 63 faults were found nearby this city; some of these faults were hundreds of kilometers in length [36]. The fault density map and earthquakes detected in the mentioned region in 2000–2003 are displayed in Figure 11(b).Figure 11
(a) Major active faults and location of Isfahan in the map of Iran; (b) satellite image for fault density in 100 km radius of Isfahan and detected earthquakes in 2000–2003.
(a)(b)In this study, as illustrated in Figure12, the underlying transportation network comprises 6 shelters and one distribution center. The distribution center is situated in Bakhtiar-Dasht, and also the shelters are in Isfahan’s disparate areas. To deliver emergency packages from the depot to the different shelters, three similar vehicles with a capacity of 2000 exist. Required information regarding the problem is provided in Table 7. The shelters’ demands are assessed according to the polls provided by area’s crisis management experts. Times of service and priorities are set based on the nodes’ demands in accordance with two formulations Si1=1.5h+0.05×∑ωDiω/2 and Si2=2h+0.05×∑ωDiω/2 and wi=10×∑ωDiω/2×∑ω∑jDjω, respectively. In order to study the case study in a larger scale, we added fourteen potential shelter locations by consultation with experts (including some existing stadiums, schools, and college locations in Isfahan), and their related demands are estimated on the interval of [500, 900]. Other required parameters are also computed based on the abovementioned formulations.Figure 12
Traffic zones with various traffic patterns and representation of the logistics network (Google Maps).Table 7
Information regarding the main shelters in the network.
SheltersLocationsDemandsPrioritiesService times (hour)Scenario 1Scenario 2Scenario 1Scenario 21Bagh-Fadak105012502.03542.272.572Pirouzi Stadium70011001.59292.152.453Enghelab Stadium6608401.32742.072.374Ghadir Park100011001.85842.222.525Vaziche Stadium7508501.41592.102.406Abrisham85011501.76992.202.50The reliability of links between nodes and the speed pattern in the disaster aftermath are set in accord with deliberation with specialists in crisis management and civil architecture fields. We classified the region based on three traffic zones of heavy, medium, and light to determine the traffic pattern, which are shown in Figure12. This classification is performed based on existing transportation statistics of Isfahan. For each category of traffic, Figure 13 depicts the speed patterns for four time intervals [0, 3], [3, 6], [6, 9], and more than 9. It is notable that in some cases a route from one node to another one may include various traffic regions. In this situation, it is required that some breakpoints be specified in the indicated route where the vehicles enter to another traffic area. Hence, by considering these breakpoints, travel times are calculated using the method of Section 3.1.Figure 13
(a) Speed changes of the light traffic zone; (b) speed changes of the medium traffic zone; (c) speed changes of the heavy traffic zone (continues line: first scenario; dash line: second scenario).
(a)(b)(c)Owing to the urban networks’ structure, it is evident that the distances among nodes are not symmetric. In this regard, two alternative links employed between each pair of nodes are obtained from Google Maps. We report the details of links’ distances and their reliabilities between the main shelter locations in Table8. For the rest of the demand points, detailed information is not presented for the sake of brevity.Table 8
Details of distances (km) for two parallel links between main shelter locations in multigraph.
Depoti1i2i3i4i5i6Depot×20.1-122.8-123.7-125.7-116.2-129.6-0.9524.9-125.5-123.9-131.4-116.4-133.1-0.95i123.3-1×11.4-120.7-0.9514.3-120.1-130.1-0.9526.9-115.9-121.1-0.9515.9-121.2-137.3-0.95i226.0-112.4-1×7.90-0.953.2-119.2-129.2-0.9531.4-115.2-113.9-0.953.9-124.0-132.3-0.90i324.8-0.9519.5-0.959.0-0.95×6.7-0.9514.0-0.9515.2-137.9-0.9522.2-0.9514.9-0.957.5-0.9020.7-0.8525.8-1i425.7-116.3-12.6-17.8-0.95×25.1-124.4-0.9531.9-120.3-13.0-113.7-0.9529.9-0.932.1-0.95i517.0-118.6-120.0-113.5-0.9025.8-1×17.2-0.9017.6-119.9-120.5-114.6-0.9526.3-123.6-0.95i617.3-0.9030.1-0.9525.9-0.9516.3-123.6-0.9518.8-0.90×35.5-0.9546.9-0.9530.7-0.9524.6-131.3-0.9524.8-0.95We solved the proposed model once for a logistics network which includes only the main shelter locations (Case I) and once for a logistics network considering both main and potential shelter locations (Case II). By solving the proposed model, a series of solutions as a Pareto front are achieved. In Figures14 and 15, we depicted the obtained results for Case I and Case II, respectively. As can be seen from the figures, increasing the reliability of routes yields to increase the first objective function. It accentuates the conflict of objectives. Herein, a decision maker should select one of the solutions in the Pareto front according to his preference. By comparing Figures 14 and 15, it becomes clear that by increasing the number of demand points in Case II, the overall reliability (objective 2) of routes decreases and cost of traveling (objective 1) increases since a vehicle should serve more demand nodes.Figure 14
Approximate Pareto optimal for the studied Case I.Figure 15
Approximate Pareto optimal for the studied Case II.


## 7. Conclusion

In this study, we developed a novel multiobjective stochastic mathematic model for urban distribution of emergency supplies named SR-TDRRPM-SSTW. The first objective strived to minimize the total operational costs and penalties of late arrivals while the second one attempted to maximize the reliability of constructed routes. The arrival times at each demand node of the designated routing scheme were supposed to adhere to semisoft time windows. Since the SR-TDRRPM-SSTW was studied in a multigraph setting, not only should it determine the sequence of serving affected areas but also should decide about effective links between nodes. In this regard, in the first stage of the model, the vehicle routing decisions were performed, and based on expected conditions, the decisions about suitable paths and arrival times were made in the second stage with a finite number of defined scenarios.To solve the problem, an improved version of NSGA-II was proposed and employed. A comparison between the proposed NSGA-II and its standard version implied the effectiveness of the proposed method in several aspects. The employed statistical tests demonstrated that the difference between these two methods was statistically significant. In further analysis, we found that consideration of nondominated parallel links as a multigraph network can remarkably improve the solution of SR-TDRRPM-SSTW. The extensions for future works can include investigating the multicommodity version of SR-TDRRPM-SSTW and considering the possibility of split deliveries.



---

*Source: 8871952-2021-09-20.xml*