home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / OSBORNE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.3 KB  |  173 lines
  1. function [mu, Ascaled, logd] = osborne(A,K)
  2. %
  3. % [MU,ASCALED,LOGD] = OSBORNE(A) 
  4. % [MU,ASCALED,LOGD] = OSBORNE(A,K) produces the scalar upper bound MU on
  5. % the structured singular value (ssv) computed via
  6. %           mu = max(svd(diag(dosb)*A/diag(dosb)))
  7. % where dosb obtained by applying the Osborne scaling technique.
  8. %
  9. % Input: 
  10. %     A  -- a pxq complex matrix whose ssv is to be computed
  11. % Optional input:
  12. %     K  -- an nx1 or nx2 matrix whose rows are the uncertainty block sizes  
  13. %           for which the ssv is to be evaluated; K must satisfy 
  14. %           sum(K) == [q,p].  If only the first column of K is given then the
  15. %           uncertainty blocks are taken 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 -- diag(dosb)*A/diag(dosb)
  19. %     LOGD    -- dosb=exp(LOGD)) is Osborne scaling vector
  20. %
  21.  
  22. % R. Y. Chiang & M. G. Safonov 8/15/88
  23. % Copyright (c) 1988-89 by the MathWorks, Inc.
  24. % All Rights Reserved.
  25. % --------------------------------------------------------------------
  26.  
  27. %
  28. % ---------------- Tolerances used in this routine -------------------
  29. %  If the matrix A is reducible, the Osborne theory does not apply 
  30. %  directly, but osborne.m compensates for this by computing the Osborne
  31. %  scaling for a slightly perturbed matrix having a slightly larger
  32. %  (i.e., slightly more conservative) SSV upper bound.  Reducibility
  33. %  is detected by examining the norm of row and column of the matrix that 
  34. %  forms the actual diagonal scaling. The following constants are used by 
  35. %  this routine; they may be changed if necessary:
  36.  
  37. tol = 1/1.e120;     % Value for variable "tol" used
  38. alfa0 = eps/10000;  % Initial Perturbation size
  39. maxiter = 20;       % Maximum allowable iterations of perturbation loop
  40. maxcounter = 10;    % Maximum allowable number of decreases in perturbation  
  41. % --------------------------------------------------------------------%- 
  42.  
  43. [p,q] = size(A);
  44. if nargin < 2,    
  45.    K = ones(q,2);   % set blocksize to default 1x1
  46. end  
  47. K = [K,K]; K=K(:,1:2);            % set K(:,2)=K(:,1) if K(:,2) is not present
  48. n = size(K)*[1;0];                % number of uncertainty blocks
  49. K2c = K(:,1);                  
  50. K1c = K(:,2);
  51. if sum(K1c)~= p | sum(K2c)~= q,
  52.     error('K and A are not dimensionally compatible')
  53. end
  54. %
  55. % ------ form matrix aa of max singular values of blocks of A
  56. onesvec = ones(n,1);
  57.  
  58. if n < p | n < q,
  59.    aa = zeros(n,n);
  60.    r=[0;conv(K1c,onesvec)];
  61.    c=[0;conv(K2c,onesvec)];
  62.    for i = 1:n
  63.       for j = 1:n
  64.           aa(i,j) = max(svd(A(r(i)+1:r(i+1),c(j)+1:c(j+1))));
  65.       end
  66.    end
  67. else
  68.    aa = abs(A);
  69. end
  70. %
  71. % ----- Initialize constants and temporary variables------
  72. % -------------used inside while loop below --------------
  73. aaa   = aa;
  74. alfa0;              % Initial perturbation of aa
  75. iter  = 0;          % Counter for while loop
  76. counter = 0;        % Counter for perturbation decreases in loop
  77. alfa = 0; lamp = 0; % This results in no perturbation in iteration 1
  78. logd=zeros(n,1);
  79. Ascaled=A;
  80. d = eye(n);         % Osborne initial scaling
  81. osbmaxiter = 20;    % Maximum Osborne iteratoins
  82. %
  83. % ------- Compute "Ascaled" and Osborne scaling "D" --------------
  84. %
  85. loop = 1;       
  86. while loop
  87.    iter = iter+1;
  88.    if (iter==maxiter)|(counter==maxcounter), % If this is maxiter-th iteration,
  89.      loop=0;                            %    make it the last
  90.      converge=0;
  91.    end
  92.    reducible=0;                 % Assume aaa is not reducible and
  93.    perttoobig=0;                %   that perturb. is not too big   
  94.                                 %   until proven otherwise
  95.    aaa = aa+alfa*ones(n,n);     %   perturb the aa matrix
  96.    % Starting Osborne iteration
  97.    for osbiter = 1:osbmaxiter
  98.       for i = 1:n
  99.           aao = aaa - diag(diag(aaa));
  100.           offrow = norm(aao(i,:));
  101.           offcol = norm(aao(:,i));
  102.           if (offcol < tol) | (offrow < tol)
  103.              reducible = 1;
  104.           end
  105.           if ~reducible
  106.               dhead(i,i) = sqrt(offrow/offcol);
  107.               aaa(i,:) = aaa(i,:)/dhead(i,i);
  108.               aaa(:,i) = aaa(:,i)*dhead(i,i);
  109.               d(i,i) = d(i,i)/dhead(i,i);
  110.           end
  111.       end  % of For loop
  112.    end     % of Osborne iteration loop
  113.    lam = eig(aaa); % checking perturbation size
  114.    [maxlam,ilam] = max(real(lam));
  115.    if iter==1,
  116.      lamp=maxlam;               % Unperturbed perron e-value 
  117.      alfa=alfa0*lamp;           % Perturbation size to be used in iter 2,
  118.    end   
  119.    if iter>1,                   % If aaa is perturbed version of aa
  120.      if maxlam>1.0001*lamp      %    determine if pert. was too big
  121.         perttoobig=1;           %    to preserve perron e-value to
  122.      else                       %    within .01 percent
  123.         perttoobig=0;
  124.      end
  125.    end
  126.    if ~reducible             % if nonreducibility confirmed by 
  127.        d = diag(d)/d(1,1);   % nomalize D scaling
  128.        dinv = onesvec./d;
  129.        dbar = d*dinv';
  130.        logd = log(d) + logd;  % overwrite logd
  131.        aa = dbar.*aa;         % overwrite aa
  132.        if n==p & n==q,
  133.             Ascaled = dbar.*Ascaled;  %     and overwrite Ascaled.
  134.        else
  135.           for i = 1:n
  136.              for j = 1:n
  137.                  Ascaled(r(i)+1:r(i+1),c(j)+1:c(j+1)) ...
  138.                = dbar(i,j)*Ascaled(r(i)+1:r(i+1),c(j)+1:c(j+1));
  139.              end
  140.           end       
  141.        end
  142.        if ~perttoobig   % If everything is ok
  143.            converge=1;  % then convergence achieved;
  144.            loop=0;      % exit while loop
  145.        else             % else decrease perturbation:
  146.            counter=counter+1;
  147.            alfa=.01*alfa;           
  148.        end
  149.    else % reducible     % else increase perturbation
  150.      alfa=10000*alfa;
  151.    end
  152. end        % of perturbation loop
  153. %
  154. % -------- Compute the osborne singular value muosv -----------
  155. %
  156. musv  = max(svd(A));
  157. if converge
  158.     muosv = max(svd(Ascaled));
  159.     if muosv > musv,          % make sure we never do worse than no scaling!
  160.       logd = zeros(n,1);
  161.       Ascaled = A;
  162.       muosv = musv;
  163.     end
  164. else                          % if numerical problem
  165.     muosv = musv; 
  166.     if nargout > 1
  167.        disp('Warning OSBORNE.M: ASCALED, LOGD values are inaccurate!');
  168.     end
  169. end
  170. mu = muosv;
  171. %
  172. % ------ End of OSBORNE.M --- RYC/MGS 12/26/90 %
  173.