#### 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.

#### 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. Section 3 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 Table 1.

Let be an acyclic direct graph, where is the set of nodes and is the set of arcs. Each timetabled trip is represented by a trip node in the graph. Denote as the set of trip nodes. Each trip node is associated with a trip start time window , where and are the scheduled start time and latest start time, respectively. The e-bus serving a trip should wait until the scheduled start time if it arrives earlier at the trip start depot. Each trip node is associated with two uncertain parameters: the trip time and trip energy consumption . is realized in the interval where is the nominal trip time; and are the upper and lower deviations from . is realized in the interval where is the nominal value; and are the upper and lower deviations from . Since the worst case will always be achieved at the right-hand side of the interval, we only consider the realization of in and in correspondingly.

Denote as the set of depots. Each depot is created with two nodes in the graph: operation start node and operation end node , indicating that the e-bus begins/ends the operation from/at depot , respectively. For nodes and , , we assign a time window , where is the earliest operation start time and is the latest operation end time. This means that the e-buses can begin and end their operation at any time within time window .

The arc set includes three kinds of arcs: (i) the pull-out arc connects a depot node , , and a trip node , representing that a vehicle begins operation from depot to carry out the first trip ; (ii) the pull-in arc connects a trip node and an end node , , representing that a vehicle finishes its last trip and ends the operation, returning to depot ; and (iii) the trip connection arc connects two trip nodes, representing that two trips are carried out consecutively. Each arc is associated with a nonservice travel time and energy consumption . To ensure that the trip start time window is respected, trip nodes and , , are connected only if the time compatible condition is satisfied. The cost of an arc includes the empty travel cost and the vehicle usage cost for the pull-out arcs. Denote , , as the unit charging cost in the time interval between the end of trip node and the start of trip node ; is calculated as the weighted average value of the time-of-use tariff in the time interval.

Figure 1 gives an example of the graph for an instance with two depots and four trips. Each depot is created with a start node and an end node , . Each trip is represented by a trip node associated with a departure time window and a scheduled trip time, .

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 of (kWh/min). Denote as the unit charging time and as the charging preparation time. The e-buses should charge integer times of .

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. Figure 2 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 , the energy consumption deviation is proportional to the trip time deviation , i.e., where is an augment constant and 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 varies within interval ; as such, we only need to give the values of and as the model inputs. Therefore, any empirical method that can give estimation on the values of and works.

**(a)**

**(b)**

**(c)**

**(d)**

We adopt a cardinality constrained uncertainty set as discussed in [24] to control the level of trip time uncertainty. is defined as where the cumulative uncertainty of the random variable is bounded by the budget . 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 set ; (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 when the trip time varies within the budgeted set . 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

Let be the binary variables that take value 1 if a vehicle housed at terminal traverses arc . At each trip node , denote as the trip start time and as the lowest battery SoC of the vehicle at the trip start time. The lower bound of is denoted as . At each operation start (end) node (), , denote as the operation start (end) time and as the lowest battery SoC of the vehicle at the operation start (end) time. On each arc , , denote as the amount of energy to be charged after completing trip before starting trip . According to our assumption, the charging place is at the ending depot of trip node . are the auxiliary variables defined corresponding to for constraint linearization. Denote as the binary variables that are equal to 1 if the e-buses can be charged on arc . Denote as the integer variables that specify the number of time units an e-bus spends charging on arc . The deterministic problem is formulated as follows (MD-EVSP):

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 is performed only if trip is succeeded by trip . Constraint (10) restricts that charging on arc can be performed only if time is available. Constraint (11) defines the upper bound of on arc according to the time available. Originally, constraint (11) is formulated as follows:

To linearize quadratic terms , auxiliary variables are introduced where is satisfied by constraints (14)–(17). Constraint (12) requires the battery SoC after charging not to exceed . 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 set . In graph , a trip chain is represented by a path starting from an operation start node to an operation end node , visiting a sequence of trip nodes . In solving the RVRP, Agra et al. [29] showed that only the travel time vector needs to be considered in the problem formulation to ensure the robustness feasibility. is the set that contains all the extreme points of set . Because of the structure of set , 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.

Denote as the earliest possible start time of node and as the lowest battery SoC of the vehicle upon the start of ; let and be the charging amount on arc , i.e., after the end of trip node and before the start of trip node , when trips from to are taking their maximum trip time, among which trip node takes the nominal trip time and maximum trip time, respectively. Denote as the lower bound of : if , ; otherwise, .

and are defined by the recursive functions (25) and (26), respectively. A trip chain is robust if and , for and . In fact, it suffices to check and for .

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

Denote as the optimal total charging amount of a trip chain. is the maximum total charging amount required when the trip times varies within uncertainty set which is calculated as follows: let be the maximum number of trips of path that can attain their maximum trip time. Choose the 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, is given by where is the total energy consumption of . The charging policy for an e-bus carrying out trip chain 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 until the total charging amount reaches .

##### 3.4. Robust Problem Formulation

We introduce an MIP formulation for the R-MD-EVSP under cardinality uncertainty set by incorporating the DP equations for trip chain robustness checking into the deterministic model.

At each node , given trips have taken their maximum trip time, , variables , , , and parameter 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 and as the binary variables that are equal to 1 if vehicle charging can be carried out on arc when trip *i* takes its nominal and maximum trip time, respectively; denote and as the integer variables that specify the number of time units an e-bus spends charging on arc when trip takes its nominal and maximum trip time, respectively. Denote and 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: 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 .(ii)Battery SoC consistency constraints:

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 .

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:s.t. constraints (2)–(5) and constraints (27)–(30).

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 is performed only if trip is succeeded by trip . Constraints (35) and (36) restrict that charging on arc can be performed if time is available. Constraints (37) and (38) define the upper bound of on arc according to the available time. Originally, constraints (37) and (38) are formulated as follows:

To linearize quadratic terms and , auxiliary variables and are introduced where is satisfied by constraints (43)–(46) and is satisfied by constraints (47)–(50).

Constraints (39) and (40) require the battery SoC after charging not to exceed . 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 is defined by , the total operation cost, and , , a binary parameter that is equal to 1 if trip is carried out in . Denote , , as a binary variable that takes value 1 if trip chain is selected to be part of the solution. The set partitioning formulation of the R-MD-EVSP is as follows:

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 , , be the dual variables associated with constraint (64). Let be the reduced cost of trip chain with respect to , , i.e., . The pricing problem is defined as follows:

A feasible trip chain corresponds to an *o-d* path in . The problem aims at finding a feasible trip chain with the minimum reduced cost. We modify the cost of each arc by a modified cost where for . The reduced cost of an *o-d* path equals the sum of , . 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 robust *o-d* paths in graph . Based on the DP equations (25) and (26) defined in Section 3.3, the robustness of path can be checked by the following procedure.

*The Robustness Checking Procedure*. On an arc of path , for each , determine the values of and by equations (67) and (68). Then, the values of and are determined by equations (25) and (26), respectively. Path is robustness feasible if and for and .

##### 4.2. Label Setting for the Pricing Problem

In this section, we describe the label setting algorithm for generating robust *o-d* paths based on the robustness checking approach proposed in Section 4.1. In solving the RCESPP on graph , 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 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.

Denote , , as a label representing a partial path from the start node to a node with three kinds of resource defined as follows: is the reduced cost of the path, is the earliest departure time of node , and is the lowest battery SoC of the vehicle upon the start of trip , when trips from node to node 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 node 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 , 1, 2, …, [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 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 to a label along an arc is performed by REFs (69)–(71). As the feasibility checking approach always ensures , the resulting label is robust feasible if and only if , .

Variables and are consistent with those defined in Section 3.3. is the estimated cost on arc consisting of the cost of vehicle usage, empty travel, trip start delay, and charging cost. is computed as follows. Let be a partial path from node to node . Denote as the maximum total charging amount required for when the trip times varies within uncertainty set . where is the vehicle usage and empty travel cost assigned to arc ; 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 Proposition 1, we introduce a dominance rule for our label setting algorithm to improve the computational efficiency.

Proposition 1. *Let and be the two labels associated with the paths ending at node . dominates if (i) ; (ii) , ; and (iii) , .*

##### 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 arcs . If not all the are integers, we choose to branch on that has the highest fractional value, creating two nodes, one with and the other with . 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, values in the RMP that are larger than 0.99 are set to 1, while 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 Table 2. 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 . Here we set . The charging preparation time and minimum charging time are set as 2 min and 5 min, respectively. Parameter takes the value of 0.5. The daily usage cost of an e-bus is set as 1000 CNY. Figure 3 shows the time-of-use tariff in Shenzhen, which is used to set the values of , .

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 Figure 4.

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.

**(a)**

**(b)**

###### 5.1.1. Results

The performance of the MIP model and BP algorithm with uncertainty budget is presented in Table 3. The deterministic solution corresponds to the case with . 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.

The 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 . This is due to the more vehicles put into use and a larger amount of planned daytime charging amount. When the level of protection 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 time varies among 0, 3, 5, and 8 minutes. 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 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.

##### 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 Table 2 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.

**(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 Table 4. 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 .

##### 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 generating random realizations of trip time from normal distribution for all trip nodes , where and are the mean and standard deviation of the trip time. Each random realization represents one day’s e-bus schedule.

Table 5 presents the statistics on the deterministic () and robust () 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 . If the schedule contains infeasible trip chains, the schedule is regarded as infeasible. Denote the total number of infeasible schedules as . IR-TC is calculated as where is the number of e-buses used; IR-S is calculated as . PoR is the percentage of increase in the operational cost of the robust solutions compared with that of the deterministic solution.

About 0.94% to 9.81% of the simulated e-bus trip chains experience delays using the deterministic solutions. With , 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 . With the increase of , the objective increase while the penalty decreases. This result coincides with that in Table 3, indicating that the schedules with and are Pareto optimal solutions.

Table 6 presents the robustness evaluation results of the multiroute scheduling cases. When , about 4.09% and 3.12% of the simulated trip chains experience delay in cases I and II, respectively. When , 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 . With the increase of , the objective increases while the penalty decreases. The results indicate that the schedules with and are Pareto optimal solutions.

**(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.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported in part by the Basic Research Program of Shenzhen Science and Technology Innovation Committee (JCYJ20180307123910003), Scientific Research Start-Up Funds of Tsinghua Shenzhen International Graduate School (QD2021007N), and National Natural Science Foundation of China (61673233).