//
//
//------------------------------------------------------
// robust ES for I beam 
disp('robust ES for I beam');
//
//-----------------------------------------------------
//
//
function y=beam(x)  // the function to optimize
n=prod(size(x));
y=2*x(2)*x(4)+x(3)*(x(1)-2*x(4));
endfunction
//
function y=penaltystress(x)
n=prod(size(x));
y1=180000*x(1)/(x(3)*(x(1)-2*x(4))^3+2*x(2)*x(4)*(4*x(4)^2+3*x(1)*(x(1)-2*x(4))));
y2=y1+15000*x(2)/((x(1)-2*x(4))*x(3)^3+2*x(4)*x(2)^3)-16
y3=5000/(1/12*x(3)*(x(1)-2*x(4))^3+2/12*x(2)*x(4)^3+2*x(2)*x(4)*((x(1)-x(4))/2)^2)-0.1;
y=max(y2,0)^2+max(y3,0)^2;
endfunction
//
function y=penaltyboud(x)  // the function to optimize
n=prod(size(x));
penalty=[max(10-x(1),0)^2,max(10-x(2),0)^2,max(9-10*x(3),0)^2,max(9-10*x(4),0)^2,max(2*x(4)-x(1),0)^2,min(80-x(1),0)^2,min(50-x(2),0)^2,min(5-x(3),0)^2,min(5-x(4),0)^2];
y=sum(penalty);
endfunction
//-----------------------------------------------------
//
//-----------------------------------------------------
function evalu=robust(x);  // evaluation of the whole population
global K std wu ws Nsamp
ech=[];
for ks=1:Nsamp
u=std*(2*rand(1,n,'u')-1);
//u=std*rand(1,n,'n'); 
   ech=[ech,beam(x+u)+K*(penaltyboud(x+u)+penaltystress(x+u))];
end
val1=mean(ech);
val2=st_deviation(ech);
evalu=wu*val1+ws*val2^2;
endfunction
//
//-----------------------------------------------------
function Xpar=evaluation(X);  // evaluation of the whole population
global K std wu ws Nsamp
Xpar=X;
[Npo,n]=size(X);n=n-1;
evalpop=zeros(Npo,1);
for i=1:Npo;
  evalpop(i,:)=robust(X(i,1:n));
end
Xpar(:,(n+1))=evalpop;
endfunction
//
//-----------------------------------------------------
function Xparnew=selection(Xpar,Xoff);
[mu,n]=size(Xpar);
[lambda,n]=size(Xoff);
n=n-1;
Xtot=[Xpar;Xoff];
[Fsort,ind]=sort(-Xtot(:,n+1));
Xtot=Xtot(ind,:);
Xparnew=Xtot(1:mu,:);
endfunction
//
//-----------------------------------------------------//
function xoff=mutate(xpar,sigma) // normal mutation
[np,n]=size(xpar);n=n-1;
xoff=xpar;
u=rand(1,n,'n');
xoff(1,1:n)=xpar(1,1:n)+sigma*u;
endfunction
//
////-------------------------------------------dynamic results displays
//
function result(Ntot,fmin,x,co)
nl=length(Ntot);
if (nl>1) then
xset('window',0);
plot2d([Ntot($-1),Ntot($)],[fmin($-1),fmin($)],co);
xtitle('convergence history','Neval','fmin');
end
xset('window',2);
xbasc();
X=[0,x(2),x(2),x(2)/2+x(3)/2,x(2)/2+x(3)/2,x(2),x(2),0,0,x(2)/2-x(3)/2,x(2)/2-x(3)/2,0,0]
Y=[0,0,x(4),x(4),x(1)-x(4),x(1)-x(4),x(1),x(1),x(1)-x(4),x(1)-x(4),x(4),x(4),0];
plot2d(X,Y,rect=[0,0,50,80],frameflag=3);
endfunction
//
function plotIbeam(x,wu,Nsamp);
//
X=[0,x(2),x(2),x(2)/2+x(3)/2,x(2)/2+x(3)/2,x(2),x(2),0,0,x(2)/2-x(3)/2,x(2)/2-x(3)/2,0,0]
Y=[0,0,x(4),x(4),x(1)-x(4),x(1)-x(4),x(1),x(1),x(1)-x(4),x(1)-x(4),x(4),x(4),0];
//
xset('window',3);
if (wu==1)&(Nsamp==1)
plot2d(X,Y,1,rect=[0,0,50,80],frameflag=3); 
//legends('non robust, prioity on mean',1,1);
end
if (wu==1)&(Nsamp>1)
plot2d(X,Y,2,rect=[0,0,50,80],frameflag=3);
//legends('robust prioity on mean',2,2);
end
if (wu==0.5)&(Nsamp>1)
plot2d(X,Y,3,rect=[0,0,50,80],frameflag=3); 
//legends('robust, equal priority mean/variance',3,3);
end
if (wu==0)&(Nsamp>1)
plot2d(X,Y,4,rect=[0,0,50,80],frameflag=3); 
//legends('robust, prioity on variance',4,4);
end
//
endfunction
//
//
//-----------------------------------------------------
//-----------------------------------------------------
//---------main program -------------------------------
//
xset('window',0);xbasc();
global K std wu ws Nsamp
//
// ----------------------------------------------------
//
K=1000;
n=4;                          // number of variables of the optimization problem
mu=30;                           //parent number       
lambda=300;                       // offspring number 
siginit=1;                       // initial mutation strength
alpha=1.1;                       // 1/5 rule coefficient
Ngen=40;
//
Nsamp=evstr(x_dialog('sampling number','300'));
std=evstr(x_dialog('variation range','0.5'));
wu=evstr(x_dialog('coefficient on the mean','1'));
ws=1-wu;
// ----------------------------------------------------
//
xmin=[10,10,0.9,0.9];
xmax=[80,50,5,5];
Xmin=ones(mu,1)*xmin;
Xmax=ones(mu,1)*xmax;
u=rand(mu,n);
pop=Xmin+(Xmax-Xmin).*u;  // random initialisation of the population
Xpar=[pop,rand(mu,1)];  Xoldpar=Xpar;
Xoff=Xpar;Ntot=[];
mineval=0;newval=zeros(mu,1);itfct=0;fmin=[];fmin2=[];
//
//-----------------------------------------------------
//-----------------------------------------------------
//
sigma=siginit;
Xpar=evaluation(Xpar);
Nevalf=mu;
//
//     MAIN LOOP : BEGINNING
//
//
//
for gen=1:(Ngen-1)
//
  for i=1:lambda
    u=int(mu*rand())+1;
    xpar=Xpar(u,:);
    xoff=mutate(xpar,sigma);             // GA mutation
    Xoff(i,:)=xoff;
  end
//
oldval=Xoff(:,n+1);
Xoff=evaluation(Xoff);  // evaluation 
Nevalf=Nevalf+Nsamp*lambda;
//
newval=Xoff(:,n+1);
ind=find(oldval>newval);
taux=length(ind)/lambda;
if (taux<0.2) then sigma=sigma/alpha; else sigma=alpha*sigma;end
//
Xpar=selection(Xpar,Xoff);  // selection
fmin=[fmin,Xpar(1,n+1)];
bestpar=Xpar(1,1:n);
Ntot=[Ntot,Nevalf];
result(Ntot,fmin,bestpar,1);       // dynamic result display
//
//     MAIN LOOP : ENDING
//
end
//-----------------------------------------------------
plotIbeam(bestpar,wu,Nsamp);
//-----------------------------------------------------
//
//
//
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(fmin($));
disp('evaluation number:');disp(Nevalf);
//
