Model and Method for Multiobjective Time-Dependent Hazardous Material Transportation

,


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 reallife traffic environment.In real life, risk on road segments is time-dependent on population density subject to time-ofday 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][3][4][5][6][7][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][19][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][22][23][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 timedependent.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-andsolve based -constraint method for the considered problem.Computational results show the effectiveness of the method.
The remainder of this paper is organized as follows.Section 2 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.

Problem Formulation
2.1.Problem Description.Consider a bidirectional graph  = (, ), where  is the set of nodes and  is the set of arcs.The graph represents the transportation network in which the vehicles carrying hazardous materials are allowed to move.Arc (, ) denotes a road segment from node  to node . kinds of hazardous materials will be transported from origin nodes  ⊂  to destination nodes  ⊂ .
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, called safety 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 Figure 1.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 timedependent because the population exposure varies with time in nature.The population exposure on each arc (, ) at time period [  ,  +1 ), denoted as   , depends on the departure time from node .Without loss of generality, set  1 = 0 as the beginning time of the first period.In a day, there are usually only several time periods [26].So travel time   on the reserved lane is far less than the length of a time period; that is,   ≪  +1 −  .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.

Problem Definition.
To model our problem, we introduce the following notation.
1 if there is a reserved lane on arc (, ) 0 otherwise

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: : Minimize subject to where  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 (, ) before lane reservation and inversely proportional to the total number of lanes on arc (, ) [13].We use the same formula in [21] to evaluate the impact   =   /(  − 1), where   is the travel time on the general lane(s) of arc (, ) and   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   = 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 , there is one and only one path departing from the source node   (resp., arriving at the destination node   ) during one and only one time period.Constraint (5) ensures the flow conservation constraint for node  in  \ {  ,   } on space and time.It represents that if shipment  arrives at a node  ( ̸ =   ,   ) via a reserved lane during time period , it must also depart from  via a reserved lane during time period .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 (, ) during any time period if no lane on the arc has been reserved.Constraints ( 7) and ( 8) mean that if shipment  passes through the reserved lane on arc (, ), then its travel time is   .Constraint (9) means that there is exactly one time period  for    on any node .In a feasible solution, shipment  passes through arc (, ) if and only if two conditions are satisfied:    > 0 and ∑  =1    = 1.Constraints ( 9) and (10) imply that departure time    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.

Solution Algorithm
In this section, an improved algorithm, called the cut-andsolve 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 Paretooptimal 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.

𝜀-Constraint
Method.Introduced by Haimes et al. [29], the -constraint method is based on a scalarization where only one objective function  (commonly, it may be the most preferred or primary one) is minimized while the others are bounded by some allowable values   ,  ∈ {1, . . ., } \ {}, 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: : Minimize   () (18) subject to where Ω is the feasible solution space.
According to [30], for a given vector of  = ( 1 , . . .,  −1 ,  +1 , . . .,   ), optimal solutions of   are weakly Paretooptimal solutions.For solving problem   , it is necessary to determine the range of   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 Section 2 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 problem  into problem  0 , which aims at minimizing objective function (1) subject to constraints (3)- (17) and Step 2. Determine the ideal point ( Step 3. Fix the values of .Determine the range of : The range is divided into  equal intervals with  + 1 points, called equidistant grid points.The value of  in constraint (20) is defined by these grid points with the following formula: Mathematical Problems in Engineering Step 4. Repeat to solve problem  0 with  and obtain  + 1 Pareto-optimal solutions.

Cut-and-Solve Method.
With the -constraint method, the considered problem is transformed into a series of singleobjective MIP problems  0 .In this section, we propose an exact and iterative approach, called cut-and-solve (CS), to solve  0 .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 th iteration of the branching tree of CS, there are only two nodes, corresponding to sparse problem (SP  ) and remaining problem (RP  ), respectively [33].The SP  can be optimally solved in a reasonable amount of computation time because it is relatively small.The optimal solution of SP  , if it exists, provides an upper bound of the IP, denoted as UB  .And if it is "good" enough, then the best upper bound of the IP found so far, denoted as UB  , is updated.Meanwhile, a lower bound of the IP, denoted as LB  , can be obtained by solving the linearly relaxed RP  .When the best upper bound found so far is smaller than or equal to LB  , a global optimal solution is obtained and the CS iteration stops.Otherwise, the current RP  is further divided into a new SP +1 and a new RP +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 Section 2,    = 1 means that shipment  passes through the reserved lane on arc (, ) when time    occurs at time period ;    = 1 means that the arriving time of shipment  at node  occurs at time period  and otherwise    = 0. We have the following property.From Property 1, constraint (23) will be added to problem  0 to obtain problem  1 , which aims at minimizing objective function (1) subject to constraints (3)-( 17), (20) and (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 SP  is too small, the optimal solution is not "good" enough to update the upper bound; if the solution space of SP  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 LB  should be provided in each iteration.For an integer programming model, LB  is usually obtained by solving a linear relaxation problem of RP  .Note that    = 1 means that shipment  passes through the reserved lane on arc (, ) during time period  and, otherwise,    = 0.If    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 LB  obtained in this case is usually very "bad." Therefore,    should be relaxed to obtain more subpaths.Similar observation can be also found for    .Another important decision variable   means whether there is a reserved lane on arc (, ).Our preliminary simulation experiments show that if   are relaxed the same as    , a considerable number of   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   are not relaxed, a better LB  can be obtained.Therefore,   are not relaxed as continous variables in this paper.Similar observation can be also found for    .In a word, in order to obtain an improved LB  , we apply the partial integrality strategy for RP  , in which only    and    are relaxed as continuous variables while the integrality of   and    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, PC  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 RP  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    may be fractional in an optimal solution of the relaxed RP  , (   ,    ,   ,    ), there may be multiple paths for some shipments.The link with the greatest value of ∑  =1    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.
Let   represent the set of shipments which contains more than one path in the fractional solution and let  V represent the first node where multipath appears for shipment V. Define   as the set of the most potential arcs for all shipments in   , which implies that the arcs in   are very likely to be selected in the final optimal solution of  1 .We have where ( V ,  V ) refers to the arc with the largest value among all the arcs outgoing from node  V .
In [32], the piercing cut is considered as a combination of some decision variables in a certain set, called   .In this paper,   refers to the set of these decision variables  V  V  V  , where V is a shipment which has multiple paths in the fractional solution and ( V ,  V ) is the most potential arcs for shipment V. Set   is defined as follows: The piercing cut (PC  ) is defined as follows: where Accordingly, the additional constraint associated with RP  can be written as follows: As mentioned above, SP  and RP  are generated by adding constraints (26) and (27) to CP  , respectively, that is, SP  : objective function (1) subject to constraints (3)-( 17), ( 20), ( 23), (26), and RP  : objective function (1) subject to constraints (3)-( 17), ( 20), ( 23), (27), and (28).
It is worth pointing out that when  = 1, the original problem  1 is considered as CP 1 .Hence, constraint (28) should be removed for RP 1 and SP 1 .
The process of the proposed -constraint is illustrated in Figure 2.

Computational Results
To evaluate the efficiency of the proposed method, 155 instances (31 sets × 5 instances) are randomly generated.The proposed algorithm is coded in .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  1 .It is allowed to run until problem  1 is solved to optimality.The transportation network (, ) 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 (, ) is generated by a probability function  (,) =  exp(−(, )/), where (, ) and  are the Euclidean distance between  and  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  is set from 1 to 5 because in real life the number of time periods is not very large [26].Let    be the accident probability of hazardous material  happening on the general lane(s) of arc (, ).Note that    >    .Our preliminary experimental results also show that the values of the input parameters    ,    ,   ,   ,   ,  int , and   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]:   = (, ) and   =   * (0.6,0.9), where  is a uniform distribution;    = (, ) *  (8,20);    =    * (0.6,0.9), whose unit is 10 −7 ;   is generated by (10, 80), whose unit is 10 4 ;  int = 10; and   is generated by  (2,5).The number of iterations of -constraint method, , is set to 20.In this paper,  = ||/|| represents the average degree of graph , where || and || are its number of nodes and arcs, respectively [35,36].The average degree of graph  is defined as its number of arcs per node, which implies the density of the graph.
Table 1 summarizes the computational results on the randomly generated instances with  = 3, || = 5, and || = 3.The total computation time for an instance represents its total running time for obtaining 21 solutions in the constraint method.Columns 1 and 2 represent the average computation time (in CPU seconds) of five instances for each set by the proposed -constraint method, in which singleobjective problem  1 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 1 and 2 are almost the same.But 2 increases with  the number of nodes more slightly than 1 in Figure 3.For example, for set 1 with || = 20, 1 is only 1.261 times as much as 2, whereas, for set 9 with || = 100, the ratio 1/2 is 2.223 times.
Table 2 summarizes the computational results on the randomly generated instances with  = 3, || = 10, and || = 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 1 and 2 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 || = 10 than for the sets with || = 5.For example, given the number of nodes and time periods, 1/2 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  0 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  0 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  0 () 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  0 () 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 || = 10, the computation time for set 16 is nearly up to 36000 s so that this set can be considered as one of the largestscale problems which can be solved in reasonable time.Table 3 shows the total computation times of the proposed method with  = 4, || = 5, and different ||.We can observe from Table 3 that, for a given number of ||, the computation time increases quickly with the value of ||.For example, the computation times for 100 nodes with || = 3 (set 30) and || = 5 (set 31) are 11.329 and 24.237 times as much as that with || = 1 (set 29), respectively.As shown in Tables 1 and 3, the computation time also increases with the average degree of network.When || = 3, given the number of nodes and shipments, the computation time with  = 4 is more than that with  = 3.For example, the computation time for set 24 is 2.645 times as much as that for set 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-andsolve 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 cutand-solve method.Finally, developing high quality heuristics is a future direction for large-scale time-dependent hazardous material transportation problems via lane reservation.
Sets and Parameters {1, 2, . . ., ||}: The set of hazardous material shipments { 1 ,  2 , . . .,  || }: The set of origin nodes { 1 ,  2 , . . .,  || }: The set of destination nodes {1, . . ., ||}: The set of time periods  int : The safety time interval between any two shipments which pass the same arc   : The travel time on the reserved lane on arc (, )   : The impact on the normal traffic due to the reserved lane on arc (, )    : The accident probability of hazardous material  on a reserved lane on arc (, )   : The population exposure along arc (, ) at time period .Decision Variables    is the arrival or departure time of shipment  at node .Note that    = 0 if shipment  does not pass node , ∀ ∈ , ∀ ∈ ,  passes the reserved lane on arc (, ) when time    occurs at time period  0 otherwise

Figure 1 :
Figure 1: Sketch of a road segment with reserved lanes.

1 Relax P 1
-solve method Determine the range of  and the S + 1 grids For each , set r := 0, UB b = +∞, and CP 1 = P with partial integrality strategy and Define set U r and PC r , and then obtain SP r and RP r Solve SP r Set UB b := UB r Relax RP r with partial integrality strategy, solve the relaxed problem, and obtain a lower bound LB r for RP r Obtain S + 1 Pareto-optimal solutions Obtain the global optimal solution UB b of P 1 UB r < UB b UB b ≤ LB r Set r := r + 1. Define CP r as RP r−1

Figure 2 :
Figure 2: Flow chart of the cut-and-solve based -constraint method.