学习笔记-Matlab算法篇-现代优化算法

现代优化算法

01遗传算法

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

为便于计算，一般来说，每一代群体的个体数目都取相等。群体规模越大、越容易找到最优解，但由于受到计算机的运算能力的限制，群体规模越大，计算所需要的时间也相应的增加。进化终止条件指的是当进化到什么时候结束，它可以设定到某一代进化结束，也可能根据找出近似最优是否满足精度要求来确定。下表列出了生物遗传概念在遗传算法中的对应关系： 

例子(e01)：已知敌方 100 个目标的经度、纬度如下表所示。

经度	纬度	经度_1	纬度_2	经度_4	纬度_5	经度_7	纬度_8
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.105	15.4562	1.9451	0.2057	26.4951	22.1221	31.4847	8.964
26.2418	18.176	44.0356	13.5401	28.9836	25.9879	38.4722	20.1731
28.2694	29.0011	32.191	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.207	5.1442	49.243	16.7044
17.1168	20.0354	34.1688	22.7571	9.4402	3.92	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.866	15.1224	3.1616	4.2428	18.5245	14.3598
58.6849	27.1485	39.5168	16.9371	56.5089	13.709	52.5211	15.7957
38.43	8.4648	51.8181	23.0159	8.9983	23.644	50.1156	23.7816
13.7909	1.951	34.0574	23.396	23.0624	8.4319	19.9857	5.7902
40.8801	14.2978	58.8289	14.5229	18.6635	6.7436	52.8423	27.288
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.649	0.398
1.3496	16.8359	49.9816	6.0828	19.3635	17.6622	36.9545	23.0265
15.732	19.5697	11.5118	17.3884	44.0398	16.2635	39.7139	28.4203
6.9909	23.1804	38.3392	19.995	24.6543	19.6057	36.998	24.3992
4.1591	3.1853	40.14	20.303	23.9876	9.403	41.1084	27.7149

 

我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

 

 5.变异操作：变异也是实现群体多样性的一种手段，同时也是全局寻优的保证。具体设计如下，按照给定的变异率，对选定变异的个体，随机地取三个整数，满足1<𝑢<𝑣<𝑤<102，1<u<v<w<102，把uv之间的基因段插入到w后面。

6.选择：采用确定性的选择策略，也就是说选择目标函数值最小的 M 个个体进化到下一代，这样可以保证父代的优良特性被保存下来。

data1.txt
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

%% 实现遗传算法
tic
clc,clear 
load data1.txt       %加载敌方 100 个目标的数据 
x=data1(:,1:2:8);x=x(:);   % 经度
y=data1(:,2:2:8);y=y(:);   % 维度
data1=[x y]; 
d1=[70,40]; 
data0=[d1;data1;d1]; 
%距离矩阵 d data0
data1=data0*pi/180; 
d=zeros(102); 
for i=1:101 
    for j=i+1:102 
        temp=cos(data1(i,1)-data1(j,1))*cos(data1(i,2))*cos(data1(j,2))+sin(data1(i,2))*sin(data1(j,2)); 
        d(i,j)=6370*acos(temp); 
    end 
end 
d=d+d';L=102;w=50;dai=100; 
%通过改良圈算法选取优良父代 A 
for k=1:w 
    c=randperm(100); 
    c1=[1,c+1,102]; 
    flag=1; 
 while flag>0 
      flag=0; 
   for m=1:L-3 
      for n=m+2:L-1 
        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)) 
           flag=1; 
           c1(m+1:n)=c1(n:-1:m+1); 
        end 
      end 
   end 
 end 
  J(k,c1)=1:102; 
end 
J=J/102;  
J(:,1)=0;J(:,102)=1; 
rng('default'); 
%遗传算法实现过程 
A=J; 
for k=1:dai  %产生 0～1间随机数列进行编码 
    B=A; 
    c=randperm(w); 
%交配产生子代 B 
    for i=1:2:w   
        F=2+floor(100*rand(1)); 
        temp=B(c(i),F:102); 
        B(c(i),F:102)=B(c(i+1),F:102); 
        B(c(i+1),F:102)=temp; 
    end 
%变异产生子代 C 
by=find(rand(1,w)<0.1); 
if length(by)==0 
    by=floor(w*rand(1))+1; 
end 
C=A(by,:); 
L3=length(by); 
for j=1:L3 
   bw=2+floor(100*rand(1,3)); 
   bw=sort(bw); 
   C(j,:)=C(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]); 
end 
   G=[A;B;C]; 
   TL=size(G,1); 
   %在父代和子代中选择优良品种作为新的父代 
   [dd,IX]=sort(G,2);temp(1:TL)=0; 
   for j=1:TL 
       for i=1:101 
           temp(j)=temp(j)+d(IX(j,i),IX(j,i+1)); 
       end 
   end 
     [DZ,IZ]=sort(temp); 
     A=G(IZ(1:w),:); 
end 
path=IX(IZ(1),:) 
long=DZ(1) 
toc 
xx=data0(path,1);yy=data0(path,2); 
plot(xx,yy,'-o') 

结果：

path =

  1 至 19 列

     1    17     3    45    67     2    92    82    48    72    14    27    10    84    18    40    79    77    31

  20 至 38 列

    97    85    65    64    11    76    69    94    70    19    63    62    26    29    34    66    90    86     8

  39 至 57 列

    39    78    47    23    58    81    25    68     7    22    71    37    32    13    24    49    28    57    88

  58 至 76 列

    61    16    91    41     4    73    33    75     5    60     9    54    53    12    55    96    89     6    56

  77 至 95 列

    21    99   101    52   100    44    38    98    51    80    50    15    42    20    30    74    83    87    59

  96 至 102 列

    46    93    43    36    35    95   102


long =

   4.0327e+04

时间已过 0.246457 秒。

 

02模拟退火算法

定义：模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态i之下的能量为E(i)，那么材料在温度T时从状态i进入状态j就遵循如下规律：

1.如果E(j)<=E(i)，接受该状态被转变

2.如果E(j)>E(i)，则状态转换以概率被接受：

式中，K是玻尔兹曼常数，T是材料温度

                                                                                             

注意事项：

（1）理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折中。

（2）要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m 次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值T𝑒T_e值，或连续几个温度下转换过程没有使状态发生变化算法就结束。

（3）选择初始温度和确定某个可行解的邻域的方法也要恰当。

例子（e02）:使用模拟退火算法求解e01的问题

 

 

 

clc,clear 
load sj.txt    %加载敌方100 个目标的数据， 数据按照表格中的位置保存在纯文本文件 sj.txt 中 
x=sj(:,1:2:8);x=x(:); 
y=sj(:,2:2:8);y=y(:); 
sj=[x y]; 
d1=[70,40]; 
sj=[d1;sj;d1]; 
sj=sj*pi/180; 
%距离矩阵 d 
d=zeros(102); 
for i=1:101 
    for j=i+1:102 
        temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)); 
        d(i,j)=6370*acos(temp); 
    end 
end 
d=d+d'; 
S0=[];Sum=inf; 
rng('default'); 
for j=1:1000 
    S=[1 1+randperm(100),102]; 
    temp=0;
    for i=1:101 
        temp=temp+d(S(i),S(i+1)); 
    end 
    if temp<Sum 
        S0=S;Sum=temp; 
    end 
end 
e=0.1^30;L=20000;at=0.999;T=1; 
%退火过程 
for k=1:L 
   %产生新解 
c=2+floor(100*rand(1,2)); 
c=sort(c); 
c1=c(1);c2=c(2); 
  %计算代价函数值 
  df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1)); 
   %接受准则 
  if df<0 
  S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];       
  Sum=Sum+df; 
  elseif exp(-df/T)>rand(1) 
  S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)]; 
  Sum=Sum+df; 
  end 
  T=T*at; 
   if T<e 
       break; 
   end 
end 
 %  输出巡航路径及路径长度 
xx=sj(S0,1);yy=sj(S0,2); 
plot(xx,yy,'-o') 
Sum

 结果：

Sum =

   4.4359e+04

03神经网络算法神

介绍：人工神经网络（artificial neural network，ANN）的基本单元的神经元模型如下：

                                     

有三个基本要素：

1.一组连接（对应于生物神经元的突触），连接强度由各连接上的权值表示，权值为正表示激活，为负表示抑制。

2.一个求和单元，用于求取各输入信号的加权和（线性组合）。

3.一个非线性激活函数，起非线性映射作用并将神经元输出幅度限制在一定范围内(一般限制在(0,1)或者(-1,1)之间)。

其计算公式如下：

                                             

也可以把阈值放到和输入一起：

Matlab中的激活函数如下表示：

函数名	功 能
purelin	线性传递函数
hardlim	硬限幅传递函数
hardlims	对称硬限幅传递函数
satlin	饱和线性传递函数
satlins	对称饱和线性传递函数
logsig	对数 S 形传递函数
tansig	正切 S 形传递函数
radbas	径向基传递函数
compet	竞争层传递函数

网络结构：

1.前馈型网络：各神经元接受前一层的输入，并输出给下一层，没有反馈。结点分为两类，即输入单元和计算单元，每一计算单元可有任意个输入，但只有一个输出（它可耦合到任意多个其它结点作为其输入）。通常前馈网络可分为不同的层，第i层的输入只与第i-1层输出相连，输入和输出结点与外界相连，而其它中间层则称为隐层。

2.反馈型网络：所有结点都是计算单元，同时也可接受输入，并向外界输出。NN 的工作过程主要分为两个阶段：第一个阶段是学习期，此时各计算单元状态不变，各连线上的权值可通过学习来修改；第二阶段是工作期，此时各连接权固定，计算单元状态变化，以达到某种稳定状态。

 从作用效果看，前馈网络主要是函数映射，可用于模式识别和函数逼近。反馈网络按对能量函数的极小点的利用来分类有两种：第一类是能量函数的所有极小点都起作用，这一类主要用作各种联想存储器；第二类只利用全局极小点，它主要用于求解最优化问题。

例子(e03)：蠓虫分类问题

生物学家试图对两种蠓虫（Af与Apf）进行鉴别，依据的资料是触角和翅膀的长度，已经测得了9只Af和6只Apf的数据如下：

翅长
Af	1.24	1.36	1.38	1.38	1.38	1.4	1.48	1.54	1.56
Apf	1.14	1.18	1.2	1.26	1.28	1.3
触角长
Af	1.27	1.74	1.64	1.82	1.9	1.7	1.82	1.9	1.7
Apf	1.82	1.9	1.7	1.82	1.82	2.08

现在的问题是：

1.根据如上资料，如何制定一种方法，正确地区分两类蠓虫。

2.对触角和翼长分别为(1.24,1.80)，(1.28,1.84)与(1.40,2.04)的 3 个标本，用所得到的方法加以识别

3.设 Af 是宝贵的传粉益虫，Apf 是某疾病的载体，是否应该修改分类方法。

如上的问题是有代表性的，它的特点是要求依据已知资料制定一种分类方法，类别是已经给定的（Af 或 Apf）。

求解(多层前馈神经网络)：为 解 决 上 述 问 题 ， 考 虑 一 个 其 结 构 如 下 图 所 示 的 人 工 神 经 网 络：

                                         

该网络的输出为：

                             

后向传播算法：最速下降法

                              

%% 使用神经网络进行蠓虫分类
clear
p1=[1.24,1.27;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90;
1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08];
p2=[1.14,1.82;1.18,1.96;1.20,1.86;1.26,2.00
1.28,2.00;1.30,1.96];
p=[p1;p2]';
pr=minmax(p);
goal=[ones(1,9),zeros(1,6);zeros(1,9),ones(1,6)];
plot(p1(:,1),p1(:,2),'h',p2(:,1),p2(:,2),'o')
%% 网络生成
net=newff(pr,[3,2],{'logsig','logsig'});
net.trainParam.show = 10;
net.trainParam.lr = 0.05;
net.trainParam.goal = 1e-10;
net.trainParam.epochs = 50000;
net = train(net,p,goal);
%% 网络使用
x=[1.24 1.80;1.28 1.84;1.40 2.04]';
y0=sim(net,p)
y=sim(net,x)
hold on
plot(x(:,1),x(:,2),'*')

结果：

y0 =

  1 至 11 列

    0.9999    0.9999    0.9999    0.9999    0.9998    0.9999    0.9999    0.9999    0.9999    0.0000    0.0000
    0.0001    0.0001    0.0001    0.0001    0.0002    0.0001    0.0001    0.0001    0.0001    1.0000    1.0000

  12 至 15 列

    0.0000    0.0000    0.0000    0.0001
    1.0000    1.0000    1.0000    0.9999


y =

    0.0010    0.0519    0.5222
    0.9990    0.9464    0.4695

 

04禁忌搜索算法

介绍：禁忌搜索算法是组合优化算法的一种，是局部搜索算法的扩展。禁忌搜索算法是人工智能在组合优化算法中的一个成功应用。禁忌搜索算法的特点是采用了禁忌技术。所谓禁忌就是禁止重复前面的工作。禁忌搜索算法用一个禁忌表记录下已经到达过的局部最优点，在下一次搜索中，利用禁忌表中的信息不再或有选择地搜索这些点。禁忌搜索算法实现的技术问题是算法的关键。禁忌搜索算法涉及侯选集合、禁忌对象、评价函数、特赦规则、记忆频率信息等概念。

 2.候选集合：侯选集合由邻域中的邻居组成。常规的方法是从邻域中选择若干个目标值或评价值最佳的邻居入选。

4.评价函数：是侯选集合元素选取的一个评价公式，侯选集合的元素通过评价函数值来选取。以目标函数作为评价函数是比较容易理解的。目标值是一个非常直观的指标，但有时为了方便或易于计算，会采用其他函数来取代目标函数。

5.特赦规则：在禁忌搜索算法的迭代过程中，会出现侯选集中的全部对象都被禁忌，或有一对象被禁，但若解禁则其目标值将有非常大的下降情况。在这样的情况下，为了达到全局最优，我们会让一些禁忌对象重新可选。这种方法称为特赦，相应的规则称为特赦规则。

6.记忆频率信息：在计算的过程中，记忆一些信息对解决问题是有利的。如一个最好的目标值出现的频率很高，这使我们有理由推测：现有参数的算法可能无法再得到更好的解。根据解决问题的需要，我们可以记忆解集合、被禁对象组、目标值集合等的出现频率。频率信息有助于进一步加强禁忌搜索的效率。我们可以根据频率信息动态控制禁忌的长度。一个最佳的目标值出现的频率很高，有理由终止计算而将此值认为是最优值。

例子(p01)：使用禁忌搜索算法求解e01中的问题。

1.解空间：解空间S可表为 {1,2,3,…,102}的所有固定起点和终点的循环排列集合。

2.目标函数：目标函数为侦察所有目标的路径长度。

5.评价函数：可以用目标函数作为评价函数，但是这样每选取一个新的路径都得去计算总时间，计算量比较大。对于上述二邻域中的邻居作为侯选集合，每一个新路径中只有两条边发生了变化，因此将目标函数的差值作为评价函数可以极大地提高算法的效率。评价函数取为：

05蚁群算法

 介绍：蚁群是自然界中常见的一种生物，人们对蚂蚁的关注大都是因为“蚁群搬家，天要下雨”之类的民谚。蚁群算法的特点是模拟自然界中蚂蚁的群体行为。科学家发现，蚁群总是能够发现从蚁巢到食物源的最短路径。经研究发现，蚂蚁在行走过的路上留下一种挥发性的激素，蚂蚁就是通过这种激素进行信息交流。蚂蚁趋向于走激素积累较多的路径。找到最短路径的蚂蚁总是最早返回巢穴，从而在路上留下了较多的激素。由于最短路径上积累了较多的激素，选择这条路径的蚂蚁就会越来越多，到最后所有的蚂蚁都会趋向于选择这条最短路径。基于蚂蚁这种行为而提出的蚁群算法具有群体合作，正反馈选择，并行计算等三大特点，并且可以根据需要为人工蚁加入前瞻、回溯等自然蚁所没有的特点。

介绍：在使用蚁群算法求解现实问题时，先生成具有一定数量蚂蚁的蚁群，让每一只蚂蚁建立一个解或解的一部分，每只人工蚁从问题的初始状态出发，根据“激素”浓度来选择下一个要转移到的状态，直到建立起一个解，每只蚂蚁根据所找到的解的好坏程度在所经过的状态上释放与解的质量成正比例的“激素”。之后，每只蚂蚁又开始新的求解过程，直到寻找到满意解。为避免停滞现象，引入了激素更新机制。

例子(p02)：使用蚁群算法解决TSP问题

1.设人工蚁群在并行地搜索 TSP 的解，并通过一种信息素做媒介相互通信，在每个结点上且和该结点相连的边上以信息素量做搜索下一结点的试探依据，直到找到一个TSP 问题的可行解。

2.在时刻 t 人工蚁 k 由位置 i 转移至位置 j 的转移概率为。其中参数α为轨迹的相对重要性，β为能见度的相对重要性；S为可行点集，即蚂蚁k下一步允许选择的城市。Α、β分别反映了蚂蚁在运动过程中所积累的信息及启发式因子在蚂蚁选择路径中所起的不同作用。