Complex equation in MATLAB (control engineering)

I want to determine at which dead time the following open loop with the transfer function (F_O) becomes unstable. Unfortunately, I cannot do this with MATLAB.

I have now determined the solution using GeoGebra. However, I do not understand what I did wrong in MATLAB. What did I type in wrong?

You can see my MATLAB code in the Attachment.

MATLAB-Code:

clc;

syms s w T real


F_O = (exp(-s*T))/(s*(1+s))

F_O_w = subs(F_O,s,(j*w))
R = real(F_O_w)
I = imag(F_O_w)

eqns = [R == -1 , I == 0]

[w_D T_t] = vpasolve(eqns,[w T]) 

Command Window:

F_O =
 
exp(-T*s)/(s*(s + 1))
 
 
F_O_w =
 
-(exp(-T*w*1i)*1i)/(w*(1 + w*1i))
 
 
R =
 
- cos(T*w)/(w^2 + 1) - sin(T*w)/(w*(w^2 + 1))
 
 
I =
 
sin(T*w)/(w^2 + 1) - cos(T*w)/(w*(w^2 + 1))
 
 
eqns =
 
[- cos(T*w)/(w^2 + 1) - sin(T*w)/(w*(w^2 + 1)) == -1, sin(T*w)/(w^2 + 1) - cos(T*w)/(w*(w^2 + 1)) == 0]
 
 
w_D =
 
-0.78615137775742328606955858584296
 
 
T_t =
 
12109.538400133501929986843400658

GeoGebra-Solution:

Realteil ... Real part

Imaginärteil ... Imaginary part

Thanks for any answer!

Answer 1

[Torsten]

Seems there are several solutions.

Use

[w_D T_t] = vpasolve(eqns,[w T],[1 1]) 

instead of

[w_D T_t] = vpasolve(eqns,[w T])

in order to get the solution from GeoGebra.

[Michael Weingartshofer]

Thank you very much!

Answer 2

[Paul]

Hi Michael,

vpasolve found a solution, even if it's not the one you want.

syms s
syms w T real
F_O = (exp(-s*T))/(s*(1+s));
F_O_w = subs(F_O,s,(j*w));
R = real(F_O_w)

R = 

I = imag(F_O_w)

I = 

eqns = [R == -1 , I == 0];
[w_D T_t] = vpasolve(eqns,[w T])

w_D = 

T_t = 

12109.538400133501929986843400658

subs(eqns,[w T],[w_D T_t])

ans = 

Use an additional input to vpasolve to specify a positive interval

[w_D T_t] = vpasolve(eqns,[w T],[0 inf;0 inf])

w_D = 

0.78615137775742328606955858584296

T_t = 

3278.008034774571436567137274222

Now, we have the correct frequency. Looking at the equations for R and I, we see that if w*T is a solution, then so is w*T + 2*k*pi. So to find the T_t we really want, we need to do mod arithmetic

T_t = vpa(mod(T_t*w_D,2*pi)/w_D)

T_t = 

1.150614143656049868422783870361

If we had some idea of the expected answer, then we can use vpasolve as

[w_D T_t] = vpasolve(eqns,[w T],1)

w_D = 

0.78615137775742328606955858584296

T_t = 

1.150614143656049868422783870361

Easier is to just use allmargin to get the delay margin and its frequency.

allmargin(tf(1,[1 1 0]))

ans = struct with fields:
     GainMargin: Inf
    GMFrequency: Inf
    PhaseMargin: 51.8273
    PMFrequency: 0.7862
    DelayMargin: 1.1506
    DMFrequency: 0.7862
         Stable: 1