遗传算法原理与 MATLAB 实现

一个具体的例子：求二元函数最大值

最近在学遗传算法，找了一道经典的例题把整个流程走一遍。

问题：求下列二元函数的最大值

个体编码

遗传算法的运算对象是字符串，所以必须先将 x1、x2 编码。这里用二进制整数来表示：x1、x2 都是 1~7 范围内的整数，分别用 3 位无符号二进制整数来表示。把 x1x2 拼在一起作为一个个体进行整体编码，例如基因型 101110 对应的表现型为 x=[5, 6]。表现型和基因型之间可以方便地互相转换。

初始群体的产生

遗传算法的本质是对群体进行进化，需要先准备一批初始个体。本例中初始群体数量取为 4，每个个体的基因型通过随机方式产生，例如：

011101, 101011, 011100, 111001

适应度计算

遗传算法以适应度作为评价个体优劣的指标，从而决定其遗传机会（淘汰还是继续遗传）。本例的优化目标是取非负数最大值，因此直接用函数值作为适应度——函数值越大，个体越好，适应度越高。

选择运算

选择运算将适应度较高的个体按照一定规则遗传到下一代群体中。适应度越高的个体获得更多遗传机会。本例采用与适应度成正比的概率来确定哪些个体进入下一代，具体步骤如下：

第一步，计算所有个体适应度的总和：

第二步，计算每个个体的相对适应度：

这就是每个个体被遗传到下一代的概率。

第三步，各概率值构成若干区间，所有概率之和为 1。

第四步，对每个个体产生一个 0~1 之间的随机数（本例共产生 4 个），根据随机数落在哪个概率区间来确定被选中的个体。

交叉运算

交叉运算是遗传算法中产生新个体的主要操作，以某一概率相互交换两个个体之间的部分染色体。本例采用单点交叉：

1. 先对群体进行随机配对（1-2 号、3-4 号）
2. 随机设置交叉点位置
3. 相互交换配对染色体之间的部分基因

可以看出，新产生的个体 111101 和 111011 的适应度相比原来两个个体的适应度都要高。

变异运算

变异运算是对个体的某一个或若干基因座上的基因值按小概率进行改变，也是产生新个体的一种方式。本例采用基本位变异：

1. 首先确定各个体的基因变异位置（下表中的数字表示变异点设置在该基因座处）
2. 按照某一概率将变异点的位取反

新一代群体 P(t+1)

对群体 P(t) 进行一轮选择、交叉、变异运算后，得到新一代群体 P(t+1)：

MATLAB 代码实现

以上流程用 Sheffield 大学遗传算法工具箱（Sheffield GA Toolbox）在 MATLAB 中实现。先自定义优化函数：

function z=test_shang(x,y)
z=x.^2+y.^2;
end

遗传算法主体代码：

NIND=40;        % 个体数目
NVAR=2;         % 变量维数为2
MAXGEN=100;     % 最大遗传代数
PRECI=20;       % 变量的二进制位数
GGAP=0.9;       % 代沟为0.9，舍掉10%适应度比较低的个体

trace=zeros(2,MAXGEN); % 算法性能跟踪，创建2*MAXGEN大小的零矩阵

Chrom=crtbp(NIND,PRECI*NVAR); % 初始化种群

%FieldD=[10;-1;2;1;0;1;1];
FieldD=[rep(length(Chrom)/2,[1,2]); rep([-1;2],[1,2]); rep([1;0;1;1],[1,2])];

% 建立区域描述器从左到右分别为
% FieldD=[len lb ub code scale lbin ubin]
% len:包含在Chrom在其中的每个子串的长度，1个变量则与变量的二进制位数相同
%     rep(length(Chrom)，2个则为rep(length(Chrom)/2
% lb,ub:分别指每个变量的下界和上界区间，为[-1,2]
% code：编码模式，取code=1，标准二进制编码
% scale：一般为算术刻度，scale=0
% lbin,ubin：是否包括边界，0则去掉边界，1则包括边界

gen=0; % 初始为0代

variable=bs2rv(Chrom,FieldD);  % 二进制字符串到表现型的转换
ObjV=test_shang(variable(:,1),variable(:,2)); % 初始化，ObjV为优化目标

while gen<MAXGEN
    FitnV=ranking(ObjV);                          % 分配适应度的值
    SelCh=select('sus',Chrom,FitnV,GGAP);         % 选择运算
    SelCh=recombin('xovsp',SelCh,0.7);            % 重组
    SelCh=mut(SelCh);                             % 变异
    variable=bs2rv(SelCh,FieldD);                 % 二进制转换为表现型
    ObjVSel=test_shang(variable(:,1),variable(:,2)); % 计算子代的目标函数值
    [Chrom, ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); % 重插入
    variable=bs2rv(Chrom,FieldD); % 初始种群转换为表现型
    % Chrom为初始种群，SelCh为过程中的值，为二进制
    gen=gen+1;
    [Y,I]=max(ObjV); hold on; % 找出群体里面的最大最小值
    trace(1,gen)=max(ObjV);
    trace(2,gen)=sum(ObjV)/length(ObjV);
end

variable=bs2rv(Chrom,FieldD);
hold on grid;

figure(1);
plot(variable(:,1),variable(:,2),'o');

figure(2);
plot(trace(1,:));
hold on;
plot(trace(2,:),'-.');grid
legend('解的变化','种群均值的变化')

trace 矩阵记录了每代的最优解和种群均值，最后可以画出收敛曲线。

Sheffield GA 工具箱函数速查

FieldD 结构

FieldD=[len  lb  ub  code  scale  lbin  ubin]

• len：包含在 Chrom 中每个子串的长度（PRECI）
• lb, ub：分别指每个变量的下界和上界
• code：编码方式，1 为标准的二进制编码，2 为格雷编码
• scale：刻度类型，0 为算术刻度，1 为对数刻度
• lbin, ubin：是否包含边界，0 去掉边界，1 包含边界

多维情形：

FieldD=[rep(len(Chrom),[1,维度]); rep([lb;ub],[1,维度]); rep([code;scale;lbin;ubin],[1,维度])];

种群初始化：crtbp、crtrp、crtbase

crtbp：创建二进制种群

[Chrom,Lind,BaseV]=crtbp(Nind, Lind*维度);  % 多维函数需乘以维度
[Chrom,Lind,BaseV]=crtbp(Nind, BaseV);
BaseV = crtbase([Nind,Lind],[a,b]);

Nind 为种群中个体的数量，Lind 指定个体或染色体的长度。

crtrp：创建实值原始种群

Chrom=crtrp(Nind, FieldDR);

FieldDR 为 2×Nvar 的矩阵，第一行为取值下界，第二行为取值上界；Nvar 向量的长度即为染色体的长度。

bs2rv：二进制串到实值的转换

Phen=bs2rv(Chrom, FieldD);

len 表示染色体的长度；lb, ub 分别为变量取值的上下界；code(i)=1 为标准的二进制编码，code(i)=0 为格雷编码；scale(i)=0 为算术刻度，scale(i)=1 为对数刻度；lbin, ubin 表示取值是否含有边界，取零去掉边界，取 1 含有边界。

适应度计算：ranking、scaling

ranking：基于排序的适应度分配

FitV=ranking(ObjV, RFun, SUBPOP);

• RFun(1)：线性排序标量在 [1, 2] 间，默认为 2；非线性排序在 [1, length(ObjV)-2]，指定选择压差
• RFun(2)：排序方法，0 为线性排序（默认），1 为非线性排序
• SUBPOP：ObjV 中子种群的数量，默认为 1

scaling：线性适应度计算

FitnV = scaling(ObjV, Smul);

注意：线性比率不适合目标函数返回负的适应度值的情形。

选择高级函数：select

SelCh=select(SEL_F, Chrom, FitnV);
SelCh=select(SEL_F, Chrom, FitnV, GGAP);
SelCh=select(SEL_F, Chrom, FitnV, GGAP, SUBPOP);

SEL_F 是一字符串，为低级选择函数名，如 rws 或 sus；GGAP 指出了代沟，默认为 1，也可大于 1，允许子代数多于父代的数量。

rws：轮盘赌选择

NewChrIx=rws(FitnV, Nsel);

使用轮盘赌选择从种群中选择 Nsel 个个体，NewChrIx 是被选个体的索引值。

sus：随机遍历抽样

NewChrIx=sus(FitnV, Nsel);

reins：重插入子群到种群（有代沟时必须使用）

Chrom=reins(Chrom, SelCh);
Chrom=reins(Chrom, SelCh, SUBPOP);
Chrom=reins(Chrom, SelCh, SUBPOP, InsOpt, ObjVch);
[Chrom, ObjVch]=reins(Chrom, SelCh, SUBPOP, InsOpt, ObjVch, ObjVSel);

• InsOpt(1)：用子代代替父代的选择方法，0 为均匀选择，1 为基于适应度的选择，默认为 0
• InsOpt(2)：每个子种群中重插入的子代个体比率，默认为 1
• ObjVch：Chrom 中个体的目标值，基于适应度的重插入时必需
• ObjVSel：SelCh 中个体的目标值，子代数量大于重插入种群子代数量时必需

交叉函数

重组准则：奇数行与它的下一个偶数行配对，如果矩阵 OldChrom 的行数是奇数，最后一个奇数行不交配并添加到 NewChrom 的最后一行。

recombin：重组个体高级函数

NewChrom=recombin(REC_F, Chrom);
NewChrom=recombin(REC_F, Chrom, RecOpt);
NewChrom=recombin(REC_F, Chrom, RecOpt, SUBPOP);

REC_F 是低级重组函数名字符串，例如 recdis, recint, reclin, xovdp, xovdprs, xovmp, xovsh, xovshrs, xovsp, xovsprs；RecOpt 为交叉概率。

常用低级交叉函数：

变异函数：mutate、mut、mutbga

mutate：高级变异函数

NewChrom=mutate(MUT_F, OldChrom, FieldDR, MutOpt, SUBPOP);

MUT_F 为低级变异函数名字符串（如 mut, mutbga）；SUBPOP 默认为 1，决定 OldChrom 中子种群的数量。

mut：离散变异算子

NewChrom=mut(OldChrom, Pm);
NewChrom=mut(OldChrom, Pm, BaseV);

Pm 为变异概率，默认为 0.7/Lind。

mutbga：实值种群的变异

NewChrom=mutbga(OldChrom, FieldDR, MutOpt);

• MutOpt(1)：变异概率，默认为 1/Nvar
• MutOpt(2)：在 [0, 1] 间压缩重组范围的标量，默认为 1（不压缩）

其他函数

rep：矩阵复制函数

MatOut=rep(MatIn, REPN);

REPN 为二值向量，REPN(1) 指明纵向复制次数，REPN(2) 指明水平复制次数。rep 是低级复制函数，通常不直接使用，被 GA 工具箱中许多函数调用。

migrate：在子种群间迁移个体

Chrom = migrate(Chrom, SUBPOP);
Chrom = migrate(Chrom, SUBPOP, MigOpt);
Chrom = migrate(Chrom, SUBPOP, MigOpt, ObjV);
[Chrom,ObjV] = migrate(Chrom, SUBPOP, MigOpt, ObjV);

MigOpt 是具有三个参数的可选向量：

• MigOpt(1)：迁移概率，默认 0.2
• MigOpt(2)：迁移选择方式，0 为均匀迁移（默认），1 为基于适应度的迁移
• MigOpt(3)：迁移种群结构，0 为完全网状结构（默认），1 为临近结构，2 为环状结构

ObjV 为所有个体对应目标函数值的列向量，当 MigOpt(2)=1 时为必选项。