遗传算法,用于求解非线性规划问题,多约束无约束都可以。计算过程:先使用matlab的GA工具箱,如果结果图形不如意,就采用修改源代码。(源代码见于无约束目标函数最大值遗传算法求解策略)
例题前知
给出函数相当于是接下来的适应度句柄函数。
当N=1时,画出图形:
ezplot('-2*pi*exp(-0.2*sqrt(x.^2))-exp(cos(2*pi*x))+2*pi',[-30,30])
当N=2时,画出图形:
[x,y]=meshgrid(-30:30);
z=-2*pi*exp(-0.2*sqrt(0.5*((x).^2+(y).^2)))-exp(0.5*(cos(2*pi*x)+cos(2*pi*y)))+2*pi;
surf(x,y,z)
xlabel('x1'),ylabel('x2'),zlabel('f(x1,x2)')
利用数据游标可以看出最小值为-2.718:
GA工具箱核心函数的用法
(1)函数GA:
[x,fval,reason]=ga(@fitnessfun,nvars,options)
其中x为经过遗传进化以后自变量最佳染色体返回值;fval为最佳染色体的适应度;reason为算法停止的原因;@fitnessfun为适应度句柄函数;nvars为目标函数自变量的个数;options为算法的属性设置,这些属性是通过函数gaoptimset赋予的。
(2)函数gaoptimset:
options=gaoptimset(‘propertyName1’,’propertyValue1’,’propertyName2’,’propertyValue2’……),常用属性如下:
由于算法程序经常重复运行多次才能得到理想结果,鉴于此,可以将前一次运行得到的最后种群作为下一次运行的初始种群,如此操作会得到更好的结果。例如:
[x,fval,reason,output,final_pop]=ga(@fitnessfcn,nvars)
最后一个输出变量final_pop返回的就是本次运行得到的最后种群,再将final_pop作为函数ga的初始种群,语法格式为:
options=gaoptimset(‘InitialPopulation’,final_pop);
[x,fval,reason,output,final_pop2]=ga(@fitnessfcn,nvars,options);
函数gaoptimset常用的11种属性:
适应度函数设计
function f=fitnessfcn(x)
f=f(x);
如果有约束条件(包括自变量的取值范围),对于求解函数的最小值问题,可以使用如下语法格式:
function f=fitnessfcn(x)
if(x<=-1|x>3)
f=inf;
else
f=f(x);
end
如果有约束条件(包括自变量的取值范围),对于求解函数的最大值问题,可以使用如下语法格式:
function f=fitnessfcn(x)
if(x<=-1|x>3)
f=inf;
else
f=-f(x); %这里是-f(x)而不是f=f(x)
end
利用GA工具箱求解例题
接下来,利用GA工具箱求解上述函数N=10时的最小值。
lbw.m:
function f=lbw(x)
if(x(1)>30|x(1)<-30|x(2)>30|x(2)<-30|x(3)>30|x(3)<-30|x(4)>30|x(4)<-30|x(5)>30|x(5)<-30|x(6)>30|x(6)<-30|x(7)>30|x(7)<-30|x(8)>30|x(8)<-30|x(9)>30|x(9)<-30|x(10)>30|x(10)<-30)
f=300;
else
f=-2*pi*exp(-0.2*sqrt(1/10*((x(1)).^2+(x(2)).^2+(x(3)).^2+(x(4)).^2+(x(5)).^2+(x(6)).^2+(x(7)).^2+(x(8)).^2+(x(9)).^2+(x(10)).^2)))-exp(1/10*(cos(2*pi*x(1))+cos(2*pi*x(2))+cos(2*pi*x(3))+cos(2*pi*x(4))+cos(2*pi*x(5))+cos(2*pi*x(6))+cos(2*pi*x(7))+cos(2*pi*x(8))+cos(2*pi*x(9))+cos(2*pi*x(10))))+2*pi;
end
command:
options=gaoptimset('Generations',800,'StallGenLimit',300,'PlotFcns',@gaplotbestf);
[x,f]=ga(@lbw,10,options)
无约束目标函数最大值遗传算法求解策略
先画出图看看:
ezplot('200*exp(-0.05*x)*sin(x)',[-2,2])
使用数据游标,得到最高点的坐标:
利用遗传算法来计算,这里没有使用现成的GA工具库,而是使用改写源代码。
GA501.m:
%主程序:用遗传算法求解y=200*exp(-0.05*x).*sin(x)在[-2 2]区间上的最大值
clc;
clear all;
close all;
global BitLength
global boundsbegin
global boundsend
bounds=[-2 2];%一维自变量的取值范围
precision=0.0001; %运算精度
boundsbegin=bounds(:,1);%=-2
boundsend=bounds(:,2);%=2
%计算如果满足求解精度至少需要多长的染色体
BitLength=ceil(log2((boundsend-boundsbegin)' ./ precision));
popsize=500; %初始种群大小
Generationnmax=12; %最大代数
pcrossover=0.90; %交配概率
pmutation=0.09; %变异概率
%产生初始种群
population=round(rand(popsize,BitLength));
%计算适应度,返回适应度Fitvalue和累积概率cumsump
[Fitvalue,cumsump]=fitnessfun(population);
Generation=1;
while Generation<Generationnmax+1
for j=1:2:popsize
%选择操作
seln=selection(population,cumsump);
%交叉操作
scro=crossover(population,seln,pcrossover);
scnew(j,:)=scro(1,:);
scnew(j+1,:)=scro(2,:);
%变异操作
smnew(j,:)=mutation(scnew(j,:),pmutation);
smnew(j+1,:)=mutation(scnew(j+1,:),pmutation);
end
population=smnew; %产生了新的种群
%计算新种群的适应度
[Fitvalue,cumsump]=fitnessfun(population);
%记录当前代最好的适应度和平均适应度
[fmax,nmax]=max(Fitvalue);
fmean=mean(Fitvalue);
ymax(Generation)=fmax;
ymean(Generation)=fmean;
%记录当前代的最佳染色体个体
x=transform2to10(population(nmax,:));
%自变量取值范围是[-2 2],需要把经过遗传运算的最佳染色体整合到[-2 2]区间
xx=boundsbegin+x*(boundsend-boundsbegin)/(power((boundsend),BitLength)-1);
xmax(Generation)=xx;
Generation=Generation+1
end
Generation=Generation-1;
Bestpopulation=xx
Besttargetfunvalue=targetfun(xx)
%绘制经过遗传运算后的适应度曲线。一般地,如果进化过程中种群的平均适应度与最大适
%应度在曲线上有相互趋同的形态,表示算法收敛进行得很顺利,没有出现震荡;在这种前
%提下,最大适应度个体连续若干代都没有发生进化表明种群已经成熟。
figure(1);
hand1=plot(1:Generation,ymax);
set(hand1,'linestyle','-','linewidth',1.8,'marker','*','markersize',6)
hold on;
hand2=plot(1:Generation,ymean);
set(hand2,'color','r','linestyle','-','linewidth',1.8,...
'marker','h','markersize',6)
xlabel('进化代数');ylabel('最大/平均适应度');xlim([1 Generationnmax]);
legend('最大适应度','平均适应度');
box off;hold off;
crosscover.m:
%子程序:新种群交叉操作,函数名称存储为crossover.m
function scro=crossover(population,seln,pc);
BitLength=size(population,2);
pcc=IfCroIfMut(pc); %根据交叉概率决定是否进行交叉操作,1则是,0则否
if pcc==1
chb=round(rand*(BitLength-2))+1; %在[1,BitLength-1]范围内随机产生一个交叉位
scro(1,:)=[population(seln(1),1:chb) population(seln(2),chb+1:BitLength)];
scro(2,:)=[population(seln(2),1:chb) population(seln(1),chb+1:BitLength)];
else
scro(1,:)=population(seln(1),:);
scro(2,:)=population(seln(2),:);
end
fitnessfun.m:
%子程序:计算适应度函数, 函数名称存储为fitnessfun
function [Fitvalue,cumsump]=fitnessfun(population);
global BitLength
global boundsbegin
global boundsend
popsize=size(population,1); %有popsize个个体
for i=1:popsize
x=transform2to10(population(i,:)); %将二进制转换为十进制
%转化为[-2,2]区间的实数
xx=boundsbegin+x*(boundsend-boundsbegin)/(power((boundsend),BitLength)-1);
Fitvalue(i)=targetfun(xx); %计算函数值,即适应度
end
%给适应度函数加上一个大小合理的数以便保证种群适应值为正数
Fitvalue=Fitvalue'+230;
%计算选择概率
fsum=sum(Fitvalue);
Pperpopulation=Fitvalue/fsum;
%计算累积概率
cumsump(1)=Pperpopulation(1);
for i=2:popsize
cumsump(i)=cumsump(i-1)+Pperpopulation(i);
end
cumsump=cumsump';
IfCroIfMut.m:
%子程序:判断遗传运算是否需要进行交叉或变异, 函数名称存储为IfCroIfMut.m
function pcc=IfCroIfMut(mutORcro);
test(1:100)=0;
l=round(100*mutORcro);
test(1:l)=1;
n=round(rand*99)+1;
pcc=test(n);
mutation.m:
%子程序:新种群变异操作,函数名称存储为mutation.m
function snnew=mutation(snew,pmutation);
BitLength=size(snew,2);
snnew=snew;
pmm=IfCroIfMut(pmutation); %根据变异概率决定是否进行变异操作,1则是,0则否
if pmm==1
chb=round(rand*(BitLength-1))+1; %在[1,BitLength]范围内随机产生一个变异位
snnew(chb)=abs(snew(chb)-1);
end
selection.m:
%子程序:新种群选择操作, 函数名称存储为selection.m
function seln=selection(population,cumsump);
%从种群中选择两个个体
for i=1:2
r=rand; %产生一个随机数
prand=cumsump-r;
j=1;
while prand(j)<0
j=j+1;
end
seln(i)=j; %选中个体的序号
end
targetfun.m:
%子程序:对于优化最大值或极大值函数问题,目标函数可以作为适应度函数
%函数名称存储为targetfun.m
function y=targetfun(x); %目标函数
y=200*exp(-0.05*x).*sin(x);
transform2to10.m:
%子程序:将2进制数转换为10进制数,函数名称存储为transform2to10.m
function x=transform2to10(Population);
BitLength=size(Population,2);
x=Population(BitLength);
for i=1:BitLength-1
x=x+Population(BitLength-i)*power(2,i);
end
多约束非线性规划问题的求解
先画出目标函数的三维图像:
[x,y]=meshgrid(-5:10);
z=exp(x).*(4*x^2+2*y^2+4*x*y+2*y+1);
surf(x,y,z)
xlabel('x1'),ylabel('x2'),zlabel('f(x1,x2)')
ch14_2f.m:
%子函数:适应度函数同时也是目标函数,函数存储名称为ch14_2f.m
function f=ch14_2f(x)
g1=1.5+x(1)*x(2)-x(1)-x(2);
g2=-x(1)*x(2);
if(g1>0|g2>10)
f=100;
else
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
end
GA502.m:
%主程序:本程序采用遗传算法接力进化,
%将上次进化结束后得到的最终种群作为下次输入的初始种群
clc;
close all;
clear all;
%进化的代数
T=100;
optionsOrigin=gaoptimset('Generations',T/2);
[x,fval,reason,output,finnal_pop]=ga(@ch14_2f,2,optionsOrigin);
%进行第二次接力进化
options1=gaoptimset('Generations',T/2,'InitialPopulation',finnal_pop,...
'PlotFcns',@gaplotbestf);
[x,fval,reason,output,finnal_pop]=ga(@ch14_2f,2,options1);
Bestx=x
BestFval=fval
这里采用了两次迭代进化,在前一次进化上再进化一次。
但是这个结果有点奇怪~