Using ode45 to create simulated data and globally optimize the simulated data to two experimental data sets?

Question

jbs427 2016-6-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/292778-using-ode45-to-create-simulated-data-and-globally-optimize-the-simulated-data-to-two-experimental-da

评论： Star Strider 2016-7-19

Research area: Chemical kinetics

Problem description: I have acquired reaction rate versus time (v-t) and photocurrent versus time (j-t) data. The ultimate goal is to simulate data and fit the experimental v-t and j-t data with the following expressions:

v(t)=k_AR*X(t), where X(t) will be determined by solving the differential equations below from the kinetic model.

The photocurrent versus time data will be fit using the following expression:

i(t)=q*A*C1[k1*z(t)*(X(t)+y(t))-k2*n(t)*X(t)], where q,A, C1 are constants, z(t), y(t), and n(t) will be determined by solving the differential equations below from the kinetic model.

I would like to use ode45 to solve the following set of ODEs and then minimize the error between the simulated results and the experimental data (using lsqcurvefit?) by optimizing k1, k2, k3, k_AR in the system of 5 differential equations below:

Eq. 1. dw/dt=k2*n*C1*X-k1*C1*w*z Eq. 2. dz/dt=G-k1*C1*w*z Eq. 3. dX/dt=k1*C1*w*z-2k3*C1*X^2-k2*n*C1*X Eq. 4. dy/dt=k3*C1*X^2-k1*z*C1*y Eq. 5. dn/dt=R-k2*n*C1*X

C1,G and R are known constants. X+y+z=1 ... (I am having trouble applying this condition in ODE solver) Initial Conditions are: X(0)=0, y(0)=0, w(0)=1, n(0)=n1, z(0)=z1, where n1 and z1 are initial numbers. I have tried the following routine to set up the ode solver:

syms w(t) z(t) X(t) y(t) n(t) 
guess=[Gen0 1 0 0 Elec0];  %%z w X y n
eqn1=diff(w,t) == (k2*n*C1*X)-(k1*C1*w*z);
eqn2=diff(z,t) == Gen-(k1*C1*w*z);
eqn3=diff(X,t) == (k1*C1*w*z)-(2*k3*C1*X^2)-(k2*n*C1*X); 
eqn4=diff(y,t) ==(k3*C1*X^2)-(k1*z*C1*y); %%dy/dt
eqn5=diff(n,t) == Elec-(k2*n*C1*X); %%dn/dt
[V,Y]=odeToVectorField([eqn1], [eqn2], [eqn3], [eqn4], [eqn5]);
M = matlabFunction(V,'vars',{'t','Y'});
sol = ode45(M,[.016 50],guess);
xaxis=[0.016:0.016:3125*0.016];
OHradconc=deval(sol,xaxis,[3]);
simulatedv=OHradconc.*k_AR;

I would then determine the "simulatedj" and I would like to minimize the error between simulatedv and experimental v-t as well as minimize the error between simulatedj and j-t. In this way, I would use two experimental data sets to determine k1, k2, k3, and k_AR.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2016-6-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/292778-using-ode45-to-create-simulated-data-and-globally-optimize-the-simulated-data-to-two-experimental-da#answer_227042

Using the Symbolic Math Toolbox is probably not the best approach to this problem. We did something similar in: Monod kinetics and curve fitting. You will likely find it helpful.

2 个评论
显示无隐藏无

jbs427 2016-7-18

在 MATLAB Online 中打开

I followed the Monod kinetics and curve fitting example above for my case (to the best of my ability!), I wrote the following "MyMonodKinetics" function below. I have the following questions:

1. The output of the solver that I am particularly interested in is X(t), eq 3. How do I "optimize" the kAR parameter? Should I put it outside or inside the differential equations?

2. The output of MyMonodKinetics is the vector rate (sim_rate), and it should have the same number of points as my experimental rate data (exp_rate). I want to write a function that will find the optimum k1,k2,k3, and kAR so that the error (exp_rate-sim_rate)^2 is minimized. I am so lost here and would really appreciate your help! I've been working on this as best as I can, without success.

     function rate = MyMonodKinetics(B, t)
    %%The goal is to fit simulated rate(t) to experimental rate(t) data, 
    %%where simulated rate(t)=kAR*X(t) and X(t) will be the third column of the ode45 solver solution. eq. 3 below.
    % Following MONODKINETICS1: codes the system of differential equations:
    %         eq 1.  dw/dt = (k2*n*C1*X)-(k1*C1*w*z);
    %         eq 2.  dz/dt = Gen-(k1*C1*w*z);
    %         eq 3.  dX/dt = (k1*C1*w*z)-(2*k3*C1*X^2)-(k2*n*C1*X); %%looking to match this to expt data.
    %         eq 4.  dy/dt = (k3*C1*X^2)-(k1*z*C1*y); %%dy/dt
    %         eq 5.  dn/dt = Elec-(k2*n*C1*X); %%dn/dt
    %         
    %   with:
    %       Constants: Gen = a number, C1 = a number, Elec = a number.
    %       Variables:  x(1) = n,  x(2) = X, x(3) = w, x(4) = z, x(5) = y 
    %       Parameters: k2 = B(1),  k1 = B(2),  k3 = B(3)
    %%the constants
    Gen0=43.3841;
    Elec0=2.8808e-08;
    x0 = [Gen0 1 0 0 Elec0]'; %%initial conditions
    [T,Sv] = ode45(@DifEq3, t, x0);
          function dS = DifEq3(t, x)
          %%some constants      
              Elec=3.4640e-06;
              C1=1e15/1e14;
              Gen=5.3190e+03;
              xdot = zeros(5,1);
              xdot(1) = (B(1) .* x(1) .* C1.*x(2)) - (B(2) .* C1.*x(3).*x(4));
              xdot(2) = Gen-(B(2) .* C1.*x(3) .* x(4));
              xdot(3) = (B(2).*C1.*x(3).*x(4))-(2.*B(3).*C1.*x(2).^2)-(B(1).*x(1).*C1.*x(2)); %%looking to match this to expt data.
              xdot(4) = (B(3).*C1.*x(2).^2)-(B(2).*x(4).*C1.*x(5));
              xdot(5) = Elec-(B(1).*x(1).*C1.*x(2)); 
              dS = xdot;
          end
    %%generate the rate data with same num of points as the experimental data  
    xaxis=[0.016:0.016:3125*0.016]; 
    OHradconc=deval(sol,xaxis,[3]);
    rate=OHradconc.*kAR;  %%how can I optimize kAR here???
    end

Star Strider 2016-7-19

在 MATLAB Online 中打开

In the Monod Kinetics example, the person was fitting the first column, corresponding to ‘[S]’, so the output returned to the objective function ‘Sv(:,1)’. If you want to fit ‘[X]’, the ‘S’ assignment in the Monod Kinetics code becomes:

S = Sv(:,3);

I admit I don’t understand this part of the code:

    %%generate the rate data with same num of points as the experimental data  
    xaxis=[0.016:0.016:3125*0.016]; 
    OHradconc=deval(sol,xaxis,[3]);
    rate=OHradconc.*kAR;  %%how can I optimize kAR here???

It’s not necessary, and will likely interfere with the correct function of the code.

All the parameters will be optimised, so you don’t have to select particular parameters for optimisation.

Also, lsqcurvefit can fit matrices (dependent data with multiple columns) if you want to fit more than one variable. You have to make appropriate changes in your code to return them as output from the objective function that lsqcurvefit uses. (In the Monod Kinetics example, ‘MonodKinetics1’ is the objective function.)

请先登录，再进行评论。