home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / MUOPT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.8 KB  |  139 lines
  1. function [mu,ascaled,log_d,x] = muopt(T,K);
  2. % [MU,ASCALED,LOGD,X] = MUOPT(A) or MUOPT(A,K) produces the scalar upper bound 
  3. % MU on the real or complex structured singular value (ssv) computed via
  4. % the generalized Popov multiplier method of Safonov and Lee (1993 IFAC 
  5. % World Congress).
  6. % Input: 
  7. %     A  -- a pxq complex matrix whose ssv is to be computed
  8. % Optional input:
  9. %     K  -- an nx1 or nx2 matrix whose rows are the uncertainty block sizes  
  10. %           for which the ssv is to be evaluated; K must satisfy 
  11. %           sum(K) == [q,p]. Real uncertainty is indicated by multiplying the 
  12. %           corresponding row of K by minus one, e.g., if the second 
  13. %           uncertainty is real then set K(2,:)=[-1,-1]. If only the first 
  14. %           column of K is given then the uncertainty blocks are taken 
  15. %           to be square, as if K(:,2)=K(:,1).
  16. % Outputs:
  17. %     MU      -- an upper bound on the structured singular value of A
  18. %     ASCALED -- mu*(I-X)/(I+X) where X=D(mu*I-A)/(mu*I_A)*inv(D')
  19. %     LOGD    -- an n-vector containing the log of the square-root of
  20. %                the diagonal elements of the optimal generalized 
  21. %                Popov multiplier multiplier M=D'*D
  22. %           x -- the normalized eigenvector associated with the smallest
  23. %                eigenvalue of X+X'>0 .
  24.  
  25. % Authors: Penghin Lee and M. G. Safonov 7/92
  26. %
  27. % R. Y. Chiang & M. G. Safonov 
  28. % Copyright (c) 1988-92 by the MathWorks, Inc.
  29. % All Rights Reserved.
  30.  
  31. [n,c] = size(T); 
  32.  
  33. [junk,T,logdpsv]=psv(T,abs(K));
  34.  
  35. if nargin<2,
  36.    K=2*ones(n,1);
  37. elseif max(abs(K))==1,
  38.    K=.5*K(:,1)+1.5*ones(n,1);
  39. else
  40.    error('The present version of muopt.m requires 1x1 uncertainties')
  41. end
  42.  
  43. k_eps = 1e-02; k_low = 1e-14; k_upp = 1e+14; 
  44. kjw = 1; zp = 6; g = 1; knogood = 0; y = 1;
  45.  
  46.     while zp > 5;
  47.          kpres(g,1) = kjw;       
  48.      Tt =  ( (eye(n) - kjw*T ) )/(eye(n) + kjw*T); 
  49. % Tt is sectf(kjw*T,[-1,1],[0,inf])
  50.  
  51. % Initial guess. Any x with norm one
  52. xr = rand(n,1); xqi = xr / norm(xr);
  53.  
  54. [mult,xqo] = mixedalg(Tt,xqi,K);
  55.  
  56. if kjw <= 1;
  57.   if knogood == 0 & norm(mult) == 0;
  58.     kjw = 0.5 * kjw;
  59.   end
  60. end;  
  61.  
  62.   if kjw < 1;
  63.         if norm(mult) > 0 ; 
  64.          k_yes(y,1) = kjw; 
  65.          knogood = kpres (g,1);
  66.            
  67.            if kpres(g,1) < kpres (g-1,1);
  68.             knogood1 = kpres(g-1,1);
  69.            end;
  70.  
  71.          kjw = 0.5 * ( knogood1 - kpres(g,1) ) + knogood; 
  72.          
  73.          mult_yes = mult; Tt_yes = Tt;
  74.          x_yes = xqo; y = y+1;
  75.         end;
  76.  
  77.     if norm(mult) == 0 & knogood >0;
  78.        kjw = 0.5 * ( kpres(g,1) - knogood ) + knogood;
  79.     end;
  80.  
  81.   end; % if kjw< ..
  82.  
  83.   if kjw > 1;
  84.         if norm(mult) ==  0; 
  85.          knogood = kpres(g,1);
  86.             if kpres(g,1) > kpres(g-1,1);
  87.              knogood1 = kpres(g-1,1);
  88.             end;
  89.  
  90.          kjw = 0.5 * (kpres(g,1)-knogood1) + knogood1;
  91.         end;
  92.          
  93.       if norm(mult) > 0 & knogood > 0;
  94.            k_yes(y,1) = kjw;
  95.            kjw = kpres(g,1) + 0.5*(knogood - kpres(g,1));
  96.  
  97.            mult_yes = mult; Tt_yes = Tt;
  98.            x_yes = xqo; y = y+1;
  99.         end;
  100.  
  101.   end; % If kjw>..
  102.  
  103. if y > 2;
  104.         if ( abs( k_yes(y-1,1)-k_yes(y-2,1) ) < k_eps * k_yes(y-2,1) ) 
  105.          zp = zp-2;
  106.         end;
  107.  
  108.     if (k_yes(y-1,1) > k_upp) | (k_yes(y-1,1) < k_low);
  109.           break,
  110.            end;
  111.  
  112. end;
  113.    
  114.   if norm(mult) > 0 & knogood == 0;
  115.    k_yes(y,1) = kjw; kjw = 2 * kjw;  mult_yes = mult; 
  116.    x_yes = xqo; Tt_yes = Tt; y = y+1;
  117.   end  
  118.  
  119. g = g+1;
  120.         end; %while z
  121.  
  122. % Outputs
  123. mu = 1 / k_yes(y-1,1) ;
  124.  
  125.  for my = 1: 2*n; 
  126.    if mult_yes(my,1) == 0; 
  127.      mult_yes(my,1) = 1e-100;
  128.    end;
  129.  end;
  130.  
  131. multip=mult_yes(n+1:2*n)+sqrt(-1)*mult_yes(1:n);
  132. log_d = 0.5 * log(multip) + logdpsv;
  133. d=diag(sqrt(multip));
  134. T_sectf=d*Tt_yes/d';
  135. ascaled=mu*(eye(n)-T_sectf)/(eye(n)+T_sectf);
  136. x = x_yes;
  137. %
  138. % ------- End of MUOPT.M % PHL/MGS
  139.