二、蜻蜓算法简介

1蜻蜓仿生寻优算法

蜻蜓算法寻优计算主要思路通过模拟蜻蜓的捕食行为来实现(图1),该飞行(进化)寻优机制可以表示为蜻蜓群体分离、对齐、聚集、食物吸引与天敌驱散五个步骤。

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_算法

图1 蜻蜓算法寻优计算主要思路

(1)分离。表示避免蜻蜓个体距离太近降低寻优效率:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_线性代数_02

式中Ek为蜻蜓k的分离度;D表示当前蜻蜓的位置;Dj代表第j个邻近蜻蜓的位置;J代表群体中第k个蜻蜓的邻近蜻蜓的数量。

(2)对齐。表示某个蜻蜓与其邻近蜻蜓个体速度的相同程度:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_线性代数_03

式中Uk为蜻蜓k的对齐度;Vj代表第j个邻近蜻蜓的飞行速度;其他符号意义同前。

(3)聚集。表示某个蜻蜓具有朝着其附近较优蜻蜓个体靠近的趋势:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_线性代数_04

式中Bk为蜻蜓k的聚集度;其他符号意义同前。

(4)食物吸引。食物是指某一次迭代计算中最优个体的位置:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_迭代_05

式中Sk为蜻蜓k的食物吸引度;D+代表当前计算中最优蜻蜓的位置。

(5)天敌驱散。目的是使个体尽可能远离最差的蜻蜓个体,提高寻优计算效率:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_算法_06

式中Tk为蜻蜓k的天敌驱散度;D-代表当前计算中最差蜻蜓的位置。

蜻蜓个体k飞行位置更新步长:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_取值_07

蜻蜓飞行位置更新:

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_线性代数_08

式中e、u、b、s、t分别表示分离度、对齐度、聚集度、食物吸引度、天敌驱散度的影响系数;β代表惯性系数;l代表反复迭代计数下标;其他符号意义同前。

2 蜻蜓算法实现

2.1 算法调整

该研究对蜻蜓算法的寻优求解机制进行如下调整。

(1)引入初始种群个体优化机制。对初始种群中随机选择的某个蜻蜓个体,通过初步寻优计算得到一个可行解,优化后得到新个体替换原有蜻蜓个体,可以保证初始种群中至少有一个可行解。

(2)增加局部合理性判定机制。在蜻蜓算法中加入了合理性审查算子,该算子通过遍历单个蜻蜓的计算维度,保留飞行后合理的局部,剔除飞行后发生不合理的局部位置,以提高算法的寻优计算效率。

2.2 计算流程

适应性调整后的蜻蜓算法的计算流程包括12个步骤。

Step1:初始化计算参数,考虑梯级电站短期优化调度模型求解规模设定蜻蜓算法种群个数N=40,迭代计算次数MAXiter=200。

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_线性代数_09

Step3:随机选择一个初始蜻蜓个体,采用POA算法,基于单站优化原则,对选定蜻蜓个体进行寻优计算,确保初始种群中至少有一个蜻蜓个体为可行解。并令迭代次数iter=1,蜻蜓维度r=0。

Step4:判断蜻蜓个体Xk是否存在邻居。如果有邻居,进入Step5,否则进入Step6。

Step5:依据式(6)-(10)进行群体进化计算:分离度Sr、对齐度Ar、聚集度Cr、食物吸引度Fr、天敌驱散度Er。依据式(11)计算蜻蜓个体位置更新步长,进入Step7。

Step6:该蜻蜓没有邻居,采用Xr=e×(ZMAX(i,t)-ZMIN(i,t))进行蜻蜓个体随机飞行,进入Step7。

Step7:根据当前维度位置更新后计算得到新的个体,计算当前维度更新前后蜻蜓个体目标函数值func1,func2。如果判断func2>func1,说明该维度的位置变化对目标函数计算有利,则保留当前维度的位置更新,ΔXkt+1=ΔXkt+1。若func2<func1,说明当前维度位置变化对目标函数不利,则ΔXkt+1=0。

Step8:使蜻蜓个体计算维度加1,即r=r+1。判断是否完成蜻蜓个体全部维度的更新,如果完成全部维度的更新r=rmax,则进入Step9,否则进入Step5循环计算。

Step9:对蜻蜓个体Xk进行位置更新,Xt+1=Xt+ΔXt+1。

Step10:令k=k+1,进入新一个蜻蜓个体的优化计算,如果k>=N,则全部蜻蜓个体位置更新计算结束,进入Step11。若k<N,则进入Stpe4进行循环计算。

Step11:当前迭代次数中,找出最优蜻蜓个体,与食物位置相比较,根据是否优于食物位置目标函数值,更新食物位置,令iter=iter+1。若iter>=MAXiter,则进入Step12,否则进入Step4。

Step12:返回当前循环迭代中食物位置作为最优计算结果,计算结束。

三、部分源代码

%__________________________________________
% fobj = @YourCostFunction
% dim = 待优化参数个数
% Max_iteration =最大迭代次数
% SearchAgents_no = 蜻蜓数量
% lb=[lb1,lb2,...,lbn] 参数取值下限
% ub=[ub1,ub2,...,ubn] 参数取值上限

% To run DA: [Best_score,Best_pos,cg_curve]=DA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj)
%__________________________________________

clear all
clc

SearchAgents_no=40; % 蜻蜓数量

Function_name='F1'; % Name of the test function that can be from F1 to F23 (Table 1,2,3 in the paper)

Max_iteration=500; % 最大迭代次数

% Load details of the selected benchmark function
[lb,ub,dim,fobj]=Get_Functions_details(Function_name);%函数相关参数初始化

[Best_score,Best_pos,cg_curve]=DA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj);

figure('Position',[400 400 560 190])

%Draw search space
subplot(1,2,1);
func_plot(Function_name);
title('Test function')
xlabel('x_1');
ylabel('x_2');
zlabel([Function_name,'( x_1 , x_2 )'])
grid off

%Draw objective space
subplot(1,2,2);
semilogy(cg_curve,'Color','r')
title('Convergence curve')
xlabel('Iteration');
ylabel('Best score obtained so far');

axis tight
grid off
box on
legend('DA')

display(['The best solution obtained by DA is : ', num2str(Best_pos')]);
display(['The best optimal value of the objective funciton found by DA is : ', num2str(Best_score)]);

% dim = 待优化参数个数
% Max_iteration =最大迭代次数
% SearchAgents_no = 蜻蜓数量
% lb=[lb1,lb2,...,lbn] 参数取值下限
% ub=[ub1,ub2,...,ubn] 参数取值上限
% To run DA: [Best_score,Best_pos,cg_curve]=DA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj)
%__________________________________________

function [Best_score,Best_pos,cg_curve]=DA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj) %function [输出变量] = 函数名称(输入变量)

display('DA is optimizing your problem');
cg_curve=zeros(1,Max_iteration); %返回一个1 x Max_iteration的零矩阵

if size(ub,2)==1 %size(ub,2)返回矩阵ub的列数,这个条件表示ub是一个数,即待优化参数个数为1
ub=ones(1,dim)*ub; %ones(1,dim)产生1xdim的全1矩阵
lb=ones(1,dim)*lb;
end

%初始化蜻蜓邻里半径
r=(ub-lb)/10; % 1*10 20
Delta_max=(ub-lb)/10; % 1*10 20

Food_fitness=inf;%正无穷
Food_pos=zeros(dim,1);%10*1 0

Enemy_fitness=-inf;%负无穷
Enemy_pos=zeros(dim,1);%10*1 0

X=initialization(SearchAgents_no,dim,ub,lb);%10*40 [-100,100] 随机数
Fitness=zeros(1,SearchAgents_no); %1*40 0

DeltaX=initialization(SearchAgents_no,dim,ub,lb);%10*40 [-100,100]

for iter=1:Max_iteration
r=(ub-lb)/4+((ub-lb)*(iter/Max_iteration)*2);%迭代次数越大。半径越大初始1*10 50.8

w=0.9-iter*((0.9-0.4)/Max_iteration);%不断减小1*1 0.899

my_c=0.1-iter*((0.1-0)/(Max_iteration/2));% 不断减小 1*1 0.0996
if my_c<0
my_c=0;
end
s=2*rand*my_c; % 分离度 0.0013
a=2*rand*my_c; % 对齐度 0.1884
c=2*rand*my_c; %内聚度 0.1791
f=2*rand; % 食物吸引力 0.8826
e=my_c; %敌排斥力 0.0996
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%找到食物和天敌
for i=1:SearchAgents_no %首先计算所有目标值
Fitness(1,i)=fobj(X(:,i)'); %X(:,i)'是取矩阵X的所有行的第i列并共轭转置
if Fitness(1,i)<Food_fitness %寻找每次迭代的最小值
Food_fitness=Fitness(1,i);%1.2728*10^4
Food_pos=X(:,i);
end

if Fitness(1,i)>Enemy_fitness %寻找每次迭代的最大值
if all(X(:,i)<ub') && all( X(:,i)>lb')
Enemy_fitness=Fitness(1,i);%5.6813*10^4
Enemy_pos=X(:,i);
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%找到每只蜻蜓的邻居
for i=1:SearchAgents_no
index=0;
neighbours_no=0;

clear Neighbours_DeltaX
clear Neighbours_X
%找到相邻邻居
for j=1:SearchAgents_no
Dist2Enemy=distance(X(:,i),X(:,j));%计算欧氏距离
if (all(Dist2Enemy<=r) && all(Dist2Enemy~=0))
index=index+1;%邻居序号
neighbours_no=neighbours_no+1;%邻居数量
Neighbours_DeltaX(:,index)=DeltaX(:,j);
Neighbours_X(:,index)=X(:,j);
end
end

% 分离 - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eq. (3.1)
S=zeros(dim,1);
if neighbours_no>1
for k=1:neighbours_no
S=S+(Neighbours_X(:,k)-X(:,i));
end
S=-S;
else
S=zeros(dim,1);
end

% 对齐%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eq. (3.2)
if neighbours_no>1
A=(sum(Neighbours_DeltaX')')/neighbours_no;
else
A=DeltaX(:,i);
end

% 内聚%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eq. (3.3)
if neighbours_no>1
C_temp=(sum(Neighbours_X')')/neighbours_no;
else
C_temp=X(:,i);
end

C=C_temp-X(:,i);

% 靠近食物%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eq. (3.4)
Dist2Food=distance(X(:,i),Food_pos(:,1));
if all(Dist2Food<=r)
F=Food_pos-X(:,i);
else
F=0;
end

% 远离天敌%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eq. (3.5)
Dist2Enemy=distance(X(:,i),Enemy_pos(:,1));
if all(Dist2Enemy<=r)
Enemy=Enemy_pos+X(:,i);
else
Enemy=zeros(dim,1);
end

for tt=1:dim
if X(tt,i)>ub(tt)%大于上限
X(tt,i)=lb(tt);
DeltaX(tt,i)=rand;
end
if X(tt,i)<lb(tt)
X(tt,i)=ub(tt);
DeltaX(tt,i)=rand;
end
end

四、运行结果

【优化算法】蜻蜓算法(DA)【含Matlab源码 1310期】_迭代_10

五、matlab版本及参考文献

1 matlab版本

2014a

2 参考文献

[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例(第2版)[M].电子工业出版社,2016.

[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社,2017.