現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
How can i solve this problem' when i am running, this message appears: Failure in initial objective function evaluation. FMINUNC cannot continue.
64 ビュー (過去 30 日間)
古いコメントを表示
Khadija
2024 年 4 月 15 日 8:44
I'm writing a matlab program to try to fit a mathematical model to real data.
採用された回答
Torsten
2024 年 4 月 15 日 9:58
編集済み: Torsten
2024 年 4 月 15 日 9:59
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 日 12:48
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 日 12:56
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 日 13:32
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 日 13:33
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 日 13:51
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
2024 年 4 月 15 日 14:39
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 日 14:39
編集済み: Torsten
2024 年 4 月 15 日 14:40
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
2024 年 4 月 15 日 15:24
編集済み: Torsten
2024 年 4 月 15 日 15:24
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 日 18:18
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 日 8:36
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 日 9:12
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 日 11:58
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 日 12:19
編集済み: Torsten
2024 年 4 月 16 日 12:21
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:
Khadija
2024 年 4 月 16 日 21:13
Thank you sir for your help, and I will take your comment into consideration.
その他の回答 (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!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の 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 (한국어)