rlsests.m

本书是电子通信类的本科、研究生辅助教材

function  [sys, x0, str, ts]  = rlsests(t,x,u,flag,nstates,lambda,dt)
%RLSESTS S-Function to perform system identification.
%
%       This M-file is designed to be used in a SIMULINK S-function block.
%       It performs parameter estimation using the Recursive Least Squares
%	Parameter Estimation Algorithm with Exponential Data Weighting
%	
%       The input arguments are
%
%		nstates:        the number of states in the states vector
%		lambda:		the exponential data weighting factor
%       	dt:             how often to sample points (secs)
%
%	The RLS estimator is defined by the following equations:
%
%		                          1      P(k-2) * phi(k-1) * [y(k) - phi(k-1)'theta(k-1)]
%		theta[k] = theta[k-1] + ------ * -------------------------------------------------
%		                        lambda         lambda + phi(k-1)' * P(k-2) * phi(k-1)
%
%		           1       P(k-2) * phi(k-1) * phi(k-1)' * P(k-2)
%		P(k-1) = ------ * ----------------------------------------
%		         lambda    lambda + phi(k-1)' * P(k-2) * phi(k-1)
%
%	where:
%
%		theta:		the parameter estimates
%		phi:		the state vector
%		P:		the covariance matrix
%		lambda:		the exponential data weighting factor
%
%	See also: SFUNC., "Adaptive Filtering, Prediction, and Control",
%	G. C. Goodwin & K. S. Sin.

%	Copyright (c) 1990-94 by The MathWorks, Inc.
%	Rick Spada 6-17-92.


if abs(flag) == 2	% flag = 2 --> real time hit
  
  % sample hit, return the next discrete states, which are the
  % next parameter estimates

  theta = x(1:nstates);					% parameter estimates
  P = zeros(nstates,nstates);				% get covariance matrix
  P(:) = x(nstates+1:nstates+nstates*nstates);
  yk = u(nstates + 1);					% system output
  phi = u(1:nstates);					% state vector
  est_err = yk - phi' * theta;				% estimation error
  den = lambda + phi' * P * phi;			% lambda + phi' * P * phi
  theta_new = theta + P * phi * (est_err / den);	% new parameter estimates
  Pnew = (P - P * phi * phi' * P / den) / lambda;	% new covariance
  sys = [theta_new', Pnew(:)']';			% return them
  
elseif flag == 4	% flag = 4 --> Return next sample hit

  sys = [];

elseif flag == 0	% flag = 0 --> Return sizes of parameters and initial conditions
  
  sys(1) = 0;				% 0 continuous states
  sys(2) = nstates+nstates*nstates;	% enough discrete states to hold the estimates
  					% and the covariance matrix
  sys(3) = nstates;			% nstate estimate outputs
  sys(4) = nstates+1;			% nstate+1 (regression vector + system output)
  					% inputs
  sys(5) = 0;				% 0 roots
  sys(6) = 0;				% no direct feedthrough
  sys(7) = 1;				% 1 sample time

  % initialize the covariance matrix and initial estimates
  P = eye(nstates, nstates) * 1e6;
  x0 = [ones(nstates,1)', P(:)']';
  ts = [dt, 0];
  
elseif flag == 3	% flag = 3 --> Return outputs, only at sample hits
  
    sys(:) = x(1:nstates);

else

  sys = [];
  
end