tfest function does not recognize nearby poles and zeros in a 2dof system

I am using tfestimate and tfest to experimentally derive the frequency response function of a 2 dof system with 1 input and 2 output.

The problem is that tfest seems to cancel or does not recognize a complex pole pair (resonance) near a complex zero pair (antiresonance) around 6 Hz in the tf from input 1 to output 1 and consequently simplifies the frequency response of the system (see image, correct plot in blue vs simplified plot in red):

I tried different types of windowing in tfestimate but the only one that allows me to obtain the correct frequency response after tfest is rectwin without overlap: this works only with mathematical model I/O data, the blue plot shown above, but with the real model data the problem persists. Other windows like hann (which would be more suitable) also leads to a pole-zero simplification after tfest, with mathematical model I/O data too.

Increasing the number of zeros-poles from 2-4 to 6-8 tfest gives me the correct transfer functions, but I absolutely need the model with 2 zeros and 4 poles for a correct dynamic of the system to control by using the identification of its FRF.

To load I/O data see attachment .mat file

The code is the following:

% Upload I/O data
load 'test_identificazione_chirp_tauGain_0,1.mat'

% Data from real model
th1abs = theta1_rad_MEAS.signals.values;   
th2rel = theta2_rad_MEAS.signals.values;
th2abs = th1abs+th2rel;
tauIn = tauInput.signals.values;

nfs = 4096*4;   % Number of samples for FFT
Fs = 1000;   % Sampling frequency

% Transfer functions estimate from input 1 (tauIn) to outputs 1 and 2 (th1abs and th2abs)
[H_11, f_H_11] = tfestimate(squeeze(tauIn), squeeze(th1abs), rectwin(floor(size(squeeze(tauIn), 1))), 0, nfs*4, Fs);
[H_21, f_H_21] = tfestimate(squeeze(tauIn), squeeze(th2abs), rectwin(floor(size(squeeze(tauIn), 1))), 0, nfs*4, Fs);
% Try also hann window with overlap to get a better curve: hann(floor(size(squeeze(tauIn), 1) / 15)), 3000

% Estimated tf smoothing, by movmean with [x,x] values
H_11 = movmean(H_11,[2,2]);
H_21 = movmean(H_21,[2,2]);

% Estimated tf plot
figure(1)
subplot(1,2,1);
semilogx(f_H_11, mag2db(abs(H_11)), 'b', 'LineWidth', 0.7);
hold on;
subplot(1,2,2);
semilogx(f_H_21, mag2db(abs(H_21)), 'b', 'LineWidth', 0.7);
hold on;

% Frequency range to keep only relevant part of tf by tfestimate
fmin = 0.1;   % min freq (Hz)
fmax = 15;   % max freq (Hz)
idx_11 = (f_H_11 >= fmin) & (f_H_11 <= fmax);   % frequency index in the selected range
idx_21 = (f_H_21 >= fmin) & (f_H_21 <= fmax);   % frequency index in the selected range

H_11_filtered = H_11(idx_11);
f_11_filtered = f_H_11(idx_11);   % corresponding frequencies
H_21_filtered = H_21(idx_21);
f_21_filtered = f_H_21(idx_21);   % corresponding frequencies

% rad/s conversion to use in `idfrd`
f_11_rad = f_11_filtered * 2 * pi;
f_21_rad = f_21_filtered * 2 * pi;

% Filtered tf plot
subplot(1,2,1);
semilogx(f_11_filtered, mag2db(abs(H_11_filtered)), 'r', 'LineWidth', 0.7);
subplot(1,2,2);
semilogx(f_21_filtered, mag2db(abs(H_21_filtered)), 'r', 'LineWidth', 0.7);

% Creating an identified model based on frequency data
H_11_idfrd = idfrd(H_11_filtered, f_11_rad, 'Ts', 0);   % continuous time
H_21_idfrd = idfrd(H_21_filtered, f_21_rad, 'Ts', 0);   % continuous time

% tf estimation in "tf" variable type
np = 4;   % system pole number
nz = 2;   % system zero number

sys_est_11 = tfest(H_11_idfrd, np, nz);
sys_est_21 = tfest(H_21_idfrd, np, nz);

sys_est = [sys_est_11; sys_est_21];   % [2x1] dimension

% Estimated FRF plot
figure(2)
P = bodeoptions; 
P.FreqUnits = 'Hz';   % Frequency axys in Hz
bodeplot(sys_est, P);
grid on;
hold on;

% Print estimated tf in the Command Window   
sys_est

sys_est =
 
  From input to output...
               90.94 s^2 + 194.4 s + 5826
   1:  ------------------------------------------
       s^4 + 2.259 s^3 + 119.8 s^2 + 135 s + 3475
 
                -196 s^2 - 998.8 s + 1.48e05
   2:  ----------------------------------------------
       s^4 + 3.027 s^3 + 1457 s^2 + 1488 s + 8.052e04
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: [4;4]   Number of zeros: [2;2]
   Number of free coefficients: 14
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                 
Created by concatenation of originally estimated models.
 

Any help? Thank you in advance!

[Mathieu NOE]

hello

In my eyes tfest is not really a full MIMO identification tools . It works fine for SISO or MISO models mostly (only ?) - have to read the doc again,, it's been a long since last time I used it.

my second comment would be that even if you have a perfect identifcation of your close by poles and zeores , it's not a very robust situation. You may be happy to have a nice model , but reality is that polse and zeoes may have slight frequency shift due to whatever operationnal parameters and therefore your controller design is not robust against those uncertainties. Do you plan to use robuts controllers (H2, Hinf ?)

do the control bandwith really requires to have a super accurate model in this frequency range ?

and last , did you try eventually ssest