遗传算法的MATLAB实现

0、前言

最近，大家听了专业老师关于遗传算法的相关介绍，我相信大家都是比较清楚的，在此我结合自己的理解以及Y老师的教授，归纳整理如下：

1、遗传算法介绍

遗传算法，顾名思义就是模拟达尔文进化论的自然选择和遗产学机理的生物进化构成的计算模型，是一种不断选择优良个体的算法。该算法的所涉及到的主要函数也就是包括自然界里“遗传”所对应的几种情况，自然界的遗传主要过程包括染色体的选择、交叉、变异等等，大自然在经过这些操作后，基本就可以保证以后的个体基本上是最优的，那么以后再继续这样下去，就可以一直最优了，当然也不排除保留了一些“不那么优秀”的个体的存在，这也就是遗传算法的一个弊端：容易陷入局部最优。

2、解决的问题

遗传算法很有名，自然能解决的问题就比较多了，在原理上不变的情况下，只要改变模型的应用环境和形式，基本上都可以。但是遗传算法主要还是解决优化类问题，尤其是那种不能直接解出来的很复杂的问题，而实际情况通常也是这样的：遗传算法也主要是被用作于求解复杂非线性拟合的系数、求取复杂函数的最大值等等

3、开始

每个算法都有自己的开始方式，遗传算法首先是选择个体了。我们知道一个种群中可能只有一个个体吗？不可能吧，肯定很多才对，这样相互结合的机会才多，产生的后代才会多种多样，才会有更好的优良基因，有利于种群的发展。那么算法也是如此，当然个体多少是个问题，一般来说20-100之间我觉得差不多了（当然了一般来说种群规模越大，所产生的不同个体就可能越多，最后的结果就可能越好，但是所带来的运算量也是十分庞大的，请慎重考虑）。那么个体究竟是什么呢？在编写代码过程中，可能就是一行或者一列矩阵了。其他情况下，个体就是所求问题的变量，这里我们假设个体数选100个。好了，现在有了100个个体组成的一个种群，那么这个种群应该怎么发展才能越来越好？说到这，我们想想，如何定义这个越来越好呢？这个应该有一个评价指标吧。我们甚至可以根据评价指标最后得分的大小来给他们排个名来决定哪些好哪些不好（不嫌麻烦的话可以利用什么综合评价啊什么的评价类算法进行处理也未尝不可）。把这个叫做对于个体的适应度，这应该算是算法的后半部分才对（其实都是遗传算法的几个函数中的一部分）。

4、编码

首先明白什么是编码？为什么要编码？如何编码？

什么是编码？其实编码就是把自变量（x）换一下形式，在这个形式下，更容易操作其他过程（比如交叉，变异什么的），通常情况下我们是进行二进制编码，也就是把十进制的自变量转变为二进制。我感觉这个编码和人体DNA，基因的排列很相似。我们高中所学的DNA编码也就是只有几个有限的“单元”不断“组合”成的双螺旋DNA链。

那么什么是二进制编码？很简单，就是1,0对应的来组合排列而已。比如：1100100010, 0011001001等等，这些都是位数长度为10的二进制编码。再想想1在计算机的二进制形式是什么？如果八位来表示的话，是不是就是0000 0001；8是不是就是0000 1000；以此类推。那么问题来了，我们需要多少位的二进制编码才能满足我们的需要呢？事实上，究竟是多少要视情况而定，一般20位左右感觉就可以了，虽然说越大越好，但是太大了消耗内存会导致运行速度很慢。

5、交叉与变异

先说变异，什么是变异？生物上基因或者染色体发生突变就叫变异。对应到我们编程上来说，变异就是一行矩阵中有“一段”或者“单个”发生了变异。

首先就讲最简单的变异，就是个体的变异。现在以10位长的编码来说，比如把x=3编码一下，假设为110 100，在变异操作时，假设第5位变异了（说一下变异就是一位或者多位1或0变成0或1），那么这个时候变成了110101，那么再反编码回去成x是多少呢？那肯定不是3了，那么把变化了之后的x带到适应度函数里面比较这两个x值对应的y值也就不同了，同时呢，可以通过这一步操作来判断该变异的“方向”：如果后面变异后的大一些就是说产生了好的变异，就可以在下一次个体选择的时候选择它了。

那么我们一起考虑一些这个事情，如果有很多很多x来一起变异会怎么样呢？肯定会生成很多的解，然后再反复这么做每次都保留最优解的话，总能找到解的。（当然了，遗传算法就是利用了电脑计算速度很快的特点，依据“适应度函数“来诱导电脑进行”好“的计算，从而更加快速的找到优解）

同时呢，部分染色体还需要进行交叉操作，对应生物界遗传过程中的“交叉”：就是把相应部分的基因交换。再以编码为例，比如现在随便从100个x值中选取两个吧，假设正好选中了x=3和4，对应的编码假设是11001 10101 和00101 01011，那么怎么交叉呢？我们知道每次交叉的染色体通常是一块一块的，恩，这里在算法设计上也来一块一块的吧。比如说就把位置在2,3,4号的编码给整体交叉了，那么x=3对应的位置是100吧，x=4对应的位置是010吧，好，交换以后x=3对应的位置就变成了010，x=4对应的位置就变成100，加回去就变成什么了？x=3是不是就是10101 10101，x=4是不是就是01001 01011了。而现在，把他们再反编码回去还是x=3和x=4吗？显然又不是了吧（当然也有小概率是一样的吧，很小）。从而又导致了新的个体产生。

同理，带到适应度函数里面去吧，再取优秀个体。同样，有些专门研究这种算法的开发出来各种各样的交叉方式，什么一个个体的前3个与后一个个体的后三个交叉，中间几位来交叉等等，总之就是生产新个体。而这样做的目的在哪呢？无非是三个字，随机性，充分保证生产新个体具有随机性，这对于算法的广阔搜索能力来说是非常好的。

6、关于选择的问题

说完了上面的部分，再说说选择吧，选择是什么？就是优胜劣汰。”好的留下来，差的走人“，在自然界中就是“人比人死，物比物扔”，不停地选择使的种群一直朝着较好的方向行走。

7、还差点什么呢

至此，在选择完后在重复上述步骤交叉，变异等等，那么什么时候可以结束呢？很简单，办法就是迭代次数。设置不同的迭代次数，比如迭代10次看一下结果，20次，看一下结果，30次、40次、100次。。。。。。当次数达到一定程度以后，优秀的个体越来越多，大都集中在最优解附近，即使变异或者交叉了也是在这个最优解附近。那么至此遗传算法就结束了。

8、代码

（1）主函数

function main()
clear;
clc;
%种群大小
popsize=100;
%二进制编码长度
chromlength=10;
%交叉概率
pc = 0.6;
%变异概率
pm = 0.001;
%初始种群
pop = initpop(popsize,chromlength);
for i = 1:100
    %计算适应度值（函数值）
    objvalue = cal_objvalue(pop);
    fitvalue = objvalue;
    %选择操作
    newpop = selection(pop,fitvalue);
    %交叉操作
    newpop = crossover(newpop,pc);
    %变异操作
    newpop = mutation(newpop,pm);
    %更新种群
    pop = newpop;
    %寻找最优解
    [bestindividual,bestfit] = best(pop,fitvalue);
    x2 = binary2decimal(bestindividual);
    x1 = binary2decimal(newpop);
    y1 = cal_objvalue(newpop);
    if mod(i,10) == 0
        figure;
        fplot('10.*sin(5.*x)+7.*abs(x-5)+10',[0 10]);
        hold on;
        plot(x1,y1,'*');
        title(['迭代次数为n=' num2str(i)]);
        %plot(x1,y1,'*');
    end
end
fprintf('The best X is --->>%5.2f\n',x2);
fprintf('The best Y is --->>%5.2f\n',bestfit);

（2）编码函数

%二进制转化成十进制函数
%输入变量：
%二进制种群
%输出变量
%十进制数值
function pop2 = binary2decimal(pop)
[px,py]=size(pop);
for i = 1:py
    pop1(:,i) = 2.^(py-i).*pop(:,i);
end
%sum(.,2)对行求和，得到列向量
temp = sum(pop1,2);
pop2 = temp*10/1023;

（3）计算目标函数值

%计算函数目标值
%输入变量：二进制数值
%输出变量：目标函数值
function [objvalue] = cal_objvalue(pop)
x = binary2decimal(pop);
%转化二进制数为x变量的变化域范围的数值
objvalue=10*sin(5*x)+7*abs(x-5)+10;

（4）交叉函数

%交叉变换
%输入变量：pop：二进制的父代种群数，pc：交叉的概率
%输出变量：newpop：交叉后的种群数
function [newpop] = crossover(pop,pc)
[px,py] = size(pop);
newpop = ones(size(pop));
for i = 1:2:px-1
    if(rand<pc)
        cpoint = round(rand*py);
        newpop(i,:) = [pop(i,1:cpoint),pop(i+1,cpoint+1:py)];
        newpop(i+1,:) = [pop(i+1,1:cpoint),pop(i,cpoint+1:py)];
    else
        newpop(i,:) = pop(i,:);
        newpop(i+1,:) = pop(i+1,:);
    end
end

（5）初始化种群

%初始化种群大小
%输入变量：
%popsize：种群大小
%chromlength：染色体长度-->>转化的二进制长度
%输出变量：
%pop：种群
function pop=initpop(popsize,chromlength)
pop = round(rand(popsize,chromlength));
%rand(3,4)生成3行4列的0-1之间的随机数
%round四舍五入
%返回的种群就是每行是一个个体，列数是染色体长度

（6）变异函数

%关于变异
%函数说明
%输入变量：pop：二进制种群，pm：变异概率
%输出变量：newpop变异以后的种群
function [newpop] = mutation(pop,pm)
[px,py] = size(pop);
newpop = ones(size(pop));
for i = 1:px
    if(rand<pm)
        mpoint = round(rand*py);
        if mpoint <= 0
            mpoint = 1;
        end
        newpop(i,:) = pop(i,:);
        if newpop(i,mpoint) == 0
            newpop(i,mpoint) = 1;
        else newpop(i,mpoint) == 1
            newpop(i,mpoint) = 0;
        end
    else newpop(i,:) = pop(i,:);
    end
end

（7）选择函数

%如何选择新的个体
%输入变量：pop二进制种群，fitvalue：适应度值
%输出变量：newpop选择以后的二进制种群
function [newpop] = selection(pop,fitvalue)
%构造轮盘
[px,py] = size(pop);
totalfit = sum(fitvalue);
p_fitvalue = fitvalue/totalfit;
p_fitvalue = cumsum(p_fitvalue);%概率求和排序
ms = sort(rand(px,1));%从小到大排列
fitin = 1;
newin = 1;
while newin<=px
    if(ms(newin))<p_fitvalue(fitin)
        newpop(newin,:)=pop(fitin,:);
        newin = newin+1;
    else
        fitin=fitin+1;
    end
end

（8）最适应函数

%求最优适应度函数
%输入变量：pop:种群，fitvalue:种群适应度
%输出变量：bestindividual:最佳个体，bestfit:最佳适应度值
function [bestindividual,bestfit] = best(pop,fitvalue)
[px,py] = size(pop);
bestindividual = pop(1,:);
bestfit = fitvalue(1);
for i = 2:px
    if fitvalue(i)>bestfit
        bestindividual = pop(i,:);
        bestfit = fitvalue(i);
    end
end

注：该代码在MATLAB 2019b编译完成（主要是为了完成Y老师作业编写，内容比较粗造而且有借鉴）

如果大家不想自己自己编写，当然MATLAB贴心的有遗传工具箱。而使用遗传工具箱主要就是编写句柄函数：

调用格式：

fun=@(para)function_0(para,t_b,c_b);

样例脚本： 

function y=ga_curfit(x)
global ydata n
t=1:n;
y=0;
for i=1:n
    y=y+(ydata(i)-(x(:,1)+x(:,2)+x(:,3).*log(t(i)+x(:,4)))).^2/n;
end
y=sqrt(y);
end

GA_RUN：

global ydata n
format long g
b=RH(:,2);
ydata=b';
n=length(ydata);
avg=sum(ydata)/n;
LB=[2/3*avg 0 0 0];
UB=[1*avg 1/3*avg 24*pi 2*pi];
nvars=4;
options=gaoptimset;
options=gaoptimset(options,'PopulationSize',300);
options=gaoptimset(options,'CrossoverFraction',0.8);
options=gaoptimset(options,'MigrationFraction',0.1);
options=gaoptimset(options,'Generations',500);
options = gaoptimset(options,'TolFun', 1e-50);
%options = gaoptimset(options,'InitialPopulation',final_pop);
options = gaoptimset(options,'Display', 'final');
options = gaoptimset(options,'PopInitRange', [LB;UB]);
options = gaoptimset(options,'PlotFcns', @gaplotbestf);
options=gaoptimset(options,'Vectorize','on');%目标函数向量化
[x,fval,exitflag,output,final_pop,scores]=ga(@ga_curfit,nvars,[],[],[],[],LB,UB,[],options);
t=1:n;
plot(t,ydata,'r*');
hold on
plot(t,(x(1)+x(2)+x(3).*log(t)+x(4)))
%plot(t,x(1)+x(2)*sin(x(3)*t+x(4)))
legend('数据','拟合')