# Optimal Differential Pricing for Intercity High-Speed Railway Services with Time-Dependent Demand and Passenger Choice Behaviors under Capacity Constraints

Authors: Huanyin Su; Shuting Peng; Lianbo Deng; Weixiang Xu; Qiongfang Zeng
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Differential pricing of trains with different departure times caters to the taste heterogeneity of the time-dependent (departure time) demand and then improves the ticket revenue of railway enterprises. This paper studies optimal differential pricing for intercity high-speed railway services. The distribution features of the passenger demand regarding departure times are analyzed, and the time-dependent demand is formulated; a passenger assignment method considering departure periods and capacity constraints is constructed to evaluate the prices by simulating the ticket-booking process. Based on these, an optimization model is constructed with the aim of maximizing the ticket revenue and the decision variables for pricing train legs. A modified direct search simulated annealing algorithm is designed to solve the optimization model, and three random generation methods of new solutions are developed to search the solution space efficiently. Experimental analysis containing dozens of trains is performed on Wuhan-Shenzhen high-speed railway in China, and price solutions with different elastic demand coefficients (ϕ) are compared. The following results are found: (i) the optimization algorithm converges stably and efficiently and (ii) differentiation is shown in the price solutions, and the optimized ticket revenue is influenced greatly by ϕ, increasing by 7%–21%.


---

## Body

## 1. Introduction

Intercity high-speed railways serve passengers that are located in densely populated metropolitan areas or city clusters. They feature high departure frequencies and high speeds, not less than 250 km/h. Since 2016, the ticket prices of part high-speed trains have been floated within certain limits on several lines (Beijing-Shanghai high-speed railway and Hangzhou-Shenzhen railway) because pricing power has been devolved to the railway transport enterprise in China. However, train prices of most high-speed lines in China are still unitary for the same train leg and fixed now, which have seriously affected the market competitiveness of intercity high-speed railways. Hence, more theoretical methods are needed to guide pricing. According to the passenger characteristics in the railway transportation market, optimal pricing methods for intercity high-speed railways are studied in this paper, which are of great significance for improving railway operation and management mechanisms and increasing the economic revenue of railway transportation enterprises.Determining the optimal pricing of intercity high-speed railways is a revenue management problem, and revenue management methods have been studied for decades and widely applied in the railway industry in Europe and Japan [1, 2]. Researches related to revenue management have mainly focused on air transport industry, and revenue management methods for intercity high-speed railways are under-researched. However, there are great differences between airlines and high-speed railways. For example, few intermediate stops are included for airline services, whereas many stop stations are dotted in train lines. Hence, the pricing problem for high-speed railways is much more complicated.Considering the complexity of pricing problem for high-speed railways, fewer researchers studied optimal pricing based on a railway network scale or multiple trains. However, all trains in an intercity high-speed railway network form an integrated passenger transport system. Multiple trains operate together to serve passengers with different departure times, and passengers make choices between trains under capacity constraints. Hence, the optimal pricing of intercity high-speed railway services should be studied on a network scale.The passenger demand is premise or basis of optimal pricing problems. Based on historical ticket reservation data, the passenger demand for different departure times in a day shows fluctuation rules, and the rules are stable over a short period. Such a demand can be called time-dependent demand. The time-dependent demand has been studied by various researchers [3–6], and a train timetable or line planning has been evaluated and optimized based on the passenger choice behaviors with different departure times. However, time-dependent demand has been under-researched in the optimal pricing studies of intercity high-speed railways.Passengers are usually characterized by a high level of taste heterogeneity considering their different preferences for scheduling and pricing [7, 8]. Time-dependent demand has heterogeneous behavior preferences for trains with different departure times. Similarly, trains with different departure times serve passengers with different preferences. Hence, train legs with different departure times can be priced differently. Differential pricing of trains caters to the taste heterogeneity of the time-dependent demand and will improve the ticket revenue of railway enterprises.In light of the above, this paper proposes an optimal differential pricing method for intercity high-speed railway services considering the time-dependent demand and passenger choice behaviors under capacity constraints. The main contributions are as follows: (i) time-dependent demand features are analyzed and formulated, and a multinomial logit model is applied to predict the passenger choice behaviors between trains considering their departure times; (ii) an optimal differential pricing model is constructed considering capacity constraints and time-dependent demand; (iii) a modified direct search simulated annealing algorithm is developed to solve the optimization model, and three random generation methods of new solutions are designed to search the solution space efficiently, and a passenger assignment method by simulating the ticket-booking process is proposed to evaluate the price solutions; and (iv) experimental analysis is performed to prove the stability and efficiency of the algorithm convergence, and price solutions with different elastic demand coefficients are compared to verify the feasibility of the optimization method.In Table1, the main notations used in the paper are listed in order of appearance. The remainder of the paper is organized as follows. A review of the literature on revenue management in the railway industry is presented in the next section. In Section 3, the time-dependent demand is analyzed and formulated. In Section 4, a travel choice model considering departure periods and capacity constraints is proposed. In Section 5, an optimal differential pricing model is constructed and a modified direct search simulated annealing algorithm is developed to solve the optimization model. In Section 6, a case study is conducted on Wuhan-Shenzhen high-speed railway in China to illustrate the optimized method and verify its feasibility. In Section 7, the conclusions and future research directions are presented.Table 1
The main notations used.
r,sAn origin-destination (O-D) pair;r is the departure station and s is the arrival station for passengersqrsThe demand of an O-D pairr,skThe departure periodfrskThe probability of departing in thekth period for the demand of an O-D pair r,sqrskThe time-dependent demand of an O-D pairr,s departing in the kth periodSNThe set of stationsRSThe set of all O-D pairsjA trainWThe set of all trainshjThe number of stop stations of trainjvxjThexth stop station of train jSTjThe sequence of stop stations of trainjvxj,vyjA leg of trainj between stop stations vxj and vyjLEGjThe set of all the legs of trainjCAjThe capacity of trainjtrsjThe running time of trainj on leg r,scrsjThe price of trainj on leg r,s, the decision variableWrsThe set of trains stopping at stationr and stopping at station strsjThe running time of trainj on leg r,scrsjThe price of trainj on leg r,s, the decision variabledjThe departure hour of trainjηmThe assignment ratio of the demand in themth stageWrsmThe set of alternative trains for O-D pairr,s in the mth stageqrsk,mThe accepted passenger demand for O-D pairr,s and period k in the mth stagePj|r,s,k,mThe probability of selecting trainj from set Wrsk,m for time-dependent demand qrsk,mFj,x,mThe number of passengers on legvxj,vx+1j at the end of the mth stageMThe number of stages of the ticket-booking processcnThe number of decision variablesCThe vector of the decision variablesCLThe vector of the lower limits of the decision variablesCUThe vector of the upper limits of the decision variablesZTicket revenue


## 2. Literature Review

In recent years, with the development of high-speed railways, more and more scholars have paid attention to the pricing problems in the high-speed railway industry [9–11]. However, most studies have considered only one train.Zheng and Liu [12] developed a ticket fare optimization model for one high-speed train and determined the number of fare grades and the price of each fare grade in the sale horizon. Hetrakul and Cirillo [8] studied the revenue optimization problem considering the pricing and seat allocation jointly. Multinomial and latent class models were used to explain the ticket purchase timing of passengers, but only one train was considered in the optimization model. Fransiscus [13] proposed a practical revenue management method of railways including three stages. Also, only one train was considered in this method.Qin et al. [14] used prospect theory to construct a differentiated pricing model under elastic demand and applied a simulated annealing algorithm to solve the model, and eight trains are tested. Some researchers [15, 16] jointly studied pricing and seat allocation optimization problems including several trains, and a multinomial model was used to analyze the passenger choice behaviors between trains. Hu et al. [17] also proposed a joint optimization method of pricing and seat allocation in high-speed rail networks and tested a large-scale instance.In the above studies, time-dependent demand was not considered deeply. Usually, ticket fares and booking days were taken as factors in the choice models, and deterministic demand using historical data and stochastic demand following a Poisson process and elastic demand were assumed [18–20]. Actually, time-dependent demand had already been studied in some passenger transportation organization optimization problems. For example, time-dependent demand was considered in the train timetable optimization by Kaspi and Raviv [3] and Niu et al. [4] and also used in the line planning problem by Su et al. [5]. Xu et al. [6] designed passenger assignment methods for high-speed railways based on time-dependent demand. Nevertheless, time-dependent demand regarding the pricing problems of high-speed railways has not been studied sufficiently.Discrete choice models were applied to analyze passenger choice behaviors for the above pricing problems. There is only once chance in the above studies when passengers were assigned to trains by the discrete choice models. Passengers were accepted to aboard trains if there were residual capacities for their choices. Otherwise, they were refused to aboard trains. However, the choices are not unique; if there are no residual capacities for a choice, passengers can be assigned to other trains with residual capacities again, which means passengers have several chances to make choice. This is a passenger assignment problem of railways. Passenger assignment problems have been widely studied for urban transit systems [21–23]. For example, Sumalee et al. [24] proposed a dynamic transit assignment model that had an explicit seat allocation process. Hamdouch et al. [25] proposed a schedule-based transit assignment model with travel strategies and loaded passengers on a first-come-first-served basis. However there are big differences between the urban transit systems and high-speed rail systems, which have been analyzed in detail by Su et al. [26]. Considering the features of high-speed rail systems, Su et al. [26] proposed a schedule-based passenger assignment method for high-speed rail networks by simulating the ticket-booking process, which can be used for a reference in this article.As previously stated, the research status of the optimal pricing of railway revenue management can be summarized as follows: (1) many researchers have focused on the optimal pricing problem of one train in the sale horizon, and optimal differential pricing of multiple trains has not been studied extensively; (2) the passenger demand is generally formulated as elastic demand or stochastic demand in the sale horizon, and the time-dependent demand deserves more attention in revenue management research; (3) discrete choice models were widely applied to assign passengers and passenger assignment methods considering the features of high-speed rail systems deserve more attention; and (4) for large-scale instances, the convergence efficiency of the solution algorithm need to be improved. The comparison of our research with the literature in pricing problems of high-speed railways is shown in Table2.Table 2
Comparison of our research with the literature in pricing problems of railways.
ReferencesTime-dependent demandTaste heterogeneityMultiple trainsSolution methodHetrakul and Cirillo [8]✕✓✕LINGO 12.0Zhang et al. [11]✕✓✕An inverse dynamic programming methodZheng and Liu [12]✕✕✕Lagrange multiplier and KKT conditionsFransiscus [13]✕✕✕EnumerationQin et al. [14]✕✓✓Simulated annealing algorithmQin et al. [15]✕✕✓Artificial bee colony algorithmWu et al. [16]✕✕✓Two-stage algorithmHu et al. [17]✕✕✓Quasi-Newton method-based solution algorithmThis paper✓✓✓Modified direct search simulated annealing algorithmConsequently, this paper aims to fill the gap by proposing an optimal differential pricing method for intercity high-speed railway services, considering the time-dependent demand and passenger choice behaviors under capacity constraints.


## 3. Time-Dependent Demand

In this section, we analyze the fluctuation rules of passenger demands with different departure times in a day and determine the formulation of the time-dependent demand. The analysis is based on the historical ticket reservation data of the Wuhan-Shenzhen railway in the first half of 2016. This railway is located in Central China and South China with a length of 1171 km, 18 stations, and a design speed of 350 km/h, as shown in Figure1. The passenger demand along this railway is great. With high departure frequency, passengers can decide their departure time according to their will. The passenger demand for different departure times in a day shows fluctuation rules.Figure 1
Wuhan-Shenzhen high-speed railway.The passenger demand is analyzed statistically in dimensions of departure times and days of a week. The average passenger demand per hour in each day of one week is calculated, and the results of some origin-destination (O-D) pairs (Wuhan-Guangzhounan, Changshanan-Shenzhenbei, Changshanan-Guangzhounan, and Guangzhounan-Shenzhenbei) are shown in Figure2.Figure 2
Passenger demand distribution. (a) From Wuhan to Guangzhounan. (b) From Changshanan to Shenzhenbei. (c) From Changshanan to Guangzhounan. (d) From Guangzhounan to Shenzhenbei.
(a)(b)(c)(d)According to Figure2, passenger demand floats along departure times and shows some rules as follows:(i)
The “peaks and valleys” of the passenger demand occur at different departure times in one day and form some variation rules.(ii)
The variation rules are similar for each day of week for an O-D pair.(iii)
The variation rules are different for different O-D pairs.The variation rules in the above analysis reflect some time-dependent features of passenger demand. Thus, passenger demand with the above variation rules can be called time-dependent demand. For the sake of description, time-dependent demand can be formulated by calculating the probability of passenger demand choosing different departure times using statistical regression methods based on historical ticket reservation data, which are defined in detail as follows.The service time horizont1,t2 is divided into periods by hours; for example, [6:00, 6:59] is called period 6, [7:00, 7:59] is called period 7, and so on. Let RS be the set of all O-D pairs. The passenger demand of O-D pair r,s in a day is denoted as qrs, r,s∈RS, and the probability of the passenger demand departing in the kth period is denoted as frsk, so the passenger demand departing in the kth period can be calculated as follows:(1)qrsk=frsk×qrs,k=t1,t1+1,…,t2−1;r,s∈RS.Therefore, equation (1) is the formulation of the time-dependent demand and period k is the expected departure hour of passenger demand qrsk.


## 4. Passenger Assignment considering Departure Periods and Capacity Constraints

The analysis of passenger assignment aims to calculate the number of passengers on each train more accurately and then to evaluate the differentiated pricing of trains. In this section, a passenger assignment method is designed by simulating the ticket-booking process to assign passengers to trains with limited capacities. In addition, a multinomial logit model is adopted to determine the choices between trains for time-dependent demand.A passenger assignment method can be constructed by simulating ticket-booking process because passengers’ choices of trains are determined in this process [26]. As the ticket-booking process continues, some train legs become fully occupied and the alternative trains for the subsequent passengers’ booking requests change. Hence, the ticket-booking process is partitioned into several stages and the passenger assignment is conducted stage by stage. In each stage, the alternative trains are constant for the unassigned passengers. When all the time-dependent demands are assigned or there is no train with residual capacity to satisfy the unassigned passengers, the passenger assignment process is over.For ease of presentation, some notations are expressed as follows. LetSN be the set of all stations. The train set is expressed as W=j. For train j, hj is the number of stop stations and STj=v1j,v2j,…,vhjj is the sequence of stop stations (vxj∈SN). The railway line section between stop station vxj and stop station vyj (x<y) for train j is called one leg, denoted as vxj,vyj. The set of all the legs of train j is denoted as LEGj. For x≤x′ andy′≤y, leg vx′j,vy′j is covered by leg vxj,vyj on the line, denoted as vx′j,vy′j⊆vxj,vyj. The capacity of train j is denoted as CAj. The travel time and price of leg r,s of train j are denoted as trsj and crsj, respectively. The departure hour of train j is denoted as dj. The set of trains stopping at stations r and s is expressed as Wrs.In the passenger assignment method, the ticket-booking process is partitioned intoM stages. In the mth stage, m=1,2,…,M, the available train set is denoted as Wrsm and the accepted passenger demand is denoted as qrsk,m for O-D pair r,s and period k. The accepted passenger demand in the whole ticket-booking process will not exceed the real passenger demand, as shown in the following equation:(2)∑m=1Mqrsk,m≤qrsk,r,s∈RS;k=t1,t1+1,…,t2−1.The seat allocation optimization problem is also quite complex for multiple trains, which is not considered in this paper. Hence, we set the same accepted ratio (the proportion of the accepted passenger demand to the total demand of each O-D pair) for all O-D pairs, denoted asηm for the mth stage in the ticket-booking process. Then, we have(3)qrsk,m=ηmqrsk,r,s∈RS;k=t1,t1+1,…,t2−1;m=1,2,…,M,(4)0≤ηm≤1,m=1,2,…,M.Passenger choices between trains are influenced by many factors, and in this paper, we mainly consider departure time deviation, travel time, and ticket price. Then, according to the multinomial logit model theory widely used in the studies [15, 17], the probability to select train j for passenger demand qrsk,m is denoted as Pj|r,s,k,m, which is calculated as follows:(5)Pj|r,s,k,m=expϑVjr,s,k−1∑j∈WrsmexpϑVjr,s,k−1,j∈Wrsm;r,s∈RS;k=t1,t1+1,…,t2−1;m=1,2,…,M,(6)Vjr,s,k=θ1×k−dj+θ2×trsj+crsj,whereVjr,s,k is the generalized travel cost; k−dj is the absolute value of departure time deviation between the expected departure time period k and the real departure time period dj; θ1 and θ2 are time value parameters; and ϑ is a positive scaling factor. If the departure time of train j is not in period k, then k−dj>0, which results in more travel cost for passenger demand qrsk,m. In this way, departure time preferences of passengers are considered in the passenger assignment method.At the end of themth stage, the number of passengers on leg vxj,vx+1j is denoted as Fj,x,m,1≤x≤hj−1, j∈Wrsm. Then, the following expressions are satisfied:(7)Fj,x,m=Fj,x,m−1+∑r,s:vxj,vx+1j⊆r,sandj∈Wrsm∑k=t1t2−1qrsk,mPj|r,s,k,m,1≤x≤hj−1;j∈W;m=1,2,…,M,(8)Fj,x,0=0,(9)Fj,x,m≤CAj,1≤x≤hj−1;j∈W;m=1,2,…,M.Equations (7) and (8) calculate the passengers on each train leg. Equation (9) guarantees the passengers on each train leg will not exceed the train capacity at the end of each phase in the ticket-booking process.Based on the above formulation, the passenger assignment method is constructed by simulating the ticket-booking process as shown below:Step 0: calculate the time-dependent demandqrsk, r,s∈RS, k=t1,t1+1,…,t2−1, by equation (1). Let phase m=0.Step 1: letm←m+1. For r,s∈RS, construct the available train set Wrsm, Wrsm=j|Fj,x,m−1<CAj,vxj,vx+1j⊆r,s,j∈Wrs. If all Wrsm=∅, r,s∈RS, the simulation is over, and let M=m−1; otherwise, go to Step 3.Step 2: forr,s∈RS, k=t1,t1+1,…,t2−1, j∈Wrsm, calculate the probability Pj|r,s,k,m by equations (5) and (6).Step 3: calculate the accepted ratioηm. Let(10)ΔFmj,x=∑r,s:vxj,vx+1j⊆r,sandj∈Wrsm∑k=t1t2−1qrskPj|r,s,k,m,1≤x≤hj−1;j∈W,whereΔFmj,x represents the number of passengers’ booking requests for leg vxj,vx+1j in the mth stage. Combining equation (3) and (7)–(9), we have(11)ηm≤CAj−Fj,x,m−1ΔFmj,x,ΔFmj,x>0;1≤x≤hj−1;j∈W.Besides, the accepted passenger demand until the current phase cannot exceed the total passenger demand; then, we have(12)∑l=1mηlqrsk≤qrsk,r,s∈RS,k=t1,t1+1,…,t2−1.It is feasible forηm satisfying equation (11) and (12), so let(13)ηm=min1−∑l=1m−1ηl,CAj−Fj,x,m−1ΔFmj,x|ΔFmj,x>0,1≤x≤hj−1,j∈W,η0=0.Step 4: calculate the accepted passenger demandqrsk,m, r,s∈RS, k=t1,t1+1,…,t2−1, by equation (3). Calculate passenger demand on each train leg Fj,x,m, 1≤x≤hj−1,j∈W, by equations (7) and (8). If ∑l=1mηl=1, the simulation is over, and let M=m; otherwise, go to Step 3.Based on Step 0 to Step 4, the ticket-booking process is simulated and passengers are assigned to each train leg. In brief, the above passenger assignment method considering departure periods and capacity constraints is denoted asPACC, which will be used to evaluate the differentiated pricing of trains in the following sections.


## 5. Optimal Differential Pricing Model and Solution Algorithm

An optimal differential pricing model for multiple trains is constructed with the objective of maximizing the ticket revenue, and a heuristic optimization algorithm is designed to solve the model.


### 5.1. Model Formulation

The price of each train leg denoted ascrsj, r,s∈LEGj, j∈W, is a decision variable and is optimized collaboratively in the optimal model. Ticket prices can influence the passengers’ choice between trains; however, if the price is too high, some passengers may consider other alternative modes, such as intercity automobiles. In other words, ticket prices can influence whether passengers choose trains; such a demand is called the elastic price demand. An exponential function is used to formulate the total passenger demand of each O-D pair regarding ticket prices as follows:(14)qrs=qrs0exp−ϕc¯rsc¯rs0−1,where c¯rs is the average price of leg r,s for all trains. If the average price c¯rsis equal to c¯rs0, the elastic demand qrs is qrs0. c¯rs0 and qrs0 can be determined based on historical ticket reservation data. ϕ is an elastic coefficient. The bigger ϕ is, the greater the demand elasticity becomes.For a given price solution, the average price of legr,s for trains departing in the kth period can be calculated simply as follows:(15)c¯rsk=∑j∈WrskcrsjWrsk.The probability of passengers departing in thekth period is frsk; then, the average price of leg r,s for all the trains perceived by passengers can be calculated as follows:(16)c¯rs=∑k=722frskc¯rsk.Based on the above description, the optimal differential pricing model is expressed as follows:(17)maxZ=∑m=1M∑r,s∈RS∑j∈Wrsm∑k=t1t2−1qrsk,mPj|r,s,k,mcrsj,s.t.(18)crsmin≤crsj≤crsmax,r,s∈LEGj,j∈W,where Z is the ticket revenue, which can be calculated by the PACC method, and crsmin and crsmax represent the lower limit and upper limit of the price of leg r,s, respectively.


### 5.2. Solution Algorithm

The number of decision variables in the above optimization model is great. Hence, we need to develop an algorithm to solve the combinational optimization problem with high dimensions. The simulated annealing (SA) algorithm converges globally and is used widely in engineering practices. Ali et al. [27] proposed the direct search simulated annealing (DSA) algorithm, which improved the calculation efficiency and accuracy of the DSA algorithm compared with the SA algorithm. To improve the convergence stability and efficiency for solving the optimization problem with high dimensions, a modified DSA (MDSA) algorithm is proposed.


#### 5.2.1. Solution Generation

Three methods of generating new solutions are included in the MDSA algorithm. The first one generates a new solution randomly in a search space based on Cauchy distribution with a variable scale parameter [28], called GNS 1. The second one generates a new solution by selecting a certain number of solutions randomly from the solution set and combining them with the current optimized solution [27], called GNS 2. The third one generates a new solution by searching randomly around the current optimized solution based on the standard normal distribution [28], called GNS 3.For ease of description, the number of decision variables is denoted ascn. C is the column vector of the decision variables, which is also called a solution vector. CL is the column vector of the lower limits of the decision variables. CU is the column vector of the upper limits of the decision variables. Then,(19)CL≤C≤CU.All the solution vectors satisfying equation (19) form a search space of solutions, denoted as S. Based on the constraints in the optimization model, the solution vectors in S are feasible solutions. ZC is used as the evaluation function in the MDSA algorithm. Let N=7cn+1. N solutions are generated randomly in search space S and form an initial solution set A. In set A, the best solution is CBest and the worst solution is CWorst. The evaluation function values are denoted as Zmax and Zmin, respectively.Three methods of generating new solutions are described in detail below.(1) GNS 1. A new solution C is generated randomly in search space S based on Cauchy distribution with a variable scale parameter. The Cauchy distribution probability density function centered at the origin is shown as follows: (20)fξ=1πφφ2+ξ2,−∞<ξ<∞,where φ is a variable scale parameter determined by equation (21) based on the number of cooling iterations in the MDSA algorithm, which leads to searching in a wide range early and on a small scale later in the optimization process.(21)φ=1.1a,a≤50,0.9a−50,a>5,where a is the number of cooling iterations.New solutionC is generated by the following:(22)C=CL+CU−CL∘absy,w∈0,0.33,CL+CU2+CU−CL2∘y,w∈0.33,0.67,CU−CU−CL∘absy,w∈0.67,1,(23)y=AYAY,where symbol “∘” in equation (22) represents the Hadamard product. AY is a column vector with cn dimensions generated randomly by equations (20) and (21), and y is a unit column vector with cn dimensions. abs∗ represents the absolute value function. w is a random number generated randomly in section 0,1 based on the uniform distribution.(2) GNS 2. RN solutions are selected randomly from the solution set, denoted as Cn1,Cn2,…,CnRN. Then, these solutions are combined randomly with the current optimized solution CBest. Each element ci (i=1,2,…,cn) in new solution C is generated by the following:(24)ci=cink,wi∈k−1RN+1,kRN+1,k=1,2,…,RN,ciBest,wi∈RNRN+1,1,where wi is a random number generated randomly in section 0,1 based on the uniform distribution. RN makes up a proportion of the number of solutions in the solution set. This method helps greatly in improving the convergence efficiency and stability.(3) GNS 3. A new solution C is generated by searching randomly around the current optimized solution CBest based on the standard normal distribution as follows:(25)C=CBest+bI∘CU−CL∘y,where symbol “∘” in equation (25) presents the Hadamard product. b is an adjustment factor of the search radius, and I is an adaptive factor. y is a unit column vector with cn dimensions. All of the above variables are calculated using equations (26)–(28).(26)b=arctanTa−0.5ma,Ta≥1,0.05,Ta<1,(27)Ii=1γi<β0γi≥β,i=1,2,…,cn,(28)y=SYSY,where Ta is the current temperature. ma and γi are random numbers generated randomly in section 0,5 and 0,1, respectively, based on the uniform distribution. β is an adaptive probability and takes 0.1 generally. SY is a column vector with cn dimensions generated randomly based on the standard normal distribution.New solutionC generated by equation (25) may not be in search space S, and, in this situation, each element ci (i=1,2,…,cn) of new solution C is amended by(29)ci=ciU−modciL−ci,ciU−ciL,ci<ciL,ci,ciL≤ci≤ciU,ciL+modci−ciU,ciU−ciL,ci>ciU.


#### 5.2.2. MDSA Algorithm Structure

In the MDSA algorithm, the initial temperatureT0 satisfies a requirement that new solutions are accepted in the early stage with a probability greater than the desired value (between 0.8 and 0.9). The refresh method of Ta refers to the paper by Ali et al. [27]. The length of Markov chain La is calculated by equation (30), which is the maximum number of inner iterations.(30)La=10cn+10cn1−eZmin−Zmax.The acceptance probability of new solutionC is determined by the following equation:(31)Pacc=1,ZC≥Zmin,eZC−Zmin/Ta,ZC<Zmin.LetΔZmax=Zmaxa+1−Zmaxa.δand τ are given positive real numbers. If ΔZmax≤δ,for successiveΔacooling timesΔa<τ, which means Zmax does not change or significantly change after a certain cooling, the MDSA algorithm is judged to enter a local extremum. The structure of the MDSA algorithm is shown in Figure 3.Figure 3
The structure of the MDSA algorithm.In the early stage (a≤50) of MDSA algorithm, new solutions are generated mainly by GNS 1, which can expand the search range and improve the coverage of solution set A in search space S. In the later stage (a>50) of MDSA algorithm, new solutions are generated mainly by GNS 2, which can improve the convergence efficiency and stability. GNS 3 can avoid the MDSA algorithm falling into a local extreme too early.


## 5.1. Model Formulation

The price of each train leg denoted ascrsj, r,s∈LEGj, j∈W, is a decision variable and is optimized collaboratively in the optimal model. Ticket prices can influence the passengers’ choice between trains; however, if the price is too high, some passengers may consider other alternative modes, such as intercity automobiles. In other words, ticket prices can influence whether passengers choose trains; such a demand is called the elastic price demand. An exponential function is used to formulate the total passenger demand of each O-D pair regarding ticket prices as follows:(14)qrs=qrs0exp−ϕc¯rsc¯rs0−1,where c¯rs is the average price of leg r,s for all trains. If the average price c¯rsis equal to c¯rs0, the elastic demand qrs is qrs0. c¯rs0 and qrs0 can be determined based on historical ticket reservation data. ϕ is an elastic coefficient. The bigger ϕ is, the greater the demand elasticity becomes.For a given price solution, the average price of legr,s for trains departing in the kth period can be calculated simply as follows:(15)c¯rsk=∑j∈WrskcrsjWrsk.The probability of passengers departing in thekth period is frsk; then, the average price of leg r,s for all the trains perceived by passengers can be calculated as follows:(16)c¯rs=∑k=722frskc¯rsk.Based on the above description, the optimal differential pricing model is expressed as follows:(17)maxZ=∑m=1M∑r,s∈RS∑j∈Wrsm∑k=t1t2−1qrsk,mPj|r,s,k,mcrsj,s.t.(18)crsmin≤crsj≤crsmax,r,s∈LEGj,j∈W,where Z is the ticket revenue, which can be calculated by the PACC method, and crsmin and crsmax represent the lower limit and upper limit of the price of leg r,s, respectively.


## 5.2. Solution Algorithm

The number of decision variables in the above optimization model is great. Hence, we need to develop an algorithm to solve the combinational optimization problem with high dimensions. The simulated annealing (SA) algorithm converges globally and is used widely in engineering practices. Ali et al. [27] proposed the direct search simulated annealing (DSA) algorithm, which improved the calculation efficiency and accuracy of the DSA algorithm compared with the SA algorithm. To improve the convergence stability and efficiency for solving the optimization problem with high dimensions, a modified DSA (MDSA) algorithm is proposed.


### 5.2.1. Solution Generation

Three methods of generating new solutions are included in the MDSA algorithm. The first one generates a new solution randomly in a search space based on Cauchy distribution with a variable scale parameter [28], called GNS 1. The second one generates a new solution by selecting a certain number of solutions randomly from the solution set and combining them with the current optimized solution [27], called GNS 2. The third one generates a new solution by searching randomly around the current optimized solution based on the standard normal distribution [28], called GNS 3.For ease of description, the number of decision variables is denoted ascn. C is the column vector of the decision variables, which is also called a solution vector. CL is the column vector of the lower limits of the decision variables. CU is the column vector of the upper limits of the decision variables. Then,(19)CL≤C≤CU.All the solution vectors satisfying equation (19) form a search space of solutions, denoted as S. Based on the constraints in the optimization model, the solution vectors in S are feasible solutions. ZC is used as the evaluation function in the MDSA algorithm. Let N=7cn+1. N solutions are generated randomly in search space S and form an initial solution set A. In set A, the best solution is CBest and the worst solution is CWorst. The evaluation function values are denoted as Zmax and Zmin, respectively.Three methods of generating new solutions are described in detail below.(1) GNS 1. A new solution C is generated randomly in search space S based on Cauchy distribution with a variable scale parameter. The Cauchy distribution probability density function centered at the origin is shown as follows: (20)fξ=1πφφ2+ξ2,−∞<ξ<∞,where φ is a variable scale parameter determined by equation (21) based on the number of cooling iterations in the MDSA algorithm, which leads to searching in a wide range early and on a small scale later in the optimization process.(21)φ=1.1a,a≤50,0.9a−50,a>5,where a is the number of cooling iterations.New solutionC is generated by the following:(22)C=CL+CU−CL∘absy,w∈0,0.33,CL+CU2+CU−CL2∘y,w∈0.33,0.67,CU−CU−CL∘absy,w∈0.67,1,(23)y=AYAY,where symbol “∘” in equation (22) represents the Hadamard product. AY is a column vector with cn dimensions generated randomly by equations (20) and (21), and y is a unit column vector with cn dimensions. abs∗ represents the absolute value function. w is a random number generated randomly in section 0,1 based on the uniform distribution.(2) GNS 2. RN solutions are selected randomly from the solution set, denoted as Cn1,Cn2,…,CnRN. Then, these solutions are combined randomly with the current optimized solution CBest. Each element ci (i=1,2,…,cn) in new solution C is generated by the following:(24)ci=cink,wi∈k−1RN+1,kRN+1,k=1,2,…,RN,ciBest,wi∈RNRN+1,1,where wi is a random number generated randomly in section 0,1 based on the uniform distribution. RN makes up a proportion of the number of solutions in the solution set. This method helps greatly in improving the convergence efficiency and stability.(3) GNS 3. A new solution C is generated by searching randomly around the current optimized solution CBest based on the standard normal distribution as follows:(25)C=CBest+bI∘CU−CL∘y,where symbol “∘” in equation (25) presents the Hadamard product. b is an adjustment factor of the search radius, and I is an adaptive factor. y is a unit column vector with cn dimensions. All of the above variables are calculated using equations (26)–(28).(26)b=arctanTa−0.5ma,Ta≥1,0.05,Ta<1,(27)Ii=1γi<β0γi≥β,i=1,2,…,cn,(28)y=SYSY,where Ta is the current temperature. ma and γi are random numbers generated randomly in section 0,5 and 0,1, respectively, based on the uniform distribution. β is an adaptive probability and takes 0.1 generally. SY is a column vector with cn dimensions generated randomly based on the standard normal distribution.New solutionC generated by equation (25) may not be in search space S, and, in this situation, each element ci (i=1,2,…,cn) of new solution C is amended by(29)ci=ciU−modciL−ci,ciU−ciL,ci<ciL,ci,ciL≤ci≤ciU,ciL+modci−ciU,ciU−ciL,ci>ciU.


### 5.2.2. MDSA Algorithm Structure

In the MDSA algorithm, the initial temperatureT0 satisfies a requirement that new solutions are accepted in the early stage with a probability greater than the desired value (between 0.8 and 0.9). The refresh method of Ta refers to the paper by Ali et al. [27]. The length of Markov chain La is calculated by equation (30), which is the maximum number of inner iterations.(30)La=10cn+10cn1−eZmin−Zmax.The acceptance probability of new solutionC is determined by the following equation:(31)Pacc=1,ZC≥Zmin,eZC−Zmin/Ta,ZC<Zmin.LetΔZmax=Zmaxa+1−Zmaxa.δand τ are given positive real numbers. If ΔZmax≤δ,for successiveΔacooling timesΔa<τ, which means Zmax does not change or significantly change after a certain cooling, the MDSA algorithm is judged to enter a local extremum. The structure of the MDSA algorithm is shown in Figure 3.Figure 3
The structure of the MDSA algorithm.In the early stage (a≤50) of MDSA algorithm, new solutions are generated mainly by GNS 1, which can expand the search range and improve the coverage of solution set A in search space S. In the later stage (a>50) of MDSA algorithm, new solutions are generated mainly by GNS 2, which can improve the convergence efficiency and stability. GNS 3 can avoid the MDSA algorithm falling into a local extreme too early.


## 5.2.1. Solution Generation

Three methods of generating new solutions are included in the MDSA algorithm. The first one generates a new solution randomly in a search space based on Cauchy distribution with a variable scale parameter [28], called GNS 1. The second one generates a new solution by selecting a certain number of solutions randomly from the solution set and combining them with the current optimized solution [27], called GNS 2. The third one generates a new solution by searching randomly around the current optimized solution based on the standard normal distribution [28], called GNS 3.For ease of description, the number of decision variables is denoted ascn. C is the column vector of the decision variables, which is also called a solution vector. CL is the column vector of the lower limits of the decision variables. CU is the column vector of the upper limits of the decision variables. Then,(19)CL≤C≤CU.All the solution vectors satisfying equation (19) form a search space of solutions, denoted as S. Based on the constraints in the optimization model, the solution vectors in S are feasible solutions. ZC is used as the evaluation function in the MDSA algorithm. Let N=7cn+1. N solutions are generated randomly in search space S and form an initial solution set A. In set A, the best solution is CBest and the worst solution is CWorst. The evaluation function values are denoted as Zmax and Zmin, respectively.Three methods of generating new solutions are described in detail below.(1) GNS 1. A new solution C is generated randomly in search space S based on Cauchy distribution with a variable scale parameter. The Cauchy distribution probability density function centered at the origin is shown as follows: (20)fξ=1πφφ2+ξ2,−∞<ξ<∞,where φ is a variable scale parameter determined by equation (21) based on the number of cooling iterations in the MDSA algorithm, which leads to searching in a wide range early and on a small scale later in the optimization process.(21)φ=1.1a,a≤50,0.9a−50,a>5,where a is the number of cooling iterations.New solutionC is generated by the following:(22)C=CL+CU−CL∘absy,w∈0,0.33,CL+CU2+CU−CL2∘y,w∈0.33,0.67,CU−CU−CL∘absy,w∈0.67,1,(23)y=AYAY,where symbol “∘” in equation (22) represents the Hadamard product. AY is a column vector with cn dimensions generated randomly by equations (20) and (21), and y is a unit column vector with cn dimensions. abs∗ represents the absolute value function. w is a random number generated randomly in section 0,1 based on the uniform distribution.(2) GNS 2. RN solutions are selected randomly from the solution set, denoted as Cn1,Cn2,…,CnRN. Then, these solutions are combined randomly with the current optimized solution CBest. Each element ci (i=1,2,…,cn) in new solution C is generated by the following:(24)ci=cink,wi∈k−1RN+1,kRN+1,k=1,2,…,RN,ciBest,wi∈RNRN+1,1,where wi is a random number generated randomly in section 0,1 based on the uniform distribution. RN makes up a proportion of the number of solutions in the solution set. This method helps greatly in improving the convergence efficiency and stability.(3) GNS 3. A new solution C is generated by searching randomly around the current optimized solution CBest based on the standard normal distribution as follows:(25)C=CBest+bI∘CU−CL∘y,where symbol “∘” in equation (25) presents the Hadamard product. b is an adjustment factor of the search radius, and I is an adaptive factor. y is a unit column vector with cn dimensions. All of the above variables are calculated using equations (26)–(28).(26)b=arctanTa−0.5ma,Ta≥1,0.05,Ta<1,(27)Ii=1γi<β0γi≥β,i=1,2,…,cn,(28)y=SYSY,where Ta is the current temperature. ma and γi are random numbers generated randomly in section 0,5 and 0,1, respectively, based on the uniform distribution. β is an adaptive probability and takes 0.1 generally. SY is a column vector with cn dimensions generated randomly based on the standard normal distribution.New solutionC generated by equation (25) may not be in search space S, and, in this situation, each element ci (i=1,2,…,cn) of new solution C is amended by(29)ci=ciU−modciL−ci,ciU−ciL,ci<ciL,ci,ciL≤ci≤ciU,ciL+modci−ciU,ciU−ciL,ci>ciU.


## 5.2.2. MDSA Algorithm Structure

In the MDSA algorithm, the initial temperatureT0 satisfies a requirement that new solutions are accepted in the early stage with a probability greater than the desired value (between 0.8 and 0.9). The refresh method of Ta refers to the paper by Ali et al. [27]. The length of Markov chain La is calculated by equation (30), which is the maximum number of inner iterations.(30)La=10cn+10cn1−eZmin−Zmax.The acceptance probability of new solutionC is determined by the following equation:(31)Pacc=1,ZC≥Zmin,eZC−Zmin/Ta,ZC<Zmin.LetΔZmax=Zmaxa+1−Zmaxa.δand τ are given positive real numbers. If ΔZmax≤δ,for successiveΔacooling timesΔa<τ, which means Zmax does not change or significantly change after a certain cooling, the MDSA algorithm is judged to enter a local extremum. The structure of the MDSA algorithm is shown in Figure 3.Figure 3
The structure of the MDSA algorithm.In the early stage (a≤50) of MDSA algorithm, new solutions are generated mainly by GNS 1, which can expand the search range and improve the coverage of solution set A in search space S. In the later stage (a>50) of MDSA algorithm, new solutions are generated mainly by GNS 2, which can improve the convergence efficiency and stability. GNS 3 can avoid the MDSA algorithm falling into a local extreme too early.


## 6. Numerical Experiments

In this section, experimental analysis is performed on Wuhan-Shenzhen high-speed railway in China, and price solutions with different elastic demand coefficients are compared.


### 6.1. Sample Setting

In the case study, 51 trains with different departure times and stop patterns are selected. 8 stations, including Wuhan, Yueyangdong, Changshanan, Hengyangdong, Chenzhouxi, Shaoguan, Guangzhounan, and Shenzhenbei (numbered from 0 to 7), are shown in Table3. Only 11 O-D pairs are involved because passenger demand of other O-D pairs is low and insufficient to study. The historical ticket reservation data in the first half of 2016 are used. The fixed price c¯rs0 and the relevant demand qrs0 of each O-D pair r,s are shown in Table 3. The total demand under fixed prices is 41,759 and ticket revenue Z0=1.1775×107. The upper and lower limits of the price of each leg are shown in Table 3.Table 3
Fixed prices and demands and price limits.
O-D0-60-71-21-62-32-42-62-73-65-66-7qrs067403786201034912352252649474565214723236872c¯rs0(¥)463.553869.5373.579.5149.5314388.5244104.574.5crsmin(¥)32437749261561052202721717352crsmax(¥)6036999048610319440850531713697The solution dimensioncn=296 and N=2079 (the number of solutions in set A). In the early stage of the algorithm, La=5920 (the length of the Markov chain). ϵ1=0.05 and ϵ2=0.05 are the convergence precisions. Time value coefficients θ1=0.8¥/min and θ2=1¥/min.


### 6.2. Optimization Results

9 different elastic coefficientsϕ are tested, and 30 experiments are performed for the MDSA algorithm with each elastic coefficient. The ticket revenue (optimization objective) of the optimized solution is denoted as Z∗. The test obtains 30 optimized solutions for each elastic coefficient, and the biggest ticket revenue of these solutions is selected, called the best Z∗. The average ticket revenue of these solutions is also calculated, called the average Z∗. The achieved proportion of the best Z∗ after 200 cooling iterations is denoted as AP.The statistical results ofZ∗ for each elastic coefficient are shown in Table 4, where the last column is the increase in the proportion of the best Z∗ compared with Z0. The relative standard deviation (RSD) of Z∗ is no more than 0.10%, illustrating the stability of the convergence of MDSA algorithm. Over 98% of the best Z∗ has been achieved after 200 cooling iterations, illustrating the efficiency of the convergence of MDSA algorithm. The average number of cooling iterations (ANCI) is less than 900 due to high convergence precisions (ϵ1=0.05 and ϵ2=0.05). The convergence details of the best optimized solutions are shown in Figure 4.Table 4
Optimization results.
ϕBestZ∗AP (%)AverageZ∗ANCIRSD (%)Increase (%)0.214,265,42098.9714,263,0355200.0121.150.413,561,21098.5113,556,5585150.0215.170.612,923,75098.8112,913,5037370.049.760.812,620,75099.2812,615,2825770.027.18112,669,03099.2212,659,6355960.057.601.212,875,18099.2812,870,7594570.029.341.413,234,92098.8513,219,3486710.1012.401.613,707,83098.5613,681,1898590.0916.411.814,132,78098.4914,093,1976590.1020.02Figure 4
Convergence of the MDSA algorithm.According to Table4, when ϕ≥1 (ϕ<1), the best ticket revenue Z∗ increases as the elastic coefficient ϕ becomes greater (smaller), which is also shown in Figure 5. The best Z∗ increases by 7%–21% compared with Z0. Hence, the optimal differential pricing method for multiple trains can improve ticket revenue, and the ticket revenue is greatly influenced by the elastic coefficients.Figure 5
The best ticket revenue (Z∗) of the optimized solutions with different elastic coefficients.For the limited space, the best optimized solutions of some OD pairs withϕ=0.6 and ϕ=1.4 are shown in Table 5. Based on Table 5, the train legs are priced differentially. Most prices are approximate to the upper limits with ϕ=0.6, illustrating that the ticket revenue can be increased by raising the prices in this situation. In contrast, most prices are approximate to the lower limits with ϕ=1.4, illustrating that the ticket revenue can be improved by reducing the prices in this situation. Some train leg prices that are influenced greatly by time-dependent demand are not consistent with this tendency.Table 5
Best optimized solutions (ϕ=0.6/ϕ=1.4).
Departure periodTrain nameLeg0-60-71-21-62-32-42-67G1001603/603699/699———191/173408/3527G1003603/603699/699————407/3278G1005602/325699/39989/63485/261102/56—380/2219G1007429/324416/37789/69365/26191/62106/105274/2209G1009603/327699/378———194/108399/22210G1011603/599699/699————408/37312G1013377/325378/37789/65311/261102/56107/105235/22013G1015352/324377/37790/49314/261101/57105/105225/22114G1017596/333624/37790/49486/263102/56—370/23815G1019603/395699/465————405/27716G1021603/331688/37789/65486/303—194/110408/2556G1101603/484———84/80150/136367/2647G1103603/451—90/49486/39980/72—359/2357G1105602/462—90/49485/39999/57173/126370/2588G1107603/468—90/62486/42089/71172/140370/2739G1109603/494———91/72171/140372/2779G1111603/324—90/70486/262—187/106375/22010G1113576/324—82/62454/261102/66181/107381/22011G1115399/324—85/65346/261—119/105284/22011G1117476/325—90/50—103/99—320/22012G1121582/324—88/64486/262—194/106371/22013G1123601/325———103/79—372/22115G1125593/324————194/110388/22315G1127600/484—————384/28515G1129595/324———99/64—367/22117G1133603/603—89/69486/44484/80168/139374/28218G1135602/600—88/66485/42381/77—371/27719G1139602/603—73/62486/486—187/130390/28819G1141603/502—76/67486/483103/56188/126385/2849G1143602/363———92/73—384/24215G6001——————385/2967G6011—————162/106372/2407G6013—————156/111380/2577G6015—————152/106375/26010G6017————92/69—371/26410G6019——————377/26813G6023——————373/26918G6025————103/56—395/28719G6027————86/76169/126375/2618G6029————81/72—361/23215G6031————101/79194/166369/28917G6033——90/49484/411—155/131399/30013G6111—————194/172371/28617G6117————102/69188/142361/28119G6119—————169/136375/27921G6121—————188/118385/27011G6131——87/67486/47695/77178/170386/2959G77603/603699/699————408/3777G6101————80/61157/111378/2557G6103————80/59-/-369/2557G6105————82/62157/110369/256The elastic demand and unsatisfied demand after passenger assignment of the best optimized solution for eachϕ are shown in Table 6. From Table 6, we can see that the elastic demand increases with greater ϕ except for ϕ=0.2. Passengers are less sensitive to the prices when ϕ decreases, and hence the prices in the optimized solutions increase. The increase of prices results in less passenger demand. However, prices cannot be higher than the upper limits. When ϕ decreases to 0.4 and 0.2, most prices are very close to the upper limits. Thus, there are no big differences for prices with ϕ=0.2 and ϕ=0.4. In this situation, according to equation (11), the elastic demand for ϕ=0.2 is greater than that for ϕ=0.4.Table 6
Elastic demand/unsatisfied demand of the best optimized solution for eachϕ.
O-Dϕ=0.2ϕ=0.4ϕ=0.6ϕ=0.8ϕ=1.0ϕ=1.2ϕ=1.4ϕ=1.6ϕ=1.80-66350/06004/06069/06509/07162/07749/08362/08826/09146/00-73567/03387/03596/03960/04137/04327/04633/4924766/6794974/8501-21900/01796/01727/01738/01984/02287/02503/02475/03180/01-63307/03133/03185/03469/03766/03998/04258/04468/04455/02-32226/02129/02118/02194/02424/02595/02930/02600/02637/02-42411/02297/02362/02456/02609/02763/03058/03234/03629/02-64689/04418/04494/04817/05251/05745/06281/06778/07450/02-74337/04088/04075/04322/04623/04981/05776/06133/06510/03-62025/01911/01951/02119/02268/02446/02678/02827/03031/05-62198/02079/02078/02313/02484/02699/02867/02884/02937/06-76484/06122/05927/05977/06719/07506/08064/08816/08243/0Sum39494/037364/037582/039874/043427/047096/051410/49253807/67956192/850


## 6.1. Sample Setting

In the case study, 51 trains with different departure times and stop patterns are selected. 8 stations, including Wuhan, Yueyangdong, Changshanan, Hengyangdong, Chenzhouxi, Shaoguan, Guangzhounan, and Shenzhenbei (numbered from 0 to 7), are shown in Table3. Only 11 O-D pairs are involved because passenger demand of other O-D pairs is low and insufficient to study. The historical ticket reservation data in the first half of 2016 are used. The fixed price c¯rs0 and the relevant demand qrs0 of each O-D pair r,s are shown in Table 3. The total demand under fixed prices is 41,759 and ticket revenue Z0=1.1775×107. The upper and lower limits of the price of each leg are shown in Table 3.Table 3
Fixed prices and demands and price limits.
O-D0-60-71-21-62-32-42-62-73-65-66-7qrs067403786201034912352252649474565214723236872c¯rs0(¥)463.553869.5373.579.5149.5314388.5244104.574.5crsmin(¥)32437749261561052202721717352crsmax(¥)6036999048610319440850531713697The solution dimensioncn=296 and N=2079 (the number of solutions in set A). In the early stage of the algorithm, La=5920 (the length of the Markov chain). ϵ1=0.05 and ϵ2=0.05 are the convergence precisions. Time value coefficients θ1=0.8¥/min and θ2=1¥/min.


## 6.2. Optimization Results

9 different elastic coefficientsϕ are tested, and 30 experiments are performed for the MDSA algorithm with each elastic coefficient. The ticket revenue (optimization objective) of the optimized solution is denoted as Z∗. The test obtains 30 optimized solutions for each elastic coefficient, and the biggest ticket revenue of these solutions is selected, called the best Z∗. The average ticket revenue of these solutions is also calculated, called the average Z∗. The achieved proportion of the best Z∗ after 200 cooling iterations is denoted as AP.The statistical results ofZ∗ for each elastic coefficient are shown in Table 4, where the last column is the increase in the proportion of the best Z∗ compared with Z0. The relative standard deviation (RSD) of Z∗ is no more than 0.10%, illustrating the stability of the convergence of MDSA algorithm. Over 98% of the best Z∗ has been achieved after 200 cooling iterations, illustrating the efficiency of the convergence of MDSA algorithm. The average number of cooling iterations (ANCI) is less than 900 due to high convergence precisions (ϵ1=0.05 and ϵ2=0.05). The convergence details of the best optimized solutions are shown in Figure 4.Table 4
Optimization results.
ϕBestZ∗AP (%)AverageZ∗ANCIRSD (%)Increase (%)0.214,265,42098.9714,263,0355200.0121.150.413,561,21098.5113,556,5585150.0215.170.612,923,75098.8112,913,5037370.049.760.812,620,75099.2812,615,2825770.027.18112,669,03099.2212,659,6355960.057.601.212,875,18099.2812,870,7594570.029.341.413,234,92098.8513,219,3486710.1012.401.613,707,83098.5613,681,1898590.0916.411.814,132,78098.4914,093,1976590.1020.02Figure 4
Convergence of the MDSA algorithm.According to Table4, when ϕ≥1 (ϕ<1), the best ticket revenue Z∗ increases as the elastic coefficient ϕ becomes greater (smaller), which is also shown in Figure 5. The best Z∗ increases by 7%–21% compared with Z0. Hence, the optimal differential pricing method for multiple trains can improve ticket revenue, and the ticket revenue is greatly influenced by the elastic coefficients.Figure 5
The best ticket revenue (Z∗) of the optimized solutions with different elastic coefficients.For the limited space, the best optimized solutions of some OD pairs withϕ=0.6 and ϕ=1.4 are shown in Table 5. Based on Table 5, the train legs are priced differentially. Most prices are approximate to the upper limits with ϕ=0.6, illustrating that the ticket revenue can be increased by raising the prices in this situation. In contrast, most prices are approximate to the lower limits with ϕ=1.4, illustrating that the ticket revenue can be improved by reducing the prices in this situation. Some train leg prices that are influenced greatly by time-dependent demand are not consistent with this tendency.Table 5
Best optimized solutions (ϕ=0.6/ϕ=1.4).
Departure periodTrain nameLeg0-60-71-21-62-32-42-67G1001603/603699/699———191/173408/3527G1003603/603699/699————407/3278G1005602/325699/39989/63485/261102/56—380/2219G1007429/324416/37789/69365/26191/62106/105274/2209G1009603/327699/378———194/108399/22210G1011603/599699/699————408/37312G1013377/325378/37789/65311/261102/56107/105235/22013G1015352/324377/37790/49314/261101/57105/105225/22114G1017596/333624/37790/49486/263102/56—370/23815G1019603/395699/465————405/27716G1021603/331688/37789/65486/303—194/110408/2556G1101603/484———84/80150/136367/2647G1103603/451—90/49486/39980/72—359/2357G1105602/462—90/49485/39999/57173/126370/2588G1107603/468—90/62486/42089/71172/140370/2739G1109603/494———91/72171/140372/2779G1111603/324—90/70486/262—187/106375/22010G1113576/324—82/62454/261102/66181/107381/22011G1115399/324—85/65346/261—119/105284/22011G1117476/325—90/50—103/99—320/22012G1121582/324—88/64486/262—194/106371/22013G1123601/325———103/79—372/22115G1125593/324————194/110388/22315G1127600/484—————384/28515G1129595/324———99/64—367/22117G1133603/603—89/69486/44484/80168/139374/28218G1135602/600—88/66485/42381/77—371/27719G1139602/603—73/62486/486—187/130390/28819G1141603/502—76/67486/483103/56188/126385/2849G1143602/363———92/73—384/24215G6001——————385/2967G6011—————162/106372/2407G6013—————156/111380/2577G6015—————152/106375/26010G6017————92/69—371/26410G6019——————377/26813G6023——————373/26918G6025————103/56—395/28719G6027————86/76169/126375/2618G6029————81/72—361/23215G6031————101/79194/166369/28917G6033——90/49484/411—155/131399/30013G6111—————194/172371/28617G6117————102/69188/142361/28119G6119—————169/136375/27921G6121—————188/118385/27011G6131——87/67486/47695/77178/170386/2959G77603/603699/699————408/3777G6101————80/61157/111378/2557G6103————80/59-/-369/2557G6105————82/62157/110369/256The elastic demand and unsatisfied demand after passenger assignment of the best optimized solution for eachϕ are shown in Table 6. From Table 6, we can see that the elastic demand increases with greater ϕ except for ϕ=0.2. Passengers are less sensitive to the prices when ϕ decreases, and hence the prices in the optimized solutions increase. The increase of prices results in less passenger demand. However, prices cannot be higher than the upper limits. When ϕ decreases to 0.4 and 0.2, most prices are very close to the upper limits. Thus, there are no big differences for prices with ϕ=0.2 and ϕ=0.4. In this situation, according to equation (11), the elastic demand for ϕ=0.2 is greater than that for ϕ=0.4.Table 6
Elastic demand/unsatisfied demand of the best optimized solution for eachϕ.
O-Dϕ=0.2ϕ=0.4ϕ=0.6ϕ=0.8ϕ=1.0ϕ=1.2ϕ=1.4ϕ=1.6ϕ=1.80-66350/06004/06069/06509/07162/07749/08362/08826/09146/00-73567/03387/03596/03960/04137/04327/04633/4924766/6794974/8501-21900/01796/01727/01738/01984/02287/02503/02475/03180/01-63307/03133/03185/03469/03766/03998/04258/04468/04455/02-32226/02129/02118/02194/02424/02595/02930/02600/02637/02-42411/02297/02362/02456/02609/02763/03058/03234/03629/02-64689/04418/04494/04817/05251/05745/06281/06778/07450/02-74337/04088/04075/04322/04623/04981/05776/06133/06510/03-62025/01911/01951/02119/02268/02446/02678/02827/03031/05-62198/02079/02078/02313/02484/02699/02867/02884/02937/06-76484/06122/05927/05977/06719/07506/08064/08816/08243/0Sum39494/037364/037582/039874/043427/047096/051410/49253807/67956192/850


## 7. Conclusion and Future Research

An optimal differential pricing method for intercity high-speed railway services is proposed in this paper. To analyze passenger choice behaviors between trains, a time-dependent demand formulation is determined based on the historical ticket reservation data, and a passenger assignment method considering departure periods and capacity constraints is constructed by simulating the ticket-booking process. Combining these methods, an optimization model is built with the aim of maximizing the ticket revenue and the decision variables for pricing train legs.A modified direct search simulated annealing algorithm is designed to solve the optimization model, and three random generation methods of new solutions are designed to search the solution space efficiently. The passenger assignment method is used to determine the number of passengers on each train and then to evaluate the price solutions.A case study containing dozens of trains is performed on Wuhan-Shenzhen high-speed railway in China. The results illustrate the stability and efficiency of the convergence of the algorithm. The ticket revenue which can be improved greatly by pricing train legs differently increases by 7%–21% for different elastic demand coefficients.The optimization method in this paper can be extended to the optimal pricing problem in the sale horizon. Then, there are two dimensions of passenger choices: the first one being to choose the ticket-booking time and the second one being to make a selection from multiple trains. To solve this problem, data about the demand distribution in the sale horizon are needed to analyze passenger choice behaviors, and a more efficient algorithm is crucial for this problem considering thousands of decision variables.



---

*Source: 8420206-2021-11-15.xml*