改進飛蛾火焰算法在多目標水資源優(yōu)化配置中的應用

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

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

0 引 言

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

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

1 改進飛蛾火焰算法

1.1 飛蛾火焰算法

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

1.1.1 種群初始化

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

(1)

(2)

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

(3)

(4)

1.1.2 位置更新機制

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

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

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

(5)

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

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

(6)

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

1.2 改進飛蛾火焰算法

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

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

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

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

(7)

因為各目標函數(shù)之間的單位不一樣，需要對其進行歸一化操作，式(7)可改寫為[28]：

(8)

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

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

i≥nj(ifirankjd)

(9)

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

1.3 運算流程

改進飛蛾火焰算法運算流程如下：

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

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

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

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

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

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

2 仿真實驗

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

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

(10)

(2)間距指標：Spacing Metric(Δ)。Δ的計算公式如下：

(11)

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

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

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

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

多目標飛蛾火焰算法計算得到的CD值與VEGA、NSGA-Ⅱ、

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

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

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

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

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

3 案例分析

3.1 區(qū)域概況

三亞市是海南省南部的政治、經(jīng)濟、文化中心，位于東經(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)化配置模型構建

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

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

(12)

(13)

(3)約束條件。

①可供水量約束：

(14)

②輸水能力約束：

(15)

③需水量約束。

(16)

3.2.2 模型參數(shù)確定

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

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

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

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

(17)

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

(4)采用式(17)計算得到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 結 果

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

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

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

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

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

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

三亞市各子區(qū)各用水戶分配水量結果見表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 結 語

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

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

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

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