How can i solve this problem' when i am running, this message appears: Failure in initial objective function evaluation. FMINUNC cannot continue.

I'm writing a matlab program to try to fit a mathematical model to real data.

Answer 1

[Torsten]

Call your function in which you define the sum of differences squared between real and modelled data for the initial values of your parameters and see what the problem is. Maybe the function returns something undefined or a vector of values instead of a single value.

[Khadija]

It returns:

>> moderCVD

Unrecognized function or variable 'k'.

Error in moderCVD>@(t,y)(model_CVD(t,y,k)) (line 25)

[T Y] = ode23s(@(t,y)(model_CVD(t,y,k)),tforward,[ 52204.0 2000.0 36728.0 68.0 1000.0]);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode23s (line 121)

= odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in moderCVD (line 25)

[T Y] = ode23s(@(t,y)(model_CVD(t,y,k)),tforward,[ 52204.0 2000.0 36728.0 68.0 1000.0]);

[Torsten]

It seems that "k" is undefined when you call ode23s in function "modercvd".

And usually "modercvd" should return a value, namely the sum of differences squared between your real data and your mode data, shouldn't it ?

[Khadija]

Your'you're right; how can i solve this problem

[Torsten]

By defining k and returning the sum of differences squared between your real data and your model data.

That's all that can be said without seeing your code and your data.

[Khadija]

can i send you my code?

[Khadija]

function dy=model_CVD(~,y,k)

delta=0.02 ;

deltaC=0.03 ;

h=1244;

xi = k(1);

sigma = k(2);

alpha= k(3);

gammaH= k(4);

gammaA= k(5);

deltaOC= k(6);

p1= k(7);

p2= k(8);

p3= k(9);

psi=k(10);

dy = zeros(5,1);

dy(1)= h-(delta+sigma+xi+gammaH)*y(1);

dy(2)=sigma*y(1) -(delta+alpha+gammaA)*y(2);

dy(3)=gammaH*y(1)+ psi*y(4) -(delta+deltaOC+p1*sigma+p3*xi)*y(3);

dy(4)=gammaA*y(2)+ p1*sigma*y(3) -(psi+delta+p2*alpha)*y(4);

dy(5)=xi*y(1)+alpha*y(2)+p2*alpha*y(4)+p3*xi*y(3) -(delta+deltaC)*y(5);

end

function error_in_data = moderCVD(k)

% Données spécifiques

specific_data = [

2014 52204 2000 36728 68;

2015 43000 2007 41899 94;

2016 30400 3124 51385 91;

2017 26300 1320 55583 97;

2018 20800 2153 59948 99;

2019 16613 2000 64288 99;

2020 13500 2523 65877 100;

2021 12700 2154 67446 10;

2022 10150 1473 68279 100

];

tforward = 2014:0.01:2022;

tmeasure = [1:100:901]';

% Utilisez les données spécifiques

tdata = specific_data(:, 1);

Hdata = specific_data(:, 2);

Adata=specific_data(:, 3);

THdata=specific_data(:, 4);

TAdata=specific_data(:, 5);

[T, Y] = ode23s(@(t,y)(model_CVD(t,y,k)),tforward,[ 52204.0 2000.0 36728.0 68.0 1000.0]);

Hq = Y(tmeasure(:),1);

Aq = Y(tmeasure(:),2);

THq = Y(tmeasure(:),3);

TAq = Y(tmeasure(:),4);% assignts the y-coordinates of ...

%the solution at

A=(Hq - Hdata).^2;

error_in_data =sum(A);

%%

end

clear all ; clf ; clc ;

% Données spécifiques

specific_data = [

2014 52204 2000 36728 68;

2015 43000 2007 41899 94;

2016 30400 3124 51385 91;

2017 26300 1320 55583 97;

2018 20800 2153 59948 99;

2019 16613 2000 64288 99;

2020 13500 2523 65877 100;

2021 12700 2154 67446 10;

2022 10150 1473 68279 100

];

tforward = 2014:0.01:2022;

tmeasure = [1:100:801]';

% Utilisez les données spécifiques

tdata = specific_data(:, 1);

Hdata = specific_data(:, 2);

Adata = specific_data(:, 3);

THdata = specific_data(:, 4);

TAdata = specific_data(:, 5);

% Paramètres initiaux

%deltaOC= 0.01;

p1= 0.5;

delta=0.016+0.0004 ;

deltaC=0.3 ;

h=1244;

sigma=0.06 ;

xi= 0.02;

alpha= 0.04;

epsilon= 0.03;%epsilon=xi*1.5

beta= 0.06;% beta=alpha*1.5

phi= 0.03;%phi=sigma*0.05

psi= 0.3;

gammaH=0.2 ;

gammaA=0.3 ;

delta=0.0016 ;

deltaA=0.0004 ;

deltaC=0.01 ;

deltaOC= 0.01;

p2=1.5;

p3=2;

k0 = [xi sigma alpha gammaH gammaA deltaOC p1 p2 p3 psi];

% Résolution de l'équation différentielle avec les paramètres initiaux

[T, Y] = ode23s(@(t, y) model_CVD(t, y, k0), tforward, [52204.0 2000.0 36728.0 68.0 1000.0]);

yintH = Y(tmeasure(:),1);

yintA = Y(tmeasure(:),2);

yintTH = Y(tmeasure(:),3);

yintTA = Y(tmeasure(:),4);

% Affichage des données spécifiques et de la solution de l'équation différentielle

% Display the results

figure(1)

%subplot(1,2,1);

plot(tdata, Hdata, 'r.');

hold on

plot(tdata, yintH, 'b-');

xlabel('time in days');

ylabel('Number of cases');

legend('Data', 'Model estimation');

axis([2014 2022 0 100000]);

% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)

% Assigns the new values of parameters to k and the error to fval

% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)

% Assigns the new values of parameters to k and the error to fval

[k,fval] = fminunc(@moderCVD,k0);

%print final values of fitted parametres

disp(k);

xi=k(1);

sigma=k(2);

alpha= k(3);

gammaH=k(4) ;

gammaA=k(5) ;

deltaOC= k(6);

p1= k(7);

p2= k(8);

p3= k(9);

psi= k(10);

%Draw the data with the final ODE

[T Y] = ode45(@(t,y)(model_CVD(t,y,k)),tforward,[52204.0 2000.0 36728.0 68.0 1000.0]);

yintH = Y(tmeasure(:),1);

yintA = Y(tmeasure(:),2);

yintTH = Y(tmeasure(:),3);

yintTA = Y(tmeasure(:),4);

residuals = Hdata - yintH;% a chercher comment calcumer les residus pour plusieurs variables

%subplot(1,2,2);

figure(2)

plot(tdata,Hdata,'r.');

hold on

plot(tdata,yintH,'g-');

xlabel('Time in days');

ylabel('Number of newly infected cases with fitted parameters');

axis([2014 2022 0 100000]);

%Graphe des résidus

%subplot(1,2,2);

figure(3)

plot(tdata, residuals, 'r+');

xlabel('Time in days');

ylabel('Residuals');

title('Residuals Plot');

axis([2014 2022 min(residuals) max(residuals)]);

% % Affichage de la légende

legend('Residuals');

%

% % Vous pouvez également ajouter une ligne horizontale à zéro pour référence

% hold on

% plot([3, 14], [0, 0], 'k--');

%

%

%Graphe des résidus

%subplot(1,2,2);

%figure(3)

%plot(tdata, residuals, 'g.');

%xlabel('Time in days');

%ylabel('Residuals');

%title('Residuals Plot');

%axis([2014 2022 min(residuals) max(residuals)]);

% Affichage de la légende

%legend('Residuals');

% Ajout de lignes verticales reliant chaque point de résidu à la ligne horizontale à zéro

%hold on

for i = 1:length(tdata)

%plot([tdata(i), tdata(i)], [0, residuals(i)], 'k--');

end

% Ajout de la ligne horizontale à zéro

%plot([2014, 2022], [0, 0], 'k--');

[Torsten]

Sorry, no, I only answer questions here in the forum.

But if you have a valid licence, you could contact MATLAB support.

[Torsten]

Most probably you mean something like this:

clear all ; clf ; clc ;


% Données spécifiques
specific_data = [
    2014   52204 2000 36728 68;
    2015   43000  2007 41899 94;
    2016    30400 3124 51385 91;
    2017  26300 1320 55583 97;
    2018   20800 2153  59948 99;
    2019  16613 2000 64288 99;
    2020   13500 2523 65877 100;
    2021   12700 2154 67446 10;
    2022   10150 1473 68279 100
    
    ];

tforward = 2014:1:2022; 
%tmeasure = [1:100:801]'; 
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
Adata = specific_data(:, 3);
THdata = specific_data(:, 4);
TAdata = specific_data(:, 5);
% Paramètres initiaux


%deltaOC= 0.01;

p1= 0.5;


delta=0.016+0.0004 ;
deltaC=0.3 ;
h=1244;

sigma=0.06 ; 
xi= 0.02;  
alpha= 0.04;
epsilon= 0.03;%epsilon=xi*1.5
beta= 0.06;% beta=alpha*1.5
phi= 0.03;%phi=sigma*0.05
psi= 0.3;
gammaH=0.2 ;
gammaA=0.3 ;
delta=0.0016 ;
deltaA=0.0004 ;
deltaC=0.01 ;
deltaOC= 0.01;
p2=1.5;
p3=2;
k0 = [xi sigma alpha gammaH gammaA deltaOC p1 p2 p3 psi];

% Résolution de l'équation différentielle avec les paramètres initiaux
[T, Y] = ode23s(@(t, y) model_CVD(t, y, k0), tforward, [52204.0 2000.0 36728.0 68.0 1000.0]);
yintH = Y(:,1); 
yintA = Y(:,2);
yintTH = Y(:,3);
yintTA = Y(:,4);


% Affichage des données spécifiques et de la solution de l'équation différentielle
% Display the results

figure(1)
%subplot(1,2,1);
plot(tdata, Hdata, 'r.');
hold on
plot(tdata,  yintH, 'b-');
xlabel('time in days');
ylabel('Number of cases');
legend('Data', 'Model estimation');
axis([2014 2022 0 100000]);

% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval



[k,fval] = fminunc(@moderCVD,k0);

Local minimum possible.

fminunc stopped because the size of the current step is less than
the value of the step size tolerance.

%print final values of fitted parametres
disp(k);

    0.0099    0.0499    0.0400    0.1899    0.3000    0.0100    0.5000    1.5000    2.0000    0.3000

xi=k(1);
sigma=k(2);
alpha= k(3);
gammaH=k(4) ;
gammaA=k(5) ;
deltaOC= k(6);
p1= k(7);
p2= k(8);
p3= k(9);
psi= k(10);


%Draw the data with the final ODE
[T Y] = ode45(@(t,y)(model_CVD(t,y,k)),tforward,[52204.0 2000.0 36728.0 68.0 1000.0]);
yintH = Y(:,1);
yintA = Y(:,2);
yintTH = Y(:,3);
yintTA = Y(:,4);

residuals = Hdata - yintH;% a chercher comment calcumer les residus pour plusieurs variables

%subplot(1,2,2);

figure(2)
plot(tdata,Hdata,'r.');
hold on
plot(tdata,yintH,'g-');
xlabel('Time in days');
ylabel('Number of newly infected cases with fitted parameters');
axis([2014 2022 0 100000]);

%Graphe des résidus
%subplot(1,2,2);
figure(3)
plot(tdata, residuals, 'r+');
xlabel('Time in days');
ylabel('Residuals');
title('Residuals Plot');
axis([2014 2022 min(residuals) max(residuals)]);

% % Affichage de la légende
legend('Residuals');

% 
% % Vous pouvez également ajouter une ligne horizontale à zéro pour référence
% hold on
% plot([3, 14], [0, 0], 'k--');
% 
% 
%Graphe des résidus
%subplot(1,2,2);
%figure(3)

%plot(tdata, residuals, 'g.');
%xlabel('Time in days');
%ylabel('Residuals');
%title('Residuals Plot');
%axis([2014 2022 min(residuals) max(residuals)]);

% Affichage de la légende
%legend('Residuals');

% Ajout de lignes verticales reliant chaque point de résidu à la ligne horizontale à zéro
%hold on
for i = 1:length(tdata)
    %plot([tdata(i), tdata(i)], [0, residuals(i)], 'k--');
end

% Ajout de la ligne horizontale à zéro
%plot([2014, 2022], [0, 0], 'k--');

function error_in_data = moderCVD(k)
% Données spécifiques
specific_data = [
    2014   52204 2000 36728 68;
    2015   43000  2007 41899 94;
    2016    30400 3124 51385 91;
    2017  26300 1320 55583 97;
    2018   20800 2153  59948 99;
    2019  16613 2000 64288 99;
    2020   13500 2523 65877 100;
    2021   12700 2154 67446 10;
    2022   10150 1473 68279 100
    ];

tforward = 2014:1:2022; 
%tmeasure = [1:100:901]'; 
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
Adata=specific_data(:, 3);
THdata=specific_data(:, 4);
TAdata=specific_data(:, 5);
[T, Y] = ode23s(@(t,y)(model_CVD(t,y,k)),tforward,[ 52204.0 2000.0 36728.0 68.0 1000.0]);

Hq = Y(:,1);
Aq = Y(:,2);
THq = Y(:,3);
TAq = Y(:,4);% assignts the y-coordinates of ...
%the solution at
A=(Hq - Hdata).^2;
error_in_data =sum(A);
%% 
end

function dy=model_CVD(~,y,k)
delta=0.02 ;
deltaC=0.03 ;
h=1244;
xi = k(1);
sigma = k(2);
alpha= k(3);
gammaH= k(4);
gammaA= k(5);
deltaOC= k(6);
p1= k(7);
p2= k(8);
p3= k(9);
psi=k(10);
dy = zeros(5,1);
dy(1)= h-(delta+sigma+xi+gammaH)*y(1);
dy(2)=sigma*y(1) -(delta+alpha+gammaA)*y(2);   
dy(3)=gammaH*y(1)+ psi*y(4) -(delta+deltaOC+p1*sigma+p3*xi)*y(3); 
dy(4)=gammaA*y(2)+ p1*sigma*y(3) -(psi+delta+p2*alpha)*y(4);
dy(5)=xi*y(1)+alpha*y(2)+p2*alpha*y(4)+p3*xi*y(3) -(delta+deltaC)*y(5);
end

[Khadija]

c'est un tres bon ajustement, c'etait quoi le probleme exectament?

[Khadija]

you have change Y(tmeasure(:),1); to Hq = Y(:,1); frankly I did not understand this change, since I want to assign to Hq the new values of the first variable at time tmeasure of the data!!

[Torsten]

You have to specify the output times of ode23s exactly as the times for which you have data to compare with.

Vous devez spécifier les heures de sortie de ode23s exactement comme les heures pour lesquelles vous disposez de données à comparer.

[Torsten]

you have change Y(tmeasure(:),1); to Hq = Y(:,1); frankly I did not understand this change, since I want to assign to Hq the new values of the first variable at time tmeasure of the data!!

Y(:,1) are the simulated data corresponding to "Hdata" at times "tdata". And these are exactly the values you have to use when computing the error between measured and simulated data for Hdata.

[Khadija]

Merci infiniment, Monsieur. j'ai essayé et j'ai la meme erreur avec une courbe mal ajusté que je vais vous envoyer !

j'ai effectuer le changement de temps que vous avez fait, je pense que la fonction fminunc ne fonctionne

[Torsten]

At least for the data you provided the green curve looks like a good fit to me.

I don't know which new data you tried, what changes you made to the code and why fminunc didn't succeed in this case.

[Khadija]

Bonjour Mr. Torsten; je vous remercie infiniment. J'ai reverfié mon code et ca marche tres bien. j'ai bouquiné a propos de la fonction fminunc, et j'ai trouvé qu'elle exige la non_linearité. vu que je travail sur un systeme lineaire(ode), est ce que ca ne pose pas de probleme, c'est a dire la non-linearité demandée depend des points geometriques de la representation si j'ai bien compris?

[Torsten]

First:

fminunc is a general optimizer that does not only work for nonlinear problems.

Second:

You try to estimate the parameter vector using "fminunc". So if your problem were linear, the solution of your ODE system had to be linear with respect to the parameters. This will not be the case. Consider the simple (linear) ODE y' = a*y with "a" being an unknown parameter. The solution is y(t) = C*exp(a*t), a function that is not linear in the parameter "a".

[Khadija]

je suis tres convaicu par votre exemple; si vous me recommand un lien ou un ouvrage pour ce type d'ajustement et d'optimisation, je tiens a vous informer que le programme que je vous ai envoyé est mon premier en sujet ;je vous remercie pour patience!!

[Torsten]

Star Strider's code under

https://uk.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

demonstrates how to use curve fitting for a system of ODEs.

But your code is ok - I doubt you can learn much from the link (except for "lsqcurvefit" as a different possible solver).

By the way: it would help if you used the google translator french -> english and post your comments in English:

https://translate.google.com/?hl=de

[Khadija]

Thank you sir for your help, and I will take your comment into consideration.