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

评论： Star Strider 2024-6-26

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.

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Yeuran 2024-6-26

@Star Strider Hello, I have another question about your code. What should we do if the experimental data(epsilon) is shorter than time and Temp? Do you have any suggestions? For instance, if 'x' contains only 3 data points (like x = dataTable1{[1,5,7],3}), and 'tT' and 't' are of equal legth(where tT = dataTable1{:,[1 2]} and t = tT(:,1)), how should we proceed to use 'lsqcurvefit' in such a scenario?

Star Strider 2024-6-26

@Yeuran — All data have to have the same number of elements (in this instance, rows). I would only use the data with corresponding with the rows in ‘x’, deleting everything else in the longer data.

I would not interpolate or extrapolate to equalise them, since in both instances that creates data where there were no data previously, and the parameter estimates are then useless.

So, if you only have three data points, you can only estimate at most two parameters. If you are estimating more parameters than that, you will need to repeat the experiment and get more data. There are simply no other optioons.

请先登录，再进行评论。

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