程亞麗, 王致杰, 江秀臣
(1. 上海電機學院 電氣學院,上海 201306; 2. 上海交通大學 電氣工程學院,上海 200240)
目前,對功率爬坡事件的研究主要集中于爬坡預測上[1],而爬坡事件對電網(wǎng)產(chǎn)生的風險的研究很少,且只是對爬坡單一特征量進行分析[2],已經(jīng)很難滿足對于極端氣候變化特征和電網(wǎng)風險管理的需求。因此,有必要對爬坡要素間的相互關系和聯(lián)合概率分布特征進行研究。Sklar[3]在1959年提出的Archimedean Copula函數(shù),是一類將變量聯(lián)合分布函數(shù)與邊緣分布函數(shù)連接在一起的函數(shù),很好地描述變量之間的相關性并計算其聯(lián)合概率分布。聯(lián)合分布里包含了所有風險要素的信息,可以全面評估事件的風險大小,被廣泛地用于水文事件、干旱、沙塵暴等方面的多變量研究中[4-10]。文獻[11-13]把Copula理論運用到了金融市場的風險分析中,分析風險因子之間的相關性;文獻[14-16]結合Copula理論對投資組合風險進行度量,能夠更準確地度量投資組合的風險值。但是將此函數(shù)運用在風電場功率爬坡風險評估方面較少,本文利用Archimedean Copula函數(shù)多維聯(lián)合的靈活性,構建風電功率爬坡事件的二維聯(lián)合分布,建立了萊州風電場發(fā)生功率爬坡的概率模型,定量分析風電場發(fā)生功率爬坡的風險,為風電場功率爬坡風險評估研究提供新的思路。
在短時間內(nèi)風電功率的劇烈波動使得風力發(fā)電具有不確定性、間斷性和隨機性,這種現(xiàn)象被稱為風電功率爬坡事件,描述一個爬坡事件的主要特征量是確定的,包括爬坡幅值、爬坡方向、爬坡率、爬坡起止時間和持續(xù)時間,如圖1所示。爬坡幅值通常是在爬坡持續(xù)時間內(nèi)最大和最小功率的差值,爬坡方向用爬坡幅值的正負表示,正代表上爬坡事件,負代表下爬坡事件。
圖1 爬坡事件特征圖
由圖1可見,爬坡幅值越大,爬坡持續(xù)時間越短,則爬坡事件越嚴重。因此,選取爬坡幅值和爬坡持續(xù)時間這兩個特征量作為風電功率爬坡事件的主要特征量,建立二維聯(lián)合分布。
Archimedean Copula函數(shù)主要包括Frank Copula、Clayton Copula和Gumbel Copula函數(shù),本文選用3種二維Archimedean Copula函數(shù)進行爬坡特征量二維聯(lián)合,具體的分布函數(shù)及參數(shù)范圍如表1所示。
對于Archimedean Copula函數(shù)擬合度的檢驗方法通常采用離差平方和(Oridinary Least Square,OLS)準則法,采用文獻[5]中的OLS準則法公式。該函數(shù)參數(shù)擬合度優(yōu)劣的判斷依據(jù)是OLS越小,表明該Archimedean Copula與其經(jīng)驗點距的誤差越小,函數(shù)擬合程度越好。
表1 二維Archimedean Copula分布函數(shù)及參數(shù)范圍
表1中,u1,u2為Archimedean Copula函數(shù)中的兩個隨機變量。
通過對目前國內(nèi)外概率分布函數(shù)的研究,本文假定爬坡持續(xù)時間服從對數(shù)正態(tài)分布,風電功率爬坡幅值服從廣義極值分布,則相應的分布函數(shù)和概率密度函數(shù)設計如下:
設隨機變量x的樣本點為x1,x2,…,xn,隨機變量指的是風電功率爬坡持續(xù)時間的實際數(shù)據(jù),對數(shù)正態(tài)分布的概率密度函數(shù)為
(1)
風電功率爬坡持續(xù)時間t的累積分布函數(shù)為
(2)
式中:μ為隨機變量的均值,也稱為位置參數(shù);σ2為標準差,也稱為形狀參數(shù)。
極值分布是指在概率論中極大值(或者極小值)的概率分布,概率密度函數(shù)為
(3)
爬坡幅值ΔP的累積分布函數(shù)為
(4)
式中,α為尺度參數(shù)。
利用極大似然函數(shù)進行邊緣分布參數(shù)估計,爬坡持續(xù)時間的參數(shù)估計計算公式如下:
(5)
對式(5)兩邊取對數(shù),得
(6)
根據(jù)式(6)建立似然函數(shù)組為
(7)
可得
式中:L(μ,σ2)為似然函數(shù);ti為爬坡持續(xù)時間變量。爬坡幅值分布函數(shù)的參數(shù)估計按照式(5)~式(7)類推。
本章算例采用的數(shù)據(jù)是2015年一年萊州風電場的實測數(shù)據(jù),風力發(fā)電機組的輸出功率每15 min采樣一次,分別采用對數(shù)正態(tài)分布、廣義極值分布對爬坡持續(xù)時間和爬坡幅值進行擬合,利用極大似然法對各邊緣分布函數(shù)進行參數(shù)估計,并通過K-S法對各邊緣分布函數(shù)進行檢驗。采用Archimedean Copula函數(shù)構建風電功率爬坡幅值和爬坡持續(xù)時間的聯(lián)合概率分布模型,并根據(jù)OLS準則法最小原則,選取最合適的Archimedean Copula函數(shù)構建爬坡幅值和爬坡持續(xù)時間的聯(lián)合概率分布函數(shù)。
運用極大似然法分別對爬坡持續(xù)時間、幅值的邊緣分布函數(shù)參數(shù)進行估計,結果如表2所示。利用K-S法對各分布函數(shù)進行檢驗,取K-S檢驗的顯著性水平α=0.05,當n=44時,查詢檢驗分位數(shù)表得到對應的分位點D0=0.200 56,當邊緣分布的檢驗統(tǒng)計量D<0.200 56時,表明K-S檢驗通過且擬合情況較好。爬坡幅值、持續(xù)時間的K-S檢驗的計算結果如表3所示,對數(shù)正態(tài)分布和廣義極值分布對于爬坡特征量的擬合檢驗值都通過了0.05的顯著性檢驗,說明上述兩種概率分布函數(shù)對于爬坡特征量的擬合效果較好。因此,選取廣義極值分布進行爬坡幅值邊緣分布的擬合,選取對數(shù)正態(tài)分布進行爬坡持續(xù)時間邊緣分布的擬合。
表2 邊緣分布函數(shù)的參數(shù)值
表3 爬坡特征變量的K-S檢驗
對爬坡特征量中的爬坡持續(xù)時間和爬坡幅值構建Archimedean Copula函數(shù)二維聯(lián)合分布,利用式(5) ~式(7) 進行參數(shù)估計,采用表1中的函數(shù)公式進行二維擬合,通過OLS準則法進行優(yōu)度檢驗,相關參數(shù)及檢驗值見表4。
表4 分布函數(shù)參數(shù)及檢驗值
從表4中可見,從K-S檢驗統(tǒng)計量D值分析,3種Archimedean Copula聯(lián)合函數(shù)的D值都小于臨界值0.200 56。因此,3種函數(shù)都可以構建萊州風電場爬坡持續(xù)時間與爬坡幅值的聯(lián)合分布模型;從擬合優(yōu)度的評價指標分析,Clayton Copula函數(shù)構建的聯(lián)合分布模型的OLS值最小,擬合效果最優(yōu)。因此,選取Clayton Copula函數(shù)模型為萊州風電場爬坡持續(xù)時間與爬坡幅值的聯(lián)合分布模型。
為了進一步驗證擬合效果的優(yōu)劣,利用Matlab軟件對萊州風電場爬坡持續(xù)時間與爬坡幅值的數(shù)據(jù)編程,F(xiàn)rank Copula函數(shù)可以很好地描述對稱的相關關系,以Frank Copula函數(shù)為例,運用Matlab求Frank Copula函數(shù)參數(shù)程序如下。
?
調(diào)用copulafit函數(shù)估計二元Frank Copula中的參數(shù)
rho_norm = copulafit('Frank',[U(:), V(:)])
[Udata,Vdata] = meshgrid(linspace(0,1,31));
Cpdf_norm = copulapdf('Frank',[Udata(:), Vdata(:)],rho_norm);
Ccdf_t = copulacdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);
繪制二元Frank Copula的密度函數(shù)圖
figure;
surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata)));
xlabel('爬坡持續(xù)時間');
ylabel('爬坡幅值');
zlabel('二元Frank Copula密度函數(shù))');
?
通過Matlab仿真,得到爬坡持續(xù)時間與爬坡幅值的概率密度分布圖如圖2~4所示。Frank Copula函數(shù)具有對稱性,尾部很厚,可用于描述具有對稱的相關關系(見圖2)。由于風電場發(fā)生嚴重爬坡事件和一般爬坡事件時不具有對稱性,因此該函數(shù)不能很好地描述風電場功率爬坡嚴重事件引起的風險。Gumbel Copula函數(shù)分布似“J”形,擁有較厚的上尾部特征,上尾部特征是描述變量間同時出現(xiàn)大的變化的概率(見圖3)。風電功率爬坡特征量的Gumbel Copula函數(shù)上尾部特征突出,不能很好地描述風電場中爬坡幅值較大且爬坡持續(xù)時間較短的極端爬坡事件。Clayton Copula不具有對稱性,分布似“L”形,具有較厚的下尾部特征(見圖4)。風電功率爬坡特征量的Clayton Copula函數(shù)具有很好的下尾部特征,能較好地描述極端爬坡事件小概率發(fā)生的現(xiàn)象,從而可以分析極端事件引起的損失。
圖2 Frank Copula函數(shù)的爬坡幅值和爬坡持續(xù)時間概率密度分布
圖3 Gumbel Copula函數(shù)的爬坡幅值和爬坡持續(xù)時間累計概率分布
圖4 Clayton Copula函數(shù)的爬坡幅值和爬坡持續(xù)時間概率密度分布
根據(jù)上述分析可知,可用Clayton Copula函數(shù)描述萊州風電場功率爬坡幅值和爬坡持續(xù)時間的聯(lián)合分布,表達式如下:
FΔP,t=Cu1,u2=
(8)
萊州風電場功率爬坡幅值和爬坡持續(xù)時間的聯(lián)合分布概率計算模型如下:
F(ΔP1? ΔP? ΔP2,t1?t?t2)=
F(ΔP2t2)-F(ΔP2,t1)-F(ΔP1,t2)+
F(ΔP1,t1)
(9)
該模型同時包含了爬坡幅值和爬坡持續(xù)時間信息,可準確地計算爬坡事件發(fā)生的概率,能定量分析風電功率爬坡事件的風險評估,方便調(diào)度人員更好地做出調(diào)度策略,減少風險引起的損失。
選取了風電功率爬坡事件的爬坡持續(xù)時間和爬坡幅值兩個特征量,對其采用對數(shù)正態(tài)分布和廣義極值分布進行邊緣分布擬合,并在建立風電功率爬坡幅值和爬坡持續(xù)時間邊緣分布模型中,推導了利用極大似然法進行參數(shù)估計的過程,最后采用K-S檢驗法對擬合優(yōu)度檢驗。在此基礎上,根據(jù)3種Archimedean Copula函數(shù)模型的表達式建立了爬坡幅值和爬坡持續(xù)時間的二維聯(lián)合分布函數(shù),并采用OLS準則法作為擬合優(yōu)度的評價指標,確定了Clayton Copula模型作為萊州風電場爬坡事件概率計算模型,從而定量地分析風電場功率爬坡的風險。本文為風電場功率爬坡風險評估方面提供一定的理論依據(jù)和支撐,拓展了Copula理論的應用前景。