Matlab Margin Display From a Bode Plot Code

w = logspace(-1,2,100);
for k = 1:100
s = 1i * w(k);
G(k) = (106.46 / (s+7.69))*(1/(s + 23.09)*((0.1446*s + 11.568)/s));

end

figure(1)
subplot(2,1,1)
semilogx(w,20*log10(abs(G)));

figure(2)
subplot(2,1,2)
semilogx(w,angle(G)*180/pi)

I have a code for bode plot although I also need to find margins from plot or check them by using 'all margin' command. However I couldn't be able to do that since I didn't wrote the G(k) as tf[ ]. How can I find the margin values (gain margin, phase margin, crossover freq.) ?

Answer 1

[Jon]

Assuming you have the Control System Toolbox, as you say there is already a function for computing the margins. If so, why not just define your transfer function and use that command like this:

% (106.46 / (s+7.69))*(1/(s + 23.09)*((0.1446*s + 11.568)/s));

% define transfer functions
tf1 = tf(106.46,[1 7.69]);
tf2 = tf(1,[1 23.09]);
tf3 = tf([0.1446 11.568],[1 0]);

% define overall transfer function
tf = tf1*tf2*tf3;

% compute gain and phase margins
S = allmargin(tf)

S = struct with fields:
     GainMargin: [7.2132 Inf]
    GMFrequency: [16.9884 Inf]
    PhaseMargin: 44.9582
    PMFrequency: 5.5006
    DelayMargin: 0.1427
    DMFrequency: 5.5006
         Stable: 1

Otherwise, you can just go to the fundamental definition of these quantitities and solve for it yourself. You could use fzero for this.

Answer 2

[Harsh Kumar]

Hi Yigit,

It is my understanding that you are interested in calculating the margins without using the "tf"command to convert the transfer function.

To resolve this issue, you may utilize the "zpk" function instead of "tf" function to achieve the same accuracy. OR

In case, you want to achieve the same results through manual method, use the basic definitions of margins where the accuracy depends on step size of your frequency w = logspace(-1,2,100) .

For your reference, the provided code demonstrates a simple method to compute the margins without using built-in "tf" function using basic definations of margins:

w = logspace(-1, 2, 100);
G = @(s) (106.46 / (s+7.69))*(1/(s + 23.09)*((0.1446*s + 11.568)/s));

magnitude_dB = zeros(size(w));
phase_deg = zeros(size(w));
gain_margin=0;
phase_margin=0;
Gain_crossover_freq=0;
Phase_crossover_freq=0;


for k = 1:length(w)
s = 1i * w(k);
Gk = G(s);

% Calculate magnitude in dB and phase in degrees
magnitude_dB(k) = 20 * log10(abs(Gk));
phase_deg(k) = angle(Gk) * 180 / pi;

%to get the most nearest value
if magnitude_dB(k)*magnitude_dB(1)>0
Gain_crossover_freq=w(k);
phase_margin=phase_deg(k)+180;
end

%to get the most nearest value
if phase_deg(k)*phase_deg(1)>0
Phase_crossover_freq=w(k);
gain_margin=-1*magnitude_dB(k);
end

end

subplot(2,1,1)
semilogx(w, magnitude_dB)
title('Magnitude Bode Plot')
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')

subplot(2,1,2)
semilogx(w, phase_deg)
title('Phase Bode Plot')
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')

% Plotting margins
Gain_crossover_freq
gain_margin
Phase_crossover_freq
phase_margin