|
|
|
|
Below is the function nyquist1.m. This function is a modified version of the nyquist command, and has all the same attributes as the original, with a few improvements. Nyquist1.m takes poles on the imaginary axis into account when creating the Nyquist plot, and plots around them. Furthermore, this function outputs the number of open-loop right hand plane poles, the number of anti-clockwise encirclements, and the number of closed-loop right half plane poles onto your screen.
Copy the following text into a file nyquist.m, and put it in the same directory as the Matlab software, or in a directory which is contained in Matlab's search path.
function [reout,imt,w] = mnyquist(a,b,c,d,iu,w)
%NYQUIST1 Nyquist frequency response for continuous-time linear systems.
%
% This Verison of the NYQUIST Command takes into account poles at the
% jw-axis and loops around them when creating the frequency vector in order
% to produce the appropriate Nyquist Diagram (The NYQUIST command does
% not do this and therefore produces an incorrect plot when we have poles in the
% jw axis). As an added feature, this function outputs the number of open loop right
% hand plane poles, the number of anti-clockwise encirclements, and the number of
% closed loop right half plane poles on your screen.
%
% NOTE: This version of NYQUIST1 does not account for pole-zero
% cancellation. Therefore, the user must simplify the transfer function before using
% this command.
%
% NYQUIST1(A,B,C,D,IU) produces a Nyquist plot from the single input
% IU to all the outputs of the system:
% . -1
% x = Ax + Bu G(s) = C(sI-A) B + D
% y = Cx + Du RE(w) = real(G(jw)), IM(w) = imag(G(jw))
%
% The frequency range and number of points are chosen automatically.
%
% NYQUIST1(NUM,DEN) produces the Nyquist plot for the polynomial
% transfer function G(s) = NUM(s)/DEN(s) where NUM and DEN contain
% the polynomial coefficients in descending powers of s.
%
% NYQUIST1(A,B,C,D,IU,W) or NYQUIST1(NUM,DEN,W) uses the user-supplied
% freq. vector W which must contain the frequencies, in radians/sec,
% at which the Nyquist response is to be evaluated. When invoked
% with left hand arguments,
% [RE,IM,W] = NYQUIST1(A,B,C,D,...)
% [RE,IM,W] = NYQUIST1(NUM,DEN,...)
% returns the frequency vector W and matrices RE and IM with as many
% columns as outputs and length(W) rows. No plot is drawn on the
% screen.
%
% Modified by Mithran Mathew at Carnegie Mellon University
% under the supervision of Dr. William Messner.
%
%
% See also: LOGSPACE,MARGIN,BODE, and NICHOLS.
% J.N. Little 10-11-85
% Revised ACWG 8-15-89, CMT 7-9-90, ACWG 2-12-91, 6-21-92,
% AFP 2-23-93
% Copyright (c) 1986-93 by the MathWorks, Inc.
%
%********************************************************************
% Modifications made to the nyquist - takes into account poles on jw axis.
% then goes around these to make up frequency vector
if nargin==0, eval('exresp(''nyquist'')'), return, end
% --- Determine which syntax is being used ---
if (nargin==1),
error('Wrong number of input arguments.');
elseif (nargin==2), % Transfer function form without frequency vector
num = a; den = b;
w = freqint2(num,den,30);
[ny,nn] = size(num); nu = 1;
elseif (nargin==3), % Transfer function form with frequency vector
num = a; den = b;
w = c;
[ny,nn] = size(num); nu = 1;
elseif (nargin==4), % State space system, w/o iu or frequency vector
error(abcdchk(a,b,c,d));
w = freqint2(a,b,c,d,30);
[iu,nargin,re,im]=mulresp('nyquist',a,b,c,d,w,nargout,0);
if ~iu, if nargout, reout = re; end, return, end
[ny,nu] = size(d);
elseif (nargin==5), % State space system, with iu but w/o freq. vector
error(abcdchk(a,b,c,d));
w = freqint2(a,b,c,d,30);
[ny,nu] = size(d);
else
error(abcdchk(a,b,c,d));
[ny,nu] = size(d);
end
if nu*ny==0, im=[]; w=[]; if nargout~=0, reout=[]; end, return, end
% *********************************************************************
% depart from the regular nyquist program here
% now we have a frequency vector, a numerator and denominator
% now we create code to go around all poles and zeroes here.
if (nargin==5) | (nargin ==4) | (nargin == 6)
[num,den]=ss2tf(a,b,c,d)
end
tol = 1e-6; %defined tolerance for finding imaginary poles
z = roots(num);
p = roots(den);
zp = [z;p]; % combine the zeros and poles of the system
nzp = length(zp); % number of zeros and poles
ones_zp=ones(nzp,1);
%index poles with zero real part + non-neg imag part
Ipo = find((abs(real(p))=0));
% **** only if we have such poles do we do the following:****
if ~isempty(Ipo)
po = p(Ipo); % poles with 0 real part and non-negative imag part
% check for distinct poles
[y,ipo] = sort(imag(po)); % sort imaginary parts
po = po(ipo);
dpo = diff(po);
idpo = find(abs(dpo)>tol);
idpo = [1;idpo+1]; % indexes of the distinct poles
po = po(idpo); % only distinct poles are used
nIpo = length(idpo); % # of such poles
originflag = find(imag(po)==0); % locate origin pole
s = []; % s is our frequency response vector
w = sqrt(-1)*w; % create a jw vector to evaluate all frequencies with
for ii=1:nIpo % for all Ipo poles
[nrows,ncolumns]=size(w);
if nrows == 1
w = w.'; % if w is a row, make it a column
end;
if nIpo == 1
r(ii) = .1;
else % check distances to other poles and zeroes
pdiff = zp-po(ii)*ones_zp % find the differences between
% poles you are checking and other poles and zeros
ipdiff = find(abs(pdiff)>tol) % ipdiff is all nonzero differences
r(ii)=0.2*min(abs(pdiff(ipdiff))) % take half this difference
r(ii)=min(r(ii),0.1) % take the minimum of this diff.and .1
end;
t = linspace(-pi/2,pi/2,25);
if (ii == originflag)
t = linspace(0,pi/2,25);
end;
% gives us a vector of points around each Ipo
s1 = po(ii)+r(ii)*(cos(t)+sqrt(-1)*sin(t)); % detour here
s1 = s1.'; % make sure it is a column
% Now here I reconstitute s - complex frequency - and
% evaluate again.
bottomvalue = po(ii)-sqrt(-1)*r(ii); % take magnitude of imag part
if (ii == originflag) % if this is an origin point
bottomvalue = 0;
end;
topvalue = po(ii)+sqrt(-1)*r(ii); % the top value where detour stops
nbegin = find(imag(w)imag(topvalue)); % find all points greater than
nnend = length(nend); % encirclement around jw root
if (nnend == 0)
send = 0
else send = w(nend);
end
w = [sbegin; s1; send]; % reconstituted half of jw axis
end
else
w = sqrt(-1)*w;
end
% this ends the loop for imaginary axis poles
% *************************************************************
% back to the regular nyquist program here
% Compute frequency response
if (nargin==2)|(nargin==3)
gt = freqresp(num,den,w);
else
gt = freqresp(a,b,c,d,iu,w);
end
ret=real(gt);
imt=imag(gt);
% If no left hand arguments then plot graph.
if nargout==0,
status = ishold;
plot(ret,imt,'r-',ret,-imt,'g--')
% plot(real(w),imag(w))
% modifications added here to count encirclements
[numc,denc] = tfchk(num,den);
gw = [(ret + j*imt); numc(1)/denc(1); (flipud(ret) - j*flipud(imt))]; % look at G(jw)
[Ngw,Mgw] = size(gw); % size of evaluated G(jw)
gwp1 = gw + ones(Ngw,Mgw);
% moves origin from 0 to -1, so we can count encirclements around -1 now
initial_angle = rem(180/pi*angle(gwp1(1)) + 360, 360); % define initial angle
angle_gwp1 = rem(180/pi*angle(gwp1) - initial_angle + 720,360); % angle from origin for all points
dagw = diff(angle_gwp1); % subtract off initial angle find degrees of encirclement
tolerance = 180; % define tolerance - where encirclement counter "flips" over
Ipd = find(dagw < -tolerance); % number of anti-clockwise encirclements
Ind = find(dagw > tolerance); % number of clockwise encirclemtents
Nacw = max(size(Ipd)) - max(size(Ind)); % number of encirclements
% Nyquist Criterion says Z = P - N
P = length(find(p>0))
disp('P = number of Open loop poles in rhp');
N = Nacw
disp('N = number of anti-clockwise encirclements')
Z = P - N %
disp('Z = number of closed loop poles in rhp')
set(gca, 'YLimMode', 'auto')
limits = axis;
% Set axis hold on because next plot command may rescale
set(gca, 'YLimMode', 'auto')
set(gca, 'XLimMode', 'manual')
hold on
% Make arrows
for k=1:size(gt,2),
g = gt(:,k);
re = ret(:,k);
im = imt(:,k);
sx = limits(2) - limits(1);
[sy,sample]=max(abs(2*im));
arrow=[-1;0;-1] + 0.75*sqrt(-1)*[1;0;-1];
sample=sample+(sample==1);
reim=diag(g(sample,:));
d=diag(g(sample+1,:)-g(sample-1,:));
% Rotate arrow taking into account scaling factors sx and sy
d = real(d)*sy + sqrt(-1)*imag(d)*sx;
rot=d./abs(d); % Use this when arrow is not horizontal
arrow = ones(3,1)*rot'.*arrow;
scalex = (max(real(arrow)) - min(real(arrow)))*sx/50;
scaley = (max(imag(arrow)) - min(imag(arrow)))*sy/50;
arrow = real(arrow)*scalex + sqrt(-1)*imag(arrow)*scaley;
xy =ones(3,1)*reim' + arrow;
xy2=ones(3,1)*reim' - arrow;
[m,n]=size(g);
hold on
plot(real(xy),-imag(xy),'r-',real(xy2),imag(xy2),'g-')
end
xlabel('Real Axis'), ylabel('Imag Axis')
limits = axis;
% Make cross at s = -1 + j0, i.e the -1 point
if limits(2) >= -1.5 & limits(1) <= -0.5 % Only plot if -1 point is not far out.
line1 = (limits(2)-limits(1))/50;
line2 = (limits(4)-limits(3))/50;
plot([-1+line1, -1-line1], [0,0], 'w-', [-1, -1], [line2, -line2], 'w-')
end
% Axis
plot([limits(1:2);0,0]',[0,0;limits(3:4)]','w:');
if ~status, hold off,end % Return hold to previous status
return % Suppress output
end
%reout = ret;
%plot(real(p),imag(p),'x',real(z),imag(z),'o');
Use your browser's "Back" button to return to the previous page.
8/29/96 LO