Multi-Objective Stochastic Synchronous Timetable Optimization Model Based on a Chance-Constrained Programming Method Combined with Augmented Epsilon Constraint Algorithm

e design of the timetable is essential to improve the service quality of the public transport system. A lot of random factors in the actual operation environment will aect the implementation of the synchronous timetable, and adjusting timetables to improve synchronization will break the order of normal service and increase the cost of operation. A multi-objective bus timetable optimization problem is characterized by considering the randomness of vehicle travel time and passenger transfer demand. A multi-objective optimization model is proposed, aiming at minimizing the total waiting time of passengers in the whole bus network and the inconsistency between the timetable after synchronous optimization and the original timetable. rough large sample analysis, it is found that the random variables in the model obey normal distribution, so the stochastic programming problem is transformed into the traditional deterministic programming problem by the chance-constrained programming method. Amodel solvingmethod based on the augmented epsilon-constraint algorithm is designed. Examples show that when the random variables are considered, the proposed algorithm can obtain multiple high-quality Pareto optimal solutions in a short time, which can provide more practical benets for decisionmakers. Ignoring the random inuence will reduce the eectiveness of the schedule optimization scheme. Finally, sensitivity analysis of random variables and constraint condence in the model is made.


Introduction
With the acceleration of social and economic development and urban construction, the travel demand of citizens increases. In order to adapt to this change, the number of routes and stops in the public transport systems is increasing, and the optimization of the public transport network is becoming more and more complicated. erefore, e ectively improving the public transport system's service quality and operational e ciency has become a key measure. In the public transport network planning problem, the bus timetable is used to determine the departure time of each line. e relevance of generating a timetable relies on the fact that inadequate and/or inaccurate timetables confuse the passengers and reinforce the bad image of public transit as a whole [1]. At the same time, the design scheme of bus timetables also a ects the operation planning of vehicle scheduling and crew arrangement. Hence, the design of timetables is essential for maximizing the quality of service of the bus system.
In the bus network of large cities, passengers often have to transfer between di erent routes because they cannot get to the destination directly, making it inevitable to spend time waiting for the bus at the transfer station. For commuters who choose to take the bus daily, on the premise that the transfer station can catch up with the following bus line in time, reducing their transfer waiting time to the maximum extent is of positive signi cance to save commuting time and improve the commuting experience. In order to minimize the transfer waiting time, it is necessary to cooperate between different bus lines so that the bus network has good connectivity and the buses on different lines can complete a seamless connection at the transfer station. However, the public transport system is different from the railway transport system. ere is no definite arrival schedule for each bus. Only the first departure time of each bus and the departure interval time of each bus can be known. e randomness of the arrival time of the bus may lead to a long waiting time for passengers and even miss the opportunity of transfer. A synchronized timetable, with good coordination between buses at transfer nodes, enables passengers to enjoy a smooth transfer service, which is a very important and attractive bus transportation service [2]. e public transit planning process is usually divided into a sequence of five steps: (1) the design of routes, (2) the setting of frequencies, (3) the timetabling, (4) the vehicle scheduling, and (5) the crew scheduling and rostering [3]. e optimization of bus timetables is essential because it is a link between the preceding and the following. As far as the existing research on bus timetable optimization is concerned, most of them are based on the assumption of the operating environment, and there are few studies on the actual operation problems. However, in the actual operating environment, there are a large number of random factors in the public transport system, such as vehicle travel time, stop time, bad weather, and passenger demand, which will affect the implementation of the synchronous timetable (Yu et al. [4]; Yao et al. [5]; Wu et al. [6]; Wu et al. [7]). Although some studies consider the uncertainty of vehicle driving in the design of the timetable, some do not consider the impact of this on operational planning. In the collaborative optimization of bus timetables, most researchers set up the condition of fixed headways without considering the effect of nonfixed headways, and most synchronous optimization methods are only a single objective. Few scholars have carried out multi-objective schedule optimization research. Recently, Tang et al. [8] proposed a robust scheduling strategy for the electric bus given the randomness of bus operating conditions. At the same time, they considered the limitation of the timetable and battery mileage and constructed an optimization model of static and dynamic combinations. e static model introduced the buffer distance strategy to solve the adverse effects caused by randomness. In contrast, the dynamic model used the real-time changing urban road traffic conditions to re-plan the bus schedule in operation within a day.
Besides the stochastic effect of operation conditions, the dynamic change of passengers' demand also increases the difficulty of vehicle synchronized optimization. In order to realize passengers' continuous transfer, bus companies often make buses between different routes arrive at the transfer point synchronously by adjusting the existing timetable. But suppose there is a big deviation between the revised timetable and the original timetable. In that case, it will not only destroy the normal service order of the bus timetable but also significantly increase the operation cost of the bus companies due to the adjustment of the next operation plan. e primary tradeoff faced in the planning and operating tasks is between the level of service faced by the user and the operating costs for agencies [9]. ere is a conflict of interest between synchronizing bus timetables between different routes to minimize passengers' waiting time to transfer and minimize the deviation from the original timetables, which requires a difficult tradeoff between improving synchronization and operating costs.
Based on the above factors, we consider the influence of various random factors and take them as the research focus. At the same time, we consider the balance between passenger convenience and operation cost. In order to minimize the total waiting time of passengers in the whole road network and minimize the inconsistency between the synchronized optimized timetable and the original timetable, a multi-objective optimization method is proposed to solve the problem. is study optimizes the original timetable to improve the synchronization of vehicles arriving at transfer stations on different routes so that more passengers can reduce the waiting time and improve the service experience. In order to achieve this goal, we propose a multi-objective optimization model, which fully considers the randomicity of vehicle travel time in public transport routes as a random variable and finds that the stochastic programming model obeys normal distribution through large sample analysis. erefore, the stochastic programming model is transformed into a new deterministic model using chance-constrained programming. In this study, we design a modelsolving algorithm based on the augmented epsilon-constraint algorithm, which can obtain several Pareto optimal solutions reasonably. Finally, the sensitivity analysis of the random variables and the confidence of constraints in the model are carried out to verify the robustness of the model. e rest of this study is as follows: Section 2 reviews the past work. Section 3 describes the multi-objective model. Section 4 analyzes the solution algorithm. Section 5 validates the model's effectiveness and makes a sensitivity analysis. Section 6 gives the conclusion.

Literature Review
e collaborative design and optimization of bus timetable are to optimize the departure time and arrival time from the planning layer and the operation layer, so as to reduce the waiting time of passengers and the operation cost of the bus company, and thus provide the service quality of the whole bus system.
From the point of view of transfer optimization, Cevallos and Zhao [10] proposed a systematic approach based on a genetic algorithm to optimize transfer time by adjusting existing timetables to find the best feasible solution for the optimization of transfer time, taking into account the randomness of bus arrivals. Shafahi and Khani [11] constructed a mixed-integer programming model intending to minimize the waiting time for passengers to transfer between different lines. e model considered the fixed headways, necessary vehicle stop time, and transfer time. e author also extended the model with other buses arriving as new variables. Parbo et al. [12] developed a heuristic algorithm, which used a tabu search algorithm to reduce the offset from the original timetable to minimize passengers' waiting time.
Abdolmaleki et al. [13] assumed that the bus on each line had fixed headways and proposed a model with the same departure constraint to minimize the vehicle travel time in the bus network. Based on the proof that the local search for the generalized problem is equivalent to the maximum-directed cut problem, the authors proposed an approximate algorithm for solving the problem of maximum-directed cut to solve the timetable synchronization problem. At last, the authors relaxed the assumption of the fixed headways and the endless number of transfer passengers between any two lines and designed a recursive quasi-linear time algorithm to solve the problem. To solve the operation problem of the BRT system with two-way lanes, Seman et al. [14] proposed a B-IHCPS strategy based on the simultaneous integrated control of the headway and sequencing of the directional bipartitions, which aims to minimize the total time for passengers to travel and wait. Based on a bimodal network composed of pedestrian and bus modes, De-Los-Santos et al. [15] constructed a mixed integer linear programming model to minimize the total travel time of users and considered service operator objectives. Finally, two new models of line network design and line planning are proposed.
While reducing waiting time for passengers, it is also necessary to reduce the operating costs of bus companies. Castelli et al. [16] proposed a heuristic algorithm for traffic network scheduling based on the Lagrange relaxation method to minimize the weighted sum of the time spent by passengers in the public transportation system and the cost of vehicles used by operators. Moreover, the model proposed by the authors is more suitable for the local network where passenger demand changes significantly within a day. Wu et al. [17] proposed a new coordination design model of stochastic bus timetables. Passengers can readjust their routes when they missed the bus transfer and established a two-level programming model. In the experiment, the authors found that the combination of vehicle scheduling coordination and dynamic diversion can reduce the vehicle size to save the cost of operation and reduce the loss caused by passengers' missed transfers. Considering that the traditional timetable adjustment is restricted by the dynamic traffic jam on the road, Zhang and Liu [18] adopted the dual dynamic method and used the Macroscopic Fundamental Diagram framework to model the traffic dynamics, thereby helping the operator reduce bus operating costs and net benefit while maintaining costs for passengers. e above work is based on minimizing the waiting time of passengers or minimizing the operation cost. However, the collaborative design of the bus timetable can also be optimized by avoiding bus congestion at common stops and maximizing the number of synchronizations.
Considering the importance of designing the maximum synchronism timetable, Ceder et al. [19] defined synchronization as the simultaneous arrival of vehicles on different routes at the transfer station and then set up a mixed-integer linear programming model to maximize the number of vehicles arriving at the transfer station simultaneously. Based on the constraint of headways, they designed a heuristic algorithm that can solve this problem in polynomial time. en, Ibarra-Rojas and Rios-Solis [20] loosened the definition of Ceder et al. [19], considering that the intervals between two routes are all synchronized in a small time window, and established a mixed-integer programming model considering the uneven headways. e flexible, collaborative design of the public transport network was adopted to maximize the number of vehicles arriving at the transfer station simultaneously and avoid bus congestion. After proving that the bus network timetable collaborative design problem is an NP-hard problem, the authors designed a multi-start iterative local search algorithm that efficiently finds high-quality solutions. In order to effectively combine global coordination and long-term operation, Wang and Sun [21] proposed a multi-agent deep reinforcement learning framework to develop a dynamic and flexible bus route retention control strategy to solve the bus bunching problem. To better explore the best strategy, the authors developed an effective headway-based reward function in the framework and used a joint action tracker. In order to train each bus within the framework, the authors also designed an efficient learning algorithm using approximate strategy optimization. In order to reduce the occurrence of bus congestion, Ma et al. [22] proposed a nonlinear optimal control model considering driving disturbance and passenger demand uncertainty. Considering the uncertainty of the bus system and the real-time control requirement of bus regulation, a robust optimal predictive control algorithm is proposed.
Collaborative optimization of bus timetables often involves multiple objectives, and the conflict of interests among them is unavoidable. erefore, it is necessary to design a multi-objective optimization method to find the balance of interests.
Kwan and Chang [23] studied the synchronization of multi-objective metro timetables, aimed at minimizing the total passenger dissatisfaction index related to transfer waiting time and the total deviation index deviating from the original timetable. e Pareto frontier between these two indexes was searched by the nondominated solution sorting genetic algorithm-II (NSGA-II), and it was improved by the differential evolution of NSGA-II. In order to improve the service quality of the bus network, Hassold and Ceder [24] studied the dual-objective timetable design problem to minimize the expected waiting time of passengers and minimize the deviation from the deviation expected vehicle seating rate. By increasing the condition of using multiple vehicle types on the same line to avoid the constraints of uniform headways and passenger load. Assuming that the possible departure time is known, they designed a multiobjective label correction algorithm to solve the multi-objective problem. Wu et al. [25] studied the multi-objective public transport network design and frequency setting problem to minimize the cost of passengers and operators and proposed an alternate objective genetic algorithm to solve the problem with large search space and multiple constraints. Liu et al. [26] developed a dual-objective integer programming model, which aims at maximizing the number of vehicles arriving at the transfer station and minimizing the required vehicle size, and designed a sequential search method solving model based on the two-stage deficit Mathematical Problems in Engineering function. Fonseca et al. [27] proposed a mathematical method to solve the two-objective mixed integer programming problem of timetable coordination and vehicle scheduling. is method allows the timetable to be modified within the departure distance limit to minimize the weighted sum of passenger transfer time cost and bus operation cos. e timetable can be coordinated by modifying the departure time and prolonging the vehicle stay time at the transit station. Tang et al. [28] constructed a bi-objective bus timetable optimization model based on a data-driven method to minimize the total waiting time of passengers and the departure time of bus companies and designed an improved nondominated sorting genetic algorithm-II (NSGA-II) with the special coding scheme, which can search the Pareto optimal solution quickly. Wu et al. [29] proposed a mixed integer linear programming model to optimize the bus schedule and bus matching scheme, which aims at minimizing passenger waiting time, passenger travel time, and operating cost. In order to solve the model, the author designs an algorithm based on improved Lagrangian relaxation and proposes a rolling level scheme for dynamic algorithm implementation. Tian et al. [30] proposed a new three-level programming model to solve the problem of public transport network design with congestion, which minimizes the total operating cost and passenger transfer cost. In order to solve this problem, two methods are proposed to transform the nonlinear nonconvex problem into mixed integer linear programming and surrogate optimization.
ere are many uncertain factors in the actual environment of bus operation, such as the random travel time of vehicles, the dynamic demand of passengers, and so on. In an uncertain environment, the above bus schedule optimization model based on the deterministic operation environment is challenging to solve such problems. In order to solve the problem of uncertainty, researchers consider adding random variables to the construction of a model.
Ting and Schonfeld [31] studied the problem of minimizing the total operating costs of multi-hub transportation networks and designed a heuristic algorithm to simultaneously optimize the headways and timetable relaxation time of transfer lines. e results show that the relaxation delay costs of vehicles and passengers increase with the relaxation time, while the costs of passenger missed transfer and vehicle scheduling delay decrease. Considering the stochastic disturbance caused by the change in passenger demand, Yan et al. [32] set up a stochastic demand scheduling model and developed two heuristic algorithms to solve the model using simulation technology and routing strategy-based online network and designed a performance evaluation method. Zhao and Zeng [33] defined a passenger cost function based on the random arrival time of passengers, line network, vehicle spacing, and timetable. ey proposed a meta-heuristic optimization method of bus line network combined with simulated annealing, taboo, and greedy search to determine a bus line network that minimizes the passenger cost function. Yan et al. [34] proposed a bus network design problem considering the randomness of vehicle driving time, constructed a robust optimization model aiming at minimizing the weighted sum of operator cost, and designed a heuristic algorithm based on the k-shortest path algorithm, simulated annealing algorithm, Monte Carlo simulation, and probabilistic discrete selection model to solve the model. Wu et al. [6] constructed a random integer programming model for the three types of passengers, namely, transfer passengers, passenger passengers, and direct passengers, to minimize the total waiting time.
e relaxation time was added to the timetable to reduce the randomness of bus travel time, and a genetic algorithm with local search was designed to solve the model. Berrebi et al. [35] regarded the bus scheduling problem as a random decision-making process. Assuming that there are random operating conditions and unstable driving distances, they proposed a strategy for scheduling buses using real-time information so that passengers can minimize the waiting time at the transfer station and maximize the bus frequency on the route. Li et al. [36] took the bus travel time as a fuzzy variable, constructed a bi-objective optimization model to minimize the total vehicle travel time and the total passenger waiting time in the bus network, and designed a genetic algorithm with variable chromosome length to solve the model. To cope with the impact of randomness, Morales et al. [37] proposed a bus injection bus operation strategy, namely the bus scheduling strategy for the situation of extremely long headway. e authors established a random model based on the second moment of interval distribution to determine whether the bus is worth injecting and developed a complete service model to determine when the bus should be injected. Based on the above literature research, a multiobjective optimization model is constructed in this study.

Assumptions.
e bus schedule optimization problem studied in this study is to find a reasonable and relatively optimal set of departure times for a certain route in a certain period from the point of view of the transfer station. Let the line set be X, and both line M and line N belong to the line set X. In order to simplify the problem, assuming that line M is the mainline and line N is another line, and line M to line N has transfer demand, it is necessary to optimize the bus departure schedule of line N. To simplify the problem, this study only considers the one-way transfer demand. For a large public transport network in real life, no matter whether the number of bus transfer lines is increased or the request for two-way transfer is considered, we only need to consider the problem by extending the model. Due to data confidentiality reasons, the GPS location information of Shenzhen bus cards is missing. It is no longer feasible to build a model considering the transfer of a single passenger in the traditional model building method. erefore, we designed an innovative method that considered from the overall perspective, ignoring the heterogeneity of stations and stations, and considering the full sample planning problem.
where l i is the drop-off time of the transfer passenger who holds the ith IC card in line M and L M is the total drop-off time of all transfer passengers online M. (c) e total boarding time of transfer passengers online N: assuming the number of trips corresponding to the record of the ith IC card online N is a(i). In general, the number of all trips in line N is counted, and it is assumed that there are A trips in total. Considering the statistical analysis of large samples, we can find out the earliest card swiping time of all card swiping times of each train in the data set, and set Q N min � q 1 , q 2 , . . . , q a , . . . , q A , and set q A and q 1 as the minimum and maximum card swiping time in the set, then we can calculate the original departure time of different trains approximately by the following formula: where f a(i) N represents the departure time of the a(i)th bus online N in the existing timetable. Similarly, based on the consideration of a large sample, it is assumed that the total travel time of line N is T N . And the travel time of line N is approximately equal to the mean of T N , which is approximately substituted by (1/2)T N and recorded as H N . e boarding time of transfer passengers online N is the arrival time of the bus online N, which can be expressed by the sum of the departure time and travel time of the bus number. We can calculate the total boarding time of the transfer passengers online N according to the following formula: where s i is the boarding time of the transfer passenger who holds the ith IC card in line N, and S N is the total boarding time of all transfer passengers in line N. (d) e total waiting time of passengers transferring from line M to line N: from (2) and (3), we can get the passengers' total drop-off time online M and total boarding time online N, respectively. e total waiting time of transfer passengers from line M to line N can be obtained by the difference calculation of the formula: where T N M represents the total waiting time of passengers transferring from line M to line N.

Notations.
e following terms have been defined and used in our model (Table 1).

Objective Function Analysis.
e objective function of bus departure schedule optimization is composed of two parts: the total waiting time of transfer passengers and the total adjustment time of the bus timetable. Passenger cost is determined by the total waiting time from line M to line N. e total bus departure adjustment time is equal to the deviation value between the optimized timetable and the original timetable of each bus.
us, the total departure adjustment time of line M to line N of the bus company during the research period can be calculated based on the following formula: where B represents the total departure adjustment time of line M to line N.

Mathematical Problems in Engineering
According to the research contents of this study, the objective of the optimization is to minimize both the waiting time for transfer passengers and the total bus departure adjustment time as much as possible. e objective of the optimization is to minimize both the objective function F 1 and F 2 , defined as follows: where x a(i) N is the departure time of the a(i)th bus online N in the optimized timetable and H N ′ represents optimized bus travel time.

Model Constraint Analysis
(a) Vehicle Headway Constraint. After a large sample analysis and in light of actual needs, assuming that the maximum departure time interval is G max and the minimum is G min , the vehicle headway constraint can be obtained as follows: we can see that the departure time of the a(i)th bus online N in the existing timetable is f a(i) N . Assuming that based on the existing timetable, the departure time of each vehicle can be fluctuated up and down K ranges. erefore, the adjustable departure time constraints are defined as follows:  6 Mathematical Problems in Engineering Since the departure time of the next bus must be later than that of the previous bus, the departure time sequence constraints are defined as follows: (d) Reasonable Departure Time Constraint. Assume that D early and D late represent the earliest and latest departure times of a bus line, respectively, and that the departure times of each bus cannot be earlier than D early and not later than D late . e reasonable departure time constraint can be obtained as follows:

Consideration and Selection of Random
Variables. e three branches of stochastic programming are the expected value model, chance-constrained programming, and related chance programming. Chance-Constrained programming, which A. Charnes and W. W. Cooper proposed in 1959, is an optimal theory in a certain sense of probability. Considering randomness is more practical and can bring higher solution effect and precision. erefore, we consider the travel time of line N as a random variable, defined as H N ′ . en, the objective function and the corresponding constraint conditions are simplified by chance-constrained programming.

Objective Function
Optimization. First, we define F 1 and F 2 as the optimal values of the two objective functions obtained by solving the deterministic programming through the above modeling process without considering random variables. F 1 and F 2 are the minimum of F 1 and the maximum of F 2 respectively. By counting the travel times of buses on route N from route M to route N on different days, it is clear that H N ′ follows a normal distribution, thus converting it from random variables to deterministic equivalents u N + σ N · Φ − 1 (λ 1 ). Since we consider the general situation and ignore the extreme situation, here we calculate based on the existing data that the random variable obeys the normal distribution, which is consistent with the general reality and is directly simplified to consider the normal distribution. For extreme situations that may exist in reality, such as holidays, etc., passenger transfer needs will change greatly. Based on 3.2 (6), the optimized objective function can be expressed in a joint expression: max λ 1 , Let D � 1 − λ 1 , convert the above union expression to the formula:

Model Constraint Optimization
(a) Constraint of Relationship between Maximum Swipe Time and Departure Time. Considering that the travel time of the buses online N is a random variable H N ′ , we change the model constraint equation (11) to the following equation (17), where α N is the confidence level of constraint satisfaction.

Mathematical Problems in Engineering
Let (N) � c max − x a(i) N − 2H N ′ , then it satisfies the following: ere is a constraint (20) and can also be expressed as equation (21): Let the inverse of the left equation of (20) be η, then η is expressed as equation (22), so n satisfies the standard normal distribution Φ(0, 1), and its probability density is expressed as equation (23): erefore, equation (21) can be transformed in form as shown in the following: According to probability theory, there must be a number K a for the confidence level α N , which makes the following formula (26) true. So the derivation can be done as shown in equations (27)- (30): en, the constraint x a N + H N ′ ≥ c max is equivalent to equation (31), where Φ − 1 (α N ) is an inverse function of α N . Further merging and simplification lead to equation (32).
(b) Variable Constraint. Standardize the range of values of ordinary parameters in the model, the specific meaning of which has been stated in the notation, and therefore is not repeated here but expressed in the formula as follows: 8

Optimized Model.
After considering random variables and changing the objective function and constraint condition, ns, as shown above, the final multi-objective chance-constrained model is proposed as follows:

Algorithm Design
e constrained programming model in this study belongs to the category of multi-objective optimization problems. Because multiple goals often conflict, tradeoffs can only be made between multiple goals, and different tradeoffs can be combined into a set of Pareto optimal solutions. e Pareto solution is defined as if there are any two solutions S 1 and S 2 , S 1 is superior to S 2 for all objectives, we call S 1 dominates S 2 . And if the solution of S 1 is not dominated by other solutions, S 1 is called a nondominated solution, also known as the Pareto solution. And the definition of the weak dominating solution is as follows: if any two solutions S 1 and S 2 have f(S 1 ) < f(S 2 ) for all targets, but at least one objective function has g(S 1 ) < g(S 2 ), S 1 is called weak dominating S 2 , that is, S 1 is the weak dominating solution.
In order to effectively find the Pareto solution to the multi-objective optimization problem, scholars have made a long effort. Epsilon-constraint is a representative constraint processing technique proposed by Takahama and Sakai [38], which can relax the constraint to a certain extent so that the algorithm can get a good solution at the edge of the feasible region. Later, on the basis of the original epsilon-constraint method, Mavrotas, and Florios [39] proposed an augmented epsilon-constraint method to improve the original algorithm, which solves the problem that when some constrained targets are not well constrained by inequalities, the original epsilon-constraint method may match to the weak dominated solution.
erefore, the augmented epsilonconstraint method can be well applied to our proposed twoobjective optimization model.

4.1.
Augmented Epsilon-Constraint Algorithm. Epsilon-Constraint is one of the multi-objective integer programming algorithms with both theoretical and computational attractions. It can loosen the constraint and get a feasible solution with low constraint violation and small objective function. Its basic principle is to retain only one objective in the multi-objective problem and add the remaining objective constraint values to the constraints to obtain the Pareto solution of the selected single objective by adjusting the values of the restricted objective.
To facilitate the algorithm description, we simplify the original problem to min z(x) � z 1 (x), z 2 (x) , x ∈ X.
Here, z 1 (x) represents the objective function D (the probability that the stochastic programming transfer time is less than the deterministic programming transfer time), and z 2 (x) represents the objective function F 2 (the variance between the optimized timetable and the original timetable), and the feasible region X includes all constraints of the optimization model proposed in Section 3.3.3.
For the optimization problem in this study, the model focuses on optimizing the total passenger transfer time and improving the passenger service experience, while minimizing changes in the timetable is less important than the former. erefore, the priority of z 1 (x) is greater than that of Mathematical Problems in Engineering z 2 (x). us, by minimizing z 1 (x), constraining z 2 (x), and setting σ 2 as the constraint value of z 2 (x), the standard epsilon-constraint model of the original problem is e Pareto solution of the original problem can be obtained by solving the above model. When all the inequality constraints of the constrained target are in effect, the original epsilon-constrained algorithm can correctly identify the nondominated solution when all the inequality constraints of the constrained target are in effect. However, it is possible to output the weak-dominated solution when the inequality constraints of the partially constrained target are not in effect [40]. e dominant solution has not reached the optimization effect of the Pareto solution, so we should avoid the weak dominant solution as far as possible. erefore, we need to improve the original epsilon-constraint algorithm.

Customized Augmented Epsilon-Constraint Algorithm.
e augmented epsilon-constraint algorithm is an improved multi-objective analytic class optimization algorithm aiming at the problem that the traditional epsilon-constraint algorithm may not get a practical Pareto solution set. By introducing relaxation variables, the inequality constraints corresponding to each constrained object are normalized as equality constraints. e range of values of the constrained object is divided into equidistant meshes and multiplied by a sufficiently small weight ρ > 0, which is then extended to the original objective function. e augmented epsilon-constraint algorithm ensures that only valid Pareto solutions are searched. Here are the steps of the augmented epsilonconstraint algorithm to solve this scenario in this study.
Step 1.Get the range of values of z 1 (x) and z 2 (x) as [z min 1 , z max 1 ] and [z min 2 , z max 2 ], where z min i and z max i are the minimum and maximum values of the objective function z i (x), respectively. When z 1 (x) obtains the minimum value z min 1 , let x 1 be the optimal solution at this time, z 2 (x) will get the maximum value z max 2 � z 2 (x 1 ) in the feasible region. Similarly, when z 2 (x) obtains the minimum value z min 2 , let x 2 be the optimal solution at this time, z 1 (x) will get the maximum value z max 2 � z 2 (x 1 ) in the feasible region. e two objective functions are optimized for lexicographic order as follows: Step 2. Mesh the range of values of z 2 (x). If we want to get g 2 + 1 Parato solutions, there will be g 2 + 1 grid points (including the grid points corresponding to the minimum and maximum of z 2 (x) ). en, the range of values will be meshed into g 2 grids, with each grid having an interval size of (z max 2 − z min 2 )/g 2 .
Step 3. Determine the constraint level for each z 2 (x) in each grid interval. To facilitate the presentation of the algorithm, we set k 2 as the index of the grid points, numbering from 0, so that the grid points k 2 � 0, 1, . . ., g 2 correspond to the constraint level X12.
Step 4.By introducing the relaxation variable u2, the inequality constraint of constrained targets is changed into an equality constraint. e ratio of u2 to z1 is calculated to eliminate the dimensional effect. en multiplied by the sufficiently small weight p, and extended to the objective function 1, the augmented epsilon-constraint method model is constructed. e final model is as follows: x ∈ X; For each lattice k 2 , a corresponding Pareto solution (z 1 (x * ), z 2 (x * )) and the current optimal solution x * , μ * 2 will be obtained for the above model. If μ * 2 ≥ (z max 2 − z min 2 )/g 2 , that is, the value of the relaxation variable is greater than the size of the grid interval, then the next few grid intervals will not be optimized and the same Pareto solution will still be obtained. e Pareto solution can be searched directly by skipping these covered grid intervals. We make ξ 2 is the jump coefficient of the total risk objective and take ξ 2 � [μ 2 g 2 /6500]. If ξ 2 ≥ 1, then let k 2 � k 2 + ξ 2 , the algorithm continues. If the desired number of Pareto solutions is not reached, increase the value of G. Each Pareto solution is the optimal solution of a multi-objective function under various tradeoffs and can be used as a candidate solution. In different situations, decisionmakers can make decisions based on demand. e algorithm flow chart is shown in Figure 1.

Overview and Visualization Analysis of Examples.
is study selects bus route 338 and bus route 615 in Shenzhen, China, as examples to verify the feasibility and effectiveness of the multi-objective optimization model of a bus timetable. In the actual operation of public transport, generally speaking, vehicles driving up and down independently do not affect each other, and the vehicle routes are used respectively according to the upward and downward directions. Bus route 338 has 40 stops, 95 minutes for the whole upward direction of the line, 101 minutes for the entire downward direction of the line, and Bus route 615 has 62 stops, 168 minutes for the whole upward direction of the line and 153 minutes for the whole downward direction of the line. is study chooses a downward direction bus to study to simplify the problem. Figure 2 shows the route of 338 and 615 buses and the location of each station in the real environment. It can be clearly seen that the two bus routes have a section of overlapping paths. e route's two endpoints are Fuyong Seafood Market Station and Hedong Station. According to Section 3.3.1, we use the mean of the total travel time in one direction of the vehicles in di erent routes to approximate the passenger travel time. To prove the reasonableness of the assumption, we selected 338 and 615 roads with multiple identical sites and all sites evenly distributed on both sides of the line midpoint. As shown in Figure 2, passengers' transfer at any station in the overlap path can be equivalent to transferring at the path's endpoints (namely, Fuyong Seafood Market Station and Hedong Station). erefore, for the convenience of the study, the selected bus network is simpli ed. e simpli ed diagram is shown in Figure 3. e arrows on each line in the diagram indicate the direction in which the bus will run.
In this study, the data include: passenger's IC card number, passenger's card transaction date and time, passenger's bus route, and vehicle's plate number. In the research period selected in this study, there are 23838 data on Bus route 338 and 6798 data on Bus route 615, where 65 pieces of data are generated by passengers swiping their cards from Bus route 338 to Bus route 615. Some of which are shown in Table 2. For passengers who take Bus route 338 but transfer to another route, use the same optimization method. Each passenger has a separate IC card number, so the IC card number can represent an independent passenger. e same IC card number record can be used to determine whether a passenger has a transfer demand. For the data shown in date and time of transaction, the rst eight digits represent the speci c transaction date, and the last six digits represent the speci c transaction time (in 24-hour format). Each bus has a unique license plate number, through the license The value range of the total time variation target is divided into several grids. For each grid, a Pareto solution will be obtained, and for all grids, the final Pareto solution set will be obtained plate number can be obtained in the number of bus routes, and each bus in the study period of the number of vehicles.

Summary and Presentation of Calculation Results.
is model is nonconvex optimization. e augmented epsilon-constraint algorithm designed for solving the model is programmed by Python 3.7, and the optimization model is solved by calling Gurobi 9.12. All calculations pass on PCs with CPU Intel Core i5-9300H 2.40 GHz and RAM 8.00 GB, with all parameters being the default, except setting the Time Limit to a maximum of 300 seconds.
In the above initial optimization model, vehicle travel time is stochastic, so this study takes bus travel time as a random variable, introduces the chance constraint based on the scene, and deals with the uncertainty of the random variable by adjusting the con dence level. We transform the  stochastic programming problem into a deterministic programming problem using chance-constrained programming. is study assumes that the travel time of each bus obeys the normal distribution N(u N , σ 2 N ), where (u N , σ N ) are the mean and variance of the travel time of each bus. We transform the objective F 1 of the initial model from solving the minimum passenger transfer waiting time without considering the random variable to solving the probability λ 1 that the random planning transfer time after considering the random variable is less than the deterministic planning transfer time. For ease of presentation, we make D the new deterministic programming objective function, where D 1 − λ 1 , which means that D is the complementary probability of λ 1 . In addition, for the convenience of calculation, we transform the objective F 2 of the initial model from the sum of the absolute values of the optimized timetable and the original timetable o set to the sum of the squares of the timetable o set, which can be e ectively substituted.
In this calculation example, Bus No.615 is selected as the route for passengers to take after transfer, and the departure timetable of 25 trains on this route during the study period (i.e., 6: 00 to 21: 30 during the operation period of the route on December 23, 2016) is optimized. e Pareto approximate optimal solution of the multi-objective optimization problem is obtained by test example. e target values of Pareto nondominated front and Pareto solutions are shown in Figure 4 and Table 3, respectively.
It can be seen from Figure 4 that the complementarity probability of the transfer time of stochastic programming less than that of deterministic programming is negatively correlated with the o set between the optimized timetable and the original timetable. It can be concluded that minimizing passenger transfer waiting time and minimizing timetable adjustment o set are two con icting objectives when vehicle travel time is considered as a random variable.
is is in line with the basic requirement that the Pareto optimal solution can be obtained only when the interests of each objective are con icting in the multi-objective optimization problem. erefore, the optimization objectives selected in this study are reasonable.
From Table 3, we get 10 Pareto approximate optimal solutions by the designed algorithm. Among them, the nondominant solution one and the nondominant solution 10 represent the two preferences, respectively: the minimum total waiting time for passenger transfer and the minimum total o set from the original timetable when considering random variables. ese two Pareto solutions correspond to the optimal adjustment range of the timetable and the optimized departure timetable, as shown in Tables 4 and 5, where di erent shades of color represent di erent timetable o sets. e larger the color depth o set, the smaller the color light o set.
By taking bus travel time as a random variable, the Pareto solution set of D is obtained between 0.40 and 0.49    approximately. When the complementarity probability of the transfer time of stochastic programming less than that of deterministic programming is the largest among all the Pareto solutions, the offset of the optimized timetable is the minimum of 173 min 48 s. With the decrease of complementary probability, the offset of the timetable increases gradually until the complementary probability is the smallest of all Pareto solutions. e offset of timetable reaches the maximum value of 210 min 30 s. When the complementary likelihood is less than the minimum Pareto solution, the offset of the optimized timetable is no longer in the range of the Pareto optimal solution set. Increasing the offset of the timetable will only increase the operation cost of the bus company, but not reduce the complementary probability, and reduce the total waiting time for passengers to transfer. It means that the service quality of the bus system cannot be improved, and the resource of transportation energy may be wasted. erefore, the bus operation department can choose one of the Pareto solutions as the final bus schedule according to operation and management needs. If the bus operation department prefers the minimum waiting time for passenger transfer, it is necessary to make the buses from different routes arrive at the transfer station as soon as possible. At this time, Pareto solution No.1 can be selected as the departure schedule. But more and irregular changes to the original timetable will affect the normal order of public transport services, make the subsequent vehicle scheduling and personnel arrangement more complicated, and lead to an increase in operation cost. If the public transport operation department prefers the minimum operation cost of the enterprise, it is necessary to maintain the original service order as far as possible and adjust the original timetable to the smallest extent or approximately the same extent. At this time, Pareto solution No.10 can be selected as the departure schedule. However, it will increase the waiting time of passengers and reduce the passenger's service experience. If the bus operation department needs to consider both the enterprise operating cost and the passenger's travel experience, it can choose the appropriate result from Pareto solution No.2-No.9 to make the bus timetable.
In addition, it can be seen from Table 4 that for the first vehicle departure time, Pareto solution one and Pareto solution ten both make smaller optimization adjustments at the same time. We find that in the other 8 Pareto solutions, the departure time of the first vehicle is the same as that of Pareto solutions 1 and 10, both of which are 45 s earlier. It can be explained that the change of timetable optimization   Figure 5 shows that during the maximum solution period, the iterative solution process of the multi-objective function reduces with the iteration time, and the optimal target value does not change obviously after the 0.2 s iteration time.
is shows that the convergence of the augmented epsilon-constraint algorithm is good, and the algorithm can nd the approximate optimal solution more stably. is solution can be regarded as a practical solution process.

Analysis and Discussion of Calculation Results.
en, sensitivity analysis method is used to evaluate the e ect of parameter K and con dence α on the objective function. According to the preliminary calculation experiment, the parameter K is taken 30 min as the base, 5 min each time as the step length, the maximum increase 10 min upward and the maximum decrease 20 min downward. e parameter α is taken 0.25 as the base, with each step of 0.01. e maximum increase is 0.05 upward, and the maximum decrease is 0.04 downward. e sensitivity analysis results of minimizing only the objective function D and minimizing only the sum of the departure time adjustments F 2 are shown in Tables 6 and 7 respectively. As can be seen from Figure 6(a), the parameter K has a signi cant e ect on the target function D. e value of the target D decreases monotonously and tends to be stable with the increase of the value of the parameter K. When the value of the parameter K reaches 15 minutes, the change amplitude of the target D caused by the rise of the value of the parameter K becomes slower. At the same time, the total adjustment time F 2 increases monotonously and tends to be stable with the increase of the value of the parameter K. e reason is that the adjustable range of the constraint becomes larger, the constraint becomes more relaxed, the allowable range of uctuations becomes larger, and the total adjustment time F 2 is correspondingly larger.
As shown in Figure 6(b), parameter K has a signi cant e ect on the total value of departure time adjustment F 2 . e total value of departure time adjustment F 2 increases rst and then stabilizes with the increase of parameter K, uctuates slightly in a very small range, and achieves the minimum value when parameter K is 10 min. At the same time, the value of target value D decreases rst and then stabilizes with the increase of parameter K, uctuates slightly in a very small range, and maximizes when parameter K is 10 min.
To sum up, we can consider taking a better value of K between K 10 min and K 15 min to meet the balance of the minimum waiting time for passengers and the minimum total adjustment of departure time.
Empathy analysis: as shown in Figure 7(a), the condence α has a signi cant e ect on the target function D, the value of the target D increases monotonously with the increase of the con dence α, and the sum of the minimum departure time adjustment F 2 decreases monotonously with the rise of the con dence α. From Figure 7(b), it can be seen that con dence α has a signi cant e ect on the sum of departure time adjustment values F 2 , and the sum of departure time adjustment values F 2 decreases monotonously with the increase of con dence α values. In contrast, the value of the target D increases monotonously with the increase of con dence α values. To sum up, the con dence α of the compromise size can be taken into account to achieve the balance of minimizing the values of the two objective functions. e 10 nondominated solutions solved show di erent application scenarios and meanings. Considering the meaning of the value of the objective function, the value of the objective function F 1 , that is, the value of D, means the complementary probability that the sum of the waiting time of the transfer of passengers is smaller than that of the sum of the waiting time of the original deterministic programming model. In other words, it is a measure of the sum of waiting time for the transfer of passengers. e smaller the value of D, the greater the probability that passengers needless waiting time. For objective function F 2 , that is, the sum of the adjustment of departure time, in order to avoid the consumption of human and material resources, the adjustment of time deviation should be as small as possible.
us, the two objective functions are considered from the passenger's point of view and the bus company's point of view. In real life, if the number of passengers is huge, we can consider the nondominated solution which is in the front of the 10 nondominated solutions; If the network of public transport routes is especially complex, we can consider the nondominated solution which is in the back of the 10 nondominated solutions; if the number of passengers or the network of public transport routes is neither large nor complex, we can consider the non-dominated solution      which is in the middle of the 10 nondominated solutions, so that we can get the equilibrium state of two objective functions.

Conclusion
e optimization of the bus departure timetable is one of the key problems in traffic planning. Based on this, this study studies the cooperative optimization of multi-objective public transport network timetables. A deterministic multiobjective programming model is proposed considering the total waiting time and the total departure time, and the randomicity of bus travel time is fully taken into account. We optimize the preliminary proposed deterministic multiobjective programming model and simplify the model using the normal distribution property of travel time random variables.
us, the chance-constrained programming method transforms random variables into deterministic variables. Finally, stochastic programming is transformed into a new deterministic multi-objective programming problem.
is study designs a model-solving algorithm based on the augmented epsilon-constraint algorithm. Numerical results show that the algorithm can obtain a high-quality Pareto solution in a short time. For the bus operation system, the bus operation plan is crucial, which not only controls the basic operation of all buses but also directly determines the efficiency and service level of the whole bus operation system. e main objective of this study is to minimize the waiting time of transfer passengers. Using the model established in this study, we can use the cumulative method to solve the problem of considering reverse bus routes and more bus routes in real life. Moreover, reducing the waiting time of passengers can optimize the passengers' experience and improve the passengers' satisfaction, thus increasing the passengers' willingness to take the bus and improving the revenue of bus companies. On the other hand, in the algorithm proposed in this study, another objective function is constrained so that the adjustment range of the timetable is not too large, which reduces the impact of the timetable adjustment on the bus operation to a certain extent. At the same time, the support guarantee of Pareto makes the model more in line with the actual situation in real life. In addition to using the cumulative method to expand the model, we can optimize the timetable of the entire bus network, the model proposed in this study can have many other applications in real life. In real life, the problems of connecting buses to subways, buses to trains, and buses to planes can be solved using our proposed model. Furthermore, the problem of transfer between subways, the problem of transfer between trains, etc. can also be solved using the model proposed in this study.
In the following work, we will analyze the performance of the augmented epsilon algorithm and the multi-objective programming model using a real large-scale example. We consider the stochastic characteristics of road section intervals congestion planning. However, due to computing power and data limitation, we ignore the heterogeneity between stations and simplify it to consider the stochastic programming with full samples. Due to the complexity of regional bus scheduling and the limitations of research conditions, we will introduce the actual bus GPS data and propose a more accurate model. In the future, we will also combine advanced artificial intelligence algorithms such as machine learning. Using the approximate optimal solution solved by it, the solution solved in this study can be further optimized, so that a more excellent optimization effect can be obtained.

Data Availability
All data, models, and code generated or used during the study appear in the submitted article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.