智能优化算法程序代码集锦

人工蚂蚁算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%

function [x,y, minvalue] = AA(func)

% Example [x, y,minvalue] = AA('Foxhole')

clc;

tic;

subplot(2,2,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot 1 draw(func);

title([func, ' Function']);

%初始化各参数

Ant=100;%蚂蚁规模

ECHO=200;%迭代次数

step=0.01*rand(1);%局部搜索时的步长

temp=[0,0];

%各子区间长度

start1=-100;

end1=100;

start2=-100;

end2=100;

Len1=(end1-start1)/Ant;

Len2=(end2-start2)/Ant;

%P = 0.2;

%初始化蚂蚁位置

for i=1:Ant

X(i,1)=(start1+(end1-start1)*rand(1));

X(i,2)=(start2+(end2-start2)*rand(1));

%func=AA_Foxhole_Func(X(i,1),X(i,2));

val=feval(func,[X(i,1),X(i,2)]);

T0(i)=exp(-val);%初始信息素,随函数值大,信息素浓度小,反之亦

然 %%%%%***************************************** ****************************

end;

%至此初始化完成

for Echo=1:ECHO %开始寻优

%P0函数定义,P0为全局转移选择因子

a1=0.9;

b1=(1/ECHO)*2*log(1/2);

f1=a1*exp(b1*Echo);

a2=0.225;

b2=(1/ECHO)*2*log(2);

f2=a2*exp(b2*Echo);

if Echo<=(ECHO/2)

P0=f1;

else

P0=f2;

end;

%P函数定义，P为信息素蒸发系数

a3=0.1; b3=(1/ECHO).*log(9);

P=a3*exp(b3*Echo);

lamda=0.10+(0.14-0.1)*rand(1);%全局转移步长参数Wmax=1.0+(1.4-1.0)*rand(1);%步长更新参数上限

Wmin=0.2+(0.8-0.2)*rand(1);%步长更新参数下限

%寻找初始最优值

T_Best=T0(1);

for j=1:Ant

if T0(j)>=T_Best

T_Best=T0(j);

BestIndex=j;

end;

W=Wmax-(Wmax-Wmin)*(Echo/ECHO); %局部搜索步长更新参数

for j_g=1:Ant %全局转移概率求取,当该蚂蚁随在位置不是bestindex时

if j_g~=BestIndex

r=T0(BestIndex)-T0(j_g);

Prob(j_g)=exp(r)/exp(T0(BestIndex));

else%当j_g=BestIndex的时候进行局部搜索

if rand(1)<0.5

temp(1,1)=X(BestIndex,1)+W*step; temp(1,2)=X(BestIndex,2)+W*step;

else

temp(1,1)=X(BestIndex,1)-W*step; temp(1,2)=X(BestIndex,2)-W*step;

end;

Prob(j_g)=0;%bestindex的蚂蚁不进行全局转移

end;

X1_T=temp(1,1);

X2_T=temp(1,2);

X1_B=X(BestIndex,1);

X2_B=X(BestIndex,2);

%func1 =

AA_Foxhole_Func(X1_T,X2_T); %%%%%%%%%%%********* ******************************************

%F1_T=func1;

F1_T=feval(func,[X(i,1),X(i,2)]);

F1_B=feval(func,[X1_B,X2_B]);

%F1_T=(X1_T-1).^2+(X2_T-2.2).^2+1;

%func2 =

AA_Foxhole_Func(X1_B,X2_B); %%%%%%%%%%%%%******** *******************************************

%F1_B=func2;

%F1_B=(X1_B-1).^2+(X2_B-2.2).^2+1;

if exp(-F1_T)>exp(-F1_B)

X(BestIndex,1)=temp(1,1);

X(BestIndex,2)=temp(1,2);

end;

for j_g_tr=1:Ant

if Prob(j_g_tr)

X(j_g_tr,1)=X(j_g_tr,1)+lamda*(X(BestIndex,1)-

X(j_g_tr,1));%Xi=Xi+lamda*(Xbest-Xi)

X(j_g_tr,2)=X(j_g_tr,2)+lamda*(X(BestIndex,2)-

X(j_g_tr,2));%Xi=Xi+lamda*(Xbest-Xi)

X(j_g_tr,1)=bound(X(j_g_tr,1),start1,end1);

X(j_g_tr,2)=bound(X(j_g_tr,2),start2,end2);

else

X(j_g_tr,1)=X(j_g_tr,1)+((-

1)+2*rand(1))*Len1;%Xi=Xi+rand(-1,1)*Len1

X(j_g_tr,2)=X(j_g_tr,2)+((-

1)+2*rand(1))*Len2;%Xi=Xi+rand(-1,1)*Len2

X(j_g_tr,1)=bound(X(j_g_tr,1),start1,end1);

X(j_g_tr,2)=bound(X(j_g_tr,2),start2,end2);

end;

%信息素更新

subplot(2,2,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Plot 1

bar([X(BestIndex,1) X(BestIndex,2)],0.25);

%colormap (cool);

axis([0 3 -40 40 ]) ;

title ({date;['Iteration ',

num2str(Echo)]});

xlabel(['Min_x = ',num2str(X(BestIndex,1)),' ', 'Min_y = ', num2str(X(BestIndex,2))]);

for t_t=1:Ant

%func=AA_Foxhole_Func(X(t_t,1),X(t_t,2)); val1=feval(func,[X(t_t,1),X(t_t,2)]);

T0(t_t)=(1-P)*T0(t_t)+(exp(-

val1)); %**************************************** *********************************

end;

[c_iter,i_iter]=max(T0); %求取每代全局最优解 minpoint_iter=[X(i_iter,1),X(i_iter,2)];

%func3 =

AA_Foxhole_Func(X(i_iter,1),X(i_iter,2)); %%%%%%% %%*********************************************** ****************************

val2=feval(func,[X(i_iter,1),X(i_iter,2)]); minvalue_iter= val2;

%minvalue_iter=(X(i_iter,1)-

1).^2+(X(i_iter,2)-2.2).^2+1;

min_local(Echo)=minvalue_iter;%保存每代局部最优解

subplot(2,2,3);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Plot 2plot(X(BestIndex,1),X(BestIndex,2),'rs','MarkerFa ceColor','r', 'MarkerSize',8), grid on;

title (['Global Min Value = ',

num2str(minvalue_iter)]);

hold on;

plot(X(:,1),X(:,2),'g.'), pause (0.02); hold off ;

axis([-100 100 -100 100]);

grid on;

%将每代全局最优解存到min_global矩阵中

if Echo >= 2

if min_local(Echo)

min_global(Echo)=min_local(Echo);

else

min_global(Echo)=min_global(Echo-1);

end;

else

min_global(Echo)=minvalue_iter;

end;

subplot(2,2,4);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 3

min_global=min_global';

index(:,1)=1:ECHO;

plot(Echo, min_global(Echo),'y*')

%axis([0 ECHO 0 10]);

hold on;

title ([func,' (X) = ',

num2str(minvalue_iter)],'Color','r');

xlabel('iteration');

ylabel('f(x)');

grid on;

end;%ECHO循环结束

[c_max,i_max]=max(T0);

minpoint=[X(i_max,1),X(i_max,2)];

%func3 =

AA_Foxhole_Func(X(i_max,1),X(i_max,2)); %%%****** ************************************************* ******************

%minvalue = func3;

minvalue=feval(func,[X(i_max,1),X(i_max,2)]);

x=X(BestIndex,1);

y=X(BestIndex,2);

runtime=toc

人工免疫算法

function [x,y,fx,vfx,vmfit,P,vpm] =

AI(func,gen,n,pm,per);

% Example [x,y,fx] = AI('Foxhole')

subplot(2,2,1);

draw(func);

title( [func, ' Function']);

if nargin == 1,

% gen = 200; n = round(size(P,1)/2); pm =

0.0005; per = 0.0; fat = 10;

%gen = 250; n = size(P,1); pm = 0.01; per = 0.0; fat = .1;

P = cadeia(200,44,0,0,0);

gen = 40; n = size(P,1); pm = 0.2; per = 0.0; fat = 0.1;

end;

while n <= 0,

n = input('n has to be at least one. Type a new value for n: ');

end;

xmin=-100;

xmax=100;

ymin=-100;

ymax=100;

x = decode(P(:,1:22),xmin,xmax); y =

decode(P(:,23:end),ymin,ymax);

%fit = eval(f);

%fit=AI_Foxhole_Func(x,y);%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%

fit=feval(func,[x' y']);

%imprime(1,vxp,vyp,vzp,x,y,fit,1,1);

% Hypermutation controlling parameters

pma = pm; itpm = gen; pmr = 0.8;

% General defintions

vpm = []; vfx = []; vmfit = []; valfx = 1; [N,L] = size(P); it = 0; PRINT = 1;

% Generations

while it <= gen & valfx <= 100,

x = decode(P(:,1:22),xmin,xmax); y =

decode(P(:,23:end),ymin,ymax); T = []; cs = [];

%fit = eval(f);

%fit=AI_Foxhole_Func(x,y);%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%

fit=feval(func,[x' y']);

[a,ind] = sort(fit);

valx = x(ind(end-n+1:end)); valy = y(ind(end-n+1:end));

fx = a(end-n+1:end); % n best individuals (maximization)

% Reproduction [T,pcs] = reprod(n,fat,N,ind,P,T);

% Hypermutation

M = rand(size(T,1),L) <= pm;

T = T - 2 .* (T.*M) + M;

T(pcs,:) = P(fliplr(ind(end-n+1:end)),:);

% New Re-Selection (Multi-peak solution)

x = decode(T(:,1:22),xmin,xmax); y =

decode(T(:,23:end),ymin,ymax);

%fit=AI_Foxhole_Func(x,y);%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%

fit=feval(func,[x' y']);

%fit = eval(f);

pcs = [0 pcs];

for i=1:n,

[out(i),bcs(i)] =

min(fit(pcs(i)+1:pcs(i+1))); % Mimimazion problem %%%*************************

bcs(i) = bcs(i) + pcs(i);

end;

P(fliplr(ind(end-n+1:end)),:) = T(bcs,:);

% Editing (Repertoire shift)

nedit = round(per*N); it = it + 1;

P(ind(1:nedit),:) = cadeia(nedit,L,0,0,0);

pm = pmcont(pm,pma,pmr,it,itpm); valfx =

min(fx); %*************************************** **********************

vpm = [vpm pm]; vfx = [vfx valfx]; vmfit = [vmfit mean(fit)];

disp(sprintf('It.: %d pm: %.4f x: %2.2f

y: %2.2f Av.: %2.2f

f(x,y): %2.3f',it,pm,valx(1),valy(1),vmfit(1),val fx));

subplot(2,2,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 1

bar([valx(1) valy(1)],0.25);

axis([0 3 -40 40 ]) ;

title (['Iteration ', num2str(it)]); pause (0.1);

subplot(2,2,3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 2

plot(valx(1),valy(1),'rs','MarkerFaceColor','r', 'MarkerSize',8)

hold on;

%plot(x(:,1),x(:,2),'k.');

set(gca,'Color','g')

hold off;

grid on;

axis([-100 100 -100 100 ]) ;

title(['Global Min =

',num2str(valfx)]);

xlabel(['Min_x= ',num2str(valx(1)),'

Min_y= ',num2str(valy(1))]);

subplot(2,2,4); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 3 plot(it,valfx ,'k.')

axis([0 gen 0 10]);

hold on;

title ([func ,'(X) = ', num2str(valfx)]);

xlabel('iteration');

ylabel('f(x)');

grid on;

end; % end while

%imprime(PRINT,vxp,vyp,vzp,x,y,fit,it,1);

x = valx(1); y = valy(1); fx =

min(fx); %***************************************

********************************

% x = P(ind(end),1:22); y = P(ind(end),23:44); fx

= max(fx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% Plot 4

% --------------------- %

% INTERNAL SUBFUNCTIONS

% --------------------- %

% Print

function [] =

imprime(PRINT,vx,vy,vz,x,y,fx,it,mit);

% x,fx -> actual values

% vxplot, vplot -> original (base) function

if PRINT == 1,

if rem(it,mit) == 0,

mesh(vx,vy,vz); hold on; axis([-100 100 -

100 100 0 500]);

xlabel('x'); ylabel('y'); zlabel('f(x,y)');

plot3(x,y,fx,'k*'); drawnow; hold off;

end;

% Reproduction

function [T,pcs] = reprod(n,fat,N,ind,P,T);

% n -> number of clones

% fat -> multiplying factor

% ind -> best individuals

% T -> temporary population

% pcs -> final position of each clone

if n == 1,

cs = N;

T = ones(N,1) * P(ind(1),:);

else,

for i=1:n,

% cs(i) = round(fat*N/i);

cs(i) = round(fat*N); pcs(i) = sum(cs);

T = [T; ones(cs(i),1) * P(ind(end-i+1),:)];

end;

% Control of pm

function [pm] = pmcont(pm,pma,pmr,it,itpm);

% pma -> initial value

% pmr -> control rate

% itpm -> iterations for restoring

if rem(it,itpm) == 0,

pm = pm * pmr;

if rem(it,10*itpm) == 0,

pm = pma;

end;

% Decodify bitstrings

function x = decode(v,min,max);

% x -> real value (precision: 6)

% v -> binary string (length: 22)

v = fliplr(v); s = size(v);

aux = 0:1:21; aux = ones(s(1),1)*aux;

x1 = sum((v.*2.^aux)');

x = min + (max-min)*x1 ./ 4194303;

function [ab,ag] = cadeia(n1,s1,n2,s2,bip)

%default parameter value seeting

if nargin == 2,

n2 = n1; s2 = s1; bip = 1;

elseif nargin == 4,

bip = 1;

end;

% Antibody (Ab) chains

ab = 2 .* rand(n1,s1) - 1;%create n1 row s1 column array, its value range is between -1 or 1

if bip == 1,

ab = hardlims(ab);

else,

ab = hardlim(ab);

end;

% Antigen (Ag) chains

ag = 2 .* rand(n2,s2) - 1;

if bip == 1,

ag = hardlims(ag);

else,

ag = hardlim(ag);

end;

% End Function CADEIA

%------免疫粒子群优化算法（Artificial Immune - Particle Swarm Optimization）

function [x,y,Result]=PSO_AI(func)

% Example [x, y,minvalue] = PSO_AI('Foxhole') clc;

subplot(2,2,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot 1 draw(func);

title([func, ' Function']);

tic

format long;

%------给定初始化条件----------------------------------------------

c1=1.4962; %学习因子1

c2=1.4962; %学习因子2

w=0.7298; %惯性权重

MaxDT=200; %最大迭代次数

D=2; %搜索空间维数（未知数个数）N=100; %初始化群体个体数目

eps=10^(-20); %设置精度(在已知最小值时候用)

DS=10; %每隔DS次循环就检查最优个体是否变优

replaceP=0.6; %粒子的概率大于replaceP将被免疫替换

minD=1e-015; %粒子间的最小距离

Psum=0; %个体最佳的和

range=100;

count = 0;

%------初始化种群的个体------------

for i=1:N

for j=1:D

x(i,j)=-range+2*range*rand; %随机初始化位置

v(i,j)=randn; %随机初始化速度

end

%------先计算各个粒子的适应度，并初始化Pi和Pg----------------------

for i=1:N

%p(i)=Foxhole(x(i,:),D); %fitness是计算各个粒子适应度的函数，见文件

fitness.m %%%%%%%%****************************** ***********************

p(i)=feval(func,x(i,:));

y(i,:)=x(i,:);

end

pg=x(1,:); %Pg为全局最优

for i=2:N

feval(func,x(i,:))

end

%------进入主要循环，按照公式依次迭代，直到满足精度要求------------

for t=1:MaxDT

for i=1:N

v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-

x(i,:))+c2*rand*(pg-x(i,:));

x(i,:)=x(i,:)+v(i,:);

feval(func,x(i,:))

p(i)=feval(func,x(i,:)); %%%%%%%%%%%%%%%******** ************************************************* *****************************

y(i,:)=x(i,:);

end

p(i)

pg=y(i,:);

subplot(2,2,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 1

bar(pg,0.25);

axis([0 3 -40 40 ]) ;

title (['Iteration ', num2str(t)]); pause (0.1);

subplot(2,2,3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 2

plot(pg(1,1),pg(1,2),'rs','MarkerFaceColor','r', 'MarkerSize',8)

hold on;

plot(x(:,1),x(:,2),'k.');

set(gca,'Color','g')

hold off;

grid on;

axis([-100 100 -100 100 ]) ;

title(['Global Min =

',num2str(p(i))]);

xlabel(['Min_x= ',num2str(pg(1,1)),' Min_y= ',num2str(pg(1,2))]);

end

Pbest(t)=feval(func,pg) ; %%%%%%%%************* ************************************************* **********************************************

% if Foxhole(pg,D)

% break;

% end

%-----------开始进行免疫－－－－－－－－－－－－－－－－

if t>DS

if mod(t,DS)==0 && (Pbest(t-DS+1)-

Pbest(t))<1e-020 %如果连续DS代数，群体中的最优没有明显变优，则进行免疫.

%在函数测试的过程中发现，经过一定代数的更新，个体最优不完全相等，但变化非常非常小，

%我认为这个时候也应用免疫了，所以我没有用“Pbest(t-DS+1)＝Pbest(t)”作为判断条件，

%不过“(Pbest(t-DS+1)-Pbest(t))<1e-020”是否合理也值得探讨。

for i=1:N %先计算出个体最优的和

Psum=Psum+p(i);

end

for i=1:N %免疫程序

for j=1:N %计算每个个体与个体i的距离

distance(j)=abs(p(j)-p(i));

end

num=0;

for j=1:N %计算与第i个个体距离小于minD的个数

if distance(j)

num=num+1;

end

PF(i)=p(N-i+1)/Psum; %计算适应度概率

PD(i)=num/N; %计算个体浓度

a=rand; %随机生成计算替换概率的因子

PR(i)=a*PF(i)+(1-a)*PD(i); %计算替换概率

end

for i=1:N

if PR(i)>replaceP

x(i,:)=-

range+2*range*rand(1,D);

subplot(2,2,4);

axi; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% Plot 4

plot(x(i,1),x(i,2),'k*'); grid on;

axis([-100 100 -100 100 ]) ;

title(['New Min =

',num2str( feval(func,x(i,:)))]); %%%%%%%%%***** ************************************************* *************************

xlabel(['Immune

',num2str(count)]); pause (0.2);

count=count+1;

end

%------最后给出计算结果

x=pg(1,1);

y=pg(1,2);

Result=feval(func,pg); %%%%%%%%%%%************** ************************************************* ************************************************* *

toc

%------算法结束-------------------------

function probabolity(N,i)

PF=p(N-i)/Psum;%适应度概率

disp(PF);

for jj=1:N

distance(jj)=abs(P(jj)-P(i));

end

num=0;

for ii=1:N

if distance(ii)

num=num+1;

end

PD=num/N; %个体浓度

PR=a*PF+(1-a)*PD; %替换概率

% result=PR;

差分进化算法

function [sol]= DE(func)

%% Example [sol]= DE('Foxhole')

tic

popsize= 100;

lb=[-100 -100];

ub=[100 100];

[sol ]= diffevolve(func, popsize, lb, ub);

toc

end

function [sol, fval, evals] = diffevolve(varargin) %DIFFEVOLVE Differential Evolution Optimization

% Usage:

% sol = DIFFEVOLVE(PROBLEM)

% sol = DIFFEVOLVE(func, popsize, lb, ub)

% sol = DIFFEVOLVE(func, popsize, lb, ub,

'option1', value1, ...)

% [sol, fval] = DIFFEVOLVE(...)

% [sol, fval, evals] = DIFFEVOLVE(...)

% DIFFEVOLVE(func, popsize, lb, ub) tries to

find the global optimum of

% the fitness-function [func], using a transversal differential evolution

% strategy. The population size set by [popsize], and the boundaries for

% each dimension are set by the vectors [lb]

and [ub], respectively.

% [sol, fval, evals] = DIFFEVOLVE(...) returns the trial vector found to

% yield the global minimum in [sol], and the corresponding function value

% by [fval]. The total amount of function evaluations that the algorithm

% performed is returned in [evals].

% The function [func] must accept a vector argument of length [N], equal to

% the number of elements in the vectors [lb] or [ub]. Also, the function

% must be vectorized so that inserting a matrix of [popsize]x[N] will return

% a vector of size [popsize]x[1] containing the corresponding function values

% for the [N] trial vectors.

% The default control parameters DIFFEVOLVE uses, are

% -1 < F < 1 scaling parameter for the difference vector.

% By default, it is randomly selected from the % range [-1, 1] for each individual in each

% iteration.

% crossconst = 0.9 Probability for crossover.

% convergence = 100 Maximum amount of iterations the best

% individual ever found should remain the best

% individual before the algorithm converges.

% These values, as well as additional options, may be given manually in

% any extra input variables. Additionally, DIFFEVOLVE may be called as

% DIFFEVOLVE(PROBLEM), where [PROBLEM] is a complete problem-structure

% constructed by HEURSET.

% Some optimizations require repeated evaluation of a function that takes

% in the order of hours per evaluation. In those cases, it might be

% convenient to interrupt the optimization

after a few days and use the

% results thus far obtained. Unfortunately, usually after you interrupt a

% function, its results are lost. For that reason, DIFFEVOLVE creates two

% global variables, [DIFFEVOLVE_bestind] and [DIFFEVOLVE_bestfval], which

% store the best individual and function value thus far found. When the

% optimization is interrupted, the intermediate results from the

% optmization are still available. Once an optimization has succesfully

% converged, these variables are deleted from the workspace again.

% See also heurset, genetic, swarm, simanneal, godlike.

% Author: Rody P.S. Oldenhuis

% Delft University of Technology

% E-mail: oldenhuis@dds.nl

% Last edited: 28/Feb/2008

%% initialize & parse input

% initial (default) values

skippop = false;

options = heurset;

[pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...

convvalue, crossconst, Flb, Fub, n] = parseprob(options);

% temporary globals

global DIFFEVOLVE_bestind DIFFEVOLVE_bestfval

% common inputs

if (nargin >= 4)

func = varargin{1};

popsize = varargin{2};

lb = varargin{3};

ub = varargin{4};

end

% with additional options

if (nargin >= 5)

if ~isstruct(varargin{5})

options = heurset(varargin{5:end});

else

options = varargin{5};

end

[dum1, dum2, dum3, dum4, dum5, grace, display, maxfevals, convmethod,...

convvalue, crossconst, Flb, Fub, n] = parseprob(options);

end

% if called from GODLIKE

if (nargin == 2)

problem = varargin{2};

% errortrap

if ~isstruct(problem)

error('PROBLEM should be a structure. Type `help heurset` for details.')

end

[pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...

convvalue, crossconst, Flb, Fub, n] = parseprob(problem);

% make sure the options are correct convmethod = 'maxiters';

grace = 0;

display = false;

skippop = true;

end

% if given a full problem structure

if (nargin == 1)

problem = varargin{1};

% errortrap

if ~isstruct(problem)

error('PROBLEM should be a structure. Type `help heurset` for details.')

end

[pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...

convvalue, crossconst, Flb, Fub, n] = parseprob(problem);

end

% initialize convergence method

if strcmpi(convmethod, 'exhaustive')

convergence = 1; maxiterinpop = convvalue;

maxiters = inf;

elseif strcmpi(convmethod, 'maxiters')

convergence = 2;

maxiterinpop = inf;

maxiters = convvalue;

elseif strcmpi(convmethod, 'achieveFval')

convergence = 3;

%errortrap

if isempty(convvalue)

error('Please define function value

to achieve.')

end

maxiterinpop = inf;

maxiters = inf;

else

convergence = 1;

maxiterinpop = convvalue;

end

% problem dimensions

dims = size(lb, 2);

% errortraps

if ( (size(lb, 2) ~= 1 || size(ub, 2) ~= 1)

&&...

(size(lb, 2) ~= dims || size(ub, 2) ~=

dims) )

error('Upper- and lower boundaries must

be the same size.')

end

if (popsize < 3)

error('DIFFEVOLVE needs a population size

of at least 3.')

end

%% initial values

% iteration-zero values

oldbestfit = inf;

improvement = maxiterinpop;

iters = 0;

evals = 0;

sizepop = [popsize, dims];

fitnew = inf(popsize, 1);

firsttime = true;

converging1 = false;

converging2 = false;

% boundaries

if (numel(lb) == 1)

mins = repmat(lb, popsize, dims);

maxs = repmat(ub, popsize, dims);

end

if (numel(lb) == numel(ub) && size(lb, 2) ==

dims)

mins = repmat(lb, popsize, 1);

maxs = repmat(ub, popsize, 1);

end

range = (maxs - mins);

% initialize population

if (~skippop)

pop = repmat(ub, popsize, 1) - repmat(ub-lb, popsize, 1).* rand(popsize, dims);

end

newpop = pop;

prevpop = pop;

% Display strings

if display

fprintf(1, ['\nNote that when DIFFEVOLVE

is interrupted, the best solution and

function\n',...

'value are saved in global variables DIFFEVOLVE_bestind and

DIFFEVOLVE_bestfval.\n\n']);

fprintf(1, 'Optimizing with ');

switch convergence

case 1

fprintf(1, 'exhaustive convergence criterion...\n');

case 2

fprintf(1, 'maximum iterations convergence criterion...\n');

case 3

fprintf(1, 'goal-attain convergence criterion...\n');

end

fprintf(1, 'Current best solution has

cost of ------------- ');

overfevals1 = '\n\nCould not find better solution within given maximum';

overfevals2 = '\nof function evaluations. Increase MaxFevals option.\n\n';

fitbackspace =

'\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b';

converge1bck = [fitbackspace,

'\b\b\b\b\b\b\b\b\b\b\b'];

converge2bck = [fitbackspace,

'\b\b\b\b\b\b\b'];

convergestr = ', possibly converging: '; itersstr = ', iterations left: ';

end

%% Differential Evolution loop

% loop while one of the convergence criteria

is not met

while (improvement > 0 && iters <= maxiters)

% evaluate fitnesses and adjust

population

fitold = fitnew;

try

fitnew = feval(func, newpop);

catch

error('diffevolve:fevalerror', ...

['Diffevolve cannot continue because the supplied cost function ',...

'gave the following

error:\n %s'], lasterr);

end

if (numel(fitnew) ~= popsize)

error('diffevolve:function_not_vectorized_or_inco rrect_dimensions', ...

['The user-supplied cost

function does not appear to get enough arguments,\n',...

'or is not vectorized properly. Make sure the dimensions of the limits [Lb]\n',...

'and [Ub] are the same as the required input vector for the function, or

that\n',...

'the function is vectorized properly.'])

end

prevpop = newpop;

pop(fitnew < fitold, :) = newpop(fitnew < fitold, :);

% increase number of function evaluations evals = evals + numel(fitnew);

% improving solutions

[bestindfit, ind] = min(fitnew);

bestind = newpop(ind, :);

if (bestindfit < oldbestfit)

% new best individual

oldbestfit = bestindfit;

oldbestind = bestind;

improvement = maxiterinpop;

% assign also the globals

DIFFEVOLVE_bestind = bestind;

DIFFEVOLVE_bestfval = bestindfit;

% display improving solution

if display

if converging1

fprintf(1, converge1bck);

pause(0.05)

end

if converging2

fprintf(1, converge2bck);

pause(0.05)

end

if (oldbestfit < 0)

fprintf('\b')

end

fprintf(1, fitbackspace);

fprintf(1, '%1.8e', oldbestfit); converging1 = false;

firsttime = true;

pause(0.05)

end

% Check for convergence

if (evals > maxfevals)

% maximum allowed function

evaluations has been superceded

fprintf(overfevals1);

fprintf(overfevals2);

break

end

switch convergence

% exhaustive

case 1

if (oldbestfit <= bestindfit) && (oldbestfit ~= inf)

if (improvement <= 100) && display

if ~converging1

fprintf(1, convergestr);

fprintf(1, '%3.0f', improvement-1);

converging1 = true;

pause(0.05)

else

fprintf(1,

'\b\b\b%3.0f', improvement-1);

pause(0.05)

end

improvement = improvement - 1;

end

% max iterations

case 2

if display && (maxiters-iters < 100)

if firsttime

fprintf(1, itersstr);

fprintf(1, '%3.0f', maxiters-iters);

pause(0.05)

else

fprintf(1, '\b\b\b%3.0f', maxiters-iters);

pause(0.05)

end

firsttime = false;

converging2 = true;

end

iters = iters + 1;

% goal-attain

case 3

if (oldbestfit <= convvalue)

break;

end

% Transversal Differential Evolution

for i = 1:popsize

for j = 1:n base = 1; d1 = 1; d2 = 1;

while base == d1 || base == d2 || d1 == d2

base = round(rand*(popsize-1))+1;

d1 = round(rand*(popsize-1))+1;

d2 = round(rand*(popsize-1))+1;

end

replace = rand < crossconst |

rand == i;

if replace

F = (Fub-Flb)*rand + Flb;

newpop(i, :) = pop(base, :) + F*(pop(d1, :) - pop(d2, :));

else

newpop(i, :) = pop(i, :);

end

% enforce boundaries

outsiders1 = (newpop < mins);

outsiders2 = (newpop > maxs);

if any(outsiders1(:)) ||

any(outsiders2(:))

reinit = rand(sizepop).*range + mins; newpop(outsiders1) =

reinit(outsiders1);

newpop(outsiders2) =

reinit(outsiders2);

end

%% (pre-) end values

% if solution has been found

if isfinite(oldbestfit)

% when called normally

if (~skippop)

fval = oldbestfit;

sol = oldbestind;

% when called from GODLIKE

else

fval = fitnew;

sol = prevpop;

end

% display final message

if display

fprintf(1, '\nDifferential Evolution optimization algorithm has converged.\n');

pause(0.05)

end

% all trials might be infeasible

else

fprintf(1,'\n');

warning('diffevolve:no_solution',...

'DIFFEVOLVE was unable to find any solution that gave finite values.')

fval = oldbestfit;

sol = NaN;

end

%% Grace function evaluations

if (grace > 0)

% display progress

if display

fprintf(1, 'Performing direct-search...');

pause(0.05)

end

% perform direct search

options = optimset('TolX', eps,

'MaxFunEvals', grace, 'TolFun', ...

eps, 'MaxIter', 1e4, 'Display', 'off');

[soltry, fvaltry] = fminsearch(func, sol, options);

% enforce boundaries

if ~any(soltry >= ub | soltry <= lb)

sol = soltry;

fval = fvaltry;

end

evals = evals + grace;

end

%% finalize

% display progress

if display

fprintf(1, 'All done.\n\n');

pause(0.05)

end

% clear the temp globals

clear global DIFFEVOLVE_bestind

DIFFEVOLVE_bestfval

end

% parser function to easily parse the input arguments

function [pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...

convvalue, crossconst, Flb, Fub, n] = parseprob(problem)

func = problem.costfun;

popsize = problem.popsize;

lb = problem.lb; ub = problem.ub;

grace = problem.grace;

display = problem.display;

maxfevals = problem.maxfevals;

convmethod = problem.conv.method;

convvalue = problem.conv.value;

pop = problem.DE.pop;

crossconst = problem.DE.crossconst;

Flb = problem.DE.Flb;

Fub = problem.DE.Fub;

n = problem.DE.n;

end

function problem = heurset(varargin)

%HEURSET Set options for the various heuristic optimizers

% Usage:

% heurset (to display a quick summary of all the options)

% options = heurset('option1', value1,

'option2', value2, ...)

% PROBLEM = heurset('option1', value1,

'option2', value2, ...)

% HEURSET is an easy way to set all options for the global optimization

% algorithms SWARM, DIFFEVOLVE, GENETIC, SIMANNEAL and GODLIKE. All

% options, and their possible values are:

================================================= =====================

% To define a full problem (PROBLEM =

heurset(...)):

================================================= =====================

% CostFunction : function handle to the N-dimensional objective function

% Lb : [1 x N] vector defining the lower bounds for each

% dimension.

% Ub : [1 x N] vector defining the upper bounds for each

% dimension.

% Popsize : the population size to be used by the heuristic metod.

================================================= =====================

% General options for all optimizers:

================================================= =====================

% Grace : positive scalar. When nonzero, this will pass the result

% found after convergence to the FMINSEARCH function,

% allowing it this many function evaluations (on top of the

% "MaxFevals" option, see two options down). This function

% is part of the optimization toolbox, so it should be set

% to zero when this is not present. The default is 0.

% Display : string, either 'on' or 'off'. Determines whether progress

% will be displayed during the optimizations. The default

% is "on".

% MaxFevals : positive scalar, defining the maximum number of

% allowable function evaluations. The default is 1e6.

% Maxiters : positive scalar. This enables the "maximum iterations"

% convergence criterion. The optimizer will be "converged"

% when this number of iterations has been performed.

% BestInPop : positive scalar. This enables the "exhaustive" convergence

% criterion--The best individual ever found needs to remain

% the best for this many iterations before the algorithm is

% "converged".

% AchieveFval : scalar. This enables the "goal-attain" convergence

% criterion--basically, the algorithm continues until

% either this function value has been reached, or the

% maximum number of iterations has been superceeded.

% *NOTE*: only the last convergence criterion given will be used. The

% default convergence criterion is the "BestInPop", with a value of 100.

================================================= =====================

% Options specific to the Differential Evolution algorithm:

================================================= =====================

% Flb : scalar. This value defines the lower bound for the range

% from which the scaling parameter will be taken. The

% default value is -1.2.

% Fub : scalar. This value defines the upper bound for the range % from which the scaling parameter will be taken. The

% default value is +1.2. These two values may be set equal

% to each other, in which case the scaling parameter F is

% simply a constant.

% Crossconst : positive scalar. It defines the probability with which a

% new trial individual will be inserted in the new

% population. The default value is 0.9.

% n : positive scalar. This defines the number of transversal

% steps that the Differential Evolution takes every

% iteration. The default value is 10. Note that setting this

% value to 1 results in the Neoteric Differential Evolution

% Algorithm.

================================================= =====================

% Options specific to the Genetic Algorithm:

================================================= =====================

% Crossprob : positive scalar, defining the probability for crossover

% for each individual. The default value is 0.25.

% Mutprob : positive scalar, defining the mutation probability for

% each individual. The default value is 0.01.

================================================= =====================

% Options specific to the Simulated Annealing Algorithm:

================================================= =====================

% T0 : scalar. This is the initial temperature for all particles. The

% default value is 1.

% minT : scalar. This is the minimum temperature before convergence. The

% default value is 1e-8.

% k : scalar. This is the (normalized) Bolzmann constant. The default

% is 1.

% cool : function handle, with "maxiters", T0, T, and minT as parameters.

% This function defines the cooling schedule to be applied each

% iteration. The default is

@(T,T0,minT,maxiters)(minT/T0)^(1/maxiters).

================================================= =====================

% Options specific to the Particle Swarm Algorithm:

================================================= =====================

% eta1 : scalar < 4. This is the

'social factor', the

% acceleration factor in front of the difference with the

% particle's position and its neighorhood-best. The

% default value is 2. Note that negative values result in

% a Repulsive Particle Swarm algorithm.

% eta2 : scalar < 4. This is the

'cooperative factor', the

% acceleration factor in front of the difference with the

% particle's position and the location of the global

% minimum found so far. The default value is 2.

% eta3 : scalar < 4. This is the

'nostalgia factor', the

% acceleration factor in front of the difference with the

% particle's position and its personal-best. The default

% value is 0.5.

% omega : scalar. This is the 'inertial constant', the tendency of

% a particle to continue its motion undisturbed. The

% default value is 0.5.

% numneighbors : positive scalar. This defines the number of neighbors

% assigned to each particle. The default value is 5.

================================================= =====================

% Options specific to the GODLIKE Algorithm:

================================================= =====================

% AlgIters : positive scalar. This sets the number of iterations

% spent in each heuristic optimizer. The default value is

% 250.

% PutbackIters : positive scalar. This defines the number if iterations % before the best individual

ever found is inserted back

% into one of the populations.

The default value is 25.

% SwapIters : positive scalar. This sets the number of iterations

% before the different

populations are shuffled and

% swapped between the algorithms. The default value is 25.

% See also genetic, diffevolve, swarm, simanneal, godlike.

% Author: Rody P.S. Oldenhuis

% Delft University of Technology

% E-mail: oldenhuis@dds.nl

% Last edited: 28/Feb/2009

% display quick summary for HEURSET in case

of zero input- and output-arguments

if (nargin == 0) && (nargout == 0)

fprintf(1, 'CostFunction: [ function

handle | {[]} ]\n')

fprintf(1, ' Lb: [ (1 x N)

vector | {[]} ] \n')

fprintf(1, ' Ub: [ (1 x N)

vector | {[]} ] \n')

fprintf(1, ' Popsize: [ positive

scalar | {[]} ]\n\n')

fprintf(1, ' Grace: [ positive

scalar | {0} ]\n')

fprintf(1, ' Display: [ {on} | off ] \n')

fprintf(1, ' MaxFevals: [ positive

scalar | {[]} ]\n')

fprintf(1, ' Maxiters: [ positive

scalar | {[]} ]\n')

fprintf(1, ' BestInPop: [ positive

scalar | {100} ]\n')

fprintf(1, ' AchieveFval: [ scalar | {[]} ]\n\n')

fprintf(1, ' Flb: [ scalar | {-

1.2} ]\n')

fprintf(1, ' Fub: [ scalar |

{1.2} ]\n')

fprintf(1, ' Crossconst: [ positive

scalar | {0.9} ]\n')

fprintf(1, ' n: [ positive

scalar | {10} ]\n\n')

fprintf(1, ' Crossprob: [ positive

scalar | {0.5} ]\n')

fprintf(1, ' Mutprob: [ positive

scalar | {0.01} ]\n\n')

fprintf(1, ' T0: [ scalar |

{1} ]\n')

fprintf(1, ' minT: [ scalar | {1e-

8} ]\n')

fprintf(1, ' k: [ scalar |

{1} ]\n')

fprintf(1,[' cool: [ function handle | {@(T, T0, minT, maxiters)',...

'(minT/T0)^(1/maxiters)} ]\n\n'])

fprintf(1, ' eta1: [ scalar |

{2} ]\n')

fprintf(1, ' eta2: [ scalar |

{2} ]\n')

fprintf(1, ' eta3: [ scalar |

{0.5} ]\n')

fprintf(1, ' omega: [ scalar |

{0.5} ]\n')

fprintf(1, 'numneighbors: [ positive

scalar | {5} ]\n\n')

fprintf(1, ' AlgIters: [ positive

scalar | {10} ]\n')

fprintf(1, 'PutbackIters: [ positive

scalar | {250} ]\n')

fprintf(1, ' SwapIters: [ positive

scalar | {25} ]\n\n')

return

end

% create structure with default values if no input is given

if (nargin == 0) && (nargout ~= 0)

problem = struct;

problem.costfun = [];

problem.lb = [];

problem.ub = [];

problem.popsize = [];

problem.grace = 0;

problem.display = true;

problem.maxfevals = 1e6;

problem.conv.method = 'exhaustive';

problem.conv.value = 150;

% differential evolution

problem.DE.pop = [];

problem.DE.Flb = -1.2;

problem.DE.Fub = 1.2;

problem.DE.crossconst = 0.9;

problem.DE.n = 10;

% genetic algorithm

problem.GA.pop = [];

problem.GA.crossprob = 0.5;

problem.GA.mutprob = 0.01;

% simulated annealing

problem.SA.pop = [];

problem.SA.T0 = 1;

problem.SA.minT = 1e-8;

problem.SA.k = 1;

problem.SA.cool = @(T, T0, minT, maxiters) (minT/T0)^(1/maxiters);

% particle swarm

problem.PS.pop = [];

problem.PS.eta1 = 2;

problem.PS.eta2 = 2; problem.PS.eta3 = 0.5;

problem.PS.omega = 0.5;

problem.PS.numneigh = 5;

% GODLIKE

problem.GODLIKE.swapiters = 10;

problem.GODLIKE.algiters = 250;

problem.GODLIKE.putbackiters = 25;

% otherwise, create structure with fields according to user input

else

% default values

problem = heurset;

% errortrap

if (mod(nargin, 2) ~= 0)

error('Please provide values for all

the options.')

end

% loop through all the inputs, and use an "if-else cancer" to

% create the problem structure

for i = 1:2:nargin

option = varargin{i};

value = varargin{i+1};

% if option is not recognized,

continue to the next argument

if ~isa(option, 'char')

throwwarning(option, [], [], []);

continue;

end

% General options

if strcmpi(option, 'costfunction')

if ~isa(value, 'function_handle')

throwwarning('costfunction',

'function_handle', value);

continue;

end

problem.costfun = value;

elseif strcmpi(option, 'Lb')

if ~isnumeric(value)

throwwarning('Lb', 'double', value);

continue;

end

problem.lb = value;

elseif strcmpi(option, 'Ub')

if ~isnumeric(value)

throwwarning('Ub', 'double', value);

continue;

end

problem.ub = value;

elseif strcmpi(option, 'Popsize')

if ~isnumeric(value)

throwwarning('Popsize', 'double', value); continue ; end problem.popsize = value; elseif strcmpi(option, 'Grace') if ~isnumeric(value) throwwarning('Grace', 'double', value); continue ; end problem.grace = value; elseif strcmpi(option, 'Display') if ~ischar(value) throwwarning('Display', 'char', value); continue ; end if strcmpi(value, 'off') problem.display = false; else problem.display = true; end elseif strcmpi(option, 'MaxFevals') if ~isnumeric(value) throwwarning('MaxFevals', 'double', value); continue ; end problem.maxfevals = value; elseif strcmpi(option, 'MaxIters') if ~isnumeric(value) throwwarning('MaxIters', 'double', value); continue ; end problem.conv.method = 'MaxIters'; problem.conv.value = value; elseif strcmpi(option, 'BestInPop') if ~isnumeric(value) throwwarning('BestInPop', 'double', value); continue ; end problem.conv.method = 'exhaustive'; problem.conv.value = value; elseif strcmpi(option, 'AchieveFval') if ~isnumeric(value) throwwarning('AchieveFval', 'double', value); continue ; end problem.conv.method = 'AchieveFval'; problem.conv.value = value; % options specific to Differential Evolution elseif strcmpi(option, 'Flb') if ~isnumeric(value) throwwarning('Flb', 'double', value); continue ; end problem.DE.Flb = value; elseif strcmpi(option, 'Fub') if ~isnumeric(value) throwwarning('Fub', 'double', value); continue ; end problem.DE.Fub = value; elseif strcmpi(option, 'crossconst') if ~isnumeric(value) throwwarning('crossconst', 'double', value); continue ; end problem.DE.crossconst = value; elseif strcmpi(option, 'n') if ~isnumeric(value) throwwarning('n', 'double', value); continue ; end problem.DE.n = value; % options specific to Genetic Algorithm elseif strcmpi(option, 'MutationProb') if ~isnumeric(value) throwwarning('mutationprob', 'double', value); continue ; end options.GA.mutprob = value; elseif strcmpi(option, 'CrossProb') if ~isnumeric(value) throwwarning('CrossProb', 'double', value); continue ; end options.GA.crossprob = value; % options specific to Simulated Annealing elseif strcmpi(option, 'T0') if ~isnumeric(value) throwwarning('T0', 'double', value); continue ; end problem.SA.T0 = value; elseif strcmpi(option, 'minT') if ~isnumeric(value) throwwarning('minT', 'double', value); continue ; end problem.SA.minT = value; elseif strcmpi(option, 'k')

if ~isnumeric(value)

throwwarning('k', 'double', value);

continue;

end

problem.SA.k = value;

elseif strcmpi(option, 'cool')

if ~isa(value, 'function_handle')

throwwarning('cool',

'function_handle', value);

continue;

end

problem.SA.cool = value;

% options specific to Particle Swarm Optimization

elseif strcmpi(option, 'eta1')

if ~isnumeric(value)

throwwarning('eta1', 'double', value);

continue;

end

problem.PS.eta1 = value;

elseif strcmpi(option, 'eta2')

if ~isnumeric(value)

throwwarning('eta2', 'double', value);

continue;

end

problem.PS.eta2 = value;

elseif strcmpi(option, 'eta3')

if ~isnumeric(value)

throwwarning('eta3', 'double', value);

continue;

end

problem.PS.eta3 = value;

elseif strcmpi(option, 'omega')

if ~isnumeric(value)

throwwarning('omega',

'double', value);

continue;

end

problem.PS.omega = value;

elseif strcmpi(option, 'numneighbors')

if ~isnumeric(value)

throwwarning('numneighbors',

'double', value);

continue;

end

problem.PS.numneigh = value;

% options specific to GODLIKE

algorithm

elseif strcmpi(option, 'swapiters')

if ~isnumeric(value)

throwwarning('swapiters',

'double', value);

continue;

end

problem.GODLIKE.swapiters = value;

elseif strcmpi(option, 'algiters')

if ~isnumeric(value)

throwwarning('algiters',

'double', value);

continue;

end

problem.GODLIKE.algiters = value;

elseif strcmpi(option, 'putbackiters')

if ~isnumeric(value)

throwwarning('putbackiters',

'double', value);

continue;

end

problem.GODLIKE.putbackiters = value;

else

error('heurset:incorrectoption', ...

'Unrecognized option: %s', option);

end

function throwwarning(option, required, given, varargin)%#ok

if nargin == 3

provided = whos('given');

provided = provided.class;

warning('heurset:incorrectvalue', ...

['Incorrect type given for option "%s";\n',...

'required type is %s,

got %s.\n',...

'Using default...'], option, required, provided);

else

warning('heurset:incorrectoption', ...

'Ignoring option `%s`',

num2str(option))

end

上帝算法

function [sol]=God(func)

%%% Example : [sol]=God('Foxhole')

tic

popsize= 100;

lb=[-100 -100];

ub=[100 100];

[sol ]= godlike(func, popsize, lb, ub);

toc

end

function [sol, fval, evals] = godlike(varargin) %GODLIKE Optimization using the GODLIKE algorithm

% Usage:

% sol = GODLIKE(PROBLEM)

% sol = GODLIKE(func, popsize, lb, ub)

% sol = GODLIKE(func, popsize, lb, ub,

'option1', value1, ...)

% [sol, fval] = GODLIKE(...)

% [sol, fval, evals] = GODLIKE(...)

% GODLIKE(func, popsize, lb, ub) tries to find the global optimum of

% the fitness-function [func], using the GODLIKE-algorithm. The

% population size set by [popsize], and the boundaries for each

% dimension are set by the vectors [lb] and [ub], respectively.

% [sol, fval, evals] = GODLIKE(...) returns the trial vector found to

% yield the global minimum in [sol], and the corresponding function value

% by [fval]. The total amount of function evaluations that the algorithm

% performed is returned in [evals].

% The function [func] must accept a vector argument of length [N], equal to

% the number of elements in the vectors [lb] or [ub]. Also, the function

% must be vectorized so that inserting a matrix of [popsize]x[N] will return

% a vector of size [popsize]x[1] containing the corresponding function values

% for the [N] trial vectors.

% GODLIKE stands for "Global Optimum Determination by Linking and

% Interchanging Kindred Evaluators". It basically links 4 different % optimization algorithms and interchanges

their findings iteratively, in

% an attempt to negate the various problems all the other algorithms

% encounter.

% The default control parameters GODLIKE uses, are:

% AlgIters = 250 The number of iterations spent in each

% heuristic optimizer.

% Swapiters = 10 The number of iterations before the

% populations are swapped between the

% algorithms.

% PutbackIters = 25 The number if iterations before the best

% individual ever found is inserted back into

% one of the populations.

% These values, as well as additional options, may be given manually in

% any extra input variables. Additionally, GODLIKE may be called as

% GODLIKE(PROBLEM), where [PROBLEM] is a complete problem-structure

% constructed by HEURSET.

% Some optimizations require repeated evaluation of a function that takes

% in the order of hours per evaluation. In those cases, it might be

% convenient to interrupt the optimization

after a few days and use the

% results thus far obtained. Unfortunately, usually after you interrupt a

% function, its results are lost. For that reason, GODLIKE creates two

% global variables, [GODLIKE_bestind] and [GODLIKE_bestfval], which store

% the best individual and function value thus far found. When the

% optimization is interrupted, the intermediate results from the

% optmization are still available. Once an optimization has succesfully

% converged, these variables are deleted from the workspace again.

% See also heurset, diffevolve, genetic, swarm, simanneal.

% Author: Rody P.S. Oldenhuis

% Delft University of Technology

% E-mail: oldenhuis@dds.nl

% Last edited: 02/Feb/2009

%% initialize & parse input

% initial values

problem = heurset;

[func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...

convvalue, swapiters, algiters, putbackiters] = parseprob(problem);

% define temporary solution globals

global GODLIKE_bestind GODLIKE_bestfval

% common inputs

if (nargin >= 4)

func = varargin{1};

popsize = varargin{2};

lb = varargin{3};

ub = varargin{4};

end

% with additional options

if (nargin >= 5)

if ~isstruct(varargin{5})

options = heurset(varargin{5:end});

else

options = varargin{5};

end

[dum1, dum2, dum3, dum4, grace, display, maxfevals, convmethod,...

convvalue, swapiters, algiters, putbackiters] = parseprob(options);

end

% if given a full problem structure

if (nargin == 1)

problem = varargin{1};

% errortrap

if ~isstruct(problem)

error('PROBLEM should be a structure. Type `help heurset` for details.')

end

[func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...

convvalue, swapiters, algiters, putbackiters] = parseprob(problem);

end

% initialize convergence method, using

default values if omitted

if strcmpi(convmethod, 'exhaustive')

convergence = 1;

maxiterinpop = convvalue;

maxiters = inf;

elseif strcmpi(convmethod, 'MaxIters')

convergence = 2;

maxiterinpop = inf;

maxiters = convvalue;

elseif strcmpi(convmethod, 'achieveFval')

convergence = 3;

%errortrap

if isempty(convvalue)

error('Please define function value

to achieve.')

end

maxiterinpop = inf;

maxiters = inf;

else

convergence = 1;

maxiterinpop = convvalue;

end

% problem dimensions

dims = size(lb, 2);

% errortraps

if ( (size(lb, 2) ~= 1 || size(ub, 2) ~= 1) &&...

(size(lb, 2) ~= dims || size(ub, 2) ~= dims) )

error('Upper- and lower boundaries must be the same size.')

end

if (popsize < 3)

error('GODLIKE needs a population size of at least 3.')

end

%% initial values

% iteration-zero values

firsttime = true;

iters = 0;

evals = 0;

globalbestfit = inf;

improvement = maxiterinpop;

converging1 = false;

converging2 = false;

% set popsize to be divisible by 4

while (rem(popsize, 4) ~= 0)

popsize = popsize + 1;

end

% initialize populations

pop = repmat(ub, popsize, 1) - repmat(ub-lb, popsize, 1).* rand(popsize, dims);

% initial fitnesses

fits = inf(popsize, 1);

% update options structure

problem.costfun = func;

problem.lb = lb;

problem.ub = ub;

problem.popsize = popsize/4;

problem.grace = 0;

problem.display = false;

problem.conv.method = 'MaxIters';

problem.conv.value = algiters; % initialize different evaluators probGA = problem; probDE = problem; probPS = problem; probSA = problem; % Display strings if display fprintf(1, ['\nNote that when GODLIKE is interrupted, the best solution and function\n',... 'value are saved in global variables GODLIKE_bestind and GODLIKE_bestfval.\n\n']); fprintf(1, 'Optimizing with '); switch convergence case 1 fprintf(1, 'exhaustive convergence criterion...\n'); case 2 fprintf(1, 'maximum iterations convergence criterion...\n'); case 3 fprintf(1, 'goal-attain convergence criterion...\n'); end fprintf(1, 'Current best solution has cost of ------------- '); overfevals1 = '\n\nCould not find better solution within given maximum'; overfevals2 = '\nof function evaluations. Increase MaxFevals option.\n\n'; fitbackspace = '\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b'; converge2bck = [fitbackspace, '\b\b\b\b\b\b\b']; converge1bck = [converge2bck, '\b\b\b\b']; convergestr = ', possibly converging: '; itersstr = ', iterations left: '; end %% optimization loop % loop while one of the convergence criteria is not met while (improvement > 0 && iters <= maxiters) % if [iters] exceeds [swapiters], swap the populations randomly if (mod(iters, swapiters) == 0) % randomize the row-ordering [ignore, inds] = sort(rand(1, popsize)); pop = pop(inds, :); % assign new populations to the different evaluators probGA.GA.pop = pop(1:4:end, :); probDE.DE.pop = pop(2:4:end, :); probPS.PS.pop = pop(3:4:end, :);

probSA.PS.pop = pop(4:4:end, :); % probSA.SA.pop = pop(4:4:end, :); % Simulated Annealing is horribly slow!! end % use all heuristic evaluators [pop(1:4:end, :), fits(1:4:end), evals1] = genetic([], probGA); [pop(2:4:end, :), fits(2:4:end), evals2] = diffevolve([], probDE); [pop(3:4:end, :), fits(3:4:end), evals3] = swarm([], probPS); [pop(4:4:end, :), fits(4:4:end), evals4] = swarm([], probSA); % update function evaluations evals = evals + evals1 + evals2 + evals3 + evals4; % find best and worst fitnesses [bestfit , bind] = min(fits); [worstfit, wind] = max(fits); % improving solutions if (bestfit < globalbestfit) % new best individual improvement = maxiterinpop; globalbestfit = bestfit; globalbestind = pop(bind, :); % store also in globals GODLIKE_bestind = globalbestind; GODLIKE_bestfval = globalbestfit; % display improving solution if display if converging1 fprintf(1, converge1bck); pause(0.05) end if converging2 fprintf(1, converge2bck); pause(0.05) end if (globalbestfit < 0) fprintf('\b') end fprintf(1, fitbackspace); fprintf(1, '%1.8e', globalbestfit); converging1 = false; firsttime = true; pause(0.05) end end % Check for convergence if (evals > maxfevals)

% maximum allowed function

evaluations has been superceded

fprintf(overfevals1);

fprintf(overfevals2);

break

end

switch convergence

% exhaustive

case 1

if (bestfit >= globalbestfit) && (globalbestfit ~= inf)

if (improvement <= 100) && display

if ~converging1

fprintf(1, convergestr);

fprintf(1, '%3.0f', improvement-1);

converging1 = true;

pause(0.05)

else

fprintf(1,

'\b\b\b%3.0f', improvement-1);

pause(0.05)

end

improvement = improvement - 1;

end

% max iterations

case 2

if display && (maxiters-iters < 100)

if firsttime

fprintf(1, itersstr);

fprintf(1, '%3.0f', maxiters-iters);

pause(0.05)

else

fprintf(1, '\b\b\b%3.0f', maxiters-iters);

pause(0.05)

end

firsttime = false;

converging2 = true;

end

iters = iters + 1;

% goal-attain

case 3

if (globalbestfit <= convvalue)

break;

end

% replace worst individual with best individual

if (mod(iters, putbackiters) == 0)

pop(wind, :) = globalbestind;

end

%% (pre-) end values

% if a solution has been found

if isfinite(globalbestfit)

sol = globalbestind;

fval = globalbestfit;

if display

fprintf(1, '\nGODLIKE-algorithm has Converged.\n');

end

% no feasible value might be found

else

fprintf(1,'\n');

warning('godlike:no_solution',...

'GODLIKE was unable to find any solution that gave finite values.')

fval = globalbestfit;

sol = NaN;

end

%% Grace function evaluations

if (grace > 0)

% display progress

if display

fprintf(1, 'Performing direct-search...');

pause(0.05)

end

% perform direct search

options = optimset('TolX', eps,

'MaxFunEvals', grace, 'TolFun', ...

eps, 'MaxIter', 1e4, 'Display', 'off');

[soltry, fvaltry] = fminsearch(func, sol, options);

% enforce boundaries

if ~any(soltry >= ub | soltry <= lb)

sol = soltry;

fval = fvaltry;

end

evals = evals + grace;

end

%% finalize

% display progress

if display

fprintf(1, 'All done.\n\n');

pause(0.05)

end

% clear temp globals

clear global GODLIKE_bestfval GODLIKE_bestind