司守奎《数学建模算法与应用》-遗传算法-保有怀疑

　　遗传算法属于是智能算法里比较好懂的一类算法，在司守奎老师的《数学建模算法与应用》中，遗传算法被应用于求解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') %画出路径
View Code

 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') %画出路径
View Code

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') %画出路径
View Code

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
View Code

 

posted @ 2022-10-15 17:38  不撞楠乔  阅读(350)  评论(0)    收藏  举报