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