Solution 9.10.2.11
The MATLAB statements
load bodeid11
topm = size(bodeid11)
top = topm(1,1)
w=bodeid11(1:top,1);;
mag = bodeid11(1:top,2);;
phase = bodeid11(1:top,3);;
semilogx(w,mag);;
grid on
axis([0.1,1000,-80,20])
print -deps 910211mag.eps
semilogx(w,phase);;
grid on
axis([0.1,1000,-180,-80])
print -deps 910211phase.eps
can be used to load the data and plot the magnitude response. Figure 1
shows the magnitude data with some asymptotes added. The transfer func-
tion clearly has a pole at the origin and a pole near s = ;80. Beyond that
the magnitude plot is not muchhelp. If welookatthe phase plot wesee
there maybeapolearound s = ;5andazerointhevicinityofs = ;10.
Atlowfrequency only the gain term and the term 1=s will be present. The
term 1=s byitself will cross the vertical line through ! =1rad/sec at 0 db.
The plot actually crosses at ;1 dB, so wecannd the gain to be
K =10
;1=20
=0:89125:
Then our rst guess at the transfer function is
G(s) =
0:891(1+ s=10)
s(1 + s=7)(1+s=80)
=
49:91(s+10)
s(s +7)(s+80)
:
Wecheckthe accuracy of the model bycomparing the actual magnitude and
phase to the phase of the derived transfer function, as shown in Figures 3
and 4 The the shape of the magnitude curveiscorrect but the gain is t0oo
low. The phase curvelooksallrightonthefrontandtheback, but is not
very close in the middle. Weconclude that a single zero will not pull the
1
10
-1
10
0
10
1
10
2
10
3
-80
-70
-60
-50
-40
-30
-20
-10
0
10
20
Figure 1: ABode magnitude plot
2
10
-1
10
0
10
1
10
2
10
3
-180
-170
-160
-150
-140
-130
-120
-110
-100
-90
-80
Figure 2: ABode phase plot
3
10
-1
10
0
10
1
10
2
10
3
-80
-70
-60
-50
-40
-30
-20
Figure 3: Comparison of actual and derived Bode magnitude plots
4
10
-2
10
-1
10
0
10
1
10
2
10
3
-180
-170
-160
-150
-140
-130
-120
-110
-100
-90
-80
Figure 4: Comparison of actual and derived Bode phase plots
5
10
-1
10
0
10
1
10
2
10
3
-80
-70
-60
-50
-40
-30
-20
-10
0
10
20
Figure 5: Comparison of actual and derived Bode magnitude plots
phase plot backtothe larger values werequre. If wehaveasecond zero
then we need another pole. So our next guess is
G(s) =
0:87(1+s=10)(1+s=12)
s(1 +s=5)(1+s=15)(1+ s=80)
=
43:55(s+10)(s+12)
s(s+5)(s+15)(s+80)
:
As shown in Figures 5 and 6, the t is muchbetterfor both the mag-
nitude and the phase. It is nowjustamatter of moving a pole and a zero
around to improvethe t. After some eort we nd Wecould stop here,
because wehaveaworkable transfer function. But consider
G(s) =
0:875(1+s=10)
2
s(1 +s=5)(1+s=20)(1+ s=80)
=
70(s+10)
2
s(s+5)(s+20)(s+80)
:
6
10
-2
10
-1
10
0
10
1
10
2
10
3
-180
-170
-160
-150
-140
-130
-120
-110
-100
-90
-80
Figure 6: Comparison of actual and derived Bode phase plots
7
10
-2
10
-1
10
0
10
1
10
2
10
3
-80
-70
-60
-50
-40
-30
-20
-10
0
10
20
Figure 7: Comparison of actual and derived Bode magnitude plots
The results are shown in Figures 7 and 8. As can be seen the t is perfect.
Thus, it maybehard to detect a near pole/zero cancellation, and wemay
just havetorootaround a bit until wegetagoodt. The MATLAB
program that does this analysis is
load bodeid11
topm = size(bodeid11)
top = topm(1,1)
w=bodeid11(1:top,1);;
mag = bodeid11(1:top,2);;
phase = bodeid11(1:top,3);;
semilogx(w,mag);;
grid on
axis([0.1,1000,-80,20])
print -deps 910211mag.eps
semilogx(w,phase);;
grid on
axis([0.1,1000,-180,-80])
8
10
-2
10
-1
10
0
10
1
10
2
10
3
-180
-170
-160
-150
-140
-130
-120
-110
-100
-90
-80
Figure 8: Comparison of actual and derived Bode phase plots
9
print -deps 910211phase.eps
z1 = 10
p1 = 0
p2 = 5
p3 = 80
Ktc = 10^(-1/20)
K=(Ktc * p2 * p3)/z1
wp = logspace(-2,3,20);;
jw = j*wp;;
mag1 = 20*log10( (K.*abs(jw+z1) )./ ( abs(jw +p1).* abs(jw + p2).*abs(jw+p3) ) );;
semilogx(w,mag,'k-',wp,mag1,'kd')
grid on
axis([0.1,1000,-80,-20])
print -deps 910211mag1.eps
phase1 = ( angle(jw + z1) -angle(jw +p1) -angle(jw + p2) - angle(jw+p3) )*180/pi;;
semilogx(w,phase,'k-',wp,phase1,'kd')
grid on
axis([0.01,1000,-180, -80])
print -deps 910211phase1.eps
z1 = 10
z2 = 12
p1 = 0
p2 = 5
p3 = 15
p4 = 80
Ktc = 10^(-1.2/20)
K=(Ktc * p2 * p3*p4)/(z1*z2)
mag2 = 20*log10( (K.*abs(jw + z1).*abs(jw+z2) )
./ ( abs(jw +p1).* abs(jw + p2).*abs(jw+p3).*abs(jw + p4) ) );;
semilogx(w,mag,'k-',wp,mag2,'kd')
grid on
axis([0.1,1000,-80,20])
print -deps 910211mag2.eps
phase2 = (angle(jw + z1) + angle(jw + z2)
-angle(jw +p1) -angle(jw + p2) - angle(jw+p3)-angle(jw + p4))*180/pi;;
semilogx(w,phase,'k-',wp,phase2,'kd')
grid on
axis([0.01,1000,-180,-80])
print -deps 910211phase2.eps
z1 = 10
z2 = 10
p1 = 0
p2 = 5
p3 = 20
p4 = 80
10
Ktc = 0.875
K=(Ktc *p2* p3 * p4)/(z1*z2)
mag3 = 20*log10( ( K*abs(jw + z1).*abs(jw + z2) )
./ ( abs(jw +p1).* abs(jw + p2).*abs(jw+p3).*abs(jw+p4) ) );;
semilogx(w,mag,'k-',wp,mag3,'kd')
grid on
axis([0.01,1000,-80,20])
print -deps 910211mag3.eps
phase3 = (angle(jw + z1) + angle(jw + z2)
-angle(jw +p1) -angle(jw + p2) - angle(jw+p3)-angle(jw + p4) )*180/pi;;
semilogx(w,phase,'k-',wp,phase3,'kd')
grid on
axis([0.01,1000,-180,-80])
print -deps 910211phase3.eps
11