Dynamic Optimization Model and Algorithm Design for Emergency Materials Dispatch

Abstract

Emergency materials dispatch (EMD) is a typical dynamic vehicle routing problem (DVRP) and it concentrates on process strategy solving, which is different from the traditional static vehicle routing problem. Based on the characteristics of emergency materials dispatch, DVRP changed the EMD into a series of static problems in time axis. A mathematical multiobjective model is established, and the corresponding improved ant colony optimization algorithm is designed to solve the problem. Finally, a numeric example is provided to demonstrate the validity and feasibility of this proposed model and algorithm.

1. Introduction

In the emergency situations, such as earthquakes, snowstorms, and wars, ensuring all kinds of materials to be efficiently dispatched is the key problem to the emergency logistics. And it is also a typical dynamic vehicle routing problem whose aim is exactly a guarantee of the smooth logistic workflow process. Meanwhile, emergency material dispatch problem is a much more complex vehicle routing problem because of its high efficiency and strict time limitation. Therefore, exploring the modeling methods and efficient solving algorithms under the dynamic condition is greatly significant [1].

The research scope of dynamic vehicle routing problem (DVRP) includes uncertain demands, uncertain network performance, uncertain service vehicles, and subjective preference of decision makers. The solving strategy can be divided into two categories. One is local optimization strategy, which deals with the dynamic adjustment of implementation process and local optimization algorithms. The insertion method and k-opt method are often used to redesign the routing lines and vehicles distribution [2]. The other one is restarting optimization strategy, which is actually a static solution for dynamic problems. Once received determined real-time information, the strategy starts from the very beginning again to find the optimal material dispatch scheme. As one successful application of restarting optimization strategy, Paraftis used dynamic programming (DP) to deal with dial-a-ride problems with one vehicle [3]. However, this method can only solve small-scale problems not exceeding ten demand nodes.

Ant colony algorithm (ACA) has the nature advantage to adapt the dynamic environment with inherent robustness. Eyckelhof [4] carried out some research of VRP under the condition of unchanging node number and changing node distance, which indicates the sudden change of traffic congestion degree. The preliminary study results show that the original version and some slightly improved versions of ACA are quite good in solving some simple test performances. Wang [5] changed the dynamic vehicle routing problem into a series of static subproblems based on time axis. The path pheromone is reinitialized and the new problem with old unmet demand nodes and new nodes in next dispatch round is modeled to be solved.

Because emergency logistics is focusing on speed instead of quantity, the dynamic effective material management is to replace the inventory materials’ lack of mobility. The key point of emergency material dispatch in this paper is how to ensure the reliability of material supply stability in dynamic environment. Based on the characteristics of emergency material dispatch problem, the multiobjective mathematic model is established, an improved ACA optimization strategy and algorithm are developed, and a numeric example proved this method is feasible to generate real-time emergency material dispatch scheme corresponding to each point of the time axis.

2. Mathematic Model of Emergency Material Dispatch in Dynamic Environment

As to dynamic events of DVRP, the normal treatments are the event trigger mechanism and rolling horizon principle. While we adopt the event trigger mechanism, the new arrived material assignment is immediately inserted into the running dispatch scheme. Meanwhile we adopt the rolling horizon principle and divide the operation time into several small horizons. In the intervals, the material dispatch scheme is adjusted. The two treatments have both advantages and disadvantages. The former one may cause the frequent scheme altering, while the latter one may decrease the service quality.

2.1. Strategic Analysis

Because sudden incidents in the emergency logistics are various and discrete, the demand of materials is changing with the development of the situation. According to the actual material assignment operation, a rolling horizon multitask decision-making emergency material dispatch model is built based on time axis, and the dispatch scheme is updated while distributing unfinished dispatch scheme and addon dispatch scheme in these intervals.

In the interval, the event occurrences trigger the simple 𝐾 exchange neighborhood optimization and service denial strategy, while guaranteeing that every undamaged vehicle finally returns to the logistic center. Based on the instant statistics information about each node demands of the moment, the dynamic model can be converted to submodels. As the time goes on, the material demands of customer nodes will accumulate because of the frequent scheme altering or the assignment finishing rate as shown in Figure 1.

Figure 1 

Accumulated material demands of customer nodes.

The restart strategy based on rolling horizon principle has three advantages. Firstly, the dispatch scheme can change with the dynamic network attributes. Because the material flow is continuous, the new dispatch scheme round can be adjusted to take full use of limited transport network resources. Secondly, it can reduce the complexity of problem and generate simpler and easier schemes, in which every vehicle starts from the logistic center and the initial capacity is the loading capacity. Thirdly, the improved ACA algorithm can speed up response time by using the information generated from the solving process of former dispatch scheme.

2.2. Variable Parameter Descriptions

The vehicles that have been the starting must try their abilities to serve all nodes and the vehicles which are in the logistic center are dispatched according to actual need. The mathematic model of emergency material dispatch is built based on some cases as follows.(i)Emergency materials are all generic class material and dispatched as standardized container units whose unit of measurement is twenty-foot equivalent units (TEU).(ii)A route can only start from the logistic center and end in the logistic center.(iii)The type of vehicles is single and the lay time is ignored. The drive speed is constant and the initial route weight matrix 𝐴 is represented by the distance between nodes.(iv)The emergency materials should be dispatched into all customer nodes safely, cheaply, and fast during the transportation.(v)The parameters of the mathematic model are all integer decision variables.

According to the current research situation, the problem to be studied is described as a vector set of 𝐺=(𝑉,𝐷,𝑊,𝑅) at the moment of 𝑡. The problem is a network topology 𝑉={𝑣₀,𝑣₁,𝑣₂,…,𝑣_𝑛} that contains a logistic center and 𝑁 customer nodes (𝐶 node for short), where 𝑁 is the maximum number of the 𝐶 nodes. Each node has an associated location [𝑥_𝑖,𝑦_𝑖], and the network is in the Euclidean plane, so the travel distance 𝑑_𝑖𝑗 between any two nodes 𝑣_𝑖 and 𝑣_𝑗 is the straightline distance between them. The distance 𝑑_𝑖𝑗 between any two nodes of vector 𝑉 constitutes the route weight matrix 𝐷, which is symmetric and satisfy the triangle inequality: 𝑑_𝑖𝑗≤𝑑_𝑖𝑙+𝑑_𝑙𝑗. And the travel time can be obtained if you divide the distance by the drive speed, which is also a constant. Each 𝐶 node denotes a service request. All requests will randomly arrive within a planning fixed horizon. Accumulated material demands which contain both unfinished dispatch scheme and addon dispatch scheme are expressed as vector 𝑊={𝑤₁,𝑤₂,𝑤₃,…,𝑤_𝑛}, which corresponds to the emergency material demand of 𝐶 node 𝑣_𝑖 at the moment of 𝑡. Every vehicle must depart from and return to the logistic center, and the transportation parameter matrix is 𝑅=(𝐾,𝑄). The logistic center has 𝐾 service vehicles whose loading capacity is 𝑄 TEU. Definebinaryvariable:𝑥_𝑖𝑗𝑘=1ifvehicle𝑘drivesfromnode𝑣_𝑖tonode𝑣_𝑗0else. (1)

2.3. Objective Function

min𝑧₁=^𝐾𝑁_𝑘=1^𝑁_𝑗=0𝑖=0𝑑_𝑖𝑗𝑥_𝑖𝑗𝑘,min𝑧₂=^𝐾𝑁_𝑘=1𝑗=1𝑥_0𝑗𝑘. (2)(3)

The model seeks the optimization routes of the current snapshot so that the total dispatch service benefit is maximized. Equation (2) is an object function minimizing the total dispatch route length, which also makes the drive cost the lowest. Equation (3) is an object function minimizing the number of vehicles that have been put into use in the logistic center, which also cuts down greatly the fixed charges. And the two targets will guarantee the maximum of emergency material dispatch safety.

In this paper the multiple objectives are integrated into a single objective with the weight coefficient method, and the weight vector of the two objectives is 𝜔=(0.354,0.646), which was obtained from importance comparison matrix in AHP. Considering the data in numeric example, the value of 𝑧₁ ranged from 800 to 1000, and the value of 𝑧₂ ranged from 9 to 22. Therefore we used weighted summation of the two data to get the single objective function 𝑍 as follows: min𝑍=0.354⋅lg𝑧₁+0.646⋅𝑧₂. (4)

2.4. Constraint Conditions

Consider ^𝐾𝑁_𝑘=1𝑗=1𝑥_0𝑗𝑘≤𝐾;_{𝑁𝑁𝑗=0𝑖=0}𝑤_𝑗𝑥_𝑖𝑗𝑘≤𝑄,𝑘∈𝐾,_𝑁𝑗=1𝑥_0𝑗𝑘=^𝑁_𝑖=1𝑥_𝑖0𝑘≤1,𝑘∈𝐾,_{𝐾𝑁𝑘=1𝑖=0}𝑥_𝑖𝑗𝑘=1;𝑗={1,2,3,…,𝑁_{𝐾𝑁𝑘=1𝑗=0}},𝑥_𝑖𝑗𝑘=1;𝑖={1,2,3…,𝑁},𝑖≠𝑗. (5)(6)(7)(8)(9)(10)

Equation (5) ensures that the total number departed from the logistic center is less than the upper bound in every interval of dispatch scheme. Equation (6) ensures that the sum of total demand of all 𝐶 nodes in a vehicle’s route is less than its loading capacity. Equation (7) can ensure that all vehicles have departed from the logistic center and finally returned to the depot. Equations (8) and (9) ensure that every 𝐶 node has one and only one predecessor node and successor node and is visited only once by a vehicle. Expression (10) can ensure the travel route of every vehicle forbid loop node, especially looping in the logistic center.

3. Ant Colony Algorithm Design of DVRP

Dealing with DVRP, the new input information could make the former optimal dispatch scheme change into the suboptimal, even the infeasible scheme. In the static environment, it is tolerable to spend more processing time for commanders getting high quality or optimal dispatch scheme. However, commanders want to obtain the dispatch scheme immediately; the fewer seconds or minutes are the better in an emergency.

The main idea of this algorithm design is to improve the performance of basic ant colony algorithm by state transformation rules using load rating function and pheromone update rules based on ant ranking system [6]. Furthermore, the route pheromone information is reserved while solving the former dispatch scheme, so the latter dispatch scheme solving process is speeded up and the emergency material dispatch scheme is generated.

3.1. State Transformation Rules of ACA

Dorigo M put forward the adaptive pseudorandom probability transformation rule in ant-Q algorithm in 1993 [6]. In general approach of basic Dorigo ant algorithm, each of 𝑚 ants constructs several routes touring through all the given 𝑁 cities in every generation, which makes up the dispatch scheme. Starting at a random city an ant selects the next city using heuristic information as well as pheromone information, which serves as a form of memory by indicating which choices were good in the past. In order to enhance the accuracy and reality of the ants while they choose travel routes, this paper proposes the improved state transformation rule which takes full use of load rating in view of constraint expression (6) as follows.

Consider 𝜏𝑆=argmax^𝛼_𝑖𝑗(𝑡)⋅𝜂^𝛽_𝑖𝑗(𝑡)⋅𝜇^𝛾_𝑖𝑗(𝑡),if𝑅≤𝑅₀𝑝^𝑘_𝑖𝑗(𝑡)𝑝,else,_𝑘𝑖𝑗(𝑡)=𝜏^𝛼_𝑖𝑗(𝑡)⋅𝜂^𝛽_𝑖𝑗(𝑡)⋅𝜇^𝛾_𝑖𝑗(𝑡)_{𝑠∈allowed𝑘}𝜏^𝛼_𝑖𝑗(𝑡)⋅𝜂^𝛽_𝑖𝑗(𝑡)⋅𝜇^𝛾_𝑖𝑗(𝑡),if𝑠∈allowed_𝑘0,else. (11)(12)

The probability for ant 𝑘 to append arc (𝑣_𝑖,𝑣_𝑗) to its partial solution is then given as 𝑝^𝑘_𝑖𝑗(𝑡); 𝜏_𝑖𝑗 indicates the pheromone concentration information; 𝜂_𝑖𝑗 represents the expectation heuristic information, which equals the reciprocal of 𝑑_𝑖𝑗 and minimizes the total dispatch route length. 𝜇_𝑖𝑗 represents the load rating heuristic information, which minimizes the number of vehicles in use. 𝛼, 𝛽, and 𝛾 are constants that determine the relative influence corresponding to the pheromone values, distance values, and load rating values on the decision of the ant 𝑘. Random variable 𝑅 follows 0-1 uniform distribution. Parameter 𝑅₀(0≤𝑅₀≤1) determines state transformation rule, which also indicates the relative importance of exploitation versus exploration. The nodes set allowed_𝑘 contains the remaining nodes that have not been visited so far by ant 𝑘.

While ant 𝑘 has 𝑄_𝑖 TEU capacity at node 𝑣_𝑖, and the material demand of successor node 𝑣_𝑗 is 𝑤_𝑗 TEU. Function 𝜇_𝑖𝑗=𝑄_𝑖+𝑤_𝑗𝑄,if𝜇_𝑖𝑗>1,then𝜇_𝑖𝑗=0. (13)

3.2. Pheromone Update Rules

Once the 𝑚 ants of the colony have completed their computation, the better dispatch routes are used to globally modify the pheromone trail. In this way a “preferred route” is memorized in the pheromone trail matrix and future ants will use this information to generate new solutions in a neighborhood of this preferred route. For pheromone initialization we set 𝜏₀=1/(𝑁−1) for every arc (𝑣_𝑖,𝑣_𝑗) at the beginning of each iteration.

When ant 𝑘 moves from 𝑣_𝑖 to 𝑣_𝑗, a local updating is performed on the pheromone matrix, according to the following rule: 𝜏_𝑖𝑗(𝑡+1)=1−𝜌₁⋅𝜏_𝑖𝑗(𝑡)+𝜌₁𝜏₀, (14) where parameter 𝜌₁ determines the evaporation rate of pheromone.

The global pheromone matrix adapts the elite ant strategy, which only uses the optimal and suboptimal dispatch route. And it is updated as follows: 𝜏_𝑖𝑗(𝑡+1)=1−𝜌₂⋅𝜏_𝑖𝑗(𝑡)+𝜌₂Δ𝜏_𝑖𝑗(𝑡),Δ𝜏_𝑖𝑗(𝑡)=𝜙₁,if𝑒_𝑖𝑗𝜙∈optimaldispatchroute₂,if𝑒_𝑖𝑗∈suboptimaldispatchroute0,otherwise, (15) where parameter 𝜌₂ determines the evaporation rate of pheromone, 𝜙₁ and 𝜙₂ are constants, and 𝜙₁>𝜙₂>0.

3.3. Dynamic Modification Rules of Pheromone Matrix

During the emergency material dispatch, service requests are received in real time, and serviceability of requesting node is immediately verified. The target of this paper is to get the updated dispatch scheme based on rapid reaction to the dynamic information about the traffic capacity of transport network and demand nodes. However, the situations of dynamic environment can be divided into three types as follows:(i)the change in route weight matrix 𝐷,(ii)the change in material demand of customer nodes,(iii)the change in optimization objectives of material dispatch.

How to modify the pheromone matrix as much as possible to adjust the change in route weight matrix is the key point. Among these three situations, types (i) and (ii) are more common in connectivity and travel time reliability along with the development of the situation. The simplest strategy is reinitializing all pheromones to 𝜏₀ and then restarting the ant algorithm while making the solving process inefficiency, so the current solving process should make full use of the previous results of the dispatch scheme calculation in order to reduce the computation time. For description, the modification rules are as follows.

(1) Rule 1: Dynamic Distance Weight between Nodes. In the case of the distance between any two nodes 𝑣_𝑖 and 𝑣_𝑗 changing, 𝑀𝑎𝑥𝐷𝑖𝑠 is the maximum of all the 𝑑_𝑖𝑗 in the former route weight matrix 𝐷; the pheromone modification rule of arc (𝑣_𝑖,𝑣_𝑗) is expressed as follows: If𝑑_𝑎𝑖<𝑃×𝑀𝑎𝑥𝐷𝑖𝑠∨𝑑_𝑏𝑖∨𝑑<𝑃×𝑀𝑎𝑥𝐷𝑖𝑠_𝑎𝑗∨𝑑<𝑃×𝑀𝑎𝑥𝐷𝑖𝑠_𝑏𝑗<𝑃×𝑀𝑎𝑥𝐷𝑖𝑠Then𝜏_𝑖𝑗=𝜏₀⋅𝜏1+log^𝑖𝑗𝜏₀, (16) where parameter 𝑃∈[0,1] determines the adjustment range. 𝑃=0 means no adjustment, while 𝑃=1 means all pheromones are changing according to this rule.

(2) Rule 2: Dynamic Regulation of Customer Nodes Demands. In the case of the demand 𝑤_𝑖 of node 𝑣_𝑖 changing, 𝑁 is the number of all customer nodes; the pheromone modification rule of arc (𝑣_𝑖,𝑣_𝑗) is expressed as follows: 𝜏_𝑖𝑗=1−𝜃_𝑖⋅𝜏_𝑖𝑗+𝜃_𝑖⋅1𝑁−1 (17)

(1) 𝜂-strategy [2], as Figure 2. Consider 𝜃_𝑖=max0,1−𝜂𝜆⋅𝜂_𝑖𝑗;𝜆∈[0,+∞),1𝜂=()⋅𝑛⋅𝑛−1^𝑛_{𝑖=1𝑗=1,𝑗≠𝑖}𝜂_𝑖𝑗. (18)

Figure 2 

𝜂-strategy.

Because parameter 𝜂_𝑖𝑗 equals the reciprocal of 𝑑_𝑖𝑗, so it can be inferred that parameter 𝜃_𝑖 is inversely proportional to the distance between node 𝑣_𝑖 and the changed node.

(2) 𝜏-strategy [2], as shown in Figure 3. Consider 𝜃_𝑖=min1,𝜆⋅𝜓;𝜆∈[0,+∞),𝜓=max_{𝑥,𝑦∈𝑃𝑖𝑘}𝜏_𝑥𝑦𝜏_max. (19)

Figure 3 

𝜏-strategy.

From the equations it can be deduced that the closer arc (𝑣_𝑖,𝑣_𝑗) is next to the changed node, the more equalized pheromone is in the route.

To illustrate the 𝜂-strategy and 𝜏-strategy operation, we used a 10×10 grid of cities and set 𝜆=1 to visualize how the distribution of reset values takes place. In Figure 2 the distribution for the 𝜂-strategy is proportionate to Euclidean distance, while Figure 3 shows that the 𝜏-strategy tends to distribute along the path. By comparing the two strategies, we choose 𝜂-strategy to deal with dynamic regulation of customer nodes demands.

(3) Rule 3: Dynamic Regulation of Nodes Number. When the increased node number equals the decreased number, replace the nodes and then deal in Rule (1).(1)When the increased node number is more than the decreased number, directly delete the row and column in which the node stays.(2)When the increased node number is less than the decreased number, insert the new node, calculate the distance weight, and initialize pheromone to 𝜏₀.

3.4. Detailed Steps of Improved ACA

Step 1. Calculate the dispatch scheme at the initial time of 𝑡₀ using improved ACA and generate the solved result.

Step 2. Overwrite the pheromone matrix into the excel document “pheromone” until the iteration stops.

Step 3. When calculating the new dispatch scheme, read the data in excel document “pheromone” and save it initialize pheromone matrix.

Step 4. According to the specific changing information of dynamic environment, update the pheromone matrix with modification rules in Section 3.2.

Step 5. Calculate the dispatch scheme at the time of 𝑡 using improved ACA, generate the new solved result, and then go to and execute Step 2.

4. Experimental Results

The emergency material dispatch system contains one logistic center and 100 demand nodes. The loading capacity of each vehicle is 200TEU. The static data of numeric example in this paper adapts the C101 and R101 from Solomon benchmark at the initial moment. As the situation develops, such as earthquake, snowstorm, and war, each problem has 5 groups of dynamic data which are used to depict 5 different dynamic degrees of 10%, 30%, 50%, 70%, and 90%, respectively. We randomly choose the proportional elements of distance matrix 𝐷 and nodes demand vector 𝑊 from 𝐺=(𝑉,𝐷,𝑊,𝑅) to simulate the dynamic environment. The distance matrix is not symmetric and does not satisfy the triangle inequality any more. The test data set is designed by modifying benchmark method, which was proposed by Lackner [7] in 2004.

The test work is completed on a desktop computer with Intel Pentium Dual-Core [email protected] GHz, 2046 MB RAM, and Windows XP SP3 Professional Operating system. The algorithm is coded and compiled using M file in Matlab 7.8.0. The simulation experimental results are shown in Table 1.

Table 1 

Calculation results of Solomon benchmark problem for various dynamic degrees.

For description convenience, we call the column (1) C1, column (2) C2, and so on. C1 represents the data type of Solomon’s 100-customer benchmark problems. C2 and C3 are the best known results obtained by the heuristics for Solomon’s static VRPTW problems reported in http://www.sintef.no/static/am/opti/projects/top/vrp/bknown.html (last updates: March 24, 2005). C4 DoD is short for the degree of dynamism, which means the changing nodes’ percentage of all the nodes. Among the data items from the table, C5, C6, and C7 are the computational results for proposed approach, and they are also the average number of used vehicles (NV), the average travel distance (length), and the average CPU time (CPU, in s), respectively. The CPU time is the time used for solving each static version problem. So are the C8, C9, and C10.

Other parameters are set as follows. The ant colony size is 50 ants, the number of iterations is 100, and the heuristic factor is g 𝛼=1,𝛽=5,𝛾=1. The evaporation rate of pheromone 𝜌=0.4. The state transformation rule parameter 𝑅₀=0.6, and the dynamic modification rule parameters 𝑃=0.1, 𝜆=2.

As shown in Table 1, the basic ACO algorithm with a simple restart strategy uses one more vehicle and generates the suboptimal total dispatch route, while its solution time is far more than other algorithms. The improved ACO algorithm with pheromone reserved strategy can meet the requirement of quick decisions, which saves an average of 50% computation time searching the optimal emergency material dispatch scheme. And the total computation time is less than 80 seconds in absolute amount dealing with 100-node problem. Because of simplified modeling, this improved ACO algorithm is superior to the LNS algorithm in [8]. So the improved ACO algorithm designed in this paper is optimal in the vehicles number and suboptimal in the total route length.

5. Conclusions

The multiobjective materials dispatch model is built by linear weighted method. The dynamic information generated in the process of emergency material dispatch is disposed in batches by rolling horizon strategy, such as traffic capacity of transport network or demand nodes. Eventually, a numeric example is designed in 5 different dynamic degrees to prove the validity of this model and algorithm. The modified ant colony optimization algorithm utilizes the pheromone reserved by old material dispatch scheme to initialize the ant colony pheromone matrix in seeking for the next new dispatch scheme. This method improves the quality of scheme solutions and obviously the solving efficiency, while generating material dispatch scheme in real time based on current emergency material dispatch network. However, there are still some problems that need to be studied such as split demand, complete set, and multimodal transport of complex dynamic emergency materials dispatch problem.