EE C128 / ME134 Problem Set 6 Solution Fall 2011

1 (20 pts) Nyquist Exercise

Consider a close loop system with unity feedback. For each G(s), hand sketch the Nyquist diagram, determine
Z = P − N , algebraically find the closed-loop pole location, and show that the closed loop pole location is
consistent with the Nyquist diagram calculation. Let controller D(s) = k = 2.

a) G(s) = s − 1
D(s)G(s) = 2(s − 1). The open-loop system has no poles, so P = 0.

Now, plot the Nyquist diagram:

Once feedback is added (k = D(s)), from the Nyquist diagram, the point − k1 = −0.5 is inside the
contour and is encircled once CW, therefore N = −1 and Z = P − N = 0 − (−1) = 1. So there is one
zero of the characteristic equation in the RHP, the closed-loop system has one pole in the RHP, and
the system is unstable.
To verify, ∆(s) = 2(s − 1) + 1 = 2s − 1. So the closed loop pole is at 12 , in the RHP.
b) G(s) = s + 1
D(s)G(s) = 2(s + 1). The open-loop system has no poles in the RHP, so P = 0.
 Once feedback is added (k = D(s)), the point − k1 = −0.5 is not inside the contour, therefore N = 0 =
Z = P . So there are no zeros of the characteristic equation in the RHP, and no poles of the system in
the RHP.
To verify, ∆(s) = 2(s + 1) + 1 = s + 3. One pole at -3. No poles or zeros in the RHP, thus N = 0.
This corresponds to the Nyquist diagram where −1 lies outside the contour.
1
c) G(s) = s+1

P = Open loop poles in the RHP = 0

N = CCW encirclements of − k1 = 0
Z = Closed loop poles in the RHP = P − N = 0
Algebraic solution of CL poles: ∆(s) = D(s) + kN (s) = s + 1 + k.
Solving ∆(s) = 0 with k = 2 yields s = −3. Thus, no CL poles in the RHP, and the Nyquist diagram
is consistent.
1
d) G(s) = s−1
 P = Open loop poles in the RHP = 1
N = CCW encirclements of − k1 = 1
Z = Closed loop poles in the RHP = P − N = 0
Algebraic solution of CL poles: ∆(s) = D(s) + kN (s) = s − 1 + k.
Solving ∆(s) = 0 with k = 2 yields s = −1. No CL poles in the RHP, so the Nyquist diagram is
consistent.

2 (20 pts) Nyquist Plot

For controller D(s) = k and
s−2
G(s) =
(s + 1)(s + 2)(s + 3)
a) Hand sketch the asymptotes of the Bode plot magnitude and phase for the open-loop transfer functions.

b) Hand sketch Nyquist diagram.

At s = 0, G(s) = − 13 .
From there, we can see that the Nyquist diagram should go through the point -0.333. We also know
the output contour should travel the same direction since there are poles on the left side and a zero on
the right side. The positive real intercept will occur when ∠G(jω) = 0, which can be estimated from
the Bode plot as ω ≈ 2, which is good enough for sketching purposes.
 To more accurately find the positive real intercept, one would in practice read off the value from the
Matlab Bode plot. The positive intercept can also be found from:
2πn = 0 =∠(jω − 2) − ∠(jω + 1) − ∠(jω + 2) − ∠(jω + 3) (1)
=π + atan(−ω/2) − atan(ω/1) − atan(ω/2) − atan(ω/3) (2)
=π − 2atan(ω/2) − atan(ω) − atan(ω/3) (3)
⇒ ω = 1.8708 (4)
⇒ G(s) = 0.1333 (5)
(Using):
x = fsolve(@(x) pi-2*atan(x/2) - atan(x) - atan(x/3),[0]);
s = j*x;
(s - 2)/((s + 1)*(s + 2)*(s + 3))

c) From Nyquist diagram, determine range of k for stability.

We need to place k such that − k1 is not encircled CW by the Nyquist diagram. A CW encirclement
would imply N = −1 and Z = 1 for the characteristic equation, and thus a closed-loop pole added to
the RHP of the system. Examining the real axis of the system, the range of − k1 is from (a) −∞ to
−0.333 in addition to from (b) 0.1333 to ∞, as these ranges lie outside of the CW encirclement.
Therefore, the range of stable k is (a) 0 to 3 as well as (b) −7.5 to 0. The union of these two sets is
−7.5 < k < 3.
d) Verify sketches with MATLAB and hand in.
3 (20 pts) Nyquist Plot
For controller D(s) = k and
s+1
G(s) =
s2 (s
+ 10)
a) Hand sketch the asymptotes of the Bode plot magnitude and phase for the open-loop transfer functions.

b) Hand sketch Nyquist diagram.

The interesting behavior happens when the Nyquist contour approaches the double pole at the origin.
We can choose to go around the left (dotted) or right (solid) side of the poles to avoid them. If we
take the left path, the outside of the diagram (at infinite magnitude) travels CCW, whereas if we take
the right, the diagram travels CW.

Case 1: Left Side of polese at origin

Case 2: Right Side of poles at origin

 c) From Nyquist diagram, determine range of k for stability.

Case 1: If we took the left (dotted) path, then there are zero crossings at −∞ and (infinitesimally
close to) 0. P = 2 since there we placed the two poles in the Nyquist contour. Between 0 and
−∞ there are two CCW encirclements where N = 2 and thus Z = 0 and the system is stable.
On the positive real axis, there is only one CCW encirclement and thus the closed-loop system
has an unstable pole for k < 0. Therefore, the stable range of k is between 0 and ∞.

Case 2: If we took the right (solid) path, then there is one zero crossing infinitesimally close to 0.
P = 0 since there are no open loop poles in the Nyquist contour. If − k1 ≤ 0, then there are
no encirclements and N = Z = 0, and the system is stable. If − k1 > 0 then there is one CW
encirclement, N = −1, Z = 1, and the system is unstable. These are the same ranges as the first
case, and the range of stable k is between 0 and ∞.
d) Verify sketches with MATLAB and hand in.

The Nyquist plot from MATLAB is misleading since it doesn’t include the behavior close to the zero
poles.

4 (20 pts) Nyquist Plot

For kD(s) = 1 and
(s + 2)2
G(s) =
(s + 1)3 (s + 4)2
a) Hand sketch the asymptotes of the Bode plot magnitude and phase for the open-loop transfer functions.

Normalize the components of the transfer function:

!2  3 !2 !2  3 !2
(s + 2)2
 
s 1 1 1 1 s 1 1
= 2 +1 4 s = 4 +1 s
(s + 1)3 (s + 4)2 2 s+1 4 +1 2 s+1 4 +1

From this we can see:

DC gain is 14 = −12 dB.
There are two zeros with breakpoint at 2.
There are three poles with breakpoint at 1 and two poles with breakpoint at 4.
All of those breakpoints are within one decade of each other, so this diagram is going to be crowded.
b) Hand sketch Nyquist diagram.
For this diagram, we will use the Bode plot information for the section of the test contour that is on the
+jω axis. Note that for this diagram we are using the technique where points that are infinitesimally
close to the origin (b, c, d) are drawn a small distance away so that the phase is visible.

c) From Nyquist diagram, determine range of k for stability.

There are two regions of interest on the negative real axis of the Nyquist diagram. In the drawing
above they are marked with a square and a triangle. We need to find the intersection with the real
axis that divides these two regions. There are several methods we can use.
◦
• Estimate from the Bode plot.  The phase plot shows that ]G(jω) = 180 when ω ≈ 4.
The magnitude plot shows that G(j4) ≈ −40 dB ≈ 0.01. Therefore, the intersection is at
approximately −0.01.
 
• Estimate ω from the Bode plot, then calculate G(jω). As before, ω ≈ 4 from the phase
plot.  
G(j4) = |−0.00885 + 0.00171j| = 0.00892

This answer shows that G(j4) isn’t really that close to the real axis. Our approximation of the
intersection is −0.00892.
• Start with ω = 4 and use trial and error to dial it in.

G(j4) = −0.00885 + 0.00171j start

G(j3) = −0.01640 − 0.00120j crossed the real axis
G(j3.5) = −0.01190 + 0.000837j crossed the real axis
G(j3.3) = −0.01350 + 0.000214j didn’t cross
G(j3.1) = −0.01536 − 0.000651j crossed, but got further away
G(j3.2) = −0.01440 − 0.000184j didn’t cross
G(j3.25) = −0.01394 + 0.000023j crossed, got closer
G(j3.24) = −0.01403 − 0.000017j crossed. close enough!
  
Take the absolute value: G(j3.24) = 0.01403: to finalize the estimate: the crossing is at −0.01403.
• Use MATLAB’s x-y picker to get these values off the Bode or Nyquist plot. Be careful:
unless you do some fancy tricks, your picked values might not be very accurate. I got ω = 3.17,
and the intersection at −0.0147, when I tried it.
• Solve G(jω) = x, with negative real x, for ω.
(s + 2)2
G(s) =
(s + 1)3 (s + 4)2
s2 + 4s + 4
= 5
s + 11s4 + 43s3 + 73s2 + 56s + 16
−ω 2 + 4jω + 4
G(jω) =
jω + 11ω − 43jω 3 − 73ω 2 + 56jω + 16
5 4

Multiply top and bottom by the complex conjugate of the denominator.

−ω 2 + 4jω + 4 −jω 5 + 11ω 4 + 43jω 3 − 73ω 2 − 56jω + 16

=
jω 5 + 11ω 4 − 43jω 3 − 73ω 2 + 56jω + 16 −jω 5 + 11ω 4 + 43jω 3 − 73ω 2 − 56jω + 16
The denominator is something huge, but it’s real, so we discard it. We want the numerator to be
real.

N = jω 7 − ω 6 − 3jω 5 − 55ω 4 − 64jω 3 − 84ω 2 − 160jω + 64

Isolate the terms with j and force them to cancel. Ignore the terms without j.
 
0 = j ω 7 − 3ω 5 − 64ω 3 − 160ω
 
0 = ω ω 6 − 3ω 4 − 64ω 2 − 160

One solution is ω = 0, but we already knew about that one (it’s on the positive real axis).
There’s another solution that we have to find! Since this polynomial is even it’s really a cubic.
Unfortunately, cubics are pretty hard. We can use MATLAB’s roots command to find:

ω = 3.2442357
G(jω) = −0.013996

In summary: depending on the method you use, your value for the intersection on the negative real
axis should be something like −0.014.
Now we need to interpret the Nyquist diagram in these two regions (marked with a square and a
triangle on the diagram above).
Triangle region: − k1 < −0.014:
P = Open loop poles in the RHP = 0
N = CCW encirclements of − k1 = 0
Z = Closed loop poles in the RHP = P − N = 0
Stable.
Square region: 0 > − k1 > −0.014:
P = Open loop poles in the RHP = 0
N = CCW encirclements of − k1 = −2
Z = Closed loop poles in the RHP = P − N = 2
Unstable.
Solving for k, we find that
k ∈ (0, 71.45) ⇒ stable
k ∈ (71.45, +∞) ⇒ unstable

Normally we only consider positive k. If you searched for stability conditions on negative k, you would
find that N = 1 on the real axis before s = 0.25, and N = 0 afterwards. Therefore:
k ∈ (−∞, −4) ⇒ unstable
k ∈ (−4, 0) ⇒ stable
However, this is not a required part of the solution.
 d) Verify sketches with MATLAB and hand in.
Bode Diagram
Nyquist Diagram
0
0.2
Magnitude (dB)

0.15
−50

0.1
−100

Imaginary Axis
0.05

−150 0
0

−0.05
Phase (deg)

−90
−0.1

−180
−0.15

−270 −0.2
−1 0 1 2
10 10 10 10 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3
Frequency (rad/s) Real Axis

5 (20 pts) Gain and phase margin

A closed loop system with unity gain has loop transfer function

125(s + 1)
G(s) =
(s + 5)(s2 + 4s + 25)

a) Plot the Bode magnitude and phase plots for the open loop system (MATLAB ok).
Bode Diagram
Gm = Inf dB (at Inf rad/s) , Pm = 42.1 deg (at 11.4 rad/s)

50
Magnitude (dB)

−50

−100
90
Phase (deg)

−90

−180
−2 −1 0 1 2 3
10 10 10 10 10 10
Frequency (rad/s)
b) Determine the gain and phase margin.
MATLAB’s margin command provides these values, or they can be determined graphically from the
Bode plot.
Gain margin: +∞
Phase margin: 42.1◦

c) Assuming a second order approximation for the closed loop system, estimate the transient response for
a step input from the phase margin and gain margin. (That is estimate ξ, overshoot, peak time, and
settling time.)

Our approach is to find an open-loop second-order system whose phase and gain margins match those
found in part (b), and then estimate the transient response parameters for the closed-loop system from
this approximation. The OL second-order system is characterized by ξol and ωn ol, so our first step is
to estimate those two parameters.

We will examine two approximation methods for the second order system.

• Example 10.13, - 6dB

By definition, the closed loop bandwidth ωn is the frequency where the magnitude of the response
of the closed loop system is -3 dB. Following example 10.13 in Nise, for unity gain feedback, (with
phase in range -135 to -225 degrees), the open loop system would have magnitude response -6 dB
at this ωn . From the Matlab plot, the -6dB frequency ωn ≈ 16.
The CL ξ may be estimated using the “small ξ” approximation discussed in class:
(deg)
ΦM
ξ≈ = 0.421
100
or via the relationship in section 10.10 of Nise:
2ξ
ΦM = tan−1 q p
−2ξ 2 + 1 + 4ξ 4

which proceeds as follows:

q p
tan ΦM −2ξ 2 + 1 + 4ξ 4 = 2ξ
 p 
(0.9022)2 −2ξ 2 + 1 + 4ξ 4 = 4ξ 2
p
1 + 4ξ 4 = (4 + 2 · 0.90222 )ξ 2
1 + 4ξ 4 = 31.674ξ 4
1
ξ4 =
27.674
ξ = 0.4360

Either method is fine; obviously the results are very similar.

So this closed loop estimate has ωn = 16 and ξ = 0.44.

• Approximate open loop 2nd order system
For a second order system, we know that the phase at ωn is 90◦ . Here phase is 90◦ at ωn ≈ 6.3
from the Bode diagram. We know that |G(jωn )| = 2ξ1ol which from the Bode plot is 10 dB ( ≈ 3),
thus ξol ≈ 0.17.
Therefore, our second order approximation of the open loop TF is:
40
s2 + 2.14s + 40
G(s)
Since we are interested in the closed loop step response, we apply the 1+G(s) feedback formula:
 ωn2
CLTF =
ωn2 + (s2 + 2ξωn s + ωn2 )
ωn2
=
s2 + 2ξωn s + 2ωn2
40
= 2
s + 2.14s + 80
√
ω̄n = 80 = 8.9
2.14
ξ¯ = = 0.12
2ω̄n

So when we close the loop around our second order approximation, we get another canonical
second order system, but the DC gain is 0.5 and the ωn and ξ parameters have been changed.
Now we can estimate the transient response parameters from this system. (values in brackets come
from the -6dB BW and Φ100 approximation).
M

ξ = 0.12 [0.44]
√ 2
OS = e−(ξπ/ 1−ξ ) × 100% = 68% [21%]
π
tpeak = p = 0.36 s [0.22 s]
ωn 1 − ξ 2
 p 
− ln 0.02 1 − ξ 2
tsettle = = 3.7 s [0.6 s]
ξωn

d) Compare the actual closed loop step response from MATLAB with the estimates from c).
Remember to find the closed loop step response of the original system.

G(s)
CLTF =
1 + G(s)
125(s + 1)
=
125(s + 1) + (s + 5)(s2 + 4s + 25)
125s + 125
= 3
s + 9s2 + 170s + 250
 Step Response

1.4
True system
Example 10.13
1.2 Second−order approx

1
Amplitude

0.8

0.6

0.4

0.2

0
0 1 2 3 4 5 6
Time (seconds)

From the simulated step response, we can make the following observations:

• The true system and the approximated system have the same steady state behavior
• The peak time is approximately 0.254 s
• The overshoot is approximately 102%
• The settling time is approximately 2.16 s
The true system seems faster than the second-order approximation system; its peak time is significantly
earlier and its overshoot is greater. The second order approximation is not going to be very accurate,
since the phase margin is not small, implying that the 2nd order poles are not that close to the jω
axis, and hence the first order pole and zero are affecting the response.