利用遗传算法找极小值_工程优化设计与Matlab实现——基本遗传算法求解一维问题...

亦凉 2023-01-06 04:26 153阅读 0赞

*从猴子到人，本来就是一个优化过程。*

## 遗传算法 ##

## Part 1 遗传算法的自然界原型 ##

遗传算法根据进化论的思想，模拟了生物界遗传选择和自然淘汰的生物进化过程。

### Point 1 从生物遗传谈起 ###

先回想一下生物遗传知识：

*  基因是遗传物质最基本的结构和功能单位
 *  基因在染色体上
 *  基因是有遗传效应的DNA片段

基因是如何在生物之间遗传的：

*  假设某生物种群有n个个体；
 *  其中第n个个体有m个染色体Cn1、Cn2，……，Cnm；
 *  第n个个体的第m个染色体上有k个基因Gnm1，Gnm2，……，GnmK，

如下图：

![23745117e752e03afa92dda04fee699c.png][]

以P1和P2产生p1的过程为例，它们在繁衍后代的时候，对应的染色体（C11和C21）如两两之间可能发生某些基因的交换，也可能发生某些基因的变异，然后产生新的后代染色体。也就是说后代的染色体是在两个父本染色体的基础上，外加交叉变异的影响得到的。

### Point 2 自然选择 ###

从父本产生了新的子代种群后，子代种群不一定能全部存活到产生下一代，环境的适应性和一些意外情况会对它们进行自然选择。总的来说，适应性差的被淘汰的概率就高，适应性强的被淘汰的概率就低。

### Point 3 回到我们要解决的问题 ###

在启发式算法中，我们是用多个迭代点并行进行迭代计算的，每一代的所有迭代点构成了一个种群，迭代一次后得到的点就是后代群体。

![2a0759cb9ffe234c3254ec70a63c770b.png][]

对于m维问题，我们可以将每个点的m个坐标分量视为其染色体。

![595c2710bba5f5d372759e6af3d24eaf.png][]

我们可以用编码的方式将m个坐标分量（视为染色体）进行离散化，使离散后的各个元素所为每个坐标分量的基因。

![40842f6c9a9264a56b8d640a786ca055.png][]

一个群里中的两个迭代点繁衍后代时，对应的坐标分量（视为染色体）中的离散元素（基因）相互交叉，基因变异后组成新的坐标分量（新染色体），所有的新坐标量可以组成新的迭代点（后代个体），所有的后代个体组成了新的后代群体。

也就是说，我们在优化问题中找到对应的元素，套用自然界中的迭代方法即可，如下图。

![d2f93cf6be9aeff579736b52a7b8705e.png][]

我们只是对第一个图的左侧的概念进行了对应替换，就得到了遗传算法的迭代思路。

所以我们完全可以模拟生物进化的方式来更新、筛选我们的迭代点。

## Part 2 基本遗传算法求解一维问题的实现 ##

本篇中我们仅对简单一维问题的遗传算法作出分析。

由于我们将迭代点的每一个坐标分量视为染色体，所以一维问题中每个迭代点只有一个坐标分量。相当于每个个体只有一个染色体，这就使得上面的问题得以简化，如下图：

![ed77fd36464f632f04405d7d6e6a4d21.png][]

对与一维问题，由于一个个体只有一个染色体，所以许多人写的博客或者专栏里都将一维问题的染色体与个体混为一谈，也没有什么问题。

### **Point 1 确定个体、种群、染色体** ###

首先我们先创建由

![Image 1][] 个迭代点组成的一个种群

![Image 1][] ，种群里的每一个体都携带一个染色体（一维问题中每一个点只有一个坐标）。由此我们可以建立一个矩阵：

![cce139afc264c99c2ba5ba6e10c79ae1.png][]

### **Point 2 确定基因：二进制编码** ###

之前我们反复提高了将坐标分量的离散化。

**为什么要离散化？**

每一个染色体是以更小的单元——基因来参与交叉和变异的，所以需要找到构成染色体的基本元素。

在实际问题中，每一个染色体（坐标分量）都是一个常数，我们很难找到一个常数中更小的元素进行交叉或者变异，所以要对染色体进行离散化（编码）。这里我们介绍二进制编码。

**如何实现？**

二进制编码的思想是利用二进制数只有0和1组成且位数众多的特点，使用

![Image 1][] 位二进制数代表更新前的实数，这样就可以对二进制数进行交叉、变异等更新操作，最后再将更新后的二进制数转化为十进制，即为更新后的点。

![be8f9788128ec0f8ce4e42897db8538c.png][]

那么就要建立实数与二进制数之间的映射关系，也就是说在解空间内找任何一个点，都能在编码空间内找到与之对应的点，这正是因为如此，编码空间里的容量一定大于等于解空间的容量。

![20b78ddd02b07a6c0fb8b40b20acbd34.png][]

一方面，如果

![Image 1][] 的取值范围是下界

![Image 1][] 到上界

![Image 1][] ，我们的计算精度至少为0.001，所以至少要从（0,1）区间内划分1001个点（包括端点0、1）。

另一方面，

![Image 1][] 位二进制数可以表示

![Image 1][] 个二进制数。如果找到

![Image 1][] 使

![Image 1][] 即可。由于

![Image 1][] ，所以我们用10位二进制数可以表示1024个点，比1001个点更多，完全可以满足精度。

如果我们让这1024个点按以下方式对应：

![2e8867d53cdb98ada1b1dea2d331b100.png][]

这个映射关系如何表达呢？

![bc056fe420d37692db4f1fedac1b7964.png][]

以上就是我们编码公式，其对应的解码公式为：

![4af3d4b6fb42e639e38ee3c8dc7f51bb.png][]

二进制编码的稳定性高，种群多样性大，但需要的存储空间大，需要解码且难以理解。

**具体实现：**

对于编码长度，我们可以根据上界、下界和精度要求推出来：

%计算由精度与上下界决定的编码长度
    l_Coding = ceil(log2((ub - lb) / E + 1));        %编码长度

编码:我们用一个个体个数行、编码长度列的矩阵X\_Coding来存储编码值，每一列代表一个染色体：

%**********************************************************************%
    %****************************二进制编码函数****************************%
    %**********************************************************************%
    function X_Coding = Coding(n_Population, lb, ub, l_Coding, X)
        for i = 1:n_Population
            temp = X(i);                                                 %十进制原值，临时变量减少重复索引
            xb_dec = round((2^l_Coding - 1) * (temp - lb) / (ub - lb));  %十进制编码值，四舍五入取整round(相当于第几个二进制数)
            xb_bin = dec2bin(xb_dec, l_Coding);                          %二进制编码值,返回的是字符串
            for ii = 1:l_Coding
                X_Coding(ii,i) = str2num(xb_bin(ii));                    %将二进制编码值转换为列向量存储到X_Coding
            end
        end
    end

解码:

function X = DeCoding(n_Population, lb, ub, l_Coding, X_Coding)
        for i = 1:n_Population
            xb_bin = num2str(X_Coding(:,i)');                             %二进制编码值
            xb_dec = bin2dec(xb_bin);                                     %十进制编码值
            X(i) = xb_dec * (ub - lb) / (2^l_Coding - 1) + lb;            %十进制无符号原值
        end
    end

### Point 3 交叉规则 ###

产生一个随机数，如果该随机数小于交叉概率则交叉计算。（轮盘赌法，轮盘只有两个分区）

基本的遗传算法采用相邻两个个体的染色体交叉。使用一个随机数确定交叉的基因片段的起始基因位置，将起始基因到末尾基因之间的所有基因全部交换。（称为随机单点交叉）

![7a8ada6aee4fa82f47aa080b3a680788.png][]

%**********************************************************************%
    %********************相邻两个染色体，随机单点交叉函数********************%
    %**********************************************************************%
    function X_Coding = Cross(n_Population, l_Coding, X_Coding, pc)
        if mod(n_Population ,2) == 0
            cross_num = n_Population;        %如果个体数为偶数，参与交叉个体数为全部个体数
        else
            cross_num = n_Population - 1;    %如果个体数为奇数，参与交叉个体数为全部个体数-1，,最后一个不参与交叉
        end
        for i = 1:2:cross_num
            hand_of_god = rand;            %上帝之手随机做选择
            if hand_of_god < pc            %如果上帝手小于交叉概率pc，就进行交叉，否则不交叉
                cross_point = 2 + round((l_Coding - 2)*rand);    %随机生成一个2到编码长度l_Coding区间内的整数为交叉点的位置（0的话没办法索引，1的话相当于整体换位置）
                Coss_Temp = X_Coding(cross_point:l_Coding, i);   %Coss_Temp为交叉临时中间列向量
                X_Coding(cross_point:l_Coding, i) = X_Coding(cross_point:l_Coding, i+1);
                X_Coding(cross_point:l_Coding, i+1) = Coss_Temp;
            end
        end
    end

### Point 4 变异规则 ###

产生一个随机数，如果该随机数小于交叉概率则变异计算。（轮盘赌法，轮盘只有两个分区）

基本的遗传算法采用随机单点变异，即使用一个随机数确定变异的基因的位置，将该基因变异。

%**********************************************************************%
    %****************************随机单点变异函数****************************%
    %**********************************************************************%
     function X_Coding = mutate(n_Population, l_Coding, X_Coding, pm)
        for i = 1:n_Population
            hand_of_god = rand;            %上帝之手随机做选择
            if hand_of_god < pm            %如果上帝手小于变异概率pc，就进行变异，否则不变异
                mutate_point = 1 + round((l_Coding - 1)*rand);    %随机生成一个1到编码长度l_Coding区间内的整数为变异点的位置（0的话没办法索引）
                temp = X_Coding(mutate_point, i);    %临时变量存储当前变异值，避免反复索引
                if temp == 0
                    temp = 1;
                else
                    temp = 0;
                end
                X_Coding(mutate_point, i) = temp;    %变异完成
            end
        end
     end

### Point 5 自然选择规则 ###

自然选择采用轮盘赌法，即适应度大的被选择的概率高，适应度小的被选择的概率低。

轮盘赌法详见：工程优化设计与Matlab实现——蚁群算法（一）求解旅行商问题

唯一注意的是，这里为了加快收敛速度，每次都将种群中最优解保留了下来，称为精英保留策略。

%**********************************************************************%
    %********************************选择函数********************************%
    %**********************************************************************%
    function X = Select(n_Population, X)
        X_Fitness = max(func(X)) * ones(size(X)) - func(X);    %由于本程序求极小值，函数值与最大函数值的差越大，适应度越大
        fitness_sum = sum(X_Fitness);                          %适应度求和
        Possibility = cumsum(X_Fitness / fitness_sum);         %根据适应度计算，X_Fitness / fitness_Sum是被选择到的概率，Possibility为累积概率
        for i = 1:(n_Population - 1)                           %n_Population选n_Population-1，可重复选（就算全选了一样的也没关系，后面还有变异）
            hand_of_god = rand;                                %“上帝之手”转动轮盘，皮一下好开心
            for ii = 1:n_Population
                if hand_of_god < Possibility(ii)
                    Temp_X(i) = X(ii);                         %Temp_X为临时存储选择到的值的矩阵
                    break
                end
            end
        end
        [~, max_x_position] = max(X_Fitness);                 %找出最大适应度
        Temp_X(n_Population) = X(max_x_position);
        X = Temp_X;                                           %得到选择后的种群
    end

### Point 6 迭代流程 ###

![7351f6af052af397546b8299a5a450fa.png][]

## Part 3 基本遗传算法求解一维问题的栗子 ##

利用基本遗传算法求解

![Image 1][] 的极大值，

![Image 1][]

目标函数定义如下，由于遗传算法的程序是求极小值点的程序，在这里求极大值点需要取相反数函数：

function f = func(x)
        f = x + 10 .* sin(5 .* x) + 7 .* cos(4 .* x);
        f = - f;

基本遗传算法求解：

%----------基本遗传算法求函数极小值----------%
    clc;
    clear all;
    close all;
    %----------设置参数----------%
    n_Population = 5;    %种群个体数目
    ub = 10;              %给定上界
    lb = -10;             %给定下界
    E = 0.01;          %给定精度要求
    pc = 0.4;            %交叉概率
    pm = 0.1;            %变异概率
    n_loop_max = 200;     %最大迭代次数
    
    %计算由精度与上下界决定的编码长度
    l_Coding = ceil(log2((ub - lb) / E + 1));        %编码长度
    
    %----------一些矩阵说明----------%
    X = zeros(1, n_Population)*nan;                       %X为种群矩阵，1行n_Population列
    X_Coding = zeros(l_Coding, n_Population)*nan;         %X_Coding为种群编码矩阵，l_Coding行n_Population列；
    X_Fitness = zeros(1, n_Population)*nan;               %X_Fitness适应度矩阵，1行n_Population列
    Possibility = zeros(1, n_Population)*nan;             %Possibility累积概率概率矩阵，1行n_Population列
    History_X = zeros(n_loop_max, n_Population)*nan;      %X为种群矩阵历史记录，最终为n_loop_max+1行n_Population列
    History_FX = zeros(n_loop_max, n_Population)*nan;      %X为种群函数值矩阵历史记录，最终为n_loop_max+1行n_Population列
    History_X_Fitness = zeros(n_loop_max, n_Population)*nan; %X为种群适应度矩阵历史记录，最终为n_loop_max行n_Population列
    
    
    for n_loop = 1:n_loop_max
        if n_loop == 1
            X = X_init(n_Population,lb,ub);           %种群初始化
        end
        
        %将当前种群状态存储进历史记录
        History_X(n_loop, :) = X;                     %存储当前种群
        History_X_Fitness(n_loop, :) = max(func(X)) * ones(size(X)) - func(X);     %存储当前适应度,由于本程序求极小值，函数值与最大函数值的差越大，适应度越
        
        fprintf('第%d次迭代n',n_loop)
        fprintf('初始种群：')
        display(X)
        fprintf('初始种群适应度：')
        display(X_Fitness)
        
        %编码
        X_Coding = Coding(n_Population, lb, ub, l_Coding, X);
        fprintf('种群编码：')
        display(X_Coding)
        
        %交叉
        X_Coding = Cross(n_Population, l_Coding, X_Coding, pc);
        fprintf('交叉：')
        display(X_Coding)
        
        %变异
        X_Coding = mutate(n_Population, l_Coding, X_Coding, pm);
        fprintf('变异：')
        display(X_Coding)
        
        %解码
        X = DeCoding(n_Population, lb, ub, l_Coding, X_Coding);
        fprintf('解码后得到新种群：')
        display(X)
        
        %选择
        X = Select(n_Population, X);
        fprintf('选择后的新种群：')
        display(X)
        
        if n_loop == n_loop_max
            History_X(n_loop_max+1,:) = X;              %将最后一次的种群存入历史记录
        end
    end 
    
    [min_f, min_x_position] = min(func(X));    %找出最小函数值
    min_x = X(min_x_position);                 %找到最小函数值对应的个体
    
    %**********************打印结果**********************%
    fprintf('二进制编码遗传算法：n')
    fprintf('****************参数设置****************n')
    fprintf('种群个体数目: %dn', n_Population)
    fprintf('给定上界: %f      给定下界：  %fn', ub, lb)
    fprintf('精度为：  %fn', E)
    fprintf('编码长度: %dn', l_Coding)
    fprintf('交叉概率: %fn', pc)
    fprintf('变异概率: %fn', pm)
    fprintf('最大迭代次数: %dn', n_loop_max)
    fprintf('****************迭代结果****************n')
    fprintf('共迭代 %d 次n', n_loop)
    fprintf('极小值点为： %fn', min_x)
    fprintf('极小值为： %fn', min_f)
    fprintf('****************************************n')
    fprintf('endn')
    
    % %**********************作图**********************%
    % History_FX = func(History_X); 
    % for i = 1 : n_loop_max+1
    %     figure(i)
    %     xx = lb:0.1:ub;
    %     yy = func(xx);
    %     plot(xx,yy);
    %     hold on
    %     plot(History_X(i,:),History_FX(i,:),'r*');
    %     saveas(figure(i),['C:UsersAirHeDesktop1',int2str(i),'.jpg'])
    %     close(figure(i))
    % end
    
                                                                               
    %**********************************************************************%
    %****************************种群X初始化函数****************************%
    %**********************************************************************%
    function X = X_init(n_Population, lb, ub)
        for i = 1:n_Population               %随机生成位于上下界内的种群
            X(i) = lb + rand * (ub - lb);
            if X(i) > ub        %如果圆整后超过上界则纠正到上界
                X(i) = ub;
            elseif X(i) < lb    %如果圆整后超过下界则纠正到下界
                X(i) = lb;
            end
        end
    end
    
    %**********************************************************************%
    %****************************二进制编码函数****************************%
    %**********************************************************************%
    function X_Coding = Coding(n_Population, lb, ub, l_Coding, X)
        for i = 1:n_Population
            temp = X(i);                                                 %十进制原值，临时变量减少重复索引
            xb_dec = round((2^l_Coding - 1) * (temp - lb) / (ub - lb));  %十进制编码值，四舍五入取整round(相当于第几个二进制数)
            xb_bin = dec2bin(xb_dec, l_Coding);                          %二进制编码值,返回的是字符串
            for ii = 1:l_Coding
                X_Coding(ii,i) = str2num(xb_bin(ii));                    %将二进制编码值转换为列向量存储到X_Coding
            end
        end
    end
    
    %**********************************************************************%
    %****************************二进制解码函数****************************%
    %**********************************************************************%
    function X = DeCoding(n_Population, lb, ub, l_Coding, X_Coding)
        for i = 1:n_Population
            xb_bin = num2str(X_Coding(:,i)');                             %二进制编码值
            xb_dec = bin2dec(xb_bin);                                     %十进制编码值
            X(i) = xb_dec * (ub - lb) / (2^l_Coding - 1) + lb;            %十进制无符号原值
        end
    end
    
    %**********************************************************************%
    %********************相邻两个染色体，随机单点交叉函数********************%
    %**********************************************************************%
    function X_Coding = Cross(n_Population, l_Coding, X_Coding, pc)
        if mod(n_Population ,2) == 0
            cross_num = n_Population;        %如果个体数为偶数，参与交叉个体数为全部个体数
        else
            cross_num = n_Population - 1;    %如果个体数为奇数，参与交叉个体数为全部个体数-1，,最后一个不参与交叉
        end
        for i = 1:2:cross_num
            hand_of_god = rand;            %上帝之手随机做选择
            if hand_of_god < pc            %如果上帝手小于交叉概率pc，就进行交叉，否则不交叉
                cross_point = 2 + round((l_Coding - 2)*rand);    %随机生成一个2到编码长度l_Coding区间内的整数为交叉点的位置（0的话没办法索引，1的话相当于整体换位置）
                Coss_Temp = X_Coding(cross_point:l_Coding, i);   %Coss_Temp为交叉临时中间列向量
                X_Coding(cross_point:l_Coding, i) = X_Coding(cross_point:l_Coding, i+1);
                X_Coding(cross_point:l_Coding, i+1) = Coss_Temp;
            end
        end
    end
    
    %**********************************************************************%
    %****************************随机单点变异函数****************************%
    %**********************************************************************%
     function X_Coding = mutate(n_Population, l_Coding, X_Coding, pm)
        for i = 1:n_Population
            hand_of_god = rand;            %上帝之手随机做选择
            if hand_of_god < pm            %如果上帝手小于变异概率pc，就进行变异，否则不变异
                mutate_point = 1 + round((l_Coding - 1)*rand);    %随机生成一个1到编码长度l_Coding区间内的整数为变异点的位置（0的话没办法索引）
                temp = X_Coding(mutate_point, i);    %临时变量存储当前变异值，避免反复索引
                if temp == 0
                    temp = 1;
                else
                    temp = 0;
                end
                X_Coding(mutate_point, i) = temp;    %变异完成
            end
        end
     end
    
    %**********************************************************************%
    %********************************选择函数********************************%
    %**********************************************************************%
    function X = Select(n_Population, X)
        X_Fitness = max(func(X)) * ones(size(X)) - func(X);    %由于本程序求极小值，函数值与最大函数值的差越大，适应度越大
        fitness_sum = sum(X_Fitness);                          %适应度求和
        Possibility = cumsum(X_Fitness / fitness_sum);         %根据适应度计算，X_Fitness / fitness_Sum是被选择到的概率，Possibility为累积概率
        for i = 1:(n_Population - 1)                           %n_Population选n_Population-1，可重复选（就算全选了一样的也没关系，后面还有变异）
            hand_of_god = rand;                                %“上帝之手”转动轮盘，皮一下好开心
            for ii = 1:n_Population
                if hand_of_god < Possibility(ii)
                    Temp_X(i) = X(ii);                         %Temp_X为临时存储选择到的值的矩阵
                    break
                end
            end
        end
        [~, max_x_position] = max(X_Fitness);                 %找出最大适应度
        Temp_X(n_Population) = X(max_x_position);
        X = Temp_X;                                           %得到选择后的种群
    end

计算结果：

![e3d88c9fa49bfe98b978933c24354569.png][]

遗传算法中种群数量、交叉概率、变异概率、精度要求都会对算法的性能造成影响。

遗传算法由于在较小的区间内搜索时，交叉对个体影响很小，变异影响较大，所以小区域内搜索能力差。

另外，由于编码的特点，区间内的点仍然有些无法完全覆盖，也会对求解精度产生影响。

[23745117e752e03afa92dda04fee699c.png]: /images/20221119/43340354244c4fa4a4976522e80e5eb2.png
[2a0759cb9ffe234c3254ec70a63c770b.png]: /images/20221119/1139f6a407a94bacb860f07c6591cced.png
[595c2710bba5f5d372759e6af3d24eaf.png]: /images/20221119/04139011ae8744709aafc010f6a9f6a8.png
[40842f6c9a9264a56b8d640a786ca055.png]: /images/20221119/83281c83c1d04762a05769b4dff1cb37.png
[d2f93cf6be9aeff579736b52a7b8705e.png]: /images/20221119/ac55edeae7364c6faf63c6014ee0a54d.png
[ed77fd36464f632f04405d7d6e6a4d21.png]: /images/20221119/b3b99ca637fc4929bc4dd2b014456b1a.png
[Image 1]: 
[cce139afc264c99c2ba5ba6e10c79ae1.png]: /images/20221119/555bd6101cd449259d6e45e734dd9324.png
[be8f9788128ec0f8ce4e42897db8538c.png]: /images/20221119/cf731b11e33b468586d4bcf012014248.png
[20b78ddd02b07a6c0fb8b40b20acbd34.png]: /images/20221119/45f7af068e6a41fcb2534a3b1669dab5.png
[2e8867d53cdb98ada1b1dea2d331b100.png]: /images/20221119/37be710205144a1cb5d90a0f63540ade.png
[bc056fe420d37692db4f1fedac1b7964.png]: /images/20221119/987d190ab4c3466f97a700e55cfa2250.png
[4af3d4b6fb42e639e38ee3c8dc7f51bb.png]: /images/20221119/116560d6a6704ef8a4b33f1b0ab88873.png
[7a8ada6aee4fa82f47aa080b3a680788.png]: /images/20221119/17d92af32bc74b15854a90468e01f806.png
[7351f6af052af397546b8299a5a450fa.png]: /images/20221119/ec33a1e0ba98495aa944973b8a068011.png
[e3d88c9fa49bfe98b978933c24354569.png]: https://img-blog.csdnimg.cn/img_convert/e3d88c9fa49bfe98b978933c24354569.png