劉龍龍，李文竹，劉 心

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

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

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

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

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

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

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

1.1 DA-GPR基本原理

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

(1)

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

(2)

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

預(yù)測(cè)值的后驗(yàn)分布為：

(3)

(4)

(5)

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

(6)

1.2 DA-GPR目標(biāo)函數(shù)設(shè)計(jì)

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

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

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

(7)

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

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

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

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

(8)

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

(9)

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

(10)

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

Fi=X+-X

(11)

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

Ei=X-+X

(12)

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

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

在蜻蜓尋優(yōu)的過(guò)程中，大量的蜻蜓群體朝著同一個(gè)方向進(jìn)行長(zhǎng)距離遷移，這個(gè)群體活動(dòng)就是對(duì)算法的全局搜索過(guò)程。

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

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

(13)

有鄰近蜻蜓時(shí)，下一代蜻蜓的位置為：

Xt+1=Xt+ΔXt+1

(14)

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

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

(15)

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

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

(16)

(17)

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

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

1.5 DA-GPR算法步驟

流程圖如圖1所示。

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

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

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

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

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

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

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

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

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

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

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

2.1 仿真實(shí)驗(yàn)平臺(tái)

本次實(shí)驗(yàn)的運(yùn)行環(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ù)來(lái)自河北工程大學(xué)在線水量檢測(cè)平臺(tái)，采用居民家屬院2016年6月17日到7月7日，共21 d每天24 h的實(shí)際社區(qū)時(shí)用水量資料，如圖2所示。將樣本集21 d/h用水量中的前20 d 480組時(shí)用水量資料作為訓(xùn)練集，最后一天24 h用水量資料作為測(cè)試集。

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

2.3 模型輸入和輸出

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

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

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

慣性權(quán)重ω、分離權(quán)重s、對(duì)齊權(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 對(duì)比模型及評(píng)價(jià)標(biāo)準(zhǔn)

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

(18)

(19)

(20)

2.6 結(jié)果與分析

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

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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