程亞麗， 王致杰， 江秀臣

(1. 上海電機(jī)學(xué)院 電氣學(xué)院，上海 201306； 2. 上海交通大學(xué) 電氣工程學(xué)院，上海 200240)

目前，對(duì)功率爬坡事件的研究主要集中于爬坡預(yù)測(cè)上[1]，而爬坡事件對(duì)電網(wǎng)產(chǎn)生的風(fēng)險(xiǎn)的研究很少，且只是對(duì)爬坡單一特征量進(jìn)行分析[2]，已經(jīng)很難滿足對(duì)于極端氣候變化特征和電網(wǎng)風(fēng)險(xiǎn)管理的需求。因此，有必要對(duì)爬坡要素間的相互關(guān)系和聯(lián)合概率分布特征進(jìn)行研究。Sklar[3]在1959年提出的Archimedean Copula函數(shù)，是一類將變量聯(lián)合分布函數(shù)與邊緣分布函數(shù)連接在一起的函數(shù)，很好地描述變量之間的相關(guān)性并計(jì)算其聯(lián)合概率分布。聯(lián)合分布里包含了所有風(fēng)險(xiǎn)要素的信息，可以全面評(píng)估事件的風(fēng)險(xiǎn)大小，被廣泛地用于水文事件、干旱、沙塵暴等方面的多變量研究中[4-10]。文獻(xiàn)[11-13]把Copula理論運(yùn)用到了金融市場(chǎng)的風(fēng)險(xiǎn)分析中，分析風(fēng)險(xiǎn)因子之間的相關(guān)性；文獻(xiàn)[14-16]結(jié)合Copula理論對(duì)投資組合風(fēng)險(xiǎn)進(jìn)行度量，能夠更準(zhǔn)確地度量投資組合的風(fēng)險(xiǎn)值。但是將此函數(shù)運(yùn)用在風(fēng)電場(chǎng)功率爬坡風(fēng)險(xiǎn)評(píng)估方面較少，本文利用Archimedean Copula函數(shù)多維聯(lián)合的靈活性，構(gòu)建風(fēng)電功率爬坡事件的二維聯(lián)合分布，建立了萊州風(fēng)電場(chǎng)發(fā)生功率爬坡的概率模型，定量分析風(fēng)電場(chǎng)發(fā)生功率爬坡的風(fēng)險(xiǎn)，為風(fēng)電場(chǎng)功率爬坡風(fēng)險(xiǎn)評(píng)估研究提供新的思路。

1 風(fēng)電場(chǎng)功率爬坡特征量的選取

在短時(shí)間內(nèi)風(fēng)電功率的劇烈波動(dòng)使得風(fēng)力發(fā)電具有不確定性、間斷性和隨機(jī)性，這種現(xiàn)象被稱為風(fēng)電功率爬坡事件，描述一個(gè)爬坡事件的主要特征量是確定的，包括爬坡幅值、爬坡方向、爬坡率、爬坡起止時(shí)間和持續(xù)時(shí)間，如圖1所示。爬坡幅值通常是在爬坡持續(xù)時(shí)間內(nèi)最大和最小功率的差值，爬坡方向用爬坡幅值的正負(fù)表示，正代表上爬坡事件，負(fù)代表下爬坡事件。

圖1 爬坡事件特征圖

由圖1可見(jiàn)，爬坡幅值越大，爬坡持續(xù)時(shí)間越短，則爬坡事件越嚴(yán)重。因此，選取爬坡幅值和爬坡持續(xù)時(shí)間這兩個(gè)特征量作為風(fēng)電功率爬坡事件的主要特征量，建立二維聯(lián)合分布。

2 Archimedean Copula函數(shù)

Archimedean Copula函數(shù)主要包括Frank Copula、Clayton Copula和Gumbel Copula函數(shù)，本文選用3種二維Archimedean Copula函數(shù)進(jìn)行爬坡特征量二維聯(lián)合，具體的分布函數(shù)及參數(shù)范圍如表1所示。

對(duì)于Archimedean Copula函數(shù)擬合度的檢驗(yàn)方法通常采用離差平方和(Oridinary Least Square,OLS)準(zhǔn)則法，采用文獻(xiàn)[5]中的OLS準(zhǔn)則法公式。該函數(shù)參數(shù)擬合度優(yōu)劣的判斷依據(jù)是OLS越小，表明該Archimedean Copula與其經(jīng)驗(yàn)點(diǎn)距的誤差越小，函數(shù)擬合程度越好。

表1 二維Archimedean Copula分布函數(shù)及參數(shù)范圍

表1中，u1，u2為Archimedean Copula函數(shù)中的兩個(gè)隨機(jī)變量。

3 基于Archimedean Copula函數(shù)的爬坡特征量聯(lián)合分布

3.1 爬坡特征量邊緣分布函數(shù)

通過(guò)對(duì)目前國(guó)內(nèi)外概率分布函數(shù)的研究，本文假定爬坡持續(xù)時(shí)間服從對(duì)數(shù)正態(tài)分布，風(fēng)電功率爬坡幅值服從廣義極值分布，則相應(yīng)的分布函數(shù)和概率密度函數(shù)設(shè)計(jì)如下：

設(shè)隨機(jī)變量x的樣本點(diǎn)為x1,x2,…,xn，隨機(jī)變量指的是風(fēng)電功率爬坡持續(xù)時(shí)間的實(shí)際數(shù)據(jù)，對(duì)數(shù)正態(tài)分布的概率密度函數(shù)為

(1)

風(fēng)電功率爬坡持續(xù)時(shí)間t的累積分布函數(shù)為

(2)

式中：μ為隨機(jī)變量的均值，也稱為位置參數(shù)；σ2為標(biāo)準(zhǔn)差，也稱為形狀參數(shù)。

極值分布是指在概率論中極大值(或者極小值)的概率分布，概率密度函數(shù)為

(3)

爬坡幅值ΔP的累積分布函數(shù)為

(4)

式中，α為尺度參數(shù)。

利用極大似然函數(shù)進(jìn)行邊緣分布參數(shù)估計(jì)，爬坡持續(xù)時(shí)間的參數(shù)估計(jì)計(jì)算公式如下：

(5)

對(duì)式(5)兩邊取對(duì)數(shù)，得

(6)

根據(jù)式(6)建立似然函數(shù)組為

(7)

可得

式中：L(μ，σ2)為似然函數(shù)；ti為爬坡持續(xù)時(shí)間變量。爬坡幅值分布函數(shù)的參數(shù)估計(jì)按照式(5)～式(7)類推。

4 實(shí)例研究

本章算例采用的數(shù)據(jù)是2015年一年萊州風(fēng)電場(chǎng)的實(shí)測(cè)數(shù)據(jù)，風(fēng)力發(fā)電機(jī)組的輸出功率每15 min采樣一次，分別采用對(duì)數(shù)正態(tài)分布、廣義極值分布對(duì)爬坡持續(xù)時(shí)間和爬坡幅值進(jìn)行擬合，利用極大似然法對(duì)各邊緣分布函數(shù)進(jìn)行參數(shù)估計(jì)，并通過(guò)K-S法對(duì)各邊緣分布函數(shù)進(jìn)行檢驗(yàn)。采用Archimedean Copula函數(shù)構(gòu)建風(fēng)電功率爬坡幅值和爬坡持續(xù)時(shí)間的聯(lián)合概率分布模型，并根據(jù)OLS準(zhǔn)則法最小原則，選取最合適的Archimedean Copula函數(shù)構(gòu)建爬坡幅值和爬坡持續(xù)時(shí)間的聯(lián)合概率分布函數(shù)。

4.1 爬坡特征量擬合分布

運(yùn)用極大似然法分別對(duì)爬坡持續(xù)時(shí)間、幅值的邊緣分布函數(shù)參數(shù)進(jìn)行估計(jì)，結(jié)果如表2所示。利用K-S法對(duì)各分布函數(shù)進(jìn)行檢驗(yàn)，取K-S檢驗(yàn)的顯著性水平α=0.05，當(dāng)n=44時(shí)，查詢檢驗(yàn)分位數(shù)表得到對(duì)應(yīng)的分位點(diǎn)D0=0.200 56，當(dāng)邊緣分布的檢驗(yàn)統(tǒng)計(jì)量D<0.200 56時(shí)，表明K-S檢驗(yàn)通過(guò)且擬合情況較好。爬坡幅值、持續(xù)時(shí)間的K-S檢驗(yàn)的計(jì)算結(jié)果如表3所示，對(duì)數(shù)正態(tài)分布和廣義極值分布對(duì)于爬坡特征量的擬合檢驗(yàn)值都通過(guò)了0.05的顯著性檢驗(yàn)，說(shuō)明上述兩種概率分布函數(shù)對(duì)于爬坡特征量的擬合效果較好。因此，選取廣義極值分布進(jìn)行爬坡幅值邊緣分布的擬合，選取對(duì)數(shù)正態(tài)分布進(jìn)行爬坡持續(xù)時(shí)間邊緣分布的擬合。

表2 邊緣分布函數(shù)的參數(shù)值

表3 爬坡特征變量的K-S檢驗(yàn)

4.2 爬坡特征量二維聯(lián)合分布

對(duì)爬坡特征量中的爬坡持續(xù)時(shí)間和爬坡幅值構(gòu)建Archimedean Copula函數(shù)二維聯(lián)合分布，利用式(5) ～式(7) 進(jìn)行參數(shù)估計(jì)，采用表1中的函數(shù)公式進(jìn)行二維擬合，通過(guò)OLS準(zhǔn)則法進(jìn)行優(yōu)度檢驗(yàn)，相關(guān)參數(shù)及檢驗(yàn)值見(jiàn)表4。

表4 分布函數(shù)參數(shù)及檢驗(yàn)值

從表4中可見(jiàn)，從K-S檢驗(yàn)統(tǒng)計(jì)量D值分析，3種Archimedean Copula聯(lián)合函數(shù)的D值都小于臨界值0.200 56。因此，3種函數(shù)都可以構(gòu)建萊州風(fēng)電場(chǎng)爬坡持續(xù)時(shí)間與爬坡幅值的聯(lián)合分布模型；從擬合優(yōu)度的評(píng)價(jià)指標(biāo)分析，Clayton Copula函數(shù)構(gòu)建的聯(lián)合分布模型的OLS值最小，擬合效果最優(yōu)。因此，選取Clayton Copula函數(shù)模型為萊州風(fēng)電場(chǎng)爬坡持續(xù)時(shí)間與爬坡幅值的聯(lián)合分布模型。

為了進(jìn)一步驗(yàn)證擬合效果的優(yōu)劣，利用Matlab軟件對(duì)萊州風(fēng)電場(chǎng)爬坡持續(xù)時(shí)間與爬坡幅值的數(shù)據(jù)編程，F(xiàn)rank Copula函數(shù)可以很好地描述對(duì)稱的相關(guān)關(guān)系，以Frank Copula函數(shù)為例，運(yùn)用Matlab求Frank Copula函數(shù)參數(shù)程序如下。

?

調(diào)用copulafit函數(shù)估計(jì)二元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ù)時(shí)間');

ylabel('爬坡幅值');

zlabel('二元Frank Copula密度函數(shù))');

?

通過(guò)Matlab仿真，得到爬坡持續(xù)時(shí)間與爬坡幅值的概率密度分布圖如圖2～4所示。Frank Copula函數(shù)具有對(duì)稱性，尾部很厚，可用于描述具有對(duì)稱的相關(guān)關(guān)系(見(jiàn)圖2)。由于風(fēng)電場(chǎng)發(fā)生嚴(yán)重爬坡事件和一般爬坡事件時(shí)不具有對(duì)稱性，因此該函數(shù)不能很好地描述風(fēng)電場(chǎng)功率爬坡嚴(yán)重事件引起的風(fēng)險(xiǎn)。Gumbel Copula函數(shù)分布似“J”形，擁有較厚的上尾部特征，上尾部特征是描述變量間同時(shí)出現(xiàn)大的變化的概率(見(jiàn)圖3)。風(fēng)電功率爬坡特征量的Gumbel Copula函數(shù)上尾部特征突出，不能很好地描述風(fēng)電場(chǎng)中爬坡幅值較大且爬坡持續(xù)時(shí)間較短的極端爬坡事件。Clayton Copula不具有對(duì)稱性，分布似“L”形，具有較厚的下尾部特征(見(jiàn)圖4)。風(fēng)電功率爬坡特征量的Clayton Copula函數(shù)具有很好的下尾部特征，能較好地描述極端爬坡事件小概率發(fā)生的現(xiàn)象，從而可以分析極端事件引起的損失。

圖2 Frank Copula函數(shù)的爬坡幅值和爬坡持續(xù)時(shí)間概率密度分布

圖3 Gumbel Copula函數(shù)的爬坡幅值和爬坡持續(xù)時(shí)間累計(jì)概率分布

圖4 Clayton Copula函數(shù)的爬坡幅值和爬坡持續(xù)時(shí)間概率密度分布

根據(jù)上述分析可知，可用Clayton Copula函數(shù)描述萊州風(fēng)電場(chǎng)功率爬坡幅值和爬坡持續(xù)時(shí)間的聯(lián)合分布，表達(dá)式如下：

FΔP,t=Cu1,u2=

(8)

萊州風(fēng)電場(chǎng)功率爬坡幅值和爬坡持續(xù)時(shí)間的聯(lián)合分布概率計(jì)算模型如下：

F(ΔP1? ΔP? ΔP2,t1?t?t2)=

F(ΔP2t2)-F(ΔP2,t1)-F(ΔP1,t2)+

F(ΔP1,t1)

(9)

該模型同時(shí)包含了爬坡幅值和爬坡持續(xù)時(shí)間信息，可準(zhǔn)確地計(jì)算爬坡事件發(fā)生的概率，能定量分析風(fēng)電功率爬坡事件的風(fēng)險(xiǎn)評(píng)估，方便調(diào)度人員更好地做出調(diào)度策略，減少風(fēng)險(xiǎn)引起的損失。

5 結(jié) 語(yǔ)

選取了風(fēng)電功率爬坡事件的爬坡持續(xù)時(shí)間和爬坡幅值兩個(gè)特征量，對(duì)其采用對(duì)數(shù)正態(tài)分布和廣義極值分布進(jìn)行邊緣分布擬合，并在建立風(fēng)電功率爬坡幅值和爬坡持續(xù)時(shí)間邊緣分布模型中，推導(dǎo)了利用極大似然法進(jìn)行參數(shù)估計(jì)的過(guò)程，最后采用K-S檢驗(yàn)法對(duì)擬合優(yōu)度檢驗(yàn)。在此基礎(chǔ)上，根據(jù)3種Archimedean Copula函數(shù)模型的表達(dá)式建立了爬坡幅值和爬坡持續(xù)時(shí)間的二維聯(lián)合分布函數(shù)，并采用OLS準(zhǔn)則法作為擬合優(yōu)度的評(píng)價(jià)指標(biāo)，確定了Clayton Copula模型作為萊州風(fēng)電場(chǎng)爬坡事件概率計(jì)算模型，從而定量地分析風(fēng)電場(chǎng)功率爬坡的風(fēng)險(xiǎn)。本文為風(fēng)電場(chǎng)功率爬坡風(fēng)險(xiǎn)評(píng)估方面提供一定的理論依據(jù)和支撐，拓展了Copula理論的應(yīng)用前景。