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.


s2fieldsofstudy:
Computer Science
Business

text:
A Look Inside 5G Standards to Support Time Synchronization for Smart Manufacturing

Connectivity has a major role in the current transformation of smart manufacturing. 5G is foreseen as an integral part of an end-to-end networking infrastructure supporting smart manufacturing operation. The integration of 5G with time-sensitive networking (TSN)is seen as a holistic communication solution for smart factories. The next generation of industrialization consists of real-time processes, automation, fixed and cellular networks, and high diversity of devices. The integration of fixed and cellular devices in a single infrastructure to enable end-to-end deterministic communications requires high accuracy of time synchronization. The transport technologies in fixed networks have been evolving toward TSN and deterministic networking to fulfill the requirements for industrial communications. The latest cellular specifications for 5G are targeting ultra-reliable and low-latency communications that enable the Industrial Internet of Things. This article describes the state of the art for integrating TSN with 5G networks based on the support given in 3GPP TS 23.501 Release 16. The requirements for time synchronization in factory automation and integration with 5G networks are presented, together with the most recent advancements in the standardization process.

I. INTRODUCTION
HE next generation of industrialization integrates fixed and cellular networking technologies to interconnect devices that require accurate timing information to provide reliable and low latency deterministic communications. The work in fixed networks has been progressing in the IEEE 802.1 Time-Sensitive Networking (TSN) Task Group [1] and the IETF Deterministic Networking (DetNet) Working Group [2]. A core set of TSN specifications are already available and applicable to multiple verticals. The first three generations of cellular networks have been designed mainly for consumer communications and specifically for downlink data transfer, i.e., asymmetric data communication. The latest specifications for 5G are targeting ultra-reliable and low-latency communications (URLLC) that enable Industrial Internet of Things (IIoT).
The industrial automation infrastructure demands reliable transport as well as seamless connectivity of sensors, actuators, and controlling devices distributed in both cellular and fixed networks. To address the needs of next generation industrial networks, high accuracy time synchronization and symmetric transport between fixed and mobile devices is required. This paper reviews the requirements and the state-of-the-art in terms of time synchronization for industrial automation in integrated 5G-TSN networks. This paper describes the TSN capabilities in fixed networks and explains in detail how to provide similar capabilities via cellular networks and devices based on the support given in 3GPP TS 23.501 Release 16 [11].
The paper is structured as follows. Chapter II provides an overview of use cases and requirements for time synchronization in factory automation. Chapter III introduces the solutions defined in TSN for synchronization. The state-ofthe-art in 5G for tackling accurate time synchronization is presented in Chapter IV. Chapter V highlights the main gaps between state-of-the-art and requirements and discusses future works in this area, while Chapter VI concludes this work.

II. REQUIREMENTS FOR TIME-SYNC IN FACTORY AUTOMATION
Industrial applications require not only robust and reliable data transport, but often accurate time information for synchronous communications as well. Over time several technologies have been applied, but Ethernet-based networks are the most widely adopted infrastructure in industrial practice today. Such a proliferation of Ethernet stems from the facts that it is standardized and is a de-facto transport for TCP/IP-based industrial wireline communication, which also facilitates integration of industrial applications with Information Technology (IT) services. The use of Ethernet in industrial automation has been realized via protocols like EtherCAT, Ethernet Powerlink and PROFINET.
Time or clock synchronization is a key requirement of many applications deployed over industrial networks. It aims to harmonize the time among independent, often physically distributed, clocks of network devices. Even if the clocks are initially set to the same time, they drift eventually and accumulate possibly significant time deviations. Many aspects of every-day operation of a network rely on its devices to have a common timescale, e.g., for fault management purposes. There are different industrial applications in the factory automation domain that depend on the clocks to be tightly synchronized.
For example, a system of Programmable Logic Controllers (PLCs) usually has a shared automation objective and needs to implement control actions in a specific sequence. Therefore, the PLCs need to have their clocks synchronized. In such a case, the required time deviation (with respect to a common time T reference) provided up to 100 PLCs should be at most around one microsecond [3]. Another example of stringent time accuracy requirement relates to closed-loop motion control. Motion control governs movement and rotation of machine parts in real-time, by regularly transmitting, e.g., values of requested speed and position to actuators, while also having sensors reporting on the actual values for these parameters.
Both the transmission to all involved actuators and the reporting from all sensors must be accomplished at a specific moment in time to enable precise motion. Building on simultaneous communication under few tens of microseconds, a motion control system requests that up to 100 actuators and sensors have guaranteed clock synchronicity in the order of one microsecond, or even below [3]. High accuracy demand is also put forth by factory deployments of mobile robots, which aid, e.g., transportation of materials or manufacturing of products. These tasks often include a tight cooperation between multiple mobile robot units, for which robot actions need to be strictly synchronized. A summary of time/clock synchronicity demands in industrial automation is given in [3] and in Table 1. These requirements were major input for 3GPP SA2 standardization work when developing time synchronization support provided by 5G [11]. The time synchronicity given in Table 1 is defined as the maximum allowed time offset within a synchronization domain between the master clock and any individual device clock.

III. SYNCHRONIZATION IN WIRELINE NETWORKS AND TSN
To achieve time synchronicity among devices in factory automation, a mechanism is required to deliver proper time from one device to another. This generally comprises the exchange of several messages to enable an estimation of the latency that the network introduces. A device that network nodes synchronize their clocks to, is usually referred as the master device, while others are called slave devices. Many realizations of a time synchronization mechanism can be found today, such as the Network Time Protocol (NTP) [4] , the Precision Time Protocol (PTP) [5], and the generalized PTP (gPTP) [6], which can be considered as specific PTP profile defined in IEEE 802.1AS specifications and requires multicast transport. While NTP is used for wide-area networks such as the Internet, PTP is used, e.g., for telecom networks, gPTP focuses on networks with a limited coverage area and provides better time precision. gPTP may be used on Ethernet-based, optical and microwave networks, and in, for example, WiFi networks. Due to their higher synchronization accuracy, PTP and gPTP are adopted both in mobile backhaul networks and in time-sensitive industrial networks.

A. Time synchronization in mobile backhaul networks
In mobile networks, standard development organizations such as 3GPP and ITU-T play a major part to study and develop synchronization solutions and architectures. 3GPP typically defines timing accuracy requirements and standardizes solutions for mobile air interface. ITU-T also develops solutions and architectures related to 3GPP requirements. Specifically, the ITU-T Recommendation G.8275.1 Precision Time Protocol Telecom Profile [7] specifies a PTP profile that can be used for synchronizing time in wireline mobile backhaul networks. Such mechanism can also help achieving synchronicity in mobile industrial networks, where the timing from the wireline backhaul could be transported and distributed to the industrial network.

B. Time synchronization in TSN
The gPTP mechanism, defined in the IEEE 802.1AS standard [6], is a building block of the TSN standardization efforts carried out by the IEEE 802.1 TSN Task Group (TG). gPTP is explained in [8] and the feature topic on TSN of the IEEE Communications Standards Magazine [9] provides in-depth tutorials of TSN. The initial goal of the TSN TG was to build audio and video bridges within a LAN environment, thus allowing audio and video devices to align in time the streams coming from various sources. This initial use case expanded towards a general addon to traditional Ethernet to offer bounded and low latency, low packet delay variation (jitter) and extremely low packet losses for many use cases, including (but not limited to) industrial automation and automotive in-vehicle communication. Nowadays, the TSN toolbox includes more than 20 standards. Recently, TSN profile specifications (e.g., IEC/IEEE 60802) are being developed for selected use cases, e.g., industrial automation [9].
Industrial automation devices in a TSN network align their activities in time through a shared clock according to the gPTP protocol, thus forming a gPTP domain. The shared clock can either be a universal clock (also known as "wall clock"), or a working clock based on an arbitrary timescale. One node in the gPTP domain acts as the source of the shared clock and it is denoted as GrandMaster (GM). Clock information is distributed from the GM to all the nodes in the domain with dedicated messages defined by the gPTP protocol. To provide an example of clock distribution, consider a typical gPTP domain illustrated in Fig. 1.
The domain includes TSN end stations, one of which is the GM, as well as TSN bridges. In each link between two nodes, the node closer to the GM is called Clock Master (CM) and the other node is called Clock Slave (CS). The gPTP protocol is implemented in a distributed way by each node in the network.
To provide an example, consider the following sequence of operations for Node B (TSN bridge) in Fig. 1: 1. Receive a gPTP Protocol Data Unit (PDU) from its direct CM (Node A) containing the time of the GM and the cumulative time delay from the GM. 2. Compute the propagation delay of the link with its direct CM. 3. Combine the information in 1 and 2 to compute the total delay from the GM and update its internal clock. 4. Compute its residence time, i.e., the time elapsed from the reception of the gPTP PDU (ingress) and its transmission (egress). 5. Send the gPTP PDU to its direct CSs (e.g., Node C), which contains the time delay from the GM (computed in step 3) and the residence time (computed in step 4).
IEC/IEEE 60802 includes time synchronization requirements for both working clock and universal clock. In the first case, the maximum deviation between a node's clock and the GM's working clock must be the range between 100 ns and 1 µs. In addition, redundant working clock domains can be setup for the same set of devices, so that zero failover time can be achieved in case a GM node experiences a failure. In practice, redundant domains can be realized through "cold standby" or "hot standby". In both cases, a node in the gPTP domain is ready to take the role as GM in case of failure of the primary GM. In the "hot standby" case, the back-up GM shares its clock through gPTP PDUs even when the primary GM is active. The current practice suggests the use of "hot standby" backup GMs for achieving the highest synchronization availability.
In case synchronization to a universal clock is implemented, the clock is provided by an external source (e.g., Global Positioning System) to the GM and shared within the plant via gPTP mechanisms. One possible use case for universal time synchronization is the need to record sequence of events by different devices within a plant, to enable root-cause analysis and plant optimization. In this use case, a more relaxed synchronization accuracy is allowed (up to 100 µs) and no redundant clock domain is required.
IV. 5G TIME SYNCHRONIZATION SUPPORTING TSN Integration of 5G with TSN to support time critical industrial applications requires end-to-end time synchronization, too. Time reference information is needed for the applications running in end devices in most of the industrial automation deployments. Time reference information is also required in the bridges of a TSN network when time-based TSN tools are used, like Scheduled Traffic (IEEE 802.1Qbv), to provide deterministic low latency for time-critical traffic.

A. 5G internal system clock
Time synchronization is an integral part of the 5G cellular radio system and is essential for its operation. Time synchronization has been common practice already for cellular networks of different generations. The radio network components themselves are also time synchronized, e.g., for advanced radio transmission, such as synchronized Time Division (TDD) operation, cooperative multipoint (CoMP) transmission, or carrier aggregation. The new capability to provide reference time delivery as a service over the 5G system is designed to exploit the existing synchronized operation of the 5G radio access network. Such a building block enables end-toend time synchronization for the communication of industrial applications running over 5G system. The 5G internal system clock can be made available to a User Equipment (UE) with signaling of time information related to absolute timing of radio frames as described in 3GPP TS 38.331 [10]. 5G NR radio transmission is done in 10 ms radio frames, and each frame is identified with System Frame Number (SFN). These frames are used to define transmission cycles and they are also used to calculate accurate time reference over radio interface to a UE.
The 5G system can use a 5G internal clock as time reference that is distributed over the transport network utilising PTP Telecom Profile [7] to its gNBs (5G base stations) as Fig. 2 shows. Based on the accurate time distributed over transport network, a gNB acquires the reference time value, further it modifies this value considering the propagation time to specific reference points in the system frame structure. Furthermore, a gNB transmits modified time reference value to UEs in either a System Broadcast (SIB) or a Radio Resource Control (RRC) message. It depends on the implementation what frequency is used to update the time reference at the end point connected to the UE, e.g., a TSN end station.  B. 5G system enhancements to support gPTP 3GPP working groups conducted studies on how 5G can support the industrial automation vertical together with TSN including time synchronization. The 5G System (5GS) needs to interwork with the gPTP of the connected TSN network as gPTP is the default time synchronization solution for TSNbased industrial automation [12].
The 5GS is considered as a virtual bridge for each User Plane Function (UPF) [11], which is connected to a TSN bridge. The UPF is a gateway that interconnects the 5G system to external networks, in this case a TSN network. The other end of the 5G virtual TSN bridge is at the UEs. Such a 5G virtual bridge is also considered as a time-aware system as per IEEE 802.1AS. Fig. 3 illustrates the new entities introduced at both edges of the 5G system to provide TSN Translation (TT) functionalities: Network side (NW-TT) at the UPF side and Device Side (DS-TT) at the UE side.
With the selected solution in the standard, the TSN timing information is delivered from NW-TT to DS-TT as gPTP messages via the 5G user plane. This approach is scalable, builds on the established 5G security model, and is flexible in various deployments, including network sharing and the support of non-public networks.
The 5G system is modelled as an IEEE 802.1AS time-aware system for supporting TSN time synchronization as specified in [11]. There are two synchronization processes running in parallel in an integrated 5G-TSN system: a 5G system synchronization process and a TSN synchronization process.
-5GS synchronization: the 5G GM provides the reference clock for 5GS internal synchronization, so the UPF and the UE are synchronized to the 5G reference clock. -TSN domain synchronization: It provides synchronization service to devices in the TSN network. The gPTP Grand Master and the devices that are synchronized are outside the 5G system. Only the TSN translators at the edge of the 5G system track the gPTP domain. That is, only the TTs need to support IEEE 802.1AS operations, e.g., gPTP, timestamping and best master clock algorithm (BMCA).
The two synchronization processes can be considered independent from each other. For example, the gNB only needs to be synchronized to the 5G GM clock, the 5G synchronization process that guarantees RAN function working properly is kept intact and independent from the external gPTP synchronization process.
The independence of the two synchronization processes brings flexibility in time synchronization deployment, e.g., when upgrading an existing deployment. For example, if an operator has already deployed a 5G system in one area, in order to further support factories in the neighborhood with gPTP, only the UPF and UE side need additional enhancement and the overall 5G time synchronization remains untouched. Also, if 5G is added to a fixed TSN network with time synchronization, the TSN time synchronization is not modified due to the introduction of 5G.

C. Downlink time synchronization
For every gPTP PDU entering the 5GS at the UPF, the NW-TT entity generates ingress Timestamps (TSi) based on the 5GS internal clock and embeds the timestamp in the gPTP message. The UPF forwards the gPTP message to the UE via the user plane (i.e., via a PDU session). UE receives the gPTP message and forwards it to the DS-TT. The DS-TT then creates an egress Timestamp (TSe) for the gPTP PDU of the external gPTP working domain. The difference between TSi and TSe is considered as the calculated residence time that the given gPTP message has spent within the 5GS. The DS-TT modifies the TSN timing information received from gPTP messages based on the calculated 5GS residence time and sends it towards to the next time-aware system.

D. Support for multiple working clock domains
In smart manufacturing, multiple gPTP domains can coexist in a single shared communication network infrastructure where each gPTP domain consists of the commonly managed industrial automation devices which can have common timescale [9]. According to the ongoing standardization in [9] the minimum number of domains that a TSN industrial equipment needs to support depends on its device class. A constrained (currently called Class B) device has to support at least two gPTP domains: one working clock domain and one universal time domain. For redundancy, a feature rich (currently called Class A) device has to support at least two of each time domain, i.e., at least four gPTP domains altogether. 3GPP requires the 5G system in Release 16 to support networks with up to 32 working clock domains [3]. A gPTP working clock domain is characterized with a specific domain number in the gPTP PDUs. To support multiple gPTP domains in the 5GS, as shown in Fig. 4, the NW-TT generates an ingress timestamp (TSi) for every gPTP message carrying a specific domain number (where the domain number indicates a specific gPTP domain).  The DS-TT generates TSes for the gPTP PDUs for every external gPTP working domain. Ingress and egress time stamping are based on the 5G system clock to which the NW-TT and DS-TT are synchronized. An end station can select timing information of interest based on the domain number in the gPTP PDU.
A special case that is supported by the 5GS is when the 5GS clock acts as a general grandmaster and provides the time reference not only within the 5GS, but also to the rest of the devices in the deployment, including connected TSN bridges and end stations where, the time domains are merged into one common single time domain.

E. Support for non-public networks and RAN sharing
3GPP defines Non-Public Networks (NPN) for industrial IoT in which the access to the NPN is restricted to a private group of authorized devices. There are two variants of NPN. A standalone NPN is a dedicated and isolated 5G network deployment for a private entity, like a factory. A public network integrated NPN is a private network that is partly sharing communication infrastructure with a public network.
It is also possible to share a common infrastructure for multiple NPNs. Normally, there are many companies in an industry park that provide common facilities to its enterprise customers. Each company can use a dedicated 5G NPN, where the core network nodes are deployed in the company premises and connected to, e.g., its local TSN infrastructure. In such case, a common 5G radio access network (RAN) can be used in the entire industry park in order to provide good coverage and capacity, meanwhile prohibit interference and coexistence challenges of multiple local RANs in a close vicinity. This RAN is shared among the multiple NPNs, but the NPNs are logically separated and capacity reservations can be provided per NPN. Each NPN can have its own working clock domains, which remain isolated, even if a shared RAN is used, as depicted in Fig 5. The proposed solution is not limited in the number of working domains that are supported, which means that each NPN can support multiple working domains independently without limitations of the 5G system.

F. Security aspects of time synchronization
In the 5G time synchronization, the gPTP messages are transmitted as user plane data. Therefore, all defined 5G security mechanisms that are used for data traffic also apply for the time synchronization data. Only authorized devices can participate in the 5G communication and all user data is encrypted. It is not possible for a device to monitor the time synchronization of a factory, even in a case like multiple NPNs sharing a common RAN, unless the device is a member of the NPN and included in the time synchronization communication.
What is available to any device within the coverage of the 5G RAN is the time used in the 5GS. However, this time does not reveal anything about the time that is used by any end device within any NPN. This functionality is commonly available in today's cellular networks, even if the timing is typically on much coarser level (for example, smart phones can receive the time-of-day from the network).

G. Overhead for conveying gPTP messages
Time synchronization is typically applied in real-time industrial control applications like motion control. Let us consider a motion control application that has a user plane data rate of 50 Byte / 500 μs = 800 kbps. The gPTP data rate is around 2 * 50 Byte / 125 ms = 6.4 kbps. The overhead for 5GS to deliver the time sync message in user plane is 0.8%, and the overhead is constant regardless to the number of UEs. The 5GS is very well capable to handle a data rate of 6.4 kb/s for a UE, which traffic is marginal compared to the industrial URLLC traffic. In conclusion, the overhead of gPTP messages is insignificant to provide a time-synchronization service.

V. GAP ANALYSIS AND FUTURE WORK
From the analysis of the state-of-the-art, a few areas emerge in which more work needs to be done. These areas, summarized in the following, will be the subject of future research and standardization efforts.

A. Arbitrary placement of synchronization master
When bringing 5G into an industrial setting, it may be of interest that a timing GM can be connected wirelessly to the system. When a gPTP GM is at UE side, the 5G system needs to support the TSN timing distribution both in uplink direction (UE-to-UPF synchronization) and to other UEs (UE-to-UE synchronization). Although the support for TSN timing distribution in the uplink is only considered for Release 17, it is important to note that the above presented 5G time synchronization solution is well prepared for this functionality, because the gPTP working clock domain is distributed in the user-plane. The UE-to-UE synchronization case introduces more stringent accuracy requirement on the synchronization performance, since it doubles radio propagation delay over the air interface. Further enhancement of 5GS is required in order to fulfill the industrial automation synchronicity requirement for the UE-to-UE synchronization case. Furthermore, the downlink and uplink time synchronization are independent of each other via timestamping and determination of 5GS residence time.
B. Synchronicity budget for 5G system Synchronization for Industry Automation applications, such as motion control, allows maximum 1 µs synchronicity [3]. The requirement is between the sync master and any device of the clock domain, which is an end-to-end (E2E) requirement [13]. A TSN network can have multiple hops in the E2E chain. When 5GS is modelled as a virtual bridge interworking with other TSN nodes, the 5GS only contributes to a fraction of the 1 µs E2E synchronicity budget, therefore, the actual synchronicity budget left for 5GS is even less. This can be partially compensated since the 5GS replaces multiple wired network hops with reliable wireless connectivity, simplifying the deployment in the factory. A reasonable synchronicity budget for 5GS in real industrial automation scenarios needs to be further investigated and activities are planned on 3GPP Rel-17 [14] in this domain.

C. UE relative timing
Once provided with absolute timing from the gNB as 5G reference time, multiple UE devices will be synchronized together, nominally to the same time.
In the 5G radio network the precision is determined by 1. The quantization of the absolute time information provided [10].
2. The difference in the propagation time between each of the UEs and their respective gNB antennas.
Uncertainty in the synchronization of the UEs is introduced by the difference in the delay resulting from the radio signal propagation to the UE, which depends on the relative locations of the UE devices and their respective gNB. If the UE devices happen to be the same distance from their gNB but in different directions, they can still experience the same propagation time of flight delay. Clearly if the UE devices are moving then their relative timing can change.
The differences in the propagation time are largely countered in 5G by the Timing Advance mechanism, a control loop which aligns the timing of the incidence of the uplink signals from the UE as they arrive at the gNB. One major contribution to the error in this mechanism is due to quantization, in 5G of 32 ns or more depending on the numerology of the actual radio configuration. This appears as a source of error in the synchronization between the UEs, whatever the relative separation of the UE devices [15].
The UE radio layer needs to be smarter to distinguish attenuation, packet loss, signal reflections and other propagation delays. Directional antennas with beam forming and other methods to assist the signal to find the best path can be explored for increased timing accuracy.
Taking these two factors into account the relative timing accuracy between two local UE devices connected to the same serving cell should still be within 1 µs.

VI. CONCLUSION
This paper explains how 5G supports time synchronization for industrial automation even in integrated deployments with wireline networks. As future wireline industrial communication networks will be based on TSN, the paper focuses on the integration between 5G and TSN and discusses several aspects that need to be considered for accurate time synchronization between the two domains. If DetNet is also used, then the time synchronization provided by 5G and TSN will apply to DetNet as well because DetNet leverages time synchronization provided by the sub-network technology. We describe how the 5G system internal synchronization is maintained and leveraged, furthermore, what specific enhancements have been added to 5G to support the time synchronization protocol in the wireline domain, i.e., gPTP. Support for multiple clock domains, time synchronization in non-public 5G networks, and security aspects are also explained. The paper also discusses aspects that will be refined by 3GPP Release 17 work, such as the possibility of connecting a gPTP time source via 5G and potentially consider relative synchronization between different UEs.