模擬退火演算法並行機調度問題

① 退火演算法的應用領域及示例

作為模擬退火演算法應用，討論旅行商問題(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演算法在模型修改、接受概率以及降溫曲線上的改進使得模擬退火演算法收斂速度大大加快。後人在此基礎上還有很多的改進，讀者可以參考相關文獻。

模擬退火法的優點:由於不需要計算偏導數矩陣，不需要解線性方程組(當然正演計算的除外)，結構簡單，易於編程;此外，由於它搜索范圍大，能接受較差模型，因此易於達到全局極小。缺點:隨機搜索，計算量巨大，往往要計算成百上千次正演，這與前面的最小二乘法十幾次的正演計算相比反演時間太長，因此一般應用在一維反演之中，在二維、三維等高維反演中應用較少。