Design Case Studies

Steady-State Design

You can design the steady-state Kalman filter described above with the function kalman. First specify the plant model with the process noise.

This is done by

% Note: set sample time to -1 to mark model as discrete
Plant = ss(A,[B B],C,0,-1,'inputname',{'u' 'w'},...
                          'outputname','y');

Assuming that , you can now design the discrete Kalman filter by

Q = 1; R = 1;
[kalmf,L,P,M] = kalman(Plant,Q,R);

This returns a state-space model kalmf of the filter as well as the innovation gain

M

M =
   3.7980e-01
   8.1732e-02
  -2.5704e-01

The inputs of kalmf are and , and its outputs are the plant output and state estimates and .

Because you are interested in the output estimate , keep only the first output of kalmf. Type

kalmf = kalmf(1,:);
kalmf
a = 
                      x1_e         x2_e         x3_e
         x1_e       0.7683       -0.494       0.1129
         x2_e       0.6202            0            0
         x3_e    -0.081732            1            0
 
 
b = 
                         u            y
         x1_e      -0.3832       0.3586
         x2_e       0.5919       0.3798
         x3_e       0.5191     0.081732
 
 
c = 
                      x1_e         x2_e         x3_e
          y_e       0.6202            0            0
 
 
d = 
                         u            y
          y_e            0       0.3798
 
I/O groups:                            
      Group name      I/O    Channel(s)
      KnownInput       I         1     
     Measurement       I         2     
    OutputEstimate     O         1     
                                       
Sampling time: unspecified
Discrete-time model.

To see how the filter works, generate some input data and random noise and compare the filtered response with the true response . You can either generate each response separately, or generate both together. To simulate each response separately, use lsim with the plant alone first, and then with the plant and filter hooked up together. The joint simulation alternative is detailed next.

The block diagram below shows how to generate both true and filtered outputs.

You can construct a state-space model of this block diagram with the functions parallel and feedback. First build a complete plant model with as inputs and and (measurements) as outputs.

a = A;
b = [B B 0*B];
c = [C;C];
d = [0 0 0;0 0 1];
P = ss(a,b,c,d,-1,'inputname',{'u' 'w' 'v'},...
                      'outputname',{'y' 'yv'});

Then use parallel to form the following parallel connection.

sys = parallel(P,kalmf,1,1,[],[])

Finally, close the sensor loop by connecting the plant output to the filter input with positive feedback.

% Close loop around input #4 and output #2
SimModel = feedback(sys,1,4,2,1)
% Delete yv from I/O list
SimModel = SimModel([1 3],[1 2 3])

The resulting simulation model has as inputs and as outputs.

SimModel.inputname

ans = 
    'w'
    'v'
    'u'

SimModel.outputname

ans = 
    'y'
    'y_e'

You are now ready to simulate the filter behavior. Generate a sinusoidal input and process and measurement noise vectors and .

t = [0:100]';
u = sin(t/5);

n = length(t)
randn('seed',0)
w = sqrt(Q)*randn(n,1);
v = sqrt(R)*randn(n,1);

Now simulate with lsim.

[out,x] = lsim(SimModel,[w,v,u]);

y = out(:,1);   % true response
ye = out(:,2);  % filtered response
yv = y + v;     % measured response

and compare the true and filtered responses graphically.

subplot(211), plot(t,y,'--',t,ye,'-'), 
xlabel('No. of samples'), ylabel('Output')
title('Kalman filter response')
subplot(212), plot(t,y-yv,'-.',t,y-ye,'-'),
xlabel('No. of samples'), ylabel('Error')

The first plot shows the true response (dashed line) and the filtered output (solid line). The second plot compares the measurement error (dash-dot) with the estimation error (solid). This plot shows that the noise level has been significantly reduced. This is confirmed by the following error covariance computations.

MeasErr = y-yv;
MeasErrCov = sum(MeasErr.*MeasErr)/length(MeasErr);
EstErr = y-ye;
EstErrCov = sum(EstErr.*EstErr)/length(EstErr);

The error covariance before filtering (measurement error) is

MeasErrCov
 
MeasErrCov =
    1.1138

while the error covariance after filtering (estimation error) is only

EstErrCov
 
EstErrCov =
    0.2722

Discrete Kalman Filter

Time-Varying Kalman Filter