Deriving the transfer function from a data set that has magnitude,phase and omega

Hello,

I have this dataset that contains magnitude, phase and omega of a transfer function. I can bodeplot it but I need to derive the transfer function from it. Please help.

Thanks

Answer 1

[Paul]

Hi George,

I plotted the data and it looks like the data comes from a transfer function where the degree of the denominator is two more than the degree of the numerator.

First, store the data in an frd object. Then try using tfest to find the transfer function that fits the data. It should work very well.

If you have any trouble, feel free to come back, show your code, and and ask more questions.

[George Isaac]

I tried this but didn't work

gain=M;

response=gain.*exp(i*P*pi*180);

Ts=0.1;

gfr=idfrd(response,w,Ts);

sys=tfest(gfr,2);

[Paul]

The equation for response needs to be corrected use pi/180 (not pi*180). Also, it's safer to use 1i instead of i

Is the data from a continuous-time or discrete-time system? I assumed it was continuous.

If it's for a discrete-time system, there seems to be a disconnect between the sample period, Ts, and the data. For a discrete-time system, the mag and phase data is typically only provided up to the Nyquist frequency pi/Ts = 31.4 rad/sec. For frequencies higher than that, the frequency response repeats because it's periodic, but that's not the case for the data in activity.mat.

I'm not familiar with idrfd objects, which is why I suggested using an frd object. idfrd might work, but I can't say for sure.

[George Isaac]

I got a very accurate transfer function but im not sure how i can analyze the poles and zeros chosen for it!

406.5 s + 40.67

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

s^3 + 32.51 s^2 + 81 s + 610.1

[Sam Chak]

Hi @George Isaac

There are a few commands you can apply to analyze the estimated transfer function.

s = tf('s');
G = (406.5*s + 40.67)/(s^3 + 32.51*s^2 + 81*s + 610.1)

G =
 
         406.5 s + 40.67
  ------------------------------
  s^3 + 32.51 s^2 + 81 s + 610.1
 
Continuous-time transfer function.

z = zero(G)     % zero

z = -0.1000

p = pole(G)     % 3 poles

p = 
 -30.5106 + 0.0000i
  -0.9997 + 4.3585i
  -0.9997 - 4.3585i

rlocus(G)

margin(G)

step(G, 10)

stepinfo(G)

ans = struct with fields:
         RiseTime: 0.0131
    TransientTime: 4.1285
     SettlingTime: NaN
      SettlingMin: -0.9759
      SettlingMax: 2.2098
        Overshoot: 3.2149e+03
       Undershoot: 1.4639e+03
             Peak: 2.2098
         PeakTime: 0.3472

[George Isaac]

Does anyone have an explanation of how zeta (damping ratio) is 0.22. The peak value gain is 6.59 and the Kdc is 0.06667. How come zeta is 0.22 in the transfer function?

[Sam Chak]

s = tf('s');
G = (406.5*s + 40.67)/(s^3 + 32.51*s^2 + 81*s + 610.1)

G =
 
         406.5 s + 40.67
  ------------------------------
  s^3 + 32.51 s^2 + 81 s + 610.1
 
Continuous-time transfer function.

If you transform the transfer function into the zero-pole form, you will notice that the characteristic polynomial (denominator) has a pair of complex-valued poles (roots):

G2 = zpk(G)

G2 =
 
          406.5 (s+0.1)
  -----------------------------
  (s+30.51) (s^2 + 1.999s + 20)
 
Continuous-time zero/pole/gain model.

.

Compare with , you can obtain the damping ratio:

zeta = 1.999/(2*sqrt(20))   % just an approximation!

zeta = 0.2235

The DC gain is essentially .

Kdc = dcgain(G)

Kdc = 0.0667

damp(G)

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
 -1.00e+00 + 4.36e+00i     2.24e-01       4.47e+00         1.00e+00    
 -1.00e+00 - 4.36e+00i     2.24e-01       4.47e+00         1.00e+00    
 -3.05e+01                 1.00e+00       3.05e+01         3.28e-02    

[George Isaac]

Yes I understand all of this isn't the formula for zeta the one i just attached?

How is Mr in this case 0.155?

[Paul]

Hi George,

Looks like you got a great fit to the data.

I suspect that formula for zeta is only applicable for an ideal second order system. But in this case we have

sys = zpk(tf([406.5  40.67],[1 32.51 81 610.1]))

sys =
 
          406.5 (s+0.1)
  -----------------------------
  (s+30.51) (s^2 + 1.999s + 20)
 
Continuous-time zero/pole/gain model.

That zero s = -0.1 will cause the magnitude response to start to increase starting at w = 0.1 rad/sec before the pole at wn = sqrt(20) = 4.6 rad/sec starts the rolloff.

bode(sys)

The resonant peak at 4.6 rad/sec is so high because of the low frequency zero.

[George Isaac]

Hi @Paul

Is there a reason why if you go back to the formula and trace back Mr, it comes up as 0.155 which is exactly the inverse of the Mr thats apparent in the graph at 6.59?

[Paul]

The formula for Mr as cited is

sys = zpk(tf([406.5  40.67],[1 32.51 81 610.1]))

sys =
 
          406.5 (s+0.1)
  -----------------------------
  (s+30.51) (s^2 + 1.999s + 20)
 
Continuous-time zero/pole/gain model.

Kdc = freqresp(sys,0);
zeta = 0.224; % damping ratio of second order poles
Mr = Kdc/2/zeta/sqrt(1-zeta^2)

Mr = 0.1527

which I guess is the 0.155 that you're referring to.

Is the 6.59 the peak of the amplitude curve?

m = bode(sys,logspace(-1,2,1000));
maxm = max(m)

maxm = 6.5942

So I guess that is it.

1/maxm

ans = 0.1516

So close, but not exact. Then again, the results in this comment aren't based on the exact transfer function you have.

Offhand, I can't see a reason why this is anything more than a coincidence, but I haven't thought about it too much. Is there any reason to think otherwise?

[George Isaac]

Yes I’m talking about the frequency where the bode plot is peaking. It peaks at 6.59 rad/sec which is the Mr? But somehow it’s 0.15

[Paul]

I was editing my comment as you were typing in yours.

In any case, the magnitude peaks at a value of Mr = 6.5948 and the frequency at the peak is w = 4.464 rad/sec.

sys = zpk(tf([406.5  40.67],[1 32.51 81 610.1]));
w = logspace(-1,2,5000);
m = bode(sys,w);
[maxm,index] = max(m);
maxm

maxm = 6.5948

w(index)

ans = 4.4641

[George Isaac]

Mr=6.59;
Wr=4.464;
Kdc=0.06667;

zeta_vec= roots([4 0 -4 0 (Kdc/Mr)^2]);
zeta_vec;

yes but isn't the peak value of Mr=6.5948 the one that should be used in the formula.. This would come out with zeta not being 0.22

zeta_vec =

-1.0000

1.0000

-0.0051

0.0051

[Paul]

Strictly speaking, the formula only applies for a second order system. For example, with zeta = 0.2 and Kdc = 1, we would expect to have

Kdc = 1;
zeta = 0.2;
Mr = 1/2/zeta/sqrt(1-zeta^2)

Mr = 2.5516

20*log10(Mr)

ans = 8.1361

And we do

bode(tf(1,[1 2*zeta*1 1]),tf(Mr))

If the transfer function is not second order, the formula might still be a very good approximation if the additional poles and zeros are at higher frequency than that of the dominant second order poles. For example, let's add a zero at s= -10 and a pole at s = -20, which don't have much impact on the low-frequency response and don't really affect the Mr

bode(tf([1/10 1],conv([1/20 1],[1 2*zeta*1 1])),tf(Mr),logspace(-1,1,500))

In your transfer function, the dominant second order poles are at ~4.45 rad/sec, but the transfer function has a zero at s = -0.1, which is smaller than the frequency of the dominant poles. That zero has a sigificant impact on the low frequency portion of the Bode plot (and the overshoot of the step response as illuatrated above by @Sam Chak), so the formula doesn't apply.