T2 = readtable('E7i_CSV.csv');   
                                                                       % Fill With Actual Sampling Frequency
FHz = T2.F;
Ts = 1/(2*(max(FHz)+10000))
for k = 1:size(T2,1)-1
    if FHz(k+1) == FHz(k)
        FHz(k+1) = FHz(k+1)+0.5;                                                % 'Brute Force' Interpolation 
    end
end
Mag = T2.G;
PhDeg = T2.P;
Response = Mag.*exp(1j*deg2rad(PhDeg));                                         % Complex Vector
sysfr = idfrd(Response, FHz, Ts, 'FrequencyUnit','Hz')
%bode(sysfr)
%tfsys = tfest(sysfr,4,3) This is much better
tfsys = tfest(sysfr,4,3) %82.69%
%tfsys = tfest(sysfr,6,3)
figure
compare(sysfr, tfsys)
is_stable = isstable(tfsys);
disp(is_stable);

I want to convert this transfer function to statespace equations, which will be used for Model Predictive Control Design. Simple tf2ss command give me TF but it doesn't look very accrurate.

[A, B, C, D] = tf2ss(num,den)

So, I was reading a paper which i will follow, they have mentioned that they had FRF (Frequency Response Function) data and they find state space model/equations using Prediction Error Method (PEM). But I am unable to understand how I can used Bode Plot data and convert it state space using PEM. I will appreciate your help.

Secondy I verified the paper StateSpace Model by first converting it to Transfer Function and then using tf2ss command. which shows inaccurate results similiar to my results. I am looking for your help. Following is the verification code for paper i am following

%Rana Paper State Space Equations Coversion to
                            %Transfer Function for Validation Reason
% Define state-space matrices
Ax = [-49.3277, -980.7648, 618.2198;
      567.3557, -120.2273, 651.7342;
      -45.8885, -415.4194, -115.7522];
Bx = [-1.8122; 3.901370; 4.3011];
Cx = [14.9368, 16.9208, -23.6827];
Dx = 0;
% Create a state-space model
sys = ss(Ax, Bx, Cx, Dx);
% Display the state-space model
disp('State-Space Model:');
disp(sys);
% Convert the state-space model to a transfer function
tf_sys = tf(sys);
% Display the transfer function
disp('Transfer Function:');
disp(tf_sys);
                        %TF to StateSpace Rana Model
num_rana = [ -62.92 3.625*10^4 -1.065*10^8];
den_rana = [ 1 285.3 8.811*10^5 1.982*10^8];
tff_rana = tf(num_rana,den_rana);
tffss_rana = tf2ss(num_rana,den_rana);
[Arana,Brana,Crana,Drana] = tf2ss(num_rana,den_rana);