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

31 views (last 30 days)

Show older comments

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 Comments

Torsten
on 27 Jan 2024

Edited: Torsten
on 27 Jan 2024

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

### Accepted Answer

Star Strider
on 27 Jan 2024

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 Comments

### More Answers (1)

William Rose
on 27 Jan 2024

Matlab has excellent tools for this.

and

explain ways to do it.

##### 2 Comments

William Rose
on 27 Jan 2024

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!