Computer Control Systems - Assignment Example

Topic: COMPUTER CONTROL SYSTEMS Instructions: This is not an essay it is a 100% assignment for a module in control using Matlab and Simulink.I am not certain of the number of pages that will be required to complete this assignment; There is no minimum word requirement. It is a programming assignment. Refrences not required Instructions files attached: 1. matlab_simulink assignments.pdf Created: 2006-02-25 01:38 Style: APA Pages: 10 Deadline: 2006-03-03 23:49 Language Style: English UK Sources: 1 Time Left: 30 hours Grade: 2:1 Classical Control Systems Based upon a mainframe digital computer, Offered the ability to use data storage and retrieval, Alarm functions and process optimization. First installed on a refinery in 1959. Had reliability limitations. 2nd-Order#1 2nd-Order#2 Bode Plots#1 Bode Plots#2 Root Locus#1 Root Locus#2 Design Example#1 Design Example#2 Design Example#3 Design Example#4 Root Locus#2 Design Example#1 Design Example#2 Design Example#3 Design Example#4 Design Example #5 Frequency Domain Phase Lag Design Procedure Frequency Domain Phase Lead Design Procedure Root Locus Transient Performance Design Procedure Root Locus Steady-State Performance Design Procedure Signals & Systems Examples Design Procedures and Examples of Course Concepts ECE 421 Routh Examples The Routh array is a tabular procedure for determining how many roots of a polynomial are in the right-half of the s-plane. We can also determine if there are any roots on the jw axis and their locations. An important use of the Routh array is to determine upper and lower limits on the value of some parameter, such as gain, so that all roots of the closed-loop characteristic equation are in the left-half plane. We will look at four examples. The first three have known numbers for all the polynomial coefficients. The last example will illustrate how to determine parameter value limits. In each example, the open-loop system represented by the transfer functions G(s)H(s) is given, and then the closed-loop characteristics equation is formed. This is the polynomial that the Routh array uses. The first system is represented by the transfer functions and the resulting Routh array. Compare the numbers in the array with the discussion in class on the structure of the array. The number of right-half plane roots of the characteristic equation (closed-loop poles) is given by the number of sign changes in the first column of the array. By inspection, there are no sign changes. Therefore, there are no roots with positive real parts. Since the array was constructed without a 0 appearing anywhere in the first column, there are no roots on the jw axis. The second system is slightly more complex, but the Routh array is formed in the same manner. Note that the number of terms in each row decreases by 1 at each odd-powered row, and that the last element in each even-powered row is the constant coefficient from the characteristic equation. Since there are no sign changes, there are no roots in the right-half plane. There are no roots on the jw axis since there were no 0s in the first column. The third system is the same as the second system except that the gain has been increased by a factor of 10. Note that several of the coefficients in the characteristic equation have changed. Also note that there is a negative coefficient in the polynomial. That guarantees that there is at least one unstable root. Since the constant coefficient is positive, there is an even number of unstable roots. Examination of the first column of the array shows that there are 2 sign changes, from +761.7 to -355.5 and from -355.5 to +120. Therefore, there are 2 roots of the characteristic equation with positive real parts and 2 closed-loop poles in the right-half plane. The location of those roots is not available from the Routh array. Since there were no 0s in the first column, there are no poles on the jw axis. The last system has its gain K left as a variable. We want to determine the upper and lower bounds on K that guarantee that all closed-loop poles are in the left-half of the s-plane. The elements in the Routh array will now be functions of K. The necessary and sufficient condition for all the roots of the characteristic equation to be in the left-half plane is for all the elements in the first column of the Routh array to have the same sign. This puts restrictions on the values of K. The first two rows have positive constants in the first column; the last three rows of the array have K in the first column. Therefore, K must be limited so that those elements are positive. We will look at those rows one at a time, starting with the easy ones, and then combine the results. The s^2 and s^0 rows have simple constraints on the value of K to make the elements positive. They both put upper limits on the allowed value of K. Note that those upper limits are not the same. In order to satisfy both limits, K is restricted to be less than the smallest upper limit. Therefore, based on these two equations, K < 12 is required. The s^1 row is messier, having a quadratic term in K. The limits on K from this term can be found by setting the expression to 0 and solving for the 2 roots. You have to determine which of the roots the upper limit is and which is the lower limit; in some cases it may not be obvious. In this example, 476.4 is the upper limit, and -26.4 is the lower limit. Now all the limits have to be combined. We have upper limits of 300, 476.4, and 12. None of those values can be exceeded by K. Therefore, the actual upper limit on K is 12. The only lower limit encountered in -26.4. Therefore, the limits on K for closed-loop stability is The values that make the s^1 row equal to 0 have a special significance. When the first column of that row is 0, all elements in that row are 0, since there is only 1 element. That means that there are roots of the characteristic equation which are located at points which are symmetric with respect to the origin of the s-plane. These symmetric roots are also roots of the polynomial formed from the elements in the Routh array 1 row above the all-zero row, in this case the s^2 row. This row is called the auxiliary equation. Either value of K that makes the s^1 row 0 can be substituted into the elements of the s^2 row, forming the auxiliary equation Note that only even-powered terms appear in this equation, since this is an even-powered row. When K = 476.4 is substituted into the auxiliary equation, it yields roots at s = +/- j7.2562. When K = -26.4, the roots of the auxiliary equation (which are also roots of the complete characteristic equation) are at s = +/- j1.5339. These are the locations where closed-loop poles cross the jw axis in going from the left-half plane to the right-half plane or vice versa. When K = 12, the coefficient in the s^0 row is 0; this means that there is at least on root of the characteristic equation at the origin of the s-plane. The root locus (studied in detail in Chapter 6) is a plot of the movement of the closed-loop poles as K is varied. Plots for this system are included so that you can see how the closed-loop poles move around as you change the value of K. The open-loop poles and zeros of G(s)H(s) are shown by x and o, respectively. The lines represent the movement of the closed-loop poles. Each line represents one pole. The values of K where a crossing of the jw axis are shown on the plots. The Routh array is often used to calculate these values of K. Root locus for K > 0 Root locus for K < 0 % ***** MATLAB Code Starts Here % function a = routh(x) %ROUTH % % a = routh(x); % % Function to create the Routh array for a nth degree % polynomial to determine the number of roots of the polynomial % with positive real parts. % % x is the polynomial to be tested % a is the Routh array for polynomial x % n = length(x)-1; % degree of x % if mod(n,2) == 0, % no. of columns in array ncol=(n/2)+1; else ncol=(n+1)/2; end % a = zeros(n+1,ncol); % initialize array % for i = 1:ncol % 1st row of array a(1,i) = x(2*i-1); end % for i = 1:ncol-1 % 2nd row of array a(2,i) = x(2*i); end % if mod(n,2) == 1, a(2,ncol) = x(n+1); end % for i = 3:n+1 % rest of the rows for j = 1:ncol-1 a(i,j) = (a(i-1,1)*a(i-2,j+1) - a(i-1,j+1)*a(i-2,1)) / a(i-1,1); end end % % ***** MATLAB Code Stops Here Nyquist Example #1 ECE 421 Nyquist Example #1 The Nyquist plot is a graph of the magnitude and phase of a transfer function evaluated along the jw axis, with the graph displayed as real part vs. imaginary part or magnitude vs. phase. The Nyquist plot contains the same magnitude and phase information as the Bode plot. In the Nyquist plot, however, there is only a single graph, and frequency is not explicitly shown in the plot; it is a parameter along the graph. Perhaps the main use of the Nyquist plot in control system analysis is the application of the Nyquist stability criterion, discussed in later examples. The Nyquist plot and stability criterion can be applied to transfer functions with poles and/or zeros in the right-half of the s-plane. The usual interpretation of a Bode plot limits its application to transfer functions having no poles or zeros in the right-half plane. In that respect, the Nyquist plot is more general than the Bode plot. A simple system is represented by the transfer function Where the gain K = 1. The MATLAB function "nyquist" can be used to compute the real and imaginary parts of the transfer function. The user can specify a set of frequency values or have MATLAB choose the frequencies at which the function is to be evaluated. With no output arguments, "nyquist" will make a plot of the data; with output arguments, no plot is made. The plot linked below is the Nyquist plot for this transfer function. The solid (red) line is for positive frequencies, and the dashed (green) line is for negative frequencies. Three specific frequencies are shown on the graph: w = 0.1 r/s, w = 1 r/s, and w = 10 r/s. The line drawn from the origin to the point on the graph at w = 1 r/s represents the transfer function G(jw) evaluated at w = 1 r/s. The length of the line is the magnitude and the angle that the line makes with respect to the positive real axis is the phase angle. The real and imaginary parts of G(j1) can be read off the axes at that point. Nyquist plot for G1(s) From the plot, notice that the graph starts at a magnitude of 1 and a phase angle of 0 degrees, which is G(j0). As frequency goes to infinity, the magnitude goes to 0 (more poles than zeros in G(s)), and the phase goes to -90 degrees (1 more pole than zero in G(s)). Also note that for negative frequencies, the graph is the mirror image about the real axis of the graph for positive frequencies. The Nyquist plot is always symmetric about the real axis. If the gain in the transfer function is changed, the Nyquist plot is changed in a very simple manner. If the gain is kept positive, then changing the gain changes the magnitude of the transfer function at each frequency, without changing the phase angle. Therefore, the graph just gets larger or smaller, corresponding to increases or decreases in gain. If the gain is made negative, 180 degrees is added to the phase at each frequency, so the entire curve is rotated by that amount. The next figure shows the Nyquist plot for the previous transfer function for gains K of {1, 2, 4, -1, -2, -4}. Only positive frequencies are shown; negative frequencies would give the mirror images of the graphs shown. Effects of gain changes on G1(s) A second system is described by the transfer function With three poles and no zeros, the phase (for w > 0) decreases monotonically from 0 degrees to -270 degrees, and the magnitude decreases monotonically from 1 to 0. For negative frequencies, the mirror image is obtained. Three frequencies are again indicated on the graph to provide a sense of scale to the plot. Notice that for frequencies greater than 10 r/s, the graph is essentially at the orgin. Without making a second plot with an expanded scale, no information is available in that part of the plot. That is one advantage of Bode plots, where both very large and very small magnitudes are visible due to the logarithmic scale used. Nyquist plot for G2(s) Another system is described by the transfer function This transfer function shows the effect of zeros on the Nyquist plot. Each zero (in the left-half plane) has a phase which goes from 0 degrees to +90 degrees as frequency varies from 0 to infinity. The smallest pole or zero in this transfer function is the pole at s=-0.4. Therefore, the Nyquist graph starts having negative phase due to that pole. The three zeros at s=-1 then start providing positive phase, with the result that the phase goes positive and the magnitude increases. The three poles at s=-5 provide negative phase and decreasing magnitude. As frequency goes to infinity, the net result is a magnitude of 0 (more poles than zeros) and a phase of -90 degrees (1 more pole than zero). The frequency of w=2.089 r/s is the frequency at which the maximum positive phase shift occurs, and was obtained from the data. Nyquist plot for G3(s) The last system is described by the transfer function This system also has a collection of poles and zeros. The phase goes negative first since the smallest term is the pole at s=-0.1. The phase then goes positive and the magnitude increases due to the zeros at s=-0.6 and s=-4. Finally, the phase goes negative and the magnitude decreases due to the poles at s = -20, s=-130, and s=-1000. The net result is a magnitude of 0 (more poles than zeros) and a phase of -180 degrees (2 more poles than zeros). Nyquist plot for G4(s) Assuming that the gain is positive and that there are no poles or zeros in right-half of the s-plane, then starting (w = 0) and ending (w = +infinity) points of the graph in a Nyquist plot are determined by two things: the number of poles at the origin the difference between the total number of poles and the number of zeros in G(s). Specifically, let N = the number of poles of G(s) at s=0, n = the total number of poles in H(s) (including those at s=0), and m = the total number of zeros in G(s). The following table gives the relationship between the value of the magnitude curve and those quantities and between the value of the phase curve and those quantities in the limits as frequency goes to 0 and goes to infinity -- under the assumptions of positive gain and no poles or zeros in the right-half of the s-plane. If the assumptions are false, the phase relationships are modified slightly, as described in the next set of examples (no changes in the magnitude relationships). Note that at low frequencies only N is a factor, and at high frequencies (n-m) is the deciding factor. Frequency Value of Magnitude Curve Value of Phase Curve 0 Inf if N > 0 M1 if N = 0 0 in N < 0 - 90*N degrees Inf 0 if n > m M2 if n = m Inf if n < m - 90*(n-m) degrees where M1 and M2 are the magnitudes of the transfer function when it is evaluated at s=0 and s=infinity, respectively, given by MATLAB Code % ***** MATLAB Code Starts Here % w = logspace(-2,3,501); % frequency vector % num1 = 1; den1 = [1 1]; % 1st transfer function [re1,im1] = nyquist(num1,den1,w); % Nyquist data % [i0,i1] = min(abs(w-.1)); % locating frequencies [i0,i2] = min(abs(w-1)); [i0,i3] = min(abs(w-10)); % figure(1),clf,plot(re1,im1,'r',re1,-im1,'g--'),grid % Nyquist plot xlabel('Real Axis'),ylabel('Imag Axis'),... title('Nyquist Plot for 1/(s+1)'),... axis('square'),axis([-0.5 1.5 -1 1]),... % hold on,plot(re1(i1),im1(i1),'wo',... % add frequency markers re1(i2),im1(i2),'wo',re1(i3),im1(i3),'wo',... % and vector [0 re1(i2)],[0 im1(i2)],'w'),hold off,... gtext('w=0.1'),gtext('w=1'),gtext('w=10') % % Show effect of gain on Nyquist plot % figure(2),clf,plot(re1,im1,2*re1,2*im1,4*re1,4*im1,... -re1,-im1,-2*re1,-2*im1,-4*re1,-4*im1),grid,... xlabel('Real Axis'),ylabel('Imag Axis'),... title('Nyquist Plot for K/(s+1)'),... axis('square'),axis([-5 5 -5 5]),... % gtext('1'),gtext('2'),gtext('4'),... gtext('-1'),gtext('-2'),gtext('-4') % num2 = 1; den2 = [1 2 3 1]; % 2nd transfer function [re2,im2] = nyquist(num2,den2,w); % Nyquist data % figure(3),clf,plot(re2,im2,'r',re2,-im2,'g--'),grid % Nyquist plot xlabel('Real Axis'),ylabel('Imag Axis'),... title('Nyquist Plot for G2(s)'),... axis('square'),axis([-0.5 1.5 -1 1]) % hold on,plot(re2(i1),im2(i1),'wo',... % add frequency markers re2(i2),im2(i2),'wo',re2(i3),im2(i3),'wo'),hold off gtext('w=0.1'),gtext('w=1'),gtext('w=10') % olz3 = [-1;-1;-1]; % 3rd transfer function olp3 = [-0.2;-5;-5;-5]; K3 = 2 * abs(prod(olp3) / prod(olz3)); num3 = K3 * real(poly(olz3)); den3 = real(poly(olp3)); % [re3,im3] = nyquist(num3,den3,w); % [i0,i4] = min(abs(w-100)); % locating frequencies [i0,i5] = max(atan2(im3,re3)); % figure(4),clf,plot(re3,im3,'r'),grid % Nyquist plot xlabel('Real Axis'),ylabel('Imag Axis'),... title('Nyquist Plot for G3(s)'),... axis('square'),axis([-1 5 -3 3]) % hold on,plot(re3(i1),im3(i1),'wo',... % add frequency markers re3(i2),im3(i2),'wo',re3(i3),im3(i3),'wo',... % and vectors re3(i4),im3(i4),'wo',re3(i5),im3(i5),'wo',... [0 re3(i1)],[0 im3(i1)],'w',... [0 re3(i2)],[0 im3(i2)],'w',... [0 re3(i5)],[0 im3(i5)],'w'),hold off gtext('w=0.1'),gtext('w=1'),gtext('w=2.089'),gtext('w=10'),gtext('w=100') % olz4 = [-0.6;-4]; % 4rd transfer function olp4 = [-0.1;-20;-130;-1000]; K4 = 2 * abs(prod(olp4) / prod(olz4)); num4 = K4 * real(poly(olz4)); den4 = real(poly(olp4)); % w1 = logspace(-3,5,401); [re4,im4] = nyquist(num4,den4,w1); % figure(5),clf,plot(re4,im4,'r'),grid % Nyquist plot xlabel('Real Axis'),ylabel('Imag Axis'),... title('Nyquist Plot for G4(s)'),... axis('square'),axis([-.5 2.5 -1.5 1.5]),... % hold on,plot(re4(i1),im4(i1),'wo',... % add frequency markers re4(i2),im4(i2),'wo',re4(i3),im4(i3),'wo',... re4(i4),im4(i4),'wo'),hold off,... % gtext('w=0.1'),gtext('w=1'),gtext('w=10'),gtext('w=100') % % ***** MATLAB Code Stops Here ECE 421 Nyquist Example #2 The Nyquist stability criterion is a basic result which can be applied to any transfer function. The stability criterion is applied to the frequency response of the open-loop system G(s)H(s). From the plot of magnitude vs. phase or real part vs. imaginary part for G(jw)H(jw) as w varies from negative infinity to positive infinity, the number of unstable closed-loop poles can be determined. Additional information is also available. We can determine how "close" the closed-loop system is to becoming unstable, and we have some guidance on how to design a compensator to improve the stability characteristics of the system. Thus, more information is available from the Nyquist plot than just the number of unstable closed-loop poles. The Nyquist stability criterion is based on Cauchy's "Principle of Argument" from complex analysis. This principle allows one to determine the difference between the number of zeros and the number of poles in a specified closed region of a complex plane. The complex plane that we are interested in is the s-plane. The items of interest to us are the closed-loop poles of the system. The region of interest in the s-plane is the jw axis and the right-half of the s-plane. This is the region of instability. No closed-loop poles can be allowed to be in this region. The region is defined by the jw axis (from minus to plus infinity) and an infinitely-large semicircle which encloses the entire right-half of the s-plane. The information that is available to us is the open-loop transfer function G(s)H(s) and its frequency response. Assuming the usual negative feedback closed-loop configuration, the closed-loop characteristic equation can be written as where we have not multiplied through by D(s) as we normally do. Thus, the characteristic equation is really a transfer function. The zeros of this transfer function are the zeros of the closed-loop characteristic polynomial and are the closed-loop poles. The poles of the characteristic transfer function are the open-loop poles, that is, the roots of D(s)=0. If we substitute values of s=jw into the transfer function G(s)H(s), we can plot the magnitude of G(jw)H(jw) vs. the phase angle or we can plot the real part of G(jw)H(jw) vs. the imaginary part. Either way, we get the same plot, which defines a complex plane, the G(s)H(s) plane. Plotting the frequency response of G(jw)H(jw) for all w, the Nyquist stability criterion states that the number of encirclements of the point s = -1 + j0 (-1 point) in the G(s)H(s) plane is equal to difference between the number of zeros of 1+G(s)H(s) (the closed-loop characteristic equation) and the number of poles of 1+G(s)H(s). That is, where N0 is the number of encirclements of the -1 point (positive being clockwise; negative being counter-clockwise), Z0 is the number of closed-loop poles in the right-half of the s-plane, and P0 is the number of open-loop poles in the right-half of the s-plane. In order for the closed-loop system to be stable, Z0 must be 0. Therefore, for the closed-loop system to be stable, the Nyquist plot of G(jw)H(jw) must make 1 counter-clockwise encirclement of the -1 point for each unstable open-loop pole of G(s)H(s). Four example transfer functions will be considered. Their Nyquist plots will be made and their closed-loop stability determined by using the Nyquist stability criterion. The four transfer functions are very similar, differing only by the movement of an open-loop pole and/or zero from the left-half to the right-half of the s-plane. Knowing this last piece of information, you should be able to look at the plots and determine which plot goes with which transfer function. The first transfer function and its Nyquist plot are given by: Nyquist plot for system #1 The solid line is for w > 0, and the dashed line is for w < 0. The arrows indicate the direction of movement along the curve for increasing frequency. All the open-loop poles and zeros are in the left-half of the s-plane. Since there are no unstable open-loop poles, P0=0. Therefore, the number of encirclements of the -1 point must be 0 if the closed-loop system is to be stable. From the plot, we can see that the -1 point is not encircled. Therefore, the number of unstable closed-loop poles is 0, and the closed-loop system is stable. N0=Z0=P0=0. Note that the smallest pole is smaller than the smallest zero, and the phase of the Nyquist plot starts out going negative. Also note that there are 4 poles and 2 zeros, and the phase angle of the curve goes to -180 degrees as w goes to infinity. The second system is described by the transfer function and Nyquist plot Nyquist plot for system #2 This system has a zero in the right-half of the s-plane. A right-half plane zero has the normal magnitude response of a zero, but it has the phase shift of a pole, with an additional 180 degrees added to all frequencies. Therefore, for this system there are 5 poles and 1 zero from a phase standpoint. Thus, the phase will go through a total change of 360 degrees as w goes from 0 to infinity. The phase shift starts at 180 degrees and goes in a clockwise direction. The -1 point is encircled one time in the positive direction by the Nyquist plot; N0 = 1. Since the number of unstable open-loop poles P0 = 0, the number of unstable closed-loop poles is Z0 = 1. The closed-loop characteristic equation can be factored to verify this; the unstable root is at s=0.139. The third system is described by the transfer function and Nyquist plot Nyquist plot for system #3 This system has both an open-loop pole and an open-loop zero in the right-half of the s-plane. The unstable pole has the phase shift of a stable zero, and the unstable zero has the phase shift of a stable pole. The magnitudes are unchanged. The 180 degree shift at all frequencies is cancelled out since there is one pole and one zero in the right-half of the s-plane. The phase shift starts at 0 degrees, at it did for the first system, but now it goes positive first due to the reversal in the directions of the first pole and first zero. There is no encirclement of the -1 point. However, there is one unstable open-loop pole. Thus, Z0=1, and there is one unstable closed-loop pole, located at s=0.4072. The last system is described by the transfer function and Nyquist plot Nyquist plot for system #4 This system has one unstable open-loop pole but no unstable open-loop zeros. In order for the closed-loop system to be stable, there must be 1 counter-clockwise encirclement of the -1 point by the Nyquist plot. The phase starts at 180 degrees since there is one unstable pole or zero. The system now has the equivalent of 3 poles and 3 zeros since the phase shift of the unstable pole will be that of a zero. Now the two smallest terms in the transfer function are zeros from a phase perspective. Therefore, the curve goes in a counter-clockwise direction. The Nyquist plot makes 1 encirclement of the -1 point in a negative direction. Therefore, the closed-loop system is stable. Knowing the relative ordering of the values of the poles and zeros, and the sign of the gain, figure out how to look at the Nyquist plots and determine which plot goes with which transfer function. Nyquist Example #2 MATLAB Code % ***** MATLAB Code Starts Here % % System 1 -- all open-loop poles & zeros in LHP % olz1 = [-0.5;-10]; olp1 = [-0.3;-5;-20;-20]; num1 = 250 * real(poly(olz1)); den1 = real(poly(olp1)); figure(1),clf,nyquist(num1,den1);title('System #1') % % System 2 -- one open-loop zero in RHP, all open-loop poles in LHP % olz2 = [0.5;-10]; olp2 = [-0.3;-5;-20;-20]; num2 = 250 * real(poly(olz2)); den2 = real(poly(olp2)); figure(2),clf,nyquist(num2,den2);title('System #2') % % System 3 -- one open-loop zero and one open-loop pole in RHP % olz3 = [0.5;-10]; olp3 = [0.3;-5;-20;-20]; num3 = 250 * real(poly(olz3)); den3 = real(poly(olp3)); figure(3),clf,nyquist(num3,den3);title('System #3') % % System 4 -- one open-loop pole in RHP, all open-loop zeros in LHP % olz4 = [-0.5;-10]; olp4 = [0.3;-5;-20;-20]; num4 = 250 * real(poly(olz4)); den4 = real(poly(olp4)); figure(4),clf,nyquist(num4,den4);title('System #4') % % Form the closed-loop characteritics equations % dcl(1,:) = den1 + [0 0 num1]; dcl(2,:) = den2 + [0 0 num3]; dcl(3,:) = den3 + [0 0 num3]; dcl(4,:) = den4 + [0 0 num4]; % % Find the closed-loop poles % for i = 1:4 clp(:,i) = roots(dcl(i,:)); end % clear i % % ***** MATLAB Code Stops Here ECE 421 Nyquist Example #3 Gain margin and phase margin are measures of relative stability. They measure how "close" a system is to crossing the boundary between stable and unstable, in one direction or the other. Gain margin is the amount of change in the value of the gain of the transfer function, from its present value, to that value that will make the Nyquist plot pass through the -1 point in the G(s)H(s) plane. The gain margin is a ratio of the value of the gain needed to have the curve pass through the -1 point to the current value of the gain. The phase margin is the amount of pure phase shift (no change in magnitude) that will rotate the Nyquist curve so that it passes through the -1 point. The pure phase shift would be due to a pure time delay in the time domain. Both gain and phase margin can be seen by looking at the closed-loop characteristic equation evaluated on the jw axis. If a closed-loop pole of the system is on the jw axis, then the magnitude equals 1 and the phase is 180 degrees at that frequency. If that occurs, the Nyquist plot will pass through the -1 point at that frequency. If there are no closed-loop poles on the jw axis, the plot will not pass through the -1 point, but we can use the gain and phase margins as measures of how much avoidance there is of the -1 point. Two special frequencies are defined. The frequency at which the magnitude of G(jw)H(jw) equals 1 is called the gain crossover frequency, w_x. The frequency at which the phase shift of G(jw)H(jw) equals 180 degrees is called the phase crossover frequency, w_phi. The gain margin is defined only at the frequency w_phi. At this frequency, the phase shift is 180 degrees. A gain change will not rotate the curve, only make it larger or smaller. From a graphical perspective, the gain margin is the reciprocal of the magnitude of G(jw)H(jw) where it crosses the negative real axis. In the first figure, the unit circle represents all points where the magnitude is 1, and the negative real axis represents all points where the phase shift is 180 degrees. The gain margin, GM, is 1 divided by the magnitude of G(jw)H(jw) at the point where the negative real axis is crossed. In the figure, the crossing occurs at -1/3, so the magnitude is 1/3, and the gain margin is 3. The phase margin is defined only at the frequency w_x. At this frequency, the magnitude is 1. This is shown graphically by the point where the Nyquist plot crosses the unit circle. The phase margin is the amount by which that point of intersection can be rotated so that it coincides with the -1 point. This represents a pure phase change, with no change in the magnitude. In the figure, the phase shift where G(jw)H(jw) crosses the unit circle is -135.5 degrees, so the phase margin is +44.5 degrees. Defining gain and phase margins An example transfer function will be considered with three values of gain. The open-loop system is given by The Nyquist plots for the three gains are shown in the next 3 figures. Note that the shape of the three plots is the same (note the different scales on the plots), but their relationships with the unit circle have changed. Nyquist Plot, Gain = 20 Nyquist Plot, Gain = 40 Nyquist Plot, Gain = 80 An expanded view of all three plots are shown in the next figure, along with the unit circle. Here it can be clearly seen that the gain margin and phase margin are both reduced as the gain is increased. The number of unstable open-loop poles is 0 for the system. For gains of K=20 and K=40, there are no encirclements of the -1 point: P0=0, N0=0, and Z0=0. For the case of K=80, the -1 point is encircled: P0=0, N0=2, Z0=2. Therefore, for the first two gains, the closed-loop system is stable; for the third gain, the closed-loop system is unstable. The maximum value of K for closed-loop stability is K=60. Expanded View of Nyquist Plots For the two smaller gains, the curves cross the negative real axis at magnitudes less than 1, and they cross the unit circle in the 3rd quadrant. Therefore, the gain margins are greater than 1, the phase margins are positive, and the gain crossover frequencies are less than the phase crossover frequencies. For these conditions, the closed-loop systems are stable. For the largest gain shown, the curve crosses the negative real axis at a magnitude greater than 1 and crosses the unit circle in the 2nd quadrant. Therefore, the gain margin is less than 1, the phase margin is negative, and the gain crossover frequency is greater than the phase crossover frequency. The closed-loop system is unstable. Gain PM (deg) GM w_x (r/s) w_phi (r/s) 20 44.5 3 1.84 3.32 40 13.9 1.5 2.73 3.32 80 -8.6 0.75 3.77 3.32 MATLAB Code Nyquist Example #3 MATLAB Code % ***** MATLAB Code Starts Here % K = [20 40 80];% system numerators den = [1 6 11 6];% system denominator % [GM(1),PM(1),wcgm(1),wcpm(1)] = margin(K(1),den); % stability margins [GM(2),PM(2),wcgm(2),wcpm(2)] = margin(K(2),den); [GM(3),PM(3),wcgm(3),wcpm(3)] = margin(K(3),den); % [re1,im1] = nyquist(K(1),den);% frequency responses [re2,im2] = nyquist(K(2),den); [re3,im3] = nyquist(K(3),den); % [xx,yy] = circle;% unit circle % figure(1),clf,plot(re1,im1,'r',xx,yy,'w',[-5 5],[0 0],'w--',... [0 0],[-5 5],'w--',[-1/GM(1) -1/GM(1)],[-0.05 0.05],'y',[0 -2],[0 -2*tan(PM(1)*pi/180)],'g'),... axis('square'),axis([-2 0 -1 1]),xlabel('Real Axis'),... ylabel('Imag Axis'),title('Definitions of Gain and Phase Margin'),... text(-.8,.5,'unit circle'),... text(-1/GM(1)-0.1,0.1,'GM'),text(-1.1,-0.5,'PM') % figure(2),clf,plot(re1,im1,'r',re1,-im1,'g--',xx,yy,'w',[-5 5],[0 0],'w--',... [0 0],[-5 5],'w--'),axis('square'),axis([-2 4 -3 3]),xlabel('Real Axis'),... ylabel('Imag Axis'),title('G(s) = 20 / [(s+1)(s+2)(s+3)]') % figure(3),clf,plot(re2,im2,'r',re2,-im2,'g--',xx,yy,'w',[-5 5],[0 0],'w--',... [0 0],[-5 5],'w--'),axis('square'),axis([-3 7 -5 5]),xlabel('Real Axis'),... ylabel('Imag Axis'),title('G(s) = 40 / [(s+1)(s+2)(s+3)]') % figure(4),clf,plot(re3,im3,'r',re3,-im3,'g--',xx,yy,'w',[-15 15],[0 0],'w--',... [0 0],[-15 15],'w--'),axis('square'),axis([-5 15 -10 10]),xlabel('Real Axis'),... ylabel('Imag Axis'),title('G(s) = 80 / [(s+1)(s+2)(s+3)]') % figure(5),clf,plot(re1,im1,'y',re2,im2,'r',re3,im3,'g',xx,yy,'w',[-2 0],[0 0],'w--'),... axis('square'),axis([-2 0 -1 1]),xlabel('Real Axis'),... ylabel('Imag Axis'),title('Comparison of Nyquist Plots') % text(-.65,-.5,'K = 20'),... text(-1.25,-.8,'K = 40'),... text(-1.75,-.35,'K = 80'),... text(-.8,.5,'unit circle') % % ***** MATLAB Code Stops Here ECE 421 Nyquist Example #4 Three sets of plots are given in this example. The first 2 sets each have four plots, and the third set has three plots. The first two sets show the effects on the Nyquist plot and on the stability of the closed-loop system of moving a pole or a zero from the left-half to the right-half of the s-plane and the effect of changing the sign of the gain K. The third set shows the effects on the Nyquist plot and on closed-loop stability of changing the value of the gain. G1a(s), all poles & zeros in LHP, K > 0, N0 = 0, P0 = 0, Z0 = 0 G1b(s), 1 zero in RHP, K < 0, N0 = 0, P0 = 0, Z0 = 0 G1c(s), 1 zero in RHP, K > 0, N0 = 1, P0 = 0, Z0 = 1 G1d(s), 1 pole in RHP, K < 0, N0 = 0, P0 = 1, Z0 = 1 G2a(s), all poles & zeros in LHP, K > 0, N0 = 0, P0 = 0, Z0 = 0 G2b(s), 1 zero in RHP, K < 0, N0 = 2, P0 = 0, Z0 = 2 G2c(s), 1 zero in RHP, K > 0, N0 = 1, P0 = 0, Z0 = 1 G2d(s), 1 pole in RHP, K < 0, N0 = 0, P0 = 1, Z0 = 1 G3a(s), K = 1000, N0 = 0, P0 = 2, Z0 = 2 G3b(s), K = 3000, N0 = -2, P0 = 2, Z0 = 0 G3c(s), K = 9000, N0 = -1, P0 = 2, Z0 = 1 ECE 421 Second-Order System Example This example examines the effects that a physical parameter, amplifier gain in this case, has on the damping ratio and undamped natural frequency of the standard 2nd-order system. The open-loop system consists of an ideal amplifier, with gain KA, in series with a DC servo motor rigidly connected to an inertial load. The inertia of the motor and load are lumped together in this case, as is the viscous friction of the motor and load. The only deviation from the standard DC motor is that the armature inductance is assumed to be 0. This yields the standard second-order system model. The following Table provides the parameter values for the open-loop system, and the relevant equations are shown below. Table of Parameter Values 2nd-Order System Example Table of Parameter Values Variable Description Units Value R Motor Armature Resistance Ohms 1 L Motor Armature Inductance Henrys 0 Kb Motor Back EMF Constant Volts / rad/sec 100 KT Motor Torque Constant ft-lbs / amp 139 J Motor/Load Rotational Inertia ft-lbs-s^2 44,200 B Motor/Load Viscous Damping ft-lbs-s 7,000 KA Amplifier Gain volts/volt Variable The output signal of the given system is the angular position of the motor and load. Unity feedback is placed around the system, with the actual position being compared with the desired position to produce the position error. This error signal is the input to the open-loop system G(s). The amplifier in the system converts this position error to a voltage signal which is applied as the armature voltage of the motor. The motor turns in a direction to reduce the position error. The open-loop transfer function G(s) is shown both in terms of the physical variables as well as the standard variables of the 2nd-order system model. The un-damped natural frequency and damping ratio can be easily computed by comparing the two forms of G(s). For this example, the value of the amplifier gain KA will be varied from 20 to approximately 6300 in a logarithmic manner. This provides a well-distributed range of values for the damping ratio from just less than 1 to just greater than 0. The first two plots are for the standard second-order system parameters, damping ratio and un-damped natural frequency, both plotted versus the amplifier gain. As the first plot illustrates, the system is under-damped for the entire range of KA values chosen, but would become critically damped and over-dampled for smaller values of KA. If KA is increased beyond 6300, the system approaches being undamped. The undamped natural frequency increases with increasing KA, as shown in the second plot. Therefore, the amplifier gain affects both of the parameters in the second-order system. Damping Ratio Undamped Natural Frequency The next plot shows percent overshoot to a step input and phase margin. These are both measures of relative stability, one from the time domain and one from the frequency domain. Overshoot and phase margin are functions of the damping ratio only. Large values of overshoot are generally considered bad, and small values of phase margin are generally considered bad. The figure illustrates that both of those measures of relative stability are improved by increasing the value of the damping ratio, which corresponds to decreasing the value of the amplifier gain KA. For larger values of KA, the closed-loop system is "less stable", although it is still stable in an absolute sense. Percent Overshoot and Phase Margin The next plot shows three time measurements of the closed-loop step response -- rise time (10% - 90%), settling time (2%), and peak time. Peak time and rise time are both seen to decrease with increasing KA. The decreases in rise and peak times are caused by the combined effects of decreasing damping ratio and increasing undamped natural frequency. Thus, if the system becomes "less stable" (smaller zeta) or develops a higher natural frequency, the times associated with responding to the discontinuity in the step input become shorter. The interesting graph in the figure is settling time. It is completely independent of the amplifier gain, depending only on the product of damping ratio and undamped natural frequency. The amplifier gain affects those two parameters in equal and opposite directions so that their product remains equal. Thus, the total time for the transient response to decay does not depend on the amplifier gain. Rise, Settling, and Peak Times The next two plots show the closed-loop bandwidth and the open-loop gain crossover frequency and their ratio. Increasing the gain is shown to increase the bandwidth and the gain crossover frequency. The wider bandwidth means that the closed-loop system will be able to respond to more rapidly changing reference input signals, in this case, desired values for the output position. The plot of the ratio of bandwidth to gain crossover frequency shows that the ratio is fairly constant. A usable rule-of-thumb for underdamped second-order systems is that the closed-loop bandwith is approximately 1.5 times the open-loop gain crossover frequency. Bandwidth and Gain Crossover Frequency Ratio of Bandwidth to Gain Crossover Frequency The last plot shows the product of the step response rise time and the closed-loop bandwidth. For the range of values for KA used in this example, the ratio of largest to smallest rise time is approximately 55. The ratio of largest to smallest value for the bandwidth is approximately 39. However, for the product of rise time and bandwidth, the ratio of largest to smallest value is only 1.5!! Therefore, we can see that there is a fundamental relationship between the range of frequencies that the closed-loop system can respond to without significant attenuation and the speed with which the closed-loop system can respond to an abrupt change in the input signal. Product of Rise Time and Bandwidth Second-Order System Example MATLAB Code % ***** MATLAB Code Starts Here % R = 1; % motor parameter values KT = 139; Kb = 100; % J = 44200; % load inertia and damping K = sqrt(3.8201e+04 * 1.6978e+05); B = 0.5 * (26*sqrt(K/10) + 130*sqrt(K/10)); % KA = 2 * logspace(1,4,19); % range of amplifier gains KA = KA(1:16); % wn = sqrt(KA*KT ./ (J*R)); % undamped natural frequency z = ((B*R+Kb*KT) ./ (J*R)) ./ (2*wn); % damping ratio % % Calculate percent overshoot (po), rise time (tr), settling time (ts), peak time (tp), % gain crossover frequency (wx), phase margin (pm), and bandwidth (wb), % for i = 1:length(z) po(i) = 100 * exp(- pi * z(i) / sqrt(1 - z(i)^2)); tr(i) = (1 + 1.1 * z(i) + 1.4 * z(i)^2) / wn(i); ts(i) = 4 / (z(i) * wn(i)); tp(i) = pi / (wn(i) * sqrt(1 - z(i)^2)); wx(i) = wn(i) * sqrt(- 2 * z(i)^2 + sqrt(4 * z(i)^4 + 1)); pm(i) = (180 / pi) * atan(2 * z(i) * wn(i) / wx(i)); wb(i) = wn(i) * sqrt(1 - 2 * z(i)^2 + sqrt(4 * z(i)^4 - 4 * z(i)^2 + 2)); end % % Plot the results versus amplifier gain % figure(1),clf,semilogx(KA,z,KA,z,'wo'),grid,... xlabel('KA (v/rad)'),ylabel('Zeta'),title('Damping Ratio vs. Amplifier Gain') % figure(2),clf,semilogx(KA,wn,KA,wn,'wo'),grid,... xlabel('KA (v/rad'),ylabel('wn (r/s)'),title(' Undamped Natural Frequency vs. Amplifier Gain') % figure(3),clf,semilogx(KA,po,KA,po,'wo',KA,pm,KA,pm,'wo'),grid,xlabel('KA (v/rad)'),... ylabel('PO (%) & PM (deg)'),title('Overshoot & Phase Margin vs. Amplifier Gain'),... gtext('Overshoot'),gtext('Phase Margin') % figure(4),clf,semilogx(KA,tr,KA,tr,'wo',KA,ts,KA,ts,'wo',KA,tp,KA,tp,'wo'),grid,... xlabel('KA (v/rad)'),ylabel('Tr, Ts, & Tp (s)'),... title('Rise, Settling, & Peak Times vs. Amplifier Gain') gtext('Rise Time'),gtext('Settling Time'),gtext('Peak Time') % figure(5),clf,semilogx(KA,wb,KA,wx,KA,wb,'wo',KA,wx,'wo'),grid,xlabel('KA (v/rad)'),... ylabel('wb & wx (r/s)'),title('Bandwidth & Gain Crossover vs. Amplifier Gain'),... gtext('Bandwidth'),gtext('Gain Crossover Freq') % figure(6),clf,semilogx(KA,wb./wx,KA,wb./wx,'wo'),grid,xlabel('KA (v/rad)'),... ylabel('wb / wx'),title('Ratio of Bandwidth to Gain Crossover vs. Amplifier Gain') % figure(7),clf,semilogx(KA,tr.*wb,KA,tr.*wb,'wo'),grid,xlabel('KA (v/rad)'),... ylabel('tr * wb'),title('Product of Rise Time & Bandwidth vs. Amplifier Gain') % % ***** MATLAB Code Stops Here Read More