Qualitative Graphical Representation of Nyquist Plots

Download as pdf or txt
Download as pdf or txt
You are on page 1of 20

Qualitative Graphical Representation of Nyquist Plots

R. Zanasia , F. Grossia,∗, L. Biagiottia


a Department of Engineering “Enzo Ferrari”,

University of Modena and Reggio Emilia, via Pietro Vivarelli 10, 41125 Modena, Italy

Abstract

In this paper, the procedure for manually drawing the Nyquist plot of a generic
transfer function is revised and, in particular, two novel parameters (∆τ , ∆p ),
which allow to simplify this process, are presented. Thanks to these parameters,
the analysis of the frequency response at low and high frequencies is considerably
enhanced, with a very little effort. These parameters allow to predict initial and
final directions of the polar curve in the vicinity of initial and final points and
consequently the sectors of the complex plane where the plot starts/ends. In
many cases it is possible to obtain a qualitative Nyquist plot, able to correctly
predict the stability properties of the closed-loop system, by simply joining the
initial and final tracts found with the proposed procedure. Moreover, the anal-
ysis based on these parameters can aid to correctly interpret the plots obtained
with computer programs which often, in particular when poles at the origin are
present, hide the behavior of the frequency response in the area close to the
origin of the complex plane.
Keywords: Nyquist plot, Frequency response, Stability analysis.

1. Introduction

In classical control theory, Nyquist stability criterion plays a central role for
discussing the stability of closed-loop systems. Its effectiveness is strictly related

∗ Corresponding author
Email addresses: [email protected] (R. Zanasi),
[email protected] (F. Grossi), [email protected] (L. Biagiotti)

Preprint submitted to Elsevier June 17, 2015


to the correct drawing of the polar plots of the loop transfer function G(jω).
However, in many applications the exact plot is not essential and a rough sketch
reproducing the shape of the polar plot of G(jω) may be sufficient, provided
that this plot correctly reproduces the encirclements of the loop transfer function
about the critical point −1+j0 and the intersections with the real axis. The rules
for drawing Nyquist plots seem well settled and are illustrated in all controls
textbooks with no substantial differences, see e.g. [1, 2, 3] among many others.
They are based on the analytical computation of the frequency response G(jω)
for ω = 0, ω = ∞ and for those frequencies where the diagram intersects the real
and imaginary axes. In more complex cases, i.e. high-order systems, polar plots
are obtained indirectly, by means of a sort of “translation” of Bode diagrams in
the polar plane. Unfortunately, standard approaches are characterized by some
drawbacks. They may require complex calculations. Moreover, an analysis
based on Bode diagrams (with the phase starting/ending at an integer multiple
of π/2) may lead to the wrong conclusion that the polar plots always start
from or end to a point located in the real or imaginary axis (see for instance
the figures reported in textbook [2], and in particular Fig. 8-33). Indeed, as
highlighted by [4], when the polar plots start at infinity the locus generally does
not approach a coordinate axis but tends to get further and further away from
these axes.
The use of CACSD (Computer Aided Control System Design) programs is the
option that provides the best results in terms of accuracy. But, also in these
cases some problems may arise. In fact, often the obtained polar plots may
result unreadable because of the large span in the magnitude over the entire
frequency range, that hides the local behavior of the curve, in particular in the
region enclosed by the unit circle. This situation is quite common when systems
owning poles on the imaginary axis are considered. In order to cope with this
problem in [5, 6] a logarithmic scaling of the magnitude of the frequency response
is proposed, that allows to magnify the parts of the polar plot close to the
origin without losing the diagram overview. This approach gives good results
in particular with respect to the problem of detecting the intersection with

2
the (negative) real axis, but it is characterized by two important drawbacks:
the proposed plotting technique needs a numerical analysis software (Matlab
functions have been developed by the authors) and it requires to arbitrarily set
the minimum value of the magnitude that can be represented in the diagram.
In this paper, the method for manually drawing polar plots is revised. The
proposed approach has a great value from an educational point of view for a
twofold reason. On the one hand, the definition of two novel parameters (∆τ ,
∆p ) contributes to greatly simplify the drawing procedure and to improve the
correctness of the Nyquist plot, in particular in the regions “far from” and “close
to” the origin. On the other hand, these two parameters and the related con-
siderations can also be the key to correctly interpret many polar plots obtained
with CACSD programs. This allows to correctly analyze the stability of the sys-
tem in a feedback configuration on the basis of the Nyquist criterion. Moreover,
the use of parameters ∆τ and ∆p can be useful to study the dynamic behavior
of a linear system in feedback configuration with a nonlinear static element.
In this case, these parameters used together with well-known methods such as
circle criterion [7], Popov criterion [7] or describing function method [8], may
be conclusive to assess the stability of the nonlinear feedback system.
The paper is organized as follows. In Sec. 2 the proposed parameters for a
qualitative evaluation of Nyquist plots are defined and their meaning explained.
In Sec. 3 the method for drawing Nyquist plots by exploiting the proposed
parameters is illustrated step by step, and some numerical examples are provided
in Sec. 4. Concluding remarks are given in the last section.

2. Qualitative graphical analysis of the frequency response in the


complex plane
Consider a transfer function G(s), without time-delays, expressed as follows
p
Y
βi,mi smi + βi,mi −1 smi −1 + . . . + βi,1 s + βi,0

K i=1
G(s) = q , (1)
sh Y ni ni −1

αi,ni s + αi,ni −1 s + . . . + αi,1 s + αi,0
i=1

3
Pp Pq
where αi,0 6= 0, βi,0 6= 0, i=1 mi = m and i=1 ni = n. Note that from (1) it
is possible to obtain all the standard expressions of rational polynomial transfer
functions, i.e.
bm sm + bm−1 sm−1 + . . . + b1 s + b0
• Polynomial form G(s) =
sh (an sn + an−1 sn−1 + . . . + a1 s + a0 )
(s − z1 ) (s − z2 ) · · · (s − zm )
• Zero-pole-gain form G(s) = ρ
sh (s − p1 ) (s − p2 ) · · · (s − pn )
(1+τ1′ s) (1+τ2′ s) · · · (1+τm

s)
• Time-constant form G(s) = µ
sh (1+τ1 s) (1+τ2 s) · · · (1+τn s)

2.1. Approximating functions

The initial and final points of the polar plot can be obtained by considering
the functions G0 (s) and G∞ (s) that approximate G(s) for |s| ≃ 0+ and |s| ≃ ∞,
respectively. The approximating function G0 (s) is obtained from G(s) in (1) by
neglecting all the s terms except for zeros and poles at the origin:
p
Y
βi,0
µ i=1
G0 (s) = lim+ G(s) = h where µ=K q . (2)
|s|≃0 s Y
αi,0
i=1

The function G∞ (s) is deduced from G(s) in (1) by considering in the numerator
and denominator polynomials only the polynomials of s of higher degree, i.e.
p
Y
βi,mi
ρ i=1
G∞ (s) = lim G(s) = r where ρ=K q (3)
|s|≃∞ s Y
αi,ni
i=1

where r = h + n − m is the relative degree of function G(s). The start-point


and the end-point of the polar plot, denoted respectively with q0 = M0 ejϕ0 and
q∞ = M∞ ejϕ∞ , can be computed by considering the magnitude and the phase
of the frequency responses G0 (jω) and G∞ (jω) for ω → 0+ and ω → ∞:

 ∞
 if h > 0
|µ|  π
M0 = lim+ h = |µ| if h = 0 , ϕ0 = arg (µ) − h
ω→0 ω  2
 0

if h < 0

4
Im
Im
h = −1 ϕ∞

r=2
r = −1
r=0
ϕ0
ϕ0
ϕ0 h=0 ϕ∞ ϕ∞

Re Re
ϕ∞
h=2

h=1
ϕ0 r=1

(a) (b)
Figure 1: Polar plots in the low frequency range when ρ > 0 for different values of h (a) and
plots in the high frequency range when µ > 0 for different values of the relative degree r (b).


 0 if r > 0


|ρ| π
M∞ = lim r = |ρ| if r = 0 , ϕ∞ = arg (ρ) − r
ω→∞ ω  2
 ∞

if r < 0
where arg (µ) and arg (ρ) are 0 or −π according to the sign, positive or negative,
of the constants µ and ρ. Therefore, as well-known, for ω → 0+ polar plots start
from the origin if h < 0, from a point of the real axis if h = 0, or from infinity
if h > 0, see Fig. 1(a). Dually, for ω → ∞ Nyquist plots end to infinity if r < 0,
in a point of the real axis if r = 0 or into the origin if r > 0, see Fig. 1(b).
The initial and final directions of the polar plot are determined by parameters
ϕ0 and ϕ∞ . However, as shown in Fig. 1 for both ω ≃ 0+ and ω ≃ ∞, the
graphical behavior of the polar plot is not univocal because system G(s) may
exhibit a phase lead or a phase lag with respect to ϕ0 and ϕ∞ .

2.2. Low frequency behavior

In order to evaluate the initial phase shift of G(s), it is necessary to take


into account the contribution of all the poles and zeros of the system. From the

5
generic transfer function G(s) given in (1), rewritten as follows
p  
Y βi,mi mi βi,mi −1 mi −1 βi,1
s + s + ... + s+1
µ βi,0 βi,0 βi,0
G(s) = h i=1 q
s Y αi,ni n
 
αi,ni −1 ni −1 αi,1
s +
i
s + ... + s+1
i=1
αi,0 αi,0 αi,0
µ
= G̃0 (s) = G0 (s) G̃0 (s)
sh
it is possible to derive a first order approximation of G(s) for low frequencies
by expanding function G̃0 (s) in Taylor series at s = 0, i.e.
! !
µ dG̃0 (s) 2 µ dG̃0 (s)
G(s) = h G̃0 (s) + s + o(s ) ≃ h 1 + s = G′0 (s)
s s=0 ds s ds
s=0 s=0

Since the derivative of function G̃0 (s) can be expressed as



β β βi,1
dG̃0 (s) Xp
mi βi,m
i,0
i
smi −1 + (mi − 1) i,m i −1
βi,0 smi −2 + . . . + βi,0
= G̃0 (s)  βi,mi m β
ds i=1
β
s i + i,mi −1 smi −1 + . . . + i,1 s + 1
βi,0 βi,0 βi,0
αi,ni ni −1 αi,ni −1 αi,1
q
sni −2 + . . . +
!
X ni αi,0 s + (ni − 1) αi,0 αi,0
− αi,ni n αi,ni −1 αi,1
i=1 αi,0 s
i +
βi,0 sni −1 + . . . + αi,0 s+1

it follows that the approximation of G(s) for low frequencies is


µ
G(s)|s≃0 ≃ G′0 (s) = (1 + ∆τ s) (4)
sh
where
p q
dG̃0 (s) X βi,1 X αi,1
∆τ = = − . (5)
ds i=1
βi,0 i=1
αi,0
s=0
When G(s) is expressed in one of the standard forms reported at the beginning
of the section, the parameter ∆τ can be easily computed as follows
b1 a1
∆τ = − , for polynomial form
b0 a0
m n
X 1 X 1
∆τ = − + , for zero-pole-gain form
i=1 i
z i=1 i
p
Xm n
X

∆τ = τi − τi , for time-constant form.
i=1 i=1

Note that symbol ∆τ comes from the last relation: when G(s) is in time-constant
form the parameter ∆τ is the difference between the time constants τi of the

6
zeros and τi′ of the poles of the transfer function G(s). However, it is also worth
noticing that when the values of the zeros/poles or the equivalent time constants
are not available the computation of ∆τ can be performed by directly using the
coefficients of the polynomial form, see (5).

Remark 1. Consider the transfer function

GT (s) = G(s) T (s) (6)

where G(s) is a rational function, as given in (1), and T (s) is a non-rational


function, e.g. a time-delay or more generally a transfer function describing
an infinite dimensional system in the Laplace domain. The parameter ∆τ for
system GT (s) can be computed as follows
p q
X βi,1 X αi,1
∆τ = − + δτ (7)
i=1
βi,0 i=1
αi,0

where
1 dT (s)
δτ = lim .
ω→0 T (s) s=jω ds s=jω

Relation (7) can be used only if parameter δτ is finite. In the case of a time-delay
system it is T (s) = e−t0 s and δτ = −t0 .

Theorem 1. For ω ≃ 0+ , the phase shift of G(jω) with respect to ϕ0 is given


by ∆ϕ0 = ∆τ ω. Therefore ∆τ > 0 implies an initial phase lead ∆ϕ0 > 0 and
∆τ < 0 implies an initial phase lag ∆ϕ0 < 0.

Proof. The frequency response of the approximating function G′0 (s) in (4) is

µ
G(jω)|ω≃0+ ≃ G′0 (jω) = [1 + jω∆τ ] (8)
(jω)h

and the argument of G(jω) for ω ≃ 0+ is

π
arg G(jω)|ω≃0+ ≃ arg(µ) − h + arctan (∆τ ω)
2
≃ ϕ0 + ∆τ ω (9)

where the approximation arctan(x) ≃ x for x ≃ 0 has been considered. 

7
If ∆τ = 0 the Nyquist plot starts tangentially to the coordinate axis defined
by ϕ0 . In this case, the initial phase shift ∆ϕ0 can be estimated by expanding
the function G̃0 (s) in a third-order Taylor series at s = 0, but, because of the
complexity of the involved computations, this result has no practical interest.

Theorem 2. For system with h > 0, the Nyquist plot starts from infinity.
When h = 1, the plot is tangent to a vertical asymptote whose abscissa is
σ0 = µ ∆τ , while when h > 1 no asymptotes exist.

Proof. From (8) it is possible to conclude that


• if h = 0, limω→0+ G(j0) = µ;
• if h = 1, limω→0+ Re{G(jω)} = µ ∆τ and limω→0+ Im{G(jω)} = ∞;
• if h > 1, limω→0+ Re{G(jω)} = ∞ and limω→0+ Im{G(jω)} = ∞ but no
linear asymptotes exist1 . 

2.3. High frequency behavior

To predict the behavior of G(s) for high frequencies, it is useful to rewrite


the transfer function G(s) as follows
p  
Y βi,0 1 βi,1 1 βi,mi −1 1
+ + . . . + + 1
ρ βi,mi smi βi,mi smi −1 βi,mi s
G(s) = r i=1q 
s Y αi,0 1

αi,1 1 αi,ni −1 1
+ + . . . + + 1
i=1
αi,ni sni αi,ni sni −1 αi,ni s
ρ
= G̃∞ (s) = G∞ (s) G̃∞ (s).
sr

Let G̃∞ (p) = G̃∞ (s)|s= p1 denote the following function

p  
Y βi,0 mi βi,1 mi −1 βi,mi −1
p + p + ... + p+1
i=1
βi,mi βi,mi βi,mi
G̃∞ (p) = q   .
Y αi,0 αi,1 ni −1 αi,ni −1
ni
p + p + ... + p+1
i=1
αi,ni αi,ni αi,ni

1 When h > 1, from (8) it follows that the cartesian expression of the polynomial curve in the
complex plane that approximates G(jω) for ω ≃ 0+ is Re{G(jω)}h = Const · Im{G(jω)}h−1
if h is odd or Im{G(jω)}h = Const · Re{G(jω)}h−1 if h is even, see [4].

8
Since the Taylor series expansion of function G̃∞ (s) at s = ∞ is equal to the
series expansion of function G̃∞ (p) at p = 0, the first order approximation of
function G(s) at s = ∞ can be written as follows
 
 
ρ dG̃∞ (p) ρ ∆p
G(s)|s≃∞ ≃ r G̃∞ (p)|p=0 + p  = r 1− = G′∞ (s) (10)
s dp s s
p=0 p= 1
s

where
q p
dG̃∞ (p) X αi,n i −1
X βi,m i −1
∆p = = − . (11)
dp i=1
αi,ni i=1
βi,mi
p=0

Expression (11) is obtained directly from (5) noting that function G̃∞ (p) shares
the same structure of function G̃0 (s) with inverted coefficients: j → (mi − j) at
numerator and j → (ni − j) at denominator. When G(s) is expressed in one of
the standard forms, the parameter ∆p can be computed as follows:

an−1 bm−1
∆p = − , for polynomial form
an bm
m
X Xn
∆p = zi − pi , for zero-pole-gain form
i=1 i=1
n m
X 1 X 1
∆p = − ′ , for time-constant form.
i=1
τi i=1 τi

In this case the symbol ∆p is due to the fact that ∆p is the difference between
zeros and poles of the transfer function G(s).

Remark 2. Consider the transfer function GT (s) given in (6) and let T (p) =
T (s)|s= p1 . In this case the parameter ∆p in (11) must be modified as follows

q p
dG̃∞ (p) X αi,n i −1
X βi,m i −1
∆p = = − + δp (12)
dp i=1
αi,ni i=1
βi,mi
p=0

where
1 dT (p)
δp = lim . (13)
ω→0 T (p) p=jω dp p=jω

Relation (12) can be used only if parameter δp is finite. For time-delay systems
T (s) = e−t0 s the parameter δp cannot be used because in (13), where T (p) =
t0 t0
dT (p)
e− p and dp = t0
p2 e− p , the limit for ω → 0 does not exist.

9
Theorem 3. For ω ≃ ∞ the phase shift of G(jω) with respect to ϕ∞ is given
1
by ∆ϕ∞ = ∆p ω. Therefore ∆p > 0 implies a final phase lead ∆ϕ∞ > 0 and
∆p < 0 implies a final phase lag ∆ϕ∞ < 0.

Proof. The frequency response of the approximating function G′∞ (s) in (10)
is  
ρ 1
G(jω)|ω≃∞ ≃ G′∞ (jω) = 1 + j ∆p (14)
(jω)r ω
and the argument of G(jω) in the high frequency range can be expressed as
 
π 1
arg G(jω)|ω≃∞ ≃ arg (ρ) − r + arctan ∆p (15)
2 ω
1
≃ ϕ∞ + ∆p . 
ω
Similarly to ∆τ , if ∆p = 0 it is not possible to determine the final phase shift
of G(jω) on the basis of a first-order Taylor series but it is necessary to consider
a third-order approximation. In this case, the polar plot ends tangentially to
the coordinate axis defined by ϕ∞ .

Theorem 4. For system with r < 0, the Nyquist plot ends to infinity. When
r = −1, the plot is tangent to a vertical asymptote with abscissa σ∞ = −ρ ∆p ,
while when r < −1 no asymptotes exist.

Proof. From (14) it is straightforward to conclude that


• if r = 0, limω→∞ G(jω) = ρ;
• if r = −1, limω→∞ Re{G(jω)} = −ρ ∆p and limω→∞ Im{G(jω)} = ∞;
• if r < −1, limω→∞ Re{G(jω)} = ∞ and limω→∞ Im{G(jω)} = ∞ but no
linear asymptotes exist2 . 

3. Qualitative drawing of the Nyquist plot

The qualitative drawing of the Nyquist plot of a generic transfer function


G(s) having the structure given in (1) can be done using the following procedure.

2 When r < −1, from (14) it follows that the cartesian expression of the polynomial curve in
the complex plane that approximates G(jω) for ω ≃ ∞ is Re{G(jω)}r = Const·Im{G(jω)}r+1
or Im{G(jω)}r = Const · Re{G(jω)}r+1 according to the value of r.

10
Im Im
II I ∆p > 0 ∆p < 0
II I
∆τ > 0 ∆τ < 0
ϕ∞ = π2
ϕ0 = π2
σ∞
∆τ < 0 ∆τ > 0 ∆p < 0 ϕ∞ = 0 ∆p > 0
ϕ0 = −π σ0
ϕ0 = 0 ∆τ < 0 Re ∆p > 0 σ∞ ∆p < 0 Re
∆τ > 0
ϕ∞ = −π
ϕ∞ = − π2
ϕ0 = − π2
∆p < 0 ∆p > 0
III IV III IV
∆τ < 0 ∆τ > 0

(a) (b)
Figure 2: Low frequency (a) and high frequency (b) behaviors of the systems shown in Fig. 1
for positive and negative values of parameters ∆τ and ∆p , respectively.

1. Initial point. The initial point of the diagram can be determined by com-
puting magnitude M0 = |G0 (jω)| and phase ϕ0 = arg{G0 (jω)} of the approxi-
mating function G0 (jω) for ω → 0+ .
2. Phase shift (lead or lag) for ω ≃ 0+ . For ω ≃ 0+ the Nyquist plot starts
with a phase shift ∆ϕ0 which is concordant with the sign of parameter ∆τ
defined in (5). Therefore, by using parameter ∆τ it is possible to deduce in
which quadrant the polar plot starts, see Fig. 2(a).
3. Final point. The final point of the diagram can be determined by computing
magnitude M∞ = |G∞ (jω)| and phase ϕ∞ = arg{G∞ (jω)} of the approximat-
ing function G∞ (jω) for ω → ∞.
4. Phase shift (lead or lag) for ω ≃ ∞. For ω ≃ ∞ the Nyquist plot ends with
a phase shift ∆ϕ∞ which is concordant with the sign of parameter ∆p defined
in (11). As shown in Fig. 2(b), by calculating ∆p it is possible to find in which
quadrant of the complex plane the ending segment of the polar plot is contained.
5. Presence of asymptotes. For systems with h = 1 the polar plot exhibits
a vertical asymptote, for ω → 0+ , whose abscissa is σ0 = µ∆τ as stated in
Theorem 2. Dually, for systems with r = −1, the abscissa of the vertical
asymptote for ω → ∞ is σ∞ = −ρ∆p , see Theorem 4. In all the other cases, no

11
asymptotes exist.
6. Total phase variation. When ω varies from 0+ to ∞, the rotation angle of
G(jω) about the origin can be calculated as

π
∆ϕ = − (np,s − np,u − nz,s + nz,u ) (16)
2

where np,s , nz,s , np,u and nz,u denote the number of “stable” and “unstable”
poles and zeros of function G(s), respectively. Note that the computation of
∆ϕ only requires the knowledge of the amount of stable and unstable poles
and zeros but not necessarily their numerical value. As a consequence, if poles
and zeros are not explicitly given, e.g. when G(s) is given in polynomial form,
it is sufficient to apply the Routh-Hurwitz criterion to both numerator and
denominator of the transfer function G(s) to find np,s , nz,s , np,u and nz,u .
7. Analytical computation of G(jω) for some frequency values. A precise draw-
ing of the Nyquist plot needs the exact computation of some peculiar points that
belong to the curve. In particular, the intersections with the real axis are of
great importance for many aims, including stability analysis of linear systems,
stability margins evaluation, stability analysis of systems with static nonlin-
earities, etc. While a direct computation3 of such intersections may lead to
complex calculations, the application of the Routh-Hurwitz stability criterion
to the characteristic equation 1+kG(s) = 0 of the feedback system as a function
of k provides a simple method for obtaining such intersections. As a matter of
fact, if k ⋆ denotes a value of the gain k which makes the system marginally
stable, i.e. it nullifies an element of the first column of the Routh table, an
1
intersection of the Nyquist plot with the real axis occurs in G(jω) = − ⋆ .
k
8. Nyquist plot drawing for ω from 0+ to ∞ . Once that initial and final seg-
ments of the polar curve are known, the Nyquist plot can be obtained by con-
necting them with a continuous curve which performs a rotation of an angle
∆ϕ about the origin and crosses the real axis in correspondence of the points

3 Usually the intersection with the real axis are deduced by computing G(jω ⋆ ) for those
angular frequencies ω ⋆ which guarantee Im{G(jω ⋆ )} = 0.

12
computed in the previous step. Obviously, the exact computation of some
points contributes to improve the precision of the polar plot, but a qualita-
tive sketch/analysis only based on initial and final directions and on the phase
shift ∆ϕ may be sufficient for many goals including stability analysis, see Sec. 4.
9. Complete Nyquist plot. The frequency response for negative values of ω is
the complex-conjugate of function G(j|ω|). Therefore, the Nyquist plot for ω
ranging from −∞ to 0− can be deduced from that of G(jω) obtained for positive
values of ω by reflecting this plot with respect to the real axis. Finally, if the
curve starts from infinity for ω ≃ 0, i.e. h > 0, the complete Nyquist plot
is obtained by connecting point G(j0− ) with point G(j0+ ) with h clockwise
semicircles of infinite radius. Dually, if the curve goes to infinity for ω ≃ ∞, i.e.
r < 0, the complete Nyquist plot is obtained by connecting point G(+j∞) with
point G(−j∞) with r clockwise semicircles of infinite radius.

4. Numerical examples

In order to highlight the significance of the two parameters ∆τ and ∆p , a


few examples, that appear in the literature4 , have been taken into account. For
the sake of simplicity, the considered transfer functions are supposed to be char-

acterized by positive static gain µ and positive time constants τi , τi , but it is
worth noticing that the correctness of the proposed procedure does not depend
on these conditions.

(1 + τ1 s)
1. Consider the transfer function G1 (s) = µ 2
.
s (1 + τ1 s) (1 + τ2 s)
The initial and final points are characterized by
 
µ M 0 = ∞ ′
µτ1 1 M ∞ = 0
G1,0 (s) = 2 → , G1,∞ (s) = →
s ϕ = −π τ1 τ2 s3 ϕ∞ = − 3 π
0
2
The parameters ∆τ and ∆p are
 
′ 1 1 1
∆τ = τ1 − (τ1 + τ2 ) ∆p = − ′ − − − .
τ1 τ1 τ2

4 In particular, some of the transfer functions reported in [3], Tab. 9.6, have been used.

13
Im Im Im

∆p> 0 ∆p> 0 ∆p< 0


ϕ∞ ϕ∞ ϕ∞
ϕ0 ∆τ < 0 ∆τ < 0

1
Re ϕ0 Re ϕ0 Re
∆τ > 0 k⋆

(a) (b) (c)


Figure 3: Nyquist plot of G1 (s) for different combinations of ∆τ and ∆p .

Being ϕ∞ = − 32 π, the value of ∆p which determines the quadrant where the


polar plot reaches the origin does not influence the stability properties of the
feedback system. On the contrary, the value of ∆τ is rather critical: ∆τ < 0
implies that the initial segment of the polar curve is located in the second quad-
rant and therefore no intersections with the negative real axis occur; ∆τ > 0
implies the polar curve starts in the third quadrant and necessarily intersects
the negative real axis, see Fig. 3(a). As a consequence, by applying the Nyquist
stability criterion, it comes out that G1 (s) in a closed-loop configuration with
a proportional regulator, whose gain is denoted by k, is unstable ∀k if ∆τ < 0
and is stable for 0 < k < k ⋆ if ∆τ > 0.

(1 + τ1 s)
2. Consider the transfer function G2 (s) = µ .
s (1 + τ1 s) (1 + τ2 s)
The approximating functions and the initial and final points are:
 
M 0 = ∞ ′ M = 0
µ µτ1 1 ∞
G2,0 (s) = → π , G2,∞ (s) = 2

s  ϕ0 = − τ1 τ2 s  ϕ∞ = −π
2
The parameters ∆τ and ∆p are
 
′ 1 1 1
∆τ = τ1 − (τ1 + τ2 ), ∆p = − ′ − − −
τ1 τ1 τ2
For ∆τ < 0, the polar plot starts from third quadrant along the vertical asymp-
tote σ0 = µ∆τ and the final segment approaches the origin: a) from the third
quadrant when ∆p > 0, see Fig. 4(a); b) from the second quadrant when ∆p < 0,
see Fig. 4(b). The two cases are completely different: in case (a) the system
G2 (s) in a feedback configuration is stable for any k > 0, while in case (b) the

14
Im Im

σ0 ϕ∞ σ0 ∆p < 0
1
∆p > 0 Re − ϕ∞ Re
k⋆

ϕ0 ϕ0
∆τ < 0 ∆τ < 0

(a) (b)
Figure 4: Nyquist plots of G2 (s) for ∆τ < 0 and ∆p > 0 (a) or ∆p < 0 (b).

stability is guaranteed only for 0 < k < k ⋆ . Note that the simple knowledge of
the two parameters ∆τ and ∆p allows to obtain a qualitative plot of the polar
curve that can be used to correctly predict the stability of the feedback system.
The numerical value of k ⋆ can be deduced by means of the Routh-Hurwitz sta-
bility criterion.
It is worth noticing that the plots obtained with the nyquist function of the
Matlab Control Toolbox of the following two transfer functions:
1 1
20( 10 s + 1) 20( 10 s + 1)
G2,a (s) = 1 , G2,b (s) =
s(s + 1)( 20 s + 1) s( 3 s + 1)( 21 s + 1)
1

may be very similar. For function G2,a (s) the values of the two parameters are
∆τ = −0.95 < 0 and ∆p = 11 > 0; for function G2,b (s) the values of the two
parameters are ∆τ = −0.73 < 0 and ∆p = −5 < 0. Only a proper magnification
of the polar plot obtained by Matlab allows to establish whether an intersection
with the negative real axis exists or not, see Fig. 5(c) and 5(d).
In case ∆τ > 0 and ∆p > 0, the vertical asymptote for ω → 0+ is located in
the fourth quadrant, while for ω → ∞ the polar curve approaches the origin
from third quadrant, tangent to the negative real axis, see Fig. 6(a). As a con-
sequence, no intersections with the real axis exist, and the application of the
Nyquist stability criterion leads us to conclude that, under these conditions, the
system G2 (s) with a feedback control is stable ∀k > 0. Note that a simple curve
drawn by joining the initial and the final segments found with the proposed

15
Nyquist Diagram Nyquist Diagram
300 300

200 200

100 100
Imaginary Axis

Imaginary Axis
0 0

−100 −100

−200 −200

−300 −300
−20 −18 −16 −14 −12 −10 −8 −6 −4 −2 0 −20 −18 −16 −14 −12 −10 −8 −6 −4 −2 0
Real Axis Real Axis

(a) (b)
Nyquist Diagram Nyquist Diagram
1 0.4

0.8
0.3

0.6

0.2

0.4

0.1
0.2
Imaginary Axis

Imaginary Axis
0 0

−0.2
−0.1

−0.4

−0.2

−0.6

−0.3
−0.8

−1 −0.4
−2.5 −2 −1.5 −1 −0.5 0 0.5 −2.5 −2 −1.5 −1 −0.5 0 0.5
Real Axis Real Axis

(c) (d)
Figure 5: Matlab plots of the polar curves G2,a (s) (∆τ < 0 and ∆p > 0) (a) and G2,b (s)
(∆τ < 0 and ∆p < 0) (b) and magnification about the origin of the polar curves (c), (d).

methods based on ∆τ and ∆p is able to well reproduce the real diagram, see
10(s+1)
Fig. 6(b) where the polar plot of G2,c (s) = s(s+2)(s+3) obtained with Matlab is
reported. The last case for system G2 (s), that is ∆τ > 0 and ∆p < 0, cannot
occur with any combination of τ1′ , τ1 and τ2 .
Finally, it is worth noticing that, as mentioned in the introduction, the use of
parameters ∆τ and ∆p , together with the well-known methods such as circle
criterion [7], Popov criterion [7] or describing function method [8], may be con-
clusive to assess the stability of the LTI system G(s) when connected in a feed-
back configuration with a static nonlinearity. Consider, for example, the system
G2 (s) in Fig. 7 connected in feedback with the saturation function y = y(x) with
slope m. Without a precise drawing of the Nyquist plot, the circle criterion can-
not be used to study the stability of the system. As a matter of fact, a sufficient
condition for the absolute stability of the feedback system is that the Nyquist
plot of G2 (s) does not touch the critic circle that in this case is the half-plane
1
to the left of the vertical line σ = − m . Unfortunately, this condition cannot be
assured only on the basis of the knowledge of two parameters ∆τ and ∆p and
of the possible intersections with the negative real axis. On the contrary, the

16
PSfrag

Nyquist Diagram

Im 15

10

ϕ∞ σ0

Imaginary Axis
∆p > 0 Re 0

−5

ϕ0 ∆ > 0 −10

τ
−15
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
Real Axis

(a) (b)
Figure 6: Nyquist plot of G2 (s) for ∆τ > 0 and ∆p > 0 (a) and Matlab plot of G2,c (s) =
10(s + 1)
(b).
s(s + 2)(s + 3)

Popov Criterion assures that when ∆p < 0 the system of Fig. 7 is absolutely sta-
ble. In fact, in this case, the Popov plot Re{G2 (jω)} + jω Im{G2 (jω)} does not
intersect the real axis and therefore it is always possible to find a line of proper
1
slope passing through the point − m + j0 which does not touch the Popov plot.
Similar conclusions can also be obtained with the describing function method.
The describing function F (X) of the saturation function y = y(x) in Fig. 7
1
is positive real, and therefore the curve − F (X) , for 0 ≤ X ≤ ∞, belongs to
the real negative axis. If ∆p > 0 the Nyquist plot of G2 (s) does not intersect
the negative real axis and therefore no oscillations can appear in the feedback
system. Conversely, if ∆p < 0 the Nyquist plot of G2 (s) intersects the nega-
tive real axis at point − k1⋆ , and therefore two different situations may occur:
if m > k ⋆ , as shown in Fig. 8(a), an intersection between curves G2,b (jω) and
−1/F (X) exists and a stable persistent oscillation is present in the system, see
Fig. 8(c); on the contrary, if m < k ⋆ no intersections and no oscillations exist
in the feedback system, see Fig. 8(b) and Fig. 8(d).
3. In order to show that the computation of the two parameters ∆τ and ∆p can
be easily performed even if the time-constants of the system, or equivalently the
poles and zeros, are unknown, the transfer function

1200(s + 1/3)(s + 1/2)


G3 (s) =
s (1 + 0.5s) (50s3 + 506s2 + 60.1s + 1)
is considered [5]. Note that G3 (s) is written in a mixed form (time-constants,

17
y(x)
m
x(t) = X sin(ωt) y(t) ≈ Y1 sin(ωt + ϕ1 )
1 G2 (s)
− x

F (X)

Figure 7: Feedback connection of the linear system G2 (s) and the static nonlinearity y(x)
(saturation function) with describing function F (X).

pole-zero and polynomial). By using (5) and (11), the values of ∆τ and ∆p are
 
1 506 1 1
∆τ = 3 + 2 − (0.5 + 60.1) = −55.6 < 0, ∆p = + − + = 11.29 > 0
0.5 50 3 2

Since ϕ0 = − π2 , ϕ∞ = − 3π
2 and ∆ϕ = −π, the polar curve starts from the third

quadrant and ends in the second one, see Fig. 9. Therefore, the simple analysis
based only on initial and final segments leads us to conclude that at least one
intersection with the negative real axis exists, and the plot of Fig. 9(a) would
be the result. Indeed a deeper analysis, based for instance on Routh-Hurwitz
stability criterion, shows that the transfer function with a proportional feedback
control is stable for 0 ≤ k ≤ k1⋆ = 0.0013 and k2⋆ = 0.0713 ≤ k ≤ k3⋆ = 2.861.
As a consequence, see Fig. 9(b), there are three intersections with the real axis:
σ1 = − k1⋆ = −769.23, σ2 = − k1⋆ = −14.02 and σ3 = − k1⋆ = −0.3495. Note
1 2 3

that the Nyquist plot obtained with the Matlab function does not highlight
the presence of intersections with the real axis, see Fig. 10(a), and only with a
progressive enlargement of the plot about the origin it is possible to detect such
points, see Fig. 10(b) and Fig. 10(c).

5. Conclusion

In this paper, the use of two novel parameters ∆τ and ∆p for improving
and simplifying the standard techniques for qualitative drawing of the Nyquist
plot is proposed. The two parameters can be computed with simple calcula-
tions involving only additions and divisions, and in many cases are sufficient to
draw qualitative plots that provide all the information needed for assessing the

18
Im Im
1 1
− ⋆ −
k m
Re 1 Re
1 1 1 −
− − − k⋆
F (X) m F (X)

G(jω) G(jω)
(a) (b)
0.03

0.25

0.2 0.02

0.15

0.1 0.01

0.05
x(t)

x(t)
0 0

−0.05

−0.1 −0.01

−0.15

−0.2 −0.02

−0.25

−0.03
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t t
(c) (d)
Figure 8: Nyquist plot of G2,b (s) and plot of −1/F (X) for m > k ⋆ (a) and m < k ⋆ (b).
Impulse response of the closed-loop system for m = 1 (c) and m = 0.1 (d).

stability of linear and nonlinear systems with a feedback control. Moreover, the
two parameters can be helpful for a correct interpretation of the Nyquist plots
obtained with numerical programs, such as Matlab. As a matter of fact, when
the range of variation of the module of the polar plot is very large the plots
obtained with numerical methods may hide some fundamental details, like the
intersections with the real axis. A simple analysis based on ∆τ and ∆p can
suggest if any intersection exists, and if it is necessary to increase the level of
magnification of the plot to show them.

[1] G. F. Franklin, D. J. Powell, A. Emami-Naeini, Feedback Control of Dy-


namic Systems, 4th Edition, Prentice Hall PTR, Upper Saddle River, NJ,
USA, 2001.

[2] B. C. Kuo, F. Golnaraghi, Automatic Control Systems, 9th Edition, John


Wiley & Sons, Inc., New York, NY, USA, 2002.

19
Im Im

∆p > 0 ∆p > 0
ϕ∞ ϕ∞
σ0 σ0
σ1 σ2 σ3
Re Re

ϕ0 ϕ0
∆τ < 0 ∆τ < 0
(a) (b)

Figure 9: Nyquist plot of G4 (s) for ∆τ < 0 and ∆p > 0

Nyquist Diagram Nyquist Diagram


5 Nyquist Diagram
x 10 1.5
2 80

1.5 60
1

1 40

0.5

0.5 20

Imaginary Axis
Imaginary Axis
Imaginary Axis

0 0 0

−0.5 −20
−0.5

−1 −40

−1
−1.5 −60

−2 −80 −1.5
−12000 −10000 −8000 −6000 −4000 −2000 0 2000 4000 −900 −800 −700 −600 −500 −400 −300 −200 −100 0 100 −20 −15 −10 −5 0
Real Axis Real Axis Real Axis

(a) (b) (c)


Figure 10: Polar plot of G4 (s) obtained with Matlab for different levels of magnification.

[3] R. C. Dorf, R. H. Bishop, Modern Control Systems, 11th Edition, Prentice


Hall PTR, Upper Saddle River, NJ, USA, 2008.

[4] P. Ferreira, Concerning the Nyquist plots of rational functions of nonzero


type, Education, IEEE Transactions on 42 (3) (1999) 228 –229.

[5] T. Andresen, A logarithmic-amplitude polar diagram, Modeling, Identifica-


tion and Control 22 (2) (2001) 65–72.

[6] R. Zanasi, F. Grossi, Open and closed logarithmic nyquist plots, in: Pro-
ceedings of 2014 European Control Conference (ECC), 2014.

[7] H. Khalil, Nonlinear Systems, Prentice Hall PTR, 2002.

[8] K. J. Astrom, R. M. Murray, Feedback Systems: An Introduction for Scien-


tists and Engineers, Princeton University Press, Princeton, NJ, USA, 2008.

20

You might also like