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

Authors: Yu Yuan; Pengcheng Wang; Minghui Wang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The 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 affect 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. Through 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. A model solving method 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 benefits for decisionmakers. Ignoring the random influence will reduce the effectiveness of the schedule optimization scheme. Finally, sensitivity analysis of random variables and constraint confidence in the model is made.


---

## Body

## 1. 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. Therefore, effectively improving the public transport system’s service quality and operational efficiency 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. The 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 affects 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 different 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 significance 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. There 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. The 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].The 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]. The 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. The 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. The 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]. There 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. This 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. Therefore, the stochastic programming model is transformed into a new deterministic model using chance-constrained programming. In this study, we design a model-solving 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.The rest of this study is as follows: Section2 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.


## 2. Literature Review

The 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. The model considered the fixed headways, necessary vehicle stop time, and transfer time. The 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.The 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. Then, 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. The 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. Therefore, 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. The 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 multi-objective 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 function. Fonseca et al. [27] proposed a mathematical method to solve the two-objective mixed integer programming problem of timetable coordination and vehicle scheduling. This 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. The 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.There 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. The 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. They 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. The 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. The 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 multi-objective optimization model is constructed in this study.


## 3. Model Formulation

### 3.1. Assumptions and Notations

#### 3.1.1. Assumptions

The 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 beX, 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. Therefore, 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.(a)
Passenger transfer demand: considering the same bus IC card, based on the analysis of large sample data, for lineM and line N, if the bus IC card has an adjacent swiping record of line M and line N and the bus IC card has to swipe record of line N after line M, the number of transfer passengers from line M to line N is assumed to be plus 1. We count the number of bus IC cards with these two properties to determine the number of interchange passengers from line M to line N and record them as PMN to represent passenger transfer demand.(b)
The total drop-off time of transfer passengers onlineM: assuming the boarding time of the transfer passenger who holds the ith IC card is ti. All IC cards, namely all transfer passengers, are anonymously processed to avoid losing generality. Considering the statistical analysis of large samples, it is assumed that the total travel time of line M is TM. And the passenger travel time online M can be approximately equal to the mean value of TM, which is approximately substituted by 1/2TM and recorded as HM. Since the passenger’s time to get off is equal to the sum of the passenger’s boarding time and the ride time, the total drop-off time of transfer passengers online M can be calculated by the following formula:(1)li=ti+HM,LM=∑i=1PMNli,whereli is the drop-off time of the transfer passenger who holds the ith IC card in line M and LM is the total drop-off time of all transfer passengers online M.(c)
The total boarding time of transfer passengers onlineN: assuming the number of trips corresponding to the record of the ith IC card online N is ai. 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 QminN=q1,q2,…,qa,…,qA, and set qA and q1 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:(2)t=qA−q1A−1,fNai=q1+a−1t,wherefNai represents the departure time of the aith 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 lineN is TN. And the travel time of line N is approximately equal to the mean of TN, which is approximately substituted by 1/2TN and recorded as HN. The 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:(3)si=fNai+HN,SN=∑i=1PMNsi,wheresi is the boarding time of the transfer passenger who holds the ith IC card in line N, and SN is the total boarding time of all transfer passengers in line N.(d)
The total waiting time of passengers transferring from lineM 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. The total waiting time of transfer passengers from line M to line N can be obtained by the difference calculation of the formula:(4)TMN=SN−LM=∑i=1PMNfNai+HN−∑i=1PMNti+HM,whereTMN represents the total waiting time of passengers transferring from line M to line N.


#### 3.1.2. Notations

The following terms have been defined and used in our model (Table1).Table 1
Compliance description of notations used in this study.
SetXRoutes set,N∈XQminNThe earliest card-brush time set of all vehicles in lineN, qi∈QminNCNaiAll card-brush time of thea(i)th bus in line N, ci∈CNaiParametersATotal number of cars in routenPMNPassenger transfer demandLMTotal drop-off time of all transfer passengers in lineMTMTotal vehicle travel time for lineMTNTotal vehicle travel time for lineNHMMean of total vehicle travel time for lineMHNMean of total vehicle travel time for lineNtiBoarding time of the transfer passenger who holds theith IC cardliDrop-off time of the transfer passenger who holds theith IC card in line MsiBoarding time of the transfer passenger who holds theith IC card in line NaiNumber of trips corresponding to the record of theith IC card in line NqiThe earliest card swiping time of theith IC bus in line NSNTotal boarding time of all transfer passengers in lineNfNaiDeparture time of theaith bus online N in the existing timetable, ∀ai∈1,2,…,ATMNTotal waiting time of passengers transferring from lineM to line NBTotal departure adjustment time of lineM to line NHN′Optimized bus travel time onlineNGmaxMaximum departure time intervalGminMinimum departure time intervalKAdjustable range of departure timetableJFloating range of a vehicle’s travel timeDearlyThe earliest departure times of a bus lineDlateThe latest departure times of a bus linecmaxMaximum card swiping time of theaith bus online NF1¯Optimal value ofF1 derived from deterministic programmingF2¯Optimal value ofF2 derived from deterministic programmingHN′˜Random travel time of vehicles onlineNαNConfidence level of constraint satisfaction for stochastic programminguNMean of random travel time in lineNσNVariance of random travel time in lineNλ1Probability that the random planning transfer time after considering the random variable is less than the deterministic planning transfer timeDecision variablesD1 −λ1xNaiDeparture time of theai-th bus online N in the optimized timetable


### 3.2. Model Analysis

#### 3.2.1. Objective Function Analysis

The 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 lineM to line N. The total bus departure adjustment time is equal to the deviation value between the optimized timetable and the original timetable of each bus. Thus, 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:(5)B=min∑i=1PMNxNai−fNai,where B represents the total departure adjustment time of line M to line N.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. The objective of the optimization is to minimize both the objective functionF1 and F2, defined as follows:(6)minF1=∑i=1PMNxNai+HN′−∑i=1PMNti+HMminF2=min∑i=1PMNxNai−fNai,,where xNai is the departure time of the aith bus online N in the optimized timetable and HN′ represents optimized bus travel time.


#### 3.2.2. 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 Gmax and the minimum is Gmin, the vehicle headway constraint can be obtained as follows:(7)Gmin≤xNai+1−xNai≤Gmax,∀ai∈1,2,…,A,N∈X.(b) Timetable Adjustment Range Constraint. From 3.1.1 (3), we can see that the departure time of the aith bus online N in the existing timetable is fNai. Assuming that based on the existing timetable, the departure time of each vehicle can be fluctuated up and down K ranges. Therefore, the adjustable departure time constraints are defined as follows:(8)fNai−K≤xNai≤fNai+K,∀ai∈1,2,…,A,N∈X.(c) Departure Time Sequence Constraint. 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:(9)xNai≤xNai+1,∀ai∈1,2,…,A,N∈X.(d) Reasonable Departure Time Constraint. Assume that Dearly and Dlate represent the earliest and latest departure times of a bus line, respectively, and that the departure times of each bus cannot be earlier than Dearly and not later than Dlate. The reasonable departure time constraint can be obtained as follows:(10)Dearly≤xNai≤Dlate,∀ai∈1,2,…,A,N∈X.(e) Constraint of Relationship between Maximum Swipe Time and Departure Time. For line N, it is necessary to ensure at least that the sum of the maximum departure time in the timetable and the travel time is greater than the maximum card swiping time for transfer passengers. Assuming that all card swiping times of the aith bus online N is set to CNai=c1,c2,…, and cmax is the maximum card swiping time in the set. The constraint of the relationship between the maximum card swiping time and the departure time is set out in (11) as follows:(11)xNai+2HN≥cmax,∀ai∈1,2,…,A,N∈X.(f) Route N Bus Travel Time Floating Upper and Lower Bound. For line N, assuming that the travel time of the bus is a floating J based on the original data, the upper and lower limits of the travel time of the line N bus can be obtained as shown in:(12)HN−J≤HN′≤HN+J,∀N∈X.(g) Constraint of Relationship between the Total Drop-Off Time and the Total Boarding Time. Obviously, given that the total drop-off time of all transfer passengers online M is ∑i=1PMNti+HM and the total optimized boarding time is ∑i=1PMNxNai+HN′, the constraint of the relationship between the total drop-off time and the total boarding time is defined as follows:(13)∑i=1PMNxNai+HN′≥∑i=1PMNti+HM,∀ai∈1,2,…,A,M,N∈X.


### 3.3. Establishment and Optimization of the Model

#### 3.3.1. Consideration and Selection of Random Variables

The 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. Therefore, we consider the travel time of line N as a random variable, defined asHN′˜. Then, the objective function and the corresponding constraint conditions are simplified by chance-constrained programming.


#### 3.3.2. Objective Function Optimization

First, we defineF1¯and F2¯ as the optimal values of the two objective functions obtained by solving the deterministic programming through the above modeling process without considering random variables. F1¯ and F2¯ are the minimum of F1 and the maximum of F2 respectively. By counting the travel times of buses on route N from route M to route N on different days, it is clear that HN′˜ follows a normal distribution, thus converting it from random variables to deterministic equivalents uN+σ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:(14)Prmin∑i=1PMNxNai+HN′˜−∑i=1PMNti+HM≤F1¯=λ1,minF2=min∑i=1PMNxNai−fNai,(15)maxλ1,minF2=min∑i=1PMNxNai−fNai.LetD=1−λ1, convert the above union expression to the formula:(16)minD=1−λ1,minF2=min∑i=1PMNxNai−fNai.


#### 3.3.3. 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 HN′˜, we change the model constraint equation (11) to the following equation (17), where αN is the confidence level of constraint satisfaction.(17)Prcmax−xNai≤2HN′˜≥αN,∀ai∈1,2,…,A,N∈X.LetN=cmax−xNai−2HN′˜, then it satisfies the following:(18)EzN=cmax−xNai−2uN,(19)DzN=4σN2=2σN2.There is a constraintcmax−xNai≤2HN′˜ equivalent to cmax−xNai−2HN′˜≤0 that satisfies equation (20) and can also be expressed as equation (21):(20)cmax−xNai−2HN′˜−cmax−xNai−2uN+cmax−xNai−2uN2σN≤0,(21)cmax−xNai−HN′˜−cmax−xNai−uN2σN≤−cmax−xNai−uN2σN.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):(22)η=HN′˜−2uN2σN,(23)Φη=12π∫−∞ηe−t2/2dt.Therefore, equation (21) can be transformed in form as shown in the following:(24)−cmax−xNai−2HN′˜−cmax−xNai−2uN2σN≥cmax−xNai−2uN2σN,(25)η≥cmax−xNai−2uN2σN.According to probability theory, there must be a numberKa for the confidence level αN, which makes the following formula (26) true. So the derivation can be done as shown in equations (27)–(30):(26)PrKa≤η=αN(27)Ka=Φ−11−αN(28)Prcmax−xNai−2uN2σN≤η≥αN(29)cmax−xNai−2uN2σN≤Ka(30)cmax−xNai−2uN2σN≤Φ−11−αN.Then, the constraintxNa+HN′˜≥cmax is equivalent to equation (31), where Φ−1αN is an inverse function of αN. Further merging and simplification lead to equation (32).(31)cmax−xNai≤2uN+2σN⋅Φ−11−αN,(32)cmax−xNai≤2uN+σN⋅Φ−11−αN.(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:(33)D∈0,1.


#### 3.3.4. 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:(34)minDminF2=min∑i=1PMNxNai−fNaisubject.to∑i=1PMNxNai+uN+σN⋅Φ−1λ1−∑i=1PMNti+HM≤F1¯xNai+2uN+σN⋅Φ−11−αN≥cmax∀ai∈1,2,…,A,N∈XGmin≤xNai+1−xNai≤Gmax,∀ai∈1,2,…,A,N∈XfNai−K≤xNai≤fNai+K,∀ai∈1,2,…,A,N∈XxNai≤xNai+1,∀ai∈1,2,…,A,N∈XDearly≤xNai≤Dlate,∀ai∈1,2,…,A,N∈X∑i=1PMNxNai+uN+σN⋅Φ−1λ1≥∑i=1PMNti+HM,∀ai∈1,2,…,A,M,N∈XD=1−λ1D∈0,1.


## 3.1. Assumptions and Notations

### 3.1.1. Assumptions

The 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 beX, 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. Therefore, 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.(a)
Passenger transfer demand: considering the same bus IC card, based on the analysis of large sample data, for lineM and line N, if the bus IC card has an adjacent swiping record of line M and line N and the bus IC card has to swipe record of line N after line M, the number of transfer passengers from line M to line N is assumed to be plus 1. We count the number of bus IC cards with these two properties to determine the number of interchange passengers from line M to line N and record them as PMN to represent passenger transfer demand.(b)
The total drop-off time of transfer passengers onlineM: assuming the boarding time of the transfer passenger who holds the ith IC card is ti. All IC cards, namely all transfer passengers, are anonymously processed to avoid losing generality. Considering the statistical analysis of large samples, it is assumed that the total travel time of line M is TM. And the passenger travel time online M can be approximately equal to the mean value of TM, which is approximately substituted by 1/2TM and recorded as HM. Since the passenger’s time to get off is equal to the sum of the passenger’s boarding time and the ride time, the total drop-off time of transfer passengers online M can be calculated by the following formula:(1)li=ti+HM,LM=∑i=1PMNli,whereli is the drop-off time of the transfer passenger who holds the ith IC card in line M and LM is the total drop-off time of all transfer passengers online M.(c)
The total boarding time of transfer passengers onlineN: assuming the number of trips corresponding to the record of the ith IC card online N is ai. 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 QminN=q1,q2,…,qa,…,qA, and set qA and q1 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:(2)t=qA−q1A−1,fNai=q1+a−1t,wherefNai represents the departure time of the aith 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 lineN is TN. And the travel time of line N is approximately equal to the mean of TN, which is approximately substituted by 1/2TN and recorded as HN. The 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:(3)si=fNai+HN,SN=∑i=1PMNsi,wheresi is the boarding time of the transfer passenger who holds the ith IC card in line N, and SN is the total boarding time of all transfer passengers in line N.(d)
The total waiting time of passengers transferring from lineM 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. The total waiting time of transfer passengers from line M to line N can be obtained by the difference calculation of the formula:(4)TMN=SN−LM=∑i=1PMNfNai+HN−∑i=1PMNti+HM,whereTMN represents the total waiting time of passengers transferring from line M to line N.


### 3.1.2. Notations

The following terms have been defined and used in our model (Table1).Table 1
Compliance description of notations used in this study.
SetXRoutes set,N∈XQminNThe earliest card-brush time set of all vehicles in lineN, qi∈QminNCNaiAll card-brush time of thea(i)th bus in line N, ci∈CNaiParametersATotal number of cars in routenPMNPassenger transfer demandLMTotal drop-off time of all transfer passengers in lineMTMTotal vehicle travel time for lineMTNTotal vehicle travel time for lineNHMMean of total vehicle travel time for lineMHNMean of total vehicle travel time for lineNtiBoarding time of the transfer passenger who holds theith IC cardliDrop-off time of the transfer passenger who holds theith IC card in line MsiBoarding time of the transfer passenger who holds theith IC card in line NaiNumber of trips corresponding to the record of theith IC card in line NqiThe earliest card swiping time of theith IC bus in line NSNTotal boarding time of all transfer passengers in lineNfNaiDeparture time of theaith bus online N in the existing timetable, ∀ai∈1,2,…,ATMNTotal waiting time of passengers transferring from lineM to line NBTotal departure adjustment time of lineM to line NHN′Optimized bus travel time onlineNGmaxMaximum departure time intervalGminMinimum departure time intervalKAdjustable range of departure timetableJFloating range of a vehicle’s travel timeDearlyThe earliest departure times of a bus lineDlateThe latest departure times of a bus linecmaxMaximum card swiping time of theaith bus online NF1¯Optimal value ofF1 derived from deterministic programmingF2¯Optimal value ofF2 derived from deterministic programmingHN′˜Random travel time of vehicles onlineNαNConfidence level of constraint satisfaction for stochastic programminguNMean of random travel time in lineNσNVariance of random travel time in lineNλ1Probability that the random planning transfer time after considering the random variable is less than the deterministic planning transfer timeDecision variablesD1 −λ1xNaiDeparture time of theai-th bus online N in the optimized timetable


## 3.1.1. Assumptions

The 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 beX, 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. Therefore, 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.(a)
Passenger transfer demand: considering the same bus IC card, based on the analysis of large sample data, for lineM and line N, if the bus IC card has an adjacent swiping record of line M and line N and the bus IC card has to swipe record of line N after line M, the number of transfer passengers from line M to line N is assumed to be plus 1. We count the number of bus IC cards with these two properties to determine the number of interchange passengers from line M to line N and record them as PMN to represent passenger transfer demand.(b)
The total drop-off time of transfer passengers onlineM: assuming the boarding time of the transfer passenger who holds the ith IC card is ti. All IC cards, namely all transfer passengers, are anonymously processed to avoid losing generality. Considering the statistical analysis of large samples, it is assumed that the total travel time of line M is TM. And the passenger travel time online M can be approximately equal to the mean value of TM, which is approximately substituted by 1/2TM and recorded as HM. Since the passenger’s time to get off is equal to the sum of the passenger’s boarding time and the ride time, the total drop-off time of transfer passengers online M can be calculated by the following formula:(1)li=ti+HM,LM=∑i=1PMNli,whereli is the drop-off time of the transfer passenger who holds the ith IC card in line M and LM is the total drop-off time of all transfer passengers online M.(c)
The total boarding time of transfer passengers onlineN: assuming the number of trips corresponding to the record of the ith IC card online N is ai. 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 QminN=q1,q2,…,qa,…,qA, and set qA and q1 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:(2)t=qA−q1A−1,fNai=q1+a−1t,wherefNai represents the departure time of the aith 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 lineN is TN. And the travel time of line N is approximately equal to the mean of TN, which is approximately substituted by 1/2TN and recorded as HN. The 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:(3)si=fNai+HN,SN=∑i=1PMNsi,wheresi is the boarding time of the transfer passenger who holds the ith IC card in line N, and SN is the total boarding time of all transfer passengers in line N.(d)
The total waiting time of passengers transferring from lineM 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. The total waiting time of transfer passengers from line M to line N can be obtained by the difference calculation of the formula:(4)TMN=SN−LM=∑i=1PMNfNai+HN−∑i=1PMNti+HM,whereTMN represents the total waiting time of passengers transferring from line M to line N.


## 3.1.2. Notations

The following terms have been defined and used in our model (Table1).Table 1
Compliance description of notations used in this study.
SetXRoutes set,N∈XQminNThe earliest card-brush time set of all vehicles in lineN, qi∈QminNCNaiAll card-brush time of thea(i)th bus in line N, ci∈CNaiParametersATotal number of cars in routenPMNPassenger transfer demandLMTotal drop-off time of all transfer passengers in lineMTMTotal vehicle travel time for lineMTNTotal vehicle travel time for lineNHMMean of total vehicle travel time for lineMHNMean of total vehicle travel time for lineNtiBoarding time of the transfer passenger who holds theith IC cardliDrop-off time of the transfer passenger who holds theith IC card in line MsiBoarding time of the transfer passenger who holds theith IC card in line NaiNumber of trips corresponding to the record of theith IC card in line NqiThe earliest card swiping time of theith IC bus in line NSNTotal boarding time of all transfer passengers in lineNfNaiDeparture time of theaith bus online N in the existing timetable, ∀ai∈1,2,…,ATMNTotal waiting time of passengers transferring from lineM to line NBTotal departure adjustment time of lineM to line NHN′Optimized bus travel time onlineNGmaxMaximum departure time intervalGminMinimum departure time intervalKAdjustable range of departure timetableJFloating range of a vehicle’s travel timeDearlyThe earliest departure times of a bus lineDlateThe latest departure times of a bus linecmaxMaximum card swiping time of theaith bus online NF1¯Optimal value ofF1 derived from deterministic programmingF2¯Optimal value ofF2 derived from deterministic programmingHN′˜Random travel time of vehicles onlineNαNConfidence level of constraint satisfaction for stochastic programminguNMean of random travel time in lineNσNVariance of random travel time in lineNλ1Probability that the random planning transfer time after considering the random variable is less than the deterministic planning transfer timeDecision variablesD1 −λ1xNaiDeparture time of theai-th bus online N in the optimized timetable


## 3.2. Model Analysis

### 3.2.1. Objective Function Analysis

The 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 lineM to line N. The total bus departure adjustment time is equal to the deviation value between the optimized timetable and the original timetable of each bus. Thus, 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:(5)B=min∑i=1PMNxNai−fNai,where B represents the total departure adjustment time of line M to line N.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. The objective of the optimization is to minimize both the objective functionF1 and F2, defined as follows:(6)minF1=∑i=1PMNxNai+HN′−∑i=1PMNti+HMminF2=min∑i=1PMNxNai−fNai,,where xNai is the departure time of the aith bus online N in the optimized timetable and HN′ represents optimized bus travel time.


### 3.2.2. 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 Gmax and the minimum is Gmin, the vehicle headway constraint can be obtained as follows:(7)Gmin≤xNai+1−xNai≤Gmax,∀ai∈1,2,…,A,N∈X.(b) Timetable Adjustment Range Constraint. From 3.1.1 (3), we can see that the departure time of the aith bus online N in the existing timetable is fNai. Assuming that based on the existing timetable, the departure time of each vehicle can be fluctuated up and down K ranges. Therefore, the adjustable departure time constraints are defined as follows:(8)fNai−K≤xNai≤fNai+K,∀ai∈1,2,…,A,N∈X.(c) Departure Time Sequence Constraint. 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:(9)xNai≤xNai+1,∀ai∈1,2,…,A,N∈X.(d) Reasonable Departure Time Constraint. Assume that Dearly and Dlate represent the earliest and latest departure times of a bus line, respectively, and that the departure times of each bus cannot be earlier than Dearly and not later than Dlate. The reasonable departure time constraint can be obtained as follows:(10)Dearly≤xNai≤Dlate,∀ai∈1,2,…,A,N∈X.(e) Constraint of Relationship between Maximum Swipe Time and Departure Time. For line N, it is necessary to ensure at least that the sum of the maximum departure time in the timetable and the travel time is greater than the maximum card swiping time for transfer passengers. Assuming that all card swiping times of the aith bus online N is set to CNai=c1,c2,…, and cmax is the maximum card swiping time in the set. The constraint of the relationship between the maximum card swiping time and the departure time is set out in (11) as follows:(11)xNai+2HN≥cmax,∀ai∈1,2,…,A,N∈X.(f) Route N Bus Travel Time Floating Upper and Lower Bound. For line N, assuming that the travel time of the bus is a floating J based on the original data, the upper and lower limits of the travel time of the line N bus can be obtained as shown in:(12)HN−J≤HN′≤HN+J,∀N∈X.(g) Constraint of Relationship between the Total Drop-Off Time and the Total Boarding Time. Obviously, given that the total drop-off time of all transfer passengers online M is ∑i=1PMNti+HM and the total optimized boarding time is ∑i=1PMNxNai+HN′, the constraint of the relationship between the total drop-off time and the total boarding time is defined as follows:(13)∑i=1PMNxNai+HN′≥∑i=1PMNti+HM,∀ai∈1,2,…,A,M,N∈X.


## 3.2.1. Objective Function Analysis

The 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 lineM to line N. The total bus departure adjustment time is equal to the deviation value between the optimized timetable and the original timetable of each bus. Thus, 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:(5)B=min∑i=1PMNxNai−fNai,where B represents the total departure adjustment time of line M to line N.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. The objective of the optimization is to minimize both the objective functionF1 and F2, defined as follows:(6)minF1=∑i=1PMNxNai+HN′−∑i=1PMNti+HMminF2=min∑i=1PMNxNai−fNai,,where xNai is the departure time of the aith bus online N in the optimized timetable and HN′ represents optimized bus travel time.


## 3.2.2. 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 Gmax and the minimum is Gmin, the vehicle headway constraint can be obtained as follows:(7)Gmin≤xNai+1−xNai≤Gmax,∀ai∈1,2,…,A,N∈X.(b) Timetable Adjustment Range Constraint. From 3.1.1 (3), we can see that the departure time of the aith bus online N in the existing timetable is fNai. Assuming that based on the existing timetable, the departure time of each vehicle can be fluctuated up and down K ranges. Therefore, the adjustable departure time constraints are defined as follows:(8)fNai−K≤xNai≤fNai+K,∀ai∈1,2,…,A,N∈X.(c) Departure Time Sequence Constraint. 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:(9)xNai≤xNai+1,∀ai∈1,2,…,A,N∈X.(d) Reasonable Departure Time Constraint. Assume that Dearly and Dlate represent the earliest and latest departure times of a bus line, respectively, and that the departure times of each bus cannot be earlier than Dearly and not later than Dlate. The reasonable departure time constraint can be obtained as follows:(10)Dearly≤xNai≤Dlate,∀ai∈1,2,…,A,N∈X.(e) Constraint of Relationship between Maximum Swipe Time and Departure Time. For line N, it is necessary to ensure at least that the sum of the maximum departure time in the timetable and the travel time is greater than the maximum card swiping time for transfer passengers. Assuming that all card swiping times of the aith bus online N is set to CNai=c1,c2,…, and cmax is the maximum card swiping time in the set. The constraint of the relationship between the maximum card swiping time and the departure time is set out in (11) as follows:(11)xNai+2HN≥cmax,∀ai∈1,2,…,A,N∈X.(f) Route N Bus Travel Time Floating Upper and Lower Bound. For line N, assuming that the travel time of the bus is a floating J based on the original data, the upper and lower limits of the travel time of the line N bus can be obtained as shown in:(12)HN−J≤HN′≤HN+J,∀N∈X.(g) Constraint of Relationship between the Total Drop-Off Time and the Total Boarding Time. Obviously, given that the total drop-off time of all transfer passengers online M is ∑i=1PMNti+HM and the total optimized boarding time is ∑i=1PMNxNai+HN′, the constraint of the relationship between the total drop-off time and the total boarding time is defined as follows:(13)∑i=1PMNxNai+HN′≥∑i=1PMNti+HM,∀ai∈1,2,…,A,M,N∈X.


## 3.3. Establishment and Optimization of the Model

### 3.3.1. Consideration and Selection of Random Variables

The 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. Therefore, we consider the travel time of line N as a random variable, defined asHN′˜. Then, the objective function and the corresponding constraint conditions are simplified by chance-constrained programming.


### 3.3.2. Objective Function Optimization

First, we defineF1¯and F2¯ as the optimal values of the two objective functions obtained by solving the deterministic programming through the above modeling process without considering random variables. F1¯ and F2¯ are the minimum of F1 and the maximum of F2 respectively. By counting the travel times of buses on route N from route M to route N on different days, it is clear that HN′˜ follows a normal distribution, thus converting it from random variables to deterministic equivalents uN+σ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:(14)Prmin∑i=1PMNxNai+HN′˜−∑i=1PMNti+HM≤F1¯=λ1,minF2=min∑i=1PMNxNai−fNai,(15)maxλ1,minF2=min∑i=1PMNxNai−fNai.LetD=1−λ1, convert the above union expression to the formula:(16)minD=1−λ1,minF2=min∑i=1PMNxNai−fNai.


### 3.3.3. 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 HN′˜, we change the model constraint equation (11) to the following equation (17), where αN is the confidence level of constraint satisfaction.(17)Prcmax−xNai≤2HN′˜≥αN,∀ai∈1,2,…,A,N∈X.LetN=cmax−xNai−2HN′˜, then it satisfies the following:(18)EzN=cmax−xNai−2uN,(19)DzN=4σN2=2σN2.There is a constraintcmax−xNai≤2HN′˜ equivalent to cmax−xNai−2HN′˜≤0 that satisfies equation (20) and can also be expressed as equation (21):(20)cmax−xNai−2HN′˜−cmax−xNai−2uN+cmax−xNai−2uN2σN≤0,(21)cmax−xNai−HN′˜−cmax−xNai−uN2σN≤−cmax−xNai−uN2σN.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):(22)η=HN′˜−2uN2σN,(23)Φη=12π∫−∞ηe−t2/2dt.Therefore, equation (21) can be transformed in form as shown in the following:(24)−cmax−xNai−2HN′˜−cmax−xNai−2uN2σN≥cmax−xNai−2uN2σN,(25)η≥cmax−xNai−2uN2σN.According to probability theory, there must be a numberKa for the confidence level αN, which makes the following formula (26) true. So the derivation can be done as shown in equations (27)–(30):(26)PrKa≤η=αN(27)Ka=Φ−11−αN(28)Prcmax−xNai−2uN2σN≤η≥αN(29)cmax−xNai−2uN2σN≤Ka(30)cmax−xNai−2uN2σN≤Φ−11−αN.Then, the constraintxNa+HN′˜≥cmax is equivalent to equation (31), where Φ−1αN is an inverse function of αN. Further merging and simplification lead to equation (32).(31)cmax−xNai≤2uN+2σN⋅Φ−11−αN,(32)cmax−xNai≤2uN+σN⋅Φ−11−αN.(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:(33)D∈0,1.


### 3.3.4. 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:(34)minDminF2=min∑i=1PMNxNai−fNaisubject.to∑i=1PMNxNai+uN+σN⋅Φ−1λ1−∑i=1PMNti+HM≤F1¯xNai+2uN+σN⋅Φ−11−αN≥cmax∀ai∈1,2,…,A,N∈XGmin≤xNai+1−xNai≤Gmax,∀ai∈1,2,…,A,N∈XfNai−K≤xNai≤fNai+K,∀ai∈1,2,…,A,N∈XxNai≤xNai+1,∀ai∈1,2,…,A,N∈XDearly≤xNai≤Dlate,∀ai∈1,2,…,A,N∈X∑i=1PMNxNai+uN+σN⋅Φ−1λ1≥∑i=1PMNti+HM,∀ai∈1,2,…,A,M,N∈XD=1−λ1D∈0,1.


## 3.3.1. Consideration and Selection of Random Variables

The 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. Therefore, we consider the travel time of line N as a random variable, defined asHN′˜. Then, the objective function and the corresponding constraint conditions are simplified by chance-constrained programming.


## 3.3.2. Objective Function Optimization

First, we defineF1¯and F2¯ as the optimal values of the two objective functions obtained by solving the deterministic programming through the above modeling process without considering random variables. F1¯ and F2¯ are the minimum of F1 and the maximum of F2 respectively. By counting the travel times of buses on route N from route M to route N on different days, it is clear that HN′˜ follows a normal distribution, thus converting it from random variables to deterministic equivalents uN+σ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:(14)Prmin∑i=1PMNxNai+HN′˜−∑i=1PMNti+HM≤F1¯=λ1,minF2=min∑i=1PMNxNai−fNai,(15)maxλ1,minF2=min∑i=1PMNxNai−fNai.LetD=1−λ1, convert the above union expression to the formula:(16)minD=1−λ1,minF2=min∑i=1PMNxNai−fNai.


## 3.3.3. 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 HN′˜, we change the model constraint equation (11) to the following equation (17), where αN is the confidence level of constraint satisfaction.(17)Prcmax−xNai≤2HN′˜≥αN,∀ai∈1,2,…,A,N∈X.LetN=cmax−xNai−2HN′˜, then it satisfies the following:(18)EzN=cmax−xNai−2uN,(19)DzN=4σN2=2σN2.There is a constraintcmax−xNai≤2HN′˜ equivalent to cmax−xNai−2HN′˜≤0 that satisfies equation (20) and can also be expressed as equation (21):(20)cmax−xNai−2HN′˜−cmax−xNai−2uN+cmax−xNai−2uN2σN≤0,(21)cmax−xNai−HN′˜−cmax−xNai−uN2σN≤−cmax−xNai−uN2σN.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):(22)η=HN′˜−2uN2σN,(23)Φη=12π∫−∞ηe−t2/2dt.Therefore, equation (21) can be transformed in form as shown in the following:(24)−cmax−xNai−2HN′˜−cmax−xNai−2uN2σN≥cmax−xNai−2uN2σN,(25)η≥cmax−xNai−2uN2σN.According to probability theory, there must be a numberKa for the confidence level αN, which makes the following formula (26) true. So the derivation can be done as shown in equations (27)–(30):(26)PrKa≤η=αN(27)Ka=Φ−11−αN(28)Prcmax−xNai−2uN2σN≤η≥αN(29)cmax−xNai−2uN2σN≤Ka(30)cmax−xNai−2uN2σN≤Φ−11−αN.Then, the constraintxNa+HN′˜≥cmax is equivalent to equation (31), where Φ−1αN is an inverse function of αN. Further merging and simplification lead to equation (32).(31)cmax−xNai≤2uN+2σN⋅Φ−11−αN,(32)cmax−xNai≤2uN+σN⋅Φ−11−αN.(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:(33)D∈0,1.


## 3.3.4. 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:(34)minDminF2=min∑i=1PMNxNai−fNaisubject.to∑i=1PMNxNai+uN+σN⋅Φ−1λ1−∑i=1PMNti+HM≤F1¯xNai+2uN+σN⋅Φ−11−αN≥cmax∀ai∈1,2,…,A,N∈XGmin≤xNai+1−xNai≤Gmax,∀ai∈1,2,…,A,N∈XfNai−K≤xNai≤fNai+K,∀ai∈1,2,…,A,N∈XxNai≤xNai+1,∀ai∈1,2,…,A,N∈XDearly≤xNai≤Dlate,∀ai∈1,2,…,A,N∈X∑i=1PMNxNai+uN+σN⋅Φ−1λ1≥∑i=1PMNti+HM,∀ai∈1,2,…,A,M,N∈XD=1−λ1D∈0,1.


## 4. Algorithm Design

The 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.The Pareto solution is defined as if there are any two solutionsS1 and S2, S1 is superior to S2 for all objectives, we call S1 dominates S2. And if the solution of S1 is not dominated by other solutions, S1 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 S1 and S2 have fS1<fS2 for all targets, but at least one objective function has gS1<gS2, S1 is called weak dominating S2, that is, S1 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. Therefore, the augmented epsilon-constraint method can be well applied to our proposed two-objective 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 tominzx=z1x,z2x,x∈X.Here,z1x represents the objective function D (the probability that the stochastic programming transfer time is less than the deterministic programming transfer time), and z2x represents the objective function F2 (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. Therefore, the priority ofz1x is greater than that of z2x. Thus, by minimizing z1x, constraining z2x, and setting σ2 as the constraint value of z2x, the standard epsilon-constraint model of the original problem is(35)minz¯x=z1x,subject. to.z2x≤σ2,x∈X.The 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]. The dominant solution has not reached the optimization effect of the Pareto solution, so we should avoid the weak dominant solution as far as possible. Therefore, we need to improve the original epsilon-constraint algorithm.


### 4.2. Customized Augmented Epsilon-Constraint Algorithm

The 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. The 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. The augmented epsilon-constraint algorithm ensures that only valid Pareto solutions are searched. Here are the steps of the augmented epsilon-constraint algorithm to solve this scenario in this study.Step 1.Get the range of values of z1x and z2x as z1min,z1max and z2min,z2max, where zimin and zimax are the minimum and maximum values of the objective function zix, respectively. When z1x obtains the minimum value z1min, let x1 be the optimal solution at this time, z2x will get the maximum value z2max=z2x1 in the feasible region. Similarly, when z2x obtains the minimum value z2min, let x2 be the optimal solution at this time, z1x will get the maximum value z2max=z2x1 in the feasible region. The two objective functions are optimized for lexicographic order as follows:(36)z1max=minz1x,s.t.z2x=z2min,x∈X,z2max=minz2x,s.z1x=z1min,x∈X.Step 2. Mesh the range of values of z2x. If we want to get g2+1 Parato solutions, there will be g2+1 grid points (including the grid points corresponding to the minimum and maximum of z2x). Then, the range of values will be meshed into g2 grids, with each grid having an interval size of z2max−z2min/g2.Step 3. Determine the constraint level for each z2x in each grid interval. To facilitate the presentation of the algorithm, we set k2 as the index of the grid points, numbering from 0, so that the grid points k2=0,1,…, g2 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. The ratio of u2 to z1 is calculated to eliminate the dimensional effect. Then multiplied by the sufficiently small weight p, and extended to the objective function 1, the augmented epsilon-constraint method model is constructed. The final model is as follows:(37)minz¯1x−ρμ2v2,s.t.z¯2x+μ2=σ2;μ2≥0;x∈X;k2=0,1,…,g2.For each latticek2, a corresponding Pareto solution z¯1x∗,z¯2x∗ and the current optimal solution x∗,μ2∗ will be obtained for the above model. If μ2∗≥z2max−z2min/g2, 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. The 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=μ2g2/6500. If ξ2≥1, then let k2=k2+ξ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.The algorithm flow chart is shown in Figure1.Figure 1
Algorithm flow chart of input data, model, algorithm, and output.


## 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 tominzx=z1x,z2x,x∈X.Here,z1x represents the objective function D (the probability that the stochastic programming transfer time is less than the deterministic programming transfer time), and z2x represents the objective function F2 (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. Therefore, the priority ofz1x is greater than that of z2x. Thus, by minimizing z1x, constraining z2x, and setting σ2 as the constraint value of z2x, the standard epsilon-constraint model of the original problem is(35)minz¯x=z1x,subject. to.z2x≤σ2,x∈X.The 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]. The dominant solution has not reached the optimization effect of the Pareto solution, so we should avoid the weak dominant solution as far as possible. Therefore, we need to improve the original epsilon-constraint algorithm.


## 4.2. Customized Augmented Epsilon-Constraint Algorithm

The 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. The 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. The augmented epsilon-constraint algorithm ensures that only valid Pareto solutions are searched. Here are the steps of the augmented epsilon-constraint algorithm to solve this scenario in this study.Step 1.Get the range of values of z1x and z2x as z1min,z1max and z2min,z2max, where zimin and zimax are the minimum and maximum values of the objective function zix, respectively. When z1x obtains the minimum value z1min, let x1 be the optimal solution at this time, z2x will get the maximum value z2max=z2x1 in the feasible region. Similarly, when z2x obtains the minimum value z2min, let x2 be the optimal solution at this time, z1x will get the maximum value z2max=z2x1 in the feasible region. The two objective functions are optimized for lexicographic order as follows:(36)z1max=minz1x,s.t.z2x=z2min,x∈X,z2max=minz2x,s.z1x=z1min,x∈X.Step 2. Mesh the range of values of z2x. If we want to get g2+1 Parato solutions, there will be g2+1 grid points (including the grid points corresponding to the minimum and maximum of z2x). Then, the range of values will be meshed into g2 grids, with each grid having an interval size of z2max−z2min/g2.Step 3. Determine the constraint level for each z2x in each grid interval. To facilitate the presentation of the algorithm, we set k2 as the index of the grid points, numbering from 0, so that the grid points k2=0,1,…, g2 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. The ratio of u2 to z1 is calculated to eliminate the dimensional effect. Then multiplied by the sufficiently small weight p, and extended to the objective function 1, the augmented epsilon-constraint method model is constructed. The final model is as follows:(37)minz¯1x−ρμ2v2,s.t.z¯2x+μ2=σ2;μ2≥0;x∈X;k2=0,1,…,g2.For each latticek2, a corresponding Pareto solution z¯1x∗,z¯2x∗ and the current optimal solution x∗,μ2∗ will be obtained for the above model. If μ2∗≥z2max−z2min/g2, 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. The 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=μ2g2/6500. If ξ2≥1, then let k2=k2+ξ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.The algorithm flow chart is shown in Figure1.Figure 1
Algorithm flow chart of input data, model, algorithm, and output.


## 5. Example Analysis

### 5.1. Overview and Visualization Analysis of Examples

This 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. This study chooses a downward direction bus to study to simplify the problem.Figure2 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. The 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 different 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). Therefore, for the convenience of the study, the selected bus network is simplified. The simplified diagram is shown in Figure 3. The arrows on each line in the diagram indicate the direction in which the bus will run.Figure 2
A real-world bus network in Shenzhen, China.Figure 3
A simplified diagrammatic illustration of bus network.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 Table2. 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. The 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 first eight digits represent the specific transaction date, and the last six digits represent the specific transaction time (in 24-hour format). Each bus has a unique license plate number, through the license plate number can be obtained in the number of bus routes, and each bus in the study period of the number of vehicles.Table 2
Passenger data samples of bus lines.
IC card numberDate and time of transactionRouteLicense plate number66548344820161223071249Bus route 338YUE B9377232921994720161223084405Bus route 338YUE B9375302081650920161223063450Bus route 338YUE BBB43968433596120161223081852Bus route 338YUE BL718066673461820161223072423Bus route 338YUE B8739368406046020161223071933Bus route 338YUE BL828766548344820161223082307Bus route 615YUE BH097232921994720161223093908Bus route 615YUE BBA94702081650920161223071808Bus route 615YUE B8690068433596120161223090336Bus route 615YUE BBB70266673461820161223073903Bus route 615YUE BBB74368406046020161223074057Bus route 615YUE BBB479


### 5.2. Summary and Presentation of Calculation Results

This model is nonconvex optimization. The 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 confidence level. We transform the stochastic programming problem into a deterministic programming problem using chance-constrained programming. This study assumes that the travel time of each bus obeys the normal distributionNuN,σN2, where uN,σN are the mean and variance of the travel time of each bus. We transform the objective F1 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 F2 of the initial model from the sum of the absolute values of the optimized timetable and the original timetable offset to the sum of the squares of the timetable offset, which can be effectively 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. The Pareto approximate optimal solution of the multi-objective optimization problem is obtained by test example. The target values of Pareto nondominated front and Pareto solutions are shown in Figure4 and Table 3, respectively.Figure 4
Pareto optimal front for the network example shown in Figure2.Table 3
Target values of Pareto optimal solutions.
Pareto optimal solution no.DF210.4005043.69420.4084871.26430.4174698.83340.4264526.40250.4354353.97260.4454181.54170.4564009.11180.4663836.68090.4773664.250100.4883491.819It can be seen from Figure4 that the complementarity probability of the transfer time of stochastic programming less than that of deterministic programming is negatively correlated with the offset between the optimized timetable and the original timetable. It can be concluded that minimizing passenger transfer waiting time and minimizing timetable adjustment offset are two conflicting objectives when vehicle travel time is considered as a random variable. This is in line with the basic requirement that the Pareto optimal solution can be obtained only when the interests of each objective are conflicting in the multi-objective optimization problem. Therefore, the optimization objectives selected in this study are reasonable.From Table3, 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 offset from the original timetable when considering random variables. These 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 different shades of color represent different timetable offsets. The larger the color depth offset, the smaller the color light offset.Table 4
Adjustment time and optimized timetable of Pareto solution no.1.Table 5
Adjustment time and optimized timetable of Pareto solution no.10.By taking bus travel time as a random variable, the Pareto solution set ofD 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. The 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.Therefore, 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 Table4 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 goal will not have an obvious influence on the optimization of the departure time of the first vehicle. Similarly, for the bus departing at 8:37:04, all Pareto solutions have an optimized adjustment time range of 4 to 5 minutes (that is, the adjustment time range of the timetable is significantly smaller than that of other calculation results except for the first vehicle), and the corresponding optimization results are all in the time range of 8:32 to 8:33. Different optimization objectives will adjust this bus to the same time range. It shows that the bus operation department can get a smaller passenger transfer waiting time through a smaller timetable adjustment during this period.Figure5 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. This shows that the convergence of the augmented epsilon-constraint algorithm is good, and the algorithm can find the approximate optimal solution more stably. This solution can be regarded as a practical solution process.Figure 5
Algorithm iteration diagram.


### 5.3. Analysis and Discussion of Calculation Results

Then, sensitivity analysis method is used to evaluate the effect of parameterK and confidence α 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. The parameter α is taken 0.25 as the base, with each step of 0.01. The maximum increase is 0.05 upward, and the maximum decrease is 0.04 downward. The sensitivity analysis results of minimizing only the objective function D and minimizing only the sum of the departure time adjustments F2 are shown in Tables 6 and 7 respectively.Table 6
Sensitivity analysis of parameterK.
K10152025303540minD0.4000.3830.3830.3830.3830.3830.383F25043.6945447.3695447.3695447.3695447.3695447.3695447.369minF23491.8193532.1873532.1873532.1873532.1873532.1873532.187D0.4880.4860.4860.4860.4860.4860.486Table 7
Sensitivity analysis of parameterα.
α0.210.220.230.240.250.260.270.280.290.30minD0.3260.3320.3440.3560.3690.3830.4000.4190.4380.459F26500.0006357.5616107.6385866.4805601.6045353.9015043.6944710.9964398.8324055.678minF23637.4503623.2063598.2133574.0983547.6103522.8403491.8193458.5493428.4753407.052D0.4790.4800.4820.4830.4850.4860.4880.4910.4930.495As can be seen from Figure6(a), the parameter K has a significant effect on the target function D. The 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 F2 increases monotonously and tends to be stable with the increase of the value of the parameter K. The reason is that the adjustable range of the constraint becomes larger, the constraint becomes more relaxed, the allowable range of fluctuations becomes larger, and the total adjustment time F2 is correspondingly larger.Figure 6
Sensitivity analysis of parameterK. (a) D takes the minimum. (b) F2 takes the minimum.
(a)(b)As shown in Figure6(b), parameter K has a significant effect on the total value of departure time adjustment F2. The total value of departure time adjustment F2 increases first and then stabilizes with the increase of parameter K, fluctuates 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 first and then stabilizes with the increase of parameter K, fluctuates slightly in a very small range, and maximizes when parameter K is 10 min.To sum up, we can consider taking a better value ofK 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 Figure7(a), the confidence α has a significant effect on the target function D, the value of the target D increases monotonously with the increase of the confidence α, and the sum of the minimum departure time adjustment F2 decreases monotonously with the rise of the confidence α. From Figure 7(b), it can be seen that confidence α has a significant effect on the sum of departure time adjustment values F2, and the sum of departure time adjustment values F2 decreases monotonously with the increase of confidence α values. In contrast, the value of the target D increases monotonously with the increase of confidence α values. To sum up, the confidence α of the compromise size can be taken into account to achieve the balance of minimizing the values of the two objective functions.Figure 7
Sensitivity analysis of confidenceα. (a) D takes the minimum. (b) F2 takes the minimum.
(a)(b)The 10 nondominated solutions solved show different application scenarios and meanings. Considering the meaning of the value of the objective function, the value of the objective functionF1, 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. The smaller the value of D, the greater the probability that passengers needless waiting time. For objective function F2, 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. Thus, 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.


## 5.1. Overview and Visualization Analysis of Examples

This 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. This study chooses a downward direction bus to study to simplify the problem.Figure2 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. The 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 different 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). Therefore, for the convenience of the study, the selected bus network is simplified. The simplified diagram is shown in Figure 3. The arrows on each line in the diagram indicate the direction in which the bus will run.Figure 2
A real-world bus network in Shenzhen, China.Figure 3
A simplified diagrammatic illustration of bus network.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 Table2. 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. The 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 first eight digits represent the specific transaction date, and the last six digits represent the specific transaction time (in 24-hour format). Each bus has a unique license plate number, through the license plate number can be obtained in the number of bus routes, and each bus in the study period of the number of vehicles.Table 2
Passenger data samples of bus lines.
IC card numberDate and time of transactionRouteLicense plate number66548344820161223071249Bus route 338YUE B9377232921994720161223084405Bus route 338YUE B9375302081650920161223063450Bus route 338YUE BBB43968433596120161223081852Bus route 338YUE BL718066673461820161223072423Bus route 338YUE B8739368406046020161223071933Bus route 338YUE BL828766548344820161223082307Bus route 615YUE BH097232921994720161223093908Bus route 615YUE BBA94702081650920161223071808Bus route 615YUE B8690068433596120161223090336Bus route 615YUE BBB70266673461820161223073903Bus route 615YUE BBB74368406046020161223074057Bus route 615YUE BBB479


## 5.2. Summary and Presentation of Calculation Results

This model is nonconvex optimization. The 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 confidence level. We transform the stochastic programming problem into a deterministic programming problem using chance-constrained programming. This study assumes that the travel time of each bus obeys the normal distributionNuN,σN2, where uN,σN are the mean and variance of the travel time of each bus. We transform the objective F1 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 F2 of the initial model from the sum of the absolute values of the optimized timetable and the original timetable offset to the sum of the squares of the timetable offset, which can be effectively 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. The Pareto approximate optimal solution of the multi-objective optimization problem is obtained by test example. The target values of Pareto nondominated front and Pareto solutions are shown in Figure4 and Table 3, respectively.Figure 4
Pareto optimal front for the network example shown in Figure2.Table 3
Target values of Pareto optimal solutions.
Pareto optimal solution no.DF210.4005043.69420.4084871.26430.4174698.83340.4264526.40250.4354353.97260.4454181.54170.4564009.11180.4663836.68090.4773664.250100.4883491.819It can be seen from Figure4 that the complementarity probability of the transfer time of stochastic programming less than that of deterministic programming is negatively correlated with the offset between the optimized timetable and the original timetable. It can be concluded that minimizing passenger transfer waiting time and minimizing timetable adjustment offset are two conflicting objectives when vehicle travel time is considered as a random variable. This is in line with the basic requirement that the Pareto optimal solution can be obtained only when the interests of each objective are conflicting in the multi-objective optimization problem. Therefore, the optimization objectives selected in this study are reasonable.From Table3, 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 offset from the original timetable when considering random variables. These 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 different shades of color represent different timetable offsets. The larger the color depth offset, the smaller the color light offset.Table 4
Adjustment time and optimized timetable of Pareto solution no.1.Table 5
Adjustment time and optimized timetable of Pareto solution no.10.By taking bus travel time as a random variable, the Pareto solution set ofD 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. The 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.Therefore, 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 Table4 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 goal will not have an obvious influence on the optimization of the departure time of the first vehicle. Similarly, for the bus departing at 8:37:04, all Pareto solutions have an optimized adjustment time range of 4 to 5 minutes (that is, the adjustment time range of the timetable is significantly smaller than that of other calculation results except for the first vehicle), and the corresponding optimization results are all in the time range of 8:32 to 8:33. Different optimization objectives will adjust this bus to the same time range. It shows that the bus operation department can get a smaller passenger transfer waiting time through a smaller timetable adjustment during this period.Figure5 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. This shows that the convergence of the augmented epsilon-constraint algorithm is good, and the algorithm can find the approximate optimal solution more stably. This solution can be regarded as a practical solution process.Figure 5
Algorithm iteration diagram.


## 5.3. Analysis and Discussion of Calculation Results

Then, sensitivity analysis method is used to evaluate the effect of parameterK and confidence α 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. The parameter α is taken 0.25 as the base, with each step of 0.01. The maximum increase is 0.05 upward, and the maximum decrease is 0.04 downward. The sensitivity analysis results of minimizing only the objective function D and minimizing only the sum of the departure time adjustments F2 are shown in Tables 6 and 7 respectively.Table 6
Sensitivity analysis of parameterK.
K10152025303540minD0.4000.3830.3830.3830.3830.3830.383F25043.6945447.3695447.3695447.3695447.3695447.3695447.369minF23491.8193532.1873532.1873532.1873532.1873532.1873532.187D0.4880.4860.4860.4860.4860.4860.486Table 7
Sensitivity analysis of parameterα.
α0.210.220.230.240.250.260.270.280.290.30minD0.3260.3320.3440.3560.3690.3830.4000.4190.4380.459F26500.0006357.5616107.6385866.4805601.6045353.9015043.6944710.9964398.8324055.678minF23637.4503623.2063598.2133574.0983547.6103522.8403491.8193458.5493428.4753407.052D0.4790.4800.4820.4830.4850.4860.4880.4910.4930.495As can be seen from Figure6(a), the parameter K has a significant effect on the target function D. The 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 F2 increases monotonously and tends to be stable with the increase of the value of the parameter K. The reason is that the adjustable range of the constraint becomes larger, the constraint becomes more relaxed, the allowable range of fluctuations becomes larger, and the total adjustment time F2 is correspondingly larger.Figure 6
Sensitivity analysis of parameterK. (a) D takes the minimum. (b) F2 takes the minimum.
(a)(b)As shown in Figure6(b), parameter K has a significant effect on the total value of departure time adjustment F2. The total value of departure time adjustment F2 increases first and then stabilizes with the increase of parameter K, fluctuates 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 first and then stabilizes with the increase of parameter K, fluctuates slightly in a very small range, and maximizes when parameter K is 10 min.To sum up, we can consider taking a better value ofK 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 Figure7(a), the confidence α has a significant effect on the target function D, the value of the target D increases monotonously with the increase of the confidence α, and the sum of the minimum departure time adjustment F2 decreases monotonously with the rise of the confidence α. From Figure 7(b), it can be seen that confidence α has a significant effect on the sum of departure time adjustment values F2, and the sum of departure time adjustment values F2 decreases monotonously with the increase of confidence α values. In contrast, the value of the target D increases monotonously with the increase of confidence α values. To sum up, the confidence α of the compromise size can be taken into account to achieve the balance of minimizing the values of the two objective functions.Figure 7
Sensitivity analysis of confidenceα. (a) D takes the minimum. (b) F2 takes the minimum.
(a)(b)The 10 nondominated solutions solved show different application scenarios and meanings. Considering the meaning of the value of the objective function, the value of the objective functionF1, 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. The smaller the value of D, the greater the probability that passengers needless waiting time. For objective function F2, 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. Thus, 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.


## 6. Conclusion

The 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 multi-objective 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 multi-objective programming model and simplify the model using the normal distribution property of travel time random variables. Thus, the chance-constrained programming method transforms random variables into deterministic variables. Finally, stochastic programming is transformed into a new deterministic multi-objective programming problem.This 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. The 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.



---

*Source: 9222636-2022-08-28.xml*