Train schedule optimization in a high-speed railway system using a hybrid simulation and meta-model approach

Computers & Industrial Engineering 138 (2019) 106110

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Train schedule optimization in a high-speed railway system using a hybrid simulation and meta-model approach

T

Erfan Hassannayebia, , Morteza Borounb, Shafagh Alaei Jordehic, Hamrah Kord ⁎

a

Industrial Engineering Department, Islamic Azad University, Central Tehran Branch, Tehran, Iran School of Industrial and Systems Engineering, College of Engineering, University of Tehran, Iran c School of Railway Engineering, Iran University of Science and Technology, Iran d Department of Industrial Engineering, College of Engineering, University of Tehran, Iran b

ARTICLE INFO

ABSTRACT

Keywords: High-speed railway Waiting time Response surface methodology Simulation-optimization Meta-model

In urban metro systems, several factors such as passenger congestion may affect the movement of trains, resulting in significant deviations from the planned timetable. This study proposes a meta-model simulation-based optimization approach for minimizing the passenger wait time under random disturbances. The modeling framework is based on a hybrid discrete-event simulation and response surface methodology (RSM), thereby capturing the benefits of both methods. The verification process is handled using the calculation of second-order derivate along with an Analysis of variance (ANOVA) test. Simulation meta-data that is most consistent with the data collected from simulation experiments are used to describe the relationship between response variables, i.e., average waiting time per passenger, and input variable, i.e., time intervals between two departures. The optimal parameters of the RSM are decided using a factorial design. The validity of the results is confirmed using the real instances of a high-speed rail transit line.

1. Introduction Rail transportation organizations seek for optimized schedules to minimize passenger dissatisfaction and improve timetable reliability (Shakibayifar, Hassannayebi, Mirzahossein, Taghikhah, & Jafarpur, 2019). During railway planning, an inconsistency is present between designers who emphasize the train operations and passengers who look for efficient and flexible services with minimum waiting time (Parbo, Nielsen, & Prato, 2016). The design procedure of a train timetable involves different constraints, e.g., capacity limitations, safety headway, periodic or non-periodic timetables, operator cost, overtaking rules, service quality and user costs (Kang et al., 2016). In high-speed rail transit services, it is crucial to reduce the waiting time for passengers by modifying the timetable robustness, which has the potential for absorbing the small disturbances and does not reduce the system performance under unexpected conditions (Shakibayifar et al., 2017b). Fundamentally, passenger waiting time is as an essential indicator of service quality in urban rail networks (Hassannayebi, Sajedinejad, & Mardani, 2016a). In high-speed rail systems, several causes may affect the traverse of trains, resulting in instability of the services (Xu, Li, Liu, Ran, & Qin, 2019). For example, an unexpected increase in travel demand that is causing crowds or frequent travel disturbances, often

⁎

occurring accidentally (Hassannayebi, Sajedinejad, & Mardani, 2016b). The occurrence of even small disturbances may result in prolongation of train dwell time at stops and make the train headway irregular (Lai & Leung, 2018). Therefore, it is vital to design an optimal train timetable by modifying train headways under random disturbances (Shakibayifar et al., 2017a). Due to complex and stochastic nature of the transportation planning problem, several studies addressed the use of hybrid simulation-based optimizations (SBO) techniques (Hassannayebi, Sajedinejad, & Mardani, 2014, Sajedinejad, Mardani, Hasannayebi, Mir Mohammadi Kabirian, 2011). The literature of SBO is extensive with a wide range of methods, e.g., discrete SBO, meta-heuristics, surrogate approaches, and meta-models (Antal et al., 2019). Among different simulation-optimization approaches, RSM is one of the most popular statistical methods in the field of quasi meta-model approaches. RSM uses the fitted simulation outputs to build an approximate regression function, which mathematically expresses the relationship between a set of independent input variables and the response variable (Dengiz & Belgin, 2014). However, the number of applications of meta-models to the design problem of the service industry is increasing because of its potential to decrease the number of required simulation replications for solution evaluation (Korytkowski, Wiśniewski, & Rymaszewski, 2013). In this

Corresponding author. E-mail addresses: [email protected] (E. Hassannayebi), [email protected] (M. Boroun), [email protected] (H. Kor).

https://doi.org/10.1016/j.cie.2019.106110 Received 18 July 2019; Accepted 1 October 2019 Available online 09 October 2019 0360-8352/ © 2019 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

paper, a train headway optimization problem is studied to illustrate the applicability of a meta-model based simulation-optimization framework. This research is motivated by a need to optimize the utilization of fleet and minimization of the operating cost by adjustment of train headways under random disturbances. In this study, the ultimate goal is to find an optimal headway time for train services running on both the inbound and outbound routes while minimizing the average waiting time for passengers concerning the capacity constraints. The contributions of the present study are twofold: First, a hybrid meta-model approach using Design of Experiments (DOE) is proposed to optimize the level of service for passengers. The proposed framework enables to build a meta-model and the decide the optimal values of the input variables, i.e., headway times, with fewer computational cost as against the pure simulation considering. Second, the present study utilizes RSM and leverages the flexibilities that RSM provides. The first and secondorder models are constructed using the Box–Beckhen design method. Several statistical tests, along with stationary point analysis, are also provided to validate the result of the simulation model. The hybridization of the event-driven simulation and RSM algorithm reduces the number of required simulation runs. The simulation model is first used to develop a meta-model and in the next step to optimize the response variable. RSM is used to estimate the pattern of response variations at a specific point with low variations in the input variable, which results in the approximation of the optimal input variables. The structure of the paper is organized as follows: next section provides a literature review of train timetabling problem in urban railway with the most emphasis on simulation-optimization applications. Section 3 describes the developed event-driven simulation model. Section 4 discusses the response surface methodology used for the simulation-optimization method. The results of the proposed SBO, as well as the experimental outcomes, are discussed in Section 5, which followed by concluding remarks in Section 6.

(Earliness-tardiness) and regularity costs. Each sub-problem related to routing a train solved optimality as an instance of the blocking, no-wait job-shop scheduling problem by branch and bound (B&B) algorithm. As a real application, the computational efficiency of the proposed approach has been proved on test instances of Milan metro. Researchers have studied the effects of uncertainty on optimizing service quality in an urban rail system. Several techniques and solution methods have been used in the literature. For example, Yalcınkaya and Bayhan (2009) applied a meta-model approach to optimize the travel time of passenger in an urban metro system. They assumed the headway as an input variable while average travel time and the train capacity utilization are regarded as response variables. The optimization process was handled by a Derringer–Suich multi-response method to decide the optimal values of the input factors. Dong, Ning, Cai, and Hou (2010) addressed the real-time train scheduling problem in an urban metro. They applied simulation model along with a dynamic control model for optimizing the train service under capacity constraints. Canca, Barrena, Zarzo, Ortega, and Algaba (2012) introduced optimal train reallocation strategies and scheduling shuttle trains under service disruptions on a double direction line to minimize the demand congestion and maintain a level of quality of service. They proposed a mixed integer programming formulation to determine the optimal specification of short-turning and deadheading strategies with the capability of generating both periodic and non-periodic train timetables. Canca et al. (2012) studied the problem of inserting special train services; namely, shuttles, to distribute transportation supply efficiently. To achieve this objective, an optimization model is proposed by incorporating deadheading and short-turning strategies, attaining lower headway at crowded locations with minimum changes imposed on the planned services. As output, the detailed timetable of shuttles includes the departure and arrival times, dwell times and speeds are achieved by solving the proposed optimization model. Rahimi Mazrae Shahi, Fallah Mehdipour, & Amiri, (2016) developed a simulation-based optimization model using response surface methodology for the design of the train schedule for subway systems. A multinomial meta-model was employed to fit the simulation data. The optimization process was carried out by a sequential quadratic programming (SQP) and a weighed metric technique. Azad, Hassini, and Verma (2016) proposed an optimization model for disruption risk management of rail networks. The designed method was implemented on real-world instances of the US rail network. Pouryousef, Lautala, and Watkins (2016) proposed a multi-objective linear formulation of train rescheduling for rail networks. The optimization model has been combined with a rail simulation tool to improve the capacity of the rail system. The model is capable of both conflict resolution and schedule compression. The validation procedure was performed by using RAILSYS simulation software. Yin, Yang, Tang, Gao, and Ran (2017) proposed an alternative formulation for passenger-oriented train service planning in the metro system. The mathematical model was solved by Lagrangian relaxation heuristic to find an optimized timetable regarding the passenger wait time and energy consumption. Liu et al. (2018b) proposed approximate dynamic programming (DP) method to minimize energy consumption in a subway system under time-dependent demand assumption. The modeling framework allows finding a tradeoff between energy consumption, capacity utilization, and passenger waiting time. Kim et al. (2019) designed a data-based framework for optimizing urban train timetables. The method involves the use of automated fare collection (AFC) data to minimize the total sum of passenger riding, waiting, and walking times. Mo, Yang, Wang, and Qi (2019) presented a mixed-integer nonlinear programming model for generating an energy-efficient train timetable under capacity constraints. The method was used to improve the service quality by explicitly modeling of passenger flows. A Tabu search algorithm solved several instances of Beijing Yizhuang subway line. Li et al. (2019a) aim to find a trade-off between effectiveness and equity metrics in high-speed urban rail system considering

2. Literature review The previous works on metro train scheduling problem are comprehensive, ranging from train timetabling (Shi, Yang, Yang, & Gao, 2018), energy-efficient scheduling (Liu et al., 2018a, Ning, Xun, Gao, & Zhang, 2014), and headway optimization (Hassannayebi & Zegordi, 2017, Hassannayebi, Zegordi, Amin-Naseri, & Yaghini, 2018) to transfer synchronization (Wong, Yuen, Fung, & Leung, 2008, Wu et al., 2015). Among early works, Bookbinder and Desilets (1992) proposed a mathematical optimization model for scheduling the transfers in a rail network. Eberlein, Wilson, Barnhart, and Bernstein (1998) studied the real-time optimization of train timetable in an urban transit line using a nonlinear integer programming model. It was assumed that the vehicle dwell times, speeds, and passenger arrival rates are constant. Eberlein, Wilson, and Bernstein (1999) presented mathematical models and a heuristic algorithm to solve train control problems. They considered the combination of skip-stop and holding policies to provide robust and efficient train timetable. Vázquez-Abad and Zubieta (2005) designed an adaptive and dynamic simulation model to optimize the service level offered in a complex urban transit system. The simulation model accounts for the implicit flow of passengers while attempts to adjust the headway times on each rail line. Flamini and Pacciarelli (2008) addressed a real-time routing and sequencing problem of a metro terminal to optimize punctuality or tardiness/earliness and regularity concerning deviation from planned headway. The problem is formulated as a bi-criteria blocking job shop scheduling problem by using the alternative graph model of Mascis and Pacciarelli (2002) with particular constraints. A successful application of optimization methods for real-time traffic control in metro stations is presented by Mannino and Mascis (2009). They developed a real-time routing and scheduling model to operate trains in metro stations minimizing a weighted sum of punctuality 2

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

the dynamic demand situation. The underlying problem is to improve the fairness of the generated train timetables regarding the wait time of different passengers. The problem was solved using a simulated annealing (SA) algorithm and a neighborhood search method. Likewise, Li et al. (2019b) presented an MINLP model to minimize the total cost of passenger wait time and operating costs in an urban rail transit system. The decision model accounts to optimize the short turning services and train headways simultaneously. Several instances of Shanghai Metro network were examined by a two-phase genetic algorithm (GA). Högdahl, Bohlin, and Fröidh (2019) addressed the train delay prediction and minimization problem in a railway system using a hybrid simulation-optimization method. They proposed a bi-objective mixedinteger programming formulation that aims to minimize the weighted sum of approximate delay time and train travel time. To the best of our knowledge, very limited studies were devoted to the meta-model based simulation-optimization approaches to train scheduling problem in high-speed railway. Our work is different from Hassannayebi et al. (2014) in that it develops a meta-model based approach to improve the efficiency of the SBO. Also, this research differs from the existing studies, i.e. (Rahimi Mazrae Shahi, Fallah Mehdipour, & Amiri, 2016, Yalcınkaya and Bayhan, 2009), in that it takes particular constraints into consideration in optimal planning the train headway times. Our research work also considers the second-order analysis of stationary point processes. The next following section presents the simulation modeling framework and the response surface methodology used to find the optimal levels of the input factors.

modules, i.e., a decision-making module, a scenario planning module, an solution assessment module, and an optimization module (Fig. 1). The event-list calendar is the core module of the designed simulation framework which controls the overall simulation processes via accurate adoption of the state variables. The notations used for defining the simulation procedures are listed in Table 1. The fleet is homogeneous with a maximum carrying capacity (C). Train capacity is considered as a hard constraint to analyze the congestion and overload effects explicitly in the simulation model. The passenger demand flows are aggregated and represented by time-dependent arrival rates and alighting fractions at stops. A total number of n train units are available. The double-track high-speed railroad infrastructure underlying in this study include 2 m stops for passenger boarding/alighting purposes (Fig. 2). The track capacity of all intermediate stations is limited. Thus, there is no permission for train overtaking. It is assumed that the passengers follow a first-come-firstserve (FCFS) policy when boarding a train. The dynamics of train traffic simulation involve the running time, dwell time, and headway control procedures. Eq. (1) defines the dynamics of train traffic in terms of the train dwell time at stops k and k + 1. The dynamic equation calculates the headway time of i + 1th train service as a function of previous headway time and the dwell time differences.

Hi (k + 1) = Hi (k ) + dwi (k )

1 (k ),

i

I {1},

k

K {2m} (1)

Eq. (2) implies the dynamic progress of passenger flows that evolve by increasing the number of passengers waiting to board the trains.

3. Simulation modeling framework

l i (k ) = l i

1 (k )

+ Hi (k )

i (k )

B i (k )

i (k ),

i

I {1},

k

K (2)

This section describes the designed event-oriented simulation model. The simulation platform is integrated with RSM optimization method to improve the response variable iteratively. The proposed computer simulation model is implemented in Enterprise Dynamics simulation tool. The designed simulation modeling framework is built upon the concepts and specifications of an object-oriented programming paradigm which facilitate the creation and simulation of the dynamic processes. The train operations are modeled using the creation and adjustment of different entities as the building blocks of the simulation model. Each stationary or moving entity, i.e., vehicles or track segments, has its attributes, a set of behaviors, and an event calendar to control the dynamic evolution of the system state variables. The simulation model is constructed in three parts and includes four basic

Likewise, the total number of passengers on ith train service, once it departs from stop k ∈ K, is calculated as follows:

ni (k ) = ni (k

1)

[Ai (k )

Bi (k )],

i

I,

k

K {1}

where the number of passengers who want to alight from i vice at stop k ∈ K is given by:

Ai (k ) =

i (k )

n i (k

1),

i

I,

k

K {1}

th

(3) train ser(4)

The evolution of passenger flows requires the calculation of the number of passengers who actually board ith train service at stop k ∈ K. Eq. (5) estimate the number of passengers boarded on the trains with respect to the remaining capacity of the arriving train.

Event-list control

Passenger flow control

dwi

Decision making module

Train flow control Scenario control module

Rail infrastructurescenario designand parameter setting

Passenger routing model Headway control module

Sensitivity analyzer for disturbances

Solution assessment module New parameters

Statistical test analyzer Data fitting module

Optimization module Surrogatefunction estimationandanalysis RSM optimizer

Adjusted headway time Fig. 1. The structure of the design event-driven simulation platform. 3

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Table 1 The definition of state variables. Symbol

Definition

Hi (k ) li (k ) ni (k ) Bi (k ) Ai (k ) i (k ) i (k ) Ri (k ) dwi (k ) i (k )

The The The The The The The The The The

headway time which separates the departure times of the i-1th and ith train services at stop k K number of passengers remains on stop k K just at the departure time of ith train service due to the capacity constraint estimated number of passengers on ith train service once it departs from stop k K number of passengers who actually board ith train service at stop k K number of passengers who want to alight from ith train service at stop k K estimated arrival rate of passengers at stop k K who want to board ith train service percentage of passengers alight from ith train service at stop k K travel time of ith train service from stop k 1 to k stop time of ith train service at stop k K number of passengers who failed to board ith train service at stop k due to the capacity restrictions

Inbound direction 1

2m

2

3

m-2

m-1

m

2m-1

2m-2

m+3

m+2

m+1

Outbound direction Fig. 2. The double-track high-speed railroad infrastructure underlying in this study.

Bi (k ) = min{li i

I {1},

1 (k )

+ Hi (k )

k

K {1}

i (k )

i (k ),

C

[1

i (k )]

ni (k

the train dwell times at stops based on the input and output flow of passengers. The evolution of the dynamic processes carries on until all remaining passengers reach their destinations. Finally, the average wait time per passenger is return as the response variable.

1)},

(5)

Given the above equations, the event-driven diagram of the dynamic evolution of passenger flows is shown in Fig. 3. The flow chart includes the system of equations associated with the passenger’s events, i.e., arrival process, boarding process, and alighting process, which evolve the passenger flow counts. For example, when new travelers arrive at the stops, they have to wait until a train arrives. Also, the flow diagram of the simulation processes associated with the train arrival/ departure events is shown in Fig. 4. The dynamic evolution of train operations involves the update of the state variables at the time when a new departure or arrival events occurs. As mentioned previously, the dynamic traffic control module verifies the departure of trains with respect to the minimum safety headway and the availability of a free track segment at a stop. In the cases when all these conditions are not satisfied, the involved train must wait at the platform. The eventhandler module controls the traffic simulation process and calculates

4. Response surface methodology Box and Wilson in 1951 developed Pseudo-response or pseudo-regression models. The use of RSM in the simulation was first discussed in 1975. RSM encounters the simulation model with the black-box approach and observes the input/output of the simulation model. RSM is an innovative method, which uses a sequence of local experiments and ultimately leads to an optimal combination of input variables. RSM also uses the DOE with linear regression analysis. Although real experiments may take a long time and it is impossible for a large sample size. RSM may be combined with white-box meta-models that provide an estimated gradient and improve its performance. In this method, it is assumed that input variables are continuous, and the output variable is

Fig. 3. Event-driven diagram of the dynamic evolution of passenger flows. 4

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Fig. 4. Flow diagram of the simulation processes associated with the train arrival/departure events.

mathematically expressed as an average response. The response variable can also be formulated as an expected value of a binary variable using the indicator function. In general, the application of the RSM method has been developed in simple optimization problems such as minimizing an objective function with continuous inputs without any constraint. Suppose E [w0 z ] as a response function that is minimized by the optimal combination of an input vector, where z = (z1, , z k )T ; j = 1, , k . Without loss of generality, vector z represents the initial input value, and vector x shows the coded input values. The objective is to minimize the expected value of the response function min z E [w0 z ]. The flowchart of the RSM method for the headway optimization problem is shown in Fig. 5. At the beginning of the algorithm, the main factors that affect the response variable are determined using factorial experimental design approaches. At this stage, the starting point and the dimensions of the first area are also determined. In the next step, the first-order regression quadratic model is fitted for the average response level as follows:

k

Y (x ) =

0

+

j zj

+

i=1

z nj

zj cin

zj = x i cin + z jn

j xj

+

(8)

The coefficients of the regression model { 0, 1, , n} are also estimated from variance analysis. The next step is to test the adequacy of the first-order model. Before applying the first-order model to improve the response, the adequacy test for the first-order model is performed to justify whether or not it can describe the response behavior. If the actual response shows the interactions between the factors or the pure curvature, the first-order model will be shown to be inadequate, which can be derived from the ANOVA. The goodness of fit test requires the factorial experiment of resolution III to estimate first-order regression coefficients that are not saturated, e.g., the total number of observations must be higher than the number of regression coefficients. If the firstorder model is accepted, the model is used to determine the direction in which the improvement in the response variable is expected. In this case, a linear search algorithm is applied to the starting point of the experimental space. The steepest descent search method is used for the maximization problem to moves towards the best points. The first-order regression model is used to find direction for the highest decrease based on the gradient information. It should be noted that the steepest gradient method is dependent on the scale, and therefore, the numbers are coded within the interval [−1, +1]. It should be noted that the step size of the movement is not known in the steepest descent search algorithm. Usually, an initial value is selected, and if the response value is higher than the current value, the step size changes to a smaller value, otherwise, the step size will be more significant. Finally, if the first-order model is not validated, there is a pure curvature or interaction between the factors in the current experimental space, in which case the value of P-value is lower

(6)

where the error term follows a standard normal probability distribution, i.e., NID (0, 2 ) . To increase the computational accuracy in estimating the parameters of the regression model, the input factors are encoded as follows:

xi =

+ i=1

k

Y (x ) =

0

(7)

The first-order model is expressed as follows 5

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Fig. 5. The flowchart of the RSM method for the headway optimization problem.

than the confidence level. Under these conditions, the first-order polynomial quadratic model is replaced with a second-order quadratic function as follows. k

y (x ) =

0

+

k j xj

j=1

k 2 j, j x j

+ j=1

to estimate the stationary point. If the estimation of the second-order model is as follows:

y =

k

+

j, i x j x i

+

NID (0,

2)

+ x 'b + x 'Bx

0

(10)

where

j=1 j< i

(9)

0

As indicated previously, regression coefficients of the second-order model are determined by regression analysis using the most central composite design (CCD) test. In the next step, similar to the first-order model, the adequacy of the second-order model in the current experimental space is tested using the variance analysis table. If the model adequacy is not confirmed, the size of the experimental space is reduced or the number of simulation replication increases. If the second-order model is sufficient, as mentioned in the previous section, a canonical analysis should be made to determine the stationary point of the second-order model. The optimal point for the second-order model is the superstructure in which the relative derivatives are zero. This point is called stationary. A stationary point may be a maximum point, a minimum point, or a saddle point. The following procedure can be used

1

b=

. . .

(11)

k

11

1k

2

B= sym .

kk

(12)

The stationary point is determined as follows:

y = b + 2Bx = 0 x

Fig. 6. Pareto chart for the standardized effects. 6

(13)

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

1 1 B b 2

xs =

estimated coefficients for the average waiting time per passenger in terms of coded values. The estimated value of the adjusted coefficient of determination, R2, shows 99.9% of the variance in the passenger wait time. The obtained value indicates that the coefficient of determination is adequately large, thereby there is a significant match between the simulated and estimated values. We observe that for all rows, the value of P-value is smaller than the α = 0.05 significance level, which indicate that all are significant. Since the curvature is significant, we need to look for a second-order regression model. Table 5 shows the result of ANOVA based on coded values. This table reports the statistical analysis of the main effects as well as two-way interactions between the inbound and outbound headway times at α = 0.05 significance level. The effects whose p-values are higher than α = 0.05 can be eliminated since they have no significant effect on the response. Table 6 reports the values of the estimated coefficients of the first-order regression model in terms of real values. SE Mean or standard error of the mean shows the variability between the average samples. Also, the standard deviation measures the variability within a sample. In fact, SE Mean shows the accuracy of the average estimate of a sample taken from a population, and the smaller it is, the better. The larger the sample, the smaller the SE Mean, and the more accurate average sample estimate. In this study, the number of replications is 10 per point. Since the single factors, i.e., “Inbound Headway” and “Outbound Headway” and their combined effects (Inbound Headway * Outbound Headway) are located behind the red reference line, they are all significant (Fig. 7). Normal probability plot of residuals is used to test the weather the points are normally distributed or not. The residuals versus order plot are also used to confirm the assumption of the independence of the residuals (observations) relative to each other. Observations will be independent of each other if no sequence diagram is shown. The existence of a pattern may be due to the correlation between the data and the consequent lack of independence between them. In general, the residuals should be randomly scattered around the center line. Residuals Versus Fitted Values must show a random pattern of residual values on both sides of the horizontal axis. In fact, this diagram is used to confirm the assumption that the variance of the random distribution of residuals is constant. The larger the fitted values increase, the distribution of residual values tends to increase. Any pattern in points may be a violation of the assumption that the variance is constant. The Cube plot is used to display the relationship between two to eight factors with or without response variables for 2 × 2 factorial design (Fig. 9). When the effect of a factor depends on the level of another factor, the interaction plot can be used to represent possible interactions (Fig. 10). Parallel lines indicate non-interconnection. The higher the difference between the slope of the lines, the higher the interaction effect. This plot is often used to visualize interactions in the variance analysis table. The main effects plot is used to check the difference between the average levels of one or more factors. Main effects are significant when the difference in the levels of a factor has an effect on the response variable. Fig. 11 shows the average response for each factor with a separate line. Based on the result obtained in all cases, the main effects, interactions, and curvature are significant. As a result, the CCD is carried out.

(14)

The value of the response variable in the stationary point is also obtained as follows:

ys =

0

+

1 ' xs b 2

(15)

After the optimum point is obtained, it should be determined whether this point is a maximal, minimum, or saddle point. For this purpose, the Eigenvector, i.e. i , is calculated.

|B

(16)

I| = 0

The type of stationary point is obtained through the sign of i . If all i values are negative, a maximal point obtained and if all i are positive, the resulting point is minimal. Otherwise, if i values have different signs, the result is a saddle point. Subsequently, the similarity of the optimum point with the output of the simulation is examined. If the equality test is rejected, we will return to the step of fitting the pseudosecond-order model with the use of CCD. Otherwise, the stationary point is accepted, and the algorithm is terminated. 5. Result and discussion This section presents the results of the factorial design experiments and the optimized headway obtained by the RSM. The simulation model is programmed using the Enterprise Dynamics simulation software platform developed by INCONTROL Simulation Solutions. All the computational experiments are implemented in MATLAB 2018b on Windows 10 system, and runs on laptop equipped with an Intel Core i7 CPU at 2.9 GHz and 16 gigabyte RAM. The designed SBO is implemented on real instances of Line No. 1 from Tehran underground metro system. This high-speed rail line is about 37 km long, which include 29 stations positioned from south to north. About 25% of stations are at surface level while the others are underground. 5.1. First-order model The process of fitting the full factorial design involves the calculation of the first-order regression meta-model as represented in Eq. (6). A full 2k factorial design considers the simulation experiments, and the pseudo-first-order model is developed. Parameter k is coded as ± 1 levels in which 0 characterizes the central points. The use of extra central points would enhance the resistance of the meta-model against curvature through estimation of the experimental error. The description of the factors used in the experiments is given in Table 2. A complete 22 factorial design is used with the two input factors, headway time for inbound and output routes, which are described in Table 2. The experimental test includes one replication at a corner and five replications at the central point for estimating the sum of squared errors (SSE) to check the validity of the model. The test points, along with the response values, are given in Table 3. To check the validity of the first-order model, we examine whether or not the fitted model reproduces the outputs of the simulation model. The validity of the meta-model is usually checked by applying a F statistics test. The simulation outputs for the 22 full factorial design are analyzed and the constants associated with the regression meta-model are calculated. The estimated effects and coefficients are obtained using MINITAB statistical package. Table 4 reports the results of effects and Table 2 The upper and lower values of the input factors. Input factors

Lower limit with the coded value of −1 (seconds)

Upper limit with the coded value of +1 (seconds)

Minimum coded value

Maximum coded value

Headway time for inbound route ( X1) Headway time for inbound route ( X2 )

240 240

1200 1200

−1 −1

+1 +1

7

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Table 3 The test points along with the response values. Simulation run

Central point

Block

X1

X2

Inbound headway time

Outbound headway time

Average wait time per passenger (seconds)

1 2 3 4 5 6 7 8 9

0 0 1 1 0 1 1 0 0

1 1 1 1 1 1 1 1 1

0 0 −1 +1 0 +1 −1 0 0

0 0 −1 −1 0 +1 +1 0 0

12 12 4 20 12 20 4 12 12

12 12 4 4 12 20 20 12 12

508.11 502.75 227.16 2231.67 521.13 1907.79 1507.50 509.96 508.03

Since the validity of the first-order model does not confirm, then the second-order model is built using a Box–Beckhen design method (see Figs. 6, 8, 12–14).

Table 4 Effects and estimated coefficients for the average waiting time per passenger in terms of coded values. Term

Effect

Coef

SE Coef

T

P

Constant X1 X2 X1 * X2 Ct Pt

1202.4 478.2 −802.1

1468.5 601.2 239.1 −401.1 −958.5

3.390 3.390 3.390 3.390 4.548

433.17 177.34 70.53 −118.30 −210.74

0.000 0.000 0.000 0.000 0.000

S = 2.99169

R-Sq = 100.00%

5.2. CDD and second-order quadratic model This section provides the result of experiments with CCD and second-order quadratic model. The experimental design is based on face centered and the value of Alpha = 1 in all experiments. The experimental points with their responses are given in Table 9. Table 10 reports the results of the fitting of the second-order regression quadratic model for the test points and the estimated regression coefficients in terms of coded values. If the p-value associated with the t-statistic is lower than the α, it can be concluded that the corresponding coefficient is significantly different from zero. If the value of p-value is greater than α, then the corresponding coefficient can be eliminated from the regression model. Based on the results, we see that for all rows, the value of P-value is smaller than the α = 0.05, so all of the effects are significant. S represents the standard deviation which shows the spread between the values of the data and the fitted values. “Lack-of-Fit” shows the data distance from the fitted model. The F statistic determines that whether or not a higher order model, which includes current predictor variables, is more appropriate. The F statistic is calculated as follows:

R-Sq(adj) = 99.99%

Table 5 The ANOVA result based on coded values. Source

DF

Seq SS

Adj SS

Adj MS

F

P

Main effects X1 X2 2-Way interactions X1 * X2 Curvature Residual error Pure error Total

2 1 1 1 1 1 4 4 8

1,674,470 1,445,766 228,704 643,380 643,380 2,041,750 184 184 4,359,784

1,674,470 1,445,766 228,704 643,380 643,380 2,041,750 184 184

837,235 1,445,766 228,704 643,380 643,380 2,041,750 46 46

18211.26 31447.82 4974.69 13994.60 13994.60 44411.47

0.000 0.000 0.000 0.000 0.000 0.000

Table 6 Estimated coefficients for the average waiting time in terms of real values. Term

Coef

Constant X1 X2 X1 * X2 Ct Pt

−694.316 150.348 105.087 −6.26648 −958.534

F=

AdjMS Residual ErrorAdjMS

(17)

Pure Error also specifies the variance in the response in the tested design points. To verify the relationship between the response and the predictive variables, the p-value for the Lack-of-Fit test should be compared with a threshold value, i.e., α. Estimated coefficients for the average waiting time in terms of real values are presented in Table 12. Given the pseudo-fitted model, the optimal solution is obtained as follows. y = 360.515

39.7354x1

Table 7 Least Squares Means for waiting time of each passenger at the station. X1

Mean

SE Mean

4 20

867.3 2069.7

4.794 4.794

X2 4 20

Mean 1229.4 1707.6

SE Mean 4.794 4.794

X1 * X2 4 4 20 4 4 20 20 20

Mean 227.2 2231.7 1507.5 1907.8

SE Mean 6.780 6.780 6.780 6.780

20.3684x2 + 5.86248x12 + 3.90633x 22

3.39954x1 x2

(18)

b=

xs =

39.7354 B= 20.3684

5.86248 1.69977

1.69977 3.90633

1 1 4.743287 B b= 4.671059 2

(19) (20)

To determine the nature of the optimum point, the Eigenvector of I mamatrix B is obtained by calculating the determinant of the B trix:

|B

I| = 0 1

2 9.76881 + 20.01156 = 0 = 6.845, 2 = 2.923

(21)

According to the obtained results, the point is a minimum point because of the Eigenvectors are positive. The value of the response variable is equal to: 8

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Table 8 The upper and lower values of the input factors for the CCD test. Input factor (headway)

Route

Lower limit with the coded value of −1 (seconds)

Upper limit with the coded value of +1 (seconds)

X1 X2

Inbound Outbound

240 240

600 600

Fig. 7. The normal probability plot for the standardized effects.

Fig. 8. The residual plot for average passenger wait time.

ys =

0

+

1 ' xs b = 218.7058 2

H0: y¯ = 218.7058 H : y¯ 218.7058

(22)

The output of the simulation model for the optimal solution is 219.983. In order to investigate the validity of the RSM, the simulation model is repeated for n = 10 replications. The hypothesis test for the difference between the average response obtained from the simulation model with the value obtained from the meta-model at 97.5% confidence level is expressed as follows:

(23)

The results obtained from the simulation model are as follows: The value of the t-test statistic is equal to:

t =

X¯ S

µ n

=

219.983 218.7058 2.607 10

t (0.975,9) = 2.262 9

= 1.55

t < t(0.975,9)

(24)

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Fig. 9. The cube plot for average wait time.

Fig. 10. The interaction plot for average wait time.

Fig. 11. The main effects plot for average wait time.

As a result, the H0 hypothesis is not rejected, and the optimal point is accepted and thus validated. Fig. 15 shows the boxplot for 97.5% tconfidence interval for average values.

Fig. 16 illustrates the comparison of the output obtained from the fitted pseudo-model and simulation model vs. a number of replications. It is expected that with the increase in the number of simulations, the 10

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Fig. 12. The residual plot.

difference between the output of the pseudo-fit model and the output of the simulation model is reduced. Based on the results (Fig. 17), it is seen that as the number of replications increase, the variance of the output, i.e., the average waiting time per passenger, is decreased and the outputs go towards convergence. The oscillation in the simulation output is also due to the random nature of the problem. Similarly, Fig. 18 shows that with increasing the number of replications, the optimal solution obtained from the RSM goes towards convergence. Based on the results, the difference between the optimal solution obtained from the RSM and the simulation model is reduced with increasing number of simulations, and the solutions go towards convergence. Fig. 19 shows the comparison of the output obtained from the quasi-fit model and the simulation model by adding points to the test space. According to the results, it is seen that the output from the quasifitted model for the average passenger waiting time goes toward convergence and the difference between the responses is decreasing by adding points to the solution space (see Tables 7 and 8, 11, 13). Figs. 20 and 21 also show the optimal headway obtained by quasi-

Fig. 13. The contour plot for the fitted RSM model.

Fig. 14. The surface plot for the fitted RSM model. 11

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Table 9 The experimental points with their responses. Design of experimental points

Type of point

X1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

1 0 0 −1 1 0 1 1 1 1 0 1 −1 −1 0 0 0 0 0 −1 −1 0 −1 −1 1 −1

−1 0 0 −1 −1 0 1 1 −1 −1 0 1 −1 0 0 0 0 0 0 0 1 0 0 0 1 1

X2

−1 0 0 0 −1 0 1 −1 1 1 0 1 0 1 0 0 0 0 0 1 0 0 −1 −1 −1 0

Headway time (inbound route)

Headway time (outbound route)

Average wait time per passenger (seconds)

4 7 7 4 4 7 10 10 4 4 7 10 4 7 7 7 7 7 7 7 10 7 7 7 10 10

4 7 7 7 4 7 10 4 10 10 7 10 7 10 7 7 7 7 7 10 7 7 4 4 4 7

231.30 243.61 242.11 251.77 230.08 252.19 375.93 388.85 351.04 328.37 245.80 375.56 243.89 327.04 243.97 249.90 248.09 249.16 245.08 326.15 381.98 240.82 251.98 247.13 390.53 389.52

Table 12 Estimated coefficients for the average waiting time in terms of real values. Term

Coef

Constant X1 X2 X1 * X1 X2 * X2 X1 * X2

368.865 −48.9479 −13.2162 6.52533 3.33477 −3.41556

Fig. 15. The boxplot for 97.5% t-confidence interval for average values.

Table 10 Estimated regression coefficients for the average waiting time per passenger in terms of coded values. Term

Coef

SE Coef

T

P

Constant X1 X2 X1 * X1 X2 * X2 X1 * X2

249.50 55.49 28.68 58.73 30.01 −30.74

3.578 3.518 3.518 5.185 5.185 4.309

69.726 15.774 8.154 11.326 5.788 −7.134

0.000 0.000 0.000 0.000 0.000 0.000

S = 12.1871 PRESS = 6037.92 R-Sq = 96.84% R-Sq(pred) = 93.59% R-Sq (adj) = 96.06%.

Fig. 16. Comparison of the output obtained from the fitted meta-model and simulation model vs. the number of replications.

Table 11 The ANOVA result based on coded values for the second-order model. Source

DF

Seq SS

Adj SS

Adj MS

F

P

Regression Linear X1 X2 Square X1*X1 X2*X2 Interaction X1*X2 Residual Error Lack-of-Fit Pure Error Total

5 2 1 1 2 1 1 1 1 20 3 17 25

91171.1 46828.1 36954.1 9874.0 36783.5 31807.7 4975.7 7559.6 7559.6 2970.5 2519.2 451.3 94141.6

91171.1 46828.1 36954.1 9874.0 36783.5 19051.5 4975.7 7559.6 7559.6 2970.5 2519.2 451.3

18234.2 23414.0 36954.1 9874.0 18391.7 19051.5 4975.7 7559.6 7559.6 148.5 839.7 26.5

122.77 157.64 248.81 66.48 123.83 128.27 33.50 50.90 50.90

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

31.63

0.000

Fig. 17. Optimal inbound headway obtained by RSM vs. the number of replications.

model fitted by adding points to the experimental space for inbound and outbound routes, respectively. It can be observed that the output from the quasi-fitted model converges by adding points to the experimental space for both inbound and outbound routes, and the difference between the responses is decreasing. 12

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

Fig. 18. Optimal outbound headway obtained by RSM vs. the number of replications.

Fig. 20. Optimal inbound headway obtained by quasi-model fitted by adding points to the experimental space.

Fig. 19. Comparing the output obtained from the quasi-fit model and the simulation model by adding points to the test space. Table 13 Results of simulation outputs for mean difference test. Replication No.

The output of the simulation

1 2 3 4 5 6 7 8 9 10

219.25 218.39 224.81 222.31 222.25 218.33 220.85 216.33 217.26 220.05

Fig. 21. Optimal outbound headway obtained by quasi-model fitted by adding points to the experimental space.

computational results indicate that the developed RSM could provide optimal values of control factors to ensure a desirable level of service for passengers in terms of average waiting time per passenger. As future work, the proposed meta-model simulation–optimization approach can be extended to generate energy-efficient train schedules. It is also valuable to estimate other metrics, i.e. robustness and reliability of the train schedule, through hybrid SBO. Also, the developed SBO method would benefit from a reduction of the computational complexity of the model. Consequently, future researches need to be directed on using neural networks or deep learning algorithms to speed up the search process, thus improving the efficiency of the proposed SBO.

6. Conclusion The daily operations in high-speed rail transit systems are sensitive to several uncertain factors, i.e., stochastic passenger flows and running time variations that may affect the baseline schedule. This study addressed the optimal design of train timetable in a high-speed railway using hybrid meta-modeling and event-driven simulation method. The modeling framework uses the response surface methodology as a regression-based simulation-based optimization approach to minimize the passenger wait time on platforms. Several complete factorial designs of DOE are built and tested to determine the significance of the main effects on the response variable. The validity of the designed simulation-optimization model was confirmed using real instances of Tehran metro networks. The results of ANOVA and residual plots specify the implication of the second-order surrogate model, thereby suggesting it tolerably characterizes an accurate relationship between the train headway times and the response variable. To sum up, the

References Antal, M., Pop, C., Petrican, T., Vesa, A. V., Cioara, T., Anghel, I., ... NiewiadomskaSzynkiewicz, E. (2019). MoSiCS: modeling, simulation and optimization of complex systems – A case study on energy efficient datacenters. Simulation Modelling Practice and Theory, 93, 21–41. Azad, N., Hassini, E., & Verma, M. (2016). Disruption risk management in railroad networks: An optimization-based methodology and a case study. Transportation Research Part B: Methodological, 85, 70–88. Bookbinder, J. H., & Desilets, A. (1992). Transfer optimization in a transit network. Transportation Science, 26, 106–118. Canca, D., Barrena, E., Zarzo, A., Ortega, F., & Algaba, E. (2012). Optimal train reallocation strategies under service disruptions. Procedia-Social and Behavioral

13

Computers & Industrial Engineering 138 (2019) 106110

E. Hassannayebi, et al.

approach. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 1–16. Mannino, C., & Mascis, A. (2009). Optimal real-time traffic control in metro stations. Operations Research, 57, 1026–1039. Mascis, A., & Pacciarelli, D. (2002). Job-shop scheduling with blocking and no-wait constraints. European Journal of Operational Research, 143, 498–517. Mo, P., Yang, L., Wang, Y., & Qi, J. (2019). A flexible metro train scheduling approach to minimize energy cost and passenger waiting time. Computers & Industrial Engineering, 132, 412–432. Ning, B., Xun, J., Gao, S., & Zhang, L. (2014). An integrated control model for headway regulation and energy saving in urban rail transit. IEEE Transactions on Intelligent Transportation Systems, 16, 1469–1478. Parbo, J., Nielsen, O. A., & Prato, C. G. (2016). Passenger perspectives in railway timetabling: A literature review. Transport Reviews, 36, 500–526. Pouryousef, H., Lautala, P., & Watkins, D. (2016). Development of hybrid optimization of train schedules model for N-track rail corridors. Transportation Research Part C: Emerging Technologies, 67, 169–192. Rahimi Mazrae Shahi, M., Fallah Mehdipour, E., & Amiri, M. (2016). Optimization using simulation and response surface methodology with an application on subway train scheduling. International Transactions in Operational Research, 23, 797–811. Sajedinejad, A., Mardani, S., Hasannayebi, E., Mir Mohammadi K, S. & Kabirian, A. (2011). Simarail: simulation based optimization software for scheduling railway network. Simulation Conference (WSC), Proceedings of the 2011 Winter, 2011. IEEE, 3730–3741. Shakibayifar, M., Hassannayebi, E., Jafary, H., & Sajedinejad, A. (2017a). Stochastic optimization of an urban rail timetable under time-dependent and uncertain demand. Applied Stochastic Models in Business and Industry, 33, 640–661. Shakibayifar, M., Hassannayebi, E., Mirzahossein, H., Taghikhah, F., & Jafarpur, A. (2019). An intelligent simulation platform for train traffic control under disturbance. International Journal of Modelling and Simulation, 39, 135–156. Shakibayifar, M., Sheikholeslami, A., Corman, F., & Hassannayebi, E. (2017b). An integrated rescheduling model for minimizing train delays in the case of line blockage. Operational Research, 1–29. Shi, J., Yang, L., Yang, J., & Gao, Z. (2018). Service-oriented train timetabling with collaborative passenger flow control on an oversaturated metro line: An integer linear optimization approach. Transportation Research Part B: Methodological, 110, 26–59. Vázquez-Abad, F. J., & Zubieta, L. (2005). Ghost simulation model for the optimization of an urban subway system. Discrete Event Dynamic Systems, 15, 207–235. Wong, R. C., Yuen, T. W., Fung, K. W., & Leung, J. M. (2008). Optimizing timetable synchronization for rail mass transit. Transportation Science, 42, 57–69. Wu, J., Liu, M., Sun, H., Li, T., Gao, Z., & Wang, D. Z. W. (2015). Equity-based timetable synchronization optimization in urban subway network. Transportation Research Part C: Emerging Technologies, 51, 1–18. Xu, X., Li, H., Liu, J., Ran, B., & Qin, L. (2019). Passenger flow control with multi-station coordination in subway networks: Algorithm development and real-world case study. Transportmetrica B: Transport Dynamics, 7, 446–472. Yalcinkaya, Ö., & Bayhan, G. M. (2009). Modelling and optimization of average travel time for a metro line by simulation and response surface methodology. European Journal of Operational Research, 196, 225–233. Yin, J., Yang, L., Tang, T., Gao, Z., & Ran, B. (2017). Dynamic passenger demand oriented metro train scheduling with energy-efficiency and waiting time minimization: Mixedinteger linear programming approaches. Transportation Research Part B: Methodological, 97, 182–213.

Sciences, 54, 402–413. Dengiz, B., & Belgin, O. (2014). Simulation optimization of a multi-stage multi-product paint shop line with Response Surface Methodology. Simulation, 90, 265–274. Dong, H., Ning, B., Cai, B., & Hou, Z. (2010). Automatic train control system development and simulation for high-speed railways. IEEE Circuits and Systems Magazine, 10, 6–18. Eberlein, X. J., Wilson, N. H., Barnhart, C., & Bernstein, D. (1998). The real-time deadheading problem in transit operations control. Transportation Research Part B: Methodological, 32, 77–100. Eberlein, X. J., Wilson, N. H. & Bernstein, D. (1999). Modeling real-time control strategies in public transit operations. Computer-Aided Transit Scheduling. Springer. Flamini, M., & Pacciarelli, D. (2008). Real time management of a metro rail terminus. European Journal of Operational Research, 189, 746–761. Hassannayebi, E., Sajedinejad, A., & Mardani, S. (2014). Urban rail transit planning using a two-stage simulation-based optimization approach. Simulation Modelling Practice and Theory, 49, 151–166. Hassannayebi, E., Sajedinejad, A. & Mardani, S. 2016a. Disruption management in urban rail transit system: A simulation based optimization approach. Handbook of Research on Emerging Innovations in Rail Transportation Engineering. IgI-Global. Hassannayebi, E., Sajedinejad, A., & Mardani, S. (2016b). Disruption management in urban rail transit system: A simulation based optimization approach. Handbook of Research on Emerging Innovations in Rail Transportation Engineering, 420–450. Hassannayebi, E., & Zegordi, S. H. (2017). Variable and adaptive neighbourhood search algorithms for rail rapid transit timetabling problem. Computers & Operations Research, 78, 439–453. Hassannayebi, E., Zegordi, S. H., Amin-Naseri, M. R., & Yaghini, M. (2018). Optimizing headways for urban rail transit services using adaptive particle swarm algorithms. Public Transport, 10, 23–62. Högdahl, J., Bohlin, M., & Fröidh, O. (2019). A combined simulation-optimization approach for minimizing travel time and delays in railway timetables. Transportation Research Part B: Methodological, 126, 192–212. Kang, L., Zhu, X., Sun, H., Puchinger, J., Ruthmair, M., & Hu, B. (2016). Modeling the first train timetabling problem with minimal missed trains and synchronization time differences in subway networks. Transportation Research Part B: Methodological, 93, 17–36. Kim, M., Yoon, H. J., Son, S. H. & Eun, Y. (2019) Data-based model of metro scheduling for passenger wait-time optimization with constraints: WIP abstract. 2019. ACM, 326–327. Korytkowski, P., Wiśniewski, T., & Rymaszewski, S. (2013). An evolutionary simulationbased optimization approach for dispatching scheduling. Simulation Modelling Practice and Theory, 35, 69–85. Lai, D. S. W., & Leung, J. M. Y. (2018). Real-time rescheduling and disruption management for public transit. Transportmetrica B: Transport Dynamics, 6, 17–33. Li, D., Zhang, T., Dong, X., Yin, Y., & Cao, J. (2019a). Trade-off between efficiency and fairness in timetabling on a single urban rail transit line under time-dependent demand condition. Transportmetrica B: Transport Dynamics, 1–29. Li, S., Xu, R., & Han, K. (2019b). Demand-oriented train services optimization for a congested urban rail line: Integrating short turning and heterogeneous headways. Transportmetrica A: Transport Science, 15, 1459–1486. Liu, P., Yang, L., Gao, Z., Huang, Y., Li, S., & Gao, Y. (2018a). Energy-efficient train timetable optimization in the subway system with energy storage devices. IEEE Transactions on Intelligent Transportation Systems, 19, 3947–3963. Liu, R., Li, S., Yang, L., & Yin, J. (2018b). Energy-efficient subway train scheduling design with time-dependent demand based on an approximate dynamic programming

14