Plots & Analysis
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.
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);
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')