mpbal.m

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

function [sysb,G,T,Ti] = mpbal(sys)
%MPBAL  Minimum-phase balancing of a minimum-phase state-space realization.
%
%   SYSb = MPBAL(SYS) returns a stable, minimum-phase balanced realization
%   of the system SYS in the sense that the controllability gramian of SYSb
%   and the observability gramian of the inverse of SYSb are equal and
%   diagonal. SYS has to be a controllable, stable, minimum-phase SS object
%   and the inverse of SYS must be observable.
%
%   [SYSb,G,T,Ti] = MPBAL(SYS) also returns a vector G containing the
%   diagonal of the Gramian of the balanced realization.  The matrix T is
%   the state transformation xb = Tx used to convert SYS to SYSb, and Ti is
%   its inverse.
%
%   See also ISMPBALANCED, BALREAL, GRAM, SS.
%
% 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.

if ~isa(sys,'ss') 
  error('SYS is not a continuous-time SS object.')
else
  if sys.Ts ~= 0
	error('SYS is not a continuous-time SS object.')
  end
end

try
  sysinv = inv(sys);
catch
  error('SYS is not invertible: The D matrix is singular.')
end

if max(real(eig(sys))) >= 0 
  error('SYS is not asymptotically stable.')
end

if max(real(eig(sysinv))) >= 0 
  error('SYS is not minimum-phase: A-B*inv(D)*C is not Hurwitz.')
end

try
  P = gram(sys,'c');
  R = chol(P);
catch
  error('SYS is not controllable.')
end

try
  Q = gram(sysinv,'o');
  [U,S,V]=svd(R*Q*R'); % Should have U = V 
  %if norm(U-V) > 1e-6
  %	error('U is not equal to V in the SVD.')
  %end
  GRAM = sqrt(S);
  G = diag(GRAM);
  T = sqrt(GRAM)*inv(U)*inv(R');
  Ti = inv(T);
  %Ti = R'*U*inv(sqrt(GRAM));
  %T = inv(Ti);
catch
  error('The inverse of SYS is not observable.')
end

if norm(U-V) > 1e-6
  error('U is not equal to V in the SVD.')
end

sysb = ss2ss(sys,T);

%Check: Should have GRAM = Phat = Qhat = Pb = Qb
%Phat = T*P*T';
%Qhat = Ti'*Q*Ti;
%sysbinv = inv(sysb);
%Pb = gram(sysb,'c');
%Qb = gram(sysbinv,'o');
%GRAM
%Pb
%Qb
%Phat
%Qhat