A GRASP Approach for Solving Large-Scale Electric Bus Scheduling Problems

: Electrifying public bus transportation is a critical step in reaching net-zero goals. In this paper, the focus is on the problem of optimal scheduling of an electric bus (EB) ﬂeet to cover a public transport timetable. The problem is modelled using a mixed integer program (MIP) in which the charging time of an EB is pertinent to the battery’s state-of-charge level. To be able to solve large problem instances corresponding to real-world applications of the model, a metaheuristic approach is investigated. To be more precise, a greedy randomized adaptive search procedure (GRASP) algorithm is developed and its performance is evaluated against optimal solutions acquired using the MIP. The GRASP algorithm is used for case studies on several public transport systems having various properties and sizes. The analysis focuses on the relation between EB ranges (battery capacity) and required charging rates (in kW) on the size of the ﬂeet needed to cover a public transport timetable. The results of the conducted computational experiments indicate that an increase in infrastructure investment through high speed chargers can signiﬁcantly decrease the size of the necessary ﬂeets. The results also show that high speed chargers have a more signiﬁcant impact than an increase in battery sizes of the EBs.


Introduction
In recent years there has been a staggering increase in the use of electric vehicles (EVs) on a global scale.Many governments have provided significant incentives for EV adoption to meet net-zero emission goals.One part of EV adoption is connected to the use of motorised individual transport that represents nearly one-third of the global carbon emissions [1].A wide range of research is dedicated to optimizing their use through developing charging infrastructure and exploiting the potential of charging flexibility of EVs to solve problems related to the increased use of renewable energy production as discussed in [2][3][4].
The use of electric buses (EBs or eBuses) in public transport has been more prolific than personal transport vehicles.For example, in China there are currently more than 400 000 EVs in use.This corresponds to more than 14% of all buses in public transport in this country.An analysis on the expansion of EBs in China, through a case study regarding Beijing, can be found in [5].The use of EBs offers a multitude of advantages such as decreased air and noise pollution within metropolitan areas.Another reason is that EBs are essential to reach net-zero goals and the adoption of EBs offers an economically viable public transportation when compared to diesel buses [6][7][8][9].This is partly due to the high daily mileage of EBs.This results in a significant decrease in operational cost, due to the use of electricity instead of diesel, which quickly manages to cover the initial capital investment.
For the same reason, a high penetration of EVs can be observed within car hailing and taxi services.In parallel, there is a growing research interest on the efficient planning of charging infrastructure for large fleets and exploiting smart charging scheduling in such systems [10][11][12][13].
The use of EBs introduces new operational challenges related to limited driving range and long charging times in comparison to diesel buses (DBs).The use of DBs in urban public transport is well represented using the vehicle scheduling problem (VSP) [14,15].In the VSP, the goal is to optimize the assignment of a set of scheduled trips to a set of vehicles, where each trip is associated with one vehicle, based on some cost function frequently related to the number of used buses.This problem becomes significantly harder to solve in case of EBs due to additional constraints related to the range and the scheduling of the charging.These new issues are frequently modeled using the electric vehicle scheduling problem (E-VSP) [16] and its variations.It should be noted that even before the expanded use of EBs, the VSP has been considered with range constraints [17][18][19].A recent survey can be found in [20].We should note in passing that beyond vehicle scheduling there is a wealth of research dealing with public transport network design.Background regarding related models and research can be found in the following surveys [21,22].
Due to the growing amount of interest in the use of EBs, different approaches have been developed to address related practical problems.One part of this research is dedicated to charging infrastructure [23,24], while another direction is the analysis of operational issues.In [25,26], the long-term use of EBs is analyzed regarding scheduling and battery degradation.It should be noted that the infrastructural and operational aspects are closely connected and have been jointly considered [27,28] resulting in highly complex models.The scheduling of the use of EBs within a public transport system is well represented using the E-VSP and variations.Some of these variations are similar to the ones for the VSP, like the use of single [29][30][31][32] and multiple depots [33,34].This problem has been considered from the point of pure EB [29][30][31] and mixed fleets [32,33].A very interesting analysis of the problem of interest is given in [29], where the robustness is considered in relation to traffic conditions.In some variations, constraints related to the amount of charging that can be provided at a single time period have also been considered [29,32].The majority of these models have been solved using mixed integer programs [29][30][31][32].Only limited work has been done on the use of heuristic methods; some examples are the use of a combination of a greedy algorithm and simulated annealing [32] as well as the use of a genetic algorithm in [33].
The developed models and corresponding optimization methods are frequently used to evaluate the use of the EBs on relatively small systems.In the general case, only a few lines of a public transport system are analyzed, which is to a large extent a consequence of the high computational cost of solving the related MIPs.In this paper, the main objective is to provide a method that can be used to analyze the use of EBs on a large scale to cover a full public transport timetable.To achieve this, a greedy randomized adaptive search procedure (GRASP) [35] is developed for solving a variation of the single-depot E-VSP.In the proposed model, the charging time depends on the battery state.With the intention of lowering the complexity of the model, the following assumptions are made: charging is only allowed at the depot and there is no limit in the charging capacity.For the proposed model, firstly a MIP is developed to evaluate the performance of the GRASP algorithm.This is achieved with the following two objectives: (i) to observe from what system size on it is reasonable to use a metaheuristic method and (ii) to be able to evaluate the quality of solutions of the GRASP.
The developed GRASP algorithm is used for case studies on the use of EBs for public transport based on several real-world cases (cities).One part of the evaluation is dedicated to the battery size of EBs in the fleet.The second part is focused on the use of different chargers at the depot whose power can range from 50 kW to 600 kW.The conducted computational experiments are used to observe the size of the fleet needed to cover a public transport timetable for different battery sizes and charger types.
The main contributions of the conducted research are the following.Firstly, a new metaheuristic method is developed that extends the limited research on their application to eBus scheduling problems.The proposed GRASP algorithm manages to find high quality solutions for problem instances of larger size than previously developed methods.The second important contribution is related to the analysis of real-world public transport systems.To be more specific, the conducted case studies provide valuable insights on the use of eBuses with different ranges and charging infrastructure.In addition, the effect of different weather conditions, observed through different levels of energy consumption, on the variation in the size of eBus fleets needed to achieve a public transport timetable over the year is also analyzed.
The paper is organized as follows.The next section provides an outline of the model.The third section gives details of the graph formulation of the base problem.In the fourth section the MIP model is presented.Section 5 introduces the proposed GRASP algorithm.In Section 6, we present the results of the conducted computational experiments with a focus on some case studies for several cities.The paper ends with concluding remarks.

Model Outline
In this section, an outline of the proposed model is presented.The goal of the model is to find the minimal number of EBs that can cover all the bus lines in a timetable.More precisely, this is done by assigning a set of trips (timetable) to a set of EBs.A trip is specified with its origin, destination, duration, and starting time.The proposed model uses the following assumptions.The scheduling (assignment of EBs to trips) is done for a fixed time window, over a set of stations, and in a single depot.Moreover, the distances between all stations and the depot are assumed to be known in advance and recharging can only be done at the depot.
All trips are specified only by using origin and destination pairs and no intermediate stops are considered.It is assumed that all the buses start and end the day at the depot.An EB has a fixed driving range and a battery that must be recharged before the start of its schedule.Therefore, an EB can only be used if there is sufficient energy stored in its battery.It is assumed that all EBs move at a constant speed, so it is possible to use a scaling between length and duration.This is done with the goal of simplifying the calculation of the battery usage which will be proportional to the duration of a trip.
For an EB to be able to perform a trip from some origin, it has to be there before the starting time.The goal of the proposed model is to compare the number of needed EBs (having different ranges) and DBs to complete a public transport timetable.To that end, for a given timetable let us specify the specifics of the model as follows.

•
There are N locations and the corresponding set of locations is represented by set L = {1, ..., N}.

•
There is a single depot D ∈ L.

•
The distance between any two locations i, j ∈ L is known and is equal to The schedule is done over P time periods and a corresponding set of time periods is denoted by P = {1, ..., P}.

•
A timetable trip is a 4-tuple {o, d, s, l}, where origin-destination pairs are in the set of existing locations, that is, {o, d} ⊂ L and other parameters are part of time periods, that is {s, l} ⊂ P. The timetable is a set T of timetable trips.For simplicity of notation, we will use e = s + l for each trip to indicate the time of its completion.

•
An EB has a battery capacity that provides a fixed range R. It is assumed that all EBs have the same range.

•
The objective is to find the minimal number of EBs that can perform all the trips in the timetable.

•
An EB can only be charged at a depot.

•
An EB is fully charged at the first time period.

•
There is a charging rate γ, i.e., how much a battery is charged in one time period.

•
Charging of an EB is always done to full capacity.

Graph Formulation
In this section, the graph formulation for the scheduling of buses is given.The basic idea of the graph formulation is similar to the one used in [29].It only considers assignment of buses to trips in the timetable.In the later sections details on how this model is extended to include battery capacity and battery recharging are discussed more specifically.
Let us define a directed graph G(V, E) to represent the problem.For each trip t ∈ T there is a node t ∈ V.In addition, the node set V contains two auxiliary nodes D s and D e as the start (departure from the depot) and the end (arrival to the depot) trip of all buses, respectively.In the following text, the notation V t will be used for the set of all nodes that correspond to timetable trips.Let us define the distance between two nodes (trips) i and j as the distance between the destination i d of trip i and origin j o of trip j as follows.
An edge (i, j) is a part of the edge set E if it is possible to perform trip j after trip i.Let us formally define the set of edges E in the following way.
The set E n contains all the trip pairs (i, j) such that the starting time of trip j (s j ) is greater or equal to the sum of the completion time e i of trip i and the travel time from the destination of trip i to the origin of trip j (d ij ).An edge connects the start-trip node D s to all other nodes (E s ).All the trip nodes are connected to the end node D e .Finally, E is the set of all edges that exist in the graph.Note that, due to the definition of the edge set, the directed graph G does not contain any cycles.Each edge (i, j) ∈ E has an additional property equal to its length.
The schedule for buses is presented as a directed graph S(V, Ẽ) which is a subgraph of G.An edge (i, j) indicates that a bus performs trip j right after trip i.Let us make some observations about the solution defined in this way for the problem of interest.The set of vertices of graph S is the same as the one of G which means that all the trips are performed by some EB.Each trip node i ∈ V t will have exactly one incident edge for a node x ∈ V of type (i, x) since it can be performed directly before/after only one other trip.The number of edges having the form (D s , i) is equal to the number of buses needed to complete the timetable.Consequently, the schedule for a single bus will be a directed path having the form (D s , . . ., D e ).An illustration on how a trip timetable with known distances between locations and bus schedule is converted to the graph formulation can be seen in Figure 1.

Mathematical Model
In this section, based on the presented formulations, the mixed integer program for the problem of interest is presented.The proposed formulation has similarities to the ones given in [29,36].The model uses the following input parameters:

•
The distance between any two trips i, j ∈ V is known and is equal to d ij which is calculated as described in the previous section.

•
For each trip i there is a starting time s i and a duration d i .For simplicity of notation, let us define e i = s i + d i as the time when trip i is completed.

•
The charging rate per time period is denoted by γ.

•
Battery usage per time period is denoted by α.

•
The battery capacity of all EBs is denoted by G.
The schedule of EBs and battery-related constraints are specified in the model using the following decision variables:

•
For each (i, j) ∈ E let us define a binary variable x ij which states whether some EB performs a trip j after i.

•
For each trip i ∈ V let us define b i as the state of the battery of the EB when it has arrived at the origin of trip i but before performing it.

•
For each i ∈ V let us define a binary variable c i which states whether the EB that performs trip i charges its battery after completing it.
The selected objective function of the model is to minimize the number of used EBs to complete all the trips in the timetable.As previously stated, the number of buses is equal to the number of edges starting from the start node D s , which can be formally represented using the following formula: The last step in specifying the proposed model is defining its constraints.This is done in ( 9)- (20).The set of constraints can be divided into three groups.The first one, in ( 9)- (11), is dedicated to the completion of all the trips in the timetable by some EB.This set of constraints does not consider the range of EBs.The constraints ( 12)-( 18) focus on the state of the battery and consequently on the range of the EBs.The constraints ( 19)-( 20) are dedicated to the charging of EBs.
Parts of the battery-related constraints are conditional.For the sake of clarity, let us define two kinds of auxiliary variables, that are strictly specified by the values of x ji and c j .In practice they are not necessary but increase the readability of the constraints.The first one is as follows: The variable v ji is used to recognize the state when an EB performs trip j before trip i and the EB is not charged after trip j.In this case the value of v ji is equal to 0. The second type of auxiliary variables is as follows: The variable w ji is used to recognize the state when an EB performs trip j before trip i and it is charged after trip j.In this case the value of w ji is equal to 0. Now, all the constraints of the proposed model can be defined as follows: The constraints in Equations ( 9) and ( 10) guarantee that only one trip is performed by an EB after and before a timetable trip i, respectively.Note that the starting at the depot and finalizing the day at the depot are also considered as trips (pull-out and pull-in).The constraint in Equation (11) states that there is an equal number of EBs that leave the starting depot and arrive at the ending depot.
The next set of constraints is related to the battery state of charge.Constraints given in (12) provide the minimal value of the battery load before a timetable trip i is performed.The battery value must be greater or equal to the charge needed to perform trip i and to reach the depot afterwards.The constraints in ( 13) and ( 14) are conditional ones implemented using a large value M.These constraints only have an effect if v ji = 0, or, in other words, if an EB performs trip i after trip j and does not charge between these trips.In this case, the state of the battery is equal to the state of the battery after trip j (b j ) minus the battery used for completing trip j and relocating to the origin of trip i.In a similar way, the constraints given in ( 15) and ( 16) are conditional constraints for the case if a bus is charged between trips j and i, specified in the variable w ji .In this case the state of the battery charge before performing trip i is equal to the full battery (G) minus the charge needed to move from the depot to the origin of timetable trip i.It should be noted that a sufficiently high value of M for ( 13) -( 16) is 2G (two times the EB battery capacity).The constraint in (17) guarantees that an EB leaves the starting depot with a full battery.The constraints given in (18) are used to enforce the battery capacity before the start of any timetable trip.
The conditional constraints given in ( 19) are used to specify when the battery can be charged between trips j and i.These constraints are related to time periods.To be more precise, charging between trips i and j can only be performed if after completing trip j, at time period e j , there is enough time to move from the destination of trip j to the depot, move from the depot to the origin of trip i and charge the battery to its full capacity before the starting time of trip j (s j ).Note that a sufficient value for large M in this constraint is 2P.

GRASP Algorithm
In practical real-world cases, the scheduling of EBs to achieve a public transport timetable generally requires solving large problem instances.As it will be seen in the results section, the use of the presented MIP is effective only for problem instances of limited size.Hence, a metaheuristic approach-a GRASP algorithm-is developed to achieve high quality solutions for large-scale problem instances.In this section, details of the proposed algorithm are presented.First, a simple greedy algorithm is introduced.Second, it is extended to a GRASP algorithm by adding a local search and randomization.

Greedy Algorithm
The basic idea of the proposed algorithm is to start with a schedule containing a minimal number of trips and gradually expanding it with new ones.In the implementation of this simple idea, a graph formulation is used.The solution has the form of a set of paths P = {P 1 , . . .P n }, each corresponding to an EB (path), as previously described.Let us use the notation T(P) as the set of all timetable trips t such that t ∈ P. In addition let us define the same function for a set of paths (EBs) as A set of paths P is a complete solution if it contains all the trips in the timetable, formally written as T ⊂ T(P ).To fully specify the greedy algorithm, we need to define the initial solution, a list of candidates and the heuristic function.The simplest initial partial solution can consist of a single path P 0 = {D s , D e }.On the other hand, it is evident that there is a set U of trips that cannot be performed after any timetable trip, or, in other words, must be the first daily trip of an EB, which corresponds to the nodes that have no entering edges in the set E n but only in E s .Using the set U, a larger initial solution, that has a path for each such trip, is defined as Note that for a partial solution P, a trip t can only be added to at most one position in each path P i ∈ P.This is due to the fact that the starting/ending times of trips in P are strictly ascending.Consequently, a candidate for expansion will have the form (t, P) where t ∈ T \ T(P ) is a trip, not already in the partial solution, and P ∈ P is a path (EB) in the partial solution.Note that it is possible to have a trip t that cannot be inserted to any of the paths P ∈ P in the current partial solution.Because of this fact, the set of paths, used to generate the candidates for the greedy algorithm, has an additional one that makes it possible to add new buses to the solution.This path will have the form (D s , D e ).
Let us analyze the valid candidates that are in the form of (t, P).The term "valid" is used in the sense that the battery charge level of an EB has to be positive and less than the battery size during the entire schedule.To achieve this, let us define a function B(P, i) as the state of the battery for a path P after i trips have been completed.Let us use the notation p i for the i-th timetable trip in path P. The function B(P, i) is defined recursively as follows The equation given in (23) states that the initial (bus at depot) level of charge of the battery is equal to its capacity.Equation ( 24) is used to calculate the battery state after trip i.The notation and ⊥ is used for the values true and false, respectively.It is assumed that if it is possible to charge the EB between trips p i−1 and p i it will be charged.In this case the battery state after trip p i is equal to the full battery minus the battery usage needed to travel from the depot to the origin of trip p i and performing the trip p i .In case the battery cannot be charged, this value is equal to the battery state after trip p i−1 minus the battery used for moving from the destination of trip p i−1 to the origin of trip p i and performing the trip p i .Let us specify when the EB can be charged between p i−1 and p i , or, in other words, the function CanCharge(P, i).This is done in the same way as in constraint (19) of the MIP formulation.For the sake of increasing the readability, the definition of the corresponding function is explicitly given in the following way: The equation given in (25) defines the function that gives the time needed to charge the battery after the i-th trip of path P is completed and the bus has reached the depot.The charging time is inversely proportional to the charging rate γ.The equation given in (26) states that an EB can be charged between trips p i−1 and p i if the starting time of trip p i (s p i ) is after the earliest time trip p i−1 can be completed (e p i−1 ) and the time the EB needs to move to the depot, and fully charge the battery and move from the depot to the origin of trip p i .
Using the function B(P, i), a Boolean function valid(P) can be defined that states if a path is valid as follows Equation (27) states that for each trip in a path, the charge of the battery is sufficient to arrive to the depot.Now, let us use the notation P +t as the path acquired by inserting trip t in path P. A candidate (t, P) is valid if valid(P +t ) is true.
The next step in defining the greedy algorithm is specifying the heuristic function.The idea of the heuristic function is to first avoid adding new paths (EBs) to the solution, unless necessary.The second part is to avoid unnecessary movements of EBs.Let us first define the length of additional movement of an EB, or more precisely for its path P, using the following formula: In ( 28), the function Move(P, i) is equal to the additional movement performed between trips i − 1 and i of path P. The additional movement depends on the case if charging is performed or not between these two trips.The total additional movement for a path P is given in (29) with the function Move(P) and is equal to the sum of additional movements for all consecutive trips.
Using function Move, the heuristic function h, for a valid expansion candidate (P, t), can be specified as follows The heuristic function h(P, t), given in (30), differentiates between the cases when a new path is added or not.In case when no new path is added, the preference is for the smallest increase in the additional movement.In case when a new path is added, the heuristic function h n (t) is used and only depends on the trip t.In Equation (31), m and M are constants satisfying m << Move(P, t) << M for any candidate (P, t).The heuristic h n (t) states that in case a trip t can be added to some path Q ∈ P producing a valid path Q +t , adding a new path (D s , t, D e ) to the solution is the least desirable.On the contrary, if the trip t cannot be added to any path Q ∈ P producing a valid one, then adding the path (D s , t, D e ) is most desirable.The reason for this is that if a path Q +t is not valid, it is not possible to expand Q to a path R in a way that R +t is a valid path.Therefore, it is preferable to add the new necessary path earlier to the solution.This is because in the later iterations of the greedy algorithm there will be a larger number of candidates for expansion, potentially of higher quality.
Using the heuristic function h, the greedy algorithm can be fully specified for which details can be seen in the pseudocode given in Algorithm 1.In this algorithm, the first stage is the generation of an initial solution P based on Equation (22).Next, the set of available trips T c is set to all the trips in the timetable T .In the main loop, at each iteration, firstly, the set of all valid candidates for expansion of the partial solution C is generated.Elements of C have the form (P, t) where P ∈ P and t ∈ T c .Out of all (P, t) ∈ C, the one having the minimal value of the heuristic function is selected for expanding the partial solution.Ties are randomly broken.This means that the timetable trip t is inserted in the path P. Further, notice that t can only be inserted at one unique place in path P based on its starting and ending times s t and e t , respectively.Next, the trip t is removed from the set of available trips T c .The main loop is completed when all timetable trips are added to the solution P, or, in other words, T c is an empty set.

Algorithm 1 Greedy Algorithm
Generate initial solution P using Equation ( 22) T c = T \ T(P ) while T c = ∅ do Generate set of valid candidates C for P, T c (P, t) = argmin (P,t)∈C h(P, t) Expand P with candidate (P, t) T c = T c \ {t} end while

GRASP Extension
In this section, details are provided on how the presented greedy algorithm is extended to the GRASP metaheuristic.To achieve this, the greedy algorithm needs to be randomized and a local search is added.In the proposed algorithm, the randomization is done using the standard restricted candidate list (RCL) approach.In case of an RCL with n elements, instead of selecting the candidate for expansion having the best value of the heuristic function, as in the deterministic greedy algorithm, one of the best n candidates is randomly selected.
The local search follows an approach commonly used in bin packing problems [37].The idea is to check if the schedule has an EB such that all of its trips can be reassigned to other already used EBs.If this is possible, the procedure is repeated for the new schedule.The application of this simple concept is given in the following text.Firstly, the procedure on how the trips performed by one EB are relocated to the other ones is explained and later the complete local search.
The goal is to reassign one by one all the trips from one EB to the remaining ones.The pseudocode corresponding to the procedure of attempting the relocation of all the trips from the EB corresponding to path P to all the remaining EBs, a subset of the set of paths P, is given in Algorithm 2. The function RelocatePath returns a new set of paths Q if all the trips in path P can be relocated or an empty set if this is not possible.

Algorithm 2 Path Removal
function RELOCATEPATH(P, P) In this function, firstly, Q is set to the set of all paths in P except P. The set Q is used for creating the new schedule.Next, a set of trips R is set to all the trips that need to be relocated.The goal of the main loop is to relocate all the trips in R to other EBs one by one.At each iteration of the main loop, one of the trips r ∈ R is selected randomly.For trip r, the set of valid candidates (EBs) for reassignment Q v is generated using the function valid(Q +r ) for all the EBs Q ∈ Q .In case Q v is an empty set, the function RelocatePath cannot reassign all the trips in path P, and the function returns an empty schedule.
Note that by randomizing the order in which the trips in the original path P are selected for reassignment increases the number of potential schedules that are explored in case the local search is applied to a same solution P. The next step is selecting to which EB the trip r is reassigned from the set of all the valid ones Q v .This is done using the same heuristic h as in the case of the greedy algorithm.To be more precise, the trip r is reassigned to the EB, corresponding to the path Q ∈ Q v , that has the minimal value of h(Q, r).The final step in the main loop is removing r from the set of trips R that need to be reassigned.The main loop is completed when all the trips have been reassigned or, in other words, when set R is empty.
Let us now define the complete local search using the function RelocatePath.The details of the local search can be seen in Algorithm 3. In this algorithm, the first step is setting the current best solution P c to the schedule P. The outer loop attempts to improve P c when possible.This is achieved by setting the set of candidate paths for removal Q to all the paths in P c .The inner loop checks if any of the paths in Q can be removed from P c .This is done by selecting a random path Q ∈ Q. Next, the function Relocate(P c , Q) checks if it is possible to remove the path Q.If this is true, the current best solution is updated, and the inner loop is exited.In case this is not true, the path Q is removed from the set of candidate paths Q.This procedure is repeated until set Q is empty, or, in other words, it is not possible to further improve the solution P c .

Algorithm 3 Local Search
function LOCALSEARCH(P)

Implementation
In this subsection, the details of the proposed GRASP algorithm and its pseudocode are discussed and details are presented in Algorithm 4.

Algorithm 4 GRASP algorithm
N = 0 while (N ≤ MaxSolutions) do Generate solution P using GreedyRCL(n) P = LocalSearch(P ) The GRASP algorithm repeatedly generates new solutions P using function GreedyRCL(n).This is done using the randomized version of the greedy algorithm given in Algorithm 1.In this algorithm, the heuristic function h(P, t) is substituted with an RCL having n elements.On each such solution, the local search is applied and is checked if it is the new best solution.This procedure is repeated until a maximal number of solutions is generated or potentially some other stopping criterion is satisfied.
It is important to point out that in the practical implementation, it is not necessary to recalculate all the values of the heuristic function at every iteration.For instance, if at one iteration, it is impossible to add a trip to an EB it will remain to be so in all the later steps.The values of functions h, CanCharge, B only need to be recalculated at iteration i + 1 for the paths that have been changed at iteration i. Exploiting this can significantly improve the computation performance of the method.

Results
In this section, the results of the conducted computational experiments are presented.The main goal is to provide a tool that can analyze the potential increase in the number of EBs needed to achieve a public transport timetable compared to the use of DBs.Because of this, the computational experiments have two objectives.The first one is to evaluate the effectiveness of the proposed MIP and GRASP for the newly proposed problem.The limits on the size of the problem instances that can be solved to optimality using the MIP within a reasonable time are found.The obtained solutions are, then, used to evaluate the performance of the proposed GRASP algorithm with respect to the quality of results.In addition, the computational cost of the GRASP algorithm is also analyzed.These tests have been conducted on synthetic problem instances.
In the final experiments, the GRASP algorithm is used to address practical issues related to the use of EBs in public transport.This is done through case studies on several real-world cities and corresponding public transport timetables.Moreover, the relation between EBs' battery size and charging rates is analyzed.
The MIP has been implemented in ILOG CPLEX for C#.NET through Concert Technology.The implementation of the GRASP algorithm has been done in C# in Visual Studio 2019.The computational experiments have been performed on a personal computer running Windows 10 having Intel(R)Xeon(R) Gold 6244 CPU @3.60 GHz with 128 GB memory.

Synthetic Data Experiments
In this subsection, the experiments conducted on synthetic data are presented.A large number of problem instances is generated having a wide range of sizes.In this way, it is possible to have an extensive evaluation of the proposed GRASP and MIP approaches.
Let us start with providing details of the methods for generating such instances.The first step is randomly selecting N points in a square (0, α) × (0, α).Each point represents a location in the problem with one being selected as the depot.The distance between any two locations a and b is equal to the Euclidean distance between the corresponding points.The next stage is generating random bus routes.The scheduling is conducted over a time window of one day, and a period corresponds to one minute, or, equivalently, the problem is solved over 1440 periods.A bus route is defined by its origin l o and its destination l d , the start time t s and the end time t e , the duration d of a trip, and the frequency f of departure times.Note that the duration of a trip with some origin and destination locations is not directly related to the distance between them.The reason for this is that bus lines generally do not have a route that corresponds to the shortest path between the origin and the destination but instead try to satisfy the transport needs of passengers by visiting points of interest.
When generating the bus lines, the following rules have been used.A single location is randomly selected as the depot.Every other location l ∈ L is an origin or destination of at least one bus line.The starting time, duration, and frequency of all bus lines are randomly selected from an interval (s min , s max ), (d min , d max ), and ( f min , f max ), respectively.The ending time for each line is defined using its relation to the start time.To be more precise, it is equal to t s + δ where δ is randomly selected from (δ min , δ max ).Individual trips for each line are generated in a natural way, having the form (l o , l d , t s + i f , d) where an integer i ≥ 0 is such that t s + i f is less than the ending time.
To evaluate the MIP and the GRASP algorithms, the size of the problem instances has been related to the total number of trips R and varies from 20 to 2000.The number of locations depends on this value and is equal to N = 1 2 R.An overview of the values of the parameters used for generating the synthetic problem instances can be seen in Table 1.The setting for the comparison of the MIP and the GRASP is as follows.The MIP has been solved in CPLEX with the default parameter setting with a time limit of 1200 s for each instance.In case of the GRASP, a total of 5000 iterations is performed, or, in other words, 5000 solutions are generated.The size of the RCL is 2. The aggregate results for each problem size are given in Table 2.The averages of solutions and computational times are given.In addition, the average values of the lower bounds found by CPLEX and the number of found optimal solutions are also provided.The first observation is that the MIP can find the proven optimal solutions for small problem instances having 20 to 40 trips, but not for all of them.The MIP could find feasible solutions and lower bounds for problem instances having up to 200 trips within the time limit of 1200 s.The GRASP algorithm managed to find 20 out of the 23 proven optimal solutions.In case of the problem sizes of 20 and 30 trips, for which the MIP found all the optimal solutions the average error of the GRASP is 0 and around 5%, respectively.The performance of the GRASP algorithm, evaluated by average solution quality, becomes better than the MIP starting from instances with 40 trips on.The difference between the average solution quality and the average lower bounds found by the MIP, in case of problem sizes where the MIP did not have a good performance, is significant and close to 50%.It is expected that the lower bounds provided by CPLEX are not tight.
The main reason for developing a metaheuristic for this problem is to provide a tool that can generate high quality feasible solutions for large problem instances within a reasonable time.The proposed GRASP algorithm achieves this by solving the instances with 2000 trips within less than an hour.In Table 2, the average computational times are included to provide insights to the scaling of the method.The computational time grows 100 times for instances of 100 and 200 trips to 1000 and 2000, respectively.This means that, in practice, the computational cost has a growth which is close to being quadratic.
With the goal of having a better understanding of the practical performance of the proposed method, in the sense of convergence speed, time-to-target plots (ttt-plots) are used to visualise runtime distributions.Such plots provide information about the probability that a method will find a solution at least as good as a given target value for a given problem instance within a given running time.This type of information is very useful for characterizing the running times of stochastic algorithms for combinatorial optimization problems.An extensive description of ttt-plots and tools for their generation can be found in [38,39].
The experimental setup for generating the ttt-plots is the following.Three different problems sizes having 400, 600, and 800 trips are evaluated.For each of them, one characteristic instance is selected and the proposed GRASP algorithm is run 75 times using different seeds for the random number generator.The number of generated solutions for each run is 10 000.A large number of generated solutions is used with the intention of having all the runs reach stagnation.The target value is the solution found with the minimal number of used eBuses.It should be noted that, in case of problem sizes 400 and 600, all the runs found the same quality solution.In case of the largest size, less than 5% of the runs did not find this solution and are excluded from generating the graph.The related ttt-plots can be seen in Figures 2-4.From them, it can be observed that the computational cost for the GRASP algorithm, to find the solutions of the targeted quality, grows faster than quadratic.In case of instances having 400, 600, and 800 trips the time needed to have 80% of the runs find the solution of the desired quality is close to 25, 50, and 370 s.This is a common behavior for GRASP methods for large-scale problem instances, since it lacks a learning mechanism, and the relatively random search of the solution space is not very efficient for instances larger than some problem size.In the final part of the analysis of the performance of the proposed GRASP algorithm, we make a comparison to other approaches for scheduling charging of eBuses.It should be noted that this is not a direct comparison, since the existing models focus on different properties of such system.In addition some of the methods find proven optimal solutions and others only near optimal ones.As previously stated, the goal of the proposed method is to find high quality solutions for large problem instances which correspond to real-world public transport systems.Because of this the comparison focuses on the size of instances such methods can solve as well as the computational cost.In [29], a more complex model is proposed, which considers robust scheduling in a stochastic traffic environment.In this research, a branch-and-price algorithm is used to find approximate solutions to problem instances of at most 96 trips with a computational cost of 117 s.In the work by Awesabi et al. [28], a wireless setting of single-depot dynamic charging systems is considered and the system is optimized using a MIP approach.This model also considers the optimization of the battery size of the used eBuses.The proposed method optimizes a small off-campus college transportation system having six routes.A cost optimization of a public transport system based on eBuses is provided in [30].In this work the proposed MIP is solved using a combination of column generation and Lagrangian relaxation for a system having four routes and 543 trips with a computational cost of 539 s.To the best of the authors knowledge, the largest problem instances for eBus scheduling have been solved using a column generation-based approach combined with a local search [36].In this work, problem instances having 941 trips are solved within a time period of 46 to 58 hours.It can be observed that the proposed GRASP algorithm manages to solve much larger instances at a significantly lower computational cost.To be more precise, the conducted computational experiments indicate that an average of less than 3000 s is needed to find high quality solutions for instances having 2000 trips.This increased computational efficiency makes it possible to conduct an extensive analysis of real-world systems as it will be seen in the following case studies.

Case Studies
In this section, the results of case studies on real-world public transport systems are presented.The computational experiments are conducted using the previously presented GRASP algorithm.The procedure for generating test instances for the selected public transport systems is as follows.The used timetables for each case study are obtained from official sources available online.The used data for bus lines include the origin and destination locations and the departure and arrival times for each of the trips performed by a bus.The departure and arrival times are mapped to time periods in the model in which a single time period had the length of one minute.Note that, some bus lines are circular and some are bi-directional.In case of bi-directional lines, a trip would correspond to one direction of travel.The distance between locations is calculated based on Googlemaps using the Google Cloud Services.Traffic conditions have not been considered in the model.
The tests are used to evaluate the effect of different battery capacities and charging powers on the number of EBs needed to achieve a public transport timetable.The battery capacity of EBs commonly used for this task is from 200 kWh to 425 kWh [40].The charging power used for EBs can greatly vary based on the installed selected infrastructure [41].In the conducted computational experiments, a charging range between 50 kW and 300 kW is considered.Contrary to DBs, the range of EBs has a much larger level of variety [42,43] which especially depends on driver habits and weather conditions.Empirical studies have shown that the consumption for 12 foot EBs averages around 1.4 kWh/km, but can range from 1 kWh/km in case of an experienced driver in ideal conditions to even 2.35 kWh/km in extremely cold weather when electric heating is used [44].In the conducted case studies the average energy consumption and the two extreme values are evaluated.In the case studies, the speed of EBs is considered constant and equal to 20km/h [45,46].The cities/counties chosen for the case studies have been selected in a way to represent public transport systems having a wide range of different properties.In the following text the relevant details are provided.
The Ocala/Marion County, Florida, USA system consists of seven bus lines.Most of the lines are bi-directional and have departures every hour in each direction.The trip in each direction has a duration of around 20-30 min.All the bus lines start their trips from a central terminal.As a second case study, the Hernando, Florida, USA system with only four bidirectional lines is chosen.The line's length and departure frequencies are similar to the case of Ocala/Marion County.The system does not have a main terminal; instead, the lines connect different locations in the city.
The third case study is the Collier, Florida, USA system with 20 lines, half of which are circular bus lines.This bus network has one main terminal and around half of the lines use this terminal as an arrival or departure station.The frequencies vary from 60 to 160 min, with the majority of lines having departures every 90 min.The duration of lines ranges from 35-80 min, with the majority being 40 min.The fourth case study, the Wilmington, North Carolina, USA system has close to 20 circular lines.The vast majority of the lines have an hourly frequency and a duration of close to an hour.The system has two main terminals from which most of the buses depart.
The fifth case is the Colorado Springs, Colorado, USA system which has around 30, mostly bi-directional, bus lines.The majority of bus lines have frequencies of 15, 30, and 60 min.The duration of most of the trips would be 20-30 min.The system has four main stations and the majority of the lines have a departure station at one of them.As a sixth case study, the Doha, Qatar system with 50 bus lines is considered.Nearly half of the bus lines are circular.The majority of the lines have 20 and 30 minute frequencies.There are also two large terminals from which most of the buses depart.In addition there is a group of lines that depart from the airport.One characteristic of the Doha system is that the lines have a long duration ranging from 60 to even 180 min.
The locations of the bus depots for all the tested cities are not available.On the other hand, currently all these transport systems are operated using DBs so there is a potential for changing the location of the EB depot due to power distribution requirements for the installation of chargers.In the generated test instances, the locations of the depots are selected by choosing empty locations that are close to the largest terminal in the system.
It should be noted that, in case of Doha and the other systems having major terminals, the arrival times and departure times are set in a way to maximize the possibility of passengers changing lines.This results in the corresponding problem instances becoming much harder to solve by the MIP compared to the synthetic ones.This is due to an increased level of symmetries in the solutions space, which does not notably effect the GRASP algorithm.
The objective of the conducted case studies is to evaluate the number of EBs needed to cover a public transport timetable.This is done through a comparison of the number of DBs needed to achieve the same task.DBs, in general, have ranges that are sufficient to perform all the daily trips on a full tank.In practice this means that in this case the range constraints are not relevant.Consequently, the number of DBs needed to achieve a timetable can be acquired using a MIP having the objective function given in Equation ( 6) and only the constraints in ( 9)- (11).This model is significantly less computationally expensive than the one for EBs.These simplified problem instances generated for the systems used in the case studies can all be solved to optimality using CPLEX within less than 10 min.
Before presenting results of the conducted computational experiments, let us point out the following.The GRASP algorithm finds high quality solutions for EB-based systems which are compared to the optimal ones found by the MIP for DBs.Consequently, the provided results are used to analyze trends and not to provide exact configurations of the systems.In case of using the proposed GRASP algorithm to practical real-world systems, with the intention of optimizing them, the method has a significant drawback of not being able to evaluate the precise distance from an optimal solution.Because of this, for such applications, it would be essential to conduct a more in-depth error analysis through finding tight lower bounds.It should be noted that this is a simplification of the real-world problem that does not consider driver-related constraints.Such constraints in many cases also impact the number of buses needed to perform all the trips in a public transport timetable.
The summary of the results is presented in Figures 5-7 for different system sizes.For each city, a subfigure is used to show the number of EBs needed to perform a timetable for different combinations of EB battery capacities, charging powers at the depot, and energy consumption per unit distance.The first thing that can be observed is that, in case of ideal conditions with 1 kWh/km energy consumption, the increase in the needed EBs compared to DBs, is generally relatively small for small and medium systems.In case of large systems, the increase is more significant and can reach even 30%.It should be noted, that GRASP algorithms generally have a worse performance for large problem instances and it is expected that part of the increase is related to this fact.In case of all tested systems, it is noticeable that the use of higher charging powers has a more significant impact than the battery capacities.An increase in charging power at the depot decreases the time needed to charge an EB.The impact of increasing charging power is the most significant for low battery capacity EBs, but is significant for all of them.The difference in the number of EBs needed in case of average energy consumption compared to the extreme case is very significant.Once more this increase is most notable in case of low charging powers and battery capacities.It should be understood that although the extreme energy consumption events are rare, it is expected that they will occur annually during periods of extreme weather.The results indicate that it can be extremely beneficial adding alternative energy sources for heating and cooling.Another observation that can be made is that, in case of larger initial investments, longer range EBs, and higher power chargers, the fluctuation in the number of needed EBs during the year can be greatly reduced.

Conclusions
In this paper, the practical of using electric buses in public transport systems has been analyzed.The addressed problem consists of scheduling EBs to perform all the trips in a public transport timetable.This has been done by firstly developing a mathematical model for a single-depot system in which EB charging times depend on the current state of the battery.Due to the fact that the developed MIP model can only solve small-scale problem instances, that do not reflect real-world systems, a metaheuristic approach has also been developed.To be more precise, a GRASP algorithm has been created that can solve problem instances that correspond to urban public transport systems.The effectiveness of the GRASP algorithm has been evaluated through a comparison with optimal solutions acquired using the MIP for small instances.The conducted computational experiments have shown that the GRASP manages to find many optimal or near optimal solutions for those instances.The computational cost of the GRASP algorithm is significantly lower than of the MIP.It manages to find better solutions than the MIP, even for relatively small instances, where the MIP terminates early and can find solutions for problem sizes on which the MIP cannot be applied.
The developed GRASP algorithm has been used for case studies on real-world public transport systems.The related computational experiments have shown that the number of EBs needed to achieve a public transport timetable can be significantly higher than DBs.The developed model indicates that by using high speed chargers the number of needed EBs can be greatly reduced and using this type of chargers is more relevant than the range of the EBs.In addition, systems using higher charging powers and having EB fleets with larger battery capacities have a much lower seasonal fluctuation in the number of needed EBs.The results of the conducted case studies provide valuable insights for the cost analysis of such systems.
In the future there are several directions in which the presented research can be extended.The first one is extending the model to include multiple-charging depots, limiting the charging capacities, considering driver-related constraints and different types of disturbances in a stochastic model.It should be noted that the used solution method, the GRASP, is very suitable for parallelization.This type of speedup would be highly relevant for such extended models.Another direction is adapting the proposed GRASP algorithm to previously developed mathematical models related to the use of EBs.This can enable the application of such model to much larger problem instances that well represent real-world systems.Another issue relates to disturbances in public transport (see [47] for a survey and [48] for some conceptional work) as they can even arise with special features related to EBs (see, e.g., [49]).Finally, the Fixed Set Search metaheuristic [50], which adds a learning mechanism to GRASP, can be applied to try to increase the quality of found solutions for the problem of interest.

Figure 1 .
Figure 1.Illustration of the conversion of a timetable with known distances between origins and destinations between trips.The directed graph (middle) has a node for each trip and edges only connect trips that can be sequentially done.The set of red edges (bottom) is a schedule of EBs.It contains the two buses with schedules (Depot, a, c, Depot) and (Depot, b, d, Depot).

Figure 2 .
Figure 2. Time-to-target plot for an instance with 400 trips based on 75 independent runs of the GRASP algorithm.

Figure 3 .
Figure 3. Time-to-target plot for an instance with 600 trips based on 75 independent runs of the GRASP algorithm.

Figure 4 .
Figure 4. Time-to-target plot for an instance with 800 trips based on 75 independent runs of the GRASP algorithm.

Figure 5 .
Figure 5. Number of EBs with different ranges dependent on depot charging power to achieve a small public bus transport system.In each of the subfigures the dashed line corresponds to the number of DBs needed for the same task.

(a) 1 Figure 6 .
Figure 6.Number of EBs with different ranges dependent on depot charging power to achieve a medium public bus transport system.In each of the subfigures the dashed line corresponds to the number of DBs needed for the same task.

(a) 1 Figure 7 .
Figure 7. Number of EBs with different ranges dependent on depot charging power to achieve a large public bus transport system.In each of the subfigures the dashed line corresponds to the number of DBs needed for the same task.

Table 1 .
Parameters used to generate the synthetic problem instances.

Table 2 .
Comparison of the proposed GRASP and MIP for synthetic problem instances of different sizes.