Timetable merging for the Periodic Event Scheduling Problem

Lindner, Niels; Liebchen, Christian Article Timetable merging for the Periodic Event Scheduling Problem EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Lindner, Niels; Liebchen, Christian (2022) : Timetable merging for the Periodic Event Scheduling Problem, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 11, Iss. 1, pp. 1-9, https://doi.org/10.1016/j.ejtl.2022.100081 This Version is available at: https://hdl.handle.net/10419/325157 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ EURO Journal on Transportation and Logistics 11 (2022) 100081 Available online 20 April 2022 2192-4376/© 2022 The Author(s). Published by Elsevier B.V. on behalf of Association of European Operational Research Societies (EURO). This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.elsevier.com/locate/ejtl Timetable merging for the Periodic Event Scheduling Problem Niels Lindnera,∗, Christian Liebchenb aZuse Institute Berlin, Takustr. 7, 14195 Berlin, Germany bTechnical University of Applied Sciences Wildau, Hochschulring 1, 15745 Wildau, Germany ARTICLE INFO Keywords: Periodic Event Scheduling Problem Periodic timetabling Railway timetabling PESPlib Benchmark solutions Mixed integer programming ABSTRACT We propose a new mixed integer programming based heuristic for computing new benchmark primal solutions for instances of the PESPlib. The PESPlib is a collection of instances for the Periodic Event Scheduling Problem (PESP), comprising periodic timetabling problems inspired by real-world railway timetabling settings, and attracting several international research teams during the last years. We describe two strategies to merge a set of good periodic timetables. These make use of the instance structure and minimum weight cycle bases, finally leading to restricted mixed integer programming formulations with tighter variable bounds. Implementing this timetable merging approach in a concurrent solver, we improve the objective values of the best known solutions for the smallest and largest PESPlib instances by 1.7 and 4.3 percent, respectively. 1. Introduction In periodic event scheduling, we consider a set of events, where each event repeats periodically with the same period time 𝑇. A periodic timetable is then an assignment of time values within the interval [0, 𝑇 ) to all the events, subject to the condition that the (periodic) time differences between certain pairs of events meet a feasibility interval. This type of problem arises in particular in the context of computing cyclic timetables in public transport (aka fixed-interval timetables), where one may think of a periodic event as either arrival or departure of a directed traffic line at some station, and in coordinating the traffic lights at several road intersections (Hassin,1996;Wünsch,2008). In the context of timetabling in public transport, the Periodic Event Scheduling Problem (PESP) introduced by Serafini and Ukovich (1989) has been the core model for computing the first optimized railway timetable in practice (Liebchen,2008), and also for the timetabling part of the company-wide success story of Operations Research at Dutch Railways (Kroon et al.,2009). In order to attract researchers in developing more powerful solution methods for periodic event scheduling, Goerigk (2012) composed a collection of instances (PESPlib), most of which are motivated from railway timetabling. This library attracted several international research teams, however, none of the PESPlib instances has been solved to proven optimality. Finding provably optimal timetables is out of reach with the currently available techniques, even for moderatelysized instances. However, there is plenty of heuristics that are able to compute good quality timetables. The current incumbent solutions in the PESPlib stem from the concurrent PESP solver by Borndörfer ∗Corresponding author. E-mail addresses: [email protected] (N. Lindner), [email protected] (C. Liebchen). et al. (2020), in which a variety of these algorithms and heuristics are used as subroutines. In particular, the most common local search procedures are implemented, such that the best solutions will constitute local minima for all of these heuristics when the solver terminates. We present a strategy to escape local minima to further improve timetables that cannot be improved anymore by the aforementioned local search procedures. We follow the spirit of Cook and Seymour (2003): Their tour merging heuristic constructs the union of ten heuristically determined tours for the Traveling Salesman Problem (TSP), applies a branch-decomposition-based dynamic program to this restricted instance, and is able to produce even better TSP tours. Our philosophy is to merge a small set of good quality timetables. But we do not copy two half timetables and simply compile their parts. Rather, we do so by restricting the bounds of the periodic tension or cycle offset variables in the cycle-based mixed integer programming formulation of PESP. We exploit the similarities of the input timetables, and take into account the special characteristics of a periodic timetabling problem in public transport. This leads to two timetable merging heuristics, one for each variable type. From the perspective of mixed integer programs, our method can be understood as a generalization of the well-known crossover heuristic (Rothberg,2007). We tighten the variable bounds such that all the solutions in the selected small set of good quality timetables remain feasible. However, our approach is more general, as we consider more than two initial solutions and allow bound restrictions of arbitrary variables. This is because merely fixing integer variables with the same solution value rarely leads to better solutions. Moreover, the crossover https://doi.org/10.1016/j.ejtl.2022.100081 Received 3 December 2021; Received in revised form 12 April 2022; Accepted 14 April 2022 EURO Journal on Transportation and Logistics 11 (2022) 100081 2 N. Lindner and C. Liebchen heuristic for general mixed integer programs is supposed to run very fast for the purpose of solution polishing, but we spend a considerable amount of time on solving the restricted problem, as solution quality is our main goal, and PESP is notoriously hard. It turns out that applying this tailored timetable merging heuristic can make better solutions accessible. Invoking the concurrent PESP solver that computed the current PESPlib incumbents, we find several timetables improving the best known primal bounds, and end up with a 1.7 and 4.3 percent lower objective value for the smallest and largest PESPlib instance, respectively. We define PESP in Section 2, where we also review some problem features, and introduce the cycle-based mixed integer programming formulation. Our timetable merging heuristic is discussed in Section 3. Section 4presents computational results, on the one hand on the structure of the restricted scenarios from the merging heuristic, and on finding better periodic timetables for the PESPlib instances on the other. We conclude this paper in Section 5. 2. The periodic event scheduling problem In this section we first define the combinatorial optimization problem for which we are about to propose our merging heuristics, provide a short review of other solution techniques that have been applied earlier, and on some of which we are going to build on, and formulate the mixed-integer linear programming (MIP) formulation that constitutes the basis of our computational study. 2.1. Problem definition The Periodic Event Scheduling Problem (PESP) dates back to Serafini and Ukovich (1989), and shows certain similarities to models that were already considered by Rüger (1986). Formally, the input is a tuple (𝐺, 𝑇 , 𝓁, 𝑢, 𝑤), where •𝐺= (𝑉 , 𝐴)is a directed graph, called event–activity network, whose vertices are called events and whose arcs are called activities, •𝑇∈Nis a period time, •𝓁∈R𝐴 ≥0is a vector of lower bounds such that 0≤𝓁< 𝑇 , •𝑢∈R𝐴 ≥0is a vector of upper bounds,0≤𝑢−𝓁< 𝑇 , and •𝑤∈R𝐴 ≥0is a vector of weights. A vector 𝜋∈ [0, 𝑇 )𝑉is called a periodic timetable if there exists a periodic tension 𝑥∈R𝐴such that 𝓁≤𝑥≤𝑢and ∀𝑎= (𝑖, 𝑗) ∈ 𝐴∶𝜋𝑗−𝜋𝑖≡𝑥𝑎mod 𝑇 . Intuitively, a periodic timetable 𝜋assigns times modulo 𝑇to each event in 𝐺, and fixes the duration of any activity 𝑎= (𝑖, 𝑗) ∈ 𝐴to 𝜋𝑗−𝜋𝑖 modulo 𝑇. The actual duration of 𝑎is then chosen to lie within the interval [𝓁𝑎, 𝑢𝑎]. Starting from a periodic timetable 𝜋, a corresponding periodic tension 𝑥can be computed by 𝑥𝑎∶= [𝜋𝑗−𝜋𝑖−𝓁𝑎]𝑇+𝓁𝑎for all 𝑎= (𝑖, 𝑗) ∈ 𝐴, where [⋅]𝑇denotes the modulo 𝑇operator with values in [0, 𝑇 ). Conversely, given a periodic tension 𝑥, a corresponding periodic timetable 𝜋can be reconstructed by a graph traversal. We further define the periodic slack as 𝑦∶= 𝑥−𝓁∈R𝐴 ≥0. In a public transport context, an event 𝑖is usually modeling either the arrival or the departure of a directed traffic line at some station, e.g., the departure of the trains from Berlin to Munich in the city of Erfurt. An arc 𝑎= (𝑖, 𝑗)models the time duration from event 𝑖to event 𝑗. If 𝑖and 𝑗are two subsequent departure and arrival events of the same directed line, then 𝑎= (𝑖, 𝑗)models the trip duration from the station of event 𝑖to the station of event 𝑗. In turn, if 𝑖and 𝑗are the arrival and departure events of the same directed line within the same station, then 𝑎= (𝑖, 𝑗)models the dwell duration within this station. To illustrate many other commercial and operational types of constraints, we refer to Liebchen and Möhring (2004). If in an hourly service (i.e., 𝑇= 60 minutes), for a dwell arc 𝑎= (𝑖, 𝑗)we require that 𝓁𝑎= 3 and 𝑢𝑎= 7, then of course 𝜋𝑖= 29 and 𝜋𝑗= 33 constitute a feasible timetable, because 33 − 29 = 4 ∈ [3,7]. However, notice that 𝜋𝑖= 58 and 𝜋𝑗= 3 constitute a feasible timetable, too, because 𝑥𝑎= [3 − 58 − 3]60 + 3 = 2 + 3 = 5 ∈ [3,7]. Definition 2.1. Given (𝐺, 𝑇 , 𝓁, 𝑢, 𝑤)as above, the Periodic Event Scheduling Problem (PESP) is to find a periodic slack 𝑦such that the weighted periodic slack ∑𝑎∈𝐴𝑤𝑎𝑦𝑎is minimum or to decide that no periodic timetable exists. Alternatively, one may seek to minimize the weighted periodic tension ∑𝑎∈𝐴𝑤𝑎𝑥𝑎, which differs from the corresponding weighted periodic slack by the constant ∑𝑎∈𝐴𝑤𝑎𝓁𝑎. Speaking in application terms, where weights may reflect the number of passengers using an activity, this amounts to minimizing the total passenger travel time, given that the routes that the passengers are taking – and thus the activities on which they are showing up – are known in advance. Notice that recently there have been made advances in relaxing this assumption (Schiewe and Schöbel,2020). But the activities’ weights may also model other practical aspects, as for instance in the context of minimizing the amount of rolling stock that is required to operate a periodic timetable (Liebchen and Möhring,2004). Notice that in the literature, sometimes PESP is regarded as the pure feasibility decision problem, thus with no weights 𝑤𝑎defined on the arcs. 2.2. Complexity and solution approaches For an arbitrary PESP instance, it is not clear at all that a periodic timetable resp. tension resp. slack exists. In fact, the feasibility problem is NP-complete even if 𝑇≥3is not regarded as part of the input, because PESP generalizes Vertex Coloring (Odijk,1994). Moreover, the feasibility can be checked in linear time if undirecting 𝐺results in a forest, but is already NP-hard for graphs of treewidth or branchwidth ≥2(Lindner and Reisch,2020). In particular, in this paper we deal with an NP-hard optimization problem. A lot of powerful tools stimulated by different viewpoints within discrete optimization have been developed to solve PESP instances: These tools comprise, e.g., mixed-integer programming (Liebchen,2006), simplex-style algorithms (Nachtigall and Opitz, 2008;Goerigk and Schöbel,2013), satisfiability methods (Großmann et al.,2012), machine learning (Matos et al.,2020), matching heuristics (Pätzold and Schöbel,2016), and graph partitioning approaches (Lindner and Liebchen,2019). A large part of these methods has been integrated into a single PESP solver based on concurrency, ConcurrentPESP (Borndörfer et al.,2020). However, up to today even medium-sized PESP instances withstand all attempts to compute optimal solutions: Since 2012, none of the 20 instances of the PESP benchmark library PESPlib (Goerigk,2012) could be solved to optimality, the current best primal–dual gap being ≈35%, where the smallest instance R1L1 even has no more than 6,386 activities. For all these instances, ConcurrentPESP is the record holder concerning both primal and dual bounds. 2.3. Mixed integer programming formulation Our merging strategy relies on the well-known cycle-based mixed integer programming formulation for PESP: Minimize ∑ 𝑎∈𝐴 𝑤𝑎(𝑥𝑎−𝓁𝑎) subject to 𝛤 𝑥 =𝑇 𝑧, 𝓁≤𝑥≤𝑢, 𝑧∈Z𝐵. (1) EURO Journal on Transportation and Logistics 11 (2022) 100081 3 N. Lindner and C. Liebchen Here, 𝐵denotes an integral cycle basis of 𝐺with cycle matrix 𝛤, see, e.g., Liebchen (2006) for details. The variables are the periodic tension 𝑥as defined above and the cycle offset 𝑧∈Z𝐵. Starting from a periodic tension 𝑥, the corresponding cycle offset is of course recovered by computing 𝑧=𝛤 𝑥∕𝑇. Conversely, given 𝑧, an optimal periodic tension 𝑥for 𝑧can be found by a minimum cost network flow computation (Nachtigall and Opitz,2008). In the case of integer input vectors 𝓁and 𝑢, the existence of a feasible solution implies the existence of an integral solution vector 𝑥, as 𝛤is the cycle matrix of an integral cycle basis. In addition, recall from Odijk (1994) that for each oriented cycle 𝛾∈ {−1,0,1}𝐴the following cycle inequalities 𝑧ODIJK 𝛾∶= ⌈𝛾𝑡 +𝓁−𝛾𝑡 −𝑢 𝑇⌉≤𝑧𝛾≤⌊𝛾𝑡 +𝑢−𝛾𝑡 −𝓁 𝑇⌋=∶ 𝑧ODIJK 𝛾(2) are valid and turned out to be useful for solving (1) in many computational studies. Here, an oriented cycle is the incidence vector of a cycle in 𝐺, where arcs traversed in forward direction have a +1 entry, and those in backward direction have a −1 entry. For such an incidence vector, we let 𝛾+∶= max(𝛾, 0) resp. 𝛾−∶= max(−𝛾, 0) denote the positive resp. negative part of 𝛾, so that 𝛾=𝛾+−𝛾−. 3. Merging timetables The incumbent solutions in the PESPlib are hard to improve further: They are locally optimal for the local heuristics implemented in the ConcurrentPESP solver, so that only the global component realized by branch-and-cut is able to yield better solutions. This is unsatisfactory, given the enormous size of the branch-and-bound trees being a result of the weak trivial linear programming relaxations (Liebchen, 2006). Bound restriction approaches are not new to periodic event scheduling. For example, the Arc Selection heuristic suggested by Roth (2019) tries to improve a periodic tension 𝑥by decreasing the upper bounds 𝑢to 𝑥and iteratively re-optimizing over subsets of activities. The optimization is carried out by calling a MaxSAT solver after transforming the PESP instance to an instance of the weighted partial maximum satisfiability problem. On large instances, this approach is reported to be superior to pure branch-and-cut. Although inspired by Roth (2019), our idea is somewhat different: For a PESP instance, we consider a set 𝑆of solutions with ‘‘good’’ objective values. We now want to identify certain variables of (1) whose values are the same or at least similar for all the solutions in 𝑆. Then, by sharpening bounds for these variables, we aim at deriving a modified mixed-integer program that •is easier to solve because of the tighter bounds, •still contains all the already good solutions of 𝑆as feasible solutions, and •has a feasible set which is a subset of the feasible set of the initial problem. Ideally, solving this sharpened problem, we could find even better solutions for the initial problem. As there are two types of variables in (1), there are two directions to pursue: Restricting the periodic tension 𝑥or restricting the cycle offset 𝑧. We will discuss the details of our strategy for these two cases separately. 3.1. Restricting periodic tensions Given a PESP instance 𝐼= (𝐺, 𝑇 , 𝓁, 𝑢, 𝑤), we first analyze the span 𝑢𝑎−𝓁𝑎of an activity 𝑎. For example, in the smallest PESPlib instance R1L1, having 𝑇= 60, it turns out that each activity either has span at most 17 or exactly 59. A similar structure is present in the other railway instances of the PESPlib. Motivated by the typical modeling of event–activity networks for periodic timetabling in public transport, we refer to these free activities in the group with the larger span as transfer activities, which is also supported by an investigation of the connectivity structure of the PESPlib instances due to Goerigk and Liebchen (2017). With integer bounds, as is the case for the PESPlib instances, a span of 59 on an activity 𝑎in conjunction with a period time 𝑇= 60 means that any integer periodic slack 𝑦𝑎∈ [0,60) is feasible. However, restricting the span of these activities in a too rigorous way might turn the PESP instance infeasible. Of course, increasing the lower bounds and decreasing the upper bounds of the activities is our main goal. However, in order to keep some flexibility, we restrict the upper bounds in a different manner than the lower bounds: Non-transfer activities as well as transfer activities with low impact on the objective function always remain at their original upper bounds, because there would be no strong empirical evidence which of these arcs will have to face large slack in the best solutions for an instance. To this end, we choose a parameter 𝛼∈ [0,1] and denote by 𝐴𝛼the largest set of transfer activities 𝑎∈𝐴that, when sorted in ascending order w.r.t. weight 𝑤𝑎, sum up to at most 𝛼∑𝑎∈𝐴𝑤𝑎. Now let 𝑆be a set of periodic tensions for 𝐼. We define the following bounds for each activity 𝑎∈𝐴: 𝓁ORIG 𝑎∶= 𝓁𝑎, 𝑢ORIG 𝑎∶= 𝑢𝑎, 𝓁FIX 𝑎∶= min{𝑥𝑎∣𝑥∈𝑆}, 𝑢FIX 𝑎∶= {max{𝑥𝑎∣𝑥∈𝑆}if 𝑎∈𝐴𝛼, 𝑢𝑎otherwise, 𝓁MED 𝑎∶= 1 2𝓁ORIG 𝑎+1 2𝓁FIX 𝑎, 𝑢MED 𝑎∶= 1 4𝑢ORIG 𝑎+3 4𝑢FIX 𝑎. Observe that for all 𝑥∈𝑆we find 𝓁=𝓁ORIG ≤𝓁MED ≤𝓁FIX ≤𝑥≤𝑢FIX ≤𝑢MED ≤𝑢ORIG =𝑢. We can now combine different lower and upper bound strategies freely by considering the PESP instances 𝐼𝑆(𝐿, 𝑈) ∶= (𝐺, 𝑇 , 𝓁𝐿, 𝑢𝑈, 𝑤) for 𝐿, 𝑈 ∈ {ORIG,MED,FIX}. Clearly, 𝐼𝑆(ORIG,ORIG)is the original PESP instance 𝐼. Observe that any solution in 𝑆is feasible for all scenarios 𝐼𝑆(𝐿, 𝑈). Moreover, restricting the changes to either lower or upper bound, a periodic tension for FIX is a periodic tension for MED, and any periodic tension for MED is a periodic tension for ORIG. However, formally, the weighted slacks might differ, as the lower bounds are potentially altered. 3.2. Restricting cycle offsets Consider a PESP instance 𝐼= (𝐺, 𝑇 , 𝓁, 𝑢, 𝑤)with a set 𝑆of periodic tensions. For every integral cycle basis 𝐵with cycle matrix 𝛤, we can compute the corresponding cycle offset vector 𝑧∈Z𝐵for each tension 𝑥∈𝑆. Denote by 𝑆𝐵∶= {𝛤 𝑥∕𝑇∣𝑥∈𝑆}the set of these cycle offsets. The idea is now to find bounds 𝑧, 𝑧 ∈Z𝐵and to solve the restricted mixed-integer program Minimize ∑ 𝑎∈𝐴 𝑤𝑎𝑦𝑎 subject to 𝛤 𝑥 =𝑇 𝑧, 𝓁≤𝑥≤𝑢, 𝑧≤𝑧≤𝑧, 𝑧∈Z𝐵. (3) We do this in two ways by defining for each oriented cycle 𝛾∈𝐵 𝑧ALL 𝛾∶= min{𝑧𝛾∣𝑧∈𝑆𝐵}, 𝑧ALL 𝛾∶= max{𝑧𝛾∣𝑧∈𝑆𝐵}, 𝑧PARTIAL 𝛾∶= {𝑧ALL 𝛾if 𝑧ALL 𝛾=𝑧ALL 𝛾, 𝑧ODIJK 𝛾otherwise, cf. (2), 𝑧PARTIAL 𝛾∶= {𝑧ALL 𝛾if 𝑧ALL 𝛾=𝑧ALL 𝛾, 𝑧ODIJK 𝛾otherwise, cf. (2). EURO Journal on Transportation and Logistics 11 (2022) 100081 4 N. Lindner and C. Liebchen In the PARTIAL version, we fix all cycle offset variables if they are the same in all solutions in 𝑆, and we do not impose any restrictions that would depend on 𝑆on the other entries of 𝑧. Hence, for any set of tensions 𝑆and any integral cycle basis 𝐵, we obtain two restricted mixed-integer programs of the form (3): MIP𝑆,𝐵(𝐼, ALL)and MIP𝑆,𝐵(𝐼, PARTIAL). Note that these do not necessarily correspond to PESP instances anymore, as we have not changed any part of the input data 𝐼= (𝐺, 𝑇 , 𝓁, 𝑢, 𝑤)– we are merely restricting the program (1). Clearly, the solutions in 𝑆are feasible for both MIP𝑆,𝐵(𝐼, ALL)and MIP𝑆,𝐵(𝐼, PARTIAL). Moreover, any feasible solution to MIP𝑆,𝐵(𝐼, ALL)is feasible for MIP𝑆,𝐵(𝐼, PARTIAL), and any feasible solution to MIP𝑆,𝐵(𝐼, PARTIAL)is feasible for the original unrestricted MIP (1). As the cycle inequalities (2) are valid inequalities, the unrestricted MIP (1) is equivalent to MIP𝑆,𝐵(𝐼, ODIJK)for any cycle basis 𝐵and any 𝑆. Referring to MIP𝑆,𝐵(𝐼, ⋅)as the feasible regions of the respective mathematical programs, we summarize 𝑆 ⊆ MIP𝑆,𝐵(𝐼, ALL)⊆MIP𝑆,𝐵(𝐼, PARTIAL)⊆MIP𝑆,𝐵(𝐼, ODIJK) ≡MIP (1). 3.3. Choosing a cycle basis So far, we have not discussed which cycle basis 𝐵to choose when restricting cycle offsets. For the purpose of solving the mixed integer program (3), it is desirable to fix as many integer variables as possible, and more generally, to minimize the possible number of values of 𝑧. Hence, we seek to minimize ∏ 𝛾∈𝐵 (𝑧ALL 𝛾−𝑧ALL 𝛾+ 1), or equivalently, log ∏ 𝛾∈𝐵 (𝑧ALL 𝛾−𝑧ALL 𝛾+ 1) = ∑ 𝛾∈𝐵 log(𝑧ALL 𝛾−𝑧ALL 𝛾+ 1).(4) We call (4) the log width of the cycle basis 𝐵. In the context of periodic timetabling, the concept of integral cycle bases of small (log) width has proven to be a valuable tool in order to speed up the branch-and-bound process for the cycle-based mixed integer programming formulation of PESP (Liebchen and Peeters,2009). Finding a so-called directed or undirected cycle basis of minimum weight is well-understood (Kavitha et al.,2009). Yet, in order to profit from these insights, we must get around the following two pitfalls: First, since a minimum-weight cycle basis in general is not integral, after having computed a minimum-weight cycle basis we need to check, whether it is even an integral cycle basis, because otherwise the computed cycle basis risks to be useless. Second, we can only compute a minimum-weight cycle basis efficiently, if the weight is given as a function on the arcs rather than on the cycles. Therefore, we will construct an arc-weight function 𝑐∶𝐴→R≥0such that for all oriented cycles 𝛾in 𝐺holds 𝑐(𝛾) ∶= ∑ 𝑎∈𝛾 𝑐(𝑎) ≈ 𝑧ALL 𝛾−𝑧ALL 𝛾. Lemma 3.1. Define 𝑐∶𝐴→R≥0via 𝑐(𝑎) ∶= max{𝑥𝑎∣𝑥∈𝑆} − min{𝑥𝑎∣𝑥∈𝑆} 𝑇. Then for all oriented cycles 𝛾∈𝐵consisting of |𝛾|arcs holds 0≤𝑐(𝛾)−(𝑧ALL 𝛾−𝑧ALL 𝛾)≤𝑐(𝛾)<|𝛾|. Proof. Let 𝛾be an oriented cycle in 𝐵. Then 𝑧ALL 𝛾−𝑧ALL 𝛾= max{𝑧𝛾∣𝑧∈𝑆𝐵} − min{𝑧𝛾∣𝑧∈𝑆𝐵} =max{(𝛤 𝑥)𝛾∣𝑥∈𝑆} − min{(𝛤 𝑥)𝛾∣𝑥∈𝑆} 𝑇 =max{𝛾𝑡𝑥∣𝑥∈𝑆} − min{𝛾𝑡𝑥∣𝑥∈𝑆} 𝑇. Recall that we can decompose 𝛾=𝛾+−𝛾−into its positive and negative part. Then for any 𝑥∈𝑆holds 𝛾𝑡 +𝑥−𝛾𝑡 −𝑥≤𝛾𝑡𝑥≤𝛾𝑡 +𝑥−𝛾𝑡 −𝑥, where 𝑥, 𝑥 are the vectors in R𝐴with 𝑥𝑎∶= max{𝑥𝑎∣𝑥∈𝑆}and 𝑥𝑎∶= min{𝑥𝑎∣𝑥∈𝑆}for all 𝑎∈𝐴. In particular 𝑧ALL 𝛾−𝑧ALL 𝛾≤ (𝛾𝑡 +𝑥−𝛾𝑡 −𝑥)−(𝛾𝑡 +𝑥−𝛾𝑡 −𝑥) 𝑇=(𝛾++𝛾−)𝑡(𝑥−𝑥) 𝑇 =∑ 𝑎∈𝛾 𝑐(𝑎) = 𝑐(𝛾). Since 𝑧ALL 𝛾−𝑧ALL 𝛾≥0, we arrive at 0≤𝑐(𝛾)−(𝑧ALL 𝛾−𝑧ALL 𝛾)≤𝑐(𝛾). Finally note that, as we require 𝑢−𝓁< 𝑇 for our PESP instances, we have 𝑐(𝛾)≤∑ 𝑎∈𝛾 𝑢𝑎−𝓁𝑎 𝑇<|𝛾|.□ □ We can now compute a minimum-weight undirected cycle basis 𝐵∗ w.r.t. 𝑐. Note that such a cycle basis is not necessarily integral. However, our empirical observation is that on non-artificial PESP instances, as the ones that can be found in particular in the PESPlib, minimum undirected cycle bases typically turn out to be integral. Being an F2vector space, the undirected cycles form a matroid (Horton,1987), so that 𝐵∗is also of minimum weight w.r.t. log(𝑐) + 1: Up to breaking ties, the greedy algorithm sorts the cycles in the same order both regarding 𝑐and log(𝑐)+1. This enables us to bound the error that we make when we are using the cycle basis 𝐵∗instead of the one that minimizes the actual log width (4). Corollary 3.2. Let 𝐵′be an undirected cycle basis of 𝐺of minimum log width (4). Then 0≤∑ 𝛾∈𝐵∗ log(𝑧ALL 𝛾−𝑧ALL 𝛾+ 1) − ∑ 𝛾∈𝐵′ log(𝑧ALL 𝛾−𝑧ALL 𝛾+ 1) <|𝐵∗|log(|𝑉|+ 1). Proof. The left inequality is clear as 𝐵′minimizes (4). As ∑𝛾∈𝐵′log(𝑧ALL 𝛾−𝑧ALL 𝛾+1) ≥0, it remains to invoke Lemma 3.1 to obtain ∑ 𝛾∈𝐵∗ log(𝑧ALL 𝛾−𝑧ALL 𝛾+ 1) <∑ 𝛾∈𝐵∗ log(|𝛾|+ 1) ≤|𝐵∗|log(|𝑉|+ 1). Note that we can assume |𝛾|≤|𝑉|as there is always a minimum undirected cycle basis composed of simple cycles, and we only compare the objectives. □ An empirical analysis of the quality of different cycle bases will be given in Section 4.2, the cycle basis 𝐵∗as defined in Lemma 3.1 being superior for our purposes. 4. Results We turn now to the computational results that we obtained by timetable merging. We describe the details of the set-up in Section 4.1. In Section 4.2, we present a structural analysis of the restricted scenarios constructed in Section 3, including an assessment of various cycle bases. Finally, the actual timetables found by our strategy are evaluated in Section 4.3. The full set of solutions is available at GitHub.1 1https://github.com/nielslindner/timetable-merging-solutions EURO Journal on Transportation and Logistics 11 (2022) 100081 5 N. Lindner and C. Liebchen Fig. 1. Workflow for a PESP instance 𝐼. 4.1. General setup We test the merging approach outlined in Section 3on the smallest and largest PESPlib instance, R1L1 and R4L4, respectively. For each of the two instances, we choose a set 𝑆of 5 solutions whose weighted slack is close to the current PESPlib record. Among these solutions are the incumbent solutions found by Lindner and Roth (see Borndörfer et al. 2020), and the best solution found by the iterative procedure by Goerigk and Liebchen (2017). Our workflow is visualized in Fig. 1. We construct the 8 scenarios with the tightened bounds as described in Section 3.1, and another 2 scenarios according to Section 3.2. For each of these 10 scenarios, we conduct 18 = 6 ⋅3runs of ConcurrentPESP with 60 min each (wall time): These arise from 6 different initial solutions (the ones from the set 𝑆, or none) in combination with 3 different parameter settings for fine-tuning the actual solution process (Round 1). For each of these runs we check whether the optimum solution of that run achieved an improvement compared to the initial solution that gave rise to the respective run. If this is the case, then the possibility arises to escape local minima and we thus launch Round 2: We feed the ConcurrentPESP solver on the original instance R1L1 resp. R4L4 with the output of Round 1 as input. Each such run of Round 2 is performed as a single solver run for 4 h (wall time), dedicating most computational capacity to MIP solving. For the scenarios of Round 1, we compute dual bounds with separated additional dedicated runs of ConcurrentPESP of 24 hours each, running exclusively a MIP solver with best bound emphasis, and a separator for violated flip inequalities (Lindner and Liebchen, 2020). For the ALL and PARTIAL scenarios with restricted cycle offsets, we have to turn off the modulo network simplex algorithm, as this too frequently violates the tighter bounds. Instead, we adjust the mixed-integer programming based maximum cut heuristic available in ConcurrentPESP to work with the cycle offset bound constraints. The ConcurrentPESP solver is run on up to 8 threads on an Intel Xeon E3-1270 CPU running at 3.80 GHz with 32 GB RAM. We use IBM CPLEX 12.10 as underlying MIP solver. 4.2. Scenario analysis Tension-restricted scenarios The approach from Section 3.1 yields PESP instances 𝐼𝑆(𝐿, 𝑈)for 𝐿, 𝑈 ∈ {FIX,MED,ORIG}, where 𝐼𝑆(ORIG,ORIG)corresponds to the original instance 𝐼∈ {R1L1,R4L4}, and is thus omitted. Our approach does neither alter the graph 𝐺nor the weights 𝑤, it merely affects the lower bounds 𝓁and the upper bounds 𝑢. We demonstrate hence the effect of our approach by analyzing the span: The span of an activity 𝑎∈𝐴is defined as 𝑢𝑎−𝓁𝑎. We will call an activity fixed if its span is 0, and free if the span is 59. The weighted span of an activity 𝑎is given by 𝑤𝑎(𝑢𝑎−𝓁𝑎). Table 1 presents an analysis of the activity spans in the tensionrestricted scenarios for the PESPlib instances R1L1 and R4L4. As parameter 𝛼, we chose 𝛼= 0.125. The ‘‘difficulty’’ of the restricted instances 𝐼(𝐿, 𝑈), in terms of the average span or the average weighted span, is dictated by the upper bound strategy: 𝑈=FIX has lower average span than 𝑈=MED, which in turn has lower average span than 𝑈=ORIG. Each change in 𝑈to a more restrictive policy reduces the weighted average span by roughly a factor of 1 2. For a given strategy for the upper bounds, there is a clear ranking for the lower bound strategies, which is also FIX →MED →ORIG with increasing difficulty. Moreover, the number of free activities is significantly lower than on the original instance. The number of fixed activities is approximately the same for all scenarios with 𝑈∈ {MED,ORIG}, but jumps up to ≈54%–55% (R1L1) resp. ≈ 43% (R4L4) when 𝑈=FIX. Cycle-offset-restricted scenarios and evaluation of cycle bases Recall from Section 3.3 that, when restricting the cycle offset vectors 𝑧, we want to choose a cycle basis 𝐵that minimizes the log width (4). In a kind of pre-test, we compare several cycle bases: 1. a fundamental cycle basis obtained from a minimum spanning tree in 𝐺w.r.t. 𝑤, 2. an undirected cycle basis minimizing the span 𝑢−𝓁, 3. an undirected cycle basis minimizing the number of arcs, EURO Journal on Transportation and Logistics 11 (2022) 100081 6 N. Lindner and C. Liebchen Table 1 Span analysis of the tension-restricted R1L1 and R4L4 scenarios. 𝐼 𝐿 𝑈 Fixed activities Free activities Avg. span Avg. wt. span R1L1 FIX FIX 55.18% 28.03% 18.98 6 511 MED FIX 53.78% 28.05% 19.76 7 608 ORIG FIX 53.78% 28.07% 20.49 8 646 FIX MED 10.34% 29.13% 21.03 15 285 MED MED 10.12% 29.24% 21.81 16 382 ORIG MED 10.12% 29.79% 22.55 17 420 FIX ORIG 10.34% 35.07% 26.19 35 391 MED ORIG 10.12% 35.61% 26.96 36 488 ORIG ORIG 10.12% 44.28% 27.70 37 526 R4L4 FIX FIX 42.78% 32.38% 24.99 3 802 MED FIX 42.66% 32.40% 25.76 4 073 ORIG FIX 42.66% 32.54% 26.46 4 320 FIX MED 8.93% 34.36% 27.02 7 643 MED MED 8.86% 34.70% 27.79 7 915 ORIG MED 8.86% 36.53% 28.50 8 162 FIX ORIG 8.93% 40.91% 31.84 16 221 MED ORIG 8.86% 42.01% 32.61 16 493 ORIG ORIG 8.86% 54.27% 33.32 16 740 Table 2 Evaluation of cycle bases for R1L1 and R4L4. The ‘‘fixed’’ column indicates the percentage of cycles 𝛾in the basis with 𝑧𝛾=𝑧𝛾. The log width is given w.r.t. decadic logarithms to make the numbers more intuitive. 𝐼Cycle basis ALL PARTIAL ODIJK Fixed Logwidth Fixed Logwidth Fixed Logwidth R1L1 Fundamental 33.28% 667 33.28% 1 516 0.00% 2107 min span 58.67% 343 58.67% 645 0.00% 1543 min #arcs 59.59% 346 59.59% 609 0.00% 1484 min 𝑐73.62% 216 73.62% 597 0.00% 2302 R4L4 Fundamental 22.01% 2 862 22.01% 6 542 0.00% 7 939 min span 37.66% 1 777 37.66% 3 399 0.00% 5 378 min #arcs 42.26% 1 722 42.26% 2 932 0.00% 5 000 min 𝑐59.46% 1 144 59.46% 2 886 0.00% 6 957 Table 3 Results for the 10 = 8 + 2 scenarios of round 1 for R1L1. 𝐿 𝑈 Best objective Better objectives Dual bound Gap FIX FIX 30415672 0 / 18 30415672 0.00% MED FIX 30415672 0 / 18 30415672 0.00% ORIG FIX 30415672 0 / 18 30415672 0.00% FIX MED 30415672 0 / 18 29673500 2.44% MED MED 30415672 0 / 18 28612511 5.93% ORIG MED 30036475 3 / 18 25883704 13.83% FIX ORIG 30415672 0 / 18 28697443 5.65% MED ORIG 30415672 0 / 18 25620519 15.77% Cycle offset strategy Best objective Better objectives Dual bound Gap ALL 30415672 0 / 18 30174317 0.79% PARTIAL 30415672 0 / 18 30118941 0.98% 4. an undirected cycle basis minimizing the function 𝑐as defined in Lemma 3.1. For our two instances, all the undirected cycle bases turn out to be integral, so that they are minimum integral cycle bases as well. Table 2 evaluates the fixed cycle offset variables and the log width of the aforementioned cycle bases. The cycle basis 𝐵⋆minimizing 𝑐 suggested in Section 3.3 comes out as a clear winner: It fixes by far the most variables and the log width is smallest for our strategies ALL and PARTIAL. It is clear by construction that the strategies ALL and PARTIAL fix the same set of cycle offset variables. However, it is remarkable that none of the cycle bases is able to fix an integer variable in the standard MIP formulation (1) just with Odijk’s cycle inequalities (2). This is due to the fact that every cycle in the R1L1 resp. R4L4 instance contains at least two transfer activities. The minimum span cycle basis is tailored to decrease the log width w.r.t. the Odijk bounds (Liebchen and Peeters,2009), but surprisingly, minimizing the number of activities in the cycle basis often produces an even smaller log width. Minimizing 𝑐is a bad choice for tightening Odijk’s bounds. In the sequel, we thus select MIP𝑆,𝐵∗(R1L1,ALL)and MIP𝑆,𝐵∗(R1L1,PARTIAL)for our computations with the two cycleoffset-restricted scenarios. 4.3. Computational results We finally turn to the timetables found by our merging approach. In this section, we omit the subscripts 𝑆and 𝐵∗, as we consider only a single 𝑆and a single 𝐵∗per instance. The results for R1L1 are summarized in Table 3 (Round 1) and Table 4 (Round 2). In Round 1, R1L1(ORIG,MED)is the only scenario where we could produce a better timetable (best weighted slack: 30036475) than the current PESPlib incumbent (weighted slack: 30415672). This underlines the ‘‘hardness’’ of the local optimum at the latter solution. The scenarios R1L1(−,FIX)can in fact be solved to EURO Journal on Transportation and Logistics 11 (2022) 100081 7 N. Lindner and C. Liebchen Fig. 2. Objective values of the solutions obtained from R1L1(ORIG,MED). After Round 1, runs #2, #3, and #5 produced better solutions than the current PESPlib incumbent (value 30415672, dotted blue line). After Round 2, in total 7 better timetables have been found. Fig. 3. Objective values of the solutions obtained from MIP(R1L1,ALL). No improvement upon the current PESPlib incumbent shows up in Round 1, but 7 better timetables are found after Round 2. Interestingly, these arise from weaker initial solutions or none at all. optimality: The current PESPlib incumbent is optimal. In particular, we cannot gain much information out of these 3 scenarios and discard them for Round 2. However, when reaching the computation time limit, several runs of Round 1 end up at a different timetable than the current PESPlib incumbent with only slightly higher slack. These timetables turn out to be ‘‘far enough’’ from the incumbent, and this is why we find improving solutions for all scenarios except 𝑈=FIX in Round 2. The best timetable has weighted slack 29894745 and is computed from one of the three better timetables produced by R1L1(ORIG,MED) in Round 1, but not from the one with lowest weighted slack. The PARTIAL cycle offset strategy produces the second best timetable with a weighted slack of 29907 781. The evolution of the objective values over both rounds for all 18 runs of the scenarios R1L1(ORIG,MED)and MIP(R1L1,ALL)is visualized in Figs. 2 and 3, respectively. These figures, as well as the subsequent ones, read as follows: The upper end of the green bar is the objective value of the initial solution from 𝑆that Round 1 has been fed with, or ∞in the case where no solution has been provided to Round 1 (∅). The border between the green and the yellow bars is the objective value of the best solution that has been found in Round 1, and is thus plugged into Round 2 as initial solution. Finally, the bottom of the yellow bar is the objective value that has been achieved as the result of Round 2. Table 4 Results of round 2 for R1L1. 𝐿 𝑈 Best objective Better objectives FIX MED 30348 574 2 / 18 MED MED 30003 486 4 / 18 ORIG MED 29 894 745 6/18 FIX ORIG 30373924 1/18 MED ORIG 30 335 565 3 / 18 Cycle offset strategy Best objective Better objectives ALL 29973362 7/18 PARTIAL 29907781 6/18 The picture for R4L4 is somewhat different (see Tables 5 and 6): Although the current PESPlib incumbent is the optimal solution to the R4L4(−,FIX)instances, we find at least 6 better timetables for each other scenario in Round 1, the best has weighted slack 37281703 and is found by R4L4(ORIG,MED), too. For Round 2, we again discard 𝑈=FIX. The remaining 7 scenarios bring plenty of better solutions, at least 12 per scenario. The best timetable is an outcome of a Round 1 timetable for R4L4(MED,ORIG), and has weighted slack 36729402. The objective values of the best timetables found in Round 2 are close for all scenarios, the second best with a weighted slack of 36753295 has again been produced by EURO Journal on Transportation and Logistics 11 (2022) 100081 8 N. Lindner and C. Liebchen Table 5 Results of round 1 for R4L4. 𝐿 𝑈 Best objective Better objectives Dual bound Gap FIX FIX 38381922 0/ 18 38381922 0.00% MED FIX 38381922 0/ 18 38381922 0.00% ORIG FIX 38381922 0/ 18 38381922 0.00% FIX MED 38095741 7/ 18 29918556 21.46% MED MED 37616308 6/ 18 24688127 34.37% ORIG MED 37281703 6/ 18 19414100 47.93% FIX ORIG 37862826 9/ 18 27834556 26.49% MED ORIG 37398748 9/ 18 23196127 37.98% Cycle offset strategy Best objective Better objectives Dual bound Gap ALL 37507260 9/ 18 30458434 18.80% PARTIAL 37499535 9/ 18 30560248 18.51% Fig. 4. Objective values of the solutions obtained from R4L4(MED,ORIG). Round 1 produces 6 timetables improving upon the current PESPlib incumbent (objective value 38 381 922, dotted blue line), and 3 more are found in Round 2. Runs #16-18 have been conducted without an initial solution, and the final objective value is larger than 41000 000. Table 6 Results of round 2 for R4L4. 𝐿 𝑈 Best objective Better objectives FIX MED 36784 153 18 / 18 MED MED 36772 886 12 / 18 ORIG MED 36757 228 13 / 18 FIX ORIG 36770965 12/18 MED ORIG 36 728 402 12 / 18 Cycle offset strategy Best objective Better objectives ALL 36775559 12/18 PARTIAL 36753295 13/18 the PARTIAL cycle offset strategy. The objective value evolution for R4L4(MED,ORIG)and MIP(R4L4,PARTIAL)is visualized in Figs. 4 and 5, respectively. To sum up, although our restricted scenarios do not always give rise to better solutions of the original instance in Round 1, the ConcurrentPESP solver is able to compute timetables of very good quality that can help overcome local optima on the original instance in Round 2. Restricting tensions produced lower objective values, but the difference to the cycle offset approach is only minor, so that we consider both approaches as fruitful. 5. Conclusion We presented two strategies to merge periodic timetables in the context of the Periodic Event Scheduling Problem. Both rely on mixed integer programming, one constructing restricted PESP instances with tighter bounds on periodic tensions, and the other one defining restricted mixed-integer programs tightening the bounds of the integer cycle offset variables. The new bounds are computed using a set of initial solutions, considering the span of the activities and using minimum cycle basis techniques. Varying the bounding schemes and the type of variables to restrict, we constructed in total 10 scenarios per original instance. These scenarios have been given as input to the PESP solver ConcurrentPESP with several combinations of initial solutions and parameter settings. On the smallest and largest PESPlib instance, this approach overcomes local optima for the classical PESP heuristics, and is therefore able to produce better periodic timetables. We are confident that the approach will turn out to be fruitful on the other PESPlib instances, and on general PESP instances as well: For public transport applications, there is often a natural distinction between activities of large span (e.g., transfers, turnarounds) and small span (e.g., driving, dwelling), so that the techniques from Section 3 apply. Although we did not experiment with different sets 𝑆of initial solutions, the outcome was satisfactory, given the hardness of periodic timetabling problems. It would be interesting to integrate the search for these sets 𝑆into a merging framework, perhaps guided by reinforcement learning. A related question is whether a statistical tuning of parameters (e.g., 𝛼, or parameters of the concurrent solver) could help improving the method. The main difficulty here is that it takes a large amount of time to produce high-quality solutions. Another direction of investigation could be to combine our two directions, i.e., restricting periodic tension and cycle offset variables simultaneously. Bearing similarities to the crossover heuristic as it is applied for MIPs, the basic principle of our method is certainly applicable to general mixed integer programs. Note however that our computation times are several hours, whereas standard crossover for solution polishing is