gradient.m

optimisation中的梯度法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         TP1 Optimisation et commande                                     %
%          Recherche d'optimum--------Methode Gradient + Moindres Carrees                  %
%          Programmation par   Nicolas MARTIN  et  Ning QU                                 %
%          Date: 11/10/2007                                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

%courbe de donnees
load('granul');

%valeur initiale
%estimation
a1=10;a2=8;
mu1=450;mu2=695;sig1=69;sig2=91;
teta=[mu1;mu2;sig1;sig2];
epsilon=0.001;
k=2;
G=[2;2;2;2];
i=0;

%tableau d evolution
evo_teta=[];
evo_T=[];

while norm(G)>epsilon
    i=i+1; %calcul le nombre d iteration
    %construction du signal
    ym=a1*exp(-(x-mu1).^2/sig1^2)+ a2*exp(-(x-mu2).^2/sig2^2);
    e=y-ym; %erreur
    %Gradient
    G1=-a1/sig1^2*e'*((x-mu1).*exp(-(x-mu1).^2/sig1^2));
    G2=-a2/sig2^2*e'*((x-mu2).*exp(-(x-mu2).^2/sig2^2));
    G3=-a1/sig1^3*e'*((x-mu1).^2.*exp(-(x-mu1).^2/sig1^2));
    G4=-a2/sig2^3*e'*((x-mu2).^2.*exp(-(x-mu2).^2/sig2^2));

    G=[G1;G2;G3;G4]; %vecteur
    teta=teta - k*G; % calcul des nouvelles valeur
    evo_teta(i,:)=teta';

    mu1=teta(1);
    mu2=teta(2);
    sig1=teta(3);
    sig2=teta(4);

    %methode moindre carree
    E1=exp(-(x-teta(1)).^2/teta(3)^2);
    E2=exp(-(x-teta(2)).^2/teta(4)^2);
    X=[E1 E2];
    T= inv(X'*X)*X'*y;
    a1=T(1);
    a2=T(2);
    evo_T(i,:)=T';
end

teta,i,T

%signal final
ym=a1*exp(-(x-mu1).^2/sig1^2)+ a2*exp(-(x-mu2).^2/sig2^2);
%signal bruite estime
ymb=a1*exp(-(x-mu1).^2/sig1^2)+ a2*exp(-(x-mu2).^2/sig2^2)+rand(size(x))/20;
%courbe bruite
plot(x,y,'*')
hold on
%courbe estime
plot(x,ym,'r')
% signal bruite estime
hold on
plot(x,ymb,'g')