Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (2024)

1]HUN-REN Wigner Research Centre for Physics, Konkoly-Thege Miklós út 29-33., Budapest, 1121,Hungary2]Faculty of Transport and Aviation Engineering, Silesian University of Technology, Akademicka 2A, Gliwice, 44-100, Poland3]Institute of Theoretical and Applied Informatics, Polish Academy of Sciences, Bałtycka 5, Gliwice, 44-100, Poland4]Joint Doctoral School, Silesian University of Technology, Akademicka 2A, Gliwice, 44-100, Poland5]“Friedrich List” Faculty of Transport and Traffic Sciences, Technical University of Dresden,Dresden, 01069, Germany

Mátyás Koniorczykkoniorczyk.matyas@wigner.hun-ren.hu Krzysztof Krawieckrzysztof.krawiec@polsl.pl Ludmila Botelholbotelho@iitis.pl Nikola Bešinovićnikola.besinovic@tu-dresden.de Krzysztof Dominokdomino@iitis.pl[[[[[

Abstract

We address the applicability of hybrid quantum-classical heuristics for practical railway rescheduling management problems.We build an integer linear model for the given problem and solve it with D-Wave’s quantum-classical hybrid solver as well as with CPLEX for comparison.The proposed approach is demonstrated on a real-life heterogeneous urban network in Poland, including both single- and multi-track segments and covers all the requirements posed by the operator of the network.The computational results demonstrate the readiness for application and benefits of quantum-classical hybrid solvers in the realistic railway scenario: they yield acceptable solutions on time, which is a critical requirement in a rescheduling situation. At the same time, the solutions that were obtained were feasible. Moreover, though they are probabilistic (heuristics) they offer a valid alternative by returning a range of possible solutions the dispatcher can choose from. And, most importantly, they outperform classical solvers in some cases.

Keywords

railway rescheduling, \sepconflict management, \sepheterogeneous urban railway network, \sepquantum annealing, \sephybrid quantum-classical heuristics

1 Introduction

Rail transport is expected to experience an increase in capacity demands due to changes in mobility needs resulting from climate policy, leading to traffic challenges in passenger and cargo rail transport. The situation is aggravated by the fact that the rolling stock and especially railway infrastructure cannot keep up with the increase in transport needs, which makes railway systems overloaded.

Rail transport, due to its technical and organizational characteristics is very sensitive to disturbances in traffic: these extraordinary events have an impact on railway operations, typically resulting in delays [23]. Examples of such disturbances include: late train departures/arrivals, extended dwell times, or (partial) track closures (also referred to as disruptions). These can last from several minutes up to hours.The impact of disturbances can propagate to multiple sections in the railway network (cf. the work of [50] and references therein).Thus, ensuring stable railway traffic and providing reliable service for passengers, rail cargo companies, and their clients is in the best common interest of railway infrastructure managers, and train operators.To limit these impacts as much as possible it is necessary to make proper rescheduling decisions quickly. The dispatchers need to reschedule and partially reroute the trains, aiming at the minimization of negative consequences. Still, in many places, dispatchers are using their own intuition or simple heuristics like FCFS (First Come First Served); resulting in decisions far from the objective of consequence minimization. In the last decades, there has been a growing research interest in mathematical optimization methods in support of rail dispatchers in decision making. While doing so, diverse objective functions were addressed, such as the (weighted) sum of delays [36], the maximal delay that cannot be avoided[14], or fuel consumption measures [29].

As a railway network is a complex non-local structure, modeling a bigger portion of it is necessary for efficient suppression of the consequences of disturbances. Hence, large-scale rescheduling problems have to be addressed, and as the time available for making decisions is limited, they have to be solved almost real time. Therefore it is vital to develop efficient algorithms. However, the current conventional rescheduling optimization models have difficulties addressing complex and large-scale instances in suitable computation time [8, 43, 35].

The railway dispatching/rescheduling problem can broadly be recognized as being equivalent to job-shop scheduling with blocking and no-wait constraints[49, 41]. A possible modeling approach is based on alternative graphs employing order and precedence variables[36], facilitating the formulation as a mixed integer linear program which is often large, hence, specialized algorithms are often used instead[35, 14]. An alternativeapproach is to use time-indexed models: discrete-time units and binary decision variables that assign events to particular time instants. Even though this results in very large problems, it is applied for both timetabling and dispatching/rescheduling[8, 39, 43, 18].

As time indexed models result in integer or 0-1 programs, they are suitable[53] for to be solved on new types of hardware: quantum annealers[2, 33]. These physical devices can be considered as stochastic heuristic solvers for Quadratic Unconstrained Binary Optimization (QUBO) problems, and there are few commercially available options, including D-Wave[31]. The exploration of their potential applications attracts a growing research interest (c.f. Section2.2), in which railway applications are not yet strongly represented.

In this paper, we address the train rescheduling problem in complex railway networks with mixed infrastructure including single, double, and multiple track railway lines with given planned train paths. Our consideration includes shunting movement of rolling stock between depots and stations followed by rolling stock connections.We apply a new hybrid rescheduling algorithm combining classical-quantum modelling and based on Quantum Annealing (QA). This work extends ona particular linear modeling strategy, partly explored on a toy model in[21].We develop an Integer Linear Programming (ILP) model, which is solved with proprietary D-Wave solvers as well as with CPLEX for comparison.The D-Wave approach is analysed in detail, with respect to its applicability and performance in a practically relevant situation.The results suggest that quantum computing and QA in particular, although an early-stage technology, are ready for tackling challenging railway rescheduling problems. During the preparation of the present work, another contribution has appeared in the literature[56], which is similar to our in that it sets up a QUBO model for a somewhat similar problem, and solves it by a simulation of (quantum) annealing. This contribution, however, considers a different railway infrastructure: high-speed trains, and considers timetable optimization, for which it employs a model different from the one described here.

The main contributions of the present work are the following:

•
A railway traffic dispatching model is introduced for rescheduling trains due to disturbances and disruptions taking into account given timetable, rolling stock connections and a network with single or multi-track segments.
•
QA-based hybrid heuristics are, for the first time, applied to real-life problems of railway rescheduling optimization on urban railway networks, these heuristics have some but still meaningfully input of real QA.
•
The performance of the hybrid (quantum classical) D-Wave solver is demonstrated on a real-life network in Poland and diverse disturbance scenarios of different sizes and types.
•
The D-Wave hybrid solver provides good quality solutions in short time, also for the complex instances in limited computational time.
•
For certain instances the hybrid (quantum classical) D-Wave solver outperforms classical solvers in terms of computational time, yielding still feasible solutions.

The paper is organized as follows. In Section2 we review state-of-the-art literature and identify scientific gaps we intend to fill, in Section3 we describe problem under investigation, in Section4 we evaluate our model, in Section5 we discuss the hybrid solver we use, in Section6 we present computational results.

2 Literature review

In this section, we present time-indexed models known from the literature, comparing them briefly with our modelling approach.For a general review of railway timetabling and rescheduling we refer to the reviews by [38] and [7].Second, we review recent QA applications in optimization in general, and also in the rail/transport domain.Finally, we summarize the existing scientific gaps.

2.1 Classical railway dispatching/rescheduling

Time-indexed modelling of railway scheduling and rescheduling is quite common in literature for routing [39, 8] and scheduling [9, 29, 43, 13] .[8] focuson determining train’s path for a complex central railway station area(called the blocking - stairway), used a discrete-time model resulting in a $0-1$ program. Their model is successfully demonstrated on an operational day at the central railway station area Berne, Switzerland.Similarly, [39] consider train movements in a single major (and thus complex) railway junction, including a freight yard and a few minor stations. They build a mixed integer programming model that can also be solved close to optimal in case of practically relevant instances. They address the train movements on the level of detail of acceleration and deceleration strategies, hence adopting a significantly finer discrete time scale than our one-minute resolution.When compared to our problem, these two models addressed a high traffic in a complex station area. Instead, we are concerned more with a bigger area, multiple stations, and mixed track, but a lower train density and complexity.Our train routes fixed within the stations, and we also model shunting movementsthat is not a subject of the cited reference.

[9] proposed a graph theoretic formulation for the problem using a directed multigraph in which nodes correspond to departures/arrivals at a certain station at a given time instant. This formulation is used to derive an integer linear programming model that is solved using a Lagrangian relaxation. [13] introduced a new pure 0-1 programming formulation, and call it Tick Formulation, to model the Deadlock Detection (DD) problem.[43] focused “on the train dispatching problem on an N-track network, with the main challenge of how to formulate specific retiming, reordering, retracking and rerouting options in combination”. Their model was demonstrated on generic networks of different structure than what we address. In our situation we are limited to retiming and reordering of trains.Structurally, their model shows similarities to ours, including the role of binary variables and also the introduction of "cells" that generalize the objects the occupancy of which can be controlled: they can be one or more blocks, stations, etc. In our model, the decision variables will be explicitly linked to a subset of stations.

2.2 Quantum Annealing and its applications

Quantum Annealing (QA)[15] is an optimization methodanalogous to Simulated Annealing (SA). Both can be used for theunconstrained minimization of a quadratic function of $\pm 1$ -valued variables, which represents the energy configuration for a set of spins of the Ising model[30, 4]; an important andcelebrated problem in physics. The Ising problem can be equivalentlyrewritten in a form of a Quadratic Unconstrained Binary Optimization(QUBO) problem with binary ( $0-1$ ) variables.The Max-Cut problem[5] is also equivalent to the Ising model or QUBO. The problem is NP-hard, has tremendous literature in classical optimization, andthere is significant ongoing progress in the development ofalgorithms to solve it[22, 28]. An extensive list of Ising problems and their formulations are presented e.g. in the work of [37].

Hardware quantum annealers, like the D-Wave machine, implement a quantum version of the Ising model, assigning a quantumbit, i.e. a two-level physical quantum system to each bit. The system employs real physical two-level quantum systems with a tunable energyoperator (Hamiltonian), the energy being the objective function of theoptimization problem.An adiabaticevolution of the quantum system is implemented: the interactionsbetween the spins are slowly changed from an initial Hamiltonian witha simple minimal-energy state to the Hamiltonian corresponding to theobjectiveaccording to the rules of quantummechanics. According to the adiabatic theorem of quantummechanics[3], under certain conditions thephysical system remains in its lowest energy state during theevolution, and an optimal, minimum-energy solution can be possibly read out atthe end.

Anideal adiabatic quantum computer would be a device with a large numberof perfect quantum bits completely decoupled from their environment,very close to zero temperature. Any pair of qubits could be coupled, and even, more generalcouplings, e.g. involving $3$ or more spins could take place. Such asystem would exhibit perfect quantum coherence and entanglement.Even in such an ideal system the evolution time to reachat a minimum energy state with high enough probability is inverse proportional to thegap between the minimal energy configuration and the one closest to it. Certainly the gap is not known in advance, hence, the annealing time has to be determined by experimenting in practice. A more severe issue is that the gap can bevery small, hence, it is not possible tosolve all possible hard problems on such a setup.In spite of that, there are problems in which this approach can be efficient.

Meanwhile, the state-of-the-art physical quantum annealers are smallersystems of a few hundred quantum bits in which not all qubits arecoupled, and the arrangement operates at a finite temperature. Thefixed topology of qubit couplings means that the problem’s graphdefined by the nonzero coupling has to be embedded as an inducedsubgraph of the topology of the system. This embedding is a hardproblem itself[59]. It often requires tocouple multiple physical qubits to represent a logical bit of theproblem. As for the finite temperature, the adiabatic evolution worksalso in systems that interact with theirenvironment[54]. This introduces noise, whoseimpact, as discussed in detail byAlbash\BOthers. [1],increases with the system size.

Operationally, from the optimizer’s point of view quantum annealerscan be viewed as probabilistic heuristic solvers returning astatistical sample[20] of configurationswhich are supposed to be optimal or close to optimal. It must bestressed that the quantum annealers are not algorithms running ondigital computers; they are analog devices implementedphysically. This latter implies that the coefficients of the problemare encoded with a limited accuracy, and the algorithmic properties ofthe particular optimization problem such as its complexity class willnot determine the solver’s behaviour.

To overcome the problem of limited size and accuracy, quantumannealers of the present state of the art are often used in hybrid(quantum classical) solvers, orchestrating classical algorithms andusing QA as a subroutine in order to address hardproblem instances more efficiently. Solvers available in D-Wave’s ’LeapHybrid Solver Service (HSS)’[11], including the one usedin the present work, belong to this family. Meanwhile, quantumtechnology keeps on developing, systems of bigger size and bettertopology are regularly announced, they are more and more affordable,and there is a growing community around them.

Currently, the applicability of QA technologies is being explored. Benchmark problems are solved[55, 42], comparisons with digital computers aremade[32], etc.The application of quantum annealers in transportation optimization is a new area with only a few contributions so far.There is an apparent interest in this research direction in the aviation domain[47, 51, 48], which relates also to the already mentioned benchmarking of annealers as done by[55]. Shipment rerouting was also addressed in the context of QA[58].There are a few examples of hybrid solvers’ applications, too. These include traffic flow optimization [44], multi-car paint shop optimization [57], tail assignment problem [40],and vehicle optimization [24].

In railway context, the first applications of QA can be found in train dispatching/rescheduling [19, 21] and rolling stock planning[27]. In particular, the first proof-of-conceptdemonstration of a (pure) quantum computing approach to railway rescheduling was presented by some of us[18]. A followup paper[21] laid down the principles of a more general modeling approach to railway rescheduling in light of QA, introducing a suitable QUBO / HOBO (Higher-Order Binary Optimization) encoding of these problems. The present work continues this line of research.

2.3 Scientific gaps

We recognize several important scientific gaps. First, existing optimization models suffer from the curse of dimensionality, and inability to solve larger real-life instances.Second, the QA models for the studied optmimzation problem that have been introduced so far [19, 21] had been designed for pure QA implementations and have been demonstrated only on simple network setups due to the size limitations of the currently available Noisy Intermediate Scale Quantum (NISQ) devices.Third, no hybrid QA-based models have been used for real-life railway, or even other schedule-based modes, like public transport and air traffic planning and/or rescheduling.

In this paper we demonstrate the quantum readiness of medium-scale railway rescheduling models: we successfully apply quantum methods in the rescheduling problem of a urban railway networks. To do so, we use a hybrid approach combining quantum and classical computing. For a fair comparison of current and future QA approaches with state-of-the-art classical approaches, we elaborate the railway model readily suitable for both classical and hybrid (quantum classical) solvers in a comparable manner. Quantum computing is believed to develop rapidly in near future. Hence, demonstration quantum readiness for railway problems right now is step forward towards the demonstration of quantum advantage in future. On future quantum devices the comparison with classical solvers is expected to be in strong favour for the quantum ones. In this way we also contribute in developing a novel general set of railway rescheduling models that has a potential to scale well with the size of the problem, aiming to overcome the curse of dimensionality.

3 Problem description

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (1)

We consider a railway rescheduling problem that includes train reordering, retiming, and shunting movements with rolling stock circulation at stations. Also, we consider an urban railway network with mixed tracks from single, double, up to quadruple tracks.

We model a railway network with edges and nodes of a graph.Figure 1 depicts an exemplary network layout composed of nodes (station or junction)and edges (single, double, multiple track lines).Each edge is composed of one or more tracks. Each track consists ofblocks sections defined by pairs of signals [39]. Each block section can be occupied by at most one train at a time.A subsequent train is allowed to enterthe block only after a minimal headway: minimal time spanbetween the trains. We assume a 2-block signalling system, meaning that two free blocks are required between the consecutive trains. We also assume green way policy (trains have the free way to move at the maximal allowed speed between stations), resulting in constant running times for each block.

Each node is composed of station tracks (blocks) and interlocking areas.One station track can be occupied by one train at a time.Routing dependencies between pairs of trains competing for the same resource in interlocking areas are considered in order toto guarantee that only one train from the pair can occupy the area at the time.

A train’s route is the sequence of blocksthe train passes during its journey.A train path is a sequence of arrival and departure times of a particular train assigned to a train route.We assume a one-minute resolution for all time parameters such as timetable time, running and dwell time minimal headway time or passing time, resulting in integer variables.This assumption of discretized time is needed to enable the use of quantum-based solvers, and, more importantly, introduces the possibility of pruning the inherently binary variables.The relaxation of integer constraints on time variables is not expected to improve computation time significantly, as the real difficulty is tied to precedence of trains, encoded later onprecedence binary variables.

In the network, two trains can follow each other, i.e. keep the given order, meet and pass (M-P) when going the opposite directions, and meet and overtake (M-O), i.e. change the order, when going the same direction. Single tracks are designed for bidirectional traffic, double-tracks for unidirectional, one for each direction, and multi-track lines can be combination of unidirectional and bidirectional.On single-track lines, M-Ps and M-Os are only possible at stations. To prevent M-P on single-track lines, we determine the set $\mathcal{J}^{2}_{\text{single}}$ as the setof all pairs of trains that can potentially meet on the single trackline heading in opposite directions.

Trains in the same direction preserve their order between stations, and keep the minimal headway time between each other. To prescribe these we determine set $\mathcal{J}^{2}_{\text{headway}}$ : the set of all pairs of trains that can potentially violate the minimal headway condition.

The train traffic is scheduled on the basis of the given timetable. The timetable contains all the train paths.In the train rescheduling problem, we are given an initial timetable that is conflict-free. However, conflicts may appear due to disturbances such as late departures and/or arrivals due to excessive passenger demand, malfunctioning rolling stock, etc.Following [10], a conflict is an inadmissible situation when at least a pair of trainsclaim the same resource (e.g., block section, switch) simultaneously. We assume that the conflict may occur either on the railway line or at the station.Possible conflicts that can occur on the line include the lack of the minimum headway between two subsequent trains heading on the same track in the same direction, or claiming a segment of a single track line by two trains heading in opposite directions at the same time.Conflicts at stations include claiming a station track by two trains at the same time, claiming a station switch (in the interlocking area) by two trains at the same time.

The conflicts have to be solved by modifying the original timetable, applying decisions on the train sequencing and retiming for trains claiming for the same resource.Such an intervention in the structure of train paths implies additional changes to maintain the feasibility of the modified timetable.These changes typically result in delays of additional trains that were not directly involved in the conflict otherwise, giving rise to secondary delays depending on the dispatching decision depending on the modification of the timetable. It is common in the literature to minimize a function of these secondary delays (see e.g. the work of [14]); we will do this in the present contriubution, too).

An important element of our approach is that from amongst all stations we first determine a subset of the decision stations; we assume that the direct decisions implied by our model affect thesestations only.Hence, as decision stations we select those stations, where routes of trains intersect, where trains start or terminate, or where selected part network is bounded. The motivation is that if the routes of trains are fixed, the decisions on modifying their train paths has to be made with respect to these stations. (If we change some of the trains’ routes in a re-routing process, new model with new parameters if developed, e.g. additional decision stations may appear.)In what follows, by a station we always mean a decision station.This also means that non-decision stations appear only through the parameter values in the model, they do not appear as indices of decision variables or parameters. Headways, for instance, are calculated between decision stations, taking into account all the line blocks and station blocks of non-decision stations in between.

4 Methods and Model

In the following we describe our model in detail.Section 4.1 defines sets, parameters and decision variables. Section 4.2 describes our Integer Linear Programming (ILP)formulation.

4.1 Sets, parameters and decisions

To formulate our decision variables, constraints and the objectivefunction, we determine sets of index tuples needed to find the actual index sets of variables, from the given infrastructuredata, timetable data, and the rolling stock circulation plan. Also, weintroduce parameters calculated from the same input.

4.1.1 Sets

Let us denote by $\hat{S}$ the set of all stations, and by $S\subset\hat{S}$ the decision stations we have determined in advance, as already described.The main objects of our model for railway rescheduling are thetrains $j\in\mathcal{J}$ and the stations $s\in\mathcal{S}_{j}$ intheir path. Set $\mathcal{S}_{j}$ includes decision stations only and it is an ordered set.In addition to $\mathcal{J}$ and $\mathcal{S}_{j}$ , the relevant sets are the following.Set $\mathcal{J}^{2\,(\text{turn})}_{s}\subset\mathcal{J}\times\mathcal{J}\left(%\forall s\in\mathcal{S}\right)$ is the set of all pairs of trains so that the firsttrain of the pair terminates at station $s$ and its rolling stockcontinues as the second train of the pair. This set is deducedfrom the rolling stock circulation plan and the timetable. Set $\mathcal{J}^{2\,(\text{close})}\subset\mathcal{J}\times\mathcal{J}$ is the set of trains which are close enough to each other in time so that precedence variables have to be defined for them. More details will be given at the description of the parameter $d_{\max}$ which this set depends on.Set $\mathcal{J}^{2\,(\text{headway})}\subset\mathcal{J}^{2\,(\text{close})}$ is the set of train pairs that can potentially violate theminimal headway time between two subsequenttrains on a line segment. For a pair of trains to be in $\mathcal{J}^{2\,(\text{headway})}$ , both their routes need to includethe same line segment so that the trains are moving inthe same direction on it and can meet there according to model parameters, i.e. maximal allowed delay.Set $\mathcal{J}^{2\,(\text{single})}\subset\mathcal{J}^{2\,(\text{close})}$ is the set of all trains that share at least one single-track linesegment as a common part of their routes so that they are heading inthe opposite direction.Set $\left(\forall s\in\mathcal{S}\right)\ \mathcal{J}^{2\,(\text{track})}_{s}%\subset\mathcal{J}^{2\,(\text{close})}$ is the set of train pairs that areplanned to occupy the same track on station $s$ anytime during theplanning horizon, thereby competing for the same track.Set $\left(\forall s\in\mathcal{S}\right)\ \mathcal{J}^{2\,(\text{switch},\text{out%})}_{s}\subset\mathcal{J}^{2\,(\text{close})}$ is the set oftrain pairs that are planned to pass the same interlocking area of $s$ upon theirdeparture from station $s$ .Set $\left(\forall(s,s^{\prime})\in\mathcal{S}^{\times 2}\right)\ \mathcal{J}^{2\,(%\text{switch},\text{out},\text{in})}_{s,s^{\prime}}\subset\mathcal{J}^{2\,(%\text{close})}$ is theset of train pairs that are planned to pass the same interlocking area of $s$ in order to have $j$ to departure $s$ while $j^{\prime}$ arrive $s$ fromthe direction of $s^{\prime}$ .Set $\left(\forall s\in\mathcal{S}^{\times 2}\right)\ \mathcal{J}^{2\,(\text{switch%},\text{in},\text{noMP})}_{s,s^{\prime}}\subset\mathcal{J}^{2\,(\text{close})}$ is theset of train pairs that are planned to pass the same interlocking area of $s$ upon their arrival at $s$ from the direction of $s^{\prime}$ , and there is no M-Ppossibility for them between $s$ and $s^{\prime}$ Set $\left(\forall s\in\mathcal{S}^{\times 2}\right)\ \mathcal{J}^{2\,(\text{switch%},\text{in},\text{MP})}_{s,s^{\prime}}\subset\mathcal{J}^{2\,(\text{close})}$ is the setof train pairs that are planned to pass the same interlocking area of $s$ upon their arrival at $s$ so that either they both come from the directionof $s^{\prime}$ but there is a M-P possibility for them between $s$ and $s^{\prime}$ , or one of them is approaching $s$ from a directionother than $s^{\prime}$ .Set $\left(\forall j\in\mathcal{J}\right)\ \mathcal{C}^{2}_{j}$ isthe set of all station pairs that are subsequent in the route of $j$ .Set $\left(\forall(j,j^{\prime})\in\mathcal{J}^{\times 2}\right)\ \mathcal{C}^{2\,(%\text{common})}_{j,j^{\prime}}$ is the set of all stationpairs that appear as subsequent stations in the route of both $j$ and $j^{\prime}$ heading in the same direction.Set $\left(\forall(j,j^{\prime})\in\mathcal{J}^{\times 2}\right)\ \mathcal{C}^{2\,(%\text{common, single})}_{j,j^{\prime}}$ is the set of allstation pairs that appear as subsequent stations in the route ofboth $j$ and $j^{\prime}$ which are connected with a single-track linesegment, and trains are heading in opposite directions. The order within the pairs is determined by $j$ .All these sets can be enumerated on the basis of the input data in astraightforward manner.

4.1.2 Parameters

The parameters that appear in our model are the following:

•
$\tau^{(\text{pass})}(j,s\rightarrow s^{\prime})$ is the running time of $j$ from $s$ to $s^{\prime}$ .
•
$\tau^{\text{(headway)}}(j,j^{\prime},s\rightarrow s^{\prime})$ is theminimal headway time for $j^{\prime}$ following $j$ from $s$ to $s^{\prime}$ .
•
$\tau^{(\text{switch})}(j,j^{\prime},s)$ is the running time of $j$ over a interlocking area of $s$ , where $j$ may be in conflict with $j^{\prime}$ . (In our examples we do not consider cases when, e.g., on bigger stations there are multiple switches with different technological times. The model could trivially be extended to cover such scenarios by adding extra indices if relevant.)
•
$\tau^{(\text{dwell})}(s,j)$ is the minimal dwell time of $j$ at $s$ .
•
$\tau^{\text{(turn)}}(s,j,j^{\prime})$ is the minimum turnaround time for therolling stock of a train $j$ terminating at $s$ to continue itsjourney as train $j^{\prime}$ .
•
$\sigma(j,s)$ is the scheduled departure time of $s$ from $j$ , given by the original timetable.
•
$\upsilon(j,s)$ is the earliest possible departure time of $j$ from $s$ when a disturbance and/or disruption in the network happens. It is the maximum of $\sigma(j,s)$ and the technically feasibleearliest departure time, and may be more constraining than the first depending on the initial delays.The latter is calculated assuming that the train follows its planned route at minimum running time, and that there are no other trains present on the network.
•
$d_{\text{max}}$ is an upper bound assumedfor thesecondary delay $t^{(\text{out})}(j,s)-\upsilon(j,s)$ . This parameter sets an upper limit for the possible secondary delays that can arise on the particular part of the network. Setting such a bound is common in the literature, see e.g. the work of [14]. We set the same value of this parameter for all our testing instances. It has to be big enough so that no secondary delay exceeds it (e.g. $d_{\max}=40$ excludes the possibility of a one-hour delay due to waiting for another delayed train).Meanwhile it is desirable to set its value as small as possible; restricting the time variables to small intervals decreases the number of binary decision variables in the model.

The choice of a small $d_{\max}$ parameter results in a decrease of the model size, which is achieved via defining the set $J^{2\,(\text{close})}$ . A pair of trains $(j,j^{\prime})$ is included in this set if and only if they can meet on any station, given the original timetable, the disturbed/disrupted timetable, and assuming that no train can have a secondary delay greater than $d_{\max}$ . The number of constraints will be proportional to the average number of trains that can meet another train at a station; this will be linear in $d_{\max}$ .

4.1.3 Decision variables

Our decision variables are the following.First, we use departure time variables:

t^{(\text{out})}(j,s)\in\mathbb{N}

(1)

defining the departure time of train $j\in\mathcal{J}$ from station $s\in\mathcal{S}_{j}$ . Such a variable is defined for all decisionstations $\mathcal{S}_{j}$ . In the case of trains that terminate within the modeled part of the network, the last station is not included in $\mathcal{S}_{j}$ .The arrival time of the trains at stations, $t^{(\text{in})}(j,s)\in\mathbb{N}$ are trivially related to the $t^{(\text{out})}(j,s)$ variables through a constant offset (c.f. Eq.(3)), given the fixed running time assumption.

Second, in addition to the time variables, we use three sets of binary precedence variables. Decision variables $y^{(\text{out})}(j,j^{\prime},s)\in\{0,1\}$ determine the order of trains to leave stations: the variable takes the value $1$ if train $j$ leaves station $s$ before train $j^{\prime}$ , and 0 otherwise. The $(j,j^{\prime},s)$ tuples for which an $y$ variable is definedwill be specified later.Similarly, the precedence variables $y^{(\text{in})}(j,j^{\prime},s)\in\{0,1\}$ prescribe the order at the entry to stations; the value is $1$ if $j$ arrives to $s$ before $j^{\prime}$ .The last set of precedence variables describes the precedence of trains at some resource located between stations $s$ and $s^{\prime}$ (e.g. single track used by $j$ and $j^{\prime}$ heading in opposite direction). The binary variable $z(j,j^{\prime},s,s^{\prime})\in\{0,1\}$ will be $1$ if $j$ uses the given resource before $j^{\prime}$ . Also in thiscase, the quadruples $(j,j^{\prime},s,s^{\prime})$ for which we have such a variablewill be specified later.

4.2 ILP formulation

Given the index sets, variables, and the parameters and decision variables of the model, now we formulate the optimization objective and the constraints.

4.2.1 Objective function

Our goal is to minimize the secondary delays that occur in the analyzed part of the railway network.Hence, a suitable objective function is the weighted sum of secondarydelays at the destination station:

f(t)=\frac{1}{d_{\text{max}}}\sum_{j\in\mathcal{J}}w(j)\left(t^{(\text{out})}(%j,s^{*})-\upsilon(j,s^{*})\right).

(2)

where $s^{*}$ is the last element of $S_{j}$ , and $w_{j}$ are weights for each train representing its priority.The reason for restricting ourselves to the secondary delays at the destination stations only is that this is the actual figure of merit chosen by the respective railway authorities; the generalization to other linear objectives such as the weighted sum of secondary delays on all stations is straightforward.Theconstant multiplier $1/d_{\text{max}}$ is optional; we use itfor better comparability of different instances.

4.2.2 Constraints

The constraints of the model are the following.

Minimal running time

Each train needs a minimal time to get to the subsequent station:

\forall j\in\mathcal{J}\ \forall(s,s^{\prime})\in\mathcal{C}_{j}\quad t^{(%\text{in})}(j,s^{\prime})=t^{(\text{out})}(j,s)+\tau^{(\text{pass})}(j,s%\rightarrow s^{\prime})

(3)

Headways

A minimal headway time is required between subsequent train pairs on the common part of their route as

	$\displaystyle\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{headway})}\ %\forall(s,s^{\prime})\in\mathcal{C}^{2\,(\text{common})}_{j,j^{\prime}}t^{(%\text{out})}(j^{\prime},s)\geq$
	$\displaystyle t^{(\text{out})}(j,s)+\tau^{(\text{headway})}(j,j^{\prime},s%\rightarrow s^{\prime})-C\cdot y^{(\text{out})}(j^{\prime},j,s),$		(4)

where $C$ is a constant big enough to make the constraint satisfiedwhenever the binary variable $y^{(\text{out})}(j^{\prime},j,s)$ takes thevalue of $1$ . (Constraints of this structure are often termed as big-M constraints in this context.) In our implementation we calculate and use thesmallest suitable value of $C$ given particular $d_{\max}$ value, i.e.

C=-\min(t^{(\text{out})}(j^{\prime},s))+\max(t^{(\text{out})}(j,s))+\tau^{(%\text{headway})}(j,j^{\prime},s\rightarrow s^{\prime})

(5)

Single-track occupancy

Trains moving in opposite directions cannot meet on the same single-track line segment:

	$\displaystyle\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{single})}\ \forall%(s,s^{\prime})\in\mathcal{C}^{2\,(\text{common, single})}_{j,j^{\prime}}$
	$\displaystyle t^{(\text{out})}(j^{\prime},s^{\prime})\geq t^{\text{(in)}}(j,s^%{\prime})-C\cdot z(j^{\prime},j,s^{\prime},s),$		(6)

where $C$ is a big enough constant chosen similarly to that in Eq.(4).

Minimal dwell time

Each train has to occupy the station node for a prescribed time duration at each station:

\forall j\in\mathcal{J}\forall s\in S_{j}\quad t^{(\text{out})}(j,s)\geq t^{(%\text{in})}(j,s)+\tau^{(\text{dwell})}(j,s).\\

(7)

Timetable

No train is allowed to depart before its scheduled departure time:

\forall j\in\mathcal{J}\forall s\in S_{j}\quad t^{(\text{out})}(j,s)\geq\sigma%(j,s).

(8)

Station track occupancy

Station tracks can be occupied by at most one train at a time:

\small{\forall s\in S\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{track})}_{%s}}t^{(\text{in})}(j^{\prime},s)\geq t^{(\text{out})}(j,s)-C\cdot y^{(\text{%out})}(j^{\prime},j,s),\\

(9)

where $C$ is chosen similarly to that inEq.(4) again. Note that this requirementmay not be needed for depot tracks; the exceptions can be handledby the proper definition of $\mathcal{J}^{2\,(\text{track})}_{s}$ .

Interlocking area occupancy

These ensure that trains cannot meet in interlocking areas:

	$\displaystyle\forall s\in S\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{%switch},\text{out})}_{s}t^{(\text{out})}(j^{\prime},s)\geq$
	$\displaystyle t^{(\text{out})}(j,s)+\tau^{(\text{switch})}(j,j^{\prime},s)-C%\cdot y^{(\text{out})}(j^{\prime},j,s),$		(10)
	$\displaystyle\forall(s,s^{\prime})\in S^{\times 2}\forall(j,j^{\prime})\in%\mathcal{J}^{2\,(\text{switch},\text{out},\text{in})}_{s,s^{\prime}}t^{(\text{%in})}(j^{\prime},s)\geq$
	$\displaystyle t^{(\text{out})}(j,s)+\tau^{(\text{switch})}(j,j^{\prime},s)-C%\cdot z(j^{\prime},j,s^{\prime},s),$		(11)
	$\displaystyle\forall(s,s^{\prime})\in S^{\times 2}\forall(j,j^{\prime})\in%\mathcal{J}^{2\,(\text{switch},\text{in},\text{noMP})}_{s,s^{\prime}}t^{(\text%{in})}(j^{\prime},s)\geq$
	$\displaystyle t^{(\text{in})}(j,s)+\tau^{(\text{switch})}(j,j^{\prime},s)-C%\cdot y^{(\text{out})}(j^{\prime},j,s^{\prime}),$		(12)
	$\displaystyle\forall(s,s^{\prime})\in S^{\times 2}\forall(j,j^{\prime})\in%\mathcal{J}^{2\,(\text{switch},\text{in},\text{MP})}_{s,s^{\prime}}t^{(\text{%in})}(j^{\prime},s)\geq$
	$\displaystyle t^{(\text{in})}(j,s)+\tau^{(\text{switch})}(j,j^{\prime},s)-C%\cdot y^{(\text{in})}(j^{\prime},j,s),$		(13)

with a choice of $C$ similar again to that inEq.(4).

Rolling stock circulation constraints

These are introduced in order to bind the train with the shunting movement called also service train. If the train set of train $j$ which terminates at $s$ is supposed to continue its trip as (service train) $j^{\prime}$ or vice versa, a precedence of these trains including a minimum turnaround time has to be ensured:

\forall s\in S\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{turn})}_{s}t^{(%\text{out})}(j^{\prime},s)\geq t^{(\text{in})}(j,s)+\tau^{\text{(turn)}}(j,j^{%\prime},s).

(14)

Order of trains

We have additional conditions on the $y$ -variables, concerning the case when M-O is not possible on the line or station.

\begin{split}\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{headway})}\ %\forall(s,s^{\prime})\in\mathcal{C}^{2\,(\text{common})}_{j,j^{\prime}}\\\forall(j,j^{\prime})\in\mathcal{J}^{2\,(\text{headway})}\cap\mathcal{J}^{2\,(%\text{track})}_{s^{\prime}}\\y^{(out)}(j,j^{\prime},s)=y^{(out)}(j,j^{\prime},s^{\prime})\\\forall(s,s^{\prime})\in S^{\times 2}\forall(j,j^{\prime})\in\mathcal{J}^{2\,(%\text{switch},\text{in},\text{MP})}_{s,s^{\prime}}\cap\mathcal{J}^{2\,(\text{%track})}_{s^{\prime}}\\y^{(in)}(j,j^{\prime},s^{\prime})=y^{(out)}(j,j^{\prime},s^{\prime})\end{split}

(15)

The objective function in Eq.(2), together with theconstraints in Eq.(3)-(14) defineour mathematical programming model.

The model yields an integer linear program, and the time variables can be constrained even into a finite range using the parameter $d_{\max}$ .Note that the binary variables have the obvious symmetry property

\begin{split}y^{(\text{out})}(j^{\prime},j,s)=1-y^{(\text{out})}(j,j^{\prime},%s)\\z(j^{\prime},j,s^{\prime},s)=1-z(j,j^{\prime},s,s^{\prime}),\end{split}

(16)

which is taken into account explicitly by using independent binary variables only.

As for the scaling of our model, the number of time variables i.e. $\#t$ is bounded by number of trains times number of stations $\#\mathcal{J}\#\mathcal{S}$ (see also Section 3.1 of the work of [21]). In order to obtain a more accurate estimate of the number of these variables we observe that in our model no trains visit all stations. We assume that in average each train visits $\alpha\in(0,1)$ fraction of the stations, we have then

\#(t)\approx\alpha\#\mathcal{J}\#\mathcal{S}.

(17)

as a tighter bound.

Suppose, that most of the network is composed mainly of double track lines. Then there will be typically a single precedence variable per train and station ( $y^{(out)}(j,j^{\prime},s)$ ), and thus their number can be estimated by:

\#y+\#z\approx\#y\approx\alpha N(d_{\max})\#\mathcal{J}\#\mathcal{S},

(18)

where $N(d_{\max})$ is number of trains the train can meet at a station on average.Under such assumptions the total number of variables grows linearly with number of trains and stations:

\#vars=\#t+\#y+\#z\approx\text{const}\#\mathcal{J}\#\mathcal{S}.

(19)

If, on the other hand, the network has dominantly single track lines, two precedence variables per train and station ( $y^{(out)}(j,j^{\prime},s)$ and $z(j,j^{\prime},s,s^{\prime})$ ) will be needed, and thus

\#y=\#z\approx\alpha N(d_{\max})\#\mathcal{J}\#\mathcal{S}.

(20)

In both cases we assume that there are not many ( $y^{(in)}(j,j^{\prime},s)$ ) variables, as they are tied certain particular interlocking area conditions.

Based on Eq.(16)there are typically $2$ constraints per each $y$ variable as headway and station track occupation constraints, and $2$ constraints per each $z$ variable as single track occupation constraints. For the interlocking area occupation constraints it is more complicated as it involves both $z$ and $y$ variables; a good assumption is to consider $2$ constraints for each $y$ variable. For the running time and minimal dwell time we expect one constraint per train and station. Thus number of constraints can be estimated by

\#constr.\approx 6\#y+2\#z+2\#\mathcal{J}\#\mathcal{S}

(21)

where the second term is tied to single track line conditions (one condition per $z$ variable).

\#constr.\approx\begin{cases}\alpha(6N(d_{\max})+2)\#\mathcal{J}\#\mathcal{S}%\text{ for double track}\\\alpha(8N(d_{\max})+2)\#\mathcal{J}\#\mathcal{S}\text{ for single track}.\end{cases}

(22)

Because of the limitations in time imposed by $d_{\max}$ the model can be viewed as linear in number of trains and number of stations.Hence, as expected, the model is more complex for the network in which single track lines dominate. Although the analysed network is mostly composed of the double track lines, we will analyse certain use-cases with the fraction of single track lines increased.

Observe that without limiting the secondary delays by $d_{\max}$ , the value of $N(d_{\max})$ would only be bounded by $\#\mathcal{J}-1$ (as each train is considered to possibly meet each other), which would result in quadratic scaling of the number of precedence variables and number of constraints in the number of trains.

5 Hybrid quantum-based approach

To overcome the limitations described in in Section II of currenthardware quantum annealers, we apply hybrid (quantum classical) solver. Inparticular, we use the ’Leap Hybrid Solver Service(HSS)’[11]; a cloud-based proprietary solver which isdeveloped by the market leader of QA hardware and isavailable as a service. It utilizes a hybrid approach that combinesclassical computational power with quantum processing.

In particular, we have used the Constrained Quadratic Model(CQM)[12] Solver. The CQM hybrid solver works as follows.

Input

A constrained linear or quadraticmodel, and the value of the solver parameter t_min, which prescribes the minimum required run time (in seconds) solver is allowed to work on any problem.

Preprocessing

Preprocessing includes subproblem identification and decomposition to smaller instances[52].Furthermore, the solver’s internal preprocessing mechanism takes constraints into accountautomatically, implementing the suitable penalty method internally.

Processing

The hybridsolver implements a workflow that combines a portfolio of variousclassical heuristics (including tabu search, simulated annealing,etc.) running on powerful CPUs (Central Processing Units) and GPUs (Graphics Processing Units). In the course of the solutionprocess the hardware quantum annealer is also invoked, in order tosolve smaller QUBOproblems that can potentially boost the classical heuristics. Integerproblems are transformed to binary using state-of-the-art methods likethe one described by [25][34], generally applicablefor integer problems. When using the quantum hardware, the solverperforms the required embedding, runs the QUBO subproblem on thephysical QA device several times, and returns a sampleof potential solutions. After these readouts, the sample isincorporated into the solution workflow. In this way, even though thephysical annealer supports problems with limited size (e.g. thePegasus machine has $5500$ physical qubits[16],and several of these may be needed to represent a binary variablebecause of the embedding), the approximate solutions of smallsubproblems contribute to the solution process. The final solution is composed by the solver’s module.

Output

The solution of the integer or quadratic constrained problem. The workflow resultsin a solution that is not guaranteed to be optimal, hence, the hybridsolver as a whole is in nature a heuristic solver. Furthermore, the solver is probabilistic, hence each at each run the final solution may be different.

While D-Wave’s proprietary solvers act as black boxes, i.e. they hide the exact details of thedescribed process process from the user, the output is supplemented withtiming parameters, including run time: the total elapsed timeincluding system overhead, and Quantum Processing Unit (QPU) access time which is is thetime spent accessing the actual quantum hardware. These parametersenable a comparison with other solvers.

As mentioned before, in addition to the actual problem to be solved, the CQM solver optionally inputs anotherparameter, t_min.This is a time limit for heuristics running in parallel[12]: unless a thread does not terminate by itself after t_min, it is allowed to run at least till t_min before it is stopped at the current best solution.We have sampled the CQM problems with various settings of t_min anduncovered the impact of this setting to the model performanceand solution quality for ourproblem.

The CQM Solver can solve problems encoded in the form of Eq.(2) – (15) (c.f. Section4.2) with up to $5000$ binary or integer variables and $100,000$ constraints.Computational results of railway rescheduling problems were obtained using classical as well as hybrid (quantum classical) solvers.As a classical solver, IBM ILOG CPLEX Interactive Optimizer (version 22.1.0.0) was used. The CPLEX computation was performed on $16$ cores of Intel(R) Core(TM) i7-10700kF CPU 3.80GHz with $64$ GBs of memory. As hybrid solver, the D-Wave CQM hybrid solver (version: 1.12) was used.

6 Computational results

In this section, we demonstrate the performance of the proposed hybrid (quantum classical) approach on a network at Polish railways - central part of the Upper Silesia Metropolis.In particular, we compare the performance of DWave’s hybrid quantum-classical solver (CQM) with CPLEX.The comparison is performedfor diverse network layouts including single-, double-, and multiple-track lines, and different level of disturbances and disruptions including various initial delays and track closures.As the objective is tied to secondary delays, we can roughly apply it to asses the degree of difficulty of the problem.

Section6.1 describes the considered network and traffic characteristics.Section6.2 presents results concerning generic examples on the core line of the network, comparing double track line scenarios with single track line scenarios and closure scenarios. Section6.3 presents analogous results for larger network with real-life timetable. Here the case studies concerns various scenarios with different degree of difficulty.

6.1 Railway network

We use the Open Railway Map¹¹1https://www.openrailwaymap.org/, visited 2022-11-25 representation of the selected part of the Polish railway network located in the central part of the Metropolis GZM (Poland), presented in Fig2. The selected part of the network combines heterogeneous network layout including single-, double-, triple-, and quadruple-track lines. It comprises 25 nodes including $11$ stations and $3$ junctions, $146$ blocks, and 2 depots.The infrastructure is managed by Polish state infrastructure manager PKP PLK (Polskie Koleje Państwowe, Polskie Linie Kolejowa eng. Polish State Railways, Polish Railway Lines).

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (2)

The decision stations include stations {KO, KO(STM), CB, KL}, stations bounding the analysed part of the network {GLC, CM, KZ, Ty, Mi}, depots {KO(KS), KO(IC)} from which regional and intercity trains are being shunted to and from {KO}, and stations on the single track line – {MJ} to allow meet and pass (M-P) or meet and overtake (M-O).As described before, the decision stations are the only ones that appear explicitly in our model. Due to operational reasons resulting from track design, the decisions on train departure times (and thus on trains’ order) are, in the current operational practice, de facto made with respect to the decision stations.All other stations and junctions are treated as line blocks as long as no rerouting or retracking is considered; they will appear in the model implicitly via parameters. These stations include:

•
stations $\{\text{ZZ},\text{CB}\}$ in which usually there are rigid assignments of the platform tracks to the traffic direction (i.e., track $1$ towards $\{\text{GLC}\}$ ; track $2$ towards $\{\text{KO}\}$ )
•
branch junctions $\{\text{KTC},\text{Bry},\text{Mc}\}$ in which usually there is a rigid assignment to the tracks.

As an example, let for a particular train $j$ , $\hat{\mathcal{S}}_{j}=\{\text{KZ, KO(STM), KO, KTC, CB, RCB, ZZ, GLC}\}$ be the ordered set of stations. Then, the dispatching decisions de facto affect directly the decision stations $\mathcal{S}_{j}=\{\text{KZ, KO(STM), KO, CB, GLC}\}$ . The set of decision stations can be extended if necessary.

In all computation priority weights in Eq.(2) we use $w_{j}=1$ for stopping (local) trains $w_{j}=1.5$ for intercity (fast) trains $w_{j}=1.75$ for express trains, and $w_{j}=0$ for shunting (if applicable).As for the solver parameter t_min, after initial experiments we have decided to prefer t_min $=5$ s as it provides a good balance between computation time and expected solution quality. Sensitivity analysis regarding this parameter is performed in AppendixA for selected instances.

6.2 Synthetic experiments on the selected railway line

As the first set of experiments we address, synthetic instances on a part of the network - railway line KO-GLC, see Fig2.The goal in these scenarios is to compare the CPLEX and CQM solvers for: the double track line ( $line1$ ), the double track line with closures ( $line2$ ), and the single track line ( $line3$ ). Shunting movements and rolling stock circulation are not considered in this set of experiments.In particular, the addressed configurations are the following:

1.
$Line1$ , the double track line with dense traffic. We use cyclic $3$ -hour timetable with $10$ trains each hour and each direction, i.e. we have $59$ trains.
2.
$Line2$ , similar to $line1$ , with the additional closure of one of the tracks between ZZ and RCB. We use a cyclic $2$ -hour timetable with $10$ trains each hour and each direction, i.e. we have $40$ trains.
3.
$Line3$ , the whole line is considered as a single track, with a timetable of $3$ hours traffic and $21$ trains.

For $line1$ and $line3$ , the original timetable is feasible, while for $line2$ it is not feasible due to the closure.For each of the lines we compute $12$ instances, each with different initial delays of trains subsets. Instance $0$ yields no initial delays. Thereafter, we use the parameter value $d_{\max}=40\ \text{min}$ , in all cases in this and subsequent subsection.Assuming that each train interacts with $N(d_{\max})=12$ trains for $line1$ and $line2$ , and with $N(d_{\max})=6$ trains for $line3$ , Eqs.(17) - (19), (22) allows us to quickly estimate the sizes of the problem.This sizes are compared with computed (comp.) actual number of variables and constraints (mean over instances) in Tab.1.Here, $line2$ is the most complex in terms of constraints, what is reflected later by computational difficulties of some of its instances. Furthermore, $line2$ is partially double track and partially single track, hence the actual number of constraints is bounded from above by the estimate for the single track line and from below by the estimate for the double track line. $Line1$ includes double-, whereas $line3$ single-track segments only. There actual number of constraints and variables appeared to be close to the estimate.

\#\mathcal{S}

\#\mathcal{J}

# int vars

# precedence vars

# constraints

comp. mean

Eq.(17)

act. mean

Eq.(18)/Eq.(20)

doub./sing.

comp. mean

Eq.(22)

doub./sing.

118

1538

1415 / -

7889

8731 / -

142

133

1394

1600/3200

8491

6867/13067

592

- / 840

3057

- / 3499

The computational results are presented in Figs.3, 4, and 5 for $line1$ , $line2$ , and $line3$ , respectively. Each figure shows objective $\times$ $d_{\text{max}}$ (top part), computation time (comp. time, middle part) and QPU access time (bottom part) for CQM hybrid solver with t_min equal $5$ and $20$ s. CPLEX results are also displayed.

For scenario $line1$ (Fig.3), CQM solver obtained optimal solutions for all cases (except for instances $6$ and $11$ ).Similarly, for $line3$ (Fig.5), only small deviations from the optimum can be observed for a number of instances.All computational times for $line1$ and $line3$ were significantly below $0.4$ min, although CPLEX performed faster.Comparing CQM with t_min $=5s$ and t_min $=20s$ , we have almost the same objectives but t_min $=20s$ requires almost $4$ times longer computational time.

For $line2$ (Fig.4), all CQM solutions were suboptimal, although feasible. More remarkably, for $9$ instances at t_min $=5s$ and $3$ instances at t_min $=20s$ (out of $12$ instances) CQM was faster than CPLEX, in several cases up to $5\times$ faster.This difference in time is clear while setting the minimal processing time of the CQM solver t_min $=5s$ (see Section5). This observation is approved by more detailed discussion on t_min parameter in AppendixA.In principle, $line1$ (Fig.3) and $line3$ (Fig.5) are less complicated in terms of number of constraints, and here solutions of CPLEX and CQM are similar in terms of objective, but CPLEX is faster. We can conclude that the utility of the CQM solver can be expected for the more complex instances of railway scheduling problems.

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (3)

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (4)

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (5)

6.3 Real-life experiments on the considered rail network

In this subsection, we will test the CQM solver on more realistic scenarios of the railway traffic on presented network.We will also focus on the track closure situation, turning a double track line segment into a single track one, under dense traffic, to elaborate our findings from $line2$ in Section6.2.Our computations address various use cases of train delays outside the network and within the network as well as the partial track closure. We constructed $9$ networks with increased level of difficulty ranging from small delays to disruptions, i.e. it is expected more trains are involved and/or disturbances are spread more over the network. These networks are tabulated in Table 2.Their details are as follows:

1.
Networks $1$ - $3$ concern only disturbances due to delayed trains and no closures.
2.
Networks with number $4$ or higher concern also disturbances due to closures.
1. (a)
  Networks $4$ and $5$ concerns rerouting trains from double track line KTC-CB to the single track line with higher passing times, see Fig.2 (trains have no initial delays in network $3$ and some initial delays in network $4$ ).
2. (b)
  Network $6$ concerns multiple closures, i.e. change of multiple line KZ-KO(STM) and double track line KO-KL to single track lines.(See Fig.2; the reason can be e.g. upcoming reconstruction works on this part of network.)
3. (c)
  Networks $7$ - $9$ concerns closures of both $5$ and $6$ , but each with different initial delays of trains. The intention of the last $3$ networks is to create a really challenging rescheduling problem.
3.
For comparison, network $0$ is the default problem with no disturbances.

network	$\#S$	place	trainsinvolved	description
A. Disturbtions from outside the network
1	10	KO-Ty	IC 14006KS 94113	Delayed IC (higher priority)is in conflict with on-time KS
2	10	arr. Ty	KS 94766KS 40518IC 41004KS 44862IC 4120	Delayed 5 trainson the arrival to Ty, randomdelays of several minutes
B. Disturbtions originated within the network
3	10	wholenetwork	14 trainsKS + IC	Random delays of 14 trainsdelays up to 35 minutes
C. Disturbtions within the network and closures
4	11	KTC-CB	10 trainsKS + IC	Double track line KTC-CB closedtrains rerouted to single track line
5	11	wholenetwork	mosttrains	Double track line KTC-CB closed as in 4+ delays as in 3
6	10	wholenetwork	all trains	Multiple track KZ - KO(STM) and double trackKO - KL changed into single track + delays as in 3
7	11	wholenetwork	all trains	Multiple track KZ - KO(STM) and double trackKO - KL changed into single track + double trackline KTC-CB closed as in 4 + delays as in 3
8	11	wholenetwork	all trains	Closures as in 7 + random delays of 13 trainsup to 30 minutes
9	11	wholenetwork	all trains	Closures as in 7 + random delays of 15 trainsup to 30 minutes

CPLEX

CQM hyb. t_min =

5

mean value over

5

realis.

comparison

CQM vs. CPLEX

#vars /

# constr.

obj.

\times

dmax

comp.

time [s]

obj.

\times

dmax

comp.

time [s]

QPU

time [s]

obj.

diff. %

comp t.

diff %

556 /1756

0.0

0.07

0.0

5.19

0.03

0 %

556 / 1756

1.0

0.08

1.0

5.24

0.02

0 %

556 / 1740

6.0

0.08

6.0

5.37

0.03

0 %

556 / 1769

7.5

0.09

7.5

5.05

0.03

0 %

662 / 2210

78.25

0.20

82.70

5.10

0.03

-5.6 %

662 / 2204

114.75

0.25

132.55

5.25

0.02

-15.5 %

711 / 2599

91.25

0.41

142.3

5.14

0.03

-55.9 %

817 / 3029

188.75

7.98

263.4

5.12

0.02

-39.5 %

35.8 %

817 / 3074

157.75

3.70

271.65

5.12

0.02

-72.2 %

-38.4 %

817 / 3081

185.5

6.51

263.85

5.11

0.02

-42.2 %

21.5 %

Table3 tabulates the results comparing CPLEX and CQM in terms of objective function value and computation time; problem sizes and QPU times are included.As CQM is probabilistic, $5$ runs have been performed for each network, yielding $5$ different realizations. This also facilitates the analysis of the statistics of the output and the differences between particular realizations.

The level of difficulty increases with respect to the number of variables and constraints.Observe that CQM always returns a feasible solution.For the networks $0-3$ with initial delays only, CQM generated optimal solutions, while in the remaining networks, the objective value was somewhat higher than the actual optima obtained with CPLEX.As for computational times, for CPLEX they grow significantly - from $0.072$ to almost $8s$ (network $9$ ). Meanwhile those of CQM remain almost constant around $5.1s$ .QPU times remain not greater than $0.032s$ .Remarkably, for the hardest instances: networks $7$ and $9$ , the CQM hybrid solver even outperforms CPLEX with respect to computational time.It does not find the actual optimum, just a feasible solution close to it.For example, in network $9$ , the computed CQM solution is around 30% worse than the optimal one.In conclusion, we expect benefits from the application of the hybrid solver on medium scale railway networks with multiple closures,as in networks 7-9. This is in accordance with our observation on synthetic data in Fig.4, where the benefits appeared for more complex instances, too.

To evaluate the extent of difficulty of the problems analysed in Tab.3, Tab.4 presents mean (over stations and realizations) delays, resulting from CQM solutions, on selected stations (i.e. stations in central part of the network in Fig.2).Meaningful delays are introduced starting from network $4$ onward, i.e. for networks with closures in a limited number of stations, see also Tab.2.From network $5$ onward, delays become more uniformly distributed among the stations when compared with results in Tab.2.However, for networks $7$ - $9$ the mean delays are the highest and almost uniformly distributed among all stations, showing that the spread of disruption impacts the whole network.These have appeared to be more challenging network disruptions.

Fig.6 shows more detailed statistical analysis for two networks $4$ and $7$ and compares CPLEX (top panel) vs CQM (bottom panel).For network $4$ (left panel of Fig.6), bigger delays occur mainly at stations KO(STM), KO and CB for both CPLEX and CQM solutions, albeit for the former are somewhat smaller.This reflects the initial disturbances, which indeed occurred mainly between KO and CB (see Tab.2), and thus the spread is limited to these 3 stations.For network $7$ (right panel of Fig.6), delays are more spread through the network due to disruption and initial disturbances. To conclude, we can expect advantage from a hybrid (quantum classical) approach particularly in case of more complex problems with extensive disruptions and disturbances occurring throughout the whole network.

	Delay [min]
network	KO(STM)	KO	CB	KL	MJ
0	3.7	1.0	1.7	1.5	0.0
1	3.1	1.3	2.1	1.5	0.0
2	3.8	1.3	2.0	1.3	2.0
3	2.6	2.1	1.4	0.7	1.5
4	6.6	4.4	7.4	1.6	0.0
5	7.1	6.2	10.7	1.0	2.7
6	6.7	6.7	3.2	8.7	4.1
7	9.8	10.7	12.5	9.5	9.9
8	10.4	10.5	10.0	10.5	7.6
9	10.0	10.3	11.5	9.8	9.6

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (6)

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (7)

In order to further evaluate the solutions, we present time-distance diagrams for network $7$ , for the CPLEX solution and two different quantum realizations; these can be seen in Figs.7, 8, and 9, respectively. We have chosen the line segment between GLC and KZ, offering a good demonstration.(We have introduced a small artificial ’distance’ between KO and KO(STM) wich are contiguous locations for sake of better visibility, hence the horizontal lines in the diagram.)The red dotted lines represent the conflicted situation, and the green solid lines represent the solution. All the three proposed solutions have similar figures of merit, they all provide a valid option to be realized.The conflicts are resolved in somewhat different ways, see e.g. the paths of trains 40150, 5312, and 26103.In details:

•
in the CPLEX solution, train 26103 passes first the KTC-CB closure heading for GLC, then 40150 followed by 5312 go in opposite directions,
•
in the CQM solution, first realization, 40150 and 5312 pass first the KTC-CB closure and 26103 waits at KTC for them to pass,
•
in the CQM solution, second realization, 26103 passes first KTC-CB closure, then 5312 followed by 40150 go in opposite directions.

This results in following objective $\times d_{\max}$ values: for CPLEX - $188.75$ , for CQM first realization - $279$ for CQM second realization - $261.5$ . Although, the first solution is optimal, it is reasonable to supply dispatchers with the range feasible solutions.Hence, these diagrams proved that CQM hybrid solver can produce usable solutions for railway traffic management.

An additional benefit of the probabilistic nature of the quantum-based method is that it produces different alternatives after each realization run without any modification. These could be offered to the dispatcher in a decision support system as possibilities. The dispatcher’s choice can be influenced by other circ*mstances not covered by the present model.

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (10)

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (11)

Solving rescheduling problems in heterogeneous urban railway networks using hybrid quantum-classical approach (12)

7 Conclusions

We have demonstrated that quantum annealers can be readily applicable in train rescheduling optimization in urban railway networks. In particular we have encoded rescheduling problems in a linear integer program, and solved these using a state-of-the-art solver (CPLEX), as well as on the QA based CQM hybrid (quantum classical) solver.While the CQM hybrid solver does not outperform classical solvers in general, we have found specific cases in which they were actually better. This supports the expectation that future quantum devices will be efficient in solving large scale-problems, and in our case, real-time railway dispatching problems on the national scale; even in the range that are beyond the scope of current exact models and heuristics.

While the quantum-based solvers would possibly return suboptimal solutions, these are still often better than the solutions obtained manually or based on smaller scale models. More importantly, all generated solutions were feasible. In each computation, the CQM hybrid solver always has taken some advantage of the QPU.Hence, we believe that such a little help from the QPU must have boosted the classical heuristics in the hybrid solver.

Recall that quantum annealer devices are subjects of many ongoing debates. Previous studies claimed that the quantum nature of their operation is limited[45], while recent works argue that their operation relies more on the thermalization with a cold environment in place of actual quantum adiabatic evolution[6]. They are often criticized for the small size of problems they can address and various classical solvers do outperform them on certain problems at the present state of the art. In spite of all these, our computational results demonstrate the readiness of hybrid (quantum classical) solvers to handle real-life railway problems. We believe that future investigations shall focus on hybrid approaches both for railway applications as well as more generally.

Finally, several future research directions can be determined. First concerns the ongoing worldwide works on the development of fault tolerant quantum computing devices [17] with proper error correction/mitigation, this willboost results of approaches supported by quantum computing (and QA in particular). Second concerns creation of custom open source hybrid solver that is dedicated for railway problems and will use noise of quantum device to our advantage e.g. by modelling some stochastic behaviour on railways. Third concerns evolution and application of quantum inspired techniques such as tensor networks technique, simulated bifurcation[26].With all these future developments, QA-based approaches tend to become more powerful and ready to address for the first time large-scale real-time challenges in railway traffic management, which are currently unsolvable with traditional techniques.

Data availability

The code and the data used for generating the numerical results can be found in https://github.com/iitis/railways_dispatching_silesia under DOI 10.5281/zenodo.10657323

Acknowledgements

This research was supported by the Ministry of Culture and Innovation and the National Research, Development and Innovation Office within the Quantum Information National Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004) (MK), and by the Silesian University of Technology Rector’s Grant no. BKM - 720/RT2/2023 12/020/BKM_2023/0252 (KK).This research was funded in part by National Science Center, Poland, under grant number 2019/33/B/ST6/0201 (LB) and 2023/07/X/ST6/00396 (KD). For the purpose of Open Access, the authors has applied a CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.

We would like to thank to Sebastian Deffner as well as to the whole scientific community of the Department of Physics, University of Maryland, Baltimore County (UMBC) for valuable discussion; and to Bartłomiej Gardas, and Zbigniew Puchała for valuable motivation; and to Katarzyna Gawlak, Akash Kundu, and Özlem Salehi for data validation and assistance with coding.We acknowledge the cooperation with Koleje Ślaskie sp. z o.o. (eng. Silesian Railways).

Compliance with ethical standards

Conflict of interest

The authors declare that they have no conflict of interest.

Ethical approval

This article does not contain any studies with human participants or animals performed byany of the authors.

References

\bibcommenthead
Albash\BOthers. [\APACyear2017]\APACinsertmetastaralbash2017temperature{APACrefauthors}Albash, T., Martin-Mayor, V.\BCBL Hen, I.\APACrefYearMonthDay2017.\BBOQ\APACrefatitleTemperature scaling law for quantum annealingoptimizers Temperature scaling law for quantum annealingoptimizers.\BBCQ\APACjournalVolNumPagesPhysical Review Letters11911110502,\PrintBackRefs\CurrentBib
Apolloni\BOthers. [\APACyear1989]\APACinsertmetastarapolloni1989quantum{APACrefauthors}Apolloni, B., Carvalho, C.\BCBL DeFalco, D.\APACrefYearMonthDay1989.\BBOQ\APACrefatitleQuantum stochastic optimization Quantum stochasticoptimization.\BBCQ\APACjournalVolNumPagesStochastic Processes and theirApplications332233–244,\PrintBackRefs\CurrentBib
Avron\BBA Elgart [\APACyear1999]\APACinsertmetastaravron.elgart.99{APACrefauthors}Avron, J.E.\BCBT\BBA Elgart, A.\APACrefYearMonthDay199901.\BBOQ\APACrefatitleAdiabatic Theorem without a Gap Condition Adiabatictheorem without a gap condition.\BBCQ\APACjournalVolNumPagesCommun. Math. Phys.2032445–463,\PrintBackRefs\CurrentBib
Bian\BOthers. [\APACyear2010]\APACinsertmetastarBian2010TheIM{APACrefauthors}Bian, Z., Chudak, F., Macready, W.G.\BCBL Rose, G.\APACrefYearMonthDay2010.\APACrefbtitleThe Ising model: teaching an old problem new tricks. TheIsing model: teaching an old problem new tricks.{APACrefURL}https://www.dwavesys.com/media/vbklsvbh/weightedmaxsat_v2.pdf\APACrefnotevisited 2022.04.29.\PrintBackRefs\CurrentBib
Boros\BBA Hammer [\APACyear1991]\APACinsertmetastarboros1991max{APACrefauthors}Boros, E.\BCBT\BBA Hammer, P.L.\APACrefYearMonthDay1991.\BBOQ\APACrefatitleThe Max-Cut problem and quadratic 0–1 optimization;polyhedral aspects, relaxations and bounds The Max-Cut problem andquadratic 0–1 optimization; polyhedral aspects, relaxations andbounds.\BBCQ\APACjournalVolNumPagesAnnals of Operations Research333151–180,\PrintBackRefs\CurrentBib
Buffoni\BBA Campisi [\APACyear2020]\APACinsertmetastarBuffoni_2020{APACrefauthors}Buffoni, L.\BCBT\BBA Campisi, M.\APACrefYearMonthDay2020jun.\BBOQ\APACrefatitleThermodynamics of a quantum annealer Thermodynamics ofa quantum annealer.\BBCQ\APACjournalVolNumPagesQuantum Sci. Technol.53035013,\PrintBackRefs\CurrentBib
Cacchiani\BOthers. [\APACyear2014]\APACinsertmetastarcacchiani2014overview{APACrefauthors}Cacchiani, V., Huisman, D., Kidd, M., Kroon, L., Toth, P., Veelenturf, L.\BCBL Wagenaar, J.\APACrefYearMonthDay2014.\BBOQ\APACrefatitleAn overview of recovery models and algorithms forreal-time railway rescheduling An overview of recovery models andalgorithms for real-time railway rescheduling.\BBCQ\APACjournalVolNumPagesTransportation Research Part B:Methodological6315–37,\PrintBackRefs\CurrentBib
Caimi\BOthers. [\APACyear2012]\APACinsertmetastarCAIMI20122578{APACrefauthors}Caimi, G., Fuchsberger, M., Laumanns, M.\BCBL Lüthi, M.\APACrefYearMonthDay2012.\BBOQ\APACrefatitleA model predictive control approach for discrete-timerescheduling in complex central railway station areas A model predictivecontrol approach for discrete-time rescheduling in complex central railwaystation areas.\BBCQ\APACjournalVolNumPagesComputers & Operations Research39112578-2593,\PrintBackRefs\CurrentBib
Caprara\BOthers. [\APACyear2002]\APACinsertmetastarcaprara2002modeling{APACrefauthors}Caprara, A., Fischetti, M.\BCBL Toth, P.\APACrefYearMonthDay2002.\BBOQ\APACrefatitleModeling and solving the train timetabling problemModeling and solving the train timetabling problem.\BBCQ\APACjournalVolNumPagesOperations Research505851–861,\PrintBackRefs\CurrentBib
Corman\BOthers. [\APACyear2012]\APACinsertmetastarCORMAN201279{APACrefauthors}Corman, F., D’Ariano, A., Pacciarelli, D.\BCBL Pranzo, M.\APACrefYearMonthDay2012.\BBOQ\APACrefatitleBi-objective conflict detection and resolution inrailway traffic management Bi-objective conflict detection and resolutionin railway traffic management.\BBCQ\APACjournalVolNumPagesTransportation Research Part C: EmergingTechnologies20179-94,\APACrefnoteSpecial issue on Optimization in Public Transport+ISTT2011\PrintBackRefs\CurrentBib
D-Wave Quantum Inc. [\APACyear2022\APACexlab\BCnt1]\APACinsertmetastarDwaveHSS{APACrefauthors}D-Wave Quantum Inc.\APACrefYearMonthDay2022\BCnt1.\APACrefbtitleD-Wave Hybrid Solver Service: An Overview [WhitePaper].D-Wave Hybrid Solver Service: An Overview [WhitePaper].{APACrefURL}https://www.dwavesys.com/media/4bnpi53x/14-1039a-b_d-wave_hybrid_solver_service_an_overview.pdf\APACrefnotevisited 2022.04.29.\PrintBackRefs\CurrentBib
D-Wave Quantum Inc. [\APACyear2022\APACexlab\BCnt2]\APACinsertmetastardwaveCQM{APACrefauthors}D-Wave Quantum Inc.\APACrefYearMonthDay2022\BCnt2.\APACrefbtitleHybrid Solver for Constrained Quadratic Model [WhitePaper].Hybrid Solver for Constrained Quadratic Model [WhitePaper].{APACrefURL}https://www.dwavesys.com/media/rldh2ghw/14-1055a-a_hybrid_solver_for_constrained_quadratic_models.pdf\APACrefnotevisited 2022.04.29.\PrintBackRefs\CurrentBib
DalSasso\BOthers. [\APACyear2021]\APACinsertmetastardal2021tick{APACrefauthors}DalSasso, V., Lamorgese, L., Mannino, C., Onofri, A.\BCBL Ventura, P.\APACrefYearMonthDay2021.\BBOQ\APACrefatitleThe tick formulation for deadlock detection andavoidance in railways traffic control The tick formulation for deadlockdetection and avoidance in railways traffic control.\BBCQ\APACjournalVolNumPagesJournal of Rail Transport Planning &Management17100239,\PrintBackRefs\CurrentBib
D’Ariano\BOthers. [\APACyear2007]\APACinsertmetastarDarianoPDPbb{APACrefauthors}D’Ariano, A., Pacciarelli, D.\BCBL Pranzo, M.\APACrefYearMonthDay2007.\BBOQ\APACrefatitleA branch and bound algorithm for scheduling trains in arailway network A branch and bound algorithm for scheduling trains in arailway network.\BBCQ\APACjournalVolNumPagesEuropean Journal of OperationalResearch1832643-657,\PrintBackRefs\CurrentBib
Das\BBA Chakrabarti [\APACyear2008]\APACinsertmetastardas_colloquium_2008{APACrefauthors}Das, A.\BCBT\BBA Chakrabarti, B.K.\APACrefYearMonthDay2008\APACmonth09.\BBOQ\APACrefatitleColloquium : Quantum annealing and analogquantum computation Colloquium : Quantum annealing and analogquantum computation.\BBCQ\APACjournalVolNumPagesReviews of Modern Physics8031061–1081,\PrintBackRefs\CurrentBib
Dattani\BOthers. [\APACyear2019]\APACinsertmetastardattani2019pegasus{APACrefauthors}Dattani, N., Szalay, S.\BCBL Chancellor, N.\APACrefYearMonthDay2019Jan.\BBOQ\APACrefatitlePegasus: The second connectivity graph for large-scalequantum annealing hardware Pegasus: The second connectivity graph forlarge-scale quantum annealing hardware.\BBCQarXiv:1901.07636[quant-ph]\PrintBackRefs\CurrentBib
DiVincenzo [\APACyear2000]\APACinsertmetastardivincenzo2000physical{APACrefauthors}DiVincenzo, D.P.\APACrefYearMonthDay2000.\BBOQ\APACrefatitleThe physical implementation of quantum computation Thephysical implementation of quantum computation.\BBCQ\APACjournalVolNumPagesFortschritte der Physik: Progress ofPhysics489-11771–783,\PrintBackRefs\CurrentBib
Domino\BOthers. [\APACyear2023]\APACinsertmetastardomino2023quantum{APACrefauthors}Domino, K., Koniorczyk, M., Krawiec, K., Jałowiecki, K., Deffner, S.\BCBL Gardas, B.\APACrefYearMonthDay2023.\BBOQ\APACrefatitleQuantum annealing in the NISQ era: railway conflictmanagement Quantum annealing in the NISQ era: railway conflictmanagement.\BBCQ\APACjournalVolNumPagesEntropy252191,\PrintBackRefs\CurrentBib
Domino\BOthers. [\APACyear2021]\APACinsertmetastardomino2021quantum{APACrefauthors}Domino, K., Koniorczyk, M., Krawiec, K., Jałowiecki, K.\BCBL Gardas, B.\APACrefYearMonthDay2021.\BBOQ\APACrefatitleQuantum computing approach to railway dispatching andconflict management optimization on single-track railway lines Quantumcomputing approach to railway dispatching and conflict managementoptimization on single-track railway lines.\BBCQarXiv:2010.08227[cs.ET]\PrintBackRefs\CurrentBib
Domino, Koniorczyk\BCBL\BBA Puchała [\APACyear2022]\APACinsertmetastardomino2022statistical{APACrefauthors}Domino, K., Koniorczyk, M.\BCBL Puchała, Z.\APACrefYearMonthDay2022.\BBOQ\APACrefatitleStatistical quality assessment of Ising-based annealeroutputs Statistical quality assessment of ising-based annealeroutputs.\BBCQ\APACjournalVolNumPagesQuantum Information Processing218288,\PrintBackRefs\CurrentBib
Domino, Kundu\BCBL\BOthers. [\APACyear2022]\APACinsertmetastarqubohobo{APACrefauthors}Domino, K., Kundu, A., Salehi, Ö.\BCBL Krawiec, K.\APACrefYearMonthDay2022.\BBOQ\APACrefatitleQuadratic and higher-order unconstrained binaryoptimization of railway rescheduling for quantum computing Quadratic andhigher-order unconstrained binary optimization of railway rescheduling forquantum computing.\BBCQ\APACjournalVolNumPagesQuantum Information Processing2191–33,\PrintBackRefs\CurrentBib
Dunning\BOthers. [\APACyear2018]\APACinsertmetastarDunning_2018{APACrefauthors}Dunning, I., Gupta, S.\BCBL Silberholz, J.\APACrefYearMonthDay2018aug.\BBOQ\APACrefatitleWhat Works Best When? A Systematic Evaluation ofHeuristics for Max-Cut and QUBO What Works Best When? A SystematicEvaluation of Heuristics for Max-Cut and QUBO.\BBCQ\APACjournalVolNumPagesINFORMS Journal on Computing303608–624,\PrintBackRefs\CurrentBib
Ge\BOthers. [\APACyear2022]\APACinsertmetastarGe2022{APACrefauthors}Ge, L., Voß, S.\BCBL Xie, L.\APACrefYearMonthDay2022Jun04.\BBOQ\APACrefatitleRobustness and disturbances in public transportRobustness and disturbances in public transport.\BBCQ\APACjournalVolNumPagesPublic Transport14191-261,\PrintBackRefs\CurrentBib
Glos\BOthers. [\APACyear2023]\APACinsertmetastarglos2022optimizing{APACrefauthors}Glos, A., Kundu, A.\BCBL Salehi, Ö.\APACrefYearMonthDay2023.\BBOQ\APACrefatitleOptimizing the production of test vehicles using hybridconstrained quantum annealing Optimizing the production of test vehiclesusing hybrid constrained quantum annealing.\BBCQ\APACjournalVolNumPagesSN Computer Science45609,\PrintBackRefs\CurrentBib
Glover\BOthers. [\APACyear2019]\APACinsertmetastarglover2019quantum{APACrefauthors}Glover, F., Kochenberger, G.\BCBL Du, Y.\APACrefYearMonthDay2019.\BBOQ\APACrefatitleQuantum Bridge Analytics I: a tutorial on formulatingand using QUBO models Quantum bridge analytics i: a tutorial on formulatingand using qubo models.\BBCQ\APACjournalVolNumPages4or17335–371,\PrintBackRefs\CurrentBib
Goto [\APACyear2019]\APACinsertmetastargoto2019quantum{APACrefauthors}Goto, H.\APACrefYearMonthDay2019.\BBOQ\APACrefatitleQuantum computation based on quantum adiabaticbifurcations of Kerr-nonlinear parametric oscillators Quantum computationbased on quantum adiabatic bifurcations of Kerr-nonlinear parametricoscillators.\BBCQ\APACjournalVolNumPagesJournal of the Physical Society ofJapan886061015,\PrintBackRefs\CurrentBib
Grozea\BOthers. [\APACyear2021]\APACinsertmetastar2109.07212v1{APACrefauthors}Grozea, C., Hans, R., Koch, M., Riehn, C.\BCBL Wolf, A.\APACrefYearMonthDay2021Sep.\BBOQ\APACrefatitleOptimising Rolling Stock Planning including Maintenancewith Constraint Programming and Quantum Annealing Optimising rolling stockplanning including maintenance with constraint programming and quantumannealing.\BBCQarXiv:2109.07212v1[cs.AI]\PrintBackRefs\CurrentBib
Gusmeroli\BOthers. [\APACyear2022]\APACinsertmetastarGusmeroli_2022{APACrefauthors}Gusmeroli, N., Hrga, T., Lužar, B., Povh, J., Siebenhofer, M.\BCBL Wiegele, A.\APACrefYearMonthDay2022jun.\BBOQ\APACrefatitleBiqBin: A Parallel Branch-and-bound Solver for BinaryQuadratic Problems with Linear Constraints BiqBin: A parallelbranch-and-bound solver for binary quadratic problems with linearconstraints.\BBCQ\APACjournalVolNumPagesACM Trans. Math. Softw.4821–31,\PrintBackRefs\CurrentBib
Harrod [\APACyear2011]\APACinsertmetastarharrod2011modeling{APACrefauthors}Harrod, S.\APACrefYearMonthDay2011.\BBOQ\APACrefatitleModeling network transition constraints withhypergraphs Modeling network transition constraints withhypergraphs.\BBCQ\APACjournalVolNumPagesTransportation Science45181–97,\PrintBackRefs\CurrentBib
Ising [\APACyear1925]\APACinsertmetastarIsing_1925{APACrefauthors}Ising, E.\APACrefYearMonthDay1925feb.\BBOQ\APACrefatitleBeitrag zur Theorie des Ferromagnetismus Beitrag zurtheorie des ferromagnetismus.\BBCQ\APACjournalVolNumPagesZ. Physik311253–258,\PrintBackRefs\CurrentBib
Johnson\BOthers. [\APACyear2011]\APACinsertmetastarjohnson2011quantum{APACrefauthors}Johnson, M.W., Amin, M.H., Gildert, S., Lanting, T., Hamze, F., Dickson, N.\BDBLothers\APACrefYearMonthDay2011.\BBOQ\APACrefatitleQuantum annealing with manufactured spins Quantumannealing with manufactured spins.\BBCQ\APACjournalVolNumPagesNature4737346194–198,\PrintBackRefs\CurrentBib
Jünger\BOthers. [\APACyear2021]\APACinsertmetastarJ_nger_2021{APACrefauthors}Jünger, M., Lobe, E., Mutzel, P., Reinelt, G., Rendl, F., Rinaldi, G.\BCBL Stollenwerk, T.\APACrefYearMonthDay2021dec.\BBOQ\APACrefatitleQuantum Annealing versus Digital Computing Quantumannealing versus digital computing.\BBCQ\APACjournalVolNumPagesACM Journal of Experimental Algorithmics261–30,\PrintBackRefs\CurrentBib
Kadowaki\BBA Nishimori [\APACyear1998]\APACinsertmetastarkadowaki1998quantum{APACrefauthors}Kadowaki, T.\BCBT\BBA Nishimori, H.\APACrefYearMonthDay1998.\BBOQ\APACrefatitleQuantum annealing in the transverse Ising modelQuantum annealing in the transverse Ising model.\BBCQ\APACjournalVolNumPagesPhysical Review E5855355,\PrintBackRefs\CurrentBib
Karimi\BBA Ronagh [\APACyear2019]\APACinsertmetastarkarimi2019practical{APACrefauthors}Karimi, S.\BCBT\BBA Ronagh, P.\APACrefYearMonthDay2019.\BBOQ\APACrefatitlePractical integer-to-binary mapping for quantumannealers Practical integer-to-binary mapping for quantumannealers.\BBCQ\APACjournalVolNumPagesQuantum Information Processing1841–24,\PrintBackRefs\CurrentBib
Lamorgese\BOthers. [\APACyear2018]\APACinsertmetastarLamorgese2018{APACrefauthors}Lamorgese, L., Mannino, C., Pacciarelli, D.\BCBL Krasemann, J.T.\APACrefYearMonthDay2018.\BBOQ\APACrefatitleTrain Dispatching Train dispatching.\BBCQ\BIn R.Borndörfer, T.Klug, L.Lamorgese, C.Mannino, M.Reuther\BCBL\BBA T.Schlechte(\BEDS), \APACrefbtitleHandbook of Optimizationin the Railway Industry Handbook of Optimization in the RailwayIndustry(\BPGS265–283).\APACaddressPublisherChamSpringer International Publishing.\PrintBackRefs\CurrentBib
Lange\BBA Werner [\APACyear2018]\APACinsertmetastarlange2018approaches{APACrefauthors}Lange, J.\BCBT\BBA Werner, F.\APACrefYearMonthDay2018.\BBOQ\APACrefatitleApproaches to modeling train scheduling problems asjob-shop problems with blocking constraints Approaches to modeling trainscheduling problems as job-shop problems with blocking constraints.\BBCQ\APACjournalVolNumPagesJournal of Scheduling212191–207,\PrintBackRefs\CurrentBib
Lucas [\APACyear2014]\APACinsertmetastarlucas2014ising{APACrefauthors}Lucas, A.\APACrefYearMonthDay2014.\BBOQ\APACrefatitleIsing formulations of many NP problems Isingformulations of many NP problems.\BBCQ\APACjournalVolNumPagesFrontiers in Physics25,\PrintBackRefs\CurrentBib
Lusby\BOthers. [\APACyear2018]\APACinsertmetastarlusby2018survey{APACrefauthors}Lusby, R.M., Larsen, J.\BCBL Bull, S.\APACrefYearMonthDay2018.\BBOQ\APACrefatitleA survey on robustness in railway planning A survey onrobustness in railway planning.\BBCQ\APACjournalVolNumPagesEuropean Journal of OperationalResearch26611–15,\PrintBackRefs\CurrentBib
Lusby\BOthers. [\APACyear2013]\APACinsertmetastarLUSBY2013713{APACrefauthors}Lusby, R.M., Larsen, J., Ehrgott, M.\BCBL Ryan, D.M.\APACrefYearMonthDay2013.\BBOQ\APACrefatitleA set packing inspired method for real-time junctiontrain routing A set packing inspired method for real-time junction trainrouting.\BBCQ\APACjournalVolNumPagesComputers & Operations Research403713-724,\PrintBackRefs\CurrentBib
Martins\BOthers. [\APACyear2021]\APACinsertmetastarmartins2021qubo{APACrefauthors}Martins, L.N., Rocha, A.P.\BCBL Castro, A.J.\APACrefYearMonthDay2021.\BBOQ\APACrefatitleA QUBO Model to the Tail Assignment Problem. A QUBOModel to the Tail Assignment Problem.\BBCQ\APACrefbtitleICAART (2) Icaart (2)(\BPGS899–906).\PrintBackRefs\CurrentBib
Mascis\BBA Pacciarelli [\APACyear2002]\APACinsertmetastarMASCIS2002498{APACrefauthors}Mascis, A.\BCBT\BBA Pacciarelli, D.\APACrefYearMonthDay2002.\BBOQ\APACrefatitleJob-shop scheduling with blocking and no-waitconstraints Job-shop scheduling with blocking and no-waitconstraints.\BBCQ\APACjournalVolNumPagesEuropean Journal of OperationalResearch1433498-517,\PrintBackRefs\CurrentBib
McLeod\BBA Sasdelli [\APACyear2022]\APACinsertmetastarMcLeod_2022{APACrefauthors}McLeod, C.R.\BCBT\BBA Sasdelli, M.\APACrefYearMonthDay2022.\BBOQ\APACrefatitleBenchmarking D-Wave Quantum Annealers: Spectral GapScaling ofMaximum Cardinality Matching Problems Benchmarking D-WaveQuantum Annealers: Spectral Gap Scaling ofMaximum Cardinality MatchingProblems.\BBCQ\APACrefbtitleComputational Science – ICCS 2022Computational science – ICCS 2022(\BPGS150–163).\APACaddressPublisherSpringer International Publishing.\PrintBackRefs\CurrentBib
Meng\BBA Zhou [\APACyear2014]\APACinsertmetastarMENG2014208{APACrefauthors}Meng, L.\BCBT\BBA Zhou, X.\APACrefYearMonthDay2014.\BBOQ\APACrefatitleSimultaneous train rerouting and rescheduling on anN-track network: A model reformulation with network-based cumulative flowvariables Simultaneous train rerouting and rescheduling on an n-tracknetwork: A model reformulation with network-based cumulative flowvariables.\BBCQ\APACjournalVolNumPagesTransportation Research Part B:Methodological67208-234,\PrintBackRefs\CurrentBib
Neukart\BOthers. [\APACyear2017]\APACinsertmetastarneukart2017traffic{APACrefauthors}Neukart, F., Compostella, G., Seidel, C., von Dollen, D., Yarkoni, S.\BCBL Parney, B.\APACrefYearMonthDay2017\APACmonth12.\BBOQ\APACrefatitleTraffic Flow Optimization Using a Quantum AnnealerTraffic flow optimization using a quantum annealer.\BBCQ\APACjournalVolNumPagesFrontiers in ICT4,{APACrefDOI} https://doi.org/10.3389/fict.2017.00029{APACrefURL} http://dx.doi.org/10.3389/fict.2017.00029 \PrintBackRefs\CurrentBib
Shin\BOthers. [\APACyear2014]\APACinsertmetastar1401.7087v2{APACrefauthors}Shin, S.W., Smith, G., Smolin, J.A.\BCBL Vazirani, U.\APACrefYearMonthDay2014Jan.\BBOQ\APACrefatitleHow "Quantum" is the D-Wave Machine? How "quantum"is the D-Wave machine?\BBCQarXiv:1401.7087v2[quant-ph]\PrintBackRefs\CurrentBib
Soriani\BOthers. [\APACyear2022]\APACinsertmetastarsoriani2022three{APACrefauthors}Soriani, A., Nazé, P., Bonança, M.V., Gardas, B.\BCBL Deffner, S.\APACrefYearMonthDay2022.\BBOQ\APACrefatitleThree phases of quantum annealing: Fast, slow, and veryslow Three phases of quantum annealing: Fast, slow, and very slow.\BBCQ\APACjournalVolNumPagesPhysical Review A1054042423,\PrintBackRefs\CurrentBib
Stollenwerk\BOthers. [\APACyear2019]\APACinsertmetastar10.1007/978-3-030-14082-3_9{APACrefauthors}Stollenwerk, T., Lobe, E.\BCBL Jung, M.\APACrefYearMonthDay2019.\BBOQ\APACrefatitleFlight Gate Assignment with a Quantum AnnealerFlight Gate Assignment with a Quantum Annealer.\BBCQS.Feld\BBA C.Linnhoff-Popien(\BEDS), \APACrefbtitleQuantumTechnology and Optimization Problems Quantum technology and optimizationproblems(\BPGS99–110).\APACaddressPublisherChamSpringer International Publishing.\PrintBackRefs\CurrentBib
Stollenwerk\BOthers. [\APACyear2020]\APACinsertmetastar8643733{APACrefauthors}Stollenwerk, T., O’Gorman, B., Venturelli, D., Mandrà, S., Rodionova, O., Ng, H.\BDBLBiswas, R.\APACrefYearMonthDay2020.\BBOQ\APACrefatitleQuantum Annealing Applied to De-Conflicting OptimalTrajectories for Air Traffic Management Quantum annealing applied tode-conflicting optimal trajectories for air traffic management.\BBCQ\APACjournalVolNumPagesIEEE Transactions on Intelligent TransportationSystems211285-297,\PrintBackRefs\CurrentBib
Szpigel [\APACyear1973]\APACinsertmetastarSzpigel1973{APACrefauthors}Szpigel, B.\APACrefYearMonthDay1973.\BBOQ\APACrefatitleOptimal train scheduling on a single line railwayOptimal train scheduling on a single line railway.\BBCQ\APACjournalVolNumPagesOperational Research72343-352,\PrintBackRefs\CurrentBib
Törnquist [\APACyear2007]\APACinsertmetastartornquist2007railway{APACrefauthors}Törnquist, J.\APACrefYearMonthDay2007.\BBOQ\APACrefatitleRailway traffic disturbance management—Anexperimental analysis of disturbance complexity, management objectives andlimitations in planning horizon Railway traffic disturbancemanagement—An experimental analysis of disturbance complexity, managementobjectives and limitations in planning horizon.\BBCQ\APACjournalVolNumPagesTransportation Research Part A: Policy andPractice413249–266,\PrintBackRefs\CurrentBib
P.N.Tran\BOthers. [\APACyear2020]\APACinsertmetastarTran_2020{APACrefauthors}Tran, P.N., Pham, D\BHBIT., Goh, S.K., Alam, S.\BCBL Duong, V.\APACrefYearMonthDay2020jun.\BBOQ\APACrefatitleAn Interactive Conflict Solver for Learning Air TrafficConflict Resolutions An Interactive Conflict Solver for Learning AirTraffic Conflict Resolutions.\BBCQ\APACjournalVolNumPagesJournal of Aerospace InformationSystems176271–277,\PrintBackRefs\CurrentBib
T.Tran\BOthers. [\APACyear2016]\APACinsertmetastartran2016hybrid{APACrefauthors}Tran, T., Do, M., Rieffel, E., Frank, J., Wang, Z., O’Gorman, B.\BDBLBeck, J.\APACrefYearMonthDay2016.\BBOQ\APACrefatitleA hybrid quantum-classical approach to solvingscheduling problems A hybrid quantum-classical approach to solvingscheduling problems.\BBCQ\APACrefbtitleProceedings of the International Symposium on CombinatorialSearch Proceedings of the international symposium on combinatorial search(\BVOL7, \BPGS98–106).\PrintBackRefs\CurrentBib
Venturelli\BOthers. [\APACyear2015]\APACinsertmetastarventurelli2016job{APACrefauthors}Venturelli, D., Marchand, D.J.J.\BCBL Rojo, G.\APACrefYearMonthDay2015.\BBOQ\APACrefatitleQuantum Annealing Implementation of Job-Shop SchedulingQuantum annealing implementation of job-shop scheduling.\BBCQarXiv:1506.08479\PrintBackRefs\CurrentBib
Venuti\BOthers. [\APACyear2016]\APACinsertmetastarPhysRevA.93.032118{APACrefauthors}Venuti, L.C., Albash, T., Lidar, D.A.\BCBL Zanardi, P.\APACrefYearMonthDay2016Mar.\BBOQ\APACrefatitleAdiabaticity in open quantum systems Adiabaticity inopen quantum systems.\BBCQ\APACjournalVolNumPagesPhys. Rev. A93032118,\PrintBackRefs\CurrentBib
Willsch\BOthers. [\APACyear2022]\APACinsertmetastarWillsch_2022{APACrefauthors}Willsch, D., Willsch, M., Calaza, C.D.G., Jin, F., Raedt, H.D., Svensson, M.\BCBL Michielsen, K.\APACrefYearMonthDay2022apr.\BBOQ\APACrefatitleBenchmarking Advantage and D-Wave 2000Q quantumannealers with exact cover problems Benchmarking advantage and D-Wave2000q quantum annealers with exact cover problems.\BBCQ\APACjournalVolNumPagesQuantum Inf Process214,\PrintBackRefs\CurrentBib
Xu\BOthers. [\APACyear2023]\APACinsertmetastarxu2023high{APACrefauthors}Xu, H\BHBIZ., Chen, J\BHBIH., Zhang, X\BHBIC., Lu, T\BHBIE., Gao, T\BHBIZ., Wen, K.\BCBL Ma, Y.\APACrefYearMonthDay2023.\BBOQ\APACrefatitleHigh-speed train timetable optimization based onspace–time network model and quantum simulator High-speed train timetableoptimization based on space–time network model and quantumsimulator.\BBCQ\APACjournalVolNumPagesQuantum Information Processing2211418,\PrintBackRefs\CurrentBib
Yarkoni, Alekseyenko\BCBL\BOthers. [\APACyear2021]\APACinsertmetastaryarkoni2021multi{APACrefauthors}Yarkoni, S., Alekseyenko, A., Streif, M., VonDollen, D., Neukart, F.\BCBL Bäck, T.\APACrefYearMonthDay2021.\BBOQ\APACrefatitleMulti-car paint shop optimization with quantumannealing Multi-car paint shop optimization with quantum annealing.\BBCQ\APACrefbtitle2021 IEEE International Conference on Quantum Computing andEngineering (QCE) 2021 ieee international conference on quantum computingand engineering (qce)(\BPGS35–41).\PrintBackRefs\CurrentBib
Yarkoni, Huck\BCBL\BOthers. [\APACyear2021]\APACinsertmetastar10.1007/978-3-030-87672-2_33{APACrefauthors}Yarkoni, S., Huck, A., Schülldorf, H., Speitkamp, B., Tabrizi, M.S., Leib, M.\BDBLNeukart, F.\APACrefYearMonthDay2021.\BBOQ\APACrefatitleSolving the Shipment Rerouting Problem with QuantumOptimization Techniques Solving the shipment rerouting problem with quantumoptimization techniques.\BBCQM.Mes, E.Lalla-Ruiz\BCBL\BBA S.Voß(\BEDS), \APACrefbtitleComputational Logistics Computational logistics(\BPGS502–517).\APACaddressPublisherChamSpringer International Publishing.\PrintBackRefs\CurrentBib
Zbinden\BOthers. [\APACyear2020]\APACinsertmetastarzbinden2020embedding{APACrefauthors}Zbinden, S., Bärtschi, A., Djidjev, H.\BCBL Eidenbenz, S.\APACrefYearMonthDay2020.\BBOQ\APACrefatitleEmbedding Algorithms for Quantum Annealers with Chimeraand Pegasus Connection Topologies Embedding Algorithms for QuantumAnnealers with Chimera and Pegasus Connection Topologies.\BBCQ\APACrefbtitleInternational Conference on High Performance ComputingInternational conference on high performance computing(\BPGS187–206).\PrintBackRefs\CurrentBib

Appendix A. Comparative analyses for different values of the t_min parameter on selected solutions

The right choice of the t_min parameter can improve results meaningfully.Hence, in this Appendix we analyse in more detail the role of the t_min solver parameter for selected problems in Section6.3.In particular, for networks $7$ and $9$ therein, we have performed runs with t_min sweeping over a range $[5,60]$ ; the results are presented in Fig.10and Fig.11. The objective value is roughly log-linear in computational time. For small t_min parameter values the computational time is shorter, whereas for large t_min it grows exponentially. Overall, the improvement in the objective is overwhelmed by the increase of the computational time. This is important from the point of view of the algorithm layout, as it may be tied to a power law scaling. Such a behaviour is plausible in the case of QA methods[46, 1, 20].