5.1.1.1 The State Space Model

The state at time i is denoted by x(i) ∈ X, where X is the state space. For finite discrete state space, X = Ω = {ω₁, …, ω_k}, where ω_k is the kth symbol (label or class) out of K possible classes. For infinite discrete state space, X = Ω = {ω₁, …, ω_k, …}, where 1 ≤ k < ∞. For real-valued vectors with dimension m, we have . In all these cases we assume that the states can be modelled as random variables.

For a causal system with input vector u(i), an mth-order state space model takes the form

(5.1)

where w(i) is the system noise (also known as process noise).

Suppose for a moment that we have observed the state of a process during its whole history, that is from the beginning of time up to the present. In other words, the sequence of states x(0), x(1), …, x(i) are observed and as such fully known (i denotes the present time). In addition, suppose that – using this sequence – we want to estimate (predict) the next state x(i + 1). We need to evaluate the conditional probability density p(x(i + 1)|x(0), x(1), …, x(i)). Once this probability density is known, the application of the theory in Chapters 3 and 4 will provide the optimal estimate of x(i + 1). For instance, if X is a real-valued vector space, the Bayes estimator from Chapter 4 provides the best prediction of the next state (the density p(x(i + 1)|x(0), x(1), …, x(i)) must be used instead of the posterior density).

Unfortunately, the evaluation of p(x(i + 1)|x(0), x(1), …, x(i)) is a nasty task because it is not clear how to establish such a density in real-world problems. Things become much easier if we succeed to define the state such that the so-called Markov condition applies:

(5.2)

The probability of x(i + 1) depends solely on x(i) and not on the past states. In order to predict x(i + 1), knowledge of the full history is not needed. It suffices to know the present state. If the Markov condition applies, the state of a physical process is a summary of the history of the process.

Example 5.1 Motion tracking

Suppose a particle P is moving on a two-dimensional plane whose position at the time i is denoted by the coordinate vector (x(i), y(i)). Also at the time moment i the speed of P is expressed as (v_x(i), v_y(i)), where v_x(i) and v_y(i) represent the velocity components along the x-axis and y-axis, respectively. The continuous time t is equidistantly sampled with period Δ. If Δ is small enough, it can be assumed that P’s moving speed is approximately uniform during a sampling interval. Thus the position of P can be recursively calculated as

(5.3)

where w_x(i) and w_y(i) are the approximation errors. Define the state variable as

Then we have

(5.4)

where

and

Equation (5.4) is in the form of (5.1). Furthermore, it is easy to see that the Markov condition (5.2) holds if {w_x(i)} and {w_y(i)} are independent white random sequences with normal distribution (see Equation (5.14) below).

Example 5.2 The density of a substance mixed with a liquid

Mixing and diluting are tasks frequently encountered in the food industries, paper industry, cement industry, and so. One of the parameters of interest during the production process is the density D(t) of some substance. It is defined as the fraction of the volume of the mix that is made up by the substance.

Accurate models for these production processes soon involve a large number of state variables. Figure 5.1 is a simplified view of the process. It is made up by two real-valued state variables and one discrete state. The volume V(t) of the liquid in the barrel is regulated by an on/off feedback control of the input flow f₁(t) of the liquid: f₁(t) = f₀x(t). The on/off switch is represented by the discrete state variable x(t) ∈ {0, 1}. A hysteresis mechanism using a level detector (LT) prevents jitter of the switch. The x(t) switches to the ‘on’ state (=1) if V(t) < V_low and switches back to the ‘off’ state (=0) if V(t) > V_high.

The rate of change of the volume is with f₂(t) the volume flow of the substance and f₃(t) the output volume flow of the mix. We assume that the output flow is governed by Torricelli's law: . The density is defined as , where V_S(t) is the volume of the substance. The rate of change of V_S(t) is . After some manipulations the following system of differential equations appears:

A discrete time approximation (notation: V(i) ≡ V(iΔ), D(i) ≡ D(iΔ), and so on) is

(5.5)

This equation is of the type x(i + 1) = f(x(i), u(i), w(i)) with x(i) = [V(i) D(i) x(i)]^T. The elements of the vector u(i) are the known input variables, which is the non-random part of f₂(i). The vector w(i) contains the random input, i.e. the random part of f₂(i). The probability density of x(i + 1) depends on the present state x(i), but not on the past states.

Figure 5.1 shows a realization of the process. Here, the substance is added to the volume in chunks with an average volume of 10 litres and at random points in time.

If the transition probability density p(x(i + 1)|x(i)) is known together with the initial probability density p(x(0)), then the probability density at an arbitrary time can be determined recursively:

(5.6)

The joint probability density of the sequence x(0), x(1), …, x(i) follows readily:

(5.7)

5.1.1.2 The Measurement Model

Unfortunately in most cases we cannot directly observe the system states. Instead we can only observe the system outputs such as the data from the sensor in relation to the state. This means, in addition to the state space model, we also need a measurement model that describes the data from the system output. Suppose that at moment i the measurement data are z(i) ∈ Z, where Z is the measurement space. For the real-valued state variables, the measurement space is often a real-valued vector space, that is for some positive integer n. For the discrete case, one often assumes that the measurement space is also discrete. The measurement model is described by

(5.8)

where v(i) is the measurement noise.

The probabilistic model of the measurement system is fully defined by the conditional probability density p(z(i)|x(0), …, x(i), z(0), …, z(i − 1)). We assume that the sequence of measurements starts at time i = 0. In order to shorten the notation, the sequence of all measurements up to the present will be denoted by

We restrict ourselves to memoryless measurement systems which is systems where z(i) depends on the value of x(i), but not on previous states nor on previous measurements. In other words:

(5.9)

The state model (5.1) together with the measurement model (5.8) forms the general discrete-time state space model for state estimation. However, this general model is not very useful until we move to more specific situations, which will be discussed later in this chapter.

5.2.1.1 Linear-Gaussian State Space Models

The state model is said to be linear if the state space model (5.1) can be expressed by a so-called linear system equation (or linear state equation, linear plant equation, linear dynamic equation):

(5.13)

where F(i) is the system matrix. It is an m × m matrix where m is the dimension of the state vector; m is called the order of the system. The vector u(i) is the control vector (input vector) of dimension l. Usually, the vector is generated by a controller according to some control law. As such the input vector is a deterministic signal that is fully known, at least up to the present. L(i) is the gain matrix of dimension m × l. Sometimes the matrix is called the distribution matrix as it distributes the control vector across the elements of the state vector.

w(i) is the process noise (system noise, plant noise). It is a sequence of random vectors of dimension m. The process noise represents the unknown influences on the system, for instance, formed by disturbances from the environment. The process noise can also represent an unknown input/control signal. Sometimes process noise is also used to take care of modelling errors. The general assumption is that the process noise is a white random sequence with normal distribution. The term ‘white’ is used here to indicate that the expectation is zero and the autocorrelation is governed by the Kronecker delta function:

(5.14)

where C_w(i) is the covariance matrix of w(i). Since w(i) is supposed to have a normal distribution with zero mean, C_w(i) defines the density of w(i) in full.

The initial condition of the state model is given in terms of the expectation E[x(0)] and the covariance matrix C_x(0). In order to find out how these parameters of the process propagate to an arbitrary time i, the state equation (5.13) must be used recursively:

(5.15)

The first equation follows from E[w(i)] = 0. The second equation uses the fact that the process noise is white, that is E[w(i)w^T(j)] = 0 for i ≠ j (see Equation (5.14)).

If E[x(0)] and C_x(0) are known, then Equation (5.15) can be used to calculate E[x(1)] and C_x(1). From that, by reapplying (5.15), the next values, E[x(2)] and C_x(2), can be found, and so on. Thus, the iterative use of Equation (5.15) gives us E[x(i)] and C_x(i) for arbitrary i > 0.

In the special case, where neither F(i) nor C(i) depends on i, the state space model is time invariant. The notation can then be shortened by dropping the index, that is F and C_w. If, in addition, F is stable (i.e. the magnitudes of the eigenvalues of F are all less than one), the sequence C_x(i), i = 0, 1, … converges to a constant matrix. The balance in Equation (5.15) is reached when the decrease of C_x(i) due to F compensates for the increase of C_x(i) due to C_w. If such is the case, then

(5.16)

which is the discrete Lyapunov equation.

5.2.1.2 Prediction

Equation (5.15) is the basis for prediction. Suppose that at time i an unbiased estimate is known together with the associated error covariance C_e(i), which was introduced in Section 4.2.2. The best predicted value (MMSE) of the state for ℓ samples ahead of i is obtained by the recursive application of Equation (5.15). The recursion starts with and terminates when E[x(i + ℓ)] is obtained. The covariance matrix C_e(i) is a measure of the magnitudes of the random fluctuations of x(i) around . As such it is also a measure of uncertainty. Therefore, the recursive usage of Equation (5.15) applied to C_e(i) gives C_e(i + ℓ), that is the uncertainty of the prediction. With that, the recursive equations for the prediction become

(5.17)

Example 5.3 Prediction of motion tracking Example 5.1 formulated a simple motion tracking problem with the state equation (5.4), where

and

According to Equation (5.17), the one step ahead prediction (ℓ = 0) can be expressed as

(5.18)

where represents the input term. The prediction error covariance at the time i is

(5.19)

where is the covariance of , which is

if is a white random sequence. Suppose the initial prediction error covariance is 0, that is

Then can be generated recursively according to Equation (5.19). By Equation (5.4), can be set to 0. Also suppose the sample interval is 1, that is

This example will be extended to a complete example of online estimation with a discrete Kalman filter by establishing the measurement model in Example 5.4.

5.2.1.3 Linear-Gaussian Measurement Models

A special form of Equation (5.8) is the linear measurement model, which takes the following form:

(5.20)

where H(i) is the so-called measurement matrix, which is an n × m matrix and v(i) is the measurement noise. It is a sequence of random vectors of dimension n. Obviously, the measurement noise represents the noise sources in the sensory system. Examples are thermal noise in a sensor and the quantization errors of an AD converter.

The general assumption is that the measurement noise is a zero mean, white random sequence with normal distribution. In addition, the sequence is supposed to have no correlation between the measurement noise and the process noise:

(5.21)

where C_v(i) is the covariance matrix of v(i). C_v(i) specifies the density of v(i) in full.

The measurement system is time invariant if neither H nor C_v depends on i.

5.2.1.4 The Discrete Kalman Filter

The concepts developed in the previous section are sufficient to transform the general scheme presented in Section 5.1 into a practical solution. In order to develop the estimator, first the initial condition valid for i = 0 must be established. In the general case, this condition is defined in terms of the probability density p(x(0)) for x(0). Assuming a normal distribution for x(0) it suffices to specify only the expectation E[x(0)] and the covariance matrix C_x(0). Hence, the assumption is that these parameters are available. If not, we can set E[x(0)] = 0 and let C_x(0) approach infinity, that is C_x(0) → ∞ I. Such a large covariance matrix represents the lack of prior knowledge.

The next step is to establish the posterior density p(x(0)|z(0)) from which the optimal estimate for x(0) follows. At this point, we enter the loop of Figure 5.2. Hence, we calculate the density p(x(1)|z(0)) of the next state and process the measurement z(1) resulting in the updated density p(x(1)|z(0), z(1)) = p(x(1)|Z(1)). From that, the optimal estimate for x(1) follows. This procedure has to be iterated for all of the next time cycles.

The representation of all the densities that are involved can be given in terms of expectations and covariances. The reason is that any linear combination of Gaussian random vectors yields a vector that is also Gaussian. Therefore, both p(x(i + 1)|Z(i)) and p(x(i)|Z(i)) are fully represented by their expectations and covariances. In order to discriminate between the two situations a new notation is needed. From now on, the conditional expectation E[x(i)|Z(j)] will be denoted by . It is the expectation associated with the conditional density p(x(i)|Z(j)). The covariance matrix associated with this density is denoted by C(i|j).

The update, that is the determination of p(x(i)|Z(i)) given p(x(i)|Z(i − 1)), follows from Section 4.1.5, where it has been shown that the unbiased linear MMSE estimate in the linear-Gaussian case equals the MMSE estimate, and that this estimate is the conditional expectation. Application of Equations (4.33) and (4.45) to (5.20) and (5.21) gives

(5.22)

The interpretation is as follows: is the predicted measurement. It is an unbiased estimate of z(i) using all information from the past. The so-called innovation matrix S(i) represents the uncertainty of the predicted measurement. The uncertainty is due to two factors: the uncertainty of x(i) as expressed by C(i|i − 1) and the uncertainty due to the measurement noise v(i) as expressed by C_v(i). The matrix K(i) is the Kalman gain matrix. This matrix is large when S(i) is small and C(i|i − 1)H^T(i) is large, that is when the measurements are relatively accurate. When this is the case, the values in the error covariance matrix C(i|i) will be much smaller than C(i|i − 1).

The prediction, that is the determination of p(x(i + 1)|Z(i)) given p(x(i)|Z(i)), boils down to finding out how the expectation and the covariance matrix C(i|i) propagate to the next state. Using equations (5.13) and (5.15) we have

(5.23)

At this point, we increment the counter and x(i + 1|i) and C(i + 1|i) become x(i|i − 1) and C(i|i − 1). These recursive equations are generally referred to as the discrete Kalman filter (DKF).

In the Gaussian case, it does not matter much which optimality criterion we select. MMSE estimation, MMAE estimation and MAP estimation yield the same result, that is the conditional mean. Hence, the final estimate is found as . It is an absolute unbiased estimate and its covariance matrix is C(i|i). Therefore, this matrix is often called the error covariance matrix.

In the time invariant case, and assuming that the Kalman filter is stable, the error covariance matrix converges to a constant matrix. In that case, the innovation matrix and the Kalman gain matrix become constant as well. The filter is said to be in the steady state. The steady state condition simply implies that C(i + 1|i) = C(i|i − 1). If the notation for this matrix is shortened to P, then Equations (5.22) and (5.23) lead to the following equation:

(5.24)

Equation (5.24) is known as the discrete algebraic Ricatti equation. Usually, it is solved numerically. Its solution implicitly defines the steady state solution for S, K and C.

Example 5.4 Application to the motion tracking In this example we reconsider the motion tracking described in Example 5.3. We model the measurements with the model: where is the constant matrix

Figure 5.3 shows the results of the motion tracking. The covariance of the measurement noise is . The filter is initiated with to be the actual position of the object. Recall that in Example 5.3, C_x(0) is set to be 1. It can be seen from Figure 5.3 that if both system error and measurement error are white noise with covariance 1, the motion tracking works quite well.

Listing 5.1 is the code of this example.

Listing 5.1 Discrete Kalman filter for motion tracking

clear
x=(1:80);
for idx = 1: length(x)
 y(idx)=(idx-22)*(idx-22)/20;
end
y=2*randn(1,length(x)) + y;
 
kalman = [];
for idx = 1: length(x)
   location = [x(idx),y(idx),0,0];
   if isempty(kalman)
       stateModel = [1 0 1 0;0 1 0 1;0 0 1 0;0 0 0 1];
       measurementModel = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];
       kalman = vision.KalmanFilter(stateModel,measurementModel,
       'ProcessNoise',1,'MeasurementNoise',1);
       kalman.State = location;
   else
       trackedLocation = predict(kalman);
       plot(location(1),location(2),'k+');
       trackedLocation = correct(kalman,location);
       pause(0.2);
       plot(trackedLocation(1),trackedLocation(2),'ro');
   end
end

For comparison, Figure 5.4 depicts the tracking results when the system noise becomes 1 × 10⁻⁴ (=0.0001) while the measurement noise becomes 4. This time the performance is much worse after time step 28.

5.2.2.1 The Linearized Kalman Filter

The simplest approximation occurs when the system equations are linearized around some fixed value of x(i). This approximation is useful if the system is time invariant and stable, and if the states swing around an equilibrium state. Such a state is the solution of

Defining , the linear approximation of the state equation (5.25) becomes . After some manipulations:

(5.26)

By interpreting the term as a constant control input and by compensating the offset term in the measurement vector, these equations become equivalent to Equations (5.13) and (5.20). This allows the direct application of the DKF as given in Equations (5.22) and (5.23).

Many practical implementations of the discrete Kalman filter are inherently linearized Kalman filters because physical processes are seldom exactly linear and often a linear model is only an approximation of the real process.

Example 5.5 A linearized model for volume density estimation In Section 5.1.1 we introduced the non-linear, non-Gaussian problem of the volume density estimation of a mix in the process industry (Example 5.2). The model included a discrete state variable to describe the on/off regulation of the input flow. We will now replace this model by a linear feedback mechanism:

The liquid input flow has now been modelled by f₁(i) = α(V₀ − V(i)) + w₁(i). The constant V₀ = V_ref + (c − f₂)/α (with realizes the correct mean value of V(i). The random part w₁(i) of f₁(i) establishes a first-order AR model of V(i), which is used as a rough approximation of the randomness of f₁(i). The substance input flow f₂(i), which in Example 5.2 appears as chunks at some discrete points of time, is now modelled by , that is a continuous flow with some randomness.

The equilibrium is found as the solution of V(i + 1) = V(i) and D(i + 1) = D(i). The results are

The expressions for the Jacobian matrices are

The considered measurement system consists of two sensors:

• A level sensor that measures the volume V(i) of the barrel.
• A radiation sensor that measures the density D(i) of the output flow.

The latter uses the radiation of some source, for example X-rays, that is absorbed by the fluid. According to Beer–Lambert's law, P_out = P_inexp ( − μD(i)), where μ is a constant depending on the path length of the ray and on the material. Using an optical detector the measurement function becomes z = Uexp ( − μD) + v with U a constant voltage. With that, the model of the two sensors becomes

(5.27)

with the Jacobian matrix:

The best fitted parameters of this model are as follows:

Volume	Substance	Output	Measurement
control	flow	flow	system
Δ = 1 (s)		Vref = 4000 (l)
V0 = 4001 (l)		c = 1 (l/s)	U = 1000 (V)
α = 0.95 (1/s)			μ = 100

Figure 5.5 shows the real states (obtained from a simulation using the model from Example 5.2), observed measurements, estimated states and estimation errors. It can be seen that:

• The density can only be estimated if the real density is close to the equilibrium. In every other region, the linearization of the measurement is not accurate enough.
• The estimator is able to estimate the mean volume but cannot keep track of the fluctuations. The estimation error of the volume is much larger than indicated by the 1σ boundaries (obtained from the error covariance matrix). The reason for the inconsistent behaviour is that the linear-Gaussian AR model does not fit well enough.

5.2.2.2 The Extended Kalman Filter

A straightforward generalization of the linearized Kalman filter occurs when the equilibrium point is replaced with a nominal trajectory , recursively defined as

Although the approach is suitable for time variant systems, it is not often used. There is another approach with almost the same computational complexity, but with better performance. That approach is the extended Kalman filter (EKF).

Again, the intention is to keep track of the conditional expectation x(i|i) and the covariance matrix C(i|i). In the linear-Gaussian case, where all distributions are Gaussian, the conditional mean is identical to the MMSE estimate (= minimum variance estimate), the MMAE estimate and the MAP estimate (see Section 4.1.3). In the present case, the distributions are not necessarily Gaussian and the solutions of the three estimators do not coincide. The extended Kalman filter provides only an approximation of the MMSE estimate.

Each cycle of the extended Kalman filter consists of a ‘one step ahead’ prediction and an update, as before. However, the tasks are much more difficult now because the calculation of, for instance, the ‘one step ahead’ expectation:

requires the probability density p(x(i + 1)|Z(i)) (see Equation (5.10)). However, as stated before, it is not clear how to represent this density. The solution of the EKF is to apply a linear approximation of the system function. With that, the ‘one step ahead’ expectation can be expressed entirely in terms of the moments of p(x(i)|Z(i)).

The EKF uses linear approximations of the system functions using the first two terms of the Taylor series expansions. Suppose that at time i we have the updated estimate and the associated approximate error covariance matrix C(i|i). The word ‘approximate’ expresses the fact that our estimates are not guaranteed to be unbiased due to the linear approximations. However, we do assume that the influence of the errors induced by the linear approximations is small. If the estimation error is denoted by e(i), then

(5.28)

In these expressions, is our estimate. It is available at time i and therefore is deterministic. Only e(i) and w(i) are random. Taking the expectation on both sides of (5.28), we obtain approximate values for the ‘one step ahead’ prediction:

We have approximations instead of equalities for two reasons. First, we neglect the possible bias of , that is a non-zero mean of e(i). Second, we ignore the higher-order terms of the Taylor series expansion.

Upon incrementing the counter, becomes and we now have to update the prediction by using a new measurement z(i) in order to get an approximation of the conditional mean . First we calculate the predicted measurement based on using a linear approximation of the measurement function, that is . Next, we calculate the innovation matrix S(i) using that same approximation. Then we apply the update according to the same equations as in the linear-Gaussian case.

The predicted measurements are

The approximation is based on the assumption that E[e(i|i − 1)]≅0 and on the Taylor series expansion of h( · ). The innovation matrix becomes

From this point on, the update continues as in the linear-Gaussian case (see Equation (5.22)):

Despite the similarity of the last equation with respect to the linear case, there is an important difference. In the linear case, the Kalman gain K(i) depends solely on deterministic parameters: H(i), F(i), C_w(i), C_v(i) and C_x(0). It does not depend on the data. Therefore, K(i) is fully deterministic. It could be calculated in advance instead of online. In the EKF, the gains depend upon the estimated states through , and thus also upon the measurements z(i). Therefore, the Kalman gains are random matrices. Two runs of the extended Kalman filter in two repeated experiments lead to two different sequences of the Kalman gains. In fact, this randomness of K(i) can cause unstable behaviour.

Example 5.6 The extended Kalman filter for volume density estimation Application of the EKF to the density estimation problem introduced in Example 5.2 and represented by a linear-Gaussian model in Example 5.5 gives the results as shown in Figure 5.6. Compared with the results of the linearized KF (Figure 5.5) the density errors are now much more consistent with the 1σ boundaries obtained from the error covariance matrix. However, the EKF is still not able to cope with the non-Gaussian disturbances of the volume.

Note also that the 1σ boundaries do not reach a steady state. The filter remains time variant, even in the long term.

5.2.2.3 The Iterated Extended Kalman Filter

A further improvement of the update step in the extended Kalman filter is within reach if the current estimate is used to get an improved linear approximation of the measurement function yielding an improved predicted measurement . In turn, such an improved predicted measurement can improve the current estimate. This suggests an iterative approach.

Let be the predicted measurement in the ℓth iteration and let be the ℓth improvement of . The iteration is initiated with . A naive approach for the calculation of simply uses a relinearization of h( · ) based on :

Hopefully, the sequence , with ℓ = 0, 1, 2, …, converges to a final solution.

A better approach is the so-called iterated extended Kalman filter (IEKF). Here, the approximation is made that both the predicted state and the measurement noise are normally distributed. With that, the posterior probability density (Equation (5.11))

being the product of two Gaussians, is also a Gaussian. Thus, the MMSE estimate coincides with the MAP estimate and the task is now to find the maximum of p(x(i)|Z(i)). Equivalently, we maximize its logarithm. After the elimination of the irrelevant constants and factors, it all boils down to minimizing the following function w.r.t. x:

(5.29)

For brevity, the following notation has been used:

The strategy to find the minimum is to use the Newton–Raphson iteration starting from . In the ℓth iteration step, we already have an estimate obtained from the previous step. We expand f(x) in a second-order Taylor series approximation:

(5.30)

where ∂f/∂x is the gradient and ∂²f/∂x² is the Hessian of f(x) (see Appendix B.4). The estimate is the minimum of the approximation. It is found by equating the gradient of the approximation to zero. Differentiation of Equation (5.30) w.r.t. x gives

(5.31)

The Jacobian and Hessian of Equation (5.29), in explicit form, are

(5.32)

where is the Jacobian matrix of h(x) evaluated at . Substitution of Equation (5.32) in (5.31) yields the following iteration scheme:

(5.33)

The result after one iteration, that is , is identical to the ordinary extended Kalman filter. The required number of further iterations depends on how fast converges. Convergence is not guaranteed, but if the algorithm converges, usually a small number of iterations suffice. Therefore, it is common practice to fix the number of iterations to some practical number L. The final result is set to the last iteration, that is .

Equation (4.44) shows that the factor (C^{− 1}_p + H_ℓ^TC^{− 1}_vH_ℓ)^{− 1} is the error covariance matrix associated with x(i|i):

This insight gives another connotation to the last term in Equation (5.33) because, in fact, the term C(i|i)H^T_ℓC_v^{− 1} can be regarded as the Kalman gain matrix K_ℓ during the ℓth iteration (see Equation (4.20)).

Example 5.7 The iterated EKF for volume density estimation In the previous example, the EKF was applied to the density estimation problem introduced in Example 5.2. The filter was initiated with the equilibrium state as prior knowledge, that is . Figure 5.7(b) shows the transient that occurs if the EKF is initiated with E[x(0)] = [2000 0]^T. It takes about 40(s) before the estimated density reaches the true densities. This slow transient is due to the fact that in the beginning the linearization is poor. The iterated EKF is of much help here. Figure 5.7(c) shows the results. From the first measurement, the estimated density is close to the real density. There is no transient.

The extended Kalman filter is widely used because for a long period of time no viable alternative solution existed. Nevertheless, it has numerous disadvantages:

• It only works well if the various random vectors are approximately Gaussian distributed. For complicated densities, the expectation–covariance representation does not suffice.
• It only works well if the non-linearities of the system are not too severe because otherwise the Taylor series approximations fail. Discontinuities are deadly for the EKF's proper functioning.
• Recalculating the Jacobian matrices at every time step is computationally expensive.
• In some applications, it is too difficult to find the Jacobian matrix analytically. In these cases, numerical approximations of the Jacobian matrix are needed. However, this introduces other types of problems because now the influence of having approximations rather than the true values comes in.
• In the EKF, the Kalman gain matrix depends on the data. With that, the stability of the filter is no longer assured. Moreover, it is very hard to analyse the behaviour of the filter.
• The EKF does not guarantee unbiased estimates. In addition, the calculated error covariance matrices do not necessarily represent the true error covariances. The analysis of these effects is also hard.

5.3.3.1 Individually Most Likely State

Here the strategy is to determine the posterior probability p(x(i)|Z(I)) and then to determine the MAP estimate: x(i|I) = arg max p(x(i)|Z(I)). As stated earlier, this method minimizes the error probabilities of the individual states. As such, it maximizes the expected number of correctly estimated states.

Section 5.3.1 discussed the forward algorithm, a recursive algorithm for the calculation of the probability p(x(i), Z(i)). We now introduce the backward algorithm, which calculates the probability p(z(i + 1), …, z(I)|x(i)). During each recursion step of the algorithm, the probability p(z(j), …, z(I)|x(j − 1)) is derived from p(z(j + 1), …, Z(I)|x(j)). The recursion proceeds as follows:

The algorithm starts with j = I and proceeds backwards in time, that s I, I − 1, I − 2, … until finally j = i + 1. In the first step, the expression p(z(I + 1)|x(I)) appears. Since that probability does not exist (because z(I + 1) is not available), it should be replaced by 1 to have the proper initialization.

The availability of the forward and backward probabilities suffices for the calculation of the posterior probability:

(5.36)

As stated earlier, the individually most likely state is the one that maximizes p(x(i)|Z(I)). The denominator of Equation (5.36) is not relevant for this maximization since it does not depend on x(i).

The complete forward–backward algorithm is as follows.

Algorithm 5.2 The forward–backward algorithm

1. Perform the forward algorithm as given in Section 5.3.1, resulting in the array F(i, k) with i = 0, …, I and k = 1, …, K.
2. Backward algorithm:

   • Initialization:

     B(I, k) = 1 for k = 1, …, K

   • Recursion:

     for i = I − 1, I − 2, …, 0 and x(i) = 1, …, K

     

3. MAP estimation of the states:
   

The forward–backward algorithm has a computational complexity that is on the order of (I + 1)K². The algorithm is therefore feasible.

5.3.3.2 The Most Likely State Sequence

A criterion that involves the whole sequence of states is the overall uniform cost function. The function is zero when the whole sequence is estimated without any error. It is unity if one or more states are estimated erroneously. Application of this cost function within a Bayesian framework leads to a solution that maximizes the overall posterior probability:

The computation of this most likely state sequence is done efficiently by means of a recursion that proceeds forwards in time. The goal of this recursion is to keep track of the following subsequences:

(5.37)

For each value of x(i), this formulation defines a particular partial sequence. Such a sequence is the most likely partial sequence from time zero and ending at a particular value x(i) at time i given the measurements z(0), …, z(i). Since x(i) can have K different values, there are K partial sequences for each value of i. Instead of using Equation (5.37), we can equivalently use

because p(X(i)|Z(i)) = p(X(i), Z(i))p(Z(i)) and Z(i) is fixed.

In each recursion step the maximal probability of the path ending in x(i) given Z(i) is transformed into the maximal probability of the path ending in x(i + 1) given Z(i + 1). For that purpose, we use the following equality:

Here the Markov condition has been used together with the assumption that the measurements are memoryless.

The maximization of the probability proceeds as follows:

(5.38)

The value of x(i) that maximizes p(x(0), …, x(i), x(i + 1), Z(i + 1)) is a function of x(i + 1):

(5.39)

The so-called Viterbi algorithm uses the recursive equation in (5.38) and the corresponding optimal state dependency expressed in (5.39) to find the optimal path. For that, we define the array

Algorithm 5.3 The Viterbi algorithm

1. Initialization:

   for x(0) = 1, …, K

   • Q(0, x(0)) = p₀(x(0))p_z(z(0)|x(0))
   • R(0, x(0)) = 0

2. Recursion:

   for i = 2, …, I and x(i) = 1, …, K

   • 
   • 
3. Termination:
   • 
   • 

4. Backtracking:

   for i = I − 1, I − 2, …, 0

   • 

The computational structure of the Viterbi algorithm is comparable to that of the forward algorithm. The computational complexity is also on the order of (i + 1)K².

Example 5.9 Offline licence plate detection in videos Figure 5.15 shows the results of the two offline state estimators applied to the video line shown in Figure 5.11. Figure 5.16 provides the results of the whole image shown in Figure 5.10.

Both methods are able to prevent the delay that is inherent in online estimation. Nevertheless, both methods show some falsely detected licence plate pixels on the right side of the plate. These errors are caused by a sticker containing some text. Apparently, the statistical properties of the image of this sticker are similar to the one on a licence plate.

A comparison between the individually estimated states and the jointly estimated states shows that the latter are more coherent and that the former are more fragmented. Clearly, such a fragmentation increases the probability of having erroneous transitions of estimated states. However, generally the resulting erroneous regions are small. The jointly estimated states do not show many of these unwanted transitions. However, if they occur, then they are more serious because they result in a larger erroneous region.

5.3.3.3 Matlab® Functions for HMM

The Matlab^® functions for the analysis of hidden Markov models are found in the Statistics toolbox. There are five functions:

• hmmgenerate: given p_t( · | · ) and p_z( · | · ), generate a sequence of states and observations.
• hmmdecode: given p_t( · | · ), p_z( · | · ) and a sequence of observations, calculate the posterior probabilities of the states.
• hmmestimate: given a sequence of states and observations, estimate p_t( · | · ) and p_z( · | · ).
• hmmtrain: given a sequence of observations, estimate p_t( · | · ) and p_z( · | · ).
• hmmviterbi: given p_t( · | · ), p_z( · | · ) and a sequence of observations, calculate the most likely state sequence.

The function hmmtrain() implements the so-called Baum–Welch algorithm.

5.4.3.1 Matlab® Functions for Particle Filtering

Many Matlab^® users have already implemented particle filters, but no formal toolbox yet exists. Section 9.4 contains a listing of Matlab^® code that implements the condensation algorithm. Details of the implementation are also given.

5.6.1.1 Structuring

The structure of the model is settled by addressing the following questions. What is considered part of the system and what is environment? What are the (controllable) input variables? What are the possible disturbances (process noise)? What are the state variables? What are the output variables? What are the physical laws that relate to the physical variables? Which parameters of these laws are known and which are unknown? Which of these parameters can be measured directly?

Usually, there is not a unique answer to all these questions. In fact, the result of structuring is a set of candidate models.

5.6.1.2 Experiment Design

The purpose of the experimentation is to enable the determination of the unknown parameters and the evaluation and selection of the candidate models. The experiment comes down to the enforcement of some input test signals to the system at hand and the acquisition of measured data. The design aspects of the experiments are the choice of input signals, the choice of the sensors and the sensor locations and the pre-processing of the data.

The input signals should be such that the relevant aspects of the system at hand are sufficiently covered. The bandwidth of the input signal must match the bandwidth of interest. The signal magnitude must be in a range as required by the application. Furthermore, the signal energy should be sufficiently large so as to achieve suitable signal-to-noise ratios at the output of the sensors. Here, a trade-off exists between the bandwidth, the registration interval and the signal magnitude.

For instance, a very short input pulse covers a large bandwidth and permits a low registration interval, but requires a tremendous (perhaps excessive) magnitude in order to have sufficient signal-to-noise ratios. At the other side of the extremes, a single sinusoidal burst with a long duration covers only a very narrow bandwidth, and the signal magnitude can be kept low yet offering enough signal energy. A signal often used is the pseudo-random binary signal.

The choice of sensors is another issue. First, the set of identifiable, unknown parameters of the model must be determined. Sometimes it might be necessary to temporarily use additional (often expensive) sensors so as to enable the estimation of all relevant parameters. As soon as the system identification is satisfactorily accomplished, these additional sensors can be removed from the system.

In order to prevent overfitting, it might be useful to split the data according to two time intervals. The data in the first interval are used for parameter estimation. The second interval is used for model evaluation. Cross-evaluation might also be useful.

5.6.1.3 Parameter Estimation

Suppose that all unknown parameters are gathered in one parameter vector α. The discrete system equation is denoted by f(x, w, α). The measurement system is modelled, as before, by z = h(x, v). We then have the sequence of measurements according to the following recursions:

(5.49)

where I is the length of the sequence, x(0) is the initial condition (which may be known or unknown) and v(i) and w(i) are the measurement noise and the process noise, respectively.

One possibility for estimating α is to process the sequence z(i) in batches. For that purpose, we stack all measurement vectors to one I × n (recall that n is the dimension of z(i)) dimensional vector Z. Equation (5.49) defines the conditional probability density p(Z|α). The stochastic nature of Z is due to the randomness of w(i), v(i) and possibly x(0). Equation (5.49) shows how this randomness propagates to Z. Once the conditional density p(Z|α) has been settled, the complete estimation machinery from Chapter 4 applies, thus providing the optimal solution of α. Especially, maximum likelihood estimation is popular since (5.49) can be used to calculate the (log-) likelihood of α. A numerical optimization procedure must provide the solution.

Working in batches soon becomes complicated due to the (often) non-linear nature of the relations involved in Equation (5.49). Many alternative techniques have been developed. For linear systems, processing in the frequency domain may be advantageous. Another possibility is to process the measurements sequentially. The trick is to regard the parameters as state vectors α(i). Static parameters do not change in time. Therefore, the corresponding state equation is α(i + 1) = α(i). Sometimes, it is useful to allow slow variations in α(i). This is helpful in order to model drift phenomena, but also to improve the convergence properties of the procedure. A simple model would be a process that is similar to a random walk: α(i + 1) = α(i) + ω(i). The white noise sequence ω(i) is the driving force for the changes. Its covariance matrix C_ω should be small in order to prevent a too wild behaviour of α(i).

Using this model, Equation (5.49) is transformed into:

The original state vector has been augmented with α(i). The new state equation can be written as ξ(i + 1) = ϕ(ξ(i), w(i), ω(i)) with . This opens the door to simultaneous online estimation of both x(i) and α(i) using the techniques discussed earlier in this chapter. However, the new state function ϕ( · ) is non-linear and for online estimation we must resort to estimators that can handle these non-linearities, for example extended Kalman filtering (Section 5.2.2) or particle filtering (Section 5.4).

Note that if ω(i) ≡ 0, then α(i) is a random constant and (hopefully) its estimate converges to a constant. If we allow ω(i) to deviate from zero by setting C_ω to some (small) non-zero diagonal matrix, the estimator becomes adaptive. It has the potential to keep track of parameters that drift.

5.6.1.4 Evaluation and Model Selection

The last step is the evaluation of the candidate models and the final selection. For that purpose, we select a quality measure and evaluate and compare the various models. Popular quality measures are the log-likelihood and the sum of squares of the residuals (the difference between measurements and model-predicted measurements).

Models with a low order are preferable because the risk of overfitting is minimized and the estimators are less sensitive to estimation errors of the parameters. Therefore, it might be beneficial to consider not only the model with the best quality measure but also the models that score slightly less, but with a lower order. In order to evaluate these models, other tests are useful. Ideally, the cross correlation between the residuals and the input signal is zero. If not, there is still part of the behaviour of the system that is not explained by the model.

5.6.2.1 Observability

We consider a deterministic linear system, which is a special case of Equations (5.13) and (5.20):

(5.50)

The system is called observable if with known F(i), L(i)u(i) and H(i) the state x(i) (with fixed i) can be solved from a sequence z(i), z(i + 1), … of measurements. The system is called completely observable if it is observable for any i. It can be assumed that L(i)u(i) = 0. This is without any loss of generality since the influence to z(i), z(i + 1), … of an L(i)u(i) not being zero can be neutralized easily. Hence, the observability of a system solely depends on F(i) and H(i).

An approach to find out whether the system is observable is to construct the observability Gramian (Bar-Shalom and Li, 1993). From Equation (5.50):

(5.51)

Equation (5.51) is of the type z = H. The least squares estimate is (see Section 4.3.1). The solution exists if and only if the inverse of H^TH exists, or in other words, if the rank of H^TH is equal to the dimension m of the state vector. Equivalent conditions are that H^TH is positive definite (i.e. y^TH^THy > 0 for every y ≠ 0) or that the eigenvalues of H^TH are all positive (see Appendix B.5).

Translated to the present case, the requirement is that for at least one K ≥ 0 the observability Gramian , defined by

has rank equal to m. Equivalently we check whether the Gramian is positive definite. For time-invariant systems, F and H do not depend on time and the Gramian simplifies to

If the system F is stable (the magnitude of all eigenvalues of F are less than one), we can set K → ∞ to check whether the system is observable.

A second approach to determine the observability of a time-invariant, deterministic system is to construct the observability matrix:

According to Equation (5.51), x(i) can be retrieved from a sequence z(i), …, z(i + m − 1) if is invertible; that is if the rank of equals m.

The advantage of using the observability Gramian instead of the observability matrix is that the former is more stable. Modelling errors and round-off errors in the coefficients in both F and H could make the difference between an invertible or and a non-invertible one. However, is less prone to small errors than .

A more quantitative measure of the observability is obtained by using the eigenvalues of the Gramian . A suitable measure is the ratio between the smallest eigenvalue and the largest eigenvalue. The system is less observable as this ratio tends to zero. A likewise result can be obtained by using the singular values of the matrix (see singular value decomposition in Appendix B.6).

5.6.2.2 Controllability

In control theory, the concept of controllability usually refers to the ability that for any state x(i) at a given time i a finite input sequence u(i), u(i + 1), …, u(i + ℓ − 1) exists that can drive the system to an arbitrary final state x(i + ℓ). If this is possible for any time, the system is called completely controllable. As with observability, the controllability of a system can be revealed by checking the rank of a Gramian. For time-invariant systems, the controllability can also be analysed by means of the controllability matrix. This matrix arises from the following equation:

The minimum number of steps, ℓ, is at most equal to m, the dimension of the state vector. Therefore, in order to test the controllability of the system (F, L) it suffices to check whether the controllability matrix (L FL … F^{m − 1}L) has rank m.

The Matlab^® functions for creating the controllability matrix and Gramian are ctrb() and gram(), respectively.

5.6.2.3 Dynamic Stability and Steady State Solutions

The term stability refers to the ability of a system to resist to and recover from disturbances acting on this system. A state estimator has to face three different causes of instabilities: sensor instability, numerical instability, and dynamic instability.

Apart from the usual sensor noise and sensor linearity errors, a sensor may produce unusual glitches and other errors caused by phenomena hard to predict, such as radio interference, magnetic interference, thermal drift, mechanical shocks, and so on. This kind of behaviour is sometimes denoted by sensor instability. Its early detection can be done using consistency checks (to be discussed later in this section).

Numerical instabilities originate from round-off errors. Particularly, the inversion of a near singular matrix may cause large errors due to its sensitivity to round-off errors. A careful implementation of the design must prevent these errors (see Section 5.6.3).

The third cause for instability lies in the dynamics of the state estimator itself. In order to study the dynamic stability of the state estimator it is necessary to consider the estimator as a dynamic system with as inputs the measurements z(i) and the control vectors u(i) (see Appendix D). The output consists of the estimates . In linear systems, the stability does not depend on the input sequences. For the stability analysis it suffices to assume zero z(i) and u(i). The equations of interest are derived from (5.22) and (5.23):

(5.52)

with

where P(i) is the covariance matrix C(i + 1|i) of the predicted state and P(i) is recursively defined by the discrete Ricatti equation:

(5.53)

The recursion starts with . The first term in Equation (5.53) represents the absorption of uncertainty due to the dynamics of the system during each time step (provided that F is stable; otherwise it represents the growth of uncertainty). The second term represents the additional uncertainty at each time step due to the process noise. The last term represents the reduction of uncertainty thanks to the measurements.

For the stability analysis of a Kalman filter it is of interest to know whether a sequence of process noise, w(i), can influence each element of the state vector independently. The answer to this question is found by writing the covariance matrix C_w(i) of the process noise as C_w(i) = G(i)G^T(i). Here, G(i) can be obtained by an eigenvalue/eigenvector diagonalization of C_w(i), that is C_w(i) = V_w(i)Λ_w(i)V^T_w(i) and G(i) = V_w(i)Λ^1/2_w(i) (see Appendix B.5 and C.3). Λ_w(i) is a K × K diagonal matrix where K is the number of nonzero eigenvalues of C_w(i). The matrices G(i) and V_w(i) are K × K matrices. The introduction of G(i) allows us to write the system equation as:

where is a K-dimensional Gaussian white noise vector with covariance matrix I. The noise sequence w(i) influences each element of x(i) independently if the system (F(i), G(i)) is controllable by the sequence . We then have the following theorem (Bar-Shalom and Li, 1993), which provides a sufficient (but not a necessary) condition for stability:

If a time-invariant system (F, H) is completely observable and the system (F, G) is completely controllable, then for any initial condition P(0) = C_x(0) the solution of the Ricatti equation converges to a unique, finite, invertible steady state covariance matrix P(∞).

The observability of the system ensures that the sequence of measurements contains information about the complete state vector. Therefore, the observability guarantees that the prediction covariance matrix P(∞) is bounded from above. Consequently, the steady state Kalman filter is BIBO stable (see Appendix D.3.2). If one or more state elements are not observable, the Kalman gains for these elements will be zero and the estimation of these elements is purely prediction driven (model based without using measurements). If the system F is unstable for these elements, then so is the Kalman filter.

The controllability with respect to ensures that P(∞) is unique (does not depend on C_x(0)) and that its inverse exists. If the system is not controllable, then P(∞) might depend on the specific choice of C_x(0). If for the non-controllable state elements the system F is stable, then the corresponding eigenvalues of P(∞) are zero. Consequently, the Kalman gains for these elements are zero and the estimation of these elements is again purely prediction driven.

The discrete algebraic Ricatti equation, P(i + 1) = P(i), mentioned earlier in this chapter, provides the steady state solution for the Ricatti equation:

The steady state is reached due to the balance between the growth of uncertainty (the second term and possibly the first term) and the reduction of uncertainty (the third term and possibly the first term).

A steady state solution of the Ricatti equation gives rise to a constant Kalman gain. The corresponding Kalman filter, derived from Equations (5.22), (5.23) and (5.52),

with

becomes time invariant.

5.6.4.1 Orthogonality Properties

We recall that in linear-Gaussian systems, the innovations (residuals) vectors in the update step of the Kalman filter (5.22) are zero mean with the innovation matrix as covariance matrix:

Furthermore, let be the estimation error of the estimate . We then have the following properties:

(5.56)

These properties follow from the principle of orthogonality. In the static case, any unbiased linear MMSE satisfies

The last step follows from K = C_xzC^{− 1}_z (see Equation (4.29)). The principle of orthogonality, , simply states that the covariance between any component of the error and any measurement is zero. Adopting an inner product definition for two random variables e_k and z_l as , the principle can be expressed as e ⊥ z.

Since x(i|j) is an unbiased linear MMSE estimate, it is a linear function of the set . According to the principle of orthogonality, we have e(i|j) ⊥ Z(j). Therefore, e(i|j) must also be orthogonal to any z(l), where l ≤ j, because z(l) is a subspace of Z(j). This proves the first statement in Equation (5.56).

In addition, x(k|l), l ≤ j, is a linear combination of Z(l). Therefore, it is also orthogonal to e(i|j), which proves the second statement.

The whiteness property of the innovations follows from the following argument. Suppose j < i. We may write . In the inner expectation, the measurements Z(j) are known. Since is a linear combination of Z(j), is non-random. It can be taken outside the inner expectation: . However, must be zero because the predicted measurements are unbiased estimates of the true measurements. If , then (unless i = j).

5.6.4.2 Normalized Errors

The NEES (normalized estimation error squared) is a test signal defined as

In the linear-Gaussian case, the NEES has a χ²_m distribution (chi-square with m degrees of freedom; see Appendix C.1.5).

The χ²_m distribution of the NEES follows from the following argument. Since the state estimator is unbiased, E[e(i|i)] = 0. The covariance matrix of e(i|i) is C(i|i). Suppose A(i) is a symmetric matrix such that A(i)A^T(i) = C^{− 1}(i|i). Application of A(i) to e(i) will give a random vector y(i) = A(i)e(i). The covariance matrix of y(i) is A(i)C(i|i)A^T(i) = A(i)(A(i)A^T(i))^{− 1}A^T(i) = I. Thus, the components of y(i) are uncorrelated and have unit variance. If both the process noise and the measurement noise are normally distributed, then so is y(i). Hence, the inner product y^T(i)y(i) = e^T(i)C^{− 1}(i|i)e(i) is the sum of m squared, independent random variables, each normally distributed with zero mean and unit variance. Such a sum has a χ²_m distribution.

The NIS (normalized innovation squared) is a test signal defined as

In the linear-Gaussian case, the NIS has a χ²_n distribution. This follows readily from the same argument as used above.

5.6.4.3 Consistency Checks

From Equation (5.56) and the properties of the NEES and the NIS, the following statement holds true.

If a state estimator for a linear-Gaussian system is optimal, then:

• The sequence Nees(i) must be χ²_m distributed.
• The sequence Nis(i) must be χ²_n distributed.
• The sequence (innovations) must be white.

Consistency checks of a ‘state estimator-under-test’ can be performed by collecting the three sequences and by applying statistical tests to see whether the three conditions are fulfilled. If one or more of these conditions are not satisfied, then we may conclude that the estimator is not optimal.

In a real design, the NEES test is only possible if the true states are known. As mentioned above, such is the case when the system is (temporarily) provided with extra measurement equipment that measures the states directly and with sufficient accuracy so that their outputs can be used as a reference. Usually such measurement devices are too expensive and the designer has to rely on the other two tests. However, the NEES test is applicable in simulations of the physical process. It can be used to see if the design is very sensitive to changes in the parameters of the model. The other two tests use the innovations. Since the innovations are always available, even in a real design, these are of more practical significance.

In order to check the hypothesis that a dataset obeys a given distribution, we have to apply a distribution test. Examples are the Kolmogorov–Smirnov test and the chi-square test (Section 9.4.3). Often, instead of these rigorous tests a quick procedure suffices to give an indication. We simply determine an interval in which, say, 95% of the samples must be. Such an interval [A, B] is defined by F(B) − F(A) = 0.95. Here, F( · ) is the probability distribution of the random variable under test. If an appreciable part of the dataset is outside this interval, then the hypothesis is rejected. For chi-square distributions there is the one-sided 95% acceptance boundary with A = 0 and B such that . Sometimes the two-sided 95% acceptance boundaries are used, defined by and . Table 5.2 gives values of A and B for various degrees of freedom (Dof).

For state estimators that have reached the steady state, S(i) is constant and the whiteness property of the innovations implies that the power spectrum of any element of must be flat. Suppose that is an element of . The variance of , denoted by σ²_l, is the lth diagonal element in S. The discrete Fourier transform of , calculated over i = 0, …, I − 1, is

The periodogram of , defined here as , is an estimate of the power spectrum (Blackman and Tukey, 1958). If the power spectrum is flat, then E[P_l(k)] = σ²_l. It can be proven that in that case the variables 2P_l(k)/σ²_l is χ²₂ distributed for all k (except for k = 0). Hence, the whiteness of is tested by checking whether 2P_l(k)/σ²_l is χ²₂ distributed.

5.6.4.4 Fudging

If one or more of the consistency checks fail, then somewhere a serious modelling error has occurred. The designer has to step back to an earlier stage of the design in order to identify the fault. The problem can be caused anywhere, from inaccurate modelling during the system identification to improper implementations. If the system is non-linear and the extended Kalman filter is applied, problems may arise due to the neglect of higher-order terms of the Taylor series expansion.

A heuristic method to catch the errors that arise due to approximations of the model is to deliberately increase the modelled process noise (Bar-Shalom and Li, 1993). One way to do so is by increasing the covariance matrix C_w of the process noise by means of a fudge factor γ. For instance, we simply replace C_w by C_w + γI. Instead of adapting C_w we can also increase the diagonal of the prediction covariance C(i + 1|i) by some factor. The fudge factor should be selected such that the consistency checks now pass as successfully as possible.

Fudging causes a decrease of faith in the model of the process and thus causes a larger impact of the measurements. However, modelling errors give rise to deviations that are autocorrelated and not independent from the states. Thus, these deviations are not accurately modelled by white noise. The designer should maintain a critical attitude with respect to fudging.