MOPSO 多目标粒子群（Matlab）

MOPSO

本文提出了一种将Pareto占优性融入粒子群优化( PSO )的方法，以允许这种启发式方法处理多个目标函数的问题。与目前其他扩展PSO来解决多目标优化问题的方案不同，我们的算法使用了一个次级的(即,外部的)粒子库，这些粒子后来被其他粒子用来指导自己的飞行。我们还加入了一个特殊的变异算子，丰富了算法的探索能力。使用几个测试函数和从标准的进化多目标优化文献中提取的度量来验证所提出的方法。结果表明该方法具有很强的竞争力，可以被认为是解决多目标优化问题的一种可行的选择。

• 主要创新

我们的方法与文献中存在的其他建议的主要区别是：

• 我们采用类似于 PAES [21] 自适应网格的外部（或辅助）存储库。

• 我们使用的变异算子既作用于群体的粒子，也作用于待解决问题的每个设计变量的范围（参见第 IV-C 节）。这不仅旨在探索搜索空间的偏远区域，而且还试图确保探索每个决策变量的全部范围。

• 我们对MOPSO 参数对其性能的影响进行了广泛的分析。我们还使用三个指标将我们的 MOPSO 与其他三种算法（代表进化多目标优化的最新算法）进行比较。

• 算法框架

1. 初始化种群POP

   FOR i=0 TO MAX /* MAX 种群大小*/

   初始化 POP[i]

2. 初始化每个粒子的速度：

   FOR　ｉ＝０　ＴＯ　ＭＡＸ

   VEL[i]＝0

3. 评估种群中每一个粒子

4. 存储存储库中表示非支配向量的粒子的位置

5. 生成迄今为止探索的搜索空间的超立方体，并使用这些超立方体作为一个坐标系来定位粒子，其中每个粒子的坐标根据其目标函数的值定义。

6. 初始化每个粒子的记忆(这个记忆充当一个穿过搜索空间的指南。该内存也存储在存储库中)：

   FOR i=0 TO MAX

   PBEST[i]=POP[i]

7. WHILE 最大循环次数尚未达到 DO

         a)使用以下表达式计算每个粒子1的速度：
   

   $$\begin{array}{rl}VEL[i]& =W\times VEL[i]+R_1\\ & \times (PBESTS[i]-POP[i])+R_2\times (REP[h]-POP[i])\end{array}$$

            其中(惯性权重)取0.4；且均为范围内的随机数；是粒子已经拥有的最佳位置；2是从存储库中取出的值；索引的选择方式如下：那些包含一个以上粒子的超立方体被赋予一个适应值，该适应值等于任意数量(我们在实验中使用)除以其包含的粒子数。这样做的目的是减少那些包含更多粒子的超立方体的适应度，可以看作是适应度共享的一种形式。然后，我们应用轮盘赌选择使用这些适应度值来选择超立方体，从中我们将采取相应的粒子。一旦超立方体被选中，我们在超立方体中随机选择一个粒子，就是该粒子的当前值。
   
            
   
            b)计算粒子的新位置加上上一步产生的速度
   

   $$POP[i]=POP[i]+VEL[i]$$

               c)保持搜索空间内的粒子，以防超出其边界(避免产生不依赖于有效搜索空间的解)。当决策变量超出其边界时，我们做两件事：
   
                      1 )决策变量取其对应边界(无论是下边界还是上边界)的值；
   
                      2 )决策变量的速度乘以( -1 )，使其向相反的方向搜索。
   
              d)对超立方体中的每个粒子进行评价。
   
              e)更新超立方体中粒子的内容和地理表示。该更新包括将当前所有非支配位置插入到存储库中。在这个过程中，从存储库中的任何占主导地位的位置都被消除。由于存储库的大小是有限的，每当存储库满时，我们应用第二个保留准则：那些位于目标空间中人口较少的区域的粒子优先于那些位于人口密集区域的粒子。                        
   
             f)当粒子的当前位置优于其内存中包含的位置时，粒子的位置使用更新
   

$$PBESTS[i]=POP[i]$$

                  决定从记忆中保留什么位置的标准是简单地应用帕累托占优(即,如果当前位置由内存中的位置控制,则保留内存中的位置;否则,当前位置替换内存中的位置;如果两个人都不被另一个人支配,那么我们随机选择其中一个)。

            g)增加循环计数器

  8. END WHILE

• 主要不同

• 外部存储库

  外部存储库(或归档文件)的主要目标是保存搜索过程中发现的非支配向量的历史记录。外部存储库由两个主要部分组成：存档控制器和网格。

  1 )存档控制器：存档控制器的功能是决定是否要在存档中添加某个解决方案。

  在我们算法的主种群中每次迭代找到的非支配向量相对于外部存储库的内容被比较(在一对一的基础上)，在搜索开始时将是空的。
  
  如果外部存档为空，则接受当前解(见图情形1)。如果这个新的解是由外部存档中的一个个体支配，那么这样的解会被自动丢弃(见图情形2 )。否则，如果外部种群所包含的解中没有一个支配将要进入的解，那么这样的解就存储在外部存档中。如果归档文件中存在由新解支配的解，则将此类解从归档文件(见图1中的情形3和情形4)中移除。最后，如果外部种群已经达到其最大允许容量，则调用自适应网格程序(见图情形5)。

  2 )网格：为了产生均匀分布的Pareto前沿，我们的方法使用了自适应网格的变化 。

  其基本思想是使用一个外部存档来存储所有相对于存档内容非支配的解。 自适应网格实际上是由超立方体组成的空间。这种超立方体具有和目标函数一样多的组件。每个超立方体可以被解释为一个包含无数个体的地理区域。自适应网格的主要优点是计算成本低于小生境。然而，网格必须在每一代更新。在这种情况下，自适应网格的计算复杂度将与小生境相同，即O(N²) .

• 变异操作

众所周知，PSO具有非常高的收敛速度。然而，在多目标优化的背景下，这种收敛速度可能是有害的，因为基于PSO的算法可能会收敛到一个虚假的Pareto前沿(即相当于全局优化中的一个局部最优)。PSO的这一缺点在一些问题中很明显，在这些问题中，我们原来的方法并没有很好地表现出来。这激发了变异算子的发展，它试图在搜索开始时与所有的粒子一起探索。然后，快速减少受变异算子(见图4)影响的粒子数(关于迭代次数)。值得注意的是，我们的变异算子不仅适用于群体的粒子，也适用于待求解问题(使用相同的变异函数)的每个设计变量的范围。这样做的目的是在搜索开始时覆盖每个设计变量的全部范围，然后使用非线性函数缩小覆盖的范围。从图4可以看出，在开始时，种群中的所有粒子都受到变异算子(以及决策变量的全范围)的影响。这就意图在算法中产生一种高度的探索行为。随着迭代次数的增加，变异算子的效果下降。

• MOPSO代码：

for it = 1:MaxIt
    
    for i = 1:nPop
        
        leader = SelectLeader(rep, beta);
        
        pop(i).Velocity = w*pop(i).Velocity ...
            +c1*rand(VarSize).*(pop(i).Best.Position-pop(i).Position) ...
            +c2*rand(VarSize).*(leader.Position-pop(i).Position);
        
        pop(i).Position = pop(i).Position + pop(i).Velocity;
        
        pop(i).Position = max(pop(i).Position, VarMin);
        pop(i).Position = min(pop(i).Position, VarMax);
        
        pop(i).Cost = CostFunction(pop(i).Position')';
        
        % Apply Mutation
        pm = (1-(it-1)/(MaxIt-1))^(1/mu);
        if rand<pm
            NewSol.Position = Mutate(pop(i).Position, pm, VarMin, VarMax);
            NewSol.Cost = CostFunction(NewSol.Position')';
            if Dominates(NewSol, pop(i))
                pop(i).Position = NewSol.Position;
                pop(i).Cost = NewSol.Cost;

            elseif Dominates(pop(i), NewSol)
                % Do Nothing

            else
                if rand<0.5
                    pop(i).Position = NewSol.Position;
                    pop(i).Cost = NewSol.Cost;
                end
            end
        end
        
        if Dominates(pop(i), pop(i).Best)
            pop(i).Best.Position = pop(i).Position;
            pop(i).Best.Cost = pop(i).Cost;
            
        elseif Dominates(pop(i).Best, pop(i))
            % Do Nothing
            
        else
            if rand<0.5
                pop(i).Best.Position = pop(i).Position;
                pop(i).Best.Cost = pop(i).Cost;
            end
        end
        
    end
    
    % Add Non-Dominated Particles to REPOSITORY
    rep = [rep
         pop(~[pop.IsDominated])]; %#ok
    
    % Determine Domination of New Resository Members
    rep = DetermineDomination(rep);
    
    % Keep only Non-Dminated Memebrs in the Repository
    rep = rep(~[rep.IsDominated]);
    
    % Update Grid
    Archive_costs=GetCosts(rep);
    Grid = CreateGrid(Archive_costs, nGrid, alpha);

    % Update Grid Indices
    for i = 1:numel(rep)
        rep(i) = FindGridIndex(rep(i), Grid);
    end
    
    % Check if Repository is Full
    if numel(rep)>nRep
        
        Extra = numel(rep)-nRep;
        for e = 1:Extra
            rep = DeleteOneRepMemebr(rep, gamma);
        end
        
    end
    
    % Plot Costs
    figure(1);
    PlotCosts(pop, rep);
    pause(0.01);
    
    % Show Iteration Information
    disp(['Iteration ' num2str(it) ': Number of Rep Members = ' num2str(numel(rep))]);
    
    % Damping Inertia Weight
    w = w*wdamp;
    
end

扫码关注公众号获取更多代码