# A Sparse Optimization Approach for Energy-Efficient Timetabling in Metro Railway Systems

Authors: Ziyan Luo; Xiaoyu Li; Naihua Xiu
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

In this paper, we propose a sparse optimization approach to maximize the utilization of regenerative energy produced by braking trains for energy-efficient timetabling in metro railway systems. By introducing the cardinality function and the square of the Euclidean norm function as the objective function, the resulting sparse optimization model can characterize the utilization of the regenerative energy appropriately. A two-stage alternating direction method of multipliers is designed to efficiently solve the convex relaxation counterpart of the original NP-hard problem and then to produce an energy-efficient timetable of trains. The resulting approach is applied to Beijing Metro Yizhuang Line with different instances of service for case study. Comparison with the existing two-step linear program approach is also conducted which illustrates the effectiveness of our proposed sparse optimization model in terms of the energy saving rate and the efficiency of our numerical optimization algorithm in terms of computational time.


---

## Body

## 1. Introduction

With the rapid increase in industry and population in cities, the metro railway system has been playing a more and more important role in urban economic and social development for metropolis due to its high loading capacity and low pollution. Despite the inherent efficiency, the energy consumption of metro railway systems is of a huge magnitude and the efficient energy management for such systems is then of global importance. It is known that the energy consumption of trains is largely influenced by the driving strategy and the arrivals and departures in the timetable, which inspires the extensive study on energy-efficient train control and energy-efficient train timetabling. A careful treatment on such topics can be found in the nice survey paper by Scheepmaker et al. in [1], with focuses on reviewing the mathematical problem formulations and solution approaches in the literature. Related work on energy optimization for rail traffic systems such as tramway systems has also been well studied; see, e.g., [2, 3].Among all the involved energy saving strategies in metro railway systems, an energy-efficient timetable with regenerative braking which schedules the arrival time and the departure time of each train to and from the platforms it visits for the minimization of energy consumption and the maximization of produced regenerative energy utilization, has gained considerable academical attention in recent years. In 2012, Peña-Alcaraz et al. [4] proposed a large-scale mixed integer programming model to maximize the utilization of regenerative braking energy by means of maximizing the total duration of all possible synchronization processes of acceleration and regenerative braking between all possible train pairs. The simulation model is applied to line three of the Madrid metro system and the energy savings of the optimized timetable are reported to be about 7% on average. By merely choosing those train pairs suitable for regenerating energy transfer other than all possible train pairs in [4], Das Gupta et al. developed a more tractable mixed integer programming model in [5]. In 2013, Yang et al. [6] presented a cooperative integer programming model to maximize the time overlaps of nearby accelerating and braking trains and applied a genetic algorithm to approximately solve the resulting optimization model. For the case study on Beijing Metro Yizhuang Line, they reported 22.1% and 15.2% increase of effective time overlaps during peak hours and off-peak hours, respectively. A two-objective timetable optimization model for maximizing the regenerative braking energy with consideration of the passenger waiting time was developed by Yang et al. in [7] where the genetic algorithm was again adopted to solve the optimization model. An integrated energy-efficient operation model was proposed by Li & Lo in [8] with the speed profile part to minimize the net energy consumption and the timetabling part with headway constraints to synchronize the accelerating and regenerative braking for the reuse of regenerative energy. About 25% energy savings were reported for Beijing Metro Yizhuang Line with the headway between trains being 90 seconds. Li & Lo further developed an integrated energy-efficient timetable and speed profile optimization model to determine the cycle time, the headway time, and the speed profiles for a metro line to minimize the energy consumption [9]. The resulting mathematical model is a strictly convex quadratic programming problem with simple bound constraints and an analytical approach based on the first-order optimality condition (i.e., the KKT condition). An integrated timetable and speed profile optimization model with Taylor approximation were also considered by Yang et al. in [10] where the active set method was employed to handle the resulting strictly convex quadratic programming with linear constraints. Under the assumption that trains are operating in the optimal speed profiles, Yang et al. [11] proposed a mixed integer programming model aiming at synchronizing the arrivals and departures of the trains such that regenerative braking energy can be used as much as possible during the whole day service. A two-step linear programming model with the first step to minimize the total energy consumed by all trains and the second step to maximize the utilization of regenerative energy produced by braking trains for a whole day service in metro railway networks was proposed by Das Gupta et al. in [12]. To ensure the transfer of maximum possible regenerative energy from braking trains to accelerating trains, a time overlap vector was introduced which requires to be as close as possible to the zero vector. To effectively quantify such a requirement, the l0-norm and the l1-norm were employed in [12]. By adopting the l1-norm heuristic for the original l0-norm (i.e., the number of nonzero components in a vector), together with epigraph approach that equivalently transforms the convex l1-norm minimization to a linear programming, the resulting relaxation problem in the second step was efficiently solved by Gurobi Optimizer. Numerical experiments on different instances of service PES2-SFM2 of line 8 of Shanghai Metro network spanning a full service period of one day were conducted and a significant reduction in effective energy consumption with the worst case being 19.27% was reported by applying the proposed approach.The originall0-norm minimization problem proposed in [12] is well-known to be NP-hard due to the discontinuity and nonconvexity of the l0-norm. Such a type of optimization problems is coined by sparse optimization in the optimization community and has a wide range of applications in compressed sensing, signal processing, statistical regression, financial portfolio, and so on. Various tractable relaxation strategies are proposed to replace the l0-norm, such as the popular convex surrogate l1-norm [13], and some weighted l1-norm variants [14]. Numerical algorithms are also developed for the resulting relaxation models, such as the greedy-type methods [15], the first-order optimization methods [16], and second-order methods for some special structured models [17].Inspired by the development of sparse optimization methods, in this paper, a sparse optimization model with thel0-norm and the squared Euclidean norm as the objective function is built to maximize the utilization of the regenerative braking energy. Motivated by the adaptive LASSO model in the high-dimensional statistical regression [14], we propose a weighted l1-norm relaxation scheme. A two-stage alternating direction method of multipliers (TSADMM) is then proposed with the first stage for calculating an appropriate weight vector and an initial point for the iterative algorithm in the second stage. To illustrate the effectiveness of our new mathematical model, the linear programming model in [12], a quadratic model merely considering the squared Euclidean norm, and our weighted-l1 companied with the squared Euclidean norm are applied to Beijing Metro Yizhuang Line for comparison. The numerical results on energy savings show that our model is the best among them. To demonstrate the efficiency of our proposed TSADMM, the academical version of CVX (http://cvxr.com/cvx/) is called to solve our model as well. The numerical results in computation time show that our algorithm outperforms that used in CVX.The rest of this paper is organized as follows. In Section2, the sparse optimization model for energy-efficient timetabling in metro railway systems is built. In Section 3, the weighted l1 norm relaxation model is proposed and a two-stage alternating direction method of multipliers (TSADMM) is designed to get an approximate optimal solution. Theoretical global convergence is established as well. By applying our proposed approach to Beijing Metro Yizhuang Line, we demonstrate the effectiveness of our proposed model and the efficiency of our proposed TSADMM by the numerical results in Section 4. Conclusions and future research are given in Section 5. For reference, the appendix provides the proof of the global convergence theorem for our algorithm.


## 2. Model Formulation

### 2.1. Notation

The notation that will be used throughout the paper is listed in Table1 for the convenience of the subsequent model formulation.Table 1
A list of notations.
a i t The arrival time of train t at platform i d i t The departure time of train t from platform i T The set of all trains ⋅ The cardinality of a set, i.e., the number of elements in a set [ t r _ i j t , t r ¯ i j t ] The trip time window for traint from platform i to j ϕ The set of all the crossing-overs O i j The set of all train pairs associated with the crossing-over(i,j) [ κ _ i j t t ′ , κ ¯ i j t t ′ ] The trip window for traint on the crossing-over (i,j) [ d w _ i t , d w ¯ i t ] The dwell time window for traint at platform i S The set of all platform pairs at the same interchange station C i j The set of all connection train pairs for a platform pair(i,j) [ c _ i j t t ′ , c ¯ i j t t ′ ] The connection window between the traint at the platform i and the train t′ at platform j H i , j The set of train-pairs that move along the track(i,j) [ h _ i t t ′ , h ¯ i t t ′ ] The headway window between traint and t′ arrival at or departure from platform i [ t t _ t , t t ¯ t ] The window of total travel time constraint for traint f i j Energy consumption associated with the trip(i,j)∈Atr P The set of all platform pairs powered by a same electrical substation t → Temporally closest train to the right of train t t ← Temporally closest train to the left of train t ∇ i t The relative distance fromait to regenerative alignment point Δ i t The relative distance fromdit to the consumption alignment point ε → A set ofε containing elements of the form (i,j,t,t→) ε ← A set ofε containing elements of the form (i,j,t,t←) N The set of all platforms K The set of all tracks K t The set of tracks visited by a traint N t The set of platforms visited by a train t K t r The set of all arcs associated with trip time constraints δ S ( ⋅ ) The indicator function with respect to the set S Π S ( ⋅ ) The projection operator onto the set S s i g n ( ⋅ ) The sign function ∘ The componentwise product, also known as the Hadamard product AT The transpose of the matrixA vec(X) The column vector generated by stacking all columns of the matrixX


### 2.2. Constraints

The set of constraints in the metro train networks including constraints of the trip time, the dwell time, the headway time, the total travel time, and the domain of event times will be described in this subsection, following from [12].(i)
Trip time constraint: it includes the trip time constraints with tracks and with crossing-overs. For the former one, when a traint∈T has a trip from platform i∈Nt to platform j∈Nt along the track (i,j)∈Kt, the departure time dit and the arrival time ajt are bounded by tr_ijt and tr¯ijt, respectively. The corresponding constraints are(1)tr_ijt≤ajt-dit≤tr¯ijt,∀t∈T,∀i,j∈Kt.For the latter one, one knows that a crossing-over connecting two train-lines is a special track, where a train-line is a path with some nonopposite platforms and tracks. The crossing-over is to connect the terminal platform of a train-line with the first platform of another. A traint turns around by a crossing-over and departs from the first platform to travel through the train-line after it arrives at the final platform at another train-line, and then the same train is labelled as t′ in order to treat it as two different trains [20]. Let ϕ represent the set of all crossing-overs. Consider a crossing-over (i,j)∈ϕ, we denote the set of all train pairs associated with a same physical train through the crossing-over (i,j) by Oij. The trip time to travelling through the crossing-over (i,j)∈ϕ is the time from the departure from the platform i (train is labelled as t) to the arrival at the platform j (train is labelled as t′), which has to stay within a window [κ_ijtt′,κ¯ijtt′]. And the corresponding constraints can be described as below:(2)κ_ijtt′≤ajt′-dit≤κ¯ijtt′,∀i,j∈ϕ,∀t,t′∈Oi,j.(ii)
Dwell time constraint: when a traint∈T arrives at a platform i∈Nt, it should dwell there in a time interval denoted by [dw_it,dw¯it] to make sure that passengers can get off and get on the train before it departures from the platform i. The dwell time constraint can be written as follows:(3)dw_it≤dit-ait≤dw¯it,∀t∈T,∀i∈Nt.(iii)
Connection constraint: an interchange station used by connecting trains is to make passengers transfer between trains on different lines. Any platform pair(i,j)∈S consists of two platforms at a same interchange station. For a connection train pair (t,t′)∈Cij, we assume that the train t arrives at platform i and another train t′ departs from platform j. The difference of both times can be in a window [c_ijtt′,c¯ijtt′] to maintain that passengers can get off from the former train and get on the latter. The constraint of connection can be described as follows:(4)c_ijtt′≤djt′-ait≤c¯ijtt′,∀i,j∈S,∀t,t′∈Cij.(iv)
Headway time constraint: similar to the trip time constraints, the headway time constraints includes those with tracks and with crossing-overs. In any subway system, a minimum interval between the departures and arrivals of sequenced trains on the same track is maintained, called headway time. To meet the needs of the passengers, railway networks in many cities designate an upper bound between the departures and arrivals of consecutive trains, in order to ensure passengers do not wait for a long time before the arrival of next train. DefineHij as the set of pairs of successive trains that pass by track (i,j)∈K. Consider (t,t′)∈Hij, and let [h_itt′,h¯itt′] and [h_jtt′,h¯jtt′] be the time intervals that maintained the departures and arrivals of consecutive trains t and t′ from and to platform i and platform j, respectively. For any (i,j)∈K between platform i and platform j, headway time constraint with tracks can be described as follows:(5)h_itt′≤dit-dit′≤h¯itt′,h_jtt′≤ajt-ajt′≤h¯jtt′,∀i,j∈K,t,t′∈Hij.In the other case, consider two consecutive trainst1 and t2 that go through a crossing-over (i,j)∈ϕ after departing from the final platform i of a train-line. They will be labelled as t1′ and t2′ when they arrive at the first platform j of some other train-line. Let O~ij be the set of all such two-train pairs (t1,t1′),(t2,t2′). There is an interval time window [h_it1t2,h¯it1t2] between the departures of trains t1 and t2 from platform i and a window [h_jt1′t2′,h¯jt1′t2′] between the arrivals of trains t1′ and t2′ at platform j. We can then write the headway constraints with crossing-overs as follows:(6)h_it1t2≤dit1-dit2≤h¯it1t2,h_jt1′t2′≤ajt1′-ajt2′≤h¯jt1′t2′,∀i,j∈ϕ,∀t1,t1′,t2,t2′∈O~ij.Besides, the headway time is related to the passenger demand and the number of trains. Letm¯′ be the number of trains that are in service per hour, ct be the train capacity, u be the train utility rate, and D be the passenger demand. Then, D=ct×u×m¯′ [8]. Let h be the headway time. Since h=3600/m¯′, we have(7)h=3600×ct×uD.Here we assume thatu and c are constant while the passenger demand varies. In other words, the headway time of trains changes at different periods in a day.(v)
Total travel time constraint: when a traint∈T has a trip, traversing all platforms and tracks in chronological order, it can have the total travel time between tt_t and tt¯t, to meet the needs of passengers in subway systems. Total travel time constraints can be written as follows:(8)tt_t≤dNtNtt-aNt1t≤tt¯t,∀t∈T,(vi)
Domain of event times: we set zero second as the time of the departure of the first train from a platform in a day when the railway is in service. By adding the maximum possible values of all trip times and dwell times, we can set an upper bound for arrival of the last train at the final platform. And we denote the upper bound bymT∈Z++. So, we can write the domain of the event time constraints as follows:(9)0≤ait≤mT,0≤dit≤mT,∀t∈T,∀i∈Nt.


### 2.3. Minimizing the Energy Consumption

The major energy consumption of a train is from the acceleration phase. To efficiently minimize the total energy consumption of trains in a metro railway network, Das Gupta [12] proposed an optimization model with linear inequality constraints to get a feasible timetable. The involved optimization model takes the form of(10)min∑i,j∈Ktr,t∈Tfijajt-dits.t.(1)–(9)where fij:R++→R++ with argument (ajt-dit) being the energy consumption function that is unknown in advance and some best possible affine approximate function in the sense of least-squares by using the measured energy in practice is then utilized. By denoting x≔vec(XT) with (11)X=a11a12⋯a1T-1a1Td11d12⋯d1T-1d1Ta21a22⋯a2T-1a2Td21d22⋯d2T-1d2T⋮⋮⋱⋮⋮aN1aN2⋯aNT-1aNTdN1dN2⋯dNT-1dNTthe approximate counterpart of (10) can be formulated by the following linear programming:(12)min∑i,j∈K¯trcijxj-xis.t.lij≤xj-xi≤uij,i,j∈K¯0≤xi≤mT,i=1,2,…,n,where n=2TN, K¯ is a subset of i,j∣i,j=1,2,…,n that collects all those indices according to constraints (1)-(8), and K¯tr is a subset of K¯ that collects all those indices according to constraints (1)-(2). By utilizing its optimal solution x¯, we can get a¯jt’s and d¯jt’s such that(13)ajt-dit=a¯jt-d¯it,∀t∈T,∀i,j∈Kt.(14)ajt′-dit=a¯jt′-d¯it,∀i,j∈ϕ,∀t,t′∈Oi,j.For all other constraints, lower and upper bounds are allowed to vary as described in (3)-(9).


### 2.4. Maximizing the Energy Regeneration: A Sparse Optimization Model

We adopt the speed profile by maximum accelerating, speed holding, maximum braking strategy from [18] with a multi-phase-speed-limit section, as described in Figures 1 and 2.Figure 1
Three working phases (accelerating-cruising-braking) of speed profile.Figure 2
The speed limit of a train with a multiphase section.Lett1, t2, and t3 be the time of maximum accelerating, speed hold, and maximum braking phases, and aacc and abra are the accelerating rate of the accelerating and braking phase. It is known that, with regenerative braking, kinetic energy can be converted into electricity which can be fed back to the power supply system to be used by other nearby trains in the same electric substation. Denote the conversion factor from electricity to kinetic energy as θacc and the conversion factor from kinetic energy to regenerative electricity as θbra and the resistance rate at the speed holding phase as γ. For any two adjacent stations, we adopt the formula from [18] to calculate the energy consumption for accelerating and the regenerative energy produced during braking, denoted by Econ and Ereg, respectively, in the way as below: (15)Econ=aacc2t122θacc+∫t1t1+t2γ2t-t1θaccdt=aacc2t122θacc+γ2t222θaccandEreg=θbraabra2t322.Here Econ and Ereg are both for the unit mass. A simple heuristic method is applied to approximate the power graph as shown in Figure 3, where c and c~ are the width of the rectangles for accelerating and braking and h and h~ are heights of the rectangles. Similar to those in [12], the midpoint of the width in the rectangle is called the regenerative or consumptive alignment point, where Δit represents the distance between dit and the regenerative alignment point and ∇it represents the distance between ait and the consumption alignment point.Figure 3
Use rectangles instead of energy consumption and regenerative energy.Here, we align the distance between regenerative and consumption alignment points to maximum the utilization. Before introducing such a distance, we recall the suitable train pairs as introduced in [12]. Assume that platform pairs are opposite to each other at the same electrical substation. Let P be the set of all platform pairs. For any pair (i,j)∈P and the set Ti⊆T contains all trains that arrives at, dwells, departs from platform i. When the train t∈Ti departs from the platform i, it is a best choice to find a train t→ at the platform j which are producing the regenerative braking energy using the following definition.Definition 1 (see [12]).
Consider any(i,j)∈P. For every train t∈Ti, the train t→∈Tj is called the temporally closest train to the right of t if(16)
t
→
=
arg
min
⁡
t
′
∈
x
∈
T
j
:
0
≤
a
¯
j
x
+
d
¯
j
x
/
2
-
a
¯
i
t
+
d
¯
i
t
/
2
≤
β
⁡
a
¯
i
t
+
d
¯
i
t
2
-
a
¯
j
t
′
+
d
¯
j
t
′
2where β is an empirical parameter determined by the timetable designer and is much smaller than the time horizon of the entire timetable [12]. Similarly, when a train t∈Ti arrives at platform i, it produces the regenerative energy. To find a suitable train t← departing from platform j for transferring the regenerative braking energy, we recall the following definition from [12].Definition 2 (see [12]).
Consider any(i,j)∈P. For every train t∈Ti, the train t←∈Tj is called the temporally closest train to the left of t if(17)t←=argmin⁡t′∈x∈Tj:0≤a¯it+d¯it/2-a¯jx+d¯jx/2≤β⁡a¯it+d¯it2-a¯jt′+d¯jt′2.
Denote the setsε← and ε→ whose components are (i,j,t,t←) and (i,j,t,t→), respectively. Consider any (i,j,t,t←)∈ε←. To maximize the transfer of the regenerative energy from the braking train t← by the accelerating train t, we attempt to make the term (dit+Δit-ajt→+∇jt→) to be zero or as close to zero as possible. Similarly, for (i,j,t,t→)∈ε→, our goal is to make the term (djt←+Δjt←-ait+∇it) to be zero or as close to zero as possible. Following from [12], we define an auxiliary variable y∈Rn with n=2n¯m¯, n¯=N and m¯=T by(18)y=dit+Δit-ajt→+∇jt→i,j,t,t→∈ε→,djt←+Δjt←-ait+∇iti,j,t,t←∈ε←.The components of y represent difference of the time between the alignment points of accelerating and braking. Once the running time is fixed, Δit and ∇it are two constants.
The ultimate energy consumption can be approximately calculated as follows (see Figure4):(19)Ei,j,t,t~=Ekt,con,ifyk≥c~jt~+cit2,Ekt,con-min⁡c~jt~,cit∗min⁡hit,h~jt~,ifyk≤cit-c~jt~2,Ekt,con-c~jt~+cit2-yk∗min⁡hit,hjt~,otherwise.where train t is accelerating and train t~ is braking.Figure 4
Three cases for net energy consumption function.By sorting all possible consecutive platforms for braking energy absorbing, a0-1 vector α∈Rn is introduced to indicate all those locations in y with regeneration energy usage when the components of α take the value 1. In order to make the utilization of the regenerative energy as much as possible, we then turn to minimizing the number of nonzero elements of α∘y which corresponds to minimizing α∘y0 and making nonzero elements in α∘y as close to zero as possible which corresponds to minimizing α∘y22. The resulting sparse optimization model turns out to be(20)minx,yλ1α∘y0+λ22α∘y22s.t.(3)–(9),13,14,18.where x,y∈Rn are decision variables. Similar to (12), (20) can be rewritten as(21)minx,yλ1α∘y0+λ22α∘y22:y-Ax-b=0,Fx-g=0,f-Ex≥0.where A∈Rn×n, E∈Rm×n, F∈Rl×n are (22)E=E1-E1E2-E2E3-E3E4-E4,F=Om¯×m¯E^1Om¯×2m¯⋯Om¯×2m¯Om¯×m¯Om¯×m¯Om¯×2m¯E^2⋯Om¯×2m¯Om¯×m¯⋮⋮⋮⋱⋮⋮Om¯×m¯Om¯×2m¯⋯Om¯×2m¯E^n¯-1Om¯×m¯with (23)E1=E^1Om¯×2m¯⋯Om¯×2m¯Om¯×2m¯E^2⋯Om¯×2m¯⋮⋮⋱⋮Om¯×2m¯Om¯×2m¯⋯E^n¯,E^i=-Im¯×m¯Im¯×m¯,E2=-Im¯×m¯Om¯×2n¯m¯-2m¯Im¯×m¯,E3=I2n¯m¯×2n¯m¯,E4=E~1Om¯-1×m¯⋯Om¯-1×m¯Om¯-1×m¯E~2⋯Om¯-1×m¯⋮⋮⋱⋮Om¯-1×m¯Om¯-1×m¯⋯E~2n¯,E~i=-110⋯00-11⋯0⋮⋮⋱⋱⋮00⋯-11,m=2(2n¯(m¯-1)+n¯m¯+m¯+n) and l=(n¯-1)m¯, b, g, and f are the corresponding column vectors from constraints (18), (13), (14), and (3)-(9). Apparently, y-Ax-b=0 is the reformulation of the constraint (18), Fx-g=0 is the condensed form of constraints (13) and (14), and f-Ex≥0 is the representation of all the inequality constraints distributed in (3)-(9).


## 2.1. Notation

The notation that will be used throughout the paper is listed in Table1 for the convenience of the subsequent model formulation.Table 1
A list of notations.
a i t The arrival time of train t at platform i d i t The departure time of train t from platform i T The set of all trains ⋅ The cardinality of a set, i.e., the number of elements in a set [ t r _ i j t , t r ¯ i j t ] The trip time window for traint from platform i to j ϕ The set of all the crossing-overs O i j The set of all train pairs associated with the crossing-over(i,j) [ κ _ i j t t ′ , κ ¯ i j t t ′ ] The trip window for traint on the crossing-over (i,j) [ d w _ i t , d w ¯ i t ] The dwell time window for traint at platform i S The set of all platform pairs at the same interchange station C i j The set of all connection train pairs for a platform pair(i,j) [ c _ i j t t ′ , c ¯ i j t t ′ ] The connection window between the traint at the platform i and the train t′ at platform j H i , j The set of train-pairs that move along the track(i,j) [ h _ i t t ′ , h ¯ i t t ′ ] The headway window between traint and t′ arrival at or departure from platform i [ t t _ t , t t ¯ t ] The window of total travel time constraint for traint f i j Energy consumption associated with the trip(i,j)∈Atr P The set of all platform pairs powered by a same electrical substation t → Temporally closest train to the right of train t t ← Temporally closest train to the left of train t ∇ i t The relative distance fromait to regenerative alignment point Δ i t The relative distance fromdit to the consumption alignment point ε → A set ofε containing elements of the form (i,j,t,t→) ε ← A set ofε containing elements of the form (i,j,t,t←) N The set of all platforms K The set of all tracks K t The set of tracks visited by a traint N t The set of platforms visited by a train t K t r The set of all arcs associated with trip time constraints δ S ( ⋅ ) The indicator function with respect to the set S Π S ( ⋅ ) The projection operator onto the set S s i g n ( ⋅ ) The sign function ∘ The componentwise product, also known as the Hadamard product AT The transpose of the matrixA vec(X) The column vector generated by stacking all columns of the matrixX


## 2.2. Constraints

The set of constraints in the metro train networks including constraints of the trip time, the dwell time, the headway time, the total travel time, and the domain of event times will be described in this subsection, following from [12].(i)
Trip time constraint: it includes the trip time constraints with tracks and with crossing-overs. For the former one, when a traint∈T has a trip from platform i∈Nt to platform j∈Nt along the track (i,j)∈Kt, the departure time dit and the arrival time ajt are bounded by tr_ijt and tr¯ijt, respectively. The corresponding constraints are(1)tr_ijt≤ajt-dit≤tr¯ijt,∀t∈T,∀i,j∈Kt.For the latter one, one knows that a crossing-over connecting two train-lines is a special track, where a train-line is a path with some nonopposite platforms and tracks. The crossing-over is to connect the terminal platform of a train-line with the first platform of another. A traint turns around by a crossing-over and departs from the first platform to travel through the train-line after it arrives at the final platform at another train-line, and then the same train is labelled as t′ in order to treat it as two different trains [20]. Let ϕ represent the set of all crossing-overs. Consider a crossing-over (i,j)∈ϕ, we denote the set of all train pairs associated with a same physical train through the crossing-over (i,j) by Oij. The trip time to travelling through the crossing-over (i,j)∈ϕ is the time from the departure from the platform i (train is labelled as t) to the arrival at the platform j (train is labelled as t′), which has to stay within a window [κ_ijtt′,κ¯ijtt′]. And the corresponding constraints can be described as below:(2)κ_ijtt′≤ajt′-dit≤κ¯ijtt′,∀i,j∈ϕ,∀t,t′∈Oi,j.(ii)
Dwell time constraint: when a traint∈T arrives at a platform i∈Nt, it should dwell there in a time interval denoted by [dw_it,dw¯it] to make sure that passengers can get off and get on the train before it departures from the platform i. The dwell time constraint can be written as follows:(3)dw_it≤dit-ait≤dw¯it,∀t∈T,∀i∈Nt.(iii)
Connection constraint: an interchange station used by connecting trains is to make passengers transfer between trains on different lines. Any platform pair(i,j)∈S consists of two platforms at a same interchange station. For a connection train pair (t,t′)∈Cij, we assume that the train t arrives at platform i and another train t′ departs from platform j. The difference of both times can be in a window [c_ijtt′,c¯ijtt′] to maintain that passengers can get off from the former train and get on the latter. The constraint of connection can be described as follows:(4)c_ijtt′≤djt′-ait≤c¯ijtt′,∀i,j∈S,∀t,t′∈Cij.(iv)
Headway time constraint: similar to the trip time constraints, the headway time constraints includes those with tracks and with crossing-overs. In any subway system, a minimum interval between the departures and arrivals of sequenced trains on the same track is maintained, called headway time. To meet the needs of the passengers, railway networks in many cities designate an upper bound between the departures and arrivals of consecutive trains, in order to ensure passengers do not wait for a long time before the arrival of next train. DefineHij as the set of pairs of successive trains that pass by track (i,j)∈K. Consider (t,t′)∈Hij, and let [h_itt′,h¯itt′] and [h_jtt′,h¯jtt′] be the time intervals that maintained the departures and arrivals of consecutive trains t and t′ from and to platform i and platform j, respectively. For any (i,j)∈K between platform i and platform j, headway time constraint with tracks can be described as follows:(5)h_itt′≤dit-dit′≤h¯itt′,h_jtt′≤ajt-ajt′≤h¯jtt′,∀i,j∈K,t,t′∈Hij.In the other case, consider two consecutive trainst1 and t2 that go through a crossing-over (i,j)∈ϕ after departing from the final platform i of a train-line. They will be labelled as t1′ and t2′ when they arrive at the first platform j of some other train-line. Let O~ij be the set of all such two-train pairs (t1,t1′),(t2,t2′). There is an interval time window [h_it1t2,h¯it1t2] between the departures of trains t1 and t2 from platform i and a window [h_jt1′t2′,h¯jt1′t2′] between the arrivals of trains t1′ and t2′ at platform j. We can then write the headway constraints with crossing-overs as follows:(6)h_it1t2≤dit1-dit2≤h¯it1t2,h_jt1′t2′≤ajt1′-ajt2′≤h¯jt1′t2′,∀i,j∈ϕ,∀t1,t1′,t2,t2′∈O~ij.Besides, the headway time is related to the passenger demand and the number of trains. Letm¯′ be the number of trains that are in service per hour, ct be the train capacity, u be the train utility rate, and D be the passenger demand. Then, D=ct×u×m¯′ [8]. Let h be the headway time. Since h=3600/m¯′, we have(7)h=3600×ct×uD.Here we assume thatu and c are constant while the passenger demand varies. In other words, the headway time of trains changes at different periods in a day.(v)
Total travel time constraint: when a traint∈T has a trip, traversing all platforms and tracks in chronological order, it can have the total travel time between tt_t and tt¯t, to meet the needs of passengers in subway systems. Total travel time constraints can be written as follows:(8)tt_t≤dNtNtt-aNt1t≤tt¯t,∀t∈T,(vi)
Domain of event times: we set zero second as the time of the departure of the first train from a platform in a day when the railway is in service. By adding the maximum possible values of all trip times and dwell times, we can set an upper bound for arrival of the last train at the final platform. And we denote the upper bound bymT∈Z++. So, we can write the domain of the event time constraints as follows:(9)0≤ait≤mT,0≤dit≤mT,∀t∈T,∀i∈Nt.


## 2.3. Minimizing the Energy Consumption

The major energy consumption of a train is from the acceleration phase. To efficiently minimize the total energy consumption of trains in a metro railway network, Das Gupta [12] proposed an optimization model with linear inequality constraints to get a feasible timetable. The involved optimization model takes the form of(10)min∑i,j∈Ktr,t∈Tfijajt-dits.t.(1)–(9)where fij:R++→R++ with argument (ajt-dit) being the energy consumption function that is unknown in advance and some best possible affine approximate function in the sense of least-squares by using the measured energy in practice is then utilized. By denoting x≔vec(XT) with (11)X=a11a12⋯a1T-1a1Td11d12⋯d1T-1d1Ta21a22⋯a2T-1a2Td21d22⋯d2T-1d2T⋮⋮⋱⋮⋮aN1aN2⋯aNT-1aNTdN1dN2⋯dNT-1dNTthe approximate counterpart of (10) can be formulated by the following linear programming:(12)min∑i,j∈K¯trcijxj-xis.t.lij≤xj-xi≤uij,i,j∈K¯0≤xi≤mT,i=1,2,…,n,where n=2TN, K¯ is a subset of i,j∣i,j=1,2,…,n that collects all those indices according to constraints (1)-(8), and K¯tr is a subset of K¯ that collects all those indices according to constraints (1)-(2). By utilizing its optimal solution x¯, we can get a¯jt’s and d¯jt’s such that(13)ajt-dit=a¯jt-d¯it,∀t∈T,∀i,j∈Kt.(14)ajt′-dit=a¯jt′-d¯it,∀i,j∈ϕ,∀t,t′∈Oi,j.For all other constraints, lower and upper bounds are allowed to vary as described in (3)-(9).


## 2.4. Maximizing the Energy Regeneration: A Sparse Optimization Model

We adopt the speed profile by maximum accelerating, speed holding, maximum braking strategy from [18] with a multi-phase-speed-limit section, as described in Figures 1 and 2.Figure 1
Three working phases (accelerating-cruising-braking) of speed profile.Figure 2
The speed limit of a train with a multiphase section.Lett1, t2, and t3 be the time of maximum accelerating, speed hold, and maximum braking phases, and aacc and abra are the accelerating rate of the accelerating and braking phase. It is known that, with regenerative braking, kinetic energy can be converted into electricity which can be fed back to the power supply system to be used by other nearby trains in the same electric substation. Denote the conversion factor from electricity to kinetic energy as θacc and the conversion factor from kinetic energy to regenerative electricity as θbra and the resistance rate at the speed holding phase as γ. For any two adjacent stations, we adopt the formula from [18] to calculate the energy consumption for accelerating and the regenerative energy produced during braking, denoted by Econ and Ereg, respectively, in the way as below: (15)Econ=aacc2t122θacc+∫t1t1+t2γ2t-t1θaccdt=aacc2t122θacc+γ2t222θaccandEreg=θbraabra2t322.Here Econ and Ereg are both for the unit mass. A simple heuristic method is applied to approximate the power graph as shown in Figure 3, where c and c~ are the width of the rectangles for accelerating and braking and h and h~ are heights of the rectangles. Similar to those in [12], the midpoint of the width in the rectangle is called the regenerative or consumptive alignment point, where Δit represents the distance between dit and the regenerative alignment point and ∇it represents the distance between ait and the consumption alignment point.Figure 3
Use rectangles instead of energy consumption and regenerative energy.Here, we align the distance between regenerative and consumption alignment points to maximum the utilization. Before introducing such a distance, we recall the suitable train pairs as introduced in [12]. Assume that platform pairs are opposite to each other at the same electrical substation. Let P be the set of all platform pairs. For any pair (i,j)∈P and the set Ti⊆T contains all trains that arrives at, dwells, departs from platform i. When the train t∈Ti departs from the platform i, it is a best choice to find a train t→ at the platform j which are producing the regenerative braking energy using the following definition.Definition 1 (see [12]).
Consider any(i,j)∈P. For every train t∈Ti, the train t→∈Tj is called the temporally closest train to the right of t if(16)
t
→
=
arg
min
⁡
t
′
∈
x
∈
T
j
:
0
≤
a
¯
j
x
+
d
¯
j
x
/
2
-
a
¯
i
t
+
d
¯
i
t
/
2
≤
β
⁡
a
¯
i
t
+
d
¯
i
t
2
-
a
¯
j
t
′
+
d
¯
j
t
′
2where β is an empirical parameter determined by the timetable designer and is much smaller than the time horizon of the entire timetable [12]. Similarly, when a train t∈Ti arrives at platform i, it produces the regenerative energy. To find a suitable train t← departing from platform j for transferring the regenerative braking energy, we recall the following definition from [12].Definition 2 (see [12]).
Consider any(i,j)∈P. For every train t∈Ti, the train t←∈Tj is called the temporally closest train to the left of t if(17)t←=argmin⁡t′∈x∈Tj:0≤a¯it+d¯it/2-a¯jx+d¯jx/2≤β⁡a¯it+d¯it2-a¯jt′+d¯jt′2.
Denote the setsε← and ε→ whose components are (i,j,t,t←) and (i,j,t,t→), respectively. Consider any (i,j,t,t←)∈ε←. To maximize the transfer of the regenerative energy from the braking train t← by the accelerating train t, we attempt to make the term (dit+Δit-ajt→+∇jt→) to be zero or as close to zero as possible. Similarly, for (i,j,t,t→)∈ε→, our goal is to make the term (djt←+Δjt←-ait+∇it) to be zero or as close to zero as possible. Following from [12], we define an auxiliary variable y∈Rn with n=2n¯m¯, n¯=N and m¯=T by(18)y=dit+Δit-ajt→+∇jt→i,j,t,t→∈ε→,djt←+Δjt←-ait+∇iti,j,t,t←∈ε←.The components of y represent difference of the time between the alignment points of accelerating and braking. Once the running time is fixed, Δit and ∇it are two constants.
The ultimate energy consumption can be approximately calculated as follows (see Figure4):(19)Ei,j,t,t~=Ekt,con,ifyk≥c~jt~+cit2,Ekt,con-min⁡c~jt~,cit∗min⁡hit,h~jt~,ifyk≤cit-c~jt~2,Ekt,con-c~jt~+cit2-yk∗min⁡hit,hjt~,otherwise.where train t is accelerating and train t~ is braking.Figure 4
Three cases for net energy consumption function.By sorting all possible consecutive platforms for braking energy absorbing, a0-1 vector α∈Rn is introduced to indicate all those locations in y with regeneration energy usage when the components of α take the value 1. In order to make the utilization of the regenerative energy as much as possible, we then turn to minimizing the number of nonzero elements of α∘y which corresponds to minimizing α∘y0 and making nonzero elements in α∘y as close to zero as possible which corresponds to minimizing α∘y22. The resulting sparse optimization model turns out to be(20)minx,yλ1α∘y0+λ22α∘y22s.t.(3)–(9),13,14,18.where x,y∈Rn are decision variables. Similar to (12), (20) can be rewritten as(21)minx,yλ1α∘y0+λ22α∘y22:y-Ax-b=0,Fx-g=0,f-Ex≥0.where A∈Rn×n, E∈Rm×n, F∈Rl×n are (22)E=E1-E1E2-E2E3-E3E4-E4,F=Om¯×m¯E^1Om¯×2m¯⋯Om¯×2m¯Om¯×m¯Om¯×m¯Om¯×2m¯E^2⋯Om¯×2m¯Om¯×m¯⋮⋮⋮⋱⋮⋮Om¯×m¯Om¯×2m¯⋯Om¯×2m¯E^n¯-1Om¯×m¯with (23)E1=E^1Om¯×2m¯⋯Om¯×2m¯Om¯×2m¯E^2⋯Om¯×2m¯⋮⋮⋱⋮Om¯×2m¯Om¯×2m¯⋯E^n¯,E^i=-Im¯×m¯Im¯×m¯,E2=-Im¯×m¯Om¯×2n¯m¯-2m¯Im¯×m¯,E3=I2n¯m¯×2n¯m¯,E4=E~1Om¯-1×m¯⋯Om¯-1×m¯Om¯-1×m¯E~2⋯Om¯-1×m¯⋮⋮⋱⋮Om¯-1×m¯Om¯-1×m¯⋯E~2n¯,E~i=-110⋯00-11⋯0⋮⋮⋱⋱⋮00⋯-11,m=2(2n¯(m¯-1)+n¯m¯+m¯+n) and l=(n¯-1)m¯, b, g, and f are the corresponding column vectors from constraints (18), (13), (14), and (3)-(9). Apparently, y-Ax-b=0 is the reformulation of the constraint (18), Fx-g=0 is the condensed form of constraints (13) and (14), and f-Ex≥0 is the representation of all the inequality constraints distributed in (3)-(9).


## 3. Model Relaxation and Solution Algorithm

### 3.1. Relaxation Model

Problem (21) is generally NP-hard due to the combinatorial property of the involved l0-norm [21]. By adopting the well-known convex surrogate l1-norm initiated in compressed sensing [13, 22], we get the following l1+l22 relaxation model:(24)minx,yλ1α∘y1+λ22α∘y22:y-Ax-b=0,Fx-g=0,f-Ex≥0.Inspired by Oracle property beneficial from the adaptive lasso strategy proposed by Zou et al. [14], we will further utilize a reweighted l1-norm term instead of the l1-norm in (24) and get the following relaxation counterpart:(25)minx,yλ1α∘w∘y1+λ22α∘y22:y-Ax-b=0,Fx-g=0,f-Ex≥0.where w is some appropriate positive column vector to encourage a low sparsity of y.


### 3.2. A Two-Stage Alternating Direction Method of Multipliers

By introducing an auxiliary vectorz∈Rm, (25) can be rewritten as the following equality constrained optimization problem:(26)minx,y,zλ1α∘w∘y1+λ22α∘y22+δR+mz:y-Ax-b=0,Fx-g=0,z-f+Ex=0,where δR+m(h) is the indicator function with respect to R+m which takes the value 0 if h∈R+m, and +∞ otherwise. The augmented Lagrangian function of (26) takes the form of (27)LσAx,y,z;u,v,s=λ1∑i=1nwiαiyi+λ22α∘y22+δR+mz-u,y-Ax-b-v,Fx-g-s,z-f+Ex+σ2y-Ax-b22+Fx-g22+z-f+Ex22,where σ>0 is the penalty parameter, and u∈Rn, v∈Rl, s∈Rm are Lagrangian multipliers corresponding to those three equality constraint systems in (26). The well-known alternating direction method of multipliers (ADMM) proposed by Glowinski and Marroco [23] and Gabay and Mercier [24] is then adopted for handling (26) with the iterative scheme for k=0,1,…,(28)yk+1,zk+1=argminLσAxk,y,z;uk,vk,sk∣y∈Rn,z∈Rm;xk+1=argminLσAx,yk+1,zk+1;uk,vk,sk∣x∈Rn;uk+1,vk+1,sk+1=uk,vk,sk-τσyk+1-Axk+1-b,Fxk+1-g,zk+1-f+Exk+1;where τ∈(0,(1+5)/2) is the dual stepsize. Note that LσA(x,y,z;u,v,s) is separable of y and z. Thus the first subproblem in (28) can be easily solved by the corresponding first-order optimality condition which leads to the following closed-form solution:(29)yik+1=Axk+b+ukσiifαi=0,signqimaxqi-λ1wiλ2+σ,0otherwise,with qi=1/(λ2+σ)(uk+σAxk+σb)i for all i=1,…,n, and(30)zk+1=ΠR+mf-Exk+1σsk.The update for xk+1 from the second subproblem can be easily obtained by solving the following linear system:(31)Hx=rhs,where rhs=AT(yk+1-b-1/σuk)+FT(g+1/σvk)+ET(f-zk+1+1/σsk) and H=ATA+FTF+ETE. Note that E is full column rank. Thus, ETE is positive definition and also is H. This indicates that (31) admits a unique solution:(32)xk+1=H-1rhs.To effectively reduce the computation cost resulting from the large size of the involved matrices in H, the preconditioned conjugate gradient (PCG) method is employed to numerically solve the linear system (31) instead of the direct way as in (32).It is worth mentioning that we can utilize the above ADMM method to solve anl1 problem which is exactly (26) by taking w=e and λ2=0 to a relatively low accuracy. The resulting numerical solution (x^,y^,z^) can not only be served as the initial point, but also produce the weight vector w in the way that w=1/|y^i|r+ϵ with r∈(0,+∞) and ϵ>0 sufficiently small for a low sparsity promotion when solving (26). The two-stage ADMM algorithmic framework then can be summarized in Algorithm 1.Algorithm 1:A two-stage alternating direction methods of multipliers (TSADMM).
Require: Choose an initial point x00,y00,z00,u00,v00,s00, accuracy parameters ε1,ε2>0, the parame- ters λ10, λ20, σ, r, ϵ. Ensure: ( x ∗ , y ∗ , z ∗ ); Stage I Compute x^,y^,z^,u^,v^,s^. Step 1.1 Set λ1=λ10, λ2=0, x0,y0,z0,u0,v0,s0=x00,y00,z00,u00,v00,s00, w=e and k=0; Step 1.2 Compute xk+1,yk+1,zk+1,uk+1,vk+1,sk+1 by (31), (29), (30) and (28); Step 1.3 Set k=k+1. If some criterion is satisfied to the accuracy ε1, then stop. Otherwise, go to Step 1.2. Set x^,y^,z^,u^,v^,s^=xk,yk,zk,uk,vk,sk and update w=1/|y^|r+ϵ; Stage II Solve Problem (26) and output x∗,y∗,z∗. Step 2.1 Set λ1=λ10, λ2=λ20, (x0,y0,z0)=(x^,y^,z^), (u0,v0,s0)=0 and k=0; Step 2.2 Compute xk+1,yk+1,zk+1,uk+1,vk+1,sk+1 by (31), (29), (30) and (28); Step 2.3 Set k=k+1. If some criterion is satisfied to the accuracy ε2, then set x∗,y∗,z∗= x k , y k , z k and stop. Otherwise, go to Step 2.2.


### 3.3. Global Convergence

Algorithm1 admits a global convergence property as stated in the following theorem whose proof will be established in the appendix.Theorem 3.
Let(xk,yk,zk,w) be generated from Algorithm 1. If τ∈(0,1+5/2), then the sequence {(xk,yk,zk)} converges to an optimal solution to (26).


## 3.1. Relaxation Model

Problem (21) is generally NP-hard due to the combinatorial property of the involved l0-norm [21]. By adopting the well-known convex surrogate l1-norm initiated in compressed sensing [13, 22], we get the following l1+l22 relaxation model:(24)minx,yλ1α∘y1+λ22α∘y22:y-Ax-b=0,Fx-g=0,f-Ex≥0.Inspired by Oracle property beneficial from the adaptive lasso strategy proposed by Zou et al. [14], we will further utilize a reweighted l1-norm term instead of the l1-norm in (24) and get the following relaxation counterpart:(25)minx,yλ1α∘w∘y1+λ22α∘y22:y-Ax-b=0,Fx-g=0,f-Ex≥0.where w is some appropriate positive column vector to encourage a low sparsity of y.


## 3.2. A Two-Stage Alternating Direction Method of Multipliers

By introducing an auxiliary vectorz∈Rm, (25) can be rewritten as the following equality constrained optimization problem:(26)minx,y,zλ1α∘w∘y1+λ22α∘y22+δR+mz:y-Ax-b=0,Fx-g=0,z-f+Ex=0,where δR+m(h) is the indicator function with respect to R+m which takes the value 0 if h∈R+m, and +∞ otherwise. The augmented Lagrangian function of (26) takes the form of (27)LσAx,y,z;u,v,s=λ1∑i=1nwiαiyi+λ22α∘y22+δR+mz-u,y-Ax-b-v,Fx-g-s,z-f+Ex+σ2y-Ax-b22+Fx-g22+z-f+Ex22,where σ>0 is the penalty parameter, and u∈Rn, v∈Rl, s∈Rm are Lagrangian multipliers corresponding to those three equality constraint systems in (26). The well-known alternating direction method of multipliers (ADMM) proposed by Glowinski and Marroco [23] and Gabay and Mercier [24] is then adopted for handling (26) with the iterative scheme for k=0,1,…,(28)yk+1,zk+1=argminLσAxk,y,z;uk,vk,sk∣y∈Rn,z∈Rm;xk+1=argminLσAx,yk+1,zk+1;uk,vk,sk∣x∈Rn;uk+1,vk+1,sk+1=uk,vk,sk-τσyk+1-Axk+1-b,Fxk+1-g,zk+1-f+Exk+1;where τ∈(0,(1+5)/2) is the dual stepsize. Note that LσA(x,y,z;u,v,s) is separable of y and z. Thus the first subproblem in (28) can be easily solved by the corresponding first-order optimality condition which leads to the following closed-form solution:(29)yik+1=Axk+b+ukσiifαi=0,signqimaxqi-λ1wiλ2+σ,0otherwise,with qi=1/(λ2+σ)(uk+σAxk+σb)i for all i=1,…,n, and(30)zk+1=ΠR+mf-Exk+1σsk.The update for xk+1 from the second subproblem can be easily obtained by solving the following linear system:(31)Hx=rhs,where rhs=AT(yk+1-b-1/σuk)+FT(g+1/σvk)+ET(f-zk+1+1/σsk) and H=ATA+FTF+ETE. Note that E is full column rank. Thus, ETE is positive definition and also is H. This indicates that (31) admits a unique solution:(32)xk+1=H-1rhs.To effectively reduce the computation cost resulting from the large size of the involved matrices in H, the preconditioned conjugate gradient (PCG) method is employed to numerically solve the linear system (31) instead of the direct way as in (32).It is worth mentioning that we can utilize the above ADMM method to solve anl1 problem which is exactly (26) by taking w=e and λ2=0 to a relatively low accuracy. The resulting numerical solution (x^,y^,z^) can not only be served as the initial point, but also produce the weight vector w in the way that w=1/|y^i|r+ϵ with r∈(0,+∞) and ϵ>0 sufficiently small for a low sparsity promotion when solving (26). The two-stage ADMM algorithmic framework then can be summarized in Algorithm 1.Algorithm 1:A two-stage alternating direction methods of multipliers (TSADMM).
Require: Choose an initial point x00,y00,z00,u00,v00,s00, accuracy parameters ε1,ε2>0, the parame- ters λ10, λ20, σ, r, ϵ. Ensure: ( x ∗ , y ∗ , z ∗ ); Stage I Compute x^,y^,z^,u^,v^,s^. Step 1.1 Set λ1=λ10, λ2=0, x0,y0,z0,u0,v0,s0=x00,y00,z00,u00,v00,s00, w=e and k=0; Step 1.2 Compute xk+1,yk+1,zk+1,uk+1,vk+1,sk+1 by (31), (29), (30) and (28); Step 1.3 Set k=k+1. If some criterion is satisfied to the accuracy ε1, then stop. Otherwise, go to Step 1.2. Set x^,y^,z^,u^,v^,s^=xk,yk,zk,uk,vk,sk and update w=1/|y^|r+ϵ; Stage II Solve Problem (26) and output x∗,y∗,z∗. Step 2.1 Set λ1=λ10, λ2=λ20, (x0,y0,z0)=(x^,y^,z^), (u0,v0,s0)=0 and k=0; Step 2.2 Compute xk+1,yk+1,zk+1,uk+1,vk+1,sk+1 by (31), (29), (30) and (28); Step 2.3 Set k=k+1. If some criterion is satisfied to the accuracy ε2, then set x∗,y∗,z∗= x k , y k , z k and stop. Otherwise, go to Step 2.2.


## 3.3. Global Convergence

Algorithm1 admits a global convergence property as stated in the following theorem whose proof will be established in the appendix.Theorem 3.
Let(xk,yk,zk,w) be generated from Algorithm 1. If τ∈(0,1+5/2), then the sequence {(xk,yk,zk)} converges to an optimal solution to (26).


## 4. Numerical Experiments

Numerical experiments on Beijing Metro Yizhuang Line will be conducted in this section to evaluate our proposed model and to illustrate the performance our designed TSADMM algorithm. The procedure of the numerical study is shown in Figure6. All the computational results are obtained by running Matlab (version 2016b) on a windows desktop (Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz 3.40GHz RAM 16.0 G).Beijing Metro Yizhuang Line, opened on December in 2010, is one of the 20 lines in Beijing Metro Networks. It has two train lines with 14 stations as shown in Figure5. The total length of Yizhuang Line is 23.23 km, with the average distance of two station being 1.8 km, the minimum distance being 1 km, and the maximum being 2.6 km. In the section, we calculate the energy conversion when a train has a trip under unit mass.Figure 5
Railway network considered for numerical experiment.Figure 6
The procedure of the numerical experiments.


### 4.1. Inputs and Algorithm Initialization

For (26), α∈Rn in the objective function is a 0-1 vector defined in Section 2.4. And the weight vector is denoted by w=1/|y^|r+ϵ, where y^ is an approximate solution generated by Stage I in Algorithm 1. Here we set r=0.5 and ϵ=10-20. In the equality constraint y-Ax-b=0, A∈Rn×n with n=2n¯m¯, n¯=|N|=28, m¯ the number of trains in one day service, b∈Rn with components (Δit+∇jt~)(i,j,t,t~) sorting in the appropriate way, where train t is accelerating and t~ is braking, and platforms i and j are opposite to each other, satisfying the equation i+j=29. Fx-g=0 in (26) is an equality constraint about trip times. F∈Rl×m is the coefficient matrix in (26), l=(n¯-1)m¯. The components of g∈Rl are the trip times. E∈Rm×n in z-f+Ex=0 is the coefficient matrix of inequality constraints, and m=2(2n¯(m¯-1)+n¯m¯+m¯+n). The elements in the vector f∈Rm are upper and lower bounds of dwell time constraints, headway time constraints, total travel time constraints, and domain of event times constraints. For testing purpose, elements of f are given in the following way. The bounds of dwell time constraints are given in Table 2. Let D be the passenger demand for one-day service (16 hours) and uD∈R16 be the distribution of D for each hour which is presented in Table 3. Let ct be the capacity of a train, which takes the value 1440, and u∈R16 be the train utility rate for different hours as stated in Table 3 as well. By utilizing formula (7), we can get the headway time for each different period of time (per hour). Specifically, denote the headway time vector by h∈R16. Then hk=(3600×1440×uk)/(D×ukD) for all k=1,2,…,16. Set the windows for the headway time constraints to be hi±8 for all possible train pairs in the ith hour. The total travel time constraints are related to the dwell time and trip time. The bounds of trip times are calculated by Table 4 as reference. For each train t, the bounds of trip time are more than ∑q=15sq/vq for all sections (ij),i=1,2,…,27,j=i+1. The bounds of trip time constraints are presented in Table 5. The following formulae are adopted for calculating the windows of total travel time constraints,(33)tt_=r1×∑i=128dw_i+∑i=1i≠1427tr_ij+∑i=14κ_ijtt¯=r2×∑i=128dw¯i+∑i=1i≠1427tr¯ij+∑i=14κ¯ijwhere r1=1.1 and r2=0.85. For the domain of event times constraints, we set mT=16×3600.Table 2
The windows of the dwell time constraints.
Platforms SJZ1 XC1 XHM1 JG1 YZQ1 YZWHY1 WYJ1 RJDJ1 RCDJ1 TJNL1 JHL1 CQN1 CQ1 YZHCZ1  Lower bounds 25 25 25 25 25 25 25 25 25 25 25 25 25 25 Upper bounds 45 45 45 45 45 45 45 45 45 45 45 45 45 45  Platforms SJZ2 XC2 XHM2 JG2 YZQ2 YZWHY2 WYJ2 RJDJ2 RCDJ2 TJNL2 JHL2 CQN2 CQ2 YZHCZ2  Lower bounds 25 25 25 25 25 25 25 25 25 25 25 25 25 25 Upper bounds 45 45 45 45 45 45 45 45 45 45 45 45 45 45Table 3
The train utility rate and distribution ofD for one-day (16 hours) services.
Hours Periods 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 u 0.3 0.8 0.8 0.6 0.4 0.4 0.5 0.5 0.3 0.6 0.6 0.8 0.8 0.4 0.3 0.3 u D 0.03 0.11 0.13 0.07 0.05 0.05 0.04 0.05 0.04 0.05 0.07 0.13 0.12 0.03 0.02 0.02Table 4
Speed limits of every section on Beijing Metro Yizhuang Line following from [18].
Segment Segment Segment Segment Segment Sections v 1(m/s) s 1(m) v 2(m/s) s 2(m) v 3(m/s) s 3(m) v 4(m/s) s 4(m) v 5(m/s) s 5(m) SJZ-XC 13.89 0-150 23.61 150-480 18.06 480-1161 23.61 1161-2501 16.67 2501-2631 XC-XHM 16.67 2631-2643 23.61 2643-2797 20.83 2797-3534 23.61 3534-3780 16.67 3780-3905 XHM-JG 16.67 3905-3918 23.61 3918-5808 20.83 5808-6141 23.61 6141-6271 - - JG-YZQ 16.67 6271-6281 23.61 6281-8122 16.67 8122-8254 - - - - YZQ-YZWHY 16.67 8254-8265 23.61 8265-9116 16.67 9116-9246 - - - - YZWHY-WYJ 16.67 9246-9256 23.61 9259-10655 16.67 10655-10785 - - - - WYJ-RJDJ 16.67 10785-10797 23.61 10797-11933 16.67 11933-12065 - - - - RJDJ-RCDJ 16.67 12065-12077 23.61 12077-13289 16.67 13289-13419 - - - - RCDJ-TJNL 16.67 13419-13431 23.61 13419-14649 20.83 14649-15426 23.61 15426-15624 16.67 15624-15756 TJNL-JHL 16.67 15756-15768 23.61 15768-17891 16.67 17891-18021 - - - - JHL-CQN 16.67 18021-18033 23.61 18033-19982 16.67 19982-20107 - - - - CQN-CQ 16.67 20107-20120 23.61 20120-21264 16.67 21264-21394 - - - - CQ-YZHCZ 16.67 21394-21406 23.61 21406-22569 16.67 22569-22728 - - - -Table 5
The windows of the trip time constraints.
Sections 1 2 3 4 5 6 7 8 9 10 11 12 13 14  Lower bounds 160 92 150 130 80 100 100 100 140 130 130 100 100 170 Upper bounds 190 110 170 150 110 125 120 120 170 160 145 120 120 200  Sections 15 16 17 18 19 20 21 22 23 24 25 26 27 –  Lower bounds 100 100 130 130 140 100 100 100 80 130 150 92 160 – Upper bounds 120 120 145 160 170 120 120 125 110 150 170 110 190 –The initial point(u0,v0,s0) in Algorithm 1 is chosen to be 0. To calculate the energy consumption and regenerative energy, we set the maximum acceleration aacc=0.5m/s2 in accelerating, abra=-0.8m/s2 for braking, the resistance acceleration γ=-0.05m/s2 during the speed holding phase, the conversion factor from electricity to kinetic energy θacc=0.7, and the ratio of kinetic energy converting to regenerative energy θbra=0.5 as used in [18]. For testing purpose, we approximately calculate the energy consumption for accelerating and the regenerative energy produced during braking in the following way. As presented in Section 2.4, the maximum accelerating, speed holding, maximum braking strategy from [18] is adopted and the parameters t1, t2, and t3 are set in Table 6. The parameters c and c~ in Figure 3 which represent the width of the rectangles for Econ and Ebra (in percentage of the trip time) are also given in Table 6 for approximately calculating the energy consumption for accelerating Econ and the regenerative energy produced during braking Ebra.Table 6
Parameters for calculatingEcon and Ebra.
t1 t 2 t 3 c c ~  20% 45% 35% 14% 21%


### 4.2. Stopping Criteria

To measure the accuracy of an approximate optimal solution(x,y,z) of Problem (26) from (x,y,z,u,v,s) generated by Stage II in Algorithm 1, the relative primal and dual infeasibilities will be adopted. Specifically, denote(34)ηp=maxy-Ax-b1+y,Fx-g1+g,z-f+Ex1+z,ηd=ATu-FTw-ETs1+ATu+FTw+ETs.For a given tolerance ε2>0, we will stop the tested algorithms when both ηp and ηd are less than ε2. The algorithms will also be terminated when they reach a given maximum number of iterations MaxIter2. Here we set ε2=10-3 and MaxIter2=104. For Stage I in Algorithm 1, it will be stopped when both ηp and ηd are less than ε1 or the maximum iteration reaches MaxIter1. Here we set ε1=10-1 and MaxIter1=103.


### 4.3. Experiment Results

To evaluate our model withl0-norm plus l22-norm as the object, the other two possible models are computed for comparison. The first model is the l0-norm plus l1-norm proposed by Das Gupta et al. in [12] which takes the form of(35)minxy0+λy1:y-Ax-b=0,Fx-g=0,Ex≤fwhich is handled by approximately solving its linear programming relaxation problem(36)minxeTs:y≤s,y≥-s,y-Ax-b=0,Fx-g=0,Ex≤f.The second possible model is the squared l2-norm related quadratic programming,(37)minx12y22:y-Ax-b=0,Fx-g=0,Ex≤f,where only the squared Euclidean norm is employed to measure the vanishing of y. Since we focus on the reduction of the effective energy consumption, we define the energy saving rate η as follows: (38)η=∑t=1TNEi,j,t,t~Σt=1TNEtcon.Energy savings by using these three different models are listed in Table 7 with different choices of m¯.Table 7
Results of numerical experiments on Beijing Metro Yizhuang Line: A1 stands for the linear programming model (36) solved by Gurobi; A2 stands for the quadratic programming (37) solved by CVX; A3 stands for our proposed model (26) solved by Algorithm 1.
η(%) m - A 1 A 2 A 3 357 25.54 24.52 29.92 368 23.67 23.44 28.36 380 25.56 23.01 28.54 391 24.50 23.80 28.49 403 24.78 24.17 27.47 414 27.02 26.57 30.05 426 27.73 27.07 30.52 437 27.59 26.94 30.46 449 28.26 27.54 30.84As shown in Table7, the energy savings by using our model are better than those by the other two in all the testing instances with at least 2.58 percentage higher in the energy saving rate η.Other choices for the initial trip times from [6, 9–11, 19] with m¯=368 are also considered as inputs. The comparison results in terms of the energy saving rate η are stated in Table 8.Table 8
Energy saving rates form-=368: A1, A2, and A3 are the same as in Table 7.
η(%) R e f . A 1 A 2 A 3 [19] 19.82 21.34 24.94 [9] 18.85 18.83 23.11 [10] 18.99 18.88 23.11 [11] 19.02 18.85 23.12 [6] 18.97 18.78 23.29As one can see in Table8, the utilization of regenerative energy of our approach is the best among these three models as well with the energy saving rate at least 3.60 percent higher than the other two. These evidently show the effectiveness of our proposed model.In order to demonstrate the efficiency of our proposed TSADMM algorithm, we call the academical CVX to solve ours model for comparison. The time comparisons for all the aforementioned instances are presented in Tables9 and 10, respectively.Table 9
Time comparison between CVX and TSADMM with different choices ofm-.
Alg. m - 357 368 380 391 403 414 426 437 449 time(s) CVX 47.17 47.06 54.33 56.90 31.09 50.63 56.18 70.72 71.14 TSADMM 24.35 12.80 31.06 21.26 25.39 19.36 27.61 30.80 23.66Table 10
Time comparison between CVX and TSADMM with different trip times from [6, 9–11, 19].
Alg. Ref. [19] [9] [10] [11] [6] time(s) CVX 55.61 57.86 55.45 51.25 47.81 TSADMM 20.36 29.22 28.21 27.64 31.84As can be seen in Tables9 and 10, our proposed TSADMM outperforms the CVX solver for the model (26) in terms of computation time. Particularly, for the cases of m¯ equaling 368 and 449 in Table 9, the computation time for TSADMM is less than one-third of those for CVX. Besides the theoretical global convergence of our TSADMM, the numerical convergence behaviors in terms of the infeasibility measure for all the above 14 instances are presented in Figures 7 and 8. One can see that for each instance, the infeasibility decreases rapidly as the number of iterations increases and it meets the required accuracy 10-3 within less than 4000 iterations.Figure 7
The numerical convergence behavior of TSADMM with different choices ofm¯.Figure 8
The numerical convergence behavior of TSADMM for different trip times: tr1,…,tr5 stand for the cases for different trip times proposed in [6, 9–11, 19].


## 4.1. Inputs and Algorithm Initialization

For (26), α∈Rn in the objective function is a 0-1 vector defined in Section 2.4. And the weight vector is denoted by w=1/|y^|r+ϵ, where y^ is an approximate solution generated by Stage I in Algorithm 1. Here we set r=0.5 and ϵ=10-20. In the equality constraint y-Ax-b=0, A∈Rn×n with n=2n¯m¯, n¯=|N|=28, m¯ the number of trains in one day service, b∈Rn with components (Δit+∇jt~)(i,j,t,t~) sorting in the appropriate way, where train t is accelerating and t~ is braking, and platforms i and j are opposite to each other, satisfying the equation i+j=29. Fx-g=0 in (26) is an equality constraint about trip times. F∈Rl×m is the coefficient matrix in (26), l=(n¯-1)m¯. The components of g∈Rl are the trip times. E∈Rm×n in z-f+Ex=0 is the coefficient matrix of inequality constraints, and m=2(2n¯(m¯-1)+n¯m¯+m¯+n). The elements in the vector f∈Rm are upper and lower bounds of dwell time constraints, headway time constraints, total travel time constraints, and domain of event times constraints. For testing purpose, elements of f are given in the following way. The bounds of dwell time constraints are given in Table 2. Let D be the passenger demand for one-day service (16 hours) and uD∈R16 be the distribution of D for each hour which is presented in Table 3. Let ct be the capacity of a train, which takes the value 1440, and u∈R16 be the train utility rate for different hours as stated in Table 3 as well. By utilizing formula (7), we can get the headway time for each different period of time (per hour). Specifically, denote the headway time vector by h∈R16. Then hk=(3600×1440×uk)/(D×ukD) for all k=1,2,…,16. Set the windows for the headway time constraints to be hi±8 for all possible train pairs in the ith hour. The total travel time constraints are related to the dwell time and trip time. The bounds of trip times are calculated by Table 4 as reference. For each train t, the bounds of trip time are more than ∑q=15sq/vq for all sections (ij),i=1,2,…,27,j=i+1. The bounds of trip time constraints are presented in Table 5. The following formulae are adopted for calculating the windows of total travel time constraints,(33)tt_=r1×∑i=128dw_i+∑i=1i≠1427tr_ij+∑i=14κ_ijtt¯=r2×∑i=128dw¯i+∑i=1i≠1427tr¯ij+∑i=14κ¯ijwhere r1=1.1 and r2=0.85. For the domain of event times constraints, we set mT=16×3600.Table 2
The windows of the dwell time constraints.
Platforms SJZ1 XC1 XHM1 JG1 YZQ1 YZWHY1 WYJ1 RJDJ1 RCDJ1 TJNL1 JHL1 CQN1 CQ1 YZHCZ1  Lower bounds 25 25 25 25 25 25 25 25 25 25 25 25 25 25 Upper bounds 45 45 45 45 45 45 45 45 45 45 45 45 45 45  Platforms SJZ2 XC2 XHM2 JG2 YZQ2 YZWHY2 WYJ2 RJDJ2 RCDJ2 TJNL2 JHL2 CQN2 CQ2 YZHCZ2  Lower bounds 25 25 25 25 25 25 25 25 25 25 25 25 25 25 Upper bounds 45 45 45 45 45 45 45 45 45 45 45 45 45 45Table 3
The train utility rate and distribution ofD for one-day (16 hours) services.
Hours Periods 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 u 0.3 0.8 0.8 0.6 0.4 0.4 0.5 0.5 0.3 0.6 0.6 0.8 0.8 0.4 0.3 0.3 u D 0.03 0.11 0.13 0.07 0.05 0.05 0.04 0.05 0.04 0.05 0.07 0.13 0.12 0.03 0.02 0.02Table 4
Speed limits of every section on Beijing Metro Yizhuang Line following from [18].
Segment Segment Segment Segment Segment Sections v 1(m/s) s 1(m) v 2(m/s) s 2(m) v 3(m/s) s 3(m) v 4(m/s) s 4(m) v 5(m/s) s 5(m) SJZ-XC 13.89 0-150 23.61 150-480 18.06 480-1161 23.61 1161-2501 16.67 2501-2631 XC-XHM 16.67 2631-2643 23.61 2643-2797 20.83 2797-3534 23.61 3534-3780 16.67 3780-3905 XHM-JG 16.67 3905-3918 23.61 3918-5808 20.83 5808-6141 23.61 6141-6271 - - JG-YZQ 16.67 6271-6281 23.61 6281-8122 16.67 8122-8254 - - - - YZQ-YZWHY 16.67 8254-8265 23.61 8265-9116 16.67 9116-9246 - - - - YZWHY-WYJ 16.67 9246-9256 23.61 9259-10655 16.67 10655-10785 - - - - WYJ-RJDJ 16.67 10785-10797 23.61 10797-11933 16.67 11933-12065 - - - - RJDJ-RCDJ 16.67 12065-12077 23.61 12077-13289 16.67 13289-13419 - - - - RCDJ-TJNL 16.67 13419-13431 23.61 13419-14649 20.83 14649-15426 23.61 15426-15624 16.67 15624-15756 TJNL-JHL 16.67 15756-15768 23.61 15768-17891 16.67 17891-18021 - - - - JHL-CQN 16.67 18021-18033 23.61 18033-19982 16.67 19982-20107 - - - - CQN-CQ 16.67 20107-20120 23.61 20120-21264 16.67 21264-21394 - - - - CQ-YZHCZ 16.67 21394-21406 23.61 21406-22569 16.67 22569-22728 - - - -Table 5
The windows of the trip time constraints.
Sections 1 2 3 4 5 6 7 8 9 10 11 12 13 14  Lower bounds 160 92 150 130 80 100 100 100 140 130 130 100 100 170 Upper bounds 190 110 170 150 110 125 120 120 170 160 145 120 120 200  Sections 15 16 17 18 19 20 21 22 23 24 25 26 27 –  Lower bounds 100 100 130 130 140 100 100 100 80 130 150 92 160 – Upper bounds 120 120 145 160 170 120 120 125 110 150 170 110 190 –The initial point(u0,v0,s0) in Algorithm 1 is chosen to be 0. To calculate the energy consumption and regenerative energy, we set the maximum acceleration aacc=0.5m/s2 in accelerating, abra=-0.8m/s2 for braking, the resistance acceleration γ=-0.05m/s2 during the speed holding phase, the conversion factor from electricity to kinetic energy θacc=0.7, and the ratio of kinetic energy converting to regenerative energy θbra=0.5 as used in [18]. For testing purpose, we approximately calculate the energy consumption for accelerating and the regenerative energy produced during braking in the following way. As presented in Section 2.4, the maximum accelerating, speed holding, maximum braking strategy from [18] is adopted and the parameters t1, t2, and t3 are set in Table 6. The parameters c and c~ in Figure 3 which represent the width of the rectangles for Econ and Ebra (in percentage of the trip time) are also given in Table 6 for approximately calculating the energy consumption for accelerating Econ and the regenerative energy produced during braking Ebra.Table 6
Parameters for calculatingEcon and Ebra.
t1 t 2 t 3 c c ~  20% 45% 35% 14% 21%


## 4.2. Stopping Criteria

To measure the accuracy of an approximate optimal solution(x,y,z) of Problem (26) from (x,y,z,u,v,s) generated by Stage II in Algorithm 1, the relative primal and dual infeasibilities will be adopted. Specifically, denote(34)ηp=maxy-Ax-b1+y,Fx-g1+g,z-f+Ex1+z,ηd=ATu-FTw-ETs1+ATu+FTw+ETs.For a given tolerance ε2>0, we will stop the tested algorithms when both ηp and ηd are less than ε2. The algorithms will also be terminated when they reach a given maximum number of iterations MaxIter2. Here we set ε2=10-3 and MaxIter2=104. For Stage I in Algorithm 1, it will be stopped when both ηp and ηd are less than ε1 or the maximum iteration reaches MaxIter1. Here we set ε1=10-1 and MaxIter1=103.


## 4.3. Experiment Results

To evaluate our model withl0-norm plus l22-norm as the object, the other two possible models are computed for comparison. The first model is the l0-norm plus l1-norm proposed by Das Gupta et al. in [12] which takes the form of(35)minxy0+λy1:y-Ax-b=0,Fx-g=0,Ex≤fwhich is handled by approximately solving its linear programming relaxation problem(36)minxeTs:y≤s,y≥-s,y-Ax-b=0,Fx-g=0,Ex≤f.The second possible model is the squared l2-norm related quadratic programming,(37)minx12y22:y-Ax-b=0,Fx-g=0,Ex≤f,where only the squared Euclidean norm is employed to measure the vanishing of y. Since we focus on the reduction of the effective energy consumption, we define the energy saving rate η as follows: (38)η=∑t=1TNEi,j,t,t~Σt=1TNEtcon.Energy savings by using these three different models are listed in Table 7 with different choices of m¯.Table 7
Results of numerical experiments on Beijing Metro Yizhuang Line: A1 stands for the linear programming model (36) solved by Gurobi; A2 stands for the quadratic programming (37) solved by CVX; A3 stands for our proposed model (26) solved by Algorithm 1.
η(%) m - A 1 A 2 A 3 357 25.54 24.52 29.92 368 23.67 23.44 28.36 380 25.56 23.01 28.54 391 24.50 23.80 28.49 403 24.78 24.17 27.47 414 27.02 26.57 30.05 426 27.73 27.07 30.52 437 27.59 26.94 30.46 449 28.26 27.54 30.84As shown in Table7, the energy savings by using our model are better than those by the other two in all the testing instances with at least 2.58 percentage higher in the energy saving rate η.Other choices for the initial trip times from [6, 9–11, 19] with m¯=368 are also considered as inputs. The comparison results in terms of the energy saving rate η are stated in Table 8.Table 8
Energy saving rates form-=368: A1, A2, and A3 are the same as in Table 7.
η(%) R e f . A 1 A 2 A 3 [19] 19.82 21.34 24.94 [9] 18.85 18.83 23.11 [10] 18.99 18.88 23.11 [11] 19.02 18.85 23.12 [6] 18.97 18.78 23.29As one can see in Table8, the utilization of regenerative energy of our approach is the best among these three models as well with the energy saving rate at least 3.60 percent higher than the other two. These evidently show the effectiveness of our proposed model.In order to demonstrate the efficiency of our proposed TSADMM algorithm, we call the academical CVX to solve ours model for comparison. The time comparisons for all the aforementioned instances are presented in Tables9 and 10, respectively.Table 9
Time comparison between CVX and TSADMM with different choices ofm-.
Alg. m - 357 368 380 391 403 414 426 437 449 time(s) CVX 47.17 47.06 54.33 56.90 31.09 50.63 56.18 70.72 71.14 TSADMM 24.35 12.80 31.06 21.26 25.39 19.36 27.61 30.80 23.66Table 10
Time comparison between CVX and TSADMM with different trip times from [6, 9–11, 19].
Alg. Ref. [19] [9] [10] [11] [6] time(s) CVX 55.61 57.86 55.45 51.25 47.81 TSADMM 20.36 29.22 28.21 27.64 31.84As can be seen in Tables9 and 10, our proposed TSADMM outperforms the CVX solver for the model (26) in terms of computation time. Particularly, for the cases of m¯ equaling 368 and 449 in Table 9, the computation time for TSADMM is less than one-third of those for CVX. Besides the theoretical global convergence of our TSADMM, the numerical convergence behaviors in terms of the infeasibility measure for all the above 14 instances are presented in Figures 7 and 8. One can see that for each instance, the infeasibility decreases rapidly as the number of iterations increases and it meets the required accuracy 10-3 within less than 4000 iterations.Figure 7
The numerical convergence behavior of TSADMM with different choices ofm¯.Figure 8
The numerical convergence behavior of TSADMM for different trip times: tr1,…,tr5 stand for the cases for different trip times proposed in [6, 9–11, 19].


## 5. Conclusion

In this paper we have proposed a sparse optimization model usingl0-norm and squared l2-norm as objective function so as to maximize the utilization of regenerative energy in metro railway system. To overcome the NP-hardness resulting from the l0-norm, the weighted l1-norm by invoking the adaptive lasso strategy in statistical regression has been introduced to make the relaxation counterpart computationally tractable and numerically effective. A two-stage alternating direction method of multipliers (TSADMM) has been designed to efficiently solve the proposed mathematical model and the global convergence of the algorithm has been established as well. Numerical experiments have been conducted on Beijing Metro Yizhuang Line for case study. The comparison results with the linear programming proposed by Das Gupta et al. and a quadratic programming with least-squares loss function have illustrated the effectiveness of our proposed model in terms of the energy saving rate, and the timing comparison results with the CVX toolbox have demonstrated the efficiency of our proposed TSADMM algorithm.Since the weightedl1-norm is a convex approximation of the l0-norm, the optimal solution to the (26) is a proximal optimal solution to the original sparse optimization model. How to design some more effective nonconvex relaxation scheme to get a better proximal optimal solution so as to get a better energy saving rate will be one of our future work. Meanwhile, the delay perturbations of the metro trains and the randomness of passenger demands have not been considered in our static timetabling. How to use advanced optimization techniques for handling such cases will also be our future research issue.



---

*Source: 1784789-2018-11-01.xml*