A timetable optimization model and an improved artificial bee colony algorithm for maximizing regenerative energy utilization in a subway system

Abstract

Maximizing regenerative energy utilization in subway systems has become a hot research topic in recent years. By coordinating traction and braking trains in a substation, regenerative energy is optimally utilized and thus energy consumption from the substation can be reduced. This article proposes a timetable optimization problem to maximize regenerative energy utilization in a subway system with headway and dwell time control. We formulate its mathematical model, and some required constraints are considered in the model. To keep the operation time duration constant, the headway time between different trains can be different. An improved artificial bee colony algorithm is designed to solve the problem. Its main procedure and some related tasks are presented. Numerical experiments based on the data from a subway line in China are conducted, and improved artificial bee colony is compared with a genetic algorithm. Experimental results prove the correctness of the mathematical model and the effectiveness of improved artificial bee colony, which improves regenerative energy utilization for the experimental line and performs better than genetic algorithm.

Keywords

Regenerative energy utilization timetable optimization headway time dwell time artificial bee colony

Introduction

Although subway systems are considered to be energy-efficient compared with other transportation, for example, private cars or buses, their energy consumption is very high. The studies on railway energy conservation have attracted many interests.¹ Around 40% of the total energy consumed in a subway system is consumed by train traction systems.^2–4 Therefore, many researchers are devoting themselves to reduce the traction energy consumption from the substations. Benefit from the development and deployment of regenerative braking systems, trains’ kinetic energy can be converted to electrical energy during braking phases.⁴ The electrical energy recovered from braking trains is called regenerative energy (RE). It can go back to power supply line, thus it can be utilized by accelerating trains and/or other electric equipment (e.g. air conditioners, lighting facilities) to reduce energy consumption in the same substation. As the energy demands of the other electric equipment are relatively small, there are two ways to utilize RE in general. One is to accelerate trains in the same substation immediately, and the other one is to store it in energy storage devices (e.g. super-capacitors) for later use. Since energy storage devices usually cost high and have limited capacities,⁵ maximizing regenerative energy utilization (REU) through the synchronization of traction and braking trains is definitely a preferential way.^2,6 Note that if RE is not utilized sufficiently and immediately, its surplus would be wasted by resistors into heat. The coordination of traction and braking trains can be fulfilled by slightly adjusting the timetable used in a subway system. It costs little and the impacts to the quality of subway service are limited. This explains why the studies on maximizing REU through timetable optimization have attracted great attention in recent years.^7–13 The objective is to make the braking and traction of different trains occur simultaneously such that RE can be optimally utilized in a subway system.

Nag and Pal¹⁴ first present a timetable model with the consideration of RE. Ramos et al.⁷ present a mixed integer optimization problem for timetabling during off-peak hours, they take the overlap time of the traction and braking phases of the trains at the same substation as the measurement of REU, and solved the problem using a CPLEX software. Nasri et al.⁸ use a simulation method to study REU, and the impacts of headway time between trains and reserve time at stations on REU are studied. Kim et al.⁹ optimize the departure time of trains at the starting station to maximize REU. A multi-criteria mixed integer programming (MIP) method is used in their work. Fournier et al.¹¹ develop an optimization model to maximize REU by subtly modifying the dwell time for trains at stations. A hybrid genetic/linear programming algorithm is implemented to solve their problem. Li and Lo¹⁵ present an integrated energy-efficient operation model to jointly optimize the timetable and speed profiles at sections. A genetic algorithm (GA) is designed to solve their problem.

Yang et al.¹² present a cooperative scheduling model to coordinate the traction and braking phases of successive trains in the same operation direction. They assume that RE from braking trains can only be utilized by its adjacent traction trains in the same direction. GA is designed to solve their problem. Both headway time and dwell time are optimized in their work. Yang et al.¹³ present a two-objective timetable optimization model to maximize REU and reduce passenger waiting time concurrently. Their model coordinates the traction and braking phases of trains in the opposite directions at the same station, that is, one train is departing from a platform, while the other train is arriving at the opposite platform at the same station. GA is designed to solve their problem with the decision variables as the headway and dwell time. Their work lay solid foundations for the timetable optimization models of maximizing REU. However, several improvements can be made from their work. For example, the surplus RE not utilized by the adjacent trains is wasted in their model, rather than being utilized by other traction trains in the same substation, which disaccords with the facts of REU. In addition, the headway time between any two successive trains is assumed to be identical and got decreased in the experiments. However, shortened headway time leads to either shortened operation time or increased number of train trips. Both of them are not desired in a subway system in fact.

Motivated by the previous studies, we propose a timetable optimization problem to maximize REU for a subway system with more realistic considerations. RE generated from braking trains is utilized by all traction trains in the same substation, and the constraint that the operation time should be constant is considered in this work, under the assumption that the number of train trips is fixed. To make the problem solvable, headway time between different successive train pairs are not required to be identical. We formulate a mathematical model of the timetable optimization problem and then design an improved artificial bee colony (IABC) algorithm to solve the problem. The main procedures of IABC are presented, and the generation/update of a solution is designed according to the characteristics of the problem. Numerical studies based on the operation data from a subway line in China are conducted to illustrate the correctness of the mathematical model and the effectiveness of IABC. In addition, IABC is also compared with and performs better than GA.

The rest of the article is organized as follows. We present the timetable optimization problem and its mathematical model in section “Timetable optimization mathematical model.” Section “IABC algorithm” introduces the principle of the proposed IABC. We present the experimental studies in section “Experimental results and analysis” and conclude this article in section “Conclusion.”

Timetable optimization mathematical model

For a better understanding of this article, the parameters and variables used are introduced first.

Parameters and variables

I: number of train trips.

$i, j$ : train trip indices, $i, j \in {1, 2, \dots, I}$ .

N: number of stations.

$n, k$ : station/platform/section indices, $n, k \in {1, 2, \dots, N}$ for stations, $n, k \in {1, 2, \dots, 2 N - 1}$ for platforms and $n, k \in {1, 2, \dots, 2 N - 2}$ for sections, respectively.

m: train mass.

$γ$ : train running phase index, $γ \in {1, 2, 3}$ where 1 to 3 represents a traction, coasting, and braking phase, respectively.

$a_{n}^{γ}$ : train acceleration in the γth phase at section n.

$r_{n}^{γ}$ : time duration of the γth phase at section n.

$r_{n}$ : total running time duration at section n, which is the time duration between a train departs from platform n and arrives at platform n + 1.

u : switch time point index of the running process at a section, $u \in {1, 2, 3, 4}$ where 1 to 4 represents the start time of the traction phase, the switch time point from the traction phase to coasting, the switch time point from coasting to braking and the end time of the braking phase, respectively.

$t_{i, n}^{u}$ : time instant of the uth switch time point for train i at section n.

$t_{i, 0}^{4}$ : time instant of the start of the ith train trip, which is the time point that train i ends its braking phase at the imaginary section 0.

$v_{i, n}^{u}$ : train speed when train i is at the uth switch time point at section n.

$Φ$ : turnaround time duration at terminal station N.

$τ_{i}$ : trip time duration of train i, which is the time duration between $t_{i, 0}^{4}$ and $t_{i, 2 N - 2}^{4}$ . Note that 2N − 2 is the index of the last platform in our model, which will be discussed in the next section.

$Λ$ : operation time duration.

$\underline{h}, \bar{h}$ : lower and upper limits of headway time, respectively.

${\underline{x}}_{n}, {\bar{x}}_{n}$ : lower and upper limits of dwell time at platform n, respectively.

$\underline{τ}, \bar{τ}$ : lower and upper limits of the train trip time, respectively.

$E_{i, n}^{a} (t)$ : required electrical energy to accelerate train i at section n at time t.

$E_{i, n}^{b} (t)$ : RE generated from train i at section n at time t.

$η_{1}$ : conversion efficiency from traction electricity to kinetic energy.

$η_{2}$ : conversion efficiency from kinetic energy to regenerative electricity.

$ζ$ : transmission loss coefficient of REU.

$f_{i, n} (t)$ : required electrical power to accelerate train i at section n at time t.

$g_{i, n} (t)$ : power of RE generated from train i at section n at time t.

Z : number of substations.

z : substation index, $z \in {1, 2, \dots, Z}$ .

$S_{z}$ : set of sections in the zth substation.

E : total REU in a subway system.

Decision variables

$h_{i}$ : headway time between trains i and $i + 1$ .

$x_{i, n}$ : dwell time for train i at platform n.

Problem statement

As shown in Figure 1, a subway line comprised N stations and the sections connecting them. Every train travels along the line, during which it stops at every platform and passes the sections connecting them. Every train starts its trip from platform 1. After the dwell time for passengers to get on/off, it departs from the platform and runs along the down direction, until it arrives at the next platform. Similar process repeats until it arrives at the terminal station N. Then, it turns around to the up direction with a short time interval. After the dwell time at platform N, it departs and runs along the up direction. Thereafter, the similar process repeats as in the down direction, until it arrives at the last platform 2N − 1, where it ends its trip. Note that for station 1 ≤ n ≤ N − 1, there are two platforms located in the opposite directions of it, labeled as platform n and platform 2N − n, respectively. The terminal station N has only one platform labeled as platform N, as the detailed turnaround process is dismissed in the analysis of REU. The path connects platform n and platform n + 1 is labeled as section n.

Figure 1.

Train trip process.

Every train travels at the same path, but with a time interval between any two successive trains. The time interval between trains i and i + 1 is named the headway time of train i, denoted as $h_{i}$ . A timetable is used in subway systems, which defines all the time instant that each train arrives at and departs from each platform. The headway time, dwell time, trip time, and operation time can all be obtained from a timetable. However, train speed profiles at a section are not specified in a timetable, while it is very important for calculating REU. The train running process at a section is divided into three phases, that is, the traction, coasting, and braking phases, respectively.^12,16 The detailed parameters are illustrated in Figure 2.

Figure 2.

General process of a train running at an inter-section.

For any train i and section n, we have

t_{i, n}^{u + 1} = t_{i, n}^{u} + r_{n}^{γ}

(1)

where $γ \in {1, 2, 3}$ and u = γ. The calculation of the departure time at platform N differs from that at other platforms, as each train turns around at station N, as shown in equation (2)

t_{i, n}^{1} = {\begin{matrix} t_{i, n - 1}^{4} + x_{i, n} & n \in [1, 2 N - 2] / {N} \\ t_{i, n - 1}^{4} + x_{i, n} + Φ & n = N \end{matrix}

(2)

For any train i, $t_{i, 0}^{4}$ can be obtained by $t_{i - 1, 0}^{4}$ and the headway time $h_{i - 1}$ . Recursively, it can be obtained by $t_{1, 0}^{4}$ and the sum of the headway time between any two successive trains before train i. It is described as follows

t_{i, 0}^{4} = t_{i - 1, 0}^{4} + h_{i - 1} = t_{1, 0}^{4} + \sum_{j = 1}^{i - 1} h_{j}

(3)

Operation time duration of a subway system is represented by the time difference between the time instant that the first train starts its trip and that the last train starts its trip in this work, that is, $Λ = t_{I, 0}^{4} - t_{1, 0}^{4}$ . By substituting it into equation (3), we can obtain the following

Λ = t_{I, 0}^{4} - t_{1, 0}^{4} = \sum_{j = 1}^{I - 1} h_{j}

(4)

Trip time duration of a train is determined by its dwell time at platforms, running time at sections and the turnaround time. It can be described as follows

τ_{i} = \sum_{n = 1}^{2 N - 2} x_{i, n} + \sum_{n = 1}^{2 N - 2} r_{n} + Φ

(5)

where $r_{n} = r_{n}^{1} + r_{n}^{2} + r_{n}^{3}$ .

The train speed is 0 when a train stops at a platform, and it can be calculated according to the kinematics equation when the train is running at sections. Thus, the speed profile of a train is as follows

v_{i} (t) = {\begin{matrix} 0 & t \in [t_{i, n - 1}^{4}, t_{i, n}^{1}] \\ a_{n}^{1} (t - t_{i, n}^{1}) & t \in [t_{i, n}^{1}, t_{i, n}^{2}] \\ a_{n}^{1} r_{n}^{1} - a_{n}^{2} (t - t_{i, n}^{2}) & t \in [t_{i, n}^{2}, t_{i, n}^{3}] \\ a_{n}^{3} (t_{i, n}^{4} - t) & t \in [t_{i, n}^{3}, t_{i, n}^{4}] \end{matrix}

(6)

where $a_{n}^{1}$ and $a_{n}^{3}$ are the maximum traction and braking acceleration of the train at section n, respectively. They both are given constants. $a_{n}^{2}$ can be obtained as $a_{n}^{2} = (a_{n}^{1} Δ t_{n}^{1} - a_{n}^{3} Δ t_{n}^{3}) / Δ t_{n}^{2}$ .

RE has to be utilized immediately by other trains in their traction phases; otherwise, it will be wasted into heat. So REU should be realized by coordinating traction and braking trains in the same substation. It is hard to know exactly how RE will spread throughout the power supply line and other trains.¹¹ The RE generated from braking trains can be shared by the traction trains in the same substation in our model. REU for each train is determined by the total amount of its traction energy demands and the amount of the RE available in this substation, that is, the RE transmission lost is assumed to be a constant. Figure 3 illustrates a case of four trains in a substation. The black arrows denote the direction of the electricity current. The sum of RE from the two braking trains is consumed by the other two traction trains together.

Figure 3.

Flow chart of electrical energy.

The required electrical energy for accelerating train i at section n is determined by the kinetic energy increased in the traction phase, which is shown as follows

E_{i, n}^{a} (t) = {\begin{matrix} m ({(v_{i, n} (t))}^{2} - {(v_{n}^{1})}^{2}) / 2 η_{1} & t \in [t_{i, n}^{1}, t_{i, n}^{2}] \\ 0 & otherwise \end{matrix}

(7)

The required electrical power $f_{i, n} (t)$ for accelerating train i at section n at time t is the derivative of the consuming electrical energy and is given as follows

f_{i, n} (t) = \frac{{dE}_{i, n}^{a} (t)}{dt} = {\begin{matrix} ψ (t - t_{i, n}^{1}) & t \in [t_{i, n}^{1}, t_{i, n}^{2}] \\ 0 & otherwise \end{matrix}

(8)

where ψ is a constant and $ψ = m (a_{n}^{1})^{2} / η_{1}$ .

The total electrical power required to accelerate trains at substation z is determined as follows

P_{z}^{a} (t) = \sum_{i = 1}^{I} \sum_{n \in S_{z}} f_{i, n} (t)

(9)

Similarly, the available RE from train i at section n is determined by the kinetic energy decreased in a braking phase, and the power of RE at time t is the derivative of $E_{i, n}^{b} (t)$ . They are determined as follows

E_{i, n}^{b} (t) = {\begin{matrix} χ_{1} ({(v_{n}^{3})}^{2} - {(v_{i, n} (t))}^{2}) & t \in [t_{i, n}^{3}, t_{i, n}^{4}] \\ 0 & otherwise \end{matrix}

(10)

g_{i, n} (t) = \frac{{dE}_{i, n}^{b} (t)}{dt} = {\begin{matrix} χ_{2} (t_{i, n}^{4} - t) & t \in [t_{i, n}^{3}, t_{i, n}^{4}] \\ 0 & otherwise \end{matrix}

(11)

where $χ_{1}$ is a coefficient, which converts kinetic energy into electrical energy, it is a constant satisfying $χ_{1} = m η_{2} (1 - ζ) / 2$ ; $χ_{2}$ is a coefficient, which determines the power of RE at time t, it is obtained from the derivation of $E_{i, n}^{b} (t)$ , it is a constant satisfying $χ_{2} = m (a_{n}^{3})^{2} η_{2} (1 - ζ)$ . Thus, the total power of RE at a substation is determined as follows

P_{z}^{b} (t) = \sum_{i = 1}^{I} \sum_{n \in A_{z}} g_{i, n} (t)

(12)

REU at the zth substation is determined by the minimal value of $P_{z}^{a} (t)$ and $P_{z}^{b} (t)$ and the operation time duration as follows

E_{z} = \int_{t = t_{1, 0}^{4}}^{t_{I, 2 N - 2}^{4}} min {P_{z}^{a} (t), P_{z}^{b} (t)} dt

(13)

Thus, total REU in a subway system is determined as follows

E = \sum_{z = 1}^{Z} E_{z}

(14)

Mathematical model

The objective of our model is to maximize the total REU in a subway system. The decision variables are headway time and dwell time. The constraints include the safety, cost, and service quality of a subway system. The timetable optimization model is shown as follows, where $H = (h_{1}, h_{2}, \dots, h_{i}, \dots, h_{I - 1})$ and $X = (X_{1}, X_{2}, \dots, X_{i}, \dots, X_{I})$ , where $X_{i} = (x_{i, 1}, x_{i, 2}, \dots, x_{i, 2 N - 2}), \forall i \in {1, 2, \dots, I} .$

max E = \sum_{z = 1}^{Z} E_{z} = f (H, X)

(15)

s.t.

\underline{h} \leq h_{i} \leq \bar{h}, i \in {1, 2, \dots, I}

(16)

{\underline{x}}_{n} \leq x_{i, n} \leq {\bar{x}}_{n}, i \in {1, 2, \dots, I}, n \in {1, 2, \dots, 2 N - 2}

(17)

\sum_{i = 1}^{I - 1} h_{i} = Δ P = C_{1}

(18)

Δ \underline{T} \leq \sum_{n = 1}^{2 N - 2} x_{i, n} + C_{2} \leq Δ \bar{T}, i \in {1, 2, \dots, I}

(19)

h_{i} \in Z^{+}, i \in {1, 2, \dots, I - 1}

(20)

x_{i, n} \in Z^{+}, i \in {1, 2, \dots, I}, n \in {1, 2, \dots, 2 N - 2}

(21)

Objective (equation (15)) aims to maximize the total REU in a subway system. Constraint (equation (16)) guarantees that the headway time between any two successive trains is within a predefined time window, where $\underline{h}$ is a constant determined by the signal system and the operation cost, and $\bar{h}$ is a constant determined by the service quality. equation (17) guarantees that there is enough time for passengers to get on and off a train at each platform and the dwell time is not too long, where ${\underline{x}}_{n}$ and ${\bar{x}}_{n}$ are constants determined by the passenger demands at platform n and they are the same for all trains. Note that the difference between $x_{i, n}$ and ${\bar{x}}_{n}$ is called the reserve time, which is not less than 0 and not larger than ${\bar{x}}_{n} - {\underline{x}}_{n}$ , that is, $0 \leq x_{i, n} - {\underline{x}}_{n} \leq {\bar{x}}_{n} - {\underline{x}}_{n}, \forall i \in (1, 2, \dots, I}$ . The upper limit of the reserve time at each platform is usually the same in reality, that is, ${\bar{x}}_{n} - {\underline{x}}_{n} = {\bar{x}}_{k} - {\underline{x}}_{k}, \forall k, n \in {1, 2, \dots, 2 N - 2}$ . Equation (18) guarantees that operation time duration of a subway system is a constant, where $C_{1} = t_{I, 0}^{4} - t_{1, 0}^{4}$ . Equation (19) guarantees that the trip time of each train is within a required time window of the planned timetable, as the trip time is an important measurement of the service quality of a subway system, where $C_{2}$ is a constant and $C_{2} = \sum_{n = 1}^{2 N - 2} r_{n} + Φ$ . Equations (20) and (21) ensure that the headway time and dwell time are positive integers.

IABC algorithm

Artificial bee colony (ABC) algorithm is a swarm intelligence algorithm. It was first proposed by Karaboga¹⁷ in 2005. Since then, it has been successfully used to handle many complicated optimization problems.^18–22 Its performance has been proven to be better than or similar to other intelligent algorithms, including differential evolution (DE),^23,24 GA,^19,25,26 particle swarm optimization (PSO),^27,28 and evolutionary algorithm (EA).^29,30,31 ABC can be efficiently employed to solve engineering problems with high dimensionality.^32,33

Motivated by the previous studies, we design an IABC algorithm to solve our timetable optimization problem. Its main improvements over the prior ABC algorithms^17,32,33 are as follows. First, a restart mechanism is introduced when the algorithm reaches a local optimum and thus enhances its global search ability. Also, the best solution found in each iteration is kept as one solution in the next iteration to avoid the results getting worse. Second, the current timetable is initialized as one solution, and specific random generation processes are specially designed for the initialization of other solutions. The processes are also used by each scout bee to make the generated solutions feasible, that is, satisfying constraints (16)–(21). Finally, for each employed bee and looker bee, a local search operator is randomly chosen from the swap, insertion, crossover, and mutation operators with the probability of $p_{s}$ , $p_{i}$ , $p_{c}$ , and $p_{m}$ , respectively, and repair is immediately applied to the newly generated solution to ensure that it is feasible.

The main tasks of IABC include encoding of solutions, initialization of solutions, and the optimization with employed bee phase, onlooker bee phase, and scout bee phase, respectively.

Encoding of solutions

Encoding of a solution directly affects the complexity of the programming and the efficiency of IABC. Note that a food source represents a solution to the optimization model in IABC. According to the characteristics of our timetable optimization problem, we encode each solution as a vector of the decision variables, that is, $S = (H, X)$ , where $S$ represents a solution; $H = (h_{1}, h_{2}, \dots, h_{i}, \dots, h_{I - 1})$ ; and $X = (X_{1}, X_{2}, \dots, X_{i}, \dots, X_{I})$ , where $X_{i} = (x_{i, 1}, x_{i, 2}, \dots, x_{i, 2 N - 2}), \forall i \in {1, 2, \dots, I}$ .

The main procedure of IABC is shown in Algorithm 1. $n_{e}$ , $n_{o}$ , and $n_{s}$ are the number of employed bees, onlooker bees, and scout bees, respectively. $n_{b}$ is the total number of all kinds of bees, where $n_{b} = n_{e} + n_{o} + n_{s}$ . $S^{*}$ is the current timetable of a subway system. $S_{n_{b} + 1}$ is used to record the best solution found in each iteration. M₁ is the maximum count of the restart mechanism. It means that IABC would not terminate until it has restarted at least M₁ times (i.e. totally executed for M₁ + 1 times). The restart mechanism ensures that all the other solutions are initialized except the one to record the best result found. M₂ is the maximum iteration count in each execution. The solutions are updated by the employed bees, onlooker bees, and scout bees in each iteration, with the same process in the prior ABC algorithms.^17,32,33

Algorithm 1. Procedure of IABC
Input: $n_{e}$ , $n_{o}$ , $n_{s}$ , $n_{b}$ , $M_{1}$ , $M_{2}$ , $S^{*}$
Output: $S_{1}, S_{2}, \dots, S_{n_{b} + 1}$
Begin
1 $S_{n_{b} + 1} = S^{*}$ ; %init the best solution found
2 for $i = 1$ to $M_{1} + 1$ do
3 for $k = 1$ to $n_{b}$ do
4 randomly generate solution $S_{k}$ via Algorithm 2;
5 end for
6 for $j = 1$ to $M_{2}$ do
7 for $k = 1$ to $n_{b} + 1$ do
8 V_k = Fitness(S_k); %evaluate the fitness value
9 end for
10 $[\tilde{V}, \tilde{d}] = sort ([V_{1}, \dots, V_{n_{b}}], ″ descend ″)$
11 $[B_{1}, \dots, B_{n_{e}}] = [S_{{\tilde{d}}_{(1)}}, \dots, S_{{\tilde{d}}_{(n_{e})}}]$ ; $% n_{e}$ best food
12 $[F_{1}, \dots, F_{n_{e}}] = [{\tilde{V}}_{(1)}, \dots, {\tilde{V}}_{(n_{e})}]$ ; %fitness values
13 $S_{n_{b} + 1} = S_{{\tilde{d}}_{(1)}}$ ; %record the best food
%employed bee phase
14 for $k = 1$ to $n_{e}$ do
15 generate $S_{k}$ from $B_{k}$ via a local search operator;
16 end for
%onlooker bee phase
17 for $k = 1$ to $n_{o}$ do
18 find $B_{l}$ in $[B_{1}, \dots, B_{n_{e}}]$ by roulette wheel selection;
19 generate $S_{n_{e} + k}$ from $B_{l}$ via a local search operator;
20 end for
%scout bee phase
21 for $k = 1$ to $n_{s}$ do
22 generate solution $S_{n_{e} + n_{o} + k}$ via Algorithm 2;
23 end for
24 end for
25 end for
End

Initialization of solutions

At the beginning of IABC, we initialize $n_{b} + 1$ solutions as the basis of optimization. Since our objective is to find a better solution than the currently used timetable on REU, we initialize one solution as the current timetable before the iteration starts. And in the initialization phase after each restart, the best solution found previously is used as one solution. The other $n_{b}$ solutions are generated randomly while ensuring that all of them are feasible, that is, satisfying equations (16)–(21). The procedure of the random generation of a feasible solution is shown in Algorithm 2.

Algorithm 2. Random generation of a feasible solution
Input: $I, N, \underline{h}, \bar{h}, {\underline{x}}_{n}, {\bar{x}}_{n}, C_{1}, C_{2}$
Output: $S$ %a feasible solution
Begin
1 repeat
2 for $i = 1$ to $I - 2$ do
3 $h_{i} = \underline{h} + (\bar{h} - \underline{h}) \times rand$ ; %randomly generate $h_{1}$ to $h_{I - 2}$
4 end for
5 $h_{I - 1} = C_{2} - sum (h_{1}, h_{2}, \dots, h_{I - 1})$ ; %calculate $h_{I - 1}$
6 until $\underline{h} \leq h_{I - 1} \leq \bar{h}$
7 $H = (h_{1}, h_{2}, \dots, h_{I - 1})$ ; %headway time vector
8 for $i = 1$ to Ido
9 Repeat
10 for $n = 1$ to $2 N - 2$ do
11 $x_{i, n} = {\underline{x}}_{n} + ({\bar{x}}_{n} - {\underline{x}}_{n}) \times rand$ ; %randomly generate $x_{i, n}$
12 end for
13 until $\underline{τ} \leq sum (x_{i, 1}, x_{i, 2}, \dots, x_{i, 2 N - 2}) + C_{1} \leq \bar{τ}$
14 $X_{i} = (x_{i, 1}, x_{i, 2}, \dots, x_{i, 2 N - 2})$ ; %dwell time of train i
15 end for
16 $X = (X_{1}, X_{2}, \dots, X_{i}, \dots, X_{I})$ ; %dwell time for all trains
17 $S = (H, X)$ ; %a feasible solution
End

Optimization with employed bee phase

In each iteration, all the solutions are evaluated and a list that contains the n_e best solutions is generated. Each employed bee is assigned to a specific solution in the list, that is, the kth employed bee is assigned to the kth solution in the list, where $k \in {1, 2, \dots, n_{e}}$ . Then, a local search operator is applied to each employed bee to find a new neighbor solution. To improve the local search ability of IABC, a local search operator is randomly chosen from the swap, insertion, mutation, and crossover operators with the probability of $p_{s}$ , $p_{i}$ , $p_{m}$ , and $p_{c}$ , respectively. Note that each solution is encoded as a vector consisted by H and each $X_{i}$ . For the convenience of description, we call H and each $X_{i}$ a basic unit of a solution in the following of this article. To ensure the newly generated solution feasible, the chosen operator is applied on each basic unit of a solution separately. The principles of the swap insertion, mutation, and crossover operators are shown in Figures 4 –7, respectively. C is a basic unit of a solution and $\hat{C}$ is the newly generated basic unit after applying the local search operator. Each basic unit consists of several elements, labeled as $c_{1}, c_{2}, \dots, c_{e}$ , respectively, where $c_{e}$ is the last element of the basic unit, that is, e = I − 1 when the basic unit is H , and e = 2 N − 2 when the basic unit is $X_{i}$ , respectively.

Figure 4.

Principle of swap operator in IABC.

Figure 5.

Principle of insertion operator in IABC.

Figure 6.

Principle of mutation operator in IABC.

Figure 7.

Principle of crossover operator in IABC.

The principle of swap operator is illustrated in Figure 4. Two positions $i, k \in {1, 2, \dots, e}$ are randomly generated and the elements at these two positions are swapped. Note that when the basic unit is $X_{i}$ , to ensure the newly generated solution feasible, the reserve time at these two positions are swapped, rather than the dwell time, that is, $c_{i} - {\underline{c}}_{i}$ is swapped with $c_{k} - {\underline{c}}_{k}$ .

The principle of insertion operator is shown in Figure 5. Two positions $i, k \in {1, 2, \dots, e}$ are randomly generated first, then the element at i is inserted into position k, and the element at k and the elements between i and k are moved by one position toward position i accordingly. Again, reserve time is moved when the basic unit is $X_{i}$ , to make the newly generated solution feasible.

The principle of mutation operator is shown in Figure 6. A position $i \in {1, 2, \dots, e}$ is randomly generated first; then, the element at that position is replaced by a new randomly generated feasible element. Note that a feasible element represents that it satisfies its constraints of the lower and upper limits, that is, equation (16) for headway time and equation (17) for dwell time. Then, repair is applied on the newly generated solution $\hat{C}$ . We calculate the sum of $\hat{C}$ first, that is, $s = \sum \hat{C}$ , then check whether it satisfies equation (18) or (19), depends on the basic unit is H or $X_{i}$ . If the new solution is infeasible, a position is randomly generated, and the element at that position is increased or decreased by one according to the difference of s and its lower/upper limits while ensuring that the newly generated element satisfying its bound limits. This repair process repeats until the newly generated solution is feasible.

The principle of crossover operator is shown in Figure 7. The crossover operator is applied on two basic units with the same meaning in two solutions, for example, two headway time vectors $H_{1}$ and $H_{2}$ in two solutions. During the operation, two positions $i, k \in {1, 2, \dots, e}, i < k$ are randomly generated first. Then, the elements between i and k are exchanged in the two basic units $C_{1}$ and $C_{2}$ . Thereafter, repair is applied on the newly generated ${\hat{C}}_{1}$ and ${\hat{C}}_{2}$ , respectively. In the repair process, we calculate the sum of the basic unit, that is, $s_{j} = \sum C_{j}, j \in {1, 2}$ , then check whether $s_{j}$ is feasible. If not, a position is randomly generated in ${\hat{C}}_{j}$ , and the element at that position is increased/decreased by one according to the difference of $s_{j}$ and its lower/upper limits, while ensuring that the newly generated element is feasible. The repair process repeats until the newly generated solution is feasible.

Optimization with onlooker bee phase

The process of onlooker bee phase is similar to that of the employed bee phase, except that each onlooker bee is assigned to a solution in the n_e best solutions list by spinning a roulette wheel. Suppose the n_e best solutions in the list have been ordered based on their fitness values. The ordered fitness values are denoted as $F_{1}, \dots, F_{n_{e}}$ , respectively, as shown in Algorithm 1. In the roulette wheel selection process, the probability of a solution been chosen is proportional to its fitness value. The detailed procedure for an onlooker bee to choose a solution in the list is shown as follows:

Let $p_{0} = 0$ ;

Calculate $p_{i} = \sum_{j = 1}^{i} F_{j}, i = 1, 2, \dots, n_{e}$ ;

Randomly generate a number $r \in (0, p_{n_{e}})$ ;

Compare r with each $p_{i}$ and select the ith solution in the list satisfying $r \in (p_{i - 1}, p_{i})$ .

For each onlooker bee, once a food source is chosen, a local search operator is chosen and it is applied on this solution to generate a new one. The processes to choose an operator and to apply it are the same with that of the employed bee phase mentioned above.

Optimization with scout bee phase

To avoid IABC from being trapped in a local optimum, each scout bee randomly generates a new solution in each iteration. The random generation of a new solution is the same with that introduced in the initialization of solutions, which is shown in Algorithm 2.

Experimental results and analysis

To illustrate the correctness of the mathematical model and the effectiveness of the proposed IABC, several numerical experiments are conducted based on the data obtained from Yanfang Line, a subway line in Beijing, China. The running time at each section is given in Table 1, the sections in each substation are given in Table 2, and the other parameters of Yanfang Line are given in Table 3.

Table 1.

Actual running time at each section.

Section index	1	2	3	4	5	6	7	8	9	10	11	12	13	14
Δt_n (s)	121	93	111	227	79	118	133	132	118	79	227	111	94	115

Table 2.

Sections in each substation.

Substation	Section index				Substation	Section index
S ₁	1	2	13	14	S ₂	3	4	11	12
S ₃	5	6	9	10	S ₄	7	8	−	−

Table 3.

Parameters for the Beijing Yanfang Line.

Parameter	Value	Unit	Parameter	Value	Unit
N	8	–	h	482	s
$a_{n}^{1}$	0.8	m/s²	$\underline{h}$	422	s
$a_{n}^{3}$	1	m/s²	$\bar{h}$	542	s
$Δ t_{n}^{1}$	27	s	$Δ \underline{T}$	2348	s
$Δ t_{n}^{3}$	21	s	$Δ \bar{T}$	2468	s
$Δ t_{t}$	230	s	$x_{n}$	30	s
$η_{1}$	0.7	–	${\underline{x}}_{n}$	25	s
$η_{2}$	0.8	–	${\bar{x}}_{n}$	35	s
$β$	0.05	–	m	287,080	kg

The parameters of IABC in the experiments are shown in Table 4. To compare IABC and GA,^12,13,16 we keep the number of solutions searched nearly equal for the two algorithms. Thus, the population size of GA is set to be 50 and its maximum iteration count is 200. The local search operators used in GA include selection, crossover, and mutation. Their detailed processes can be found in the work by Yang et al.¹⁶

Table 4.

Parameters of IABC.

$n_{e}$	$n_{o}$	$n_{s}$	$M_{1}$	$M_{2}$	$p_{s}$	$p_{i}$	$p_{m}$	$p_{c}$
10	10	20	4	50	0.1	0.1	0.2	0.6

IABC: improved artificial bee colony.

The timetable optimization problem is solved by IABC and GA, respectively. We compare their results with different train trip numbers (i.e. different values for I). Our previous work has proved that REU can be improved by different dwell time for different trains, but the increment is little, while it increases the solution space dramatically.³⁴ Thus, we keep identical dwell time for different trains in the experiments in this work, that is, $x_{i, n} = x_{j, n}, \forall i, j \in {1, 2, \dots, I}$ .

Performance indicators to evaluate the effectiveness of IABC

In order to evaluate the effectiveness of IABC in solving our timetable optimization problem, we apply it and GA, respectively, and obtain the optimized REU with different train trip numbers. The following performance indicators are used to analyze the experimental results and evaluate the effectiveness of IABC in the following sections:

REU. The values of REU for the optimized timetable and currently used timetable are shown, respectively. They can be compared directly. A greater REU is better than a smaller one, as it represents more energy can be saved.

Improvement ratio of optimized REU over the current timetable. It shows the relative improvement of the optimized REU over the current timetable. Suppose the optimized REU obtained by IABC is denoted as REU^I and REU for the current timetable is denoted as REU^C, then the improvement ratio of optimized REU by IABC over the current timetable is calculated as $I R^{I} = R E U^{I} / R E U^{c} \times 100 %$ . The improvement ratio of optimized REU by GA is obtained in the same way.

Money saving with the optimized timetable. To show the effect of the optimized timetable in energy saving intuitively, we calculate the amount of energy saved by the optimized timetable and obtain the corresponding money saving by considering the price of electricity in US$. It can be used as a reference to judge the economic effect of the timetable optimization model.

Experimental results analysis

The experimental results are shown in Table 5.

Table 5.

Optimized REU by IABC and GA.

I	131	100	50	40	30	20	10
REU^c (kJ/kg)	30.55	23.16	11.23	8.84	6.46	4.07	1.68
REU^I (kJ/kg)	54.44	38.63	21.60	19.34	15.66	10.44	4.61
IR^I (%)	78	67	92	119	143	156	174
REU ^G (kJ/kg)	38.71	31.61	16.48	14.48	11.35	7.95	4.17
IR^G (%)	27	36	47	64	76	95	147

GA: genetic algorithm; IABC: improved artificial bee colony; REU: regenerative energy utilization.

I: train trip number; REU^c: REU of the currently used timetable; REU^I: REU of the timetable optimized by IABC; REU^G: REU of the timetable optimized by GA; IR^I: REU improvement ratio of the timetable optimized by IABC over the current one, that is, $I R^{I} = R E U^{I} / R E U^{c} \times 100 %$ ; IR^G:REU improvement ratio of the timetable optimized by GA over the current one, that is, $I R^{G} = R E U^{G} / R E U^{c} \times 100 %$ .

It is clear that REU improves dramatically after optimization, and IABC performs better than GA in solving our timetable optimization problem. IR^G decreases with I; IR^I first decreases with I and increases a little when I is bigger than 100. The actual number of train trips is 131 in this subway line. It is seen that REU^I is 78% more than REU^c when I is 131, which means a great amount of energy can be saved if we optimized the timetable by IABC.

To compare the effectiveness of IABC and GA intuitively, Figure 8 shows the optimization process with I = 10 and I = 50 as two examples. Note that the conclusion is the same for all train trip numbers. The horizontal axis is the iteration count, and the vertical axis is the REU improvement ratio over the current timetable, that is, IR^I and IR^G. We can see that IABC obtains better results than GA, almost at any iteration count.

Figure 8.

Optimization process of IABC versus GA.

The near-optimal solution obtained by applying IABC when the number of train trips is 131 is shown in Figure 9 as an example. There are two subfigures in Figure 9. To see the changes of the optimized timetable with the current one clearly, the upper part of Figure 9 shows the differences between the optimized headway time and the headway time in the current timetable, and the lower part shows the optimized dwell time and that in the current timetable, respectively.

Figure 9.

Optimal solution obtained by IABC.

From Figure 9 and Table 5, it is clear that by slightly adjusting the headway time and dwell time in a subway system, REU can be dramatically improved. This represents tremendous cost saving for the subway operators, and significant contribution to the environment, at the expense of very little impact on the service quality of the subway system.

Total money saving from the electrical energy by the optimized timetable for a subway system is determined as equation (22)

ϒ = \frac{(RE U^{I} - RE U^{c}) \times 1000 \times m \times I \times 365 \times μ}{(3600 \times 1000)}

(22)

where $ϒ$ is the money saved (US$); µ is the unit price of the electricity (US$/kWh). Take I = 131 and µ = 0.112 for an example, it is easy to obtain that the total money saved in this subway in a year is (38.63 − 31.61) × 1000 × 287,080 × 131 × 365 × 0.112/(3600 × 1000) = 2,997,917.88 (US$).

From above, we can see that this work is very useful for the decision makers to improve the timetable in a subway system on energy saving. Based on the proposed mathematical model and IABC, we can easily obtain an optimized timetable and thus save a great amount of money by maximizing REU in a subway system.

Conclusion

In this work, we propose a timetable optimization problem to maximize REU in a subway system, and its mathematical model is formulated with the decision variables of headway time and dwell time. RE generated from braking trains can be utilized by all traction trains in the same substation in the model. Some required constrains are considered, for example, the constant operation time constraint with fixed number of train trips in a day. Then, IABC is designed to solve our problem. A restart mechanism is adopted to reduce its chance of being trapped in a local optimum. The procedure that generates a new solution is specially designed to make every new solution feasible, including the random generation of a solution and the repair of a solution after a local search operator has been applied. Numerical experiments are conducted based on the data from a subway line in China. Experiment results prove the correctness of the model and the effectiveness of IABC in solving timetable optimization problems. IABC is also compared with and performs better than GA in the experiments.

Footnotes

Handling Editor: Jiangchen Li

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article. This paper was partially supported by the National Key Research and Development Program of China under grant 2018YFB1201501, Beijing municipal natural science foundation under grant L161008, Fundamental Research Funds for the Central Universities under grant 2016JBZ004, TCT Funding Program under grant 9907006510, and Chinese Railway Certification Center funding program under grant 1852ZJ1303, Beijing Laboratory of Urban Rail Transit, and China Scholarship Council.

ORCID iD

Hongjie Liu

References

Tian

Weston

Zhao

et al . System energy optimisation strategies for metros with regeneration. Transport Res C: Emer 2017; 75: 120–135.

González-Gil

Palacin

Batty

et al . A systems approach to reduce urban rail energy consumption. Energ Convers Manage 2014; 80: 509–524.

González-Gil

Palacin

Batty

Optimal energy management of urban rail systems: key performance indicators. Energ Convers Manage 2015; 90: 282–291.

Yang

Ning

et al . A survey on energy-efficient train operation for urban rail transit. IEEE T Intell Transp 2016; 17: 2–13.

Kampeerawat

Koseki

A strategy for utilization of regenerative energy in urban railway system by application of smart train scheduling and wayside energy storage system. Energy Procedia 2017; 138: 795–800.

González-Gil

Palacin

Batty

Sustainable urban rail systems: strategies and technologies for optimal management of regenerative braking energy. Energ Convers Manage 2013; 75: 374–388.

Ramos

Pena

Fernandez

et al . Mathematical programming approach to underground timetabling problem for maximizing time synchronization. In: Proceedings of the international conference on industrial engineering and industrial management, Madrid, 5–7 September 2007, pp.88–95. IEEE.

Nasri

Moghadam

Mokhtari

. Timetable optimization for maximum usage of regenerative energy of braking in electrical railway systems. In: Proceedings of the international symposium on power electronics electrical drives automation and motion, Pisa, Italy, 14–16 June 2010, pp.1218–1221. New York: IEEE.

Kim

Han

. A model and approaches for synchronized energy saving in timetabling. In: Proceedings of the 9th world congress on railway research, Lille, France, 20–26 May 2011. IEEE.

10.

Peña-Alcaraz

Fernández

Cucala

et al . Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy. Proc IMechE, Part F: J Rail and Rapid Transit 2012; 226: 397–408.

11.

Fournier

Mulard

Fages

. Energy optimization of metro timetables: a hybrid approach. In: Proceedings of the 18th international conference on principles and practice of constraint programming, Quebec, Canada, 8–12 October 2012, pp.7–12. Springer.

12.

Yang

Gao

et al . A cooperative scheduling model for timetable optimization in subway systems. IEEE T Intell Transp 2013; 14: 438–447.

13.

Yang

Ning

et al . A two-objective timetable optimization model in subway systems. IEEE T Intell Transp 2014; 15: 1913–1921.

14.

Nag

Pal

. Optimal Design of Timetables to Maximize Schedule Reliability & Minimize Energy Consumption, Rolling Stock and Crew Deployment (September 10, 2003). Proceedings of 2nd UIC(International Congress of Railways) Energy Efficiency Conference, February 2–4, 2004, Paris. Available at SSRN: https://ssrn.com/abstract=2323254

15.

HK.

An energy-efficient scheduling and speed control approach for metro rail operations. Transport Res B: Meth 2014; 64: 73–89.

16.

Yang

Chen

et al . An energy-efficient scheduling approach to improve the utilization of regenerative energy for metro systems. Transport Res C: Emer 2015; 57: 13–29.

17.

Karaboga

An idea based on honey bee swarm for numerical optimization. Technical report-TR06, October 2005. Kaiseri, Turkey: Department of Computer Engineering, Faculty of Engineering, Erciyes University.

18.

Deng

A discrete artificial bee colony algorithm for minimizing the total flow time in the blocking flow shop scheduling. Chin J Chem Eng 2012; 20: 1067–1073.

19.

Pan

Gao

Pareto-based discrete artificial bee colony algorithm for multi-objective flexible job shop scheduling problems. Int J Adv Manuf Technol 2011; 55: 1159–1169.

20.

Artificial bee colony algorithm-based parameter estimation of fractional-order chaotic system with time delay. IEEE/CAA J Automat Sinica 2017; 4: 107–113.

21.

Singh

An artificial bee colony algorithm for the leaf-constrained minimum spanning tree problem. Appl Soft Comput 2009; 9: 625–631.

22.

Q-K

Pan

Tasgetiren

Suganthan

et al . A discrete artificial bee colony algorithm for the lot-streaming flow shop scheduling problem. Inform Sciences 2011; 181: 2455–2468.

23.

Tian

Ren

Zhou

Dual-objective scheduling of rescue vehicles to distinguish forest fires via differential evolution and particle swarm optimization combined algorithm. IEEE T Intell Transp 2016; 17: 3009–3021.

24.

Wang

Zheng

A hybrid estimation of distribution algorithm for unrelated parallel machine scheduling with sequence-dependent setup times. IEEE/CAA J Automat Sinica 2016; 3: 235–246.

25.

Mareda

Gaudard

Romerio

A parametric genetic algorithm approach to assess complementary options of large scale windsolar coupling. IEEE/CAA J Automat Sinica 2017; 4: 260–272.

26.

Hou

Zhou

et al . Pareto optimization for scheduling of crude oil operations in refinery via genetic algorithm. IEEE T Syst Man Cy: S 2017; 47: 517–530.

27.

Dong

Zhou

A supervised learning and control method to improve particle swarm optimization algorithms. IEEE T Syst Man Cy: S 2017; 47: 1135–1148.

28.

Tang

Luo

Yang

et al . A chance constrained optimal reserve scheduling approach for economic dispatch considering wind penetration. IEEE/CAA J Automat Sinica 2017; 4: 186–194.

29.

Kang

Feng

Zhou

et al . Optimal load scheduling of plug-in hybrid electric vehicles via weight-aggregation multi-objective evolutionary algorithms. IEEE T Intell Transp 2017; 18: 2557–2568.

30.

Guo

Liu

Zhou

et al . Disassembly sequence optimization for large-scale products with multiresource constraints using scatter search and Petri nets. IEEE T Cybernetics 2016; 46: 2435–2446.

31.

Guo

Liu

Zhou

et al . Dual-objective program and scatter search for the optimization of disassembly sequences subject to multiresource constraints. IEEE T Automat Sci Eng 2018; 15: 1091–1103.

32.

Karaboga

Basturk

On the performance of artificial bee colony (ABC) algorithm. Appl Soft Comput 2008; 8: 687–697.

33.

Karaboga

Akay

A comparative study of artificial bee colony algorithm. Appl Math Comput 2009; 214: 108–132.

34.

Liu

Zhou

Guo

et al . An ABC-based subway timetable optimization model for regenerative energy utilization. In: Proceedings of the international conference on systems, man, and cybernetics Miyazaki, Japan, 7–10 October 2018. IEEE.