home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / OHKAPP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.1 KB  |  168 lines
  1. function [ax,bx,cx,dx,ay,by,cy,dy,aug] = ohkapp(Z1,Z2,Z3,Z4,Z5,Z6)
  2. % [SS_X,SS_Y,AUG] = OHKAPP(SS_,TYPE,IN) or
  3. % [AX,BX,CX,DX,AY,BY,CY,DY,AUG] = OHKAPP(A,B,C,D,TYPE,IN) produces 
  4. %   an optimal Hankel norm approximation via descriptor system.
  5. %   (AX,BX,CX,DX) is the reduced model, (AY,BY,CY,DY) is the 
  6. %   anticausal part of the solution.
  7. %   The infinity norm of (G - Ghed) <= sum of the Kth Hankel s.v. to 
  8. %   nth Hankel s.v. times 2. NO BALANCING IS REQUIRED.
  9. %
  10. %   Type = 1, in = order "k" of the reduced model.
  11. %
  12. %   Type = 2, find a reduced order model such that the total error is 
  13. %             less than "in".
  14. %
  15. %   Type = 3, display Hankel singular values, prompt for order "k" (in 
  16. %             this case, no need to specify "in").
  17. %
  18. %   AUG = [max. Hank SV, state(s) removed, Error Bound, Hankel SV's].
  19.  
  20. % R. Y. Chiang & M. G. Safonov 4/11/87
  21. % Copyright (c) 1988 by the MathWorks, Inc.
  22. % All Rights Reserved.
  23. %--------------------------------------------------------------------
  24. %
  25.  
  26. inargs = '(a,b,c,d,Type,in)';
  27. eval(mkargs(inargs,nargin,'ss'))
  28.  
  29. [ma,na] = size(a);
  30. [md,nd] = size(d);
  31. [hsv,p,q] = hksv(a,b,c);
  32. %
  33. if Type == 1
  34.    kk = in+1;
  35.    r = 1;
  36. end
  37. %
  38. if Type == 2
  39.    tails = 0;
  40.    kk = 1;
  41.    r = 1;
  42.    for i = ma:-1:1
  43.        tails = tails + hsv(i);
  44.        if 2*tails > in
  45.           kk = i+1;
  46.           break
  47.        end
  48.    end
  49. end
  50. %
  51. if Type == 3
  52.    format short e
  53.    format compact
  54.    [mhsv,nhsv] = size(hsv);
  55.    if mhsv < 60
  56.       disp('    Hankel Singular Values:')
  57.       hsv'
  58.       for i = 1:mhsv
  59.         if hsv(i) == 0
  60.            hsvp(i) = eps;
  61.         else
  62.            hsvp(i) = hsv(i);
  63.         end
  64.       end
  65.       disp(' ')
  66.       disp('                              (strike a key to see the plot ...)')
  67.       pause
  68.       subplot
  69.       plot(20*log10(hsvp),'--');hold
  70.       plot(20*log10(hsvp),'*');grid;hold
  71.       ylabel('DB')
  72.       title('Hankel Singular Values')
  73.       pause
  74.       in = input('Please assign the order of the K_th Hankel MDA: ');
  75.       r = input('Please input the Multiplicity of the (K+1)th Hankel SV: ');
  76.    else
  77.       disp('    Hankel Singular Values:')
  78.       hsv(1:60,:)'
  79.       disp('                              (strike a key for more ...)')
  80.       pause
  81.       hsv(61:mhsv,:)'
  82.       for i = 1:mhsv
  83.         if hsv(i) == 0
  84.            hsvp(i) = eps;
  85.         else
  86.            hsvp(i) = hsv(i);
  87.         end
  88.       end
  89.       disp(' ')
  90.       disp('                              (strike a key to see the plot ...)')
  91.       pause
  92.       subplot
  93.       plot(20*log10(hsvp),'--');hold
  94.       plot(20*log10(hsvp),'*');grid;hold
  95.       ylabel('DB')
  96.       title('Hankel Singular Values')
  97.       pause
  98.       in = input('Please assign the order of the K_th Hankel MDA: ');
  99.       r = input('Please input the Multiplicity of the (K+1)th Hankel SV: ');
  100.     end
  101.     kk = in+1;
  102.     format loose
  103. end
  104. %      
  105. if kk > na 
  106.    ax = a; bx = b;
  107.    cx = c; dx = d;
  108.    ay = zeros(ma,na);   by = zeros(ma,nd);
  109.    cy = zeros(md,na);   dy = zeros(md,nd);
  110.    aug = [0 0 0 hsv'];
  111.    return 
  112. end
  113. if kk < 1
  114.    ax = zeros(ma,na); bx = zeros(ma,nd);
  115.    cx = zeros(md,na); dx = zeros(md,nd);
  116.    ay = a; by = b;
  117.    cy = c; dy = d;
  118.    bnd = 2*sum(hsv); 
  119.    aug = [hsv(1,1) na bnd hsv'];
  120.    return
  121. end 
  122. %
  123. ro = hsv(kk,1);
  124. bnd = ro*r;
  125. kp1 = 0;
  126. for i = 1 : na
  127.     if hsv(i,1) < ro
  128.        kp1 = kp1 + 1;
  129.        bnd = bnd + hsv(i,1);
  130.     end
  131. end
  132. strm = r + kp1;
  133. aug = [hsv(1,1) strm 2*bnd hsv'];
  134. %
  135. aa = ro^2*a' + q*a*p;
  136. bb = q*b;
  137. cc = c*p;
  138. dd = d;
  139. ee = q*p - ro*ro*eye(ma);
  140. [u,s,v] = svd(ee);
  141. %
  142. u1 = u(:,1:(na-r));
  143. u2 = u(:,(na-r+1):na);
  144. v1 = v(:,1:(na-r));
  145. v2 = v(:,(na-r+1):na);
  146. %
  147. sigi = inv(u1'*ee*v1);
  148. a11 = u1'*aa*v1;
  149. a12 = u1'*aa*v2;
  150. a21 = u2'*aa*v1;
  151. a22i = pinv(u2'*aa*v2);
  152. b1 = u1'*bb; b2 = u2'*bb; 
  153. c1 = cc* v1; c2 = cc*v2;
  154. %
  155. axy = sigi * (a11 - a12*a22i*a21);
  156. bxy = sigi * (b1  - a12*a22i*b2);
  157. cxy = c1 - c2 * a22i * a21;
  158. dxy = dd - c2 * a22i * b2; 
  159. %
  160. [ax,bx,cx,dx,ay,by,cy,dy,msat] = stabproj(axy,bxy,cxy,dxy);
  161. %
  162. if xsflag
  163.    ax = mksys(ax,bx,cx,dx);
  164.    bx = mksys(ay,by,cy,dy);
  165.    cx = aug;
  166. end
  167. %
  168. % ------ End of OHKAPP.M --- RYC/MGS %