Integrated optimization of electric bus scheduling and charging planning incorporating flexible charging and timetable shifting strategies

In a battery electric bus (BEB) network, buses are scheduled to perform timetabled trips while satisfying time, energy consumption, charging, and operational constraints. Increasing research efforts have been dedicated to the integrated optimization of multiple planning tasks to reduce system costs. At a high integration level, this study determines the BEB scheduling and charging planning with flexible charging and timetable shifting strategies. We first formulate an integrated arc-based model to minimize the total costs considering the power grid pressure cost and subsequently reformulate it into a two-stage model, for which we develop an effective solution method. The first stage minimizes the total operational costs including the fleet, charging, and battery degradation costs based on the column generation technique, and the second stage minimizes the peak power demand through two timetable shifting strategies. It is found through numerical experiments that the proposed integrated optimization model and solution method can achieve significant improvement in the utilization rate and reductions in the fleet size, operational costs, and peak power demand compared to the two baseline models.


Introduction
It has been reported that the transportation sector accounts for about one-quarter of the global carbon dioxide emissions due to the wide use of carbon-fueled vehicles (IEA, 2019), which is one of the main causes of climate change. This proportion is even higher in developed economies, for example, one-third in the European Union (IEA, 2022). Therefore, it is necessary to find substitute mobility options to reduce emissions and secure a more sustainable environment. Electric mobility has proven to be instrumental in decarbonizing the transportation sector (Qu et al., 2022;Razeghi and Samuelsen, 2016). Amongst, electric public bus systems are considered an important pillar that offers affordable and environmentally friendly mobility. Local governments in many countries implement a series of policies, such as financial incentives, customer subsidies, and taxation on petroleum, to promote the deployment of battery electric buses (BEBs) (Nie et al., 2016). Despite the advantages of BEBs, they still face challenges in large-scale adoption (Ji et al., 2022;Zhang et al., 2021a). The high capital investment, including the costs of embedded battery packs and auxiliary charging facilities, makes BEBs less appealing than traditional diesel buses. The battery capacity of the BEBs gradually decreases during the charging and discharging cycles, resulting in the degradation of battery performance (Schoch et al., 2018;Zhang et al., 2019). Furthermore, due to the complex planning processes for a BEB system, inefficient charging strategies and fleet schedules are usually adopted by BEB operators (Yao et al., 2020;Zhang et al., 2021b).

Literature review
In general, the planning tasks of bus systems include line designing, timetabling, bus scheduling, crew rostering, and charging planning (if the buses are BEBs). In the existing studies of bus systems, timetabling and bus scheduling are two primary research subjects. Timetabling is to determine the departure and arrival times of each trip in a bus network constrained by the given frequency and level of service of each line. Bus scheduling is to assign buses to serve the timetabled trips on the condition that each trip is covered by one bus. Because of the limited travel range of BEBs, the timetable needs reconstructions and more idle time for the charging process, where charging planning makes decisions on when, how long, and how much electricity the bus is charged. As reviewed below, these planning tasks have often been handled separately due to the high model complexity.
Typically, a timetable is determined based on the line characteristics and the passenger demands between origins and destinations. With an emphasis on passenger services, the common objective of timetabling is to minimize passenger travel or transfer times (Wu et al., 2016), while ignoring the subsequent planning tasks. For the BEB scheduling problem, route duration, route distance, and the duration for replenishing energy at the charging points should be constrained due to the limitation of the driving range per charge. Some studies assigned the BEBs to meet the condition that each timetabled trip is covered once, aiming to yield the least operational costs. For example, given fixed charging duration, Tang et al. (2019) proposed robust scheduling strategies for BEBs under stochastic traffic conditions, and Rinaldi et al. (2020) focused on the optimal scheduling of a mixed BEB fleet. These studies did not show any charging flexibility since the charging duration or volume is fixed. Specifically, each charging operation starts immediately upon arrival at the charging station and has a determined charging duration or volume according to the "first-in-first-out" principle. Such an inflexible charging strategy cannot fully coordinate the demand of BEBs and the availability of charging resources. Assuming predetermined fleet schedules, some studies were dedicated to the orderly and flexible charging strategy. For example, Wang et al. (2017) arranged the charging demand of BEBs to ensure recharging without any delays. He et al. (2020) optimized charging schedules for a fast-charging system and determined when to charge a BEB. The obtained charging strategies are called "semi-flexible charging", which indicates that the buses are charged to a pre-given electricity level while optimizing the sequence of charging operations. Unlike their studies, Abdelwahed et al. (2020) optimized the charging schedules and determined where, when, and how long each BEB should be charged while considering the time-dependent electricity prices in time-of-use (TOU) plans. This charging strategy is called "fullflexible charging", as all the charging-related components need to be determined endogenously, such as the charging start time, duration, and volume. However, the studies mentioned above adopted separated optimization for single planning tasks, which tends to cause inefficient BEB planning outcomes (Schmid and Ehmke, 2015). Therefore, it is desirable to optimize these planning tasks integratively to exploit the system's capacity to the greatest extent and maximize its productivity and efficiency (Ceder, 2001;Liu and Ceder, 2018).
Integrated optimization has recently been a hotspot topic in the study of timetabling, bus scheduling, and charging planning. Commonly, two interdependent planning tasks are addressed simultaneously, such as integrating timetabling with bus scheduling or integrating bus scheduling with charging planning. The integrated optimization of timetabling and bus scheduling is only targeted at conventional buses. For instance, Ibarra-Rojas et al. (2014) combined two integer linear programming models of timetabling and vehicle scheduling to form a bi-objective integrated model. One objective is to maximize the number of passengers benefiting from well-coordinated transfers, and the other is to minimize the operational costs related to fleet management. Schmid and Ehmke (2015) aimed to minimize operational costs by balancing consecutive departures on a line according to predefined departure time intervals. In addition to the cost-efficiency objective, they also considered the quality of a timetable from a passenger's point of view. Liu and Ceder (2018) developed a bi-level integer programming model to attain the optimality of timetable synchronization with vehicle scheduling. The integration of three or more planning tasks is even more desirable but seems infeasible for most practical problems (Teng et al., 2020).
Increasing research has integrated charging planning for BEBs into the aforementioned integrated optimization models developed for non-BEBs. For instance, considering the limited driving range and charging requirement constraints, Liu and Ceder (2020) proposed two formulations to minimize the total number of BEBs and chargers. One is based on the deficit function theory, and the other is an equivalent bi-objective integer programming model. Under travel range uncertainty,  addressed the multi-depot vehicle location-routing-scheduling problem with two-stage stochastic programming. Some realistic operational conditions, such as battery degradation, nonlinear charging profile, and peak power demand, were incorporated in different studies. To be specific, the capacity and power of batteries in BEBs gradually fade due to the chemical and mechanical processes, and this degradation mechanism always occurs owing to cyclic charging and discharging (Barré et al., 2013;Pelletier et al., 2018); the nonlinear charging profile implies varying charging efficiency responding to different charging conditions (Liu et al., 2022;Zhang et al., 2021a); peak power demand has a massive impact on the existing power system (Chen et al., 2018;Lopes et al., 2011;Vicini et al., 2012), which is triggered by a onetime occurrence of charging demand and usually measured as the highest average power over a time interval (e.g., varying from 15 to 60 min) (He et al., 2020;Hledik, 2014;Qin et al., 2016). Pelletier et al. (2018) developed a comprehensive mathematical model that incorporates a realistic charging process, time-dependent energy costs, battery degradation, grid restrictions, and facility-related demand charges. They pointed out that optimizing the charging schedule contributes to keeping the facility-related charging demand low at the expense of spreading out the charging activities throughout the day. Zhang et al. (2021) determined the BEB charging strategies to minimize the total operational costs of a transit system. These studies reported that the incorporation of these charging conditions enhances realism and accounts for feasible fleet size reductions and charging strategies. The integrated optimization problems of any two of the above planning tasks have been demonstrated to be NP-hard (Ibarra-Rojas et al., 2014;Rinaldi et al., 2020). The integrated optimization of bus scheduling, charging planning, and timetabling with the consideration of realistic charging conditions would be even more challenging.

Focus of this study
As stated above, to the best of our knowledge, no study has investigated the bus scheduling and charging planning with timetable shifting problem while considering the battery degradation, non-linear charging profile, and TOU plans. For comparative convenience, we summarize the detailed characteristics of the relevant studies on the operational planning of BEBs in Table 1. As seen, the literature is evolving to take into account realistic operational conditions such as battery degradation, non-linear charging profile, and TOU plans. However, three notable literature gaps can be identified. First, the integrated optimization models of bus scheduling and charging planning have not considered timetabling. As such, the existing studies fall short of investigating the influence of charging strategies and moderating the peak power demand by timetabling-related strategies, such as timetable shifting. Second, most existing studies on bus scheduling and charging planning are based on fixed or semi-flexible charging strategies to obtain the sequence of service trips. This treatment of charging simplifies the modeling complexity but cannot capture the real operational conditions. Abdelwahed et al. (2020) was the only one that considered full-flexible charging, but their model focused on this single planning task. Third, existing approaches for solving the scheduling problem focused on heuristics, branch-and-price, and optimization solver. Although their problems can be solved well by these methods, they cannot solve the proposed integrated optimization problem well. To fill these gaps, this study aims to coordinate bus scheduling, charging planning, and timetabling with full-flexible charging as well as incorporate realistic operational conditions in a BEB system. Therefore, the main contributions of this study are presented as follows. First of all, we investigate an integrated optimization model for timetable modifications and joint bus and charging scheduling. This study thus far has the highest model integrity in the literature. Second, with the consideration of full-flexible charging, battery degradation, nonlinear charging profile and TOU plan, the model aims to generate a joint fleet and charging schedule to minimize the total costs, considering the power grid pressure cost. These considerations can capture the real operational conditions. Third, we reformulate the arc-based model and develop a two-stage solution method. This treatment separates the integrated objective into the operational costs and the peak power demand, which simplify the complexity of the solution method.
The remainder of this paper is organized as follows. Section 2 presents the problem description. In Section 3, an optimization model is proposed to characterize the integrated timetabling, BEB scheduling, and charging planning problem. Section 4 discusses the twostage solution method based on an arc-based model reformulation. Numerical experiments are conducted in Section 5 to verify the effectiveness and computational efficiency of the proposed model and solution method. Finally, Section 6 concludes the main contributions and provides potential extensions for future work.

Problem description
In this section, we first describe the setting of the research problem and then introduce the timetable shifting strategies in detail.  Liu et al. (2021) N

Research problem
The research problem focuses on BEB scheduling and charging planning with timetable shifting. The studied transit network can be illustrated in Fig. 1, which is composed of a single depot where all buses are housed at night, one central charging station with a limited number of candidate chargers located in the depot, several transit lines served by BEBs, and bus stops. The central depot charging requires that a BEB returns to the depot to replenish their electricity during the idle time between two service trips. According to the practice of BEB operations, the BEBs run round trips along the lines, and the start and end locations of the bus line, i.e., the bus stops of departure and arrival, are called terminals. For example, Line1 is composed of a sequence of stops; once a daily operation begins, a BEB departs from the depot and then performs the round trip based on terminal A. Due to the limited space capacity, a depot is usually located in a different location from the terminals. Therefore, there is certain distance between the depot and the terminals. To ensure that BEBs have enough electricity to perform the tasks in chronological order, they are not allowed to return to the depot except for replenishing their electricity. This network structure with one central depot and multiple neighboring terminals applies broadly to small-sized and medium-sized cities (Perumal et al., 2021) and has been widely used in the BEB literature.
The planning horizon T is discretized into equal time intervals as T = {0,1,2,⋯,|T|}, where | • | gives the cardinality of a set. The sets of terminals and directed lines are denoted by S and L, respectively. Each line l ∈ L has its original timetable T l 0 , and the set of timetabled trips is represented as I l = {1, ..., |I l |} in ascending order by their departure times. Thus, we have I = ∪ l∈L I l and I l ∩ I l ′ = ∅, l, l ′ ∈ L, l ∕ = l ′ . Trip i is associated with a line l i ∈ L, terminal s i ∈ S, departure time d i ∈ T l 0 , travel time tt i , and electricity e i (in percentage) consumed by the trip. For consistency, the volume of electricity replenished at the charging station is also expressed in percentage. The objective is to minimize the total costs, considering the power grid pressure cost. The three key components, including modified timetable, fleet schedule, and the corresponding charging schedule are elaborated below.

Timetable shifting strategy
The limited number of chargers at the central depot often leads to a high peak power demand (Liu and Ceder, 2018;Lin et al., 2019). We adopt the timetable shifting strategy, which has been applied for timetable synchronization in public transport (Fonseca et al., 2018), to adjust the start charging time to moderate the peak power demand. Large deviations from the original timetable would not only perturb the regular timetabled services but also confuse passengers with accustomed trip plans. More importantly, they would cause massive adjustments to the following scheduling tasks. To coordinate with the fleet and charging schedules, the original timetable should only be shifted within an acceptable range. We consider two types of shifting strategies, including block shifting and differentiated shifting. Block shifting indicates that the timetabled trips of a line are shifted equally while differentiated shifting means that the departure time of each timetabled trip makes a small shift within a time window. Essentially, block shifting is a special case of differentiated shifting. In a mathematical way, a generic form is proposed to represent the shifting strategy. Specifically, let d i denote the original departure time of trip i and Δ i denote the shift of trip i. Note that Δ i may be positive or negative, where a positive value indicates the time is moved backward and a negative value indicates the time is moved forward. Shift Δ i should be constrainted by a time window, that is where d − i and d + i denote the earliest and latest departure times at the terminal, respectively. When all trips in the timetable of one line are shifted equally, it indicates a block shifting strategy. To distinguish it from differentiated shifting, Δ l is introduced to denote the block shift of line l, and d + l represent the pre-determined minimum and maximum time shifts of line l, respectively. After timetable shifting, the actual departure time of trip i is represented as d i + Δ l + Δ i ,∀i ∈ T l 0 ,l ∈ L. Note that the departures of the trips in the same line can never overtake each other, meaning that the sequence of departures remains. Similar to the formulations in Fonseca et al.(2018), the shift window for trip i can be created based on the original timetables as if i refers to the first or the last trip of a line, we set d − i = d + i . For illustration purposes, Fig. 2 shows the two timetable shifting strategies given the original timetable. Fig. 2(a) presents the block shift Δ l < 0 for line l, indicating that all the timetabled trips of line l are moved forward |Δ l |; Fig. 2(b) presents the shift time windows , indicating that the departure times can be shifted within the time windows.

Model formulation
This section presents the formulations of the research problem. First, we show the primary assumptions and notations of this study. Second, the network structure and components used to model the problem are presented. Next, we discuss the constraints concerning the flow conservation of BEBs, time window, energy consumption and replenishment, timetable shifting, and charging station capacity. Finally, we consider three parts of the operational costs and the peak power demand in the objective function. The total operational costs include the fleet cost, the charging cost, and the cost incurred by battery degradation.

Notations and assumptions
The primary notations used in the formulations are listed in Table 2, and the assumptions used throughout this paper are summarized as follows.
Assumption 1. BEBs are homogeneous in terms of battery size and capacity. Each bus can only serve one line on an operation day, and the battery capacity is sufficient to cover at least one service trip.

Assumption 2.
There is only one charging station located in the depot where the physical distances between the chargers are neglected. A charging operation of a BEB can only occur at the charging station and is carried out with one of the chargers.

Assumption 3. A BEB can be charged only after finishing a round trip at the terminal.
Assumption 4. The routes of the buses for all the origin-destination pairs are known and fixed. Stochasticity in the road network is not considered, and the running times of BEBs between any two consecutive stops remain the same as those in the original timetable.
Assumption 5. A BEB with an initial level of battery electricity leaves the depot and begins its service, and the electricity of buses running on lines should always fall within the feasible range.
Assumptions 1 and 5 are consistent with the common practices of bus systems An, 2020). Assumptions 2, 3 and 4 are adopted to simplify the problem, which are widely adopted in the literature related to the planning and operation of BEBs (Pelletier et al., 2018;Chen et al., 2018;Liu et al., 2021;Uslu and Kaya, 2021;Zhang et al., 2021a). Although the fast-charging method occurs at terminals in many studies Wang et al., 2017), we make an assumption that the charging operations occur at the charging station since they can be transformed in the same formulations by the method of reduction (i.e., setting the distances between the depot and the terminals to 0). In addition, the consideration of time and energy consumption stochasticity is beyond the scope in this paper.

Network of fleet schedule and charging strategy
To address the research problem rigorously, we follow Li (2014)   also include deadhead trips without transporting passengers. The arc set A o = {(i, j)|i ∈ O o , j ∈ I} denotes the bus moving from the depot to the departure terminal of trip j when its daily operation starts, called a pull-out arc shown in Fig. 3.
called pull-in arc, denotes the bus moving from the arrival terminal of trip i to the depot when its daily operation ends. Arc set A I = {(i, j)|i, j ∈ I, i ∕ = j} denotes the bus moving from the arrival terminal of trip i to the departure terminal of trip j, called a service arc for connecting two trips. Suppose the remaining energy is not sufficient for the next trip. In that case, the bus needs to replenish its electricity at the charging station and then perform the next trip, which can be represented as {(i, c)→(c, j)|i, j ∈ I, i ∕ = j} and included in the service arcs. The trip chain indicates that a BEB starts from the depot, then performs several timetabled trips, replenishes its electricity between the idle time of trips if the battery electricity is insufficient for the next service trip, and finally ends at the depot. The fleet schedule refers to the generation of trip chains, subject to the constraint that each timetabled trip is only performed once by one bus. One trip chain consists of some nodes and arcs, and the trips are in ascending order regarding their departure times, satisfying energy consumption and charging time constraints. We illustrate the feasible fleet and charging schedules in Fig. 3(a). The time interval associated with each node represents the period of the corresponding nodes, such as the trip node with the scheduled departure and arrival times, the charging node with the start and end charging times, the state of charge (SOC) range before and after charging. This feasible fleet schedule R ⊆ R contains three trip chains, denoted as R = {r 1 , r 2 , r 3 }, namely, where c denotes charging. As for the charging schedule in Fig. 3

(a), it includes three charging operations expressed by a tuple
where c j is an indicator variable equal to 1 if a charging operation happens after trip j and 0 otherwise, ct j denotes the start charging time, et j denotes the end charging time, [soc o ,soc t ] indicates the SOC range before and after charging (%). This representation also implies the charging duration and volume of the charging operation. The original charging schedule requires 2 chargers; by timetable shifting, the charging demand can be replenished by 1 charger; the adjusted timetable and the optimized charging time are represented below nodes in Fig. 3(a). In addition, Fig. 3 indicates the power demand comparison between without and with timetable shifting for the illustrative example. Setting 15 min as a time interval, the peak power demand is 106.6 kW without shifting and 74.6 kW with shifting; the difference is 30% in this example, which indicates an advantage of the timetable shifting strategy. To guarantee the physical turnaround of the bus, we consider an average turnaround time ρ.
Related to the charging schedule, the charging volume needs to be determined, which is impacted by the charging duration and other factors such as the TOU plans and the nonlinear charging profile. Let c e t denote the electricity price in $/kWh associated with timestep t ∈ T. As for the nonlinear charging profile, we use a constant current-constant voltage (CC-CV) scheme shown in Fig. 4(a), commonly adopted by electric vehicles (Lam and Bauer, 2013;Schoch et al., 2018;Xu and Meng, 2019). During the constant current (CC) phase, the charging current is held constant, and the battery's terminal voltage increases until it reaches a certain maximum value, which is specified by the manufacturer. The constant voltage (CV) phase starts subsequently, and the terminal voltage must be maintained at the maximum value to avoid overcharging the battery. This scheme indicates that the variation of the battery's SOC cannot be assumed linear with the charging time (Pelletier et al., 2016). Similar to Montoya et al. (2017), we make a piecewise linear function to formulate the variation of SOC and the charging duration. Let f(χ) be the charging function when soc 0 = 0 and the battery is charged for χ time intervals. As shown in Fig. 4(b), soc 1 , soc 2 , and soc 3 represent the switch points of SOC in each piecewise part. Based on soc 0 , we have χ 0 = f − 1 (soc 0 ) and know in which piecewise part χ 0 +Δχ is located. To put it simply, soc 0 and Δχ are the input parameters of the piecewise function to obtain soc χ . Therefore, the final SOC can be calculated given the piecewise function of a charging profile, the SOC before visiting the charging station, and the charging duration.

Flow conservation constraints
A fleet schedule consists of several trip chains, and each trip chain performed by one bus is connected by the arcs from the pull-out arc, then to the service arcs, and finally to the pull-in arc as illustrated in Fig. 3(a). To capture the relationship, x ij is introduced as a binary decision variable, where x ij = 1 means trips i and j are successively served by one bus and x ij = 0 otherwise. To reduce the search space for the trip chain, infeasible or incompatible service arcs need to be removed from A I beforehand. For service trip arc (i,j), let d i and d j denote the scheduled departure time, respectively, where i, j ∈ T l 0 , l ∈ L. Considering the shifting strategy, the departure time of trips will be modified possibly and shifted in a window s j is greater than the latest departure time from the same terminal, service arc (i, j) should be eliminated. This condition indicates that trip j can never be performed on time even if trip i arrives at the earliest, that is, i +tt i denotes the earliest arrival time of trip i, d + j denotes the latest departure time of trip j, and ρ is the bus turnaround time.
In practice, each timetabled trip is exactly included once in the fleet schedule, as formulated in Eq. (2). To guarantee the continuity of a trip chain, flow conservation between the connecting arcs needs to be satisfied as Eq.

Time window constraints
As stated above, this study considers flexible charging with varying start charging time and duration. With minor shifts from the scheduled departure time, we use τ d i and τ a i as non-negative integer variables to denote the actual departure and arrival times of trip i respectively, which are limited by a shift time window of the original timetable. If the charging operation happens in the trip chain, the continuity should also be guaranteed. We define binary decision variable μ i as a charging indicator, where μ i = 1 if a charging operation happens after trip i and μ i = 0 otherwise. In addition, given the flexible charging and TOU plans, the charging schedule is related to two binary cumulative time decision variables. If the bus is charged at a timestep t * after trip i, a it = 1 at t ≥ t * and a it = 0 at other timesteps. Similarly, if the bus finishes charging at the timestep during t * +Δ t after trip i, b it = 1 at t ≥ t * +Δ t and b it = 0 at other timesteps. Fig. 5 illustrates the relationship between a it and b it .
For service arc (i,j), τ a i is affected by the timetable shifts and the turnaround time ρ, and thus should be equal to or earlier than τ d j , as formulated in Eq. (4). The start and end charging times are constrained by Eqs. (5a) and (5b) where M indicates a large positive number. The charging duration is bounded by the minimum and maximum charging durations . The minimum charging duration is determined due to economic efficiency, while the maximum charging duration is constrained by the battery capacity of the BEB. Let τ(i) be the charging duration immediately after trip i, which can be calculated by a it and b it in Eq. (6a). The first term on the right side of the equation represents the end charging timestep, and the second term represents the start charging timestep.
Furthermore, let τ min denote the minimum charging duration and the actual charging duration is constrained by Eq. (6b), indicating that τ(i) is positive if the charging operation happens after trip i, and τ(i) = 0 otherwise. Given the maximum charging duration |T|, the upper charging duration is limited by the latter constraint.
Finally, a it and b it delimit the flexible charging duration as Eqs. (7a)-(7c). The coupling constraints between the charging operation and the start and end charging times are shown in Eqs. (8a) and (8b).

Energy consumption and replenishment constraints
The energy consumption constraints are used to guarantee the continuity of the trip chain. If the remaining electricity of a BEB does not satisfy the next trip j, the bus needs to replenish energy at the charging station and then perform trip j. Eq. (9a) captures the discharging dynamics between trips i and j, where H i denotes the remaining SOC after completing trip i and H c i denotes the remaining SOC of the bus after trip i but before charging happens. Similarly, if a charging operation happens after trip i, the energy consumption and the replenishment should be considered. As described above, the charging electricity increases nonlinearly with the charging duration and the initial electricity. Let w(i) be the proportion of charging electricity, which is determined by τ(i) and H c i . The proportion of charging electricity is calculated as w − H c i . The energy consumption and replenishment are constrained by Eqs. (9b) and (9c), where e(s i , O c ) denotes the proportion of deadhead energy consumption from the terminal to the charging station, and e(O c , s j ) denotes the energy consumption from the charging station to the arrival terminal. Noted that if a charging operation immediately happens after trip i, Let H min and H max be the pre-defined minimum and maximum electricity and H init be the initial SOC of a BEB at the beginning of the daily operations. To mitigate the range anxiety of the driver, the remaining electricity after trip i should be constrained as Eq. (10a) following the natural process. Also, to avoid a BEB to visit the charging station frequently, the minimum level of electricity after charging is constrained by Eq. (10b), indicating that the sum of H c i and w(i) is non-negative if and only if the charging operation happens; H c i = 0, otherwise.
Finally, to ensure the integrity of the energy consumption of each bus when a daily operation starts or ends, the following constraints should be met. Specifically, Eq. (11a) indicates that each bus departs with battery electricity H init from the depot, and Eq. (11b) indicates that when the bus finishes its operation, the remaining SOC should be larger than the minimum value.

Timetable shifting constraints
Based on the original timetable, the timetable shifting strategies adjust slightly departure times of the trips for better synchronization. The block shift Δ l of timetable T l 0 for line l is constrained by Eq. (12a), and the adjusted departure time should fall within the bounds defined by the lower and upper time shifts, formulated as Eq. (12b). According to Section 2.2, the actual departure time after adjustment is formulated as Eq. (13).
In addition, the adjusted timetable should meet the headway constraints. The minimum headway is determined by the operator for safety consideration and the maximum headway is usually related to the service level. The headway between sequenced trips i and i +1 is constrained by Eq. (14), where h − i and h + i represent the minimum and maximum headways, respectively.

Capacity constraints
The number of chargers is limited by the available space at the charging location (Abdelwahed et al., 2020). The number of buses charged at any timestep should not exceed the number of chargers at the charging station. As shown in Fig. 5, a it − b it indicates whether the bus is charging at timestep t. μ i (a it − b it ) = 1 denotes that a charging operation occurs for trip i at timestep t. Each bus can serve a maximum of one trip at any timestep, and a trip can only be served by one bus. Thus, the constraint of charging station capacity is formulated as Eq. (15), where κ denotes the number of chargers at the depot.

Fleet cost
Fleet cost is dependent on the number of used BEBs. In general, it is related to the amortized cost, the depreciation cost, the cost of insurance and maintenance, etc. (Zhang et al., 2021a). In the BEB network, the number of buses corresponds to the number of trip chains in the fleet schedule. With the definition of the pull-out arc indicator x oi , the fleet cost associated with the number of buses can be calculated by where c b is the fixed cost of one BEB.

Charging cost
The charging cost is dependent on the charging time, charging volume, and unit electricity price. As formulated above, the charging volume is represented as To calculate the charging cost, w it is introduced as the charging electricity at timestep t after trip i. If a charging operation happens after trip i, the proportion of charging electricity is formulated as where Δt denotes the unit time interval, t indicates the charging timestep, and H i,t− 1 denotes the electricity level at timestep t − 1. Because of the implement of TOU plans, c e t denotes the unit electricity price at timestep t. And the charging cost can be calculated as Eq. (17), where E denotes the battery capacity of a BEB.

Battery degradation cost
The battery degradation cost is incurred by each cyclic discharging and charging process of the battery (Lam and Bauer, 2013). Let ∇ϑ = {(i, j)|i, j ∈ I} be the set of consecutively charging operations, which can be obtained from the charging schedule along the trip chain. Let ϑ ij be a consecutively charging operation indicator, where ϑ ij = 1 means two consecutively charging operations happen after trips i and j served by the same bus, and ϑ ij = 0 otherwise. According to Zhang et al. (2021), we show the discharging and charging of electricity between the consecutively charging operations in Fig. 6. If ϑ ij = 1, the SOC changes of one cyclic discharging and charging Fig. 6. Illustration of discharging and charging cycle. process are represented as where H c i and H c j denote the remaining SOC at the charging station before charging, H c i +w(i) and H c j +w(j) denote the SOC after charging, respectively. Let d(*) denote the charging or discharging degradation cost corresponding with the SOC changes, which is calculated by an empirical model (Lam and Bauer, 2013) developed to describe the aging of batteries over each charge and discharge cycle. We provide detailed calculations in the supplementary document. The total battery degradation cost is formulated as Eq. (18)

Peak power demand
The peak power demand reflects the stress of the power grid, which is a prescribed measurement over time intervals during the billing period (Perumal et al., 2021). Let p t denote the unit power of the charging station at timestep t, represented as p t = ∑ i∈I w it E, ∀t ∈ T. Let α denote the length (in the unit of h) of a time interval for calculating the peak power demand; thus, a day (24 h) is divided into 24/α demand periods. The set of demand periods is denoted by Q = {1,2,...,24/α}. For a demand period indexed by q ∈ Q, t 1 q and t 2 q denote the start and end times, respectively. The average power demand θ q of q is calculated as Eq. (19a) and the peak power demand θ is defined as Eq. (19b).
Taken together, the objective function includes costs related to the fleet, charging, battery degradation, and peak power demand.
Because of the different units, σ is introduced as a positive weight factor for θ by converting the pressure of the power grid to the economic consequence. The BEB scheduling and charging planning problem can be formulated as an integrated optimization model as Eq. (20). In fact, σ is a critical parameter indicating how the pressure of the power grid is integrated into the model. In a relevant study with a similar objective function, He et. al (2020) set the value of σ fixed as 12.765 $/kW. However, since it is hard to calibrate the value for applications, we do not set a fixed value of σ but develop an efficient method to solve it in the next section.  (2022). It should be pointed out that our model is an extension of Zhang et al. (2021) and the primary differences are that the full-flexible charging and timetable shifting strategies as well as the objective to reduce peak power demand are considered in our study. These features undoubtedly add up to model complexity. Given the high model complexity, an efficient solution method is needed to address a sizable BEB network.

Solution algorithm
To solve the above-formulated problem, we design a two-stage solution method, which is commonly adopted in the literature (Carosi et al., 2019;Huang et al., 2021;. The objective function of the arc-based model is decomposed to separate the operational costs and the peak power demand. This is done by assuming that for transit companies, minimizing the operational costs is the primary goal and reducing the peak power demand is the secondary goal.Section 4.1 describes the solution method of the first stage, minimizing the total operational costs based on the column generation (CG) technique. Section 4.2 describes the timetable shifting strategies for minimizing the peak power demand.

First-stage method based on CG
The problem of BEB scheduling and charging planning in the first stage is NP-hard as the number of variables and constraints in the formulation increases exponentially with the network size. We reconstruct a path-based formulation and use the CG technique to solve the first-stage problem. The CG technique has been demonstrated as effective in solving large-scale optimization problems (Steiner and Irnich, 2018;Xu and Meng, 2019;Liu et al., 2021). The CG starts with a restricted master problem (RMP) with a subset of variables and associated columns. After optimally solving the RMP, the dual solutions are used to find the non-basic variables with the minimum reduced costs implicitly, referred to as a pricing problem, which can be solved by a labeling method (Gao et al., 2022). If the reduced minimum cost is negative, the variable and column corresponding to the reduced minimum cost are added to the RMP, and the next iteration is triggered; otherwise, the optimal solution for this stage is claimed to be found in the current iteration. We use a larger timestep in the first-stage problem, while a smaller one in the second-stage problem ensures computational efficiency without lossing solution fidelity. The reasons are as follows. First, the size of the timestep affects the number of labels in the pricing problem, and as the step size increases, the labels will grow exponentially (Xu and Meng, 2019;Lin et al., 2020). Second, the purpose of the first stage is to obtain the fleet schedule and the charging volume, for which the step size has less effect on the objective value, as demonstrated in Section 5.1.

Path-based formulation
We define a binary variable λ r , which equals 1 if path r is selected, and 0 otherwise. The path-based model can be formulated as where C and C r denote the total operational costs of path r, respectively; incidence variable K r i equals 1 if trip i is covered by path r, and 0 otherwise; incidence variable U r t equals 1 if the BEB is charged at timestep t in path r, and 0 otherwise.
where path r is associated with one BEB; FC r denotes the cost of the BEB; CC r denotes the charging cost for the BEB; DC r denotes the battery degradation cost, and ∇ϑ r denotes the set of consecutively charging operations along path r.
For notational convenience, δ = {δ i |i ∈ I} and β = {β t |t ∈ T} are used as the dual variables associated with Eqs. (21b) and (21c). By relaxing Eq. (21d), we can use the CG technique to tackle the RMP and pricing problem iteratively. The CG technique starts with an initial feasible solution R to solve the RMP and determine the dual variables. Next, we construct and minimize the reduced cost Θ(λ r ) to solve the pricing problem. If Θ(λ r ) ≥ 0, no column can be added to R and the RMP has been solved with a proven optimal solution in this stage; otherwise, variable λ r obtained from the pricing problem is added to R for future iterations. The reduced cost function of variables λ r is expressed as As the CG technique starts with an initial solution, we introduce a simple heuristic method to construct a feasible path set, which is usually adopted by operators (Abdelwahed et al., 2020;Liu et al., 2021). In the beginning, the departure time τ d i is fixed to the original timetable T l 0 for all i ∈ I l ,l ∈ L, and the initial SOC of a BEB is equal to H init when departing from depot O. At the outset, a BEB serves the trip with the earliest departure time. After this trip is completed, we check if the BEB can still serve another trip satisfying the minimal battery level constraint. If so, we choose the trip with the earliest departure time from the remaining trip set; otherwise, we choose available timesteps for charging the bus until the battery electricity reaches to H init . If no trip can be served anymore, we direct the current trip to the depot, and the constructed trip chain constitutes an initial feasible RMP solution. This heuristic method is repeated until all trips in the trip set I are served. As a consequence, a feasible solution is generated.

Pricing problem
It can be seen from Eq. (22) that checking Θ(λ r ) ≤ 0 corresponds to the problem of finding the minimum-cost paths in the network G (V, A). The pricing problem aims to find a minimum-cost path from source node O o to destination node O d subject to time and electricity constraints. Essentially, the constraints of the pricing problem are similar to the original arc-based model. The pricing problem is to find a feasible trip chain with a negative reduced cost, while the original arc-based model is to find the optimal fleet schedules composed of several trip chains. The pricing problem can be formulated as where Eq. (23a) minimizes the reduced cost of path r; the pull-in and pull-out from the depot are constrained by Eqs. (23b) and (23c).
The pricing problem is an extension of the resource-constrained shortest path problem (Li, 2014) by allowing replenishment along a path, for which the label correcting method is suitable for solving the problem. To implement the labeling method, each node i ∈ I is assumed to be associated with multiple labels representing partial paths. The label extension at each node is due to the different charging strategies, including whether a charging operation happens after this node, the start charging times, and the charging duration. According to Liao (2016Liao ( , 2019, we initially construct a three-dimensional matrix of 2 × |T| × |T| units for each node, where each unit represents the potential label of this node. G is reconstructed as a space-time network with acyclic paths. Let i→j be a link of G(i, j ∈ I), and we code any label k at node i as l k where c k and w k are the cost and current SOC of the corresponding path, s k denotes the charging state with a value of 0 or 1, t → k and t k are the start charging time and charging duration, respectively. Also, p k , sp k , and cs k are used to record the current path, the SOC change at each node of the path, and the corresponding charging strategies of this partial path, respectively. All these labels at node i are grouped into a set denoted by L(i). Taking arc (i, j) for example to illustrate the label correcting process, suppose we have label k at node i represented by l k at node j is generated if the feasibility and dominance checks are performed for link i→j, where c u = c k + c j + c ij , w u = w k + H ij + w(i), c ij and H ij denote the arc cost and energy consumption respectively, s u = 1 or 0 is determined by the electricity feasibility check, the values of t → u and t u are determined by the time feasibility check. The dominance check is to discard the labels at node j that are dominated by any other labels in terms of the cost and SOC.
The CG-based method is typically incurred with the tailing-off effect, indicating that the method reaches a plateau and becomes slow to approach the optimum (Wang et al., 2018). To overcome this problem, we use the following stabilization technique (Addis et al., 2012;L. Zhang et al., 2021) to speed up the convergence. Let Φ = {δ, β} denote the vector of a dual solution of the RMP and Φ represent the update of Φ. Initially, we set Φ as vector 0, and a modified dual vector Φ can be computed by Φ = ε 1 • Φ +(1 − ε 1 )Φ in the iterative process, where ε 1 ∈ [0, 1] is a weighting parameter. Given ε 1 , the CG procedure is executed until no columns can be added to the RMP or the objective function reduction of RMP is no greater than a priori defined threshold ε 2 after a pre-defined number of iterations. Then, we update Φ as Φ and increase ε 1 gradually for a new iteration of the above process. The CG procedure terminates when ε 1 = 1 and no columns can be added.

Second-stage method of timetable shifting
Timetable shifting affects the charging schedule but has little impact on the trip sequences served by the buses. As described in Section 2.2, two timetable shifting strategies are applied. The second stage aims to minimize the peak power demand while controlling the cost increase of the first stage and modifications to the timetable. The outputs of this stage include the adjusted departure times and the start and end charging times. We formulate the model of this stage as Eqs. (4)-(6a), (12)-(15). Different time shifts may induce the same peak power demand, and modifications to the timetable are also presented in the objective function. The first term of Eq. (24) denotes the weighted peak power demand and the second term denotes the weighted value of time shifts, where ω 1 is a weight factor. Eq. (25) indicates that the total operational costs increase influenced by the shifting strategy should be within a certain range, where ω 2 is the controlling factor. In addition, it should be noted that the fleet and battery degradation costs do not change as the fleet schedule and the SOC changes are determined in the first stage. Therefore,FC and DC in Eq. (25) can be neglected. In other words, the timetable shifting strategy simply changes the charging cost CC.
In this stage, a smaller timestep size (e.g., 1 min) is adopted because a larger timestep makes no sense for the timetable modifications. Nevertheless, a small timestep size significantly increases the solution space, and it is difficult to solve directly by a solver (e.g., CPLEX, GUROBI). We design an iterative timetable shifting algorithm involving two shifting strategies, in which block shifting is first applied for all trips, followed by the differentiated shifting implemented many times for the subset trips with minor modifications. In other words, once the set of shifts in the blocking shifting is determined, the differentiated shifting strategy is implemented several times until the peak power demand does not decrease in an acceptable time; then the next iteration of block shifting begins. The timetable shifting algorithm involves a restricted model of Eq. (24), where variable Δ l is pre-determined and variable τ d i for the subset trips are fixed. With this variable fixed strategy, the restricted model could be directly solved by a solver.
As for the selection of the subset of timetabled trips, based on Fonseca et al. (2018), we designed a random selection strategy, denoted as Z(u), where u denotes the size of the selected trip set. We assume that any trip i ∈ I has an equal probability of being selected. We set the selected size u as 0.3|I|. In each iteration of the second stage, the selected trips are denoted as subset I ′ . For each i ∈ I ′ , the differentiated shifts Δ i are set as 0, i.e., the departure time τ d i is fixed as d i + Δ l . For unselected trip i ∈ I\I ′ , Δ i is the decision variable and need to be solved in the second stage. Note that the selecting strategy are applied multiple times until the objective value of the second stage remains the same comparing with the last selection. Obviously, u is related to the solution precision of the second stage and this value is associated with the solution efficiency.
The pseudo-code of the timetable shifting algorithm is described as follows. The input consists of a set of trips I, the original timetable T 0 , fleet schedule R, and charging strategy C 1 R with specified charging durations, two stopping criteria (stopCriterion1 and stopCriterion2), and a selection strategy Z(u). Both the two stopping criteria are based on the gaps between the solutions obtained in two adjacent iterations. Specifically, given Φ 1 and Φ 2 as the solution values of two adjacent iterations and gap threshold ε, if (Φ 2 − Φ 1 )/Φ 1 ≤ ε, the iterative process ends and otherwise continues. The algorithm starts with solving Eq. (24) without timetable modifications; thus, departure time τ d i is fixed to T 0 for all i ∈ I. An initial solution S 0 is generated composed of the initial timetable T 0 and charging strategy C 2 R including all start charging times. The iterative procedure is described in Lines 3-10, which runs until the stopCriterion2 is met.
The iteration of block shifting starts in Line 4 by random selecting Δ l satisfying Eq. (12a) and a new timetable T η 1 is generated, where η 1 records the number of iterations. The iteration of differentiated shifting starts in Line 6 by applying trip selection strategy Z(u) . A new solution is calculated in Line 7 by solving Eq. (24), with τ d i fixed to T η 1 for all i ∈ I ′ . The obtained solution is at least as good as the current best solution due to the stopping criteria. S * is constantly updated as the best-found solution and returned once the stopCriterion2 is met.
Iterative Timetable Shifting Algorithm Input: I, T 0 , R, C 1 R , stopCriterion1, stopCriterion2, Z Initialization: 1: solve Eq. (24) by a solver with Δ l = 0 for all l ∈ L, τ d i fixed to T 0 for all i ∈ I 2: While stopCriterion1 not satisfied do 4: T η 1 ← random selected Δ l for all ∀l ∈ L satisfying Eq. (12a) 5: While stopCriterion2 not satisfied do 6: I ′ ← selecting trips with strategy Z 7: (24) by a solver with τ d i fixed to T η 1 for all i ∈ I ′ 8: end while 9: The pseudo-code of the two-stage solution method is outlined below, and the flowchart is depicted in Fig. 7. First, the method starts by obtaining the initial solution with simple heuristic method. Subsequently, the two stages including CG technique and timetable shifting performs until the stopCriterion0 is met. The stopCriterion0 of the whole algorithm is set as (C ′ − C)/C < ε, where C ′ and C are the total operational costs of two adjacent iterations in the CG procedure, and ε is the relative gap threshold.
Remark 2. Note that the solutions obtained by the CG technique may not be optimal due to the lack of an optimality procedure (e.g., branch and bound) (Xu and Meng, 2019). However, this is acceptable because the result of the first stage is further optimized in the second stage through an iterative process. As the possible combinations of time shifts constitute a sizeable set, it is difficult to solve Eq. Obtain F1 = {R1, C R1 } and C1 using the CG procedure with input F0, presented in Section 4.1.1;

7:
Obtain F η 2 = {R η 2 , C R η2 } and C η 2 using the CG procedure with inputs R0, C R η2 − 1 and C 2 R η2 , presented in Section 4.1.1; 8: end while return F η 2 , C η 2 Remark 3. We present the complexity analyses for the proposed arc-based model and the two-stage solution method. In general, the complexity of the models is closely related to the number of independent variables, such as the number of service trips, lines, and power demand periods as well as the number of timesteps. The arc based-model is a MIP optimization problem, which is a NP-hard problem per se. The order of magnitude is dictated by variables a, b, θ q and x whose sizes are respectively |I| • |T|, |Q|, and |I| 2 , and by constraints (4), (5b), (7), (9a), (9c), and (19a) whose sizes are |I| 2 , |I| • |T| and |Q|. To optimize the operational costs and reduce the peak power demand, we reduce the computational complexity by reformulating the arc-based model into two-stage model. The first stage model involves the column generation  technique, resulting in a runtime complexity of O( by the recursive formulation in each iteration; the second stage model with the timetable shifting strategy has a heavily reduced scale of variables, which is an integral programming problem and can be easily solved by the heuristic rules and an optimization solver.

Numerical experiments
In this section, small instances are first generated to evaluate the performance of the proposed two-stage solution method. To further assess the effectiveness, a real-life case study with a similar setting in Zhang et al. (2021) is conducted. We explore the impact of flexible charging and timetable shifting strategies on optimal solutions to highlight their roles in reducing the total operational costs and moderating peak power demand. All computations are performed on a PC with i7-9850U @ 2.60 GHz CPU and 8 GB RAM. The programming language is Python.

Small instances
We generate a set of illustrative small-sized instances in a planning time horizon from 6:00 to 12:00 according to Zhou et al (2022). These instances have different headways ranging from 5 to 20 min to assess the efficiency of the proposed solution method as shown in Table 3. We set the number of chargers to 1 for the convince of studying the utilization rate in such small cases. Supposed that the energy consumption is 1.35 kWh/km (Gao et al., 2017) per round trip of 20 km in length. A BEB has a battery capacity of 120 kWh, with which a full charge can empower around three trips. The timestep is set as 10 min, meaning that the start charging time and charging duration should be an integer multiple of the timestep. Considering the battery discharges when not in use after a full charge overnight, each BEB is assumed to have a SOC of 0.95 when it starts daily operations. The minimum charging duration is 10 min (Duan et al., 2021); the minimum headway between any two timetabled trips is 5 min. We use the time-dependent electricity prices (Wang et al., 2021) shown in Table 4. For safety concerns, the SOC of the BEB battery should be within the range [H min = 0.2,H max = 0.95]. The fixed cost of employing one BEB is 16.5 $/day, a linear charging profile is adopted for slowing charging at night, and the parameters for the battery degradation are adopted from Zhang et al. (2021). The peak power demand is in two or three order of magnitude while the sum of time shifts is in one order of magnitude. Thus, we set ω 1 to 0.001, with the aim of optimizing the peak power demand while choosing the solution with the smallest time shifts. As shown in Fig. 4, the nonlinear charging profile in CV phase is approximated by a piecewise linear function under fast charging for daily operations. ω 2 is the controlling factor, and set as 1.05.
Thus, we adopt the charging profile of fast charging, which is expressed as Eq. (26) where χ is the charging duration (h) and f(χ) denotes the charging SOC.
We implement the simple heuristic method (Section 4.1.1) as the first baseline to compare the optimal schedule obtained by the two-stage method. According to the battery level, the simple heuristic method always chooses to replenish energy or serve the trip with the earliest departure time. To further validate the effectiveness of our model and method, the model suggested by Zhang et al. (2021) is also solved by the CG technique for comparison as the second baseline. In their model, the charging strategy is semi-flexible. Once a bus visits the charging station, it will be charged to a pre-given electricity level while optimizing the sequence of charging operations at the depot.
For the sake of convenience, we denote these three methods respectively as P1 (simple heuristic method and fixed charging strategy), P2 (CG method with semi-flexible charging strategy), and P3 (our model and two-stage method). For the charging in P1, the start charging time is not flexible and the batteries are charged to a specified electricity level; for P2, the start charging time is flexible and the charging duration is fixed given the remaining electricity volume; and for P3, both the start charging time and charging duration (related to electricity volume) are flexible. The comparations of the results are presented in Fig. 8. The figure shows that the total operational costs and the cost components including fleet, charging, and battery degradation costs increase as the number of trips increases, as expected. The average utilization rate of BEBs fluctuates with the number of trips, where the average utilization rate is derived by the service time of transporting passengers divided by the total operational time. In addition, this figure shows the result of the peak power demand, which is determined by the use of the single charger with an upper bound value of 96 kW.
Compared with the results of P1 and P2, P3 obtains additional total operational costs savings by 17-21% and 7-14% for the six Table 4 Price of the electricity of the day.
Electricity price ($/kWh) Time of the day 0.049 22 pm -7 am 0.121 7 pm -9 am, 11 am -13 pm, 16 pm -18 pm, 20 pm -22 pm 0.172 18 pm -20 pm 0.173 9 am − 11 am, 13 pm -16 pm instances ( Fig. 8(a)), respectively, indicating that the proposed model and solution method can produce better fleet and charging schedules. Compared to the fleet cost obtained by P1 and P2, the costs of P3 are reduced by approximately 31% and 16% respectively ( Fig. 8(b)), implying a more compact fleet schedule. This result is corroborated by a 45% and 20% increase in the utilization rate of BEBs (Fig. 8(c)), respectively. Nevertheless, the charging cost associated with P3 is higher than the results of P1 and P2 (Fig. 8(d)). The results can be explained by the reduction in fleet size but increases in the number of service trips per BEB and electricity consumption for the operations. As for battery degradation, we know that larger SOC variations cause faster battery degradations (L. Zhang et al., 2021;Zhou et al., 2022). As expected, the battery degradation cost from P3 is less than those of P1 and P2 (Fig. 8(e)) because the fullflexible charging strategy in P3 avoids larger SOC variations. Similar to the charging cost, we notice that the battery degradation costs are approximately the same for IB5 and IB6 for the three methods because of the setting of only one charger. Finally, in Fig. 8(f), we see that P1 has the highest peak power demands in the six instances, whereas the peak power demands of P3 are not always lower than those of P2. To demonstrate that the number of used chargers impacts on peak power demand, we further investigate a large case study in Section 5.2.
To evaluate the effects of the length of the timestep, we test the computation time and the efficiency using the six instances. The results are presented in Table 5. Columns 1 and 2 present the problem instance; columns 3 and 4 show the computational results of the two-stage solution method, including the total operational costs and the peak power demand; columns 5 and 6 are the running times to obtain the optimal solution (T_CPU Time) and spending on solving the pricing problem (PP_CPU Time), respectively. With different timesteps for the same instance, the optimal costs appear approximately stable, indicating that the optimal costs are to a small extent influenced by timesteps in these cases. According to the CPU times, we observe that the computation efficiency is largely and negatively affected by the number of trips and positively affected by the timesteps. In other words, longer computation times are associated with more trips and shorter timesteps. In addition, it is found that the size of the timestep also has a significant impact on the peak power demand in the instances. When the timestep is larger, it means more bias in the peak power demand. Summarily, when the timestep is 10, there is no significant difference in the accuracy from a smaller timestep, but the computation time is significantly less.
To balance computation efficiency and accuracy, the timestep of 10 min is adopted in the large case study below.

Real BEB network
In this subsection, we further demonstrate the effectiveness of the proposed model and solution method with a real-world BEB network from Zhang et al. (2021), of which the original timetable has 210 trips as listed in the supplementary document (Table S1). The network has one depot, two terminals, and six lines (numbered from 1 to 6), as depicted in Fig. 9. The lengths of the six lines are 10 km, 12.5 km, 30 km, 17.5 km, 20 km, and 22.5 km, respectively. Lines 1, 2, 4, 5 and 6 are unidirectional, while line 3 is circular. We consider the round trip travel times (terminal to terminal) are 69 min, 86 min, 103 min, 120 min, 138 min, and 155 min, respectively. That is to say, for lines 1, 2, 4, 5 and 6, the travel time of a trip doubles the one-way travel time; for line 3, the travel time of a trip is equal to the one-way travel time. The BEBs are medium-duty buses equipped with 162 kWh batteries and the planning time horizon is set for the whole day including daily operations, electricity replenishment, and overnight charging. The number of fast chargers for electricity replenishment is 15, and a sufficient number of slow chargers with a charging power of 50 kW is provided for overnight charging. Other parameters are consistent with those for the small instances.  Fig. 9. Geographical distribution of selected BEB lines.
The optimal solution can be obtained by the two-stage solution method within 5 h, 2.5 h for each stage. As this is a deterministic and planning problem that does not require real-time responses, we consider it feasible. The optimized fleet schedule and charging operations are shown in Table 6, which includes 45 BEBs and 104 charging operations. The charging details including the start charging time, end charging time and charging volume, and the optimized timetable is shown in the supplementary document (Table S2). The total operational costs are $2207.09, including a fleet cost of $726, a charging cost of $979.28, and a battery degradation cost of $501.81. The optimized peak power demand is 855.36 kW, and the utilization rate is 0.571. The block shifts are [-5, − 5, 1, 5, 2, 1] in minutes for lines 1 to 6, where a negative and a positive number represent moving forward and backward for the whole line, respectively. Also, some differentiated shifts in the timetable are made for 17 trips, as shown in the supplementary document (Table S3).
We compare the results of P3 with those obtained from P1 and P2, as reported in Fig. 10. Compared to P1, P3 has a 22% reduction in the total operational costs, 20% reduction in fleet cost, 1% savings in charging cost, and 45% savings in the battery degradation cost. Thus, it can be concluded that our model and solution method have significant advantages over P1, which is usually adopted by the BEB operators. Compared to P2, it is worth noting that P3 has a 7.1% of total operationals cost savings but a 17% increase in the battery degradation cost. The result indicates that the full-flexible charging strategy involves less total operational costs, while it accelerates the battery aging process and results in higher costs for battery capacity fading. Further, we can see in Fig. 11 that the peak power demand is reduced by 17% and the used number of chargers is reduced by 2. This makes sense because more buses being charged at the Table 6 Fleet schedule and the charging operations of the six lines.   Fig. 13. Comparison of the peak power demand and number of used chargers for the three charging strategies. same time create a higher peak power demand; in other words, reducing the peak power demand is associated with decreasing the number of used chargers. These results in Figs. 10 and 11 confirm that the integrated optimization model and the two-stage solution method together lead to a more efficient fleet and charging schedules, and also have a significant advantage in further reducing the peak power demand and decreasing the number of chargers.
The following analysis further illustrates the advantages of the full-flexible charging strategy. Fig. 12 shows the comparison derived from the proposed two-stage solution method with the three charging strategies. As expected, the full-flexible charging strategy results in the lowest total operational costs, the lowest fleet size, and the highest utilization rate of BEBs. A flexible charging strategy increases the charging opportunity, allows for more trips to be performed, and requires fewer BEBs to finish the trip service, thereby resulting in lower total operational costs and an increased utilization rate. In addition, Fig. 13 shows that the flexible charging strategy uses fewer chargers as well as a smaller peak power demand. The disadvantage is that the full-flexible charging strategy results in a greater battery degradation cost compared to the semi-flexible charging strategy because the higher SOC variation or average SOC causes faster battery capacity degradation .
For illustration purposes, the SOC curves of the three charging strategies of line 2 are presented in Fig. 14, where we set the initial and maximum battery electricity as 95%. Fig. 14(a) shows that the BEBs visit the charging station if the remaining electricity level is too low for the next trip. Obviously, the fixed charging strategy leads to high SOC between 30 and 95% and hence more battery degradations compared with the other two charging strategies. Thus, it is an inferior charging strategy because it does not only occupy more vehicles but also has a serious impact on battery life. Fig. 14(b) and (c) show the SOC curves for the semi-flexible charging strategy and the full-flexible charging strategy, respectively. They both take into account the number of available chargers, TOU plans, the current electricity level, and the power load on the grid. Comparing these two figures, we can see that the semi-flexible charging strategy maintains high SOCs between 70 and 95% in most cases, while the full-flexible charging strategy keeps the SOC between 20 and 65%. Based on the finding that a higher average SOC level results in a faster battery capacity degradation , the full-flexible charging strategy with an average of 42.5% SOC seems to outperform the semi-flexible counterpart with 82.5% average SOC. However, it is not true in terms of the higher battery degradation cost associated with the full-flexible charging strategy. It is because, with the full-flexible charging strategy, the battery SOC varies considerably, including an initial discharge approximately from 95% to 30% or a final charging from 30% to 100% at the end of the operation as seen in Fig. 14(c). Finally, it can also be seen that because of the limited resources (e.g., chargers and power grid) or the search for better charging timing due to the TOU plans, BEBs in the full-flexible charging strategy are enabled to wait for a short time before starting to charge and choose an appropriate charging duration. Overall, the full-flexible charging strategy leads to the more efficient fleet and charging schedules.
To examine the effects of the timetable shifting strategy, we set up a scenario without the consideration of timetable shifting for comparison. As seen in Table 7, regardless of the solution method, the timetable shifting strategy leads to a 17-26% decrease in peak power demand and a reduction in the number of chargers, although the total operational costs increase slightly by 0.3-0.7%. The increase in the total operational costs is mainly due to the increase in charging costs as a result of the TOU plans, while the charging volumes remain the same. To exclude the influence of the TOU plans on the total operational costs, we set a fixed electricity price of 0.15$/kWh in the test. It is found that the total operational costs remain unchanged after applying the timetable shifting strategy but the peak power demand decreased considerably by 26% compared with not applying the strategy. Using the full-flexible charging strategy and keeping other setups the same, Fig. 15 shows the daily power demand for central depot charging before and after applying the strategy, demonstrating that the timetable shifting strategy has a significant impact on reducing the grid load.
Through the above numerical analyses, we conclude that the proposed model and two-stage solution method perform well in generating the fleet and charging schedules. The results show that the full-flexible charging and timetable shifting strategies significantly influence the total operational costs and the peak power demand. For the operators of BEB networks, the design of the timetable, the fleet schedule, and the charging strategy are not only crucial for reductions in the fleet size and number of chargers but also important for the relief of the power grid pressure. We recommend that the operators should apply the flexible charging strategy and adopt the timetable shifting strategies including block shifting and differentiated shifting strategies. In this way, the total operational costs can be saved considerably and the power grid pressure can be relieved, associated with a decrease in number of used chargers.

Conclusions and future work
This paper investigates the optimal fleet and charging schedules with the consideration of timetable shifting, battery degradation, nonlinear charging profile, and TOU plans. We propose an integrated optimization model with high complexity and reformulate the arc-based model into a two-stage model. We further develop a two-stage solution method involving the CG technique and iterative adjustments of timetable shifts to find optimal solutions. Numerical experiments considering small instances and a real-world case are performed to demonstrate the efficiency and effectiveness of the model and solution method. The results show that the proposed model and method can achieve significant reductions in the total operational costs and peak power demand compared to two baselines. To highlight the influence of the full-flexible charging and timetable shifting strategies, sensitivity analyses are conducted. The results show that the strategies contribute to the peak power demand reduction and the increased utilization rate of BEBs. These results confirm the outperformances of the proposed fleet schedule and charging planning model that incorporates the two useful strategies.
In future work, we would make more efforts to improve our work. First, we will further investigate the conjecture that using a larger capacity battery could reduce the battery degradation cost caused by the full-flexible charging strategy. Second, the influence of timetable shifting on passenger demand and preferences (Liao et al., 2013(Liao et al., , 2020 should be considered because timetable shifts may lead to a change in waiting time for passengers and further decrease their travel accessibility Liao, 2021, 2022). Third, some M. Duan et al. realistic conditions, such as the more bus types, variable charging power, wireless charging strategy, and multiple depots in the BEB network should be considerd. Fourth, the coefficients weighing the importance of different components in the objective functions should be calibrated before real-world applications.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
No data was used for the research described in the article. Notes: C = total operational costs; CC = charging cost; θ = peak power demand; ↑ = increase; ↓ = decrease. Fig. 15. Comparison of the power demand with and without the timetable shifting strategy.