第二十三章 现代优化算法
现代优化算法是80年代初兴起的启发式算法。这些算法包括禁忌搜索（tabu search），模拟退火（simulated annealing），遗传算法（genetic algorithms），人工神经网络（neural networks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标－求NP-hard组合优化问题的全局最优解。虽然有这些目标，但NP-hard理论限制它们只能以启发式的算法去求解实际问题。
启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（Ant Colony Algorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。
现代优化算法解决组合优化问题，如TSP（Traveling Salesman Problem）问题，QAP（Quadratic Assignment Problem）问题，JSP（Job-shop Scheduling Problem）问题等效果很好。
§1 模拟退火算法
1.1 算法简介
模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。
如果用粒子的能量定义材料的状态，Metropolis算法用一个简单的数学模型描述了退火过程。假设材料在状态i之下的能量为E(i)，那么材料在温度T时从状态i进入状态j就遵循如下规律：
（1）如果E(j)≤E(i)，接受该状态被转换。
（2）如果E(j)>E(i)，则状态转换以如下概率被接受：
e
其中KT是材料温度。
在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态的概率满足波尔兹曼分布：
PT(x=i)=
E(i)−E(j)
KT
e
j∈S
−
E(i)KTE(j)−KT
∑e
−E(i)KT−E(j)KT
其中x表示材料当前状态的随机变量，S表示状态空间集合。
显然
lim
e
j∈S
T→∞
=
∑e
1 |S|
-450-
其中|S|表示集合S中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，
lim
T→0
e
j∈S
−
E(i)−Emin
KT−
E(j)−Emin
KT
=lim
T→0
j∈Smin
e
−
E(i)−Emin
KT
∑e
e
−
∑e
−
E(j)−Emin
KT
+
j∉Smin
∑e
−
E(j)−Emin
KT
=lim
T→0
E(i)−Emin
KT−
E(j)−Emin
KT
j∈Smin
∑e
⎧1
若i∈Smin⎪||S=⎨min ⎪0 其它⎩
其中Emin=minE(j)且Smin={i|E(i)=Emin}。
j∈S
上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。
想应用于优化问题就可以得到模拟退火寻优方法。
考虑这样一个组合优化问题：优化函数为f:x→R，其中x∈S，它表示优化问题的一个可行解，R={y|y∈R,y>0}，S表示函数的定义域。N(x)⊆S表示x的一个邻域集合。
首先给定一个初始温度T0和该优化问题的一个初始解x(0)，并由x(0)生成下一个解x'∈N(x(0))，是否接受x'作为一个新解x(1)依赖于下面概率：
+
+
⎧1 若f(x')
⎪
P(x(0)→x')=⎨−f(x')−f(x(0))K呢丠
T0
⎪其它⎩e
换句话说，如果生成的解x'的函数值比前一个解的函数值更小，则接受x(1)=x'作为
接受x'作为一个新解。 一个新解。否则以概率e
泛泛地说，对于某一个温度Ti和该优化问题的一个解x(k)，可以生成x'。接受x'作为下一个新解x(k+1)的概率为：
−
f(x')−f(x(0))
T0
⎧1 若f(x')
P(x(k)→x')=⎨−f(x')−f(x(k)) （1）
Ti
⎪其它⎩e
在温度Ti下，经过很多次的转移之后，降低温度Ti，得到Ti+1
过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。
我们注意到，在每个Ti下，所得到的一个新状态x(k+1)完全依赖于前一个状态
x(k)，x(0),L,x(k−1)无关，因此这是一个马尔可夫过程。使用马尔可夫过程对上述模拟退火的步骤进行分析，结果表明：从任何一个状态x(k)生成x'的概率，在N(x(k))中是均匀分布的，且新状态x'被接受的概率满足式（1），那么经过有限次的转换，在温度Ti下的平衡态xi的分布由下式给出：
-451-
Pi(Ti)=
e
−
f(xi)T−f(xi)Ti
（2）
∑e
j∈S
当温度T降为0时，xi的分布为：
⎧1
若xi∈Smin⎪*
Pi=⎨|Smin|
⎪0 其它⎩
并且
xi∈Smin
∑P
*
i
=1
这说明如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率1被找到。因此可以说模拟退火算法可以找到全局最优解。
在模拟退火算法中应注意以下问题：
（1）理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。
（2）要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续m次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值Te，或连续几个温度下转换过程没有使状态发生变化算法就结束。
（3）选择初始温度和确定某个可行解的邻域的方法也要恰当。 1.2 应用举例
已知敌方100个目标的经度、纬度如表1所示。
表1 经度和纬度数据表
经度 纬度 经度 纬度 经度 纬度 53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 56.5432 21.4188 10.8198 16.2529 22.7891 23.1045 20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 26.2418 18.1760 44.0356 13.5401 28.9836 25.9879 28.2694 29.0011 32.1910 5.8699 36.4863 29.7284 8.9586 24.6635 16.5618 23.6143 10.5597 15.1178 8.1519 9.5325 22.1075 18.5569 0.1215 18.8726 31.9499 17.6309 0.7732 0.4656 47.4134 23.7783 43.5474 3.9061 53.3524 26.7256 30.8165 13.4595 23.9222 7.6306 51.9612 22.8511 12.7938 15.7307 21.5051 24.0909 15.2548 27.2111 6.2070 5.1442 17.1168 20.0354 34.1688 22.7571 9.4402 3.9200 52.1181 0.4088 9.5559 11.4219 24.4509 6.5634 37.5848 16.8474 35.6619 9.9333 24.4654 3.1644 14.4703 13.6368 19.8660 15.1224 3.1616 4.2428 58.6849 27.1485 39.5168 16.9371 56.5089 13.7090 -452-
经度 纬度
30.3313 6.9348 10.1584 12.4819 31.4847 8.9640 38.4722 20.1731 0.9718 28.1477 50.2111 10.2944 48.2077 16.8889 41.8671 3.5667 27.7133 5.0706 4.9568 8.3669 49.2430 16.7044 11.5812 14.5677 26.7213 28.5667 0.7775 6.9576 18.5245 14.3598 52.5211 15.7957
38.4300 8.4648 13.7909 1.9510 40.8801 14.2978 39.9494 29.5114 8.0831 27.6705 1.3496 16.8359 15.7320 19.5697 6.9909 23.1804 4.1591 3.1853 51.8181 23.0159 34.0574 23.3960 58.8289 14.5229 47.5099 24.0664 9.1556 14.1304 49.9816 6.0828 11.5118 17.3884 38.3392 19.9950 40.1400 20.3030 8.9983 23.6440 23.0624 8.4319 18.6635 6.7436 10.1121 27.2662 53.7989 0.2199 19.3635 17.6622 44.0398 16.2635 24.6543 19.6057 23.9876 9.4030 50.1156 23.7816 19.9857 5.7902 52.8423 27.2880 28.7812 27.6659 33.6490 0.3980 36.9545 23.0265 39.7139 28.4203 36.9980 24.3992 41.1084 27.7149
我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。
这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵D=(dij)102×102，其中dij表示表示i,j两点的距离，i,j=1,2,L,102，这里D为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。
上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设A,B两点的地理坐标分别为(x1,y1)，(x2,y2)，过A,B两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0
度经线圈所在的平面为XOZ平面建立三维直角坐标系。则A,B两点的直角坐标分别为：
A(Rcosx1cosy1,Rsinx1cosy1,Rsiny1) B(Rcosx2cosy2,Rsinx2cosy2,Rsiny2) 其中R=6370为地球半径。
A,B两点的实际距离
⎛
⎜ d=Rarccos⎞⎟， 化简得
d=Rarccos[cos(x1−x2)cosy1cosy2+siny1siny2]。 求解的模拟退火算法描述如下： （1）解空间
解空间S可表为{1,2,L,101,102}的所有固定起点和终点的循环排列集合，即 S={(π1,L,π102)|π1=1,(π2,L,π101)为{2,3,L,101}的循环排列,π102=102} 其中每一个循环排列表示侦察100个目标的一个回路，πi=j表示在第i−1次侦察目标j，初始解可选为(1,2,L,102)，本文中我们先使用Monte Carlo方法求得一个较好的初始解。
（2）目标函数
此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求
minf(π1,π2,L,π102)=
∑dππ
i=1
101
ii+1
-453-
而一次迭代由下列三步构成：
（3）新解的产生 ① 2变换法
任选序号 u,v（u
π1Lπu−1πvπv−1Lπu+1πuπv+1Lπ102
② 3变换法
任选序号u,v和w，将u和v之间的路径插到w之后，对应的新路径为（设u
π1Lπu−1πv+1LπwπuLπvπw+1Lπ102 （4）代价函数差
对于2变换法，路径差可表示为
Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1)
（5）接受准则
Δf
P=⎨
exp(/)0fTf−ΔΔ≥⎩
如果Δf
α
（6）降温
利用选定的降温系数α进行降温即：T←αT，得到新的温度，这里我们取=0.999。
（7）结束条件
−30
用选定的终止温度e=10，判断退火过程是否结束。若T
我们编写如下的matlab程序如下： 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=zeros(102); %距离矩阵d 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;
rand('state',sum(clock)); for j=1:1000 S=[1 1+randperm(100),102]; temp=0;
-454-
for i=1:101
temp=temp+d(S(i),S(i+1)); end if temp
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
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
break; end end
% 输出巡航路径及路径长度 S0,Sum
计算结果为44小时左右。其中的一个巡航路径如图1所示。
4332211 图1 模拟退火算法求得的巡航路径示意图
§2 遗传算法
-455-
2.1 遗传算法简介
遗传算法（Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色体的基因并随机变异某些染色体的基因后生成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。其实现方法如下：
（1） 根据具体问题确定可行解域，确定一种编码方法，能用数值串或字符串表示可行解域的每一解。
（2） 对每一解应有一个度量好坏的依据，它用一函数表示，叫做适应度函数。 （3） M、交叉概率p、变异概率pm、进化终止条件。 为便于计算，一般来说，每一代群体的个体数目都取相等。群体规模越大、越容易找到最优解，但由于受到计算机的运算能力的限制，群体规模越大，计算所需要的时间也相应的增加。进化终止条件指的是当进化到什么时候结束，它可以设定到某一代进化结束，也可以根据找出近似最优是否满足精度要求来确定。表2列出了生物遗传概念在遗传算法中的对应关系。
表2 生物遗传概念在遗传算法中的对应关系
生物遗传概念 适者生存 个体 染色体 基因 适应性 种群 交配 变异
遗传算法中的作用
算法停止时，最优目标值的解有最大的可能被留住 解 解的编码
解中每一分量的特征 适应度函数值
根据适应度函数值选取的一组解 通过交配原则产生一组新解的过程 编码的某一分量发生变化的过程
2.2模型及算法
我们用遗传算法研究1.2中的问题。 求解的遗传算法的参数设定如下： 种群大小：M=50 最大代数：G=1000
交叉率：pc=1，交叉概率为1能保证种群的充分进化。 变异率：pm=0.1， 一般而言，变异发生的可能性较小。
（1） 编码策略
采用十进制编码，用随机数列ω1ω2Kω102作为染色体，其中0
-456-
应，例如一个9目标问题的一个染色体为
[0.23，0.82，0.45，0.74，0.87，0.11，0.56，0.69，0.78]
其中编码位置i代表目标i，位置i的随机数表示目标i在巡回中的顺序，我们将这些随机数按升序排列得到如下巡回：
6－1－3－7－8－4－9－2－5
（2） 初始种群
本文中我们先利用经典的近似算法—改良圈算法求得一个较好的初始种群。即对于初始圈
C=π1Lπu−1πuπu+1Lπv−1πvπv+1Lπ102，2≤u
π1Lπu−1πvπv−1Lπu+1πuπv+1Lπ102
记Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1)，若Δf
目标函数为侦察所有目标的路径长度，适应度函数就取为目标函数。我们要求
minf(π1,π2,L,π102)=
∑dππ
i=1
101
ii+1
（4） 交叉操作
我们的交叉操作采用单点交叉。设计如下，对于选定的两个父代个体
f1=ω1ω2Kω102，f2=ω'1ω'2Kω'102，我们随机地选取第t个基因处为交叉点，则经过交叉运算后得到的子代编码为s1和s2，s1的基因由f1的前t个基因和f2的后102−t个基因构成，s2的基因由f2的前t个基因和f1的后102−t个基因构成，例如：
f1=[0,0.14,0.25,0.27,|0.29,0.54,L,0.19,1] f2=[0,0.23,0.44,0.56,|0.74,0.21,L,0.24,1]
设交叉点为第四个基因处，则
s1=[0,0.14,0.25,0.27,|0.74,0.21,L,0.24,1] s2=[0,0.23,0.44,0.56,|0.29,0.54,L,0.19,1]
交叉操作的方式有很多种选择，我们应该尽可能选取好的交叉方式，保证子代能继承父代的优良特性。同时这里的交叉操作也蕴含了变异操作。 （5） 变异操作
变异也是实现群体多样性的一种手段，同时也是全局寻优的保证。具体设计如下，按照给定的变异率，对选定变异的个体，随机地取三个整数，满足1
（6） 选择
采用确定性的选择策略，也就是说选择目标函数值最小的M个个体进化到下一代，这样可以保证父代的优良特性被保存下来。 2.3 模型求解及结论
编写MATLAB程序如下： tic
clc,clear
load sj.txt %加载敌方100个目标的数据
-457-
x=sj(:,1:2:8); x=x(:); y=sj(:,2:2:8); y=y(:); sj=[x y]; d1=[70,40];
sj0=[d1;sj;d1]; sj=sj0*pi/180; d=zeros(102); %距离矩阵d 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';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))
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; rand('state',sum(clock)); %遗传算法实现过程 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)
by=floor(w*rand(1))+1;
-458-
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=sj0(path,1);yy=sj0(path,2); plot(xx,yy,'-o')
计算结果为40小时左右。其中的一个巡航路径如图2所示。
图2 遗传算法求得的巡航路径示意图
§3 禁忌搜索算法
3.1 禁忌搜索算法简介
禁忌搜索算法是组合优化算法的一种，是局部搜索算法的扩展。禁忌搜索算法是人工智能在组合优化算法中的一个成功应用。禁忌搜索算法的特点是采用了禁忌技术。所谓禁忌就是禁止重复前面的工作。禁忌搜索算法用一个禁忌表记录下已经到达过的局部最优点，在下一次搜索中，利用禁忌表中的信息不再或有选择地搜索这些点。
禁忌搜索算法实现的技术问题是算法的关键。禁忌搜索算法涉及侯选集合、禁忌对象、评价函数、特赦规则、记忆频率信息等概念。
（1）邻域
-459-
在组合优化中，距离的概念通常不再适用，但是在一点附近搜索另一个下降的点仍然是组合优化数值求解的基本思想。因此，需要重新定义邻域的概念。
一个组合优化问题可用三参数(D,F,f)表示，其中D表示决策变量的定义域，F表示可行解区域，F中的任何一个元素称为该问题的可行解，f表示目标函数。满足
f(x*)=min{f(x)|x∈F}的可行解x*称为该问题的最优解。 定义1 对于组合优化问题(D,F,f)，D上的一个映射：
D
N:x∈D→N(x)∈2
D
称为一个邻域映射，其中2表示D的所有子集组成的集合，N(x)称为x的邻域，y∈N(x)称为x的一个邻居。
（2）侯选集合
侯选集合由邻域中的邻居组成。常规的方法是从邻域中选择若干个目标值或评价值最佳的邻居入选。
（3）禁忌对象和禁忌长度
禁忌表中的两个主要指标是禁忌对象和禁忌长度。禁忌算法中，由于我们要避免一些操作的重复进行，就要将一些元素放到禁忌表中以禁止对这些元素进行操作，这些元素就是我们指的禁忌对象。禁忌长度是被禁对象不允许选取的迭代次数。一般是给被禁对象x一个数（禁忌长度）t，要求对象x在t步迭代内被禁，在禁忌表中采用tabu(x)=t记忆，每迭代一步，该项指标做运算tabu(x)=t−1，直到tabu(x)=0时解禁。于是，我们可将所有元素分成两类，被禁元素和自由元素。禁忌长度t的选取可以有多种方法，例如t=常数，或t=[n]，其中n为邻域中邻居的个数；这种规则容易在算法中实现。
（4）评价函数
评价函数是侯选集合元素选取的一个评价公式，侯选集合的元素通过评价函数值来选取。以目标函数作为评价函数是比较容易理解的。目标值是一个非常直观的指标，但有时为了方便或易于计算，会采用其他函数来取代目标函数。
（5）特赦规则
在禁忌搜索算法的迭代过程中，会出现侯选集中的全部对象都被禁忌，或有一对象被禁，但若解禁则其目标值将有非常大的下降情况。在这样的情况下，为了达到全局最优，我们会让一些禁忌对象重新可选。这种方法称为特赦，相应的规则称为特赦规则。
（6）记忆频率信息
在计算的过程中，记忆一些信息对解决问题是有利的。如一个最好的目标值出现的频率很高，这使我们有理由推测：现有参数的算法可能无法再得到更好的解。根据解决问题的需要，我们可以记忆解集合、被禁对象组、目标值集合等的出现频率。
频率信息有助于进一步加强禁忌搜索的效率。我们可以根据频率信息动态控制禁忌的长度。一个最佳的目标值出现的频率很高，有理由终止计算而将此值认为是最优值。
3.2 模型及求解
我们用禁忌搜索算法研究如下的两个问题： （1）研究1.2中同样的问题。
（2）我方有三个基地，经度、纬度分别为（70,40），（72,45），（68,48）。假设我方所有无人侦察机的速度都为1000公里/小时。三个基地各派出一架飞机侦察敌方目标，怎样划分任务，才能使时间最短，且任务比较均衡。
3.2.1 问题（1）的求解
-460-
求解的禁忌搜索算法描述如下： （1）解空间
解空间S可表为{1,2,L,101,102}的所有固定起点和终点的循环排列集合，即
S={(π1,L,π102)|π1=1,(π2,L,π101)为{2,3,L,101}的循环排列,π102=102} 其中每一个循环排列表示侦察100个目标的一个回路，πi=j表示第i−1次侦察目标j。
（2）目标函数
目标函数为侦察所有目标的路径长度。我们要求
minf(π1,π2,L,π102)=
∑dππ
i=1
101
ii+1
（3）候选集合
定义2 对于路径π1Lπu−1πuπu+1Lπv−1πvπv+1Lπ102，交换u与v之间的顺序，得到的新路径为
π1Lπu−1πvπv−1Lπu+1πuπv+1Lπ102，
称为原路径二邻域的一个邻居。
定义3 对于路径π1Lπuπu+1Lπv−1πvπv+1Lπwπw+1Lπ102，将u和v之间的路径插到w之后，得到的新路径为
π1Lπu−1πv+1LπwπuLπvπw+1Lπ102，
称为原路径三邻域的一个邻居。
如果要考虑当前解的全部二邻域（或三邻域）的邻居，将面临着太大的工作量。因此我们用随机选取的方法每次选取50个邻居组成的集合作为侯选集合。而将省下的时间作更多次搜索，这样做同样可以保证较高的精确度，同时可以大大提高算法的效率。 （4）禁忌长度及禁忌对象
对于上述定义的二邻域中的邻居总数为C100，我们的禁忌长度取为
2
t=70≈C100。
2
我们把禁忌表设计成一个循环队列，初始化禁忌表H=Φ。从候选集合C中选出一个向量x，如果x∉H，并且H不满，则把向量x添加到禁忌表中；如果H已满，则最早进入禁忌表的向量出列，向量x进入到队列。
（5）评价函数
可以用目标函数作为评价函数，但是这样每选取一个新的路径都得去计算总时间，计算量比较大。对于上述二邻域中的邻居作为侯选集合，每一个新路径中只有两条边发生了变化，因此将目标函数的差值作为评价函数可以极大地提高算法的效率。评价函数取为
Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1)
禁忌搜索算法的流程如下： STEP1 选定一个初始解x
now
及给以禁忌表H=Φ；
now
STEP2 若满足停止条件，停止计算；否则，在x足禁忌要求的侯选集Can_N(x
now
的邻域N(H,x
now
)中选出满
)，在Can_N(xnow)中选一个评价值最佳的解
-461-
xnext，xnow:=xnext，更新禁忌表H，重复STEP2。
利用Matlab程序求得，我们的巡航时间大约在41小时左右，其中的一个巡航路径如图3所示。
4332211
图3 禁忌搜索算法求得的巡航路径示意图
3.2.2 问题（2）的求解 对于这个问题，我们的基本想法是，先根据敌方基地的分布特点将敌方的基地大体划分在三个区域之内，并使三架侦察机分别对这三个区域的敌方基地进行侦察，求取各自的最短时间。然后对任务不均衡区域之中的点做适当调整。
我们解决问题的步骤如下：
（1）划分子图。要达到比较均衡，应使每架飞机的巡航时间基本相同，由于敌方目标的分布较均匀，可以将敌方目标的地理位置图分成面积基本相同的三部分。如过点
443
作一条斜线，过点（68，48）以斜率k2=作一条斜线，568
把基地所在的地区划分成三部分G1,G2,G3（见图4）。
（70，40）以斜率k1=
544332211
图4 任务划分示意图
（2）再对每个子图Gi(i=1,2,3)分别运用禁忌搜索算法求得其最短侦察路径ψi，和最短侦察时间t(ψi)(i=1,2,3)。
（3）均衡任务。定义均衡率 η=
min{t(ψ1),t(ψ2),t(ψ3)}
×100%
max{t(ψ1),t(ψ2),t(ψ3)}
若η接近于1，则上面划分的任务就可以接受。否则的话，根据t(ψi)（i=1,2,3）的大小用局部搜索算法，调整k1,k2的斜率，从而调整各分区内点的个数，直至任务达到
-462-
均衡。
最后计算结果如下
两直线的斜率分别为k1=0.68，k2=0.6484，均衡率η=97.57%。 三架飞机侦察完所有目标所用的时间分别为
t(ψ1)=20.401小时，t(ψ2)=20.910小时，t(ψ3)=20.729小时。
§4 改进的遗传算法
4.1 引言
无人机航路规划问题实际上是一个组合优化问题，是优化理论中的NP—hard问题。因为其解空间不连续，解邻域表达困难，所以难以用通常的算法求解。遗传算法作为现代优化算法之一，其主要特点是对非线性极值问题能以概率1跳出局部最优解，找到全局最优解。而遗传算法这种跳出局部最优寻找全局最优特性都基于算法中的交叉和变异。在传统遗传算法的结构中，变异操作在交叉操作基础上进行，强调的是交叉作用，认为变异只是一个生物学背景机制。在具体交叉操作中，人们通常采用断点交叉法（段交叉）多点交叉与均匀交叉，其中断点交叉是指随机地在基因序列中选择一个断点，然后交换双亲上断点右端的所有染色体。在变异操作中，变异算子一般是用Guassian 分布的随机变异来实现[46,47]。近年来，也有学者尝试用Cauchy 分布的随机序列来实现变异[48], 希望通过Cauchy 分布宽大的两翼特性实现更大范围的变异，以利于找到全局最优解。[49]从理论上分析了采用Cauchy分布随机变异进化算法的局部收敛性。[50]进一步把二者结合起来, 采用两种分布的线性叠加，但仿真结果显示，算法改进效果并不十分明显。文献[51]将生物进化看成是随机性加上反馈，并指出其中的随机性主要是由系统的内在因素所引起，而不是由外部环境的随机扰动所造成。而混沌系统在其混沌域中表现为随机性，它是确定系统内部随机性的反映，不同于外在的随机特性。本节根据以上特点对基于求解航路规划的遗传算法进行改进，首先将变异操作从交叉操作中分离出来，使其成为独立的并列于交叉的寻优操作，在具体遗传操作中，混沌与遗传操作联系在一起，在交叉操作中，以“门当户对”原则进行个体的配对，利用混沌序列确定交叉点，实行强度最弱的单点交叉，以确保算法收敛精度，削弱和避免寻优抖振问题；在变异操作中，利用混沌序列对染色体中多个基因进行变异，以避免算法早熟。
下面我们研究1.2中同样的问题。 4.2 模型及算法
与标准的遗传算法相比，我们做了如下的两点改进。 （1）交叉操作
我们的交叉操作采用改进型交叉。具体设计如下：首先以“门当户对”原则，对父代个体进行配对，即对父代以适应度函数（目标函数）值进行排序，目标函数值小的与小的配对，目标函数值大的与大的配对。然后利用混沌序列确定交叉点的位置，最后对确定的交叉项进行交叉。例如(Ω1,Ω2)配对，他们的染色体分别是Ω1=ω1ω2Kω102，
22
Ω2=ω12ω2Kω102，采用Logistic混沌序列x(n+1)=4x(n)(1−x(n))产生一个2到101之间
111
的正整数，具体步骤如下：
取一个(0,1)随机初始值，然后利用x(n+1)=4x(n)(1−x(n))迭代一次产生1个
(0,1)上的混沌值，保存以上混沌值作为产生下一代交叉项的混沌迭代初值，再把这个
值分别乘以100并加上2，最后取整即可。假如这个数为33，那么我们对(Ω1,Ω2)染色
'
,Ω'2) 体中相应的基因进行交叉，得到新的染色体(Ω1
-463-
'111112111Ω1=ω1ω2ω3ω4ω5Lω33ω34Lω60ω61L
22221222Ω'2=ω12ω2ω3ω4ω5Lω33ω34Lω60ω61L
很明显这种单点交叉对原来的解改动很小，这可以削弱避免遗传算法在组合优化应用中产生的寻优抖振问题，可以提高算法收敛精度。
（2）变异操作
变异也是实现群体多样性的一种手段，是跳出局部最优，全局寻优的重要保证。在本文具体变异算子设计如下，首先根据给定的变异率(本文选为0.02)，随机地取两个在2到101之间的整数，对这两个数对应位置的基因进行变异，具体变异以当前的基因值为初值利用混沌序列x(n+1)=4x(n)(1−x(n))进行适当次数的迭代，得到变异后新的基因值，从而得到新的染色体。
4.3 仿真结果对比及算法性能分析 我们研究1.2中的同样问题。 计算的MATLAB程序如下： tic
clc,clear
load sj.txt %加载敌方100个目标的数据 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=zeros(102); %距离矩阵d 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';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))
c1(m+1:n)=c1(n:-1:m+1); end end
-464-
end end
J(k,c1)=1:102; end
J=J/102;
J(:,1)=0;J(:,102)=1;
rand('state',sum(clock)); %遗传算法实现过程 A=J;
for k=1:dai %产生0～1间随机数列进行编码 B=A;
%交配产生子代B for i=1:2:w
ch0=rand;ch(1)=4*ch0*(1-ch0); for j=2:50
ch(j)=4*ch(j-1)*(1-ch(j-1)); end
ch=2+floor(100*ch); temp=B(i,ch);
B(i,ch)=B(i+1,ch); B(i+1,ch)=temp; end
%变异产生子代C
by=find(rand(1,w)
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),:);
-465-
end
path=IX(IZ(1),:) long=DZ(1) toc
在仿真试验中，我们对文中航路规划问题分别利用断点交叉和换位变异结合的遗传算法，多点交叉和移位变异结合的遗传算法和文中提出的改进算法进行求解比较。表3是各种算法种群规模（M=50）和迭代次数（G=100）都相同时连续20次求解的平均值（公里），算法平均运算时间（秒）。
指标
平均航路距离 算法执行时间
表3 算法性能比较表 断点交叉算法 多点交叉算法 文中改进算法
本文从算法结构到具体的遗传操作都进行了改进，其中变异操作从交叉操作中分离出来，使得遗传算法也可以通过并行计算实现，提高算法实现效率。其次改进后的算法，分别采用变化强度不同的交叉操作和变异操作，其中交叉操作采用强度最弱的单点交叉，保证了算法收敛精度，削弱和避免算法因交叉强度大而产生的寻优抖振问题。当然单一的单点交叉很容易使算法早熟，文中采用较大强度的多个基因变异正好解决早熟问题。从仿真结果可以看到改进后的算法效果较为明显。
§5 Matlab遗传算法工具
5.1 遗传算法与直接搜索工具箱概述
遗传算法与直接搜索（Genetic Algorithm and Direct Search，简称GADS）工具箱是一系列函数的集合，它们扩展了优化工具箱和Matlab数值计算环境。遗传算法与直接搜索工具箱包含了要使用遗传算法和直接搜索算法来求解优化问题的一些例程。这些算法使我们能够求解那些标准优化工具箱范围之外的各种优化问题。所有工具箱函数都是Matlab的M文件，这些文件由实现特定优化算法的Matlab语句所写成。
使用语句
Type function_name
就可以看到这些函数的Matlab代码。我们也可以通过编写自己的M文件来实现和扩展遗传算法和直接搜索工具箱的性能，也可以将该工具箱与Matlab的其它工具箱或Simulink结合使用来求解优化问题。
工具箱函数可以通过图形界面或Matlab命令行来访问，它们是用Matlab语言编写的，对用户开放，因此可以查看算法，修改源代码或生成用户函数。
遗传算法与直接搜索工具箱可以帮助我们求解那些不易用传统方法解决的问题，譬如旅行商问题等。
遗传算法与直接搜索工具箱有一个精心设计的图形用户界面，可以帮助我们直观、方便、快速地求解最优化问题。 1．功能特点
遗传算法与直接搜索工具箱的功能特点如下：
（1）图形用户界面和命令行函数可用来快速地描述问题、设置算法选项以及监控进程。
（2）具有多个选项的遗传算法工具可用于问题创建、适应度计算、选择、交叉和变异。
（3）直接搜索工具实现了一种模式搜索方法，其选项可用于定义网格尺寸、表决
-466-
方法和搜索方法。
（4）遗传算法与直接搜索工具箱函数可与Matlab的优化工具箱或其它的Matlab程序结合使用。
（5）支持自动的M代码生成。 2．图形用户界面和命令行函数
遗传算法工具函数可以通过命令行和图形用户界面来使用遗传算法。直接搜索工具函数也可以通过命令行和图形用户界面来进行访问。图形用户界面可用来快速地定义问题，设置算法选项，对优化问题进行详细定义。
遗传算法和直接搜索工具箱还同时提供了用于优化管理、性能监控及终止准则定义的工具，同时还提供了大量的标准算法选项。
在优化运行的过程中，可以通过修改选项来细化最优解，更新性能结果。用户也可以提供自己的算法选项来定制工具箱。
3．使用其它函数和求解器
遗传算法与直接搜索工具箱和Matlab及优化工具箱是紧密结合在一起的。用户可以用遗传算法或直接搜索算法来寻找最佳起始点，然后利用优化工具箱或用Matlab程序来进一步寻找最优解。通过结合不同的算法，可以充分地发挥Matlab和工具箱的功能以提高求解的质量。对于某些特定问题，使用这种方法还可以得到全局（最优）解。
4．显示、监控和输出结果
遗传算法与直接搜索工具箱还包括一系列绘图函数，用来可视化优化结果。这些可视化功能直观地显示了优化的过程，并且允许在执行过程中进行修改。
使用输出函数可以将结果写入文件，产生用户自己的终止准则，也可以写入用户自己的图形界面来运行工具箱求解器。除此之外，还可以将问题的算法选项导出，以便日后再将它们导入到图形界面中去。
5.2 使用遗传算法工具初步 5.2.1 遗传算法使用规则
遗传算法是一种基于自然选择、生物进化过程来求解问题的方法。在每一步中，遗传算法随机地从当前种群中选择若干个体作为父辈，并且使用它们产生下一代的子种群。在连续若干代之后，种群朝着优化解的方向进化。我们可以用遗传算法来求解各种不适宜于用标准算法求解的优化问题，包括目标函数不连续、不可微、随机或高度非线性的问题。
遗传算法在每一步使用下列三类规则从当前种群来创建下一代： （1）选择规则（Selection Rules）：选择对下一代种群有贡献的个体（称为父辈）。 （2）交叉规则（Crossover Rules）：将两个父辈结合起来构成下一代的子辈种群。 （3）变异规则（Mutation Rules）：施加随机变化给父辈个体来构成子辈。
遗传算法与标准优化算法主要在两个方面有所不同，它们的比较情况归纳于表4中。
表4 遗传算法与标准优化算法比较 标准算法 遗传算法
每次迭代产生一个单点，点的序列逼近一个优每次迭代产生一个种群，种群逼近一个优化解 化解
通过确定性的计算在该序列中选择下一个点 通过随机进化选择计算来选择下一代种群
5.2.2 遗传算法使用方式
遗传算法工具有两种使用方式：
（1）以命令行方式调用遗传算法函数ga。
-467-
（2）通过图形用户界面使用遗传算法工具。
1．在命令行使用遗传算法，可以用下列语法调用遗传算法函数ga： [x, fval]= ga(@fitnessfun,nvars,A,b,Aeq,beq,LB,UB,@nonlcon,options)
其中@fitnessfun是适应度函数句柄，nvars是适应度函数的独立变量的个数，options是一个包含遗传算法选项参数的数据结构，其它参数的含义与非线性规划fmincon中的参数相同。
函数所给出的结果为：x最终值到达的点，这里x为行向量，fval适应度函数的最终值。
例1 求下列问题的解
maxf(x)=2x1+3x1+3x2+x2+x3
2
⎧x1+2x12+x2+2x2+x3≤10⎪22
⎪x1+x1+x2+x2−x3≤50⎪2x+x2+2x+x≤40⎪123
s.t. ⎨1
2
⎪x1+x3=2⎪x+2x≥1
2⎪1
⎪⎩x1≥0,x2,x3不约束
22
解：（1）编写适应度函数（文件名为fun1.m） function y=fun1(x); %x为行向量 c1=[2 3 1]'; c2=[3 1 0]'; y=x*c1+x.^2*c2; y=-y;
（2）编写非线性约束函数（文件名为fun2.m） function [f,g]=fun2(x);
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10 x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50 2*x(1)+x(1)^2+2*x(2)+x(3)-40]; g=x(1)^2+x(3)-2; （3）主函数
a=[-1 -2 0;-1 0 0];b=[-1;0];
[x,y]=ga(@fun1,3,a,b,[],[],[],[],@fun2); x,y=-y
遗传算法程序的运行结果每一次都是不一样的，我们要运行多次，找一个最好的结果。
2．通过GUI使用遗传算法
遗传算法工具有一个图形用户界面GUI，它使我们可以使用遗传算法而不用工作在命令行方式。打开遗传算法工具，可键入以下命令：
gatool
使用遗传算法工具首先必须输入下列信息：
（1）Fitness function（适应度函数）－欲求最小值的目标函数。输入适应度函数的形式为@fitnessfun，其中fitnessfun.m是计算适应度函数的M文件。符号@产生一个对于函数fitnessfun的函数句柄。
（2）Number of variables（变量个数）－适应度函数输入向量的长度。
其它参数的含义，我们就不一一介绍了，与优化工具的图形用户界面解法中的参数一样。如果某个参数的值为空，可以输入空矩阵[]，或者什么都不输入。
-468-
例2 用图形用户界面求解例1。
在Matlab命令窗口运行gatool，打开图形用户界面，如图5所示，填入有关的参数，未填入的参数取值为空或者为默认值，然后用鼠标点一下“start”按钮，就得到求解结果，再使用“file”菜单下的“Export to Workspace…”选项，把计算结果输出到Matlab工作空间中去。
图5 遗传算法图形用户界面解法
注意：Matlab工作空间中必须存在变量a，b。或者编一个小程序，对a,b进行赋值，然后运行一下该小程序。
5.3 直接搜索工具 直接搜索的命令为
[x,fval] = patternsearch(@fun,x0,A,b,Aeq,beq,LB,UB, @nonlcon,options) 打开直接搜索的图形用户界面的命令为 psearchtool
例3 用直接搜索算法求解例1。
解 编写程序如下（所有程序放在一个文件中）： function example23_3 a=[-1 -2 0;-1 0 0];b=[-1;0];
[x,y]=patternsearch(@fun1,rand(1,3),a,b,[],[],[],[],@fun2); %初始值必须为行向量 x,y=-y
%定义目标函数 function y=fun1(x); %x为行向量 c1=[2 3 1]'; c2=[3 1 0]'; y=x*c1+x.^2*c2; y=-y; %定义非线性约束函数 function [f,g]=fun2(x);
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10 x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50 2*x(1)+x(1)^2+2*x(2)+x(3)-40]; g=x(1)^2+x(3)-2;
直接搜索算法每次的计算结果也是不一样的。
-469-
§6 蚁群算法
6.1 蚁群算法简介
蚁群是自然界中常见的一种生物，人们对蚂蚁的关注大都是因为“蚁群搬家，天要下雨”之类的民谚。然而随着近代仿生学的发展，这种似乎微不足道的小东西越来越多地受到学者们地关注。1991年意大利学者M. Dorigo等人首先提出了蚁群算法，人们开始了对蚁群的研究：相对弱小，功能并不强大的个体是如何完成复杂的工作的（如寻找到食物的最佳路径并返回等）。在此基础上一种很好的优化算法逐步发展起来。
蚁群算法的特点是模拟自然界中蚂蚁的群体行为。科学家发现，蚁群总是能够发现从蚁巢到食物源的最短路径。经研究发现，蚂蚁在行走过的路上留下一种挥发性的激素，蚂蚁就是通过这种激素进行信息交流。蚂蚁趋向于走激素积累较多的路径。找到最短路径的蚂蚁总是最早返回巢穴，从而在路上留下了较多的激素。由于最短路径上积累了较多的激素，选择这条路径的蚂蚁就会越来越多，到最后所有的蚂蚁都会趋向于选择这条最短路径。基于蚂蚁这种行为而提出的蚁群算法具有群体合作，正反馈选择，并行计算等三大特点，并且可以根据需要为人工蚁加入前瞻、回溯等自然蚁所没有的特点。
在使用蚁群算法求解现实问题时，先生成具有一定数量蚂蚁的蚁群，让每一只蚂蚁建立一个解或解的一部分，每只人工蚁从问题的初始状态出发，根据“激素”浓度来选择下一个要转移到的状态，直到建立起一个解，每只蚂蚁根据所找到的解的好坏程度在所经过的状态上释放与解的质量成正比例的“激素”。之后，每只蚂蚁又开始新的求解过程，直到寻找到满意解。为避免停滞现象，引入了激素更新机制。
6.2 解决TSP问题的蚁群算法描述
现以TSP问题的求解为例说明蚁群系统模型。首先引进如下记号：n为城市的个数；m为蚁群中蚂蚁的数量；dij为两城市i和j之间距离；bi(t)为t时刻位于城市i的蚂蚁的个数，m=∑b(t)；τi
i=1nij(t)为t时刻边弧(i,j)的轨迹强度（即ij连线上残留的
，i,j=1,2,L,n，i≠j；ηij(t)为t时刻边弧(i,j)信息量），且设τij(0)=c（c为常数）
的能见度，反映由城市i转移到城市j的期望程度。
根据上述原理，蚂蚁k(k=1,2,L,m)在运动过程中根据各条路径上的信息量决定转移方向。与真实蚁群系统不同，人工蚁群系统具有一定的记忆功能。随着时间的推移，以前留下的信息逐渐消逝，经n个时刻，蚂蚁完成一次循环，各路径上信息量要作调整。由此得到下述的人工蚁群系统模型：
1）设人工蚁群在并行地搜索TSP的解，并通过一种信息素做媒介相互通信，在每个结点上且和该结点相连的边上以信息素量做搜索下一结点的试探依据，直到找到一个TSP问题的可行解。
2）在时刻t人工蚁k由位置i转移至位置j的转移概率为
αβ⎧τij(t)ηij(t)⎪⎪τα(t)ηβ(t), j∈S kpij(t)=⎨iv （3） iv
v∈S⎪⎪⎩0, j∉S
；β为能见度的相对重要性（β≥0）；S为其中参数α为轨迹的相对重要性（α≥0）
可行点集，即蚂蚁k下一步允许选择的城市。α,β分别反映了蚂蚁在运动过程中所积累的信息及启发式因子在蚂蚁选择路径中所起的不同作用。
-470-
3）当m个人工蚁按（3）式找到了可行解，则将各边的信息量用下式修改。即调整信息量的轨迹强度更新方程为
τij(t+1)=ρτij(t)+Δτij，ρ∈(0,1) （4）
kΔτij=∑Δτij
k=1m
其中Δτij为第k只蚂蚁在本次循环中留在路径(i,j)上的信息量；Δτij为本次循环中路径(i,j)上的信息量的增量；参数ρ为轨迹的持久性；1−ρ为轨迹衰减度，表示信息消逝程度。
对上述系统模型，采用人工蚁群方法求解的算法步骤可归结为：
step 1：NC←0（NC为迭代步数或搜索次数）；各τij和Δτij的初始化；将m个蚂蚁置于n个顶点上。
step 2：将各蚂蚁的初始出发点置于当前解集中；对每个蚂蚁k（k=1,2,L,m）按概率pij转移至下一顶点j；将顶点j置于当前解集。
step 3：计算各蚂蚁的目标函数值zk(k=1,L,m)，记录当前的最好解。
step 4：按更新方程修改轨迹强度。
step 5：NC←NC+1，若NC
若为了简化计算，增加处理较大规模的TSP问题的能力，则可将（4）式修改为： τij(t+1)=ρτij(t)+(1−ρ)Δτij，ρ∈(0,1)
其中 kk
⎧1⎪d,(i,j)∈BEk Δτij=⎨ij
⎪0, 其它⎩
此处BE为本次最优路线上的边集。
6.3 人工蚁群算法性能的讨论
人工蚁群算法是一种基于种群的进化算法。作为一个新兴的研究领域，虽它还远未像GA、SA等算法那样形成系统的分析方法和坚实的数学基础，但目前已有一些基本结果。
在M. Dorigo 三种不同的模型中，循环路径(i,j)上信息量的增量Δτij不同。
1）Ant-quantity system 模型中，
⎧Q⎪d,若第k只蚂蚁在时刻t和t+1之间经过ijk Δτij=⎨ij
⎪0, 其它⎩
2）在Ant-density system 模型中，
⎧Q,若第k只蚂蚁在时刻t和t+1之间经过ijk Δτij=⎨⎩0, 其它
3）在Ant-cycle system 模型中，
-471-
⎧Q⎪,若第k只蚂蚁在本次循环中经过ijkΔτij=⎨Lk
⎪0, 其它⎩
其中Q是反映蚂蚁所留轨迹数量的常数，Lk表示第k只蚂蚁在本次循环中所走路径的长度；且t=0时，τij(0)=c，Δτij=0。算法中模型1）、2）利用的是局部信息，模型3）利用的是整体信息。
人工蚁群算法中，α,β,Q等参数对算法性能也有很大的影响。α值的大小表明留在每个结点上的信息量受重视的程度，α值越大，蚂蚁选择以前选过的点的可能性越大，但过大会使搜索过早陷于局部极小点；β的大小表明启发式信息受重视的程度；Q值会影响算法的收敛速度，Q过大会使算法收敛于局部极小值，过小又会影响算法的收敛速度，随问题规模的增大Q的值也需要随之变化；蚂蚁的数目越多，算法的全局搜索能力越强，但数目加大将使算法的收敛速度减慢。
习题二十三
1．用遗传算法求解下列非线性规划问题：
kkf(x)=(x1−2)2+(x2−1)2 s.t. x1−2x2+1≥0 min x122−x2+1≥0 4
2．学生面试问题
高校自主招生是高考改革中的一项新生事物，现在仍处于探索阶段。某高校拟在全面衡量考生的高中学习成绩及综合表现后再采用专家面试的方式决定录取与否。该校在今年自主招生中，经过初选合格进入面试的考生有N人，拟聘请老师M人。每位学生要分别接受4位老师（简称该学生的“面试组”）的单独面试。面试时，各位老师独立地对考生提问并根据其回答问题的情况给出评分。由于这是一项主观性很强的评价工作，老师的专业可能不同，他们的提问内容、提问方式以及评分习惯也会有较大差异，因此面试同一位考生的“面试组”的具体组成不同会对录取结果产生一定影响。为了保证面试工作的公平性，组织者提出如下要求：
（1）每位老师面试的学生数量应尽量均衡；
（2）面试不同考生的“面试组”成员不能完全相同；
（3）两个考生的“面试组”中有两位或三位老师相同的情形尽量的少；
（4）被任意两位老师面试的两个学生集合中出现相同学生的人数尽量的少。 请回答如下问题：
问题一：设考生数N已知，在满足条件（2）下，说明聘请老师数M至少分别应为多大，才能做到任两位学生的“面试组”都没有两位以及三位面试老师相同的情形。
问题二：请根据（1）～（4）的要求建立学生与面试老师之间合理的分配模型，并就N=379，M=24的情形给出具体的分配方案（每位老师面试哪些学生）及该方案满足（1）～（4）这些要求的情况。
问题三：假设面试老师中理科与文科的老师各占一半，并且要求每位学生接受两位文科与两位理科老师的面试，请在此假设下分别回答问题一与问题二。
问题四：请讨论考生与面试老师之间分配的均匀性和面试公平性的关系。为了保证-472-
面试的公平性，除了组织者提出的要求外，还有哪些重要因素需要考虑，试给出新的分配方案或建议。
3．编写用禁忌搜索算法求解TSP问题的Matlab程序。
4．编写用蚁群算法求解TSP问题的Matlab程序。
5．用遗传算法求解下列非线性整数规划
Max z=x2224x22
1+x2+3x3+4+2x5−8x1−2x2−3x3−x4−2x5
⎧0≤xi≤99(i=1,L,5)⎪⎪⎪x1+x2+x3+x4+x5≤400
⎨x1+2x2+2x3+x
⎪4+6x5≤800
⎪2x1+x2+6x3≤200
⎪⎩x3+x4+5x5≤200
-473-