余 佳，焦 錚，蘇 哲，任炳昱，王國浩，肖 堯

(天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072)

引水隧洞施工過程面臨著活動時間隨機(jī)變化、施工機(jī)械故障以及塌方、涌水等風(fēng)險因素。然而，現(xiàn)有的施工方案優(yōu)選研究缺乏考慮風(fēng)險因素的影響，這導(dǎo)致無法獲得合理的施工方案。因此，在考慮風(fēng)險因素影響下進(jìn)行引水隧洞施工方案優(yōu)選具有重要的意義。考慮到隧洞施工待選方案數(shù)量龐大，有必要首先采用多目標(biāo)優(yōu)化方法篩選出Pareto方案集，然后綜合考慮專家意見從中優(yōu)選出最終方案。

在進(jìn)行多目標(biāo)優(yōu)化問題[1]求解時，粒子群優(yōu)化(PSO)算法由于結(jié)構(gòu)簡單、收斂速度快、參數(shù)設(shè)置少等優(yōu)點，已被廣泛應(yīng)用[2-3]。針對傳統(tǒng)PSO算法易陷入早熟收斂等缺點，國內(nèi)外學(xué)者通過改進(jìn)算法慣性權(quán)重和加速因子[4]，與人工魚群算法[5]、模擬退火算法[6]、蟻群算法[7]等結(jié)合來優(yōu)化粒子的搜索精度和效率。本文提出Levy飛行自適應(yīng)混沌粒子群(LFACPSO)算法，利用Levy飛行的強(qiáng)跳躍特性來改進(jìn)粒子的更新方式，提高粒子跳出局部最優(yōu)的能力；此外，采用混沌算法來初始化粒子，并自適應(yīng)地調(diào)整粒子慣性權(quán)重系數(shù)以提高粒子的局部和全局搜索能力。

國內(nèi)外學(xué)者對水電地下工程和大壩工程施工方案優(yōu)選進(jìn)行了研究，如畢磊等[8]提出了基于數(shù)據(jù)包絡(luò)分析的隧洞施工方案評價與優(yōu)選方法，趙延喜等[9]利用層次分析法和模糊集法構(gòu)建了TBM施工風(fēng)險綜合評判計算模型，何偉[10]利用灰色系統(tǒng)方案決策理論進(jìn)行隧洞施工方案優(yōu)選，張孝遠(yuǎn)等[11-13]采用模糊多屬性決策法、逼近理想解法、數(shù)據(jù)驅(qū)動的D-AHP方法對大壩工程施工方案進(jìn)行了優(yōu)選。雖然這些研究考慮了模糊、灰色等不確定性，但是現(xiàn)有研究缺乏從活動時間隨機(jī)變化、施工機(jī)械故障以及塌方、涌水等風(fēng)險因素入手分析其對優(yōu)化過程的影響，且在方案優(yōu)選過程中缺乏考慮決策者的猶豫度。因此，本文提出考慮風(fēng)險因素影響的引水隧洞施工方案優(yōu)選方法，建立綜合考慮風(fēng)險因素影響的施工方案優(yōu)化模型，采用LFACPSO算法求解模型獲得Pareto方案集，并采用多屬性群決策直覺模糊熵權(quán)冪平均(IFEWPA)方法，在考慮專家猶豫度條件下從Pareto方案集中優(yōu)選出最佳方案。

1 考慮風(fēng)險因素影響的引水隧洞施工方案優(yōu)選方法

1.1 考慮風(fēng)險因素影響的引水隧洞施工方案優(yōu)化數(shù)學(xué)模型

建立考慮風(fēng)險因素影響的引水隧洞施工方案優(yōu)化數(shù)學(xué)模型。以考慮風(fēng)險因素影響的工期和成本均值最小為目標(biāo)建立目標(biāo)函數(shù)。由于引水隧洞工程可能涉及多個施工標(biāo)段，因此優(yōu)化目標(biāo)中常包含多個工期目標(biāo)，目標(biāo)函數(shù)g為

(1)

約束條件為

(2)

1.2 考慮風(fēng)險因素影響的引水隧洞施工仿真方法

引水隧洞施工中面臨著活動時間隨機(jī)波動、施工機(jī)械故障以及塌方、涌水等風(fēng)險的影響。采用網(wǎng)絡(luò)計劃技術(shù)(CPM)和循環(huán)網(wǎng)絡(luò)仿真技術(shù)(CYCLONE)耦合方法[14]來模擬各類風(fēng)險綜合影響下的隧洞施工過程。其中，活動時間隨機(jī)波動通過概率分布來表示[15]；施工機(jī)械故障風(fēng)險的影響通過為鉆機(jī)、自卸汽車、裝載機(jī)等機(jī)械賦予故障間隔時間及故障處理時間參數(shù)，并將其考慮到鉆孔和出渣活動模擬過程中；塌方、涌水等突發(fā)地質(zhì)風(fēng)險事件的影響通過采用貝葉斯網(wǎng)絡(luò)[16]計算風(fēng)險發(fā)生概率并將其嵌入CYCLONE模型來實現(xiàn)。采用貝葉斯網(wǎng)絡(luò)計算風(fēng)險事件發(fā)生概率的過程詳見文獻(xiàn)[17]。

1.3 基于LFACPSO算法的引水隧洞施工方案多目標(biāo)優(yōu)化方法

1.3.1 LFACPSO算法

1.3.1.1 基于Levy飛行的粒子更新方法

Levy飛行是一種非高斯隨機(jī)過程，是一種短距離搜索與偶爾長距離行走相間的行走方式，其步長S服從Levy分布。目前常采用Mantegna等[18]提出的算法模擬Levy分布：

(3)

(4)

在優(yōu)化方案求解過程中，當(dāng)粒子停留在某一局部位置的次數(shù)大于預(yù)設(shè)值時(本文取10次[20])，則在粒子速度更新公式中引入Levy飛行，以使其跳出局部最優(yōu)。基于Levy飛行的粒子速度和位置更新公式如下:

υi,t+1=ωL(xi,t)+c1ζ(pi,t-xi,t)+c2η(pt-xi,t)

(5)

xi,t+1=xi,t+υi,t+1

(6)

其中L(xi,t)=0.01Sxi,trandom(size(xi,t))

式中：xi,t——第i個粒子在第t代進(jìn)化過程中的位置；pi,t——第i個粒子搜索到的歷史最優(yōu)位置；pt——整個粒子群搜索到的歷史最優(yōu)位置；υi,t+1——第i個粒子在t+1代進(jìn)化過程中的速度；ω——慣性權(quán)重；c1、c2——學(xué)習(xí)因子，取c1=c2=2.0；ζ、η——符合[0,1]均勻分布的隨機(jī)數(shù)；size(xi,t)——xi,t的維度；random(size(xi,t))——在xi,t各維度上生成0～1隨機(jī)數(shù)。

1.3.1.2 粒子混沌初始化

混沌運動具有遍歷性、隨機(jī)性、對初值的敏感性等特點[21]。采用混沌序列來初始化粒子有助于提高種群多樣性及粒子搜索的遍歷性，克服隨機(jī)初始化的不足。采用如下的Logistic混沌方程：

Zq+1=dZq+(1-Zq)

(7)

式中：Zq——第q次迭代后的混沌變量取值；d——混沌系統(tǒng)的控制參數(shù)，取d=4。初始值Z0∈(0,1)且Z0?{0.25,0.50,0.75}時，系統(tǒng)處于完全的混沌狀態(tài)。

1.3.1.3 自適應(yīng)慣性權(quán)重

采用自適應(yīng)調(diào)整慣性權(quán)重的策略，使慣性權(quán)重隨迭代次數(shù)的改變而變化。慣性權(quán)重ω的調(diào)整方式如下：

(8)

式中：ωt——第t次迭代的慣性權(quán)重值；ωmax、ωmin——最大、最小慣性權(quán)重；tmax——最大迭代次數(shù)。

1.3.2 基于LFACPSO算法的施工方案多目標(biāo)優(yōu)化流程

基于LFACPSO算法的施工方案多目標(biāo)優(yōu)化流程如下：(a) 確定每個工序的待選方案，采用考慮風(fēng)險因素影響的仿真方法計算各方案下的目標(biāo)值；(b) 初始化算法參數(shù)，包括粒子數(shù)量、最大迭代次數(shù)、學(xué)習(xí)因子c1和c2、最大和最小慣性權(quán)重等；(c)采用混沌理論初始化粒子，其中一個粒子對應(yīng)一種施工方案組合；(d) 計算初始粒子的適應(yīng)度值，并將粒子當(dāng)前的適應(yīng)度值設(shè)為其初始最優(yōu)適應(yīng)度值；(e)選擇領(lǐng)導(dǎo)粒子，判斷粒子位置未更新的次數(shù)是否大于10次，若未大于10次，則采用傳統(tǒng)的粒子更新公式[22]來更新速度和位置，否則，采用式(5)(6)來更新粒子的速度和位置；(f)選取非支配解，并根據(jù)式(8)更新慣性權(quán)重；(g)判斷算法是否滿足終止條件，若滿足，則輸出Pareto方案集合，否則，返回步驟(e)繼續(xù)迭代。

1.4 施工方案多屬性群決策IFEWPA方法

施工方案多屬性群決策IFEWPA方法首先采用IFPWA算子，在考慮專家猶豫度的條件下進(jìn)行專家意見融合，獲得集成矩陣；然后采用直覺模糊交叉熵來計算屬性權(quán)重，以優(yōu)選出最終方案。

1.4.1 基于IFPWA算子的專家意見融合方法

將專家評價用直覺模糊數(shù)來表示，形成直覺模糊評價矩陣。第k個專家的直覺模糊矩陣Ak為

(9)

式中：αi,j,k——專家k對方案i中目標(biāo)j的直覺模糊評價，αi,j,k=(μi,j,k,vi,j,k,γi,j,k)；μi,j,k、vi,j,k、γi,j,k——專家評價的隸屬度、非隸屬度、猶豫度，μi,j,k,vi,j,k,γi,j,k∈[0,1]；m——待選方案數(shù)量；n——優(yōu)化目標(biāo)數(shù)量。

采用Song等[23]提出的方法計算直覺模糊集之間的相似度，并替代P-A算子中的支持度函數(shù)，得到IFPWA算子，該算子的支持度函數(shù)Sup計算公式為

(10)

(11)

利用IFPWA算子將p個直覺模糊評價矩陣Ak(k=1,2,…,p)融合為集成評價矩陣A：

(12)

式中：ξi,j,k——各專家評價直覺模糊集之間的相互支持程度，ξ越大則相互支持程度越高。

1.4.2 基于直覺模糊交叉熵的屬性權(quán)重計算方法

存在直覺模糊集A={〈x,μA(x),vA(x)〉|x∈X}和B={〈x,μB(x),vB(x)〉|x∈X}，其中μA(x)、vA(x)分別為有限非空集X中元素x對集合A的隸屬度和非隸屬度，μB(x)、vB(x)分別為元素x對集合B的隸屬度和非隸屬度。記a=1+μA(x)-vA(x)，b=1+μB(x)-vB(x)，c=1-μA(x)+vA(x)，d=1-μB(x)+vB(x)，則A與B的直覺模糊交叉熵為

(13)

由于上述直覺模糊交叉熵不滿足對稱性，故直覺模糊集的對稱交叉熵為

C*(A,B)=C(A,B)+C(B,A)

(14)

在多屬性決策過程中，屬性Gj的相似度F(Gj)可表示為

(15)

F(Gj)越高，表明屬性Gj越不重要，一般賦予較低的權(quán)重值。因此，對于m個方案n個屬性的多屬性決策問題，屬性Gj(j=1,2,…,n)的權(quán)重wj可由下式計算得到：

(16)

式中：ri,j——集成評價矩陣中第i個方案下第j個目標(biāo)的直覺模糊數(shù)。

確定所有屬性權(quán)重wj(j=1,2，…,n)后，將m個待選方案的n個指標(biāo)綜合起來，得到全局偏好值Ri(式(17))，根據(jù)該值將所有待選方案排序即可獲得最優(yōu)方案。

(17)

式中：μij、vij——集成矩陣A中第i個方案下第j個屬性的直覺模糊評價的隸屬度與非隸屬度。

方案i的得分函數(shù)值ti為

(18)

2 工 程 應(yīng) 用

某引水隧洞C2標(biāo)段采用鉆爆法施工，包括1號隧洞和2號隧洞的西端洞段，1號隧洞長4 700 m，2號隧洞長6 000 m。1號、2號隧洞均分成上部和下部開挖，且兩條隧洞之間每間隔500 m布置一條橫通道，整個標(biāo)段共設(shè)置8條橫通道。采用考慮風(fēng)險因素影響的施工方案優(yōu)選方法進(jìn)行C2標(biāo)段施工方案優(yōu)選。

2.1 基于LFACPSO算法的施工方案多目標(biāo)優(yōu)化

表1 C2標(biāo)段可選機(jī)械參數(shù)

采用考慮風(fēng)險因素影響的引水隧洞施工仿真方法獲得各工序在不同施工機(jī)械方案下的工期及成本，然后采用LFACPSO算法進(jìn)行施工方案多目標(biāo)求解。對200個粒子進(jìn)行混沌初始化，設(shè)置最大迭代次數(shù)為1 000次，算法學(xué)習(xí)因子c1=c2=2.0，最大慣性權(quán)重ωmax=1.0，最小慣性權(quán)重ωmin=0.2。通過算法求解獲得27個Pareto方案如圖1所示，圖中圓的半徑越大，表示成本越高。

圖1 Pareto方案集Fig.1 Pareto scheme set

2.2 基于IFEWPA方法的施工方案多屬性群決策

圖2 專家集成評價結(jié)果Fig.2 Integrated evaluation results of experts

根據(jù)獲得的各目標(biāo)權(quán)重及集成的專家評價結(jié)果，由式(17)計算得到各施工方案的全局偏好值，并根據(jù)式(18)獲得相應(yīng)的得分函數(shù)值。由圖3可知，各方案的得分函數(shù)值均不相同，且方案17的得分函數(shù)值最大，為綜合最優(yōu)方案。在方案17中，1號隧洞上部配置21臺手風(fēng)鉆、10臺自卸汽車、2臺裝載機(jī)，2號隧洞上部配置21臺手風(fēng)鉆、12臺自卸汽車、2臺裝載機(jī)，1號隧洞下部配置6臺手風(fēng)鉆、2臺自卸汽車、1臺裝載機(jī)，2號隧洞下部配置10臺手風(fēng)鉆、3臺自卸汽車、1臺裝載機(jī)，橫通道配置6臺手風(fēng)鉆、4臺自卸汽車、1臺裝載機(jī)。

圖3 施工方案的全局偏好值及得分函數(shù)值Fig.3 Global preference values and score function values of construction scheme

2.3 對比分析與討論

為了證明本文所提出方法的有效性，以2號隧洞上部800 m開挖段為對象，將優(yōu)選方案17下考慮風(fēng)險綜合影響的仿真進(jìn)度與不考慮風(fēng)險影響的仿真進(jìn)度、實際進(jìn)度、計劃進(jìn)度進(jìn)行對比，結(jié)果如圖4所示。由圖4可知，在800 m段開挖過程中，相比于不考慮風(fēng)險影響的仿真進(jìn)度和計劃進(jìn)度，考慮風(fēng)險綜合影響的仿真進(jìn)度與實際進(jìn)度更加貼近。在完成800 m開挖時，實際工期為175 d，考慮風(fēng)險綜合影響的平均仿真工期為169 d，不考慮風(fēng)險影響的平均仿真工期為155 d，計劃工期為159 d，考慮風(fēng)險綜合影響的仿真工期更接近實際工期。

圖4 2號隧洞上部施工進(jìn)度對比Fig.4 Comparison of construction schedule of number 2 tunnel upper section

為了證明提出的LFACPSO算法在進(jìn)行隧洞施工方案多目標(biāo)優(yōu)化方面的優(yōu)越性，選取不帶Levy飛行的自適應(yīng)混沌粒子群(ACPSO)算法、傳統(tǒng)的PSO算法和帶精英策略的非支配排序遺傳(NSGA-Ⅱ)算法分別進(jìn)行優(yōu)化模型的求解。各算法的初始種群均設(shè)置為200個，迭代次數(shù)為1 000次，每種算法分別重復(fù)計算10次，并對10次計算的結(jié)果進(jìn)行對比分析，結(jié)果如圖5所示。

圖5 方案優(yōu)化模型的不同算法求解結(jié)果對比Fig.5 Comparison of solution results of different algorithms for scheme optimization model

從算法的優(yōu)化結(jié)果來看，LFACPSO算法所得的平均Pareto優(yōu)化方案數(shù)量最大，為24.6個，而ACPSO算法、PSO算法、NSGA-Ⅱ算法所得的平均方案數(shù)量分別為20.6個、21個、20.7個。由此可知，相比其他算法，LFACPSO算法求解方案多目標(biāo)優(yōu)化模型的能力更強(qiáng)，所得Pareto方案集更加完備。

從算法的魯棒性來看，LFACPSO、ACPSO、PSO、NSGA-Ⅱ算法10次計算所得Pareto方案數(shù)量的標(biāo)準(zhǔn)差分別為2.61、5.12、5.10、3.95。對比其他算法，LFACPSO算法求解結(jié)果更穩(wěn)定，魯棒性更高。

3 結(jié) 語

本文提出了考慮風(fēng)險因素影響的引水隧洞施工方案優(yōu)選方法。該方法首先在綜合考慮活動時間隨機(jī)變化、施工機(jī)械故障以及塌方、涌水等風(fēng)險因素影響下建立施工方案優(yōu)化模型，克服了現(xiàn)有施工方案優(yōu)選研究缺乏考慮風(fēng)險因素影響的不足。其次，針對傳統(tǒng)PSO算法局部搜索能力差、易陷入早熟收斂的不足，提出了LFACPSO算法來求解優(yōu)化模型，該算法采用基于Levy飛行的粒子更新方式來提高粒子跳出局部最優(yōu)的能力，同時采用混沌算法來初始化粒子，并自適應(yīng)地調(diào)整粒子慣性權(quán)重系數(shù)以提高算法的全局與局部搜索能力。最后，針對現(xiàn)有施工方案優(yōu)選研究未能體現(xiàn)專家猶豫度的不足，提出了施工方案多屬性群決策的IFEWPA方法，在考慮專家猶豫度的條件下優(yōu)選出最終方案。相比計劃進(jìn)度和不考慮風(fēng)險影響的仿真進(jìn)度，綜合考慮風(fēng)險影響的仿真進(jìn)度與實際進(jìn)度更加吻合，證明了本文提出方法的有效性。將LFACPSO算法與ACPSO、PSO和NSGA-Ⅱ算法進(jìn)行對比，結(jié)果表明，LFACPSO算法所得的Pareto方案集更加完備，結(jié)果更穩(wěn)定，證明了該算法的有效性與魯棒性。