DISASTER RELIEF ROUTING IN LIMITED CAPACITY ROAD NETWORKS WITH HETEROGENEOUS FLOWS

. In the aftermath of a major earthquake, delivery of essential ser-vices to survivors is of utmost importance and in urban areas it is conducted using road networks that are already stressed by road damages, other urban traﬃc and evacuation. Relief distribution eﬀorts should be planned carefully in order to create minimal additional traﬃc congestion. We propose a dynamic relief distribution model where relief trucks share limited capacity road net- works with counterﬂows resulting from car traﬃc. We develop a MIP model for this problem and solve it by decomposing the road network geographically and solving each subnetwork iteratively using the Relax and Fix method.


1.
Introduction. Traffic congestion occurs in the aftermath of earthquakes and other disasters due to demolition debris, relief distribution, evacuation traffic and other urban traffic. The problem of post-disaster relief distribution over limited capacity road networks has not received much attention in the Disaster Operations Management (DOM) literature [13,18,4,33,23]. To our knowledge, the only citation that mentions post-earthquake traffic management is that of Feng and Wen [21] who emphasize traffic control given certain origin-destination (o-d ) relief trips. However, in this article the vehicle routing problem is not explicitly considered, rather a path selection model for satisfying traffic flow between o-d pairs and other incoming traffic is presented.
In this study, we consider optimizing the routes of relief trucks that provide service to handicapped, old, sickly or otherwise affected people who cannot walk or drive to the nearest aid distribution center. Each street has a given number of needy people determined by census data. Once a truck arrives, these people fetch their goods within a few minutes and return to their homes. The truck then moves on to the next street on its tour. A truck is also allowed to transit a street without stopping. In our model, a set of relief trucks enters a neighborhood and visits a number of streets, i.e., links in a road network, during a given time span (maximum tour length) while sharing the road capacity with car traffic flows. We assume that relief distribution to needy persons is repeated daily until consumer supply chains are restored. Therefore, relief traffic is inevitably combined with car traffic in urban areas. Congestion becomes critical especially in old city neighborhoods where single lane alley clusters exist. The objective function of the proposed model maximizes the number of persons served over all visited streets as well as the number of people leaving the neighborhood in private cars over the planning horizon.
This problem resembles the "arc routing problem with profits" mentioned by Archetti et al. [7] where capacitated vehicles collect profit from arcs while serving arc demand. The arc routing problem is also considered by Aroz et al. [5,6] who consider the rural postman problem with the objective of maximizing total profit minus total travel cost, and by Feillet et al. [20] where uncapacitated vehicles can collect profit from the same arc multiple times. None of the so-called arc routing articles consider "arc capacity", or, in other words, restricted traffic flow. In the disaster operations management context, the arc routing problem is recently discussed in debris clearance [2] where a typical uncapacitated vehicle routing formulation is adopted to clear blocked roads with the goal of minimizing total clearing time. The authors show that only small networks can be solved to optimality with this formulation, larger instances are solved using a local search procedure that might fail to find feasible solutions in larger scale networks. Note that, in these arc routing formulations, there are no counterflows that represent another type of traffic flow (e.g., private car traffic in our study).
Our problem, and, therefore, the proposed mathematical model, differ from the above-mentioned arc routing problems and other relief distribution models in the following respects: a) There exists counterflow traffic in the road system. Outbound private cars traveling towards highways/boulevards share narrow road networks with inbound relief trucks; b) Traffic flow over each arc during a given time is limited by arc flow capacity; c) A relief truck blocks all other traffic while serving a street; d) The objective function maximizes both the profit collected (number of people served) over all visited edges and the number of people leaving the neighborhood in outbound private cars. In other words, the goal here is integrated, the obstruction caused by the relief truck traffic is minimized while maximizing the amount of people served. To our knowledge, the relief distribution problem has not been investigated in an environment where traffic counterflows exist. This results in a novel mathematical programming model with two different types of interacting flow variables. Thereby, the model requires a different solution approach that takes advantage of the particular model structure.
To solve this problem, we develop a dynamic arc routing model with relief truck and other kinds of traffic flows over limited capacity road networks. The outputs of the model consist of a set of vehicle routes and their street service schedules. Since the arc routing model is NP-Hard, we decompose the road network into smaller scale subnetworks and solve resulting smaller scale MIP models by the iterative Relax and Fix (RF) method. RF method generically reduces the number of integer variables in a model by discretizing only a subset of the variables. In our implementation, we demonstrate that selecting the vehicle flow variables related to a few consecutive time intervals and relaxing the rest produces an effective solution strategy. In each iteration, the optimal values of discrete variables obtained in the previous iteration are fixed. Then, a new subset of discrete variables is created and the problem is solved again.
In the literature, RF heuristics have been successfully used in different applications, such as large-scale lot-sizing problems [10,22,9], facility location and stochastic supply chain problems [19,3], among others. Therefore, the implementation of RF to the vehicle routing problem introduced here is a novel contribution. The motivation for using an RF heuristic here is to take advantage of the problem structure that involves route construction for all relief trucks over multiple periods. In fact, using the RF method, relief truck routes are constructed gradually adding one or more new links to a truck's partial route in each iteration. Hence, the selected solution method agrees with the problem structure and we demonstrate that the RF heuristic reduces CPU times.
In our numerical experiments, we test RF results against optimal solutions for a number of smaller scale road networks. We observe that the solutions generated by RF are as good as the optimal ones and their CPU times are significantly lower. However, experiments with larger scale problems indicate that RF is not capable of solving them within reasonable time. In order to solve large scale problems, we combine smaller scale disjoint adjacent networks into aggregate networks. Each small network is then solved with RF and the solutions of disjoint networks are combined to represent the solution of the aggregate network. We also attempt to solve each aggregate network optimally and compare the objective function with the sum of the objective functions of its disjoint subsets. This comparison illustrates the degree of suboptimality resulting from dissecting large regions into smaller disjoint regions and solving the latter with RF. Our numerical results illustrate that there is a limit to the size of the networks that can be solved optimally. Therefore, the approach of using RF on a set of smaller scale networks appears to be a feasible and computationally efficient way to deal with this problem.
2. Literature survey. The DOM literature usually treats the traffic flow problem separately from the relief distribution problem. To deal with the traffic problem, Daganzo [17] presents the Cell Transmission Model (CTM) that partitions the car traffic network into grids and controls traffic flow advancement over adjacent grids (cells). Liu et al. [28] use combined optimization-simulation method to solve a bi-level CTM that maximizes traffic flow and minimizes travel time. The maximum flow model is another model category that authors utilize to represent outgoing (evacuation) traffic flows in DOM. The quickest flow model given by Miller-Hooks and Patterson [29] minimizes network clearance time in a capacitated evacuation network. Road intersections and lane changes are considered by Cova and Johnson [15]. Lim et al. [26] propose the capacitated maximum flow model and Campos et al. [12] minimize intersection delays and traffic congestion. Among traffic flow management articles, Cui et al. [16] are the only authors so far to consider evacuation flows and the counterflows of rescuers over a capacitated road network, but they do not consider relief distribution or routing.
An area where the closure of roadways is explicitly considered in the relief distribution problem is the road repair problem [36,14,25]. Here, some links in the road network are closed and relief trucks can traverse only repaired or healthy links. In these models the road network is usually uncapacitated, that is, the only restriction on truck traversal is binary and traffic flow restrictions on arcs are otherwise absent. These models aim to schedule link repair jobs such that relief can be distributed. The latter models are different from the one proposed here both in context and problem structure.
Relief distribution models without traffic flow constraints are abundant in the DOM literature. All related vehicle routing models in the literature are based on the "node routing problem" where demand exists on nodes rather than arcs. These models can be classified according to their type of routing model structure [33]. Examples of dynamic network flow routing models are provided by Haghani and Oh [24] who minimize operation costs, by Ozdamar et al. [32] who minimize cumulative unmet demand, and by Afshar and Haghani [1], Ozdamar and Demir [31] and Najafi et al. [30] all of whom consider the same goal. Route enumeration model examples are provided by Balk et al. [8] where the goal is to minimize transport, backorder, and lost sales costs, and, by Lin et al. [27] who propose a multi objective model with the goals of minimizing cumulative unmet demand, minimizing travel time, and minimizing the pairwise difference between demand fill rates at nodes. Classical vehicle routing models are proposed by Vitoriano et al. [34,35] who propose a multi-objective model where the goals are minimizing transport and operating costs, minimizing the maximum looting probability, and minimizing the maximum percentage of unmet demand over nodes. Tour length limitation is considered by Berkoune et al. [11] who minimize travel distance. Solution methods in these articles are usually exact methods, and some propose heuristics such as Lagrangean relaxation [24,32], Genetic Algorithms [27,11], and, hierarchical network disaggregation and optimization [31]. Except for the article by Ozdamar and Demir [31], the size of the problems solved are quite modest.
The literature survey summarized above indicates that the model proposed here has novel properties such that it considers an arc routing problem rather than a node routing problem over a network with two types of interacting flows (truck traffic and private car traffic) sharing the same road capacity. Furthermore, the selected RF solution method with network decomposition has not been applied previously to the routing problems in the literature.
3. The mathematical model. The mathematical model (Model R) for the problem at hand is a MIP. The model assumptions are as follows: a) There is a planning horizon consisting of T time periods. b) Residents in each street are served with relief materials by at most one truck during T . The number of residents in need of relief materials on street (i, j) is denoted by pop ij . c) If a truck is on a street in any time period, it can either be traversing it in transit or serving it. d) When a truck is on a street, it blocks the street to all other traffic. e) Each street junction j generates an expected amount of outbound car traffic Q jt in each time period t. Q jt arise from outbound traffic on streets that are incident to junction j. f) A street (i, j) has a traversal capacity a ij (number of cars per period). g) Trucks and cars traverse a street with limited speed. The notation used in the model is provided below. x ijt : number of cars traversing link (i, j) in period t.
The model R is provided below. Objective Function The objective function maximizes the time weighted number of served residents on each street and the number of cars reaching the exit nodes of the neighborhood (region) connected to major highways. The time weights favor early servicing of trucks and early exit of cars to highways. Constraints Constraints (1) state that a link (i, j) can only be serviced once by one truck during the planning horizon.
Constraints (2) ensure truck flow balance over each road junction in each period.
Constraints (3) state that a truck that enters the neighborhood from any entry/exit node in E has to return to any entry/exit node in E by the end of the planning horizon. Note that trucks traverse links (i, j) ∈ D only in transit.
Constraints (4) and (5) allow each truck at most one entry and one exit to/from the neighborhood.
Constraints (6) do not allow transit trucks to traverse a link that is being served at a given time period. x Constraints (7) impose link capacities per period on car traffic in transit. These constraints ensure that car traffic is blocked while a link is in service or while it is being traversed by a truck. t q=1 i∈N Constraints (8) enable the balance of car traffic flow over each road junction j in each period.
The model R presented here describes an arc routing problem with two categories of traffic flows and various flow restrictions imposed on bidirectional links. 4. Solution methodology. Since routing problems are known to be NP-Hard, we choose to use an iterative MIP heuristic that solves smaller and more tractable MIP problems in each iteration. This Relax and Fix (RF) heuristic is based on partitioning a set of discrete variables into P disjoint sets according to a predefined criterion. In our implementation, we partition vehicle related discrete variables y and v in Model R into two sets. In each iteration, only those variables pertaining to a small time interval ∆t are defined as discrete variables (∆t < T ) and the remaining y and v variables outside this interval (t / ∈ ∆t) are relaxed. Therefore, a subproblem of much lower complexity is solved in each iteration. Once a subproblem is solved, the discrete variables defined in the incumbent iteration are fixed at their current optimal values and the process is repeated by re-solving the problem for all the remaining sets. A pseudo code for RF is given in Algorithm 1. Line 1.4 indicates that only the variables that belong to the corresponding partition (∆t) are requested to be discrete.
The iterations of the RF approach are applied in a forward scheme performed from the first ∆t periods where truck related y and v variables are defined as k + +

5:
Solve Model R

6:
Fix y ijcq and v ijcq for q = (k − 1)∆t + 1, ..., k∆t at their optimal integer values 7: end while 8: Report final solution 9: End (Algorithm1) integers only for ∆t periods and the rest of the variables are defined as positive float variables. Once an integer solution is found, these variables are fixed at their optimal values as either 0 or 1, and the next set of variables ∆t + 1, ∆t + 2, ..., 2∆t are defined as integers while the remaining variables are still float. As the length of ∆t increases, the number of iterations is reduced, because there are T /∆t iterations to be completed. However, the computation time of each iteration may increase. Hence, it is important to reach a balance between the number of iterations and the length of ∆t in order to reduce the total computation time. 5. Numerical experiments. We performed our computational experiment on a portion of the road network of Fatih County in Istanbul depicted in Figure 1. Every intersection in this road network is represented by a node and every road segment between two intersections corresponds to a link in our problem. The resulting network depicted in Figure 1 contains 260 nodes and 390 links in total. The numbers on the figure indicate the nodes in the network. The nodes marked with a circle reside on wide boulevards and they are linked to the dummy node set E (dummy exit nodes are excluded from Figure 1).
Due to the large problem size, Model R cannot be utilized to obtain a solution for the entire Fatih network within reasonable time. To remedy this, we have divided this large network into seven smaller subnetworks or regions such that each region has a distinct set of links (however two regions may possibly have common nodes). The boundaries of these regions were determined such that all regions have similar numbers of links (between 46-64 non-dummy links). The boundaries of the seven regions are marked on Figure 1 with a dashed line. Each region is allocated a fleet of 3 relief trucks, resulting in 21 trucks for the entire network.
All of the runs in this section were performed using the MIP solver of CPLEX version 12.1 on a Intel Pentium Core 2 Duo 2.4 GHz computer. For each run, we used a planning horizon of 200 minutes, which was discretized using 5 minute intervals, resulting in a total of T = 40 time periods. The population to be served on each street (pop ij ) was calculated by allocating the total population of this portion of the Fatih County to each street in proportion to its length. Truck service time for each street is proportionate to its pop ij value. The total number of cars originating from a road junction over the time span T was calculated by summing the car owners on the links incident to that road junction and taking about 10 to 20% of the sum. Then, this number was distributed probabilistically over T using the discrete uniform distribution to generate Q jt . On the average, each car carries 3 persons and car ownership ratio to population is 1:3 in Istanbul. Arc flow capacities for car traffic were set by the following formula: (rush hour velocity per minute * number of lanes on street * time interval length) / (average car length). Transit times for both cars and trucks take 1 time interval. In order to explore the performance of the RF approach described in Algorithm 1, we ran the algorithm on all seven regions for a range of ∆t values (∆t=2, 4,8) and compared these results with the solution of the Model R in (1)-(8) as a single MIP problem without using the RF approach. To provide a fair comparison, we allocated a total of 3600 seconds of CPU time and 3600/ 40/∆t seconds for each iteration of the RF algorithm. (Here * indicates the smallest integer greater than or equal to *.) This ensures that the time spent to run the RF algorithm will also not exceed 3600s in total. Table 1 summarizes the solution quality and computational performance of the RF approach compared to the optimal or near-optimal solution of the single MIP (indicated by "No RF" rows in column 2). In this table, the first column displays the region of the problem solved (as indicated in Figure 1) and the number of links in that region. The second column indicates the solution method used while the third one displays the objective function value obtained for this method. Columns

Region Method
Obj. CPU(secs.) Rel. 4 and 5 display the total CPU secs. taken to solve the problem and the relative gap at the end of the CPLEX run for rows that correspond to a single MIP run ("No RF" runs). Finally, the last column (indicated as % Dev in Table 1) gives the percentage deviation of the RF solution from the solution obtained by the single MIP run displayed in the rows titled "No RF". Here, a positive value corresponds to a superior solution by the RF approach. In terms of the solution quality, the results in Table 1 indicate that the RF approach provides high quality results that are within around 1% of the single MIP (No RF) objective function value. Note that the single MIP runs often terminate with a high relative gap value at the end of the allowed 3600 CPU secs. Therefore, this approach does not guarantee an optimal solution even for the smaller subnetworks formed by decomposing the original Fatih network. In fact, RF approach is capable of producing slightly better results for three out of the seven regions. While the performance of the approaches with and without RF is comparable in terms of the objective function value, the difference in CPU times is very significant. Except for Region 1 where the CPU time of the run without RF is slightly less than that of the version with RF only for ∆t=2, all runs with RF require much less CPU time than those without RF. Moreover, the computational time requirement of the RF approach decreases as ∆t increases while also improving the solution quality. The latter result is expected because a larger ∆t implies a longer look-ahead horizon for RF and a lower number of CPLEX iterations. The results indicate that ∆t=8 is the best setting in terms of both CPU time and solution quality. For ∆t=8, the RF algorithm requires on the average 19.96% of the single MIP run's CPU time and yields comparable results. In fact, the average percentage deviation of RF over the single MIP run is only 0.11%. Overall, the computational results in Table 1 demonstrate that the RF algorithm is advantageous in producing high quality results in limited time.
In addition to solving each of the seven regions separately, we also experimented with combining some of the regions that have common borders to build larger networks that are solved as a single MIP. The objective function values of the larger regions were then compared with the sum of the individual smaller region objective function values in order to observe the degree of suboptimality caused by solving the smaller regions instead. We also compared these solutions to the results of the single region runs with RF from Table 1 to evaluate the performance of the RF approach for larger size networks. To allow a fair comparison, we again limited the CPU time for each run to 3600 seconds per region, i.e. for a single region problem, the CPU limit is 3600 CPU secs., whereas for problems with two, three and five regions, the limit is set at respectively 7200, 10800, and 18000 CPU secs. Similarly, the number of trucks allocated to each problem is 3 times the number of regions in that problem.
In Table 2, we summarize the computational results of Model R for each region as well as various combinations of regions. The first column in this table describes the regions included in the run and how the run was performed. The regions written consecutively (without a plus sign in between e.g. 12) indicate that those regions (regions 1 and 2) are combined (merged) into one single MIP problem. On the other hand, if regions are listed with a "+" sign in between (e.g. 1+2), a separate problem is solved for each region in the list and the summation of the statistics from all involved regions are displayed in that row. For instance, the second column in Table 2 titled "Obj." for subnetwork "34+5" displays the summation of the objective function values of two problems: one combined problem involving merged regions 3 and 4 and another problem with only region 5. The abbreviation RF next to the subnetwork indicators indicates that each region in the list is solved using the RF approach with ∆t=8, which provides the best performance among all three ∆t values. Rows that do not include an RF designation indicate that the listed problems are solved as a single MIP whether combined or single. In Table  2, the 3 rd and 4 th columns display the total CPU time used by CPLEX to solve all problems listed in the first column and the relative gap at the end of this time for rows corresponding to a single run in merged subnetworks. Finally, the total number of links involved in the merged problems is provided in the last column.
The results in Table 2 demonstrate the impact of decomposition with or without RF on solution quality and computational performance. The runs involving the same set of regions are grouped together and shaded in the same color. The results indicate that solving Model R directly using CPLEX is not a viable option as problem size (as reflected by the number of links in the network) increases within the imposed CPU time limits. While smaller problems with up to around 120 links can be solved within a reasonable relative gap, larger problems yield large relative gaps and significantly inferior solutions. For instance, subnetworks "123", "12345", "34567" could not be solved as a single MIP problem, that is, truck routes could not be constructed, only the flow of outbound private cars was enabled. Based on these results, we can conclude that decomposition is crucial in obtaining good quality solutions within reasonable time. For each set of regions, we report results of decomposition results both with and without RF. Comparison of these results suggests that runs that employ RF yield a solution quality that is comparable to their counterparts without RF. Though the "solved" single MIP solutions for the combined regions outperform the decomposed RF solutions, we observe that the decomposed RF solutions require an average of only 9.29% of the computational time required by the combined (first row in each group) subnetwork without RF. Based on these results, we can conclude that using decomposition along with RF offers an efficient alternative for the solution of larger size problems. 6. Conclusion. In this study, we focus on an arc routing problem that arises in disaster relief delivery operations. The model proposed here is novel because it acknowledges the existence of counterflow traffic during relief distribution. The resulting model describes an arc routing problem on a capacitated road network with two types of traffic flows: relief trucks and other traffic such as private cars. The proposed model strikes a balance between delivering relief to the highest number of affected residents and allowing the most private car traffic flow over the same road network. To our knowledge, this capacitated arc routing model with multiple types of interacting vehicle flows has not been studied in the vehicle routing literature. Furthermore, it is also different than the relief distribution problems in disaster management.
Given the complexity of the resulting MIP model, we propose an iterative RF algorithm that relaxes some of the discrete variables and solves a MIP problem with reduced number of discrete variables in each iteration. While the RF approach has been applied to lot-sizing, facility location, stochastic supply chain and other complex MIP problems, to our knowledge it has not been used to solve vehicle routing problems. Over all, we can say that this study contributes to the literature by both defining a realistic arc routing problem in disaster management and proposing a novel approach for its solution.
Computational results on a realistic size network with 390 links and 260 nodes indicate that solving a problem of realistic size without decomposition is not feasible and a network decomposition approach is inevitable. The proposed approach, where the network is decomposed into a number of smaller subnetworks each of which is solved using the RF algorithm produces encouraging results. As we compare the performance of the decomposition approach with and without RF, we also see that the decomposition approach that employs the RF algorithm is capable of solving subproblems within significantly reduced CPU times without a significant loss in solution quality. These results suggest that RF is a promising approach in solving dynamic vehicle routing problems that can be decomposed over the time dimension.