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