# Optimal Electric Bus Scheduling under Travel Time Uncertainty: A Robust Model and Solution Method

Authors: Mengyan Jiang; Yi Zhang; Yi Zhang
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

With the increasing adoption of electric buses (e-buses), e-bus scheduling problem has become an essential part of transit operation planning. As e-buses have a limited battery capacity, e-bus scheduling problem aims to assign vehicles to timetabled service trips on the bus routes considering their charging demand. Affected by the dynamic operation environment, the travel time and energy consumption of the e-buses often display considerable randomness, resulting in unexpected trip start delays and battery energy shortages. In this paper, we addressed the e-bus scheduling problem under travel time uncertainty by robust optimization approaches. We consider the cardinality constrained uncertainty set to formulate a robust multidepot EVSP model considering trip time uncertainty and partial recharging. The model is developed based on the dynamic programming equations that we formulated for trip chain robustness checking. A branch-and-price (BP) algorithm is devised to generate provably high-quality solutions for large-scale instances. In the BP algorithm, an efficient label setting algorithm is developed to solve the robust resource-constrained shortest path subproblem. Comprehensive numerical experiments are conducted based on the bus routes in Shenzhen to demonstrate the effectiveness of the suggested methodology. The robustness of the schedules was evaluated through Monte Carlo simulation. The results show that the trip start delay and battery energy shortage caused by the travel time uncertainty can be effectively reduced at the expense of an increase in the operational cost. A trade-off should be made between the reduction in infeasibility rate and increase in operational cost to choose a proper uncertainty budget.


---

## Body

## 1. Introduction

In recent years, an increasing number of electric buses (e-buses) have been introduced into the public transit systems because of their environmental and social benefits such as reducing on-road pollution, energy saving, and better onboard experience [1]. However, as the driving range of the e-buses per charge is limited and recharging the battery is time-consuming, the vehicles need extra time for battery recharging during the operation. To ensure operational efficiency, a charging plan for the e-bus fleet is required which specifies the optimal time, location, and amount to charge. In the operation planning phase, the electric vehicle scheduling problem (EVSP) aims to assign timetabled trips to the e-buses while making optimal charging plans to minimize the operational cost.Two essential parameters taken as input to the EVSP are the expected travel time and energy consumption of the bus trips which have a great impact on the reliability of the schedule. In real-life cases, getting accurate estimations of the trip travel time is difficult because of the substantial variability of traffic condition, passenger demand, and driving conditions. Consequently, the energy consumption of a trip can also deviate from the estimated value considerably. To mitigate possible variations, on the operation planning side, the operators usually introduce a buffer time between consecutive trips to absorb small delays and a safe range for the battery state of charge (SoC). However, buffer time is usually very short and not able to protect the schedule against trip time variations effectively. If some trips experience longer travel time than the expected value, the delay will propagate along the trip chain and influence the punctuality of the consecutive trips. The original charging plan can also become infeasible as the charging time and vehicle battery SoC are also influenced by the trip travel time.With the aim to generate robust schedules against trip travel time variation and reduce the delays that cannot be absorbed by the buffer time, the stochastic vehicle scheduling model and dynamic rescheduling strategies are proposed in the literature [2–7]. Different from the exiting studies, we aim to tackle the travel time uncertainty in the EVSP by robust optimization methods. The optimized schedule can remain feasible regarding the trip start time and vehicle battery SoC when trip time varies in the budgeted uncertainty set. The main contribution of this study is as follows:(i)
We consider the cardinality constrained uncertainty set to formulate a robust multidepot EVSP (R-MD-EVSP) model considering trip time uncertainty and partial recharging. The R-MD-EVSP model is developed based on the dynamic programming (DP) equations that we formulated for individual trip chain robustness checking.(ii)
A branch-and-price (BP) algorithm is devised where the pricing subproblem, a robust resource-constrained shortest path problem, is solved by an efficient labeling algorithm based on the DP equations.(iii)
Comprehensive numerical experiments are conducted based on the cases in Shenzhen with hundreds of trips on multiple bus routes. The robustness of the robust solutions is verified through Monte Carlo simulations.The remainder of the paper is organized as follows. In the next section, a literature review on the EVSP is provided. Section3 proposes the DP equations for trip chain robustness checking and the robust EVSP formulation. The BP algorithm is described in Section 4. The results of the numerical experiments are presented in Section 5. We conclude the paper in Section 6.


## 2. Literature Review

In the operation planning stage of public transportation, the vehicle scheduling problem (VSP) aims at allocating transit vehicles to carry out timetabled trips with the objective to minimize the total number of vehicles used. Single-depot VSP is polynomial time solvable while multidepot VSP has proven to be NP-hard by Bertossi et al. [8]. With the wide adoption of e-buses in the transit system in recent years, EVSP arises under different scenarios regarding the charging technology, charging policy, and fleet composition. Li [9] and Yang et al. [10] developed models and algorithms for the EVSP under battery swapping modes. Considering the plug-in charging mode, Liu and Ceder [11] proposed a model based on the deficit function theory aimed at minimizing the number of vehicles and chargers. Mixed fleet EVSP was considered in [12–14]. Zhang et al. [15] addressed an EVSP considering the degradation of battery and nonlinear charging process. A tailored BP algorithm was devised to solve the problem. Their computational experiments showed that the optimized schedule can reduce the cost considerably which is mainly achieved by the substantial extension of battery life. Considering nonlinear charging process and multivehicle type, Zhang et al. [16] developed an MIP model for the EVSP with linear approximation of the nonlinear charging function. An ALNS heuristic was devised to solve the problem.The above studies assume a full charging policy, while allowing for partial charging increases the operational flexibility and also the problem complexity. Considering partial charging, Wen et al. [17] developed a multidepot EVSP and an adaptive large neighborhood search heuristic to solve instances with customized bus trips. van Kooten Niekerk et al. [18] proposed two models for single-depot scenarios with a different level of detail resembling the actual charging processes. A column generation algorithm was developed to solve the problem. Koháni and Kohánia [19] proposed a linear MIP model for a single-depot EVSP considering the location of the charging stations and the assignment of chargers to the vehicles. The model was solved by the standard solver. Li et al. [20] addressed EVSP along with the charger deployment problem. An adaptive genetic algorithm is designed to solve the problem. Yıldırım and Yıldız [21] considered an optimal fleet composition and scheduling problem. An IP-column-generation algorithm is proposed to solve the pricing subproblem.Charging scheduling problem aims to determine the time, amount, and specific charger for each e-bus to get recharged. In large-sized transit systems, due to the problem scale and complexity, charging scheduling often assumes that the e-bus operation schedule is given beforehand. Qin et al. [22] conducted a simulation analysis on the daily charging patterns and demand charges of a fleet of e-buses. An optimal charging strategy is identified to minimize demand charges. Wang et al. [23] developed a mixed-integer programming (MIP) model for the e-bus recharging scheduling based on a real-world case. The model aims to assign the charging station, chargers, and charging time to the e-buses to minimize the system operating costs, assuming given timetable and operation schedule. With the goal to minimize the electricity demand charges and energy charges, He et al. [24] proposed a modeling framework to determine the time to charge a vehicle and the actual charging power. Liu et al. [25] considered the e-bus charging scheduling problem under limited charging resources in the charging station. A column-generation-based algorithm is developed to solve the large-scale problem instances. The results show that the optimal charging plan can greatly reduce the charging cost compared with the uncontrolled charging.Most of the studies take deterministic trip time as the model input and assume that the energy discharge is proportional to the travel distance. To tackle the trip time uncertainty in the real world, on the input side, some methods were proposed to improve the accuracy of the trip time estimation based on historical running data [26, 27]. On the planning side, dynamic rescheduling strategies and stochastic VSP models were developed. Huisman et al. [3] introduced a dynamic vehicle scheduling approach to solve the VSP periodically with renewed estimation on the future trip time. He et al. [5] formulated a stochastic dynamic VSP and adopted an approximate dynamic programming approach where the objective function is approximated by a feed-forward neural network. Huisman et al. [3] developed a stochastic VSP model using typical disruption scenarios to minimize the expected total planned costs and costs caused by disruptions. Shen et al. [2] proposed a probabilistic network flow model assuming certain trip time distributions and developed a hybrid heuristic solution method. Considering the stochastic travel time, Tang et al. [6] proposed a stochastic model and a dynamic rescheduling paradigm for a single-depot EVSP with the full charging policy. A speed-energy-consumption relationship was incorporated in their models. Bie et al. [7] proposed a multiobjective stochastic e-bus scheduling model considering the variability of travel time and energy consumption. The numerical study showed that the optimized schedule can effectively protect against the accumulation of stochastic volatilities.Compared with the EVSP, more studies have considered the travel time uncertainty in the vehicle routing problem (VRP). VRP arises in the logistics context which aims to plan for vehicle routes to serve customers at different locations. To tackle the uncertainty related to customer demand and vehicle travel time, stochastic and robust models as well as dynamic reoptimization approaches were proposed in the literature. The advantage of robust optimization is that it only requires defined bounds on the data rather than the underlying distribution while maintaining the tractability of the problem. Different kinds of uncertainty sets have been considered in robust VRP (RVRP) with the typical ones including the hypercube, budgeted uncertainty sets, and ellipsoidal sets. The cardinality constrained set proposed by Bertsimas and Sim [28] is widely used which specifies an upper bound of the number of customer demand or the travel links that can attain their maximum value.Motivated by the travel time uncertainty in maritime transportation, Agra et al. [29] proposed two RVSP formulations. Dynamic programming recursive equations were proposed for path feasibility checking based on which path inequalities can be generated. They showed that it suffices to consider a subset of the extreme points of the uncertainty polytope and developed scenario reduction solution techniques. Salicru et al. [26] proposed a compact formulation for the RVRP under demand uncertainty with dynamic programming equations integrated. To solve the RVRP and its variants, exact algorithms were developed based on the branch-price-and-cut method by reformulating the robust problem into its deterministic counterpart [30, 31]. Local search-based heuristic frameworks with solution robustness checking procedure embedded in were devised in [32, 33] to solve large-sized instances. Pelletier et al. [34] considered a capacitated electric VRP with energy consumption uncertainty. They considered different kinds of uncertainty sets and proposed an LNS-based heuristic to solve the problem.From the above discussion, we can see that although many studies have focused on the EVSP, few of them have considered the uncertainty related to the travel time and energy consumption in the problem. To the best of our knowledge, there has been no research to cope with the trip time uncertainty of the EVSP using the robust optimization method. Our study proposed a robust EVSP model and solution method considering the travel time uncertainty under budgeted uncertainty set.


## 3. Robust Multidepot EVSP

### 3.1. Problem Description

A bus route in the transit network is defined by two end stops and a series of intermediate stops. On a bus route, a service trip starts from one end stop at a scheduled time and ends at the other end stop to carry passengers. The timetable of the route includes all the trips that should be carried out in one day with their scheduled start time from the start stop and expected end time at the end stop. On the electrified bus routes, these trips are carried out by a fleet of e-buses. Each e-bus undertakes a sequence of trips in a day which is called a trip chain. An operation schedule consists of the trip chains for all the e-buses. For the convenience of operation, bus depots are established close to the end stops of the routes for vehicle parking, maintenance, and charging if charging facilities are established. Based on this setting, the R-MD-EVSP is described as follows. The list of notations is provided in Table1.Table 1
Definitions of the sets and parameters for the R-MD-EVSP model.
SetsKSet of depots with indexkVtSet of trip nodesVSet of nodesASet of arcsParametersa0Earliest operation start timeb0Latest operation end timeaiScheduled departure time of trip nodeibiLatest departure time of trip nodeiτiTrip time of trip nodei, an uncertain parametertiNominal trip time of trip nodeit^iUpper trip time deviation fromtit˜iLower trip time deviation fromtiδitTrip time deviation of trip nodei, δit∈0,t^iεiTrip energy consumption of trip nodei, an uncertain parametereiNominal energy consumption of trip nodeie^iUpper energy consumption deviation fromeie˜iLower energy consumption deviation fromeiδieEnergy consumption deviation of trip nodei, δie∈0,e^isijEmpty travel time from nodei to node jeijEnergy consumption of the empty travel from nodei to node jcijCost of arci,j∈AcdDeparture delay cost per minuteceUnit charging costtpCharging preparation timersEnergy charging rateruEnergy consumption rateΓtUncertainty budget, maximum number of trips that can take their upper bound trip time, indexed byγumaxUpper bound of the vehicle battery SoCuminLower bound of the vehicle battery SoCuiγminLower bound of variableYiγ, γ≤ΓtLetG=V,A be an acyclic direct graph, where V is the set of nodes and A is the set of arcs. Each timetabled trip is represented by a trip node in the graph. Denote T⊆V as the set of trip nodes. Each trip node i∈T is associated with a trip start time window ai,bi, where ai and bi are the scheduled start time and latest start time, respectively. The e-bus serving a trip should wait until the scheduled start time ai if it arrives earlier at the trip start depot. Each trip node i∈T is associated with two uncertain parameters: the trip time τi and trip energy consumption εi. τi is realized in the interval ti−t˜i,ti+t^i where ti is the nominal trip time; t^i≥0 and t˜i≥0 are the upper and lower deviations from ti. εi is realized in the interval ei−e˜i,ei+e^i where ei is the nominal value; e^i≥0 and e˜i≥0 are the upper and lower deviations from ei. Since the worst case will always be achieved at the right-hand side of the interval, we only consider the realization of τi in ti,ti+t^i and εi in ei,ei+e^i correspondingly.DenoteK as the set of depots. Each depot k∈K is created with two nodes in the graph: operation start node ok and operation end node dk, indicating that the e-bus begins/ends the operation from/at depot k, respectively. For nodes ok and dk, k∈K, we assign a time window a0,b0, where ao is the earliest operation start time and b0 is the latest operation end time. This means that the e-buses can begin and end their operation at any time within time window a0,b0.The arc setA includes three kinds of arcs: (i) the pull-out arc connects a depot node ok, k∈K, and a trip node i∈T, representing that a vehicle begins operation from depot k to carry out the first trip i; (ii) the pull-in arc connects a trip node i∈T and an end node dk, k∈K, representing that a vehicle finishes its last trip i and ends the operation, returning to depot k; and (iii) the trip connection arc connects two trip nodes, representing that two trips are carried out consecutively. Each arc i,j∈A is associated with a nonservice travel time sij and energy consumption eij. To ensure that the trip start time window is respected, trip nodes i and j, i,j∈T, are connected only if the time compatible condition ai+ti+sij≤bj is satisfied. The cost of an arc cij includes the empty travel cost and the vehicle usage cost f for the pull-out arcs. Denote gij, i,j∈A, as the unit charging cost in the time interval between the end of trip node i and the start of trip node j; gij is calculated as the weighted average value of the time-of-use tariff in the time interval.Figure1 gives an example of the graph for an instance with two depots and four trips. Each depot is created with a start node ok and an end node dk, k=1,2. Each trip is represented by a trip node vi associated with a departure time window and a scheduled trip time, i=1,2,3,4.Figure 1
An illustrative example of the graph for the robust multidepot EVSP.We assume that the e-buses begin the operation with a fully charged battery and get recharged at the destination depot after finishing a trip if needed. The charging amount is proportional to the charging time with a rate ofrs (kWh/min). Denote ts as the unit charging time and tp as the charging preparation time. The e-buses should charge integer times of ts.The historical running data of 150 bus routes in five months in Shenzhen were used to analyze the relationship between the trip time and trip energy consumption. Most of them display similar relationship. Figure2 shows the relationship between the average trip time and energy consumption of the selected four bus routes. The solid lines are the linear regression results of the sample average data with the R-square values given next to it. The R-square values are both above 0.9, meaning that the trip energy consumption can be described as a linear function of the trip time, although the data of Route 3 and Route 4 display some nonlinear trends. The trip time and trip energy consumption may display different relationships in different cities and routes. Based on the cases of Shenzhen, we assume that for a trip i, the energy consumption deviation δie∈0,e^i is proportional to the trip time deviation δit∈0,t^i, i.e., δie=αδitru where α≥1 is an augment constant and ru is the energy consumption rate (kWh/min). Note that the robust EVSP model developed in this study does not depend on the assumption of a linear relationship between trip time and trip energy consumption. We assume that the energy consumption of a trip i varies within interval ei,ei+e^i; as such, we only need to give the values of ei and e^i as the model inputs. Therefore, any empirical method that can give estimation on the values of ei and e^i works.Figure 2
Relationship between the trip time and energy consumption for the selected bus routes. (a) Route 1. (b) Route 2. (c) Route 3. (d) Route 4.
(a)(b)(c)(d)We adopt a cardinality constrained uncertainty setU as discussed in [24] to control the level of trip time uncertainty. U is defined as U=τ∈ℝT|τi=ti+θit^i,0≤θi≤1,∀i∈T,∑i∈Tθi≤Γt where the cumulative uncertainty of the random variable θi is bounded by the budget Γt. We choose this uncertainty set because it is commonly used in the RVRP, the solution structure of which is similar to our problem: the trip chain assigned to an e-bus is equivalent to the customer route assigned to a vehicle in the VRP, where a bus trip can be regarded as a customer.Based on the above settings, the problem looks for the minimum cost schedule that satisfies the following constraints: (i) each trip is carried out exactly once; (ii) the departure time window of each trip is respected when the trip time varies within the budgeted setU; (iii) each vehicle begins and ends its operation at the same depot; and (iv) the vehicle battery SoC is always kept within the safety range umin,umax when the trip time varies within the budgeted set U. Following the robust optimization paradigm, a trip chain is robust feasible if it satisfies conditions (ii)–(iv) for all possible trip time realizations. A feasible schedule of the problem is given as set of feasible trip chains that together also satisfies condition (i).


### 3.2. Deterministic Problem Formulation

LetXijk be the binary variables that take value 1 if a vehicle housed at terminal k traverses arc i,j∈A. At each trip node i∈T, denote Zi as the trip start time and Yi as the lowest battery SoC of the vehicle at the trip start time. The lower bound of Yi is denoted as uimin=umin+ei. At each operation start (end) node ok (dk), k∈K, denote Zi as the operation start (end) time and Yi as the lowest battery SoC of the vehicle at the operation start (end) time. On each arc i,j∈A, i,j∈T, denote Wij as the amount of energy to be charged after completing trip i before starting trip j. According to our assumption, the charging place is at the ending depot of trip node i. Rij are the auxiliary variables defined corresponding to Zi for constraint linearization. Denote Eij as the binary variables that are equal to 1 if the e-buses can be charged on arc i,j∈A. Denote Uij as the integer variables that specify the number of time units an e-bus spends charging on arc i,j∈A. The deterministic problem is formulated as follows (MD-EVSP):(1)min∑k∈K∑i,j∈AcijXijk+∑i,j∈AgijWij,(2)s.t.∑k∈K∑j:i,j∈AXijk=1,i∈T,(3)∑j:i,j∈AXijk−∑j:j,i∈AXjik=0,i∈T,k∈K,(4)∑β∈K\k∑i∈VXoβik=0,k∈K,(5)∑β∈K\k∑i∈VXidβk=0,k∈K,(6)ai≤Zi≤bi,i∈V,(7)Zj−Zi−ti−sij+M1−∑k∈KXijk≥0,i,j∈A,(8)Yi−ei+Wij−eij+M1−∑k∈KXijk≥Yj,i,j∈A,(9)Eij−∑k∈KXijk≤0,i,j∈A,(10)−aj−Zi−ti−sij−tp−M1−Eij≤0,i,j∈A,(11)Wij+Rijrs−aj−ti−sij−tprsEij≤0,i,j∈A,(12)Wij+Yi−ei−umax≤0,i,j∈A,(13)Wij−rstsUij=0,i,j∈A,(14)Rij−biEij≤0,i,j∈A,(15)aiEij−Rij≤0,i,j∈A,(16)Rij−Zi+ai1−Eij≤0,i,j∈A,(17)Zi−bi1−Eij−Rij≤0,i,j∈A,(18)uimin≤Yi≤umax,i∈V,(19)Wij≥0,i,j∈A,(20)Xijk∈0,1,i,j∈A,k∈K,(21)Eij∈0,1,i,j∈A,(22)Rij≥0,i,j∈A,(23)Uij∈Z,i,j∈A.Objective (1) minimizes the sum of vehicle usage, empty travel, and charging cost. Constraint (2) ensures that each trip is carried out only once. Constraint (3) is the flow conservation constraint. Constraints (4) and (5) require that each vehicle begin and end its operation at its base depot. Constraint (6) limits the departure time of each trip to be within the departure time window. Constraint (7) ensures the time consistency of two successive trips. Constraint (8) ensures the vehicle battery energy consistency of two successive trips. Constraint (9) requires that charging on arc i,j is performed only if trip j is succeeded by trip i. Constraint (10) restricts that charging on arc i,j can be performed only if time is available. Constraint (11) defines the upper bound of Wij on arc i,j∈A according to the time available. Originally, constraint (11) is formulated as follows:(24)Wij+Zi−aj−ti−sij−tprsEij≤0,i,j∈A.To linearize quadratic termsZiEij, auxiliary variables Rij are introduced where Rij=ZiEij is satisfied by constraints (14)–(17). Constraint (12) requires the battery SoC after charging not to exceed umax. Constraint (13) ensures that the charging time can only be integer times of the minimum charging time unit. Constraint (18) keeps the vehicle battery SoC within the safety range. Constraints (19)–(23) define the domains of the variables.


### 3.3. DP Equations to Define the Robustness of a Trip Chain

In this section, we introduce the DP equations to define the robustness of a trip chain under uncertainty setU. In graph G, a trip chain is represented by a path P=ok,v1,v2,…,vn,dk starting from an operation start node ok to an operation end node dk, visiting a sequence of trip nodes vi,i=1,2,…,n. In solving the RVRP, Agra et al. [29] showed that only the travel time vector τ∈extU needs to be considered in the problem formulation to ensure the robustness feasibility. extU is the set that contains all the extreme points of set U. Because of the structure of set U, they defined DP recursive equations to define the robustness of a vehicle route. As our R-MD-EVSP problem shares a similar solution structure with the RVRP, we can define the robustness of a trip chain by the DP recursive equations (25) and (26) as follows.DenoteZiγ as the earliest possible start time of node vi and Yiγ as the lowest battery SoC of the vehicle upon the start of vi; let Wijγ and W˜ijγ be the charging amount on arc i,j, i.e., after the end of trip node vi and before the start of trip node vj, when γ≤Γt trips from v1 to vi are taking their maximum trip time, among which trip node vi takes the nominal trip time and maximum trip time, respectively. Denote uiγmin as the lower bound of Yiγ: if γ<Γt, uiγmin=umin+ei+e^i; otherwise, uiγmin=umin+ei.Ziγ and Yiγ are defined by the recursive functions (25) and (26), respectively. A trip chain is robust if Ziγ∈ai,bi and Yiγ∈uiγmin,umax, for ∀vi∈P and γ≤Γt. In fact, it suffices to check ZiΓt and YiΓt for vi∈P.(25)Zjγ=a0maxaj,Ziγ+ti+tijmaxaj,Ziγ+ti+Wijγrs+tij,Zi,γ−1+ti+t^i+W˜ijγrs+tijj=0γ=0otherwise,(26)Yjγ=umaxj=ominumax,Yiγ−ei+Wijγ−eijγ=0minminumax,Yiγ−ei+Wijγ−eij,otherwiseminumax,Yi,γ−1−ei−e^i+W˜ijγ−eij.


#### 3.3.1. The Optimal Total Charging Amount and Charging Policy

DenoteH as the optimal total charging amount of a trip chain. H is the maximum total charging amount required when the trip times varies within uncertainty set U which is calculated as follows: let b=minn,Γt be the maximum number of trips of path P that can attain their maximum trip time. Choose the b trips with the maximum worst case trip time and take their maximum trip times. The rest of the trips take their nominal trip times. Then, H is given by H=umin+E−umax where E is the total energy consumption of P. The charging policy for an e-bus carrying out trip chain P is as follows: the e-buses can get recharged during the layover time between two successive trips for any flexible amount that is integer time of the unit charging time ts until the total charging amount reaches H.


### 3.4. Robust Problem Formulation

We introduce an MIP formulation for the R-MD-EVSP under cardinality uncertainty setU by incorporating the DP equations for trip chain robustness checking into the deterministic model.At each nodei∈V, given γ trips have taken their maximum trip time, γ=0,1,2,…,Γt, variables Ziγ, Yiγ, Wijγ, W˜ijγ and parameter uiγmin are consistent with those defined in Section 3.3. Given that γ trips on the trip chain from the first trip to trip j (trip j not included) are taking their maximum trip time, denote Eijγ and E˜ijγ as the binary variables that are equal to 1 if vehicle charging can be carried out on arc i,j∈A when trip i takes its nominal and maximum trip time, respectively; denote Uijγ and U˜ijγ as the integer variables that specify the number of time units an e-bus spends charging on arc i,j∈A when trip i takes its nominal and maximum trip time, respectively. Denote Rijγ and R˜ijγ as the auxiliary variables defined to linearize the quadratic terms in the constraints.In the robust model, the robustness of a trip chain is ensured by linearizing of the DP equations (25) and (26) and adding them into the problem constraints. The following constraints ((27)–(30)) are the linearized forms of the DP equations (25) and (26).(i)
Time consistency constraints:(27)Zjγ−Ziγ−ti−sij+M1−∑k∈KXijk≥0,i,j∈Ap,γ=0,1,2,…,Γt,(28)Zjγ−Zi,γ−1−ti−t^i−sij+M1−∑k∈KXijk≥0,i,j∈Ap,γ=1,2,…,Γt.Constraints (27) and (28) are the robust constraints corresponding to constraint (7). It is developed based on the linearization of equation (25), which means the time consistency of two successive trips should be satisfied for ∀γ=0,1,2,…,Γt.(ii)
Battery SoC consistency constraints:(29)Yiγ−ei+Wijγ−eij+M1−∑k∈KXijk≥Yjγ,i,j∈A,γ=0,1,2,…,Γt,(30)Yi,γ−1−ei−e^i+W˜ijγ−eij+M1−∑k∈KXijk≥Yjγ,i,j∈A,γ=1,2,…,Γt.Constraints (29) and (30) are the robust constraints corresponding to constraint (8). It is developed based on the linearization of equation (26), which means the battery SoC consistency of two successive trips should be satisfied for ∀γ=0,1,2,…,Γt.Other constraints are defined corresponding to the deterministic models with the variables substituted by the ones defined in the robust model. The R-MD-EVSP model is formulated as follows:(31)min∑k∈K∑i,j∈AcijXijk+∑γ=0Γt∑i,j∈AgijWijγ+∑γ=1Γt∑i,j∈AgijW˜ijγ,s.t. constraints (2)–(5) and constraints (27)–(30).(32)ai≤Ziγ≤bi,i∈V,γ=0,1,2,…,Γt,(33)Eijγ−∑k∈KXijk≤0,i,j∈A,γ=0,1,2,…,Γt,(34)E˜ijγ−∑k∈KXijk≤0,i,j∈A,γ=1,2,…,Γt,(35)−aj−Ziγ−ti−sij−tp−M1−Eijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(36)−aj−Zi,γ−1−ti−t^i−sij−tp−M1−E˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(37)Wijγ+Rijγrs−aj−ti−sij−tprsEijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(38)W˜ijγ+R˜ijγrs−aj−ti−t^i−sij−tprsE˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(39)Wijγ+Yiγ−ei−umax≤0,i,j∈A,γ=0,1,2,…,Γt,(40)W˜ijγ+Yi,γ−1−ei−e^i−umax≤0,i,j∈A,γ=1,2,…,Γt,(41)Wij−rstsUijγ=0,i,j∈A,γ=0,1,2,…,Γt,(42)Wij−rstsU˜ijγ=0,i,j∈A,γ=0,1,2,…,Γt,(43)Rijγ−biEijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(44)aiEijγ−Rijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(45)Rijγ−Ziγ+ai1−Eijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(46)Ziγ−bi1−Eijγ−Rijγ≤0,i,j∈Ap,γ=0,1,2,…,Γt,(47)R˜ijγ−biE˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(48)aiE˜ijγ−R˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(49)R˜ijγ−Zi,γ−1+ai1−E˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(50)Zi,γ−1−bi1−E˜ijγ−R˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(51)uiγmin≤Yiγ≤umax,i∈V,γ=0,1,2,…,Γt,(52)Wijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(53)W˜ijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(54)Xijk∈0,1,i,j∈A,k∈K,(55)Eijγ∈0,1,i,j∈A,γ=0,1,2,…,Γt,(56)Rijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(57)Uijγ∈Z,i,j∈A,γ=0,1,2,…,Γt,(58)E˜ijγ∈0,1,i,j∈A,γ=0,1,2,…,Γt,(59)R˜ijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(60)U˜ijγ∈Z,i,j∈A,γ=0,1,2,…,Γt.Objective (31) consists of the vehicle usage, empty travel, and the sum of charging cost in all possible cases regarding the realizations of the trip times. The worst-case charging cost cannot be expressed explicitly; therefore, it is not included in the objective but calculated after the optimal solution is obtained. Constraint (32) limits the departure time of each trip to be within the departure time window. Constraints (33) and (34) require that charging on arc i,j is performed only if trip j is succeeded by trip i. Constraints (35) and (36) restrict that charging on arc i,j can be performed if time is available. Constraints (37) and (38) define the upper bound of Wij on arc i,j∈A according to the available time. Originally, constraints (37) and (38) are formulated as follows:(61)Wijγ+Ziγ−aj−ti−sij−tprsEijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(62)W˜ijγ+Zi,γ−1−aj−ti−t^i−sij−tprsE˜ijγ≤0,i,j∈A,γ=1,2,…,Γt.To linearize quadratic termsZiγEijγ and Zi,γ−1E˜ijγ, auxiliary variables Rijγ and R˜ijγ are introduced where Rijγ=ZiγEijγ is satisfied by constraints (43)–(46) and R˜ijγ=Zi,γ−1E˜ijγ is satisfied by constraints (47)–(50).Constraints (39) and (40) require the battery SoC after charging not to exceed umax. Constraints (41) and (42) ensure that the charging time can only be integer times of the minimum charging time unit. Constraint (49) keeps the vehicle battery SoC within the safety range. Constraints (52)–(60) define the domains of the variables.


## 3.1. Problem Description

A bus route in the transit network is defined by two end stops and a series of intermediate stops. On a bus route, a service trip starts from one end stop at a scheduled time and ends at the other end stop to carry passengers. The timetable of the route includes all the trips that should be carried out in one day with their scheduled start time from the start stop and expected end time at the end stop. On the electrified bus routes, these trips are carried out by a fleet of e-buses. Each e-bus undertakes a sequence of trips in a day which is called a trip chain. An operation schedule consists of the trip chains for all the e-buses. For the convenience of operation, bus depots are established close to the end stops of the routes for vehicle parking, maintenance, and charging if charging facilities are established. Based on this setting, the R-MD-EVSP is described as follows. The list of notations is provided in Table1.Table 1
Definitions of the sets and parameters for the R-MD-EVSP model.
SetsKSet of depots with indexkVtSet of trip nodesVSet of nodesASet of arcsParametersa0Earliest operation start timeb0Latest operation end timeaiScheduled departure time of trip nodeibiLatest departure time of trip nodeiτiTrip time of trip nodei, an uncertain parametertiNominal trip time of trip nodeit^iUpper trip time deviation fromtit˜iLower trip time deviation fromtiδitTrip time deviation of trip nodei, δit∈0,t^iεiTrip energy consumption of trip nodei, an uncertain parametereiNominal energy consumption of trip nodeie^iUpper energy consumption deviation fromeie˜iLower energy consumption deviation fromeiδieEnergy consumption deviation of trip nodei, δie∈0,e^isijEmpty travel time from nodei to node jeijEnergy consumption of the empty travel from nodei to node jcijCost of arci,j∈AcdDeparture delay cost per minuteceUnit charging costtpCharging preparation timersEnergy charging rateruEnergy consumption rateΓtUncertainty budget, maximum number of trips that can take their upper bound trip time, indexed byγumaxUpper bound of the vehicle battery SoCuminLower bound of the vehicle battery SoCuiγminLower bound of variableYiγ, γ≤ΓtLetG=V,A be an acyclic direct graph, where V is the set of nodes and A is the set of arcs. Each timetabled trip is represented by a trip node in the graph. Denote T⊆V as the set of trip nodes. Each trip node i∈T is associated with a trip start time window ai,bi, where ai and bi are the scheduled start time and latest start time, respectively. The e-bus serving a trip should wait until the scheduled start time ai if it arrives earlier at the trip start depot. Each trip node i∈T is associated with two uncertain parameters: the trip time τi and trip energy consumption εi. τi is realized in the interval ti−t˜i,ti+t^i where ti is the nominal trip time; t^i≥0 and t˜i≥0 are the upper and lower deviations from ti. εi is realized in the interval ei−e˜i,ei+e^i where ei is the nominal value; e^i≥0 and e˜i≥0 are the upper and lower deviations from ei. Since the worst case will always be achieved at the right-hand side of the interval, we only consider the realization of τi in ti,ti+t^i and εi in ei,ei+e^i correspondingly.DenoteK as the set of depots. Each depot k∈K is created with two nodes in the graph: operation start node ok and operation end node dk, indicating that the e-bus begins/ends the operation from/at depot k, respectively. For nodes ok and dk, k∈K, we assign a time window a0,b0, where ao is the earliest operation start time and b0 is the latest operation end time. This means that the e-buses can begin and end their operation at any time within time window a0,b0.The arc setA includes three kinds of arcs: (i) the pull-out arc connects a depot node ok, k∈K, and a trip node i∈T, representing that a vehicle begins operation from depot k to carry out the first trip i; (ii) the pull-in arc connects a trip node i∈T and an end node dk, k∈K, representing that a vehicle finishes its last trip i and ends the operation, returning to depot k; and (iii) the trip connection arc connects two trip nodes, representing that two trips are carried out consecutively. Each arc i,j∈A is associated with a nonservice travel time sij and energy consumption eij. To ensure that the trip start time window is respected, trip nodes i and j, i,j∈T, are connected only if the time compatible condition ai+ti+sij≤bj is satisfied. The cost of an arc cij includes the empty travel cost and the vehicle usage cost f for the pull-out arcs. Denote gij, i,j∈A, as the unit charging cost in the time interval between the end of trip node i and the start of trip node j; gij is calculated as the weighted average value of the time-of-use tariff in the time interval.Figure1 gives an example of the graph for an instance with two depots and four trips. Each depot is created with a start node ok and an end node dk, k=1,2. Each trip is represented by a trip node vi associated with a departure time window and a scheduled trip time, i=1,2,3,4.Figure 1
An illustrative example of the graph for the robust multidepot EVSP.We assume that the e-buses begin the operation with a fully charged battery and get recharged at the destination depot after finishing a trip if needed. The charging amount is proportional to the charging time with a rate ofrs (kWh/min). Denote ts as the unit charging time and tp as the charging preparation time. The e-buses should charge integer times of ts.The historical running data of 150 bus routes in five months in Shenzhen were used to analyze the relationship between the trip time and trip energy consumption. Most of them display similar relationship. Figure2 shows the relationship between the average trip time and energy consumption of the selected four bus routes. The solid lines are the linear regression results of the sample average data with the R-square values given next to it. The R-square values are both above 0.9, meaning that the trip energy consumption can be described as a linear function of the trip time, although the data of Route 3 and Route 4 display some nonlinear trends. The trip time and trip energy consumption may display different relationships in different cities and routes. Based on the cases of Shenzhen, we assume that for a trip i, the energy consumption deviation δie∈0,e^i is proportional to the trip time deviation δit∈0,t^i, i.e., δie=αδitru where α≥1 is an augment constant and ru is the energy consumption rate (kWh/min). Note that the robust EVSP model developed in this study does not depend on the assumption of a linear relationship between trip time and trip energy consumption. We assume that the energy consumption of a trip i varies within interval ei,ei+e^i; as such, we only need to give the values of ei and e^i as the model inputs. Therefore, any empirical method that can give estimation on the values of ei and e^i works.Figure 2
Relationship between the trip time and energy consumption for the selected bus routes. (a) Route 1. (b) Route 2. (c) Route 3. (d) Route 4.
(a)(b)(c)(d)We adopt a cardinality constrained uncertainty setU as discussed in [24] to control the level of trip time uncertainty. U is defined as U=τ∈ℝT|τi=ti+θit^i,0≤θi≤1,∀i∈T,∑i∈Tθi≤Γt where the cumulative uncertainty of the random variable θi is bounded by the budget Γt. We choose this uncertainty set because it is commonly used in the RVRP, the solution structure of which is similar to our problem: the trip chain assigned to an e-bus is equivalent to the customer route assigned to a vehicle in the VRP, where a bus trip can be regarded as a customer.Based on the above settings, the problem looks for the minimum cost schedule that satisfies the following constraints: (i) each trip is carried out exactly once; (ii) the departure time window of each trip is respected when the trip time varies within the budgeted setU; (iii) each vehicle begins and ends its operation at the same depot; and (iv) the vehicle battery SoC is always kept within the safety range umin,umax when the trip time varies within the budgeted set U. Following the robust optimization paradigm, a trip chain is robust feasible if it satisfies conditions (ii)–(iv) for all possible trip time realizations. A feasible schedule of the problem is given as set of feasible trip chains that together also satisfies condition (i).


## 3.2. Deterministic Problem Formulation

LetXijk be the binary variables that take value 1 if a vehicle housed at terminal k traverses arc i,j∈A. At each trip node i∈T, denote Zi as the trip start time and Yi as the lowest battery SoC of the vehicle at the trip start time. The lower bound of Yi is denoted as uimin=umin+ei. At each operation start (end) node ok (dk), k∈K, denote Zi as the operation start (end) time and Yi as the lowest battery SoC of the vehicle at the operation start (end) time. On each arc i,j∈A, i,j∈T, denote Wij as the amount of energy to be charged after completing trip i before starting trip j. According to our assumption, the charging place is at the ending depot of trip node i. Rij are the auxiliary variables defined corresponding to Zi for constraint linearization. Denote Eij as the binary variables that are equal to 1 if the e-buses can be charged on arc i,j∈A. Denote Uij as the integer variables that specify the number of time units an e-bus spends charging on arc i,j∈A. The deterministic problem is formulated as follows (MD-EVSP):(1)min∑k∈K∑i,j∈AcijXijk+∑i,j∈AgijWij,(2)s.t.∑k∈K∑j:i,j∈AXijk=1,i∈T,(3)∑j:i,j∈AXijk−∑j:j,i∈AXjik=0,i∈T,k∈K,(4)∑β∈K\k∑i∈VXoβik=0,k∈K,(5)∑β∈K\k∑i∈VXidβk=0,k∈K,(6)ai≤Zi≤bi,i∈V,(7)Zj−Zi−ti−sij+M1−∑k∈KXijk≥0,i,j∈A,(8)Yi−ei+Wij−eij+M1−∑k∈KXijk≥Yj,i,j∈A,(9)Eij−∑k∈KXijk≤0,i,j∈A,(10)−aj−Zi−ti−sij−tp−M1−Eij≤0,i,j∈A,(11)Wij+Rijrs−aj−ti−sij−tprsEij≤0,i,j∈A,(12)Wij+Yi−ei−umax≤0,i,j∈A,(13)Wij−rstsUij=0,i,j∈A,(14)Rij−biEij≤0,i,j∈A,(15)aiEij−Rij≤0,i,j∈A,(16)Rij−Zi+ai1−Eij≤0,i,j∈A,(17)Zi−bi1−Eij−Rij≤0,i,j∈A,(18)uimin≤Yi≤umax,i∈V,(19)Wij≥0,i,j∈A,(20)Xijk∈0,1,i,j∈A,k∈K,(21)Eij∈0,1,i,j∈A,(22)Rij≥0,i,j∈A,(23)Uij∈Z,i,j∈A.Objective (1) minimizes the sum of vehicle usage, empty travel, and charging cost. Constraint (2) ensures that each trip is carried out only once. Constraint (3) is the flow conservation constraint. Constraints (4) and (5) require that each vehicle begin and end its operation at its base depot. Constraint (6) limits the departure time of each trip to be within the departure time window. Constraint (7) ensures the time consistency of two successive trips. Constraint (8) ensures the vehicle battery energy consistency of two successive trips. Constraint (9) requires that charging on arc i,j is performed only if trip j is succeeded by trip i. Constraint (10) restricts that charging on arc i,j can be performed only if time is available. Constraint (11) defines the upper bound of Wij on arc i,j∈A according to the time available. Originally, constraint (11) is formulated as follows:(24)Wij+Zi−aj−ti−sij−tprsEij≤0,i,j∈A.To linearize quadratic termsZiEij, auxiliary variables Rij are introduced where Rij=ZiEij is satisfied by constraints (14)–(17). Constraint (12) requires the battery SoC after charging not to exceed umax. Constraint (13) ensures that the charging time can only be integer times of the minimum charging time unit. Constraint (18) keeps the vehicle battery SoC within the safety range. Constraints (19)–(23) define the domains of the variables.


## 3.3. DP Equations to Define the Robustness of a Trip Chain

In this section, we introduce the DP equations to define the robustness of a trip chain under uncertainty setU. In graph G, a trip chain is represented by a path P=ok,v1,v2,…,vn,dk starting from an operation start node ok to an operation end node dk, visiting a sequence of trip nodes vi,i=1,2,…,n. In solving the RVRP, Agra et al. [29] showed that only the travel time vector τ∈extU needs to be considered in the problem formulation to ensure the robustness feasibility. extU is the set that contains all the extreme points of set U. Because of the structure of set U, they defined DP recursive equations to define the robustness of a vehicle route. As our R-MD-EVSP problem shares a similar solution structure with the RVRP, we can define the robustness of a trip chain by the DP recursive equations (25) and (26) as follows.DenoteZiγ as the earliest possible start time of node vi and Yiγ as the lowest battery SoC of the vehicle upon the start of vi; let Wijγ and W˜ijγ be the charging amount on arc i,j, i.e., after the end of trip node vi and before the start of trip node vj, when γ≤Γt trips from v1 to vi are taking their maximum trip time, among which trip node vi takes the nominal trip time and maximum trip time, respectively. Denote uiγmin as the lower bound of Yiγ: if γ<Γt, uiγmin=umin+ei+e^i; otherwise, uiγmin=umin+ei.Ziγ and Yiγ are defined by the recursive functions (25) and (26), respectively. A trip chain is robust if Ziγ∈ai,bi and Yiγ∈uiγmin,umax, for ∀vi∈P and γ≤Γt. In fact, it suffices to check ZiΓt and YiΓt for vi∈P.(25)Zjγ=a0maxaj,Ziγ+ti+tijmaxaj,Ziγ+ti+Wijγrs+tij,Zi,γ−1+ti+t^i+W˜ijγrs+tijj=0γ=0otherwise,(26)Yjγ=umaxj=ominumax,Yiγ−ei+Wijγ−eijγ=0minminumax,Yiγ−ei+Wijγ−eij,otherwiseminumax,Yi,γ−1−ei−e^i+W˜ijγ−eij.


### 3.3.1. The Optimal Total Charging Amount and Charging Policy

DenoteH as the optimal total charging amount of a trip chain. H is the maximum total charging amount required when the trip times varies within uncertainty set U which is calculated as follows: let b=minn,Γt be the maximum number of trips of path P that can attain their maximum trip time. Choose the b trips with the maximum worst case trip time and take their maximum trip times. The rest of the trips take their nominal trip times. Then, H is given by H=umin+E−umax where E is the total energy consumption of P. The charging policy for an e-bus carrying out trip chain P is as follows: the e-buses can get recharged during the layover time between two successive trips for any flexible amount that is integer time of the unit charging time ts until the total charging amount reaches H.


## 3.3.1. The Optimal Total Charging Amount and Charging Policy

DenoteH as the optimal total charging amount of a trip chain. H is the maximum total charging amount required when the trip times varies within uncertainty set U which is calculated as follows: let b=minn,Γt be the maximum number of trips of path P that can attain their maximum trip time. Choose the b trips with the maximum worst case trip time and take their maximum trip times. The rest of the trips take their nominal trip times. Then, H is given by H=umin+E−umax where E is the total energy consumption of P. The charging policy for an e-bus carrying out trip chain P is as follows: the e-buses can get recharged during the layover time between two successive trips for any flexible amount that is integer time of the unit charging time ts until the total charging amount reaches H.


## 3.4. Robust Problem Formulation

We introduce an MIP formulation for the R-MD-EVSP under cardinality uncertainty setU by incorporating the DP equations for trip chain robustness checking into the deterministic model.At each nodei∈V, given γ trips have taken their maximum trip time, γ=0,1,2,…,Γt, variables Ziγ, Yiγ, Wijγ, W˜ijγ and parameter uiγmin are consistent with those defined in Section 3.3. Given that γ trips on the trip chain from the first trip to trip j (trip j not included) are taking their maximum trip time, denote Eijγ and E˜ijγ as the binary variables that are equal to 1 if vehicle charging can be carried out on arc i,j∈A when trip i takes its nominal and maximum trip time, respectively; denote Uijγ and U˜ijγ as the integer variables that specify the number of time units an e-bus spends charging on arc i,j∈A when trip i takes its nominal and maximum trip time, respectively. Denote Rijγ and R˜ijγ as the auxiliary variables defined to linearize the quadratic terms in the constraints.In the robust model, the robustness of a trip chain is ensured by linearizing of the DP equations (25) and (26) and adding them into the problem constraints. The following constraints ((27)–(30)) are the linearized forms of the DP equations (25) and (26).(i)
Time consistency constraints:(27)Zjγ−Ziγ−ti−sij+M1−∑k∈KXijk≥0,i,j∈Ap,γ=0,1,2,…,Γt,(28)Zjγ−Zi,γ−1−ti−t^i−sij+M1−∑k∈KXijk≥0,i,j∈Ap,γ=1,2,…,Γt.Constraints (27) and (28) are the robust constraints corresponding to constraint (7). It is developed based on the linearization of equation (25), which means the time consistency of two successive trips should be satisfied for ∀γ=0,1,2,…,Γt.(ii)
Battery SoC consistency constraints:(29)Yiγ−ei+Wijγ−eij+M1−∑k∈KXijk≥Yjγ,i,j∈A,γ=0,1,2,…,Γt,(30)Yi,γ−1−ei−e^i+W˜ijγ−eij+M1−∑k∈KXijk≥Yjγ,i,j∈A,γ=1,2,…,Γt.Constraints (29) and (30) are the robust constraints corresponding to constraint (8). It is developed based on the linearization of equation (26), which means the battery SoC consistency of two successive trips should be satisfied for ∀γ=0,1,2,…,Γt.Other constraints are defined corresponding to the deterministic models with the variables substituted by the ones defined in the robust model. The R-MD-EVSP model is formulated as follows:(31)min∑k∈K∑i,j∈AcijXijk+∑γ=0Γt∑i,j∈AgijWijγ+∑γ=1Γt∑i,j∈AgijW˜ijγ,s.t. constraints (2)–(5) and constraints (27)–(30).(32)ai≤Ziγ≤bi,i∈V,γ=0,1,2,…,Γt,(33)Eijγ−∑k∈KXijk≤0,i,j∈A,γ=0,1,2,…,Γt,(34)E˜ijγ−∑k∈KXijk≤0,i,j∈A,γ=1,2,…,Γt,(35)−aj−Ziγ−ti−sij−tp−M1−Eijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(36)−aj−Zi,γ−1−ti−t^i−sij−tp−M1−E˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(37)Wijγ+Rijγrs−aj−ti−sij−tprsEijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(38)W˜ijγ+R˜ijγrs−aj−ti−t^i−sij−tprsE˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(39)Wijγ+Yiγ−ei−umax≤0,i,j∈A,γ=0,1,2,…,Γt,(40)W˜ijγ+Yi,γ−1−ei−e^i−umax≤0,i,j∈A,γ=1,2,…,Γt,(41)Wij−rstsUijγ=0,i,j∈A,γ=0,1,2,…,Γt,(42)Wij−rstsU˜ijγ=0,i,j∈A,γ=0,1,2,…,Γt,(43)Rijγ−biEijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(44)aiEijγ−Rijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(45)Rijγ−Ziγ+ai1−Eijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(46)Ziγ−bi1−Eijγ−Rijγ≤0,i,j∈Ap,γ=0,1,2,…,Γt,(47)R˜ijγ−biE˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(48)aiE˜ijγ−R˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(49)R˜ijγ−Zi,γ−1+ai1−E˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(50)Zi,γ−1−bi1−E˜ijγ−R˜ijγ≤0,i,j∈A,γ=1,2,…,Γt,(51)uiγmin≤Yiγ≤umax,i∈V,γ=0,1,2,…,Γt,(52)Wijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(53)W˜ijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(54)Xijk∈0,1,i,j∈A,k∈K,(55)Eijγ∈0,1,i,j∈A,γ=0,1,2,…,Γt,(56)Rijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(57)Uijγ∈Z,i,j∈A,γ=0,1,2,…,Γt,(58)E˜ijγ∈0,1,i,j∈A,γ=0,1,2,…,Γt,(59)R˜ijγ≥0,i,j∈A,γ=0,1,2,…,Γt,(60)U˜ijγ∈Z,i,j∈A,γ=0,1,2,…,Γt.Objective (31) consists of the vehicle usage, empty travel, and the sum of charging cost in all possible cases regarding the realizations of the trip times. The worst-case charging cost cannot be expressed explicitly; therefore, it is not included in the objective but calculated after the optimal solution is obtained. Constraint (32) limits the departure time of each trip to be within the departure time window. Constraints (33) and (34) require that charging on arc i,j is performed only if trip j is succeeded by trip i. Constraints (35) and (36) restrict that charging on arc i,j can be performed if time is available. Constraints (37) and (38) define the upper bound of Wij on arc i,j∈A according to the available time. Originally, constraints (37) and (38) are formulated as follows:(61)Wijγ+Ziγ−aj−ti−sij−tprsEijγ≤0,i,j∈A,γ=0,1,2,…,Γt,(62)W˜ijγ+Zi,γ−1−aj−ti−t^i−sij−tprsE˜ijγ≤0,i,j∈A,γ=1,2,…,Γt.To linearize quadratic termsZiγEijγ and Zi,γ−1E˜ijγ, auxiliary variables Rijγ and R˜ijγ are introduced where Rijγ=ZiγEijγ is satisfied by constraints (43)–(46) and R˜ijγ=Zi,γ−1E˜ijγ is satisfied by constraints (47)–(50).Constraints (39) and (40) require the battery SoC after charging not to exceed umax. Constraints (41) and (42) ensure that the charging time can only be integer times of the minimum charging time unit. Constraint (49) keeps the vehicle battery SoC within the safety range. Constraints (52)–(60) define the domains of the variables.


## 4. Branch-and-Price Method

In this section, we introduce a BP algorithm to solve the R-MD-EVSP based on the set partitioning formulation. LetΩ be the set of feasible trip chains. Each trip chain p∈Ω is defined by cp, the total operation cost, and api, i∈T, a binary parameter that is equal to 1 if trip i is carried out in p. Denote θp, p∈Ω, as a binary variable that takes value 1 if trip chain p is selected to be part of the solution. The set partitioning formulation of the R-MD-EVSP is as follows:(63)min∑p∈Ωcpθp,(64)s.t.∑k∈Kapiθp=1,i∈T,(65)θp∈0,1,p∈Ω.Objective (63) minimizes the total operation cost of the schedule. Constraint (64) ensures that each trip is carried out once by a vehicle. As set Ω includes a large number of feasible trip chains, it is impractical to solve the model directly. Therefore, we developed a column generation (CG) algorithm to solve the linear relaxation of models (63)–(65) called master problem (MP). The CG algorithm starts by solving the linear relaxation of models (63)–(65) with an initial set of feasible trip chains Ω¯, called the restricted master problem (RMP). In each iteration, a pricing problem is solved to generate columns with negative reduced cost to be added to the RMP. Let πi, i∈T, be the dual variables associated with constraint (64). Let c¯p be the reduced cost of trip chain p∈Ω with respect to πi, i∈T, i.e., c¯p=cp−∑i∈Tapiπi. The pricing problem is defined as follows:(66)min∑p∈Ωc¯pθp.A feasible trip chainp∈Ω corresponds to an o-d path in G. The problem aims at finding a feasible trip chain with the minimum reduced cost. We modify the cost of each arc i,j∈A by a modified cost c¯ij=cij−πi where πi=0 for i∈o,d. The reduced cost of an o-d path p equals the sum of c¯ij, i,j∈p. The problem is a resource-constrained elementary shortest path problem (RCESPP), and we solve it by developing a label setting algorithm [35]. The CG alternates between the optimization of the RMP and the pricing problems until no more columns with negative reduced cost are generated, implying that the MP has been solved to optimality.


### 4.1. Robustness Checking of a Trip Chain

In solving the pricing problem, we need to generate robusto-d paths in graph G. Based on the DP equations (25) and (26) defined in Section 3.3, the robustness of path P can be checked by the following procedure.The Robustness Checking Procedure. On an arc vi,vj of path P, for each γ≤Γt, determine the values of Wijγ and W˜ijγ by equations (67) and (68). Then, the values of Zjγ and Yjγ are determined by equations (25) and (26), respectively. Path P is robustness feasible if Ziγ∈ai,bi and Yiγ∈uiγmin,umax for vi∈P and γ≤Γt.(67)Wijγ=max0,minumax,aj−Ziγ−ti−tij−tprs,(68)W˜ijγ=max0,minumax,aj−Zi,γ−1−ti−t^i−tij−tprs.


### 4.2. Label Setting for the Pricing Problem

In this section, we describe the label setting algorithm for generating robusto-d paths based on the robustness checking approach proposed in Section 4.1. In solving the RCESPP on graph G, the labels are used to represent the partial paths starting from the start node and the resources accumulated along the path. New labels are generated by extending the existing labels on graph G following the resource extension functions (REFs). New labels are checked for resource feasibility and infeasible labels are discarded. At each node, dominance rule is then applied to eliminate the labels which cannot lead to a path with a better reduced cost than the dominate labels.DenoteLi=Ci,Riγt,Riγe, γ≤Γt, as a label representing a partial path from the start node o to a node i∈V with three kinds of resource defined as follows: Ci is the reduced cost of the path, Riγt is the earliest departure time of node i, and Riγe is the lowest battery SoC of the vehicle upon the start of trip i, when γ trips from node o to node i are attaining their maximum trip time.In the VRP, when the customer time windows are wide, two customers can be served in different orders so that the graph is cyclic. As such, in solving the pricing problem, at a given nodei in the graph, an elementary resource is defined for each customer to record the number of times the customer has been visited on the partial path o, 1, 2, …, i [31]. In the R-MD-EVSP, because the departure time window of the trips is small, any two trips cannot be served in different orders and graph G is acyclic. Therefore, in solving the pricing problem, the elementary resource is not needed.Based on the recursive equations (25) and (26) and the feasibility checking approach introduced in Section 4.1, the forward extension of a label Li to a label Lj along an arc i,j∈A is performed by REFs (69)–(71). As the feasibility checking approach always ensures Rjγe∈ujγmin,umax, the resulting label Lj is robust feasible if and only if Rjγt∈aj,bj, γ=0,1,2,…,Γt.(69)Cj=Ci+cij'−πi,(70)Rjγe=minumax,Riγe−ei+Wijγ−eij,Ri,γ−1e−ei−e^i+W˜ijγ−eij,(71)Rjγt=maxai,Riγt+ti+Wijγrs+tij,Riγ−1t+ti+t^i+W˜ijγrs+tij.VariablesWijγ and W˜ijγ are consistent with those defined in Section 3.3. cij' is the estimated cost on arc i,j∈A consisting of the cost of vehicle usage, empty travel, trip start delay, and charging cost. cij' is computed as follows. Let Pj=o,1,2,…,i,j be a partial path from node o to node j. Denote Hj as the maximum total charging amount required for Pj when the trip times varies within uncertainty set U. cij'=cij+gijHj−Hi where cij is the vehicle usage and empty travel cost assigned to arc i,j; gijHj−Hi is the increased charging cost.In the label setting algorithm, dominance rule plays a critical role in reducing the number of extensions required as the labeling process goes on. In Proposition1, we introduce a dominance rule for our label setting algorithm to improve the computational efficiency.Proposition 1.
LetLi1 and Li2 be the two labels associated with the paths ending at node i. Li1 dominates Li2 if (i) Ci1≤Ci2; (ii) Riγe1≥Riγe2, γ=0,1,2,…,Γt; and (iii) Riγt1≤Riγt2, γ=0,1,2,…,Γt.


### 4.3. Heuristic Pricing

The label setting algorithm can be time-consuming, especially in solving large-scale instances. As we know, it is not necessary to solve the pricing problem to its optimal and return a column with the most negative reduced cost at each iteration of the CG. The computational efficiency of the GC procedure can be improved by only finding suboptimal solutions in the pricing problem at the expense of an increase in the number of iterations. As such, we adopt heuristic decisions before invoking the exact pricing algorithm. Instead of storing all the nondominant columns in the label setting, we maintain only a prespecified number of them to improve the algorithm efficiency. When a new nondominant label is added to the list and the total number of labels has reached the limit, the label with the largest objective value in the list will be discarded. In this way, fewer labels are extended at each node, speeding up the process of negative path detection. Note that this heuristic may lead to more CG iterations and the number of labels maintained should be carefully set through preliminary experiment.


### 4.4. Branching

An initial solution is formulated by assigning each e-bus with one trip. After each CG procedure ends, to obtain the integer solution, branches are created in the branch-and-bound tree by branching on the total flow of the arcsxiji,j∈A. If not all the xiji,j∈A are integers, we choose to branch on xij that has the highest fractional value, creating two nodes, one with xij=1 and the other with xij=0. Otherwise, an integer solution is obtained. A depth-first search strategy is used to explore the branch-and-bound tree. In order to obtain integer solutions more efficiently, after each CG procedure ends, xij values in the RMP that are larger than 0.99 are set to 1, while xij values that are less than 0.01 are set to 0. The preliminary experiments show that this setting works well on speeding up the computational time and has minor impact on the solution quality.


### 4.5. Obtaining Optimal Charging Plan

After the optimal schedule is obtained by the BP algorithm, we further generate the optimal charging plan for each individual trip chain, considering the time-of-use tariff. The problem can be defined as a robust fixed-route electric vehicle refueling problem (R-FRVRP) which aims to determine the optimal charging time and amount for an e-bus on a given trip chain.


## 4.1. Robustness Checking of a Trip Chain

In solving the pricing problem, we need to generate robusto-d paths in graph G. Based on the DP equations (25) and (26) defined in Section 3.3, the robustness of path P can be checked by the following procedure.The Robustness Checking Procedure. On an arc vi,vj of path P, for each γ≤Γt, determine the values of Wijγ and W˜ijγ by equations (67) and (68). Then, the values of Zjγ and Yjγ are determined by equations (25) and (26), respectively. Path P is robustness feasible if Ziγ∈ai,bi and Yiγ∈uiγmin,umax for vi∈P and γ≤Γt.(67)Wijγ=max0,minumax,aj−Ziγ−ti−tij−tprs,(68)W˜ijγ=max0,minumax,aj−Zi,γ−1−ti−t^i−tij−tprs.


## 4.2. Label Setting for the Pricing Problem

In this section, we describe the label setting algorithm for generating robusto-d paths based on the robustness checking approach proposed in Section 4.1. In solving the RCESPP on graph G, the labels are used to represent the partial paths starting from the start node and the resources accumulated along the path. New labels are generated by extending the existing labels on graph G following the resource extension functions (REFs). New labels are checked for resource feasibility and infeasible labels are discarded. At each node, dominance rule is then applied to eliminate the labels which cannot lead to a path with a better reduced cost than the dominate labels.DenoteLi=Ci,Riγt,Riγe, γ≤Γt, as a label representing a partial path from the start node o to a node i∈V with three kinds of resource defined as follows: Ci is the reduced cost of the path, Riγt is the earliest departure time of node i, and Riγe is the lowest battery SoC of the vehicle upon the start of trip i, when γ trips from node o to node i are attaining their maximum trip time.In the VRP, when the customer time windows are wide, two customers can be served in different orders so that the graph is cyclic. As such, in solving the pricing problem, at a given nodei in the graph, an elementary resource is defined for each customer to record the number of times the customer has been visited on the partial path o, 1, 2, …, i [31]. In the R-MD-EVSP, because the departure time window of the trips is small, any two trips cannot be served in different orders and graph G is acyclic. Therefore, in solving the pricing problem, the elementary resource is not needed.Based on the recursive equations (25) and (26) and the feasibility checking approach introduced in Section 4.1, the forward extension of a label Li to a label Lj along an arc i,j∈A is performed by REFs (69)–(71). As the feasibility checking approach always ensures Rjγe∈ujγmin,umax, the resulting label Lj is robust feasible if and only if Rjγt∈aj,bj, γ=0,1,2,…,Γt.(69)Cj=Ci+cij'−πi,(70)Rjγe=minumax,Riγe−ei+Wijγ−eij,Ri,γ−1e−ei−e^i+W˜ijγ−eij,(71)Rjγt=maxai,Riγt+ti+Wijγrs+tij,Riγ−1t+ti+t^i+W˜ijγrs+tij.VariablesWijγ and W˜ijγ are consistent with those defined in Section 3.3. cij' is the estimated cost on arc i,j∈A consisting of the cost of vehicle usage, empty travel, trip start delay, and charging cost. cij' is computed as follows. Let Pj=o,1,2,…,i,j be a partial path from node o to node j. Denote Hj as the maximum total charging amount required for Pj when the trip times varies within uncertainty set U. cij'=cij+gijHj−Hi where cij is the vehicle usage and empty travel cost assigned to arc i,j; gijHj−Hi is the increased charging cost.In the label setting algorithm, dominance rule plays a critical role in reducing the number of extensions required as the labeling process goes on. In Proposition1, we introduce a dominance rule for our label setting algorithm to improve the computational efficiency.Proposition 1.
LetLi1 and Li2 be the two labels associated with the paths ending at node i. Li1 dominates Li2 if (i) Ci1≤Ci2; (ii) Riγe1≥Riγe2, γ=0,1,2,…,Γt; and (iii) Riγt1≤Riγt2, γ=0,1,2,…,Γt.


## 4.3. Heuristic Pricing

The label setting algorithm can be time-consuming, especially in solving large-scale instances. As we know, it is not necessary to solve the pricing problem to its optimal and return a column with the most negative reduced cost at each iteration of the CG. The computational efficiency of the GC procedure can be improved by only finding suboptimal solutions in the pricing problem at the expense of an increase in the number of iterations. As such, we adopt heuristic decisions before invoking the exact pricing algorithm. Instead of storing all the nondominant columns in the label setting, we maintain only a prespecified number of them to improve the algorithm efficiency. When a new nondominant label is added to the list and the total number of labels has reached the limit, the label with the largest objective value in the list will be discarded. In this way, fewer labels are extended at each node, speeding up the process of negative path detection. Note that this heuristic may lead to more CG iterations and the number of labels maintained should be carefully set through preliminary experiment.


## 4.4. Branching

An initial solution is formulated by assigning each e-bus with one trip. After each CG procedure ends, to obtain the integer solution, branches are created in the branch-and-bound tree by branching on the total flow of the arcsxiji,j∈A. If not all the xiji,j∈A are integers, we choose to branch on xij that has the highest fractional value, creating two nodes, one with xij=1 and the other with xij=0. Otherwise, an integer solution is obtained. A depth-first search strategy is used to explore the branch-and-bound tree. In order to obtain integer solutions more efficiently, after each CG procedure ends, xij values in the RMP that are larger than 0.99 are set to 1, while xij values that are less than 0.01 are set to 0. The preliminary experiments show that this setting works well on speeding up the computational time and has minor impact on the solution quality.


## 4.5. Obtaining Optimal Charging Plan

After the optimal schedule is obtained by the BP algorithm, we further generate the optimal charging plan for each individual trip chain, considering the time-of-use tariff. The problem can be defined as a robust fixed-route electric vehicle refueling problem (R-FRVRP) which aims to determine the optimal charging time and amount for an e-bus on a given trip chain.


## 5. Numerical Experiments

We conducted numerical experiments based on the case of bus routes in Shenzhen. We first conducted experiments on single-route scheduling and then on multiroute scheduling where e-buses are allowed to carry out trips on different routes. In the following cases, plug-in DC chargers are established in the depots, the majority of which have the maximum charging power of 100 kW. BYD e-buses with a battery capacity of 260 kWh, length of 10.5 m, weight of 18000 kg, and maximum driving distance of 220 km operate on the routes. The characteristics of the bus routes and charging facilities are presented in Table2. According to the historical running data, the charging rate is set to be 1.6 kWh per minute, and the battery energy consumption rate is estimated to be 0.6 kWh per minute. As defined in Section 3.1, the energy consumption of a trip depends on the trip time following the equation δie=αδitru. Here we set α=1. The charging preparation time tp and minimum charging time ts are set as 2 min and 5 min, respectively. Parameter α takes the value of 0.5. The daily usage cost of an e-bus cf is set as 1000 CNY. Figure 3 shows the time-of-use tariff in Shenzhen, which is used to set the values of gij, i,j∈A.Table 2
Characteristics of the bus routes in the single and multiroute scheduling cases.
RouteLength (km)Number of stopsOperation periodHeadway (min)Ave. trip time (min)Number of tripsEnd depot IEnd depot IIM13337.5446 : 10–22 : 305–121102761 (16)2 (8)E739.7157 : 00–22 : 0010–6070843 (20)4 (5)E1175.5186 : 00–21 : 306–1012019435 (10)E1441.9227 : 30–22 : 0015–501107036 (10)E2257.1196 : 30–20 : 00601203437 (5)4217.3236 : 30–22 : 3010–15451328 (10)9 (8)4317.6256 : 30–22 : 3010–1560130398131.6446 : 30–21 : 305–158020438Note. End depots are numbered from 1 to 9 as follows: 1, Changlingdong depot; 2, Shekou bus depot; 3, Shenzhen North Railway Station; 4, Longgangyiwu bus depot; 5, Nanaoxinda bus depot; 6, Tantou West; 7, Hangzi bus depot; 8, Moon Bay; 9, Window of the World. The values in the brackets next to the depot index are the number of charging piles established at the depots.Figure 3
Time-of-use tariff in Shenzhen.The MIP model was solved by the standard optimization solver Cplex 12.6, and the BP framework was coded in Java using Cplex 12.6 to solve the RMP. All the experiments were run on a PC with Windows 10, Intel Core i5-8250U, 1.80 GHz, and 8 GB RAM.


### 5.1. Single-Route Scheduling

We first analyze the performance of the MIP model and BP algorithm based on Route M133 in Shenzhen. A display of bus Route M133 is shown in Figure4.Figure 4
A display of bus Route M133 in Shenzhen.In generating a timetable, estimation of the scheduled trip times is critical. We adopted the method suggested in [27] to generate homogeneous running time (HRT) periods, the trip time within each of which follows the same distribution. We generate the HRT periods for Route M133 based on the trip time data of 100 weekdays from August 1 to December 31 in 2017. Figure 5 presents the distribution of the trip time of Route M133 during different periods of the day. Within each HRT period, the scheduled trip time is set according to a rule-of-thumb: the average trip time plus the standard deviation of the trip time. Figure 6 shows the scheduled trip time and maximum deviation from the scheduled trip time within each running time period. We generate the timetable including 276 trips. Except for the original timetable, we also create timetables with the number of scheduled trips ranging from 64 to 154 by changing the headway requirements. These timetables account for the small-scale instances.Figure 5
Trip time distribution of Route M133 in a day based on the historical running data.Figure 6
HRT periods of Route M133 generated based on the historical running data. (a) Outbound. (b) Inbound.
(a)(b)


#### 5.1.1. Results

The performance of the MIP model and BP algorithm with uncertainty budgetΓt=0,1,2,3 is presented in Table 3. The deterministic solution corresponds to the case with Γt=0. The table shows the objective of the MIP model obtained by Cplex and BP algorithm (Obj), number of vehicles used (#V) in the solution obtained by the BP, optimal gap of the MIP and BP, and the computational time. When the MIP gap is higher than 90%, we report the root relaxation solution obtained by Cplex in the bracket. If the root relaxation is not solved by Cplex within 3600 s, an em dash (“—”) is used.Table 3
Computational results of the MIP model and BP on single-route instances.
Number of tripsΓtObj#VGap (%)Run time (s)MIPBPBPMIPBPCplexBP62016916.516912.8165.4036000118107.316400.01611.64036001219876.816800.01619.50036001321763.317000.01726.480.123600196022900.722676.7218.300.1436008125664.725504.52218.180.2136005230395.027279.22230.910.2136006335098.333331.92240.170.2636008112026294.325711.12312.530.42360010130508.726714.72424.610.6036007236842.230202.22437.570.5836008351072.533858.52454.970.71360010154038082.336578.73313.350.20360024169344.337041.83352.410.223600342(33000)42180.033—0.68—563(33000)45606.033—0.84—132760(49000)56443.349—0.82—1291(49000)67309.149—0.86—1692—82386.749—1.02—1913—99134.449—0.92—113The results indicate that Cplex is not able to solve the model to near optimal in a short time. Some instances are not solved by Cplex in 3600 s; however, the BP algorithm is able to generate high-quality solutions with small BP gaps within a reasonable computational time. The objective obtained by the BP algorithm is consistently lower than that obtained by the MIP model. For a given instance, the total operational cost, including the cost of vehicle usage and charging, becomes higher with the increase of the uncertainty budgetΓt. This is due to the more vehicles put into use and a larger amount of planned daytime charging amount. When the level of protection Γt increases, the planned daytime charging amount becomes higher if the number of vehicles used stays the same; the daytime charging amount can be reduced when the number of vehicles used increases.


#### 5.1.2. Impact of the Trip Start Buffer Time

We investigated the impact of trip start buffer time on the optimal schedule. Each case is named as “number of trips—budget value.” The trip start delay timetd varies among 0, 3, 5, and 8 minutes. td=0 means that trip start delay is not permitted. As shown in Figure 7, the number of e-buses used is reduced with the increase of td in cases “96-0,” “112-0,” and “276-0.” The results indicate that allowing for trip start buffer time can reduce the number of e-buses put into operation in some cases. Trip start buffer time can increase time flexibility of the schedule and saves the operational cost in some scenarios. With the increase of the uncertainty budget, the planned charging amount increases, protecting the schedule against the variation of energy consumption.Figure 7
The influence of departure delay time on the number of vehicles used under different values of uncertainty budget.


### 5.2. Multiroute Scheduling

We carried out numerical experiments on two cases with multiple bus routes and depots. The multiroute scheduling is applied in the cases. The characteristics and layouts of the transit system in the case are shown in Table2 and Figure 8. Case I includes five depots and four express bus routes: Route Nos. E7, E11, E14, and E22, that all have one end stop at the Shenzhen North Railway Station, formulating a hub-and-spoke topology. The timetable of the four routes includes a total number of 382 trips. Case II includes three depots and three bus routes: Route Nos. 42, 42, and 81, that share end stops with each other, formulating a circle topology. All the depots are equipped with charging facilities. The timetable of case I and case II includes a total number of 382 and 466 trips.Figure 8
A display of the route layout of (a) case I and (b) case II.
(a)(b)


#### 5.2.1. Results

We obtained the optimal schedule for cases I and II by the BP algorithm. The computational results are shown in Table4. The table shows the objective of the MIP model obtained by the BP algorithm (Obj), number of vehicles used (#V), total charging amount, computational time, and BP gap with uncertainty budget Γt=0,1,2,3.Table 4
Computational results of the BP on multiroute instances.
CaseΓtObj.#VCharging amount (kWh)Run time (s)BP gap (%)Case I081638.37311141.33921.02192313.47616313.43980.962114236.07617522.46761.103132200.57618247.85901.32Case II069279.0646895.662481.89176556.26411178.489791.848261286540.86413533.282610.12398607.36413444.163501.06


### 5.3. Schedule Robustness

We investigated the schedule robustness through Monte Carlo simulation based on the randomly generated trip time data. The simulation is performed by generatingNrand=10000 random realizations of trip time from normal distribution Ntia,sdi for all trip nodes i∈T, where tia and sdi are the mean and standard deviation of the trip time. Each random realization represents one day’s e-bus schedule.Table5 presents the statistics on the deterministic (Γt=0) and robust (Γt=1,2,3) solutions for single-route instances. The number of vehicles used (#V), infeasibility rate of the trip chain (IR-TC), schedule (IR-S), and price of robustness (PoR) are reported. Denote the total number of infeasible trip chains as NtcIF. If the schedule contains infeasible trip chains, the schedule is regarded as infeasible. Denote the total number of infeasible schedules as NsIF. IR-TC is calculated as NtcIF/Nrand∗Nveh where Nveh is the number of e-buses used; IR-S is calculated as NsIF/Nrand. PoR is the percentage of increase in the operational cost of the robust solutions compared with that of the deterministic solution.Table 5
Robustness of the single-route schedules.
Number of tripsΓtInfeasibility rate (%)PoR (%)Trip chainSchedule9600.9419.27—10.000.009.5220.000.0015.4130.000.0030.7711207.2087.31—10.000.053.7620.000.0014.8730.000.0024.0615405.6592.08—10.010.261.2520.020.7013.2830.031.0419.7927609.8199.86—10.000.0016.1420.000.0031.4930.000.0043.06About 0.94% to 9.81% of the simulated e-bus trip chains experience delays using the deterministic solutions. WithΓt=1, the infeasibility rate drops to almost zero while the total operational cost increases by 1.25% to 16.14%. We calculated the penalty of the trips chains incurred by the violation of the trip start time window and battery SoC safety range. Figure 9 shows the changes of the objective and penalty per 100 trip chains under different values of Γt. With the increase of Γt, the objective increase while the penalty decreases. This result coincides with that in Table 3, indicating that the schedules with Γt=0 and Γt=1 are Pareto optimal solutions.Figure 9
Objective and penalty under different values ofΓt in single-route cases.Table6 presents the robustness evaluation results of the multiroute scheduling cases. When Γt=0, about 4.09% and 3.12% of the simulated trip chains experience delay in cases I and II, respectively. When Γt=1, the infeasibility rate drops to zero while the objective increases by 10.49% and 8.48% in cases I and II, respectively. We also calculated the penalty of the trips chains incurred by the violation of the trip start time window and battery SoC safety range. Figure 10 shows the changes of the objective and penalty per 100 trip chains under different values of Γt. With the increase of Γt, the objective increases while the penalty decreases. The results indicate that the schedules with Γt=0 and Γt=1 are Pareto optimal solutions.Table 6
Robustness of the multiroute schedules.
CaseΓtInfeasibility rate (%)PoR (%)Trip chainScheduleCase I04.0998.84—10.000.0010.4920.000.0019.1930.000.0013.59Case II03.1295.36—10.000.008.4820.000.0011.5430.000.0012.24Figure 10
Objective and penalty under different values ofΓt in multiroute cases. (a) Case I. (b) Case II.
(a)(b)Based on the above evaluation results, we can make the following conclusions:(i)
The infeasibility rate of the deterministic solution is significantly higher than that of the robust solutions in the simulation.(ii)
The schedule robustness can be greatly improved at the expense of an increase in the total operational cost.(iii)
Trade-off between the infeasibility rate and PoR is needed to choose a proper value of uncertainty budget when generating robust schedules.


## 5.1. Single-Route Scheduling

We first analyze the performance of the MIP model and BP algorithm based on Route M133 in Shenzhen. A display of bus Route M133 is shown in Figure4.Figure 4
A display of bus Route M133 in Shenzhen.In generating a timetable, estimation of the scheduled trip times is critical. We adopted the method suggested in [27] to generate homogeneous running time (HRT) periods, the trip time within each of which follows the same distribution. We generate the HRT periods for Route M133 based on the trip time data of 100 weekdays from August 1 to December 31 in 2017. Figure 5 presents the distribution of the trip time of Route M133 during different periods of the day. Within each HRT period, the scheduled trip time is set according to a rule-of-thumb: the average trip time plus the standard deviation of the trip time. Figure 6 shows the scheduled trip time and maximum deviation from the scheduled trip time within each running time period. We generate the timetable including 276 trips. Except for the original timetable, we also create timetables with the number of scheduled trips ranging from 64 to 154 by changing the headway requirements. These timetables account for the small-scale instances.Figure 5
Trip time distribution of Route M133 in a day based on the historical running data.Figure 6
HRT periods of Route M133 generated based on the historical running data. (a) Outbound. (b) Inbound.
(a)(b)


### 5.1.1. Results

The performance of the MIP model and BP algorithm with uncertainty budgetΓt=0,1,2,3 is presented in Table 3. The deterministic solution corresponds to the case with Γt=0. The table shows the objective of the MIP model obtained by Cplex and BP algorithm (Obj), number of vehicles used (#V) in the solution obtained by the BP, optimal gap of the MIP and BP, and the computational time. When the MIP gap is higher than 90%, we report the root relaxation solution obtained by Cplex in the bracket. If the root relaxation is not solved by Cplex within 3600 s, an em dash (“—”) is used.Table 3
Computational results of the MIP model and BP on single-route instances.
Number of tripsΓtObj#VGap (%)Run time (s)MIPBPBPMIPBPCplexBP62016916.516912.8165.4036000118107.316400.01611.64036001219876.816800.01619.50036001321763.317000.01726.480.123600196022900.722676.7218.300.1436008125664.725504.52218.180.2136005230395.027279.22230.910.2136006335098.333331.92240.170.2636008112026294.325711.12312.530.42360010130508.726714.72424.610.6036007236842.230202.22437.570.5836008351072.533858.52454.970.71360010154038082.336578.73313.350.20360024169344.337041.83352.410.223600342(33000)42180.033—0.68—563(33000)45606.033—0.84—132760(49000)56443.349—0.82—1291(49000)67309.149—0.86—1692—82386.749—1.02—1913—99134.449—0.92—113The results indicate that Cplex is not able to solve the model to near optimal in a short time. Some instances are not solved by Cplex in 3600 s; however, the BP algorithm is able to generate high-quality solutions with small BP gaps within a reasonable computational time. The objective obtained by the BP algorithm is consistently lower than that obtained by the MIP model. For a given instance, the total operational cost, including the cost of vehicle usage and charging, becomes higher with the increase of the uncertainty budgetΓt. This is due to the more vehicles put into use and a larger amount of planned daytime charging amount. When the level of protection Γt increases, the planned daytime charging amount becomes higher if the number of vehicles used stays the same; the daytime charging amount can be reduced when the number of vehicles used increases.


### 5.1.2. Impact of the Trip Start Buffer Time

We investigated the impact of trip start buffer time on the optimal schedule. Each case is named as “number of trips—budget value.” The trip start delay timetd varies among 0, 3, 5, and 8 minutes. td=0 means that trip start delay is not permitted. As shown in Figure 7, the number of e-buses used is reduced with the increase of td in cases “96-0,” “112-0,” and “276-0.” The results indicate that allowing for trip start buffer time can reduce the number of e-buses put into operation in some cases. Trip start buffer time can increase time flexibility of the schedule and saves the operational cost in some scenarios. With the increase of the uncertainty budget, the planned charging amount increases, protecting the schedule against the variation of energy consumption.Figure 7
The influence of departure delay time on the number of vehicles used under different values of uncertainty budget.


## 5.1.1. Results

The performance of the MIP model and BP algorithm with uncertainty budgetΓt=0,1,2,3 is presented in Table 3. The deterministic solution corresponds to the case with Γt=0. The table shows the objective of the MIP model obtained by Cplex and BP algorithm (Obj), number of vehicles used (#V) in the solution obtained by the BP, optimal gap of the MIP and BP, and the computational time. When the MIP gap is higher than 90%, we report the root relaxation solution obtained by Cplex in the bracket. If the root relaxation is not solved by Cplex within 3600 s, an em dash (“—”) is used.Table 3
Computational results of the MIP model and BP on single-route instances.
Number of tripsΓtObj#VGap (%)Run time (s)MIPBPBPMIPBPCplexBP62016916.516912.8165.4036000118107.316400.01611.64036001219876.816800.01619.50036001321763.317000.01726.480.123600196022900.722676.7218.300.1436008125664.725504.52218.180.2136005230395.027279.22230.910.2136006335098.333331.92240.170.2636008112026294.325711.12312.530.42360010130508.726714.72424.610.6036007236842.230202.22437.570.5836008351072.533858.52454.970.71360010154038082.336578.73313.350.20360024169344.337041.83352.410.223600342(33000)42180.033—0.68—563(33000)45606.033—0.84—132760(49000)56443.349—0.82—1291(49000)67309.149—0.86—1692—82386.749—1.02—1913—99134.449—0.92—113The results indicate that Cplex is not able to solve the model to near optimal in a short time. Some instances are not solved by Cplex in 3600 s; however, the BP algorithm is able to generate high-quality solutions with small BP gaps within a reasonable computational time. The objective obtained by the BP algorithm is consistently lower than that obtained by the MIP model. For a given instance, the total operational cost, including the cost of vehicle usage and charging, becomes higher with the increase of the uncertainty budgetΓt. This is due to the more vehicles put into use and a larger amount of planned daytime charging amount. When the level of protection Γt increases, the planned daytime charging amount becomes higher if the number of vehicles used stays the same; the daytime charging amount can be reduced when the number of vehicles used increases.


## 5.1.2. Impact of the Trip Start Buffer Time

We investigated the impact of trip start buffer time on the optimal schedule. Each case is named as “number of trips—budget value.” The trip start delay timetd varies among 0, 3, 5, and 8 minutes. td=0 means that trip start delay is not permitted. As shown in Figure 7, the number of e-buses used is reduced with the increase of td in cases “96-0,” “112-0,” and “276-0.” The results indicate that allowing for trip start buffer time can reduce the number of e-buses put into operation in some cases. Trip start buffer time can increase time flexibility of the schedule and saves the operational cost in some scenarios. With the increase of the uncertainty budget, the planned charging amount increases, protecting the schedule against the variation of energy consumption.Figure 7
The influence of departure delay time on the number of vehicles used under different values of uncertainty budget.


## 5.2. Multiroute Scheduling

We carried out numerical experiments on two cases with multiple bus routes and depots. The multiroute scheduling is applied in the cases. The characteristics and layouts of the transit system in the case are shown in Table2 and Figure 8. Case I includes five depots and four express bus routes: Route Nos. E7, E11, E14, and E22, that all have one end stop at the Shenzhen North Railway Station, formulating a hub-and-spoke topology. The timetable of the four routes includes a total number of 382 trips. Case II includes three depots and three bus routes: Route Nos. 42, 42, and 81, that share end stops with each other, formulating a circle topology. All the depots are equipped with charging facilities. The timetable of case I and case II includes a total number of 382 and 466 trips.Figure 8
A display of the route layout of (a) case I and (b) case II.
(a)(b)


### 5.2.1. Results

We obtained the optimal schedule for cases I and II by the BP algorithm. The computational results are shown in Table4. The table shows the objective of the MIP model obtained by the BP algorithm (Obj), number of vehicles used (#V), total charging amount, computational time, and BP gap with uncertainty budget Γt=0,1,2,3.Table 4
Computational results of the BP on multiroute instances.
CaseΓtObj.#VCharging amount (kWh)Run time (s)BP gap (%)Case I081638.37311141.33921.02192313.47616313.43980.962114236.07617522.46761.103132200.57618247.85901.32Case II069279.0646895.662481.89176556.26411178.489791.848261286540.86413533.282610.12398607.36413444.163501.06


## 5.2.1. Results

We obtained the optimal schedule for cases I and II by the BP algorithm. The computational results are shown in Table4. The table shows the objective of the MIP model obtained by the BP algorithm (Obj), number of vehicles used (#V), total charging amount, computational time, and BP gap with uncertainty budget Γt=0,1,2,3.Table 4
Computational results of the BP on multiroute instances.
CaseΓtObj.#VCharging amount (kWh)Run time (s)BP gap (%)Case I081638.37311141.33921.02192313.47616313.43980.962114236.07617522.46761.103132200.57618247.85901.32Case II069279.0646895.662481.89176556.26411178.489791.848261286540.86413533.282610.12398607.36413444.163501.06


## 5.3. Schedule Robustness

We investigated the schedule robustness through Monte Carlo simulation based on the randomly generated trip time data. The simulation is performed by generatingNrand=10000 random realizations of trip time from normal distribution Ntia,sdi for all trip nodes i∈T, where tia and sdi are the mean and standard deviation of the trip time. Each random realization represents one day’s e-bus schedule.Table5 presents the statistics on the deterministic (Γt=0) and robust (Γt=1,2,3) solutions for single-route instances. The number of vehicles used (#V), infeasibility rate of the trip chain (IR-TC), schedule (IR-S), and price of robustness (PoR) are reported. Denote the total number of infeasible trip chains as NtcIF. If the schedule contains infeasible trip chains, the schedule is regarded as infeasible. Denote the total number of infeasible schedules as NsIF. IR-TC is calculated as NtcIF/Nrand∗Nveh where Nveh is the number of e-buses used; IR-S is calculated as NsIF/Nrand. PoR is the percentage of increase in the operational cost of the robust solutions compared with that of the deterministic solution.Table 5
Robustness of the single-route schedules.
Number of tripsΓtInfeasibility rate (%)PoR (%)Trip chainSchedule9600.9419.27—10.000.009.5220.000.0015.4130.000.0030.7711207.2087.31—10.000.053.7620.000.0014.8730.000.0024.0615405.6592.08—10.010.261.2520.020.7013.2830.031.0419.7927609.8199.86—10.000.0016.1420.000.0031.4930.000.0043.06About 0.94% to 9.81% of the simulated e-bus trip chains experience delays using the deterministic solutions. WithΓt=1, the infeasibility rate drops to almost zero while the total operational cost increases by 1.25% to 16.14%. We calculated the penalty of the trips chains incurred by the violation of the trip start time window and battery SoC safety range. Figure 9 shows the changes of the objective and penalty per 100 trip chains under different values of Γt. With the increase of Γt, the objective increase while the penalty decreases. This result coincides with that in Table 3, indicating that the schedules with Γt=0 and Γt=1 are Pareto optimal solutions.Figure 9
Objective and penalty under different values ofΓt in single-route cases.Table6 presents the robustness evaluation results of the multiroute scheduling cases. When Γt=0, about 4.09% and 3.12% of the simulated trip chains experience delay in cases I and II, respectively. When Γt=1, the infeasibility rate drops to zero while the objective increases by 10.49% and 8.48% in cases I and II, respectively. We also calculated the penalty of the trips chains incurred by the violation of the trip start time window and battery SoC safety range. Figure 10 shows the changes of the objective and penalty per 100 trip chains under different values of Γt. With the increase of Γt, the objective increases while the penalty decreases. The results indicate that the schedules with Γt=0 and Γt=1 are Pareto optimal solutions.Table 6
Robustness of the multiroute schedules.
CaseΓtInfeasibility rate (%)PoR (%)Trip chainScheduleCase I04.0998.84—10.000.0010.4920.000.0019.1930.000.0013.59Case II03.1295.36—10.000.008.4820.000.0011.5430.000.0012.24Figure 10
Objective and penalty under different values ofΓt in multiroute cases. (a) Case I. (b) Case II.
(a)(b)Based on the above evaluation results, we can make the following conclusions:(i)
The infeasibility rate of the deterministic solution is significantly higher than that of the robust solutions in the simulation.(ii)
The schedule robustness can be greatly improved at the expense of an increase in the total operational cost.(iii)
Trade-off between the infeasibility rate and PoR is needed to choose a proper value of uncertainty budget when generating robust schedules.


## 6. Conclusion

In this paper, an R-MD-EVSP model and BP algorithm are introduced as a novel approach for e-bus scheduling under travel time uncertainty. Partial charging and trip start time window are taken into consideration. The cardinality constrained uncertainty set was adopted to control the level of uncertainty. The proposed methodology was tested on the real-world transit operation cases in Shenzhen including single-route and multiroute cases. The size of the instances varies from 62 to 466 timetabled trips. The results showed that the BP algorithm is able to generate provably high-quality solutions efficiently for large-scale instances. With the increase of the uncertainty budget, the planned charging amount increases to protect the schedule against the variation of energy consumption when the number of e-buses used stays the same; the charging amount can be reduced when the number of vehicles used increases. Additionally, a sensitivity analysis of the trip start buffer time was carried out. The analysis showed that allowing for trip start buffer time can increase time flexibility of the schedule and save the operational cost in some scenarios. The robustness of the schedules was evaluated through Monte Carlo simulation based on the randomly generated trip time data. The results showed that the infeasibility rate of the deterministic solution is significantly higher than that of the robust solutions. Applying the robust schedules can effectively protect against trip time and energy consumption variation at the expense of an increase in the operational cost. Therefore, transit operators need to trade off between the infeasibility rate and PoR to choose a proper value of uncertainty budget when adopting robust schedules.In future research studies, more realistic battery charging and energy consumption process should to be considered in the robust e-bus scheduling model. Besides, heuristic methods are needed to solve the problem instances of larger scale. Furthermore, different kinds of uncertainty sets can be considered to control the trip time variation. Trip time and energy consumption can be described by independent uncertainty sets, accounting for more general e-bus operation cases.



---

*Source: 1191443-2021-12-20.xml*