本人是完全自学GEP的,由于时间紧张,又无人指导,所以用Matlab编写的GEP程序。下面,首先,我将具体给出GEP的流程及实现算法;其次,将我在编写过程中遇到的一些实现方式的Matlab代码给出,希望能对GEP初学者有帮助。本人采用的数值常量控制方式为GEP-RNC,所以很多问题是针对这个来展开的。时间有限,陆续添加~
注意:本文中所述文字及图标,均为自己按理解描述,可能会与现有论文表述不是完全一致。如果要放入论文或公开发表的刊物上,尤其是大段落直接摘录引用时,请注明本人的原出处,否则可能引起版权纠纷。
[1] 雷七峰. 改进的基因表达式变成算法及应用. 硕士学位论文,西安电子科技大学, 2010.
§·GEP的流程及实现算法
◆种群
GEP算法是由种群初始化开始的。由于初始化是随机的,因而有可能初始种群中所有个体(与染色体对应,本文中不严格区分二者)对应的解都是非法的。Ferreira指出,此时只需要反复重新初始化种群,直到产生至少有一个合法个体为止。而初始化种群中只要有一个个体对应的解是合法的,算法都可以开始,并产生合法的解。
GEP算法采用了简单精英策略(simple elitism)(即复制最佳个体直接进入新一代种群),以防止通过遗传操作而导致种群中的最优个体消失。又因为算法中使用了轮盘赌选择,所以GEP算法的选择策略又称为“保留精英的轮盘赌策略(roulette-wheel sampling with simple elitism)”。
◆染色体编/解码
编码:一条就是一个染色体是线性字符串,这就是染色体的编码形式;具体到单基因、多基因会有关于编码的定义和组织形式。由于是字符串,因而在matlab中应该考虑用1×n或者n×1的单元格结构来存储一条染色体。参考下文中“函数符集F(functionset)和终结符集T(terminalset)的设置”部分代码,该代码中的Fall就是采用单元格结构来存储的。
解码:对代表染色体的线性字符串进行解码。一般有两种做法,但实现起来多以较简单的线性解码法来实现。具体请参考文献《基因表达式编程的2 种解码方法》,作者:谢大同,陈巧云。
只要确定了染色体的MATLAB编码结构,解码参考上文,很容易实现。
◆适应度计算
适应度计算规则根据应用有不同的表述,一般是染色体解码结果的一个函数式。可以参考鼻祖ferreira的GEP著作fitness部分。
需要注意的是,染色体在解码中有可能会遭遇累如1/0的情况,需要对这部分染色体的适应度单独强制赋值为0或其他值,以使其在后续的进化中不占优势。
◆遗传操作
无论如何修改基因中的任何一位或几位,都可以保证该基因是合法的,即该基因一定能被编码为合法的表达式树。这个特点称为GEP对遗传操作的封闭性。Ferreira提出了以下几种遗传操作算子:选择算子,复制算子,变异算子,反转算子,插串算子和重组算子。Ferreira指出,GEP算法中遗传算子的先后顺序对最终的结果不是十分重要。
这个在matlab中实现是很容易的。
◆轮盘赌选择
参考本人的文章《GEP中的轮盘赌选择程序(MATLAB)》:
◆新一代种群
新种群是上一代种群按照一定的挑选规则挑选出来的,具体的规则可以自己设计。但总体来讲,需要遵循“优胜略汰”的进化原则,使得上一代种群中适应度大的个体进入新种群的概率占优。
§·GEP的matlab实现重要细节
1.自变量的设置。
函数发现,就是找出一个因变量与若干个自变量的函数关系,这也是GEP的一个重要应用。GEP中需要用字符或字符串来表示每个自变量,以及函数操作符号。有人用a,b,c来表示这些自变量。这样,是可以的,但更好的做法是手工输入3,自动分配三个字符或字符串。这样做,可以避免,如果有50个自变量,手动去输入,会非常麻烦。我们不妨用x1,x2,x3,x4,...来表示这些自变量。
代码:
Nx=10;%自变量个数
for p=1:Nx,
T{p}=strcat('x',int2str(p));
end
T
运行结果:
T =
'x1' 'x2' 'x3' 'x4' 'x5' 'x6' 'x7' 'x8' 'x9' 'x10'
2.函数符集F(functionset)和终结符集T(terminalset)的设置。
GEP中需要设置函数符集F和终结符集T,用于进行个体编码。最简单的GEP的单个基因分为头部和尾部,头部元素是从F或者T中随机获取,而尾部是从T中获取。
代码:
Fall={'+','-','*','/','sqrt','log','exp','sin','cos','tan'};%所有的可用的函数符号
Ftype=[1 2 3 8];%程序运行时需要手动输入的参数,表示选取Fall中的哪些元素作为F
F=Fall(Ftype);
Nx=10;%自变量个数
for p=1:Nx,
T{p}=strcat('x',int2str(p));
end
T=[T {'?'}];%对于GEP-RNC来讲,需要用一个'?'用来表示数值常量
fandt=[T F];%F和T的合集,头部元素从其中选取
T
F
fandt
运行结果:
T =
'x1' 'x2' 'x3' 'x4' 'x5' 'x6' 'x7' 'x8' 'x9' 'x10' '?'
F =
'+' '-' '*' 'sin'
fandt =
'x1' 'x2' 'x3' 'x4' 'x5' 'x6' 'x7' 'x8' 'x9' 'x10' '?' '+' '-' '*' 'sin'
3.个体的数据结构
对于GEP-RNC来讲,单基因分为3部分:头部,尾部和Dc域。此外,还需要一个随机常量数组C。Dc域和C用来控制数值常量。Dc域的长度等于尾部的长度,而尾部的长度是由头长和函数符集F共同决定的。Dc域的基因为一些数字,假设数组C长度为Lc,则Dc域的取值范围为0~9。通常,Lc=10。Dc域的0表示C的第0个元素。
在编写程序时,我们会发现,Dc域为数字,而头部和尾部为字符,要严格按照算法所述,将三者放在一起,需要先将Dc域转换为字符类型。另一方面,GEP的很多操作算子是针对头部和尾部的,而针对Dc域也有专门的操作算子,也就是说,头部和尾部是放在一起进化的,而Dc域是单独进化的。那么,我们不妨在编程时,利用结构体,下设三个域,分别存储:头部+尾部,Dc域,数组C。这样,也方便操作算子的操作。C++语言里面存在第0个元素,而matlab里面不存在第0个元素,为了清楚起见,不妨令Dc域的取值范围为1~10。
代码:
Lc=10;%随机常量数组C长度,一般取10即可
LowerBound=-3;UpperBound=4;%随机数值常量的取值范围
HeadSize=3;TailSize=4;%头长为3,尾长为4(尾长其实是根据头长计算出来的,为方便这里不给出计算公式而直接给出)
chrom{1}.genes={'+' '*' 'x1' 'x4' 'x5' 'x6' 'x5'}
chrom{1}.DC=ceil(rand(1,TailSize)*Lc);
chrom{1}.RNC=rand(1,Lc)*(UpperBound-LowerBound)+LowerBound;
chrom{1}.genes
chrom{1}.DC
chrom{1}.RNC
运行结果:
ans =
'+' '*' 'x1' 'x4' 'x5' 'x6' 'x5'
ans =
3 2 1 8
ans =
0.1157 3.5227 0.2620 -0.0695 2.9235 0.6761 -1.5815 1.7050 2.8668 -2.8625
4.