FitnV=ranking(-ObjV);%分配适应度值(Assignfitnessvalues)
SelCh=select('sus',Chrom,FitnV,GGAP);%选择
SelCh=recombin('xovsp',SelCh,0.7);%重组
SelCh=mut(SelCh);%变异
variable=bs2rv(SelCh,FieldD);%子代个体的十进制转换
ObjVSel=variable.*sin(10*pi*variable)+2.0;%计算子代的目标函数值
[ChromObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);%重插入子代的新种群
gen=gen+1;%代计数器增加
%输出最优解及其序号,并在目标函数图象中标出,Y为最优解,I为种群的序号
[Y,I]=max(ObjV),holdon;
plot(variable(I),Y,'bo');
trace(1,gen)=max(ObjV);%遗传算法性能跟踪
trace(2,gen)=sum(ObjV)/length(ObjV);
end
variable=bs2rv(Chrom,FieldD);%最优个体的十进制转换
holdon,grid;
plot(variable',ObjV','b*');
figure
(2);
plot(trace(1,:
)');
holdon;
plot(trace(2,:
)','-.');grid;
legend('解的变化','种群均值的变化')
使用基于适应度的重插入确保四个最适应的个体总是被连续传播到下一代。
这样在每一代中有36(NIND*GGAP)个新个体产生。
区域描述器FieldD描述染色体的表示和解释,每个格雷码采用20位二进制,变量区间为[-1,2]。
程序段Chrom=crtbp(NIND,PRECI)表示一个初始种群Chrom被函数crtbp创建,它是由NIND个均匀分布长度为PRECI的二进制串矩阵构成。
基于排序的适应度分配计算由程序段FitnV=ranking(-ObjV)实现。
对这个等级评定算法的缺省设置是选择等差为2和使用线性评估,给最适应个体的适应度值为2,最差个体的适应度值为0,这里的评定算法假设目标函数是最小化的,所以ObjV乘了一个负号,使目标函数为最大化。
适应度值结果被向量FitnV返回。
选择层使用高级函数选择调用低级函数随机遍历抽样例程sus,SelCh包含来自原始染色体的GGAP*NIND个个体,这些个体将使用高级函数recombin进行重组,recombin使个体通过SelCh被选择再生产,并使用单点交叉例程xovsp,使用交叉概率Px=0.7执行并叉。
交叉后产生的子代被同一个矩阵SelCh返回,实际使用的交叉例程通过支持使用不同函数名字串传递给recombin而改变。
为了产生一组子代,变异使用变异函数mut。
子代再次由矩阵SelCh返回,变异概率缺省值PM=0.7/Lind=0.0017,这里Lind是假定的个体长度。
再次使用bs2rv,将个体的二进制编码转换为十进制编码,计算子代的目标函数值ObjVSel。
由于使用了代沟,所以子代的数量比当前种群数量要小,因此需要使用恢复函数reins。
这里Chrom和SelCh是矩阵,包含原始种群和子代结果。
这两个事件的第一个被使用单个种群和采用基于适应度的恢复,基于适应度的恢复用SelCh中的个体代替Chrom中最不适应的个体。
新种群中的个体是由原始种群中的优良个体和子代中新产生的个体组成。
原始种群中个体的目标函数值ObjV随后又作为函数reins的输入参数,子代中个体的目标函数值由ObjVSel提供。
Reins返回具有插入子代的新种群Chrom和该种群中个体的目标函数值ObjV。
每次迭代后的最优解和解的均值存放在trace中。
这个遗传优化的结果包含在矩阵ObjV中。
决策变量的值为variable(I)。
画出迭代后个体的目标函数值分布图和遗传算法性能跟踪图。
遗传算法的运行结果如下:
(1)图7.1为目标函数
的图象。
图7.1目标函数图像
(2)图7.2为目标函数的图像和初始随机种群个体分布图。
图7.2初始种群分布图
(3)经过1次遗传迭代后,寻优结果如图7.3所示。
x=1.6357,f(x)=3.4729。
图7.3一次遗传迭代后的结果
(4)经过10次遗传迭代后,寻优结果如图7.4所示。
x=1.8518,f(x)=3.8489。
图7.4经过10次遗传迭代后的结果
(5)经过25次遗传迭代后,寻优结果如图7.5所示。
x=1.8505,f(x)=3.8503。
图7.5经过25次遗传迭代后的结果
(6)经过25次迭代后最优解的变化和种群均值的变化见图7.6。
图7.6经过25次迭代后最优解的变化和种群均值的变化
7.2多元单峰函数的优化实例
目标函数是DeJong函数,是一个连续、凸起的单峰函数,它的M文件objfun1包含在GA工具箱软件中。
DeJong函数的表达式为
求解
这里n是定义问题维数的一个值。
这个例子中选取n=20。
由DeJong函数的表达式可以看出,DeJong函数是一个简单的平方和函数,只有一个极小点(0,0,…,0),理论最小值为f(0,0,…,0)=0。
程序的主要变量:
个体的数量NIND为40,最大遗传代数为MAXGEN=300,变量维数为NVAR=20,每个变量使用20位表示,即PRECI=20,使用代沟GGAP=0.9。
下面为求解DeJong函数最小值的MATLAB代码。
%定义遗传算法参数
NIND=40;%个体数目(Numberofindividuals)
MAXGEN=500;%最大遗传代数(Maximumnumberofgenerations)
NVAR=20;%变量的维数
PRECI=20;%变量的二进制位数(Precisionofvariables)
GGAP=0.9;%代沟(Generationgap)
trace=zeros(MAXGEN,2);
%建立区域描述器(Buildfielddescriptor)
FieldD=[rep([PRECI],[1,NVAR]);rep([-512;512],[1,NVAR]);rep([1;0;1;1],[1,NVAR])];
Chrom=crtbp(NIND,NVAR*PRECI);%创建初始种群
gen=0;%代计数器
ObjV=objfun1(bs2rv(Chrom,FieldD));%计算初始种群个体的目标函数值
whilegenFitnV=ranking(ObjV);%分配适应度值(Assignfitnessvalues)
SelCh=select('sus',Chrom,FitnV,GGAP);%选择
SelCh=recombin('xovsp',SelCh,0.7);%重组
SelCh=mut(SelCh);%变异
ObjVSel=objfun1(bs2rv(SelCh,FieldD));%计算子代目标函数值
[ChromObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);%重插入
gen=gen+1;%代计数器增加
%输出最优解及其对应的20个自变量的十进制值,Y为最优解,I为种群的序号
trace(gen,1)=min(ObjV);%遗传算法性能跟踪
trace(gen,2)=sum(ObjV)/length(ObjV);
end
plot(trace(:
1));holdon;
plot(trace(:
2),'-.');grid;
legend('种群均值的变化','解的变化')
区域描述器的构建采用矩阵复制函数rep建立矩阵FieldD,描述染色体的表示和解释。
一个初始种群被函数crtbp创建,随后产生一个矩阵Chrom,它由NIND个均匀分布长度为NVAR*PRECI的二进制串构成。
使用函数objfun1重新计算目标函数,初始种群中的所有个体的目标函数值由下面程序段计算:
ObjV=objfun1(bs2rv(Chrom,FieldD));
函数bs2rv根据域描述器FieldD转换矩阵Chrom的二进制串为实值,返回一实值表现型的矩阵。
这个bs2rv返回值矩阵通过直接作为目标函数objfun1的输入变量,目标函数结果被返回在矩阵ObjV中。
这个例子中基于排序的适应度分配计算由下面程序段实现:
FitnV=ranking(ObjV);
对这个等级评定算法的缺省设置是选择等差2,使用线性评估,给最适应个体的适应度值为2和最差个体的适应度值为0。
适应度值结果被向量FitnV返回。
使用高级函数选择调用低级函数随机遍历抽样例程sus,程序段为:
SelCh=select(’sus’,Chrom,FitnV,GGAP);
后面的选择中,SelCh包含来自原始染色体的GGAP*NIND个个体,这些个体将使用高级函数recombin进行重组,程序段为:
SelCh=recombin('xovsp',SelCh,0.7);
函数recombin使个体通过SelCh被选择再生产,并使用单点交叉例程函数xovsp,使用交叉概率Px=0.7执行并叉。
个体作为矩阵SelCh的输入被排序,以便使奇数位置的个体与它相邻的个体进行交叉,如果SelCh个体的数量是奇数个,则最后一个个体不进行交叉而返回。
交叉后产生的子代被同一个矩阵SelCh返回,实际使用的交叉例程通过支持使用不同函数名字串传递给recombin而改变。
为了产生一组子代,变异现在使用变异函数mut:
SelCh=mut(SelCh);
子代再次由矩阵SelCh返回,函数调用中变异概率没有被指定,这个缺省值PM=0.7/Lind=0.0017,这里Lind是个体的长度。
子代的目标函数值ObjVSel由程序段计算:
ObjVSel=objfun1(bs2rv(SelCh,FieldD));
因为使用了代沟,子代的数量比当前种群数量要小,完成使用恢复函数reins:
[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);
这里,Chrom和SelCh是矩阵包含原始种群和子代结果。
这两个事件的第一个被使用单个种群和采用基于适应度的恢复,基于适应度的恢复用SelCh中的个体代替Chrom中最不适应的个体。
原始种群目标函数值ObjV随后被要求作为reins的参数,另外新种群的目标函数值不用重新计算整个种群的目标函数而返回,子代的目标值ObjVSel被提供。
Reins返回具有插入子代的新种群Chrom和这个种群的目标函数值ObjV。
这个遗传优化的结果包含在矩阵ObjV中,决策变量的值可能被包含在下面程序段中:
Phen=bs2rv(Chrom,FieldD);
遗传算法的运行结果如下:
(1)图7.7为二维DeJong函数图形。
图7.7二维DeJong函数图形
(2)初始种群中个体的目标函数值分布图如图7.8所示。
图7.8初始种群中个体的目标函数值分布图
(3)经过20次遗传迭代后的结果如图7.9所示。
此时
=2.7412e+005。
图7.9经过20次遗传迭代后的结果
(4)经过100次遗传迭代后的结果如图7.10所示。
此时
=6.9666e+003。
图7.10经过100次遗传迭代后的结果
(5)经过250次遗传迭代后的结果如图7.11所示。
此时
=409.4516。
图7.11经过250次遗传迭代后的结果
(6)经过500次遗传迭代后的结果如图7.12所示。
此时
=4.2550。
图7.12经过500次遗传迭代后的结果
(7)经过250次迭代后种群目标函数均值的变化和最优解的变化如图7.13所示。
=4.2550时20个变量的值见表7.1。
图7.13经过250次迭代后种群目标函数均值的变化和最优解的变化
表7.1目标函数取最优解时20个变量的值
-2.0015
0.1362
-0.0190
-0.0171
-0.0659
-0.1382
-0.0181
-0.0610
-0.0112
0.1274
-0.0142
0.0640
0.0474
0.0122
-0.0962
0.2524
-0.0708
0.2690
-0.1411
-0.0952
7.3多元多峰函数的优化实例
Shubert函数为:
求
。
图7.14为Shubert函数图像。
图7.14Shubert函数图
下面为Shubert函数寻优的MATLAB代码。
functionz=Shubert(x,y)%Shubert函数
z=((1*cos((1+1)*x+1))+(2*cos((2+1)*x+2))+(3*cos((3+1)*x+3))+…
(4*cos((4+1)*x+4))+(5*cos((5+1)*x+5))).*((1*cos((1+1)*y+1))+…
(2*cos((2+1)*y+2))+(3*cos((3+1)*y+3))+(4*cos((4+1)*y+4))+(5*cos((5+1)*y+5)));
[x1,x2]=meshgrid(-10:
.1:
10);
Figure
(1);mesh(x1,x2,Shubert(x1,x2));%画出Shubert函数图像
%定义遗传算法参数
NIND=40;%个体数目(Numberofindividuals)
MAXGEN=50;%最大遗传代数(Maximumnumberofgenerations)
NVAR=2;%变量数目
PRECI=25;%变量的二进制位数(Precisionofvariables)
GGAP=0.9;%代沟(Generationgap)
%建立区域描述器(Buildfielddescriptor)
FieldD=[rep([PRECI],[1,NVAR]);rep([-3;3],[1,NVAR]);rep([1;0;1;1],[1,NVAR])];
Chrom=crtbp(NIND,NVAR*PRECI);%创建初始种群
gen=0;
trace=zeros(MAXGEN,2);%遗传算法性能跟踪初始值
x=bs2rv(Chrom,FieldD);%初始种群十进制转换
ObjV=Shubert(x(:
1),x(:
2));%计算初始种群的目标函数值
whilegenFitnV=ranking(ObjV);%分配适应度值(Assignfitnessvalues)
SelCh=select('sus',Chrom,FitnV,GGAP);%选择
SelCh=recombin('xovsp',SelCh,0.7);%重组
SelCh=mut(SelCh);%变异
x=bs2rv(SelCh,FieldD);%子代十进制转换
ObjVSel=Shubert(x(:
1),x(:
2));%重插入
[ChromObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);
gen=gen+1;
[Y,I]=min(ObjVSel);
Y,bs2rv(Chrom(I,:
),FieldD));%输出最优解及其对应的自变量值
trace(gen,1)=min(ObjV);
trace(gen,2)=sum(ObjV)/length(ObjV);%遗传算法性能跟踪
if(gen==50)%迭代数为50时画出目标函数图
figure
(2);
plot(ObjV);holdon;
plot(ObjV,'b*');grid;
end
end
figure(3);clf;
plot(trace(:
1));holdon;
plot(trace(:
2),'-.');grid;
legend('解的变化','种群均值的变化')
遗传算法的显示结果如下:
(1)初始种群的目标函数值分布如图7.15所示。
图7.15初始种群的目标函数值
(2)经过一次迭代后的目标函数值如图7.16所示。
此时
-1.3900,
-2.9601,
=-44.5224。
图7.16经过一次迭代后的结果
(3)经过10次迭代后的目标函数值如图7.17所示。
此时
-1.3811,
0.8980,
=-186.0739。
图7.17经过10次迭代后的结果
(4)经过50次迭代后的目标函数值如图7.18所示。
此时
-1.4342,
-0.8003,
=-186.7309。
图7.18经过50次迭代后的结果
(5)经过50次迭代后,种群目标函数均值的变化和最优解的变化如图7.19所示。
图7.19经过50次迭代后种群目标函数均值的变化和最优解的变化
7.4收获系统最优控制
收获系统(Harvest)是一个一阶的离散方程,表达式为:
这里
是一个初始的状态条件,
是一个刻度常量,
是状态和非常控制,
是解决问题使用的步骤数。
目标函数为:
这个问题的精确优化解答可由下式确定:
GA工具箱提供一个M文件实现这个函数objharv。
注意,这里是一个最大化问题,工具箱例程实现是最小化问题,这个目标函数objharv,
乘–1即产生为最小化问题。
这个最初条件x(0)设为100,
取为1.1。
这个问题的控制步骤
=20,因此使用的决策变量个数NVAR=20,作为每个控制输入Uk。
决策变量被限制在RANG=[0,200]范围,限制最大的控制输入在任意步为200,这个区域描述器FieldD描述决策变量可以使用矩阵复制函数rep构造。
这个GA的参数可用MATLAB变量来指定。
除了传统GA参数例如代沟(GGAP)和交叉概率(XOVR)外,还有大量与多种群GAS有关的其它参数被定义。
这里INSR=0.9,说明每一代中只有90%产生个体被复制到种群中,SUBPOP=8,说明8个子种群使用迁移概率为MIGR=0.2或20%,每20代(MIGGEN=20)在子种群与当前迁移之间。
每个子种群包含20个个体,NIND=20。
脚本文件使用的函数使用MATLAB串指定。
用串SEL_F、XOV_F、MUT_F、OBJ_F指定分别代表:
选择函数名、重组函数名、变异函数名、目标函数名。
由于使用离散重组函数recdis进行子代的培育,交叉概率没有使用,XOVR=1。
初始种群的创建使用函数crtrp,代计数器gen设置为0。
在由FieldD指定的范围内使用一致性随机挑选个体决策变量组成SUBPOP*NIND个个体。
矩阵Chrom包含了所有子种群并且子种群中所有个体的目标函数值能使用MATLAB的feval命令直接计算。
ObjV=feval(OBJ_F,Chrom);Feval执行函数计算获得第一个输入变量,这里是目标函数名objharv,包含在OBJ_F中,这个函数将被计算并用剩余的参数作为输入调用这个函数。
这里函数调用是:
ObjV=objharv(Chrom);
由于使用实值编码,不需要将染色体转换为表现型表示。
GA现在将进入代循环。
下面为收获系统最优控制的MATLAB代码。
%定义遗传算法参数
NVAR=20;%变量维数
RANGE=[0;200];%变量范围
GGAP=0.8;%代沟(Generationgap)
XOVR=1;%交叉率
MUTR=1/NVAR;%变异率
MAXGEN=500;%最大遗传代数(Maximumnumberofgenerations)
INSR=0.9;%插入率
SUBPOP=8;%子种群数
MIGR=0.2;%迁移率
MIGGEN=20;%在子种群与迁移之间20代
NIND=20;%个体数目(Numberofindividuals)
SEL_F='sus';%