基于子集模擬的建設(shè)項(xiàng)目離散型工期 - 成本優(yōu)化算法研究

王 家， 韓淙吉， 陳雅含， 李炎彪

(湖南大學(xué) 土木工程學(xué)院， 湖南 長(zhǎng)沙 410082)

工期和成本是建設(shè)項(xiàng)目管理的重要目標(biāo)，且兩者相互聯(lián)系、相互制約。項(xiàng)目工期的壓縮，往往伴隨著勞動(dòng)力的增加、更高效率施工設(shè)備和施工方法的采用等，進(jìn)而導(dǎo)致項(xiàng)目成本的增加[1]。因此，在進(jìn)度計(jì)劃管理決策時(shí)，需綜合考慮工期與成本間的合理權(quán)衡，考察項(xiàng)目的工期 - 成本優(yōu)化問題。同時(shí)，建設(shè)項(xiàng)目的資源(如施工設(shè)備、勞動(dòng)力)分配通常具有離散型的特征，項(xiàng)目中各工序的實(shí)施方式只能從有限個(gè)執(zhí)行模式中選取[2～5]。因此，離散型工期 - 成本優(yōu)化問題更為貼合工程實(shí)際，受到項(xiàng)目管理者和研究人員的廣泛關(guān)注。

離散型工期 - 成本優(yōu)化問題屬于組合優(yōu)化問題，其求解方法主要分為精確算法和啟發(fā)式算法。其中，精確算法以整數(shù)規(guī)劃為代表，旨在尋求問題的準(zhǔn)確最優(yōu)解[6]。但離散型工期 - 成本優(yōu)化問題為NP困難問題，其復(fù)雜程度隨項(xiàng)目中工序數(shù)量的增加而急劇上升[7]。當(dāng)優(yōu)化問題中項(xiàng)目的工序數(shù)量較多時(shí)，精確算法并不適用。此時(shí)，適宜采用啟發(fā)式算法，如遺傳算法[8～11]、粒子群算法[12]、蟻群算法[13,14]等，來尋找問題的最優(yōu)解或近似最優(yōu)解。但是，啟發(fā)式算法的求解具有隨機(jī)性，每次運(yùn)行獲得的最優(yōu)解不一定相同(不穩(wěn)定)，現(xiàn)有的啟發(fā)式算法在最優(yōu)解獲取穩(wěn)定性上仍有較大的改進(jìn)空間。

為此，本文針對(duì)建設(shè)項(xiàng)目的離散型工期 - 成本優(yōu)化問題，提出一種基于子集模擬法的新型啟發(fā)式算法。同時(shí)，為解決算法中可行域及收縮區(qū)域中隨機(jī)樣本點(diǎn)的高效產(chǎn)生難題，本文建議采用“帶有反射壁的隨機(jī)游動(dòng)”的馬爾科夫鏈蒙特卡羅方法。通過算例驗(yàn)證，與應(yīng)用較廣的遺傳算法相比，本文建議的優(yōu)化算法性能更好，在最優(yōu)解的獲取穩(wěn)定性上有較大改進(jìn)。

1 建設(shè)項(xiàng)目離散型工期 - 成本優(yōu)化問題模型

建設(shè)項(xiàng)目離散型工期 - 成本優(yōu)化問題中，涉及的項(xiàng)目工序執(zhí)行模式為有限個(gè)。因此，考慮一包含N項(xiàng)工序的建設(shè)項(xiàng)目，并設(shè)任一工序i的執(zhí)行模式總數(shù)為Ki，i=1,…,N。項(xiàng)目工序的不同執(zhí)行模式，可能對(duì)應(yīng)于不同的施工工藝，也可能對(duì)應(yīng)于不同的的人力、機(jī)械投入數(shù)量，進(jìn)而導(dǎo)致不同的工序成本和工序工期。例如，表1為某橋梁工程項(xiàng)目中大孔板梁安裝工序的3種備選執(zhí)行模式，對(duì)應(yīng)于3種不同的施工工藝，其對(duì)應(yīng)的成本和工期亦各不相同。

表1 大孔板梁安裝工序的備選執(zhí)行模式

假設(shè)工序i在第j種執(zhí)行模式下的持續(xù)時(shí)間和直接成本分別為ti,j和ci,j。優(yōu)化模型中，一般采用xi,j作為執(zhí)行模式的指示變量，即當(dāng)工序i選擇執(zhí)行模式j(luò)時(shí)，xi,j=1，否則xi,j=0。因任一工序同時(shí)只能執(zhí)行一種模式，所以指示變量xi,j需滿足：

(1)

項(xiàng)目的實(shí)施方案，由項(xiàng)目中所有工序的執(zhí)行模式所確定，可表示為包含所有指示變量的矩陣X=[xij,i=1,…,N,j=1,…,Ki]。例如，表2為某建設(shè)項(xiàng)目的實(shí)施方案，該項(xiàng)目包含6項(xiàng)工序，每項(xiàng)工序的執(zhí)行模式總數(shù)分別為5,4,3,4,4,5。由指示變量xi,j的含義可知，該項(xiàng)目實(shí)施方案中，工序1～6分別選擇第3、第2、第2、第4、第1及第4種執(zhí)行模式。

表2 某建設(shè)項(xiàng)目實(shí)施方案矩陣X

對(duì)應(yīng)項(xiàng)目的某一實(shí)施方案X，工序i的持續(xù)時(shí)間ti和直接成本ci可表示為：

(2)

(3)

項(xiàng)目的總工期T(X)，可采用關(guān)鍵路徑法，由關(guān)鍵路徑上所有工序的持續(xù)時(shí)間求和得到[15]：

(4)

式中:A為關(guān)鍵路徑上的工序集合?？紤]項(xiàng)目的要求完工工期TD時(shí)，項(xiàng)目的實(shí)施方案X需滿足約束條件：

T(X)≤TD

(5)

項(xiàng)目的總成本C(X)包含直接總成本Cd(X)和間接總成本Ci(X)。其中，直接總成本Cd(X)對(duì)應(yīng)于所有工序的直接成本之和。間接總成本Ci(X)一般假定為項(xiàng)目總工期的線性函數(shù)，通過引入單位工期的間接成本l來計(jì)算[16]。因此，對(duì)應(yīng)于某一項(xiàng)目實(shí)施方案，項(xiàng)目的總成本C(X)可表示為：

C(X)=Cd(X)+Ci(X)

=Cd(X)+l×T(X)

(6)

在此基礎(chǔ)上，建設(shè)項(xiàng)目的離散型工期 - 成本優(yōu)化模型可描述為如下的0-1型整數(shù)規(guī)劃問題[17]：

(7)

建設(shè)項(xiàng)目的實(shí)施方案，也可采用向量的表示形式，寫為Z=[Z1,Z2,…,ZN]。此時(shí)，向量Z中的元素Zi，i=1,…,N代表工序i的執(zhí)行模式編號(hào)，為[1,Ki]內(nèi)的自然數(shù)。由項(xiàng)目實(shí)施方案的矩陣X與向量Z的含義可知，X與Z間存在一一對(duì)應(yīng)關(guān)系。例如，表2的矩陣X對(duì)應(yīng)的實(shí)施方案中，工序1～6分別選擇第3、第2、第2、第4、第1及第4種執(zhí)行模式，因此該方案的向量表示為Z=[3,2,2,4,1,4]。由向量Z的含義可知，其對(duì)應(yīng)的矩陣X自動(dòng)滿足公式(1)對(duì)應(yīng)約束條件。因此，在項(xiàng)目實(shí)施方案的向量表示下，建設(shè)項(xiàng)目的離散型工期 - 成本優(yōu)化模型可寫為：

(8)

其中，為得到實(shí)施方案Z下的項(xiàng)目總工期T(Z)和總成本C(Z)，可首先求得Z對(duì)應(yīng)的實(shí)施方案矩陣X(Z)，進(jìn)而利用公式(4)和公式(6)計(jì)算。

2 基于子集模擬的離散型工期 - 成本優(yōu)化算法

子集模擬法(Subset Simulation)是由Au和Beck提出的、針對(duì)高維空間中可靠度估算分析的高效數(shù)值模擬算法[18]。該方法通過引入一組中間失效事件，將偶發(fā)事件的小概率表示為一組較大的條件概率的乘積，進(jìn)而將小概率估算問題轉(zhuǎn)化為較大條件概率的估計(jì)問題，以提高數(shù)值模擬算法的計(jì)算效率。通過建立可靠度問題與優(yōu)化問題間的內(nèi)在聯(lián)系，子集模擬法在結(jié)構(gòu)優(yōu)化問題中取得了較好應(yīng)用[19]。但是，在建設(shè)項(xiàng)目工期 - 成本優(yōu)化問題中，尚未見到基于子集模擬法的相關(guān)應(yīng)用。

2.1 基于子集模擬的優(yōu)化框架

針對(duì)公式(8)的工期 - 成本優(yōu)化問題，其可行域?yàn)棣?={Z:T(Z)≤TD,1≤Zi≤Ki,i=1,…,N}?；谧蛹M的優(yōu)化算法中，通過引入一系列與目標(biāo)函數(shù)相關(guān)的、遞減的邊界值c1>c2>…>cJ，在可行域Ω0內(nèi)定義了一組逐漸縮小的區(qū)域Ωj={Z:Z∈Ω0,C(Z)≤cj}，j=1,…,J。其中，定義的區(qū)域間存在子集關(guān)系，即Ω0?Ω1?Ω2?…?ΩJ?；谧蛹M的優(yōu)化算法中，通過逐次在可行域Ω0及逐漸縮窄的區(qū)域Ωj(j=1,…,J)中按均勻分布隨機(jī)取樣，將搜索范圍逐漸縮窄到最優(yōu)解附近的較小區(qū)域，以得到問題的最優(yōu)解。

基于子集模擬的優(yōu)化算法中有兩個(gè)重要參數(shù)：指定區(qū)域中的抽樣樣本數(shù)M和條件概率參數(shù)p0。參數(shù)M和p0的取值，需保證Mp0和1/p0均為正整數(shù)。條件概率參數(shù)p0決定了迭代中考察區(qū)域的收縮快慢，如p0過大，迭代中考察區(qū)域的收縮變慢，優(yōu)化算法需要較多的迭代次數(shù)。同時(shí)，Mp0對(duì)應(yīng)于每一個(gè)迭代中收縮區(qū)域的初始樣本點(diǎn)數(shù)量，如p0過小，則需要較大的樣本數(shù)量M才能保證Mp0的數(shù)值不至過小，以達(dá)到對(duì)收縮區(qū)域的有效搜索。綜合目前研究成果，p0一般在0.05～0.3區(qū)間取值，進(jìn)而可根據(jù)計(jì)算資源的承受能力確定抽樣樣本數(shù)量M。

圖1～4描述了參數(shù)M=25和p0=0.2時(shí)的上述流程。針對(duì)圖1中隨機(jī)產(chǎn)生(滿足可行域Ω0內(nèi)的均勻分布)的25個(gè)抽樣樣本點(diǎn)，圖2將它們的目標(biāo)函數(shù)值按升序排列。此時(shí)，定義下一收縮區(qū)域Ω1的臨界值c1，可由這25個(gè)目標(biāo)函數(shù)值中第Mp0=5小的數(shù)值確定，見圖2。對(duì)應(yīng)上述確定的c1，可行域Ω0內(nèi)隨機(jī)產(chǎn)生的M=25個(gè)抽樣樣本點(diǎn)，有Mp0=5個(gè)落在收縮區(qū)域Ω1={Z:Z∈Ω0,C(Z)≤c1}內(nèi)，見圖3。為產(chǎn)生收縮區(qū)域Ω1內(nèi)的M=25個(gè)抽樣樣本點(diǎn)，以目前的Mp0=5個(gè)樣本點(diǎn)為種子，以Ω1內(nèi)的均勻分布為目標(biāo)概率分布，采用MCMCS方法產(chǎn)生5條馬爾科夫鏈鏈條，見圖4。每條鏈條以其中一個(gè)種子為起點(diǎn)，并額外產(chǎn)生(1/p0-1)=4個(gè)狀態(tài)點(diǎn)，以確?？倶颖緮?shù)為Mp0(1/p0-1+1)=M=25。

圖1 可行域Ω0內(nèi)隨機(jī)產(chǎn)生的抽樣樣本點(diǎn)

圖2 定義收縮區(qū)域Ω1的邊界值c1

圖3 收縮區(qū)域Ω1及落入其中的初始樣本點(diǎn)

圖4 收縮區(qū)域Ω1的抽樣樣本點(diǎn)

2.2 收縮區(qū)域Ωj={Z:Z∈Ω0,C(Z)≤cj}內(nèi)隨機(jī)樣本點(diǎn)的產(chǎn)生

如2.1節(jié)所述，基于子集模擬的優(yōu)化算法中，需根據(jù)收縮區(qū)域Ωj內(nèi)的少量樣本點(diǎn)，利用MCMCS方法產(chǎn)生要求數(shù)量的更多樣本點(diǎn)。本文采用“帶有反射壁的隨機(jī)游動(dòng)”方法[21]來實(shí)現(xiàn)馬爾科夫鏈的產(chǎn)生。具體而言，為產(chǎn)生收縮區(qū)域Ωj內(nèi)Z1為起點(diǎn)的馬爾科夫鏈{Z1,Z2,…}，采用如下Zk=[Zk,1,Zk,2,…,Zk,N]到Zk+1=[Zk+1,1,Zk+1,2,…,Zk+1,N]的迭代操作：

(1)從樣本點(diǎn)Zk=[Zk,1,Zk,2,…,Zk,N]中，隨機(jī)選擇一可變?cè)豘k,r(對(duì)應(yīng)工序r的執(zhí)行模式總數(shù)Kr>1)。

2.3 可行域Ω0內(nèi)隨機(jī)樣本點(diǎn)的產(chǎn)生

2.4 基于子集模擬的離散型工期 - 成本優(yōu)化算法總結(jié)

針對(duì)建設(shè)項(xiàng)目的離散型工期 - 成本優(yōu)化問題，本文基于子集模擬法，建議的優(yōu)化算法的步驟總結(jié)如下：

(1)確定優(yōu)化算法中采用的參數(shù)，即在指定區(qū)域中隨機(jī)抽樣的樣本數(shù)量M、條件概率參數(shù)p0。

(4)算法采用的迭代終止原則為：消耗的計(jì)算資源達(dá)到設(shè)定的上限值。如不滿足迭代終止條件，令j=j+1，并重復(fù)步驟3。

3 案例分析

為檢驗(yàn)基于子集模擬的優(yōu)化算法的性能，本文考慮文獻(xiàn)[22]中的兩個(gè)離散型工期 - 成本優(yōu)化問題。兩個(gè)問題分別包含18項(xiàng)工序和63項(xiàng)工序，為項(xiàng)目離散型工期-成本優(yōu)化問題的經(jīng)典算例。同時(shí)，為增加算例的復(fù)雜程度，本文將18項(xiàng)工序的算例網(wǎng)絡(luò)首尾相接、重復(fù)3次，以構(gòu)成包含54項(xiàng)工序的新算例?？紤]到版面限制，項(xiàng)目中各工序在不同執(zhí)行模式下的工期和成本信息在此不再列出，請(qǐng)參考文獻(xiàn)[22]中的表1及表3。

針對(duì)算例問題的算法實(shí)現(xiàn)，筆者采用MatLab2020a自主編程完成，使用的計(jì)算機(jī)硬件配置為處理器Intel Core i7-7700、內(nèi)存16 G。

3.1 算例一(包含54項(xiàng)工序，真實(shí)最優(yōu)解已知)

圖5 算例一項(xiàng)目的單代號(hào)網(wǎng)絡(luò)圖

圖6 建議優(yōu)化算法中最優(yōu)目標(biāo)函數(shù)值隨迭代階段的變化

為檢驗(yàn)基于子集模擬的建議優(yōu)化算法的求解穩(wěn)定性，表3給出了建議優(yōu)化算法100次獨(dú)立運(yùn)行求解后的統(tǒng)計(jì)結(jié)果。由第2節(jié)優(yōu)化算法的闡述可知，算法中的計(jì)算資源主要用于給定實(shí)施方案Z下的項(xiàng)目總工期T(Z)和總成本C(Z)的分析。因此，迭代終止條件中的計(jì)算資源限制取60000次的項(xiàng)目工期和成本分析。目前，針對(duì)工期 - 成本優(yōu)化問題的啟發(fā)式算法間的對(duì)比研究較少，學(xué)界對(duì)各種啟發(fā)式算法的優(yōu)劣未達(dá)成共識(shí)。同時(shí)，遺傳算法因其自行概率搜索、運(yùn)算并行性、應(yīng)用不依賴問題種類的強(qiáng)魯棒性等特點(diǎn)，在工期 - 成本優(yōu)化問題中的應(yīng)用更為廣泛[8～11]。因此， 本文選擇遺傳算法與建議算法進(jìn)行對(duì)比分析?？紤]到種群規(guī)模m、交叉概率參數(shù)pc和變異概率參數(shù)pm對(duì)遺傳算法性能的影響，本文依據(jù)這些參數(shù)的一般取值范圍，進(jìn)行了大量不同參數(shù)組合下的遺傳算法性能檢驗(yàn)。檢驗(yàn)發(fā)現(xiàn)，針對(duì)本算例的離散型工期 - 成本優(yōu)化問題，遺傳算法在參數(shù)取值m=100，pc=0.1和pm=0.01時(shí)的求解性能最優(yōu)。表3提供了上述參數(shù)取值下遺傳算法100次獨(dú)立運(yùn)行求解后的統(tǒng)計(jì)結(jié)果。其中，遺傳算法求解采用相同的計(jì)算資源限制，即60000次的項(xiàng)目工期和成本分析。

表3 基于子集模擬的建議優(yōu)化算法與遺傳算法獨(dú)立運(yùn)行100次后的統(tǒng)計(jì)結(jié)果

由前人的研究成果可知，該問題的真實(shí)最優(yōu)解對(duì)應(yīng)目標(biāo)函數(shù)值為813810[22]。由表3可知，基于子集模擬的建議優(yōu)化算法100次運(yùn)行求解中，均找到了問題的真實(shí)最優(yōu)解(對(duì)應(yīng)目標(biāo)函數(shù)值813810)。相比之下，遺傳算法的100次運(yùn)行求解中，僅有78次獲得了問題的真實(shí)最優(yōu)解。由此可見，基于子集模擬的建議優(yōu)化算法獲取最優(yōu)解的穩(wěn)定性更高。

3.2 算例二(包含63項(xiàng)工序，真實(shí)最優(yōu)解未知)

圖7 算例二項(xiàng)目的單代號(hào)網(wǎng)絡(luò)圖

表4提供了上述參數(shù)取值下建議優(yōu)化算法與遺傳算法100次獨(dú)立運(yùn)行求解后的統(tǒng)計(jì)結(jié)果。對(duì)比可知，建議優(yōu)化算法獲得最優(yōu)目標(biāo)函數(shù)值的平均值為5468054，小于遺傳算法的平均值5487892。同時(shí)，建議優(yōu)化算法獲得最優(yōu)目標(biāo)函數(shù)值的最大值(5493430)、最小值(5453770)均優(yōu)于遺傳算法獲得最優(yōu)目標(biāo)函數(shù)值的最大值(5537150)、最小值(5457580)，且最優(yōu)目標(biāo)函數(shù)值的標(biāo)準(zhǔn)差更小。圖8給出了兩種算法100次獨(dú)立運(yùn)行獲得的最優(yōu)目標(biāo)函數(shù)值的分布情況。由圖8可見，建議優(yōu)化算法獲得的最優(yōu)目標(biāo)函數(shù)值更小，且分布更為集中，有91%(高于遺傳算法的47%)的最優(yōu)目標(biāo)函數(shù)值分布在區(qū)間[5450000，5483000)中，表明基于子集模擬的建議優(yōu)化算法獲取最優(yōu)解的穩(wěn)定性更高。

表4 基于子集模擬的建議優(yōu)化算法與遺傳算法獨(dú)立運(yùn)行100次獲得的最優(yōu)目標(biāo)函數(shù)值的統(tǒng)計(jì)結(jié)果

圖8 兩種優(yōu)化算法獨(dú)立運(yùn)行100次獲得的最優(yōu)目標(biāo)函數(shù)值的分布

4 結(jié) 論

針對(duì)進(jìn)度計(jì)劃管理決策中常見的離散型工期 - 成本優(yōu)化問題，本文提出了基于子集模擬法的優(yōu)化算法。同時(shí)，本文采用“帶有反射壁的隨機(jī)游動(dòng)”的馬爾科夫鏈蒙特卡羅方法，重點(diǎn)解決了優(yōu)化算法中可行域及收縮區(qū)域中隨機(jī)樣本點(diǎn)的產(chǎn)生問題。此外，本文利用包含54項(xiàng)工序和63項(xiàng)工序的兩個(gè)算例，檢驗(yàn)了建議優(yōu)化算法的性能，并與應(yīng)用較廣的遺傳算法進(jìn)行了對(duì)比。在包含54項(xiàng)工序，真實(shí)最優(yōu)解已知的算例中，建議優(yōu)化算法獲取真實(shí)最優(yōu)解的頻率為100%，高于遺傳算法的78%。而在包含63項(xiàng)工序，真實(shí)最優(yōu)解未知的算例中，建議優(yōu)化算法獲取的最優(yōu)目標(biāo)函數(shù)值更優(yōu)，且更為集中。研究表明，相較于遺傳算法，基于子集模擬的建議優(yōu)化算法在最優(yōu)解的獲取穩(wěn)定性上有較大改進(jìn)。