闲鱼技术01 发表于 2022-8-19 14:28

非线性优化算法(二)多目标粒子群算法-matlab实现

许多优化问题的求解需要优化算法,因为各种语言都有很多工具包、工具箱、开源算法等等,典型如matlab全局优化工具箱。基本上是不需要再自己开发的。我自己研究和实现过两种,其他的就也是用工具箱。这里分享一下第二种。
粒子群算法是一种启发式多目标全局优化算法,启发自鸟群觅食的自然现象。启发式算法有很多,如遗传算法、蚁群算法、模拟退火算法、蒙特卡洛算法、免疫算法、鱼群算法等。这些算法是从不同的自然现象中受到启发而发展的,具有全局最优的搜索能力,不需要获得目标函数的梯度,针对复杂非线性系统比较有效。
(一)
通过多次应用,逐渐完善了一种自适应网格策略的多目标粒子群优化算法(AG-MOPSO)。基本原理如下:


首先,粒子群算法首先建立初始粒子种群。每个粒子表示在设计参数空间内的一个点,表示一种设计方案。每个粒子的目标函数值可以被计算作为粒子的适应度。在多目标优化问题中,粒子的相互优劣以支配的概念进行辨别。如果A粒子的所有目标函数值都优于B粒子,称A粒子支配B粒子,即A比B更优。如果A粒子只有部分目标函数优于B粒子,则两个粒子为非支配关系,并无优劣之分。粒子群中总有一部分粒子不被任何其他粒子支配,这部分粒子称为非劣集,也即粒子群的最优解集。通过对所有粒子进行目标函数计算可以得到初始粒子群的非劣解集。然后进入迭代。
在每一个迭代步中,每个粒子将获得速度


其中,X_i^k和V_i^k分别是第i个粒子在第k迭代步的位置和速度,omega是粒子保持其前一步速度的惯性权重系数,P_i^k和P_g^k分别是粒子的个体历史最优位置和整个粒子群的历史最优位置,c_1和c_2是指向最优位置的加速度因子,r_1和r_2是0和1之间的随机数。
为了增大开始时的全局搜索能力和后期的局部搜索能力,惯性权重设为随迭代步线性减小,


每一步迭代完成后,计算粒子的目标函数值,获得当前步的非劣解集,与原有非劣解集对比,得出整体最新的非劣解集。随着迭代的进行非劣解集的粒子数量逐渐增长。如果非劣解集的数量过大,容易导致局部收敛,同时也降低了计算效率。因此,多目标粒子群算法中,通常会采用自适应网格技术,对非劣解集进行删减。自适应网格是指在目标函数空间内,进行均匀网格划分,计算每个网格中的非劣解集数量,这个数量也成为网格密度。当非劣解集的数量超过给定值时,每个网格中随机删去一定数量的非劣解。第j个网格的删除数量为


其中N_k是当前非劣解集的总数量,N_lim是给定的非劣解集数量限制,d_j是第j个网格内的粒子数量。
为了计算每个粒子的速度,需要获得个体历史最优位置和整个粒子群的历史最优位置。
其中,整个粒子群的历史最优位置是从非劣解集中选择。为了避免聚集效应,粒子数量较少的网格内的粒子有较高的概率被选中


而个体历史最优位置,只需将当前粒子位置和个体历史最优位置进行比较,如果相互为非支配关系,则任选其一作为新的个体历史最优位置,否则选择两者中的支配者作为新的个体历史最优位置。
通常对限制条件的处理是,如果粒子不满足限制条件,则给予其一个极大的目标函数值。但针对一些比较复杂的问题,可能粒子不满足限制条件的情况有很大概率出现,如果采用这一方法,种群中满足条件的粒子很少,将极大降低最优解的搜索能力。因此,可以采用新的限制条件策略,即如果粒子移动后不满足限制条件,则重新产生速度,在当前位置重新移动,直至满足条件。这样做虽然会引入额外的计算量,但能够保证算法寻优的有效性。同时,为了避免部分情况下搜索陷入僵局,如果重新移动的次数过多,则给定粒子全局变异能力。变异的概率随不成功的重复次数的增加而增大,但最大不超过0.1。变异概率为


其中k_mu是基础变异概率,k_c是不成功的迭代次数,N_c是一个常数。全局变异策略也能够增加算法的全局搜索能力。
(二)matlab实现
为了便于使用,把粒子群法写成一个通用的matlab函数,方便调用。使用时首先要定义好优化问题,输入所需的参数(详见代码注释)。同时需要定义两个计算限制条件和目标函数的子函数供调用。
代码如下
function pop = PSO_MO_main(p_range,p_discrete,Ini_gene_method,p0_ini_para,NO)
%多目标粒子群算法主函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%输入参数:
%p_range,参数取值范围,用于随机产生初始种群,2行n列,第一行是参数下限,第二行是参数上限,
%p_discrete,共n项,表示设计变量是连续的(0),只能取整(1),或指定取值精度(n>0)
%Ini_gene_method初始种群的定义方式,0表示完全随机产生,1表示给定满足限制条件的初始种群,2计算中断后可继续计算(需要有输出文件pop.mat);
%p0_ini_para定义初值所需参数,初值定义方法为0或2时,本参数无意义,初值方法为1时,本参数为初始种群数据,共size_pop行n+NO列,需满足限制条件,前n列为设计变量值,后NO列为目标函数值;
%NO优化目标的个数;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%输出参数:
%pop优化过程详细信息,包含每一代种群的:
% 设计变量(size_pop行n列)、移动速度(size_pop行n列)、目标函数值(size_pop行NO列),
% 群体非劣解集(x行n列),群体非劣解的目标函数值(x行NO列),群体最优解(1*n),群体最优解目标函数值(1*NO)
% 个体最优解(size_pop行n列)、个体最优解目标函数值(size_pop行NO列)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%需定义以下子函数,以供调用
%限制条件J_lim=estimate_limit(p,p_range);满足限制条件返回1,否则返回0
%目标函数Obj=estimate_Obj(p);计算目标函数值,返回目标值1*NO,以目标值更小为优

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%关于限制条件的处理方式
%1——如不满足限制条件,重新移动,直至满足,为防止陷入局部死循环,设定不满足的次数越多,发生变异的概率越高,此种策略适用于限制条件较强的情况
%2——如不满足限制条件,目标函数取一极大值,以便完全阻止其成为非劣解,此种策略适用于限制条件较弱,即很容易满足的情况
limit_method=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%自适应网格尺寸:自适应网格技术是为了实现种群的多样性,将非劣解集进行网格划分,选择密度小的作为全局最优。同时删减非劣解集
grid_num=20; %对非劣解集进行网格划分的数量
noninf_num=200; %限定非劣解集的数量上限

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%算法参数定义
%种群数量
size_pop=100;
%最大代数
max_gen=200;
%加速度因子,一般取0-4
c1=1.495;
c2=1.495;
%初始惯性权重和最终惯性权重
w_start=0.8;
w_end=0.2;
%种群总变异概率
km=0.99/size_pop;


%变量维度
n=length(p_range(1,:));
%初始化输出参数,以结构体的形式储存优化过程信息。
pop=struct;

%定义速度最大值和最小值:参数一般取取值范围的(0.05-0.5)
V_max=(p_range(2,:)-p_range(1,:))*0.5;
V_min=-V_max;

%初始化种群
if Ini_gene_method==0
    p=zeros(size_pop,n);%参数
    V=zeros(size_pop,n);%速度
    f=zeros(size_pop,NO); %适应度
    for i=1:size_pop
      while 1==1
            p(i,:)=rand(1,n).*(p_range(2,:)-p_range(1,:))+p_range(1,:); %随机产生初值
            p(i,:)=get_discrete(p(i,:),p_discrete);
            lim=estimate_limit(p(i,:),p_range);
            if lim==1 %满足限制条件
                f(i,:)=estimate_Obj(p(i,:));
                break
            end
      end
      %V(i,:)=rand(1,n).*(V_max-V_min)+V_min; %可以随机定义初始速度,也可以定义为0
      gen_ini_pop=i   %实时显示初始化进度
    end
elseif Ini_gene_method==1
    p=p0_ini_para(:,1:n);
    f=p0_ini_para(:,n+1:n+NO);
    V=zeros(size_pop,n);
%   for i=1:size_pop
%       V(i,:)=rand(1,n).*(V_max-V_min)+V_min;   %可以随机定义初始速度,也可以定义为0
%   end
elseif Ini_gene_method==2
    %不产生初始种群,读取种群数据继续计算
end

%根据初始种群,计算非劣解集,种群最优解和个体最优解
if Ini_gene_method~=2
    %获得初始非劣解集
    =get_noninferior(p,f);
    %初始极值和适应度
    ibestp=p; %个体极值
    ibestf=f; %个体极值适应度
    %群体极值按照网络密度进行选择
    grid_inform=grid_noninferior(pnonf,grid_num);
    =select_pbest(pnonp,pnonf,grid_inform,size_pop);
else
    pop=load('pop.mat','pop'); %读取已有的数据,继续计算
    pop=pop.pop;
    s=length(pop);
    p=pop(s).pop;
    V=pop(s).V;
    f=pop(s).fitness;
    pnonp=pop(s).noninferp;
    pnonf=pop(s).noninferf;
    grid_inform=pop(s).grid_inform;
    pbestp=pop(s).pbestpop;
    pbestf=pop(s).pbestfitness;
    ibestp=pop(s).ibestpop;
    ibestf=pop(s).ibestfitness;
end


%开始迭代
if Ini_gene_method~=2
    t=1;
else
    t=s;
end

while t<=max_gen
    %记录当前代的结果并实时保存
    pop(t).pop=p;
    pop(t).V=V;
    pop(t).fitness=f;
    pop(t).noninferp=pnonp;
    pop(t).noninferf=pnonf;
    pop(t).pbestpop=pbestp;
    pop(t).pbestfitness=pbestf;
    pop(t).ibestpop=ibestp;
    pop(t).ibestfitness=ibestf;
    pop(t).grid_inform=grid_inform;
    save pop
   
    %画出当前代的非劣解集,暂只适用于2个或3个目标
    if NO==2
      clf
      scatter(pnonf(:,1),pnonf(:,2));
      pause(0.01)   %短暂暂停以保证图成功显示
    elseif NO==3
      clf
      scatter3(pnonf(:,1),pnonf(:,2),pnonf(:,3));
      pause(0.01)   %短暂暂停以保证图成功显示
    end
      
   
    %惯性权重,可以选择如下几种形式
    %w=0.5;
    w=w_start-(w_start-w_end)*(t/max_gen);
    %w=w_start-(w_start-w_end)*(t/max_gen)^2;
    %w=w_start-(w_start-w_end)*(2*t/max_gen-(t/max_gen)^2);
    %w=w_end*(w_start/w_end)^(1/(1+10*t/max_gen));
   
   
    %更新粒子种群的位置
   
    if limit_method==1
      %%限制条件策略:如不满足限制条件,重新产生移动,直至满足,为防止陷入局部死循环,设定不满足的次数越多,发生变异的概率越高,但有上限。
      for i=1:size_pop
            kc=1;
            while 1==1
                %更新速度
                V(i,:)=w*V(i,:)+c1*rand.*(ibestp(i,:)-p(i,:))+c2*rand.*(pbestp(i,:)-p(i,:));
                for j=1:n
                  if V(i,j)>V_max(j)
                        V(i,j)=V_max(j);
                  elseif V(i,j)<V_min(j)
                        V(i,j)=V_min(j);
                  end
                end
               
                %更新位置
                p(i,:)=p(i,:)+V(i,:);
                p(i,:)=get_discrete(p(i,:),p_discrete);
               
                %自适应变异,如果多次循环均未获得满足限制条件的粒子位置,增加变异率
                if rand<min(km*kc/10,0.1)
                  p(i,:)=rand(1,n).*(p_range(2,:)-p_range(1,:))+p_range(1,:);
                  p(i,:)=get_discrete(p(i,:),p_discrete);
                end

                %判断是否满足限制条件,如不满足则重新移动
                J_lim=estimate_limit(p(i,:),p_range);
                if J_lim==1
                  Obj=estimate_Obj(p(i,:));
                  f(i,:)=Obj;
                  break
                end
                kc=kc+1;
            end
            
            step1=t+i/1000%实时显示当前计算进度
      end
      
    elseif limit_method==2
      
      %限制条件策略:如不满足限制条件,目标函数取极大值,以便完全阻止其成为非劣解。
      for i=1:size_pop
            %更新速度
            V(i,:)=w*V(i,:)+c1*rand.*(ibestp(i,:)-p(i,:))+c2*rand.*(pbestp(i,:)-p(i,:));
            for j=1:n
                if V(i,j)>V_max(j)
                  V(i,j)=V_max(j);
                elseif V(i,j)<V_min(j)
                  V(i,j)=V_min(j);
                end
            end
            
            %更新位置
            p(i,:)=p(i,:)+V(i,:);
            p(i,:)=get_discrete(p(i,:),p_discrete);
            
            %参数上下限修正
            for j=1:n
               if p(i,j)>p_range(2,j)
                   p(i,j)=p_range(2,j);
               elseif p(i,j)<p_range(1,j)
                   p(i,j)=p_range(1,j);
               end
            end
            
            %自适应变异,如果多次循环均未获得满足限制条件的粒子位置,增加变异率
            if rand<km
                p(i,:)=rand(1,n).*(p_range(2,:)-p_range(1,:))+p_range(1,:);
                p(i,:)=get_discrete(p(i,:),p_discrete);
            end

            %判断是否满足限制条件,如满足,计算目标函数,否则目标函数取极大值
            J_lim=estimate_limit(p(i,:),p_range);
            if J_lim==1
                f(i,:)=estimate_Obj(p(i,:));
            else
                f(i,:)=ones(1,NO)*10^10;
            end
            step=t+i/1000
      end
    end
                  
    %更新个体极值
    for i=1:size_pop
      if sum(f(i,:)<ibestf(i,:))==NO
            ibestf(i,:)=f(i,:);
            ibestp(i,:)=p(i,:);
      elseif sum(f(i,:)<ibestf(i,:))>0&&sum(f(i,:)<ibestf(i,:))<NO
            if rand>=0.5
                ibestf(i,:)=f(i,:);
                ibestp(i,:)=p(i,:);
            end
      end
    end
   
    %更新群体非劣解集
    %首先求当前代的非劣解集
    =get_noninferior(p,f);
    %将历史非劣解集和当前代非劣解集合并,求出新的非劣解集
    pnonp_new_candidate=;
    pnonf_new_candidate=;
    =get_noninferior(pnonp_new_candidate,pnonf_new_candidate);
    %非劣解集网格划分
    grid_inform=grid_noninferior(pnonf,grid_num);
    %非劣解集删减
    =cut_noninferior(pnonp,pnonf,grid_inform,noninf_num);
    %重新网格划分
    grid_inform=grid_noninferior(pnonf,grid_num);
    %更新群体极值和群体适应度
    =select_pbest(pnonp,pnonf,grid_inform,size_pop);

    t=t+1;
end

end


function =get_noninferior(p,f)
    %由解集p及其目标值f,求出解集的非劣解集及目标值,并划分网格
    %判断解之间是否支配,第(i,j)项表示第i个解是否被第j个解支配,0表示不支配,1表示支配
    NO=length(f(1,:));
    np=length(p(:,1));
    dom=zeros(np,np);
    for i=1:(np-1)
      for j=(i+1):np
            ifdom=sum(f(i,:)<=f(j,:));
            if ifdom==0 %i被j支配
                dom(i,j)=1;
                dom(j,i)=0;
            elseif ifdom==NO %j被i支配
                dom(i,j)=0;
                dom(j,i)=1;
            else %互不支配
                dom(i,j)=0;
                dom(j,i)=0;
            end
      end
    end
    %非劣解集
    sdom=sum(dom,2);
    nond=find(sdom==0);
    pnonp=p(nond,:);
    pnonf=f(nond,:);
end

function grid_inform=grid_noninferior(pnonf,grid_num)
    %对非劣解集划分网格,返回网格密度信息
    %以结构体形式储存密度信息,包括有grid(所属网格),density(点的数量),non_num(属于该网格的点的编号,对应pnonp的行数)

    %非劣解集的size
    Nnon=length(pnonf(:,1));
    NO=length(pnonf(1,:));
    %网格的范围
    grid_limit=zeros(2,NO);
    grid_size=zeros(1,NO);
    for i=1:NO
      maxf=max(pnonf(:,i));
      minf=min(pnonf(:,i));
      grid_size(i)=(maxf-minf)/(grid_num-1);
      grid_limit(1,i)=minf-grid_size(i)/2;
      grid_limit(2,i)=maxf+grid_size(i)/2;
    end
    %所属网格
    dgrid=[]; %用于保存密度大于1的网格的号码
    grid_inform=struct;%存储网格密度和属于该网格的非劣解的序号
    for i=1:Nnon
      gridi=zeros(1,NO);
      for j=1:NO
            gridi(j)=ceil((pnonf(i,j)-grid_limit(1,j))/grid_size(j));
      end
      s=size(dgrid);
      k=1;
      while k<=s(1)
            if sum(gridi==dgrid(k,:))==NO
                grid_inform(k).grid=gridi;
                grid_inform(k).density=grid_inform(k).density+1;
                grid_inform(k).non_num=;
                break
            end
            k=k+1;
      end
      if k>s(1)
            dgrid=;
            grid_inform(k).grid=gridi;
            grid_inform(s(1)+1).density=1;
            grid_inform(s(1)+1).non_num=i;
      end
    end
end

function =cut_noninferior(pnonp,pnonf,grid_inform,noninf_num)
%对非劣解集进行删减
Nnon=length(pnonf(:,1));
if Nnon>noninf_num
    Ngrid=length(grid_inform);
    Nun_non=[];
    for i=1:Ngrid
      delN=round((Nnon-noninf_num)/Nnon*grid_inform(i).density);
      if delN<grid_inform(i).density*0.8
            nonnum=grid_inform(i).non_num;
            for j=1:delN
                dN=ceil(rand*(grid_inform(i).density+1-j));
                nonnum(dN)=[];
            end
            grid_inform(i).non_num=nonnum;
      end
      Nun_non=;
    end
    pnonp=pnonp(Nun_non,:);
    pnonf=pnonf(Nun_non,:);
end
end

function p=get_discrete(p,p_discrete)
    %按照离散度定义,更新参数
    n=length(p);
    for j=1:n
      if p_discrete(j)==0
            p(j)=p(j);
      elseif p_discrete(j)==1
            p(j)=round(p(j));%取整
      else
            dis=mod(p(j),p_discrete(j));
            if dis/p_discrete(j)<0.5
                p(j)=p(j)-dis;
            else
                p(j)=p(j)+p_discrete(j)-dis;
            end
      end
    end
end

function =select_pbest(pnonp,pnonf,grid_inform,size_pop)
%采用网格法选择群体极值
Ngrid=length(grid_inform);
dens=zeros(1,Ngrid);
pp=zeros(1,Ngrid+1);
for i=1:Ngrid
    dens(i)=grid_inform(i).density;
end
for i=1:Ngrid
    pp(i+1)=sum(1./dens(1:i))/sum(1./dens);
end

n=length(pnonp(1,:));
NO=length(pnonf(1,:));
pbestp=zeros(size_pop,n);
pbestf=zeros(size_pop,NO);
for i=1:size_pop
    kk=rand;
    for j=1:Ngrid
      if kk>=pp(j)&&kk<=pp(j+1)
            idgrid=j;
            break
      end
    end
    idp=ceil(rand*grid_inform(idgrid).density);
    gridp=grid_inform(idgrid).non_num;
    sele_n=gridp(idp);
    pbestp(i,:)=pnonp(sele_n,:);
    pbestf(i,:)=pnonf(sele_n,:);
end
end

HuldaGnodim 发表于 2022-8-19 14:35

博主好高产[大笑]

stonstad 发表于 2022-8-19 14:36

图文并茂[赞同][赞同][赞同][谢邀][谢邀][谢邀]
页: [1]
查看完整版本: 非线性优化算法(二)多目标粒子群算法-matlab实现