|
| 1 | +## Kalman Filtering |
| 2 | +# Goal: Estimate the true state of a system when measurements are noisy. |
| 3 | +# |
| 4 | +# In many practical systems, we cannot directly observe the full state x_t. |
| 5 | +# Instead, we have noisy measurements y_t. The Kalman filter provides the |
| 6 | +# optimal linear estimate of the hidden state, assuming linear dynamics |
| 7 | +# and Gaussian noise. |
| 8 | + |
| 9 | +using LinearAlgebra |
| 10 | + |
| 11 | +# The Kalman filter operates in two steps: prediction and update. |
| 12 | + |
| 13 | +## Prediction Step: |
| 14 | +# Forecast the next state based on the previous estimate and the system model. |
| 15 | +# x_pred = A * x_prev + B * u_prev |
| 16 | +# - This uses the known dynamics of the system (A, B) and the previous control input u_prev |
| 17 | +# P_pred = A * P_prev * A' + Q |
| 18 | +# - The uncertainty (covariance) is propagated forward |
| 19 | +# - Q represents process noise: uncertainty in how the system evolves |
| 20 | +# - Larger P_pred means we are less confident in the prediction |
| 21 | + |
| 22 | +## Update Step: |
| 23 | +# Correct the prediction using the latest measurement y |
| 24 | +# K = P_pred * C' * inv(C * P_pred * C' + R) |
| 25 | +# - The Kalman gain K determines how much we trust the measurement vs. prediction |
| 26 | +# x_hat = x_pred + K * (y - C * x_pred) |
| 27 | +# - Innovation (y - C*x_pred) adjusts our estimate based on new information |
| 28 | +# P = (I - K * C) * P_pred |
| 29 | +# - Covariance decreases after incorporating measurement, increasing confidence |
| 30 | + |
| 31 | +function kalman_filter(A, B, C, Q, R, y, u=nothing) |
| 32 | + n = size(A,1) |
| 33 | + N = length(y) |
| 34 | + x_hat = zeros(n, N) # estimated states |
| 35 | + P = zeros(n,n) # covariance |
| 36 | + I = Matrix(I, n, n) |
| 37 | + |
| 38 | + # Handle optional control input |
| 39 | + u = isnothing(u) ? zeros(size(B,2), N) : u |
| 40 | + |
| 41 | + for k in 1:N |
| 42 | + # Prediction |
| 43 | + if k == 1 |
| 44 | + x_pred = zeros(n) # initial guess for first step |
| 45 | + P_pred = P + Q |
| 46 | + else |
| 47 | + x_pred = A * x_hat[:,k-1] + B * u[:,k-1] |
| 48 | + P_pred = A * P * A' + Q |
| 49 | + end |
| 50 | + |
| 51 | + # Update |
| 52 | + K = P_pred * C' * inv(C * P_pred * C' + R) |
| 53 | + x_hat[:,k] = x_pred + K * (y[k] - C * x_pred) |
| 54 | + P = (I - K * C) * P_pred |
| 55 | + end |
| 56 | + |
| 57 | + return x_hat |
| 58 | +end |
| 59 | + |
| 60 | +## Example usage: |
| 61 | +# A simple demonstration with a small system. This shows how the filter |
| 62 | +# reconstructs the true state from noisy measurements. |
| 63 | + |
| 64 | +if abspath(PROGRAM_FILE) == @__FILE__ |
| 65 | + # Define system matrices (example: 2-state thermal system) |
| 66 | + A = [1.0 1.0; 0.0 1.0] # state transition |
| 67 | + B = [0.0; 1.0] # control input |
| 68 | + C = [1.0 0.0] # measurement matrix |
| 69 | + Q = 0.01 * I(2) # process noise covariance |
| 70 | + R = 0.1 # measurement noise covariance |
| 71 | + |
| 72 | + # Simulated noisy measurements |
| 73 | + true_states = [20.0, 21.0, 23.0, 22.0] |
| 74 | + y = true_states .+ randn(length(true_states)) * sqrt(R) |
| 75 | + |
| 76 | + # Run Kalman filter |
| 77 | + x_hat = kalman_filter(A, B, C, Q, R, y) |
| 78 | + |
| 79 | + println("Noisy measurements: ", y) |
| 80 | + println("Estimated states: ", x_hat) |
| 81 | +end |
0 commit comments