遗传算法求数值函数的最值
0. 引言
设有函数:
f(x) = x + 10*sin(5*x) + 7*cos(4*x);
其图像容易画出,如下所示:
先要求求该函数的最大值,读者可能已经有了很多种思路,本文介绍遗传算法是如何解决此类问题的。
1. 遗传算法简介
如果不关心算法的实现细节的话,遗传算法可以使用如下的流程描述。
这基本是借鉴生物种群的自然演化规律而抽象得到的流程图。下面分别解释流程图中的各个步骤。
编码
众所周知的,生物中都有保存其遗传信息的物质——染色体,染色体中的记录的遗传信息——基因决定了一个生物个体的表现性状,可以简单地认为一个染色体唯一确定一个生物个体。
染色体的本质是一种蛋白质,其结构非常复杂,用计算机模拟是不现实的。因此遗传算法简单地把染色体视为一串二进制位组成的序列(即该序列中的所有元素取自集合{0,1}
)。
假设有命名为A的生物,规定其染色体序列长度为8,这时可以取得一个个体:
[0,1,0,0,1,0,1,1]
初始化种群
显然一个生物种群不可能只包含一个个体,因此可以规定种群的大小,为方便说明原理,这是规定其大小为4。初始种群可以随机生成。如下:
[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
[1,0,1,0,0,0,0,1]
[0,0,0,1,0,1,0,1]
评估种群中的个体适应度
生物种群的存续与兴衰决定于其对于环境的适应能力。对于使用二进制位序列表示的个体,总可以找出一个函数来评估这个个体对环境的适应能力。可以规定函数:
f(x) = int(x)
# x是输入的一个个体,如[0,1,0,0,1,0,1,1]
# 返回值是该二进制序列表示的整数
此外还可以规定上述函数的返回值(一个整数)越大,个体的适应度越高。显然适应度最高的个体是:
[1,1,1,1,1,1,1,1]
选择
选择操作用于选出一个种群中存活概率更高的个体,这些个体组成一个新的种群。显然个体适应度越高的个体存活概率更高。因此可以使用赌轮盘算法选出新种群。赌轮盘算法的原理不是本文的重点,因此略去,其实现代码见后文。
交叉
交叉操作也是对生物遗传过程的抽象。设有两个生物个体,这两个个体的生殖繁衍过程中,其染色体会进行交叉交换,生成新的染色体,繁衍得来的新的个体携带新的染色体。
设有两个生物个体如下:
[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
设对这两个生物个体有一定的概率发生交叉操作,而交叉操作必然发生在二进制位序列中的某一位置。这一位置的选择也是随机的。假设从第5位(从1开始计数)开始交叉,其步骤如下:
[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
->
[1,1,1,1,1],[1,0,0]
[1,1,1,1,0],[1,1,0]
->
[1,1,1,1,1],[1,1,0]
[1,1,1,1,0],[1,0,0]
->
[1,1,1,1,1,1,1,0]
[1,1,1,1,0,1,0,0]
变异
生物在自然环境中受环境的影响,有可能发生遗传信息的永久改变,也就是变异。遗传算法规定,二进制位序列的变异指的是在一定的概率下,随机选择该序列中的一位取反。假设在第2位发生了变异,步骤如下:
[0,1,0,0,1,0,1,1]
->
[0],1,[0,0,1,0,1,1]
->
[0],0,[0,0,1,0,1,1]
->
[0,0,0,0,1,0,1,1]
小结
上述的各个步骤按照流程图所示依次执行,直到一个种群中所有的个体都相同(即收敛),如:
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
或者主循环的迭代次数超过规定次数,算法停止,输出结果。
Matlab入门
本文中的算法使用Matlab描述,因此需要对Matlab有一定的了解。但是帮读者入门Matlab不是我们的任务,因此推荐以下手册,可以帮助读者快速掌握Matlab:
https://ww2.mathworks.cn/help/pdf_doc/matlab/getstart_zh_CN.pdf
Matlab是收费软件,而且对于运行该软件的性能有一定要求,因此可以使用开源免费且占用资源较少的Octave作为替代品,其官网如下:
https://www.gnu.org/software/octave/
Octave使用的语法规则与Matlab基本相同,不兼容之处见下:
https://en.wikibooks.org/wiki/MATLAB_Programming/Differences_between_Octave_and_MATLAB
本文所附的代码在Octave上运行通过,并未在Matlab上进行测试,如果有需要使用Matlab运行代码,需要自己检查兼容性。
算法实现
function retval = main (input1, input2)
% 主函数,input1,input2,以及retval都没有实际作用
clear all;
close all;
clc;
NP = 50; % 种群大小
L = 16; % 染色体序列长度
Pc = 0.8; % 交叉概率
Pm = 0.1; % 变异概率
G = 50; % 主循环迭代次数
Xs = 10; % 所求函数的定义域的最大值
Xx = 0; % 所求函数的定义域的最小值
S = 10; % 什么意思自己悟吧,我只是不想硬编码成10
pop = initPop(NP,L);
for i = 1:G
dpop = decode(pop);
xpop = getX(dpop,L,Xx,Xs);
fit = fitValue(xpop);
if all(fit == fit(1)) % 如果收敛,停止主循环
break;
endif
[maxFit,maxIndex] = max(fit);
maxX(i) = getX(decode(pop(maxIndex,:)),L,Xx,Xs);
maxY(i) = maxFit;
pop = select(pop,fit,xpop);
pop = crossOver(pop,Pc);
pop = mutation(pop,Pm);
endfor
x = Xx:(Xs-Xx)/(2^L-1)/S:Xs;
y = targetFunc(x);
subplot(2,1,1);
plot(x,y); % 绘制函数图像
hold on;
plot(maxX,maxY,'r*'); % 绘制历次迭代生成的最大值点
[px,py] = size(maxY);
x = 1:py;
y = maxY;
subplot(2,1,2);
plot(x,y); % 绘制历次生成的最大值点关于迭代次数的折线图
hold off;
endfunction
function pop = initPop(NP,L)
pop = round(rand(NP,L));
endfunction
function res = targetFunc(x)
res = x + 10*sin(5*x) + 7*cos(4*x);
endfunction
function res = decode(pop)
% 用于将二进制位序列转化成整数
[px,py] = size(pop);
res = ones(px,1);
for j = 1:py
res(:,j) = 2.^(py-j)*pop(:,j);
endfor
res = sum(res,2);
endfunction
function res = getX(dpop,L,Xx,Xs)
% 缩放二进制位序列表示成的整数到
% 定义域中
res = dpop./(2^L-1).*(Xs-Xx);
endfunction
function res = fitValue(xpop)
res = targetFunc(xpop);
endfunction
function newpop = select(pop,fit,xpop)
% 赌轮盘算法的实现
[maxFit,maxIndex] = max(fit);
[minFit,minIndex] = min(fit);
[px,py] = size(pop);
newpop = ones(px,py);
fit = (fit-minFit)./(maxFit-minFit);
fit = fit./sum(fit);
cumfit = cumsum(fit);
ms = sort(rand(px,1));
fiti = 1;
newi = 1;
while newi <= px
if ms(newi) < cumfit(fiti)
newpop(newi,:) = pop(fiti,:);
newi = newi + 1;
else
fiti = fiti + 1;
endif
endwhile
endfunction
function newpop = crossOver(pop,Pc)
% 交叉操作
[px,py] = size(pop);
newpop = ones(size(pop));
for i = 1:2:px-1
if rand<Pc
cpoint = ceil(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,:);
endif
endfor
endfunction
function newpop = mutation(pop,Pm)
% 变异操作
[px,py] = size(pop);
newpop = ones(px,py);
for i = 1:px
if rand<Pm
mpoint=ceil(rand*py);
newpop(i,mpoint) = 1 - pop(i,mpoint);
newpop(i,:) = pop(i,:);
else
newpop(i,:) = pop(i,:);
endif
endfor
endfunction
绘图结果如下: