Veronica Dal Sasso

Luigi De Giovanni

Martine Labbé

Strengthened Formulations and Valid Inequalities for Single Delay Management in Public Transportation

Keywords: Delay Management, Mixed Integer Linear Programming, Polyhedral Study, Branch-and-cut

The Delay Management Problem arises in Public Transportation networks, often characterized by the necessity of connections between different vehicles. The attractiveness of Public Transportation networks is strongly related to the reliability of connections, which can be missed when delays or other unpredictable events occur. Given a single initial delay at one node of the network, the Delay Management Problem is to determine which vehicles have to wait for the delayed ones, with the aim of minimizing the dissatisfaction of the passengers. In this paper, we present strengthened mixed integer linear programming formulations and new families of valid inequalities. The implementation of branch-and-cut methods and tests on a benchmark of instances taken from real networks show the potential of the proposed formulations and cuts.

Introduction

Intermodal public transportation networks are composed of different means of transport such as trains, buses and metros. In order to encourage the use of public transports, it is essential to offer a good service. It is however impossible to provide direct routes for all origin-destination pairs, thus passengers often have to change vehicles. To minimize the inconvenience for the users, that we can associate to the total amount of time spent travelling, usually connections between vehicles of different lines are considered, and timetables are designed in such a way that the following vehicle in a connection departs shortly after the arrival of the preceding one. This choice has a drawback, that comes out when any vehicle is delayed for some reason. To overcome this issue, it is possible to 2 Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!) allow some extra time in the schedule so that small delays can be recovered on the way. However, some longer delays are unavoidable, resulting in these two alternatives: either the vehicle that follows in the route waits for the delayed one (the connection is maintained), or it departs on time (the connection is dropped). In the former case the connection is not missed but the delay spreads over the network. In the latter case, only the passengers on the feeding vehicle will arrive at their destination later than planned. It can be easily understood that the best solution to minimize the dissatisfaction of passengers will be to maintain a subset of connections and to drop the remaining ones.

In this paper, we deal with a particular delay management problem as initially stated in [START_REF] Schöbel | A model for the Delay Management Problem based on mixed-integer-programming[END_REF], where a single initial delay is considered and known in advance (offline delay management), and passenger re-routing is not allowed, meaning that passengers missing a connection will wait for the next vehicle of the same line, instead of looking for alternative paths to their destinations. Moreover, vehicles in the next time period are assumed to be on time, so that the next vehicle will arrive after a known time interval, which is the same for all lines. Finally, delay propagation does not take capacity constraints into account. Under these assumptions, the Delay Management Problem (DMP) consists in finding the wait-depart decisions that minimize the total inconvenience over customers, measured by the sum of all passenger delays when they reach their destinations. This problem has been formulated by [START_REF] Schöbel | A model for the Delay Management Problem based on mixed-integer-programming[END_REF] as a Mixed Integer Linear Programming (MILP) model based on an event-activity network representation. Further formulations and improvements are presented in [START_REF] Schöbel | Integer programming approaches for solving the delay management problem[END_REF]. In [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] two equivalent models are presented, that progressively reduce the number of variables. The first one, called MILP1, is based on a simplified event-activity network called Arr-arr graph, while the second one, called MILP2, is obtained from MILP1 by Fourier-Motzkin elimination. MILP2 reduces the number of variables at the cost of a substantial increase in the number of constraints, so that a row-generation approach is proposed to solve the model. Some improvements to MILP1 and a preliminary polyhedral study are presented in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF], where valid inequalities for both MILP1 and MILP2 are deduced from the results on the mixed 0-1 knapsack polytope by [START_REF] Marchand | The 0-1 knapsack problem with a single continuous variable[END_REF], and used to strengthen the linear relaxation.

Among other works on DMP, we cite some initial study based on simulation (Suhl andMellouli 1999, Suhl et al. 2001) or max-plus algebra [START_REF] De Vries | On the max-algebraic models for transportation networks[END_REF][START_REF] Goverde | Computer in railways VI , chap. The max-plus algebra approach to railway timetable[END_REF]. The computational complexity of DMP is discussed in [START_REF] Gatto | The computational complexity of delay management[END_REF]. On-line versions of the problem have been considered by [START_REF] Gatto | On-line delay management on a single train line[END_REF], [START_REF] Berger | Online railway delay management: hardness, simulation and computation[END_REF] and by [START_REF] Adenzo-Diaz | On-line timetable re-scheduling in regional train services[END_REF]. In the latter the wait/depart strategy aims at minimizing the impact of delays on vehicle schedules, rather than passengers' discomfort. Further extensions to more general DMPs have been also presented in literature. The effect of network capacity restrictions in delay propagation are considered, e.g., in [START_REF] Dollevoet | Delay management including capacities of stations[END_REF], where stations have a limited number of platforms with impact on wait/depart decisions, or in [START_REF] Schöbel | Capacity constraints in delay management[END_REF], where DMP is integrated with decisions on optimal train-to-platform reallocation. Problems arising from a limited number of tracks are addressed in [START_REF] Schachtebeck | To wait or not to wait -and who goes first? delay management with priority decisions[END_REF]. A bicriterial approach, in which both the total vehicles' delay and the number of passengers that miss a connection are minimized, is presented in [START_REF] Ginkel | To wait or not to wait? the bicriteria delay management problem in public transportation[END_REF]. Another extension to DMP considers the re-routing of passengers that may choose different routes to reach their destination in case of delays [START_REF] Dollevoet | Delay management with re-routing of passengers[END_REF][START_REF] Berger | Online railway delay management: hardness, simulation and computation[END_REF][START_REF] Schmidt | Simultaneous optimization of delay management decisions and passenger routes[END_REF]. A different approach is investigated in [START_REF] Liebchen | Computing delay resistant railway timetables[END_REF] and consists in building timetables that achieve a good trade-off between robustness against delays and increase of the total travel time. [START_REF] Corman | Integrating train scheduling and delay management in real-time railway traffic control[END_REF] integrate delay management and train scheduling, and develop a model and fast heuristics to minimize the time passengers spend in the system. [START_REF] Schön | A stochastic dynamic programming approach for delay management of a single train line[END_REF] introduce a stochastic version of DMP on a single line where delays are affected by uncertainty, and propose a dynamic programming approach to minimize the total passengers' delay at destination. The aim of this paper is to improve the formulations and extend the theoretical results in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF], as well as to integrate them into a branch-and-cut algorithm to solve DMP.

Results from previous works are summarized in Section 1, where we report the models MILP1 and MILP2 and the valid inequalities stated in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF]. In Section 2 we improve the formulations for DMP by providing tighter values for models' coefficients and by reducing the number of variables and constraints for MILP2. In Section 3 we propose new facets for the mixed 0-1 knapsack polytope associated to each constraint of MILP2. Related inequalities represent a new class of cuts for MILP1 and MILP2 that generalizes and dominates the one proposed in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF]. In order to verify the effectiveness of the proposed models and cuts, in addition to the linear relaxation which was briefly investigated in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF], two alternative branch-and-cut algorithms for DMP have been implemented, based on the improved MILP1 and MILP2: Section 4 describes exact and heuristic separation procedures, and Section 5 reports on the computational results obtained for real instances taken from literature: a first set representing a portion of the Bruxelles' intermodal public transportation network, and other instances inspired by the railway system in the Randstad, in the Netherlands [START_REF] Dollevoet | Delay management including capacities of stations[END_REF]. Some final remarks and hints for future research are given in the concluding Section 6.

Basic formulations and valid inequalities

We summarize the two MILP models and the valid inequalities presented in [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] and De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF].

The transportation network is represented by an Arr-Arr graph (K, A R ∪ A C ), that is an eventactivity network [START_REF] Elmaghraby | Activity Networks[END_REF] in which each node k ∈ K represents the arrival of one vehicle Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!) at one station (events) and there are two types of arcs (activities): direct-ride arcs (A R ), linking the arrivals of the same vehicle at two consecutive stations of the same line, and connection arcs (A C ), linking the arrivals of two vehicles of different lines among which a connection exists (a passenger transfer activity is included). Let slack ij be the time that can be recovered along the arc (i, j), that is computed as the difference between the scheduled arrival times (according to timetable) and the minimum necessary time related to the riding and/or transfer time corresponding to the arc itself. Passenger routes correspond to paths on the Arr-Arr graph, and the corresponding set is denoted by P . Each path p ∈ P is a sequence of direct-ride and connection arcs from an origin a p to a destination v p . Let A p be the set of arcs of path p ∈ P ,

A C p ⊆ A C (resp. A R p ⊆ A R
) the set of connection (resp. direct-ride) arcs in p, and c p the number of passengers traveling on p. Finally, let T be the inter-arrival time of two vehicles of the same line at the same station, which we suppose to be the same for all the lines, and D the initial delay.

Model MILP1 considers the following decision variables:

• d k that represents the delay at node k ∈ K,

• z p = 1 if all connections on path p ∈ P are maintained 0 otherwise,

• u p = d vp if all connections on path p ∈ P are maintained 0 otherwise. Assuming, without loss of generality, that the initial delay occurs at node 1, the MILP1 formulation is as follows [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF][START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF]):

(MILP1) min p∈P c p [(1 -z p )T + u p ] (1)
s.t. d 1 ≥ D (2) d j ≥ d i -slack ij ∀ (i, j) ∈ A R (3) d j ≥ d i -slack ij -M ij (1 -z p ) ∀ p ∈ P, (i, j) ∈ A C p (4) u p ≥ d vp -D(1 -z p ) ∀ p ∈ P (5) d i ≥ 0 ∀ i ∈ K (6) u p ≥ 0 ∀ p ∈ P (7) z p ∈ {0, 1} ∀ p ∈ P (8)
where

M ij = D -slack ij .
The objective function (1) minimizes the total delay suffered by all passengers on all the paths.

Each delay is weighted by u p if all connections on path p are maintained, T otherwise. In fact, since we do not allow here for passenger re-routing, people missing a connection will have to wait for the next vehicle of the same line. The initial delay is set by constraint (2). Constraints (3) describe how the delay spreads along direct-ride arcs, while constraints (4) involve connection arcs: if z p = 1 they describe how the delay spreads, otherwise they are redundant (the delay does not spread).

Constraints (5) set the value of u p according to its definition. The last constraints set variables' bounds and domains.

Observing that variables d k do not appear in the objective function, Fourier-Motzkin elimination is used to project the feasible region on the space of the variables z p and u p . To this end, the set Π of all the paths in the Arr-Arr graph is defined. Notice that P ⊂ Π and that c p = 0 for all p ∈ Π \ P .

This yields the following formulation [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF][START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF]):

(MILP2*) min p∈ Π c p [(1 -z p )T + u p ] (9) s.t. u p ≥ Dz p - (i,j)∈Ap slack ij - (i,j)∈A C p M ij (1 -zi,j ) ∀ p, p ∈ Π : a p = 1, v p = v p (10) u p ≥ 0 ∀ p ∈ Π (11) z p ∈ {0, 1} ∀ p ∈ Π. ( 12
)
where zi,j = max p∈ Π:(i,j)∈A C p z p . Constraints (10) allow determining the final delay on a path p by linking it to the delay that spreads through any of the paths p ∈ Π starting at the delayed node 1 and ending in the same node as p . In fact, these constraints state that if all connections on path p and on at least a path for each connection in p are maintained, then the delay u p at the end of path p will be at least the initial delay minus the time that can be recovered along p (in this case, by the objective function, u p = u p ). If any of the connections is not maintained, these constraints are redundant.

Due to the definition of variables z, constraints (10) are not linear. The constraint is linearized by considering, in the last term, the sum over each possible ordered set of paths (p 1 , p 2 . . . p

|A C p | )
where the t-th path must contain the t-th connection of path p. Since we add a constraint for each possible such subset of paths, the constraint corresponding to the maximum value of z for each connection is also included in the formulation, as required by (10), together with other redundant constraints. This is formalized by the following model, corresponding to MILP2 [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF][START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF], where P p t ⊆ Π denotes the set of all the paths containing the t-th connection arc of p, and α(p, i, j) = t if and only if (i, j) is the t-th connection arc of p:

(MILP2) min p∈ Π c p [(1 -z p )T + u p ] (13)
Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

s.t. u p ≥ Dz p - (i,j)∈Ap slack ij - (i,j)∈A C p M ij 1 -z p α(p,i,j) ∀ p , p ∈ Π : a p = 1, v p = v p , p t ∈ P p t , t = 1 . . . |A C p | (14) u p ≥ 0 ∀ p ∈ Π (15) z p ∈ {0, 1} ∀ p ∈ Π. ( 16
)
These formulations can be improved, as shown in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF], by special classes of cover cuts derived from the ones defined for the Mixed 0-1 Knapsack problem [START_REF] Marchand | The 0-1 knapsack problem with a single continuous variable[END_REF]. We recall that a cover cut identifies a set of binary variables that cannot be all equal to 1 at the same time. The Mixed 0-1 Knapsack set is defined as follows:

Y = {(y, s) ∈ B n × R + : j∈J a j y j ≤ b + s}. ( 17
)
where J = {1, . . . , n}, a j > 0, b > 0.

A single constraint of type ( 14) is obtained by fixing p and p and the combination of paths

Y f DM P = (z p , z p 1 , z p 2 . . . z p |A C p | , u p ) ∈ {0, 1} |A C p |+1 × R + s.t. Dz p + (i,j)∈A C p M ij z p α(p,i,j) ≤ (i,j)∈Ap slack ij + (i,j)∈A C p M ij + u p
and, among other polyhedral results on the polytope conv(Y f DM P ), the following is proven.

Proposition 1 (De Giovanni et al. ( 2014)). Let p, p , p 1 , . . . , p |A C p | be paths defined as in ( 14) and let p i = p j and p i = p for each i, j ∈ {1, . . . , |A C p |}. Then:

• inequality (14) is a facet of conv(Y f DM P ) if and only if slack ij = 0, ∀ (i, j) ∈ A p .
• the following inequality defines a facet of conv(Y f DM P ) if and only if D > (i,j)∈Ap slack ij :

z p + (i,j)∈A C p z p α(p,i,j) ≤ |A C p | + u p D -(i,j)∈Ap slack ij . ( 18
)
Inequalities ( 18) say that, if we consider two paths p and p towards the same final node v p and we want to keep the final delay d v p = u p = 0, we need to drop at least a connection along p: indeed, by definition, v p = v p so that u p = u p = 0; moreover, the second term on the right hand side is non negative, as we are not able to recover all the delay along p, having D > (i,j)∈Ap slack ij . If u p > 0, in order to satisfy ( 18), either we drop some connection on path p (z p = 0), or we still need to drop at least one connection on p, unless we consider for p and p a large enough delay For the case when less than |A C p | + 1 distinct paths are involved by ( 14), De Giovanni et al. ( 2014) observe that inequalities (18) are still valid for the related Mixed 0-1 Knapsack set and, consequently, MILP1 and MILP2. In this paper, we will focus on this case, showing that (18) are not facet defining and introducing a class of more general facets.

d v p = u p ≥ D -(i,j)∈Ap slack ij ,

Formulation refinements

The formulations MILP1 and MILP2 can be further improved by choosing better values for slack ij and M ij constants and by refining the set of paths to be taken into account for constraints ( 10) and ( 14).

Improving model coefficients

Constraints (4) use the bound M ij = Dslack ij to force the RHS to a non positive value in case at least one connection on p is dropped (z p = 0), thus making the constraints redundant. In fact, we can assume that the delay at each node (variables d k ) is not greater than the initial delay D.

The constraints can be tightened by computing, for each node k ∈ K, a better bound M k for the maximum delay that we can observe at k provided that all the connections are maintained, taking advantage of the extra time allowed for each activity, represented by the parameters slack ij on the arcs of the Arr-arr graph. We define

M k = max{0, D -L 1k }
where L 1k is the length of the shortest path in the Arr-Arr graph from the delayed node 1 to k, using slack ij as arc weights. Thus, M k represents the maximum delay that can be recorded at node k. Let us now consider an arc (i, j): since the delay on i is bounded by M i , and since vehicles cannot depart before the scheduled time (d j ≥ 0), the maximum delay that can be recovered on (i, j) is

σ ij = min{slack ij , M i } 1 2 3 4 Figure 1
Example: necessary paths to correctly spread the delay.

We thus obtain the following valid formulation MILP1' for DMP:

(MILP1') min p∈P c p [(1 -z p )T + u p ] (19) s.t. d 1 ≥ D (20) d j ≥ d i -σ ij ∀ (i, j) ∈ A R (21) d j ≥ d i -σ ij z p -M i (1 -z p ) ∀ p ∈ P, (i, j) ∈ A C p (22) u p ≥ d vp -M vp (1 -z p ) ∀ p ∈ P (23)
d i ≥ 0 ∀ i ∈ K (24) u p ≥ 0 ∀ p ∈ P (25) z p ∈ {0, 1} ∀ p ∈ P ( 26
)
where constraints ( 22) and ( 23) replace, respectively, (4) and (5) of MILP1. Note that, again, constraints ( 22) and ( 23) are redundant if z p = 0 (the RHS is non positive), while they correctly set the value of the node delays in case the initial delay spreads through maintained connections (z p = 1).

Choosing the paths

Model MILP2 has been defined by [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] after extending the set P of passenger paths to Π, including all the possible paths on the Arr-Arr graph: this was exploited during the Fourier-Motzkin elimination steps leading to MILP2, to prove that the delay correctly spreads. We remark that MILP1 formulation, which is defined on P ⊂ Π, is still valid, as the delay correctly spreads thanks to constraints (3) and ( 4), which directly involve the arcs of the Arr-Arr graph. On the contrary, using P instead of Π in MILP2 is not valid, as shown by the following example.

Let us consider the Arr-Arr graph in Figure 1, with four nodes 1, 2, 3, and 4, two direct-ride arcs

(1,2) and (2,3), and a connection arc (3,4). Let 1 be the delayed node and let P contain only the path 2-3-4. Since no path starts at the delayed node, the feasible region of MILP2 would result in simply the non-negativity constraints: the optimal solution would set z variable to 1 with objective value 0, which denotes that the delay does not correctly spread.

The example suggests that additional paths starting at the delayed node have to be taken into account. Instead of using the set Π of all the paths as in [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF], we define a smaller set as follows. Let pv be the path from the delayed node 1 to node v without connection arcs (i.e.

A C pv = ∅). Note that pv may not exist and, in case it exists, it is unique. Let Θ v be defined as follows

Θ v = (i,j)∈ pv σ ij if pv exists +∞ otherwise
We define, for each node k ∈ K, the following subset:

Π k = {p k } ∪ {q ∈ Π | a q = 1 ∧ v q = k ∧ A C q ⊆ A C P ∧ (i,j)∈A R q σ ij < Θ k } (27)
where A C P = p∈P A C p is the set of all connection arcs included in at least one passenger path. Each Π k includes a subset of relevant paths starting at the delayed node, ending at some actual passenger's destination node k and propagating a large-enough delay: the following proposition shows that a valid formulation for DMP can be obtained by writing constraints (10) or ( 14) with reference to the sets Π k , instead of Π.

Proposition 2. The following formulation MILP2' is valid for DMP and the value of the linear relaxation is the same as MILP1':

(MILP2') min p∈P c p [(1 -z p )T + u p ] (28) s.t. u p ≥ D - (i,j)∈A R p σ ij - (i,j)∈A C p σ ij z p α(p,i,j) + M i 1 -z p α(p,i,j) -M v p (1 -z p ) ∀p ∈ P, p ∈ Π v p , (p 1 , p 2 . . . p |A C p | ) : p 1 ∈ P p 1 , p 2 ∈ P p 2 . . . p |A C p | ∈ P p |A C p | (29)
u p ≥ 0 ∀ p ∈ P (30) z p ∈ {0, 1} ∀ p ∈ P. ( 31
)
The proof is quite technical and it is reported in Appendix A: it reviews the steps of the Fourier-Motzkin elimination provided by [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] in the light of making minimal the set of paths that guarantees that the initial delay correctly spreads through the Arr-Arr graph, which leads to definition (27).

Strengthened valid inequalities

Formulations MILP1' and MILP2' can be tightened using the following inequalities:

z p + (i,j)∈A C p z p α(p,i,j) ≤ |A C p | + u p D -(i,j)∈Ap σ ij . ( 32
)
Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

They directly correspond to (18) and, by the same argument used in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF] to prove Proposition 1, they represent facets of the mixed 0-1 knapsack set associated to each constraint [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] under the conditions D > (i,j)∈Ap σ ij , and, for each i, j ∈ {1, . . . , |A C p |}, p i = p j and p i = p . In this case, (4) involves |A C p | + 1 distinct variables, that is, the paths selected for each connection of p are different from each others and from p .

In this section, we will propose a new class of valid inequalities that dominates (32) in the more general case when just the condition D > (i,j)∈Ap σ ij holds, that is, constraints [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] involve

|A C p | + 1 distinct
u p ≥ D -M vp (1 -z p ) - (i,j)∈A R p σ ij -M a -M c + (M a -σ ab )z p + (M c -σ cd )z p that is (M vp + M a -σ ab + M c -σ cd )z p ≤ (i,j)∈A R p σ ij + (M vp + M a + M c ) -D + u p
The corresponding cut of type (32) is

z p + z p + z p ≤ 2 + u p D -(i,j)∈Ap σ ij .
which is valid for MILP1' and MILP2'. However, we will see in this section that the cut is dominated by the inequality:

z p ≤ u p D -(i,j)∈Ap σ ij .
This scenario occurs every time p ∈ P p i , or P p i ∩ P p j = ∅, for some i, j = 1 . . . |A C p |, i = j. Let us fix one constraint of type [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] and consider the related paths p ∈ P , p ∈ Π v p and tuple

(z p 1 , z p 2 , . . . , z p |A C p |
), where p t ∈ P p t , ∀ t = 1, . . . , |A C p |. Let W be the set of distinct passenger paths p whose associated variable z p appears in the fixed constraint [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF], path p excluded. We denote with w 1 , . . . , w |W | the elements of W . We can define the following Mixed 0-1 Knapsack set

Y f DM P = (z p , z w 1 , z w 2 , . . . ,z w |W | , u p ) ∈ {0, 1} |W |+1 × R + , w∈W ∪{p } γ w z w ≤ (i,j)∈A R p σ ij + (i,j)∈A C p M i + M v p -D + u p ( 33
)
11 where coefficients γ w are defined as the sum of the coefficients in ( 29) related to the same variable z w . More formally, recalling the definition of α(p, i, j), we denote, for each w ∈ W ∪ {p }, the set

A C p (w) = {(i, j) ∈ A C p | p α(p,i,j)
= w}, and we have

γ p = (i,j)∈A C p (p ) (M i -σ ij ) + M v p γ w = (i,j)∈A C p (w) (M i -σ ij ) ∀ w ∈ W.
The mixed cover cut related to (33) is

z p + w∈W z w ≤ |W | + u p D -(i,j)∈Ap σ ij (34) It reduces to (32) if |W | = |A C p |
, which is the case studied by De Giovanni et al. ( 2014). In the following, we will prove that (34) is facet defining for conv(Y f DM P ). The proof is based on the following lemmata, whose rather technical proofs can be found in Appendix B.

Lemma 1. For any path p ∈ Π such that a p = 1, M vp ≥ D -(i,j)∈Ap σ ij . Lemma 2. Assuming that D > (i,j)∈Ap σ ij , W ∪ {p } is a minimal cover of Y f DM P and the inequality z p + w∈W z w ≤ |W | (35) defines a facet of conv(Y f DM P | u p =0 ). Lemma 3. If W = ∅, for any w ∈ W ∪ {p }, γ w ≤ (i,j)∈A R p σ ij + (i,j)∈A C p M i + M v p -D.
By Lemmata 1, 2 and 3, we can exploit the following result related to general Mixed 0-1 Knapsack sets Y defined in (17).

Proposition 3 [START_REF] Marchand | The 0-1 knapsack problem with a single continuous variable[END_REF]). Let a j ≤ b, ∀ j ∈ J. If j∈J π j y j ≤ π 0 defines a facet of conv(Y | s=0 ), having π ≥ 0 and j∈J π j > π 0 > 0, then the inequality Proof. We distinguish two cases: (i) W = ∅ and (ii) W = ∅.

In case (i), by Lemma 3 and being |W | > 0, the conditions of Proposition 3 hold. Hence we have to determine 36) is strictly bounded from above by

β = min s>0 s η(s) -|W | η(s) = max z p + w∈W z w s.t. w∈W ∪{p } γ w z w ≤ M vp -D + (i,j)∈Ap σ ij + (i,j)∈A C p (M i -σ ij ) + s (36) z w ∈ {0, 1} ∀w ∈ W ∪ {p }. If 0 < s < D -(i,j)∈Ap σ ij , then the RHS of (
M vp + (i,j)∈A C p (M i -σ ij ) = w∈W ∪{p } γ w .
This means that we cannot consider all variables z w equal to 1, that is, η(s) < |W | + 1. Also, by Lemma 2, W ∪ {p } is a minimal cover for Y f DM P 0 and, since we are interested in s > 0, we have

η(s) = |W |.
If s ≥ D -(i,j)∈Ap σ ij , then the RHS of ( 36) is bounded from below by

M vp + (i,j)∈Ap (M i -σ ij ) = w∈W ∪{p }
γ w so that we can consider all z w variables equal to 1 and η(s

) = |W | + 1. Summarizing η(s) = |W | if 0 < s < D -(i,j)∈Ap σ ij |W | + 1 if s ≥ D -(i,j)∈Ap σ ij and the minimum value β is obtained for s = D -(i,j)∈Ap σ ij , that is, β = D -(i,j)∈Ap σ ij .
Let us now consider case (ii). In this case, and since D > i,j)∈Ap σ ij , inequality (34) becomes

u p ≥   D - (i,j)∈Ap σ ij   z p .
It is valid for Y f DM P , since it always holds if z p = 0 and, if z p = 1, noticing that in case (ii)

A C p (p ) = A C
p , it reduces to the inequality that defines Y f DM P in (33). Moreover, Y f DM P is in R 2 and, since the two affinely independent points z p = u p = 0 and z p = 1, u p = D -(i,j)∈Ap σ ij are tight, (34) defines a facet.

We remark here that inequality ( 32 We have already observed that inequalities (18) can be used as cuts to strengthen both formulations MILP1 and MILP2: by similar arguments, under the condition D > (i,j)∈Ap σ ij , inequalities (34) are valid cuts for MILP1' and MILP2'.

A branch-and-cut algorithm

Models MILP1, MILP2, MILP1' and MILP2' are solved with a branch-and-cut algorithm, that is, a branch-and-bound algorithm where additional cuts are used to strengthen the linear relaxation. In particular, we solve the linear relaxation and, in case the optimal solution has binary variables at a fractional value, we recursively branch by creating two subproblems where one of these variables takes either value 0 or 1. The value of the linear relaxation provides a lower bound that can be used to stop branching whenever the bound is larger than the incumbent integer-feasible solution.

MILP1 and MILP1' can be directly handled by on-the-shelf solvers, whereas MILP2 and MILP2' present a huge number of constraints ( 14) and ( 29) respectively, which suggests a row generation approach to solve their linear relaxation. We start from a partial model with no constraints ( 14) or ( 29), but the one corresponding to the trivial path starting and ending at the delayed node, and we iterate through the following steps:

1. solve the linear relaxation of the partial model; 2. determine one (or more) constraints ( 14) or ( 29) violated by the current solution;

3. if no violated constraints exist, stop; 4. otherwise add them to the partial model and iterate from Step 1.

Step 2 is performed by the constraint separation procedure.

As observed above, cuts (18), resp. (34), can be used to strengthen the linear relaxation of MILP1 and MILP2, resp. MILP1' and MILP2', to obtain tighter lower bounds for each subproblem in the branch-and-bound procedure: this gives rise to the branch-and-cut approach we propose to solve the models. However, the number of the proposed cuts is huge, and it is convenient, from a computational point of view, to generate and add them dynamically. To this end, cut separation procedures are required. Given a solution of the linear relaxation including just a subset of cuts ( 18) or (34), a cut separation procedure determines one or more further cuts violated by the solution itself: they are added to the linear relaxation and the process iterates. Notice that, for MILP2 and MILP2', both constraint and cut separation procedures are used to dynamically add inequalities to the partial models. In the following, we present constraint and cut separation procedures with reference to models MILP1' and MILP2', with improved coefficients: the corresponding procedures for MILP1 and MILP2 are obtained by substituting σ ij with slack i,j .

Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

Constraint separation

As detailed in Appendix A, constraints (29) of MILP2' come from linearizing inequalities

u p ≥ D - (i,j)∈Ap l ij -M v p (1 -z p ) ∀p ∈ P, p ∈ Π v p , (37)
which link the delay at the end of a passenger path p to the delay propagated on a path p from 1 to v p = v p . Inequalities (37) use the terms l ij , representing the amount of delay that can be saved on arc (i, j), computed as a non-linear function of the z variables:

l ij =    min p∈P :(i,j)∈A C p {σ ij z p + M i (1 -z p )} if (i, j) ∈ A C P σ ij if (i, j) ∈ A R . ( 38
)
Since in the separation procedure the solution to the current partial MILP2' is available, we can use the values related to variables z to compute the terms l ij and directly separate (37): if such a violated constraint exists, a violating constraint ( 29) is determined by the same pair of paths p and p , and by the paths obtained choosing, for each connection arc (i, j) ∈ A C p , p α(p,i,j) = arg max p∈P :(i,j)∈A C p z p , according to (38); otherwise, no constraint ( 29) is violated.

Given the optimal solution (z * , u * ) to the linear relaxation of the current partial MILP2', we want to determine two paths p and p such that (37), and hence [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF], is violated, that is,

(i,j)∈Ap l ij < D -u * p -M v p (1 -z * p ).
Noticing that the right term depends only on p , once p is fixed the maximum constraint violation is obtained by selecting the shortest path p ∈ Π v p with l ij as arc weights: this is done by following Procedure A.

Procedure A:

1. For each (i, j) ∈ A R ∪ A C P , compute l ij according to (38); for each (i, j) ∈ A C \ A C P , set l ij = +∞ 2. For each node k ∈ q∈P {v q }, determine d(k) = min p∈Π k (i,j)∈Ap l ij ; let p(k) be the related path 3. For each path p ∈ P :

4. if d(v p ) > u * p + M v p (1 -z * p )
, then inequality (37) is violated for the paths p , p = p(v p ).

Step 1 computes the terms l ij corresponding to the solution at hand, setting to an infinite cost the connection arcs not included in any passenger path (we are not interested in paths including any such arcs, since the delay is not propagated on them). Step 2 computes the shortest path on the Arr-Arr graph from the delayed node 1 to any possible passenger path destination, according to the l ij metric. We notice that this corresponds to propagating the initial delay and determining, for each destination node k, the path propagating the maximum delay, which is at least the one on the direct path pk . As a consequence, p(k) ∈ Π k , since, by ( 27), all paths from 1 to k not included in Π k propagate a delay not greater than the one on pk .

Observing that the Arr-Arr graph is acyclic, Step 2 can be done by a breadth first visit, and the overall computational complexity of Procedure A is

O(|A C ||P | + |A R | + |A C ∪ A R | + |P |) = O(|A C ||P | + |A C ∪ A R |),
hence polynomial in the length of the DMP input.

Cut separation

Given a fractional solution (z * , u * ) of the LP relaxation of MILP1' or MILP2', we want to determine if it violates any of the inequalities (34) that is, we search for paths p ∈ P , p ∈ Π v p and a suitable set W ⊆ P \ p made of paths including the connection arcs of p, such that

|W | - w∈W z * w < z * p - u * p D -(i,j)∈Ap σ ij . ( 39
)
In the following, we present cut separation procedures that fix the pair of paths p and p, where p starts at the delayed node and v p = v p , and determine the set W minimizing the LHS of ( 39), which corresponds to the maximum cut violation.

We remark that inequalities (34) have been derived under the condition D > (i,j)∈Ap σ ij , as otherwise they are not valid for Y f DM P : e.g., the violating point z p = z w = 1, ∀ w ∈ W , u p = 0 belongs to Y f DM P , as can be easily checked. As a consequence, paths p with D ≤ (i,j)∈Ap σ ij are not taken into account.

Preliminary checks. Before applying a cut separation procedure to a pair of paths p ∈ P , p ∈ Π v p , we can check the following conditions under which no violated cuts involving p and p exist. We present this checking by increasing order of computational complexity.

1. A C p = ∅: under this condition, inequality (34) defines a facet of Y f DM P , but it is redundant for the optimal solutions of relaxed MILP1' and MILP2', for any p ∈ P . In fact, in this case, the constraints of the two formulations imply u * p ≥ D -(i,j)∈Ap σ ij and, from optimality,

u * p = u * p . It follows that z * p - u * p D -(i,j)∈Ap σ ij ≤ 0 and (39) cannot hold, as |W | -w∈W z * w ≥ 0. 2. u * p ≥ M v p z * p : in fact, since M v p = M vp ≥ D -(i,j)∈Ap σ ij , in this case, we have z * p - u * p D -(i,j)∈Ap σ ij ≤ z * p - u * p M v p ≤ z * p -z * p = 0
and (39) cannot hold. Notice the this condition includes the case z * p = 0. 3. u * p ≥ D -(i,j)∈Ap σ ij z * p : as above, in this case we have z * p -

u * p D -(i,j)∈Ap σ ij ≤ 0,
guaranteeing that no cut is violated.

If none of the above conditions holds, the following procedures allow separating inequalities (34)

for a given pair of paths p ∈ P , p ∈ Π v p .

Set covering approach. The problem of finding the maximum cut violation, that is, for fixed p and p , minimizing the LHS of (39), can be solved by the following set-covering formulation:

min w∈P \{p } (1 -z * w )y w s.t. w∈P :(i,j)∈A C w y w ≥ 1 ∀ (i, j) ∈ A C p y w ∈ {0, 1} ∀ w ∈ P

Computational results

The branch-and-cut approach described above has been implemented using the solver Cplex 12.4

(IBM 2014) and Python 2.7 (and related Cplex interface) as programming language. In order to separate cuts for both MILP1' and MILP2', we use the heuristic procedures HP1 and HP2, applied to any pair of paths p ∈ P , p ∈ Π v p . The cut separation procedure (implemented using the Cplex featured user-cuts callback ) works at each branch-and-bound node and iterates until no violated inequalities are found. Concerning the solution of MILP2', at each branch-and-bound node we apply both the constraint separation procedure A and the cut separation procedure. In particular, constraints ( 29) are separated only if the current linear relaxation is integer (using the Cplex featured lazy-constraints callback ), since computational results show a remarkable increase in the running times if we look for violated constraints also at fractional solutions. By combining different models and cut separation procedures, we have obtained and tested several algorithm configurations for DMP. In the following, we report on the ones with better performance, which are described in Table 1: each column is devoted to a configuration and shows its name in the first row and, in the following rows, the model solved by the configuration, the constraint defining the feasible region and the related separation procedure (for models MILP2 and MILP2'), and, in the last row, if cuts (34) are added to strengthen the formulation together with, in this case, the used separation procedure (HP1 or HP2).

The test benchmark is made of two groups: Group I is made of 236 instances derived from different subnetworks of the intermodal public transportation network around Bruxelles [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF]; Group II includes 64 instances inspired by the railway system in the Randstad, in the Netherlands [START_REF] Dollevoet | Delay management including capacities of stations[END_REF]. Notice that we needed to adapt the latter set of instances as follows. Data from [START_REF] Dollevoet | Delay management including capacities of stations[END_REF] are related to an event-activity network of type Dep-Arr, so we shrink it to an Arr-Arr graph (see [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF]. Moreover, the data describe the path and times of each single vehicle and vehicles on different lines run according to different time periods: we consider instead a fixed time period T coherent with the first group of instances, and we consider the lines along which vehicles move. The number of passengers provided is related to their origin and destination stations and the time in which people start their journey, without information on which path is used in between, while we need passengers associated to each route, thus we aggregate the number of passengers from each origin to each destination within the provided time window and re-distribute them on the routes we take into account. In particular, we create instances in which passengers use routes with at most one, two or three connection arcs (with no restriction in the total number of arcs) or at most a total of five arcs (with no restriction in the number of connections). Different instances are obtained considering a different network or a different initially delayed node. As from Table 2, instances of Group I have been divided into two size classes S and L, according to the number of variables and constraints needed by model MILP1: S includes the instances where the number of variables times the number of constraints is not larger than 100000, L the remaining instances. Group II only has instances in class L. The initial delay D is equal to 4, 8, 18 or 28 minutes and the period T between two vehicles of the same line has been set to 30 minutes. In the following, we report the results obtained on a 2 GHz Intel Core i7 processor, 64 bit operating system with 4GB RAM, with Cplex presolve, heuristics and general purpose cuts switched off, and a time limit of 300 seconds per instance.

A first set of experiments performed on Group I instances aims at verifying the impact of the proposed model and cuts improvements on the strength of the formulations for DMP, in terms of bound provided by the linear relaxation. Results are reported in Tables 3 and4. Rows aggregate the instances by size class (S or L) and by initial delay D ("all" refers to the average of all instances within the same size class), and, further, the last row reports average results for all the instances. The columns of Table 3 are divided into two parts, related respectively to models MILP1 and MILP2, and to MILP1' and MILP2'. It presents results on the effect of using tighter constants in the new models. In column %int the number of instances that are solved to integrality by the linear relaxation are reported, while column Gap% shows the average percent integrality gap for the remaining instances, that is the difference between the optimal integer value and the value of the linear relaxation, divided by the optimal integer value and multiplied by 100. We remark that, since MILP1 and MILP2 have the same linear relaxation, they share columns %int and Gap%.

The same holds for MILP1' and MILP2'. As we can easily see, the percentage of instances solved to integrality increases thanks to the stronger constants, and the overall Gap for the remaining instances is reduced. In general, we observe that the improvements to the model coefficients are quite effective: the number of instances that are solved to integrality by the linear relaxation of MILP1' and MILP2' is always remarkably larger (from 54.24% to 73.73% on average), and even the gap for the remaining instances is slightly reduced (from 5.36% to 4.88%). The third column of each part of Table 3 reports the running times of MILP1 and MILP1' models, whose linear relaxations have been always solved within the time limit. In the last two columns of each part, we show results for MILP2 and MILP2', namely the percentage of instance whose linear relaxation is solved within the time limit and the related computing time (average on solved instances). We notice that MILP1 and MILP1' not only solve all the instances within the time limit, but the running time is always fairly smaller than MILP2 and MILP2'. As a matter of fact, the proposed constraint separation procedure shows rather slow, so that, in the following, we will focus on results from MILP1 and MILP1' models.

Table 4 shows the impact of the new cuts, generated using the two different heuristic approaches HP1 and HP2. To this end, from now on, we just consider those instances such that the linear relaxation of MILP1' and MILP2' does not provide the integer optimal solution without applying cuts (the same instances considered for the column Gap% of the second half of Table 3): column #inst2 shows their number. With reference to these second round instances, the column %solved2 reports the percentage of instances solved within the time limit. Columns %int2 and Gap2% show respectively the percentage of instances solved to integrality and the percent integrality gap for the remaining ones. Column #cuts shows the number of cuts separated and time2 is the computing time. Results show that separating inequalities (34) yields further sensible improvements, both in terms of percentage of instances closed at the root node, and integrality gap in the remaining cases.

It can be seen how procedure HP2 has better performance on MILP1' model, as it solves the same number of instances as HP1 in smaller running times for all instance subgroups but I L with initial delay of 8 minutes: in this case, we are solving one extra instance that needs a large number of cuts and more running time than the average to be solved. For this reason, test on the branch-and-cut procedure considers the algorithm configuration MILP1'+HP2.

The branch-and-cut procedure has been tested on instances of both Group I and Group II and results are presented in Table 5, where we consider only the second round instances, that is the ones that are not solved to integrality at the root node without adding any cut: they correspond to instances of Table 4 and to all the instances of Group II with delay 18 and 28 minutes. The 3 rd and 4 th columns of this table show results regarding the pure branch-and-bound procedure, while the next 4 columns show how the new cuts (configuration MILP1'+HP2) impact on the solution process in terms of number of instances solved within the time limit and, with reference to these instances, average number of branch-and-bound nodes, average number of applied cuts and average running times. We remark that the row "average" refers to the instances of Group I. In general, the number of nodes in the branching tree is reduced, due to the addition of cuts, which confirms the effectiveness of the proposed cuts. However, the separation procedures seem to lack in efficiency, so that the running times are definitively larger than the ones needed by MILP1'. In fact, HP2 enumerates all possible pairs of paths p ∈ P , p ∈ Π v p , and their number may be prohibitively large, especially for the larger instances. In particular, as soon as the cut separation procedure is triggered and starts providing violated inequalities, no instance of Group II can be solved, so that just instances with integer linear relaxation for MILP1' (among which all the ones with smaller initial delays) terminate within the time limit.

The running times of the branch-and-cut procedures may be reduced by generating a small subset of cuts, in order to limit the computational effort devoted to separation. To this end, we performed some preliminary experiments in which we selected just a few cuts to be added at the root node. We recall that HP2 derives a cut (34) from a violated inequality (32) and we observe that a stronger cut may be obtained by choosing the paths p in (32) equal to p whenever possible.

Thus, in the attempt of selecting more effective cuts, we modify HP2 by choosing 5 report the results of this cut-and-branch method where, after a preliminary tuning on the number of cuts to be added to the root node, a maximum of 11 cuts are added at the root node. We report, for each group, the average number of cuts, the average number of processed nodes, the average running time and, in the last column, the percentage of instances where the cut-and-branch procedure is faster than MILP1'. In general, the cut-and-branch approach solves a larger percentage of instances within the time limit. Moreover, as expected, the average number of nodes is smaller than MILP1' for every group but for instances of type I-L with 28 minutes of initial delay: this is due to a single instance showing an unusual large number of nodes (731 for MILP1' and 1887 for cut-and-branch, which also affect the average computing time in this group). In fact, about 58% of the instances of Group I are solved faster, in particular the small ones, as well as instances with small or close to the period T initial delay. Within the same group, the cut-and-branch procedure seems to be faster when the maximum number of connection arcs per passenger path is small (typically 1): in this case, added cuts involve less variables and seem to be more effective.

Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

Concerning instances of Group II, adding a limited number of cuts yields more instances solved within the time limit than MILP1'+HP2, at the cost of larger running times, mostly spent by the separation procedure, so that MILP1' alone is always faster on Group II instances.

The overall results for instances of Groups I and II show that, in order to combine the advantages of a tightened formulation together with reduced running times, a way to select, among all the violated cuts, the most effective ones, as well as separation procedures able to more efficiently deal with large number of paths should be investigated.

Conclusions

In this paper we have presented strengthened formulations for the Delay Management Problem, starting from two formulations MILP1 and MILP2 taken from literature. A first improvement has been obtained by tightening the bounds for the model coefficients, leading to a new formulation MILP1'. By projection, we have obtained a further new formulation MILP2': thanks to the selection of restricted sets of paths, the number of constraints is sensibly reduced with respect to MILP2, even if it remains exponential in the worst case.

We have also proposed new classes of valid inequalities for the new formulations, obtained by studying one single MILP2' constraint and deriving facets for the related polytope: we have shown that the new inequalities generalize the ones presented in literature and provide stronger cuts.

In order to analyze the potential of these theoretical results in computational terms, we have implemented different branch-and-cut algorithms for DMP. To this end, preliminary separation procedures have been provided for both MILP2' constraints and the new family of cuts. Results on a benchmark of instances derived from real transportation networks show that the use of the new formulations and cuts can lead to more efficient solution methods.

In fact, the performance of the models with improved constants are sensibly better, in terms of instances solved at the root node, integrality gap, number of branch-and-bound nodes and efficiency.

In particular MILP1' is sensibly faster than MILP1 and MILP2. Concerning MILP2', a better implementation may be obtained by investigating more efficient constraint separation procedures that do not enumerate all passenger paths. Also, computational results show that the cuts are quite effective in tightening the linear relaxation and reducing the number of branch-and-cut nodes. Even in this case, the cut separation procedures should be fastened, by investigating cut selection procedures that avoid the enumeration of all possible pairs of paths. In this direction, some preliminary experiments, where a restricted number of selected cuts is applied, often provide sensible improvements in the overall efficiency and confirm the potential of the proposed family of cuts.

  which makes (18) redundant, being the second term of the right hand side at least 1. Since conv(Y f DM P ) is defined by one of MILP2 model's constraints, inequalities (18) are valid for MILP2. Furthermore, they are valid for MILP1 model, as they do not involve the variables d k and MILP2 is the projection of MILP1 on the space of the variables z p and u p . As a consequence, (18) can be used as cuts to strenghten both MILP1 and MILP2 formulations. De Giovanni et al. (2014) solve the linear relaxation of MILP1 by dinamically adding inequalities (18) separated by a heuristic procedure, showing sensible improvements in terms of better dual bounds and number of instances closed at the root node.

  variables or less. We demonstrate this case by the following example. We consider a path p starting at the delayed node 1 and containing two connection arcs (a, b) and (c, d), such that α(p, a, b) = 1 and α(p, c, d) = 2. Since p ∈ P p 1 ∩ P p 2 , by considering p = p and the combination with p 1 = p 2 = p, one of the constraints (29) of MILP2' is given by

  ) -π 0 , η(s) = max j∈J π j y j j∈J a j y j ≤ b + s, y ∈ {0, 1} |J| .Proposition 4. Assuming that D > (i,j)∈Ap σ ij , inequality (34) defines a facet of conv(Y f DM P ).Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

  ) is still valid for Y f DM P and reduces to (34) in the specialcase |W | = |A C p |. In the case |W | < |A C p |,(32) are not facet defining, since they can be obtained by combining the corresponding inequality (34) with suitably chosen valid inequalities of type z w ≤ 1.

  13

p

  α(p,i,j) = p if p ∈ A C p arg max p∈P :(i,j)∈A C p z * p otherwise where z * p is the current fractional solution. The last five columns of Table

  p 1 , p 2 . . . p |A C p | . De Giovanni et al. (2014) focus on the case when p i = p j and p i = p for each

	i, j ∈ {1 . . . |A C p |}, i.e., |A C p | + 1 distinct paths are involved. In this case, the associated Mixed 0-1
	Knapsack set is

Table 1

 1 Branch-and-cut algorithms for DMP.

	Algorithm:	MILP1 MILP2 MILP1' MILP2' MILP1'+HP1 MILP1'+HP2
	Model	MILP1 MILP2 MILP1' MILP2'	MILP1'	MILP1'
	Constraints / Procedure	-	(14) / A	-	(29) / A	-	-
	Cuts / Procedure	-	-	-	-	(34) / HP1	(34) / HP2

-: not used

  Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

Table 2

 2 Instance summary: size groups.

	Group		min #var. max #var. min #constr. max #constr. #instances T period
	I	S	136	336	99	246	184	30
		L	382	12833	309	19023	52	30
	II	L	610	4522	508	7762	64	30

Table 3

 3 Linear relaxation without cuts.

	Group I D		MILP1	MILP2	MILP1'	MILP2'
		%int Gap % time %solved time %int Gap%	time	%solved time
	I S	all 50.00	5.14	0.05	100.00 0.55 72.28 4.64	0.05	100.00 0.41
		4 91.30	6.31	0.05	100.00 0.10 95.65 1.71	0.05	100.00 0.06
		8 41.30	6.23	0.05	100.00 0.44 67.39 5.76	0.05	100.00 0.28
		18 32.61	5.84	0.04	100.00 0.91 56.52 4.74	0.05	100.00 0.67
		28 34.78	3.30	0.05	100.00 0.73 69.57 3.69	0.04	100.00 0.65
	I L	all 69.23	6.62	0.23	92.31 3.14 78.85 6.02	0.30	94.23	8.05
		4 100.00 0.00	0.22	92.31 1.32 100.00 0.00	0.27	100.00 21.73
		8 69.23	3.33	0.22	92.31 2.26 76.92 2.01	0.31	92.31	1.36
		18 53.85	7.32	0.23	92.31 4.16 69.23 5.66	0.31	92.31	4.05
		28 53.85	8.12	0.24	92.31 4.82 69.23 9.39	0.31	92.31	3.93
	average	54.24	5.36	0.09	98.31 1.08 73.73 4.88	0.10	98.73	2.02

Table 4

 4 Linear relaxation with cuts (MILP1'+HP1 and +HP2).

	Group I D #inst2		MILP1'+HP1				MILP1'+HP2
				%solved2 %int2 Gap2% #cuts time2 %solved2 %int2 Gap2% #cuts time2
	I S	all	51	100.00	43.14	6.51	69.76 0.49	100.00	74.51	3.05	80.18	0.27
		4	2	100.00	0.00	1.71	0.00	0.12	100.00	0.00	0.99	1.00	0.12
		8	15	100.00	46.67	7.23	43.00 0.35	100.00	80.00	5.26	47.07	0.16
		18	20	100.00	35.00	5.37	88.65 0.56	100.00	75.00	2.00	104.10 0.34
		28	14	100.00	57.14 10.89 81.43 0.58	100.00	78.57	3.97	92.79	0.32
	I L	all	11	72.73	50.00	6.24		1.34	81.82	77.78	1.92	697.33 33.68
		4	0	-	-	-	-	-	-	-	-	-	-
		8	3	66.67	100.00 0.00	29.00 0.62	100.00 100.00 0.00 1790.33 99.24
		18	4	75.00	33.33	1.19	74.00 0.94	75.00	66.67	0.90	73.67	0.48
		28	4	75.00	33.33 11..29 133.67 2.24	75.00	66.67	2.94	228.00 1.31
	average			95.16	44.07	6.47	71.85 0.60	96.77	75.00	2.90	172.75 5.28

-: no instance is taken into account in this group.

  Dal Sasso, De Giovanni, and Labbé: Strengthened Formulations and Cuts for Delay Management 20 Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

Table 5

 5 Branch-and-bound: results for MILP1' and MILP1'+HP2.

	Group D	MILP1'		MILP1'+HP2		Cut-and-branch MILP1'+HP2
			#nodes time2 %solved2 #nodes #cuts time2 %solved2 #nodes #cuts time2 %fast
	I S	all	8.29	0.14	100.00	2.10	80.49	0.33	100.00	3.98	10.01 0.11 62.75
		4	3.00	0.12	100.00	3.00	1.00	0.27	100.00	1.00	2.00	0.09 100.00
		8	4.40	0.10	100.00	1.87	46.87	0.23	100.00	3.40	11.00 0.11 26.66
		18	9.75	0.14	100.00	1.75	104.15 0.37	100.00	3.95	10.10 0.11 70.00
		28 11.14	0.16	100.00	2.71	94.07	0.39	100.00	5.07	10.21 0.13 85.71
	I L	all 123.27 2.59	81.82	2.00	703.56 33.97	100.00	218.00 2.54	2.54 36.36
		4	-	-	-	-	-	-	-	-	-	-	-
		8	4.00	0.39	100.00	1.00	1790.33 99.75	100.00	2.33	11.00 0.37	0.00
		18 140.75 3.14	75.00	1.67	73.67	0.49	100.00	117.75 10.50 1.88 50.00
		28 195.25 3.70	75.00	3.33	246.67 1.66	100.00	480.00 11.00 4.83 50.00
	average		28.69	0.57	96.77	2.10	86.61	0.40	100.00	41.95 10.21 0.54 58.06
	II L	all	2.31	0.36	34.38	1.00	0.00	2.73	50.00	1.75	3.19 45.95 0.00
		4	-	-	-	-	-	-	-	-	-	-	-
		8	-	-	-	-	-	-	-	-	-	-	-
		18	2.81	0.40	50.00	1.00	0.00	2.75	56.25	1.22	1.22 16.07 0.00
		28	4.44	0.42	18.75	1.00	0.00	2.67	43.75	2.43	5.71 84.38 0.00

-: no instance is taken into account in this group.

Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

† The research of the third author is supported by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office. Dal Sasso, De Giovanni, and Labbé: Strengthened Formulations and Cuts for Delay Management Article submitted to Transportation Science; manuscript no.

Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

where y w are the variables that characterize the set W . If the value of the optimal solution is strictly negative, the related W , together with p and p, identifies a violated cut, otherwise the current solution (z * , u * ) violates no cuts. We remark that p is not included in the objective function (no associated penalty), but it is included in the constraints (it can be used to cover the connection arcs), so that if p = p , the model has an optimal solution with all variables set to 0 but y p = 1. In fact, if z * p > u * p D-(i,j)∈Ap σ ij , a violated inequality (34) exists with W = ∅. Solving the set covering formulation is computationally expensive and we propose the following heuristic procedure HP1, based on the well-known approach by [START_REF] Chvàtal | A greedy heuristic for the set-covering problem[END_REF]:

Procedure HP1: , where n is the number of connection arcs of p included in w and not yet covered by other paths in W . The procedure stops as soon as all connection arcs of p are covered. Since the output W may be not optimal, the procedure may fail in identifying a violated cut, even if it exists, that is, the cut separation procedure HP1 is heuristic.

A further heuristic solution to the set covering problem comes from the observation that, given p ∈ P and p ∈ Π v p , it is easy to separate inequalities (32), since the maximum violation comes from setting p α(p,i,j) = arg max p∈P :(i,j)∈A C p z * p . We can thus consider W as the set of involved path p (different from p ), and check if the corresponding (34) is violated. To perform this step we can apply the same separation procedure used in De [START_REF] De Giovanni | A polyhedral study for delay management in public transportation[END_REF], based on enumerating all possible pairs (p , p) and checking the corresponding inequality.

We denote this procedure as Procedure HP2 and we remark that it determines the inequality that maximizes the violation of (32), whereas the corresponding W is a heuristic solution to the set covering problem and may not maximize the violation of (34). As a consequence, Procedure HP2 is heuristic, as shown by the following example. Let p and p be paths such that z * p = 0.3, u * p = 0, p includes 3 connection arcs (a, b), (c, d) and (e, f ). Suppose we have two paths q and r with d), (e, f )}. Let z * q = 0.85 and let z * r = 0.8. In this case, the maximum violation for (32) yields the inequality z p + z q + z q + z r ≤ 3 + 0, and then the inequality (34) z p + z q + z r ≤ 2 + 0 (corresponding to W = {q, r}) would be checked and no violated inequality would be found. Actually, the optimal solution to the set covering problem is W = {r}, yielding the inequality z p + z r ≤ 1 + 0, which is violated.

The investigation of more efficient separation procedures for MILP2' constraints and for the new improved cuts as well as their integration into branch-and-cut methods for DMP are the object of ongoing research.

Further research could also consider to extend these models to deal with more that one initial delay: model MILP1' is easily extended, as it suffices to add a suitable constraint for each distinct initial delay. An extension of model MILP2' has to take into account the propagation of the different delays, so that it requires different families of extra constraints and related cuts should be studied.

Appendix A: Proof of Proposition 2

Proposition 2. The following formulation MILP2' is valid for DMP and the value of the linear relaxation is the same as MILP1':

Proof. The result is obtained by showing that the feasible domain of MILP2' is the projection of the feasible domain of MILP1' onto the space of u and z variables. This can be done by applying Fourier-Motzkin elimination to d variables with steps similar to the ones used by [START_REF] Heilporn | Optimization models for the single delay management problem in public transportation[END_REF] to obtain MILP2 from MILP1. First, we rewrite the constraints ( 21) and ( 22): from ( 22) we have

and, by defining for each arc (i, j)

we can write

Now, let p be one of the paths starting at the delayed node 1. For ease of notation, we denote with 1 the first node in path p, with 2 the second and so on, until the last node v p . By constraints (20), ( 21) and ( 22)

Article submitted to Transportation Science; manuscript no. (Please, provide the mansucript number!)

Similarly, for any path p whose first node is k = 1, starting from d k ≥ 0, we have

Hence, we have

Notice that, since l ij is nonnegative, the last inequalities are redundant and can be omitted.

We remark also that paths p ∈ Π with at least one connection not included in any passenger path are not considered, since the delay does not spread on them, according to constraints ( 22).

We recall that, by ( 23),

As it holds for every p ∈ P , it also holds for a p ∈ P such that v p = v p :

Noticing that

and recalling that just variables z and u appear in the objective function, we can apply Fourier-Motzkin elimination and replace constraints (20), ( 21), ( 22) and ( 23) by

Let us now consider a passenger path p ∈ P such that pv p , i.e. the path from the delayed node 1 to v p with no connection arcs, exists. Let us define, for ease of notation, p = pv p . In this case

and one of the constraints (40) is

Consider also another path p ∈ Π from 1 to v p = v p such that

The resulting constraints of type ( 40) is

25 which is redundant. In fact, as

Summarizing, (i,j)∈Ap l ij ≥ Θ v p , and hence D -

meaning that ( 41) is redundant.

The discussion above shows that MILP1' is equivalent to

In case (i, j) ∈ A C , the term l ij is not linear and can be written as follows:

where zij = max p∈P :(i,j)∈A C p z p . Then MILP1' is equivalent to

This model is still not linear, because of the term zij . In order to get a linear formulation, we write the same inequality involving z for each combination of paths in P such that the t-th path in the combination includes the t-th connection arc of p. The resulting constraints are

yielding formulation MILP2'.

Appendix B: Proofs of technical lemmas

Lemma 1. For any path p ∈ Π such that a p = 1, M vp ≥ D -(i,j)∈Ap σ ij .

Proof. Let k be one of the nodes of p and denote by p(k) the sub-path of p from 1 to k. By defining Mk = max{0, D -

we have, for each node k included in the path p, Mk = D -(i,j)∈A p(k) σij . By the definition of M k , we have Mk ≤ M k and, by definition of

and hence

Lemma 2. Assuming that D > (i,j)∈Ap σ ij , W ∪ {p } is a minimal cover of Y f DM P and the inequality

Proof. First, we prove that the set is a cover. In fact, as by hypothesis D > (i,j)∈Ap σ ij ,

In order to prove that the cover is minimal, we show that W ∪ {p } \ w is not a cover for each w ∈ W ∪ {p }.

We have two cases: w = p and w ∈ W .

If w = p , recalling that A C p (w) denotes the set {(h, l) ∈ A C p |p α(p,h,l) = w},

If w ∈ W , we have

σ hl .

We have to show that

Now, let (h 1 , j) be the first connection arc of A C p ( w) found on path p moving from node 1. Denoting with p(h 1 ) the sub-path of p from 1 to h 1 , we have

Since W = ∅, the set A C p \ A C p (p ) = ∅, and let (k, j) be the first connection arc in this set on the way from 1 to v p . Recalling that p(k) is the sub-path of p limited to node k, for M k , by Lemma 1, we have

Summarizing RHS -γ w ≥ 0, for any w ∈ W ∪ {p }.