Bi-objective decision support system for routing and scheduling of hazardous materials

Socio-Economic Planning Sciences xxx (2014) 1e14

Contents lists available at ScienceDirect

Socio-Economic Planning Sciences journal homepage: www.elsevier.com/locate/seps

Bi-objective decision support system for routing and scheduling of hazardous materials Rojee Pradhananga a, *, Eiichi Taniguchi b,1, Tadashi Yamada b, 2, Ali Gul Qureshi b, 3 a b

Department of Mechanical and Industrial Engineering (MIE), College of Engineering, Qatar University, P.O. Box: 2713, Doha, Qatar Graduate School of Engineering, Kyoto University, C-1 Kyotodaigaku Katsura, Nishikyo-ku, Kyoto 615-8540, Japan

a r t i c l e i n f o

a b s t r a c t

Article history: Available online xxx

This paper presents a Pareto-based bi-objective optimization of hazardous materials vehicle routing and scheduling problem with time windows and shows its application to a realistic hazardous material logistics instance. A meta-heuristic solution algorithm is also proposed, which returns a set of routing solutions that approximate the frontier of the Pareto optimal solutions based on total scheduled travel time and total risk of whole transportation process. It works in a single-step fashion simultaneously constructing the vehicle route and selecting the optimal paths connecting the routed locations from a set of non-dominated paths obtained in terms of travel time and risk value.
2014 Elsevier Ltd. All rights reserved.

Keywords: Hazardous material Vehicle routing and scheduling problem with time windows Multi-objective optimization Realistic logistics instance

1. Introduction The nature of Hazardous Material (HazMat) is such that its production, storage, and transportation activity inherit many risks to both society and environment. Of the various modes of the HazMat transportation, large proportion of the HazMat shipment is carried out in-land by trucks. About 94% of total HazMat shipments in US are reported to be carried out by means of truck [17]. Despite of the efforts to mitigate the adverse effects of HazMat, accidents do happen [18] and the consequences in most cases are severe. Amount of HazMat transported by all transport modes is increasing year after year [24]. Without doubt, the HazMat incidents in highways are the most tremendous in number than in any other modes. Of the 15,007 HazMat incidents that occurred in US in 2011, 12,795 occurred in-land in highways causing an estimated damage of 64 million USD [43]. In regard to this fact, land based transportation of the HazMat is the main concern in this paper. At beginning, the line haul transport of HazMats is basically a shortest path problem involving large shipments between two defined points, for example from a refinery to the central depot. Once it reaches to the depot, it must be distributed to the end users.

* Corresponding author. Tel.: þ974 55958476. E-mail addresses: [email protected] (R. Pradhananga), [email protected]. kyoto-u.ac.jp (E. Taniguchi), [email protected] (T. Yamada), [email protected] (A.G. Qureshi). 1 Tel.: þ81 75 383 3229; fax: þ81 7 5950 3800. 2 Tel.: þ81 75 383 3230; fax: þ81 75 950 3800. 3 Tel.: þ81 75 383 3422; fax: þ81 75 950 3800.

The problem therefore is to determine efficient routes for a fleet of vehicles transporting the HazMats from depot to a set of customers. A similar problem to the latter case is commonly known as Vehicle Routing and scheduling Problem with Time Windows (VRPTW) [12,42]. The VRPTW is a complex problem than the shortest path problem. In HazMat transportation context, the VRPTW is extend to the Hazardous materials Vehicle Routing and scheduling Problem with Time Windows (HVRPTW) with additional objective of minimizing risk in addition to the conventional objective of minimizing the cost. Such multi-objective nature further increases the complexity of the HVRPTW. Therefore, despite of the practical significance, the HVRPTW-related research is very scant. Present study is an effort to fill in the gap of the literature and aims for development of an efficient optimization model for the HVRPTW addressing the practical aspects of HazMat transportation: a Presenting an optimization model that explicitly considers the multi-objective characteristics of the HazMat distribution problem. b Demonstrating the model’s application to a realistic Intelligent Transportation System (ITS) data-based HazMat logistics instance. Travel time data used for the HazMat logistics instance is the average of the historical data accumulated through Vehicle Information and Communication System (VICS), an ITS application in Japan. The HVRPTW involves a wide variety of stakeholders with often conflicting objectives. While the primary objective for the carrier is to minimize the transportation cost by optimizing the vehicle

http://dx.doi.org/10.1016/j.seps.2014.02.003 0038-0121/
2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

2

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

Fig. 1. The hypothetical urban road network and its transformed HVRPTW instance.

number (fixed cost of vehicle) and the total scheduled travel time (operation cost of vehicle), that of the local government and the people in the area is to reduce the risk. It is quite difficult to exactly determine how important the particular objective to the particular group is. Therefore, multi-objective analysis is an important tool in HazMat planning. It serves stakeholders of various needs for decision making in the HazMat transportation [17,25] and more importantly, it provides a basis for evaluation of various city planning policies. Although there exist some literature on multiobjective HVRPTW [27,46,47], the problem was solved using weight-based scalar approach. Scalar approaches are easier but major drawback with these approaches is that the optimal routing is highly influenced by the assigned weight values. Therefore, it requires precise determination of weight values which needs extreme analysis of the field data. More importantly, the approach results a single optimal solution and to examine trade-off among the objectives, the problem must be solved several times, which on the whole takes longer computation time. This paper takes an initiative to solve HVRPTW applying concept of Pareto optimization. Pareto optimization works simultaneously to determine a set of non-dominated solutions applying conditions of Pareto dominance, and is the best approach to deal with multi-objective problems. Typically, path choice and routing are two separate steps in any classical VRPTW. During the path choice, shortest paths from each customer to all other customers including the depot are determined. The pre-determined paths are then used in routing step to determine the optimal order of customers to be visited by the fleet of vehicles satisfying the demand and time windows constraints. Risk and time (cost) are the two indispensable objective functions of the HVRPTW that governs both path choice and routing. Scalar approach such as weighted sum approach combines the multiple objectives of path choice into a single objective and result a single

path between the customers. Therefore, it allows the two steps to be processed separately, similar to the case in the classical VRPTW. Pareto optimization, in HVRPTW results in a set of non-dominated paths for path choice between each customer pair. Therefore, the final non-dominated set of routing solutions must be based on all these sets of shortest paths. Consequently, it becomes necessary to process the two steps together as a single-step process. Fig. 1 along with Table 1 illustrates the necessity of single-step procedure in Pareto-based optimization of HVRPTW. Route optimization is shown comparing the classical two-step approach with the proposed single-step approach on a small hypothetical HVRPTW instance derived from an urban road network with 4 nodes and 10 links. Nodes b and d are customers 1 and 2 in the network and have demand values of 200 units each. Node a is a depot (customer 0) and node 3 is a non-customer node. Capacity of vehicle is 1000 units. Time windows at the depot and the customer nodes are given in the boxes. The first and the second values in the parenthesis represent the travel time and the risk value associated with the links, respectively. These values are same for the links in both directions connecting same pair of nodes. The travel time and the time windows are in the same units. Unloading time of 10 in the

Table 1 Routing results in single-step and two-step optimization. Optimal routing

Customer order

Single-step optimization 1 0/2/1/0 2 0/2/1/0 3 0/2/1/0 Two-step optimization (Step 1 0/2/1/0

Arc order

p4 p3 p3 2) p4

Total scheduled travel time

Total risk

/ p10 / p6 / p10 / p6 / p10 / p7

120 110 100

130 150 180

/ p10 / p6

120

130

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

same units as that of travel time has been used. For simplicity, the available number of vehicles is assumed to be 1. The urban road network has been transformed into a directed graph of vertices (customers and depot) and arcs (paths) as shown in Fig. 1. Choice of path for movement between the vertices is independent of time windows and unloading process. Hence, minimizing total travel time and total risk of all the associated links are the objectives used, to determine the arcs. Based on these objectives, ten arcs (non-dominated paths) were obtained for path choice in the instance. The single-step approach uses all these arcs to determine the optimal non-dominated set of routing solutions. On the contrary, the two-step approach requires a single arc to be identified from each vertex to another vertex in the first step. The two-step approach considered in this illustration utilizes a scalar method (weighted sum method) to select a single arc among the set of non-dominated paths for path choice between a pair of vertices. The weight values are taken equal to 1 for both objectives. Therefore, as shown in the figure, the final non-dominated routing solutions in the two-step approach are based on only six arcs among the ten arcs, used in single-step approach. Time value in routing stage depends on both time windows constraint and unloading process. Therefore, total scheduled travel time, which is sum of the waiting time, unloading time and travel time of the arcs in use is defined specifically for routing. Minimizing the total risk and the total scheduled travel time are the objectives considered for optimal routing. Table 1 shows the final routing solutions obtained in the single-step and two-step optimizations. Neglecting the feasible routing solutions with avoidable waiting times, there are three non-dominated routing solutions obtained with single-step optimization. While due to time windows restriction, there is only one final optimal routing solution for twostep optimization. The total scheduled travel times given in Table 1 for both cases of optimizations are obtained with zero waiting time values. These zero values are obtained for this particular illustration by either re-ordering of appropriate customer nodes or by delaying start of the vehicle at depot, yet ensuring the time windows constraints are satisfied. While the first solution in the single-step optimization is the same as that obtained in the two-step optimization, two more Pareto optimal solutions are obtained with the single-step optimization. Therefore, it can be concluded that all the possible arcs (non-dominated paths for path choice) must be used to obtain full range of Pareto optimal solutions in the HVRPTW. This ultimately illustrates the need of single-step optimization to investigate use of several nondominated paths for obtaining several non-dominated routing solutions. Therefore, a new formulation to the multi-objective HVRPTW is presented in this paper. A decision support system based on the meta-heuristics Multi-Objective Ant Colony System (MOACS) was also developed, which allows analysis of Pareto optimal solutions in the HVRPTW. VRPTW instances including those in HazMat transportation in general are validated on benchmark instances such as Solomon benchmark instances [38]. Although the instances cover several aspects of VRPTW, the instances are mainly based on Euclidean distances and fail to represent the actual traffic characteristics of the road network [35]. ITS and its incorporation with Geographical Information System (GIS) has enabled acquisition of real information of the road network, enabling more realistic analysis of logistics studies. Basically, three types of information: a) geographic information from GIS b) real static traffic information and c) on-line (real time) traffic information have been used. Although in the HazMat transportation literature, there are several GIS based studies that focus mainly on geographical characteristics of road network, less importance has been placed on historical and on-line real traffic information. ITS data-based HazMat instances with real

3

road traffic information provide excellent base for evaluating Pareto fronts in the proposed research, as it becomes feasible to graphically represent the realistic trade-off occurring between the cost and risk objectives and their potential contribution to social and economic effects on urban areas. Therefore, the proposed algorithm for HVRPTW has been tested in an application to a realistic HazMat logistics instance derived from Osaka road network for distribution of Liquefied Petroleum Gas (LPG). Travel times used for links in the instance are based on historical traffic information provided by VICS. Risk is obtained based on the worst case death consequence tested for various LPG release incidents (boiling liquid expanding vapor explosion, flash fire, jet fire and vapor cloud explosion) for the area. GIS is used for the necessary geographical data and for the final graphical displays, highlighting the trade-off between various Pareto optimal HVRPTW solutions. 2. Literature review All available models in the HazMat routing, both single objective and multi-objective, can be basically classified into two categories: the shortest path models and the Vehicle Routing and scheduling Problem (VRP) models. While the first class is dominated with large volume of the literature, research relating to second class is very limited. Shortest path models are designed to find shortest path between a given origin and destination node. Most multi-objective shortest path models considered both cost and risk attributes in the HazMat routing. A few models basically focused on risk attributes considering various risk factors. Both deterministic and static [28,11,6,7,27,45], and stochastic and dynamic [9,16,30,32] studies have been presented. The VRP arises naturally as a central problem in the fields of transportation, distribution and logistics [10]. As mentioned earlier, the complex characteristics of the HVRPTW have limited its studies. To the author’s knowledge, in past, only few routing studies with single objective [33,41] and multi-objective [27,46,47] have been presented in this direction. All the existing literature in the field is limited to simplification of the problem either by performing single objective optimization [33,41] or by transforming the multiple objectives into a single objective using weighted sum approach [27,46,47] for both path choice and routing. Tarantilis and Kiranoudis [41] proposed a risk based technique to solve a Capacitated Vehicle Routing Problem (CVRP) in the HazMat distribution. Routing decisions were made in a new space, which they called risk space. A List Based Threshold Accepting (LBTA) algorithm was used to solve the routing problem. Pradhananga et al. [33] applied single objective Genetic Algorithm (GA) to solve a typical HazMat routing problem in Thailand. The routing problem was solved for a single cost based objective. Considering hard time windows constraints, Zografos and Androutsopoulos [47] proposed a multi-objective HVRPTW model. However, a weighted sum approach (i.e. a linear combination of risk and time) was used to solve the problem. Dijkstra’s shortest path algorithm was used to find shortest paths between all vertices based on the combined (single) objective function. The resulting HVRPTW instance was then solved for optimal routing using an insertion-based heuristic approach. Varying the weight values for each objective, so called non-dominated solutions were obtained. In Ref. [46], they extended the model and the heuristic algorithm for developing a GIS-based decision support system for the integrated HazMat routing and emergency response decisions. Multicriteria analysis of a real case VRP for routing gasoline in Mexico City is presented in Ref. [27]. Vehicle routes were analyzed for two criteria: the length and the exposure population, which includes both inhabitants’ exposure and exposed travelers. Weight sum approach was used to obtain the routes. Three typical routes: two

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

4

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

extreme routes with minimized length and minimized exposure population and one with equal weight values to the two criteria were analyzed. So far known, none of the studies has used Pareto optimization to deal with these problems. Contrary to the HazMat VRP presented in Refs. [46,47,27], the routing model presented in this paper is based on the concept of Pareto optimization and explicitly reflects the multi-objective nature of the problem during the whole process of transportation. Both path choice and routing in the proposed algorithm are solved for non-dominated solutions and are processed in single step. Furthermore, a GIS-based decision support system that can be used by various stakeholders based on their own perceived importance of risk and travel time is developed in this research. Among the three types of information mentioned in earlier section, used for analysis of realistic logistics instances, studies using geographic information from GIS [39,36,37] and real traffic information [8,40,2] are prevalent in literature. While use of real time traffic information [3,20,26] is a relatively current advancement in the field. All the HazMat VRP researches mentioned earlier use such information to carry out analysis in applications to actual network-based instances. On the contrary, use of real, on-line traffic information from ITS data are very limited in HazMat transportation and Lozano et al. [27] is the only HazMat VRP study that used such traffic information. The traffic data from ITS system directly affect the cost term and have potential to represent realistic trade-off occurring between the cost and risk objectives. VICS is an ITS application in Japan which provides traffic data such as travel time for every 5-min slot. Although it has been used in several traffic related studies, its application in logistics sector is not common. Recently, studies on probabilistic VRPTW [2] and VRP with soft time windows [35] were carried out on test instances based on VICS data to evaluate the results based on actual traffic conditions. The ITS data-based HazMat logistics instance considered in this paper utilizes the static real traffic information based on historical records from VICS and data on geographical distribution from GIS. Therefore, it provides better and realistic characteristics of both objectives of the HVRPTW (i.e. time and risk). 3. Model formulation 3.1. Problem definition Based on discussion in Section 1, HVRPTW in specific to this study can be defined as a problem of determining a set of Pareto optimal routes for a fleet of vehicles carrying the HazMat in order to serve a given set of customers satisfying the following conditions: a. Both path choice and routing must be multi-objective and should be carried out as a single-step process. b. Minimizing the total scheduled travel time (sum of the travel times between customers, unloading times and waiting times at customer locations) and the total risk value of the routes are the two main objectives to be minimized in routing. On the other hand, minimizing the total travel time and the total risk value are the objectives for path choice. c. All vehicles must start and end their routes at depot. The time of starting and ending the operations of vehicle must be within the time windows specified at depot. d. Demands of customers must be serviced within pre-specified time windows. e. Waiting due to early arrival is possible, while late arrival is not allowed at all. The final set of routes are however expected to shorten waiting, since minimizing the total scheduled travel time that includes waiting at customers is one of the objective used in the routing process.

3.2. Risk calculation The risk term in the HazMat transportation differentiates it from rest of the other transportation problems. Thus, minimizing risk is one of the primary objectives in the HVRPTW. Ample research with efforts to obtain better qualitative and quantitative risk models have been carried out in the past. Considering occurrence of a HazMat accident event H on a road link l, the accident can associates O types of releases, U types of incidents and Y types of consequences. Therefore, the total HazMat transport risk associated with a road link l ðRTl Þ is expressed as:

RTl ¼

XXX

dl ðYy =H; Oo ; Uu Þdl ðUu =H; Oo Þdl ðOo =HÞdl ðHÞcyl

o˛O u˛U y˛Y

(1) where,

dl (H): Probability of a HazMat accident event H on link l dl ðOo =HÞ: Conditional probability of release event Oo on link l dl ðUu =H; Oo Þ: Conditional probability of incident Uu on link l dl ðYy =H; Oo ; Uu Þ: Conditional probability of consequence Yy on link l cyl : Possible y-type consequence Several consequences of a HazMat accident can occur. Since safety to human life has the top priority, the focus here is on the death consequence (D). Equation (1) can be re-written as Equation (2) to obtain HazMat risk of death associated with link l ðRD Þ. l

RD l ¼

XX

dl ðD=H; Oo ; Uu Þdl ðUu =H; Oo Þdl ðOo =HÞdl ðHÞcDl

(2)

o˛O u˛U

where,

dl ðD=H; Oo ; Uu Þ: Conditional probability of death consequence D on link l cDl : Population in neighborhood of link l Setting dl ðD=H; Oo ; Uu Þ as 1, dl ðD=H; Oo ; Uu Þ cD can be defined l as the exposure population on link l, where all the population suffer the same consequence of life loss on occurrence of a HazMat release incident. For the routing problem in this paper, we use the expected (worst case) loss measure of risk. It assumes occurrence of a worst HazMat release incident and uses the exposure population associated to it ðfl Þ. This allows the risk expression in Equation (2) to neglect the conditional probabilities. Therefore, HazMat transport risk associated with link l ðRl Þ can be expressed as:

Rl ¼ dl fl

(3)

where,

dl : Probability of a HazMat accident on link l Number of HazMat transport accidents in a given time period and the total distance traveled by HazMat trucks in the same period is used to determine HazMat accident rate on a unit road segment [31]. dl is then obtained multiplying the accident rate by the length of the link. fl is obtained moving a danger circle of radius l [5] along the entire link between the two nodes. As fl corresponds to the worst HazMat release incident, l value is defined comparing consequences of all the possible HazMat release incidents. The value of l is dependent upon the particular HazMat class under consideration.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

3.3. Mathematical formulation

XXX

Let us consider an urban road network (N, L) of nodes and links. b includes the depot node, a The node set N ¼ f1; 2; 3; ::::::::::; Ng b set of customer nodes and/or some non-customer nodes, where N

j˛V p˛P k˛K

denotes the total number of nodes. The link set L ¼ f1; 2; 3; ::::::::::; b Lg includes all available connections between nodes in N with b L as the total number of links. To each link l˛L connecting a pair of nodes, two attributes tl and Rl which represents its average travel time and risk, respectively, are defined. The HVRPTW is defined in a directed graph G(V, A) of vertices and arcs. The vertex set V includes the depot (vertex 0) and a set of b g. C b is the total number customer vertices C ¼ f 1; 2; 3; ::::::::::; C of customers. Unlike the two-step optimization where there is only one arc between each pair of vertices (i, j) ˛ V, the single-step model preserves various arcs (non-dominated paths) between the pairs. Therefore, each arc in this formulation is represented by a three index notation namely (i, j, p), where p is the unique identity associated with each non-dominated path commencing from i to j. Therefore, A is the set of all feasible arcs (i, j, p), (i, j) ˛ V and p ˛ P, where P is the set of all non-dominated paths with each path p b The nonhaving its unique identity in P ¼ f1; 2; 3; ::::::::::; Pg. dominated paths are assumed to be determined beforehand using labeling algorithm to be described in Section 4.2. For each arc (i, j, p), the following attributes are defined:

tijp : Travel time P of arcði; j; pÞ tl ¼

Rijp : Risk of Parcði; j; pÞ Rl ¼ l˛ði;j;pÞ

Similarly, let: Di : Demand at vertex i; i ˛C si : Unloading time at vertex i; i ˛V ¼ 0; for i ¼ 0 bi : Start of time window at vertex i; i ˛V ei : End of time window at vertex i; i ˛V K: Set of vehicles at depot Wk : Capacity of vehicle k; k ˛K Z1 : Total scheduled travel time of all the vehicles in operation Z2 : Total risk exposure associated with the transportation process xijpk : ¼1; if vehicle k uses arc (i, j, p) ¼ 0; otherwise cijp : Service travel time for arc (i, j, p) ¼ si þ tijp tiw : Waiting time at customer i; i ˛C Tik : Time at which vehicle k begins servicing at customer i; i ˛C TkA : Arrival time of vehicle k at depot TkD : Departure time of vehicle k from depot In the following, mathematical formulation to the HVRPTW is provided:

Z1 ¼

XXXX i˛V j˛V p˛P k˛K

Z2 ¼

XXXX

xijpk cijp þ

X

tiw

(4)

i˛C

xijpk Rijp

(5)

j˛C

(6)

i˛V j˛V p˛P k˛K

subject to

XXX i˛V p˛P k˛K

xijpk ¼ 1;

i˛C

(7)

x0jpk ¼ 1;

k˛K

(8)

xi0pk ¼ 1;

k˛K

(9)

j˛C p˛P

XX i˛C p˛P

XX

xihpk 

i˛V p˛P

X i˛C

0 Di @

XX

xhjpk ¼ 0;

k˛K; h˛C

(10)

j˛V p˛P

XX

1 xijpk A  Wk ;

k˛K

(11)

j˛V p˛P

  Tik þ si þ tijp þ tjw  Tjk  1  xijpk M;

i; j˛C; p˛P

k˛K (12)

  Tik þ si þ ti0p  TkA  1  xi0pk M;   TkD þ t0jp  Tjk  1  x0jpk M; bi  Tik  ei

l˛ði;j;pÞ

Minimize

XX

xijpk ¼ 1;

5

i˛C

i˛C;

j˛C;

p˛P

p˛P

k˛K

k˛K

(13) (14) (15)

b0  TkA  e0

k˛K

(16)

b0  TkD  e0

k˛K

(17)

xijpk ˛f0; 1g tiw ; Tik ; TkA ; TkD

i; j˛V; k˛K ˛<þ

i˛C; k˛K

(18) (19)

Equations (4) and (5) express objective functions for minimizing the total scheduled travel time and the total risk objectives of the transportation process, respectively. The total scheduled travel time in Equation (4) includes waiting times of the vehicles at all customers. xijpk is a binary decision variable (Equation (18)) that determines whether an arc (i, j, p) is used in the solution (xijpk ¼1) or not (xijpk ¼0). The objective function of the HVRPTW is subjected to several constraints similar to that of the VRPTW. Constraints (6) and (7) enforce each customer to be serviced once by a unique vehicle and a unique arc. Equations (8)e(10) are the flow conservation constraints to the problem while Equation (11) represents the capacity constraint to all vehicles and indicates that the total demand of all customers to be satisfied with a particular vehicle should not exceed its capacity. The HVRPTW formulation presented in this paper impose hard time windows constraints at each customer. While servicing by vehicles that arrives after the end of the time window is not allowed, vehicles that arrive earlier are allowed to wait to the start of the time windows. Equations (12)e(14) assure that if a vehicle passes from vertex i to j using arc (i, j, p), the time it starts servicing at j is the service start time at i plus sum of the unloading time at i, the travel time of arc (i, j, p) and waiting time (if any) at j, where M is a big constant. Waiting time at the customer j would be a positive value if vehicle k reaches customer j earlier than its start of the time window and would be 0 otherwise. Equations (15)e(17) assure fulfillment of time windows constraints

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

6

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

Step 1: Perform labeling algorithm at all vertices (customers and depot) to find set P and define G(V, A) Step 2: Set parameters. Initialize trail pheromone value and Pareto optimal set S based on initial solution using nearest neighborhood search. For Iteration = 1, initialize pheromone values for all arcs in A based on the trail pheromone value. Step 3: For ant = 1 to number of ants, construct solution. Construct solution: A: depotno =1, place ant at depot, wherein i = 0. B: Identify feasible set of customer vertices from i that is V’ and the set of arcs A’ connecting i to V’. If V’ is empty, add all arcs from i to depot vertex to A’. C: Select an arc (i, j, p) from A’. If exploitation, choose an arc (i, j, p) with maximum heuristic and pheromone product value. Else, choose randomly using probability of selection. D: Local update of pheromone value of the chosen arc (i, j, p). E: If the vertex corresponding to selected arc (i, j, p) is depot, set depotno = depotno +1, i = 0. Else, i = customer vertex corresponding to selected (i, j, p). F: If more customers to be served, go to step B. G: Set Total vehicles = depotno-1 Step 4: If ant < number of ants, ant = ant +1 and go to Step 3. Step 5: Update S based on dominance rule. Step 6: Apply insertion local search to each solution in S. Step 7: Find new trail pheromone value based on average objective values of current S. Step 8: If previous trail pheromone value < new trail pheromone, set trail pheromone = new trail pheromone and initialize pheromone values of all arcs in A based on pheromone value of this trail. Else, carry out global update of pheromone for all arcs (i, j, p) that belongs to solutions in S. Step 9: Iteration = Iteration+1, go to Step 3 until Iteration <10000

Fig. 2. MOACS-based meta-heuristic algorithms for HVRPTW.

at customers and depot. Service beginning time at customers which must be a positive value (Equation (19)) must satisfy the specified time windows at the corresponding customer locations. Similarly, the arrivals and departures of the vehicles from depot must be done within time window specified at the depot vertex. 4. A meta-heuristic approach for HVRPTW HVRPTW as formulated in Section 3 is much more difficult to solve than other multi-criteria VRPTW due to necessity of singlestep optimization in addition to requirement of determining a set of equally valid optimal solutions (i.e., Pareto set). Meta-heuristics are the most practical alternatives to solve such problems. GarciaMartinez et al. [22] showed competency of Multi-Objective Ant Colony Optimization (MOACO) over well-known evolutionary algorithms, Non-dominated Sorting Genetic Algorithm (NSGA) and Strength Pareto Evolutionary Algorithm (SPEA) to bi-criteria traveling salesman problem. This characterizes suitability of MOACO algorithms for solving complex multi-objective combinatorial optimization problems such as the HVRPTW. Hence, a new MOACO algorithm based on the MOACS has been presented in this section for solving the HVRPTW. Ant Colony Optimization (ACO) is a meta-heuristic inspired by the evolutionary mechanism of artificial ants to search for good paths. Ant System (AS) by Dorigo et al. [15] was the first ACO algorithm. Ant Colony System (ACS) by Dorigo and Gambardella [14] is a variant of AS that has been preferably used in the literature of multi-objective optimization. A number of efforts (e.g., Refs. [4,13,21,29]) have been carried out for developing MOACOs to solve multi-objective transportation problems. Among the various MOACO approaches, Multi-Objective Ant Colony System for VRPTW (MOACS-VRPTW) proposed by Baran and Schaerer [4] is the only approach that solved VRPTW for non-dominated set of solutions. Moreover, the survey carried out by Garcia-Martinez et al. [22] specified that the algorithm shows better convergence to the Pareto front than most other MOACO approaches. Hence, the proposed algorithm in this research adopts its basic structure

from MOACS-VRPTW, however, for utilizing all the non-dominated paths based on travel time and risk objectives for the path choice and carrying out single-step optimization, the proposed algorithm shows some significant changes over MOACS-VRPTW. Another major difference is that the proposed algorithm also includes a local search step to improve the quality of Pareto optimal solutions. Before proceeding to the solution construction step, labeling algorithms are carried out at all customers and depot vertices to obtain set P and thereby to obtain the set A. Detail on labeling algorithm is presented in Section 4.1. The set P includes all the non-dominated paths for path choice. Ants in the proposed MOACS use these paths to construct their routes (solution). A local search method is employed to improve quality of solutions constructed in each iteration. Details on the solution construction process and the local search method are provided in Sections 4.3 and 4.4, respectively. Fig. 2 presents flowchart of the proposed MOACS algorithm for HVRPTW.

4.1. Labeling algorithm Labeling algorithms serve as an iterative approach for solving shortest path problems. These algorithms assign labels to nodes at each step. Label at a node stores information about a path leading to it by maintaining a connection with the label at the predecessor node. Basic knowledge on the labeling algorithm for shortest path problems is given in detail in Ref. [1]. Labeling algorithm used in this paper is based on the template labeling algorithm for shortest b extra path problem with resource constraint [23]. To each label, N resources are applied as in Ref. [19]. These resources indicate connectivity of a node with the resident node in the label as well as provide information on whether it has been already visited or not in the label. Labeling algorithm with travel time and risk value as b þ 1 times to obtain the the resource constraints is performed C non-dominated paths for path choice between all the vertices. At each time, the algorithm starts at a different origin, which could be the depot vertex or either one of the customer vertices.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

b þ 5 dimensions Label U created at a node is a row vector of N and is represented as:

U ¼ fRESðUÞ; tðUÞ; RðUÞ; VISðUÞ; PREDðUÞ; LABEL NO:g where, RES(U): Resident node of the label U tðUÞ: Travel time resource of the label, which shows the cumulative travel time along the path, p(U) ending at RES(U) R(U): Risk resource of the label, which shows cumulative risk value along the path, p(U) ending at RES(U) b entries of 2 for all VIS(U): Visited node vector of the label with N nodes already visited, 1 for all unreachable nodes due to lack of connectivity and 0 for all reachable nodes PRED(U): Immediate predecessor node of RES(U) on p(U), represented by the label number of the label with RES(U) ¼ predecessor node LABEL NO.: A label counter assigning a unique number to each label. The algorithms work on dynamic manipulation of two sets, namely set of unprocessed labels (L) and set of useful labels (U). Label selection, dominance rule and path extension are repetitively carried out until all labels in L are processed (i.e. L ¼ Q, where Q is an empty set). The set U, finally obtained is then used to trace the paths connecting origin vertex with the other vertices in label filtration process to obtain set P. 4.2. Initialization An initial routing solution is obtained using nearest neighborhood (nn) heuristic. The solution obtained itself represents the first member of the Pareto optimal set S. s0 defined as the initial trail pheromone value for the first iteration in ACS can be any heuristics value but it is good to start with some prior knowledge for better convergence. Equation (20) is the expression to obtain it, which is similar to the initialization used in MOACS algorithms for VRPTW [4].

s0 ¼

1 ðjNjÞnn $ðZ1 Þnn $ðZ2 Þnn

(20)

where, ðjNjÞnn : Sum of total number of customers and number of vehicles used in nearest neighborhood solution ðZ1 Þnn : Total scheduled travel time in nearest neighborhood solution ðZ2 Þnn : Total risk value in nearest neighborhood solution For computation of the initial pheromone values for subsequent iterations ðs00 Þ, all these three terms are obtained based on the average number of vehicles used and the average objective values for the Pareto optimal set of the previous iteration.

and time windows allow. Additional vehicles are used to service all the customers. All vehicles in operation end their routes at depot. Each ant m at vertex i˛V selects an arc (i, j, p) among A’ based on a pseudorandom proportional rule. A’ is a subset of A containing all arcs from vertex i to all the feasible set of customer vertices (V’). In case when V’ is empty, the depot vertex is added in it. A random number q ð0  q  1Þ is generated, and the route is extended using exploitation process (if q  q0 ) or exploration process (if q > q0 ). Here, q0 is a parameter that defines the relative importance of exploration and exploitation. For increased value of q0, the probability to use exploration is reduced. The maximum value of pheromone and heuristics combination represented by variable (D) in Equation (21) is used for exploitation of previous best choice of arc at vertex i. Several choices of other feasible arcs are explored using probability of choosing arc(i, j, p) (prði; j; pÞ ) in Equation (22).



h

D ¼ maxði;j;pÞ˛A0 sijp $ hijp

iqh b   n  qm $ nijp

iqh b   n qm $ nijp h i qh b   n qm sijp $ hijp $ nijp

(21)

h

prði; j; pÞ ¼

sijp $ hijp P ði;j;pÞ˛A0

(22)

where,

sijp : Pheromone value of arc (i, j, p) hijp : Heuristic value relating to scheduled travel time objective of arc (i, j, p) nijp : Heuristic value relating to risk objective of arc (i, j, p) b: Parameter that define the relative influence of scheduled travel time objective with respect to the pheromone m: Parameter that define the relative influence of risk objective with respect to the pheromone qh : Ant specific weight/preference for normalizing scheduled travel time objective; ¼ ðm=ajm ¼ 1; ::: ; aÞ a: Parameter that defines number of ants qn : Ant specific weight/preference for normalizing risk objecn h tive; q ¼ 1  q

hijp in the proposed algorithm is obtained using concept proposed earlier by Gambardella et al. [21]; which allows ants to construct their solutions with maximum utilization of the time windows. Equation (23) is the expression for calculation of the time-related heuristic value of arc (i, j, p). In addition to the scheduled travel time necessary for vehicle k in passing arc(i, j, p), the second expression in the equation also considers urgency of servicing customer j given by the time interval between the current time at customer i and the one in which the time window at j closes. For case with i ¼ 0 that is when ant is at depot, Tik is taken as the start time window at depot and si as zero.   Tjk ¼ max Tik þ si þ tijp ; bj

hijp ¼

4.3. Solution construction Set P, obtained from labeling algorithm is used for constructing solutions in the proposed MOACS algorithm for HVRPTW. Contrary to the algorithm proposed by Baran and Schaerer [4] in which ants construct their solution with insertion of feasible customers, ants in the proposed algorithm construct their solution with insertion of feasible arcs in A. Each ant m begins its route with a vehicle k stationed at depot and goes on satisfying customers until the capacity

7

1     max 1; Tjk  Tik ej  Tik

(23)

The risk-related heuristic value relies on the risk value associated with the arc (i, j, p) and is evaluated as given in Equation (24).

nijp ¼

1 Rijp

(24)

Once an arc (i, j, p) is selected for insertion, the next vertex j to be visited is recognized. Also, the service beginning time Tjk and the

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

8

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

waiting time tjw (if any) at the customer j are known. Objective values for the portion of the route (cijp þ tjw and Rijp ) can then be determined. TkD and TkA are determined once service beginning times Tjk for all customer vertices j˛C are known. Finally, the total number of vehicles used in the solution and the total objective values, Z1 (sum of cijp þ tjw of all used vehicles, for all inserted arcs (i, j, p)) and Z2 (sum of Rijp of all used vehicles, for all inserted arcs (i, j, p)) are determined. h n It should be noted that b, m, q and q are MOACS specific parameters that determine relative influence of pheromone and h n heuristics information. q b and q m therefore, should not be confused with the objective weights used in weight sum approach. h n Moreover, different ants in the algorithm have different q b and q m values therefore, exploit or explore the solutions for different pheromone and heuristics combinations while the non-dominated set of solutions is decided based on the dominance rule. Once all ants in the iteration form their routes, Pareto optimal set S is updated checking its solutions for dominance to the newly constructed solutions. Dominance rule used in this algorithm separately compares both objective values of the routing solutions and removes the solutions that are inferior in terms of both objectives. Solution with lower number of vehicle is preferable from cost aspect even if they are inferior in terms of both objectives. Hence, apart from checking dominance in terms of the two objectives, number of vehicles used is also considered as a criterion for selection of solutions in S. Set S therefore is obtained discarding solutions that are inferior in terms of both the objectives while using equal or higher number of vehicles. 4.4. Local search Local search procedures are often used to improve the quality of solutions in meta-heuristics. Therefore, unlike the MOACS-VRPTW, MOACS for the HVRPTW presented in this research also employs a local search method. After updating the Pareto optimal solution in the iteration, an insertion-based local search is applied to all solutions belonging to S. All vertices in solution j that belongs to set S are given chances to be inserted to the same vehicle route at different position or to the route of other vehicles without violating feasibility requirements. The newly obtained solutions are then checked for improved objective values as well as for reduction in use of vehicles. Pareto optimal set S is revised accordingly.

pheromone to increase pheromone of each arc (i, j, p) of j˛S as shown in Equation (26) before proceeding to the next iteration. The increase r=½ðZ1 Þj $ðZ2 Þj  which is of range s0 0 is also used by other similar Pareto based MOACS such as presented by Baran and Schaerer [4].

h

s0ijp ¼ ð1  rÞs0ijp þ r= ðZ1 Þj $ðZ2 Þj

i

ði; j; pÞ˛j

(26)

where, ðZ1 Þj : Total scheduled travel time of solution j˛S ðZ2 Þj : Total risk value of solution j˛S

5. An application of the multi-objective HVRPTW 5.1. Test instance The test instance is a realistic HVRPTW logistics instance, which is derived from the road network data in Osaka, Japan, and deals with distribution of LPG. There are 225 nodes and 781 links in the network. Each one way road is modeled with a unidirectional link and each two way road is modeled by two unidirectional links. Out of 225 nodes in the network, 24 are customer vertices represented as customer 1 to customer 25 except customer 17, which is the depot vertex. Fig. 3 shows the road network and the locations of customers and the depot in the instance. Time windows at the customers are randomly generated, and a demand value of 3300 L is assumed at each customer location. All vehicles stationed at the depot are assumed to have a uniform capacity of 20,000 L. Service time to unload the material is assumed to be 9 min at each customer location. All the links in the instance are equipped with the VICS technology, which provides traffic information such as travel time for each 5-min slot. Since the HVRPTW presented in this paper is a static model, the travel time for each link is obtained as the static average value of all 24 h slots in weekdays over a month. All the links are identified in a GIS map of Osaka as shown in Fig. 3 to identify the geographical distribution of population in the area. This population data is used to calculate risk associated with each link in

4.5. Pheromone update To each arc (i, j, p) in the HVRPTW complete graph, used by an ant for solution construction, local pheromone update is carried out as given in Equation (25).

s0ijp ¼ ð1  rÞs0ijp þ rs0

(25)

where,

s0ijp : Pheromone value of arc (i, j, p) for current iteration s0ijp : Pheromone value of arc (i, j, p) for next iteration r: Evaporation coefficient which powers exploration process by evaporating trail pheromone values of used arcs. Based on the revised Pareto optimal set S after the local search procedure, new initial trail pheromone value s0 0 is determined as defined in Section 4.2. If set S is obtained with better s0 0 than the previously initialized trail pheromone value s0 , pheromone value equal to s0 0 is assigned to all arcs in the set A to improve exploration. Otherwise exploitation is favored by global updating

Fig. 3. Road network and customer locations in realistic HVRPTW instance.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

the instance, which is described in the following subsection in detail. 5.2. Risk calculation in test instance The test network is mainly spread within twenty four wards in Osaka city while a few links are located in five other cities in the Osaka prefecture. Statistics of HazMat freight vehicles accidents and the vehicle kilometers traveled are unavailable for Osaka Prefecture. Therefore, the HazMat accident rate on a unit road segment in Osaka Prefecture is assumed to be equal to that of a normal freight vehicle’s accident rate, and is obtained using the freight vehicle accident statistics and vehicle kilometers traveled. In reality, each link may have a different probability of the HazMat accident depending upon the road and vehicle conditions, driver’s behavior and the surrounding environment. Historical records of accidents in the neighborhood area to some extent represent these differences. Therefore, the HazMat accident rate on unit road segment of Osaka Prefecture is modified with a coefficient that depends upon the specific number of accidents in a ward in the Prefecture. dl for each link l in the instance is then obtained using these modified rates and the corresponding link length with an assumption that all links within a ward are similar (i.e. pose same probability to a HazMat accident). These assumptions are useful for illustration purposes, but may not match the exact actual situations. The traffic accident statistics are obtained from the Institute for Traffic Accident Research and Data Analysis (ITARDA), Japan while the vehicle kilometer traveled is obtained from and the Ministry of Land, Infrastructure, Transport and Tourism (MLIT), Japan. The impact radius l corresponds to the death consequence of the worst HazMat release incident. Therefore, l for the instance is obtained testing threat distances for various LPG release incidents such as Boiling Liquid Expanding Vapor Explosion (BLEVE), flash fire, jet fire and vapor cloud explosion, that are modeled for Osaka, Japan using Areal Locations of Hazardous Atmospheres (ALOHA) [44]. The value of l for the instance is determined to be 275 m in all direction, from the above mentioned worst-case analysis. Impact area created on each link l in the network using this value is then multiplied with the population density of the associated ward to obtain fl for the link. Finally Rl for the link is obtained using Equation (3). 5.3. Results and discussions Although our algorithm is basically designed to solve multiobjective problems, multi-objective benchmark instances relating to the problem dealt with in this paper are unavailable in the literature. In order to investigate the quality of the objective values obtained with the proposed algorithm, with respect to those obtained with the best known solutions in the literature, numerical tests were carried out on Solomon benchmark instances using single objective of time only, similar to VRPTW. It was found that the algorithm works satisfactorily with average discrepancy of 2.43% and worst gap of 10.66% from the best known optimal solutions to the instances (for detail, please refer Ref. [34]). Furthermore, a sensitivity analysis of the parameters in specific to the multi-objective HVRPTW was also carried out performing several numerical tests on the realistic HVRPTW instance discussed in Section 5.1. The algorithms were executed using Borland C compiler in Core 2 Duo desktop PC of 2.67 GHz with 2 GB RAM. In previous studies [4,21], the best computational results of MOACS for VRPTW were obtained with parameter setting of a ¼ 10, b ¼ 1, q0 ¼ 0.9 and r ¼ 0.1. These values along with a m ¼ 1 and a stopping criterion of

9

number of iterations, u ¼ 10,000 are used as the basic parameter setting for the instance. Section 5.3.1 shows the result of the sensitivity analysis of the parameters. Detail analysis of the routing solutions is presented in Section 5.3.2. 5.3.1. Sensitivity analysis Sensitivity of the Pareto front obtained for the realistic HVRPTW instance is tested against several parameter values. The basic parameter setting is used as the reference setting. Therefore, while testing different values of a parameter, the values of all other parameters are held constant as in the reference setting. Fig. 4(a)e(f) shows effects of parameters: u, a, b, m, q0, and r to the quality of routing solutions of the instance when varied between [1000, 15,000], [2, 20], [0.5, 2], [0.5, 2], [0.1, 0.9] and [0.1, 0.9], respectively. Changes of u and a parameters are also observed to significantly affect computation times of the Pareto fronts in the instance. On the contrary, computation time is found be less sensitive to the changes in parameters b, m, q0, and r. The Pareto front of the instance in the basic setting is observed to have two groups of routing solutions: using 5 vehicles and using 4 vehicles. The Pareto fronts obtained for various parameter settings are therefore compared for the two cases (i.e. cases with 4 vehicles and 5 vehicles used) of the routing solutions. Fig. 4(a) shows the Pareto fronts obtained for u of 1000, 5000, 10,000 and 15,000 iterations. The corresponding computation times in seconds are 197, 312, 1412 and 1857, respectively. Although, the Pareto fronts obtained for the instance with lower values of u (1000 and 5000) have computation times much lower than that for higher values of u (10,000 and 15,000), quality of the Pareto fronts obtained for the higher values are found to be significantly improved. Furthermore, comparison of the Pareto fronts obtained for u ¼ 10,000 and 15,000 shows that the increase of u beyond 10,000 iterations in the present instance is insignificant as it does not have any enhancement in the solution quality while the computation time is still increasing. Increase in number of ants improves solution while it also increases computation time. Fig. 4(b) shows the Pareto fronts obtained for a value of 2, 5, 10 and 20 with computation times in seconds of 1141, 967, 1412 and 2343, respectively. The figure shows that solutions of the instance mainly those obtained for the case of 4 vehicles used are inferior for lower values of a (2, 5). For larger values of a, better Pareto fronts are obtained. However significant increase in computation time is observed mainly while increasing a from 10 to 20. For increased value of b, ant favors larger savings in the scheduled travel time objective. Therefore, the chances to have solutions with better scheduled travel time increases. Nevertheless, higher the value of b, more is the chance that the algorithm will be trapped in a local optimum. Fig. 4(c) shows sensitivity of the solutions with b parameter. Increase of b from 0.5 to 1 for the test instance is found to be beneficial as the part of the Pareto front corresponding to the case of 4 vehicles used is observed to have significant improvement resulting solutions with lower scheduled travel time values. The further increase of b however is found to deteriorate the Pareto front for both cases of solutions due to early convergence. Fig. 4(d) shows sensitivity of the solutions with m parameter. Increase of m tends to provide better solution mainly in terms of risk value since the relative importance of the risk saving is increased. However, there is added possibility of early convergence. For the test instance, increase of m from 0.5 to 1 is observed to improve the Pareto front for both cases of the solutions. Further increase of m continues slight improvements in the solutions corresponding to case of 5 vehicles used however solutions corresponding to the case of 4 vehicles used are observed to be deteriorated due to the early convergence with the increase.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

10

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

Fig. 4. Sensitivity of parameters: (a) u (b) a (c) b (d) m (e) q0, and (f) r.

Large q0 favors exploitation of the information from previous best solution. Fig. 4(e) shows Pareto fronts obtained for the instance for various q0 values. The Pareto fronts shows that the solutions of the instance mainly those corresponding to case of 4 vehicles used are very sensitive to the increase of q0 parameter and are significantly better for high value of q0 than for its lower values. Fig. 4(f) shows sensitivity of the solutions with parameter r. Solutions of the instance corresponding to the case of 4 vehicles

used are found very sensitive to the increase of parameter r. The solutions however are observed to be deteriorated with the increase. Coefficient r helps exploration of new arcs. However, high value of it can result a fast focus on the new information and some previously very attractive choices can be forgotten quickly. Therefore, the algorithm can be trapped into a local optimum as is observed for higher values of r for the test instance.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

Case of 5 vehicles used Solution with minimized scheduled travel time (single objective) Case of 4 vehicles used Solution with minimized risk value (single objective) 550

Risk value (x10-2)

Overall, among the two cases of the solutions, solutions corresponding to the case of 4 vehicles used are found to be more sensitive to the parameters changes. Comparison of the Pareto fronts obtained for various parameters setting shows that the Pareto fronts obtained with two settings: the basic setting and the setting with u ¼ 10,000, a ¼ 20, b ¼ 1, m ¼ 1, q0 ¼ 0.9, and r ¼ 0.1 are comparatively better than those obtained for the rest of the settings. Pareto front obtained with the later setting have a few better solutions compared to that with the basic setting, however its computation time is also 1.5 times of that of the basic setting. Therefore, parameter setting of the study is kept fixed to the basic setting u ¼ 10,000, a ¼ 10, b ¼ 1, m ¼ 1, q0 ¼ 0.9, and r ¼ 0.1. Solution analysis in the next section corresponds to the numerical tests carried out with the basic setting.

11

500 450 400 350

5.3.2. Routing results This section provides detailed analysis of the multi-objective routing solutions obtained for the realistic instance. A total of 3041 dominant paths (arcs) are obtained as the result of the labeling algorithm for the instance. This is significantly higher than 600 arcs in the typical two-step/weighted sum approach, where each pair of vertices is connected by one arc only. Therefore, the complexity of the single-step/Pareto optimal approach is significantly higher than the typical case, as up to 32 arcs are observed connecting customer 18 to customer 20. Routing solutions for the instance presented both in Sections 5.3.1 and 5.3.2 are obtained solving the HVRPTW using all these arcs. To analyze the Pareto front and the consistency of the results in different test runs, the problem is solved for ten consecutive runs with basic parameter setting. The Pareto fronts obtained in all test runs are similar in terms of distribution of the Pareto optimal solutions and are found to have solutions using 4 and 5 numbers of vehicles. This means that all solutions for the instance using more than 5 vehicles are dominated by solutions using 5 and 4 vehicles in terms of both the scheduled travel time and risk objectives. Fig. 5 shows the typical Pareto front obtained in test run 1. To demonstrate how effective the algorithm is to capture the wide range of Pareto optimal set, the problem is also solved for minimized total scheduled travel time and minimized total risk values, respectively, using single objective ant colony system. The triangle and the diamond in Fig. 5 are respectively, the resulting solutions which serve as a benchmark for minimum total scheduled travel time and total risk objectives for all test runs. It is very important to ensure that the algorithms can provide Pareto optimal sets that have large coverage, contain solutions that are equally distributed and have good (dominant) objective values. Pareto fronts in most test runs have Pareto optimal solutions as close as to the benchmark solution obtained in relation to risk minimization as in Fig. 5. In all Pareto fronts, solutions with time values of range 1030e1100 and with risk values of range 510(102) to 540(102) are missing. These slight gaps in the Pareto fronts can be argued due to the shortage of feasible solutions at these locations in specific to this problem. Moreover, the Pareto fronts capture a large portion of well-distributed central Pareto optimized solutions. In all test runs, it is observed that the solutions in the top-left section of the Pareto front (with lesser time values and higher risk values) are abundant with those using 5 vehicles. On the other hand, solutions in the bottom right section (with higher time values and lesser risk values) are mostly using 4 vehicles. Reduction in the number of vehicles in most cases is achieved by accommodating more number of customers by each vehicle. This often leads to the requirement of waiting for earliest possible service beginning time defined by the time windows. This explains the lower values of the total scheduled travel time observed with the 5 vehicle cases. We

300 730

830

930

1030

1130

Scheduled travel time (min) Fig. 5. Pareto front in realistic HVRPTW instance.

may get results same as given for the case of 4 vehicles used and 5 vehicles used if we keep fixed cost of vehicle as a separate (third) objective. Finally it is up to the choice of decision maker to make decision between lower vehicles and related travel time and risk values. In future, fixed cost of vehicle can be included in the objective to emphasize the above mentioned trade-off. However, such optimization when the fixed cost of vehicle is very high can drive the whole optimization process towards minimizing the number of vehicles used and may result in both longer and high risk HazMat routes, which is not desirable. Because of the conflicting natures of the objectives in the problem, the solutions with less value of one objective often cause larger value of the other objective. For both cases of solutions (case of 4 vehicles used and case of 5 vehicles used), the Pareto front shows large number of Pareto optimal solutions representing trade-offs between scheduled travel time and risk objectives. The decrease of risk values, with increment in scheduled travel time is significant in case of 5 vehicles used. The trade-off is stabilized as the risk value reached around 400  102 units. For case of 4 vehicles used, the reduction in risk value with increase of scheduled travel time is not so large. It is up to the need and preference of the decision maker to choose the particular case of vehicles used. But in both cases, the central portions of the Pareto fronts are the areas that the decision makers may be more interested in, as it provides the better compromise between the time and risk values. Computation time is the other point that most decision makers are interested in. Column 10 in Table 2 provides computation times for the ten consecutive test runs carried out on the realistic instance. The computation time for the instance is found in the range between 700 s and 1400 s. Considering the large network size and the complex nature of the instance, computation time of range 700 se1400 s is within a reasonable range. In fact, most exact approaches are unable to solve such problems, and for heuristics approaches, the above range is quite acceptable. To check consistency of the final Pareto optimal solutions, minimum objective values of the Pareto optimal sets obtained in each test run are also provided. The values are provided for both cases of vehicles used. Columns 2 and 4 are the minimum scheduled travel time and risk values for the case of 4 vehicles used while Columns 3 and 5 are their corresponding risk and scheduled travel time values, respectively. Columns 6 to 9 give similar values for the case of 5 vehicles used. The Pareto optimal

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

12

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

Table 2 Minimum and maximum objective values for 10 test runs in realistic instance. Test run

1 2 3 4 5 6 7 8 9 10

Case of 4 vehicles used

Case of 5 vehicles used

Min. scheduled travel time (a) (min)

Risk of (a) (102)

Min. risk (b) (102)

Scheduled travel time of (b) (min)

Min. scheduled travel time (c) (min)

Risk of (c) (102)

Min. risk (d) (102)

Scheduled travel time of (d) (min)

919.5 972.9 934.2 918.8 912.7 933.3 929.3 924.3 931.7 925.0

392.3 388.5 445.7 427.1 402.3 400.6 385.1 387.0 387.0 420.5

339.2 358.4 338.2 339.4 362.8 349.8 340.9 339.8 359.1 349.9

1108.3 1080.7 1111.7 1108.2 1081.1 989.8 1107.6 1014.5 1029.1 1153.3

801.2 825.1 827.3 775.2 799.8 772.5 768.4 762.7 794.8 773.5

509.2 463.7 491.9 511.2 460.1 522.5 488.4 493.5 455.2 499.9

369.4 374.7 378.0 368.2 390.0 377.8 382.3 368.7 392.1 385.9

906.4 904.4 899.3 907.3 867.6 883.2 927.7 907.8 893.4 890.4

sets found in all test runs are found consistent in terms of the objective values. In all sets, the minimized risk values are found in association with solutions using 4 vehicles. Similarly, the minimized scheduled travel times in all sets are found with solutions using 5 vehicles. These results are also consistent with the benchmark solutions obtained for scheduled travel time objective (obtained for the case of 5 vehicles used) and risk objective (obtained for the case of 4 vehicles used). Pareto optimal sets in all test runs comprise of two categories of solutions. Category 1 includes solutions with different customers visiting order (sequence) and Category 2 includes solutions with exactly same visiting order but with different connecting arcs (i, j, p) between them. While the reasons of Category 1 to be in Pareto set is quite obvious, Category 2 is the result of use of all nondominated paths in single-step simultaneous optimization of path choice and routing, which is the novel idea presented in this paper. To further elaborate the results, Solution I, II and III are

Computation time (s)

1412 790 1193 986 782 1137 833 976 718 740

chosen as the three typical Pareto optimal solutions from Pareto optimal set, obtained in test run 1. Fig. 6 shows comparison of vehicle routes in the three solutions. All of them used 5 vehicles to deliver to 24 customers from the depot. 5.3.2.1. Category 1 e different customer sequence. Solutions I and II are the example solutions belonging to Category 1. As mentioned earlier, these solutions have different customer visiting order, hence providing different objective values. The total scheduled travel time value and the total risk value for Solution I are obtained to be 827.12 min and 430.12  102 units respectively while the values for Solution II are 852.55 min and 406.57  102 units respectively. Comparison of routes in Fig. 6(a, b, c and e) shows that more than 50% links are different on the routes of Solution I and II. Yet, both of these solutions are Pareto optimal with slight difference in objective values. Such Pareto optimal solutions, therefore, can

Fig. 6. Comparison of vehicle routes of the three Pareto optimal solutions: (a) Vehicle 1 (b) Vehicle 2 (c) Vehicle 3 (d) Vehicle 4 (e) Vehicle 5.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

provide suitable alternatives for decision making mainly in situation with failure of some of the links in the network or to alternate between different areas for avoiding continuous risk to the same population. However, it must be made clear that the model does not capture the equity aspect, and not much insight can be drawn from these results on the risk equity of the vehicle routes. Therefore, it is left to the decision maker to make choice of equitable routes. 5.3.2.2. Category 2 e same customer sequence different (i, j, p). Solution I and III on the other hand are the example solutions belonging to Category 2. These solutions have vehicles with exactly same customer visiting order, however, uses different arcs (i, j, p) to connect to same pair of vertices i and j, which results difference in objective values. As given in the previous section, Solution I has the total scheduled travel time value and the total risk value of 827.12 min and 430.12  102 units, respectively, while the values for Solution III are obtained to be 838.11 min, 421.34  102 units, respectively. Comparing routes of solution I and III, the paths followed are exactly same for Vehicle 3 to 5. On the other hand, Vehicle 1 and 2 in the two solutions follow different paths in the network despite of the fact that both vehicles serve to same customers in same order. Again, as Fig. 6(a, b) clearly shows about 50% of the routes use different links (given by the blue and green colored routes). As mentioned previously, this routing pattern is an outcome of using all the non-dominated paths in path choice. Because the difference of objective values for this category of solutions is lesser than that for the first category, these provide a better alternative for routing with differing link usage. 6. Conclusions This paper provides advancement over the literature on HazMat transportation studies by considering for the first time the concept of Pareto optimization to deal with multi-objective routing and multi-objective path choice in HVRPTW. Pareto optimization on contrary to the prevalently used weighted sum approach solves the HVRPTW routing problem to determine a set of non-dominated solutions. Practically, Pareto optimal solutions can provide a significant base to decision making process in HazMat transportation. Firstly it provides the decision maker with clear idea on the cost and risk trade-off occurring for a particular instance without generating need for defining weight values to the objectives beforehand. Secondly, decision makers of varying preference to the objectives can be benefited. Use of ITSbased HazMat logistics instance enhances the quality of analysis in realistic terms. Analysis of results signifies followings:  Numerical results on the HVRPTW test instance shows the importance of carrying out multi-objective path choice for determination of complete Pareto front for routing in the HVRPTW. Therefore, a new four index formulation to the HVRPTW is presented based on the multiple non-dominated connections i.e. arcs (i, j, p) between two vertices i, j.  Sensitivity analysis of the parameters shows best results for the HVRPTW instance for basic setting of u ¼ 10,000, a ¼ 10, b ¼ 1, m ¼ 1, q0 ¼ 0.9, and r ¼ 0.1. Despite the large size and complex nature of the realistic instance in the HVRPTW, high quality Pareto fronts (with large extent and centrally well distributed solutions) within reasonable computation time are obtained with the algorithm in the setting, which shows the algorithm’s applicability for solving practical problems in HazMat transportation.

13

The model has potential to guide alternative routes for HazMat transportation during restricted situations of road network such as during accidents, construction, disasters, and to alternate routes pass through different areas so that same population is not subjected to continuous risk. Noticing the high complexity of singlestep routing, the exposure population and the travel time information used in this study are static values obtained from historical data and the stochastic and time varying nature of the HazMat routing has not been taken into consideration at this stage. Extensions of the model considering such characteristics of the HazMat routing problem utilizing the real time ITS-based traffic information and considering the effects of land use characteristics of the road network are the future tasks to the study. References [1] Ahuja R, Magnanti T, Orlin J. Network flows: theory, algorithms, and applications. New Jersey: Prentice Hall; 1993. [2] Ando N, Taniguchi E. Travel time reliability in vehicle routing and scheduling with time windows. Networks Spat Econ 2006;6:293e311. [3] Ausiello G, Feuerstein E, Leonardi S, Stougie L, Talamo M. Algorithms for the on-line traveling salesman. Algorithmica 2001;29(4):560e81. [4] Baran B, Schaerer M. A multiobjective ant colony system for vehicle routing problem with time windows. In: Proceeding of 21st IASTED international conference on applied informatics. Austria; 2003. pp. 97e102. [5] Batta R, Chiu SS. Optimal obnoxious paths on a network: transportation of hazardous materials. Operations Res 1988;36(1):84e92. [6] Bonvicini S, Spadoni G. A hazmat multi-commodity routing model satisfying risk criteria: a case study. J Loss Prev Process Ind 2008;21:345e58. [7] Caramia M, Giordani S, Iovanella A. On the selection of k routes in multiobjective hazmat route planning. IMA J Manag Math 2010;21(3):239e51. [8] Crainic TG, Ricciardi N, Storchi G. Advanced freight transportation systems for congested urban areas. Transp Res Part C 2004;12:119e37. [9] Chang TS, Nozick LK, Turnquist MA. Multi-objective path finding in stochastic dynamic networks, with application to routing hazardous materials shipments. Transp Sci 2005;39(3):383e99. [10] Dantzig GB, Ramser JH. The truck dispatching problem. Manag Sci 1959;6:80e 91. [11] Dell’Olmo P, Gentili M, Scozzari A. On finding dissimilar Pareto-optimal paths. Eur J Oper Res 2005;162:70e82. [12] Desrosiers J, Dumas Y, Solomon MM, Sournis E. Time constrained routing and scheduling. In: Ball MO, Magnanti TL, Monma CL, Nemhauser GL, editors. Handbooks in operations research and management science: network routing, vol. 8. Amsterdam: North-Holland; 1995. pp. 35e139. [13] Doerner K, Hartl RF, Reimann M. CompetAnts for problem solving e the case of full truckload transportation. Cent Eur J Oper Res 2003;11(2):115e41. [14] Dorigo M, Gambardella L. Ant colony system: a cooperative learning approach to the travelling salesman problem. IEEE Trans Evol Comput 1997;1(1):53e66. [15] Dorigo M, Maniezzo V, Colorni A. The ant system: optimization by a colony of cooperating agents. IEEE Trans Syst Man Cybern Part B 1996;26(1):29e41. [16] Erkut E, Alp O. Integrated routing and scheduling of hazmat trucks with stops en route. Transp Sci 2007;41(1):107e22. [17] Erkut E, Tjandra SA, Verter V. Hazardous material transportation. In: Laporte G, Barnhart C, editors. Handbooks in operations research and management science: transportation, vol. 14. Amsterdam: North-Holland; 2007. pp. 539e621. [18] Erkut E, Verter V. Hazardous materials logistics. In: Drezner Z, editor. Facility location: a survey of application and methods. New York; 1995. pp. 467e506. [19] Feillet D, Dejax P, Gendreau M, Guenguen C. An exact algorithm for the elementary shortest path problem with resource constraints: application to some vehicle routing problems. Network 2004;44(3):216e29. [20] Fleischmann B, Gnutzmann S, Sandvoß E. Dynamic vehicle routing based on online traffic information. Transp Sci 2004;38(4):420e33. [21] Gambardella LM, Taillard E, Agazzi G. MACS-VRPTW: a multiple ant colony system for vehicle routing problems with time windows. In: Corne D, Dorigo M, Glover F, editors. New ideas in optimization. London: McGraw-Hill; 1999. pp. 63e79. [22] Garcia-Martinez C, Cordon O, Herrera F. A taxonomy and an empirical analysis of multiple objective ant colony optimization algorithms for the bi-criteria TSP. Eur J Oper Res 2007;180:116e48. [23] Irnich S, Villeneuve D. The shortest-path problem with resource constraints and k-cycle elimination for k  3. Informs J Comput 2006;18(3):391e406. [24] Kuncyte R, Laberge-Nadeau C, Crainic TG, Read JA. Organisation of truckdriver training for the transportation of dangerous goods in Europe and North America. Accid Anal Prev 2003;35:191e200. [25] List GF, Mirchandani PB, Turnquist MA, Zografos KG. Modeling and analysis for hazardous materials transportation: risk analysis, routing/scheduling and facility location. Transp Sci 1991;25:100e14. [26] Lois A. On the online dial-ride-problem [PhD thesis]. Greece: University of Thessaly; 2010.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003

14

R. Pradhananga et al. / Socio-Economic Planning Sciences xxx (2014) 1e14

[27] Lozano A, Muñoz Á, Macías L, Antún JP. Hazardous materials transportation in Mexico City: chlorine and gasoline cases. Transp Res Part C 2011;19:779e89. [28] Marianov V, Revelle C. Linear, non-approximated models for optimal routing in hazardous environments. J Oper Res Soc 1998;49:157e64. [29] McMullen PR. An ant colony optimization approach to addressing a JIT sequencing problem with multiple objectives. Artif Intell Eng 2001;15(3): 309e17. [30] Miller-Hooks ED, Mahmassani HS. Optimal routing of hazardous materials in stochastic, time-varying transportation networks. Transp Res Rec 1998;1645: 143e51. [31] Nicolet-Monnier M, Gheorghe AV. In: Keller AZ, editor. Topics in safety, risk, reliability and quality: quantitative risk assessment of hazardous materials transport systems, vol. 5; 1996. Dordrecht, Netherland. [32] Nozick LK, List GF, Turnquist MA. Integrated routing and scheduling in hazardous materials transportation. Transp Sci 1997;31:200e15. [33] Pradhananga R, Hanaoka S, Sattayaprasert W. Optimization model for hazardous material transport routing in Thailand. Int J Logist Syst Manag 2011;9(1):22e42. [34] Pradhananga R, Taniguchi E, Yamada T. Ant colony system based routing and scheduling of hazardous material transportation. Procedia Soc Behav Sci 2010;2:6097e108. [35] Qureshi AG, Taniguchi E, Yamada T. An analysis of exact VRPTW solutions on ITS data-based logistics instances. Int J Intelligent Transp Syst Res 2012;10(1): 34e46. [36] Ruiz R, Moroto C, Alcaraz J. A decision support system for a real vehicle routing problem. Eur J Oper Res 2004;153:593e606. [37] Santos L, Coutinho-Rodrigues J, Current JR. Transp Res Part A 2008;42: 922e34. [38] Solomon M. Algorithms for the vehicle routing and scheduling problems with the time window constraints. Operations Res 1987;35:254e65. [39] Southworth F, Peterson B. Intermodal and international freight network modeling. Transp Res Part C 2000;8:147e66. [40] Taniguchi E, Noritake M, Yamada T, Izumitani T. Optimal size and location planning of public logistics terminals. Transp Res Part E 1999;35:207e22. [41] Tarantilis CD, Kiranoudis CT. Using the vehicle routing problem for the transportation of hazardous materials. Operational Res Int J 2001;1(1):67e78. [42] Taniguchi E, Thompson RG, Yamada T, Duin RV. City logistics: network modeling and intelligent transport systems. Oxford: Pergamon; 2001. [43] US DOT (Department of Transportation). Hazardous materials safety ; 2012. [44] US EPA (Environmental Protection Agency). ALOHA, CAMEO ; 2012.

[45] Xie C, Waller T. Optimal routing with multiple objectives: efficient algorithm and application to the hazardous materials transportation problem. Comput Aided Civ Infrastruct Eng 2012;27:77e94. [46] Zografos KG, Androutsopoulos KN. A decision support system for integrated hazardous materials routing and emergency response decisions. Transp Res Part C 2008;16:684e703. [47] Zografos KG, Androutsopoulos KN. A heuristic algorithm for solving hazardous materials distribution problems. Eur J Oper Res 2004;152(2):507e19.

Dr. Rojee Pradhananga is a post-doctoral fellow of Department of Mechanical and Industrial Engineering, Qatar University. She holds a PhD in Civil Engineering from Kyoto University, Japan and Masters in Civil Engineering from Asian Institute of Technology (AIT), Thailand. Dr. Pradhananga is an active young researcher in the field of logistics and transportation. Her research interests include but are not limited to optimization, meta-heuristics, emergency planning, city logistics and hazardous material transportation.

Professor Eiichi Taniguchi is a professor of Department of Urban Management, Kyoto University with a PhD in Civil Engineering. His research area covers urban freight transport, traffic planning and intelligent transport systems. He published more than 200 academic papers and has edited and authored several books including “City Logistics d Network Modelling and Intelligent Transport Systems”, “Innovations in Freight Transport” and “Logistics Systems for Sustainable Cities” and “Urban Transportation and Logistics: Health, Safety, and Security Concerns”.

Tadashi Yamada is an associate professor at Graduate School of Engineering, Kyoto University, Japan, with a PhD in Engineering. His research focuses mainly on Freight Transport Network Design, Supply Chain-Transport Supernetwork Modelling, Vehicle Routing and Scheduling, and Combinatorial Optimisation with Metaheuristics. His articles have been published in academic journals including Transportation Science, Transportation Research Part E, Transportation Research Record, and Transport Reviews.

Dr. Qureshi is a civil engineer by training with specialization in transportation engineering, logistics and operations research. His research interests lie in city logistics policies, such as urban consolidation centers (UCC), environmental road pricing, truck ban, etc. He is also involved in researches related with exact and heuristics optimization of different variants of vehicle routing and location problems including the hazardous material transportation.

Please cite this article in press as: Pradhananga R, et al., Bi-objective decision support system for routing and scheduling of hazardous materials, Socio-Economic Planning Sciences (2014), http://dx.doi.org/10.1016/j.seps.2014.02.003