I am looking for an example (or some guidance) on how to fit a Kalman smoother to data using log-likelihood maximization using LowLevelParticleFilters.jl.
I already have a model structure and have collected a dataset. It is a fairly simple 3-state LTI system, of which I can measure 2. It models the (very slow) thermal dynamics of a building in the presence of several large disturbances. I have decently accurate data on a few of these, like the ambient temperature and the solar irradiance. But others I am still missing, like people and heat gains from electronics. (it is actually very similar to the content described in the tutorials by Baggepinnen)
Now my goals are as follows:
- Calibrate the model parameters for 24 hour simulation performance. I will use this in an MPC context.
- Calibrate the process and measurement noise covariance matrices.
I augmented the state-space model with an integrating disturbance, I tried to follow: Disturbance gallery · LowLevelParticleFilters Documentation
I already followed the tutorials on youtube by Baggepinnen. I also read through the documentation, and followed the tutorials on JuliaHub.
First, I understand I need to disretize the linear model. As described in How to Tune a Kalman Filter: Step-by-Step Guide | JuliaHub, we can use the matrix exponential for this. I implemented the discrete dynamics as follows:
function setup_filter(p::Vector{T}) where {T}
# unpack parameters
C_i, C_e, C_h, R_ia, R_ie, R_ea, R_ih, A_e, A_w = p[1:9]
σi, σe, σh, σq = p[10:13] # process noise std dev
σmi, σmh = p[14:15] # measurement noise std dev (interior thermostat, supply temp)
Ti0, Te0, Th0 = p[16:18] # initial state
# continuous-time state-space matrices (3x3 system)
A_cont = @SMatrix [
(-1/R_ia-1/R_ie-1/R_ih)/C_i 1/(C_i*R_ie) 1/(C_i*R_ih);
1/(C_e*R_ie) (-1/R_ea-1/R_ie)/C_e 0;
1/(C_h*R_ih) 0 -1/(C_h*R_ih)
]
B_cont = @SMatrix [
1/(C_i*R_ia) 0 A_w/C_i;
1/(C_e*R_ea) 0 A_e/C_e;
0 1/C_h 0
]
G = @SMatrix [1 / C_i; 0; 0]
# x_aug = [Ti, Te, Th, q]' where q is the integrating disturbance
A_cont_aug = @SMatrix [
A_cont[1, 1] A_cont[1, 2] A_cont[1, 3] G[1];
A_cont[2, 1] A_cont[2, 2] A_cont[2, 3] G[2];
A_cont[3, 1] A_cont[3, 2] A_cont[3, 3] G[3];
0.0 0.0 0.0 0.0
]
B_cont_aug = @SMatrix [
B_cont[1, 1] B_cont[1, 2] B_cont[1, 3];
B_cont[2, 1] B_cont[2, 2] B_cont[2, 3];
B_cont[3, 1] B_cont[3, 2] B_cont[3, 3];
0.0 0.0 0.0
]
# augmented AB matrix for discretization
AB_cont = zeros(7, 7)
AB_cont[1:4, 1:4] = Matrix(A_cont_aug)
AB_cont[1:4, 5:7] = Matrix(B_cont_aug)
# discretize
AB_disc = exp(AB_cont * Ts)
A_disc = SMatrix{4,4}(AB_disc[1:4, 1:4])
B_disc = SMatrix{4,3}(AB_disc[1:4, 5:7])
# y = C * x + 0 * u --> output is [Ti, Th]
C = @SMatrix [1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0]
# process noise covariance R1
R1 = Diagonal([σi^2 / C_i^2, σe^2 / C_e^2, σh^2 / C_h^2, σq^2 / C_i^2])
# measurement noise covariance R2
R2 = Diagonal([σmi^2, σmh^2])
# initial state distribution d0
d0 = MvNormal(T.([Ti0, Te0, Th0, 0.0]), Diagonal([1.0, 5.0, 5.0, 1.0]))
return KalmanFilter(
A_disc, B_disc, C, zeros(ny, nu),
R1, R2, d0; Ts, p
)
end
# Re-initialize with corrected filter
p = [p_init[Param] for Param in
["C_i", "C_e", "C_h", "R_ia", "R_ie", "R_ea", "R_ih", "A_e", "A_w",
"σi", "σe", "σh", "σq", "σmi", "σmh", "Ti0", "Te0", "Th0"]]
kf = setup_filter(p);
measure = LowLevelParticleFilters.measurement(kf)
The problem here is that it doesn’t seem to be possible to autodiff through this, because of the matrix exponential. When I replace the matrix exponential by a simple forward euler step (not ideal for accuracy), I still get autodiff complaints that KalmanFilter constructor in LowLevelParticleFilters.jl is internally calling eigvals() which then fails.
Should I use the c2d(ss) as described here instead? Will that work with autodiff?
Discretization · LowLevelParticleFilters Documentation
I would appreciate some guidance on how to approach this problem.
