Bayesian Filtering and Smoothing to Measure Damper Characteristics

(1)

DEGREE PROJECT, IN MECHATRONICS , SECOND LEVEL

STOCKHOLM, SWEDEN 2014

Bayesian Filtering and Smoothing to

Measure Damper Characteristics

CHRISTIAN SVENSSON

(2)

(3)

Bayesian Filtering and Smoothing to Measure

Damper Characteristics

Christian Svensson

Master of Science Thesis MMK 2014:57 MDA 488 KTH Industrial Engineering and Management

(4)

(5)

Master of Science Thesis MMK 2014:57 MDA 488 Bayesian Filtering and Smoothing to Measure

Damper Characteristics Christian Svensson

Approved Examiner Supervisor

2014-06-30 Jan Wikander Bengt Eriksson

Commissioner Contact person

Öhlins Racing AB Janne Färm

Abstract

During the development of high performance damper products for vehicle ap-plications, it is essential to measure the dampers’ characteristics to ensure their performance. Some of the measurements are performed in dynamometers. The dynamometer actuates the damper with a predefined position signal and measures the actual position and the resulting force. These measurements su↵er from noise; the problem is especially bad if the signals are studied with respect to velocity. This is due to that the velocity is not directly measurable; it must be calculated from position by di↵erentiation, which amplifies the noise.

The purpose of this work was to improve these measurements by doing a comparison of di↵erent Bayesian filters and smoothers. The comparison inclu-ded the Kalman filter, the exteninclu-ded Kalman filter, the unscented Kalman filter, the Rauch-Tung-Striebel smoother, the extended Rauch-Tung-Striebel smoother and the unscented Rauch-Tung-Striebel smoother. The smoothers are only applicable in o✏ine applications, while the filters could be used in real time.

(6)

(7)

Examensarbete MMK 2014:57 MDA 488 Bayesisk filtrering och glättning för mätning av

stötdämparkaraktäristik Christian Svensson

Godkänt Examinator Handledare

2014-06-30 Jan Wikander Bengt Eriksson

Uppdragsgivare Kontaktperson

Öhlins Racing AB Janne Färm

Sammanfattning

Under utvecklingen av högpesterande dämparprodukter för fordonstillämpningar är det ett måste att mäta dämparnas karakteristik för att säkerställa deras prestanda. Vissa av mätningarna görs i dynamometrar. Dynamometern aktuerar dämparen med en fördefinierad positionssignal och mäter den faktiska positionen samt den resulterande kraften. Dessa mätningar påverkas av brus, problemet är extra tydligt om signalerna studeras med avseende på hastighet. Detta beror på att hastigheten inte är direkt mätbar, utan måste räknas ut genom derivering, vilken förstärker bruset.

Syftet med detta arbete var att förbättra dessa mätningar genom att göra en jämförelse mellan olika Bayesiska filter och glättare. I jämförelsen ingick Kal-manfiltret, det utökade KalKal-manfiltret, det oparfymerade KalKal-manfiltret, Rauch-Tung-Stribelglättaren, den utökade Rauch-Tung-Striebelglättaren och den oparfymera-de Rauch-Tung-Stribelglättaren. Glättarna kan endast användas o✏ine, medan fil-terna kan användas i realtid.

(8)

(9)

Chapter 1 Introduction

This chapter explains the background of why this work was needed and its pur-pose. The chapter also gives an introduction to Bayesian filtering and smoothing.

1.1 Background

This work was done in collaboration with the company Öhlins Racing AB in Upplands Väsby, Sweden. Öhlins Racing AB develops and manufactures high performance damping products for vehicle applications. The products are found in vehicles such as super sport cars, sport bikes, snowmobiles, downhill bicycles and various types of racing vehicles.

During the development of high performance damper products, it is crucial to test the products to ensure that they meet the engineers’ and the customers’ re-quirements. Some of the tests are performed with purpose-built dynamometers, which actuate the damper with a predefined position signal and measures position and resulting force. In some cases, the accuracy of the measurements is not sat-isfactory. This is due to measurement noise and quantities that are not directly measurable, e.g. velocity, which today is calculated from the position signal.

(12)

CHAPTER 1. INTRODUCTION −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 −2000 −1500 −1000 −500 0 500 1000 1500 2000 2500 3000 Velocity/(m/s) Force/N Low damping Medium damping High damping

Figure 1.1: Example damping force curves from dynamometer measurements. These curves consist of mean valued data points at constant velocities. This measurement method results in a quasi-stationary relationship between velocity and force. This is also the damping curves for the three di↵erent settings that were used in this work.

1.1.1 The Dynamometer

The dynamometer could be divided into four sub-systems: Actuator, mounting assembly for the test object, sensors and data acquisition system. The actuator in this case is a hydraulic cylinder, which is attached to the mounting assembly. The actuator is regulated with a closed loop controller, usually with position feedback. Today, primary two sensors are in use; one load cell and one linear variable dif-ferential transformer (LVDT) position sensor. The load cell is placed at the top of the mounting assembly and its purpose is to measure force. The LVDT position sensor is found inside of the hydraulic cylinder and is measuring cylinder posi-tion. The data acquisition system consists of hardware to interface the sensors and a computer to log the measurements. See Figure 1.2 for a schematic illustration of the dynamometer.

(13)

1.2. PURPOSE AND DELIMITATION

Figure 1.2: A schematic illustration of the dynamometer.

1.2 Purpose and Delimitation

The purpose of this work was to improve the dynamometer measurements by comparing di↵erent Bayesian filters and smoothers.

This work was delimited to Bayesian filtering and smoothing with focus on techniques suitable for unimodal distributed noise. This delimitation is motivated in section 2.1.

1.3 Bayesian Filtering and Smoothing

(14)

CHAPTER 1. INTRODUCTION

Figure 1.3: An illustration of the filtering process.

how noisy the prediction and the measurements are. A noisy prediction or mea-surement is trusted less than a more accurate one. The principle behind Bayesian filtering is illustrated in Figure 1.3.

One of the most famous implementation of the Bayes filter algorithm is the Kalman filter (KF) (Kalman, 1960). The KF assumes that the underlying system is linear and that model and sensor noise are Gaussian distributed with zero mean, i.e. white Gaussian noise. There exist several solutions to extend the KF for nonlinear systems. One common method is the extended Kalman filter (EKF) (see e.g. Thrun et al., 2005; Anderson & Moore, 1979), which uses a first order Taylor expansion to linearize the underlying system. Another approach presented by Julier & Uhlmann (1997) linearize the system with so called sigma points. This method is called the unscented Kalman filter (UKF) and is more accurate than the EKF and have the same order of computational complexity.

The filters introduced above treats current data and data collected in the past. However, if all the data is already collected, i.e. the application is not a real time problem, then it is possible to get better results by considering all measurements. One way of doing this is to run the filter from the oldest measurement to the newest and then run the filter backwards. This results in a much better estimate and is called smoothing. The di↵erences of prediction, filtering and smoothing are illustrated in Figure 1.4. All the filters mentioned above have a corresponding smoother.

The KF and the EKF assumes that the noise is white Gaussian, the UKF can compensate for some deviations from the standard Gaussian distribution depend-ing on how the sigma points are chosen. All three methods are limited to unimodal distributions, which are assumed in this case. An analysis of the measurement noise is given in section 3.2.

(15)

1.3. BAYESIAN FILTERING AND SMOOTHING

(16)

(17)

Chapter 2 Frame of Reference

This chapter gives a theoretical description of the methods and theories that form the foundation of this work. The chapter starts with a description of the models that are considered. Then each filter method is presented, followed by its corre-sponding smoother. The methods are presented with relevant assumptions in mind and not as generalized as possible. E.g. only time invariant models are considered, even though all filters and smoothers can handle time variant systems.

2.1 A Motivation to the Chosen Methods

The three chosen filtering methods and their corresponding smoothers were cho-sen since they probably are the most common solutions for similar filtering prob-lems. However, there are other methods in the literature that deserves to be men-tioned. A method fairly similar to the UKF is the central di↵erence Kalman filter (CDKF), which was developed by Ito & Xiong (2000) and Nørgaard et al. (2000). According to Van Der Merwe (2004), this filter’s performance is almost identical to the UKF. For this reason and due to that it seems less popular, it is not included in this work.

Another method is the Gauss-Hermite Kalman filter (Ito & Xiong, 2000). Ac-cording to Särrkä (2013, p. 205) this filter seems to be more robust and consistent in errors compared to the UKF. But it is su↵ering from high computational com-plexity and it does not seem to introduce any significant improvements over the UKF, therefore it is not a part of this comparison.

(18)

CHAPTER 2. FRAME OF REFERENCE An option that is completely di↵erent from the Kalman filtering methods is the particle filter (PF). The PF was introduced by Gordon et al. (1993) and represents the state space with a large number of samples, called particles. The PF is not limited to unimodal or Gaussian posterior distributions; this makes it popular for localization methods and problems with posteriors that cannot be approximated as a Gaussian with good accuracy. However, the nature of the problem in this work makes it suitable for Kalman filtering methods. The posterior is assumed to be approximately Gaussian and therefore, the PF is not included in the comparison. A good introduction to particle filters and smoothing could be found in a tutorial by Doucet & Johansen (2009).

2.2 Basic Probability Theory

The probability that a random variable X takes the value x is denoted p(X = x) .

A common abbreviation is to omit the random variable as p(x) = p(X = x) .

The probability of that X = x and the random variable Y = y is p(X = x and Y = y) = p(x, y) .

This is called the joint distribution between X and Y. Another important concept is conditional probability. A conditional probability is if an event Y = y happens given that X = x already has happened. This is written as

p_{(y | x) .}

2.3 Model Basics

Consider the following system of nonlinear first order di↵erential equations: ˙x(t) = gc x(t),u(t) , (2.1a)

y(t) = hc x(t) , (2.1b)

(19)

2.3. MODEL BASICS

signals, y(t) 2 Rp _{is a vector with the measured quantities and} _h

c(·) 2 Rp is a nonlinear, time invariant vector valued function called the measurement model. The subscript "c" indicate continuous time.

A special case of Equation 2.1 is if the functions g_c(·) and hc(·) are linear, then the system could be described as

˙x(t) = Acx(t) + Bcu(t), (2.2a)

y(t) = Cx(t), (2.2b)

whereAc 2 Rn⇥nis a matrix describing the system behavior,Bc 2 Rn⇥m is a ma-trix describing the control signals’ impact on the system andC 2 Rp⇥n_describes which states that are observed, with proper scaling.

If the input u(t) is constant during equidistantly spaced time intervals, with the length Ts, such that

u (k + ⌧)Ts =u(kTs), 0  ⌧ < 1, k = 1,2,. . . , (2.3) then the time continuous system in Equation 2.1 could be written in discrete time as

x(kTs) =g⇣x (k 1)Ts ,u (k 1)Ts ⌘, (2.4a)

y(kTs) =h x(kTs) , (2.4b)

whereg(·) and h(·) are the discrete counterparts of gc(·) and hc(·) (Ljung & Glad, 2004). These functions could be approximated with e.g. Euler forward, Euler backward or Tustins method. The linear state space in Equation 2.2 is discretized with zero-order hold as

x(kTs) =Ax (k 1)Ts +Bu (k 1)Ts , (2.5a)

y(kTs) =Cx(kTs), (2.5b)

where the matricesA and B are given by A B 0 I ! = e 0 BB@A₀c B₀c 1 CCATs , (2.6)

whereI is the identity matrix (DeCarlo, 1989, p. 215). For practical reasons, the following notation will be used:

xk :=x(kTs), with this notation Equation 2.5 becomes

xk = Axk 1+Buk 1, (2.7a)

(20)

CHAPTER 2. FRAME OF REFERENCE

Figure 2.1: A visualization of a hidden Markov model. A state only depends on the previous state and the control that form the transition between the states. The states, which are hidden, are only observable through noisy measurements.

In reality, the models are seldom perfect and sensor readings are a↵ected by noise. These imperfections could be modeled as additive noise with

xk =Axk 1+Buk 1+ k, (2.8a)

yk =Cxk+"k, (2.8b)

where k 2 Rnis the process noise and "k 2 Rpis the measurement noise. Equation 2.8 could also be written in probabilistic form as

xk ⇠ p(xk| xk 1,uk 1) , (2.9a)

y_k ⇠ p yk| xk . (2.9b)

As seen in Equation 2.9a, the statexkonly depends on the previous statexk 1 and controluk 1. This mean that all states beforexk 1 are irrelevant for the current state. The true state is also unknown, it could only be observed through noisy measurements. This is called a hidden Markov model (HMM) and is illustrated in Figure 2.1.

2.4 Derivation of Models

This section describes the di↵erent models that have been used. A motivation to the exclusion of models for the mechanical parts of the dynamometer is given in Appendix A.

2.4.1 Damping Force

The simplest model of a damper is a linear relationship between the damping force Fd and velocity ˙z:

Fd =d ˙z, (2.10)

(21)

2.4. DERIVATION OF MODELS

where d is the constant damping coefficient. This model is suitable and used for the Kalman filter. However, most dampers are better described with the nonlinear model

Fd = d( ˙z) ˙z. (2.11)

The nonlinear damping coefficient d( ˙z) could preferably be implemented from measurements as a piecewise linear function.

2.4.2 Spring Force

A pressurized non-through rod damper will have a position dependent force due to the displacement of the piston rod. This position dependent force could be described as

Fk = kz, (2.12)

where k is the spring constant. A detailed analysis of the gas spring force is found in Appendix B.

2.4.3 Static Friction

The static friction force Ff is modeled as

Ff = Fcsgn( ˙z), (2.13)

where Fcis the absolute friction force and sgn(·) is the signum function. In com-puter simulations, this is commonly implemented as Karnopp’s stick-slip model (Karnopp, 1985). However, if the implementation mentioned in subsection 2.4.1 with a piecewise linear damping force model, the static friction is already incor-porated from the measurements.

2.4.4 Total Damper Force

The total damper force is given by the sum of the forces from the modeled damper components:

F = Fd+Fk +Ff. (2.14)

2.4.5 Linear Model for the Kalman Filter

(22)

CHAPTER 2. FRAME OF REFERENCE to estimate the position z, velocity ˙z and force F, these will be the states. This gives the state vector

x(t) = 0 BB@xx12(t)(t) x₃(t) 1 CCA = 0 BB@z(t)˙z(t) F(t) 1 CCA . (2.15)

Now, the matrices in Equation 2.2 needs to be derived. The time derivative of the first state is

˙x1(t) = x2(t). (2.16)

The time derivative for the second state is acceleration, which is given as an input. This gives

˙x2(t) = u(t). (2.17)

The equation for the third and last state is given by taking the time derivative of Equation 2.14. Since this is a linear model, the nonlinear static friction is neglected and the linear damping coefficient is used. This gives

˙F(t) = k ˙z(t) + d ¨z(t) = kx2(t) + du(t). (2.18) Finally, the model could be written as

˙ x(t) = 0 BB@0 1 00 0 0 0 k 0 1 CCA | {z } Ac x(t) + 0 BB@01 d 1 CCA |{z} Bc u(t), (2.19) y(t) = 1 0 0_{0 0 1} ! | {z } C x(t). (2.20)

An illustration describing the system model can be seen in Figure 2.2.

(23)

2.5. THE KALMAN FILTER

Figure 2.2: A figure describing the system model with sensors and actuator.

2.5 The Kalman Filter

The KF was derived by Kalman (1960) with orthogonal projections as a recursive solution to the linear filtering problem. The KF is an optimal estimator if the system is linear and the noises are white Gaussians (Anderson & Moore, 1979, pp. 46–49), i.e. on the form in Equation 2.8 with

k ⇠ N (0,Q) (2.21a)

"_k _{⇠ N (0,R)} (2.21b)

where N (·) is the multivariate normal distribution function with the mean as first argument and the covariance matrix as second. Q 2 Rn⇥n _and_{R 2 R}p⇥p _{are the} covariance matrices for the process and the measurement noise, respectively. Note that the noises are assumed to be stationary, i.e. the statistics does not change over time. This is not an assumption claimed by the KF, but an assumption made in this work.

Since the dynamics are linear and the disturbances are Gaussians, the posterior distribution will be Gaussian as well. This is formulated as

p xk| y1:k, u1:k =N µk, ⌃k (2.22) where µ_k is the mean vector and ⌃k is the covariance matrix for the belief. y1:k are all measurements from time step 1 to k andu1:kare all controls from time step 1 to k.

(24)

CHAPTER 2. FRAME OF REFERENCE the prior distribution, i.e. before the measurement is taken into account. This is denoted by the minus sign superscript.

The next and last step is called the update step and is performed in the lines 4 to 6. In this step the measurement is taken into account and statistics for the posterior distribution are calculated. The matrix calculated in line 4,Kk, is called the Kalman gain and determines how much the filter should trust the measurement or the process.

Algorithm 1 Kalman Filter

1: function KF(µ_{k 1}, ⌃k 1,uk 1,yk) 2: µ_k = Aµ_{k 1}+Bu_{k 1} 3: ⌃_k =A⌃_{k 1}AT+Q 4: K_k = ⌃_kCT ⇣C⌃_kCT+R⌘ 1 5: µ_k = µ_k +K_k ⇣y_k Cµ_k⌘ 6: ⌃_k = (I K_kC) ⌃_k 7: return µ_k, ⌃k 8: end function

2.6 The Rauch-Tung-Striebel Smoother

The KF only considers measurements before and at time step k, which is well suited for online implementations. However, if all data already is collected and to be treated o✏ine, it is possible to use the information from all measurements up to time step N. The belief for this problem is

p xk| y1:N,u1:N =N⇣µs_k, ⌃s_k⌘, 1  k < N, (2.23) where µs

k is the smoothed mean and ⌃sk is the smoothed covariance matrix. A solution to this problem is the Rauch-Tung-Striebel smoother (RTSS), which was presented by Rauch et al. (1965).

This algorithm requires the means and covariances at all time steps from the Kalman filter. So first the data is Kalman filtered and then the smoothing algorithm is started with µs

N = µN and ⌃sN = ⌃N and goes backwards in time. The algorithm is similar to the KF and is presented in Algorithm 2.

(25)

2.7. THE EXTENDED KALMAN FILTER Algorithm 2 Rauch-Tung-Striebel Smoother

1: function RTSS(µs_k+1, ⌃s_k+1,µ_k, ⌃k,uk) 2: µ_k+1=Aµ_k+Bu_k 3: ⌃_k+1 =A⌃kAT+Q 4: L_k = ⌃_kAT⇣⌃_k+1⌘ 1 5: µs_k = µ_k+L_k ⇣µs_k+1 µ_k+1⌘ 6: ⌃s_k =⌃_k +L_k ⇣⌃s_k+1 ⌃_k+1⌘LT_k 7: return µs_k, ⌃s_k 8: end function

2.7 The Extended Kalman Filter

As seen in section 2.5, the KF requires linear models. However, if a system is nonlinear, it could be approximated as a linear system with a first order truncated Taylor expansion. This yields the EKF and is probably one of the most used methods for filtering problems.

The EKF is presented in Algorithm 3. The predict step takes place in line 2 to 4 and the update is done in line 4 to 8. Since the Jacobians of the models are needed in the EKF, the functions g(·) and h(·) need to be di↵erentiable at least once with respect to the state vector.

Algorithm 3 Extended Kalman Filter

(26)

2.8 The Extended Rauch-Tung-Striebel Smoother

The extended Rauch-Tung-Striebel smoother (ERTSS) is the corresponding smoother to the EKF and is presented in Algorithm 4 (Särrkä, 2013). The di↵erences com-pared to the RTSS are corresponding to the di↵erences between the KF and the EKF. Just like the RTSS, the ERTSS is started with µs

N = µN and ⌃sN = ⌃N. Algorithm 4 Extended Rauch-Tung-Striebel Smoother

1: function ERTSS(µs_k+1, ⌃s_k+1,µ_k, ⌃_k,u_k) 2: µ_k+1 =g µ_k,uk 3: Gk = @g(_@µµk,uk) k 4: ⌃_k+1= G_k⌃_kGT_k+Q 5: Lk = ⌃kGT_k ⇣⌃_k+1⌘ 1 6: µs_k = µ_k +L_k ⇣µs_k+1 µ_k+1⌘ 7: ⌃s_k = ⌃k+Lk ⇣⌃s_k+1 ⌃_k+1⌘LT_k 8: return µs_k, ⌃s_k 9: end function

2.9 The Unscented Kalman Filter

The UKF uses a set of discrete points, called sigma points, to parameterize the mean and covariance of the posterior. The UKF captures the mean and covariance with an accuracy of at least a second order truncated Taylor expansion (Julier & Uhlmann, 1997). The sigma points are passed directly through the nonlinear func-tion, which means that no Jacobian matrices need to be evaluated. A result of this is that the models do not need to be di↵erentiable. This benefits in that the models could be treated as “black boxes”. The algorithm is presented in Algorithm 5 and is for systems of the form in Equation 2.7 with additive noise.

A state vector with mean µ_k 2 Rnneeds 2n + 1 sigma points. These are

(27)

2.9. THE UNSCENTED KALMAN FILTER matrix and 2 R is a scaling factor given by

= ↵2(n + ) n, (2.25)

0  ↵  1, 0  .

↵controls the spread of the sigma points and  0 to guarantee a positive definite covariance matrix (Van Der Merwe, 2004).

To perform the prediction step of the filter, the sigma points are first passed through the nonlinear process function in line 4, this gives a new set of sigma points, Sk. From this set, the statistics of the prediction is calculated in line 5 and 6, where the weights for the mean is given by

w(i)_m = 8>_<> : n+ for i = 0, 1 2(n+ ) for i , 0 (2.26) and the weights for the covariance by

w(i)_c = 8>_<> : n+ + ⇣ 1 ↵2+ ⌘ for i = 0, 1 2(n+ ) for i , 0. (2.27) 0 is a constant that could be used to include information about the fourth order term in the Taylor expansion. For a Gaussian distribution, = 2 is optimal (Julier, 2002).

A new set of sigma points are now calculated in line 7 with the predicted mean and covariance matrix. This set with predicted sigma points are passed through the measurement model in line 9, which gives a set of sigma points, Yk, for the predicted measurement. The statistics of this set are calculated in line 10 to 12, where ˆy_k is the mean,Sk is the covariance matrix andTk is the cross-covariance matrix for the predicted measurement.

(28)

Algorithm 5 Unscented Kalman Filter

1: function UKF(µ_{k 1}, ⌃_{k 1},u_{k 1},y_k) 2: _X(i)_{k 1} = 8>>>> >>> <>>>> >>> : µ_{k 1} for i = 0, µ_{k 1}+ ✓q(n + )⌃_{k 1} ◆(i) for i = 1, . . . , n µ_{k 1} ✓q(n + )⌃_{k 1}◆(i n) for i = n + 1, . . . ,2n 3: _X_{k 1} = ⇣X(0)_{k 1} X(1)_{k 1} . . . _X(2n)_{k 1}⌘ 4: _S_k =g(X_{k 1},u_{k 1})

5: µ_k = P2n_i=0w_m(i)_S(i)_k

6: ⌃_k =P2n_i=0w_c(i) ⇣S_k(i) µ_k⌘ ⇣_S(i)_k µ_k⌘T+Q

7: _X(i)_k = 8>>>> >>> <>>>> >>> : µ_k for i = 0, µ_k + ✓q(n + )⌃_k ◆(i) for i = 1, . . . , n, µ_k ✓q(n + )⌃_k◆(i n) for i = n + 1, . . . ,2n 8: _X_k = ⇣X(0)_k X(1)_k . . . _X(2n)_k ⌘ 9: _Y_k = h⇣X_k⌘

10: yˆ_k =P2n_i=0w_m(i)Y(i)_k

11: S_k =P2n_i=0w_c(i) ⇣_Y(i)_k yˆ_k⌘ ⇣_Y(i)_k yˆ_k⌘T+R

12: T_k =P2n_i=0w_c(i) ⇣_X(i)_k µ_k⌘ ⇣_Y(i)_k yˆ_k⌘T

(29)

2.10. THE UNSCENTED RAUCH-TUNG-STRIEBEL SMOOTHER

2.10 The Unscented Rauch-Tung-Striebel Smoother

The unscented Rauch-Tung-Striebel smoother (URTSS) is the corresponding smoother to the UKF and was introduced by Särkka (2008). The algorithm is presented in Algorithm 6.

Algorithm 6 Unscented Rauch-Tung-Striebel Smoother

1: function URTSS(µs_k+1, ⌃s_k+1,µ_k, ⌃k,uk) 2: _X(i)_k = 8>>>>> <>>>> > : µ_k for i = 0, µ_k + ⇣p(n + )⌃_k⌘(i) for i = 1, . . . , n µ_k ⇣p(n + )⌃_k⌘(i n) for i = n + 1, . . . ,2n 3: _X_k = ⇣X(0)_k X(1)_k . . . _X(2n)_k ⌘ 4: _Sk+1=g(Xk,uk)

5: µ_k+1=P2n_i=0w_m(i)_S(i)_k+1

6: ⌃_k+1 =P_i=02n w(i)_c ⇣_S_k+1(i) µ_k+1⌘ ⇣_S(i)_k+1 µ_k+1⌘T+Q

(30)

(31)

Chapter 3 Implementation

This chapter describes how the models, filters and smoothers were implemented and which methods that were used.

All software and simulations were implemented in MATLAB (2014) and Simulink with the signal processing toolbox.

The filter performance depends on, to some degree, how nonlinear the process function is and the nonlinearity of the process function depends on the damper setting. Therefore, three di↵erent damping settings were used to compare filter performance. These settings could be seen in Figure 1.1.

3.1 Parameter Identification

The damper models include a damping coefficient and a spring constant that needs to be measured. This section describes how this was done.

3.1.1 Damping Coefficient

The damping coefficients for the linear models were identified with linear interpo-lation of the damping curves in Figure 1.1. The slope of the curve is the damping coefficient and this resulted in the coefficients presented in Table 3.1.

Table 3.1: Linear damping coefficients

Setting Damping coefficient Low damping 4289 N s/m

(32)

CHAPTER 3. IMPLEMENTATION The velocity dependent damping coefficient Equation 2.11 was implemented as a lookup table with data points from Figure 1.1. Values outside and in between the data points were found with linear extrapolation and interpolation. These mea-surements also include the static friction force and therefore no separate term was needed.

3.1.2 Spring Constant

The spring constant determines how strong the position dependent force is. The obvious way to identify its value is to run the damper from one position to another and measure how much the force changes. This is preferably done at low speed to minimize the damping force. This were done by collecting a force measure-ment data series at a velocity of 0.02 m/s and a position span of 50 mm. The relationship between force and position were linear interpolated and this resulted in a spring constant of 628.9 N/m

3.2 Noise Analysis

The filters and smoothers considered in this work assumes that the noise are Gaus-sian, therefore, it were of interest to make an analysis of the noise to see how well the noise could be approximated as Gaussian. To do this, a 30 s long data series were collected at 3184 Hz with the actuator at a static position. No damper was mounted to prevent eventual vibrations to propagate to the load cell. The mean of the data series were subtracted to render the noise unbiased.

A frequency spectrum was generated for each sensor to reveal eventual fre-quency peaks. The frefre-quency spectrum for the position sensor could be seen in Figure 3.1. Three prominent peaks are visible in the figure at 50 Hz, 336 Hz and 1424 Hz. The peak at 50 Hz does most likely come from the switching frequency in the mains electricity. The sources of the other two are unknown.

The spectrum of the force sensor noise reveals one high peak at 50 Hz and a smaller one at 550 Hz. This concludes that the measurement system seems to su↵er from disturbances from the mains electricity.

To see how well the noise matches a Gaussian distribution, a histogram over the measurement values were drawn. For the position signal, this show an approx-imately Gaussian distribution with a negative skewness, see Figure 3.3. The high frequency peaks does not seem to have a↵ected the distribution too much. The probability distribution of the force signal seems to be approximately Gaussian as well, see Figure 3.4. However, the histogram shows some negative kurtosis. The frequency peaks does not seem to have much influence on this distribution either.

(33)

3.2. NOISE ANALYSIS 0 200 400 600 800 1000 1200 1400 1600 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 −6 Frequency/Hz |Position|/m

Figure 3.1: A single sided frequency spectrum of the absolute value of the position noise. Three frequencies stand out with high amplitudes, 50 Hz, 336 Hz and 1424 Hz.

0 200 400 600 800 1000 1200 1400 1600 0 1 2 3 4 5 6 Frequency/Hz |Force|/N

(34)

CHAPTER 3. IMPLEMENTATION −2 −1.5 −1 −0.5 0 0.5 1 1.5 x 10−5 0 2 4 6 8 10 12 x 10 4 Position/m Probability Density/(1/m) Histogram Gaussian approximation

Figure 3.3: A histogram of the position noise and a white Gaussian with the same variance as the noise for reference. The distribution seems to be approximately Gaussian with a negative skewness. −40 −30 −20 −10 0 10 20 30 0 0.01 0.02 0.03 0.04 0.05 0.06 Force/N Probability Density/(1/N) Histogram Gaussian approximation

Figure 3.4: A histogram of the force noise and a white Gaussian with the same variance as the noise for reference. The distribution seems to be approximately Gaussian, but with a lower kurto-sis.

(35)

3.3. DERIVATION OF INPUT SIGNAL

Table 3.2: Statistics of the sensors’ noise

Sensor Variance Skewness Kurtosis Position 2.233 ⇥ 10 5_m2 _-0.1949 _2.349 Force 83.53 N2 _-0.01907 _2.329

3.3 Derivation of Input Signal

As described in subsection 2.4.5, the input to the model was acceleration. Since the reference position signal was known, it was possible to calculate the reference acceleration. However, since the control loop and the time constants in the sys-tem a↵ected the output, the dynamometer could not follow the reference exactly. Hence, a model of the dynamometer was needed. The reason for using a simulated position signal and not a measured signal to generate the input is noise. Noise is amplified with derivation and an input signal that was calculated from measured position would not have been as useful as a simulated signal.

By inspection and after some simulations, it was found that the dynamometer transfer function from reference position to output position were on the form

G(s) = K (s + a)

(s + b)(s + c), (3.1)

where s is the Laplace variable, a, b and c are constants that determine the loca-tions of the zero and poles. K is a constant that controls the static gain. The static gain is given by the final value theorem (Ljung & Glad, 2006) and a step input to the transfer function:

lim s!0sG(s) 1 s = K a bc. (3.2)

The static gain were one, which gives K a

bc = 1 ) K = bc

a . (3.3)

The parameters were found by trial and error, this method resulted in the fol-lowing values:

a = 3.525, b = 3.6, c = 500.

(36)

CHAPTER 3. IMPLEMENTATION 53.5 54 54.5 55 55.5 56 1.04 1.06 1.08 1.1 1.12 1.14 1.16 x 10 −3 Time/s Position/m Measured position Simulated position

Figure 3.5: A comparison between model output and measured output for a selected part of a drive file.

3.4 Comparison Method

The filters’ and smoothers’ performance were evaluated with the mean squared error (MSE) between the estimate and the true value. This is a convenient measure since it could be compared directly with the variance of the sensor noise. The mean squared error is defined as

MSE = 1 N N X k=1 ( ˆxk xk)2, (3.4)

where ˆxk is the estimate and xkis the true value. As usual, N is the total number of samples. It should be mentioned that these comparisons were done with simulated data since it is impossible to know the true value xk in the real system.

3.5 Filter Tuning

(37)

3.6. UKF PARAMETERS

measured sensor variance. The first element is the position variance and the sec-ond is the force variance. The cross-covariance between the noise from the force and position sensor are assumed to be zero, which gives

R = 2.233 ⇥ 10₀ 5 _83.530 !

. (3.5)

This left the covariances for the process noise as tuning parameters. The idea was to make the estimate look “good” with the real data and use the same covariances for the simulated process noise. This was done with the UKF since it possesses the best linearization method and hence should give the lowest model error in the force state. The EKF and the KF must be compensated for their flaws by increas-ing the process noise covariance for the force model. This strategy also gave the covariance for the velocity model. However, the position model has no uncer-tainty since it is a pure integration from position, which is a correct model. But setting this value to zero is not possible since it would make the matrix singular and non-invertible. The solution was to choose a low number, compared to the other covariances, but as high as possible to get a well-conditioned matrix. This gave one covariance matrix for each filter and for the UKF it was

Q = 0 BB@1 ⇥ 10 13 ₀ ₀ 0 1 ⇥ 10 5 ₀ 0 0 3500 1 CCA , (3.6)

for the EKF

Q = 0 BB@1 ⇥ 10 13 ₀ ₀ 0 1 ⇥ 10 5 ₀ 0 0 4000 1 CCA , (3.7)

and for the KF

Q = 0 BB@1 ⇥ 10 13 ₀ ₀ 0 1 ⇥ 10 5 ₀ 0 0 7000 1 CCA . (3.8)

The order of the variances in the covariance matrices are the same as in the state vector; position, velocity and force. The same covariance matrix was used for each filter’s corresponding smoother.

3.6 UKF Parameters

(38)

CHAPTER 3. IMPLEMENTATION and  (see section 2.9 for more information about the parameters). The starting point for tuning these parameters were

↵ =1, = 2,  =0.

The parameters were tuned by trial and error with the goal of achieving a lower MSE. But it turned out that the starting point parameters were as good as it gets.

(39)

Chapter 4 Results and Analysis

This chapter presents, describe and analyze the outcome of this work. First, a performance comparison with simulated data is presented. The simulated results are followed by figures showing the e↵ect of the filters and smoothers on real measured data. The results are presented in figures that are relevant for a real application. Some results were so similar that figures that would have appeared identical are left out; instead a representing case is presented.

4.1 Performance Comparison

The results for each filter for the position state are presented in Figure 4.1. The KF and EKF have the same MSE, while the UKF is a bit lower. The reason for this interesting result is that the UKF is better at incorporating the information from the force sensor to estimate the position. However, if the force signal is left out, or extremely noisy, all filters will be equally good. This is as expected since the position model is linear and there is no benefit from using a filter that could handle nonlinear models. The corresponding results hold for the smoothers. There is a small visible di↵erence between the low damping setting compared to the other settings with the UKF and the URTSS.

The results for the velocity state are shown in Figure 4.2. As seen, since the velocity state is not measured directly, but derived from the position signal, it is very beneficial to use some sort of filter to use information from other states to form a velocity estimate. The reason for the discrepancy between KF, EKF and UKF is the same as described for the position state.

(40)

CHAPTER 4. RESULTS AND ANALYSIS

Unfiltered KF EKF UKF RTSS ERTSS URTSS 0 0.5 1 1.5 2 2.5 x 10 −11 MSE/m 2 High damping Medium damping Low damping

Figure 4.1: A comparison of the position mean squared error for each filter type and each damper setting, including the variance of the position sensor noise as a reference.

Unfiltered KF EKF UKF RTSS ERTSS URTSS 0 1 2 3 x 10 −4 MSE/(m/s) 2 High damping Medium damping Low damping

Figure 4.2: A comparison of the velocity mean squared error for each filter type and each damper setting, including the variance of the time derivative of the position sensor noise as a reference.

(41)

4.2. RESULTS WITH REAL DATA

Unfiltered KF EKF UKF RTSS ERTSS URTSS 0 20 40 60 80 100 120 MSE/N 2 High damping Medium damping Low damping

Figure 4.3: A comparison of the force mean squared error for each filter type and each damper setting, including the variance of the force sensor as a reference.

4.2 Results with Real Data

This section present and illustrate the filtered results with real data. Only the UKF and the URTSS were used to form these results since the di↵erence between other filters and smoothers are so small that they basically look the same as the UKF and the URTSS. The UKF and the URTSS gives a representable view of how the filters and smoothers can improve the measurements.

A comparison between an UKF estimate and a measurement is shown in Fig-ure 4.4 as the velocity with respect to position during one constant velocity stroke. As seen, there is a significant improvement of the signal quality.

The corresponding smoothed result for the velocity with respect to position is shown in Figure 4.5. This figure should be compared with Figure 4.4. The result for the URTSS is much "smoother" than the UKF estimate, hence the name smoother.

The filtered results for the velocity with respect to time during the same stroke as shown before is presented in Figure 4.6. It is clear how the filter improves the measurement and this result could be put in relation to the variance comparison in Figure 4.2 to give an idea of what the improvement looks like in reality.

(42)

CHAPTER 4. RESULTS AND ANALYSIS −0.03 −0.02 −0.01 0 0.01 0.02 0.03 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Position/m Velocity/(m/s) Measurement UKF estimate

Figure 4.4: Filtered velocity with respect to filtered position during a 0.2 m/s stroke.

−0.03 −0.02 −0.01 0 0.01 0.02 0.03 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Position/m Velocity/(m/s) Measurement URTSS estimate

Figure 4.5: Smoothed velocity with respect to smoothed position during a 0.2 m/s stroke.

(43)

4.2. RESULTS WITH REAL DATA 45.5 45.6 45.7 45.8 45.9 46 46.1 46.2 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Time/s Velocity/(m/s) Measurement UKF estimate

Figure 4.6: Filtered velocity with respect to time during a 0.2 m/s stroke.

45.5 45.6 45.7 45.8 45.9 46 46.1 46.2 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Time/s Velocity/(m/s) Measurement URTSS estimate

(44)

CHAPTER 4. RESULTS AND ANALYSIS The UKF estimate of the dynamic damping force with respect to velocity is presented in Figure 4.8. The estimate contains less noise than the measurement, however, the estimate does not seem to conform with the measurement at 0.2 m/s. Figure 4.9 shows the result for the smoothed relationship between damping force and velocity. Here, the estimate seems to have better conformity with the measurement at 0.2 m/s compared to the UKF estimate.

(45)

4.2. RESULTS WITH REAL DATA −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 −4000 −3000 −2000 −1000 0 1000 2000 3000 4000 Velocity/(m/s) Force/N Measurement UKF estimate

Figure 4.8: Filtered force with respect to filtered velocity. The damper were actuated with a 8 Hz sine wave during 1 s. The sine wave had and an amplitude of 0.5 m/s.

−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 −4000 −3000 −2000 −1000 0 1000 2000 3000 4000 Velocity/(m/s) Force/N Measurement URTSS estimate

(46)

(47)

Chapter 5 Discussion and Recommendations

There was almost no improvement of the force signal; this might be solved in various ways. One way is to develop improved models and more accurate input signals to let the filters trust the process more. Another way might be to add more sensors. Depending on hardware and resources, it could be possible to add load cells in series with each other. This would give more information to the filter about the force signal and a lower variance as a result.

The measurements could be improved by investigate the noise sources, this is out of the scope of this work and is only discussed briefly. As seen in the noise analysis in section 3.2, the measurement system is probably a↵ected by the mains supply. This might be difficult or easy to resolve, depending on how the noise is coupled to the system. Even if the noise sources are eliminated and the measurement system is improved, old measurements remain noisy. The filters however, are able to improve data that is already collected.

The filtering techniques could be used to improve the dynamometer perfor-mance. With a filter in the control loop, the system benefits from the availability of a velocity estimate. I.e. the filter acts as a state observer. With a proper de-signed feed forward controller, this could result in higher controller performance.

5.1 Rejected Ideas

This section present and discuss ideas that were rejected, either because they did not work or because a better solution were found.

(48)

CHAPTER 5. DISCUSSION AND RECOMMENDATIONS Another theory was that the mounting assembly of the dynamometer was not sti↵ enough. Therefore, measurements were made with an accelerometer mounted on the cross head as well. But as Appendix A shows, including a model of the mounting assembly would have resulted in a very sti↵ problem with numerical is-sues. Besides, it would have been necessary to develop some sort of measurement method to verify the models, maybe with strain gauges.

5.2 Recommendations

If the application is o✏ine and a good model of the damper is available, it is recommended to use the URTSS since it is the most accurate in this case. If the same conditions are met, but the estimate is needed at real time, the UKF is the first choice. If there is no damper model or if no force measurements are available, the system is reduced and the model becomes linear. In this case, there is no need for a filter that is capable of handling nonlinear models. Therefore the RTSS is a suitable choice for o✏ine applications and the KF is the appropriate filter for online implementations.

(49)

Chapter 6 Conclusion

(50)

(51)

Appendix A

Dynamometer Eigenfrequencies

This appendix analyze the dynamic properties of the mounting assembly of the dynamometer and will act as a motivation for the exclusion of these dynamics in the system model derived in section 2.4. The mounting assembly of the dy-namometer could be seen as three flexible beams, one horizontal beam supported by two vertical columns.

A.1 Cross-Head

The horizontal beam is called cross-head and is modeled as a simply supported beam with the length L, the density ⇢, the flexural rigidity EI and the cross sec-tional area A. The eigenfrequencies !n is given by the bending wave equation (Bodén et al., 2001): d4 n d x4 !2n ⇢A EI n= 0, n = 1,2, . . . , (A.1) where nis the mode shape for eigenfrequency n. With the given boundary con-ditions, the mode shape is

n=sin(knx) , (A.2)

where kn is the wave number for mode n and x denotes the position where the deflection is measured. (Bodén et al., 2001). Only the first mode is considered relevant and this gives k1 = ⇡_L (Bodén et al., 2001). Equation A.2 inserted in Equation A.1 gives

k_n4sin(knx) !2n⇢A_EI sin(knx) = 0. (A.3) The force F is acting on the beam at half its length, which gives

x = L

(52)

APPENDIX A. DYNAMOMETER EIGENFREQUENCIES When Equation A.3 is evaluated at this position and with the relevant wave num-ber, it becomes ⇡4 L4 !21 ⇢A EI =0. (A.5) This gives !2₁= ⇡ 4_EI L4_⇢A. (A.6)

The eigenfrequency fb of the beam is given by plugging in physical data of the cross-head into Equation A.6. This gives

!₁2_{⇡ 1.4 ⇥ 10}4_{rad/s ) f}_b_{⇡ 2200 Hz,} (A.7) If included in the system model, this would have resulted in a very fast pole and a sti↵ system. Besides, the eigenfrequency is out of the bandwidth of the measure-ment system and therefore not included in the model.

A.2 Vertical Columns

The two vertical columns are assumed to vibrate in phase with each other and therefore modeled as one beam with twice the cross sectional area. The boundary conditions are considered to be fixed in the lower end and with a point mass at the other (the horizontal beam). Only quasi-longitudinal waves are considered and this gives the transcendental equation

cot = _⇢SLm , (A.8)

where m is the mass of the horizontal beam, ⇢ is the density of the vertical beam, S is the cross sectional area and L is the length (Bodén et al., 2001). is used in the following equation to calculate the eigenfrequency:

fn = _2⇡LcL n, (A.9)

where cL is the phase velocity for quasi-longitudinal waves and fnis the eigenfre-quency for the nth modal shape (Bodén et al., 2001). This gives

f1⇡ 1200 Hz,

which will a↵ect the system more than the horizontal beam. However, even this is considered irrelevant to include in the system model.

(53)

Appendix B

Gas Reservoir Spring Force

This appendix presents an analysis of the gas spring force from a thermodynam-ical point of view. Denote the pressure in the gas reservoir as p, the volume of the gas as V and the position coordinate z (origin at fully extension). According to Dixon (2007, p. 269), the process could be seen as adiabatic for fast dynamics and isothermal for slow dynamics. A more accurate model could be an polytropic process with an identified polytropic index, as discussed by Warner (1996). How-ever, an adiabatic process is considered accurate enough for this application. The following holds for an adiabatic process:

pV = constant, (B.1)

where is the heat capacity ratio, in this case = 1.4 since the gas reservoir is filled with nitrogen. The gas volume is

V = V0 Arz, (B.2)

where V0is the initial volume and Ar is the cross sectional area of the piston rod. Equation B.1 and Equation B.2 yields

p0V₀ = p(V0 Arz) ) p = _V V0 0 Arz

!

p0. (B.3)

The gas force contribution Fk is given by

F_k = A_r(p p_a), (B.4)

where pa is the atmospheric pressure. The nonlinear spring force could be lin-earized around a stationary point z⇤_as

(54)

APPENDIX B. GAS RESERVOIR SPRING FORCE For readability reasons, Equation B.5 is rewritten as

F_k_{(z) ⇡ k(z + Z}_p), (B.6)

where the spring constant is

k = p0V0 A 2 r

(V0 Arz⇤) +1. (B.7)

The equivalent preload Zpis derived as

k Zp =Fk(z⇤) kz⇤ ) Zp = Fk(z ⇤₎

k z⇤. (B.8)

If the physical data of the considered damper is inserted into Equation B.7, the calculated spring sti↵ness is

k ⇡ 250 N/m. (B.9)

This result is consistent with Dixon (2007, p. 269), who claims that the sti↵-ness usually is about 200 N/m. However, this is not even close to the measured value of 628.9 N/m from section 3.1. This discrepancy is left unexplained and the measured value was used in the models.

(55)

Bibliography

Anderson, B. D. O. & Moore, J. B. (1979). Optimal Filtering. Prentice-Hall, Inc. Arasaratnam, I. & Haykin, S. (2009). Cubature kalman filters. Automatic Control,

IEEE Transactions on, 54(6), 1254–1269.

Bodén, H., Carlsson, U., Glav, R., Wallin, H., & Åbom, M. (2001). Ljud och Vi-brationer. Marcus Wallenberg Laboratoriet för Ljud- och Vibrationsforskning. DeCarlo, R. A. (1989). Linear systems: A State Variable Approach with

Numeri-cal Implementation. Prentice Hall.

Dixon, J. C. (2007). The Shock Absorber Handbook. John Wiley & Sons, Ltd. Doucet, A. & Johansen, A. M. (2009). A tutorial on particle filtering and

smooth-ing: Fifteen years later. In D. Crisan & B. Rozovskii (Eds.), The Oxford Hand-book of Nonlinear Filtering, volume 12 (pp. 656–704). Oxford: Oxford Uni-versity Press.

Gordon, N., Salmond, D., & Smith, A. F. M. (1993). Novel approach to nonlinear/non-gaussian bayesian state estimation. Radar and Signal Process-ing, IEE Proceedings F, 140(2), 107–113.

Ito, K. & Xiong, K. (2000). Gaussian filters for nonlinear filtering problems. Automatic Control, IEEE Transactions on, 45(5), 910–927.

Julier, S. (2002). The scaled unscented transformation. In American Control Conference, 2002. Proceedings of the 2002, volume 6, (pp. 4555–4559 vol.6). Julier, S. J. & Uhlmann, J. K. (1997). New extension of the kalman filter to

nonlinear systems.

(56)

BIBLIOGRAPHY Karnopp, D. (1985). Computer simulation of stick-slip friction in mechanical dy-namic systems. Journal of dydy-namic systems, measurement, and control, 107(1), 100–103.

Ljung, L. & Glad, T. (2004). Modellbygge och simulering. Studentlitteratur. Ljung, L. & Glad, T. (2006). Reglerteknik. Studentlitteratur.

MATLAB (2014). MATLAB and Signal Processing Toolbox Release R2014a. Natick, Massachusetts: The MathWorks Inc.

Nørgaard, M., Poulsen, N. K., & Ravn, O. (2000). New developments in state estimation for nonlinear systems. Automatica, 36(11), 1627 – 1638.

Rauch, H. E., Striebel, C., & Tung, F. (1965). Maximum likelihood estimates of linear dynamic systems. AIAA journal, 3(8), 1445–1450.

Särkka, S. (2008). Unscented rauch–tung–striebel smoother. Automatic Control, IEEE Transactions on, 53(3), 845–849.

Särrkä, S. (2013). Bayesian Filtering and Smoothing. Cambridge University Press.

Thrun, S., Burgard, W., & Fox, D. (2005). Probabilistic Robotics. MIT Press. Van Der Merwe, R. (2004). Sigma-point Kalman filters for probabilistic inference

in dynamic state-space models. PhD thesis, University of Stellenbosch.

Warner, B. (1996). An analytical and experimental investigation of high perfor-mance suspension dampers. PhD thesis, Concordia University.

(57)

(58)