遗传算法解决函数优化问题

  • 作者: Cukor丘克
  • 环境: MatlabR2020a + vscode

为什么要学习遗传算法

为什么要学习遗传算法，或者说遗传算法有什么厉害的地方。例如求解以下函数优化问题：
\(min f(x_1, x_2)=x^2_1+x^2_1+25*(sin^2x_1+sin^2x_2), -10 \le x_1 \le 10, -10 \le x_2 \le 10 .\)

方法一：

这个时候你可能想从设一个变量\(i\)从-10到10，一个变量\(j\)从-10到10，每个涉及到的点都遍历，然后比较前后数值取最大的那个，当遍历完成了就确定了最大值。

方法二：

对函数求导，然后令导为零，取最小值。

两种方法在一般理论上都可行。但是有各自的缺点。

方法一的缺点：

每次都从头开始遍历，直到把所有的都遍历完才能确定最值，太耗时间，而且大多数做的都是没有意义的比较。而且计算机对浮点数的比较有精度的要求，有时候两个数只差一点点，计算机就会将它们当作是相等的。

方法二的缺点：

不是所有函数的导函数都是容易求解的，而且当出现很多个峰值或谷值的时候也是比较麻烦，就展现不出求导的优势了。

这个时候就需要像遗传算法这样的智能优化算法来求解。

目标函数的图像

遗传算法的基本认识演化计算

生物进化是指一个种群经过漫长的时间所发生的累积变化，这些变化是由于生物体的基因变异或在繁殖期间以不同方式重组基因所产生的，而且这些变化可以被遗传到生物体的后代。

生物的进化可以看成是一个优化过程，而优化的结果是产生能够很好地适应环境的生物体。现在地球上的种类繁多‘结构复杂的生物都是通过漫长的由简单到复杂、由低级到高级的进化过程而得到的优化结果。生物的进化过程也可以看成是在众多可能性中搜索“解”的一种方法。

演化算法

演化计算所涉及的算法称为演化算法。所有的演化算法都有一个共同点：求解问题的过程也就是模拟大自然生物进化的过程。

演化算法模仿自然进化过程，在求解问题的过程中，保持一个个体的种群，每个个体表示问题的一个可能解。个体适应环境的程度用一个适应函数判断，每个个体安装适应函数来度量该个体作为问题解的好坏程度。

遗传算法是应用最为广泛的一种演化算法。遗传算法是美国密歇根大学的J.H.Holland教授在研究自然界自适应现象的过程中提出来的。

算法流程图

参数初始化主要影响遗传算法的四个参数

  1. 种群大小（20~200）
  2. 交叉概率（0.7~0.9）
  3. 变异概率（0.02~0.2）
  4. 迭代次数（20~500）

次要的参数

  1. 下限
  2. 上限
  3. 保留小数点后t位
  4. 编码长度数组
  5. 个体编码长度

个体编码

设计演化算法的第一步是对问题的可能解进行编码，其目的是为了能够有效地执行遗传操作。

演化算法不是直接作用在问题的解空间上的，而是交替地作用在编码空间和解空间上。

编码是一个从问题的解空间到编码空间的映射。

编码可以是二进制编码也可以是Gray码编码，一般就使用二进制编码。

一个编码就对应着一个实数解，一个实数解对应一个编码。这就需要编码的长度一定要足够大，这样才能把所有的可能解都使用编码表示 。

实例：

设\(a_j \le x_j \le b_j\)所要求的精度为小数点后t位，这要求将区间\([a_j, b_j]\)划分为至少\((b_j-a_j)10^t\)份。假设表示变量\(x_j\)的位串长度用\(l_j\)表示，则\(l_j\)可取为满足下列不等式的最小正整数\(m\):$$(b_j-a_j)10^t \le 2^m-1,$$

将\(x_j\)的二进制表示转换为十进制表示可按下式计算：$$x_j=a_j+decimal(substring_j)*(b_j-a_j)/(2^k-1), k=l_j,$$

编码

现有\(-3.0 \le x_1 \le 12.1\),精度为小数点后t位，\(t=4\),则带入上述公式$$(12.1-(-3.0))*10^4=151000,$$

\[2^n-1<151000 \le 2^m-1，n=17,m=18,\]

所以表示\(x_j\)的二进制位串的长度为\(l_1=18,\)

同理\(4.1 \le x_2 \le 5.8\)可以使用\(l_2=15\)的二进制位串来表示。

而一个个体是由\(x_1、x_2\)共同表示的，所以种群中的每个个体应该使用\(l_1+l_2=33\)的二进制位串表示。

解码

给定下列33位二进制位串：$$011011100110001010110011011001101,$$

那么表示前18位所表示的变量\(x_1\)的值为$$x_1=-3.0+decimal(011011100110001010)*(12.1-(-3.0))/(2^k-1), k=18,$$

求得\(x_1=-3.0+113034*15.1/262143=3.5110\)

对\(x_2\)解码同理。

计算个体的编码长度的Matlab代码如下：

% 主函数中如下调用L = getcodelen(a, b, t);    % 编码长度数组len = sum(L);               % 编码长度% 返回的是一个编码长度数组，后续代码需要用到它，所以把它和个体的编码长度分开function L = getcodelen(a, b, t)%getcodelen - 计算a到b的十进制数的最大二进制编码长度%% Syntax: L = getcodelen(a, b, t)%% Long description    n = length(a);    L = zeros(n, 1);    for i = 1:n       while ((b(i) - a(i))*10.^t) > ((2^L(i)) - 1)          L(i) = L(i) + 1;       end    endend

初始化种群

因为使用的是二进制编码，所以只需要确定种群的大小和编码的长度即可初始化一个种群。然后不断迭代产生最优的个体。

初始化种群的Matlab代码如下:

function pop = initpop(popsize, len)    pop = round(rand(popsize, len));end

疑问环节：

问：因为种群是随机生成的，有没有一种可能，初始化的个体编码不能和变量的定义域对应。

答：没有这种可能。因为在设计编码的时候就是按照定义域a到b的范围设定的，所以不管怎么随机得到的二进制位串，解码之后都是在a到b的区间内。并且在解码的公式中可以看出不会将二进制位串解码到a到b的区间之外。x=a+….就表示了最小值一定是a，如果\(decimal(substring)/(2^l-1)=1\)，那么整个式子就剩下\(x=a+b-a=b\)，所以最大也就只能到b.

计算适应值

适应函数是区分种群个体好坏的唯一方式，是进行选择的基础。

在设计适应函数时，应遵守以下原则：

  1. 最优解与具有最大适应值的个体相对应
  2. 适应值能够反映个体质量的差异
  3. 计算量应尽可能小

此外，有些选择策略，如轮盘赌选择，还要求适应值非负。

常见的适应函数

  • 原始适应函数。原始适应函数是直接由目标函数变化而来。当优化问题为\(max f(x)\)时，适应函数为\(f(x)\)，当优化问题为\(min f(x)\)时，适应函数为\(-f(x)\).
  • 简单适应函数。为了防止适应值为负的情形，常常需要对目标函数作简单变换。对于优化问题\(max f(x)\)，适应函数可以定义为$$1/(cmax-f(x)),$$

其中，\(cmax>f(x)\).对应优化问题\(min f(x)\),适应函数可以定义为$$1/(f(x)-cmin),$$

其中，\(f(x)>cmin\).

function fitvalue = fitness(pop, f, a, b, L, cmin)%fitness - 计算适应值%% Syntax: fitvalue = fitness(pop, f, a, b, L, cmin)%% Long description    [px, ~] = size(pop);    fitvalue = zeros(px, 1);    x = zeros(px, length(L));    for i = 1:px        x(i, :)  = decode(pop(i, :), a, b, L);        fitvalue(i) = 1./(f(x(i, 1), x(i, 2)) - cmin);    endend

父体选择

选择或复制操作是决定哪些个体可以进入下一代。选择哪些适应值比较高的个体。我们采用选择策略是轮盘赌选择，因为这个方法比较简单，也比较直观。

轮盘的设计

  1. 累加适应值，得到一个常数F
  2. 每个适应值除以F，得到了每个个体的适应值所占整体的份数\(p_i\)
  3. 累加每个\(p_i\)，实现摆盘操作

转动轮盘

随机地产生一些随机数，相当于转动了轮盘。可以使用rand函数来实现。

function newpop = parent_select(pop, fitvalue)%parent_select - 父体选择%% Syntax: newpop = parent_select(pop, fitvalue)%% Long description    [px, py] = size(pop);    newpop = zeros(px, py);    %% 轮盘的设计    p = cumsum(fitvalue./sum(fitvalue));    %% 转动轮盘    r = sort(rand(px, 1));  % 排序是为了减少频繁变动轮盘的指针    j = 1;    for i = 1:px        while r(i) > p(j)            j = j + 1;        end        % r(i) <= p(j)        newpop(i, :) = pop(j, :);    endend

交叉

遗传算子的设计是与编码密切相关的。因为之前采用的编码是二进制编码，所以后面提到的遗传算子的交叉、变异都是基于二进制位串的。

  1. 点式杂交。
    1. 单点杂交
    2. 多点杂交

本次案例采用的杂交方式是单点杂交。

单点杂交算法过程：

设二进制位串的长度是\(L\)，首先随机地产生一个整数\(pos\)作为杂交点的位置，\(1 \le pos \le L\)，然后将两个父体在该杂交点右边的子串进行交换，产生两个后代个体。

例如：给定两个父体如下：

$v_1=(100111000101\ |\ 01011001110)， $

$v_2=(100010111100\ |\ 11001100101). $

假设杂交点的位置是13，也就是上面画 |之后的二进制位串进行交叉。结果如下：

$v’ _1 =(100111000101\ |\ 11001100101)， $

$v’ _2 =(100010111100\ |\ 01011001110)， $

多点杂交算法过程：

设二进制位串的长度是\(L\)，多点杂交在$1~ ~\ ~L-1 $之间随机地选择多个杂交点，然后再保持第一个杂交点左边的对应字串不交换的情形下，间隔地交换两个父体在杂交点之间的对应字串，生成两个后代 。

例如，给定两个父体如下：

\(f_1=(1000111110~|~10101110101~|~001001),\)

\(f_2=(1011011101~|~00100110100~|~110101).\)

假设所选择的两个杂交点分别为11和22，那么经两点杂交后，所得到的两个后代个体如下：

\(s_1=(1000111110~|~00100110100~|~001001),\)

\(s_2=(1011011101~|~10101110101~|~110101).\)

在多点杂交中，有三段杂交的也是类似的操作。

  1. 均匀杂交。

均匀杂交是依概率交换两个父体位串的每一位，其过程如下：

先随机地产生一个与父体等长的二进制位串，其中，0表示不交换，1表示交换，这个二进制位串称为杂交模板，然后根据所产生的杂交模板对两个父体进行杂交。

例如，给定两个父体如下：

\(f_1=(01001),\)

\(f_2=(10101).\)

假设所产生的杂交模板为：

\(（00101）\)

所得到的两个后代如下：

\(s_1=(01101),\)

\(s_2=(10001).\)

下面给出单点杂交的Matlab代码，其他的两种杂交算法由你自己实现，也很好实现的。

function newpop = crossover(pop, L, pc)%crossover - 交叉%% Syntax: newpop = crossover(pop, L, pc)%% Long description    [px, ~] = size(pop);    newpop = pop;    start = 3;  % 奇数的情况，从第3个个体开始    if mod(px, 2) == 0        start = 2;    end    n = length(L);    for i = start:2:px        if rand <= pc        rear = 0;            for j = 1:n                front = rear + 1;                rear = front + L(j) - 1;                [newpop(i, front:rear), newpop(i-1, front:rear)] ...                  = sub_crossover(pop(i, front:rear), pop(i-1, front:rear));            end        end    endendfunction [newcode1, newcode2] = sub_crossover(code1, code2)%sub_crossover - 子交叉函数%% Syntax: [newcode1, newcode2] = sub_crossover(code1, code2)%% Long description    n = length(code1);    pos = -1;           % 交叉点    while pos < 1        pos = ceil(rand * n);    end    newcode1 = code1;    newcode2 = code2;    newcode1(pos:n) = code2(pos:n);    newcode2(pos:n) = code1(pos:n);end

变异

二进制编码时的变异算子非常简单，只是依一定概率\(p_m\)(变异概率)将所选个体的位串取反，即若是1则取0；若是0，则取1.

变异的Matlab代码

function newpop = mutation(pop, pm)%mutation - 变异%% Syntax: newpop = mutation(pop, pm)%% Long description    [px, py] = size(pop);    newpop = pop;    cnt = ceil(rand*py/2);      % 确定变异点的个数    pos = mod(ceil(rand(px,cnt) * py), py) + 1;  % 变异点的位置    for i = 1:px        if rand <= pm            for j = 1:cnt                newpop(i, pos(i, j)) = ~newpop(i, pos(i, j));            end        end    endend

到此为止，已经将遗传算法的主体部分介绍完毕，接下来就是回到实际问题中，使用遗传算法解决它，即绘图分析、结果分析之类的。

回到问题

原问题是求解以下函数优化问题：
\(min f(x_1, x_2)=x^2_1+x^2_1+25*(sin^2x_1+sin^2x_2), -10 \le x_1 \le 10, -10 \le x_2 \le 10 .\)

那就把图画出来，用遗传算法把每一代的最优个体所对应图像上的点绘制出来，最后在历代最优个体中选出适应值最高的个体，这个个体（可能解）就当作是函数\(f(x_1,x_2)\)的全局最优解。

下面给出本次案例的主函数的matlab代码。所涉及到的其他调用子函数已经在上面实现，在自己编写的时候记得自己另创建一个对应的*.m文件。

%% 遗传算法解决优化问题%% 清屏clear; clc;%% 目标函数f = @(x1, x2)(x1.^2 + x2.^2 + 25*(sin(x1).^2 + sin(x2).^2));%% 参数初始化% 主要影响遗传算法的四个参数popsize = 100;   % 种群大小pc = 0.8;        % 交叉概率pm = 0.05;       % 变异概率gen = 200;       % 迭代次数% 次要的参数a = [-10.0 -10.0];  % 下限b = [10.0 10.0];    % 上限t = 4;              % 保留小数点后t位L = getcodelen(a, b, t);    % 编码长度数组len = sum(L);               % 个体编码长度cmin = -100;% 与画图有关的变量x = zeros(gen, length(L));  % 每一代的最优解的x分量z = zeros(gen, 1);  % 每一代的最优解的适应值%% 初始化种群pop = initpop(popsize, len);%% 开始迭代for it = 1:gen    fitvalue = fitness(pop, f, a, b, L, cmin);      % 计算适应值    [z(it), best_index] = max(fitvalue);            % 记录每一代的最优个体（最高适应值、下标）    x(it, :) = decode(pop(best_index, :), a, b, L);    newpop = parent_select(pop, fitvalue);          % 父体选择    newpop = crossover(newpop, L, pc);              % 交叉    newpop = mutation(newpop, pm);                  % 变异    pop = newpop;end%% 打印最优解[best_fitvalue, best_index] = max(z);best_solution = x(best_index, :);best_value = f(best_solution(1), best_solution(2));disp(['最优解为: ' num2str(best_solution)]);disp(['最优适应值为：' num2str(best_fitvalue)]);disp(['最小值为：' num2str(best_value)]);%% 画图figure(1)[X, Y] = meshgrid(a(1):0.05:b(1), a(2):0.05:b(2));Z = f(X, Y);mesh(X, Y, Z);text(x(:, 1), x(:, 2),f(x(:,1), x(:, 2)), 'o', 'Color','blue', 'Fontsize', 8);text(best_solution(1), best_solution(2), best_value, 'O', 'Color', 'red', 'FontSize', 15)hold on;best = z;for i = 2:gen    best(i) = max(best(i), best(i-1));endfigure(2)plot(z);figure(3);plot(best);

结果图

图1 原函数图像及各代最优个体散点

图2 适应值变化曲线

图3 历代适应值最大的适应值变化曲线

从图1可以看出，通过遗传算法计算出的最优解基本贴近原函数的的最小值的位置，说明这个最优解可行。从图2和图3可以看到遗传算法是可以收敛到最优解的。本次运行程序大概在第145代左右确定全局最优解，虽然不是标准最优解\((0,0)\)，但是这个最优解已经贴近标准最优解了。

同时本文章也有配套的视频讲解，下面是笔者在bilibli发布的遗传算法的讲解视频链接:

Cukor丘克-《通过例子学习遗传算法（遗传算法解决函数优化问题）》

参考文献

[1] 黄竞伟, 朱福喜, 康立山. 《计算智能》. 科学出版社,2010.6.

[2] 温正. 《智能优化算法》.