Variable and adaptive neighbourhood search algorithms for rail rapid transit timetabling problem

Variable and adaptive neighbourhood search algorithms for rail rapid transit timetabling problem

Computers & Operations Research 78 (2017) 439–453 Contents lists available at ScienceDirect Computers & Operations Research journal homepage: www.el...

1MB Sizes 22 Downloads 74 Views

Computers & Operations Research 78 (2017) 439–453

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

Variable and adaptive neighbourhood search algorithms for rail rapid transit timetabling problem Erfan Hassannayebi, Seyed Hessameddin Zegordi n Industrial Engineering Department, Tarbiat Modares University, Tehran, Iran

art ic l e i nf o

a b s t r a c t

Available online 31 December 2015

Supplying affordable and efficient transportation services to the users is one of the main tasks of the public transport systems. In this study, the objective is the minimization of the total and maximum waiting time of the passengers through optimization of the train timetables for urban rail transit systems. For this purpose, mixed-integer linear and non-linear programming models are developed which could solve the small to medium-sized test instances optimally. In order to tackle large instances, adaptive and variable neighbourhood search algorithms are designed based on different novel solution encoding schemes and decoding approaches. The effectiveness of the proposed models and solution methods are illustrated through the application to the Tehran intercity underground rail lines in IRAN. The outcomes demonstrate that the variable neighbourhood search algorithm outperforms the adaptive step-size neighbourhood search method in the different scenarios of the real case. Furthermore, the generated headway for the period of study result in a significant reduction in total waiting time of the passengers compared with the current baseline timetables. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Urban rail Waiting time Adaptive step-size Variable neighbourhood search

1. Introduction An important aspect of the public transportation systems is the passenger waiting times. Specifically, inefficient transportation systems incur additional waiting time to the passengers. Since it is impossible to avoid waiting time completely, it is an essential concern in the daily operation of a public transit system to minimize the passenger waiting times. On the other hands, managing the growing travel demand for public transport is another challenge that urban rail transit systems currently face. To deal with these challenges, the railway companies should try to optimize the rail operations. The optimization of the urban public transportation systems has been intensively studied from both cost minimization and profit maximization viewpoints. In recent researches, passenger waiting time has appeared as an important performance measure to be minimized by optimization methods [50]. Train timetabling problem (TTP) in an urban rail system refers to the procedure of determining the optimal departure times of transit vehicles according to the dynamic behavior of the travel demand in order to minimize the operator and passenger performance measures [31]. In this regard, prediction models aim to analyze the main influencing factors of passenger waiting time. n

Corresponding author. Tel.: þ 98 2182883394; fax: þ98 2188005040. E-mail addresses: [email protected] (E. Hassannayebi), [email protected] (S.H. Zegordi). http://dx.doi.org/10.1016/j.cor.2015.12.011 0305-0548/& 2015 Elsevier Ltd. All rights reserved.

Here, as in early studies in this area, Welding [47], and Osuna and Newell [34] explored the characteristics of passenger waiting time and developed an estimation model for average waiting time per passenger (AWT) at station stops. They defined the expected waiting time as a function of the average and variance of the headway. The rail managers generally offer a limited number of vehicle services during the daily normal operation, due to the operational and maintenance costs. In this regard, the train headways are affected by the total number of services, fleet size, and vehicle capacity and therefore, better distribution of the train services result in decreasing the travel time and waiting time of passengers. However, urban rail transit is a complex dynamic traffic system, making it challenging to be optimized [51]. Furth and Muller [21] presented methods for generating the waiting time distribution from headway data in order to measure the service reliability. The problem of estimating passenger waiting times and measuring the headway variance at bus stops with incomplete input data was studied by Mcleod [30]. Usually, the establishment of the regular timetables results in minimum waiting time at stops. However, from the application point of view, the evenheadway or regular schedules suffer from the lack of enough agility to handle variations of travel demand over the study period [13]. Moreover, the total waiting time during periods of peak demand is expected to be affected by fleet size and train capacity constraints. Very few academic efforts have been devoted to assessment of these factors on waiting times. Islam et al. [24]

440

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

analyzed the effect of the headway variability and vehicle capacity on service performance and reliability of a high-frequency bus transit system through a Markov chain model. Amin-Naseri and Baradaran [6] extended the Welding's formulations and developed a more accurate formula to estimate the AWT at a bus stop with independent headways. A discrete-event simulation method has been applied to evaluate the accuracy of the proposed formulations. Ceder and Philibert [12] presented a method for the generation of timetables in order to eliminate the overcrowding through adjustment of departure times. The developed method produces schedule for vehicles reaching to an even maximum load. The results of this scheduling policy demonstrate the significant improvement over even-headway schedules or with headway plan constructed based on hourly even maximum loads. Although there is a rigorous literature dealing with train scheduling in the context of off-line planning, there is still a lack of efficient solution methods that can efficiently solve the demandoriented TTP. Furthermore, a limited number of studies have been found that addressed urban train timetabling problem under dynamic demand and train capacity constraint. The motivation for this study is the development of effective meta-heuristic algorithms for providing solutions to large-sized instances of the demand-oriented TTP. Also, this study aims to contribute to establish new solution encoding–decoding methods for the urban rail transit scheduling. The outline of the study is organized as follows. A review of the related literature is presented in the next section. The research objectives and contributions are provided in Section 3. The train timetabling models are described and mathematically formulated in Section 4. The proposed neighbourhood search meta-heuristics are presented in Section 5. The specification of the real case is given in Section 6. Computational results and the corresponding discussion are presented in Section 7. Finally, in Section 8 conclusions are drawn and recommendations for further research are given.

2. Literature review The optimization approach to timetables design aims to maintain the desired level of service to passengers. Moreover, the minimization of the passenger waiting time in rail transport services is a challenging scientific field that is devoted to optimization models and solution methods. The studies reviewed in this article are in this context a more effective way of improving the quality of service compared with the prediction models. The first optimization model for underground rail transit system is developed by Cury et al. [17] taking into account the train and passenger actions, service levels for passengers, the fleet size, and the trains performance. They developed a Lagrangian-based multilevel hierarchical optimization algorithm to solve the schedule optimization problem. Chierici et al. [15] developed a mixed-integer non-linear optimization model to solve the real instances of rail network in Italy. The model generates regular train timetable using a branchand-bound as well as a heuristic algorithm in order to maximize the total demand supplied by the vehicles. The service frequency directly influences the waiting time of passengers. Albrecht [4] designed a two-level method for optimizing services in the railway train timetabling problem. In the first phase, the optimal frequency of trains per line and train capacity will be determined. In the second phase, a detailed schedule of trains will be produced. Zhao et al. [53] addressed real-time train dispatching problem for urban rail transportation systems capturing the passenger flow data. The objective of the optimization model is minimization of the departure time interval subject to the minimum departure interval, the maximum acceptable waiting time, and the capacity of trains.

Niu and Zhang [32] addressed train timetabling problem under the multi-period and uniform demand condition. The train capacity is formulated as a soft constraint and the objective function includes minimizing the waiting times of passengers and the penalty costs. A hybrid genetic algorithm is proposed for solving the non-linear integer programming model. Peña-Alcaraz et al. [35] developed mathematical optimization models to design timetables in underground rail system with the aim of synchronizing the trains movements reducing the energy consumption. Canca et al. [10] addressed the demand-oriented train scheduling problem. They developed a mixed-integer nonlinear formulation to determine the number of train services and their departure and stop times. The quality of train schedule is examined through average waiting time and train occupancy ratio. Sato et al. [36] presented rescheduling approach in order to minimize the total traveling time and the waiting time at stations. Niu and Zhou [33] developed a mixed integer nonlinear programming model for the capacitated train timetabling problem in an urban transit system, subject to fleet size constraint. A heuristic search approach based on local improvement and dynamic programming approaches are proposed to solve train timetabling problem optimally for singlestation cases. The heuristic algorithm uses the gradient information when adjusting the departure time of a train during the peak and off peak periods. In order to solve the practically-sized instances of the problem, a genetic algorithm is suggested. Sun et al. [41] adopted Lagrangian duality theory to optimize train schedule for metro line taking into account the robustness and energy saving aspects. The proposed train scheduling model aims to minimize the waiting time of passengers and operational cost subject to the passenger demand, headway, and dwell time constrains. Mao and Li [28] presented three schedule adjustment strategies for the train dispatching problem. A novel algorithm is proposed to calculate the maximum traffic capacity in railway lines and a rescheduling method is employed to adjust the freighttrain timetables. Wang et al. [45] proposed mathematical models to minimize the total passenger waiting time plus passenger traveling time. Since they assumed that the passenger arrival rate is static over the study period, the problem is formulated as a nonlinear non-convex programming model and it has been solved by SNOPT solver of GAMS. The decision variables are the optimal departure times, running times, and dwell times of trains along a single-way railway line. Sun et al. [39] formulated optimization models to design demand-oriented train schedules. The mathematical programs include models without capacity constraints and models with limited train capacity. It is assumed that trains are restricted to depart at the equivalent discrete time intervals. Furthermore, the operation costs are given and the number of daily train services is constant. The proposed mathematical models are solved by CPLEX and a quantitative analysis is conducted to investigate the sensitivity of optimal total waiting time of passengers by changing the important model parameters on a metro line in Singapore. Tong et al. [42] developed a method to determine transit frequencies considering the passenger flow along the transit line. The presented approach uses the information of the un-served passengers due to capacity constraints in peak hours to adjust the frequencies. In order to generate train schedules with minimum passenger average waiting time, Barrena et al. [7] proposed a branch-and-cut algorithm with a set of valid inequalities which tightened the formulation. The proposed exact method could not yields optimal solutions for test instances in a reasonable running time but it reduces the average waiting time comparing with the evenheadway timetable of around 30% on average. According to the literature classification of the related studies (Table 1), the most of the articles optimally solved small and medium-sized instances of the problem. Also, due to the dynamic

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

441

Table 1 Summary of current literature available on demand-adapted train timetabling problem. Reference

Objective functions

Zhao et al. [53]

Minimization of the departure time interval



Dynamic



Niu and Zhang [32]

Waiting times of passengers and the in-train crowded costs

Genetic Algorithm (GA)

Dynamic



Niu and Zhou [33]

Total weighted waiting time

Local improvement and GA algorithms

Dynamic



Wang et al. [45]

Total waiting time and travel time

SNOPT solver of GAMS

Static



Canca et al. [10]

Total average waiting time

Epsilon-constraint method

Dynamic



Niu and Tian [31]

Heuristic algorithms

Dynamic



Barrena et al. [7] Sun et al. [39]

Bi-objectives (total waiting time and the number of required vehicles) Multi-objectives (energy consumption, total waiting time and travel time, total number of left behind passengers) Average waiting time Average waiting time

Barrena et al. [8]

Average waiting time

Wang et al. [46]

Total traveling time of passengers and the energy consumption

Wang et al. [44]

Solution approach

Travel demand Train capacity

Two-level optimization method using GA Static



Branch and Cut (B&C) CPLEX solver

Dynamic Dynamic



Adaptive large neighbourhood search (ALNS) Bi-level optimization method

Dynamic

√ 

Static



Hassannayebi et al. [23] Expected waiting time per passenger

GA

Dynamic

Sun et al. [40] This study

Lagrangian decomposition Adaptive and variable neighbourhood search algorithms

Static Dynamic

√ 

Waiting time of passengers and train operation cost Total and maximum waiting time of passengers

nature of the travel demand, most of the mathematical formulations were mixed-integer nonlinear and therefore intractable in practice. Although, we have found studies attempting to develop linear formulations. Moreover, the demand-driven train timetabling problem with capacity constraint is proven to be NP-hard and computationally difficult to be solved [39]. Therefore, heuristic and meta-heuristic approaches are suitable when solving NP-hard problems. A wide range of meta-heuristic algorithms were used to solve the train scheduling problem. Abril et al. [1] presented distributed and asynchronous search algorithms with application to the train scheduling problem based on an innovative meta-tree Constraint Satisfaction Problem (CSP) structure. Cucala et al. [16] addressed the optimal train operation model with the aim of minimizing the energy consumption subject to the uncertain delays caused by the driver response in high speed rail. The train timetabling model was formulated as a fuzzy linear programming model and a Genetic Algorithm (GA) was proposed to solve the problem. Jamili et al. [25] proposed a hybridized meta-heuristic algorithm based on integration of simulated annealing (SA) and particle swarm optimization (PSO) methods for a periodic singletrack train timetabling problem. Domínguez et al. [18] applied a multi-objective particle swarm optimization algorithm to produce the optimal speed profiles in Madrid's underground rail system. Barrena et al. [8] presented an adaptive large neighbourhood search (ALNS) method to minimize the average passenger waiting time at the stations. In this study, efficient variants of neighbourhood search algorithms are proposed for the rapid rail transit scheduling problem with the total and maximum waiting time criteria. Meanwhile, new decoding–encoding approaches are employed to improve the performance of the search procedure. Moreover, nonlinear and linear mathematical optimization models are presented to find global optimal solution under the given operational constraints. The results demonstrate the successful application of the proposed variable and adaptive neighbourhood search algorithms to the train timetabling problem.

3. Research objectives and contributions This paper aims to achieve two specific objectives by proposing effective neighbourhood search algorithms. First, the proposed approach improves on current solution methods of the demandoriented train scheduling problem by developing different efficient



neighbourhood search algorithms that can generate train timetables with minimum passenger waiting time in a reasonable computational time. The current study aims to find an optimized daily train timetable for a metro system through minimizing the total waiting time of passengers. Second, this paper explores the advantages of different solution representation schemes and decoding methods. For this purpose, real-valued and partial discrete frequency-based representations of the solutions are suggested and tested on a real-world underground rail system. In the real-valued solution encoding procedure, a variable step-size neighbourhood search algorithm (called adaptive neighbourhood search) is suggested to solve the train timetabling problem. The suggested encoding scheme could effectively prevent obtaining infeasible and sub-optimal solutions through adaptive moves in the neighborhoods, thus improving the search efficiency and solution quality. The discrete representation incorporates biased random search using demand information at periods to enhance search performance. Furthermore, the solution decoding is handled through solving a mathematical optimization model.

4. Problem statement and formulation Demand for rail travel in both time and space is variable. In a subway system, train services required to operate within the specified safety time and the separation time between two successive departures of the trains known as headway. In this study, the aim is to optimize the operation of train services through adjustment of the headway times in accordance with the dynamic passenger demand. As stated, the passenger waiting times depend on the number of trains, fleet capacity, operational characteristics of the route and the headway times. A typical viewpoint is that the greater the number of trains and their capacity, the shorter the passenger waiting time will be. In another stand point, a cost effective approach to better utilization of the fleet is the adjustment of headway times according to the travel demand patterns. With a regular schedule or even-headway timetable, vehicles arrive and depart at fixed intervals. Under conditions of dynamic demand, an even-headway policy may result in the long passenger waiting time during peak hours, due to the capacity constraints, or inefficient use of the capacity of the vehicles during the off-peaks. In order to measure the desirability of the public transportation system, this study addressed both the average and maximum waiting time of passengers.

442

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

Fig. 1. A time-station diagram of train services.

4.1. Total waiting time minimization approach In this section, mixed-integer nonlinear and linear formulations are presented for the capacitated train timetabling problem under multi-period demand. Before presenting mathematical formulations, some definitions and notations are introduced. Consider a set of train services (i A I) that are required to be scheduled during the period of service ½0; T . The period of service refers to the duration of time that trains provide transport service to the passengers. The total number of available train services is denoted by D. Let H i be the headway between i-th and (iþ 1)-th train departures. Thus, H 0 and H D represent the dispatch time of the first train from vehicle depot and the distance between the end of the period (T) and the last departure, respectively (Fig. 1). Let bi and wi denote the number of on-board passengers of i-th train and the number of waiting passengers after the departure of i-th train, respectively. The Eq. (4) states that the technological constraint and signaling system ensure that vehicles are operating at their most efficient as well as their most safe condition. To control the maximum waiting time of passengers, the maximum headway time between train services are incorporated in the formulation. A common and important characteristic of travel demand profile is the existence of some peaks point (local maximum) over time. Demand peaks are related to the rush hours and normally condensed to two points per day. In our approach to the passengers demand, a piece-wise linear approximation of the demand function is proposed which enables capturing demand profiles with different shapes and patterns. In this regard, the scheduling cycle is divided into a number of demand periods (p A P) that ðpÞ correspond to different demand patterns, where λ denotes the arrival rate of passenger (passenger per minute) during period p A P. Total number of demand periods is denoted by NT. Furthermore, the length of each demand period and the start of the pth period are denoted by θp and t p , respectively. The arrival process of passengers is assumed to be a homogeneous uniform for each period. The decision variables are the departure times of all train services with the aim of minimizing total waiting time of passengers. The binary variables (yðpÞ i ) are included in the formulation to represent the assignment of train services to demand periods. The variable yðpÞ i takes value 1 if i-th train runs at period p, and takes value 0 otherwise. This variable is linked with departure time variables (di ) through constraint (5). The total number of passengers arriving between two consecutive trains i and i þ1 (δi ) is dependent on not only the headway time but also the departure time of train services. Eq. (6) ensures that each train service must be scheduled in a specific time period. The assignment of train services to the demand periods should be respected in chronological order which is expressed in Eq. (7). The input flow of passengers between t ¼ 0 and the departure time of i-th train is denoted by Di . With this notation, Eq. (8) calculates the passenger flows between two consecutive departures. The flow conservation constraint is presented in Eq. (9). Eq. (10) states the train capacity constraint. Assume that i-th train departs at period t. In order to calculate Di we count the number of passengers arrived at P the periods p0 ¼ 1; 2; …; p  1. As a result, Di is equal to t0 A P   0 Pt 0 ðt Þ ðt Þ 1  t ¼ 1 yi :θt0 :λ plus the number of passengers arrived at

the beginning of the period p. Eqs. (11) and (12) are presented to calculate Di . For special cases, Eqs. (13)–(16) refer to the number of passengers arriving to station before the first and after the last train services. In addition, Eq. (17) states the initial conditions. Finally, the capacitated train timetabling model under dynamic demand is written as a mixed-integer nonlinear programming (MINLP) model with the following equations:  X  1 ½MINLP: min z ¼ w þ δ ð1Þ :H i i i iAI 2 Subject to: X H ¼T iAI i

ð2Þ

H i ¼ di þ 1 –di

i A I=fDg

hmin r H i r hmax 

ð3Þ

i A I=f0; Dg 

ð4Þ 

t p  1  M U 1  yði pÞ r di o t p þ M U 1  yði pÞ X

yðpÞ pAP i

¼1



iA I; p A P

iAI

ð5Þ ð6Þ

  0 yiðpþÞ1 þ yði pÞ r 1p; p0 A P; p ¼ 2; 3; …P ; p0 r p  1;

iAI

ð7Þ

δi ¼ Di þ 1  Di i A I=f0; Dg

ð8Þ

wi ¼ wi  1 þ δi  1  bi

ð9Þ

bi r C Di r

iAI 

X t0 A P



ð10Þ 1



X

δ0 r

δ0 Z

δD r

δD Z

t0 A P

1





t0 A P

1



t0 A P

1

X



t0 A P

1

t0 A P



1

t¼1

  ðpÞ þ di –t p  1 U λ

pAP

pAP

t¼1

ð13Þ

ð14Þ

   ðpÞ ðt 0 Þ yðDt Þ U θt0 U λ þ t p  dD U λ

pAP

t ¼ t0

ð12Þ

   ðpÞ ðt 0 Þ yð1t Þ U θt0 U λ þ d1 –t p  1 U λ

pAP

t ¼ t0

ð11Þ

   ðpÞ ðt 0 Þ yð1t Þ U θt0 U λ þ d1 –t p  1 U λ

pAP

XNT

  þM U 1  yðDpÞ

ðt 0 Þ

   ðpÞ ðt 0 Þ yði tÞ U θt0 U λ þ di –t p  1 U λ

iA I=f0g;

XNT

  þM U 1  yðDpÞ X

t¼1

Xt 0

  þM U 1  yði pÞ

yði tÞ U θt0 U λ

iA I=f0g;

Xt 0

  þM U 1  yði pÞ X

t¼1

Xt 0

  þM U 1  yði pÞ X



Xt 0

þM U 1  yði pÞ Di Z

iAI

ð15Þ

   ðpÞ ðt 0 Þ yðDt Þ U θt0 U λ þ t p  dD U λ

pAP

ð16Þ

w0 ¼ b0 ¼ 0

ð17Þ

di ; δi ; Di ; wi ; bi ; H i A R þ 8 i A I; yði pÞ A f0; 1g; 8 i A I; p A P

ð18Þ

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

The above non-linear mathematical model is intractable for the practically-sized instances of the train timetabling problem. As a result, an alternative linear model is presented to overcome the difficulty of solving the problem. Although, the nonlinear model is used later in the decoding procedure applied in the neighbourhood search algorithm. The mixed integer linear programming models (MILP) offers a time-indexed definition of the decision variables for scheduling trains with the aim of minimizing the waiting time of passengers at stops. Assume a capacitated train timetabling problem under multi-period travel demand where the inter-arrival time of passengers is constant. In the proposed MILP model, the index t A T is used to display the time intervals that trains are authorized to dispatch. Thus, the length of the planning horizon is divided into equal-sized small time intervals with length α (e.g., 1 min) and trains dispatch only at discrete time points. The arrival rate of passengers at interval t is denoted by λt . For modelling purposes, the binary decision variables xðtÞ is defined where the value of i 1 means the dispatch of i-th train at the beginning of the interval [t,tþ1]. With above notation of the decision variables, the passenger flow variables are defined as follows: bt denotes the number of boarding passengers on the train that depart at the beginning of the interval [t,tþ 1] and wt denotes the number of waiting passengers at the beginning of the interval [t,tþ 1]. The objective can be expressed as a linear function because the floworiented variables are expressed at any time instances. The Eq. (19) represents the total waiting time of passengers assuming the constant inter-arrival time of passenger. The set of constraints to be satisfied are the following: The minimum and maximum allowed headways are stated in Eq. (20). Constraints (21) and (22) impose that each train dispatches at only one time slot and in each time slot only one train allowed to be dispatched, respectively. Eq. (23) defines the relationship between the passenger flow variables. Eq. (24) also states that if a train dispatched at timeslot t, then the number of passengers boarded the train is up to the maximum capacity. that the total passenger  be 1remarked  P It should waiting time is t A T α U wt þ 2λt . Although, the term λt is a constant and it is removed from the objective function. X ½MILP: minz ¼ wt ð19Þ tAT Subject to: X X αðt  1Þ Uxði tþÞ 1 α  α U ðt  1Þ Uxði tÞ r hmax hmin r tAT

tAT

tAT

xði tÞ r1

X

xðtÞ r 1 iAI i

iAI

ð21Þ

t AT

ð22Þ

wt ¼ wt  1 þ λt  1  bt bt r

X

xðtÞ UC iAI i

bt ; wt A R þ

t AT

ð23Þ

t AT

xði tÞ A f0; 1g 8 i A I;

ð24Þ t AT

presented mixed-integer linear mathematical optimization model when the number of train services and time slots increase. Consequently, heuristic and meta-heuristic approaches are more promising in practice. Another point is that the optimal solution of the presented time-indexed formulation gives an upper bound for the main objective function. The value of the optimality gap depends on the level of discretization (amount of α). In other words, decreasing the value of α increases the CPU time exponentially. Thus, a trade-off is needed between the solution quality and the computational effort. 4.2. Mini–max approach In real situations, passengers may leave the train subway system after a certain amount of waiting time and choose another mode of transport. In this situation, setting an upper bound for waiting time is crucial from the passengers' viewpoint. As noticed by Shen and Wilson [38], minimizing the maximum waiting time of passengers over the period of analysis or the number of passengers who cannot board the first incoming train service are alternative measures of service quality. Also Wang et al. [43] argued that the average passenger waiting time is dependent on the number of train services and maximum waiting time. To the best of author's knowledge, the first attempt to achieve the minimization of maximum passenger waiting time was made by Ceder [11]. He introduced the maximum headway as a criterion to provide a passenger level of service. In aforementioned study, a complete timetable was given and the aim is to decrease the number of departures through minimizing the maximum headway time. In this study, a reference schedule with minimum of maximum waiting time will be constructed subject to an upper bound for the accumulative waiting time of passengers. The time-indexed formulation gives us the flexibility of defining an alternative objective function easily. The mini–max approach to train timetabling problem is written in the following mixed-integer programming model:

ð26Þ MiniMax : min z ¼ max wt tAT

Subject to: i A I=f0; Dg ð20Þ

X

443

ð25Þ

The demand-oriented TTP as described above is known to be NP-hard [39] which necessitates applying heuristic or metaheuristic techniques to solve the practical instances of the problem. In other words, the above presented model can be solved up to certain sizes by direct use of commercial optimization solvers. It should be noted that the number of binary decision variables in the presented MILP model depends on the amount of α and the number of train services, usually with integer linear programming models presented in this article more binary variables needs to be defined compared with the MINLP model. Therefore, it is still computationally difficult to find an optimum solution for the

AWT r β Uhmin

ð27Þ

Constraints (20)–(25) remain unchanged. According to the Eq. (27), a desired level of service is maintained through imposing an upper bound for the average waiting time per passengers (AWT). AWT  is defined  P as total waiting time P ðt Þ ðt Þ to total demand ratio t A T α U wt þ 12λ = t A T α U λ . In the best case, the minimum AWT is half of minimum headway (0:5  hmin ). Thus, parameter β Z 0:5 determines the maximum level of service. 4.3. Numerical examples In this section, the impact of changing patterns of demand for travel on the optimal headway is studied through two numerical examples. The first example is adopted from Niu and Zhou [33]. The passenger demands profiles are presented in Fig. 2. In example 2, demand pattern includes two peak points which made it more relevant to the real-world situation. The total demand for travel in example 1 and 2 are 70 and 80 passengers, respectively. Passengers arrive to station with constant interval during the study period [0, 20] min. The length of each demand period is θ ¼2 min. The minimum and maximum headways (hmin ; hmax ) are 1 and 13 min, respectively. In order to obtain average waiting time per passenger (AWT), the total waiting time experienced by all

444

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

Arrival Rate (passenger per minutes)

9 Example 1

Example 2

8 7 6 5 4 3 2 1 0 1

2

3

4

5

6

7

8

9

10

Demand Periods Fig. 2. The demand profiles for example 1 and 2.

Table 2 The computational results of the MINLP and MILP mathematical models in example 1 (C¼ 5). The number of services (D)

5 6 7 8 9 10 11 12 13 14 15 16 17 18 Average

MINLP model

MILP model (α ¼ 1s)

AWT (min) CPU time (s)

AWT (min)

CPU time (s)

4.699 4.157 3.667 3.181 2.728 2.207 1.744 1.309 1.0657 0.981 0.940 0.9019 0.879 0.864

4.699 4.157 3.667 3.181 2.728 2.207 1.744 1.309 1.066 0.981 0.940 0.902 0.881 0.868

31 30 15 15 18 14 17 14 73 990 39640 3736 43875 14469 7352.6

133 464 1104 673 2692 3963 4807 5723 6693 4400 7518 9270 6265 5176 4205.7

Table 3 The computational results of the MINLP and MILP mathematical models in example 2 (C¼5). The number of services (D)

5 6 7 8 9 10 11 12 13 14 15 16 17 18 Average a

 P  passengers during the planning horizon ( i A I wi þ 12δi U H i ) is P ðt Þ divided by the total number of passengers (θ U t A T λ ). The optimal solutions of the MINLP and MILP mathematical models in two numerical examples are reported in Tables 2 and 3, respectively. The MINLP model could solved the instances of the example 1 using KNITRO/GAMS solver and the result indicates that it could find an optimal solutions in a reasonable CPU time. Moreover, the MILP model is solved using CPLEX solver of GAMS. In example 1, the average computation time is nearly 75% more than the MINLP model. This is mainly because of the convexity of the demand function and the fewer binary variables compared with the mixed-integer linear model. However, in example 2, the difference between the average CPU time of the linear and nonlinear models is significant and the state-of-the-art MINLP solvers need extensive computational efforts to find an optimal solution. As shown, the complexity of the MINLP model depends on the convexity of demand function, the number of train services, and the number of demand periods. The result indicates that the shape of the demand profile has as significant effect on the computational effort needed to solve the nonlinear model. On the other

MINLP model

MILP model (α ¼ 1s)

AWT (min) CPU time (s)

AWT (min)

CPU time (s)

4.70 3.948 3.2119 2.55 2.0928 1.72 1.4807 1.2948 1.1474 1.05 1.031a 0.9385a 0.8857 0.86

5.23 4.392 3.303 2.56 2.095 1.73 1.482 1.296 1.148 1.06 0.969 0.918 0.889 0.86

9.5 22 51 261 1152 491 6040 2693 642 269 5349 1926 694 7159 1911.3

104 192 707 637 2453 1199 4711 6602 14546 42430 82440 70986 11839 12911 17982.6

Upper bound.

hand, the CPU time of the linear model is rather independent of the shape of the demand function. The average optimality gaps for the best solutions obtained in examples 1 and 2 using MILP model are 0.05% and 2.33%, respectively. Although, the nonlinear model failed to find optimal solutions for the cases when D¼ 15 and D ¼16 in example 2. In both cases, the MILP could find better solutions compared with the MINLP model. In the investigated small to medium-sized instances of the example 2, the non-convexity of the demand function leads to a significant increase in computational effort. As a conclusion, the results in example 2 show the computational efficiency of the MILP model compared with the non-linear model. This is mostly as a result of the considered non-convex demand function. The influence of the parameter β on the optimality gap   AWTmini  max  AWToptimum =AWToptimum ) of the solutions obtained by the mini–max approach is analyzed and the results are summarized in Tables 4 and 5. The result shows that the problems are solved in a reasonable CPU time while the optimality gaps are highly sensitive to the value of parameter β .

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

Table 4 The computational results of the mini–max approach in example 2 (C¼ 5, D ¼14). β

1.5 1.3 1.1 1.05 1.04

Mini–max approach Gap%

AWT (min)

Max waiting time (min)

CPU time (s)

12.57 14.86 5.90 1.52 1.71

1.182 1.206 1.112 1.066 1.068

9.7 9.7 9.7 9.7 9.7

90 87 82 165 231

Table 5 The computational results of the mini–max approach in example 2 (C¼ 5, D ¼18). β

0.820 0.824 0.825 0.826 0.830 0.850 1

Mini–max approach Gap%

AWT (min)

Max waiting time (min)

CPU time (s)

3.60 1.74 1.51 3.84 3.49 5.81 7.67

0.891 0.875 0.873 0.893 0.890 0.910 0.926

9.7 9.7 9.7 9.7 9.7 9.7 9.7

1144 261 228 1412 2028 3179 33

The presented mathematical formulations are characterized with large numbers of constrains and integer variables for real-sized instances, which make it difficult to obtain optimal solutions in a reasonable time. Also, the outcomes indicate that the CPU time grows exponentially with the number of binary variables. Due to this fact, effective meta-heuristic methods are customized to solve the problem. In the next section, the details of the proposed neighbourhood search algorithms are presented.

5. Neighbourhood search variants The demand-oriented train timetabling is a NP-hard combinatorial optimization problem which is particularly challenging to be solved optimally for practical cases [39]. Consequently, heuristic and meta-heuristic methods are decidedly worthwhile in practice. This section presents different variants of neighbourhood search algorithm to solve train timetabling problem. First, a general summary of the method and the framework on which the neighbourhood search heuristic is applied will be provided. Then, the proposed discrete and real-valued solution encoding schemes are detailed. Then, the basic neighbourhood search method is elaborated in order to reduce the computational time, and the quality of the heuristic solution delivered by the method. For these reasons, adaptive and variable neighbourhood mechanisms are proposed. The proposed neighbourhood search algorithms start with generating random feasible solutions. The initial solutions are evaluated and the best one is selected as a basic solution for neighbourhood search. The main factors in the efficiency and success of neighbourhood search methods are the solution encoding strategies, appropriate definition of the neighbourhood structures, and the efficiency of searching neighbourhood. A neighbourhood search algorithm gradually generates improved sequence of solutions in the neighbourhood of the starting point. Generally speaking, three strategies for choosing a best neighbourhood which include (a) choosing the best improvement, (b) first improvement, and (c) random selection strategies. In the first strategy, all the neighbours are evaluated and the best solution will be chosen. For

445

the large size neighbourhood structures, it may be very time consuming to examine all possible moves to find the best neighbouring solution. In first-improvement strategy, once a better solution is found, the solution is replaced with the original one. This study adopts the first improvement strategy and proposes two different types of solution encoding–decoding procedures. 5.1. Partial frequency-based solution encoding: a discrete scheme The selection of the encoding scheme is of significant importance, as it affects the performance of the algorithm in terms of the quality of solutions and speed of convergence. In this study, the main decision variables are the departure time of trains from the station. An approach to solution encoding is to decide the frequency of train services at each interval. In this discrete type of solution representation, the departure times are still undetermined. Hence, the solution encoding is partial and the decoding step needs a post-optimization procedure. A main advantage of the proposed frequency-based solution representation is that fewer input data are required. In addition, the decoding process can be handled by an exact method or a commercial solver. In this study, the decoding procedure is controlled by the use of MINLP solvers of the GAMS system. Likewise, Niu and Zhang [32] presented an integer nonlinear formulation to determine the number of scheduled trains (frequencies) during each demand period. A genetic algorithm was applied to solve the problem. They assumed that the passenger demands at each period should be supplied by the train services scheduled at the same period. However, no proof of optimality was given. In our implementation, the train frequency is obtained at the first step and the complete timetable will be generated according to the given number of services at the demand intervals. Then, the proposed optimization model is used to decode the partial solution and determine the departure time of train services with respect to the given frequencies at each interval. In frequency based encoding procedure, a solution is represented by an array having size equal to the number of demand periods as shown in Fig. 3. The elements of this array denote the number of train services or frequency (F t ) for the corresponding period (t). Obviously, the sum of frequencies equals to the total number of trains scheduled during the study period: X F ¼D ð28Þ tAT t The relation between the binary decision variables (yðtÞ i ) and the frequency variable (F t ) can be written as: X yðtÞ ¼ F t t A T ð29Þ iAI i

5.1.1. Decoding approach As stated previously, determining the number of trains dispatched within each period does not necessarily mean finding the complete train schedule and it will need to decode the partial solution through solving an integer nonlinear programming model with known frequency at each time interval. Fig. 4 illustrates the proposed MINLP-based decoding–encoding approach employed in discrete solution representation scheme. It should be noted that the basic MINLP model cannot be solved for the large-sized instances due to the complexity of the problem. However, if the frequency variables (F t ) take values, then the resulting problem is a much more easy model to be solved by

F1

F2

F3

.

.

Ft

.

.

Fig. 3. Frequency-based solution encoding.

FNT-1

FNT

446

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

Encoding Discrete representation of service frequency

Decoding using GAMS Solve a MINLP optimization model for each partial solution Fig. 4. The flowchart of the proposed MINLP-based encoding–decoding procedure.

commercial solvers. On the other hand, one of the disadvantages of discrete representation scheme is that the solution is partial and in order to produce the complete solution still an integer nonlinear model should be solved. In order to encode each chromosome, the following model ½MINLPd  should be solved:  X  1 ½MINLPd : min z ¼ w þ δ ð30Þ U Hi i i iAI 2 Subject to: (2)–(17), and (29). di ; δi ; Di ; wi ; bi ; H i ; F t A R þ yðtÞ i A 0; 1

ð31Þ

5.1.2. Initial random feasible solutions We can limit the service frequency at each period according to the minimum and maximum headways. The maximum frequency (F max ) can be achieved only when the trains run with minimum headway. In order to calculate the minimum j k and maximum j k freθ θ and F min ¼ hmax are quencies in each period, Eqs. F max ¼ hmin considered. In each demand period, the number of trains dispatched (F t ) is defined within the minimum and maximum frequencies where the frequencies of the previous periods are known. The following equation determines the feasible range for the frequency of the current period: n n oo Xt  1 r Ft Min F max ; max F min ; D  ðNT  t Þ  F max  F0 t0 ¼ 1 t Xt  1 rD F 0  max ðNT  t Þ  F min ; D t0 ¼ 1 t o Xt  1 0  F max  F t A T =f1; NT g ð32Þ t 0 t ¼1 Likewise, for the first and the last periods, the following equations are derived: Min F max ; max F min ; D  ðNT  1Þ  F max r F 1 r D  max ðNT 1Þ  F min ; D  F max ð33Þ F NT ¼ D 

XNT  1 t0 ¼ 1

F t0

ð34Þ

Random solutions are generated according to the number of trains that have been scheduled in each time period. The feasibility of the randomly generated solutions is controlled through Eqs. (32) and (33) which help the search method to proceed efficiently. It should be noted that the Eqs. (32) and (33) are used in a completely random solution generation method. But, the frequency in periods with higher arrival rate is probably more than the frequency in periods with lower arrival rate. With this consideration, a biased random solution generation is provided in the next section. 5.1.3. Biased random feasible solution generation The generation of the initial solutions is theoretically of slight impact on lengthy runs however its quality can be critical when solutions of good quality need to be found promptly. In order to

speed up the convergence of the algorithm, a randomly biased solution generation procedure is adopted. The simple principle of biased random search is that the frequency of the period with high arrival rate should be more than those established in the periods with low arrival rates. We take the advantages of the discrete representation of the solution which enables the incorporation of travel demand information to the frequency determination procedure at each period. The generation of biased random solutions can increase the possibility of finding the local and global optima. Therefore, in this section, the modified versions of the Eqs. (32) and (33) are presented to generate biased random solutions. This modification is performed by tightening the lower and upper bounds of the train frequencies (F t ) at each period. At first, using the demand information at each period, the normalized arrival rate at period t is defined as follows: pt ¼ P

λðtÞ t0 A T

0

λðt Þ

t AT

ð35Þ

Approximately, we can consider that it is suitable to schedule at

 least min F max ; pt :D number of train services in period t. The term min F max ; ⌈pt :D⌉ is an upper-bound on the maximum number of the scheduled train services in period t. In this regard, the approximate lower and upper bounds for the minimum and maximum frequencies during period t are defined as follows:

 f min ðt Þ ¼ max F min ; min F max ; pt U D t AT ð36Þ f max ðt Þ ¼ max F min ; min F max ; ⌈pt U D⌉

t AT

ð37Þ

The floor and ceiling notations bxc and ⌈x⌉, mean the largest integer less than or equal to x (round down) and means the smallest integer greater than or equal to x (roundup), respectively. As a result, according to what is described in this section, the Eqs. (32) and (33) are modified as follows: n o XNT max f min ð1Þ; D  f ðt 0 Þ r F 1 r D t 0 ¼ 2 max nXNT o  max f ðt 0 Þ; D  f max ð1Þ ð38Þ t 0 ¼ 2 min n o XNT Xt  1 f ðt 0 Þ  F 0 rF t max f min ðt Þ; D  t 0 ¼ t þ 1 max t0 ¼ 1 t nXNT Xt  1 rD F 0  max f ðt 0 Þ; D t0 ¼ 1 t t 0 ¼ t þ 1 min o Xt  1 F 0  f max ðt Þ t A T =f1; NTg  t0 ¼ 1 t F NT ¼ D 

XNT  1 t0 ¼ 1

F t0

ð39Þ ð40Þ

Eqs. (38) and (40) state the feasible range to allocate train frequencies into demand periods. Correspondingly, the feasible range for determining the train frequency at periods t A 2; 3; ::; NT 1 is presented in Eq. (39). In each period, the uniform distribution is used to determine the number of train services according to the above recursive equations. 5.1.4. Neighbourhood structures In this section, a set of neighbourhood structures are proposed for the discrete encoding scheme. The neighbourhood structures include deterministic and probabilistic shift, swap, exchange, and insert mechanisms. The neighbourhood structures are described as follows: Deterministic shift neighbourhood refers to the decreasing the frequency of service in an interval (t 1 ) and increasing it at interval (t 2 ). Moreover, in probabilistic shift neighbourhood the probability of decreasing the service frequency at period (t 1 ) and increasing it at period (t 2 ) depends on the relative arrival rates. For more explanation, consider two time periods t 1 and t 2 in the current solution x. The frequency in these two periods will be changed with probability λðt2 Þ ðt 1 Þ ðt 2 Þ through the probabilistic operators. The ideas underlying the λ

þλ

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

definition of the probabilistic neighbourhood structures are the following: First, it avoids the search procedure from being trapped in local optimum. Second, the purpose is at maintaining the promising properties of the current solution. In deterministic swap, the frequencies of the adjacent periods will be swapped with each other while in probabilistic swap the replacement will be made by the proposed probability. Likewise, in the exchange neighbourhood, a neighbouring solution is obtained by replacing the number of services performed in a period with frequency in another period. Algorithm 1. Neighbourhood search (Pop, MaxIter). begin Generate Pop random feasible solutions; Solve ½MINLPd  to decode the initial solutions xbest ¼ SelectBest(); k: ¼1; while k o ¼ MaxIter do x':¼ NeighbourhoodSearch (xbest); Solve ½MINLPd to decode x' and find the f(x') if f(x') o f(xbest) then xbest: ¼x'; k: ¼k þ1; return (xbest); End. In Algorithm 1, the termination condition is the maximum number of iterations. The algorithm starts with the random solution generated through Eqs. (32) and (33) and then the proposed random biased search runs, using Eqs. (38) and (39), in order to obtain alternative good solutions. 5.1.5. Variable neighbourhood search (VNS) The Variable neighbourhood search (VNS) method combines local search with systematic changes of neighbourhood structures in order to escape from local optima [22]. It has been applied to a wide-range of combinatorial optimization problems including supply chain network design [19] and capacitated vehicle routing problems [49]. However, the applications of the VNS to train timetabling problem was limited. To the best of the authors' knowledge, this is the first article that applies the variable neighbourhood search method for demand-oriented train timetabling problem. In VNS methods, a set of nested neighbourhood structures N k ; 1 r kmin r k r kmax are given. In different variant of the proposed VNS, deterministic and stochastic changes of the neighbourhood structures are performed. For example, Variable Neighbourhood Descent (VND) method changes the neighbourhoods deterministically. Reduced variable neighbourhood search (RVNS) is a stochastic variant of VNS in which solutions from the predefined neighbourhood structures are selected randomly [27]. The pseudo-code of the proposed VND is given in Algorithm 2. The considered neighbourhood structures are in an ordered list consist of the shift, swap, exchange, and insert operators. Algorithm 2. Variable neighbourhood descent VND (P, x, kvnd). 1 Select a set of neighbourhood structures Nk:S - P(S), 1 rk rkvnd; 2 Set stop ¼false; 3 repeat 4 Set k ¼1; 5 repeat 6 x0 ¼ BestImprovement (P, x, Nk(x)); 7 if (f(x0 )of(x)) then 8 x ¼x0 , k ¼1; //Make a move. 9 else k¼ kþ 1; //Next neighbourhood. 10 endif 11 Update stop;

12 13 14

447

until (k ¼ ¼kvnd or stop); until (f(x0 )Zf(x) or stop); return x.

The main steps of the proposed general variable neighbourhood search are given in Algorithm 3. GVNS works with a finite set of pre-specified neighbourhood structures for diversification purposes. The algorithm includes shaking step in order to generate random solution in the k-th neighbourhood of the current solution to be used as an initial solution for VND. At the next step, the produced random solution (x0 ) is improved using a variable neighbourhood descent method, generating a new improved solution (x00 ). Then, the algorithm compares the new solution (x00 ) and the best known solution (x). If an improvement is achieved, the search returns to the first neighbourhood structure (k ¼kmin) and it updates the best known solution (x¼ x00 ); otherwise, the GVNS method attempts to obtain a better solution using a different neighbourhood structure (k¼k þkstep). The GVNS algorithm terminates when it reaches to the maximum number of iterations (Iter4MaxIter) or all the neighbourhood structures are visited (k Zkmax). Algorithm 3. General variable neighbourhood search GVNS (P, x, MaxIter, kmin, kstep, kmax, kvnd). 1 Select a set of neighbourhood structures Nk:S-P(S), kmin rk rkmax; 2 Set stop¼false; 3 Set k ¼kmin; 4 repeat 5 Select x0 ANk(x) at random; //Shaking. 6 x00 ¼VND(P, x0 , kvnd); 7 if (f(x00 ) of(x)) then 8 x ¼x00 , k ¼kmin; //Make a move. 9 else k ¼k þkstep; //Next neighbourhood. 10 Iter¼ Iterþ 1; 11 endif 12 if Iter4 MaxIter, stop¼true; 13 until (k Zkmax or stop); 14 return x.

5.2. Real-valued solution encoding In the proposed real-valued encoding scheme, a candidate solution is represented as a vector of length D, each element corresponding to a train departure time. In this schedule-based solution representation, the feasibility of the solutions only depends on the minimum and maximum headways. In our implementation, the random feasible solutions will be generated using the uniform probability density function taking into account the minimum and maximum allowed headways based on the following recursive equation: max T  ðD  iÞ  hmax ; di  1 þ hmin rdi r min di  1 þhmax ; T  ðD  iÞ:hmin iAI ð41Þ Fig. 5 illustrates the earliest and latest departure times of i-th train with respect to the minimum and maximum headway times. The quality of each candidate solution is assessed by means of the objective function defined in Eq. (1). 5.2.1. Neighbourhood structures The real-valued encoding of the search space is usually suitable for solving optimization problems in which the decision variables takes values in continues domains. Due to the both discrete and continuous nature of the decision variables, a discretization

448

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

Fig. 5. Earliest and latest departure time of trains.

. . . Train i-1

. . . Train i

Step-size

Train i+1

Fig. 6. The illustration of the forward/backward neighbourhood move using different step-sizes.

procedure is suggested. As described at the beginning of this section, in the real-valued encoding scheme, the departure time of each train service is determined. In this encoding scheme, neighbourhood structure is defined as follows: due to the continuity of the decision variables (di ), the continuous nature of the search space is discretized by introducing move steps (step-size). In other words, the step-size parameter is used to determine the number of neighbours. For each train service, the possible moves are performed according to the step-size and the minimum and maximum headways (Fig. 6). This solution coding scheme is proposed to avoid constraints violation during the search procedure. The step-size parameter determines the percentage of the whole neighbourhood scanned at every iteration. Hence, the accuracy of the neighbourhood search algorithm increases by lowering the step-size. On the other hand, in this case, more neighbours are examined and the computational time increases. Here, the objective is to find a trade-off between the obtained quality and the required computational effort. Consider an initial solution (x) in which the departure time of train services are given d1 ; d2 ; …; dD . The earliest and latest departure times of i-th train are denoted by ei and li , respectively, where the departure times of the prior services are fixed. For the first train service, the feasible departure interval is [e1 ; l1 ] where: e1 ¼ max T  ðD  1Þ  hmax ; hmin ; d2  hmax ð42Þ l1 ¼ min T  ðD 1Þ  hmin ; hmax ; d2  hmin

ð43Þ

For the last train service, the earliest and latest departure times are given by eD ¼ dD  1 þ hmin and lD ¼ min di  1 þ hmax ; T , respectively. In the same way, the earliest and latest departure intervals for the other trains are derived with respect to the previously determined departure time: 2 r i rD  1 ei ¼ max T  ðD iÞ  hmax ; di þ 1  hmax ; di  1 þ hmin ð44Þ li ¼ min T  ðD  iÞ  hmin ; di  1 þ hmax ; di þ 1 hmin

2 r i rD 1 ð45Þ

Therefore, in order to search the neighbourhood solutions the j k li  ei . Observably, it is step

number of feasible moves for i-th train is

not practical to search all possible moves and pick the best ones, due to the huge number of possible moves for each candidate solution. Therefore, this study designed an algorithm which allows adaptive search distances causing an improvement in solution quality and search convergence behavior.

5.2.2. Neighbourhood search with adaptive step-size A well-known approach to tuning parameters is by adjusting them throughout the search which called adaptive or self-tuning strategy. Although, the application of the adaptive strategies in the local search methods were limited. Brunato and Battiti [9] developed a novel adaptive random search algorithm for the continuous optimization problems. Alabas-Uslu and Dengiz [3] proposed a self-adaptive mechanism embedded in local search heuristic algorithm for solving the classic vehicle routing problem. AminNaseri and Afshari [5] applied a hybrid genetic algorithm and adaptive local search for process planning and scheduling problem in manufacturing systems. Kim and Liou [26] proposed a novel adaptive local search method for hybrid evolutionary multiobjective algorithms to improve the convergence speed. In our proposed algorithm, adaptive strategies are adopted to achieve a balance between solution quality and convergence speed. In order to escape from local optima, the proposed neighbourhood search algorithm changes the step-size adaptively which is similar to the variable step-size random search algorithm [29,37,48]. The pseudo-code of the proposed neighbourhood search with adaptive step-size is given in Algorithm 4. The algorithm parameters include initial population of the random solutions, the maximum number of iterations of the algorithm (MaxIter), as well as minimum and maximum step-sizes (stepmax and stepmin). The neighbourhood search algorithm starts with generating a set of random solutions (as Pop) then the best obtained solution is used as a starting point for the local search algorithm. In the phase of random solution generation, departure time of train services are determined randomly using uniform distribution di  U½ei ; li . Beginning with an initial step-size value, the neighbourhood search continues with an adapted step-size during the search process. The neighbourhood search algorithm presented here to decide how step-size should be adapted to attain the desired accuracy. Algorithm 4. Neighbourhood Search with Adaptive step-size (Pop, MaxIter, stepmax, stepmin). 1 begin 2 Generate Pop random feasible solutions; 3 xbest ¼ SelectBest(); 4 k:¼1; 5 step:¼stepmax; 6 while ko ¼MaxIter do 7 x': ¼NeighbourhoodSearch (xbest, step); 8 if f(x') of(xbest) then xbest: ¼x'; 9 k:¼ kþ 1; 10 Update step; 11 EndWhile

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

12 return (xbest); 13 End.

stepðkÞ ≔z  stepmax 

In what follows, several suggestions to dynamic modification of the step-size parameter in each iteration (k ¼1, 2,…, MaxIter) of the neighbourhood search algorithm are given. A basic linear adoption mechanism is presented in Eq. (46). In this self-adaptive method, the search starts with the maximum step-size (stepmax) and the step-size gradually is reduced through a linear function of the number of iterations. Similarly, the search may start from a minimum step-size (stepmin ) and decreases gradually with the number of iterations as provided in Eq. (47).   stepmax  stepmin k ð46Þ stepðkÞ ≔ stepmax  MaxIter  stepðkÞ ≔ stepmin þ

stepmax  stepmin MaxIter

 k

ð47Þ

Chatterjee and Siarry [14] questioned the capability to escape from the local optimal with linear functions and they proposed a nonlinear time-varying adapting mechanism as follows:  n   MaxIter k ð48Þ  stepmax  stepmin stepðkÞ ≔ stepmax  MaxIter  stepðkÞ ≔ stepmin þ

MaxIter k MaxIter

n

   stepmax  stepmin

ð49Þ

With n ¼ 1, the above formula becomes linearly adaptive mechanisms. The search starts with a high (or low) step-size which allows the diverse (or intense) exploration of the search space and then decreases (or increase) the step-size progressively. Stochastic adaptive search is an alternative approach which behaves between pure random and pure adaptive search methods [52]. In a pure random adoption of the step-size, the uniform distribution is employed as follows:   ð50Þ stepðkÞ ≔ uniform stepmin ; stepmax Linear time-varying with random changes area also suggested by Feng et al. [20]:     MaxIter k  stepmax  stepmin ð51Þ stepðkÞ ≔z  stepmin þ MaxIter

    MaxIter  k  stepmax  stepmin MaxIter

449

ð52Þ

Another efficient way of escaping from the local optima is to change the step-size using the search history. Precisely, if the neighbourhood search algorithm fails to improve the upper bound after a specified number of consecutive iterations (L), the step-size return to the maximum (stepmax ) or minimum values (stepmin ). In such a self-tuning mechanism, different linear or nonlinear stepsize functions can be utilized. We suggest the nonlinear increasing time-varying adapting mechanisms. In what follows, the selfadaption parameter tuning procedures are compared through numerical examples. The initial solution used in all algorithms is an even-headway schedule. Fig. 7 illustrates the convergence graph of the adaptive stepsize neighbourhood search algorithm assuming stepmax ¼ 1 min and stepmin ¼ 1 s. The outcomes indicate that the pure random adaptive approach has the highest convergence rate as a result of its perturbation mechanism which helps to escape from the local optima. Due to the small sizes of the instances, it could reach to near-optimal solutions in a very short time. Although, the pure random adoption may fails to find global solution for large-scale instances in a reasonable time. The gradually linear increasing and non-linear decreasing (n ¼2, L¼5) outperforms the other complementary adaptive methods. Fig. 8 illustrates the comparison between the pure-random adoption and history-enabled mechanisms using convergence graph. As can be seen from the graph, the search efficiency improves through increasing the number of consecutive non-improving iterations (L). Although, L¼ 5 is a good choice in terms of solution quality and convergence criteria. Fig. 9 shows the convergence graph of the pure-random adoption mechanism with different values of the maximum stepsize. The results suggest increasing the value of the maximum step-size in order to improve the convergence rate. With stepmax ¼ 2:5 the neighbourhood search algorithm could find the optimal solution. As a conclusion, adaptive parameter tuning strategies, as shown in this section, could significantly improve the performance of the basic neighbourhood search algorithm.

6. Application to real cases Tehran urban and suburban railway is a rapid system transit rail with five operating lines serving Tehran. Currently, the Tehran

Fig. 7. Convergence graph of the pure-random adoption (D ¼10).

450

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

Fig. 8. Convergence graph of the pure-random adoption and history-enabled mechanisms (D ¼10).

Fig. 9. Convergence graph of the pure-random adoption mechanism with different step-max (D¼ 10).

Fig. 10. Line 5 of Tehran metro.

metro network is 122 km long and on average two million passengers use it daily. Tehran–Karaj subway line is investigated that connect two major stations (Golshahr and Tehran). It is a 41.5 km long route starting from Tehran terminus and ending at Golshahr. In public transportation systems, the travel demand is characterized by peak and off-peak fluctuations. In the investigated case, the transport demand is characterized by arrival rate of passengers. Different scenarios of passenger arrival rate at Tehran and Golshahr stations are shown in Figs. 11 and 12, respectively. Each one of these demand scenarios represents the arrival rate of the passengers in a specific month of the year according to the collected historical data. The following real-world examples illustrate the analysis of the mathematical models derived above and demonstrate the advantages and disadvantages of both the linear and the nonlinear

Fig. 11. Demand profiles at Tehran station.

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

models. Consider two main stations of the line 5 of Tehran Underground Metro in Fig. 10. The operational information of the Line 5 (Tehran Metro) is given (Table 6). The baseline daily headways used in terminal stations are illustrated in Fig. 13. It can be observed that the headway reaches to minimum value (hmin ) during the peak hours, which indicates that the train service plan fits the passenger flow demand relatively. The baseline headways were obtained from the rail experts that are produced manually with trial and error procedure.

7. Results and discussion For a better comparison of the solution quality, a simple lower bound of the average waiting time per passenger is obtained. The basic idea is designing a regular timetable. Welding [47] estimated the expected passenger waiting time AWT where both the vehicle's and passenger's arrivals at stop are random. The proposed formulation comprises of the average headway H and the

Fig. 12. Demand profiles at Golshahr station.

451

headway variance σ 2 h :   1 σ2 h AWT ¼ Hþ 2 H

ð53Þ

From the above equation, it can be concluded that the average passengers waiting time is at minimum value when the headway variance is zero and vehicles run with minimum headway. On the other hand, the average headway H can be regarded as a scheduled constant headway [2]. Therefore, the lower bound for of the average waiting time is half of minimum headway (0:5  hmin ). This lower bound is used to assess the optimally gap of the best solutions obtained from optimization methods. The performance of the adaptive neighbourhood search and VNS meta-heuristic are assessed on real instances with different demand scenarios. All developed mathematical formulations and algorithms were written in GAMS modelling language. The computational results of the MILP mathematical model are given in Table 7. This table provides the best obtained solutions using MILP model (where α ¼ 3) in term of average waiting cost and CPU time required by GAMS/ CPLEX solver. The column "Gap" indicates the gap (upper bound  lower bound)/lower bound. This gap is calculated according to the proposed lower bound. The column "Baseline" specifies the percentage of reduction in AWT compared with the baseline timetable. The results indicate that the solution quality of CPLEX was rather poor, with a maximum average relative MILP gap of 52.64% and 45.52% at Tehran and Golshahr stations, respectively. However, in all of the solutions, the MILP model found the better timetables with fewer waiting cost (5.14% on average) compared with the baseline schedule. The average value of the total CPU time required by CPLEX to solve six instances of the problem at Tehran and Golshahr stations are 2.94 and 1.77 h, respectively. Table 8 provides the comparison between the adaptive and variable neighbourhood search algorithms in both terms of the solution quality and computational time. For the result reported in Table 8, the ANS and VNS parameters were set as presented in Tables 9 and 10, respectively. The column "Imp" shows the percentage reduction of the upper bounds obtained from the neighbourhood search algorithms over GAMS/CPLEX's performance which is calculated through equations

Table 6 The operational information of the Line 5 (Tehran Metro). Station name

C (passengers)

θ (min)

NTðperiodsÞ

hmin ðminÞ

hmax ðminÞ

D (train services)

Tehran Golshahr

2000 2000

60 60

19 19

7 7

25 20

97 99

Fig. 13. Baseline headway time between trains at terminal stations of the Tehran–Karaj metro line.

452

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

ðUBGVNS  UBGAMS Þ=UBGAMS % and ðUBANS  UBGAMS Þ=UBGAMS %. Furthermore, column "Baseline" indicates the percentage reduction of the upper bounds obtained from the neighbourhood search algorithms compared with the baseline schedules which is calculated through equations ðUBGVNS  AWTBaseline Þ=AWTBaseline % and ðUBANS  AWTBaseline Þ=AWTBaseline %. The average reduction of the AWT resulted from the GVNS and ANS algorithms in all demand scenarios are 2.11% and 1.43%, respectively, compared with the GAMS/CPLEX's performance. The average run time for the GVNS was much smaller compared to ANS and CPLEX solver (Table 8). The average value of the total CPU time required by GVNS to solve all instances of the problem is on average 1.15 h. The results demonstrate that the proposed GVNS method provides significant reduction in computational effort (50.97% on average) compared with the GAMS/CPLEX solver. Computational results comparing the performances of GVNS and ANS algorithms with the CPLEX solvers demonstrate the effectiveness of the frequency-based solution encoding combined with a MINLP-based decoding procedure. Moreover, due to the huge size of the solution space, the realvalued coding approach was not suitable for the large-sized instances. This was due to the fact that the proposed frequencybased encoding–decoding scheme provides a better search property during the optimization compared with the real-valued encoding method. The results obtained are also promising enough to suggest that using biased random search is an effective pre-optimization step.

Baseline schedule

MIP model (α ¼ 3)

In railway networks with high rate of congestion, the minimization of the passenger waiting time becomes more and more important. In order to gain higher efficiency and better utilization of the infrastructure and fleet capacity in railway traffic systems, it is required to use optimization approaches. In this study, the focus was on the comfort of the passengers through presenting effective solution methods for minimizing the total waiting time in metro system. Firstly, a mixed integer nonlinear formulation was presented. A linear mixed-integer program was extended based on discretization of the planning horizon. Mixed-integer linear model could find an optimal solution for almost all of the small- and medium-sized instances. From the result of the numerical examples, a conclusion can be drawn that the convexity of the demand function is an important factor that influence the solution time of the MINLP model. Due to the non-convexity of the real-world demand functions, the MINLP models are intractable and linear models are alternative options. Although, it was still difficult to solve the real-world instances by the direct use of the commercial optimization packages. In order to cope with this problem, two versions of the neighbourhood search algorithms (including the adaptive and variable neighbourhood search) were employed to find the near-optimal solutions that minimizes total waiting time of passengers. The adaptive neighbourhood search algorithms utilized the complete real-valued encoding scheme while the variable neighbourhood search method employed the partial frequency-based solution representation where the decoding was made using the mixed-integer nonlinear solvers of the GAMS. The proposed novel solution encoding–decoding procedures were tested on real-sized examples of the Tehran metropolitan network. The outcomes indicated that the proposed meta-heuristic

AWT

UB

Table 9 The parameters of the ANS algorithm.

Table 7 The computational results of the MILP mathematical model. Station name

Tehran

Golshahr

Demand scenario

1 2 3 4 5 6 1 2 3 4 5 6

8. Conclusion

5.533 5.561 5.559 5.595 5.617 5.587 5.571 5.404 5.390 5.410 5.357 5.412

5.301 5.343 5.377 5.394 5.315 5.325 5.318 5.041 5.046 5.058 5.047 5.050

Gap (%) Baseline (%) CPU time (h) 51.46 52.66 53.63 54.11 51.86 52.14 51.94 44.03 44.17 44.51 44.20 44.29

 4.19  3.92  3.27  3.59  5.38  4.69  4.54  6.72  6.38  6.51  5.79  6.69

9.65 1.15 1.32 1.04 0.52 3.98 3.02 1.72 2.01 1.95 0.47 1.50

MaxIter

stepmax (min)

stepmin (min)

Pop (solutions)

4000

5

0.5

200

Table 10 The parameters of the GVNS algorithm. MaxIter

kmin

kmax

kstep

kvnd

3000

1

4

1

2

Table 8 The computational results of the neighbourhood search algorithms. Station name

Tehran

Golshahr

Demand scenario

1 2 3 4 5 6 1 2 3 4 5 6

Adaptive step-size neighbourhood search (ANS)

General variable neighbourhood search algorithm (GVNS)

Upper bound (UB)

CPU time (h)

Baseline (%)

Imp (%)

Upper bound (UB)

CPU time (h)

Baseline (%)

Imp (%)

5.271 5.263 5.344 5.337 5.270 5.319 5.230 5.006 4.879 5.033 4.779 5.002

1.65 1.41 2.74 1.61 1.74 1.34 1.47 1.37 1.35 1.3 1.53 1.46

 4.74  5.36  3.87  4.61  6.18  4.80  6.12  7.36  9.48  6.97  10.79  7.58

 0.57  1.50  0.61  1.06  0.85  0.11  1.65  0.69  3.31  0.49  5.31  0.95

5.245 5.259 5.297 5.303 5.263 5.268 5.228 4.936 4.879 4.918 4.766 4.948

3.01 2.25 1.15 0.12 0.08 0.37 0.11 1.42 1.26 1.09 1.27 1.76

 5.21  5.43  4.71  5.22  6.30  5.71  6.16  8.66  9.48  9.09  11.03  8.57

 1.06  1.57  1.49  1.69  0.98  1.07  1.69  2.08  3.31  2.77  5.57  2.02

E. Hassannayebi, S.H. Zegordi / Computers & Operations Research 78 (2017) 439–453

methods could find the near-optimal timetables in a reasonable amount of time. The disadvantage of the complete solution representation was the moderately extensive computational effort needed to find the near-optimal solution. However, the proposed frequency-based solution representation could efficiently search the solution space aimed to compact the solution information. An important topic for further research is to combine the proposed real-valued and discrete solution encoding methods. Such hybrid algorithms may help to obtain better upper bounds. It can take the advantages of both discrete and real-valued encoding procedures. The hybrid algorithm can start with the solutions that are generated by the frequency-based encoding structure. The best solution obtained can be regarded as an initial solution for the adaptive neighbourhood search algorithm.

References [1] Abril M, Salido MA, Barber F. Distributed search in railway scheduling problems. Eng Appl Artif Intell 2008;21:744–55. [2] Adebisi O. A mathematical model for headway variance of fixed-route buses. Transp Res B Methodol 1986;20:59–70. [3] Alabas-Uslu C, Dengiz B. A self-adaptive local search algorithm for the classical vehicle routing problem. Expert Syst Appl 2011;38:8990–8. [4] Albrecht T. Automated timetable design for demand-oriented service on suburban railways. Public Transp 2009;1:5–20. [5] Amin-Naseri M, Afshari AJ. A hybrid genetic algorithm for integrated process planning and scheduling problem with precedence constraints. Int J Adv Manuf Technol 2012;59:273–87. [6] Amin- Naseri, Baradaran Vahid. Accurate estimation of average waiting time in public transportation systems. Transp Sci 2014;49(2):213–22. [7] Barrena E, Canca D, Coelho LC, Laporte G. Exact formulations and algorithm for the train timetabling problem with dynamic demand. Comput Oper Res 2014;44:66–74. [8] Barrena E, Canca D, Coelho LC, Laporte G. Single-line rail rapid transit timetabling under dynamic passenger demand. Transp Res B Methodol 2014;70:134–50. [9] Brunato M, Battiti R. Rash: a self-adaptive random search method. Adaptive and Multilevel Metaheuristics. Springer; 2008;136:95–117. [10] Canca D, Barrena E, Algaba E, Zarzo A. Design and analysis of demand-adapted railway timetables. arXiv preprint arXiv 2013;1307:0970. [11] Ceder A. A procedure to adjust transit trip departure times through minimizing the maximum headway. Comput Oper Res 1991;18:417–31. [12] Ceder A, Philibert L. Transit Timetables Resulting in Even Maximum Load on Individual Vehicles. 2014. [13] Ceder AA, Hassold S, Dano B. Approaching even-load and even-headway transit timetables using different bus sizes. Public Transp 2013;5:193–217. [14] Chatterjee A, Siarry P. Nonlinear inertia weight variation for dynamic adaptation in particle swarm optimization. Comput Oper Res 2006;33:859–71. [15] Chierici A, Cordone R, Maja R. The demand-dependent optimization of regular train timetables. Electron Notes Discret Math 2004;17:99–104. [16] Cucala A, Fernández A, Sicre C, Domínguez M. Fuzzy optimal schedule of high speed train operation to minimize energy consumption with uncertain delays and driver's behavioral response. Eng Appl Artif Intell 2012;25:1548–57. [17] Cury J, Gomide F, Mendes M. A methodology for generation of optimal schedules for an underground railway system. Autom Control IEEE Trans 1980;25:217–22. [18] Domínguez M, Fernández-Cardador A, Cucala AP, Gonsalves T, Fernández A. Multi objective particle swarm optimization algorithm for the design of efficient ATO speed profiles in metro lines. Eng Appl Artif Intell 2014;29:43–53. [19] Eskandarpour M, Nikbakhsh E, Zegordi SH. Variable neighbourhood search for the bi-objective post-sales network design problem: A fitness landscape analysis approach. Comput Oper Res 2014;52:300–14. [20] Feng, Y, Teng, G-F, Wang, A-X, Yao, Y-M. Chaotic inertia weight in particle swarm optimization. In: Proceedings of the IEEE second international conference on innovative computing, information and control, ICICIC'07. 2007. p. 475. [21] Furth PG, Muller TH. Part 4: capacity and quality of service: Service reliability and hidden waiting time: insights from Automatic Vehicle Location data. Transp Res Rec J Transp Res Board 2006;1955:79–87. [22] Hansen P, Mladenović N, Perez JAM. Variable neighbourhood search: methods and applications. Ann Oper Res 2010;175:367–407. [23] Hassannayebi E, Sajedinejad A, Mardani S. Urban rail transit planning using a two-stage simulation-based optimization approach. Simul Model Pract Theory 2014;49:151–66.

453

[24] Islam MK, Vandebona U, Dixit VV, Sharma, A. A bulk queue model for the evaluation of impact of headway variations and passenger waiting behavior on public transit performance. 2014. [25] Jamili A, Shafia MA, Sadjadi SJ, Tavakkoli-Moghaddam R. Solving a periodic single-track train timetabling problem by an efficient hybrid algorithm. Eng Appl Artif Intell 2012;25:793–800. [26] Kim H, Liou M-S. Adaptive directional local search strategy for hybrid evolutionary multiobjective optimization. Appl Soft Comput 2014;19:290–311. [27] Lazic J. New variants of variable neighbourhood search for 0–1 mixed integer programming and clustering. Brunel University, School of Information Systems, Computing and Mathematics; 2010. [28] Mao H, Li Z. Train schedule adjustment strategies for train dispatch. TELKOMNIKA Indones J Electr Eng 2013;11:2526–34. [29] Masri S, Bekey G, Safford F. A global optimization algorithm using adaptive random search. Appl Math Comput 1980;7:353–75. [30] Mcleod F. Estimating bus passenger waiting times from incomplete bus arrivals data. J Oper Res Soc 2007;58:1518–25. [31] Niu H, Tian X. An approach to optimize the departure times of transit vehicles with strict capacity constraints. Math Probl Eng 2013:2013. [32] Niu Huimin, Zhang Minghui. An optimization to schedule train operations with phase-regular framework for intercity rail lines, Discret Dyn Nat Soc, 2012, (Article ID 549374, 13 pp.), http://dx.doi.org/10.1155/2012/549374. [33] Niu H, Zhou X. Optimizing urban rail timetable under time-dependent demand and oversaturated conditions. Transp Res C Emerg Technol 2013;36:212–30. [34] Osuna E, Newell G. Control strategies for an idealized public transportation system. Transp Sci 1972;6:52–72. [35] Peña-Alcaraz M, Fernández A, Cucala AP, Ramos A, Pecharromán RR. Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy. Proc Inst Mech Eng F J Rail Rapid Transit 2012;226:397–408. [36] Sato K, Tamura K, Tomii N. A MIP-based timetable rescheduling formulation and algorithm minimizing further inconvenience to passengers. J Rail Transp Plan Manag 2013;3:38–53. [37] Schumer M, Steiglitz K. Adaptive step size random search. Autom Control IEEE Trans 1968;13:270–6. [38] Computer-aided scheduling of public transport. Vol. 505 of the series Lecture notes in economics and mathematical systems p. 335–63. [39] Sun L, Jin JG, Lee D-H, Axhausen KW, Erath A. Demand-driven timetable design for metro services. Transp Res C Emerg Technol 2014;46:284–99. [40] Sun X, Zhang S, Dong H, Chen Y, Zhu H. Optimization of metro train schedules with a dwell time model using the lagrangian duality theory. Intell Transp Syst IEEE Trans 2014;16:1–9. [41] Sun X, Zhang S, Dong H, Zhu H. Optimal train schedule with headway and passenger flow dynamic models. In: Proceedings of the 2013 IEEE International Conference on Intelligent Rail Transportation (ICIRT). 2013. p. 307–12. [42] Tong L, Zhou L, Tang J. Frequency Setting of Urban Mass Transit Based on Passenger Flow Along the Entire Line. In: Proceedings of the 2nd International Symposium on Rail Transit Comprehensive Development (ISRTCD). Springer: 2014. p. 59–65. [43] Wang J, Zhai G, Tan Z. Lost Time of Passenger Transfer and Waiting for Trains. In: Proceedings of third international conference on transportation engineering (ICTE). 2011. [44] Wang Y, de Schutter B, Delft C, van den Boom T, Ning B, Tang T. OriginDestination Dependent Train Scheduling Problem with Stop-Skipping for Urban Rail Transit Systems. In: Proceedings of the transportation research board 93rd annual meeting. 2014a. [45] Wang Y, de Schutter B, van den Boom T, Ning B, Tang T. Real-time scheduling for single lines in urban rail transit systems. In: Proceedings of the IEEE international conference on intelligent rail transportation (ICIRT). 2013. p. 1–6. [46] Wang Y, de Schutter B, van den Boom TJJ, Ning B, Tang T. Efficient bilevel approach for urban rail transit operation with stop-skipping. Intell Transp Syst IEEE Trans 2014;15:1–13. [47] Welding P. The instability of a close-interval service. Oper Res 1957;8:133–42. [48] White L, Day R. An evaluation of adaptive step-size random search. Autom Control IEEE Trans 1971;16:475–8. [49] Xiao Y, Zhao Q, Kaku I, Mladenovic N. Variable neighbourhood simulated annealing algorithm for capacitated vehicle routing problems. Eng Optim 2014;46:562–79. [50] Yalcinkaya Ö, Mirac Bayhan G. Modelling and optimization of average travel time for a metro line by simulation and response surface methodology. Eur J Oper Res 2009;196:225–33. [51] Yao X, Zhao P, Qiao K. Simulation and evaluation of urban rail transit network based on multi-agent approach. J Ind Eng Manag 2013;6:367–79. [52] Zabinsky Zelda B. Stochastic adaptive search for global optimization, vol. 72. Springer Science & Business Media; 2013. [53] Zhao S, Cao Y, Tian Q. The real-time dispatch methods for urban rail transportation. In: Proceedings of IEEE 2011 international conference on transportation, mechanical, and electrical engineering (TMEE). 2011. p. 12–15.