Combined Rolling Resistance and Road Grade Estimation Based on EV Powertrain Data

Energy consumption prediction is increasingly important for eco-driving, energy management, and charging scheduling of electric vehicles. Detailed knowledge of the rolling resistance and road grade, combined here in a road-resistance profile, improves the accuracy of these predictions. This paper presents a recursive method to identify the position-dependent road-resistance coefficient using GPS position, powertrain power, and vehicle speed. The calculations make explicit assumptions regarding the spatial continuity of both road gradient and rolling resistance by defining road segments. A recursive least-squares method with Gaussian basis functions allows the estimates to be updated whenever a route segment is traversed anew. The method is tested on data gathered by a 12 m battery electric bus. The resulting road-resistance profile shows a strong resemblance to the road slope and captures changes in rolling resistance well, including a dependency on ambient temperature, which is in accordance with literature on tire rolling resistance. Including the resistance profile in a vehicle model reduces the error of the predicted powertrain power by 1.7 percent point compared to a conventional method, without the limitation of requiring a high-resolution digital elevation model.


I. INTRODUCTION
T O REDUCE global climate change and local air pollution, the transportation sector is transitioning to electric mobility [1]. Electric Vehicles (EVs) are already a key solution in this transition, and in the next decade, the total share of EVs is expected to keep growing exponentially [1, p. 75]. The shift towards electric propulsion is also taking place in the public transport sector, where inner-city transport is electrified, mostly by introducing Battery Electric Buses (BEBs) [2]. The shorter routes with relatively low average velocities, together with the need to reduce local air pollution, make the city centers a suitable use case for BEBs [3].
The limited driving range is an important technical challenge still stalling the dominance of EVs with respect to conventional Internal Combustion Engine (ICE) vehicles. For EVs, the available energy stored in the battery is generally far less than is available in the fuel tank of an ICE vehicle. Even though electric powertrains operate more efficiently, the resulting net driving range is lower. A secondary result of the relatively high and constant efficiency of the electric powertrain is that the road loads, i.e., aerodynamic resistance, rolling resistance, and the longitudinal force due to road gradient, dictate a large part of the vehicle's energy consumption. Therefore, the driving range can vary significantly when these resistance forces change from route to route or day to day. Even though progress is made both in increasing the available battery capacity and reducing the energy consumed per driven distance, the variability of the driving range is inherent to the efficient powertrain of the EV. Accurately modeling the energy consumed by an EV powertrain remains an essential requirement to predict the driving range and apply energy-saving control algorithms. Examples of these control strategies include eco-driving [4], look-ahead cruise controllers [5], on-board energy management [6], and eco-routing algorithms [7]. For BEBs, additional uses of energy consumption models are found in solving fleet scheduling problems [2], [8] and exploring the vehicle design space [9], [10], often to minimize Total Cost of Ownership (TCO).
There are generally two types of EV energy consumption models. The first type is a data-driven approach where historical measurements are used to predict the energy consumption of a future trip under similar conditions [11], [12]. Evermore often, machine learning-type methods are used for this purpose. Alternatively, a physics-based approach can be employed, where the longitudinal dynamics of the vehicle are reconstructed to determine the energy that will be required to traverse a route at a certain forward velocity [13]. This method typically relies on models of the different road load forces; aerodynamic resistance, rolling resistance, the longitudinal component of gravity due to road gradient, the acceleration force, and the powertrain losses. Compared to data-based models, physics-based models offer an increased understanding of the energy losses and better extrapolatability to different routes and operating conditions. However, the parameters of the physics-based model can be challenging to obtain.
The rolling resistance is an important physical effect that is challenging to model accurately. Rolling resistance is a function of tire-design parameters, tire inflation pressure, tire temperature, road wetness, and road roughness [14], [15]. A second component strongly influencing the EV power request is the road gradient. In a physics-based prediction, the quality of road This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ gradient data is detrimental to the accuracy of the energy prediction [16]. Road gradient profiles can be obtained from GPS [17], [18], Inertial Measurement Units (IMUs) [19] or Digital Elevation Models (DEMs) [20]. However, each of these methods brings its own disadvantages; GPS is considered inaccurate for elevation measurements, IMU measurements require additional sensors, and DEMs are often of limited spatial resolution or only locally available. Furthermore, DEMs provide no data for parts of the route consisting of tunnels or obscured by overpasses.

A. Literature Overview: Online Rolling Resistance and Road Grade Estimation
Because of the increasing number of sensors in road vehicles, interest arose in the 1990s to identify vehicle parameters online [17], [21], [22]. Most of these earlier studies, both in simulation and experimentally, use an ICE vehicle as a case study [23], [24], [25], [26]. In the presented estimator algorithms, engine torque is often required to be measured. Because direct measurement of this torque is challenging, most studies employ static engine maps to express this signal as a function of speed and fuel rate. In the last decade, more EVs have appeared as a subject in parameter estimation studies [27], [28], [29]. These vehicles offer the advantage of measurement of the motor current, which directly relates to motor torque. Alternatively, the powertrain power can be used [28]. Except for engine or motor torque, additional sensors are often employed to base the estimation on. These include GPS [17], [18], [24], [27], which is often used to assess the road grade, or IMUs [19], [23], [26]. In most studies, the vehicle speed signal is also required and measured using wheel-speed sensors or a tachograph.
Most studies focus on the estimation of either road grade or rolling resistance. Road grade estimators are often combined with simultaneous estimation of the vehicle mass [23], [25], [26], [27], [30]. This is possible because the vehicle mass affects the longitudinal acceleration of the vehicle, which can be measured directly with an IMU or indirectly via the vehicle speed. Online rolling resistance estimators are also described in literature [29], [31], [32] and are sometimes combined with estimations for the aerodynamic coefficient [33], [34] or the vehicle mass [28]. Combined estimation of rolling resistance and road grade is seldomly encountered [24].
While the applied methods can give a good indication of the rolling resistance or road grade of past a route, this gives no direct information on any future conditions the vehicle might encounter. This aspect is left undiscussed in most studies, where results are presented as a function of time or traveled distance with respect to the start of a measurement. Some studies, e.g., [19] and [29], make the first step towards combining data from multiple vehicles by visualizing the estimated parameter as function of position on a map. One other representative study, [35], proposes a detailed method to combine data from a connected fleet of vehicles to estimate rolling resistance and wind influence, yet is only applied to simulated data. To the best of the authors' knowledge, only in [5] multiple measurements from multiple vehicles are combined to arrive at a single estimate for the road gradient. Multiple measurements are averaged using a spatial Kalman filter resulting in a road gradient that compares favorably to a reference profile.

B. Contribution
Although previous studies exist that estimate the rolling resistance or road gradient online, only a few account for the location-dependency of both these parameters and use this to combine data from multiple measurements. Additionally, other studies rarely estimate rolling resistance and road grade combined. This paper presents a novel approach to the problem by defining route segments to include both the continuous nature of road gradient and the spatial discontinuities introduced by the rolling resistance. Recursive least-squares with Gaussian basis functions is applied to iteratively evaluate the measured powertrain power of a battery electric bus that is repeatedly driving the same route. The resulting parameter estimates are subsequently used in a physics-based longitudinal vehicle model to predict the power request of a future vehicle traversing the same route.
This paper is organized as follows. In Section II the physicsbased longitudinal dynamics model that is at the basis of this work is explained. In Section III the theory behind the proposed recursive identification method is detailed. This method is evaluated with measurement data obtained from a 12 m electric bus, as described in Section IV. The results of this experiment are presented in Section V. Lastly, conclusions and future work are described in Section VI.

II. VEHICLE MODEL
This paper follows a physics-based approach to energy consumption prediction. This implies that the longitudinal dynamics of the vehicle are considered, as described in Section II-A. Notably, this model includes an empirically determined map to represent the powertrain losses as discussed in Section II-B and accounts for the temperature dependency of the aerodynamic drag force as discussed in Section II-C. Only the powertrain power will be considered, without auxiliary energy consumption, because the powertrain power directly relates to both rolling resistance and road slope.

A. Longitudinal Vehicle Dynamics
The energy consumption of the electric vehicle is assumed to adhere to the road-load equation [13], which describes the longitudinal dynamics of the vehicle, as seen in Fig. 1, according to: where constants m ef f and m are respectively the effective vehicle mass, which includes rotational inertia of wheels and driveline, and the vehicle mass, c a (T amb ) is the aerodynamic resistance coefficient as function of ambient temperature T amb , c r (s) is the rolling resistance coefficient, α(s) is the road gradient, and g is the gravitational acceleration. F pt (t) represents the driving force applied at the driven wheels, which can be determined according to where P pt (t) is the electric powertrain power supplied by the high-voltage battery, and P loss (T mot (t), v(t)) is the energy lost between wheel and battery as function of motor torque T mot (t) and forward vehicle velocity v(t). In these equations, the forward velocity v(t), powertrain losses P loss (t), and the powertrain power P pt (t) are functions of time t. In contrast, the road gradient α(s), and the rolling resistance c r (s) are assumed functions of the vehicle position s(t), which by itself is a function of time. This time-dependency of s, will be dropped hereafter for brevity. Next, (1) is rewritten as From this equation, we define the road-specific resistance coefficient c road (s) as c road (s) := cos α(s) c r (s) + sin α(s) , which represents the combined effect of rolling resistance and road grade. Using (1) and (2), the variables in left term of (3) are either measured (P pt (t), T mot (t), v(t)) or known constants (g, m, m ef f , c a ). Therefore, it is possible to calculate c road (s) when measurements on these signals are available. There are two conditions under which direct calculation of the road resistance is not possible. Firstly, because (3) does not contain any term for the dissipative friction brakes of the vehicle, c road (s) can only be calculated if these brakes are not engaged. Secondly, because (2) is badly conditioned for low velocities, c road (s) is only calculated in case v(t) > 5 km/h. The constants in (3) are subject to several assumptions. In practice, the vehicle mass m of an electric bus can change due to differences in passenger occupancy. However, modern vehicles are equipped with an Electronically Controlled Air Suspension (ECAS) that enables estimates of the vehicle mass by measuring the air bellow pressure. Given the availability of these sensors, the vehicle mass is assumed to be a known piecewise constant parameter that changes when the vehicle is halted at a bus stop. Furthermore, while the aerodynamic resistance coefficient c a depends on the ambient temperature, as described in Section II-C, these changes are typically slow. Therefore, during a single trip, c a is treated as a constant.

B. Empirical Powertrain Losses Model
The term P loss (T mot (t), v(t)) in (3) represents the combined losses of the inverter, motor, possible driveshafts, the final drive, and bearings. Because modeling these components in detail requires extensive knowledge about the powertrain, an empirical method is followed. To this end, the losses between battery and wheels are measured on a heavy-duty chassis dynamometer for different velocities and torques. The resulting power losses map is visualized in Fig. 2, and is described in more detail in [36]. The map is a function of the wheel torque T whl (t), which here is determined from the motor torque T mot (t) that is reported by the motor controller during the measurements and the known final-drive gear ratio.

C. Temperature-Dependent Aerodynamics
The aerodynamic resistance coefficient c a consist of where c d and A f describe respectively the aerodynamic drag coefficient and the frontal area of the vehicle. Both these parameters are assumed constant, indicating that no side-wind effects are considered here. The expression for the air density ρ is based on [37], which describes a static relation as function of temperature T amb , relative humidity h rel , and the atmospheric pressure p atm . The values for these three parameters are obtained from openly accessible databases of nearby weather stations [38]. The time-resolution of this weather data is 1 hr, and is matched to the start time of each trip.

III. RECURSIVE LEAST-SQUARES IDENTIFICATION
Having defined the road-resistance coefficient c road (s), this section describes the system identification method used to combine multiple observations of this coefficient in a single positiondependent road-resistance profile. Assumptions regarding the position-dependency are summarized in Section III-A and the used recursive least-squares method is detailed in Section III-B. This section thereby describes a baseline method designed to estimate c road (s) from high-frequency data, specifically with the powertrain power-request prediction for electric buses in mind. However, the method might be adapted depending on the available data and the assumption surrounding other intended applications.

A. Route Segmentation
The road resistance c road (s), as defined in (4), combines the influence of the road gradient α(s) and the rolling resistance coefficient c r (s). While α(s) can typically be considered a continuous function of position s, the rolling resistance depends strongly on the road surface [14] and can thus be discontinuous as function of s. Therefore, it is assumed here that roads consist of shorter road segments that each contain a specific road-surface type. Road segments are defined based on the following criteria that can be obtained for any road via OpenStreetMap: r The road name, e.g., 'Station Road' [39]. r The road type, e.g., 'motorway' or 'primary' [40]. r Whether or not the road is part of a junction [41].
Sections of the route with the above three properties in common are defined as part of a single road segment, which is assumed to have a constant rolling resistance c r . This implies continuity of c road (s) within every route segment. Discontinuities in c road (s) are assumed to be possible on points between road segments to represent a possible change in the road surface. If variations in the rolling resistance are present within a single road segment, these might still be captured but only in a continuous way. In Fig. 3 a highway overpass is shown, with colors indicating the different road segments according to the above criteria.
Besides realizing specific assumptions regarding the positiondependency of c road (s), the above segment-based definition of the route offers the advantage that at any location two roads meet, a new road segment is defined. This opens up the ability to combine resistance profiles of vehicles that have several segments, but not the entire route, in common.

B. Recursive Least-Squares With Gaussian Basis Functions
A vehicle is considered that traverses several route segments i = 1, . . ., N s , while measuring P pt (t), T mot (t), v(t), and GPS position. Thereby, c road (s(t)) can be determined through (3) and (2), which is referred to as an observation of c road (s(t)). Whenever a similar vehicle drives the route the next time, a second set of observations of c road (s(t)) can be made. These sequential batches of observations are indicated by k = 1, . . ., N t , where N t is the number of trips. In this section, a single road segment i, defined by the criteria in Section III-A, is considered. To this end, the subscript i will be dropped in the subsequent explanation.
Every time a vehicle traverses this route segment for the k th time, a batch of measurements of the powertrain power where N m,k is the total number of measurements in the k th trip. The total number of measurements N m,k can vary and is a function of the route segment length d, the vehicle speed v(t), and the sampling rate of the measured signals. The measured signals are used to calculate N o,k observations of c road,k (t 1 , . . ., t N o,k ) using (3) and (2). Note that this calculation is only valid if the dissipative brakes of the vehicle are not engaged and v(t) > 5 km/h. These conditions are checked by monitoring the brake-pedal position and forward vehicle velocity, and invalid data is excluded from further processing. Therefore, N o,k ≤ N m,k . The resulting observations for c road (s k ) are gathered in a column c road,k , along with the positions s k at which these are captured; . . .
If the route segment is traversed the next time (k + 1), a new set of observations: c road,k+1 at possibly different locations s i,k+1 is generated. Therefore, the challenge arises to determine an average c road -profile from the spatially non-equidistant observations of c road , and subsequently add new observations to this profile in a recursive manner.

1) Gaussian Basis Functions:
The positions at which the road resistance is measured s k are not equidistant due to the non-constant forward velocity of the vehicle. Additionally, in the next set of observations, the positions s k+1 are likely to differ. Gaussian basis functions are proposed in the position-domain, as shown in Fig. 4, to parameterize the road-resistance profile and achieve a definition for the road resistance that does not depend on sensor sampling time or forward vehicle velocity. Gaussian functions are chosen because these are non-periodic, C 1 continuous, and vanish at infinity, which are all properties that apply to a local road gradient measurement. Additionally, by using the amplitudes of the functions as the unknown parameters, physical interpretation of the parameter values is maintained. Fig. 4 shows that the basis functions are spatially distributed along the length of the road segment. The number of basis functions per road segment N φ depends on the road segment length d and a pre-defined inter-basis function distance d μ ; For the results generated in this paper, d μ = 66 m is taken. The midpoints of the Gaussian curves μ j for the basis functions j = 1, 2, . . ., N φ are chosen to be uniformly distributed over the length of the road segment: where the latter term ensures symmetry with respect to the center of the segment, which is convenient if the segment is traversed in both directions. Lastly, the basis function width σ φ , is constant and is taken σ φ = 44 m in this paper. The values described here result in approximately 1000/d μ ≈ 15 basis functions per kilometer and are sufficiently small to capture any gradient changes in the measurement data, such as those that will be described in Section V-B. Additionally, the fact that σ φ = 2/3d μ ensures that the different basis functions partially overlap, enabling even relatively constant profiles to be approximated.
The Gaussian basis functions are defined as [42] The lower condition in (9) ensures that if d < d μ , e.g., when only one basis function fits the road segment, a matrix of ones J is used and thus c road (s) is assumed constant over the entire road segment. If d ≥ d μ , multiple basis functions are required. A regressor matrix Φ(s k ) is constructed by ordening the N φ basis functions as The regressor matrix is used in a linear regression model defined as where c road,k (s k ) is the data-vector containing the observed values for c road , Φ(s k ) is the previously defined regressor matrix, containing basis functions that are a function of measured positions s k , θ is the column containing the individual amplitudes of the basis functions, and ε k (s k ) is the error term. The positions s k at which the c road -observations are made are based on GPS-measurements, as is explained in Section V-A.
The assumption is made that any error in s k is negligibly small compared to the width of the basis functions, thereby making Φ(s k ) completely known. Also, assuming that the errors ε k (s k ) are zero-mean, the parameter estimates in this least-squares problem can shown to be unbiased.

2) Combining Multiple Observations:
Having defined the regressor matrix Φ, next, parameter estimatesθ are sought that minimize the errors ε k (s k ) in (11) in a least-squares sense for all available N t trips, i.e., A recursive least-squares algorithm is applied [43], by repeatedly evaluating for k = 1, . . ., N t . Whenever new data becomes available in the form of c road,N t +1 (s N t +1 ), (13), (14), and (15) can be re-evaluated to arrive at updated parameter estimatesθ N t +1 . In order to fulfill the excitation condition of the RLS algorithm defined by (13), (14), and (15), the matrix Φ(s k )Φ(s k ) T should be full rank. According to (10), the regressor matrix Φ contains equidistantly spaced Gaussian basis functions. Therefore, insufficient excitation can occur if the positions of the observations s k are such that one or more of the basis functions is not excited. Therefore, the average distance between measurements should be of the same order or smaller than the basis function width σ φ . During the first observation, i.e, k = 1, no prior information is available. Therefore, the following initial conditions are assumed: with a a constant that should be sufficiently large compared to the covariance introduced by measurement data, and I the identity matrix, thereby assuming no prior knowledge. In these equations, the matrix P k can be used to estimate the covariance of the parameter estimates according to [44, p. 74] where the estimated error variance is given bŷ There are several arguments for choosing the above method over other function-estimation methodologies. First of all, the use of a parametric approach with basis functions allows the estimated profile to be stored efficiently, with only 15 parameters per km of road. This computational efficiency allows the method to be applied to a large road network. In contrast, nonparametric approaches can be computationally less efficient. An exemplary study of a nonparametric method is [45] where Gaussian Progress (GP) regression is applied to estimate a spatial disturbance, similar to the spatial road-resistance discussed here, in a repetitive control problem. Even though computationally less attractive, GP regression does not require the basis functions to be explicitly defined, which can be advantageous if the shape of the estimated profile is not known in advance. Secondly, the use of a recursive algorithm gives the ability to conveniently update c road -profile any time new observations become available, i.e., when the same road segment is traversed again, by simply re-evaluating the covariance matrix P k and the parameter estimatesθ according to (13), (14), and (15). This is illustrated in Fig. 5. The figure shows the road-resistance profile, first shown in Fig. 4, gradually evolving as function of the number of sets of observations k.
The above explanation describes the calculations for a single route segment i and is repeated for all segments constituting the route. In the end, the profiles of the individual road segments are aggregated into a total-route profile. Note that the resistance profile is defined as function of segment position s and not as function of traveled distance. As a result, data from vehicles with different start positions can be combined.

IV. BATTERY ELECTRIC BUS EXPERIMENT
The method described in Section III is tested on real-world data. A dedicated experiment is performed by driving a battery electric bus, as described in Section IV-A, repeatedly over a predefined route, as described in Section IV-B. The resulting dataset is summarized in Section IV-C.

A. Vehicle Details
The used vehicle is a 12 m battery electric bus with a curb weight of 14 tonnes and a 160 kW central electric motor. The vehicle was fitted with the factory-default set of sensors which are logged via the CAN-bus. These measurement signals include the powertrain power P pt (t) and motor torque T mot (t) as reported by the powertrain inverter sampled at 20 Hz and the forward vehicle velocity v(t) as measured by ABS-sensors sampled at 50 Hz. The resolution of the sensors reporting these signals is respectively 0.01 kW, 125 Nm, and 0.05 km/h. Additionally, for this test, a GPS sensor was added that reports the GPS coordinates of the vehicle at a sampling frequency of 10 Hz. Further vehicle details are provided in Table I. To simulate different passenger loading conditions, the vehicle was loaded with sandbags resulting in vehicle weights also mentioned in Table I.

B. Route Details
The considered route is 9.5 km long and mainly includes a primary road and a highway. Fig. 6 shows the route details. The route contains three longer, straight sections and is relatively flat except for two motorway links; an elevated highway on-ramp at 4-5 km and a highway off-ramp at 10 km. The figure also shows the road type and junctions as indicated by OpenStreetMap. Based on these, the route is divided into 7 route segments, indicated by i = 1 ,..., 7 . Even though 2 and 3 are both primary roads, these are considered separate route segments because they bear different names in OpenStreetMap. This is a result of the fact that these two roads are maintained by different municipalities, resulting in a difference in road surface between the two. The figure also shows that, based on the OpenStreetMap data, two roundabouts are identified and the two highway on/off-ramps are considered separately from the highway.

C. Resulting Dataset
The instrumented vehicle described in Section IV-A is driven on the route multiple times in the span of approximately one month. The vehicle adheres to standard traffic rules during driving without making any special maneuvers. The resulting dataset describes 12 trips, as detailed in Table II. The trips vary in driver, vehicle load, and weather conditions. For each trip, the calculated air density is included in the last column [37]. Due to  TABLE II  RESULTING REAL-WORLD DATASET INDICATED BY DATE, DRIVER, AVERAGE FORWARD VELOCITY, VEHICLE LOAD INCL. NR. OF PEOPLE ON BOARD, WEATHER  CONDITIONS, AMBIENT TEMPERATURE AND AIR DENSITY faulty measurement equipment, a part of the data is missing in trip #5. Nevertheless, this trip is still used because the presented method accepts non-equidistant data. The signals measured continuously during each of the trips are the forward vehicle velocity v(t), the powertrain power P pt (t), the motor torque T mot (t), the brake pedal position, and the GPS position. The measurements of the powertrain power also serve as a reference for the validation presented in Section V-C. The brake pedal position is measured to distinguish when the dissipative brakes of the vehicle are used.

V. RESULTS
This section summarizes the results of the experiments. Preprocessing of the measurement signals is explained in Section V-A and the first results are highlighted in Section V-B. The resulting road-resistance profile is used in a power-request prediction in Section V-C. Lastly, temperature effects are investigated in Section V-D.

A. Data Pre-Processing and Map Matching
When driving over a road segment, the vehicle provides data in the form of measurements of velocity v(t), powertrain power P pt(t), motor torque T mot (t), and GPS coordinates. The position s i of each measurement is determined by finding a point on the OpenStreetMap route with the smallest Euclidian distance to the respective GPS measurement. Hereafter, the other measurement signals v, P pt and T mot , are interpolated to the GPS time-grid, which has a temporal frequency of 10 Hz. The signals can thereafter be defined as function of position; v(s i ), P pt (s i ), and T mot (s i ), where the positions s i can be non-equidistant. The vehicle acceleration dv(t) dt , which is also required for (3), is calculated before this step, in the time-domain. To this end, the high-sampling rate forward velocity signal is filtered using a Savitzky-Golay filter with an approximate cut-off frequency of 0.6 Hz [46]. Advantages of the Savitzky-Golay filter are that it preserves high-frequency signal content well, at the cost of a limited noise reduction. Furthermore, it can be applied to non-equidistantly sampled data. Note that the filter can be non-causal, as the measurement information is processed in batches after the vehicle has completed the road segment i.
During pre-processing, all recorded CAN data is matched to the position s i along the route segment. As a result, data from all trips can be visualized on a common position-axis, as seen in Fig. 7. These results show that the 12 measured forward velocity profiles are similar for the majority of the route, except route segment i = 3 , which is a primary road where more traffic might be encountered. This is also visible in the acceleration of route segment 3 . In contrast, on route segment 6 , which is a highway, the forward velocity is nearly constant and the vehicle can drive at its maximum speed of 80 km/h. Over the entire route, the powertrain power shows some correlation with the acceleration and often reaches the limits of ± 160 kW.

B. Resulting Route Resistance Profile
The measured signals P pt (s), T mot (s), and v(s) of trips #1,...,#11 are used to calculate separate observations of c road (s) according to (3) and (2). Based on these, the method described in Section III-B is used to determine the parameters estimateŝ θ N t =11 for every road segment i.
After 11 iterations, i.e., N t = 11, the estimated roadresistance profileĉ road is reconstructed according tô and is displayed for the entire route in Fig. 8. Also included in this figure is the road grade. The road grade is obtained from a high-resolution DEM [47], which has been low-pass filtered using a Savitzky-Golay filter with a spatial cut-off frequency of approximately 1/50 m -1 and is differentiated with respect to position using the forward difference method. The results in Fig. 8 show that between 5.5 km and 9.5 km, where the road is nearly level, α(s) ≈ 0, the c road -profile is also nearly constant atĉ road ≈ĉ r ≈ 0.56%, which is a plausible value for the rolling resistance coefficient c r of bus tires on good asphalt [48]. In contrast, in regions of the route that are sloped, e.g., the highway on-ramp at 4.0-5.5 km and the highway  off-ramp at 9.5-10 km, the profile ofĉ road (s) deviates from this average value and shows strong similarity to the road gradient.
The gray dots in Fig. 8 represent the individual observations of c road (s) of all 11 trips. Because the measurement signals P pt (t), T mot (t), and v(t) are measured at a constant sampling frequency, there are more observations of c road (s) on sections of the route where the vehicle speed is low, for instance at the beginning of the route, around s = 4 km, and near the end of the route. However, because of the 1/v(t) term in (2) the accuracy of c road (s) decreases with velocity, resulting in a wider spread of the observations here. This again motivates the condition to not take into account measurements where v(t) < 5 km/h. Additionally, it is worth noting that if a complete route segment is to be analyzed, persistency of excitation is required; the vehicle needs to be driven at a sufficiently high speed without the usage of the dissipative brakes at every position at least once.
Route segment 3 is of special interest. This section of road, which visually would be considered 'flat,' displays gradient oscillations in the order of 0.5% (α ≈ 0.3 • ). These small oscillations are reflected by the c road -profile in this region between 0.8 Fig. 9. Road coefficientĉ road (s) of road segment 3 after N t = 11 trips. and 3.9 km, which shows a strong correlation to the displayed road gradient. This is further highlighted in Fig. 9. In the bottom part of this figure, the resulting estimated c road -profile is visualized with a 0.8% offset and can be seen to coincide with the road gradient. This offset can be interpreted as an estimated rolling resistance coefficientĉ r , 3 ≈ 0.8%, because for small road grades c road ≈ c r + α(s). Note that this rolling resistance value is larger than the earlier established average on road segment 6 ;ĉ r , 6 ;≈ 0.56%. This difference of (0.8 − 0.56)/0.56 = 43% is plausible, because rolling resistances obtained from previously conducted coast-down tests [36] showed a similar difference between smooth asphalt and bad asphalt.

C. Leave-One-Out Cross-Validation
As described in Section V-B, 11 out of 12 trips are used to generate an estimate of the location-dependent road load coefficientĉ road (s), as displayed in Figs. 8 and 9. The 12th trip of the route is reserved for validation purposes. In the following validation,ĉ road (s) will be compared to a reference profile defined byĉ where α(s) is the road gradient based on a high-resolution DEM [47], as also visualised in Fig. 8, andĉ r,ref = 0.72% is an optimal constant rolling resistance, that is determined by minimizing the powertrain-power error. Taking the rolling resistance constant is a reasonable assumption because the road surface of the entire route is indicated as 'asphalt' in OpenStreetMap. Thus, no distinction can be made between the seven road segments based on this information. Bothĉ road (s) andĉ road,ref (s) are visualised in Fig. 10 and are shown to be similar. However, small offsets can be seen, for instance between 5.5 and 9.5 km, whereĉ road,ref (s) is consistently higher. Next, bothĉ road (s) based on trips #1,...,#11 and the reference profileĉ road,ref (s) are used to predict the powertrain power request of trip #12. By substituting (1) in (2) the powertrain power requestP pt (t) can be calculated according tô where v(t) and T mot (t) are the measured velocity and motor torque of the validation trip. Furthermore, γ 1 and γ 2 are scaling factors equal to 1, unless stated otherwise. This results in P pt (ĉ road (s)) andP pt (ĉ road,ref (s)), respectively. Both predictions are visualised in Fig. 11 and are compared to the measured where N m,k is the total number of data points in the k = 12th trip, and the factor 1/160 kW is used to scale the error with respect to the maximum power of the vehicle. The resulting errors for the power predictions of trip #12 are P E RM S (ĉ road (s)) = 6.0% P E RM S (ĉ road,ref (s)) = 8.2% .
Note that the value forĉ r,ref = 0.72% inĉ road,ref (s) has been obtained by minimizing P E RM S (ĉ road,ref (s)) and it is therefore the optimal average constant rolling resistance coefficient for this approach. Nevertheless, the results show that using the estimated road-resistance profileĉ road (s) results in a lower power error, indicating that the rolling resistance is not constant over the entire route. For this particular route, the rolling resistance seems to be slightly higher on route segments 2 and 3 , which are primary roads, and slightly lower on route segment 6 , which is a highway with relatively new, smooth asphalt. The proposed method captures these effects accurately and also includes gradient estimation without requiring a high-resolution DEM.
A leave-one-out cross-validation is performed by repeating the above validation procedure for different estimation and validation dataset combinations. In all cases, 11 trips are used to estimate the position-dependent road load profileĉ road (s) and the remaining trip is used as validation. The order of the   Table III. These clearly indicate that, irrespective of the validation dataset, the error consistently decreases when using the c roadprofile, compared toĉ road,ref (s). On average, the proposed method reduces the RMS power error P E RM S from 8.3% to 6.6%. This is a strong indication that using a constant rolling resistance results in a less accurate power-request prediction even when the optimized value is known. In contrast, the proposed RLS method with Gaussian basis functions captures the position-dependency of the rolling resistance. Furthermore, the estimatedĉ road (s) also includes the effect of road gradient, which is, per definition, position-dependent.
Likely,ĉ road (s) includes other physical effects or modeling errors that are correlated with position apart from rolling resistance and road gradient. One example is the energy lost due to additional tire slip when cornering, which only occurs in the corners of the route [49]. Even though these additional effects are technically not rolling resistance, the fact that additional unmodeled dynamics are included inĉ road (s) can improve the power-request prediction because future vehicles are likely to experience these same effects.
1) Sensitivity to Aerodynamic Errors: The aerodynamic model described in Section II-C does not account for changes in wind direction, wind magnitude, or effects due to the presence of other vehicles. While the former two effects are assumed to be random, the effect of a leading vehicle will generally result in a reduction of aerodynamic drag force experienced by the following vehicle. For buses, a drag force reduction of as much as 30% can be expected, depending on inter-vehicle distance [50]. Therefore, the effects of a deviating aerodynamic model in the validation dataset are evaluated by varying γ 1 in (21). While the procedure to estimateĉ road (s) remains unchanged, changing the parameter γ 1 implies that the evaluation of the validation data is carried out with a slightly deviating aerodynamic model.
The RMS power errors are displayed as function of γ 1 in Fig. 12. The results show that for 0.8 < γ 1 < 1.33 the method using the road resistance profile outperforms the reference method. Therefore, some errors in the aerodynamic model are considered acceptable. However, additional information would be required to correctly model other aerodynamic effects, such as those caused by a preceding vehicle.
2) Sensitivity to Changes in Powertrain Losses: By varying γ 2 in (21) the robustness of the results with respect to changes in the powertrain losses model P loss can be evaluated. From a practical point of view, deviations in P loss could be caused by a vehicle with a different powertrain configuration or a worn powertrain.
The results are displayed in Fig. 12 and show that at γ 2 = 1, the error of the road resistance profileĉ road (s) is in accordance with Table III and that the error increases for values of γ 2 that are much smaller or larger. The error of the reference method, according to (20), remains almost constant, because its rolling resistanceĉ r,ref is based on minimization of the error, thereby largely compensating for the effect of a changing γ 2 . Nevertheless, the proposed method shows a smaller error for 0.66 < γ 2 < 1.63, indicating that the validation is valid for a large range of deviations of the powertrain losses model.

3) Sensitivity to Inter-Basis Function Distance:
The interbasis function distance d μ determines the smoothness of the road resistance profile within a road segment. A sensitivity study is performed to investigate the effect of this parameter on the validation results. Therefore, the RMS power error is visualised in Fig. 12 as function of varying d μ . The result shows that the error of the road resistance profileĉ road (s) generally increases with decreasing d μ . Up until d μ = 310 m, the result is relatively insensitive to changes in d μ . For this specific route, the error increases at d μ ≥ 320 m, because the number of basis functions N φ in segment 7 decreases to 1. According to (9),ĉ road (s) is therefore assumed constant for the entire route segment 7 , which results in a larger error. In contrast, choosing a high number of basis functions, i.e., d μ < 20 m, results in long computation times, and larger uncertainty on the individual parameter estimates θ. The error of the reference methodĉ road,ref (s) does not depend on the basis function definition and is constant.

D. Temperature Effects
Becauseĉ road (s) includes the rolling resistance coefficient c r , it is expected to have similar physical dependencies. After the road-surface type, the temperature has the second largest effect on the rolling resistance of pneumatic tires [28]. To evaluate if this effect is captured by the proposed method, the trips in Table II, are divided into separate sets. Three trips with the lowest ambient temperatures, trips #1, #3, and #11 are combined, as well as the three trips with the highest ambient temperatures; #5, #7, and #8. This results in a low temperature dataset with an average ambient temperature T amb = 14.7 • C and a high temperature dataset with an average ambient temperature of T amb = 22.3 • C. Both these datasets are each used to estimate aĉ road (s) for the entire route, the result of which is visualized in Fig. 13.
The results shown in Fig. 13 indicate that the colder dataset results in an overall higher c road -profile. This is also evident when the profiles are averaged over s, thereby averaging out the influence of the road gradient, resulting in: The results indicate a decrease of 20% between the cold-weather and warm-weather datasets. This is more than the expected decrease based on literature [13], which dictates a rolling resistance decrease of 7% for passenger car tires over the same temperature interval. Note that the difference cannot be ascribed to a change in air density because this is included in the model as described in Section II-C.
Several reasons can explain this exaggerated temperature effect. Firstly, two of the three cold-weather trips are driven on a wet road, which has an increasing effect on the rolling resistance c r . Secondly, two of the three warm-weather trips were recorded as the second trip of that day, implying that the powertrain and tires are likely to have already reached operating temperature. Most other trips were driven at the beginning of the day, without prior heating of the vehicle. This difference in heating of the powertrain components, especially the drive-axle oil, can cause a significant difference in the powertrain losses, which is reflected here as a change inĉ road (s). Nevertheless, the results indicate that there is a strong temperature dependency in c road (s), which reflects the change of the rolling resistance and also that of other physical effects. Including a forgetting factor allows for tracking changing parameters in time. However, the parameters here represent rolling resistance and road grade and are assumed constant for the duration of the 12 tests shown. Including a forgetting factor could allow for capturing long-term temperature effects that affect the rolling resistance. Also, it would allow the algorithm to adapt if ever a road is replaced, resulting in updated rolling resistance estimates.

E. Separating Rolling Resistance and Road Grade
Throughout this work, the rolling resistance and road grade are considered in a combined parameter c road . A possibility exists to separate this parameter into a direction-dependent and a direction-independent part if a route segment is traversed in the opposite direction as well. In this case, the sign of the resistance due to road slope will change with respect to the rolling resistance. This is demonstrated in Fig. 14, where thê c road -profile of route segment 3 , which was driven from west to east, is repeated. An additional road resistance profile is shown based on two trips where the vehicle traversed the segment in the reverse direction; from east to west. Assuming that both lanes of the road are of the same surface type and have the same road slope, the average of these two profilesĉ average (s) represents the direction-independent resistance or the rolling resistance, according to: c average (s) = 0.5 (ĉ road (s) W →E +ĉ road (s) E→W ) ≈ 0.5 (ĉ r (s) +α(s) +ĉ r (s) −α(s)) ≈ĉ r (s) .
The resulting direction-dependent and direction-independent resistance are both displayed in Fig. 14 and indicate that the direction-independent resistance is indeed constant over almost the entire segment, as is assumed for the rolling resistance in Section III-A. Additionally, the direction-dependent resistance shows a strong similarity to the actual road slope. Deviations at the beginning and end of the segment are likely due to a broader median strip at these locations, causing the road slopes of both lanes to no longer be equal.

VI. CONCLUSION
A recursive least-squares method with Gaussian basis functions is proposed to identify the position-dependent road resistance coefficient to accurately predict the powertrain power request of electric vehicles. Route segments are defined to include both the continuous nature of road gradient and the possible discontinuities due to changes in rolling resistance. The method is tested on data gathered by a 12 m battery electric bus, driving the same 9.5 km route 12 times, under varying conditions. The results indicate that including the estimated road-resistance profile improves the accuracy of the power-request prediction from 8.3% to 6.6%, i.e., by 1.7 percent point, with respect to using a high-resolution elevation model and an optimized constant rolling resistance.
Based on the estimated road-resistance profile from 11 trips, it is shown that road gradient features, such as highway on-ramps and off-ramps, are easily recognizable. Furthermore, even small gradient oscillations in the order of 0.5% are captured well. Also, the effect of a changing rolling resistance is recognizable in the results, ranging from 0.8% for route segment 3 to 0.56% for route segment 6 . The road-resistance profile is shown to be strongly temperature dependent, displaying a 20% decrease for a 7.5 • C ambient temperature increase. This effect can partially be explained by the decreased rolling resistance at higher temperatures, yet also other unmodelled effects might be captured here, such as the temperature-dependent powertrain losses.
The used RLS parameter estimation method is computationally cheap and, because of its recursive aspect, allows for the easy addition of new observations of the same route segments. While the results in this paper are based on data from electric buses, the underlying method can, in theory, be used to combine data from electric vehicles in general. In this case, vehicles from which data is combined should have the same expected rolling resistance coefficient, e.g., similar type of tires and suspension. Additionally, the algorithm assumes a mass measurement to be possible and therefore allows for the combination of data from vehicles with different masses. Due to the presence of the mass m in (1), any deviation of the estimated vehicle mass will result in a bias of the estimated road-resistance profile.
Future work remains to prove that the method can be applied to data originating from different vehicles with possibly different powertrain characteristics. Assuming this is feasible, the method described in this paper could be used to estimate the road-specific resistance for an entire route network based on data from a fleet of similar electric vehicles. Secondly, a forgetting factor can be introduced in the recursive-least squares method to empathize the importance of more recent trips, thereby enabling tracking of temporal changes due to weather influence or road maintenance in the road resistance.