【优化算法】多目标蜻蜓优化算法（MODA）【含Matlab源码 1350期】

二、蜻蜓算法简介

1蜻蜓仿生寻优算法

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

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

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

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

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

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

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

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

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

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

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

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

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

蜻蜓飞行位置更新:

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

2 蜻蜓算法实现

2.1 算法调整

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

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

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

2.2 计算流程

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

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

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:返回当前循环迭代中食物位置作为最优计算结果,计算结束。

三、部分源代码

%___________________________________________________________________%
%  Multi-Objective Dragonfly Algorithm (MODA) source codes demo     %
%                           version 1.0                             %
%                                                                   %
           %
%___________________________________________________________________%

clc;
clear;
close all;

% Change these details with respect to your problem%%%%%%%%%%%%%%
ObjectiveFunction=@ZDT1;
dim=5;
lb=0;
ub=1;
obj_no=2;

if size(ub,2)==1
    ub=ones(1,dim)*ub;
    lb=ones(1,dim)*lb;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initial parameters of the MODA algorithm
max_iter=100;
N=100;
ArchiveMaxSize=100;

Archive_X=zeros(100,dim);
Archive_F=ones(100,obj_no)*inf;

Archive_member_no=0;

r=(ub-lb)/2;
V_max=(ub(1)-lb(1))/10;

Food_fitness=inf*ones(1,obj_no);
Food_pos=zeros(dim,1);

Enemy_fitness=-inf*ones(1,obj_no);
Enemy_pos=zeros(dim,1);

X=initialization(N,dim,ub,lb);
fitness=zeros(N,2);

DeltaX=initialization(N,dim,ub,lb);
iter=0;

position_history=zeros(N,max_iter,dim);

for iter=1:max_iter
    
    r=(ub-lb)/4+((ub-lb)*(iter/max_iter)*2);
    
    w=0.9-iter*((0.9-0.2)/max_iter);
    
    my_c=0.1-iter*((0.1-0)/(max_iter/2));
    if my_c<0
        my_c=0;
    end
    
    if iter<(3*max_iter/4)
        s=my_c;             % Seperation weight
        a=my_c;             % Alignment weight
        c=my_c;             % Cohesion weight
        f=2*rand;           % Food attraction weight
        e=my_c;             % Enemy distraction weight
    else
        s=my_c/iter;        % Seperation weight
        a=my_c/iter;        % Alignment weight
        c=my_c/iter;        % Cohesion weight
        f=2*rand;           % Food attraction weight
        e=my_c/iter;        % Enemy distraction weight
    end
    
    for i=1:N %Calculate all the objective values first
        Particles_F(i,:)=ObjectiveFunction(X(:,i)');
        if dominates(Particles_F(i,:),Food_fitness)
            Food_fitness=Particles_F(i,:);
            Food_pos=X(:,i);
        end
        
        if dominates(Enemy_fitness,Particles_F(i,:))
            if all(X(:,i)<ub') && all( X(:,i)>lb')
                Enemy_fitness=Particles_F(i,:);
                Enemy_pos=X(:,i);
            end
        end
    end
    
    [Archive_X, Archive_F, Archive_member_no]=UpdateArchive(Archive_X, Archive_F, X, Particles_F, Archive_member_no);
    
    if Archive_member_no>ArchiveMaxSize
        Archive_mem_ranks=RankingProcess(Archive_F, ArchiveMaxSize, obj_no);
        [Archive_X, Archive_F, Archive_mem_ranks, Archive_member_no]=HandleFullArchive(Archive_X, Archive_F, Archive_member_no, Archive_mem_ranks, ArchiveMaxSize);
    else
        Archive_mem_ranks=RankingProcess(Archive_F, ArchiveMaxSize, obj_no);
    end
    
    Archive_mem_ranks=RankingProcess(Archive_F, ArchiveMaxSize, obj_no);
    
    % Chose the archive member in the least population area as foods
    % to improve coverage
    index=RouletteWheelSelection(1./Archive_mem_ranks);
    if index==-1
        index=1;
    end
    Food_fitness=Archive_F(index,:);
    Food_pos=Archive_X(index,:)';
       
    % Chose the archive member in the most population area as enemies
    % to improve coverage
    index=RouletteWheelSelection(Archive_mem_ranks);
    if index==-1
        index=1;
    end
    Enemy_fitness=Archive_F(index,:);
    Enemy_pos=Archive_X(index,:)';
    
    for i=1:N
        index=0;
        neighbours_no=0;
        
        clear Neighbours_V
        clear Neighbours_X
        % Find the neighbouring solutions
        for j=1:N
            Dist=distance(X(:,i),X(:,j));
            if (all(Dist<=r) && all(Dist~=0))
                index=index+1;
                neighbours_no=neighbours_no+1;
                Neighbours_V(:,index)=DeltaX(:,j);
                Neighbours_X(:,index)=X(:,j);
            end
        end
        
        % Seperation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 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
        
        % Alignment%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Eq. (3.2)
        if neighbours_no>1
            A=(sum(Neighbours_V')')/neighbours_no;
        else
            A=DeltaX(:,i);
        end

四、运行结果

五、matlab版本及参考文献

1 matlab版本

2014a

2 参考文献

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

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