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