How to find fitting optimized parameters to fit a system of non-linear ODEs to experiment.

Question

Continuum 2024-1-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2074946-how-to-find-fitting-optimized-parameters-to-fit-a-system-of-non-linear-odes-to-experiment

编辑： Torsten 2024-1-27

Hi

I have a set of ODEs (attached), I have been able to solve them using ode45, however, my issue now is my experimental results don't match the integrated values of the equations. So, I am looking to fit only the solution for epsilon with it's experimental results to find the best parameters A, B, (A0/alpha), k0, Q, and QG. Attached is my code based on an answer from another thread, but it just runs continuously but I couln't figure out what the problem is. Could it be that the there are too many parameters to fit? Any help is greatly appreciated. Thank you.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Continuum 2024-1-27

If you observe the data, there is a period where the temp stayed almost constant for some time.

Torsten 2024-1-27

编辑：Torsten 2024-1-27

But it's one long experiment with a dynamic development in time. So your measurements have a history. It's usually necessary that you start at t = 0 with a fixed temperature which is kept constant over time until the experiment has finished.

But we are in a MATLAB forum here ...

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2024-1-27

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2074946-how-to-find-fitting-optimized-parameters-to-fit-a-system-of-non-linear-odes-to-experiment#answer_1398266

在 MATLAB Online 中打开

You will need to post the ‘modified_data.xlsx’ file as well to run this here.

My edit of your code —

Model_Exp_Fit
function Model_Exp_Fit
% 2016 12 03
% NOTES: 
% 
%   1. The theta (parameter) argument has to be first in your 
%       kinetics funciton,
%   2. You need to return ALL the values from DifEq since you are fitting
%       all the values
        dataTable1=readtable('modified_data.xlsx','Sheet','Sheet1', 'VariableNamigRule','preserve');
        dataTable1 = dataTable1(:,4:end);
        tempFunc = griddedInterpolant(dataTable1.time1_min_,dataTable1.T1_K_);
        % Temp = tempFunc(t);
        % tT = [t Temp];
        tT = dataTable1{:,[1 2]};                       % Time And Temperature
        t = tT(:,1);
        x = dataTable1{:,3};                            % Epsilon
        % return
    function X = evolution(theta,tT)
    T = tT(:,2);
    tv = tT(:,1);
    x0=[20;0.4363;0.0000001];
    [T,Xv] = ode15s(@DifEq,t,x0);
    k0 = theta(1);
    QG = theta(2);
    A0_alpha = theta(3);
    Q = theta(4);
    A = theta(5);
    B = theta(6);
        function dX = DifEq(t,x)
        %thet_c = 0.508; 
        rho_c = 0.948; 
        R = 8.314;
        Temp = interp1(tv, T, t);                   % Interpolates 'T' For Each Value Of 't'
        % dataTable1=readtable('modified_data.xlsx','Sheet','Sheet1')
        % tempFunc = griddedInterpolant(dataTable1.time1_min_,dataTable1.T1_K_);
        % Temp = tempFunc(t);
        % tT = [t Temp];
        dxdt = zeros(3,1);
        phi1 = (theta(1)/3)*(1-rho_c)^1.5;
        phi2 = 9/(2*theta(3));
        phi3 = 3/(2*theta(3));
        dxdt(1) = (exp(-theta(2)./(R*Temp))).*(phi1./x(1).^2).*1./((1-rho_c+x(2)).^1.5);
        dxdt(2) = (-phi2./(x(1).*Temp)).*(1./exp(theta(4)./(R*Temp))).*x(2).^theta(6).*(((1-x(2)).^3)./(1-x(2)).^theta(5));
        dxdt(3) = (-phi3./(x(1).*Temp)).*(1./exp(theta(4)./(R*Temp))).*x(2).^theta(6).*(((1-x(2)).^2)./(1-x(2)).^theta(5));
        dX = dxdt;
        end
    X=Xv(:,3);
    end
% return
% eps_exp = load("exp_epsilon.txt");
% t = eps_exp(:,1);
% t = [0;t];
% % x = eps_exp(:,2)/100;
% x = [0; x];
theta0 = [(29.65e-5)/(1/60);164.8e3;2.03*((1/60)/(10^6));217.2e3;11.35;0.49];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@evolution,theta0,tT,x);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
% tv = linspace(min(t), max(t));
Xfit = evolution(theta, tT);
figure(1)
plot(t, x, 'p')
hold on
hlp = plot(t, Xfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
% legend(hlp, 'G(t)', 'theta(t)', 'epsilon(t)', 'Location','N')
end

I left in the sections that I commented-out so that you can understand my edits. I will give this a shot with some standardized code that I use to run genetic algorithm parameter estimations to see if I can improve on these results. If I get good results with it, I will post the results and the code here. For the time being, this code runs, and you can experiment with it.

.

2 个评论
显示无隐藏无

Continuum 2024-1-27

Thanks very much, it works now, I will play around with it a bit and see if I can get some relatively good results. Thanks very much for your input.

Star Strider 2024-1-27

As always, my pleasure!

请先登录，再进行评论。

Answer 2

William Rose 2024-1-27

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2074946-how-to-find-fitting-optimized-parameters-to-fit-a-system-of-non-linear-odes-to-experiment#answer_1398261

@Buhari Ibrahim

Matlab has excellent tools for this.

https://www.mathworks.com/help/optim/ug/optimizing-a-simulation-or-ordinary-differential-equation.html

and

https://www.mathworks.com/help/optim/ug/fit-differential-equation-ode.html

explain ways to do it.

2 个评论
显示无隐藏无

William Rose 2024-1-27

@Buhari Ibrahim, You are fitting six parameters. That is not too many. I have fitted ODE models with more parameters. You say you have attached code, but none is attached.

Continuum 2024-1-27