多目标蜻蜓优化算法（Matlab代码实现）

目录

​​1 蜻蜓算法概述​​

​​2 蜻蜓算法数学模型​​

​​2.1 分离​​

​​2.2 对齐​​

​​2.3 聚集​​

​​2.4 食物吸引​​

​​2.5 天敌驱散​​

​​3 完整matlab代码实现 ​​

---

 此算法在电力系统中也有很多应用 

1 蜻蜓算法概述

自然充满了执行不同任务的社会行为。尽管所有个体和集体行为的最终目标是生存，但生物在群体、畜群、学校、殖民地和羊群中合作和互动有几个原因：狩猎、防御、导航和觅食。例如，狼群拥有组织最完善的狩猎社交互动之一。狼倾向于遵循社会领导以不同的步骤捕猎猎物：追逐猎物、盘旋猎物、骚扰猎物和攻击猎物 。集体防御的一个例子是海洋中的鱼群。数以千计的鱼类形成了一个鱼群，并通过相互警告来避开捕食者，使得捕食者的捕食变得非常困难。

蜻蜓算法（DA）。 DA算法的主要灵感来源于自然界蜻蜓静态和动态的成群行为。优化、探索和开发的两个基本阶段是通过对蜻蜓在导航、寻找食物和在动态或统计上成群时避开敌人的社会互动进行建模来设计的。接下来文章还考虑了 DA 的二进制和多目标版本的方法，分别称为二进制 DA (BDA) 和多目标 (MODA)。所提出的算法由几个数学测试函数和一个定性和定量的真实案例研究进行了基准测试。该算法能够改善给定问题的初始随机种群，收敛到全局最优，并提供非常有竞争力的结果。 MODA 的结果还表明，该算法倾向于为多目标问题找到具有高度均匀分布的帕累托最优解的非常精确的近似值。潜艇螺旋桨设计问题得到的一组设计证明了MOD的优点。

2 蜻蜓算法数学模型

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

2.1 分离

表示避免蜻蜓个体距离太近降低寻优效率:

                       

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

2.2 对齐

表示某个蜻蜓与其邻近蜻蜓个体速度的相同程度:

                        

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

2.3 聚集

表示某个蜻蜓具有朝着其附近较优蜻蜓个体靠近的趋势:

                            

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

2.4 食物吸引

食物是指某一次迭代计算中最优个体的位置:

                           

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

2.5 天敌驱散

目的是使个体尽可能远离最差的蜻蜓个体,提高寻优计算效率:

                         

式中Tk为蜻蜓k的天敌驱散度;D-代表当前计算中最差蜻蜓的位置。
蜻蜓个体k飞行位置更新步长:

蜻蜓飞行位置更新:

                  

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

3 完整matlab代码实现 

关注微信公众号：​​电力系统与算法之美​​

%===============多目标蜻蜓优化算法MODA==============
%%%欢迎关注公众号：电力系统与算法之美%%%%%
%%%回复关键字：多目标蜻蜓即可获得全部代码%%%

clc;
clear;
close all;

%% 多目标目标函数
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

%% 初始化多目标蜻蜓优化算法参数
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;             % 分离权重
        a=my_c;             %对齐权重
        c=my_c;             % 聚集权重
        f=2*rand;           % 食物吸引权重
        e=my_c;             % 天敌驱散权重
    else
        s=my_c/iter;       
        a=my_c/iter;        
        c=my_c/iter;        
        f=2*rand;           
        e=my_c/iter;        
    end
    %% 首先计算所有的目标值
    for i=1:N 
        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);
    
    %% 提高复盖率的食物（选择种群最少地区的归档成员以提高复盖率的食物）
    index=RouletteWheelSelection(1./Archive_mem_ranks);
    if index==-1
        index=1;
    end
    Food_fitness=Archive_F(index,:);
    Food_pos=Archive_X(index,:)';
       
    %% （提高复盖率）选择种群最多的地区的档案成员作为敌人来提高复盖率
    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
        %% 找到附近的解决方案
        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
        
        %% 分离
        
        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
        
        %% 对齐
        
        if neighbours_no>1
            A=(sum(Neighbours_V')')/neighbours_no;
        else
            A=DeltaX(:,i);
        end
        
        %% 聚集
       
        if neighbours_no>1
            C_temp=(sum(Neighbours_X')')/neighbours_no;
        else
            C_temp=X(:,i);
        end
        
        C=C_temp-X(:,i);
        
        %% 食物吸引
       
        Dist2Attraction=distance(X(:,i),Food_pos(:,1));
        if all(Dist2Attraction<=r)
            F=Food_pos-X(:,i);
            iter;
        else
            F=0;
        end
        
        %% 天敌驱散
       
        Dist=distance(X(:,i),Enemy_pos(:,1));
        if all(Dist<=r)
            E=Enemy_pos+X(:,i);
        else
            E=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
        
        
        if any(Dist2Attraction>r)
            if neighbours_no>1
                for j=1:dim
                    DeltaX(j,i)=w*DeltaX(j,i)+rand*A(j,1)+rand*C(j,1)+rand*S(j,1);
                    if DeltaX(j,i)>V_max
                        DeltaX(j,i)=V_max;
                    end
                    if DeltaX(j,i)<-V_max
                        DeltaX(j,i)=-V_max;
                    end
                    X(j,i)=X(j,i)+DeltaX(j,i);
                end
                
            else
                X(:,i)=X(:,i)+Levy(dim)'.*X(:,i);
                DeltaX(:,i)=0;
            end
        else    
            for j=1:dim
                DeltaX(j,i)=s*S(j,1)+a*A(j,1)+c*C(j,1)+f*F(j,1)+e*E(j,1) + w*DeltaX(j,i);
                if DeltaX(j,i)>V_max
                    DeltaX(j,i)=V_max;
                end
                if DeltaX(j,i)<-V_max
                    DeltaX(j,i)=-V_max;
                end
                X(j,i)=X(j,i)+DeltaX(j,i);
            end
        end
        
        Flag4ub=X(:,i)>ub';
        Flag4lb=X(:,i)<lb';
        X(:,i)=(X(:,i).*(~(Flag4ub+Flag4lb)))+ub'.*Flag4ub+lb'.*Flag4lb;
        
    end
    
    display(['在迭代时 ', num2str(iter), '有 ', num2str(Archive_member_no), ' 非支配解决方案']);
end


figure

Draw_ZDT1();

hold on
if obj_no==2
    plot(Archive_F(:,1),Archive_F(:,2),'ko','MarkerSize',8,'markerfacecolor','k');
else
    plot3(Archive_F(:,1),Archive_F(:,2),Archive_F(:,3),'ko','MarkerSize',8,'markerfacecolor','k');
end
legend('True PF','Obtained PF');
title('MODA');