bilinid.m

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

function [model,extra] = bilinid(data,n,i,Dzero,algorithm,lsmethod)
%BILINID   Identify a bilinear model from input-output data.
%
%   [M,EXTRA] = bilinid(DATA)
%   [M,EXTRA] = bilinid(DATA,n)
%   [M,EXTRA] = bilinid(DATA,n,k)
%   [M,EXTRA] = bilinid(DATA,n,k,DMAT)
%   [M,EXTRA] = bilinid(DATA,n,k,DMAT,ALG)
%   [M,EXTRA] = bilinid(DATA,n,k,DMAT,ALG,LS)
%
%   DATA is the input-output data given as an IDDATA object. M is the
%   estimated model returned as a discrete-time BILIN state-space system.
%
%   n is the system order (optional):
%    (1) If n = [] then the algorithm will automatically select the order
%        such that n <= 10 (default). 
%    (2) If n is a scalar, then n is the system order. 
%    (3) If entered as a vector (e.g. [1 2 3 4 5]), a plot of singular
%        values will be given and the user will be prompted to select an order.
%
%   k is the block size given as a positive integer (optional). If not
%   specified, k is chosen to be as large as possible while still trying to
%   be compatible with the size of DATA and choice of ALG. It is recommended
%   that k >= max(n). 
%
%   DMAT determines whether the D matrix is estimated or set to zero (optional).
%     'Estimate': Estimate D (default).
%     'Zero'    : Set D = 0.
%
%   ALG determines the choice of system identification algorithm and can be
%   one of {'general' | 'fast' | 'accurate'} (optional).  ALG='fast' is only
%   allowed if the number of outputs < min(n), and ALG='accurate' is only
%   allowed if the number of outputs >= max(n). If ALG is not specified,
%   then the most appropriate algorithm is automatically chosen based on the
%   number of outputs, i.e. if number of outputs < max(n), then
%   ALG='general'. If number of outputs >= max(n), then ALG='accurate'.
%
%   LS determines the choice of least squares method (optional) used in
%   estimating the system matrices and can be one of 'cls' (default) or
%   'ols'. LS='cls' will result in using a constrained least squares method
%   and EXTRA.ls=0. LS='ols' will result in the use of an ordinary least
%   squares method. If LS='ols', then EXTRA.ls is a matrix containing
%   information about the accuracy of the estimation and should be close to
%   0 for a good estimate.
%
%   EXTRA is a structure containing additional information about the model
%   and data. EXTRA.SV is a vector containing the singular values resulting
%   from the SVD decomposition. EXTRA.Q, EXTRA.R and EXTRA.S are the
%   matrices that form the noise covariance matrix EXTRA.COV = [Q S; S' R].
%
%   See also BILIN, IDDATA, BILIN/SIM, BILIN/COMPARE, BILIN/PE.
%
% 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.

% Check the arguments

if nargin < 1
  error('BILINID needs at least one input argument.')
end

% Turn the data into row vectors
%    Y=[y(1) y(2) ...  y(N)];  U=[u(1) u(2) ...  u(N)]

if ~isa(data,'iddata')
  error('DATA should be an IDDATA object.')
end

Y = get(data,'OutputData');
U = get(data,'InputData');
Ts = get(data,'Ts');
%clear data

Y = Y';
[p,cy] = size(Y);

if p < 1
  error('Need a non-empty output vector in DATA.')
end

U = U';
m = size(U,1);

if m < 1
  error('Need a non-empty input vector in DATA.')
end

if nargin > 1
  if ~isempty(n)
	n = sort(n(:));
	if any(n < 1 | round(n) ~= n)
	  error('Only positive integers are allowed for the choice of system order.')
	end
%	if max(n) > i*p
%	  warning('The block size might be too small to support this choice of order.')
%	end
  end
else
  n = [];
end

if nargin > 3
  if ~isempty(Dzero)
	switch lower(Dzero)
	 case {'e','estimate','estimated'}
	  isDzero=0;
	 case {'z','zero','zeros'}
	  isDzero=1;
	 otherwise
	  error('Invalid choice for DMAT. Select ''Estimate'' or ''Zero''.')
	end
	clear Dzero;
  else
	isDzero=0;
  end
else
  isDzero=0;
end

if nargin > 4
  if isempty(algorithm)
	if isempty(n) 
	  if p >= 10 % Maximum n will be 10
		algorithm = 42;
	  else
		algorithm = 3;
	  end
	else
	  if p < max(n)
		algorithm = 3;
	  else
		algorithm = 42;
	  end
	end
  else
	algorithm = lower(algorithm);
	switch algorithm
	 case {'general','g'}
	  algorithm = 3;
	 case {'fast','f'}
	  if ~isempty(n)
		if p < min(n)
		  algorithm = 41;
		else
		  error('ALG = ''fast'' is not allowed if the number of outputs is not less than min(n).')
		end
	  else
		error('ALG = ''fast'' is not allowed if n is empty.')
	  end
	 case {'accurate','a'}
	  if ~isempty(n)
		if p >= max(n)
		  algorithm = 42;
		else
		  error('ALG = ''accurate'' is not allowed if the number of outputs is less than max(n).')
		end
	  else
		error('ALG = ''accurate'' is not allowed if n is empty.')
	  end
	 otherwise
	  error('Invalid choice for ALG - choose ''general'', ''fast'' or ''accurate''.')
	end
	%	 case {'3','iii'}
	%	  algorithm = 3;
	%	 case {'4.1','iv.i','iv.1'}
	%	  algorithm = 41;
	%	 case {'4.2','iv.ii','iv.2'}
	%	  algorithm = 42;
	%	 otherwise
	%	  error('Invalid choice for ALG - choose ''III'', ''IV.I'' or ''IV.II''.')
	%  end
  end
else
  if isempty(n) 
	if p >= 10 % maximum n is 10
	  algorithm = 42;
	else
	  algorithm = 3;
	end
  else
	if p < max(n)
	  algorithm = 3;
	else
	  algorithm = 42;
	end
  end
end

switch algorithm
 case 3
  disp(' ')
  disp('Using general four-block, deterministic-stochastic algorithm.')
 case 41
  disp(' ')
  disp('Number of outputs < system order.')
  disp('Using fast four-block, deterministic-stochastic algorithm.')
 case 42
  disp(' ')
  disp('Number of outputs >= system order.')
  disp('Using accurate four-block, deterministic-stochastic algorithm.')
end


if nargin > 5
  if isempty(lsmethod)
	lsmethod = 'cls';
  else
	lsmethod = lower(lsmethod);
	switch lsmethod
	 case {'ols','o','ordinaryleastsquares'}
	  lsmethod = 'ols';
	 case {'cls','c','constrainedleastsquares'}
	  lsmethod = 'cls';
	 otherwise
	  error('Invalid choice for LSMETHOD - choose ''ols'' or ''cls''.')
	end
  end
else
  lsmethod = 'cls';
end

switch lsmethod
 case 'ols'
  disp('Ordinary least squares will be used when estimating system matrices.')
 case 'cls'
  disp('Constrained least squares will be used when estimating system matrices.')
end
disp(' ')


if nargin > 2
  if isempty(i)
	i=1;
	ei=0;
	for j=0:i-1
	  ei=ei+m*(m+1)^(j);
	end 
	if algorithm ~= 3
	  di=0;
	  for j=0:i-1
		di=di+p*(m+1)^(j);
	  end  
	end
	switch algorithm
	 case 3
	  while cy >= ei^(2)+ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1)
		i=i+1;
		ei=0;
		for j=0:i-1
		  ei=ei+m*(m+1)^(j);
		end 
	  end
	 case 41
	  while cy >= ei^(2)+ei+di 
		i=i+1;
		ei=0;
		for j=0:i-1
		  ei=ei+m*(m+1)^(j);
		end 
		di=0;
		for j=0:i-1
		  di=di+p*(m+1)^(j);
		end 
	  end
	 case 42
	  while cy >= di+3*(ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1))
		i=i+1;
		ei=0;
		for j=0:i-1
		  ei=ei+m*(m+1)^(j);
		end 
		di=0;
		for j=0:i-1
		  di=di+p*(m+1)^(j);
		end 
	  end
	end
	i=i-1;
	disp(sprintf('Block size k = %d.',i))
  else
	if i < 1 | round(i) ~= i
	  error('Block size k should be a positive integer.')
	end
	ei=0;
	for j=0:i-1
	  ei=ei+m*(m+1)^(j);
	end 
	if algorithm ~= 3
	  di=0;
	  for j=0:i-1
		di=di+p*(m+1)^(j);
	  end  
	end
	switch algorithm
	 case 3
	  if (cy < ei^(2)+ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1))
		error('Not enough data points for this choice of block size.')
	  end
	 case 41
	  if (cy < ei^(2)+ei+di )
		error('Not enough data points for this choice of block size.')
	  end
	 case 42
	  if (cy < di+3*(ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1)))
		error('Not enough data points for this choice of block size.')
	  end
	end
   end
else
  i=1;
  ei=0;
  for j=0:i-1
	ei=ei+m*(m+1)^(j);
  end 
  if algorithm ~= 3
	di=0;
	for j=0:i-1
	  di=di+p*(m+1)^(j);
	end  
  end
  switch algorithm
   case 3
	while cy >= ei^(2)+ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1)
	  i=i+1;
	  ei=0;
	  for j=0:i-1
		ei=ei+m*(m+1)^(j);
	  end 
	end
   case 41
	while cy >= ei^(2)+ei+di 
	  i=i+1;
	  ei=0;
	  for j=0:i-1
		ei=ei+m*(m+1)^(j);
	  end 
	  di=0;
	  for j=0:i-1
		di=di+p*(m+1)^(j);
	  end 
	end
   case 42
	while cy >= di+3*(ei+(m/2)*(m+1)^(i)+p*((m+1)^(i)-1))
	  i=i+1;
	  ei=0;
	  for j=0:i-1
		ei=ei+m*(m+1)^(j);
	  end 
	  di=0;
	  for j=0:i-1
		di=di+p*(m+1)^(j);
	  end 
	end
  end
  i=i-1;
  disp(sprintf('Block size k = %d.',i))
end

clear ei
if algorithm ~= 3
  clear di
end


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

% ************************************************** %
%                   STEP 1                           %
%                                                    %
%        Construct Y^{f} and U^{f,u,y}               %
%                                                    %
% ************************************************** %

% Call Blochank
% Call Bloctoep

% Add some some noise to the input in order to make UUU full row rank if
% the rank condition in the following is not met.

% MOD ECK 19/06/02 - following code is incorrect
if algorithm == 42
  
  jj = cy-3*i+1;  % jj=N-3i+1, N=cy; 
  UUU=blochank(U(:),3*i,cy); % UUU=[u(1),...,u(3i);...;u(N-3i+1),...,u(N)]
  if rank(UUU)~=3*i*m
	error('The input block Hankel matrix is not full rank. Add some noise to the input.')
  end	
  YYY=blochank(Y(:),3*i,cy);  % YYY=[y(1),...,y(3i);...;y(N-3i+1),...,y(N)]
  
else
  
  jj = cy-4*i;  % jj=N-4i+1, N=cy; 
  UUU=blochank(U(:),4*i+1,cy); % UUU=[u(1),...,u(3i);...;u(N-3i+1),...,u(N)]
  if rank(UUU)~=(4*i+1)*m
	error('The input block Hankel matrix is not full rank. Add some noise to the input.')
  end	
  YYY=blochank(Y(:),4*i+1,cy);  %  YYY=[y(1),...,y(3i);...;y(N-3i+1),...,y(N)]
  
end					

clear U Y cy
%size(UUU)