# Instability of identified model in time domain with the step function

6 views (last 30 days)
Romain Liechti on 15 Mar 2019
Edited: Romain Liechti on 20 May 2019
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);

Sulaymon Eshkabilov on 19 May 2019
Can you post your data out_sweep to play with your code and come up with a solution to your problem?
Romain Liechti on 20 May 2019
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.
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);

R2018b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!