李承龍，楊 侃，蔡玉棽，薛 晴，陳 靜，楊晶晶

（1.河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098； 2.宜興市水利農(nóng)機(jī)局，江蘇 宜興 214207）

1 引言

忻州-陽泉地區(qū)位于我國嚴(yán)重缺水的山西省，水資源供需矛盾十分突出，如何將該地區(qū)緊缺的水資源進(jìn)行合理高效調(diào)配是一個亟待解決的問題［1］。 忻州-陽泉地區(qū)作為山西大水網(wǎng)規(guī)劃的第二橫區(qū)域，對整個山西省水網(wǎng)的連通有至關(guān)重要的作用。 目前對于忻州-陽泉地區(qū)的研究較少，對該區(qū)域進(jìn)行水資源優(yōu)化調(diào)配研究有重要的現(xiàn)實(shí)意義。

目前，水資源優(yōu)化調(diào)配模型的求解多采用遺傳算法［2］、模擬退火算法［3］、蟻群算法［4］、煙花算法［5］等智能優(yōu)化算法，但以上算法存在參數(shù)過多導(dǎo)致求解過程過于煩瑣、易早熟、局部尋優(yōu)能力差等問題［6］。 而粒子群算法具有抗干擾性好、參數(shù)少、容易實(shí)現(xiàn)等特點(diǎn)，在圖像處理、模式識別、大地測量等方面得到了廣泛應(yīng)用。 筆者對粒子群算法慣性權(quán)重系數(shù)和適應(yīng)度函數(shù)進(jìn)行了改進(jìn)，使其更加適用于以區(qū)域供水、發(fā)電和生態(tài)為目標(biāo)建立的水資源優(yōu)化調(diào)配模型的求解。 本文以忻州-陽泉地區(qū)為研究對象，建立多目標(biāo)水資源優(yōu)化調(diào)配模型，并使用改進(jìn)的粒子群算法求解，以期為研究區(qū)水資源優(yōu)化調(diào)配提供參考。

2 多目標(biāo)水資源優(yōu)化調(diào)配模型

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

（1）供水目標(biāo)：區(qū)域年供水量最大。 表達(dá)式為

式中：G為區(qū)域年供水量；Gj，t為j單元在t時段的供水量；J為供水單元數(shù)；T為時段數(shù)。

（2）生態(tài)目標(biāo)：水庫下游河道生態(tài)流量滿足程度最大。 表達(dá)式為

式中：E為年內(nèi)下游河道生態(tài)流量之和；ei，t為i斷面t時段的適宜生態(tài)流量；qi，t為i斷面t時段的生態(tài)流量；I為生態(tài)流量監(jiān)測斷面數(shù)。

（3）發(fā)電目標(biāo)：年發(fā)電出力最大。 表達(dá)式為

式中：N為水庫群年出力；Ak為第k個水庫的出力系數(shù)；Qk，t為第k個水庫t時段的發(fā)電流量；Hk，t為第k個水庫t時段的發(fā)電水頭；K為水庫數(shù)量。

2.2 約束條件

（1）水量平衡約束：

式中：Vk，t、Ik，t、qk，t、Ok，t分別為第k個水庫在t時段的蓄水量、上游來水流量、區(qū)間來水流量、出庫流量；Vk，t-1為第k個水庫在t -1 時段的蓄水量；Δt為時段差；Wk，t為第k個水庫在t時段的供水量。

（2）水庫水位約束：

式中：Zk，t為第k個水庫在t時段的水位；Zk，min和Zk，max分別為第k個水庫的最低水位（取死水位）和最高水位（非汛期取正常高水位、汛期取汛限水位）。

（3）下泄能力約束：

式中：q（k，t）為第k個水庫在t時段的下泄流量；qmin（k，t）和qmax（k，t）分別為第k個水庫在t時段的最小和最大下泄流量。

（4）出力約束：

式中：Nk為第k個水庫的出力；Nk，p和Nk，max分別為第k個水庫的保證出力和裝機(jī)容量。

（5）水庫蓄水量約束：

式中：Vk，t為第k個水庫t時段的蓄水量；Vk，min和Vk，max分別為第k個水庫的最小和最大蓄水量。

（6）供水限制約束：

式中：Xj，t和Gj，t分別為j單元在t時段的需水量和供水量。

3 改進(jìn)的粒子群算法

3.1 粒子群算法

粒子群算法是美國學(xué)者Eberhart 和Kennedy 在20世紀(jì)90年代通過觀察鳥群飛行覓食提出的仿生智能算法［7］。 其尋優(yōu)思想：在一個多維的空間中生成鳥群，鳥群擁有空間位置屬性和移動速度屬性，同時鳥與鳥之間能夠相互交流并且有記憶能力，每一次飛行都能夠通過信息共享和記憶的最優(yōu)位置來確定下一次飛行的位置。 在算法當(dāng)中，鳥群（即算法中的粒子群）能夠通過記錄的自身最優(yōu)解和群體最優(yōu)解來迭代更新自己的位置，從而找到空間內(nèi)的最優(yōu)解。

在d維空間中，第n個粒子的位置、速度、個體最優(yōu)值和群體最優(yōu)值表達(dá)式分別為

式中：Xn為第n個粒子的位置；Xn，d為第n個粒子在d維的位置；Vn為第n個粒子的速度；Vn，d為第n個粒子在d維的速度；pn，d為第n個粒子的d維適應(yīng)度值；gn為第n個粒子的個體最優(yōu)值；gbest為群體最優(yōu)值。

在粒子更新過程中，其在d維的位置、速度變化會限制在一定的空間內(nèi)，即當(dāng)粒子在d維內(nèi)超出位置邊界或速度邊界時，該粒子的位置或速度將更新至該維的位置邊界或速度邊界。

在算法迭代過程中，第n個粒子在d維的速度更新公式為

第n個粒子在d維的位置更新公式為

式中：l為迭代次數(shù)；為時間步長，分別為第l次和l -1次迭代時第n個粒子在d維的速度；分別為第l次和l-1次迭代時第n個粒子在d維的位置；c1、c2為學(xué)習(xí)因子，分別調(diào)節(jié)個體最優(yōu)值和群體最優(yōu)值的最大步長，控制粒子在迭代過程中尋找最優(yōu)位置的方向；r1、r2為在［0，1］區(qū)間的隨機(jī)數(shù)，表示搜索的隨機(jī)性；Ω為慣性權(quán)重，控制對解空間的搜索范圍，通常為定值或線性變化。

3.2 改進(jìn)的慣性權(quán)重

慣性權(quán)重為定值或線性變化時，易導(dǎo)致粒子過早陷入局部最優(yōu)解。 為了解決該問題，筆者使用S 型遞減函數(shù)［8］代替?zhèn)鹘y(tǒng)的定值和線性遞減函數(shù)計(jì)算慣性權(quán)重，這樣可使粒子在搜索前期以較快速度對全局進(jìn)行搜尋而避免過早陷入局部最優(yōu)解、全局搜索之后以較慢的局部搜索速度收斂至最優(yōu)值。 改進(jìn)的慣性權(quán)重計(jì)算公式為

式中：Ωmax和Ωmin分別為慣性權(quán)重的最大值和最小值，一般取Ωmax=0.9、Ωmin=0.4；lmax為粒子的最大迭代次數(shù)；α為控制系數(shù)，用于調(diào)節(jié)慣性權(quán)重變化速度的快慢，一般取10。

3.3 定權(quán)重目標(biāo)比適應(yīng)度函數(shù)

區(qū)域年總需水量XY、水庫下游河道年生態(tài)流量總和EY、各水庫按裝機(jī)容量運(yùn)行的出力之和NY計(jì)算公式分別為

式中：Nk，t為第k個水庫t時段按裝機(jī)容量運(yùn)行的出力。

為解決所構(gòu)建的模型各個目標(biāo)的計(jì)算結(jié)果數(shù)量級相差較大的問題，將供水目標(biāo)、生態(tài)目標(biāo)和發(fā)電目標(biāo)相關(guān)指標(biāo)數(shù)值進(jìn)行歸一化處理，使其取值范圍均限制在［0，1］區(qū)間：

式中：F1、F2、F3分別為供水目標(biāo)適應(yīng)度、生態(tài)目標(biāo)適應(yīng)度、發(fā)電目標(biāo)適應(yīng)度。

構(gòu)建的定權(quán)重目標(biāo)比適應(yīng)度計(jì)算公式為

式中：a、b、c分別為F1、F2、F3占適應(yīng)度F的比例，取a =0.4、b =0.3、c =0.3。

4 粒子群算法參數(shù)設(shè)置及求解步驟

4.1 參數(shù)設(shè)置

在本研究中，筆者編碼了100 個維度為26 的粒子，每個粒子都由代表水位的向量X和代表水位變幅的向量V兩個參數(shù)構(gòu)成，其中：X向量由26 個水位組成，前13 個水位為第一個水庫的，后13 個水位為第二個水庫的。 各粒子的X和V編碼為

經(jīng)過多次嘗試和調(diào)整，式（14）中學(xué)習(xí)因子c1、c2均取2，式（16）中最大迭代次數(shù)lmax取300。

4.2 改進(jìn)的粒子群算法求解步驟

（1）根據(jù)模型實(shí)例需要，確定粒子群的規(guī)模、維度、范圍和速度區(qū)間等和粒子群相關(guān)的參數(shù)，以及結(jié)束迭代的條件。

（2）在步驟（1）確定的范圍內(nèi)生成初代粒子群，其中的粒子都有起始速度和起始位置。

（3）根據(jù)目標(biāo)函數(shù)值計(jì)算粒子群中各粒子的適應(yīng)度，同時記錄各粒子的個體最優(yōu)值gn、群體最優(yōu)值gbest，并更新個體最優(yōu)值對應(yīng)的d維適應(yīng)度值pn，d和群體最優(yōu)值gbest。

（4）利用式（14）和式（15）更新下一代粒子的速度和位置。

（5）根據(jù)所設(shè)置的結(jié)束迭代條件判斷是否結(jié)束迭代，若已達(dá)到結(jié)束條件，則輸出結(jié)果；否則返回步驟（3），再次進(jìn)行適應(yīng)度計(jì)算及搜索。

粒子群算法流程見圖1。

圖1 粒子群算法流程

5 實(shí)例應(yīng)用

5.1 研究區(qū)概況與模型概化

忻州-陽泉地區(qū)屬滹沱河流域，包括忻州市忻府區(qū)、定襄縣、原平市、代縣、繁峙縣、五臺縣以及陽泉市、晉中市昔陽縣。 選擇研究區(qū)的坪上水庫、龍華口水庫作為梯級水庫調(diào)度的研究對象。 坪上水庫位于滹沱河干流上，壩高109 m，壩頂長494 m，控制流域面積11996 km2，正常蓄水位678.0 m，總庫容3.11 億m3，是一座以防洪、供水、灌溉為主的大（2）型水庫；龍華口水庫位于陽泉市盂縣北部滹沱河支流龍華河上，大壩為碾壓混凝土重力壩，壩高66 m，正常蓄水位571.6 m，總庫容3001 萬m3，其中發(fā)電庫容2018 萬m3，是一座以發(fā)電、供水、灌溉為主的中型水庫。

為方便計(jì)算，根據(jù)河流水系和實(shí)際供水情況將研究區(qū)概化為定襄縣、盂縣、郊區(qū)、城區(qū)、礦區(qū)、平定縣6個受水區(qū)，上游來水先經(jīng)過定襄縣受水區(qū)后進(jìn)入坪上水庫、龍華口水庫，經(jīng)過調(diào)度之后再向下游盂縣、郊區(qū)、城區(qū)、礦區(qū)、平定縣受水區(qū)供水，供水時考慮發(fā)電效益和河道生態(tài)保護(hù)。

5.2 研究區(qū)需水量及適宜生態(tài)流量

鑒于預(yù)測周期較長，選擇定額法進(jìn)行需水量預(yù)測。以2020年為現(xiàn)狀水平年，根據(jù)忻州-陽泉地區(qū)各產(chǎn)業(yè)和生態(tài)用水量，對規(guī)劃水平年2030年在豐水年、平水年、枯水年來水頻率下6 個受水區(qū)的需水量進(jìn)行預(yù)測，結(jié)果見表1。

表1 忻州-陽泉地區(qū)2030年需水量預(yù)測結(jié)果 萬m3

由表1可以看出，定襄縣、郊區(qū)、礦區(qū)、平定縣4 個受水區(qū)農(nóng)業(yè)用水量較大，原因是這4 個受水區(qū)農(nóng)作物種植面積較大，產(chǎn)業(yè)多以農(nóng)業(yè)為主。 農(nóng)業(yè)用水量隨來水頻率的變化而變化，總體表現(xiàn)為枯水年＞平水年＞豐水年，原因是枯水年降水量小且蒸發(fā)量大，需從水庫向農(nóng)業(yè)供給更多的水量。 生產(chǎn)、生活和河道外生態(tài)用水量受降水和天然來水的影響較小。 2030年，忻州-陽泉地區(qū)豐水年、平水年、枯水年總需水量分別為27974萬、35217 萬、36735 萬m3。

河道內(nèi)生態(tài)需水量計(jì)算方法有很多種［4］，本文選擇最為簡便的90%保證率最枯月平均流量法，計(jì)算兩個生態(tài)監(jiān)測斷面的河道適宜生態(tài)流量，結(jié)果見表2。

表2 生態(tài)監(jiān)測斷面的河道適宜生態(tài)流量計(jì)算結(jié)果 m3／s

5.3 優(yōu)化調(diào)度結(jié)果及分析

分別使用傳統(tǒng)粒子群算法和改進(jìn)粒子群算法，對忻州-陽泉地區(qū)水資源優(yōu)化調(diào)配模型進(jìn)行多次尋優(yōu)迭代求解，結(jié)果見表3?？梢钥闯觯涸?0次求解過程中，傳統(tǒng)粒子群算法在第1次、第3次、第4次、第5次、第8次求解時，在迭代60次以內(nèi)就收斂到局部最優(yōu)值；而改進(jìn)的粒子群算法均在迭代100次左右找到最優(yōu)解，且除了第4次求解陷入局部最優(yōu)值0.819 外，最優(yōu)解均分布在0.851 ～0.856 之間，說明改進(jìn)的粒子群算法全局搜索能力更強(qiáng)，不易陷入局部最優(yōu)解。

表3 尋優(yōu)迭代次數(shù)和收斂值對比

從多次求解結(jié)果中選出最優(yōu)調(diào)配方案的各目標(biāo)優(yōu)化調(diào)度結(jié)果，見表4。可以看出，豐水年受水區(qū)2030年總供水量為27974 萬m3，各地區(qū)的需水滿足率達(dá)100%；平水年定襄縣、礦區(qū)、平定縣出現(xiàn)不同程度缺水，原因是平水年降水減少導(dǎo)致需從水庫、河道供水的農(nóng)業(yè)灌溉用水增加，同時上游來水減少，進(jìn)而使得這3個農(nóng)業(yè)灌溉需水量較大的受水區(qū)出現(xiàn)不同程度的缺水。 因此，農(nóng)業(yè)需水量較高的區(qū)域會出現(xiàn)平水年缺水情況，需水滿足率相對較低，枯水年各區(qū)域需水滿足程度有差異也是同樣原因；枯水年需水量為36735 萬m3，而可供水量僅為24079 萬m3，缺水量為12656 萬m3，出現(xiàn)了較為嚴(yán)重的供水不足情況；需水量豐水年＜平水年＜枯水年。

表4 忻州-陽泉地區(qū)2030年豐水、平水、枯水年供、需水量對比

由表5、表6可知，研究區(qū)豐水年、平水年、枯水年發(fā)電量分別為25373 萬、22443 萬、20238 萬kW·h。除6月外豐水年因上游來水較多且供水能得到保障，水庫蓄水量充盈，發(fā)電水頭高和發(fā)電下泄流量充足，故生態(tài)流量基本可以滿足、發(fā)電目標(biāo)基本可以實(shí)現(xiàn)。 平水年隨著受水區(qū)出現(xiàn)一定的缺水，存在一定程度的河道生態(tài)缺水和發(fā)電量降低的問題。 枯水年來水量最少，河道生態(tài)滿足率和發(fā)電量最低。

表5 忻州-陽泉地區(qū)2030年豐水、平水、枯水年生態(tài)流量優(yōu)化調(diào)度結(jié)果

表6 忻州-陽泉地區(qū)2030年豐水、平水、枯水年發(fā)電量優(yōu)化調(diào)度結(jié)果

6 結(jié)論

（1）利用本文建立的多目標(biāo)水資源優(yōu)化調(diào)配模型，進(jìn)行規(guī)劃水平年2030年忻州-陽泉地區(qū)水資源優(yōu)化調(diào)配，得出研究區(qū)豐水年供水、發(fā)電、生態(tài)滿足率均最大，平水年次之，枯水年最小。 平水年定襄縣、礦區(qū)、平定縣出現(xiàn)不同程度缺水，枯水年出現(xiàn)較為嚴(yán)重的供水不足。

（2）改進(jìn)粒子群算法模型求解時不易陷入局部最優(yōu)解，能夠收斂到全局最優(yōu)的結(jié)果，說明使用S 型遞減慣性權(quán)重和定權(quán)重比適應(yīng)度函數(shù)相結(jié)合的方法能夠改善傳統(tǒng)粒子群算法的求解能力。

雖然水資源優(yōu)化調(diào)配模型能夠一定程度上優(yōu)化研究區(qū)水資源時空分布，協(xié)調(diào)各目標(biāo)之間的滿足情況，但難以從根本上解決地區(qū)缺水問題，應(yīng)當(dāng)通過建設(shè)調(diào)水引水工程、開展節(jié)水減排、提高水資源利用效率等措施解決水資源短缺問題。