The information age has radically shifted every aspect of our lives, perhaps first and foremost our scientific ones.

Before the information age, 99% of humanity's interaction with technology was limited to a car, airplane, telephone, or television. While technology had had a dramatic impact on life-from medicine to cheap manufacturing-for the populace at large direct everyday interaction with technology was rare.

Internet has changed all that. It has brought high technology into everyday lives, to the point that it is unimaginable to do without. Science is no longer the exclusive domain of the scientist. Now we enter into the age of the makers.

Makers don't invent from whole cloth, they assemble. They take bits bits and pieces gathered from the far-flung corners of the world and combine them to do something new, something different. Oftentimes, something scientific.

Who is a scientist? Historically, science was largely performed by the nobles and well-todo. In the 19 th century the natural sciences-geology, astronomy, etc...-were regarded as social activities. Discovering a new plant, or planet for that matter, was as much a way to impress ones friends as it was to push back the frontiers of science.

It wasn't until the early-to-mid 20 th century that science become a career instead of a hobby. WWI and especially WWII were times of tremendous scientific progress, pushed by the various governments' fervent desires to discover and invent the secret weapons that would win the war. With the governments' desire came funding, stability, and jobs. Now in the early 21 st Century, science is undergoing another transformation. Scientific principles are no longer hidden in an ivory tower, papers are no longer locked behind closed doors. Individuals are now able to quickly and efficiently access the scientific world. In a way, science has gone full loop. In past times, only the rich had the time to indulge in idle fantasies. Now, we are perhaps all rich enough for such explorations. Now, we are becoming makers.

--------------Most hard science will still be done by professional scientists, just like most hard journalism is done by professional journalists. But just as the journalistic world is supplemented by the bloggers, so is the scientific world supplemented by the makers. We must not forget this burgeoning world of the amateur scientist.

To this end, our scientific output must not be hidden. It must be obvious, available, and accessible. In this spirit, this dissertation is an experiment, much the same that a PhD is an experiment. This dissertation is designed such that the results and application are easy to understand for those with a less-than-formal background.

It is difficult to express in words the journey of science.

It is an open-ended trip where the consequences of decisions are years in the making. An uphill trip littered with the refuse of cast-off ideas and rejected works. A sublime trip that is filled with the joy of accomplishment, of exploration, and of knowledge.

It is said that research is "spending six months of work in order to find out which two weeks made it all worthwhile". My thanks go to all those who helped me recognize my two weeks.

First and foremost, I'd like to thank my advisor, Professor Juergen Sachau for giving me an invaluable opportunity to work on challenging and extremely interesting projects over the past four years. He has always made himself available for help and advice and there has never been an occasion when I've knocked on his door and he hasn't given me time.

I would also like to thank my co-advisor, Professor Eric Busvelle. Without his extraordinary theoretical ideas and computational expertise, this thesis would have been a distant dream. Seldom have I encountered anyone with such clarity of vision.

Thanks are due to Professors Ugo Boscain and Jean-Regis Hadji-Minaglou for agreeing to serve on my thesis committee, and to Professors Thierry-Marie Guerra and Hassan Hammouri for taking the time to review the manuscript as my reviewers.

Special thanks to Professor Jean-Paul Gauthier, for whom nothing is impossible; and to Dr. Nicolas Boizot, who walked side by side with me on our common path, in Burgundy and Luxembourg alike.

iii Introduction Notation conventions

Introduction

Every age has its theme, and ours is energy. Our civilization is founded on the stuff. It warms our houses, keeps the lights on, runs the computers, and is present in every single aspect of our lives. It is no coincidence that, according to the Economist, three of the ten largest companies in the world are energy companies, and eight of the largest ten are involved either in transportation and/or energy.

Without energy, the author would not have written this paper, and the reader would not be reading it.

It is impossible to overestimate the impact energy has had on our lives, and perhaps only the most savage of cultures lost in the deepest jungles of the Amazon or the African heartland can even truly relate to what it means to go without.

Civilization as we know it would disappear overnight if we were to lose our energy source. Thus, it would seem wise to make every drop, so to speak, count. Unfortunately, progress along this front has been difficult. Some would have us believe that we can drill our way out of our problem. With the greatest respect [24], and in light of the ongoing Deep Horizons catastrophe in the Gulf of Mexico, the author disagrees.

The case has been laid out and it is clear that there is a pressing need to reduce energy waste and consumption.

--------------Everyone talks about saving gas, but do we actually know how to do it? Of course, we can always say that "better technology" will find a solution, but that begs the question "What exactly is better technology?" Do we need better technology in order to improve fuel efficiency, or can we do that now, with the existing tools we have at hand?

We might learn something about the current state of technology vs. the current state of applied controls by examining the Grand DARPA Challenge. The Grand DARPA Challenge was a competition to build an autonomously driven vehicle that would race across the Mojave Desert, on a 240km course driving at average speeds in excess of 35km/hr. The first year, 2004, the farthest any of the 21 teams' vehicles went was a little less than 12km. Contrast this to 2005, when the challenge was not only successfully met, it was completed by a total of 5 out of 23 teams. Indeed, only one team failed to best the previous year's record. While technology evolves quickly, it does not evolve that quickly.

No, we have to look elsewhere for an explanation of the dramatic improvement in only 365 days.

We have to look to controls.

--------------

The field of modern optimal control is a relatively young field of mathematics. Fiftyone years ago, Lev Semenovich Pontryagin published The Mathematical Theory of Optimal Processes, founding the modern study of optimal control.

(Of course, like all mathematics, even if optimal control is a relatively recent branch of mathematics, the roots go back much further, all the way to Bernoulli, who first mentioned the Principle of Optimality in connection with the Brachistochrone Problem1 , and Newton, who solved it.)

Pontryagin inspired many of his contemporaries. Indeed, Rudolph Kalman used his work to develop the Linear Quadratic Regulator in 1960, followed later on in the year by what is now known as the discrete-discrete Kalman Filter and then in 1961 the continous-continuous Kalman Filter.

Shortly before Pontryagin published his seminal paper, Bellman introduced the idea of Dynamic Programming, a sort of brute force approach to optimal control of discrete-time systems. Optimal control was found to express a fundamental property of motion in natural systems. Bellman demonstrated that the natural direction for solving optimal control problems is backwards in time.

Since then, much fundamental research has been performed in the field of optimal control. Clarke [START_REF]The maximum principle under minimal hypotheses[END_REF] and Sussmann [START_REF] Sussmann | A strong version of the maximum principle under weak hypotheses[END_REF] have relaxed and generalized its principles to overcome some of the original limitations of the PMP. Today, the advent of high-power computers in low-power packages is opening the door to optimizing processes that just a few years ago would have been too costly or time consuming.

--------------According to Benjamin Dessus, roughly half of the installed power base in the West is made up of cars and trucks. Obviously, there is a wealth of opportunities to use optimal control to make reductions in emissions and fuel consumption. However, all optimality studies depend on an engine efficiency model (or alternatively an engine specific fuel consumption (SFC) model). In previous studies, these models were assumed to be available ahead of time. [25,[START_REF] Van Keulen | Optimal energy management in hybrid electric trucks using route information[END_REF][START_REF] Rousseau | Véhicule Hybride et Commande Optimale[END_REF] Unfortunately, this is not the case in reality. Car manufacturers are notoriously tightlipped and guard their data jealously. Furthermore, engine data alone does not tell the whole story, as it is ideal data from an ideal engine under ideal conditions. As any driver knows, the manufacturer reported fuel economy is terribly imprecise in real-world situations.

For an optimality study, it is essential to know the true system efficiency, which cannot be predicted from engine data alone. Many other car-specific factors come into play, such as transmission losses, rolling resistance, air drag, etc... Even when these values can be found in literature or from the manufacturer, they are oftentimes unreliable (e.g. coefficients of air resistance) and have limited usefulness.

A better optimality approach would take into account the real-world parameters of the individual vehicle in question. This approach should be able to adapt to changing circumstances, such as mass change when a passenger gets out of a car or an air resistance change when the driver puts down a convertible top. This approach also would even allow monitoring motor efficiency with time, opening the door to early fault warning. (A tire low-pressure warning comes to mind.) --------------Can we improve fuel efficiency for every single car in the world, without modifying a single one of them? Can we do this without imposing burdensome-and dangerous-driving styles on the drivers? Can efficiency be improved without sacrificing air quality?

The author believes the answer is yes. But first we have to develop some tools along the way.

********************

Dissertation overview

Goal: The goal of this PhD is to develop a system by which we can optimize fuel efficiency in individual vehicles, based on data gathered from these same. We offer targeted, specific controls trajectories for each individual car and individual trip.

• Chapter 1 is a survey of the state of the art in modeling and optimal control.

• Chapter 2 details current state of the art in sensoring and processing components in addition to describing the sensor package and data logger developed during the research project.

• Chapter 3 shows the application to the fuel efficiency mapping problem. A case study is presented demonstrating the observer. There are also tips and techniques that are of interest to those who would like to apply this work in a real-world environment.

• Chapter 4 illustrates the process for applying Pontryagin's Maximum Principle to an arbitrary vehicle. Several case studies are presented, demonstrating the controller.

• Appendix A is devoted to the subject of non-linear observers, with a special focus on the continuous-discrete time domain and fusion of asynchronous data.

• Appendix B is devoted to the subject of optimal controllers, with a special focus on Pontryagin's Maximum Principle.

In addition to the ones cited above, additional appendices contain many of the tools that had to be created and/or invented in order to turn the theory into practice. Specifically: data protocols, schematics, a short description of the CAN bus and a working CAN logger, a description of the data files, an example of using Google Earth to visualize data, Quad coordinates, etc... Lastly, there is a worked example of Pontryagin's Maximum Principle, designed to be especially useful for engineers and scientists who have not yet been exposed to optimal control.

Notation conventions

Throughout this dissertation, the following notation conventions apply: Optimal control theory has only recently broken into the automotive industry. An important motivation for this progress is the increasing number of electronic sensors in modern cars, opening the door to more and better controls, for instance in the automatic control of the air-fuel mixture. Sophisticated control methods are currently used in a variety of applications, however, in the industry most such problems typically model very short time intervals and do not require optimal control techniques. Optimizing a passenger car's efficiency along its trajectory is still confined to the world of research.

The car efficiency optimization problem has been considered in [START_REF] Rousseau | Véhicule Hybride et Commande Optimale[END_REF]53]. Most of the literature has chosen to use Dynamic Programming (Sec. B.2.3) [START_REF] Bonnans | Parametric optimization of hybrid car engines[END_REF]25] or Non-Linear Programming (Sec. B.2.4), although recently some studies have examined the use of Pontryagin's Maximum Principle as it applies to hybrid cars. [START_REF] Rousseau | Véhicule Hybride et Commande Optimale[END_REF]53] The authors of [25] found that optimal control applied to semi-trucks could both increase fuel efficiency and decrease driving time, all while reducing wear and tear on the transmission.

This chapter is devoted to establishing the necessary models for (1) generating a efficiency map for an individual car and (2) creating the optimal control from this map.

State models

The optimization problem presented in this dissertation is split by necessity into two parts: observation of holistic system efficiency and application of optimal control to realworld examples. To this end, in the following two sections we present state-space models for 1) the observation problem and 2) the control problem.

The observer model represents the link between the sensors and the state, whereas the controller model needs only the state variables necessary for optimization. The upshot is that the 17-state observer model is reduced to a 2-state controller model, or even 1-state if the optimal target is only a one-dimensional speed target, and not the two dimensional pairing of position and speed.

Observer model

The observer model is designed to give an estimation of holistic fuel efficiency, based on data collected only from real-world driving. [START_REF] Sebesta | On mapping engine efficiency through normal driving[END_REF] It is built on a simplified model of car kinematics, ignoring complex effects such as increased rolling resistance during turns, etc... If additional measurements are available, the model can be extended in order to give an increase in accuracy, but with the associated increase in computation time.

Modeling

We begin by establishing the system model. Two reference frames are used:

Model variables

• An inertial reference frame ijk linked to the car, where + î points in the direction of forward movement, + ĵ points 90 degrees left, and + k points straight up from the car. See Fig. 1.1a

• A global reference frame N ED (North-East-Down) which is the local mercator projection of geodetic latitude-longitude-altitude coordinates. Note that in order to maintain right-handed orthogonality, Down represents negative altitude. Also, note that proper care must be taken while projecting geodetic coordinates into a mercator projection.

[52] provides a good explanation of the process, along with Matlab code. See

M a vi = ΣF = F thrust -F drag -F roll -F gravity (1.1)
F drag = 1 2 ρC d Av 2 i (1.
2)

F roll = C rr (M g cos(θ) + S k v 2 i ) (1.3) F gravity = M g sin(θ) (1.4) F drag F roll F thrust F gravity v i , a i ω, α θ Figure 1.2: Free-body car diagram where M a = J r 2 w + M
M a is the equivalent accelerated mass, and takes into account the angular momentum of the various spinning parts of the car, where J is the sum of all polar moment of inertias, and r w the wheel radius. This constant represents all the various bits and pieces that must accelerate in order for the car itself to accelerate.

Note that this model does not separate braking forces (both engine braking and friction brakes) from thrust forces. For non-hybrid vehicles, this is an unnecessary distinction, as only data gathered while the engine is propelling the car forward is interesting.

The car is assumed fixed to the road, and there is a no-slip condition in both longitudinal and latitudinal directions. Furthermore, it is assumed that the car undergoes no rolling rotation about the î-axis (φ, φ ≡ 0). Thus, v j , v k =0 and any acceleration Σa j or Σa k represents a centripetal acceleration, i.e. the car is on a curving path:

a j,b = S j v 2 i a k,b = S k v 2 i (1.5)
Note 1 S j and S k , the road straightness coefficients, represent the inverse of road curvature. Straightness coefficients are preferable to road curvature coefficients, since straightness coefficients are finite when on a straight section of road. Realistic values for S j seem to lie on the interval [-0.05, 0.05] and for S k in [-0.005, 0.005] Likewise, the rate of change of ψ and θ is dependent on the road straightness coefficients and forward velocity:

ψ = S j v i θ = S k v i (1.6)
Note 2 In the N ED coordinate system, +ψ represents a yaw to the right, and +θ represents a pitch upwards.

The N, E, D dynamic equations are written in spherical coordinates centered about the body reference frame:

Ṅ = cos(ψ) cos(θ)v i Ė = sin(ψ) cos(θ)v i Ḋ = -sin(θ)v i (1.7)
The engine turns at a fixed rate with respect to the wheels, neutral gear and differential withstanding. In a normal car transmission, (See 4.4.2), this ratio is different between gears but otherwise constant. ωe = vi N spd (1.8)

As we cannot directly measure N spd , we suppose that the engine and wheels are not connected, which is the same as taking N spd ≡ 0. While this is obviously physically incorrect, since we have very good independent sensors for estimating ω e and v i , this approach yields satisfactory results.

Lastly, a measured acceleration will always have slight biases in each direction, not only due to imperfections in fabrication and in mounting, but also due to changes of load distribution in the car. Note also that an accelerometer in a gravity field measures the pull of gravity in addition to the net body acceleration. These considerations lead us to the following equations:

a i,m = a i,b + M Ma g sin(θ) + b i a j,m = a j,b + b j a k,m = a k,b -g cos(θ) + b k (1.9)

Energy analysis

Instantaneous power-in is the volume of fuel consumed per second times the amount of energy per volume.

P in = k Vfuel (1.10)
Instantaneous power-out is the thrust force times velocity:

P out = F thrust v i (1.11)
Efficiency is defined as

η = P out P in (1.12)
Power-out is also equal to torque times angular velocity:

P out = T e ω e (1.13)
Thus, by rearranging eqns. (1.11) and (1.13) the torque is

T e = F thrust v i ω e (1.14) Note 3 Result 1.14 is completely independent of transmission gear ratios.      x i,m a i,m a j,m a k,m DEM ψ DEM alt GP S lat GP S lon GP S alt IN J duration IN J time                   (1.15)
where the state equations are where the observation equations are

ẋi = v i vi = 1 M a F thrust -sign(v i ) 1 2 ρC d Av 2 i + M C rr g cos(θ) + S k v 2 i -M g sin(θ) ψ = S j v i θ = S k v i Ṅ = v i cos(ψ) cos(θ) Ė = v i sin(ψ) cos(θ) Ḋ = -v i sin(θ) Ėin = P in θe = w e (1.
x i,m = x i a i,m = 1 M a F thrust -sign(v i ) 1 2 ρC d Av 2 i + M C rr g cos(θ) + S k v 2 i + b i a j,m = S j v 2 i + b j a k,m = -g cos(θ) -S k v 2 i + b k DEM ψ = ψ DEM alt = -D GP S lat = N GP S lon = E GP S alt = -D IN J duration = E in IN J time = θ e (1.
18) and where

• x i,m is the measured distance at time t k
• a i,m , a j,m , a k,m are the measured accelerations in the longitudinal, latitudinal, and vertical body-frame axes at time t k

• DEM ψ , DEM alt are road headings and road altitude, respectively, at time t k . This data is read from a DEM (Digital Elevation Map), when available. Note that it is entirely possible to get reliable heading data from the DEM without having the associated altitude data, and that furthermore these two measurements are decoupled in time.

• GP S lat , GP S lon , GP S alt are, respectively, the GPS latitude, longitude, and altitude measurements taken at time t k and geodetically transformed into the local N ED projection system.

• IN J duration is the duration the fuel injector is open, times the quantity of fuel per second of opening, times the amount of energy per quantity of fuel. In other words, ∆E in ∝ t, and so this is equivalent to the measured energy input at time t k

• IN J time is the time t k at which the injector opens, which represents a multiple of π radians (See Sec. 2.3.4.1 for more details.) of angular distance traveled by the engine since the previous measurement.

Note 4

The GPS is not used as a speed measurement, since it was impossible to know if our particular GPS unit measured speed directly through doppler shift or if it simply reported speed as a filtered change of position with respect to time. In the former case, it would indeed have been a good measurement. If it is possible to find a GPS that uses doppler shift, then it would improve system accuracy to integrate it into the measurement model.

Observability normal form

We present here a normal form of the observer, as described in [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF]. The system is identical to the previous model, but transformed such that the structure follows certain rules (such as each state in a block can only depend on the states above it, and at the very least the topmost state in a block must be an output). The observability normal form is useful for demonstrating observability (cf. Sec. A.1.2), amongst other things.

We would like to write the normal form system as

ż = Az + b(z) (1.19)
First we rewrite the variable names:

                                                                 x 1 = x i x 2 = v i x 3 = F T x 4 = ψ x 5 = θ x 6 = N x 7 = E x 8 = D x 9 = S j x 10 = S k x 11 = E in x 12 = P i n x 13 = θ e x 14 = w e x 15 = b i x 16 = b j x 17 = b k ẋ1 = x 2 ẋ2 = 1 Ma x 3 + sign(x 2 ) 1 2 ρC d Ax 2 2 + M C rr g cos(x 5 ) + x 10 x 2 2 -M g sin(x 5 )) ẋ3 = 0 ẋ4 = x 9 x 2 ẋ5 = x 10 x 2 ẋ6 = x 2 cos(x 4 ) cos(x 5 ) ẋ7 = x 2 sin(x 4 ) cos(x 5 ) ẋ8 = -x 2 sin x 5 ẋ9 = 0 ẋ10 = 0 ẋ11 = x 12 ẋ12 = 0 ẋ13 = x 14 ẋ14 = ẋ2 N spd ẋ15 = 0 ẋ16 = 0 ẋ17 = 0
(1.20) where the measurements are

y 1 = x 1 y 2 = x 4 y 3 = x 6 y 4 = x 7 y 5 = -x 8 y 6 = x 11 y 7 = x 13 y 8 = 1 Ma x 3 -sign(x 2 ) 1 2 ρC d Ax 2 2 + M C rr g cos(x 5 ) + x 10 x 2 2 + x 15 y 9 = x 9 x 2 2 + x 16 y 10 = -g cos(x 5 ) -x 10 x 2 2 + x 17 (1.21)
In order to transform system (1.16) into the observability normal form, we propose the following change of coordinates Φ(x) = z and its inverse Φ

-1 (z) = x                                                                                                    z 1 = x 1 z 2 = x 2 z 3 = x 3 -M g sin(x 5 ) -M C rr g cos(x 5 ) + x 10 x 2 2 z 4 = x 4 z 5 = x 9 x 2 z 6 = x 6 z 7 = x 7 z 8 = -x 8 z 9 = x 2 sin(x 5 ) z 10 = x 2 2 x 10 cos(x 5 ) + sin(x 5 ) Ma (z 3 -1 2 ρC d Az 2 2 ) z 11 = x 11 z 12 = x 12 z 13 = x 13 z 14 = x 14 z 15 = 1 Ma x 3 -1 2 ρC d Ax 2 2 -M C rr (g cos(x 5 ) +x 10 x 2 2 + x 15 z 16 = x 9 x 2 2 + x 16 z 17 = -g cos(x 5 ) -x 10 x 2 2 + x 17 x 1 = z 1 x 2 = z 2 x 3 = z 3 + M C rr g 1 -(z 9 /z 2 ) 2 + M a z 2 z 10 -z 9 z 3 -1 2 ρC d Az 2 2 M a z 2 1 -(z 9 /z 2 ) 2 +M g z 9 z 2 x 4 = z 4 x 5 = arcsin (z 9 /z 2 ) x 6 = z 6 x 7 = z 7 x 8 = -z 8 x 9 = z 5 /z 2 x 10 = M a z 2 z 10 -z 9 z 3 -1 2 ρC d Az 2 2 M a z 3 2 1 -(z 9 /z 2 ) 2 x 11 = z 11 x 12 = z 12 x 13 = z 13 x 14 = z 14 x 15 = z 15 -1 Ma z 3 -1 2 ρC d Az 2 2 +M g(z 9 /z 2 )) x 16 = z 16 -z 5 z 2 x 17 = z 17 + g 1 -(z 9 /z 2 ) 2 + M a z 2 z 10 -z 9 z 3 -1 2 ρC d Az 2 2 M a z 2 1 -(z 9 /z 2 ) 2 which gives ż =                 A 1 0 0 0 0 0 0 0 0 0 A 2 0 0 0 0 0 0 0 0 0 A 0 0 0 0 0 0 0 0 0 A 4 0 0 0 0 0 0 0 0 0 A 5 0 0 0 0 0 0 0 0 0 A 6 0 0 0 0 0 0 0 0 0 A 7 0 0 0 0 0 0 0 0 0 A 8 0 0 0 0 0 0 0 0 0 A 9 0 0 0 0 0 0 0 0 0 A 10      (1.22)
where

A 1 =   0 1 0 0 0 1 Ma 0 0 0   , A 5 =   0 1 0 0 0 1 0 0 0   A 2 = A 6 = A 7 = 0 1 0 0 A 3 = A 4 = A 8 = A 9 = A 10 = 0
In addition, we have

y t k =                
1 0| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1 0| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1 0 0| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1 0| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1 0| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1| 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0| 1

                z (1.23)
Eqns. (1.22) and (1.23) are the model's normal form. Note that

x 2 > v i,min > 0, that θ = x 5 ∈ [-π 4 ; π 4 ], and that cos(arcsin(α)) = √ 1 -α 2 (because |α| < π 2 )
. Therefore the function is well defined. Again, it bears repeating that the system is no longer observable when x 2 = 0, and certainly has numerical stability problems as x 2 becomes tiny. This will be taken into account in the following sections.

This change of variables brings the model into the form

ż =                 A 1 0 0 0 0 0 0 0 0 0 0 A 2 0 0 0 0 0 0 0 0 0 0 A 3 0 0 0 0 0 0 0 0 0 0 A 4 0 0 0 0 0 0 0 0 0 0 A 5 0 0 0 0 0 0 0 0 0 0 A 6 0 0 0 0 0 0 0 0 0 0 A 7 0 0 0 0 0 0 0 0 0 0 A 8 0 0 0 0 0 0 0 0 0 0 A 9 0 0 0 0 0 0 0 0 0 0 A 10 b 3 (z) 0 b 5 (z) b 6 (z) b 7 (z) 0 0 b 10 (z) 0 0 0 b 14 (z) b 15 (z) b 16 (z) b 17 (z)                                (1.24)
where

A 1 =   0 1 0 0 0 1 Ma 0 0 0   , A 5 =   0 1 0 0 0 1 0 0 0   A 2 = A 6 = A 7 = 0 1 0 0 A 3 = A 4 = A 8 = A 9 = A 10 = 0

Controller model

The controller model is very basic, as nothing more complicated is needed to model distance traveled and energy consumed. Thus, the only controls are the engine speed, the transmission gear ratio, and the desired torque. Calculating power input requires only the addition of the engine efficiency to the model. The rest of the variables in eq. (1.15) are not necessary for minimizing the energy usage.

In addition, we must define the cost, J, to be minimized across the controller trajectory. While for this work absolute energy consumption was chosen, this cost could easily include other factors, such as particulate output, C0 2 emissions, etc... (In which case, it goes without saying that the controller model would have to be remade in order to model these new costs.) A minimizing controller must have a cost, J, to minimize. A heat engine transforms the fuel energy into work energy and waste heat. In the case of a car engine, a part of the work energy is lost due to various frictions and thermodynamic inefficiencies in a car (air resistance, rolling resistance, transmission losses, etc...) a part of the energy is used to run the various car accessories, and the balance is transformed into kinetic and/or potential energy.

Modeling

We seek to minimize fuel consumption, that is to say, minimize E in

J = E in = T 0 P in dt (1.25)
T can be fixed-horizon or not, i.e. if we wish we can apply a time criteria to the optimization or not, depending on the goals of a particular driving phase. See Sec. 4.2 for further details.

Kinematics

We start with the basic net forces equation developed in eq. (1.1), given in its refined form:

vi = 1 M a F thrust -sign(v i ) 1 2 ρC d Av 2 i + M C rr g cos(θ) + S k v 2 i -M g sin(θ)
By combining eqns. (1.8) and (1.14) into the above, we arrive at

vi = 1 M a T e N spd -sign(v i ) 1 2 ρC d Av 2 i + M C rr g cos(θ) + S k v 2 i -M g sin(θ) (1.26)
We keep v i in the model, instead of transforming it into ω N spd because it is a continuous variable, whereas N spd is not.

Complete controller model

The state vector is

x = x i v i (1.27)
and the dynamic equations are given by:

ẋ =   v i 1 M a T e N spd -sign(v i ) 1 2 ρC d Av 2 i + M C rr g cos(θ) + S k v 2 i -M g sin(θ)   (1.28)
While the engine efficiency, η, is a function of many parameters, most of these are intrinsic to the engine construction. Most of the remainder are a function of the engine control system. Practically, the only two parameters that are left to the user to be controlled are the engine speed and torque. As luck would have it, these are the most important in determining engine efficiency. As such, we suppose η to be only a function of ω e and T e , and define a function φ and its polynomial approximation (to be used in Chap. 4).

φ(ω e , T e ) = 1 η(ω e , T e ) ≈ a 0 + a 1 ω e + a 2 T e + a 3 ω e T e + a 4 ω 2 e + a 5 T 2 e + . . . (1.29)
and thus we can calculate

P in = ω e T e η(ω e , T e )
= ω e T e φ(ω e , T e ) (1.30)

Discussion

The key parameter in the minimization problem is η. This is the parameter that links the controls with the cost. The model for η is difficult to estimate, as it depends on several other estimated parameters, one of which, F thrust is very noisy and itself difficult to measure/estimate. Thus the observation problem hinges on the best possible estimation of F thrust possible. This is in itself provokes a cascade of dependencies. A close examination of F thrust shows that it depends heavily on θ, and to a less extent several constants (C r r, C d , ρ, etc...) and the car's velocity (v i ). While v i is easy to measure and estimate (from GPS and odometry data), θ is somewhat difficult to measure in a dynamic system without a gyroscope. Of course, at a complete stop, θ can be inferred from the direction of the gravity vector, but when the car is accelerating it becomes difficult to have a clear measurement.

Other state parameters, such as N, E, are important only insofar as they contribute to improving the estimation of F thrust . The only exceptions to this are the engine speed and energy input states, (θ e , ω e , E in , P in ). The engine speed is important as it is a control. The input energy is important because it is the quantity to be minimized. However, for both of these the model and measurements are excellent, so their impact on the model is not as critical as F thrust . As frequently seen in this dissertation, controlling and observing have different needs, and thus onboard equipment can be split along these two lines.

Monitoring equipment is composed of a sensor package, a µC (microcontroller), and an SD card for data storage.

The controlling equipment is a µP (microprocessor) that calculates the optimal control trajectory for a given driving route. It bases its controls on information from previously observed sensor data and GPS routing software.

In an industrial setting, a car manufacturer would already have access to all this sensor data, and indeed in a more global sense to the efficiency map. However, the challenge is just that: to recreate the efficiency map with other techniques than depending on the manufacturer. This project's sensor suite was inspired by the MPGuino [17], a basic mpg (mile per gallon) calculator based on the arduino [1], that has inspired a large following due to its simplicity and ease of use. The MPGuino calculates fuel consumption based on a calibrated fuel flow and car speed. Unfortunately, the arduino has several structural defects which prevent it from being used for our needs, not least of all the choice of an basic model ATMega chip which has insufficient memory and peripherals.

There are several commercially available data loggers available, but these have the disadvantage of:

• being closed source

• being expensive

• having limited expansion

Microprocessors

A microprocessor puts most or all of a computer's CPU onto a single integrated circuit. The microprocessor has seen incredible progress since its introduction forty years ago. They see use in every single facet of life, from portable music players to computer peripherals to aerospace control systems.

This project uses a microprocessor to compute the optimal control.

Microcontrollers

Microprocessors are much like microcontrollers. In fact, a microcontroller can be seen as a microprocessor core with peripheral support, such as analog-to-digital converters, serial port protocols, etc... Today, a microcontroller has an impressive array of sensors and speaks a number of bus protocols.

Due to the addition of the peripherals, in general they do not perform at quite such dizzying speeds as pure microprocessors. However, they can throttle back to the point that they are only consuming nW of power.

While there are still speed gains to be made, and floating point math is to be avoided due to the lack of a floating processing unit, it is fully sufficient to use a microcontroller as the foundation of our sensor suite.

As such, an Atmel ATMega644p forms the core of the system. Chip speed is set by an 8Mhz internal RC oscillator. The microcontroller was programmed in C. The chip was chosen as it is relatively easy to develop, has a low price, requires no external components (thus saving development time), and has an active user community (www.avrfreaks.net).

Future versions will be built around the Atmel XMEGA, which employs the same instruction set as the ATMega series, but has much more advanced peripherals and a much faster (32MHz) and more stable (<.1%) RC oscillator. This allows true 32-bit timers, along with a much more accurate timestamp. 1Project source code can be found online at www.eissq.com [START_REF] Kenneth | [END_REF] 2

.3 Sensor package

The sensor package is made up of 4 sensors:

1. Three-axis accelerometer 2. GPS 3. Fuel injector timer (for measuring fuel flow and engine position)

Odometer

Note 5 While neither a gyroscope nor a compass is used in the sensor package due to the requirement for a low-cost design, the observer detailed below was nonetheless designed with the possibility of adding these sensors. If gyroscope or compass data is available, it can be directly used without modification.

GPS

GPS, first started by the American military in the 1970s, has in the last 10 years revolutionized geolocalization. Its weaknesses are largely due to a non-gaussian error (due in part to the fact that only satellites overhead can be seen) and the inability to penetrate buildings and earth, meaning that a GPS fix is oftentimes unavailable in urban canyons, and always unavailable in tunnels. GPS also suffers from a lack of (relative) precision, only being accurate to sub-meter scales.

GPS works on the principle of triangulation. The GPS receiver is based on the very simple idea that it can calculate absolute distance from signal propagation time, i.e. the difference in time between when a signal was sent from the satellite and when it was received by the receiver. With the signal propagation time the receiver can calculate distance to the the satellite's known position. Geometrically, the intersection of any three lines is sufficient to place the receiver in 3-dimensional space so the receiver seeks out three known positions in order to have a fix. (In theory. In reality, the GPS receiver must solve an equation with 4 unknowns, as an absolute time solution must also be found before the receiver can calculate a signal's ∆t. Thus, in practice a GPS receiver must have 4 satellites visible in the constellation before it can return a 3-D fix.)

Trivial as the idea may seem, the realities of compensating for time compression due to a satellite's relativistic speed with respect to the Earth, and time dialation due to proximity to Earth's gravity well are, to greatly understate, complicated.

The US is constantly upgrading the GPS network with new and better satellites, and many of GPS's inherent problems will be done away with by the promising new EU-sponsored Galileo system, but at the time of this dissertation Galileo was still many years away due to inter-EU squabbling. If it ever flies, it will herald in a new age of global positioning, as the accuracy and precision are far higher than GPS.

For our research, we used the San Jose Technology FV-M8 GPS. This particular model was very popular with hobbyists and researchers because of its 5Hz update rate and low price (<$50). While better GPSes have appeared on the market since, this one fulfilled needs.

However, if a new GPS unit were to be selected, we would take care to use one that output the raw GPS data, instead of only the filtered data. The filtered data poses problems as by its very nature the filtered data always lags behind the real-time state. If the raw GPS data were available, the position could be calculated as part of the observer model, allowing the Kalman Filter to do the work of filtering noise. This would help capture the true dynamics of the system, as manufacturers are loath to release the proprietary data of how they regulate their filters.

Care must be taken with GPS data to properly interpret and filter the measurement data. As show in Fig. 2.2, the GPS can have significant error, even after it has had complete 3-D fixes with > 4 satellites. Ostensibly, the GPS delivers an estimation of the tracking accuracy, the HDOP, but in our experience this value is to be taken skeptically, as we have already seen the GPS get "dynamically stuck" at a clearly wrong series of coordinates for extended periods of time, all while the GPS reports a very low HDOP.

Accelerometer

Current accelerometer technology is MEMS. From Wikipedia: "Conceptually, an accelerometer behaves as a damped mass on a spring. When the accelerometer experiences an acceleration, the mass is displaced to the point that the spring is able to accelerate the mass at the same rate as the casing. The displacement is then measured to give the acceleration."

Most accelerometers work in-plane, so they can only measure acceleration in one dimension. A 3-D accelerometer is made of 3 in-plane accelerometers mounted in one package, each pointing along a cartesian axis.

Consumer-grade accelerometers are of the piezoelectric type. When a force is applied to a piezoelectric crystal, a voltage is induced. By transducing this voltage the sensor can output an acceleration value. Currently, more and more accelerometers have direct digital outputs, sparing the designer the difficulty of creating a precision analog circuit.

The switch from analog to digital signals saves considerable time in the design and implementation of the circuit. Digital signals reduce-or completely eliminate-the need for intensive noise-rejecting circuit simulation and design. Indeed, at low data transmission speeds, the circuit pathway is not important. In fact, the only significant consideration is proper capacitance and grounding of the Vcc and GND lines, respectively. We used the STMicro LIS3LV02D three-axis digital accelerometer. We chose this accelerometer in particular because it has digital outputs, in addition to having class-leading resolution (12-bits).

The accelerometer uses the TWI bus (sometimes called by its trademarked name I2C) to communicate with the microcontroller at 500kbps, more than sufficient for the 40Hz refresh rate. Data is timestamped for future use.

Note 6

The accelerometer can also be used to measure road quality, which would give a better estimation of rolling resistance. [3]

Odometer

The total distance the car has traveled is essential to high-quality observability. The GPS does not have sufficient reliability-not accuracy, nor resolution, nor speed-for this. Odometry is the use of data from sensors to measure a change of position. All cars have some form of odometer or another. The oldest cars have mechanical odometers that are driven by a spinning cable attached to a transmission gear, but modern cars use electric signals, usually derived from the ABS system.

The odometer sensor in cars is usually a toothed wheel that rotates past a hall sensor. This contactless sensoring system is very robust, does not wear out, is unaffected by dirt and grime, and gives high-quality signals. Its largest drawback is that the sensor teeth must go past the hall sensor quickly enough to trigger it. This speed threshold tends to be around 3-5 kph.

Many, but not all, modern cars have a Vss (vehicle speed sensor) signal, which is a PWM signal sent over a dedicated line. This signal is calibrated to give x pulses per y kilometers, and is often used by a stereo's speed-sensitive volume control feature or a GPS navigation system's dead-reckoning module, although cars can also use it to drive the speedometer.

Note 7

The aspiring tinkerer who would like to install a data logger such as is described here is advised not to use the ABS signal directly, as it is a very high impedance signal and there is great risk of perturbing the ABS system, with unknown consequences to safety.

In this case, we take the 12V Vss PWM signal and use a voltage divider to bring signals down to acceptable levels for the µC.

While this is in general a low-noise signal, care must be taken. Data is as accurate as the wheel radius model. For simplicity, we assume that wheel radius is a constant, however this is a simplification, as not only does it decrease with increasing car weights (flattening of the wheel), but it also varies based on air pressure (function of atmospheric pressure, tire temperature, and centripetal acceleration).

In addition to these inaccuracies, jitter in the µC real-time clock can influence speed observations.

In one test case, the car's Vss was non-functional, so we were forced to develop an interface to the car's CAN bus. This is detailed in Appendix D.

Fuel injector timer

The µC measures the time that the fuel injector opens, and for how long it stays open. The result is that we have two measurements with one sensor: engine angular distance (explained in Sec. 2.3.4.1) and fuel volume injected (explained in Sec. 2.3.4.2).

We timestamp the rising and falling edges of the signal, adding an opto-isolator in order to ensure that the µC can never have an adverse effect on the vehicle's ECU. All negative voltage components are blocked by a diode, thus protecting the circuitry from reverse voltage spikes. In order to protect the car ECU from short circuits, a current limiting resistor is placed immediately at the voltage shunt.

Engine angular distance/velocity

Properly controlling the engine demands measuring its position with a great deal of accuracy. The engine position is measured through multiple sensors in order for the ECU to properly time fuel injection, spark, etc... Sensors such as the camshaft position sensor are robust and accurate, and can be excellent choices for measuring angular position.

However, the engine angular position can also be measured by sensing voltage pulses in the fuel injection or in the spark plug wires. This method lacks somewhat in resolution, as it can only give position in multiples of π (or thereabouts, as explained in the following).

The position can also be measured optically with a colored dot on a rotating element of the engine, for instance the main engine pulley. However, installing the optical reader is probably more trouble than it's worth since there are so many other, excellent sources.

While we cannot regularly sample the engine's angular displacement, we can know at each fuel injection impulse that the engine has returned to its previous orientation. How many revolutions that is depends on the car. Contrary to expectations, in a 4-cycle engine a single fuel injector can inject multiple times during one complete combustion cycle. • Grouped (Fig 2 .3b): In this scheme, injectors fire in groups. Depending on the number of cylinders, this could be any factor of the total number of cylinders. For instance, some Jaguars V-12 engines use four groups of three injectors, whereas some older Citroen engines have two groups of two. Each injector pulse indicates that the engine has rotated 4π # of groups radians.

• Sequential (Fig 2 .3c): Each fuel injector fires independently of the others. This is the most common in modern cars, but by no means ubiquitous. Each injector pulse indicates that the engine has rotated 4π radians.

Note 8 Some engines can switch between injector grouping styles, depending on engine conditions. For instance, sequential injection may be used for low engine speeds, switching to simultaneous injection at high speeds. The same principle is used in changing from light loads to heavy loads. Care must be taken to identify if an engine is capable of switching injection schemes. If so, additional data will be required in order to identify the injector scheme. This information can come from the observer or from the OBD-II interface, to name two possibilities.

Fuel flow meter

There are varying ways of measuring fluid flow. It can be directly measured by impeller blade, turbine, or even ultrasound. Except for the last one, all these methods require direct contact with the fluid being measured, which in the case of gasoline is difficult due to its reactivity, and the last one itself is prohibitively expensive. In any case, each of these requires installing extra equipment in the vehicle, which would understandably be unacceptable to many owners.

Having already encountered these problems, the online community explored the technique of using the fuel injector pulse width to measure fuel economy. This technique is accurate to ±1%. [17], and is thus well within the bounds of reasonable error. Note 9 All these techniques suffer from the fact that they measure volume, and not mass. This is a non-negligible difference, but unfortunately short of measuring temperature in the fuel tank and compensating based on a particular blend of gasoline's coefficient of thermal expansion, there is little that can be done about this. Fortunately, this can be disregarded as the relative efficiencies are of interest here. They will stay the same and so the impact will not be major.

This method infers fuel flow by measuring fuel injector pulse duration. A fuel injector is effectively a needle (pintle) that plugs a small hole (discharge orifice). When current flows through an electromagnetic coil, the needle is pulled back from the hole, allowing fuel to flow through. When current stops flowing, a spring pushes the needle snug with the hole, stopping fuel flow. (See Fig. 2.4)

While measuring the voltage pulse applied to the injector neglects some factors, such as opening and closing speed, opening and closing hysteresis, variations in fuel pressure, temperature, density, etc..., the reality is that simply measuring the beginning and end of the positive voltage purposes are quite consistent. There is certainly an unknown bias, but these can be measured and other errors are absorbed into the model and have no effect on the final control. Furthermore, this technique does not work with direct injection diesel engines, as the fuel pressure is no longer constant, violating the principle that fuel quantity injected is proportional to injector opening duration.

A possible solution to both problems could be to install an additional sensor installed in order to measure fuel flow in the fuel lines.

Another very strong possibility is to measure the oxygen sensor and MAF( mass airflow sensor) instead of the fuel injector timing pulse, as this method can compensate for changes in density, changes in fuel blend (95/98 octane, summer/winter, ethanol/gasoline, etc...) , and works as well with diesel fuel as gasoline.

Other sensors 2.4.1 Gyroscope

Inertial navigation systems [START_REF] Draper | Inertial Guidance[END_REF], as they pertain to car navigation, have been studied in countless research papers. [START_REF] Titterton | Strapdown Inertial Navigation Technology[END_REF] In most of these cases, a gyroscope is used, however we opted not to use one due to costs. Commercially available gyroscopes fall into three categories:

1. Mechanical gyroscopes-These are built around a spinning (gyrating) mass. They are large, bulky, require separate power supply, expensive, and prone to failure.

2. MEMS-Micro Electro-Mechanical Systems. These are not true gyroscopes. They are composed of two three-axis accelerometers mounted to a digital signal processor or microcontroller of some kind. These gyroscopes integrate angular acceleration to arrive at angular position. The reference vector is given by gravity.

This is generally effective, as average acceleration can be assumed to be zero in the long run, but has serious shortcomings when instantaneous acceleration is large. Furthermore, accelerometers tend to be noisy and have low resolution, so noise on the level of several milli-gravity is normal. Inertial orientation as resulting from the double integration is therefore particularly prone to error.

3. Ring laser-The nec plus ultra of gyroscopes, these are extremely expensive, and thus out of range a project that hopes to have wide-scale acceptance.

Compass

Accurate digital compasses are starting to see wider acceptance, finding their way into low-end consumer electronics. While in the end, they were not necessary for the desired level of accuracy, they are still a promising technology and can only serve to improve results, once the price comes down a little more.

Note 11

The author feels that it should be possible to geolocate a car with only a DEM (Digital Elevation Map), compass, and odometer. This exercise in pattern matching is left as a difficult, but challenging, exercise to the reader.

Additional equipment 2.5.1 Data storage

All data is stored on an SD card. MMC and SD cards have an SPI compatibility mode. (Note, the higher capacity (>2GB) SDHC cards do not.) This compatibility mode, well documented online, gives a basic access to the SD card over the SPI bus, albeit at reduced speeds.

In order to reduce computational complexity, some liberties are taken with the SD card's FAT16 file system. Specifically, a large file of several GB (2GB is the maximum size permitted in FAT16), initialized to zeros, is saved on the card. Then, data is logged sequentially in 512 byte blocks, starting at the beginning of the file. Because this file is not fragmented, the sequential writes always write within the memory-space allocated to the file.

Data retrieval is simple. The original file simply needs to be copied to the computer, and then the unwritten space pared off. Data accumulates at several megabytes a minute, so a 2GB card can be used for several months for daily commuter usage.

Further details can be found in Appendix C.1.4.

Digital Elevation Map

The DEM (Digital Elevation Map) is a mapping of points along the road. In this particular case, a shapefile formatted document was provided by the Luxembourgish mapping agency. This document included all roads in Luxembourg, with an altitude precision of 1m. Unfortunately, today DEM data is difficult to come by, to the point of being impractical. In most cases, this data is sold by mapping companies at exorbitant prices, precluding widescale usage. (The authors were fortunate in Luxembourg to have the support of the national mapping administration, and thus have free access to quality DEM data.)

While the NASA developed SRTM data [START_REF]SRTM[END_REF] is widely available, it is of limited resolution and, by its nature, misses nuances such as the difference between bridges and tunnels.

As free, high-quality mapping data [18] becomes more and more available, it can be used to better the model.

Figure 2.6 shows the advantages that the DEM can bring to correcting the GPS. Notice in specific the 20m discrepency at the 7000m mark. Clearly this sort of GPS error can wreak havoc with the slope estimation and as a consequence in the force estimation if not accounted for. Efficiency mapping A cornerstone of any vehicle optimal control study is an engine efficiency model (or alternatively an engine specific fuel consumption (SFC) model). Unfortunately, this information is very difficult to come by. In previous studies, these models were assumed to be available ahead of time. [25,[START_REF] Van Keulen | Optimal energy management in hybrid electric trucks using route information[END_REF][START_REF] Rousseau | Véhicule Hybride et Commande Optimale[END_REF] In reality, this is not the case. Car manufacturers are notoriously tight-lipped and guard their data jealously. Furthermore, engine data alone does not tell the whole story, as it is ideal data from an ideal engine under ideal conditions. As any driver knows, the manufacturer reported fuel economy are vary considerably in real-world situations.

For an optimality study, it is essential to know the true system efficiency, which cannot be predicted from engine data alone. Many other car-specific factors come into play, such as transmission losses, rolling resistance, air resistance, etc... The values, when they can be found, are oftentimes unreliable (e.g. coefficients of air resistance) and have limited usefulness.

A better optimality approach would take into account the real-world parameters of the individual vehicle in question. This approach should be able to adapt to changing circumstances, such as a mass change when a passenger gets out of a car or an air resistance change when a convertible top being is put down. This approach also would allow monitoring motor efficiency with time, opening the door to an early indication of fault development.

Model

We recall the model developed in Section 1.2.2, using the variables highlighted in Table 1 It is reasonable to suppose that certain measurements will not be available, depending on sensor state and the current environment. For instance, the GPS fix can be lost in a tunnel or in urban canyons. We present here an observability study for the various loss of sensors. These results summarize a theoretical study of the differentiable system observability property, but can also be inferred from intuitive considerations.

.1. ẋi = v i vi = 1 Ma F thrust -sign(v i ) 1 2 ρC d Av 2 i +M C rr g cos(θ) + S k v 2 i -M g sin(θ)) Ḟthrust = 0 ψ = S j v i θ = S k v i Ṅ = v i cos(ψ) cos(θ) Ė = v i sin(ψ) cos(θ) Ḋ = -v i sin(θ) Ṡj = 0 Ṡk = 0 Ėin = P in Ṗin = 0 θe = w e ẇe = 0 ḃi = 0 ḃj = 0 ḃk = 0 x i,m = x i a i,m = 1 Ma F thrust -sign(v i ) 1 2 ρC d Av 2 i +M C rr g cos(θ) + S k v 2 i + b i a j,j = S j v 2 i + b j a k,m = -g cos(θ) -S k v 2 i + b k DEM ψ = ψ DEM alt = -D θ m = θ GP S lat = N GP S lon = E GP S alt = -D IN J duration = E in IN J time = θ e (3.1)

Values used

• Loss of GPS: Clearly, without the GPS we lose observability for N, E, D1 , ψ. However, observability is preserved for the rest of the model.

• Loss of accelerometer: As the model is completely observable with only the GPS and the fuel injector sensor, this only impacts the observability of the accelerometer bias. By inference, we see that the model works with only a two-axis accelerometer available.

• Loss of Vss: No impact on observability.

• Loss of DEM: No impact on observability.

• Loss of engine speed sensor: (Even though one sensor measures both engine speed and fuel flow, for questions of observability we separate the two, as it is possible to have independent sensors for both.) No impact on observability only if the transmission ratio is known at all times.

• Loss of fuel flow sensor: Clearly, observability of energy flow is lost. Without energy flow, efficiency cannot be calculated, and so this is the most important sensor for our study.

Of course, it goes without saying that as we lose sensors we lose accuracy. For obvious reasons, we only consider losing one sensor at a time. The above observability study should not be used as a rationale for decreasing the number of sensors, but instead only for validating an approach which does not, for whatever reason, use all desired sensors.

Choice of observer

The choice of observer is paramount. The model is of a large dimension (17 states) and the final result, the efficiency η, involves dividing by a derivative of a measured state.

At first attempt, the UKF (Cf. A.3.2) was explored as a possible solution, but there were grave difficulties integrating the measurements into the model. In [START_REF] Mourllion | Kalman filters predictive steps comparison for vehicle localization[END_REF], Moullion provides a reasonable explanation of why this behavior is seen. He explains that the UKF's correction compares a physical mean, the measured value, with a sigma-space mean, the value calculated by the model non-linearity. Depending on this nonlinearity, there can be a strong deviation from the true mean, and it is difficult to transform the physical mean into the sigma-space. If the true value equals the estimated value, we would expect that the system would apply no correction. However, as can be seen in Fig. 3.1, the sigma mean value μ is not equal to the estimated value, and so there will be an error ǫ and in consequence a correction.

Sigma point mean

Specifically, the problem with the UKF came from the need to estimate θ as accurately as possible in order to likewise estimate F thrust as accurately as possible. (Cf. Sec. 1.4) Since the acceleration measurements are the most important measurements for estimating θ, this difference between the sigma space mean and the cartesian space mean provokes an estimation error that reverberates through the entire model and rapidly makes the results useless.

It is outside the scope of this dissertation to develop a generalized approach to converting measurements into the appropriate coordinate system in order to fully apply UKF. Thus in the end, in spite of the UKF's superior stability, its computational costs and sensitivity to the measurement model rules it out.

The classical EKF (Cf. A.3.1), although theoretically not as robust as the UKF, still gives satisfactory performance. As such, henceforth we base our observer on this observer. (In Sec. A.4, we include a discussion of the Adaptive High-gain Extended Kalman Filter [START_REF] Sebesta | Using an adaptive high-gain extended kalman filter with a car efficiency model[END_REF][START_REF] Boizot | An adaptive high-gain observer for nonlinear systems[END_REF] as a possible continuation of the EKF approach.)

We recall that this model is in continuous-discrete time, as described in Section A.2.

3 ẋ = f (x(t), u(t), t) y k = h(x k , u k , k)
We recall that this is a CD-EKF, as described in Section A.3.1.3, and using the shorthand described in (A.11) Prediction:

x-

k = xk-1 + t k t k-1 dx dt dt P- k = Pk-1 + t k t k-1 d P dt dt where dx dt = f (x, u) and d P dt = (A c P + PA T c + Q c ) P| t k-1 = Pk-1 (3.2)
and Q c is a positive-definite matrix in R nx×nx Correction:

K k = P- k H T k (H k P- k H T k + R k ) -1 xk = x- k + K k y t k -h(x - k , u) Pk = (I -K k H k ) P- k (3.3)
where R k is a positive-definite matrix in R ny×ny

Jacobian

The CD-EKF requires the Jacobian A c of the transition equations, and the Jacobian H k of the measurement equations.

The Jacobian A c can be numerically calculated, but finding the formal Jacobian for the system described by eqns. (3.1) is not complicated. Even if Jacobians seem complicated to some, they are straightforward to calculate, can be validated numerically and/or formally, and form the foundation of the EKF, one of the most widely-used observers in modern science.

A c is a (17 × 17) matrix with zeros everywhere, except at the following indices:

A c =                                0 A 1,2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A 2,2 A 2,3 0 A 2,5 0 0 0 0 A 2,10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A 4,2 0 0 0 0 0 0 A 4,9 0 0 0 0 0 0 0 0 0 A 5,2 0 0 0 0 0 0 0 A 4,10 0 0 0 0 0 0 0 0 A 6,2 0 A 6,4 A 6,5 0 0 0 0 0 0 0 0 0 0 0 0 0 A 7,2 0 A 7,4 A 7,5 0 0 0 0 0 0 0 0 0 0 0 0 0 A 8,2 0 0 A 8,5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A 11,12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A 13,14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                                (3.4) where A 1,2 = 1 A 2,2 = -sign(v i ) 1 Ma (ρC d Av i + 2C rr M S k v i ) A 2,3 = 1 Ma A 2,5 = M Ma g(C rr sin(θ) -cos(θ)) A 2,10 = -M Ma C rr v 2 i A 4,2 = S j A 4,9 = v i A 5,2 = S k A 5,10 = v i A 6,2 = cos(ψ) cos(θ) A 6,4 = -v i sin(ψ) cos(θ) A 6,5 = -v i cos(ψ) sin(θ) A 7,2 = sin(ψ) cos(θ) A 7,4 = v i cos(ψ) cos(θ) A 7,5 = -v i sin(ψ) sin(θ) A 8,2 = -sin(θ) A 8,5 = -v i cos(θ) A 11,12 = 1 A 13,14 = 1
Likewise for the output equations Jacobian, H:

H dist =
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

H acc =   0 H acc 1,2 H acc 1,3 0 H acc 1,5 0 0 0 0 H acc 1,10 0 0 0 0 1 0 0 0 H acc 2,2 0 0 0 0 0 0 H acc 2,9 0 0 0 0 0 0 1 0 0 H acc 3,2 0 0 H acc 3,5 0 0 0 0 H acc 3,10 0 0 0 0 0 0 1   H ψ = 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 H θ = 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 H GP S =
  0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0   H S j = 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 H IN J = 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 H DEM = 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 where

H acc 1,2 = -1 Ma sign(v i )(ρC d Av i + 2C rr M S k v i ) H acc 1,3 = 1 Ma H acc 1,5 = M Ma C rr g sin(θ) H acc 1,10 = -M Ma C rr v 2 i H acc 1,15 = 1 H acc 2,2 = 2S j v i H acc 2,9 = v 2 i H acc 2,16 = 1 H acc 3,2 = -2S k v i H acc 3,5 = g sin(θ) H acc 3,10 = -v 2 i H acc 3,17 = 1
See Section A.5 for an explanation of how to create H k . See Algorithm 1 for pseudo-code of how H k is used in practice.

Observer pseudo-code

The original code was written in Matlab, but is easily adaptable to any language.

← x(t k ) P- k ← P(t k ) Measurement y t k ← event(s) at time t k ∆t ← t k -t k-1 Correction construct new measurement Jacobian, H k construct new measurement covariance, R k K k ← P- k H T k (H k P- k H T k + R k ) -1 Pk ← (I -K k H k ) P- k xk ← x- k + K k (y t k -h(x - k , u) Update loop x(t k ) ← xk P(t k ) ← Pk k ← k + 1 end
We use a Runge-Kutta-4 integrator for integrating the continuous model and covariance differential equation systems. Note, however, that the interval itself is ∆t, and as such changes every iteration, not for numerical reasons as in classical adaptive step versions of RK, but because of the asynchronous measurements.

There is not much to gain by using a variable-step integrator, such as Runge-Kutta-4-5. Due to the nature of the observer loop, there tend to be multiple integration intervals between each update of a particular variable. For example, while awaiting a new acceleration measurement, the system will have been integrated multiple times across smaller intervals as other measurements, such as fuel burn, arrive.

Tuning

It can be quite complicated to tune a model with 17 states. In the following, several strategies are provided for tuning the observer.

Tuning state model covariances

• x i : the kinematic model is exact, the model noise is thus low.

• v i : the kinematic model is exact, the model noise is thus low.

• F thrust : the model is obviously incorrect-Ḟthrust is supposed equal to 0-and moreover is also subjected to very rapid swings. The model error is quite high.

Note 12 It is possible to improve the frequency response of Ḟthrust by using a local polynomial model, of the type Ḟthrust = ξ 1 , ξ1 = ξ 2 , ξ2 = ξ 3 , ... This type of model better represents variables that are likely to have rapid swings in state, as the modeled variable can vary like an n-order polynomial.

The cost of this improved response is the addition of dimensions to the state model, with a commensurate increase in processing costs. In this particular case, satisfactory results were obtained with the basic model.

• ψ: The car heading is dependent on a kinematic process, and in the long run the model is exact. (The car will always be on the road, always pointing in the direction the road goes.) The error is low.

• θ: The car pitch is dependent on a kinematic process, and in the long run the model is exact. (The car will always be sitting with all four wheels on the ground, always parallel to the road surface). The error is low.

• N, E, D: Northing, Easting, and Down are dependent on a simplified kinematic model (Coriolis effects are neglected), and thus the model is quite good. In the author's experience, a very low covariance noise value gives good results, as will be shown below.

• S j : the horizontal road straightness coefficient has a poor model, Ṡj = 0, but due to the very small nature of S j itself the tuning value will be low on an absolute scale. See Note 12 for suggestions on how the model can be improved.

• S k : the vertical road straightness coefficient has a poor model, Ṡk = 0, but due to the very small nature of S k itself the tuning will be low on an absolute scale. Moreover, we can be certain that S k changes slowly (otherwise driving a car would be like riding on a roller-coaster) so S k should be smaller than S j . See Note 12 for suggestions on how the model can be improved.

• E in : The model for E in is based on the conservation of energy, and is thus exact. The model noise is low.

• P in : The model for P in is poor: Ṗin = 0. Moreover, P in is capable of high frequency swings (at the speed at which the driver depresses the pedal) and thus model noise is quite high. See Note 12 for suggestions on how the model can be improved.

• θ e : The engine angular distance model is very good, indeed it is kinematically exact.

As a consequence the model noise is very low.

• ω e : The engine angular velocity should be considered to be a very noisy model, as 1) the engine angular distance measurement is very good so there is little risk to placing a high degree of confidence in the measurement and not in the model, and 2) when the transmission shifts or is put into neutral the engine changes speed very rapidly. See Note 12 for suggestions on how the model can be improved.

• b i , b j , b k :
The accelerometer biases should be considered to react glacially slowly, since the accelerometer bias is estimated consistently across multiple trips. There will be slight variations with different loading patterns, however. I.e. a car will tilt when passengers get in or out, thus provoking slightly different biases depending on passenger weight and distribution, so the biases must be allowed some room for movement.

Tuning measurement model covariances

• Accelerometer noise: in the k direction (upwards-downwards from the car), the accelerometer is far noisier than in the ĵ (side-to-side) and î (forward-backwards) directions. This is due to road noise from bumps, uneven surfaces, etc...

• GPS noise: GPS error is internally estimated by the GPS receiver, and expressed as HDOP (Horizontal Degree Of Precision), a value that is part of the data sent by the GPS unit. If using the NMEA protocol, this data is sent in the $GPGGA and $GPGSA packets.

Unfortunately, HDOP was not reliable in absolute terms, as the estimated error was frequently underestimated. After losing a GPS fix and then regaining it at a later time, the GPS position error could be orders of magnitude greater than that predicted by HDOP.

We suggest using a relatively high base value for the GPS noise covariance, and increasing it as a function of increasing HDOP.

• Vss noise: the Vss is a very clean signal, but the quantization of this signal can lead to errors. (See Fig. 3.2.) These quantization errors decrease with increasing time resolution, but unfortunately even at a time resolution of 125µs the error was higher than the author would have liked. Real-world results can be seen in Fig. 3.3. While this error can somewhat be smoothed out, it would be better to increase the sampling speed. Quantization errors aside, it seems that a 1% odometer error is reasonable.

• Engine angular distance noise: this has the same quantization error problem as the Vss, although there is for all intents and purposes no error whatsoever in angular distance traveled. (Obvious, since the injector only fires at the same point in the cycle every time.) This value can be very low and still yield good results.

• If a DEM is present, the noise at the precise node is very low, but determining if one is at a node, and if not how far away one is from said node, can be difficult. We had good experience with a value roughly 1000 times smaller than the GPS.

Values used

We present here the values used in the above-described observer.

Normalization

Another approach to observer tuning is to normalize all variables. States and measurements that are considered noisier, or less reliable, are given higher relative values, and viceversa for states and measurements that are more reliable.

Variable Value

Q(1, 1) 0.1 Q(2, 2) 0.01 Q(3, 3) 90000 Q(4, 4) 0.01 Q(5, 5) 1.5625 × 10 -8 Q(6, 6) 0.0001 Q(7, 7) 0.0001 Q(8, 8) 0.00001 Q(9, 9) 10 -6 Q(10, 10) 1/v 2 i × 10 -9 Q(11,
R dist 0.05 R acc (1, 1) 1 R acc (2, 2) 4 R acc (3, 3) 50 R ψ 0.03 R θ 0.03 R GP S (1, 1) (GPS HDOP) 2 × 100 R GP S (2, 2) (GPS HDOP) 2 × 100 R GP S (3, 3) (GPS HDOP) 2 × 1000 R Sj 10 6 R IN J (1, 1) 1 R IN J (2, 2) 0.016 R DEM 1 Table 3.3: Measurement covariance values
We used a non-normalized approach as it gives an absolute feel to the problem physics, but the reader might prefer the normalized approach as it is more abstract and focuses more on the computational math, and less on the application.

Data gathering

Data was gathered through normal car driving. Not requiring special actions on the part of drivers is an essential facet of the approach.

However, due to this constraint, a large part of the rpm vs. torque relationship is not explored. Specifically, low torque at high rpm is a relatively rare occurrence in normal driving. In fact, most normal drivers do not drive at high rpms (>4000rpm) no matter what the circumstances, and when they do it is only for brief periods of heavy acceleration, i.e. high torque. Nonetheless, these points are not strictly necessary, as the peak efficiency will not occur at these high rpm-low torque relationships.

Moreover, the optimal procedure will push the driver to explore high efficiency areas, at the expense of other zones, improving our model in the area that most needs it. In short, the most important area to characterize first is the zone of greatest efficiency.

In our experience, an hour of highway driving, along with several starts, was sufficient to fill a significant part of the curve, and have a reasonable curve fitting.

Results

Several cars were selected for sampling. The specific car represented here is a Smart Roadster Coupé, with a 60kW engine that has been chip remapped at some previous pointwhen, how and with what was unknown to the owner-to a higher output. This provides a good case study as even the manufacturer does not know the performance of this particular engine. The additional cars give similar results (Sec. 3.6), although a detailed study is not presented.

Each result considers three cases: 1) GPS and DEM data were provided to the observer, 2) only GPS data is provided, and 3) neither GPS nor DEM data is provided. Let us point out that in this third degraded mode observability is no longer preserved for variables N, E, D, ψ. However, observability for the rest of the system is preserved, and thus we can still calculate efficiency.

Note that efficiency is only plotted for moments when the car's engine is propelling the car. Remark that efficiency grows with increasing torque, which is to be expected. 2 Again, remark that efficiency grows with increasing torque, which is to be expected. Also remark that the system is noisier, with transient efficiencies reaching unreasonable values. This is certainly due to the fact that in-town driving is more chaotic, with rapid shifts from acceleration to deceleration phases.

Highway

In-town

Nonetheless, the results are consistent with the previous highway example.

Maximum acceleration

.6 shows results when the car is undergoing maximum acceleration. Efficiency is globally higher in maximum acceleration than the previous two examples. This is because there are less pressure losses behind the throttle body and thus the engine is not working as hard to intake air. However, efficiency here is already slightly suffering versus peak efficiency because at maximum power the car's engine management system injects additional fuel in order to help cool the cylinders. This additional fuel is not completely burned, and thus the engine has a lower fuel efficiency.

Efficiency map

An engine efficiency map, Fig. 3.7, is made by averaging together efficiency values as functions of rpm and torque. Note that the efficiency graph shows a continuous and smooth progression from low to high efficiency, and shows a sharp drop off to the right as the rpm . This is raw data from hours of driving, in urban and highway environments, and has not been post-processed in any way, except to average the data in each grid point, and that a grid point is not drawn if the standard deviation is too high (±5%). 3A sample efficiency map, Fig. 3.8, is given for reference.

Discussion

As can be seen, the data graphed in the three test cases coincide so well it is often difficult to distinguish between results with all sensors and results without GPS and/or DEM. Where these sensors make the biggest difference is in urban environments, with lots of sharp turns and acceleration/deceleration phases.

This shows that the efficiency model is very robust, in spite of the fact that the altitude estimation becomes unsatisfactory. Therefore, DEM and GPS data do not seem to be strictly necessary. This is most likely due to the fact that the absolute altitude data is not necessary. As the odometer permits estimating the car's acceleration, θ remains observable. More to the point, the altitude data is heavily filtered, as it is less reliable, and thus the estimation of F thrust depends most heavily on the acceleration measurement.

If needed, one can eliminate these two sensors, reducing hardware costs and simplifying the programming task. However, despite the robustness of the observer, it is still preferable to use GPS and DEM when possible. While not dramatic, the results are clearly worse without the altitude data, and so as a consequence the efficiency map will be less precise.

It is interesting to note in these plots that the theoretical property of observability is clearly confirmed in practice. That is to say, despite the noise the whole state necessary for calculating η can be retrieved from the measurements and in consequence the observer can be said to perform to satisfaction. Even in the more complicated urban cycle, error remains It must be noted that the efficiency values are globally too high. (Thermal efficiency is generally in the 25-30% range, and this is before taking drivetrain losses into account.) This is not a major concern, as it indicates a bias in one or another of the car parameters. For instance, if C rr , the coefficient of rolling resistance, is estimated too high then as a result F thrust will be overestimated, which in turn leads to a higher calculation of P out and finally in consequence a globally higher efficiency. Or, if the amount of fuel injected per cycle is underestimated, then there is a percentage bias that affects the entire curve. However, as long as any optimality study uses these same values, it will lead to a similar optimal control.

Curve fitting

Analytic optimal control approaches require an efficiency or specific fuel consumption model represented by an equation that depends on the system state. We used the mesh presented in Fig. 3.7 as the raw data for the curve fitting. The curve fitting is not very sensitive to the mesh granularity, and so a mesh resolution that tracks relatively closely to a 100x100 grid is a good compromise between data point averaging and resolution. This led to grid squares 50[rpm] wide by 2[Nm] tall. Averages were only included as explained in the previous section: there must be a minimum number of points making up the average, and the standard deviation must be sufficiently small (±5% efficiency).

It is clear that that the empty areas seen in the unexplored torque-rpm-efficiency relationship will have a tendency to allow potentially undesirable fitting effects. The lack of data at these points can lead to a fitting with impossibly high efficiencies, especially in the case of a polynomial fitting.

These effects can be handled by 1) forcing the graph to conform to expectations (a car engine will always have one or two efficiency islands, and in modern passenger cars designed to run efficiently at 2000-4000rpm it is reasonable to suppose that there will never be an additional efficiency island to the extreme right), 2) cutting the fitting off at any minimum and setting all η outside this minimum to 0, or 3) seeking more data by performing additional experiments to better map this area.

While #3 would seem to make the most sense, in fact it would be seen as unacceptable by the average driver, and so this approach should be overlooked as impractical. #1 is interesting because it is easy programatically. It suffices to line the righthand border with a column of grid points with zero efficiency. #2 is most likely the best, though, as it allows the curve fitting to most faithfully model the area of greatest interest.

Efficiency fitting

Here we seek to fit the curve η(ω e , T e ). Keeping in mind that the goal is to create an optimal control for minimizing fuel consumption, it is desirable to keep the model as simple as possible, by limiting the fitting polynomial order to the greatest practical extent. The resulting curve parameters are not very sensitive to the mesh granularity, nor to polynomial class.

Curves were fit by a classical least-squares algorithm, and the coefficient of determination, R 2 is given. In least-squares fittings, the problem is to find the coefficient vector c in

F(ω e , T e )c = η
where F is a matrix with the rows made by the polynomial function we are looking to fit (in the cases below, that can be η(ω e , T e ), φ(ω e , T e ) = 1 η , or P in (ω e , T e )) evaluated at the corresponding grid point. In order to ensure numerical stability, it is always good practice to first normalize and center the system:

F n (ω e , T e )c n = η n where m F = mean of columns of F (as a row vector) σ F = std. deviation of columns of F (as a row vector)

F n,(i,j) = F i,j -m F,i σ F,i
and m η = mean(η) σ η = std. deviation(η)

η n,i = η i -m η σ η
If the original coefficients are needed, they are given by

c = (F T F) -1 F T (F n c n σ η + m η )
where m η is a column vector where every element has the value m η The classes of polynomials presented below are chosen based on the highest degree of an individual variable, not the polynomial's degree itself. This approach was chosen because the results naturally grouped together in such a way. .10a shows a slightly or somewhat better residual than the 2nd order, but not significantly so. Figure 3.10b has a higher fitting coefficient, but the efficiency island turns into a ridge. While this might indeed be the case, it is difficult to tell from the real data, so with all other things being equal we prefer the simpler model. As seen in Fig. 3.11 there is a high degree of correlation between the polynomial model and the data. It would seem that Fig. 3.11a is good candidate for an efficiency map, as with only 6 parameters it seems to capture the essential efficiency island.

Conclusion

In the end, it would seem that simple polynomials with highest degree 2 are already sufficient versus higher order fittings. In the interests of simplicity, we will choose this fitting model. It has the added advantage of being very robust in having a well defined efficiency peak. However, if increased resolution is necessary, Fig. 3.11 shows promise.

Another possible fitting approach would be to estimate the polynomial coefficients directly using the EKF. However, this is not as reliable in this case, as one could miss special features in the efficiency graph. It appears preferable to build a non-parametric efficiency map, and then curve-fit it. Another reason for this approach is explained at the end of [START_REF] Busvelle | Observation and identification tools for nonlinear systems. application to a fluid catalytic cracker[END_REF].

It must be added that these curves cannot be completely used as published, as they do not show the engine's limitations. While these can only be mapped by a dynamometer test, for all intents and purposes they can be inferred by performing a series of WOT (wide open throttle) tests. Figure 3.7, for instance, is bounded quite nicely by the results from several WOT tests that have been averaged into the data. Likewise for Fig. 3.20. Figure 3.12 shows how the efficiency curve could look like with limitations on ω e and T e .

Specific fuel consumption fitting

Here we seek to identify the inverse, φ(ω e , T e ) = 1/η. This is a useful fitting because, as explained in Sec. 1.3.2, the controller model uses the inverse of η so it is natural to try to fit η -1 The residual for these second order polynomials is quite poor. This is certainly because we are trying to fit the inverse of a function that was quite well described by a polynomial with positive coefficients. It follows, then that this function will resist fitting by a low number of coefficients. Even with 4 th order variables, the polynomial shown in Fig. 3.15 still has a great deal of difficulty capturing the essential dynamics.

Polynomials of highest degree 3

Conclusion

In fact, the poor fit is not as drastic as it might seem. By examining 1/φ we can see how well the curve fitting corresponds to the original data shape. If the overall shape of the function is correct, then we can use φ in the optimal calculations and expect to see similar results.

P in fitting

Here we seek to identify the entire power input function, P in (ω e , T e ). This could be a very useful fitting since this describes the entire efficiency expression as it appears in the Hamiltonian (Cf. Ch. 4). The residual for these second order polynomials is quite good. We see that there is little to be gained by adding additional terms, and so a simple polynomial seems to fit well.

Polynomials of highest degree 3

Figure 3.17 does not show a significant improvement from the previous examples.

Polynomials of highest degree 4

These higher order polynomials seem to show interesting behavior, but it is difficult to tell which is the "true" model. In these cases, it is important to avoid over-fitting the data.

Conclusions

Higher order models do not seem to add anything intrinsic to the general form of the fitting, and so we conclude that lower order models are sufficient.

Validation

The Smart was tested on a dynamometer in order to validate the proposed method. The test was performed by putting the dynamometer into constant speed operation, and varying The constant speed steps were taken in units of 100rpm, over a range from 2500rpm to 6000rpm. The torque was varied between 25Nm and 80Nm.

As the reader can see in Fig. 3.19, the curve fitted results correspond nicely to the observed results in Fig. 3.9b. While there is a scaling effect versus the experimental data results, this is easily corrected by calculating some of the parameters rather than using the reported manufacturer values, e.g. C rr , C d , fuel injection calibration, etc... It is also possible that the dynamometer calibration be somewhat off, as a dynamometer equipment is usually used for showing a relative change (useful when hot-rodding a car) instead of absolute accuracy.

What is more important is that the maximum efficiency occurs at roughly the same point in the two graphs, and the overall form is identical.

Additional results

As evidence that the above experiment is in no way specific nor tuned to the Smart Roadster, we present the efficiency graphs for a 2005 Peugeot 206RC. The engine is of a completely different construction, has a different number of cylinders, different displacement, and most of all is naturally aspirated, whereas the Smart Roadster is turbocharged.

It is interesting to note that in these graphs the same basic patterns emerge again, as expected. There is a maximum efficiency island at medium torque and rpm settings, although in the Peugeot's naturally aspirated engine it seems to occur at a slightly lower rpm than in the Smart's turbocharged engine.

Note that results are repeatable enough that once the efficiency map is established, the engine model can be used to estimate other car parameters. For instance, it would be possible to estimate rolling resistance and thus give a sort of low-pressure tire warning. It is also possible to estimate mass, or to see the change that aerodynamic modifications make in the

Defining optimality

On an internal combustion engine's torque-speed map, the zone of maximum efficiency does not necessarily coincide with the zones of least emissions. [START_REF] Johnson | Hev control strategy for real-time optimization of fuel economy and emissions[END_REF] In certain cases, these zones can be far apart, and any optimal study destined for mass utilization must choose the optimal cost with care. MPG, for the non-native reader, refers to miles per gallon. 1In particular, there is a definite tradeoff between NOx emissions and fuel consumption. This is the consequence of the intersection of chemistry with thermodynamics. One of the most basic elements of thermodynamics explains that the higher the compression, the higher the Carnot efficiency. On the other hand, basic chemical analysis shows that the higher the compression, the higher the quantity of NOx formed as a result.

The challenge for a holistic control strategy is thus to simultaneously reduce emissions while increasing fuel economy.

Fixed time versus non-fixed time optimization

Optimal control approaches, such as PMP, can treat the optimal dynamic problem in fixed-time or non-fixed time. In the former case, the problem consists of minimizing the cost for a trajectory between 0 and T , T being fixed. In the second case, T is a parameter to be optimized (and thus the minimum cost will be less than or equal to the fixed-time case). The PMP conditions are different in the two cases, which justifies treating them separately.

From a practical point of view it also makes sense to differentiate between the two problems. If the goal is to get to a certain point (e.g. drive a certain distance, accelerate to a certain speed), the problem is clear when T is fixed. As long as a reasonable value for T is given, it is likely that the solution is also reasonable, even without taking into account the state constraints (maximum engine speed, road speed limit, etc.) If T is non-fixed, it is more difficult to imagine the solution ahead of time. For instance, minimizing the cost can lead to T being small (the cost being the integral of a positive value between 0 and T ), but inversely the smaller T is the more time spent at the upper limits of the motor's performance envelope, with the associated efficiency penalty, and at higher speeds, with the associated aerodynamic losses. (After all, the most fuel efficient way to get anywhere is not to use fuel: pushing the car, for example.) Clearly, it comes down to a sort of mathematic compromise.

In the case where T is non-fixed, the problem's meaning depends to a great deal on the goal (position or speed), as will be seen in the case studies, Sec. 4.4, 4.5.

Global vs. local optimization

It makes sense to think about practical optimization as a series of slices which can be joined together to produce a "more optimal" trajectory. While this trajectory is no longer properly called the optimal trajectory for an entire car trip, the practicality of actually applying the trajectory make the theoretical considerations somewhat less imperative. The practical aspects lead us to an approach where we try to establish a set of rules that drivers can follow naturally for different phases of the trip. Once the system has identified η as a function of T, ω e and optimized the trajectory, the objective is to transmit the information in some way that the driver can understand. ("Accelerate briskly rather than gently-but don't floor it either.")

A natural separation that springs to mind is to distinguish between acceleration phases, such as merging into traffic, and cruise phases, such as going long distances on highways and autoroutes.

Fuel economy pptimization in cities is a more arduous task, and has less impact in absolute terms, as it first requires a city traffic model, which is outside the scope of this PhD. 2Moreover, fuel optimization has less to say about in-city driving, as the number of kilometers driven in-city are far fewer than those driven in more rural environments. Note, however, that optimization that takes into account emissions, such as NOx, SOx, particles, etc..., is still very interesting in urban environments, and with the advent of hybrid cars these optimizations become much more practical.

Case studies

We will now study several cases with multiple goals and then calculate the corresponding optimal trajectory. This allows us to validate our approach, knowing that the real-world application requires a more exhaustive study.

The problem is that optimal controls are very restrictive: the slightest slip, the tinniest error and optimality is lost. Our goal, then, is to optimize driving styles in a human way, in a way that does not require that the driver have unearthly powers of prediction and calculation. With the right rules, the driver will be able to drive more economically while only having to look a short distance into the future.

Note 13

The "pulse and glide" technique is not treated here. This approach to high-efficiency driving involves a period of acceleration (the pulse) followed by a period in which the engine is shut off (the glide). While this approach is possibly more fuel efficient, it is too dangerous to be considered here, as in modern cars shutting off the engine removes the source of vacuum for the brake boost, leading to a situation where the brakes are no longer effective after one or two uses.

Fitting polynomial: For the case study, we have chosen to use the fitting φ = 1/η. This choice is motivated by the fact that it is an easily differentiable polynomial, which eases the analytic task. While the results presented in Sec. 3.4.2 show that this is not the best fit for the ensemble of points, it is nonetheless a good fit for the area of peak efficiency, and has the necessary behavior, i.e. a pronounced peak efficiency followed by clear dropoffs to the sides.

The fitting polynomial and coefficients are presented in Tab. 4

One dimension: acceleration on the straight and level

We suppose that the ground is perfectly flat, and that the car is only allowed to have a forward velocity. As such

   ẋi = v i vi = 1 M a T e N spd - 1 2 ρC d Av 2 i + M C rr g (4.1)
where ω e , T e , x i , v i are all ≥ 03 The target is v(T ) = v * , where t f = T is supposed non-fixed. We make no assumption about the distance traveled. The problem is thus to accelerate to the given speed while using the least amount of energy necessary. The problem of T fixed, interesting from a practical point of view as explained above, will be considered later.

We now apply use the PMP theorem (see Annex B.2.5). Combining the cost eq. (1.30) with eq. (B.13), we can write the Hamiltonian, H, as

H = λ ẋ + λ 0 J = λ x v i + λ v 1 M a T e N spd - 1 2 ρC d Av 2 i + M C rr g + λ 0 ω e T e φ(ω e , T e ) (4.2)
Since x i does not appear in the target, as here the problem is only in one dimension, the costate transversality conditions impose λ x (T ) = 0, and since λx = -∂H ∂x i = 0, λ x = 0 for all t.

In the maximal version of the PMP, λ 0 ≤ 0. However, when λ 0 = 0 we can normalize and take λ 0 = -1, as the costate vector's magnitude is unimportant. This leaves the case of λ 0 = 0. If λ 0 = 0 then in order for the Hamiltonian H to be constant vi must be 0, since the vector < λ 0 , λ > = 0, but a constant v i clearly cannot be an optimal trajectory. Thus we can conclude that λ 0 = 0 and thus we take λ 0 = -1. This yields the reduced equation

H = λ v 1 M a T e N spd - 1 2 ρC d Av 2 i + M C rr g
-ω e T e φ(ω e , T e ) (

From the Hamiltonian transversality conditions (Sec. B.2.5.1), we know that H(T ) = 0 and for an optimal control H(t) =constant.

Furthermore, from the fact that along the trajectory, vi ≥ 0, and ω e T e φ(ω e , T e ) ≥ 0, then for H(T ) = 0 to hold, the costate λ v ≥ 0. We can deduce that λ v = 0 is an absurd solution, as it requires that ω e T e φ(ω e , T e ) = 0, so therefore λ v > 0.

CVT case

A CVT (Continuously Variable Transmission) is a special type of transmission that is seeing increasing use, especially on hybrid and high efficiency vehicles. While the exact nature of the transmission depends on the particular manufacturer, the objective is to be able to have an arbitrary output speed for any given input speed. This allows the engine to always run as close as possible to peak efficiency for a given power demand.

The CVT model is particularly interesting as it appears that this is a major trend in modern cars. In the CVT case, we suppose that ω e and T e are the controls, u.

Examining the Hamiltonian more closely, and recalling that N spd = ω e v i is non-constant in a CVT equipped car, we can see that λv > 0 since λv = - which is a 3rd degree polynomial in (ω e , T e ) and only depends on v i and λ v The problem is reduced, thus, to a shooting problem. (Remember that φ = 1/η, so the numerical shooting problem can depend on φ or 1/η, depending on the approach.) As discussed in Sec. B.3, we can resolve the shooting problem with numerical methods.

∂H ∂v i = λ v 1 M a ω e T e v 2 i + (ρC d Av i ) > 0

T, non-fixed horizon

As can be seen in Fig. 4.3, the results are somewhat realistic-putting aside for the moment that the Smart Roadster has a 6-spd gearbox, and not a CVT. It takes around 10 seconds to reach 100km/hr, which is consistent with the car's real-world performance. 4 The car attains the maximum velocity of 180km/hr after about 1 min. Again, consistent with real-world performance. Nonetheless, due to the CVT assumption this demonstration is not very realistic.

Note that λ v (0) was determined by resolving: max (Remember that here the Hamiltonian by definition = 0, but also that along an optimal trajectory the control u must maximize H. In other words, any control which is not optimal will not only not maximize the Hamiltonian, but it will allow it to be other than 0.) Note that the Hamiltonian was not completely preserved, due to the numerical integration. There are two ways to solve this problem. The first is to use a numerical algorithm that preserves the Hamiltonian. [START_REF] Sanz-Serna | Symplectic integrators for hamiltonian problems: an overview[END_REF] This type of integration method is effective but difficult to put into practice. Another method, based on the theory of observers, was developed by Busvelle et al. and gives highly satisfactory results for a low numerical cost. [START_REF] Busvelle | Numerical integration of differential equations in presence of first integrals: Observer method[END_REF] This problem has not been considered as critical since the errors are within tolerance.

The trajectory in the control/efficiency space can be seen in Fig. 4.3d. Note that at t ≈ 10s the control reaches its limits and is no longer singular. Conceptually, this results makes sense, since initially the car is at a low speed, where efficiency losses are greater than drag losses. As the car accelerates, however, it becomes more important to accelerate harder, leaving the area of maximum efficiency, in order to minimize the time spent with the higher drag losses.

Multi-speed transmission case

The multi-speed case is interesting as it represents the vast majority of cars that have a limited number of fixed gears. In this case, our case study uses a six-speed gear-box. The Hamiltonian is written as before:

H = λ v 1 M a T e N spd - 1 2 ρC d Av 2 i + M C rr g
-ω e T e φ(ω e , T e ) but this time the control is made of the pairing (N spd , T e ) where N spd is the ratio ω e /v which in the case of a six-speed transmission can have six distinct values. 5 The two controls are not equivalent to a static feedback since ω e = N spd v

5 Not taking reverse into account, of course. Indeed, there is a theoretically interesting problem which is not examined in this PhD. Imagine a car on a perfectly sinusoidal road (in the vertical sense). The car starts off pointing "uphill" and the optimal target is to reach the top while minimizing fuel consumption. No one would doubt that the car's most efficient trajectory would never require rolling backwards in neutral to a "downhill" part before reengaging the transmission in order to benefit from some small forward velocity and the associated higher engine efficiency at moderate speeds.

However, this is only because there is a critical point at which the car's rolling resistance and aerodynamic drag losses become so large that to roll backwards is de facto less efficient. Obviously, if there were no drag Therefore, the Hamiltonian as a function of N spd is

H = λ v 1 M a T e N spd - 1 2 ρC d Av 2 i + M C rr g -N spd vT e φ(N spd v, T e )
and so using the chain rule λv = -

∂H ∂v i = λ v 1 M a ρC d Av i + N spd T e φ(N spd v, T e ) + N 2 spd vT e ∂ ∂ω e [φ(ω e , T e )]
While passing from the domain of the physical controls in the (ω e , T e ) space to the domain in (N spd , T e ), the constraints become dependent on the state. The maximum principle thus must be one which takes into account the state constraints. This is a more complicated theorem because the costate vector can have discontinuities that must be calculated. In this case, however, it is not clear that this will be necessary. In effect, the intuition is that the controls will gravitate around the zones of maximum efficiency, at the expense of the zones where the efficiency drops off rapidly, which occurs especially at the state constraints. Thus, as discussed in Sec. B.2.5.3, the optimal solution can be found without considering the state constraints, and then compared with the true constraints. If the solution does not violate the physical bounds in (ω e , T e ), then this is the optimal solution sought after. If this is not the case, however, then there is no other choice than to apply the Pontryagin Maximum Principle with state limits, such as described in [10] and in Pontryagin's book [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF].

T, non-fixed horizon

Figures 4.5 show the results from PMP with the same target as the previous CVT section. The car goes from 36km/hr to 100km/hr in 9.3s, which represents a medium-to-high acceleration of around 0.2g. The gears are changed at around 4500rpm, which is consistent with this kind of spirited acceleration. Note, however, that these trajectories correspond to a problem that is somewhat non-intuitive: reach a given speed while using the least amount of fuel possible, with time, T , on a non-fixed horizon. Since the cost is "paid" the instant the car reaches the final speed, this implies a relatively short total acceleration time.

In Fig. 4.5d, the control is superimposed on the control/efficiency graph. The control trajectories as a function of the gears can be clearly seen. At 120km/hr, the optimal solution is close to leaving the state constraint space and thus becomes no longer optimal. It is relatively easy to guess what the optimal trajectory will be after this point (by the Bellman dynamic programming principle as elaborated in Sec. B.3.3, but the PMP as it is applied cannot be used to calculate an optimal trajectory greater than 120km/hr. However, in most cases throughout Europe and the world the speed is limited to 120km/hr and so this is inconsequential.

losses whatsoever the optimal trajectory is to coast back and forth, only applying the control torque when the energy can be "bought" at the cheapest price, that is to say the point where vi gives the highest efficiency. Defining this critical point would be an interesting study, and might have real-world applications for systems that have very low drag losses and very steep efficiency curves. Trains come to mind, such as perhaps an electric train in a coal mine.

T, fixed horizon

This time, the goal is to reach the target speed v * in a time T f , while minimizing fuel consumption. To be specific, the goals are 140km/hr in 30s, which corresponds to a relatively normal acceleration while merging onto an autoroute.

The Hamiltonian is identical to the non-fixed horizon case, but this time the transversality criteria no longer give that H = 0 for all time t, although H constant along any optimal trajectory still (and always) holds. This is then a classical shooting problem, where λ v (0) must be determined such that v(T f ) = v * . Since this is a shooting problem in only one dimension, it is straightforward to solve, for instance by the secant or bisection methods.

Results are shown in Fig. 4.6. At first glance, they seem surprising.

In the beginning, the control clearly stays at maximum N spd , minimum T e , which could be seen as almost the same thing as putting the car into neutral. Indeed, if the model had had The analogy taken to the extreme would be the case of a driver who must drive at 140km/hr at noon tomorrow. Obviously, the most fuel efficient strategy will be to wait until the last possible moment, using not a drop of gas, and then accelerate with the optimal profile in non-fixed horizon.

Real-world tests

As with the earlier sections, the real-word tests were performed on a Smart Roadster Coupé with a 6-spd manual transmission. Since the fixed-horizon case gives an identical • The improvement was seen on a driving style that would normally have been considered to be quite efficient: slow starts, easy application of throttle, no sudden movements, etc...

• These results are in 1-D only, so it makes a certain amount of sense that there is such a stark difference between optimal control and "douce" driving. The "douce" driving style takes the car much further before reaching the target speed. So while the acceleration phase is not optimal, a soft driving approach will not be as bad in a 2-D case as this data might imply.

• It is surprisingly difficult to follow the optimal trajectory, if for no other reason then that the gas pedal is a very rudimentary interface. The optimal PMP controls require relatively constant torque, which is difficult-if not impossible-to manage without instrumentation.

Two dimensions: distance on the straight and level

In this section, we search the optimal trajectory for a target speed and distance, which perhaps reflects reality better as a more common desire is to go a certain distance then it is to simply reach a certain speed. We define the target N = {x i (T ) = x * i , v i (T ) = v * i } and suppose that x * i is big enough, on the order of several kilometers, in order to determine the long term optimal strategy. In this case, we can no longer ignore the relation ẋi = v in the now two-dimensional system (implying that the augmented system, as described in App. B.2.5, is dimension 4).

CVT

Returning to the CVT case-it allows for easier computations and/or formal calculations, and above all facilitates the theoretical work and taking into account the constraints-and recalling that N spd = ω e v i is non-constant in a CVT equipped car, we see that the Hamiltonian is now written as in eq. ( 4.2)

H = λ x v i + λ v 1 M a T e ω e v i - 1 2 ρC d Av 2 i + M C rr g -ω e T e φ(ω e , T e )
The adjoint vector system is thus

λx = - ∂H ∂x = 0 λv = - ∂H ∂v = λ v 1 M a ω e T e v 2 i + (ρC d Av i ) -λ x (4.4)
from which we see that λ x is constant. Suppose λ v (0) is fixed. Since the initial system conditions (but only system, not aug-mented) are now fixed, that is to say v(0) = v i,0 and x(0) = 0, we can determine λ x (0) with the following condition: max

N spd ,Te H(v i,0 , 0, λ v , λ x , N spd , T e ) = 0
We now have all the initial conditions for the augmented system (v i (0), x(0), λ v (0), λ x (0)) and can now integrate this system. If v i stays positive, we thus obtain a time T f for which x(T f ) = x * . For this time T f , we will also have calculated λ v (T ). Thus for a given λ v (0), we can calculate λ v (T ) at the final time T f

The transversality conditions implies that λ v (T f ) = 0. We thus come, again, at a shooting problem, where the objective is to find λ v (0) such that λ v (T ) = 0

Results are seen in Fig. 4.10. The controls stay at their lower limits for the entire trajectory, and thus the car starts off at v i (0) and advances under its own inertia, eventually slowing down to a steady state speed for this minimal power condition. We can explain this sensitivity by numerically calculating λ v (0). For this, we will calculate the speed that corresponds to the minimum power and show its optimality.

The minimal power P corresponds to the lowest possible ω e , T e (the lower-left point of the control space graph). We note that P = ω e T e η( ω e , T e ) (4.5)

By realizing that the force on the car is found by P = F v, we can solve for the corresponding steady-state speed v i :

1 2 ρC d Av 3 i + C rr M gv i = ω e T e (4.6) Let p = C rr M g 1 2 C d A q = - ω e T e 1 2 C d A (4.7)
Then by finding the cubic root of eq. 4.6, we obtain

v i = 3 -q + q 2 + 4 27 p 3 2 + 3 -q -q 2 + 4 27 p 3 2 (4.8)
With v i in hand, we can look for the initial costate vector values. Since H = 0 and v = 0, we get P = λ x v i which gives

λ x (0) = λ x = ω e T e v i η( ω e , T e ) (4.9)
To find λ v (0), recall that by using the costate equation and the transversality condition

λv = λ v 1 M ω e T e v 2 i + ρC d A v i -λ x λ v ( T ) = 0 (4.10)
T is the final time, which can be calculated from the target x * i since the speed is constant, and thus T = x * v i By integrating in reverse time, that is to say time-to-go, from s=0 to s = T we have

λ v (0) = λ x α 1 -e -α T (4.11) where α = 1 M ωe Te v 2 i + ρC d A v i
Eq. (4.11) shows the extreme sensitivity of the shooting problem when T is large (which will generally happen when x * i is large). Nonetheless, the control problem itself is numerically very stable, since the control is constant and on the limit of the control space. This calculation shows the optimality of the trajectory, at the condition of verifying afterwards that ( ω e , T e ) do indeed satisfy the conditions for maximizing the Hamiltonian. This can be done numerically, or even formally, since the considered polynomials are particular and of a known degree. It is also necessary to show that this extremal is optimal, which is clearly the case here, by uniqueness.

Chapter 5

Conclusion and Discussion

As cars become more and more advanced, we will have unheard of access to sensors. Thanks to these improvements, the principles of optimality as described in this paper are increasingly possible in a real-world scenario This dissertation has shown that it is possible to create a car efficiency model from an observer based on cheap and readily available sensors. This efficiency model can then be used to create an optimal control profile.

Impact and future work

Sensor suite

The sensor suite proved to be very reliable and accurate. Since the beginning of the dissertation work, the total cost manufacturing cost has dropped significantly, by almost 50%. It is currently being extended to measure data in all cars, especially hybrid and diesel cars. These future models will be able to directly read the CAN bus, reducing the amount of additional sensors necessary.

The third version is undergoing testing, and the sensor suite is already seeing usage in other laboratories. Aside from measuring data directly related to a car's functioning, there are many future applications for sensor suite/observer combination, such as more accurate road mapping, traffic monitoring, etc...

Observer

The observer as presented in this PhD provides a basis for estimating a car's efficiency with a high degree of confidence. While there are always additional model parameters that can be taken into account for an even higher precision, the observer already gives consistent, repeatable results. The observer results are in line with results taken from a test bench dynamometer.

The hope is that this observer can be used to give researchers options for studying car efficiency and related field, other than depending on car manufacturers. Furthermore, this observer gives results that are certainly better than the original manufacturer's data, as the results are a function of a specific car, with all the variation that can occur during a car's long life.

Likewise, individual engineers and makers will be able to pursue their projects without having to have a large organization backing them in order to convince automotive manufacturers to divulge information.

Lastly, the results are repeatable enough that once the efficiency map is established, the engine model can be used to estimate other car parameters. For instance, it would be possible to estimate rolling resistance and thus give a sort of low-pressure tire warning. It is also possible to estimate mass, or to see the change that aerodynamic modifications make in the drag coefficient. It could furthermore be used to track slow deteriorations in the motor's general health, serving as an early warning system for engine maintenance.

Future work is to continue developing the AEKF for estimating the system state, as it promises to solve some of the observation problems, e.g. that the noise frequency can actually be slower than the perturbation frequency and that the model is in constant perturbations since there is no control input. See App. A.4 and [START_REF] Sebesta | Using an adaptive high-gain extended kalman filter with a car efficiency model[END_REF].

Controller

As it is difficult to verify directly the fuel efficiency as predicted from the model in a non-lab environment, the optimal control provides a further validation of the efficiency model. Furthermore, while more cases must be studied for a complete application, the results presented here show that the optimal control approach is sound.

Nonetheless, while the PMP controller is easy to understand, it is difficult to implement onboard a car. Significant advances in ergonomics must be made before such a controller can be used in everyday scenarios.

Possible implementation strategies are:

• A gas pedal that when pushed to a certain position, the floor for instance, triggers the optimal control mode.

• An in-car meter that displays the costate vector

• A display that uses an up arrow or down error to indicate speeds which are too fast or too slow. The intensity of the arrow can indicate whether the setpoint error is small or large.

• A meter that displays engine efficiency, not fuel efficiency. Such a meter would take into account the kinetic and potential energy stored in a car.

It is hoped that any future studies explore first the implementation side, as it is becoming increasingly clear that the optimal control and the "average" control are radically different. The "average" control has developed naturally after a century of driving, so it is likely that the optimal control will in some way upset the traditional way of driving.

However, there are still possibilities to improve control while allowing drivers to continue driving as they always have. For instance, the car alternator can be switched on and off at appropriate times, as could the air conditioning. In hybrid cars, the battery SOC (state-ofcharge) can be managed intelligently as a function of road profile.

A.1 Observers and observability

Alone, measured data is of little use. It is often too noisy to be directly added to the model, and so some amount of filtering and processing is first necessary. The filter of choice is a model-based parametric filter, one that combines the system model predictions and the real-world measurements in order to estimate the internal state. In a deterministic case, these filters are known as observers.

An observer's purpose is to:

1. Reconstruct a system's state based on measurements.

2. Suppress system and measurement noise.

Reconciliate data

In many, if not most, cases, it is unfeasible to measure all system states. This can be for reasons of cost, space, impracticality, or simply impossibility. Thus we employ the observer to estimate these "hidden" states.

Generally, system models can be expressed in the following terms1 (although we will see in Section A.2 that they can be written slightly differently, depending on the time domain):

dx(t) = f (x, u)dt + dw(t) dy(t) = h(x, u)dt + dv(t)
where x is the state, y the measured output, u the control input, and dw, dv random variables representing process and measurement noise, respectively.

Note 14

The above model expresses the system in the most general sense. Traditionally, the noises dw(t), dv(t) are struck from the equation as they are by definition non-derivable. Furthermore, controls texts and publications oftentimes present the affine case, dx dt = A(x) + g(x)u or even the linear case, dx dt = Ax + Bu. In the dissertation, except when noted the most general case f (x, u) is always presented.

The challenge is to filter the noises dw and dv without filtering real-world dynamics, all while still responding rapidly to perturbations. Except in simple linear cases, we will always be forced to make some compromises, and this is before considering other issues, such as numerical stability, etc... However, some very good methods already exist and we will examine them in Section A. 3. This chapter presents an overview of observers and describes the three major time-domains used in sytem modeling, along with a discussion of techniques for handling data fusion, such as redundant and asynchronous data. Algorithms for two current filters, the Extended Kalman Filter and the Unscented Kalman Filter, are given, along with a discussion of the new AEKF.

A.1.1 Observer

In rigorous terms, an observer is a mathematical algorithm that reconstructs the entire state of a deterministic system based only on the system model and its outputs. One common family of observer is the linear asymptotic observer, a model of the form:

xk+1 = Ax k + Bu k + L(y k -h(x k )) (A.1)
where the matrix L is such that x approaches x as k → ∞ for all initial conditions x(0) and x(0). This type of matrix is called a Luenberger observer.

The Kalman Filter can be seen as a special variant of a Luenberger observer for an observable time-invariant system at T ∞ . However, the Kalman method is not the only way to create a Luenberger observer. Indeed, any matrix L that places the eigenvalues of (A-LC) in such a way that their real components are all negative is of this family.

A.1.2 Observability

Observability describes whether or not there is sufficient information to observe all the system states. For instance, imagine that we wanted to know a car's speed between point A and B. As speed is just distance over time, if we have a distance measurement (odometer) and a time measurement (clock), we can calculate speed. This simple example system is observable.

Now, imagine we define a new system, one with states speed, distance, and path. This time we want to know the path the car took. If we still only have two instruments, the odometer and the clock, clearly we do not have enough information to observe the car's path. We would need additional information, such as a GPS or overhead photos of the car along its route. Thus, with only an odometer and clock this extended system is not observable.

Definition: Roughly speaking, a system is observable if, by knowing the controls leading up to a state and knowing that outputs of that state, it is possible to backtrack in time in order to find an earlier state, or apply the dynamic equations in order to predict where the system will be in the future.

In more rigorous terms, a system is said to be observable if for every vector x d in R n there is a finite time T such that if x(0) = x d and one applies an arbitrary control u, the knowledge of the trajectory on y during the time interval [0, T ] allows one to determine x d .

This definition of observability is the basic definition that corresponds to our intuition. There are, however, many more nuanced definitions possible, such as for local, infinitesimal, etc..., that are sometimes for useful. For supplemental definitions, see [START_REF] Busvelle | Observation and identification tools for nonlinear systems. application to a fluid catalytic cracker[END_REF].

Linear observability: In the linear case, the observability test is straightforward. For a linear system of form:

   dx dt = Ax + Bu y(t) = Cx + Du the system is observable if the Kalman Criterion Matrix Q o =      C CA . . . CA n-1      (A.2)
is of full rank n x . [START_REF] Novel | Commande linéaire des systèmes dynamiques[END_REF] Non-linear observability: In the nonlinear case, it can happen that a system is observable at certain states and non-observable at others. (For instance, in the system described by eqns (1.16), (1.17), (1.18), when v i = 0 then the system becomes non-observable.) It is difficult to talk about global observability. Instead, we are interested in locally observable systems. First, however, we must make a quick detour to define distinguishability.

Definition: distinguishability-Two states x 0 (0) and x 1 (0) are distinguishable if there exists a bounded input u such that y 0 (t) = y 1 (t) for some t ∈ [0, T ] Definition: local observability-A non-linear system is locally observable at x 0 if there exists a neighborhood of x 0 such that every x in that neighborhood other than x 0 is distinguishable from x 0 .

A.1.2.1 Observability test by Lie derivative

One test for local observability depends on the Lie derivative. (The derivation and purpose of Lie derivatives is far outside the scope of this PhD. If the reader would like further information, s/he is encouraged to read Boothby's excellent description. [START_REF] William | An Introduction to Differentiable Manifolds and Riemannian Geometry[END_REF])

The Lie derivative of h with respect to f is

L f h = ∇hf = ∂h ∂x f
By convention, L 0 f (h) = h. Higher-order Lie derivatives are give by:

L n f (h) = ∂ ∂x L n-1 f (h) f Thus, let l(x, u) ≡    L 0 f (h) . . . L nx-1 f (h)    (A.3)
where f, h are the dynamic and measurement models, as defined in eq. (A.1).

Note 15

If ẋ = f (x, u), then ḣ = ∂h ∂x ∂x ∂t = ∂h ∂x f = L 1 f (h). Similarly, ḧ = L 2 f (h), ... h = L 3 f (h)
, and so on and so forth. Thus the vector l(x, u) is nothing more than a series of derivatives of the output model.

Finally the test for local observability is that

Q(x 0 , u) ≡ ∂l(x 0 , u) ∂x (A.4) is rank n x A.1.2.

Observability test by normal form

A second test for observability is whether the model can be put into the normal form of observability. Indeed, all methods by which the model is brought into the normal form are by definition methods that show observability (no matter the change of coordinates, observability is an intrinsic property). [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF] Therefore, we can verify observability and compute the change of coordinates in one single action, as below.

Note 16 Both the Lie Derivative and the Observability Normal Form approaches are valid, but in Sec. 1.2.3 we present the normal form of eq. (1.15), as it has further uses beyond demonstrating observability.

The normal structure is unique in the single output case, and thus in this special case we can call it the canonical form, but not in the more general multiple output case the normal structure is not unique. We use the most common generalization, but there has been no study as of yet to create a sort of "zoology" of normal forms. [START_REF] Boizot | Adaptive High-gain Extended Kalman Filter and Applications[END_REF]

A.2 Time domains

Presented here is a short discussion of the three time domains into which systems and their observations can fall, following whether or not the dynamics and/or the outputs are discretized.

A.2.1 Continuous-continuous

These models are of the form

ẋ = f (x(t), u(t), t) y = h(x(t), u(t), t) (A.5)
These models assume that the system is continuous in time, and that the measurements are also continuous. In consequence, the prediction of the state and state covariance matrices are continuous.

Note 17 These cases are frequent in mathematical analyses, but are unlikely to be analyzed by an observer in real-world situations. Observers are almost uniformly implemented on digital computers, and so the following two cases that have a discrete component to the outputs are of more interest in real-world controls problems.

A.2.2 Discrete-discrete

These models are of the form

x k = f (x k-1 , u k-1 , k) y k = h(x k , u k , k) (A.6)
Observers based on this model rely on a discrete prediction of the state and the state covariance matrix, and a discrete correction of these two. This is the most commonly used, as it is the easiest and cheapest to program.

For linear systems, there is little reason to stay in continuous time when designing an observer. Time-invariant linear systems can be converted to discrete-time systems by [START_REF] Novel | Commande linéaire des systèmes dynamiques[END_REF]

ẋ = Ax + Bu ⇒ x(t) = e (t-t 0 )A x 0 + t t 0 e (t-τ )A Bu(τ )dτ (A.7)
which, with the condition that u k is constant along the interval [t k , t k+1 ], yields

x k+1 = A d x k + B d u k y k = h(x k , u k , k) (A.8)
Furthermore, if A is invertible:

A d = e AT , B d = A -1 (e AT -I )B (A.9)
In the more general case, the continuous time-domain system (A.5) can be discretized by

x k+1 = Φ(x k ) y k = h(x k , u k , k)
where Φ is the resolvent2 of the ODE in (A.5) between t k and t k+1 and where

A d = ∂Φ(x k ) ∂x k
which can be calculated by using the Jacobian of A(x), denoted A c = ∂A(x) ∂x

Ȧd = A c A d A d (t k ) = I (=Identity)
However, the most common and expedient approach to discretizing continuous models is to first linearize the model about its current state and then discretize as in eq. (A.7). This implies taking the Jacobian, either numerically or analytically. For finite steps in time, this is inherently less accurate than using a continuous model for prediction.

A.2.3 Continuous-discrete

These models are of the form ẋ = f (x(t), u(t), t)

y k = h(x k , u k , k) (A.10)
The continuous-discrete variant is a hybrid, best-of-both-worlds solution for physical problems treated on digital computers. It models a dynamic system that changes in continuous time, but has measurements arriving at discrete intervals.

In fact, the measurements do not even need to be synchronized. The ease with which asynchronous measurements is handled is a great advantage of the continuous-discrete system over the discrete-discrete one.

For the vast majority of physical systems, continuous-discrete time is the logical choice, since they fall neatly into the two time domains. Such systems change continuously with time t, but the data measurements are at discrete times t k . Nonetheless, discrete-discrete can be interesting in terms of numerical computation time.

Note 18 Continuous-discrete systems pose an interesting theoretical problem for the definition of observability. The attentive reader will notice that there is no guarantee that there is ever a time when all measurements arrive simultaneously. This leads to a good question: "Is this system observable?" In this specific case, the answer is clearly "yes". However, to the best of the authors' knowledge there is no general study of observability for asynchronous measurements, as in [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF] for the continuous case and [2] for the adaption to the discrete case. The problem is clearly not obvious since one can build an academic system for a continuoustime model which is observable for continuous measurements and for synchronous discrete measurements, but not for asynchronous discrete measurements.

3 3 Example: Note Π = {(a, b) ∈ R 2 ; a 2 + b 2 = 1} the unit circle. Let X = Π × R * + × R * + and the system Σ on X.        ẋ1 = -x2x4 ẋ2 = x1x4 ẋ3 = 0 ẋ4 = 0 ,    y1 = x1x3 y2 = x2x3 y3 = x4
and x3, x4 > 0 1. The system Σ is observable, whether it is continuous or discrete. Specifcally,

y 2 1 + y 2 2 = x 2 3 (x 2 1 + x 2 2 ) = x 2 3
give x3 and thus y1 and y2 suffice to give x1 and x2

2. The asynchronous system is non-observable because there are certain values for x4 which makes x3 indistinguishable.

Suppose that y1 and y2 are measured alternatively at a step time ∆t. The solution to Σ is

       x1 = cos(x4t + C1) x2 = sin(x4t + C1) x3 = C2 x4 = C3
If ∆t is poorly chosen such that x4∆t = π 2 + 2kπ then

A.3 Kalman Filters

Kalman Filters refer to a wide class of model-based stochastic filters. While originally developed as a discrete-time linear filter [START_REF] Emil | A new approach to linear filtering and prediction problems[END_REF], the Kalman Filter was quickly extended (no pun intended) to continuous time and non-linear systems [START_REF] Jazwinski | Stochastic Process and Filtering Theory[END_REF]. In the 50 years since Kalman's seminal paper, the term "Kalman Filter" has grown to encompass a much larger series of filters.

It goes without saying that different filters will have different advantages and disadvantages. For instance, the Extended Kalman Filter has good noise rejection properties, but can be slow to respond to perturbations, and global convergence has never been shown. Indeed, to the contrary, it has been shown that it does not globally converge.

The High-gain Kalman Filter [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF] demonstrates global convergence, but it overly sensitive to noise. The Adaptive High-gain Extended Kalman Filter combines the EKFs noise rejection properties with the HGKF's responsiveness, but adds a number of additional parameters to be tuned. [START_REF] Boizot | Adaptive High-gain Extended Kalman Filter and Applications[END_REF] The Unscented Kalman Filter presents interesting possibilities for better predictions of non-linear systems, but adds additional tuning parameters that are poorly understood, greatly increases computation time, and does not cope well with a non-linear measurement model.

We present in the following two Kalman Filter variants that are seeing wide application, the EKF and the UKF. These filters were chosen as likely candidates for modeling η (Chap. 1). We also present a new Kalman filter that presents an interesting direction for further study. [START_REF] Sebesta | Using an adaptive high-gain extended kalman filter with a car efficiency model[END_REF]

A.3.1 Extended Kalman Filter

The Extended Kalman Filter (EKF) is a tried-and-true approach used extensively in science and engineering. While there are other filters available, such as the Unscented Kalman Filter [START_REF] Julier | A new extension of the kalman filter to nonlinear systems[END_REF], the DD1 Filter [START_REF] Nørgaard | Advances in derivative-free state estimation for nonlinear systems[END_REF], the High Gain Kalman Filter [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF], etc..., the EKF has proven to be a very suitable approach for a wide class of problems, including vehicle localization [START_REF] Mourllion | Kalman filters predictive steps comparison for vehicle localization[END_REF]. Some of EKF's advantages are simplicity, computational costs, tuneability, and years of experience in real-world applications.

The EKF requires calculating the Jacobians of various model functions. 4 We use the

             y1(0) = x3 cos(C1) y2(∆t) = x3 sin(C1 + π 2 ) = x3 cos(C1) y1(2∆t) = x3 cos(C1 + 2π 2 ) = -x3 cos(C1) y2(3∆t) = x3 sin(C1 + 3π
2 ) = -x3 cos(C1) . . . and so it becomes impossible to find both x3 and C1, with the obvious conclusion that the system is non-observable.

following shorthand to represent commonly used Jacobians:

A c = ∂A ∂x i B c = ∂A ∂u i H c = ∂H ∂x i (A.11) A.3.1.1 Continuous-continuous EKF
In the CC-EKF, both the model and the measurements are continuous in time. This is the most classic form of the EKF, and yet the hardest to use in a strictly correct sense.

5 ẋ = A(x, u) + K(t) (y(t) -h(x, u)) Ṗ = A c P + PA T c + Q c -K(t)R c K T (t) (A.12)
where

K(t) = PH c R -1 c (A.13)
and where Q c is a positive semi-definite matrix in R nx×nx and R c is a positive-definite matrix of dimension R ny×ny

A.3.1.2 Discrete-discrete EKF

The DD-EKF has a discrete model and discrete measurements. Common engineering and scientific usage is to discretize linear and non-linear systems as this form is particularly cheap to calculate.

Prediction:

x-

k = f (x k-1 , u k-1 , t k-1 ) P- k = A k Pk-1 A T k + Q d (A.14)
where Q d is a positive semi-definite matrix in R nx×nx Correction:

K k = P- k H T k (H k P- k H T k + R k ) -1 xk = x- k + K k y k -h(x - k , u k ) Pk = (I -K k H k ) P- k (A.15)
where R k is a positive-definite matrix of dimension R ny×ny

A.3.1.3 Continuous-discrete EKF

The CD-EKF is used for systems with continuous dynamics, but discrete measurements. Prediction:

x-

k = xk-1 + t k t k-1 dx dt dt P- k = Pk-1 + t k t k-1 d P dt dt (A.16)
where dx dt = f (x, u) and

d P dt = (A c P + PA T c + Q c ) P| t k-1 = Pk-1 (A.17)
and Q c is a positive semi-definite matrix in R nx×nx Correction:

K k = P- k H T k (H k P- k H T k + R k ) -1 xk = x- k + K k y t k -h(x - k , u k ) Pk = (I -K k H k ) P- k (A.18)
where R k is a positive-definite matrix of dimension R ny×ny

The continuous-discrete Kalman Filter has found applications in diverse industries. One of the advantages of the CD-EKF is that it allows for continuous integration until a new measurement arrives, making it, when seen under a certain light, the fastest filter, i.e. it has the shortest computational time between measurement input and calculation of the state estimate. Thus in spite of additional integration costs vs. the discrete-discrete filter, the continuous-discrete can be the best solution for a real-time observation problem.

A.3.2 Unscented Kalman Filter

The Unscented Kalman Filter is a recent innovation [START_REF] Julier | A new extension of the kalman filter to nonlinear systems[END_REF] that is based on the assumption that it is easier to approximate a gaussian distribution than it is to approximate a non-linear function.

The UKF works by creating an n-dimensional cloud of sigma points around the expected mean. These points are precisely determined in order to exactly approximate the first, second, third, and fourth moments (kurtosis). The exact number of sigma points required for such an approximation is 2n + 1, where n is the number of states in the dynamic system. [START_REF] Julier | The spherical simplex unscented transformation[END_REF] The properties of the UKF have been studied in detail elsewhere [START_REF] Julier | A general method for approximating nonlinear transformations of probability distribution[END_REF][START_REF] Julier | A new approach for the nonlinear transformation of means and covariances in linear filters[END_REF], but we sum up the primary advantages here:

• As the mean and covariance of x are reproduced to the second order, the calculated values of the mean and covariance of y are correct to the second order as well. This implies that the mean is calculated to a higher order of accuracy than the EKF, whereas the covariance is calculated to the same order of accuracy. [START_REF] Julier | A new extension of the kalman filter to nonlinear systems[END_REF] • UKF requires no calculation of derivatives or Jacobians, operations that are generally non-trivial, human error-prone, and can lead to significant implementation difficulties.

The UKF was originally designed as a discrete-discrete filter. In additional to the covariance parameters-similar to the EKF case-there are three additional parameters that must be chosen: α, β, κ Note: if the true distribution of x is Gaussian, β = 2 is optimal.

Prediction:

X k-1 = [x k-1 . . . xk-1 . . . xk-1 ] + √ c 0 Pk-1 -Pk-1 X k = f d (X k-1 , k -1) x- k = X k w m P- k = X k WX T k + Q k-1 (A.19)
where

c = α 2 (n + κ) w m = [W (0) m . . . W (2n) m ] T W = (I -[w m . . . w m ]) × diag(W (0) c . . . W (2n) c ) × (I -[w m . . . w m ]) T
and where

W (0) m = λ/(n + λ) W (0) c = λ/(n + λ) + (1 -α 2 + β) W (i) m = λ/(2(n + λ)), i = 1, . . . , 2n W (i) c = λ/(2(n + λ)), i = 1, . . . , 2n
Correction:

X - k = [x - k . . . x- k . . . x- k ] + √ c 0 P- k -P- k Y - k = h(X - k , k) μk = Y - k w m S k = Y - k W[Y - k ] T + R k C k = X - k W[Y - k ] T K k = C k S -1 k xk = x- k + K k (y k -μk ) Pk = P- k -K k S k K T k (A.20)
However, in [START_REF] Sarkka | On unscented kalman filtering for state estimation of continuous-time nonlinear systems[END_REF], the UKF was extended to the continuous-discrete domain, and called the Unscented Kalman-Bucy Filter. In [START_REF] Van Der Merwe | The square-root unscented kalman filter for state and parameter estimation[END_REF], a more efficient square root form was given, which we detail below.

Y - k = h(X - k , k) μk = Y - m w m q, S - k = qr decomposition W (c) 1 (Y i=1:2L -μk ) R k S k = cholesky update{S - k , Y i=0 -μk , W (c) 0 } C k = X - k W[Y - k ] T K k = (C k (S T k ) -1 )S -1 k xk = x - k + K k (y k -μk ) U = K k S k A k = cholesky update{A - k , U, -1} (A.32)
The algorithm, which is not particularly obvious in the literature, is as follows:

Algorithm 2: SR-UKF-CD while t < T do form X t k-1 from m t k-1 and A t k-1 while no new measurement do integrate d X(t) dt end X- k ← X(t k ) Correction construct new measurement covariance, R k extract A - k from X- k S - k ← the square, lower-diagonal matrix S - k from QR decomposition S k ← cholesky factor of S - k + sign W (c) 0 W (c) 0 (Y i=0 -μk )(Y i=0 -μk ) T U ← K k S k for i = 1 to columns(U) do A - k ← cholesky factor of (A - k -U i U T i ) end Update loop m(t k ) ← mk Â(t k ) ← Â- k k ← k + 1 end Note 19
Taking the cholesky factor at each iteration instead of the cholesky update is not efficient. It is written in the above way only for algorithmic clarity. An efficiency SR-UKBF should use a cholseky update. Correction: 5. G k ≥ 1 is the high-gain parameter, explained in Sec. A.4.4. In addition, G 0 = 1 6. z 0 and P 0 are initial guesses for state and state covariance, and furthermore P 0 is symmetric, definite, positive matrix.

K k = P- k C T k (C k P- k C T k + R G (k)) -1 ẑk = ẑ- k + K k (y t k -C(ẑ - k )) Pk = (I -K k C k ) P- k I d,k = k i=m+1 ∆t i ||γ(y i - ∼ y i )|| 2 R ny G k = F(G k-1 , I d,k ) (A.35) where 1. Q G , R G (k)

A.4.3 AEKF Theory

Since, unlike the single-output case, there is no unique normal form for multiple output models [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF] the observer must be designed accordingly. This representation is described immediately below.

Definition 1 As earlier, z ∈ R nx denotes the state vector. z is divided into sub-vectors denoted z i for i ∈ {1, . . . , p}, where each one exists in a compact subset of R n x,i . Hence the dimension of z i is n x,i with p i=1 n x,i = n x . The dynamics of a given sub-vector z i = (z i,1 . . . z i,n x,i ) T are of the form:

żi =         0 α i,2 0 . . . 0 α i,3 . . . . . . . . . . . . . . . 0 α i,n x,i 0 . . . 0         z i +        b i,1 (z i,1 ) b i,2 (z i,1 , z i,2 ) . . . b i,n x,i -1 (z i,1 , . . . , z i,n x,i -1 ) b i,n x,i (z)       
While the other b i,j (•) do not, the sub-vector field components of the form b i,n x,i (•) do have the right to depend on the full state.

The full system is of the form eq. (A.33). In matrix C k , the row correponding to the i th output is composed of zeros; except for index (n x,1 + n x,2 + ... + n x,i-1 + 1), which is equal to α i,1

The hypotheses are:

1. Although omitted, α i,j can be time dependent. All the α i,j (t) functions are bounded and cannot vanish (which would lead to a loss of observability).

The vector field b(z) is supposed compactly supported and Lipschitz.

For such systems we can design an AEKF as in Sec. A.4.2 which has the following convergence property:

Theorem 1 For any time T * > 0 and any ǫ * > 0, there exist

1. 0 < d < T * that define an innovation 2. G 0 > 1 and G max > G 0 3. a function G k = F(G k-1 , I d ) such that G k ∈ [1, G max ]
such that for any time t ≥ T * :

z(t) -ẑ(t) 2 R nx ≤ ǫ * e -a(t-T * )
where a > 0 is a constant independent from ǫ *

A.4.4 AEKF considerations

Understanding the AEKF requires introducing certain concepts. We start with the highgain matrices Q G and R G . For reasons of Lipschitz constants (explained in [START_REF] Boizot | Adaptive High-gain Extended Kalman Filter and Applications[END_REF], Chapter 5):

Q G = G∆ G Q c ∆ G (A.36)
where Q c is a symmetric positive semi-definite matrix, and ∆ G is a diagonal matrix with elements determined by 1. The normal form leads to a block-wise diagonal matrix A, as in eq. (1.22). The matrix A is split into p square submatrices. Define n * x as the width of the largest submatrix.

2. In ∆ G , for each square block ∆ G,i (of width n x,i ) corresponding to one of the p blocks in A (also of width n x,i ), the elements are determined by

for k=1...n x,i ∆ G,i (k, k) = G n * x -n x,i +k-1
For similar reasons as for

Q G R G (k) = G -1 δ G R k δ G (A.37)
where R k is a symmetric positive semi-definite matrix, and δ G is a diagonal matrix with elements determined by taking the first element of each subdiagonal of ∆

G for i=1...p δ G (i, i) = ∆ G,i (1, 1)
The application of R G (k) to an asynchronous measurement model is detailed in Data fusion, Sec. A.5. The innovation at time t is calculated in a moving window from (t -d) to t, and as such d is the characterizing parameter. Innovation for a variable-step asynchronous continuousdiscrete system is written as

I d,k (t) ≡ k i=m+1 ∆t i ||γ(y i - ∼ y i )|| 2 R ny (A.38) where ∆t i = t i -t i-1 (A.39)
and m = min

tm∈[t 0 ,t k ] (t m ≥ t -d) (A.40)
and γ is a normalizing matrix.

γ = diag(γ 1 , γ 2 , . . . , γ ny ) (A.41)
and ∼ y i is found by integrating the IVP presented in eq. (A.33)

       ż = Az + b(z) ∼ y i = C i z(t i ) and z(t m ) = ẑ(t m ) (A.42)
In other words, the innovation is the sum of the length of the time interval, ∆t i , times the difference between the measured state and the innovative predicted state, γ||y - Innovation can be thought of as an experiment to see what happens if the observer correction is turned off at time t -d and the system is allowed to continue naturally. The amount of innovation is, roughly speaking, proportional to the square of the norm of the error between the observed trajectory and the non-observed trajectory.

A.4.4.2 Understanding "adaptation"

In a classical High-Gain Kalman Filter, the high-gain parameter, G is fixed. In an AEKF, a function called adaptation varies the value of G on the interval [1, G max ]. G = 1 and G = G max correspond to the unity-gain mode, otherwise known as the classical EKF, and the high-gain EKF mode, respectively. The adaptation, F, is driven by the estimatation quality measurement, the innovation.

In the CD-AEKF, the adaptation can be used as a threshold-based switch which incrementally steps G up and down between 1 and its maximum value. As explained in [START_REF] Boizot | Adaptive High-gain Extended Kalman Filter and Applications[END_REF], an incremental function is suitable in the continuous-discrete case.

F is used in order to limit the overshoot. It is known that for high-gain filters, there is a peak phenomenon [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF] when the adaptation switches. This phenomenon creates a perturbation to the covariance matrix P , which resembles a poorly initialized P 0 .

Note 21

The model described by (1.16) is driven ONLY by perturbations. Thus the highgain triggering frequency will be high, and as such the goal of the adaptation is to reduce to a maximum the peak phenomena.

An example adaptation function

F(G k-1 , I d,k ) would be If (I d,k > upper threshold) AND (G < G max ) ⇒ (G = G + 1) If (I d,k < lower threshold) AND (G > 1) ⇒ (G = G -1)
Preliminary results with the AEKF are promising as the high-gain component allows the observer to converge much more quickly when there is an unmodeled perturbation, such as a sudden braking or acceleration.

A.5 Data fusion

One of the particular problems that can occur in filtering is treating asynchronous and redundant data. In a multiple-output system, diverse sensor measurements are combined together in order to form the state estimation. Some sensor units might give multiple sensor readings, such as a GPS that in addition to the desired latitude and longitude outputs gives altitude, speed, bearing, etc... Furthermore, in a discrete-time sampled system, it is often the case that various sensor data arrives asynchronously, and this asynchronous data need not necessarily have a constant sampling rate.

These situations are both handled well by parametric model-based filters (such as the Kalman Filters presented above).

A.5.1 Redundancy

In general, sensors are chosen such that each sensor gives different information from the others, however this is not always the case. When more than one sensor measures the same data, this data becomes redundant. One way to handle the redundant data measurements is to simply ignore them, but this is less than ideal, especially as the redundant case is easy to treat.

Each redundant data source acts like a new output in the measurement model H(x, u). Each source has an associated covariance value (in the R matrix) that reflects the confidence in each particular measurement, exactly as would be done if the data were not redundant.

For instance, take the case of a system receiving altitude information from both a GPS and a DEM (Digital Elevation Map). A GPS is not very accurate in altitude terms, and tends to get a "stuck" on a particular bias, so the measurement covariance parameter for the GPS would be tuned to a relatively high value. On the other hand, a DEM is often accurate to sub-cm resolution, so it would have a relatively low covariance value.

The advantage of data redundancy is that if one sensor or another is lost, the model still retains its good qualities, and, more importantly, cannot lose observability. Again, drawing on the above example, there are times when a GPS cuts out-urban canyons, tunnels, etc..., and times when a DEM does not have any relevant data. The combination of the two "sensors" gives us the best possible results.

Redundancy can occur in any time domain. Furthermore, with certain sensors the information is complementary, e.g. digital map data is more precise, but GPS data has a finer grid resolution.

A.5.2 Asynchronous data

Asynchronous measurements, by their very nature, can only occur in discrete update time domains, e.g. continuous-discrete and discrete-discrete. There have been many papers written about data fusion of asynchronous signals, some which recommend interpolating lower frequency signals in order to form a sort of "oversampling", or resampling higher frequency signals in order to form a sort of "undersampling". The author recommends against this approach, as not only can it introduce unwanted effects, it completely breaks down if some of the measurements are not only asynchronous but also do not have a consistent sampling time.

The asynchronous case is not as straightforward as the redundancy case. Some measurements can arrive as a function of time, some as a function of a state variable, and some as externally triggered events. (See Note 18 for a discussion of asynchronous system observability.)

Such a set of measurement data necessitates an approach with a variable time step, where a new update/correction is triggered at each moment there is a new data value. This implies a constantly changing measurement matrix H k .

As shown in Algorithm 1, we integrate the differential system equations until we encounter an event at time t k . At this moment, we start by determining what measurements are present and then stack them into the output vector, y k . Then we assemble the measurement Jacobian matrix H k in the same fashion, again by stacking the individual row vectors associated with the measurements.

Once we have assembled H k , we can proceed as normal with the discrete update/correction of the observer.

Additionally, the measurement covariance matrix R k must change as a function of the time step and of the available data. It is always of dimension R ny×ny Dynamic optimization problems can be solved with different approaches, depending on the particular nature of the problem. Globally, optimization problems can be tackled along three major lines, presented in descending order of computational costs:

• Dynamic programming: uses the principle of optimality of subtrajectories in order to recursively compute a control for any x 0 and any t. The solution for a continuous problem stems from resolving the Hamilton-Jacobi-Bellman (HJB) equation, a partial differential equation in state space.

• Direct methods: referred to as "first discretize, then optimize". These methods require discretizing the commands u along multiple subsections of time, and then seeking the optimal trajectory across each subspace. The direct shooting and direct multiple shooting methods are examples of this approach.

• Indirect methods: referred to as "first optimize, then discretize". They have the advantage of being able to present a complete optimal synthesis of the problem, dramatically reducing computational costs. The shooting and multiple shooting methods are examples of this approach.

B.1 Controllers and controllability B.1.1 Controllers

Roughly speaking, the control is the vector u that brings the state from x(0) to x(t f ) in finite time t f < T .

Thus, the controller is the mathematical algorithm that calculates the control u.

B.1.2 Controllability

Controllability, like observability, is an important property for any system we wish to completely control. Controllability refers to the ability to act on all the variables required to control the system. The world is full of systems that we can observe, but not act upon: a bouncing ball, a paper airplane, a waterfall, etc... What is missing in all these cases is that we have no external actor, and thus they are not controllable.

Controllability is a mathematical way of expressing the idea that if a system starts at A then there is a sequence of commands, out of all the possible commands, that will bring the system to state B.

In more rigorous terms, a system is said to be controllable if for all vector pairs (x a , x b ) in R n there exists a finite time T and a command u defined on [0, T ] such that, when the command u is applied, the solution x(t) starting from x(0) = x a gives x(T ) = x b . Non-linear controllability: In the nonlinear case, it can happen that the system is controllable at certain points, but not at others. We can speak about a larger condition of global controllability. For an affine system of the form

ẋ = f (x) + n u i=1 g i (x)u i
where g i is a vector valued function.

We have a test for local controllability, but similar to the observability test in Sec. A.1.2, we must use Lie algebra to apply it.

First we define the Lie bracket operator (The derivation and purpose of Lie brackets is far outside the scope of this PhD. If the reader would like further information, s/he is encouraged to read Boothby's excellent description. [START_REF] William | An Introduction to Differentiable Manifolds and Riemannian Geometry[END_REF]):

[f, g] ≡ ∂g ∂x f - ∂f ∂x g
and higher order Lie brackets:

(ad 1 f , g) ≡ [f, g] (ad 2 f , g) ≡ [f, [f, g]] . . . (ad k f , g) ≡ [f, (ad k-1 f , g)]
Note: ad is read "adjoint".

The system is locally controllable about a point if and only if the controllability matrix

g 1 . . . g nu [ad 1 g i , g j ] . . . [ad nu g i , g j ] . . . [ad 1 f , g i ] . . . [ad nx-1 f , g i ] (B.2)
has rank n x .

B.2 Optimal control

We might ask ourselves the question, 'If a system is controllable, is there a such thing as a "perfect" control?' First, we have to define what a "perfect" control is. One popular choice is to define "perfect" as "minimizing a cost". (Or "maximizing a benefit". The opposite of one is mathematically identical to the other.) Thus, among all the possible solutions that bring x(0) to x(T ), we look to find a trajectory u * that minimizes a cost, J(T, x, u). Such a trajectory, if it exists, is called an "optimal control" for this cost.

The existence of an optimal trajectory is dependent on the system regularity and the cost function. It is possible that an optimal control does not exist in the class of considered controls, but that it exists in a larger space of controls. [START_REF] Trélat | Contrôle Optimal: théorie et applications[END_REF]

B.2.1 Concept of time in optimal control

In [START_REF] Bellman | Dynamic Programming[END_REF] Bellman showed that the natural approach to optimality problems is to work backwards in time from the solution. Call this time to go.

Time to go makes a certain amount of sense when it is considered that the objective is known, but not necessarily the initial point. While we can easily calculate the finite number of ways to optimally approach the final point, there are an infinite number of ways to leave the initial point, and we have no way to immediately distinguish those that are possibly optimal from those that are not.

Thus, we talk about integrating from the final target x(T ) to the initial point x(0). One confusing facet of inverse time is how to refer to it. If one talks of time T in an inverse time problem, does this refer to time T with respect to forward or backward time? For this reason, we prefer the phrase "time to go", which, when scrupulously used, removes all ambiguity. We represent "time to go" by s and thus

t = 0 → s = T t = T → s = 0
Inversing the dynamic equations is straightforward. It suffices to take the negative of the dynamic equation ẋ * =ẋ (B.3)

B.2.2 Linear quadratic control

When the system is linear, or is linearized about an operation point, a special general solution can be found in the case where the cost is of the form

J = 1 2 t f 0 (x T Q lqr x + u T R lqr u)dt (B.4)
where Q lqr and R lqr are positive semi-definite and positive definite constant matrices, respectively.

Note 23 Since the solution has already been found for all linear systems with the cost as described in eq. (B.4), the principle of time to go is not used here.

B.2.2.1 Continuous LQR

The optimal control is then

u * = -K lqr x (B.5)
where

K lqr = R -1 lqr B T S lqr (B.6)
and where S lqr is the convergence of the differential Riccati equation Ṡlqr = -S lqr A -A T S lqr + S lqr BR -1 lqr B T S lqr Q lqr (B.7)

B.2.2.2 Discrete LQR

Perhaps a more practical use for the LQR is optimal control in discrete-time. We present here the slightly modified approach for the discrete LQR.

S 0 = Q dlqr + A T k (B k R -1 dlqr B T k ) -1 A k S k = Q dlqr + A T k (S -1 k-1 + B k R -1 dlqr B T k ) -1 A k K dlqr,k = (B T k S k B k + R dlqr )(B T k ) -1 S k A k (B.8) (B.9)
where Q lqr and R lqr are as in the above case for the continuous control. The control is as before:

u * = -K dlqr,k x k
The LQR is particularly interesting for basic control, as it provides a first point of entry to optimal control. If the model is known, then the Jacobian can be derived (or calculated) and as a consequence it is computationally cheap to create an optimal control that has sufficient degrees of freedom for minimizing the inputs or the state values.

The LQR control is easy to apply, and far outperforms rudimentary controllers such as PIDs. If, however, the cost that we seek to minimize cannot be effectively approximated by eq. (B.4), we are forced to explore higher complexity optimization solutions.

B.2.3 Dynamic Programming

Dynamic programming is a design technique similar to "divide-and-conquer". Divideand-conquer algorithms partition the problem into independent subproblems, solve the subproblems recursively, and then combine their solutions in order to solve the original problem.

Thus, dynamic programming is an optimality approach where a complex optimization problem is broken down into a sequence of decision steps over time. It uses the princple of optimality of subarcs in order to recursively compute feedback control for all time t and all x 0

In the continuous-time case, this leads to the Hamilton-Jacobi-Bellman (HJB) equation.

∂S ∂t

+ H 1 x, ∂S ∂x = 0
Dynamic programming, while "optimal" in the sense that it eventually finds an optimal solution, is resource intensive. It is also not well adapted to high accuracy, and it is possible that it miss singular optimal paths. The allure of dynamic programming is its programming simplicity for complex models with arbitrary initial states. where the finite dimensional vector w represents the optimization degrees of freedom, and with the differentiable scalar function a and vector valued functions b, c NLP, while widely used in industry, does not give any feel to the problem dynamics. Moreover, while the computation technique is more graceful than dynamic programming, it still is a process that requires intense usage of computer memory and processor. Moreover, there is no way with NLP to be sure that the optimal solution as given by the solver is the true optimal solution.

B.2.4 Non-linear Programming

B.2.5 Pontryagin's Maximum Principle

Note 24 The PMP is often used as Pontryagin's Minimum Principle, which simply inverses Pontryagin's original "benefit function" to be maximized into a "cost function" to be minimized. The mathematics are the same, only the sign of the cost J is different.

PMP states a necessary condition that must hold on the entire length of an optimum trajectory, starting from a fixed initial state, x(0). [START_REF] Trélat | Contrôle Optimal: théorie et applications[END_REF] Consider the system

ẋ = f (x(t), u(t)) (B.11)
with associated cost function

J(x, u, t) = K(x(T )) + T 0 L(x(t), u(t), t)dt (B.12)
where L the instantaneous cost, and K the terminal cost.

We can define the Hamiltonian as

H(x, u, λ) = λ T f (x, u) + λ 0 L(x, u, t) (B. 13
)
where λ is the costate vector, also known as the adjoint vector, and λ 0 ≤ 0 (since this is Pontryagin's Maximum Principle. In the case of the Minimum Principle, λ 0 ≥ 0) such that the vector < λ 0 , λ > = 0

Note 25 λ 0 is generally normalized to -1 by dividing the costate vector λ by λ 0 . However, in the case where λ 0 = 0 this obviously cannot be done. This case arises when the optimal trajectory is independent of the cost, which can occurs, for instance, when there is only one possible trajectory.

Suppose that x * (t), u * (t), λ * (t) represent the optimal state, optimal control, and optimal costate trajectories, respectively. I.e. u * (t) is the control that that brings the state from an initial state x(0) to a final state x(T ) and minimizes the cost J, and x(t) is the associated state trajectory. Then there exists a costate trajectory λ(t) such that ẋ = ∂H ∂λ

λ = - ∂H ∂x (B.14)
and for almost all t,

0 ≤ t ≤ T H(x * , u, λ) ≤ H(x * , u * , λ) (B.15)
that is to say that the optimal control u * is the control that maximizes H amongst all the possible controls u that can be applied to bring x(0) to x(T ). By using the system and costate vector, we can define the augmented system z:

ż = ẋ λ (B. 16
)
The PMP is particularly adept at solving minimal time problems (when the cost function depends only on time), but can also be solved for fixed time, or can be solved independent of time. If the fixed time limit is set far enough in the future, this will yield identical results to the solving independently of time.

PMP can give a set of global rules that can be expressed in analytic terms, called an optimal synthesis. Such an optimal synthesis can easily be programmed into a microprocessor, greatly reducing the real-world resources necessary to compute the optimal control. These rules read like an instruction manual-they dictate what control should be applied and when, much like a numerical recipe. See Appendix G for a worked example.

B.2.5.1 Transversality conditions

Transversality on the costate vector If not all states at the final state x(T ) are fixed, then the costate vector can be such that the transversality conditions at one or both endpoints satisfy for all σ such that x + ǫσ is within o(ǫ) of the termination point of an optimal trajectory for all sufficiently small ǫ > 0.

λ T 0 x(0)σ = 0 (B.
In other words, these two conditions verify that the costate vector is perpendicular to the state vector at the endpoints.

Note 26

If there is no value-added cost K, this condition has an immediate consequence when a state x j is not present in the optimal target. The transversality condition gives λ j = 0.

B.2.5.2 Limits on the control

Limits on the control pose no particular problem. Indeed, this sort of control limit oftentimes presents itself in the form of bang-bang controls, controls that are applied only at the upper or lower bound of the possible control.

B.2.5.3 Limits on the state

Limits on the state are much more complicated to deal with. While there are several versions of the Maximum Principle that allow for state limits-see for instance [START_REF]The maximum principle under minimal hypotheses[END_REF], the theory and application is sufficiently complicated that satisfactory results can be more easily obtained with other methods. The problem with state limits comes from the fact that the costate vector can become discontinuous at the state boundaries.

Oftentimes, a state limit is avoided by the simple nature of the problem. For instance, it is unlikely that a fixed-time energy optimal trajectory for a car trip will push the car to break the speed limit. If it does, the fixed-time horizon can simply be elongated, which will almost certainly find a solution at a slower speed. It is in general easier to examine the PMP optimal solution and modify the input parameters in case of a state limit, than it is to treat the state limits analytically.

Penalization methods: Another method to manipulate the optimal solution so as to impose state limits it to us penalization methods. By weighting certain limits, the optimal solution can be forced to obey the state limits. For instance, suppose that we want to impose a limit on the state. We can add a cost of the form α T 0 g(x(t))dt to the original cost function J(x, u, t), where α is the weighting value, and g is a function that is 0 when the state is within its limits, and strictly positive everywhere else. The hope is that by adding this cost the the optimal solution will force the optimal trajectory to stay within the prescribed limits.

The mathematical justification for this approach can be found in [START_REF] Trélat | Contrôle Optimal: théorie et applications[END_REF].

B.3 Numerical methods for optimal control

Choosing the appropriate method is based on questions of numerical stability, computation time, and desired accuracy. While this dissertation focuses on Pontryagin, we also present here strategies for Non-Linear Programming and Dynamic Programming problems.

B.3.1 Indirect methods

Indirect methods are used for Pontryagin problems where the costate vector λ is also sought. Their advantages are a higher level of precision, converge toward the true optimal solution without risk of getting stuck at local minima, and are memory efficient. Lastly, the multiple shooting method is parallelizable, for increased computational speed.

B.3.1.1 Shooting method

The shooting method is a way of reducing a boundary value problem (BVP) to an initial value problem (IVP). Its name is inspired from the problem of shooting at a known target. It is not known in advance how high to aim the shot, so one simply aims, fires, and sees where the shot lands. If it lands short, the next is aimed higher. Conversely, if it lands far, the next is aimed lower. The principle is to take the first boundary condition as an initial point and then create a second condition stating the gradient of the function at the initial point. The IVP is now completely well-posed and the system can be integrated. If the system arrives at a final value other than the target, the initial gradient condition is adjusted. This loop is iterated until the integrated value converges on the true target to within an acceptable amount of error. The solution to the IVP is now also the solution to the BVP.

In mathematical terms, we have a BVP of the form

ż(t) = f (t, z(t)) r(z(0), z(t f )) = 0 (B.20)
where r(x, y) stands for "boundary conditions at x and boundary conditions at y", and where

z = x λ
is the augmented state.

We transform this into an IVP by supposing z 0 = z(0) and

ż(t) = f (t, z(t)) r(z 0 , z(t f , z 0 )) = 0 (B.21)
We now have a classic IVP where we seek z 0 such that the second half of eq. (B.21) holds true. Effectively, we are looking for the roots, and this can be solved by a classical zero-crossing method, for example Newton or bisection. With respect to the simple shooting method, the multiple shooting method breaks the interval [0, t f ] into N intervals [t i , t i+1 ] and assigns an initial guess values z(t i ) to the beginning of each interval. The interest is to improve the stability of the shooting method.

The application of the maximum principle as applied to the optimal control problem ż as described in eq. B.20 reduces to

ż(t) = f (t, z(t)) =          f 0 (t, z(t)) if t 0 ≤ t 1 f 1 (t, z(t)) if t 1 ≤ t 2 . . . f s (t, z(t)) if t s ≤ t f (B.22)
where t 1 , t 2 , . . . , t s ∈ [t 0 , t f ] can be commutations of the optimal control and, in the case of state constraints, contact with the state frontier. Continuity is be preserved for z at the commutations, but the costate can jump at the frontiers in state constraint cases. We begin by defining a function Φ(t) that is valid on t = [0, T ] (a line connecting each state initial value with its final value is a good first guess, but in fact, it can be any Lipschitz continuous function). An initial guess of unknowns z 0 is made. Then system (B.22) is integrated until ||z(t, z 0 ) -Φ(t)|| > ε for some ε > 0. We designate the time variable at this point as t 1 . Now the integration of the system continues from t 1 using Φ(t 1 ) as the initial 'guess' for the solution. This process continues until time t = b is reached. Now the error function e

(s) = ||z(t i ) -z i || is formulated where s = z 0 z 1 . . . z N -1 T ,
where z i is the initial state for the trajectory on the interval [z i , z i+1 . After this is accomplished, a correction is made to z 0 using a classical zero-crossing method, for example Newton or bisection. This process is repeated, only now the number of iterations N is fixed, and the function Φ(t i ) is replaced by z i . The iterations stop once the global convergence is satisfactory.

Figure B.3 shows the multiple shooting method for the simple case of a one-dimensional system.

Note 28

The multiple shooting method improves stability by increasing the number of nodes over the simple shooting method. While there are many more unknowns in the multiple shooting method, the method is parallelizable. [START_REF] Burrage | Parallel and Sequential Methods for Ordinary Differential Equations[END_REF]

B.3.2 Direct methods

Direct methods are used for Non-Linear Programming models. These solutions are less accurate than the indirect methods, but are oftentimes easier to program, trading computer time for programming time. They also easily handle state constraints and are not very sensitive to the initial value guess.

B.3.2.1 Direct single shooting

Generally, direct shooting methods are the the easiest when a controls problem is first approached. The strong points of direct single shooting are (1) that it can use fully adaptive, error controlled, state-of-the-art ODE or DAE solvers; (2) that it has only a few optimization degrees of freedom even for large ODE or DAE systems; and (3) that only initial guesses for the control degrees of freedom are needed.

The weak points are (1) that the ODE solution x(t; q), where the arguments after the semi-colon (;) denote the dependence on the interval's initial value and control, can depend (2) that unstable systems are difficult to treat; and (3) it tells us little about the nature of the costate vector. [START_REF] Diehl | Fast direct multiple shooting algorithms for optimal robot control[END_REF] Due to its simplicity and general efficiency, the direct single shooting approach is often used in engineering applications, e.g. gOPT. [START_REF] Pse | gPROMS Advanced User's Manual[END_REF] The direct single shooting approach starts by discretizing the controls into N subintervals. We choose grid points on the interval t

0 = τ 0 < τ 1 < τ 2 < • • • < τ N -1 < τ N = t f
and chose the controls u piecewise constant, u(t) = q i for t ∈ [τ i , τ i+1 ] so that u(t) depends only on the finitely many control parameters q = {q 0 , q 1 , . . . , q N -1 , T } and can be denoted by u(t; q), again where the arguments after the semi-colon (;) denote the dependence on the interval's initial value and control. If that problem is for a fixed time T , the last element of q disappears, as T is no longer an optimization variable.

Path constraints are discretized in order to avoid a semi-infinite problem, for example by requiring that the constraints hold only at grid points τ i . Thus the optimal problem reduces to the following finite dimensional nonlinear programming problem (NLP): min q T 0 L(x(t; q), u(t; q))dt + K(x(T ; q) (B.23) subject to path and terminal constraints and to

s 0 -x 0 = 0
This problem can be solved by a finite dimensional optimization solver, e.g. sequential quadratic programming (SQP). [START_REF] Nocedal | Numerical Optimization[END_REF]

B.3.2.2 Direct multiple shooting

The direct multiple shooting, as developed in [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF], combines advantages of simultaneous methods with the direct single shooting approach. In direct multiple shooting, we proceed as follows. First, as in the case for direct single shooting, we discretize the controls piecewise on a coarse grid.

u(t) = q i for t ∈ [τ i , τ i+1 ]
where the intervals can be as large as in the direct single shooting case. Second, we solve the ODE on each interval [τ i , τ i+1 ] independently, starting with an artificial initial value

s i ẋi = f (x i (t), q i ) for t ∈ [τ i , τ i+1 ] x i (τ i ) = s i
By numerical solution of these IVPs, we obtain trajectory arcs x i (t; s i , q i ) Simultaneously with the decoupled ODE solution, we also numerically compute the cost function integrals

l i (s i , q i ) = τ i+1 τ i L(x i (τ i ; s i , q i ), q i )dt
In order to constrain the artificial degrees of freedom s i to physically meaningful values, we impose continuity conditions s i+1 = x i (τ i+1 ; s i , q i ) Thus we come to the following NLP formulation that is completely equivalent to the direct single shooting NLP, but contains the extra variables s i and has a block sparse structure. min s,q

N -1 i=0 l i (s i , q i ) + K(s N ) (B.24)
subject to path and terminal constraints and to

s 0 -x 0 = 0 s i+1 -x i (τ i+1 ; s i , q i ) = 0
If we summarize all variables as w = (s 0 , q 0 , s 1 , q 1 , . . . , s N , q N ) we obtain an NLP which can be exploited in the tailored SQP solution procedure. [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF]

B.3.3 Dynamic programming methods

The Hamilton-Jacobi-Bellman equation (HJB) is

∂S ∂t + H 1 x, ∂S ∂x = 0 (B.25)
where

H 1 = max u∈U λ T ẋ -L(x, u)
There are a number of discretization methods. The scheme for finite difference discretization proposed below is one of the simplest, but higher order schemes can be used. We start by writing a finite difference method simple discretization of eq. (B.25). Remark that in order to ensure stability, the discretization scheme must be decentered. Thus, in order to discretize Let us consider a mesh in space (x i ), where i = (i 1 , . . . , i n ) ∈ Z, which we will suppose to be regular in order to simplify matters, and as consequence a regular discretization of time, (t j ). Let h = (h 1 , . . . , h n ) be the spatial step size, and k = t j+1 -t j be the temporal step size.

Let S i be the approached value of S(t j , x i ). It suffices to approach ∂S ∂x p (x i ) by a difference divided to the left if ẋp (x i , u) > 0, and to the right if ẋp (

x i , u) < 0 For a ∈ R a + = max(a, 0) = a+|a| 2 a -= min(a, 0) = a-|a| 2 (B.26)
For all p ∈ {1, . . . , n} we note êp = (0, . . . , 1, . . . , 0) where the 1 is at the p th place. We thus obtain the explicit scheme

0 = S i,k+1 -S i,k k +max u∈U   n p=1 S i,k -S i-êp,k h p ẋp (x, u) + + S i+êp,k -S i,k h p ẋp (x, u) --L(x i , u)   (B.

C.1 Overview

This appendix details the various data formats used in this dissertation. Specifically

• Quad-coordinates: database specification for quad-tile coordinate system (Sec. C.1.1)

• Configuration: the XML formats used to characterize cars and trips (Sec. C.1.3)

• Observer files manifest: the files generated by the observer program (Sec. C.1.2)

• Logging: the protocol in which the microcontroller code saves the data (Sec. C.1.4)

C.1.1 Quad-tiles and quad-coordinates

While latitude-longitude is the traditional way to geolocate objects, it is not a format that leads itself to easy and rapid search algorithms with fixed-point microcontrollers. A better system is quad-tile coordinates.

Quad-tiles split the world into 4 tiles at each zoom level. Zoom level 1, the entire Earth, has four tiles made by two splits: one along the equator and another along the prime meridian (assuming a Mercator projection where the world is already split in half along the 180 • meridian). See Subsequent tiles can be cut into four smaller tiles, ad infinium. The tile information can be represented as a series of pairs of 1's and 0's. In each pair, the first bit, the MSB, represents the left/right side(0=left, 1=right) and the second bit, the LSB, represents the top/bottom (0=top, 1=bottom). Thus each zoom level requires 2 bits.

As such, the entire Australian continent is contained in the quad-tile 11, all of Greenland is contained in the quad-tile 00, and the South American continent is split between quad-tiles 01 and 00. See With a 32-bit integer (zoom level 16), the tiles are < 600m on a vertical side (where obviously the horizontal width will be a function of latitude), which is already a quite small bounding box.

With a 64-bit integer (zoom level 32), the tiles are < 1cm on a side, which gives a very high degree of precision for geolocalization.

Thus quad-tiles lead themselves to geolocalization because these base integer sizes are (luckily) already well-suited to the problem.

Note 29

The attentive reader will remark that quadcoordinates are incapable of representing an exact point. Furthermore, any point falling on the edge of one of the two bisecting axes at zoom level 1 cannot even be contained in a quad-tile. Thus, both the Equator and North Pole pose representational problems, even if from an engineering perspective this is not relevant.

C.1.1.1 Quad-tiles database entry format

A system of roads can be seen as an ordered list of interconnected nodes. Nodes can be shared by multiple roads. The database storage and search algorithm used in this dissertation are detailed below:

Node entry: Each node ID is the same as the quad-tile coordinate for that same node. Trading storage space (cheap) for processor speed (expensive), each node also has its associated latitude-longitude coordinates. The number of intersecting ways is provided, and lastly the 32-bit ID for the n ways that intersect the node are given. Quad-tile bounding box entry: The quad-tile bounding box entry indicates all the ways that cross the bounding box. This is a useful representation for finding all the ways close to a given point. As in the case for the node entries, the quad-tile bounding box's ID is also its coordinate. The zoom level is also given (by necessity, as the 64-bit field 1100000000000000000000000000000000000000000000000000000000000000 at zoom level 1 would give one fourth the planet, whereas at zoom level 32 it would represent a square roughly one cm on a side almost exactly at the intersection of the equator and the prime meridian.), followed by the number of ways and finally by the IDs of all the contained ways. <ty pe ty pe=" 698 c c I n l i n e 3 Turbo Suprex " /> <d i s p l a c e m e n t d i s p l a c e m e n t=" 698 c c " /> <power power=" 60kW @5250rpm" /> <t o r q u e t o r q u e=" 110Nm @4500rpm" /> <c o m p r e s s i o n comp=" 9 : 1 " /> < i d l e rpm=" 850 " /> < i n e r t i a lumped=" 0 . 0 6 1 " /> </ e n g i n e> <w h e e l s> < t i r e summer=" 205/45 R16" /> < i n e r t i a lumped=" 0 . 6 2 " /> < !--I n e r t i a PER w h e e l / t i r e a s s e m b l y--> </ w h e e l s> <a c c e l e r o m e t e r ty p e="LIS3LV02DL"> < c a l i b r a t i o n> < r o l l p h i=" 2 . 5 8 2 8 7 2 3 7 2 3 7 8 4 2 " /> <p i t c h t h e t a=" 0 . 0 2 7 7 7 4 4 0 1 4 7 5 6 9 7 5 " /> <yaw p s i=" 1 . 6 6 3 4 9 3 0 5 1 0 8 4 1 3 " /> < s c a l e s c a l e=" 0 . 9 8 9 6 5 6 4 0 7 1 6 5 8 3 9 " /> </ c a l i b r a t i o n> </ a c c e l e r o m e t e r> <t r a n s m i s s i o n t yp e="6-spd s e q u e n t i a l "> < !--These a r e t h e f i n a l d r i v e r a t i o s b e t w e e n t h e e n g i n e and t h e w h e e l a x e l--> <r e v e r s e g e a r=" " /> < f i r s t g e a r=" 1 3 . 2 5 " /> <s e c o n d g e a r=" 9 . 6 0 4 " /> <t h i r d g e a r=" 7 . 3 7 " /> <f o u r t h g e a r=" 5 . 7 9 6 " /> < f i f t h g e a r=" 4 . 2 0 8 " /> <s i x t h g e a r=" 3 . 2 3 " /> </ t r a n s m i s s i o n> <drag> <a e r o Cd=" 0 . 3 9 " /> < r o l l i n g Crr=" 0 . 0 1 2 5 " /> </ drag> <a r e a f r o n t a l=" 1 . 9 3 " /> <mass wet=" 844 " /> < i n j> < c a l i b r a t i o n f u e l P e r S e c o n d=" 0 . 0 0 8 " /> < c a l i b r a t i o n l a g=" -0.0004 " /> < !--The d i f f e r e n c e b e t w e e n t h e a c t u a l i n j e c t o r o p e n i n g t i m e and t h e i n j e c t o r v o l t a g e p u l s e l e n g t h , i n [ s ] .--> <t i m i n g p a t t e r n=" s e q u e n t i a l " /> </ i n j> <v s s> < c a l i b r a t i o n m e t e r s P e r T i c k=" 0 . The CAN bus (a.k.a. CANbus), is a fault-tolerant, bi-directional bus that was originally designed for automotive use, but has since become a quasi-standard in industry and noisy environments. CAN is mandatory in all European vehicles since 2004.

While CAN is technically only a bit transmission protocol, in the automotive world the term is used interchangeably to mean ISO 11898-2, which completely defines the physical layers, including signal levels, wire count, etc...

In the case of ISO 11898-2, there are two wires, both transmitting the same bit. The wires pass between dominant, 0, and recessive, 1, states. The dominant is so-named since if both a dominant and recessive bit are sent simultaneously, only the dominant bit is transmitted on the bus.

Technically, this is achieved by having an open-collector style bus, where the line is pulled HI (recessive) by a pull-up resistor, and then is pulled LO (dominant) by a transistor. The transistor obviously being far stronger than the pull-up resistor, when both a dominant and recessive bit are sent at the same time by two different devices, the transistor signal squashes the resistor signal.

A CAN message is composed of an arbitration field, a control field, a data field, and a CRC field. The arbitration field is designed such that if two CAN devices start speaking simultaneously, then the one with lower arbitration has priority.

(This occurs by the physics of the communication: if two arbitration fields are different, at some moment one device will transmit a 0 (dominant) when the other is transmitting a 1 (recessive). Since the 1 is recessive, all other devices on the bus will see a 0. As each CAN device is monitoring the bus as well as transmitting on it, as soon as a device recognizes that its recessive bit has been overwritten by a dominant, it stops transmitting. As such, for an 11bit arbitration field, 0b00000000000 is obviously always highest priority and 0b11111111111 is always lowest.) Convention states that the arbitration field is also the message ID. While this is by no means a requirement of the protocol, unless there is strong reason to suspect otherwise it can be assumed that an arbitration field is the same as the message ID.

The data field is always between 0 and 8 bytes long, as indicated by the control field.

D.2 Case study

Inspired by [START_REF] Schroeder | Smart roadster project: Setting up driveby-wire or how to remote-control your car[END_REF], and armed with the above information, a Smart Roadster's CAN bus was reverse engineering by logging all CAN messages, and then comparing with logged data from the sensor suite. There are many messages for which the meaning was never discovered, but they were unnecessary for our purposes. Some certainly indicate if headlights are on, if the convertible top is open, etc... (These could be later interesting for trying to gauge electrical load as a part of the overall engine output.)

A simple circuit was designed to interface to the automotive CAN bus. The circuit needs

• A CAN transceiver to turn dominant/recessive bits into 0's and 1's and vice-versa

• A CAN controller in order to read the bit stream, perform error checking, bit destuffing, etc...

In this case, we chose the STM32F103VET6, an ARM Cortex M3 microcontroller with onboard CAN controller, and coupled it to a MAX3051 CAN transceiver.

The CAN messages as decoded are below:

E.3 Improvements

If the reader would like to use these schematics directly, s/he would be well advised to take into account the following bugs and worthwhile upgrades:

• The ATMega has been superseded by the XMEGA line. These chips are in all ways superior; it is cheaper, faster, more accurate, and has more peripherals.

• In the functioning board, the transistor Q1 had to be replaced by an optorelay because the FV-M8 GPS had trouble initializing correctly. The optorelay has a faster rise time. It's primary disadvantage is cost, so with a different GPS the transistor might not need replacing.

• The CK trace on the accelerometer can be removed. It is not important that the accelerometer have an external syncing signal.

• The reset button should be reconfigured to use an RC filter so that a reset is delayed by several seconds, in order that the button can also be used for other purposes.

• The PPS input should have an RC filter A future version could/should include:

• A CAN bus and LIN bus transceiver, in order to be compatible with all foreseeable cars.

• An analog input that can be used for fuel pressure sensing (diesel cars)

• A DC-DC converter in order to draw less current from the car's battery (important when the car is off)

As a small bonus, a separate circuit trace was made under the SD card. This is specifically designed to be able to test the STMicro LIS331DL. The relevant area must first be cut out from the circuit board. The pad layout is designed to be soldered on top of the existing accelerometer area. In fact, any digital accelerometer can be installed into this area, as long as the output pads correspond.

In addition, there are a number of PA* pads on the board. These are unconnected ADC inputs, and can be used for future expansion, depending on needs. The Ball and Wheel system is made of a spherical ball with angular momentum I b rolling on the rim of a wheel with angular momentum I w . The wheel is driven by a torque source, τ w that is limited by τ w,min ≤ τ w ≤ τ w,max . The objective is to balance the ball at the unstable equilibrium point at the very top of the wheel. In order for the system to be controllable (see Sec. B.1.2), τ w,min < 0 and τ w,max > 0 must hold. The ball is acted upon by gravity and by friction with the wheel rim. We assume a no-slip condition between the wheel and the ball.

E.4 Sensor suite parts manifest

This appendix presents a two-variable system, where only the ball's position and translational velocity with respect to the wheel are modeled. In consequence, the wheel and the ball can have arbitrary angular velocities, as long as the ball does not physically move from the set-point position. In this example, we develop a full synthesis solution.

Optimal target: In this example, we are looking for the time optimal control. In other words, we want the sequence of controls that minimize the time required to go from an initial point to a final set-point. Note that this is a different optimization problem from the one we wish to solve for the car, which minimizes fuel; and is also different from the LQR (See Sec. B.2.2), which seeks to minimize the state and control deviations.

The time optimal cost is expressed as J(x, u, t) = G.2 Two-state example (We rapidly present here the Ball and Wheel dynamic equations, without modeling them. Interested readers are encouraged to read [START_REF] Ho | Controlling a ball and wheel system using full-state-feedback linearization[END_REF] for a detailed derivation using Euler-Lagrange.)

The state vector is

G.2.1 Linearization by feedback

Linearization by feedback is a technique which greatly facilitates designing the control. The idea behind this approach is to find a diffeomorphism and a state-feedback control law that transform the nonlinear system into a linear time-invariant (LTI) system, as these LTI systems have certain desirable properties of convergence and stability. A complicated nonlinear control can be easily designed on the linear system. We give here the equations for a full-state-feedback linearization.

G.2.2 Pontryagin's Maximum Principle

The primordial goal is to find the equations for the costate vector (a.k.a. adjoint vector ), λ T = (λ 1 . . . λ nx ), as these describe the control.

First, recall the Hamiltonian, eq. (B.13) H(x, u, λ) = λ T f (x, u) -L(x, u, t) (G.9) then combining (G.2), (G.7) (where f (x, u) = ẋ = ( θwb γu f ) T ), and (G.9)

H = λ 1 θwb + λ 2 γu f -1 (G.10)
En suite, we want the control that maximizes the Hamiltonian across the trajectory. We find this classically, by 1) testing for absolute extrema at the control boundaries and 2) looking for local extrema by taking the partial derivative of H with respect to u and setting equal to 0 ∂H ∂u f = 0 = λ 2 γ (G.11)

This reveals that there are no local extrema as a function of u f . Therefore, the absolute extrema are the only extrema. 1Recalling that the torque is limited by τ w,min ≤ τ w ≤ τ w,max (as described in the initial problem statement), we deduce that the u f that maximizes H is

u * f =
u f,min = u f (τ w,min , 0) when λ 2 < 0 u f,max = u f (τ w,max , 0) when λ 2 > 0 (G.12)

Note 33

We choose θ = 0 as it has a negligible effect on u * f . In fact, this step is not necessary if we were finding a numerical solution, but in order to have a simple synthesis we make this very minor simplification in calculating the control saturations.

Note 34 This type of control is called bang-bang. The control is either maximum or minimum, depending only on the sign of λ 2 . As such, it becomes paramount to know when λ 2 crosses zero. This bears repeating: here the most important information we seek from the costate vector is the zero-crossing time. Once we can express this analytically for any (controllable) state, then we are done and we have only to check that the trajectory is in fact optimal.

Next, recall the relationship between the Hamiltonian and the costate vector (eq. (B. 1. H = 0 (From eq. (B.19))

2. u f is bang-bang, i.e. it is either u f,min or u f,max 3. λ 1 is constant, and λ 2 is a line dependent on time, with slope -λ 1

Point (1): We know that the Hamiltonian is always 0, and we know that at t = T that u f (T ) = u * f . We also know that θwb (T ) = 0 and thus from (G.10) we can calculate that

λ 2 (T ) = 1 u * f
Unfortunately we cannot further use the transversality conditions, since-from Note 27 (pg. 112)-due to the fact that our final point is x = 0 we cannot conclude anything more from the condition.

Point (2): While we know what controls we will give, it still remains to be seen when we will give these controls.

Point (3): Unfortunately, we have no additional equations with which to fix C 1 , C 2 . We must therefore perform an exhaustive analysis of λ 2 in order to understand when to switch between controls.

Note 35

Keep in mind that Pontryagin's Maximum Principle is a necessary but not sufficient condition for optimality. The fact that we have found trajectories that maximize H does not, therefore, mean that they are necessarily optimal.

Inverse time and "time to go": Optimization problems are naturally considered in inverse time. This makes a certain amount of sense when it is considered that the objective is known, but not necessarily the initial point. While we can easily calculate the finite number of ways to optimally approach the final point, there are an infinite number of ways to leave the initial point, and we have no way to distinguish those that are possibly optimal from those that are not. Thus, we talk about integrating from the final target x(T ) to the initial point x(0). The author suggests that if the above is not clear, that it be reread, and kept in mind throughout the following analysis. It bears repeating: in optimal control, we integrate starting at the final target and only stop once we reach the initial target.

One confusing facet of inverse time is how to refer to it. If one talks of time T in an inverse time problem, does this refer to time T with respect to forward or backward time? For this reason, we prefer the phrase "time to go", which, when scrupulously used, removes all ambiguity. We represent "time to go" by s and thus Knowing that u f (T ) can be either positive or negative, we can directly infer the nature of λ * 2 -it is a line that has at most 1 zero-crossing, depending only on whether α and β have the same sign or not. As the two cases are symmetrical in the linear case, we will only treat the case where β > 0.2 If β > 0 and α ≥ 0, then from (G.17 And thus by rearranging (G. [START_REF] Novel | Commande linéaire des systèmes dynamiques[END_REF] we see that in this case the optimal trajectory lies on the parabola x * 1 = 1 2γu f,max (x * 2 ) 2 , x * 1 ≥ 0, x * 2 ≤ 0. This is half of the switching locus, the other half occurring when β < 0.

If β > 0 and α < 0, then from (G.17) λ * 2 has one zero-crossing at 0 = αs 0 + 1 u f,max

s 0 = - 1 αu f,max (G.20)
and thus by integrating the IVP (and using the final values of the segment from 0 to s 0 as the initial values for the segment from s ≥ s 0 )

u f = u f,max x * 1 (s) = 1 2 γu f,max s 2 x * 2 (s) = -γu f,max s    for 0 ≤ s ≤ s 0 u f = u f,min x *
1 (s) = γ( 1 2 u f,min s 2 + (u f,min -u f,max )s 0 s -1 2 s 2 0 (u f,min -u f,max )) x * 2 (s) = γ(-u f,min s + (u f,min -u f,max )s 0 )

   for s 0 ≤ s (G.21) Now we must perform the same analysis for β < 0, but thanks to the system's symmetry we find equivalent results, only negative.

It might not seem it, but we're almost done. Note that the first part of eq. (G.21) is identical to eq. (G. [START_REF] Novel | Commande linéaire des systèmes dynamiques[END_REF]. From this we can gather that the optimal control from any space intersects the switching locus and then commutes to follow the parabola to the final destination. (see Fig G .2)

G.3 Interpreting the synthesis

In other words, there are four possibilities:

• The initial point is on the switching locus to the right of the origin. In this case, there are no commutations to be made and the optimal control is to apply u * f until the system reaches 0.

• The initial point is above the switching locus. In this case, the optimal strategy is to apply u * f,min until the switching locus is reached, at which point the bang-bang control commutes to u * f,max , which is applied until the system reaches 0.

• The initial point is on the switching locus to the left of the origin. The exact same as when to the right, only the opposite control is applied.

• The initial point is above the switching locus. The exact same as when the point is above the switching locus, only the opposite sequence of controls is applied.

Thus we have the time optimal control. Incidentally, the optimal time T can be found by simply integrating from the final point x f = 0 to the initial point x 0 Note 36 What happens in the case where the optimal control is not perfectly applied? In this case, one can see that if the optimal control is misapplied, or the system has a perturbation, that this is equivalent to an IVP with a new initial value at the current system state. The system still converges to the set point, only an additional commutation will be required.

Fig. 1 Figure 1 . 1 :

 111 Figure 1.1: Reference frames

  [START_REF]The maximum principle under minimal hypotheses[END_REF]) and Ḟthrust = Ṡj = Ṡk = Ṗin = ωe = ḃi = ḃj = ḃk = 0 (1.17)

Figure 2

 2 Figure 2.1: MPGuino

Figure 2 . 2 :

 22 Figure 2.2: GPS repeatability errors (Plotted with GET, Appendix F)

Figure 2 . 3 :

 23 Figure 2.3: Fuel injection schemes

Figure 2 . 4 :

 24 Figure 2.4: Fuel injector cutaway

Figure 2 . 5 :

 25 Figure 2.5: DEM vs. GPS along a trajectory

Figure 2 . 6 :

 26 Figure 2.6: Disagreement between DEM vs. GPS along a trajectory

Figure 3 . 1 :

 31 Figure 3.1: Sigma point averaging

Figure 3 . 2 :

 32 Figure 3.2: Quantized signal

Figure 3 .

 3 Figure 3.4 shows results gathered from data logged on a highway.

  DEM with GPS, without DEM without GPS and DEM (a) Efficiency vs. time DEM with GPS, without DEM without GPS and DEM (b) Torque vs. time

Figure 3 . 4 :

 34 Figure 3.4: Highway data

Figure 3 .

 3 Figure 3.5 shows results gathered from data logged in a city.

  with GPS and DEM with GPS, without DEM without GPS and DEM (a) Efficiency vs. time with GPS and DEM with GPS, without DEM without GPS and DEM (b) Torque vs. time

Figure 3 . 5 :

 35 Figure 3.5: In-town data

  with GPS and DEM with GPS, without DEM without GPS and DEM (a) Efficiency vs. time with GPS and DEM with GPS, without DEM without GPS and DEM (b) Torque vs. time

Figure 3 . 6 :

 36 Figure 3.6: Maximum acceleration data

Figure 3 . 7 :

 37 Figure 3.7: Efficiency map for Smart Roadster Coupé

Figure 3 . 8 :

 38 Figure 3.8: Sample efficiency map

2 eFigure 3 . 9 :

 239 Figure 3.9: 2nd order η fittings

e +α4ω 3 e +α5ωeT 3 eFigure 3 .

 33 Figure 3.10: 3rd order η fittings

5 ( 2 e + α5T 4 e 5 ( 4 eFigure 3 .

 524543 Figure 3.11: 4th order η fittings

Figure 3 . 12 :

 312 Figure 3.12: Example efficiency map with bounds

2 eFigure 3 . 13 :

 2313 Figure 3.13: 2nd order η -1 fittings

3 e 3 eFigure 3 .Figure 3 .Figure 3 .

 33333 Figure 3.14: 3rd order η -1 fittings

2 e 2 eFigure 3 . 16 :

 22316 Figure 3.16: 2rd order P in fittings

3 eFigure 3 . 17 :

 3317 Figure 3.17: 3rd order P in fittings

Figure 3 . 4 e

 34 Figure 3.18: 4th order fitted curve, P in = α 0 + α 1 ω e T e + α 2 ω 2 e T e + α 3 ω e T 2 e + α 4 ω 2 e T 2 e + α 5 T 3 e + α 6 ω 3 e + α 7 ω e T 3 e + α 8 T 4 e

Figure 3 . 19 :Figure 3 . 20 :Figure 3 . 21 :

 319320321 Figure 3.19: Efficiency map as calculated from dynamometer data. η = α 0 + α 1 ω e T e + α 2 ω 2 e T e + α 3 ω e T 2 e + α 4 ω 3 e + α 5 ω e T 3 e

Figure 4 . 1 :

 41 Figure 4.1: Engine operation zones

Figure 4 . 2 :

 42 Figure 4.2: Prototype Lexus RX 400h CVT

(

  ωe,Te)∈U adm H(x, u, λ) = 0 (a) Speed and Costate (b) Rpm and Torque (c) Hamiltonian (d) Control trajectory

Figure 4 . 3 :

 43 Figure 4.3: Optimization on the straight and level: CVT, T f non-fixed

Figure 4 . 4 :

 44 Figure 4.4: Classic multi-speed transmission

  (a) Speed and Costate (b) Transmission gear, Rpm and Torque (c) Hamiltonian (d) Control trajectory vs. efficiency

Figure 4 . 5 :

 45 Figure 4.5: Optimization on the straight and level: 6-speed transmission, T f non-fixed

  (a) Speed and Costate (b) Transmission gear, Rpm, and Torque (c) Hamiltonian (d) Control trajectory vs. efficiency

Figure 4 . 6 :

 46 Figure 4.6: Optimization on the straight and level: 6-speed transmission, T f fixed

Figure 4 . 8 :

 48 Figure 4.8: Slow acceleration test

Figure 4 . 9 :

 49 Figure 4.9: PMP optimal control acceleration test

Figure 4 .

 4 Figure 4.10: Optimization in two dimensions: CVT, T f non-fixed

  are the high-gain modified state noise covariance matrix and measurement noise covariance matrix, both as described in Sec. A.4.4. 2. I d,k is the innovation and F is the adaptation, both explained in Sec. A.4.4.3. b * (ẑ) is the Jacobian matrix (∂b(z)/∂z) | ẑ 4. γ is a diagonal matrix of normalizing constants, as explained in Sec. A.4.4

Figure A. 1 :

 1 Figure A.1: INNOVATION TIMESCALE

∼ y|| 2

 2 R ny , of each time interval that happened after t -d seconds. See Fig. A.1. (This point is justified by a Lemma in [7])

  Linear controllability: The Kalman rank condition states that a linear system of the form    dx dt = Ax + Bu y(t) = Cx + Du is controllable if and only if the controllability matrix B AB . . . A n-1 B (B.1) is of rank n x .[19]

  Non-linear programming (NLP) reformulates the optimal control problem into a finite dimensional model of the form min w a(w) subject to b(w) = 0, c(w) ≥ 0 (B.10)

17

 17

Note 27

 27 If the target x = 0, then we can only take σ = 0 and the transversality conditions do not add any additional information to the problem.Transversality on the Hamiltonian If the final time for arriving at the final state x(T ) is not fixed, there is a transversality condition at t = T H(x * , u * , λ * ) = -λ 0 ∂K ∂t (B.[START_REF] Novel | Commande linéaire des systèmes dynamiques[END_REF]

Figure B. 1 :

 1 Figure B.1: Shooting method

Figure B. 2 :

 2 Figure B.2: Shooting method

Figure B. 3 :

 3 Figure B.3: Multiple shooting method

  divided to the right or left, depending on the sign of ẋp

"

  Everything is practice."-Pelé

  Fig. C.1.

Figure C. 1 :

 1 Figure C.1: Quad-tile, zoom level 1

  Fig. C.1.

Figure C. 3 :

 3 Figure C.3: Quad search tree

Figure C. 5 :Figure C. 6 :

 56 Figure C.5: Accelerometer data format

Figure D. 2 :

 2 Figure D.2: CAN 2.0A ("standard CAN") Data Frame

Figure D. 4 :Figure E. 2 :Figure E. 3 : 2 Figure E. 5 :

 42325 Figure D.4: CAN Messages table

Figure G. 1 :

 1 Figure G.1: Diagram of Ball and Wheel

  Rτ w + mg sin(θ wb )(I w r 2 + I b R 2 ) I b I w + I w mr 2 + I b mR 2

1 I

 1 Let u f (τ w , θ wb ) = I b Rτ w + mg sin(θ wb )(I w r 2 + I b R 2 ) b I w + I w mr 2 + I b mR 2 (u f -mg sin(θ wb )(I w r 2 + I b R 2 ) I b R (G.8)

1 (G. 14 )

 114 Thus we can deduce that λ 1 = C 1 , and λ 2 = -C 1 t + C 2 Let's pause for a moment and take stock of what we know:

t = 0 16 ) 2 ( 0 ) = λ 2 (

 016202 → s = T t = T → s = 0 Inversing the dynamic equations is straightforward. It suffices to take the negative of the dynamic equation ẋ * =ẋ (G.15)This gives us the new extended systemAs we know x * (0) and x * (s = T ) this can be seen as a BVP (Boundary Value Problem) which can be solved, for example with shooting or multiple-shooting methods. (See Sec. B.3.) However, in this particular case we can analytically turn the BVP into an IVP (Initial Value Problem) in order to have a synthesis of what to do across the entire space of attainable points (the control convex).In other words, we create the synthesis by finding the solution for the exhaustive list of all possible IVPs. Consider that λ 2 * , like λ 2 , is a linear function of time-to-go of the formλ * 2 (s) = λ * 1 s + β = αs + β (G.17)where λ * 1 = Constant = α. We know that at time-to-go s = 0 that λ *

  ) λ * 2 (s) > 0 and has no zero-crossing for all time-to-go s. Therefor u * f (t) = u f,max and thus by integrating eq. (G.16), we can solve the IVP which yields ) = -(γu f,max )s (G.[START_REF] Novel | Commande linéaire des systèmes dynamiques[END_REF]

2 Figure G. 2 :

 22 Figure G.2: Optimal control synthesis for Ball and Wheel

  .3 Interpreting the synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Notation conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	List of Tables
	1	Abbr. µC	Meaning microcontroller
		µP	microprocessor
		CAN	Controller Area Network
		CVT	Continuously Variable Transmission
		CD-EKF Continuous-Discrete Extended Kalman Filter
		DD-EKF Discrete-Discrete Extended Kalman Filter
		DEM	Digital Elevation Map
		DP	Dynamic Programming
		ECU	Engine Control Unit
		EKF	Extended Kalman Filter
		PMP	Pontryagin's Maximum Principle
		UKF	Unscented Kalman Filter
		NED	North-East-Down coordinate system
	References
			xv

1 State of the (Theoretical) Art 1.1 State models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D CAN bus D.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.2 Case study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E Data logger circuit E.1 Circuit-board schematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.2 Circuit-board layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.3 Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.4 Sensor suite parts manifest . . . . . . . . . . . . . . . . . . . . . . . . . . . . F Google Earth G The Ball and Wheel: a worked example for Pontryagin's Maximum Principle G.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.2 Two-state example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.2.1 Linearization by feedback . . . . . . . . . . . . . . . . . . . . . . . . . G.2.2 Pontryagin's Maximum Principle . . . . . . . . . . . . . . . . . . . . . G2 Variable conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Observer model variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Controller model variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 System covariance values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Measurement covariance values . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Coefficients for φ(ω e , T e ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Classic driving style vs. optimal control . . . . . . . . . . . . . . . . . . . . . C.1 Node entry data format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Way entry data format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.3 Quad-tile bounding box data format . . . . . . . . . . . . . . . . . . . . . . . C.4 Nodes binary search tree data format . . . . . . . . . . . . . . . . . . . . . . . C.5 Ways binary search tree data format . . . . . . . . . . . . . . . . . . . . . . . C.6 Quad-tile bounding box quad search tree data format . . . . . . . . . . . . . xiii List of Abbreviations

Table 1 :

 1 Notation conventions

	Description	Example Meaning
	Bold, lowercase letter	x	vector
	Bold, upppercase letter A	matrix
	c, subscripted	A c	continuous time
	d, subscripted	A d	discrete time
	i, j, subscripted	A i,j	i th , j th element
	k, subscripted	x k	k th time step
	T , superscripted	A T	matrix transpose
	Hat, over letter	x	estimated value (as opposed to real)
	Star, superscripted	u *	optimal
	Minus, superscripted	P -k	intermediate value between time k and time k + 1
	Dot, superscripted	ẋ	time deivative
	Variable conventions
	Throughout this dissertation, the following variable conventions apply:
	Variable Meaning
	f	system model vector valued function
	h	measurement model vector valued function
	x	state vector
	y	measurement vector
	u	control vector
	A	Jacobian of f with respect to x
	B	Jacobian of f with respect to u
	H	Jacobian of h with respect to x
	J	cost
	K	gain matrix
	P	state noise covariance
	n x	number of state variables
	n y	number of measurement variables
	n u	number of control variables
	∆t	time step
	η	engine efficiency
	φ	inverse of engine efficiency, i.e. 1/η

Table 2 :

 2 Variable conventions Complete observer model . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Observability normal form . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Controller model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 1
	State of the (Theoretical) Art

Contents 1.1 State models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Observer model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1.1 Model variables . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1.2 Reference frames . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1.3 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1.4 Energy analysis . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 1.3.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1.1 Model variables . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1.2 Cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1.3 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Complete controller model . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 .

 1 1.

	Variable	Description	Units
	x i	car displacement in the î axis	[m]
	v i a i , a j , a k F thrust	car speed in the î axis car acceleration in the ijk axes thrust force along î axis	[ m s ] [ m s 2 ] [N]
	ψ θ S j	car heading ∈ [0, 2π) car pitch ∈ [-π 2 , π 2 ) road straightness coefficient	[rad] [rad] [ rad m ]
		(horizontal plane)
	S k	road straightness coefficient	[ rad m ]
		(vertical plane)
	N, E, D	Northing, Easting, and Down	[m]
		(local projection plane)
	E in	energy content input, fuel	[J]
	P in	power input, fuel	[W]
	b i , b j , b k accelerometer biases in the ijk axes [ m s 2 ] θ e engine angular displacement [rad]
	ω e C rr	engine angular speed rolling resistance coeffecient	[ rad s ] [-]
	C d	air drag coefficient	[-]
	A	car frontal cross-section area	[m 2 ]
	ρ M	air density total vehicle mass	[ kg m 3 ] [kg]
	M a g	total accelerated mass acceleration due to gravity	[kg] [ m s 2 ]

Table 1 .

 1

1: Observer model variables 1.2.1.2 Reference frames

  Microprocessors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Microcontrollers . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 2
	State of the (Practical) Art
	Contents
	2.1 . . . . . . . . . . . . . . . . . . . . . . .
	2.5.2 Digital Elevation Map . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3 Sensor package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 GPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Accelerometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Odometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Fuel injector timer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4.1 Engine angular distance/velocity . . . . . . . . . . . . . . . 2.3.4.2 Fuel flow meter . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Other sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Gyroscope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Compass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Additional equipment . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Data storage . . . . . . . .

  Polynomials of highest degree 4 . . . . . . . . . . . . . . . 3.4.3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Additional results . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Contents 3.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Values used . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Car model observer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2.1 Car model observability . . . . . . . . . . . . . . . . . . . . 3.1.2.2 Choice of observer . . . . . . . . . . . . . . . . . . . . . . . 3.1.2.3 Jacobian . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Observer pseudo-code . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4.1 Tuning state model covariances . . . . . . . . . . . . . . . . 3.1.4.2 Tuning measurement model covariances . . . . . . . . . . . 3.1.4.3 Values used . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4.4 Normalization . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Data gathering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Highway . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 In-town . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Maximum acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Efficiency map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Curve fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Efficiency fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1.1 Polynomials of highest degree 2 . . . . . . . . . . . . . . . 3.4.1.2 Polynomials of highest degree 3 . . . . . . . . . . . . . . . 3.4.1.3 Polynomials of highest degree 4 . . . . . . . . . . . . . . . 3.4.1.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Specific fuel consumption fitting . . . . . . . . . . . . . . . . . . . . 3.4.2.1 Polynomials of highest degree 2 . . . . . . . . . . . . . . . 3.4.2.2 Polynomials of highest degree 3 . . . . . . . . . . . . . . . 3.4.2.3 Polynomials of highest degree 4 . . . . . . . . . . . . . . . 3.4.2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 P in fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3.1 Polynomials of highest degree 2 . . . . . . . . . . . . . . . 3.4.3.2 Polynomials of highest degree 3 . . . . . . . . . . . . . . . 3.4.3.3

Table 3 .

 3 1 gives typical model parameter values. The car parameters are taken from published values for the 2004 Smart Roadster Coupé. Note that mass is considered to include an average 70kg driver with 5kg of clothing and personal material.

	Parameter Value
	g	9.805
	ρ	1.2
	C rr	0.010
	C d	0.39
	M	920
	M a	948
	A	1.93
	r w	0.2954

Table 3 .

 3

1: Model parameters 3.1.2 Car model observer 3.1.2.1 Car model observability

Table 3 .

 3 2: System covariance values

	Variable	Value

Table 4 .

 4 .1. 1: Coefficients for φ = α 0 + α 1 ω e + α 2 T e + α 3 ω e T e + α 4 ω 2 e + α 5 T 2

	Coefficient Value
	α 0	9.4360
	α 1	-0.018136
	α 2	-0.068797
	α 3	3.7201e-06
	α 4	2.1965e-05
	α 5	3.6204e-04
		e

Table C .

 C Each way ID is represented by a 32-bit unsigned integer. The number of nodes that a way includes is given (in the case of a way that has more than 256 nodes, it can be split into two ways with no loss of information), followed by the node IDs themselves. Lastly, a way is split into a given number of subpoints, which do not necessarily correspond to the nodes making up the way, and the elevation for each node is given.

	Node entry Latitude 32-bit signed of nodes Node ID 32-bit unsigned Way ID # n Nodes Longitude # of Intersecting ways n Intersecting ways ID 32-bit signed 8-bit unsigned n * 32-bit unsigned Table C.1: Node entry data format # of Elevation points n Elevation points unsigned unsigned unsigned Way entry: Way entry 32-bit 8-bit n * 32-bit 8-bit unsigned n * 32-bit unsigned

2: Way entry data format

Note 30 If it is not immediately obvious why simply using nodes for elevation is unacceptable, consider that most mapping agencies only use nodes when necessary to represent the curvature of a road. The road is considered to be straight between two nodes. If the distance between two nodes is several kilometers (imagine a road in the Australian Outback, or the desert in western America), there will certainly be a number of hills that are hidden by this oversimplification.

Table C .

 C Note 31There is nor preferred sense to a quad-tile's binary representation. That is to say, a level 1 tile could as well put level one completely at the right, and the left-shift as the zoom increases, or a level 1 zoom could be associated only with the first two left-most bits (otherwise known as MSB) followed by level 2 at the third and fourth bits, etc... The choice depends on the needs of a particular system.Quad-tile bounding box search algorithm: A somewhat more complex quad search tree. The way value data points to the appropriate information in the data table. The center, i.e. starting node, is the very first file entry in the index file. All searches are based on this starting point. In a quadtree, subgroups are divided into four leaves (See Fig.C.3), not two as in a binary search tree. As in the case for quad-tiles, 00 represents the NW quadrant, 01 the SW, 10 the NE, and 11 the SE one.

	Quad tile quad search tree
	Quad-	Pointer to	Pointer to	Pointer to	Pointer to	Pointer to en-
	tile	00 tile	01 tile	10 tile	11 tile		try in data file
	ID
	64-bit	32-bit	32-bit	32-bit	32-bit		64-bit
	unsigned	unsigned	unsigned	unsigned	unsigned	unsigned
	Quad-tile bounding box entry
	Quad tile coordinate ID	Zoom level	# of contained ways n Contained ways
	64-bit unsigned	8-bit unsigned	8-bit unsigned	n * 32-bit unsigned
		Table C.3: Quad-tile bounding box data format

6: Quad-tile bounding box quad search tree data format

Literally, "shortest time". The Brachistochrone Problem was one of the earliest problems posed in the calculus of variations. The problem was to "find the shape of the curve down which a bead sliding from rest and accelerated by gravity will slip (without friction) from one point to another in the least time". Newton solved it in one day. The answer: a cycloid.

This version exists today, but testing has not yet been completed.

D is a nuanced case. Without the GPS, it is impossible to fix the car's position sufficiently well on the map in order to be able to use a DEM.

The reason for this might not be immediately obvious to the lay eye. Efficiency is a function of many variables, but one important variable is pumping losses. As the throttle plate restricts the airflow into the engine, it creates a zone of low pressure behind the throttle body. The engine must work harder to create this vacuum, and so there is more energy lost due to this pressure differential.

A too high standard deviation means a too large variation of ωe, Te in the grid point, which means that it might not accurately reflect reality. As there are already many valid grid points, it does not harm the overall fitting to discard these points.

While the measurement is expressed in archaic units, it is a superior way of expressing fuel economy to its pseudo-inverse "L per 100km". Due to the human propensity to understand, and favor, large numbers over small, it is sociologically better to link a greater good with a greater quantity. The difference between

2.9L/100km and 3.0L/100km appears miniscule to a potential buyer, but it represents years of research and refinement. Whereas the difference between 78mpg and 81mpg is clear to even the most mathematically challenged among us. The author hopes to one day see car efficiency given in kpl.

However, the sensor suite and observer developed earlier in the dissertation are already being used by other researchers in order to gather data for modeling traffic flow.

In the case of ωe, xi, vi, this is immediately obvious, as they can only increase. Te, on the other hand, is more subtle. If Te were negative, this would imply a braking phase. In a pure acceleration phase, clearly braking would be less efficient, hence Te ≥ 0

Manufacturer's data indicates 0-100 in 10.9s, but this is not a straight-up apples-to-apples comparison. In the maximum acceleration test, the real Smart loses a significant amount of time in gear shifts (≈ 3 s), but makes up by accelerating harder between shifts.

This representation assumes a Gaussian noise distribution. This is not strictly true, and in fact is likely not true for most models. However, the Gaussian assumption leads to mathematically proven observers that have proven to be robust even when the Gaussian assumption is patently false. The most general model would be f (x, u, w) = Φ, but this is so abstract as to lose applicability to many real-world problems.

The resolvent is the integral of Ac along the path x(t). It has the interesting property that Φ t k+1 /t k Φ t k /t k-1 = Φ t k+1 /t k-1 and so x k+1 = Φ t k+1 /t 0 x(t0)

Indeed, it is its principle disadvantage as the Jacobian is a first order approximation of a non-linear system at a specific point, a point which is only an estimation of the real state, and so by its very nature is suboptimal, and can affect the filter accuracy or even lead to divergence of the filter. Many other filters have been developed that seek to address this shortcoming. See Section A.3.2 for an example of a popular one.

Instead of using a true analog computer for the correction, a daunting task to say the least, digital computers running at extremely high sample rates are used instead.

Of course, this goes without saying as it is a necessary result for a linear system. The approach is presented for completeness in the case of a non-linear PMP optimization.

Note that B → 0 is not treated as this would imply u * f → ∞. Clearly if the control can be infinite, then in the case of a linear system T → 0.

Acknowledgments

acceleration profile to the non-fixed horizon case, we test only the latter.

A convenient section of road was chosen as the test route. The 2km length was relatively flat and straight, and allowed unimpeded acceleration to 90kph. While not ideal, as there were cars parked on either side of the road, we were able to perform conclusive tests that show that Pontryagin's Maximum Principle does indeed lead to a lower fuel consumption. Results can be seen in Figs. 4.8 and 4.9 and are summed up in Tab. 4

A.3.3 Square root filters

An algorithmic variant of classic filters is created by taking the square root of the state noise covariance matrix, P. In these cases, the mathematical characteristics are preserved, for instance positive-definiteness is guaranteed, and the numerical properties are improved, as computation is done more accurately because very small numbers become larger, and very large numbers smaller. [START_REF] Chui | Kalman Filtering with Real-Time Applications[END_REF] A.3.3.1 Square-root continuous-discrete EKF

The following algorithm is detailed in [4,[START_REF] Chui | Kalman Filtering with Real-Time Applications[END_REF], but it is repeated below for ease of reference. Assume x 0 , P 0 given. (If P 0 is not given, it can generally be initialized with the identity matrix.) Let

Define the following Cholesky square roots:

Prediction: Form matrix L 1 as:

and calculate M 1 by taking the Q factor of

where l is not further used. Determine S - k by taking the first n x columns of M 1

Correction: Compute the measurement model Jacobian H k = H c (x(t k )) as in eq. (A.11). Form matrix L 2 as:

and calculate M 2 by taking the Q factor of

where l is not further used.

Extract the matrices

Called the Square Root Unscented Kalman-Bucy Filter, it provides results mathematically identical to the Unscented Kalman-Bucy Filter, only for reasons of numerical stability and speed it avoids recalculating the cholesky decomposition at each time step.

The following algorithm is detailed in [START_REF] Sarkka | On unscented kalman filtering for state estimation of continuous-time nonlinear systems[END_REF][START_REF] Van Der Merwe | The square-root unscented kalman filter for state and parameter estimation[END_REF], but it is repeated below for ease of reference.

Prediction:

where

31)

Correction:

A.4 Adaptive High-gain EKF

The AEKF is a recent innovation in the field of non-linear observers. Proposed in [START_REF] Boizot | An adaptive high-gain observer for nonlinear systems[END_REF] as a solution that combines the best features of EKF noise filtering and high-gain EKF perturbation sensitivity, the AEKF is based on the concept of innovation, a metric that indicates the deviation between the historical predicted state and the historical measured state. When the innovation becomes arbitrarily high, it triggers an increase in gain that switches the AEKF from unity-gain mode to high-gain mode.

Furthermore, the AEKF is a mathematically proved observer. Its stability is grounded by theory, as opposed to the simple EKF. This can have important consequences for safety, if an observer is used for closed-loop control purposes.

A.4.1 Observability normal form

(In the case of a multiple-input/multiple-output system, a change of variables is not necessarily unique and thus we talk of normal form instead of canonical form.)

High-gain observers, developed within the framework of [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF], have been proven to achieve global exponential convergence. (Recall that the EKF has only been proven to converge locally, i.e. when the estimated state is in a neighborhood ε of the true state.) This exponential convergence property requires that the model be in a very specific form, called observability normal form, which characterizes the observability of the system. The general multiple-output normal form is described in Sec. A.4.3.

Note 20

At risk of stating the obvious, before building any kind of observer it is necessary to first verify that the system is observable. Obviously, it is impossible to design an observer if the system itself is not observable (no matter the change of coordinates, observability is an intrinsic property). However, all methods by which the model is brought into the normal form are by definition methods that show observability. Therefore, we can verify observability and compute the change of coordinates in one single action.

A.4.2 Continuous-discrete AEKF

The multiple output CD-AEKF [START_REF] Boizot | Adaptive High-gain Extended Kalman Filter and Applications[END_REF] is repeated here for completeness. We write the continuous-discrete time normal form as

For a system with n x states and n y measurements:

where P| t k-1 = Pk-1 and ẑ| t k-1 = ẑk-1

Note 22

The above fusion techniques are equally valid for the discrete-discrete approach. However, if the DD-EKF had been used instead of the CD-EKF, Q would have varied at each iteration as a function of time since the covariance of a gaussian noise of covariance q after a time ∆t is q∆t. I.e., the longer the observer waits without any measurement data, the higher the uncertainty, which is represented by a growing q. Thus in the discrete-discrete case: This file can be adapted to any new car. For the most part, it is self explanatory. However, included are explanations for the few that are not:

• engine-> inertia lumped: the total lumped rotational inertia for all of the drivetrain except the tires and wheels.

• accelerometer-> calibration: the rotation matrix and scale that must be applied to the accelerations measured by the sensor in order to bring them into the body frame.

• drag-> Cd: Coefficient of air resistance

• drag-> Crr: Coefficient of rolling resistance

• area-> frontal: maximum cross-sectional area as seen from the front

• mass-> wet: car mass including all fluids (coolant, oil, windshield wiper fluid, etc...)

• inj-> calibration -> lag: the difference between the actual injector opening time and the injector voltage pulse length, in [s]

• inj-> calibration -> fuelPerSecond: the quantity of fuel delivered by all injectors in one second

• vss-> calibration: how far the car rolls forward per vss pulse

• file Name: the individual trip data files that have been imported for a given car.

C.1.3 Observer file generation overview

When data is imported, a certain number of files are generated:

• <file name>, data.mat Data extracted from all sensors, except GPS

• <file name>, extract.mat Data extracted from the GPS NMEA strings

• <file name>, results.mat The Kalman Filter results

• <file name>, DEM.mat Optional. An ordered list of the ways taken during a trip.

• <file name>, gps.txt All the GPS NMEA strings in text format

• <file name>, spd.kmz A Google Earth file in KML format that shows the trip path and speed. See App. F for details. Here we discuss the data format used by the sensor suite for saving data to the SD card. Logged data is saved sequentially onto the SD card, with all data sources interwoven together.

Each type of data is prefaced by a preamble and ended with an error checking byte (modified CRC). (See 4. The message itself depends on the data type. See the following section, C.1.5.

5.

Finally, the message is terminated by an error-checking byte. This byte is calculated as the sum modulo 256 of a subset of the message bytes, starting with the data length byte and continuing to the last message byte (therefore not including the error-checking byte itself).

Note 32 Classically, CRC algorithms use XOR operations, instead of addition modulo 256. Future versions of the logging program will move to the true CRC algorithm.

C.1.5 Data types

C.1.5.1 Error messages

Error messages are saved in ASCII human-readable format, followed by a newline.

E.2 Circuit-board layout

The circuit board was designed in Eagle. Top layer: The Google Earth Toolbox [START_REF]Google Earth Toolbox[END_REF], is a useful collection of utilities that aid in the presentation, understanding, and analysis of geodetic data. Certainly, this project could not have gone as far as it did if it were not for the excellent capabilities of Google Earth and Google Earth Toolbox.

For this project, Google Earth Toolbox code was modified: cleaned, streamlined, and upgraded. Specific needs that were not present in Google Earth Toolbox were addressed, such as the (in)ability to display properly formatted XML code, or to output arbitrary XML code to the Google Earth KML format.

The Google Earth files generated by the Matlab code allow one to see, at a glance, car path, car speed, starting and ending points, and trip time. Take for instance Fig. • Speed, as represented by color. The colormap was intentionally chosen to imitate Matlab's "jet". Dark blue represents a stopped car, whereas deep red represents high speed. All outputs are normalized to 180km/hr, a value that seems quite reasonable for the upper end, and yet still leaves room for definition at slower speeds

• Waypoints, as represented by the blue boxes with floating numbers. This indicates the trip time that the car was at the particular spot.

Trip time is perhaps the most important feature of the Google Toolbox Extensions, as it allows easy cross correlation between logged data time (absolute ms since the trip started) and GPS data time (UTC).

Code modifications have been submitted to Google Earth Toolbox for integration into future versions of the Toolbox.

Appendix G

The Ball and Wheel: a worked example for Pontryagin The PhD presents a combined approach to improving individual car efficiency. An optimal observer, the Extended Kalman Filter, is used to create an efficiency model for the car. Particular attention was paid to handling the asynchronous and redundant nature of the measurement data.

A low-cost sensor suite developed to measure data is described. This sensor suite was installed on multiple vehicles to good success. It employs an accelerometer, gps, fuel injector timer, and Vss input to measure all the data necessary to reconstruct the car's state. This observer and sensor suite can be used as the base for any study which requires car efficiency maps, allowing research to proceed without manufacturer supplied data.

Once the efficiency map is found, it is then curve-fitted in order to reduce model complexity. The simplified model is then used as a basis for optimal control through Pontryagin's Maximum Principle.

Real-world test results are given, both for efficiency mapping, and for optimal control. Detailed discussion of the observer and controller is presented, in order to ease understanding and save implementation time.