A Dual-Objective Substation Energy Consumption Optimization Problem in Subway Systems

: Maximizing regenerative energy utilization is an important way to reduce substation energy consumption in subway systems. Timetable optimization and energy storage systems are two main ways to improve improve regenerative energy utilization, but they were studied separately in the past. To further improve energy conservation while maintaining a low cost, this paper presents a strategy to improve regenerative energy utilization by an integration of them, which determines the capacity of each Wayside Energy Storage System (WESS) and correspondingly optimizes the timetable at the same time. We ﬁrst propose a dual-objective optimization problem to simultaneously minimize substation energy consumption and the total cost of WESS. Then, a mathematical model is formulated with the decision variables as the conﬁguration of WESS and timetable. Afterwards, we design an (cid:101) -constraint method to transform the dual-objective optimization problem into several single-objective optimization problems, and accordingly design an improved artiﬁcial bee colony algorithm to solve them sequentially. Finally, numerical examples based on the actual data from a subway system in China are conducted to show the effectiveness of the proposed method. Experimental results indicate that substation energy consumption is effectively reduced by using WESS together with a correspondingly optimized timetable. Note that substation energy consumption becomes lower when the total size of WESS is larger, and timetable optimization further reduces it. A set of Pareto optimal solutions is obtained for the experimental subway line—based on which, decision makers can make a sensible trade-off between energy conservation and WESS investment accordingly to their preferences.


Introduction
Energy conservation in subway systems has attracted great attention in recent years. As a great proportion of total energy is consumed by train traction systems [1][2][3][4], many measures have been taken to reduce traction energy consumption. The adoption of energy-efficient train operations has been a focus in the past decades [5]. It aims to find a driving strategy which consumes the least energy. However, Regenerative Energy (RE) is not considered in these studies, which limits the effect of this method on energy conservation [3,6]. Regenerative energy is the electrical energy converted from kinetic energy by a regenerative braking system during the braking of trains. As a potential energy supply, regenerative energy can be used to accelerate trains. Regenerative energy takes a great amount of the total energy 1. An integration of timetable optimization and WESS is proposed to maximize regenerative energy utilization, thus to minimize substation energy consumption in a subway system. 2. To maximize energy saving with the least cost, a dual-objective optimization model is thereby formulated to simultaneously minimize substation energy consumption and the total investment cost of WESSs. Note that a subway line is divided into several electricity supply intervals. One WESS is installed in each electricity supply interval. The size of each WESS can be different from each other, and both of them are greater than or equal to 0. Their total size determines the total cost of WESSs. 3. To solve the proposed dual-objective optimization problem, an -constraint method is first designed to transform it into several single-objective optimization problems. Then, an improved artificial bee colony algorithm is designed to solve these single-objective optimization problems sequentially. 4. Numerical examples are constructed based on the actual data obtained from a subway system in China to show the effectiveness of the proposed resolution methods. A set of Pareto optimal solutions is obtained. i.e., for each value of the total size of WESSs, the minimal substation energy consumption is determined, and the optimal configuration of each WESS and the correspondingly optimized timetable are obtained to reach the maximum energy saving.
The remainder of this paper is organized as follows. We review the related work on improvement of regenerative energy utilization in Section 2. In Section 3, a dual-objective optimization problem is proposed to simultaneously minimize substation energy consumption and the corresponding WESS costs, and a mathematical model for the proposed problem is formulated. To solve the dual-objective optimization problem, we design both an -constraint method and an improved artificial bee colony algorithm in Section 4. Then, numerical experiments are conducted in Section 5 to show the effectiveness of the proposed method. Finally, Section 6 concludes this paper and points out some future research directions.

Literature Review
In this section, we review the studies on improvement of regenerative energy utilization. As two main measures, timetable optimization and the application of energy storage systems are studied separately in the literature. However, their integration was not considered until recently. Consequently, a dual-objective optimization problem to minimize energy consumption and investment cost was also seldom studied. The main related publications are reviewed in the following.
An optimized timetable can improve regenerative energy utilization between traction and braking trains, hence reduce substation energy consumption in a subway system. In addition, the cost of timetable optimization is relatively low. Therefore, timetable optimization has been studied by many researchers to save energy [7,8,[13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]. However, ESS was seldom considered in these studies, which limits their effects on energy saving. Furthermore, most of the studies are single-objective optimization problems and only a few of them are dual-objective optimization ones, which limits the resolution methods to be applied in our work.
In a single-objective timetable optimization problem to save energy, the direct objective is to maximize regenerative energy utilization, which may be represented by energy or the overlap time between trains' traction and braking phases. Mathematical models of timetable and regenerative energy utilization are formulated in these studies, and the problems are usually solved by an analytical method or an intelligent algorithm, which are helpful for further research. For example, Ramos et al. [13] present a timetabling problem to maximize the overlap time of the speed-up and slow-down actions of trains. Kim et al. [15] propose a multi-criteria mixed integer program to minimize the peak energy consumed and to maximize regenerative energy utilization through timetable optimization. Pena et al. [16] propose a timetable optimization model for an underground rail system to maximize regenerative energy utilization. Fournier et al. [17] develop an optimization model to maximize regenerative energy utilization by subtly modifying dwell time for trains at stations, and a hybrid genetic/linear programming algorithm is implemented to tackle this problem. Yang et al. [18] propose a cooperative scheduling model to maximize the overlap time of accelerating and braking processes of adjacent trains. Yang et al. [21] formulate an integer programming model with real-world speed profiles to minimize traction energy consumption by adjusting dwell time. They coordinate the arrivals and departures of trains in the same electricity supply interval (ESI), such that regenerative energy is effectively utilized. A genetic algorithm (GA) is designed to solve their problems in both [18,21]. Zhao et al. [20] develop a nonlinear integer program to maximize regenerative energy utilization, which searches for the optimal headway and dwell time at each station. Gong et al. [23] present a timetable optimization model to maximize regenerative energy utilization with dwell time control, and GA is used to find a near-optimal solution.
Few studies on timetable optimization are dual-objective optimization problems, where energy conservation and passenger time are usually minimized at the same time. The way to formulate a dual-objective optimization model and the possible resolving method are referrable. e.g., Yang et al. [8] propose a dual-objective timetable optimization model to coordinate up and down trains at the same station to improve regenerative energy utilization and reduce passenger waiting time. Two objectives are combined into one by a weighted-sum method, and it is solved by GA. Zhao et al. [19] propose a dual-objective optimization problem to maximize regenerative energy utilization measured by the overlap time and to shorten total passenger time. Two objectives are combined into one through weighting, and a simulated annealing (SA) method is designed to solve it. Xu et al. [24] propose a dual-objective optimization problem to minimize both passenger time and traction energy, by controlling running time at each section and dwell time at each platform. They adopt a weighted-sum method to combine two objectives into one, and designed a genetic algorithm to solve it.
Installing ESS also can improve regenerative energy utilization, thereby reducing substation energy consumption in subway systems. Therefore, few researchers have studied the application of ESS in subway systems. Although timetable optimization is seldom considered in these studies, the modelling method of ESS is referable to our work. Ceraolo et al. [11] develop models to analyze the impact of regenerative braking in high-speed railway systems, where ESS is used. Feasibility of using wayside and on-board ESS is analyzed, respectively. They evaluate the cost-effectiveness of different solutions by taking into account the capital cost of the investment and annual energy saving. The results prove the effectiveness on improving regenerative energy utilization through ESS application. Ciccarelli et al. [12] propose a control strategy for on-board super-capacitors integrated with motor drive control. Simulation results show its effectiveness on energy saving and reducing the voltage surge at the overhead contact line during train braking. Liu et al. [28] propose a single-objective timetable optimization problem with application of WESS to minimize total energy consumption. An algorithm integrating tabu search and simulation is designed to solve it. Experimental results prove the effectiveness of WESS on energy saving. The model of WESS energy management is relatively simple in their work. Gao et al. [29] propose a control strategy to use super-capacitors as WESS, which is helpful for WESS energy management. Huang et al. [30] propose an energy-saving model to optimize trains' speed profiles in a subway system. On-board ESS was used as a basis of their optimization problem and the running time is allowed to be optimized within a predefined window. A dynamic programming method is used to solve their problems. Kampeerawat and Koseki [10] propose a single-objective optimization problem to reduce substation energy consumption, which raised the studies on integration of timetable optimization and WESS. A mathematical model is formulated to minimize a linear weighted sum of substation energy consumption and the energy capacity of WESS; then, GA is designed to solve their problem. Although different weight factors can be assigned in theory, it is hard to find all the non-dominated solutions for the original dual-objective optimization problem. Ahmadi et al. [31] propose to reduce substation energy consumption in subway systems by simultaneous application of WESS and speed profile optimization. To demonstrate the validity of the proposed method, they first optimized the configuration of WESS together with the real world (i.e., non-optimized) speed profiles to reduce substation energy consumption. Then, the speed profiles and the configuration of WESS were simultaneously optimized. GA was used to solve their problems. Experimental results prove that substantial reduction in substation energy consumption was achieved and total size of WESS is decreased when speed profiles were optimized in comparison with the non-optimized speed profiles.
From the above we can see that the integration of timetable optimization and WESS is very few in the existing studies, which highlights the first feature of our work. i.e., the energy saving strategy proposed in this work, which combines timetable optimization and WESS, is challenging and different from the existing studies. Consequent on the integration of the two different methods, to simultaneously minimize substation energy consumption and the financial cost of WESS, a dual-objective optimization problem is proposed, which is also innovative. Furthermore, to solve a dual-objective optimization problem, a weighted-sum method is often adopted in the previous studies to transform it into a single-objective optimization problem, where the optimal solutions are limited by the weights adopted. While in this work we design an -constraint method to transform the original dual-objective optimization problem into several single-objective ones, then the Pareto front is able to be obtained in theory [32,33]. Finally, instead of applying the most frequently used genetic algorithm to solve the single-objective optimization problems, an improved artificial bee colony algorithm is designed and experimental results prove its better performance over a genetic algorithm.

Problem Assumptions
According to the characteristics of the dual-objective substation energy-consumption optimization problem, we formulate a mathematical model based on the following assumptions.
1. There is reserve time for dwell time at each platform. 2. The subway line is divided into several electricity supply intervals. Each electricity supply interval provides energy for several stations and sections. Electricity can be transmitted within an an electricity supply interval. 3. Each electricity supply interval includes one substation, one WESS and a bank of resistors. Each WESS can be installed near the place of a substation. Each WESS consists of several parallel Basic Energy Storage Modules (BESMs). The total cost of all WESSs grows linearly with their total energy capacity (represented by the number of BESMs in this work).
4. The transmission loss coefficient of regenerative energy is a constant. 5. The other electrical utilities (e.g., lights and air conditioners) are not considered during the simulation, i.e., electrical energy can only be consumed by traction trains, resistors and charged into WESS. 6. Regenerative energy is fed back into the power supply line and can be immediately used by traction trains in the same electricity supply interval. If regenerative energy cannot be used fully by traction trains, the surplus is charged into WESS and consumed by resistors if needed. 7. The condition of a charge/discharge state transition for a WESS is defined as follows: • When trains are braking and the power of surplus regenerative energy exceeds a predefined limit, WESS in the same electricity supply interval begins to charge; When the power of surplus regenerative energy gets smaller, the corresponding WESS stops charging.

•
When trains are in a traction phase and the power of surplus traction energy exceeds a predefined limit, WESS begins to discharge; when the surplus traction energy gets smaller, WESS stops discharging (the surplus traction energy demand is satisfied by a substation).

Train Movement and Timetable Modeling
As shown in Figure 2, we label each station, platform and section with indices y ∈ {1, 2, · · · , Y}, n ∈ {1, 2, · · · , N} and l ∈ {1, 2, · · · , L} in this work, respectively. Note that total numbers of stations, platforms and sections are Y, N and L, respectively. There are two platforms at each station, except the terminal station Y, since we dismiss the detailed turnaround process in this work. Thus, we have N = 2Y − 1. The platforms in station y ∈ {1, 2, · · · , Y − 1} are labeled as platforms n and N − n + 1 respectively, by noting y = n. In addition, the platform in terminal station Y is labeled as platform N + 1 2 . A section is used to connect any two successive platforms, thus L = N − 1 = 2Y − 2. The section connecting platforms n and n + 1 is labeled as section l, by noting l = n. In a subway system, every train travels at the same closed path in a subway system with a time interval. Every train i ∈ {1, 2, · · · , I} begins its travel from platform 1. It runs along the down direction, through each section and platform sequentially, until it arrives at the terminal station. Then, it turns around to the up direction and repeats a similar process until it arrives the last platform N and finishes its travel. In the process, it stays at each platform n ∈ {1, 2, · · · , N − 1} for some time period x i,n for passengers to get on and off the train.  A timetable is often used to describe the key time points of trains in a subway system. Given a timetable, we can easily obtain all the key time points of trains, including the time instant when each train arrives at and departs from any platform. Dwell time and headway time are two of the most important factors in a timetable, and they are to be optimized in this work. The running process at a section of each train is also important for calculating energy consumption. It is often divided into three phases, i.e., traction, coasting and braking phases, in the timetable optimization studies [18,21]. Detailed process at any section l and the related parameters are illustrated in Figure 3. A red dotted line stands for a traction phase, a green solid line for a braking phase and a blue dashed line for a coasting phase. t u i,l is the time instant when train i is at time point u and section l, where u ∈ {1, 2, 3, 4} is the key time point index of the running phases at a section, and its value from 1 to 4 represents the start time of a traction phase (train leaves a platform), switch time from traction to a coasting phase, switch time from coasting to a braking phase and the end time of a braking phase (train arrives at a platform), respectively. v u i,l is the train speed when train i is at time point u and section l. r γ l is the time duration of phase γ at section l, where γ ∈ {1, 2, 3} is the train running phase index at a section, and its value from 1 to 3 represents a traction, coasting and braking phase, respectively. a γ l is the train acceleration in phase γ at section l. x i,n is the dwell time for train i ∈ {1, 2, · · · , I} at platform n ∈ {1, 2, · · · , N − 1}, which is the time interval from train i arrives at to it leaves from platform n. Based on the above, the time of all the key points for all trains is determined by the time instant when the first train starts its travel (t 4 1,0 ) via the following procedures: 1. Determine the time of the other key points for train i at section l by that train's start time at that section: 2. Determine t 1 i,l by the time instant when train i starts its travel: 3. Determine t 4 i,0 by using t 4 1,0 : where t 4 1,0 is a given constant which represents the time instant that the first train ends its braking phase at a virtual section 0 (also means arrival time at platform 1); h i is the headway time between trains i and i + 1 and it is represented by the time difference of their arrival time at platform 1 in this work, i.e., Travel time of train i is obtained as: where τ i is the travel time duration of train i, i.e., Operation time of a subway system is obtained as: where Λ is the operation time duration of a subway system. It is represented by the first platform in this work. Train speed is 0 when a train stops at a platform, and it obeys uniformly accelerated/decelerated motion in every phase at a section. Thus, train speed is determined as follows: where i ∈ {1, 2, · · · , I} and l ∈ {1, 2, · · · , L}; v i (t) is the train speed at any time t. Note that a 1 l and a 3 l are the maximum traction and braking acceleration of a train at section l respectively, which are both given constants. In addition, it is easy to obtain a 2 l as a 2 l = a 1 l r 1 l − a 3 l r 3 l r 2 l . Thus, the kinetic energy of each train at any time is determined actually.

Energy Consumption Calculation
Each electricity supply interval z ∈ {1, 2, · · · , Z} has several energy suppliers and consumers, where the energy suppliers include a substation, braking trains and WESS in a discharging state, and energy consumers include traction trains, a charging WESS and resistors installed in this electricity supply interval, as shown in Figure 4. P a z (t) is the power of traction energy demand in electricity supply interval z at time t; P b z (t) is the power of available regenerative energy in electricity supply interval z at time t; P c z (t) is the charging power of WESS z at time t; P d z (t) is the discharging power of WESS z at time t; P r z (t) is the power of energy consumed by the resistors in electricity supply interval z at time t; P s z (t) is the power of energy supplied by substation z at time t.
According to the energy conservation law, the total power of energy suppliers should be the same as that of energy consumers, i.e., Note that traction trains consume electrical energy and convert it into kinetic energy, and braking trains supply regenerative energy to the power supply line by converting from kinetic energy. The electrical energy demand from traction trains are satisfied by regenerative energy first. If regenerative energy is not enough, the surplus traction energy demand are satisfied by the discharging energy of WESS and the energy from a substation. Otherwise, if there is surplus regenerative energy, it can be consumed by resistors into heat and be charged into WESS for later use.

Total Traction Energy Demand
The required electrical energy for accelerating train i at section l (denoted as E a i,l ) is determined by the kinetic energy increased in its traction phase, and the required electrical power for accelerating train i at section l at time t (denoted as P a i,l (t)) is its derivative. They are determined as follows: where ψ 1 and ψ 2 are constants satisfying ψ 1 = m/2η 1 and ψ 2 = m a 1 l 2 /η 1 , respectively, where m is train mass and η 1 is the conversion efficiency from electrical energy to a train's kinetic energy. Thus, the total traction energy demand of all trains (denoted as E a ) in a subway system is determined as The power of traction energy demand in electricity supply interval z at time t is determined by all the traction trains in that electricity supply interval, i.e., Note that the total traction energy demand E a is a constant in this work, as it is not affected by the change of headway time and dwell time.

Available Regenerative Energy from Braking Trains
The available regenerative energy from train i at section l (denoted as E b i,l ) is determined by the kinetic energy decreased in its braking phase, and the power of regenerative energy at time t (denoted as P b i,l (t)) is the derivative of regenerative energy. They are determined as follows: where ψ 3 and ψ 4 are constants satisfying ψ 3 = mη 2 (1 − ζ 1 ) /2 and ψ 4 = m(a 3 l ) 2 η 2 (1 − ζ 1 ), respectively; η 2 is the conversion efficiency from kinetic energy to regenerative energy; ζ 1 is the transmission loss coefficient of regenerative energy.
The total available regenerative energy from all braking trains in a subway system (E b ) is The total power of regenerative energy from braking trains in electricity supply interval z is

Energy Dynamics of Wess
As described in Section 3.1, one WESS is installed in each electricity supply interval, and it is comprised of several BESMs in parallel. Thus, the energy capacity and maximum power of WESS z are determined by the number of BESMs κ z and BESM specification. i.e., where C e z and P e z is the maximum power and energy capacity of WESS z respectively; C M and P M are given constants, represents the maximum power and energy capacity of a BESM respectively; κ z is the number of BESMs in WESS z ∈ {1, 2, · · · , Z}. Thus, the total energy capacity of all WESSs (denoted as C e ) in a subway system is: where K is the total number of BESMs in a subway system, satisfying K = ∑ Z z=1 κ z . As for the control strategy design and modelling of ESS, there are several good references in smart home domains. For example, Carli and Dotoli [34,35] present a distributed technique to control the charging/discharging of ESS, controllable electrical appliances and renewable energy sources in a smart city scenario. Sperstad et al. [36] formulate the charging/discharging of ESS as a multi-period optimal power flow problem, and propose a framework of methods and models to handle the uncertainties in energy exchange due to distributed wind and solar photovoltaic power generation. These studies are very valuable for the control strategy and explicit modelling of ESS. Nevertheless, besides timetable optimization, this paper mainly focuses on the optimization of WESS configuration in a subway system, i.e., to determine the capacity of the WESS in each electrical supply interval. Thus, the control strategies of WESS are reasonably simplified as in [11,29]. For example, we assume that each WESS can only be in a charging, discharging or idle state at a time, the abnormal case that charging and discharging commands are simultaneously sent to a ESS is not considered.
Note that, according to [11,29], the maximum charging and discharging power of WESS are both variables changing with its state-of-charge, as shown in Figure 5. The horizontal axes are state-of-charge, and the vertical axes are maximum charging/discharging power of a WESS. Note that WESS cannot be thoroughly discharged. For example, the discharging process of WESS stops when its state-of-charge decreases to β 0 . The state-of-charge of WESS z at time t is defined as: where s z (t) and Q e z (t) are the state-of-charge and energy stored in WESS z at time t, respectively. The maximum charging and discharging power of WESS z at time t (denoted as P c z (t) and P d z (t) respectively) are determined as follows: where β 2 is the point where the maximum charging power of WESS begins to decrease; β 1 is the point where the maximum discharging power of WESS begins to decrease, and β 0 is the point where the maximum discharging power of WESS decreases to 0. The energy stored in WESS z at time t is determined by its initial state, charging and discharging energy, i.e., where Q e z (t = 0) is the energy stored in WESS z at the start time; η 3 and η 4 are the energy conversion efficiencies when a WESS is in a charging and discharging state, respectively.

Real-time charging power of a WESS
WESS charges when the regenerative energy from braking trains is more than the energy demand from traction trains, and the surplus regenerative energy exceeds the threshold of WESS charging. When WESS charges, its real-time charging power (denoted at P c z (t)) is determined by the surplus regenerative energy and maximum charging power of WESS, i.e., where P + z (t) is the power of surplus regenerative energy in electricity supply interval z , i.e., P + z (t) = max P b z (t) − P a z (t), 0 ; ζ 2 is the proportion coefficient of surplus regenerative energy charging into WESS. Note that the surplus regenerative energy is partially charged into WESS and partially consumed via resistors into heat if needed.

Real-time discharging power of a WESS
WESS discharges when the power demand from traction trains is more than regenerative energy from braking trains, and the surplus traction energy demand is greater than the threshold of WESS discharging. When WESS discharges, its real-time discharging power (denoted at P d z (t)) is determined by the surplus traction energy and maximum discharging power of WESS. i.e., where P − z (t) is the power of surplus traction energy in electricity supply interval z, i.e., P − z (t) = max P a z (t) − P b z (t), 0 ; ζ 3 is the proportion coefficient of electrical energy supplied by WESS when it is discharging. Note that the surplus traction energy is partially provided by the discharging power of WESS and partially by a substation if needed.

Substation Energy Consumption
The realtime power of energy consumption from substation z is determined by the surplus traction energy in electricity supply interval z and discharging power of WESS z: Thus, the whole day total energy consumption from all substations (denoted as E s ) in a subway system is obtained as:

Dual-Objective Substation Energy Consumption Optimization Model
For energy optimization problems, it is necessary to consider both profit and costs paid for that. Note that total traction energy demand is a constant in this work, thereby the profit is maximized if substation energy consumption is minimized. The cost is determined by the energy capacity of WESS in this work. It is correspondingly determined by the total number of BESMs. Thus, a mathematical model of the dual-objective substation energy-consumption optimization problem is formulated as follows:

•
Problem P: x n ≤ x i,n ≤ x n , i ∈ {1, 2, · · · , I}, n ∈ {1, 2, · · · , N − 1}, x i,n ∈ Z + , i ∈ {1, 2, · · · , I}, n ∈ {1, 2, · · · , N − 1}, where (26)  Objective (26) aims to minimize substation energy consumption. Objective (27) aims to minimize the total number of BESMs, which is a measurement of the cost paid for energy conservation. Constraint (28) guarantees the lower and upper limits of headway time. The lower limit is determined by safety requirements of a train signaling system and operation cost, and the upper one is determined by required service quality. Constraint (29) ensures the lower and upper limits of dwell time for a train at a platform. They are determined by the passenger flow and its potential variance. Constraint (30) guarantees that the operation time duration of a subway system is fixed. (31) guarantees the punctuality of each train travel, which is an important measurement of service quality in a subway system. (32) and (33) ensure that headway time and dwell time are positive integers, respectively. (34) ensures that the number of BESMs in a WESS is non-negative.

Resolution Method
There are several techniques to solve a multi-objective optimization problem [32,[37][38][39][40]. The most popular and straightforward one is a weighted-sum method. It converts the former into a single-objective optimization problem by using a linear weighted sum formulation that combines all the objectives. Then, the single-objective optimization problems can be solved by using an analytical method, or commercial software (e.g., CPLEX), or an intelligent optimization algorithm. A weighted-sum method together with GA are frequently used in solving dual-objective problems in subway systems [8,24,41]. However, the weights are hard to decide sometimes. Moreover, this method is inappropriate if not all objectives can be represented via a linear combination, and it is ill-suited for a multi-objective problem with non-convex objective space [32,42].
Another well-known technique to solve a multi-objective optimization problem is an ε-constraint method. This method is first introduced by Haimes et al. [43], some good application examples can also be found in [38,44]. In solving a multi-objective optimization problem, it aims to optimize only one primary objective each time, and the other objectives are transformed into constraints. By gradually changing the constraints with a step length of ε, a series of optimal results can be obtained. Thus, a set of Pareto optimal solutions can be found for the original multi-objective optimization problem. It is very convenient to be applied in the resolution of a dual-objective optimization problem. Since its first use to solve a dual-objective shortest path problem [45], it has been successfully applied to solve many problems [32,38,[45][46][47][48][49]. Note that the dual-objective substation energy-consumption optimization problem presented in this work is a dual-objective problem with discrete decision variables (integers), and objective (27) is obviously linear. Thus, its Pareto front is able to be obtained by using an -constraint method in theory [32,33]. Therefore, an -constraint method is designed to transform the dual-objective substation energy-consumption optimization problem into several single-objective optimization problems.
To solve these single-objective optimization problems, an improved artificial bee colony algorithm is designed. By combining the -constraint method and the improved artificial bee colony algorithm, the non-dominated solutions of the dual-objective substation energy-consumption optimization problem are obtained.

-Constraint Method
As our main objective is to minimize substation energy consumption in a subway system, we treat objective (26) in problem P as the primary objective, and objective (27) is converted into a constraint. In this way, the initial dual-objective substation energy-consumption optimization problem can be transformed into the following single-objective optimization problems by applying an -constraint method. The detailed procedure of applying the proposed -constraint method is shown in Algorithm 1, and the obtained single-objective optimization problems are listed as follows. Note that only one objective in the original dual-objective optimization problem is kept in each single-objective optimization problem, the other objective is either dismissed or converted into a constraint. The order to solve the single-objective optimization problems is also important, as the solution of the previous problem may be used in the following problems.

Constraints (28)-(34).
Note that the objective f 1 () in this problem is corresponding to objective (26) in problem P, and objective (27) in problem P is not considered in this problem. By solving problem P 1 , we obtain the minimal value of E s denoted as f 0 1 . It represents the minimal substation energy consumption when sufficient WESS is installed.

Constraints (28)-(34).
Note that the objective f 2 () in this problem is corresponding to objective (27) in problem P, and objective (26) in problem P is not considered in this problem. By solving problem P 2 , we obtain the minimal value of K denoted as Ω . By using an analytical method, it is easy to obtain that Ω = 0. It represents the minimal number of BESM installed in a subway system, when substation energy consumption is not considered. Note that objective vector f 0 1 , Ω is the ideal point of problem P. It represents the lower limit of the Pareto front, which is unreachable.
Note that objective (26) in problem P is transformed into constraint (35) in this problem. By solving this problem, we obtain its optimal result denoted as 0 . It represents the minimal number of BESMs needed while substation energy consumption is optimally minimized. Thus, ( f 0 1 , 0 ) is a non-dominated point of problem P. Substation energy consumption cannot be reduced by further increasing the number of BESMs.
Note that objective (27) in problem P is transformed into constraint (36) in this problem. By solving this problem, we obtain the optimal result of E s denoted as f Ω 1 . Note that ( f Ω 1 , Ω ) is also a non-dominated point of problem P. It represents the minimal substation energy consumption when a minimal number of WESS is installed. As Ω = 0, f Ω 1 is the minimal substation energy consumption when a timetable is optimized and no WESS is installed.

•
Problem P 5 : H, X, κ) , Constraints (28)- (34). Note that objective (27) in problem P is transformed into constraint (37) in this problem. Problem P 5 aims to minimize substation energy consumption with a limited number of BESMs. By solving it, we obtain the optimal result denoted f ω 1 . As ω can be assigned with different values, problem P 5 represents a series of problems. With ω decreasing from 0 to Ω , we obtain a series of optimal results of problem P 5 . Each pair of ( f ω 1 , ω ) is a non-dominated point of problem P. Thus, f ω 1 , ω |ω ∈ {0, 1, · · · , Ω} constitutes the Pareto front of problem P. Note that Ω = 0 and K is a non-negative integer in this work.
Algorithm 1 Procedure of applying the -constraint method Input: Mathematical model of problem P

Output:
f ω 1 , ω %Non-dominated points on Pareto front 1: Transform problem P into problems P 1 , P 2 , P 3 , P 4 and P 5 ; 2: Solve problems P 1 , P 2 , P 3 and P 4 sequentially to obtain the optimal results f 0 1 , Ω , 0 and f Ω 1 , respectively; 3: Set Ω = 0 ; 4: for ω = 1 to Ω − 1 do 5: Set ω = ω−1 − 1; 6: Solve problem P 5 to obtain the optimal result f ω 1 ; 7: end for 8: Remove dominated points from f ω 1 , ω |ω ∈ {0, 1, · · · , Ω} and output the remaining as the Pareto front of problem P; According to [50], one drawback of an -constraint method is that an improper selection of may result in a formulation with no feasible solution in general. Fortunately, the second constraint in our dual-objective optimization problem, i.e., (27) in problem P is a linear function, and its decision variables are all non-negative integers. Thus, its objective value K is definitely a non-negative integer too, i.e., the changing step length of ω can be set to 1 naturally and the complete set of Pareto optimal solutions is able to be obtained correspondingly, providing the exact solution for each single-objective optimization problem can be obtained. In addition, for any non-negative value of K, there is always at least one feasible solution for each κ z satisfying K = ∑ Z z=1 κ z . e.g., one possible solution is κ 1 = K and κ z = 0, ∀z ∈ {2, 3, · · · , Z}. Furthermore, whatever value is K, the timetable currently used in a subway system is always a feasible solution to (26). Thus, the drawback of potentially no feasible solution in an -constraint method does not exist for our dual-objective optimization problem.
The results obtained by applying an -constraint method are illustrated in Figure 6. Point A is the ideal point and B is the nadir one of the dual-objective optimization problem. Points C and E are the end points on the Pareto front. Note that point D represents a typical point on the Pareto front. With ω changes from 0 to Ω , f ω 1 is changed from f 0 1 to f Ω 1 correspondingly, thus point D can be any point on the green curve in Figure 6. The procedure of applying an -constraint method to transform a dual-objective optimization problem P into several single-objective ones (i.e., problems P 1 to P 5 ) is detailed as follows. Note that the first objective in problem P (i.e., substation energy consumption) is treated as the primary objective.
1. Obtain problems P 1 and P 2 respectively, where problem P 1 dismisses the second objective in problem P and problem P 1 dismisses the other one. The constraints in problem P are kept unchanged in both problems P 1 and P 2 . 2. By solving problem P 1 , we obtain the minimal value of the first objective (i.e., substation energy consumption) without consideration of the other objective (WESS cost). The result is illustrated as f 0 1 in Figure 6. Similarly, by solving problem P 2 , we obtain the minimal value of the second objective without consideration of the other one. The result is illustrated as Ω . 3. Obtain problem P 3 by adding constraint (35) into obtaining problem P 2 . Note that the constraint is the first objective function in problem P and f 0 1 is the optimized result of problem P 1 . 4. By solving problem P 3 , we obtain the minimal value of WESS size needed (denoted as 0 ) whilst maintaining minimum substation energy consumption. Note that ( f 0 1 , 0 ) is a non-dominated point of problem P, shown as point C in Figure 6. 5. Obtain problem P 4 by adding constraint constraint (36) into obtain problem P 1 . Note that the constraint is the second objective function in problem P and Ω is the optimized result of problem P 2 . 6. By solving problem P 4 , we obtain the minimal value of substation energy consumption (denoted as f Ω 1 ) whilst the least size of WESS is installed. Note that ( f Ω 1 , Ω ) is also a non-dominated point of problem P, shown as point E in Figure 6. 7. Obtain problem P 5 by adding constraint constraint (37) into obtain problem P 1 . Note that the constraint is the second objective function in problem P and ω is a variable which changes from 0 to Ω with a predefined step length. For different values of ω , it becomes a different problem. Thus, problem P 5 represents a series of optimization problems actually. 8. By solving problem P 5 with each value of ω , we obtain a set of minimal values of substation energy consumption (denoted as f ω 1 ) whilst certain size of WESS (i.e., totally ω number of BESMs) is installed. Note that each ( f ω 1 , ω ) is a non-dominated point of problem P, shown as point D in Figure 6. Also note that point D actually represents a set of points in the green curve, as ω is a changing from 0 to Ω .
The proposed IABC starts with an initial set of feasible solutions, where one of them is the currently used timetable with no WESS, and the others are randomly generated. There are three kinds of bees, i.e., employed, onlooker and scout bees, used in an optimization process. In each iteration, each employed bee is employed at a particular solution, and finds a neighbor one via a local search operator, which is randomly chosen from swap, insertion, mutation and crossover operators; each onlooker bee chooses a solution by spinning a roulette wheel, and then finds a neighbor one with the same procedure as an employed bee; each scout bee randomly generates a feasible solution, to enhance IABC's global search ability. The best solution found by all bees is kept as an initial one in the next iteration. When the iteration count reaches a predefined threshold, the algorithm restarts. Note that all the solutions, except the best one found, are replaced by randomly generated ones. IABC terminates when the restart count reaches a predefined limit. Then, the best solution found is output as an optimal result. The main procedure of IABC is shown in Algorithm 2. The detailed procedure of IABC can also be found in our previous work in [7]. n 1 , n 2 and n 3 are the numbers of employed, onlooker and scout bees, respectively; n 4 is the total number of all bees, i.e., n 4 = n 1 + n 2 + n 3 ; p 1 , p 2 , p 3 and p 4 are the probability of swap, insertion, mutation and crossover operators been chosen, respectively; M 1 and M 2 are predefined iteration count limits; S * is a specified initial solution; S 1 , S 2 , · · · , and S n 4 +1 are one of n 4 + 1 solutions, respectively; S and F are the optimal solution found and its fitness value, respectively. for k = 1 to n 4 + 1 do 9: Calculate the fitness value F k for each solution S k ; 10: end for 11: Sort F 1 , F 2 , · · · , F n 4 +1 in 'descend' order and pick out the first n 1 number of elements to generate a list F; 12: Generate a listS containing the n 1 number of solutions corresponding to F; 13: Set S n 4 +1 =S (1) ; %record the best solution 14: Set F n 4 +1 = F (1) ; %Fitness value of the best solution %employed bee phase 15: for k = 1 to n 1 do 16: Set solution B =S (k) ; 17: Generate S k from B via a local search operator randomly chosen from swap, insertion, mutation and crossover; 18: end for %onlooker bee phase 19: for k = 1 to n 2 do 20: Set B as a solution inS by roulette wheel selection; 21: Generate S n 1 +k from B via a local search operator randomly chosen from swap, insertion, mutation and crossover; 22: end for %scout bee phase 23: for k = 1 to n 3 do 24: Randomly generate a feasible solution S n 1 +n 2 +k ; 25: end for 26: end for 27: Output S = S n 4 +1 and F = F n 4 +1 ; The currently used timetable without WESS is used as an initial solution and the idea of elitism is adopted in IABC, to ensure the optimized solution is not getting worse than it. In each iteration, all solutions are evaluated and sorted based on their fitness values. A list containing n 1 best solutions is generated accordingly. Then, each employed bee is assigned to a specific solution in the list sequentially, and a new solution is generated by applying a local search operator. The local search operator is randomly chosen from swap, insertion, mutation and crossover operators with the probability of p 1 , p 2 , p 3 and p 4 , respectively. To ensure each newly generated solution is feasible, i.e, satisfying constraints (28)-(34), its feasibility is checked immediately and repair is applied if it is infeasible. Similarly, each onlooker bee generates a new solution based on an assigned solution in the same way. The difference between employed and onlooker bees is as follows: the former is associated with a particular solution in the best solution list, and the latter randomly chooses one in the list by spinning a roulette wheel. A scout bee randomly generates a feasible solution in each iteration to ensure the global search ability of IABC. IABC converges to a local optimum in M 2 iterations. Then, the best solution found so far is kept and the other solutions are replaced by randomly generated feasible ones, to reduce the chance of IABC being trapped in a local optimum. Thereafter, a similar process repeats until the total iteration count reaches M 1 × M 2 . Finally, a near-optimal solution and its fitness value are obtained.

Experimental Results and Analysis
To show the effectiveness of the proposed method, numerical examples are conducted based on the actual data obtained from Yanfang Line in Beijing, China. The actual dwell time at each platform and running time at each section are shown in Table 1. Note that x n is dwell time at a departure platform with index n, x n and x n being its lower and upper bounds in the optimization process, and r l is the running time from platform x n to x n+1 . They are the same for all trains. The headway time and other parameters of Yanfang Line are listed in Table 2. Note that headway time between any two successive trains is identical in the current timetable. The sections in each electricity supply interval are given in Table 3. A super-capacitor is selected as a BESM in this work. Parameters about WESS are shown in Table 4.  1  30  25  35  2  129  2  30  25  35  3  98  3  30  25  35  4  117  4  30  25  35  5  135  5  25  20  30  6  139  6  30  25  35  7  84  7  30  25  35  8  128  8  30  25  35  9  141  9  30  25  35  10  136  10  30  25  35  11  124  11  30  25  35  12  83  12  30  25  35  13  140  13  25  20  30  14  132  14  30  25  35  15  117  15  30  25  35  16  96  16  30  25  35 17 119  Table 4. Parameters about WESS.

Parameter Value Parameter Value Parameter Value
Parameters of IABC are set as shown in Table 5. GA is also applied to solve the single-objective optimization problems as a baseline method. The detailed processes of GA can be found in [21]. To compare IABC and GA fairly, the population size of GA is set to be 40 and its maximum iteration count is 300. The local search operators used in GA include selection, crossover and mutation. The probability of crossover and mutation are 0.8 and 0.2, respectively.  Compared with identical dwell time, more energy can be saved when dwell time at a platform varied for different trains. But regenerative energy utilization improvement is relatively limited, and it increases the solution space of an optimization problem dramatically [25]. Thus, in order to simplify the problem and keep consistent with the actual timetable, we keep the dwell time for different trains i and j at platform n identical in the experiments. i.e., x i,n = x j,n , ∀i, j ∈ {1, 2, ..., I} and n ∈ {1, 2, ..., N − 1}.
The experiments are implemented in MATLAB 2014 and runs on a notebook with Intel(R) Core(TM) i5-3210M CPU @2.50 GHz, 12 GB RAM and a Windows 7 Operating System.
The simulation results are are shown in Figure 7. Note that there are several WESS in a subway line, different WESS may have different number of BESMs. The horizontal axis in Figure 7 is the sum of the BESM numbers in each WESS. The number of BESMs equals 0 means there are no WESS installed at all. When it is greater than 0, the detailed configuration in each WESS should be obtained in the optimized solution. The vertical axis is the substation energy consumption. A black asterisk in Figure 7 denotes substation energy consumption with the current timetable (i.e., the timetable is not optimized). Note that the configuration of each WESS is set as the optimized result, as there does not exist an actual configuration of WESS for each case. A blue cross denotes an optimized result obtained by GA where timetable and WESS are simultaneously optimized, a green circle denotes an optimized result of our dual-objective optimization problem obtained by IABC, and a red "+ " in a green circle denotes a non-dominated point on the Pareto front of problem P. From the experimental results, it is easy to see that: (1) Substation energy consumption decreases with the total size of WESS increases, and the decreasing speed gradually slows down. (2) With the same size of WESS installed, timetable optimization further reduces substation energy consumption, which shows the effectiveness of the integration optimization of timetable and WESS. (3) IABC and GA can both to solve the single-objective optimization problems, and IABC performs better than GA, as the results obtained by IABC dominate those obtained by GA. (4) A set of Pareto optimal solutions of the dual-objective optimization problem is obtained by applying the proposed -constraint method and IABC. The non-dominated points are diverse and well distributed over the Pareto front.
Set the current timetable with no WESS installed, which is the actual conditions in the experimental subway line, as a basis for comparison, substation energy consumption reduces by (56.293 − 52.180)/56.293 × 100% = 7.31% when timetable is optimized without WESS invested, optimization of WESS with total size as 37 BESMs reduces substation energy consumption by (56.293 − 50.663)/56.293 × 100% = 10.00% when timetable is optimized, and the decrement comes to (56.293 − 46.331)/56.293 × 100% = 17.70% when timetable is optimized together with 37 BESMs installed as WESS. Thus, timetable optimization and WESS installation can both improve energy conservation, and the integration of them reaches the greatest extent. Note that the results we obtained is a set of Pareto optimal solutions, each element of them reduces substation energy consumption to different degrees with different WESS investment cost. The experimental results are helpful for optimal decision making. Based on relationship between energy conservation profit and WESS cost, decision makers can easily choose a preferred solution according to their particular needs, e.g., the solution with the lowest cost, the one with the greatest profit, or that with the best cost-effectiveness.
As an illustration of one optimized solution, the optimized timetable and configuration of each WESS when the total number of BESMs as K = 37 are given in Figure 8. The upper part shows the difference between the optimized headway time and that in the current timetable; the lower left part shows the optimized dwell time and the current dwell time; in addition, the lower right part shows the numbers of BESMs in each WESS, note that there are four electricity supply intervals in this subway line, which means four WESSs are installed. It is seen that the headway time and dwell time are slightly changed from the current timetable, and the configuration for each WESS is proposed. Note that there are several electricity supply intervals in a subway system, and one WESS is installed in an electricity supply interval. Thus, there are several WESSs in a subway system (e.g., WESS 1 to z). As each WESS is independent from others, they can have different capacities (i.e., number of BESMs) in theory. The profiles of the state-of-charge of each WESS are shown in Figure 9. The peak values of all WESS are between 90-100%. There is some residual capacity, but it is very small. Thus, this kind of WESS configuration is reasonable for regenerative energy utilization, i.e., if the total number of BESMs is reduced, some regenerative energy may not be able to be utilized; and if the number was increased, the increased energy capacity would be a waste of money. A real-world subway line is often subject to minor real-time perturbations which affects the adherence to the timetable [7,17]. The presence of unexpected disturbances brings uncertainties to scheduling optimization problems [65,66], which may affect the effectiveness of the optimized results. To make it clear, the problem and method proposed in this work are intended to be used offline indeed, which means the optimization problem is solved only once before trains go into operation and the potential uncertainties are not particularly considered. However, the offline method is still of great significance, as the size of each WESS is hard to be changed over time; it should be decided at the first stage of design.
To solve the potential uncertainties in real-world application, there are two alternative ways. One is to adopt a real-time optimization method [65,66], e.g., using a receding horizon principle to optimize the timetable iteratively; the other is to evaluate the robustness of the optimized results on disturbances and subject to it if acceptable (i.e., the optimized results are not too sensitive to disturbances). It is relatively complicated to introduce an online optimization method in this work, so we leave it as our future work direction and adopt the second method here.
To evaluate the robustness of the optimized timetable, we conduct the following experiments by adding random noise to it. In the experiment, we traverse the headway time and dwell time vector and add a random number to each element in it. The random number is randomly chosen from set {−δ, 0, δ} (s), which represents an actual noise. As a comparison, we add this king of noise both on the current timetable and the optimized one, and δ is assigned with different values to compare their effects on energy conservation. To be justified, we run the experiment 100 times with each value of δ and their average substation energy consumption is obtained, respectively. The experimental results with δ varies from 0 to 3 are given in Table 6. Note that δ = 0 represents no noise added on the timetable. It is seen that the energy saving ratio slightly decreases when the disturbance on train operations becomes larger. However, even with 3-s of noise, the energy saving ratio is still around 17% for the optimized timetable over the current one. Thus, the optimized result is robust enough for energy conservation under disturbances. The proposed offline optimization method is obviously acceptable to this particular problem.

Conclusions
We propose an integration of timetable optimization and WESS to improve regenerative energy utilization, thus to reduce substation energy consumption in a subway system. A dual-objective optimization model is formulated to simultaneously minimize substation energy consumption and WESS investment cost accordingly. To solve the dual-objective optimization problem, an -constraint method is first designed to transform it into several single-objective optimization problems; then, an improved artificial bee colony algorithm is designed to solve them sequentially. By combining the proposed -constraint method and improved artificial bee colony algorithm, a set of Pareto optimized solutions for the original dual-objective optimization can be obtained. Experiments based on the actual data from a subway system in China are conducted to illustrate the procedure to use them and their effectiveness. The results are useful for decision makers in a subway system, based on which they can easily make a sensible trade-off between energy saving profit and WESS investment cost according to their practical needs. In addition, the proposed improved artificial bee colony algorithm is also compared with a commonly used genetic algorithm during the experiments to prove its effectiveness.
The proposed model and resolution method can be used in any subway system which is equipped with regenerative braking systems, to maximize energy saving with the minimum investment cost paid for that. Although the optimized results are robust enough to random disturbances, one limit of the proposed method is that it is mainly applicable for offline optimization now. Thus, part of our future work will be the study of an online optimization method, to better handle the parameter uncertainty problems caused by unexpected disturbances to train operations. In addition, train speed profiles at sections also affect regenerative energy utilization and substation energy consumption. Thus, in our future research, we also plan to integrate energy-efficient train operations into our model to further reduce substation energy consumption.

Conflicts of Interest:
The authors declare no conflict of interest.