杜佰林，張建豐，高澤海，李濤，黃子奇，張娜

(西安理工大學(xué)省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安 710048)

當(dāng)前，隨著中國區(qū)域經(jīng)濟的快速發(fā)展和城市化進程的加快，水資源開發(fā)不平衡、地下水開采量大、水資源利用效率低等問題日趨嚴峻，造成水資源短缺現(xiàn)象逐漸凸顯[1].同時，區(qū)域降雨量時空分布不均、用水管理水平不高，導(dǎo)致水資源供需矛盾有所惡化，制約了社會經(jīng)濟快速發(fā)展[2].為了緩解上述危機，研究區(qū)域水資源優(yōu)化配置顯得尤為必要，而傳統(tǒng)的配置方法面對水資源系統(tǒng)的復(fù)雜性和不確定性時，往往只進行經(jīng)驗估計或者忽略，致使配置結(jié)果存在偏差.因此，研究如何實現(xiàn)區(qū)域水資源可持續(xù)開發(fā)利用和水資源優(yōu)化配置，具有十分重要的現(xiàn)實意義.

國內(nèi)外眾多學(xué)者針對水資源優(yōu)化配置中模型的建立開展了大量研究工作.國際上利用計算機技術(shù)對水資源系統(tǒng)的復(fù)雜關(guān)系進行了模擬仿真，逐步成為一種被普遍認可的配置思路[3].MASSE[4]針對水資源配置中的單目標(biāo)問題進行了研究，該成果的特點是將配置理論應(yīng)用到了實踐中.ABDULBAKI等[5]提出并解決了水資源配置的多目標(biāo)問題.WADA等[6]建立了不同水質(zhì)下的區(qū)域多水源、多用水戶水資源配置模型.EL-NAQA等[7]在綜合考慮地下水超采、水質(zhì)惡化的基礎(chǔ)上建立了多水源聯(lián)合調(diào)度模型.

中國的水資源分配研究工作以水庫優(yōu)化調(diào)度為基點，經(jīng)歷了水資源的供需平衡配置、經(jīng)濟效益最大化配置、面向生態(tài)效益合理配置這3個階段.粟曉玲等[8]提出了干旱地區(qū)補給水的同時考慮生態(tài)效益的配置供需理念.劉玒玒等[9]應(yīng)用蟻群和粒子群的混合智能優(yōu)化算法解決了同時考慮水質(zhì)和水量的多目標(biāo)水資源優(yōu)化配置模型.張倩等[10]應(yīng)用改進的粒子群算法解決了當(dāng)前灌區(qū)有限農(nóng)業(yè)水資源在不同作物中合理分配的問題.

可以發(fā)現(xiàn)，水資源優(yōu)化配置模型由于自身的多元化呈指數(shù)增長，簡單的線性規(guī)劃求解已不能滿足配置模型的發(fā)展趨勢，甚至一些智能優(yōu)化算法，例如遺傳算法、蟻群算法、粒子群算法等往往由于自身的局限，同樣導(dǎo)致配置結(jié)果不理想.因此，在水資源優(yōu)化配置過程中選擇適宜的模型，克服算法自身尋優(yōu)的缺陷，提高算法收斂精度及效率，是目前配置過程中亟待解決的問題.

鑒于此，文中建立綜合效益相對最優(yōu)的大荔縣水資源優(yōu)化配置模型，以實現(xiàn)該區(qū)域水資源的可持續(xù)開發(fā)利用；通過對粒子群算法(particle swarm optimization, PSO)和模擬退火算法(simulated annea-ling, SA)深入分析，將模擬退火算法在一定概率控制下，暫時接受一些劣質(zhì)解的特性，用于改善基本粒子群算法的局部尋優(yōu)能力，從而獲得具有較好全局尋優(yōu)能力的模擬退火粒子群算法(simulated annealing-particle swarm optimization, SA-PSO)，以期為今后水資源規(guī)劃模型提供新的求解方法，為區(qū)域水資源管理者提供更多的決策依據(jù).

1 模型的建立與求解算法

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

水資源優(yōu)化配置是以公平性、持續(xù)性和高效性為原則，在一些特定的區(qū)域內(nèi)運用一些技術(shù)方法，將有限的水資源科學(xué)合理地分配給各用水部門，以滿足社會、經(jīng)濟和環(huán)境的可持續(xù)發(fā)展，尋求綜合效益的最優(yōu)值.因此構(gòu)建以社會、經(jīng)濟、生態(tài)效益為目標(biāo)的綜合評價函數(shù)，表達式為

Z=opt[Fecon(X),Fsoci(X),Fecol(X)]，

(1)

式中：X為不同賦存形式、不同質(zhì)量和數(shù)量的水資源所構(gòu)成的決策變量；Fecon(X)，F(xiàn)soci(X),Fecol(X)分別為本模型求解的經(jīng)濟效益、社會效益和生態(tài)效益.

1) 經(jīng)濟效益目標(biāo).以水資源開發(fā)利用中不同水平年下各用水戶的直接經(jīng)濟效益最大化作為經(jīng)濟效益目標(biāo)，表達式為

(2)

式中：i為供水水源類別；j為用水戶類別；Qij為水源i向用水戶j的供水量，萬m3；bij為水源i向用水戶j的效益系數(shù)，元/m3；cij為水源i向用水戶j的費用系數(shù)，元/m3；αi為水源i的供水次序系數(shù)；βj為用水戶j的用水公平系數(shù).

2) 社會效益目標(biāo).考慮到社會效益很難量化，而社會的發(fā)展和穩(wěn)定受區(qū)域缺水程度影響較大，且區(qū)域缺水率又可以度量區(qū)域的缺水程度，故以水資源開發(fā)利用中區(qū)域綜合缺水率最小化作為社會效益目標(biāo)，表達式為

(3)

式中：Dj為用水戶j的需水量，萬m3.

3) 生態(tài)效益目標(biāo).以各用水戶在對應(yīng)規(guī)劃水平年中所排放的重要污染因子(chemical oxygen demand, COD)的總量最小化作為生態(tài)效益目標(biāo)，表達式為

(4)

式中：dj為用水戶j排放單位體積廢水中所含重要污染因子的量，mg/L；pij為用水戶j的污水排放系數(shù).

1.2 約束條件

文中研究目的旨在配置各用水戶在不同水源地下的供水量.該供水量在滿足配置目標(biāo)的同時，仍受供水水源、需水能力、環(huán)境能力及各變量非負等約束.約束條件表示如下.

1) 供水能力約束.水源i向用水戶j的供水總量應(yīng)小于其最大可供水量，表達式為

(5)

式中：Wimax為要求水源i的最大可供水量，元/m3.

2) 需水能力約束.配水過程中各用水戶從水源i所得到的水量應(yīng)該結(jié)合自身需水能力，既不小于最低約束也不得超過額定約束，表達式為

(6)

3) 環(huán)境能力約束.各用水戶所排放單位體積廢水中所含重要污染因子的含量應(yīng)在國家允許的排放指標(biāo)內(nèi)；且排放的重要污染因子總量應(yīng)不超過該區(qū)域的最大允許排放量，表達式為

(7)

式中：m0j為在各行業(yè)中國家標(biāo)準(zhǔn)規(guī)定所排放的重要污染物COD質(zhì)量濃度，mg/L；Mmax為區(qū)域內(nèi)最大允許排放重要污染物的總量.

4) 變量非負約束.水源i向用水戶j的供水量要滿足非負約束，表達式為

Qij≥0,i=1,2,3,…;j=1,2,3,….

(8)

1.3 模擬退火粒子群算法

1.3.1 模擬退火算法

模擬退火算法是由METROPOLIS等[11]根據(jù)物理學(xué)中固體物質(zhì)的冷卻退火過程和一般組合優(yōu)化問題求解過程之間的相似性，提出的一種隨機全局優(yōu)化算法.該求解步驟如下.

Step1 初始一個足夠大的溫度T0>0，降溫次數(shù)k=0，隨初始化x(0)，使得x(i)=x(0)，并計算其能量值E(x(0)).

Step3 根據(jù)降溫公式d(Tk)，計算Tk+1=d(Tk),k=k+1，判斷是否滿足Metropolis準(zhǔn)則，若滿足則計算結(jié)束，輸出結(jié)果；否則返回Step2.

1.3.2 粒子群算法

粒子群算法是由KENNDY等[12]受到鳥群覓食行為的啟發(fā)后，提出的一種涉及參數(shù)少、易實現(xiàn)的高效搜索算法，與模擬退火算法相同屬于進化算法，均以迭代的方式探尋最優(yōu)解.該求解步驟如下.

Step1 初始化含有種群規(guī)模N，初始位置popi和速度vi的粒子群體.

Step2 計算每個粒子的適應(yīng)度值fitness(i).

Step3 對每個粒子，比較個體極值Pbest(i)與適應(yīng)度值fitness(i)，若fitness(i)>Pbest(i)，則令Pbest(i)=fitness(i).

Step4 對每個粒子，比較全局極值gbest(i)與適應(yīng)度值fitness(i)，若fitness(i)>gbest(i),則令gbest(i)=fitness(i).

Step5 更新粒子的位置popi和速度vi，即

vid(t+1)=wvid(t)+c1r1[pid-popid(t)]+c2r2[pgd-popid(t)]，

(9)

popid(t+1)=popid(t)+vid(t+1)，

(10)

式中：w為慣性權(quán)重；c1，c2為學(xué)習(xí)因子，分別代表i的個體認知和社會加速常數(shù)；r1，r2為[0,1]相互獨立的隨機數(shù).

Step6 判斷是否滿足終止條件，若符合其精度要求或達到最大迭代次數(shù)，則計算結(jié)束，輸出結(jié)果；否則返回Step2.

1.3.3 模擬退火粒子群算法

粒子群算法雖然能夠有效地優(yōu)化各種函數(shù)，但對離散函數(shù)、多目標(biāo)優(yōu)化等問題，往往容易陷入局部極小值；模擬退火算法則可通過在一定概率下接受一些劣質(zhì)解，從而跳出原狀態(tài)并繼續(xù)搜索.文中在求解過程中將兩者融合，構(gòu)建模擬退火粒子群算法，使算法在搜索過程中有較強的全局尋優(yōu)能力，有效避免搜索陷入局部極小解.其求解步驟如下.

Step1 初始化含有種群規(guī)模N，初始位置popi和速度vi的粒子群體.

Step2 計算每個粒子的適應(yīng)度fitness(i)，比較個體極值Pbest(i)與適應(yīng)度值fitness(i)，若fitness(i)>Pbest(i),則令Pbest(i)=fitness(i)；同理比較全局極值gbest(i)與適應(yīng)度值fitness(i)，若fitness(i)>gbest(i),則令gbest(i)=fitness(i).

Step3 設(shè)定足夠大的初始溫度T0>0并進行相應(yīng)的退溫操作，即

T0=f(pi)/ln 5，

(11)

Tk+1=λTk，

(12)

式中：λ為退火過程中的衰減參數(shù).

Step4 確定當(dāng)前溫度T下每個粒子pi的適應(yīng)值，即

(13)

式中：TF為粒子適應(yīng)值.

Step5 從所有pi中確定其全局最優(yōu)的代替值p′i，并采用式(9)和(10)更新其位置和速度.

Step6 計算確定粒子目標(biāo)值，并更新Pbest(i)和gbest(i)，隨后進行退溫操作.

Step7 判斷是否滿足終止條件，若符合條件，則計算結(jié)束并輸出結(jié)果；否則返回Step4.

該算法流程如圖1所示.

圖1 模擬退火粒子群算法計算流程

2 實例研究

2.1 研究區(qū)域概況

大荔縣(109°43′～110°19′E，34°36′～35°02′N)位于陜西省關(guān)中東部平原區(qū)，地處黃、洛、渭河的交匯帶，并以“三秦通衢，三輔重鎮(zhèn)”而出名，是一座典型的陜晉豫承接產(chǎn)業(yè)轉(zhuǎn)移示范區(qū)，總縣域面積為1 776 km2，屬暖溫帶半干旱大陸性季風(fēng)氣候.全縣年平均氣溫為14.4 ℃，年均降水量為514 mm，年均蒸發(fā)量為968.3 mm，降水量遠小于蒸發(fā)量，很難滿足蒸發(fā)量需求，時有水資源短缺現(xiàn)象發(fā)生.大荔縣水資源供給系統(tǒng)依據(jù)《全國水資源綜合規(guī)劃》概要[13]的有關(guān)要求，具體描述如圖2所示.

圖2 大荔縣水資源供給系統(tǒng)

2.2 數(shù)據(jù)來源

考慮到配置的合理性，在2016現(xiàn)狀年的基礎(chǔ)上選擇在規(guī)劃水平年2020和2030年來水保證率P={50%，75%，95%} 的3種不同來水水平情況下進行水量優(yōu)化配置.統(tǒng)計大荔縣多年外調(diào)水量資料和渭南市“十三五”水利發(fā)展規(guī)劃得出縣域目前水資源總量為67 938.67萬m3.結(jié)合大荔縣現(xiàn)狀庫壩工程、調(diào)水工程及規(guī)劃新建工程的設(shè)計用途、供水對象、工程布局和調(diào)整方案等實際情況，將可供水量按地表水、地下水、其他水等3個供水水源進行統(tǒng)計，即得到不同水平年不同來水保證率下的可供水量Qor；根據(jù)1996—2016年的《大荔縣統(tǒng)計年鑒》、《大荔縣志》、《洛惠渠志》等相關(guān)資料，在獲得區(qū)域社會經(jīng)濟發(fā)展指標(biāo)的基礎(chǔ)上，采用水資源規(guī)劃中定額法、趨勢法進行預(yù)測，從而得到在2020和2030年中P={50%，75%，95%}保證率時，不同用水行業(yè)的需水量Qre具體數(shù)值見表1.

表1 大荔縣不同水平年不同保證率下需水量預(yù)測

2.3 模型參數(shù)確定

2.3.1 模型參數(shù)

依據(jù)上述水資源優(yōu)化配置模型并結(jié)合大荔縣實際用水情況，確定本次水資源配置的供水次序：地表水、地下水、其他水；參考式(14)將供水次序轉(zhuǎn)化為0～1的系數(shù)，即供水次序系數(shù)αi分別為0.50，0.33，0.17.

(14)

式中：I為供水水源的總數(shù)；hi為水源i的供水次序序號；hmax為其供水序號的最大值.

配置過程以“先生活、后生產(chǎn)”為前提條件，在堅持以生活和生態(tài)用水為優(yōu)先原則的基礎(chǔ)上，設(shè)定需水量部門得到用水的順序依次為生活、生態(tài)、工業(yè)、第三產(chǎn)業(yè)和農(nóng)業(yè)用水.參數(shù)設(shè)置參考式(14)，即各部門用水公平系數(shù)βj分別為0.33，0.27，0.20，0.13，0.07.用水效益系數(shù)bij參考《陜西省行業(yè)用水定額》(DB 61/T 943—2013)得出，即生活、生態(tài)、工業(yè)、第三產(chǎn)業(yè)和農(nóng)業(yè)的bij分別為500，500，580，620，2元/m3.據(jù)大荔縣現(xiàn)狀年2016年的水費征收標(biāo)準(zhǔn)確定各部門的費用系數(shù)cij，即生活、生態(tài)、工業(yè)、第三產(chǎn)業(yè)和農(nóng)業(yè)的cij分別為2.58，2.58，3.98，5.76，0.45元/m3.需水系數(shù)rj的設(shè)置參照文獻[14]，采用概率分布函數(shù)和置信區(qū)間，設(shè)置置信水平為0.90時，rj～N(0.85,0.382)，即生活和生態(tài)的rj為1，工業(yè)、第三產(chǎn)業(yè)和農(nóng)業(yè)的rj均服從rj～N(0.85,0.382).大荔縣污水排放主要表現(xiàn)為城鎮(zhèn)生活和工業(yè)廢水排放，故污水排放系數(shù)pij統(tǒng)一參照《城市排水工程規(guī)劃規(guī)范》(DB 50318—2000)為0.75.目標(biāo)函數(shù)權(quán)重的設(shè)置，結(jié)合水資源系統(tǒng)中不確定性的特點，在采用概率分布的基礎(chǔ)上參照文獻[15]中熵權(quán)法的相關(guān)知識，計算得社會、經(jīng)濟、生態(tài)效益的目標(biāo)函數(shù)權(quán)重依次為0.5，0.3，0.2.

2.3.2 算法參數(shù)

以Matlab R2018b為開發(fā)工具編寫模擬退火粒子群算法程序，依據(jù)前人經(jīng)驗[16]并經(jīng)過多次試算求解，設(shè)定粒子種群數(shù)目N為1 000，慣性權(quán)重ω為0.729 8，學(xué)習(xí)因子c1=c2=2，r1和r2均為[0，1]的隨機數(shù)；設(shè)定模擬退火初始溫度T0為1 000，終止溫度Tf為1，衰減參數(shù)λ為0.6，最大迭代次數(shù)maxTit為1 000.

3 結(jié)果與分析

3.1 最優(yōu)配水量

利用所建立的基于模擬退火粒子群算法的優(yōu)化配置模型對大荔縣進行水資源優(yōu)化配置，得出大荔縣在不同規(guī)劃水平年、不同來水保證率下的水資源優(yōu)化配置結(jié)果Qop，結(jié)合水資源供需平衡特征，對水資源配置結(jié)果進行分析，以P=75%為例，分析結(jié)果見表2，表中Qor，Qre，rws分別為供水量、需水量、缺水率.

由表2可知，到2020年和2030年，在“優(yōu)先保證生活和生態(tài)用水”的前提下，大荔縣生活和生態(tài)用水均得到了完全滿足，工業(yè)用水缺水情況極少，缺水情況主要表現(xiàn)在農(nóng)業(yè)用水和第三產(chǎn)業(yè)用水，符合縣域現(xiàn)狀.優(yōu)化配置結(jié)果表明，2020年大荔縣農(nóng)業(yè)缺水率為3.24%，第三產(chǎn)業(yè)缺水率為0.07%，整體缺水率為2.00%；2030年大荔縣的農(nóng)業(yè)缺水率為3.57%，第三產(chǎn)業(yè)缺水率為0.13%，整體缺水率為1.67%，缺水情況仍然存在.但以上缺水率均未超過4%，可通過采用一定農(nóng)藝措施或調(diào)整農(nóng)業(yè)種植結(jié)構(gòu)，促使該行業(yè)用水供需平衡.配置結(jié)果符合規(guī)劃模型中約束條件的設(shè)置，滿足5類不同用水戶的用水需求，配置結(jié)果合理.同時，研究區(qū)整體缺水率降低，一定程度上提高了各行業(yè)用水需求，這也表明基于模擬退火粒子群算法的優(yōu)化配置模型具有可行性.

表2 大荔縣不同水平年P(guān)=75%保證率下的水資源優(yōu)化配置結(jié)果

3.2 配置結(jié)果

為了說明優(yōu)化配置后的水量Qop與原始供水量Qor之間的差異，在模型求解得到地表水、地下水和其他水源在不同規(guī)劃水平年、不同來水保證率下配置水量的基礎(chǔ)上，將大荔縣不同水平年、不同保證率下供水量對比結(jié)果進行分析，具體分析結(jié)果如表3所示，表中Q為原始供水量與優(yōu)化配置后水量的差值.

表3 大荔縣不同水平年不同保證率下供水量對比結(jié)果

由表3可知，優(yōu)化配置后的水量較原始供水量均有遞減的趨勢，其中2020年在P=50%，75%和95%的3種不同來水水平情況下，優(yōu)化配置后水量較原始供水量的總節(jié)水量分別為372.66，333.11，255.14萬m3；2030年在P=50%，75%和95%的3種不同來水水平情況下，優(yōu)化配置后水量較原始供水量的總節(jié)水量分別為426.63，404.76和352.85萬m3.總節(jié)水量均超過250萬m3，節(jié)水效果明顯，達到了節(jié)約用水的目的，可滿足渭南市水資源綜合規(guī)劃要求.

對大荔縣進行水資源優(yōu)化配置后的3個目標(biāo)函數(shù)效益值進行計算和分析，以P=75%為例，分析結(jié)果見表4.與2016現(xiàn)狀年相比，大荔縣在2020和2030年2個規(guī)劃水平年中的經(jīng)濟效益fecon(指各用水戶的直接經(jīng)濟效益)、社會效益fsoci(以區(qū)域綜合缺水率表示，即各用水戶的缺水量/需水量的百分數(shù))和生態(tài)效益fecol(指各用水戶所排放的重要污染因子COD的總量)改善均比較明顯，配置結(jié)果合理，驗證了基于模擬退火粒子群算法優(yōu)化配置模型具有可行性.

表4 大荔縣不同水平年P(guān)=75%保證率下的效益值

3.3 模型精度

為了檢驗?zāi)M退火粒子群算法的可靠性，引入經(jīng)典測試函數(shù)——Ackley函數(shù)[17]：

-8≤xi≤8.

(15)

以Matlab R2018b為開發(fā)工具，編寫Ackley函數(shù)程序并繪制三維圖像，獲取Ackley函數(shù)的空間幾何特性.圖3為Ackley函數(shù)圖.該函數(shù)較為復(fù)雜，是一個n維函數(shù)并有多個局部極小值點，且fmin(0,0)=0，可用于檢測一個算法的全局收斂速度.分別采用SA，PSO和SA-PSO 的3種進化算法依次迭代計算測試函數(shù)，從而評價比較3種進化算法的性能，整理對比后結(jié)果見表5，表中Titf和Tits分別為最快、最慢迭代次數(shù).由表可知3種進化算法都能在有限的Tit內(nèi)找到其最優(yōu)解，其中SA的收斂速度最慢，為327次；SA-PSO的收斂速度最快，僅為15次，優(yōu)化效率顯著提升.

圖3 Ackley函數(shù)圖

表5 3種進化算法的性能對比

為了進一步表明模擬退火粒子群算法全局尋優(yōu)的優(yōu)越性，在Matlab R2018b中對性能差異較小的基本PSO和SA-PSO的個體最優(yōu)適應(yīng)值Fva和迭代次數(shù)Tit進行了仿真, 仿真結(jié)果如圖4所示.

根據(jù)圖4可以得出：SA-PSO的性能較基本PSO有了一定程度的提高，在初始條件基本相同的情況下，SA-PSO由于結(jié)合了SA思想中一定概率控制下暫時接受一些劣質(zhì)解的特性，克服了基本PSO容易陷入局部最優(yōu)解的缺點，使其全局尋優(yōu)能力增強，進而在精度和進化速度兩方面都要優(yōu)于基本PSO，且SA-PSO原理通俗易懂、容易實現(xiàn)，可為今后水資源規(guī)劃模型提供一個新的求解思路.

圖4 大荔縣2020，2030年優(yōu)化配置的進化過程

4 結(jié) 論

構(gòu)建以社會、經(jīng)濟、生態(tài)效益為目標(biāo)的綜合評價函數(shù)，建立基于模擬退火粒子群算法的水資源優(yōu)化配置模型；結(jié)合大荔縣進行實例驗證，證明算法的可行性，實現(xiàn)大荔縣的不同行業(yè)用水量的合理優(yōu)化配置.

1) 通過將模擬退火算法在一定概率控制下，暫時接受一些劣質(zhì)解的特性，用于改善基本粒子群算法的局部尋優(yōu)能力，有效克服了單一粒子群算法和模擬退火算法相應(yīng)的缺點，避免粒子陷入局部極小值，提高了算法的精度，加快了算法的進化速度.

2) 構(gòu)建大荔縣水資源優(yōu)化配置模型，并采用模擬退火粒子群算法進行求解，通過粒子編碼、設(shè)置約束條件等，求解模型得到大荔縣不同行業(yè)用水量的優(yōu)化配置方案.優(yōu)化后節(jié)水效果明顯，經(jīng)濟效益增長顯著，可滿足大荔縣的可持續(xù)發(fā)展要求.

3) 通過對比3種進化算法計算測試函數(shù)的結(jié)果，比較個體最優(yōu)適應(yīng)值的求解結(jié)果，發(fā)現(xiàn)模擬退火粒子群算法能夠獲得更好的個體適應(yīng)度值，收斂速度更快，整體運行更穩(wěn)定，為區(qū)域水資源優(yōu)化配置提供了一種新的解決途徑.