home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / SIGMA.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.6 KB  |  139 lines
  1. function [svout,w] = sigma(Z1,Z2,Z3,Z4,Z5,Z6)
  2. %SIGMA Singular value frequency response of continuous linear systems.
  3. %    SIGMA(A,B,C,D) or SIGMA(SS_) produces a singular value plot of 
  4. %       the complex matrix      -1
  5. %                G(jw) = C(jwI-A) B + D 
  6. %    as a function of frequency.  The singular values are an extension
  7. %    of Bode magnitude response for MIMO systems.  The frequency range
  8. %    and number of points are chosen automatically.  For square 
  9. %    systems, SIGMA(A,B,C,D,'inv') produces the singular values of the
  10. %    inverse complex matrix
  11. %         -1               -1      -1
  12. %           G (jw) = [ C(jwI-A) B + D ]
  13. %
  14. %    SIGMA(A,B,C,D,W) or SIGMA(A,B,C,D,W,'inv') uses the user-supplied
  15. %    frequency vector W which must contain the frequencies, in 
  16. %    radians/sec, at which the singular value response is to be 
  17. %    evaluated. When invoked with left hand arguments,
  18. %        [SV,W] = SIGMA(A,B,C,D,...) or [SV,W] = SIGMA(SS_,...)
  19. %    returns the frequency vector W and the matrix SV with as many 
  20. %    columns    as MIN(NU,NY) and length(W) rows, where NU is the number
  21. %    of inputs and NY is the number of outputs.  No plot is drawn on 
  22. %    the screen.  The singular values are returned in descending order.
  23. %
  24.  
  25. %    See also: LOGSPACE,SEMILOGX,NICHOLS,NYQUIST and BODE.
  26.  
  27. %    Andy Grace  7-10-90, 
  28. %    R. Chiang   Modified to include System data type.
  29. %    Revised W. Wang 7-22-92
  30. %    copyright (c) 1990 by the MathWorks, Inc.
  31.  
  32. if nargin==0, eval('exresp(''sigma'',1)'), return, end
  33.  
  34. if exist('mkargs') == 2, %If RCT installed
  35.   inargs='(a,b,c,d,w,invflag)';
  36.   eval(mkargs(inargs,nargin,'ss'))
  37. else
  38.    if nargin<4,
  39.      error('Too few input arguments')
  40.    else
  41.      a=Z1; b=Z2; c=Z3; d=Z4;
  42.      if nargin > 5
  43.        invflag = Z6;
  44.      elseif nargin > 4,
  45.        w=Z5;
  46.      end;
  47.    end;
  48. end;
  49.  
  50. if nargin==6,      % Trap call to RCT function
  51.   if ~isstr(invflag),
  52.     eval('svout = sigma2(a,b,c,d,w,invflag);')
  53.     return
  54.   end
  55.   if ~length(invflag)
  56.     nargin = nargin - 1;
  57.   end
  58. end
  59.  
  60. error(nargchk(4,6,nargin));
  61. error(abcdchk(a,b,c,d));
  62.  
  63. % Detect null systems
  64. if ~(length(d) | (length(b) &  length(c)))
  65.        return;
  66. end
  67.  
  68. % Determine status of invflag
  69. if nargin==4, 
  70.   invflag = [];
  71.   w = [];
  72. elseif (nargin==5)
  73.   if (isstr(w)),
  74.     invflag = w;
  75.     w = [];
  76.     [ny,nu] = size(d);
  77.     if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
  78.   else
  79.     invflag = [];
  80.   end
  81. else
  82.   [ny,nu] = size(d);
  83.     if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
  84. end
  85.  
  86. % Generate frequency range if one is not specified.
  87.  
  88. % If frequency vector supplied then use Auto-selection algorithm
  89. % Fifth argument determines precision of the plot.
  90. if ~length(w)
  91.   w=freqint(a,b,c,d,30);
  92. end
  93.  
  94. [nx,na] = size(a);
  95. [no,ns] = size(c);
  96. nw = max(size(w));
  97.  
  98. % Balance A
  99. [t,a] = balance(a);
  100. b = t \ b;
  101. c = c * t;
  102.  
  103. % Reduce A to Hessenberg form:
  104. [p,a] = hess(a);
  105.  
  106. % Apply similarity transformations from Hessenberg
  107. % reduction to B and C:
  108. b = p' * b;
  109. c = c * p;
  110.  
  111. s = w * sqrt(-1);
  112. I=eye(length(a));
  113. if nx > 0,
  114.   for i=1:length(s)
  115.     if isempty(invflag),
  116.       sv(:,i)=svd(c*((s(i)*I-a)\b) + d);
  117.     else
  118.       sv(:,i)=svd(inv(c*((s(i)*I-a)\b) + d));
  119.     end
  120.   end
  121. else
  122.   for i=1:length(s)
  123.     if isempty(invflag),
  124.       sv(:,i)=svd(d);
  125.     else
  126.       sv(:,i)=svd(inv(d));
  127.     end
  128.   end
  129. end
  130.   
  131. % If no left hand arguments then plot graph.
  132. if nargout==0
  133.   semilogx(w,20*log10(sv),w,zeros(1,length(w)),'w:')
  134.   xlabel('Frequency (rad/sec)')
  135.   ylabel('Singular Values dB')
  136.   return % Suppress output
  137. end
  138. svout = sv; 
  139.