Instability of identified model in time domain with the step function

I'm trying to fit the frequency response of a system, with a high order model. I need good precision in a wide frequency range. I'm getting good results with high order models (using ARX or tfest) in the frequency domain, but the time response of the system, using the step function is unstable. My physical system is stable.

If I reduce the order of the fit, the step response became stable but the frequency fit is not good enough for what I have to do.

Does it comes from a limit of the step response ? Is my frequency fit wrong ?

t1 = linspace(0,length(out_sweep),length(out_sweep));

in_sweep = cell2mat(struct2cell(in_sweep));
out_sweep = cell2mat(struct2cell(out_sweep));

response_sweep=out_sweep*2e-5./in_sweep;
f_ax = linspace(1,fs,length(in_sweep)+1).';
f_ax(end)=[];

% Get a frd model from the raw frequency data 
Gfrd = idfrd(permute(response_sweep,[2 3 1]),f_ax,0,'FrequencyUnit','Hz');

% First fit of the transfer function with a 32 poles 32 zeros model
Gfit = tfest(Gfrd,32,32);
bode(Gfrd, Gfit)
xlim([20,20e3])
legend('Measured','Estimated')

w = Gfrd.Frequency; %#ok<NASGU>

% Crop the frequency response 
GfrdProcessed = Gfrd;
GfrdProcessed = fselect(GfrdProcessed,20,20e3);

% create a weight vector 
w = GfrdProcessed.Frequency;
f = squeeze(GfrdProcessed.ResponseData);
idx = w<40;
f(idx) = movmean(f(idx),3);

% Low weight for high frequencies
Weight = ones(size(f));
idx = w>100;
Weight(idx) = Weight(idx)/10;

% High weight for low frequencies
idx = w<1e3;
Weight(idx) = Weight(idx)*30;

% Use the weighting option for the tfest identification
tfestOpt = tfestOptions('WeightingFilter',Weight);
% GfitNew is the sys estimation of my frequency data
GfitNew = tfest(GfrdProcessed,10,tfestOpt);

% Convert the sys to plotable data 
[mag_gfrd,pha_gfrd,wout]=bode(Gfrd,2*pi*f_ax); %#ok<ASGLU>
mag_gfrd=squeeze(mag_gfrd);
pha_gfrd=squeeze(pha_gfrd);

[mag_gfitnew,pha_gfitnew,wout]=bode(GfitNew,2*pi*f_ax);
mag_gfitnew=squeeze(mag_gfitnew);
pha_gfitnew=squeeze(pha_gfitnew);

Answer 1

[Sulaymon Eshkabilov]

Can you post your data out_sweep to play with your code and come up with a solution to your problem?

[Romain Liechti]

Hi, thanks for answering. I should've edited or answered my own post since I found the answer in the documention of the options of the tfest function.

https://fr.mathworks.com/help/ident/ref/tfestoptions.html

The option 'EnforceStability' just needed to be set to true.

It doesn't explain how a model of a physical system which is stable can be unstable and fitting perfectly the frequency response but it gives a stable model.

% Use the weighting option for the tfest identification
% Stability of the identified function is forced with 'EnforceStability', true
tfestOpt = tfestOptions('WeightingFilter',Weight,'EnforceStability',true);