閆志宏，王樹謙，劉 彬，徐 丹，李 蘇

(河北工程大學(xué)水利水電學(xué)院，河北 邯鄲 056038)

0 引 言

水資源優(yōu)化配置一般是多目標(biāo)優(yōu)化問題，各個(gè)目標(biāo)一般不可比較，改善其中某個(gè)目標(biāo)往往會造成其他目標(biāo)變劣[1]，不可能使所有目標(biāo)同時(shí)達(dá)到最優(yōu)。最初，多目標(biāo)優(yōu)化問題一般通過加權(quán)法、目標(biāo)規(guī)劃法和約束法等方式將其轉(zhuǎn)化為單目標(biāo)問題，然后利用現(xiàn)有較為成熟的單目標(biāo)算法求解，每次運(yùn)算往往僅能得到一組局部最優(yōu)解[2]。自Rosenberg[3]于1967年提出利用基于進(jìn)化的搜索方法來求解多目標(biāo)優(yōu)化問題以來，進(jìn)化多目標(biāo)優(yōu)化算法(EMO)引起了很多學(xué)者的關(guān)注[4-10]。近年來，一些基于自然元啟發(fā)式優(yōu)化算法(Nature-inspired meta-heuristic algorithms)相繼被提出用于求解多目標(biāo)優(yōu)化問題[11-15]，并已經(jīng)成功應(yīng)用到水資源優(yōu)化配置問題上[16-19]。

2015年，Mirjalili[20]提出了飛蛾火焰優(yōu)化算法(Moth-Flame Optimization，簡稱MFO)，其求解思路來源于飛蛾橫向定位導(dǎo)航機(jī)制，是一種新型的群智能算法。MFO已被成功應(yīng)用于電力系統(tǒng)[21-23]、圖像分割[24,25]和網(wǎng)絡(luò)入侵檢測[26]等工程實(shí)踐中。MFO算法具有局部搜索能力強(qiáng)的特點(diǎn)，但是該算法全局收斂能力較弱，運(yùn)行時(shí)容易收斂到局部最優(yōu)。針對MFO算法存在的不足，本文提出一種改進(jìn)MFO算法，為人工飛蛾的捕焰行為引入帶精英策略的快速非支配排序方法、擁擠度和擁擠度比較算子，并通過仿真實(shí)驗(yàn)，與VEGA、NSGA-II、MODE、BEES和SPEA等算法對比來檢驗(yàn)改進(jìn)的飛蛾火焰算法性能，仿真實(shí)驗(yàn)結(jié)果證明改進(jìn)MFO算法在求解多目標(biāo)優(yōu)化問題時(shí)，能較好的收斂到全局最優(yōu)解，求解精度較高。然后將其應(yīng)用到水資源優(yōu)化配置模型中，旨在為求解水資源優(yōu)化配置問題提供新的求解方法。

1 改進(jìn)飛蛾火焰算法

1.1 飛蛾火焰算法

MFO算法[20]是模仿飛蛾種群追逐火焰做螺旋運(yùn)動的行為，之后將該種行為模型化的新型群智能算法。

1.1.1 種群初始化

在飛蛾火焰算法中，假設(shè)優(yōu)化問題的解是飛蛾，并且用于描述優(yōu)化問題的變量為飛蛾所處空間位置。其數(shù)學(xué)模型描述如下：飛蛾種群規(guī)模為n，待尋優(yōu)變量個(gè)數(shù)為d，矩陣M存儲飛蛾所處的空間位置，OM存儲飛蛾個(gè)體的適應(yīng)度值。例如，矩陣 的第一行每個(gè)飛蛾個(gè)體輸出的目標(biāo)函數(shù)值存儲在矩陣OM的OM1中。

(1)

(2)

MFO算法中另一關(guān)鍵部分是火焰。設(shè)火焰矩陣規(guī)模F同樣為n*d，如式(3)所示，飛蛾矩陣M根據(jù)其適應(yīng)度值排序所得的矩陣存儲到火焰矩陣F中。并利用式(4) 矩陣存儲火焰的適應(yīng)度值。

(3)

(4)

1.1.2 位置更新機(jī)制

人工飛蛾圍繞火焰螺旋運(yùn)動的行為可以概化為捕焰及棄焰過程。

(1)捕焰。人工飛蛾Mi會朝著離自身距離最近的火焰Fi做對數(shù)螺旋運(yùn)動來捕獲火焰。其計(jì)算公式為：

S(Mi,Fj)=Di·ebt·cos(2πt)+Fj

(5)

式中：S(Mi,Fj)表示飛蛾Mi繞火焰Fj運(yùn)動更新后的位置；第i只飛蛾用Mi表示；第j個(gè)火焰用Fj表示；Mi與Fj之間的距離用Di表示：Di=|Fj-Mi|；b為常數(shù)；t為隨機(jī)數(shù)，大小為[-1，1]，用來表示飛蛾下一個(gè)位置距離火焰遠(yuǎn)近程度，當(dāng)參數(shù)t=-1時(shí)，該位置距離火焰最近，t=1時(shí)則表示最遠(yuǎn)。

(2)棄焰。在迭代計(jì)算時(shí)，飛蛾會不斷舍棄適應(yīng)度值低的火焰，朝向較優(yōu)的火焰做螺旋運(yùn)動，從而自適應(yīng)減少火焰的數(shù)量，使火焰數(shù)量越來越少最終收斂于1。其計(jì)算公式為：

(6)

式中：l為當(dāng)前迭代次數(shù)；N為最大火焰數(shù)量；T為最大迭代次數(shù)。

1.2 改進(jìn)飛蛾火焰算法

本文在采用飛蛾火焰算法處理多目標(biāo)優(yōu)化問題時(shí)，引入Deb[27]等人在帶精英策略的快速非支配排序遺傳算法(A Fast Elitist Non-dominated Sorting Genetic Algorithm, NSGA-Ⅱ)中提出的帶精英策略的快速非支配排序方法、擁擠度和擁擠度比較算子?？焖俜侵渑判蚍椒ê蛽頂D度及擁擠度比較算子在算法運(yùn)行時(shí)能夠在Pareto前沿面均勻地選擇個(gè)體，防止局部收斂。

(1)快速非支配排序策略?？焖俜侵渑判虿呗允且罁?jù)支配關(guān)系把所有個(gè)體劃分為不同的等級，旨在從種群中選擇出相對優(yōu)秀的個(gè)體。假設(shè)某飛蛾種群的大小為M，對于飛蛾種群的某個(gè)個(gè)體i，需計(jì)算其ni和Si值。其中，ni表示能支配i的個(gè)體數(shù)，Si表示能夠被i所支配的個(gè)體數(shù)。將所有ni=0的個(gè)體存儲至當(dāng)前集合H1中；被H1中的個(gè)體k支配的所有個(gè)體存放到集合Fk中，將Fk中的每個(gè)個(gè)體進(jìn)行nk-1操作，將nk-1=0的個(gè)體存放到集合S中；集合H1中的個(gè)體記為第一級，然后將集合S記為當(dāng)前集合，重復(fù)以上過程，直至種群中所有個(gè)體均被分級[27]。

(2)擁擠度和擁擠度比較算子。

①擁擠度。擁擠度[27]是指種群中圍繞個(gè)體i的其余個(gè)體的密度，其計(jì)算公式為：

(7)

因?yàn)楦髂繕?biāo)函數(shù)之間的單位不一樣，需要對其進(jìn)行歸一化操作，式(7)可改寫為[28]：

(8)

式中：n為目標(biāo)函數(shù)的個(gè)數(shù)(n=1,2,…,N)；fn(i)為個(gè)體i在第n個(gè)目標(biāo)函數(shù)上的值；fnmax為所有個(gè)體在第n個(gè)目標(biāo)函數(shù)所得到的最大值，fnmin則為其最小值。

②擁擠度比較算子。經(jīng)過快速非支配排序和擁擠度計(jì)算，每個(gè)個(gè)體都會得到irank和id兩個(gè)參數(shù)，其中irank為非支配排序，id為擁擠度。irank和id用以區(qū)分任意兩個(gè)個(gè)體的支配關(guān)系。

i≥nj(ifirankjd)

(9)

式中：≥n表示擁擠度比較算子。

1.3 運(yùn)算流程

改進(jìn)飛蛾火焰算法運(yùn)算流程如下：

(1)初始化，隨機(jī)產(chǎn)生飛蛾種群Ml，最大迭代次數(shù)為T，飛蛾種群大小為n，最大火焰數(shù)量為N，維數(shù)為d，當(dāng)前迭代次數(shù)用l表示。

(2)當(dāng)種群所有個(gè)體執(zhí)行快速非支配排序操作被分級且每個(gè)個(gè)體擁擠度計(jì)算完成后，依據(jù)式(9)判斷所有個(gè)體的支配關(guān)系，將支配關(guān)系結(jié)果保存為火焰矩陣。

(3)采用式(6)對火焰數(shù)量進(jìn)行更新；計(jì)算飛蛾與火焰之間的距離，并利用式(5)對飛蛾-火焰位置進(jìn)行更新，得到當(dāng)前迭代飛蛾種群Ml+1。

(4)將種群M1和種群Ml+1進(jìn)行合并(Rl+1=Ml∪Ml+1)，對Rl+1執(zhí)行快速非支配排序操作并計(jì)算其擁擠度，將前n個(gè)個(gè)體組成新種群Pl+1，保存為火焰矩陣，利用式(2)存儲飛蛾位置，式(4)存儲火焰空間位置。

(5)找出最優(yōu)飛蛾個(gè)體的當(dāng)前位置。如該位置優(yōu)于先前所保留的位置，則將此位置保存為最佳空間位置。如果滿足算法迭代停止條件執(zhí)行步驟(6)，否則轉(zhuǎn)至步驟(3)～(5)。

(6)輸出結(jié)果，算法結(jié)束。

2 仿真實(shí)驗(yàn)

本文采用3個(gè)多目標(biāo)測試函數(shù)來驗(yàn)證改進(jìn)飛蛾火焰算法的性能。這三個(gè)測試函數(shù)是Zitzler和Deb[29]所列舉的3個(gè)典型的測試函數(shù)ZDT1、ZDT2和ZDT3，每個(gè)測試函數(shù)都有兩個(gè)目標(biāo)函數(shù)。ZDT1、ZDT2和ZDT3的帕累托最優(yōu)前沿分別是為凸的、非凸的和不連續(xù)的。與單目標(biāo)優(yōu)化不同的是，多目標(biāo)優(yōu)化問題的解集收斂到帕累托最優(yōu)解集，并需保持多樣性，這兩個(gè)目標(biāo)不能用一個(gè)性能指標(biāo)來評價(jià)，不同文獻(xiàn)[30-32]提出了一系列評價(jià)指標(biāo)，本文選取兩個(gè)性能指標(biāo)來評價(jià)采用多目標(biāo)優(yōu)化算法計(jì)算得到的帕累托最優(yōu)解集。

(1)收斂性指標(biāo)：Convergence Distance(CD)。CD表示計(jì)算得到的帕累托最優(yōu)前沿的某個(gè)解與其距離最近的理想帕累托最優(yōu)前沿的某個(gè)解之間的最小歐幾里得距離的平均值。CD是收斂性度量指標(biāo)，其值越小表示計(jì)算得到的帕累托最優(yōu)前沿越接近于理想帕累托最優(yōu)前沿。如果采用某個(gè)多目標(biāo)優(yōu)化算法計(jì)算得到的帕累托最優(yōu)前沿全部收斂于理想帕累托最優(yōu)前沿，那么其CD值等于0。計(jì)算公式如下：

(10)

(2)間距指標(biāo)：Spacing Metric(Δ)。Δ的計(jì)算公式如下：

(11)

本文采用的多目標(biāo)優(yōu)化算法種群規(guī)模均設(shè)置為1 000，迭代次數(shù)為1 000次。測試函數(shù)ZDT1、ZDT2和ZDT3計(jì)算得到的帕累托前沿和理想帕累托前沿見圖1-圖3所示。在圖中可以看到對于所有的測試函數(shù)，多目標(biāo)飛蛾火焰算法計(jì)算得到的帕累托前沿近似于理想帕累托前沿。采用CD和Δ兩個(gè)指標(biāo)來評價(jià)計(jì)算得到的帕累托前沿的優(yōu)劣。ZDT1～ZDT3的CD值見表1所示，ZDT1～ZDT3的Δ值見表2所示。

圖1 ZTD1計(jì)算得到與理想帕累托前沿關(guān)系圖Fig.1 The relationship of the obtained Pareto front and the true Pareto front for ZTD1

圖2 ZTD2計(jì)算得到與理想帕累托前沿關(guān)系圖Fig.2 The relationship of the obtained Pareto front and the true Pareto front for ZTD2

圖3 ZTD3計(jì)算得到與理想帕累托前沿關(guān)系圖Fig.3 The relationship of the obtained Pareto front and the true Pareto front for ZTD3

多目標(biāo)飛蛾火焰算法計(jì)算得到的CD值與VEGA、NSGA-Ⅱ、

表1 多目標(biāo)優(yōu)化算法計(jì)算得到的CD值Tab.1 The CD value calculated by multi-objective optimization algorithm

表2 多目標(biāo)優(yōu)化算法計(jì)算得到的Δ值Tab.2 The Δ value calculated by multi-objective optimization algorithm

MODE、BEES和SPEA等其他多目標(biāo)算法進(jìn)行比較[32]。從表1中，對于測試函數(shù)ZDT1～ZDT2，多目標(biāo)飛蛾算法計(jì)算得到的CD均優(yōu)于其他多目標(biāo)算法；對于ZDT3，除DEMO外，多目標(biāo)飛蛾優(yōu)化算法計(jì)算得到的CD均優(yōu)于其他多目標(biāo)算法。這說明多目標(biāo)飛蛾火焰算法求解凸、非凸和離散的多目標(biāo)優(yōu)化問題時(shí)能夠找到真正的帕累托前沿。

采用多目標(biāo)飛蛾火焰算法計(jì)算得到的Δ值與NSGA-Ⅱ(RC)、NSGA-Ⅱ(BC)、SPEA、PAES、PAES-Ⅱ、σ-MOPSO、NSPSO和MOPSO等其他多目標(biāo)算法進(jìn)行了比較[28]。從表2可知，對于ZDT1而言，多目標(biāo)飛蛾火焰算法的Δ僅次于NSGA-Ⅱ(RC)和σ-MOPSO；對于ZDT2而言，多目標(biāo)飛蛾火焰算法的Δ僅次于NSGA-Ⅱ(RC)、NSGA-Ⅱ(BC)和σ-MOPSO；對于ZDT3而言，多目標(biāo)飛蛾火焰算法的Δ僅次于NSGA-Ⅱ(BC)、SPEA和NSPSO。

綜上，多目標(biāo)飛蛾火焰算法在解決多目標(biāo)優(yōu)化問題時(shí)能夠得到真正的帕累托前沿，其CD和Δ要優(yōu)于本文所列的大部分多目標(biāo)算法，能將改進(jìn)的飛蛾火焰算法應(yīng)用于工程實(shí)踐中。

3 案例分析

3.1 區(qū)域概況

三亞市是海南省南部的政治、經(jīng)濟(jì)、文化中心，位于東經(jīng)108°56′30″～109°48′28″，北緯18°09′34″～18°37′27″，東鄰陵水，北依保亭，西毗樂東，南臨南海及三沙市，東西長91.6 km，南北寬51 km。行政分區(qū)主要包括城區(qū)、海棠灣鎮(zhèn)、吉陽鎮(zhèn)、鳳凰鎮(zhèn)、天涯鎮(zhèn)、崖城鎮(zhèn)和育才鎮(zhèn)。

3.2 三亞市水資源優(yōu)化配置模型構(gòu)建

3.2.1 水資源優(yōu)化配置數(shù)學(xué)模型

多目標(biāo)水資源優(yōu)化配置的目標(biāo)是達(dá)到社會、經(jīng)濟(jì)及生態(tài)環(huán)境綜合效益最大。對水資源進(jìn)行優(yōu)化配置的最終目的是高效利用水資源，促進(jìn)水資源與經(jīng)濟(jì)社會的協(xié)調(diào)可持續(xù)發(fā)展。本文選取社會及經(jīng)濟(jì)效益作為目標(biāo)進(jìn)行求解。

(12)

(13)

(3)約束條件。

①可供水量約束：

(14)

②輸水能力約束：

(15)

③需水量約束。

(16)

3.2.2 模型參數(shù)確定

三亞市水資源優(yōu)化配置模型概化為7個(gè)子區(qū)，3種水源，6類用水戶。以2030年為規(guī)劃水平年，在P=50%保證率下進(jìn)行水資源的優(yōu)化配置，以達(dá)到高效利用水資源，社會效益和經(jīng)濟(jì)效益最佳的目標(biāo)。根據(jù)三亞市的發(fā)展規(guī)劃，需水量采用定額法對6類用水戶進(jìn)行需水預(yù)測，各分區(qū)各用水戶需水量預(yù)測[34]見表3所示。在對三亞市水資源進(jìn)行評價(jià)的基礎(chǔ)上，對地表水資源、地下水資源和再生水資源進(jìn)行了預(yù)測，預(yù)測結(jié)果見表4所示[34]。

表3 三亞市2030年(P=50%)需水預(yù)測結(jié)果 萬m3

表4 三亞市2030年(P=50%)可供水量預(yù)測結(jié)果 萬m3

注：√表示該種水源可為該子區(qū)供水。

(17)

表5 供水次序系數(shù)表Tab.5 The Water supply order coefficient table

(4)采用式(17)計(jì)算得到6類用水戶的用水公平系數(shù)如下：城鎮(zhèn)生活用水戶0.29，農(nóng)村生活用水戶為0.24，生態(tài)環(huán)境用水戶0.19，第二產(chǎn)業(yè)用水戶0.14，第三產(chǎn)業(yè)用水戶0.10，第一產(chǎn)業(yè)用水戶0.05。

3.3 結(jié) 果

采用多目標(biāo)飛蛾火焰算法求解三亞市水資源優(yōu)化配置模型。采用MATLAB編寫程序進(jìn)行計(jì)算，參數(shù)設(shè)置如下：種群大小500，火焰最大數(shù)目500，迭代次數(shù)1 000 次。得到三亞市水資源優(yōu)化配置的目標(biāo)函數(shù)值的帕累托前沿解集見表6，共求解得到23組解。決策者可以根據(jù)實(shí)際需求選擇與之相適應(yīng)的方案：如對社會效益有特殊偏好可以選擇表6中標(biāo)志為a的方案；如對經(jīng)濟(jì)效益有特殊偏好可以選擇表6中標(biāo)志為b的方案；如對兩個(gè)目標(biāo)函數(shù)沒有特殊偏好，可以選擇余下的方案。三亞市屬于旅游城市，第三產(chǎn)業(yè)用水量較大，為了提高水資源利用效率，本文擬選擇對社會效益有特殊偏好的方案8進(jìn)行詳細(xì)分析。

表6 改進(jìn)飛蛾火焰算法求解帕累托前沿解集結(jié)果Tab.6 Results of Pareto front under improved moth flame algorithm

飛蛾火焰算法種群個(gè)體平均值迭代過程如圖4。由圖4可知，社會效益目標(biāo)函數(shù)種群個(gè)體平均值迭代曲線在迭代800次左右時(shí)趨于收斂，800次以后雖有小范圍波動但以穩(wěn)定的趨勢運(yùn)行；經(jīng)濟(jì)效益目標(biāo)函數(shù)種群個(gè)體平均值迭代曲線在迭代600次左右時(shí)趨于收斂，600次以后雖有小范圍波動但以穩(wěn)定的趨勢運(yùn)行。

圖4 改進(jìn)飛蛾火焰算法種群個(gè)體平均值迭代過程圖Fig.4 Iterative process of the mean value under improved moth flame algorithm

三亞市各子區(qū)水量分配結(jié)果詳見圖5所示。從圖5可知，三亞市各用水戶總需水量為39 015 萬m3，各用水戶總分配水量為39 015 萬m3，缺水量為0。地表水、地下水和再生水分配給三亞市各用水戶的水量分別為：34 779、954和3 282 萬m3，占總分配水量的比例分別為：89.14%、2.45%和8.41%，可知地表水是三亞市的主要水源。

圖5 三亞市水資源優(yōu)化配置結(jié)果Fig.5 The water allocation results of Sanya City

三亞市各子區(qū)各用水戶分配水量結(jié)果見表7示。由表7可知，三亞市各子區(qū)各用水戶總分配水量為39 015 萬m3，6類用水戶分配水量分別為：6 456、1 402、20 517、2 397、7 525和718 萬m3，其占總分配水量的比例分別為：16.55%、3.59%、52.59%、6.14%、19.29%和1.84%，由此可知第一產(chǎn)業(yè)是用水大戶，所占比例最大。

4 結(jié) 語

飛蛾火焰算法在處理復(fù)雜優(yōu)化問題上收斂性較差和易陷入局部最優(yōu)，本文結(jié)合帕累托最優(yōu)策略，快速非支配排序策略，擁擠度及擁擠度比較算子，將種群所有個(gè)體進(jìn)行排序分級，優(yōu)選出最優(yōu)個(gè)體。通過仿真實(shí)驗(yàn)，得出改進(jìn)MFO算法在求解復(fù)雜多目標(biāo)優(yōu)化問題時(shí)其收斂性能和求解精度優(yōu)于本文所列大部分多目標(biāo)優(yōu)化算法，繼而將改進(jìn)MFO算法應(yīng)用到多目標(biāo)水資源優(yōu)化配置模型中，并得到了成功應(yīng)用。

將改進(jìn)MFO算法應(yīng)用到三亞市水資源優(yōu)化配置模型中，得到了該多目標(biāo)優(yōu)化問題的帕累托最優(yōu)前沿，共包含23組帕累托解。本文選擇對社會效益有特殊偏好的方案作為最終決策方案，結(jié)果顯示三亞市不同用水戶總需水量39 015 萬m3，各用水戶總分配水量為39 015 萬m3，缺水量為0，產(chǎn)生的經(jīng)濟(jì)效益為203.91 億元。該配置結(jié)果傾向于達(dá)到社會效益最佳的目標(biāo)，兼顧產(chǎn)生的經(jīng)濟(jì)效益，符合三亞市水資源開發(fā)利用原則，為求解多目標(biāo)水資源優(yōu)化配置問題提供了新的方法。

表7 三亞市水資源優(yōu)化配置結(jié)果 萬m3

本文在確定供水效益和供水費(fèi)用系數(shù)時(shí)未考慮分?jǐn)傁禂?shù)，有待進(jìn)一步研究。由于環(huán)境污染數(shù)據(jù)難以收集，待數(shù)據(jù)收集完整后應(yīng)采用改進(jìn)飛蛾火焰算法對包含社會、經(jīng)濟(jì)及生態(tài)環(huán)境效益的多目標(biāo)水資源優(yōu)化配置模型進(jìn)行求解。