4
$\begingroup$

I struggle to understand the Kalman filter and the documentation for KalmanEstimator. Please help me understand the first Application on that documentation page. The code in question is copied here:

antenna=StateSpaceModel[{{{0.5, 0.07869}, {0, -0.60653}}, {{0.0042, 0.0104}, {0.0786, 0.00786}}, {{1, 0}}, {{0, 0}}}, SamplingPeriod -> 0.1, SystemsModelLabels -> None];
{w,v}={{{0.01}},{{0.001}}};
kalmanFilter=SystemsModelExtract[KalmanEstimator[{antenna,All,1},{w,v}],All,{3}];
u= Table[Sin[(2*Pi*i)/20.0],{i,100}];
processNoise=RandomReal[NormalDistribution[0,Sqrt[w[[1,1]]]],{100}];
measurementNoise=RandomReal[NormalDistribution[0,Sqrt[v[[1,1]]]],{100}];
y=Flatten[OutputResponse[antenna,{u,processNoise}]]+measurementNoise;
ListLinePlot[OutputResponse[kalmanFilter,{u,y}]]

kalmanFilter output

In the last line above we provided OutputResponse with u and the plot above looks very much like the plot below.

ListLinePlot[0.015*u]

plot of u

So, all the above kalmanFilter did is give a slightly modified version of u, and that is not very impressive. Given how famous the KalmanFilter is I am expecting it to do something more impressive. What am I missing? Also what does process noise represent and what does measurement noise represent? (** UPDATE **) After watching the videos Daniel Huber recommended, I think I understand process noise and measurement noise. Process noise is errors introduced due to the mathematical approach to model the problem. Process noise is often from factors that we ignore to simplify the model (e.g. friction, wind-resistnace, assuming a flat earth over short distances, ...). Measurement noise is errors due the limitations of the sensor(s) making measurements.

$\endgroup$
1
  • 2
    $\begingroup$ Kalman filter is not easy to understand, it needs some effort. See e.g. the videos of Michel van Biezen (there are a lot of them.). The first you find here: youtube.com/… $\endgroup$ Commented Mar 25, 2022 at 9:47

1 Answer 1

6
$\begingroup$

First consider the response of the antenna to a sinusoid. The response is attenuated and has a delay. (This can also be seen from the BodePlot, if need be.)

antenna = StateSpaceModel[{{{0.5, 0.07869}, {0, -0.60653}}, {{0.0042, 0.0104}, 
 {0.0786, 0.00786}}, {{1, 0}}, {{0, 0}}}, SamplingPeriod -> 0.1, SystemsModelLabels -> None];
u = Table[Sin[(2*Pi*i)/20.0], {i, 100}];

y = OutputResponse[antenna, {u, ConstantArray[0, 100]}][[1]];
ListLinePlot[{0.015 u, %}, PlotLegends -> {"scaled input", "output"}]

enter image description here

Now introduce noise both to the process and the measurement, and as expected the output is also noisy - actually a noisy sinusoid.

{w, v} = {{{0.01}}, {{0.001}}}; 
processNoise = RandomReal[NormalDistribution[0, Sqrt[w[[1, 1]]]], {100}];
measurementNoise = RandomReal[NormalDistribution[0, Sqrt[v[[1, 1]]]], {100}];
ynoisy = Flatten[OutputResponse[antenna, {u, processNoise}]] + measurementNoise;
ListLinePlot[{y, ynoisy}, PlotLegends -> {"output", "noisy output"}]

enter image description here

What KalmanEstimator does is to clean up the noisy output using the input that produced that noisy output and the noisy output itself. The mistake you are making is comparing u and the yestim. If you compare ynoisy and yestim, the estimator has done a really good job.

kalmanEstim = SystemsModelExtract[KalmanEstimator[{antenna, All, 1}, {w, v}], All, {3}];
yestim = OutputResponse[kalmanEstim, {u, ynoisy}][[1]];
ListLinePlot[{ynoisy, yestim}, PlotLegends -> {"noisy output", "estimated output"}]

enter image description here

If you compare y and yestim, the estimator seems to have done too good of a job!

ListLinePlot[{y, yestim}]

enter image description here

In reality the model and covariances do not perfectly model the physical system and actual noise. For example, if we vary the covariances, we can see something closer to be the actual performance of the estimator.

processNoise = RandomReal[NormalDistribution[0, Sqrt[w[[1, 1]]] + 1], {100}];
measurementNoise = 
  RandomReal[NormalDistribution[0, Sqrt[v[[1, 1]]] + 1], {100}];
Flatten[OutputResponse[antenna, {u, processNoise}]] + measurementNoise;
OutputResponse[kalmanEstim, {u, %}][[1]];
ListLinePlot[{y, %}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.