Developing a Kalman filter algorithm to build my experience with them before working on avionics for spaceshot this spring.
As I go through the process of developing this algo, I am following this tutorial.
This type of filter is used in the tutorial as a lead-up to the Kalman filter. While the Kalman filter treats each variable as random and uses standard distributions to create an estimate, the
The filter makes use of two important equations for each variable (position, velocity, and acceleration):
The point of these equations are to update our estimate of the true quantity considering a new measurement and that time has passed. The coefficients can be changed based on sensor accuracy.
For POSITION:
For VELOCITY:
For ACCELERATION:
-
$example_{est}$ represents the Kalman filter's guess on the state. -
$example_{pred}$ is given by the state prediction made in the previous time step. -
$\alpha$ ,$\beta$ , and$\gamma$ are the gains. If the measurement accuracy is high, alpha/beta should approach 1, as this will make the measurement influence the equation more. If the measurement accuracy is low, alpha/beta should approach 0 for the same reason. -
$measurement_n$ is the position measurement taken by the sensor at time n. In my code, this is represented by a time series array.
This equation is given by the system-- in this case, we are using 1d kinematics. This gives us the equation:
where
This filter is very useful if the system is relatively simple, and the equations governing the system dynamics are known. However, most real life systems have more unpredictability and extremely complicated equations, which is why the Kalman filter is widely used.
This filter uses the state update equations but also accounts for the variance in the various measurements. This allows the algorithm to converge close to the true value.
It does this by not only predicting the various system states (position, velocity, acceleration etc.), but also the VARIANCE of those system states. In other words, it uses the measured states and respective variances, as well as the estimated states and variances.
There are two types of noise that can affect measurements that the Kalman filter can account for.
- Measurement Noise: This is noise caused by the error in the measurement due to sensor inaccuracy.
- Process Noise: This is noise caused by external impacts that change the physical state of the system. For example, if we are measuring the position of a plane, wind might create process noise.
We still follow the general algorithm:
- Initialize
- Measure
- Update
- Predict
and then loop over 1, 2, and 3 for all the data points in the set.
We get the following 5 equations:
- State Equation:
(where
- Covariance Update Equation:
- Kalman Gain Equation:
(where
-
State Extrapolation Equation: This equation updates the estimate for the variable according to the system dynamics. In the examples here this is either constant or in accordance with simple kinematics.
-
Covariance Extrapolation Equation: This equation does the same thing as the State Extrapolation Equation, but instead of updating the estimate for the state, it updates the estimate for the covariance of the estimated state. With process noise added, whatever covariance is of the highest degree/constant should be += with process noise covariance.
If we use the constant dynamic system to measure something with a constant increase, for example the temperature of a pot of water heating up linearly using a small process noise covariance, we will get a lag error.
- This means that while the Kalman Filter will eventually converge on the correct slope, the actual values will lag behind the true values by a potentially significant amount.
We can fix this, even if our system dynamics are not well defined. Increasing the process noise variance is the key.
- Once we do this, the Kalman filter is able to accurately follow and estimate the true temperature of the water. YAY!
Many systems take place in more than one dimension. In order to deal with this, we can use matrix formulas and linear algebra. As a result, this gets kind of mathy.
We can represent the space state with two equations:
-
$\mathbf{\dot x(t)} = \mathbf{Ax(t)} + \mathbf{Bu(t)}$ -
$\mathbf{y(t)} = \mathbf{Cx(t)} + \mathbf{Du(t)}$
in which $\mathbf{\dot x(t)} = $ the derivative of the state vector, $\mathbf{y(t)} = $ the output vector,
In higher order systems, we can solve for
$\mathbf{\hat x_{n+1, n}} = \mathbf{F \hat x_{n,n}} + \mathbf{G \hat u_{n,n}} + \mathbf{w_n}$
in which the hat represents estimation,
In order to find
$x_{n+1} = x_n e^{\mathbf{A} \Delta t}$
If we extend this to the whole state vector, it makes sense that
$\mathbf{F} = e^{\mathbf{A} \Delta t}$
Taking the exponential of a matrix can be done with a Taylor Polynomial. If we replace
$e^{\mathbf{X}} = \sum_{k=0}^{\infty} \frac{1}{k!} \mathbf{X}^k$
which is kind of neat. In most cases,
For zero-order hold sampling (which is what we have in a real life rocket scenario-- sensors sample at 10 Hz on our avionics board I believe), we have a consistent
$x_{n+1, n} = x(t + \Delta t)$
which means that
$\mathbf{x(t + \Delta t)} = e^{\mathbf{A} \Delta t} \mathbf{x(t)} + \int_{0}^{\Delta t} e^{\mathbf{A}t} dt \mathbf{B} \mathbf{u(t)}$
We notice that this equation looks similar in form to (3.), giving us
$\mathbf{F} = e^{\mathbf{A} \Delta t}$
and
$\mathbf{G} = \int_{0}^{\Delta t} e^{\mathbf{A}t} dt \mathbf{B}$
As a result, once we solve for
The general form of this equation is:
$\mathbf{P_{n+1,n}} = \mathbf{FP_{n,n}} \mathbf{F}^T + \mathbf{Q}$
where $\mathbf{P}{n,n}$ is the covariance matrix of the state estimate, $\mathbf{P}{n+1,n}$ is the covariance matrix of the state prediction,
Generalize measurement equation:
where
All of the multivariate covariance equations are in the form
$E(\mathbf{ee^T})$
where E denotes the expected value, and
This then leads to:
$\mathbf{R_n} = E(\mathbf{v_n v_n^T})$
where
$\mathbf{Q_n} = E(\mathbf{w_n w_n^T})$
where
- $\mathbf{P_{n,n}} = E(\mathbf{e_n e_n^T}) = E((\mathbf{x_n - \hat{x}{n,n}})(\mathbf{x_n - \hat{x}{n,n}}))$
$\mathbf{\hat x_{n,n} = \hat x_{n,n-1} + K_n (z_n - H \hat x_{n,n-1})}$
where the
Sources: michaelscheinfeild, Kalman Filter from the Ground Up: Alex Becker,