www.gusucode.com > rctobsolete 工具箱 matlab源码程序 > rctobsolete/robust/@lti/ohkapp.m

    function [ax,bx,cx,dx,ay,by,cy,dy,aug] = ohkapp(varargin)
%OHKAPP Optimal Hankel norm approximation (stable plant).
%
% [SS_X,SS_Y,AUG] = OHKAPP(SS_,MRTYPE,IN) or
% [AX,BX,CX,DX,AY,BY,CY,DY,AUG] = OHKAPP(A,B,C,D,MRTYPE,IN) produces
%   an optimal Hankel norm approximation via descriptor system.
%   (AX,BX,CX,DX) is the reduced model, (AY,BY,CY,DY) is the
%   anticausal part of the solution.
%   The infinity norm of (G - Ghed) <= sum of the Kth Hankel s.v. to
%   nth Hankel s.v. times 2. NO BALANCING IS REQUIRED.
%
%   mrtype = 1, in = order "k" of the reduced model.
%
%   mrtype = 2, find a reduced order model such that the total error is
%             less than "in".
%
%   mrtype = 3, display Hankel singular values, prompt for order "k" (in
%             this case, no need to specify "in").
%
%   AUG = [max. Hank SV, state(s) removed, Error Bound, Hankel SV's].

% R. Y. Chiang & M. G. Safonov 4/11/87
% Copyright 1988-2004 The MathWorks, Inc.
%--------------------------------------------------------------------
%

nag1=nargin;
[emsg,nag1,xsflag,Ts,a,b,c,d,mrtype,in]=mkargs5x('ss',varargin); error(emsg);
if Ts, error('LTI inputs must be continuous time (Ts=0)'), end

[ma,na] = size(a);
[md,nd] = size(d);
[hsv,p,q] = hksv(a,b,c);
%
if mrtype == 1
   kk = in+1;
   r = 1;
end
%
if mrtype == 2
   tails = 0;
   kk = 1;
   r = 1;
   for i = ma:-1:1
       tails = tails + hsv(i);
       if 2*tails > in
          kk = i+1;
          break
       end
   end
end
%
if mrtype == 3
   format short e
   format compact
   [mhsv,nhsv] = size(hsv);
   if mhsv < 60
      disp('    Hankel Singular Values:')
      hsv'
      for i = 1:mhsv
        if hsv(i) == 0
           hsvp(i) = eps;
        else
           hsvp(i) = hsv(i);
        end
      end
      disp(' ')
      disp('                              (strike a key to see the plot ...)')
      pause
      subplot
      plot(20*log10(hsvp),'--');hold
      plot(20*log10(hsvp),'*');grid on;hold
      ylabel('DB')
      title('Hankel Singular Values')
      pause
      in = input('Please assign the order of the K_th Hankel MDA: ');
      r = input('Please input the Multiplicity of the (K+1)th Hankel SV: ');
   else
      disp('    Hankel Singular Values:')
      hsv(1:60,:)'
      disp('                              (strike a key for more ...)')
      pause
      hsv(61:mhsv,:)'
      for i = 1:mhsv
        if hsv(i) == 0
           hsvp(i) = eps;
        else
           hsvp(i) = hsv(i);
        end
      end
      disp(' ')
      disp('                              (strike a key to see the plot ...)')
      pause
      subplot
      plot(20*log10(hsvp),'--');hold
      plot(20*log10(hsvp),'*');grid on;hold
      ylabel('DB')
      title('Hankel Singular Values')
      pause
      in = input('Please assign the order of the K_th Hankel MDA: ');
      r = input('Please input the Multiplicity of the (K+1)th Hankel SV: ');
    end
    kk = in+1;
    format loose
end
%
% ---- Keep all the state:
%
if kk > na
   ax = a; bx = b;
   cx = c; dx = d;
   ay = zeros(ma,na);   by = zeros(ma,nd);
   cy = zeros(md,na);   dy = zeros(md,nd);
   aug = [0 0 0 hsv'];
   if xsflag
      ax = mksys(ax,bx,cx,dx);
      bx = mksys(ay,by,cy,dy);
      cx = aug;
   end
   return
end
%
ro = hsv(kk,1);
bnd = ro*r;
kp1 = 0;
for i = 1 : na
    if hsv(i,1) < ro
       kp1 = kp1 + 1;
       bnd = bnd + hsv(i,1);
    end
end
strm = r + kp1;
aug = [hsv(1,1) strm 2*bnd hsv'];
%
aa = ro^2*a' + q*a*p;
bb = q*b;
cc = c*p;
dd = d;
ee = q*p - ro*ro*eye(ma);
[u,s,v] = svd(ee);
%
u1 = u(:,1:(na-r));
u2 = u(:,(na-r+1):na);
v1 = v(:,1:(na-r));
v2 = v(:,(na-r+1):na);
%
sigi = inv(u1'*ee*v1);
a11 = u1'*aa*v1;
a12 = u1'*aa*v2;
a21 = u2'*aa*v1;
a22i = pinv(u2'*aa*v2);
b1 = u1'*bb; b2 = u2'*bb;
c1 = cc* v1; c2 = cc*v2;
%
axy = sigi * (a11 - a12*a22i*a21);
bxy = sigi * (b1  - a12*a22i*b2);
cxy = c1 - c2 * a22i * a21;
dxy = dd - c2 * a22i * b2;
%
[ax,bx,cx,dx,ay,by,cy,dy,msat] = stabproj(axy,bxy,cxy,dxy);
%
if xsflag
   ax = mksys(ax,bx,cx,dx);
   bx = mksys(ay,by,cy,dy);
   cx = aug;
end
%
% ------ End of OHKAPP.M --- RYC/MGS %