【学习笔记】基于遗传算法和非线性规划的函数寻优算法

本文是《MATLAB智能算法30个案例分析（第二版）》一书第二章的学习笔记。

一、背景介绍

标准遗传算法的主要本质特征，在于群体搜索策略和简单的遗传算子，这使遗传算法获得了强大的全局最优解能力、问题域的独立性、信息处理的并行性、应用的鲁棒性和操作的简明性。但大量的实践和研究表明，标准的遗传算法存在局部搜索能力差和“早熟”等缺陷，不能保证算法收敛。

现有的许多文献中，针对遗传算法存在的上述问题，改进算法主要集中在对遗传算法的编码机制、选择策略、交叉算子、变异算子、特殊算子和参数设计（包括群体规模、交叉概率、变异概率）等。

经典的非线性规划算法大多采用梯度下降的方法求解，局部搜索能力强而全局搜索能力弱，容易陷入局部最优解。本文结合两种算法的优点，一方面使用遗传算法进行全局搜索，另一方面使用非线性规划算法进行局部搜索，以得到问题的全局最优解。

二、算法流程

1. 种群初始化
2. while  当前迭代次数小于最大迭代次数
3.  计算适应度;
4.  选出最优个体的适应度及其染色体
5.  选择算子;
6.  交叉算子;
7.  变异算子;
8.  if   当前迭代次数满足某个关系
9.   对每个个体进行非线性规划得到局部最优解，并以此结果为新的染色体继续迭代
10.  end
11. end

三、算法实现

以计算函数$$f\left(x\right)=\sum _{i=1}^n{x}_i^2\left(-20⩽x_i⩽20\right)$$
为例，其中，个体$x$的维度$n$为10.这是一个简单的平方和函数，只有一个极小点$x=(0,0,\cdots ,0)$，理论最小值$f\left(0,0,\cdots ,0\right)=0$。

3.1 主函数 main.m

%% 基于遗传算法和非线性规划的函数寻优算法
% 一、结合思想：
% 传统的非线性规划大多采用梯度下降法，局部搜索能力强、全局搜索能力弱
% 遗传算法的全局搜索能力强、局部搜索能力弱

% 二、算法流程
% 种群初始化;
% 计算适应度;
% while 当前进化次数 < 最大进化次数
%     选择;
%     交叉;
%     变异;
%     if 当前进化次数满足某个关系
%         对每个个体，进行非线性规划得到局部最优解，并以此为新的染色体继续进化
%     end
% end

% 三、代码
clear;clc
% (1)种群初始化
% 1.1 初始参数列表
Nind = 20;                                                          % 种群规模
Lind = 10;                                                          % 染色体个数
maxGen = 200;                                                       % 进化次数（迭代次数）
Pc = 0.6;                                                           % 交叉概率
Pm = 0.01;                                                          % 变异概率
Xmin = [-20, -20, -20, -20, -20, -20, -20, -20, -20, -20 ];         % 下界
Xmax = [20, 20, 20, 20, 20, 20, 20, 20, 20, 20];                    % 上界

% 1.2 初始化个体
individuals = struct('fitness', zeros(Nind, 1), 'chrom', []);       % 种群结构体
avg_Fit = [];           % 平均适应度
best_fit = [];          % 最佳适应度
best_chrom = [];        % 最佳适应度对应的个体

% 1.3 初始化种群
for i = 1:Nind
    individuals.chrom(i,:) = rand(1, Lind).*(Xmax - Xmin)+Xmin;     % 随机产生染色体
    temp = individuals.chrom(i,:);
    individuals.fitness(i) = -1*obj_fun(temp);                      % 计算适应度
end
best_chrom = [];                                                    % 历代最优个体
avg_Fit = [];                                                       % 历代适应度平均值
trace = [];                                                         % 历代最优适应度

% (2)开始迭代
for gen = 1:maxGen
    % 计算适应度
    for i = 1:Nind
        temp = individuals.chrom(i,:);
        individuals.fitness(i) = -1*obj_fun(temp);                  % 计算适应度
    end
    
    % 找出最优适应度及其对应的个体
    [best_fit best_index] = max(individuals.fitness);
    best_chrom = [best_chrom;individuals.chrom(best_index,:)];
    avg_Fit = [avg_Fit; mean(individuals.fitness)];
    trace = [trace;best_fit];
    
    % 选择算子
    individuals = Select(individuals,Nind);
    
    % 交叉操作
    individuals = Cross(individuals,Pc,Nind,Lind);
    
    % 变异操作
    individuals = Mutation(individuals,Pm,Nind,Lind,Xmax,Xmin);
    
    % 非线性规划
    if do_NLP(gen) == 1
        A = []; b = []; Aeq = []; Beq = [];
        for i = 1:Nind
            x0 = individuals.chrom(i,:);
            individuals.chrom(i,:) = fmincon(@obj_fun, x0, A, b, Aeq, Beq, Xmin, Xmax);
        end
    end
end
trace = -trace;
avg_Fit = -avg_Fit;
clc
disp(['最优适应度为：',num2str(trace(end))])
disp('对应的种群为：')
disp(best_chrom(find(trace == min(trace))-1,:))
plot(trace, 'LineWidth', 1)
title(['适应度进化曲线（最大迭代次数maxGen=',num2str(maxGen),'）'])
xlabel('进化次数')
ylabel('适应度')
hold on
plot(avg_Fit, 'LineWidth', 1)
legend('历代最佳适应度', '历代平均适应度')

3.2 选择算子 Select.m

本文使用随机联赛选择作为选择算子。

随机联赛选择也是一种基于个体适应度之间大小关系的选择方法。其基本思想是每次随机选取$N$个个体，选出其中适应度最高的一个个体遗传到下一代群体中。

在随机联赛选择操作中，只有个体适应度之间的大小比较运算，而无个体适应度之间的算术运算，所以它对个体适应度去正值还是取负值无特别影响。

随机联赛选择的具体操作过程为：
（1）从群体中随机选取$N$个个体进行适应度大小的比较，将其中适应度最高的个体遗传到下一代群体中。
（2）将上述过程重复$M$次（$M$为种群规模），就可得到下一代群体中的$M$个个体。

function new_individuals = Select(individuals, Nind)
% select 选择算子，基于随机联赛策略选择较优的个体
% individuals        input:     种群
% Nind               input:     种群规模
% new_chrom          output:    经选择算子生成的新种群

N = 2;          % 每次随机选择的个体数
index = [];     % 存储每次选择的个体的索引
for i = 1:Nind
    race_indi_index = randi(Nind, [2,1]);       % 被选中的个体的索引
    race_indi_fit = individuals.fitness(race_indi_index);   % 被选中的个体的适应度
    bigger_index = find(individuals.fitness == max(race_indi_fit));
    index = [index; bigger_index(1)];
end
individuals.chrom = individuals.chrom(index,:);
individuals.fitness = individuals.fitness(index,:);
new_individuals = individuals;

end

3.3 交叉算子 Cross.m

本文使用模拟二进制交叉算子（Simulated Binary Crossover，SBX），主要思想是在使用浮点数编码/实数编码时，可以模拟二进制编码的操作，计算公式如下：
$$\begin{gathered} \widetilde{x_{1j}}=\frac{1}{2}\left[\left(1+r_j\right)x_{1j}+\left(1-r_j\right)x_{2j}\right] \\ \widetilde{x_{2j}}=\frac{1}{2}\left[\left(1-r_j\right)x_{2j}+\left(1+r_j\right)x_{2j}\right] \end{gathered}$$
其中，$x_{1j}$、$x_{2j}$分别表示两个父代第$j$个基因，$r_j$为参数，由下面的公式确定：
$$r_j=\{\begin{array}{l}{\left(2u_j\right)}^{\frac{1}{\eta +1}},u_j⩽0.5\\ {\left(\frac{1}{2\left(1-u_j\right)}\right)}^{\frac{1}{\eta +1}},0.5⩽u_j<1\end{array}$$
其中，$u_j$为$[0,1]$上的随机数，$\eta$为分布指数，Deb和Agrawal建议 $\eta =1$。

function new_individuals = Cross(individuals, Pc, Nind, Lind)
% Cross 使用模拟二进制交叉算子(SBX)，随机选择两个个体进行染色体的交叉
% individuals           input:      种群
% Pc                    input:      交叉概率
% Nind                  input:      种群规模
% Lind                  input:      染色体长度
% new_individuals       output:     经交叉操作后产生的新种群

new_individuals = individuals;
for i = 1:Nind
    % 随机选择两个个体
    index = randi(Nind, [2,1]);
    cross_indi = individuals.chrom(index, :);
    temp = zeros(2, Lind);
    
    % 进行交叉操作
    for j = 1:Lind
        Pick = rand;
        if Pick < Pc
            eta = 0.5;          % 分布因子
            % 分布因子越大，其产生的后代个体逼近父代个体的概率越大，故SBX算子在局部搜索上表现较强
            u = rand;
            if u <= 0.5
                beta = (2*u)^(1/(1+eta));
            else
                beta = (2*(1-u))^(1/(1+eta));
            end
            
            % 交叉
            temp(1,j) = 0.5*((1+beta)*individuals.chrom(index(1),j) + (1-beta)*individuals.chrom(index(2),j));
            temp(2,j) = 0.5*((1-beta)*individuals.chrom(index(1),j) + (1+beta)*individuals.chrom(index(2),j));
        else
            temp(1,j) = individuals.chrom(index(1),j);
            temp(2,j) = individuals.chrom(index(2),j);
        end
    end
    new_individuals.chrom(index(1),:) = temp(1,:);
    new_individuals.chrom(index(2),:) = temp(2,:);
end

end

3.4 变异算子 Mutation.m

本文使用高斯变异算子。高斯变异是改进遗传算法对重点搜索区域的局部搜索性能的另一种变异操作方法。所谓高斯变异操作，是指在进行变异操作时，用符合均值为$\mu$、方差为${\sigma }^2$的正态分布的一个随机数来替换原有基因值。

具体实现高斯变异时，一个符合正态分布的随机数$Q$可由一些符合均匀分布的随机数利用公式来近似产生。假定有12个在$[0,1]$范围内均匀分布的随机数$r_i(i=1,2,\cdots ,12)$，则符合$N(\mu ,{\sigma }^2)$正态分布的一个随机数可由下式求得：
$$Q=\mu +\sigma \cdot \left(\sum _{i=1}^{12}{r}_i-6\right)$$
在进行由$X=x_1{x}_2\cdots x_k\cdots x_l$向$X=x_1{x}_2\cdots x_k^{{}^{\prime }}\cdots x_l$的高斯变异操作时，若变异点$x_k^{{}^{\prime }}$处的基因值取值范围为$[U_{min}^k,U_{max}^k]$，并假设：
$$\begin{gathered} \mu =\frac{U_{\min }^k+U_{\max }^k}{2} \\ \sigma =\frac{U_{\max }^k-U_{\min }^k}{6} \end{gathered}$$
则新的基因值$x_k^{{}^{\prime }}$可由以下公式确定：
$$x_k^{{}^{\prime }}=\frac{U_{\min }^k+U_{\max }^k}{2}+\frac{U_{\max }^k-U_{\min }^k}{6}\left(\sum _{i=1}^{12}{r}_i-6\right)$$

function new_individuals = Mutation(individuals, Pm, Nind, Lind, Xmax, Xmin)
% Mutation 对种群中每个个体，使用高斯变异算子，随机选取一点进行变异
% individuals           input:      种群
% Pm                    input:      变异概率
% Nind                  input:      种群规模
% Lind                  input:      染色体长度
% Xmax                  input:      自变量的上界
% Xmin                  input:      自变量的下界
% new_individuals       output:     经交叉操作后产生的新种群

new_individuals = individuals;
for i = 1:Nind
    for j = 1:Lind
        Pick = rand;
        if Pick < Pm
            u = (Xmax(j)+Xmin(j))/2;
            sigma = (Xmax(j)-Xmin(j))/6;
            new_individuals.chrom(i,j) = u+sigma*(sum(rand([1,12]))-6); 
        end
    end
end

end

3.5 非线性规划函数 do_NLP.m

本文的处理方法比较简单：若当前迭代次数$gen$为20的整数倍，则返回1，表示使用MATLAB内置的$fmincon$函数进行非线性规划；反之则返回0，不进行非线性规划。

function flag = do_NLP(gen)
% do_NLP 判断当前迭代次数是否满足进行非线性规划的条件
% gen                       input:      当前迭代次数
% flag                      output:     是否可行，1为可进行，0位不可进行
if mod(gen,20) == 0
    flag = 1;
else
    flag = 0;
end
end

3.6 目标函数 obj_fun.m

function fitness = obj_fun(chrom)
% OBJ_FUN 计算染色体序列Chrom对应的染色度
% chrom        input:   染色体序列
% fitness      output:  适应度
x = chrom;
fitness = sum(chrom.^2);
end

四、运行结果