Dynamic Modelling and Optimal Control of Autonomous Heavy-duty Vehicles

SECOND CYCLE, 30 CREDITS STOCKHOLM, SWEDEN 2020

Dynamic Modelling and Optimal Control of

Autonomous Heavy­duty Vehicles

KTH Thesis Report

KARTIK SESHADRI CHARI

KTH ROYAL INSTITUTE OF TECHNOLOGY

Kartik Seshadri Chari <kartikc@kth.se>

Systems, Control and Robotics,

School of Electrical and Computer Science KTH Royal Institute of Technology

Place for Project

Department for Autonomous Transport Solutions, EADM ­ SCANIA AB

Södertalje, Sweden

Examiner

Dimos V. Dimarogonas

Division of Decision and Control Systems, School of Electrical and Computer Science KTH Royal Institute of Technology

Internal Supervisor

Lars Lindemann

Division of Decision and Control Systems, KTH Royal Institute of Technology

External Supervisor

Linus Bergfors

Product Owner ­ EADM Scania AB

External Co­Supervisor

Pedro Lima

Developer ­ EADM

Scania AB

Autonomous vehicles have gained much importance over the last decade owing to their promising capabilities like improvement in overall traffic flow, reduction in pollution and elimination of human errors. However, when it comes to long­distance transportation or working in complex isolated environments like mines, various factors such as safety, fuel efficiency, transportation cost, robustness, and accuracy become very critical. This thesis, developed at the Connected and Autonomous Systems department of Scania AB in association with KTH, focuses on addressing the issues related to fuel efficiency, robustness and accuracy of an autonomous heavy­duty truck used for mining applications.

First, in order to improve the state prediction capabilities of the simulation model, a

comparative analysis of two dynamic bicycle models was performed. The first model

used the empirical PAC2002 Magic Formula (MF) tyre model to generate the tyre

forces, and the latter used a piece­wise Linear approximation of the former. On

top of that, in order to account for the non­linearities and time delays in the lateral

direction, the steering dynamic equations were empirically derived and cascaded to

the vehicle model. The fidelity of these models was tested against real experimental

logs, and the best vehicle model was selected by striking a balance between accuracy

and computational efficiency. The Dynamic bicycle model with piece­wise Linear

approximation of tyre forces proved to tick­all­the­boxes by providing accurate state

predictions within the acceptable error range and handling lateral accelerations up to

4 m/s

. Also, this model proved to be six times more computationally efficient than

the industry­standard PAC2002 tyre model.

Double Lane Change (DLC) and Truncated Slalom trajectories with added disturbances in the initial position, heading and velocities. A linear time­varying Spatial error MPC is proposed, which provides a link between spatial­domain and time­domain analysis.

This proposed controller proved to be a perfect balance between fuel efficiency which was achieved by minimising braking and acceleration sequences and offset­free tracking along with ensuring that the truck reached its destination within the stipulated time irrespective of the added disturbances. Lastly, a comparative analysis between various Prediction­Simulation model pairs was made, and the best pair was selected in terms of its robustness to parameter changes, simplicity, computational efficiency and accuracy.

Keywords

Autonomous driving, Dynamic bicycle model, Magic Formula, Model matching, MPC,

Clothoids

Under det senaste årtiondet har utveckling av autonoma fordon blivit allt viktigare på grund av de stora möjligheterna till förbättringar av trafikflöden, minskade utsläpp av föroreningar och eliminering av mänskliga fel. När det gäller långdistanstransporter eller komplexa isolerade miljöer så som gruvor blir faktorer som bränsleeffektivitet, transportkostnad, robusthet och noggrannhet mycket viktiga. Detta examensarbete utvecklat vid avdelningen Connected and Autonomous Systems på Scania i samarbete med KTH fokuserar på frågor gällande bränsleeffektivitet, robusthet och exakthet hos en autonom tung lastbil i gruvmiljö.

För att förbättra simuleringsmodellens tillståndsprediktioner, genomfördes en jämförande analys av två dynamiska fordonsmodeller. Den första modellen använde den empiriska däckmodellen PAC2002 Magic Formula (MF) för att approximera däckkrafterna, och den andra använde en stegvis linjär approximation av samma däckmodell. För att ta hänsyn till ickelinjäriteter och laterala tidsfördröjningar inkluderades empiriskt identifierade styrdynamiksekvationer i fordonsmodellen.

Modellerna verifierades mot verkliga mätdata från fordon. Den bästa fordonsmodellen

valdes genom att hitta en balans mellan noggrannhet och beräkningseffektivitet. Den

Dynamiska fordonsmodellen med stegvis linjär approximation av däckkrafter visade

goda resultat genom att ge noggranna tillståndsprediktioner inom det acceptabla

felområdet och hantera sidoacceleration upp till 4 m/s

. Den här modellen visade

sig också vara sex gånger effektivare än PAC2002­däckmodellen.

i position, riktining och hastighet lades till startpositionen. En MPC med straff på rumslig avvikelse föreslås, vilket ger en länk mellan rumsdomän och tidsdomän. Den föreslagna regleringen visade sig vara en perfekt balans mellan bränsleeffektivitet, genom att minimering av broms­ och accelerationssekvenser, och felminimering samtidigt som lastbilen nådde sin destination inom den föreskrivna tiden oberoende av de extra störningarna. Slutligen gjordes en jämförande analys mellan olika kombinationer av simulerings­ och prediktionsmodell och den bästa kombinationen valdes med avseende på dess robusthet mot parameterändringar, enkelhet, beräkningseffektivitet och noggrannhet.

Nyckelord

Autonom körning, Fordonsmodell, Modellanpassning, MPC, Klothoider

I want to extend my gratitude to the following people who have helped me throughout this work. First, I would like to thank Scania and the Scania Student Intro (SSI) programme for providing me with an opportunity to work on such a challenging yet exciting topic. Special thanks goes to my manager Cirillo Marcello and supervisors Linus Bergfors and Pedro Lima, who supported me in every aspect of my thesis and ensured that I never strayed away from the right path.

Furthermore, I would like to thank my KTH examiner Dimos Dimarogonas for supporting me throughout and helping take the pressure off my shoulder during this trying period. I would not have been able to complete my thesis had it not been for the guidance and useful insights from my KTH supervisor, Lars Lindemann. So, many thanks to him.

I would also like to thank my friends and colleagues for making my life easier and engaging with all those stupid banters, outings and parties. Last but not least, a big thank you to my parents and family for their invaluable support and love, which helped me throughout my journey.

Kartik Seshadri Chari

Nov 2020.

COG Centre of Gravity DOF Degree of Freedom

ICR Instantaneous Centre of Rotation FOS First Order System

SOS Second Order System

ODE Ordinary Differential Equations RMS Root Mean Squared

PID Proportional, Integral and Derivative LQR Linear Quadratic Regulator

MPC Model Predictive Control LTV Linear Time Varying

LTVR Linear Time Varying Robust RHC Receding Horizon Control ARE Algebraic Riccati Equation SLC Single Lane Change

DLC Double Lane Change

1 Introduction 1

1.1 Research Problem . . . . 2

1.2 Purpose . . . . 3

1.3 Objective . . . . 3

1.4 Benefits, Ethics and Sustainability . . . . 3

1.5 Stakeholders . . . . 4

1.6 Delimitations . . . . 4

1.7 Outline . . . . 5

2 Modelling 6 2.1 Literature Review . . . . 6

2.1.1 Related Work . . . . 7

2.1.2 Theoretical Background . . . . 8

2.2 Methods and Methodology . . . . 26

2.2.1 Model Choice . . . . 26

2.2.2 Tyre Models . . . . 28

2.2.3 Reference Pivot Point . . . . 29

2.2.4 Test Scenarios . . . . 29

2.2.5 Model Matching and Parameter Tuning . . . . 31

2.2.6 Interpolation . . . . 32

2.2.7 State Evolution . . . . 33

2.2.8 Evaluation Metrics . . . . 34

2.3 Results and Discussion . . . . 35

2.3.1 Simulation Framework . . . . 35

2.3.2 Results . . . . 39

3 Optimal Control 49

3.1 Literature Review . . . . 49

3.1.1 Related Work . . . . 50

3.1.2 Background Theory . . . . 51

3.2 Methods and Methodology . . . . 56

3.2.1 Reference Generation . . . . 56

3.2.2 Controller Design . . . . 60

3.3 Results and Discussion . . . . 63

3.3.1 Basic Controller Performance . . . . 64

3.3.2 Spatial Error MPC Performance . . . . 67

3.3.3 Spatial Error MPC with Rate constraints . . . . 71

4 Conclusions 72

Introduction

The last decade has been at the forefront of development of technologies for design of autonomous trucks. Coming from small scale in­lab testing to full out on­road testing, autonomous trucks have surely come a long way. Unlike the concept of autonomous cars which is dealing with its entirely different set of problems, the question with autonomous trucks is not “ if ” but “ when ” will we see them on the roads.

For quite a few years, there has been a very fine line between an autonomous and a non­autonomous vehicle, with each vehicle manufacturer setting their own precedent.

In 2018, the Society of Automotive Engineers (SAE) established the J3016TM “Levels of Driving Automation” standard that defines the six levels of driving automation, from no automation to full automation (see fig. 1.0.1).

An interesting observation can be made concerning the current global status of autonomous trucks. It is quite apparent that it is a market who’s importance and scope of benefits are understood by almost everyone in the transportation industry.

This realisation has stimulated a high­speed chase pitting the usual suspects ­ Google and Tesla along with other global tech firms ­ against small startups. Moreover, the winners of this derby, they may be poised to make untold billions; they’ll change the global transportation grid, and they will emerge as the new kings of the road [2].

The levels of automation give us a good standard to assess the situation of a particular

organisation’s product. Both the incumbent automakers and high­profile startups are

Figure 1.0.1: Levels of Driving Automation [1]

pushing pilot vehicles on public roads and testing their technology to improve self­

driving capabilities [3]. SAE Level two features are now commonplace for almost all top competitors in the autonomous vehicle market. The real game begins with level three autonomy wherein all major manufacturers are already in the prototyping or testing phase. Some autonomous trucking firms such as TuSimple, have even gone to the extent of promising first driver­out demonstrations on public highways by 2021 [2].

Naturally, the most vulnerable people to the exponential developments in autonomous trucking are the truck drivers. On one side, the labour departments argue on the loss of hundreds of thousands of jobs globally. On the other hand, the companies project the elimination of human errors, efficiency and sustainability as their point of argument.

The reality is that these jobs certainly could change, especially as companies that employ truck drivers start to look more closely at self­driving technology. While it’s not clear how fast they will get there or what the rules of the road maybe, but the positive technological change is just another sign that another revolution in transportation is on its way [4].

1.1 Research Problem

Unlike cars, trucks are essentially confined to operate in specific geographies

like mines, warehouses or farms [5]. Some of these geographies like mines

are extremely complex and hazardous. Thus, many companies are investing on autonomous technology to improve efficiency and reduce human casualties. However, this comes with its own set of challenges ­ both technological and legislative.

The technological challenges span regular wear­and­tear, unpredictable road and environment conditions, high braking distances and overall cost and efficiency whereas legislative challenges include safety, data protection, liability risks, etc.

This thesis will try to answer the following questions:

• Which is the vehicle model that guarantees high­fidelity in state prediction and can handle lateral accelerations upto 4 m/s

making it suitable to accurately predict the behaviour of the truck?

• Which is the most optimal controller that can ensure smooth and offset­free tracking of the desired trajectory along with robustness to parameter changes and external disturbances?

1.2 Purpose

This thesis report presents the work done on Dynamic Modelling and Optimal Control of autonomous heavy­duty vehicles. It illustrates a comparison between various vehicle models and controller designs and helps identify the most suitable alternative with the help of simulation results.

1.3 Objective

The goal of this thesis is to design a high­fidelity vehicle model for a truck that handles lateral accelerations upto 4 m/s

for better state prediction and simulation. This model would then be used to design an optimal controller that guarantees offset­free trajectory tracking while ensuring smooth and stable driving.

1.4 Benefits, Ethics and Sustainability

Autonomous trucking is a hot topic in the transportation industry for several reasons.

Right from improved safety to increased efficiency to decreased labour costs, the

benefits are revolutionary. A recent study revealed that in 2016, 3,326 people were killed, and around 110,000 people were seriously injured due to accidents involving a truck. Out of these, driver error caused 90% of these accidents [6]. Research suggests that autonomous truck technology can reduce these accidents by 80% and save around

$36 billion annually on casualty costs.

Also, advanced cooperative strategies like Platooning help avoid congestions, regulate the traffic and reduce the fuel consumption by 10% [7]. Reduced braking and accelerations along with decrease in the effect of air drag can be a few factors that increase the overall efficiency of the fleet.

When it comes to ethical considerations, the classical dilemma pertaining to the

”trolley problem” generally pops up. Given a lack of reaction time and limited options, the autonomous system is tasked with either hitting the pedestrians and killing five of them or swerving into the brick wall and killing the passenger. Even though many experts and philosophers have divided opinion on this problem, there is a common consensus that protecting five pedestrians while sacrificing a single passenger is the right choice [8]. This does not please the potential customers who give more priority to saving their own life over others. Another issue of worry is the vulnerability of the system to hacking and what precautions are being taken in order to prevent such problems.

1.5 Stakeholders

Automobile companies venturing into the world of autonomous driving are the main stakeholders for this project. This work is also of great interest for companies and government entities involved with hazardous activities like mining, where safety and labour costs are very critical.

1.6 Delimitations

The second­half of the thesis focuses on implementing an optimal controller for

trajectory tracking. However, due to time limitations, the results were obtained

through a computer simulation rather than real­time implementation on a physical

truck. This could have a significant effect on the robustness and accuracy of the

controller.

1.7 Outline

This report is divided into two parts ­ Modelling (chapter 2) and Optimal Control (chapter 3). Chapter 2 deals with designing and comparing various vehicle models to determine the most accurate model for state prediction and simulation purposes. In chapter 3, various Model Predictive Control formulations are compared and a Spatial­

error MPC is proposed for smooth, accurate and offset­free trajectory tracking.

Each of these chapters is further divided into three sections. The Literature Review section presents the relevant research in the domain and some background theory to support the work. The various methods and the justifications behind each choice are described in the Methodology section. The Results section describes the implementation architecture and presents important results of the work.

Finally, the overall conclusions, limitations, and future scope of this work are discussed

in chapter 4.

Modelling

As discussed earlier, any autonomous system can be broken down into four significant modules viz. Perception, Estimation, Motion Planning, and Control. Of these, all except Perception rely heavily on a system model. As this work deals with autonomous trucks, the system in focus is the vehicle model.

Within the Estimation module, the vehicle model helps estimate the vehicle parameters and predict future states and model disturbances. Motion Planning primarily deals with avoiding obstacles and generating a smooth and continuous trajectory. Here, knowing how the vehicle reacts to specific scenarios helps determine the optimal path [9]. Now, needless to say, a vehicle model is the cornerstone of stability analysis and reference tracking. Thus, right from basic control techniques like PID control or feed­forward control to more advanced control techniques like Model Predictive Control (MPC), the vehicle model hugely influences the tracking performance.

Additionally, in most applications, it is easier and safer to test the system in simulation before implementing it on the physical system. This demands a simulation model with very high fidelity, i.e., it is as closest to the physical model as possible. Thus, Modelling forms a crucial part of this work.

2.1 Literature Review

The domain of vehicle modelling is very vast and intricate. Depending upon the application area, there are several model types and modelling methods to choose from.

Section 2.1.1 provides a brief review of the work being carried out in this area and

section 2.1.2.1 defines the coordinate frames of reference that are followed throughout the text. Finally, a short description of the various concepts used in this work has been presented in section 2.1.2.

2.1.1 Related Work

In this work, a vehicle model with high fidelity, i.e. a prediction or simulation model that matched the experimental data, had to be derived for use in Motion Planning and Control. On these lines, in [10], Hebib et al. described the Volvo Transport Model (VTM) in which the chassis, wheels, tyres, and suspension were modelled in Simulink and simulated on a multi­body simulation environment Simscape Multibody™. This model was a highly complex and accurate model that enabled system dynamic analysis with bodies, joints, and forces.

In [11], the graphical object­oriented programming written in Dymola using Modelica was studied against Modelica’s code­based approach by the author. He found that the code based approach was more readable and structured but involved a very high workload as each subsystem was modelled via the mathematical equations of motion.

As the application domain of my work is real­time, the model should not be computationally too expensive (i.e., its update frequency should be faster than the system’s refresh rate). Furthermore, the focus should be on the vehicle state analysis instead of the system analysis.

Along the same lines, in [12], Rajamani described the simplest vehicle model that could be used for the vehicle state analysis ­ the Kinematic Bicycle model. It was highly useful for control applications dealing with low velocities but faded out when the vehicle velocity was high as there was a significant slip between the road and tyre. An improved version of this model, called the Dynamic Bicycle model, was explained in detail in [12, p. 27]. This model considered the external forces and yaw­moments that acted on the vehicle and, thus, proved to be better at calculating the slip. In [13], Keviczky et al.

showed the equations of motion w.r.t. the centre­of­gravity (COG) and transformed all the forces in terms of the body­fixed frame. (Refer to subsection 2.1.2.1 for more info.)

Out of all the forces, the tyre forces are the most significant but highly non­linear [12,

13]. Thus, various attempts have been made to develop a tyre model that can reflect the non­linearities in the simplest and most acceptable manner. In [12], Rajamani observed that the tyre forces were approximately linear for smaller values of slip and thus suggested a linear tyre force model. In [14], Jacob presented the Brush model approach, which described the physics behind the generation of forces using the elastic deformations in the rubber of the tyre contact patch. Finally, around 1988, Hans Pacejka proposed an empirical tyre model which went on to become the most widely used tyre model [15]. He observed the sigmoidal shape of the steady­state tyre force response and developed a sinusoidal function to imitate the same. Over the years, this model has been modified multiple times to incorporate the transients and effects of combined slip to form semi­empirical models. PAC2002 is the latest improvement and has become the industry favourite [16].

In [17], Mia Palmqvist proposed the use of two different types of models ­ one for the state prediction and the other for simulation. The former was a basic linear­tyre Dynamic bicycle model, and the latter was a modified version where instead of the bicycle model, all the four wheels were modelled separately which helped to account for the load transfer due to acceleration in a better manner. On top of that, the power­

train had also been modelled which derived the actual torque request that was sent to the wheels. This, in turn, helped to estimate the slip ratio in a better way. In [18], a four­wheel vehicle model was used, and the Longitudinal and lateral tyre forces were estimated via an Extended Kalman Filter (EKF) observer. On the same lines, [19] was inclined towards developing an observer that could reduce the error in the Lateral force prediction using simple measurements and EKF.

2.1.2 Theoretical Background

In order to derive a mathematical model that is as close to the physical truck as possible,

various vehicle models were analysed. As the focus was inclined towards the lateral

motion of the vehicle, the forces governing the yaw moment or steering were given

a priority. Furthermore, several concepts from Control Theory and Linear Algebra

were borrowed to supplement these models. A brief account of the above concepts is

provided below.

2.1.2.1 Co­ordinate Frames of Reference

Whenever position, velocity or acceleration vectors are talked about, they are defined with respect to some axes commonly referred to as Co­ordinate Frames of Reference. These reference frames can be classified into Inertial Frame or Non­

inertial Frame.

• An inertial frame is a frame of reference which is at rest or moving at a constant velocity. According to the Galilean Relativity, if a body with no net force acting on it is considered as the frame of reference, the laws of nature are valid [20]. This means that Newton’s laws of motion are valid in an inertial frame of reference.

• Non­inertial frames of reference are the ones which are accelerating i.e a non­

zero net force is acting on them. In such cases, the Newton’s Laws of Motion are no longer valid due to the introduction of Fictitious forces ­ Coriolis force, Centrifugal force and Euler force.

Y

X

x y x

x

y

l

l

Figure 2.1.1: Reference Co­ordinate Frames

In this work, three different reference frames were considered (see fig.2.1.1). The fixed

inertial frame was denoted by the XY­ axis. The body­fixed frame was a local frame with

origin at the Centre of Gravity (COG). The x­axis denoted the Longitudinal direction

of the vehicle and the y­axis denoted the Lateral direction. The third local frame of

reference was the front wheel axis x

y

.

2.1.2.2 Kinematic Bicycle Model

The simplest model to study is the Kinematic Bicycle Model. As the name suggests, this model has to do with the Kinematics of the vehicle rather than the dynamics. In other words, this model mathematically describes the motion without considering the forces or moments acting on the vehicle.

Y

X

δ

δ

V β l

ψ

l

Figure 2.1.2: Kinematic Bicycle Model (w.r.t COG)

The kinematic equations of motion for the above bicycle model (see fig. 2.1.2) were derived as follows: [12]

X(t) = V cos(β + ψ) ˙ Y (t) = V sin(β + ψ) ˙

ψ(t) = ˙ V cos(β)

l

+ l

[tan(δ

) − tan(δ

)]

β = tan

( l

tan(δ

) + l

tan(δ

) l

+ l

)

(2.1)

For the study case, the vehicle was assumed to be front­wheel steering (i.e. δ

= 0).

Thus, the equation for yaw rate and β simplifies to:

ψ(t) = ˙ V cos(β) tan(δ

) l

+ l

β = tan

( l

tan(δ

) l

+ l

)

(2.2)

Assumptions

• The main assumption of the bicycle model is to represent both the front wheels as a single front wheel and both the rear wheels as a single rear wheel. Also, this means that the inner wheel and outer wheel steering angles are assumed to be the same [12].

• The slip angle is 0. This means that the velocity vectors of the front and rear wheels are in the same direction as the corresponding wheels i.e. they make δ

and δ

angle with the longitudinal axis. This assumption is reasonable for low lateral accelerations, which typically occurs at low speeds [12].

• The inputs to the system are Steering angle (δ

) and the net velocity (V) in the inertial (global) frame. Practically if an acceleration request is sent to the vehicle, another equation has to be included. = ⇒ ˙V = a.

2.1.2.3 Dynamic Bicycle Model

At lower speeds, as it has been already established, there is negligible slip. This makes the Kinematic Bicycle model a good estimate of the vehicle. But, as the lateral acceleration of the vehicle increases, the previously made assumptions don’t hold true.

Thus, a bicycle model which describes the dynamics of the vehicle (see fig. 2.1.3)was derived and analysed.

Y

X

˙x˙

˙ x y

δ

α

V

ψ V

β

ψ, I ˙

, m F

F

F

F

R

R

F

F

R

R

F

F

l

l

Figure 2.1.3: Dynamic Bicycle Model (w.r.t COG)

Longitudinal Vehicle Model

The Longitudinal vehicle model considers all the forces acting along the longitudinal

axis (x­axis). Typically, this model is derived considering a slope with inclination θ as

follows:

m¨ x = −mg sin(θ) + 2F

+ 2F

− 2R

− 2R

− F

+ U

∀

 

 

 

 



 

 

 

 

F

= F

cos(δ) − F

sin(δ)

R

= N ( ˆ C

+ ˆ C

| ˙x| + ˆ C

˙x

) ≈ NC

| ˙x|

F

= 1

2 C

ρA ˙x

≈ C

˙x

(2.3)

where F

is the longitudinal tyre force, R

is the rolling resistance, F

is the air resistance and U

is the input throttle.

For the required application, the road was assumed to be flat, thus making the inclination θ = 0. Also, the eqn 2.3 does not include the lateral force acting on the vehicle due to the yaw rotation. Thus, the final longitudinal equation becomes:

m¨ x = m ˙ y ˙ ψ + 2F

+ 2F

− 2R

− 2R

− F

+ U

∀

 

 

 

 

 



F

= F

cos(δ) − F

sin(δ) R

≈ C

| ˙x|

F

≈ C

˙x

(2.4)

Lateral Vehicle Model

The Lateral vehicle model considers all the forces acting along the lateral axis (y­axis).

The equation of motion at the COG governing the lateral dynamics are as follows:

m¨ y = −m ˙x ˙ψ + 2F

+ 2F

I

ψ = 2l ¨

F

− 2l

F

(2.5)

∀ ß

F

= F

sin(δ) + F

cos(δ)

Thus, the consolidated Dynamic Bicycle model in the state space form is as follows [12,

13]:

X(t) = ˙x cos(ψ) ˙ − ˙y sin(ψ) Y (t) = ˙x sin(ψ) + ˙ ˙ y cos(ψ)

ψ(t) = ˙ d dt ( ¨ ψ)

¨

x = ˙ y ˙ ψ + 2

m F

+ 2

m F

− 2

m R

− 2

m R

− 1

m F

+ a

¨

y = − ˙x ˙ψ + 2

m F

+ 2 m F

ψ = 2 ¨ l

I

F

− 2 l

I

F

(2.6)

Assumptions

• The forces and effects of Suspension, road inclination and load transfer were neglected

• The factor of 2 was included for all longitudinal and lateral tyre forces which accounts for 2 tyres both at the front and the rear.

• The power train was not modelled as it is handled by the low­level controller and the use cases of this model has to do with high­level controller simulation.

2.1.2.4 Tyre Forces

While the vehicle is in motion, the tyres are its only components that are in contact with the road surface. Thus, its natural that most of the forces and moments associated with the road­vehicle interaction are transferred through the tyres [21].

Out of all the forces talked about in the subsection (2.1.2.3), only the aerodynamic

force (F

) is not related to the tyres. Moreover, the order of F

is considered to

be smaller than other forces and thus, is usually neglected. This clearly means that the

tyre forces play the most significant role in modelling the dynamics of the vehicle.

Figure 2.1.4: Forces and Moments acting on the tyre [SAE Recommended Practice J670e, 1976]

The above figure (2.1.4) depicts all the forces and moments acting on the tyre as per the SAE norms. Out of these, the most important ones are discussed as follows:

• Longitudinal Force (F

)

– The longitudinal tyre forces are generated mainly due to the acceleration and braking of the vehicle. These forces are essentially friction forces that act on the tyres [14].

– These forces are affected by the friction coefficient of the road (µ), the Normal force (N) as well as the tyre slip ratio (κ) [12].

– Slip Ratio (κ) : Slip ratio is defined as the ratio of the Longitudinal slip to the longitudinal velocity of the vehicle. (where, Longitudinal slip is the difference between the actual rotational velocity of the wheel and the velocity of the vehicle)

κ = V

− V

max(V

, V

) = ωR − V

max(ωR, V

) ( ∀ − 1 < κ < 1)

– If the wheel is spinning i.e. the vehicle is accelerating, κ is positive as ωR > V

.

– If the wheel is skidding i.e the vehicle is braking, κ is negative as ωR < V

.

– But, if the wheel is locked i.e. the vehicle is in the state of emergency braking, κ is ­1.

– It is intuitive that as F

is a friction force it is dependent on µ and N, but its dependence on the slip ratio is not so intuitive. Refer to Appendix for a detailed explanation.

Figure 2.1.5: Tyre Slip angle [12, 22]

• Lateral Force (F

)

– The lateral tyre forces are generated during cornering and govern the yaw rate and lateral accelerations.

– Experimentally it has been proven that F

depends on the Tyre Side Slip angle (α).

– Side Slip angle (α) : The side slip angle of a tyre is defined as the angle between the orientation of the contact patch and the orientation of wheel velocity vector. (see fig. 2.1.5)

– Mathematically, the side slip angle can be expressed as follows:

α

= δ − θ

∀θ

= tan

( V

+ L

ψ ˙ V

) α

= −θ

∀θ

= tan

( V

− L

ψ ˙

V

)

– The slip angle is formed due to a centrifugal force acting on the tyre tread in

the direction opposite to the wheel yaw and hence, the angles are used with

a negative sign.

Figure 2.1.6: Contact Patch during a Left turn [23]

• Self­aligning Moment (M

)

– While modelling the dynamic bicycle model (2.1.2.3), it was assumed that the tyre forces act at the contact point (centre of the contact patch) but this is not true. The lateral force actually acts at some distance from the centre known as the Pneumatic Trail (t) (see fig. 2.1.6).

– This results in a moment known as the Self­aligning Moment (M

) which is directly proportional to the Pneumatic trail.

– This is a resisting moment which can be experienced at the steering wheel and can prove to be helpful during cornering.

• Rolling Resistance (R

)

– Rolling resistance arises from a Rolling Resistance moment (M

) that acts against the vehicle motion.

– It is mainly generated due to the effect of hysteresis. In other words, every time a tyre rolls, it is deformed due to the vertical load but while coming back to its normal state, a small amount of energy is lost resulting in a moment M

.

– Thus, R

depends linearly on the tyre load and is typically smaller in

magnitude as compared to Static or Kinetic friction [14].

Tyre Models

As the tyre forces play such a significant role in the vehicle dynamics, it is of utmost importance to develop an accurate model describing the tyre behaviours and properties. Depending upon the application area and desired level of accuracy, tyre models have been classified into various categories as follows [21].

• Linear / Simple Tyre Model

– This model mainly focuses on the linear relationship between the tyre slip and the resulting forces. It does not take combined slip into consideration.

– In a simple tyre model, longitudinal tyre force is expressed as

F

= C

∗ κ

where C

is the longitudinal stiffness coefficient and κ is the tyre slip ratio.

– Similarly, the lateral tyre force is expressed as

F

= C

∗ α

where C

is the cornering stiffness and α is the tyre side slip angle.

– It is usually a good approximation for low slip applications or design of vehicle control systems.

• Analytical / Physical Tyre Model

– This model dwells into the physics involved at the tyre­road interaction.

– It provides a mathematical basis for describing the kinematics and dynamics of the tyre contact patch in detail.

– It may vary significantly depending on the area of application and can become very complex.

• Empirical Tyre Model

– This model does not have any physical basis but is rather based on non­

linear mathematical approximations and curve fitting of experimentally generated test data that is inclusive of effects of tyre forces and moments.

– These models require extensive tyre measurements, processing and

parameter optimisations, and tend to be very accurate for vehicle dynamic analyses and simulations [21].

– Typically, a tyre property file (.tir) is generated which is further used to derive the coefficients of the nonlinear force equations.

– Magic Formula, developed by Hans B. Pacejka is a well known model which is widely used for accurate vehicle dynamic simulations.

Figure 2.1.7: Magic Formula Curve Shape in Pure Slip condition [16]

Magic Formula

As mentioned above, the Magic Formula was developed by Hans B. Pacejka in 1980s as part of a cooperation between TU Delft and Volvo Car Corporation. It is a mathematical formula that can describe the basic tyre characteristics for steady state operating conditions.

It has the following 2 forms that depict the basic tyre characteristics as described in [15]:

y(x) = D sin[C tan

{Bx − E(Bx − tan

{Bx})}] (2.7)

y(x) = D cos[C tan

{Bx − E(Bx − tan

{Bx})}] (2.8)

Eqn.(2.7) is usually used for force equations and eqn.(2.8) is used for Pneumatic trail as it is has a cosine shape. For y = longitudinal force (F

), x is the slip ratio (κ) and for y = lateral force (F

), x is the tyre slip angle (α).

The experimental analysis of the tyre force characteristics produced a sigmoidal­

shaped curve. In order to replicate it in the form of an equation, a sine ­ arc­tangent combination was used which eventually turned out to be the best fit (see fig.(2.1.7)).

The combination has a number of parameters viz. B, C, D and E. Each of these parameters has a particular influence on the shape of the Magic Formula curve.

• B is the stiffness factor. It stretches the curve to fit the linear range of the tyres.

• C is the shape factor which determines what part of the sine is used and thus influences the shape of the curve.

• D is the Peak factor that, as the name suggests, determines the peak value of the curve.

• E is the curvature factor that influences the characteristics around the peak.

The slope of the curve is a very important factor. The stiffness coefficients (C

and C

) discussed in the Simple tyre model are nothing but the slopes of the Magic Formula curve at origin.

D is usually tweaked to place the peak at the desired position. C controls the shape of the curve according to the relation [21]:

C = 2 − 2

π sin

( y

D )

B is the only unknown to be determined for the slope and is typically adjusted to meet the slope of the measured curve. Now, once the peak location (y

) and the peak input (x

) are determined, the curvature factor can be easily determined as follows:

E =

Bx

− tan( π 2C ) Bx

− tan

(Bx

)

Magic Formula has been revised several times to include turn slip, transient characteristics as well as combined slip conditions into a Semi­empirical tyre model.

The latest version PAC2002 [16, 21] has been studied and implemented in this thesis.

(Refer Appendix B for a detailed analysis and relevant equations.)

2.1.2.5 Identification of Systems from Step Response

Any physical system is a combination of various sub­systems ­ electrical, mechanical, electro­mechanical or various combinations of them. Different methods can be adopted to model each subsystem depending upon the accuracy desired, technical competence or available information at hand.

In the sections above, physical, empirical and semi­empirical methods were discussed.

Here, a different approach is discussed wherein the time­response of a dynamic system is analysed to generate a linear differential equation (esp. an ordinary differential equation (ODE)). Initially, the time response characteristics are converted into a single transfer function (in case of SISO systems) or multiple transfer function (in case of MIMO systems). Finally these transfer functions are converted into ODEs using Inverse Laplace Transform.

This method is usually adopted in control theory where the linear system transients are analysed for performance and stability. Note that the methods discussed below are referred from [24] and are intended particularly for Linear Time Invariant (LTI) Single­Input Single­Output (SISO) systems

First Order System (FOS)

A general first order system can be mathematically described by:

a dy(t)

dt + b y(t) = c u(t) (2.9)

where a,b,and c are the coefficients and u(t) is the input to the system. Thus, the first order transfer function is given by:

L{a dy(t)

dt + b y(t) } = L{c u(t)} (2.10)

= ⇒ (a s + b)Y (s) = c U(s) (2.11)

= ⇒ G(s) = Y (s)

U (s) = c

a s + b (2.12)

= ⇒ G(s) = c a b

b s + 1

= K

τ s + 1 (2.13)

where, K is the system gain and τ is the time constant.

Figure 2.1.8: Step Response [25]

Fig. 2.1.8 shows a typical FOS step response. Thus, in order to derive the transfer function eq(2.13), we need to determine 2 parameters ­ K and τ .

The system gain (K) can be determined according to the equation:

K = y

u

(2.14)

where, y

is the total change in the output signal at t = ∞ and u

is the magnitude of the input step signal.

Determination of the time constant is not very straight forward. There are many different graphical methods which try to determine as exact time constant as possible.

Here, the 63.2 % method has been discussed.

63.2% Method

In general, the time constant of a FOS (τ ) is the intersection of the line drawn from the meeting point of the tangent at the origin and the asymptotic value, to the point of inflection (point where the step response begins to change slope) (see fig 2.1.8). This was found to be the time coordinate of the point when the response reached 63.2 % of the final value, thus giving the method its name.Thus,

τ = t

(2.15)

Usually, in practice, the system consists of some non­linearities like time delay (L)

which are caused by the transport latencies. Thus, in this case, the FOS transfer

function looks as follows:

G(s) = K e

τ

s + 1 (2.16)

where, τ

is the effective time constant which is given by τ ­ L. This method is a good approximation of a linear FOS. But, in practice, a physical system is never a pure first order system. This means that either the system response does not have its steepest slope at the beginning or there are some oscillations in the response. Thus, in order to make a better first­order approximations of a higher order system, several other methods like The Tangent method, Sunderasan ­ Krishnaswamy method, etc. are often used.

Second Order System (SOS)

As we now know that most of the physical systems are higher order systems, it is a good option to investigate the second order transfer function of the step response. Typically, most of the higher order systems are approximated to either a first or a second order for ease of analysis and calculation.

Having said that, a general second order differential equation looks like:

d

y(t) dt

+ a

dy(t)

dt + a

y(t) = b u(t) (2.17) and the transfer function is:

G(s) = Y (s)

U (s) = b

s

+ a

s + a

(2.18)

In general, the second order transfer function is defined in terms of the damping ration (ζ), the natural frequency (ω

) and system gain (K) as follows:

G(s) = K ω

s

+ 2ζω

s + ω

(2.19)

One of the major advantages of representing the system in this form is because the roots

of the characteristic equation (denominator of eqn. 2.19) decide the overall behaviour

of the system. On these lines, the behaviour of the system can be categorised into 5

major parts:

Figure 2.1.9: Step response of an SOS with ω

= 2rad/s and varying damping ratio

• Critically Damped System

– If ζ = 1, the roots of the quadratic characteristic equation are s

= −ζω

i.e both the poles are real and equal. This generates a monotonic response without an oscillation and the output may reach the desired input level.

• Under Damped System

– If 0 < ζ < 1, the roots of the quadratic characteristic equation are s

=

−ζω

±iω

p 1 − ζ

i.e both the poles are imaginary with negative real parts.

This generates a damping oscillation in the output time­response.

– This response is characterised by various properties like Rise time (t

), Peak Overshoot (%M), Peak time (t

), and Settling time (t

). These characteristics help to determine the values of ζ and ω

or vice­versa.

• Over Damped System

– If ζ > 1, the roots of the quadratic characteristic equation are s

=

−ζω

± ω

p ζ

− 1 i.e both the poles are real and distinct. This generates a monotonic response without an oscillation.

– Over damped systems, as the name suggests, are very sluggish and may not reach the desired input levels.

• Un­Damped System

– If ζ = 0, the system oscillations are not damped and the time­response is sinusoidal in nature

• Unstable System

– If ζ < 0, the system oscillations increase in amplitude rendering the system unstable.

The identification of an under­damped system is trivial and is often taught in the Basic Course. Its derivation is out of this thesis’ scope and is thus not discussed.

Rather, the derivation of an over­damped system transfer function using Modified Harriott’s Method with Iterative Improvements is discussed in the following sub­

section [24].

Modified Harriot’s Method (for over­damped systems)

A general second­order transfer function with time delays can be written as:

G(s) = K · e

(τ

s + 1)(τ

s + 1) (2.20) where, the system gain (K) is calculated as in eqn. 2.14 and the main problem is to determine the time constants τ

and τ

. N.B.: If τ

= 0, the system becomes FOS and if τ

= τ

, it becomes critically damped.

In 1964, Harriot developed a graphic method to determine the transfer function (eqn.

2.20) from the step response. This method is typically based on a special parameter z = τ

τ

+ τ

. But, as it is easier to analyse w.r.t the actual time parameter t, it is transformed into the following procedure [24]:

• Determine the time t

at which the response reaches 72% of the total output change.

• Estimate the sum of time constants P

τ

using X τ

= 1

τ

(t

− L) ∀τ

= 1.25 (2.21)

• where L is the time delay which can be observed from the time response at hand.

If the delay is not significantly visible, assume it to be zero. If, at later stages, this approximation becomes invalid, some other value could be assigned.

• Calculate the time t

in the actual time scale where the step responses from

different types of systems deviate the most from each other by using:

t

= 0.4t

+ 0.6L (2.22)

• Note the output y

corresponding to the time t

.

• Normalize the output by calculating y

y

.

• Determine the parameter z = τ

τ

+ τ

corresponding to the y

y

value obtained in the previous step. This is basically read from a graph which is obtained for an over­damped system by varying the value of z (see fig. 2.1.10).

Figure 2.1.10: Normalized y

for different values of z [24]

– if y

y

< 0.27, L has to be increased.

– if y

y

> 0.4, L has to be decreased.

– when the value of L is modified, the calculation for t

and all the next steps have to be repeated.

• With the values calculated from the above steps, determine the time constants as follows:

τ

= z X

τ

, τ

= X

τ

− τ

(2.23)

The Modified Harriott’s method gives almost satisfactory output, but it is based on an

assumption that τ

is 1.25 but it slightly varies with z. In order to compensate this

error, an Iterative Improvement was proposed as follows.

Iterative Improvement

Figure 2.1.11: τ

for different values of z

• Read the corresponding value of τ

for the calculated z from the figure 2.1.11.

• If its same as 1.25, we are good to go. Else, calculate the new time constant sum as in equation 2.21.

• Calculate the new values of time constants from the equation 2.23 and obtain the desired SOS transfer function from eqn.2.20.

2.2 Methods and Methodology

The main aim of the Modelling part of this work was to derive a Dynamic model for the truck with high fidelity and better prediction capabilities as compared to the existing Kinematic bicycle model. This meant that apart from deriving the equations of motion, the vehicle model had to be matched against the physical truck to ensure similar behaviours at various velocities. Along the way, there were some choices made and methods adopted which are motivated briefly in the following section.

2.2.1 Model Choice

The first task was to choose a relevant vehicle model for which various factors like

complexity, scope, accuracy, application domain, computational cost, etc. had to

be considered. In order to make an informed choice, the thesis scope and general

requirements were revisited. According to the requirements, a lateral vehicle model

for state estimation and prediction had to be derived. This model would be used not only as a Virtual vehicle for Simulation purposes but also as a state predictor for Motion Planning and Lateral Control. Thus, the longitudinal dynamics were kept as simple as possible by neglecting the modelling of power­train. Furthermore, as lateral control is mostly related to yaw dynamics, the roll and pitch dynamics were neglected, and the model was narrowed down to a 6 Degree of Freedom (DOF) model.

At this juncture, a choice had to be made between a two­track model or a bicycle model.

[17, 18] showed that a bicycle model is a better choice for state prediction, and that the two­track model is computationally expensive and could be used as the simulation model if necessary. Thus, a Dynamic bicycle model was finalised (see section 2.1.2.3) with greater importance given to lateral dynamics.

Steering Dynamics

Steering Ratio

Dynamic Bicycle

Model

Controller

Planner

SWA SA

Vehicle states

ref SWA request

Figure 2.2.1: Cascaded Structure of the Vehicle Model (SWA: Steering Wheel Angle, SA: Steering Angle)

A Dynamic bicycle model may be typically of 2 types depending upon the tyre model used. The Linear Tyre model is computationally less expensive and performs well in the Linear range of side slip. On the other hand, the Nonlinear (Pacejka’s Magic formula) model is more accurate and works well even outside the linear range (see section 2.1.2.4). One downside is that as this formula includes a lot of parameters and calculations, it may be computationally heavy on a real­time system. Thus, both these models were pitched against each other to determine the best model for the target application.

Figure 2.2.1 shows a block diagram of the cascaded vehicle model. One aspect of the

vehicle that is usually neglected is its steering dynamics. The steering system is usually

assumed to be ideal without time delays and steady state errors which in reality is far from truth. Thus, a separate second­order model was derived using the Modified Harriot Method with Iterative Improvement (see section 2.1.2.5). This modelling method was adopted because the steering system was experimentally found out to be over­damped and the time response plot made it evident that the system is not a first order system as discussed in section 2.2.4.1. The output of this model was then converted into the Steering angle via the Steering ratio look­up table.

2.2.2 Tyre Models

Section 2.1.2.4 summarises the basic concepts behind the various tyre models. In the case of the Pacejka (PAC2002) Magic Formula, a tyre property file (.tir) was obtained from the Vehicle Dynamics Group at Scania AB and the equations were referred from [16]. On the other hand, the Linear tyre model was modified in order to mimic the experimentally obtained Sigmoidal curve.

The Piece­wise Linear Approximation of the sigmoidal curve is shown by dashed lines in the fig. 2.2.2 which was mathematically devised as follows:

F =

 



 

Cx, ∀ |x| < x

F

, ∀ |x| > x

(2.24)

where, C denotes the Stiffness coefficient and |x

| denotes the Linear slip limit for the truck. This value was experimentally determined to be ≈ 0.2 rad or 11.46°.

N.B.: Over the course of the Results and Discussion, the Dynamic bicycle model with

PAC2002 Nonlinear tyre model will be expressed as Dynamic Model (MF) and the

Dynamic Model with Piecewise Linear Approximation of tyre forces will be expressed

as Dynamic Model (Linear).

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

slip( )[ ] / slip angle ( )[rad]

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Force (Fx or Fy) [N]

10⁴

Tyre Force Characteristics under pure slip (MF)

Lateral Force Longitudinal Force

Figure 2.2.2: Tyre Force Characteristics under Pure Slip

2.2.3 Reference Pivot Point

The reference pivot point is the point on the vehicle about which the vehicle rotates.

The distance between this point and the Instantaneous Centre of Rotation (ICR) gives the Radius of Curvature (R).

The Scania Autonomous Framework considers the centre of the first rear axle as the Pivot point and thus the feedback and all other calculations are based on this reference.

But, as the COG is a better point of reference for investigating the effects of lateral acceleration and load transfer, it was decided to translate the net feedback velocity to the COG throughout this work. Even though this assumption gave rise to some systematic modelling errors, it did not create a quantitative bias in the comparisons as this assumption was kept uniform across all the vehicle models.

2.2.4 Test Scenarios

In order to ensure that the derived Dynamic Bicycle Model mimicked the physical

truck, some experiments were devised. In the process, logs were generated against

which the performances were compared. These experiments were classified into 2

different categories as follows:

7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12 time [sec]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

SWA [rad]

Step response for Amplitude pi/3

Expt Inp

(a) Step Response for step change of magnitude π/3

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

time [sec]

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

SWA [rad]

Step response for Amplitude pi/2

Expt Inp

(b) Step Response for step change of magnitude π/2

14 15 16 17 18 19

time [sec]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

SWA [rad]

Step response for Amplitude pi/4

Inp Expt

(c) Step Response for step change of magnitude π/4

Figure 2.2.3: Experiments and Observations for Identification of Steering Model

2.2.4.1 Identification of Steering Dynamics

As mentioned above, the steering dynamics of the truck had to be derived in order to account for the non­linearities in the steering system. For this, a non­parametric method of system identification was chosen as there was no reference parametric model available. Non­parametric methods involve exciting the system with external stimuli and then observe the system response to derive the model function. There are a number of standard inputs available to choose from ­ impulse (Dirac Delta function) input, step input, ramp input, parabolic input, etc. Each input corresponds to identification of a specific property or parameter of the model and thus, the stimulus had to be carefully chosen to cater to the needs of the modelling.

In this case, as the transient response of the system and its position errors are of specific

significance, the step input was selected. Thus, multiple experiments were designed

wherein a rectangular pulse with 71.43% duty cycle (5 sec On­time and 2 sec Off­time)

was given as the Steering Wheel Angle request to the steering controller. 3 different experiments were performed with step change amplitudes π/4, π/3, and π/2 and the steering responses were recorded in the form of a DDS log (refer figure 2.2.3).

From the above figure it was deduced that the steering system was an over­damped system as the responses never reached the input reference. Also, it was clear that the system is not a First Order System (FOS) as a typical FOS has the largest slope at the origin. Thus, a second­order linear approximation was performed in the form of a transfer function as follows:

G(s) = K · e

(T

s + 1)(T

s + 1) (2.25) Furthermore, sinusoidal inputs with varying frequency (0.1 rad/s to 10 rad/s) were used to determine the maximum input frequency after which no significant change in the output was observed. This value was later compared with the Bode plot analysis’

Bandwidth frequency and was found to be overlapping (Check section 2.3 for a detailed analysis.).

2.2.5 Model Matching and Parameter Tuning

After the entire system is modelled, the next and the most important step is Model Matching. This is done to ensure that the derived vehicle model has a behaviour similar to the physical truck. This meant that several experiments had to be designed such that the dynamic range of the truck was exploited and enough lateral acceleration was recorded. On this note, the scenarios shown in fig. 2.2.4 were tested.

The scenarios in sub­figures 2.2.4a and 2.2.4b were designed at a low constant velocity

of 5 m/s due to a velocity limit in the truck framework. This meant that the tyre forces

operated in the linear range itself. Thus, a third experiment was later devised wherein

the truck was manually driven at a velocity between 10 m/s and 15 m/s. In this scenario,

the lateral acceleration reached around 4.5 m/s

and the tyre forces went past the linear

range as required.