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.

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 Section 3, 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. where is the acceleration of train operation, is the speed of train operation, is the time of train operation, is the maximum unit traction of the train, is the coefficient of the traction output ratio, and 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:where is the optimal maximum speed of the train in section ; it is not equal to the maximum speed limit of the train in section ; the maximum speed limit (the maximum speed limit of subway trains is equal to 80 km/h in China); 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].More simplified formulations can be derived from (3) as follows: there is , because, in practice, is much smaller than and , when the train runs at the accelerating phase in a section [24]. Therefore, (2) and (4) can be reformulated into (5), given .

In (4) and (5), is the running distance of the train at the accelerating phase in section . We know that , where is the average running speed of the train in section ; is the quantity of the state that satisfies . It refers to the time period when the speed of the train increases from to . After some calculation and simplification, (5) can be further developed into

Owing to the theoretical derivation of the mathematical formula, we assume the correlation coefficient between and , to be

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 represents the average operating speed while 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.

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 by , the headway between successive trains in the same direction by , the operation time of a train in a full day by , and the running time of a train in section by . Also, we set as the dwell time of a train at station , the passenger carrying capacity of a train as , the occupied rate of passenger carrying capacity as , the total length of the Metro line as , and, finally, the length of section as . 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 indicates different ticket prices, where to the number of sections that passengers pass through.

Assuming that the number of passengers arriving at station during the time period is , 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 ; we can express it by (11).

As the headway between two adjacent trains is , the waiting time of most passengers would be less than , and the number of passengers waiting for a specific approaching train at station is defined as . 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 . We let the number of passengers traveling from station to station to be and the number of passengers traveling from station to station to be . The variation in the number of passengers in the train at station , as illustrated in Figure 2, can be calculated by (12).

If the energy consumption of a train in section is and the total energy consumption of the train from the original station to the terminal station is , then the relationship between and can be expressed in (13). In this model, the power conversion rate is denoted by , and the electricity rate is denoted by .

The total number of passengers orientated from station in a day is . Passengers’ value of time is , 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:where is the generalized travel expenses of the Metro; is the travel time, which is defined as the beginning of passenger arrive at station to the end of departure (the passenger arrive at station and it is the time that the train arrives at the station ), including the waiting time and running time; is the cost (ticket fares) for passengers; is a set of random variables that related to the travel cost, where , and 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

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), is the number of passengers traveling from station to station . is running time of passengers from station to station; for example, if , , is equal to , the represents the running time of passengers from station to station. Therefore, 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 express the total waiting time of all passengers at station within . Therefore, 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 , since the ticket fare costs can be expressed by (18), in (18), represents the fare from station to station (see (10)). is a set of random variables related to the travel cost (see (14)). represents the total fare of passengers from station to station. 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).

To 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 is , 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:

(2) The carrying capacity of the trains should meet the following constraint as (21). In (21), represents the train arriving at station and the number of passengers getting on the train minus the number of passengers getting off the train. is the maximum carrying capacity of train, so means the total number of passengers on the train can not exceed the maximum carrying capacity of the train at any station.

(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 , is the terminal velocity, is the Initial velocity, is the acceleration, and is the distance. Assuming that the speed of the following train is and the speed of the preceding train is zero, the maximum braking force of train is , so the distance between the adjacent trains should be satisfied to , ensuring that there is no conflict between the adjacent trains.

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 function into a standard form, and let , so that the objective function can be converted to a minimization objective function .

Step 2. Construct minimized ideal value vector and maximized inverse ideal value vector of two objective functions and while satisfying all of the constraints, as shown in where and are the minimum value of the objective functions and , respectively, while and are the maximum value of the objective functions and .

Step 3. Build membership functions and of two objective functions and , as shown in the following.

Step 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 where is the new objective function combining and ; with is the weight coefficient of the original objective functions and . In (26), is the distance parameter.

Step 5. Change the value of according to their preferences; change the weight coefficients of the original objective functions and ; select different distance models by changing the value of ; the decision makers can specify the value of according to different preferences which are corresponding to different distance calculation methods, as shown in Table 1.

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 individual p are not worse than individual q and there is at least one target of the individual p which is better than that of the individual q, then we say that p dominates q, and it also implies that p is noninferior to q.

(2) Rank and Front. If p dominates q, then the order value of p is lower than q; otherwise, if p and q do not dominate each other, or if p and q are noninferior to each other, then p and q 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 , we know that if p and q are noninferior to each other, then p and q have the same order value; it means p and q in the same front. The crowding distance is used to calculate the distance between the individuals (p and q) 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 and q) 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 = , so the Pareto front =100.

In summary, the framework of the proposed NSGA-II algorithm is shown in Figure 3, 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.

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 Table 2. 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.

Also, the empirical train traction characteristic curve and running characteristic resistance curve corresponding to (27), (28), and (29) are shown in Figure 4 [26].

The daily passenger flow of the metro line is depicted in Figure 5. 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”.

Finally, given the number of sections that the passengers pass through, the price of tickets in Chengdu Metro is listed as follows:

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 Table 4. 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., ), the optimal speed () is 18.3m/s. In that case, the maximum revenue of the Metro company is 491.50 CNY per train, and the average travel costs of the passengers are 9.13 CNY with the distance parameter being equal to 1 or 2.

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., , ; , ), the optimal train operating speed () is approximately 11.5m/s. In that case, the maximum revenue of the Metro company is about 578.60 CNY per train, and the average travel costs of the passengers 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. Figure 6 depicts the Pareto-frontier solutions, in which objective 1 means and objective 2 represents . It presents the Pareto solutions based on NSGA-II algorithm, showing the 100 sets of data of , , and , including two extreme optimization results. It can be concluded that is the best optimal solution from the viewpoint of the company while 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.

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 Table 4. 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 Figure 7, which indicates that operation energy consumption is increasing with the increasing of average speed.

The optimization results of the NSGA-II algorithm, as shown in Figure 8, indicate that the higher operating speed benefits the passengers more, and lower operating speed benefits the metro operating company more.

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 Figure 9. Tables 5 and 6 present the timetable of OTC and OTP.

The optimization results indicate that the operating revenue () in CUT, OTC, and OTP is 536.52 CNY, 544.32 CNY, and 527.73 CNY, respectively. The average travel costs of the passengers () 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 variables in the model, we analyzed the sensitivity of the main parameters (the electricity rate and passengers’ value of time ) affecting train running costs and passengers’ average travel cost in the model. Assuming the sensitivity range of is 0.5-1.5 CNY/(), the step length is 0.1 CNY/(), and the sensitivity range of is 18-30 CNY/h, the step length is 1 CNY/h. Therefore, we get the relationship between the two parameters and decision variables under different weights at , as shown in Figure 10.

The train average operating speed in the sections decreases with the increase of the electricity rate and decreases with the increase of when the electricity rate is greater than 1.3 CNY/(). Also, the train average operating speed in the sections increases with the increase of passengers’ value of time and increases with the decrease of when passengers’ value of time 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.

Appendix

See Table 7.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Additional Points

Highlights. (1) We present a biobjective timetable optimization model to make a balance between the operating revenue of the railway company and average travel cost of passengers. (2) We apply a fuzzy multiobjective optimization and a nondominated sorting genetic algorithm II to solve the optimization problem. (3) We characterize the trade-off between the conflicting objective functions under different types of travel distances.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Key R&D Program [Grant no. 2017YFB1200700] and the National Nature Science Foundation of China [Grants nos. 71871188 and 61503311]. We acknowledge the support of the Science & Technology Department of Sichuan Province [Grant no. 2018JY0567] and the China Scholarship Council. We are grateful for the useful contributions made by our project partners.