重庆工商大学数学建模算法讲义第23章 现代优化算法

-271-第二十三章现代优化算法现代优化算法是80年代初兴起的启发式算法。这些算法包括禁忌搜索（tabusearch），模拟退火（simulatedannealing），遗传算法（geneticalgorithms），人工神经网络（neuralnetworks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标－求NP-hard组合优化问题的全局昀优解。虽然有这些目标，但NP-hard理论限制它们只能以启发式的算法去求解实际问题。启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（AntColonyAlgorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。现代优化算法解决组合优化问题，如TSP（TravelingSalesmanProblem）问题，QAP（QuadraticAssignmentProblem）问题，JSP（Job-shopSchedulingProblem）问题等效果很好。§1模拟退火算法1.1算法简介模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，昀终形成处于低能状态的晶体。如果用粒子的能量定义材料的状态，Metropolis算法用一个简单的数学模型描述了退火过程。假设材料在状态i之下的能量为)(iE，那么材料在温度T时从状态i进入状态j就遵循如下规律：（1）如果)()(iEjE≤，接受该状态被转换。（2）如果)()(iEjE，则状态转换以如下概率被接受：KTjEiEe)()(−其中K是物理学中的波尔兹曼常数，T是材料温度。在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态i的概率满足波尔兹曼分布：∑∈−−==SjKTjEKTiETeeixP)()()(其中x表示材料当前状态的随机变量，S表示状态空间集合。显然||1lim)()(SeeSjKTjEKTiET=∑∈−−∞→-272-其中||S表示集合S中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，∑∑∑∉−−∈−−−−→∈−−−−→+=minminminminminminmin)()()(0)()(0limlimSjKTEjESjKTEjEKTEiETSjKTEjEKTEiETeeeee⎪⎩⎪⎨⎧∈==∑∈−−−−→其它若0||1limminmin)()(0minminminSiSeeSjKTEjEKTEiET其中)(minminjEESj∈=且})(|{minminEiEiS==。上式表明当温度降至很低时，材料会以很大概率进入昀小能量状态。假定我们要解决的问题是一个寻找昀小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。考虑这样一个组合优化问题：优化函数为+→Rxf:，其中Sx∈，它表示优化问题的一个可行解，}0,|{∈=+yRyyR，S表示函数的定义域。SxN⊆)(表示x的一个邻域集合。首先给定一个初始温度0T和该优化问题的一个初始解)0(x，并由)0(x生成下一个解))0(('xNx∈，是否接受'x作为一个新解)1(x依赖于下面概率：⎪⎩⎪⎨⎧=→−−其它若))0(()'(0))0(()'(1)')0((TxfxfexfxfxxP换句话说，如果生成的解'x的函数值比前一个解的函数值更小，则接受')1(xx=作为一个新解。否则以概率0))0(()'(Txfxfe−−接受'x作为一个新解。泛泛地说，对于某一个温度iT和该优化问题的一个解)(kx，可以生成'x。接受'x作为下一个新解)1(+kx的概率为：⎪⎩⎪⎨⎧=→−−其它若))(()'(0))(()'(1)')((TkxfxfekxfxfxkxP（1）在温度iT下，经过很多次的转移之后，降低温度iT，得到iiTT+1。在1+iT下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。昀终的解是对该问题寻优的结果。我们注意到，在每个iT下，所得到的一个新状态)1(+kx完全依赖于前一个状态)(kx，可以和前面的状态)1(,),0(−kxxL无关，因此这是一个马尔可夫过程。使用马尔可夫过程对上述模拟退火的步骤进行分析，结果表明：从任何一个状态)(kx生成'x的概率，在))((kxN中是均匀分布的，且新状态'x被接受的概率满足式（1），那么经过有限次的转换，在温度iT下的平衡态ix的分布由下式给出：-273-∑∈−−=SjTxfTxfiiiiieeTP)()()(（2）当温度T降为0时，ix的分布为：⎪⎩⎪⎨⎧∈=其它若0||1minmin*SxSPii并且1min*=∑∈SxiiP这说明如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局昀优解将以概率1被找到。因此可以说模拟退火算法可以找到全局昀优解。在模拟退火算法中应注意以下问题：（1）理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能昀终得不到全局昀优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。（2）要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续m次的转换过程没有使状态发生变化时结束该温度下的状态转换。昀终温度的确定可以提前定为一个较小的值eT，或连续几个温度下转换过程没有使状态发生变化算法就结束。（3）选择初始温度和确定某个可行解的邻域的方法也要恰当。1.2应用举例例已知敌方100个目标的经度、纬度如表1所示。表1经度和纬度数据表经度纬度经度纬度经度纬度经度纬度53.712115.304651.17580.032246.325328.275330.33136.934856.543221.418810.819816.252922.789123.104510.158412.481920.105015.45621.94510.205726.495122.122131.48478.964026.241818.176044.035613.540128.983625.987938.472220.173128.269429.001132.19105.869936.486329.72840.971828.14778.958624.663516.561823.614310.559715.117850.211110.29448.15199.532522.107518.55690.121518.872648.207716.888931.949917.63090.77320.465647.413423.778341.86713.566743.54743.906153.352426.725630.816513.459527.71335.070623.92227.630651.961222.851112.793815.73074.95688.366921.505124.090915.254827.21116.20705.144249.243016.704417.116820.035434.168822.75719.44023.920011.581214.567752.11810.40889.555911.421924.45096.563426.721328.566737.584816.847435.66199.933324.46543.16440.77756.957614.470313.636819.866015.12243.16164.242818.524514.359858.684927.148539.516816.937156.508913.709052.521115.795738.43008.464851.818123.01598.998323.644050.115623.781613.79091.951034.057423.396023.06248.431919.98575.7902-274-40.880114.297858.828914.522918.66356.743652.842327.288039.949429.511447.509924.066410.112127.266228.781227.66598.083127.67059.155614.130453.79890.219933.64900.39801.349616.835949.98166.082819.363517.662236.954523.026515.732019.569711.511817.388444.039816.263539.713928.42036.990923.180438.339219.995024.654319.605736.998024.39924.15913.185340.140020.303023.98769.403041.108427.7149我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，昀后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵102102)(×=ijdD，其中ijd表示表示ji,两点的距离，102,,2,1,L=ji，这里D为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个昀短路径。上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设BA,两点的地理坐标分别为),(11yx，),(22yx，过BA,两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则BA,两点的直角坐标分别为：)sin,cossin,coscos(11111yRyxRyxRA)sin,cossin,coscos(22222yRyxRyxRB其中6370=R为地球半径。BA,两点的实际距离⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛⋅⋅=OBOBRdOAOAarccos，化简得]sinsincoscos)(arccos[cos212121yyyyxxRd+−=。求解的模拟退火算法描述如下：（1）解空间解空间S可表为{102,101,,2,1L}的所有固定起点和终点的循环排列集合，即}102,}101,,3,2{),,(,1|),,{(102101211021===ππππππ的循环排列为LLLS其中每一个循环排列表示侦察100个目标的一个回路，ji=π表示在第i次侦察j点，初始解可选为)102,,2,1(L，本文中我们使用MonteCarlo方法求得一个较好的初始解。（2）目标函数此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求∑=+=1011211102),,,(miniiidfπππππL而一次迭代由下列三步构成：（3）新解的产生①2变换法-275-任选序号vu,（vu）交换u与v之间的顺序，此时的新路径为：10211111ππππππππLLL+++−vuuvvu②3变换法任选序号vu,和w，将u和v之间的路径插到w之后，对应的新路径为（设wvu）1021111ππππππππLLLL++−wvuwvu（4）代价函数差对于2变换法，路径差可表示为)()(1111+−+−+−+=Δvvuuvuvuddddfππππππππ（5