Matlab应用笔记--遗传算法

注：本篇随笔依据《Matlab在数学建模上的应用》中第4章介绍来写，主要介绍简单遗传算法的思想及其Matlab实现
（博客以及Matlab小白，若有不当欢迎指出）

遗传算法（GA）简介

模拟达尔文生物进化论的自然选择和孟德尔遗传学机理的生物进化过程的计算模型，一种通过模拟自然进化过程搜索最优解的方法。
遗传算法本质是启发式随机搜索算法，通过遗传算法得到的解多是全局最优解。
简单的遗传算法注重模拟竞争关系，而自然界还有其它作用（协作、寄生关系等），当问题的优化异常复杂难以处理时，传统的简单遗传算法无能为力，这时引入协同关系的“协同进化遗传算法”就发挥作用了。

遗传算法的实现

（1）编码

遗传算法编码主要有浮点编码和二进制编码两种，这里只介绍二进制编码。

设某一参数的取值范围是$\left ( L,U \right )$，精度是$p$，则有效取值长度为$length=\frac{U-L}{p}$
现在，要用一个二进制数串能表示所有可能数值，则该数串的长度$k$应该满足
$2^{k-1}-1< length\leq 2^{k}-1$ （数串0代表L，$2^{k}-1$代表U）
（要保证二进制数串表示的数精度不得低于$p$，即$\delta =\frac{length}{2^{k}-1}\leq 1$）
若待编码的数值为$m$，则$m$与$L$间的长度为$length'=\frac{m-L}{p}$
对应的“二进制编码大小”为“$m$与$L$间的偏移长度”除以“二进制数的精度”，
即$val=length'\div \delta =length'\div \frac{length}{2^{k}-1}=length'\times \frac{2^{k}-1}{length}$
最后对直接对val进行转换即可（普通的十进制转二进制）

例如：取值范围为$\left ( 3,34 \right )$，精度是0.01，则$length=\frac{34-3}{0.01}=3100$
令$k=12$，则满足$2^{12-1}-1< 3100\leq 2^{12}-1$
若待编码的数为$17$，则$17$与$3$间的偏移长度为$length'=\frac{17-3}{0.01}=1400$
所以$val=1400\div \frac{3100}{2^{12}-1}=1400\times \frac{2^{12}-1}{3100}\approx 1849$
直接对1849转成二进制数串：$0111\; 0011\; 1001$

（2）解码

解码是编码的逆过程。

设某一参数的取值范围是$\left ( L,U \right )$，精度是$p$，则有效取值长度为$length=\frac{U-L}{p}$
二进制数串的长度直接读出为$k$
二进制数的精度是$\delta =\frac{length}{2^{k}-1}$
若待解码的数串为$b_{k}b_{k-1}b_{k-2}...b_{3}b_{2}b_{1}$，先将其直接转为二进制编码大小$val$
则“待求十进制数$m$与$L$间的偏移长度”为“二进制编码大小”乘以“二进制数的精度”
即$length'=val\times \delta =val\times \frac{length}{2^{k}-1}$
最后$L$加上“偏移长度乘以十进制精度”即为待求十进制数：$m=L+lenght'\times p$

例如：取值范围为$\left ( 3,34 \right )$，精度是0.01，则$length=\frac{34-3}{0.01}=3100$
待解码数串为$0111\; 0011\; 1001$
长度$k=12$，二进制数的精度是$\delta =\frac{3100}{2^{12}-1}$
转成十进制数：$val=1849$
偏移长度：$length'=1849\times \frac{3100}{2^{12}-1}\approx 1400$
待求十进制数：$m=3+1400\times 0.01=17$

（3）个体适应度评估

根据一定的筛选机制算出染色体的适应度，以便后继的复制和选择操作。
通常求目标函数最大值的问题可以直接把目标函数作为检测个体适应度大小的函数。

先将染色体串解码为十进制真值$m$
根据真值$m$评价目标函数$f(m)$
将目标函数值转为适应度$eval(U)=G(f(m))$。
（对于极大值问题，目标函数值可直接作为适应度）

例如：给定染色体串$0111\; 0011\; 1001$
解码后：$m=17$
设目标函数：$f(x)=e^{\frac{x}{2}}sin(x)$，代入得$f(17)\approx 1437$
若对于求极大值问题，令$G(x)=x$则$eval(U)=G((f(17))=G(1437)=1437$
即染色体$0111\; 0011\; 1001$在当前环境的适应度为1437

（4）染色体复制

复制运算是根据个体适应度大小决定下代遗传的可能性。
适应度越高（低），复制的几率越大（小），遗传基因就会在种群中扩散（淘汰）。

设种群中个体总数为$N$，个体$i$的适应度为$f_{i}$，则个体i被选取的几率为
$P_{i}=\frac{f_{i}}{\sum_{k=1}^{N}f_{k}}$
个体复制的几率决定后，再产生$[0,1]$区间的均匀随机数决定哪些个体参加复制

例如：种群中个体总数为100，第29个个体的适应度为$f_{29}=5140$
所有个体的适应度之和为$\sum_{k=1}^{100}f_{k}=114514$
则第29个个体被选取的几率为$P_{29}=\frac{5140}{114514}\approx4.5％$
将所有个体的复制几率按序号依次填充到区间$[0,1]$中
假设第29个个体的覆盖区间是$％33\sim％37.5$
再产生$[0,1]$区间的均匀随机数，若该随机数在$％33\sim％37.5$内，则第29号个体被选中进行复制。

（4）染色体交叉（配）

交叉运算是使用单点或多点进行交叉的算子。
首先用随机数产生一个或多个交配点位置，然后两个个体在交配点位置互换部分基因码，形成两个子个体。

设染色体总量为$N$，交配概率为$P_{c}$，则参与交配的染色体数量大概为$n=[N\times P_{c}]$（取整）
产生在$[0,1]$区间内的$N$个随机数$R_{1}、R_{2}、...、R_{N}$，挑选出最低的$n$个随机数作为交配对象
（或者相邻序号为一对，从$1$号开始到$N$号依次产生$\frac{N}{2}$对，每一对内产生一随机数$R$，若$R<P_{c}$则交配，否则不交配）
交配时以两个染色体位单位进行，随机产生一个$[0,k-1]$区间内的数$t$（k是染色体长度）
然后将两个染色体的第t个位点右边的位段（$[t+1,k-1]$部分）进行互换，得到两个新的子辈染色体。

例如：染色体总量为$100$，交配概率为$25％$，则参与交配的染色体数量大概为$n=25$
以相邻序号为一对，从$1$号开始到$100$号依次产生$50$对，每一对内产生一随机数
假设第7对（第13、14号染色体）产生的随机数$0.17<0.25$，它们进行交配
假设染色体位长为12，再在$[0,11]$区间内产生一随机数$t=5$为交配位点
假设原来的13、14号染色体分别为
$U_{13}=0110\; 00{\color{Red} {11\; 0111}}$
$U_{14}=1010\; 11{\color{Blue} {01\; 1001}}$
将它们的$[6,11]$位段进行互换，得子代染色体分别为
$U'_{13}=0110\; 00{\color{Blue} {01\; 1001}}$
$U'_{14}=1010\; 11{\color{Red} {11\; 0111}}$

（5）染色体变异（基因突变）

依照突变运算规则，种群内所有的基因（每个二进制位）都有概率进行突变（取反）。

设染色体总量为$N$，每个染色体的长度为$k$，突变概率为$P_{m}$，则每一代突变的基因数大概为$n=[N\times k\times P_{m}]$（取整）
产生$N\times k$个位于$[0,1]$区间内的随机数，随机数从$1$开始到$N\times k$依次编号
若某号$num$随机数小于突变概率$P_{m}$，则对应染色体内的某个二进制位取反
（$\left \lfloor \frac{num}{k}\right \rfloor$号染色体上的$num％k$号基因）

例如：设染色体总量为$10$，每个染色体的长度为$12$，突变概率为$P_{m}=0.01$，则每一代突变的基因数大概为$n=1$（取整）
产生$120$个位于$[0,1]$区间内的随机数，随机数从$1$开始到$120$依次编号
假设第75号随机数$num=0.005<0.01$
对应第$\left \lfloor \frac{75}{12}\right \rfloor=6$染色体上的$75％12=3$号基因位进行突变，即
$U_{6}=10{\color{Red} 1}1\; 0110\; 0001$ （原6号染色体）
$U'_{6}=10{\color{Red} 0}1\; 0110\; 0001$ （基因突变后6号染色体）

（6）染色体倒位

一般简单的问题不用考虑染色体倒位，当问题比较复杂时则可能需要倒位运算。
这里只简单介绍倒位的效果，具体实现不细展开（可自行实现）。

$U=110{\color{Red} {1\; 0100\; 1011\; 00}}01$（倒位前）
$U=110{\color{Red} {0\; 0110\; 1001\; 01}}01$（选择$[4,14]$位段进行180°翻转）

遗传算法的Matlab实现

（1）Matlab遗传算法工具箱

本篇随笔只基于Matlab7.0 Version版本的工具箱讨论，其他版本的工具箱可自行尝试或找其它方法。
Matlab的遗传算法工具箱（GA Toolbox）对常用的遗传运算进行集成，用户能方便调用，但是过度的封装也让该工具箱缺乏灵活性，某些问题最好还是编写算法源程序去解决。
Matlab7.0 Version的遗传算法工具箱有自带的遗传算法与直接搜索工具箱。

使用方法：直接在命令窗口中输入optimtool，在弹出的界面上填写相关参数和函数句柄或使用M文件进行程序设计。（这里只细说程序设计部分，工具箱的其它具体使用部分略）
遗传算法与直接搜索工具箱有ga和gaoptimset两个核心函数。

函数ga的一般语法格式

[x,fval,reason] = ga(@fitnessfun,nvars,options)

x为遗传进化后自变量最佳染色体返回值
fval为最佳染色体的适应度
reason为算法停止的原因
@fitnessfun为适应度句柄函数
nvars为目标函数自变量的个数
options为算法的属性设置（通过函数gaoptimset赋予）

函数gaoptimset的语法格式

options = gaoptimset('PropertyName1',PropertyValue1',
'PropertyName2',PropertyValue2','PropertyName3',PropertyValue3',...,...)

通过不同的“属性名”和“属性值”设置遗传算法的参数，返回句柄options用于函数ga中

（函数gaoptimset属性）

一般程序设计

function f = fitnessfun(x)
%适应性函数设计
 
if(...)      %符合约束条件
      f = f(x);      %目标函数的表达式     
else         %不符合约束条件
      f = ...;      %相应处理（比如要求最大值时可以令f = -inf）
      end

%主函数
clc;
clear;
close all;
 
options = gaoptimset(...);      %设置好属性
[x,fval,reason] = ga(@fitnessfun,nvars,options)  %根据设置好的属性和目标函数求解

若采用接力进化，则主函数还可以这样写

%主函数
clc;
clear;
close all;
 
%第一轮遗传运算
options = gaoptimset(...);
[x,fval,reason,output,finnal_pop] = ga(@fitnessfun,nvars,options);
 
%第二轮遗传运算，初始化种群为第一趟结束时的种群
options = gaoptimset('InitialPopulation',finnal_pop,...);
[x,fval,reason,output,finnal_pop2] = ga(@fitnessfun,nvars,finnal_pop);
%若需要，可以往后继续接力

（2）Matlab遗传算法的源程序模版

工具箱并非万能，很多时候工具箱不够灵活，需要我们自己写出源程序以针对复杂的多约束非线性规划问题。

一般遗传算法的设计（伪代码）

t = 0;      %遗传代数
初始化P(t);      %初始化种群
whlie(不满足停止准则）
begin
      t = t+1;
      从P(t-1)中选择P(t);      %选择
      重组P(t);      %交叉和变异
      计算P(t)的适应值；
end

一个简单遗传算法的详细程序设计模版

主函数

%主函数
clc;
clear;
close all;
 
%全局变量
global BitLength      %染色体位长
global boundsbegin      %变量取值范围中的下限
global boundsend        %变量取值范围中的上限
 
bounds = [-2 2];    %一维自变量取值范围S
precision = 0.0001; %运算精度
boundsbegin = bounds(:,1);
boundsend = bounds(:,2);
BitLength = ceil(log2((boundsend - boundsbegin)'./precision));
 
popsize = 50;           %种群大小
Generationmax = 12;     %遗传代数
pcrossover = 0.90;      %交配概率
pmutation = 0.09;       %变异概率
 
%产生初始种群
population = round(rand(popsize,BitLength));
%计算适应度，返回适应度Fitvalue和累积概率cumsump
[Fitvalue,cumsump] = fitnessfun(population);
 
%开始迭代模拟遗传
Generation = 1;
xmax(1:Generationmax) = 0;
ymax(1:Generationmax) = 0;
ymean(1:Generationmax) = 0;
while Generation <= Generationmax
    for j = 1:2:popsize
        %复制操作略，可添加
        seln = selection(cumsump);%选择
        scro = crossover(population,seln,pcrossover);%交配
        %变异
        smnew(j,:) = mutation(scro(1,:),pmutation);
        smnew(j+1,:) = mutation(scro(2,:),pmutation);
    end
    %产生新种群
    population = smnew;
    %计算新种群的适应度
    [Fitvalue,cumsump] = fitnessfun(population);
    %记录当前代最好的适应度和平均适应度
    [fmax,nmax] = max(Fitvalue);
    ymax(Generation) = fmax;
    ymean(Generation) = mean(Fitvalue);
    %记录当前代的最佳染色体个体
    x = transform2to10(population(nmax,:));%译码
    xx = boundsbegin + ...      %转到区间范围内对应的值
        x*(boundsend-boundsbegin)/(power(2,BitLength)-1);
    xmax(Generation) = xx;
    Generation = Generation + 1;
end
 
%记录所有代中最好的适应度及其对应的染色体
Generation = Generation - 1;
[fmax,nmax] = max((ymax'));
Bestpopulation = xmax(nmax)
Besttargetfunvalue = targetfun(Bestpopulation)
 
%绘制曲线
%若发现平均适应度与最大适应度在曲线上有趋同的形态，则表明算法收敛进行顺利
%最大适应度个体连续若干代都没有发生进化表明种群已经成熟
figure(1);
hand1 = plot(1:Generation,ymax);
set(hand1,'linestyle','-','linewidth',1.8,'marker','*',...
    'markersize',6)
hold on;
hand2 = plot(1:Generation,ymean);
set(hand2,'color','r','linestyle','-','linewidth',1.8,...
'marker','h','markersize',6)
xlabel('进化代数');ylabel('最大/平均适应度');xlim([1 Generationmax]);
legend('最大适应度','平均适应度');
box off;
hold off;

计算适应度的函数

function [Fitvalue,cumsump] = fitnessfun(population)
%计算适应度的函数
 
global BitLength
global boundsbegin
global boundsend
 
popsize = size(population,1);   %求个体数
Fitvalue(1:popsize) = 0;    %初始化
cumsump(1:popsize) = 0;     %初始化
 
%计算适应值
for i = 1:popsize
    x = transform2to10(population(i,:));    %译码
    xx = boundsbegin + ...      %转到区间范围内对应的值
        x*(boundsend-boundsbegin)/(power(2,BitLength)-1);
    Fitvalue(i) = targetfun(xx);    %计算适应度
end
Fitvalue = Fitvalue' + 230;     %保证适应值为正数
 
%计算选择概率
fsum = sum(Fitvalue);
Pperpopulation = Fitvalue/fsum;
%计算累积概率
cumsump(1) = Pperpopulation(1);
for i = 2:popsize
    cumsump(i) = cumsump(i-1)+Pperpopulation(i);
end
cumsump = cumsump';

新种群交叉操作

function scro = crossover(population,seln,pc)
%新种群交叉操作
 
global BitLength
pcc = IfCroIfMut(pc);   %根据交叉概率决定是否交叉
if pcc == 1
    chb = round(rand*(BitLength-2))+1;  %随机产生交叉位
    scro(1,:) = [population(seln(1),1:chb) population(seln(2),chb+1:BitLength)];
    scro(2,:) = [population(seln(2),1:chb) population(seln(1),chb+1:BitLength)];
else
    scro(1,:) = population(seln(1),:);
    scro(2,:) = population(seln(2),:);
end

新种群变异操作

function snnew = mutation(snew,pmutation)
%新种群变异操作
 
global BitLength
 
snnew = snew;
pmm = IfCroIfMut(pmutation);    %根据变异概率决定是否变异
if pmm == 1
    chb = round(rand*(BitLength-1))+1;  %随机产生交叉位
    snnew(chb) = abs(snew(chb)-1);
end

判断遗传运算是否需要进行交叉或变异

function pcc = IfCroIfMut(mutORcro)
%判断遗传运算是否需要进行交叉或变异
 
%设立长度为100的池，在里面抽选
test(1:100) = 0;
length = round(100*mutORcro);
test(1:length) = 1;
n = round(rand*99)+1;
pcc = test(n);

从种群中选中两个个体

function seln = selection(cumsump)
%从种群中选中两个个体
 
seln(1:2) = 0;
for i = 1:2
    r = rand;   %产生一个随机数
    j = 1;
    while cumsump(j)<r
        j = j+1;
    end
    seln(i) = j;    %选中个体的序号
end

将二进制数转换为十进制数

function x = transform2to10(Population)
%将二进制数转换为十进制数
 
global BitLength
 
x = Population(BitLength);
for i = 1:BitLength-1
    x = x+Population(BitLength-i)*power(2,i);
end

目标函数

function y = targetfun(x)
%目标函数
 
y = 200*exp(-0.05*x).*sin(x);

书上例题以及实机操作结果：

例1

结果

例2

遗传500代

遗传100代

遗传36代

遗传12代

原解决算法中的变异概率设置为0.09，若改为0.01，则在多代遗传的效果会更好

例3

获得的解不稳定
（某次跑的结果）

（另一次跑的结果）