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

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.

[Torsten]

Since you want to estimate activation energies, I hope you have measurement data for different temperatures, but under adiabatic conditions.

https://www.tutorchase.com/answers/ib/chemistry/how-can-the-activation-energy-be-determined-experimentally

[Continuum]

In this study, we are taking Q to be a constant at a specific temperature of interest. Same applies to QG. Thanks.

[Torsten]

Your data seem to indicate that you are performing non-adiabatic experiments where temperature increases over time because of the heat of reaction. That's inappropriate to estimate activation energies and preexponential factors.

[Continuum]

Yes, you are right that there are different activation energies at the different temperatures, but for this study, the relevant process is considered only at a single temperature. So, what I am looking for are estimates of these parameters to get a 'feel' of the model vs my current data . I will then compare the obtained model with data from another experiment of the same material (think of it like a different specimen). Thanks.

[Torsten]

but for this study, the relevant process is considered only at a single temperature.

Then you must refer to data different from those that you included as modified_data.xlsx here:

https://uk.mathworks.com/matlabcentral/answers/314965-parameter-estimation-for-a-system-of-differential-equations?s_tid=srchtitle

[Continuum]

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

[Torsten]

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]

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.

.

[Continuum]

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]

As always, my pleasure!

Answer 2

[William Rose]

@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.

[William Rose]

@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]

Thanks for the pointers, I reattached the code, sorry about that.