Why do I not get a perfect fit while estimating a frequency response data with a higher natural frequency using the System Identification Toolbox 7.3 (R2009a)?

Question

MathWorks Support Team 2009-6-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/95507-why-do-i-not-get-a-perfect-fit-while-estimating-a-frequency-response-data-with-a-higher-natural-freq

I was trying a identify transfer function from a frequency response dataset ( magnitude- phase-response) using IDPROC method from the System Identification Toolbox.

So I used an example of a simple second order underdamped transfer function with natural frequency of 4 rad/sec and Zeta of 0.25. This transfer function was identified perfectly with a fit of 100%.

However, when I changed the natural frequency to 2.5133e+004 rad/sec (corresponding to a frequency of 4KHz), I did not get a good fit. Here is my reproduction code:

f = 4000;
wn = 2*pi*f; %rad/sec
zeta = 0.25;
tf('s');
num = [1];
sec= 2*wn*zeta;
third = wn^2;
den = [1 sec third];
csys = tf(num,den); % continuous transfer function
Ts = 0.00001; % Ts = 0.00001 sec;
dsys = c2d(csys,Ts); 
[P,PHA,W] = bode(dsys); % bode of discrete sys, W is in rad/sec, P is in Db, PHA is in degrees
zfr = P.*exp(i*PHA*pi/180);
gfr = idfrd(zfr,W,Ts);
figure(1);
bode(gfr), legend('gfr')
mproc = pem(gfr,'p2u') % 2nd-order, continuous-time model with underdamped poles
m = pem(gfr, mproc)

请先登录，再回答此问题。

Answer 1

MathWorks Support Team 2009-6-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/95507-why-do-i-not-get-a-perfect-fit-while-estimating-a-frequency-response-data-with-a-higher-natural-freq#answer_104860

This is an expected behavior while using System Identification Toolbox 7.3 (R2009a). The reason for this behavior is as follows:

IDPROC imposes a lower limit on the value of Tw = 1/wn = 1e-3 as a default. For 4 kHz, wn is 2.5133e4 which leads to Tw value of the order of 1e-5.

Hence one is unlikely to get anything good using default settings. One must configure the parameter value bounds before estimating, as in:

f = 4000;
wn = 2*pi*f; %rad/sec
zeta = 0.25;
tf('s');
num = [1];
sec= 2*wn*zeta;
third = wn^2;
den = [1 sec third];
csys = tf(num,den); % continuous transfer function
Ts = 0.00001; % Ts = 0.00001 sec;
dsys = c2d(csys,Ts); 
[P,PHA,W] = bode(dsys); % bode of discrete sys, W is in rad/sec, P is in Db, PHA is in degrees
zfr = P.*exp(i*PHA*pi/180);
gfr = idfrd(zfr,W,Ts);
figure(1);
bode(gfr), legend('gfr')
mproc = idproc('p2u');
mproc.Tw.min = 1e-6; % lower the min bound to be smaller than expected value
m = pem(gfr, mproc);

In general, it is advisable to follow this multi-step process rather than directly doing "pem(data, 'p2u')". One should carefully inspect the min and max values of each parameter and preferably also set the initial guess values in the 'value' field of parameter (wherever possible).