司守奎《数学建模算法与应用》-遗传算法-保有怀疑
遗传算法属于是智能算法里比较好懂的一类算法,在司守奎老师的《数学建模算法与应用》中,遗传算法被应用于求解TSP旅行商问题,展现出优良的求解性质,仅仅用秒到毫秒级别的时间求解了含有101个地点的TSP问题解决方案,求解问题的路线长的平均值明显优于本书上一章讲解的模拟退火算法。
但在研读遗传算法的过程中,我发现了几个很疑惑的事情。
染色体的冲突检测
%完整代码及附件见文末GA_pre.txt
F=2+floor(100*rand(1)); %产生交叉操作的地址,floor代表取整
temp=A(c(i),[F:102]); %中间变量的保存值,这应该是单点交叉,交叉后半段
A(c(i),[F:102])=A(c(i+1),[F:102]); %交叉操作
A(c(i+1),F:102)=temp; %交换染色体段,需要temp中转
司老师这里直接用亲本的染色体交叉段给到子代,却没有进行冲突检测,由于亲本的基因序列一般为乱序,这意味着子代直接接受交叉段后,可能会与本身未交叉段有相同的基因,同时也注定会遗漏别的基因。
通过length(unique(G(i,1:102))),求子代、父代、变异代共计100+个个体,的非重复点,发现仅仅50+个个体有完整的长度102(尾部同首部一个点,故为101+1);具体检查的代码内容写入GA_pre文件了。
但是父代已知有50个,长度是完整的,这意味着有大量夭折的子代没有利用,遗传算法没有发挥真正的作用。
故此,我已经采用setdiff函数求得差集,并按乱序把非交叉段给到子代了,让子代也能拥有完整路径;代码文本见附件my_trial.txt
案例中遗传算法徒有其表
本着探究遗传算法的能效的初心,于是我将“改良圈算法”的结果和“在这之后的再进行遗传操作”的结果进行对比,看看遗传算法在这基础有多大优化,结果令我大跌眼镜。
遗传算法的代码占了大约一半篇幅,却经过对比,我发现遗传操作前和操作后的解的路径和时间都没变,也就是在改良圈后,遗传操作本身没有发挥优化作用。具体检验操作紧挨在在my_trial.txt的改良圈算法后,如下:
%多列注释按ctrl+R,多列取消注释按ctrl+T
[SG,ind1]=sort(J,2);
num=size(J,1); long=zeros(1,num); %路径长度的初始值
for j=1:num%j就是对于每个染色体
for i=1:101
long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度
end
end
[slong]=sort(long); %对路径长度按照从小到大排
path_pre=ind1(1,:);flong_pre=slong(1)%改良圈的最优解用_pre
讽刺的是,当我把遗传操作全部注释掉后,算法运算速度提高到不到0.1秒的级别;而若取消注释,多次运行,发现遗传前后遗传后的路径path完全相同,具体检验代码写入了my_trail.txt的后边,具体如下:
if path==path_pre
disp("path equals to GA before")
else
disp("path different to GA before")
end
改良后的遗传算法颇有鸡肋
本书的下一节对遗传算法进行了改进,采用了什么混沌序列的交叉操作,经过检测,交叉冲突的个体确实是很少了,但耗时也变多了,差不多是3秒左右,相对于原来的遗传算法翻了近10倍
但是在改良圈之后,它的优化作用似乎也没怎么显现出来。多次运行,发现path_pre和path确实不同了,但是flong_pre和flong还是一样的,也就是路线不同,耗时还是一样的。属实给我看的一脸懵逼了。代码见附件GA.txt。
写在最后总的来说,我属于是一脸懵逼的,作为智能算法小白,仅仅通过个人的一点验证,并不能有很好的把握确定自己的对错,因此记录与此。
GA_pre.txt

tic clc,clear, close all sj0=load('data12_1.txt'); x=sj0(:,1:2:8); x=x(:); y=sj0(:,2:2:8); y=y(:); sj=[x y]; d1=[70,40]; xy=[d1;sj;d1]; sj=xy*pi/180; %单位化成弧度 d=zeros(102); %距离矩阵d的初始值 for i=1:101 for j=i+1:102 d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*... cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2))); end end d=d+d'; w=50; g=100; %w为种群的个数,g为进化的代数 for k=1:w %通过改良圈算法选取初始种群 c=randperm(100); %产生1,...,100的一个随机排列 c1=[1,c+1,102]; %生成初始解 for t=1:102 %该层循环是修改圈 %至于为什么改良上限是102次,作者应该有他的理论依据 flag=0; %修改圈退出标志 for m=1:100 for n=m+2:101 if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<... d(c1(m),c1(m+1))+d(c1(n),c1(n+1)) c1(m+1:n)=c1(n:-1:m+1); flag=1; %若产生改良一次,flag就为1。修改圈 end end end %对于单个随机圈,可以多次改良(有待思考),一旦无法改良,flag=0,则可以退出,进行另一圈的改良 if flag==0 %J(k,c)=1:102是指第k行,第c列进行赋值[1,2,3...102],由于c1是乱序的,刚好可以把1:102乱序为自己本身 J(k,c1)=1:102; break %记录下较好的解并退出当前层循环 end end end J(:,1)=0; J=J/102; %把整数序列转换成[0,1]区间上的实数,即转换成染色体编码 for k=1:g %该层循环进行交配的操作 A=J; %交配产生子代A的初始染色体 c=randperm(w); %产生下面交叉操作的染色体对 for i=1:2:w %跨度为2,是为了每两个凑一对亲本 F=2+floor(100*rand(1)); %产生交叉操作的地址,floor代表取整 temp=A(c(i),[F:102]); %中间变量的保存值,这应该是单点交叉,交叉后半段 %一次交配产生两个子代,先揪出配对亲本,再对交叉部分进行交换 A(c(i),[F:102])=A(c(i+1),[F:102]); %交叉操作 A(c(i+1),F:102)=temp; %交换染色体段,需要temp中转 % 这里没有做冲突检测。导致有重复遍历 % 发现G中通常只有50+个染色体是可行解,这其中显然有50个是最初改良圈的50个父本的,效率很低 end by=[]; %为了防止下面产生空地址,这里先初始化 while isempty(by) by=find(rand(1,w)<0.1); %产生w个对象,选择<0.1的进行变异操作(即变异概率十分之一) end B=A(by,:); %产生变异操作的初始染色体(初始子代) for j=1:length(by) bw=sort(2+floor(100*rand(1,3))); %产生变异操作的3个地址 B(j,:)=B(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]); %分成4段,交换中间两段 end G=[J;A;B]; %父代和子代种群合在一起,总计50+50+length(by) [SG,ind1]=sort(G,2); %把染色体翻译成1,...,102的序列ind1,ind1是单纯的index,是索引 num=size(G,1); long=zeros(1,num); %路径长度的初始值 for j=1:num%j就是对于每个染色体 for i=1:101 long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度 end end [slong]=sort(long); %对路径长度按照从小到大排 J=G(1:w,:); %精选前w个较短的路径对应的染色体 end path=ind1(1,:), flong=slong(1) %解的路径及路径长度 %但是多次运行发现 解的长度总为102,某种程度上说明:不可行解的路径长度一般是大于可行解的 toc count=0; for i=1:length(G) if length(unique(G(i,1:102)))==102 count=count+1; end end fprintf("%d",count) xx=xy(path,1);yy=xy(path,2); plot(xx,yy,'-o') %画出路径
my_trial.txt

1 tic 2 clc,clear, close all 3 sj0=load('data12_1.txt'); 4 x=sj0(:,1:2:8); x=x(:); 5 y=sj0(:,2:2:8); y=y(:); 6 sj=[x y]; d1=[70,40]; 7 xy=[d1;sj;d1]; sj=xy*pi/180; %单位化成弧度 8 d=zeros(102); %距离矩阵d的初始值 9 for i=1:101 10 for j=i+1:102 11 d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*... 12 cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2))); 13 end 14 end 15 d=d+d'; w=50; g=100; %w为种群的个数,g为进化的代数 16 for k=1:w %通过改良圈算法选取初始种群 17 c=randperm(100); %产生1,...,100的一个随机排列 18 c1=[1,c+1,102]; %生成初始解 19 for t=1:102 %该层循环是修改圈 %至于为什么改良上限是102次,作者应该有他的理论依据 20 flag=0; %修改圈退出标志 21 for m=1:100 22 for n=m+2:101 23 if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<... 24 d(c1(m),c1(m+1))+d(c1(n),c1(n+1)) 25 c1(m+1:n)=c1(n:-1:m+1); flag=1; %若产生改良一次,flag就为1。修改圈 26 end 27 end 28 end 29 %对于单个随机圈,可以多次改良(有待思考),一旦无法改良,flag=0,则可以退出,进行另一圈的改良 30 if flag==0 31 %J(k,c)=1:102是指第k行,第c列进行赋值[1,2,3...102],由于c1是乱序的,刚好可以把1:102乱序为自己本身 32 J(k,c1)=1:102; break %记录下较好的解并退出当前层循环 33 end 34 end 35 end 36 37 [SG,ind1]=sort(J,2); 38 num=size(J,1); long=zeros(1,num); %路径长度的初始值 39 for j=1:num%j就是对于每个染色体 40 for i=1:101 41 long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度 42 end 43 end 44 [slong]=sort(long); %对路径长度按照从小到大排 45 path_pre=ind1(1,:);flong_pre=slong(1) 46 47 48 J(:,1)=0; J=J/102; %把整数序列转换成[0,1]区间上的实数,即转换成染色体编码 49 for k=1:g %该层循环进行交配的操作 50 %A=J; %交配产生子代A的初始染色体 51 c=randperm(w); %产生下面交叉操作的染色体对 52 for i=1:2:w %跨度为2,是为了每两个凑一对亲本 53 F=2+floor(100*rand(1)); %产生交叉操作的地址,floor代表取整 54 55 %一次交配产生两个子代,先揪出配对亲本,再对交叉部分进行交换 56 57 A(c(i),F:102)=J(c(i+1),F:102); %将母本交叉段给孩子1后面 58 muji=J(c(i),1:102); 59 ziji=J(c(i+1),F:102); 60 chaji=setdiff(muji,ziji); 61 [~,loc]=ismember(chaji,muji); 62 A(c(i),1:F-1)=muji(sort(loc));%将父本中没有母本交叉段的基因输入孩子1前面 63 A(c(i+1),F:102)=J(c(i),F:102);%同理,构造孩子2 64 muji=J(c(i+1),1:102); 65 ziji=J(c(i),F:102); 66 chaji=setdiff(muji,ziji); 67 [~,loc]=ismember(chaji,muji); 68 A(c(i+1),1:F-1)=muji(sort(loc)); 69 70 71 end 72 by=[]; %为了防止下面产生空地址,这里先初始化 73 while isempty(by) 74 by=find(rand(1,w)<0.1); %产生w个对象,选择<0.1的进行变异操作(即变异概率十分之一) 75 end 76 B=A(by,:); %产生变异操作的初始染色体(初始子代) 77 78 for j=1:length(by) 79 bw=sort(2+floor(100*rand(1,3))); %产生变异操作的3个地址 80 B(j,:)=B(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]); %分成4段,交换中间两段 81 end 82 G=[J;A;B]; %父代和子代种群合在一起,总计50+50+length(by) 83 [SG,ind1]=sort(G,2); %把染色体翻译成1,...,102的序列ind1,ind1是单纯的index,是索引 84 num=size(G,1); long=zeros(1,num); %路径长度的初始值 85 for j=1:num%j就是对于每个染色体 86 for i=1:101 87 long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度 88 end 89 end 90 [slong]=sort(long); %对路径长度按照从小到大排 91 J=G(1:w,:); %精选前w个较短的路径对应的染色体 92 end 93 path=ind1(1,:);flong=slong(1) %解的路径及路径长度 94 % %但是多次运行发现 解的长度总为102,某种程度上说明:不可行解的路径长度一般是大于可行解的 95 if path==path_pre 96 disp("path equals to GA before") 97 else 98 disp("path different to GA before") 99 end 100 toc 101 xx=xy(path_pre,1);yy=xy(path_pre,2); 102 plot(xx,yy,'-o') %画出路径
GA.txt

1 tic %计时开始 2 clc,clear, close all 3 sj0=load('data12_1.txt'); 4 x=sj0(:,1:2:8); x=x(:); 5 y=sj0(:,2:2:8); y=y(:); 6 sj=[x y]; d1=[70,40]; 7 xy=[d1;sj;d1]; sj=xy*pi/180; %单位化成弧度 8 d=zeros(102); %距离矩阵d的初始值 9 for i=1:101 10 for j=i+1:102 11 d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*... 12 cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2))); 13 end 14 end 15 d=d+d'; w=500; g=1000; %w为种群的个数,g为进化的代数 16 rand('state',sum(clock)); %初始化随机数发生器 17 for k=1:w %通过改良圈算法选取初始种群 18 c=randperm(100); %产生1,...,100的一个全排列 19 c1=[1,c+1,102]; %生成初始解 20 for t=1:102 %该层循环是修改圈 21 flag=0; %修改圈退出标志 22 for m=1:100 23 for n=m+2:101 24 if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<... 25 d(c1(m),c1(m+1))+d(c1(n),c1(n+1)) 26 c1(m+1:n)=c1(n:-1:m+1); flag=1; %修改圈 27 end 28 end 29 end 30 if flag==0 31 J(k,c1)=1:102; break %记录下较好的解并退出当前层循环 32 end 33 end 34 end 35 36 [SG,ind1]=sort(J,2); 37 num=size(J,1); long=zeros(1,num); %路径长度的初始值 38 for j=1:num%j就是对于每个染色体 39 for i=1:101 40 long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度 41 end 42 end 43 [slong]=sort(long); %对路径长度按照从小到大排 44 path_pre=ind1(1,:);flong_pre=slong(1) 45 46 47 J(:,1)=0; J=J/102; %把整数序列转换成[0,1]区间上的实数,即转换成染色体编码 48 for k=1:g %该层循环进行遗传算法的操作 49 A=J; %交配产生子代A的初始染色体 50 for i=1:2:w 51 ch1(1)=rand; %混沌序列的初始值 52 for j=2:50 53 ch1(j)=4*ch1(j-1)*(1-ch1(j-1)); %产生混沌序列 54 end 55 ch1=2+floor(100*ch1); %产生交叉操作的地址 56 temp=A(i,ch1); %中间变量的保存值 57 A(i,ch1)=A(i+1,ch1); %交叉操作 58 A(i+1,ch1)=temp; 59 end 60 by=[]; %为了防止下面产生空地址,这里先初始化 61 while ~length(by) 62 by=find(rand(1,w)<0.1); %产生变异操作的地址 63 end 64 num1=length(by); B=J(by,:); %产生变异操作的初始染色体 65 ch2=rand; %产生混沌序列的初始值 66 for t=2:2*num1 67 ch2(t)=4*ch2(t-1)*(1-ch2(t-1)); %产生混沌序列 68 end 69 for j=1:num1 70 bw=sort(2+floor(100*rand(1,2))); %产生变异操作的2个地址 71 B(j,bw)=ch2([j,j+1]); %bw处的两个基因发生了变异 72 end 73 G=[J;A;B]; %父代和子代种群合在一起 74 [SG,ind1]=sort(G,2); %把染色体翻译成1,...,102的序列ind1 75 num2=size(G,1); long=zeros(1,num2); %路径长度的初始值 76 for j=1:num2 77 for i=1:101 78 long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度 79 end 80 end 81 [slong,ind2]=sort(long); %对路径长度按照从小到大排序 82 J=G(ind2(1:w),:); %精选前w个较短的路径对应的染色体 83 end 84 path=ind1(ind2(1),:);flong=slong(1) %解的路径及路径长度 85 toc %计时结束 86 if path==path_pre 87 disp("path equals to GA before") 88 else 89 disp("path different to GA before") 90 end 91 92 xx=xy(path,1);yy=xy(path,2); 93 plot(xx,yy,'-o') %画出路径
data.txt

1 53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 30.3313 6.9348 2 56.5432 21.4188 10.8198 16.2529 22.7891 23.1045 10.1584 12.4819 3 20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 31.4847 8.9640 4 26.2418 18.1760 44.0356 13.5401 28.9836 25.9879 38.4722 20.1731 5 28.2694 29.0011 32.1910 5.8699 36.4863 29.7284 0.9718 28.1477 6 8.9586 24.6635 16.5618 23.6143 10.5597 15.1178 50.2111 10.2944 7 8.1519 9.5325 22.1075 18.5569 0.1215 18.8726 48.2077 16.8889 8 31.9499 17.6309 0.7732 0.4656 47.4134 23.7783 41.8671 3.5667 9 43.5474 3.9061 53.3524 26.7256 30.8165 13.4595 27.7133 5.0706 10 23.9222 7.6306 51.9612 22.8511 12.7938 15.7307 4.9568 8.3669 11 21.5051 24.0909 15.2548 27.2111 6.2070 5.1442 49.2430 16.7044 12 17.1168 20.0354 34.1688 22.7571 9.4402 3.9200 11.5812 14.5677 13 52.1181 0.4088 9.5559 11.4219 24.4509 6.5634 26.7213 28.5667 14 37.5848 16.8474 35.6619 9.9333 24.4654 3.1644 0.7775 6.9576 15 14.4703 13.6368 19.8660 15.1224 3.1616 4.2428 18.5245 14.3598 16 58.6849 27.1485 39.5168 16.9371 56.5089 13.7090 52.5211 15.7957 17 38.4300 8.4648 51.8181 23.0159 8.9983 23.6440 50.1156 23.7816 18 13.7909 1.9510 34.0574 23.3960 23.0624 8.4319 19.9857 5.7902 19 40.8801 14.2978 58.8289 14.5229 18.6635 6.7436 52.8423 27.2880 20 39.9494 29.5114 47.5099 24.0664 10.1121 27.2662 28.7812 27.6659 21 8.0831 27.6705 9.1556 14.1304 53.7989 0.2199 33.6490 0.3980 22 1.3496 16.8359 49.9816 6.0828 19.3635 17.6622 36.9545 23.0265 23 15.7320 19.5697 11.5118 17.3884 44.0398 16.2635 39.7139 28.4203 24 6.9909 23.1804 38.3392 19.9950 24.6543 19.6057 36.9980 24.3992 25 4.1591 3.1853 40.1400 20.3030 23.9876 9.4030 41.1084 27.7149