模拟退火算法并行机调度问题

① 退火算法的应用领域及示例

作为模拟退火算法应用，讨论旅行商问题(Travelling Salesman Problem,简记为TSP)：设有n个城市，用数码1,…,n代表。城市i和城市j之间的距离为d(i,j) i,j=1,…,n．TSP问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最短.。
求解TSP的模拟退火算法模型可描述如下：
解空间 解空间S是遍访每个城市恰好一次的所有回路，是{1,……,n}的所有循环排列的集合,S中的成员记为(w1,w2,……,wn)，并记wn+1= w1。初始解可选为(1,……,n)
目标函数 此时的目标函数即为访问所有城市的路径总长度或称为代价函数：
我们要求此代价函数的最小值。
新解的产生 随机产生1和n之间的两相异数k和m,
若k<m,则将
(w1,w2,…,wk,wk+1,…,wm,…,wn)
变为：
(w1,w2,…,wm,wm-1,…,wk+1,wk,…,wn).
如果是k>m,则将
(w1,w2,…,wm,wm+1,…,wk,…,wn)
变为：
(wm,wm-1,…,w1,wm+1,…,wk-1,wn,wn-1,…,wk).
上述变换方法可简单说成是“逆转中间或者逆转两端”。
也可以采用其他的变换方法，有些变换有独特的优越性，有时也将它们交替使用，得到一种更好方法。
代价函数差 设将(w1,w2,……,wn)变换为(u1,u2,……,un)，则代价函数差为：
根据上述分析，可写出用模拟退火算法求解TSP问题的伪程序：
Procere TSPSA:
begin
init-of-T; { T为初始温度}
S={1,……,n}; {S为初始值}
termination=false;
while termination=false
begin
for i=1 to L do
begin
generate(S′form S); { 从当前回路S产生新回路S′}
Δt:=f(S′））-f(S);{f(S)为路径总长}
IF（Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1])
S=S′；
IF the-halt-condition-is-TRUE THEN
termination=true;
End;
T_lower;
End;
End
模拟退火算法的应用很广泛，可以较高的效率求解最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)、调度问题(Scheling Problem)等等。 模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点：
⑴ 温度T的初始值设置问题。
温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。
⑵ 退火速度问题。
模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索（退火）是相当必要的，但这需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。
⑶ 温度管理问题。
温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：
T(t+1)=k×T(t)
式中k为正的略小于1.00的常数,t为降温的次数 优点：计算过程简单，通用，鲁棒性强，适用于并行处理，可用于求解复杂的非线性优化问题。
缺点：收敛速度慢，执行时间长，算法性能与初始值有关及参数敏感等缺点。
经典模拟退火算法的缺点：
⑴如果降温过程足够缓慢，多得到的解的性能会比较好，但与此相对的是收敛速度太慢；
⑵如果降温过程过快，很可能得不到全局最优解。
􀂄 模拟退火算法的改进
⑴ 设计合适的状态产生函数，使其根据搜索进程的需要
表现出状态的全空间分散性或局部区域性。
⑵ 设计高效的退火策略。
⑶ 避免状态的迂回搜索。
⑷ 采用并行搜索结构。
⑸ 为避免陷入局部极小，改进对温度的控制方式
⑹ 选择合适的初始状态。
⑺ 设计合适的算法终止准则。
也可通过增加某些环节而实现对模拟退火算法的改进。
主要的改进方式包括：
⑴ 增加升温或重升温过程。在算法进程的适当时机，将温度适当提高，从而可激活各状态的接受概率，以调整搜索进程中的当前状态，避免算法在局部极小解处停滞不前。
⑵ 增加记忆功能。为避免搜索过程中由于执行概率接受环节而遗失当前遇到的最优解，可通过增加存储环节，将一些在这之前好的态记忆下来。
⑶ 增加补充搜索过程。即在退火过程结束后，以搜索到的最优解为初始状态，再次执行模拟退火过程或局部性搜索。
⑷ 对每一当前状态，采用多次搜索策略，以概率接受区域内的最优状态，而非标准SA的单次比较方式。
⑸ 结合其他搜索机制的算法，如遗传算法、混沌搜索等。
⑹上述各方法的综合应用。

② 2011数学建模国赛B题 求解答

一 问题的重述
110警车在街道上巡逻，既能够对违法犯罪分子起到震慑作用，降低犯罪率，又能够增加市民的安全感，同时也加快了接处警时间，提高了反应时效，为社会和谐提供了有力的保障。
现给出某城市内一区域，其道路数据和地图数据已知，该区域内三个重点部位的坐标分别为：（5112，4806），（9126， 4266），（7434 ，1332）。该区域内共有307个道路交叉口，为简化问题，相邻两个交叉路口之间的道路近似认为是直线，且所有事发现场均在下图的道路上。
该市拟增加一批配备有GPS卫星定位系统及先进通讯设备的110警车。设110警车的平均巡逻速度为20km/h，接警后的平均行驶速度为40km/h。警车配置及巡逻方案要尽量满足以下要求：
D1. 警车在接警后三分钟内赶到现场的比例不低于90％；而赶到重点部位的时间必须在两分钟之内。
D2. 使巡逻效果更显着；
D3. 警车巡逻规律应有一定的隐蔽性。
现在我们需要解决以下几个问题：
一. 若要求满足D1,该区最少需要配置多少辆警车巡逻？
二. 请给出评价巡逻效果显着程度的有关指标。
三．请给出满足D1且尽量满足D2条件的警车巡逻方案及其评价指标值。
四. 在第三问的基础上，再考虑D3条件，给出你们的警车巡逻方案及其评价指标值。
五．如果该区域仅配置10辆警车，应如何制定巡逻方案，使D1、D2尽量得到满足？
六. 若警车接警后的平均行驶速度提高到50km/h，回答问题三。
七. 你们认为还有哪些因素、哪些情况需要考虑？给出你们相应的解决方案。

二 问题分析
本题为城区道路网络中警车配置及巡逻问题。在进行警车配置时，首先要考虑警车在接警后在规定时间内赶到现场的比例，在此条件下，以车数最少为目标，建模、求解；在制定巡逻方案时，要考虑巡逻的效果及隐蔽性问题。
问题一只要求满足D1，求最少的警车配置数，可以认为警车是不动的，在三分钟或两分钟内它能到达的区域就是它的覆盖范围。据此，在满足所有街道的覆盖率不低于90%的条件下，寻找最优解。
问题二要评价巡逻效果，有两个方面需要考虑：一是巡逻的全面性，即经过一段时间后警车走过的街道数占总街道数的比例；二是巡逻的不均匀性，即经过一段时间后警车经过每一条街道的次数相差不大，用方差来衡量。
问题三是在满足D1的条件上尽量满足问题二所给的指标，并给出评价方案的指标。首先找到一组满足D1的各警车位置，然后在和各警车位置相连的点中随机寻找一个点，判断新的点是否满足D1，如果满足则警车行驶到该点，否则重新寻找，直到满足为止。一段时间后统计所有车走过的点数及每个点被走过的次数，用问题二给出的两个指标进行评价。综合两个指标，可判断此路径的好坏，重复这个过程，直到综合评价指标达到一个满意的值为止。
问题四增加了隐蔽性要求，首先给出评价隐蔽性的指标，隐蔽性可用路线的随机性来评价，将它加入到问题三的模型中去进行求解。
问题五限制警车数量为10，要综合考虑D1、D2，先分配这10辆车使道路的覆盖率最高，然后按照问题三的步骤进行求解，其中每一步对D1的判断只需使道路的覆盖率尽量高即可。
问题六同问题三，只需将车速改为50km/h即可。

三 模型的假设

1. 警车都在路上巡逻，巡警去处理案件的时间不考虑；
2. 所有事发现场都在道路上，案件在道路上任一点是等概率发生的；
3. 警车初始停靠点是随机的，但尽量让它们分散分布，一辆警车管辖一个分区；
4. 假定各个划分区域内，较短时间内，最多会发生一个案件；
5. 假设区域内的每条道路都是双行线，不考虑转弯对结果造成的影响；
6. 如果重点部位不在道路上的，假设这些重点部位在离它们最近的道路上；
7. 图中水域对巡逻方案没有影响。

四 符号说明

表示警车数目
表示警车初始停靠点到各道路的最短距离
表示整个区域的总道路长度
表示不能在3分钟内到达的区域的道路的长度
表示非重点部位的警车在3分钟内不能到达现场的比例
表示三分钟内能从接警位置赶到事发现场的最大距离是
表示整个区域总的离散点个数
表示第 区内的节点个数
表示区内调整函数
表示模拟退火的时间，表征温度值
表示区间调整函数
表示全面性指标
表示不均匀性指标
表示综合评价指标
表示第 辆车经过每条道路的次数
表示整个区域每条道路经过的平均次数

五 模型的建立与算法的设计
5.1 满足D1时，该区所需要配置的最少警车数目和巡逻方案
5.1.1 满足D1条件时，区域最少警车的规律
题目要求警车的配置和巡逻方案满足D1要求时，整个区域所需要配置的警车数目最少。由假设可知警车都在道路上，且所有事发现场也都在道路上，但区域内总的道路长度是个定值的；警车在接警后赶到事发现场有时间限制和概率限制：三分钟内赶到普通区域案发现场的比例不低于90％，而赶到重点部位的时间必须控制在两分钟之内。由此可知每辆警车的管辖范围不会很大，于是考虑将整个区域分成若干个分区，每辆警车管辖一个分区域。
由上面的分析，求解整个区域的警车数目最少这个问题可转化为求解每一辆警车所能管辖的街道范围尽量的大。于是我们寻找出使每辆警车管辖的范围尽量大的规律。为了简化问题，我们不考虑赶到现场的90%的几率的限制，仅对警车能在三分钟内赶到事发现场的情况作定性分析，其分析示意图如图1所示。警车的初始停靠位置是随机的分布在道路上的任一节点上，我们假设一辆警车停靠在A点上。

图1 一辆警车管辖范围分析示意图

由于警车的平均巡逻速度为20km/h，接警后的平均行驶速度为40km/h，由于距离信息比较容易得到，于是我们将时间限制转化为距离限制，这样便于分析和求解。当警车接警后，在三分钟内能从接警位置赶到事发现场的最大距离是 ，其中 。
如图1所示，我们设警车初始停靠位置在A点，A点是道路1，2，3，4的道路交叉口。我们仅以警车在道路1巡逻为例来进行分析，警车以 的速度在道路1上A到 点之间巡逻， 与初始停靠点A的距离为 。由于案件有可能在道路上任一点发生，当警车巡逻到A点时，若案发现场在道路2，3，4上发生时，警车以40km/h的速度向事发现场行驶，警车能在三分钟内从 点赶到现场的最大距离为 。如果警车在道路1上继续向前行驶，则该警车能在三分钟内赶到现场的距离继续缩小，当警车从初始点向A点行驶但没有达到 点时，此时该警车的最大管辖范围比警车到达 点时的最大管辖范围大。为了使警车的管辖范围尽量大，警车的巡逻范围越小越好，当 时，即警车在初始停靠点静止不动时，警车的管辖范围达到最大值 。
图1所分析的是特殊的情况，道路1，2，3，4对称分布，现在我们来对一般的情况进行分析，如图2所示。

图2.1 图2.2
图2 一辆警车最大管辖范围分析示意图

图2.1所示的情况是道路分布不对称，与图1相比，图2.1所示的道路方向和角度都发生了改变，图2.3中的情形更为复杂。参照对图1的分析方法，我们分析这两种情形下，警车巡逻时能在三分钟内赶到现场的最大距离的规律，我们只分析图2.2的情况，道路1，2，3，4，5相交于点C,同时道路1与道路6也有个道路交叉口D， 由于警车巡逻时是在道路上行驶的，行走的路线是分段直线,并不影响路径的长度，所以当警车巡逻到距离初始停靠点C点 远处的D，此时若有案件发生时，该警车要在三分钟内能赶到现场处理案件，最大行驶距离在 之内，如果警车在道路1上继续向前行驶，则该警车能在三分钟内赶到现场的距离继续缩小，当警车没有行驶到D点时，此时该警车的最大管辖范围比 大，为了使警车的管辖范围尽量大，警车的巡逻范围越小越好。当 时，即警车静止不动时，一辆警车的管辖范围能达到最大值。
以上分析的仅作定性的分析，对于三个重点部位也可以同理分析，所得的结论是一致的，以上的分析没有考虑到90%的到达几率限制，但在设计算法需要充分考虑。
综上所述，当警车静止在初始停靠点时，在三分钟时间限制内，警车能从初始停靠点赶到事发现场的最大距离为 。

5.1.2 将道路离散化
由于事发现场是等概率地分布在道路上的，由区域地图可以发现，整个区域中的道路长度不均，为了使计算结果更加精确，可将这些道路离散化。只要选取合适的离散方案，就能使警车在经过道路上的离散的点时就相当于经过了这条道路。这样，不论是求解警车初始停靠点还求解警车赶到事发现场所经过的道路时，所计算得的的结果显然比仅考虑整条道路的叉路口要精确得多。
区域中共有307个道路交叉口，458条道路。我们采用线性插值方法对道路进行离散化，以 的速度行走一分钟的距离作为步长，一分钟时间的选择是参照问题三的结果要求来设定的，步长 。用线性插值的方法，从道路的一个方向进行线性插值，实现将每条道路离散化的目标，考虑到有些道路不是 的整数倍，我们就一般情况进行讨论，其分析示意图如图3所示。道路AB长度为 个 与 长度的和，为了更精确处理CB段道路，那么就要考虑在CB之间是否要插入一个新的点， 根据 的长度不同，其对应的处理方式也有所不同。

图3 道路离散化分析示意图

引进临界指数 ，选取 大小的准则是使尽量离散化后警车等效的平均巡逻速度和题目给定的速度（ ）的差值尽量小，经过计算得 时，不再插入新的坐标点时能使整个区域的道路离散效果较好。此时，将CB段长度设定为 处理，于是离散后的AB道路长度会比实际长度短些；当 时，需要在两个点之间再插入一点，因为这样处理能使整个区域的整体道路的离散化效果比较理想。如图3所示，在C与B间再插入新的坐标点，插入的位置在距C点 的D点处,这样处理后所得的道路长度比实际长度长了 。采用这样的方法进行线性插值，我们使用MATLAB编程实现对整个区域道路的离散，所得的离散结果如图4所示，离散后共得到762个节点，比原始数据多了455个节点，离散后的节点数据见附件中的“newpoint.txt”。

图4 整个区域离散结果图

采用这种插值方法道路离散后，将直线上的无穷多个点转化有限个点，便于分析问题和实现相应的算法，由图4可知，所取得的整体离散效果还是比较理想的。

5.1.3 分区域求解警车数目的算法设计
考虑到警车配置和巡逻方案需要满足：警车在接警后三分钟内赶到普通部位案发现场的比例不低于90%，赶到重点部位必须控制在两分钟之内的要求。设计算法的目标就是求解出在满足D1情况下，总的警车数目最小，即每个区域都尽可能多地覆盖道路节点。由于警车的初始位置是未知的，我们可设警车初始停靠点在道路上的任一点，即分布在图4所示的762个离散点中的某些点节点上，总体思路是让每两辆车之间尽量分散地分布，一辆警车管辖一个分区，用这些分区覆盖整个区域。
于是我们设计算法1，步骤如下所示：

Step1：将整个区域预分配为 个分区，每个分区分配一辆警车，警车的初始停靠位置设在预分配区中心的道路节点上，若区域的中心不在道路节点上，则将警车放在离中心最近的道路节点上；
Step2：统计分区不能覆盖的节点，调整警车的初始停靠点，使分区覆盖尽可能多的道路节点，调整分为区内调整和区间调整方案：（1）区内调整按照模拟退火思想构造的函数，在区间调整调整车辆初始点的位置（后文中有详细说明），当分区内节点数较多时，调整的概率小些，分区内节点数较少时，调整的概率大些，（2）当区域中存在未被覆盖的节点或节点群（大于等于三个节点集中在一个范围内）时，将警车初始位置的调整方向为朝着这些未被覆盖的节点按一定的规则（在算法说明中有详细叙述）移动，同时要保证 3个重点部位能在2分钟之内100%到达；
Step3：用Floyd算法计算出警车初始停靠点到周边各道路节点的最短距离 ；
Step4：以 个划分区域未覆盖的总的道路长度 与整个区域的道路总长度 的比值 来表示警车不能3分钟内到达现场的概率；
Step5：模拟足够多的次数，若 ，将车辆数 减1，跳转到Step1；
Step6：计算结束后，比较当 时所对应的 值， 当 取得最小值时，记录此时的区域划分方案， 即为最少的警车数。

对算法的几点说明：
（1）该算法所取的车辆数 是由多到少进行计算的， 初始值设为20，这个值的选取是根据区域图估算的。
(2)预分区的优点在于使警车的初始位置尽可能均匀地分散分布，警车的初始停靠点在一个分区的中心点附近寻找得到，比起在整个区域随机生成停靠点，计算效率明显得到提高。
预分配之后，需要对整个区域不断地进行调整，调整时需要考虑调整方向和 调整概率。
警车调整借鉴的是模拟退火算法的方法，为了使分区内包含道路节点数较多的分区的初始停车点调整的概率小些，而分区内包含道路节点数的少的分区内的初始停车点调整的概率大些，我们构造了一个调整概率函数 ,
（1）
（1）式中， 均为常数， 为整个区域车辆数， 为第 分区内覆盖的节点数， 为时间，同时 也能表征模拟退火的温度变化情况：初始温度较高，区域调整速度较快，随着时间的增加，温度不断下降，区域调整速度逐渐变慢，这个调整速度变化也是比较符合实际情况的。
由式（1）可以得出调整概率函数 ，假设在相同的温度 （时间）的条件下，由于总的车辆数目 是定值，当 时，即第 分区内的节点数大于第 分区的节点数时，分区 调整的概率大些，分区 的调整概率小些。分析其原因：当分区内包含了较多的节点个数时，该分区的警车初始停靠位置选取地比较合适了，而当分区内包含的道路节点数较少时，说明警车的初始停靠位置没有选好，需要更大概率的调整，这样的结论也是比较客观的。
对于所有分区外未被覆盖的道路节点和很多节点（称之为节点群），用来调整警车位置迁移的方向，其分析示意图如图5所示。调整方案目标是使未被覆盖的节点数尽量的少。在设计调整方向函数时，需要考虑：（1）节点群内节点的数目；（2）警车距离节点群的位置。优先考虑距离，所以在公式（2）中，用距离的平方来描述调整方向函数。
由于某一个区域范围内的未被覆盖节点数，整个区域未被覆盖的节点总数，分区域与未被覆盖的节点或节点群的距离等几个因素会影响到调整的方案，所以要综合考虑这些因素。于是设计了区间调整函数 ，
式中， 表示第 个分区内未被覆盖的节点数， 表示第 分区域与未被覆盖的节点或节点群的距离， 表示未被覆盖的节点和节点群个数。
现在简要分析第 分区按区间调整函数的调整方案，当某两节点群 的节点数目相等，但是距离不等时，如 ，由区间调整公式可知，该区间向节点群 方向调整。当某个分区与两个节点群的距离相等，但节点群的内节点个数不相等，如 时，由（4）可知，该分区域会想节点群 方向调整。
注意在整个调整过程中，调整几率控制是否调整，调整方向函数控制调整的方向，寻找在这种调整方案下的最优结果。

图5 调整分区域示意图

（3）在step3中，使用Floyd算法计算出警车初始停靠点到周边各节点的最短距离 ,目的是当区域内有情况发生时，警车能在要求的时间限制内到达现场。
（4）为求出较优的警车停靠点，采用模拟退火算法，算出局部最优的方案。
5.1.4 警车的配置和巡逻方案

使用MATLAB编程实现算法1得到，整个区域配备13辆警车，这些警车静止在初始停靠点时，能满足D1要求。警车的初始停靠位置分别为道路交叉节点6，25，30，37，82，84，110，111，126，214，253，258，278处。每个警车所管辖的交叉点（原始的交叉节点）如图6所示，求解的分区结果见附录所示。
图6 满足D1条件下的区分划分图

13个分区共覆盖了252个交叉点，另外的55个原始交叉点没有被这些分区域覆盖：137，138，151，159，167，168，170，174，175，186，188，189，211，215，226，242，255，260，261，262，263，267，270，271，272，275，282，283 ，284，287，288，289，292，296，297，299，304，305，307。在这种分区方案下，这些点中，每两个相连的点间的道路离散值长度占整个区域总的长度的比值为 。因此，在整个区域配置13辆警车，每个警车在初始停靠点静止不动，当有案件发生时，离案发现场最近的警车从初始停靠点赶到现场。

5.2 评价巡逻效果显着的指标
110警车在街道上巡逻是目的是为了对违法犯罪分子起到震慑作用，降低犯罪率，又能够增加市民的安全感，同时还加快了接处警（接受报警并赶往现场处理事件）时间，提高了反应时效，为社会和谐提供了有力的保障。巡警在城市繁华街道、公共场所执行巡逻任务, 维护治安, 服务群众, 可以得良好的社会效应[1]。
在整个区域中，由于案发现场都在道路上，道路上的每一点都是等概率发生的，因此警车巡逻的面越广，所巡逻的街道数目越多，警车的巡逻效果就越好，对违法犯罪分子就越有威慑力，警车也能更及时地处理案件。
我们采用全面性 来衡量巡逻的效果显着性，即用警车巡逻所经过的街道节点数占区域总节点数的比值。当警车重复经过同一条街道同一个离散点时， 仅记录一次。
（3）
式中， 表示警车经过的离散点数， 代表整个区域总的离散点数。 值越大，表明警车所经过的街道数目越多，所取得的效果越显着。
同时考虑到在巡逻过程中可能会出现这样的情况：在相同的时段内，警车会多次巡逻部分街道，而一些街道却很少巡逻甚至没有警车到达，这样会造成一些巡逻盲区。分布很不均衡。这样就可能出现巡逻密度大的街道上的违法犯罪分子不敢在街道上作案，而流窜到巡逻密度稀疏的街道上作案，因此在相同的警车数目条件下，密度不均衡的巡逻方式的巡逻效果的效果较差，而密度较均衡的巡逻方式所取得的巡逻效果会更好些。我们引入一个巡逻的不均匀度 来衡量巡逻效果的显着性，考虑到方差能表示不均衡度，于是我们用方差的大小来表征不均衡，方差越大，巡逻密度越不均衡，所取得的巡逻效果越差。
（4）

问题1所给出的满足D1条件下的警车数目为13辆，这时每辆警车在初始停靠点静止不动，只有该管辖区域内发生了案件时，警车才从初始停靠点赶到案发现场处理案件。当警车在巡逻状态时，所需要考虑的问题就更复杂一些，如当节点运动时，警车还能否达到D1的要求，警车的运动方向如何等问题，但基本算法思想与问题1类似，所得的算法2的框图如图7所示，
为了简化问题，我们假设各分区警车的巡逻时候，尽量保证所有的警车的行驶方向相一致，且警车都走双行道，即当警车走到某个节点后，它们又同时返回初始停靠点，警车的行驶方向有四种方式，如6所示。
在图6中，数字1代表走巡逻走的第一步，2表示朝1的巡逻方向相反的方向巡逻。在具体程序实现时，四种巡逻方向任意选择，但是尽量保证所有的警车向同一个方向巡逻。

图6 各警车巡逻方向图

我们用MATLAB编程对这种巡逻方式进行计算，所得的车辆数目为18辆，综合评价指标为 ，其结果巡逻方案见附件中的“1193402-Result3.txt”所示。

5.4 在满足问题三的基础上讨论D3条件，警车的巡逻方案和评价指标
巡逻的隐蔽性体现在警车的巡逻路线和时间没有明显的规律，主要目的是让违法犯罪分子无可乘之机，防止他们在非巡逻时间实施违法犯罪活动，危害人民的生命和财产安全。
为了使巡逻的规律具有隐蔽性，这就需要警车在巡逻时至少具有两条不同的路线，时间最好也是不相同的。因此，考虑到隐蔽性时，只需要在问题2的基础上加上一个随机过程即可。对于其评价指标，由于警车有几条可选的巡逻路线，当相同的路线在同一时间内重复出现时，重新将所设定的方案再执行一遍，我们用这个时间间隔来衡量隐蔽性的程度，当循环周期 越大，表明可选的巡逻方案越多，其规律就越具有隐蔽性，而循环周期 越小时，表明巡逻方案比较少，其隐蔽性较差。在巡逻状态时，最差的隐蔽性巡逻方案是巡逻方案只有一个，并且时间固定，这样的巡逻方案没有任何隐蔽性可言。

5.5 整个区域为10辆车时的巡逻方案
由第三问的结果可知，10辆车的数量是不能把整个区域完全覆盖的，其算法与算法2类似，不同的是此时车的数目已经固定了，要求使D1,D2尽量大的满足，我们求得的评价指标值为 ，所得的巡逻方案见附件中的“1193402-Result5.txt”所示。

5.6 平均行驶速度提高到 时的巡逻方式和评价指标值
问题六的分析方法与具体实现与问题三一致，但是警车的接警后的平均速度由原来的 提高到 ，于是各分区的覆盖范围也增大了，将数值带入问题3的算法中求解， 计算得的指标值为 ，其巡逻方案见附件中的“1193402-Result6.txt”所示。

图7 算法2框图

六 模型的分析和评价

在求解满足D1的条件下，整个区域需要配备多少辆警车问题中，采用分区巡逻的思想，先分析能使各区管辖范围达到最大值时的规律，由特殊到一般层层进行分析，逻辑严密，结果合理。
在求解区域和警车数目时，在初步设定警车停靠点位置的基础上，用模拟退火算法思路构造函数 来确定调整的概率大小，综合考虑了影响区间调整的因素后构造了 函数来确定分区的调整方向，当分区按照这两个调整函数进行调整时，各分区能管辖尽可能多的道路节点，所取得效果也比较理想。

参 考 文 献
[1]中小城市警察巡逻勤务方式的探讨，俞详，江苏公安专科学校学报，1998年第1期
[2]Matlab7.0从入门到精通，求是科技，人民邮电出版社；
[3]不确定车数的随机车辆路径问题模型及算法，运怀立等，工业工程，第10卷第3期，2005年5月；
[4]随机交通分配中的有效路径的确定方法，李志纯等，交通运输系统工程与信息，第3卷第1期，2003年2月。

③ 模拟退火算法的模拟退火算法的原理

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schele)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。 1模拟退火算法可以分解为解空间、目标函数和初始解三部分。
2模拟退火的基本思想:
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L
(2) 对k=1，……，L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数
(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。 模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。 s:=s0;e:=E(s)//设定目前状态为s0，其能量E(s0)k:=0//评估次数kwhilek<kmaxande>emax//若还有时间（评估次数k还不到kmax）且结果还不够好（能量e不够低）则：sn:=neighbour(s)//随机选取一临近状态snen:=Esn)//sn的能量为E(sn)ifrandom()<P(e,en,temp(k/kmax))then//决定是否移至临近状态sns:=sn;e:=en//移至临近状态snk:=k+1//评估完成，次数k加一returns//回转状态s

④ 模拟退火法^[1，]

模拟退火算法最早在1953年由 Metropolis等人提出。在地球物理中的最早应用是Rothman在1983年利用模拟退火算法处理地震资料的剩余静校正。模拟退火法也是类似于蒙特卡洛法的随机搜索方法。但是在产生模型的过程中引入一些规则，能有效地加快搜索速度，有时又称这类方法为启发式蒙特卡洛法。

模拟退火法概念源于统计物理学，是模拟固体熔化状态逐渐缓慢冷却最终达到能量最小的结晶状态的物理过程。对于一个熔化的金属，当处于某个温度的热平衡状态时，它的每一个分子都有它可能所处的状态，有些分子可能能量高一些，有些分子可能能量低一些，分子处于何种状态的概率由分子所具有的能量决定。设分子所有可能的能级总数为n(微观粒子的能量都是量子化的，不连续的)，则分子处于某种状态的概率满足玻尔兹曼概率分布:

地球物理反演教程

其中:Ei为第i个分子的能量;K为玻尔兹曼常数;T为绝对温度;n为分子所有可能的能级总数，分母称为配分因子;p_i为第i个分子处于能量E_i的概率。

如果把地球物理反演的模型向量看作分子，把目标函数看作分子的能量，把目标函数的极小值看成分子冷却结晶的最小能量，反演问题(最优化问题)可以模拟式(8.11)金属退火的过程，通过缓慢地减小温度进行反演，使目标函数(能量)逐渐达到极小值，这时所对应的模型(分子状态)就是反演结果。

为了改善于蒙特卡洛法的随机搜索方法，1953年 Metropolis等人在产生模型的过程中引入Metropolis接受准则，模型产生并不是完全随机，而是以前一个模型为基础随机产生。对能量减小的模型完全接受，对能量增加的模型按一定的概率接受，这样能有效地加快搜索速度，同时又有可能跳出局部极小值。具体如下:

设原来模型向量为m_i，新的模型为m_i+1(在m_i基础上随机修改产生)，各自的能量(目标函数)为E(m_i)和E(m_i+1)。如果E(m_i+1)＜E(m_i)，则目标函数在减小，新模型可以接受。如果E(m_i+1)＞E(m_i)，则目标函数在增加，按照一定概率来确定是否接受新的模型。具体规则见式(8.12):

E(m_i+1)＜E(m_i) 完全接受m_i+1为新模型

地球物理反演教程

式(8.12)就是Metropolis接受准则。它使得反演过程可以接受使目标函数增加的模型，因此也就使得模拟退火法有可能跳出局部极小，收敛于全局极小值点。由于玻尔兹曼常数K只是起到尺度因子的作用，在实际计算中K可取为1来简化公式。从式(8.12)可以看出，当温度较低时，p_i+1/p_i较小，因此接受使能量增加的新模型的可能性较小。而一般温度较低时，目标函数较小，模型比较靠近真实模型，这时基本上只接受使目标函数减小的模型，使模型尽快收敛于极小值点。

在模拟退火反演中，要求温度T随着迭代次数的增加而缓慢降温。常用的温度函数有两种。

(1)指数下降型:

T_k=T₀·exp(-ck^1/N) (8.13)

式中:k为迭代次数;c为衰减因子;N为模型参数的个数;T₀为初始温度。上式也可以改写为

地球物理反演教程

通常选择0.7≤α≤1。在实际应用中可采用0.5或1代替式(8.14)的1/N。图8.4(a)为指数降温曲线。采用参数为:T₀=200℃，α=0.99，1/N=0.9。

(2)双曲线下降型:

T=T₀α^k (8.15)

式中:T₀为初始温度;k为迭代次数;α为衰减因子，通常取0.99。初始温度T₀不能取得太高，否则增加计算时间浪费机时;T₀也不能太低，否则模型选取不能遍及整个模型空间，只是在初始模型附近选取，不能进行全局寻优。所以T₀的确定只有通过实验计算得到。图8.4(b)为双曲线降温曲线。采用参数为:T₀=200℃，α=0.99。从图8.4可以看出通过对不同温度曲线和相关参数进行选择，可以控制温度下降的方式和速度。

图8.4 模拟退火法降温曲线

模拟退火法主要有三种:

(1)MSA算法(Metropolis Simulated Annealing);

(2)HBSA算法(Heat Bath Simulated Annealing);

(3)VFSA算法(Very Fast Simulated Annealing)。

图8.5 模拟退火MSA算法程序流程图

前面介绍的利用 Metropolis接受准则的算法就是经典的模拟退火法。图8.5为模拟退火 MSA算法的程序流程图。从中可以看出 MSA算法有一套模型修改准则，依次改变模型参数，每次改变都是在原来模型基础上改变一个参数，因此容易保持已有搜索成果，持续不断地向目标函数最小值点接近，因此搜索效率比蒙特卡洛法高。此外，MSA算法允许接受使目标函数增加的模型，这样又易于跳出局部极小，达到全局极小。但 MSA算法在任何温度下和蒙特卡洛法一样都是在模型全空间进行搜索，不能根据当前温度和模型减小搜索空间，此外由于模型的修改全凭运气，所以不可能像前面介绍的最小二乘法那样目标函数基本上持续减小，而是呈不规则振荡在宏观上逐渐减小，因此效率较低。

HBSA算法与 MSA算法的不同之处是在模型的修改上。也是首先随机选择一个初始M维模型向量m₀(它具有M个参数);然后限制各个模型参数可能的取值范围，对取值离散化。假设每个模型参数都有N个可能的值，首先固定模型第2个参数m₀(2)直到第M个参数m₀(M)保持不变，只修改第1个参数m₀(1);计算m₀(1)的所有取值时的目标函数，然后按式(8.16)计算“概率”，它就是式(8.11)配分因子取1的公式。即

地球物理反演教程

选择“概率”最大的为模型第1个参数的修改值。照此依次对所有模型参数进行修改完成依次迭代计算。在每次迭代计算中保持温度不变。随着迭代次数增加，温度降低，最终达到稳定状态，获得最小能量解。这种方法的计算由于要计算某个参数的所有可能值，所以计算量也是很大的。

1989年Ingber提出了VFSA算法，由于速度较快，最为常用。它使得模拟退火法从理论走向了实际应用。VFSA算法在流程上与传统的模拟退火法相同，但是在模型修改、接受概率以及降温曲线上有所改进。

(1)模型修改:常规模拟退火法采用高斯随机分布修改模型，在任何温度下都是在模型全空间进行搜索。而Ingber提出采用依赖于温度的似cauchy分布产生新的模型。即

地球物理反演教程

y_i=Tsgn(u-0.5)[(1+1/T^|2u-1|-1](8.18)

其中:m_i为当前模型第i个参数，m'_i为修改后的模型参数;u为[0，1]的随机数;[A_i，B_i]为m_i和m'_i的取值范围;sgn( )为符号函数。

采用以上方式能在高温下进行大范围的搜索，低温时在当前模型附近搜索，而且由于似cauchy分布具有平坦的“尾巴”，使其易于迅速跳出局部极值。这一改进大大加快了模拟退火法的收敛速度。

(2)接收概率:当E(m_i+1)＞E(m_i)时，VFSA算法采用如下概率接受公式:

地球物理反演教程

上式当h→1时变为式(8.12)。h通过实验获得。

(3)降温曲线(退火计划):Ingber在1989年采用式(8.13)得出指数降温曲线。从图8.4可知，温度下降较快。

总之，VFSA算法在模型修改、接受概率以及降温曲线上的改进使得模拟退火算法收敛速度大大加快。后人在此基础上还有很多的改进，读者可以参考相关文献。

模拟退火法的优点:由于不需要计算偏导数矩阵，不需要解线性方程组(当然正演计算的除外)，结构简单，易于编程;此外，由于它搜索范围大，能接受较差模型，因此易于达到全局极小。缺点:随机搜索，计算量巨大，往往要计算成百上千次正演，这与前面的最小二乘法十几次的正演计算相比反演时间太长，因此一般应用在一维反演之中，在二维、三维等高维反演中应用较少。