Licence
This worksheet and the package of procedures entitled "Bode2", may be used, adapted, or copied, provided that copyright is acknowledged and the whole of this Licence section is included in any copies, or adaptations, of the worksheet, or of the package.
This worksheet and the "Bode2" package may be used for educational purposes free of charge.
Commercial users are not required to pay a licence fee for this worksheet, or for the "Bode2" package. Instead, they are requested to make a donation to Cancer Research UK, National Office, 61 Lincoln's Inn Fields, LONDON, WC2A 3PX, UK. The donations line telephone number is +44 (0)20 7009 8820; the Fax number is +44 (0)20 7269 3100; and the Email address is "customerservices@cancer.org.uk". The size of donation is left to the user's discretion.
This worksheet and associated "Bode2" package of procedures are offered without warranty of any kind whatsoever.
Introduction
The servo described in this worksheet formed part of an automatic aircrafttracking radar in World War II that was used by C.H. Dowker in his work on minimizing the effect of noise. The servo's input was a signal proportional to the angular error between the centre of the radar beam and the aircraft. Naturally, the signal contained a great deal of noise generated by changes in the radar signal caused by random changes in radar propagation and reflection from the aircraft's surfaces. By appropriate choice of the servo's constants, Dowker found that it was possible to minimize the effect of noise on the servo while obtaining good performance from it. Chapter 6 of the reference gives the method used by Dowker for obtaining the best choice of servo constants.
This worksheet obtains the loop transfer function of the proportional control servo with tacho feedback compensation, as shown in the block diagram. The servo's constants, derived by C.H. Dowker at equations 87, 89 and 90 in Chapter 8 of the reference, are then inserted in the loop transfer function. These are the constants needed to minimize the effect of noise while producing the best response from the servo. The loop transfer function is then plotted in Nyquist diagrams and in Bode attenuation and phase margin diagrams, which yield parameters describing the servo's closed loop stability and performance. Plots showing the servo's frequency response, peak magnification, damped natural frequency and bandwidth are produced, together with transient responses to unit step functions of input position and of velocity. Important operational parameters are obtained from these plots.
The nomenclature used in the reference has been preserved in this worksheet as far as possible, so that the reader can move easily from one to the other.
The block diagram shown at Figure 8.4 of the reference, in slightly different form, has been inserted into this worksheet as an OLE2 object (or as an image in the Standard interface).
Loop Transfer Function
Initialisation
The "Bode2" package will be required for plotting the Bode attenuation and phase margin diagrams. It is also necessary to call the "plots" package and to change the imaginary unit to j, instead of I. It is also convenient to define some settings that will be used in the various plots.
> 
with(plots): with(inttrans): with(Bode2): interface(imaginaryunit=j): 
Warning, the name changecoords has been redefined
Loop Transfer Function
By inspection, the inner loop transfer function (from amplifier output to the servo's output) is
where Yf(s) is the forward loop transfer function =
, and
Yfb(s) is the feedback loop transfer function =
(also referred to as H).
By definition, the loop transfer function for the servo is,
=
It is convenient to abbreviate
to Y, and similarly for the other loop transfer functions.
> 
Yf := Km/(s*(Tm*s + 1)); 
> 
Yfb := (Kt*s)*(Kf*Ta*s/(Ta*s + 1)); 
From the block diagram, the gain of the tacho loop is seen to be
. Substituting this in Y, and collecting terms in s,
> 
Y := collect(algsubs(Kt*Kf=A,Y),s); 
This is the final symbolic form of the loop transfer function, before substituting values for the gains and time constants.
Steady State Error
At this point, it is convenient to determine the steady state error in response to a constant velocity input,
. From the block diagram,
> 
epsilon(s) := 'theta[i](s)/(1+Y)'; 
The Laplace transform of
is
. Hence,
> 
epsilon(s) := eval(subs(theta[i](s)= Omega/s^2,epsilon(s))); 
The steady state error is obtained by using the final value theorem of the Laplace transformation,
ss error =
=
.
> 
`ss error` := limit(s*epsilon(s),s=0); 
System Constants & Loop Transfer Function
The system constants were chosen as follows (see equations 87, 89 and 90 in Chapter 8 of the reference)
A = Tacho loop gain = Kt Kf = 14,
where Kt = Tacho Gain and Kf = Filter Gain (Attenuation),
Ta = Tacho filter time constant = 0.252 second,
Tm = Motor (including load) time constant = Jm/fm = 0.316 second,
where Jm = motor and load inertias and fm = viscous friction and back e.m.f.
Ga = Amplifier Gain Error Measuring Device Gain = 158 volts/radian,
Km = Motor Gain = radian/sec/volt,
Kv = Ga x Km = Velocity error coefficient for the servo = 158/second.
(Note. In this case the motor gain Km is unity, so Kv is numerically equal to Ga, the amplifier gain.)
Assigning these values,
> 
A := 14: Ga := 158: Km := 1: Ta := 0.252: Tm := 0.316: Kv := Ga*Km: 
The loop transfer function in the frequency domain, Y(f), is
> 
`Y(f)` := subs(s=j*omega,Y); 
Nyquist Diagram
Nyquist Diagram
The loop transfer function Y(f) is plotted on the Nyquist diagram as follows,
> 
p1 := complexplot(`Y(f)`,omega=3..1000,color=BLUE): 
> 
p2 := plot([[1,0]],style=POINT,symbol=circle,symbolsize=14,color=BLACK): 
> 
p3 := coordplot(polar,[0..1,Pi/2..3*Pi/2],grid=[5,13],color=[RED,RED],linestyle=[1,2]): 
> 
display([p1,p2,p3],title="Nyquist Diagram");

The Nyquist point,
, is shown on the diagram surrounded by a small black circle, for ease of identification. The loop transfer function shows a phase margin of slightly less than 60 degrees at feedback cutoff (unity magnitude), which indicates a very stable  possibly overdamped  system.
Closed Loop Frequency Magnification
The peak magnification, M, of the closed loop frequency response is obtained by plotting a circle of radius
, with its centre at
on the real axis of the Nyquist diagram, so that it just touches the loop transfer function. The loop transfer function's phase margin of nearly 60 degrees suggests that the peak magnification of the closed loop frequency will be in the region of 1.2  1.3, so four M circles are plotted on the Nyquist diagram thus
> 
M := [1.2,1.25,1.3,1.35]: 
> 
x := seq((M[i]^2/(M[i]^21) + M[i]/(M[i]^21)*cos(theta)),i=1..4): 
> 
y := seq((j*M[i]/(M[i]^21)*sin(theta)),i=1..4): 
> 
p := seq(complexplot(x[i]+y[i],theta=0..2*Pi,color=BLACK),i=1..4): 
> 
display(p1,p2,p, scaling=constrained, title="Nyquist Diagram with Magnification Circles\nM = 1.2, 1.25, 1.3, 1.35" ); 
As
, the radius of the M circle increases and its centre approaches
, with the circle becoming a vertical line through the Nyquist point. Noting that the loop transfer function runs between the second and third largest circles, it will touch a circle with M between 1.25 and 1.3. Choosing M = 1.27, and plotting again,
> 
M := 1.27:x := M^2/(M^21) + M/(M^21)*cos(theta): 
> 
y := j*M/(M^21)*sin(theta): 
> 
p := complexplot(x+y,theta=0..2*Pi,color=BLACK): 
> 
display( [p1,p2,p], scaling=constrained, title="Nyquist Diagram with Magnification Circle M = 1.27"); 
The peak magnification,
, of the closed loop frequency response is seen to be 1.27. This agrees with equation 109 in Section 8.8 of the reference.
> 
unassign('p1','p2','p3','p','x','y','M'); 
Bode Diagrams
Attenuation & Phase Margin Diagrams
The Bode diagrams, i.e. the loop gain of the transfer function, Y(f), together with the asymptotic straightline approximations, and the phase margin plotted against normalised frequency,
, are produced by procedures in the "Bode2" package. These procedure accept the loop transfer function expressed as a function of the Laplace operator, s, or as a function of the complex frequency. The latter is used in this worksheet for ease of comparison with the transfer functions as shown in the reference.
To obtain the straightline asymptotic approximations, the denominator of the loop transfer function must be factored, thus
> 
`factored Y(f)` := numer(`Y(f)`)/factor(denom(`Y(f)`),complex); 
The procedure for the straightline approximation assumes that the loop transfer function is expressed in canonical form with gain K in the numerator,
in the denominator, and only terms of the form
in the numerator and denominator. In this case, N = 1 and the only time constant in the numerator is 0.252. To reduce the denominator to the canonical form, it must be divided by the product 0.2453105557 51.19129784, yielding two time constants
and
. The numerator must also be divided by this product, which results in K = 158.
This is the value of Kv shown in equations 87, 89 and 90 of the reference. It is worth noting that for a Type 1 servo, provided its loop transfer function is expressed in canonical form, K is the velocity error coefficient. It is normally shown as Kv and has the units of 1/sec.
The straightline approximations (shown in red) can be superimposed on the loop gain (shown in blue) by using the "display" command.
> 
display(attensemilog(`Y(f)`,.01,100),
linesemilog(158,1,[.252],[1/.2453105557,1/51.19129784],.01,100)); 
The gain Kv = 158 is shown by the intercept of the initial 6dB/octave slope of the loop gain with the
axis at 42dB = 158. Feedback cutoff (also known as gain crossover) occurs about 1/3 the way down a 6dB/octave slope, at about
= 10. (Or, since the normalizing time constant, T, is 1 second,
= 10/sec.)
The phase margin is plotted thus,
The big increase in phase margin between
= 4/second and
= 44/second is caused by the tacho feedback loop producing a long 6dB/octave attenuation between these frequencies.
Improved Phase Margin Diagram
The shapes of the Bode diagrams above with semilog scales, particularly that for phase margin, are rather distorted below about
= 0.1. (This is a known problem with Maple's semilog plots  see Robert Israel's Maple Advisor.) The diagrams can be improved by using the "attenlog" and "phaselog" procedures in the Bode2 package, which plot the loop gain and phase margin against
. For example,
> 
phaselog(`Y(f)`,.01,100); 
Closed Loop Response
Closed Loop Transfer Function
The closed loop transfer function, or system response, is given by
=
Again abbreviating
to Y and similarly for
,
> 
`C/R` := factor(`C/R`); 
The closed loop transfer function in the frequency domain, C/R(f), is
> 
`C/R(f)` := subs(s=j*omega,`C/R`); 
Closed Loop Frequency Response
The closed loop frequency response is usually plotted as magnitude in dB against
thus,
> 
Freq_resp := 20*log10(simplify(abs(`C/R(f)`))): 
> 
p := plot(Freq_resp,w=1..1.5,10..3,labels=["log10 (omega T )",""],title="Closed Loop Servo Frequency Response \n Magnitude dB"): 
> 
q := pointplot([[1,3],[1.5,3]],connect=true,colour=black,linestyle=2): 
Using the cursor, the plot shows that the peak magnification is 2.05 dB, or
= 1.27, occurring at the servo's damped angular natural frequency. This is given by
, i.e.
= 5.4/second (the normalizing constant, T, having been chosen as 1 second). The plot also shows that feedback cutoff, or gain crossover, frequency is given by
, i.e. w = 10.5/second. The peak magnification and feedback cutoff (or gain crossover) frequency agree well with those found from the Nyquist diagram and the Bode attenuation diagram respectively.
The peak magnification and the frequency at which it occurs can be found algebraically as follows. A function C/R(a) is obtained by substituting
in C/R(f). From this function, the maximum value of C/R(f) can then be found thus,
> 
`C/R(alpha)` := map(abs, subs(omega=alpha,`C/R(f)`) ): 
> 
`C/R max` = maximize(`C/R(alpha)`,location,alpha=3..10); 
This shows that the maximum value of C/R(f) is 1.262 and is located at
, demonstrating that using the cursor is sufficiently accurate for practical purposes.
As measured with the cursor, the magnitude of the servo's frequency response falls to 3 dB at
= 15.5/sec. This frequency is known as the servo's cutoff frequency (not to be confused with the feedback cutoff frequency). In this case, the servo's bandwidth extends from the very lowest input frequency  an input of "zero frequency"  up to the cutoff frequency. That is, the bandwidth is 15.5/sec. The precise value can be extracted from C/R(a), thus
> 
omega[c]*T = fsolve( `C/R(alpha)`=10^(3/20), alpha,0..40); 
The greater a servo's bandwidth, the faster its response, but the less the attenuation of any noise in the input.
> 
unassign('omega','alpha'); 
Response to Unit Step Displacement
The servo's response to a unit step function of displacement is given by Output(s) =
. Abbreviating C/R(s) to C/R, the servo's response in the time domain to a unit step function of displacement,
=
, is obtained using the inverse Laplace transform. Since the input is 1 radian, the servo's normalized response is obtained by dividing by 1 radian. The time axis is normalized by dividing by T = 1 sec, so that the plot is completely free of dimensions.
> 
`C/R step disp` := invlaplace(`C/R`/s,s,t): 
> 
p := plot(`C/R step disp`/1,t=0..1,labels=["t/T","C/R"],
title="Normalized Response to Unit Step Displacement", colour=blue): 
> 
q := plot([[0.02,0.08],[0.13,1.0]],style=LINE, colour=red): 
> 
r := plot(1,t=0..1,linestyle=DOT, colour=black): 
From the transient response plot above:
a. the peak of the transient response is 1.21 (i.e. 21% overshoot), which is roughly the same as the peak magnification of the frequency response (1.27).
b. the response time  i.e. time to reach peak of the transient  is 0.3 second.
These can be obtained analytically
> 
zero := diff(`C/R step disp`,t): 
> 
`Max overshoot at t/T` := fsolve(zero=0,t,0.1..0.4); 
> 
`Max overshoot` := eval(`C/R step disp`,t=`Max overshoot at t/T`); 
This demonstrates that using the cursor is accurate enough for most practical purposes.
The maximum slope of the response, judged by eye, is shown by the red line on the plot. The reciprocal of its slope gives an indication of the time for the response, i.e. C/R, to reach 1.0, and is known as the buildup, or rise, time.
> 
`buildup time` := evalf((0.13  0.02)/(1.0  0.09),3); 
This can be obtained analytically,
> 
inflect := diff(`C/R step disp`,t,t): 
> 
tim := fsolve(inflect=0,t,0..0.2): 
> 
`buildup time` := 1/(eval( diff(`C/R step disp`,t), t=tim)); 
Again, using the cursor is accurate enough.
Response to Unit Step of Velocity
The servo's response to a unit step input of velocity is given by Output(s) =
. Abbreviating,
=
. Whence the servo's response in the time domain is obtained using the inverse Laplace transform. Again, the time axis is normalized by dividing it by T = 1 sec, while the response is normalized by dividing it by the input's (angular) position, R. That is,
is divided by
= velocity t = 1 radian/sec t sec = t radian, thus,
> 
`vel step resp` := (invlaplace(`C/R`/s^2,s,t))/ t: 
> 
p := plot(`vel step resp`,t=0..3/2,labels=["t/T","C/R"],
title="Normalized Response to Unit Step of Input Velocity" ): 
> 
q := plot(1,t=0..3/2,linestyle=DOT,colour=black): 
The plot shows that the servo's normalized (position) response to a unit step input of velocity approaches unity smoothly and without any overshoot, reaching 0.99 at t/T = 0.8. This fast, "pursuit" response is ideal for a tracking radar's servo that has to follow targets with rapidly changing velocities.
The response,
, to a unit step input of velocity is simply the normalized response shown above, multiplied by the normalizing factor R = 1 radian/sec t sec, and is shown below.
> 
p := plot(`vel step resp`*t,t=0..1, colour=blue ): 
> 
q := plot(1*t,t=0..1,colour=red): 
> 
display(p,q,title = "Response (blue) to Unit Step Input (red) of Velocity",labels=["t/T","C & R"]);

The input =
= 1 radian/sec t is shown by the straight line on the plot.
The error when responding to a unit step input of velocity is
=
and is shown below
> 
plot( 1000*(1  `vel step resp`)*t, t=0..3/2,title="Error of Response to Unit Step Input of Velocity", labels=["t/T","1000 Error"],labeldirections=[horizontal,vertical] );

The error falls to a minimum of 5.8 milliradians (mils) between t/T = 0.9 and t/T = 1.1, and then increases slowly to its steady state value. Earlier this was shown to be
=
. In this case, W = 1 radian/sec.
> 
`ss error` := evalf(eval(`ss error`,Omega=1)); 
The design requirements included a maximum angular velocity input of 0.5 radian/second, generated by an aircraft flying at 150 yards/second (307 mph), which corresponds to the aircraft passing 300 yards away from a radar antenna controlled by the servo. (See Section 6.12 in Chapter 6 of the Reference). As noted above, the error increases slightly after t/T = 0.9, so the worst case tracking error will be obtained by using the steady state tracking error and the maximum input of 0.5 radian/sec. Since the steady state error is the response to a unit step input of 1 radian/sec, the worst case tracking error with an input of 0.5 radian/second would have been:
> 
`Tracking Error` := evalf(`ss error`/2*(300*yards),3); 
In fact, the maximum tracking error at any range, r(t), is given by
. Noting that the maximum angular input velocity occurs when the aircraft flies on a course perpendicular to the radar's beam,
, whence maximum tracking error =
.