劉龍龍，李文竹，劉 心

(1.河北工程大學(xué)，河北 邯鄲 056038；2.邯鄲河務(wù)局，河北 邯鄲 056001)

我國水資源存在的時空分布不均、用水矛盾突出等現(xiàn)實問題，越來越成為新時期經(jīng)濟社會發(fā)展的制約性因素[1]。水資源日益短缺，對區(qū)域水資源進行合理規(guī)劃、利用和調(diào)度勢在必行。居民社區(qū)用水是水資源消耗的重要組成部分[2]?？茖W(xué)的管網(wǎng)調(diào)度能節(jié)省大量的供水消耗，全面提高管網(wǎng)的安全性和可靠性，但調(diào)度方案是否可行很大程度上取決于用水量的預(yù)測精度[3]。

對用水量預(yù)測方法的研究一直是國內(nèi)外學(xué)者研究的熱點，主要有以時間序列法為代表的傳統(tǒng)方法、以人工神經(jīng)網(wǎng)絡(luò)(ANN)為代表的人工智能方法和以支持向量機(SVM)為代表的機器學(xué)習(xí)算法[4]。然而這些方法都只能得到確定的點預(yù)測結(jié)果和未來某一時刻的預(yù)測結(jié)果，而居民實際用水量蘊含各種不確定性因素，點預(yù)測結(jié)果必然存在不同程度的誤差，進而使調(diào)度決策工作面臨一定程度的風(fēng)險。另外，如果預(yù)測的結(jié)果不是實時的，僅是未來某一時刻的值，那么調(diào)度人員無法根據(jù)準確的時用水量值進行實時調(diào)度。如果能夠給出實時預(yù)測結(jié)果和其變化區(qū)間，使調(diào)度決策人員隨時了解未來用水量的變化波動范圍，有利于其做出更合理的決策。

高斯過程回歸(GPR)是通過尋找訓(xùn)練數(shù)據(jù)之間的關(guān)系來進行系統(tǒng)辨識的一種非參數(shù)黑箱模型,是用于概率問題預(yù)測建模的一種強大工具[5]。與ANN和SVM方法相比，GPR不僅能夠預(yù)測未知量的期望值，還能給出其分布狀況,同時也能進行實時預(yù)測。此外，GPR模型參數(shù)較少、易實現(xiàn)和泛化能力強等優(yōu)點[6]，可以直接方便的用于區(qū)間預(yù)測和實時預(yù)測[7]。通常采用共軛梯度法求解最優(yōu)超參數(shù)，但存在優(yōu)化效果對初始值依賴性強和容易陷入局部最優(yōu)的缺點。

蜻蜓算法是美國學(xué)者MIRJALILI S在2015年提出的一種基于種群的啟發(fā)式智能優(yōu)化算法[8]。自然界中的蜻蜓生活方式包括靜態(tài)和動態(tài)群體行為，這和啟發(fā)式優(yōu)化算法中的兩個重要階段探索和開發(fā)非常類似[9]。蜻蜓算法通過模擬蜻蜓群體航行、捕食及躲避外敵等行為進行全局和局部搜索，尋找獵物的過程就是算法尋優(yōu)的過程[10,11]。本文提出用擅長全局搜索和對初始值沒有依賴的蜻蜓算法進行超參數(shù)尋優(yōu)，形成蜻蜓-高斯過程回歸耦合算法(DA-GPR)，對居民社區(qū)時用水量進行動態(tài)實時區(qū)間預(yù)測。

1 蜻蜓-高斯過程回歸耦合預(yù)測模型

本文通過高斯過程回歸與蜻蜓算法的有機結(jié)合，利用擅長全局搜索和對初始值沒有依賴的蜻蜓算法進行最優(yōu)超參數(shù)尋優(yōu)，并利用高斯過程回歸(GPR)尋找時用水量訓(xùn)練數(shù)據(jù)之間的關(guān)系來進行系統(tǒng)辨識,給出其分布狀況,對居民社區(qū)時用水量進行動態(tài)實時區(qū)間預(yù)測。

1.1 DA-GPR基本原理

從函數(shù)空間角度來看，高斯過程回歸其性質(zhì)由均值函數(shù)m(x)和協(xié)方差函數(shù)k(x,x′)確定[12]，定義為：

(1)

式中：函數(shù)f(x)的高斯過程(Gaussian process,GP)數(shù)學(xué)表達式為f(x)～GP[m(x),k(x,x′)],x,x′∈Rd是任意的隨機變量。

(2)

式中：K(X,X)=Kn=(kij)為n×n階對稱正定協(xié)方差矩陣。

預(yù)測值的后驗分布為：

(3)

(4)

(5)

本文的協(xié)方差函數(shù)采用平方指數(shù)協(xié)方差函數(shù)，為：

(6)

1.2 DA-GPR目標函數(shù)設(shè)計

蜻蜓-高斯過程回歸耦合時用水量預(yù)測模型的最終目標是對未來的時用水量進行較準確的動態(tài)實時區(qū)間預(yù)測，并具有較好的泛化能力。因此在使用蜻蜓算法優(yōu)化高斯過程回歸的參數(shù)時，選取能充分代表泛化能力的目標函數(shù)是非常重要的。常規(guī)的蜻蜓優(yōu)化算法將訓(xùn)練集上獲得最小平均相對誤差(MAPE)的參數(shù)作為最優(yōu)解。即，將訓(xùn)練集上獲得最小平均相對誤差(MAPE)作為目標函數(shù)，對蜻蜓個體進行選著，目標函數(shù)值越小，蜻蜓個體越優(yōu)，對應(yīng)的參數(shù)越優(yōu)，目標函數(shù)的選著至關(guān)重要。

居民社區(qū)每小時的用水量比較少，基本上在幾立方米到十幾立方米，僅僅采用平均相對誤差(MAPE)作為目標函數(shù)值進行參數(shù)尋優(yōu)的效果不理想。在此基礎(chǔ)上，本文對DA-GPR進行改進，將體現(xiàn)絕對誤差大小的平均絕對誤差(MAE)和均方誤差(MSE)引入到目標函數(shù)中，將平均相對誤差(MAPE)、平均絕對誤差(MAE)和均方誤差(MSE)進行線性組合，構(gòu)造出新的目標函數(shù)，如下式所示。

F(x)=3MAPE+2MAE+MSE

(7)

1.3 DA-GPR局部搜索公式設(shè)計

蜻蜓群體在一個小的區(qū)域內(nèi)來回飛翔尋找獵物以及躲避外敵。通過模擬蜻蜓群體之間尋找獵物以及躲避外敵的社會互動行為來實現(xiàn)優(yōu)化的過程，就是算法的局部搜索過程。

蜻蜓個體在群體運動中可以分為分離、對齊、內(nèi)聚、食物吸引和天敵排斥等5種行為模式。該行為的具體意義和數(shù)學(xué)表達方法如下[13]。

(1)分離：指蜻蜓與相鄰個體之間避免碰撞。該行為的計算公式如下:

(8)

(2)對齊：指相鄰個體之間傾向于保持相同的速度。該行為的計算公式如下:

(9)

(3)內(nèi)聚：指蜻蜓個體試圖向自己認為所屬的群體靠近。該行為的計算公式如下:

(10)

(4)食物吸引力：指食物對蜻蜓的吸引力。該行為的計算公式如下:

Fi=X+-X

(11)

(5)天敵排斥力：指蜻蜓對天敵的排斥力。該行為的計算公式如下:

Ei=X-+X

(12)

式中：Si為第i個個體的分離量；X為當前個體的位置；Xj為相鄰個體j的位置；N為相鄰個體的數(shù)量；Ai為第i個個體的對齊量；Vj為第j個鄰近個體的速度；Ci為第i個個體的內(nèi)聚量；Fi為第i個個體對獵物的吸引力；X+為獵物所在的位置；Ei為第i個個體需逃離天敵的距離；X-為天敵所在的位置。

1.4 DA-GPR全局搜索公式設(shè)計

在蜻蜓尋優(yōu)的過程中，大量的蜻蜓群體朝著同一個方向進行長距離遷移，這個群體活動就是對算法的全局搜索過程。

根據(jù)5種蜻蜓行為，下一代蜻蜓的步長計算如下：

ΔXt+1=(sSi+aAi+cCi+fFi+eEi)+ωΔXt

(13)

有鄰近蜻蜓時，下一代蜻蜓的位置為：

Xt+1=Xt+ΔXt+1

(14)

無鄰近蜻蜓時，設(shè)置為隨機游走(Le′xy飛行)行為，下一代蜻蜓的個體位置為:

Xt+1=Xt+Le′xy(d)×Xt

(15)

式中：t為當前迭代次數(shù)；i為第i個蜻蜓個體；Xi為當前第t代種群個體位置；ΔXt+1為下一代種群位置更新步長；Xt+1為下一代種群個體位置；s，a，c，f，e分別對應(yīng)于5種行為的權(quán)重；ω為慣性權(quán)重；d為個體位置向量的維數(shù)。

Le′vy函數(shù)計算公式如下：

(16)

(17)

式中：r1，r2為[0,1]隨機數(shù)；Γ(x)=(x-1)!；β是一個常數(shù)(這里取為0.5)。

蜻蜓算法在尋優(yōu)的過程中，對每個個體的鄰近個體數(shù)量的計算是非常重要的，因此這里假定一個鄰域半徑，該半徑隨迭代次數(shù)的增加而成比例地增長。同時，為了達到算法的局部搜索和全局搜索的平衡，權(quán)重(s，a，c，f，e和ω)在優(yōu)化過程中自適應(yīng)地調(diào)整。

1.5 DA-GPR算法步驟

流程圖如圖1所示。

圖1 DA-GPR模型區(qū)間預(yù)測流程圖Fig.1 DA-GPR model interval prediction flow chart

具體實現(xiàn)步驟如下：

Step1：讀取樣本數(shù)據(jù)，產(chǎn)生訓(xùn)練集和測試集。

Step2：初始參數(shù)設(shè)置：設(shè)置蜻蜓種群規(guī)模、最大迭代次數(shù)和參數(shù)l,σf,σn的取值范圍。

Step3：初始化蜻蜓：隨機初始化蜻蜓的位置H、步長ΔH、鄰域半徑R，慣性權(quán)重ω，分離權(quán)重s，對齊權(quán)重a，內(nèi)聚權(quán)重c，食物吸引權(quán)重f，天敵排斥權(quán)重e。

Step4：將蜻蜓個體位置H的信息依次賦值給l,σf,σn。其中位置H矩陣的第一行存放參數(shù)l的值，第二行存放參數(shù)σf的值，第三行存放參數(shù)σn的值，每個蜻蜓個體對應(yīng)一組參數(shù)值。

Step5：創(chuàng)建GPR模型，對訓(xùn)練樣本進行訓(xùn)練，求出每組參數(shù)對應(yīng)的目標函數(shù)值，并將其作為蜻蜓算法的適應(yīng)度函數(shù)值，判斷當前的適應(yīng)度函數(shù)值是否為最佳適應(yīng)度值，若是則將對應(yīng)的超參數(shù)保存為最優(yōu)超參數(shù)值,否則仍保存原適應(yīng)度值及其所對應(yīng)的參數(shù)值。

Step6：判斷是否達到最大迭代次數(shù)，如果達到最大迭代次數(shù)，則輸出最優(yōu)超參數(shù)并創(chuàng)建最優(yōu)的GPR模型；如果達不到最大迭代次數(shù)，則依次更新蜻蜓最優(yōu)個體和最差個體、更新鄰域半徑、更新個體位置，然后返回步驟Step4繼續(xù)迭代。

Step7：將測試集數(shù)據(jù)輸入創(chuàng)建的最優(yōu)GPR模型，輸出預(yù)測值的均值和方差。

Step8：根據(jù)給定置信度，得到時用水量實時區(qū)間預(yù)測結(jié)果。

2 仿真實驗

2.1 仿真實驗平臺

本次實驗的運行環(huán)境為：操作系統(tǒng)：Windows 7；處理器：Intel(R) Core(TM)i7-2640M；主頻：2.80GHz；內(nèi)存：4.00GB；編程工具：MATLAB R2015(b)。

2.2 樣本數(shù)據(jù)

本文所用的數(shù)據(jù)來自河北工程大學(xué)在線水量檢測平臺，采用居民家屬院2016年6月17日到7月7日，共21 d每天24 h的實際社區(qū)時用水量資料，如圖2所示。將樣本集21 d/h用水量中的前20 d 480組時用水量資料作為訓(xùn)練集，最后一天24 h用水量資料作為測試集。

圖2 單日時用水量變化趨勢Fig.2 Change in water consumption on a single day

2.3 模型輸入和輸出

從圖3可以明顯看出居民社區(qū)時用水量具有明顯的24 h周期性，為預(yù)測下一時刻用水量，模型的輸入取該時刻前24 h的時用水量。模型輸出是居民社區(qū)下一時刻用水量，輸入是該時刻前24 h每小時用水量(數(shù)據(jù)來源于2016年6月17日到7月6日居民社區(qū)的用水數(shù)據(jù))。

圖3 不同模型時用水量預(yù)測結(jié)果Fig.3 Water consumption prediction results of different models

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

慣性權(quán)重ω、分離權(quán)重s、對齊權(quán)重a、內(nèi)聚權(quán)重c、食物吸引權(quán)重f和天敵排斥權(quán)重e均為自適應(yīng)線性遞減權(quán)重，其最小值為0.4，最大值為0.9；設(shè)種群規(guī)模為10，最大迭代次數(shù)為20，參數(shù)M的取值范圍為[-1，1]，參數(shù)σf的取值范圍為[-1，1]，參數(shù)σn的取值范圍為[-1，1]。

2.5 對比模型及評價標準

為了使DA-GPR模型的預(yù)測結(jié)果具有可比性，選著BP神經(jīng)網(wǎng)絡(luò)、最小二乘支持向量機和高斯過程回歸模型作為對比模型，模型性能采用平均相對誤差(MAPE)、平均絕對誤差(MAE)和均方誤差(MSE)進行衡量，它們的定義為：

(18)

(19)

(20)

2.6 結(jié)果與分析

為了驗證本文創(chuàng)建的蜻蜓-高斯過程回歸耦合模型的預(yù)測可靠性，本文將其與BP神經(jīng)網(wǎng)絡(luò)、最小二乘支持向量機和高斯過程回歸預(yù)測模型進行試驗對比，預(yù)測結(jié)果如圖3所示。在95%置信度下，高斯過程回歸和蜻蜓-高斯過程回歸兩個模型的區(qū)間預(yù)測結(jié)果如圖4和圖5所示。從圖3可知，BP神經(jīng)網(wǎng)絡(luò)和最小二乘支持向量機模型預(yù)測結(jié)果誤差較大，蜻蜓-高斯過程回歸模型預(yù)測的結(jié)果誤差較小。從圖4可知，高斯過程回歸模型預(yù)測的結(jié)果誤差較大，極個別點落在了預(yù)測區(qū)間外。從圖5可知，蜻蜓-高斯過程回歸預(yù)測模型預(yù)測的結(jié)果走勢與實測值一致，預(yù)測值與實測值基本上重合，誤差較小。

圖4 高斯過程回歸(GPR)時用水量預(yù)測結(jié)果Fig.4 Water consumption prediction results in Gaussian process regression (GPR)

圖5 蜻蜓-高斯過程回歸(DA-GPR)時用水量預(yù)測結(jié)果Fig.5 Prediction of water consumption in the dragonfly Gauss process regression (DA-GPR)

不同預(yù)測模型的預(yù)測相對誤差如圖6所示。從圖6可以看出，BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型相對誤差最大為0.5，最小二乘支持向量機預(yù)測模型相對誤差最大為0.9，高斯過程回歸預(yù)測模型相對誤差最大為1.653，蜻蜓-高斯過程回歸預(yù)測模型相對誤差最大為0.019。顯然，蜻蜓-高斯過程回歸預(yù)測模型具有更高的預(yù)測精度。

圖6 不同預(yù)測模型預(yù)測的相對誤差Fig.6 Relative error of prediction model of different prediction models

采用平均相對誤差(MAPE)、平均絕對誤差(MAE)和均方誤差(MSE)對不同模型的預(yù)測結(jié)果進行分析，結(jié)果如表1所示。通過表1可以看出，蜻蜓-高斯過程回歸模型的居民社區(qū)時用水量預(yù)測精度優(yōu)于BP神經(jīng)網(wǎng)絡(luò)、最小二乘支持向量機和高斯過程回歸模型，取得了較好的預(yù)測效果。

表1 不同預(yù)測模型誤差分析Tab.1 Error analysis of different prediction models

3 結(jié) 語

居民社區(qū)時用水量具有較強的不確定性和隨機性，傳統(tǒng)的ANN、SVM等預(yù)測方法只能得到確定的點預(yù)測結(jié)果和未來某一時刻的預(yù)測結(jié)果，無法給出預(yù)測的區(qū)間，也不能進行實時預(yù)測。為了克服這些缺點，本文提出了一種基于蜻蜓-高斯過程回歸耦合的居民社區(qū)時用水量動態(tài)實時區(qū)間預(yù)測方法。為進一步提高預(yù)測精度，進行了改進，最終得到了一定置信水平下的區(qū)間預(yù)測結(jié)果。仿真結(jié)果表明，本文構(gòu)建的區(qū)間預(yù)測方法與常規(guī)方法相比，不僅能夠預(yù)測未知量的期望值，還能給出其分布狀況,同時也能進行實時預(yù)測。而且預(yù)測精度較高，最大的相對誤差為僅0.019，具有較強的實用價值，為未來水資源實時調(diào)度提供理論依據(jù)。