Abstract: Emergency materials scheduling is one of the key links in the process and also is a complicated system engineering that is influenced by multitudinous factors. In order to improve the efficiency of emergency rescue, dynamic change of demand function for emergency materials was analyzed according to the characters of emergency materials scheduling. On this basis, optimization model of emergency materials scheduling have been built to minimize the total casualty losses when the demand doesn’t meet and to minimize the total cost which is from distribution center of emergency materials to the rescue point. Considering that the encoder mode of the different evolution and particle swarm optimization algorithm is similarity and the search capability of two types of algorithm is very strong, the Different Evolution (DE) optimization algorithm and Particle Swarm Optimization algorithm (PSO) are combined into an efficient hybrid algorithm (H-PD) to solve this problem. Finally, a case study has been carried out in order to testify validity of this model and its algorithm by using calculating and comparing analysis. The results show that the hybrid optimization algorithm above mentioned can gain smaller loss and larger cost than greedy algorithm, the improved simulated annealing optimization algorithm is reasonable. This model and its algorithm can well reflect the laws of things development and can be of practical and effective value in actual project.
INTRODUCTION
With the rapid development of social economy, modern metropolitan areas is increasing more and more. The unexpected events were also increased. If there are no scientific scheduling schemes of emergency materials to deal with the unexpected events at any moment such as natural disasters, SARS, terrorist and so on, the people’s lives and property may be greatly harmed (Larson et al., 2006; Altay and Green, 2006).
Emergency materials scheduling is a complicated system engineering that is influenced by multitudinous factors, such as the type of relief materials, the degree of disaster, demand for materials and so on (Green and Kolesar, 2004). At present, world governments, sociology and relevant scholars have paid highly attention to this problem and the relevant scholars have done much study in this area. The number of retrieval depots for massive earthquake disaster are studied (Liu et al., 2001). The process model form arterial dispatching is designed and explained from the preparation (Tang and Zhang, 2009). Two groups of decision variables are defined to indicate the nodes and arcs on the path selected and an integer programming optimization model of the problem is constructed (He and Tao, 2009). The ELS3 theory was applied to analyze emergency materials scheduling and an optimization model for emergency materials scheduling was proposed considering the timeliness, economy and reliability (Han et al., 2009). A management and scheduling model for emergency logistics under the situation that supply points supplies are less than the supplies of emergency points is set up (Chen and Wang, 2010). Optimization scheduling model for emergency materials has been proposed based on rough set theory (Song and Liu, 2010). Model of emergency materials scheduling has been provided proposed based on the variable traffic network (Chen et al., 2010). A multi-objective mixed integer model was built for studying emergency road repairing and emergency supplies distributing after earthquake disasters (Chen and Shuai, 2012). Emergency cooperative distribution model under conditions of multiple periodic emergencies and the transverse and longitudinal joint distribution patterns are analyzed and discussed (Wang and Wang, 2012). A multi-object model was established considering that of uncertain (Tian et al., 2011).
Although, there are many studies on emergency materials scheduling problems, in general, the existing study on emergency material scheduling including static and dynamic scheduling.
Static scheduling concentrates the combinatorial optimization of the emergency supply point. The existing study supposes that the supply of emergency materials is infinite. Nevertheless, demand of emergency material is not the same in different stage of emergency rescue and emergency supplies are often in short at the early period of emergency rescue. Consequently, static scheduling cannot be applied to the large-scale emergency.
Dynamic scheduling considering emergency task and material characteristics at different stages of emergency rescue process, but the materials scheduling theory and method of large-scale emergency still needs to be explored. The existing study did not grasp the characteristics of large-scale emergency demand and process understanding of large-scale emergency scheduling is still absent and the corresponding algorithm cannot be applied to large-scale emergency scheduling, the application of the existing model and its algorithm need further improvement and perfection.
For this purpose, firstly, the demand function of emergency materials is analyzed considering that the influence of seasonal, regional and so on. Secondly, models have been built considering the multi-state of emergency rescue, multiple supply point, multiple demand point and limited resources. On this basis, DE and PSO algorithm are combined into an efficient hybrid algorithm to solve this problem.
DEMAND FUNCTION
There is a close relationship between the demand of the disaster points for emergency materials and survival rate without considering rescue. Therefore the initial demand for Emergency materials changes little, but the demand is at its peak next time, such as within 48 h after the disaster, due to victims’ physiological condition gradually deteriorated, such as hunger, pain and so on, at the same time, mental endurance gradually decreased too, therefore the sharp decline in the survival rate lead to a sharp decline in demand.
After the golden period of emergency rescue, the survival rate will gradually approach zero. Therefore, the demand for materials is almost close to zero too and this moment is defined as rescue termination time tS. The normal distribution function curve agree with this trend, so it is used to describe the demand of the disaster points for emergency materials, namely N~(0, σ2). According to the knowledge of probability distribution, the value of tS is 3σ. The value of σ is inversely proportional to the level of the disaster, the affected area and the population density. That is, the higher the level, the greater the scope of the impact, the greater the population density, the smaller the value of σ, meantime, the steeper the curve, the faster change in demand for the materials, the more urgent the rescue time. The function is introduced here, σ = f (α, β, γ), where α, β, γ measure the disaster level, the impact of the scope and the population density, respectively.
The initial demand function of the kth class materials has been identified by the above analysis and is denote as:
| (1) |
where, the initial demand function of the disaster points for materials kth is expressed as gk(1)t. After a period of response time t1, a number of materials kth (pk) is delivered to the disaster points. At this time, due to shortage of materials and capacity constraints of transport capacity, supply of emergency materials pk is less than demand gk(1)t of kth materials. Therefore, the next stage of the demand function will no longer be gk(1)t, but gk(2)t. The values of σ is determined by the initial demand of every stage. Every time when the materials is delivered to that point, the demand function of it is continuous on left and not continuous on right and casualty loss for the entire process is the area which is surrounded by the demand function curve and the axis. The dynamic change of demand function curve under conditions of single rescue point, single disaster point and multi-variety emergency material is shown as Fig. 1. The dynamic change of demand function curve under conditions of multi-rescue point, multi-disaster point and multi-variety emergency material is shown as Fig. 2.
According to the 3σ principle of probability and statistics that is to say, although the range of normal variables is (-∞, +∞), but it is almost certainly which the value falls within (μ-3σ, μ+3σ). The distribution of the demand curve is (0, +∞), but it is almost certainly which the value falls within (3σ, +∞).
| Fig. 1: | Dynamic change of demand function under conditions of single rescue point, single disaster point and multi-variety emergency material |
| Fig. 2: | Dynamic change of demand function under conditions of multi-rescue point, multi-disaster point and multi-variety emergency material |
Associated with the rescue process, 3δ is defined as tS, the best rescue time is (0, tS).
gk(2)t satisfy recursive conditions and is denote as:
| (2) |
Demand function gk(t) with uncertainty is a discontinuous piecewise recursive function and is determined by the initial function and emergency supplies distribution decision. Piecewise function recursive relationship is a good explanation for this problem.
ESTABLISH OF OPTIMAL MODEL
wkj denotes the preference coefficient of the region jth for the emergency materials kth. Characterizing the demand characteristics of the region jth, for instance, the ethnic minorities have different needs for different materials:
| wj | = | Vulnerability coefficient of the region jth. Quantify the degree of difficulty to lose, when the emergency materials does not meet the demand |
| wkf | = | Materials coefficient of materials kth in season f, the material demand varies with seasons |
| Q | = | Total cost of the emergency rescue |
| Z | = | Casualty losses |
| cij | = | Unit distance transportation cost which is from distribution center of emergency materials Ai to the rescue point Dj |
| dij | = | Transportation distance which is from distribution center of emergency materials Ai to the rescue point Dj |
| pijk | = | Supply of emergency materials kth which is from distribution center of emergency materials Ai to the rescue point Dj at the time tkij |
| xijk | = | 0-1 integer variables, as follows |
| sk | = | Production cost of materials k |
Setting conditions:
| • | Emergency material kinds of each emergency rescue point is enough, but number of emergency materials is limited |
| • | Initial demand of each disaster point for emergency materials have been predicted by experts according to the emergency level, population density and the scope of the influence and dynamic change of demand function |
| • | Transportation from each emergency rescue point to each disaster point has been determined and thus the transportation time is definite, as well as the transportation costs is definite |
| • | This model does not consider the secondary disasters (such as after shocks). |
Analysis of constraints:
| • | Under the conditions of large-scale emergencies, emergency materials are always scarce. So, the supply of emergency materials kth from Ai to Dj should be non-negative and can’t exceed the storage of materials kth of the rescue point Ai , namely: |
| (3) |
| • | Response time from the last rescue point Am to the disaster point Dj must be within effective rescue time tsj, Otherwise the rescue is invalid: |
| (4) |
| • | Supply of emergency materials k from Ai to Dj is subject to the limitation of the capacity of the path: |
| (5) |
| • | Supply of emergency materials k from Ai to Dj shouldn’t be exceeded the demand of the disaster point this moment, else that would be a waste of transport capacity and materials: |
| (6) |
Defined 0-1 variable:
Analysis of the objective function: First objective function is to minimize the total casualty losses when the demand doesn’t meet. Meanwhile, the second objective function is to minimize the total cost which is from distribution center of emergency materials Ai to the rescue point Dj.
The objective function to minimum the casualty losses is denote as:
The objective function to minimum the cost which is from distribution center of emergency materials Ai to the rescue point Dj is denote as:
is the total transportation costs:
is the total production cost. The sum of the two is the total cost:
SOLUTION BASED ON H-PD OPTIMIZATION ALGORITHM
When dealing with the problem of high-dimensional complex, PSO algorithm will inevitably cease at some local minimum point with a certain probability, how to avoid premature stagnation and enhance the searching ability is the key to solve complex problems. Considering that the encoder mode of the different evolution and particle swarm optimization algorithm is similarity and the search capability of two types of algorithm is very strong. For this purpose, the two kinds of optimization algorithm are combined into an efficient hybrid algorithm to improve the search efficiency and strengthen the search ability.
Meanwhile, in order to strengthen constraint boundary search and to obtain the global optimal solution with higher probability, more particles are pulled to constrained boundary by introducing appropriate infeasible solutions.
General principles of particle swarm optimization algorithm: The PSO algorithm is based on population, but it find the optimal solution through the collaboration between individuals.
Suppose that location of ith particle in N dimensional space is represented as a vector, flight speed is expressed as vector. Each particle has a fitness function decided by optimized objective and the best position and the current position are discovered, it can be regarded as the flying experience of particle:
where, i = 1, 2,…, M, is the total number of particles in the group. vidk+1 is the dth dimensional component of flight speed vector of the ith particle by kth iteration. xidk+1 is the dth dimensional component of position vector of the ith particle by kth iteration. pid is the dth dimensional component of position vector pbest of therefore ith particle.
| Fig. 3: | Diagram of particle motion |
pgd is the dth dimensional component of the best position gbest for the group. c1 , c2 are the weighting factors. Rand() is a random function. w is inertia weight function. Diagram of particle motion is shown as Fig. 3.
General principles of the different evolution: According to the Different Evolution (DE), in the N dimensional space, selecting NP vectors of N dimensional as a group, population size does not change in the whole operation process. The mutation operation is definite as follows:
| (7) |
where, r1, r2, r3ε{1, 2,…, NP} F is mutagenic factors which is belong to [1, 2].
The interlace operation of different evolution is definite as follows:
Each component is determined as follows:
where, ranb(j) is the random number which belongs to [0, 1]. CR is cross constants which belong to [0, 1]. ranb(i) is the random integer which belong to [1, N], in order to ensure that has at least one component from vi,G+1.
The greedy strategy is taken in selecting operation of the different evolution, namely, only when the offspring ui,G+1 is better than that of parent xi,G, then ui,G+1 is retained. Let, xi,G = ui,G+1, otherwise, xi,G is retained.
Variation principle of DE algorithm is shown as Fig. 4.
| Fig. 4: | Variation principle of DE algorithm |
Design of H-PD optimization algorithm: Give consideration to the efficiency of the algorithm and the groups search strategy of PSO algorithm, a part infeasible solutions closed boundary are introduced to compare their fitness value. Moreover, a certain proportion of better fitness infeasible solution are retained to attract the group moving to the constraint boundary and help groups to find better solutions.
According to the calculating size of this model above mentioned, only 15% individual selected from the particle swarm is used for DE search in this study. Some particles will deviate from its trajectory when groups gather to gbest. Search the neighborhood of its trajectory, not only can effectively avoid the premature stagnation particles, but also can strengthen the search range (Yang et al., 2009).
Defined function as:
where, v(x) which denote close degree between the particle and the boundary.
Death penalty function method is adopted in the process of the group of search to give directly poor fitness value for the infeasible solution. But every G generation, q particles are selected from a group of infeasible solutions according to the function v(x) and then only maintain q particles position of infeasible solution in the subsequent G generation computing, if there are other infeasible solution, poor fitness value is directly given. q particles are remove still using death penalty function method in the subsequent G generation computing and a part of infeasible solutions are retain. The iterative process of H-PD is shown as Fig. 5.
| Fig. 5: | Iterative process of H-PD |
The algorithm process of H-PD optimization algorithm is designed as follows:
| Step 1: | Initialization of parameters |
| Step 2: | If not meeting the stop conditions, then performance step 3~6. Otherwise, turn to step 7 |
| Step 3: | Iter ++ |
| Step 4: | Determining whether it should retain the infeasible solution. If retaining the infeasible solution, then turn to step 7. Otherwise, turn to step 5 |
| Step 5: | Searching by using the death penalty function method |
| Step 6: | Calculating V(x) and retain the infeasible solution |
| Step 7: | Output of the optimum results |
| Step 8: | Designed to push particle swarm left the infeasible region and search only in the feasible region. Step 6 is designed to attract part of the particle in the feasible region through the constraint boundary to flight the infeasible region. The search ability for near the boundary is strengthened through step 5 and 6 continuous interleaved execution |
EMPIRICAL ANALYSES
To verify the precision of the model in the previous sections, disaster which occurs in a region is chosen to verify the effectiveness of model and algorithm. This disaster includes three rescue point and two disaster point. Three rescue points supply emergency materials for two disaster points. The materials storage of the three rescue points has been given as Ai = (100, 110, 130; 120, 140, 110). The demand of each disaster point for two kinds of materials gj = (200, 100; 150, 180), Tij = (2, 1.5; 1, 2.5; 3, 3), the termination rescue time of two disaster points are 6 and 8 days, the vulnerability coefficients of the disaster points are 0.8 and 0.6 and the material coefficients are (0.5, 0.8), respectively. The first objective is to minimize the casualties loss caused by unmet the demand and the second is to minimize the total rescue cost. The topology of the rescue points and disaster points is shown as Fig. 6.
Design of parameters are as follows: M = 600, wmax = 0.9, wmin = 0.4, itermax = 3000, c1 = c2 = 2, F = 0.8, q = 6, CR = 0.3, G = 10. Independent operation 10 times, the optimal solution is selected and then the average optimal solution and variance are calculated.
| Fig. 6: | Topology of rescue points and disaster points |
| Fig. 7: | Variation relationship between CR and the evolutional generation |
Variation relationship between CR and the evolutional generation is shown, respectively as Fig. 7. Variation relationship between F and the evolutional generation is shown, respectively as Fig. 8.
As can be seen from Fig. 7 and 8, Crm (average value of CR) continues to increase with the evolution generation while Fm (average value of F) continues to decreases with increase of evolution generation.
Variation relationship between fitness and iterations under condition of different w is shown as Fig. 9.
As can be seen from Fig. 9, the H-PD algorithm is easy to fall into the local optimum when w is given as 0, yet convergence of the H-PD algorithm has certain difficulty when w is given as 0.8.
Thus, it can be seen, for the calculating scale of this model above mentioned, it is more appropriate when w is given as about 0.5.
The optimal solution based on H-PD algorithm is as the following Table 1.
As can be seen from Table 1, demand of various types of emergency materials for each disaster point has certain difference.
| Fig. 8: | Variation relationship between F and the evolutional generation |
The relationship between ideal solutions, unmet cumulative loss and learning factor is shown as Fig. 10. The relationship between ideal solutions, unmet cumulative loss and iterations is shown as Fig. 11.
| Fig. 9: | Variation relationship between fitness and iterations under condition of different w |
| Fig. 10: | The relationship between ideal solutions, unmet cumulative loss and learning factor |
| Fig. 11: | The relationship between ideal solutions, unmet cumulative loss and iterations |
| Table 1: | Optimal solution based on H-PD algorithm |
| Table 2: | Results based on the greedy algorithm |
In order to further illustrate the superiority of the algorithm proposed in this study, the greedy algorithm (Yan et al., 2010) is employed to compute, the results based on the greedy algorithm is as follows in Table 2.
As can be seen from Table 2, despite demand of various types of emergency materials for each disaster point has certain difference, but emergency materials which rescue point 3 supplied to disaster point 1 is 0. It seems inconsistent with the actual, also shows that the defect of the greedy algorithm.
Thus it can be seen that H-PD optimization algorithm can gain smaller loss and larger cost than greedy algorithm. In the emergency process, a smaller loss program is chosen. Therefore, the H-PD optimization algorithm above mentioned is reasonable.
ANALYSIS AND DISCUSSION
| • | Emergency materials scheduling is a complicated system engineering that is influenced by multitudinous factors, such as the type of relief materials, degree of disaster, demand for materials and so on. Demand function gk(t) is determined by the initial function and emergency supplies distribution decision. Piecewise function recursive relationship is a good explanation for this problem |
| • | Exploring ability of particle swarm is expected in the time of flying start. However, with the increasing of iterations, the developing ability of particle swarm is expected in the time of flying late. So the dynamically inertia weight is adjusted by the varying weight |
| • | Optimization model of emergency materials scheduling was built considering dynamic of emergency material demand and adapting to characteristics at different stages of emergency rescue process |
| • | Considering that encoder mode of DE and PSO is similarity and search capability of two types of algorithm is very strong, DE and PSO are combined into H-PD to solve this problem. From the result, individual diversity and global search capability should been keep in early stage, with increase of generations, algorithms need enhance local search ability. The precision of solution can be improved |
| • | Furthermore, comparing H-PD algorithm and the greedy algorithm applied to this problem, the minimum unmet loss is 350 and the total cost is 3041 based on H-PD optimization algorithm. While the total unmet losses accumulated by the greedy algorithm are 391.1 and the total cost is 2841. This show that the hybrid optimization algorithm above mentioned can gain smaller loss and larger cost than greedy algorithm, the improved simulated annealing optimization algorithm is reasonable |
| • | In addition, w plays a role to adjust the balance between individual in the global search and local search, therefore, for different objective functions, the ideal value of w is different and the selection of w should be dynamically adapted. Meanwhile, considering that w decreases linearly, global search and local search can be effectively combined with,but the decline rate and control of upper and lower limit need to be further tested |
| • | For the prevalent drawback of local minimum in potential field based navigation algorithm, certain heuristic techniques are required to be appended in some unfavorable cases. The mending methods will be presented in the future study |
| • | Because secondary disasters are likely to occur frequently and also the uncertainty of decision environment are likely to increase. Evolution mechanism of traffic network under the sudden disaster should be paid attention to. These are the mainly subject in the next phase |
CONCLUSION
Emergency materials scheduling is a complicated system engineering that is influenced by multitudinous factors, such as the type of relief materials, the degree of disaster and demand for materials and so on. At present, world governments, sociology and relevant scholars have paid highly attention to this problem and the relevant scholars have done much study in this area. In this study a class of multi-objective optimization model and its algorithm is studied taking the system losses into consideration, as well as the cost which has a strong practical application background. Comparing with the tradition method, this model above mentioned comprehensive considers the various influencing factors, it has certain practical significance.
ACKNOWLEDGMENT
This study is supported by the general planning project of humanities and Social Sciences from Ministry of education of China (No. 11YJCZH170). The authors are grateful for the anonymous reviewers who made constructive comments.