遗传算法Matlab代码实现及算法函数封装

---

前言

遗传算法（GA）作为求解单目标优化问题的有效算法，自提出以来，便被广泛采用。该算法主要是模仿达尔文进化论，通过种群不断的进行自然选择、繁衍交叉变异，最终得到最优个体的过程。本文将详细介绍遗传算法的流程及各步骤的代码实现，并给出3个实例进行展示，最后给出已封装函数的参数介绍与使用示例。

---

一、遗传算法介绍

遗传算法最早是美国的Holland教授于1975 年在其专著中提出，是一类借
鉴生物界自然选择和自然遗传机制的随机搜索算法。

遗传算法大致可以分为以下数个步骤：种群初始化、适应度计算、种群初步筛选、筛选后交叉操作、筛选后变异操作、种群合并、适应度计算、种群筛选、算法终止。算法流程详情如下图所示：

二、遗传算法算子

从上图很明显能够发现，遗传算法主要是由初始化、选择、交叉、变异数个步骤构成，接下来本文将结合matlab代码对各个步骤进行详细介绍。

1.种群初始化

依据算法流程，首先需要对种群进行初始化，该过程依据随机采样的原则，在目标函数的定义域范围内随机选取数个点，以此作为第一代种群（F1）。

该步骤的关键的在于种群初始化时，编码方式的选择，不同的编码方式所产生的效果不同，且各有利弊。常见的编码方式有：二进制数编码、格雷码编码、浮点数编码。本文主要介绍二进制数编码与浮点数编码。

1.1 二进制数编码

顾名思义，该种编码方式是利用二进制数0、1进行对种群个体的刻画。以函数f(x,y)=sin(x+y)为例，分别取x、y定义域为[a,b]、[c,d]，并令编码长度为4与6（可依据实际情况进行自定义），对于自变量x来说，二进制数0000就代表定义域内a的值，二进制数1111就代表b，定义域内其他的数，便由其余4位二进制数代替，另一自变量类似。代码实现如下所示：

%% 双自变量初始化如下所示
% x的定义域为[a,b]或者(a,b)
long_chromo=[4,6];
% 更多自变量时可以令：
% long_chromo=[4,3,7];
% 具体取值可以根据具体问题可以进行自主更改
chromo_x=zeros(1,4+6);
% 代表  全部自变量的二进制编码组成的最终基因编码，
% 多个自变量时需按照各个自变量的基因长度，编写对应的解码函数
long_chromo1=[0,long_chromo];
for i=1:length(long_chromo)
	for j=1:long_chromo(i)
		chromo_x(sum(long_chromo1(1:i))+j)=round(rand(1));
	end
end
% 随机生成二进制数组成两个自变量的基因编码

%% 双自变量解码如下所示：
x1_min=a;x1_max=b;
x2_min=c;x2_max=d;
% 自变量上下限
t1=3:-1:0;
t2=5:-1:0;
idx1=2.^t1;
idx2=2.^t2;
x1_10=chromo_x(1:4).*idx1;
x2_10=chromo_x(5:10).*idx2;
% 位权法计算四位二进制数对应的十进制数
x1=(x1_max-x1_min)*(x1_10)/(2^3+2^2+2^1+2^0)+x1_min;
x2=(x2_max-x2_min)*(x2_10)/(2^5+2^4+2^3+2^2+2^1+2^0)+x2_min;
% 将十进制数映射到自变量的区间上，之后才可带入目标函数进行计算

很明显，由二进制编码构成的基因在后续使用前，需要一个解码过程，该过程用的是位权法及区间映射方法。易知，位数确定的二进制数显然是有限个，而自变量定义域是连续时，是有无限个，因此很难构建两者的一个双射或者满射，即二进制编码会漏算不少的自变量。我们只能通过加长单个自变量的基因编码长度而减小该问题的影响。
但是，二进制编码在处理离散的自变量时，尤其是面对0-1规划模型时，有着较好的求解结果。

1.2 浮点数编码

相较于二进制数编码，该方法极其直接、方便，依旧以函数f(x,y)=sin(x+y)为例，分别取x、y定义域为[a,b]、[c,d]，由于仅有两个自变量，故基因编码长度为2，具体代码如下所示：

% 双自变量，浮点数基因编码初始化
chromo_x=zeros(1,2);
x_min=[a,c];
x_max=[b,d];
for i=1:2
	chromo_x(i)=(x_max(i)-x_min(i))*rand(1)+x_min(i);
end

浮点数编码相较于二进制数编码、格雷码编码最大的优点在于：无需解码，使用方便、应对连续自变量的优化模型时求解效果极好。但是在处理整数规划问题时，上述代码还需进行进一步的修改。

1.3 小结

不同的编码方式应对不同种情况，各有利弊。除去上文两种示例编码方式之外，常用的还有格雷码编码方式，该格雷码本质上是类二进制数，但是该编码方式相较于二进制编码更为连续、更不容易“犯错”，读者感兴趣可以自行查阅相关文献进行了解。
上文中的示例代码均仅展示如何初始化单个个体，如果想要进行种群初始化，仅需将上述代码套入一个循环即可（或者运用矩阵运算化简掉循环），详见下例：

% pop为种群大小
pop=100;
% x_num为自变量个数
x_num=2;
% 浮点数编码种群初始化
chromo=zeros(pop,x_num);
for i=1:pop
    for j=1:x_num
        chromo(i,j)=x_min(j)+(x_max(j)-x_min(j))*rand(1);
    end
end
% 二进制数编码种群初始化
chromo=zeros(pop,sum(long_chromo));
long_chromo=[0,long_chromo];
for i=1:pop
    for j=1:x_num
        chromo(i,long_chromo(j)+1:sum(long_chromo(1:j))+long_chromo(j+1))=round(rand(1,long_chromo(j+1)));
    end
end

种群大小决定遗传算法的搜索规模，种群越大，搜索到最优解的自变量的可能就越大，但是种群过大会降低代码运行效率（尤其是在相对复杂的问题中），一般可以针对问题，多次运行代码尝试出合适的种群大小。

至此，完成遗传算法中的种群初始化。

2.选择算子

依据流程图，种群初始化后，需要计算适应度，然后依据适应度进行种群的选择操作。简单来说，适应度的计算完全可以采用目标函数来计算，但是相应的，一定要记清楚优化方向：是最小化，还是最大化。这个直接影响到选择算子的代码编写。

选择算子主要是实现对种群的筛选，即：剔除掉劣个体（劣解），保留下优个体（优解）的过程。

常用的选择方法有：轮盘赌、锦标赛选取法等等，详细介绍可以参考文章：

https://blog.csdn.net/hba646333407/article/details/103251008

本文以锦标赛选取法为例，进行代码实现。
锦标赛选取法具体指的是，我们从种群中随机选择数个个体（随机抽样且放回），再从中选出适应度最好的个体（数个个体中最优的个体），放入下一代中。以求目标函数最大值为例：

function pop_choosed=choose_pop(chromo_value,num)
% chromo_value为种群中各个个体的基因编码与适应度（目标函数值）
% 其大小为一个矩阵，行数为个体数，列数为基因长度+1，即最后一列为各个个体的适应度
pop_choosed=zeros(size(chromo_value));
pop=length(chromo_value(:,1));
for i=1:pop
    index=randperm(pop,num);
    % 随机选取num个个体
    r=chromo_value(index,end);
    [~,index_max]=max(r);
    % 找出适应度最优的个体编号
    pop_choosed(i,:)=chromo_value(index(index_max),:);
    % 保留该个体的基因编码及适应度
end
end

上述代码中num为随机抽样的规模大小，需要自行依据种群大小手动设定，其值越大，则说明选出来的解具有一个较大的适应范围（解为最优解的可能越大），其值越小，则说明选出来的解只有较小的适应范围（解为最优解的可能越小）。

但该值并不是越大越好，值太大，算法容易陷入局部最优解，值太小，降低搜索最优解的效率，因此该参数的具体取值可以通过多次运行来寻找一个合适的值。考虑到该参数的特点，本文构建了一种自适应参数锦标赛选取法作为封装函数中的选择算子，详见后文。

3.交叉算子

经过选择操作之后，筛选后的种群开始了繁衍，也就是进行交叉操作，但是也有可能不发生。由此我们能够得到交叉算子的一个重要参数：交叉概率。

通常情况下，交叉概率较大，故经常选取0.8、0.9这样的与1相近的数作为交叉概率。

与选择算子类似，交叉算子也有多种方法，具体可以见下链接文章：

https://blog.csdn.net/hba646333407/article/details/103349279

由于本文侧重浮点数编码，故在此介绍模拟二进制交叉方法，即采用SBX算子[2]。该方法本质上是模仿二进制编码中的单点交叉，但是在其基础上又引入另一参数yita_c，该参数值越大，则交叉后的子代基因与父代基因越相似。

（该算子代码与变异算子的代码一同展示）

4.变异算子

种群在交叉繁衍的过程中，基因难免会出现基因突变的情况，但是概率较低。由此能够得到变异算子的一个重要参数：变异概率。

该参数通常取值为0.1、0.05或者种群个体数的倒数。

变异算子也有多种方法，最常见的即单点变异，即对每一个基因生成一个随机数，若小于变异概率，则发生变异（比如说对于二进制数编码，变异操作即，值从0变为1，或者从0变为1，若为浮点数编码，则可以在原数的某邻域内进行随机选择数字作为变异后的结果）。

本文主要采用的方法是多项式变异，具体原理可以从第二篇、第三篇参考文献中找到。该方法最早是Deb（多目标规划经典算法NSGA2的作者）提出的一种变异策略。多项式变异较之于上文提到的变异方法多出了一个分布指数，分布指数的大小一定程度决定了变异后的基因与原基因的差异大小。

模拟二进制交叉与多项式变异代码如下：

function new_pop=cross_mutation(pop_chromo,pc,pm,yita_c,yita_m,x_min,x_max,temp,long_chromo,p_eq,Aeq)
%% 参数说明
% pop_chromo  行数为个体数，列数为基因数，其中不包含价值
% pc          交叉概率
% pm          变异概率
% yita_c      模拟二进制单点交叉中的参数
% yita_m      多项式变异中的参数，为分布参数
% x_min       自变量x的下界，为行向量
% temp        不同编码形式
% long_chromo 每个基因的二进制编码长度

off=zeros(size(pop_chromo));
% 采用的交叉方法为：模拟二进制单点交叉，即SBX算子
num=length(pop_chromo(:,1));
for i=1:num/2
    father_idx=randperm(num,2);
    % 选取两个不同的个体作为父代
    temp_chromo=pop_chromo(father_idx,:);
    off_chromo=zeros(size(temp));
    % 提取基因
    if rand(1)<pc
        % 进行交叉操作
        for j=1:length(pop_chromo(1,:))
            % 对每个基因使用SBX算子进行计算
            r=rand(1);
            if r<=0.5
                beta1=(2*r)^(1/(1+yita_c));
            else
                beta1=(1/(2-2*r))^(1/(1+yita_c));
            end
            off_chromo(1,j)=0.5*((1+beta1)*temp_chromo(1,j)+(1-beta1)*temp_chromo(2,j));
            off_chromo(2,j)=0.5*((1-beta1)*temp_chromo(1,j)+(1+beta1)*temp_chromo(2,j));
        end
        %% 越界处理，该步骤可以更加严格
        if off_chromo(1,j)>x_max(j)||off_chromo(1,j)<x_min(j)
            off_chromo(1,j)=x_min(j)+rand(1)*(x_max(j)-x_min(j));
        end
        if off_chromo(2,j)>x_max(j)||off_chromo(2,j)<x_min(j)
            off_chromo(2,j)=x_min(j)+rand(1)*(x_max(j)-x_min(j));
        end
    else
        % 不进行交叉操作
        off_chromo=temp_chromo;
    end
    % 采用的变异方法为多项式变异
    if rand(1)<pm
        for j=1:length(pop_chromo(1,:))
            r=rand(1);
            if r<=0.5
                deta1=(off_chromo(1,j)-x_min(j))/(x_max(j)-x_min(j));
                beta2=(2*r+(1-2*r)*(1-deta1)^(yita_m+1))^(1/(yita_m+1))-1;
            else
                deta2=(x_max(j)-off_chromo(1,j))/(x_max(j)-x_min(j));
                beta2=1-(2-2*r+(2*r-1)*(1-deta2)^(yita_m+1))^(1/(yita_m+1));
            end
            off_chromo(1,j)=off_chromo(1,j)+beta2*(x_max(j)-x_min(j));
            if off_chromo(1,j)>x_max(j)||off_chromo(1,j)<x_min(j)
                off_chromo(1,j)=x_min(j)+rand(1)*(x_max(j)-x_min(j));
            end
        end
    end
    if rand(1)<pm
        for j=1:length(pop_chromo(1,:))
            r=rand(1);
            if r<=0.5
                deta1=(off_chromo(2,j)-x_min(j))/(x_max(j)-x_min(j));
                beta2=(2*r+(1-2*r)*(1-deta1)^(yita_m+1))^(1/(yita_m+1))-1;
            else
                deta2=(x_max(j)-off_chromo(2,j))/(x_max(j)-x_min(j));
                beta2=1-(2-2*r+(2*r-1)*(1-deta2)^(yita_m+1))^(1/(yita_m+1));
            end
            off_chromo(2,j)=off_chromo(2,j)+beta2*(x_max(j)-x_min(j));
            %% 越界处理，该步骤可以更加严苛
            if off_chromo(2,j)>x_max(j)||off_chromo(2,j)<x_min(j)
                off_chromo(2,j)=x_min(j)+rand(1)*(x_max(j)-x_min(j));
            end
        end
    end
    % 对两个后代分别进行基因变异操作
    off(2*i-1:2*i,:)=off_chromo;
end
if p_eq
    new_pop=deal_eq5(off,Aeq,x_min);
else
    new_pop=off;
end

5.小结

完成上述部分的代码之后，就可以对遗传算法进行构建了，按照本文文初流程图搭建即可。种群选择这一步骤虽然在流程图中出现两次，但其实本质上一样，代码大差不差，可以用同一个函数进行选择操作。

另外，循坏外的第一次选择操作很明显可以加入到循环中，但是若加入循环，则在最后一次迭代中，种群交叉变异合并之后，并未进行筛选，很有可能在种群中存在劣解。故考虑按照次流程进行遗传算法的设计。

三、算法实例

本文共展示三个实例，题目来源于《最优化方法（第二版）》以及文末参考文献中的例题。

1.例一

目标函数为3*x1+6*x2

要求求解目标函数的最大值。

问题的约束条件为：
3*x1+x2>=9
x1+2*x2>=6
7*x1+5*x2<=35
x1,x2>=0

该问题为有不等式约束的线性优化问题，遗传算法求解该问题主函数代码如下所示：

clear
%% 初始化参数
pop=300;
% 种群大小，建议是偶数，后面有个地方默认偶数进行计算
gen=100;
% 遗传代数
x=zeros(2,gen);
% 记录轨迹
pc=0.9;
% 交叉概率
pm=0.1;
% 变异概率，建议小于0.1.
yita_c=10;
% 交叉中的分布参数
yita_m=20;
% 变异中的分布参数
x_min=[0,0];
% 自变量的最小值，无下界时，可用一个极小值替代
x_max=[5,7];
% 自变量的最大值，无上界时，可用一个极大值替代

long_chromo=[9,10];
temp="浮点数编码";

chromo=Initialize_pop(pop,2,x_min,x_max,temp,long_chromo);
% 种群初始化
chromo_value=object_fun(chromo,1,temp,long_chromo,x_min,x_max);
% 计算各个个体适应度
aary=zeros(1,gen);

% 迭代求解
for i=1:gen
    % 竞技选择法进行个体选取
    father=choose_pop(chromo_value,10);
    % 剔除适应度
    father=father(:,1:end-1);
    % 交叉变异结束
    son=cross_mutation(father,pc,pm,yita_c,yita_m,x_min,x_max,temp,long_chromo);
    % 父代与子代合并
    new_temp=[father;son];
    % 选择
    temp_value=object_fun(new_temp,i,temp,long_chromo,x_min,x_max);
    temp_value=choose_pop(temp_value,20);
    temp_value=sort(temp_value,3);
    chromo_value=temp_value(1:pop,:);
    aary(i)=mean(chromo_value(:,end));
    x(:,i)=chromo_value(1,1:2)';
    % 进度提示
    if mod(i,10)==0
        fprintf('%d gen has completed!\n',i);
    end
end

plot(aary,'o-')

求解结果如下图所示：

正确答案为35.25，算法计算结果为35.2455相对误差限很小，能够接受。除此之外，我们很明显能够发现，遗传算法收敛速度快的特点。

上述主函数代码中各个步骤的函数均在上文出现，唯有计算适应度的函数object_fun未出现代码，该函数代码如下：

function chromo_value=object_fun(chromo,w,temp,long_chromo,x_min,x_max)
% 目标函数为3*x1+6*x2的最大值
% 约束条件为
% x1+2*x2>=6
% 3*x1+x2>=9
% 7*x1+5*x2<=35
% x1,x2>=0
if temp=="浮点数编码"
    % 未加入约束条件
    value_first=chromo(:,1).*3+chromo(:,2).*6;
    % 采用罚函数的方法
    gama1=zeros(length(chromo(:,1)),1);
    gama2=zeros(length(chromo(:,1)),1);
    gama3=zeros(length(chromo(:,1)),1);
    for i=1:length(chromo(:,1))
        % 第一个约束条件
        if chromo(i,1)+chromo(i,2)*2>=6
            gama1(i)=0;
        else
            %gama1(i)=-w*(chromo(i,1)+chromo(i,2)*2-6)^2;
            gama1(i)=-inf;
        end
        % 第二个约束条件
        if 3*chromo(i,1)+chromo(i,2)>=9
            gama2(i)=0;
        else
            %gama2(i)=-w*(3*chromo(i,1)+chromo(i,2)-9)^2;
            gama2(i)=-inf;
        end
        % 第三个约束条件
        if 7*chromo(i,1)+5*chromo(i,2)<=35
            gama3(i)=0;
        else
            %gama3(i)=-w*(7*chromo(i,1)+5*chromo(i,2)-35)^2;
            gama3(i)=-inf;
        end
        % 第四个自适应满足
    end
    % 计算最终适应度
    value=value_first+w.*gama1+w.*gama2+w.*gama3;
    chromo_value=[chromo,value];
elseif temp=="二进制编码"
    %%
    x_num=length(long_chromo);
    pop=length(chromo(:,1));
    x_10=zeros(pop,x_num);
    x_2=x_10;
    long_chromo1=[0,long_chromo];
    for i=1:x_num
        t=chromo(:,long_chromo1(i)+1:sum(long_chromo1(1:i))+long_chromo1(i+1));
        idx=0:long_chromo(i)-1;
        idx=2.^idx;
        idx=sort(idx,'descend');
        x_2(:,i)=t*idx';
    end
    idx=2.^long_chromo-1;
    x_long=x_max-x_min;
    r=x_long./idx;
    x_10=x_2.*r;
    % 未加入约束条件
    value_first=x_10(:,1).*3+x_10(:,2).*6;
    % 采用罚函数的方法
    gama1=zeros(length(x_10(:,1)),1);
    gama2=zeros(length(x_10(:,1)),1);
    gama3=zeros(length(x_10(:,1)),1);
    for i=1:length(x_10(:,1))
        % 第一个约束条件
        if x_10(i,1)+x_10(i,2)*2>=6
            gama1(i)=0;
        else
            %gama1(i)=-w*(x_10(i,1)+x_10(i,2)*2-6)^2;
            gama1(i)=-inf;
        end
        % 第二个约束条件
        if 3*x_10(i,1)+x_10(i,2)>=9
            gama2(i)=0;
        else
            %gama2(i)=-w*(3*x_10(i,1)+x_10(i,2)-9)^2;
            gama2(i)=-inf;
        end
        % 第三个约束条件
        if 7*x_10(i,1)+5*x_10(i,2)<=35
            gama3(i)=0;
        else
            %gama3(i)=-w*(7*x_10(i,1)+5*x_10(i,2)-35)^2;
            gama3(i)=-inf;
        end
        % 第四个自适应满足
    end
    % 计算最终适应度
    value=value_first+w.*gama1+w.*gama2+w.*gama3;
    chromo_value=[chromo,value];
else
    fprintf("wrong in object_fun\n");
end

该部分代码仅作参考，不同的题目，该函数的代码一定是不同的，均需要针对问题进行修改，本文后续例子的object_fun函数的代码与此代码也是存在不小的差异的。

但本文在第四部分，将遗传算法进行了一个函数封装，该函数能够应对大多数的连续型优化问题，离散型优化问题的封装函数正在跟进，关注作者，等待下一版本的函数即可。

2.例二

题目如下图所示：

遗传算法求解该问题的主函数代码如下所示：

% 目标函数为最大值
% x1^2*x2*x3^2/(2*x1^3*x3^2+3*x1^2*x2^2+2*x2^2*x3^2+x1^2*x2^2*x3^2)
% 约束条件为
% x1^2+x2^2+x3^2>=1
% x1^2+x2^2+x3^2<=4
% x1,x2,x3>=0

clear
tic
%% 初始化参数
pop=300;
% 种群大小，建议是偶数，后面有个地方默认偶数进行计算
gen=100;
% 遗传代数
arry=zeros(2,gen);
% 记录轨迹
pc=1;
% 交叉概率
pm=0.1;
% 变异概率，建议小于0.1.
yita_c=10;
% 交叉中的分布参数
yita_m=50;
% 变异中的分布参数
x_min=[0,0,0];
% 自变量的最小值，无下界时，可用一个极小值替代
x_max=[3,3,3];
% 自变量的最大值，无上界时，可用一个极大值替代

%temp="二进制编码";
long_chromo=[9,10,9];
temp="浮点数编码";
fprintf("此次运行采用%s\n",temp)
chromo=Initialize_pop(pop,3,x_min,x_max,temp,long_chromo);
chromo_value=object_fun(chromo,1,temp,long_chromo,x_min,x_max);
% 计算各个个体适应度
arry(1,:)=1:gen;

% 迭代求解
for i=1:gen
    % 竞技选择法进行个体选取
    father=choose_pop(chromo_value,10);
    % 剔除适应度
    father=father(:,1:end-1);
    % 交叉变异结束
    son=cross_mutation(father,pc,pm,yita_c,yita_m,x_min,x_max,temp,long_chromo);
    % 父代与子代合并
    new_temp=[father;son];
    % 精英保留策略
    temp_value=object_fun(new_temp,i,temp,long_chromo,x_min,x_max);
    temp_value=choose_pop(temp_value,20);
    temp_value=sort(temp_value,3);
    chromo_value=temp_value(1:pop,:);
    arry(2,i)=mean(chromo_value(:,end));
    % 进度提示
    if mod(i,10)==0
        fprintf('%d gen has completed!\n',i);
    end
end
toc;
plot(arry(1,:),arry(2,:),'o-','LineWidth',1.5)

求解结果如下图所示：

按照书中给出的遗传算法计算结果为0.1537，显然，本文的遗传算法求解结果更优，约为0.1751，且完全符合约束条件。

3.例三

题目如下所示：

题目为包含等式、不等式线性约束的优化问题，该问题的遗传算法主函数代码如下：

%% 初始化参数
pop=50;
% 种群大小，建议是偶数，后面有个地方默认偶数进行计算
gen=200;
% 遗传代数
arry=zeros(2,gen);
% 记录轨迹
pc=1;
% 交叉概率
pm=0.3;
% 变异概率，建议小于0.1.
yita_c=2;
% 交叉中的分布参数
yita_m=10;
% 变异中的分布参数
x_min=[0,0];
% 自变量的最小值，无下界时，可用一个极小值替代
x_max=[5,5];
% 自变量的最大值，无上界时，可用一个极大值替代

long_chromo=[7,7];
temp="浮点数编码";
fprintf("此次运行采用%s\n",temp)

type='min';
fprintf("求解目标函数的%s值\n",type)

type_0="强初始化";
fprintf("初始种群生成采用：%s\n",type_0)

chromo=Initialize_pop(pop,length(x_min),x_min,x_max,temp,long_chromo,type_0);
chromo_value=object_fun(chromo,1,temp,long_chromo,x_min,x_max);
% 计算各个个体适应度
arry(1,:)=1:gen;

% 迭代求解
for i=1:gen
    % 竞技选择法进行个体选取
    father=choose_pop(chromo_value,10,type);
    % 剔除适应度
    father=father(:,1:end-1);
    % 交叉变异结束
    son=cross_mutation(father,pc,pm,yita_c,yita_m,x_min,x_max,temp,long_chromo);
    % 父代与子代合并
    new_temp=[father;son];
    % 精英保留策略
    temp_value=object_fun(new_temp,i,temp,long_chromo,x_min,x_max);
    temp_value=choose_pop(temp_value,20,type);
    if strcmp(type,'min')
        temp_value=sort(temp_value,3,"ascend");
    elseif strcmp(type,'max')
        temp_value=sort(temp_value,3,"descend");
    else
        fprintf("type is wrong")
    end
    chromo_value=temp_value(1:pop,:);
    arry(2,i)=mean(chromo_value(:,end));
    % 进度提示
    if mod(i,10)==0
        fprintf('%d gen has completed!\n',i);
    end
end
toc;
plot(arry(1,:),arry(2,:),'o-','LineWidth',1.5)

%% 问题结果输出
x=chromo_value(1,:);
x1=x(1);
x2=x(2);
x3=2-x2/2-x1;
x4=1-x2*5/2-2*x1;
fprintf("目标函数最终值为：%f\n",x(end))
for i=1:4
    fprintf("x%d的值为：%f\n",i,eval(['x',num2str(i)]))
end

该题正确答案为2.2，代码运行结果如下图所示：

很明显，求解效果很好，但是受限于等式约束，若采用与前两例题相同的遗传算法，则很容易得到无解或者很大的一个解，显然不符合要求。其主要原因在于，等式约束本质上是将可行解的空间降维，也就是可行解空间的维度与遗传算法的搜索空间维度不一致（小于遗传算法的搜索空间维度），导致普通的初始化难以搜索到可行解。

故在该题中，引入强初始化，逼迫初始种群必须全为符合各个约束条件的个体，但也因为这个原因，种群初始化耗费大量时间，但也因为强初始化，遗传算法的求解结果很好。最终用时90秒左右。

4.小结

在这一章节，本文举出三个例题，进行遗传算法的实战演示，针对问题的求解结果均表现出不错的结果，下面将简述一下本文是如何针对题目对遗传算法进行修改：

（1）目标函数在object_fun函数中变量value_first进行计算，读者可以仿照其格式进行修改；
（2）对不等式约束条件，本文采用的是罚函数法与内点法两种方法进行处理，展示代码中采用的是内点法，代码注释区域有注释错误——注释区域标注的是罚函数法，但是用的是内点法。
（3）对等式线性约束条件，本文上述代码采用的是第四篇参考文献中的处理方法，同时在种群初始化、目标函数等多个地方进行了或多或少的修改。
（4）对无约束的优化问题，只用删去代码中关于约束的部分即可。

四、算法函数封装

考虑到上述代码针对不同类型的问题，均需做出不小的改动，操作实属麻烦，且容易出现各种奇怪的报错，故本文在这一部分，将上述代码进行整合重构，将遗传算法整体封装成一个函数，并将上文中提到的两种编码方式、两种初始化方式均内嵌入函数中，使用者仅需调整即可。

1.示例一

问题同本文第三章的例三，使用封装函数的代码如下：

object="x1-2*x2+3";
% 字符串表示目标函数
Aeq=[0,0];
Beq=0;
% 等式约束
A=[-2,-8;
    4,5;
    2,1];
B=[-3;2;4];
% 不等式约束，必须设定为A<=B或者A<B
temp="浮点数编码";
%temp="二进制编码";
type='min';
%type='max';
type_0="强初始化";
%type_0="弱初始化";

%% 算法参数
pop=50;
% 种群大小，建议是偶数，后面有个地方默认偶数进行计算
gen=200;
% 遗传代数
pc=1;
% 交叉概率
pm=0.3;
% 变异概率，建议小于0.1.
yita_c=2;
% 交叉中的分布参数
yita_m=10;
% 变异中的分布参数
x_min=[0,0];
% 自变量的最小值，无下界时，可用一个极小值替代
x_max=[3,3];
% 自变量的最大值，无上界时，可用一个极大值替代

[chromo_best,best_value,convergence_gen]=GA_main_fun(object,Aeq,Beq,A,B,pop,gen,pc,pm,yita_c,yita_m,x_min,x_max,temp,type,type_0);

等式约束与不等式约束仅需输入系数矩阵与常数矩阵即可，遗传算法的算子选择通过消除注释与注释即可实现，运行结果如下所示:

相较于前文的代码，这个遗传算法的封装函数使用简单，效率也提升不少，但是聪明的小伙伴很容易就会发现，它仅支持线性的约束条件。但其并未解决线性等式约束的问题，仍旧是采用第四篇参考文献中的方法，使用者需提前对目标函数、约束条件进行化简，才可使用。

2.示例二

题目见下图：

封装函数代码为：

clear 
tic

object="-3.*x1-x2-4.*x3";
% 字符串表示目标函数

Aeq=["6.*x1+3.*x2+5.*x3+x4-45";
     "3.*x1+4.*x2+5.*x3+x5-30"];
Beq=zeros(length(Aeq),1);
% 等式约束，化简为Aeq=0

%A="x1-5.*x2+x3+x4";
A=[];
B=zeros(length(A),1);
% 不等式约束，必须设定为A<=B或者A<B

temp="浮点数编码";
%temp="二进制编码";
type='min';
%type='max';
%type_0="强初始化";
type_0="弱初始化";

%% 算法参数
pop=100;
% 种群大小，建议是偶数，后面有个地方默认偶数进行计算
gen=200;
% 遗传代数
pc=0.9;
% 交叉概率
pm=0.1;
% 变异概率，建议小于0.1.
yita_c=10;
% 交叉中的分布参数
yita_m=20;
% 变异中的分布参数
x_min=[0,0,0,0,0];
% 自变量的最小值，无下界时，可用一个极小值替代
x_max=[6,2,6,2,2];
% 自变量的最大值，无上界时，可用一个极大值替代


x_min=x_min-0.00001;
x_max=x_max+0.00001;

[chromo_best,best_value,convergence_gen]=GA_main_fun(object,Aeq,Beq,A,B,pop,gen,pc,pm,yita_c,yita_m,x_min,x_max,temp,type,type_0);
toc

运行结果如下所示：

正确答案为-27，很明显，该代码在处理更为复杂的等式约束条件时，代码运行时间增长幅度很大，好在运行效果十分理想，相较于上个示例中封装好的函数，本示例中函数支持以字符串的形式进行输入约束条件，且使用者无需手动化简等式约束条件，也能针对非线性等式约束条件进行求解，使用起来更为方便。

3.示例三

题目同上一示例，但在初始化、选择算子、处理等式约束等多个部分进行进一步的优化，求解结果如下所示（使用代码与上一示例相同，仅将迭代次数从200改为100）：

初始化部分采用了普通初始化与强初始化相结合的办法；选择算子中的随机抽样规模采用了自适应参数，不再是之前固定的20或者10；处理等式约束部分采用了更高效率的求解办法。在示例二中，迭代次数为100时运行时长约为200秒左右，此时仅需55秒左右（多次运行的结果）。

需要封装函数的请私信获取代码压缩包

五、参考文献

[1] 施光燕,钱伟懿,庞丽萍.最优化方法(第二版)[M].北京:高等教育出版社,2007:141-145.
[2]闫泽远.多目标进化算法的算子选择研究[D].山东师范大学,2022.
[3]王之仓,李和成.采用多项式变异策略和分解方法的多目标进化算法[J].微电子学与计算机,2021,38(01):95-100.
[4]胡宽,常新龙,宋笔锋,等.求解含等式约束优化问题的遗传算法[J].上海交通大学学报,2011,45(07):966-969+974.
[5]刘伟,蔡前凤,王振友.一种新的遗传算法求解有等式约束的优化问题[J].计算机工程与设计,2007,(13):3184-3185+3194.