이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
How can i solve this problem' when i am running, this message appears: Failure in initial objective function evaluation. FMINUNC cannot continue.
조회 수: 1 (최근 30일)
이전 댓글 표시
채택된 답변
Torsten
2024년 4월 15일
편집: Torsten
2024년 4월 15일
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.
댓글 수: 19
Khadija
2024년 4월 15일
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
2024년 4월 15일
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
2024년 4월 15일
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
2024년 4월 15일
Sorry, no, I only answer questions here in the forum.
But if you have a valid licence, you could contact MATLAB support.
Torsten
2024년 4월 15일
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]);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1669221/image.png)
% 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]);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1669226/image.png)
%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');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1669231/image.png)
%
% % 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
2024년 4월 15일
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
2024년 4월 15일
편집: Torsten
2024년 4월 15일
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
2024년 4월 15일
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
Khadija
2024년 4월 16일
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
2024년 4월 16일
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
2024년 4월 16일
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
2024년 4월 16일
편집: Torsten
2024년 4월 16일
Star Strider's code under
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:
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)