# Optimizing Vehicle Scheduling Based on Variable Timetable by Benders-and-Price Approach

Authors: Zekang Lan; Shiwei He; Rui Song; Sijia Hao
Journal: Journal of Advanced Transportation
 (2019)

---

## Abstract

In practice, vehicle scheduling is planned on a variable timetable so that the departure times of trips can be shifted in tolerable ranges, rather than on a fixed timetable, to decrease the required fleet size. This paper investigates the vehicle scheduling problem on a variable timetable with the constraint that each vehicle can perform limited trips. Since the connection-based model is difficult to solve by optimization software for a medium-scale or large-scale instance, a designed path-based model is developed. A Benders-and-Price algorithm by combining the Benders decomposition and column generation is proposed to solve the LP-relaxation of the path-based model, and a bespoke Branch-and-Price is used to obtain the integer solution. Numerical experiments indicate that a variable timetable approach can reduce the required fleet size with a tolerable timetable deviation in comparison with a fixed timetable approach. Moreover, the proposed algorithm is greatly superior to GUROBI in terms of computational efficiency and guarantees the quality of the solution.


---

## Body

## 1. Introduction

In the planning process of public bus transport, timetabling problem, and vehicle scheduling problem (VSP) are dealt with by a sequential approach where the solution of previous subproblem is taken as input of the following subproblem, because of the complexity of solving the integrated model of the two planning processes [1]. Both the timetabling problem and VSP have been well-studied. The former one determines the departure and arrival times at each station for each trip, aiming to offer a high level of service for passengers [2–4]. Based on the timetabled trips, the latter one focuses on minimizing the operational cost by arranging the vehicles to cover each trip exactly once and turning back to their corresponding origin terminals or depots after finishing a defined work period, e.g., one day. The connections of trips depend on the arrival/departure times at the terminals. Hence, the process of vehicle scheduling is highly dependent on the timetable. Thus, in order to lower the number of costly vehicles used in vehicle scheduling, a variable timetable with trips shifting in acceptable tolerances is adopted instead of a fixed timetable. However, the process of shifting is a very time-consuming mission and usually not done in a systematic manner [2]. Ceder [2] and Liu and Ceder [5] have presented deficit-function procedures with and without insertions of deadheading trips to assist vehicle schedulers in achieving the most efficient shifts. However, when taking into account fuel restrictions, maintenance considerations, etc., each vehicle has a limited daily workload that implies that the number of trips taken by one vehicle in one day cannot exceed an upper bound (limited trips, LT). In this case, the deficit-function method is not applicable and a new solution method should be proposed. To the best of our knowledge, the VSP based on a variable timetable with limited trips (VSPVT-LT) has not been studied.Previous studies typically researched into vehicle scheduling based on a fixed timetable. The usual modeling framework of the VSP is built on a connection-based network, where trips and depots are expressed as vertices, and compatible connections between any two vertices are connected using arcs [6–8]. Considering that the connection-based method could bring an immense number of deadhead connection arcs when the number of trips grows, a space-time network modeling concept has been commonly used [8, 9]. Bertossi et al. [10] have proved that single depot vehicle scheduling problem (SDVSP) can be solved in polynomial time; however multiple depot vehicle scheduling problem (MDVSP) is NP-hard even for two depots. Generally, the connection-based model is hard to solve by optimization software when thousands of trips have to be scheduled. Costa et al. [11] and Hadjar et al. [12] formulated the VSP as a Set Portioning Problem which can be solved successfully by column generation. Bodin et al. [13] and Freling and Paixão [14] discussed the vehicle scheduling problem with time constraint, and they constructed a special connection-based network with inserting additional arcs or alternative graphs with different levels to handle the time constraint. Desaulniers et al. [15] and Hadjar and Soumis [16] researched into the vehicle scheduling problem with time window (VSPTW). However, considering that the time windows of departure times may overlap, the minimum headway requirement, which is usually not considered in the VSPTW, has to be imposed in the VSPVT-LT.In recent years, integrated optimization models for two consecutive phases of public transit have attracted more and more attention. Ibarra-Rojas et al. [17] proposed a completed biobjective integration model that combines the SDVSP and the timetabling problem, and a ε-constraint method was implemented to obtain Pareto-optimal solutions. On the other hand, most of the literature addresses a sequential integration of the timetabling problem with the VSP. A sequential integration is an approach that considers the characteristics of one subproblem while another subproblem is being optimized. Van den Heuvel et al. [18] applied an iterated heuristic method between a local search strategy to alter the timetable slightly and a network-flow model for vehicle scheduling. In each iteration, the modification of the timetable with an improvement is always adopted, while with a degradation it is accepted with a probability which decreases in the wake of the algorithm processing. Guihaire and Hao [19] also presented an iterated local search algorithm where each iteration implements trip shifting to obtain a neighboring timetable and then the VSP is optimized by an efficient auction algorithm on the neighboring timetable. Liu et al. [20] provided a new biobjective, bilevel mathematical programming model, and a novel deficit-function-based sequential search approach by combining a network-flow technique and a shifting departure time procedure, presented to solve the problem to achieve a set of Pareto-efficient solutions. Fonseca et al. [21] proposed a metaheuristic that iteratively optimized the mathematical formulation of the integrated timetabling and VSP that permitted a subset of timetabled trips shifting only, whereas solving the full VSP.The remainder of this paper is organized as follows. In Section2, the overall problem description and notations are elaborated, and a directed network is built to make provision for the connection-based model in the Appendix and the shortest path algorithm in Section 4.2. The path-based model is put forward in Section 3. A computationally efficient solution method that consists of two phases to solve the path-based model is presented in Section 4. The efficiency of the proposed algorithm and the benefits of the vehicle scheduling in a variable timetable are systematically evaluated with numerical experiments in Section 5. Finally, conclusion and future work extensions are discussed in Section 6.


## 2. The Problem Description and Notations

This study considers a bus line or a set of bus lines. At the ends of each line are terminals. An initial feasible daily timetable generated by timetabling specifies the origin terminal, destination terminal, departure time at origin terminal, and arrival time at destination terminal of each trip. A timetable is feasible if several constraints such as the minimum headway requirement and the minimum travel time of trip are satisfied, where the headway time is the time interval between the departure times of two adjacent trips with the same origin terminal. Considering that the departure time of each trip can be shifted backward or forward in a given range and respecting the minimum headway requirement of operation, there are a lot of possible feasible timetables. In our case, fuel restrictions and maintenance plans are taken into account; that is, the running mileage or operational time of each vehicle in one day cannot exceed a defined upper bound for the vehicle. Assuming that there are a few differences between the running mileages or travel times of any two trips, the restriction is equivalent to limiting the number of trips carried out by each vehicle. The VSPVT-LT has to find out a feasible timetable with the least fleet size in vehicle scheduling with LT constraints, but also needs to balance the timetable deviation in relation to the initial timetable. We do not consider the inclusion of deadheading trips insertions, and we assume that the travel time of each trip is fixed in timetable modification. For illustrative purpose, some notations for the entire article are firstly introduced in Table1.Table 1
Some notations for the entire article.
Sets T Set of trips,T=i:i=1,⋯,n,where n is the total number of trips Ω Set of all the feasible timetables K Set of terminals R Set of vehicles,R=r:r=1,2,⋯,u, where u is the total number of vehicles R k Set of vehicles located in terminalk∈K, Rk⊂R T d k T a k Sets of trips that start from (end at) the terminalk∈K Parameters t 0 The given initial daily timetable D i A i The departure (arrival) time of tripi∈T in t0 l i The travel time of tripi∈T λ i The penalty cost for shifting the departure time of tripi∈T per unit time, e.g., one minute w i The penalty cost for canceling tripi∈T c r The cost for vehicler∈R if it is used c R c R = min ⁡ c r ∣ r ∈ R o i d i The origin (destination) terminal of tripi∈T h The minimum headway of operational requirement Δ d i - Δ d i + The maximum allowed earliness (tardiness) for the departure time of tripi∈T χ m i n χ m a x The minimum (maximum) time for a turnaround operation N r The maximum allowed number of trips executed by vehicler∈R Decision variables x i j k Binary variables, 1 if vehicler∈R travels from vertex (trip) i directly to vertex (trip) j, and 0 otherwise d i a i Non-negative integer variables, indicating the departure (arrival) time of tripi∈T in the final output feasible timetable m i Binary variables, 1 if tripi∈T is canceled, and 0 otherwiseThe maximum turnaround timeχmax is used to reduce nonproductive time and the number of connection arcs in the connection-based model to improve computation efficiency. Two trips i,j∈T are compatible in a given timetable t∈Ω if they can be performed by the same vehicle in sequence. The route is a sequence of trips, and any two consecutive trips in a route are compatible. For any two trips i,j∈T in timetable t0 such that (i) d(i)=o(j), (ii) (Dj+Δdj+)-(Ai-Δdi-)≥χmin, and (iii) (Dj-Δdj-)-(Ai+Δdi+)≤χmax, we call them restrictively compatible, since i and j will be compatible by adjusting their departure times in a given range even if i and j are incompatible in t0. If i and j are restrictively compatible, this is denoted as i↔j.To facilitate explanation for VSPVT-LT, consider a bus line with two terminals, which has 14 trips in its initial timetable as depicted in Figure1. Assuming that the minimum turnaround time is 10 min, the minimum headway time is 3 min, and the maximum turnaround time and vehicle daily workload are not limited, we seek to obtain the vehicle schedule solution based on the fixed timetable with the least required fleet size, which has four routes and they are revealed by the different color trajectories. For instance, the route as the red trajectory shows starts from depot 1, travels trip 1, trip 10, trip 5, and trip 14, sequentially, and end at depot 1. Additionally, because a feasible route should always have the same starting and ending depots, trip 2 and trip 9 cannot put into service, so they have to be canceled in the final timetable. The result will certainly be unsatisfactory to the planning staffs.Figure 1
The initial timetable and the corresponding vehicle schedule solution.Considering the case on a variable timetable, e.g., we setΔdi-=Δdi+ = 10 min for each trip i∈T, a modified timetable shown in Figure 2 is obtained from the initial timetable by shifting the departure times of trip 8 and trip 9 backward 10 min, trip 2 and trip 3 forward 5 min, respectively. Similarly, the corresponding vehicle schedule solution is generated and reported in Figure 2. Despite the same number of vehicles used as before adjustment, all trips are performed by exactly one vehicle. In the example, at first sight, it is a relatively simple structure. However, one should not draw a conclusion that these connections can be accomplished in manual planning, because the example exhibits an isolated part of a complicate structure, especially when LT constraints are taken into account.Figure 2
The modified timetable and the corresponding vehicle schedule solution.We associate each terminalk∈K with a directed network Gk=(Vk,Ak), where Vk and Ak denote the sets of vertices and arcs in this network, respectively. Each trip i∈T is represented by a trip-vertex, and each terminal k∈K by a source-vertex n+2k-1 and sink-vertex n+2k (every terminal is duplicated into a source-vertex and a sink-vertex). The vertex set Vk contains all the trip-vertices, corresponding source-vertex n+2k-1, and sink-vertex n+2k. For simplicity, the set of trip-vertices is also represented by T. Thus Vk=T∪n+2k-1,n+2k. The set Ak contains the start arcs n+2k-1×Td(k), the end arcs Ta(k)×n+2k, and the inner arcs I=(i,j)∣i,j∈T:i↔j. Thus, Ak is defined as Ak = (n+2k-1×Td(k)) ∪(Ta(k)×n+2k)∪I. For any vertex i∈Vk in each network Gk(k∈K), Ik(i) and Ok(i) specify the set of ingoing and outgoing vertices of vertex i, respectively.


## 3. The Path-Based Model for VSPVT-LT

In order to evaluate the benefits of the path-based model, we also develop a connection-based model, which is detailed in the Appendix. We will verify that the connection-based model is hard to solve by GUROBI in the medium-scale or large-scale instances in Section5.2. To overcome the handicap, the column generation technology should be employed. Embedded in Branch-and-bound (B&B) algorithm, column generation has been applied to successfully solve real-world instances of routing and scheduling problems [16].Let us firstly define the term “path” that is distinct from the “route” defined before. The only difference is that any two consecutive trips in a path are restrictively compatible in the initial timetablet0, without requiring being compatible. However, a path will be a route in a certain timetable t∈Ω. Considering the example mentioned in Section 2 again, the sequence of trips (8,2,11,6) is a path in the initial timetable rather than a route, because trip 8 and trip 2 are incompatible but restrictively compatible. However, (8,2,11,6) is also a route in the modified timetable. The notations for the path-based model are presented in Table 2.Table 2
Notations for the path-based model.
Sets S Set of vehicle types. Different types vehicles can execute different maximum allowed number of trips S k Set of vehicle types in terminalk, Sk⊂S ϑ , ϑ k ϑ = max ⁡ s ∣ s ∈ S, ϑk=max⁡s∣s∈Sk P Set of all possible paths P k , s Set of paths that can be covered by a vehicle from terminalk∈K of type s∈Sk Parameters c s The cost of pathp∈P covered by a vehicle of type s∈S. We suppose that two different vehicles with the same type have the same cost ς k , s The number of vehicles of types∈Sk in terminal k∈K, ∑k∈K∑s∈Skςk,s=u δ p i A binary constant taking the value 1 if the pathp includes the trip i, 0, otherwise θ p i j A binary constant taking the value 1 if the pathp includes trips i and j in sequence, 0, otherwise Decision variables x s p A binary variable which takes value 1 if the pathp is selected and covered by a vehicle of type s∈S in an optimal solution, 0, otherwise y i j A binary variable which takes value 1 if tripsi and j are compatible, 0, otherwiseThe path-based model is now given as follows.(1)min∑i∈Tλidi-Di+∑k∈K∑s∈Sk∑p∈Pk,scsxsp+∑i∈Twi1-∑k∈K∑s∈Sk∑p∈Pk,sδpixsp=∑i∈Tλidi-Di+∑i∈Twi+∑k∈K∑s∈Sk∑p∈Pk,scs-∑i∈Twiδpixsp(2)s.t.constraintsA.2-A.4See  the  Appendix(3)dj-ai≥χmin-1-yijM∀i,j∈I(4)dj-ai≤χmax+1-yijM∀i,j∈I(5)∑p∈Pk,sxsp≤ςk,s∀k∈K,s∈Sk(6)∑k∈K∑s∈Sk∑p∈Pk,sδpixsp≤1∀i∈T(7)∑k∈K∑s∈Sk∑p∈Pk,sθpijxsp≤yij∀i,j∈I(8)xsp∈0,1∀k∈K,s∈Sk,p∈Pk,sThe objective function (1) is to minimize the total cost as (A.1) (see the Appendix) expresses. Constraints (3) and (4) ensure that if trips i and j are compatible, then the difference between the departure time of j and the arrival time of i has to be greater than or equal to χmin, and smaller than or equal to χmax. Constraints (5) make sure that there are adequate vehicles to execute the paths being selected. Constraints (6) require that for each trip i∈T, at most one path is selected among those paths containing trip i. Constraints (7) state that if trips i and j are incompatible, the paths containing i and j in sequence are forbidden to be selected. Constraints (8) define the domains of the variables.


## 4. Solution Method

### 4.1. Benders-and-Price Solution Framework

In the path-based model, we observe that the variables associated with the timetabling problem, i.e.,ai and di, and the variables associated with the VSP, i.e., xs(p), are linked together by yij. For a given timetable and corresponding values of yij satisfying constraints (3) and (4), the LP-relaxation of path-based model only involves the variables xs(p) and can be solved easily by column generation. Therefore, it is reasonable to use Benders decomposition technology to divide the original problem into two subproblems as follows. The Benders master problem (BMP) is the timetabling problem that determines the timetable and corresponding variables of yij. The Benders subproblem (BSP) is a path-based MDVSP model which can be solved by column generation. The approach that combines Benders decomposition and column generation has been applied to other simultaneous problems in transportation, such as locomotives and cars assignment to passenger trains problem [22], and aircraft routing and crew scheduling problem [23]. For more recent literature about this approach, we refer the reader to Restrepo et al. [24].Introducing an additional variablez0(z0≥0), the Benders master problem can be formulated as(9)BMPmin∑i∈Tλidi-Di+z0(10)s.t.constraintsA.2-A.4,3-4,and(11)z0≥z0downwhere z0down is the lower bound of z0. Because relative large penalty costs are set for canceling trips as mentioned in the Appendix, it is more cost-efficient to use a vehicle than to cancel a trip. Therefore, the lower bound of the objective function (1) can be the least cost of using vehicles, which can be expressed as (n/ϑ)∗cR. Moreover, remember that the objective value of BMP is the lower bound of the objective function (1). By comparing the objective function (1) and (9), it is reasonable to set z0down as (n/ϑ)∗cR. After solving BMP, the departure times of trips and the variables of yij((i,j)∈I) are fixed to particular values d-i(i∈T) and y-ij((i,j)∈I). The path-based model reduces to the following primal Benders subproblem (BSP).(12)BSP∑i∈Tλid-i-Di+∑i∈Twi+min⁡∑k∈K∑s∈Sk∑p∈Pk,scs-∑i∈Twiδpixsp(13)∑p∈Pk,sxsp≤ςk,s∀k∈K,s∈Sk(14)∑k∈K∑s∈Sk∑p∈Pk,sδpixsp≤1∀i∈T(15)∑k∈K∑s∈Sk∑p∈Pk,sθpijxsp≤y-ij∀i,j∈I(16)xsp∈0,1∀k∈K,s∈Sk,p∈Pk,sObviously, it is impractical to enumerate the huge number of paths in the set ofP and corresponding P(k,s)(k∈K,s∈Sk). Therefore, a column generation method should be applied to generate the more profitable paths. Considering a LP-relaxation of restricted BSP (denoted by RBSP) obtained from the BSP by replacing P(k,s) with its subset P-(k,s) and replacing (16) with xs(p)∈[0,1]. Let π=(πk,s≤0), μ=(μi≤0), and σ=(σij≤0) be the dual variables associated with constraints (13), (14), and (15), respectively. For each k∈K and s∈Sk, the reduced cost of path p for the dual solution (π, μ, σ) is given by cs-∑i∈Twiδpi-∑k∈K∑s∈Skπk,s-∑i∈Tδpiμi-∑(i,j)∈Iθpijσij. Thus, identifying a column with minimum reduced cost of the pricing step for vehicle type s∈Sk(k∈K) is equivalent to the pricing subproblem.(17)PSPζk,s=minp∈P-k,s⁡cs-∑i∈Twiδpi-∑k∈K∑s∈Skπk,s-∑i∈Tδpiμi-∑i,j∈IθpijσijIf the objective value of PSPζk,s<0, we can enter the corresponding path into the path set P-(k,s) to improve the solution. If for each k∈K and s∈Sk, ζk,s≥0, the optimal solution of RBSP at the current iteration has been obtained. In order to solve the pricing subproblems faster, we present a label correcting shortest path algorithm, which will be talked about in the next section.The RBSP is always feasible since the null vector0 satisfies constraints (13)-(16), then artificial variables and Benders feasibility cuts are not essential to be generated in RBSP. The objective values of RBSP and BMP are the upper bounds and lower bounds at Benders-and-Price approach iterations, denoted as ηUB and ηLB, respectively. Note that the last convergence value of the Benders decomposition is actually the lower bound of the optimal objective value of the original path-based model, because we have applied LP-relaxation to RBSP. If ηUB-ηLB/ηUB≤ε, where ε is a predefined tolerance, Benders-and-Price can be terminated. Otherwise, a Benders cut as presented below has to be imposed into BMP.(18)z0≥∑i∈Twi+∑k∈K∑s∈Skςk,sπk,s+∑i∈Tμi+∑i,j∈Iσijyij


### 4.2. Shortest Path Algorithm for Pricing Problem

In order to solve the problem PSP faster, a shortest path algorithm with vertices constrained, i.e., the shortest path consisting of a limited number of vertices, is proposed. In Section2, we have associated each terminal k∈K with a directed network Gk=(Vk,Ak). Each arc (i,j)∈Ak has a cost cost(i,j), which can be determined by function (19) below. More specifically, for each start arc (i,j), the cost is mapped as -wj-μj. For each inter-arc (i,j)∈I, the cost includes not only -wj-μj, but also the connection cost -σij between vertices i and j. For each end arc, the cost is 0. The reason that the negative sign “-” remains before the dual values μ and σ is to comply with (17) and to keep the problem as a minimization problem.(19)costi,j=-wj-μji,j∈n+2k-1×Tdk0i,j∈Tak×n+2k-σij-wj-μji,j∈IThe shortest path algorithm is to find a least cost path from the source-vertexn+2k-1 to the sink-vertex n+2k for each k∈K in the corresponding graph Gk. Since the cost of arcs may be negative numbers, the Bellman-Ford algorithm is adopted to get the shortest path. Remember that there is an upper bound of intermediate vertices requirement on the shortest path, and although the standard version of Bellman-Ford algorithm does not comprise a constraint on the number of vertices used, it is not difficult to modify the algorithm by making full use of the methodology in which each iteration of the Bellman-Ford algorithm finds the shortest path containing one vertex more than previous iteration to add such a restriction. Obviously, with different upper bounds of intermediate vertices, the output shortest paths may be different. Before detailing the shortest path, some notations are firstly given. For each k∈K, let w(m,i) be the least cost (i.e., label) from source-vertex n+2k-1 to vertex i with at most m(m≤ϑk+2) vertices and pv(m,i) be the preceding vertex of vertex i on the shortest path that uses at most m(m≤ϑk+2) vertices. Pseudocode for the shortest path algorithm is given in Algorithm 1.Algorithm 1:The shortest path algorithm with vertices constrained in Gk.
Step 1. Initialization. Setwm,j=∞, ∀j∈Vk∖n+2k-1, m=1,2,⋯,ϑk+2; wm,n+2k-1=0, ∀m=1,2,⋯,ϑk+2; p v m , j = n o n e, ∀j∈Vk, m=1,2,⋯,ϑk+2. Step 2. Label updating Form=2 to ϑk+2 For each arc i,j∈Ak If wm,j>wm-1,i+costi,j Then w m , j = w m - 1 , j + c o s t i , j p v m , j = i End if End for End for Step 3. Fetch the shortest paths For eachs∈S, do Set the sink-vertex n+2k as the current vertex i; For m=s+2 to 1 If vertex i is not the source node n+2k-1 Find the preceding vertex pvm,i of the current vertex i, and update pvm,i as the current vertex i; End if End for Reverse the backward path and output the least cost path from n+2k-1 to n+2k with at most s intermedium vertices. End forInStep 3 of Algorithm 1, a shortest path is fetched for each s∈S; thus the number of output paths may be more than one. Note that the same paths may be generated with different values of s∈S. If this happens, what we need to do is just to eliminate the duplicates path(s) in the output result.


### 4.3. Overall Algorithm

The overall algorithm comprises two phases. In the first stage, the Benders-and-Price approach alternates between RBSP and BMP with Benders cuts (18) generation until the gap between ηLB and ηUB is small enough. The BMP can be solved by IP solvers, like GUROBI, and RBSP is solved by column generation method (not embedded in B&B). If the final solution is integer valued in the end of the first phase algorithm, then the second phase can be eliminated. If not, the algorithm then enters the second phase, where it seeks integer solution of RBSP by performing a Branch-and-Price (BAP) algorithm. Note that, after the first phase, the timetable has been fixed. Figure 3 gives a flowchart for the overall algorithm.Figure 3
Flowchart of the overall algorithm.As we can see in the BAP algorithm, the branching rule which branches on the variablesxs(p) with fractional values near 1 and Depth-First node picking strategy with focus on “x=1” cuts are adopted. Actually, an exact branching scheme is to branch decisions on two compatible trips to be covered by the same vehicle [15, 16]. However, for large instances, this scheme is likely to require excessive computing time because of the large number of possible connections [23]. On the other hand, in the way of branching on the path variables, once one path (also a route) is fixed to be selected, all the paths which are conflicted with the fixed path will be eliminated. Moreover, recall that canceling trips is allowed, so infeasibility never occurs during the search. Then, it is hopeful to obtain a satisfactory upper bound. For a CSP imposed by an “x=0” cut from its parent subproblem, we refrain from solving it by column generation, but dual simplex method, because the type of “x=0” cut changes the solution space slightly, and using column generation to solve will result in the same paths being generated again, which hinders the convergence of algorithm. Although the final solution can be a suboptimal solution without solving each subproblem by column generation, the final solution is always high-quality and even optimal in our test cases. Last but not least, if xp∗ is a fractional solution, it is not always necessary to branch on the fractional variable, only when it is likely to save at least one vehicle, i.e., z(xp∗)≤UB-cR.


## 4.1. Benders-and-Price Solution Framework

In the path-based model, we observe that the variables associated with the timetabling problem, i.e.,ai and di, and the variables associated with the VSP, i.e., xs(p), are linked together by yij. For a given timetable and corresponding values of yij satisfying constraints (3) and (4), the LP-relaxation of path-based model only involves the variables xs(p) and can be solved easily by column generation. Therefore, it is reasonable to use Benders decomposition technology to divide the original problem into two subproblems as follows. The Benders master problem (BMP) is the timetabling problem that determines the timetable and corresponding variables of yij. The Benders subproblem (BSP) is a path-based MDVSP model which can be solved by column generation. The approach that combines Benders decomposition and column generation has been applied to other simultaneous problems in transportation, such as locomotives and cars assignment to passenger trains problem [22], and aircraft routing and crew scheduling problem [23]. For more recent literature about this approach, we refer the reader to Restrepo et al. [24].Introducing an additional variablez0(z0≥0), the Benders master problem can be formulated as(9)BMPmin∑i∈Tλidi-Di+z0(10)s.t.constraintsA.2-A.4,3-4,and(11)z0≥z0downwhere z0down is the lower bound of z0. Because relative large penalty costs are set for canceling trips as mentioned in the Appendix, it is more cost-efficient to use a vehicle than to cancel a trip. Therefore, the lower bound of the objective function (1) can be the least cost of using vehicles, which can be expressed as (n/ϑ)∗cR. Moreover, remember that the objective value of BMP is the lower bound of the objective function (1). By comparing the objective function (1) and (9), it is reasonable to set z0down as (n/ϑ)∗cR. After solving BMP, the departure times of trips and the variables of yij((i,j)∈I) are fixed to particular values d-i(i∈T) and y-ij((i,j)∈I). The path-based model reduces to the following primal Benders subproblem (BSP).(12)BSP∑i∈Tλid-i-Di+∑i∈Twi+min⁡∑k∈K∑s∈Sk∑p∈Pk,scs-∑i∈Twiδpixsp(13)∑p∈Pk,sxsp≤ςk,s∀k∈K,s∈Sk(14)∑k∈K∑s∈Sk∑p∈Pk,sδpixsp≤1∀i∈T(15)∑k∈K∑s∈Sk∑p∈Pk,sθpijxsp≤y-ij∀i,j∈I(16)xsp∈0,1∀k∈K,s∈Sk,p∈Pk,sObviously, it is impractical to enumerate the huge number of paths in the set ofP and corresponding P(k,s)(k∈K,s∈Sk). Therefore, a column generation method should be applied to generate the more profitable paths. Considering a LP-relaxation of restricted BSP (denoted by RBSP) obtained from the BSP by replacing P(k,s) with its subset P-(k,s) and replacing (16) with xs(p)∈[0,1]. Let π=(πk,s≤0), μ=(μi≤0), and σ=(σij≤0) be the dual variables associated with constraints (13), (14), and (15), respectively. For each k∈K and s∈Sk, the reduced cost of path p for the dual solution (π, μ, σ) is given by cs-∑i∈Twiδpi-∑k∈K∑s∈Skπk,s-∑i∈Tδpiμi-∑(i,j)∈Iθpijσij. Thus, identifying a column with minimum reduced cost of the pricing step for vehicle type s∈Sk(k∈K) is equivalent to the pricing subproblem.(17)PSPζk,s=minp∈P-k,s⁡cs-∑i∈Twiδpi-∑k∈K∑s∈Skπk,s-∑i∈Tδpiμi-∑i,j∈IθpijσijIf the objective value of PSPζk,s<0, we can enter the corresponding path into the path set P-(k,s) to improve the solution. If for each k∈K and s∈Sk, ζk,s≥0, the optimal solution of RBSP at the current iteration has been obtained. In order to solve the pricing subproblems faster, we present a label correcting shortest path algorithm, which will be talked about in the next section.The RBSP is always feasible since the null vector0 satisfies constraints (13)-(16), then artificial variables and Benders feasibility cuts are not essential to be generated in RBSP. The objective values of RBSP and BMP are the upper bounds and lower bounds at Benders-and-Price approach iterations, denoted as ηUB and ηLB, respectively. Note that the last convergence value of the Benders decomposition is actually the lower bound of the optimal objective value of the original path-based model, because we have applied LP-relaxation to RBSP. If ηUB-ηLB/ηUB≤ε, where ε is a predefined tolerance, Benders-and-Price can be terminated. Otherwise, a Benders cut as presented below has to be imposed into BMP.(18)z0≥∑i∈Twi+∑k∈K∑s∈Skςk,sπk,s+∑i∈Tμi+∑i,j∈Iσijyij


## 4.2. Shortest Path Algorithm for Pricing Problem

In order to solve the problem PSP faster, a shortest path algorithm with vertices constrained, i.e., the shortest path consisting of a limited number of vertices, is proposed. In Section2, we have associated each terminal k∈K with a directed network Gk=(Vk,Ak). Each arc (i,j)∈Ak has a cost cost(i,j), which can be determined by function (19) below. More specifically, for each start arc (i,j), the cost is mapped as -wj-μj. For each inter-arc (i,j)∈I, the cost includes not only -wj-μj, but also the connection cost -σij between vertices i and j. For each end arc, the cost is 0. The reason that the negative sign “-” remains before the dual values μ and σ is to comply with (17) and to keep the problem as a minimization problem.(19)costi,j=-wj-μji,j∈n+2k-1×Tdk0i,j∈Tak×n+2k-σij-wj-μji,j∈IThe shortest path algorithm is to find a least cost path from the source-vertexn+2k-1 to the sink-vertex n+2k for each k∈K in the corresponding graph Gk. Since the cost of arcs may be negative numbers, the Bellman-Ford algorithm is adopted to get the shortest path. Remember that there is an upper bound of intermediate vertices requirement on the shortest path, and although the standard version of Bellman-Ford algorithm does not comprise a constraint on the number of vertices used, it is not difficult to modify the algorithm by making full use of the methodology in which each iteration of the Bellman-Ford algorithm finds the shortest path containing one vertex more than previous iteration to add such a restriction. Obviously, with different upper bounds of intermediate vertices, the output shortest paths may be different. Before detailing the shortest path, some notations are firstly given. For each k∈K, let w(m,i) be the least cost (i.e., label) from source-vertex n+2k-1 to vertex i with at most m(m≤ϑk+2) vertices and pv(m,i) be the preceding vertex of vertex i on the shortest path that uses at most m(m≤ϑk+2) vertices. Pseudocode for the shortest path algorithm is given in Algorithm 1.Algorithm 1:The shortest path algorithm with vertices constrained in Gk.
Step 1. Initialization. Setwm,j=∞, ∀j∈Vk∖n+2k-1, m=1,2,⋯,ϑk+2; wm,n+2k-1=0, ∀m=1,2,⋯,ϑk+2; p v m , j = n o n e, ∀j∈Vk, m=1,2,⋯,ϑk+2. Step 2. Label updating Form=2 to ϑk+2 For each arc i,j∈Ak If wm,j>wm-1,i+costi,j Then w m , j = w m - 1 , j + c o s t i , j p v m , j = i End if End for End for Step 3. Fetch the shortest paths For eachs∈S, do Set the sink-vertex n+2k as the current vertex i; For m=s+2 to 1 If vertex i is not the source node n+2k-1 Find the preceding vertex pvm,i of the current vertex i, and update pvm,i as the current vertex i; End if End for Reverse the backward path and output the least cost path from n+2k-1 to n+2k with at most s intermedium vertices. End forInStep 3 of Algorithm 1, a shortest path is fetched for each s∈S; thus the number of output paths may be more than one. Note that the same paths may be generated with different values of s∈S. If this happens, what we need to do is just to eliminate the duplicates path(s) in the output result.


## 4.3. Overall Algorithm

The overall algorithm comprises two phases. In the first stage, the Benders-and-Price approach alternates between RBSP and BMP with Benders cuts (18) generation until the gap between ηLB and ηUB is small enough. The BMP can be solved by IP solvers, like GUROBI, and RBSP is solved by column generation method (not embedded in B&B). If the final solution is integer valued in the end of the first phase algorithm, then the second phase can be eliminated. If not, the algorithm then enters the second phase, where it seeks integer solution of RBSP by performing a Branch-and-Price (BAP) algorithm. Note that, after the first phase, the timetable has been fixed. Figure 3 gives a flowchart for the overall algorithm.Figure 3
Flowchart of the overall algorithm.As we can see in the BAP algorithm, the branching rule which branches on the variablesxs(p) with fractional values near 1 and Depth-First node picking strategy with focus on “x=1” cuts are adopted. Actually, an exact branching scheme is to branch decisions on two compatible trips to be covered by the same vehicle [15, 16]. However, for large instances, this scheme is likely to require excessive computing time because of the large number of possible connections [23]. On the other hand, in the way of branching on the path variables, once one path (also a route) is fixed to be selected, all the paths which are conflicted with the fixed path will be eliminated. Moreover, recall that canceling trips is allowed, so infeasibility never occurs during the search. Then, it is hopeful to obtain a satisfactory upper bound. For a CSP imposed by an “x=0” cut from its parent subproblem, we refrain from solving it by column generation, but dual simplex method, because the type of “x=0” cut changes the solution space slightly, and using column generation to solve will result in the same paths being generated again, which hinders the convergence of algorithm. Although the final solution can be a suboptimal solution without solving each subproblem by column generation, the final solution is always high-quality and even optimal in our test cases. Last but not least, if xp∗ is a fractional solution, it is not always necessary to branch on the fractional variable, only when it is likely to save at least one vehicle, i.e., z(xp∗)≤UB-cR.


## 5. Numerical Experiments

### 5.1. Input Data and Parameter Settings

For simplicity, the proposed models and algorithms are verified on one of the bus lines in Beijing city of China. There are two terminals; thusK=1,2. The total number of vehicles u=30, and 20 of them can execute no more than 10 trips, 10 of which can execute no more than 8 trips. All the vehicles are equally assigned to two terminals; i.e., each terminal has 10 vehicles which can perform no more than 10 trips. For each trip i∈T, let Δdi-=Δdi+=ϕ = 5 min, the travel time li = 45 min. Additionally, we set h = 3 min, χmin = 10 min, and χmax = 60 min. The optimality gap ε is 0.1%. The cost factors λi, cr, and wk are set as 1, 100, and 200, respectively. All experiments are performed on a laptop equipped with 2.20 GHz Intel(R) Pentium(R) CPU and 8 GB memory, and all the algorithms are coded by C#. GUROBI 7.0.2 Optimizer is adopted as both the MIP and LP solvers.


### 5.2. Comparisons between GUROBI and the Proposed Approach

In this case study, we will compare the solving efficiency between using GUROBI to solve the connection-based model and using the proposed approach to solve the path-based model. Table3 shows the experiment results on five different scale instances. Each ordered triple from the first column represents the headway time in off-peak periods, the headway time in peak periods, and total number of trips, respectively. For example, the instance (10,8,228) denotes an initial timetable t0, for which the headway times are 10 min in off-peak periods and 8 min in peak periods, and it has 228 trips in total. As we can see, the connection-based model is difficult to solve by GUROBI. For example, in the last instance, the best feasible solution obtained by GUROBI in 3 h has an optimality gap 25.4%. Those demonstrate the limited capability of GUROBI solver on finding satisfactory solutions for somewhat large instances. Instead, the proposed approach can deal with the problem more efficiently in terms of solving time and solution quality.Table 3
Results and comparisons for different approaches.
Instance GUROBI The proposed approach Time Obj. Gap Time Obj. (80,60,30) 4 s 800 0 3 800 (60,40,42) 21 s 800 0 6 800 (40,20,68) 487 s 1000 0 9 1000 (20,10,132) 3 h 1610 17.1% 193 1610 (10,8,228) 3 h 3056 25.4% 3481 2470


### 5.3. Comparisons between the Fixed Timetable and Variable Timetable Approaches

In this section, we will solve all the instances by the designed approach which has two main procedures, namely, the Benders-and-Price approach and the BAP algorithm, with different values ofϕ. If ϕ=0, the fixed timetable (FT) approach is implemented, whereas ϕ>0 means the variable timetable (VT) approach is adopted. If ϕ=0 is set, i.e., timetable modification is not allowed, the path-based model is actually solved by the BAP algorithm without involving Benders decomposition. Table 4 reports the experimental results for FT approach and VT approach.Table 4
The solutions for FT approach and VT approach.
Instance FT approach VT approach Used vehicles Canceled trips Timetable deviation (min) Solving time (s) Used vehicles Canceled trips Timetable deviation (min) Solving time (s) (80,60,30) 4 2 0 3 4 2 0 3 (60,40,42) 8 0 0 5 8 0 0 6 (40,20,68) 10 0 0 9 10 0 0 9 (20,10,132) 17 0 0 690 16 0 10 193 (10,8,228) 26 0 0 2581 24 0 70 3481For the first three instances, the integer optimal solution can be found within 9 s without implementing the BAP algorithm. That is, after Benders-and-Price algorithm is terminated, the solutions obtained are integer valued. The solution quality in terms of the number of used vehicles does not have any improvement. For the last two instances, the final solutions are fractional in the first phase, so we have to perform the BAP algorithm to obtain integer solutions, which consumes most of the solving time. The number of used vehicles reduces by 1 and 2 in the instances (20,10,132) and (10,8,228), respectively. However, it should be borne in mind that not all the nodes in the BAP trees are solved by column generation in order to accelerate the convergence of BAP algorithm, which may make the final solution not optimal. That is, one should not come to the conclusion that the solution quality is improved only according to the required fleet size. Coincidently, both instances solved by VT approach have no branching procedures. Therefore, the solutions obtained by VT approach are optimal. As for FT approach, the lower bound of fleet size for FT approach (denoted asφLBFZ(FT), φLBFZ(VT) for VT approach correspondingly) is available, and if the actually required fleet size for FT approach (denoted as φARFZ(FT), φARFZ(VT) for VT approach correspondingly) in the final result is equal to φLBFZ(FT), namely, φARFZ(FT)=φLBFZ(FT), then the corresponding solution of FT approach is optimal. To derive φLBFZ(FT), let us see Figure 4, which shows the lower and upper bound of objective values at iterative process for VT approach, and timetable deviations are shown as the scatter plot. At first iteration, from the upper bound value 2533 (the optimal objective value of the root node of BAP tree in FT approach is the same as this value), we can deduce the φLBFZ(FT)=2533/cR=26, where · denotes the arithmetic that rounds up to the nearest integer. Furthermore, φARFZ(FT) = 26 (see Table 4) is the same as φLBFZ(FT), so the final solutions obtained by FT approach are also optimal. Thus, compared to FT approach, the solution quality is certainly improved by VT approach concerning the required fleet size and also the total cost. In addition, we can deduce that φLBFZ(VT)=(2470-70)/cR=24, which is equal to φARFZ(VT)=24 (see Table 4).Figure 4
Upper bounds and lower bounds of VT approach with timetable deviations at Benders-and-Price iterations.Figure5 displays the timetable with vehicle schedule solution of VT approach for the instance (10,8,228). In this figure, the routes are divided into two parts according to the origin terminal of the first trip in the corresponding route for clarity. In order to display the alterations of the departure times of trips, the trips whose departure times are changed in t0 are also depicted as dotted lines show.Figure 5
The schedule solution by VT approach.The local enlarged picture of the right upper corner (the hatched rectangle area) of Figure5 is shown in Figure 6, which gives more details about the alterations. For example, let us see a part of one route as the thick line shows. For the trips in the route, the departure time of trip 199 is 1 min ahead of schedule to satisfy the minimum turnaround time between trips 199 and 91, trip 211 is 5 min backward to satisfy the minimum turnaround time between trips 211 and 102, and trip 221 is 5 min forward to fulfill the minimum turnaround time between trips 102 and 221.Figure 6
The more details of modifications for the timetable in Figure5.


## 5.1. Input Data and Parameter Settings

For simplicity, the proposed models and algorithms are verified on one of the bus lines in Beijing city of China. There are two terminals; thusK=1,2. The total number of vehicles u=30, and 20 of them can execute no more than 10 trips, 10 of which can execute no more than 8 trips. All the vehicles are equally assigned to two terminals; i.e., each terminal has 10 vehicles which can perform no more than 10 trips. For each trip i∈T, let Δdi-=Δdi+=ϕ = 5 min, the travel time li = 45 min. Additionally, we set h = 3 min, χmin = 10 min, and χmax = 60 min. The optimality gap ε is 0.1%. The cost factors λi, cr, and wk are set as 1, 100, and 200, respectively. All experiments are performed on a laptop equipped with 2.20 GHz Intel(R) Pentium(R) CPU and 8 GB memory, and all the algorithms are coded by C#. GUROBI 7.0.2 Optimizer is adopted as both the MIP and LP solvers.


## 5.2. Comparisons between GUROBI and the Proposed Approach

In this case study, we will compare the solving efficiency between using GUROBI to solve the connection-based model and using the proposed approach to solve the path-based model. Table3 shows the experiment results on five different scale instances. Each ordered triple from the first column represents the headway time in off-peak periods, the headway time in peak periods, and total number of trips, respectively. For example, the instance (10,8,228) denotes an initial timetable t0, for which the headway times are 10 min in off-peak periods and 8 min in peak periods, and it has 228 trips in total. As we can see, the connection-based model is difficult to solve by GUROBI. For example, in the last instance, the best feasible solution obtained by GUROBI in 3 h has an optimality gap 25.4%. Those demonstrate the limited capability of GUROBI solver on finding satisfactory solutions for somewhat large instances. Instead, the proposed approach can deal with the problem more efficiently in terms of solving time and solution quality.Table 3
Results and comparisons for different approaches.
Instance GUROBI The proposed approach Time Obj. Gap Time Obj. (80,60,30) 4 s 800 0 3 800 (60,40,42) 21 s 800 0 6 800 (40,20,68) 487 s 1000 0 9 1000 (20,10,132) 3 h 1610 17.1% 193 1610 (10,8,228) 3 h 3056 25.4% 3481 2470


## 5.3. Comparisons between the Fixed Timetable and Variable Timetable Approaches

In this section, we will solve all the instances by the designed approach which has two main procedures, namely, the Benders-and-Price approach and the BAP algorithm, with different values ofϕ. If ϕ=0, the fixed timetable (FT) approach is implemented, whereas ϕ>0 means the variable timetable (VT) approach is adopted. If ϕ=0 is set, i.e., timetable modification is not allowed, the path-based model is actually solved by the BAP algorithm without involving Benders decomposition. Table 4 reports the experimental results for FT approach and VT approach.Table 4
The solutions for FT approach and VT approach.
Instance FT approach VT approach Used vehicles Canceled trips Timetable deviation (min) Solving time (s) Used vehicles Canceled trips Timetable deviation (min) Solving time (s) (80,60,30) 4 2 0 3 4 2 0 3 (60,40,42) 8 0 0 5 8 0 0 6 (40,20,68) 10 0 0 9 10 0 0 9 (20,10,132) 17 0 0 690 16 0 10 193 (10,8,228) 26 0 0 2581 24 0 70 3481For the first three instances, the integer optimal solution can be found within 9 s without implementing the BAP algorithm. That is, after Benders-and-Price algorithm is terminated, the solutions obtained are integer valued. The solution quality in terms of the number of used vehicles does not have any improvement. For the last two instances, the final solutions are fractional in the first phase, so we have to perform the BAP algorithm to obtain integer solutions, which consumes most of the solving time. The number of used vehicles reduces by 1 and 2 in the instances (20,10,132) and (10,8,228), respectively. However, it should be borne in mind that not all the nodes in the BAP trees are solved by column generation in order to accelerate the convergence of BAP algorithm, which may make the final solution not optimal. That is, one should not come to the conclusion that the solution quality is improved only according to the required fleet size. Coincidently, both instances solved by VT approach have no branching procedures. Therefore, the solutions obtained by VT approach are optimal. As for FT approach, the lower bound of fleet size for FT approach (denoted asφLBFZ(FT), φLBFZ(VT) for VT approach correspondingly) is available, and if the actually required fleet size for FT approach (denoted as φARFZ(FT), φARFZ(VT) for VT approach correspondingly) in the final result is equal to φLBFZ(FT), namely, φARFZ(FT)=φLBFZ(FT), then the corresponding solution of FT approach is optimal. To derive φLBFZ(FT), let us see Figure 4, which shows the lower and upper bound of objective values at iterative process for VT approach, and timetable deviations are shown as the scatter plot. At first iteration, from the upper bound value 2533 (the optimal objective value of the root node of BAP tree in FT approach is the same as this value), we can deduce the φLBFZ(FT)=2533/cR=26, where · denotes the arithmetic that rounds up to the nearest integer. Furthermore, φARFZ(FT) = 26 (see Table 4) is the same as φLBFZ(FT), so the final solutions obtained by FT approach are also optimal. Thus, compared to FT approach, the solution quality is certainly improved by VT approach concerning the required fleet size and also the total cost. In addition, we can deduce that φLBFZ(VT)=(2470-70)/cR=24, which is equal to φARFZ(VT)=24 (see Table 4).Figure 4
Upper bounds and lower bounds of VT approach with timetable deviations at Benders-and-Price iterations.Figure5 displays the timetable with vehicle schedule solution of VT approach for the instance (10,8,228). In this figure, the routes are divided into two parts according to the origin terminal of the first trip in the corresponding route for clarity. In order to display the alterations of the departure times of trips, the trips whose departure times are changed in t0 are also depicted as dotted lines show.Figure 5
The schedule solution by VT approach.The local enlarged picture of the right upper corner (the hatched rectangle area) of Figure5 is shown in Figure 6, which gives more details about the alterations. For example, let us see a part of one route as the thick line shows. For the trips in the route, the departure time of trip 199 is 1 min ahead of schedule to satisfy the minimum turnaround time between trips 199 and 91, trip 211 is 5 min backward to satisfy the minimum turnaround time between trips 211 and 102, and trip 221 is 5 min forward to fulfill the minimum turnaround time between trips 102 and 221.Figure 6
The more details of modifications for the timetable in Figure5.


## 6. Conclusions and Future Works

Considering the inefficiency of using GUROBI to solve the connection-based model for VSPVT-LT when the size of instances is relatively large, this paper formulated a path-based model, which can be solved by a bespoke approach working in two phases. In the first phase, Benders decomposition is used to decompose the original model into a timetabling problem in BMP and a VSP in BSP. Moreover, the LP-relaxation of RBSP is solved by column generation. After the Benders-and-Price algorithm is finished, a BAP algorithm is utilized to obtain an integer solution. Additionally, a vertices-constrained shortest path algorithm modified from Bellman-Ford algorithm is proposed to deal with the pricing problems of column generation. The experimental results show that the VT approach can reduce the needed number of vehicles in comparison with FT approach on the test cases, and the proposed algorithm performs well to handle the VSPVT-LT. Our future researches will concentrate on the following main extensions:(1)
In our case studies, with the Benders cuts added to BMP, the model becomes more and more difficult to handle by optimization software. An efficient algorithm should be further developed to solve BMP, such as the three-phase method [24, 25], which has one more phase because it relaxes all integrality constraints in the BMP in first phase and reintroduces the integrality constraints to the BMP in the next phase.(2)
Regarding the decrease in the number of iterations of Benders decomposition, Pareto-optimal cut method [25] can be employed. The approach exploits the fact that for a degenerate RBSP there exists multiple optimal dual solutions. Thus, it is likely to select the dual solution that is the nearest to the interior of the BMP polyhedron.



---

*Source: 2781590-2019-01-06.xml*