Optimal charging plan for electric bus considering time-of-day electricity tariff

Abstract

Purpose

The purpose of this study is to develop an optimization method for charging plans with the implementation of time-of-day (TOD) electricity tariff, to reduce electricity bill.

Design/methodology/approach

Two optimization models for charging plans respectively with fixed and stochastic trip travel times are developed, to minimize the electricity costs of daily operation of an electric bus. The charging time is taken as the optimization variable. The TOD electricity tariff is considered, and the energy consumption model is developed based on real operation data. An optimal charging plan provides charging times at bus idle times in operation hours during the whole day (charging time is 0 if the bus is not get charged at idle time) which ensure the regular operation of every trip served by this bus.

Findings

The electricity costs of the bus route can be reduced by applying the optimal charging plans.

Originality/value

This paper produces a viable option for transit agencies to reduce their operation costs.

Citation

Liu, Y., Wang, L., Zeng, Z. and Bie, Y. (2022), "Optimal charging plan for electric bus considering time-of-day electricity tariff", Journal of Intelligent and Connected Vehicles, Vol. 5 No. 2, pp. 123-137. https://doi.org/10.1108/JICV-04-2022-0008

1. Introduction

The energy crisis and environment pollution caused by transportation has exacerbated in the past few decades (Guo et al., 2019), and the importance of developing clean transportation options has become the general awareness of governments around the world. As an advanced clean transportation alternative, electric buses (EBs) have the features of zero-emissions, low noise level and high driving stability. They can reduce the urban air pollution emissions while also contributing to the energy structure diversity, alleviating the urban environmental pressure and energy crisis accordingly (Li et al., 2009). Therefore, the deployment of EBs has been expanded quickly in recent years, leading to the rapid increase in the numbers of the global EB sales. With China as an example, the numbers of EBs have been over 10 thousand in Beijing, Guangzhou, Shenzhen and some other cities. Until the end of 2020, total 400 thousand EBs operate in China, accounting for 60% of all urban buses, which makes electric bus the type with the largest scale among all bus types.

Electric buses need to get charged during daily operation due to the limitation of short driving range caused by the hurdle of onboard battery capacity, which makes the electricity costs an important component of the operation costs of transit agencies. Hence, creating favorable charging plans to reduce the electricity costs is an effective way to decrease the operation costs. The nonlinear increase of electricity consumption rates observed from the statistical analysis of the real data should be considered when transit agencies make bus scheduling and charge plans. This feature refers to the lower electricity consumption rates per unit of miles at relatively high state of charge (SOC) level and the higher consumption rates per unit of miles at low SOC level.

Besides, Time-of-Day (TOD) electricity tariff is applicable to corporations in some cities. Specifically, 24 hours per day are divided into several periods, and electricity consumption will be billed differently with the periods, being higher in on-peak load periods which are usually daylight hours and lower in off-peak load periods which are usually night hours. It is implemented to incentivize corporations to shift loads to off-peak periods, reducing the operation costs of the electricity system.

Under the TOD electricity tariff, charging buses only at night in nonoperation hours (i.e. 22:00 pm–04:00 a.m.) would lead to the reduction on electricity costs due to the lower electricity rate. However, the battery SOC might not be enough to support all trips per day, causing the service interrupt. In addition, without charging in operation hours, an electric bus has to run at low SOC levels for some later trips with a higher energy consumption rate due to its nonlinear increasing feature, increasing the total electricity consumption and costing more on electricity bills. On the other side, charging buses in daylight operation hours would maintain the SOC around a high level, so the bus would run with a low energy consumption rate, but the electricity rate is higher in the daytime. Therefore, several factors should be considered comprehensively including the varying of battery SOC during the bus operation, the corresponding energy consumption rates and the different electricity rates, so as to generate the optimal charging plans.

More studies of electric buses have been conducted recently, mainly on the optimization of charging infrastructure/electric bus infrastructure planning (He et al., 2019; An, 2020), lifecycle costs evaluation (Lajunen, 2018; Ritari et al., 2020; Ma et al., 2021), the electric bus fleet transition problem (Pelletier et al., 2019; Logan et al., 2020) and the energy management system (Zhang and Yao, 2015; Du et al., 2016; Yang et al., 2017; Liu et al., 2017; Du et al., 2018; Vepsäläinen et al., 2018). Little research effort has been put on the optimization of charging plans. Qin et al. (2016) simulated daily charging patterns and demand charges of a fleet of electric buses in Tallahassee, FL, and identified an optimal charging strategy to minimize demand charges. Chen et al. (2018) proposed an optimal real-time coordinated charging and discharging strategy for a Plug-in electric bus fast charging station with an energy storage system to achieve maximum economic benefits. Chen et al. (2020) proposed a configuration-control integrated strategy to optimize the number of second-used electric bus batteries to minimize the overall equivalent annual cost of EBCS (electric bus charging station). The coordinated charging-discharging power of the echelon battery system and EBBs (electric bus batteries) was optimized to minimize the daily electricity cost of EBCS. Rogge et al. (2018) pointed out that the transition process from conventional diesel to electric buses faces major hurdles caused by range limitations and required charging times of battery buses. They addressed these constraints and provided a methodology for the cost-optimized planning of depot charging battery bus fleets and their corresponding charging infrastructure. Considering TOD electricity tariff, Li et al. (2017) established the load model of the bus FCS (fast charging station), took maximize the net income of the bus FCS system as the objective and then established the economic model of the FCS. Rupp et al. (2020) presented a new methodology for optimizing the vehicles’ charging time as a function of the parameters CO_2eq emissions and electricity costs. Different charging scenarios were defined to analyze the influence of the temporal variability of CO_2eq intensity and electricity price on the environmental impact and economy of the bus.

Two issues from previous studies can be summarized as follows:

1. In most studies, the energy consumption rate of bus batteries during electric bus operation is assumed as a constant, which is determined by electricity consumption models obtained from fitting the laboratory data. Those data, however, have significant differences from real operation data. From the analysis of real collected data, the authors notice the nonlinear discharge characteristic of electric bus batteries, which indicates the nonlinear increase of energy consumption rate with the decrease of battery SOC. This characteristic has an important impact on the optimization of charging plans.

2. The influence of the stochastic trip travel times on charging plans has not been considered in previous studies. Actual bus operation is affected by multiple stochastic factors including traffic flow conditions, queuing of vehicles at intersections and passenger flows at bus stations, which cause the stochastic volatility in trip travel times. This volatility directly affects bus electricity consumption and the formulation of charging plans.

To address the above issues, this paper develops two optimization models for charging plans respectively with fixed and stochastic trip travel times, to minimize the electricity costs of daily operation of an electric bus. The charging time at each bus idle time is taken as the optimization variable. The TOD electricity tariff is considered, and the energy consumption model is developed based on real operation data. An optimal charging plan provides charging times at bus idle times in operation hours during the whole day (charging time is 0 if the bus is not get charged at idle time) which ensure the regular operation of every trip served by this bus. This paper produces a viable option for transit agencies to reduce their operation costs.

The remainder of this paper is organized as follows. Section 2 introduces the formulation and solving algorithm of the charging plan optimization model with fixed trip travel times. Section 3 articulates the formulation and solving algorithm of the charging plan optimization model considering the stochastic volatility in trip travel times. Section 4 presents the verification of the proposed optimization models based on the case of a real bus route and the results of the sensitivity analysis. Section 5 provides some concluding comments and future directions.

2. Model development with fixed trip travel times

2.1 Problem statement

The notations and descriptions of all parameters and variables used in this section are listed in Table 1.

An electric bus is arranged to serve I trips per day according to the transit schedules. Here, one trip is defined as a bus running from the starting station to the terminal station or inversely along the one-way route and the schedule departure time and travel time of each trip are known. Charging stations are deployed at both starting station and terminal station of the bus route and an electric bus can choose to get charged there or not at idle time after each trip. The duration of idle time after trip i is equal to the departure time of trip i + 1 minus the arrival time of trip i, denoted by t_i+1 − (t_i + T_i). The bus can get charged in nonoperation hours after finishing all trips to get prepared for the operation of the next day. An integer variable θ_i (min) represents the charging time after trip i and θ_i = 0 when no charging is arranged.

Twenty-four hours per day are divided into J time periods for the implementation of the TOD electricity tariff, and the electricity rates between adjacent time periods are different. The division moment between the j-th period and the j + 1-th period is tj,j+1c (j = 1, 2, …, J-1) and the boundary of the j-th period [tj−1,jc,    tj,j+1c) is denoted by λ_j. The electricity rate for the charging arrangement after trip i is determined by the period where idle time between the arrival time of trip i and the departure time of trip i + 1 falls into. The bus may idle after trip i crossing two time periods with various electricity rates. Under this circumstance, θ_i equals to the summation of θiF and θiS, which are the charging times respectively in the former and latter period within idle time after trip i.

2.2 Model formulation

To formulate the optimization model of charging plans, we take it as the objective to minimize the total electricity costs ω_total of daily operation of an electric bus and utilize the charging time at each idle time θ_i (or θiF, θiS) as the optimization variable.

No matter the battery gets charged or not after each trip, the requirement of remaining SOC to complete the next trip should be satisfied and manifested by the constraints of the model. To be specific, SOC_i should be sufficient to ensure the whole operation of trip i. Besides, if the charging strategy is taken that a bus gets charged using idle time in operation hours, it is better to maintain the SOC of its on-board battery within a judicious range [SOC_min, SOC_max] to keep the battery in the half-charge/half-discharge status and avoid over-charge and over-discharge, so as to extend the battery life (Jana et al., 2019). Usually, SOC_min takes a value of 20%–30% and SOC_max of 70%–80%.

Derived from the above description, the optimization model of charging plans (Model 1) is established as follows:

Model 1:

(1) min⁡ωtotal=ω1+ω2

(2) s. t.  SOCmin⁡+QiB≤SOCi≤SOCmax⁡, i=1,2,3,…,I

The constraints denoted by equation (2) ensure that the SOC after completing each trip is within the range [SOC_min, SOC_max].

(1) Calculation of ω₁:

(3) ω1=∑i=1I−1Ci

(4) Ci={θi⋅P⋅Mj/60,                    (ti+Ti  )  ∈λj     and    ti+1∈λjθiF⋅P⋅Mj/60+θiS⋅P⋅Mj+1/60,                 (ti+Ti  )  ∈λj     and    ti+1∉λj

where (t_i+T_i) is the arrival time of trip i, t_i+1 is the departure time of trip i + 1, i = 1,2,3,…, I − 1, j = 1,2,3,…, J. The piecewise function illustrated by equation (4) indicates the calculation of the electricity cost for charging after trip i based on the judgement whether (t_i+T_i) and t_i+1 fall into the same period. If they fall into the same period, the energy cost is calculated using one electricity rate (M_j) as shown in the first case of equation (4). Otherwise, the costs in two periods need to be calculated independently with different rates (M_j and M_j+1) and summed together as shown in the second case of equation (4).

The maximum charging time can be either the time needed for the battery to be charged to SOC_max or the duration of the whole idle time. A minimum charging time T_min is set to guarantee the effectiveness of each charging arrangement. Specifically, when the idle time is shorter than T_min, namely, t_i+1 − (t_i+T_i)≤ T_min, the bus will not get charged, which means that θ_i (or θiF, θiS) is 0. Therefore, with regard to the feasibility of charging plans, the charging time variables θ_i, θiF and θiS can take only nonnegative integers and their value intervals each are θ_i ∈{0}U[T_min,min{t_Ei, t_Li}], θiF∈{0}∪[Tmin⁡,min⁡{tEFi,tLFi}] and θiS∈{0}∪[Tmin⁡,min⁡{tESi,tLSi}]. Since the bus stops operating after completing trip I and get charged then in nonoperation hours during off-peak period, i only takes up to I − 1 for θ_i (or θiF, θiS). The calculations of t_Ei, t_Li, t_EFi, t_LFi, t_ESi and t_LSi are expressed as equations (5) – (10):

(5) tEi=(SOCmax⁡−SOCEi)BP×60

(6) tLi=ti+1−(ti+Ti)

(7) tEFi=(SOCmax⁡−SOCEi)BP×60

(8) tLFi=tj,j+1c−(ti+Ti)

(9) tESi=(SOCmax⁡−SOCEi)BP×60−θiF

(10) tLSi=ti+1−tj,j+1c

(2) Calculation of ω₂

(11) ω2=(SOC1−SOCEI)⋅B⋅Mn

(3) Calculation of SOC_i and SOCE_i

SOC_i and SOCE_i are calculated as follows:

(12) SOC1=SOCmax⁡

(13) S O C i = { S O C i − 1 + θ i − 1 ⋅ P / 60 − Q i - 1 B × 100 % ,                         ( t i − 1 + T i − 1     )     ∈ λ j         a n d       t i ∈ λ j   S O C i − 1 + ( θ i − 1 F + θ i − 1 S ) ⋅ P / 60 − Q i - 1 B × 100 % �                 ( t i − 1 + T i − 1     )     ∈ λ j         a n d       t i ∉ λ j �       2 ≤ i ≤ I

(14) SOCEi=SOCi−QiB×100%            1≤i≤I

SOCE_I is expressed as equation (15), deriving from equations (13) – (14):

(15) SOCEI=SOC1+∑i=1I−1θi⋅P/60−∑i=1IQiB×100%

where θ_i is replaced by θiF+θiS if (t_i+T_i)∈λ_j and t_i+1∉λ_j, 1 ≤ i ≤ I − 1.

The electricity consumption of trip i (Q_i), which is requisite in the above calculation, is affected by multiple factors. To model Q_i, we analyzed the collected data and selected three independent variables SOC_i, T_i and µ_i which have the strongest correlation to the consumption. Accordingly, the electricity consumption of trip i is estimated as follows:

(16) Qi=β^1SOCi+β^2Ti+β^3μi+β^0             i=1, 2, …, I

where β^1, β^2, β^3 andβ^0 are regression coefficients from fitting the actual collected data. Due to limitations of space, the data collection and the regression of equation (16) are not presented in this paper and the details will be demonstrated in other papers.

2.3 Solving algorithm

Since Model 1 is a constraint nonlinear optimization model, its accurate optimal solution cannot be obtained from the normal derivation. Therefore, we propose an ant colony optimization algorithm with mutation to solve the model.

Ant colony optimization (ACO) is the simulation of the food-searching behaviors of ant colonies in the natural world. Ants in a colony initially walk through all paths to search for the food and deposit pheromones on the paths after finding the food source. It takes less time for ants to walk between the food source and their nest through the shorter paths, leading to higher frequency and then more travels of ants within a unit of time. Consequently, more pheromones are deposited on shorter paths, on which more subsequent ants are attracted to search. Mapping the abstraction of this process into the ACO algorithm, each path refers to the corresponding feasible solution of the optimization problem and the pheromone on each path represents the assessment of the quality of the corresponding solution, which indicates the degree of optimality of this solution. The higher pheromone on the path, the higher quality of the corresponding solution will be demonstrated. The optimal solution can be achieved by locating the feasible solution related to the path with the highest pheromone. ACO has features of high efficiency and robustness as other parallel optimization approaches (Dorigo and Gambardella, 1997).

To apply ACO to solve Model 1, the feasible solutions from Model 1 are represented by the decision paths of ants, and the solution space is formed by the set of decision paths of the ant colony. Specifically, the decision vectors consisted of θ_i (or θiF, θiS) representing charging times are the decision paths of ants, denoted by D_θ. Since the departure times and travel times of trips are scheduled and known, it is elementary to determine whether each charging time at idle time takes θ_i or the summation of θiF and θiS, which is represented by an auxiliary variable φ_i expressed by equation (17). When φ_i = 1, the charging time is θ_i, and when φ_i = 2, the charging time consists of θiF and θiS. The dimension of D_θ is ∑i=1I−1φi. Each ant constructs one D_θ and all D_θ from ants in the colony comprise the solution space. Since the objective is to minimize ω_total, the lower objective value calculated with the D_θ satisfying all constraints indicates the higher quality of the associated solution, which is closer to the optimal solution. A greater amount of pheromone accumulates on those D_θ associated with high-quality solutions. The decision path with the highest pheromone value will be identified with a positive feedback effect, in turn achieving the optimal solution with the minimal ω_total.

(17) φ i = { 1           ,   ( t i + T i     )     ∈ λ j         a n d       t i + 1 ∈ λ j   2         ,   ( t i + T i     )     ∈ λ j         a n d       t i + 1 ∉ λ j  

However, when the number of decision stages and the number of alternative states in each stage are large, it is difficult to determine the solutions both satisfying constraints and having high quality from a large number of decision paths, since it takes long time to accumulate pheromones and show the divergence on pheromone values among decision paths. To increase the computational efficiency and reduce the searching time, we adapt the mutation procedure from Genetic Algorithm (GA) to the basic ACO and form a new ACO algorithm with mutation. The steps of utilizing the new algorithm are listed as follows:

Step 1: Randomly generate multiple qualified decision vectors D_θ. A decision vector which satisfies all constraints of the optimization model is defined as a qualified decision vector. In this step, no less than (N + N × G) qualified decision vectors D_θ are generated. Here N is the number of ants in the colony (20 ≤ N ≤ 50 usually) and G denotes the maximum number of iterations. Among all generated D_θ, N vectors are provided for all ants in the colony as initial paths and N × G vectors to ensure sufficient solutions in the global search for searching during iterations.

For instance, if I = 10 and φ_i(i = 1,2,3…, I − 1) are all assumed to be 1, then the dimension of a decision vector D_θ like [0, 0, 0, 5, 0, 10, 0, 6, 0] is nine. When each θ_i of it satisfies SOCmin⁡+QiB≤SOCi≤SOCmax⁡, this vector is a qualified D_θ.

Step 2: Initialize decision paths for the ant colony. N vectors from generated D_θ are taken as the initial decision paths, denoted as Dθn(1)  , 1≤ n ≤ N. The values of ω_total associated with Dθn(1) are set as initial pheromone values, denoted by Rθn(1) , 1≤ n ≤ N. Since the optimization in this paper is minimization problem, the lower the value of Rθn(1)  is, the higher quality of the corresponding solution has.

For example, when N = 20, 20 D_θ (e.g. [0, 0, 0, 5, 0, 10, 0, 6, 0]) are taken as the initial decision paths for all 20 ants, and the values of ω_total associated with Dθn(1)   are defined as their initial pheromone values Rθn(1) , consisting of the initialization of ACO (here 1≤ n ≤ 20).

Step 3: Calculating the state transition probability. The Dθn(g)  , 1≤g≤G with the minimum value of Rθn(g)   is selected as the baseline, of which the pheromone value is denoted by R(g), 1≤ g ≤ G. The state transition probability of each Dθn(g) can be obtained by comparing Rθn(g)   with R(g), as expressed in equation (18):

(18) Pn(g)=Rθn(g)−R(g)Rθn(g)   ,       1≤n≤N

Step 4: Conducting the local and global search. To avoid the local optimization, two search algorithms are essential to be utilized simultaneously including local search to find a better solution within a small neighborhood of the current solution of high quality, and global search to restart the search in the global scope when it is of less potential to find a better solution within current neighborhood. Here we set a threshold of the state transition probability P₀ to determine whether to conduct local or global search. Specifically, when P_n(g) ≤ P₀, the local search is carried out since Rθn(g) of the solution is close to R(g), and it is of great potential to find a better solution. When P_n(g) > P₀, the global search is carried out since Rθn(g) of the solution is too large and it is not necessary to search in its neighborhood. Under this circumstance, another qualified vector among D_θ generated in Step 1 is selected for update. The updated Dθn(g) after the local or global search is denoted as Dθn(g)′.

In the local search, a mutation procedure is applied and the basic logic is explained for the following cases with a randomly generated value rand within [0,1]. When the value of θ_i (or θiF,θiS) of the original vector Dθn(g) is 0, the threshold is set as pm₁ where the value of θ_i mutates to T_min. To be specific, if rand ≤ pm₁, the value of θ_i (or θiF,θiS) will mutate to T_min; otherwise, it remains consistent. When the value of θ_i (or θiF,θiS) of the original vector Dθn(g) is T_min, the thresholds are set respectively as pm₂ and pm₃ where the value of θ_i drops to 0 and where the value rises to (T_min+1). If rand ≤ pm₂, the value of θ_i drops to 0; if pm₂ < rand ≤ pm₃, the value rises to (T_min+1); otherwise, the value remains the same. When the value of θ_i (or θiF,θiS) of the original vector Dθn(g) is greater than T_min, the thresholds are set respectively as pm₄ and pm₅ where the value of θ_i rises to (θ_i+1) and where the value drops to (θ_i−1). Similarly, if rand ≤ pm₄, the value of θ_i rises to (θ_i+1); if pm₄ < rand ≤ pm₅, the value drops to (θ_i−1); otherwise, the value remains the same. The value of every variable θ_i (or θiF,θiS) will mutate according to the above logic and update Dθn(g) with Dθn(g)′. The mutation logic is summarized in Table 2.

In addition, the values after mutation should be checked by the boundary conditions of corresponding variables. If the value of θ_i (or θiF,θiS) after mutation surpasses its upper boundary, its value will restore to the original one.

In the global search, a new vector D_θ is randomly selected from all (N + N × G) qualified decision vectors generated in Step 1 and this selected vector is denoted by Dθn(g)′.

Step 5: Determining whether ants move. After the local and global search, the value of ω_total associated with every Dθn(g)′, denoted by Rθn(g)′, is calculated and compared with Rθn(g) of the corresponding original solution Dθn(g). If the new solution is better than the original one ( Rθn(g)′<Rθn(g)), then Dθn(g+1)=Dθn(g)′, indicating that the ant moves to a new path; if the new solution is worse than the original one ( Rθn(g)′≥Rθn(g)), then Dθn(g+1)=Dθn(g) and Rθn(g)′ restores to Rθn(g), indicating no movement of the ant.

Step 6: Updating the pheromone values. The pheromone evaporation rate is denoted by ρ. (1−ρ) in turn represents the residue of the pheromone. The updated pheromone value is achieved by decreasing the pheromone value of the original solution through evaporation and increasing the value associated with the updated solution, as expressed as equation (19):

(19) Rθn(g+1)=(1−ρ)×Rθn(g)+Rθn(g)  ′                     1≤g≤G

Step 7: Conducting the iteration. A threshold for iteration numbers is set as G₀. When g does not reach to the maximum number of iteration G or the optimum varies within the consecutive G₀ times of iteration, go back to Step 3 and let g =  g + 1. When g reaches to G or the optimum remains consistent for the consecutive G₀ times of iteration, stop the iteration and go to Step 8.

Step 8: In the last iteration, the Dθn(g) of the minimum Rθn(g) is the solution with the highest quality, namely, the optimal solution. The corresponding ω_total value is the minimum electricity costs for charging.

3. Model development with stochastic trip travel times

3.1 Model formulation

In the real transit operation, due to the interference from multiple stochastic factors such as road traffic flow states, queuing status at intersections and passenger flows at stations, instead of fixed travel times, the actual travel times are usually reflected by the fluctuations around scheduled travel times. Therefore, the stochasticity of the actual travel time T_i of trip i should be considered and T_i is assumed to follow a normal distribution of which the mean is β_i and the variance is σi2, denoted by Ti∼N(βi,σ)i2.

With the account of the stochasticity of travel times, the objective function and constraints will contain a random variable representing travel times, which leads to the formulation of a typical Stochastic Chance-constrained Programming Model (SCCPM) (Charnes and Cooper, 1959) in this paper. In SCCPM, the solution is allowed to break the constraints to some extent, but the probability that the solution satisfies the constraints should not be lower than a certain confidence level. Thus, Model 2 is achieved by modifying Model 1 according to the above description:

Model 2:

(20) min⁡ ω¯total

(21) s. t.  Pr{ωtotal≤ω¯total}≥β

(22) Pr{SOCmin⁡+QiB≤SOCi≤SOCmax⁡}≥α,    i=1,2,3,…,I

where α and β are the preset confidence levels, ω_total is the objective function, ω¯total is the minimum value of ω_total when the confidence level is greater than or equal to β; Pr{ωtotal≤ω¯total} represents the probability measure for the condition ωtotal≤ω¯total and Pr{SOCmin⁡+QiB≤SOCi≤SOCmax⁡} represents the probability measure for the condition SOCmin⁡+QiB≤SOCi≤SOCmax⁡.

3.2 Solving algorithm

There are two main ways to solve the SCCPM problems. One is to transfer the stochastic programming to the deterministic programming before solving the problems; the other way is to directly utilize intelligent optimization algorithms, such as genetic algorithm, simulated annealing algorithm and neural network algorithm. We choose the second way in terms of the flexibility of intelligent optimization algorithms. A mixed intelligent optimization algorithm incorporating stochastic simulation and immune algorithm (IA) are applied in particular to solve Model 2.

IA is an artificial intelligent optimization algorithm inspired by the principles and processes of biological immune systems as well as those of genetic evolution. With the characteristics of self-adaptation, global convergence and population diversity, IA keeps the diversity of the population and inhibits precocity phenomena comparing to other algorithms, ensuring the achievement of the global optimum and robust solution (Lin et al., 2016; Souza et al., 2016). How to deal with the random variables and random functions is the crucial issue in solving SCCPM problems. To address this issue, stochastic simulation provides an alternative to generate random variables satisfying the constraints and conduct the valuation of random functions. Therefore, we propose a mixed intelligent optimization algorithm incorporating stochastic simulation and immune algorithm. The former is used to control the probability that the solutions satisfy all constraints and the latter to search for the optimal solution.

3.2.1 Stochastic simulation

Based on the requirement of solving Model 2, the estimation algorithm of stochastic simulation is displayed below:

Step 1: Randomly generate values of T_i which follows the distribution N(βi,   σ)i2. Generated values of T_i for total I trips per day form a I-dimension vector of trip travel times by the trip sequence. Generate sT such vectors and comprise a sT × I matrix of random travel times. It is advised to set sT as 100, which is generally sufficient to describe the stochasticity of trip travel times.

To illustrate, consider the case where I = 10 and sT = 100. [30, 36, 32, 34, 30, 35, 29, 31, 40, 35] is an example of total 100 trip travel time vectors in this case, consisting of a 100 × 10 trip travel time matrix.

Step 2: Randomly generate values of charging time variable θ_i (or θiF,θiS) satisfying constraints based on how long the bus idles after trip i, consisting of sX samples for decision vectors D_θ of charging times. The dimension of D_θ is ∑i=1I−1φi, as illustrated in Section 2.4. sX vectors form a sX×∑i=1I−1φi matrix of charging times. It is advised to set sT from 50,000 to 90,000, to ensure that the number of qualified decision vectors D_θ meets the requirement for the number of D_θ in the following immune algorithm implementation. The value of sT can also be modified based on real cases.

Step 3: Calculate the probability α₀ that SOCmin⁡+QiB≤SOCi≤SOCmax⁡ is satisfied under the stochastic trip travel time matrix for each D_θ. If α₀ > α, D_θ is qualified and will be saved.

Step 4: Calculate ω_total under the stochastic trip travel time matrix for every qualified D_θ and obtain sT values of ω_total for every D_θ. For each D_θ, sort its ω_total values in ascending order and take the β×sT-th value as its value of ω¯total, so as to make it satisfy the constraint Pr{ωtotal≤ω¯total}≥β under all random generated trip travel times.

Step 5: Deliver no less than (NP+NP×H2) qualified D_θ, the values of α₀ and ω¯total under each D_θ for seeking the optimal solution via immune algorithm later. Here, NP is the population size in the algorithm (50 ≤ NP ≤500 usually), H is the maximum number of iterations in the IA. Among all qualified D_θ, NP of them are to satisfy the demand of initial population on the number of D_θ and NP×H2 are to satisfy the need of the following population mutation on the number of D_θ. Besides, it is an option to adjust the value of sX in Step 2 to control the number of delivered qualified D_θ.

3.2.2 Immune algorithm

In IA, the objective problem to be solved is indicated by antigen, the feasible solution of the objective problem is indicated by antibody and the quality of the feasible solution is indicated by affinity. Cloning, mutation and selection are conducted on individuals with higher affinity to form new populations. The optimal solution will be achieved by repeating this process. In our optimization problem, the antigen indicates minimizing the total electricity costs for a whole day while antibodies indicate the decision vectors D_θ consisting of variables θ_i (or θiF, θiS) which represent charging times at idle times in the whole day. The detailed process of implementing immune algorithm is listed as follows (Hong et al., 2000):

Step 1: Generate an initial population of antibodies based on the qualified vectors D_θ output through the stochastic simulation. Specifically, NP vectors are recruited from the total (NP+NP×H2) samples of D_θ to form the initial population, denoted by Dθ,1np, 1≤np≤NP. The ω¯total value of each antibody is the affinity which measures the matching between the antibody and the antigen, denoted by ω¯total,1np,      1≤np≤NP.

Step 2: Calculate similarity between all antibodies. The affinity between antibodies should be obtained first by calculating the Hamming distance as illustrated in equations (23)–(24), since the charging times θ_i (or θiF, θiS) are discrete variables:

(23) HMh(Dθ,hnp,    Dθ,hj)=∑k=1∑i=1I−1φi∂k,    1≤np≤NP, 1≤j≤NP, 1≤h≤H

(24) ∂k={1,         if       Dθ,hnp,k=Dθ,hj,k0,       otherwise

where Dθ,hnp and Dθ,hj represent antibody np and antibody j of the h-th generation, HMh(Dθ,hnp,    Dθ,hj) is the affinity between antibody np and antibody j of the h-th generation and Dθ,hnp,k and Dθ,hj,k are, respectively, the k-th dimension of antibody np and antibody j. ∂_k is a binary variable specifying whether the k-th element of antibody np and antibody j are equal to each other.

The similarity between antibodies is then determined by comparing the affinity and a predefined similarity threshold, as demonstrated in equation (25):

(25) Sh(Dθ,hnp,    Dθ,hj)={1,         HMh(Dθ,hnp,    Dθ,hj)≥ζ0,       HMh(Dθ,hnp,    Dθ,hj)<ζ  

where Sh(Dθ,hnp,    Dθ,hj) is a binary variable representing the similarity between Dθ,hnp and Dθ,hj, and ζ is the predefined similarity threshold. Two antibodies are determined to be similar ( Sh(Dθ,hnp,  Dθ,hj)=1) if the number of identical elements of two antibodies is higher than or level with the threshold, as stated by the first clause of equation (25). This means, in the perspective of our problem, that if the number of charging time variables which take the identical values by the two decision vectors is large enough, then these two vectors are determined to be similar.

Step 3: Calculate the antibody density denh(Dθ,hnp) by equation (26):

(26) denh(Dθ,hnp)=1NP∑j=1NPSh(Dθ,hnp,    Dθ,hj)

Step 4: Calculate the antibody incentive simh(Dθ,hnp). In this algorithm, the antibody incentive is defined to measure the antibody quality. Better incentive is achieved with higher affinity with the antigen and lower antibody density. In our problem, a lower ω¯total,hnp  value indicates the higher affinity with the antigen. Low incentive reflects the high quality of antibodies, which should be motivated to get cloned and mutated. Thus, the incentive can be calculated by equation (27):

(27) simh(Dθ,hnp)=a×ω¯total,hnp +b×denh(Dθ,hnp)

where a and b are the positive coefficients of the incentive.

Step 5: Select antibodies based on antibody incentive and clone. Select the first NP/2 antibodies among the population sorted by the antibody incentive in the descending order and get each of them cloned for the number of Nc1.

Step 6: Mutate the cloned antibodies except the parent ones. The identical mutation procedure is implemented in this step with the one used in solving Model 1, which will not be repeated here.

Step 7: Reselect antibodies from the mutated cloned ones. Calculate ω¯total,hnp , HMh(Dθ,hnp,    Dθ,hj), Sh(Dθ,hnp,    Dθ,hj), denh(Dθ,hnp) and simh(Dθ,hnp) of the cloned antibodies after mutation. Based on the simh(Dθ,hnp) values, the antibody with the highest quality among the cloned antibodies and the parent antibody is selected to remain. That is, totally NP/2 antibodies remain.

Step 8: Recruit new antibodies. Recruit another NP/2 antibodies from (NP+NP×H2) qualified samples via stochastic simulation.

Step 9: Form a new population. Combine the NP/2 antibodies remaining after cloning, mutation and selection and the other NP/2 antibodies newly recruited to form a new population of which size is NP.

Step 10: Repeat the iterative process. A threshold for iteration numbers is set as H₀. When h does not reach to the maximum number of iteration H or the optimum varies within the consecutive H₀ times of iteration, go back to Step 2 and let h = h +1. When h reaches to H or the optimum remains consistent for the consecutive H₀ times of iteration, stop the iteration and go to Step 11.

Step 11: In the last iteration, the solution with the minimal simh(Dθ,hnp) is the one with highest quality, which is the optimal solution for the charging plans. The corresponding ω¯total,hnp value is the minimal electricity costs for charging.

4. Case study

4.1 Data investigation

We collected actual operation data of a bus route in a city and conducted the numerical analysis of one electric bus in this route to verify the effectiveness of the proposed models. Basic information of bus operation is displayed here. The operation time of this route is 5:00–23:00 and the bus runs a total of 19 trips per day (I = 19). The travel length of the route is 14 km per direction. The battery capacity (B) of the bus is 162 kWh. The charging power (P) of a charging station is 120 kW. The departure timetable is illustrated in Table 3, which also provides the average ambient temperatures of each trip on November 2^nd, 2020, as well as the coefficients of the distribution function that the actual trip travel times T_i follow obtained by fitting the collected data. Table 4 lists the three electricity rates of four periods per day. Periods are, namely, based on their electricity rates as off-peak, mid-peak and on-peak periods from low to high rate.

With the collected operation data, the energy consumption model of trip i run by the bus is obtained using the weighted least squares method, which is expressed as follows:

(28) Qi=−3.0SOCi+0.27Ti+0.09μi+0.85          i=1, 2,…, I

4.2 Result analysis and assessment

4.2.1 Verification of model 1

To prolong the battery lifespan, the SOC_min is set as 30% and the SOC_max as 80%, which means that the battery SOC should always be maintained within 30%–80%. The minimum charging time (T_min) is set to be 5 min.

The ant colony optimization algorithm with mutation is applied to solving Model 1. In this numerical case, the dimension of qualified decision vectors D_θ is 18, namely, ∑i=1I−1φi=18. Other parameters are taken as follows: the number of ants N = 50, the maximum number of iterations G = 500, the threshold for iteration numbers G₀ = 15, the threshold of state transition probability p₀ = 0.2, thresholds used in the mutation procedure in the local search pm₁ = 0.5, pm₂ = pm₄ = 0.4, pm₃ = pm₅ = 0.8 and the pheromone evaporation rate ρ = 0.9. The optimal charging plan was achieved after 212 times of iteration with the minimal ω_total for the whole-day bus charging, which is demonstrated in Table 5. SOC_i, SOCE_i and the electricity cost of every trip are also provided by Table 5.

It is observed from Table 5 that the bus gets charged four times in the operation during 5:00–23:00, which are 8-min charging after trip 3, 9-min charging after trip 8, 20-min charging after trip 17 and 20-min charging after trip 18. Besides, the bus is arranged to charge for another 6 min after completing all trips of the day, so as to increase the SOC to SOC_max for the operation of the next day.

In this case, the scheduled arrival time of trip 17 is 21:10 and the off-peak period is from 21:00 to 8:00 of the next day, which indicates that the charging after trip 17, trip 18 and that in the nonoperation hours are arranged in the off-peak period. It is particularly noted here that the optimal solution is not unique when the charging time after trip 17 (θ₁₇), the charging time after trip 18 (θ₁₈) and the charging time in the nonoperation hours (θ_f) are only allowed to take nonnegative integers. As long as θ₁₇, θ₁₈ and θ_f satisfying constraints [equations (29)–(31)], all possible value combinations of these three variables can lead to the electricity cost of 53.28 CNY.

(29) 11≤θ17≤20,        θ17∈Z

(30) 0≤θ18≤20 ,       θ18∈Z

(31) θ17+θ18+θf=46,     θ17∈Z,   θ18∈Z,   θf∈Z

Two solutions that θ₁₇ = 11, θ₁₈ = 0, θ_f = 35 and θ₁₇ = 20, θ₁₈ = 20, θ_f = 6, which satisfy equations (29)–(31), both lead to the minimum electricity cost at 53.28 CNY. The nonlinear discharge characteristic of batteries suggests lower electricity consumption rates at higher SOC levels. The second solution taking full use of every idle time for charging, which helps the battery keep high SOC levels to the maximum extent, could theoretically reduce the total energy consumption, bringing the lower electricity bill. However, this advantage is not reflected by the electricity bill, as revealed by the identical objective values of the two solutions. The reason is that the potential values of all charging time variables are restricted to integers considering the feasibility of the optimized charging plan. To be specific, to achieve the minor decrement in the electricity consumption by charging the bus for the whole duration of idle time after trip 17 and trip 18 (i.e. θ₁₇ and θ₁₈ take their maximum values), the corresponding reduction in the value of another variable (i.e. θ_f) is less than 1 min in this case. The reduction being rounded to the nearest integer due to the integer constraints of variables eliminates the minor improvement in the optimization objective.

To better illustrate, we relax the integer constraint of θ_f. Consider θ₁₇ = 11 and θ₁₈ = 0, where θ_f = 35 can be obtained. With this solution, the total charging time is 63 min and the electricity cost is 53.28 CNY. If θ₁₇ = 20 and θ₁₈ = 20, which indicates the full utilization of idle time after trip 17 and 18, θ_f = 5.3 can be achieved. With this solution, the total charging time is 62.3 min and the electricity cost is 52.76 CNY. Comparing to the integer solution where θ₁₇ = 20, θ₁₈ = 20, and θ_f = 6, the charging time in the nonoperation hours (θ_f) decreases by 0.7 min and the electricity cost reduces by 0.52 CNY.

The charging plan provided by Table 5 is labeled as Plan A. Without the comprehensive optimization model considering multiple factors, the following two charging plans are usually applied for electric buses:

Plan B: A normal charging strategy is to charge the bus at the idle time in off-peak operation hours. The bus is usually arranged to get charged once when it idles longer in off-peak operation hours. Otherwise, the energy requirements for smoothly completing all trips or keeping the battery SOC at its best status might not be fulfilled. Based on this strategy, a 26-min charging time will be arranged at a 40-min idle time after trip 8, which is observed in Table 3. Besides, it is required to charge the bus in nonoperation hours for 37 min. Under this charging plan, the total electricity bill is 65.86 CNY.

Plan C: Another normal charging strategy is to charge the bus mainly in the off-peak period during which the electricity rate is low. To be specific, idle time in the off-peak period is utilized as fully as possible for the bus charging. The bus will not get charged in any mid-peak or on-peak periods unless the SOC is not high enough to complete the next trip or drops out of the best range of SOC.

In this numerical case, idle times after the first three trips and the last three trips (i = 17. 18, 19) are all in the off-peak period. According to this strategy, Plan C starts with an 8-min charging time at the idle time after trip 3. No charging is scheduled in the following mid-peak/on-peak periods until the end of trip 14 when the SOC is going to be lower than SOC_min. The bus will get charged for 15 min after trip 14. Next charging is arranged for 20 min after trip 17 when it is off-peak period. Finally, the bus will get charged for another 21 min in the nonoperation hours. The total electricity cost for the whole day is 75.56 CNY under Plan C.

The comparison among the above three charging plans is displayed in Table 6.

As observed from Table 6, comparing to Plan B, Plan A helps to save 12.58 CNY, which is around 19.1%, on the electricity cost. Similar result is revealed comparing to Plan C that Plan A saves 22.88 CNY, which is 29.5%, on the electricity cost. Therefore, Plan A effectively reduces the electricity cost and the total charging time for a whole day, verifying the effectiveness of Model 1.

4.2.2 Verification of model 2

Model 2 is solved by the proposed mixed intelligent algorithm incorporating stochastic simulation and immune algorithm. In the stochastic simulation phase, sizeT is set as 100 and sizeX as 80,000, under which the number of generated D_θ is sufficient for the implementation of immune algorithm later. α and β in constraints (20) and (21) both take 0.9. When applying immune algorithm, the size of the population NP is set as 100, the maximum number of iteration H is 300, the threshold for iteration numbers H₀ is 15, the dimension of vectors in the population is 18, the coefficients of the incentive are set as a = 1, b = 10 and the number of cloned antibodies from each parent antibody Nc1 is 10. Besides, the threshold of similarity ζ is set as 17, which means that if no less than 17 elements of two antibodies are identical, then these two antibodies are determined to be similar. The values of all thresholds used in the mutation procedure are the same with those in solving Model 1.

Via the mixed intelligent algorithm, the optimal solution and the corresponding minimum electricity cost ω¯total were achieved after 160 times of iteration, as illustrated in Table 7.

With modeling the stochasticity of trip travel times in Model 2, the optimal charging plan consists of 8-min charging after trip 3, 11-min charging after trip 8, 11-min charging after trip 17 and 34-min charging in nonoperation hours. Comparing to the optimal solution of Model 1, the electricity bill increases from 53.28 CNY to 55.50 CNY. Under the circumstance with stochastic travel times, the optimal solution in terms of θ₁₇, θ₁₈ and θ_f is still not unique. When θ₁₇, θ₁₈ and θ_f satisfy equations (32)–(34), the electricity bill remains at the minimum value of 55.50 CNY.

(32) 11≤θ17≤20 ,        θ17∈Z

(33) 0≤θ18≤18  ,      θ18∈Z

(34) θ17+θ18+θf=45 , θ17∈Z, θ18∈Z, θf∈Z

Comparing to one optimal solution of Model 1 that θ₃ = 8, θ₈ = 9, θ₁₇ = 11 and θ_f = 35, the largest difference occurs in the increase of θ₈ from 9 min to 11 min. In the setup of this case study, with the charge power of 120 kW, the energy obtained from 2-min increased charging can increase 14.8-min traveling of an electric bus, of which the electricity consumption rate is around 0.27 kWh/min. Due to the tradeoff between the positive and negative fluctuations of trip travel times around their scheduled values, the extra electricity power supporting 14.8-min traveling which is achieved from the extra 2-min charging from the charging plan from Model 2, should be sufficient for actual bus operation considering the effects of the stochasticity of travel times. Therefore, the optimal solution from Model 2 is practicable.

To test how well the optimal solution satisfies the constraints, the simulation runs 50 times with random travel times based on the solution. 100 units of travel times for each trip per day are randomly generated for each simulation run. The probability value α₀ is estimated by counting the number of units which satisfy the constraint SOCmin⁡+QiB≤SOCi≤SOCmax⁡ in each run. Among the total 50 runs, the minimum α₀ is 0.97 and greater than 0.9, demonstrating that Pr{SOCmin⁡+QiB≤SOCi≤SOCmax⁡}≥α (α = 0.9) is well satisfied and validating the optimality of the solution from Model 2.

We also conduct the simulation 10 times based on the solution from Model 1, which is formulated with fixed travel times, to test whether this solution is applicable with stochastic travel times. The results are demonstrated in Table 8.

It is noticed from Table 8, that the maximum probability among 10 runs that SOCmin⁡+QiB≤SOCi≤SOCmax⁡ is satisfied is merely 0.64 based on the charging plan obtained from Model 1, specifying no applicability of the plan under the circumstance with stochastic travel times.

4.2.3 Sensitivity analysis

(1) The effect of total trip number I on charging plans

A bus in urban areas usually serves 10–20 trips per day. In the above case, we take I = 19, and cases where I =11, 12, …, 18 are supplemented here to investigate the effects of different trip numbers on charging plans. Normally the departure time of the first trip is no later than 7:00 while the arrival time of the last trip is no earlier than 17:00. The earliest and latest trips in Table 3 are gradually removed based on this condition to form the cases where I =11, 12, …, 18. The optimal charging plans of all cases are listed in Table 9 and the comparison among them are displayed in Figure 1.

Obviously, the accumulated operation mileage reduces with the decrease of the total trip number, as well as the total electricity consumption. The charging times at idle time gradually decrease consequently, as illustrated in Table 9. Besides, the idle times when the bus is arranged to get charged are identical, which are the last idle time (7:30) in the morning off-peak period, the first idle time (12:20) in the mid-peak period and the idle times (21:10, 22:05) in the night off-peak period. When the total trip number per day is lower than 15, the bus only needs to get charged once (12:20) in the operation hours. When the trip number is lower than 11, under fixed travel time condition, the bus only gets charged in the nonoperation hours during the off-peak period, without the requirement of charging in the operation hours, whereas considering the stochasticity of travel times, the bus still needs to get charged for 5 min in the operation hours.

Comparing the optimal solutions of Model 1 and Model 2, the effect of stochasticity of travel times on the optimal charging plans is reflected by the longer charging time of the latter solution, which considers the stochasticity, than that of the former one, which satisfies all constraints as well as possible, ensuring the whole-day operation while leading to the higher electricity bill. To be more detailed based on Table 5, the amount that the electricity cost of the latter solution exceeds that of the former one varies from 1.48 to 4.44 CNY, which accounts for 3.33%–15.00% of the cost under each case. Large gaps in the cost occurring in cases where the trip number of 17 and 16, are, respectively, 2.96 CNY and 3.70 CNY. The peak value shows in the case of 11 trips which is 4.44 CNY, accounting up for 15% of the cost. The least gap is 1.48 CNY in cases where the trip number is 12, 13 and 14.

Therefore, the effect of the stochastic of travel times on the electricity cost varies with different trip numbers. Generally speaking, with more trips served per day (i.e. I >15 in this study), the effect on the total electricity consumption is significate enough to influence the charging scheduling, demonstrating the necessity of using Model 2 to achieve the optimal charging plan under the condition of larger trip numbers. On the other side, with less trips served per day (i.e. I <15 in this study), the minor effect on the charging scheduling is revealed by the similarity between two solutions from Model 1 and Model 2, which indicates that Model 1 is eligible to generate the optimal charging plan under the condition of smaller trip numbers. An exception occurs when the trip number is 11, where large differences show in both electricity cost and optimal charging plan. When the bus does not need to get charged in operation hours according to the solution from Model 1, it is necessary to verify the result further through Model 2. Under this case, the stochasticity of travel times should be considered, which means that Model 2 should be applied for charging scheduling to ensure regular operation of the bus.

(2) The effect of battery capacity B on the charging plans

The battery capacity of the bus has a significant impact on its driving range, affecting charging plans consequently. A large battery capacity ensures a long driving range without multiple times of charging, while a small capacity can only support a short driving range of the bus, which needs to get charged in a high frequency to guarantee the normal transit operation. The capacity of a normal bus battery is around 140–220 kWh. Within this range, we take the step of 10 kWh to form multiple cases to explore the effects of battery capacity on the charging plans. The achieved optimal charging plans as well as the electricity costs are demonstrated in Table 10, and the variation of the electricity cost in Figure 2.

As observed from Table 10, when B is within 140–190 kWh, the charging time after trip 8 declines with the increase of B, as well as the total electricity cost. The charging time after trip 8 reduces by 2–3 min with the 10 kWh increment of the B. When B is within 210–220 kWh, the bus only needs to get charged after trip 3, trip 17 and trip 18 in off-peak operation hours, avoiding the charging scheduled after trip 8 in the mid-peak period, which indicates that the battery capacity of 210–220 kWh is large enough for regular bus operation.

Another comment suggested by Table 10 is that the charging plan from Model 2 requires longer charging time and leads to high electricity cost comparing to that from Model 1. The gaps between the costs of two plans are between 0.74–2.96 CNY, where the peak occurs with the battery capacity being 140 kWh while the valley occurs with the capacity being 210 kWh. Besides, a decline of the gaps with the increase of the battery capacity is displayed in Figure 2, indicating that more attention should be paid upon the effect of stochasticity of travel times on formulating charging plans when the capacity is lower.

(3) The effect of electricity rate gaps of different load periods on charging plans

With implementing the TOD electricity tariff, there is a high correlation between the optimal charging plans and the electricity rate gaps of different load periods. In our study, three electricity rates are taken including off-peak rate, mid-peak rate and on-peak rate from low to high. Due to the nonlinear discharge characteristic of batteries, the bus is more likely to be arranged to get charged at idle time in operation hours to reduce the total energy consumption and the electricity cost, if the gaps are narrowed between off-peak rate (applying to the electricity use in nonoperation hours) and mid-peak rate or on-peak rate (applying to the electricity use in operation hours). We conduct the sensitivity analysis of different gaps by keeping the off-peak rate and reducing other rates. To be specific, as shown in the first column in Table 11, three cases are formed by decrementing the gap between off-peak and mid-peak rate, and the gap between mid-peak and on-peak rate by the step of 0.1 CNY/(kWh) based on the original case (the first row). One more case is generated by setting a 0.1 CNY/(kWh) gap between mid-peak and off-peak rate and the same gap between on-peak and mid-peak rate. Table 11 lists these cases and corresponding optimal charging plans.

The charging time after trip 8 increases by two minutes in both optimal charging plans from Model 1 and Model 2 when the mid-peak rate is adjusted from 0.64 to 0.54 CNY/(kWh) and the gap reduces to 0.17 CNY/(kWh) between mid-peak and off-peak rate, indicating the intension and feasibility to reduce the total energy consumption and the total electricity cost by prolonging charging times in mid-peak operation hours. By testing with smaller steps from 0.64 to 0.54 CNY/(kWh), we obtain the thresholds of the gap and mid-peak rate where the charging time after trip 8 starts to change, which are 0.55 CNY/(kWh) for mid-peak rate and 0.18 CNY/(kWh) for the gap. The thresholds remain the same for two models. Therefore, when the gap between mid-peak and off-peak rate drops to 0.18 CNY/(kWh), the optimal charging plan will change no matter the stochasticity of travel times is considered or not by the optimization model.

The charging time after trip 8 increases sharply from 11 min to 16 min in optimal charging plan from Model 1 when the mid-peak rate is adjusted from 0.54 to 0.44 CNY/(kWh) and the gap reduces to 0.07 CNY/(kWh) between mid-peak and off-peak rate. By testing with smaller steps from 0.54 to 0.44 CNY/(kWh), it is revealed that the charging time after trip 8 changes once more exactly at the point where the mid-peak rate is 0.44 CNY/(kWh) and the gap is 0.07 CNY/(kWh). For the optimal charging plan from Model 2 which models the stochasticity of travel times, the charging time after trip 8 remains consistent although the mid-peak rate reduces to 0.44 CNY/(kWh) and the gap reduces to 0.07 CNY/(kWh).

With regard to the modification of on-peak rate, it is observed that until the gap between off-peak and on-peak period drops to 0.2 CNY/(kWh) as the last case shown, any charging time in on-peak period is still not included in either optimal charging plans of Model 1 or Model 2. Since the gap less than 0.2 CNY/(kWh) between off-peak and on-peak period is not realistic in the real TOD electricity tariff application, a conclusion can be drawn that bus charging will avoid being arranged in on-peak periods in terms of reducing the electricity cost although the battery has the nonlinear discharge characteristic.

From the analysis of the above three factors, it is concluded that trip number, battery capacity and electricity rate gaps of different load periods all have significant impacts on the formulation of charging plans for electric buses. The differences between the electricity costs of optimal charging plans from Model 1 and Model 2 are, respectively, 1.48–4.44 CNY, 1.06–2.22 CNY and 0.74–2.96 CNY, with one factor varying and the other two remaining consistent. This indicates that the stochasticity of travel times performs a higher influence on the charging plans under various trip numbers than that under various battery capacities and electricity rate gaps. Therefore, when the number of trips served per day is changed, it is suggested to modify the charging plan based on the optimization model considering the stochasticity of travel times.

5. Conclusions

With the implementation of the TOD electricity tariff, we propose two optimization models respectively with fixed trip travel times and stochastic trip travel times, of which the objectives both are minimizing the electricity costs of charging a bus for the whole-day operation and the optimization variables are the charging times at idle time after each trip. The numerical test based on a real electric bus is conducted to verify our models, followed by the sensitivity analysis of three factors including trip number, battery capacity and electricity rate gaps of different load periods. The main conclusions are summarized as follows:

• From the model without considering the stochasticity of travel times, the optimal charging plan costs 12.58 CNY and 22.28 CNY less in electricity use than the two plans formulated based on two normal charging strategies, saving 19.1% and 29.5% on the electricity bill, verifying the effectiveness of the proposed model.

• The stochasticity of trip travel times affects charging plans. In summary, when considering this stochasticity, the total charging time and the electricity cost will increase. The more trips a bus serves per day and the smaller its battery capacity is, the greater influence the stochasticity has on charging plans, the more attention should be paid on the stochasticity.

• Even with the nonlinear discharge characteristic of batteries, it will not be arranged for a bus to get charged in on-peak periods to reduce the total energy consumption and the electricity cost. This feature, however, has a significant impact on the charging schedules in mid-peak periods.

In this study, the optimization model addresses the charging plan for only one electric bus of a transit route. More research efforts will be focused on the studies of the optimization of coordinated charging plans for buses of one route or multiple routes.

Figures

Figure 1

The minimal electricity cost varying with different trip numbers

Figure 2

The minimum electricity cost with different battery capacities

Table 1

Parameter and variable descriptions

Notations	Descriptions	Unit
I	the number of trips running in one-day operation
i	the serial number of trips, 1 ≤ i≤ I
ti	the scheduled departure time of trip i
SOCi	the battery SOC at the start of trip i	%
SOCEi	the battery SOC at the end of trip i	%
Ti	the travel time of trip i	min
Qi	the electricity consumption of trip i	kWh
µi	the average ambient temperature during trip i	°F
P	charge power	kW
θi	the charging time at idle time after trip i without crossing various tariff time periods	min
θiF	the charging time at idle time after trip i in the former period when crossing various tariff time periods	min
θiS	the charging time at idle time after trip i in the latter period when crossing various tariff time periods	min
J	the total number of electrical tariff time periods per day
j	the serial number of electrical tariff time periods
Mj	the electricity rate in the j-th tariff time period	CNY/(kWh)
Mn	the electricity rate in nonoperation hours	CNY/(kWh)
tj,j+1c	the division moment between the j-th period and the j + 1-th period
λj	the time boundary of the j-th period
Tmin	the minimum charging time	min
SOCmin	the lower bound of SOC	%
SOCmax	the upper bound of SOC	%
ωtotal	electricity costs for charging per day	CNY
ω1	electricity costs for charging at idle time in operation hours per day	CNY
ω2	electricity costs for charging in nonoperation hours per day	CNY
B	battery capacity	kWh
tEi	the time needed for the battery being charged from SOCEi to SOCmax without crossing various tariff time periods	min
tLi	the duration of idle time between trip i and trip i + 1 without crossing various tariff time periods	min
tEFi	the time needed for the battery being charged from SOCEi to SOCmax in the former period when crossing various tariff time periods	min
tLFi	the duration of idle time in the former period when crossing various tariff time periods, which equals tj,j+1c−(ti+Ti)	min
tESi	the time needed for the battery being charged from the SOC at tj,j+1c to SOCmax in the latter period when crossing various tariff time periods	min
tLSi	the duration of idle time in the latter period when crossing various tariff time periods, which equals ti+1−tj,j+1c	min
Ci	the electricity cost for charging after trip i	CNY

Table 2

The mutation logic of the solving algorithm

Original value of θi (or θiF,θiS)	Value of rand	Value of θi (or θiF,θiS) after mutation
0	rand ≤ pm1	Tmin
Tmin	rand ≤ pm2	0
Tmin	pm2 < rand ≤ pm3	Tmin+1
>Tmin	rand ≤ pm4	θi+1
>Tmin	pm4 < rand ≤ pm5	θi−1

Table 3

Departure timetable of the electric bus

Trip	Scheduled departure time	Scheduled travel time/min	Scheduled idle time/min	Ambient temperature/°F	Coefficient β	Coefficient σ
1	5:00	35	15	32.0	35	1.4
2	5:50	35	20	35.6	35	1.4
3	6:45	45	10	35.6	45	1.7
4	7:40	50	5	35.6	50	1.4
5	8:35	45	20	33.8	45	2.1
6	9:40	40	20	37.4	40	1.5
7	10:40	40	20	46.4	40	1.4
8	11:40	40	40	51.8	40	1.6
9	13:00	40	15	53.6	40	1.3
10	13:55	40	15	51.8	40	1.7
11	14:50	45	15	50.0	45	1.9
12	15:50	50	10	50.0	50	2.3
13	16:50	50	5	44.6	50	1.6
14	17:45	50	15	44.6	50	2.0
15	18:50	40	15	42.8	40	1.5
16	19:45	35	15	42.8	35	1.4
17	20:35	35	20	39.2	35	1.4
18	21:30	35	20	39.2	35	1.3
19	22:25	35	—	37.4	35	1.3

Table 4

Time of day electricity rates

Time period	Electricity rate/ (CNY/kWh)
8:00–12:00	1.31
12:00–17:00	0.74
17:00–21:00	1.31
21:00–the next 8:00	0.37

Table 5

Optimal charging plan under model 1

Trip	SOCi/%	SOCEi/%	Charging time/min	Electricity cost/CNY
1	80.00	76.80	0	0
2	76.80	73.73	0	0
3	73.73	68.94	8	5.92
4	78.81	73.28	0	0
5	73.28	68.39	0	0
6	68.39	64.42	0	0
7	64.42	60.86	0	0
8	60.86	57.51	9	13.32
9	68.62	65.51	0	0
10	65.51	62.25	0	0
11	62.25	58.00	0	0
12	58.00	52.83	0	0
13	52.83	47.29	0	0
14	47.29	41.65	0	0
15	41.65	37.47	0	0
16	37.47	34.05	0	0
17	34.05	30.38	20	14.80
18	55.07	51.79	20	14.80
19	76.49	73.50	0	0
Nonoperation hours	73.50%→80.00%		6	4.44
Total charging time is 63 min and ωtotal =53.28CNY

Table 6

Comparison among two normal charging plans and proposed charging plan

Charging plan	Plan A	Plan B	Plan C
θ1/min	0	0	0
θ2/min	0	0	0
θ3/min	8	0	8
θ4/min	0	0	0
θ5/min	0	0	0
θ6/min	0	0	0
θ7/min	0	0	0
θ8/min	9	26	0
θ9/min	0	0	0
θ10/min	0	0	0
θ11/min	0	0	0
θ12/min	0	0	0
θ13/min	0	0	0
θ14/min	0	0	15
θ15/min	0	0	0
θ16/min	0	0	0
θ17/min	20	0	20
θ18/min	20	0	0
θf/min	6	37	21
Total charging time/min	63	63	64
ωtotal/CNY	53.28	65.86	75.56

Table 7

Optimal charging plan under model 2

Trip	Charging time after each trip/min
1	0
2	0
3	8
4	0
5	0
6	0
7	0
8	11
9	0
10	0
11	0
12	0
13	0
14	0
15	0
16	0
17	11
18	0
19	0
Nonoperation hours	34
The total charging time is 64 min and ω¯total=55.50 CNY

Table 8

Applicability of the optimal solution of model 1 for stochastic travel times

Simulation run	Probability that the constraint is satisfied
1	0.60
2	0.64
3	0.57
4	0.61
5	0.58
6	0.59
7	0.57
8	0.60
9	0.64
10	0.59

Table 9

Optimal charging plans with different trip numbers

Trip number	Operation hours	Model 1		Model 2		Relative changes(%)
		Optimal charging plan	ωtotal	Optimal charging plan	ω¯total
19	5:00–23:00	[0 0 8 0 0 0 0 9 0 0 0 0 0 0 0 0 20 20 0; 6]	53.28	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 11 0 0; 34]	55.50	2.22/4.17
18	5:00–22:05	[0 0 8 0 0 0 0 9 0 0 0 0 0 0 0 0 13 0; 30]	51.06	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 13 0; 29]	53.28	2.22/4.35
17	5:50–22:05	[0 6 0 0 0 0 8 0 0 0 0 0 0 0 0 5 0; 39]	48.84	[0 6 0 0 0 0 11 0 0 0 0 0 0 0 0 8 0; 34]	51.80	2.96/6.01
16	5:50–21:10	[0 6 0 0 0 0 8 0 0 0 0 0 0 0 0 0; 41]	46.62	[0 5 0 0 0 0 11 0 0 0 0 0 0 0 0 0; 41]	50.32	3.70/7.94
15	6:45–21:10	[0 0 0 0 0 12 0 0 0 0 0 0 0 0 0; 41]	48.10	[0 0 0 0 0 14 0 0 0 0 0 0 0 0 0; 40]	50.32	2.22/4.62
14	6:45–20:20	[0 0 0 0 0 10 0 0 0 0 0 0 0 0; 40]	44.40	[0 0 0 0 0 11 0 0 0 0 0 0 0 0; 40]	45.88	1.48/3.33
13	6:45–19:30	[0 0 0 0 0 7 0 0 0 0 0 0 0; 41]	40.70	[0 0 0 0 0 8 0 0 0 0 0 0 0; 41]	42.18	1.48/3.64
12	6:45–18:35	[0 0 0 0 0 5 0 0 0 0 0 0; 39]	36.26	[0 0 0 0 0 5 0 0 0 0 0 0; 41]	37.74	1.48/4.08
11	6:45–17:40	[0 0 0 0 0 0 0 0 0 0 0; 40]	29.60	[0 0 0 0 0 5 0 0 0 0 0; 36]	34.04	4.44/15.0

Notes:

In each cell of the columns “Optimal charging plan”, the numbers before semicolons refer to charging times at idle time in operation hours and the number after semicolons refers to the charging time in nonoperation hours. In each cell of the column “Relative changes”, the former value is the difference between ω¯total and ω_total, and the latter is calculated by (ω¯total-ω_total)/ω_total×100%

Table 10

Optimal charging plans with different battery capacities

B/kWh	Model 1		Model 2		Relativechanges(%)
	Optimal charging plan	ωtotal	Optimal charging plan	ω¯total
140	[0 0 8 0 0 0 0 14 0 0 0 0 0 0 0 0 18 17 0; 6]	56.98	[0 0 8 0 0 0 0 16 0 0 0 0 0 0 0 0 15 15 0; 10]	59.20	2.96/5.26
150	[0 0 8 0 0 0 0 12 0 0 0 0 0 0 0 0 5 14 0; 24]	55.50	[0 0 8 0 0 0 0 14 0 0 0 0 0 0 0 0 15 15 0; 12]	57.72	2.22/4.00
162	[0 0 8 0 0 0 0 9 0 0 0 0 0 0 0 0 20 20 0; 6]	53.28	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 11 0 0; 34]	55.50	2.22/4.17
170	[0 0 8 0 0 0 0 7 0 0 0 0 0 0 0 0 14 20 0; 14]	51.80	[0 0 8 0 0 0 0 9 0 0 0 0 0 0 0 0 15 15 0; 17]	54.02	2.22/4.29
180	[0 0 8 0 0 0 0 5 0 0 0 0 0 0 0 0 6 18 0; 26]	50.32	[0 0 8 0 0 0 0 6 0 0 0 0 0 0 0 0 20 20 0; 10]	51.80	1.48/2.94
190	[0 0 6 0 0 0 0 5 0 0 0 0 0 0 0 0 13 19 0; 20]	50.32	[0 0 8 0 0 0 0 5 0 0 0 0 0 0 0 0 12 12 0; 27]	51.06	0.74/1.47
200	[0 0 8 0 0 0 0 5 0 0 0 0 0 0 0 0 18 20 0; 11]	49.58	[0 0 8 0 0 0 0 5 0 0 0 0 0 0 0 0 15 15 0; 21]	51.06	1.48/2.99
210	[0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 6 14 0; 35]	46.62	[0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 18 18 0; 20]	47.36	0.741.59
220	[0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 12 19 0; 26]	46.62	[0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 18 18 0; 20]	47.36	0.74/1.59

Table 11.

Optimal charging plans under different price gaps

Electricity rate/CNY/(kWh)	Model 1		Model 2		Relative changes(%)
	Optimal charging plan	ωtotal	Optimal charging plan	ω¯total
[0.37,0.74,1.31]	[0 0 8 0 0 0 0 9 0 0 0 0 0 0 0 0 20 20 0; 6]	53.28	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 11 0 0; 34]	55.50	2.22/4.17
[0.37,0.64,1.11]	[0 0 8 0 0 0 0 9 0 0 0 0 0 0 0 0 5 15 0; 26]	51.48	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 19 19 0; 7]	53.30	1.82/3.54
[0.37,0.54,0.91]	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 16 20 0; 7]	49.62	[0 0 8 0 0 0 0 13 0 0 0 0 0 0 0 0 14 14 0; 15]	51.78	2.16/4.35
[0.37,0.44,0.71]	[0 0 8 0 0 0 0 16 0 0 0 0 0 0 0 0 14 17 0; 6]	47.38	[0 0 8 0 0 0 0 13 0 0 0 0 0 0 0 0 17 17 0; 8]	48.44	1.06/2.24
[0.37,0.47,0.57]	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 19 13 0; 11]	48.08	[0 0 8 0 0 0 0 11 0 0 0 0 0 0 0 0 19 19 0; 7]	49.56	1.48/3.08