Untitled

Housekeeping

close all; clear; clc; s = tf('s');

The first step in our analysis is to create a sample system based on chosen poles & zeros to examine the effects of the Root Locus Plot.

Below we choose our plant's poles, zeros, and gain before building our closed-loop system based on a unity-gain controller & feedback.

gain = 1;
 
plant = gain;
for z = [-1+0j]
plant = plant*(s-z);
end
for p = [-2+4j, -2-4j]
plant = plant/(s-p);
end
 
Plant = tf(plant);
Controller = 1;
System = feedback(Controller*plant,1);

Typically one of the first tests we perform on a plant or control system is it's response to a step input, and that gives us a pretty good idea of the transient response, but in many cases the input to a system may not be a static value and may vary with time. A good example would be a self-driving car. The input to such a system would never be as simple as "go 30 mph", but would include speeding up/down, turning, obstacle avoidance, etc. A good way to examine such a system would be to test a signal that changes constantly.

Luckily, sinusoidal inputs provide a convenient template to test such an event because the signal is constantly accelerating. By substituting a pure sine wave into a transfer function, we will arrive at a ratio of two complex numbers, the numerator and the denominator.

For example, our Plant would simplify to:

Because an LTI system can scale & shift a sinusoid, but not change its frequency, the steaty-state response to an input

will be

The magnitude the response will be scaled by is:

The phase-shift the response will be scaled by is:

Examining our given plant, we can see that the magnitude and phase will simplify to the following:

Magnitude: and Phase:

With these equations in hand, let's test our system's response to a range of frequencies!

w = [];
mag = [];
phase = [];
 
for i = 0:0.01:100
w = [w i];
mag = [mag sqrt(1+i*i)/sqrt(400-24*i*i+i*i*i*i)];
phase = [phase atan(i/1)-atan(4*i/(20-i*i))];
end
 
 
f = figure();
f.Position = [0 0 1300 800];
sgtitle('Step Response Based on Pole');
 
subplot(2,1,1)
plot(w,mag,'DisplayName','magnitude');
xlabel('frequency (rad/s)');
ylabel('magnitude');
legend('FontSize',15, 'Location', 'northeast')
 
subplot(2,1,2)
plot(w, phase,'DisplayName','phase');
xlabel('frequency (rad/s)');
ylabel('phase (rad)');
legend('FontSize',15, 'Location', 'northeast')

At it's most basic, The "Bode Plot" is nothing more than this pair of plots!

Unfortunately, while the magnitude plot is reasonably easy to interperet, with a clear peak at a given frequency, the phase plot is more difficult ot interperet with most of the "action" happehing at lower frequencies, as well as a steep "jump" in phase at a given frequency.

Typically, a few modifications are made to make the plots more intuitive:

- Rather than plotting directly, the frequency is plotted logarithmically, and the magnitude is expressed in decibels.

- π is subtracted from frequencies over to better line up with the plot; essentially expressing -270° as 90° for a smoother plot.

- Phase shifts are converted to degrees. In this case, we also plot vertical lines showing the frequencies of maximum gain and phase.

mag_adj = 20*log10(mag);
 
phase_adj = phase;
phase_adj(phase_adj>pi/2) = phase_adj(phase_adj>pi/2)- pi;
phase_adj = phase_adj*180/pi;
 
f = figure();
f.Position = [0 0 1300 800];
sgtitle('Step Response Based on Pole');
 
subplot(2,1,1)
plot(w, mag_adj,'DisplayName','magnitude');
xlabel('frequency (rad/s)');
ylabel('magnitude (dB)');
legend('FontSize',15, 'Location', 'northeast')
set(gca, 'XScale', 'log')
 
subplot(2,1,2)
plot(w, phase_adj,'DisplayName','phase');
xlabel('frequency (rad/s)');
ylabel('phase (rad)');
legend('FontSize',15, 'Location', 'northeast')
set(gca, 'XScale', 'log')

Now let's compare our manual plot to the built-in Bode Plot:

f = figure();
f.Position = [0 0 1300 800];
bode(Plant)

As we can see, the Bode Plot is nothing more that the magnitude and phase-shift of the system response as the frequency of the input is varied.

Compared to the root locus, this plot (plots-plural?) are more suited to frequency-based design rather than transient-response-based design. If we wanted our system to respond most arressively to a particular frequency for example, applying a control system with the intent of shifting the magnitude peak left or right would be much more straightforward than doing the same in the s-plane.

Another excellent use for the Bode Plot[s] are in observing, or engineering around, gain and phase margin.

Firstly, there are two so-called "crossover frequencies", the gain-crossover-frequency and the phase-crossover-frequency. The gain-crossover-frequency is the frequency where the gain crosses the 0dB line, and the phase-crossover-frequency is the frequency where the phase crosses the -180° (or -180 + n*360) line.

BodeMargins.png

The Gain Margin is defined as the vertical distance (gain) we can shift the magnitude plot before the system becomes unstable.

The Phase Margin, is defined as the vertical distance (phase) we can shift the phase plot before the system becomes unstable.

NOTE: GAIN margin is defined by the PHASE crossover frequency.

NOTE: PHASE margin is defined by the GAIN crossover frequency.

In our chosen example, the gain and phase margins are considered infinite because the crossover frequencies do not occur.

For another example, consider the transfer function:

gain2 = 5;
 
plant2 = gain2;
for z = []
plant2 = plant2*(s-z);
end
for p = [-2, -2, 0]
plant2 = plant2/(s-p);
end
 
Plant2 = tf(plant2, display_format='zpk');
f = figure();
f.Position = [0 0 1300 800];
bode(Plant2)

By observation, we can see that the magnitude plot crosses the 1dB line at approximately 1 rad/s and the phase plot crosses the -180° line at approximately 2 rad/s. We can confirm this by using MATLAB's built-in function.

[gmargin2, pmargin2, gcrossover2, pcrossover2] = margin(Plant2);
gmargin2_reg = 10^(gmargin2/20);
pmargin2_reg = pmargin2*pi/180;
 
fprintf(['The Gain Margin is: %f dB (%f) at a frequency of %f rad/s\n' ...
'The Phase Margin is: %f deg (%f rad) at a frequency of %f rad/s'], ...
gmargin2, gmargin2_reg, pcrossover2,pmargin2, pmargin2_reg, gcrossover2);
The Gain Margin is: 3.200000 dB (1.445440) at a frequency of 0.999968 rad/s The Phase Margin is: 36.871359 deg (0.643527 rad) at a frequency of 2.000000 rad/s

As an illustration, let's show the gain and phase margin's graphically.

First we solve for when yielding

The equations for gain and phase margin then give us:

w = [];
mag2 = [];
phase2 = [];
 
for i = 0:0.01:100
w = [w i];
mag2 = [mag2 5/sqrt(i*i*i*i*i*i+8*i*i*i*i+16*i*i)];
phase2 = [phase2 -atan((i*i-4)/(4*i))];
end
 
mag2_adj = 20*log10(mag2);
phase2_adj = phase2;
phase2_adj(phase2_adj>pi/2) = phase2_adj(phase2_adj>pi/2)- pi;
phase2_adj = phase2_adj*180/pi;
 
mag2_adj_limit = mag2_adj + gmargin2;
phase2_adj_limit = phase2_adj + pmargin2;
 
f = figure();
f.Position = [0 0 1300 800];
sgtitle('Step Response Based on Pole');
 
subplot(2,1,1)
plot(w, mag2_adj,'DisplayName','magnitude'); hold on;
plot(w, mag2_adj_limit,'DisplayName','magnitude limit'); hold off;
xlabel('frequency (rad/s)');
ylabel('magnitude (dB)');
legend('FontSize',15, 'Location', 'northeast')
set(gca, 'XScale', 'log')
 
subplot(2,1,2)
plot(w, phase2_adj,'DisplayName','phase'); hold on;
plot(w, phase2_adj_limit,'DisplayName','phase limit'); hold off;
xlabel('frequency (rad/s)');
ylabel('phase (rad)');
legend('FontSize',15, 'Location', 'northeast')
set(gca, 'XScale', 'log')

As we can see, the system has approximately 3.2dB (or 1.45) gain margin, and 36.87° (or 0.64 rad) phase margin. This means the maximum additional gain we could add to the system would be 1.45, and the maximum additional phase we could add to the system would be 36.87°, before the system becomes unstable.