balpem.m

剑桥大学用于线性和双线性动力学系统辨识的工具箱

function mfdb = balpem(data,mid,bd,method,varargin)
%BALPEM   Identifies a balanced state-space model from input-output data.
%
%   BALPEM uses a balanced parameterization to improve on an initial
%   estimate of a system using the prediction error algorithm PEM and the
%   IDGREY class. The estimated model is stable, minimal and balanced in
%   some sense. 
%
%   Usage:
%      M = balpem(DATA,Mi)
%      M = balpem(DATA,Mi,BD)
%      M = balpem(DATA,Mi,BD,ALG)
%      M = balpem(DATA,Mi,BD,ALG,Property_1,Value_1,...,Property_n,Value_n)
%
%   Inputs:  
%      DATA : Input-output data given as an IDDATA object.
%      Mi   : Initial estimate of the system given as an IDSS object. 
%      BD   : Options for estimating the B and D  matrices (optional):
%               'Estimate' : Estimate B and D matrix (default).
%               'ZeroD'    : Estimate B and set D=0 (valid if ALG='sb' or 'mp').
%               'ZeroB'    : Set B=0 and estimate D (valid if ALG='mp').
%               'ZeroBD'   : Set B=0 and D=0 (valid if ALG='mp').
%      ALG  : Choice of balanced parameterization (optional):
%               'sb'  : Stable, balanced parameterization (default).
%               'mp'  : Minimum-phase balanced parameterization.
%      Property_n,Value_n: See IDPROPS ALGORITHMS or IDPROPS IDGREY for a
%               list of possible Property/Value pairs.
%
%   Outputs:
%      M : The estimated state-space model given as an IDSS model object:
%               x[k+1] = A x[k] + B u[k] + K e[k];     x[0] = X0
%               y[k]   = C x[k] + D u[k] + e[k]
%
%   See also IDDATA, IDSS, SUBID3B, N4SID, PEM, IDGREY.
%
% CUED System Identification Toolbox.
% Cambridge University Engineering Department.
% Copyright (C) 1998-2002. All Rights Reserved.

% Version 1.00, Date: 01/06/2002
% Created by H. Chen and E.C. Kerrigan.

% The algorithm is described in C.T. Chou and J.M. Maciejowski, "System
% Identification Using Balanced Parameterizations", IEEE Transactions on
% Automatic Control 42:7 (1997) pp. 956-974.

%%%%%%%%%%%%%%%%%%%%%%%
% Check the arguments %
%%%%%%%%%%%%%%%%%%%%%%%

warning on

if nargin < 2
  error('BALPEM needs at least two input arguments.')
end

if nargin > 1
  if isempty(mid)
	error('Mi has not been given.')
  else
	if ~isa(mid,'idss') | mid.Ts == 0
	  error('Mi must be a discrete-time IDSS object.')
	end
  end
end

if nargin > 3
  if isempty(method)
	method = 'sb';
  else
	method = lower(method);
	if ~(strcmp(method,'sb') | strcmp(method,'mp'))
	  error('Invalid choice for ALG. Must be ''sb'' or ''mp'.')
	end	
  end
else
  method = 'sb';
end

% aux.Best == 1 if B is to be estimated, 0 if B == 0
% aux.Dest == 1 if D is to be estimated, 0 if D == 0
if nargin > 2
  switch method
   case {'sb'}
	if isempty(bd)
	  aux.Dest = 1; 
	else 
	  switch lower(bd)
	   case {'zerod','zd','z','zero','ebzd'}
		aux.Dest = 0;
	   case {'estimate','e','est','ebd','ebed'}
		aux.Dest = 1;
	   otherwise
		error('Invalid choice for BD. Choose ''Estimate'' or ''ZeroD''.')
	  end
	end
   case 'mp'
	if isempty(bd)
	  aux.Best = 1;
	  aux.Dest = 1; 
	else
	  switch lower(bd)
	   case {'zerob','zb','zbed'}
		aux.Best = 0;
		aux.Dest = 1;
	   case {'zerod','zd','ebzd'}
		aux.Best = 1;
		aux.Dest = 0;
	   case {'estimate','e','est','ebd','ebed'}
		aux.Best = 1;
		aux.Dest = 1;
	   case {'zerobd','zbd','zero','z','zbzd'}
		aux.Best = 0;
		aux.Dest = 0;
	   otherwise
		error('BD must be one of ''Estimate'', ''ZeroB'', ''ZeroD'' or ''ZeroBD''.')
	  end	
	end
   case 'ncf'
	aux.Dest = 1;
	if ~isempty(bd) 
	  if ~strcmpi(bd,'Estimate')
		error('BD can only be ''Estimate''.')
	  end
	end
  end
else
  switch method
   case 'sb'
	aux.Dest = 1;
   case 'mp'
	aux.Best = 1;
	aux.Dest = 1;
   case 'ncf'
	aux.Best = 1;
	aux.Dest = 1;
  end
end

[aux.n,aux.m] = size(mid.b); % Number of states n and inputs m
aux.p = size(mid.c,1);       % Number of outputs p
%if strcmp(method,'ncf')
%  aux.m = aux.m-aux.p; % Make correction if ALG is 'ncf'
%end

if ~all(all(mid.k == 0)) & ~strcmp(method,'mp')
  warning(sprintf('Mi.K is not allowed to be non-zero if ALG is ''sb''. Setting Mi.K = 0.'))
  mid.k = mid.k*0;
end

%nk = mid.nk;
if aux.Dest == 1
  if ~all(mid.nk == 0)
	warning('Mi.nk is not equal to [0 ... 0]. Estimated model M might be poor.')
  end
else
  if ~all(mid.nk == 1)
	warning('Mi.nk is not equal to [1 ... 1]. Estimated model M might be poor.')
  end
  mid.d = 0*mid.d;
end
%mid.nk = nk;

Ts = mid.Ts; % Sample time

if strcmp(method,'mp') | strcmp(method,'sb')
  if max(abs(eig(mid.a))) >= 1
	error('Mi is not stable. The eigenvalues of A do not lie inside the unit circle.')
  end 
end

if strcmp(method,'mp') 
  if max(abs(eig(mid.a-mid.k*mid.c)) >= 1)
	error('Mi is not minimum-phase. The eigenvalues of A-K*C do not lie inside the unit circle.')
  end 
  if aux.Best == 0
	if ~all(all(mid.b == 0))
	  mid.b = 0*mid.b;
	  warning('Mi.B is not equal to zero. Estimated model M might be poor.')
	end
  end
  % Replace B with K and ensure that error feedthrough term is the identity matrix
  Dmp = mid.d; % Discrete-time estimate of D
  Bmp = mid.b; % Discrete-time estimate of B
  mid = ss(mid.a,mid.k,mid.c,eye(aux.p),Ts);
else
  mid = ss(mid.a,mid.b,mid.c,mid.d,Ts); % Convert from IDSS to SS format
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                              BEGIN ALGORITHM                              %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert to a continuous-time, balanced SS system using the bilinear
% transform, followed by calling BALREAL. If ALG is 'mp', then the next
% command actually balances (A,K,C,I).

try
  mic = d2c(mid,'tustin'); % Use the bilinear transform
  switch method
   case 'mp'
	%ismpbalanced(mic)
	[micb,s,Tbal] = mpbal(mic);
	%ismpbalanced(micb)
   case 'sb'
	[micb,s,Tbal] = balreal(mic);
	%isbalanced(micb)
   case 'ncf'
	micb=mic;
  end
catch
  switch method
   case 'sb'
	error('(A,B,C,D) of Mi has to be a minimal, stable system.')
   case 'mp'
	error('(A,K,C,I) of Mi has to be a reachable, stable, minimum-phase system with observable inverse.')
   case 'ncf'
	error('(A,B,C,D) of Mi has to be a minimal system.')
  end
end

% Correct the sign of the elements of B (K) and extract the data from B (K)
[micb,aux.Bind,Bdat,Tsig] = signextb(micb);

% Extract directional information about the columns of C
aux.U = extu(micb,s,method);

% Convert Hankel singular values to a vector of estimated parameters
par = conhsv(s,method);

if aux.Dest == 1 % If D needs to be estimated
  switch method
   case 'sb'
	par = [par; reshape(mid.d,aux.p*aux.m,1)]; % d.t. estimate
   case 'mp'
	par = [par; reshape(Dmp,aux.p*aux.m,1)];   % d.t. estimate
  end
end

if strcmpi(method,'ncf')
  par = [par; reshape(micb.d,aux.p*aux.m,1)]; % c.t. estimate
end

if strcmp(method,'mp') 
  aux.Kind = aux.Bind; % Recall that B has been replaced by K, so Bind
  aux.Bind = [];       % actually relates to K , hence Bind is redundant
  if aux.Best == 1     % If B needs to be estimated
	% Next two lines implement the bilinear transform on Bmp followed with
    % the state transform Tsig*Tbal, as done above to (A,K,C,I) of Mi.
	sysc = d2c(ss(mid.a,Bmp,mid.c,Dmp,Ts),'tustin');
	sysc = ss2ss(sysc,Tsig*Tbal);
	Bmp = sysc.b; % Bmp now represents the continuous-time estimate of B
	par = [par; reshape(Bmp,aux.n*aux.m,1)] ;
  end
end

par = [par; Bdat];  % Add estimated B (K) data

% If method = 'sd', then PAR = [Sc; Dd; Bc]
% If method = 'mp'  then PAR = [Sc; Dd; Bc; Kc]
% If merhod = 'ncf' then PAR = [Sc; Dc; Bc]

mgrey = idgrey(method,par,'d',aux,Ts,varargin{:});
mfdb = pem(data,mgrey,varargin{:}); % MfDB in IDGREY form
mfdb = idss(mfdb); % Convert to IDSS form
switch method
 case{'sb'}
  mfdb.EstimationInfo.Method='BALPEM - Standard balanced parameterization';
 case{'mp'}
  mfdb.EstimationInfo.Method='BALPEM - Minimum-phase balanced parameterization';
end

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                          END ALGORITHM                                    %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% *** last line of balpem.m ***