基于改進模擬退火算法的非線性土方調配方案優(yōu)化方法

范 斌 李遲典 周 誠

(1.武漢地鐵集團有限公司 質量安全監(jiān)察部，武漢 430070；2.華中科技大學 國家數(shù)字建造技術創(chuàng)新中心，武漢 430074；3.華中科技大學 土木與水利工程學院，武漢 430074)

引言

土方工程是建設項目的重要一環(huán)，其順利、穩(wěn)定完工是整個工程按時完成的前提[1]。由于土方工程的復雜性，其成本通常占整個工程成本的較大部分[2]。土方調配方案是影響土方工程成本的重要因素，針對土方調配方案的優(yōu)化對成本的控制有很大作用。

將土方調配過程簡化并歸納出目標函數(shù)與約束函數(shù)，再利用數(shù)學模型對其求解是土方調配優(yōu)化的一般思路。當前研究多集中于線性規(guī)劃模型[3]。其中，Hare考慮路徑可訪問性并利用混合整數(shù)線性規(guī)劃模型進行求解[4]；邢志國利用綜合運距替換平均運距進行模型求解，使得土方調配優(yōu)化結果與實際工程成本相適應[5]；閆若鈺利用Revit對土方量進行計算并利用Dynamo實現(xiàn)土方調配線性規(guī)劃自動求解[6]。除此之外，多目標規(guī)劃與動態(tài)規(guī)劃的研究也較多，Parente依據(jù)成本最低和持續(xù)時間最短構建多目標動態(tài)規(guī)劃模型并利用多目標遺傳算法進行土方調配方案求解[7]；黎天勝以成本最低和環(huán)保為目標構建土方調配模型來實現(xiàn)成本節(jié)約與環(huán)境保護[8]；常峻從實時監(jiān)控系統(tǒng)獲取信息來建立動態(tài)規(guī)劃模型并據(jù)此對土方調配模型進行實時規(guī)劃[9]。

但當前針對非線性土方調配的研究較少。Easa考慮實際土方施工中借方和棄方的單位成本隨借方和棄方量的變化而變化，構建了一個非線性土方調配模型[10]。其他研究大多則是先利用啟發(fā)式算法求解線性調配問題，再提出其用于求解非線性調配問題的可能性。

現(xiàn)實土方調配中存在借方階梯費率、調配量過大導致交通擁堵、調配量不同而采用不同運輸工具導致費率不同等多種非線性因素，單純的線性規(guī)劃模型已不能滿足實際需要，亟需對非線性規(guī)劃模型進行研究。

1 土方調配非線性規(guī)劃模型

土方調配非線性規(guī)劃模型最突出的特點是目標函數(shù)是一個非線性函數(shù)，即土方成本與土方量的函數(shù)關系是非線性的。土方成本的非線性主要是由于在考慮土方調配非線性因素的情況下，土方工程費率不再是一個固定值，而是會隨著土方量的增加而發(fā)生變化的。分析各因素對土方成本的影響原理，在線性成本的基礎上提出非線性土方調配的成本函數(shù)，如式1所示：

式1中，Vij、dTij為由挖方區(qū)i調配到填方區(qū)j的調配體積與轉運運距，VFi、VCi、dCi、dDi為挖方區(qū)i的填方部分的體積、挖方部分體積、收集運距及攤鋪運距；VFj、VCj、dCj、dDj為填方區(qū)j的填方部分的體積、挖方部分體積、收集運距及攤鋪運距；fC、fT、fD、fW、fB為土方調配中的收集費率函數(shù)、轉運費率函數(shù)、攤鋪費率函數(shù)、棄方費率函數(shù)及借方費率函數(shù)；Nc、Nf則分別表示挖方區(qū)數(shù)量、填方區(qū)數(shù)量。費率函數(shù)是非線性土方成本的重要計算因素，一般為調配土方量的函數(shù)。如果費率函數(shù)是一個常數(shù)，則式1可以表示為線性調配成本；如果費率函數(shù)不是常數(shù)，則式1為非線性規(guī)劃的目標函數(shù)。

對于土方調配問題，其約束一般為土方量約束?，F(xiàn)假設土方施工分區(qū)整體土方量為填方大于挖方，需要借土且所有土方材料都能滿足填方需求，不需要進行土方換填等操作。針對上述工程假設，土方量約束函數(shù)用式2到式6表示：

(2)

(3)

(4)

VF=0

(5)

Vij≥0

(6)

其中，VB表示借方體積，VF表示棄方體積，其余字母含義與式1相同。式2到式6的約束條件是假設填方量多余挖方量的情況下得到的，每個公式代表不同的含義。式2表示填方區(qū)收到的土方材料體積要小于或等于自己需求的材料體積，小于部分的體積由之后的借方體積補充； 式3表示挖方區(qū)外運的體積要等于自己多余的體積； 式4表示借方體積的計算，因為假設總挖方量要小于總填方量且所有挖方土壤都可以滿足回填要求，因此借方體積為總填方體積與總挖方體積的差； 式5表示棄方體積的計算，在假設條件下的棄方體積即為零； 式6表示調配土方量要滿足非負性約束。每個棄方區(qū)的棄方量約束和每個借方區(qū)的借方量約束在本文中不做詳細考慮。

通過將約束條件整合到目標函數(shù)中來實現(xiàn)對約束的處理。等式約束通過消除一個參數(shù)來轉化為不等式約束，例如對于式3所示的約束，需要將Vi1作為換元主體換到等式一邊，如式7中的第二個式子所示；由于Vi1是一個非負參數(shù)，因此可以將其轉換為一個不等式，如下列第三個式子所示：

Vi1=Vci-VFi-(Vi2+Vi3+…+ViNf)

(7)

Vci-VFi-(Vi2+Vi3+…+ViNf)≥0

對于不等式約束則通過構建罰函數(shù)來進行轉化。如果解滿足不等式要求則罰函數(shù)應該為0，否則通過罰函數(shù)來對解的越界進行懲罰。本文以二次罰函數(shù)來構建罰函數(shù)：以式2為例，此約束是不等式約束，其罰函數(shù)構造如式8所示，其中σ為常數(shù)，控制罰函數(shù)F的放縮，對于不同的問題有不同的大小，需要根據(jù)實際情況考慮，本文通過多次嘗試后確定σ取值為1 000。

(8)

為了利用啟發(fā)式算法對有約束的非線性規(guī)劃問題求解，需要將成本與約束整合為一個統(tǒng)一的目標函數(shù)，如式9所示。式9一共由兩部分組成：第一部分即為土方成本，第二部分是罰函數(shù)。罰函數(shù)由式2、式3通過處理得到，式4和式5是借方和棄方體積，其成本不隨調配方案改變而改變，不進行轉換； 式6為非負約束，不用轉化，可以直接代入啟發(fā)式算法進行計算。

通過上述論述，可以將非線性土方調配問題轉化為式9所示的非負約束下的函數(shù)最小值問題，然后利用啟發(fā)式算法進行求解。

2 改進模擬退火優(yōu)化算法

模擬退火算法是一種通過模擬固體的退火過程進行優(yōu)化的常用啟發(fā)式算法，其算法流程可以用圖1表示。首先，在解空間種隨機產生一個新解，將這個新解作為一個初始解帶入后續(xù)算法； 通過一定的規(guī)律在舊解的基礎上擾動產生新解并計算新舊兩解目標函數(shù)值的變化量； 如果變化量≤0，則接受新解； 如果變化量>0，則按照式10所示的Metropolis準則判斷是否接受，其中ΔE即為目標函數(shù)值的變化量； 根據(jù)內循環(huán)次數(shù)要求執(zhí)行內循環(huán)并判斷當前溫度下是否滿足了終止條件，如果滿足則輸出最終解。

(10)

模擬退火算法可以跳出局部最優(yōu)獲得全局最優(yōu)解，過程簡單，魯棒性強，在非線性優(yōu)化問題的求解上有較大優(yōu)勢。但該算法也有一些缺點，例如算法運行時間較長，算法運行中可能會導致最優(yōu)解丟失，不同參數(shù)對于算法結果影響較大等問題。

圖1 模擬退火算法流程

針對算法運行時間較長的問題，本文參考L. Ingber提出的方法對降溫函數(shù)與擾動方式進行改進[11]。其降溫函數(shù)改進為式11所示的指數(shù)函數(shù)，其中，T0為初始溫度，k為降溫次數(shù)，α為降溫系數(shù)，一般為[0.7， 1]之間取值，D一般取1，此降溫函數(shù)可顯著提高降溫速度。改進的擾動方式如式12所示，舊解x的i維分量xi的上下限為[mini，maxi]，則其i維分量xi經(jīng)過擾動后產生的新解為xi’，sgn為階躍函數(shù)，T為當前情況下的溫度，ui為[0， 1]上均勻分布的隨機數(shù)。

T=T0αk1/D

(11)

針對算法運行中的最優(yōu)解丟失問題，本文在算法中加入最優(yōu)解保留策略，最優(yōu)解僅記錄整個過程中的最優(yōu)目標函數(shù)值而不參與退火過程。

對于不同參數(shù)對算法結果影響較大，本文嘗試通過對比調參來保證結果的可靠性。

3 算法驗證與對比

3.1 案例簡介

為了驗證本文方法的效果，選用武漢某土方工程進行計算對比與驗證。該工程填挖方量大，土方工期長，土方工程成本大，具有代表性。本文選擇三個不同的算例進行算法分析：根據(jù)整體工程分包分區(qū)情況和施工段劃分情況進行土方工程分區(qū)，具體分區(qū)數(shù)量為6，如圖2所示； 根據(jù)現(xiàn)場實際的土方工程分區(qū)進行分區(qū)調配，具體分區(qū)數(shù)量為13，如圖3所示； 根據(jù)區(qū)域聯(lián)通性和區(qū)域面積進行分區(qū)，具體分區(qū)數(shù)量為29，如圖4所示。其中13分區(qū)情況為實際的分區(qū)情況，故將其作為主要研究場景進行研究并進行算法比選。為了解不同問題規(guī)模下改進算法的求解性能，本文利用改進算法對6分區(qū)問題和29分區(qū)問題進行求解，分析不同的問題規(guī)模下算法的變化。參考實際情況，假設土方調配過程中土方轉運量不大于20 000m3時，收集費率、轉運費率和攤鋪費率為6.65元/(km·m3)； 如果土方轉運量大于20 000m3，土方費率降低為原費率的0.9，即為5.985元/(km·m3)，不考慮其他影響土方成本的線性與非線性因素。

圖2 6分區(qū)方案劃分示意圖

圖3 13分區(qū)方案劃分示意圖

圖4 29分區(qū)方案劃分示意圖

表1 各算法可選參數(shù)

3.2 算法調參及結果展示

分析前人的研究經(jīng)驗，對于本文的土方調配優(yōu)化問題，擬采用改進模擬退火、遺傳算法、粒子群算法和差分進化算法對比求解，各算法的可選參數(shù)組合如表1所示。其中改進模擬退火的冷卻進度表還有一個參數(shù)為終止溫度，這里將終止溫度固定為10-7來減少調參。為了減小算法的隨機性，每個參數(shù)組合均計算三次來減小誤差，參數(shù)組合利用“算法種類(SA/GA/PSO/DE)-參數(shù)一—參數(shù)二—參數(shù)三—參數(shù)四”來表示。

將6分區(qū)調配問題帶入改進模擬退火計算得到1 050個運行結果。因為數(shù)據(jù)量過多，將所有數(shù)據(jù)按照馬爾可夫鏈長度數(shù)N進行分類，以初始溫度和降溫系數(shù)的參數(shù)組合為橫軸，土方成本為縱軸繪制圖形，如圖5所示。由圖可知，最優(yōu)值為1 960 892元，當N大于500 000時，基本所有參數(shù)組合均接近這個最優(yōu)值。例如起始溫度=100，馬爾可夫鏈長度=500 000，降溫系數(shù)=0.75得到的最優(yōu)值就是1 960 892元，此參數(shù)下目標函數(shù)隨時間變化如圖6所示。

圖5 6分區(qū)改進模擬退火算法結果匯總

按照同樣的方法對分區(qū)數(shù)量為13、29的兩種情況求解，得到如下結果：對于13分區(qū)問題，調參結果如圖7所示，改進算法得到的最小成本為1 765 672元，參數(shù)組合SA-1000000-300000-0.85，此參數(shù)組合下的算法隨時間收斂如圖8所示； 對于29分區(qū)，結果如圖9所示，成本的最小值為1 922 424元，是在參數(shù)組合SA-1000000-600000-0.8處得到的，此參數(shù)組合下的算法隨時間收斂如圖10所示。

圖6 6分區(qū)改進模擬退火最優(yōu)解之一的收斂情況

圖7 13分區(qū)改進模擬退火算法結果匯總

圖8 13分區(qū)下改進模擬退火最優(yōu)解收斂過程

圖9 29分區(qū)改進模擬退火算法結果匯總

圖10 29分區(qū)下改進模擬退火最優(yōu)解變化過程

圖11 13分區(qū)下遺傳算法種群數(shù)量和終止代數(shù)對比圖

對于遺傳算法，種群數(shù)量和終止代數(shù)是影響算法運行時間的主要參數(shù)，因此根據(jù)種群數(shù)量和終止代數(shù)的組合可以將算法結果分為18類，選取每類最優(yōu)的參數(shù)組合，得到表3所示的表格。將表3得到的18個數(shù)據(jù)按照種群數(shù)量N分類匯總可以得到如圖11所示的不同參數(shù)對比圖。從圖11可以得出針對13分區(qū)問題，遺傳算法最優(yōu)的參數(shù)組合為GA-3000-3000-0.99-0.001，得到的最優(yōu)解為1 772 355元，最優(yōu)參數(shù)組合的目標函數(shù)隨算法運行時間的圖像，如圖12所示。

圖12 13分區(qū)下遺傳算法最優(yōu)解收斂過程

表3 13分區(qū)遺傳算法交叉概率和變異概率對比表

對于粒子群算法，種群數(shù)量和終止代數(shù)是影響算法運行時間的主要參數(shù)，根據(jù)種群數(shù)量和終止代數(shù)的組合可以將算法結果分為9類，得到表4和圖13所示的不同參數(shù)對比圖。從圖13可以得出，針對13分區(qū)問題，粒子群算法最優(yōu)的參數(shù)組合為GA-3000-12000-1.1-1，得到的最優(yōu)解為1 858 936元，最優(yōu)參數(shù)組合的目標函數(shù)隨算法運行時間的圖像如圖14所示。

表4 13分區(qū)粒子群算法慣性權重和學習因子對比表

圖13 13分區(qū)下粒子群算法種群數(shù)量和終止代數(shù)對比圖

圖14 13分區(qū)粒子群算法最優(yōu)解收斂過程

同理，對于差分進化算法，可以得到表5和圖15。最優(yōu)參數(shù)組合的目標函數(shù)隨算法運行時間變化的圖像如圖16所示。

表5 13分區(qū)差分進化算法交叉概率和放縮因子對比表

圖15 13分區(qū)下差分進化算法種群數(shù)量和終止代數(shù)對比圖

圖16 13分區(qū)差分進化算法最優(yōu)解收斂過程

3.3 算法結果分析

以13分區(qū)的不同算法收斂情況進行對比分析，得到算法對非線性調配的模型的求解效果，匯總四種算法的最優(yōu)參數(shù)組合，最優(yōu)目標函數(shù)和算法運行時間，得到表6所示的各種算法效果對比表。

表6 13分區(qū)非線性調配的各算法效果對比

由表6可以得出，對于13分區(qū)的情況，改進模擬退火得到最優(yōu)解； 遺傳算法與改進模擬退火的計算時間大致相同，得到的解比改進算法的最優(yōu)解大0.4%； 粒子群算法計算時間最長，得到成本解比改進算法大5.3%； 差分進化算法的計算時間最短，但是其成本相較于改進算法求得的最優(yōu)解大1.2%。綜上，改進模擬退火的解最優(yōu)且運行時間較短，選取改進模擬退火作為非線性土方調配的最優(yōu)求解方法。

在此基礎上，為了解不同問題規(guī)模對改進模擬退火的求解性能和參數(shù)組合的影響，利用改進模擬退火對該項目6分區(qū)、13分區(qū)和29分區(qū)情況下的非線性土方調配求解，得到不同分區(qū)數(shù)量下最優(yōu)參數(shù)組合、最優(yōu)目標函數(shù)和最優(yōu)組合的函數(shù)圖像，如表7所示。

表7 不同分區(qū)下改進模擬退火算法效果對比

由前文可知，分區(qū)數(shù)量為6時，有多個參數(shù)組合都能得到最優(yōu)解，此情況下的求解難度較??； 分區(qū)數(shù)量為13時，僅有一個最優(yōu)參數(shù)組合且得到的成本解是收斂的，比6分區(qū)得到的成本小，求解難度適中； 分區(qū)數(shù)量為29的情況下也僅有一個最優(yōu)參數(shù)組合，此參數(shù)組合下的成本仍有下降可能，求解難度最大。通過上述分析可知，隨著分區(qū)數(shù)量的增加，求解難度逐步增加，對于不同規(guī)模的問題，最優(yōu)參數(shù)組合也不同。

4 結論

本文將非線性規(guī)劃問題轉化為函數(shù)極值問題并利用多種啟發(fā)式算法對其求解，得到了以下結論：

(1)對于非線性規(guī)劃模型，可以將問題轉化為無約束的函數(shù)優(yōu)化問題并進行求解；在本案例中，改進模擬退火是求解非線性土方調配的最優(yōu)方法，其求得最優(yōu)成本為1 765 672元；

(2)隨著非線性土方調配問題中分區(qū)數(shù)量的增加，待求解問題的規(guī)模變大，改進模擬退火算法的求解難度會變大，所以需要采用不同的參數(shù)組合來求得最優(yōu)解。