遗传算法就是在一定的自变量有限的取值范围内,随机取若干个个体,每个个体相当于自变量范围内的一个取值,若干个个体共同组成一个种群,个体对于环境的适应能力体现为该个体对应的因变量,不同的个体得到的结果不同,对于结果较好的个体,其下一代在种群中的占比更高,对于结果不好的个体,其下一代在种群中的占比会更少,简单来说,就是好的个体被保留,坏的个体被淘汰。经过不断的更新换代,最后结果会不断逼近最优的结果。
同时参考自然界种群的进化过程,除了物竞天择这一思想外,还有基因的突变和基因的交叉互换来探索更多的优秀的形状,从而提高种群的适应性。所以在遗传算法中也参考了这种思想,在遗传算法中,一个个体由能代表该个体属性的多个二进制编码构成。具有相同属性的个体二进制编码相同。在每一轮的种群更新中,都会进行一次交叉互换和一次基因变异。
在进行交叉互换时,根据交叉互换的概率,对随机的个体与另外一个随机的个体进行一次基因的交叉互换,即随机取基因的某一处断开,将断开后的基因片段进行交换,再组合成新的片段。
在进行变异的时候,根据变异的概率,对随机的个体的基因中的随机某个基因进行取反,得到新的基因。
那么一个遗传算法的过程就是以上的流程。下面具体谈谈怎么实现。
一个个体的性状的表达,是根据个体的基因属性所对应的环境的适应度。以求某函数的最大值为例,个体的不同基因属性代表了不同的自变量取值,该自变量对应的函数结果就是该个体在整个环境中的性状的表达。环境对该性状的选择就是通过数学方法来将结果与最理想的结果进行一个评估,例如当求最大值时,就可以让函数结果大的自变量取值尽可能的保留,结果小的取值尽可能抛弃,这就是自然选择的过程。
这里在进行自然选择的时候,用到了轮盘赌的数学模型。轮盘赌是一个有着不同面积区块的圆形转盘,每次转动停到的位置都是随机的,那么最后的结果就是根据转盘中的区块的面积占比来决定每次转动出现的各个结果的概率。这是一个很简单的数学模型。在自然选择过程中,就可以通过这样的方式来淘汰形状较差的个体。具体方法就是,每个个体都有一个数值作为其对环境的适应能力的体现,所有个体的适应力数值加在一起就是一整个转盘的面积1,各个个体的适应力比上所有的适应力之和就是该个体区域面积占整个转盘面积的大小。数值越大,对应的就是转盘上的面积越大,随机取一个数,落在该区域的概率就越大,因此适应性强的将更有概率出现在下一代。这就是轮盘赌。
以上一个遗传算法的基本原理过程就搞清楚了。
然后就是代码部分
主程序:
clc;clear;close all;
%%
% f= @(a,b)(1/2*a .* sin(a) .* cos(1/2 * a) ).*(1/2*b .* sin(b) .* cos(1/2 * b) ); % 函数表达式
f= @(a,b)(10*sin(a)+7*abs(a-5)+10).*(10*sin(b)+7*abs(b-5)+10); % 函数表达式
figure(1);
[x0_1, x0_2]=meshgrid(0:.2:20);
y0=f(x0_1,x0_2);
mesh(x0_1, x0_2, y0);
hold on;
N = 50; % 初始种群个数
chromosome = 2; % 空间维数
ger = 50; % 最大迭代次数
chromlength =[16,16];
limit = [0, 20;0, 20];
pc = 0.3; % 交叉概率
pm = 0.008; % 变异概率
all_best=0;
all_x=0;
%%
% for i=1:chromosome
% eval(['POP.chrom',num2str(i),'=round(rand(N,chromlength(i)));']);
% for j=1:N
% eval(['POP.x',num2str(i),'(j)=(binary2decimal(POP.chrom1(j,:)))/(2^chromlength(i)-1)*(limit(i,2)-limit(i,1)) ;'])
% % POP.x'i'(j)=(binary2decimal(POP.chrom1(j,:)))/(2^chromlength(i)-1)*(limit(i,2)-limit(i,1)) ;
% end
% end
POP.chrom1=round(rand(N,chromlength(1)));
POP.chrom2=round(rand(N,chromlength(2)));
%%
% 开始迭代
for n=1:ger
%%
% 交叉互换
for i = 1:2:N-1
if(rand<pc)
cpoint = round(rand*chromlength);
POP.NEWchrom1(i,:) = [POP.chrom1(i,1:cpoint(1)),POP.chrom1(i+1,cpoint(1)+1:chromlength(1))];
POP.NEWchrom1(i+1,:) = [POP.chrom1(i+1,1:cpoint(1)),POP.chrom1(i,cpoint(1)+1:chromlength(1))];
POP.NEWchrom2(i,:) = [POP.chrom2(i,1:cpoint(2)),POP.chrom2(i+1,cpoint(2)+1:chromlength(2))];
POP.NEWchrom2(i+1,:) = [POP.chrom2(i+1,1:cpoint(2)),POP.chrom2(i,cpoint(2)+1:chromlength(2))];
else
POP.NEWchrom1(i,:) = POP.chrom1(i,:);
POP.NEWchrom1(i+1,:) = POP.chrom1(i+1,:);
POP.NEWchrom2(i,:) = POP.chrom2(i,:);
POP.NEWchrom2(i+1,:) = POP.chrom2(i+1,:);
end
end
% 根据交叉互换的结果 更新种群的基因
POP.chrom1=POP.NEWchrom1;
POP.chrom2=POP.NEWchrom2;
%%
% 基因变异并更新种群
for i=1:N
for j1=1:chromlength(1)
if(rand<pm)
POP.chrom1(i,j1)=~POP.chrom1(i,j1);
end
end
for j2=1:chromlength(2)
if(rand<pm)
POP.chrom2(i,j2)=~POP.chrom2(i,j2);
end
end
end
%%
% 形状表达与选择
% 将基因(二进制编码)转化为自变量的取值(10进制的数)
for i=1:N
POP.x1(i)=(binary2decimal(POP.chrom1(i,:)))/(2^chromlength(1)-1)*(limit(1,2)-limit(1,1)) ;
end
for i=1:N
POP.x2(i)=(binary2decimal(POP.chrom2(i,:)))/(2^chromlength(2)-1)*(limit(2,2)-limit(2,1)) ;
end
% 根据自变量的取值得到函数的输出
for i=1:N
POP.y(i)=f(POP.x1(i),POP.x2(i));
end
%将输出的结果单位化,转化为0-1之间的数值长度(相当于轮盘赌的各个区域的面积)
a=min(POP.y);
b=sum(POP.y)+N*(-a);
for i=1:N
POP.adapt(i)=(POP.y(i)-a)/b;
end
%数值长度转换为0-1之间的区间的节点(相当于把面积转化为了轮盘赌上各个区域的边界线)
POP.NWEadapt(1)=POP.adapt(1);
for i=2:N
POP.NWEadapt(i)=POP.adapt(i)+POP.NWEadapt(i-1);
end
%进行轮盘赌,任取一个随机数cs,求这个随机数在轮盘赌中的位置区域
%到达某个区域,就代表下一个种群在第i个个体就拥有该区域所表示的基因,从而得到新种群
for i=1:N
cs=rand;
[a,b]=POP_erfen(1,N,POP.NWEadapt,cs);
POP.NEWchrom1(i,:)=POP.chrom1(b,:);
POP.NEWchrom2(i,:)=POP.chrom2(b,:);
end
%更新种群
POP.chrom1=POP.NEWchrom1;
POP.chrom2=POP.NEWchrom2;
%%
% 获得最优位置
[this_x this_best]=best(POP)
if all_best<this_best
all_best=this_best;
all_x=this_x;
end
cla;
mesh(x0_1, x0_2, y0);
plot3(POP.x1,POP.x2,POP.y, 'ro');title('状态位置变化');
pause(0.1);
end
figure();
mesh(x0_1, x0_2, y0);
hold on;
plot3(all_x(1),all_x(2),all_best, 'ro');title('最优位置图');
disp(['最大值:',num2str(all_best)]);
disp(['变量取值:',num2str(all_x)]);
另外在代码中用了两个自己写的函数,一个是二分法求随机数在轮盘赌中的位置,一个是求种群中的最优个体:
二分法:
function [Nmin,Nmax]=POP_erfen(min,max,fun,tal)
if tal<fun(1)
Nmin=0;
Nmax=1;
else
a=min;
b=max;
if (b-a)<=1
Nmin=a;
Nmax=b;
else
c=round((a+b)/2);
if tal>=fun(c)
a=c;
b=b;
[Nmin,Nmax]=POP_erfen(a,b,fun,tal);
else
a=a;
b=c;
[Nmin,Nmax]=POP_erfen(a,b,fun,tal);
end
end
end
end
最优解:
function [bestindividual bestfit] = best(POP)
[px,py] = size(POP.y);
bestindividual = [POP.x1(1) , POP.x2(1)];
bestfit = POP.y(1);
for i = 2:px
if POP.y(i)>bestfit
bestindividual = [POP.x1(i) , POP.x2(i)];
bestfit = POP.y(i);
end
end
二进制转十进制函数:
%二进制转化成十进制函数
%输入变量:
%二进制种群
%输出变量
%十进制数值
function pop2 = binary2decimal(pop)
[px,py]=size(pop);
for i = 1:py
pop1(:,i) = 2.^(py-i).*pop(:,i);
end
%sum(.,2)对行求和,得到列向量
pop2 = sum(pop1,2);
(后面打算把整个过程用GIF图表现出来,未完待续)