Open mfbolus opened 3 years ago
The steady state matlab Kalman gain identified by n4sid
is often apparently the wrong sign.
Besides this larger issue, the Kalman gain for the 'innovations form' of the Kalman filter is slightly different from the normal Kalman gain:
K_matlab = A*K
This particular change has been accounted for already. But the apparent sign mismatch is persistent. Matlab defines the error as (z - yHat)
, so it is consistent with this implementation. To make things more confusing, if n4sid
is used to fit 2nd order model with CVA weighting on ED example L6CT dataset, the identified Kalman gain appears fine after accounting for A
.
It is not just the wrong sign. The forward innovations form conversion to the filter update form above (K_innovations = A*K_filter
) assumes that the process and measurement noise are independent, i.e., that S=0
where S
is cross-cov of the state-output residuals.
From Katayama 2005, p. 120:
Kinn_t = (A*P_t*C' + S) * inv(C*P_t*C' + R)
Kfilt_t = (P_t*C') * inv(C*P_t*C' + R)
I did some testing and found that if I used van Overschee's matlab toolbox for SSID that accompanies van Overschee et de Moor 1996, but manually zero out this cross term, the Kalman gain is consistent with Kinn = A * Kfilt
(though not exact for whatever reason). Upon inspection, the van Overschee identified cross term S
is non-negligibly large. I just ignore this term when I fit models and do Kalman filtering because I assume the noises to be independent. This is a common assumption.
I now need to circle back to matlab's n4sid results but there is actually not an option to force S
to zeros and the function does not return S, so we have no way of back-calculating the filter update form of the Kalman gain. It's becoming pretty clear at this point that I either need to write the library to do control using predicted state (x_(t+1|t)
) rather than the filtered estimate of the current state (x_(t|t)
), or we need everyone to use my GLDS ssid fit code which is not ideal.
I have used Matlab's kalman
function to design filter gains, given Q
and R
(for how-to see this note). In Matlab's language, what we are looking for they call M_x
. The innovation form gain they call L
. They also often call L
K
, unless I am wrong, which is confusing.
Using this matlab function with Q
and R
identified from either the van Overschee toolbox or my fit code, I get consistent answers for the two forms of the Kalman gains, so it appears it is not a strange inconsistency between matlab gain definitions and other literature. I would not have thought so, but wanted to check.
What remains, however, is that the Kalman gains fit by ssest
or n4sid
matlab functions are usually the wrong sign at the very least. I will have to leave this issue open for now.
This is an issue carried forward from when this was hosted on the github.gatech.edu domain. In the interest of resolving it, I will past my comments below.