//
////////////////////////////////////////////////////////////
//          probleme de la canette
// 
//   méthode de pénalisation  avec Newton plus recherche linéaire
// contrainte inégalité: volume cannette >330 cm3
///////////////////////////////////////////////////////////////
//
rmin=(165/%pi)^(1/3)
hmin=330/(%pi*rmin^2)
//
function s=surfcan(X)   
    r=X(1,:);h=X(2,:);
    s=2*%pi*r.^2+2*%pi*r.*h
endfunction
//
function s=volcan(X)   
    r=X(1,:);h=X(2,:);
    s=%pi*r.^2.*h
endfunction
//
function s=gradsurfcan(X)
    r=X(1);h=X(2);
    s=[4*%pi*r+2*%pi*h;2*%pi*r]
endfunction
//
function s=Hsurfcan(X)
    r=X(1);h=X(2);
    s=[4*%pi,2*%pi;2*%pi,0];
endfunction
//
function v=cE2(X)
    r=X(1);h=X(2);
    if(%pi*r^2*h-330<0) then 
      v=(%pi*r^2*h-330)^2
    else
       v=0
    end   
endfunction
//
function v=gradcE2(X)
    r=X(1);h=X(2);
    if(%pi*r^2*h-330<0) then 
       v=[4*%pi^2*r^3*h^2-4*330*%pi*r*h;2*%pi^2*r^4*h-2*%pi*330*r^2]
    else
       v=[0;0]
    end 
endfunction
//
function v=HcE2(X)
    r=X(1);h=X(2);
    if(%pi*r^2*h-330<0) then 
       v=[12*%pi^2*r^2*h^2-4*330*%pi*h,8*%pi^2*r^3*h-4*330*%pi*r;8*%pi^2*r^3*h-4*330*%pi*r,2*%pi^2*r^4]
   else
       v=zeros(2,2)
   end 
endfunction
//
function Xt=NewtonCont(J,gradJ,HJ,cEcarre,gradcEcarre,HcEcarre,X0,N,sigma)
   XN=X0;Xt=XN;
   for i=1:N
      alpha=alphain
      Jsig=J(XN)+sigma*cEcarre(XN)
      Jsig_grad=gradJ(XN)+sigma*gradcEcarre(XN) 
      Jsig_H=HJ(XN)+sigma*HcEcarre(XN)
      d=-inv(Jsig_H)*Jsig_grad  
      while(J(XN+alpha*d)+sigma*cEcarre(XN+alpha*d)>Jsig+bet*alpha*d'*(Jsig_grad))
      alpha=tau*alpha
     end
     XN=XN+alpha*d
   Xt=[Xt,XN];          
   end      
endfunction
//
X0=[10;6];  // point initial
sigma=100;  // coefficient de pénalisation
bet=0.01    // constante d'Armijo
alphain=1;  // pas initial
tau=0.3;    // constante de backtracking
N=100;       // nombre d'itérations
Xt=NewtonCont(surfcan,gradsurfcan,Hsurfcan,cE2,gradcE2,HcE2,X0,N,sigma)
disp(Xt)
//
//tracé
//
figure(2)
r=linspace(2,4,100)
h=linspace(6,9,100)
z=zeros(100,100);hvol=zeros(100,1)
for i=1:100
    hvol(i)=330/(%pi*r(i)^2);
    for j=1:100
      z(i,j)=surfcan([r(i);h(j)]) 
    end   
end
clf()
contour(r,h,z,20)
plot2d([rmin,rmin],[hmin,hmin],-1)
plot2d(r,hvol,1,rect=[2,6,4,9])
//
plot(Xt(1,:),Xt(1,:),'b-')
//
//
figure(3)
clf()
plot(0:N,Xt(1,:),'r-')
plot(0:N,rmin*ones(N+1,1),'b-')
plot(0:N,Xt(2,:),'r-')
plot(0:N,hmin*ones(N+1,1),'b-')
xtitle('évolution du rayon et de la hauteur (bleu: valeur optimale')
//
//
figure(4)
clf()
plot(0:N,surfcan(Xt),'r-')
plot(0:N,volcan(Xt),'r-')
plot(0:N,330*ones(N+1,1),'b-')
xtitle('évolution du volume et de la surface (bleu: volume 33cl)')
//
//
// on peut constater que la convergence vers le point (rmin,hmin) est très sensible en
// fonction du point de départ

