MATLAB Session 7 Results

First a reminder following Session 6 about the availability of the nice SISO graphics design tool:

>> sisotool

The default view includes both the root locus and Bode plots. An illustration after importing a third order function is shown below. As you vary the position of the closed-loop poles (the red squares in the root locus plot), the Bode plot will be updated instantaneously (adjusting to the corresponding loop gain and reporting the new gain and phase margins).

We don't need to say more as the graphics tool is well designed and intuitive. Use the on-line Help if in case you have the need to.


Now back to the collection of MATLAB statements and screen display:

>> %Polar (Nyquist) plots
>> p=poly([-1; -2; -3]);
>> G=tf(10,p);
>> nyquist(G);

>> %We'll do our own plot, mapping only the positive imaginary axis
>> [re,im,w]=nyquist(G);
>> plot(re(1,:),im(1,:))
>> w=logspace(-1,1);
>> [re,im]=nyquist(G,w);
>> plot(re(1,:),im(1,:))
>> hold
Current plot held
>> x=-1; y=0;  %Add the (-1,0) point and the axes
>> xh=[-2 2]; yh=[0 0];
>> xv=[0 0]; yv=[-2 1];
>> plot(x,y,'o',xh,yh,'-',xv,yv,'-')

>> %Add the curves of much larger gains
>> G=tf(50,p);
>> [re,im]=nyquist(G,w);
>> plot(re(1,:),im(1,:))
>> G=tf(60,p);
>> [re,im]=nyquist(G,w);
>> plot(re(1,:),im(1,:))
>> %Just a sample Nichols plot
>> p=poly([-1; -2; -3]);
>> G=tf(10,p);
>> ngrid    %New MATLAB requires generating the grid first
>> nichols(G)
>> zoom    %Instead of using axis(), zoom is another way to cheat
>> %Magnitude and Phase Angle (Bode) Plots
>> G=tf(1,[1 0.4 1]);
>> bode(G)
>> %If we do not use dB as the unit for magnitude, we'll need to do our own plots:
>> %This is not so bad if you put all these statements in an M-file
>> %You actually must learn how to do this for problems with dead time
>> w=logspace(-1,1);
>> [mag,phase]=bode(G,w);
>> mag=mag(1,:);
>> phase=phase(1,:);
>>
>> subplot(211),loglog(w,mag),ylabel('Magnitude'),grid
>> subplot(212),semilogx(w,phase),ylabel('Phase, deg'),grid
>> xlabel('Frequency (rad/time)')

>> %Calculation of gain and phase margins
>> %Option 1
Note: We are still doing calculations based on G=tf(1,[1 0.4 1]). This is a simple second order system that will not become unstable. Yes, the MATLAB result is meaningless. We do it anyway just to underscore the point that MATLAB won't even give you a warning if you are doing calculations that are meaningless.

>> [Gm,Pm, Wcg,Wcp]=margin(mag,phase,w)

Gm =

  652.8328    %<--meaningless results, and unit is absolute, not dB

Pm =

   33.1675

Wcg =

   25.4255

Wcp =

    1.3573

>> %Option 2
>> [Gm,Pm, Wcg,Wcp]=margin(G)

Gm =

   Inf    %<--This is a better telltale sign that the system is always stable

Pm =

   32.8599

Wcg =

   NaN

Wcp =

    1.3565

>> %Option 3
>> margin(G)

>> %Here, the gain margin results are in the plot, just underneath the "Bode Diagrams" title.
>> %If the result is valid, the value of Gm in the plot is in dB, which is NOT the case using
>> %options 1 and 2.
 

>> %To handle systems with dead time:
>> G=tf(1,[1 0.4 1]);
>> freq=logspace(-1,1);
>> [mag,phase]=bode(G,freq);
>> mag=mag(1,:);
>> phase=phase(1,:);
>> tdead=0.2;
>> phase = phase - ((180/pi)*tdead*freq); % phase is in degrees
>>
>> subplot(211), loglog(freq,mag)
>> ylabel('Magnitude'),title('Bode Plot')
>> grid
>> subplot(212), semilogx(freq,phase)
>> ylabel('Phase (degree)'),xlabel('Frequency')
>> grid

>> %With the "mag" and the "new" phase lag, we can use Option 1 as explained earlier
>> %to do gain margin calculations. Yes, we could have skipped the plotting steps.
>> [Gm,Pm,Wcg,Wcp]=margin(mag,phase,freq)

Gm =

    2.0432    %<--The unit is absolute, not dB

Pm =

   17.6017    %<--The phase margin is a bit small, but this system is stable

Wcg =

    1.7112

Wcp =

    1.3573

>> %Optional reading based on Web Supplement-- redundant calculations to see that
>> %                    there is no magic behind those MATLAB functions
>>
>> G=tf(1,[1 0.4 1]);
>> w=logspace(-1,1,200);   %Deviation from text. Use 200 points for smoother curve
>> gjw=freqresp(G,w); % does the s=jw calculation for each w
>> [re,im]=nyquist(G,w);
>>
>> plot(re(1,:),im(1,:))
>> hold
Current plot held
>> R=real(gjw(1,:)); % Real and imaginary parts of gjw are identical to the nyquist() calculation
>> I=imag(gjw(1,:));
>> plot(R,I,'x')
>> hold off

>> %Now we show that the magnitude and phase angle of gjw are identical to the bode() calculation
>> [mag,phase]=bode(G,w);
>> subplot(211), loglog(w,mag(1,:)), ylabel('Magnitude')
>> subplot(212), semilogx(w,phase(1,:)), ylabel('Phase')
>> M=abs(gjw(1,:));
>> P=angle(gjw(1,:));
>> P=P*180/pi;
>>
>> subplot(211), hold on, loglog(w,M,'x')
>> subplot(212), hold on, semilogx(w,P,'x')
>> hold off

>> %Finally, we show that nichols() is just a semilog plot of the magnitude versus phase angle
>> plot(phase(1,:),20*log10(mag(1,:)),'x')
>> title('Nichols Plot'),ylabel('Mag'),xlabel('Phase')
>> hold
Current plot held
>> nichols(G)