# Model and Method for Multiobjective Time-Dependent Hazardous Material Transportation

Authors: Zhen Zhou; Ada Che; Feng Chu; Chengbin Chu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

In most of the hazardous material transportation problems, risk factors are assumed to be constant, which ignores the fact that they can vary with time throughout the day. In this paper, we deal with a novel time-dependent hazardous material transportation problem via lane reservation, in which the dynamic nature of transportation risk in the real-life traffic environment is taken into account. We first develop a multiobjective mixed integer programming (MIP) model with two conflicting objectives: minimizing the impact on the normal traffic resulting from lane reservation and minimizing the total transportation risk. We then present a cut-and-solve basedε-constraint method to solve this model. Computational results indicate that our method outperforms the ε-constraint method based on optimization software package CPLEX.


---

## Body

## 1. Introduction

Hazardous materials are a kind of goods with physical, chemical, and biological properties, which could cause lots of accidents, such as flame, explosion, and leak, in the course of production, storage, and transportation. With the development of industry, more and more hazardous substances, including raw materials, intermediate, final products, and wastes, are produced and moved daily through the transportation network in different transportation modes, such as road, rail, water, air, and pipeline. Their transportation risk imposed on environment and human is widening and deepening year by year. Hazardous material transportation has become a serious problem worldwide and attracts many researchers’ attentions in the related field.It has been pointed out in much relevant literature that the essential objective of hazardous material transportation is to minimize the transportation risk due to the nature of this problem. As we know, the selection of routes in a network for hazardous material transportation can affect its risk factors, such as the probability of hazardous material accidents and the risk exposure to the surrounding population and environment. Therefore, appropriate routing decisions are very important for hazardous material transportation management. In the last couple of decades, various applications of operations research models to hazardous material transportation have focused on risk reduction and fruitful achievements in this area have been published; please see [1] for details. However, in almost all of the hazardous material transportation problems, the transportation risk is considered to be time-invariant. That is to say, the risk of a road segment in a transportation network is assumed to be constant, which fails to capture the dynamic nature of the real-life traffic environment. In real life, risk on road segments is time-dependent on population density subject to time-of-day variation, peak and off-peak periods, various weather conditions, and so on. The time-dependent risk is one of the important features of hazardous material transportation. The time-dependent transportation problem is to decide the path for each shipment and its starting time so as to minimize the transportation risk.An important branch of the time-dependent transportation concerns hazardous materials. Time-dependent transportation problems can be distinguished into deterministic and stochastic settings. The dynamic characteristic of transportation networks usually depends on one or more traffic quantities [2], such as travel time, link volume, and queue length. For deterministic time-dependent transportation problems, part or all of the traffic quantities are assumed to be variant but all of the traffic quantities are known for a road segment. For stochastic time-dependent transportation problems, the traffic quantities are considered as random variables with time-dependent distribution functions, such as in [2–8]. For time-dependent hazardous material transportation, traffic quantities are usually travel time and transportation risk. According to optimization criteria, time-dependent transportation problems can be divided into single-objective ones and multiobjective ones.For the hazardous material transportation routing problem, a lot of works have concentrated on time-invariant risk and travel time. However, time-dependent hazardous material transportation problems have not been widely studied and only a few related publications can be found in the literature. Jia et al. [9] investigated a hazardous material transportation problem with deterministic time-dependent risk for minimizing the transportation risk. The proposed model guaranteed the minimum distance between hazardous material shipments at any time. They transformed the considered problem into a number of time-dependent shortest path problems for each truck and proposed an iterative heuristic. Erkut and Alp [10] proposed an integrated routing and scheduling problem for hazardous material transportation in a network with stochastic time-dependent accident probability, population exposure, and travel time. The model aimed to minimize transportation risk while imposing a constraint on the total travel time of the shipment. Meng et al. [11] examined a similar problem with multiple objectives. They transformed it into a time-dependent multiobjective shortest path problem subject to three types of time constraints. A dynamic programming approach was constructed to solve the problem. Both of the methods of references [10, 11] were pseudopolynomial. Nozick et al. [12] developed an approach for a hazardous material routing and scheduling problem with deterministic time-dependent risk. But their approach could not guarantee generating all Pareto-optimal paths. Chang et al. [13] proposed an effective algorithm for finding a path in stochastic time-dependent networks that could address multiple optimization criteria. In their work, travel time, transportation risk, and other traffic quantities along paths were random variables. However, the performance of the algorithm was sensitive to some parameters and the computational burden increased with the number of the dominated paths. Miller-Hooks and Mahmassani [14] proposed an optimal routing algorithm for a single hazardous material shipment in a stochastic time-dependent network, where the travel time and risk followed time-dependent normal distribution functions. They presented several procedures for determining the best compromise path to minimize the total travel time and risk (population exposure).Lane reservation is considered as a flexible and economic strategy for special events or situations, like sport games and emergencies. As stated in [15], it has been successfully applied to the Olympic Games in 2000 in Sydney [16] and in 2004 in Athens [17]. The principle of lane reservation is to reserve lanes on some road segments and/or in some specific time periods in a transportation network. Only special transportation tasks are allowed to pass through them. A probable consequence of this strategy is that it may worsen traffic congestions on general lanes, which is also an important feature of the lane reservation problem. Research works about lane reservation involving simulation tools and statistical methods can be found in the literature [18–20]. Some recent works focused on minimizing the traffic impact by linear programming models and optimal methods for various applications, such as automated truck, large sportive events, and hazardous material transportation [15, 21–24]. Our previous work [15] investigated a hazardous material transportation problem via lane reservation, in which the risk along a road segment was assumed to be constant. The goal was to minimize the transportation risk and the impact on normal traffic due to lane reservation. Remarkably, it has been shown in [15] that the lane reservation strategy can greatly reduce the transportation risk at a reasonable cost of its traffic impact. The problem proposed in this paper, in which the transportation risk is considered to be deterministic and time-dependent, is more realistic. To the best of our knowledge, this is the first work for time-dependent hazardous material transportation via lane reservation.In this paper, we investigate a novel multiobjective hazardous material transportation problem via lane reservation with a deterministic time-dependent transportation risk. As we know, the factors of the transportation risk generally include the hazardous material accident probability and the population exposure to the accidents. The accident probability estimation is influenced by the nature of roads, characteristics of the trucks, transportation environment, driver conditions [25], and so forth. Estimating the accident probability is a complicated and difficult work. For simplification, the probability of an accident is regarded to be time-invariant in this paper. Population exposure is determined by population density and area. In real life, the population density along a road segment strongly depends on time and space. As we know, the population density in hospitals, schools, factories, and so on in day time is greater than that at night, and the opposite happens in residential areas. In this work, the accident probability on reserved lanes is assumed to be constant and the population exposure is assumed to be time-dependent. Therefore, the transportation risk varies with time and space. This work is motivated by the dynamic characteristic of risk and it is a natural extension of our previous work.The contributions of this work can be summarized as follows. Firstly, we propose a novel multiobjective MIP model for the time-dependent hazardous material transportation problem via lane reservation. Secondly, we develop a cut-and-solve basedε-constraint method for the considered problem. Computational results show the effectiveness of the method.The remainder of this paper is organized as follows. Section2 describes the time-dependent hazardous material transportation problem via lane reservation and a new multiobjective model is formulated with the objectives of minimizing the total impact on the normal traffic and the total transportation risk. In Section 3, some properties are analyzed to reduce the search space of solutions and a cut-and-solve based ε-constraint method is developed. Section 4 reports the computational results of experiments. Section 5 concludes the paper and discusses future research directions.


## 2. Problem Formulation

### 2.1. Problem Description

Consider a bidirectional graphG=(V,A), where V is the set of nodes and A is the set of arcs. The graph represents the transportation network in which the vehicles carrying hazardous materials are allowed to move. Arc (i,j) denotes a road segment from node i to node j. W kinds of hazardous materials will be transported from origin nodes O⊂V to destination nodes D⊂V.Before formulation of the problem, we make the following nonrestrictive assumptions.(1) There are at least two lanes on a road segment such that one lane is allowed to be reserved. (2) From the point of view of transport safety, any path for hazardous material shipments consists of only reserved lanes. (3) Vehicles with hazardous materials travel on the reserved lanes without congestion. Consequently, travel time on reserved lanes is time-invariant. (4) Any two hazardous material shipments on the same road segment must maintain a minimum time interval, calledsafety time interval. (5) The accident probability on a reserved lane is constant and hazardous material accidents happen independently.In the network, nonreserved lanes are called general lanes, as shown in Figure1. From assumption (3), the travel time on a reserved lane is time-invariant throughout the day. Nevertheless, the risk on a road segment should be time-dependent because the population exposure varies with time in nature. The population exposure on each arc (i,j) attime period[Ik,Ik+1), denoted as Eijk, depends on the departure time from node i. Without loss of generality, set I1=0 as the beginning time of the first period. In a day, there are usually only several time periods [26]. So travel time Tij on the reserved lane is far less than the length of a time period; that is, Tij≪Ik+1-Ik. The problem is to choose lanes on the existing network to be reserved, select the path for each hazardous material shipment, and decide the travel time period for each shipment on its selected path. The objective of the problem is to seek the best trade-off for minimizing the total traffic impact on the normal traffic and the total transportation risk.Figure 1
Sketch of a road segment with reserved lanes.


### 2.2. Problem Definition

To model our problem, we introduce the following notation.Sets and Parameters
W
{
1,2
,
…
,
|
W
|
}: The set of hazardous material shipments
O
{
o
1
,
o
2
,
…
,
o
|
W
|
}: The set of origin nodes
D
{
d
1
,
d
2
,
…
,
d
|
W
|
}: The set of destination nodes
K
{
1
,
…
,
|
K
|
}: The set of time periods
T
int
⁡: The safety time interval between any two shipments which pass the same arc
T
i
j: The travel time on the reserved lane on arc (i,j)
C
i
j: The impact on the normal traffic due to the reserved lane on arc (i,j)
P
i
j
w: The accident probability of hazardous material w on a reserved lane on arc (i,j)
E
i
j
k: The population exposure along arc (i,j) at time period k.Decision Variables
t
i
w is the arrival or departure time of shipment w at node i. Note that tiw=0 if shipment w does not pass node i, ∀w∈W, ∀i∈N,
x
i
j
k
w
=
1
if
shipment
w
passes
the
reserved
lane
on
⁡
arc
i
,
j
when
time
t
i
w
occurs
at
time
period
k
0
otherwise
y
i
j
=
1
if
there
is
a
reserved
lane
on
⁡
arc
(
i
,
j
)
0
otherwise
β
i
k
w
=
1
if
time
t
i
w
occurs
at
time
period
k
0
otherwise
z
i
w
w
′
=
1
if
t
i
w
<
t
i
w
′
,
w
<
w
′
0
otherwise
.


### 2.3. Mathematical Model

The mathematical model for the hazardous material transportation problem via lane reservation in a time-dependent network is presented by constraints (3)–(17). Consider the following:P:
(1)Minimizef1=∑i,j∈ACijyij,(2)Minimizef2=∑w∈W∑i,j∈A∑k∈KPijwEijkxijkw,
subject to
(3)∑j:(ow,j)∈A∑k=1Kxowjkw=1,∀w∈W,(4)∑i:(i,dw)∈A∑k=1Kxidwkw=1,∀w∈W,(5)∑j:(i,j)∈A∑k=1Kxijkw=∑j:(j,i)∈A∑k=1Kxjikw,∀w∈W,∀i∈V∖ow,dw,(6)∑k=1Kxijkw≤yij,∀i,j∈A,∀w∈W,(7)tjw-tiw≤Tij+M1-∑k=1Kxijkw,∀i,j∈A,∀w∈W,∀i≠dw,∀j≠ow,(8)tjw-tiw≥Tij+M∑k=1Kxijkw-1,∀i,j∈A,∀w∈W,∀i≠dw,∀j≠ow,(9)∑k=1Kβikw=1,∀w∈W,∀i≠dw,(10)∑k=1KβikwIk≤tiw<∑k=1KβikwIk+1,∀w∈W,∀i≠dw,(11)tiw′-tiw≥Tint⁡∑k=1Kxijkw+∑k=1Kxijkw′-1-M1-ziww′,∀i,j∈A,∀i≠ow,ow′,dw,dw′,∀w,w′∈W,w<w′,(12)tiw-tiw′≥Tint⁡∑k=1Kxijkw+∑k=1Kxijkw′-1-Mziww′,∀i,j∈A,∀i≠ow,ow′,dw,dw′,∀w,w′∈W,w<w′,(13)tiw≥0,∀i∈N,∀w∈W,(14)yij∈0,1,∀i,j∈A,(15)xijkw∈0,1,∀i,j∈A,∀w∈W,∀k∈K,(16)ziww′∈0,1,∀i∈N,∀w,w′∈W,w<w′,(17)βikw∈0,1,∀i∈N,∀w∈W,∀k∈K,
where M is a very large positive number.Objective (1) is to minimize the total impact on the normal traffic. The impact can be considered as the increase in travel time on the general lane(s) because of the lane reservation strategy, which is proportional to the travel time on the general lane(s) of arc (i,j) before lane reservation and inversely proportional to the total number of lanes on arc (i,j) [13]. We use the same formula in [21] to evaluate the impact Cij=τij/(Mij-1), where τij is the travel time on the general lane(s) of arc (i,j) and Mij is the total number of lanes in the road segment. Statistical results from [27] revealed that the travel time on the general lanes increased by about 53% after reserving one of three lanes on A1 highway in Paris. This figure is very close to the theoretical test result (50%) obtained by the computational experiment done in [22] using the above formula to estimate the impact. If no lane is reserved, then Cij=0. Objective (2) is to minimize the total transportation risk. In this paper, the risk is measured in a traditional way, and it is measured by multiplying the probability of the hazardous material accident by its consequence. See [1] for more details. Constraint (3) (resp., (4)) means that, for shipment w, there is one and only one path departing from the source node ow (resp., arriving at the destination node dw) during one and only one time period. Constraint (5) ensures the flow conservation constraint for node i in V∖{ow,dw} on space and time. It represents that if shipment w arrives at a node i  (i≠ow,dw) via a reserved lane during time period k, it must also depart from i via a reserved lane during time period k. Constraints (3), (4), and (5) together ensure that there is one and only one path for each shipment from its origin to its corresponding destination during one and only one time period. Constraint (6) guarantees that no shipment would pass through arc (i,j) during any time period if no lane on the arc has been reserved. Constraints (7) and (8) mean that if shipment w passes through the reserved lane on arc (i,j), then its travel time is Tij. Constraint (9) means that there is exactly one time period k for tiw on any node i. In a feasible solution, shipment w passes through arc (i,j) if and only if two conditions are satisfied: tiw>0 and ∑k=1Kxijkw=1. Constraints (9) and (10) imply that departure time tiw should be located within one and only one time period and on exactly one arc. Constraints (11) and (12) guarantee that if two or more shipments pass the same reserved lane, then the safety time interval between any two shipments must be satisfied. Constraints (13)–(17) specify the restriction on the decision variables.


## 2.1. Problem Description

Consider a bidirectional graphG=(V,A), where V is the set of nodes and A is the set of arcs. The graph represents the transportation network in which the vehicles carrying hazardous materials are allowed to move. Arc (i,j) denotes a road segment from node i to node j. W kinds of hazardous materials will be transported from origin nodes O⊂V to destination nodes D⊂V.Before formulation of the problem, we make the following nonrestrictive assumptions.(1) There are at least two lanes on a road segment such that one lane is allowed to be reserved. (2) From the point of view of transport safety, any path for hazardous material shipments consists of only reserved lanes. (3) Vehicles with hazardous materials travel on the reserved lanes without congestion. Consequently, travel time on reserved lanes is time-invariant. (4) Any two hazardous material shipments on the same road segment must maintain a minimum time interval, calledsafety time interval. (5) The accident probability on a reserved lane is constant and hazardous material accidents happen independently.In the network, nonreserved lanes are called general lanes, as shown in Figure1. From assumption (3), the travel time on a reserved lane is time-invariant throughout the day. Nevertheless, the risk on a road segment should be time-dependent because the population exposure varies with time in nature. The population exposure on each arc (i,j) attime period[Ik,Ik+1), denoted as Eijk, depends on the departure time from node i. Without loss of generality, set I1=0 as the beginning time of the first period. In a day, there are usually only several time periods [26]. So travel time Tij on the reserved lane is far less than the length of a time period; that is, Tij≪Ik+1-Ik. The problem is to choose lanes on the existing network to be reserved, select the path for each hazardous material shipment, and decide the travel time period for each shipment on its selected path. The objective of the problem is to seek the best trade-off for minimizing the total traffic impact on the normal traffic and the total transportation risk.Figure 1
Sketch of a road segment with reserved lanes.


## 2.2. Problem Definition

To model our problem, we introduce the following notation.Sets and Parameters
W
{
1,2
,
…
,
|
W
|
}: The set of hazardous material shipments
O
{
o
1
,
o
2
,
…
,
o
|
W
|
}: The set of origin nodes
D
{
d
1
,
d
2
,
…
,
d
|
W
|
}: The set of destination nodes
K
{
1
,
…
,
|
K
|
}: The set of time periods
T
int
⁡: The safety time interval between any two shipments which pass the same arc
T
i
j: The travel time on the reserved lane on arc (i,j)
C
i
j: The impact on the normal traffic due to the reserved lane on arc (i,j)
P
i
j
w: The accident probability of hazardous material w on a reserved lane on arc (i,j)
E
i
j
k: The population exposure along arc (i,j) at time period k.Decision Variables
t
i
w is the arrival or departure time of shipment w at node i. Note that tiw=0 if shipment w does not pass node i, ∀w∈W, ∀i∈N,
x
i
j
k
w
=
1
if
shipment
w
passes
the
reserved
lane
on
⁡
arc
i
,
j
when
time
t
i
w
occurs
at
time
period
k
0
otherwise
y
i
j
=
1
if
there
is
a
reserved
lane
on
⁡
arc
(
i
,
j
)
0
otherwise
β
i
k
w
=
1
if
time
t
i
w
occurs
at
time
period
k
0
otherwise
z
i
w
w
′
=
1
if
t
i
w
<
t
i
w
′
,
w
<
w
′
0
otherwise
.


## 2.3. Mathematical Model

The mathematical model for the hazardous material transportation problem via lane reservation in a time-dependent network is presented by constraints (3)–(17). Consider the following:P:
(1)Minimizef1=∑i,j∈ACijyij,(2)Minimizef2=∑w∈W∑i,j∈A∑k∈KPijwEijkxijkw,
subject to
(3)∑j:(ow,j)∈A∑k=1Kxowjkw=1,∀w∈W,(4)∑i:(i,dw)∈A∑k=1Kxidwkw=1,∀w∈W,(5)∑j:(i,j)∈A∑k=1Kxijkw=∑j:(j,i)∈A∑k=1Kxjikw,∀w∈W,∀i∈V∖ow,dw,(6)∑k=1Kxijkw≤yij,∀i,j∈A,∀w∈W,(7)tjw-tiw≤Tij+M1-∑k=1Kxijkw,∀i,j∈A,∀w∈W,∀i≠dw,∀j≠ow,(8)tjw-tiw≥Tij+M∑k=1Kxijkw-1,∀i,j∈A,∀w∈W,∀i≠dw,∀j≠ow,(9)∑k=1Kβikw=1,∀w∈W,∀i≠dw,(10)∑k=1KβikwIk≤tiw<∑k=1KβikwIk+1,∀w∈W,∀i≠dw,(11)tiw′-tiw≥Tint⁡∑k=1Kxijkw+∑k=1Kxijkw′-1-M1-ziww′,∀i,j∈A,∀i≠ow,ow′,dw,dw′,∀w,w′∈W,w<w′,(12)tiw-tiw′≥Tint⁡∑k=1Kxijkw+∑k=1Kxijkw′-1-Mziww′,∀i,j∈A,∀i≠ow,ow′,dw,dw′,∀w,w′∈W,w<w′,(13)tiw≥0,∀i∈N,∀w∈W,(14)yij∈0,1,∀i,j∈A,(15)xijkw∈0,1,∀i,j∈A,∀w∈W,∀k∈K,(16)ziww′∈0,1,∀i∈N,∀w,w′∈W,w<w′,(17)βikw∈0,1,∀i∈N,∀w∈W,∀k∈K,
where M is a very large positive number.Objective (1) is to minimize the total impact on the normal traffic. The impact can be considered as the increase in travel time on the general lane(s) because of the lane reservation strategy, which is proportional to the travel time on the general lane(s) of arc (i,j) before lane reservation and inversely proportional to the total number of lanes on arc (i,j) [13]. We use the same formula in [21] to evaluate the impact Cij=τij/(Mij-1), where τij is the travel time on the general lane(s) of arc (i,j) and Mij is the total number of lanes in the road segment. Statistical results from [27] revealed that the travel time on the general lanes increased by about 53% after reserving one of three lanes on A1 highway in Paris. This figure is very close to the theoretical test result (50%) obtained by the computational experiment done in [22] using the above formula to estimate the impact. If no lane is reserved, then Cij=0. Objective (2) is to minimize the total transportation risk. In this paper, the risk is measured in a traditional way, and it is measured by multiplying the probability of the hazardous material accident by its consequence. See [1] for more details. Constraint (3) (resp., (4)) means that, for shipment w, there is one and only one path departing from the source node ow (resp., arriving at the destination node dw) during one and only one time period. Constraint (5) ensures the flow conservation constraint for node i in V∖{ow,dw} on space and time. It represents that if shipment w arrives at a node i  (i≠ow,dw) via a reserved lane during time period k, it must also depart from i via a reserved lane during time period k. Constraints (3), (4), and (5) together ensure that there is one and only one path for each shipment from its origin to its corresponding destination during one and only one time period. Constraint (6) guarantees that no shipment would pass through arc (i,j) during any time period if no lane on the arc has been reserved. Constraints (7) and (8) mean that if shipment w passes through the reserved lane on arc (i,j), then its travel time is Tij. Constraint (9) means that there is exactly one time period k for tiw on any node i. In a feasible solution, shipment w passes through arc (i,j) if and only if two conditions are satisfied: tiw>0 and ∑k=1Kxijkw=1. Constraints (9) and (10) imply that departure time tiw should be located within one and only one time period and on exactly one arc. Constraints (11) and (12) guarantee that if two or more shipments pass the same reserved lane, then the safety time interval between any two shipments must be satisfied. Constraints (13)–(17) specify the restriction on the decision variables.


## 3. Solution Algorithm

In this section, an improved algorithm, called the cut-and-solve basedε-constraint method, is developed for solving the multiobjective model. There are several techniques to solve a multiobjective problem in the literature, such as the weighted-sum method, the ε-constraint method, the goal attainment approach, and metaheuristics [28]. With the ε-constraint method, a multiobjective problem can be converted into a series of single-objective problems to obtain Pareto-optimal solutions. Moreover, this method can alleviate the difficulties in setting up an appropriate weight vector faced by the weighted sum method and the goal attainment approach. Therefore, in this paper, the ε-constraint method is used to solve the considered problem.


### 3.1.ε-Constraint Method

Introduced by Haimes et al. [29], the ε-constraint method is based on a scalarization where only one objective function n (commonly, it may be the most preferred or primary one) is minimized while the others are bounded by some allowable values εi, i∈1,…,m∖{n}, and added to the original model as constraints. The original multiobjective optimization problem is converted into a series of single-objective optimization problems as given below:P
ε:
(18)Minimizefnx
subject to
(19)fix≤εi,∀i∈1,…,m∖nx∈Ω,
where Ω is the feasible solution space.According to [30], for a given vector of ε=(ε1,…,εn-1,εn+1,…,εm), optimal solutions of Pε are weakly Pareto-optimal solutions. For solving problem Pε, it is necessary to determine the range of εi that can be defined by ideal and nadir points of the original problem [31]. The ideal and nadir points can be obtained from the solutions of subproblems using the procedure described below.The proposed mathematical model in Section2 has conflicting objective functions (1) and (2) subject to constraints (3)–(17). From assumption (2), any path for hazardous material shipments contains only exclusive reserved lanes. Therefore, in this paper, the objective of minimizing the traffic impact is considered as the main objective function. In this way, the original multiobjective model can be converted into a single-objective one, in which the traffic impact is minimized, while the transportation risk subject to ε forms a new constraint and is added to the original model.The procedure of theε-constraint method is described as follows [15, 29].Step 1.
Transform problemP into problem P0, which aims at minimizing objective function (1) subject to constraints (3)–(17) and
(20)∑w∈W∑(i,j)∈A∑k∈KPijwEijkδijkw≤ε.Step 2.
Determine the ideal point (f1I,f2I) and the nadir point (f1N,f2N).
The ideal point corresponds to lower limits offl  (l=1,2), where flI=min⁡fl subject to constraints (3)–(17). Thenadir point indicates upper limits of f1 and f2, respectively, where f1N=min⁡f1 subject to constraints (3)–(17) and f2=f2I. And f2N=min⁡f2 subject to constraints (3)–(17) and f1=f1I. That is, the values of the objective function f2 are bounded by [f2I,f2N].Step 3.
Fix the values ofε.
Determine the range ofε:
(21)range2=f2N-f2I.
The range is divided intoS equal intervals with S+1 points, calledequidistant grid points. The value of ε in constraint (20) is defined by these grid points with the following formula:
(22)εs=f2N-range2S*s,s=0,1,…,S.Step 4.
Repeat to solve problemP0 with ε and obtain S+1 Pareto-optimal solutions.


### 3.2. Cut-and-Solve Method

With theε-constraint method, the considered problem is transformed into a series of single-objective MIP problems P0. In this section, we propose an exact and iterative approach, called cut-and-solve (CS), to solve P0. The CS method was first presented by Climer and Zhang and was proved to be efficient for the traveling salesman problem [32]. It is a special branch-and-bound search strategy but can avoid making wrong choices in depth-first branch-and-bound. Given an integer programming problem (IP) with the minimization of an objective function, at the rth iteration of the branching tree of CS, there are only two nodes, corresponding tosparse problem (SPr) andremaining problem (RPr), respectively [33]. The SPr can be optimally solved in a reasonable amount of computation time because it is relatively small. The optimal solution of SPr, if it exists, provides an upper bound of the IP, denoted as UBr. And if it is “good” enough, then the best upper bound of the IP found so far, denoted as UBb, is updated. Meanwhile, a lower bound of the IP, denoted as LBr, can be obtained by solving the linearly relaxed RPr. When the best upper bound found so far is smaller than or equal to LBr, a global optimal solution is obtained and the CS iteration stops. Otherwise, the current RPr is further divided into a new SPr+1 and a new RPr+1 by adding a piercing cut for the next iteration. The process continues until an optimal solution to the original IP is obtained.(
1) Preprocessing. Before using the cut-and-solve method, we analyze properties of the model in order to reduce the search space. If some variables can be fixed in advance, the reduction in the number of decision variables will possibly speed up the CS process.As defined in Section2, xijkw=1 means that shipment w passes through the reserved lane on arc (i,j) when time tiw occurs at time period k; βikw=1 means that the arriving time of shipment w at node i occurs at time period k and otherwise βikw=0. We have the following property.Property 1.
Ifxijkw=1, then βjkw+βj(k+1)w=1, ∀k∈K, w∈W, ∀(i,j)∈A, i≠dw.Proof.
Note that ifxijkw=1, then βikw=1. It implies that tiw∈[Ik,Ik+1). That is, Ik≤tiw<Ik+1. For arc (i,j), we can deduce that Ik+Tij≤tjw=tiw+Tij<Ik+1+Tij. Note that Tij≪Ik+1-Ik. We thus have Ik≤Ik+Tij≤tjw<Ik+1+Tij≤Ik+2, which implies two cases. In one case, Ik≤Ik+Tij≤tjw<Ik+1, and, in the other case, Ik+1≤tjw<Ik+1+Tij≤Ik+2. That is to say, either tjw∈[Ik,Ik+1) or tjw∈[Ik+1,Ik+2). It follows that βjkw+βj(k+1)w=1.From Property1, constraint (23) will be added to problem P0 to obtain problem P1, which aims at minimizing objective function (1) subject to constraints (3)–(17), (20) and
(23)βikw+βik+1w≥xijkw,∀i,j∈A,∀k∈K,∀w∈W.(
2) Definition of Piercing Cut, Remaining Problem, and Sparse Problem. A key issue for the cut-and-solve method is to find an appropriate piercing cut that separates the current remaining problem into a new sparse problem and a new remaining problem. The reasons are the following. If the solution space of SPr is too small, the optimal solution is not “good” enough to update the upper bound; if the solution space of SPr is too large, it will take too much time to obtain an optimal solution. Another important factor which influences the efficiency of the cut-and-solve method is that tight LBr should be provided in each iteration. For an integer programming model, LBr is usually obtained by solving a linear relaxation problem of RPr. Note that xijkw=1 means that shipment w passes through the reserved lane on arc (i,j) during time period k and, otherwise, xijkw=0. If xijkw are not relaxed, our preliminary simulation experiments show that most of them take the values of zero, which means that very few reserved lanes are passed by the shipments at any time period. The LBr obtained in this case is usually very “bad.” Therefore, xijkw should be relaxed to obtain more subpaths. Similar observation can be also found for βikw. Another important decision variable yij means whether there is a reserved lane on arc (i,j). Our preliminary simulation experiments show that if yij are relaxed the same as xijkw, a considerable number of yij take the values greater than zero. It means that there may exist a reserved lane on the corresponding arc. However, the experiments also show that if yij are not relaxed, a better LBr can be obtained. Therefore, yij are not relaxed as continous variables in this paper. Similar observation can be also found for ziww′. In a word, in order to obtain an improved LBr, we apply the partial integrality strategy for RPr, in which only xijkw and βikw are relaxed as continuous variables while the integrality of yij and ziww′ is maintained.Climer and Zhang introduced a general procedure for generating piercing cuts based on reduced cost from an optimal solution of linear relaxation problem [32]. In their work, PCr was defined as a set including the decision variables with large reduced cost. But the general procedure is not appropriate to MIP because it has been shown by our preliminary experimental results that the lower bound of the proposed problem obtained by the linearly relaxed RPr is not good enough and the reduced costs of decision variables are often missing. Reference [24] proposed a new piercing cut technique for MIP using the number of “critical links.” Since the value of xijkw may be fractional in an optimal solution of the relaxed RPr, (t¯iw,x¯ijkw,y¯ij,z¯iww′), there may be multiple paths for some shipments. The link with the greatest value of ∑k=1Kxijkw is called “critical link.” For more details of the piercing cut based on “critical link,” please see [24]. The considered problem in this paper is different from that in [24]. The piercing cut in [24] cannot be used directly for our problem. A new piercing cut based on “critical link” is proposed as follows.LetLr represent the set of shipments which contains more than one path in the fractional solution and let iv represent the first node where multipath appears for shipment v. Define ar as the set of the most potential arcs for all shipments in Lr, which implies that the arcs in ar are very likely to be selected in the final optimal solution of P1. We have
(24)ar=i*,j*∣i*,j*=argmax⁡iv,jv∈A∑k=1Kx¯ivjvkv,v∈Lr,
where (iv,jv) refers to the arc with the largest value among all the arcs outgoing from node iv.In [32], the piercing cut is considered as a combination of some decision variables in a certain set, called Ur. In this paper, Ur refers to the set of these decision variables xivjvkv, where v is a shipment which has multiple paths in the fractional solution and (iv,jv) is the most potential arcs for shipment v. Set Ur is defined as follows:
(25)Ur=xivjvkv∣iv,jv∈ar.The piercing cut (PCr) is defined as follows:
(26)∑xivjvkv∈Ur∑k=1Kxivjvkv≥hr,
where hr is a given integer in [1,Lr].Accordingly, the additional constraint associated withRPr can be written as follows:
(27)∑xivjvkv∈Ur∑k=1Kxivjvkv≤hr-1.As mentioned above,SPr and RPr are generated by adding constraints (26) and (27) to CPr, respectively, that is,SPr: objective function (1) subject to constraints (3)–(17), (20), (23), (26), and
(28)∑xivjvkv∈Ul∑k=1Kxivjvkv≤hl-1,l=1,…,r-1,RPr: objective function (1) subject to constraints (3)–(17), (20), (23), (27), and (28).It is worth pointing out that whenr=1, the original problem P1 is considered as CP1. Hence, constraint (28) should be removed for RP1 and SP1.The process of the proposedε-constraint is illustrated in Figure 2.Figure 2
Flow chart of the cut-and-solve basedε-constraint method.


## 3.1.ε-Constraint Method

Introduced by Haimes et al. [29], the ε-constraint method is based on a scalarization where only one objective function n (commonly, it may be the most preferred or primary one) is minimized while the others are bounded by some allowable values εi, i∈1,…,m∖{n}, and added to the original model as constraints. The original multiobjective optimization problem is converted into a series of single-objective optimization problems as given below:P
ε:
(18)Minimizefnx
subject to
(19)fix≤εi,∀i∈1,…,m∖nx∈Ω,
where Ω is the feasible solution space.According to [30], for a given vector of ε=(ε1,…,εn-1,εn+1,…,εm), optimal solutions of Pε are weakly Pareto-optimal solutions. For solving problem Pε, it is necessary to determine the range of εi that can be defined by ideal and nadir points of the original problem [31]. The ideal and nadir points can be obtained from the solutions of subproblems using the procedure described below.The proposed mathematical model in Section2 has conflicting objective functions (1) and (2) subject to constraints (3)–(17). From assumption (2), any path for hazardous material shipments contains only exclusive reserved lanes. Therefore, in this paper, the objective of minimizing the traffic impact is considered as the main objective function. In this way, the original multiobjective model can be converted into a single-objective one, in which the traffic impact is minimized, while the transportation risk subject to ε forms a new constraint and is added to the original model.The procedure of theε-constraint method is described as follows [15, 29].Step 1.
Transform problemP into problem P0, which aims at minimizing objective function (1) subject to constraints (3)–(17) and
(20)∑w∈W∑(i,j)∈A∑k∈KPijwEijkδijkw≤ε.Step 2.
Determine the ideal point (f1I,f2I) and the nadir point (f1N,f2N).
The ideal point corresponds to lower limits offl  (l=1,2), where flI=min⁡fl subject to constraints (3)–(17). Thenadir point indicates upper limits of f1 and f2, respectively, where f1N=min⁡f1 subject to constraints (3)–(17) and f2=f2I. And f2N=min⁡f2 subject to constraints (3)–(17) and f1=f1I. That is, the values of the objective function f2 are bounded by [f2I,f2N].Step 3.
Fix the values ofε.
Determine the range ofε:
(21)range2=f2N-f2I.
The range is divided intoS equal intervals with S+1 points, calledequidistant grid points. The value of ε in constraint (20) is defined by these grid points with the following formula:
(22)εs=f2N-range2S*s,s=0,1,…,S.Step 4.
Repeat to solve problemP0 with ε and obtain S+1 Pareto-optimal solutions.


## 3.2. Cut-and-Solve Method

With theε-constraint method, the considered problem is transformed into a series of single-objective MIP problems P0. In this section, we propose an exact and iterative approach, called cut-and-solve (CS), to solve P0. The CS method was first presented by Climer and Zhang and was proved to be efficient for the traveling salesman problem [32]. It is a special branch-and-bound search strategy but can avoid making wrong choices in depth-first branch-and-bound. Given an integer programming problem (IP) with the minimization of an objective function, at the rth iteration of the branching tree of CS, there are only two nodes, corresponding tosparse problem (SPr) andremaining problem (RPr), respectively [33]. The SPr can be optimally solved in a reasonable amount of computation time because it is relatively small. The optimal solution of SPr, if it exists, provides an upper bound of the IP, denoted as UBr. And if it is “good” enough, then the best upper bound of the IP found so far, denoted as UBb, is updated. Meanwhile, a lower bound of the IP, denoted as LBr, can be obtained by solving the linearly relaxed RPr. When the best upper bound found so far is smaller than or equal to LBr, a global optimal solution is obtained and the CS iteration stops. Otherwise, the current RPr is further divided into a new SPr+1 and a new RPr+1 by adding a piercing cut for the next iteration. The process continues until an optimal solution to the original IP is obtained.(
1) Preprocessing. Before using the cut-and-solve method, we analyze properties of the model in order to reduce the search space. If some variables can be fixed in advance, the reduction in the number of decision variables will possibly speed up the CS process.As defined in Section2, xijkw=1 means that shipment w passes through the reserved lane on arc (i,j) when time tiw occurs at time period k; βikw=1 means that the arriving time of shipment w at node i occurs at time period k and otherwise βikw=0. We have the following property.Property 1.
Ifxijkw=1, then βjkw+βj(k+1)w=1, ∀k∈K, w∈W, ∀(i,j)∈A, i≠dw.Proof.
Note that ifxijkw=1, then βikw=1. It implies that tiw∈[Ik,Ik+1). That is, Ik≤tiw<Ik+1. For arc (i,j), we can deduce that Ik+Tij≤tjw=tiw+Tij<Ik+1+Tij. Note that Tij≪Ik+1-Ik. We thus have Ik≤Ik+Tij≤tjw<Ik+1+Tij≤Ik+2, which implies two cases. In one case, Ik≤Ik+Tij≤tjw<Ik+1, and, in the other case, Ik+1≤tjw<Ik+1+Tij≤Ik+2. That is to say, either tjw∈[Ik,Ik+1) or tjw∈[Ik+1,Ik+2). It follows that βjkw+βj(k+1)w=1.From Property1, constraint (23) will be added to problem P0 to obtain problem P1, which aims at minimizing objective function (1) subject to constraints (3)–(17), (20) and
(23)βikw+βik+1w≥xijkw,∀i,j∈A,∀k∈K,∀w∈W.(
2) Definition of Piercing Cut, Remaining Problem, and Sparse Problem. A key issue for the cut-and-solve method is to find an appropriate piercing cut that separates the current remaining problem into a new sparse problem and a new remaining problem. The reasons are the following. If the solution space of SPr is too small, the optimal solution is not “good” enough to update the upper bound; if the solution space of SPr is too large, it will take too much time to obtain an optimal solution. Another important factor which influences the efficiency of the cut-and-solve method is that tight LBr should be provided in each iteration. For an integer programming model, LBr is usually obtained by solving a linear relaxation problem of RPr. Note that xijkw=1 means that shipment w passes through the reserved lane on arc (i,j) during time period k and, otherwise, xijkw=0. If xijkw are not relaxed, our preliminary simulation experiments show that most of them take the values of zero, which means that very few reserved lanes are passed by the shipments at any time period. The LBr obtained in this case is usually very “bad.” Therefore, xijkw should be relaxed to obtain more subpaths. Similar observation can be also found for βikw. Another important decision variable yij means whether there is a reserved lane on arc (i,j). Our preliminary simulation experiments show that if yij are relaxed the same as xijkw, a considerable number of yij take the values greater than zero. It means that there may exist a reserved lane on the corresponding arc. However, the experiments also show that if yij are not relaxed, a better LBr can be obtained. Therefore, yij are not relaxed as continous variables in this paper. Similar observation can be also found for ziww′. In a word, in order to obtain an improved LBr, we apply the partial integrality strategy for RPr, in which only xijkw and βikw are relaxed as continuous variables while the integrality of yij and ziww′ is maintained.Climer and Zhang introduced a general procedure for generating piercing cuts based on reduced cost from an optimal solution of linear relaxation problem [32]. In their work, PCr was defined as a set including the decision variables with large reduced cost. But the general procedure is not appropriate to MIP because it has been shown by our preliminary experimental results that the lower bound of the proposed problem obtained by the linearly relaxed RPr is not good enough and the reduced costs of decision variables are often missing. Reference [24] proposed a new piercing cut technique for MIP using the number of “critical links.” Since the value of xijkw may be fractional in an optimal solution of the relaxed RPr, (t¯iw,x¯ijkw,y¯ij,z¯iww′), there may be multiple paths for some shipments. The link with the greatest value of ∑k=1Kxijkw is called “critical link.” For more details of the piercing cut based on “critical link,” please see [24]. The considered problem in this paper is different from that in [24]. The piercing cut in [24] cannot be used directly for our problem. A new piercing cut based on “critical link” is proposed as follows.LetLr represent the set of shipments which contains more than one path in the fractional solution and let iv represent the first node where multipath appears for shipment v. Define ar as the set of the most potential arcs for all shipments in Lr, which implies that the arcs in ar are very likely to be selected in the final optimal solution of P1. We have
(24)ar=i*,j*∣i*,j*=argmax⁡iv,jv∈A∑k=1Kx¯ivjvkv,v∈Lr,
where (iv,jv) refers to the arc with the largest value among all the arcs outgoing from node iv.In [32], the piercing cut is considered as a combination of some decision variables in a certain set, called Ur. In this paper, Ur refers to the set of these decision variables xivjvkv, where v is a shipment which has multiple paths in the fractional solution and (iv,jv) is the most potential arcs for shipment v. Set Ur is defined as follows:
(25)Ur=xivjvkv∣iv,jv∈ar.The piercing cut (PCr) is defined as follows:
(26)∑xivjvkv∈Ur∑k=1Kxivjvkv≥hr,
where hr is a given integer in [1,Lr].Accordingly, the additional constraint associated withRPr can be written as follows:
(27)∑xivjvkv∈Ur∑k=1Kxivjvkv≤hr-1.As mentioned above,SPr and RPr are generated by adding constraints (26) and (27) to CPr, respectively, that is,SPr: objective function (1) subject to constraints (3)–(17), (20), (23), (26), and
(28)∑xivjvkv∈Ul∑k=1Kxivjvkv≤hl-1,l=1,…,r-1,RPr: objective function (1) subject to constraints (3)–(17), (20), (23), (27), and (28).It is worth pointing out that whenr=1, the original problem P1 is considered as CP1. Hence, constraint (28) should be removed for RP1 and SP1.The process of the proposedε-constraint is illustrated in Figure 2.Figure 2
Flow chart of the cut-and-solve basedε-constraint method.


## 4. Computational Results

To evaluate the efficiency of the proposed method, 155 instances (31 sets × 5 instances) are randomly generated. The proposed algorithm is coded inC. The computational experiments are carried out on an HP PC with 3.10 GHz Intel Core processor and 4 GB RAM under Windows 7 environment. The CPLEX MIP solver (version 12.5) under default settings is used to solve P1. It is allowed to run until problem P1 is solved to optimality.The transportation networkG(V,A) in this work is generated according to the random network topology generator proposed by Waxman [34]. The nodes are randomly and uniformly generated in the plane [0,100]×[0,100], while arc (i,j) is generated by a probability function p(i,j)=βexp⁡⁡(-di,j/αL), where d(i,j) and L are the Euclidean distance between i and j and the maximum distance between any two nodes in the graph, respectively, and 0<α,β≤1. The origin and destination pairs are randomly selected from the set of nodes. Parameter K is set from 1 to 5 because in real life the number of time periods is not very large [26]. Let πijw be the accident probability of hazardous material w happening on the general lane(s) of arc (i,j). Note that πijw>Pijw. Our preliminary experimental results also show that the values of the input parameters πijw, Pijw, Tij, τij, Eijk, Tint⁡, and Mij have little effect on the performance of the proposed algorithm for a considerable number of instances and they are generated according to the following ways. The data of parameters are set with the similar way of [15]: τij=d(i,j) and Tij=τij*U(0.6,0.9), where U is a uniform distribution; πijw=d(i,j)*U(8,20); Pijw=πijw*U(0.6,0.9), whose unit is 10−7; Eijk is generated by U(10,80), whose unit is 104; Tint⁡=10; and Mij is generated by U(2,5). The number of iterations of ε-constraint method, S, is set to 20.In this paper,N=|A|/|V| represents the average degree of graph G, where |V| and |A| are its number of nodes and arcs, respectively [35, 36]. The average degree of graph G is defined as its number of arcs per node, which implies the density of the graph.Table1 summarizes the computational results on the randomly generated instances with N=3, |W|=5, and |K|=3. The total computation time for an instance represents its total running time for obtaining 21 solutions in the ε-constraint method. Columns T1 and T2 represent the average computation time (in CPU seconds) of five instances for each set by the proposed ε-constraint method, in which single-objective problem P1 is solved by the optimization software package CPLEX and the CS based method, respectively. We can observe from Table 1 that the proposed CS based method is more efficient than CPLEX and both of the total computation times moderately increase with the number of nodes. It is worth pointing out that the trends of two curves of T1 and T2 are almost the same. But T2 increases with the number of nodes more slightly than T1 in Figure 3. For example, for set 1 with |V|=20, T1 is only 1.261 times as much as T2, whereas, for set 9 with |V|=100, the ratio T1/T2 is 2.223 times.Table 1
Computational results withN=3, |W|=5, and |K|=3.
Set
|
V
|
T
1
T
2
T
1
/
T
2
1
20
42.492
33.689
1.261
2
30
159.133
101.038
1.575
3
40
229.150
155.282
1.476
4
50
328.191
174.990
1.875
5
60
319.409
200.318
1.593
6
70
865.940
252.524
3.429
7
80
1375.141
465.900
2.956
8
90
2089.713
830.784
2.515
9
100
2353.960
1058.735
2.223Figure 3
Comparison ofT1 and T2.Table2 summarizes the computational results on the randomly generated instances with N=3, |W|=10, and |K|=3. It can be observed from Table 2 that the computation times of CPLEX and the proposed method drastically increase with the number of nodes, but the latter increases more slightly than the former. Take sets 10 and 11 for example; given the number of shipments and time periods, the computation times T1 and T2 for set 11 are 6.636 and 4.440 times as much as those for set 10, respectively. From Tables 1 and 2, we can also observe that the CS based method is more efficient for the sets with |W|=10 than for the sets with |W|=5. For example, given the number of nodes and time periods, T1/T2 for set 3 is only 1.476, while it increases to 3.573 for set 12. When the number of nodes increases to 60 in Table 2, CPLEX cannot find 21 optimal solutions of P0 within 36000 s, but the proposed method can do it well for sets with up to 80 nodes. We set a threshold for the computation time for each problem P0 to 36000 s/21. When the threshold is reached and a problem is not solved to optimality, CPLEX is terminated. For sets 14–16, although CPLEX fails to find all optimal solutions, it can provide lower bounds and upper bounds of P0(s) which are not solved to optimality when CPLEX is terminated. The gap between the lower bound and upper bound is denoted by (upper bound − lower bound)/upper bound and it implies the extent of optimality. GAP in Table 2 denotes the average gap for problems P0(s) that are not solved to optimality. The average gaps of sets 14–16 are 5.309%, 3.217%, and 3.949%. The gaps are relatively small, which means that the obtained solutions might be relatively close to optimality for sets 14–16. In Table 2, given that |W|=10, the computation time for set 16 is nearly up to 36000 s so that this set can be considered as one of the largest-scale problems which can be solved in reasonable time.Table 2
Computational results withN=3, |W|=10, and |K|=3.
Set
V
T
1
T
2
T
1
/
T
2
GAP
10
20
529.566
385.258
1.375
0
11
30
3514.437
1710.633
2.055
0
12
40
11489.591
3215.313
3.573
0
13
50
28094.791
6698.321
4.194
0
14
60
>36000.000
5703.657
—
5.309%
15
70
>36000.000
7568.321
—
3.217%
16
80
>36000.000
35434.740
—
3.949%Table3 shows the total computation times of the proposed method with N=4, |W|=5, and different |K|. We can observe from Table 3 that, for a given number of |V|, the computation time increases quickly with the value of |K|. For example, the computation times for 100 nodes with |K|=3 (set 30) and |K|=5 (set 31) are 11.329 and 24.237 times as much as that with |K|=1 (set 29), respectively. As shown in Tables 1 and 3, the computation time also increases with the average degree of network. When |K|=3, given the number of nodes and shipments, the computation time with N=4 is more than that with N=3. For example, the computation time for set 24 is 2.645 times as much as that for set 5.Table 3
Computational results withN=4, |W|=5, and |K| = 1–5.
Set
|
V
|
K
T
1
T
2
T
1
/
T
2
17
20
1
26.923
23.804
1.131
18
3
63.177
41.690
1.515
19
5
160.727
141.198
1.138
20
40
1
45.259
31.647
1.430
21
3
202.990
165.206
1.229
22
5
311.470
251.380
1.239
23
60
1
53.883
37.412
1.440
24
3
625.306
529.896
1.180
25
5
1289.057
1141.264
1.129
26
80
1
83.724
51.417
1.628
27
3
2920.514
831.624
3.512
28
5
4002.699
1575.911
2.540
29
100
1
149.272
101.348
1.473
30
3
3885.786
1148.121
3.384
31
5
9610.322
2456.343
3.912


## 5. Conclusion

This paper investigated the time-dependent lane reservation problem for hazardous material transportation, in which the transportation risk varies with time throughout the day. A new multiobjective mixed integer programming model was first proposed for the considered problem. Then a cut-and-solve method basedε-constraint method was developed. For improving the cut-and-solve method, a specific property of the considered problem was explored. Computational results showed that our method outperforms the optimization software package CPLEX.For the considered problem, one of the future research directions is to study properties of the model and attempt to reduce the search space of solution. Exploiting efficient piercing cuts is another way to improve the efficiency of the cut-and-solve method. Finally, developing high quality heuristics is a future direction for large-scale time-dependent hazardous material transportation problems via lane reservation.



---

*Source: 832814-2014-12-14.xml*