遗传算法的matlab实现（续 ）

       上一篇博客中代码部分有人应该会有疑惑，在这里我先做一些介绍：

• 自然选择部分使用了排名法，适应度越高越难被淘汰，但是适应度最高的不会被淘汰；
• 里面出现了很多个个体，其中父代是一个种群N，交叉之后有了两种子代2N，变异之后也会出现子代N，共4N，因此自然选择需要淘汰到只剩N个个体，作为下一代的父代。 

       不过关于多维函数极值的遗传算法写法，很多人直接将每一个变量单独作为一个种群，以此并行遗传，虽然这样可行，但是从代码编写上很不美观，而且当变量很多时，这种策略显然不可行。下面我也给出了一种多维函数极值的遗传算法程序：

        问题背景是求z=2-exp(-x^2-y^2)的极值（最大值和最小值）。

clc;clear;close all;%% 参数设定N = 50;                    %种群数量ger = 100;                 %迭代次数L = [5 5];                 %基因长度pc = 0.8;                  %交叉概率pm = 0.1;                  %变异概率R = [-5 5;-5 5];           %变量范围%% 初始化%初始化解码器decode = cell(length(L),1);for k = 1 : length(L)    decode{k}= zeros(L(k),1);    for i = 1 : L(k)        decode{k}(i) = 10^(L(k)-i);    endend%初始化种群x = randi([0 9],[N,sum(L)]);     %初始化父代种群x1 = x;                          %初始化交叉子代种群1x2 = x;                          %初始化交叉子代种群2x3 = x;                          %初始化变异子代种群a = zeros(4*N,length(L));        %中间种群%% 开始遗传for i = 1 : ger    for j = 1 : N        if rand < pc            p = randi(N);        %确定交叉对象            start = 1;            for k = 1 : length(L)                q = randi(L(k)-1);%确定交叉节点                x1(j,start:sum(L(1:k))) = [x(j,start:sum(L(1:k-1))+q),x(p,start+q:sum(L(1:k)))];                x2(j,start:sum(L(1:k))) = [x(p,start:sum(L(1:k-1))+q),x(j,start+q:sum(L(1:k)))];                start = start + L(k);            end        end        if rand < pm            start = 0;            q = randi(L(k));%确定变异节点            for k = 1 : length(L)                x3(j,start+q) = randi([0 9]);                start = start + L(k);            end        end    end    x = [x;x1;x2;x3];%合并父代子代    start = 1;    for k = 1 : length(L)        a(:,k) = x(:,start:sum(L(1:k)))*decode{k}*(R(k,2)-R(k,1))/(10^L(k)-1)+R(k,1);%解码        start = start + L(k);    end    y = -f(a);%计算适应度    z = [x,y];    z = sortrows(z,sum(L)+1);%对适应度排序    while size(z,1) > N        n = randi(size(z,1));        if rand > n/size(z,1)%自然选择，排名法            z(n,:) = [];        end    end    z(:,sum(L)+1) = [];%去除最后一列（适应度）    x = z;    t(i)=-max(y);endstart = 1;xm = zeros(1,length(L));for k = 1 : length(L)      xm(k) = x(end,start:sum(L(1:k)))*decode{k}*(R(k,2)-R(k,1))/(10^L(k)-1)+R(k,1);%取最后一代适应度最高的个体      start = start + L(k);endxm(abs(xm)<1e-3)=0;ym =f(xm);disp(['最优解为：',num2str(xm)]);disp(['最优值为：',num2str(ym)]);plot(t);title('适应度曲线')

子函数：

function y = f(x)    y = 2 - exp(-x(:,1).^2-x(:,2).^2);end

      上面的代码求的是最小值，因此在求适应度时加上了负号（记录适应度的地方也加上了负号），如果要求最大值，自然只需要把负号去掉即可，结果如下：

     在这里我们发现了一个问题，最大值处对应的最优值不稳定，为何呢？我们先来做一个实验，编一个自定义穷举函数，大家可以自己试试看：

clc;clear;close all;x = -5 : 0.01 : 5;y = -5 : 0.01 : 5;x1 = repmat(x,length(x),1);x2 = x1(:);y1 = repmat(y',1,length(y));y2 = y1(:);z = f([x2,y2]);n = find(z == max(z));xm = [x2(n),y2(n)];ym = max(z)

     最后发现，最优值的确为2，但是对应的最优解却多达4w多种，从函数形式来看，由函数单调性可知当x,y为-5/5时，函数取最大值，这是什么原因呢？我们先画出函数的深度图：

       可以发现，除了中间部分，绝大部分的值都差不多，这是因为函数中有指数函数，并且matlab默认误差值为2e-16,二者综合作用下，很难找到真正的极值点，那么我们是否没办法呢？这类问题可以通过改造适应度函数来实现，当然这只能在可以明确适应度函数单调性的情况下使用，可以将适应度函数改为z0=x^2+y^2,那么我们最终就只需要做一下转换z=2-exp(-z0)即可，再套用上面的代码可轻松实现。