Calculation of Phase margin

Please, I would like to ask a control system toolbox expert a question. Solving an exercise in which we are asked to calculate the value of the gain k which causes a transfer function G (s) = k / [s (s ^ 2 + s + 4)] to have a phase margin of 50°. The analytical calculation by hand gave me the following answer:

gain plot crossing frequency : wcf = 1.491 rad/s

The desired gain k = 3.4585

fase plot crossing frequency : wcg = 2 rad/s

Gain margin : GM=1.2634 dB

Instead the Matlab script to calculate the phase margin with this data k, gives an answer that is not the answer I expected and I can't understand why Matlab makes a different choice.

__________________________________________

% Problem_B_7_26_Modern_Control_Engineering_Katsuhiko_Ogata 
clear all
clc
fprintf(' \n')

fprintf('--------------------- \n')

--------------------- 

w=1.4910;
k=sqrt((w^4)+(w^2)*((4-w^2)^2));
fprintf('Desired gain : k=%g° \n',k)

Desired gain : k=3.45852° 

fprintf('--------------------- \n')

--------------------- 

fprintf('Open loop transfer function \n')

Open loop transfer function 

s=tf('s');
G=zpk((k)/(s*(s^2+s+4)))

G =
 
      3.4585
  ---------------
  s (s^2 + s + 4)
 
Continuous-time zero/pole/gain model.

fprintf('--------------------- \n')

--------------------- 

fprintf('Figure 1 : Bode diagram \n')

Figure 1 : Bode diagram 

figure(1)
margin(G)
grid on

fprintf('--------------------- \n')

--------------------- 

[GM,PM,wcg,wcf]=margin(G);
fprintf('Phase margin: \t\t\t\t PM=%g° \n',PM)

Phase margin: 				 PM=31.2163° 

fprintf('Gain margine: \t\t\t\t GM=%g dB \n',db(GM))

Gain margine: 				 GM=1.2634 dB 

fprintf('Gain crossing frequency : \t %scf=%g rad/s \n',char(969),wcf)

Gain crossing frequency : 	 ωcf=1.71982 rad/s 

fprintf('Phase crossing frequency : \t %scg=%g rad/s \n',char(969),wcg)

Phase crossing frequency : 	 ωcg=2 rad/s 

fprintf('--------------------- \n')

--------------------- 

Answer 1

[Paul]

Hi T, it's a tricky problem becase of the shape of the gain curve around the phase cross over frequency.

w=1.4910;
k=sqrt((w^4)+(w^2)*((4-w^2)^2));
s = tf('s');
G=zpk((k)/(s*(s^2+s+4)));

Let's check the stability margins of G

allmargin(G)

ans = struct with fields:
     GainMargin: [1.1566 1.1566]
    GMFrequency: [2.0000 2.0000]
    PhaseMargin: [58.2662 50.0064 31.2163]
    PMFrequency: [1.3487 1.4909 1.7198]
    DelayMargin: [0.7540 0.5854 0.3168]
    DMFrequency: [1.3487 1.4909 1.7198]
         Stable: 1

Sure enough, we see a 50 deg phase margin at the expected frequency of 1.491 (I don't know why it shows two gain margins at the same frequency, that seems odd to say the least). However, we see that there are two additional gain crossovers and so two additional phase margins. I think that margin() returns only the minimum of the gain/phase margins

[GM,PM,wcg,wcf]=margin(G);
PM

PM = 31.2163

The Bode plot reveals reveals the three gain crossovers at 1.3487, 1.491, and 1.7198 rad/sec

bode(G,linspace(1.2,1.8,1000)),grid

[T Hafid]

My version of Matlab 2015a gives me a single value for the gain margin and not two identical values as in your case.

I must say that I would have preferred Matlab to have immediately provided me with the 3 different values of the phase margin, rather than a single value chosen by Matlab on the basis of its own logic (minimum value as you say or other reason).

In fact, in that case I would have immediately seen that the solution calculated by me, calculated by hand, was among the solutions and would have been obliged me to investigate the correctness of the other 2 solutions. Instead, having received only the solution with the minimum phase magin, forced me only to redo the calculations by hand more than three times to no avail.

Thanks a lot Mr. Paul for your help.

[T Hafid]

Thanks Mr Paul for your help