//
// optimisation d'une fonction J avec une contrainte de type égalité
//
////////////////////////////////////////////////
// Uzawa
///////////////////////////////////////////////
function [X,lamb]=uzawa(J, nablaJ,cE,nablacE,X0,lamb0,alpha1,alpha2,N)
    X=X0;lamb=lamb0;
    for i=1:N
        X=X-alpha1*(nablaJ(X)+lamb*nablacE(X))
        lamb=lamb+alpha2*cE(X)
    end
endfunction
////////////////////////////////////////////////////
// exemple: fonctionnelle quadratique, contrainte affine
/////////////////////////////////////////////////////
//
n=10;
A=toeplitz([2,-1,zeros(1,n-2)]);
b=ones(n,1)
//
function y=quad(x)
    y=1/2*(A*x)'*x-b'*x
endfunction
function y=nablaquad(x)
    y=A*x-b
endfunction
function y=cEex(x)
    y=sum(x)-1
endfunction        
function y=nablacEex(x)
    y=ones(n,1)
endfunction  
////////////////// 
alpha1=0.01;
alpha2=0.01;
X0=zeros(n,1);
lamb0=3;
N=10000;
[X,lamb]=uzawa(quad,nablaquad,cEex,nablacEex,X0,lamb0,alpha1,alpha2,N)
disp('solution Uzawa ');
disp(X)
disp('somme de X (égale à 1 en théorie):')
disp(sum(X))
/////////////////////////////////////////////////////////////////
// penalisation
/////////////////////////////////////////////////////////
function X=penal(J, nablaJ,cE,nablacE,X0,rho,N)
    X=X0;alphain=1;bet=0.001;tau=0.3;
    for i=1:N
     alpha=alphain;d=-(nablaJ(X)+2*rho*cE(X)*nablacE(X));
     while(J(X+alpha*d)+rho*cE(X+alpha*d)^2>J(X)+rho*cE(X)^2+bet*alpha*d'*d)
      alpha=tau*alpha
     end
     X=X+alpha*d
    end
endfunction
///////////////////////////////////////////////////////
// exemple: fonctionnelle quadratique, contrainte affine
/////////////////////////////////////////////////////
//
rho=250;
N=4000;
X=penal(quad,nablaquad,cEex,nablacEex,X0,rho,N)
disp('solution penalisation');
disp(X)
disp('somme de X(égale à 1 en théorie):')
disp(sum(X))
