Stability Analysis of Bode Plot

Dear Control Community,

I have found open loop frequent response of my plant and then estimated transfer function from it. When I do stability analysis of my plant bode plot, it shows GM -17dB and Phase Margin -54 degree, and shows close loop unstable.

I am very confused and worried about it how can I make this system stable. This device is currently working with PI controller in close loop but I have taken plant open loop response without PI controller.

Is it possible to make this system close loop stable by PI controller or model predictive controller?

My main goal is to design Kalman Estimator and Model Predictive Control for this plant to improve its performance.

%Experimental Data coversion to Bode Plot
T2 = readtable('EIDEL_CSV.csv');
% T2 = readtable('EIDEL2_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
% Remove duplicate entries from the frequency vector
FHz = unique(FHz);
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,5,4); %82.69%
%tfsys = tfest(sysfr,6,3)
figure
compare(sysfr, tfsys);
is_stable = isstable(tfsys);
disp(is_stable);
% Assuming you have a tfsys transfer function object named 'tf'
num = tfsys.Numerator;  % Extract numerator coefficients
den = tfsys.Denominator;  % Extract denominator coefficients
% Plot poles and zeros
figure;
pzmap(tfsys);
title('Poles and Zeros of the Estimated Transfer Function');
ss_model = ss(tfsys);
[A, B, C, D] = tf2ss(num,den);
figure
compare(sysfr, tfsys, ss(A,B,C,D), ss_model);
%compare(sysfr, tfsys, ss_model, ss_model2)
step(tfsys);
% Define the frequency range
% Define the frequency range
% frequencies = logspace(0, log10(20000), 1000); % Creates a logarithmic frequency range from 1 Hz to 20 kHz
% 
% % Create a Bode plot
% bode(tfsys, frequencies);
Y = feedback(tfsys,1);
pole (Y);
bode(tfsys);
margin(tfsys)