| Dr. Boris J. Lurie |
|
|
|
These are the listings of all MATLAB and SPICE scripts found in Classical Feedback Control by B. Lurie and P. Enright, exclusive of the functions in Appendix 14. The scripts can be copied freely for nonprofit use. The author bears no responsibility for the results of the scripts' applications.
| Chapter 1 | Chapter 2 | Chapter 3 | Chapter 4 | Chapter 5 | Chapter 6 |
| Chapter 7 | Chapter 10 | Chapter 11 | Appendix 2 | Appendix 5 | Appendix 13 |
Example 3. The frequency responses can be plotted with MATLAB ® script:
w = logspace(-1, 2); % log scale of angular % frequency w num = 5000; den = [1 55 250 0]; bode(num, den, w) a = [1 0]; b = [1 5]; c = [1 50]; ab = conv(a,b); den = conv(ab,c)
Example 1. The L-plane Nyquist diagram is charted with MATLAB script
num = [20 10]; den = [1 10 20 1 0];
[mag, phase] = bode(num, den);
plot(phase, 20*log10(mag), 'r', -180, 0, 'wo')
title('L-plane Nyquist diagram')
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
The frequency response can be plotted using MATLAB with:
w = logspace(-2,1); % frequency range % 0.01 to 10 rad/sec num = 1; den = [1 2.1 0.2 0 0]; bode(num, den, w)
The plots can be made in MATLAB with:
w = logspace(-1,1); % freq range 0.1 to 10 rad/sec den = [1 15 50 0 0 0]; num = [0 0 0 50 27.5 0.25]; % equal length of the vectors g = num + den; % makes the addition allowable bode(num, den, w) % for T hold on bode(g, den, w) % for F bode(num, g, w) % for M hold off
Example 1. The transient response is found with:
num = [0 0 0 50 27.5 0.25]; den = [1 15 50 0 0 0]; g = num + den; step(num, g) grid
22 (a) Exact solution found with MATLAB script:
A = [2 0.2 0.3; 0.1 2.1 0.1; 0.04 0.1 1.9]; inv(A)
ans = 0.5038 -0.0443 -0.0772
-0.0235 0.4795 -0.0215
-0.0094 -0.0243 0.5291
Example 1. The open-loop and closed-loop poles can be calculated with the MATLAB commands
n = [0 0 10 20]; d = [1 1 1 0]; roots(d) %open-loop poles ans = 0 -0.5000 ± 0.8660i roots(n + d) %closed-loop poles ans = -1.6551 0.3275 ± 3.4608i
The root loci in Fig. 3.11 are plotted with:
n = [10 20]; d = [1 1 1 0];
rlocus(n,d)
hold on
k = [0.05 0.1 0.2 0.5 1];
rlocus(n, d, k)
title('k = [0.05 0.1 0.2 0.5 1]');
hold off
The weight function is charted with the script
u = linspace(-3,3,200); b = log(coth(abs(u/2))); plot(u,b,'w'); grid
14 The following is a SPICE input file:
** unst_pl.cir for feedback loop with unstable plant
* compensator, inverting, zero s=10, pole s=40,
* crossover s=20
GC 2 0 1 0 1
RC1 1 0 1MEG ; for SPICE, to make node (1) non-floating
RC2 2 0 1
RC3 2 3 0.33333
LC 3 0 0.033333 ; Z2=(s+10)/(s+4)
* plant loop, non inverting, closed loop 0.1s/(s^2+0.2s+4)
EP 5 0 2 4 10
GP1 4 0 5 0 1
RP 4 0 0.5
CP1 4 0 10
LP 4 0 0.025 ; Z5=0.1s/(s^2+0.05s+4)
* plant integrators, (4) to (6) to (7)
GP2 6 0 0 5 80
GP3 7 0 0 6 1
RP2 5 0 1MEG ; for SPICE
RP3 6 0 1MEG ; for SPICE
RP4 7 0 1MEG ; for SPICE
CP2 6 0 1
CP3 7 0 1
* loop closing resistor
RL 7 1 1MEG ; to close the loop, place semicolon
; after M to reduce RL
*
VIN 8 0 AC 1
RIN 8 1 1
.AC DEC 100 0.1 20
.PROBE ; (if using PSPICE)
.END
Example 1. The loop transfer function is plotted with
n = conv([11 55 110 36], [-1 10]); d = conv([1 7.7 34 97 83 0 0], [1 10]); w = logspace(-1,1,200); bode(n,d,w)
The L-plane Nyquist diagram is plotted with
[mag, phase] = bode(n,d,w); plot(phase,20*log10(mag),'w', -180,0,'wo'); grid
The closed-loop transfer function M = n/g where g = n + d is plotted with
n = conv([0 0 0 1],n); g = n + d; w = logspace(-1,1,200); bode(n,g,w)
The prefilter response is plotted with
nr = [1 1 0.81]; dr = [1 2 0.81]; bode(nr, dr, w)
and the closed-loop response with the prefilter is plotted with
nc = conv(nr, n); gc = conv(dr, g); bode(nc, gc, w)
By compensating for the hump with a small, 1.6 dB additional notch with
nr2 = [1 0.5 1]; dr2 = [1 0.6 1]; nc2 = conv(nr2, nc); gc2 = conv(dr2, gc); bode(nc2, gc2, w) step(nc2, gc2)
The effects of variations in k can be calculated as follows:
n = conv([11 55 110 36], [-1,10]); d = conv([1 7.7 34 97 83 0], [1 10]); k = 1; n = conv(n, k); % specify k w = logspace(-1,1,200); bode(n, d, w) % open-loop response [mag, phase] = bode(n, d, w); plot(phase, 20*log10(mag)) % Nyquist diagram grid n = conv([0 0 1], n); g = n + d; bode(n, g, w) % closed-loop response nr = [1 1 0.81]; dr = [1 2 0.81]; bode(nr, dr, w) % prefilter notch response nc = conv(nr, n); gc = conv(dr, g); bode(nc, gc, w) % closed-loop response with notch nr2 = [1 0.5 1]; dr2 = [1 0.6 1]; bode(nr2 ,dr2, w) % second notch response nc2 = conv(nr2, nc); gc2 = conv(dr2, gc); bode(nc2, gc2, w) % closed-loop response with notches step(nc2, gc2) % step response with notches
Example 2. The MATLAB code for this transfer function's numerator and denominator is
k = 0.7; res = 2; ja = 0.027; a = k*k/(res*ja); nc1 = conv(res*ja/k,[1 a]); dc1 = [1 0];
Example 3. The numerator ns and the
denominator ds can be found with
n = conv([11 55 110 36], [-1,10]); wb = 75.4; d = conv([1 7.7 34 97 83 0 0], [1 10]); format short e; [ns, ds] = lp2lp(n, d, wb)
and the numerators and denominators of the prefilter with
nr = [1 1 0.81]; dr = [1 2 0.81]; [nrs, drs] = lp2lp(nr, dr, wb) nr2 = [1 0.5 1]; dr2 = [1 0.6 1]; [nr2s, dr2s] = lp2lp(nr2, dr2, wb)
Example 2. The gain and phase responses are plotted in Fig. 4.36 with
n = conv([1 0 10],[1 0 50]); d = conv([1 0],[1 0 11]); d = conv(d,[1 0 54]); w = logspace(0, 1, 1000); bode(n,d,w)
4
fst = 120; Q = 10; fb = fst/2^[20*log10(Q)/18 + 2] fb = 13.8881
Example 1. The polynomial coefficients can be found from the roots with
rn = [-0.5 -2 -20]; rd = [-0.2 -0.3 -10]; num = poly(rn); denum = poly(rd);
and the Bode plot can be obtained with
bode(num, denum)
or by using the command freqs.
After multiplication of the polynomials in the numerator and denominator:
n = 4 * 2.8 * conv([1 0.42],[1 2 4]) n = 11.2000 27.1040 54.2080 18.8160 d = conv(conv([1 0.06],[1 1.4]), conv([1 2.8],[1 2.4 9 0])) d = 1.0000 6.6600 23.3960 48.5880 38.1125 2.1168 0
The Nyquist diagram on the L-plane is plotted with
w = logspace(-1, 1); [mag, phase] = bode(n, d, w);
plot(phase,20*log10(mag),'r', -180, 0,'wo')
title('L-plane Nyquist diagram')
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
The SPICE input file is shown below.
*** PID example Figs. 5.19, 5.20 ***
ES 3 0 1 2 1 ; input signal summer
RSR1 1 0 1MEG ; leakage resistor
RSP2 2 0 1MEG ; leakage resistor
***
GSAT 5 0 0 3 0.001 ; saturation in I-path
RSAT 5 0 1K ; threshold = (0.7+VT1)*GI1/GSAT
D1 5 6 DIODE
D2 7 5 DIODE
.MODEL DIODE D
VT1 6 0 1V
VT2 0 7 1V
***
GI1 8 0 0 5 1 ; I-path
CI2 9 0 0 8 10 ; integral coefficient
RSP8 8 0 1MEG ; leakage resistor
***
GP 9 0 0 3 2 ; proportional coefficient
***
GD1 4 0 0 3 3 ; differential coefficient
LD 4 0 1
GD2 9 0 0 4 1
***
RS 9 0 1 ; summing resistor
***
GA1 10 0 9 0 1 ; actuator gain coefficient
***
RSATA 10 0 1K ; saturation in actuator
D3 10 11 DIODE
D4 12 10 DIODE
VT3 11 0 1V
VT4 0 12 1V
GA2 13 0 0 10 1 ; force source
***
CP1 13 0 5 ; mass of the main body
RSP13 13 0 1MEG ; leakage resistor
LP2 13 14 0.1 ; spring of flexible mode
CP2 12 0 0.5 ; mass of second body
RP 14 15 0.02 ; losses in the flexible mode
GINT 2 0 0 13 1
CINT 2 0 1 ; integrator to generate position
*** V10 is force, V13 is velocity, V2 is position
***
VTEST 1 0 AC 1 ; use only when frequency responses
; are tested
.AC DEC 20 0.01 10 ; use only when frequency responses tested
** Pulse (Vmin Vmax delay rise fall width period)
* VPULSE 10 PULSE ( 0V 10V 0S 0S 0S 500 500 )
; when transient responses tested
* .TRAN 0.1 10 ; when transient responses tested
.PROBE ; or other graphical postprocessor
.END
Example 2. The same problem is solved with MATLAB by
n = [1 5]; d = [1 10]; fs = 100; % sampling frequency in Hz [nd,dd] = bilinear(n,d,fs) nd = 0.9762 -0.9286 dd = 1.0000 -0.9048
15 MATLAB function conv is used to multiply
the polynomials in the numerator:
a = [601]; b = [1 0.88]; c = [1 3.5 19.4]; ab = conv(a,b); num = conv(ab,c)
and in the denominator:
d = [1 0]; e = [1 0.22]; f = [1 5.5]; g = [1 10.7 157]; de = conv(d,e); def = conv(de,f); deff = conv(def,f); den = conv(deff,g)
16 (1) This response plotted with MATLAB commands
n = [10.8 37.9 52.8]; d = [1 6.4 25.6 64 0 0]; w = logspace(-1,1); bode(n,d,w)
With:
n1 = conv(3.3,[ 1 0.3]); % during iterations, adjust the zero % (try also zero 0.75, pole 1.5) d1 = conv([1 4 0], [1 1]); % during iterations, adjust the pole d2 = [1 2.4 16]; d1d2 = conv(d1,d2); d = conv(d1d2,[1 0]) n = conv(n1,d2) + conv(7.5,d1) % vectors have equal length since % both polynomials are cubic w = logspace(-1,1); bode(n,d,w)
(2) With MATLAB, the function can be decomposed:
num = [601 2632 13510 10260];
den = [1 21.9 309.7 2117.8 5200.4 1044.8];
[Res,Pol,K] = residue(num, den)
Res = Pol = K = []
1.0e+002 *
0.1598 - 0.2042i -5.3460 +11.3430i
0.1598 + 0.2042i -5.3460 -11.3430i
-0.1684 - 4.4410i -5.4940 + 0.1396i
-0.1684 + 4.4410i -5.4940 - 0.1396i
0.0172 - 0.0000i -0.2200
The polynomials can be found as follows:
num_prod = [2*a (-2*(a*c + b*d))] den_prod = [1 (-2*c) (c*c + d*d)]
When d is the result of calculation inaccuracy, then, approximately,
num_prod = [2*a (-2*a*c)] den_prod = [1 (-2*c) (c*c)]
(3) With:
num = conv([1 0.88],[ 1 3.5 19.4]) d1 = conv([1 0.22], [1 5.5]); den = conv(d1, [1 10.7 157])
With MATLAB, the function can be decomposed as follows:
num = [1 4.38 22.48 17.072];
den = [1 16.42 219.414 910.987 189.97];
[Res,Pol,K] = residue(num, den)
Res = Pol = K = []
0.3889 + 0.2973i -5.3500 +11.3304i
0.3889 - 0.2973i -5.3500 -11.3304i
0.2072 -5.5000
0.0151 -0.2200
The sum of the two complex pole fractions can be found as follows:
a = 0.3889; b = 0.2973; c = -5.35; d = 11.3304; prod_num = [2*a (-2*(a*c + b*d))] prod_den = [1 (-2*c) (c*c + d*d)]
In MATLAB, the response can be found using
r = [r1 r2 .. ri ..]; p = [p1 p2 .. pi ..]; [num, den] = residue(r,p); bode(num,den)
Example 2. The MATLAB code calculates the residues:
num = [0.000001 0.0005]; den = [1 2100 200000]; [r, p] = residue(num, den)
Example 1.
[A,B,C,D]=linmod('file_name'); bode(A,B,C,D,1)
Example 1. The logarithmic Nyquist diagram for
TP is plotted in Fig. 10.24(b) with the
nyqlog function from the Bodestep toolbox described in Appendix 14:
np = [2 1]; dp = [1 2 0 0]; nyqlog(1,np,dp)
(A) The L-plane Nyquist diagrams are plotted with the following script:
np = [2 1]; dp = [1 2 0 0]; nyqlog(1,np,dp); hold on
ne = 2; de = [1 2 0]; nyqlog(1,ne,de); hold on
ng = 1; dg = [1 2 2 0]; nyqlog(0.5,ng,dg); zoom off
gtext('TP'); gtext('TE'); gtext('G') %place labels with mouse
Example 2.
(C)
npg = [0 0 2 1.2 0.1]; dpg = [1 2 1.6 0.16 0]; gpg = npg + dpg; sstep(npg,gpg); grid ne = [0.4 1.04 0.1]; de = [1 2 1.6 0.16 0]; ge = n + d; sstep(n,g); grid
Example 1. The SPICE simulation input file is shown below.
**** ch9ex1.cir for iso-w simulation of NDC Figs. 11.21(b), **** 11.22 *** input integrator G2 2 0 0 1 1 C2 2 0 1 R2 2 0 1MEG *** feedback summer G3 3 0 7 2 1 R3 3 0 1 *** kDF path: G5 =.1, 1, or 10 G5 5 0 0 3 10 R5 5 0 1 *** forward path, inverting G4 4 0 3 0 1 C4 4 0 1 R4 4 0 1MEG *** forward summer, the output is VDB(6) G6 6 0 4 5 1 R6 6 0 1 *** feedback path G7 7 0 0 5 1 C7 7 0 1 R7 7 0 1MEG *** VIN 1 0 AC 1 RIN 1 0 1MEG .AC DEC 20 .001 10 .PROBE .END
1(a) When employing SPICE, the input file is
*ch9oc1.cir for determining oscillation condition *T = 200(s+300)(s+600)/[s(s+20)(s+50)(s+100)], open loop G2 2 0 0 1 200 ; gain, zero at 300 L2 2 8 1 G3 3 0 0 2 1 ; zero at 600 L3 3 9 1 R3 9 0 600 G4 4 0 0 3 1 ; pole at 20 R4 4 0 0.05 C4 4 0 1 G5 5 0 0 4 1 ; pole at 50 R5 5 0 0.02 C5 5 0 1 G6 6 0 0 5 1 ; pole at 100 R6 6 0 0.01 C6 6 0 1 G7 7 0 0 6 1 ; integrator R7 7 0 1MEG ; de-floating resistor C7 7 0 1 VIN 1 0 AC 1 RIN 1 0 1MEG ; source loading resistor .AC DEC 20 1 100 .PROBE .END
Example 2. The time-function can be plotted either with
num = [5 7]; den = [1 3 2 0]; impulse(num, den)
or with
num = [5 7]; den = [1 3 2]; step(num, den)
The programs' listing follows:
Function find_phase2
function [phase] = find_phase2(magdb, natfreq)
% function [phase] = find_phase2(magdb ,freq)
% This routine uses the Bode Integral to generate phase data.
% from a magnitude response of a m.p. function.
% magdb: row gain vector given in dB
% freq: row frequency vector given in rad/sec
% The Magnitude and Frequency vectors must be the same length.
% Before running this function prepare table with table_maker.
[row,col] = size(magdb); if col == 1, magdb = magdb'; end;
[row,col] = size(natfreq); if col == 1, natfreq = natfreq'; end;
%%% table_maker % calls the function creating the table
load table % load data needed for toe and tail calculations:
% table con hilimit lolimit numentries numintstep
points = length(natfreq);
numsteps = points - 1;
% The following variables are for the lookuptable (u domain)
ilnfreq = log(natfreq(1));
flnfreq = log(natfreq(points));
toeslope = (magdb(2) - magdb(1))/(log(natfreq(2))-ilnfreq);
tailslope = (magdb(points) - magdb(numsteps))/(flnfreq -log(natfreq(numsteps)));
dmagdb = magdb(2:points) - magdb(1:numsteps);
dnatfreq = natfreq(2:points) - natfreq(1:numsteps);
deriv = dmagdb./dnatfreq;
nnfreq = natfreq(1:numsteps);
w1 = nnfreq';w2 = natfreq(2:points)';
for I = 1:points,
% The next lines perform the integration...
wc = natfreq(i)*ones(w1);
weights = (w2+ wc).*log(w2 + wc) - (w1 + wc).*log(w1 + wc) +...
(wc-w2).*log(abs(wc-w2+eps)) - (wc-w1).*log(abs(wc-w1+eps));
looplnfreq = log(natfreq(i));
u = [flnfreq ilnfreq] - [looplnfreq]*[1 1] + [eps eps];
ind = (log10(abs(u))+[-log10(lolimit)-log10(lolimit)])/con+[1 1];
ind = max([1 1;ind]);ind=min([numentries numentries;ind]);
tailtoe = abs(pi^2/4*[1;1]-table(ind',2));
phase(i) = (deriv*weights + [tailslope toeslope]*tailtoe)/pi;
end
Function table_maker.m
function table_maker % table_maker.m, called by function find_phase2. This routine % creates the lookup table for calculation the toe and tail % 1e-15 is used for zero, 650 is used for infinity % reasonable limits: 1e-9 to 100 % This contains nearly all of the area (to 6 places). lolimit = 1e-9; hilimit = 100; numentries = 4001; % number of table entries, should be odd numintstep = 100; % number of steps for each integration vector = logspace(log10(lolimit),log10(hilimit),numentries); table(:,1) = vector'; [table(1,2)] = integral_u(1e-15,lolimit, numintstep); for k = 2:length(vector), [table(k,2)]=table(k-1,2)+integral_u(table(k-1,1),table(k,1), numintstep); end clear vector k con = log10(table(7,1))-log10(table(6,1)); save table con hilimit lolimit numentries numintstep table clear con hilimit lolimit numentries numintstep table
Function integral_u
function [trap]=integral_u(u1,u2,numsteps)
% function [trap]=integral_u(u1,u2,numsteps)
% This routeens integrates in u domain.
% It may not handle some special cases, however.
% vector=linspace calculation(u1,u2,numsteps+1);
if (u1==0),
u1=1e-15;
end
if (u2==0),
u2=1e-15;
end
vector = logspace(log10(u1),log10(u2),numsteps+1);
% This next calculation is the logu_function
valvector = log(abs((exp(vector)+1)./(exp(vector)-1)));
delta = vector(2:numsteps+1)-vector(1:numsteps);
%up = valvector(1:numsteps)*delta';
%down = valvector(2:numsteps+1)*delta';
trap=(.5*(valvector(1:numsteps)+valvector(2:numsteps+1)))*delta';
if trap==NaN,
trap=0;
disp('Warning: NaN found in integration, result set to 0')
end
The Bode diagram is plotted with
w = logspace(1,3);
[A,B,C,D] = linmod('mlsplan');
bode(A,B,C,D,1,w)
hold on
% remove the connection to the input of block Z2
bode(A,B,C,D,1,w)
The Bode diagrams are obtained with the following script:
npzt = [2.33e-10 6.7e-6 -5.14e11 1.26e17 1.5e21 9.5e24 1.2e28]; dpzt = [1 2.83e5 7e9 9.18e13 6.5e17 1.684e21 0 0]; nvc = [2.5e9]; dvc = [1 40 400 0]; ntot1 = conv(npzt,dvc); ntot2 = conv(dpzt,nvc); adz = zeros(1,length(ntot1)-length(ntot2)); ntot = ntot1 + [adz ntot2]; dtot = conv(dvc,dpzt); w = logspace(1,4.5); bode(npzt,dpzt,w); hold on bode(nvc,dvc,w); bode(ntot,dtot,w); hold off zoom on
The logarithmic plane Nyquist diagrams are plotted with
nyqlog:
nyqlog(3700, nvc,dvc) hold on nyqlog(3700, npzt,dpzt) hold on nyqlog(3700, ntot,dtot) hold off
The Bode diagram is obtained with the script
npzt = [2.33e-10 6.7e-6 -5.14e11 1.26e17 1.5e21 9.5e24 1.2e28]; dpzt = [1 2.83e5 7e9 9.18e13 6.5e17 1.684e21 0 0]; nvc = [2.5e9]; dvc = [1 40 400 0]; n1 = conv(npzt,dvc); d1 = conv(dvc,dpzt); d11 = [0 n1] + d1; d2 = conv(dpzt,nvc); deq = d11 + [0 0 0 d2]; w = logspace(0,4); % nyqlog(1260,-n1,deq) bbode(n1,deq,w)
The Bode diagram of this transfer function found with the script
npzt = [2.33e-10 6.7e-6 -5.14e11 1.26e17 1.5e21 9.5e24 1.2e28]; dpzt = [1 2.83e5 7e9 9.18e13 6.5e17 1.684e21 0 0]; nvc = [2.5e9]; dvc = [1 40 400 0]; n1 = conv(nvc,dpzt); d1 = conv(dvc,dpzt); d2 = conv(npzt,dvc); deq = d1 + [0 d2]; w = logspace(1,4); bode(n1,deq,w)