///////////////////////////////////////
//
// 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=c(X)   //contrainte
    d=X(1);h=X(2);
    z=%pi*d^2*h/4-300;
endfunction
//////////////////////////////////////////////////
function z=penal(X)  // fonction de penalisation externe
    d=X(1);h=X(2);
    z=max(0,-c([d,h]));
endfunction
//
function z=penal2(X)  // fonction de penalisation interne
    d=X(1);h=X(2);
        if (c(X)>0) then
            z=log(c(X));
        else
            z=0;
        end        
endfunction
//
function z=gradpenal(X)  // gradient fonction de penalisation externe
    d=X(1);h=X(2);
    if (c(X)>0) then
      z=[0;0]
    else
      z=-[%pi*d*h/2;%pi*d^2/4]      
    end
endfunction     
//
function z=gradpenal2(X)  // gradient fonction de penalisation interne
    d=X(1);h=X(2);
    if (c(X)<0) then
      z=[0;0]
    else
      z=[%pi*d*h/2/c(X);%pi*d^2/4/c(X)]      
    end    
endfunction
//
//////////////////////////////////////////////////////
// solution exacte
//
d0=(2400/%pi)^(1/3);h0=1200/(%pi*d0^2);
///////////////////////////////////////////////////////
Ntype=evstr(x_dialog('type de pénalisation (1: externe, 2: interne)','2'));
mu=evstr(x_dialog('taux de pénalisation (mu)','10'));
X0=evstr(x_dialog('point initial','[10;12]'));
//
deff('z=ftot(X)','z=f(X)+mu*penal(X)');
deff('z=ftot2(X)','z=f(X)-mu*penal2(X)');
deff('z=gradftot(X)','z=gradf(X)+mu*gradpenal(X)');
deff('z=gradftot2(X)','z=gradf(X)-mu*gradpenal2(X)');
/////////////////////////////////////////////////////////
//
///////////////// tracés ////////////////////
//
Nt=200;
x=linspace(2,14,Nt);
y=linspace(2,14,Nt);
z=zeros(Nt,Nt);zpenal=zeros(Nt,Nt);
for i=1:Nt
    for j=1:Nt
       z(i,j)=f([x(i),y(j)]) 
       zpenal(i,j)=f([x(i),y(j)])+mu*penal([x(i),y(j)]); 
       zpenal2(i,j)=f([x(i),y(j)])-mu*penal2([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(d0,h0,-1)
 xtitle('pas de penalisation')
//
/////////////////////////////////////////////////////////
//
alphainit=1
beta=0.001
tau=0.3
//
Niter=40;
X=X0;
Xtot=X;

if (Ntype==1) then

 for i=1:Niter
  d=-gradftot(X);
  alpha=alphainit
  while ftot(X+alpha*d)>ftot(X)+beta*alpha*(d'*gradftot(X))
  alpha=tau*alpha;
  end    
 X=X+alpha*d;Xtot=[Xtot,X];      
 end  
 figure(2)
 clf()
 contour(x,y,zpenal,20)
 d=linspace(5,12,Nt);
 h=1200/%pi*ones(d)./(d.^2);
 plot2d(dplot,hplot)
 plot2d(d0,h0,-1)
 plot2d(Xtot(1,:),Xtot(2,:))
 xtitle('penalisation externe')
//    
else      
 for i=1:Niter
  d=-gradftot2(X); 
  alpha=alphainit
  while (ftot2(X+alpha*d)>ftot2(X)+beta*alpha*(d'*gradftot2(X)))|(c(X+alpha*d)<0)
  alpha=tau*alpha;
  end    
 X=X+alpha*d;Xtot=[Xtot,X];      
 end 
 figure(3)
 clf()
 contour(x,y,zpenal2,20)
 plot2d(dplot,hplot)
 plot2d(d0,h0,-1)
 plot2d(Xtot(1,:),Xtot(2,:))
 xtitle('penalisation interne')
//     
end
disp('valeur exacte:');disp([d0;h0])
disp('surface minimale canette:');disp(f([d0,h0]))

disp('valeur optimale obtenue:');disp(X)
disp('contrainte:');disp(c(X))

disp('surface canette:');disp(f(X))
disp('volume canette:');disp(300+c(X))
//



