extu.m

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

function U = extu(sys,s,method)
% EXTU Extracts direction information of the columns of C from the system.
%
%   U = EXTU(SYS,S,METHOD) extracts the directions of the column vectors of
%   the C matrix of the state-space system SYS, where METHOD indicates which
%   method will be used later on to identify the system. U(:,j) should have
%   unity 2-norm and is the direction of the column vector C(:,j).
%
% 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.

n = length(s);

if strcmp(method,'mp')
  K = sys.b; % Recall that B has been replaced with K
  try
	Ct = inv(sys.d)*sys.c;
	U = Ct*0; % Assign memory space for U to speed up computations
  catch
	% There should be no problem with inverting D if the original system
	% is minimum-phase and the error feedthrough term is the identity matrix
	error('The system is not minimum-phase.');	
  end
else
  B = sys.b;
  C = sys.c;
  D = sys.d;
  U = C*0; % Assign memory space for U to speed up computations
end

if strcmp(method,'pr')
  Ds = sys.d + sys.d';
  if any(real(eig(Ds)) <= 0)
	% Ds has to be positive definite for the results to hold - see
    % CM97. Unfortunately this constraint cannot be implemented in a simple
    % fashion. If an error occurs, please try a different estimate.
	error('Dc + Dc'' is not positive definite - please try a different estimate.')
  end
end

clear sys
warning on
for j = 1:n
  switch method
   case 'sb'
	% CM97, Equation (9)
	U(:,j) = C(:,j)/sqrt(B(j,:)*B(j,:)');
   case 'mp'	
	% CM97, Equation (29)
	U(:,j) = (Ct(:,j)-s(j)*K(j,:)')/sqrt(1+s(j)^2)/sqrt(K(j,:)*K(j,:)'); 
   case 'pr'	
	% CM97, Equation (44)
	U(:,j) = inv(sqrtm(Ds))*C(:,j)/sqrt(B(j,:)*inv(Ds)*B(j,:)'); 
   case 'ncf'
	% CM97, Equation (55)
	U(:,j) = C(:,j)*sqrt(1-s(j)^2)/sqrt(B(j,:)*B(j,:)');
   otherwise
	error('Unknown method.')
  end
  warning on
  Unorm = U(:,j)'*U(:,j);
  if Unorm == 0
	error('One of the columns of U == 0. This is not allowed.')
  end
  if abs(1-Unorm) > 1e-6
	warning(sprintf('Warning: Solution might be inaccurate. 1-||U(:,%d)||_2 = %4.3e.',j,1-sqrt(Unorm)))
  end
end

if any(any(isnan(U) | isinf(U)))
  error('The system estimates are ill-conditioned - cannot extract U.')
end