Matlab hacks
This is a collection of modified Matlab files, with slight
modifications to make them work correctly for Problem 6 of
Homework 9.
rlocus.m
function [rout,kout] = rlocus(a,b,c,d,k)
%RLOCUS Evans root locus.
%
% For a SISO system SYS with numerator N(s) and denominator D(s),
% RLOCUS(SYS) calculates and plots the locus of the roots of
%
% H(s) = D(s) + k * N(s) = 0
%
% for a set of positive gains K which are adaptively calculated to
% produce a smooth plot.
%
% RLOCUS(SYS,K) uses a user-specified vector K of gains.
%
% [R,K] = RLOCUS(SYS) or R = RLOCUS(SYS,K) returns the matrix R
% of complex root locations for the gains K. R has LENGTH(K) columns
% and its j-th column lists the closed-loop roots for the gain K(j).
%
% See also RLOCFIND, PZMAP, EIG.
% Old help
%RLOCUS Evans root locus.
% RLOCUS(NUM,DEN) calculates and plots the locus of the roots of:
%
% H(s) = DEN(s) + k*NUM(s) = 0
%
% for a set of gains K which are adaptively calculated to produce a
% smooth plot. Alternatively the vector K can be specified with an
% optional right-hand argument RLOCUS(NUM,DEN,K). Vectors NUM and
% DEN must contain the numerator and denominator coefficients in
% descending powers of s or z. When invoked with left hand arguments
% R = RLOCUS(NUM,DEN,K) or [R,K] = RLOCUS(NUM,DEN)
% returns the matrix R with LENGTH(K) rows and LENGTH(DEN)-1
% columns containing the complex root locations. Each row of the
% matrix corresponds to a gain from vector K. If a second left hand
% argument is included, the gains are returned in K.
%
% RLOCUS(A,B,C,D), R = RLOCUS(A,B,C,D,K), or [R,K] = RLOCUS(A,B,C,D)
% finds the root-locus from the equivalent SISO state-space system
% (A,B,C,D).
% dx/dt = Ax + Bu u = -k*y
% y = Cx + Du
%
% See also: PZMAP and RLOCTOOL.
% J.N. Little 10-11-85
% Revised A.C.W.Grace 7-8-89, 6-21-92
% Revised A.Potvin 6-1-94
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 1.5 $ $Date: 1997/12/01 22:15:25 $
ni = nargin;
no = nargout;
if ni==0,
if no~=0, error('Missing input argument(s).'), end
eval('exresp(''rlocus'');')
return
end
error(nargchk(2,5,ni));
if ni<4,
% Transfer function
if ni==3,
k = c;
else
k = [];
end
[num,den] = tfchk(a,b);
p = size(num,1);
if p>1,
disp('Must call rlocus with a SISO system. e.g. rlocus(num(i,:),den)')
return
end
sys = tf(num,den);
else
% State space
[msg,a,b,c,d] = abcdchk(a,b,c,d); error(msg);
if ni==4,
k = [];
end
[p,m] = size(d);
if p*m>1,
disp('Must call rlocus with a SISO system. e.g. rlocus(a,b(:,i),c(j,:),d(i,j))')
return
end
sys = ss(a,b,c,d);
end
if no,
[rout,kout] = rlocus2(sys,k);
% Make output consistent with old syntax
rout = rout.'; % length(k) x length(Poles)
else
rlocus2(sys,k)
end
% end rlocus
rlocus2.m
function [rout,kout] = rlocus(sys,k)
%RLOCUS Evans root locus.
%
% For a SISO system SYS with numerator N(s) and denominator D(s),
% RLOCUS(SYS) calculates and plots the locus of the roots of
%
% H(s) = D(s) + k * N(s) = 0
%
% for a set of positive gains K which are adaptively calculated to
% produce a smooth plot.
%
% RLOCUS(SYS,K) uses a user-specified vector K of gains.
%
% [R,K] = RLOCUS(SYS) or R = RLOCUS(SYS,K) returns the matrix R
% of complex root locations for the gains K. R has LENGTH(K) columns
% and its j-th column lists the closed-loop roots for the gain K(j).
%
% See also RLOCFIND, PZMAP, EIG.
% J.N. Little 10-11-85
% Revised A.C.W.Grace 7-8-89, 6-21-92
% Revised A.Potvin 6-1-94, P. Gahinet 7-19-96
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 1.9 $ $Date: 1997/12/01 22:05:15 $
ni = nargin;
no = nargout;
if ni==0,
if no~=0, error('Missing input argument(s).'), end
eval('exresp(''rlocus'');')
return
elseif ni==1
k = [];
end
error(nargchk(1,2,ni));
% Extract ZPK data
[Zeros,Poles,Gain,Ts,Td] = zpkdata(sys,'v');
if prod(size(Gain))>1,
error('System SYS must be a SISO system.')
elseif Td>0,
error('RLOCUS does not support delay systems.')
end
% Quick exit for empty models
if isempty(Gain),
% Empty model: exit
if no~=0, rout = []; kout = []; end
return
end
% For locus computation, convert model to SS if proper
% and to TF otherwise
if length(Zeros)<=length(Poles),
sys = ss(sys);
else
sys = tf(sys);
end
% Generate root locus
kgiven = ~isempty(k);
if kgiven,
% When k is given: make sure it's a row vector
k = k(:).';
% Remember that genrloc returns r which is length(Poles) by length(k)
r = genrloc(sys,Zeros,Poles,k);
r = sortrloc(r);
elseif Gain==0,
warning('System has zero gain')
r = zeros(0,1); k = zeros(1,0);
elseif isempty(Zeros) & isempty(Poles),
warning('System is a static gain')
r = zeros(0,1); k = zeros(1,0);
else
% Adaptively search for values of gain if they are not specified
[k,r] = gainrloc(sys,Zeros,Poles,Gain);
end
% Draw plot
if no==0,
% Find the "right" axis,
ax = newplot;
OldNextPlot = get(ax,'NextPlot');
if Gain==0 | (isempty(Zeros) & isempty(Poles))
pzmap(sys);
return
end
%-Get new xlimits
newxmin=min([-0.5, min([real(Zeros);real(Poles)])]);
newxmax=max([0.5, max([real(Zeros);real(Poles)])]);
set(ax,'Box','on','Nextplot','add')
% Make X and Y axes, Poles, and Zeros
PlantPoleH = plot(real(Poles),imag(Poles),'x');
PlantZeroH = plot(real(Zeros),imag(Zeros),'o');
% Plot root locus data in "right" axis
rplot=r.';
LocusH = line(real(rplot),imag(rplot));
% Shrink axis limits, if necessary
[valmin,indmin]=min(real(rplot));
[valmax,indmax]=max(real(rplot));
[rr,rc] = size(r);
belowXmin=find(valminnewxmax);
aboveflag=0;
xminnew=newxmin;
xmaxnew=newxmax;
MaxPole=max([1;real(Zeros);real(Poles)]);
tol=1e7;
for ctmin = 1:length(belowXmin),
if indmin(belowXmin(ctmin)) < rc & ( ~(valmin(belowXmin(ctmin)) < -MaxPole*tol ) ),
onedown = max([1,indmin(belowXmin(ctmin))-1]);
oneup = min([rc,indmin(belowXmin(ctmin))+1]);
if ( ~(rplot(onedown,belowXmin(ctmin)) > MaxPole*tol) ) & ...
( ~(rplot(oneup,belowXmin(ctmin)) > MaxPole*tol) ),;
xminnew=min([xminnew, valmin(belowXmin(ctmin))]);
else
belowflag=1;
end
else
belowflag=1;
end
end
if belowflag
xminnew = (xminnew+newxmax)/2 - (1.5*((newxmax-xminnew)/2));
else
xminnew = xminnew*1.05;
end
for ctmax = 1:length(aboveXmax),
if indmax(aboveXmax(ctmax)) < (rc-1) & ( ~(valmax(aboveXmax(ctmax)) > MaxPole*tol ) ),
xmaxnew = max([xmaxnew, valmax(aboveXmax(ctmax))]);
else
aboveflag=1;
end
end
if aboveflag
xmaxnew = (xmaxnew+newxmin)/2 + (1.5*((xmaxnew-newxmin)/2));
else
xmaxnew = xmaxnew*1.05;
end
set(gca,'xlim',[floor(xminnew),ceil(xmaxnew)]);
% Plot the XY axis,
xlim=[floor(xminnew),ceil(xmaxnew)];
linexlim=[min([xlim(1),-1]); max([xlim(2),1])];
ylim=get(gca,'Ylim');
lineylim=[min([ylim(1),-1]); max([ylim(2),1])];
XYaxisH = line([linexlim [0; 0]],[[0; 0] lineylim],'color','k','linestyle',':');
if kgiven,
% Remember that when k is given, only plot crosses
set(LocusH,'LineStyle','none','Marker','x');
end
xlabel('Real Axis')
ylabel('Imag Axis')
set(ax,'NextPlot',OldNextPlot);
else
rout = r;
kout = k;
end
% end @lti/rlocus
genrloc.m
function r = genrloc(sys,z,p,k)
%GENRLOC Generates points along root locus
%
% R = GENRLOC(SYS,Z,P,K) computes the poles R of the
% negative feedback loop
%
% ---->o---->| SYS |---+--->
% | |
% +<-----| G |----+
%
% for the values of the gain G specified in the vector K.
% The vectors Z and P contain the zeros and poles of SYS
% and the matrix R is N-by-length(K) where N is the number
% of closed-loop poles.
%
% See also ROCUS and RLOCFIND.
% Author(s): A. Potvin, 12-1-93, PG 7-9-97
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 1.4 $ $Date: 1997/12/01 22:06:09 $
% Extract transfer function data
[num,den] = tfdata(sys,'v');
n = max(length(num),length(den))-1;
% Pre-allocate space for R
lk = length(k);
r = zeros(n,lk);
r(:) = Inf;
% For all k, determine the roots of den+k*num
for i=1:(n>0)*lk,
ki = k(i);
if ~isfinite(ki),
clr = z;
elseif ki==0,
clr = p;
else
clr = roots(den+k(i)*num);
end
nclr = length(clr);
r(n-nclr+1:n,i) = clr;
end
% end tf/genrloc
A few other files, namely
- axesrloc.m
- matchlsq.m
- zpkdata.m
- gainrloc.m
- sortrloc.m
had to be comped to the working directory for this to work.