朱思峰,李子胥
(天津城建大學計算機與信息工程學院,天津 300384)
我國水資源呈現時空分布不均,人均占有量小等特點[1]。由于水資源的自然因素、經濟因素以及全球氣候變化等自然問題的綜合影響,我國面臨著水資源供應不足、水環境污染、水生態退化以及極端與突發事件頻繁等突出的問題。為了實現可持續發展,需要對區域水資源進行科學合理的優化配置[2]。
在水資源總量有限的前提下,通過對水資源進行優化配置,可以更好地滿足社會用水需求。目前,水資源優化配置成為了一個研究熱點,引起了許多學者的關注。對于復雜的水資源優化配置問題,傳統的數學方法難以求解,出現了許多基于啟發式智能算法的水資源優化配置方案:文獻[3]通過使用快速非支配排序策略和增加擁擠度比較算子來改進飛蛾火焰算法,并在三亞市的水資源配置模型中得到了缺水量最小的決策方案;文獻[4]通過增加Logistic映射和慣性權重來改進鯨魚算法,收斂速度和精度有明顯提升,并應用于邯鄲市的水資源優化配置模型中,得到了以缺水量最小為特殊偏好的解集;文獻[5]采用模擬退火粒子群算法從各用戶配水量、缺水程度、配置目標及模型精度等多方面進行分析,得出了陜西省大荔縣各子區域在不同條件下水資源優化配置結果;文獻[6]采用基于改進MFO算法求解水資源優化配置,通過自適應權重表示改進MFO算法中的燭火,以自適應權重數值越小為模型最優解,得到了凈收益率更高的水資源優化配置解集;文獻[7]通過使用有并行機制和全局優化特性的遺傳算法,將水資源優化配置問題建模為生物進化問題,通過種群間的優勝劣汰生成最優解集,以完成水資源的優化配置;文獻[8]通過使用多智能體Q-學習算法將城市不同區域的用水戶抽象為基于agent的模型,同時采用了一種自適應的獎勵值函數來提高多智能體Q學習算法的性能以處理城市水資源配置中的各種情況;文獻[9]從慣性因子及學習因子選擇、外部檔案維護和全局最優選取策略3個方面改進了多目標粒子群算法,并對宿遷市水資源優化調度進行實例研究,得到了分布性更好的解集。
上述這些文獻大多是通過使用智能優化算法對水資源優化配置模型進行求解,為水資源優化配置提供了全新的思路。第二代非支配排序遺傳算法(Non-Dominated Sorted Genetic Algorithm-II,NSGA-Ⅱ)是當前解決多目標優化問題最有效的方法之一[10]。NSGA-Ⅱ算法是 Srinivas和 Deb于2000年在NSGA 的基礎上提出的[11],它具有良好的收斂性和分布性,在許多領域被廣泛應用,并有很好的表現。但NSGA-Ⅱ沒有考慮個體擁擠度相同時的情況,同時未考慮淘汰某個個體后相鄰個體擁擠度變化造成的影響,因此提出了采用一種新型的擁擠度計算策略和一種動態的擁擠度計算方式的改進型NSGA-Ⅱ算法(簡稱為INSGA-II算法),并把其應用于求解建立的水資源優化配置模型。得到的非支配解集可以為多目標優化配置提供可選的解決方案,同時可為區域水資源的可持續發展提供決策依據。
山西水資源十分貧乏。據數據統計,山西省人均水資源在全國列為倒數第二位。以1984年所統計的人口耕地面積計算,山西地區人均占水量466 m3,占全國人均占水量的18.8%,山西地區單位面積平均水量占全國的8%~10.9%,與南方各省相比差距較大,這足以說明山西水資源問題極其嚴峻[12]。
研究區域為山西省晉中市南部,行政區域劃分為3個子區,分別為G1子區(平遙縣)、G2子區(榆社縣)、G3子區(左權縣),依次對應k=1,2,3;共有3種水源,分別為地表水、地下水及引調水,依次對應i=1,2,3;共有5種用戶,分別為生活用水、農業用水、第二產業用水、第三農業用水和生態用水,依次對應j=1,2,3,4,5。以最小區域內的缺水量和最大化供水經濟效益為目標,建立了晉中南部供水區水資源優化配置模型。晉中市南部水資源配置網絡圖如圖1。

圖1 晉中市南部水資源配置網絡圖Fig.1 Water resources allocation network in the south of Jinzhong City
以2025年為規劃水平年,在來水頻率P=75%下進行水資源的優化配置,以滿足缺水量最小化和經濟效益最大化的目標。依據晉中市的發展規劃,環境最小需水量通過參考10年來枯水月的平均流量來計算,然后對不同子區內的5類用水戶的需水量進行預測,不同子區內的用水戶的需水量預測結果[13]如表1所示。

表1 晉中市子區2025年需水量預測 萬m3Tab.1 Water demand forecast of Jinzhong sub district in 2025
考慮到實際情況中晉中市的可供水總量的不確定性、不同年份來水頻率不同和晉中市近期的實際調水工程布置政策,對2025年晉中市不同子區在來水頻率P=75%的情況下的可供水量進行預測[13],預測結果如表2所示。

表2 晉中市子區2025年供水預測結果 萬m3Tab.2 Forecast results of water supply in 2025 for southern Jinzhong water supply area
多目標水資源優化配置模型是綜合考慮社會效益和經濟效益,對區域內水資源進行合理有效的配置,提高區域內水資源的利用率的模型。模型將社會效益和經濟效益作為目標函數、對用水戶的供水量作為決策變量進行優化,其中社會效益函數用所有子區K用水戶J的需水量減去所有子區K的水源對用水戶J的供水量之和來表示,目標函數取缺水量(下同)最小值。經濟效益函數為用水戶J用水產值和用水費用之差,目標函數取最大值。
水資源優化配置的目標是使區域內的水資源配置的缺水量最小、經濟效益最大。將晉中市南部區域分為3個子區,設第k(1≤k≤3)個子區有3個水源,5個用水戶。決策向量其中表示第i(1≤i≤3)個水源可為子區k的j(1≤j≤5)用水戶提供的水量。
(1)社會效益目標函數。以缺水量最小來滿足社會效益目標函數,設第k(1≤k≤3)個子區的第j(1≤j≤5)個用水戶的需水量用來表示。
社會效益目標函數計算公式如下式(1)。
(2)經濟效益目標函數。以某區域用水戶能夠產生的最大經濟效益來表示,bkj和ckj為第k(1≤k≤3)個子區的第j(1≤j≤5)個用水戶的供水產值系數和供水費用系數,第i(1≤i≤3)個水源的供水次序系數用來表示,第k(1≤k≤3)個子區的第j(1≤j≤5)個用水戶的用水公平系數為。經濟效益目標函數如式(2)。
(1)供水量約束。不同水源對不同子區的不同用戶的供水量不超過其最大供水能力,供水量約束如式(3)、(4)和(5)所示。
式中:W1、W2、W3是三類水源分別對五類用水戶的可供水量。
(2)需水量約束。區域內的水源要滿足各用戶的需水量要求,但不能供水過多導致浪費,所以各水源供給量既要大于或者等于各用戶對各水源的最小需要量,也要小于或者等于各用戶對各水源的最大需水量,需水量約束如式(6)。
(3)非負約束。所有變量都應滿足大于等于零的要求。在該模型中均為已知量。
(1)供水效益系數bkj。通過分析DB14/T 1049.1-2015《山西省用水定額》和產值較高的一些用水部門的用水量,得出了晉中市不同行業的生產總值和用水量的比例系數,以確定生活用水、農業用水、第二用水、第三產業用水和生態環境用水的供水效益系數。供水效益系數取值如下:生活用水為600 元/m3,第二產業用水為480 元/m3,第三產業用水為450 元/m3,生態環境用水 300元/m3,農業用水為15 元/m3。
(2)供水費用系數ckj。通過分析晉中市2016年水費征收標準和相關水價政策,以確定生活用水、農業用水、第二用水、第三產業用水和生態環境用水的供水費用系數。供水費用系數取值如下:生活用水3.90 元/m3,農業用水0.25 元/m3,第二產業用水4.58 元/m3,第三產業用水5.76 元/m3,生態環境用水2.59元/m3。
(3)供水次序系數。由于晉中南部的第二產業較多,為了防止煤炭開采和金屬冶煉等產業對地下水的過度使用和污染,在水資源分配中應減少對地下水的使用,增加獲取和保護較容易的地表水的使用,因此本模型的供水次序為地表水、引調水、地下水。采用式(7)計算各類水源的供水次序系數,得到結果為0.50,0.33,0.17。
式中:子區k中水源i的供水次序號用表示表示其最大值。
(4)供水公平系數。本模型通過分析晉中南部各類用水戶的需求和影響,得到各用水戶獲水次序為生活用水、生態用水、農業用水、第二產業用水、第三產業用水,根據公式(7),得到生活、生態、農業、第二產業、第三產業對應的供水公平系數為0.33,0.27,0.20,0.13,0.07。
傳統NSGA-II算法依據個體間的支配關系將種群分為若干層,將進化群體的非支配個體集合設置為第一層,將去掉第一層個體后的種群中所求得的非支配個體集合設置為第二層,將去掉第一層和第二層個體后的種群中所求得的非支配個體集合設置為第三層,以此類推,然后計算每個個體的擁擠距離。NSGA-II算法的選擇操作首先對個體所在層的次序和擁擠距離進行比較,從中選擇出最優個體形成新的父代種群;接著使用遺傳算法產生新的子代種群,最后將新的父代種群與子代種群合并,進行基本的遺傳操作,反復執行,在滿足新進化種群的大小要求后停止操作[11]。
對NSGA-II算法進行了改進,并把改進后的NSGA-II算法稱為INSGA-II算法。改進思路是:通過改進擁擠度的計算方式和使用動態擁擠度來解決擁擠度相同時應該刪除哪個個體的問題,同時在刪除擁擠度最低的解后更新相鄰個體的擁擠度,以此來保證種群的多樣性。提出的INSGA-II算法從以下兩個方面對NSGA-II算法進行了改進。
(1)左右擁擠度。在NSGA-II算法中,若兩個個體的擁擠距離相等,在依據擁擠距離刪除個體時無法判斷該刪除哪一個個體。舉例,如表3。

表3 擁擠距離表Tab.3 Congestion distance table
表3為6個個體在單目標函數下的擁擠度值,按照NSGA-II算法的規則,根據公式(8)可得個體3的擁擠距離為P[4]-P[2],個體4的擁擠距離為P[5]-P[3],這兩個個體的擁擠距離相等。
式中:P[i]distance為個體i的擁擠距離;P[i]fk為個體i在子目標上fk上的函數值。
在依據擁擠距離刪除個體時就無法判斷該刪除哪一個個體,NSGA-II算法將會隨機選擇一個個體刪除[14]。為了避免這種問題的發生,設計的INSGA-II算法采用了一種新的左右擁擠度指標:首先初始化每個個體的擁擠度,給邊界節點的擁擠度賦予最大值保證其能入選下一代,對非邊界節點給其添加左擁擠度P[i]distance-front和右擁擠度P[i]distance-after并初始化為0,在對子目標m的函數值排序后,依次計算非邊界節點的左擁擠度和右擁擠度,最后根據左擁擠度和右擁擠度計算節點的擁擠度。具體計算方法如下:

(2)動態擁擠度排序策略。原始NSGA-II算法通過以個體的Pareto分層數和個體擁擠度為依據來選擇優秀個體,在同一Pareto分層中,擁擠度是選擇和淘汰個體的關鍵因素。但實驗分析,通過使用固定擁擠度所得到的個體排序方法存在不足,當某個區域聚集了擁擠度低的個體時,固定擁擠的排序策略可能會將該區域的全部個體被淘汰,使最終得到的Pareto解集的多樣性變差[15]。
在實際應用中,當某個個體被淘汰后,與之相鄰的個體擁擠度將發生變化,而用固定擁擠度排序法將會忽略這一點變化,所以會導致解集多樣性變差。
為了克服固定擁擠度排序法在實際應用中的不足,設計的INSGA-II算法使用一種動態擁擠度排序策略:當算法需要從種群第Fn層的h個個體中淘汰擁擠度最小的個體時,使用動態擁擠度排序策略,在淘汰之后記錄該解的位置i,并且重新計算位于i+1和i-1的解的擁擠度,根據更新后的擁擠度排序來刪除擁擠度最小的解,直到解集中解的數量滿足要求為止。
以不同擁擠度策略在標準測試函數MOP2上的表現為例,其中種群個數為8,使用固定擁擠度排序策略的解的分布情況如圖2,使用動態擁擠度排序策略解的分布情況如圖3。

圖2 使用固定擁擠度排序策略的解分布情況Fig.2 Distribution of solutions using fixed congestion sorting strategy

圖3 使用動態擁擠度排序策略的解分布情況Fig.3 Distribution of solutions using state congestion ranking strategy
對比圖2與圖3可以看出:采用固定擁擠度排序策略進行非支配解篩選得到的解疏密不一,采用動態擁擠度排序策略進行非支配解篩選得到的解分布性更好,證明通過使用動態擁擠度的排序策略得到的種群的多樣性更優。
INSGA-II算法的流程如下:
步驟1:初始化種群。根據約束條件隨機產生規模為N的初始種群Pt(t=0),Qt=?。對每個個體作函數評價,并依次對種群中的個體進行快速非支配排序、動態擁擠度計算。
步驟2:選擇父代種群。對P0中的個體所在層的次序和擁擠距離進行比較,從中選擇出最優個體形成新的父代種群。
步驟3:生成子代種群。對父代進行交叉和變異,產生新的子代種群Qt。
步驟4:混合父代和子代種群,形成新種群Rt=Pt∪Qt。
步驟5:選擇。依據中間種群快速非支配排序的分層數和本文使用的動態擁擠度排序的方式選擇優秀個體加入種群Pt+1。
步驟6:終止條件。令t=t+1,若t大于最大迭代次數MaxGen,則算法終止,同時令Pt中非支配解集作為Pareto最優解集; 若t小于最大迭代次數MaxGen,則返回步驟2。
采用Zitzler 和 Deb[16]所列舉的6 個測試函數ZDT1、ZDT2 、ZDT3 、ZDT4、ZDT5和ZDT6來驗證INSGA-II算法的性能。這6個測試函數見表4,它們的Pareto前沿包含了連續、非連續、凸函數、凹函數等不同情況,這些測試函數可以科學有效地反映算法的性能,得到的測試結果也更加全面。

表4 測試函數Tab.4 Test functions
多目標優化問題通過使用收斂帕累托最優解集來獲得解,并需要解保持多樣性。選取兩個性能指標來評價多目標優化算法得到的帕累托最優解集。
(1)反世代距離評價指標(Inverted Generational Distance,IGD)。IGD主要通過將Pareto前沿上的個體到算法獲取的Pareto最優解集之間的最小距離進行求和,來對算法的綜合性能包括收斂性和分布性進行評價,其計算公式如式(9)。IGD越小,表示算法的收斂性和分布性越好。
式中:P為分布在Pareto面上的個體;|P|為分布在Pareto 面上的個體的數量;Q為算法獲取的Pareto 最優解集;d(v,Q)為P中個體v到種群Q的最小歐幾里得距離。
(2)間距指標(Spacing Metric) 。Spacing的計算公式如式(10)。
式中:dp表示兩個連續帕累托前沿之間的歐幾里得距離;表示dp的平均值;表示算法所得Pareto前沿與理想Pareto前沿之間的歐幾里得距離;Spacing值的大小表示Pareto前沿分布情況,Spacing值越小說明算法所得的Pareto前沿分布更加均勻。
采用的測試函數的種群大小設置為200,算法最多迭代1 000次,采用IGD和Spacing兩個指標來評價計算得到的帕累托前沿的優劣[17]。各個對比算法在測試函數上的IGD值如表5所示。

表5 不同優化算法的IGD值Tab.5 IGD values for different optimization algorithms
從表5可以看出:通過將INSGA-II算法和NSGA-II、NS‐GA-II(SBX)、MOPSO、SPEA2、dMOPSO算法進行比較[18],在測試函數ZDT1和ZDT2中,INSGA-II算法計算得到的IGD值比其他對比算法更小,在測試函數ZDT5中表現不如原始的NSGAII算法,對于其余測試函數,INSGA-II算法計算得到的IGD值處于次優。這說明本文所提出的INSGA-II算法在求解凸、非凸和離散的多目標優化問題時能夠找到真正的帕累托前沿,并且具有良好的收斂性和分布性。將提出的改進算法運用于多目標水資源優化配置模型的求解,有利于高效精準地獲得最優調配方案
各個對比算法在測試函數上的Spacing值如表6所示。

表6 不同優化算法的Spacing值Tab.6 SpacingValues of different optimization algorithms
從表6可以看出:通過將INSGA-II算法計算得到的Spacing值與 NSGA-II(SBX)、MOPSO、SPEA2、dMOPSO算法進行比較[18],在測試函數 ZDT1、ZDT2和 ZDT3中,INSGA-II算法計算得到的Spacing值比其他對比算法更小,在ZDT4、ZDT5和ZDT6中,INSGA-II算法計算得到的Spacing值均較優,這說明所提出的INSGA-II算法在求解凸、非凸和離散的多目標優化問題時所得到的帕累托前沿分布較均勻。在求解多目標水資源優化配置模型時,使用本文提出的改進算法有利于得到分布更均勻的優化調配方案。
綜上所述,INSGA-II算法在解決多目標優化問題時能夠有效得到帕累托前沿,其IGD和Spacing值要優于其他對比算法。
使用INSGA-II算法求解晉中市水資源優化配置模型,并通過MATLAB進行編程,算法參數設置如下:種群大小為50,最大迭代次數為400,交叉概率為0.7,個體間變異概率為0.4,共得到非支配解集的個數為50。迭代次數為400時,INSGA-II算法和NSGA-II算法非支配解集對比如圖4。因為提出模型中的子目標相互沖突,子目標同時達到最優值相當困難,只能通過合理分配使其向最優方向優化,即多目標問題得到的是一個非支配解集(其中一個解代表一種配置方案)。

圖4 INSGA-II算法和NSGA-II算法非支配解集對比Fig.4 Comparison of non dominated solution set between INSGA-II and NSGA-II
由圖4可知:在迭代次數為400時,基于INSGA-II算法的水資源優化配置方案得到的非支配解集更符合要求,即解集向缺水量減小和經濟效益增大的方向收斂,其中部分解集得到的經濟效益更高、缺水量更小,證明使用的INSGA-II算法在求解晉中市水資源優化模型時有更好的收斂性和分布性。
在確定社會效益和經濟效益兩目標時,采用主觀賦權法,構造公式如式(11)。
式中:x表示不同的水資源配置構成的決策變量;h表示目標個數;f1(x)表示社會效益函數,f2(x)表示經濟效益函數;ωh表示目標對應的權重系數,權重之和為1。
根據對兩目標需求的不同賦予不同的權重,最終得到3種方案如下。
方案1:優先最大程度滿足供水,再考慮經濟效益,得到該方案的多目標綜合權重Ω=(ω1,ω2) =(0.8,0.2)。水資源分配方案見表7。

表7 方案1的水資源分配方案Tab.7 Water resources allocation scheme of scheme 1
方案2:同時考慮經濟效益和社會效益,得到該方案的多目標綜合權重Ω=(ω1,ω2)=(0.5,0.5)。水資源分配方案見表8。

表8 方案2的水資源分配方案Tab.8 Water resources allocation scheme of scheme 2
方案3:優先考慮經濟效益,再滿足社會效益,得到該方案的多目標綜合權重Ω=(ω1,ω2) =(0.2,0.8)。水資源分配方案見表9。

表9 方案3的水資源分配方案Tab.9 Water resources allocation scheme of scheme 3
3種方案的社會效益(缺水量)和經濟效益比較見表10。

表10 3種水資源配置方案的效益比較Tab.10 Benefit comparison of three water resources allocation schemes
從表10可看出:方案1是偏好社會效益的方案,其缺水量最低(社會效益最高);方案2是均衡方案,該方案同時兼顧了社會效益和經濟效益,其社會效益和經濟效益值都居中;方案3是偏好經濟效益的方案,其經濟效益值最高。決策者可以在實際情況中根據不同的情況選擇不同的方案,為晉中市水資源優化配置提供了全新的思路。
提出了一種使用新型擁擠度的計算策略和使用動態擁擠度的INSGA-II算法。同時分析晉中南部水資源狀況,綜合考慮社會效益和經濟效益,構建了晉中南部水資源優化配置模型,并將INSGA-II算法應用于晉中南部水資源優化配置模型中,得到了缺水量最小和經濟效益最大的非支配解集。該解集可以為晉中南部地區水資源的優化配置提供多種可選的解決方案,同時可為晉中市南部地區水資源的可持續發展提供決策依據。該配置結果同時考慮了社會效益和產生的經濟效益,在解決晉中市南部水資源優化配置模型的同時,為其他地區的水資源優化配置模型的求解提供了新的思路。本文未收集污染數據,同時在經濟效益函數中確定供水效益和費用系數時沒有考慮不同用水戶的分攤情況,待相關數據收集完整后應在水資源優化配置模型中加入生態環境效益函數,并且在經濟效益函數中考慮分攤情況。通過INSGA-II算法將水質和水量分配相聯系、以及在更大區域內多約束條件下的水資源優化配置問題,同時對以社會、經濟及生態環境效益為目標的多目標水資源優化配置模型進行求解。