# A MultiObjective Optimization Approach for Integrated Timetabling and Vehicle Scheduling with Uncertainty

Authors: Yindong Shen; Wenliang Xie; Jingpeng Li
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

The timetabling problem (TTP) and vehicle scheduling problem (VSP) are two indispensable problems in public transit planning process. They used to be solved in sequence; hence, optimality of resulting solutions is compromised. To get better results, some integrated approaches emerge to solve the TTP and VSP as an integrated problem. In the existing integrated approaches, the passenger comfort on bus and the uncertainty in the real world are rarely considered. To provide better service for passengers and enhance the robustness of the schedule to be compiled, we study the integrated optimization of TTP and VSP with uncertainty. In this paper, a novel multiobjective optimization approach with the objectives of minimizing the passenger travel cost, the vehicle scheduling cost, and the incompatible trip-link cost is proposed. Meanwhile, a multiobjective hybrid algorithm, which is a combination of the self-adjust genetic algorithm (SGA), large neighborhood search (LNS) algorithm, and Pareto separation operator (PSO), is applied to solve the integrated optimization problem. The experimental results show that the approach outperforms existing approaches in terms of service level and robustness.


---

## Body

## 1. Introduction

The public transit planning process is usually divided in sequential steps: network route design, timetabling [1], vehicle scheduling [2], and crew scheduling and rostering [3]. The output of each activity positioned higher in the sequence becomes an important input for lower-level decisions [4]. Timetabling and vehicle scheduling are two important activities in the transit operational planning. The timetabling problem (TTP) is to find the best dispatching policy for transit vehicles with the objectives of minimizing passenger waiting time at stops and minimizing the headway unevenness governed by the headway smoothing rule [5]. These objectives refer to the service level provided to passengers. The vehicle scheduling problem (VSP) is concerned with determining the optimal allocation of vehicles to carry out all the trips in a given transit timetable. The objective is to minimize the operation cost related with the usage of vehicles and fuel consumption [6]. The TTP and VSP are usually solved individually in sequence. This may lead to a compromised solution for the optimization problem of timetabling and vehicle scheduling as a whole [7]. In order to meet the passenger demand with high social benefits and lower operators’ costs, it is necessary to coordinate the timetabling and vehicle scheduling. However, building a proper model and developing an integrated approach are challenging.The mostly used integrated approach to the TTP and VSP is the recursive approach. The basic idea can be described as follows: (1) generate an initial timetable based on passenger demand; (2) determine the minimum number of vehicles that is needed to execute the timetable; (3) change the trip chain to reduce the vehicle scheduling cost; and (4) adjust the initial timetable and return to step (2) until the generated vehicle schedule is acceptable [8]. For example, Ceder [9] proposed a solution approach to combine timetabling and vehicle scheduling so as to improve the correspondence of vehicle departure times with passenger demand while minimizing the fleet size required. The author stated that the proposed approach could obtain efficient schedules to the test instances from both the passenger and operator perspectives. Ceder [10] proposed a method to construct timetables based on even-headways and even-load with different vehicle types. The timetable was used as an input to solve the vehicle scheduling problem using the deficit function theory. Schmid and Ehmke [11] developed hybrid mathematical models and solution approaches to solve the vehicle scheduling problem with time windows as well as balanced departure times based on the research of Kliewer et al. [12], which discussed the multiple-depot, multiple-vehicle-type bus VSP while using time-space-based instead of connection-based networks for modeling.The recursive approach is hard to find the global optimum solution. In recent years, the approaches to build an integrated model of TTP and VSP have been widely studied. For instance, Petersen et al. [13] formulated a partial integrated model of TTP and VSP, in which the objectives were to minimize the costs of vehicles usage and passenger transfers. Meanwhile, a large neighborhood search metaheuristic approach was proposed. Carosi et al. [14] studied the integrated timetabling and vehicle scheduling and proposed a mixed-integer linear programming (MILP) multicommodity flow-type model to obtain the vehicle schedule with desired headways and minimum deviation from the ideal frequency of service. Yue et al. [15] proposed an optimization methodology that integrated both train timetabling and rolling stock scheduling based on time-dependent passenger flow demands and employed a simulated-annealing (SA)-based heuristic algorithm to solve the proposed model.Obviously, the integrated problem has multiple objectives, which conflict with each other. More and more researchers tend to establish multiobjective models and build intelligent optimization algorithms for the integrated problem. For example, Ibarra-Rojas et al. [16] proposed a biobjective optimization problem to jointly solve the single-depot VSP and the synchronization bus timetabling problem in which time windows for departure times were considered and the passenger demand were assumed to be constant. The objectives were to maximize the number of passengers benefited by well-timed transfers and to minimize the fleet size. Teng and Chen [17] focused on the TTP and VSP for electric buses and developed a multiobjective optimization model under single-line mode. The objectives included smoothing the vehicle departure intervals and minimizing the number of vehicles and total charging costs. A multiobjective particle swarm optimization (MOPSO) algorithm was developed to get the Pareto optimal solution set. Wang et al. [18] studied the integration of train scheduling and rolling stock circulation planning under time-varying passenger demand for an urban rail transit line. A multiobjective approximated mixed-integer linear programming approach was presented to provide better passenger services and reduce the operational complexity and costs of the scheduling.Table1 summarizes the relevant studies on the integration of TTP and VSP, in terms of the approach type of integration, objectives, solution method, and data type. Among the existing studies on the integrated problem, the passenger crowding on the bus is rarely considered. The timetable obtained without considering crowding on the bus might bring passengers unsatisfied or uncomfortable and increase the cost of the transport system. Moreover, these researches do not consider the uncertainty of the passenger flow and trip times, which exists commonly in the real-world bus operation. Evaluating the passenger travel cost that consists of waiting cost and comfort cost would be inaccurate regardless of the stochastic passenger flow. Furthermore, when the stochastic trip times are employed, some infeasible trip links in the traditional VSP with fixed and tight trip times might become compatible, and the resulting schedule might be more robust (i.e., with high on-time performance) [19]. Considering the stochastic passenger flow and trip time, the integrated problem would be more complex and hard to optimize. Therefore, to have better performance, an improved optimization strategy for the integrated problem with uncertainty should be proposed.Table 1
Summary of relevant studies on the integration of TTP and VSP.
Approach typePublicationObjectivesSolution methodData typeRecursive approachCeder [9]Maximize the correspondence of vehicle departure times, the fleet sizeDeficit function approachFixed passenger load at each stop, fixed trip timeCeder [4]Even-headways and even-load with different vehicle types, minimize vehicle scheduling costHeuristicsFixed passenger load at each stop, fixed trip timeSchmid and Ehmke [11]Minimize costs of operation, maximize quality of timetablesHybrid metaheuristicFixed time windows, fixed trip timeIntegrated modelPetersen et al. [13]Minimize the costs of vehicles usage and passenger transfersLNS approachFixed passenger volume estimates, fixed trip timeCarosi et al. [14]Minimum deviation from the ideal frequency of service, vehicle schedule costDiving-type approachChangeable time window, fixed trip timeYue et al. [15]Minimizing infeasible trains, waiting times for passengersSA algorithmFixed passenger volume, fixed trip timeMulti-objective modelIbarra-Rojas et al. [16]Maximize the number of passengers benefited, minimize the fleet sizeε-constraint methodFixed passenger volume estimates, fixed trip timeTeng and Chen [17]Smooth the vehicle departure intervals, minimize the number of vehicles and total charging costsMOPSO algorithmFixed trip timeWang et al. [18]Minimize the load factor variation and the headway variation of trains, the number of entering and existing depot operationsMINLP approachFixed passenger traveling rates, fixed trip time limitsTo enhance the passenger satisfaction and the robustness of the vehicle schedule generated, we proposed a multiobjective model of the integrated TTP and VSP with stochastic passenger flow and trip times. And a hybrid algorithm is applied to optimize the multiobjective model. The remainder of the paper is structured as follows: Section2 proposes a multiobjective optimization model for the integrated TTP and VSP, which aims at lowering the passenger travel cost, vehicle operating cost, and incompatible trip-link chain cost while satisfying the passenger and bus corporation demands. Section 3 develops a solution method to solve the integrated problem of timetabling and vehicle scheduling. Section 4 displays experimental results. Some concluding remarks and possible future research directions are given in Section 5.


## 2. Mathematical Model

The aim of the coordinated study about timetabling and vehicle scheduling is the assignment of departure intervals and the allocation of a fleet of vehicles to carry out all the trips in the timetables minimizing the costs at the same time. In this section, an integrated multiobjective optimization model of timetabling and vehicle scheduling is presented to solve the problem efficiently. This model considers the passenger waiting cost and comfort cost and scheduling scheme costs with the stochastic passenger flow and trip time.


### 2.1. Problem Description

#### 2.1.1. Solution Structure of the Integrated Problem

The integrated problem of timetabling and vehicle scheduling is concerned with the generation of a timetable and a vehicle schedule. Considering the departure intervals are different in different time periods, the whole day trips should be divided into different periods. Hence, a timetable contains a set of headways within each time period. The departure time of each trip could be calculated based on the set of headways in the timetable. An example of the timetable based on headway is presented as Figure1. Based on the timetable, the vehicle schedule is acquired with a network flow model. The network flow model is shown in Figure 2, where the letter (A or B) in the cycle before a solid arrow denotes the departure time/point of a trip, while the letter in the cycle after a solid arrow denotes the arrival time/point of a trip. The daily work of a vehicle starts from a pull-out from a depot followed by a sequence of service trips and ends at a pull-in to the depot. Each trip has to be assigned to exactly one vehicle. The connection of two consecutive trips in a vehicle is called a trip-link arc.Figure 1
The example of the timetable based on headway.Figure 2
A network flow model for the VSP (Figure2 is reproduced from [19]).A solution that contains the set of headway and allocation of the vehicle trip-link chain can refer to the scheduling scheme of the integrated problem. Suppose there aref time periods, the solution structure of the integrated problem is shown in Figure 3, where T0 is the beginning of the timetable and N is the fleet size of the vehicle schedule.Figure 3
Schematic diagram of solution structure.


#### 2.1.2. Stochastic Passenger Flow and Trip Time Problem

The passenger boarding and alighting flow would be different under different timetables. A small example of the changeable passenger flow is presented in Figure4. From the illustration, we can observe that the passenger flow varies at same time under different timetables, and the distribution of passenger flow follows some certain rules. Therefore, evaluating the cost based on a certain value of the passenger flow would be inaccurate. To estimate the passenger flow under different timetables feasibly, the passenger flow may be calculated as an integral of the distribution.Figure 4
An example of the passenger flow of same stop under different timetables.
(a)(b)Considering the stochastic trip time, the compatibility of any pair of trips becomes uncertain, and the compatibility probability is defined to determine the existence of the arcs between trips. A compatibility probability of two tripsi and j can be defined as illustrated in Figure 5, where the curve denotes the arrival time distribution fie (t) of trip i; Tie,ear and Tie,lat denote the earliest and latest arrival times of trip i, respectively; Tj denotes the scheduled departure time of trip j; and DHij denotes the deadhead time [19]. Different from the VSP with certain trip time, the compatibility probability of the pair of trips is defined as an integral of the probability distribution. Evaluating the cost of trip-link arc with uncertainty is more complicated since the incompatibility probability of the pair of trips should be penalized. Dividing the cost of VSP into the vehicle scheduling cost and the incompatible trip-link penalty cost would simplify the evaluating of the cost and make the influence of the stochastic trip time more obvious.Figure 5
Compatibility probability of any two tripsi and j (Figure 5 is reproduced from [19].).


#### 2.1.3. Multiobjective Optimization Problem

As mentioned before, there are three optimization objectives in the integrated problem, which are minimizing the passenger travel cost, minimizing the vehicle scheduling cost, and minimizing the incompatible trip-link costs. The passenger travel cost reflects the passenger waiting time spent on each stop and passenger crowding kilometers on the bus. The costs of the available vehicles and the links between trips consist of the cost of vehicle scheduling and reflect the vehicle cost and fuel consumption. The incompatible trip-link cost penalizes the infeasible trip link in the trip chain, which is attributed to the stochastic trip time. Obviously, these optimization objectives conflict with each other. A shorter headway would reduce the passenger cost but require more bus vehicles, which would increase the vehicle scheduling cost. A lower incompatibility trip-link penalty would require a longer headway and shorter trip-link chain, which would increase the passenger cost and vehicle scheduling cost. Therefore, establishing a multiobjective optimization model is a beneficial way to solve the integrated problem. The relationship between the objectives in the integrated problem is shown as Figure6.Figure 6
The conflict objectives of the integrated model.


### 2.2. A Multiobjective Optimization Model for the Integrated Problem

To formulate a mathematical model for the integrated timetabling and vehicle scheduling problem, some notations are provided in Tables2 and 3.Table 2
Subscripts and parameters of the problem formulation.
NotationsDefinitionSNumber of stopsNNumber of tripsFNumber of time periodsts,nArrival time of the tripn on stop sfsb (t)Passenger boarding flow distribution functionfso (t)Passenger alighting flow distribution functionQsbExpectation of total passenger boarding flowQsoExpectation of total passenger alighting flowCvehVehicle available costCijLink cost of arc (i, j)PijPenalty of the infeasible time of arc (i, j)TfTrip time in periodfhfminLower bound of the headway at periodfhfmaxUpper bound of the headway at periodfCmaxMaximum capacity of the bus vehiclesTable 3
Sets and variables of the problem formulation.
NotationsDefinitionPSet of pull-out arcsQSet of pull-in arcsRSet of the trip-link arcsASet of all arcs,A = P∪Q∪RDSet of depotsTSet of tripshfHeadway during the periodfxij1 if arc (i, j) is selected; 0 otherwiseThe objectives of the multiobjective optimization model are to minimize the passenger travel cost, the bus vehicle scheduling cost, and the incompatible trip-link cost. The solution of the integrated model would optimize the objectives and satisfy the constraints at the same time. In the integrated model, the decision variables are departure interval and arc link, which refer to the headway and vehicle trip-link chain, respectively. The integrated model can be expressed as(1)minZ=Z1,Z2,Z3,whereZ1=∑s=1S∑n=1Nts,n−ts,n−1⋅∫ts,n−1ts,nfsbtdt⋅Qsb+α⋅∑n=1N∑s=1Sts+1,n−ts,n⋅Fcn,s,Z2=∑d,j∈PCveh⋅xdj+∑i,j∈ACij⋅xij,Z3=∑i,j∈Rpij⋅xij,(2)s.t.ts+1,n−ts,n=TfS,s=1,…,S,n=1,…,N,(3)t1,n−t1,n−1=hf,n=1,…,N,(4)hminf≤hf≤hmaxf,(5)∑s=1s∫ts,n−1ts,nfsbtdt⋅Qsb−∫ts,n−1ts,nfsotdt⋅Qso≤Cmax,(6)∑i:i,j∈Axij−∑i:i,j∈Axji=0,∀j∈T,(7)∑i:d,j∈Axdj−∑i:i,d∈Axid=0,∀d∈D,(8)∑i:i,j∈Axij=1,∀j∈T,(9)xij∈0,1,∀i,j∈A.Formula (1) is the objective function, where Z1 is the objective function to minimize the passenger travel cost, Z2 is the objective function to minimize the vehicle scheduling cost, and Z3 is the objective function to minimize the incompatible trip-link cost. Formula (2) is the assumption of vehicle running at constant speed in one period. Formula (3) is the assumption of vehicle leaving on schedule. Formula (4) is the departure interval constraint for decision variables hf. Formula (5) is the passenger capacity constraint, which requires that the number of passengers on the bus does not exceed the capacity of the vehicle. Formulae (6) and (7) are the flow-conservation constraints. Formula (8) is the cover constraint requiring that each trip must be covered exactly by one vehicle, and formula (9) is the 0–1 constraint for trip-link decision variables xij.Considering the uncertainty, the passenger flow in functionZ1 is assigned as the time integral of passenger flow distribution to get more accurate estimation. In general, the passenger flow distribution is fitted with a certain function, such as normal distribution [20]. In objective function Z1, ∑s=1S∑n=1Nts,n−ts,n−1⋅∫ts,n−1ts,nfsbtdt⋅Qsb is the passenger waiting cost, α is a non-negative weight to adjust the penalty of crowding, and Fc (n, s) is the passenger crowding function expressed as formula (10). This formula reflects that the passenger crowding cost increases with the passenger volume on the bus at different periods.(10)Fcn,s=∑s=1s∫ts,n−1ts,nfsbtdt⋅Qsb−∫ts,n−1ts,nfsotdt⋅Qso.Objective functionZ2 is defined as the sum of the costs of the available vehicles and the links between trips. The costs of the available vehicles is assigned as ∑d,j∈PCveh⋅xdj, in which xdj is defined as the pull-out trip links. The costs of the links between trips is expressed as ∑i,j∈ACij⋅xij.Objective functionZ3 reflects the gap between the theoretical trip time and the realistic trip time, which would reduce the feasibility of the scheduling scheme. As suggested by Shen et al. [20], Pij in function Z3 denotes the penalty of the infeasible time of arc (i, j) and can be defined as(11)Pij=EIFij2=∫Tj−DHij+∞tie−Tj+DHij2fietdt=∫Tj−Ti−DHij+∞Ti+t−Tj+DHij2fitdt,where Ti and Tj denote the scheduled departure time of trip i and j, respectively; DHij denotes the deadhead time from arrival point of trip i to departure point of trip j; IFij=tie−Tj+DHij is the infeasible time of an arc (i, j); and fit is the probability density function about the duration of trip i.


## 2.1. Problem Description

### 2.1.1. Solution Structure of the Integrated Problem

The integrated problem of timetabling and vehicle scheduling is concerned with the generation of a timetable and a vehicle schedule. Considering the departure intervals are different in different time periods, the whole day trips should be divided into different periods. Hence, a timetable contains a set of headways within each time period. The departure time of each trip could be calculated based on the set of headways in the timetable. An example of the timetable based on headway is presented as Figure1. Based on the timetable, the vehicle schedule is acquired with a network flow model. The network flow model is shown in Figure 2, where the letter (A or B) in the cycle before a solid arrow denotes the departure time/point of a trip, while the letter in the cycle after a solid arrow denotes the arrival time/point of a trip. The daily work of a vehicle starts from a pull-out from a depot followed by a sequence of service trips and ends at a pull-in to the depot. Each trip has to be assigned to exactly one vehicle. The connection of two consecutive trips in a vehicle is called a trip-link arc.Figure 1
The example of the timetable based on headway.Figure 2
A network flow model for the VSP (Figure2 is reproduced from [19]).A solution that contains the set of headway and allocation of the vehicle trip-link chain can refer to the scheduling scheme of the integrated problem. Suppose there aref time periods, the solution structure of the integrated problem is shown in Figure 3, where T0 is the beginning of the timetable and N is the fleet size of the vehicle schedule.Figure 3
Schematic diagram of solution structure.


### 2.1.2. Stochastic Passenger Flow and Trip Time Problem

The passenger boarding and alighting flow would be different under different timetables. A small example of the changeable passenger flow is presented in Figure4. From the illustration, we can observe that the passenger flow varies at same time under different timetables, and the distribution of passenger flow follows some certain rules. Therefore, evaluating the cost based on a certain value of the passenger flow would be inaccurate. To estimate the passenger flow under different timetables feasibly, the passenger flow may be calculated as an integral of the distribution.Figure 4
An example of the passenger flow of same stop under different timetables.
(a)(b)Considering the stochastic trip time, the compatibility of any pair of trips becomes uncertain, and the compatibility probability is defined to determine the existence of the arcs between trips. A compatibility probability of two tripsi and j can be defined as illustrated in Figure 5, where the curve denotes the arrival time distribution fie (t) of trip i; Tie,ear and Tie,lat denote the earliest and latest arrival times of trip i, respectively; Tj denotes the scheduled departure time of trip j; and DHij denotes the deadhead time [19]. Different from the VSP with certain trip time, the compatibility probability of the pair of trips is defined as an integral of the probability distribution. Evaluating the cost of trip-link arc with uncertainty is more complicated since the incompatibility probability of the pair of trips should be penalized. Dividing the cost of VSP into the vehicle scheduling cost and the incompatible trip-link penalty cost would simplify the evaluating of the cost and make the influence of the stochastic trip time more obvious.Figure 5
Compatibility probability of any two tripsi and j (Figure 5 is reproduced from [19].).


### 2.1.3. Multiobjective Optimization Problem

As mentioned before, there are three optimization objectives in the integrated problem, which are minimizing the passenger travel cost, minimizing the vehicle scheduling cost, and minimizing the incompatible trip-link costs. The passenger travel cost reflects the passenger waiting time spent on each stop and passenger crowding kilometers on the bus. The costs of the available vehicles and the links between trips consist of the cost of vehicle scheduling and reflect the vehicle cost and fuel consumption. The incompatible trip-link cost penalizes the infeasible trip link in the trip chain, which is attributed to the stochastic trip time. Obviously, these optimization objectives conflict with each other. A shorter headway would reduce the passenger cost but require more bus vehicles, which would increase the vehicle scheduling cost. A lower incompatibility trip-link penalty would require a longer headway and shorter trip-link chain, which would increase the passenger cost and vehicle scheduling cost. Therefore, establishing a multiobjective optimization model is a beneficial way to solve the integrated problem. The relationship between the objectives in the integrated problem is shown as Figure6.Figure 6
The conflict objectives of the integrated model.


## 2.1.1. Solution Structure of the Integrated Problem

The integrated problem of timetabling and vehicle scheduling is concerned with the generation of a timetable and a vehicle schedule. Considering the departure intervals are different in different time periods, the whole day trips should be divided into different periods. Hence, a timetable contains a set of headways within each time period. The departure time of each trip could be calculated based on the set of headways in the timetable. An example of the timetable based on headway is presented as Figure1. Based on the timetable, the vehicle schedule is acquired with a network flow model. The network flow model is shown in Figure 2, where the letter (A or B) in the cycle before a solid arrow denotes the departure time/point of a trip, while the letter in the cycle after a solid arrow denotes the arrival time/point of a trip. The daily work of a vehicle starts from a pull-out from a depot followed by a sequence of service trips and ends at a pull-in to the depot. Each trip has to be assigned to exactly one vehicle. The connection of two consecutive trips in a vehicle is called a trip-link arc.Figure 1
The example of the timetable based on headway.Figure 2
A network flow model for the VSP (Figure2 is reproduced from [19]).A solution that contains the set of headway and allocation of the vehicle trip-link chain can refer to the scheduling scheme of the integrated problem. Suppose there aref time periods, the solution structure of the integrated problem is shown in Figure 3, where T0 is the beginning of the timetable and N is the fleet size of the vehicle schedule.Figure 3
Schematic diagram of solution structure.


## 2.1.2. Stochastic Passenger Flow and Trip Time Problem

The passenger boarding and alighting flow would be different under different timetables. A small example of the changeable passenger flow is presented in Figure4. From the illustration, we can observe that the passenger flow varies at same time under different timetables, and the distribution of passenger flow follows some certain rules. Therefore, evaluating the cost based on a certain value of the passenger flow would be inaccurate. To estimate the passenger flow under different timetables feasibly, the passenger flow may be calculated as an integral of the distribution.Figure 4
An example of the passenger flow of same stop under different timetables.
(a)(b)Considering the stochastic trip time, the compatibility of any pair of trips becomes uncertain, and the compatibility probability is defined to determine the existence of the arcs between trips. A compatibility probability of two tripsi and j can be defined as illustrated in Figure 5, where the curve denotes the arrival time distribution fie (t) of trip i; Tie,ear and Tie,lat denote the earliest and latest arrival times of trip i, respectively; Tj denotes the scheduled departure time of trip j; and DHij denotes the deadhead time [19]. Different from the VSP with certain trip time, the compatibility probability of the pair of trips is defined as an integral of the probability distribution. Evaluating the cost of trip-link arc with uncertainty is more complicated since the incompatibility probability of the pair of trips should be penalized. Dividing the cost of VSP into the vehicle scheduling cost and the incompatible trip-link penalty cost would simplify the evaluating of the cost and make the influence of the stochastic trip time more obvious.Figure 5
Compatibility probability of any two tripsi and j (Figure 5 is reproduced from [19].).


## 2.1.3. Multiobjective Optimization Problem

As mentioned before, there are three optimization objectives in the integrated problem, which are minimizing the passenger travel cost, minimizing the vehicle scheduling cost, and minimizing the incompatible trip-link costs. The passenger travel cost reflects the passenger waiting time spent on each stop and passenger crowding kilometers on the bus. The costs of the available vehicles and the links between trips consist of the cost of vehicle scheduling and reflect the vehicle cost and fuel consumption. The incompatible trip-link cost penalizes the infeasible trip link in the trip chain, which is attributed to the stochastic trip time. Obviously, these optimization objectives conflict with each other. A shorter headway would reduce the passenger cost but require more bus vehicles, which would increase the vehicle scheduling cost. A lower incompatibility trip-link penalty would require a longer headway and shorter trip-link chain, which would increase the passenger cost and vehicle scheduling cost. Therefore, establishing a multiobjective optimization model is a beneficial way to solve the integrated problem. The relationship between the objectives in the integrated problem is shown as Figure6.Figure 6
The conflict objectives of the integrated model.


## 2.2. A Multiobjective Optimization Model for the Integrated Problem

To formulate a mathematical model for the integrated timetabling and vehicle scheduling problem, some notations are provided in Tables2 and 3.Table 2
Subscripts and parameters of the problem formulation.
NotationsDefinitionSNumber of stopsNNumber of tripsFNumber of time periodsts,nArrival time of the tripn on stop sfsb (t)Passenger boarding flow distribution functionfso (t)Passenger alighting flow distribution functionQsbExpectation of total passenger boarding flowQsoExpectation of total passenger alighting flowCvehVehicle available costCijLink cost of arc (i, j)PijPenalty of the infeasible time of arc (i, j)TfTrip time in periodfhfminLower bound of the headway at periodfhfmaxUpper bound of the headway at periodfCmaxMaximum capacity of the bus vehiclesTable 3
Sets and variables of the problem formulation.
NotationsDefinitionPSet of pull-out arcsQSet of pull-in arcsRSet of the trip-link arcsASet of all arcs,A = P∪Q∪RDSet of depotsTSet of tripshfHeadway during the periodfxij1 if arc (i, j) is selected; 0 otherwiseThe objectives of the multiobjective optimization model are to minimize the passenger travel cost, the bus vehicle scheduling cost, and the incompatible trip-link cost. The solution of the integrated model would optimize the objectives and satisfy the constraints at the same time. In the integrated model, the decision variables are departure interval and arc link, which refer to the headway and vehicle trip-link chain, respectively. The integrated model can be expressed as(1)minZ=Z1,Z2,Z3,whereZ1=∑s=1S∑n=1Nts,n−ts,n−1⋅∫ts,n−1ts,nfsbtdt⋅Qsb+α⋅∑n=1N∑s=1Sts+1,n−ts,n⋅Fcn,s,Z2=∑d,j∈PCveh⋅xdj+∑i,j∈ACij⋅xij,Z3=∑i,j∈Rpij⋅xij,(2)s.t.ts+1,n−ts,n=TfS,s=1,…,S,n=1,…,N,(3)t1,n−t1,n−1=hf,n=1,…,N,(4)hminf≤hf≤hmaxf,(5)∑s=1s∫ts,n−1ts,nfsbtdt⋅Qsb−∫ts,n−1ts,nfsotdt⋅Qso≤Cmax,(6)∑i:i,j∈Axij−∑i:i,j∈Axji=0,∀j∈T,(7)∑i:d,j∈Axdj−∑i:i,d∈Axid=0,∀d∈D,(8)∑i:i,j∈Axij=1,∀j∈T,(9)xij∈0,1,∀i,j∈A.Formula (1) is the objective function, where Z1 is the objective function to minimize the passenger travel cost, Z2 is the objective function to minimize the vehicle scheduling cost, and Z3 is the objective function to minimize the incompatible trip-link cost. Formula (2) is the assumption of vehicle running at constant speed in one period. Formula (3) is the assumption of vehicle leaving on schedule. Formula (4) is the departure interval constraint for decision variables hf. Formula (5) is the passenger capacity constraint, which requires that the number of passengers on the bus does not exceed the capacity of the vehicle. Formulae (6) and (7) are the flow-conservation constraints. Formula (8) is the cover constraint requiring that each trip must be covered exactly by one vehicle, and formula (9) is the 0–1 constraint for trip-link decision variables xij.Considering the uncertainty, the passenger flow in functionZ1 is assigned as the time integral of passenger flow distribution to get more accurate estimation. In general, the passenger flow distribution is fitted with a certain function, such as normal distribution [20]. In objective function Z1, ∑s=1S∑n=1Nts,n−ts,n−1⋅∫ts,n−1ts,nfsbtdt⋅Qsb is the passenger waiting cost, α is a non-negative weight to adjust the penalty of crowding, and Fc (n, s) is the passenger crowding function expressed as formula (10). This formula reflects that the passenger crowding cost increases with the passenger volume on the bus at different periods.(10)Fcn,s=∑s=1s∫ts,n−1ts,nfsbtdt⋅Qsb−∫ts,n−1ts,nfsotdt⋅Qso.Objective functionZ2 is defined as the sum of the costs of the available vehicles and the links between trips. The costs of the available vehicles is assigned as ∑d,j∈PCveh⋅xdj, in which xdj is defined as the pull-out trip links. The costs of the links between trips is expressed as ∑i,j∈ACij⋅xij.Objective functionZ3 reflects the gap between the theoretical trip time and the realistic trip time, which would reduce the feasibility of the scheduling scheme. As suggested by Shen et al. [20], Pij in function Z3 denotes the penalty of the infeasible time of arc (i, j) and can be defined as(11)Pij=EIFij2=∫Tj−DHij+∞tie−Tj+DHij2fietdt=∫Tj−Ti−DHij+∞Ti+t−Tj+DHij2fitdt,where Ti and Tj denote the scheduled departure time of trip i and j, respectively; DHij denotes the deadhead time from arrival point of trip i to departure point of trip j; IFij=tie−Tj+DHij is the infeasible time of an arc (i, j); and fit is the probability density function about the duration of trip i.


## 3. The Multiobjective Hybrid Algorithm

Due to the NP-hard nature of the multiobjective optimization problem, the exact approach is difficult to solve the problem. Therefore, we develop a multiobjective hybrid algorithm to optimize the multiobjective integrated optimization problem. The algorithm is a combination of a self-adjust genetic algorithm (SGA), large neighborhood search (LNS) algorithm, and Pareto separation operator (PSO). The SGA is applied to construct a new headway, the LNS is employed to find the optimal trip-link chain under the headway, and the PSO is used to select the multiobjective optimal solution. The solution obtained by the hybrid algorithm is the scheduling scheme that consists of the headway and trip-link chain. Basic steps of the multiobjective hybrid algorithm are shown in Figure7. Some definitions and the constructive approaches for the solution structure, the SGA, the LNS, and the PSO are elaborated in Sections 3.1–3.4.Figure 7
The framework of the multiobjective hybrid algorithm.


### 3.1. The Structure of Solution Set

As mentioned in Section2, the solution of the multiobjective model consists of two parts: headway vector Tg = (t1, t2,…, tn) and trip chain set Dg = {X1, X2,…, XM}. The vector element tn in Tg represents the headway value of the n-th characteristic periods. The element XM in Dg represents a feasible trip chain operated by vehicle M and can be expressed as XM = (xi, …, xj), in which xi denotes the i-th trip in the trip set.For example, suppose there are three trips and two characteristic periods in the model. The headway vectorTg = (10, 8) refers to the headway being 10 minutes in period 1 and 8 minutes in period 2, respectively. The trip-link chain set Dg\ = {X1} refers to that vehicle 1 is to carry out all the trips based on Tg, and the trip-link chain set element X1 = (1, 3, 2) refers to the trips 1 and 3 and 3 and 2 that are linked as one chain.


### 3.2. The Self-Adjust GA for Headway Vector

#### 3.2.1. Basic Frame of the MethodG (s)

The GA is the most commonly used metaheuristic algorithm and performs well in operation research fields such as crew scheduling [21] and vehicle scheduling [22]. The GA search is based on a population of solutions and some operators used to perform the evolution in search process. The main operators in this method are selection, crossover, and mutation. The solution in GA is usually coded in finite length string, which is called chromosome. Starting from a population of chromosomes, all the chromosomes are evaluated using fitness function to rank the individuals. The crossover operator and mutation operator are responsible for the creation of new individuals.In our self-adjust genetic algorithmG (s), a chromosome represents a headway vector Tg. The new headway vector Tg at each iteration g is acquired from a pair of parent chromosomes, which are selected by the roulette wheel. Each gene in the chromosome Tg represents a headway value in one period. The process of the crossover operator and mutation operator is illustrated as Figure 8. The location of the chromosome in the crossover operator and mutation operator are selected randomly at each iteration. The offspring individual created in each iteration should satisfy the headway constraints in the model to ensure the solutions obtained from the algorithm are feasible in the real world. The basic steps of the G (s) are shown in Table 4.Figure 8
The crossover and mutation process of methodG (S).Table 4
Framework and steps of methodG (s).
MethodG (s)1Initialize population setPg=Ts1,Ts2,Ts32Forg=1toG3Select a pair of parent chromosomesTi, Tj from Pg using roulette wheel4Create a new chromosomeTg using crossover and mutation operators based on parent chromosomes5IfTg satisfy the headway constraints, go to step 6; otherwise, go to step 36AddTg to the population set as Pg+1=Pg,Tg7End for8Return the headway vector setPg=Ts1,Ts2,Ts3,T1,…,Tn


#### 3.2.2. Self-Adjust Strategy ofG(s)

Different from the normal GA approach in which population size is assigned as a fixed value on initialization, the population size in our method is adjusted with iterationg. That is, the population size is initially defined as three and is adjusted with the new feasible solution added. In general, the median value of the headway constraints is supposed to be the expected headway. Therefore, we define that there are three individuals in the initial population set Pg, one is the headway vector consisting of the median value of the headway constraints in each period, and others are the headway vectors plus one and reduce one, respectively, based on the median value headway vector. That is, the three initial chromosomes are Ts1 = (t1, …, tn), Ts2 = (t1 + 1, …, tn + 1), and Ts3 = (t1 − 1, …, tn − 1), where tf = (hfmin + hfmax)/2, and hfmin and hfmax are the lower and upper bounds of the headway in period f. After each iteration, the new headway vector is added to the population set, and the population size is enlarged by one.The parent chromosomes are selected using the roulette wheel at each iteration. Different from the normal roulette wheel method in which the selection probability of each chromosome in the population is calculated by fitness function, the selection probability of the three initial individuals is adjusted with the population size change. The selection probability of the initial individuals is calculated as1/n, where n is the number of individuals in the population set at each iteration. The fitness function value in the roulette wheel is the sum values of three objectives in formula (1).


### 3.3. The LNS Method for Trip-Link Chain

The LNS is originally proposed by Shaw [23] and has been successfully applied to the VSP and some extensions [24]. The basic idea of the LNS is to search in large neighborhoods, which may contain more potentially better solutions compared with small neighborhoods. The search neighborhood of the solution is defined implicitly by a destroy method and a repair method. A destroy method disrupts part of the current solution, while a repair method rebuilds the destroyed solution. In the LNS, one or more destroy and repair methods would be employed. And in each iteration, the destroy and repair methods to be applied are selected according to their weights, which are defined by users before the iterations. The resultant solution is accepted based on an acceptance criterion, and the heuristic stops when the stop criterion is met.The solution search of methodL (S) begins from an initial trip chain set Dg. This trip set is received by a greedy constructive heuristic method based on the timetable acquired from headway Tg. First, the trips are sorted in a nondecreasing order of their start times. Then, these trips are added to the solution during iteration. In each iteration, the selected trip is inserted to the end of the existing vehicle tasks that is compatible in time and leads to the minimum increase in the objective value of vehicle scheduling cost. If the trip is unfeasible to be added to any existing trip-link arc, it would be set as a new trip-link chain. The procedure stops until all the trips are served. The basic steps of the LNS method L (S) are shown in Table 5.Table 5
Framework and steps of methodL (S).
MethodL (S)1Initialize trip chain setDg=X1,…,XM. Define the destroy operator set α−=α1−,…,αi− and repair operator set β+=β1+,…,βj+2Forl=1toL3Select a destroy methodαi− and a repair method βj+ randomly4Removeq trips from Dg using αi−.5Insert the removed trips onDg usingβj+7End for8Return the final solutionDgThe operators used in methodL (S) at each iteration are selected randomly. The iteration of method L (s) is to search trip link combination in solution neighborhood. The destroy operator is to remove some trips, and the repair operator is to insert these trips back to find new trip links. The definitions of the destroy operators and repair operators are described in the following sections.


#### 3.3.1. Destroy Operators

A selection of destroy operators are implemented in our LNS, including random remove operator, short trip-link remove operator, and incompatibility trip-link remove operator. These removal operators are given different weights before iterationL and selected by the roulette wheel. The removal methods used in our LNS are defined as follows:(1)
Random remove operator. It simply removes a number of trips randomly from the solution to diversify the search into different parts of the solution space.(2)
Short trip-link remove operator. It would remove the shortest trip-link chain to reduce the vehicle scheduling costs.(3)
Incompatibility trip-link remove operator. It would remove some maximal incompatibility penalty trip links from the solution to reduce the penalty costs.


#### 3.3.2. Repair Operators

A class of regret insertion heuristics is implemented as the repair operators in methodL (S). These operators insert one trip to the best position at a time. They differ in their ways of selecting the trip to be inserted [25]. These approaches could be defined as below:LetNr denote the set of trips to be inserted, and Δfi,p denotes the cost of inserting trip i∈Nr to position p in the current partial solution. If trip i is inserted to position p, Δfi,p=cprep,i+ci,sucp−cprep,sucp, where pre (p) and suc (p) are the predecessor trip and successor trip corresponding to position p, respectively, and ci,j denotes the cost of trip-link arc (i, j). If the insertion leads to an infeasible solution, we set Δfi,p=∞.Letci1 denote the insertion cost of trip i at its best position, ci1=minp∈Pcip, where P is the set of all the positions, and cik is the insertion cost at the k-th best position. In a regret-k heuristic, the trip i∗=argmaxi∈Nr∑j=1j=kcij−ci1 is selected as the next trip to be inserted.In our LNS, the regret-2, regret-3, and regret-4 heuristic methods are implemented and selected randomly with equal probability at each iteration. The trip insertion cost is calculated as the weighted sum of the vehicle scheduling cost and the penalty cost in formula (1).


### 3.4. The PSO to Simplify Computational Complexity

#### 3.4.1. Pareto Dominance Relation

Multiobjective optimization deals with optimizing multiple objectives. As mentioned in Section2, the integrated model formulates a problem with three objectives. Different from the optimization with one objective, a solution cannot be compared with another using relational operators in a multiobjective search space because objectives can be in conflict and require special considerations. This is due to the existence of more than one criterion for comparison [26]. Therefore, we need other operators to measure and find out how much a solution is better than another. The Pareto dominance is the most widely used operator to measure the solution of multiobjective model, and a series of multiobjective optimization algorithms are developed based on it [27]. The operator can be mathematically defined in a minimization optimization problem as follows:(12)∀i∈1,2,…,k:fix⟶≤fiy⟶∧∃i∈1,2,…,k:fix⟶≤fiy⟶,where x⟶=x1,x2,…,xk and y⟶=y1,y2,…,yk. This equation shows that a solution (vector x) is better than another (vector y) if it has equal and at least one better value on all objectives and can be denoted as x≺y. Suppose there are two objective functions f1 and f2 in a minimization optimization problem, the dominance relation could be presented from an example in Figure 9. This figure shows that the circles are better than squares because they provide a lower value in both objectives.Figure 9
The Pareto dominance relation.For the multiobjective optimization problem, there is a set of nondominated solutions. Consequently, the projection of Pareto optimal solutions in the objective space are stored in a set called Pareto optimal front. The optimization of multiobjective problem aims to find the Pareto optimal front of the problem.


#### 3.4.2. Pareto Separation Operator

The Pareto dominance relation is applied in our multiobjective hybrid method to find the optimal operation scheme of the integrated model. Suppose there arem objectives and n acquired solutions in the multiobjective optimization problem, the number of comparison is Cn2⋅m in general, which would increase rapidly with the addition of the objectives and solutions. In order to reduce the number of comparison in our multiobjective algorithm and find the Pareto optimal front quickly, a separation solution comparison method which is called the Pareto separation operator is applied. The process of the PSO is presented as below.It is easy to find that if one solution is nondominated in the multiobjective problem, it would be nondominated in at least one subproblem of the multiobjective optimization problem. That is, if one solution is nondominated in a two-objective subproblem of the multiobjective model, it would be a Pareto optimal solution. The relation between the multiobjective problem and the subproblem is illustrated in Figure10, in which the point in the Pareto front is also the nondominated point in the projection. Therefore, in our comparison method, we only need to find the nondominated solutions in the two-objective subproblem of the multiobjective model, and the union of these solutions is the Pareto optimal solution set.Figure 10
The relation between multiobjective problem and subproblem.The comparison between two objectives is easier than the comparison between more objectives. Suppose there aren solutions in the two-objective subproblem, after we find the optimal solution Sm on the objective 1, we just need to find the solution better than the solution Sm on the objective 2, and these solutions are the Pareto optimal solution of the two-objective subproblem. The number of comparison in our method would be no more than Cm2⋅2n−1, where Cm2 denotes the number of the two-objective subproblem. It is easy to prove that, if n is larger than 2 (m −1), Cn2⋅m would be larger than Cm2⋅2n−1. In the most multiobjective problem, the number of solutions n is far larger than the number of objectives m. Hence, our method would effectively reduce the computational complexity. An example of the comparison method is shown in Figure 11.Figure 11
The example of the Pareto separation operator.


## 3.1. The Structure of Solution Set

As mentioned in Section2, the solution of the multiobjective model consists of two parts: headway vector Tg = (t1, t2,…, tn) and trip chain set Dg = {X1, X2,…, XM}. The vector element tn in Tg represents the headway value of the n-th characteristic periods. The element XM in Dg represents a feasible trip chain operated by vehicle M and can be expressed as XM = (xi, …, xj), in which xi denotes the i-th trip in the trip set.For example, suppose there are three trips and two characteristic periods in the model. The headway vectorTg = (10, 8) refers to the headway being 10 minutes in period 1 and 8 minutes in period 2, respectively. The trip-link chain set Dg\ = {X1} refers to that vehicle 1 is to carry out all the trips based on Tg, and the trip-link chain set element X1 = (1, 3, 2) refers to the trips 1 and 3 and 3 and 2 that are linked as one chain.


## 3.2. The Self-Adjust GA for Headway Vector

### 3.2.1. Basic Frame of the MethodG (s)

The GA is the most commonly used metaheuristic algorithm and performs well in operation research fields such as crew scheduling [21] and vehicle scheduling [22]. The GA search is based on a population of solutions and some operators used to perform the evolution in search process. The main operators in this method are selection, crossover, and mutation. The solution in GA is usually coded in finite length string, which is called chromosome. Starting from a population of chromosomes, all the chromosomes are evaluated using fitness function to rank the individuals. The crossover operator and mutation operator are responsible for the creation of new individuals.In our self-adjust genetic algorithmG (s), a chromosome represents a headway vector Tg. The new headway vector Tg at each iteration g is acquired from a pair of parent chromosomes, which are selected by the roulette wheel. Each gene in the chromosome Tg represents a headway value in one period. The process of the crossover operator and mutation operator is illustrated as Figure 8. The location of the chromosome in the crossover operator and mutation operator are selected randomly at each iteration. The offspring individual created in each iteration should satisfy the headway constraints in the model to ensure the solutions obtained from the algorithm are feasible in the real world. The basic steps of the G (s) are shown in Table 4.Figure 8
The crossover and mutation process of methodG (S).Table 4
Framework and steps of methodG (s).
MethodG (s)1Initialize population setPg=Ts1,Ts2,Ts32Forg=1toG3Select a pair of parent chromosomesTi, Tj from Pg using roulette wheel4Create a new chromosomeTg using crossover and mutation operators based on parent chromosomes5IfTg satisfy the headway constraints, go to step 6; otherwise, go to step 36AddTg to the population set as Pg+1=Pg,Tg7End for8Return the headway vector setPg=Ts1,Ts2,Ts3,T1,…,Tn


### 3.2.2. Self-Adjust Strategy ofG(s)

Different from the normal GA approach in which population size is assigned as a fixed value on initialization, the population size in our method is adjusted with iterationg. That is, the population size is initially defined as three and is adjusted with the new feasible solution added. In general, the median value of the headway constraints is supposed to be the expected headway. Therefore, we define that there are three individuals in the initial population set Pg, one is the headway vector consisting of the median value of the headway constraints in each period, and others are the headway vectors plus one and reduce one, respectively, based on the median value headway vector. That is, the three initial chromosomes are Ts1 = (t1, …, tn), Ts2 = (t1 + 1, …, tn + 1), and Ts3 = (t1 − 1, …, tn − 1), where tf = (hfmin + hfmax)/2, and hfmin and hfmax are the lower and upper bounds of the headway in period f. After each iteration, the new headway vector is added to the population set, and the population size is enlarged by one.The parent chromosomes are selected using the roulette wheel at each iteration. Different from the normal roulette wheel method in which the selection probability of each chromosome in the population is calculated by fitness function, the selection probability of the three initial individuals is adjusted with the population size change. The selection probability of the initial individuals is calculated as1/n, where n is the number of individuals in the population set at each iteration. The fitness function value in the roulette wheel is the sum values of three objectives in formula (1).


## 3.2.1. Basic Frame of the MethodG (s)

The GA is the most commonly used metaheuristic algorithm and performs well in operation research fields such as crew scheduling [21] and vehicle scheduling [22]. The GA search is based on a population of solutions and some operators used to perform the evolution in search process. The main operators in this method are selection, crossover, and mutation. The solution in GA is usually coded in finite length string, which is called chromosome. Starting from a population of chromosomes, all the chromosomes are evaluated using fitness function to rank the individuals. The crossover operator and mutation operator are responsible for the creation of new individuals.In our self-adjust genetic algorithmG (s), a chromosome represents a headway vector Tg. The new headway vector Tg at each iteration g is acquired from a pair of parent chromosomes, which are selected by the roulette wheel. Each gene in the chromosome Tg represents a headway value in one period. The process of the crossover operator and mutation operator is illustrated as Figure 8. The location of the chromosome in the crossover operator and mutation operator are selected randomly at each iteration. The offspring individual created in each iteration should satisfy the headway constraints in the model to ensure the solutions obtained from the algorithm are feasible in the real world. The basic steps of the G (s) are shown in Table 4.Figure 8
The crossover and mutation process of methodG (S).Table 4
Framework and steps of methodG (s).
MethodG (s)1Initialize population setPg=Ts1,Ts2,Ts32Forg=1toG3Select a pair of parent chromosomesTi, Tj from Pg using roulette wheel4Create a new chromosomeTg using crossover and mutation operators based on parent chromosomes5IfTg satisfy the headway constraints, go to step 6; otherwise, go to step 36AddTg to the population set as Pg+1=Pg,Tg7End for8Return the headway vector setPg=Ts1,Ts2,Ts3,T1,…,Tn


## 3.2.2. Self-Adjust Strategy ofG(s)

Different from the normal GA approach in which population size is assigned as a fixed value on initialization, the population size in our method is adjusted with iterationg. That is, the population size is initially defined as three and is adjusted with the new feasible solution added. In general, the median value of the headway constraints is supposed to be the expected headway. Therefore, we define that there are three individuals in the initial population set Pg, one is the headway vector consisting of the median value of the headway constraints in each period, and others are the headway vectors plus one and reduce one, respectively, based on the median value headway vector. That is, the three initial chromosomes are Ts1 = (t1, …, tn), Ts2 = (t1 + 1, …, tn + 1), and Ts3 = (t1 − 1, …, tn − 1), where tf = (hfmin + hfmax)/2, and hfmin and hfmax are the lower and upper bounds of the headway in period f. After each iteration, the new headway vector is added to the population set, and the population size is enlarged by one.The parent chromosomes are selected using the roulette wheel at each iteration. Different from the normal roulette wheel method in which the selection probability of each chromosome in the population is calculated by fitness function, the selection probability of the three initial individuals is adjusted with the population size change. The selection probability of the initial individuals is calculated as1/n, where n is the number of individuals in the population set at each iteration. The fitness function value in the roulette wheel is the sum values of three objectives in formula (1).


## 3.3. The LNS Method for Trip-Link Chain

The LNS is originally proposed by Shaw [23] and has been successfully applied to the VSP and some extensions [24]. The basic idea of the LNS is to search in large neighborhoods, which may contain more potentially better solutions compared with small neighborhoods. The search neighborhood of the solution is defined implicitly by a destroy method and a repair method. A destroy method disrupts part of the current solution, while a repair method rebuilds the destroyed solution. In the LNS, one or more destroy and repair methods would be employed. And in each iteration, the destroy and repair methods to be applied are selected according to their weights, which are defined by users before the iterations. The resultant solution is accepted based on an acceptance criterion, and the heuristic stops when the stop criterion is met.The solution search of methodL (S) begins from an initial trip chain set Dg. This trip set is received by a greedy constructive heuristic method based on the timetable acquired from headway Tg. First, the trips are sorted in a nondecreasing order of their start times. Then, these trips are added to the solution during iteration. In each iteration, the selected trip is inserted to the end of the existing vehicle tasks that is compatible in time and leads to the minimum increase in the objective value of vehicle scheduling cost. If the trip is unfeasible to be added to any existing trip-link arc, it would be set as a new trip-link chain. The procedure stops until all the trips are served. The basic steps of the LNS method L (S) are shown in Table 5.Table 5
Framework and steps of methodL (S).
MethodL (S)1Initialize trip chain setDg=X1,…,XM. Define the destroy operator set α−=α1−,…,αi− and repair operator set β+=β1+,…,βj+2Forl=1toL3Select a destroy methodαi− and a repair method βj+ randomly4Removeq trips from Dg using αi−.5Insert the removed trips onDg usingβj+7End for8Return the final solutionDgThe operators used in methodL (S) at each iteration are selected randomly. The iteration of method L (s) is to search trip link combination in solution neighborhood. The destroy operator is to remove some trips, and the repair operator is to insert these trips back to find new trip links. The definitions of the destroy operators and repair operators are described in the following sections.


### 3.3.1. Destroy Operators

A selection of destroy operators are implemented in our LNS, including random remove operator, short trip-link remove operator, and incompatibility trip-link remove operator. These removal operators are given different weights before iterationL and selected by the roulette wheel. The removal methods used in our LNS are defined as follows:(1)
Random remove operator. It simply removes a number of trips randomly from the solution to diversify the search into different parts of the solution space.(2)
Short trip-link remove operator. It would remove the shortest trip-link chain to reduce the vehicle scheduling costs.(3)
Incompatibility trip-link remove operator. It would remove some maximal incompatibility penalty trip links from the solution to reduce the penalty costs.


### 3.3.2. Repair Operators

A class of regret insertion heuristics is implemented as the repair operators in methodL (S). These operators insert one trip to the best position at a time. They differ in their ways of selecting the trip to be inserted [25]. These approaches could be defined as below:LetNr denote the set of trips to be inserted, and Δfi,p denotes the cost of inserting trip i∈Nr to position p in the current partial solution. If trip i is inserted to position p, Δfi,p=cprep,i+ci,sucp−cprep,sucp, where pre (p) and suc (p) are the predecessor trip and successor trip corresponding to position p, respectively, and ci,j denotes the cost of trip-link arc (i, j). If the insertion leads to an infeasible solution, we set Δfi,p=∞.Letci1 denote the insertion cost of trip i at its best position, ci1=minp∈Pcip, where P is the set of all the positions, and cik is the insertion cost at the k-th best position. In a regret-k heuristic, the trip i∗=argmaxi∈Nr∑j=1j=kcij−ci1 is selected as the next trip to be inserted.In our LNS, the regret-2, regret-3, and regret-4 heuristic methods are implemented and selected randomly with equal probability at each iteration. The trip insertion cost is calculated as the weighted sum of the vehicle scheduling cost and the penalty cost in formula (1).


## 3.3.1. Destroy Operators

A selection of destroy operators are implemented in our LNS, including random remove operator, short trip-link remove operator, and incompatibility trip-link remove operator. These removal operators are given different weights before iterationL and selected by the roulette wheel. The removal methods used in our LNS are defined as follows:(1)
Random remove operator. It simply removes a number of trips randomly from the solution to diversify the search into different parts of the solution space.(2)
Short trip-link remove operator. It would remove the shortest trip-link chain to reduce the vehicle scheduling costs.(3)
Incompatibility trip-link remove operator. It would remove some maximal incompatibility penalty trip links from the solution to reduce the penalty costs.


## 3.3.2. Repair Operators

A class of regret insertion heuristics is implemented as the repair operators in methodL (S). These operators insert one trip to the best position at a time. They differ in their ways of selecting the trip to be inserted [25]. These approaches could be defined as below:LetNr denote the set of trips to be inserted, and Δfi,p denotes the cost of inserting trip i∈Nr to position p in the current partial solution. If trip i is inserted to position p, Δfi,p=cprep,i+ci,sucp−cprep,sucp, where pre (p) and suc (p) are the predecessor trip and successor trip corresponding to position p, respectively, and ci,j denotes the cost of trip-link arc (i, j). If the insertion leads to an infeasible solution, we set Δfi,p=∞.Letci1 denote the insertion cost of trip i at its best position, ci1=minp∈Pcip, where P is the set of all the positions, and cik is the insertion cost at the k-th best position. In a regret-k heuristic, the trip i∗=argmaxi∈Nr∑j=1j=kcij−ci1 is selected as the next trip to be inserted.In our LNS, the regret-2, regret-3, and regret-4 heuristic methods are implemented and selected randomly with equal probability at each iteration. The trip insertion cost is calculated as the weighted sum of the vehicle scheduling cost and the penalty cost in formula (1).


## 3.4. The PSO to Simplify Computational Complexity

### 3.4.1. Pareto Dominance Relation

Multiobjective optimization deals with optimizing multiple objectives. As mentioned in Section2, the integrated model formulates a problem with three objectives. Different from the optimization with one objective, a solution cannot be compared with another using relational operators in a multiobjective search space because objectives can be in conflict and require special considerations. This is due to the existence of more than one criterion for comparison [26]. Therefore, we need other operators to measure and find out how much a solution is better than another. The Pareto dominance is the most widely used operator to measure the solution of multiobjective model, and a series of multiobjective optimization algorithms are developed based on it [27]. The operator can be mathematically defined in a minimization optimization problem as follows:(12)∀i∈1,2,…,k:fix⟶≤fiy⟶∧∃i∈1,2,…,k:fix⟶≤fiy⟶,where x⟶=x1,x2,…,xk and y⟶=y1,y2,…,yk. This equation shows that a solution (vector x) is better than another (vector y) if it has equal and at least one better value on all objectives and can be denoted as x≺y. Suppose there are two objective functions f1 and f2 in a minimization optimization problem, the dominance relation could be presented from an example in Figure 9. This figure shows that the circles are better than squares because they provide a lower value in both objectives.Figure 9
The Pareto dominance relation.For the multiobjective optimization problem, there is a set of nondominated solutions. Consequently, the projection of Pareto optimal solutions in the objective space are stored in a set called Pareto optimal front. The optimization of multiobjective problem aims to find the Pareto optimal front of the problem.


### 3.4.2. Pareto Separation Operator

The Pareto dominance relation is applied in our multiobjective hybrid method to find the optimal operation scheme of the integrated model. Suppose there arem objectives and n acquired solutions in the multiobjective optimization problem, the number of comparison is Cn2⋅m in general, which would increase rapidly with the addition of the objectives and solutions. In order to reduce the number of comparison in our multiobjective algorithm and find the Pareto optimal front quickly, a separation solution comparison method which is called the Pareto separation operator is applied. The process of the PSO is presented as below.It is easy to find that if one solution is nondominated in the multiobjective problem, it would be nondominated in at least one subproblem of the multiobjective optimization problem. That is, if one solution is nondominated in a two-objective subproblem of the multiobjective model, it would be a Pareto optimal solution. The relation between the multiobjective problem and the subproblem is illustrated in Figure10, in which the point in the Pareto front is also the nondominated point in the projection. Therefore, in our comparison method, we only need to find the nondominated solutions in the two-objective subproblem of the multiobjective model, and the union of these solutions is the Pareto optimal solution set.Figure 10
The relation between multiobjective problem and subproblem.The comparison between two objectives is easier than the comparison between more objectives. Suppose there aren solutions in the two-objective subproblem, after we find the optimal solution Sm on the objective 1, we just need to find the solution better than the solution Sm on the objective 2, and these solutions are the Pareto optimal solution of the two-objective subproblem. The number of comparison in our method would be no more than Cm2⋅2n−1, where Cm2 denotes the number of the two-objective subproblem. It is easy to prove that, if n is larger than 2 (m −1), Cn2⋅m would be larger than Cm2⋅2n−1. In the most multiobjective problem, the number of solutions n is far larger than the number of objectives m. Hence, our method would effectively reduce the computational complexity. An example of the comparison method is shown in Figure 11.Figure 11
The example of the Pareto separation operator.


## 3.4.1. Pareto Dominance Relation

Multiobjective optimization deals with optimizing multiple objectives. As mentioned in Section2, the integrated model formulates a problem with three objectives. Different from the optimization with one objective, a solution cannot be compared with another using relational operators in a multiobjective search space because objectives can be in conflict and require special considerations. This is due to the existence of more than one criterion for comparison [26]. Therefore, we need other operators to measure and find out how much a solution is better than another. The Pareto dominance is the most widely used operator to measure the solution of multiobjective model, and a series of multiobjective optimization algorithms are developed based on it [27]. The operator can be mathematically defined in a minimization optimization problem as follows:(12)∀i∈1,2,…,k:fix⟶≤fiy⟶∧∃i∈1,2,…,k:fix⟶≤fiy⟶,where x⟶=x1,x2,…,xk and y⟶=y1,y2,…,yk. This equation shows that a solution (vector x) is better than another (vector y) if it has equal and at least one better value on all objectives and can be denoted as x≺y. Suppose there are two objective functions f1 and f2 in a minimization optimization problem, the dominance relation could be presented from an example in Figure 9. This figure shows that the circles are better than squares because they provide a lower value in both objectives.Figure 9
The Pareto dominance relation.For the multiobjective optimization problem, there is a set of nondominated solutions. Consequently, the projection of Pareto optimal solutions in the objective space are stored in a set called Pareto optimal front. The optimization of multiobjective problem aims to find the Pareto optimal front of the problem.


## 3.4.2. Pareto Separation Operator

The Pareto dominance relation is applied in our multiobjective hybrid method to find the optimal operation scheme of the integrated model. Suppose there arem objectives and n acquired solutions in the multiobjective optimization problem, the number of comparison is Cn2⋅m in general, which would increase rapidly with the addition of the objectives and solutions. In order to reduce the number of comparison in our multiobjective algorithm and find the Pareto optimal front quickly, a separation solution comparison method which is called the Pareto separation operator is applied. The process of the PSO is presented as below.It is easy to find that if one solution is nondominated in the multiobjective problem, it would be nondominated in at least one subproblem of the multiobjective optimization problem. That is, if one solution is nondominated in a two-objective subproblem of the multiobjective model, it would be a Pareto optimal solution. The relation between the multiobjective problem and the subproblem is illustrated in Figure10, in which the point in the Pareto front is also the nondominated point in the projection. Therefore, in our comparison method, we only need to find the nondominated solutions in the two-objective subproblem of the multiobjective model, and the union of these solutions is the Pareto optimal solution set.Figure 10
The relation between multiobjective problem and subproblem.The comparison between two objectives is easier than the comparison between more objectives. Suppose there aren solutions in the two-objective subproblem, after we find the optimal solution Sm on the objective 1, we just need to find the solution better than the solution Sm on the objective 2, and these solutions are the Pareto optimal solution of the two-objective subproblem. The number of comparison in our method would be no more than Cm2⋅2n−1, where Cm2 denotes the number of the two-objective subproblem. It is easy to prove that, if n is larger than 2 (m −1), Cn2⋅m would be larger than Cm2⋅2n−1. In the most multiobjective problem, the number of solutions n is far larger than the number of objectives m. Hence, our method would effectively reduce the computational complexity. An example of the comparison method is shown in Figure 11.Figure 11
The example of the Pareto separation operator.


## 4. Experiments and Results

To verify the effectiveness of the multiobjective optimization model, the experiments are carried out based on a case study on Route 3 of Jingmen Bus (JMB3) in China, which operates 128 trips and carries around 800 passengers per day. 23 buses are currently deployed on this route. The scheduled headway is about 5 to 15 minutes at most time throughout the day, with a maximum of 20 minutes in the afternoon off-peak period.The experimental data preprocessing is presented in Section4.1. This section proposes the approach based on history data to get the passenger flow distribution function, expectation of total passenger flow, and characteristic period division. The parameter setting of the multiobjective heuristic algorithm proposed in Section 3 is presented in Section 4.2. After that, the experimental results are presented in Section 4.3.


### 4.1. Passenger Flow Data and Trip Time Data Processing

#### 4.1.1. Distribution Curve Fitting

The passenger flow and trip time would be changed over time. According to clustering method proposed by Fisher [28], the bus operating period could be divided to f periods with each period called a homogeneous running time (HRT) period [29]. In the same period, the passenger flow has the same distribution function, and the trips with departure time in same period have common time probability distribution. In this paper, the distribution is fitted with normal distribution. An example of the distribution curve fitting is shown as Figure 12, where the blue bar is the passenger flow and the red line is the fitting curve.Figure 12
The normal distribution curve fitting of passenger flow distribution.


#### 4.1.2. HRT Period Division

The HRT period is calculated on loss function. According to the method proposed by Fisher [28], the loss function values of different classification numbers are shown in Figure 13, in which the loss function values corresponding to the classification numbers of 5 and 6 are quite different, and there is not much difference from the loss function with the number of classifications being 6. Therefore, it is appropriate to divide the whole day into six characteristic periods in order to balance the efficiency and complexity of the method.Figure 13
Loss function value of different classification numbers.Based on the method proposed before, the HRT period division and trip time during each period are shown as Table6.Table 6
Trip time during each period on Route 3 of Jingmen Bus.
Characteristic periodsRange of departure interval (min)Expected value of trip time (min)6 : 30–7 : 0010–20627 : 01–8 : 005–15618 : 01–10 : 005–156610 : 01–13 : 0010–207513 : 01–15 : 005–158115 : 01–17 : 305–1578


### 4.2. Parameter Setting of the Multiobjective Hybrid Algorithm

In the study case experiments, the multiobjective algorithm executes 1000 iterations, that is,g = 1000 in the algorithm. Meanwhile, the method L (S) executes 1000 iterations. For the destroy operator 1 removes more trips and searches a larger neighborhood than other two destroy operators, the weight of destroy operator 1 is larger than others. The weights of three destroy operators in L (S) are (0.5, 0.3, 0.2), and the random remove operator removes 10 trips randomly at each iteration L. In the trip insertion, the weights of the vehicle scheduling costs and the penalty costs are assigned as 0.6 and 0.4 separately.


### 4.3. Experimental Results

Using the above data of the study case and the algorithm parameters, the MATLAB version R2017b is used to solve the integrated model. It runs on a personal computer (PC) with Windows 10, Intel Core i5, and 16 GB random access memory (RAM).Figure14 presents a solution of the integrated problem acquired through the multiobjective hybrid algorithm. Figure 15 illustrates the Pareto optimal solution set and initial solution set of the integrated model obtained by the multiobjective method. From Figure 15, we can see that our multiobjective hybrid algorithm is effective to find the Pareto optimal front. Our multiobjective method focuses on finding the solutions with lower vehicle scheduling cost and incompatibility penalty.Figure 14
An example of the integrated problem solution.Figure 15
Pareto optimal solution set and initial solution set of the integrated model.Figure16 illustrates the comparison between the Pareto optimal front of the integrated model, the cost of recursive approach schedule, and the cost of the current implementation schedule on JMB3. The figure shows that the scheduling scheme of the integrated model is partially better than the recursive approach schedule and better than the current implementation schedule.Figure 16
Solution costs of different schedules.The fleet size and objective function value comparison among the Pareto optimal front, the recursive approach schedule, and the current schedule are shown in Table7. Based on Euclidean distance, the closest point of the Pareto optimal solution to the recursive approach schedule and the current implementation schedule are shown in Figure 17.Table 7
Results of the different schedules.
Fleet sizePassenger travel costIncompatible trip-link costVehicle scheduling costPareto optimal solution20–231.2 × 106∼1.4 × 1060–1.2 × 1063.6 × 106∼7.8 × 106Recursive approach schedule231.3 × 1061.5 × 1066.3 × 106Current implementation schedule231.3 × 1061.2 × 1068.7 × 106Figure 17
Closest Pareto optimal point to the different schedules.


## 4.1. Passenger Flow Data and Trip Time Data Processing

### 4.1.1. Distribution Curve Fitting

The passenger flow and trip time would be changed over time. According to clustering method proposed by Fisher [28], the bus operating period could be divided to f periods with each period called a homogeneous running time (HRT) period [29]. In the same period, the passenger flow has the same distribution function, and the trips with departure time in same period have common time probability distribution. In this paper, the distribution is fitted with normal distribution. An example of the distribution curve fitting is shown as Figure 12, where the blue bar is the passenger flow and the red line is the fitting curve.Figure 12
The normal distribution curve fitting of passenger flow distribution.


### 4.1.2. HRT Period Division

The HRT period is calculated on loss function. According to the method proposed by Fisher [28], the loss function values of different classification numbers are shown in Figure 13, in which the loss function values corresponding to the classification numbers of 5 and 6 are quite different, and there is not much difference from the loss function with the number of classifications being 6. Therefore, it is appropriate to divide the whole day into six characteristic periods in order to balance the efficiency and complexity of the method.Figure 13
Loss function value of different classification numbers.Based on the method proposed before, the HRT period division and trip time during each period are shown as Table6.Table 6
Trip time during each period on Route 3 of Jingmen Bus.
Characteristic periodsRange of departure interval (min)Expected value of trip time (min)6 : 30–7 : 0010–20627 : 01–8 : 005–15618 : 01–10 : 005–156610 : 01–13 : 0010–207513 : 01–15 : 005–158115 : 01–17 : 305–1578


## 4.1.1. Distribution Curve Fitting

The passenger flow and trip time would be changed over time. According to clustering method proposed by Fisher [28], the bus operating period could be divided to f periods with each period called a homogeneous running time (HRT) period [29]. In the same period, the passenger flow has the same distribution function, and the trips with departure time in same period have common time probability distribution. In this paper, the distribution is fitted with normal distribution. An example of the distribution curve fitting is shown as Figure 12, where the blue bar is the passenger flow and the red line is the fitting curve.Figure 12
The normal distribution curve fitting of passenger flow distribution.


## 4.1.2. HRT Period Division

The HRT period is calculated on loss function. According to the method proposed by Fisher [28], the loss function values of different classification numbers are shown in Figure 13, in which the loss function values corresponding to the classification numbers of 5 and 6 are quite different, and there is not much difference from the loss function with the number of classifications being 6. Therefore, it is appropriate to divide the whole day into six characteristic periods in order to balance the efficiency and complexity of the method.Figure 13
Loss function value of different classification numbers.Based on the method proposed before, the HRT period division and trip time during each period are shown as Table6.Table 6
Trip time during each period on Route 3 of Jingmen Bus.
Characteristic periodsRange of departure interval (min)Expected value of trip time (min)6 : 30–7 : 0010–20627 : 01–8 : 005–15618 : 01–10 : 005–156610 : 01–13 : 0010–207513 : 01–15 : 005–158115 : 01–17 : 305–1578


## 4.2. Parameter Setting of the Multiobjective Hybrid Algorithm

In the study case experiments, the multiobjective algorithm executes 1000 iterations, that is,g = 1000 in the algorithm. Meanwhile, the method L (S) executes 1000 iterations. For the destroy operator 1 removes more trips and searches a larger neighborhood than other two destroy operators, the weight of destroy operator 1 is larger than others. The weights of three destroy operators in L (S) are (0.5, 0.3, 0.2), and the random remove operator removes 10 trips randomly at each iteration L. In the trip insertion, the weights of the vehicle scheduling costs and the penalty costs are assigned as 0.6 and 0.4 separately.


## 4.3. Experimental Results

Using the above data of the study case and the algorithm parameters, the MATLAB version R2017b is used to solve the integrated model. It runs on a personal computer (PC) with Windows 10, Intel Core i5, and 16 GB random access memory (RAM).Figure14 presents a solution of the integrated problem acquired through the multiobjective hybrid algorithm. Figure 15 illustrates the Pareto optimal solution set and initial solution set of the integrated model obtained by the multiobjective method. From Figure 15, we can see that our multiobjective hybrid algorithm is effective to find the Pareto optimal front. Our multiobjective method focuses on finding the solutions with lower vehicle scheduling cost and incompatibility penalty.Figure 14
An example of the integrated problem solution.Figure 15
Pareto optimal solution set and initial solution set of the integrated model.Figure16 illustrates the comparison between the Pareto optimal front of the integrated model, the cost of recursive approach schedule, and the cost of the current implementation schedule on JMB3. The figure shows that the scheduling scheme of the integrated model is partially better than the recursive approach schedule and better than the current implementation schedule.Figure 16
Solution costs of different schedules.The fleet size and objective function value comparison among the Pareto optimal front, the recursive approach schedule, and the current schedule are shown in Table7. Based on Euclidean distance, the closest point of the Pareto optimal solution to the recursive approach schedule and the current implementation schedule are shown in Figure 17.Table 7
Results of the different schedules.
Fleet sizePassenger travel costIncompatible trip-link costVehicle scheduling costPareto optimal solution20–231.2 × 106∼1.4 × 1060–1.2 × 1063.6 × 106∼7.8 × 106Recursive approach schedule231.3 × 1061.5 × 1066.3 × 106Current implementation schedule231.3 × 1061.2 × 1068.7 × 106Figure 17
Closest Pareto optimal point to the different schedules.


## 5. Conclusion

This paper has proposed a multiobjective optimization approach for timetabling and vehicle scheduling with stochastic trip time and passenger flow. A multiobjective optimization model has been established comprehensively with the objectives of minimizing the passenger travel costs, vehicle scheduling costs, and incompatible trip-link costs while satisfying the passenger and bus corporation demands and improving the scheduling robustness by considering the uncertainty. In this paper, the objective of minimizing the passenger travel costs is different from the previous research projects, which optimize the timetabling based on headway smoothing rule by considering the passenger waiting cost and comfort cost. And the objective of minimizing the incompatible trip-link costs is to penalize the incompatible trip-link chain attributed to the stochastic trip time and enhance the robustness of the vehicle schedule. Meanwhile, the probability distribution of the stochastic trip time and passenger flow following are obtained by the measurement history data in advance. After the model establishment, a multiobjective hybrid algorithm, which is combined with self-adjust GA, LNS, and Pareto separation operator, is proposed to solve the problem. A case study on real instance of Route 3 of Jingmen Bus (JMB3) is reported.The experimental results on the problem instance derived from the JMB3 case show that our multiobjective optimization approach can produce the scheduling schemes that perform better than the current implementation schedule at one or more objectives. Compared with the recursive approach, the multiobjective hybrid method proposed in this paper can generate the scheduling scheme that balances the different costs.Since the scheduling schemes generated by the multiobjective optimization method are a series of feasible schemes, our proposed approach in this paper would contribute to the transit network planning. Although the method in this paper performs well, there is still room for further improvement. For example, the model does not consider the stochastic deadhead time, and the probability distribution model of passenger flow and trip time is derived from limited data; hence, the probability distribution model is simple. These problems will require more work in the future.



---

*Source: 3529984-2021-08-05.xml*