- Email: [email protected]
Freight transportation in railway networks with automated terminals: A mathematical model and MIP heuristic approaches
European Journal of Operational Research 214 (2011) 588–594
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Production, Manufacturing and Logistics
Freight transportation in railway networks with automated terminals: A mathematical model and MIP heuristic approaches D. Anghinolfi, M. Paolucci, S. Sacone, S. Siri ⇑ Department of Communications, Computer and System Sciences, University of Genova, Italy
a r t i c l e
i n f o
Article history: Received 10 December 2010 Accepted 11 May 2011 Available online 17 May 2011 Keywords: Freight transportation Optimal planning Mathematical programming MIP heuristics
a b s t r a c t In this paper we propose a planning procedure for serving freight transportation requests in a railway network with fast transfer equipment at terminals. We consider a transportation system where different customers make their requests (orders) for moving boxes, i.e., either containers or swap bodies, between different origins and destinations, with specific requirements on delivery times. The decisions to be taken concern the route (and the corresponding sequence of trains) that each box follows in the network and the assignment of boxes to train wagons, taking into account that boxes can change more than one train and that train timetables are fixed. The planning procedure includes a pre-analysis step to determine all the possible sequences of trains for serving each order, followed by the solution of a 0–1 linear programming problem to find the optimal assignment of each box to a train sequence and to a specific wagon for each train in the sequence. This latter is a generalized assignment problem which is NP-hard. Hence, in order to find good solutions in acceptable computation times, two MIP heuristic approaches are proposed and tested through an experimental analysis considering realistic problem instances. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction Freight transportation is a very important activity in our society and it is also a complex domain to be studied. Freight transportation planning has been often faced by mathematical programming and optimization approaches and dealt with in many research works (see for instance the reviews [1,2]). In this paper, a kind of long haul transportation of containerized goods [3] is considered, characterized by a railway network whose terminals are automated, i.e., equipped with transhipment devices allowing fast train–train transfer operations. The study of this transportation service derives from a research project carried out with an Italian company that is developing a new automated container handling system named Metrocargo [4]. The adoption of Metrocargo in the national railway network would provide a very innovative transportation service for the Italian logistic context. Our main interest in this project is the definition of optimization methods for the relevant off-line planning. We consider the point of view of a central decision maker who provides the transportation service in the railway network and has to satisfy transportation requests (orders) from customers. The timetable and stop sequence of trains in the railway network are fixed by the provider of the transportation service and are not ⇑ Corresponding author. Tel.: +39 019 21945137; fax: +39 019 21945212. E-mail address: [email protected] (S. Siri). 0377-2217/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2011.05.013
matter of decision. The orders coming from customers regard the movement of a set of boxes (containers or swap bodies) from an origin to a destination within a deadline. We aim at generating transportation plans for orders considering the boxes similarly to passengers planning a train trip: for each box of each order we determine the sequence of trains and train changes at terminals so that the box can reach its destination in time. The planning problem faced in this paper combines two main decisions. The former concerns the routing of boxes in the network and the selection of train changes at terminals; the latter decision deals with the assignment of boxes to wagons of the selected trains. These two decisions have been generally considered separately in the literature. The box routing and train selection are often related to aggregate flows of goods represented as a set of commodities, as in the general formulation provided in the review [5]. In this paper, instead, we determine the routes and the sequences of trains so that the boxes of a same order can reach their common destination independently from each other, i.e., even through different routes and with different sequences of trains. The decision on wagon assignment has already been treated in the literature usually considering a single terminal. In [6] the authors assume a rail–rail container terminal with a rapid transshipment yard and they define different models and algorithms for determining the optimal placement of containers on arriving trains and on departing trains, in order to minimize the movements of containers in the terminal. Analogously, in [7] the
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
so-called load planning problem (assignment of containers to train slots) is faced for an intermodal terminal where containers are transferred between trucks and trains. Differently, our purpose in this work is optimally matching the whole transportation demand for the railway network with the overall available transportation offer (the fixed train timetables and stop sequences). Specifically, we determine the minimum cost assignment of the boxes to wagons of the different trains in order to satisfy train and wagon capacity constraints. Since the train timetables and stop sequences are fixed input data, this paper does not deal with train routing and scheduling, as done for instance in [8,9], or with other classical railway problems, such as the definition of railroad blocking plans [10]. Even though box routing and train selection, as well as wagon assignment, have been already treated in the literature, the main novelty of this paper stands in considering them jointly. As detailed in the following, the considered planning problem (a preliminary version can be found in [11,12]) basically consists in a generalized assignment problem (which is known to be NP-hard [13]) that we formulate as a 0–1 linear programming problem. This kind of problems can be faced by means of general purpose mixed integer programming (MIP) solvers that seek for the optimal solution by exploiting branch-and-bound or branch-and-cut algorithms. However, most MIP problems are very difficult to solve when their dimensions become too large. For this reason, complex combinatorial optimization problems from both academic research and real world applications have been tackled by specialized heuristics or metaheuristics. Recently, new approaches have been proposed to combine ideas from metaheuristics with MIP solver algorithms [14–16], and they are generally referred to as matheuristics [17]. In many cases, matheuristics have been adopted for large scale problems coming from real applications in order to reduce the solution search space and to speed up the computation (e.g., in production planning [18], in production scheduling [19], in supply chain management [20] and in vehicle routing and transportation logistics [21]). In this paper, we propose two matheuristic approaches for the considered planning problem; the first one is an ad hoc procedure integrated with the MIP solver to drive and restrict the solution space exploration, whereas the second one is a general purpose heuristic method to solve MIP problems by iteratively calling the MIP solver as a black box tool. This latter approach shares some concepts with the MIP heuristics recently appeared in literature as Local Branching [14], Relaxation Induced Neighborhood Search [15] and Variable Neighborhood Decomposition Search [22], but, differently from them, it exploits only a randomization mechanism to guide the MIP solver. This paper is organized as follows. In Section 2 the main aspects of the considered problem are introduced; in Section 3 the preanalysis and the mathematical programming formulation adopted in the planning procedure are described in detail. Then, in Section 4 the two proposed MIP heuristic approaches are introduced. The experimental results are discussed in Section 5 and some conclusive remarks are finally reported in Section 6. 2. Problem description The Metrocargo system is a new technology, now at a prototype level, allowing to load and unload cargo units from trains in a horizontal way and under the electric feeding line, thanks to some shuttles and a properly devised storage area where all the operations are realized automatically. With this technology rail–rail transfer operations at terminals become fast. The planning procedure proposed in this paper refers to a railway network with automated terminals equipped with Metrocargo system or any other innovative system allowing fast transfers. With such systems, the boxes can be considered as passengers that change trains in their
589
route to destination by using the available (i.e., scheduled) train services. The transportation demand to be satisfied is given by a set of orders. For each order we consider the origin and destination, the release time (i.e., the earliest time from which the boxes are available at the origin), the deadline (i.e., the latest time by which the boxes must be delivered at destination), and the number of boxes, each one characterized by length and weight. The railway network is composed of railway terminals and railway links connecting the terminals. We assume that for each pair of terminals in the network, i.e., for each origin–destination pair, the available routes (defined as sequences of railway links) are specified a priori. Such routes are identified by the transportation service provider among all the possible link sequences connecting an origin–destination pair, as the most effective ones to provide the transportation service. Note that this is usually done in order to eliminate possible but not reasonable/acceptable routes. Moreover, the available routes connecting each origin–destination pair are a priori ordered by the transportation service provider according to a given priority. The goal of the proposed planning procedure is to make the planning of transportation operations for all the orders within a specified time horizon. The problem is to determine for each box of each order (a) the route that it must follow from its origin to its destination, (b) the sequence of trains that it must use along such route so that the order time constraints are satisfied, and (c) the wagons used to transport it for each train in the sequence. The problem data are listed in the following.
N R RN n dn qn hn
L
nlr Sr
kr gr Wr f r;w gr;w O t so t do npo
nlo,p lo,p
po,p
set of railway terminals set of trains set of trains passing in terminal n; RN n # R fixed unitary cost for handling boxes at terminal n hourly unitary cost for storing boxes at terminal n maximum number of handling operations (loading and unloading) allowed for each train at terminal n set of railway links, L ¼ fði; jÞ : i; j 2 N ; i–jg; we denote also a railway link as l = (i, j), where i is the tail of l and j is its head number of railway links covered by train r schedule for train r, specifying the nlr links covered by r and the arrival and departure times at the visited terminals cost related to the use of train r maximum bearable weight for train r set of wagons of train r length of wagon w of train r maximum bearable weight for wagon w of train r set of orders release time of order o at the origin terminal deadline for order o at the destination terminal number of alternative available routes connecting the origin and destination terminals of order o number of railway links of route p of order o vector 1 nlo,p indicating the sequence of railway links of route p of order o cost of route p of order o (this cost is weighted according to the route priority)
590
B ob fb gb
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
set of boxes order including box b; ob 2 O length of box b weight of box b
3. The planning approach The proposed planning procedure is divided in two sequential phases, corresponding to a pre-analysis and to the solution of a mathematical programming problem, as described in the following. 3.1. The pre-analysis The pre-analysis consists of an algorithm that computes all the sequences of trains available for serving each order, taking into account the network structure, the timetables and the stop sequences of trains, the origin, the destination and the time requirements of the order. Note that we refer to sequences of trains because each box is generally transported by more than one train from its origin to its destination (the maximum number of trains that a box can change is given by the number of railway links of the route). The pre-analysis algorithm is applied to one order at a time without taking into account capacity constraints but only considering feasibility in terms of time delivery and train connections. More in detail, this algorithm allows to compute all the feasible train sequences (for each order o and for each route p a priori specified) as those available to connect the order origin and destination. The algorithm proceeds backward analysing the railway links in lo,p, starting from the last link of the route, then verifying which trains, among the ones scheduled to cover the links, can be used to serve the order according to their timetable. In particular this is done for the last railway link in lo,p by finding all the trains arriving at the link head not after the deadline t do and leaving from the link tail not before t so . Then, the algorithm proceeds backward along the route, selecting among the trains travelling on the considered link the ones allowing a feasible connection with the feasible trains already determined for the subsequent railway links. For example, in a path including three links we may find that trains 1, 2 and 3 cover the last link in a feasible way for an order. Proceeding backward, we can find that train 1 is also feasible for the second link and it allows a feasible connection with train 2 but not with train 3. Considering the first link we can found that train 1 is still feasible and also train 4 allows a feasible connection with train 1. After the completion of the backward analysis, a forward procedure is applied to determine all the feasible train sequences for an order o on a route p. In the considered example, we find that the order can be served on that route by four feasible train sequences, i.e., 1-1-1 (train 1 on the three links), 1-1-2 (train 1 on the first two links and train 2 on the third link), 4-1-1 (train 4 on the first link and train 1 on the following two links) and 4-1-2 (train 4 on the first link, train 1 on the second link and train 2 on the last link). The following quantities are obtained from the pre-analysis and are used as input data in the mathematical programming problem: So Ro;s S Uo;n;r S Lo;n;r S To;n;r
set of train sequences found for order o set of trains in sequence s of order o set of indices of sequences for order o which involve the unloading of boxes at terminal n from train r set of indices of sequences for order o which involve the loading of boxes at terminal n to train r set of indices of sequences for order o which involve a transfer operation at terminal n with train r (i.e., neither loading nor unloading operations are executed)
po,s co,s
route associated with train sequence s of order o cost of sequence s of order o
In addition, the following two quantities are computed:
uno;s 2 f0; 1g denoting whether the boxes of o following sequence s involve a train change at terminal n2N storage time (in hours) required for boxes of o following sequence s at terminal n 2 N
Hno;s
The cost co,s (cost of serving an order o with sequence s) is determined as the priority po;po;s of order o on route po,s multiplied by the sum of unloading and storage costs over all the involved terminals, i.e.:
co;s ¼ po;po;s
X
n2N
uno;s dn þ Hno;s qn ;
ð1Þ
3.2. The mathematical programming formulation The mathematical programming problem is defined to optimally assign each box to a train sequence and to a wagon of the trains in the sequence, by considering the sequences of trains for each order found in the pre-analysis and taking into account the physical characteristics of boxes and wagons, i.e., considering capacity and other operational constraints. The decision variables are the following. yb;s 2 f0; 1g; b 2 B; s 2 S ob , equal to 1 if box b is assigned to sequence s, 0 otherwise xb;s;r;w 2 f0; 1g; b 2 B; s 2 S ob ; r 2 Rob ;s ; w 2 W r , equal to 1 if box b is assigned to wagon w of train r in sequence s, 0 otherwise v b 2 f0; 1g; b 2 B, equal to 1 if box b is not served, 0 otherwise zr 2 f0; 1g; r 2 R, equal to 1 if train r is used, 0 otherwise The planning problem is stated with the following 0–1 linear programming formulation. Problem 1.
min
X X
cob ;s yb;s þ
X
kr zr þ M
r2R
b2B s2S ob
X
vb;
ð2Þ
b2B
subject to
X
yb;s þ v b ¼ 1;
b 2 B:
xb;s;r;w ¼ yb;s ;
b2B
ð3Þ
s2S ob
X
s 2 S ob
r 2 Rob ;s :
ð4Þ
w2W r
X
X b2B
s2 S Lo
b
X
T ;n;r [S o
g b xb;s;r;w 6 gr zr ;
n2N
r 2 RN n :
w2W r
b ;n;r
ð5Þ X
X b2B
s2
S Lo ;n;r [S To ;n;r b b
g b xb;s;r;w 6 gr;w ;
n2N
r 2 RN n
w 2 Wr :
ð6Þ X
X b2B
s2 S Lo
b ;n;r
[S To
fb xb;s;r;w 6 f r;w ;
n2N
r 2 RN n
w 2 Wr :
b ;n;r
ð7Þ
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
X b2B
X s2
X
S Lo ;n;r [S U ob ;n;r b
yb;s 2 f0; 1g xb;s;r;w 2 f0; 1g
xb;s;r;w 6 hn ;
n2N
r 2 RN n :
ð8Þ
w2W r
b2B
s 2 S ob ;
b2B
s 2 S ob
ð9Þ r 2 Rob ;s
w 2 Wr ;
591
The partially fixed MIP sub-problems are solved by imposing both a maximum computation time and a maximum number of nodes for the MIP solver. Other heuristic details are here omitted for the sake of brevity, as, for example, the use of a tabu list in order to avoid variables yb,s entering and leaving set F too frequently.
ð10Þ
4.2. The randomized neighborhood search heuristic The second method used to face Problem 1 is called Randomized Neighborhood Search (RANS) and it is a general purpose MIP heuristic that simply operates similarly to an iterated local search using randomization to generate the explored neighborhood. The RANS algorithm starts from a feasible first incumbent solution for the original MIP problem and proceeds iteratively following three main steps:
zr 2 f0; 1g
r 2 R;
ð11Þ
v b 2 f0; 1g
b 2 B:
ð12Þ
The objective function (2) is the minimization of the costs associated with train sequences (computed in the pre-analysis), train costs, and the penalty for not serving orders (here M denotes a very large constant). Constraints (3) impose that, if served, each box is assigned to one and only one train sequence, while (4) establish that if a box is assigned to a train sequence, then it must be assigned to one wagon of each train belonging to that sequence. Constraints (5) impose for each train that the maximum bearable weight is not exceeded. Similarly, constraints (6) and (7) ensure that boxes assigned to wagons are compatible with the wagon length and weight limitations. Constraints (8) impose that the maximum number of handling operations to be performed for each train at a given terminal is not exceeded. Finally, constraints (9)–(12) are the definitions of the variables.
4. Two MIP heuristic approaches In order to find good quality solutions to Problem 1 in an acceptable computation time, we devised and experimented two MIP heuristic approaches: an ad hoc procedure integrated with the MIP solver and a general purpose heuristic method that iteratively calls the MIP solver as a black box tool. 4.1. The ad hoc heuristic Problem 1 is a generalized assignment which involves two main nested assignment decisions: the assignment of boxes to sequences, through yb,s variables, and the assignment of boxes to train wagons according to the selected sequence, through xb,s,r,w variables. The method focuses on the yb,s binary variables, as they determine the highest level assignment decisions, so mainly influencing the dimension of the solution space. Let Y denote the set of all yb,s variables. The heuristic iteratively solves a sequence of MIP sub-problems obtained from the original Problem 1 having fixed a subset F # Y of yb,s binary variables by imposing both their lower and upper bound equal to 0. The method then analyzes the reduced costs for the fixed variables from the solution of the branch-andcut node producing the optimal integer solution at an iteration, and determines for which variables it is convenient to reset the original upper bound to 1. This method is called Ad hoc MIP Heuristic (AMH). More in details, this method starts including all the yb,s variables in F. Then, at any iteration k, the subset YnF of variables to be explored is determined, for each box b, as follows: by removing from F (i.e., resetting the upper bound to 1) the k1 k1 k1 variable yb;s0 such that rb;s is the 0 ¼ mini fr b;i g, where r b;s reduced cost of yb,s at iteration k 1; by inserting in F (i.e., fixing the upper bound to 0) the variable yb;s00 in YnF such that yb;s00 ¼ 0 at iteration k 1 and n o k1 . r k1 b;s00 ¼ mini r b;i
1. Solution destruction: a subset V # G of k binary and integer variables is randomly selected and the partially fixed MIP subproblem (obtained from the original one by fixing the variables in GnV to their values in the incumbent solution) is defined. Note that k is a parameter initialized as k = 0.1 jGj and automatically tuned at each iteration (step 3). 2. Solution construction: this is a local search step consisting in solving the partially fixed MIP sub-problem by calling a MIP solver, having imposed a maximum computation time (note that also this time is a parameter automatically determined by the algorithm depending on the time needed to solve the linear relaxation of the MIP problem). Then, if a new best solution is found, the incumbent is updated. 3. Parameter adjustment and iteration: if the MIP sub-problem is optimally solved in the available time, then k = k 1.1, otherwise k = k 0.9, and a new iteration is started. In this simple way the algorithm controls the dimension of the explored neighborhood (i.e., the number of binary/integer variables), depending on the experienced difficulty in solving sub-problems. Similarly to the ad hoc heuristic, RANS rationale is to face the solution of a MIP problem by iteratively solving problems of smaller dimensions. RANS generates the MIP sub-problems by hard fixing a subset of binary and integer variables whose cardinality is automatically adjusted in order to find an optimal solution in the maximum time allotted. Note that this corresponds to perform a local search in a neighborhood of the incumbent solution of variable dimension. 5. Experimental results We ran a set of experimental tests on a set of randomly generated problem instances in order to evaluate the performance of the planning procedure and the solution approaches presented above. Both the pre-analysis algorithm and the solution schemes for Problem 1 were implemented in C]; in particular, the 0–1 linear optimization problem was solved using Cplex 11.0 and the ILOG Concert technology for building the model from the C] language. We introduced two instance ‘‘difficulty level’’ indexes, called AverageLoad and PeakLoad, that we computed before solving the 0–1 linear problem. The AverageLoad is the ratio between the required train slots and the available ones, defined as:
Av erageLoad ¼
P
o2O nbo
nlo
g nw nlr
ð13Þ
The denominator of (13) is an estimation of the number of available slots for serving boxes. This is obtained multiplying g, that is the a priori computed average number of boxes that can be carried on a full wagon, nw, the average number of wagons for trains, and nlr,
592
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
the number of different link-train pairs included in all the train sequences found in the pre-analysis for the orders. The numerator of (13) indicates the average number of boxes to be served on the different railway links, since nbo is the number of boxes of order o and nlo is the average number of railway links composing the routes of order o. The PeakLoad characterizes the most critical link-train pair and is defined as:
PeakLoad ¼
maxðl;rÞ2LR
P
o2O nbo ðl; rÞ
g nw
ð14Þ
:
The denominator in (14) is the estimation of the available slots on a train, whereas the numerator estimates the load of the most critical link-train pair, as the maximum number of boxes that must be transported on any railway link by any train. In (14) nbo ðl; rÞ is the average number of boxes of order o associated with link-train pair (l, r), obtained as the number of train sequences of order o including the pair (l, r) times the quantity nbo =jS o j, i.e., the average number of boxes of order o for sequences serving o. The difficulty in solving an instance should increase for higher values of these two indexes; in particular, PeakLoad values larger than 1 should characterize instances for which some boxes could be unserved. We built a set of problem instances taking into account some realistic data, when available, relevant to the type of boxes and wagons, the railway network and train timetables, and randomly generating the transportation demand. Moreover, the cost terms adopted in this experimental campaign have been provided by logistic operators. As already pointed out in the Introduction, Metrocargo system is at a prototype level and the considered railway network with fast transfer systems at terminals is not applied yet in Italy. This is the reason why realistic (and not real) data are used to build the problem instances. We analysed three different scenarios, described in Table 1, which will be referred to as Small, Medium and Large scenario, characterized by different dimensions of the railway network and the available trains. The considered planning horizon is one week for the Small and Medium scenario, whereas two weeks are considered for the Large one. The train routes and timetables are supposed to be the same on each day, the available routes for each origin-destination pair range between 1 and 3 and each train covers either 2 or 3 railway links. We generated six groups of instances, one for the Small scenario, two for the Medium and three for the Large one. The type of boxes was selected among 13 considered types from the discrete uniform distribution U[1, 13]; the origin and destination terminals were generated from U[1, nt] where nt is the number of terminals in the network for a scenario; the order release time was chosen
Table 1 Dimensions of the three scenarios. Scenario
Number of nodes
Number of arcs
Number of trains
Small Medium Large
10 13 13
34 50 50
98 154 308
as t so U½1; 48 (expressed in hours) and the deadline as tdo ¼ tso þ d where d U[96, 120]. In this way, we assume that boxes must reach their destination from their origin in 4–5 days. Finally, Table 2 reports specific data and average difficulty level indexes (including the average number of 0-1 variables) for the groups of instances. The computational tests were executed on a 2.8 GHz Pentium 4 computer with 2 GB of RAM. For each group we generated 5 instances: we firstly applied the pre-analysis procedure, then the three solution methods, corresponding to solving Problem 1 with Cplex 11.0 and applying the two proposed MIP heuristics (AMH and RANS), again using Cplex 11.0 as MIP solver. For the first method, we solved Problem 1 giving a 2 hours time limit to the Cplex 11.0 solver. Moreover, we imposed for AMH a time limit of 600 seconds for solving the MIP sub-problems and we fixed the node limit for the branch-and-cut exploration to 1000 nodes for the instances with a number of variables not greater than 30000, whereas we fixed this limit to a single node for larger instances. The maximum time limit for both the MIP heuristics was set to 1 hour. In order to evaluate the effectiveness of the AMH and RANS approaches, we report the obtained results in Table 3 together with the ones produced by Cplex. Moreover, the characteristics of each instance are specified, in terms of AverageLoad and PeakLoad indexes, and number of variables in the 0–1 LP problem. The computational results for the three methods are reported, as regards the value of the objective function and the number of not assigned P boxes (NotAss) that corresponds to b2B v b . Note that Cplex was not able to improve the first trivial solution (no boxes assigned) for one instance in groups A and D, for three instances in group B and for all the instances in groups C, E and F. In other words, only few instances are solved in a satisfactory way by Cplex. In addition, note that the solver always reached the limit of 7200 seconds. Moreover, from the results in Table 3 we can clearly appreciate the effectiveness of both the MIP heuristics that, in half the time given to the MIP solver, were able to find solutions that are largely better. Only for the instance A2 Cplex found a slightly better solution, whereas the average percentage deviations, obtained as (Heuristic Objective Cplex Objective)/Cplex Objective, are about 67% for AMH and 77% for RANS. In Table 4 the two MIP heuristics are compared, showing the percentage deviation of RANS results from AMH, computed as (RANS Objective AMH Objective)/AMH Objective, and the RANS deviation in the number of not assigned boxes, computed as RANS NotAss AMH NotAss. Comparing the two proposed MIP heuristics, we can observe the prevalence of RANS over AMH: this latter, in fact, produced a better objective solution only for 3 instances over 35 and only in one case it was able to assign one box more than RANS. The reason of the better performance of RANS with respect to AMH could be explained by more carefully analysing the behaviour of the two methods. AMH iterates the solution of MIP sub-problems where a subset of yb,s variables are fixed (bounded to zero); this in principle was supposed to generate simpler MIP sub-problems from whose solutions the method should gain guidance on the most convenient assignment decisions. However, the simplification due to
Table 2 Characteristics of the groups of instances. Group
Scenario
Orders
Boxes per order
AverageLoad
PeakLoad
Variables
A B C D E F
Small Medium Medium Large Large Large
60 80 100 100 120 140
U[6, 8] U[6, 8] U[8, 10] U[4, 6] U[5, 7] U[6, 8]
0.38 0.38 0.59 0.18 0.26 0.32
2.45 1.95 3.47 1.10 1.72 2.05
142139 218164 322087 437293 669931 859769
593
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594 Table 3 The computational results with Cplex, AMH and RANS. Inst. ID
AverageLoad
PeakLoad
No. of variables
Cplex objective
Cplex NotAss
AMH objective
AMH NotAss
RANS objective
RANS NotAss
A1 A2 A3 A4 A5
0.36 0.38 0.41 0.38 0.38
2.22 2.01 2.54 2.99 2.49
133251 136996 163414 136600 140435
3391895 9379194 424000000 19403106 12395334
3 9 424 19 12
2386326 9382326 3386911 389755 5388565
2 9 3 0 5
1389165 9379687 1388865 393552 4391352
1 9 1 0 4
B1 B2 B3 B4 B5
0.39 0.38 0.38 0.39 0.38
2.41 2.04 2 1.9 1.38
185589 197942 242570 234597 230124
73516161 50578470 563000000 552000000 578000000
73 50 563 552 578
42523507 24575480 14572205 5560182 4561566
42 24 14 5 4
41514833 17563662 10567689 560189 553609
41 17 10 5 5
C1 C2 C3 C4 C5
0.56 0.58 0.59 0.6 0.6
3.19 2.98 3.81 3.69 3.7
307070 301130 331078 335544 335615
910000000 901000000 895000000 890000000 895000000
910 901 895 890 895
113737485 129748775 153680582 124696449 132723643
113 129 153 124 132
101746837 102770643 124704947 106715182 101750611
101 102 124 106 101
D1 D2 D3 D4 D5
0.17 0.18 0.18 0.18 0.17
1.41 0.98 1.07 0.88 1.15
430006 428932 436317 468674 422538
8705542 33617886 1712172 504000000 15638562
8 33 1 504 15
575054 591723 572312 590978 578566
0 0 0 0 0
573427 585129 564187 580400 566335
0 0 0 0 0
E1 E2 E3 E4 E5
0.26 0.26 0.26 0.25 0.25
1.76 1.62 1.52 1.95 1.76
722530 732747 687204 658116 549056
734000000 707000000 719000000 711000000 720000000
734 707 719 711 720
2879459 5847630 2784957 2819057 796262
2 5 2 2 0
858728 835393 774731 810395 803483
0 0 0 0 0
F1 F2 F3 F4 F5
0.31 0.33 0.33 0.3 0.33
2.54 2.08 2.12 1.44 2.06
864289 890972 927211 730560 885812
981000000 983000000 994000000 998000000 979000000
981 983 994 998 979
6055668 10094255 31081102 14038610 18107729
6 10 31 14 18
1074804 1087820 10148714 1062927 5122048
1 1 10 1 5
Table 4 The comparison between AMH and RANS results. Inst. ID
RANS objective % dev. from AMH
A1 A2 A3 A4 A5
41,79 0,03 58,99 0,97 18,51
RANS NotAss dev. from AMH 1 0 2 0 1
B1 B2 B3 B4 B5
2,37 28,53 27,48 89,92 87,86
1 7 4 0 1
C1 C2 C3 C4 C5
10,54 20,79 18,85 14,42 23,34
12 27 29 18 31
D1 D2 D3 D4 D5
0,28 1,11 1,42 1,79 2,11
0 0 0 0 0
E1 E2 E3 E4 E5
70,18 85,71 72,18 71,25 0,91
2 5 2 2 0
F1 F2 F3 F4 F5
82,25 89,22 67,35 92,43 71,71
5 9 21 13 13
variable fixing was not enough to allow the solution of the MIP subproblems in short times; if, for example, we consider instance A1,
which has 133251 binary variables, we can observe that 94.3% are xb,s,r,w variables that are never fixed and only 5.7% are yb,s. Therefore the number of iterations performed by the AMH in one hour for instance A1 was only 49. With analogous considerations, the AMH ability to explore the solution space seems to be not very extended, since changes in variable fixing from one iteration to another one are rather limited. On the other hand, RANS does not privilege any subsets of variables for hard fixing (it is a general purpose MIP heuristic) and it uses an adaptation mechanism to control the dimension of the MIP sub-problems (i.e., the number of not fixed variables) in order to optimally solve them most of the times. In addition, this provides RANS with a good solution space exploration/diversification capability and allows RANS to execute a larger number of iterations (203 in case of instance A1) that greatly compensate for the use of a random fixing rule. As already highlighted, the tested instances can be considered representative of real cases, in particular the largest ones, corresponding to a network with 13 terminals and 50 railway links, 308 trains and an average request of about 1000 boxes to transport. Also the time horizon used for the tested instances can be considered representative of realistic applications. Of course, the proposed approach is suitable for an off-line planning, for instance one day in advance, for which the considered computation times are acceptable. 6. Conclusions In this paper we have proposed a planning procedure for serving some transportation requests (expressed in terms of origin, destination, release time, deadline, number of boxes, and so on) by using a railway network with rapid transhipment terminals. In such a system the boxes can change different trains from origin to destination and the boxes of a same order can be served in
594
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
different ways (different routes and sequences of trains). The proposed procedure is composed of a pre-analysis phase and, then, the statement and solution of a 0–1 linear programming problem. Two MIP heuristic approaches have been proposed for solving large problem instances. The performed experimental tests have shown the effectiveness of the proposed planning approach and its possible application to real cases. Finally, it can be noted that the proposed planning approach can be easily integrated in planning procedures of an overall railway network in which rapid transhipment terminals coexist with traditional ones. As a matter of fact, the planning procedure proposed in this work considers only terminals with fast transfer equipment and the relevant trains. Other trains and terminals of the same network can keep being treated with traditional planning methods. Acknowledgments This work was partially supported by the Regional Authority of Liguria in Italy within the project ‘‘ModLog – Modelling and Optimization of Dynamic Logistic Systems’’, developed at the Italian Centre of Excellence for Integrated Logistics and funded by ‘‘Parco Scientifico e Tecnologico della Liguria’’. References [1] Y.M. Bontekoning, C. Macharis, J.J. Trip, Is a new applied transportation research field emerging? – A review of intermodal rail-truck freight transport literature, Transportation research part A 38 (2004) 1–34. [2] T.G. Crainic, K.H. Kim, Intermodal transportation, in: C. Barnhart, G. Laporte (Eds.), Transportation, North Holland, 2007. [3] T.G. Crainic, Long-haul freight transportation, in: R.W. Hall (Ed.), Handbook of Transportation Science, Kluwer Academic publishers., 2003, pp. 451–516. [4]. [5] N. Wieberneit, Service network design for freight transportation: A review, OR Spectrum 30 (2008) 77–112.
[6] N. Bostel, P. Dejax, Models and algorithms for container allocation problems on trains in a rapid transshipment shunting yard, Transportation Science 32 (1998) 370–379. [7] P. Corry, E. Kozan, An assignment model for dynamic load planning of intermodal trains, Computers and Operations Research 33 (2006) 1–17. [8] J.-F. Cordeau, P. Toth, D. Vigo, A survey of optimization models for train routing and scheduling, Transportation Science 32 (1998) 380–404. [9] C.A. Yano, A.M. Newman, Scheduling trains and containers with due dates and dynamic arrivals, Transportation Science 35 (2001) 181–191. [10] H.N. Newton, C. Barnhart, P.H. Vance, Constructing railroad blocking plans to minimize handling costs, Transportation Science 32 (1998) 330–345. [11] S. Sacone, S. Siri. ‘‘A planning approach for freight transportation operations in railway networks. In Proceedings of the IEEE International Conference on Automation and Logistics, Qingdao, China (2008). [12] S. Sacone, S. Siri. ‘‘Optimal planning of order-based operations for intermodal freight transportation. In: Proceedings of the International Trade and Freight Transportation Conference, Cyprus (2008). [13] M.L. Fisher, R. Jaikamur, L.N. VanWassenhove, A multiplier adjustment method for the generalized assignment problem, Management Science 32 (1986) 1095–1103. [14] M. Fischetti, A. Lodi, Local branching, Mathematical Programming 98 (2003) 23–47. [15] E. Danna, E. Rothberg, C. Le Pape, Exploring relaxation induced neighborhoods to improve MIP solutions, Mathematical Programming 102 (2005) 71–90. [16] P. Hansen, N. Mladenovic´, D. Urosevic´, Variable neighborhood search and local branching, Computers and Operations Research 33 (2010) 3034–3045. [17] V. Maniezzo, T. Stützle, S. Voß, (Eds.). Matheuristics – Hybridizing Metaheuristics and Mathematical Programming, Annals of Information Systems 10 (2010). [18] A.R. Clark, Hybrid heuristics for planning lot setups and sizes, Computers and Industrial Engineering 45 (2003) 545–562. [19] A.M. Newman, M. Kuchta, Using aggregation to optimize long-term production planning at an underground mine, European Journal of Operational Research 176 (2007) 1205–1218. [20] M.A. Lejeune, A variable neighborhood decomposition search method for supply chain management planning problems, European Journal of Operational Research 175 (2006) 959–976. [21] K.F. Doerner, V. Schmid, Survey: Matheuristics for rich vehicle routing problems Hybrid Metaheuristics 2010, Lecture Notes in Computer Science 6373 (2010) 206–221. [22] J. Lazic´, S. Hanafi, N. Mladenovic´, D. Urosevic´, Variable neighbourhood decomposition search for 0–1 mixed integer programs, Computers and Operations Research 37 (2010) 1055–1067.
Contents lists available at ScienceDirect
European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor
Production, Manufacturing and Logistics
Freight transportation in railway networks with automated terminals: A mathematical model and MIP heuristic approaches D. Anghinolfi, M. Paolucci, S. Sacone, S. Siri ⇑ Department of Communications, Computer and System Sciences, University of Genova, Italy
a r t i c l e
i n f o
Article history: Received 10 December 2010 Accepted 11 May 2011 Available online 17 May 2011 Keywords: Freight transportation Optimal planning Mathematical programming MIP heuristics
a b s t r a c t In this paper we propose a planning procedure for serving freight transportation requests in a railway network with fast transfer equipment at terminals. We consider a transportation system where different customers make their requests (orders) for moving boxes, i.e., either containers or swap bodies, between different origins and destinations, with specific requirements on delivery times. The decisions to be taken concern the route (and the corresponding sequence of trains) that each box follows in the network and the assignment of boxes to train wagons, taking into account that boxes can change more than one train and that train timetables are fixed. The planning procedure includes a pre-analysis step to determine all the possible sequences of trains for serving each order, followed by the solution of a 0–1 linear programming problem to find the optimal assignment of each box to a train sequence and to a specific wagon for each train in the sequence. This latter is a generalized assignment problem which is NP-hard. Hence, in order to find good solutions in acceptable computation times, two MIP heuristic approaches are proposed and tested through an experimental analysis considering realistic problem instances. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction Freight transportation is a very important activity in our society and it is also a complex domain to be studied. Freight transportation planning has been often faced by mathematical programming and optimization approaches and dealt with in many research works (see for instance the reviews [1,2]). In this paper, a kind of long haul transportation of containerized goods [3] is considered, characterized by a railway network whose terminals are automated, i.e., equipped with transhipment devices allowing fast train–train transfer operations. The study of this transportation service derives from a research project carried out with an Italian company that is developing a new automated container handling system named Metrocargo [4]. The adoption of Metrocargo in the national railway network would provide a very innovative transportation service for the Italian logistic context. Our main interest in this project is the definition of optimization methods for the relevant off-line planning. We consider the point of view of a central decision maker who provides the transportation service in the railway network and has to satisfy transportation requests (orders) from customers. The timetable and stop sequence of trains in the railway network are fixed by the provider of the transportation service and are not ⇑ Corresponding author. Tel.: +39 019 21945137; fax: +39 019 21945212. E-mail address: [email protected] (S. Siri). 0377-2217/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2011.05.013
matter of decision. The orders coming from customers regard the movement of a set of boxes (containers or swap bodies) from an origin to a destination within a deadline. We aim at generating transportation plans for orders considering the boxes similarly to passengers planning a train trip: for each box of each order we determine the sequence of trains and train changes at terminals so that the box can reach its destination in time. The planning problem faced in this paper combines two main decisions. The former concerns the routing of boxes in the network and the selection of train changes at terminals; the latter decision deals with the assignment of boxes to wagons of the selected trains. These two decisions have been generally considered separately in the literature. The box routing and train selection are often related to aggregate flows of goods represented as a set of commodities, as in the general formulation provided in the review [5]. In this paper, instead, we determine the routes and the sequences of trains so that the boxes of a same order can reach their common destination independently from each other, i.e., even through different routes and with different sequences of trains. The decision on wagon assignment has already been treated in the literature usually considering a single terminal. In [6] the authors assume a rail–rail container terminal with a rapid transshipment yard and they define different models and algorithms for determining the optimal placement of containers on arriving trains and on departing trains, in order to minimize the movements of containers in the terminal. Analogously, in [7] the
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
so-called load planning problem (assignment of containers to train slots) is faced for an intermodal terminal where containers are transferred between trucks and trains. Differently, our purpose in this work is optimally matching the whole transportation demand for the railway network with the overall available transportation offer (the fixed train timetables and stop sequences). Specifically, we determine the minimum cost assignment of the boxes to wagons of the different trains in order to satisfy train and wagon capacity constraints. Since the train timetables and stop sequences are fixed input data, this paper does not deal with train routing and scheduling, as done for instance in [8,9], or with other classical railway problems, such as the definition of railroad blocking plans [10]. Even though box routing and train selection, as well as wagon assignment, have been already treated in the literature, the main novelty of this paper stands in considering them jointly. As detailed in the following, the considered planning problem (a preliminary version can be found in [11,12]) basically consists in a generalized assignment problem (which is known to be NP-hard [13]) that we formulate as a 0–1 linear programming problem. This kind of problems can be faced by means of general purpose mixed integer programming (MIP) solvers that seek for the optimal solution by exploiting branch-and-bound or branch-and-cut algorithms. However, most MIP problems are very difficult to solve when their dimensions become too large. For this reason, complex combinatorial optimization problems from both academic research and real world applications have been tackled by specialized heuristics or metaheuristics. Recently, new approaches have been proposed to combine ideas from metaheuristics with MIP solver algorithms [14–16], and they are generally referred to as matheuristics [17]. In many cases, matheuristics have been adopted for large scale problems coming from real applications in order to reduce the solution search space and to speed up the computation (e.g., in production planning [18], in production scheduling [19], in supply chain management [20] and in vehicle routing and transportation logistics [21]). In this paper, we propose two matheuristic approaches for the considered planning problem; the first one is an ad hoc procedure integrated with the MIP solver to drive and restrict the solution space exploration, whereas the second one is a general purpose heuristic method to solve MIP problems by iteratively calling the MIP solver as a black box tool. This latter approach shares some concepts with the MIP heuristics recently appeared in literature as Local Branching [14], Relaxation Induced Neighborhood Search [15] and Variable Neighborhood Decomposition Search [22], but, differently from them, it exploits only a randomization mechanism to guide the MIP solver. This paper is organized as follows. In Section 2 the main aspects of the considered problem are introduced; in Section 3 the preanalysis and the mathematical programming formulation adopted in the planning procedure are described in detail. Then, in Section 4 the two proposed MIP heuristic approaches are introduced. The experimental results are discussed in Section 5 and some conclusive remarks are finally reported in Section 6. 2. Problem description The Metrocargo system is a new technology, now at a prototype level, allowing to load and unload cargo units from trains in a horizontal way and under the electric feeding line, thanks to some shuttles and a properly devised storage area where all the operations are realized automatically. With this technology rail–rail transfer operations at terminals become fast. The planning procedure proposed in this paper refers to a railway network with automated terminals equipped with Metrocargo system or any other innovative system allowing fast transfers. With such systems, the boxes can be considered as passengers that change trains in their
589
route to destination by using the available (i.e., scheduled) train services. The transportation demand to be satisfied is given by a set of orders. For each order we consider the origin and destination, the release time (i.e., the earliest time from which the boxes are available at the origin), the deadline (i.e., the latest time by which the boxes must be delivered at destination), and the number of boxes, each one characterized by length and weight. The railway network is composed of railway terminals and railway links connecting the terminals. We assume that for each pair of terminals in the network, i.e., for each origin–destination pair, the available routes (defined as sequences of railway links) are specified a priori. Such routes are identified by the transportation service provider among all the possible link sequences connecting an origin–destination pair, as the most effective ones to provide the transportation service. Note that this is usually done in order to eliminate possible but not reasonable/acceptable routes. Moreover, the available routes connecting each origin–destination pair are a priori ordered by the transportation service provider according to a given priority. The goal of the proposed planning procedure is to make the planning of transportation operations for all the orders within a specified time horizon. The problem is to determine for each box of each order (a) the route that it must follow from its origin to its destination, (b) the sequence of trains that it must use along such route so that the order time constraints are satisfied, and (c) the wagons used to transport it for each train in the sequence. The problem data are listed in the following.
N R RN n dn qn hn
L
nlr Sr
kr gr Wr f r;w gr;w O t so t do npo
nlo,p lo,p
po,p
set of railway terminals set of trains set of trains passing in terminal n; RN n # R fixed unitary cost for handling boxes at terminal n hourly unitary cost for storing boxes at terminal n maximum number of handling operations (loading and unloading) allowed for each train at terminal n set of railway links, L ¼ fði; jÞ : i; j 2 N ; i–jg; we denote also a railway link as l = (i, j), where i is the tail of l and j is its head number of railway links covered by train r schedule for train r, specifying the nlr links covered by r and the arrival and departure times at the visited terminals cost related to the use of train r maximum bearable weight for train r set of wagons of train r length of wagon w of train r maximum bearable weight for wagon w of train r set of orders release time of order o at the origin terminal deadline for order o at the destination terminal number of alternative available routes connecting the origin and destination terminals of order o number of railway links of route p of order o vector 1 nlo,p indicating the sequence of railway links of route p of order o cost of route p of order o (this cost is weighted according to the route priority)
590
B ob fb gb
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
set of boxes order including box b; ob 2 O length of box b weight of box b
3. The planning approach The proposed planning procedure is divided in two sequential phases, corresponding to a pre-analysis and to the solution of a mathematical programming problem, as described in the following. 3.1. The pre-analysis The pre-analysis consists of an algorithm that computes all the sequences of trains available for serving each order, taking into account the network structure, the timetables and the stop sequences of trains, the origin, the destination and the time requirements of the order. Note that we refer to sequences of trains because each box is generally transported by more than one train from its origin to its destination (the maximum number of trains that a box can change is given by the number of railway links of the route). The pre-analysis algorithm is applied to one order at a time without taking into account capacity constraints but only considering feasibility in terms of time delivery and train connections. More in detail, this algorithm allows to compute all the feasible train sequences (for each order o and for each route p a priori specified) as those available to connect the order origin and destination. The algorithm proceeds backward analysing the railway links in lo,p, starting from the last link of the route, then verifying which trains, among the ones scheduled to cover the links, can be used to serve the order according to their timetable. In particular this is done for the last railway link in lo,p by finding all the trains arriving at the link head not after the deadline t do and leaving from the link tail not before t so . Then, the algorithm proceeds backward along the route, selecting among the trains travelling on the considered link the ones allowing a feasible connection with the feasible trains already determined for the subsequent railway links. For example, in a path including three links we may find that trains 1, 2 and 3 cover the last link in a feasible way for an order. Proceeding backward, we can find that train 1 is also feasible for the second link and it allows a feasible connection with train 2 but not with train 3. Considering the first link we can found that train 1 is still feasible and also train 4 allows a feasible connection with train 1. After the completion of the backward analysis, a forward procedure is applied to determine all the feasible train sequences for an order o on a route p. In the considered example, we find that the order can be served on that route by four feasible train sequences, i.e., 1-1-1 (train 1 on the three links), 1-1-2 (train 1 on the first two links and train 2 on the third link), 4-1-1 (train 4 on the first link and train 1 on the following two links) and 4-1-2 (train 4 on the first link, train 1 on the second link and train 2 on the last link). The following quantities are obtained from the pre-analysis and are used as input data in the mathematical programming problem: So Ro;s S Uo;n;r S Lo;n;r S To;n;r
set of train sequences found for order o set of trains in sequence s of order o set of indices of sequences for order o which involve the unloading of boxes at terminal n from train r set of indices of sequences for order o which involve the loading of boxes at terminal n to train r set of indices of sequences for order o which involve a transfer operation at terminal n with train r (i.e., neither loading nor unloading operations are executed)
po,s co,s
route associated with train sequence s of order o cost of sequence s of order o
In addition, the following two quantities are computed:
uno;s 2 f0; 1g denoting whether the boxes of o following sequence s involve a train change at terminal n2N storage time (in hours) required for boxes of o following sequence s at terminal n 2 N
Hno;s
The cost co,s (cost of serving an order o with sequence s) is determined as the priority po;po;s of order o on route po,s multiplied by the sum of unloading and storage costs over all the involved terminals, i.e.:
co;s ¼ po;po;s
X
n2N
uno;s dn þ Hno;s qn ;
ð1Þ
3.2. The mathematical programming formulation The mathematical programming problem is defined to optimally assign each box to a train sequence and to a wagon of the trains in the sequence, by considering the sequences of trains for each order found in the pre-analysis and taking into account the physical characteristics of boxes and wagons, i.e., considering capacity and other operational constraints. The decision variables are the following. yb;s 2 f0; 1g; b 2 B; s 2 S ob , equal to 1 if box b is assigned to sequence s, 0 otherwise xb;s;r;w 2 f0; 1g; b 2 B; s 2 S ob ; r 2 Rob ;s ; w 2 W r , equal to 1 if box b is assigned to wagon w of train r in sequence s, 0 otherwise v b 2 f0; 1g; b 2 B, equal to 1 if box b is not served, 0 otherwise zr 2 f0; 1g; r 2 R, equal to 1 if train r is used, 0 otherwise The planning problem is stated with the following 0–1 linear programming formulation. Problem 1.
min
X X
cob ;s yb;s þ
X
kr zr þ M
r2R
b2B s2S ob
X
vb;
ð2Þ
b2B
subject to
X
yb;s þ v b ¼ 1;
b 2 B:
xb;s;r;w ¼ yb;s ;
b2B
ð3Þ
s2S ob
X
s 2 S ob
r 2 Rob ;s :
ð4Þ
w2W r
X
X b2B
s2 S Lo
b
X
T ;n;r [S o
g b xb;s;r;w 6 gr zr ;
n2N
r 2 RN n :
w2W r
b ;n;r
ð5Þ X
X b2B
s2
S Lo ;n;r [S To ;n;r b b
g b xb;s;r;w 6 gr;w ;
n2N
r 2 RN n
w 2 Wr :
ð6Þ X
X b2B
s2 S Lo
b ;n;r
[S To
fb xb;s;r;w 6 f r;w ;
n2N
r 2 RN n
w 2 Wr :
b ;n;r
ð7Þ
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
X b2B
X s2
X
S Lo ;n;r [S U ob ;n;r b
yb;s 2 f0; 1g xb;s;r;w 2 f0; 1g
xb;s;r;w 6 hn ;
n2N
r 2 RN n :
ð8Þ
w2W r
b2B
s 2 S ob ;
b2B
s 2 S ob
ð9Þ r 2 Rob ;s
w 2 Wr ;
591
The partially fixed MIP sub-problems are solved by imposing both a maximum computation time and a maximum number of nodes for the MIP solver. Other heuristic details are here omitted for the sake of brevity, as, for example, the use of a tabu list in order to avoid variables yb,s entering and leaving set F too frequently.
ð10Þ
4.2. The randomized neighborhood search heuristic The second method used to face Problem 1 is called Randomized Neighborhood Search (RANS) and it is a general purpose MIP heuristic that simply operates similarly to an iterated local search using randomization to generate the explored neighborhood. The RANS algorithm starts from a feasible first incumbent solution for the original MIP problem and proceeds iteratively following three main steps:
zr 2 f0; 1g
r 2 R;
ð11Þ
v b 2 f0; 1g
b 2 B:
ð12Þ
The objective function (2) is the minimization of the costs associated with train sequences (computed in the pre-analysis), train costs, and the penalty for not serving orders (here M denotes a very large constant). Constraints (3) impose that, if served, each box is assigned to one and only one train sequence, while (4) establish that if a box is assigned to a train sequence, then it must be assigned to one wagon of each train belonging to that sequence. Constraints (5) impose for each train that the maximum bearable weight is not exceeded. Similarly, constraints (6) and (7) ensure that boxes assigned to wagons are compatible with the wagon length and weight limitations. Constraints (8) impose that the maximum number of handling operations to be performed for each train at a given terminal is not exceeded. Finally, constraints (9)–(12) are the definitions of the variables.
4. Two MIP heuristic approaches In order to find good quality solutions to Problem 1 in an acceptable computation time, we devised and experimented two MIP heuristic approaches: an ad hoc procedure integrated with the MIP solver and a general purpose heuristic method that iteratively calls the MIP solver as a black box tool. 4.1. The ad hoc heuristic Problem 1 is a generalized assignment which involves two main nested assignment decisions: the assignment of boxes to sequences, through yb,s variables, and the assignment of boxes to train wagons according to the selected sequence, through xb,s,r,w variables. The method focuses on the yb,s binary variables, as they determine the highest level assignment decisions, so mainly influencing the dimension of the solution space. Let Y denote the set of all yb,s variables. The heuristic iteratively solves a sequence of MIP sub-problems obtained from the original Problem 1 having fixed a subset F # Y of yb,s binary variables by imposing both their lower and upper bound equal to 0. The method then analyzes the reduced costs for the fixed variables from the solution of the branch-andcut node producing the optimal integer solution at an iteration, and determines for which variables it is convenient to reset the original upper bound to 1. This method is called Ad hoc MIP Heuristic (AMH). More in details, this method starts including all the yb,s variables in F. Then, at any iteration k, the subset YnF of variables to be explored is determined, for each box b, as follows: by removing from F (i.e., resetting the upper bound to 1) the k1 k1 k1 variable yb;s0 such that rb;s is the 0 ¼ mini fr b;i g, where r b;s reduced cost of yb,s at iteration k 1; by inserting in F (i.e., fixing the upper bound to 0) the variable yb;s00 in YnF such that yb;s00 ¼ 0 at iteration k 1 and n o k1 . r k1 b;s00 ¼ mini r b;i
1. Solution destruction: a subset V # G of k binary and integer variables is randomly selected and the partially fixed MIP subproblem (obtained from the original one by fixing the variables in GnV to their values in the incumbent solution) is defined. Note that k is a parameter initialized as k = 0.1 jGj and automatically tuned at each iteration (step 3). 2. Solution construction: this is a local search step consisting in solving the partially fixed MIP sub-problem by calling a MIP solver, having imposed a maximum computation time (note that also this time is a parameter automatically determined by the algorithm depending on the time needed to solve the linear relaxation of the MIP problem). Then, if a new best solution is found, the incumbent is updated. 3. Parameter adjustment and iteration: if the MIP sub-problem is optimally solved in the available time, then k = k 1.1, otherwise k = k 0.9, and a new iteration is started. In this simple way the algorithm controls the dimension of the explored neighborhood (i.e., the number of binary/integer variables), depending on the experienced difficulty in solving sub-problems. Similarly to the ad hoc heuristic, RANS rationale is to face the solution of a MIP problem by iteratively solving problems of smaller dimensions. RANS generates the MIP sub-problems by hard fixing a subset of binary and integer variables whose cardinality is automatically adjusted in order to find an optimal solution in the maximum time allotted. Note that this corresponds to perform a local search in a neighborhood of the incumbent solution of variable dimension. 5. Experimental results We ran a set of experimental tests on a set of randomly generated problem instances in order to evaluate the performance of the planning procedure and the solution approaches presented above. Both the pre-analysis algorithm and the solution schemes for Problem 1 were implemented in C]; in particular, the 0–1 linear optimization problem was solved using Cplex 11.0 and the ILOG Concert technology for building the model from the C] language. We introduced two instance ‘‘difficulty level’’ indexes, called AverageLoad and PeakLoad, that we computed before solving the 0–1 linear problem. The AverageLoad is the ratio between the required train slots and the available ones, defined as:
Av erageLoad ¼
P
o2O nbo
nlo
g nw nlr
ð13Þ
The denominator of (13) is an estimation of the number of available slots for serving boxes. This is obtained multiplying g, that is the a priori computed average number of boxes that can be carried on a full wagon, nw, the average number of wagons for trains, and nlr,
592
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
the number of different link-train pairs included in all the train sequences found in the pre-analysis for the orders. The numerator of (13) indicates the average number of boxes to be served on the different railway links, since nbo is the number of boxes of order o and nlo is the average number of railway links composing the routes of order o. The PeakLoad characterizes the most critical link-train pair and is defined as:
PeakLoad ¼
maxðl;rÞ2LR
P
o2O nbo ðl; rÞ
g nw
ð14Þ
:
The denominator in (14) is the estimation of the available slots on a train, whereas the numerator estimates the load of the most critical link-train pair, as the maximum number of boxes that must be transported on any railway link by any train. In (14) nbo ðl; rÞ is the average number of boxes of order o associated with link-train pair (l, r), obtained as the number of train sequences of order o including the pair (l, r) times the quantity nbo =jS o j, i.e., the average number of boxes of order o for sequences serving o. The difficulty in solving an instance should increase for higher values of these two indexes; in particular, PeakLoad values larger than 1 should characterize instances for which some boxes could be unserved. We built a set of problem instances taking into account some realistic data, when available, relevant to the type of boxes and wagons, the railway network and train timetables, and randomly generating the transportation demand. Moreover, the cost terms adopted in this experimental campaign have been provided by logistic operators. As already pointed out in the Introduction, Metrocargo system is at a prototype level and the considered railway network with fast transfer systems at terminals is not applied yet in Italy. This is the reason why realistic (and not real) data are used to build the problem instances. We analysed three different scenarios, described in Table 1, which will be referred to as Small, Medium and Large scenario, characterized by different dimensions of the railway network and the available trains. The considered planning horizon is one week for the Small and Medium scenario, whereas two weeks are considered for the Large one. The train routes and timetables are supposed to be the same on each day, the available routes for each origin-destination pair range between 1 and 3 and each train covers either 2 or 3 railway links. We generated six groups of instances, one for the Small scenario, two for the Medium and three for the Large one. The type of boxes was selected among 13 considered types from the discrete uniform distribution U[1, 13]; the origin and destination terminals were generated from U[1, nt] where nt is the number of terminals in the network for a scenario; the order release time was chosen
Table 1 Dimensions of the three scenarios. Scenario
Number of nodes
Number of arcs
Number of trains
Small Medium Large
10 13 13
34 50 50
98 154 308
as t so U½1; 48 (expressed in hours) and the deadline as tdo ¼ tso þ d where d U[96, 120]. In this way, we assume that boxes must reach their destination from their origin in 4–5 days. Finally, Table 2 reports specific data and average difficulty level indexes (including the average number of 0-1 variables) for the groups of instances. The computational tests were executed on a 2.8 GHz Pentium 4 computer with 2 GB of RAM. For each group we generated 5 instances: we firstly applied the pre-analysis procedure, then the three solution methods, corresponding to solving Problem 1 with Cplex 11.0 and applying the two proposed MIP heuristics (AMH and RANS), again using Cplex 11.0 as MIP solver. For the first method, we solved Problem 1 giving a 2 hours time limit to the Cplex 11.0 solver. Moreover, we imposed for AMH a time limit of 600 seconds for solving the MIP sub-problems and we fixed the node limit for the branch-and-cut exploration to 1000 nodes for the instances with a number of variables not greater than 30000, whereas we fixed this limit to a single node for larger instances. The maximum time limit for both the MIP heuristics was set to 1 hour. In order to evaluate the effectiveness of the AMH and RANS approaches, we report the obtained results in Table 3 together with the ones produced by Cplex. Moreover, the characteristics of each instance are specified, in terms of AverageLoad and PeakLoad indexes, and number of variables in the 0–1 LP problem. The computational results for the three methods are reported, as regards the value of the objective function and the number of not assigned P boxes (NotAss) that corresponds to b2B v b . Note that Cplex was not able to improve the first trivial solution (no boxes assigned) for one instance in groups A and D, for three instances in group B and for all the instances in groups C, E and F. In other words, only few instances are solved in a satisfactory way by Cplex. In addition, note that the solver always reached the limit of 7200 seconds. Moreover, from the results in Table 3 we can clearly appreciate the effectiveness of both the MIP heuristics that, in half the time given to the MIP solver, were able to find solutions that are largely better. Only for the instance A2 Cplex found a slightly better solution, whereas the average percentage deviations, obtained as (Heuristic Objective Cplex Objective)/Cplex Objective, are about 67% for AMH and 77% for RANS. In Table 4 the two MIP heuristics are compared, showing the percentage deviation of RANS results from AMH, computed as (RANS Objective AMH Objective)/AMH Objective, and the RANS deviation in the number of not assigned boxes, computed as RANS NotAss AMH NotAss. Comparing the two proposed MIP heuristics, we can observe the prevalence of RANS over AMH: this latter, in fact, produced a better objective solution only for 3 instances over 35 and only in one case it was able to assign one box more than RANS. The reason of the better performance of RANS with respect to AMH could be explained by more carefully analysing the behaviour of the two methods. AMH iterates the solution of MIP sub-problems where a subset of yb,s variables are fixed (bounded to zero); this in principle was supposed to generate simpler MIP sub-problems from whose solutions the method should gain guidance on the most convenient assignment decisions. However, the simplification due to
Table 2 Characteristics of the groups of instances. Group
Scenario
Orders
Boxes per order
AverageLoad
PeakLoad
Variables
A B C D E F
Small Medium Medium Large Large Large
60 80 100 100 120 140
U[6, 8] U[6, 8] U[8, 10] U[4, 6] U[5, 7] U[6, 8]
0.38 0.38 0.59 0.18 0.26 0.32
2.45 1.95 3.47 1.10 1.72 2.05
142139 218164 322087 437293 669931 859769
593
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594 Table 3 The computational results with Cplex, AMH and RANS. Inst. ID
AverageLoad
PeakLoad
No. of variables
Cplex objective
Cplex NotAss
AMH objective
AMH NotAss
RANS objective
RANS NotAss
A1 A2 A3 A4 A5
0.36 0.38 0.41 0.38 0.38
2.22 2.01 2.54 2.99 2.49
133251 136996 163414 136600 140435
3391895 9379194 424000000 19403106 12395334
3 9 424 19 12
2386326 9382326 3386911 389755 5388565
2 9 3 0 5
1389165 9379687 1388865 393552 4391352
1 9 1 0 4
B1 B2 B3 B4 B5
0.39 0.38 0.38 0.39 0.38
2.41 2.04 2 1.9 1.38
185589 197942 242570 234597 230124
73516161 50578470 563000000 552000000 578000000
73 50 563 552 578
42523507 24575480 14572205 5560182 4561566
42 24 14 5 4
41514833 17563662 10567689 560189 553609
41 17 10 5 5
C1 C2 C3 C4 C5
0.56 0.58 0.59 0.6 0.6
3.19 2.98 3.81 3.69 3.7
307070 301130 331078 335544 335615
910000000 901000000 895000000 890000000 895000000
910 901 895 890 895
113737485 129748775 153680582 124696449 132723643
113 129 153 124 132
101746837 102770643 124704947 106715182 101750611
101 102 124 106 101
D1 D2 D3 D4 D5
0.17 0.18 0.18 0.18 0.17
1.41 0.98 1.07 0.88 1.15
430006 428932 436317 468674 422538
8705542 33617886 1712172 504000000 15638562
8 33 1 504 15
575054 591723 572312 590978 578566
0 0 0 0 0
573427 585129 564187 580400 566335
0 0 0 0 0
E1 E2 E3 E4 E5
0.26 0.26 0.26 0.25 0.25
1.76 1.62 1.52 1.95 1.76
722530 732747 687204 658116 549056
734000000 707000000 719000000 711000000 720000000
734 707 719 711 720
2879459 5847630 2784957 2819057 796262
2 5 2 2 0
858728 835393 774731 810395 803483
0 0 0 0 0
F1 F2 F3 F4 F5
0.31 0.33 0.33 0.3 0.33
2.54 2.08 2.12 1.44 2.06
864289 890972 927211 730560 885812
981000000 983000000 994000000 998000000 979000000
981 983 994 998 979
6055668 10094255 31081102 14038610 18107729
6 10 31 14 18
1074804 1087820 10148714 1062927 5122048
1 1 10 1 5
Table 4 The comparison between AMH and RANS results. Inst. ID
RANS objective % dev. from AMH
A1 A2 A3 A4 A5
41,79 0,03 58,99 0,97 18,51
RANS NotAss dev. from AMH 1 0 2 0 1
B1 B2 B3 B4 B5
2,37 28,53 27,48 89,92 87,86
1 7 4 0 1
C1 C2 C3 C4 C5
10,54 20,79 18,85 14,42 23,34
12 27 29 18 31
D1 D2 D3 D4 D5
0,28 1,11 1,42 1,79 2,11
0 0 0 0 0
E1 E2 E3 E4 E5
70,18 85,71 72,18 71,25 0,91
2 5 2 2 0
F1 F2 F3 F4 F5
82,25 89,22 67,35 92,43 71,71
5 9 21 13 13
variable fixing was not enough to allow the solution of the MIP subproblems in short times; if, for example, we consider instance A1,
which has 133251 binary variables, we can observe that 94.3% are xb,s,r,w variables that are never fixed and only 5.7% are yb,s. Therefore the number of iterations performed by the AMH in one hour for instance A1 was only 49. With analogous considerations, the AMH ability to explore the solution space seems to be not very extended, since changes in variable fixing from one iteration to another one are rather limited. On the other hand, RANS does not privilege any subsets of variables for hard fixing (it is a general purpose MIP heuristic) and it uses an adaptation mechanism to control the dimension of the MIP sub-problems (i.e., the number of not fixed variables) in order to optimally solve them most of the times. In addition, this provides RANS with a good solution space exploration/diversification capability and allows RANS to execute a larger number of iterations (203 in case of instance A1) that greatly compensate for the use of a random fixing rule. As already highlighted, the tested instances can be considered representative of real cases, in particular the largest ones, corresponding to a network with 13 terminals and 50 railway links, 308 trains and an average request of about 1000 boxes to transport. Also the time horizon used for the tested instances can be considered representative of realistic applications. Of course, the proposed approach is suitable for an off-line planning, for instance one day in advance, for which the considered computation times are acceptable. 6. Conclusions In this paper we have proposed a planning procedure for serving some transportation requests (expressed in terms of origin, destination, release time, deadline, number of boxes, and so on) by using a railway network with rapid transhipment terminals. In such a system the boxes can change different trains from origin to destination and the boxes of a same order can be served in
594
D. Anghinolfi et al. / European Journal of Operational Research 214 (2011) 588–594
different ways (different routes and sequences of trains). The proposed procedure is composed of a pre-analysis phase and, then, the statement and solution of a 0–1 linear programming problem. Two MIP heuristic approaches have been proposed for solving large problem instances. The performed experimental tests have shown the effectiveness of the proposed planning approach and its possible application to real cases. Finally, it can be noted that the proposed planning approach can be easily integrated in planning procedures of an overall railway network in which rapid transhipment terminals coexist with traditional ones. As a matter of fact, the planning procedure proposed in this work considers only terminals with fast transfer equipment and the relevant trains. Other trains and terminals of the same network can keep being treated with traditional planning methods. Acknowledgments This work was partially supported by the Regional Authority of Liguria in Italy within the project ‘‘ModLog – Modelling and Optimization of Dynamic Logistic Systems’’, developed at the Italian Centre of Excellence for Integrated Logistics and funded by ‘‘Parco Scientifico e Tecnologico della Liguria’’. References [1] Y.M. Bontekoning, C. Macharis, J.J. Trip, Is a new applied transportation research field emerging? – A review of intermodal rail-truck freight transport literature, Transportation research part A 38 (2004) 1–34. [2] T.G. Crainic, K.H. Kim, Intermodal transportation, in: C. Barnhart, G. Laporte (Eds.), Transportation, North Holland, 2007. [3] T.G. Crainic, Long-haul freight transportation, in: R.W. Hall (Ed.), Handbook of Transportation Science, Kluwer Academic publishers., 2003, pp. 451–516. [4]
[6] N. Bostel, P. Dejax, Models and algorithms for container allocation problems on trains in a rapid transshipment shunting yard, Transportation Science 32 (1998) 370–379. [7] P. Corry, E. Kozan, An assignment model for dynamic load planning of intermodal trains, Computers and Operations Research 33 (2006) 1–17. [8] J.-F. Cordeau, P. Toth, D. Vigo, A survey of optimization models for train routing and scheduling, Transportation Science 32 (1998) 380–404. [9] C.A. Yano, A.M. Newman, Scheduling trains and containers with due dates and dynamic arrivals, Transportation Science 35 (2001) 181–191. [10] H.N. Newton, C. Barnhart, P.H. Vance, Constructing railroad blocking plans to minimize handling costs, Transportation Science 32 (1998) 330–345. [11] S. Sacone, S. Siri. ‘‘A planning approach for freight transportation operations in railway networks. In Proceedings of the IEEE International Conference on Automation and Logistics, Qingdao, China (2008). [12] S. Sacone, S. Siri. ‘‘Optimal planning of order-based operations for intermodal freight transportation. In: Proceedings of the International Trade and Freight Transportation Conference, Cyprus (2008). [13] M.L. Fisher, R. Jaikamur, L.N. VanWassenhove, A multiplier adjustment method for the generalized assignment problem, Management Science 32 (1986) 1095–1103. [14] M. Fischetti, A. Lodi, Local branching, Mathematical Programming 98 (2003) 23–47. [15] E. Danna, E. Rothberg, C. Le Pape, Exploring relaxation induced neighborhoods to improve MIP solutions, Mathematical Programming 102 (2005) 71–90. [16] P. Hansen, N. Mladenovic´, D. Urosevic´, Variable neighborhood search and local branching, Computers and Operations Research 33 (2010) 3034–3045. [17] V. Maniezzo, T. Stützle, S. Voß, (Eds.). Matheuristics – Hybridizing Metaheuristics and Mathematical Programming, Annals of Information Systems 10 (2010). [18] A.R. Clark, Hybrid heuristics for planning lot setups and sizes, Computers and Industrial Engineering 45 (2003) 545–562. [19] A.M. Newman, M. Kuchta, Using aggregation to optimize long-term production planning at an underground mine, European Journal of Operational Research 176 (2007) 1205–1218. [20] M.A. Lejeune, A variable neighborhood decomposition search method for supply chain management planning problems, European Journal of Operational Research 175 (2006) 959–976. [21] K.F. Doerner, V. Schmid, Survey: Matheuristics for rich vehicle routing problems Hybrid Metaheuristics 2010, Lecture Notes in Computer Science 6373 (2010) 206–221. [22] J. Lazic´, S. Hanafi, N. Mladenovic´, D. Urosevic´, Variable neighbourhood decomposition search for 0–1 mixed integer programs, Computers and Operations Research 37 (2010) 1055–1067.