# Joint Operating Revenue and Passenger Travel Cost Optimization in Urban Rail Transit

Authors: Wenxin Li; Qiyuan Peng; Qinlin Li; Chao Wen; Yongxiang Zhang; Javad Lessan
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

Urban rail transit (URT) scheduling requires designing efficient timetables that can meet passengers’ expectations about the lower travel cost while attaining revenue management objectives of the train operators. This paper presents a biobjective timetable optimization model that seeks maximizing the operating revenue of the railway company while lowering passengers’ average travel cost. We apply a fuzzy multiobjective optimization and a nondominated sorting genetic algorithm II to solve the optimization problem and characterize the trade-off between the conflicting objective functions under different types of distances. To illustrate the model and solution methodology, the proposed model and solution algorithms are validated against train operation record from a URT line of Chengdu metro in China. The results show that significant improvements can be achieved in terms of the travel cost and revenue return criteria when implementing the solutions obtained by the proposed model.


---

## Body

## 1. Introduction

Being fast, reliable, safe, and convenient, URT has been able to provide satisfactory trip services and thus mitigate urban traffic congestion in Metropolitan cities, e.g., Tokyo [1], Beijing [2], and New York [3]. In Beijing, for example, seventeen URT lines are currently operating to serve over 10 million passengers each day [4]. Satisfying passengers’ expectation about service quality and operators’ expectations regarding economically viable return have been among the main operational challenges of managing this massive transit mode. These issues have been the focus of many research efforts and many sophisticated solutions, with a special focus on timetable optimization, have been proposed over decades. Since the price of URT tickets and thus the passenger fares do not change for short-term periods, the managerial decisions are restricted and there is more focus on reducing the operating costs. As studies on URT system operations show, train moving operations consumes more than 50% of total electrical energy during train operations [5]. However, reducing energy consumption alone may lead to a timetable with long travel times and thus diminish the service quality. Therefore, an operating timetable should consider both passengers’ point of view and operators’ objectives. In this regard, the New Haven line of the Metro-North Commuter Railroad can be mentioned as a real-world case where minimizing energy consumption in track alignment, speed limit, and schedule adherence objectives are considered to satisfy passengers and operators expectations simultaneously [6]. In another effort, to obtain energy-efficient train operations and distribute the total trip time among different sections, a numerical algorithm is proposed in [7].Since the impact of train speed on energy consumption is significant, train energy consumption mainly depends on the driving strategies and the operational timetable. Also, any change in the train timetable can affect the operating costs, including passenger travel times (costs) and ticket fares, and thus the quality of transport services. The operating revenue and travel cost are sensitive to the duration of operations defined in the timetable. This necessitates an integrated model to plan train operations such that the operating revenue and travel cost objectives are optimized. In this paper, we present a biobjective timetable optimization model that maximizes the operating revenue and minimizes passengers’ average travel cost, to optimize the (average) operating speed of trains. To solve the optimization problem, we apply a fuzzy multiobjective optimization and a nondominated sorting genetic algorithm II. Next, using a case study from Chengdu metro in China, we perform a numerical analysis to characterize the behavior of the incompatible objectives at different scenarios, to figure out dominating operational strategies for improved and efficient train services.The rest of this paper is organized as follows. The next section provides a brief review of the relevant literature on train speed and timetable optimization problems. In Section3, we analyze the passenger boarding behavior and present our biobjective optimization model. In Section 4, we apply a fuzzy multiobjective optimization method and nondominated sorting genetic algorithm II to identify the relationship between different decision variables and objective functions. In Section 5, a case study from Chengdu metro in China is presented to verify the performance of the proposed model. Finally, Section 6 comes with conclusions and directions for future research.


## 2. Literature Review

To meet the unpredictable and varying operational requirements, timetable rescheduling is the most common practice in URTs. This problem has been tackled from various theoretical and operational perspectives by practitioners. From saving energy perspective, Zhou and Xu proposed a multitrain coordinated operation optimization algorithm that considers both the buffer time and safety constraints [8]. Miyatake and Ko used three different methods, including dynamic programming, gradient method, and sequential quadratic programming, to solve the URT timetable rescheduling problem for URT operations by optimizing energy consumption [9]. Some studies have also investigated improving multiple factors such as train movement profile, passenger comfort, safety, and operation stability. For instance, Su et al. considered timetable optimization and speed profiles among successive stations for energy-efficient and optimized train operations [7]. A stochastic optimization model is proposed in [10] that redistributes the time supplements and buffer times in a given timetable, to improve the safety and operation stability of URT system. Assis and Milani analyzed the evolution of train headways and train passenger loads along URT lines and presented a methodology to optimize train timetables in URT lines [11].Regarding train operation costs and total passenger travel time, Ghoseiri et al. and Chang et al. developed multiobjective optimization models for railways, to minimize fuel costs and total travel times [12, 13]. Li et al. proposed a multiobjective train scheduling model by minimizing the train energy and carbon emission costs as well as the total travel time of passengers [14]. They applied a fuzzy multiobjective optimization algorithm to solve the model. Corman et al. proposed an optimization model to generate timetables and to effectively manage the traffic in real-time, which illustrated the effects of changing trains’ speed profile in open corridors [15]. Chevrier et al. introduced the speed profiles found by the evolutionary algorithm produced by a set of solutions optimizing both the running time and energy consumption, which can be used to optimize the running of the trains [16].The literature listed above mainly has focused on the long-distance railways. In the URT system, there are similar models to deal with the fuel costs reduction and travel-time saving. A biobjective integer programming model with headway time and dwell time control and a genetic algorithm with binary encoding to find the optimal solution were conducted in [17]. The idea of optimizing metro train speed profiles was also applied to reduce energy consumption [18]. A bilevel train scheduling optimization is proposed in [19] that takes into account stop-skipping strategies and the passenger travel time, and energy consumption gave the origin-destination-dependent passenger demand.To improve the quality of metro service and reduce passenger costs, demand-sensitive orientation timetable models were presented in [20]. Moreover, a bilevel demand-oriented approach was applied to obtain a timetable for a suburban railway [21]. Considering the passenger demands, transfers, and passenger flow splitting, an event-driven train scheduling model for a URT network was proposed in [22] that concludes the nonfixed headway train schedules have a better performance.Though some researchers have focused on the train timetable scheduling problem, few of them considered the timetable scheduling problem from an operational efficiency perspective and the concerns about the lower travel cost from passengers’ side. This paper tries to fill this gap, proposing a multiobjective optimization model considering operating revenue and travel costs simultaneously.


## 3. Biobjective Optimization Model

### 3.1. Train Energy Consumption and Speed Analysis

Accelerating, coasting, and braking are the main phases of a train movement when performing running activities between successive track sections [23]. One can ignore small variations in train speed as the preventive maintenance of infrastructures and facilities keeps the operational conditions at the required level. With this simplifying assumption, train motion formula can be described as shown in (1). In this equation, we consider the basic line resistance, the track gauge, the maximum speed, and the signal systems, for instance, but we do not consider the curve and tunnel resistances. (1)a=dvdt=μffv-rvdtdx=1vwhere a is the acceleration of train operation, v is the speed of train operation, t is the time of train operation, f(v) is the maximum unit traction of the train, μf is the coefficient of the traction output ratio, and r(v) is the unit basic resistance of the train; it indicates the basic resistance of trains under unit weight, N/kN.Using this formula, when a train runs in a section, the relationship between energy consumption and the maximum speed can be derived as follows:(2)Ei,i+1′=∫0vi,i+1maxμffv·M0·vμffv-rvdvwhere vi,i+1max is the optimal maximum speed of the train in section i; it is not equal to the maximum speed limit of the train in section i; vi,i+1max≤ the maximum speed limit (the maximum speed limit of subway trains is equal to 80 km/h in China); M0 is the weight of the train. Equation (3) is another expression proposed by Gu et al., to calculate the energy consumption of a train in a section from the perspective of traction work [24].(3)Ei,i+1′′=∫0li,i+11μffv·M0·x·dxMore simplified formulations can be derived from (3) as follows: (4)Ei,i+1′′=∫0li,i+11μffv·M0·x·dx=∫0v-i,i+1t1μffv·M0·x·dx=∫0v-i,i+1t1μffv·M0·vt·v·dt=∫0v-i,i+1t1μffv·M0·vt·va·dv=∫0v-i,i+1t1μffv·M0·vt·vμffv-rv·dvthere is r(v)/f(v)≈0, because, in practice, r(v) is much smaller than f(v) and μf=1, when the train runs at the accelerating phase in a section [24]. Therefore, (2) and (4) can be reformulated into (5), given μffv-r(v)≈f(v). (5)Ei,i+1′=∫0vi,i+1maxvM0dvEi,i+1′′=∫0v-i,i+1t1vM0·vt·dvIn (4) and (5), li,i+11 is the running distance of the train at the accelerating phase in section i. We know that li,i+11=v-i,i+1t1, where v-i,i+1 is the average running speed of the train in section i; t is the quantity of the state that satisfies oΔt≈0. It refers to the time period when the speed of the train increases from v to v+oΔv. After some calculation and simplification, (5) can be further developed into (6)Ei,i+1′=∫0vi,i+1maxvM0dv=12vi,i+1max2M0Ei,i+1′′=∫0v-i,i+1t1vM0·vt·dv=13v-i,i+1t13M0tt∈0,ti,i+1Ei,i+1′=Ei,i+1′′⇒vi,i+1maxv-i,i+1=23v-i,i+1t13t=23li,i+11t·t1li,i+11=constant;t1=constant;t=oΔtOwing to the theoretical derivation of the mathematical formula, we assume the correlation coefficientε between vi,i+1max and v-i,i+1, to be (7)vi,i+1maxv-i,i+1=1ε=vi,i+1maxv-i,i+1=1ε.We approximate the value ofε using train movement data in track sections, including the optimal maximum speed, average operating speed, and the relationship between the optimal maximum speed and the average operating speed [23]. The instance is shown in Figure 1, in which y1 represents the average operating speed while y2 represents the maximum operating speed and the correlation value ε can be calculated by (8). Given this, (2) can be transformed into (9). Therefore, we get the relation between average speed and energy consumption, but this relation has the range of validity as follows.(1)
The object of the study must be urban rail transit system; the maximum speed of trains should not exceed 80km/h. And considering the comfort of passengers, the maximum train acceleration is restricted to 1 m/s2, and the maximum train deceleration is constrained as -1 m/s2.(2)
The distance between the two stations is very short (the distance between the two stations should be less than 3 kilometers in principle); the train running in the section only includes three phases: accelerating, coasting, and braking, without cruising phase.(8)ε=113×∑x=113y1xy2x(9)Ei,i+1=∫0v-i,i+1/εμffv·M0·vμffv-rv·dvFigure 1
The trend of polyliney1 and polyline y2.


### 3.2. Description of Variables and Model Assumptions

We first provide the definitions and assumptions used in this study and then describe ticket pricing schemes in URT. We denote the number of stations along a line byn, the headway between successive trains in the same direction by Tk, the operation time of a train in a full day by T0, and the running time of a train in section i by ti,i+1. Also, we set si as the dwell time of a train at station i, the passenger carrying capacity of a train as D, the occupied rate of passenger carrying capacity as δ, the total length of the Metro line as L, and, finally, the length of section i as li,i+1. Generally, two different kinds of ticket pricing schemes are used in URT. The first one depends on the train operating distance, and the second one depends on the number of pass-through sections. In Chengdu Metro, the latter method is applied for which the price of tickets can be expressed by (10). In this equation (c1,c2,…,cn) indicates different ticket prices, where j-i to the number of sections that passengers pass through.(10)cij=c10<j-i≤i1c2i1<j-i≤i2c3i2<j-i≤i3⋮⋮cnin-1<j-i≤inAssuming that the number of passengers arriving at stationi during the time period [0,t) is Ni(t), because the passenger flow of line 4, phase 1 in Chengdu Metro is small, and the stations are equipped with passenger volume control measures to avoid congestion, that means we can use Poisson distribution to describe the arrival of passengers. Therefore, it satisfies the Poisson distribution with a distribution intensity λi; we can express it by (11).As the headway between two adjacent trains isTk, the waiting time of most passengers would be less than Tk, and the number of passengers waiting for a specific approaching train at station i is defined as Ni(Tk). In addition, some of the passengers may choose to wait for the next train instead, if they find the first arriving train too much crowded. Therefore, we consider a passenger boarding rate κ that equals the ratio calculated by the number of boarding passengers divided by Ni(Tk). We let the number of passengers traveling from station i to station j(j>i) to be dij(Tk) and the number of passengers traveling from station i′(0<i′<i) to station i to be dii(Tk). The variation in the number of passengers Ai in the train at station i, as illustrated in Figure 2, can be calculated by (12).(11)Nit=λit(12)Ai=κNit-∑i′=1i-1di′iTkFigure 2
Schematic diagram forAi.If the energy consumption of a train in sectioni is Ei,i+1 and the total energy consumption of the train from the original station to the terminal station is Et, then the relationship between Et and Ei,i+1 can be expressed in (13). In this model, the power conversion rate is denoted by ξ, and the electricity rate is denoted by η.(13)Et=∑i=1n-1Ei,i+1The total number of passengers orientated from stationi in a day is Ni(T0). Passengers’ value of time is VOT, which can be represented by the local hourly average wage since the Metro is running in the same city every day. The average travel cost can be represented as follows:(14)Um=VOTTm+Cm+εmwhere Um is the generalized travel expenses of the Metro; Tm is the travel time, which is defined as the beginning of passenger arrive at station i to the end of departure (the passenger arrive at station j(j>i) and it is the time that the train arrives at the station j(j>i)), including the waiting time and running time; Cm is the cost (ticket fares) for passengers; εm is a set of random variables that related to the travel cost, where εm=ϕci, and ϕ∈-10%,-5% is the influence coefficient.


### 3.3. Objective Functions and Constraints

The proposed model aims to maximize the revenue while minimizing the average travel costs of passengers. We assume that fixed costs such as air conditioning and lighting are not considered. In addition, we consider that the expenditure is mainly borne by revenue and do not consider the government subsidy for the cost of the subway, though some lines are subsidized by the government, but, considering the generality, we assume that the operating revenue is only derived from ticket sales; the tickets sold to passengers can serve to represent the revenue function. Therefore, the income equals the fare income minus the sum of the operating costs and the train running energy consumption. As a result, the operating revenue can be represented by(15)max⁡Z=∑i=1n-1∑j=i+1ncij·dijTk-EtξηIn the same way, the average travel costs of passengers can be expressed as the average travel costs per trip that include the travel time and the ticket fare costs. We determine the travel time costs based on the waiting time costs and the transit time costs. Equation (16) represents the transit time costs, In (16), dij(Tk) is the number of passengers traveling from station i to station j(j>i). ∑i,jdijTk·∑ij-1li,i+1/v-i,i+1+si-s1 is running time of passengers from i station toj(j>i) station; for example, if i=1, j=3, ∑i,jdijTk·∑ij-1li,i+1/v-i,i+1+si-s1 is equal to d13(Tk)·(l1,2/v-1,2+l2,3/v-2,3+s2), the d13(Tk)·(l1,2/v-1,2+l2,3/v-2,3+s2) represents the running time of passengers dij(Tk) from i station toj(j>i) station. Therefore, VOT·∑i=1n-1∑j=i+1ndijTk·∑ij-1li,i+1/v-i,i+1+si-s1 represents the transit time costs.Equation (17) represents the waiting time costs, in (17), because the arrival of passengers conforms to the Poisson distribution, so ∫0Tkκλitdt express the total waiting time of all passengers at station i within [0,Tk]. Therefore, VOT·∑i=1n-1∫0Tkκλitdt is the total waiting time of all waiting passengers when the train is from starting station to terminal station.Considering the effects of the random variablesεm, since the ticket fare costs can be expressed by (18), in (18), cij represents the fare from i station to j(j>i) station (see (10)). ϕ is a set of random variables related to the travel cost (see (14)). (1+ϕ)cijdij(Tk) represents the total fare of passengers dij(Tk) from i station to j(j>i) station. ∑i=1n-1∑j=i+1n(1+ϕ)cijdij(Tk) represents the total fare of all passengers when the train is moving from starting station to terminal station.Therefore, the average travel costs of passengers per trip can be represented by (19).(16)Umt=VOT·∑i=1n-1∑j=i+1ndijTk·∑ij-1li,i+1v-i,i+1+si-s1(17)Umw=VOT·∑i=1n-1∫0Tkκλitdt(18)Umf=∑i=1n-1∑j=i+1n1+ϕcijdijTk(19)Um=VOT∑i=1n-1∑j=i+1ndijTk·∑ij-1li,i+1/v-i,i+1+si-s1+∑i=1n-1∫0Tkκλitdt∑i=1n-1κNiTk+∑i=1n-1∑j=i+1n1+ϕcijdijTk∑i=1n-1κNiTkTo guarantee the QoS and prevent the train from extreme cases, the proposed multiobjective optimization model should satisfy the following constraints.(1) When the number of running trains per day isN, T0 is the operation time of a train in a full day and the operating speed of the train in the sections should satisfy the following constraint:(20)N-1Tk+∑i=1nli,i+1v-i,i+1+∑i=1n-1si=T0(2) The carrying capacity of the trains should meet the following constraint as (21). In (21), Ai represents the train arriving at i station and the number of passengers getting on the train minus the number of passengers getting off the train. Dδmax is the maximum carrying capacity of train, so ∑ii′Ai≤Dδmax∀i′=(1,2,3,…,n-1) means the total number of passengers on the train can not exceed the maximum carrying capacity of the train at any station.(21)∑ii′Ai≤Dδmax∀i′=1,2,3,…,n-1(3) The minimum following distance and the following maximum speed between two adjacent trains should satisfy the relations defined in (22). In (22), according to the kinematics formula, we know that v22-v12=2as, v2 is the terminal velocity, v1 is the Initial velocity, a is the acceleration, and s is the distance. Assuming that the speed of the following train is vi,i+1max and the speed of the preceding train is zero, the maximum braking force of train is bmax, so the distance Lf between the adjacent trains should be satisfied to Lf≥(vi,i+1max)2/2bmax, ensuring that there is no conflict between the adjacent trains.(22)Lf≥vi,i+1max22bmax


## 3.1. Train Energy Consumption and Speed Analysis

Accelerating, coasting, and braking are the main phases of a train movement when performing running activities between successive track sections [23]. One can ignore small variations in train speed as the preventive maintenance of infrastructures and facilities keeps the operational conditions at the required level. With this simplifying assumption, train motion formula can be described as shown in (1). In this equation, we consider the basic line resistance, the track gauge, the maximum speed, and the signal systems, for instance, but we do not consider the curve and tunnel resistances. (1)a=dvdt=μffv-rvdtdx=1vwhere a is the acceleration of train operation, v is the speed of train operation, t is the time of train operation, f(v) is the maximum unit traction of the train, μf is the coefficient of the traction output ratio, and r(v) is the unit basic resistance of the train; it indicates the basic resistance of trains under unit weight, N/kN.Using this formula, when a train runs in a section, the relationship between energy consumption and the maximum speed can be derived as follows:(2)Ei,i+1′=∫0vi,i+1maxμffv·M0·vμffv-rvdvwhere vi,i+1max is the optimal maximum speed of the train in section i; it is not equal to the maximum speed limit of the train in section i; vi,i+1max≤ the maximum speed limit (the maximum speed limit of subway trains is equal to 80 km/h in China); M0 is the weight of the train. Equation (3) is another expression proposed by Gu et al., to calculate the energy consumption of a train in a section from the perspective of traction work [24].(3)Ei,i+1′′=∫0li,i+11μffv·M0·x·dxMore simplified formulations can be derived from (3) as follows: (4)Ei,i+1′′=∫0li,i+11μffv·M0·x·dx=∫0v-i,i+1t1μffv·M0·x·dx=∫0v-i,i+1t1μffv·M0·vt·v·dt=∫0v-i,i+1t1μffv·M0·vt·va·dv=∫0v-i,i+1t1μffv·M0·vt·vμffv-rv·dvthere is r(v)/f(v)≈0, because, in practice, r(v) is much smaller than f(v) and μf=1, when the train runs at the accelerating phase in a section [24]. Therefore, (2) and (4) can be reformulated into (5), given μffv-r(v)≈f(v). (5)Ei,i+1′=∫0vi,i+1maxvM0dvEi,i+1′′=∫0v-i,i+1t1vM0·vt·dvIn (4) and (5), li,i+11 is the running distance of the train at the accelerating phase in section i. We know that li,i+11=v-i,i+1t1, where v-i,i+1 is the average running speed of the train in section i; t is the quantity of the state that satisfies oΔt≈0. It refers to the time period when the speed of the train increases from v to v+oΔv. After some calculation and simplification, (5) can be further developed into (6)Ei,i+1′=∫0vi,i+1maxvM0dv=12vi,i+1max2M0Ei,i+1′′=∫0v-i,i+1t1vM0·vt·dv=13v-i,i+1t13M0tt∈0,ti,i+1Ei,i+1′=Ei,i+1′′⇒vi,i+1maxv-i,i+1=23v-i,i+1t13t=23li,i+11t·t1li,i+11=constant;t1=constant;t=oΔtOwing to the theoretical derivation of the mathematical formula, we assume the correlation coefficientε between vi,i+1max and v-i,i+1, to be (7)vi,i+1maxv-i,i+1=1ε=vi,i+1maxv-i,i+1=1ε.We approximate the value ofε using train movement data in track sections, including the optimal maximum speed, average operating speed, and the relationship between the optimal maximum speed and the average operating speed [23]. The instance is shown in Figure 1, in which y1 represents the average operating speed while y2 represents the maximum operating speed and the correlation value ε can be calculated by (8). Given this, (2) can be transformed into (9). Therefore, we get the relation between average speed and energy consumption, but this relation has the range of validity as follows.(1)
The object of the study must be urban rail transit system; the maximum speed of trains should not exceed 80km/h. And considering the comfort of passengers, the maximum train acceleration is restricted to 1 m/s2, and the maximum train deceleration is constrained as -1 m/s2.(2)
The distance between the two stations is very short (the distance between the two stations should be less than 3 kilometers in principle); the train running in the section only includes three phases: accelerating, coasting, and braking, without cruising phase.(8)ε=113×∑x=113y1xy2x(9)Ei,i+1=∫0v-i,i+1/εμffv·M0·vμffv-rv·dvFigure 1
The trend of polyliney1 and polyline y2.


## 3.2. Description of Variables and Model Assumptions

We first provide the definitions and assumptions used in this study and then describe ticket pricing schemes in URT. We denote the number of stations along a line byn, the headway between successive trains in the same direction by Tk, the operation time of a train in a full day by T0, and the running time of a train in section i by ti,i+1. Also, we set si as the dwell time of a train at station i, the passenger carrying capacity of a train as D, the occupied rate of passenger carrying capacity as δ, the total length of the Metro line as L, and, finally, the length of section i as li,i+1. Generally, two different kinds of ticket pricing schemes are used in URT. The first one depends on the train operating distance, and the second one depends on the number of pass-through sections. In Chengdu Metro, the latter method is applied for which the price of tickets can be expressed by (10). In this equation (c1,c2,…,cn) indicates different ticket prices, where j-i to the number of sections that passengers pass through.(10)cij=c10<j-i≤i1c2i1<j-i≤i2c3i2<j-i≤i3⋮⋮cnin-1<j-i≤inAssuming that the number of passengers arriving at stationi during the time period [0,t) is Ni(t), because the passenger flow of line 4, phase 1 in Chengdu Metro is small, and the stations are equipped with passenger volume control measures to avoid congestion, that means we can use Poisson distribution to describe the arrival of passengers. Therefore, it satisfies the Poisson distribution with a distribution intensity λi; we can express it by (11).As the headway between two adjacent trains isTk, the waiting time of most passengers would be less than Tk, and the number of passengers waiting for a specific approaching train at station i is defined as Ni(Tk). In addition, some of the passengers may choose to wait for the next train instead, if they find the first arriving train too much crowded. Therefore, we consider a passenger boarding rate κ that equals the ratio calculated by the number of boarding passengers divided by Ni(Tk). We let the number of passengers traveling from station i to station j(j>i) to be dij(Tk) and the number of passengers traveling from station i′(0<i′<i) to station i to be dii(Tk). The variation in the number of passengers Ai in the train at station i, as illustrated in Figure 2, can be calculated by (12).(11)Nit=λit(12)Ai=κNit-∑i′=1i-1di′iTkFigure 2
Schematic diagram forAi.If the energy consumption of a train in sectioni is Ei,i+1 and the total energy consumption of the train from the original station to the terminal station is Et, then the relationship between Et and Ei,i+1 can be expressed in (13). In this model, the power conversion rate is denoted by ξ, and the electricity rate is denoted by η.(13)Et=∑i=1n-1Ei,i+1The total number of passengers orientated from stationi in a day is Ni(T0). Passengers’ value of time is VOT, which can be represented by the local hourly average wage since the Metro is running in the same city every day. The average travel cost can be represented as follows:(14)Um=VOTTm+Cm+εmwhere Um is the generalized travel expenses of the Metro; Tm is the travel time, which is defined as the beginning of passenger arrive at station i to the end of departure (the passenger arrive at station j(j>i) and it is the time that the train arrives at the station j(j>i)), including the waiting time and running time; Cm is the cost (ticket fares) for passengers; εm is a set of random variables that related to the travel cost, where εm=ϕci, and ϕ∈-10%,-5% is the influence coefficient.


## 3.3. Objective Functions and Constraints

The proposed model aims to maximize the revenue while minimizing the average travel costs of passengers. We assume that fixed costs such as air conditioning and lighting are not considered. In addition, we consider that the expenditure is mainly borne by revenue and do not consider the government subsidy for the cost of the subway, though some lines are subsidized by the government, but, considering the generality, we assume that the operating revenue is only derived from ticket sales; the tickets sold to passengers can serve to represent the revenue function. Therefore, the income equals the fare income minus the sum of the operating costs and the train running energy consumption. As a result, the operating revenue can be represented by(15)max⁡Z=∑i=1n-1∑j=i+1ncij·dijTk-EtξηIn the same way, the average travel costs of passengers can be expressed as the average travel costs per trip that include the travel time and the ticket fare costs. We determine the travel time costs based on the waiting time costs and the transit time costs. Equation (16) represents the transit time costs, In (16), dij(Tk) is the number of passengers traveling from station i to station j(j>i). ∑i,jdijTk·∑ij-1li,i+1/v-i,i+1+si-s1 is running time of passengers from i station toj(j>i) station; for example, if i=1, j=3, ∑i,jdijTk·∑ij-1li,i+1/v-i,i+1+si-s1 is equal to d13(Tk)·(l1,2/v-1,2+l2,3/v-2,3+s2), the d13(Tk)·(l1,2/v-1,2+l2,3/v-2,3+s2) represents the running time of passengers dij(Tk) from i station toj(j>i) station. Therefore, VOT·∑i=1n-1∑j=i+1ndijTk·∑ij-1li,i+1/v-i,i+1+si-s1 represents the transit time costs.Equation (17) represents the waiting time costs, in (17), because the arrival of passengers conforms to the Poisson distribution, so ∫0Tkκλitdt express the total waiting time of all passengers at station i within [0,Tk]. Therefore, VOT·∑i=1n-1∫0Tkκλitdt is the total waiting time of all waiting passengers when the train is from starting station to terminal station.Considering the effects of the random variablesεm, since the ticket fare costs can be expressed by (18), in (18), cij represents the fare from i station to j(j>i) station (see (10)). ϕ is a set of random variables related to the travel cost (see (14)). (1+ϕ)cijdij(Tk) represents the total fare of passengers dij(Tk) from i station to j(j>i) station. ∑i=1n-1∑j=i+1n(1+ϕ)cijdij(Tk) represents the total fare of all passengers when the train is moving from starting station to terminal station.Therefore, the average travel costs of passengers per trip can be represented by (19).(16)Umt=VOT·∑i=1n-1∑j=i+1ndijTk·∑ij-1li,i+1v-i,i+1+si-s1(17)Umw=VOT·∑i=1n-1∫0Tkκλitdt(18)Umf=∑i=1n-1∑j=i+1n1+ϕcijdijTk(19)Um=VOT∑i=1n-1∑j=i+1ndijTk·∑ij-1li,i+1/v-i,i+1+si-s1+∑i=1n-1∫0Tkκλitdt∑i=1n-1κNiTk+∑i=1n-1∑j=i+1n1+ϕcijdijTk∑i=1n-1κNiTkTo guarantee the QoS and prevent the train from extreme cases, the proposed multiobjective optimization model should satisfy the following constraints.(1) When the number of running trains per day isN, T0 is the operation time of a train in a full day and the operating speed of the train in the sections should satisfy the following constraint:(20)N-1Tk+∑i=1nli,i+1v-i,i+1+∑i=1n-1si=T0(2) The carrying capacity of the trains should meet the following constraint as (21). In (21), Ai represents the train arriving at i station and the number of passengers getting on the train minus the number of passengers getting off the train. Dδmax is the maximum carrying capacity of train, so ∑ii′Ai≤Dδmax∀i′=(1,2,3,…,n-1) means the total number of passengers on the train can not exceed the maximum carrying capacity of the train at any station.(21)∑ii′Ai≤Dδmax∀i′=1,2,3,…,n-1(3) The minimum following distance and the following maximum speed between two adjacent trains should satisfy the relations defined in (22). In (22), according to the kinematics formula, we know that v22-v12=2as, v2 is the terminal velocity, v1 is the Initial velocity, a is the acceleration, and s is the distance. Assuming that the speed of the following train is vi,i+1max and the speed of the preceding train is zero, the maximum braking force of train is bmax, so the distance Lf between the adjacent trains should be satisfied to Lf≥(vi,i+1max)2/2bmax, ensuring that there is no conflict between the adjacent trains.(22)Lf≥vi,i+1max22bmax


## 4. Solution Approaches

### 4.1. Fuzzy Multiobjective Optimization Algorithm

A fuzzy multiobjective optimization algorithm is used to optimize trains’ operating speed. To reflect the preferences of decision-makers regarding the objective functions, we have used different weight coefficients to obtain a trade-off between the expectations of decision-makers. The steps of the algorithm are listed as follows.Step 1.
Transform the objective functionZ into a standard form, and let X=-Z, so that the objective function Z can be converted to a minimization objective function X.Step 2.
Construct minimized ideal value vectorGmin and maximized inverse ideal value vector Gmax of two objective functions X and Um while satisfying all of the constraints, as shown in (23)Gmin=Xmin,UmminGmax=Xmax,Ummaxwhere Xmin and Ummin are the minimum value of the objective functions X and Um, respectively, while Xmax and Ummax are the maximum value of the objective functions X and Um.Step 3.
Build membership functionsψ1 and ψ2 of two objective functions X and Um, as shown in the following.(24)ψ1=1X≤XminXmax-XXmax-XminXmin<X<Xmax0X≥Xmax(25)ψ2=1U≤UmminUmmax-UUmmax-UmminUmmin<U<Ummax0U≥UmmaxStep 4.
According to the principle of the fuzzy multiobjective optimization algorithm, the multiobjective optimization problem is transformed into the following single-objective optimization problem, as shown in(26)L=max⁡λ1ψ1p+λ2ψ2p1/pwhere L is the new objective function combining X and Um; λi(i=1,2)≥0 with ∑i=12λi=1 is the weight coefficient of the original objective functions X and Um. In (26), p is the distance parameter.Step 5.
Change the value ofλ according to their preferences; change the weight coefficients of the original objective functions X and Um; select different distance models by changing the value of p; the decision makers can specify the value of p according to different preferences which are corresponding to different distance calculation methods, as shown in Table 1.Table 1
Value ofp corresponding to different distance calculation methods.
Value ofp p = 1 p = 2 p = + ∞ The Meaning of Objective Function Manhattan distance Euclidean distance Chebyshev distance Objective function model L = max ⁡ ( λ 1 ψ 1 + λ 2 ψ 2 ) L = max ⁡ ( λ 1 ψ 1 ) 2 + ( λ 2 ψ 2 ) 2 L = max ⁡ max ⁡ λ 1 ψ 1 , λ 2 ψ 2


### 4.2. Nondominated Sorting Genetic Algorithm II (NSGA-II)

The NSGA-II algorithm is improved based on the original NSGA algorithm, and a fast nondominated sorting method is proposed to cope with the complexity issue [25]. The NSGA algorithm uses the congestion comparison operator and does require specifying the shared parameters. Also, the introduction of elite strategy and the expanding of the sampling space allows the parents and their offspring participate in the competition to produce the next generation of the population to generate better offspring.(1) Dominance and Noninferior. In the multiobjective optimization problem, if all the targets of the individualp are not worse than individualq and there is at least one target of the individualp which is better than that of the individualq, then we say thatp dominatesq, and it also implies thatp is noninferior toq.(2) Rank and Front. Ifp dominatesq, then the order value ofp is lower thanq; otherwise, ifp andq do not dominate each other, or ifp andq are noninferior to each other, thenp andq have the same order value. An individual with the order value of 1 belongs to the first front, while the individual with the order value of 2 belongs to the second front. With the current population, the first front is completely unobstructed, and the second front is dominated by the individuals which are in the first front. In this way, all the individuals in the population are assigned to different fronts.(3) Crowding Distance. Through step (2), we know that ifp andq are noninferior to each other, thenp andq have the same order value; it meansp andq in the same front. The crowding distance is used to calculate the distance between the individuals (p andq) in the same front, and it can be used to characterize the degree of crowding between individuals. The greater the crowded distance is, the less crowded is, and the better the diversity of the population is. As the crowded distance is used to calculate the distance between the individuals (p andq) in the same front, the greater the crowded distance is, the smaller the number of individuals distributed at the same front is and the greater the total fronts are because the population size is fixed. In this case, we call the better the diversity of the population is. For example, we assume that the population size is 500; case 1 is as follows: all individuals are distributed at three fronts; case 2 is as follows: all individuals are distributed at four fronts, so we can say that the population diversity of case 2 is better than that of case 1.(4) Optimal Front Individual ParetoFraction. The optimal front individual ParetoFraction is defined as the proportion (ranges from 0 to 1) of the individuals in the optimal front (Pareto front) in the population. In other words, it is the Pareto front that is equal to the minimal value between ParetoFraction multiplied by the population size and the number of existing individuals in the front. For example, we assume that the population size is 500 and the generations is 100, because the optimal front individual ParetoFraction is defined as the proportion (ranges from 0 to 1); we assume it equals 0.3, so the optimal front equals 150. But considering that the population has to iterate to generate new populations, after 100 iterations, the number of existing individuals in the front is 100. According to the rules, we know that the Pareto front = min⁡(150,100), so the Pareto front =100.In summary, the framework of the proposed NSGA-II algorithm is shown in Figure3, in which “Gen” is the counter of the number of generations, population consolidation is that selection, crossover, and mutation of the previous population; we can get the new population, through the nondominant relationship and the crowding degree of individuals; we can get suitable individuals of the new population, and the suitable individuals are selected to form a new paternal population and then continue to produce new offspring population. Congestion calculation is the calculation of the density of a given individual and its surrounding individuals in a population.Figure 3
The framework of the NSGA-II algorithm.


## 4.1. Fuzzy Multiobjective Optimization Algorithm

A fuzzy multiobjective optimization algorithm is used to optimize trains’ operating speed. To reflect the preferences of decision-makers regarding the objective functions, we have used different weight coefficients to obtain a trade-off between the expectations of decision-makers. The steps of the algorithm are listed as follows.Step 1.
Transform the objective functionZ into a standard form, and let X=-Z, so that the objective function Z can be converted to a minimization objective function X.Step 2.
Construct minimized ideal value vectorGmin and maximized inverse ideal value vector Gmax of two objective functions X and Um while satisfying all of the constraints, as shown in (23)Gmin=Xmin,UmminGmax=Xmax,Ummaxwhere Xmin and Ummin are the minimum value of the objective functions X and Um, respectively, while Xmax and Ummax are the maximum value of the objective functions X and Um.Step 3.
Build membership functionsψ1 and ψ2 of two objective functions X and Um, as shown in the following.(24)ψ1=1X≤XminXmax-XXmax-XminXmin<X<Xmax0X≥Xmax(25)ψ2=1U≤UmminUmmax-UUmmax-UmminUmmin<U<Ummax0U≥UmmaxStep 4.
According to the principle of the fuzzy multiobjective optimization algorithm, the multiobjective optimization problem is transformed into the following single-objective optimization problem, as shown in(26)L=max⁡λ1ψ1p+λ2ψ2p1/pwhere L is the new objective function combining X and Um; λi(i=1,2)≥0 with ∑i=12λi=1 is the weight coefficient of the original objective functions X and Um. In (26), p is the distance parameter.Step 5.
Change the value ofλ according to their preferences; change the weight coefficients of the original objective functions X and Um; select different distance models by changing the value of p; the decision makers can specify the value of p according to different preferences which are corresponding to different distance calculation methods, as shown in Table 1.Table 1
Value ofp corresponding to different distance calculation methods.
Value ofp p = 1 p = 2 p = + ∞ The Meaning of Objective Function Manhattan distance Euclidean distance Chebyshev distance Objective function model L = max ⁡ ( λ 1 ψ 1 + λ 2 ψ 2 ) L = max ⁡ ( λ 1 ψ 1 ) 2 + ( λ 2 ψ 2 ) 2 L = max ⁡ max ⁡ λ 1 ψ 1 , λ 2 ψ 2


## 4.2. Nondominated Sorting Genetic Algorithm II (NSGA-II)

The NSGA-II algorithm is improved based on the original NSGA algorithm, and a fast nondominated sorting method is proposed to cope with the complexity issue [25]. The NSGA algorithm uses the congestion comparison operator and does require specifying the shared parameters. Also, the introduction of elite strategy and the expanding of the sampling space allows the parents and their offspring participate in the competition to produce the next generation of the population to generate better offspring.(1) Dominance and Noninferior. In the multiobjective optimization problem, if all the targets of the individualp are not worse than individualq and there is at least one target of the individualp which is better than that of the individualq, then we say thatp dominatesq, and it also implies thatp is noninferior toq.(2) Rank and Front. Ifp dominatesq, then the order value ofp is lower thanq; otherwise, ifp andq do not dominate each other, or ifp andq are noninferior to each other, thenp andq have the same order value. An individual with the order value of 1 belongs to the first front, while the individual with the order value of 2 belongs to the second front. With the current population, the first front is completely unobstructed, and the second front is dominated by the individuals which are in the first front. In this way, all the individuals in the population are assigned to different fronts.(3) Crowding Distance. Through step (2), we know that ifp andq are noninferior to each other, thenp andq have the same order value; it meansp andq in the same front. The crowding distance is used to calculate the distance between the individuals (p andq) in the same front, and it can be used to characterize the degree of crowding between individuals. The greater the crowded distance is, the less crowded is, and the better the diversity of the population is. As the crowded distance is used to calculate the distance between the individuals (p andq) in the same front, the greater the crowded distance is, the smaller the number of individuals distributed at the same front is and the greater the total fronts are because the population size is fixed. In this case, we call the better the diversity of the population is. For example, we assume that the population size is 500; case 1 is as follows: all individuals are distributed at three fronts; case 2 is as follows: all individuals are distributed at four fronts, so we can say that the population diversity of case 2 is better than that of case 1.(4) Optimal Front Individual ParetoFraction. The optimal front individual ParetoFraction is defined as the proportion (ranges from 0 to 1) of the individuals in the optimal front (Pareto front) in the population. In other words, it is the Pareto front that is equal to the minimal value between ParetoFraction multiplied by the population size and the number of existing individuals in the front. For example, we assume that the population size is 500 and the generations is 100, because the optimal front individual ParetoFraction is defined as the proportion (ranges from 0 to 1); we assume it equals 0.3, so the optimal front equals 150. But considering that the population has to iterate to generate new populations, after 100 iterations, the number of existing individuals in the front is 100. According to the rules, we know that the Pareto front = min⁡(150,100), so the Pareto front =100.In summary, the framework of the proposed NSGA-II algorithm is shown in Figure3, in which “Gen” is the counter of the number of generations, population consolidation is that selection, crossover, and mutation of the previous population; we can get the new population, through the nondominant relationship and the crowding degree of individuals; we can get suitable individuals of the new population, and the suitable individuals are selected to form a new paternal population and then continue to produce new offspring population. Congestion calculation is the calculation of the density of a given individual and its surrounding individuals in a population.Figure 3
The framework of the NSGA-II algorithm.


## 5. Numerical Experiments

### 5.1. Numerical Experiments Setup

For the numerical experiment, the URT line 4 of phase 1 in Chengdu Metro is studied to show the effect of using optimal train operating speed obtained by the proposed model. It has 16 stations, which are numbered from M1 to M16; the basic information of this line, including the length of the sections, the travel time, and dwell times, respectively, for each section and at each station, is provided in Table2. The operation hour of the trains is from 6:30 to 22:30, the headway is 3 minutes in peak hours and it is 5 minutes in nonpeak hours; the maximum operating speed is 80km/h. Considering the comfort of passengers, the maximum train acceleration is restricted to 1 m/s2 and the maximum train deceleration is constrained as -1 m/s2. Finally, both the headway between two adjacent trains and the turn-back-time are set as 300s, and other parameters are listed in Table 3.Table 2
Basic information of line 4, phase 1 in Chengdu Metro.
Station ID M1 M2 M3 M4 M5 M6 M7 M8 Dwelling Time/s 25 25 25 25 25 30 25 25 Section Length/m 2585 2140 1887 2150 1171 985 900 1166 Operation Time/s 154 138 108 139 74 72 67 84  Station ID M9 M10 M11 M12 M13 M14 M15 M16  Dwelling Time/s 30 30 30 30 30 30 25 25 Section Length/m 972 1537 1140 750 1630 1560 820 Operation Time/s 71 113 83 61 105 101 64Table 3
Other parameters used in the model.
κ ϕ V O T/(CNY/h) M 0 /kg ξ η / C N Y /(kw⋅h) δ m a x 0.98 -10% 25.63 196000 0.95 1 1.2Also, the empirical train traction characteristic curve and running characteristic resistance curve corresponding to (27), (28), and (29) are shown in Figure 4 [26].(27)fmax=2030≤v≤51.5km/h-0.002032v3+0.4928v2-42.13v+134351.5<v≤80km/h(28)bmax=1660≤v≤77km/h0.1343v2-25.07v+130077<v≤80km/h(29)w0=2.031+0.0622v+0.001807v2Figure 4
The traction characteristic curve and running characteristic resistance curve.The daily passenger flow of the metro line is depicted in Figure5. The node of the polyline “passenger flow of M1” represents the passenger flow from M1 to any other succeeding stations, and the other 14 polylines have the same functions with the polyline “passenger flow of M1”.Figure 5
The daily passenger flow of line 4, phase 1 of Chengdu Metro (Source: The ACC office of Chengdu Metro).Finally, given the number of sections that the passengers pass through, the price of tickets in Chengdu Metro is listed as follows:(30)cij=20<j-i≤636<j-i≤10410<j-i≤16516<j-i≤24


### 5.2. Results of the Fuzzy Multi-Objective Optimization Algorithm

We use MATLAB to implement the proposed optimization algorithm to obtain the optimal train operating speed in every section under different input weight and distance parameters, as shown in Table4. The optimization results show that when the weight of the operating revenue is as important as the average travel cost of the passengers in the objective function (i.e., λ1=λ2=0.5), the optimal speed (v-) is 18.3m/s. In that case, the maximum revenue of the Metro company (Z) is 491.50 CNY per train, and the average travel costs of the passengers (Um) are 9.13 CNY with the distance parameter (p) being equal to 1 or 2.Table 4
The objective function values and the respective optimal train operating speed obtained by the fuzzy multiobjective optimization algorithm.
Value ofp Value ofλ1,λ2 v - (m/s) Z (CNY/train) U m (CNY/train) p = 1 (0.0,1.0) 18.3 491.50 9.13 (0.2,0.8) 18.3 491.50 9.13 (0.4,0.6) 18.3 491.50 9.13 (0.5,0.5) 18.3 491.50 9.13 (0.6,0.4) 18.3 491.50 9.13 (0.8,0.2) 18.3 491.50 9.13 (1.0,0.0) 11.5 578.60 10.77  p = 2 (0.0,1.0) 18.3 491.50 9.13 (0.2,0.8) 18.3 491.50 9.13 (0.4,0.6) 18.3 491.50 9.13 (0.5,0.5) 18.3 491.50 9.13 (0.6,0.4) 18.3 491.50 9.13 (0.8,0.2) 11.5 578.60 10.77 (1.0,0.0) 11.5 578.60 10.77  p = + ∞ (0.0,1.0) 18.3 491.50 9.13 (0.2,0.8) 11.5 578.60 10.77 (0.4,0.6) 11.5 578.60 10.77 (0.5,0.5) 11.5 578.60 10.77 (0.6,0.4) 11.5 578.60 10.77 (0.8,0.2) 11.5 578.60 10.77 (1.0,0.0) 11.5 578.60 10.77 Note: because the optimization result shows thatv-1,2=v-2,3=v-i,i+1…=v-15,16, then we use v- to express average train operating speed in every section.When the operating revenue of the Metro company is more important than the average travel cost of the passengers in the objective function (i.e.,λ1=0.8, λ2=0.2; λ1=1.0, λ2=0), the optimal train operating speed (v-) is approximately 11.5m/s. In that case, the maximum revenue of the Metro company (Z) is about 578.60 CNY per train, and the average travel costs of the passengers (Um) are almost 10.77 CNY. These results imply that the maximum operating revenue is positively related to the optimal speed while the average travel costs of the passengers are negatively correlated with the optimal train operating speed.


### 5.3. Results of the Nondominated Sorting Genetic Algorithm II

We use the “gamultiobj” built-in function in MATLAB to solve our biobjective model with the NSGA-II algorithm. For this case, our parameters are tuned as follows: ParetoFraction = 0.3, population size = 100, generations = 200, “stall Gen Limit” = 200, and “Tol Fun” = 0.01 in the options of function gaoptimset. Here, “ParetoFraction” is the optimal individual coefficient, “stall Gen Limit” is the generation of stopping iteration, and “Tol Fun” is the error of fitness function. Figure6 depicts the Pareto-frontier solutions, in which objective 1 means -Z and objective 2 represents Um. It presents the Pareto solutions based on NSGA-II algorithm, showing the 100 sets of data of v-, min⁡-Z, and min⁡Um, including two extreme optimization results. It can be concluded that (v-=11.50,Z=578.60,Um=10.77) is the best optimal solution from the viewpoint of the company while (v-=18.30,Z=491.54,Um=9.13) is the best optimal solution from the viewpoint of the passengers. The specific optimization results based on NSGA-II algorithm are shown in Table 7.Figure 6
The Pareto solutions obtained by the NSGA-II algorithm.


### 5.4. A Comparative Analysis of the Optimization Algorithms

The fuzzy multiobjective optimization algorithm is to obtain the optimal operating speed by changing the weight and distance parameters in the objective function, so the subjective intention of this algorithm is very significant, and the result will be affected easily by the coefficient values of objective functions, as shown in Table4. Because of the significant differences in the values of those two objective functions, the optimization results thus do not change after increasing the weight of the first objective function.The relationship between average speed and energy consumption is as shown in Figure7, which indicates that operation energy consumption is increasing with the increasing of average speed.Figure 7
The relationship between the average speed and energy consumption.The optimization results of the NSGA-II algorithm, as shown in Figure8, indicate that the higher operating speed benefits the passengers more, and lower operating speed benefits the metro operating company more.Figure 8
The trade-off between the objective functions and the train average operating speed in the sections.


### 5.5. Comprehensive Analysis of the Optimized Timetable

The duration of the running cycle in the current timetable is 3,938 seconds, and the absolute error of total travel time variations between the current timetable and the optimized timetable is restricted to ±120 seconds. The travel time of the current timetable (CUT), the optimized timetable with the operating speed benefitting the company (OTC), and the optimized timetable with the operating speed benefitting the passengers (OTP) respectively are presented in Figure9. Tables 5 and 6 present the timetable of OTC and OTP.Table 5
The optimized timetable of OTC.
Down M1 M2 M3 M4 M5 M6 M7 M8  Dwell (s) - 25 25 25 25 30 25 25  Arrival(s) - 161 330 469 639 741 846 941  Departure(s) 0 186 355 494 664 771 871 966  M9 M10 M11 M12 M13 M14 M15 M16  Dwell (s) 30 30 30 30 30 30 25 -  Arrival(s) 1053 1158 1304 1420 1513 1652 1788 1879  Departure(s) 1083 1188 1334 1450 1543 1682 1813 -  UP M16 M15 M14 M13 M12 M11 M10 M9  Dwell time(s) - 25 30 30 30 30 30 30  Arrival(s) - 2245 2376 2515 2608 2724 2870 2975  Departure(s) 2179 2270 2406 2545 2638 2754 2900 3005  M8 M7 M6 M5 M4 M3 M2 M1  Dwell (s) 25 25 30 25 25 25 25 -  Arrival(s) 3092 3187 3287 3394 3564 3703 3872 4058  Departure(s) 3117 3212 3317 3419 3589 3728 3897 -Table 6
The optimized timetable of OTP.
Down M1 M2 M3 M4 M5 M6 M7 M8  Dwell (s) - 25 25 25 25 30 25 25  Arrival(s) - 147 304 431 589 685 784 873  Departure(s) 0 172 329 456 614 715 809 898  M9 M10 M11 M12 M13 M14 M15 M16  Dwell (s) 30 30 30 30 30 30 25 -  Arrival(s) 979 1076 1216 1326 1415 1546 1672 1759  Departure(s) 1009 1106 1246 1356 1445 1576 1697 -  UP M16 M15 M14 M13 M12 M11 M10 M9  Dwell time(s) - 25 30 30 30 30 30 30  Arrival(s) - 2121 2242 2373 2462 2572 2712 2809  Departure(s) 2059 2146 2272 2403 2492 2602 2742 2839  M8 M7 M6 M5 M4 M3 M2 M1  Dwell (s) 25 25 30 25 25 25 25 -  Arrival(s) 2920 3009 3103 3204 3362 3489 3646 3818  Departure(s) 2945 3034 3133 3229 3387 3514 3671 -Table 7
The optimization results based on NSGA-II algorithm.
v - min ⁡ ( - Z ) min ⁡ U m v - min ⁡ ( - Z ) min ⁡ U m v - m i n ( - Z ) min ⁡ U m v - min ⁡ ( - Z ) min ⁡ U m 11.50 -578.60 10.77 16.67 -516.00 9.40 17.29 -506.88 9.29 15.76 -528.74 9.58 18.30 -491.54 9.13 13.42 -558.03 10.14 11.73 -576.25 10.68 13.53 -556.77 10.11 18.30 -491.54 9.13 15.45 -532.82 9.64 14.03 -550.82 9.98 12.48 -568.53 10.43 12.70 -566.12 10.36 11.78 -575.84 10.67 16.29 -521.35 9.47 16.11 -523.93 9.51 14.02 -550.90 9.98 17.85 -498.48 9.20 14.53 -544.73 9.85 13.80 -553.56 10.04 16.95 -511.95 9.35 18.22 -492.76 9.14 15.24 -535.67 9.69 14.85 -540.71 9.78 11.50 -578.60 10.77 15.09 -537.54 9.72 12.09 -572.61 10.56 16.22 -522.29 9.49 11.86 -575.00 10.64 16.40 -519.83 9.45 12.30 -570.44 10.49 16.92 -512.42 9.36 15.30 -534.80 9.68 12.19 -571.59 10.52 12.35 -569.87 10.47 13.96 -551.73 10.00 18.09 -494.74 9.16 17.22 -507.98 9.31 15.16 -536.70 9.71 17.60 -502.37 9.24 16.33 -520.82 9.47 12.81 -564.87 10.32 14.10 -549.96 9.96 15.58 -531.07 9.62 12.02 -573.29 10.58 12.51 -568.17 10.42 11.98 -573.79 10.60 16.06 -524.62 9.52 13.90 -552.45 10.01 17.77 -499.67 9.21 14.60 -543.87 9.84 17.53 -503.37 9.25 14.26 -548.05 9.92 11.63 -577.30 10.72 12.86 -564.40 10.31 13.84 -553.12 10.03 15.71 -529.34 9.59 16.80 -514.11 9.38 17.45 -504.56 9.27 15.39 -533.69 9.66 17.34 -506.26 9.29 17.96 -496.86 9.18 13.58 -556.17 10.10 13.65 -555.37 10.08 14.81 -541.15 9.79 15.80 -528.17 9.57 14.75 -541.91 9.80 11.89 -574.69 10.63 16.88 -512.92 9.36 14.92 -539.80 9.76 11.55 -578.15 10.75 13.24 -560.09 10.19 12.73 -565.75 10.34 16.58 -517.35 9.42 15.49 -532.25 9.63 13.01 -562.66 10.26 14.70 -542.51 9.81 14.43 -545.89 9.87 11.67 -576.93 10.71 16.80 -514.11 9.38 17.07 -510.20 9.33 14.20 -548.74 9.93 14.33 -547.20 9.90 15.03 -538.32 9.73 13.06 -562.17 10.25 16.48 -518.74 9.44 12.40 -569.32 10.45 18.29 -491.66 9.13 11.59 -577.75 10.74 17.17 -508.79 9.32 13.15 -561.09 10.22 15.16 -536.60 9.71 14.33 -547.18 9.90 17.66 -501.37 9.23 13.75 -554.23 10.05 17.08 -510.08 9.33 12.78 -565.24 10.33 12.26 -570.86 10.50 12.12 -572.28 10.55 14.47 -545.49 9.87Figure 9
Comparison between the optimized timetable and the current timetable.The optimization results indicate that the operating revenue (Z) in CUT, OTC, and OTP is 536.52 CNY, 544.32 CNY, and 527.73 CNY, respectively. The average travel costs of the passengers (Um) in CUT, OTC, and OTP are 9.70 CNY, 9.84 CNY, and 9.56 CNY, respectively. Comparing with the OTC and CUT, the revenue is increased by 1.45%, and the average energy consumption of each train in sections is reduced by 7.88%. Considering maximizing operating revenue, we conclude that the OTC performs better than the current one CUT and OTP.


### 5.6. The Sensitivity of the Main Parameters

In order to research the relationship between the main parameters and decision variablesv-i,i+1 in the model, we analyzed the sensitivity of the main parameters (the electricity rate η and passengers’ value of time VOT) affecting train running costs and passengers’ average travel cost in the model. Assuming the sensitivity range of η is 0.5-1.5 CNY/(kw·h), the step length is 0.1 CNY/(kw·h), and the sensitivity range of VOT is 18-30 CNY/h, the step length is 1 CNY/h. Therefore, we get the relationship between the two parameters and decision variables v-i,i+1 under different weights at p=1, as shown in Figure 10.Figure 10
The trade-off between the main parameters and the train average operating speed in the sections.The train average operating speed in the sections decreases with the increase of the electricity rateη and decreases with the increase of λ1 when the electricity rate η is greater than 1.3 CNY/(kw·h). Also, the train average operating speed in the sections increases with the increase of passengers’ value of time VOT and increases with the decrease of λ1 when passengers’ value of time VOT is less than 23 CNY/h.


## 5.1. Numerical Experiments Setup

For the numerical experiment, the URT line 4 of phase 1 in Chengdu Metro is studied to show the effect of using optimal train operating speed obtained by the proposed model. It has 16 stations, which are numbered from M1 to M16; the basic information of this line, including the length of the sections, the travel time, and dwell times, respectively, for each section and at each station, is provided in Table2. The operation hour of the trains is from 6:30 to 22:30, the headway is 3 minutes in peak hours and it is 5 minutes in nonpeak hours; the maximum operating speed is 80km/h. Considering the comfort of passengers, the maximum train acceleration is restricted to 1 m/s2 and the maximum train deceleration is constrained as -1 m/s2. Finally, both the headway between two adjacent trains and the turn-back-time are set as 300s, and other parameters are listed in Table 3.Table 2
Basic information of line 4, phase 1 in Chengdu Metro.
Station ID M1 M2 M3 M4 M5 M6 M7 M8 Dwelling Time/s 25 25 25 25 25 30 25 25 Section Length/m 2585 2140 1887 2150 1171 985 900 1166 Operation Time/s 154 138 108 139 74 72 67 84  Station ID M9 M10 M11 M12 M13 M14 M15 M16  Dwelling Time/s 30 30 30 30 30 30 25 25 Section Length/m 972 1537 1140 750 1630 1560 820 Operation Time/s 71 113 83 61 105 101 64Table 3
Other parameters used in the model.
κ ϕ V O T/(CNY/h) M 0 /kg ξ η / C N Y /(kw⋅h) δ m a x 0.98 -10% 25.63 196000 0.95 1 1.2Also, the empirical train traction characteristic curve and running characteristic resistance curve corresponding to (27), (28), and (29) are shown in Figure 4 [26].(27)fmax=2030≤v≤51.5km/h-0.002032v3+0.4928v2-42.13v+134351.5<v≤80km/h(28)bmax=1660≤v≤77km/h0.1343v2-25.07v+130077<v≤80km/h(29)w0=2.031+0.0622v+0.001807v2Figure 4
The traction characteristic curve and running characteristic resistance curve.The daily passenger flow of the metro line is depicted in Figure5. The node of the polyline “passenger flow of M1” represents the passenger flow from M1 to any other succeeding stations, and the other 14 polylines have the same functions with the polyline “passenger flow of M1”.Figure 5
The daily passenger flow of line 4, phase 1 of Chengdu Metro (Source: The ACC office of Chengdu Metro).Finally, given the number of sections that the passengers pass through, the price of tickets in Chengdu Metro is listed as follows:(30)cij=20<j-i≤636<j-i≤10410<j-i≤16516<j-i≤24


## 5.2. Results of the Fuzzy Multi-Objective Optimization Algorithm

We use MATLAB to implement the proposed optimization algorithm to obtain the optimal train operating speed in every section under different input weight and distance parameters, as shown in Table4. The optimization results show that when the weight of the operating revenue is as important as the average travel cost of the passengers in the objective function (i.e., λ1=λ2=0.5), the optimal speed (v-) is 18.3m/s. In that case, the maximum revenue of the Metro company (Z) is 491.50 CNY per train, and the average travel costs of the passengers (Um) are 9.13 CNY with the distance parameter (p) being equal to 1 or 2.Table 4
The objective function values and the respective optimal train operating speed obtained by the fuzzy multiobjective optimization algorithm.
Value ofp Value ofλ1,λ2 v - (m/s) Z (CNY/train) U m (CNY/train) p = 1 (0.0,1.0) 18.3 491.50 9.13 (0.2,0.8) 18.3 491.50 9.13 (0.4,0.6) 18.3 491.50 9.13 (0.5,0.5) 18.3 491.50 9.13 (0.6,0.4) 18.3 491.50 9.13 (0.8,0.2) 18.3 491.50 9.13 (1.0,0.0) 11.5 578.60 10.77  p = 2 (0.0,1.0) 18.3 491.50 9.13 (0.2,0.8) 18.3 491.50 9.13 (0.4,0.6) 18.3 491.50 9.13 (0.5,0.5) 18.3 491.50 9.13 (0.6,0.4) 18.3 491.50 9.13 (0.8,0.2) 11.5 578.60 10.77 (1.0,0.0) 11.5 578.60 10.77  p = + ∞ (0.0,1.0) 18.3 491.50 9.13 (0.2,0.8) 11.5 578.60 10.77 (0.4,0.6) 11.5 578.60 10.77 (0.5,0.5) 11.5 578.60 10.77 (0.6,0.4) 11.5 578.60 10.77 (0.8,0.2) 11.5 578.60 10.77 (1.0,0.0) 11.5 578.60 10.77 Note: because the optimization result shows thatv-1,2=v-2,3=v-i,i+1…=v-15,16, then we use v- to express average train operating speed in every section.When the operating revenue of the Metro company is more important than the average travel cost of the passengers in the objective function (i.e.,λ1=0.8, λ2=0.2; λ1=1.0, λ2=0), the optimal train operating speed (v-) is approximately 11.5m/s. In that case, the maximum revenue of the Metro company (Z) is about 578.60 CNY per train, and the average travel costs of the passengers (Um) are almost 10.77 CNY. These results imply that the maximum operating revenue is positively related to the optimal speed while the average travel costs of the passengers are negatively correlated with the optimal train operating speed.


## 5.3. Results of the Nondominated Sorting Genetic Algorithm II

We use the “gamultiobj” built-in function in MATLAB to solve our biobjective model with the NSGA-II algorithm. For this case, our parameters are tuned as follows: ParetoFraction = 0.3, population size = 100, generations = 200, “stall Gen Limit” = 200, and “Tol Fun” = 0.01 in the options of function gaoptimset. Here, “ParetoFraction” is the optimal individual coefficient, “stall Gen Limit” is the generation of stopping iteration, and “Tol Fun” is the error of fitness function. Figure6 depicts the Pareto-frontier solutions, in which objective 1 means -Z and objective 2 represents Um. It presents the Pareto solutions based on NSGA-II algorithm, showing the 100 sets of data of v-, min⁡-Z, and min⁡Um, including two extreme optimization results. It can be concluded that (v-=11.50,Z=578.60,Um=10.77) is the best optimal solution from the viewpoint of the company while (v-=18.30,Z=491.54,Um=9.13) is the best optimal solution from the viewpoint of the passengers. The specific optimization results based on NSGA-II algorithm are shown in Table 7.Figure 6
The Pareto solutions obtained by the NSGA-II algorithm.


## 5.4. A Comparative Analysis of the Optimization Algorithms

The fuzzy multiobjective optimization algorithm is to obtain the optimal operating speed by changing the weight and distance parameters in the objective function, so the subjective intention of this algorithm is very significant, and the result will be affected easily by the coefficient values of objective functions, as shown in Table4. Because of the significant differences in the values of those two objective functions, the optimization results thus do not change after increasing the weight of the first objective function.The relationship between average speed and energy consumption is as shown in Figure7, which indicates that operation energy consumption is increasing with the increasing of average speed.Figure 7
The relationship between the average speed and energy consumption.The optimization results of the NSGA-II algorithm, as shown in Figure8, indicate that the higher operating speed benefits the passengers more, and lower operating speed benefits the metro operating company more.Figure 8
The trade-off between the objective functions and the train average operating speed in the sections.


## 5.5. Comprehensive Analysis of the Optimized Timetable

The duration of the running cycle in the current timetable is 3,938 seconds, and the absolute error of total travel time variations between the current timetable and the optimized timetable is restricted to ±120 seconds. The travel time of the current timetable (CUT), the optimized timetable with the operating speed benefitting the company (OTC), and the optimized timetable with the operating speed benefitting the passengers (OTP) respectively are presented in Figure9. Tables 5 and 6 present the timetable of OTC and OTP.Table 5
The optimized timetable of OTC.
Down M1 M2 M3 M4 M5 M6 M7 M8  Dwell (s) - 25 25 25 25 30 25 25  Arrival(s) - 161 330 469 639 741 846 941  Departure(s) 0 186 355 494 664 771 871 966  M9 M10 M11 M12 M13 M14 M15 M16  Dwell (s) 30 30 30 30 30 30 25 -  Arrival(s) 1053 1158 1304 1420 1513 1652 1788 1879  Departure(s) 1083 1188 1334 1450 1543 1682 1813 -  UP M16 M15 M14 M13 M12 M11 M10 M9  Dwell time(s) - 25 30 30 30 30 30 30  Arrival(s) - 2245 2376 2515 2608 2724 2870 2975  Departure(s) 2179 2270 2406 2545 2638 2754 2900 3005  M8 M7 M6 M5 M4 M3 M2 M1  Dwell (s) 25 25 30 25 25 25 25 -  Arrival(s) 3092 3187 3287 3394 3564 3703 3872 4058  Departure(s) 3117 3212 3317 3419 3589 3728 3897 -Table 6
The optimized timetable of OTP.
Down M1 M2 M3 M4 M5 M6 M7 M8  Dwell (s) - 25 25 25 25 30 25 25  Arrival(s) - 147 304 431 589 685 784 873  Departure(s) 0 172 329 456 614 715 809 898  M9 M10 M11 M12 M13 M14 M15 M16  Dwell (s) 30 30 30 30 30 30 25 -  Arrival(s) 979 1076 1216 1326 1415 1546 1672 1759  Departure(s) 1009 1106 1246 1356 1445 1576 1697 -  UP M16 M15 M14 M13 M12 M11 M10 M9  Dwell time(s) - 25 30 30 30 30 30 30  Arrival(s) - 2121 2242 2373 2462 2572 2712 2809  Departure(s) 2059 2146 2272 2403 2492 2602 2742 2839  M8 M7 M6 M5 M4 M3 M2 M1  Dwell (s) 25 25 30 25 25 25 25 -  Arrival(s) 2920 3009 3103 3204 3362 3489 3646 3818  Departure(s) 2945 3034 3133 3229 3387 3514 3671 -Table 7
The optimization results based on NSGA-II algorithm.
v - min ⁡ ( - Z ) min ⁡ U m v - min ⁡ ( - Z ) min ⁡ U m v - m i n ( - Z ) min ⁡ U m v - min ⁡ ( - Z ) min ⁡ U m 11.50 -578.60 10.77 16.67 -516.00 9.40 17.29 -506.88 9.29 15.76 -528.74 9.58 18.30 -491.54 9.13 13.42 -558.03 10.14 11.73 -576.25 10.68 13.53 -556.77 10.11 18.30 -491.54 9.13 15.45 -532.82 9.64 14.03 -550.82 9.98 12.48 -568.53 10.43 12.70 -566.12 10.36 11.78 -575.84 10.67 16.29 -521.35 9.47 16.11 -523.93 9.51 14.02 -550.90 9.98 17.85 -498.48 9.20 14.53 -544.73 9.85 13.80 -553.56 10.04 16.95 -511.95 9.35 18.22 -492.76 9.14 15.24 -535.67 9.69 14.85 -540.71 9.78 11.50 -578.60 10.77 15.09 -537.54 9.72 12.09 -572.61 10.56 16.22 -522.29 9.49 11.86 -575.00 10.64 16.40 -519.83 9.45 12.30 -570.44 10.49 16.92 -512.42 9.36 15.30 -534.80 9.68 12.19 -571.59 10.52 12.35 -569.87 10.47 13.96 -551.73 10.00 18.09 -494.74 9.16 17.22 -507.98 9.31 15.16 -536.70 9.71 17.60 -502.37 9.24 16.33 -520.82 9.47 12.81 -564.87 10.32 14.10 -549.96 9.96 15.58 -531.07 9.62 12.02 -573.29 10.58 12.51 -568.17 10.42 11.98 -573.79 10.60 16.06 -524.62 9.52 13.90 -552.45 10.01 17.77 -499.67 9.21 14.60 -543.87 9.84 17.53 -503.37 9.25 14.26 -548.05 9.92 11.63 -577.30 10.72 12.86 -564.40 10.31 13.84 -553.12 10.03 15.71 -529.34 9.59 16.80 -514.11 9.38 17.45 -504.56 9.27 15.39 -533.69 9.66 17.34 -506.26 9.29 17.96 -496.86 9.18 13.58 -556.17 10.10 13.65 -555.37 10.08 14.81 -541.15 9.79 15.80 -528.17 9.57 14.75 -541.91 9.80 11.89 -574.69 10.63 16.88 -512.92 9.36 14.92 -539.80 9.76 11.55 -578.15 10.75 13.24 -560.09 10.19 12.73 -565.75 10.34 16.58 -517.35 9.42 15.49 -532.25 9.63 13.01 -562.66 10.26 14.70 -542.51 9.81 14.43 -545.89 9.87 11.67 -576.93 10.71 16.80 -514.11 9.38 17.07 -510.20 9.33 14.20 -548.74 9.93 14.33 -547.20 9.90 15.03 -538.32 9.73 13.06 -562.17 10.25 16.48 -518.74 9.44 12.40 -569.32 10.45 18.29 -491.66 9.13 11.59 -577.75 10.74 17.17 -508.79 9.32 13.15 -561.09 10.22 15.16 -536.60 9.71 14.33 -547.18 9.90 17.66 -501.37 9.23 13.75 -554.23 10.05 17.08 -510.08 9.33 12.78 -565.24 10.33 12.26 -570.86 10.50 12.12 -572.28 10.55 14.47 -545.49 9.87Figure 9
Comparison between the optimized timetable and the current timetable.The optimization results indicate that the operating revenue (Z) in CUT, OTC, and OTP is 536.52 CNY, 544.32 CNY, and 527.73 CNY, respectively. The average travel costs of the passengers (Um) in CUT, OTC, and OTP are 9.70 CNY, 9.84 CNY, and 9.56 CNY, respectively. Comparing with the OTC and CUT, the revenue is increased by 1.45%, and the average energy consumption of each train in sections is reduced by 7.88%. Considering maximizing operating revenue, we conclude that the OTC performs better than the current one CUT and OTP.


## 5.6. The Sensitivity of the Main Parameters

In order to research the relationship between the main parameters and decision variablesv-i,i+1 in the model, we analyzed the sensitivity of the main parameters (the electricity rate η and passengers’ value of time VOT) affecting train running costs and passengers’ average travel cost in the model. Assuming the sensitivity range of η is 0.5-1.5 CNY/(kw·h), the step length is 0.1 CNY/(kw·h), and the sensitivity range of VOT is 18-30 CNY/h, the step length is 1 CNY/h. Therefore, we get the relationship between the two parameters and decision variables v-i,i+1 under different weights at p=1, as shown in Figure 10.Figure 10
The trade-off between the main parameters and the train average operating speed in the sections.The train average operating speed in the sections decreases with the increase of the electricity rateη and decreases with the increase of λ1 when the electricity rate η is greater than 1.3 CNY/(kw·h). Also, the train average operating speed in the sections increases with the increase of passengers’ value of time VOT and increases with the decrease of λ1 when passengers’ value of time VOT is less than 23 CNY/h.


## 6. Conclusions

In this paper, we proposed a biobjective mathematical model to find an improved timetable by optimizing the average train operating speed in sections. The objectives are to maximize the revenue of the system planner who operates trains while trying to minimize the average travel costs of the passengers. The optimization model observes the actual passengers boarding rate to satisfy the operational requirements. Also, the paper studies the trade-off between the maximum running speed and the average train operating speed in sections based on the identified parameters and decision variables. In the solution process, the fuzzy multiobjective optimization algorithm is adopted as an effective method for the optimization model. Moreover, we also presented the nondominated sorting genetic algorithm II that seeks the relationship between decision variables and objective functions. Both algorithms can effectively reduce the effort to obtain acceptable solutions. As a case study, the numerical experiments of Chengdu Metro showed that higher operating speed benefits passengers more, while lower operating speed benefits the operating company more. Comparing with the OTC that obtained by the proposed model and algorithms and CUT, the OTC performs better than the current one CUT as the revenue is increased by 1.45% and the average operating energy consumption of each train in sections is reduced by 7.88%. To satisfy the more operational requirements, we will further strengthen the proposed model to integrate the timetable optimization and other practical constraints.



---

*Source: 7805168-2018-12-16.xml*