吳小濤，袁曉輝，袁艷斌，易凡茹，朱婧巍，吳育聯

（1.黃岡師范學院 數學與統(tǒng)計學院，湖北 黃岡 438000；2.華中科技大學 土木與水利工程學院，湖北 武漢 430074；3.武漢理工大學 資源與環(huán)境工程學院，湖北 武漢 430070）

0 引言

影響光伏發(fā)電系統(tǒng)輸出功率的主要因素為太陽輻照度。由于受溫度、濕度和云層等因素的影響，太陽輻照度具有非線性，這導致光伏發(fā)電系統(tǒng)輸出功率具有不穩(wěn)定性、間歇性等特點。建立太陽輻照度預測模型并獲得模擬數據，有利于電力調度部門根據模擬數據提前調整電能的存儲和釋放，實現光伏發(fā)電系統(tǒng)輸出功率平滑輸出，保持電網穩(wěn)定工作，并且能夠實現更多的光伏發(fā)電系統(tǒng)并網運行。

目前，太陽輻照度的預測模型主要包括基于太陽輻射傳遞的物理模型、基于衛(wèi)星圖像技術的預測模型和基于大數據驅動的統(tǒng)計預測模型。基于太陽輻射傳遞的物理模型建模過程復雜、抗干擾能力較差、魯棒性不強；基于衛(wèi)星圖像技術的預測模型技術性較強、建設成本較高；基于大數據驅動的統(tǒng)計預測模型能夠通過挖掘太陽輻照度及其影響因素的歷史數據之間的內在規(guī)律，利用機器學習等算法建立太陽輻照度預測模型（以下簡稱為機器學習模型），與基于太陽輻射傳遞的物理模型、基于衛(wèi)星圖像技術的預測模型相比，機器學習模型具有投入成本較低、預測速度較快等優(yōu)點，是目前主要的研究方向。

機器學習模型包括時間序列模型、神經網絡模型和支持向量機模型等，其中，時間序列模型適用于太陽輻照度變化較小的情形，當太陽輻照度變化較大時，由于神經網絡模型能解決非線性問題，其預測性能較好。神經網絡模型的經驗風險最小，但存在收斂速度慢、泛化能力弱且易陷入局部極小等問題；支持向量機模型的結構風險最小，具有很強的泛化能力，并且能夠避免陷入局部極小等問題，其預測性能往往優(yōu)于神經網絡模型[1]～[3]。文獻[4]對于同一組數據分別建立了自回歸模型、神經網絡模型和支持向量機模型，通過模擬計算發(fā)現，支持向量機模型的預測精度較高。最小二乘支持向量機（LSSVM）是一種改進的支持向量機，其求解速度優(yōu)于支持向量機。文獻[5]建立了計算太陽輻照度的LSSVM和人工神經網絡預測模型，模擬結果表明，LSSVM模型的預測速度、精度均優(yōu)于人工神經網絡模型。LSSVM模型中的核函數和核函數中的參數決定著LSSVM模型的預測性能。

本文將泛化能力強、推廣能力弱的高斯核函數和推廣能力強、泛化能力弱的多項式核函數相結合，構建了混合核函數，并采用蝗蟲優(yōu)化算法（GOA）優(yōu)化混合核函數中的參數，建立基于GOA優(yōu)化混合核LSSVM的太陽輻照度預測模型，用來預測太陽輻照度 由于太陽輻照度時間序列不穩(wěn)定，基于單一機器學習算法建立的太陽輻照度預測模型的預測精度較低。利用數據分解算法如小波分解（WD）算法、集合經驗模態(tài)分解（EEMD）算法和變分模態(tài)分解（VMD）算法等對太陽輻照度時間序列進行分解，可以有效地降低太陽輻照度時間序列的不穩(wěn)定性，因而，組合了機器學習和數據分解算法的太陽輻照度預測模型的預測精度高于基于單一機器學習算法的預測模型。

文獻[6]～[8]分別建立了基于WD和神經網絡、EEMD和神經網絡、WD和支持向量機預測模型，模擬結果表明，組合預測模型的預測精度高于單一神經網絡模型或支持向量機模型。不同的數據分解算法具有不同的特點，EEMD算法根據時間序列的極值分解太陽輻照度時間序列，雖然能有效分解出太陽輻照度時間序列的趨勢分量和細節(jié)分量，但是其隨機分量的頻率高于原始太陽輻照度時間序列，預測誤差會較大，從而降低了原始太陽輻照度時間序列的預測精度。VMD算法根據中心頻率分解太陽輻照度時間序列，可將太陽輻照度時間序列有效地分解為多個不同、相對穩(wěn)定的頻率段序列。

本文根據EEMD算法和VMD算法的特點，設計了基于EEMD和VMD的二階分解算法以分解太陽輻照度時間序列，具體的分解步驟:首先，采用EEMD算法分解原始太陽輻照度時間序列得到若干個分量；然后，利用VMD算法進一步分解頻率最高的分量，得到若干個相對穩(wěn)定的分量；接著，建立基于GOA優(yōu)化的混合核LSSVM的太陽輻照度預測模型，對所有分量進行預測，累加所有分量的預測結果，即得到原始太陽輻照度時間序列的預測結果。本文提出的預測模型簡稱為EEMD-VMD-GOA-LSSVM模型。

1 理論與方法

1.1 EEMD算法和VMD算法

EEMD算法為Huang在經驗模態(tài)分解（EMD）算法的基礎上提出的改進數據分解算法，該算法解決了EMD算法分解過程中極易出現的模態(tài)混疊問題[9]。

VMD算法為Dragomiretskiy K在2014年提出的一種數據分解算法[10]。VMD算法中假設原始序列y（t）由K個帶寬有限、頻率不同的分量uk（t）（k=1，2，…，K）組成，從而將數據分解問題轉化為變分問題，VMD算法包含變分模型的建立和求解2個步驟[11]。VMD算法在分解數據前須要確定模態(tài)個數K，若K小于真實值，會產生模態(tài)混疊現象；若K大于真實值，又會產生過分解現象；若K等于真實值，分解得到的殘差與原始序列的相關系數的絕對值較小，因此，本文提出設定某個閾值，即設定K的初始值為2，當分解得到的殘差與原始序列的相關系數不小于閾值時，令K=K+1，直到分解得到的殘差與原始序列的相關系數小于閾值時，停止分解，此時K的值即為最優(yōu)的模態(tài)個數。

1.2 GOA

GOA為Shahrzad Saremi在2017年提出的一種智能優(yōu)化算法[12]。在蝗蟲集群過程中，任意2個蝗蟲個體之間有引力和斥力2種社會行為，蝗蟲個體之間距離越近，斥力越大，引力越??；距離越遠，引力越大，斥力越小。GOA模擬的斥力能夠實現全局探索，模擬的引力能夠實現局部開發(fā)。所有蝗蟲個體之間共享位置信息，通過不斷縮小搜索空間獲得最優(yōu)解。

令第i，j個蝗蟲之間的距離為dij。dij的計算式為

式中:qi為第i個蝗蟲的位置；qj為第j個蝗蟲的位置。

令第i，j個蝗蟲之間的社會行為強度為s（dij）。s（dij）的計算式為

式中:ζ，ξ均為常數，ζ=0.5，ξ=1.5。

令第i個蝗蟲個體的社會行為為Si。Si的計算式為

GOA的尋優(yōu)步驟如下。

①將含有n0個蝗蟲的群體初始化為一群隨機粒子（隨機解），令Yi（i=1，2，…，n0）為第i個粒子的當前位置，根據適應度函數計算每個粒子的適應度值，得到所有粒子當前的最優(yōu)位置T*。

②每個粒子利用式（5）更新后的位置Yi。Yi的計算式為

式中:ub，lb分別為搜索空間的上界和下界；c為平衡算法全局和局部搜索能力的線性遞減系數。

c的計算式為

式中:L，l分別為最大迭代次數和當前迭代數；cmin，cmax分別為c的最小值、最大值，取cmin=10-5，cmax=1。

③迭代數l=l+1，并利用式（6）更新c的值。

④計算更新位置后的粒子的適應度值及最優(yōu)位置，判斷最優(yōu)位置是否優(yōu)于T*。若優(yōu)于則用其替代T*，否則保持T*不變。在更新位置時，如果粒子的位置小于搜索空間的下界lb，則將粒子的位置置于lb；如果粒子的位置大于搜索空間的上界ub，則將粒子的位置置于ub。

⑤根據迭代終止條件（最大迭代次數或誤差值）判斷是否結束尋優(yōu)，若結束則返回T*，否則返回③繼續(xù)進行迭代。

在GOA尋優(yōu)時，c的取值對算法性能的影響較大。在迭代初期，智能優(yōu)化算法往往須要變化幅度較大的c值來實現全局搜索，確定最優(yōu)解的大致范圍；在后期須要變化幅度較小的c值實現局部精細開發(fā)確定最優(yōu)解。顯然線性遞減的c不能滿足上述要求，為此本文提出基于非線性的逆不完全伽馬函數的c值調整方法，即將式（6）表達式更換為

式中:gammaincinv（τ，1-l/L）為逆不完全伽馬函數在Matlab中的調用函數，返回在τ和1-l/L處計算的下不完全伽馬函數的逆函數，其中τ為隨機參數，取0.2。

與式（6）相比，式（7）在迭代初期具有變化幅度較大的c，利于GOA實現全局探索；在迭代后期，具有變化幅度較小的c，利于GOA實現局部精細開發(fā)，從而使GOA具有更好地平衡全局探索和局部精細開發(fā)的能力。

1.3 LSSVM

LSSVM采用映射將非線性問題轉化為高維特征空間中的線性問題[13]。假設樣本集為{（xi，yi）|i=1，2，3，…，N}，其中:xi為多維的輸入樣本數據；yi為與xi對應的一維輸出樣本數據；N為樣本個數。

構造線性回歸函數為

式中:w為高維特征空間中權系數列向量；b為偏差，b∈R；φ（x）為非線性映射，φ（x）=[φ（x1），φ（x2），…，φ（xN）]T。

w和b通過最小化結構風險函數R來確定。R的計算式為

式中:εi為誤差；γ為懲罰系數，γ＞0。

利用拉格朗日乘數法求解式（9）。拉格朗日函數L（w，b，ε，μ）的計算式為

式中:μi為拉格朗日乘子。

通過KKT條件，將式（10）轉化為求解線性方程組問題，轉化結果見式（11）。

本文選擇泛化能力強、推廣能力弱的高斯核函數和推廣能力強、泛化能力弱的多項式核函數構建混合核函數。

式中:Ker（xi，xj）為核函數，其中i，j=1，2，…，N。

混合核函數的表達式為

式中:β為權值，0≤β≤1；σ為高斯核函數的核寬，σ＞0；r，p為多項式核函數中的參數，r＞0，p為正整數。

求解式（11）得到μi和b，從而得到LSSVM的最優(yōu)回歸函數表達式為

根據回歸函數式（13）即可對未知點進行預測。

混合核函數中的參數β，σ，r，p會直接影響LSSVM模型的預測性能，本文采用GOA對上述參數進行優(yōu)化，從而建立最優(yōu)的混合核LSSVM的太陽輻照度預測模型。

2 EEMD-VMD-GOA-LSSVM預測模型

為了降低原始太陽輻照度時間序列的不穩(wěn)定性，本文利用EEMD-VMD二階分解算法對原始太陽輻照度時間序列進行分解，即先利用EEMD算法分解原始太陽輻照度時間序列，得到n個分量IMF1，IMF2，…，IMFn，然后利用VMD算法對頻率最高的分量（本文假設IMF1的頻率最高）進行分解，得到K個分量IMF11，IMF12，…，IMF1K，則原始太陽輻照度時間序列被分解成n-1+K個分量。

建立基于GOA優(yōu)化的混合核LSSVM的太陽輻照度預測模型對n-1+K個分量進行預測，累加所有預測結果即可得到原始太陽輻照度時間序列的預測結果。具體預測步驟如下。

①利用EEMD-VMD二階分解算法分解原始太陽輻照度時間序列，得到若干個分量。

②對每個分量，選擇最后若干個數據作為測試數據集，其余數據作為訓練數據集。

③建立混合核LSSVM的太陽輻照度初始預測模型，在訓練數據集中利用GOA對混合核函數的參數進行優(yōu)化，利用優(yōu)化后的混合核LSSVM的太陽輻照度預測模型對每個分量的測試數據集進行預測，得到每個分量的預測結果。

④將③中每個分量的預測結果進行疊加，得到原始太陽輻照度時間序列的預測值。

EEMD-VMD二階分解算法流程圖見圖1。

圖1 EEMD-VMD二階分解算法流程圖Fig.1 Two-stage decomposition algorithm flow diagramwith EEMD-VMD

3 實驗分析

3.1 數據來源及數據整理

由于太陽輻照度具有較強的季節(jié)周期性和晝夜周期性，本文選取某個太陽能監(jiān)控平臺收集的夏季每日8:00-16:00各整點時刻的太陽輻照度和部分氣象數據（溫度、相對濕度和氣壓），共810組數據。

太陽輻照度隨時間的變化情況如圖2所示。

圖2 太陽輻照度隨時間的變化情況Fig.2 Variation of solar irradiance with time

由圖2可知，太陽輻照度時間序列整體上表現出非線性的特點，這是由于太陽輻照度受到溫度、相對濕度和氣壓等因素的影響。

溫度、相對濕度和氣壓隨時間的變化情況如圖3所示。

圖3 溫度、相對濕度和氣壓隨時間的變化情況Fig.3 Variation of temperature，relative humidity and air pressure with time

由圖2、圖3（a）可知，溫度隨時間的變化情況與太陽輻照度隨時間的變化情況基本一致，即當溫度較低時，太陽輻照度較??；當溫度較高時，太陽輻照度較大。由圖2、圖3（b）可知，相對濕度隨時間的變化情況與太陽輻照度隨時間的變化情況基本相反。由圖2、圖3（c）可知，氣壓隨時間的變化情況與太陽輻照度隨時間的變化情況無關。進一步分析發(fā)現，對原始太陽輻照度時間序列進行異常值和歸一化處理后，利用皮爾遜相關系數計算太陽輻照度與溫度、相對濕度和氣壓的相關系數，得到相關系數分別為0.7，-0.7，0.1，由此可知，太陽輻照度與溫度呈中度正相關，與相對濕度呈中度負相關，與氣壓不相關。因此，本文中的氣象因素只須考慮溫度和相對濕度，在建立混合核LSSVM的太陽輻照度預測模型時，將溫度和相對濕度作為輸入，以提高模型的預測精度。

3.2 EEMD-VMD二階分解算法

利用EEMD算法對整理后的太陽輻照度時間序列進行分解，其中，噪聲方差Nstd取0.01，噪聲組數NE取100，得到9個分量IMF1～IMF9。利用EEMD算法分解原始太陽輻照度時間序列的結果如圖4所示。

圖4利用EEMD算法分解原始太陽輻照度時間序列的結果Fig.4 Decomposition results of original solar irradiance time series using EEMD algorithm

由圖4可知:IMF9反映了原始太陽輻照度時間序列的整體變化趨勢，變化頻率最小；IMF2～IMF8反映了原始太陽輻照度時間序列的細節(jié)變化情況，變化頻率雖然較大，但小于原始序列；IMF1的變化頻率最高。

利用VMD算法對IMF1進一步分解，其中，帶寬限制alpha取2 000，噪聲容限tau取1。對于模態(tài)個數K，按照前文提出的方法，設定相關系數的閾值為0.005，計算得到K=2，3，4，5，6，7時，IMF1與分解殘差的相關系數分別為-0.01，0.307，0.1 99，0.0 02，0.0 01，0.0 01。因此，取K=5，此時殘差與原始太陽輻照度時間序列的相關系數僅為0.0 02。當K取6，7或更大時，殘差與原始太陽輻照度時間序列的相關系數變化不大，再分解會產生過分解現象。利用VMD算法對IMF1進行分解得到5個分量和1個殘差分量，記為IMF11～IMF16。

綜上所述，利用EEMD-VMD二階分解算法分解原始太陽輻照度時間序列可得到14個分量，分別為IMF11～IMF16，IMF2～IMF9。利用VMD算法分解IMF1得到的結果如圖5所示。

圖5 利用VMD算法分解IMF 1得到的結果Fig.5 Decomposition result of IMF1 using VMD algorithm

3.3 LSSVM預測模型的輸入和輸出

本文擬對未來2 d中8:00-16:00各整點時刻的太陽輻照度進行短期預測。對于利用EEMDVMD二階分解算法分解得到的各分量，選擇最后18個數據（793～810 h）作為測試數據，前792個數據作為訓練數據，并建立基于GOA優(yōu)化的混合核LSSVM的太陽輻照度預測模型，該模型的輸入根據太陽輻照度時間序列的偏自相關系數的截尾特性來確定。具體地，如果太陽輻照度時間序列為J階截尾，那么將t時刻的太陽輻照度作為輸出，其對應的t-1，t-2，……，t-J這J個時刻的太陽輻照度、溫度和相對濕度作為輸入。對原始太陽輻照度時間序列建立基于GOA優(yōu)化的混合核LSSVM的太陽輻照度預測模型，通過對比輸入中是否包含溫度和相對濕度的實驗結果發(fā)現，包含溫度和相對濕度的預測結果的平均絕對誤差為13.0%，而不包含時預測結果的平均絕對誤差為16.5%，因此，輸入變量中包含溫度和相對濕度時，預測精度更高。

3.4 預測模型評價指標

本文采用決定系數R2、平均絕對誤差MAPE、平均相對誤差MRE和均方根誤差RMSE這4個指標對EEMD-VMD-GOA-LSSVM模型進行評價分析。其中:R2用來檢驗模型的擬合度，R2值越大，說明模型的擬合程度越好；MAPE，MRE和RMSE用來檢驗模型的精度，它們的值越小，說明模型的預測精度越高。R2，MAPE，MRE和RMSE的計算式分別為

式中:testi，prei分別為i時刻太陽輻照度的實測值、預測值；Num為預測樣本總數；testmax，testave分別為測試樣本中太陽輻照度實測值的最大值和平均值。

4 實驗結果分析

為了檢驗EEMD-VMD-GOA-LSSVM模型的預測性能，本文分別建立了ARIMA（p，d，q）模型、GOA-LSSVM模型、GOA-BP神經網絡模型和EEMD-GOA-LSSVM模型用于比較。ARIMA（p，d，q）模型中p，d，q分別為自回歸（AR）的項數，差分（I）的系數，移動平均（MA）的項數。在ARIMA（p，d，q）模型中，通過反復驗證得知，當p=3，d=1，q=3時，該模型的預測精度最高；在GOA-LSSVM模型中，核函數采用本文提出的混合核函數，GOA用來優(yōu)化核參數，適應度函數選擇平均絕對誤差函數，初始蝗蟲群體個數取10，最大迭代次數為10；在GOA-BP神經網絡模型中，GOA用來優(yōu)化BP，實驗得出隱含層、輸出層中神經元的個數分別為22，12時，該模型的預測精度最高。在EEMD-GOA-LSSVM模型中，利用EEMD算法分解原始太陽輻照度時間序列得到9個分量，即IMF1～IMF9，建立基于GOA優(yōu)化的混合核LSSVM預測模型，對9個分量進行預測，再累加這9個分量的預測值，得到原始太陽輻照度時間序列的預測結果。

5種預測模型的預測曲線圖如圖6所示。由圖可知，5種預測模型中，ARIMA（p，d，q）模型的預測精度最差，其次為EEMD-GOA-LSSVM模型，GOA-LSSVM和GOA-BP模型的預測精度稍好，其預測曲線整體上追隨實測值曲線波動，但在某些點的誤差較大，EEMD-VMD-OA-LSSVM模型的預測精度最高，其預測曲線與實測值曲線幾乎重合。

圖6 5種預測模型的預測曲線圖Fig.6 Prediction curves of five forecasting models

為了便于比較，本文將801～803 h，5種預測模型的預測曲線和實測值曲線局部放大，具體如圖7所示。

圖7 5種預測模型的預測曲線局部放大圖Fig.7 Local enlarged map of prediction curve of five prediction models

由圖7可知，在801，802 h和803 h這3個時刻，EEMD-VMD-GOA-LSSVM模型的預測值均小于實測值，但它們的差值較小，而其余4種模型的預測值與實測值之間的差值較大，在802 h時最大。

5種預測模型的誤差值如表1所示。

表1 5種預測模型的誤差值Table 1 Error values of five prediction models

由表1可知，EEMD-VMD-GOA-LSSVM模型的MAPE為2.3%，MRE為1.3%，這說明該模型的預測值與實測值間的誤差非常小，預測精度非常高，其余4種預測模型的MAPE均大于10%，MRE均大于5%，這說明它們的預測值與實測值間的誤差大，預測精度較低；EEMD-VMDGOA-LSSVM模型的RMSE為12.2，遠小于其余4種模型，這說明EEMD-VMD-GOA-LSSVM模型的預測值與實測值的離散程度更?。辉谀Ｐ蛿M合度方面，EEMD-VMD-GOA-LSSVM模型的R2值為0.99，而其余4種模型的R2值均小于0.9，這說明EEMD-VMD-GOA-LSSVM模型能更好地反映原始太陽輻照度時間序列的變化情況。綜合上述預測曲線圖和誤差的分析可知，EEMD-VMDGOA-LSSVM模型具有較好的預測性能。

5 結論

本文針對太陽輻照度時間序列不穩(wěn)定的特點，提出了一種基于EEMD-VMD二階數據分解算法和GOA優(yōu)化混合核LSSVM的太陽輻照度預測模型，并利用該模型對太陽輻照度進行了預測，得到如下結論。

①利用EEMD-VMD二階分解算法得到的分量不僅能反映原始太陽輻照度時間序列的變化情況，而且較原始太陽輻照度時間序列更加穩(wěn)定。沒有采用數據分解算法的GOA-BP模型、ARIMA模型和GOA-LSSVM模型的MAPE均高于13%，而采用EEMD-VMD二階分解算法的EEMDVMD-GOA-LSSVM模型的MAPE僅為2.3%，這說明組合了EEMD-VMD二階分解算法的太陽輻照度預測模型比未使用數據分解算法的太陽輻照度預測模型的預測精度更高。

②針對VMD算法中模態(tài)個數K難以確定的問題，本文提出了基于原始太陽輻照度時間序列和分解殘差的相關系數來確定的方法。

③與ARIMA模型、GOA-LSSVM模型、GOABP神經網絡模型和EEMD-GOA-LSSVM模型相比，EEMD-VMD-GOA-LSSVM模型的預測精度更高，能有效地反映太陽輻照度的變化規(guī)律。