MATLAB: Help needed on ODE solving coupled with fmincon

fminconode

Hi, I need help to solve the problem with the code for my ODE solving. As I have many parameters and I have used lsqcurvefit to solve my problem, but there was an error message of:
"Solver stopped prematurely.
lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 800 (the default value)."
function ODE
% 2019 12 4
    function X=kinetics(k,t)
x0=[7.525;13.9;3;0;0;0];
[T,Xv]=ode45(@DifEq,t,x0);
%
    function dX=DifEq(t,x)
    dxdt=zeros(6,1);
    dxdt(1)= -(k(1)+k(6)).*x(1);
    dxdt(2)= -(k(2)+k(7)).*x(2);
    dxdt(3)= -k(8).*x(3);
    dxdt(4)= k(1).*x(1)-k(3).*x(4)-k(4).*x(4);
    dxdt(5)= k(2).*x(2)+k(3).*x(4)-k(5).*x(5);
    dxdt(6)= k(6).*x(1)+k(7).*x(2)+k(8).*x(3)+k(4).*x(4)+k(5).*x(5);
    dX=dxdt;
    end
X=Xv;
end
t=[0
5
10
15
20];
x=[7.525 13.9 3 0 0 0
5.585 9.51 8 1.44 2.89 2.0
5.275 8.27 7 1.65 3.88 2.35
4.925 7.93 6 1.90 3.97 2.7
3.885 6.00 5.5 2.89 5.70 2.95]; 
k0=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,x);
fprintf(1,'\tRate Constants:\n')
for n1 = 1:length(k)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', n1, k(n1))
end
tv = linspace(min(t), max(t));
Xfit = kinetics(k, tv);
figure(1)
plot(t, x, 'p')
hold on
hlp = plot(tv, Xfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'X_1(t)', 'X_2(t)', 'X_3(t)', 'X_4(t)', 'X_5(t)', 'X_6(t)', 'Location','N')
end
Instead, I followed advice of another expert to use fmincon to substitute lsqcurvefit to evaluate the function to find k 
k0=[0.25;0.25;0.25;0.25;0.25;0.25;0.25;0.25];
lb=[0;0;0;0;0;0;0;0];
ub=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5];
A=[];
b=[];
Aeq=[];
beq=[];
k = fmincon(@kinetics,k0,A,b,Aeq,beq,lb,ub);
But there was also an error message of:
"Not enough input arguments.
Error in kinetics
[T,Xv]=ode45(@DifEq,t,x0);
Error in fmincon
initVals.f = feval(funfcn{3},X,varargin{:});
Error in call_fmincon
k = fmincon(@kinetics,k0,A,b,Aeq,beq,lb,ub);
Caused by:
Failure in initial objective function evaluation. FMINCON cannot continue."

Best Answer

It is best to fit the initial conditions as well as the kinetics parameters. I would not put an upper bound on the kinetics parameters. A lower bound to be certain they are all is all that is necessary.
This is something of an improvement in my original code (I upgraded it later, including a change so the fitted lines and markers are the same colours): 
    function X=kinetics(k,t)
% x0=[7.525;13.9;3;0;0;0];
x0 = k(9:14);                                                           % Fit The Initial Conditions As Well
[T,Xv]=ode45(@DifEq,t,x0);
%
    function dX=DifEq(t,x)
    dxdt=zeros(6,1);
    dxdt(1)= -(k(1)+k(6)).*x(1);
    dxdt(2)= -(k(2)+k(7)).*x(2);
    dxdt(3)= -k(8).*x(3);
    dxdt(4)= k(1).*x(1)-k(3).*x(4)-k(4).*x(4);
    dxdt(5)= k(2).*x(2)+k(3).*x(4)-k(5).*x(5);
    dxdt(6)= k(6).*x(1)+k(7).*x(2)+k(8).*x(3)+k(4).*x(4)+k(5).*x(5);
    dX=dxdt;
    end
X=Xv;
end
t=[0
5
10
15
20];
x=[7.525 13.9 3 0 0 0
5.585 9.51 8 1.44 2.89 2.0
5.275 8.27 7 1.65 3.88 2.35
4.925 7.93 6 1.90 3.97 2.7
3.885 6.00 5.5 2.89 5.70 2.95]; 
k0=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;rand(6,1)];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,x,zeros(size(k0)));
fprintf(1,'\tRate Constants:\n')
for n1 = 1:length(k)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', n1, k(n1))
end
tv = linspace(min(t), max(t));
Xfit = kinetics(k, tv);
figure(1)
hd = plot(t, x, 'p');
for k1 = 1:size(x,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Xfit);
for k1 = 1:size(x,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'X_1(t)', 'X_2(t)', 'X_3(t)', 'X_4(t)', 'X_5(t)', 'X_6(t)', 'Location','N')
This takes a while to run, however it produces a better result.
If you have the Global Optimization Toolbox, I also wrote a version of this code using the ga (genetic algorithm) function that is generally much more robust at finding the optimal parameter set, although it can take a long time to run (as ga optiomisations usually do). I will post it as an edit to my Answer here if you request it.
Related Question