Hi everyone,
I am currently trying out MATLAB fitting to model a series chemical kinetics, which is based on an existing question in this forum seen here.
In my case, I add parameter of temperature to determine the reaction constant ki. The problem is intended to be solved using lsqnonlin.
Thus, the objective is to fit the data and parameters to find the values of Arrhenius coefficients A1-4 and Ea1-4.
Currently, I am having a code (seen at the end of this post) that generates some errors and can't be executed, so there I must've done some mistakes in either coding or my way of thinking. I looked up and refer to similar problems in the forum as well, but it seems it hasn't been working for my case so far.
I'm still pretty new in Matlab, especially in fitting nonlinear models. I hope someone can help me find the mistakes and give a solution to the problem.
Thanks in advance.
Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf. > In checkbounds (line 33)In lsqnsetup (line 77)In lsqnonlin (line 207)In Untitled2_trial3 (line 53) Warning: Length of upper bounds is > length(x); ignoring extra bounds. > In checkbounds (line 41)In lsqnsetup (line 77)In lsqnonlin (line 207)In Untitled2_trial3 (line 53) Exiting due to infeasibility: at least one lower bound exceeds the corresponding upper bound.>>
%Reaction system
%c1 --> c2
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The data for p1 and p2 is unknown.
%Data
time = [10 20 40 60 100 140 180];c1 = [0.508 0.442 0.385 0.354 0.299 0.267 0.258];c2 = [0.246 0.301 0.359 0.398 0.422 0.467 0.489];T_o = [400 423 433 456 467 468 463];c1_0 = 0.711;c2_0 = 0;T_o_0 = 400;R = 8.314;%Plot
time_plot = [0,time];c1_plot = [c1_0,c1];c2_plot = [c2_0,c2];T_o_plot = [0,T_o];%figure 1, C1 and C2 vs time
figure(1)plot(time_plot,c1_plot,'bo')grid onhold onplot(time_plot,c2_plot,'ro')hold offxlabel('Time [min]')ylabel('Concentration [mol/L]')legend('c1','c2')%figure 2, C1 and C2 vs T
figure(2)plot(T_o_plot,c1_plot,'bo')grid onhold onplot(T_o_plot,c2_plot,'ro')hold offxlabel('Temperature [K]')ylabel('Concentration [mol/L]')legend('c1','c2')%Data processing
x_in = [time',T_o'];y_in = [c1',c2'];%Lsqnonlin
%[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqnonlin(@diff1,p0,x_in,y_in,lb,ub);
%FitData = diff1(p,time,1);
%c1_out = FitData(:,1);
%c2_out = FitData(:,2);
% Options for fitting
options = optimoptions('lsqnonlin', ... 'Display', 'iter', ... 'Algorithm', 'levenberg-marquardt', ... 'UseParallel', false, ... 'StepTolerance', 1e-6);%initial value for iteration
B0 = [0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,c1(1),c2(1),T_o(1)];% lsqnonlin
optimizedParams = lsqnonlin(@diff1, B0, x_in, y_in, [], [], options);function C = diff1(p,time,~)%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%where
%k1 = A1*exp(-E1/RT_o)
%k2 = A2*exp(-E2/RT_o)
%k3 = A3*exp(-E3/RT_o)
%k4 = A4*exp(-E4/RT_o)
%
%variables: c1 = x(1), c2 = x(2), T_o = x(3)
%parameters:
%A1 = B(1), E1 = B(2), n1 = B(3)
%A2 = B(4), E2 = B(5), n2 = B(6)
%A3 = B(7), E3 = B(8), n3 = B(9)
%A4 = B(10), E4 = B(11), n4 = B(12)
x0 = B(11:12);[t,Cout] = ode45(@DifEq, time, x0); function dC = DifEq(t, x) xdot = zeros(3,1); xdot(1) = -(B(1)*exp(-B(2)/R/x(3)).*x(1).^B(3)-B(7)*exp(-B(8)/R/x(3)).*x(1)^B(9)); xdot(2) = B(4)*exp(-B(5)/R/x(3)).*x(1).^B(6)-B(10)*exp(-B(11)/R/x(3)).*x(2).^B(12); dC = xdot; endC = Cout;end
Best Answer