home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / SCHMR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.4 KB  |  134 lines
  1. function [am,bm,cm,dm,totbnd,hsv] = schmr(Z1,Z2,Z3,Z4,Z5,Z6)
  2. % [SS_M,TOTBND,HSV] = SCHMR(SS_,TYPE,NO) or
  3. % [AM,BM,CM,DM,TOTBND,HSV]=SCHMR(A,B,C,D,TYPE,NO) performs Schur method
  4. %    model reduction on G(s):=(a,b,c,d) such that the infinity-norm
  5. %    of the error (Ghed(s) - G(s)) <= sum of the 2(n-k) smaller Hankel
  6. %    singular values. For unstable G(s) the algorithm works by first 
  7. %    splitting G(s) into a sum of stable and antistable part. 
  8. % 
  9. %   Based on the "TYPE" selected, you have the following options:
  10. %
  11. %    1). TYPE = 1  --- no: size "k" of the reduced order model.
  12. %
  13. %    2). TYPE = 2  --- find a k-th order model such that the total error 
  14. %                      is less than "no".
  15. %
  16. %    3). TYPE = 3  --- disply all the Hankel SV prompt for "k" (in this
  17. %                      case, no need to specify "no").
  18. %
  19. %    TOTBND = Error bound, HSV = Hankel Singular Values.
  20.  
  21. % R. Y. Chiang & M. G. Safonov 9/13/87
  22. % Copyright (c) 1988 by the MathWorks, Inc.
  23. % All Rights Reserved.
  24. % ------------------------------------------------------------------------
  25. %
  26. disp('  ');
  27. disp('        - - Working on Schur Balanced Model Reduction - -');
  28.  
  29. inargs = '(a,b,c,d,Type,no)';
  30. eval(mkargs(inargs,nargin,'ss'))
  31.  
  32. if Type == 3
  33.    no = 0;
  34. end
  35. [ra,ca] = size(a);
  36. [ee,dd] = eig(a); 
  37. %
  38. % ------ Find the no. of stable roots :
  39. %
  40. m = 0;
  41. for k = 1 : ra
  42.     if real(dd(k,k)) < 0
  43.         m = m + 1;
  44.     end
  45. end
  46. %
  47. % ------ If completely stable :
  48. %
  49. kk = -1;
  50. if m == ra
  51.    [am,bm,cm,dm,augl,hsv,slbig,srbig] = schbal(a,b,c,d,Type,no);
  52.    kk = 1;
  53.    augr = [0 0];
  54.    totbnd = augl(1,2);
  55. end
  56. %
  57. % ------ If completely unstable :
  58. %
  59. if m == 0.
  60.    [am,bm,cm,dm,augr,hsv,slbig,srbig] = schbal(-a,-b,c,d,Type,no);
  61.    am = -am; bm = -bm; 
  62.    kk = 1;
  63.    augl = [0 0];
  64.    totbnd = augr(1,2);
  65. end
  66. %
  67. % ------ If having both stable & unstable parts :
  68. %
  69. if kk < 0
  70.   [al,bl,cl,dl,ar,br,cr,dr,msat] = stabproj(a,b,c,d);
  71.   [hsvl,pl,ql] = hksv(al,bl,cl);
  72.   [hsvr,pr,qr] = hksv(-ar,-br,cr);
  73.   hsv = [hsvl;hsvr];
  74.   [hsv,index] = sort(hsv);
  75.   hsv = hsv(ra:-1:1,:);
  76.   index = index(ra:-1:1,:);
  77.   if Type == 1
  78.       nor = 0;
  79.       for i = 1:no
  80.           if index(i) > msat
  81.              nor = nor + 1;           
  82.           end
  83.       end
  84.      nol = no-nor;
  85.    end
  86.    if Type == 2
  87.       tol = no;
  88.       tails = 0;
  89.       for i = ra:-1:1
  90.           tails = tails + hsv(i);
  91.           if 2*tails > tol
  92.              no = i;
  93.              break    
  94.           end
  95.       end
  96.       nor = 0;
  97.       for i = 1:no
  98.           if index(i) > msat
  99.              nor = nor + 1;           
  100.           end
  101.       end
  102.       nol = no-nor;
  103.       Type = 1;
  104.    end
  105.    if Type == 3
  106.       nol = no; nor = no;
  107.    end
  108.    [alf,blf,clf,dlf,augl,hsvl,tll,tlr] = schbal(al,bl,cl,dl,Type,nol);
  109.    [arf,brf,crf,drf,augr,hsvr,trl,trr] = schbal(-ar,-br,cr,dr,Type,nor);
  110.    arf = -arf;
  111.    brf = -brf;
  112.    totbnd = augl(1,2)+augr(1,2);
  113.    hsv = [hsvl;hsvr];
  114.    [am,bm,cm,dm] = addss(alf,blf,clf,dlf,arf,brf,crf,drf);
  115.    [ral,cal] = size(al); [rar,car] = size(ar);
  116.    if augl(1,1) == ral
  117.       am = arf; bm = brf; cm = crf; dm = d;
  118.    end
  119.    if augr(1,1) == rar
  120.       am = alf; bm = blf; cm = clf; dm = d;
  121.    end
  122. end
  123. %
  124. if xsflag
  125.    am = mksys(am,bm,cm,dm);
  126.    bm = totbnd;
  127.    cm = hsv;
  128. end
  129. %
  130. nn = augl(1,1) + augr(1,1);
  131. disp(' ')
  132. disp(['               ' int2str(nn), '    states removed !!'])
  133. %
  134. % ------- End of SCHMR.M --- RYC/MGS 9/13/87 %