Chapter 8: Introduction To Systems Control: 8.1 System Stability From Pole-Zero Locations (S-Domain)
Chapter 8: Introduction To Systems Control: 8.1 System Stability From Pole-Zero Locations (S-Domain)
Chapter 8: Introduction To Systems Control: 8.1 System Stability From Pole-Zero Locations (S-Domain)
N (s)
C1
C2
CN
C1
CN
=
+
+ L+
=
+L+
D ( s) s p1 s p 2
s p N s ( 1 + jw1 )
s ( N + jw N )
(8.1)
Rules:
a. If all p i are real then it is obvious that: w1 = w 2 = L = w N 0 and impulse response is simply a sum of real
exponentials:
h(t ) = ( C1e 1t + C 2 e 2t + L + C N e N t ).u (t )
(8.2)
If all i < 0 then the roots (poles) are on negative real axis and each exponential term decays as t
increases. Thus, the system is unconditionally stable for bounded inputs.
If any i > 0 then some of the poles are on the positive real-axis and as a result the term for this
particular pole will diverge as t increases. The system will be generally unstable for any input.
However, it is unconditionally unstable for bounded inputs.
b. If the system has complex roots then they must occur in complex conjugate pairs, i.e., if 1 + jw1 is a root
with a coefficient C1 , then 1 jw1 is also a root with constant C1* . For a simple pole-pair at location:
s k = k jw k the impulse response due to these two terms will be:
h k (t ) = C k .e k t .e jw kt .u (t ) + C k* .e k t .e jwk t .u (t )
= 2 C k .e k t .Cos( w k t + k ).u (t )
(8.3)
Im{C k }
)
Re{C k }
If k < 0 then the system is unconditionally stable for bounded inputs since the pole-zero map is on
the left half plane as before.
where k = a tan(
8.1
If k > 0 then the system is generally unstable since the pole-zero map is sometimes on the right
half plane and other times on the LHP.
Case 2. Simple Poles on jw-axis:
In this case, pole-zero map passes through jw-axis; at that instant k = 0 and the pair of roots are purely
imaginary and hence, the system is oscillatory with an impulse response:
h k (t ) = 2 C k .Cos ( wk t + k )
(8.4)
Rules:
For bounded input signals, the output signal will be bounded oscillatory or simply bounded. But if the
input is also sinusoidal with the same wk , the output will be of the form:
B.s
B
Y (s) = 2
y
(
t
)
=
.t.Sin( wk t )
(8.5)
2 wk
( s + w2 ) 2
k
for k < 0
(8.6)
Rules:
Since e k t will decay faster than the term t m1 increases and the response in (8.6) will be BIBO stable if the
input is bounded.
Case 4. Multiple -Order (Repeated) Poles on jw-axis:
Again let us assume the order of the terms is m, then we will have an impulse response:
hk (t ) = C k .t m1.Cos ( wk t + k )
for k < 0
8.2
(8.7)
Rules:
Here t m1 increases as t grows and the system will be unstable.
Case 5. Multiple -Order (Repeated) Poles in RHP of s-Plane:
The response for these poles will be:
h k (t ) = C k .t m1.e k t .Cos ( wk t + k ) but k > 0
(8.8)
Rules:
Again, both e k t and t m1 will grow as t increases. Expectedly, the system will be unstable.
Example 8.1: Determine if the following feedback system is stable.
Input X(s)
E(s)
Z(s)
s/(s+3)
Y(s)
1/(s+4)
Output
10/(s+10)
8.3
%Step-Input Response
time=0:0.1:10;
response=step(b,a,time);
plot(time,response,'rO');
axis;
As we can see observe from the above pole-zero map, both zeros and poles are all real with z1=0; z2=-10,
p1=-1.388; p2=-4; and p3=-21,6. Therefore, the system is stable for bounded inputs, which is clearly
demonstrated in the step-response plot.
8.4
Example 8.1: The system block diagram and the dynamic equations and the open-loop transfer function of
DC Motor Speed are:
s ( Js + b).( s ) = K .I ( s );
( Ls + R ). I ( s ) = V K .s.( s)
K
=
V ( Js + b).(Ls + R) + K 2
(8.9)
As you can see, the system is overdamped and the settling time is about one second, so the overshoot and
settling time requirements are satisfied. The only problem we can see from this plot is the steady- state error
8.6
of about 50%. If we increase the gain to reduce the steady-state error, the overshoot becomes too large (Try
this yourself). We need to add a lag controller to reduce the steady-state error.
Adding a lag controller
From the plot we see that this is a very simple root locus. The damping and settling time criteria were met
with the proportional controller. The steady-state error is the only criterion not met with the proportional
controller. A lag compensator can reduce the steady-state error. By doing this, we might however increase
our settling time. Try the following lag controller first:
( s + 1) /( s + 0 .1)
(8.10)
This can be done by changing the m-file to look like the following:
J=0.01; b=0.1; K=0.01; R=1; L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];
z1=1; p1=0.01;
numa = [1 z1]; dena = [1 p1];
numb=conv(num,numa);
rlocus(numb,denb)
sgrid(.8,0)
sigrid(2.3)
title('Root Locus with a lag controller')
We get the following root locus, which looks very similar to the original one. Now let's close the loop and see
the closed-loop step response by entering the following code at the end of our m-file:
[k,poles]=rlocfind(numb,denb)
[numc,denc]=cloop(k*numb,denb,-1);
t=0:0.01:3;
step(numc,denc,t)
title('Step response with a lag controller')
Rerun this m-file in the Matlab command window. When prompted to select a point, pick one that is near the damping
requirement (diagonal dotted line). You should get the a plot similar to the following:
8.7
8.8
Our gain should be about 20. As you can see the response is not quite satisfactory.
You may also note that even though the gain was selected to correlate with a position close to the damping
criterion, the overshoot is not even close to five percent. This is due to the effect of the lag controller kicking
in at a later time than the plant. (its pole is slower). What this means is that we can go beyond the dotted
lines that represent the limit, and get the higher gains without worrying about the overshoot
Now we rerunI our m-file, place the gain just above the white, dotted line. Keep trying until you get a
satisfactory response. It should look similar to the following (we used a gain of around 50):
The steady-state error is smaller than 1%, and the settling time and overshoot requirements have been met.
As you can see, the design process for root locus is very much a trial and error process. That is why it is nice
to plot the root locus, pick the gain, and plot the response all in one step. If we had not been able to get a
satisfactory response by choosing the gains, we could have tried a different lag controller, or even added a
lead controller.
8.9
Example: 8.3 Bode Plots: Matlab has a tool for plotting Bode plots, where the frequency is on a logarithmic
scale, the phase is given in degrees, and the magnitude is given as the gain in decibels.
50
H (s) = 3
s + 9 s 2 + 30 s + 40
num=[50]
denom=[1 9 30 40]
bode(num,denom)
8.10
8.11
The time delay can be thought of as an extra block in the forward path of the block diagram that adds
phase to the system but has no effect the gain. That is, a time delay can be represented as a block with
magnitude of 1 and phase w*time_delay (in radians/second).
One nice thing about the phase margin is that you don't need to replot the Bode in order to find the new
phase margin when changing the gains. If you recall, adding gain only shifts the magnitude plot up. This is
the equivalent of changing the y-axis on the magnitude plot. Finding the phase margin is simply the matter of
finding the new cross-over frequency and reading off the phase margin. For example, suppose we use the
above case bode(50,[1 9 30 40]). We will get the following bode plot:
8.12
You should see that the phase margin is about 100 degrees.
Suppose we add a gain of 100, by entering the command bode(100*50,[1 9 30 40]). We should get the
following plot (note I changed the axis so the scale would be the same as the plot above, your bode
plot may not be exactly the same shape, depending on the scale used):
8.13
As you can see the phase plot is exactly the same as before, and the magnitude plot is shifted up by 40dB
(gain of 100). The phase margin is now about -60 degrees. This same result could be achieved if the y-axis
of the magnitude plot was shifted down 40dB. Try this, look at the first Bode plot, find where the curve
crosses the -40dB line, and read off the phase margin. It should be about -60 degrees, the same as the
second Bode plot.
We can find the gain and phase margins for a system directly, by using Matlab. The command margin returns
the gain and phase margins, the gain and phase cross over frequencies, and a graphical representation of
these on the Bode plot. In our problem:
margin(50,[1 9 30 40])
8.14
Bandwidth Frequency
The bandwidth frequency is defined as the frequency at which the closed-loop magnitude response is equal
to -3 dB. However, when we design via frequency response, we are interested in predicting the closed-loop
behavior from the open-loop response. Therefore, we will use a second-order system approximation and say
that the bandwidth frequency equals the frequency at which the open-loop magnitude response is between 6 and - 7.5dB, assuming the open loop phase response is between -135 deg and -225 deg.
In order to illustrate the importance of the bandwidth frequency, we will show how the output changes with
different input frequencies. We will find that sinusoidal inputs with frequency less than Wbw (the bandwidth
frequency) are tracked "reasonably well" by the system. Sinusoidal inputs with frequency greater than Wbw
are attenuated (in magnitude) by a factor of 0.707 or greater (and are also shifted in phase).
8.15
Example: 8.3 Bandwidth from Bode Plots: Consider a second order system transfer function:
H (s) =
1
s 2 + 0.5 * s + 1
Since this is the closed-loop transfer function, our bandwidth frequency will be the frequency corresponding
to a gain of -3 dB.
Looking at the plot, we find that it is approximately 1.4 rad/s.
We can also read off the plot that for an input frequency of 0.3 radians, the output sinusoid should have
a magnitude about one and the phase should be shifted by perhaps a few degrees (behind the input).
For an input frequency of 3 rad/sec, the output magnitude should be about -20dB (or 1/10 as large as
the input) and the phase should be nearly -180 (almost exactly out-of-phase).
8.16
We can use the lsim command to simulate the response of the system to sinusoidal inputs.
1. Consider a sinusoidal input with a frequency lower than Wbw. We must also keep in mind that we want to
view the steady state response. Therefore, we will modify the axes in order to see the steady state response
clearly (ignoring the transient response).
w= 0.3;
num = 1; den = [1 0.5 1 ];
t=0:0.1:100;
u = sin(w*t);
[y,x] = lsim(num,den,u,t);
plot(t,y,t,u) ; axis([50,100,-2,2])
Note that the output (blue) tracks the input (purple) fairly well; it is perhaps a few degrees behind the input as
expected.
2. If we set the frequency of the input higher than the bandwidth frequency for the system, we get a very
distorted response (with respect to the input):
w = 3;
num = 1; den = [1 0.5 1 ];
t=0:0.1:100;
u = sin(w*t);
8.17
[y,x] = lsim(num,den,u,t);
plot(t,y,t,u) ; axis([90, 100, -1, 1])
Note that the magnitude is about 1/10 that of the input, as predicted, and that it is almost exactly out of phase
(180 degrees behind) the input.
Closed-loop performance
In order to predict closed-loop performance from open-loop frequency response, we need to have several
concepts clear:
1. The system must be stable in open loop if we are going to design via Bode plots.
2. If the gain cross over frequency is less than the phase cross over f requency(i.e. Wgc < Wpc), then the
closed-loop system will be stable.
3. For second-order systems, the closed-loop damping ratio is approximately equal to the phase margin
divided by 100 if the phase margin is between 0 and 60 deg. We can use this concept with caution if the
phase margin is greater than 60 deg.
4. For second-order systems, a relationship between damping ratio, bandwidth frequency and settling time
is given by:
In order to make our life easier, we usually plot these equations instead of plugging in the values:
8.18
Similarly, a relationship between damping ratio, bandwidth frequency and rise time is given by:
However, a very rough estimate that you can use is that the bandwidth is approximately equal to the natural
frequency.
8.19
Example: 8.4 Let's use these concepts to design a controller for the following system:
G( s ) =
10
1.25 * s + 1
The design must meet the following specifications: (1) zero steady state error, (2) maximum overshoot must
be less than 40% and (3) settling time must be less than 2 seconds.
There are two ways of solving this problem: one is graphical and the other is numerical. Within Matlab, the
graphical approach is best, so that is the approach we will use. First, let's look at the Bode plot. Create an mfile with the following code:
num = 10;
den = [1.25,1];
bode(num, den)
8.20
Observations: There are several several characteristics of the system that can be read directly from this
Bode plot.
First of all, we can see that the bandwidth frequency is around 10 rad/sec.
Since the bandwidth frequency is roughly the same as the natural frequency, the rise time is
1.8/BW=1.8/10=1.8 seconds. This is a rough estimate, so we can say the rise time is about 2 seconds.
The phase margin for this system is approximately 95 degrees. This corresponds to a damping of
PM/100=95/100=0.95. Plugging in this value into the equation relating overshoot and the damping ratio
(or consulting a plot of this relation), we find that the damping ratio corresponding to this overshoot is
approximately 1%. The system will be close to being overdamped.
The steady-state error can be read directly off the Bode plot as well. The constant (Kp, Kv, or Ka) are
located at the intersection of the low frequency asymptote with the w=1 line. (Just extend the low
frequency line to the w=1 line.) The magnitude at this point is the constant.
Since the Bode plot of this system is a horizontal line at low frequencies (slope = 0), we know this
system is of type zero. Therefore, the intersection is easy to find. The gain is 20dB (magnitude 10).
What this means is that the constant for the error function is 10.
The steady-state error is 1/(1+Kp)=1/(1+10)=0.091. If our system was type one instead of type zero,
the constant for the steady-state error would be found in a manner similar to the following
8.21
Let's check our predictions by looking at a step response plot. This can be done by adding the following two
lines of code into the Matlab command window.
[numc,denc] = cloop(num,den,-1);
step(numc,denc)
Observations:
Our predictions were very good.
The system has a rise time of about 2 s, is overdamped, and has a steady-state error of about 9%.
Now we need to choose a controller that will allow us to meet the design criteria.
We choose a PI controller because it will yield zero steady state error for a step input.
K * (s + a)
GC ( s ) =
s
PI controller has a zero at a , which we can choose. This gives us additional design flexibility to help us meet
our criteria.
1. The first thing we need to find is the damping ratio corresponding to a percent overshoot of 40%:
8.22
Plugging in this value into the equation relating overshoot and damping ratio (or consulting a plot of this
relation), we find that the damping ratio corresponding to this overshoot is approximately 0.28. Therefore, our
phase margin should be approximately 30 degrees. From our Ts*Wbw vs damping ratio plot, we find that
Ts*Wbw ~ 21. We must have a bandwidth frequency greater than or equal to 12 if we want our settling time
to be less than 1.75 seconds which meets the design specs.
2. Now that we know our desired phase margin and bandwidth frequency, we can start our design.
Remember that we are looking at the open-loop Bode plots. Therefore, our bandwidth frequency will be the
frequency corresponding to a gain of approximately -7 dB.
num = [10];
den = [1.25, 1];
numPI = [1];
denPI = [1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);
bode(newnum, newden, logspace(0,2))
8.23
Our phase margin and bandwidth frequency are too small. We will add gain and phase with a zero. Let's
place the zero at 1 for now and see what happens. Change your m-file to look like the following:
num = [10];
den = [1.25, 1];
numPI = [1 1];
denPI = [1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);
bode(newnum, newden, logspace(0,2))
8.24
Now let's try to get a higher bandwidth frequency without changing the phase margin too much. Let's
try to increase the gain to 5 and see what happens .This will make the gain shift and the phase will
remain the same.
num = [10];
den = [1.25, 1];
numPI = 5*[1 1];
denPI = [1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);
bode(newnum, newden, logspace(0,2))
That looks really good. Let's look at our step response and verify our results with the following lines in the
code:
[clnum,clden] =cloop(newnum,newden,-1);
step(clnum,clden)
8.25
As you can see, our response is better than we had hoped for. However, we are not always quite as lucky
and usually have to play around with the gain and the position of the poles and/or zeros in order to achieve
our design requirements .
8.3 System Analysis & design based on Nyquist Diagrams
The Nyquist plot allows us to predict the stability and performance of a closed-loop system by observing its
open-loop behavior. The Nyquist criterion can be used for design purposes regardless of open-loop stability
(remember that the Bode design methods assume that the system is stable in open loop). Therefore, we use
this criterion to determine closed-loop stability when the Bode plots display confusing information.
Note: The Matlab nyquist command does not provide an adequate representation for systems that have
open-loop poles in the jw-axis. Therefore, we suggest that you copy the nyquist1.m file as a new m-file. This
m-file creates more accurate Nyquist plots, since it take into account poles and zeros on the jw-axis.
The Nyquist diagram is basically a plot of G( jw ) , where G(s ) is the open-loop transfer function and w is a
vector of frequencies which encloses the entire right-half plane.
In drawing the Nyquist diagram, both positive and negative frequencies (from zero to infinity) are taken
into account. We will represent positive frequencies in red and negative frequencies in green.
8.26
The frequency vector used in plotting the Nyquist diagram usually looks like this (if you can imagine the
plot stretching out to infinity):
If we have open-loop poles or zeros on the jw-axis, G(s ) will not be defined at those points, and we
must loop around them when we are plotting the contour. Such a contour would look as follows:
Note that the contour loops around the pole on the jw axis. As we mentioned before, the Matlab nyquist
command does not take poles or zeros on the jw axis into account and therefore produces an incorrect
plot. Therefore, If we have a pole on the jw axis, we have to use nyquist1. If there are no poles or zeros
on the jw-axis, or if we have pole-zero cancellation, we can use either the nyquist command or
nyquist1.m.
8.27
G( s )
1 + G( s )
If 1+ G(s ) encircles the origin, then G(s ) will enclose the point -1. Since we are interested in the closed-loop
stability, we want to know if there are any closed-loop poles (zeros of 1+ G(s ) ) in the right-half plane.
Therefore, the lnyquist1.m command plots the Nyquist diagram using a logarithmic scale and preserves the
characteristics of the -1 point.
Example: 8.5 Consider the following transfer functions:
Case 1:
G( s ) =
0.5
s 0.5
8.28
Case 2:
G( s ) =
s+2
s2
Note: this function has a pole at the origin. Let us see the difference between using the nyquist, nyquist1, and
lnyquist1 commands with this particular function.
nyquist([1 2], [1 0 0])
8.29
Note that the nyquist plot is not the correct one, the nyquist1 plot is correct, but it's hard to see what happens close to the -1 point, and
the lnyquist1 plot is correct and has an appropriate scale.
8.30
From the Cauchy criterion that the number N of times that the plot of G(s)H(s) encircles -1 is equal to:
Number Z of zeros of 1 + G(s)H(s) enclosed by the frequency contour
Minus the number P of poles of 1 + G(s)H(s) enclosed by the frequency contour (N = Z - P).
Keeping careful track of open- and closed-loop transfer functions, as well as numerators and denominators,
you should convince yourself that:
the zeros of 1 + G(s)H(s) are the poles of the closed-loop transfer function
the poles of 1 + G(s)H(s) are the poles of the open-loop transfer function.
Note: This is only one convention for the Nyquist criterion. Another convention states that: A positive N
counts the counter-clockwise or anti-clockwise encirclements of -1. The P and Z variables remain the same.
In this case the equation becomes Z = P - N.
It is very important to learn how to count the number of times that the diagram encircles -1. Therefore, we will
go into some detail to help you visualize this. You can view this movie as an example.
Another way of looking at it is to imagine you are standing on top of the -1 point and are following the
diagram from beginning to end. Now ask yourself: How many times did I turn a full 360 degrees? Again, if the
motion was clockwise, N is positive, and if the motion is anti-clockwise, N is negative.
Knowing the number of right-half plane (unstable) poles in open loop (P), and the number of encirclements of
-1 made by the Nyquist diagram (N), we can determine the closed-loop stability of the system. If Z = P + N is
a positive, nonzero number, the closed-loop system is unstable.
Example: 8.6 We can use the Nyquist diagram to find the range of gains for a closed-loop unity feedback
system to be stable.
G( s ) =
s 2 + 10 .s + 24
s 2 8 .s + 15
This system has a gain K which can be varied in order to modify the response of the closed-loop system.
However, we will see that we can only vary this gain within certain limits, since we have to make sure that our
8.32
closed-loop system will be stable. The first thing we need to do is find the number of positive real poles in our
open-loop transfer function:
roots([1 -8 15])
ans =5; 3
The poles of the open-loop transfer function are both positive. Therefore, we need two anti-clockwise (N = -2)
encirclements of the Nyquist diagram in order to have a stable closed-loop system (Z = P + N). If the number
of encirclements is less than two or the encirclements are not anti-clockwise, our system will be unstable.
Let's look at our Nyquist diagram for a gain of 1:
nyquist([ 1 10 24], [ 1 -8 15])
There are two anti-clockwise encirclements of -1. Therefore, the system is stable for a gain of 1. Now we will
see how the system behaves if we increase the gain to 20:
nyquist(20*[ 1 10 24], [ 1 -8 15])
8.33
The diagram expanded. Therefore, we know that the system will be stable no matter how much we increase
the gain. However, if we decrease the gain to 0.5, the diagram will contract and the system might become
unstable.
nyquist(0.5*[ 1 10 24], [ 1 -8 15])
The system is now unstable. By trial and error we find that this system will become unstable for gains less
than 0.80. We can verify our answers by zooming in on the Nyquist plots as well as by looking at the closedloop steps responses for gains of 0.79, 0.80, and 0.81.
8.34
Gain Margin
We already defined the gain margin as the change in open-loop gain expressed in decibels (dB), required at
180 degrees of phase shift to make the system unstable. Now we are going to find out where this comes
from. Suppose that we have a system that is stable if there are no Nyquist encirclements of -1, such as :
50
G( s ) = 3
s + 9.s 2 + 30 .s + 40
Looking at the roots, we find that we have no open loop poles in the RHP and therefore no closed-loop poles
in the RHP if there are no Nyquist encirclements of -1. Now, how much can we vary the gain before this
system becomes unstable in closed loop?
Observations:
1. The open-loop system represented by this plot will become unstable in closed loop if the gain is
increased past a certain boundary.
2. The negative real axis area between -1/a (defined as the point where the 180 degree phase shift
occurs, where the diagram crosses the real axis) and -1 represents the amount of increase in gain that
can be tolerated before closed-loop instability.
3. If the gain is equal to a, the diagram will touch the -1 point:
8.35
G( jw ) =
1
a * .G( jw ) = 1
a
Therefore, we say that the gain margin is 'a' units. Since the gain margin is normally measured in decibels.
Hence, the gain margin in dB is :
GM = 20 * log 10 (a ) dB
Gain margin of the stable, open-loop transfer function we viewed before. Recall that the function is:
50
G( s ) = 3
2
s + 9.s + 30 .s + 40
and that the Nyquist diagram can be viewed by typing:
nyquist (50, [1 9 30 40 ])
To find the gain margin we must find 'a', as defined above. To do this, we need to find the point where there
is exactly 180 degrees of phase shift. This means that the transfer function at this point is real (has no
imaginary part). The numerator is already real, so we just need to look at the denominator. When s = j*w,
the only terms in the denominator that will have imaginary parts are those which are odd powers of s.
Therefore, for G( jw ) to be real, we must have:
8.36
jw 3 + 30. jw = 0
which means w=0 (this is the rightmost point in the Nyquist diagram) or w = 30. We can then find the value
of G( jw ) at this point using a matlab tool called: polyval:
polyval(50,j*w)/polyval([1 9 30 40],j*w)
We now have our gain margin of a = 4.6 and by zooming-in on the Nyquist plot: a = 4.6
nyquist(a*50,[1 9 30 40])
8.37
Phase Margin
We have defined the phase margin as the change in open-loop phase shift required at unity
gain to make a closed-loop system unstable.
From our previous example we know that this particular system will be unstable in closed loop if the Nyquist
diagram encircles the -1 point. However, we must also realize that if the diagram is shifted by theta degrees,
it will then touch the -1 point at the negative real axis, making the system marginally stable in closed loop.
Therefore, the angle required to make this system marginally stable in closed loop is called the phase margin
(measured in degrees). In order to find the point we measure this angle from, we draw a circle with radius of
1, find the point in the Nyquist diagram with a magnitude of 1 (gain of zero dB), and measure the phase shift
needed for this point to be at an angle of 180 deg.
8.38
APPENDIX:
Function nyquist1:
This function is a modified version of the nyquist command of matlab and has been jointly prepared by the
good folks at Carnegie-Mellon and University of Michigans CTM Group 1. It has all the same attributes as the
original, with a few improvements. Nyquist1.m takes poles on the imaginary axis into account when creating
the Nyquist plot, and plots around them. Furthermore, this function outputs the number of open-loop right
hand plane poles, the number of anti-clockwise encirclements, and the number of closed-loop right half plane
poles onto your screen.
Copy the following text into a file nyquis1t.m, and put it in the same directory as the Matlab software, or in a
directory which is contained in Matlab's search path. Be sure to copy from the document source. We have
observed some very peculiar bugs appearing when copying from the text appearing in the browser window.
function [reout,imt,w] = nyquist1(a,b,c,d,iu,w)
%NYQUIST1 Nyquist frequency response for continuous-time linear systems.
%
%
This Verison of the NYQUIST Command takes into account poles at the
%
jw-axis and loops around them when creating the frequency vector in order
%
to produce the appropriate Nyquist Diagram (The NYQUIST command does
%
not do this and therefore produces an incorrect plot when we have poles in the
%
jw axis). As an added feature, this function outputs the number of open loop right
%
hand plane poles, the number of anti-clockwise encirclements, and the number of
%
closed loop right half plane poles on your screen.
%
%
NOTE: This version of NYQUIST1 does not account for pole- zero
%
cancellation. Therefore, the user must simplify the transfer function before using
%
this command.
%NYQUIST1 Nyquist frequency response for continuous-time linear systems.
%
NYQUIST(A,B,C,D,IU) produces a Nyquist plot from the single input
%
IU to all the outputs of the system:
1
This program and many other nice control systems tutorials can be found at: http://www.engin.umich.edu/group/ctm/basic/basic.html
8.39
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
.
x = Ax + Bu
y = Cx + Du
-1
G(s) = C(sI-A) B + D
RE(w) = real(G(jw)), IM(w) = imag(G(jw))
%
J.N. Little 10-11-85
%
Revised ACWG 8-15-89, CMT 7-9-90, ACWG 2-12-91, 6-21-92,
%
AFP 2-23-93
%
WCM 8-31-97
%
% ********************************************************************
% Modifications made to the nyquist - takes into account poles on jw axis.
% then goes around these to make up frequency vector
%
sbegin = 0
else sbegin = w(nbegin);
end;
nend = find(imag(w) > imag(topvalue)); % find all points greater than
nnend = length(nend); % encirclement around jw root
if (nnend == 0)
send = 0
else send = w(nend);
end
w = [sbegin; s1; send]; % reconstituted half of jw axis
end
else
w = sqrt(-1)*w;
end
%end % this ends the loop for imaginary axis poles
% *************************************************************
% back to the regular nyquist program here
% Compute frequency response
if (nargin==2)|(nargin==3)
gt = freqresp(num,den,w);
else
gt = freqresp(a,b,c,d,iu,w);
end
% ***********************************************************
%
%
%
%
nw = length(gt);
mag = abs(gt); % scaling factor added
ang = angle(gt);
mag = log2(mag+1); % scale by log2(mag) throughout
%
%
%
for n = 1:nw
h(n,1) = mag(n,1)*(cos(ang(n,1))+sqrt(-1)*sin(ang(n,1)));
end; % recalculate G(jw) with scaling factor
%
gt = h;
% ***********************************************************
8.44
ret=real(gt);
imt=imag(gt);
% If no left hand arguments then plot graph.
if nargout==0,
status = ishold;
plot(ret,imt,'r-',ret,- imt,'g--')
% plot(real(w),imag(w))
% modifications added here
% ****************************************
% modifications added here to count encirclements
[numc,denc] = tfchk(num,den);
% create the + and - reflection of G(jw) on imag axis
gw = [(ret + j*imt); numc(1)/denc(1); (flipud(ret) - j*flipud(imt))]; %
%look at G(jw)
[Ngw,Mgw] = size(gw); % size of evaluated G(jw)
gwp1 = gw + ones(Ngw,Mgw);
% moves origin from 0 to -1, so we
% can count encirclements around -1 now
initial_angle = rem(180/pi*angle(gwp1(1)) + 360, 360);
% define initial angle
angle_gwp1 = rem(180/pi*angle(gwp1) - initial_angle + 720,360);
% angle from origin for all points
dagw = diff(angle_gwp1);
% subtract off initial angle find degrees of encirclement
tolerance = 180;
% define tolerance - where encirclement counter "flips" over
Ipd = find(dagw < -tolerance);
% number of anti-clockwise encirclements
Ind = find(dagw > tolerance);
% number of clockwise encirclemtents
Nacw = max(size(Ipd)) - max(size(Ind)); % number of encirclements
% Nyquist Criterion says Z = P - N
P = length(find(p>0))
8.45
hold on
plot(real(xy),- imag(xy),'r-',real(xy2),imag(xy2),'g-')
end
xlabel('Real Axis'), ylabel('Imag Axis')
limits = axis;
% Make cross at s = -1 + j0, i.e the -1 point
if limits(2) >= -1.5 & limits(1) <= -0.5 % Only plot if -1 point is not far out.
line1 = (limits(2)-limits(1))/50;
line2 = (limits(4)-limits(3))/50;
plot([-1+line1, -1-line1], [0,0], 'w-', [-1, -1], [line2, - line2], 'w-')
end
% Axis
plot([limits(1:2);0,0]',[0,0;limits(3:4)]','w:');
if ~status, hold off, end % Return hold to previous status
return % Suppress output
end
%reout = ret;
% plot(real(p),imag(p),'x',real(z),imag(z),'o');
8.47