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

• 自然选择部分使用了排名法，适应度越高越难被淘汰，但是适应度最高的不会被淘汰；
• 里面出现了很多个个体，其中父代是一个种群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);
    end
end
%初始化种群
x = randi([0 9],[N,sum(L)]);     %初始化父代种群
x1 = x;                          %初始化交叉子代种群1
x2 = x;                          %初始化交叉子代种群2
x3 = 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);
end
start = 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);
end
xm(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)即可，再套用上面的代码可轻松实现。