function Xeval=evaluation(X);
 mu=1;
 [Npop,n2]=size(X);
 Xeval=X;
  for i=1:Npop;
   penal=min(%pi*X(i,1)^2*X(i,2)/4-300,0)^2;  
   Xeval(i,$)=%pi*X(i,1)^2/4+%pi*X(i,1)*X(i,2)+mu*penal;
  end
endfunction
///////////////////////////////
function bestX=best(X);
[Npop,n2]=size(X);
[y,k]=gsort(-X(:,$));
bestX=X(k(1),:);
endfunction
/////////////////////////////// 
function Xsel=selection(X) 
[Npop,n2]=size(X);
[y,k]=gsort(X(:,$))  
X=X(k,:);  
p=1:Npop;
roulette=cumsum(p)/sum(p);
Xsel=[];
for i=1:Npop 
  u=rand();isel=1;
  while (u>roulette(isel)) 
     isel=isel+1;
  end
  Xsel=[Xsel;X(isel,:)];
end  
endfunction
////////////////////////////////// 
function xcrois=croisement(x)
xcrois=x;
u=rand();
xcrois(1,:)=u*x(1,:)+(1-u)*x(2,:);
xcrois(2,:)=(1-u)*x(1,:)+u*x(2,:);
endfunction
//////////////////////////////////
function xmut=mutation(x,gen,Ngen,xmin,xmax)
bet=5;
if (rand()<1/2) then
   xmut=x+rand()*(xmax-x)*(1-(gen-1)/Ngen)^bet;
else
   xmut=x-rand()*(x-xmin)*(1-(gen-1)/Ngen)^bet;
end
endfunction
//////////MAIN //////////////////////////
Npop=300;Ngen=40;
n=2;
pc=0.2;pm=0.3;
xmin=1;xmax=20;
X=xmin+(xmax-xmin)*rand(Npop,n+1);
bestX=[];
for gen=1:Ngen
X=evaluation(X);   // evaluation
X=selection(X);    // selection
bestX=[bestX;best(X)];
Xnew=[];
for j=1:Npop/2
   u1=int(Npop*rand())+1;
   u2=int(Npop*rand())+1;
   if (rand()<pc) then 
     Xnew=[Xnew;croisement(X([u1,u2],:))]; //croisement  
   else
     Xnew=[Xnew;X([u1,u2],:)];
   end  
end
for j=1:Npop
  if (rand()<pm) then
    Xnew(j,:)=mutation(Xnew(j,:),gen,Ngen,xmin,xmax);  // mutation
  end 
end  
X=Xnew;
end
//////////TRACE //////////////////////////
function area = can(x)
d=x(1);h=x(2);
area=%pi*d^2/4+%pi*d*h
endfunction
///////////////////////////////
xset('window',0);clf
plot2d(Npop*(1:Ngen),bestX(:,$))
disp('meilleur individu (diameter,height):');disp(bestX($,1:n));
disp('volume:');disp(%pi*bestX($,1)^2*bestX($,2)/4); 
disp('surface:');disp(can(X($,1:n)));
Xoptim=[(2400/%pi)^(1/3) ,1200/(%pi*(2400/%pi)^(2/3))];
disp('solution exacte:');disp(Xoptim);
/////////////////////////////////////////
xset('window',1);clf;
xmin=1;xmax=12;ymin=1;ymax=12;
xcan=xmin:((xmax-xmin)/100):xmax;
ycan=ymin:((ymax-ymin)/100):ymax;
for i=1:101;
  for j=1:101;
    zcan(i,j)=can([xcan(i),ycan(j)]);
  end
  hcan(i)=300*4/(%pi*xcan(i)^2);
end
plot2d(bestX(:,1),bestX(:,2),rect=[xmin,ymin,xmax,ymax]);
plot2d(bestX($,1),bestX($,2),-2);
plot2d(Xoptim(1),Xoptim(2),-3);
plot2d(xcan,hcan,2,rect=[xmin,ymin,xmax,ymax]);
contour2d(xcan,ycan,zcan,10,rect=[xmin,ymin,xmax,ymax]);