///////////////////////////////////////
//
// probleme de la canette:
// résolution par pénalisation externe ou interne
// et algorithme de gradient
//
////////////////////////////////////////
function z=f(X)  // fonction cout non penalisee
    d=X(1);h=X(2);
    z=%pi*d^2/4+%pi*d*h;
endfunction
//
function z=gradf(X)  // gradient fonction cout non penalisee
    d=X(1);h=X(2);
    z=[%pi*d/2+%pi*h;%pi*d];
endfunction
//
function z=g(X)   //contrainte g(X)<0
    d=X(1);h=X(2);
    z=-%pi*d^2*h/4+300;
endfunction
//
function z=gradg(X)  // gradient contrainte
          d=X(1);h=X(2);
          z=-[%pi*d*h/2;%pi*d^2/4]      
endfunction  
// 
function z=lagrangien(X,lamb)  // lagrangien
    z=f(X)+lamb*g(X);
endfunction   
//
//
//////////////////////////////////////////////////////
// solution exacte
//
d0=(2400/%pi)^(1/3);h0=1200/(%pi*d0^2);
lambexact=gradf([d0;h0])./-gradg([d0;h0]);
//
///////////////////////////////////////////////////////
//
X0=evstr(x_dialog('point initial','[9;9]'));
lamb0=evstr(x_dialog('lambda initial','2'));
//
///////////////////////
alpha=0.001;
rho=0.1;
//
Niter=4000;
X=X0;
lamb=lamb0;
Xtot=[X;lamb];
////////////////////////
 for i=1:Niter
  dgrad=gradf(X)+lamb*gradg(X);
  d=-dgrad;      
  X=X+alpha*d;
 lamb=max(0,lamb+rho*g(X));
 Xtot=[Xtot,[X;lamb]];     
end 
//
disp('valeur exacte:');disp([d0;h0])
disp('surface minimale canette:');disp(f([d0,h0]))

disp('valeur optimale obtenue:');disp(X)
//
disp('surface canette:');disp(f(X))
disp('volume canette:');disp(300-g(X))
//
//
///////////////// tracés ////////////////////
//
Nt=200;
x=linspace(2,14,Nt);
y=linspace(2,14,Nt);
z=zeros(Nt,Nt);
for i=1:Nt
    for j=1:Nt
       z(i,j)=f([x(i),y(j)]) 
    end
end
figure(1)
clf()
contour(x,y,z,20)
dplot=linspace(5,12,Nt);
hplot=1200/%pi*ones(dplot)./(dplot.^2);
plot2d(dplot,hplot)
plot2d(Xtot(1,:),Xtot(2,:))
plot2d(d0,h0,-1)
//
///////////////////////////////////////////////////////



