diff --git a/Optimisation/TP/TP-C14/Algo_Gauss_Newton.m b/Optimisation/TP/TP-C14/Algo_Gauss_Newton.m new file mode 100644 index 0000000000000000000000000000000000000000..2573a0dadd4bc908ff2325ceadda51bd79d7ec99 --- /dev/null +++ b/Optimisation/TP/TP-C14/Algo_Gauss_Newton.m @@ -0,0 +1,60 @@ +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_gradf_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 +% +% --------------------------------------------------------------------------------- + +% TO DO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + beta = [0;0]; + norm_grad_f_beta = 0; + f_beta = 0; + norm_delta = 0; + nb_it = 0; + exitflag = 0; + +end diff --git a/Optimisation/TP/TP-C14/Algo_Newton.m b/Optimisation/TP/TP-C14/Algo_Newton.m new file mode 100644 index 0000000000000000000000000000000000000000..e7ab662969535225f16c01e203790a1e401e4145 --- /dev/null +++ b/Optimisation/TP/TP-C14/Algo_Newton.m @@ -0,0 +1,59 @@ +function [beta, norm_grad_f_beta, f_beta, norm_delta, nb_it, exitflag] ... + = Algo_Newton(Hess_f_C14, 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + beta = [0;0]; + norm_grad_f_beta = 0; + f_beta = 0; + norm_delta = 0; + nb_it = 0; + exitflag = 0; + +end diff --git a/Optimisation/TP/TP-C14/C14_figures.png b/Optimisation/TP/TP-C14/C14_figures.png new file mode 100644 index 0000000000000000000000000000000000000000..906349d48b1a0d704e8413dbf1ef56970b2760ca Binary files /dev/null and b/Optimisation/TP/TP-C14/C14_figures.png differ diff --git a/Optimisation/TP/TP-C14/C14_results.txt b/Optimisation/TP/TP-C14/C14_results.txt new file mode 100644 index 0000000000000000000000000000000000000000..0109f4f050726bf98c20aa8223760edc262d783e --- /dev/null +++ b/Optimisation/TP/TP-C14/C14_results.txt @@ -0,0 +1,65 @@ +[Warning: Cannot set Position while WindowStyle is 'docked'] +[> In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('Modelisation_C14', '/Users/gergaud/git-ENS/Optimisation-et-EDP-1A/enseignants/Optimisation/TP/TP-C14/src_etudiants/Modelisation_C14.m', 59)" style="font-weight:bold">Modelisation_C14</a> (<a href="matlab: opentoline('/Users/gergaud/git-ENS/Optimisation-et-EDP-1A/enseignants/Optimisation/TP/TP-C14/src_etudiants/Modelisation_C14.m',59,0)">line 59</a>)] +[Warning: Contour not rendered for constant ZData] +[> In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('contour', '/Applications/MATLAB_R2022b.app/toolbox/matlab/specgraph/contour.m', 51)" style="font-weight:bold">contour</a> (<a href="matlab: opentoline('/Applications/MATLAB_R2022b.app/toolbox/matlab/specgraph/contour.m',51,0)">line 51</a>) +In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('Modelisation_C14', '/Users/gergaud/git-ENS/Optimisation-et-EDP-1A/enseignants/Optimisation/TP/TP-C14/src_etudiants/Modelisation_C14.m', 89)" style="font-weight:bold">Modelisation_C14</a> (<a href="matlab: opentoline('/Users/gergaud/git-ENS/Optimisation-et-EDP-1A/enseignants/Optimisation/TP/TP-C14/src_etudiants/Modelisation_C14.m',89,0)">line 89</a>)] +[Warning: Contour not rendered for constant ZData] +[> In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('contour', '/Applications/MATLAB_R2022b.app/toolbox/matlab/specgraph/contour.m', 51)" style="font-weight:bold">contour</a> (<a href="matlab: opentoline('/Applications/MATLAB_R2022b.app/toolbox/matlab/specgraph/contour.m',51,0)">line 51</a>) +In <a href="matlab:matlab.internal.language.introspective.errorDocCallback('Modelisation_C14', '/Users/gergaud/git-ENS/Optimisation-et-EDP-1A/enseignants/Optimisation/TP/TP-C14/src_etudiants/Modelisation_C14.m', 109)" style="font-weight:bold">Modelisation_C14</a> (<a href="matlab: opentoline('/Users/gergaud/git-ENS/Optimisation-et-EDP-1A/enseignants/Optimisation/TP/TP-C14/src_etudiants/Modelisation_C14.m',109,0)">line 109</a>)] +Algorithme de Gauss-Newton +-------------------------- +residu_C14(beta0, Donnees) + 0 + +J_residu_C14(beta0, Donnees) + 0 + +-------------------------------------------------------------------------------------------- + nb_iter A0 lambda ||f'(beta)|| f(beta) ||delta|| exitflag +-------------------------------------------------------------------------------------------- + 0 10 0.0001 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + +-------------------------------------------------------------------------------------------- +Algorithme de Newton +-------------------- +Hessienne f(beta^{(0)}) + 0 + +-------------------------------------------------------------------------------------------- + nb_iter A0 lambda ||f'(beta)|| f(beta) ||delta|| exitflag +-------------------------------------------------------------------------------------------- + 0 10 0.0001 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + + 0 0 0 0 0 0 0 + +-------------------------------------------------------------------------------------------- diff --git a/Optimisation/TP/TP-C14/Modelisation_C14.m b/Optimisation/TP/TP-C14/Modelisation_C14.m new file mode 100644 index 0000000000000000000000000000000000000000..d1844bbfed5118bec37683b49145de86eea864ee --- /dev/null +++ b/Optimisation/TP/TP-C14/Modelisation_C14.m @@ -0,0 +1,316 @@ +%-------------------------------------------------------------------------% +% 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); + % TO DO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + f_plot(i,j) = 0.5*(res_plot.')*res_plot; + 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('residu_C14(beta0, Donnees)') +disp(residu_C14(beta0, Donnees)) +disp('J_residu_C14(beta0, Donnees)') +disp(J_residu_C14(beta0, Donnees)) +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} = ['iterat� ' 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('--------------------') +% Initialisation de l'affichage +[H, ~, ~] = Hess_f_C14(beta0, Donnees, ... + @(beta) residu_C14(beta, Donnees), ... + @(beta) J_residu_C14(beta, Donnees)); +disp('Hessienne f(beta^{(0)})') +disp(H) + +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 + +%% Fonctions a coder : residu, Jacobienne, Hessienne + +%%------------------------------------------------------------------------% % +% 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) +% + % TO DO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + res = 0; + +end + +%-------------------------------------------------------------------------% % +% 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) +% + % TO DO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + J_res = 0; + +end + +%-------------------------------------------------------------------------% +% Fonction de calcul de la Hessienne des moindres carres de la fonction % +% de desintegration radioactive du carbone 14 % +% (residu et J_residu sont des fonctions en entree) % +%-------------------------------------------------------------------------% + +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) +% + % TO DO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + res = 0; + J_res = 0; + H_f = 0; + +end + + + + + + + + + + + diff --git a/Optimisation/TP/TP-C14/TP_Opti1A.pdf b/Optimisation/TP/TP-C14/TP_Opti1A.pdf new file mode 100644 index 0000000000000000000000000000000000000000..7b90ba960361a63645362178d2ffffb0b8a2403b Binary files /dev/null and b/Optimisation/TP/TP-C14/TP_Opti1A.pdf differ