王永濤，劉 堅(jiān)，李家春，蔡家斌，索鑫宇，李繼學(xué)

（1.湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙410082；2.貴州省水利科學(xué)研究院，貴陽(yáng)550002；3.貴州大學(xué)機(jī)械工程學(xué)院，貴陽(yáng)550025）

0 引 言

李道西博士等[1]指出經(jīng)濟(jì)、合理、實(shí)用是設(shè)計(jì)灌溉管網(wǎng)的最終目的，管網(wǎng)優(yōu)化方案制定應(yīng)綜合考慮降低能耗、節(jié)省投資和提高灌水質(zhì)量。灌溉管網(wǎng)優(yōu)化布置可視為網(wǎng)絡(luò)圖優(yōu)化問題，最佳管網(wǎng)設(shè)計(jì)具有理論和實(shí)際雙重意義。傳統(tǒng)的灌溉管網(wǎng)優(yōu)化主要為微分法、非線性規(guī)劃法、線性規(guī)劃法、動(dòng)態(tài)規(guī)劃法與枚舉法等。當(dāng)求解問題規(guī)模較大時(shí)，傳統(tǒng)方法存在求解效率低等局限性。近年來，智能優(yōu)化算法在解決規(guī)模限制問題方面，具有明顯的優(yōu)勢(shì)。李援農(nóng)教授等[2]從控制灌溉面積最大、簡(jiǎn)化灌溉管網(wǎng)系統(tǒng)、減少管網(wǎng)建設(shè)及運(yùn)行管理成本等方面入手，基于最大控制面積和單位面積管網(wǎng)投資最低分別對(duì)毛管單向布置和雙向布置的微灌小區(qū)建立了目標(biāo)函數(shù)，設(shè)定灌水均勻度與壓力為約束條件，采用遺傳（Genetic Algorithm，GA）算法進(jìn)行求解，得到最優(yōu)化的灌溉管網(wǎng)布置與管徑組合方式。該研究從灌區(qū)整體出發(fā)，通過壓差（坡度或地形高差）與水頭損失（局部水頭損失、沿程水頭損失）計(jì)算得出各灌水節(jié)點(diǎn)的相對(duì)壓力值，具有方法簡(jiǎn)單且易于獲得最優(yōu)解的特點(diǎn)。馬孝義等[3]對(duì)樹狀輸配水管網(wǎng)進(jìn)行兩級(jí)優(yōu)化，第一級(jí)采用單點(diǎn)供水的原則，構(gòu)建GA 算法優(yōu)化布置模型；第二級(jí)根據(jù)整數(shù)編碼的原則，構(gòu)建投資最低的GA 算法優(yōu)化模型。胡杰華等[4]為解決灌溉管網(wǎng)布置和管徑優(yōu)化問題，建立管網(wǎng)投資最低的目標(biāo)函數(shù)，采用結(jié)合設(shè)計(jì)者實(shí)際工程經(jīng)驗(yàn)的Kruskal 算法求解灌溉管網(wǎng)最小生成樹。實(shí)踐表明，該方法能較好解決局部最優(yōu)解的問題，高效求解出符合工程實(shí)際的樹狀管網(wǎng)布置和管徑。范興業(yè)[5]、馬雪琴[6-7]等采用GA 算法開展樹狀灌溉管網(wǎng)的優(yōu)化設(shè)計(jì)，具有較強(qiáng)的指導(dǎo)意義，但受遺傳算子的影響，求解結(jié)果穩(wěn)定性較差。朱成立等[8]采用蟻群優(yōu)化（Ant Clony Optimization，ACO）算法對(duì)樹狀灌溉管網(wǎng)的管徑和布置進(jìn)行了優(yōu)化設(shè)計(jì)，實(shí)例研究表明，ACO 算法解決了搜索方向不明確、效率低等問題，對(duì)實(shí)際工作具有借鑒意義，但易于陷入局部最優(yōu)的問題。付玉娟等[9,10]認(rèn)為灌溉制度、灌水方式會(huì)直接影響灌溉配水管網(wǎng)的流量分配模式，進(jìn)而影響工程的總投資。同時(shí)通過對(duì)輪灌工作制度下灌溉管網(wǎng)優(yōu)化設(shè)計(jì)研究，提出了將輪灌組的劃分也進(jìn)行優(yōu)化，用基于整數(shù)編碼的遺傳算法進(jìn)行優(yōu)化，通過VB 和MATLAB 混合編程來實(shí)現(xiàn)。計(jì)算結(jié)果表明提出的方法效果明顯，可以應(yīng)用于生產(chǎn)。

然而，以上單一算法在解決收斂速度過快引起早熟現(xiàn)象或無法找到符合條件的解以及算法后期的局部搜索能力較差等問題局限性較大。牛寅[11]提出了模擬退火遺傳（Simulated Annealing Genetic Algorithm，SAGA）混合算法，針對(duì)江蘇省南通市港閘區(qū)花魚地村灌溉工程項(xiàng)目，分別對(duì)傳統(tǒng)經(jīng)驗(yàn)設(shè)計(jì)方法和SAGA 混合算法進(jìn)行了對(duì)比分析，結(jié)果表明優(yōu)化后的灌溉管網(wǎng)投資降低了19.77%，但其忽略了坡度對(duì)優(yōu)化結(jié)果的影響，在山丘區(qū)地勢(shì)起伏，坡度影響較大場(chǎng)合下，還需開展大量的研究。布谷鳥優(yōu)化算法模型簡(jiǎn)單、參數(shù)少及前期搜索較快的優(yōu)點(diǎn)，但存在搜索效率低、局部搜索效率不高等問題。同時(shí)，模擬退火算法優(yōu)點(diǎn)是局部搜索能力較強(qiáng)，缺點(diǎn)是具有較低的搜索效率。為此，馬燦等[12,13]將布谷鳥（Cuckoo Search，CS）算法與模擬退火（Simulated Annealing，SA）算法結(jié)合構(gòu)建（Simulated Annealing Cuckoo Search，SACS）混合算法，充分利用布谷鳥優(yōu)化算法全局搜索能力把握搜索方向，模擬退火算法的局部搜索能力提高搜索效率，以期得到更好的求解質(zhì)量。通過對(duì)經(jīng)典測(cè)試函數(shù)和旅行商問題進(jìn)行測(cè)試，結(jié)果表明，改進(jìn)后的混合模擬退火的布谷鳥算法在收斂速度以及尋優(yōu)精度方面有明顯提升，對(duì)于函數(shù)優(yōu)化問題和組合優(yōu)化問題都具有一定的優(yōu)勢(shì)。但其僅對(duì)經(jīng)典測(cè)試函數(shù)和旅行商問題進(jìn)行了測(cè)試，將算法應(yīng)用到工程優(yōu)化領(lǐng)域還需要大量研究。本文的創(chuàng)新點(diǎn)是在馬燦等人研究得基礎(chǔ)上，進(jìn)一步研究SACS 混合算法在樹狀灌溉管網(wǎng)優(yōu)化設(shè)計(jì)等方面得應(yīng)用，將SACS 混合算法應(yīng)用于具體灌溉工程，為解決灌溉工程管網(wǎng)優(yōu)化設(shè)計(jì)實(shí)際問題提供了新的思路。

1 樹狀灌溉管網(wǎng)優(yōu)化模型設(shè)計(jì)

對(duì)石阡縣龍?zhí)链笸筒枞~園區(qū)216 hm2現(xiàn)代水利試點(diǎn)區(qū)工程實(shí)地調(diào)查的基礎(chǔ)上，在等高線地形圖上明確水源及各灌溉節(jié)點(diǎn)的位置，采用緊鄰規(guī)劃的原則并結(jié)合實(shí)際情況初步確定灌溉管網(wǎng)連接圖，并確定水源至各需水節(jié)點(diǎn)的長(zhǎng)度、高程及供水流量等工程特征信息。

1.1 樹狀灌溉管網(wǎng)數(shù)學(xué)建模

對(duì)灌溉管網(wǎng)系統(tǒng)的優(yōu)化設(shè)計(jì)問題是指在管網(wǎng)布置一定，滿足灌水均勻度、灌溉施肥制度、輪灌制度等一系列設(shè)計(jì)準(zhǔn)則的前提下，實(shí)現(xiàn)建設(shè)投資、管理成本、運(yùn)行維護(hù)成本等綜合成本最低的管網(wǎng)最優(yōu)配置。灌溉管網(wǎng)系統(tǒng)涉及的參數(shù)較多，且各參數(shù)間包含復(fù)雜的非線性關(guān)系和耦合關(guān)系，因此該設(shè)計(jì)優(yōu)化問題可以看做是一個(gè)以綜合成本最小為優(yōu)化目標(biāo)，同時(shí)滿足多個(gè)約束條件的復(fù)雜的非線性問題的優(yōu)化。山區(qū)地形地勢(shì)起伏，多采用低壓灌溉管網(wǎng)，構(gòu)建非線性規(guī)劃模型，并以管網(wǎng)的一次性建設(shè)投資最小為優(yōu)化目標(biāo)[14]。為簡(jiǎn)化模型求解過程，本文采用懲罰函數(shù)將多約束優(yōu)化問題轉(zhuǎn)化為無約束優(yōu)化問題，進(jìn)一步提高模型的求解效率。

目標(biāo)函數(shù)：

式中：F′為管網(wǎng)一次性投入，元；Di為向第i節(jié)點(diǎn)供水管道的管徑，mm；Li為向第i節(jié)點(diǎn)供水的管道長(zhǎng)度，m；N為管網(wǎng)需水節(jié)點(diǎn)數(shù)；a、B、b為管道造價(jià)系數(shù)和指數(shù)。

（1）壓力約束：

式中：I(k)為從水源到管網(wǎng)第k個(gè)節(jié)點(diǎn)處所經(jīng)過的管段個(gè)數(shù)；M為標(biāo)準(zhǔn)管徑數(shù)；α為局部水頭損失的放大系數(shù)，通常取1.1；f、m、n為管道水頭損失中與管材相關(guān)的參數(shù)，見表1；E0為自由水面高程，m；Qij為從水源到第i需水節(jié)點(diǎn)供水管網(wǎng)中的第j段管道流量，m3/h；Dij為從水源到第i需水節(jié)點(diǎn)供水管網(wǎng)中的第j段管道直徑，mm；Lij為從水源到第i需水節(jié)點(diǎn)供水管網(wǎng)中的第j段管道管長(zhǎng)，mm；Ek為管網(wǎng)第k個(gè)節(jié)點(diǎn)的地面高程，m；Hkmin為滿足第k個(gè)節(jié)點(diǎn)處流量所需要的最小水頭，m。

表1 塑料管道常用的f、m、n值

（2）管徑約束：

①管徑應(yīng)滿足標(biāo)準(zhǔn)管徑：

式中：DW為可選標(biāo)準(zhǔn)管徑的集合。

②管徑沿水流方向應(yīng)滿足非增條件：

式中：Dij-1為從水源到第i需水節(jié)點(diǎn)供水路徑中的第j- 1段管道直徑，mm。

（3）流速約束：

式中：vj(k)為第j段管道流速，m/s；vmax為系統(tǒng)允許通過的最大流速，m/s；vmin為系統(tǒng)允許通過的最小流速，m/s。

1.2 適應(yīng)度函數(shù)的構(gòu)造

采用罰函數(shù)法將有約束條件的優(yōu)化問題轉(zhuǎn)化為無約束的優(yōu)化問題。若迭代點(diǎn)滿足約束，則不予懲罰；否則，若不滿足約束條件，懲罰較嚴(yán)重，懲罰值以違反約束條件的平方增大。

式中：Mp、Mv分別為壓力、流速懲罰因子。

2 SACS混合算法設(shè)計(jì)

2.1 參數(shù)編碼

（1）管網(wǎng)布置編碼：按照近鄰規(guī)劃原則，逆時(shí)針對(duì)管段進(jìn)行編號(hào)，根據(jù)節(jié)點(diǎn)和供水管段一一對(duì)應(yīng)的關(guān)系，各個(gè)基因位上的數(shù)字代表的就是向節(jié)點(diǎn)供水的管段編號(hào)。

（2）管徑編碼：編碼組合采用市場(chǎng)上標(biāo)準(zhǔn)商用管材的管徑，鑒于管徑組合均為有限的離散值，常用的二進(jìn)制編碼會(huì)導(dǎo)致冗余和非可行解等不足，本文采用整數(shù)編碼的方式。不同管道的管徑進(jìn)行從N+ 1 到M（M為正整數(shù)）進(jìn)行編碼，見表2。

表2 管徑編碼

2.2 布谷鳥搜索操作

步驟1：種群初始化，在D維解空間中隨機(jī)生成N個(gè)鳥巢，計(jì)算所有鳥巢的適應(yīng)度值，并保留適應(yīng)度值最優(yōu)的鳥巢位置。

步驟2：設(shè)鳥巢i的位置為xi=(xi1,xi2,…xiD) ，1 ≤i≤N，計(jì)算第t次迭代的所有鳥巢的適應(yīng)度值，最優(yōu)的鳥巢位置保留到下一代，按式（9）進(jìn)行所有鳥巢位置的更新：

式中：為鳥巢i第t次迭代的位置；N(0,1)為D維的標(biāo)準(zhǔn)正態(tài)分布；?為點(diǎn)乘；α為步長(zhǎng)因子；step為由Lévy 分布產(chǎn)生的隨機(jī)步長(zhǎng)；為第t次迭代的最優(yōu)鳥巢。

步驟3：設(shè)鳥巢宿主發(fā)現(xiàn)外來鳥蛋的概率為Pa，產(chǎn)生均勻分布隨機(jī)數(shù)r，r∈(0,1)若r＞Pa，則對(duì)所發(fā)現(xiàn)的鳥巢進(jìn)行隨機(jī)變異，以新鳥巢替換舊鳥巢，計(jì)算所有鳥巢的適應(yīng)度值，且保留適應(yīng)度值最優(yōu)的鳥巢位置。

步驟4：記錄歷史最優(yōu)解，若未達(dá)到迭代終止條件，則返回步驟2。

2.3 模擬退火操作

步驟1：設(shè)初始退火溫度為T0、溫度冷卻系數(shù)為C，建立冷卻進(jìn)度表Tk=C T0以及馬爾可夫鏈的長(zhǎng)度L。

步驟2：令T=T0，根據(jù)當(dāng)前解的位置xk，按式（9）迭代產(chǎn)生下一代的解的位置xk+1。產(chǎn)生一個(gè)在(0，1)區(qū)間上均勻分布的隨機(jī)數(shù)β，按式（10）計(jì)算出在給定當(dāng)前解xk和溫度Tk下與Metropolis接受準(zhǔn)則相對(duì)應(yīng)的轉(zhuǎn)移概率P：

若β＜P則，接收新解，令xk=xk+1，否則xk不變。

步驟3：對(duì)步驟2 重復(fù)k次，按照冷卻進(jìn)度表更新當(dāng)前退火溫度Tk+1，進(jìn)行此溫度下的平衡點(diǎn)尋優(yōu)。

步驟4：判斷迭代次數(shù)是否到達(dá)馬爾可夫鏈的長(zhǎng)度L，若次數(shù)到達(dá)則終止算法，否則轉(zhuǎn)至步驟2繼續(xù)執(zhí)行。

2.4 SACS混合算法流程

SACS 混合算法的基本思路是當(dāng)基本布谷鳥（CS）搜索算法陷入局部最優(yōu)時(shí)，隨機(jī)選擇部分鳥巢按式（9）進(jìn)行更新，并進(jìn)行模擬退火（SA）搜索，產(chǎn)生新的鳥巢，然后繼續(xù)采用基本布谷鳥（CS）算法進(jìn)行搜索。SACS 混合算法流程圖見圖1，主要步驟如下：

圖1 SACS 混合算法流程圖

步驟1：相關(guān)參數(shù)初始化。設(shè)種群數(shù)量為N，最大迭代次數(shù)為M，步長(zhǎng)因子為α，發(fā)現(xiàn)概率為Pa。

步驟2：種群初始化。隨機(jī)生成N個(gè)鳥巢，計(jì)算每個(gè)鳥巢的適應(yīng)度值，并保留最優(yōu)鳥巢的位置。

步驟3：對(duì)鳥巢位置進(jìn)行Lévy 迭代。采用式（9）對(duì)鳥巢位置進(jìn)行更新，計(jì)算新一代鳥巢的適應(yīng)度值，并與上次鳥巢進(jìn)行比較，保留更好的鳥巢。

步驟4：模擬退火過程。對(duì)除最優(yōu)鳥巢外剩余的鳥巢全部進(jìn)行模擬退火操作，計(jì)算所有鳥巢的適應(yīng)度值，且保留適應(yīng)度值最優(yōu)的鳥巢。

步驟5：鳥巢宿主以概率Pa發(fā)現(xiàn)外來鳥蛋。產(chǎn)生均勻分布的隨機(jī)數(shù)r，r∈(0,1)，若r＞Pa，則對(duì)所發(fā)現(xiàn)的鳥巢進(jìn)行隨機(jī)變異，以新鳥巢替換舊鳥巢，計(jì)算所有鳥巢的適應(yīng)度值，且保留適應(yīng)度值最優(yōu)的鳥巢位置。

步驟6：記錄歷史最優(yōu)解，若未達(dá)到迭代終止條件，則返回步驟2。

3 模型驗(yàn)證及分析

采用文獻(xiàn)[14]資料數(shù)據(jù)，管道造價(jià)系數(shù)和指數(shù)為a= 0、B= 0.007 5、b= 1.438 4。設(shè)定種群規(guī)模為N= 40、最大迭代次數(shù)M= 200、步長(zhǎng)因子α= 1、發(fā)現(xiàn)概率Pa= 0.1；初始溫度T= 50、溫度冷卻系數(shù)C= 0.96、馬爾可夫鏈長(zhǎng)度L= 50；壓力懲罰因子Mp=5、流速懲罰因子Mv=5，利用MATLAB 編程計(jì)算，得出符合實(shí)際且造價(jià)最小的方案，其中管網(wǎng)布置和管徑優(yōu)化結(jié)果如圖2和圖3所示。最優(yōu)解的管網(wǎng)長(zhǎng)度為3 600 m，投資額為1.79 萬元，與文獻(xiàn)[14]GA 算法的管網(wǎng)優(yōu)化布置圖4 相比，本文SACS 混合算法管長(zhǎng)減少了16.27%，投資節(jié)約了2.13%；與SAGA 混合算法方法、CS 算法相比，雖管長(zhǎng)相等，但管徑優(yōu)化較大，投資分別節(jié)約了3.76%、17.28%，驗(yàn)證了SACS混合算法的有效性，具體結(jié)果見表3。

圖2 SACS混合算法的管網(wǎng)優(yōu)化結(jié)果

圖3 SACS混合算法的管網(wǎng)優(yōu)化布置圖

圖4 文獻(xiàn)[14]SA算法的管網(wǎng)優(yōu)化布置圖

表3 優(yōu)化結(jié)果對(duì)比表

4 石阡縣茶葉園區(qū)灌溉管網(wǎng)優(yōu)化設(shè)計(jì)

4.1 灌區(qū)概況

本文以貴州省石阡縣龍?zhí)链笸筒枞~園區(qū)216 hm2現(xiàn)代水利試點(diǎn)區(qū)噴灌工程項(xiàng)目為例，開展灌溉管網(wǎng)優(yōu)化設(shè)計(jì)。項(xiàng)目區(qū)平均海拔高程857 m，年平均氣溫16.8 ℃，年平均降水量1 105 mm，年平均日照時(shí)數(shù)1 233.2 h，全年無霜期303 d。根據(jù)《灌溉與排水工程設(shè)計(jì)標(biāo)準(zhǔn)》（GB 50288—2018）及項(xiàng)目區(qū)實(shí)際情況，確定灌區(qū)設(shè)計(jì)保證率為80%[15]。項(xiàng)目區(qū)年降水量時(shí)空分布不均，本次灌溉工程為補(bǔ)充性灌溉，根據(jù)《貴州省行業(yè)用水定額》（DB52T—2011）確定茶園灌溉定額為1 800 m3/hm2[16]。

4.2 管網(wǎng)布置

貴州省石阡縣龍?zhí)链笸筒枞~園區(qū)共有12 個(gè)灌水節(jié)點(diǎn)，傳統(tǒng)設(shè)計(jì)方法管網(wǎng)管網(wǎng)布置方案見圖5。采用緊鄰規(guī)劃的原則初步確定共有42 條可能的連接方案見圖6。1～12 個(gè)節(jié)點(diǎn)最小壓力水頭為10 m，管網(wǎng)最低流速為0.5 m/s，最高流速為3 m/s。茶葉園區(qū)灌溉管網(wǎng)各管段長(zhǎng)度見表4，各節(jié)點(diǎn)特征見表5，不同管徑的管道單價(jià)見表6。

表5 茶葉園區(qū)灌溉管網(wǎng)各節(jié)點(diǎn)特征

圖6 茶葉園區(qū)灌溉管網(wǎng)初步連接圖

表4 茶葉園區(qū)灌溉管網(wǎng)各管段長(zhǎng)度

4.3 計(jì)算結(jié)果及分析

對(duì)表6灌溉管道單價(jià)回歸擬合，得到管道造價(jià)系數(shù)和指數(shù)為a= 0、B= 0.004 8、b= 1.968 3。種群規(guī)模、最大迭代次數(shù)、步長(zhǎng)因子、發(fā)現(xiàn)概率；初始溫度、溫度冷卻系數(shù)、馬爾可夫鏈長(zhǎng)度；壓力懲罰因子、流速懲罰因子等參數(shù)設(shè)定與“3模型驗(yàn)證與分析”節(jié)均相同。

表6 不同管徑的管道單價(jià)（1.25 MPa）

采用MATLAB 編程對(duì)本文設(shè)計(jì)的SACS 混合算法進(jìn)行計(jì)算[17]，得出管網(wǎng)布置和管徑優(yōu)化結(jié)果如圖7 和圖8 所示。不同方案下管網(wǎng)投資比較見表7，傳統(tǒng)設(shè)計(jì)方法管長(zhǎng)1 661.35 m，管網(wǎng)投資為14.460 506 萬元；SACS 混合算法優(yōu)化后的管網(wǎng)長(zhǎng)度為1 649.27 m，投資為11.867 810 萬元。管網(wǎng)總長(zhǎng)度減少了12.08 m；管網(wǎng)投資節(jié)約了17.93%，灌溉工程優(yōu)化效果明顯。

圖7 茶葉園區(qū)SACS混合算法管網(wǎng)布置優(yōu)化結(jié)果

圖8 茶葉園區(qū)SACS混合算法管網(wǎng)優(yōu)化布置圖

表7 茶葉園區(qū)不同方案下管網(wǎng)投資

5 結(jié) 論

本文針對(duì)樹狀灌溉管網(wǎng)投資節(jié)約的設(shè)計(jì)需求，在滿足灌水均勻度、灌溉施肥制度、輪灌制度等一系列設(shè)計(jì)準(zhǔn)則的前提下，設(shè)計(jì)了基于整數(shù)編碼的SACS 混合算法，對(duì)樹狀管網(wǎng)管徑及布置進(jìn)行了優(yōu)化，并結(jié)合文獻(xiàn)[14]資料數(shù)據(jù)與GA 算法、SAGA 混合算法、CS 算法進(jìn)行了對(duì)比分析，表明本文SACS 混合算法比GA算法管長(zhǎng)減少了16.27%，投資節(jié)約了2.13%；比SAGA混合算法、CS算法相比，雖管長(zhǎng)相等，但管徑得到較大的優(yōu)化，投資分別節(jié)約了3.76%、17.28%，驗(yàn)證了SACS 混合算法的有效性。

應(yīng)用本文設(shè)計(jì)的SACS 混合算法，對(duì)石阡縣龍?zhí)链笸筒枞~園區(qū)216 hm2現(xiàn)代水利試點(diǎn)區(qū)工程輸配水管網(wǎng)進(jìn)行優(yōu)化，并與傳統(tǒng)設(shè)計(jì)方法進(jìn)行了對(duì)比，結(jié)果表明：優(yōu)化后的管網(wǎng)長(zhǎng)度為1 649.27 m，管網(wǎng)總長(zhǎng)度減少了12.08 m。由于管徑優(yōu)化效果明顯，故管網(wǎng)總投資由14.460 506 萬元降低為11.867 810 萬元，投資節(jié)約了17.93%，實(shí)現(xiàn)了灌溉管網(wǎng)建設(shè)投資降低的目標(biāo)，為灌溉工程優(yōu)化設(shè)計(jì)提供了新的思路。