用Matlab编写GEP程序

用Matlab编写GEP程序  

2010-12-24 21:25:49|  分类： 我眼看GEP |  标签：基因表达  轮盘赌选择  基因表达式编程  gep  expression   |字号大中小 订阅

本人是完全自学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）》：

       http://winsonman.blog.163.com/blog/static/10400910820113210117471/

 ◆新一代种群

        新种群是上一代种群按照一定的挑选规则挑选出来的，具体的规则可以自己设计。但总体来讲，需要遵循“优胜略汰”的进化原则，使得上一代种群中适应度大的个体进入新种群的概率占优。

§·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.

 

来源： http://winsonman.blog.163.com/blog/static/1040091082010112492549840/