Differential equation fit to my data: incorrect minimization

So I have been trying to fit some experimental data that I have in my differential equations. The code I am using is as follows:

Time = table2array(ValuesMaxPeaks(:,1));
x0 = table2array([ValuesMaxPeaks(1,2),ValuesMaxPeaks(1,3),ValuesMaxPeaks(1,4),ValuesMaxPeaks(1,5)]);
Sdata = table2array([ValuesMaxPeaks(:,2),ValuesMaxPeaks(:,3),ValuesMaxPeaks(:,4),ValuesMaxPeaks(:,5)]);

thau = 1/26.61;

model_funct = @(k) ode45(@(t,y) [k(3)*y(3)+k(1)*y(2)+k(5)*y(4)-(k(2)+k(4)+k(6)+thau)*y(1);
             k(2)*y(1)-(k(1)+thau)*y(2);
             k(4)*y(1)-(k(3)+thau)*y(3);
             k(6)*y(1)-(k(5)+thau)*y(4)], Time, x0);

error = @(k) sum(sum((model_funct(k).y - Sdata').^2));


options.StepTolerance = 0.0000000000000000000000000001;
options.MaxFunctionEvaluations = 12000000000000000000000;
options.FunctionTolerance = 0.000000000000000000000000000000000000000001;
k0 = [0.1, 0.01, 0.1, 0.1, 0.1, 0.1]; % Initial guess for k
k_opt = lsqnonlin(@(k) errorfunc(k,Sdata, model_funct),k0,zeros(1,6),[0.2,0.2,0.2,0.2,0.2,0.2], options);

As you can see, I use ode45 to solve the differential equations and lsqnolin to minimize the error of the solutions compared to my experimental data. The code does not return any error, but when I run it, MATLAB returns me the following:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point
is less than 1e-4 times the value of the function tolerance.

<stopping criteria details>

and the result of my minimization stays as my initial values k0. I have already tried to lower the function tolerance and step tolerance, but do not manage to make it solve my constants properly. I also attach some example experimental data in case anyone needs to see it. Do you have any suggestions?

Answer 1

[Alan Weiss]

I think that you probably did not set up your problem correctly. There are relevant documentation examples:

https://www.mathworks.com/help/optim/ug/fit-differential-equation-ode.html

https://www.mathworks.com/help/optim/ug/fit-ode-problem-based-least-squares.html

Here is a little script I wrote based on yours. It doesn't work very well, in that it produces a poor fit.

Alan Weiss

MATLAB mathematical toolbox documentation

load ValuesMaxPeaks.csv % I created a comma-separated array, attached
Time = ValuesMaxPeaks(:,1);
x0 = [ValuesMaxPeaks(1,2),ValuesMaxPeaks(1,3),ValuesMaxPeaks(1,4),ValuesMaxPeaks(1,5)];
Sdata = [ValuesMaxPeaks(:,2),ValuesMaxPeaks(:,3),ValuesMaxPeaks(:,4),ValuesMaxPeaks(:,5)];

thau = 1/26.61;

%options.StepTolerance = 0.0000000000000000000000000001;
options = optimoptions("lsqnonlin",'FiniteDifferenceStepSize',1e-4,...
    'FiniteDifferenceType','central');
%options.MaxFunctionEvaluations = 12000000000000000000000;
%options.FunctionTolerance = 0.000000000000000000000000000000000000000001;
k0 = [0.1, 0.01, 0.1, 0.1, 0.1, 0.1]; % Initial guess for k
% pos = paramfun(k0,Time,x0,thau)
fun = @(k)paramfun(k,Time,x0,thau) - Sdata;
k_opt = lsqnonlin(fun,100*rand(1,6),zeros(1,6),100*ones(1,6), options);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

plot(paramfun(k_opt,Time,x0,thau))
hold on
plot(Sdata,'.')
hold off

function pos = paramfun(k,Time,x0,thau)
f = @(t,y) [k(3)*y(3)+k(1)*y(2)+k(5)*y(4)-(k(2)+k(4)+k(6)+thau)*y(1);
             k(2)*y(1)-(k(1)+thau)*y(2);
             k(4)*y(1)-(k(3)+thau)*y(3);
             k(6)*y(1)-(k(5)+thau)*y(4)];
[~,pos] = ode45(f,Time,x0);
end