7.7 MATLAB Control Toolbox: Bode and Nyquist Functions

The MATLAB Control Toolbox has special functions to generate Bode and Nyquist plots, and to calculate gain and phase margins. Here we use Example 7.3 to illustrate the use of the MATLAB bode, nyquist, and imargin functions that are available in the Control Toolbox. For information, simply enter help bode, and so forth, in the MATLAB command window.

First, define the process transfer function (g_p).

» vdvtf = tf(0.5848*[-0.3549 1],[0.1858 0.8627 1])

Transfer function:
   -0.2075 s + 0.5848
-------------------------
0.1858 s^2 + 0.8627 s + 1

Next, generate the open-loop transfer function (g_cg_p) for a controller gain (k_c) of 2.5.

» kc = 2.5;

» gcgp = kc*vdvtf


Transfer function:
    -0.5189 s + 1.462
-------------------------
0.1858 s^2 + 0.8627 s + 1

Now, use the bode function to generate magnitude and phase angles as a function of frequency. Also, use imargin to perform the gain margin (Gm) and phase margin (Pm) calculations at their respective frequencies (Wcg and Wcp).

»[mag,phase,w] = bode(gcgp);

[Gm,Pm,Wcg,Wcp] = imargin(squeeze(mag),squeeze(phase),w)

Gm =
    1.6680

Pm =
   59.0905

Wcg =
    4.3064

Wcp =
    2.0956

Figure 7-14 is generated using the following sequence of commands.

subplot(2,1,1),loglog(w,squeeze(mag),[0.1 Wcg Wcg],
   [1/Gm 1/Gm min(mag)],'--')
subplot(2,1,2),semilogx(w,squeeze(phase),[min(w) Wcg Wcg],
   [-180 -180 0],'--')

Figure 7-15 is generated using the following sequence of commands.

subplot(2,1,1),loglog(w,squeeze(mag),[0.1 Wcp Wcp],
   [1 1 min(mag)],'--')
subplot(2,1,2),semilogx(w,squeeze(phase),[0.1 Wcp Wcp],
   [-180+Pm -180+Pm 0],'--')

The Nyquist plot, Figure 7-16, is generated using the following sequence of steps. First, define an appropriate frequency range.

w1 = logspace(-2,2,300);

Then use the nyquist function to generate real and imaginary vectors.

[regcgp,imgcgp] = nyquist(gcgp,w1);

Generate circles for convenience when interpreting the Nyquist plot.

circlex = cos(0:0.2:2*pi);
circley = sin(0:0.2:2*pi);
phaseptx = cos((-180+Pm)*pi/180);
phasepty = sin((-180+Pm)*pi/180);
plot(squeeze(regcgp),squeeze(imgcgp),circlex,circley,'--',
   [-1/Gm],[0],'o', phaseptx,phasepty,'*')
axis('square')
xlabel('real')
ylabel('imaginary')