The Emergency Vehicle Routing Problem with Uncertain Demand under Sustainability Environments

The reasonable utilization of limited resources is critical to realize the sustainable developments. In the initial 72-h crucial rescue period after the disaster, emergency supplies have always been insufficient and the demands in the affected area have always been uncertain. In order to improve timeliness, utilization and sustainability of emergency service, the allocation of the emergency supplies and the emergency vehicle routes should be determined simultaneously. Assuming the uncertain demands follow normal distribution, an optimization model for the emergency vehicle routing, by considering the insufficient supplies and the uncertain demands, is developed. The objective function is applied to minimize the total costs, including the penalty costs induced by more or less supplies than the actual demands at all demand points, as well as the constraints of the time windows and vehicle load capacity taken into account. In more details, a solution method for the model, based on the genetic algorithm, is proposed, which solves the problem in two stages. A numerical example is presented to demonstrate the efficiency and validity of the proposed model and algorithm.


Introduction
In recent years, around the world, natural disasters, such as earthquake, hurricane, volcanic eruption and flood, are occurring increasingly frequently.More than 500 disasters are estimated to strike our planet every year, killing around 75,000 people and impacting more than 200 million others [1].For example, the Wenchuan Earthquake in China left more than 110,000 casualties in 2008, including 69,197 confirmed dead and 374,176 injured, with 18,222 listed as missing.The Haiti earthquake in 2010 caused 230,000 death and more than 300,000 injuries.
The enormous scale of these disasters has called attention to the need for effective and sustainable management of the relief activities after disasters.To mitigate damage and loss in disasters, issues in pre-disaster, during-disaster and post-disaster have been widely studied in the past few years.For relief activities in post-disaster, to transport sufficient emergency supplies to affected areas in order to support basic needs for those trapped in disaster-affected areas is a critical challenge [2].Emergency supplies that are essential for human survivals include water, food and medicines.After the disaster, the supplies must be delivered to the demand points immediately, especially in the initial rescue period, which is the optimal rescue period for the people in the affected areas.Unfortunately, emergency resources are always restricted in this period.Therefore, it is necessary to enhance the sustainability of the emergency service to make full use of the resources.Thus, the effectiveness and timeliness of the delivery are significant to the humanitarian response and the sustainability of the rescue operation, which provide the motivation for this work.
Disasters always result in massive demands that often outstrip available resources.The process of planning, managing, and controlling the flow of those resources to provide relief to affected people is called emergency logistics.Different from general business logistics, emergency logistics is unique in some aspects that may increase the relative complexity in solving the induced problems.First, the demand-related information, e.g., the density of demand in the affected points, is quite limited in the initial rescue period, and is intuitively unpredictable using historical data.The lack of relevant information will lead to difficulties to determine the actual demands of the areas.Second, relative to general business logistics, the emergency resources may be insufficient to decision makers in the initial rescue period, since the disaster is unpredictable.This fact adds more challenging issues to the management of emergency logistics.Consequently, in order to improve the utilization-efficiency and delivery-timeliness of the limited resources in the initial rescue period, the distribution of emergency supplies and emergency vehicles routing has to be determined synchronously under the uncertainty conditions.Emergency logistics has emerged as a worldwide noticeable theme as disasters, either artificial or natural, may occur anytime around the world with enormous consequences [3].
The crucial rescue period refers to the initial 72 h following the onset of a disaster, which is regarded as the most critical period to search and rescue the trapped survivals [4].In the crucial rescue period, the available emergency resources, including vehicles and supplies, are always limited.How to establish an effective emergency logistics system after emergency events, and launch the emergency plan to distribute the limited resources to the affect areas with a good planning and scheduling transportation, has been a research hotspot in recent decades.
The operations of emergency logistics are based on the timely delivery services and supplies allocation.In most cases, the vehicles are the most frequently used tools to implement the delivery work.Therefore, how to assign the work to existed vehicles and how to select the vehicle routes, namely, the emergency vehicle routing problem (VRP), are the key to the entire emergency logistics operations.Thus far, there have been many relevant studies on the VRP.Lei et al. [5], Zhu et al. [6], Errico et al. [7] and Zhang et al. [8] considered the VRP with time windows and stochastic service times.Yin and Chuang [9] first addressed a previously studied problem, then proposed a new problem version for planning the least cost routing, in which the deployed green vehicles transport the final products from suppliers to customers through a cross-dock subject to a CO 2 intensity constraint, and developed an adaptive memory artificial bee colony algorithm to tackle both problems.In their recent study [10], they proposed a bi-objective mathematical formulation for the cross-docking with interrelated operations (vehicle routing and vehicle scheduling), and developed a cooperative coevolution approach consisting of hyper-heuristics and hybrid-heuristics for achieving continuous improvement in alternating objectives.Marinakis et al. [11] introduced a new hybrid algorithmic approach based on particle swarm optimization for the vehicle routing problem with stochastic demands that are known only when the vehicle arrives to the customers.This algorithm is a combination of the particle swarm optimization algorithm with the 2-opt and 3-opt local search algorithms and with the path relinking strategy.Euchi et al. [12] described hybrid iterated density estimation evolutionary algorithm with 2-opt local search to determine the specific assignment of each tour to a private vehicle (internal fleet) or an outside carrier (external fleet).
Furthermore, Parragh et al. [13,14] wrote a survey about pickup and delivery problems, which are distinguished by two problem classes: the first class deals with the transportation of goods from the depot to line haul customers and from backhaul customers to the depot, and the second class refers to all those problems where goods are transported between pickup and delivery locations.Berghida and Boukra [15] dealt with the vehicle routing problem with simultaneous pickup and delivery Sustainability 2017, 9, 288 3 of 24 by a new technique based on a cooperative approach reaping the benefits of three metaheuristics.Maquera et al. [16] presented an application of scatter search to the problem of routing vehicles that are required to deliver and pickup goods, which generates all pairs of solutions that contain at least one new element.
However, as we mentioned above, the emergency VRP is different from general VRP in demand-related information and supply of emergency resources, which bring more complexity to the problem.Wohlgemuth et al. [17] considered the dynamic vehicle routing problem with anticipation in disaster relief.Pillac et al. [18] classified routing problems from the perspective of information quality and evolution, introduced the notion of degree of dynamism, and presented a comprehensive review of applications and solution methods for dynamic vehicle routing problems.While Cao et al. [19] proposed four robust strategies to cope with the uncertain demand of an open VRP, Moghaddama et al. [20] investigated a variant of an uncertain VRP that the customer's demands are supposed to be uncertain with unknown distributions.In addition, an advanced particle swarm optimization (PSO) algorithm was proposed to solve such VRP with a novel decoding scheme.
Despite the critical importance of emergency logistics distribution, only a limited amount of related studies have been carried out.Sheu [21] presented relief-demand management model for emergency logistics operations under imperfect information conditions.Lin et al. [2] proposed a logistics model for delivery of prioritized items in disaster relief operations.Chang et al. [22] gave a novel algorithm capability of regulating the distribution of available resources and automatically generating a variety of feasible emergency logistics schedules for decision-makers.In addition, Wang et al. [23] constructed a nonlinear integer open location-routing model for relief distribution problem in emergency logistics, considering travel time, the total cost, and reliability with split delivery.
The VRP is a problem with NP-difficulty, and this nature makes the exact algorithms only for small problems.For larger instances, heuristics is a natural choice.Therefore, much attention has been focused on designing heuristics algorithms with good performance.Following the first heuristics algorithm presented by Hochbaum and Shmoys [24], there has been a long list of work on designing heuristics algorithms for this problem over the years.Several meta-heuristic approaches were reported, including Genetic Algorithm (GA) [25][26][27], Tabu Search (TS) [28][29][30], Particle Swarm Optimization (PSO) [11,20,31], hybrid algorithm [32][33][34][35], etc., which were used to solve VRP and its variant problems.With reference to Genetic Algorithm, Ombuki et al. [25] proposed an application of genetic algorithms approach for multi-depot vehicle routing problem(MDVRP) using the Pareto ranking technique.An advantage of this approach is that it is unnecessary to derive weights for a weighted sum scoring formula.Furthermore, Ursani et al. [26] introduced the Localized Optimization Framework (LOF) which is an iterative procedure between two phases, Optimization and De-optimization.They have also chosen a genetic algorithm as an optimization methodology.It was demonstrated that the LGA is, on average, able to produce better solutions than most of the other heuristics on small scale problems of VRPTW.Moreover, the LGA has attained several new best solutions on popular datasets.
In the existing research results, the goal is all about minimizing the number of vehicles, vehicle travel distance and some other economic factors.However, this approach does not work with a full consideration of the distinguishing features in the crucial rescue period of relief operation, which may be in short supply and with uncertain demands.The main research of this paper is to fully consider the supply/demand features in crucial rescue period, and the optimization model for the emergency vehicle routing problem is developed.As for the solution methods, the particularity and practicality of the emergency vehicle routing problem are taken into consideration.This fact divides the problem into two stages, in which the first stage is to determine the serviced demand point set of each center, and the second stage is to optimize the used vehicle number, the actual supplies and vehicle routes at each center.
This paper is organized as follows.The mathematical model for the problem and the solution methods for the model are proposed in Section 2. In Section 3, numerical examples are presented to demonstrate the benefit of our new logistics model.Finally, the conclusions and suggestions for future work are given in Section 4.

Model Foundation and Notations
A well-known challenge to efficiently execute disaster relief operations is to obtain sufficient and accurate demand information from the disaster-affected area.Due to the chaotic situation in the critical period, the information is always incomplete and inaccurate.Therefore, it becomes very difficult to correctly estimate useful information to deliver the emergency supplies.As a result, the precise emergency demands are uncertain due to the unpredictability of the actual demands at the demand points.In this paper, we assume that the demand at each point is uncertain with a known probability distribution.For example, we could define the actual demand at demand point k as d k , which is a random variable with probability density function given by f k (x).The probability distribution of demand is assumed to be available, which could derive from population data and/or information gathered during the disaster prevention phase [36].
Similar to the general vehicle routing problem, we also assume that all the emergency vehicles are the same type.In the initial period of rescue, a variety of rescue facilities and resources are very limited, so it is necessary to take as many emergency demand points as possible into account in this period, provide some necessities for them to ensure that they could wait for the subsequent rescue.In this case, we assume that the delivery for each demand point can be only served by one vehicle one time, which should return back to the starting point (service central) after finishing the current rescue mission.
For the sake of simplicity, the service time of vehicles on the demand point is included in travel time.The affected points considered in this paper are connected by existing roads, which can be serviced by vehicles.Other isolated points could be serviced by other methods, such as helicopters.
In formulating the mathematical model, the following notations in Table 1 are used in this paper.The decision variable is defined as follows: 1 if node j is the next node to node i on the route of vehicle v ; 0 otherwise; i, j ∈ N, v ∈ V

Mathematical Model
In this section, the mathematical model for the emergency vehicle routing problem and insufficient supplies is constructed, based on the above assumptions and notations.
As mentioned previously, emergency supplies are often insufficient in the 72-h crucial rescue period after disasters.Due to the significance of the availability of emergency resource at the demand points, a relatively large penalty τ − k is associated with the shortage of a unit of supply at demand point k, and this penalty corresponds to the social cost of a death or a severe injury of a patient because of a supply shortage.On the contrary, since the total amount of emergency supplies are limited, some demand points may suffer from a shortage, while others are assigned even more supplies.Henceforth, the penalty τ + k is also assigned to a possible supply surplus, in terms of units.In other words, the amount of emergency supplies distributed to each demand point should be consistent with that of actual demands as soon as possible, and both the excessive and insufficient distribution could effect on rescue.Similar examples of penalty costs, due to excessive supplies, as well as shortages, can be found in the literature [37,38].These penalties should be considered in the objective function of the model to measure the utilization efficiency of the limited supplies.
Specifically, for any point k ∈ DN, when the delivered emergency supply is not equal with its actual demand, the amount of shortage ∆d − k and surplus ∆d + k can be calculated, respectively, as follows: According to the definition of notations, the probability of the supply shortage (P − k (z k )) and supply surplus (P + k (z k )) at demand point k can be expressed, respectively, as: Hence, the expected amount of the shortage (E ∆d − k ) or surplus (E ∆d + k ) at demand point k can be calculated by: Based on the above quantities, the expected penalties incurred by the distribution decision due to the shortage and surplus of the relief supplies at each demand point can be easily calculated.
In addition, the available rescue vehicles are always limited in the 72-h crucial rescue period, so it is better to finish the rescue work with fewer vehicles.For this purpose, the operating costs of vehicles should be measured in the objective function of the model, which could be described as: According to the above analysis, the optimization model of the emergency vehicle route considering insufficient supplies can be established as follows: In this model, the objective function (Equation ( 8)) is to minimize the total expected cost, in which is the total expected shortage penalty incurred by the insufficient supply at demand points, is the total expected surplus penalty incurred by the excessive supply at demand points, ij is the total operating costs of vehicles.In addition to the objective, Constraint (9) represents the load capacity constraint of a vehicle.Constraint (10) ensures that the vehicles which depart from a service central are no more than the maximum available vehicles in the center.Constraint (11) determines the time that the vehicle arrives at the demand points.Constraint (12) gives the time that the vehicle arrives at the demand point, in which the service time window must be satisfied.Constraint (13) ensures the emergency vehicles cannot travel between the service centers.Constraint (14) indicates that all the vehicles should start the delivery work from a service center and return to the same center after these missions are finished.Constraint (15) ensures that one demand point can be serviced only by one vehicle.Constraint (16) represents the domain of the decision variable.

Solution Method
When implementing emergency rescue after the disaster, all kinds of emergency resources are required to arrive within a specified period of time to the emergency demand points.Considering the particular situation of the traffic organization after the disaster, it is reasonable to transform the time constraint into the distance constraint.In other words, in order to ensure a timely rescue, distance between emergency supply center and demand point should be inside a reasonable radius.In turn, we could determine the served demand point set of each supply center on the basis of time constraint of the problem, and multi-depot vehicle routing problem of emergency resources can be transformed into single depot vehicle routing problem of emergency resources.

A Combination with Nearest Assignment and Average Distance
According to the particularity and practicality of emergency rescue after disaster, the problem can be divided into two stages, in which the first stage is to determine the served demand point set of each center, and the second stage is to optimize the used vehicle number and vehicle routes at each center.
With the assumption that the speed of emergency vehicles is a constant, the delivery time could be calculated by the delivery distance.In addition, the constraint of service time windows can be transformed to the constraint of distance.Base on this assumption, the methods, including Nearest Assignment and Average Distance, are used to transform the above multi-depot (center) vehicle routing problem to multiple single-depot (center) vehicle routing problem, so that the service sets of each center are determined.
For the Nearest Assignment method, the ratios and differences between the distance from the demand point to the nearest center (depot) and that from the one to the second-nearest center (depot) are calculated in the first step.The half of the distance between the two nearest centers is selected as metric (or determined by many experimental results) for the above differences, which is used to allocate the demand points to the centers.The other demand points (not including the boundary points) are serviced by the nearest center, respectively.Based on the initial assignment results of the Nearest Assignment, the Average Distance method is used to assign the boundary points into the service set of service center.

The Improved Genetic Algorithm of the Proposed Problem
To determine the used vehicle numbers and vehicle routes at each center, the heuristics algorithm is used.In general, all the heuristics algorithms are based on the local search method, while with different characteristics.TS method traverses the solution space and finds the optimal solution by neighborhood search, which leads to a relatively faster computation, while the quality of solution depends heavily on the starting point or the given initial solution.SA is similar to TS; it traverses the search space by testing random mutations on an individual solution, which has the advantage of robustness on the initial solution, while the searching efficiency is low.GA is a related global optimization technique based on the stochastic search, which is easy to handle arbitrary kinds of constraints and objectives.Because of these strengths, GA always ensures both solution quality and computational efficiency.As a result, the solution method is developed based on the GA in this work.
GA, a highly parallel, random and adaptive optimization algorithm, was proposed initially by Holland in 1975.The idea was derived from the nature of biological genetics "survival of the fittest", opened and enlightened by biological evolutionism.The detailed GA steps can be found in References [25][26][27].Here, we focus on the designing of the key parameters.
(1) Coding and decoding of the solution The encoding method of chromosome is a very significant work, which directly determines the degree of difficulty and the algorithm performance.The emergency vehicle routing problem, which was proposed in Section 2, is no longer a simple vehicle routing problem, for its uncertain demands/supplies.Therefore, the chromosome contains two parts of information, the driving path of the vehicle and the actual supplies.The more difficult part is how to represent the actual supplies.
The chromosome encoding of vehicle path is a random arrangement of natural number from 1 to L, its length.Each gene corresponds to a natural number, the gene value represents demand point, and the chromosome sequence represents the order of demand points which it is served.The actual number of emergency supplies uses binary 0-1 coding, and the decoding method of this encoding is described as follows: according to the capacity constraints of the vehicle, following the emergency service order, combined with the distribution of the number and the actual arrival time windows, put the demand point into the distribution path successively.For example, if the supplies exceed the capacity requirement of the vehicle or the time window at the demand point l is not satisfied, then we end the vehicle path, and a new vehicle routing from the demand point l is established.Repeat the above operation, until the distribution of all demand points is achieved; finally, we get a complete distribution plan.

A. Chromosome coding method
This paper encodes vehicle routes with natural numbers (the number of the demand points), while adopting binary 0-1 code for the actual supplies in each demand point.The coding of vehicle routes is similar to the present study [25], and the following will focus on the coding of the actual supplies.
For the point i, we assume that its uncertain demand is distributed in the interval (a i , b i ), then the coding for the quantity of actual supplies (setting accuracy to be 0.1) could be obtained by the following steps.
In turn, we determine the number of binary digits through the following formula, , then the actual supplies adopts n-bit binary encoding.

B. Chromosome decoding method
Denoting B as the decimal numeral of the above n-bit binary code, the decoding formula therefore becomes: For example, we assume there are 5 demand points in the problem, and it is determined to adopt 5, 4, 4, 5, 5 bit binary code to indicate the actual supplies, respectively.In turn, using the above decoding formula, the chromosome could be 43251|01011011001111010111010, which implies that the vehicle route is 4→3→2→5→1, and the actual supplies equal 6.7, 5.4, 4.5, 6.4 and 5.7, respectively.
The specific steps of this decoding method are outlined as follows: Step 1. Regard the first emergency demand point on the chromosome as the first service object to be added to the distribution path, with the two constraints, namely capacity constraint and time window constraint, to determine whether the hybrid algorithm actual supplies of the demand point do not exceed the vehicle capacity and meets the request of time windows.If both of them meet the requirements, then enter Step 2; otherwise, at least one does not meet the requirements, then enter Step 3.
Step 2. Add its next demand point into the path, continue to implement the above operation, namely to determine whether the sum of the actual supplies of the new demand point and the previous demand point exceeds the vehicle capacity, as well as whether the time window is satisfied.If both of them meet the requirements, then repeat this step; otherwise, at least one does not meet the requirements, then enter Step 3.
Step 3. From this demand point, reopen a new path, and regard the demand point as the first demand point on the new path, repeat Step 1.
Step 4. Repeat Step 1 to Step 3, until the last demand point is added to the path.A feasible solution to vehicle routing problem is obtained.
(2) Initial population In order to maintain the diversity of individuals in the initial population and prevent the GA into the local optimum, this paper generates the initial population by random methods.
(3) Selection, crossover and mutation strategy The function for individual fitness measurement is called the fitness function, which is the standard quantity to measure the quality of individual in population, and is the only basis to determine the further search direction and search range.Similar to the previous work, the reciprocal of the objective function value is used as the fitness function in this paper, as follow: where f i is the fitness value of chromosome i, and ξ i is its objective function value.
Sustainability 2017, 9, 288 9 of 24 In the genetic operation, the selection probability determines the probability of individual inheritance to the next generation.The calculation formula proposed by Michalewicz [39] is used here: where c is the selection probability of the individual which sort first in the population, c = 0.1 1−0.9 s , s denotes the population size, and p i is the selection probability of chromosome i.
We use the method based on the combination of the roulette wheel selection mechanism and the best individual preservation strategy to implement the selection operation.
The type of crossover we choose relies heavily on the type of genotype representation at hand; according to the actual situation of the problem itself, we use partially matched crossover for crossover operation.Specific operations are similar to the one as follows: where is the fitness value of chromosome , and is its objective function value.In the genetic operation, the selection probability determines the probability of individual inheritance to the next generation.The calculation formula proposed by Michalewicz [39] is used here: where is the selection probability of the individual which sort first in the population, = .
. , denotes the population size, and is the selection probability of chromosome .We use the method based on the combination of the roulette wheel selection mechanism and the best individual preservation strategy to implement the selection operation.
The type of crossover we choose relies heavily on the type of genotype representation at hand; according to the actual situation of the problem itself, we use partially matched crossover for crossover operation.Specific operations are similar to the one as follows: Before crossover: Parents: 42351||01011011001111010111010 53142||10001010101001010101110 After crossover: Children: 42141||01011010101001010111010 53352||10001011001111010101110 Mapping set: 5→4, 3→1 Replacement: Children: 52143||01011010101001010111010 41352||10001011001111010101110 As far as mutation operation is concerned, this paper adopts inversion mutation, combined with the characteristics of practical problems, to operate chromosomes.
The crossover probability and the mutation probability have a direct effect on the performance of GA.The parametric adaptation mechanism is introduced in this paper.In other words, and are shifted with the variation of individual fitness and the operation procedure, which are listed as follows: where , , , are the initial parameters (in this paper, = 0.5, = 0.05, = 1, and = 0.1); is the fitness value of best individual in current population; ̅ is the averaged fitness value of the population; is the larger fitness value of the two individuals participating in intersection; and is the fitness value of varied individuals.It is observed that higher individual fitness values lead to lower and values, and, when the fitness value comes close to , the value of and are close to zero.

Example 1
A numerical example is presented in this section.The demand points and service centers are shown in Figure 1, in which there are three service centers (A, B, and C) and 35 demand points.The positions (represented as two-dimensional coordinate) of centers and the number of available vehicles in each center are given in Table 2.
As far as mutation operation is concerned, this paper adopts inversion mutation, combined with the characteristics of practical problems, to operate chromosomes.
The crossover probability p c and the mutation probability p m have a direct effect on the performance of GA.The parametric adaptation mechanism is introduced in this paper.In other words, p c and p m are shifted with the variation of individual fitness and the operation procedure, which are listed as follows: where k 1 , k 2 , k 3 , k 4 are the initial parameters (in this paper, k 1 = 0.5, k 2 = 0.05, k 3 = 1, and k 4 = 0.1); f 0 is the fitness value of best individual in current population; f is the averaged fitness value of the population; f is the larger fitness value of the two individuals participating in intersection; and f is the fitness value of varied individuals.It is observed that higher individual fitness values lead to lower p c and p m values, and, when the fitness value comes close to f 0 , the value of p c and p m are close to zero.

Numerical Example 1
A numerical example is presented in this section.The demand points and service centers are shown in Figure 1, in which there are three service centers (A, B, and C) and 35 demand points.The positions (represented as two-dimensional coordinate) of centers and the number of available vehicles in each center are given in Table 2.In this example, we set the size of the population and maximum generation to be 50 and 1000, respectively.In addition, we assume that the emergency demand at each demand point follows normal distribution ( , ) within the given value interval, in which is the mean or expectation of the distribution, and is its standard deviation with its variance as .Then, we know that: These parameters, including , and value interval, are shown in Table 3.The positions and the time windows that must be serviced for 35 demand points are also given in Table 3.
The number of available vehicles in centers is given by = = = 5, with the maximum load = 27 ton and the average speed as 30 km/h.The fixed cost of vehicle per delivery is = 200 Yuan, ∀ .The transportation cost per kilometer is = 5 Yuan/km, ∀ , .When the delivered supplies do not meet the demand at the demand points, the penalty cost per supply unit for shortage at demand point is set as = 500 Yuan, ∀ , and the penalty cost per supply unit for surplus at demand point is given by = 300 Yuan, ∀ .
To transform the multi-depot problem into single-depot problem, the metrics of ratios and differences are valued as 0.65 and 22, respectively, in the Nearest Assignment method.The demand points with ratios less than 0.65 and differences more than 22 are chosen as boundary points.The results of the allocation according to the Nearest Allocation method are reported in Table 4.
Based on the initial demand allocation results in Table 4, the Average Distance method is used to reallocate the boundary points, including demand points 9, 11, 15, 21, 24 and 31.The average distance between the each boundary point and the demand points in each set is calculated.The point should be inserted into the set in which the average distance is shortest.As a result, the final allocation is captured.In this example, we set the size of the population and maximum generation to be 50 and 1000, respectively.In addition, we assume that the emergency demand at each demand point follows normal distribution N µ, σ 2 within the given value interval, in which µ is the mean or expectation of the distribution, and σ is its standard deviation with its variance as σ 2 .Then, we know that: These parameters, including µ, σ 2 and value interval, are shown in Table 3.The positions and the time windows that must be serviced for 35 demand points are also given in Table 3.
The number of available vehicles in centers is given by ϑ 1 = ϑ 2 = ϑ 3 = 5, with the maximum load Q = 27 ton and the average speed as 30 km/h.The fixed cost of vehicle per delivery is ω v = 200 Yuan, ∀v.The transportation cost per kilometer is c ij = 5 Yuan/km, ∀i, j.When the delivered supplies do not meet the demand at the demand points, the penalty cost per supply unit for shortage at demand point is set as ϕ − k = 500 Yuan, ∀k, and the penalty cost per supply unit for surplus at demand point is given by ϕ + k = 300 Yuan, ∀k.To transform the multi-depot problem into single-depot problem, the metrics of ratios and differences are valued as 0.65 and 22, respectively, in the Nearest Assignment method.The demand points with ratios less than 0.65 and differences more than 22 are chosen as boundary points.The results of the allocation according to the Nearest Allocation method are reported in Table 4.
Based on the initial demand allocation results in Table 4, the Average Distance method is used to reallocate the boundary points, including demand points 9, 11, 15, 21, 24 and 31.The average distance between the each boundary point and the demand points in each set is calculated.The point should be inserted into the set in which the average distance is shortest.As a result, the final allocation is captured.The initial solution is captured from random generation of initial population within each allocation set.Since the solution quality depends on the initial population in heuristic algorithm, ten computation tests are performed with different initial solutions.The proposed algorithm has been implemented in MATLAB and all the experimental tests were carried out on an Intel Pentium 3.20 GHz processor with 8GB of RAM (the same below).The maximum running time was limited to 12 s.The objective function values (OFV) of the corresponding ten optimal solutions are shown in Figure 2, where ϕ represents the OFV.
As demonstrated in Figure 2, the OFVs of the ten optimal solutions with different initial solutions are very close.In particular, the deviation between the maximum value (9789) and minimum value (9393) is only 4.21% different.Consequently, the solution method can be viewed to be convergent.As demonstrated in Figure 2, the OFVs of the ten optimal solutions with different initial solutions are very close.In particular, the deviation between the maximum value (9789) and minimum value (9393) is only 4.21% different.Consequently, the solution method can be viewed to be convergent.The convergence rate of the method is further described in Figure 3.We select the computation test with minimal OFV (9393) for the analysis.It can be found that the OFV decreases quickly with respect to the iteration of genetic algorithm.The optimal solution ( = 9393) can be obtained at the 197th iteration, in which there are nine vehicles in three centers for the emergency supply delivery.The optimal allocation scheme and vehicle routes are reported/described in Table 5 and Figure 4.The convergence rate of the method is further described in Figure 3.We select the computation test with minimal OFV (9393) for the analysis.It can be found that the OFV decreases quickly with respect to the iteration of genetic algorithm.The optimal solution (ϕ = 9393) can be obtained at the 197th iteration, in which there are nine vehicles in three centers for the emergency supply delivery.The optimal allocation scheme and vehicle routes are reported/described in Table 5 and Figure 4.As demonstrated in Figure 2, the OFVs of the ten optimal solutions with different initial solutions are very close.In particular, the deviation between the maximum value (9789) and minimum value (9393) is only 4.21% different.Consequently, the solution method can be viewed to be convergent.The convergence rate of the method is further described in Figure 3.We select the computation test with minimal OFV (9393) for the analysis.It can be found that the OFV decreases quickly with respect to the iteration of genetic algorithm.The optimal solution ( = 9393) can be obtained at the 197th iteration, in which there are nine vehicles in three centers for the emergency supply delivery.The optimal allocation scheme and vehicle routes are reported/described in Table 5 and Figure 4.As demonstrated in Figure 2, the OFVs of the ten optimal solutions with different initial solutions are very close.In particular, the deviation between the maximum value (9789) and minimum value (9393) is only 4.21% different.Consequently, the solution method can be viewed to be convergent.The convergence rate of the method is further described in Figure 3.We select the computation test with minimal OFV (9393) for the analysis.It can be found that the OFV decreases quickly with respect to the iteration of genetic algorithm.The optimal solution ( = 9393) can be obtained at the 197th iteration, in which there are nine vehicles in three centers for the emergency supply delivery.The optimal allocation scheme and vehicle routes are reported/described in Table 5 and Figure 4.It is observed in the optimal solution that the number of used vehicles at each center and the serviced demand points of each center are almost equal respectively.For each center, there are 11 or 12 demand points which should be served and three vehicles are used to deliver supplies.For each vehicle, there are three or four demand points which should be delivered.In addition, the demand points with bigger time window, which have more flexibility, are always in the back of the vehicle routes.These features demonstrate that the results are rational.

Numerical Example 2
Taking the waterlogging disaster in Changsha City of Hunan Province as an example, in Figure 5, the red dots and the blue dots represent the demand points and the supply points, respectively.The demand points and supply points in urban area are marked into coordinate axes.The specific data can be found in the Appendix A Tables A1 and A2.The maximum load of vehicles is Q = 30 ton and the average speed as 40 km/h.The fixed cost of vehicle per delivery is ω v = 200 Yuan, ∀v.The transportation cost per kilometer is c ij = 3 Yuan/km, ∀i, j.When the delivered supplies do not meet the demand at the demand points, the penalty cost per supply unit for shortage at demand point is set as ϕ − k = 100 Yuan, ∀k, and the penalty cost per supply unit for surplus at demand point is given by ϕ + k = 60 Yuan, ∀k.It is observed in the optimal solution that the number of used vehicles at each center and the serviced demand points of each center are almost equal respectively.For each center, there are 11 or 12 demand points which should be served and three vehicles are used to deliver supplies.For each vehicle, there are three or four demand points which should be delivered.In addition, the demand points with bigger time window, which have more flexibility, are always in the back of the vehicle routes.These features demonstrate that the results are rational.

Example 2
Taking the waterlogging disaster in Changsha City of Hunan Province as an example, in Figure 5, the red dots and the blue dots represent the demand points and the supply points, respectively.The demand points and supply points in urban area are marked into coordinate axes.The specific data can be found in the Appendix A Tables A1 and A2.The maximum load of vehicles is = 30 ton and the average speed as 40 km/h.The fixed cost of vehicle per delivery is = 200 Yuan, ∀ .The transportation cost per kilometer is = 3 Yuan/km, ∀ , .When the delivered supplies do not meet the demand at the demand points, the penalty cost per supply unit for shortage at demand point is set as = 100 Yuan, ∀ , and the penalty cost per supply unit for surplus at demand point is given by = 60 Yuan, ∀ .According to the method described in Section 2.3, we obtain the radius of the emergency service of each supply point.Then we determine the served demand point set of each supply center, as shown in Figure 6, which depicts the covered demand point within the scope of the emergency service by each demand point.The results of the allocation according to the Nearest Allocation method are reported as Appendix A Table A3.According to the method described in Section 2.3, we obtain the radius of the emergency service of each supply point.Then we determine the served demand point set of each supply center, as shown in Figure 6, which depicts the covered demand point within the scope of the emergency service by each demand point.The results of the allocation according to the Nearest Allocation method are reported as Appendix A Table A3.In this example, ten computation tests are performed with different initial solutions.The computation of the proposed algorithm is implemented in MATLAB.The maximum running time was less than 18 s.The optimal allocation scheme and vehicle routes are described in Figure 7.It is observed that the total traveled distance is 3713.9, and the total cost equals 18,217.7.From the data in Table 6, we can see, in the premise of different initial solutions, the total costs of 10 calculation optimal solutions is very close through MRGA algorithm, the relative deviation is only 4.7% to the maximum value (19,128.4) and the minimum value (18,217.7); in comparison with the initial solution, the use of MRGA algorithm to solve the problem has saved 16.01%-19.52%,and the computation time, mostly being within 24 s, also proved the efficiency of the proposed optimization model and MRGA algorithm.
Similar to example 1, the number of used vehicles at each center and the serviced demand points of each center are almost equal respectively.In most cases, there are 8-14 demand points which should be delivered by three vehicles.For each vehicle, there are almost three or four demand points which should be delivered except for a few vehicles services for two demand points and only one vehicle service for six demand points.In demand intensive areas, one vehicle could serve for more demand points, while a vehicle could provide emergency resources for two demand points at most in the remote area, where the demand is not concentrated.In general, the demand In this example, ten computation tests are performed with different initial solutions.The computation of the proposed algorithm is implemented in MATLAB.The maximum running time was less than 18 s.The optimal allocation scheme and vehicle routes are described in Figure 7.It is observed that the total traveled distance is 3713.9, and the total cost equals 18,217.7.In this example, ten computation tests are performed with different initial solutions.The computation of the proposed algorithm is implemented in MATLAB.The maximum running time was less than 18 s.The optimal allocation scheme and vehicle routes are described in Figure 7.It is observed that the total traveled distance is 3713.9, and the total cost equals 18,217.7.From the data in Table 6, we can see, in the premise of different initial solutions, the total costs of 10 calculation optimal solutions is very close through MRGA algorithm, the relative deviation is only 4.7% to the maximum value (19,128.4) and the minimum value (18,217.7); in comparison with the initial solution, the use of MRGA algorithm to solve the problem has saved 16.01%-19.52%,and the computation time, mostly being within 24 s, also proved the efficiency of the proposed optimization model and MRGA algorithm.
Similar to example 1, the number of used vehicles at each center and the serviced demand points of each center are almost equal respectively.In most cases, there are 8-14 demand points which should be delivered by three vehicles.For each vehicle, there are almost three or four demand points which should be delivered except for a few vehicles services for two demand points and only one vehicle service for six demand points.In demand intensive areas, one vehicle could serve for more demand points, while a vehicle could provide emergency resources for two demand points at most in the remote area, where the demand is not concentrated.In general, the demand From the data in Table 6, we can see, in the premise of different initial solutions, the total costs of 10 calculation optimal solutions is very close through MRGA algorithm, the relative deviation is only 4.7% to the maximum value (19,128.4) and the minimum value (18,217.7); in comparison with the initial solution, the use of MRGA algorithm to solve the problem has saved 16.01%-19.52%,and the computation time, mostly being within 24 s, also proved the efficiency of the proposed optimization model and MRGA algorithm.
Similar to example 1, the number of used vehicles at each center and the serviced demand points of each center are almost equal respectively.In most cases, there are 8-14 demand points which should be delivered by three vehicles.For each vehicle, there are almost three or four demand points which should be delivered except for a few vehicles services for two demand points and only one vehicle service for six demand points.In demand intensive areas, one vehicle could serve for more demand points, while a vehicle could provide emergency resources for two demand points at most in the remote area, where the demand is not concentrated.In general, the demand points with bigger time window are always behind the vehicle routes.The instance indicates that the proposed algorithm for the emergency vehicle routing problem is reasonable and efficient.In Table 7, we change the penalty costs ϕ − k and ϕ + k , and keep other parameters unchanged; five experiments are carried out for each degree and we take the average as a result.With the increase of penalty cost, the total cost of the system increases, on the contrary; the total cost of the system decreases with the decrease of penalty cost.In the range of more than 20% increase, the number of used vehicles on several demands points increase as well; on the other hand, the results remain unchanged in the case of reducing the penalty costs in the range of 30%; with the increase rate up to 40%, the number of used vehicles on several demands points decrease.As depicted in the Figure 8, the evolution trend between total cost and distance is very closely.In the range [−30%, 20%], the distance remain unchanged, as the total cost close to linear change and change slowly; while the change degree reaches −40% and 30%, the distance present a variation with a jump, as does the total cost.In Table 7, we change the penalty costs and , and keep other parameters unchanged; five experiments are carried out for each degree and we take the average as a result.With the increase of penalty cost, the total cost of the system increases, on the contrary; the total cost of the system decreases with the decrease of penalty cost.In the range of more than 20% increase, the number of used vehicles on several demands points increase as well; on the other hand, the results remain unchanged in the case of reducing the penalty costs in the range of 30%; with the increase rate up to 40%, the number of used vehicles on several demands points decrease.As depicted in the Figure 8, the evolution trend between total cost and distance is very closely.In the range [−30%, 20%], the distance remain unchanged, as the total cost close to linear change and change slowly; while the change degree reaches −40% and 30%, the distance present a variation with a jump, as does the total cost.Next, we change the time windows while keeping other parameters unchanged (it is assumed that all time windows in the system are changed at the same rate).With the increase of time windows, the results remain unchanged; on the contrary, the total cost and distance of the system increase with the decrease of time windows.In the range of more than 10% decrease, the number of used vehicles on several demands points increase; on the other hand, the results remain unchanged in the case of raising the time windows; when the change rate runs up to 50%, there are no feasible solutions in several demands points.The specific data are in Table 8.In addition, we can see from Figure 9 that the trend of change in cost and distance are consistent.Next, we change the time windows while keeping other parameters unchanged (it is assumed that all time windows in the system are changed at the same rate).With the increase of time windows, the results remain unchanged; on the contrary, the total cost and distance of the system increase with the decrease of time windows.In the range of more than 10% decrease, the number of used vehicles on several demands points increase; on the other hand, the results remain unchanged in the case of raising the time windows; when the change rate runs up to 50%, there are no feasible solutions in several demands points.The specific data are in Table 8.In addition, we can see from Figure 9 that the trend of change in cost and distance are consistent.considered.In addition, the number of vehicles is demonstrated in the results, so that the sustainability of emergency service is satisfied.The genetic algorithm described above is compared in Table 10 with some of the other classical metaheuristics [40,41] (BD,2003).The rightmost column describes the results of our proposed method, including total traveled distance and the number of vehicles, the computer and the CPU time are also described as reported by the authors.According to Table 10, the results of those experiments do not indicate any conclusive evidence to support a dominating heuristic over the others.Due to the difference and uncertainty of the emergency demand interval, the problem relative to the previous examples has increased its complexity, but, in general, the results by the proposed method prove to be competitive and effective, since it is the closest to the performance of best-known heuristic routing procedures.Moreover, it is demonstrated that the emergency vehicle routing problem can be tested through the VRPTW benchmark problems of Solomon in an acceptable way.

Conclusions
To improve sustainability of rescue operation after disaster, utilizing limited emergency resource reasonably was required.In this paper, the emergency VRP with insufficient supplies and uncertain demands during the crucial rescue period (the initial 72-h) after the disaster were investigated, hoping to enhance timeliness, utilization and sustainability of the emergency service.Combined with the characteristics of the crucial rescue period, such as the emergency resources are always limited, the uncertain emergence demands following the probability distribution is assumed at first.Based on this assumption, the expected penalty costs are incurred by the decision of resource distribution, due to the shortage and surplus of the relief supplies at each demand point.Then, the optimization model for emergency vehicles routing problem with insufficient supplies is proposed, in which the objective function was to minimize the transportation costs, vehicle operation costs and the penalty costs at all demand points.
The problem was divided into two steps.In first step, the Nearest Assignment method and Average Distance method were used to transform the multi-depot vehicle routing problem into multiple single-depot vehicle routing problem; meanwhile, the service point set of each service center was determined.Then, the heuristics algorithm, based on Genetic Algorithm (GA), was developed to solve the single-depot vehicle routing problems.In order to improve solution quality and efficiency, the hybrid encoding mode was used to construct genetic operator in the algorithm.
The numerical example demonstrated that the presented method would address the problem quickly and efficiently.The results illuminated that the model performed well in the disaster rescue operation, where appropriate amount delivery was momentous because only correct deliveries of supplies could make sure to keep the people in the affected area alive as long as possible and to wait for the subsequent rescue.In addition, our model was sensitive to the penalty cost settings, and this fact suggests that rational determinations of the penalty cost would make our model more efficient.
Many variants of the proposed model could be considered in the future research; for instance, a penalty for missing the time window could be introduced, different types of resources could be appropriately defined, etc.Further research could be carried out related to the proposed parameters and their performance.Moreover, applications of the proposed concept of penalty to other related problems are also possible in the sustainability environments.

Figure 1 .
Figure 1.The network of affected areas.

Figure 4 .
Figure 4. Optimal routings of the vehicles.

Figure 4 .
Figure 4. Optimal routings of the vehicles.

Figure 4 .
Figure 4. Optimal routings of the vehicles.

Figure 4 .
Figure 4. Optimal routings of the vehicles.

Figure 5 .
Figure 5. Coverage provided by supply points.

Figure 5 .
Figure 5. Coverage provided by supply points.

Figure 6 .
Figure 6.Assignments of the points.

Figure 7 .
Figure 7. Optimal routings of the vehicles.

Figure 6 .
Figure 6.Assignments of the points.

Figure 7 .
Figure 7. Optimal routings of the vehicles.

Figure 7 .
Figure 7. Optimal routings of the vehicles.

Figure 8 .
Figure 8.Total cost and Distance vs.The penalty costs.

Figure 8 .
Figure 8.Total cost and Distance vs.The penalty costs.

Figure 9 .
Figure 9.Total cost and Distance vs.The time windows.

Figure 9 .
Figure 9.Total cost and Distance vs.The time windows.

Table 2 .
Information of service centers.

Table 2 .
Information of service centers.
Figure 1.The network of affected areas.

Table 3 .
Information of demand points.

Table 4 .
Initial allocation of the demand points to service centers.

Table 5 .
Vehicle routings and the actual service quantity.

Table 5 .
Vehicle routings and the actual service quantity.

Table 6 .
The initial solutions vs. the optimal solutions.
with bigger time window are always behind the vehicle routes.The instance indicates that the proposed algorithm for the emergency vehicle routing problem is reasonable and efficient. points

Table 6 .
The initial solutions vs. the optimal solutions.

Table 7 .
Total cost and Distance vs.The penalty costs.

Table 7 .
Total cost and Distance vs.The penalty costs.

Table 8 .
Total cost and Distance vs.The time windows.

Table 9 .
The performance of the proposed method.

Table 10 .
Average performance comparison among other metaheuristics.

Table A3 .
Allocation of the demand points to service centers.

Table A4 .
The performance of the proposed method to the problem (Solomon Benchmarks with narrow time windows).

Table A5 .
The performance of the proposed method to the problem (Solomon Benchmarks with wide time windows).