FUNCTION LISTINGS
See Appendix 14 of
Classical Feedback Control
for descriptions of these MATLAB funtions.
The functions can be copied freely for nonprofit use.
The author bears no responsibility for the results of the functions'
applications.
Download BodeStepToolbox.tar
(MATLAB and SPICE scripts for the examples and problems in
Classical Feedback Control are available
here.)
NYQLOG
function NYQLOG(wb,numl,denl)
%NYQLOG Nyquist diagram on logarithmic plane
%function NYQLOG(wb,numl,denl)
% The loop transfer function is numl(s)/denl(s), where numl and
% denl contain the polynomial coefficients in descending powers of s.
% wb is the crossover frequency in rad/sec, i.e. frequency at which
% loop gain is 0 dB.
% Default: wb = 1, numl = [-26 178 654 1269 429],
% denl = [1 18.8 132 543 1164 860 0 0]
%
% Copyright (c) B. J. Lurie 4-1-98, La Crescenta.
% See also BOSTEP, BOCLOS, BONYQAS, BOINTEGR, TFSHIFT, STEPLOG from the
% Bode toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 3,
denl = [1 18.8 132 543 1164 860 0 0];
end
if nargin < 2,
numl = [-26 178 654 1269 429];
end
if nargin == 0,
wb = 1;
end
[mag, phase] = bode(numl, denl); plot(phase, 20*log10(mag),'w')
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
axis([-270 -90 -20 70])
hold on
[a,b] = bode (numl,denl,wb); plot(b,20*log10(a),'wx')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode(numl,denl,wb*w); plot(b,20*log10(a),'w+');
end
plot(-180, 0, 'wo');
hold off
title('Nyquist diagram, x marks w = wb, + marks octaves')
xlabel('loop phase shift in degrees')
ylabel('loop gain in dB')
zoom on
Back
BOLAGNYQ
function BOLAGNYQ(wb,numl,denl)
%BOLAGNYQ Gain and -(phase+180) (i.e. lag margin) responses, and Nyquist
% diagram on logarithmic plane
%function BOLAGNYQ(wb,numl,denl)
% The loop transfer function is numl(s)/denl(s), where numl and
% denl contain the polynomial coefficients in descending powers of s.
% Frequency wb in rad/sec specifies the position of cursor x; it is
% convenient to put it at the frequency where the loop gain is 0 dB.
%
% Default/demo: wb = 1, numl = [-26 178 654 1269 429],
% denl = [1 18.8 132 543 1164 860 0 0]
%
% Copyright (c) B. J. Lurie 6-16-98, La Crescenta.
% See also BOSTEP, BOCLOS, BONYQAS, BOINTEGR, TFSHIFT, STEPLOG, NYQLOG from
% the Bodestep toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 3,
denl = [1 18.8 132 543 1164 860 0 0];
end
if nargin < 2,
numl = [-26 178 654 1269 429];
end
if nargin == 0,
wb = 1;
end
[mag, phase, w] = bode(numl, denl);
subplot(2,1,1)
semilogx(w, 20*log10(mag),'w', w,(-180-phase),'w--');
hold on
[mmb,ppb] = bode(numl, denl,wb);
semilogx(wb,mmb,'wx');
wmin = min(w); wmax = max(w);
axis([wmin wmax -50 100])
xlabel('rad/sec')
ylabel('gain, dB; lag margin, degr')
grid
hold off
subplot(2,1,2)
plot(phase, 20*log10(mag),'w')
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
axis([-270 -90 -20 50])
hold on
[a,b] = bode (numl,denl,wb); plot(b,20*log10(a),'wx')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode(numl,denl,wb*w); plot(b,20*log10(a),'w+');
end
plot(-180, 0, 'wo');
hold off
xlabel('phase, degr; x marks wb, + mark octaves from wb')
ylabel('loop gain, dB')
zoom on
wb
gain_at_wb = mmb
lag_margin = 180 + ppb
Back
TFSHIFT
function [numsh, densh] = tfshift(num,den,q,p)
%TFSHIFT Transfer function with shifted frequency response.
% [numsh, densh] = tfshift(num,den,q,p)
% TFSHIFT produces numerator and denominator of transfer function
% NUMSH(s)/DENSH(s) shifted in frequency q times from NUM(s)/DEN(s)
% where NUM and DEN contain the polynomial coefficients in
% descending powers of s. NUMSH(s)/DENSH(s) has the same value
% at frequency qw as the function NUM(s)/DEN(s) at frequency w.
% The numerator and denumerator are both multiplied by 10^p for
% the coefficients to have convenient values. Default is p = 0.
%
% Copyright (c) 1998 A. Ahmed and B. J. Lurie 2-27-97.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
ln = length(num)-1;
ld = length(den)-1;
pn = [ln: -1: 0];
pd = [ld: -1: 0];
numsh = num./(q.^pn);
densh = den./(q.^pd);
if nargin == 4,
numsh = numsh*10^p;
densh = densh*10^p;
end
Back
BONYQAS
function [gain1,phase1] = bonyqas(u,g,lslope,hslope,nml);
% function [gain1,phase1] = bonyqas(u,g,lslope,hslope,bnm);
% BONYQAS Asymptotic Bode and Nyquist diagrams corresponding in minimum
% phase manner to a piece-linear gain response in dB specified by the
% gain g(h) at k corner frequencies u(h) in rad/sec. The lengths of
% vectors g and u must be equal. Coefficient lslope defines the
% asymptotic slope at lower frequencies in integrator units, and
% hslope defines the asymptotic slope at higher frequencies. nml is
% nonminimum phase lag in degrees at u(h); the lag is proportional
% to the frequency, default is nml = 0.
% Smoothed Bode diagram is plotted in white. Phase is plotted in green.
% Nyquist diagram is plotted on logarithmic plane.
% The function is applicable for reasonably smooth responses. The phase
% accuracy for typical feedback loop shaping problem without sharp
% peaks or notches is better than 3 degr.
%
% Default: bonyqas without arguments plots the band-pass response:
% bonyqas([0.2 2 8 32],[-5 10 6 -12],-1,3);
% Low-pass loop Example 1: bonyqas([3 4 20 30],[16 13 -9 -10],2, 3, 0)
% Loop Example 2: bonyqas([3 3.8 22 45],[16 13 -10 -11.5],2, 3, 50)
%
% Copyright (c) B. Lurie 6-13-98.
% See also BOSTEP, BOCLOS, BOCOMP, BOINTEGR, NYQLOG, TFSHIFT from the
% Bode toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
% default
if nargin == 0,
u = [0.2 2 8 32]; % 4 corner frequencies
g = [-5 10 6 -12]; % gain in dB; band-pass system
lslope = -1; % -1 integrator, i.e. 6.02 dB/oct up-slope
hslope = 3; % 3 integrators, i.e. -18.06 dB/oct
end
if nargin < 5,
nml = 0;
end
k = length(u); h = [1:k];
% frequency range of calculations and plots
wmin = floor(log10(u(1)) - 0.33);
wmax = ceil(log10(u(k)) + 0.33);
w = logspace(wmin,wmax);
% low-frequency asymptote
gain1 = (-lslope * (20 * log10(w/u(1))) + g(1))';
phase1 = -lslope * 90;
% phase1 = -lslope * 90 * ones(length(w),1); % can be used instead
% semilogx(w, gain1, 'r--') % can be activated for purpose of teaching
% hold on
% segment slopes; segment no.h ends at u(h)
for h = [2:k]
sslope(h) = (g(h) - g(h-1))/log2(u(h)/u(h-1));
end % [0 15 -2 -9]
% sslope % can be activated for purpose of teaching
% ray slopes
rslope(1) = sslope(2) + 6.02 * lslope;
rslope(k) = - 6.02 * hslope - sslope(k);
for h = [2:(k-1)]
rslope(h) = sslope(h + 1) - sslope(h);
end
% rayslope_round = round(rslope) % [9 -17 -7 -9], positive means up
% frequency responses for k rays
nray = 1; dray = [1 1 1]; % normalized ray with damping coefficient 0.5;
for h = [1:k]
[nrayf, drayf] = lp2lp(nray,dray,u(h)); % specifying corner frequency
tfr = freqs(nrayf, drayf, w);
rgain = 20*real(log10(tfr));
rphase = (180/pi) * angle(tfr);
rasis = (-rslope(h)/12.041) * rgain'; % minus makes positive ray go up
% semilogx(w, rasis,'w--'); % can be activated for purpose of teaching
% semilogx(w, rphase,'g:'); % can be activated for purpose of teaching
gain1 = (-rslope(h)/12.041) * rgain' + gain1; % slope denormalization
phase1 = phase1 - (rslope(h)/12.041) * rphase';
end
% nonminimum phase lag addition
phase1 = phase1 - nml * (w/u(h))';
% Bode diagram and phase response
subplot(1,2,1)
semilogx(w,gain1,'w', w,phase1,'g--', u,g,'wo')
grid;
xlabel('frequency, rad/sec');
ylabel('gain, dB, and phase shift, degr');
title('Bode diagram')
hold off;
zoom on;
% Nyquist diagram on log plane
subplot(1,2,2)
plot(phase1, gain1,'w', -180,0,'wo')
hold on
gain2 = (-lslope * (20 * log10(u/u(1))) + g(1))'; % marks
phase2 = -lslope * 90;
for h = [1:k]
[nrayf, drayf] = lp2lp(nray,dray,u(h));
tfr = freqs(nrayf, drayf, u);
rgain = 20*real(log10(tfr));
rgain = real(rgain);
rphase = (180/pi) * angle(tfr);
gain2 = (-rslope(h)/12.041) * rgain' + gain2;
phase2 = phase2 - (rslope(h)/12.041) * rphase';
end
phase2 = phase2 - nml * (u/u(h))';
plot(phase2,gain2,'wo');
hold off;
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
axis([-270 -90 -20 70])
title('Nyquist diagram')
xlabel('phase shift, degr')
ylabel('gain in dB')
zoom on
Back
BOSTEP
function [wd,wc,width,numc,denc,numl,denl] = bostep(x,y,n,bnc,zetaz,zetap,typ)
% BOSTEP Normalized loop transfer function with Bode step.
% function [wd,wc,width,numc,denc,numl,denl] = bostep(x,y,n,bnc,zetaz,zetap,typ)
% Defaults: [wd,wc,width,numc,denc,numl,denl] = bostep(10,0.1667,3,0.1,0.6,0.4,2),
% x is gain stability margin in dB, default x = 10;
% y is phase stability margin y*pi rad, default y = 0.1667;
% n defines asymptotic slope -6n dB/oct, default n = 3;
% bnc (0 to 1) is n.p. shift in rad at the end of Bode step (at wc),
% default bnc = 1;
% zetaz is damping of zeros at beginning of Bode step (at wd),
% default zetaz = 0.6;
% zetap is damping of poles at the end of Bode step (at wc), default
% zetap = 0.4;
% typ (1 or 2) is the number of poles at the origin, default typ = 2.
% BOSTEP calculates the following parameters of the asymptotic Bode diagram
% withcrossover frequency 1 rad/sec: frequencies of beginning and end of
% the step wd, wc, and the width of the step width.
% BOSTEP also produces numerator and denominator of rational loop transfer
% function numl(s)/denl(s) approximating the asymptotic loop Bode diagram,
% where numl and denl contain the polynomial coefficients in descending
% powers of s. The typ is the number, 1 or 2, of the poles at the origin.
% BOSTEP also calculates the minimum phase compensator transfer function
% numc(s)/denc(s) for a single integrator plant with n.p. lag
% (-s + 2wc/bnc)/(s + 2wc/bnc).
%
% Copyright (c) B. J. Lurie 3-1-98, La Crescenta.
% See also BOCLOS, BONYQAS, BOCOMP, NYQLOG, TFSHIFT from the Bodestep toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7,
typ = 2;
end
if nargin < 6,
zetap = 0.4;
end
if nargin < 5,
zetaz = 0.6;
end
if nargin < 4,
bnc = 0.1;
end
if nargin < 3,
n = 3;
end
if nargin < 2,
y = 0.167;
end
if nargin == 0,
x = 10;
end
% ** Asymptotic diagram parameters
width = 0.9*(n + (pi/2)*bnc)/(2*(1 - y));
wd = 2^(x/12/(1 - y));
wc = wd*width;
% ** Poles and zeros for typ = 1
z1 = 0.5*(1-y);
z23 = [1 2*zetaz*wd wd^2]; % a pair of complex zeros
k = 0.35*wc^2*sqrt(0.17/y);
p1 = 0.07*(1-y); p2 = 0.7*wd; p3 = (wd + wc)/2.7;
p45 = [1 2*zetap*wc 0.81*wc^2]; ;% a pair of complex poles
% ** Poles and zeros for typ = 2
if typ == 2,
p1 = 0;
p3 = (wd + wc)/2.2;
z23 = 0.95*[1 2.2*zetaz*wd*sqrt(sqrt(y/0.17)) wd^2];
p45 = [1 2.2*zetap*wc 0.81*wc^2];
end
% ** Minimum phase compensator
numc = k * p3 * conv([1 z1], z23);
denc = conv(conv([1 p1],[1 p2]),conv([1 p3],p45));
% ** Loop with 1/s plant and n.p. lag (-s + 2wc/bnc)/(s + 2wc/bnc):
numl = conv(numc,[-1 2*wc/bnc]);
denl = conv(conv(denc,[1 0]),[1 2*wc/bnc]);
format short e
numl
denl
% numl_mp = numc; ne nada!
% denl_mp = conv(denc,[1 0]);ne nada!
% ** L-plane Nyquist diagram
[mag, phase] = bode(numl, denl); plot(phase, 20*log10(mag),'w')
set(gca,'XTick',[-240 -210 -180 -150 -120 -90])
grid
axis([-240 -90 -20 70])
hold on
nl = numl;
dl = denl;
[a,b] = bode (nl,dl,1); plot(b,20*log10(a),'wx')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode (nl,dl,w); plot(b,20*log10(a),'w+');
end
plot(-180, 0, 'wo')
hold off
title('Nyquist diagram, x marks w = 1, + marks octaves')
xlabel('loop phase shift in degrees')
ylabel('loop gain in dB')
zoom on
Back
BOCLOS
function [dclos, nprf1, dprf1, nprf2, dprf2, numtot, dentot] = boclos(numl,denl,y)
% BOCLOS Normalized closed loop transfer function without and with
% prefilter, and frequency response of 4th order Bessel filter. Also,
% step response, optional.
% function [dclos, nprf1, dprf1, nprf2, dprf2, numtot, dentot] =
% boclos(numl,denl,y)
% BOCLOS calculates denominator dclos of the closed loop transfer
% function without prefilter (the numerator is numl);
% prefilter1 and prefilter2 transfer functions nprf1(s)/dprf1(s)
% and nprf2(s)/dprf2(s); and closed loop transfer functions:
% numl(s)/dclos(s) without prefilter and numtot(s)/dentot(s)
% with prefilter; the coefficients are in descending powers of s.
% The defaults follow the default outputs of BOSTEP:
% numl = [-9.6082 626.39 1947.8 3326.0 1093]
% denl = [1 75.102 492.80 1716.3 3222.6 2217 0 0]; y = 0.1667;
%
% Copyright (c) B. J. Lurie 3-15-98.
% See also BOSTEP, BOCOMP, BONYQAS, NYQLOG, BOLAGNYQ, TFSHIFT from the
% Bodestep toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
% default:
if nargin < 3
y = 0.1667;
end
if nargin < 2
denl = [1 75.102 492.80 1716.3 3222.6 2217 0 0];
end
if nargin == 0
numl = [-9.6082 626.39 1947.8 3326.0 1093];
end
format short e
% Normalized closed loop function denominator (without prefilter)
if length(denl) - length(numl) == 2
dclos = [0 0 numl] + denl;
end
if length(denl) - length(numl) == 3
dclos = [0 0 0 numl] + denl;
end
if length(denl) - length(numl) == 4
dclos = [0 0 0 0 numl] + denl;
end
% First and second prefilter notches
nprf1 = [1 y/0.1667 1]; dprf1 = [1 2 1]; % prefilter 1
nprf2 = [1 0.47*y/0.1667 1]; dprf2 = [1 0.6 1]; % prefilter 2
% Prefilter 3 provides gain hump
nprf3 = [1 2.5 3.6]; dprf3 = [1 1.5 3.6]; % prefilter 3
% Closed loop transfer function with prefilters 2 and 3
numtot = conv(conv(nprf3,numl),conv(nprf1,nprf2));
dentot = conv(conv(dprf3,dclos),conv(dprf1,dprf2));
% Closed loop function with extra 3rd order Bessel lowpass
% prefilter,w = 0.25
numtot = conv([15],numtot);
dentot = conv([64 96 60 15],dentot);
% Open- and closed-loop and Bessel filter frequency responses
w = logspace(-1,1);
[magop,phaseop] = bode(numl, denl, w) % open-loop
[magcl,phasecl] = bode(numl, dclos, w) % clos-loop without prefilter
[magtot,phasetot] = bode(numtot, dentot, w) % clos-loop with prefilter
[magbf,phasebf] = bode(105, [1 10 45 105 105], w) % Bessel 4th order
% [magbf,phasebf] = bode(15, [1 6 15 15], w) % Bessel 3th, option
% plotting
subplot(2,1,1)
semilogx(w, 20*log10(magop),'w', w,(phaseop),'w--');
hold on
semilogx(w, 20*log10(magcl),'b', w,(phasecl),'b--');
semilogx(w, 20*log10(magtot),'g', w,(phasetot),'g--');
semilogx(w, 20*log10(magbf),'r', w,(phasebf),'r--');
wmin = min(w); wmax = max(w);
axis([wmin wmax -40 20])
xlabel('rad/sec')
ylabel('gain, dB; lag margin, degr')
grid
hold off
subplot(2,1,2)
sstep(numtot, dentot)
xlabel('time, sec')
ylabel('output')
grid
zoom on
Back
BOINTEGR
function [ao, p1, numr, denr] = bointegr(numc,denc,wb)
% BOINTEGR expresses the compensator transfer function as a sum of the
% term of the partial fraction expansion ao/(s + p1) which is dominant
% a lower frequencies and the rest of the terms combined into transfer
% function numr/denr, where the coefficients in the vectors of
% coefficients are in the descending powers of s. numc is up to
% 9th order. wb is the crossover frequency or anything close enough,
% it is used here only to scale the frequency axis; default is 1.
% Bode diagrams for the compensator and its two parallel paths are
% plotted.
% The default/demo is
% numc = [9.6082 29.262 48.975 16.016];
% denc = [1 6.8631 24.466 46.749 32.488 0]; wb = 1;
% function [ao, p1, numr, denr] = bointegr(numc,denc);
%
% Copyright (c) B. Lurie 6-19-98.
% See also BOLAGNYQ, BONYQAS, BOSTEP, BOCLOS, BOCOMP, NYQLOG, TFSHIFT from
% the Bodestep toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
% default
if nargin < 3
wb = 1;
end
if nargin < 2
numc = [9.6082 29.262 48.975 16.016];
denc = [1 6.8631 24.466 46.749 32.488 0];
end
format short e
% Partial fraction expansion
ld = length(denc);
[r,p,k] = residue(numc,denc);
ao = r(ld-1); p1 = p(ld-1);
% Combining the rest of the terms
if ld == 2
rr = 0;
pr = 1;
end
if ld == 3
rr = [r(ld-2)];
pr = [p(ld-2)];
end
if ld == 4
rr = [r(ld-3) r(ld-2)];
pr = [p(ld-3) p(ld-2)];
end
if ld == 5
rr = [r(ld-4) r(ld-3) r(ld-2)];
pr = [p(ld-4) p(ld-3) p(ld-2)];
end
if ld == 6
rr = [r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
pr = [p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
end
if ld == 7
rr = [r(ld-6) r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
pr = [p(ld-6) p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
end
if ld == 8
rr = [r(ld-7) r(ld-6) r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
pr = [p(ld-7) p(ld-6) p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
end
if ld == 9
rr = [r(ld-8) r(ld-7) r(ld-6) r(ld-5) r(ld-4) r(ld-3) r(ld-2)];
pr = [p(ld-8) p(ld-7) p(ld-6) p(ld-5) p(ld-4) p(ld-3) p(ld-2)];
end
[numr,denr] = residue(rr,pr,k);
numr = real(numr); denr = real(denr); % elimination of imaginary
% parts caused by numerical errors
% Print the answers
ao
p1
numr
denr
% if further breaking into parallel channels needed, rr and pr
wmin = floor(log10(wb)-3);
wmax = ceil(log10(wb)+2);
w = logspace(wmin,wmax);
[mc,pc] = bode(numc,denc,w);
[mr,pr] = bode(numr,denr,w);
[mi,pi] = bode(ao,[1 p1],w);
% plotting
wmin = min(w); wmax = max(w);
subplot(211)
semilogx(w, 20*log10(mc),'w', w, 20*log10(mr),'g', w, 20*log10(mi),'r');
axis([wmin wmax -40 80])
xlabel('rad/sec')
ylabel('gain, dB')
grid
subplot(212)
semilogx(w,(pc),'w', w,(pr),'g', w,(pi),'r');
axis([wmin wmax -200 20])
ylabel('phase shift, degr')
grid
hold off
zoom on
Back
BOCOMP
function [numc_pl,denc_pl] = bocomp(wb,numc,denc,kmot,resist,mom_inert);
% BOCOMP Compensator transfer function numc_pl(s)/denc_pl(s)for a dc
% motor system with loop response generated by bostep. The vectors
% NUMC_PL and DENC_PL contain polynomial coefficients in descending
% powers of s,
% The driver voltage gain is 1. The sum of its output resistance with
% the motor winding resistance is resist.The motor constant is kmot.
% The rigid body load moment of inertia is mom_inert.
% The loop crossover frequency in rad/sec is wb.
% The function numc(s)/denc(s) is the normalized compensator obtained
% with function bostep for the nominal plant (1/s)(-s+wc/bnc)/(s+wc/bnc).
% BOCOMP calculates: the compensator polynomial coefficients,
% the gain coefficient, the poles and the zeros,
% the biquad coefficients for complex poles and zeros,
% and plots: Bode diagrams for the plant (in green) and
% compensator (in read).
% Defaults/demo initiated by typing bocomp:
% mom_inert = 0.027; resist = 2; kmot = 0.7; wb = 100;
% denc = [ 1 6.8631 24.466 46.749 32.488 0];
% numc = [9.6082 29.262 48.975 16.016];
% [numc_pl,denc_pl] = bocomp(wb,numc,denc,kmot,resist,mom_inert)
%
% Copyright (c) B. J. Lurie 6-18-98.
% See also BOSTEP, BOCLOS, NYQLOG, TFSCHIFT, BOLAGNYQ, BONYQAS, BOINTEGR
% from the Bodestep toolbox.
% Further explanations are given in Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 6,
mom_inert = 0.027; ;
end
if nargin < 5,
resist = 2;
end
if nargin < 4,
kmot = 0.7;
end
if nargin < 3,
denc = [1 7.9547 33.536 74.518 58.275 0];
end
if nargin < 2,
numc = [13.046 45.855 76.234 22.827]; % ratio of 4th to 6th
end
if nargin == 0,
wb = 100;
end
% plant transfer function kmot/(resist*mom_inert)/[s(s+a)] where
a = kmot^2/(resist*mom_inert); %real pole of the plant
npm = kmot/(resist*mom_inert); dpm = [1 a 0]; % minimum phase plant
% without flex modes
% denormalized m.p. compensator for n.p. single integrator plant
[numc_wb, denc_wb] = tfshift(numc,denc,wb);
% denormalized m.p. compensator for motor plant
numc_pl = wb*deconv(conv(numc_wb,dpm),[1 0]);
denc_pl = conv(denc_wb,npm);
format short e
gain_coeff = numc_pl(1)/denc_pl(1);
roots_num = roots(numc_pl);
biquad_num = [1 -2*real(roots_num(2)) abs(roots_num(2))^2];
roots_den = roots(denc_pl);
biquad_den = [1 -2*real(roots_den(2)) abs(roots_den(2))^2];
% printout
gain_coeff
roots_num
biquad_num
roots_den
biquad_den
numc_pl
denc_pl
% plotting
wmin = floor(log10(wb)-3);
wmax = ceil(log10(wb)+1);
w = logspace(wmin,wmax);
[magc,phasc] = bode(numc_pl,denc_pl,w);
[magnpm, phasenpm] = bode(npm,dpm,w);
semilogx(w,20*log10(magc),'r', w,phasc,'r--'); % compensator
hold on
semilogx(w,20*log10(magnpm),'g', w,phasenpm,'g--') % plant
hold off
xlabel('rad/sec; phase plotted with dashed lines')
ylabel('dB, degr')
title('gain and phase of plant and compensator')
grid
zoom on
Back
NDCP
function NDCP(w1,n1,d1,n2,d2,n3,d3)
%NDCP Nyquist diagram on logarithmic plane for describing function of
% nonlinear dynamic compensator with parallel channels.
% -->df-->link1--
% --< >(+)-->link3--
% ------->link2--
%function NDCP(w1,n1,d1,n2,d2,n3,d3)
% The first channel incorporates a link with describing function
% df=0:0.2:1 and a linear link with numerator n1 and denominator d1
% of the transfer function. The second channel is a linear link with
% numerator n2 and denominator d2.
% The third linear link with transfer function numerator n3 and
% denominator d3 is connected in series to the parallel connection.
% The transfer functions are of the type n(s)/d(s), where n and d
% contain polynomial coefficients in descending powers of s.
% w1 is a frequency in rad/sec to be marked by both x and o.
% The default is: w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0]; n3=d3=1;
%
% Copyright (c) B. Lurie 6-30-99.
% For further explanations see Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = 1;
end
if nargin < 6
n3 = 1;
end
if nargin < 1
w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0];
end
denl = conv(conv(d1,d2),d3);
for df = 1 : -0.2 : 0,
n1 = df * n1;
v1 = conv(n1,d2);
v2 = conv(d1,n2);
ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
num1 = [adz v1] + v2;
else
adz = zeros(1,ll);
num1 = v1 + [adz v2];
end
numl = conv(num1,n3);
[mag, phase] = bode(numl, denl);
plot(phase, 20*log10(mag),'w',-180, 0, 'ro')
hold on
[a,b] = bode (numl,denl,w1);
plot(b,20*log10(a),'wx',b,20*log10(a),'wo')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode(numl,denl,w1*w); plot(b,20*log10(a),'w+');
end
end
hold off
grid
axis([-270 0 -20 70])
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90 -60 -30 0])
title('Nyquist diagram, x marks w = w1, + marks octaves')
xlabel('loop phase shift in degrees')
ylabel('loop gain in dB')
zoom on
Back
BNDCP
function BNDCP(w1,n1,d1,n2,d2,n3,d3)
%BNDCP Bode diagram for describing function of nonlinear dynamic
% compensator with parallel channels.
% ->df-->link1-
% --< >(+)-->link3--
% ->--->link2--
%function BNDCP(w1,n1,d1,n2,d2,n3,d3)
% The first channel incorporates a link with describing function
% df=0:0.2:1 and a linear link with numerator n1 and denominator d1
% of the transfer function. The second channel is a linear link with
% numerator n2 and denominator d2.
% The third linear link with transfer function numerator n3 and
% denominator d3 is connected in series to the parallel connection.
% The transfer functions are of the type n(s)/d(s), where n and d
% contain polynomial coefficients in descending powers of s.
% w1 is a frequency in rad/sec to be marked by both x and o.
%
% The default is: w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0]; n3=d3=1;
% Copyright (c) B. Lurie 6-30-99.
% For further explanations see Classical % Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = 1;
end
if nargin < 6
n3 = 1;
end
if nargin < 1
w1=10; n1=200; d1=[1 0 0]; n2=10; d2=[1 0];
end
denl = conv(conv(d1,d2),d3);
for df = 1 : -0.2 : 0,
n1 = df * n1;
v1 = conv(n1,d2);
v2 = conv(d1,n2);
ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
num1 = [adz v1] + v2;
else
adz = zeros(1,ll);
num1 = v1 + [adz v2];
end
numl = conv(num1,n3);
bode(numl, denl);
hold on
end
hold off
zoom on
Back
NDCB
function NDCB(w1,n1,d1,n2,d2,n3,d3)
%NDCB Nyquist diagram on logarithmic plane for describing function of
% nonlinear dynamic compensator with nonlinear feedback channel.
% -->(+)----->link1----->link3--
% (-) |
% ^--link2<--df<--
%function NDCB(w1,n1,d1,n2,d2,n3,d3)
% The first link is linear with numerator n1 and denominator d1
% of the transfer function. The feedback channel incorporates a
% link with describing function df=0:0.2:1 an a second linear
% link with numerator n2 and denominator d2.
% The third linear link with transfer function numerator n3 and
% denominator d3 is connected in series to the first link.
% The transfer functions are of the type n(s)/d(s), where n and d
% contain polynomial coefficients in descending powers of s.
% w1 is a frequency in rad/sec to be marked by both x and o.
%
% The default is: w1=10; n1=100; d1=[1 0]; n2=1; d2=[1];
% n3=[1 4]; d3=[1 15 0];
%
% Copyright (c) B. Lurie 6-30-99.
% For further explanations see Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = [1 14 0];
end
if nargin < 5
n3 = [1 2];
end
if nargin < 1
w1=8; n1=[100 100]; d1=[1 0 0]; n2=1; d2=[1];
end
numl = conv(conv(n1,d2),n3);
v1 = conv(d1,d2);
v2 = conv(n1,n2);
for df = 1 : -0.2 : 0,
v2 = df * v2;
ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
den1 = [adz v1] + v2;
else
adz = zeros(1,ll);
den1 = v1 + [adz v2];
end
denl = conv(den1,d3);
[mag, phase] = bode(numl, denl);
plot(phase, 20*log10(mag),'w',-180, 0, 'ro')
hold on
[a,b] = bode (numl,denl,w1);
plot(b,20*log10(a),'wx',b,20*log10(a),'wo')
for w = [0.0156 0.03125 0.0625 0.125 0.25 0.5 2 4],
[a,b] = bode(numl,denl,w1*w); plot(b,20*log10(a),'w+');
end
end
hold off
grid
axis([-270 0 -20 70])
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90 -60 -30 0])
title('Nyquist diagram, x marks w = w1, + marks octaves')
xlabel('loop phase shift in degrees')
ylabel('loop gain in dB')
zoom on
Back
BNDCB
function BNDCB(w1,n1,d1,n2,d2,n3,d3)
%BNDCB Bode diagrams for describing functions of nonlinear dynamic
% compensator with nonlinear feedback channel.
% -->(+)----->link1----->link3--
% (-) |
% ^--link2<--df<--
%function BNDCB(w1,n1,d1,n2,d2,n3,d3)
% The first link is linear with numerator n1 and denominator d1
% of the transfer function. The feedback channel incorporates a
% link with describing function df = 0:0.2:1 and a second linear
% link with numerator n2 and denominator d2.
% The third linear link with transfer function numerator n3 and
% denominator d3 is connected in series to the first link.
% The transfer functions are of the type n(s)/d(s), where n and d
% contain polynomial coefficients in descending powers of s.
% w1 is a frequency in rad/sec to be marked by both x and o.
%
% The default is: w1=10; n1=100; d1=[1 0]; n2=1; d2=[1];
% n3=[1 4]; d3=[1 15 0];
%
% Copyright (c) B. Lurie 6-30-99.
% For further explanations see Classical Feedback Control by
% B. J. Lurie and P. J. Enright, Marcel Dekker, 2000.
if nargin < 7
d3 = [1 14 0];
end
if nargin < 5
n3 = [1 2];
end
if nargin < 1
w1=8; n1=[100 100]; d1=[1 0 0]; n2=1; d2=[1];
end
numl = conv(conv(n1,d2),n3);
v1 = conv(d1,d2);
v2 = conv(n1,n2);
for df = 1 : -0.2 : 0,
v2 = df * v2;
ll = length(v1) - length(v2);
if ll < 0
adz = zeros(1,-ll);
den1 = [adz v1] + v2;
else
adz = zeros(1,ll);
den1 = v1 + [adz v2];
end
denl = conv(den1,d3);
bode(numl, denl);
hold on
end
hold off
zoom on
Back