Battery Electric Bus Network: Efﬁcient Design and Cost Comparison of Different Powertrains

: Recent advances in the electromobility of bus ﬂeets improve transit service sustainability but require the development of vehicle recharging facilities. The paper deﬁnes a methodology to design an efﬁcient transit network operated by battery electric buses in cities with grid-shaped road network, based on continuous approximations. An analytical model deﬁnes the optimal network conﬁguration that minimizes the agency cost, the monetization of emissions and the travel time of transit users. The analytical model allows the comparison of total cost, emissions and bus performance to other fuel powertrains. The methodology is tested in Guadalajara (Mexico) to propose an alternative bus conﬁguration, outperforming the current bus service and reducing the agency cost and environmental impact. The analytical model justiﬁes the network rationalization in fewer routes to reduce the total cost of the system. The deployment of standard battery electric buses with opportunity charging scheme obtains the lowest total cost of the system.


Introduction
This paper aims to design an optimal bus network configuration operated by full electric powertrain, estimating the total number of electric buses, charging resources and the total cost incurred by transit agency and users. The proposed methodology allows a comparative analysis of the resulting bus design with internal combustion engine bus powertrains. The paper also checks the potential benefits and cost savings of the optimized bus network design obtained by diesel powertrains with the current bus network configuration in Guadalajara, Mexico. A sensitivity analysis is also performed with regard to the unit cost parameters in order to estimate under which domains of distance cost, resource depreciation cost, emission monetization and charging infrastructure, the batteryelectric buses are more cost-efficient than conventional diesel buses.
Mass transit systems must provide a competitive performance to achieve a positive modal redistribution from private cars. To that end, several cities have undergone a significant improvement in the transit service, together with access restrictions to private vehicles, land use planning and redensification strategies. A crucial issue is to provide an efficient network with an outstanding performance at a limited operating cost. [1]. The optimal balance is achieved by the efficient provision of the spatial configuration of routes and the temporal coverage of services (headways). Therefore, the transit network design problem is aimed at minimizing the total cost incurred by users and transit agencies. This problem has received the attention of many researchers, proposing discrete or continuousbased optimization models.
Discrete-based optimization models need the construction of a graph to resemble the existing road configuration with which the bus network may operate. Given an origin-destination matrix (O-D), different mathematical optimization techniques have been proposed to solve the network design problem. This modeling approach allows representing with a good detail the street configuration of the city and the demand distribution. Nevertheless, as this problem is known to be Np-Hard in computational time [2], mathematical models can only obtain exact solutions in small problems [3]. For medium and large networks, researchers proposed heuristic algorithms [4][5][6][7][8][9] to obtain a suboptimal solution in affordable computational time. Later, metaheuristics algorithms [10][11][12][13][14] have been developed to fine tune a given suboptimal solution.
On the other hand, continuous or parsimonious models propose ideal conceptualizations of transit network whose physical layout and temporal performance are controlled by a small set of spatial and temporal decision variables. They estimate the network performance in terms of user travel time and operating costs. Each system cost component incurred by users (in-vehicle, waiting and access times) or transit agency is formulated as a function of decision variables.
A set of compact formulas is developed for each cost term using geometric probability. Then, the combination of optimal decision variable values that minimize the total cost is obtained by a continuous optimization, following a grid search procedure. This set of optimal values allows the creation of the ideal skeleton of routes and stops as well as the frequency scheme at which the system needs to operate. Later, decision makers create the final map of the routes based on the available street network, resembling the ideal scheme as much as possible. The main advantage of these models is the low computational effort to obtain the solution and the wide transferability to numerous case instances where the city shape is similar to the street network pattern considered in the model.
Nevertheless, assumptions such as the uniform demand distributed must be considered to simplify the calculations in the models. Several contributions in this field were developed for grid [15], radial network configurations [16]. Daganzo [17] proposed hybrid configurations to serve square-shaped cities. In this approach, the potentialities of both grid and radial schemes were combined in the same network pattern. This scheme was extended to a hybrid concept to cities where the road network follows a circular and ring mesh of streets [18,19]. Other works have been done comparing the performance of door-to-door and transfer-based grid and hybrid networks, stating that the latter provides lower cost designs to serve the mobility of the city when demand is consolidated in the city center [20].
However, few authors have integrated the environmental impacts derived from the bus service in the network optimization. Chester et al. and Griswold [21,22] consider the greenhouse gas (GHG) contribution of the transit vehicle manufacturing, road, station construction and the fuel consumption in the network design problem. These emissions are integrated into the objective function and the new optimal values of network decision variables are evaluated for different transport modes (bus, BRT, LRT and metro). This approach has also been analyzed in Cheng et al. and Griswold et al. [23,24] in San Francisco and in the new Barcelona bus network proposed in Estrada et al. [25] respectively. In Griswold et al. [26], travel time elasticities are considered in transit and private car systems to evaluate the modal split changes in scenarios, where the total amount of GHG emissions for the transit operator or the total system is constrained.
In Cheng et al. [27], the optimal design and emissions of a more sophisticated transit system, composed by trunk and feeder services, is evaluated under an elastic demand assumption. It implies a competition between the private car and Internal Combustion Engine (ICE) transit services. Lee and Madanat [28] proposed a methodology to define the optimal network of the charging infrastructure to minimize the emissions in a given city, promoting full electric vehicle operation. Nevertheless, the proposed charging network configuration is developed for passenger cars moving in different areas of the city, without any link with existing bus corridors. The location of the charging facility is not incorporated to the transit network design problem.
One of the crucial issues in the former references addressing the transit network design was the estimation of the amount of GHG emissions generated in the operation, construction and installation of transit facilities and auxiliary resources. The emissions Sustainability 2021, 13, 4745 3 of 28 generated during the vehicle operation due to fuel combustion (exhaust emissions), tire rolling, braking or fuel evaporation have been fully analyzed in the EEA model and Argone GREET model [29,30]. Moreover, several studies have conducted life cycle cost analysis (LCCA) of the batteries and vehicles including energy and environmental aspects, and real energy consumption of electric buses and its benefits in urban areas [31][32][33][34][35][36][37][38].
Nevertheless, the LCCA of elements not related to the bus operations (lane infrastructure and charging infrastructure) has been barely addressed in the literature. These parameters have also been complemented by the energy consumption and emission data gathered from demonstrations conducted in the ZeEUS project, one of the main studies devoted to test electric bus performance on real conditions [39]. Unfortunately, the vast majority of LCCA analysis highlighting the lower emission contribution of electric rolling resources, do not consider the extra number of resources that electric bus operations will cause. Therefore, the main motivation of this paper is to fill this gap. Hence, in this work an estimation for the number of electric resources and capital and operating cost that transit agencies will incur. Section 2 describes the methodology and row data used to estimate the agency, user and environmental costs developed as a function of bus service design variables. These formulas are merged into an optimization model to identify the best bus network operation under different scenarios. Section 3 shows the result of the model presenting the optimal parameters for the design of a door to door orthogonal network. Additionally, a sensitive analysis is representing integrating the change of variables such as demand, city size and cost of infrastructure and electric vehicles. On the other hand, Section 4 indicates the discussion of the results. Finally, general conclusions, Guadalajara achievements and most important remarks are drawn in Section 5.

Materials and Methods
This section will describe the methodology followed and all the assumptions adopted in the model.

Methodology
It was assumed a rectangular shaped city of area D x ·D y where the street network resembles a perfect grid layout. The bus network design problem is aimed at minimizing the total cost of the system, composed by the infrastructure and the operating cost incurred by the transit agency, the emission cost caused by vehicles and the time cost by transit users. The modeling approach to design the bus network will be based on the previous work of Daganzo and Estrada et al. [25,40] considering the following features and assumptions: • Street network and bus corridor structure: It was assumed that buses follow a perfect grid configuration. Another assumption is that the city is covered by orthogonal bus corridors, traveling in the west-east and north-south axis and operated in two directions of service (see Figure 1). The headways of the corridors running along the horizontal and vertical directions (x and y axes) are denoted by H x and H y respectively. Moreover, the route spacing between the vertical or horizontal lines are respectively s x = p x ·s and s y = p y ·s, where s is the stop spacing and p x , p y two integer variables. Finally, there is a transfer stop at each intersection between vertical and horizontal routes. Hence, users can move in the network from any point to any other point with one transfer. • User behavior and service characteristics: The network is designed to minimize the total cost at the peak hour, serving a total of Λ passengers with origins and destinations uniformly distributed in the service area D x and D y (in Estrada [25] it is stated that the agreement between the cost and performance estimations of continuous approximation models with OD simulation techniques is within 85%). However, the bus service will also be provided during Ω hours along the day, with an average hourly demand λ(pax/h). The door-to-door travel time of users integrates: the access and egress time spent at a v w walking speed to/from the closest stop, waiting time at stops, riding time and the transfer time to overcome the distance δ between the loading platforms of orthogonal routes.

•
Kinematic characteristics of buses: The maximal cruising speed on both horizontal and vertical lines is assumed to be identical and equal to v. Boarding and alighting time per passenger at stops is τ'. On the other hand, vehicles have a constant acceleration rate a and a maximal instantaneous speed of v max . Therefore, the additional driving time to perform one stop with regard to one vehicle cruising at maximal speed (without stopping) is τ = v max /a. • Vehicle powertrains and energy provision: It was assumed that ICE-powered buses (diesel) or battery electric buses (BEB) would operate the service. In the case of diesel vehicles, the auxiliary facility to provide energy will be a set of fuel dispensers (fuel station) to be located in the bus garage.
and vertical lines is assumed to be identical and equal to v. Boarding and alighting time per passenger at stops is τ'. On the other hand, vehicles have a constant acceler ation rate a and a maximal instantaneous speed of vmax. Therefore, the additional driv ing time to perform one stop with regard to one vehicle cruising at maximal speed (without stopping) is τ = vmax/a. • Vehicle powertrains and energy provision: It was assumed that ICE-powered buse (diesel) or battery electric buses (BEB) would operate the service. In the case of diese vehicles, the auxiliary facility to provide energy will be a set of fuel dispensers (fue station) to be located in the bus garage.
On the other hand, if the service is provided by BEB vehicles, one of the crucial deci sion variables is the energy capacity of the battery pack (E, in kWh). In fact, it was assumed that vehicles must have, at any point of the service, enough energy to run a distance to the bus garage or maintenance location. Despite the fact that the city may present several bu garages scattered in the city periphery, it was assumed that the maximal distance to be run from terminal stop to the bus garage is lG. This can be conceived as a distance margin that every bus may overcome in the necessary case. Assuming an average vehicle con sumption factor of fc (kWh/veh-km), the effective energy that can be used is (E − fc·lG) Therefore, it has to be check that the effective energy capacity (E − fc·lG) provides a suffi cient autonomy range to operate the service all day long. If so, this vehicle will be charged at the bus garage during the night period. This operation is denoted as overnight charging The total number of chargers to be deployed at the bus garage for overnight charging is given by = / / , where M is the fleet size (see Equation (12)) and [x] − a math ematical operator that gives the lower integer of x. The term into brackets is, indeed, the number of buses that a single charging station is able to charge completely their batterie at full capacity (E) at an electric power SD during the available night time. Otherwise, i On the other hand, if the service is provided by BEB vehicles, one of the crucial decision variables is the energy capacity of the battery pack (E, in kWh). In fact, it was assumed that vehicles must have, at any point of the service, enough energy to run a distance to the bus garage or maintenance location. Despite the fact that the city may present several bus garages scattered in the city periphery, it was assumed that the maximal distance to be run from terminal stop to the bus garage is l G . This can be conceived as a distance margin that every bus may overcome in the necessary case. Assuming an average vehicle consumption factor of f c (kWh/veh-km), the effective energy that can be used is (E − f c ·l G ). Therefore, it has to be check that the effective energy capacity (E − f c ·l G ) provides a sufficient autonomy range to operate the service all day long. If so, this vehicle will be charged at the bus garage during the night period. This operation is denoted as overnight charging.
The total number of chargers to be deployed at the bus garage for overnight charging is given by , where M is the fleet size (see Equation (12)) and [x] − a mathematical operator that gives the lower integer of x. The term into brackets is, indeed, the number of buses that a single charging station is able to charge completely their batteries at full capacity (E) at an electric power S D during the available night time. Otherwise, if BEB vehicles are equipped with small battery packs, it may be charged on-route. This operation is known as opportunity charging. From now on, it will be assumed that opportunity fast chargers are only located at terminal stops. The regular case considers that there are N y charging stations evenly distributed along distance D x in the North or the South periphery quadrant of the city.
Due to urban space scarcity, it was considered that all these stations are located a vertical distance l oy far from the city boundary. Similarly, there are N x stations located in the east or west periphery quadrant evenly distributed along D y distance. All these N x stations are situated l ox distance away from the vertical boundary of the city. Therefore, vertical and horizontal corridors will have one terminal stop in each round trip where the battery charging operation is available. Nevertheless, it can be also considered that N x and N y additional charging facilities can also be deployed symmetrically in the opposite periphery quadrants, allowing buses to be charged at each terminal stop of each bus route direction. Doing this, vehicles may be equipped with smaller battery packs since bus charging operation may be more frequent than before. The E-W and N-S periphery quadrants with charging stations are controlled by decision variables φ x and φ y . When these variables are equal to φ x = 2; φ y = 2, it means that charging stations are located in both periphery quadrants E-W and N-S.
Depending on the location of the chargers in W-E quadrants, the vehicles operating the last part of a horizontal corridor will run along a vertical detour distance d v to the charger. This distance is, in fact, a function of the number of chargers located at W-E quadrants, . Similarly, vehicles operating vertical corridors will incur an extra horizontal distance d h = d h N y based on the number of available chargers. Therefore, when φ x = 2; φ y = 2, the distance range between two charging operations in vertical and horizontal corridors are respectively denoted by d c,S = d c,N = D y + 2 l oy + 2d h , d c,E = d c,W = D x + 2 l ox + 2d v . They include respectively the longitudinal distance run within the city boundaries, the longitudinal distance travelled out of the city boundaries, and the detour distance in the transversal direction.
On the contrary, the situation described by φ x = 1; φ y = 1 is equivalent to the regular case, where there is a charging station in just one terminal stop. In that case, the distance run is d c,S = 2D y + 2l oy + 2d h and d c,E = 2D x + 2l ox + 2d v . Hence, the total number of charging stations to provide service in vertical and horizontal corridors is respectively θ y N y and θ x N x . Note that when the charging time of a single bus at a terminal stop is higher than the target time headway H x or H y , vehicles cannot be dispatched at the required frequency unless multiple charging areas are deployed. Therefore, several buses must be charged simultaneously at the terminal stop. Hence, it was considered that each charging station y (y = 1, . . . , N y ) located in the N-S quadrants will be deployed with m y charging areas. Each charging area has a parking slot for one bus and is equipped with an independent charging connector. Therefore, a total of m y buses can be charged simultaneously at each charging station of the set N y . Similarly, each charging station x (x = 1, . . . , N x ) located in E-W periphery quadrants will be equipped with m x charging areas.
The number of charging areas in each charger station (m x or m y ) is defined to be the minimal value that satisfies the required charging time and the dispatchment of vehicles at the target headways H x or H y . To do so, it was considered that the charging time of one vehicle at charging station y is T c,y (y = 1, . . . , N y ). In addition to that, it was considered that the expected number of routes served at each charging station is D x N y p x s . If the average headway between two consecutive bus departures from each charger y is h y , therefore the minimal number of charging areas that each terminal stop y in north and/or south peripheral quadrants must present is given by m y = T c,y h y In the W-E periphery quadrants, the number of charging areas is similarly calculated Therefore, the total number of opportunity charging areas is Another available charging alternative is the daily charging operation at the bus garage, where vehicles are charged at this facility during the service. As a difference from the opportunity and overnight charging, vehicles present very different levels of SOC during the service. The modeling approach in the following lines is developed for a horizontal route. Let µ x = 2D x /(H x ·v cx ) be the number of vehicles needed to operate a single horizontal corridor to fulfill the target headway H x . It was assumed that after ∆t 0 units of time from the beginning of the service, a subset of ∆µ x vehicles (∆µ x ≤ µ x ) will go to the bus garage to be charged. These vehicles will be replaced in Sustainability 2021, 13, 4745 6 of 28 the service by the same amount of extra vehicles that were previously parked at this facility, with batteries completely charged. These ∆µ x vehicles will spend l G /v units of time traveling to the bus garage, ∆t 0,x ·v c,x · f c S d units of time being charged (they have only provided service during ∆t 0,x ), and l G /v units of time for traveling again to the bus route to continue the service. Therefore, the time needed during this first cycle of charging is At time ∆t 0,x + ∆t 1,x , we can start the second cycle of charging, replacing those ∆µ x vehicles pending to be recharged, by the ones coming from the bus garage. Note that the charging operation of these second group of vehicles needs ∆t 2,x = 2l G v + (∆t 0,x + ∆t 1,x )·v c,x · f c /S d units of time to be reintroduced again into the service. Generally, the ∆µ x vehicles replaced at the n-th cycle of charging will consume ∆t n, units of time for being available again. Therefore, this charging scheme may charge vehicles during the whole service at the expense of deploying ∆µ x additional vehicles in the operation, and performing the unproductive trips between routes and bus garage.
The variables that control the charging operations are the number of charging cycles performed (π x ), the number of extra vehicles (∆µ x ) and the elapsed time when the first charging operation is performed (∆t 0,x ). In order to ensure a feasible charging scheme, it has to fulfill on one hand that all charging cycles are performed within the service period (∑ π x i=0 ∆t i,x ≤ Ω) and, secondly, the added vehicles and the number of charging cycles are enough to charge all vehicles operating the route (π x ·∆µ x ≥ µ x ). Similar equations must be fulfilled in the vertical corridors (y coordinate). Consequently, the bus network design problem is aimed at minimizing the total costs of the system (Z), defined by Equation (1).
The first term in brackets in Equation (1) captures the infrastructure cost that the transit agency will incur, where L is the total length of bus corridors (two ways of service) and n ch the number of electric charging areas (BEB vehicles) or the number of refueling stations (ICE vehicles). The corresponding parameters c L and c N determine the depreciation cost per unit of bus lane distance, and per unit of charging area and hour respectively. In the case of diesel vehicles, c N is calculated as the depreciation cost of one fuel station per hour, divided by the number of refueling operations that a single station can perform along the night period.
The second term in brackets estimates the agency's operating cost. The total distance run by the fleet in one hour (V) and the fleet needed (M) are multiplied by the corresponding Sustainability 2021, 13, 4745 7 of 28 unit distance cost (c V ) and unit temporal cost (c M ). The former considers all operating costs depending on the distance run (mainly energy), while c M embraces vehicle depreciation, driver salary, insurances and other fixed expenses. Typically, the depreciation cost of electric buses does not encompass the depreciation cost of batteries (c B ), since the battery lifetime is significantly lower than the vehicles. Indeed, the battery capacity (E) has to be determined, according to the route length and energy consumption factors.
Moreover, the user cost contribution is calculated in the third bracketed term by the product of the average hourly demand rate (λ in pax/h), the value of time (β T in USD/pax-h) and the total door-to-door travel time (T, in hours). This travel time component encompasses access (A), waiting (W), transfer (Tr) and in-vehicle travel time (IVTT). Finally, it was also included in the objective function the monetary cost of emissions per hour of service (Z E ). Therefore, the total cost function Z (USD/h), in terms of monetary units per hour of service.
Equations (2a) and (2b) constraint that the occupancy of vehicles running along x and y directions (O x and O y respectively) should be lower than the maximal passenger capacity of buses (C in pax/bus). The set of Equation (3) is aimed at verifying that vehicles present a sufficient range under three charging schemes. If we opted for an overnight charging scheme at the bus garage, Equation (3a) states that the total distance run along the whole day by BEB vehicles, plus a safety margin distance l G , should be lower than the maximal allowable distance (E/ f C ). In the former equation, v c,x and v c,y represent the number of kilometers run per hour in the x and y direction by a single bus (net commercial speed calculated through Equations (11a) and (11b). Equation (3b) ensures a sufficient range for the group of vehicles charged at the first (i = 1) and last cycle (i = π x or π y ) of charging at the bus garage during the service. Similarly, Equation (3c) should be only satisfied by BEB technology with the opportunity charging scheme. Equation (4) constraints the minimal threshold of stop spacing (s min ) and headways (H min ). It was assumed that below these two thresholds, buses would stop too frequently along the route or they would not be able to maintain an acceptable headway adherence. Finally, constraints (5)-(7) specify the feasible domains and the nature of the decision variables. Indeed, it was restricted that routes will be spaced, at most, twice the stop distance.
In the following subsections, it will calculate each component of the objective function dependent on the decision variables of the problem: s, H x , H y , p x , p y , N x , N y , φ x and φ y .

Bus agency metrics
The different cost components incurred by transit agencies are estimated, considering similar approaches presented in Daganzo and Estrada [17,25]. Equation (8) captures the total bidirectional length of the network (L), while Equation (9) estimates the distance run in one hour by the whole fleet (V). It was assumed that in the case of BEB vehicles with opportunity charging scheme, each vehicle operating a horizontal corridor runs an extra distance l ox + d v from the terminal stop to the closest charging station. This extra distance is equal to l oy + d h for vehicles running along vertical corridors. In diesel and BEB powertrains, l ox = d v = l oy = d h = 0. In the case of day charging at the bus garage, the terms ∑ π x i=0 ∆t i,x and ∑ π y i=0 ∆t i,y represent the time needed to charge all vehicles starting the service (one charging cycle) in the horizontal and vertical corridors. In other charging schemes and diesel powertrains, these terms must be ∑ Moreover, it was calculated the expected number of transfers p(1) per user through Equation (10). This variable is needed to compute the net commercial speed v c in x and y coordinates, addressed in Equation (11). The net commercial speed is the average speed that buses will experience to perform a complete round trip of length 2D x (horizontal) or 2D y (vertical). In the case of BEB with opportunity charging scheme, the total length is computed by 2(D x + l ox + d v ) (horizontal) or 2 D y + l oy + d h (vertical). The net commercial speed takes into account the lay-over time at terminal stops (θ E , θ W , θ S , θ N ) and the charging time T c,z spent at the corresponding terminal station z (z = W, E, S or N) only by BEB vehicles. Finally, the fleet size needed to run the service is addressed by means of Equation (12).
In the case of BEB with the opportunity charging scheme, the charging time T c,z spent at terminal station z is calculated in Equation (13), where z represents the charging stations located on the north, south, east and west part of the city. The first term is the minimal charging time whose estimation is based on the energy consumed in the distance run from the last charging operation d c,z , the energy consumption factor f c and the charging speed S D . The second term captures the positioning and clearance time spent at the charging area per bus in the terminal stop, t pos . This time is neglected in alternative solutions that do not visit the on-street opportunity charging facility (T c,z = 0).

User metrics
The expected door-to-door travel of a single user T is calculated as the summation of access (A), waiting (W), transfer (Tr) and in-vehicle travel time (IVTT) through Equations (14)-(17) respectively. Equation (15) the estimation of the waiting time specifies the different headway provision along directions x and y. Moreover, the IVTT term must take into account the gross commercial speed of buses v c , i.e., the average commercial speed experienced by users, neglecting the idle time spent at the terminal stops. Finally, Equation (18) defines the expected passenger occupancy in the vehicles operating routes in the x and y directions.

Emissions
The estimation of the local emissions and greenhouse gases (GHGs) generated within all the processes needed to run the bus service is addressed in this section. These emissions have been calculated by means of a lifecycle analysis of the vehicle manufacturing, infrastructure construction and the well to wheel analysis of the energy required to operate the vehicles, given a bus powertrain and bus size. These contributions are estimated by means of the corresponding emission factor of a specific pollutant generated by one process (vehicle operation, manufacturing and infrastructure construction). These factors will be later multiplied by the corresponding size fleet (M), infrastructures (L and n ch ) and distance run (V) in one hour of service.
The local pollutants considered are PM 10 , SO x , CO, NO x , NH 3 and VOC, while CO 2 , CH 4 and N 2 O are transformed into CO 2 eq to account for GHG, taking the corresponding warming potentials. Both local pollutants and CO 2 eq constitute the set P of different emission typologies considered in this paper. Finally, the mass of each pollutant generated per hour of service has been monetized to be integrated in the objective function. To do so, a proxy parameter ε p has been defined to monetize the effect of producing one unit of pollutant p ∈ P. The values of these proxy parameters (in USD/g of pollutant p) were determined according to van Essen et al. [40]. The processes considered that have an associated economic impact in terms of emissions and pollutants in the bus system are the following: Monetization of the emissions caused by the circulation of vehicles: Z e,F : Each ICE vehicle produces E F,p (g/veh-km) grams of local pollutant p (p ∈ P) per kilometer run. In fact, these emission factors capture the pollutants generated in the tank-to-wheel phase. These emission factors, in reality, depend on the type of vehicle, speed regime, driving style and fuel consumption factors [33]. The emission monetization associated to the circulation of vehicles in one hour can be estimated by Equation (19), taking into account the total distance run by the whole fleet (V) in the peak hour. It is considered that BEBs do not present any pollutant contribution.
Monetization of the emissions caused by the energy generation in production plans: Z e,E : The generation of electricity from primary energy sources for bus battery charging also lead to emissions when the service is operated by BEBs. These emissions, in both BEB and ICE vehicles, are estimated by means of the well-to-tank analysis. The emission factor E E,p of pollutant p (g of pollutant p/kWh) is estimated as the mass of this pollutant caused in the generation and transportation of one unit of energy to move vehicles. These emission factors depend on the energy mix of the country. Equation (20) captures the estimation of the Z e,E variable, where f c (kWh/veh-km) represents the energy consumption factor of the bus fleet.
Monetization of the emissions caused by vehicle manufacturing and maintenance: Z e,M : The manufacturing of public transport vehicles consumes and transforms a wide number of raw and semiprocessed materials. The emission factor of pollutant p, E M,p (g of pollutant p/veh-h) is defined as the ratio among the total amount of emissions of this pollutant p in the vehicle production and maintenance and the number of working hours along the vehicle lifetime. The emission monetization of this effect is obtained by Equation (21), multiplying emission factors by the proxy parameters (ε p ) and the fleet size (M).
Monetization of the emissions caused by the row infrastructure construction: Z e,I : The manufacturing of the infrastructure needed to build that transit lane also contribute to the GHG and local emissions. In that case, the emission factor for pollutant p is considered per unit of infrastructure length and unit of time used, E I,p (g/km-h). It was also accounted for the contribution of public work machinery used in the whole construction phase of the bus network. The total amount of emissions produced in one kilometer is split along the target lifetime of the infrastructure. Then, the monetization can be obtained by Equation (22).
Monetization of the emissions caused by the stop/station construction: Z e,S : This term encompasses the economic effects of the emissions generated in the construction of stops, stations and charging stations. Traditionally, the material and energy consumption in the construction of the bus stop was significantly lower than the bus lane infrastructure. It was assumed that E S,p (g/st-h) and E C,p (g/st-h) are the emission factors of pollutant p per hour of stop and electric charging area respectively. The emissions caused by the construction of refueling stations for diesel vehicles at bus garage are considered in the E E,p term, and therefore they were not addressed here. Consequently, the term Z e,S can be calculated by Equation (23).
Optimization procedure The analytical model presented before allows identifying the optimal design of the public transport network from 5 decision variables for ICE powertrains, or for BEB powertrain when overnight charging is chosen: headways in horizontal and vertical routes, stop spacing and route spacing between vertical and horizontal lines defined by the auxiliary variables p x and p y . Nevertheless, in the case of the deployment of electric bus technologies with opportunity charging, four additional decision variables are considered to account for the number of electric charging stations and charging areas in the N-S periphery (θ y , N y ) or in the E-W periphery (θ x , N x ). The optimization process followed in this paper was based on a grid enumeration of the objective function, in relation to the decision variables. For each of the variables with a continuous domain, it was identified an enumeration interval pace (∆s and ∆H) to discretize the enumeration domain. The discretized domain of these variables is therefore enumerated by s = s min + k s ∆ s and H ., k max. Input parameters s min , H min define the minimal reasonable value of stop spacing and time headway respectively, while k max is defined based on the maximal value of the former decision variables. In the case of the discrete variables, the enumeration domain of p x , p y and φ x , φ y was {1, 2}, while the corresponding to the number of chargers was N x = {1, 2, . . . , D y s·p x } and N y = {1, 2, . . . , D x s·p y }. Therefore, the total cost of the collective transport system Z (s', H', p x ', p y ', N x ', N y ', φ x ', φ y ') is evaluated by Equation (1) for each potential combination of decision variables defined in the grid search. Hence, the feasible solution of the optimization problem stated in Equation (1) is the combination of decision variable values that produce the minimal cost of the objective function (Z), satisfying the constraints of Equations (2)-(7).

Materials
The applicability of the model presented has been tested under various scenarios of bus powertrains in Guadalajara, México ( Figure 2). The city of Guadalajara, located in the Jalisco state, is the second most populated city in Mexico, with near to 5 million inhabitants within its metropolitan area spread in an area of 151.6 km 2 [41].

Materials
The applicability of the model presented has been tested under various scenarios of bus powertrains in Guadalajara, México ( Figure 2). The city of Guadalajara, located in the Jalisco state, is the second most populated city in Mexico, with near to 5 million inhabitants within its metropolitan area spread in an area of 151.6 km 2 [41]. Public transport services the Guadalajara Metropolitan Area by means of a door-todoor bus network, operated by 12 m-long diesel buses. The bus service area has been modeled by a rectangle of sides Dx× Dy = 18·15 km 2 . The peak hour demand of bus passengers in a working day was estimated to be Λ = 333,613 pax/h, while the average hourly demand considering the variability along the whole day was λ = 236,605 pax/h. Both figures have been calculated from the data provided in Parada [42]. Moreover, it was assumed that the value of time of passengers using the bus service in Guadalajara is μ = 3.14 USD/pax-h, considering that the average monthly salary in 2018 was 514 USD [43]. It was assumed that these users show an average walking speed of vw = 4.5 km/h and the average transfer distance between loading areas at stops is δ = 0.3 km [25]. The service is provided along Ω = 16 h per day. Public transport services the Guadalajara Metropolitan Area by means of a door-todoor bus network, operated by 12 m-long diesel buses. The bus service area has been modeled by a rectangle of sides D x × D y = 18·15 km 2 . The peak hour demand of bus passengers in a working day was estimated to be Λ = 333,613 pax/h, while the average hourly demand considering the variability along the whole day was λ = 236,605 pax/h. Both figures have been calculated from the data provided in Parada [42]. Moreover, it was assumed that the value of time of passengers using the bus service in Guadalajara is µ = 3.14 USD/pax-h, considering that the average monthly salary in 2018 was 514 USD [43]. It was assumed that these users show an average walking speed of v w = 4.5 km/h and the average transfer distance between loading areas at stops is δ = 0.3 km [25]. The service is provided along Ω = 16 h per day.
The infrastructure cost to build the bus lane along the transit corridors was supposed to be c L = 84.36 USD/km-h [23]. Eight scenarios regarding different vehicle capacities, powertrains and charging operation are considered in the analysis: The scenario C-12 represents the current standard size (12 m long) and diesel vehicle powertrain use in Guadalajara, without significant treatments of emissions (conventional diesel). A similar scenario entitled C-18 was created for articulated diesel buses (18 m long) with the same conventional environmental label but higher passenger capacity. Scenarios EuroVI-12 and EuroVI-18 will consider a situation with the most recent diesel-hybrid technology. They differ in terms of the vehicle typology, and consequently, in the maximal passenger occupancy that they can accommodate. Finally, the last four scenarios suppose that charging facilities and BEB vehicles would be deployed. For standard buses of 12 m of length, it was considered that overnight charging at the bus garage, day charging at the bus garage during service and opportunity charging at terminal stops are feasible.
For articulated buses of 18 m, only the opportunity charging operation was considered in the analysis.
The vehicle capacities of all standard and articulated buses are respectively C = 70 pax and C = 120 pax. It was assumed that buses of all the scenarios present the same cruising speed v = 30 km/h, unit boarding and alighting time per passenger τ = 3 s [18], and the additional running time to perform one stop due to acceleration and deceleration phase is τ = v a = 35 s. On the other hand, the maximal acceleration/deceleration rate is a = 0.85 m/s 2 [44]. The parameters that depend on the vehicle typology are summarized in Table 1. The driver cost component (Table 1, row c) was reported to be 9.45 USD/veh-h in Guadalajara, considering the current salaries [45]. The unit distance cost of the ICE vehicles in the first four scenarios was calculated considering the cost of the diesel consumption [32], and an additional cost term that accounts for vehicle maintenance staff and spare parts.
The latter cost component (row b), refuel workload (row e) and engineering staff (row f in Table 1) were calculated considering the real cost reported in Barcelona by Transportes Metropolitanos de Barcelona (TMB) [46] and multiplied by the ratio η GB = 0.2505. This ratio, calculated as the driver cost in Guadalajara divided by the corresponding value in Barcelona (37.72 Euro/veh-h), accounts for the different labor regulations, salaries and monetary parity between these cities. The vehicle purchasing costs, charger and refueling station investments are supposed to be equivalent to the figures reported in TMB [46]. The corresponding cost per hour of vehicle or hour of facility were calculated through a linear amortization along 79,200 h (15 years·330 day/year·16 h/day). The cost of each refueling station is divided by the number of buses whose fuel tank can be refilled in the same day (350 buses), and expressed in terms of USD/veh-h. It is worth to mention that the investment in the charging stations for standard buses at the bus garage is significantly lower than the corresponding to buses charged on the street (opportunity charging).
The battery cost presents a significant variability depending on the different available chemistries, lifespans and acquisition conditions. Organization for Economic Co-operation and Development and International Energy Agency (OCDE-IEA) [45] suggests that the battery cost should range between 120 and 400 USD/kWh in all BEB analysis. Indeed, they suggest a normal value of 260 USD/kWh. Nevertheless, in our study, the resulting battery cost was calculated as a straight-line depreciation of the maximal purchasing cost threshold 400 USD/kWh, along a lifetime of 4 years.
The unit emission costs (ε p ) and emission factors (E F,p , E E,p , E M,p , E I,p , E C,p ) for each pollutant p ∈ P, effect and vehicle scenario in the case applicability in Guadalajara are presented in detail in the Table A1 in the Appendix B. The unit emission costs parameters were considered the average values of contingency cost proposed by van Essen et al. [40] for the 28 European Union countries. In the case of CO 2 eq, the monetization factor was considered ε CO2 = 111 USD/Ton CO 2 based on the impact and wears on the territory. Reasonably, it is one order of magnitude greater than the figure proposed in SENDECO, since the latter is based on the market equilibrium price of CO 2 emission rights among countries.
The terms E F,p accounting for the tank-to-wheel emissions are calculated considering EEA using the Tier 3 estimation for the urban buses standard 15-18 t and urban buses articulated >18 t (diesel conventional and diesel EURO VI diesel powertrains), when the road slope is 0 and it is fully loaded. The non-exhaust emissions like fuel evaporation or emissions generated in the brake operations are not considered. Therefore, the corresponding emissions factors of this category for electric buses are considered to be 0. The emissions E E,p produced in the energy generation and transportation for diesel and battery electric bus powertrains within the well-to-tank phase were calculated using the GREET© fuel cycle model v2019 [30]. It considered the electricity mix for Mexico in 2014 published in IEA [47] based on the following primary energy sources: coal 11.3%, oil 11%, gas 57.1%, nuclear 3.3%, hydro 13%, bioenergy 0.3% and others clean sources (4%).
The emissions E M,p caused by the vehicle manufacturing process are approximated to the results presented in Nordelöf et al. [35] for standard diesel and electric buses, based on the data of Volvo's 7900 model series. The corresponding values for articulated buses were supposed to be 1.35 times higher than the values for standard buses. The emissions were split by a total vehicle lifetime of 15 years (diesel and electric), assuming 330 days of work along the year and 16 working hours per day. On the other hand, the emissions E I,p caused by the infrastructure construction were calculated based on Wayman et al. [48] to build 1 square meter of pavement surface. Assuming a lifespan of 60 years (including maintenance works, 360 days·16 h of service/day) and a bus lane width of 3.5 m, the CO 2 and SO 2 emissions were calculated per hour and kilometer of bus lane. Emissions factors of other pollutants were not available. Finally, the emissions factors associated to each stop construction were negligible in our analysis due to the minor requirements of bus shelters in Guadalajara.
Nevertheless, the emission factors for charging facilities were calculated based on Nansai et al. [49] regarding the installation and manufacturing phase of one charger (14.55 Ton CO 2 , 16 kg NO x , 15 kg SO x and 13 kg CO). Additionally, it was also considered the emissions of the charger transportation along 40 km from the assembling plant to the final location with a heavy diesel duty truck 14-20 t (EURO VI). For each effect, the summation for all pollutants p of the corresponding emission factor, times the unit emission cost is summarized in Table A1. In the case of pollutant x = CO 2 eq, the emission factors considering the vehicle manufacturing, operation and the well to tank phase in Scenarios C-12, Euro VI 12 and BEB 12 Ov. are 3153, 2763 and 2306 g CO 2 /veh-km (see Appendix B). In the case of the electric vehicle (Scenario 5), the well to tank emission factor is 514 g/kWh, similar to the 582 g/kWh for México [50].

Results
This section shows the final results derived from the model. Figure 3 shows the total cost of the bus system obtained by the optimization model in different powertrain scenarios, and in the current Guadalajara bus network configuration, operated by conventional standard (12 m) diesel buses.

Sensitivity Analysis
In this section, it was analyzed how the total cost of the system and the service performance metrics evolve in the different powertrain scenarios according to the steady changes in the demand, city size and in the unit electric vehicle and unit charging infrastructure costs.

Hourly Demand
In Figure 5a, the optimal stop spacing (s), time headway variable (only in horizontal routes, Hx) and the total cost of the system (Z) are depicted for different hourly demand values, ranging from 25,000 to 500,000 pax/h for each of the powertrain scenarios. Figure 5a shows the variation of the optimal network configuration in the total cost (z), the spacing among stops (s) and the headway (H). In general it is noted that all the powertrain scenarios have a similar performance in terms of a linear increase of the total costs as demand increases. A detail discussion of the results is presented in Section 4.
In order to zoom up in which subdomains each scenario is more efficient, it is depicted in Figure 5b the difference between the total system cost Z in each scenario and the function Zm. The Zm function was calculated through a linear regression model, considering the total cost of the scenario (C-12) in the corresponding demand data domain of anal-

Sensitivity Analysis
In this section, it was analyzed how the total cost of the system and the service performance metrics evolve in the different powertrain scenarios according to the steady changes in the demand, city size and in the unit electric vehicle and unit charging infrastructure costs.

Hourly Demand
In Figure 5a, the optimal stop spacing (s), time headway variable (only in horizontal

Sensitivity Analysis
In this section, it was analyzed how the total cost of the system and the service performance metrics evolve in the different powertrain scenarios according to the steady changes in the demand, city size and in the unit electric vehicle and unit charging infrastructure costs.

Hourly Demand
In Figure 5a, the optimal stop spacing (s), time headway variable (only in horizontal routes, H x ) and the total cost of the system (Z) are depicted for different hourly demand values, ranging from 25,000 to 500,000 pax/h for each of the powertrain scenarios.   Figure 5a shows the variation of the optimal network configuration in the total cost (z), the spacing among stops (s) and the headway (H). In general it is noted that all the powertrain scenarios have a similar performance in terms of a linear increase of the total costs as demand increases. A detail discussion of the results is presented in Section 4.
In order to zoom up in which subdomains each scenario is more efficient, it is depicted in Figure 5b the difference between the total system cost Z in each scenario and the function Z m . The Z m function was calculated through a linear regression model, considering the total cost of the scenario (C-12) in the corresponding demand data domain of analysis. The resulting analytical expression of this function was Z m = 2.3728Λ + 69,215.
Finally, Figure 5c shows the results of the variations in the size of the fleet and charging infrastructure required varying the demand.

City Size
The same analysis of the bus electric powertrain resources and cost has been repeated when the city size is varied in the domain D x ∈ [3; 30] km. The basic assumption was that the city shape and demand density is maintained with regard to the Guadalajara current network, therefore D y (km) = 15/18D x and Λ pax h = 336000 15·18 15 18 D 2 x . This sensitive analysis indicates how the size fleet (M), the energy capacity of the battery (E) and the charging numbers changes in functions to the variation of the city size. Figure 6 indicates these results. Finally, Figure 5c shows the results of the variations in the size of the fleet and charging infrastructure required varying the demand.

City Size
The same analysis of the bus electric powertrain resources and cost has been repeated when the city size is varied in the domain ∈ [3; 30] km. The basic assumption was that the city shape and demand density is maintained with regard to the Guadalajara current network, therefore ( ) = 15/18 and Λ( ) = · . This sensitive analysis indicates how the size fleet (M), the energy capacity of the battery (E) and the charging numbers changes in functions to the variation of the city size. Figure 6 indicates these results.

Electric Vehicle and Charging Infrastructure Costs
The former efficiency comparison among bus scenarios mainly depends on the cost parameters related to the energy provision, vehicle purchase and the charging/refill facility investment and their depreciation models. Although the acquisition cost of diesel vehicles is quite stabilized, the competition between electric vehicles is evolving rapidly. As a result, it may be expected that the current depreciation cost of electric vehicles cM is going to decline in the future years. A similar variability may be supposed for the unit vehicle distance cost cV and the depreciation cost of charging stations cN for electric powertrains.
Therefore, we will analyze the potential efficiency of electric powertrains in the Guadalajara city for different unit distances and temporal-based vehicle cost domains. Indeed, we will calculate the unit distance cost (USD/veh-km) for electric vehicle of typology x by ( ) , = ( ) , (1 + ), x = 12 or 18 m, where ( ) , is the corresponding unit distance cost of diesel Euro VI powertrain and a reduction factor (−1 ≤ ≤ ∞). Similarly, it was assumed that the temporal-based cost of electric vehicles (USD/veh-h) can be calculated by ( ) , = ( ) , (1 + ), x = 12 or 18 m, where ( ) , is the corresponding unit temporal cost of diesel Euro VI powertrain and an increment factor (−1 ≤ ≤ ∞). It should be noted that the new temporal based cost of electric vehicles does not take into account the hourly depreciation cost of the battery pack.
In Figure 7a, it was identified the domain of the factor pair ( , ) where the utilization of BEB-12 Ov, Opp and BEB-18 Opp reduces the total cost of the system with regard to the deployment of EVI-12 and EVI-18 respectively. In this case, it was plot the boundary

Electric Vehicle and Charging Infrastructure Costs
The former efficiency comparison among bus scenarios mainly depends on the cost parameters related to the energy provision, vehicle purchase and the charging/refill facility investment and their depreciation models. Although the acquisition cost of diesel vehicles is quite stabilized, the competition between electric vehicles is evolving rapidly. As a result, it may be expected that the current depreciation cost of electric vehicles c M is going to decline in the future years. A similar variability may be supposed for the unit vehicle distance cost c V and the depreciation cost of charging stations c N for electric powertrains.
Therefore, we will analyze the potential efficiency of electric powertrains in the Guadalajara city for different unit distances and temporal-based vehicle cost domains. Indeed, we will calculate the unit distance cost (USD/veh-km) for electric vehicle of typology x by (c V ) BEB,x = (c V ) EV I,x (1 + β V ), x = 12 or 18 m, where (c V ) EV I,x is the corresponding unit distance cost of diesel Euro VI powertrain and β V a reduction factor (−1 ≤ β V ≤ ∞). Similarly, it was assumed that the temporal-based cost of electric vehicles (USD/veh-h) can be calculated by (c M ) BEB,x = (c M ) EV I,x (1 + β M ), x = 12 or 18 m, where (c M ) EV I,x is the corresponding unit temporal cost of diesel Euro VI powertrain and β M an increment factor (−1 ≤ β M ≤ ∞). It should be noted that the new temporal based cost of electric vehicles does not take into account the hourly depreciation cost of the battery pack.
In Figure 7a, it was identified the domain of the factor pair (β M , β V ) where the utilization of BEB-12 Ov, Opp and BEB-18 Opp reduces the total cost of the system with regard to the deployment of EVI-12 and EVI-18 respectively. In this case, it was plot the boundary that limits the region where BEB was preferable (bottom-left) for different values of the emission monetization. Hence, we generated the potential monetization values in the future by ε p, f = ε p,a (1 + β ε ), where ε p,a is the current monetization value of pollutant p and β ε an increment factor (−1 ≤ β ε ≤ ∞). In Figure 7b, a similar analysis was done where the boundary limit among the BEB or EVI effective areas was provided for different depreciation values of the charging station. Finally, in Figure 7c, it was analyzed how much the battery depreciation (cB) and charger cost (cn) should be reduced to economically justify the overnight charging, instead of the opportunity charging scheme.

Discussions
The results presented in Section 3 can be dived in two main groups of results, which are the performance metrics for all the scenarios proposed, which were presented in Fig

Discussion of the Performance Metrics
The first set of results that summarize the main performance metrics of all the scenarios can be visualized in Figures 3 and 4. Figure 3 contains the total cost of the bus system obtained by the optimization model in different powertrain scenarios and the current system. It can be observed that the main effect is that the modification of the bus network design is the main responsible to lessen the total cost of the system. Only by varying the spatial and temporal configuration to a transfer-based network, the total cost estimated by the model was reduced by 23% with regard to the current network operated by the same vehicles (Scenario C-12). The total system cost barely varied in different bus powertrains under analysis, once the configuration was based on a transfer-based network, instead of a door-to-door network.
The minimal cost of the system was obtained when electric battery standard buses were deployed with an opportunity charging scheme (Scenario BEB-12 Opp). The incremental cost saving of the previous vehicles typology due to the higher powertrain efficiency was 6.02% (net saving between C-12 and BEB-12 Opp). The second-best vehicle is a diesel Euro VI standard bus (Scenario EVI-12), followed by articulated buses charged in an opportunity operation mode (Scenario BEB-18). The day charging at the bus garage during the day (Scenario BEB-12 Day) was not as competitive as other scenarios.
The charging scheme of standard battery electric buses plays a crucial role in the cost of the system. The opportunity charging scheme of standard BEB vehicles implies more vehicles (ΔM = 264 veh) and runs more kilometers per hour (ΔV = 284 veh-km/h) than In Figure 7b, a similar analysis was done where the boundary limit among the BEB or EVI effective areas was provided for different depreciation values of the charging station. Finally, in Figure 7c, it was analyzed how much the battery depreciation (c B ) and charger cost (c n ) should be reduced to economically justify the overnight charging, instead of the opportunity charging scheme.

Discussions
The results presented in Section 3 can be dived in two main groups of results, which are the performance metrics for all the scenarios proposed, which were presented in Figures 3 and 4 and the sensitive analysis, which were presented in Figures 5-7.

Discussion of the Performance Metrics
The first set of results that summarize the main performance metrics of all the scenarios can be visualized in Figures 3 and 4. Figure 3 contains the total cost of the bus system obtained by the optimization model in different powertrain scenarios and the current system. It can be observed that the main effect is that the modification of the bus network design is the main responsible to lessen the total cost of the system. Only by varying the spatial and temporal configuration to a transfer-based network, the total cost estimated by the model was reduced by 23% with regard to the current network operated by the same vehicles (Scenario C-12). The total system cost barely varied in different bus powertrains under analysis, once the configuration was based on a transfer-based network, instead of a door-to-door network.
The minimal cost of the system was obtained when electric battery standard buses were deployed with an opportunity charging scheme (Scenario BEB-12 Opp). The incremental cost saving of the previous vehicles typology due to the higher powertrain efficiency was 6.02% (net saving between C-12 and BEB-12 Opp). The second-best vehicle is a diesel Euro VI standard bus (Scenario EVI-12), followed by articulated buses charged in an opportunity operation mode (Scenario BEB-18). The day charging at the bus garage during the day (Scenario BEB-12 Day) was not as competitive as other scenarios.
The charging scheme of standard battery electric buses plays a crucial role in the cost of the system. The opportunity charging scheme of standard BEB vehicles implies more vehicles (∆M = 264 veh) and runs more kilometers per hour (∆V = 284 veh-km/h) than diesel and overnight-charged electric counterparts. Nevertheless, the cost reduction of BEB-12 Opp with regard to diesel vehicles (EVI-12) was due to the lower unit distance cost parameter (c V ). This fact reduced drastically the agency cost component. Other cost components were balanced between diesel and BEB-12 Opp technology. Indeed, the emission cost in Scenario EVI-12 was equivalent to the sum of emission and battery depreciation cost in BEB-12 Opp.
However, the BEB-12 Opp. significantly outperformed the overnight charging scheme due to the lower number of chargers and battery packs. BEB-12 Opp involves n ch = 224 fast charging areas (N x = 18, m x = 3, φ x = 2, N y = 29, m y = 2, φ y = 2) and vehicles equipped with a battery pack size of E = 50.96 kWh. On the contrary, the overnight charging scheme will require n ch = 394 areas and a battery pack of E = 463 kWh. Only the depreciation cost of batteries rose to 14% of the total cost, excluding user costs. BEB vehicles with overnight charging scheme would be competitive only if the battery and the garage charging facility cost were reduced. However, due to the battery lifespan (5 years maximum) and the current manufacturing cost, the batteries played a differential role in the competitive analysis. Day charging at the bus garage (BEB-12 Day) requires a similar number of chargers and less battery capacity with regard to overnight charging (n ch = 333 fast ch. areas and E = 367 kWh/veh). Nevertheless, it needs 330 extra vehicles than BEB-12 Ov. due to the death headings to the bus garage and charging time during the service. Moreover, the emission cost becomes a small fraction of the total system expenses in all scenarios. In the current network, C-12 and C-1 scenarios, the monetary value of emissions represents 4-6% of the total cost with the current proxy parameters.
On the other hand, analyzing Figure 4 it can be observed that the optimal time headways and stop spacing vary around H = 2-2.5 min and s = 310-350 m respectively. The existing service in operation outperforms the network accessibility due to the multiple number of routes, which can only be operated at low frequencies (huge waiting times). Hence, the new optimized design reduces waiting times at the expense of rationalizing the number of routes. The optimal bus network configuration is obtained when the stop lattice is p x = p y = 2, i.e., one out of two consecutive stops allows transfer operation to orthogonal routes. In terms of the amount of total vehicles (M) this shows a variation in the proposed scenarios from the smallest fleet requirement, which was 2330 vehicles in the C-18 scenario to the biggest size requirement, which was 3089 vehicles in the BEB-12 Day scenario. The size variation among the scenario is bonded mainly to the size of the vehicles (some of 12 m and other of 18 m see in Section 2) and the charging scheme.

Discussion of Sensitivity Analysis of Hourly Demand
Within the hourly demand (Λ) sensitivity analysis it was observed how total cost of the system (Z), the headway (H) and stop spacing (s) was modified in function to the hourly demand. This result was presented in Figure 5a. The total cost of the system presents an increasing behavior with regard to Λ, since the user cost component is a linear function of Λ. While the demand gets more consolidate, the bus route design tends to be operated at lower headways and line spacings. The different behavior in the stop spacing function around Λ = 40,000 pax/h is because the optimal stop lattice changes from p x = p y = 2 (Λ < 350,00 pax/h) to be p x = 1 and p y = 1 when Λ > 350,00 pax/h. In this domain, the lattice p x = p y = 2 will need a small stop spacing to accommodate the total demand in the number of corridors, causing a speed reduction. The total cost presents barely the same numerical value in all powertrain scenarios.
In the case of Figure 5b this showed the cost saving for each scenario compared with scenario C-12 when the demand is varying. It can be highlighted from Figure 5b that the powertrain that minimizes the total cost of the system in the whole demand domain is still the BEB-12 Opp vehicle, followed by EVI-12 and electric buses (BEB 18 Opp, 12 m-long overnight and day charging), EVI-18 and C-12.
Moreover, it was also calculated the extra cost that the system will incur when it maintains the spatial bus route configuration during the whole day of service. It was assumed that the bus route configuration was designed for the peak hour in Guadalajara (Λ' = 333,613 pax/h) where p x = p y = 2 and s = 310 m. In the whole domain of potential peak demand Λ ∈ [25; 500] kpax/h, it was optimized the bus network design problem again, maintaining the values of spatial variables s, p x and p y equal to the former figures and only varying the time headways H x and H y . This strategy is usually undertaken by transit agencies. The resulting value of the objective function was referred as Z' and it was compared to the value Z opt , corresponding to the design problem when all decision variables could be varied. The convex functions Z'-Z opt for all scenarios are depicted in Figure 5b by dotted lines. It can be observed that the extra cost of the system was negligible around Λ' = 330,000 pax/h (target demand value for the spatial design). It highlights the outstanding robustness of the network design method proposed in this paper. Decision makers may adapt the optimal values of decision variables to the real and available configuration of the street network to draw the final master plan, without incurring a significant extra cost. As the hourly passenger flow from the target demand Λ' is increased or decreased, the extra cost increases. When Λ > 400,000 pax/h, the operation of the prefixed network configuration at the lowest headways in standard buses (12 m long) was not feasible.
Finally, it was analyzed in Figure 5c how the number electric resources (vehicles, battery capacity and chargers) evolved with regard to the hourly demand, maintaining the rest of city and transit agency parameters constant. The battery pack in scenario BEB-12 Ov. varied around E = 500 kWh. Although there are some bus prototypes in the market equipped with batteries of this capacity, this figure exceeds the normal capacity of buses deployed in real service and would become a huge challenge for vehicle manufacturers. On the contrary, the battery pack needed in scenario BEB-12 where vehicles are charged at terminal stops presents a constant value of E = 51.0 kWh. The battery pack in the BEB-12 Day scenario, remained stable around 350 kWh.
The number of charging stations increased as the demand consolidates and as a result, more vehicles were needed. Nevertheless, the ratio M/n ch slightly increased as the demand growths. It remained stable in each charging scheme around 6-12 buses per charger for the whole domain of service.

Discussion of Sensitivity Analysis of City Size
In the case of the variation of the city size this approach was made to analysis the variation in the fleet size (M), the charging infrastructure required (N ch ) and the battery capacity (E) for each electric scenario. This result was presented in Figure 6.
The battery pack in scenario BEB-12 Opp presents in Figure 6 a linear increment with regard to D x , ranging from E bat = 34 to 68 kWh in the domain of analysis. This was an expected tendency since the consumption between two consecutive chargers in horizontal and vertical corridors are respectively D x · f c and 15 18 D x · f c . However, in all scenarios, the number of charging stations, the number of vehicles and consequently, the depreciation cost of batteries and charging facilities, do not present the linear growth with regard to D x variable as happened in Figure 5c. This fact is the result of the hourly demand that grows quadratically with D x .

Discussion of Sensitivity Analysis of Electric Vehicle and Charging Infrastructure Costs
Finally, in the last set of results of Section 3 a sensitivity analysis was made varying cost in charging infrastructure and vehicles. These results were presented in Figure 7a-c.
In Figure 7a, when the c V parameter (units distance cost) in BEB and Euro-VI vehicles was the same (β V = 0) and current emission monetization (β ε = 0), BEB-12 Ov. and BEB-12 Opp. vehicles were only cost competitive if its vehicle depreciation cost was reduced by 60% and 15% respectively, with regard to the corresponding EVI-12 cost. If the emission monetization values are doubled (β ε = 1), the domain of (β MB , β V ) where BEB is costcompetitive grows. Nevertheless, the domain increment was marginal. If externalities are monetized 10 times higher, the BEB technology is cost efficient even when the depreciation cost of electric vehicles is higher than diesel counterparts. On the contrary, the situation when emissions were not considered in the objective function (β ε = −1) was analyzed, the domain of successful (β MB , β V ) for the deployment of electric vehicles was reduced.
In the case of the 7b case, the charging station cost for BEB-x (x = 12 m and 18 m) was given by (c n ) BEB,x = (c n ) EV I (1 + β ch ), where (c n ) EV I = 0.104 USD/charger-h was the depreciation cost of the fuel dispenser station per hour and diesel vehicle and β ch the increment factor. As it can be noticed, any increment of β ch reduced the effective area where the BEB outperformed the diesel Euro-VI counterpart. Nevertheless, the (β M , β V ) boundary was practically insensitive to any change in β ch in BEB-12 Opp.
Finally the results showed in Figure 7c and given the current parameters defined in Table A1 for Guadalajara city, the battery cost reduction was calculated by (c B ) new = (c B ) R (1 + β B ) for both overnight and opportunity charging schemes, where (c V ) R = 0.019 USD/kWh-h is the current battery depreciation cost assuming a purchasing cost of USD 400/kWh and β B the corresponding reduction factor (−1 ≤ β B ≤ ∞). The reduction in the charging infrastructure cost was modeled by the reduction factor β ch defined before. It can be noticed from Figure 7c that the recommended domain for BEB-12 Overnight technology was really sensitive to the depreciation cost of batteries. This cost must be reduced by more than 75% (β B > 0.75) to justify the implementation of overnight charged vehicles with battery packs higher than E = 450 kWh. It means that the purchasing cost of batteries should be even lower than the most favorable scenario of 120 USD/kWh [45].

Conclusions
This paper presents a robust methodology to design efficient bus networks operated by battery electric buses or combustion engine buses. This aspect of the paper is the main contribution of the authors to the knowledge of the field. Only defining a set of parameters addressing the shape of the city, demand attributes and capital-operating costs of the vehicle technology, the proposed model directly calculates the optimal bus service design and the costs incurred by agency, users and environment in terms of emissions. Although the model applicability and charging performance was tested in Guadalajara (México), the methodology is easily transferable to any city with a grid-shaped street network.
In the Guadalajara test instance, the optimal design based on a grid network is responsible for the reduction of the current total cost of the existing door-to-door network by 27%, properly balancing the spatial and temporal variables. It is important to highlight that this reduction corresponds to an analytical model based on continuous approximations. The outstanding total cost savings obtained are similar to the corresponding other cities where the bus network design was switched from direct to transfer-based trips [28]. The magnitude of the cost savings, resulted from a new network design, is higher than the ones that can be obtained when the vehicle technology is replaced by more efficient powertrains. This is due to the lower agency cost (less infrastructure length and reduced fleet size) required by transfer-based schemes, compared to the current system. Moreover, the low pollutant proxies caused the external cost contribution in the system cost to be marginal. Indeed, the external cost represented less than 7% of the total cost in fossil fuel powertrains, while in electric powertrains this figure was less than 1%. This potential saving in the externality cost component in favor of BEB was neutralized by the current expensive purchasing price of electric vehicles and batteries. As a result, the total cost of the system did not vary significantly among the considered scenarios in the optimized network, reproducing different bus powertrains. In addition to that, the network design (spatial and temporal coverage) was quite stable in all powertrains.
However, the results obtained by 12 m-long BEB with the opportunity charging scheme outperformed the others, although it usually implies more vehicles due to the charging time spent on-route. The second-best powertrains were EVI-12, followed by BEB-18 and BEB-12 Ov and BEB-12 Day. This ranking is maintained for a wide domain of demand, considering the current purchasing and maintenance cost of resources. However, considering the sensitivity analysis performed, the increment of the emission monetization would allow a greater justification of BEB technology, even if the vehicle temporal cost parameter was currently higher than the alternative EVI technology counterpart. The variation of the facility cost of charging stations practically did not have a high impact on the domains of other cost parameters. Another thing to be aware is that the monetization costs for local pollutants are given as average values for countries in terms of Euros/vehiclekilometer, and do not consider the population density in the areas where transport projects are going to be deployed. Therefore, the temporal cost component of the objective function is affected by the high demand (as in the case of Guadalajara), while the negative health effects on the population only depends on the fleet mileage. The authors recommend the estimation of more realistic and case-sensitive external cost parameters to be balanced with other cost components.
From the vehicle design perspective, the performance of BEB-12 overnight technology is questionable due to the thermal management, weight and space needed on-board by the huge battery pack. The cost of this extremely large battery pack reduces the efficiency of this technology in favor of the opportunity charging alternative. The day charging at bus garage in service is only cost-efficient if the unproductive dead-heading travel time to bus garages are negligible. Indeed, in the BEB-12 Opp. scenario, the battery pack presents a treatable size that minimizes the time spent at chargers.
Nevertheless, the major challenge of the opportunity charging scheme is the construction of on-street charging stations. The best solution for these transit services has been obtained when routes are operated with low time headways. In order to completely charge the battery pack at the terminal station and allow some headway disruptions, each charging station will need two or more bus areas (platforms) to charge simultaneously multiple vehicles at the same time. This proposal is questionable from the city perspective, due to the public space consumption. Due to the space scarcity in the city areas, several authorities are reluctant to allocate areas in the sidewalks to install charging facilities. Moreover, the same authorities also argue that the on-street charging stations reduce the flexibility of the bus system, since routes cannot be extended or modified. In spite of the outstanding cost-efficiency of the opportunity charging scheme, there are still several controversial issues to be solved to clearly promote this charging operation.
Therefore, the adoption of BEB technology presents several challenges due to the higher depreciation cost of batteries and vehicles, and additional number of charging resources. Decision makers and transit agencies must evaluate the new incremental resources and even mileage caused by electrification in order to maintain the target service in terms of headways and the objectives of urban mobility plans.
Author Contributions: All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.
Distance travelled by the fleet in one hour of service, V: In each horizontal line, there is a vehicle that will run a complete round trip in every time headway H x . If all horizontal lines are considered, the total distance travelled by the fleet in the Cartesian component x will be equal to twice (double direction) the total length of the horizontal corridors L x divided by their step interval H x (Equation (A2a)). Nevertheless, in the opportunity charging scheme, we have to compute the deadheading distance of all vehicles operating horizontal corridors that must stop at the electric charging station, L ch,x . The resulting formula is given in Equation (A2b). The horizontal component of deadheading distance in one corridor along x-axis is given by l ox . The vertical component d v depends on the relative spacing between horizontal corridors and the spatial distribution of charging stations in the E-W periphery sectors of the city (see Figure 1). If charging is evenly distributed and centered along the city side D x , this distance can be approximated by In the former expression, d h ∼ = D x 4N y when N y < D x s x . In the case of day charging at the bus garage during the service, each vehicle of the fleet operating horizontal corridors will travel a distance equivalent to 2l G every charging cycle, i.e., every ∑ π x i=0 ∆t i,x units of time. The fleet size required to operate the horizontal routes without charging operations is given by 2L x v c,x ·H x , as it is proved in Equation (A5). Therefore, the extra distance run by the horizontal fleet in one hour due to the dead heading movements to the bus garage is given by 2L Similar expression can be found for vertical corridors.
Finally, the estimation of the total distance run by the whole fleet in one hour in the BEB-12 day scenario is made by Equation (A2c).
Average number of transfers, p(1): Let us assume that there is a user's trip that originated at point O1 of Figure A1a. The user should not make transfers if its destination is located within a longitudinal band of width p y s or a vertical band of width p x s, both centered at point O1. The corresponding area of these bands is D x p y s + D y p x s. Thus, the probability of making zero transfers is equivalent to the probability that the travel destination D1 is contained within these bands. Assuming that origins and destinations of trips are uniformly distributed in the network, the probability of making 0 transfers is equivalent to the ratio between the area of the aforementioned bands and the total area of the city. However, the area resulted from the intersection of the two bands will be subtracted from the numerator to avoid double-counting. Finally, the probability of making a transfer is easily determined knowing that p(0) + p(1) = 1. In this way, the expected number of transfers is equivalent to the probability of making 1 transfer, p(1), defined in Equation (A3).
Net commercial speed, v c : The pace of vehicles, i.e., the inverse of the net commercial speed, is calculated as the summation of the time needed to overcome one kilometer at the maximal cruising speed v, the additional time for braking and accelerating the vehicle before/after all stops located in one kilometer, and the boarding and alighting time of users in one kilometer. Note that in the latter component, the number of boarding passengers along horizontal corridors also embraces those passengers that transferred from vertical to this horizontal route p(1)Λ/2. This figure is divided by the distance run in service per hour, V x (the one given by Equation (A2a)) to calculate the boarding passengers on one vehicle along one kilometer. These three components correspond to the first three terms of Equation (A4) and are also calculated in Estrada [25]. Nevertheless, the fourth term consists of the layover time at the two terminal stops of the horizontal (θ W , θ E ) and vertical routes (θ N , θ S ). Moreover, the fifth term is the positive difference between the time spent at charging vehicles at the terminal stops and the previous lay-over time in the case of opportunity charging. If the charging time is higher than the layover time, vehicles are held for a longer time at terminal stops, becoming less productive and shrinking the net commercial speeds.
Number of vehicles required to operate the service in one hour, M: The size of the fleet can be easily calculated by means of Equation (A5), breaking it down into its horizontal (M x ) and vertical component (M y ). The number of kilometers operating the horizontal corridors can be directly calculated by the quotient between the number of kilometers run along horizontal corridors in one hour (V x ) and the net commercial speed (v c,x ). A similar approach can be done to calculate the fleet size in vertical corridors. In the case of opportunity charging we have to consider the time spent travelling the unproductive distance to access the chargers L ch,x and L ch,y . In the case of the day charging at the bus garage, we consider the extra fleet needed in each horizontal and vertical route (∆µ x and ∆µ y ) multiplied by the number of horizontal and vertical corridors respectively.
Access time, A: This component will depend on the spacing of the stops and the transport lines (p x , p y ). We assume that the probability to board on a vertical and horizontal corridor is equal. Therefore, access distance is calculated (see Figure A1b) as one-half of the summation of the average distance to a horizontal corridor a H (from brown coloured station) and to a vertical corridor a V (from green colored stations). This distance is considered twice to reflect the access and egress leg in the whole trip and it is divided by the average walking speed to finally estimate the access time through Equation (A6). Access time, A: This component will depend on the spacing of the stops and the transport lines (px, py). We assume that the probability to board on a vertical and horizontal corridor is equal. Therefore, access distance is calculated (see Figure A1b) as one-half of the summation of the average distance to a horizontal corridor (from brown coloured station) and to a vertical corridor (from green colored stations). This distance is considered twice to reflect the access and egress leg in the whole trip and it is divided by the average walking speed to finally estimate the access time through Equation (A6). Average waiting time at stop, W: This variable considers the average total time spent by a user on the stops, waiting for the arrival of the transport vehicle. This network allows covering all the trips of the city by means of a transfer, at most. Therefore, the waiting time is computed by W = p(0)·w(0) + p(1)·w (1), where w(x) and p(x) are the waiting time and the probability to make x transfers respectively. Assuming an average headway in routes, the waiting time for users who do not make transfers is w(0) = /2, while the waiting time for trips with one transfer includes the waiting time in the first stop ( /2) and the corresponding at the transfer stop ( /2). In this last case, we assume that the arrivals of the vehicles of horizontal and vertical lines are not coordinated in the transfer stops. In this way, the average waiting time is defined by Equation (A7), considering = .
= + − × + 4 Transfer time, TTR: By transfer time we will understand the time needed to travel on foot from the loading area where the user gets off the first vehicle to the new loading area where the departing bus will stop, to continue his/her journey. The waiting time at the transfer stop is already considered in variable W. The transfer stops will be considered to have a proper design so that the user must walk a distance δ at speed In this way, the transfer time can be easily calculated by Equation (A8), considering that only a fraction p(1) of all users will make a transfer to complete their trip.
In-vehicle travel time, IVTT: As we assume the L1 distance metric, the expected distance between the origin and the destination point (O1 and D1 respectively) of an average trip will be calculated independently in the two Cartesian components: (dIVVT)x and (dIVVT)y. Average waiting time at stop, W: This variable considers the average total time spent by a user on the stops, waiting for the arrival of the transport vehicle. This network allows covering all the trips of the city by means of a transfer, at most. Therefore, the waiting time is computed by W = p(0)·w(0) + p(1)·w(1), where w(x) and p(x) are the waiting time and the probability to make x transfers respectively. Assuming an average headway H in routes, the waiting time for users who do not make transfers is w(0) = H/2, while the waiting time for trips with one transfer includes the waiting time in the first stop ( H/2) and the corresponding at the transfer stop ( H/2). In this last case, we assume that the arrivals of the vehicles of horizontal and vertical lines are not coordinated in the transfer stops. In this way, the average waiting time is defined by Equation (A7) Transfer time, T TR : By transfer time we will understand the time needed to travel on foot from the loading area where the user gets off the first vehicle to the new loading area where the departing bus will stop, to continue his/her journey. The waiting time at the transfer stop is already considered in variable W. The transfer stops will be considered to have a proper design so that the user must walk a distance δ at speed v w In this way, the transfer time can be easily calculated by Equation (A8), considering that only a fraction p(1) of all users will make a transfer to complete their trip.
− p x sD y + p y sD x − p x p y s 2 D x D y (A8) In-vehicle travel time, IVTT: As we assume the L1 distance metric, the expected distance between the origin and the destination point (O1 and D1 respectively) of an average trip will be calculated independently in the two Cartesian components: (d IVVT ) x and (d IVVT ) y . For the sake of simplicity, we will only focus on the x-component (d IVVT ) x . First of all, we assume that the location of point O1 is fixed at x = x O and its distance to the destination D1 is denoted by (d IVVT ) x,O . The distance between points O1 and D1, given that the point O is located at is the probability that point D1 falls on the left side of the origin O1 and d xD<xO the distance between these two points when x D < x O . As origins and destinations of users are uniformly distributed throughout the city, the probability p(x D < x O ) can be calculated as the ratio between the effective length of the segment in which D1 can be located on the left side of O1 (x O /D x ), while the distance d xD<xO is calculated as one half of the effective length of the segment in which D1 can be located (x O /2). The probability corresponding to the situation that D1 is located on the right side of O can be estimated in a similar way, obtaining the following results: With these results we can obtain (d IVTT ) x,o = However, the origin point O1 can be fixed at any point on the x-axis and the destination point can be located on both, its right and left sides. The expected distance in the component x of the evenly distributed points O1 and D1 can be estimated by Similarly, the distance for the vertical component is estimated by (d IVTT ) y = D y /3. Finally, the travel time inside the vehicle can be estimated by means of Equation (A9), where v c ' is the gross commercial speed of public transport vehicles. The gross commercial speed is equal for all bus powertrains and can be calculated by Equation (11a,b) letting θ S = θ N = θ E = θ W = T c,S = T c,N = T c,E = T c,W =l ox = d v = l oy = d h = 0.
Occupancy, O: We aimed at calculating the section y = y that presents a maximum vehicle occupancy in a vertical route like the one depicted in Figure A1. It is assumed that one half of the users are going to complete their trip by accessing a vertical line (although later they can make a transfer to a horizontal line), and the other half will do in a downstream direction. Therefore, the total flow of passengers downwards throughout the width of the city is Q V = Λ(1 + p(1))/4. The total downward total vertical flow can be calculated by the former flow Q V divided by the total number of lines, i.e., q V,sx = Q V ·s x /D x , Finally, it is necessary to determine the passenger flow crossing a generic section y = y (discontinuous red line in Figure A1). For a user crossing the section y = y , the coordinate of its origin (y O ) must satisfy y O > y , while the coordinate of its destination (y D ) must verify y D < y . Since the origins and destinations of users are evenly distributed along D y , the probability that the origin is located in sections y O > y is equivalent to p(y O > y ) = D y −y D y . Similarly, the probability that its destination fulfills y D < y is p(y D < y ) = y D y . Finally, the probability that a vertical trip downwards will cross the section y = y is given by p( y O > y |y D < y ) = D y −y D y y D y . Therefore, the passenger flow at one route at section y = y is defined by q V,s x (y ) = Q V s x D y −y D x D y y D y . It can be demonstrated that the maximal flow is obtained when y = D y /2, obtaining y O > D y /2 y D < D y /2 = 1 4 . In this way, the value of the occupancy at section y = D y /2 in a vertical line O y can be estimated by Equation (A10), as the product of the flow q V,s x (y ) and the time headway H y . This occupancy must