Pareto optimal train scheduling for urban rail transit using generalized particle swarm optimization

Abstract

In urban rail transport, train timetable plays a crucial role, whose quality determines the whole system’s performance to a large extent. In practical urban rail operation, two contradictive aspects—service quality and operation cost—should be considered during train scheduling. A good train timetable should achieve considerable service quality with as little operation cost as possible. Previously, many studies have been conducted specific to urban rail train scheduling, although most of them do not put enough emphasis on its multi-objective nature. In this article, therefore, Pareto optimal urban rail train scheduling which can give more instruction to practical operation is studied. First, referring to some existing studies, the problem is reasonably defined, which takes time-dependent origin–destination demand as the input and aims at minimizing the passengers’ total travel time and the number of used train stocks. Then, an efficient iteration algorithm and a valid train stock assignment procedure are designed to calculate the passengers’ total travel time and required train stock number, respectively. On that basis, the studied problem is reasonably formulated as a bi-objective optimization model and a Pareto-based particle swarm optimization procedure is designed to solve it. Finally, with two different scaled urban rail lines, the whole methodology is illustrated and the algorithm is tested.

Introduction

Public transport has been recognized as one of the most effective methods to reduce urban traffic congestion. Among the common public transport modes, rail transport is more favored due to its high speed and capacity. Nowadays, in almost all metropolises around the world, rail transit is applied and plays an important role in the daily trips of residents.

Compared with other transport modes, the operation of rail transport is more dependent on advance planning, which generally consists of five phases: ¹ demand analysis, line planning, train timetabling, rolling stock planning, and crew scheduling. In this planning system, train timetable plays a crucial role, which provides the foundation for the following phases, that is, rolling stock and crew scheduling. The quality of train timetable, therefore, determines the performance of the whole system to a large extent. In early times, when train timetabling can only be conducted manually, the quality of timetable was mainly determined by designer’s experience. With the development of computer technology, automatic train timetabling draws more and more attention and many studies have been conducted on it. Some studies are focused on the description of designer’s experience with computer language,^2,3 while more other studies are conducted from the perspective of mathematical optimization.^4–11

As a particular form of rail transport, the operation of urban rail transit also requires well-predefined train timetable. Compared with intercity railways, however, urban rail systems generally have the following features: ¹²

1. Detailed origin–destination (OD) demand information can be easily collected through the automatic fare collection (AFC) system.

2. Urban rail line is always designed as double-track and each direction has a separate rail track. Hence, trains in different directions do not interact and can be scheduled separately.

3. Skip-stop pattern is not allowed and trains make stop at every station they pass.

4. Trains get similar performance (speed), so trains do not overtake each other during operation.

With these features, the train timetabling problem of urban rail transport is remarkably different with that of intercity rail system. On one hand, due to the absence of train crossing and overtaking, the decision variables in timetabling are much reduced. At the meantime, with the systematic travel demand information, the timetabling can be conducted in a more refined and demand-oriented manner.

Previously, quite a number of studies have been conducted specific to urban rail train timetabling problem and most of them are from mathematical optimization perspective. Cury et al. ¹³ proposed a nonlinear optimization model for the train timetabling problem of São Paulo Underground Railway System. In that model, a specific number of trains are given and the headways between consecutive trains need to be determined to minimize the total cost which combines the passengers’ total waiting time, the number of trains, the deviation between the actual number and the desired number of the passengers in the train, and the deviation between the actual and reference control sequence. Also on the train timetabling problem of São Paulo Underground Railway System, Wânderson and Basílio ¹⁴ proposed a new linear programming approach, which can be solved more efficiently. This approach, however, is only suited for periods of time where the passenger demand can be considered constant. Kwan and Chang ¹⁵ studied the scheduling of rapid transit systems from Pareto optimization perspective, in which operating costs and passenger dissatisfaction level are considered as optimization objectives and average train frequency, average dwell time, average coast level, and average layover time are taken as variables, and proposed a heuristic-based evolutionary algorithm to achieve that. Niu and Zhou ¹⁶ studied the urban rail timetabling under oversaturated conditions, taking into consideration that the passengers may not be able to board the next arrival train due to the limited train capacity. With a reasonable description of passenger loading and departure events, they proposed a nonlinear train scheduling model, whose objective is to minimize the total number of waiting passengers and weighted remaining passengers. Zhu et al. ¹⁷ also studied the urban rail timetabling with capacity constraints. In their study, the capacity constraint of station is also considered and a simple simulation model is proposed to calculate the passenger waiting cost. On this basis, a single-objective optimization model is proposed and generic algorithm is applied to solve it. Wang et al. ¹² studied the train scheduling problem for an urban rail transit network. With an event-driven model describing network status, a nonlinear non-convex optimization model is designed, in which the objective function is a weighted sum of the passengers’ total travel time (PTTT) and trains’ energy consumption.

As can be seen, the multi-objective nature of urban rail timetabling has been widely recognized and the objectives used in the existing studies can be divided into two main aspects—service quality and operation cost. However, the existing studies do not cope with the multiple objectives very well. In most of the above studies, simple weighted aggregation is used to integrate the multiple objectives.^{12–14,16,17} However, the practicability of these weighted aggregation methods cannot be guaranteed, since there is still no definite unified opinion on how to determine the weight proportion of each aspect. In contrast, Pareto optimality method is more practical, which can provide more information to operators when determining an optimal train timetable, that is, list more compromise solutions to let the operator choose according to the circumstances of urban rail system. Kwan and Chang ¹⁵ introduced Pareto optimality, but the variables in their model are relatively simplified, which limits its application. In this article, therefore, a Pareto optimal urban rail train scheduling problem, which is more realistic and helpful, is proposed, in which the PTTT and the number of used train stocks are chosen as objectives to represent the service quality and operation cost, respectively.

The rest of this article is organized as follows. In the next section, the studied problem is elaborated and the notations used in this article are defined. Then, in section “Performance assessment of urban rail train timetable,” an efficient iteration algorithm and a valid train stock assignment procedure are designed to calculate the PTTT and the required train stock number, respectively, based on which a bi-objective mathematical optimization model is proposed in section “Model formulation.” In section “Pareto optimal solution using generalized PSO,” a Pareto-based particle swarm optimization (PSO) procedure is designed to solve the proposed model. After that, in section “Numerical example,” two example urban rail lines are introduced to illustrate the proposed methodology and test the algorithm. Finally, conclusions and discussions are given in section “Conclusion and discussions.”

Problem statement

With reference to the existing studies,^13,14,16,17 a bi-directional urban rail line with only one rail yard located at one end is considered in our studies, whose layout is shown in Figure 1. The station on the end where the yard is located is called the start terminal, while the station on the other end is called return terminal. The direction from the start terminal to the return terminal is defined as down-bound and the opposite direction is up-bound. The stations are numbered as 1, 2,…, m sequentially along the down-bound direction. Among all the stations, only the two terminal stations are equipped with crossovers and can accommodate train turning-back.

Figure 1.

Layout of the studied urban rail line.

During operation, the trains depart from the yard and enter the line from the start terminal. After a down-bound trip with making stop at each station, trains reach the return terminal, in which they are emptied and then conduct turning-back to begin an all-stop up-bound trip. When the trains return to the start terminal, they are faced with two choices—going back into the rail yard or turning back to start a new round. In this process, as analyzed in section “Introduction,” the trains have the similar performance so that they do not overtake each other and trains are assumed to follow the same running time in each section and the same dwell time at each station.^16,17 Besides, the headway between the arrivals of two consecutive trains at the start terminal from the rail yard should ensure that no conflict happens between them in the following operation procedures. More specifically, the headway between the arrivals of two consecutive trains at the start terminal must be not less than the minimum safety headway between two consecutive running trains in the section, two consecutive arrivals at each station, or two turning-backs in the return terminal.

As for the passengers, their arrivals at each station follow a non-homogeneous Poisson process with arrival rate varying with time. In order to describe that, the scheduled time horizon is typically divided into several time periods and the arrival rate or the number of arrivals at each station in each time period is given. On the platform, the passengers board the train on a first-come-first-served (FCFS) basis, which has been proven a valid assumption in urban transit-related study.

On the above conditions, the train timetable can be considered to be composed of a series of train round trips (abbreviated as TRTs), each of which is made up of a down-bound trip and an up-bound trip. Once the start time of each TRT is determined, a definite train timetable can be derived. Therefore, the decision variables of the train timetabling problem here are the start times of these TRTs. As for the objective, PTTT and the required number of train stocks, which represent transit service quality and operation cost, respectively, are chosen as the main performance criteria of train timetable, as analyzed in section “Introduction.” To summarize, the train timetabling problem studied here can be described as determining the start time of each TRT to achieve a minimum PTTT with as few train stocks as possible.

Based on the above statements, the parameters and variables used in this study are identified. The corresponding notations are listed in Appendix 1.

Performance assessment of urban rail train timetable

Calculating PTTT

As elaborated in section “Introduction,” the key point in calculating PTTT is the reasonable consideration of train capacity constraint. Here, in order to achieve efficient calculation of PTTT under capacity constraints, an iteration algorithm is proposed.

As analyzed in section “Problem statement,” once the decision variables $x_1,x_2,\dots ,x_i,\dots ,x_n$ are determined, a definite train timetable can be derived. The departure times from station in each round trip, that is, $t_{i,u}^{\mathrm{down}}(\forall i,u)$ and $t_{i,u}^{\mathrm{up}}(\forall i,u)$ , can be calculated, as shown in equations (1) and (2)

$$t_{i,u}^{\mathrm{down}}=x_i+{\displaystyle \sum _{v=1}^{u-1}(e_{v,v+1}^{\mathrm{down}}+d_v^{\mathrm{down}})},\forall i,u$$

(1)

$$\begin{array}{l}{t}_{i,u}^{\mathrm{up}}=x_i+{\displaystyle \sum _{v=1}^p(e_{v,v+1}^{\mathrm{down}}+d_v^{\mathrm{down}})}+{\delta }_{\mathrm{return}}\\ +{\displaystyle \sum _{v=u+1}^p(e_{v-1,v}^{\mathrm{up}}+d_v^{\mathrm{up}})},\forall i,u\end{array}$$

(2)

Since the trains need to be emptied before turning back at the return terminal, the up-bound and down-bound passenger transit processes are relatively separated. Therefore, given a definite train timetable, the PTTT of each direction can be calculated separately. In the following, the calculation of down-bound PTTT (i.e. $W_{\mathrm{down}}$ ) is elaborated:

• Step 0: initialization. Set $i=1$ , $u=1$ , $W_{\mathrm{down}}=0$ , $r_{u,v}^k=0(\forall u,v,k)$ , $c_i=C(\forall i)$ , and $o_i^u=0(\forall i,u)$ . Number the TRTs in the timetable sequentially according to their start time.

• Step 1: calculate the composition of down-bound passengers who are trying to get on the train of round trip i from the platform of station $u$ , that is, updating $r_{u,v}^k(\forall k,v=u+1,\dots ,m)$ at time $t_{i,u}^{\mathrm{down}}$ . Since the passengers’ arrivals at station follow Poisson distribution, the expected number of the passengers waiting on the platform arriving in each period can be calculated according to a time proportion method.

• 1. If the TRT i is the first one in the timetable, that is, $i=1$ , the calculation of $r_{u,v}^k$ need only consider the departure time of this train $t_{i,u}^{\mathrm{down}}$ . If the train of round trip i departs from station u after time period k, then all the passengers arriving in period k are waiting on the platform for this train, that is, $r_{u,v}^k=q_{u,v}^k$ . If the train of round trip i departs from station u in time period k, then the expected number of the passengers on the platform arriving in this period can be calculated according to time proportion $(t_{i,u}^{\mathrm{down}}-{\epsilon }_{k-1})/({\epsilon }_k-{\epsilon }_{k-1})$ . If the train of round trip i departs from the station before time period k, then none of the passengers on the platform is arriving in period k. In conclusion, the $r_{u,v}^k$ should be updated as

$$r_{u,v}^k=\{\begin{array}{cc}{q}_{u,v}^k\hfill & {\epsilon }_k\le t_{i,u}^{\mathrm{down}}\hfill \\ \multicolumn{1}{c}{\frac{t_{i,u}^{\mathrm{down}}-{\epsilon }_{k-1}}{{\epsilon }_k-{\epsilon }_{k-1}}{q}_{u,v}^k}& {\epsilon }_{k-1}<t_{i,u}^{\mathrm{down}}<{\epsilon }_k\hfill \\ \multicolumn{1}{c}{0}& t_{i,u}^{\mathrm{down}}\le {\epsilon }_{k-1}\hfill \end{array},\forall k,v=u+1,\dots ,m$$

(3)

• 2. If the TRT i is not the first one, that is, $i>1$ , then $r_{u,v}^k$ should be updated on the basis of the status when the previous down-bound train departs from the same station, that is, the departure time of previous train $t_{i-1,u}^{\mathrm{down}}$ should also be considered. If the previous train and this train both depart after time period k, then $r_{u,v}^k$ stays the same. If the previous train and this train both depart in time period k, then $r_{u,v}^k$ should be calculated according to time proportion $(t_{i,u}^{\mathrm{down}}-t_{i-1,u}^{\mathrm{down}})/({\epsilon }_k-{\epsilon }_{k-1})$ . If only the previous train departs in time period k, $r_{u,v}^k$ should be calculated according to time proportion $({\epsilon }_k-t_{i-1,u}^{\mathrm{down}})/({\epsilon }_k-{\epsilon }_{k-1})$ . If only this train departs in time period k, $r_{u,v}^k$ should be calculated according to time proportion $(t_{i,u}^{\mathrm{down}}-{\epsilon }_{k-1})/({\epsilon }_k-{\epsilon }_{k-1})$ . If the previous train departs before time period k and this train departs after time period k, then all the passengers arriving in time period k should be included. If this train arrives before the time period k, then $r_{u,v}^k$ remains at 0. In conclusion, $r_{u,v}^k$ should be updated as

$$r_{u,v}^k=\{\begin{array}{cc}{r}_{u,v}^k\hfill & {\epsilon }_k\le t_{i-1,u}^{\mathrm{down}}\hfill \\ \multicolumn{1}{c}{r_{u,v}^k+\frac{t_{i,u}^{\mathrm{down}}-t_{i-1,u}^{\mathrm{down}}}{{\epsilon }_k-{\epsilon }_{k-1}}{q}_{u,v}^k}& {\epsilon }_{k-1}\le t_{i-1,u}^{\mathrm{down}}<t_{i,u}^{\mathrm{down}}\le {\epsilon }_k\hfill \\ \multicolumn{1}{c}{r_{u,v}^k+\frac{{\epsilon }_k-t_{i-1,u}^{\mathrm{down}}}{{\epsilon }_k-{\epsilon }_{k-1}}{q}_{u,v}^k}& {\epsilon }_{k-1}\le t_{i-1,u}^{\mathrm{down}}\le {\epsilon }_k<t_{i,u}^{\mathrm{down}}\hfill \\ \multicolumn{1}{c}{r_{u,v}^k+\frac{t_{i,u}^{\mathrm{down}}-{\epsilon }_{k-1}}{{\epsilon }_k-{\epsilon }_{k-1}}{q}_{u,v}^k}& t_{i-1,u}^{\mathrm{down}}<{\epsilon }_{k-1}\le t_{i,u}^{\mathrm{down}}\le {\epsilon }_k\hfill \\ \multicolumn{1}{c}{r_{u,v}^k+q_{u,v}^k}& t_{i-1,u}^{\mathrm{down}}<{\epsilon }_{k-1}\le {\epsilon }_k<t_{i,u}^{\mathrm{down}}\hfill \\ \multicolumn{1}{c}{0}& t_{i,u}^{\mathrm{down}}\le {\epsilon }_{k-1}\hfill \end{array},\forall k,v=u+1,\dots ,m$$

(4)

• Step 2: when the train of round trip i arrives at station $u$ , there are still passengers destined for the following stations staying on the train, so the available capacity of the train is

$$c_i=C-\sum _{v=u+1}^m{o}_i^v$$

(5)

Update $c_i$ accordingly and set $o_i^u=0$ .

• Step 3: update $o_i^v(\forall v=u+1,\dots ,m)$ , $r_{u,v}^k(\forall v=u+1,\dots ,m)$ , and $W_{\mathrm{down}}$ at time of the train of round trip i departing from station u, that is, $t_{i,u}^{\mathrm{down}}$ . The procedure is as follows:

•  Step 3.1: set $k=1$

•  Step 3.2: if $r_{u,v}^k=0(\forall v=u+1,\dots ,m)$ , then go to Step 3.4. Otherwise, go to Step 3.3.

•  Step 3.3: if $\sum _{v=u+1}^m{r}_{u,v}^k\le c_i$ , all the passengers waiting on the platform can board the train, so update $o_i^v$ , $r_{u,v}^k$ , and $W_{\mathrm{down}}$ according to equations (6)–(8), respectively, and go to Step 3.4

$$o_i^v=o_i^v+r_{u,v}^k,\forall v=u+1,\dots ,m$$

(6)

$$W_{\mathrm{down}}=W_{\mathrm{down}}+\sum _{v=u+1}^m{r}_{u,v}^k(t_{i,u}^{\mathrm{down}}-{\epsilon }_{k-1})$$

(7)

$$r_{u,v}^k=0,\forall v=u+1,\dots ,m$$

(8)

If $\sum _{v=u+1}^m{r}_{u,v}^k>c_i$ , only a part of the passengers waiting on the platform can board the train, so update $o_i^v$ , $r_{u,v}^k$ , and $W_{\mathrm{down}}$ according to equations (9)–(11), respectively, and go to Step 4

$$o_i^v=o_i^v+\lfloor \frac{c_i{r}_{u,v}^k}{\sum _{v=u+1}^m{r}_{u,v}^k}\rfloor ,\forall v=u+1,\dots ,m$$

(9)

$$W_{\mathrm{down}}=W_{\mathrm{down}}+\sum _{v=u+1}^m\lfloor \frac{c_i{r}_{u,v}^k}{\sum _{v=u+1}^m{r}_{u,v}^k}\rfloor (t_{i,u}^{\mathrm{down}}-{\epsilon }_{k-1})$$

(10)

$$r_{u,v}^k=r_{u,v}^k-\lfloor \frac{c_i{r}_{u,v}^k}{\sum _{v=u+1}^m{r}_{u,v}^k}\rfloor ,\forall v=u+1,\dots ,m$$

(11)

where $\lfloor \rfloor$ is the sign of floor function.

•  Step 3.4: if ${\epsilon }_k<t_{i,u}^{\mathrm{down}}$ , then set $k=k+1$ and go back to Step 3.3. Otherwise, go to Step 4.

• Step 4: if $u<m$ , then set $u=u+1$ and go back to Step 1. Otherwise, set $u=0$ and go to Step 5.

• Step 5: if $i<n$ , then set $i=i+1$ and go back to Step 1. Otherwise, stop and output $W_{\mathrm{down}}$ and $\mathbf{R}=[r_{u,v}^k{]}_{m\times m\times K}$ .

The above procedure is further illustrated in Figure 2. Through this procedure, the expected total travel time of down-bound passengers can be efficiently calculated. In a similar way, the up-bound PTTT $W_{\mathrm{up}}$ can also be calculated and then the PTTT of the whole line under the given timetable can be obtained by

$$W=W_{\mathrm{up}}+W_{\mathrm{down}}$$

(12)

Figure 2.

Calculation of down-bound PTTT.

Determining number of train stocks required

In practice, rolling stock assignment needs to be conducted after train timetabling so that the timetable can be actually implemented. Here, a valid rolling stock assignment procedure is proposed to determine the number of train stocks required for a specific train timetable.

First, an n-dimensional square matrix $\mathbf{H}=[h_{i,j}{]}_{n\times n}$ is defined to describe the connection among the scheduled TRTs. The elements of this matrix are all binary. $h_{i,j}=1$ represents that the train stock of round trip i is also used in round trip j after turning back at the start terminal or returning to the rail yard while $h_{i,j}=0$ represents the opposite. Besides, the values of elements in this matrix should meet the constraints shown in equations (13) and (14), which specifies that a train stock can only undertake one round trip at a time. As can be seen, rolling stock assignment is to determine a reasonable matrix $\mathbf{H}$

$$\sum _{j=1}^n{h}_{i,j}\le 1$$

(13)

$$\sum _{i=1}^n{h}_{i,j}\le 1$$

(14)

Given a specific train timetable, the start times of all the TRTs $\mathbf{x}=[x_1,x_2,\dots ,x_n]$ are specified, as well as the times the round trips returning to start terminal, that is, $a_{1,1}^{\mathrm{up}},a_{2,1}^{\mathrm{up}},\dots ,a_{n,1}^{\mathrm{up}}$ . Besides, an integer record variable $\lambda$ is introduced to record the number of on-going turning-backs. On this basis, the rolling stock assignment can be conducted as follows:

• Step 0: set $\mathbf{H}=0$ .

• Step 1: determine which round trips’ train stock can undertake a new round trip through turning back at the start terminal.

•  Step 1.1: check whether the train stock of round trip i can conduct turning-back at the start terminal.

•   Step 1.1.1: if $i\prime =i$ , then the train stock of round trip i can conduct turning-back at the start terminal and go to Step 1.2. Otherwise, go to Step 1.1.2.

•   Step 1.1.2: check whether the train stock of round trip $i\prime$ has been assigned to certain subsequent round trip. If $\sum _{j=1}^m{h}_{i\prime ,j}=0$ , then go to Step 1.1.4. Otherwise, extract the TRT $\overline{j}$ that follows $h_{i\prime ,\overline{j}}=1$ and go to Step 1.1.3.

•   Step 1.1.3: if $a_{i,1}^{\mathrm{up}}+d_1^{\mathrm{up}}<x_{j^{*}}-d_1^{\mathrm{down}}$ , set $\lambda =\lambda +1$ .

•   Step 1.1.4: if $\lambda <2$ , then set $i\prime =i\prime +1$ and return to Step 1.1.1. Otherwise, the train stock of round trip i cannot conduct turning-back at the start terminal and go to Step 1.3.

Step 1.2: determine which subsequent round trip this train stock undertakes after turning back at the start terminal. Check whether there is TRT starting from the start terminal after the time $a_{i,1}^{\mathrm{up}}+{\delta }_1^{\mathrm{up}}+{\delta }_{\mathrm{start}}$ . If yes, then extract the first one $\tilde{j}$ and set $h_{i,\tilde{j}}=1$ . Otherwise, go to Step 2.

Step 1.3: if $i<n$ , then set $i=i+1$ and return to Step 1. Otherwise, go to Step 2.

Step 2: determine which round trips can be undertaken by the train stock returning to the rail yard from certain previous round trip.

 Step 2.0: set $i=n$ .

 Step 2.1: if $\sum _{j=1}^m{h}_{j,i}=0$ , then go to Step 2.2. Otherwise, go to Step 2.3.

 Step 2.2: if there is TRT j returning to start terminal before the time $x_i-{\delta }_1^{\mathrm{down}}-{\delta }_1^{\mathrm{up}}-\Delta$ and the corresponding train stocks are not used afterward, that is, $\sum _{i=1}^m{h}_{j,i}=0$ , extract the last one $\hat{j}$ and set $h_{\hat{j},i}=1$ . Otherwise, go to Step 3.

 Step 2.3: if $i>1$ , then set $i=i-1$ and return to Step 2.1. Otherwise, go to Step 3.

Step 3: extract the columns whose elements are all 0 from the matrix $\mathbf{H}$ . The number of these columns is the required number of train stocks for the given train timetable.

Model formulation

On the basis of the theories proposed in section “Performance assessment of urban rail train timetable,” the multi-objective optimization problem proposed in section “Problem statement” can be formulated. In this section, the objective functions and constraints of the mathematical model are designed.

Objective function

Given a specific train timetable, the corresponding PTTT can be calculated using the iteration algorithm proposed in section “Calculating PTTT” and the number of train stocks required can be determined through the rolling stock assignment procedure proposed in section “Determining number of train stocks required.” Here, in order to facilitate description, these two procedures are abstracted as two functions

$$W=f(\mathbf{x})$$

(15)

$$N=g(\mathbf{x})$$

(16)

where $\mathbf{x}$ is an n-dimensional vector composed of the starting times of all TRTs in the train timetable, that is, $\mathbf{x}=[x_1,x_2,\dots ,x_n]$ .’W is the PTTT and N is the number of train stocks required for implementing the timetable. Then the objective function of the multi-objective optimization model can be formulated as

$$min\mathbf{Z}(\mathbf{x})={[f(\mathbf{x}),g(\mathbf{x})]}^{\mathrm{T}}$$

(17)

Constraints

Three kinds of constraints are considered in this model, which are formulated, respectively, in the following.

Operation time constraints

In practice, the operation time of urban rail line is always clearly defined and all the transit services should be conducted within the specified time horizon. In this study, therefore, the start times of TRTs are all limited predefined time horizon, that is

$$x_i\in [T_s,T_e],\forall i=1,2,\dots ,n$$

(18)

where $T_s$ and $T_e$ are the start time and end time of the scheduled time horizon, respectively.

Safety headway constraints

As analyzed in section “Problem statement,” the start times of TRTs should ensure no conflict happens in the following train operations. Therefore, the headway between the start times of two consecutive round trips should be larger than the minimum safety headway between the two consecutive running trains in the section, two consecutive arrivals at each station, and two consecutive turning-backs at the return terminal. Synthesizing these safety requirements, a minimum safety headway $\tau$ is given to restrict the start times of round trips, that is

$$x_{i+1}-x_i\ge \tau ,\forall i=1,2,\dots ,n-1$$

(19)

Unmet demand constraint

When calculating the PTTT using the iteration algorithm proposed in section “Calculating PTTT,” the unmet travel demands $\tilde{\mathbf{R}}=[{\tilde{r}}_{u,v}^k{]}_{m\times m\times p}$ can also be obtained, which can be abstracted as a function, that is

$$\tilde{\mathbf{R}}=\tilde{f}(\mathbf{x})$$

(20)

Generally, unmet travel demand has an adverse effect on the transit service quality of urban rail transit, which should be avoided as much as possible. Therefore, in order to guarantee the service quality, the amount of unmet demands should be limited in a reasonable range, that is

$$\Gamma (\tilde{\mathbf{R}})=\sum _{k=1}^K\sum _{u=1}^m\sum _{v=1}^m{\tilde{r}}_{u,v}^k\le \Re$$

(21)

where $\Gamma ()$ is a function to calculate the summation of a matrix’s elements, and ℜ is the maximum allowable amount of unmet travel demand.

Pareto optimal solution using generalized PSO

As can be seen, the problem we study here is a quite complex combinatorial optimization problem, which is usually solved by heuristic algorithms.^15–19 Besides, as analyzed in section “Introduction,” many existing studies formulate the urban rail train scheduling as multi-objective optimization and most of them choose weighted aggregation method to handle the multiple objectives.^14,15,17 In this study, Pareto optimization, which is considered more reasonable in handling contradictive objectives, is introduced to solve the proposed multi-objective train scheduling problem and a generalized PSO procedure is applied to achieve this solution.

Model adjustment and Pareto solution

In order to facilitate the PSO operation, the model proposed in section “Model formulation” is adjusted. The unmet demand constraint is introduced into the first objective with a large number M as penalty coefficient, that is

$$F(\mathbf{x})=f(\mathbf{x})+Mmax\{\Gamma (\tilde{f}(\mathbf{x}))-\Re ,0\}$$

(22)

Then the whole model can be described as

$$min\mathbf{Z}(\mathbf{x})={[F(\mathbf{x}),g(\mathbf{x})]}^{\mathrm{T}}$$

(23)

$$\mathbf{x}\in X=\{\mathbf{x}\in {\mathbb{R}}^n|\mathrm{Constraints}(18)\mathrm{and}(19)\}$$

(24)

where $X$ is the set of all feasible solutions which is determined by constraints (18) and (19).

According to the general theory of Pareto optimality, for two arbitrary solutions $\mathbf{x},\mathbf{x}\prime \in X$ , if

$$\{\begin{array}{c}F(\mathbf{x})\le F(\mathbf{x}\prime )\hfill \\ \multicolumn{1}{c}{g(\mathbf{x})<g(\mathbf{x}\prime )}\end{array}\mathrm{or}\{\begin{array}{c}F(\mathbf{x})<F(\mathbf{x}\prime )\hfill \\ \multicolumn{1}{c}{g(\mathbf{x})\le g(\mathbf{x}\prime )}\end{array}$$

it can be said that the solution $\mathbf{x}$ dominates the solution $\mathbf{x}\prime$ . If a solution ${\mathbf{x}}^{*}$ is not dominated by any other feasible solution, then it is called a non-dominated solution. All the non-dominated solutions in the solution space constitute the non-dominated set and their corresponding objective function values constitute the Pareto front. Therefore, from the perspective of Pareto optimality, solving multi-objective optimization problem is to extract a proper non-dominated set from the solution space.

Pareto-based PSO

PSO is a widely used population-based stochastic search algorithm, which is proposed by Kennedy and Eberhart. ²⁰ In PSO, each particle is endowed with three attributes—position, fitness, and velocity. Generally, the position of a particle corresponds to a particular solution in the solution space and the fitness corresponds to the objective function value of this solution. During iteration, the position of each particle is updated constantly, through which the solution space is searched. In classic PSO, the updates of particles’ positions are typically controlled by their velocities, that is

$${\mathbf{pos}}_p^{(I)}={\mathbf{pos}}_p^{(I-1)}+{\mathbf{vel}}_p^{(I-1)}$$

(25)

where I is the iteration index and p is the particle index. ${\mathbf{vel}}_p^{(I)}$ and ${\mathbf{pos}}_p^{(I)}$ represent, respectively, the velocity and position of particle p at iteration I.

In addition, the velocities of the particles are also iteration-varying. On each iteration, the velocity of each particle needs to be adjusted according to some visited positions that achieve a good fitness, that is

$$ve{l}_p^{(I)}=\omega ve{l}_p^{(I-1)}+{\mu }_1{R}_1\left(po{s}_p^{\mathrm{best}}-po{s}_p^{(I-1)}\right)+{\mu }_2{R}_2\left(po{s}_{\mathrm{global}}^{\mathrm{best}}-po{s}_s^{(I-1)}\right)$$

(26)

where ${\mathbf{pos}}_p^{\mathrm{best}}$ is the personal best position of particle p and ${\mathbf{pos}}_{\mathrm{global}}^{\mathrm{best}}$ is the global best position. $\omega$ is the inertia weight coefficient. ${\mu }_1$ and ${\mu }_2$ are the acceleration constants, which are typically set to be 2. $R_1$ and $R_2$ are the random numbers uniformly distributed over $[0,1]$ , that is, $R_1,R_2~U[0,1]$ .

Compared with other evolutionary algorithms, PSO achieves considerable computational accuracy with a much simplifier procedure and faster convergence rate. ²¹ Since being proposed, PSO has been widely used in various research areas, especially in traffic and transportation.^22–26 In this process, the classic PSO has often been extended to accommodate different situations. One successful example of these extensions is the multi-objective PSO, that is, MOPSO. ²⁷ In the past decade, a number of studies have been conducted on MOPSO and many useful theories have been proposed. Based on these existing theories, a Pareto MOPSO is designed to solve the proposed multi-objective train timetabling problem.

Main algorithm

Compared with single-objective PSO, MOPSO has some significant differences, one of which is the determination of ${\mathbf{pos}}_p^{\mathrm{best}}$ and ${\mathbf{pos}}_{\mathrm{global}}^{\mathrm{best}}$ . In MOPSO, non-dominated solutions can be obtained on each iteration, but the performance of these non-dominated solutions cannot be compared directly. A common approach is to collect all the non-dominated solutions and determine the best one according to their distribution in the objective space. In order to achieve that, an external archive is usually introduced to store the obtained non-dominated solutions. On this basis, the general procedure of MOPSO is as follows:

• Step 0: set $I=0$ , initialize the particle swarm.

• Step 1: calculate the objective vectors of each particle. On this basis, update the external archive.

• Step 2: update ${\mathbf{pos}}_p^{\mathrm{best}}$ for each particle and determine global best position ${\mathbf{pos}}_{\mathrm{global}}^{\mathrm{best}}$ .

• Step 3: if $I=I_{max}$ , stop and output the current external archive. Otherwise, continue.

• Step 4: set $I=I+1$ , update the velocity and position of each particle and return to Step 1.

The key points in above procedure include the following:

1. Initializing the particle swarm;

2. Updating the external archive;

3. Determining each particle’s personal best position and the global best position;

4. Updating velocity and position of each particle.

In the following parts, with the existing related theories, these key points are reasonably designed according to the features of the studied problem.

Swarm initialization

As can be seen, the main work during initial swarm generation is to set initial position and velocity for each particle, which influences efficiency of the whole algorithm to a large extent. Here, to solve the proposed model, the initial position and the velocity of the particles are generated in the following way.

Initial position

For each particle, the initial position is obtained through the following procedure:

• Step 0: set $i=0$ .

• Step 1: set $i=i+1$ .

• Step 2: if $i\le n$ , then assign a random integer to $\mathrm{po}{s}_i$ according to

$$\mathrm{po}{s}_j=\{\begin{array}{cc}\mathrm{rand}[T_s,T_e-\tau (n-1)]\hfill & i=1\hfill \\ \multicolumn{1}{c}{\mathrm{rand}[\mathrm{po}{s}_{i-1}+\tau ,T_e-\tau (n-i)]}& 1<i\le n\hfill \end{array}$$

(27)

where rand[] is a function which returns a random integer within the specific interval, and go back to Step 1. Otherwise, stop.

Initial velocity

In order to guarantee the feasibility of particles’ position after first iteration, the initial velocity of each particle is generated as follows:

• Step 0: set $i=0$ .

• Step 1: set $i=i+1$ .

• Step 2: if $i\le n$ , then assign a random integer to $\mathrm{ve}{l}_i$ according to

$$\mathrm{ve}{l}_i=\{\begin{array}{cc}\mathrm{rand}[-\mathrm{po}{s}_i,T_e-\tau (n-1)-\mathrm{po}{s}_i]\hfill & i=1\hfill \\ \multicolumn{1}{c}{\mathrm{rand}[\mathrm{po}{s}_{i-1}+\mathrm{ve}{l}_{i-1}+\tau -\mathrm{po}{s}_i,T_e-\tau (n-i)-\mathrm{po}{s}_i]}& 1<i\le n\hfill \end{array}$$

(28)

and go back to Step 1. Otherwise, stop.

Archive maintenance

In order to reduce the computation cost, the size of external archive is typically limited. When a new non-dominated solution is obtained, the archive is updated through the following procedure:

• Step 1: check whether these solutions in the archive are dominated by the new obtained solution. If yes, remove them and insert this new solution into the archive. If no, continue.

• Step 2: check whether the size of archive reaches the ceiling. If no, then directly insert this new solution into the archive. If yes, continue.

• Step 3: compare the performances of solutions in the archive and the new solution. If the worst solution in the archive is still better than the new one, then drop the new solution. Otherwise, replace the worst solution with the new obtained one.

To achieve this procedure, a definite approach for assessing the performance of non-dominated solution must be given. In our study, considering the penalty item in the first objective, a modified density-based approach is proposed. Given a set of non-dominated solutions of the studied problem ${\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_S$ , the fitness of each solution can be calculated through the following procedure:

• Step 1: select the solutions that do not satisfy constraint (21), set their fitness as $-F({\mathbf{x}}_s)$ , and remove them from the set.

• Step 2: for the remaining solutions, normalize their coordinates in the objective space according to

$$\hat{F}({\mathbf{x}}_s)=\frac{F({\mathbf{x}}_s)-\underset{s\prime }{min}\{F({\mathbf{x}}_{s\prime })\}}{\underset{s\prime }{max}\{F({\mathbf{x}}_{s\prime })\}-\underset{s\prime }{min}\{F({\mathbf{x}}_{s\prime })\}}$$

(29)

$$\hat{g}({\mathbf{x}}_s)=\frac{g({\mathbf{x}}_s)-\underset{s\prime }{min}\{g({\mathbf{x}}_{s\prime })\}}{\underset{s\prime }{max}\{g({\mathbf{x}}_{s\prime })\}-\underset{s\prime }{min}\{g({\mathbf{x}}_{s\prime })\}}$$

(30)

where $[\hat{F}({\mathbf{x}}_s),\hat{g}({\mathbf{x}}_s)]$ is the normalized coordinates of solution ${\mathbf{x}}_s$ .

• Step 3: for each remaining solution ${\mathbf{x}}_s$ , find its nearest solution ${\mathbf{x}}_{s\prime }$ according to the normalized coordinates and set its fitness as $|\hat{F}({\mathbf{x}}_s)-\hat{F}(v)|+|\hat{g}({\mathbf{x}}_s)-\hat{g}({\mathbf{x}}_{s\prime })|$ .

Determining personal and global best position

On each iteration, the personal best position of each particle is determined simply according to Pareto dominance, that is, if the new position ${\mathbf{pos}}_p^{(I)}$ is dominated by ${\mathbf{pos}}_p^{\mathrm{best}}$ , then ${\mathbf{pos}}_p^{\mathrm{best}}$ is kept; if ${\mathbf{pos}}_p^{(I)}$ dominates current ${\mathbf{pos}}_p^{\mathrm{best}}$ , then set ${\mathbf{pos}}_p^{\mathrm{best}}={\mathbf{pos}}_p^{(I)}$ ; if neither of them is dominated by the other, then select one of them randomly to be the new ${\mathbf{pos}}_p^{\mathrm{best}}$ . ²⁸

As for global best position, the fitness calculation procedure proposed in the last section needs to be conducted on the archive again, and the non-dominated solution with maximum fitness is chosen to be the new ${\mathbf{pos}}_{\mathrm{global}}^{\mathrm{best}}$ .

Updating position and velocity

On each iteration, the position and velocity of each particle are updated according to equations (25) and (26). Here, the control parameters in equation (26) are set as follows:

• 1. The inertia weight coefficient $\omega$ is set to follow an S-type variation during iteration, as formulated in equation (29)

$$\omega =\frac{{\omega }_{max}-{\omega }_{min}}{1+{\vartheta }^{(2\vartheta \times I/I_{max}-\vartheta )}}+{\omega }_{min}$$

(31)

where ${\omega }_{max}$ and ${\omega }_{min}$ are the maximum and minimum inertia weights, respectively. $I$ and $I_{max}$ are the current and maximum iterations, respectively. $\vartheta$ is the speed variation coefficient and set to be 10.

• 2. The acceleration constants ${\mu }_1$ and ${\mu }_2$ are set to be adaptive according to an evolution factor. For each particle $p$ , the summation of Euclidean distances between it and other particles $S{D}_p$ is calculated as

$$S{D}_p=\frac{1}{P-1}\sqrt{\sum _{p\prime =1}^P{(\mathbf{po}{\mathbf{s}}_p-\mathbf{po}{\mathbf{s}}_{p\prime })}^2}$$

(32)

Then the evolutionary factor is calculated as

$$\mathrm{EF}=\frac{S{D}_{p^{*}}-\underset{p}{min}\{S{D}_p\}}{\underset{p}{max}\{S{D}_p\}-\underset{p}{min}\{S{D}_p\}}$$

(33)

where $p^{*}$ represents the global best position. With the evolutionary factor, the acceleration constants can be determined as

$${\mu }_1=\frac{1}{0.5+e^{-2f}}$$

(34)

$${\mu }_2=\frac{1}{0.2+0.3e^{-f}}$$

(35)

Furthermore, in order to guarantee the feasibility of particles’ positions in next iteration, the velocity of each particle needs to be examined and revised if not feasible. Relative procedure is as follows:

• Step 0: set $i=1$ .

• Step 1: check whether $\mathrm{ve}{l}_i+\mathrm{po}{s}_i+\tau \le \mathrm{ve}{l}_{i+1}+\mathrm{po}{s}_{i+1}$ . If no, set $\mathrm{ve}{l}_i=\mathrm{ve}{l}_{i+1}+\mathrm{po}{s}_{i+1}-\mathrm{po}{s}_i-\tau$ .

• Step 2: set $i=i+1$ . If $i=n$ , then stop. Otherwise, return to Step 1.

Numerical example

Small-scaled example

First, the proposed methodology is tested on a fictitious small-scaled urban rail line, whose layout is shown in Figure 3. As can be seen, the line contains four stations. The numbers labeled between two adjacent stations specify the running time of train in this section along the corresponding direction, and the numbers labeled at the start and return terminal represent the required times for turning-back at these two stations. For this example line, a 4-hour OD travel demand data are collected, which are not listed here due to the limited space. Besides, on this specific line, the train dwell time at each station is 2 min, the minimum safety headway between the start times of two consecutive TRTs, that is, $\tau$ , is 3 min, and the minimum time a train spends from leaving the start terminal for the yard to arriving at the start terminal again form the yard is 10 min. All the trains have a capacity of 600 persons.

Figure 3.

Layout of the small-scaled example urban rail line.

Under a MATLAB 9.0 environment, the proposed methodology is implemented to conduct Pareto optimal train scheduling for this specific line. Based on that, the designed algorithm is tested.

First, different control parameter settings are compared. For the inertia weight coefficient $\omega$ , the proposed adaptive method is compared with three different fixed values, that is, $\omega =1.2,1.0,0.8$ . The satisfaction of unmet demand constraint during iteration is first compared, as shown in Figure 4(a). As can be seen, using the adaptive method, all the infeasible solutions are eliminated after 17 iterations, which are obviously faster than $\omega$ fixed at 1.0 and 1.2 but slightly slower than $\omega =0.8$ . Besides, the number of non-dominated solutions during iteration is also compared, as shown in Figure 4(b). We can see, with the adaptive $\omega$ , the algorithm can achieve a faster convergence than the other three situations.

Figure 4.

Comparison between different settings of inertia weight coefficient: (a) Variation of the number of infeasible solutions and (b) Variation of the number of non-dominated solutions.

For the acceleration constants ${\mu }_1$ and ${\mu }_2$ , the proposed adaptive method is compared with a recommended valuation scheme, that is, ${\mu }_1={\mu }_2=2$ , as shown in Figure 5. Although the adaptive acceleration constants do not show much advantage on removing the infeasible solutions during iteration, it provides an obvious improvement on the convergence speed.

Figure 5.

Comparison between different settings of acceleration constants: (a) Variation of the number of infeasible solutions and (b) Variation of the number of non-dominated solutions.

Moreover, based on the optimization toolbox provided by MATLAB, this example is further calculated by a branch and bound (BAB) procedure, and the result is compared with that of the proposed algorithm, as shown in Figure 6. As can be seen, the Pareto front obtained by PSO is almost as good as the result obtained by BAB. However, the PSO is much more efficient than BAB. When solving this small urban rail scheduling problem, PSO spends about 14 min while BAB spends more than 55 min. Therefore, the designed PSO procedure can achieve an adequate train scheduling with more reasonable time consumption.

Figure 6.

Comparison between PSO and BAB.

Practical-sized example

Setup

Using the proposed methodology, train scheduling for a practical-sized urban rail line is conducted. This line contains 10 stations. Its layout, section running times, and turning-back times are shown in Figure 7, and the corresponding train dwell times are listed in Table 1. For this rail line, the OD demand statistics are collected for every 30 min of a whole day. Due to the limited space, only the distribution of passenger arrivals and the total cross-sectional passenger flow volume in each section are presented here, as shown in Figures 8 and 9. Besides, the minimum safety headway between the start times of two consecutive TRTs, that is, $\tau$ , is 3 min and the minimum time a train spends from leaving the start terminal for the yard to arriving at the start terminal again form the yard is 30 min. All the trains have a capacity of 1500 persons.

Figure 7.

Layout of the practical-sized example urban rail line.

Table 1.

Train dwell time at each station on the second example line.

Station	Dwell time		Station	Dwell time
	Up-bound	Down-bound		Up-bound	Down-bound
1	1	1	6	2	2
2	1	1	7	1	2
3	2	1	8	1	1
4	1	1	9	1	1
5	2	2	10	1	1

Figure 8.

Distribution of passenger arrivals on the second example line.

Figure 9.

Total cross-sectional passenger flow volume in each section on the second example line.

Pareto optimal train scheduling

With the TRT number defined as 90, train timetabling for the example rail line is conducted using the designed MOPSO. The related parameters are set as shown in Table 2. After 100 iterations, a reasonable Pareto front is obtained as shown in Figure 10.

Table 2.

Parameter setting.

Parameter	P	$I_{max}$	ℜ	M	A
Value	60	100	100	500	80

Figure 10.

Pareto front.

As can been seen, as the number of train stocks is small ( $L<10$ ), the PTTT is rather large and deceases fast with the increase in $L$ , which implies that at this stage the number of train stocks is the main factor limiting the train scheduling. As the number of train stocks increases, the descent rate of PTTT gradually reduces. Especially when $L$ exceeds 16, the PTTT almost becomes constant. Considering the emphasis people put on service quality in practical operation, the solution whose $L=16$ is chosen as the best train timetable, as shown in Table 3.

Table 3.

Pareto optimal train timetable.

TRT	Start time	TRT	Start time	TRT	Start time	TRT	Start time	TRT	Start time
1	5:48	19	8:12	37	11:51	55	16:51	73	18:30
2	6:07	20	8:18	38	12:09	56	17:02	74	18:37
3	6:24	21	8:25	39	12:26	57	17:12	75	18:47
4	6:42	22	8:34	40	12:44	58	17:18	76	18:55
5	6:58	23	8:43	41	13:02	59	17:25	77	19:07
6	7:10	24	8:53	42	13:21	60	17:32	78	19:18
7	7:20	25	9:03	43	13:39	61	17:40	79	19:31
8	7:27	26	9:11	44	13:56	62	17:46	80	19:44
9	7:31	27	9:19	45	14:15	63	17:49	81	19:59
10	7:35	28	9:33	46	14:34	64	17:53	82	20:14
11	7:41	29	9:48	47	14:50	65	17:56	83	20:31
12	7:44	30	10:03	48	15:06	66	17:59	84	20:46
13	7:48	31	10:18	49	15:22	67	18:03	85	21:01
14	7:52	32	10:33	50	15:38	68	18:07	86	21:16
15	7:56	33	10:48	51	15:53	69	18:11	87	21:37
16	8:00	34	11:04	52	16:07	70	18:16	88	21:58
17	8:04	35	11:19	53	16:23	71	18:21	89	22:19
18	8:08	36	11:33	54	16:38	72	18:26	90	22:39

TRT: train round trip.

Conclusion and discussions

In this article, a Pareto optimal urban rail train scheduling problem is studied. Based on a representative layout of urban rail line, the studied problem is well defined as bi-objective optimization. Two innovative theories are proposed: an effective iteration algorithm to calculate PTTT and a valid train stock assignment procedure to obtain the number of used train stocks. Based on these theories, the problem is reasonably formulated. Then, based on general framework of MOPSO, a specific MOPSO is designed for the studied problem. Through the numerical experiment on two different scaled examples, the effectiveness of the designed algorithm is illustrated.

As for further research, several extensions can be considered. First of all, to further improve the practicability of obtained timetable, the dwell times of trains at stations should be determined according to the corresponding alighting and boarding process instead of predefined. Second, in the initial swarm generation of PSO, the unmet demand constraint can be further taken into account, so that the procedure can be simplified and the convergence rate of the whole algorithm can be improved. Third, the methodology should be tested under real context to verify its effectiveness in solving practical urban rail train scheduling problem.