第二十三章现代优化算法 现代优化算法是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模拟退火算法 算法简介 模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不 同结构对应于粒子的不同能量水平。在高温条件下,粒子的能量较高,可以自由运动和 重新排列。在低温条件下,粒子能量较低。如果从高温开始,非常缓慢地降温(这个过 称为退火),粒子就可以在每个温度下达到热平衡。当系统完全被冷却时,最终 成处于低能状态的晶体 如果用粒子的能量定义材料的状态, Metropolis算法用一个简单的数学模型描述了 退火过程。假设材料在状态i之下的能量为E(),那么材料在温度T时从状态i进入状 态j就遵循如下规律: (1)如果E()≤E(),接受该状态被转换。 (2)如果E()>E(i),则状态转换以如下概率被接受 E()=E 其中K是物理学中的波尔兹曼常数,T是材料温度。 在某一个特定温度下,进行了充分的转换之后,材料将达到热平衡。这时材料处 于状态i的概率满足波尔兹曼分布: P 其中x表示材料当前状态的随机变量,S表示状态空间集合。 显然 E(i) lim
-271- 第二十三章 现代优化算法 现代优化算法是 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) ，则状态转换以如下概率被接受： KT E i E j e ( )− ( ) 其中 K 是物理学中的波尔兹曼常数，T 是材料温度。 在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处 于状态i 的概率满足波尔兹曼分布： ∑∈ − − = = j S KT E j KT E i T e e P x i ( ) ( ) ( ) 其中 x 表示材料当前状态的随机变量， S 表示状态空间集合。 显然 | | 1 lim ( ) ( ) S e e j S KT E j KT E i T = ∑∈ − − →∞
其中S|表示集合S中状态的数量。这表明所有状态在高温下具有相同的概率。而当温 度下降时 lim =lim T→0 m组 = lime A7 若 I Smin I i∈Smin -E(l-Emin ∑e 其它 其中Emn=minE()且Smn={iE(1)=Emn}。 上式表明当温度降至很低时,材料会以很大概率进入最小能量状态, 假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思 想应用于优化问题就可以得到模拟退火寻优方法。 考虑这样一个组合优化问题:优化函数为f:x→R,其中x∈S,它表示优化 问题的一个可行解,R={yy∈R,y>0},S表示函数的定义域。N(x)三S表示x 的一个邻域集合 首先给定一个初始温度T。和该优化问题的一个初始解x(0),并由x(0)生成下一个 解x∈N(x(O),是否接受x作为一个新解x(1)依赖于下面概率: 若f(x)<f(x(0) P(x(0)→x)={x-xo) 其它 换句话说,如果生成的解x的函数值比前一个解的函数值更小,则接受x(1)=x作为 个新解。否则以概率e 接受x作为一个新解。 泛泛地说,对于某一个温度T和该优化问题的一个解x(k),可以生成x。接受x 作为下一个新解x(k+1)的概率为: 若f(x)<f(x(k) P(x(k)→x)={_/x)=xk (1) 其它 在温度T下,经过很多次的转移之后,降低温度T,得到T1<T。在T下重复上述 过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问 题寻优的结果。 我们注意到,在每个T下,所得到的一个新状态x(k+1)完全依赖于前一个状态 x(k),可以和前面的状态x(0),…,x(k-1)无关,因此这是一个马尔可夫过程。使用马 尔可夫过程对上述模拟退火的步骤进行分析,结果表明:从任何一个状态x(k)生成x的 概率,在N(x(k)中是均匀分布的,且新状态x被接受的概率满足式(1),那么经过 有限次的转换,在温度T下的平衡态x的分布由下式给出 -272-
-272- 其中| S | 表示集合 S 中状态的数量。这表明所有状态在高温下具有相同的概率。而当温 度下降时， ∑ ∑ ∑∉ − − ∈ − − − − → ∈ − − − − → + = min min min min min min min ( ) ( ) ( ) 0 ( ) ( ) 0 lim lim j S KT E j E j S KT E j E KT E i E T j S KT E j E KT E i E T e e e e e ⎪ ⎩ ⎪ ⎨ ⎧ ∈ = = ∑∈ − − − − → 其它 若 0 | | 1 lim min min ( ) ( ) 0 min min min i S S e e j S KT E j E KT E i E T 其中 min ( ) min E E j j∈S = 且 { | ( ) } min Emin S = i E i = 。 上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。 假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思 想应用于优化问题就可以得到模拟退火寻优方法。 考虑这样一个组合优化问题：优化函数为 → + f : x R ，其中 x ∈ S ，它表示优化 问题的一个可行解， = { | ∈ , > 0} + R y y R y ，S 表示函数的定义域。N(x) ⊆ S 表示 x 的一个邻域集合。 首先给定一个初始温度T0 和该优化问题的一个初始解 x(0) ，并由 x(0) 生成下一个 解 x'∈ N(x(0)) ，是否接受 x' 作为一个新解 x(1) 依赖于下面概率： ⎪ ⎩ ⎪ ⎨ ⎧ < → = − − 其它 若 ( ') ( (0)) 0 1 ( ') ( (0)) ( (0) ') T f x f x e f x f x P x x 换句话说，如果生成的解 x' 的函数值比前一个解的函数值更小，则接受 x(1) = x' 作为 一个新解。否则以概率 0 ( ') ( (0)) T f x f x e − − 接受 x' 作为一个新解。 泛泛地说，对于某一个温度Ti 和该优化问题的一个解 x(k) ，可以生成 x' 。接受 x' 作为下一个新解 x(k +1) 的概率为： ⎪ ⎩ ⎪ ⎨ ⎧ < → = − − 其它 若 ( ') ( ( )) 0 1 ( ') ( ( )) ( ( ) ') T f x f x k e f x f x k P x k x （1） 在温度Ti 下，经过很多次的转移之后，降低温度Ti ，得到Ti+1 < Ti 。在Ti+1 下重复上述 过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问 题寻优的结果。 我们注意到，在每个Ti 下，所得到的一个新状态 x(k +1) 完全依赖于前一个状态 x(k) ，可以和前面的状态 x(0),L, x(k −1) 无关，因此这是一个马尔可夫过程。使用马 尔可夫过程对上述模拟退火的步骤进行分析，结果表明：从任何一个状态 x(k) 生成 x' 的 概率，在 N(x(k)) 中是均匀分布的，且新状态 x' 被接受的概率满足式（1），那么经过 有限次的转换，在温度Ti 下的平衡态 i x 的分布由下式给出：
PC/、c (2) f(r) 当温度T降为0时,x,的分布为 P={|S 若x,∈Sm 0 其它 并且 P"=1 这说明如果温度下降十分缓慢,而在每个温度都有足够多次的状态转移,使之在每一个 度下达到热平衡,则全局最优解将以概率1被找到。因此可以说模拟退火算法可以找 到全局最优解。 在模拟退火算法中应注意以下问题 (1)理论上,降温过程要足够缓慢,要使得在每一温度下达到热平衡。但在计算 机实现中,如果降温速度过缓,所得到的解的性能会较为令人满意,但是算法会太慢, 相对于简单的搜索算法不具有明显优势。如果降温速度过快,很可能最终得不到全局最 优解。因此使用时要综合考虑解的性能和算法速度,在两者之间采取一种折衷 (2)要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续m次的 转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定 为一个较小的值冖’或连续几个温度下转换过程没有使状态发生变化算法就结束 (3)选择初始温度和确定某个可行解的邻域的方法也要恰当 12应用举例 例已知敌方100个目标的经度、纬度如表1所示。 经度和纬度数据表 经度纬度 经度 经度纬度 纬度 30.33136.9348 5654322141881081981625292278912310 101584124819 □20105015456219 264951 2131484 60 440356135401 289836 38.472 282694290011 32.1910586993648632972840971828477 89586246635165618361431059715.11850211102944 8.151995325 0.1215188726 48.207716.8889 31949917.6309 474134237783 41.86713.56 53.352426.7256 308165134595277133507 23922276306519612228511127938157307 4956883669 21505124090915254827211 6207051442 492430167044 171682003543416882275719440239200115812145677 5211810408 9.555914219 24450965634 584816.8474 35661999333 2446543.16440.777569576 14.470313.6368 19866015.1224 586849 3951681693 7957 〖38430084648 5l8181230159 89983236440 50.115623.7816 137909195103405742339602306248431919.985757902 273
-273- ∑∈ − − = j S T f x T f x i i i i i e e P T ( ) ( ) ( ) （2） 当温度T 降为 0 时， i x 的分布为： ⎪ ⎩ ⎪ ⎨ ⎧ ∈ = 其它 若 0 | | 1 min min * x S P S i i 并且 1 min * ∑ = x ∈S i i P 这说明如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个 温度下达到热平衡，则全局最优解将以概率 1 被找到。因此可以说模拟退火算法可以找 到全局最优解。 在模拟退火算法中应注意以下问题： （1）理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算 机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢， 相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最 优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。 （2）要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m 次的 转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定 为一个较小的值Te ，或连续几个温度下转换过程没有使状态发生变化算法就结束。 （3）选择初始温度和确定某个可行解的邻域的方法也要恰当。 1.2 应用举例 例 已知敌方 100 个目标的经度、纬度如表 1 所示。 表 1 经度和纬度数据表 经度 纬度 经度 纬度 经度 纬度 经度 纬度 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
408801142978 14.5229 18663567436 528423272880 39949429.5114 24.0664 10.1121272662 28.781227.6659 349616835949981660828193631762369545230265 15732019569715118173884440398162635397139284203 383392199950246543196057369980243992 4.15913.1853 40.140020.3030 23.987694030 41.108427.7149 我方有一个基地,经度和纬度为(70,40)。假设我方飞机的速度为1000公里/小时。 我方派一架飞机从基地出发,侦察完敌方所有目标,再返回原来的基地。在敌方每一目 标点的侦察时间不计,求该架飞机所花费的时间(假设我方飞机巡航时间可以充分长) 这是一个旅行商问题。我们依次给基地编号为1,敌方目标依次编号为2,3, 101,最后我方基地再重复编号为102(这样便于程序中计算)距离矩阵D=(d1)x02, 其中d表示表示,j两点的距离,j=12,…102,这里D为实对称矩阵。则问题是 求一个从点1出发,走遍所有中间点,到达点102的一个最短路径。 上面问题中给定的是地理坐标(经度和纬度),我们必须求两点间的实际距离。设 A,B两点的地理坐标分别为(x1,y1),(x2,y2),过A,B两点的大圆的劣弧长即为两点 的实际距离。以地心为坐标原点O,以赤道平面为XOy平面,以0度经线圈所在的平 面为YOZ平面建立三维直角坐标系。则A,B两点的直角坐标分别为 A(Rcos x cos y,, Rsin x, cos y,, Rsin yu) B(RcOS x2 cos y2, Rsin x2 cos y2, Rsin y2) 其中R=6370为地球半径 A,B两点的实际距离 )A·OB d= Marco 化简得 d=Rarccos[cos(x, -x2)cos y, cos y2+sin y, sin y2] 求解的模拟退火算法描述如下 (1)解空间 解空间S可表为{1,2,,101102}的所有固定起点和终点的循环排列集合,即 S={(x1…,m1a)|x1=1,(x2…,ro)为{23…101)的循环排列,r1a2=102} 其中每一个循环排列表示侦察100个目标的一个回路,丌=j表示在第i次侦察j点 初始解可选为(1,2,…102),本文中我们使用 Monte carlo方法求得一个较好的初始解。 (2)目标函数 此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求 minf(x,z2…,zm)=∑dx 而一次迭代由下列三步构成 (3)新解的产生 ①2变换法
-274- 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）。假设我方飞机的速度为 1000 公里/小时。 我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目 标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。 这是一个旅行商问题。我们依次给基地编号为 1，敌方目标依次编号为 2，3，…， 101，最后我方基地再重复编号为 102（这样便于程序中计算）。距离矩阵 102 102 ( ) D = dij × ， 其中 dij 表示表示i, j 两点的距离，i, j = 1,2,L,102 ，这里 D 为实对称矩阵。则问题是 求一个从点 1 出发，走遍所有中间点，到达点 102 的一个最短路径。 上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设 A, B 两点的地理坐标分别为( , ) 1 1 x y ，( , ) 2 2 x y ，过 A, B 两点的大圆的劣弧长即为两点 的实际距离。以地心为坐标原点O ，以赤道平面为 XOY 平面，以 0 度经线圈所在的平 面为 XOZ 平面建立三维直角坐标系。则 A, B 两点的直角坐标分别为： ( cos cos , sin cos , sin ) 1 1 1 1 1 A R x y R x y R y ( cos cos , sin cos , sin ) 2 2 2 2 2 B R x y R x y R y 其中 R = 6370为地球半径。 A, B 两点的实际距离 ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ ⋅ = OB OB d R OA OA arccos ， 化简得 arccos[cos( ) cos cos sin sin ] 1 2 1 2 1 2 d = R x − x y y + y y 。 求解的模拟退火算法描述如下： （1）解空间 解空间 S 可表为{1,2,L,101,102 }的所有固定起点和终点的循环排列集合，即 {( , , ) | 1,( , , ) {2,3, ,101} , 102} S = π 1 L π 102 π 1 = π 2 L π 101 为 L 的循环排列 π 102 = 其中每一个循环排列表示侦察 100 个目标的一个回路， j π i = 表示在第i 次侦察 j 点， 初始解可选为(1,2,L,102) ，本文中我们使用 Monte Carlo 方法求得一个较好的初始解。 （2）目标函数 此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求 ∑= + = 101 1 1 2 102 1 min ( , , , ) i i i f π π L π dπ π 而一次迭代由下列三步构成： （3）新解的产生 ① 2 变换法
任选序号u,v(u<v)交换l与v之间的顺序,此时的新路径为 ②3变换法 任选序号u,和w,将l和v之间的路径插到w之后,对应的新路径为(设 u<y<w) 丌1…丌u-1+…丌丌u…n,n+1 (4)代价函数差 对于2变换法,路径差可表示为 Af =(d r, +d m)-(d, (5)接受准则 △f<0 p(-A/7)4f≥0 如果Δ∫<0,则接受新的路径。否则,以概率eXp(-Δ/T)接受新的路径,即若 p(-A/7)大于0到1之间的随机数则接受。 (6)降温 利用选定的降温系数α进行降温即:T←-α,得到新的温度,这里我们取 (7)结束条件 用选定的终止温度e=10-30,判断退火过程是否结束。若T<e,算法结束,输出 当前状态。 我们编写如下的 matlab程序如下 load sj.txt%加载敌方100个目标的数据,数据按照表格中的位置保存在纯文本 文件stxt中 sj(:1:28)x=x(:) s(:228)y=y() dl=70,40] %距离矩阵d Zeros(102) for i=1: 101 for j=i+1: 102 temp=cos(sj(i, 1)-sjG, 1))*cos(sj(i, 2)*cos(sj(, 2))+sin(sj (1, 2)*sin(sj(, 2)); nd So=U; Sum=inf and(state, sum(clock) forj=1:1000 S=[l 1+randperm100), 102
-275- 任选序号 u,v （u < v ）交换u 与v 之间的顺序，此时的新路径为： π 1Lπ u−1 π vπ v+1Lπ u+1 π uπ v+1Lπ 102 ② 3 变换法 任选序号 u,v 和 w ，将 u 和 v 之间的路径插到 w 之后，对应的新路径为（设 u < v < w ） π 1Lπ u−1 π v+1Lπ wπ u Lπ vπ w+1Lπ 102 （4）代价函数差 对于 2 变换法，路径差可表示为 ( ) ( ) −1 +1 −1 +1 Δ = + − + u v u v u u v v f dπ π dπ π dπ π dπ π （5）接受准则 ⎩ ⎨ ⎧ −Δ Δ ≥ Δ < = exp( / ) 0 1 0 f T f f P 如果Δf < 0 ，则接受新的路径。否则，以概率exp(−Δf /T) 接受新的路径，即若 exp(−Δf /T)大于 0 到 1 之间的随机数则接受。 （6）降温 利用选定的降温系数α 进行降温即： T ←αT ，得到新的温度，这里我们取 α = 0.999 。 （7）结束条件 用选定的终止温度 30 10− e = ，判断退火过程是否结束。若T < e ，算法结束，输出 当前状态。 我们编写如下的 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 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; rand('state',sum(clock)); for j=1:1000 S=[1 1+randperm(100),102]; temp=0;
temp=temp+d(s(i), S(i+D)); SO=S: SI e=0.1^30,L=20000at=0.999;T=1 %退火过程 for k=1- L %产生新解 C=2+floor(100 rand(1, 2)) %计算代价函数值 dfd(So(cl-1), SO(c2))+d(so(cl),So(c2+1)d(So(cl-1), SO(cl))-d(so(c2), SO(c2+1)) %接受准则 if df<o S0=[S0(lcl-1)S0(c2-1:c1)SO(c2+1:102) Sum=Sum+df. elseif exp(-df/rprand(1) 0=S0(lcl-1)0(2-1cl)S0(c2+1:102) fT<e %输出巡航路径及路径长度 SO Sum 计算结果为44小时左右。其中的一个巡航路径如图1所示。 图1模拟退火算法求得的巡航路径示意图
-276- for i=1:101 temp=temp+d(S(i),S(i+1)); end if temprand(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 % 输出巡航路径及路径长度 S0,Sum 计算结果为 44 小时左右。其中的一个巡航路径如图 1 所示。 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 0 5 1 0 1 5 2 0 2 5 3 0 3 5 4 0 图 1 模拟退火算法求得的巡航路径示意图
§2遗传算法 2.1遗传算法简介 遗传算法( Genetic Algorithms,简称GA)是一种基于自然选择原理和自然遗传机 制的搜索(寻优)算法,它是模拟自然界中的生命进化机制,在人工系统中实现特定目 标的优化。遗传算法的实质是通过群体搜索技术,根据适者生存的原则逐代进化,最终 得到最优解或准最优解。它必须做以下操作:初始群体的产生、求每一个体的适应度、 根据适者生存的原则选择优良个体、被选出的优良个体两两配对,通过随机交叉其染色 体的基因并随机变异某些染色体的基因后生成下一代群体,按此方法使群体逐代进化 直到满足进化终止条件。其实现方法如下 (1)根据具体问题确定可行解域,确定一种编码方法,能用数值串或字符串表示 可行解域的每一解。 (2)对每一解应有一个度量好坏的依据,它用一函数表示,叫做适应度函数,适 应度函数应为非负函数。 (3)确定进化参数群体规模M、交叉概率P。、变异概率Pm、进化终止条件。 为便于计算,一般来说,每一代群体的个体数目都取相等。群体规模越大、越容 易找到最优解,但由于受到计算机的运算能力的限制,群体规模越大,计算所需要的时 间也相应的增加。进化终止条件指的是当进化到什么时候结束,它可以设定到某一代进 化结束,也可能根据找出近似最优是否满足精度要求来确定。表1列出了生物遗传概念 在遗传算法中的对应关系。 2生物遗传概念在遗传算法中的对应关系 生物遗传概念 遗传算法中的作用 适者生存 算法停止时,最优目标值的解有最大的可能被留住 个体 解的编码 基因 解中每一分量的特征 适应性 适应度函数值 种群 根据适应度函数值选取的一组解 交配 通过交配原则产生一组新解的过程 变异 编码的某一分量发生变化的过程 22模型及算法 我们用遗传算法研究3.1.2中的问题 求解的遗传算法的参数设定如下 种群大小:M=50 最大代数:G=1000 交叉率:P=1,交叉概率为1能保证种群的充分进化 变异率:Pm=0.1,一般而言,变异发生的可能性较小 (1)编码策略
-277- §2 遗传算法 2.1 遗传算法简介 遗传算法（Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机 制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目 标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终 得到最优解或准最优解。它必须做以下操作：初始群体的产生、求每一个体的适应度、 根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色 体的基因并随机变异某些染色体的基因后生成下一代群体，按此方法使群体逐代进化， 直到满足进化终止条件。其实现方法如下： （1） 根据具体问题确定可行解域，确定一种编码方法，能用数值串或字符串表示 可行解域的每一解。 （2） 对每一解应有一个度量好坏的依据，它用一函数表示，叫做适应度函数，适 应度函数应为非负函数。 （3） 确定进化参数群体规模 M 、交叉概率 pc 、变异概率 pm 、进化终止条件。 为便于计算，一般来说，每一代群体的个体数目都取相等。群体规模越大、越容 易找到最优解，但由于受到计算机的运算能力的限制，群体规模越大，计算所需要的时 间也相应的增加。进化终止条件指的是当进化到什么时候结束，它可以设定到某一代进 化结束，也可能根据找出近似最优是否满足精度要求来确定。表 1 列出了生物遗传概念 在遗传算法中的对应关系。 表 2 生物遗传概念在遗传算法中的对应关系 生物遗传概念 遗传算法中的作用 适者生存 算法停止时，最优目标值的解有最大的可能被留住 个体 解 染色体 解的编码 基因 解中每一分量的特征 适应性 适应度函数值 种群 根据适应度函数值选取的一组解 交配 通过交配原则产生一组新解的过程 变异 编码的某一分量发生变化的过程 2.2 模型及算法 我们用遗传算法研究 3.1.2 中的问题。 求解的遗传算法的参数设定如下： 种群大小： M = 50 最大代数：G = 1000 交叉率： pc = 1，交叉概率为 1 能保证种群的充分进化。 变异率： pm = 0.1， 一般而言，变异发生的可能性较小。 （1） 编码策略
采用十进制编码,用随机数列O1O2…O1a2作为染色体,其中0<O1<1 (i=2,3…101),O1=0,O12=1;每一个随机序列都和种群中的一个个体相对 应,例如一个9城市问题的一个染色体为 [0.23,0.82,0.45,0.74,0.87,0.11,0.56,0.69,0.78] 其中编码位置i代表城市i,位置i的随机数表示城市i在巡回中的顺序,我们将这些随 机数按升序排列得到如下巡回 (2)初始种群 本文中我们先利用经典的近似算法一改良圈算法求得一个较好的初始种群。即对于 初始圈C=x…xnn…-,x…m, 2≤t<v≤101 2≤丌,<丌,≤101,交换l与ν之间的顺序,此时的新路径为 n1∵丌u-1n,y-1…丌+17ay+1…丌102 记=(dx-,+dx.x,)-(d-x.+dx,x),若4<0,则以新的路经修改旧的路 经,直到不能修改为止 (3)目标函数 目标函数为侦察所有目标的路径长度,适应度函数就取为目标函数。我们要求 minf(x1,x2,…,zn)=∑x (4)交叉操作 我们的交叉操作采用单点交叉。设计如下,对于选定的两个父代个体 f1=12…Ol,f2=ω2…,lm2,我们随机地选取第t个基因处为交叉点,则 经过交叉运算后得到的子代编码为S和S2,S1的基因由f的前t个基因和f2的后 02-1个基因构成,S2的基因由/2的前t个基因和f的后102-1个基因构成,例如 f1=[0,0.14,0.25,027,0.29,0.54,…0.19,1 f/2=[0,0.23,0.44,0.56,0.74,021,…0.24,1 设交叉点为第四个基因处,则 s1=[0,0.14,0.25,0.27,0.74,0.21,…0.24,1 s2=[0,0.23,0440.56,029,0.54,…0.19,1 交叉操作的方式有很多种选择,我们应该尽可能选取好的交叉方式,保证子代能继 承父代的优良特性。同时这里的交叉操作也蕴含了变异操作。 (5)变异操作 变异也是实现群体多样性的一种手段,同时也是全局寻优的保证。具体设计如下, 按照给定的变异率,对选定变异的个体,随机地取三个整数,满足1<u<v<w<102, 把u,v之间(包括和)的基因段插到w后面 (6)选择 采用确定性的选择策略,也就是说选择目标函数值最小的M个个体进化到下一代,这 样可以保证父代的优良特性被保存下来 2、3模型求解及结论 编写 MATLAB程序如下 -278-
-278- 采用十进制编码，用随机数列 ω1 ω2Kω102 作为染色体，其中 0 < ωi < 1 （i = 2,3,L,101）， 0 ω1 = ，ω102 = 1；每一个随机序列都和种群中的一个个体相对 应，例如一个 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 < v ≤ 101 ， 2 ≤ π u < π v ≤ 101，交换u 与v 之间的顺序，此时的新路径为： π 1Lπ u−1π vπ v−1Lπ u+1 π uπ v+1Lπ 102 记 ( ) ( ) −1 +1 −1 +1 Δ = + − + u v u v u u v v f dπ π dπ π dπ π dπ π ，若 Δf < 0 ，则以新的路经修改旧的路 经，直到不能修改为止。 （3） 目标函数 目标函数为侦察所有目标的路径长度，适应度函数就取为目标函数。我们要求 ∑= + = 101 1 1 2 102 1 min ( , , , ) i i i f π π L π dπ π （4） 交叉操作 我们的交叉操作采用单点交叉。设计如下，对于选定的两个父代个体 1 = ω1ω 2 Kω102 f ， 102 ' 2 ' 1 ' f 2 = ω ω Kω ，我们随机地选取第 t 个基因处为交叉点，则 经过交叉运算后得到的子代编码为 1 s 和 2 s ， 1 s 的基因由 1f 的前 t 个基因和 2 f 的后 102 − t 个基因构成， 2 s 的基因由 2 f 的前 t 个基因和 1f 的后102 − t 个基因构成，例如： [0, 0.14, 0.25, 0.27,| 0.29, 0.54, ,0.19, 1] f1 = L [0, 0.23, 0.44, 0.56,| 0.74, 0.21, ,0.24, 1] f 2 = L 设交叉点为第四个基因处，则 [0, 0.14, 0.25, 0.27,| 0.74, 0.21, ,0.24, 1] s1 = L [0, 0.23, 0.44, 0.56,| 0.29, 0.54, ,0.19, 1] s2 = L 交叉操作的方式有很多种选择，我们应该尽可能选取好的交叉方式，保证子代能继 承父代的优良特性。同时这里的交叉操作也蕴含了变异操作。 （5） 变异操作 变异也是实现群体多样性的一种手段，同时也是全局寻优的保证。具体设计如下， 按照给定的变异率，对选定变异的个体，随机地取三个整数，满足1 < u < v < w < 102 ， 把u,v 之间（包括u 和v ）的基因段插到 w 后面。 （6） 选择 采用确定性的选择策略，也就是说选择目标函数值最小的 M 个个体进化到下一代，这 样可以保证父代的优良特性被保存下来。 2.3 模型求解及结论 编写 MATLAB 程序如下： tic
load sj. txt %加载敌方100个目标的数据 x=S(,1:2:8)x=x() y=s(:22:8)y=y(:) dl=70,40} j0=ldl; sj; d1]; %距离矩阵d s=sj0°pi/180, d=zeros(102) for i= 1: 101 temp=cos(sj(i, I)-sj(, 1))*cos(sj(1, 2))*cos(sj(, 2))+sin(sj (1, 2))*sin(sj(, 2) d(i,j6370*acos(temp d=d+d;L=102;w=50,dai=100 6通过改良圈算法选取优良父代A for k=l:w (100); l=1,c+1,102] while flag>0 flag=0 for n=m+2. L-1 ifd(cl(m,cl(n)+d(cl(m+1)cl(n+1)<d(cl(m)cl(m+1)+d(cl(n)cl(n+1) (n-1:m+1); J(k,cl=1:102; J=J/102 m(clock); %遗传算法实现过程 A=; fork=l:dai%产生0~1间随机数列进行编码 B=A. %交配产生子代B =12: F=2+floor(100*rand(1)); temp=B(c(i), F: 102) B(c(,F:102)=B(c(+1),F:102) B(c(+1),F:102)=ter
-279- 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]; sj0=[d1;sj;d1]; %距离矩阵 d sj=sj0*pi/180; 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';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; 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 -floor(w*rand(1))+ C=A(by L3=length(by) for =1: L3 bw=2+floor(100*rand(1, 3)) bw-sort(bw); C(-)=C(lbw1}1,bw(2)+1:bw(3).bw(1):bw(2),bw(3)+1:102]) G=[A B; CI TL=Size(gI %在父代和子代中选择优良品种作为新的父代 , IXFsort(G2); temp(I: TLFO for =l: TL for i=1: 101 temp(Ftemp()+d(( 1), IX(, i+l)); d [DZ, IZFsort(temp) end G(IZ(1: w). ) path=IX(IZ(1) xx=sjO(path, 1); yy=sjO(path, 2); plot(xx,yy, -o') 计算结果为40小时左右。其中的一个巡航路径如图2所示。 图2遗传算法求得的巡航路径示意图 3禁忌搜索算法 3.1禁忌搜索算法简介 禁忌搜索算法是组合优化算法的一种,是局部搜索算法的扩展。禁忌搜索算法是人 工智能在组合优化算法中的一个成功应用。禁忌搜索算法的特点是采用了禁忌技术。所 谓禁忌就是禁止重复前面的工作。禁忌搜索算法用一个禁忌表记录下已经到达过的局部
-280- %变异产生子代 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=sj0(path,1);yy=sj0(path,2); plot(xx,yy,'-o') 计算结果为 40 小时左右。其中的一个巡航路径如图 2 所示。 0 10 20 30 40 50 60 70 0 5 10 15 20 25 30 35 40 图 2 遗传算法求得的巡航路径示意图 3 禁忌搜索算法 3.1 禁忌搜索算法简介 禁忌搜索算法是组合优化算法的一种，是局部搜索算法的扩展。禁忌搜索算法是人 工智能在组合优化算法中的一个成功应用。禁忌搜索算法的特点是采用了禁忌技术。所 谓禁忌就是禁止重复前面的工作。禁忌搜索算法用一个禁忌表记录下已经到达过的局部