lsqnonlin fitting parameters with multiple equations (ODE) and dependency between parameters

For an assignment, I'm trying to fit a set of four ODEs with experimental data. I think I need to use lsqnonlin, but I'm completely lost on how to implement this properly with ODEs and on top of that, the parameters are dependent on each other. I've worked with lsqnonlin before, but this is outside my abilities. The goal is as follows:

I have 4 conditions:

• The concentration of O2- (red line) should be 0 (or close to 0 like 0.001) after 1 sec
• the ratio CO/C2H4 should be 0.62 after 1 sec
• the ratio H2/C2H4 should be 0.64 after 1 sec
• the conversion of CH4 (CH4 (at t=0) - CH4 (at t=1))/CH4 (at t=0)) should be 0.405

as you can see, these are all boundary conditions at 1 sec time. I have 5 parameters that should be fitted, which are constants (k1 to k4) and the initial condition for CH4. Here is the code:

% Parameters to fit
k1=1;
k2=1;
k3=1;
k4=1;
CH4IN=0.076;



% Time interval and initial conditions
tmin = 0;
tmax = 100;
t_interval = [tmin,tmax];
init_cond = [CH4IN,0.076,0,0]'; %initial amount of respectively CH4, O2-, C2H4, H2
% Solution
[t,y] = ode45(@(t,Y) SOEfcn(t,Y,k1,k2,k3,k4) , t_interval , init_cond);
% Calculate and plot fluxes
% Ethane, CO en H2O cant be calculated directly, thus this must be
% Calculated in the following manner:
rEthane=k2.*y(:,3).*y(:,4).^2.*y(:,2);
Ethane = cumtrapz(rEthane);
rCO = 2*k3.*y(:,3).*y(:,2).^4;
CO = cumtrapz(rCO);
rH2O = k1.*y(:,1).^2.*y(:,2)+k2.*y(:,3).*y(:,4).^2.*y(:,2)+2.*k3.*y(:,3).*y(:,2).^4+k4.*y(:,4).*y(:,2);
H2O = cumtrapz(rH2O);

plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'g',t,y(:,4),'y',t,Ethane,'v',t,CO,'p',t,H2O,'x');
legend('CH4','O2-','C2H4','H2','C2H6','CO','H2O')
xlabel('time (s)')
ylabel('concentration')

with SOEfcn (in ode45):

function dYdt = SOEfcn(t,Y,k,k2,k3,k4)
%Y(1) is CH4, Y(2) is oxide (O2-), Y(3) is ethylene (C2H4), Y(4) is H2
dYdt = [-2*k(1)*Y(1)^2*Y(2); %methane rate
         -k(1)*Y(1)^2*Y(2)-k2*Y(3)*Y(4)^2*Y(2)-4*k3*Y(3)*Y(2)^4-k4*Y(4)*Y(2); %oxide rate
         k(1)*Y(1)^2*Y(2)-k2*Y(3)*Y(4)^2*Y(2)-k3*Y(3)*Y(2)^4; %ethylene rate
         %k2*Y(3)*Y(4)^2*Y(2); %ethane rate
         %2*k3*Y(3)*Y(2)^4; %CO rate
         k(1)*Y(1)^2*Y(2)-2*k2*Y(3)*Y(4)^2*Y(2)-k4*Y(4)*Y(2);]; %H2 rate
         %k1*Y(1)^2*Y(2)+k2*Y(3)*Y(4)^2*Y(2)+2*k3*Y(3)*Y(2)^4+k4*Y(4)*Y(2);];%H2O rate
end

Just to be clear, the ethane, CO and H2 rate couldn't be solved within the ode45, so I had to calculate it differently.

Does anyone know how to fit the parameters to the given boundary conditions (I think with lsqnonlin, but if you know another method that's also fine)?

[Torsten]

You forgot to include t in your calls to cumtrapz:

H2O = cumtrapz(t,rH2O)

e.g.

Further, for reaction kinetics equations, you should use a stiff integrator, i.e. ODE15S. I'm confident that with this integrator, you will be able to calculate H2O, CO and ethane directly by the integrator.

And do you really have only those four conditions for optimizing the four parameters ? Or do you have a bunch of measurement data that you did not mention yet ? With only those four conditions, you won't be able do senseful fitting.

[maurits]

Thank you for your remarks, I will alter it accordingly. As for the data, I only have those four conditions which are the final concentrations of the compounds. Basically I need to prove that my reaction rates can achieve the final results. So yes, the fitting will not be accurate which is acceptable for my case.

[Torsten]

Take a look at the examples provided under

de.mathworks.com/help/optim/ug/fit-differential-equation-ode.html?lang=en),

[maurits]

yes i think we solved the problem, thanks for the ode15s tip.

Answer 1

[maurits]

We fixed it, but it turns out our reaction set isn't properly set up so we will have to make a new one but that is outside the scope of this question...

function ke = SOEmodel()
% Parameters to fit
% k1=1;
% k2=1;
% k3=1;
% k4=1;
% CH4IN=0.076;

% Time interval and initial conditions
tmin = 0;
tmax = 1;
t_interval = [tmin,tmax];

%initial values
k=[1,1,1,1,0.076];
fitOpts=optimset('MaxIter',10e7,'MaxFunEvals',10e7,'TolFun',1e-12);
ke = lsqnonlin(@objective, k,[],[],fitOpts);

%----------
    function F = objective(k)
    init_cond = [k(5),0.076,0,0,0,0,0]'; %initial amount of respectively CH4, O2-, C2H4, H2
    % Solution
    [t,y] = ode15s(@(t,Y) SOEfcn(t,Y,k) , t_interval , init_cond);
    F=[(y(end,5)/y(end,3)-0.62),(y(end,6)/y(end,3)-0.64),y(end,2)-0.0001,((k(5)-y(end,1))/k(5))-0.405];
    end
%----------
% Plot fluxes

plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'g',t,y(:,4),'y',t,y(:,5),'v',t,y(:,6),'p',t,y(:,7),'x');
legend('CH4','O2-','C2H4','C2H6','CO','H2','H2O')
xlabel('time (s)')
ylabel('concentration')
end

and

function dYdt = SOEfcn(t,Y,k)
%Y(1) is CH4, Y(2) is oxide (O2-), Y(3) is ethylene, Y(4) is H2
%k1=k(1) en k4=k(4) etc
%Y(2)=1;
dYdt = [-2*k(1)*Y(1)^2*Y(2); %methane rate
         -k(1)*Y(1)^2*Y(2)-k(2)*Y(3)*Y(4)^2*Y(2)-4*k(3)*Y(3)*Y(2)^4-k(4)*Y(4)*Y(2); %oxide rate
         k(1)*Y(1)^2*Y(2)-k(2)*Y(3)*Y(4)^2*Y(2)-k(3)*Y(3)*Y(2)^4; %ethylene rate
         k(2)*Y(3)*Y(4)^2*Y(2); %ethane rate
         2*k(3)*Y(3)*Y(2)^4; %CO rate
         k(1)*Y(1)^2*Y(2)-2*k(2)*Y(3)*Y(4)^2*Y(2)-k(4)*Y(4)*Y(2); %H2 rate
         k(1)*Y(1)^2*Y(2)+k(2)*Y(3)*Y(4)^2*Y(2)+2*k(3)*Y(3)*Y(2)^4+k(4)*Y(4)*Y(2);];%H2O rate
end