基于大规模连续多目标优化的共轭梯度-进化集成算法

声明：文章题目字数有限，翻译水平有限，仅供参考！

原题目: Integrating Conjugate Gradients Into Evolutionary Algorithms for Large-Scale Continuous Multi-Objective Optimization

引：这么久了，又捡起来写博客，是最近本人突然明白了很多。每每遇到学习、科研方面的重大阻碍时，记录往往可以让自己的内心更加踏实。这种踏实感在告诉你，不要计较所谓的成功与失败，只要不断进步就好。心灵的平复可以让自己头脑更加清醒，所谓的挫折、失败都是需要直视的，解决也好、解决不了也好，努力过、充实过就可以了，是的，每一次进步就是成功！

大过年的，画一幅画给大家乐呵乐呵~~🤗

进步从学习一篇文章开始！

参考资料

论文：
1、结合进化计算和数学规划求解大规模多目标优化问题
2、Integrating Conjugate Gradients Into EvolutionaryAlgorithms for Large-Scale ContinuousMulti-Objective Optimization

博客：
1、共轭梯度法的简单直观理解
2、一文了解差分进化算法的前世今生
3、从为什么梯度方向是函数变化率最快方向详谈梯度下降算法
4、梯度下降详解（主观理解+推导证明+例题）
5、共轭梯度法的简单直观理解
6、差分进化算法原理及应用
Code:
1、Differential Evolution (DE)
2、PlatEMO

在讲解文章之前，我们先来讲一下文章所需的一些基础知识，比较熟悉共轭梯度和差分进化的可自行省略一、二部分。*

1、数学优化

1.1、梯度下降

在说明共轭梯度之前，我们先来讲一下什么是梯度下降（让与我一样数学知识空白的童鞋们充一下电）。为了简便起见，我直接附上了我的笔记，大家可自行学习。

这里我们再扩充一点：为什么梯度下降的方法可以寻找到最小值？（最大值问题可以转换成最小值问题，我们统称为最小值问题）

很多人将梯度下降问题形象的描述为下山问题。如下图：

那么如何从山峰走到山谷呢？其实，证明思想比较简单，我们不妨简单的讲一下。
我们用$f(x)$表示我们的目标函数，假设该函数是连续可微函数，下山的初始点为$x_0$，我们用一阶泰勒展开式表示下山的下一个点：
$$f(x)\approx f(x_0)+f(x_0{)}^{{}^{\prime }}(x-x_0)$$
我们用$dv$表示两个点距离的差值，则上述式子可以简化为：
$$f(x)\approx f(x_0)+f(x_0{)}^{{}^{\prime }}dv$$
我们知道$f(x_0){}^{\prime }$是函数在$x_0$点处的导数，而导数可以用梯度的形式表示。要想$f(x)$变小，那么$dv$的方向应该与梯度的方向相反才行，即梯度的反方向，也就是说，要想获得$f(x)$最优值，我们应该沿着 梯度下降的方向寻找可行的解。

1.2、共轭梯度

有关为什么叫共轭梯度，为什么共轭梯度可以求解最优化问题，本文将不再赘述，请参照博文：共轭梯度法的简单直观理解。这里我们仅讲解一下共轭梯度如何求解。
步骤（我这里就偷个懒复制粘贴了🤗）：

如果这种方式比较晦涩的话，我们不妨直接来看F-R共轭梯度的求解过程，如下：

1.3、梯度实现

对于单值多元函数，有些函数难以获得解析导数，因此编程实现了导数的数值算法，基本思想就是用一阶差分代替导数：
$$f(x_0{)}^{{}^{\prime }}\approx li{m}_{dx\to 0}\frac{f(x_0+dx)-f(x_0)}{dx}$$
式中，$dx$可以用$\Delta x$来表示，表示在$x$方向的微小增量；$x_0$表示初值，一般可以设$dx=10^{-6}$。
当$X$为多维变量、目标有多个时，上式可转换成：
$$f(x_0{)}^{{}^{\prime }}\approx li{m}_{\epsilon \to 0}\frac{f_i(X+\epsilon )-f_i(X)}{\epsilon },i=1,2,...,M$$
式中，$\epsilon$表示在$X$方向的微小增量，这样貌似已经完成了梯度的近似求解。即，将微分方程转换为代数方程式的形式，用有限差分来近似求解。
但，上述方程仍存在以下问题：

“需要注意的是，这种近似是不太可靠的，因为一阶偏导数$\frac{\partial f_i(x)}{\partial x_j}$ 还没有得到，而且雅可比矩阵仍然是未知的。如果计算一个标量$f(x{)}^{{}^{\prime }}$而不是一个向量$(\frac{\partial f_i(x_1)}{\partial x_1},\frac{\partial f_1(x)}{\partial x_2},...,\frac{\partial f_1(x)}{\partial x_D})$ ，所有的决策变量只能以相同的方向和相同的步长进行更新，这在达到全局最优方面显然是不太适用。”

如何理解上边这段话呢？这里我附上我的理解（理解不一定合理，有更具备数学依据的请附在评论中），我们不妨将其放入到二元方程中，如下图：

在二元方程中，两个变量$x_1、x_2$从初始（黑线）到加入$\epsilon$（红线），这两个变量的方向是一致的，即变量的方向是一致的。三元的可以看成一个面，如下图。

那么如何体现全局最优呢？

文中，给出了一种近似雅可比矩阵的方式：
$$f(x_0{)}^{{}^{\prime }}\approx li{m}_{\epsilon \to 0}\frac{f_i(x_{j+\epsilon })-f_i(X)}{\epsilon },i=1,2,...,M；j=1,2,...,D$$
其中，$x_{j+\epsilon }=(x_1,x_2,...,x_{j-1},x_j+\epsilon ,x_{j+1},...,x_D)$。这样，目标函数的一阶偏导数就可以通过将$\epsilon$设为一个很小的值来正确地近似。

需要强调的是在对j维求偏导时，仅仅是变化的j维参量，其他维度是不变的。这一点，可以在代码中体现。详见3.2.3计算梯度

拓展: 关于微分、导数、差分之间的关系，我比较赞同博主写的这篇文章，这里我仅粘贴文章中一幅图，具体请参考：微分、导数与差分

2、进化计算—差分进化

之前一直觉得差分进化是个很神秘的算法，会比较复杂，其实当真正想要了解它时，才发现差分进化只不过是遗传算法的一种 “变种”。
为什么说是变种呢？因为差分进化的过程也包含：初始化、变异、交叉和选择四个主要过程。两者不同之处在于变异与交叉的实现过程。有关遗传算法变异与交叉过程的具体实现过程，我不再赘述，不了解的请参考我转载的另一篇博客：遗传算法的交叉变异详解。我们这里仅说明一下差分进化变异与交叉过程。

2.1、变异

假设种群个数为$N$，决策变量为$D$，为了方便，我们设定优化问题为单目标优化问题，则种群个体可以表示为：
$$\overrightarrow{x_i}=（x_i^1,x_i^2,...,x_i^D）,i=1,2,...,N$$
定义$x_{i,G}$为第$G$代种群的第$i$个个体。$x_{i,G}$ 通过以下的变异算子生成与之对应的变异向量 $v_{i,G}$:
​$$v_{i,G}=x_{r_1^i,G}+F\cdot (x_{r_2^i,G}-x_{r_3^i,G})$$
其中，$r_1^i，r_2^i，r_3^i$是范围$1$~$N$之间的随机数,且三者不等于$i$; $x_{i,G}$ 称为$v_{i,G}$的目标向量；$x_{r_1^i,G}$为基向量；$(x_{r_2^i,G}-x_{r_3^i,G})$为差分向量，$F$为缩放因子。

实现：

        x=pop(i).Position;
        
        A=randperm(nPop);
        
        A(A==i)=[];
        
        a=A(1);
        b=A(2);
        c=A(3);
        
        % Mutation
        %beta=unifrnd(beta_min,beta_max);
        beta=unifrnd(beta_min,beta_max,VarSize);
        y=pop(a).Position+beta.*(pop(b).Position-pop(c).Position);
        y = max(y, VarMin);
		y = min(y, VarMax);

2.2、交叉

交叉算子如下所示，其作用是通过组合$x_{i,G}$ 和$v_{i,G}$ 产生对应的试验向量$u_{i,G}$ :
$$u_{i,G}=\{\begin{array}{l}{x}_{i,G}\\ v_{i,G},if(j=j_{rand}or{r}_i^j<Cr)\end{array}$$
其中$Cr$为交叉概率。

实现：

	    z=zeros(size(x));
        j0=randi([1 numel(x)]);
        for j=1:numel(x)
            if j==j0 || rand<=pCR
                z(j)=y(j);
            else
                z(j)=x(j);
            end
        end

2.3 几何上的差分进化

具体的可参考：差分进化算法原理及应用，为了更加便于学习，我们偷个懒，对该博文进行了简单的复制，侵权则删。
首先我们附上一差分进化过程完整的几何变化过程，截自：一文了解差分进化算法的前世今生

1.产生初始种群的空间几何分布如下：

可见此时的种群分布是杂乱无章的，完全离最优区域的领域很远。因此，需要对种群进行变异，以增加子代种群的多样性。

2. 变异过程的几何表示：
   

可见通过变异的操作，使着部分种群经过差分改变其运动方向。

3. 交叉操作的几何表示：
   
   交叉操作，通过随机比较，产生子代U，以方便通过贪婪选择进行比较选择最优子代，达到最优化的目的。从上图 可以看出子代的选择方式的几何分布特征，最终朝着最优区域前进。

4.选择操作

可以看出最终的子代都停留在3.85这个位置，说明此时的种群已经收敛到最优值上了。

3、文章实现

3.1、Why？

首先，我们先说明一下为什么要采用数学优化与进化计算结合的方式，其原因我直接粘贴原文：
简单来概括：数学优化方法难以求解大规模多目标优化问题，而进化计算方法因其计算量大的问题，其求解效率仍存在缺陷。
(大规模多目标问题：决策变量上百个+目标函数 $\ge$ 2)
那么当前求解大规模多目标问题，现在研究如何做的呢？文中列出了三种主流解决方法（数学优化方法存在解多样性缺乏的问题，仅考虑进化算法）：

🐉 基于不同策略的决策变量分组方法

“主要包括随机分组、差分分组和变量分析。随机分组是非常有效的，因为它将决策变量随机划分为具有相同大小的预定义数量的组，但由于变量之间的相互作用被完全忽略，它可能会驱动解走向局部最优。差分分组可以检测到变量之间的相互作用，对于逼近全局最优解更有效，但由于每两个变量之间的相互作用应该独立检测，因此它具有相当高的计算复杂度。变量分作为差分分组的补充，它可以区分促进收敛性和多样性的变量，但仍然很耗时。”

🐉 问题转换或降维的方法

“问题变换策略不是优化决策变量，而是通过优化作用于现有解决方案的权向量来探索巨大的搜索空间，由于权向量比决策向量短得多，可以大大减少优化变量的数量。但这是以牺牲有效性为代价的，局部最优解可以被有效地找到，而全局最优解很可能在权值向量空间中丢失。相比之下，降维策略旨在通过机器学习、数据挖掘或其他技术来减少搜索空间，这些技术考虑变量交互，可以在减少的搜索空间中保留全局最优解。然而，降维策略只适用于特定的LSMOPs，如那些具有低内在维数或稀疏最优解的LSMOPs。”
🐉 新的变异算子或概率模型

“以有效地近似于原始搜索空间中的全局最优解。自定义变异算子是对现有的遗传算子、竞争群优化器和协方差矩阵自适应进化策的增强。定制的概率模型通过描述迄今为止发现的有希望的解来生成新的解决方案，这是由高斯过程、生成对抗网络、基于镜像划分的求解知识和其他技术实现的。这些搜索策略可以提高现有MOEAs在LSMOPs上的收敛速度，但仍需要大量的函数评估才能达到全局帕累托前沿。”

3.2、How？

那如何解决进化算法与数学优化所存在的问题呢？一种少量研究采取的进化计算与数学优化方法结合的方式有望成为解决这一困境的有力工具。当前，基于进化计算与梯度下降的混合方法还未得到充分的关注，少量学者已进行了这种混合算法的探究，但 “现有工作通过基于梯度的局部搜索来简单地跟踪变异算子，这加快了收敛速度，但不利于目标空间中的多样性保存。其次，现有的工作集中于优化单一目标，而很少涉及多个目标的梯度的聚合。总的来说，这对现有的混合算法提出了挑战。”

本文为了解决上述问题，提出了进化计+差分进化结合的方式，构建了一种求解大规模多目标问题的混合算法—MOCGDE，该算法流程图如下：

3.2.1 算法总体实现

为了更加详细说明算法步骤，我们在这里直接粘贴上算法的伪代码：

为了更清楚地理解算法，我这里对算法的伪代码进行解释，如下：

🐉🐉🐉🐉🐉🐉🐉🐉🐉🐉 算法1： 算法流程🐉🐉🐉🐉🐉🐉🐉 🐉🐉🐉
🐉行号	🐉解释
Input	决策变量、种群个数；Archive大小，用于存储最优解，也就是进化中所应用的Pareto解中Rank=1的种群个体
Output	最优解集：Archive
2	随机生成均匀分布的点，这里主要用于生成均匀分布的向量，可以采用normal-boundaryintersection method等方法
3	更新存档的详细过程在算法2，保证存档不超过其最大大小。
4~6	K:计数器；G：梯度集，用于保存共轭梯度求解的梯度值；S：搜索方向集，即梯度下降方向，也就是梯度反方向。
9~16	在每次迭代中，种群中的每个解依次通过共轭梯度(算法3)和差分进化进行更新，并使用新生成的解来更新种群和存档。（共轭梯度与差分进化的知识详看：1、数学优化和2、进化计算—差分进化）
	解在梯度的引导下可以螺旋收敛到最优，但需要搜索方向周期性地恢复到梯度的负值，以确保共轭。这是通过一个计数器k来实现的，如果k是决策变量数量D的倍数，则S被设置为-g，并且k在每次迭代结束时增加1
17	利用得到的搜索方向(梯度下降方向)和差分进化更新解，实则是共轭梯度下降方法与差分进化交替使用，为了确定一个合适的步长，文章使用了一个简单的线搜索策略，步长随着m从1开始指数级地减小，直到找到一个更好的解决方案，详见算法4
18~23	若成功，则子代替换父代，更新Archive;否则，从Archive中随机选择一个解替换当前父代个体，详见算法4（需要说明的是，在算法4中还采用了约束违背的思想处理含约束的大规模多目标优化问题，约束违背的思想详见我的另一篇博文: 基于参考点的非支配遗传算法-NSGA-III（二）

3.2.2 更新Archive

更新Archive的主要目的就是：保存不大于最大Archive大小的种群，且该种群的Rank=1。此部分伪代码如下：

此伪代码部分有些难理解，其实就是寻找两两最小距离，每个种群个体x均能找到一个种群个体y, <$\vec{x}$,$\vec{y}$>距离最小，那么到底是x还是y呢？我们让两个都进行对比一下，也就是让除了{x,y}之外的所有个体种群分别对x和y求距离，即 $<x,\vec{z}>$,$<y,\vec{z}>$均是向量。那么代码实现其实就比较好容易了，我只要找到所有欧式距离最小的那个种群对即可，即$<x,z^{\ast }>$或者$<y,z^{\ast }>$，可通过求每一行中最小，然后再求所有行中最小中的最小值。

🐉🐉🐉🐉🐉🐉🐉🐉 Algorithm 2：Update Archive 🐉🐉🐉🐉🐉🐉🐉 🐉
🐉行号	🐉解释
Input	这里的A= [上一迭代Archive，子代OffSpring]
Output	更新后的Archive
1	移除种群A中非支配解（这些非支配解肯定来自于OffSpring）
2~9	逐个删除与其他解之间欧氏距离最短的解，直到存档不超过其最大大小

2~9具体实现我们通过来代码能更容易理解，这里我附上代码，代码中我已附上注释。

function P = UpdateArchive(P,N)
    P = P(NDSort(P.objs,P.cons,1)==1);    %算法第1行
    if length(P) > N
        Choose = true(1,length(P));  %判定是否被选择，已选择了给定false
        Dis    = pdist2(P.objs,P.objs);  %计算两两欧式距离
        Dis(logical(eye(length(Dis)))) = inf;  %自身距离设定inf
        while sum(Choose) > N
            Remain   = find(Choose);   %还未选择的保存在Remain
            Temp     = sort(Dis(Remain,Remain),2); %每行从小到大排序
            [~,Rank] = sortrows(Temp);%按照行大小排序，每行从第一列开始对比，依次进行，Rank记录行索引
            Choose(Remain(Rank(1))) = false; %最小的行设定被选择，即将要被删除
        end
        P = P(Choose); %删除距离较小的种群个体，保证存档大小不超过最大值
    end
end

3.2.3 计算共轭梯度

关于共轭梯度计算方法，我们在前边中已经介绍，这里仍然先粘贴上此部分的伪代码：

🐉🐉🐉🐉🐉🐉🐉🐉 Algorithm 3：Calculate Gradient 🐉🐉🐉🐉🐉🐉🐉 🐉
🐉行号	🐉解释
Input	种群 x, 权重 w (需要注意的是文章对于约束也进行了处理，用于聚合约束和目标的权重向量是不同的。由于必须满足所有的约束条件，所以用于聚合约束的权重向量被设置为一个量为1的向量；由于解应该在目标之间做出不同的权衡，因此用于聚合目标的权向量是预定义的（也就是我们前边采用NBI生成的均分分布的权重值），其中每个解都与一个不同的权向量以及收敛方向相关联。)
Output	梯度g（含约束的问题返回的是约束的梯度，否则，返回目标函数的梯度）和diff（标志符：用于表示种群个体在各目标函数梯度符号是否相同）
1	h(x)>0表示种群不满足目标约束，即所得解为不可行解，针对此种情况执行2~5，仅需处理约束梯度即可，因为目标梯度计算此时是没有意义的；否则的话执行 7~9，即所得解为可行解，仅需求解目标梯度即可。
2~5	获取约束梯度值，构建雅可比矩阵；同时乘上权重1，获得含权重的梯度。
7~9	获取目标函数梯度值，构建雅可比矩阵；同时乘上目标权重，获得含权重的梯度。

前边我们提到“需要强调的是在对j维求偏导时，仅仅是变化的j维参量，其他维度是不变的。”代码中目标与约束求导过程则分别用代码表示为：

%%目标函数求导
 function ObjGrad = CalObjGrad(obj,Dec)
	%CalObjGrad: Calculate the gradients of the objectives of a solution.
    % Grad = obj.CalObjGrad(Dec) returns the gradients of objectives
	% of Dec, i.e., a Jacobian matrix.
    % This function is usually called by gradient-based algorithms.
	%   Example:
    %       ObjGrad = Problem.CalObjGrad(Dec)

            Dec(Dec==0) = 1e-12;
            X           = repmat(Dec,length(Dec),1).*(1+eye(length(Dec))*1e-6);
            ObjGrad     = (obj.CalObj(X)-repmat(obj.CalObj(Dec),size(X,1),1))'./Dec./1e-6;
 end

即将$X$设为$D\ast D$的矩阵
1、$X$：$f(x_j+dx)$，第 j 行变化，仅第 j 列变化，其他不变；
2、$f(x_j+dx)-f(x_{j_0}):$(obj.CalObj(X)-repmat(obj.CalObj(Dec),size(X,1),1))',$x_{j_0}$表示上一迭代的值；
3、 $dx:$ Dec.*1e-6

$X$变量：

$Dec$变量：

%% 约束求导
function ConGrad = CalConGrad(obj,Dec)
        %CalConGrad - Calculate the gradients of constraints of a solution.
        %
        %   Grad = obj.CalConGrad(Dec) returns the gradients of constraints
        %   of Dec, i.e., a Jacobian matrix.
        %
        %   This function is usually called by gradient-based algorithms.
        %
        %   Example:
        %       ConGrad = Problem.CalConGrad(Dec)
        
            Dec(Dec==0) = 1e-12;
            X           = repmat(Dec,length(Dec),1).*(1+eye(length(Dec))*1e-6);
            ConGrad     = (obj.CalCon(X)-repmat(obj.CalCon(Dec),size(X,1),1))'./Dec./1e-6;
        end

注：考虑到搜索空间的第 j 维的范围，在实验中将$\epsilon$设为$10^{-6}(u_j-l_j)$，其中$u_j$和$l_j$分别表示第 j 个决策变量的下界和上界。

%%共轭梯度计算
  [gk,site] = FiniteDifference(Problem,Population(i),W(i,:));
                    if K(i) == 1
                        dk = -gk;
                    else
                        beta = (gk'*gk)/(g0{i}'*g0{i});
                        dk   = -gk + beta*d0{i};
                        if gk'*dk >= 0
                            dk = -gk;
                        end
                    end   

function [df,site] = FiniteDifference(Problem,X,W)
    if any(X.con>0)
        df   = Problem.CalConGrad(X.dec)';
        site = false(1,length(X.dec));
        df   = sum(df,2);
    else
        df   = Problem.CalObjGrad(X.dec)';
        site = (any(df<0,2)&any(df>0,2))';
        df   = df*W';
    end
end   

3.2.4 子代生成

生成子代过程实则就是差分进化中变异过程与共轭梯度修改决策变量的过程，本文中对差分进化的变异进行了部分修改，原基向量并未采用差分进化中的方式，而是采用了当前个体，此外，差分进化中交叉也省略了（个人认为差分进化的修改均是为了能与共轭梯度能相对应）。首先，还是让我们来看下该部分的伪代码：

🐉🐉🐉🐉🐉🐉🐉🐉 Algorithm 4：Generate Solutions 🐉🐉🐉🐉🐉🐉🐉 🐉
🐉行号	🐉解释
Input	当前父代个体x，Archive集合，搜寻方向（梯度下降方向）s，diff（梯度标志符）
Output	子代y，是否产生更优解标识符success
2	m是为了确定一个合适的步长而设定的，文章使用了一个简单的线搜索策略，步长随着m从1开始指数级地减小，直到找到一个更好的解决方案
3~9	其中s是FRCG中使用的搜索方向，x'，x''是从差分进化中使用的存档中随机选择的两个解；若种群个体搜索方向不一，即diff=1，则采用差分进化方式，否则，则采用共轭梯度下降方法。各目标搜索方向相同的决策变量通过共轭梯度更新以增强种群收敛性(第7行)，而diff中的决策变量通过差分进化更新以增强种群多样性(第9行)。
10~14	如果y的约束违反程度小于解x，则解y优于解x；若两者约束违背相同，则看两者的目标函数值，即f(y)≺f(x)（注：f(y)≺f(x)是y中所有目标函数值均不大于x,即两者存在弱支配关系，具体内容可参照我的一篇博文：多目标非支配排序遗传算法-NSGA-II（一)）

附上代码：

 for step = 0 : 9
    mu        = rand(1,Problem.D) < 1/sum(site);
    OffDec    = Population(i).dec + ~site.*0.5^step.*dk' + mu.*site*0.5^step.*(Archive(randi(end)).dec-Archive(randi(end)).dec);
    Offspring = Problem.Evaluation(OffDec);
    OffPop    = [OffPop,Offspring];
   if sum(max(Offspring.con,0))<sum(max(Population(i).con,0)) || sum(max(Offspring.con,0))==sum(max(Population(i).con,0)) && all(Offspring.obj<Population(i).obj)
         success = true;
          break;
     end
 end

终于完成了！祝大家龙年快乐！