commit cc759343611db75235926b94a00142605563d936 Author: Laureηt Date: Sat Jun 10 20:58:09 2023 +0200 init diff --git a/Algo_Gauss_Newton.m b/Algo_Gauss_Newton.m new file mode 100755 index 0000000..9c5d8f8 --- /dev/null +++ b/Algo_Gauss_Newton.m @@ -0,0 +1,87 @@ +function [beta, norm_grad_f_beta, f_beta, norm_delta, nb_it, exitflag] ... + = Algo_Gauss_Newton(residu, J_residu, beta0, option) + +%***************************************************************** +% Fichier ~gergaud/ENS/Optim1a/TP-optim-20-21/TP-ref/GN_ref.m * +% Novembre 2020 * +% Université de Toulouse, INP-ENSEEIHT * +%***************************************************************** +% +% GN resout par l'algorithme de Gauss-Newton les problemes aux moindres carres +% Min 0.5||r(beta)||^2 +% beta \in \IR^p +% +% Paramètres en entrés +% -------------------- +% residu : fonction qui code les résidus +% r : \IR^p --> \IR^n +% J_residu : fonction qui code la matrice jacobienne +% Jr : \IR^p --> real(n,p) +% beta0 : point de départ +% real(p) +% option(1) : Tol_abs, tolérance absolue +% real +% option(2) : Tol_rel, tolérance relative +% real +% option(3) : n_itmax, nombre d'itérations maximum +% integer +% +% Paramètres en sortie +% -------------------- +% beta : beta +% real(p) +% norm_grad_f_beta : ||gradient f(beta)|| +% real +% f_beta : f(beta) +% real +% r_beta : r(beta) +% real(n) +% norm_delta : ||delta|| +% real +% nbit : nombre d'itérations +% integer +% exitflag : indicateur de sortie +% integer entre 1 et 4 +% exitflag = 1 : ||gradient f(beta)|| < max(Tol_rel||gradient f(beta0)||,Tol_abs) +% exitflag = 2 : |f(beta^{k+1})-f(beta^k)| < max(Tol_rel|f(beta^k)|,Tol_abs) +% exitflag = 3 : ||delta)|| < max(Tol_rel delta^k),Tol_abs) +% exitflag = 4 : nombre maximum d'itérations atteint +% +% --------------------------------------------------------------------------------- + + norm_grad_f_beta0 = norm(J_residu(beta0).'*residu(beta0)); + + nb_it = 0; + beta = beta0; + + continuer = 1; + + while continuer + + nb_it = nb_it + 1; + + beta_last = beta; + beta = beta - (J_residu(beta).'*J_residu(beta))\(J_residu(beta).'*residu(beta)); + + norm_grad_f_beta = norm(J_residu(beta).'*residu(beta)); + f_beta = norm(residu(beta))/2; + delta = beta - beta_last; + norm_delta = norm(delta); + + if norm(J_residu(beta).'*residu(beta)) <= max(option(2)*norm_grad_f_beta0, option(1)) + exitflag = 1; + continuer = 0; + elseif abs(norm(residu(beta))/2 - norm(residu(beta_last))/2) <= max(option(2)*norm(residu(beta))/2, option(1)) + exitflag = 2; + continuer = 0; + elseif norm_delta <= max(option(2)*norm(beta), option(1)) + exitflag = 3; + continuer = 0; + elseif nb_it >= option(3) + exitflag = 4; + continuer = 0; + end + + end + +end \ No newline at end of file diff --git a/Algo_Newton.m b/Algo_Newton.m new file mode 100755 index 0000000..565389a --- /dev/null +++ b/Algo_Newton.m @@ -0,0 +1,91 @@ +function [beta, norm_grad_f_beta, f_beta, norm_delta, nb_it, exitflag] ... + = Algo_Newton(Hess_f, beta0, option) +%************************************************************ +% Fichier ~gergaud/ENS/Optim1a/TP-optim-20-21/Newton_ref.m * +% Novembre 2020 * +% Université de Toulouse, INP-ENSEEIHT * +%************************************************************ +% +% Newton résout par l'algorithme de Newton les problemes aux moindres carres +% Min 0.5||r(beta)||^2 +% beta \in R^p +% +% Parametres en entrees +% -------------------- +% Hess_f_C14 : fonction qui code la hessiennne de f +% Hess_f_C14 : R^p --> matrice (p,p) +% (la fonction retourne aussi le residu et la jacobienne) +% beta0 : point de départ +% real(p) +% option(1) : Tol_abs, tolérance absolue +% real +% option(2) : Tol_rel, tolérance relative +% real +% option(3) : nitimax, nombre d'itérations maximum +% integer +% +% Parametres en sortie +% -------------------- +% beta : beta +% real(p) +% norm_gradf_beta : ||gradient f(beta)|| +% real +% f_beta : f(beta) +% real +% res : r(beta) +% real(n) +% norm_delta : ||delta|| +% real +% nbit : nombre d'itérations +% integer +% exitflag : indicateur de sortie +% integer entre 1 et 4 +% exitflag = 1 : ||gradient f(beta)|| < max(Tol_rel||gradient f(beta0)||,Tol_abs) +% exitflag = 2 : |f(beta^{k+1})-f(beta^k)| < max(Tol_rel|f(beta^k)|,Tol_abs) +% exitflag = 3 : ||delta)|| < max(Tol_rel delta^k),Tol_abs) +% exitflag = 4 : nombre maximum d'itérations atteint +% +% --------------------------------------------------------------------------------- + +% TO DO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + [H_f0, res0, J_res0] = Hess_f(beta0); + norm_grad_f_beta0 = norm(J_res0.'*res0); + + nb_it = 0; + beta = beta0; + + continuer = 1; + + while continuer + + nb_it = nb_it + 1; + + beta_last = beta; + [H_f_last, res_last, J_res_last] = Hess_f(beta_last); + + beta = beta - H_f_last\(J_res_last.'*res_last); + [H_f, res, J_res] = Hess_f(beta); + + norm_grad_f_beta = norm(J_res.'*res); + f_beta = norm(res)/2; + delta = beta - beta_last; + norm_delta = norm(delta); + + if norm_grad_f_beta <= max(option(2)*norm_grad_f_beta0, option(1)) + exitflag = 1; + continuer = 0; + elseif abs(norm(res)/2 - norm(res_last)/2) <= max(option(2)*norm(res)/2, option(1)) + exitflag = 2; + continuer = 0; + elseif norm_delta <= max(option(2)*norm(beta), option(1)) + exitflag = 3; + continuer = 0; + elseif nb_it >= option(3) + exitflag = 4; + continuer = 0; + end + + end + +end diff --git a/C14_figures.png b/C14_figures.png new file mode 100755 index 0000000..dad950e Binary files /dev/null and b/C14_figures.png differ diff --git a/C14_results.txt b/C14_results.txt new file mode 100755 index 0000000..04cd00a --- /dev/null +++ b/C14_results.txt @@ -0,0 +1,46 @@ +Algorithme de Gauss-Newton +-------------------------------------------------------------------------------------------- + nb_iter A0 lambda ||f'(beta)|| f(beta) ||delta|| exitflag +-------------------------------------------------------------------------------------------- + 0 10 0.0001 4.6322e+05 48.07 + + 1 15.022 0.00010633 15913 0.22921 5.0219 4 + + 2 15.025 0.00010433 5.9024 0.2105 0.0032964 4 + + 3 15.025 0.00010432 0.39911 0.2105 0.00068766 4 + + 4 15.024 0.00010432 0.004769 0.2105 4.9165e-06 2 + + 4 15.024 0.00010432 0.004769 0.2105 4.9165e-06 2 + + 4 15.024 0.00010432 0.004769 0.2105 4.9165e-06 2 + + 4 15.024 0.00010432 0.004769 0.2105 4.9165e-06 2 + + 4 15.024 0.00010432 0.004769 0.2105 4.9165e-06 2 + +-------------------------------------------------------------------------------------------- +Algorithme de Newton +-------------------------------------------------------------------------------------------- + nb_iter A0 lambda ||f'(beta)|| f(beta) ||delta|| exitflag +-------------------------------------------------------------------------------------------- + 0 10 0.0001 4.6322e+05 48.07 + + 1 13.954 5.5076e-05 3.0226e+05 1.6021 3.9536 4 + + 2 14.774 0.00010067 6753.7 0.26139 0.82072 4 + + 3 14.987 0.00010359 450.57 0.2117 0.21288 4 + + 4 15.019 0.00010424 80.175 0.21053 0.031965 4 + + 5 15.024 0.00010431 5.8177 0.2105 0.004593 4 + + 6 15.024 0.00010432 1.184 0.2105 0.0006653 4 + + 7 15.024 0.00010432 0.066427 0.2105 9.519e-05 2 + + 7 15.024 0.00010432 0.066427 0.2105 9.519e-05 2 + +-------------------------------------------------------------------------------------------- diff --git a/Hess_f_C14.m b/Hess_f_C14.m new file mode 100755 index 0000000..b0bf0fb --- /dev/null +++ b/Hess_f_C14.m @@ -0,0 +1,32 @@ +function [H_f, res, J_res] = Hess_f_C14(beta, donnees, residu, J_residu) +% +% Paramètres en entrés +% -------------------- +% beta : vecteur des paramètres +% real(p) +% donnees : Données +% real(n,2) +% residu : fonction qui code les résidus +% res_beta = residus(beta) +% J_residu : fonction qui code la matrice jacobienne +% J_res_beta = J_residu(beta); +% +% Paramètres en sortie +% -------------------- +% H_f : Matrice hessienne +% real(p,p) +% res : vecteur des résidus +% real(n) +% J_res : Matrice jacobienne des résiduis +% real(n,p) +% + + res = residu(beta); + J_res = J_residu(beta); + % cf formule du cours + S = [ sum(0*donnees(:,1)) , sum(donnees(:,1).*exp(-beta(2)*donnees(:, 1))); + sum(donnees(:, 1).*exp(-beta(2)*donnees(:, 1))), sum(-beta(1)*(donnees(:, 1).^2).*exp(-beta(2)*donnees(:, 1))) ]; + H_f = S + J_res.'*J_res; + +end + diff --git a/J_residu_C14.m b/J_residu_C14.m new file mode 100755 index 0000000..25a895a --- /dev/null +++ b/J_residu_C14.m @@ -0,0 +1,28 @@ +%-------------------------------------------------------------------------% % +% Fonction de calcul de la Jacobienne du residu de la fonction de % +% desintegration radioactive du carbone 14 % +%-------------------------------------------------------------------------% + +function J_res = J_residu_C14(beta, donnees) +% +% Paramètres en entrés +% -------------------- +% beta : vecteur des paramètres +% real(p) +% donnees : Données +% real(n,2) +% +% Paramètres en sortie +% -------------------- +% J_res : Matrice jacobienne des résidus +% real(n,p) +% + + % res = donnees(:, 2) - beta(1)*exp(-beta(2)*donnees(:, 1)); + + v1 = -exp(-beta(2)*donnees(:, 1)); + v2 = beta(1)*donnees(:, 1).*exp( -beta(2)*donnees(:, 1) ); + + J_res = [ v1, v2 ]; + +end \ No newline at end of file diff --git a/Modelisation_C14.m b/Modelisation_C14.m new file mode 100755 index 0000000..94c4435 --- /dev/null +++ b/Modelisation_C14.m @@ -0,0 +1,206 @@ +%-------------------------------------------------------------------------% +% 1SN - TP Optimisation % +% INP Toulouse - ENSEEIHT % +% Novembre 2020 % +% % +% Ce fichier contient le programme principal permettant l'estimation % +% des parametres de la fonction de desintegration radioactive du % +% carbone 14 par une approche des moindres carres. % % +% Modele : A(t)= A0*exp(-lambda*t) % +% Les algorithmes utilises pour la minimisation sont : % +% - l'algorithme de Gauss-Newton % +% - l'algorithme de Newton. % +%-------------------------------------------------------------------------% + +clear +close all +clc +format shortG +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); + +% Initialisation + +% Donnees +if exist('C14_results.txt','file') + delete('C14_results.txt'); +end +diary C14_results.txt +Ti = [ 500; 1000; 2000; 3000; 4000; 5000; 6300]; +Ai = [14.5; 13.5; 12.0; 10.8; 9.9; 8.9; 8.0]; +Donnees = [Ti, Ai]; + +% Estimation a priori des parametres du modele : beta0 = [A0, lambda] +beta0 = [10; 0.0001]; % Newton, Gauus-Newton, fminunc et leastsq converge +% beta0 = [15; 0.001]; % Newton, Gauss-Newton, fminunc et leastsq divergent +% beta0 = [15; 0.0005]; % Newton diverge, Gauss-Newton, fminunc et leastsq convergent +% beta0 = [10; 0.0005]; % Gauss-Newton converge + +%% Calcul et affichage du modele initial ---------------------------------- +xmin = 9; xmax = 20; +x = linspace(xmin, xmax, 100); + +ymin = -0.0001; ymax = 0.0005; +y = linspace(ymin, ymax, 100); + +[A0_plot, lambda_plot] = meshgrid(x, y); +[m , n] = size(A0_plot); +f_plot = zeros(m, n); + +for i = 1:m + for j = 1:n + res_plot = residu_C14([A0_plot(i,j); lambda_plot(i,j)], Donnees); + f_plot(i,j) = res_plot.'*res_plot /2; + end +end + +figure('Position',[0.1*L,0.1*H,0.8*L,0.8*H]); +subplot(1,3,1) % Fonction f en 3D + mesh(A0_plot, lambda_plot, f_plot) + axis('square') + title('Representation de la fonction f des moindres carres') + xlabel('A_0'); + ylabel('\lambda') + zlabel('f(A_0,\lambda)') + +T = linspace(0,6500,100); +A = beta0(1)*exp(-beta0(2)*T); +txt_legend{1} = 'donnees'; +txt_legend{2} = 'depart'; + +subplot(2, 3, 2) % Donnees + modele initial pour Gauss-Newton + title('Desintegration radioactive du carbone 14 (Gauss-Newton)') + axis([0 max(T) 0 18]) + xlabel('dur�e T') + ylabel('radioactivit� A') + hold on + grid on + plot(Ti, Ai, 'ok') + plot(T, A) + legend(txt_legend{1:2}, 'Location', 'SouthWest') + +subplot(2, 3, 5) % Courbes de niveaux de f pour Gauss-Newton + title('Recherche des parametres (Gauss-Newton)') + xlabel('A_0') + ylabel('\lambda') + hold on + contour(A0_plot, lambda_plot, f_plot, 100); + plot(beta0(1),beta0(2),'ok') + text(beta0(1),beta0(2),' depart \beta^{(0)}') + +subplot(2, 3, 3) % Donnees + modele initial pour Newton + title('Desintegration radioactive du carbone 14 (Newton)') + axis([0 max(T) 0 18]) + xlabel('dur�e T') + ylabel('radioactivit� A') + hold on + grid on + plot(Ti, Ai, 'ok') + plot(T, A) + legend(txt_legend{1:2}, 'Location', 'SouthWest') + +subplot(2, 3, 6) % Courbes de niveaux de f pour Newton + title('Recherche des parametres (Newton)') + xlabel('A_0') + ylabel('\lambda') + hold on + contour(A0_plot, lambda_plot, f_plot, 100); + plot(beta0(1),beta0(2),'ok') + text(beta0(1),beta0(2),' depart \beta^{(0)}') + +pause(0.5) + +%% Algorithmes +% Choix du nombre d'iteration maximal + +nb_iterations_max = 8; +txt_legend = cell(1,nb_iterations_max+2); +txt_legend{1} = 'donnees'; +txt_legend{2} = 'depart'; + + +%% Gauss-Newton +% ------------- +% Initialisation de l'affichage +disp('Algorithme de Gauss-Newton') +disp('--------------------------------------------------------------------------------------------') +disp(' nb_iter A0 lambda ||f''(beta)|| f(beta) ||delta|| exitflag ') +disp('--------------------------------------------------------------------------------------------') + +% Calcul et affichage des valeurs initiales +res_beta = residu_C14(beta0, Donnees); +f_beta = 0.5*(res_beta.')*res_beta; +J_res_beta = J_residu_C14(beta0, Donnees); +norm_grad_f_beta = norm((J_res_beta.')*res_beta); +disp([0 beta0(1) beta0(2) norm_grad_f_beta f_beta]); + +options = [sqrt(eps) 1.e-12 0]; + +for i = 1:nb_iterations_max + txt_legend{i+2} = ['iteration ' num2str(i)]; + options(3) = i; + [beta, norm_grad_f_beta, f_beta, norm_delta, k, exitflag] = ... + Algo_Gauss_Newton(@(beta) residu_C14(beta, Donnees), ... + @(beta) J_residu_C14(beta, Donnees), ... + beta0, options); + disp([k beta(1) beta(2) norm_grad_f_beta f_beta norm_delta exitflag]) + % Visualisation + A = beta(1)*exp(-beta(2)*T); + subplot(2, 3, 2) + plot(T,A) + legend(txt_legend{1:i+2}) + % eval(['print -depsc fig_GN_courbe' int2str(i) '_C14']) + subplot(2, 3, 5) + plot(beta(1),beta(2),'ok') + text(beta(1),beta(2),[' \beta^{(' num2str(i) ')}']) + + pause(0.5) +end + +% Affichage des itérés de beta et sauvegarde des graphique +disp('--------------------------------------------------------------------------------------------') + +%% Newton +% ------- +% Initialisation de l'affichage +disp('Algorithme de Newton') +disp('--------------------------------------------------------------------------------------------') +disp(' nb_iter A0 lambda ||f''(beta)|| f(beta) ||delta|| exitflag ') +disp('--------------------------------------------------------------------------------------------') + +% Calcul et affichage des valeurs initiales +res_beta = residu_C14(beta0, Donnees); +f_beta = 0.5*(res_beta.')*res_beta; +J_res_beta = J_residu_C14(beta0, Donnees); +norm_grad_f_beta = norm((J_res_beta.')*res_beta); +disp([0 beta0(1) beta0(2) norm_grad_f_beta f_beta]); + +options = [sqrt(eps) 1.e-12 0]; + +for i = 1:nb_iterations_max + options(3) = i; + % Algorithme de Newton + [beta, norm_grad_f_beta, f_beta, norm_delta, k, exitflag] = ... + Algo_Newton(@(beta) Hess_f_C14(beta, Donnees, ... + @(beta) residu_C14(beta, Donnees), ... + @(beta) J_residu_C14(beta, Donnees)), ... + beta0, options); + % Affichage des valeurs + disp([k beta(1) beta(2) norm_grad_f_beta f_beta norm_delta exitflag]) + % Visualisation + A = beta(1)*exp(-beta(2)*T); + subplot(2, 3, 3) + plot(T,A) + legend(txt_legend{1:i+2}) + subplot(2, 3, 6) + plot(beta(1),beta(2),'ok') + text(beta(1),beta(2),[' \beta^{(' num2str(i) ')}']) + + pause(0.5) +end + +% Affichage des iteres de beta et sauvegarde des courbes +disp('--------------------------------------------------------------------------------------------') +print('C14_figures','-dpng') +diary diff --git a/TP_Opti1A.pdf b/TP_Opti1A.pdf new file mode 100755 index 0000000..20814f4 Binary files /dev/null and b/TP_Opti1A.pdf differ diff --git a/residu_C14.m b/residu_C14.m new file mode 100755 index 0000000..63ad9db --- /dev/null +++ b/residu_C14.m @@ -0,0 +1,23 @@ +%%------------------------------------------------------------------------% % +% Fonction de calcul du residu de la fonction de desintegration % +% radioactive du carbone 14 % +%-------------------------------------------------------------------------% + +function res = residu_C14(beta, donnees) +% +% Paramètres en entrés +% -------------------- +% beta : vecteur des paramètres +% real(p) +% donnees : Données +% real(n,2) +% +% Paramètres en sortie +% -------------------- +% res : vecteur des résidus +% real(n) +% + + res = donnees(:, 2) - beta(1)*exp(-beta(2)*donnees(:, 1)); + +end