遗传算法实例及MATLAB程序解析

       遗传算法Genetic Algorithms，GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终得到最优解或准最优解。它必须做以下操作∶初始群体的产生、求每一个体的适应度、根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色体的基因并随机变异某些染色体的基因生 成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。其实现方法如下∶

实例解析

已知100个目标的经纬度信息如下（第一列为经度，第二例为纬度，以此类推）：

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

       我方有一个基地，经度和纬度为（70，40）。假设我方飞机的速度为1000km/h。我方派一架飞机从基地出发，侦察完所有目标，再返回原来的基地。在每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。
        这是一个旅行商问题。给我方基地编号为1，目标依次编号为2，，…，101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵$D=(d_{ij}{)}_{102\times 102}$，其中$d_{ij}$表示i，j两点的距离，i、j=1，2，…，102，这里$D$为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。
       上面问题中给定的是地理坐标（经度和纬度），必须求两点间的实际距离。设A，B两点的地理坐标分别为$（x_1，y_1），（x_2，y_2）$，过A，B两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点$O$，以赤道平面为$XOY$平面，以0度经线圈所在的平面为$XOZ$平面建立三维直角坐标系。则A，B两点的直角坐标分别为：

$$A(Rcos{x}_{1}cos{y}_1,Rsin{x}_{1}cos{y}_1,Rsin{y}_1),$$$$B(Rcos{x}_{2}cos{y}_2,Rsin{x}_{2}cos{y}_2,Rsin{y}_2),$$

式中∶R=6370为地球半径。

A，B两点的实际距离:
$$d=Rarccos(\frac{OA\cdot OB}{|OA|\cdot |OB|}),$$

用MATLAB求解程序如下：

%遗传算法
clc,clear
sj0=load('sj.txt');       %加载100个目标的数据
x=sj0(:,1:2:8); x=x(:);   %取经度
y=sj0(:,2:2:8); y=y(:);   %取纬度
sj=[x y]; d1=[70,40];     %基地经纬度（70,40）
sj=[d1;sj;d1]; sj=sj*pi/180;  %单位化成弧度
d=zeros(102); %距离矩阵d的初始值 100个目标+两次经过起点

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';                    %d为一个实对称矩阵，表示两点间的双向距离
w=50; g=100;               %w为种群的个数，g为进化的代数
rand('state',sum(clock));  %初始化随机数发生器

for k=1:w                  %通过改良圈算法选取初始种群
    c=randperm(100);       %产生1，...，100的一个随机排列  
    c1=[1,c+1,102];        %生成初始解
    for t=1: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; %修改圈
        end
      end
    end
   if flag==0
      J(k,c1)=1:102; break %记录下较好的解并退出当前层循环
   end
   end
end

J(:,1)=0; J=J/102;   %把整数序列转换成[0,1]区间上的实数，即转换成染色体编码
for k=1:g            %该层循环进行遗传算法的操作 
    A=J;             %交配产生子代B的初始染色体
    c=randperm(w);   %产生下面交叉操作的染色体对 
    for i=1:2:w  
        F=2+floor(100*rand(1));   %产生交叉操作的地址
        temp=A(c(i),[F:102]);     %中间变量的保存值
        A(c(i),[F:102])=A(c(i+1),[F:102]); %交叉操作
        A(c(i+1),F:102)=temp;  
    end
    by=[];                  %为了防止下面产生空地址，这里先初始化
while ~length(by)
    by=find(rand(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]); %交换位置
end
   G=[J;A;B];           %父代和子代种群合在一起
   [SG,ind1]=sort(G,2); %把染色体翻译成1，...,102的序列ind1
   num=size(G,1); long=zeros(1,num); %路径长度的初始值
   for j=1:num
       for i=1:101
           long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度
       end
   end
     [slong,ind2]=sort(long); %对路径长度按照从小到大排序
     J=G(ind2(1:w),:);       %精选前w个较短的路径对应的染色体
end
path=ind1(ind2(1),:), flong=slong(1)  %解的路径及路径长度
xx=sj(path,1);yy=sj(path,2);
plot(xx,yy,'-V')                       %画出路径

运行结果如下：

path =

  1 至 28 列

     1    17     3    45    67     2    92    87    83    82    48    72    14    27    10    84    18    40    20    30    74    42    15     9     5    60    79    77

  29 至 56 列

    31    97    85    65    64    11    76    69    94    70    19    63    62    26    29    34    66    90    86     8    39    78    47    23    58    81    25    68

  57 至 84 列

     7    22    71    37    32    13    24    61    49    28    57    88    16    91    41     4    73    33    75    54    53    12    89     6    96    55    44    38

  85 至 102 列

    50    80    51    98   100    56    21    99   101    52    46    59    93    43    36    35    95   102


flong =

   4.0096e+04

注：这是我用Markdown写的第一篇文章，之中难免有所纰漏，还望指出.