A new formulation to the shortest path problem with time windows and capacity constraints

R. Dondo

INTEC (Universidad Nacional del Litoral - CONICET). Güemes 3450 - (3000) Santa Fe - República Argentina Tel. +54 342 4559174/77 ; Fax: +54 342 4550944. E-mail: rdondo@santafe-conicet.gov.ar

Abstract— The Shortest-path problem with time-windows and capacity constraints (SPPTWCC) is a problem used for solving vehicle-routing and crew-scheduling applications. The SPPTWCC occurs as a sub-problem used to implicitly generate the set of all feasible routes and schedules in the column-generation formulation of the vehicle routing problem with time windows (VRPTW) and its variations. The problem is NP-hard in the strong sense. Classical solution approaches are based on a non-elementary shortest-path problem with resource constraints using dynamic-programming labeling algorithms. In this way, numerous label-setting algorithms have been developed. Contrarily to this approach and with the aim to obtain elemental and optimal solutions, we propose a new mixed integer-linear formulation to the SPPTWCC. Some valid inequalities that can be used to strengthen the linear relaxation of the SPPTWCC are also proposed. Numerical experiments on some VRPTW instances taken from Solomon's benchmark problems show that (near) optimal solutions are easily obtained in spite of the considerable problem size. Also the number of generated columns is kept at a very low level.

Keywords— Shortest Path Problem; MILP Formulation; Column Generation; Vehicle Routing.

I. INTRODUCTION

The shortest-path problem with time-windows and capacity constraints is a problem widely used for formulating vehicle-routing and crew-scheduling applications (Desaulniers et al., 1998). The SPPTWCC consists of finding the shortest path from a source-node to some nodes of a network (while fulfilling timing and capacity constraints) that ends in a sink node. The term "shortest path" should be carefully interpreted: given costs associated to arcs and prices associated to the nodes, the aim is to find the least-cost path from the source node to the sink node. The SPPTWCC occurs as a sub-problem used to implicitly generate the set of all feasible routes in the column-generation formulation of the vehicle routing problem with time windows (VRPTW) and its variants (Cordeau et al., 2002). It is NP-hard in the strong sense. In the VRPTW, the source and sink nodes are usually located in the same place. This place is commonly named "the depot". For n-depot routing problems, n source/sink pairs placed in the same location are usually defined. We may relax the problem to consider variants with source and sink nodes placed on different locations. This work deals just with the single-depot case but the more general cases are straight forward.

The SPPTWCC is also a problem with an economic meaning. I.e. given a set of profits associated to the nodes we aim to choose, at a non-zero cost, a subset of nodes that maximizes our net profit.

Classical solution approaches are based on the non-elementary shortest-path problem with resource constraints, using pseudo-polynomial dynamic programming labeling algorithms. Very refined and complex algorithms of this type have been developed. (See e.g. Houck et al., 1980; Irnich and Desaulniers, 2005 and Irnich and Villenueve, 2006). These algorithms are very effective in generating, in addition to the best route, many solutions per iteration. On the contrary, our purpose is just to obtain the optimal solution to the SPPTWCC. Consequently, we propose a new mixed integer-linear formulation (MILP) to the problem.

This work is organized as follows: Section 2 describes the problem and presents its conventional MILP formulation. Section 3 presents a novel MILP formulation based on global precedence relationships. Advantages and weaknesses of this model are also discussed in this section. In Section 4 several pre-processing and polyhedral techniques are applied to the new formulation in order to improve its numerical resolubility. Numerical examples that arise from the well known Solomon benchmark collection are presented and discussed in Section 5. Finally, the conclusions are outlined in section 6.

II. PROBLEM DEFINITION AND ITS USUAL MATHEMATICAL MODEL

Consider a route-network represented by an undirected graph G{I ∪ p, A } with I = {i₁, i₂, ..., i_n} denoting the set of nodes or customers and p representing a source /sink node called "the depot". Nodes and the depot are connected by a set of arcs A = {(i, j) / i,j ∈ I ∪ p}. Known load and price vectors L = [l₁ , l₂, ..., l_n] and Π = [π₁, π_{2, ...}π_n] are associated to the customer set I. Loads l_i must be collected within a time window [a_i, b_i] on each node i ∈ I. The parameters a_i stands for the earliest possible start-time of the service and the parameters b_i states the latest possible start-time of the service at the node. In addition, travel-costs C = {c_ij} and travel times = {t_ij} are given data for any route segment (i,j) ∈ A. Moreover, the service time on node i is denoted st_i. For each cargo l_i collected node i∈ I, an associated price π_i is accumulated. It is assumed that the triangle inequality is satisfied by the travel costs and travel times, i.e. and . The solution to the SPPTWCC problem must: (1) Maximize the net profit collected from the selected subset of nodes I ^opt⊆ I. This profit is defined as the sum of collected prices minus the cumulated cost incurred by traveling arcs to pick them. (2) The resulting route must start and end on the depot p. (3) The selected nodes must be visited once, so an elemental path is designed. (4) The total collected load must never exceed a given capacity q. (5)The time-length used to collect loads and premiums must be shorter than the maximum allowed working time t^max. (6)The service at every customer site i must start within the specified time window [a_i, b_i]. This problem is usually formulated as follows:

(1)

Subject to

	(2)
	(3)
	(4)
	(5)
	(6)
	(7)
	(8)

while Eq. (1) states the objective function above mentioned, Eq. (2) states a capacity constraint. Constraints (3), (4) and (5) are flow constraints resulting in a path from the depot p to the subset of chosen nodes and back to the depot. Constraints (6) and (7) are timing constraints and constraint (8) limits the routing time to a maximum value t^max. The binary variable x_ij indicates if arc (i,j) ∈ A is used (x_ij = 1) or not (x_ij = 0). The parameter γ_ij is the reduced cost (c_ij - π_i) of using the arc (i,j). While c_ij is a non-negative number, γ_ij can be any real value. As the SPPTWCC is NP-hard, no algorithm with a worst- case running-time bounded by a polynomial in the size of the problem is known. To optimally solve this problem, we have to use an enumerative algorithm such as branch and bound with a worst-case solution-time that is exponential in the size of the input. Note that the formulation given by Eqs. (1) to (8) is, in practice, unsolvable because of the high number of binary variables x_ij. Furthermore, the weak linear relaxations of constraints (6) lead to enormous search trees. Therefore, in the context of column generation algorithms, the most common solution approach to the SPPTWCC is a dynamic-programming label-setting algorithm (See Desrochers et al., 1992). While searching for the optimal route, these procedures generate many non-optimal but useful routes. Very refined and complex algorithms of this type have been developed (Houck et al., 1980; Irnich and Desaulniers, 2005; and Irnich and Villenueve 2006). These methods often lead to problems with thousands of columns that are very hard to solve. Ropke and Cordeau (2009) claim that although the SPPTWCC can initially be solved with heuristic algorithms, whenever the heuristic can no longer produce columns with negative reduced costs, it is necessary to switch to exact methods and always the last pricing problem must be solved to guaranteed optimality in order to prove that the lowest bound is valid. Consequently, optimizing formulations for the SPPTWCC are of crucial importance to prove the optimality of a given solution. This constitutes the main motivation for developing the formulation to the SPPTWCC to be next presented.

III. A REFORMULATION TO THE SPPTWCC

The computational hardness of most combinatorial problems has inspired researchers to develop good formulations that are expected to reduce the size of the enumeration tree and the computation times of these problems. In such a way, we can reformulate the SPPTWCC as follows:

(9)

subject to:

	(10)
	(11)
	(12)
	(13)
	(14)
	(15)

The objective function (9) is expressed as minimization of the difference between the overall travelled distance (CV) and the total quantity of collected prices (Σ _i∈I π_iY_i). Eq. (10) is a capacity constraint equivalent to eq. (2). Eq. (11) is like eq. (3) but states flow constraints using continuous variables. It computes the least travelling costs and times (C_i and T_i) from the depot p to a given node i. Eq. (12) combines and reformulates the information from constraints (4) and (6) in order to sequence nodes. In this way, let us assume that nodes i and j are both in the optimal path (Y_i = Y_j = 1). Then, the relative ordering of nodes i and j becomes determined by the sequencing variable S_ij. In such a case, node j can be a direct/indirect predecessor of node i or viceversa. If node i is visited before j (S_ij = 1), the travel cost from node i to node j (C_j) must always be larger than C_i by at least c_ij. Furthermore, the arrival time at node j (T_j) should be larger than T_i by at least the sum of the traveling time t_ij and the service time (st_i) at the node i. In case node j is visited earlier (S_ij = 0) , the reverse statements hold-on. Eqs. (13) states that the overall traveling cost (CV) must always be larger than the traveling expenses from the depot to any node i (C_i) along the tour by at least the amount c_ip. Also, the total time (TV) required to complete the tour is found by adding the sum of both the service time st_i at node i and the travel time t_ipalong the edge (i,p) to the initial service time at the node last visited i. Since the node last visited is not known beforehand, the eq. (13) must be written for every node i ∈ I. Eqs. (14) and (15) are time-windows and maximum routing time constraints. This formulation differs from the usual formulation to the SPPTWCC, given by Eqs. (1)-(8) in the following main issues:

• It uses a continuous representation for both the time domain (variables TV and T_i for all i ∈ I) and the cost domain (variables CV and C_i for all i ∈ I).
• The new formulation handles node-visit and node-sequencing decisions through different sets of binary variables. The variable Y_i indicates if node i ∈ I belongs to the optimal path (Y_i = 1) or not (Y_i = 0) while the variable S_ijset the generalized precedence relationship between nodes (i, j) ∈ I: i < j if both are in the optimal path (Y_i + Y_j = 1). Constraints (12) become redundant whenever nodes i and/or j are not in the optimal path (Y_i + Y_j ≤ 1). In such a case, the constraint will state that C_i, T_i, C_j and T_j are all larger than a large negative number.
• The reformulation uses the concept of generalized predecessors for sequencing variables. Variable S_ij indicates that node i is visited before (S_ij = 1) or after (S_ij = 0) node j just in case both nodes belong to the optimal tour. This variable indicates both direct and indirect precedence relationships. I.e. if i is a direct predecessor of j, the term (c_i ≥ c_i + c_ij) will be satisfied as equality. If i is a indirect predecessor of j, the inequality sign ">" will become active.

Our formulation is aimed at limiting one flaw of the original one; the high number of 0-1 variables. More precisely we have

O[|I| (|I|-1)] binary variables;

O[|I|] continuous variables;

O[3 + 3|I| + |I| (|I|-1)] constraints

on the conventional formulation and

O[|I| + ½|I| (|I|-1)] binary variables;

O[2|I| + 2] continuous variables;

O[6|I| + 4|I| (|I|-1) + 2] constraints

on the new formulation. So, a considerable saving of 0-1 variables is achieved. These can be reduced to almost a half (if |I| is large) respect to the conventional formulation, at the cost of increasing continuous variables and constraints. Consequently, we have a formulation with fewer binary variables and more constraints. This should be favourable from a resolution point of view and should lead to smaller branch and bound trees. Although the new model uses fewer binary variables, greatly exploits "big-M" type constraints. Thus, it is expected to have a loose continuous relaxation. It is well known that big-M constraints can be difficult for branch and bound solvers since they generate un-tight lower bounds that are crucial for the efficiency of the solver. In general, the larger M in the big-M constraint, the less tight the formulation. Then, the big-M value should be substituted by an upper bound in the terms that activates. I.e, in the equation:

the parameter M_T should be replaced by an upper bound on the term (T_i + st_i + t_ij). In this case the bound is (b_i + st_i + t_ij). In this way, "customized" and as small as possible big-M parameters should be introduced into each constraint. This should lead to a higher tightness of constraints and therefore, the problem might be faster solved. Nevertheless, even with this "strongest" version of inequalities, they are very unfavourable for linear relaxations. Consequently, we have achieved one objective (to reduce the number of 0-1 variables) but we are still unable to provide tight lower bounds to linear relaxations. Since the number of constraints with big-M structure are higher than in the new formulation, the original formulation (Eqs. 1-8) presents even worse bounds. Ways to improve the linear relaxations on the new formulation are presented in the next section.

IV. PRE-PROCESSING AND POLYHEDRAL TECHNIQUES IN THE NEW MODEL

The size and complexity of solved integer problems can be increased when the polyhedral theory is applied. The underlying idea of polyhedral combinatronics is to replace the constraints-set of the integer-programming problem by an alternative convexification of the feasible points and rays of the convex hull conv(S) of the problem (von Hoesel and Aardal, 1996). Then, if we can list the whole set of linear inequalities that defines conv(S) we can solve the integer problem by linear programming. Nevertheless, for most integer problems the minimal number of inequalities necessary to describe the polyhedron is exponential in the number of variables. Thus, it is unrealistic to search for the complete set of inequalities. In addition, if generated within a branch-and-bound tree they can be not valid throughout the entire tree since cuts are generated assuming that certain values are fixed. Therefore, we are interested in a set of constraints conv(S) as small as possible and with no a-priory value-assumptions. Then, one should consider inequalities that define a facet of conv(S). They are the "best possible" in the sense that they cannot be "stronger" without losing some feasible mixed-integer solutions of the original problem. Frequently only a partial set of valid inequalities "located" in the neighbourhood of the optimal solution is useful to reduce solution times (Hoffmann, 2000). It is worth noting that, to make polyhedral methods work well, one important issue is pre-processing. Important elements of pre-processing are to reduce the size of the initial formulation and to reduce the range of constraint coefficients to make instances numerically more efficient.

A. Pre-processing: In pre-processing, the formulation is tightened before the actual optimization. This is done by fixing some variables or by reducing the interval of values a variable can take. This leads to a more compact solution space and consequently to shorter solution times. In this way, to pre-set some sequencing constraints, we define the following useful sets:

Nodes-compatibility relationship sets:

Set of nodes compatible with i ∈ I : A node j is said to be compatible with a reference node i if can be visited either before or after i. This compatibility condition is stated by the following set:

(16)

Set of predecessors of node i ∈ I : A pair of nodes (i,j) is said to be pre-ordered if they must be visited in a certain pre-determined order when time-window constraints are satisfied. For instance, node j is said to be a predecessor of node i if j must be visited before node i. This condition is defined by the following set:

(17)

Set of successors of node i ∈ I : Node j is said to be a successor of node i if j must be visited after node i. Successors of node i are specified by the following set:

(18)

Set of nodes incompatibles with i ∈ I : Nodes (i,j) that cannot be assigned to the same path are called incompatible. The incompatibility condition for nodes j ≠ i is stated by the following set:

(19)

Stronger compatibility relationships can be defined if we use the concept of "immediacy" for defining the following more restrictive sets:

Set of immediate predecessors of node i ∈ I: A node j is said to be an immediate predecessor of node i if it is a predecessor of such a node and in addition, no other node k ∈ I: k ≠ j can be inserted in a path from i toward j. In such a case, the only feasible path connecting i and j is the arc (i, j) ∈ N. Immediate precedence relationships are stated by the following set:

(20)

Immediate successors are defined in the same way:

(21)

Immediately compatible nodes are nodes i and j ∈ I that are compatible and, in addition, no a third node k can be inserted between them, no matter its relative ordering. Immediate compatibility is defined by the following set:

(22)

B. Domain reduction: The narrowing of the range of a given variable is known as domain-reduction. The reduction of the domain of a continuous variable refines the information known about this variable. In the context of the VRPTW and the SPPTWCC, as the triangle inequality holds, the earliest arrival time to a customer can be strengthened by the time used to arrive straight from the depot and the latest arrival time can be strengthened by driving the fastest way to the depot. So, the customer i time-window can initially be strengthened from [a_i, b_i] to [max (a_p + t_pi, a_i), min (b_p - st_i - t_ip, b_i)]. A further reduction can be achieved, in the context of the formulation (9)-(15), by using two of the Desrochers rules (Desrochers et al., 1992). Considering that node i is assigned to the optimal tour, the beginning of the time windows for all successor j ∈ Suc(i) of customer i, and the closing of time windows for all predecessors j ∈ Pre(i) of node i, can be strengthened as follows:

(23.a)

(23.b)

These rules are useful in case a node is assigned to a partial tour before the actual optimisation as i.e. in not-root nodes of a branch and bound tree for the VRPTW. Rules (23) not only narrow time windows. They also allow to move some nodes from sets Com(i) to Pre(i) ∪ Suc(i) and from Pre(i) ∪ Suc(i) to Inc(i). So, the re-use of narrowing rules and the re-definition of compatibility sets until no further changes are achieved lead to easier problems. In such a way, the sequencing constraints (12) now can be split as follows:

(12.b)

(12.c)

(12.d)

(12.e)

So, the use of narrowing rules and precedence relationships can lead to a considerable saving of sequencing variables and sequencing constraints. Furthermore, the remaining sequencing constraints can be greatly simplified.

C. Valid inequalities: The use of information about the structure of the convex hull of feasible solutions has been so far one of the most useful approaches for solving combinatorial problems (Hoffmann, 2000) and will be fully exploited to strength the new formulation in order to tight the lower bounds of a branch and bound procedure. Consequently, we aim to develop specific inequalities necessary in the description of the convex hull of the solution space of the SPPTWCC formulated according Eqs. (9) to (15). If a problem is NP-hard, we cannot expect to have a concise description of the convex hull conv(S) of the feasible solutions. This does not have necessarily a negative impact since what we need is just a good description of the "area around" the optimal solution (Hoffmann, 2000). As the number of valid inequalities grows exponentially on the number of nodes, we should focus only in the most attractive inequalities. Those are the ones that eliminate the greatest possible number of suboptimal configurations within a given path and are easy to compute. Sometimes simple exact rules are enough to provide a number of valid inequalities and we identified the following types:

Valid inequalities based on two-vertexes: These are constraints aimed at exploit information from every pair of nodes i,j ∈ I: i < j. Therefore its computation burden is in the order O(n²). For a given couple of nodes i,j ∈ I: i < j, first let us assume that node j is a successor of node i. Then, the only feasible path involving both nodes to the depot is (i → j → p). In addition, if the earliest possible returning time to the depot is larger than t' = t_v^max - min_{j k∈ I: j ≠ k} (st_k +t_jk), no visit to another node k ∈ I: k ≠ j ≠ i is possible. In such a case, if the cost of returning via j (c_ij - π_j + c_jp) is higher than the cost of returning directly to the depot (c_ip), the partial path (i → j → p) is sub-optimal and does not belong to the optimal path. Consequently, the inequality (Y_i + Y_j) ≤ 1 becomes valid. A "mirror inequality" is valid if the only feasible path involving both nodes to the depot is (j → i → p) and if the cost of returning via i (i.e. c_ji - p_i + c_ip) is higher than the cost of returning directly to the depot (c_ip). Finally, if i,j ∈ I: i < j are compatible nodes, two partial paths would be feasible; (i → j → p) and (j → i → p). Consequently, both the direct way and the indirect way to the depot are feasible. If both indirect ways are more expensive than the direct ones, the inequality (Y_i + Y_j) ≤ 1 becomes valid. These above assertions are mathematically expressed by Eq. (26).

(24)

where b_k' = (b_k - st_k - min _{ij: i ≠ j} t_ij). Similar inequalities can be designed considering the depot as the start node for the partial paths involving nodes i and j. This lead to the following additional valid inequalities:

(25)

Valid inequalities based on three-vertexes: These inequalities could be derived considering every triplet i, j, k ∈ I: i < j < k with a computational burden in the order O(n³). In such a case, the information about time-windows and about node-prices can be exploited.

Valid inequalities based on time-windows information: Let us assume, without loss of generality, that the only feasible partial path between nodes i, j and k is (i → j → k). If, in addition, the minimum path time (i.e. a_i + st_i + t_ij + st_j + t_jk) for going to k via j is bigger than b_k, the three nodes cannot be simultaneously on the same path. Then the inequality (Y_i + Y_j + Y_k) ≤ 2, becomes binding. This is just one case on which the nodes are pre-ordered. For the two remaining cases (i.e. i,j,k ∈ I: i < j < k, j ∈ Pre(i), k ∈ Pre(j) and i,j,k ∈ I: i < j < k, j ∈ Pre(i), k ∈ Suc(i)), the same inequality is binding. They are expressed as follows:

(26)

Furthermore, if two of the three vertexes of a given triplet are compatible, two partial paths must be verified. If both paths are suboptimal, the inequality (Y_i + Y_j + Y_k) ≤ 2 would be binding. I.e. if nodes k and j are compatible and both are successors of node i, two minimum time paths are to be verified. If (a_i + st_i + t_ij + st_j + t_jk) is higher than b_k and if (a_i + st_i + t_ik + st_k + t_kj) is larger than b_j, then inequality (Y_i + Y_j + Y_k) ≤ 2 becomes valid. The same inequality must be written if compatible nodes k and j are both predecessors of node i and if both three-node paths are infeasible. The whole statement is expressed by the following equation:

(27)

Price-based inequalities: Value of prices π_i also can be exploited for identifying suboptimal paths. I.e. let us consider that the only feasible partial path between nodes i, j, k ∈ I is (i → j → k) and that, in addition, no other node can be inserted between them, so j is an immediate successor of node i and k is an immediate successor of node j. Let us now compare the costs of both feasible paths. If the direct path from i to k is cheaper than the indirect path (i → j → k) because (c_ij - π_j + c_jk) > c_ik, then the three nodes must not be in the optimal path. Therefore the inequality (Y_i + Y_j + Y_k) ≤ 2 becomes valid. For the two remaining cases (i.e. i,j,k ∈ I: i < j < k, j ∈ Pre(i), k ∈ Pre(j) and i,j,k ∈ I: i < j < k, j ∈ Pre(i), k ∈ Suc(i)), mirror inequalities are binding. They are derived in the same way that time-window based inequalities and are stated as follows:

(28)

Furthermore, if two nodes are compatible, inequalities like the ones of Eq. (27) would also be binding. They are stated by Eq. (29).

(29)

Inequalities based on immediate compatibility relationships: Let us consider that node i ∈ I belongs to the optimal path. In such a case, just one of all nodes j ∈ IPre(i) may be a predecessor of vertex i and just one node k ∈ ISuc(i) may be its successor. This raises the following two valid inequalities:

(30.a)

(30.b)

V. LIFTING THE INEQUALITIES

Lifting refers to extending valid inequalities from a low dimensional polyhedron to polyhedrons that are valid in higher dimensions. The concept of lifting has been introduced by Gomory (1958). Padberg (1979) developed a sequential lifting procedure for binary programming. After those pioneering works, numerous lifting techniques for a variety of constraint-structures have been proposed (See e.g. Hosel and Aardal, 1996). We will use some established results in order to lift the above presented inequalities. First, let us consider the incompatibility relationship (12.b). As the node i is incompatible with any node j ∈ Inc(i), if i belong to the optimal path, no just the node j ∈ Inc(i) can be inserted on it. Any other node k ∈ Inc(i): k # j incompatible with i will also be forbidden. Then, the following will be a stronger valid constraint:

(31)

This constraint can replace the arrangement of incompatible-couple constraints related to node i that are stated by Eq. (12.b). Eq. (31) is a stronger constraint because it extends inequality (12.b) from a {0, 1}² space to a {0, 1}^|Inc(i)| dimension. It also reduces the total number of constraints of the formulation.

The first constraint of Eq. (24) can be lifted by a similar procedure and the following inequality can be obtained:

(32)

where Ξ = { i,j ∈ I : j ∈ Suc(i), i < j, a_i + st_i + t_ij + st_j > t', c_ij - π_j + c_jp > c_ip}

.

This means that if node i is in the optimal path, not just one but all nodes j that fulfil the condition stated by the first constraint of Eq. (24) must be excluded from the tour. Now, this stronger inequality can replace the former constraint. A mirror equation can be written for the second constraint of Eq. (24). In the same way, the last restriction can be lifted to obtain its following enforced replacement:

(33)

where Ψ = { i,j ∈ I : j ∈ Suc(i), i < j, a_i + st_i + t_ij + st_j > t', c_ij - π_j + c_jp > c_ip, a_j + st_j + t_ij + st_i > t', c_ji - π_i + c_ip > c_jp }. In addition, inequalities that consider the depot as a tour start-point are lifted in the same way and enforced versions of eqs. (25) can be straightforward derived and written. Now, considering the first constraint of eq. (26), we can see that if both nodes i and j are in the optimal path, the least arrival time to node j moves from a_j to max[a_j ; a_i + st_i + t_ij]. Consequently some nodes k that are successors of of node j may become incompatible to the couple (i, j) and the following will be a facet-defining-constraint:

(34)

where Ω = { i,j,k ∈ I: k ∈ Suc(i), a_i + st_i + t_ij + st_j + t_kj > b_k}. This inequality states that successors of node j that cannot be visited later than (a_i + st_i + t_ij + st_j + t_k) are to be excluded from the optimal path just in case nodes i and j belong to it. Conversely, if j is predecessor of node i, the following mirror inequality would become valid:

(35)

where Ω' = {k ∈ I: k ∈ Suc(i): a_j + st_j + t_ji + st_i + t_kj > b_k}.

VI. SOME NUMERICAL EXAMPLES AND DISCUSSION

One of the most prominent vehicle routing problems with side constraints is the VRPTW. This section illustrates the use of the new MILP formulation on a simple branch-and-price algorithm aimed at finding optimal solutions to VRPTW instances. Solomon (1987) benchmark VRPTW instances have attracted numerous researchers to develop exact and heuristic solution procedures and are commonly used as a test-bed for these procedures. The collection of Solomon's 56 problems has been grouped into three categories: C, R and RC. C-class problems feature clustered customers whose time windows have been generated based on known solutions. Locations in R-class problems were randomly generated over a square while RC-class problems comprise a combination of clustered and randomly generated customers. The data set for every category comprises 100 nodes, a central depot, similar vehicle capacities but different time-window distributions. Euclidean distances among customers and traveling times are numerically identical. Furthermore, time windows are hard constraints, service times are independent of customer requirements and the tour duration cannot exceed a maximum value t^max. The objective is the minimization of the total distance. Smaller problems can be generated by selecting the first 25 or 50 nodes of each instance. Benchmark problems of each class are further classified into types "1" and "2", like C1 and C2. Type-1 problems have narrow time windows and small vehicle capacities while type-2 problems feature wider time windows and larger vehicle capacities. In order to evaluate the performance of our SPPTWCC formulation we first solved all R1-type instances with just the first 25 nodes. We selected this group because the different time-windows lead to solutions involving a wide span of solution-shapes. I.e. problem R101 have a solution with numerous trips involving a few nodes per trip while problems R104, R108 and R112 have solutions with fewer tours and many nodes per tour.

In order to "translate" these benchmark problems to the SPPTWCC, we included into the data a price vector Π = [π_1, π_{2, ...}π₂₅].

The vector was obtained by generating columns in a column generation procedure until reaching the optimal lower bound to the problem. Then, the price vector was "frozen" into values obtained in such a way. Also inter-node distances were rounded to the nearest first decimal value. Price-values are reported in Table 1. To test the proposed SPPTWCC model, different configurations of enforced formulations were coded using ILOG OPL Studio 3.7 and all R1-problems with 25 nodes were solved in a 2.8 GHz 1.0 GB RAM Pentium IV PC.

Table 1: Optimal-price values on each Solomon R1-instances comprising the first 25 nodes

We first compared three basic configurations. Configuration 1 uses the new model after applying pre-processing and domain reduction but no valid inequalities at all. Configuration 2 adds the valid inequalities and Configuration 3 replaces them by the "lifted" equations presented in section V. Table 2 summarizes objective function values of relaxed SPPTWCC. These relaxed problems were solved considering variables Y_i and S_ij as continuous variables bounded by the interval [0, 1]. As prices reported on Table 1 were computed considering variables Y_i and S_ij as binary variables, the re laxation of their integrality leaded to negative objective function values. It can be observed here the effect of valid inequalities in the relaxation of the problems. The shrinkage of the objective function value shows that these additional constraints might lead to smaller branch and price trees and therefore to shorter CPU times. Nevertheless, the relative reduction of the (negative) objective function value seems to depend on the problem morphology. These claims were tested by solving the R1-type VRPTW instances using the three different configurations of the SPPTWCC as slave problems within a basic branch and bound procedure (Barnhart et al., 2000). Results are presented in Table 3. We can note the following patterns:

• The effect of valid inequalities in problems constrained by tight time-windows (i.e. R101, R105 and R109) is null or slightly negative.
• The positive effect of valid inequalities is mostly seen in problems constrained by weak time-windows and can be considerable. See solution times for the R104, R108 and R-112 problems.
• The global effect of valid inequalities is positive because tightly constrained problems (where the effect of inequalities is slightly negative) are much easier to solve than almost no-constrained problems on which the accelerating effect of valid inequalities can be substantially positive.
• The accelerating effect of lifted inequalities respect to un-lifted valid inequalities is minor.

Table 2: Objective function values for different SPPTWCC formulations with relaxed assignment (Y_i) and sequencing (S_ij) desicion-variables

Table 3: Resolution times (in seconds) and objective function values for the Solomon R1-instances comprising the first 25 nodes.

The Table 4 shows the number of columns generated while using the MILP formulations of the SPPTWCC. Interestingly, this number remains more or less constant for the whole problem series and is dramatically smaller than in algorithms that use conventional label setting procedures as routes-generators (See Larsen, 2000). In spite of this, CPU times are larger in our algorithms. This is because label-setting algorithms generate many columns per iteration. In spite of this initial disadvantage, our procedure shows the following properties:

Table 4: Columns generated while solving R1-instances (25 nodes) with the three configurations.

• It is very flexible since it can handle almost any structure of constraints while providing good performances on known benchmark problems. Then, this MILP model can be a cornerstone to model more complex problems, like i.e. the pick-up and delivery problem and its realistic variations as multi-commodity vehicle routing problems, routing problems with multiple time-windows, etc.
• It should lead to smaller branch-and-price trees, thus erasing its disadvantages respect to label setting algorithms in larger problems.

To test the last claim we solved the R1 series involving the first 50 customers. The Table 5 shows resolution data for this problem set. Interestingly, if we compare this data with data from problems involving just 25 nodes, we can't detect an explosion of CPU times and of generated columns.

Table 5: Overview of solution parameters for the R1-Solomon examples with 50 nodes

The Table 6 shows resolution data en the R1-series involving 100 nodes. Although the number of columns grow with the problem size, this growing seems to be more or less bounded. I.e. the average number of columns is around 200 for the 25 node-series, 400 for the 50 nodes series and around 1000 for the 100 nodes series.

Table 6: Overview of solution parameters for the R1-Solomon examples with 100 nodes

VII. CONCLUSIONS

In this work, we developed a new MILP formulation for the SPPTWCC that is useful in the context of column generation methods designed to solve vehicle routing problems. It is flexible since it can handle almost any structure of constraints while providing good performances on known benchmark problems. Thus, this MILP model can be used as a cornerstone to model more complex and more realistic problems. We have also developed new valid inequalities tailored for such a model and introduced its "lifted" equivalences. The numerical research demonstrates that these inequalities are useful in substantially reducing CPU times.

We think that those developments constitute a starting point for the development of slave route-generator problems that can model a wide arrange of realistic and complex problems. We can see a remarkably constant level of generated columns on problems of a given size but with different solution-morphologies. This obviously means that columns of more complex problems demand higher generation times. Nevertheless, CPU times used to solve small problems with our formulation are higher than CPU times reported in the literature for algorithms based on label-setting procedures. This is because these procedures are able to generate many columns per iteration. Nevertheless, optimality for these columns maybe hard to prove. It is important to highlight that the objective of the proposal is not to compete with algorithms based on label setting procedures but to complement them with this alternative model. It seems that a column generation algorithm involving calling to both types of "routes generators" seems to be the most efficient option. We should first use the label-setting algorithm and then, whenever the branching mechanism demands a few but hard to find columns we should switch to the MILP formulation. This, of course, is an issue for future research. Another area of continuing research should include the design of branching methods tailored to exploit the structure of the new formulation. It seems that branching on assignment variables Y_i is an effective method to keep the search tree bounded. The continuing development of valid inequalities also deserves additional work.

REFERENCES
1. Barnhart, C., E. Johnson, G. Nemhauser, M. Savelsbergh and P. Vance, "Branch and Price. Column Generation for Solving Huge Integer Programs," Op. Res., 48, 316-329 (2000).         [ Links ]
2. Cordeau, J.., G. Desaulniers, J. Desrosiers, M. Solomon, and F. Soumis, "VRP with time windows," The Vehicle Routing Problem, P. Toth, D. Vigo, eds. SIAM, Philadelphia, PA., 7, 155-194 (2002).         [ Links ]
3. Desaulniers, G., J. Desrosiers, I. Ioachim, M. Solomon, F. Soumis and D. Villeneuve, "A unified framework for deterministic time constrained vehicle routing and crew scheduling problems," Fleet Management and Logistics, T. Crainic, G. Laporte, eds. Kluwer Academic Publisher, Boston, MA., 3, 57-93 (1998).         [ Links ]
4. Desrochers, M., J. Desrosiers and M. Solomon, "A new optimization algorithm for the vehicle routing problem with time windows," Op. Res., 40, 342-354 (1992).         [ Links ]
5. Gomory, R., "Outline of an algorithm for integer solutions to linear programs," Bull. AMS, 64, 275-278 (1958).         [ Links ]
6. Hoffmann, K., "Combinatorial Optimization: History and Future Challenges," J of Appl. and Comp. Math., 124, 341-360 (2000).         [ Links ]
7. Houck, D, J. Picard, M. Queyranne and R. Vemuganti, "The travelling salesman problem as a constrained shortest path problem: theory and computational experience," Op. Res., 17, 93-109 (1980).         [ Links ]
8. Irnich, S and G. Desaulniers, "Shortest path problems with resource constraints," Column Generation, In G. Desaulniers, J. Desrosiers and M. Solomon, editors. Chapter 2, 33-65 (2005).         [ Links ]
9. Irnich, S and D. Villenueve, "The Shortest-Path Problem with Resource Constraints and k-Cycle Elimination for k > 3. Improving a branch and price algorithm for the VRPTW. INFORMS," J. of Computing, 18, 391-406 (2006).         [ Links ]
10. Larsen, J., The Dynamic Vehicle Routing Problem, Ph.D. Thesis. Technical University of Denmark (2001).         [ Links ]
11. Padberg, M., "A note on 0-1 programming," Operations. Research, 23, 833-877 (1979).         [ Links ]
12. Ropke, S and J. Cordeau, "Branch-and-Cut-and-Price for the Pickup and Delivery Problem with Time Windows," Transp. Sci., 43, 267-286 (2009).         [ Links ]
13. Solomon, M., "Algorithms for the vehicle routing and scheduling with time window constraints," Op. Res., 15, 254-265 (1987).         [ Links ]
14. van Hoesel, C.P.M. and K.I. Aardal, "Polyhedral techniques in combinatorial optimization I: theory," Statistica neerlandica: Orgaan van de Vereniging voor Statistiek, 50, 3-26 (1996).         [ Links ]

Received: June 27, 2011.
Accepted: October 24, 2011.
Recommended by subject editor: Pedro de Alcântara Pessôa