智能优化算法源代码

更新时间：2024-05-23 11:29:01 阅读量：综合文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

人工蚂蚁算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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; 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;?stindex的蚂蚁不进行全局转移 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); %%%%%%%%%%%***************************************************

?_T=func1;

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

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

?_T=(X1_T-1).^2+(X2_T-2.2).^2+1;

%func2 =

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

?_B=func2;

?_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;

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; 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 2

plot(X(BestIndex,1),X(BestIndex,2),'rs','MarkerFaceColor','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;ìHO循环结束

[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),valfx));

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; 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; 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; 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)

Tfault 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 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 if

feval(func,x(i,:))

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

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,:); if

feval(func,x(i,:))

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

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

p(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 end

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

% if

Foxhole(pg,D)

%-----------开始进行免疫－－－－－－－－－－－－－－－－ 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 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 end end end 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)

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_incorrect_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 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 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 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 end 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 end

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 \ % 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 end end 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 \,...

'required type is %s, got %s.\\n',...

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

warning('heurset:incorrectoption', ... 'Ignoring option `%s`', num2str(option)) end 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 \

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 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 end

% replace worst individual with best individual if (mod(iters, putbackiters) == 0) pop(wind, :) =

globalbestind;

end 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 end

% parser function to easily parse the input arguments

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

convvalue, swapiters, algiters, putbackiters] = 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; swapiters =

problem.GODLIKE.swapiters; algiters =

problem.GODLIKE.algiters; putbackiters =

problem.GODLIKE.putbackiters; end

粒子群算法

% pso_Trelea_vectorized.m

% a generic particle swarm optimizer % to find the minimum or maximum of any

% MISO matlab function %

% Implements Common, Trelea type 1 and 2, and Clerc's class 1\ % also automatically try to track to a changing environment (with varied % success - BKB 3/18/05) %

% This vectorized version removes the for loop associated with particle % number. It also *requires* that the cost function have a single input % that represents all dimensions of search (i.e., for a function that has 2

% inputs then make a wrapper that

passes a matrix of ps x 2 as a single% variable) %

% Usage:

% [optOUT]=PSO(functname,D) % or:

% [optOUT,tr,te]=...

PSO(functname,D,mv,VarRange,minmax,PSOparams,plotfcn,PSOseedValue) %

% Inputs:

% functname - string of matlab function to optimize

% D - # of inputs to the function (dimension of problem) %

% Optional Inputs:

% mv - max particle velocity,

either a scalar or a vector of length D

% (this allows each component to have it's own max velocity),

% default = 4, set if not input or input as NaN %

% VarRange - matrix of ranges for each input variable,

% default -100 to 100, of form: % [ min1 max1 % min2 max2 % ...

% minD maxD ] %

% minmax = 0, funct minimized (default)

% = 1, funct maximized

% = 2, funct is targeted to P(12) (minimizes distance to errgoal) %

% PSOparams - PSO parameters

% P(1) - Epochs between updating display, default = 100. if 0, % no display

% P(2) - Maximum number of

iterations (epochs) to train, default = 2000.

% P(3) - population size, default = 24 %

% P(4) - acceleration const 1 (local best influence), default = 2

% P(5) - acceleration const 2 (global best influence), default = 2 % P(6) - Initial inertia weight, default = 0.9

% P(7) - Final inertia weight, default = 0.4

% P(8) - Epoch when inertial

weight at final value, default = 1500 % P(9)- minimum global error gradient,

% if abs(Gbest(i+1)-Gbest(i)) < gradient over

% certain length of epochs, terminate run, default = 1e-25

% P(10)- epochs before error gradient criterion terminates run, % default = 150, if the SSE does not change over 250 epochs

% then exit

% P(11)- error goal, if NaN then unconstrained min or max, default=NaN % P(12)- type flag (which kind of PSO to use)

% 0 = Common PSO w/intertia (default)

% 1,2 = Trelea types 1,2

% 3 = Clerc's Constricted PSO, Type 1\

% P(13)- PSOseed, default=0 % = 0 for initial positions all random

% = 1 for initial particles as user input %

% plotfcn - optional name of plotting function, default 'goplotpso',

% make your own and put here %

% PSOseedValue - initial particle position, depends on P(13), must be

% set if P(13) is 1 or 2, not used for P(13)=0, needs to % be nXm where n<=ps, and m<=D

% If n

% on Varrange % Outputs:

% optOUT - optimal inputs and

associated min/max output of function, of form:

% [ bestin1 % bestin2 % ... % bestinD % bestOUT ] %

% Optional Outputs:

% tr - Gbest at every iteration, traces flight of swarm

% te - epochs to train, returned as a vector 1:endepoch %

% Example: out=pso('Ackley')

% Brian Birge % Rev 3.3 % 2/18/06

function

[OUT,varargout]=PSO(functname,D,varargin)

rand('state',sum(100*clock)); if nargin < 1

error('Not enough arguments.'); end

% PSO PARAMETERS

if nargin == 1 % only specified functname and D D=2

VRmin=ones(D,1)*-100; VRmax=ones(D,1)*100;

VR=[VRmin,VRmax]; minmax = 0; P = []; mv = 4;

plotfcn='goplotpso';

elseif nargin == 3 % specified functname, D, and mv

VRmin=ones(D,1)*-100; VRmax=ones(D,1)*100; VR=[VRmin,VRmax]; minmax = 0;

mv=varargin{1}; if isnan(mv) mv=4; end

P = [];

plotfcn='goplotpso';

elseif nargin == 4 % specified functname, D, mv, Varrange mv=varargin{1}; if isnan(mv) mv=4; end

VR=varargin{2}; minmax = 0; P = [];

plotfcn='goplotpso';

elseif nargin == 5 % Functname, D, mv, Varrange, and minmax mv=varargin{1}; if isnan(mv) mv=4; end

VR=varargin{2};

minmax=varargin{3}; P = [];

plotfcn='goplotpso';

elseif nargin == 6 % Functname, D, mv, Varrange, minmax, and psoparams mv=varargin{1}; if isnan(mv) mv=4; end

VR=varargin{2};

minmax=varargin{3};

P = varargin{4}; % psoparams

plotfcn='goplotpso';

elseif nargin == 7 % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn

mv=varargin{1}; if isnan(mv) mv=4; end

VR=varargin{2};

minmax=varargin{3};

P = varargin{4}; % psoparams plotfcn = varargin{5};

elseif nargin == 8 % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn, PSOseedValue mv=varargin{1}; if isnan(mv) mv=4; end

VR=varargin{2};

minmax=varargin{3};

P = varargin{4}; % psoparams plotfcn = varargin{5}; PSOseedValue = varargin{6}; else

error('Wrong # of input arguments.'); end

% sets up default pso params

Pdef = [100 2000 100 2 2 0.9 0.4 1500 1e-25 250 NaN 0 0]; Plen = length(P);

P = [P,Pdef(Plen+1:end)];

df = P(1); me = P(2); ps = P(3); ac1 = P(4); ac2 = P(5); iw1 = P(6); iw2 = P(7); iwe = P(8); ergrd = P(9); ergrdep = P(10); errgoal = P(11);

trelea = P(12); PSOseed = P(13);

% used with trainpso, for neural net training

if strcmp(functname,'pso_neteval') net = evalin('caller','net'); Pd = evalin('caller','Pd'); Tl = evalin('caller','Tl'); Ai = evalin('caller','Ai'); Q = evalin('caller','Q'); TS = evalin('caller','TS'); end

% error checking

if ((minmax==2) & isnan(errgoal)) error('minmax= 2, errgoal= NaN: choose an error goal or set minmax to 0 or 1'); end

if ( (PSOseed==1) &

~exist('PSOseedValue') )

error('PSOseed flag set but no PSOseedValue was input'); end

if exist('PSOseedValue')

tmpsz=size(PSOseedValue); if D < tmpsz(2)

error('PSOseedValue column size must be D or less'); end

if ps < tmpsz(1)

error('PSOseedValue row length must be # of particles or less'); end end

% set plotting flag if (P(1))~=0 plotflg=1; else

plotflg=0;