Differential equation fit to my data: incorrect minimization

1 view (last 30 days)
Marina Batlló RIus on 29 Mar 2023
Answered: Alan Weiss on 29 Mar 2023
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?

Alan Weiss on 29 Mar 2023
I think that you probably did not set up your problem correctly. There are relevant documentation examples:
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