多目标正余弦算法多目标优化算法 python代码多目标规划求解算法

转载

mob64ca1407d5aa 2024-08-09 15:55:30

文章标签 matlab 算法开发语言 ci 迭代 文章分类 Python 后端开发

1 算法介绍

1.1 关于速度和位置

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。

鸟被抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子i在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)。每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)，这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。

2.2 速度和位置的更新

PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。

多目标正余弦算法多目标优化算法 python代码多目标规划求解算法_matlab

对于公式(1)：

公式(1)的第①部分称为【记忆项】，表示上次速度大小和方向的影响；

公式(1)的第②部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；

公式(1)的第③部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。

以上面两个公式为基础，再来看一个公式：

多目标正余弦算法多目标优化算法 python代码多目标规划求解算法_开发语言_02

公式(2)和公式(3)被视为标准PSO算法。

1.3 标准PSO算法的流程

1）初始化一群微粒(群体规模为N)，包括随机位置和速度；

2）评价每个微粒的适应度；

3）对每个微粒，将其适应值与其经过的最好位置pbest作比较，如果较好，则将其作为当前的最好位置pbest；

4）对每个微粒，将其适应值与其经过的最好位置gbest作比较，如果较好，则将其作为当前的最好位置gbest；

5）根据公式(2)、(3)调整微粒速度和位置；

6）未达到结束条件则转第2）步。

迭代终止条件根据具体问题一般选为最大迭代次数Gk或(和)微粒群迄今为止搜索到的最优位置满足预定最小适应阈值。

## 4.2 PSO流程图解

多目标正余弦算法多目标优化算法 python代码多目标规划求解算法_算法_03

2 部分代码

clc;
clear;
close all;

%% 问题定义，这里可以删除

TestProblem=1;   % Set to 1, 2, or 3

switch TestProblem
   case 1
       CostFunction=@(x) MyCost1(x);
       nVar=50;
       VarMin=0;
       VarMax=1;
       
   case 2
       CostFunction=@(x) MyCost2(x);
       nVar=3;
       VarMin=-5;
       VarMax=5;
       
   case 3
       CostFunction=@(x) MyCost3(x);
       nVar=2;
       VarMin=0;
       VarMax=1;
end

VarSize=[1 nVar];

VelMax=(VarMax-VarMin)/10;

%% MOPSO 设置
nPop=500;   % Population Size

nRep=100;   % Repository Size

MaxIt=500;  % Maximum Number of Iterations

phi1=2.05;%学习因子
phi2=2.05;
phi=phi1+phi2;
chi=2/(phi-2+sqrt(phi^2-4*phi));

w=chi;              % Inertia Weight
wdamp=1;            % Inertia Weight Damping Ratio
c1=chi*phi1;        % Personal Learning Coefficient
c2=chi*phi2;        % Global Learning Coefficient

alpha=0.1;  % Grid Inflation Parameter

nGrid=10;   % Number of Grids per each Dimension

beta=4;     % Leader Selection Pressure Parameter

gamma=2;    % Extra (to be deleted) Repository Member Selection Pressure

%% 初始化

particle=CreateEmptyParticle(nPop);

for i=1:nPop
   particle(i).Velocity=0;
   particle(i).Position=unifrnd(VarMin,VarMax,VarSize); %unifrnd在[VarMin,VarMax]区间模拟VarSize的数组
   % 初始化了自变量X的位置
   particle(i).Cost=CostFunction(particle(i).Position); % 求解目标函数
   particle(i).Best.Position=particle(i).Position;
   particle(i).Best.Cost=particle(i).Cost;
end

particle=DetermineDomination(particle);

rep=GetNonDominatedParticles(particle);

rep_costs=GetCosts(rep);
G=CreateHypercubes(rep_costs,nGrid,alpha);

for i=1:numel(rep)
  [rep(i).GridIndex rep(i).GridSubIndex]=GetGridIndex(rep(i),G);
end
   
%% MOPSO 主循环

for it=1:MaxIt
   for i=1:nPop
       rep_h=SelectLeader(rep,beta);

       particle(i).Velocity=w*particle(i).Velocity ...
                            +c1*rand*(particle(i).Best.Position - particle(i).Position) ...
                            +c2*rand*(rep_h.Position -  particle(i).Position);

       particle(i).Velocity=min(max(particle(i).Velocity,-VelMax),+VelMax);

       particle(i).Position=particle(i).Position + particle(i).Velocity;
   
       flag=(particle(i).Position<VarMin | particle(i).Position>VarMax);    %非劣检查:flag
       particle(i).Velocity(flag)=-particle(i).Velocity(flag);
       
       particle(i).Position=min(max(particle(i).Position,VarMin),VarMax);

       particle(i).Cost=CostFunction(particle(i).Position);

       if Dominates(particle(i),particle(i).Best)
           particle(i).Best.Position=particle(i).Position;
           particle(i).Best.Cost=particle(i).Cost;
           
       elseif ~Dominates(particle(i).Best,particle(i))
           if rand<0.5
               particle(i).Best.Position=particle(i).Position; 
               particle(i).Best.Cost=particle(i).Cost;
           end
       end

   end
   
   particle=DetermineDomination(particle);
   nd_particle=GetNonDominatedParticles(particle);
   
   rep=[rep
        nd_particle];
   
   rep=DetermineDomination(rep);
   rep=GetNonDominatedParticles(rep);
   
   for i=1:numel(rep)
      [rep(i).GridIndex rep(i).GridSubIndex]=GetGridIndex(rep(i),G);
   end
   
   if numel(rep)>nRep
       EXTRA=numel(rep)-nRep;
       rep=DeleteFromRep(rep,EXTRA,gamma);
       
       rep_costs=GetCosts(rep);
       G=CreateHypercubes(rep_costs,nGrid,alpha);
       
   end
  
%     disp(['Iteration ' num2str(it) ': Number of Repository Particles = ' num2str(numel(rep))]);
   
   w=w*wdamp;
end

%% 结果

costs=GetCosts(particle);
rep_costs=GetCosts(rep);

figure;

plot(costs(1,:),costs(2,:),'b.');
hold on;
plot(rep_costs(1,:),rep_costs(2,:),'rx');
legend('Main Population','Repository');