System Identification - Frequency Domain

I am trying to model the transfer function for the plant of a magnetic actuator. The input is in current (A) and output is distance moved x (mm). A sine sweep is run as disturbance and the input output plot looks like this:

Since it is a sine sweep, i convert the data into frequency domain as follows:

I am trying to identify the transfer function for this frequency domain plot. I use this as an idfrd data. I do manage to find a 6th order transfer function that fits the data. The corresponding code is as follows:

f_data=idfrd(H3c,fc,t_loop,'FrequencyUnit','Hz','Name', 'Openloop Plant', ...
    'InputName', 'currrent', 'OutputName', 'xdist',...
    'InputUnit', 'A', 'OutputUnit', 'mm');
bode(f_data)

f_data.InterSample = 'zoh';
model_1=tfest(f_data,6)
bode(model_1)
hold on
bode(f_data)
figure()
compare(f_data,model_1)

From this, I get the following match:

The problem is, I can't get this back in the time domain. Or atleast, it is becoming unstable in the time domain like follows:

y=sim(model_1,Signal);
plot(t,y)

Can anybody give any suggestions as to why this is happening and what to do to make it right?

Answer 1

[Rajiv Singh]

Make an attempt with stability enforced.

opt = tfestOptions('EnforceStability', true);

model=tfest(f_data,6,opt)

Also, you don't necessarily need to generate FRD representation for frequency domain identification (conversion of time domaoin data to FRD is itself an estimation problem, subject to bias/variance considerations). You could try something like this too:

Time_Domain_Data = iddata(y,t,Ts);

Frequency_Domain_Data = fft(Time_Domain_Data);

model=tfest(Frequency_Domain_Data,6,opt)

[Neha Chohan]

Yes, switching on enforce stability in options does the job. Thanks!

I will also take into consideration your second comment. :)

[Neha Chohan]

The phase angle is not matching yet. Any suggestions for that?

Thanks in advance.

[Rajiv Singh]

The algorithm fits the complex frequency response, not magnitude and phase separately. I can't think of a way of improving phase only while keeping the magnitude curve unchanged, with the exception of rotations induced by delays. If there is a delay present from input to output, try finding out its value (e.g., using impulse or step test, or using delayest) and adding that value to the estimated model. But I doubt that is the case here.. you are probably forcing a stable model on the data from a (seemingly) unstable process.