Estimating fuel consumption using regression and machine learning

STOCKHOLM SWEDEN 2018 ,

Estimating fuel consumption using regression and machine learning

LUKAS EKSTRÖM

KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF ENGINEERING SCIENCES

using regression and machine learning

LUKAS EKSTRÖM

Degree Projects Systems Engineering (30 ECTS credits) Degree Programme in Aerospace Engineering (120 credits) KTH Royal Institute of Technology year 2018

Supervisor at Scania: Henrik Wentzel

Supervisor at KTH: Johan Karlsson

Examiner at KTH: Johan Karlsson

TRITA-SCI-GRU 2018:378 MAT-E 2018:80

Royal Institute of Technology School of Engineering Sciences KTH SCI

SE-100 44 Stockholm, Sweden

URL: www.kth.se/sci

This thesis focuses on investigating the usage of statistical models for esti- mating fuel consumption of heavy duty vehicles. Several statistical models are assessed, along with machine learning using artificial neural networks.

Data recorded by sensors on board trucks in the EU describe the oper- ational usage of the vehicle. The usage of this data for estimating the fuel consumption is assessed, and several variables originating from the opera- tional data is modelled and tested as possible input parameters.

The estimation model for real world fuel consumption uses 8 parameters describing the operational usage of the vehicles, and 8 parameters describing the vehicles themselves. The operational parameters describe the average speed, topography, variation of speed, idling, and more. This model has an average relative error of 5.75%, with a prediction error less than 11.14% for 95% of all tested vehicles.

When only vehicle parameters are considered, it is possible to make predictions with an average relative error of 9.30%, with a prediction error less than 19.50% for 95% of all tested vehicles.

Furthermore, a computer software called Vehicle Energy Consumption Calculation tool(VECTO) must be used to simulate the fuel consumption for all heavy duty vehicles, according to legislation by the EU. Running VECTO is a slow process, and this thesis also investigates how well statistical models can be used to quickly estimate the VECTO fuel consumption. The model estimates VECTO fuel consumption with an average relative error of 0.32%

and with a prediction error less than 0.65% for 95% of all tested vehicles.

I

gression och maskininl¨ arning

Sammanfattning

Denna rapport fokuserar p˚ a att unders¨ oka anv¨ andningen av statistiska mod- eller f¨ or att uppskatta br¨ anslef¨ orbrukningen hos tunga fordon. Flera statis- tiska modeller utv¨ arderas, tillsammans med maskinl¨ arning med artificiella neurala n¨ atverk.

Data som registreras av sensorer ombord p˚ a Scania-lastbilar i EU beskriver fordonets drift. Anv¨ andningen av dessa data f¨ or att uppskatta br¨ anslef¨ orbrukningen unders¨ oks och flera variabler som kommer fr˚ an operativa data modelleras och testas som m¨ ojliga in-parametrar.

Uppskattningsmodellen f¨ or den verkliga br¨ anslef¨ orbrukningen anv¨ ander 8 parametrar som beskriver anv¨ andningen av fordonet och 8 parametrar som beskriver sj¨ alva fordonet. Bland annat beskrivs medelhastighet, topografi, hastighetsvariation, andel tomg˚ ang. Denna modell har ett genomsnittligt relativfel p˚ a 5.75 %, med ett skattningsfel mindre ¨ an 11.14% f¨ or 95% av de de fordon som testats.

Om endast fordonsparametrar beaktas som in-parametrar ¨ ar det m¨ ojligt att g¨ ora skattningar med ett genomsnittligt relativfel p˚ a 9.30 %, med ett skattningsfel mindre ¨ an 19.50% f¨ or 95% av de de fordon som testats. Ett da- torprogram kallat VECTO m˚ aste anv¨ andas f¨ or att simulera br¨ anslef¨ orbrukningen f¨ or alla tunga fordon enligt EU-lagstiftning. Att k¨ ora VECTO ¨ ar en tid- skr¨ avande process, och denna rapport unders¨ oker ocks˚ a hur v¨ al statistiska modeller kan anv¨ andas f¨ or att snabbt uppskatta VECTO-br¨ anslef¨ orbrukningen.

Modellen uppskattar VECTO-br¨ anslef¨ orbrukningen med ett genomsnittligt relativfel p˚ a 0.32% och med ett skattningsfel mindre ¨ an 0.65% f¨ or 95% av de de fordon som testats.

II

I would like to express my gratitude to Henrik Wentzel, Johan Karlsson, Marcus Forslin, Antonius Kies, and G¨ oran Svensson, for the help and guid- ance they have given me.

III

Contents IV

1 Introduction 1

2 Theory 3

2.1 Preliminary math . . . . 3

2.1.1 Basic statistical measures . . . . 3

2.1.2 Normalizing . . . . 4

2.1.3 Standardized euclidean distance . . . . 4

2.1.4 Prediction error . . . . 5

2.1.5 Categorical variables . . . . 6

2.2 Regression . . . . 7

2.2.1 Linear regression . . . . 8

2.2.2 KNN regression . . . . 9

2.2.3 Multivariate adaptive regression splines . . . . 10

2.3 Artificial neural network . . . . 16

2.3.1 Artificial neuron . . . . 17

2.3.2 Artificial neural network . . . . 18

2.3.3 Input . . . . 19

2.3.4 Output . . . . 19

2.3.5 Learning . . . . 20

2.3.6 Batch learning . . . . 25

2.3.7 Evaluating observations . . . . 25

3 Data 26 3.1 Vehicle parameters . . . . 26

3.2 Operational parameters . . . . 28

3.2.1 Bin data . . . . 29

3.2.2 Considered variables . . . . 30

3.3 Output data . . . . 32

3.4 Data gathering . . . . 32

3.4.1 Estimation scenario A . . . . 33

3.4.2 Estimation scenario B and C . . . . 33

IV

4.2 Models . . . . 35

4.2.1 Linear regression . . . . 35

4.2.2 KNN regression . . . . 35

4.2.3 KNN Regression with inverse distance weighting . . . 36

4.2.4 Multivariate adaptive regression splines (MARS) . . . 36

4.2.5 Artificial neural network . . . . 37

4.3 Choice of input parameters . . . . 38

4.3.1 Estimation scenario A - VECTO fuel consumption . . 38

4.4 Estimation scenario B and C - Real world fuel consumption . 40 4.5 Selection of model . . . . 42

5 Results 43 5.1 Estimation scenario A - VECTO fuel consumption . . . . 43

5.1.1 Selected model . . . . 43

5.2 Estimation scenario B - Real world fuel consumption using only vehicle parameters . . . . 46

5.2.1 Selected model . . . . 46

5.3 Estimation scenario C - Real world fuel consumption using vehicle and operational parameters . . . . 48

5.3.1 Selected model . . . . 48

5.4 Results summary . . . . 50

5.4.1 Comparison of estimation scenario B and C . . . . 50

5.4.2 Variable influence . . . . 51

6 Discussion and conclusion 52 6.1 Discussion . . . . 52

6.2 Conclusion . . . . 53

References 54 A Glossary 55 A.1 Acronyms . . . . 55

A.2 Terms . . . . 56

B Fundamental statistics 59

C Assumed values for bin data 60

D Additional tables 63

E Activation functions 66

V

F.2 Estimation scenario B . . . . 84 F.3 Estimation scenario C . . . . 96

G Additional figures 109

H Code 110

VI

Chapter 1

Introduction

This report is the result of a thesis made at Scania AB, in cooperation with Royal Institute of Technology (KTH). The thesis is being conducted in the group YDMC at Scania in S¨ odert¨ alje, Sweden, which is specialized in analy- sis of fuel consumption, and by extent, CO

emissions. It is of great interest to be able to quickly and accurately estimate the fuel consumption (FC) and CO

emissions for vehicles. Current methods for estimation relies on sim- ulation, which is too slow and too inaccurate. More advanced simulations require several minutes to complete, which is too slow for being implemented into a sales process, and are still not very accurate. Furthermore, a new leg- islation by the European union require all new heavy duty vehicles (HDVs) to be simulated using a software called Vehicle Energy Consumption Cal- culation Tool (VECTO). The simulation tool is used for certification of the vehicles. Therefore, it is also of great interest to be able to quickly deliver and estimate the VECTO certified CO

to haulage contractors. This paper assesses statistical models for predicting fuel consumption. Several regres- sional models are created and tested. Furthermore, the use of Artificial intelligence (AI) using artificial neural networks are also tested. The paper also investigates the usage of Operational data - data uploaded from vehicles in use - and their potential role in estimating the real world FC.

When VECTO is used for certification purposes, simulations are per-

formed for a set of predefined transport missions, each with a set of prede-

fined transport payloads. Since the driving patterns are predefined, the FC

obtained through VECTO, for a given predefined transport mission, can be

regarded as a function only of the vehicle itself. This thesis focuses on the

transport mission called long haulage with reference payload, which simu-

lates transportation across a large distance, with few stops and dominantly

motorway velocity. A previous thesis, [1] by Max af Klintberg, has assessed

how several regressional models can be used to evaluate the VECTO FC

for the long haulage cycle. That thesis has heavily influenced the methods

assessed in this thesis for estimating VECTO FC. New prediction models

Figure 1.1: The three estimation scenarios

are tested and their performance compared.

The main objective of this thesis is to investigate how well prediction models can be used to estimate fuel consumption. Estimation of two differ- ent values of fuel consumption - VECTO and real world - are to be assessed.

The VECTO fuel consumption describes the average fuel consumption dur- ing one simulated transport mission. The real world fuel consumption de- scribes the average fuel consumption in a longer term perspective.

The input for the models should be a set of variables that can be de- rived from vehicle specification and/or operational data. The estimation of VECTO is done solely based on vehicle data, and does not consider any operational data. Estimation of real world fuel consumption considers both vehicle and operational data. A scenario when the real world FC is estimated solely based on vehicle parameters is also assessed.

Consequentially, three estimation scenarios can be introduced, according to

• Estimation scenario A - Estimation of VECTO fuel consumption with vehicle data as input.

• Estimation scenario B - Estimation of real world fuel consumption with vehicle data as input.

• Estimation scenario C - Estimation of real world fuel consumption with vehicle data and operational data as input.

The three estimation scenarios are also illustrated in Figure 1.1.

Chapter 2

Theory

2.1 Preliminary math

This section introduces mathematical concepts which form a basis for fol- lowing sections. Methods for describing properties in a quantifiable manner, measuring similarity between observations, and for modelling estimation er- ror, are explained.

2.1.1 Basic statistical measures

If X

, X

, ..., X

are random variables, and they form a vector

X = X

X

... X

 . (2.1)

Let x

be the ith independent observation on the jth random variable.

Assume that m ∈ N independent observations have been drawn for all n ∈ N random variables. These observations then form a data set X = {x

∈ R

, i = 1, ..., m}. Further, let ˆ x

be the sample mean and s

be the sample standard deviation of the jth random variable. Obtained as described in Appendix B.

The sample covariance q

of the j:th and k:th random variables is de- fined as

q

= 1 m − 1

m

X

i=1

(x

− ˆ x

)(x

− ˆ x

). (2.2) The pearson correlation coefficient r can then be computed using

r

= q

s

s

. (2.3)

Pearson correlation coefficient takes a value in the range [−1, 1] where 1

is positive linear correlation, -1 is negative linear correlation and 0 is no

correlation. The value r

is called the coefficient of determination, which

takes a value in the range [0, 1] where 1 is linear correlation and 0 is no

correlation.

2.1.2 Normalizing

Consider a random variable x

in a sample, that has a sample mean ˆ x

and a sample standard deviation s

, obtained as in subsection 2.1.1. The sample can be normalized to have a new sample mean ˆ x

and sample standard deviation s

. This is done by modifying all observations x

according to

x

:= x

− s

ˆ

x

s

+ ˆ x

. (2.4) 2.1.3 Standardized euclidean distance

For comparing similarity of quantitative components of observations, a mea- sure of distance between the observations are to be obtained. A formulation of the distance is derived below.

According to [2], a metric on the set X d(x

, x

) is a function such that for all x

, x

∈ X , the following holds:

1. d(x

, x

) ≥ 0

2. d(x

, x

) = 0 if, and only if x

= x

3. d(x

, x

) = d(x

, x

)

4. d(x

, x

) ≤ d(x

, x

) + d(x

, x

).

A metric can then be considered to be the distance between observations.

There are a lot of ways to define a metric. However, the following metrics are relevant for this thesis.

The euclidean distance between x

∈ R

and x

∈ R

is the L

-norm (also called the euclidean norm) of the difference x

− x

, i.e. ||x

− x

||

. The L

-norm of a vector v ∈ R

is defined as

||v||

= v u u t

n

X

j=1

v

. (2.5)

Therefore, the euclidean distance d

(x

, x

) is defined as:

d

(x

, x

) = ||x

− x

||

= v u u t

n

X

j=1

(x

− x

)

. (2.6)

Problems can arise when the euclidean distance is used and the com-

ponents of the vectors vary in size, and larger components can come to

dominate the value of the norm. The standardized euclidean distance nor-

malizes the components in X and makes them dimensionless. If ˆx

and s

are the sample mean and sample standard deviation of component j in the data set X , the standardized euclidean distance d

(x

, x

) is defined as

d

(x

, x

) = v u u t

n

X

j=1

 x

− s

ˆ

x

− x

− s

ˆ x



. (2.7)

2.1.4 Prediction error

Consider a set of m observations for one random variable y = {y

, i = 1, ..., m} ,

and a set of predictions for the same variable ˆ

y = {ˆ y

, i = 1, ..., m} . The prediction error for each observation is defined as

ˆ

e

= ˆ y

− y

(2.8)

and the relative prediction error is ˆ

e

= e ˆ

y

= y ˆ

− y

y

. (2.9)

The root-mean-squared error (RMSE) is defined as

RMSE(y) = v u u t

1 m

m

X

i=1

ˆ

e

. (2.10)

The root-mean-squared for the relative errors is

RMSE

(y) = v u u t

1 m

m

X

i=1

ˆ

e

. (2.11)

Quantiles

The performance of an estimation model can also be measures by quantiles of the prediction error. The quantiles are associated with a probability p, which creates a threshold for the prediction error, given that probability p.

This is useful for measuring the performance regarding outliers. This thesis

uses quantiles for p = 95%, 99%, 99.5%.

2.1.5 Categorical variables

Categorical variables are variables that is set to one of a finite set of possible values. They are not trivially quantifiable, i.e. it is not clear how values for one categorical variable relate to each other, other than being equal or not.

There are several different ways to implement categorical variables into a quantitative model.

Dummy variables

Assume that a categorical variable x

is given which can take on one value from {v

, ..., v

}. This categorical variable can be represented by a series of binary variables y

, j = 1, ..., n − 1 [3], such that

y

(x

) = 0 if x

= v

or x 6= v

1 if x

= v

. (2.12)

The result is a n − 1 binary variables. For example, if x is set to the first value v

, the resulting vector is

v|

= 0 0 ... 0 . If x

is set to the second or third value:

v|

= 1 0 ... 0

v|

= 0 1 ... 0 . Exclusivity

Assume that a categorical variable x

is given, which can take on one value from {v

, ..., v

}, where x

is one component of vectors x ∈ X . Subsets can be created for every value {v

, ..., v

} and X can then be split into these subsets depending on the categorical variable. This would result in n different data sets.

This effectively removes the categorical variable x

from the data sets, and the resulting data sets would be on the form Y = {x

∈ R

}. The drawback of this method is that some of the resulting data sets may be small, and using this method on several categorical variables will partition the subsets into even smaller subsets.

Penalized distance

Using penalized distance, d

(·, ·), is a method of measuring distance between

categorical observations. Assume that x

and x

are two observations for

the categorical variable x

. The penalized distance weight w

> 0 is defined for x

. The penalized distance for the two observations are then

d

(x

, x

) = w

I(x

= x

) = w

if x

= x

0 otherwise. (2.13) The concept can be expanded to vectors of categorical variable. As- sume the vectors x

∈ R

and x

∈ R

are observations of the same set of categorical variable, and that penalized distance weights w

, i = 1, ..., n are defined for all variables. The penalized distance between the two categorical observations are then

d

(x

, x

) = d

(x

, x

) + d

(x

, x

) + ... + d

(x

, x

)

= w

I(x

= x

) + w

I(x

= x

) + ... + w

I(x

= x

).

(2.14)

2.2 Regression

Regressional analysis is a branch within statistic that aims to estimate the relationship between variables. A set of independet variables X (input vari- ables) are used to estimate one or more dependent variables Y (output vari- ables).

A regressional f (X) model is constructed using a set o training data, which is a collection of previous observations consisting of both independent and one dependent variable where x

∈ X and y

∈ Y forms a pair (x

, y

).

The regressional model will then be constructed to make estimates ˆ Y of the dependent variable, using the independent variables:

ˆ

y

= E(y

|x

) = f (x

). (2.15) The above equation combined with Equation 2.8 yields the following func- tion for the prediction error

ˆ

e

= f (x

) − y

. (2.16) The aim of regressional analysis is to construct the prediction function f so to minimize the prediction error. Assume a set of m previous observations, and therefore i = 1, ..., m. A measure on the total prediction error can be formulated as

E = 1 m

m

X

i=1

|ˆ e

| (2.17)

or

E = 1 m

m

X

i=1

ˆ

e

. (2.18)

The process of constructing a regressional model can then be formulated as minimize

f

N

X

i=1

(f (x

) − y

)

, (2.19) which is called the method of least squares.

There are different methods for regressional analysis. The ones used in this thesis are linear regression, K nearest neighbors-regression, Multivariate adaptive regression splines and Artificial neural network.

2.2.1 Linear regression

Assume that there are k independent variables. Linear regression assumes that the relationship between X and Y is linear, and that f (X) is a linear function, according to:

f (x

, β) =β

+ β

x

+ β

x

+ ... + +β

x

=β

a

(2.20)

where there a m independent variables and a

= 1 x



∈ R

. The parameter vector β ∈ R

is obtained through

β = arg min

β m

X

i=1

(f (x

, β) − y

)

. (2.21)

Let A ∈ R

be the design matrix, and y a column vector with all the observations for the dependent variable, according to

A =









 a

a

.. . a











=









 1 x

1 x

.. . .. . 1 x











, y =









 y

y

.. . y









 .

The optimization problem in Equation 2.21 can then be written as β = arg min

β

(Aβ − y)

(Aβ − y) = arg min

β

Φ(β) (2.22)

. Simplifications of the objective function Φ yield Φ(β) = (Aβ − y)

(Aβ − y)

=(Aβ)

Aβ − (Aβ)

y − y

(Aβ) + y

y

=β

A

Aβ − 2(Aβ)

y + y

y. (2.23) Differentiate twice:

∂Φ

∂β = 2A

Aβ − 2A

y (2.24)

∂

Φ

∂β

= 2A

A (2.25)

The second derivative is positive semi-definite and therefore Φ is convex.

The optimum to Equation 2.21 and Equation 2.22 is then obtained from

∂Φ

∂β =2A

Aβ − 2A

y = 0

⇒ A

Aβ = A

y. (2.26)

If A

A is invertible, there exists a unique solution

β

= (A

A)

A

y. (2.27) Thus, predictions for new observations x

can be computed with

ˆ

y

= f (x

) = β

a

= β

1 x



.. (2.28)

2.2.2 KNN regression

K nearest neighbors (KNN) is a non-parametric regressional model, where K ∈ N is a parameter for the model. An observation x

is evaluated by measuring the standardized euclidean distance to observations in the training data x

, i = 1, ..., m. The value of the dependent variables are then equal to the average of the K nearest observations, i.e, the K items with the shortest distance.

d

(x

, x

) (2.29) is obtained as explained in Equation 2.7. Penalized distance is used to evaluate the distance over the categorical input variables x

, according to Equation 2.14, which yields

d

(x

, x

). (2.30) The total distance between the observations is then

d

(x

x

) = d

(x

, x

) + d

(x

, x

). (2.31) Inverse distance weighting

KNN regression may also implement Inverse distance weighting (IDW).

The value of a dependent variable o for new observations x

will then be a weighted average for the K nearest observations in the training data:

f (x

) = (

k=1w_k(x0,x_k)o_k PK

k=1wk(x0,xk)

if d

(x

x

) 6= 0∀i o

if d

(x

x

) = 0 for some i.

(2.32)

In the above equation x

, x

, ..., x

are the K nearest observations to x

. The weights will be inversely proportional to the distance, according to

w

(x

, x

) = 1

d

(x

x

)

(2.33) where u ∈ R

is a parameter that can be used to control the magnitude of the distance weighting.

2.2.3 Multivariate adaptive regression splines

Multivariate adaptive regression splines (MARS) is a nonlinear regression algorithm introduced in [4]. The algorithm results in an estimation function that is a finite sum of basis functions B

(·):

f (x

) =

Mf inal

X

M =1

c

B

(x

) (2.34)

where c

, for M = 1, ..., M

, are coefficients. An example of a function fit to data by MARS can be seen in Figure 2.1.

Figure 2.1: The function f (x) = 0.10968+1.1890∗max(0, x−5.5102) created

by MARS to fit the data.

The MARS algorithm consists of two sequential parts

1. The forward stepwise part finds a set of M

basis functions.

2. The backwards stepwise part tries to reduce the set of basis function provided from the forward stepwise part, to increase generalizability.

The number of basis functions in the final set of basis functions M

, where 1 ≤ M

≤ M max, is selected as a trade-of of generalizability and fit to training data.

Basis functions

The original MARS algorithm implements methods to form basis functions of quantitative input variables, as explained in [4]. Futhermore, an approach for managing categorical input variables is given in [5]. These methods are explained here.

The most basic form of basis function considered by MARS is

B

(x

) = 1. (2.35)

All new basis functions are created as a product of a previously created basis function and some function b(x

) of one variable x

. Assume there are M basis functions already created, basis function M + 1 will then be on the form

B

(x

) = B

(x

) · b(x

). (2.36) The function b(x

) can assume two different forms, depending on whether x

is quantitative or categorical.

For quantitative input variables, MARS uses hinge functions with one quantitative variable as a parameter. Assume x

is a quantitative input variable. The hinge functions are on the form

b

(x

) = h(x

) = max(0, x

− c) (2.37) or

b

(x

) = h(x

) = max(0, c − x

) (2.38) where c is a constant, called the knot. The two above equations are consid- ered a mirrored pair.

Assume now that x

is a categorical input variable, which assumes K

different values so that x

∈ C

where C

= C

, ..., C

. Let A

, A

, ...

be subsets of C

. According to the binomial theorem, there are

2

(2.39)

subsets of a set with K

elements, as explained in [6]. MARS considers functions with a binary output value, on the form

b

(x

) = I(x

∈ A

) (2.40)

or

b

(x

) = I(x

∈ A /

) (2.41) where A

is some of the subsets of C

and

I(x

∈ A

) = 1 if x

∈ A

0 otherwise I(x

∈ A /

) = 1 if x

∈ A /

0 otherwise .

Thus, the basis functions created by MARS are a product of 1 and a number of hinge functions defined in Equation 2.37 and Equation 2.38, and categorical functions defined in Equation 2.40 and Equation 2.41. The in- teraction of a basis function is the number of hinge and categorical functions present in a basis function. Some examples of basis functions are

• B(x

) = 1, the most basic basis function considered, which has inter- action = 0.

• B(x

) = h

(x

), where h

(x

) = max(0, x

− c

), where c

is the knot of h

, This basis function has interaction = 1.

• B(x

) = h

(x

) · I(x

∈ A

), where A

is a subset of the the values that the categorical variable x

can assume, h

(x

) = max(0, x

−c

), where c

is the knot of h

. This basis function has interaction = 2.

• B(x

) = h

(x

) · h

(x

), where h

(x

) = max(0, x

− c

), h

(x

) = max(0, x

− c

) where c

and c

are the knots of h

and h

. This basis function has interaction = 2.

The function in Figure 2.1 have two basis functions, B

(x

) = 1 and B

(x

) = max(0, x

− 5.5102). B

implements a hinge function with a knot at 5.5102.

Finding coefficients

For a set of basis functions B = {B

(·), B

(·), ..., B

(·)}. The coefficients c = c

c

...

 present in Equation 2.34 are found using linear regres- sion. The design matrix is obtained by evaluating the basis functions in B for every observation x

in the training data.

A =











B

(x

) B

(x

) ... B

(x

) B

(x

) B

(x

) ... B

(x

)

.. . .. . . .. .. . B

(x

) B

(x

) ... B

(x

)











. (2.42)

The output vector is y = y

y

... y

. This results in a linear problem.

If A

A is invertible, the solution c

, for B is obtained using Equation 2.27.

Forward stepwise part

The forward stepwise part of MARS creates a basis function set B containing M

basis functions. The set B of basis functions increases in size as the algorithm proceeds. At the start of the algorithm the set contains the basis function B

(x

) = 1. New basis function created for MARS are products of one of the basis functions already in B and a hinge function or categorical function. Furthermore, a pair of basis functions are added at a time, i.e.

two basis functions are added simultaneously.

Consider B = {B

(·), B

(·), ..., B

(·)}, i.e. there are currently M basis functions. Assume M < M

. A pair of new basis functions

B

(·) and B

(·)

should now be added. The new pair of basis function are based on the same basis function B

(·) ∈ B and appends a mirrored pair b

(x

) and b

(x

). So that

B

(x

) = B

(x

)b

(x

) B

(x

) = B

(x

)b

(x

)

where x

is one of the variables x

, j = 0, ..., n. If x

is quantitative, the new basis functions take the form

B

(x

) = B

(x

)max(0, x

− c

) B

(x

) = B

(x

)max(0, c

− x

)

for some knot c used for this pair. If x

is categorical, the new basis functions take the form

B

(x

) = B

(x

)I(x

∈ A

) B

(x

) = B

(x

)I(x

∈ A /

) where A

is some of the subsets A

, A

, ... of C

.

Thus, for every new pair B

(·) and B

(·), the forward stepwise algorithm should determine:

• B

(x

) - One of the basis function already in B.

• x

- The variable to use in the mirrored pair b

and b

.

• b

(x

) and b

(x

) - The mirror paired of b(·)-functions to use. If x

is

– quantitative - b

and b

are hinge-functions as described

in Equation 2.37 and Equation 2.38. The knot location c

has

to be determined.

– categorical - b

and b

are categorical functions as de- scribed in Equation 2.40 and Equation 2.41. The set A

has to be determined.

Determining the above values is an iterative process. MARS considers every combination of B

(x

) ∈ B and x

. Furthermore, the algo- rithm considers a range of different ways to formulate b

(x

) and b

(x

):

• If x

is quantitative, then b

(·) are a pair of hinge functions as described in Equation 2.37 and Equation 2.38. There is an infinite number of possible knot locations. However, MARS considers every observation x

, i = 1, ..., n in the training data and as a possible location for the knot.

• If x

is categorical, then b

(·) are a pair of of categorical func- tions as described in Equation 2.40 and Equation 2.41. All possible subsets A

, A

, ... of C

are candidates to form A

.

For all of these combinations, two candidate basis functions ˆ B

and ˆ B

are created:

B ˆ

(x

) = B

(x

)b

(x

) B ˆ

(x

) = B

(x

)b

(x

).

MARS then considers a set consisting of B, ˆ B

and ˆ B

:

B ˆ = {B

, B

, ..., B

, B

, B

} (2.43) and thereafter solves the linear optimization problem with the design matrix provided by Equation 2.42 given ˆ B, which yields the coefficients ˆ

c ∈ R

corresponding to ˆ B. Thereafter the RM SE, with respect to the training data, is calculated for the current solution.

When RMSE has been computed for every combination of the elements in B and b , the pair ˆ B

, ˆ B

associated with the smallest RMSE is selected and added as basis functions to B. Hence, M := M + 2. The procedure repeats untill M ≥ M

.

Backward stepwise part

The algorithm featured in the forward stepwise part is a greedy algorithm and resulting in overfitting of the model given by B, meaning the that model has adapted to noise in the training data and does not generalize well.

Therefore, the backward stepwise part attempts to remove basis functions

from B. Removing these functions will increase the prediction error for the

training data collection. However, the basis function that gives the smallest

increase in RMSE are removed one at a time. The final basis function set B

is selected by a trade-off of generalizability and fit to training data.

The algorithm proceeds as following:

1. The initial data set is B

:= B, obtained from the forward stepwise algorithm.

2. The basis function that can be removed with the smallest increase in LOF is removed from the previous B

. The resulting set of basis function is B

. The process is repeated until there is a set that only contains one basis function. The result is a range of basis function sets B

, B

, ..., B

where N (M ) = M

−M +1 is the number of basis functions in B

.

3. The final set B

is determined using generalized cross-validation (GCV). Let RM SE

be the root-mean-squared error obtained from performing regression with the data set B

on the training data.

Furthermore, MARS defines effective number of parameters (EN OF ) of B

as

EN OF

= N (M ) + ψ N (M ) − 1

2 (2.44)

where ψ = 2.3 is a penalty, which is fixed in this thesis. The GCV for B

is then

GCV

= RMSE

N (M ) (1 − EN OF

/N (M ))

. (2.45) The final set of basis functions is the set that has the lowest GCV:

M

= arg min

M ∈{1,...,Mmax}

GCV

, (2.46)

B

:= B

M

:= N (M

).

Given B

, the corresponding optimal coefficients, c

= c

c

... C



is obtained using linear least squares fitting with the design matrix provided in Equation 2.42.

Optimization

The forward stepwise algorithm may be very computationally demanding.

For every new basis function to add to B, the algorithm will consider every

combination of

• Basis functions already in B.

• Every input variable x

.

• Every point x

in the training data if x

is quantitative, or every subset of C

if x

is categorical.

For combination of these variables, the MARS algorithm solves the linear problem resulting from Equation 2.42. Solving the problem using the normal equation solution as given in Equation 2.27 is computationally fast. The number of columns in the design matrix are proportional to the amount of basis functions in ˆ B, computation time for solving every linear problem increases as B grows.

The computational time for the entire forward stepwise part grows rapidly with M

and the training data size, as described in [7]. Furthermore, the computation time grows exponentially with the length of C

for every cat- egorical variable x

, according to Equation 2.39.

Two methods of reducing the amount of computations for the forward stepwise part are presented in [7]. These include:

• Introducing an interaction limit. This limits the interaction of basis functions to a certain positive integer. Basis functions in B that have interaction equal to the interaction limit will not be considered to form a factor in new basis functions.

• Introducing a cutoff count and a priority queue for which basis func- tions B to consider. The basis functions in B are ordered in a de- scending order according to their latest computed improval of fit - the reduction in LOF. The algorithm then only considers the basis func- tions at the top of the priority queue. The cutoff count is the number of basis functions considered.

• Ignoring categorical parameters with many subsets.

• Only considering some of the points x

for quantitative parameters x

as possible knot location for hinge-functions. The points should be selected so that they are evenly distributed over the range of values of the variable.

2.3 Artificial neural network

Artificial neural networks (ANNs) are systems of computation in the field of

AI that can be used for regression and classification. Inspired by biological

neural networks, which exists in brains, artificial neural networks have a

set of artificial neurons distributed in different layers. There is one input

layer and one output layer. Between them, there can also be a number of

hidden layers. Signals are passed from the input layers through the net- work to the output layer. The neurons are connected by synapses, which all have synaptic weights that signifies how the transported value is amplified between the two neurons connected by the synapse. The neural network can learn by adjusting the synaptic weights to produce certain outputs de- pending on the inputs. This section describes the type of Artificial neural network considered in this thesis.

2.3.1 Artificial neuron

A neuron that has n input synapses receives input values x

, x

, ..., x

and produces one output value. The procedure is described in Figure 2.2.

Figure 2.2: An artificial neuron The following components are used by the neuron:

1. Inputs x

, x

, ..., x

. There is one input from input synapse - one input from every neuron in the previous neurons layer.

2. Synapse weights w

, w

, ..., w

. The weights of the input synapses.

3. Total net input is the weighted sum of all the inputs, according to u =

n

X

i=1

w

x

. (2.47)

4. Output value y. The total net input is transformed using an activation function g to produce the output, according to

y = g(u). (2.48)

The output of the neuron can then be defined as

y(x

, x

, ..., x

) = g

n

X

i=1

w

x

!

. (2.49)

Activation function

The activation function g : R 7→ R, is a function that transforms the total net input. Determining the activation function of an Artificial neural network (ANN) is directly linked to the performance of the network. Non-linear ac- tivation functions give ANN:s their ability to work with non-linear patterns.

Some common activation functions are

• Logistic function g(u) = σ(u) =

• Hyperbolic tangent g(u) = tanh(u) =

• Identity function g(u) = u

• Rectified linear unit (ReLU) g(u) = u for x ≥ 0 0 for x < 0

• Leaky rectified linear unit (L-ReLU) g(u) =

 u for x ≥ 0 0.01u for x < 0 More information about activation functions can be found in Appendix E.

2.3.2 Artificial neural network

The type of ANNs considered in this thesis are called Feed forward artificial neural network. The neurons are organized into layers. All neurons in the hidden layers and output layer are called receiving neurons. All receiving neurons are connected through synapses to every neuron in the previous layer so that the outputs of a layer are the inputs of the next layers, as can be seen in Figure 2.3. Every synapse has a specific weight associated with it. Furthermore, all neurons in the hidden layers and output layer have a transfer function associated with it. Is is possible to have different transfer functions for different neurons in the same ANN, however this thesis only considers one transfer function for all hidden layer neurons and one transfer function for the output layer neurons. This concept is explained further in subsection 2.3.4.

Bias neuron

It may be beneficial to add a constant term to the ANN. This is done by

adding a neuron with a constant output, called bias neuron. All neurons in

hidden layers and output layers are connected to the bias neuron by regular

synapses. This allows the weight of the constant term to be adjusted in the

same way as other weights in the ANN. The output value of the bias neuron

is usually set to y

= 1. An example of an ANN with a bias neuron can be

seen in Figure 2.4.

Figure 2.3: Example of an artificial neural network with 3 input values, 1 output value and one hidden layer with 4 neurons. The arrows represent synapses. Every synapses has an individual weight w associated with it.

2.3.3 Input

The input values of an ANN can be normalized as explained in subsec- tion 2.1.2. According to [8], this increases the computational efficiency of the learning process.

Categorical input variables are managed using the dummy variable ap- proach explained in subsection 2.1.5.

2.3.4 Output

Many activation functions have a limited codomain. For example:

• The codomain of the logistic function σ(u) =

is described by 0 < σ(u) < 1 ∀ u ∈ R.

• The codomain of the hyperbolic tangent function tanh(u) is described by −1 < tanh(u) < 1 ∀ u ∈ R.

The output values of a pattern that the ANN is to learn must be in the

codomain of the transfer functions of the output layer. To achieve this, 2

methods are considered:

Figure 2.4: Example of an artificial neural network that features a bias neuron.

1. Normalizing the output data so it is distributed in a range compatible with the codomain of the output transfer functions. The normalizing is performed as described in subsection 2.1.2.

2. Setting the transfer function of the output layer neurons to some func- tion g(x) that has a codomain that includes the range of the output values. The identity function g(u) = u can conveniently be used, as its codomain is R.

2.3.5 Learning

The learning process of an ANN is done using the backpropagation algorithm (BPA). Using this algorithm, the network adapts to the training data. This is done by adjusting the synaptic weights

Assume an ANN created for a set of data consisting of n

input variables

and n

output variables, and which contains m observations. Let v

∈ R

and o

∈ R

be the input and output vectors for the i:th observation.

There are therefore n

neurons in the input layer and n

neurons in the output layer. There are N

hidden layers. The number of neurons in the hidden layers are n

, 1 ≥ i ≥ N

. Total number of layers is N = N

+ 2.

The number of neurons in all layers are described by the set

n = {n

, n

, n

, ..., n

, n

} . (2.50) The activation functions for the hidden layer neurons and output layer neu- rons are g

(u) and g

(u) respectively. Furthermore, assume g

(u) and g

(u) are C

or a higher level of differentiability.

The BPA can be summarized by the following stages:

1. Forward pass - The procedure of sending values through the ANN to obtain an estimate of the output values and a metric of the prediction error.

2. Backward pass - The process of finding the influence of the synaptic weights on the prediction error.

3. Updating synaptic weights using the Gradient descent algorithm.

For this subsection, the following notation is introduced:

• L

= {2, 3, ..., N − 1} is the collection of hidden layers. L = 1 repre- sents the input layer and L = N represents the output layer.

• n

is the number of neurons in layer L, indexed as the set in Equa- tion 2.50.

• y

is a vector consisting of the current output values of the neurons in layer L. Consequentially, y

is the output value of the j:th neuron in layer L.

• ˆo

are the output values generated by the ANN for the i:th observation.

That is, the estimations of o

made with the input v

.

• u

is the total net input of the j:th neuron in layer L.

• w

is the synaptic weight of the synapse from neuron j in layer L − 1 to neuron k in layer L.

• If the ANN features a bias neuron, w

is the synaptic weight of the

synapse from the bias neuron to neuron k in layer L.

Cost function

A cost function is defined to represent the error of prediction. The BPA is used to minimize the cost function. The cost function used in this thesis is the Quadratic cost function.

Assume that the ANN has been used to predict the i:th observation where the target values are o

∈ R

and the predicted values are ˆ o

∈ R

. The quadratic cost function is defined as

φ(o

, ˆ o

) = 1 2

nout

X

k=1

(ˆ o

− o

)

. (2.51)

Forward pass

The forward pass is the process of sending values through the ANN form the input layer to the output layer. Before the forward pass, an input vector v

is assigned to the input neurons, so that

y

:= v

.

Thereafter, all receiving neuron layers compute their output values y

according to Equation 2.49, sequentially, starting with the first receiving layer, i.e. layer L = 2. At this point, the following holds for the input and output vectors:

ˆ

o

= y

. (2.52)

The estimations for the i:th observation is thereby obtained, and prediction error φ(o

, ˆ o

) can be calculated using Equation 2.51.

Backward pass

The objective of the backward pass is to find the partial derivative

∂φ

∂w

for all synaptic weights. To find these partial derivatives, the chain rule is applied. The derivatives are computed sequentially, starting with the final layer. The formulation of the chain rule product used for all synaptic weights is

∂φ

∂w

= ∂φ

∂y

· ∂y

∂u

· ∂u

∂w

. (2.53)

For the synapses leading into the output layer, L = N , Equation 2.53 yields

∂φ

∂w

= ∂φ

∂y

· ∂y

∂u

· ∂u

∂w

. (2.54)

For the first partial derivative, the following holds, as explained in Equa- tion 2.52:

∂φ

∂y

= ∂φ

∂ ˆ o

. (2.55)

Differentiation yields

∂φ

∂ ˆ o

= ∂

∂ ˆ o

"

1 2

nout

X

k=1

(ˆ o

− o

)

#

= ∂

∂ ˆ o

 1

2 (ˆ o

− o

)



= ˆ o

− o

⇒ ∂φ

∂y

= y

− o

. (2.56)

The second partial derivative is given by differentiating Equation 2.48, using the transfer function g

(u). This yields:

∂y

∂u

= dg

du    

. (2.57)

The third partial derivative

(N ) k

∂w^N_jk

can be obtained by differentiating Equa- tion 2.47, according to

∂u

∂w

= ∂

∂w

"

X

i=1

w

y

#

= y

(2.58)

where y

is the output of the neuron in the preceding layer connected sending its output through the current synapse.

By considering Equation 2.56, Equation 2.57 and Equation 2.58, all par- tial derivatives of Equation 2.54 are known and the procedure can be used to calculate

∂φ

∂w

for L = N .

For partial derivatives of synaptic weights leading into a hidden layer, i.e. for L ∈ L

, Equation 2.59 becomes

∂φ

∂w

= ∂φ

∂y

· ∂y

∂u

· ∂u

∂w

∀ L ∈ L

. (2.59) The second and third partial derivatives of the above product are obtained in a manner consistent with the output layer,

∂y

∂u

= dg

du    

(2.60)

and

∂u

∂w

= ∂

∂w

"

X

i=1

w

y

#

= y

(2.61)

for all L ∈ L

. Because the backward pass starts with the last layer, all partial derivatives for layer L + 1 are known. This can be used to calculate the partial derivative

∂y_k^(L)

∀ L ∈ L

:

∂φ

∂y

=

n^(L+1)

X

l=1

w

· ∂φ

∂u

=

n^(L+1)

X

l=1

w

· ∂φ

∂y

· ∂y

∂u

∀ L ∈ L

, (2.62)

where n

is the number of neurons in layer L + 1. The partial derivatives

∂φ

∂y_l^(L+1)

and

(L+1) l

∂u^(L+1)_l

have previously been obtained and Equation 2.62 can be computed.

By considering Equation 2.62, Equation 2.60 and Equation 2.61, all par- tial derivatives of Equation 2.59 are known and the procedure can be used to calculate

∂φ

∂w

for all L ∈ L

.

If the ANN has a bias neuron, the synaptic weights w

are updated in the same way as the synaptic weights w

.

Gradient descent

After the backward pass is completed, all partial derivatives are known.

Gradient descent is used to update the synaptic weights. The method is explained in detail in [9]. The weights are updated according to

w

:= w

− η∆w

, where ∆w

= ∂φ

∂w

(2.63)

for all L, j and k. The value η is called step size, and is central to the

gradient descent method. Specifically for ANNs, η is usually called learning

rate. The choice of learning rate greatly affects training performance - A

high value may result in divergence and failure of the algorithm, whereas a

low value may result in entrapment in local minima or slow learning.

2.3.6 Batch learning

During training, the observations in the training data may be used multiple times. Each complete cycle of the entire training data through the BPA is called a training epoch.

The BPA can also be defined so that the forward pass and backward pass for several observations in the training data are performed before the synaptic weights are updated. This is called batch learning. Let γ ≥ 1 be the batch size. During batch learning, a sample of γ observations are randomly selected from the training data. The forward pass and backward pass are performed for the observations. The weights are then updated using gradient descent as defined in Equation 2.63, but instead letting ∆w

be the average of

jk

for observations in the batch.

Furthermore, two special cases of batch learning are described in [9].

These are

• Online learning. γ = 1 - The weights are updated after every iteration of forward pass and backward. pass, i.e. m times per training epoch.

• Offline learning. γ = m, where m is the number of samples in the training data - The weights are updated once per training epoch.

2.3.7 Evaluating observations

When an ANN has been trained, observations can quickly be evaluated.

Consider a new observations with input vector v

. An estimation ˆ o

of the

output can be obtained using the forward pass, explained above.

Chapter 3

Data

This chapter describes parameters considered as variables in the models.

The variables introduced here are candidates for being in the final estimation models. Therefore, a variable introduced here will not necessarily be used in the final models.

The data gathering shall result in a database where the considered vari- ables can be obtained using vehicle Chassis number .

The first two sections of this chapter describes variables considered as input for the models, i.e. variables that might be used to predict the fuel consumption. The third section describes the two output variables.

Two categories of variables are considered. Vehicle data and operational data. Operational data variables are not present in estimation of VECTO FC.

3.1 Vehicle parameters

Vehicle parameters are variables that are part of, or can be derived from, the vehicle specification. The variables can be determined from the vehicle itself, regardless of driving patterns and other operational data. Most of the variables can be directly obtained from vehicle specification and vehicle data.

• Wheel configuration c

- A string representing the total number of wheels, number and location of driven wheels and the number of location of supporting wheel axles.

• Chassis adaption c

- Whether the vehicle is a tractor or rigid.

• Rolling resistance coefficient RRC - Ratio between the rolling

resistance force and normal force. Assumed to be temperature- and

speed invariant. Two different values of rolling resistance are consid-

ered: