Planning rolling stock maintenance: Optimization of train arrival dates at a maintenance center

A railway network is an indispensable part of the public transportation system in many major cities around the world. In order to provide a safe and reliable service, a fleet of passenger trains must undergo regular maintenance. These maintenance operations are lengthy procedures, which are planned for one year or a longer period. The planning specifies the dates of trains' arrival at the maintenance center and should take into account the uncertain duration of maintenance operations, the periods of validity of the previous maintenance, the desired number of trains in service, and the capacity of the maintenance center. The paper presents a nonlinear programming formulation of the considered problem and several optimization procedures which were compared by computational experiments using real world data. The results of these experiments indicate that the presented approach is capable to be used in real world planning process.


1.
Introduction. Passenger trains provide one of the major means of public transport in many cities around the world. For example, the suburban passenger rail network in Sydney, Australia connects central Sydney with northern, southern, western, and eastern suburbs with 174 stations. The network spans over 813 kilometers of track and delivered about 359.2 million passenger journeys in 2017/18 [16].
There is a risk of sudden breakdowns or even a derailment if the trains (rolling stock) do not undergo maintenance regularly. There are several mandatory levels of maintenance that differ from each other by their scope and periodicity. Readers are referred to [9] for a detailed description of the various maintenance levels. This paper focuses on the high-level heavy maintenance which is the most involved and time consuming procedure and is performed at a specialized maintenance center. The rolling stock arrives at the maintenance center in groups. Each group is comprised of several cars coupled together and is referred to as a set or a train-set [9]. All cars in a train-set undergo maintenance in the maintenance center simultaneously.
The reliable functioning of a passenger transportation system is not possible without a plan, specifying the availability of the rolling stock. The dates when the train-sets should be withdrawn from service and sent to the maintenance center are crucial for the rolling stock operator as well as for the planning at the maintenance center. Given that the heavy maintenance of a train-set is a long process, both, rolling stock operator and maintenance center, need to specify at least for a year when the heavy maintenance of each train-set should commence [10], [15]. Furthermore, often the transportation of passengers and the rolling stock heavy maintenance are carried out by two separate organizations which operate on the basis of a long term contract that specifies for each train-set the precise date of the commencement of its maintenance.
This paper is concerned with the development of the plan that specifies the dates when the train-sets should arrive at the maintenance center. The presented optimization procedures take into account the specifics of the heavy maintenance procedure, including its uncertain duration and the restriction on the period between consecutive arrivals of the train-sets at the maintenance center, as well as the information provided by the rolling stock operator, including the engineering restrictions on the permissible period between heavy maintenance procedures and the demand for transportation.
The heavy maintenance has a validity period after which another heavy maintenance procedure must take place according to the engineering restrictions. Consequently, the rolling stock operator determines for every train-set a time window within which the maintenance of this train-set should commence (see, for example, [10]). The length of such time window depends on the practice of the rolling stock operator and in some cases can be zero (see, for example, [15]). Although it is equally possible to start maintenance at any point in time within the corresponding time window, for the purpose of the optimization procedures below, the middle of a time window will be considered as a preferred date and the time window will be viewed as a specification of the permissible deviation from this preferred date.
The dates when the train-sets should arrive at the maintenance center, specified by the mentioned above plan, may violate the time windows of the train-sets which are dictated only by the practice of the rolling stock operator and do not take into account the capacity of the maintenance center, the uncertain duration of maintenance, and the demand for transportation. Such violation is highly undesirable which often is modeled by introducing a penalty for the violation of the time windows.
A train-set that arrives at the maintenance center before its time window is referred to as "early", while a train-set that arrives after its time window is referred to as "tardy". If despite the penalties for the violation of the time windows the plan still contains some early and tardy train-sets, this situation is resolved by the consultations between the rolling stock operator and the maintenance center. Similar to the majority of publications on this topic, the consultation phase is beyond the scope of this paper.
Upon arriving at the maintenance center, a train-set is shunted to the first operation line, where thorough inspections and replacement of some components and parts such as air-conditioning units are performed. A train-set can arrive only if the first operation line is not occupied by the previous train-set. Normally, there are several types of trains and each type may require different time at the first operation line.
After the first operation line, the train-set undergoes various maintenance operations such as bogies replacement, system testing, refurbishment, etc. In order to perform these operations, the train-set has to be shunted between several lines. The duration of each operation depends on the condition of the train-set; the availability and composition of the workforce; the availability of spare parts; and many other factors. Consequently, the dwelling time of a train-set at the maintenance center (also referred to as a cycle time) is uncertain at the time of planning.
On the arrival at the maintenance center, a train-set is completely withdrawn from service and must stay at the maintenance center for at least one month. This long cycle time directly impacts the number of train-sets available in active service. Therefore, as part of the input data for the heavy maintenance planning, a permissible number of train-sets that can be out of service simultaneously are specified for each type of train-sets. Furthermore, the capacity of the maintenance center also imposes the restriction on the number of train-sets which can undergo maintenance simultaneously. Any violation of all these restrictions causes serious problems and therefore must be taken into consideration at the planning stage.
The goal of the planning process is to determine an arrival plan, specifying the arrival dates for all train-sets. The discussion above suggests that it is reasonable to have the objective function as a weighted sum of two components: the total penalty for the deviation (earliness and tardiness) from the arrival time windows, and the expected total penalty for violating the center capacity as well as for violating the permissible number of out-of-service train-sets of all types.
The body of literature on planning the rolling stock maintenance falls into two broad categories: (i) planning the rolling stock low-level maintenance subject to the trains timetable, and (ii) planning the rolling stock high-level maintenance for a long planning horizon. The first category considers only low-level maintenance, such as daily and monthly inspections, which is often combined with the decisions on the rolling stock utilization. Accordingly, the low-level maintenance is often considered as constraints in the planning process of rolling stock utilization (see for example, [9], [7]).
The second category is concerned with the high-level maintenance which has a long cycle time. In authors' opinion, the high-level maintenance planning has not attracted literature which it deserves. To the best of the authors' knowledge, only the publications [15], [5], and [10] are closely related to the planning problem considered in this paper.
The first of these three publications, [15], is concerned with maintenance scheduling at the Hong Kong Mass Transit Railway Corporation. In contrast to our paper, the authors of [15] postulate that the duration of maintenance is given. This simplification allows them to approach the minimization of the earliness and tardiness with respect to the completion of maintenance as a deterministic scheduling problem. The authors of [15] also assume that the given permissible number of trains which can dwell at the maintenance center simultaneously can not be violated, which may not be possible to ensure when the duration of maintenance is uncertain. The resultant deterministic scheduling problem is solved by a genetic algorithm. The authors reported that the proposed approach produced near optimal solutions for randomly generated instances with linear earliness and tardiness objective.
As in [15], [5] assumes that the duration of maintenance is known and that the limit on the number of trains that can dwell at the maintenance center simultaneously cannot be violated. In addition, [5] is concerned with planning under the condition that the maintenance operations must commence prior to their due dates. The planning problem was formulated as a mixed integer linear programming model and was solved using IBM ILOG CPLEX 11.2. The objective function is a weighted sum of the maintenance cost, the cost of shunting activities, the cost related to the spare parts, and the penalty for early maintenance. The authors reported that the inclusion of the cost related to the spare parts positively affected the quality of the plan.
The publication [10] is concerned with planning the heavy maintenance of the electric multiple units at the Shanghai Railway Bureau. As in [15] and [5], it is assumed that the duration of maintenance is known and that the limit on the number of trains that can dwell at the maintenance center simultaneously cannot be violated. As in our paper, each train has a time window when its maintenance should commence, but in contrast to our paper, [10] postulates that this time window cannot be violated. The problem is formulated as an integer linear program with the objective of minimizing the total penalty for early maintenance. It was reported that small instances can be solved exactly. For large instances, the authors propose a simulated annealing algorithm and report solving instances with up to 124 trainsets and a planning horizon of 607 days.
In summary, the contribution of our paper is that, in contrast to the previous publications, our paper considers the planning problem where • for each train-set, the duration of maintenance is uncertain; • the violation of the time windows within which the maintenance should commence is undesirable but possible; • there are several types of train-sets; • the violation of the given limit on the number of train-sets that can dwell at the maintenance center simultaneously is undesirable but possible; • for each type of train-sets, the violation of the given limit on the number of train-sets that can dwell at the maintenance center simultaneously is undesirable but possible.
These differences to the previous publications lead to new optimization procedures designed to be used in a complex planning process involving the negotiations between the rolling stock operator and the maintenance center. More specifically, the paper presents • a new nonlinear programming formulation of the considered planning problem, which takes into account the uncertain duration of maintenance; • a method that permits to evaluate the objective function of the introduced nonlinear mathematical program for any feasible solution without resorting to computationally expensive Monte-Carlo simulation; • a new mixed integer linear programming relaxation that is based on Jensen's Inequality [3] and therefore provides a lower bound to the optimal value of the objective function for the nonlinear mathematical program and hence allows to assess the quality of approximate solutions; • an Iterated Local Search (ILS) metaheuristic that, in contrast to the previous publications, takes into account the uncertain duration of maintenance; • a hybrid two-stage optimization procedure that combines the Jensen's Inequality based relaxation (or another mixed integer linear program also introduced in the paper) with either a local search subroutine or the presented ILS subroutine; • a fast method for evaluation of the neighborhoods for the local search and ILS subroutines; • results of the comparison, by means of computational experiments on realworld data, of the presented optimization procedures taking into account two main characteristics: the solution quality and the time needed for obtaining the solution.
The remainder of the paper is organized as follows. Section 2 presents a nonlinear mathematical programming formulation of the considered problem, an efficient algorithm for evaluation of the objective function, and a Jensen's Inequality based mixed integer linear programming relaxation. Section 3 presents an alternative mixed integer linear program, local search and iterated local search subroutines, and a fast method for the neighborhood evaluation. Section 4 presents the results of computational experiments that use real-world data provided by a major maintenance center in Australia. The conclusion can be found in Section 5.

2.
Mathematical programming formulation. This section presents a Nonlinear Integer Programming Model (NIPM) for the considered problem and a Mixed Integer Programming Model (MIPM) that provides a lower bound for NIPM. The considered planning problem can be stated as follows. A set N = {1, · · · , n} of train-sets is to undergo maintenance at a maintenance center. The planning period is T days which are numbered from 0 to T − 1. An arrival plan specifies for each train-set j ∈ N the day s j of its arrival at the maintenance center. No two train-sets can arrive on the same day. Furthermore, the train-sets are of m different types and, for each type k, there exists a restriction when the next train-set can arrive after the arrival of a train-set of type k. This restriction is given by the number of days τ k . If a train-set j of type k arrives on day s j , then, regardless of the type of the next train-set, it can arrive only on the day s j + τ k or later.
It is convenient to consider the partition F 1 , ..., F m of N , where each F k is comprised of all train-sets of the same type k. Observe that arrival days that satisfy the restriction on the time between two consecutive arrivals exist if and only if In what follows, it is assumed that this inequality holds. For each j ∈ N , the time that a train-set j spends at the maintenance center (its cycle time) is a discrete random variable D j which assumes integer values. All these random variables are independently distributed and, for each 1 ≤ k ≤ m, the cycle times of all train-sets in F k are identically distributed between a k -the minimal possible duration of a cycle for a train-set of type k, and b k -the maximal possible duration of a cycle for a train-set of type k. For each 1 ≤ k ≤ m, τ k < a k .
Each train-set j has the associated time window [θ j − ∆, θ j + ∆], where θ j is the preferred day of the commencement of maintenance, i.e. the preferred day of the arrival at the maintenance center, and ∆ is the permissible deviation from this preferred day of arrival. The penalty for the violation of time window is defined by the function: where λ 1 > 0 and λ 2 > 0. For any arrival plan σ = [s 1 , · · · , s n ], the total penalty for the violation of time windows will be denoted by G 1 (σ) = G 1 (s 1 , · · · , s n ): For any arrival plan σ, any 1 ≤ k ≤ m, and any day t, denote by W k t (σ) the number of train-sets in F k that are at the maintenance center on day t. Since all cycle times are random variables, W k t (σ) is a random variable. If W k t (σ) exceeds the given limit C kt , this attracts the penalty δ kt (W k t (σ) − C kt ), where δ kt > 0. The total number of train-sets at the maintenance center on day t is If this number exceeds the given limit C t , this attracts the penalty where δ t > 0. Let σ be an arbitrary arrival plan and G 2 (σ) be the mathematical expectation of the total penalty for the violation of the limits on the number of train-sets that can dwell at the maintenance center simultaneously. Then where E[·] is the mathematical expectation. The goal is to minimize where α and β are two positive weights. The heavy maintenance planning problem can be formulated as follows.
(NIPM) Z N IP M = min αG 1 (s 1 , · · · , s n ) + βG 2 (s 1 , · · · , s n ) subject to (1), (2), (3) In the above formulation, the objective function (6) is the weighted sum of two components: the total penalty for the violation of time windows; and the expected penalties for the violation of the limits C t and C kt . Constraint set (7) ensures that each train-set must arrive for maintenance within the planning horizon. Constraint set (8) enforces that at most one train-set can occupy the first operation line on any given day. The arrival day of each train-set can be calculated as (9). Constraint set (10) states that the decision variables are binary.

2.2.
Evaluation of the objective function. For optimization problems involving random variables, Monte-Carlo simulation is a commonly used approach for approximately evaluating the objective function [2,8] unless the model can be written in closed form by assuming some special distributions [6,1]. The objective function (6) can be evaluated exactly using the method below.
Given an arrival plan σ = [s 1 , . . . , s n ], let which is a Bernoulli random variable that takes value 1 if the train-set j is at the maintenance center on day t and 0 otherwise. Then, the probability of Y j (s j , t) = 1 can be computed as For any arrival plan σ, and any day t, W t (σ) = j∈N Y j (s j , t) is a sum of Bernoulli random variables with success probabilities according to (12). Therefore, W t (σ) is a random variable that follows Poisson Binomial distribution [4]. For any arrival plan σ = [s 1 , · · · , s n ], any 1 ≤ i ≤ n, and any day t, denote p i = Prob(Y i (s i , t) = 1). Let Prob(W t (σ l ) = i) be the probability that i train-sets from the partial arrival plan σ l = [s 1 , . . . , s l ] dwell at the maintenance center on day t. Then, Let P M F be the probability mass function of the Poisson binomial distributed variable. The entire procedure is outlined in Algorithm 1 [4].
Similarly, for any arrival plan σ, any 1 ≤ k ≤ m, and any day t, W k t (σ) = j∈F k Y j (s j , t) is a sum of Bernoulli random variables with success probabilities according to (12). Then, all probabilities Prob(W k t (σ)) can be obtained using Algorithm 1.

Algorithm 1 Direct Convolution (DC)
1: Input: total number of Bernoulli random variables n; success probability p i of the i-th random variable 2: Output: P M F 3: procedure DC(p 1 , · · · , p n ) 4: for i from 2 to n do 6: end for 14: return P M F 15: end procedure As a result, the objective function (4) can be computed as 2.3. Integer linear programming relaxation based on Jensen's Inequality. As in (3), denoting the mathematical expectation by E[·], the right-hand side in (14) can be written as and, taking into account that max{0, ·} is convex, by Jensen's Inequality [3], Using the variables x jt , introduced in (5), for any t ∈ {0, ..., T − 1}, and, for any t ∈ {0, ..., T − 1} and any 1 ≤ k ≤ m, This leads to the following Mixed Integer Programming Model (MIPM): subject to The nonlinear programming problem NIPM and the mixed integer linear programming problem MIPM use the same variables x jt , introduced in (5), that satisfy the same set of constraints (7), (8), (10). Therefore, for any arrival plan σ = [s 1 , ..., s n ] that is feasible for NIPM, there exists a feasible solution for MIPM with variables x jt , w t and w k t satisfying, for all j ∈ N , t ∈ {0, ..., T − 1} and 1 ≤ k ≤ m, or, by taking into account (16) and (17), In other words, for any arrival plan σ = [s 1 , ..., s n ] that is feasible for NIPM, there exists a solution for the mixed integer linear programming problem MIPM which variables x jt , w t and w k t satisfy This, by virtue of (15), implies that the optimal value Z M IP M of the objective function for MIPM is a lower bound on the optimal value Z N IP M of the objective function for NIPM.
3. Construction of arrival plans. As has been discussed in Section 1, the existing literature often ignores the stochastic nature of cycle times and focuses on the optimization procedures that construct arrival plans assuming that the duration of maintenance is known. This also reflects the practice often encountered in industry. Under the assumption of deterministic cycle times, the arrival plan is constructed either by solving a mathematical programming problem [10], [5] or by some metaheuristic [10], [15]. The Mixed Integer Linear Programming Problem (MILP) below also assumes that the cycle times are known constants, but in contrast to the previous publications the resulting arrival plan is evaluated as described in Section 2.2.
As far as metaheuristics are concerned, this paper presents an Iterated Local Search (ILS) algorithm which is embedded into the multi-start framework. In contrast to the previous publications, the presented ILS does not assume that the cycle times are deterministic and evaluates each element in a neighborhood taking into account the stochastic nature of cycle times.
Furthermore, this paper presents a hybrid optimization procedure which generates a starting solution by solving either the mixed integer linear programming problem MILP or the mixed integer linear programming problem MIPM and then enhances the obtained solution using the mentioned above ILS algorithm. As has been discussed above, when the hybrid optimization procedure generates the starting solution using the mixed integer linear programming problem MIPM, the optimal value of the objective function is a lower bound on the optimal value of the original nonlinear programming problem. Section 4 reports the results of computational experiments which were aimed at the comparison of the quality of the solutions produced by the presented optimization procedures and the times needed to generate these solutions. Since the considered optimization problem is a component of the complex planning process which normally involves a lot of negotiations between the rolling stock operator and the maintenance center, both characteristics, quality and time, are important.

Mixed integer linear program MILP.
It is reasonable to model random cycle times using a beta-PERT distribution [12], which is commonly used in project management [13]. In this case, the constant duration is chosen as the most likely value of the cycle time according to the beta-PERT distribution. This approach is adopted in the computational experiments below. Let q k be the most likely value of the cycle time for the train-sets in F k . Then, the arrival day for each train-set j can be determined by where variables x jt are obtained by solving the following mixed integer linear program 3.2. Local search subroutines. As has been discussed above, the considered planning problem can be solved either by some local search based metaheuristic or by a hybrid algorithm that produces a starting solution by solving one of the two mixed integer linear programs, MIPM or MILP, and then enhances this solution by some local search based optimization procedure. Four neighborhood operators, N 1 , N 2 , N 1 , and N 2 , as described below, will be used for this purpose.
3.2.1. Neighborhood operators. Consider an arbitrary arrival plan σ = [s 1 , ..., s n ] and a sequence t 1 , ..., t n of arrival days which are listed in a nondecreasing order and which are obtained by changing a single arrival day in σ, say by changing the arrival day s g of some train-set g. For each 1 ≤ j ≤ n, let n(j) be the train-set that arrives at t j , and let k(j) be the type of n(j). Let t i be the new arrival day assigned to the train-set g. The replacement of s g by t i results in a new feasible arrival plan if The neighborhood explored by the operator N 1 is comprised of all feasible arrival plans (solutions) that can be obtained from the input arrival plan σ by assigning a different arrival day to a single train-set. In other words, the neighborhood explored by the operator N 1 is comprised of all arrival plans that are obtained when the change of a single arrival day in σ results either in (p1), or (p2), or (p3). The benefit of using N 1 is in the fast evaluation of each solution in the neighborhood because each solution is obtained by changing only a single arrival day and this change does not affect any other arrival days. The operator N 2 explores all solutions that can be obtained by changing any two arrival days of the input arrival plan. Similar to N 1 , these two changes must not affect any other arrival days in the input arrival plan.
Other two operators N 1 and N 2 are similar to N 1 and N 2 but their neighborhoods are constructed without the restriction that the change of an arrival day in σ does not affect any unchanged days in σ. More specifically, the operator N 1 explores the neighborhood comprised of the solutions that result from a change of a single arrival day in the input arrival plan, but in contrast to N 1 , this change is allowed to violate all three feasibility conditions (p1), (p2), and (p3). In such case, the set of arrival days is transformed into a feasible arrival plan by the algorithm TRANSFORMATION and its two subroutines LEFT and RIGHT. Similarly, the operator N 2 explores all solutions that can be obtained by changing any two arrival days in the input arrival plan, but in contrast to N 2 , these two changes may affect some other arrival days in the input solution. In such case, analogously to N 1 , the resultant set of arrival days is transformed into a feasible arrival plan using the same algorithm TRANSFORMATION.  for j = i + 1; j ≤ n; j++ do Let N be one of the four operators, N 1 , or N 1 , or N 2 , or N 2 . For each solutioñ σ = [s 1 , ...,s n ] in the neighborhood of an input arrival plan σ, the operator N computes the value of the objective function Letσ be a solution in the neighborhood of σ with the smallest value of the objective function. Denote by N(σ) the output solution produced by the operator N. Then The Algorithm 5 below outlines the local search procedure for an input arrival plan σ. If N is N i , where i ∈ {1, 2}, then this procedure will be referred to as LSi. If N is N i , where i ∈ {1, 2}, this procedure will be referred to as LSi .

Algorithm 5 Local Search
1: repeat 2:σ = σ 3: σ = N(σ) 4: until αG 1 (σ) + βG 2 (σ) == αG 1 (σ) + βG 2 (σ) 5: returnσ The four search operators, mentioned above, can be combined in different ways. In this paper, four options are considered: (1) LS1 is a local search that uses only the operator N 1 ; (2) LS1 is a local search that uses only the operator N 1 ; (3) Sequential Local Search (SLS) applies LS1 to the input arrival plan and when LS1 terminates applies LS2 to the output of LS1; and (4) SLS applies LS1 to the input arrival plan and when LS1 terminates applies LS2 to the output of LS1 . Each option has its own merit as demonstrated in the computational study.

Evaluation of solutions in a neighborhood.
The main computational burden in many local search algorithms is the evaluation of the elements constituting a neighborhood. The local search procedures LS1, LS2, LS1 and LS2 are no exception. Therefore, an algorithm for computing (26) for each element of the neighborhood of an input solution is the key factor that determines the efficiency of these optimization procedures.
As (14) indicates, in order to recompute (26) efficiently, one requires a fast algorithm for recomputing the probability mass functions of the random variables W t , W k t for all t ∈ {0, ..., T − 1} and 1 ≤ k ≤ m. All these random variables have the same nature -they are sums of Bernoulli random variables. Given the structure of neighborhoods explored by the operators N 1 and N 2 , each element in the neighborhood of the input arrival plan can be obtained by changing the arrival day of one or two trains-sets in this input solution. For each such train-set, the change in the arrival day results in the change of one of the Bernoulli random variables in the sums defining the random variables W t , W k t for all t ∈ {0, ..., T − 1} and 1 ≤ k ≤ m.
This replacement can be viewed as the elimination of one of the Bernoulli random variables from the sum followed by the addition of another Bernoulli random variable to the result of the elimination. The Algorithm 6 and Algorithm 7 below do this efficiently for any probability mass function of a random variable that is a sum of Bernoulli random variables, and therefore, are used for recomputing probability mass functions of the random variables W t , W k t for all t ∈ {0, ..., T − 1} and 1 ≤ k ≤ m. Let p be the success probability of the Bernoulli random variable that is to be eliminated, P M F be the original probability mass function, and new P M F be the probability mass function resulted from the elimination of this Bernoulli random variable. The probability mass function new P M F is defined for all integers 0 ≤ i ≤ n − 1, whereas the probability mass function P M F is defined for all integers 0 ≤ i ≤ n. end for 14: end if 15: return new P M F Let p be the success probability of the Bernoulli random variable that is to be added to a sum of Bernoulli random variables which has the probability mass function P M F . Let new P M F be the probability mass function of the sum after this addition. The probability mass function new P M F is defined for all integers 0 ≤ i ≤ n, whereas the probability mass function P M F is defined for all integers 0 ≤ i ≤ n − 1. The reasoning similar to the above lead to the Algorithm 7.
3.3. Iterated local search. Iterated local search (ILS) is one of the commonly used metaheuristics which was successful in solving a wide range of optimization problems [11]. The Algorithm 8 below outlines this metaheuristic as it was implemented and used in the computational experiments, the results of which are reported in Section 4. In this pseudocode, G is the objective function (4) and the parameter U specifies the maximum permissible number of consecutive unsuccessful attempts to improve the current best known arrival plan σ * .
The Algorithm 8 interchangeably invokes two subroutines, SEARCH and PER-TURB. In some computational experiments, reported in Section 4, the subroutine SEARCH is a local search procedure LS1 or LS1 (see Algorithm 5), whereas in the others, the subroutine SEARCH is the sequential local search SLS or SLS . The subroutine PERTURB randomly chooses three train-sets and one by one assigns to them new arrival days without violating the feasibility (that is, without violating the restrictions on the duration of time intervals between any two consecutive arrivals of train-sets). In the process of assigning the arrival days to the three selected train-sets, the subroutine PERTURB does not take into account the new value of the objective function G. In what follows, for any arrival planσ, SEARCH(σ) and PERTURB(σ) denote the output of SEARCH and PERTURB, respectively, resulted from their application toσ. The arrival plan σ in SEARCH(σ) in line 1 is a feasible (in terms of the time between the consecutive arrivals) input solution.
An input arrival plan for the ILS is generated either by solving one of the two mixed integer linear programs, MIPM or MILP, or by the heuristic INITIAL described below. The heuristic INITIAL chooses randomly n different arrival days t 1 < ... < t n and associates randomly with each t i one of m types of train-sets in such a manner that each type 1 ≤ k ≤ m receives |F k | arrival days. Denote by k(t i ) the type associated with t i . If, for each 1 ≤ i < n, then the arrival days are feasible. If they are not feasible, then the heuristic INITIAL transforms the generated arrival days into feasible arrival days, using the end for 9: end if 10: return t 1 , ..., t n 4. Computational results. The proposed solution approaches are tested and evaluated on data provided by one of the leading maintenance centers in Australia. The planning horizon is one year. There are 35 train-sets of 3 different types. The parameters of train-sets in F k , 1 ≤ k ≤ 3 are given in Table 1. The column '|F k |' reports the total number of train-sets of the same type k; the column 'τ k ' gives the number of days a train-set of the corresponding type spends on the first operation line. Since an increase in transport demand is often expected during public holidays, the number of train-sets of type k which can dwell at the maintenance center simultaneously during these days is normally less than that on any other days of the year. The out-of-service limit 'C kt ' for type k on day t is reported in column 'non-PH' if day t is not a public holiday, and in column 'PH' if day t is a public holiday, where PH is an abbreviation of public holiday. For any day t, the capacity of the maintenance center imposes a limit C t = 5, and the daily penalty factor δ t is fixed at 1. For any 1 ≤ k ≤ 3, and any day t, the daily penalty factor δ kt = 1 if t is not a public holiday, and δ kt = 10 if t is a public holiday. A larger penalty is applied for public holidays because it is more undesirable to violate the given limit C kt during these periods. With a larger penalty factor δ kt , the violation of the permissible number of train-sets of type k that can dwell at the maintenance center simultaneously on public holidays is still possible in the optimal solution.
The permissible deviation from the preferred day of the commencement of maintenance is 14 days, i.e. ∆ = 14. The penalty factors for the violation of time windows are chosen as λ 1 = λ 2 = 1.
As has been discussed in Section 3.1, the random cycle times are modeled using beta-PERT distribution with the minimum, most likely, and maximum values as given in Table 2 for each train type k. Note that the beta-PERT distribution is a continuous probability distribution. So for the computational experiments, beta-PERT distribution is discretized into days.  1  20  25  40  beta-PERT  2  27  30  46  beta-PERT  3  29  30  52 beta-PERT Extensive experiments were performed on numerous settings for the weights α and β. The choice of values presented in Table 3 corresponds to the appropriate weight coefficients that can generate solutions representing typical scenarios. For case 1, a large relative weight is assigned to the second component of the objective function G 2 (σ), and the optimization procedures aim at minimizing the expected penalties for the violation of the limits C t and C kt . The importance of the objective G 1 (σ) increases when proceeding from case 1 to case 9. All algorithms are implemented in Python 2.7. The mixed integer linear programs are solved with IBM ILOG CPLEX 12.7 via the mathematical programming modeling language PuLP [14]. All tests are run on a computer with Intel i5-6300U 2.4GHz processor and 8GB of RAM.

4.1.
Comparison of the performance of MIPM and MILP. Table 4 compares the results of the two mixed integer linear programs. As has been discussed in Section 2.3, solving the MIPM produces a lower bound which is reported under the column titled 'LB'. The column 'Time' gives the running time (in seconds) by CPLEX to obtain an optimal solution. For all cases, the objective values of the solutions are computed as described in Section 2.2 and reported under the column titled 'Obj.'. The relative gap is calculated as %Rel = (Obj. − LB)/LB × 100.
The results in Table 4 shows that the MILP can be solved in short computation time, yet the MIPM has the advantage of providing a lower bound. It is observed that the MIPM gives better solutions in 6 out of 9 cases, and on average 0.54% better than the MILP. For all cases, CPLEX obtains an optimal solution to MILP in less than 11 seconds, whereas more time is needed to solve the MIPM to optimality. By investigating the output of CPLEX in case 1 which takes the longest time, it was found that the optimal solution in fact was obtained in less than 2 minutes. The remaining time was taken by CPLEX for proving optimality, which is a common behavior of this software. As the relative weight of β decreases, the computation time of MIPM reduces substantially. This observation suggests that the second component of the objective function is harder to optimize for the MIPM. The improvements in solution quality by the hybrid two-stage optimization procedure that combines the mixed integer linear program with the local search LS1 and the Sequential Local Search (SLS) are reported in Tables 5 and 6, respectively. In both tables, the column 'LS Time' gives the computation time (in seconds) of the local search procedure. Results in Table 5 show that the local search LS1 is effective in improving the solutions of both models for all cases. On average, the relative gap is reduced by 30.56% for MIPM and 29.71% for MILP. The behavior of the local search is consistent for both models: the best performance is achieved in case 5 with a change in the objective value of 11.48% for MIPM, and 11.76% for MILP; whereas the worst performance is observed in case 8 at 0.004% and 0.23%, respectively. The   Table 6 show that the MIPM with SLS yields better solutions in shorter running times. For MIPM, the search in the neighborhood explored by the operator N 2 finds better solution in only one of the nine cases, i.e. case 7. For MILP, the SLS is seen to be especially useful on cases with large relative weight β. Since the MIPM produces better results over the MILP for most cases, MIPM is used in the remainder of the computational experiments.

4.2.
Comparison of hybrid ILS and multi-start ILS. In this section, eight algorithms, namely hybrid ILS with LS1, hybrid ILS with LS1 , hybrid ILS with SLS, hybrid ILS with SLS , multi-start ILS with LS1, multi-start ILS with LS1 , multi-start ILS with SLS, and multi-start ILS with SLS , are compared by means of computational experiments. For all eight algorithms, a maximum permissible number of iterations without improvement U = 20 is used. The computational results are reported in Tables 7 -10. In the preliminary testing, performing an exhaustive search on the neighborhood by N 2 and N 2 is too time consuming, which significantly reduces the number of perturbation in the ILS procedure due to the time limit of 1 hour. As a result, the ability of ILS to escape the local optimum is severely impaired. Therefore instead of allowing the arrival days of two train-sets to be changed to any feasible days from {0, ..., T − 1}, the newly assigned arrival day t of each such train-set j is selected from the range s j − 36 ≤ t ≤ s j + 36, where s j is obtained from the input arrival plan σ. Such restriction substantially reduces the size of the search space so that the neighborhood can be explored in a reasonable time. The reduced structure of neighborhood explored by N 2 and N 2 is employed in the computational experiments of hybrid ILS with SLS, hybrid ILS with SLS , multi-start ILS with SLS, and multi-start ILS with SLS . We will discuss the results of hybrid ILS first.
The hybrid ILS is implemented by solving the MIPM, which provides an input arrival plan to the iterated local search in Algorithm 8. Because of the faster evaluation of solutions in a neighborhood and the smaller neighborhood size, the versions of the hybrid iterated local search with the operators N 1 and N 2 are faster than their counterparts with the operators N 1 and N 2 often at a cost of inferior solution quality. The computational experiments took this into account and ran versions with the operators N 1 and N 2 with one hour limit on the permissible computation time, recorded for each such optimization procedure the average of the actual computation times for ten runs, and then set this recorded average time as the limit on the computation time for the version with the corresponding operators N 1 and/or N 2 .
For each case, i.e. for each choice of the parameters α and β, Tables 7 and 8 present the average computation time in seconds (Time), the average value of the objective function (Obj.), and the average relative gap (%Rel) obtained for ten runs of the hybrid ILS. The column In. Obj. displays the value of the objective function obtained by solving the MIPM. Table 7 indicates that the version with N 1 obtains better quality solutions in six of the nine cases, with the average relative gap improving from 9.12% to 5.84%. Table 8 also indicates the superior performance of the version with N 1 and N 2 in comparison with the version with N 1 and N 2 .
The output of the multi-start ILS is the best output obtained by the application of the Algorithm 8 to five different input arrival plans generated by the heuristic INITIAL. These five applications of the Algorithm 8 constitute one run of the multistart ILS. In the course of the computational experiments, the duration of each run of the multi-start ILS was limited by one hour. For each case, i.e. for each choice of the parameters α and β, Tables 9 and 10 present the average required time (Time), average value of the objective function (Obj.), and average relative gap (%Rel) obtained for ten runs of the multi-start ILS. The column In. Obj. contains the average value of the objective function for the input arrival plans that resulted in the output of the multi-start ILS.
In Table 9, the multi-start ILS with LS1 is superior to the multi-start ILS with LS1 in both time and solution quality. The same observation can be seen in Table  10, in which the multi-start ILS with SLS showing better performance over the version with SLS . It is worth noting that although the multi-start ILS algorithms begin     with poor initial solutions, significantly better results are achieved after the local search based enhancement procedures, with the best reported average improvement of 72% observed in the multi-start ILS with LS1. In summary, the comparison of the eight optimization procedures indicates that the best solution quality is obtained by the combination of the mixed integer linear program MIPM and the iterated local search with the operators N 1 and N 2 . The comparison of the different neighborhoods is summarized in Table 11. The column Tables contains references to the tables that present the results of computational experiments. The column N contain the number of cases for which the neighborhood structure that does not require the subroutine TRANSFORMATION yields a better solution quality. The column N contain the number of cases for which the neighborhood structure obtained, using the subroutine TRANSFORMATION, yields a better solution quality. The column EQUAL contains the number of cases when both neighborhood types produced the same solution quality.

4.3.
Visualization of quality of arrival plan. During the negotiations between the rolling stock operator and the maintenance center, it is useful to have information about the risk of violating the center capacity which may occur as a result of the uncertain duration of maintenance. For this reason, a powerful visualization tool, based on the idea of heat map [17], is developed to provide insights into the risk over the planning horizon. Figures 1A and 1B show examples of the risk heat map associated with the arrival plans of cases 1 and 9, respectively. The horizontal axis indicates the month and year (for example, 2018-06 stands for June 2018), while the vertical axis indicates the day of the month. Each cell of the heat map  Figure 1. Heat maps displaying the probability of having more than 5 train-sets residing in the maintenance center for each day across the planning horizon of one year for cases (A) α = 1, β = 1000; and (B) α = 1, β = 1.
corresponds to a particular day in the planning horizon, and the probability of violating the limit is clearly stated in each cell. The color intensity reflects the level of risk whereby the darker the color, the higher the risk.
The resulting heat map in Figures 1A and 1B show a trade-off example in which the constructed arrival plan must prioritize either the technological restrictions of the maintenance center or the arrival time windows. Figure 1A considers the perspective of the maintenance center who is more concerned about keeping the number of train-sets below the capacity of the maintenance center. As a result, the total penalty for the violation of the capacity limitation is insignificant and it can be seen on Figure 1A that there are few days which have high probability of exceeding the center capacity. However, the total penalty for the violation of time windows, G 1 (σ), is 23,487. On the other hand, the heat map in Figure 1B is associated with an arrival plan σ that is constructed considering the perspective of the rolling stock operator, the main concern of whom is to satisfy the arrival time windows. In this case, the total penalty G 1 (σ ) is only 6,218 but the maintenance center has a high risk of violating the capacity, i.e. it is harder to have an efficient operational plan.

5.
Conclusion. The paper contributes to the existing body of literature on train maintenance by introducing a nonlinear programming problem that determines the arrival days of train-sets to a maintenance center, taking into account the uncertain duration of maintenance and the requirements specified by the rolling stock operator as well as the technological restrictions of the maintenance center. A fast method of evaluation of the objective function for any feasible solution of the nonlinear program is presented together with a mixed integer programming relaxation based on Jensen's inequality. This relaxation provides a lower bound on the optimal value of the objective function of the nonlinear program and generates a starting solution for the hybrid optimization procedure which enhances this solution by using iterated local search. This hybrid optimization procedure is compared with iterated local search metaheuristic diversified by the multi-start framework. The results of the computational experiments on real-world data warrant the implementation of the presented approach in the process of maintenance planning. Further research should be focused on the operational level of the maintenance planning for a shorter planning horizon and more detailed information about maintenance procedures.