6
$\begingroup$

The first "Application" in the documentation for KalmanEstimator is copied below with additional variables for clarity.

Clear["Global`*"];
a={{0.5,0.07869},{0,-0.60653}};
b={{0.0042,0.0104},{0.0786,0.00786}};
c={{1,0}};
d={{0,0}};
antenna=StateSpaceModel[{a,b,c,d},SamplingPeriod->0.1];
processVariance=0.01;
measurementVariance=0.001;
processCovariance={{processVariance}};
measurementCovariance={{measurementVariance}};
kalmanEstimate=KalmanEstimator[{antenna,All,1},
   {processCovariance,measurementCovariance}];
kalmanFilter=SystemsModelExtract[kalmanEstimate,All,{3}];

The documentation uses (u, y) where I use (inputSignal, noisySignal) respectively.

inputSignal= Table[Sin[0.1*Pi*i],{i,100}];
processNoise=RandomReal[NormalDistribution[0,Sqrt[processVariance]],{100}];
measurementNoise=RandomReal[
   NormalDistribution[0,Sqrt[measurementVariance]],{100}];
noisySignal=Flatten[OutputResponse[antenna,{inputSignal,processNoise}]]+
   measurementNoise;

A Kalman filter should take (noisySignal), and give an estimate of (inputSignal), and I suppose the next line is supposed to do that for us. However, in an actual application I have (noisySignal), but I don't know what (inputSignal) is. If I know what (inputSignal) is, I don't need to estimate it! Please explain where I am confused, or show how KalmanEstimator in Mathematica can be used to implement a Kalman filter.

estimatedSignal=OutputResponse[kalmanFilter,{inputSignal,noisySignal}];

I include some plots below.

ListLinePlot[noisySignal]

ListLinePlot[{inputSignal, 67.41*First@estimatedSignal}]

$\endgroup$
  • $\begingroup$ KalmanEstimator is estimating noisySignal not inputSignal. For your application have you got the state-space representation, and identified which inputs/outputs are noisy? Is KalmanFilter more suited for your purpose? $\endgroup$ – Suba Thomas Sep 2 '15 at 19:59
  • $\begingroup$ Suba Thomas, Really? I added plots at the end of the original question. You will see that noisySignal is hardly like estimated Signal. One the other hand 67.41*estimatedSignal is close to a delayed version of inputSignal. $\endgroup$ – Ted Ersek Sep 2 '15 at 22:37
  • $\begingroup$ Yes, really, and it can be readily verified. Set the process and measurement noise to zero and get the actual output signal. Now compare actual and estimatedSignal and see the performance of the estimator. $\endgroup$ – Suba Thomas Sep 3 '15 at 13:56
  • $\begingroup$ It is very evident from your edit that the estimator is not estimating the input signal. If it was, it is doing a horribly poor job because it has delay and is off by a factor of 67.41! $\endgroup$ – Suba Thomas Sep 3 '15 at 13:58
6
$\begingroup$

The short answer is yes, we need to know the input signal when using the Kalman filter.

Let's consider an example. Let's say we have a robot and we want to estimate it's position as it drives around in the real world. We can change the location of the robot by applying a torque we can control on its wheels. In theory, we can calculate the location of the robot by measuring the torque we apply and figuring out how many rounds the wheels have turned, which can then be mapped into a location. However, there are all sorts of uncertainties in the calculation. For example, the wheel may slip, and the location may change the by the wind, etc. Luckily, the robot is equipped with a GPS sensor which can measure its location in the world. However, the GPS has its own error in measuring the location of the robot. A better way of estimating the position of the robot is to combine the measurements from the GPS system and the measurements of the torque on the wheels, using the Kalman filter.

So in this example, the state of the system is the location of the robot. We know a model of the system, which is the mapping between the torque on the wheels and the location of the robot. But we know that this model alone is not an accurate description of the state of the system since there are noises in this process. And our aim is to estimate the state of the system, using the input signal (the torque), the model of the system(mapping between torque and location), and the measured noised signal (GPS measurements).

Now back to the example you provide (I've changed the notation slightly for clarity).

Mathematica graphics

We have an input signal of the system u that we know, and we have a model for the system that maps the input signal to the output signal y. And there are two kinds of noises in this process. There is the process noise (because our model doesn't capture all the perturbations in the world), and the measurement noise (because the measurement instrument is not accurate). And our aim is to estimate the real signal from the system, output from the response to our input together with all the perturbations. For example, in the figure above, u is the input we know. The perturbation that we don't know can be thought as another input to the system. We can denote the deterministic input signal we know as $u_d$ and the perturbation we don't know as $u_w$. The direct output from the system under $u_d + u_w$ is y, which is the output we want to estimate. But due to the measurement noise v, what we can get out from the system is $y_s = y + v$. We can use the Kalman estimator to get an estimate of the state y as $\hat{y}_s$. We use the deterministic input $u_d$ and the noised output $y_s$ as the input to the Kalman estimator, and the output of the estimator is $\hat{y}_s$.

Mathematica graphics

This is the model of our system.

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 and the v are the covariances of the process and measurement noise. They describe how large our noise can vary.

{w, v} = {{{0.01}}, {{0.001}}};

This describes how the antenna system is connected to the Kalman estimator

kal = KalmanEstimator[{antenna, 1, 1}, {w, v}];

This is our deterministic input ($u_d$ in the figure)

u = Table[Sin[(2 π i)/20.0], {i, 100}];

and the process and measurement noise

processNoise = RandomReal[NormalDistribution[0, Sqrt[w[[1, 1]]]], {100}];
measurementNoise = RandomReal[NormalDistribution[0, Sqrt[v[[1, 1]]]], {100}];

This is the output from the system, under the deterministic input $u_d$ and a noise input $u_w$. This signal is also what we want to estimate.

y = Flatten[OutputResponse[antenna, {u, processNoise}]];

Due to the measurement noise, we don't have access to y. Instead, we only have access to a noised output $y_s$

ys = ytrue + measurementNoise;

We can get an estimate of y as $\hat y_s$

yEstimate = Last@OutputResponse[{kal}, {u, ys}];

We can see that it is a good estimation, compared to the true output y

ListPlot[{y, yEstimate}, Joined -> True, 
 PlotLegends -> Placed[{"true y", "estimated y"}, Top]]

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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