基于Copula函數(shù)的風(fēng)速相關(guān)性建模及概率最優(yōu)潮流分析

徐玉琴，張林浩

(華北電力大學(xué) 電氣與電子工程學(xué)院，河北 保定 071003)

0 引言

針對日趨嚴(yán)重的環(huán)境問題與能源危機(jī)，風(fēng)電產(chǎn)業(yè)的開發(fā)利用受到了國家的高度重視，成為發(fā)展速度最快的新能源發(fā)電方式［1］。風(fēng)能具有隨機(jī)性、間歇性的特點(diǎn)，同時(shí)由于風(fēng)力資源的地域分布特征，風(fēng)電場往往以集群方式接入電網(wǎng)，這些特性對電力系統(tǒng)的穩(wěn)定運(yùn)行與經(jīng)濟(jì)調(diào)度產(chǎn)生了一定影響。

在現(xiàn)有的最優(yōu)潮流分析方法中，一般假定各風(fēng)電場之間的風(fēng)速是相互獨(dú)立的［2］，事實(shí)上，風(fēng)電場都選址在風(fēng)能資源充裕的地區(qū)，地理位置比較接近，其風(fēng)速具有較強(qiáng)的相關(guān)性，從而各風(fēng)電的出力也存在較強(qiáng)的相關(guān)性［3］。因此，在分析最優(yōu)潮流的過程中，應(yīng)該考慮到風(fēng)速相關(guān)性因素。由于多采用威布爾(Weibull)概率分布函數(shù)來描述風(fēng)速的分布特性，研究風(fēng)速相關(guān)性問題的關(guān)鍵在于正確處理非正態(tài)相關(guān)變量之間的相關(guān)性。文獻(xiàn)［4］采用了三階多項(xiàng)式正態(tài)變換技術(shù)來產(chǎn)生具有指定邊際概率分布及相關(guān)性的多變量隨機(jī)數(shù)，并應(yīng)用于概率最優(yōu)潮流分析。文獻(xiàn)［5］采用Nataf逆變換技術(shù)建立不同風(fēng)電場之間具有相關(guān)性的風(fēng)速分布樣本空間，進(jìn)而得到具有相關(guān)性的風(fēng)電場出力。在最近的研究中，Copula函數(shù)逐漸應(yīng)用于隨機(jī)變量相關(guān)性問題的處理。文獻(xiàn)［6］基于多個(gè)實(shí)測負(fù)荷水平的相關(guān)性，采用Copula函數(shù)建立了具有相關(guān)性的負(fù)荷需求樣本。文獻(xiàn)［7］介紹了基于Copula函數(shù)構(gòu)建電力系統(tǒng)隨機(jī)變量之間的聯(lián)合概率計(jì)算方法，并利用正態(tài)Copula函數(shù)描述離岸風(fēng)電場與近岸風(fēng)電場出力之間的相關(guān)性。與其他方法相比，利用Copula函數(shù)描述隨機(jī)變量相關(guān)性的優(yōu)點(diǎn)主要有兩點(diǎn):第一點(diǎn)是邊緣分布的選擇不受限制，即便是服從不同分布的隨機(jī)變量，其相關(guān)性也可以用Copula函數(shù)刻畫;第二點(diǎn)是Copula函數(shù)可以描述隨機(jī)變量之間的非線性相關(guān)性，并且對變量進(jìn)行非線性單調(diào)增變換時(shí)，由Copula函數(shù)導(dǎo)出的相關(guān)性測度不會發(fā)生改變。

本文利用Copula函數(shù)構(gòu)建多風(fēng)電場的風(fēng)速相關(guān)性模型，針對風(fēng)電的不確定性，采用基于機(jī)會約束規(guī)劃(Chance Constrained Programming，CCP)的電力系統(tǒng)概率最優(yōu)潮流(Probabilistically Optimal Power Flow，P－OPF)模型，利用一種基于隨機(jī)模擬技術(shù)的粒子群優(yōu)化(Particle Swarm Optimization，PSO)算法對上述模型進(jìn)行求解。最后，以實(shí)際風(fēng)電場的歷史風(fēng)速數(shù)據(jù)和IEEE－30節(jié)點(diǎn)測試系統(tǒng)為算例，驗(yàn)證本文所提方法的合理性與有效性。

1 基于Copula函數(shù)的風(fēng)電場相關(guān)風(fēng)速建模

1.1 Copula函數(shù)的概念

Copula是一個(gè)定義在N維［0，1］空間內(nèi)的多元同一分布函數(shù)［8］，其實(shí)際意義是一類將聯(lián)合分布函數(shù)與它們各自的邊緣分布函數(shù)連接在一起的函數(shù)，因此也稱為連接函數(shù)，其理論基礎(chǔ)是Sklar定理。

Sklar定理:令 F(·，…，·)為具有邊緣分布F1(·)，F(xiàn)2(·)，…，F(xiàn)N(·) 的聯(lián)合分布函數(shù)，那么存在一個(gè)Copula函數(shù)C(·，…，·)，滿足:

若 F1(·)，F(xiàn)2(·)，…，F(xiàn)N(·) 連續(xù)，則 C(·，…，·)唯一確定;反之，若 F1(·)，F(xiàn)2(·)，…，F(xiàn)N(·) 為一元分布，C(·，…，·)為相應(yīng)的 Copula函數(shù)，那么由式(1)定義的函數(shù)F(·，…，·)是具有邊緣分布F1(·)，F(xiàn)2(·)，…，F(xiàn)N(·) 的聯(lián)合分布函數(shù)。

可以看出，Copula函數(shù)建立了邊緣分布到聯(lián)合分布的映射關(guān)系，利用Copula函數(shù)，可以將邊緣分布與變量之間的相關(guān)結(jié)構(gòu)分離，從而減小含非正態(tài)相關(guān)變量模型分析的難度。

1.2 Copula函數(shù)的分類

常用的Copula函數(shù)分為正態(tài)Copula，t－Copula和阿基米德 Copula，其中阿基米德 Copula包括Gumbel－Copula，Clayton－Copula 以及 Frank－Copula，下面主要介紹阿基米德Copula函數(shù)的相關(guān)概念和特性。

阿基米德Copula(Archimedean Copula)分布函數(shù)的表達(dá)式為

式中:φ(·)稱為阿基米德Copula函數(shù)的生成元(generator)，由表達(dá)式可知，阿基米德Copula函數(shù)由其生成元唯一確定，三種常見阿基米德Copula函數(shù)的生成元如表1所示，其中t∈［0，1］，表示Copula函數(shù)所描述的邊緣分布函數(shù);θ代表生成元相關(guān)參數(shù)，與隨機(jī)變量的相關(guān)性程度有關(guān)。對應(yīng)于建立風(fēng)電場風(fēng)速聯(lián)合概率分布而言，t即為各風(fēng)電場風(fēng)速的邊緣分布;θ表征各風(fēng)電場風(fēng)速的相關(guān)性程度。

表1 阿基米德函數(shù)的生成元Tab.1 Archimedean Copulas generators

其中，Gumbel－Copula的密度函數(shù)圖像呈“J”字形，上尾高，下尾低，對變量在分布的上尾部的變化比較敏感;Clayton－Copula的密度函數(shù)圖像呈“L”字形，下尾高，上尾低，能夠捕捉到下尾相關(guān)的變化;Frank－Copula具有對稱性，尾部相關(guān)系數(shù)為0，適用于描述具有對稱相關(guān)結(jié)構(gòu)的隨機(jī)變量。

1.3 Copula函數(shù)的選取

選取適當(dāng)?shù)腃opula來描述多風(fēng)電場風(fēng)速的相關(guān)結(jié)構(gòu)，是得到相關(guān)性風(fēng)速樣本空間的首要步驟。有多種Copula函數(shù)的選取方法，本文采用基于經(jīng)驗(yàn)Copula(Empirical Copula，EMC)和理論Copula(Theoretical Copula，THC)的最短歐氏距離法［8］來對Copula函數(shù)進(jìn)行選取。

設(shè)N個(gè)風(fēng)電場的風(fēng)速向量為V=(V1，V2，…，VN)，其中 Vi(i=1，2，…，N) 的概率密度函數(shù)fi(vi)以及累積分布函數(shù)Fi(vi)已知，風(fēng)速向量的EMC可由下式求得:

式中:I(·)為示性函數(shù)，括號中的條件滿足時(shí)，I=1，反之為0;{(v1，j，v2，j，…，vN，j)，j=1，2，…，n}為容量為n的風(fēng)速觀測樣本;為順序統(tǒng)計(jì)量，1≤ i1，…，iN≤ n 。

得到EMC后，可由下式計(jì)算EMC與THC之間的歐式距離:

得到計(jì)算結(jié)果后，選取歐式距離最小的Copula函數(shù)作為描述風(fēng)電場間風(fēng)速相關(guān)性的連接函數(shù)。

1.4 基于Copula函數(shù)的相關(guān)性風(fēng)速

利用選取的Copula函數(shù)并結(jié)合式 (1)，可以得到各風(fēng)電場風(fēng)速的聯(lián)合概率分布函數(shù):

例如:選取Gumbel－Copula作為描述風(fēng)速相關(guān)性的最佳函數(shù)，將邊緣分布函數(shù)Fi(vi)代入Gumbel－Copula的生成元，可以得到具體的風(fēng)速聯(lián)合概率分布函數(shù):

式中:參數(shù)θ可用最大似然法進(jìn)行估計(jì)。

利用上述Copula函數(shù)相關(guān)理論，計(jì)算具有相關(guān)性的風(fēng)電場風(fēng)速樣本空間和出力，算法流程如圖1所示。

圖1 基于Copula函數(shù)的風(fēng)速和出力計(jì)算框圖Fig.1 Calculation diagram of wind speed and output power based on Copula function

2 基于機(jī)會約束規(guī)劃的最優(yōu)潮流模型

風(fēng)電場的出力具有隨機(jī)性，此時(shí)宜引入機(jī)會約束規(guī)劃理論，以解決約束條件中含有隨機(jī)變量，且必須在觀察到隨機(jī)變量的實(shí)現(xiàn)前做出決策的情況。考慮到所做決策在不利情況下可能使約束條件不滿足，而采取一種原則:即允許所做決策在一定程度上不滿足約束條件，但該決策應(yīng)使約束條件成立的概率不小于某一置信水平［9］。本文建立了風(fēng)速相關(guān)情況下基于機(jī)會約束規(guī)劃的含風(fēng)電場電力系統(tǒng)概率最優(yōu)潮流模型。

2.1 目標(biāo)函數(shù)

模型以常規(guī)機(jī)組的發(fā)電成本最低作為優(yōu)化目標(biāo)，表達(dá)式如下:

式中:Ng為常規(guī)機(jī)組數(shù);S為其中平衡機(jī)組數(shù);C(P)=αcP2+βcP+γc，為常規(guī)機(jī)組的發(fā)電成本函數(shù)，αc，βc，γc為相應(yīng)系數(shù);Pi為非平衡機(jī)組有功出力;為平衡機(jī)組有功出力。

由于風(fēng)電出力的不確定性，平衡機(jī)組為了配合風(fēng)電場的運(yùn)行，其發(fā)電功率也是隨機(jī)量。對于n個(gè)隨機(jī)風(fēng)速樣本，可計(jì)算出n個(gè)風(fēng)電場出力，忽略系統(tǒng)損耗，用負(fù)荷減去風(fēng)電場的發(fā)電功率和非平衡機(jī)組的發(fā)電功率，即可得到n個(gè)平衡機(jī)組出力，從而計(jì)算出n個(gè)目標(biāo)函數(shù)值，取其平均值作為對目標(biāo)函數(shù)的估計(jì)。

2.2 等式約束

等式約束為系統(tǒng)功率平衡約束，表達(dá)式如下:

式中:n 為系統(tǒng)總節(jié)點(diǎn)數(shù);Pg，i，Pw，i，Pl，i分別為節(jié)點(diǎn)i上的常規(guī)機(jī)組發(fā)電功率、風(fēng)電場發(fā)電功率以及負(fù)荷需求。

2.3 不等式約束

模型中的控制變量為常規(guī)機(jī)組的發(fā)電功率，其不等式約束為

本文考慮的狀態(tài)變量有節(jié)點(diǎn)電壓幅值、常規(guī)機(jī)組無功出力、線路傳輸功率以及系統(tǒng)的旋轉(zhuǎn)備用。由于風(fēng)電出力的不確定性，潮流計(jì)算后得到的狀態(tài)變量也都是隨機(jī)變量，利用機(jī)會約束規(guī)劃，將不等式約束條件的表達(dá)式表述為概率形式，使之在一定的置信水平上成立，具體如下:

式中:Pr{·}為 {·}中事件成立的概率;上標(biāo)max和min分別表示相應(yīng)狀態(tài)變量的上、下限;Ui為節(jié)點(diǎn)i的電壓幅值;Qg，i為機(jī)組i的無功出力;Pli為支路i的傳輸功率;SPQ為PQ節(jié)點(diǎn)集合，Sl為支路集合;Psr為系統(tǒng)的旋轉(zhuǎn)備用，可取系統(tǒng)總負(fù)荷的5%;β1～β4分別為相應(yīng)不等式約束的置信水平。

3 基于隨機(jī)模擬的粒子群優(yōu)化算法

3.1 隨機(jī)模擬技術(shù)

隨機(jī)模擬技術(shù)以大量試驗(yàn)作為基礎(chǔ)，對計(jì)算結(jié)果進(jìn)行概率形式的驗(yàn)證，為處理機(jī)會約束規(guī)劃問題提供了一條有效途徑，具體實(shí)現(xiàn)的方法如下。

考慮機(jī)會約束:

首先，設(shè)置N1=0，根據(jù)ξ的概率分布Φ(ξ)生成N個(gè)隨機(jī)變量，將控制變量x與生成的隨機(jī)變量代入式 (11)，若滿足不等式g(x，ξ)≤0，則N1=N1+1，根據(jù)大數(shù)定律，當(dāng)N足夠大時(shí)，可以用N1/N的數(shù)值來估計(jì)不等式成立的概率，若N1/N≥β，則機(jī)會約束成立，否則不成立。

3.2 基于隨機(jī)模擬的PSO算法

PSO算法由Eberhart和Kennedy提出，源于對鳥群覓食行為的研究。PSO算法的特點(diǎn)和優(yōu)勢體現(xiàn)在簡單易實(shí)現(xiàn)，能夠記憶個(gè)體最優(yōu)和全局最優(yōu)信息，并沒有太多的參數(shù)需要調(diào)整，非常適用于工程應(yīng)用。

本文將隨機(jī)模擬技術(shù)與粒子群算法結(jié)合，用以求解含機(jī)會約束規(guī)劃的概率最優(yōu)潮流問題，具體算法流程如下:

(1)讀取系統(tǒng)信息，輸入風(fēng)電場容量、風(fēng)電機(jī)組運(yùn)行參數(shù)及用以建立相關(guān)性風(fēng)速模型的Copula函數(shù)的參數(shù);

(2)輸入基于隨機(jī)模擬的粒子群算法參數(shù)、約束條件的置信水平以及迭代的終止條件;

(3)隨機(jī)生成粒子的初始位置，即常規(guī)可調(diào)機(jī)組的有功出力，生成具有相關(guān)性的隨機(jī)風(fēng)速并計(jì)算相應(yīng)的風(fēng)電場出力，計(jì)算系統(tǒng)潮流，利用隨機(jī)模擬技術(shù)檢驗(yàn)粒子的可行性，即不等式約束條件是否達(dá)到指定的置信水平，直至生成滿足粒子規(guī)模n的n個(gè)可行初始粒子群位置及速度;

(4)計(jì)算各粒子的目標(biāo)函數(shù)值;

(5)獲取粒子的個(gè)體極值、個(gè)體最優(yōu)位置以及全局極值、全局最優(yōu)位置;

(6)更新粒子的位置與速度;

(7)計(jì)算系統(tǒng)潮流，利用隨機(jī)模擬技術(shù)檢驗(yàn)更新后粒子的可行性，若不可行，則重新生成粒子的速度繼而更新其位置，直至滿足約束條件，如果重復(fù)更新的次數(shù)超過某一規(guī)定的數(shù)值，則用原來可行的粒子代替;

(8)重復(fù) (4)至 (7)，直至符合終止條件;

(9)輸出最優(yōu)粒子位置及所對應(yīng)的最優(yōu)解。

4 算例分析

本文算例分為兩部分，第一部分以云南某地區(qū)兩個(gè)風(fēng)電場2011年1至12月的實(shí)際歷史風(fēng)速數(shù)據(jù)作為樣本，選取Copula函數(shù)并估計(jì)其相關(guān)參數(shù)，用以生成具有相關(guān)性的風(fēng)速樣本;第二部分將兩風(fēng)電場接入IEEE30標(biāo)準(zhǔn)測試系統(tǒng)，應(yīng)用基于隨機(jī)模擬技術(shù)的粒子群優(yōu)化算法計(jì)算第2節(jié)所述的最優(yōu)潮流問題，并對結(jié)果進(jìn)行分析。

首先對歷史風(fēng)速數(shù)據(jù)進(jìn)行處理，利用最大似然估計(jì)法估計(jì)出兩風(fēng)電場風(fēng)速的威布爾分布參數(shù)如表2所示。

表2 兩風(fēng)電場風(fēng)速的威布爾分布參數(shù)Tab.2 Weibull distribution parameters of the two wind farms wind speed

利用威布爾分布參數(shù)可得到兩個(gè)風(fēng)電場風(fēng)速的邊緣分布，為便于觀察，圖2繪制出了兩個(gè)邊緣分布的二元頻率直方圖。

圖2 邊緣分布的二元頻率直方圖Fig.2 Frequency histogram for two marginal distributions

從圖2可以看出，頻率直方圖的上尾高，下尾低，也就是說兩個(gè)風(fēng)電場風(fēng)速的聯(lián)合密度函數(shù)具有在上尾部比較集中的特點(diǎn)，因此可以選取Gumbel－Copula函數(shù)描述原始數(shù)據(jù)的相關(guān)結(jié)構(gòu)。

上述分析較為直觀，為了得到更客觀的評價(jià)，利用5種Copula函數(shù)分別構(gòu)建兩風(fēng)電場風(fēng)速的聯(lián)合概率分布，通過最大似然估計(jì)法估計(jì)相關(guān)參數(shù)，并利用基于EMC與THC的最短歐氏距離法計(jì)算5種模型的歐氏距離，結(jié)果如表3所示。

表3 兩風(fēng)電場Copula函數(shù)參數(shù)估計(jì)和歐氏距離Tab.3 Parameter estimation and Euclidean distance of Copula function for two wind farms

可以看出，Gumbel－Copula具有最小的歐式距離，能更好地描述兩風(fēng)電場風(fēng)速的相關(guān)性，與利用頻率直方圖得到的結(jié)論是一致的。

下面以IEEE30節(jié)點(diǎn)測試系統(tǒng)［10］作為算例，驗(yàn)證上述基于Gumbel－Copula的相關(guān)性風(fēng)速模型、考慮機(jī)會約束規(guī)劃的最優(yōu)潮流模型以及基于隨機(jī)模擬技術(shù)的粒子群算法的可行性與有效性。該系統(tǒng)有6臺常規(guī)發(fā)電機(jī)組和41條支路，假設(shè)6號機(jī)組不參與優(yōu)化，其發(fā)電功率保持不變，其他機(jī)組參與優(yōu)化，1號機(jī)組為平衡機(jī)組，配合風(fēng)電場的運(yùn)行，可調(diào)機(jī)組對應(yīng)于發(fā)電成本函數(shù)的經(jīng)濟(jì)參數(shù)及出力上下限如表4所示，系統(tǒng)總負(fù)荷為283.4 MW。風(fēng)電場1，2分別接入節(jié)點(diǎn)19和節(jié)點(diǎn)20，額定裝機(jī)容量分別為18 MW和9 MW，風(fēng)電機(jī)組的型號相同，切入風(fēng)速為3.5 m/s，額定風(fēng)速為15 m/s，切出風(fēng)速為25 m/s。

表4 常規(guī)發(fā)電機(jī)組經(jīng)濟(jì)、技術(shù)參數(shù)Tab.4 Economic and technical parameters for conventional units

在置信水平為0.98的條件下，考慮風(fēng)速具有相關(guān)性和風(fēng)速獨(dú)立這兩種情況，分別進(jìn)行10次最優(yōu)潮流計(jì)算，計(jì)算結(jié)果如圖3所示。

圖3 可調(diào)常規(guī)機(jī)組發(fā)電總成本Fig.3 Total fuel cost of adjustable generators

由圖3可以看出，當(dāng)考慮風(fēng)速相關(guān)性時(shí)，可調(diào)常規(guī)機(jī)組的發(fā)電成本均值要比不考慮風(fēng)速相關(guān)性時(shí)略高。分析可知，當(dāng)風(fēng)電場間風(fēng)速相關(guān)性較強(qiáng)時(shí)，風(fēng)速同時(shí)增大與減小的概率較大，導(dǎo)致風(fēng)電場出力波動變大，為了配合風(fēng)電場的運(yùn)行，常規(guī)機(jī)組需要更大幅度地調(diào)節(jié)其出力，從而使發(fā)電成本增加;反之，當(dāng)風(fēng)速相對獨(dú)立時(shí)，風(fēng)電場之間的風(fēng)速呈現(xiàn)出一定的互補(bǔ)性，降低了出力的變化率，較之前一種情況，相對減輕了常規(guī)機(jī)組的調(diào)節(jié)負(fù)擔(dān)，發(fā)電成本也就隨之降低。

在上述模型的基礎(chǔ)上設(shè)置不同機(jī)會約束的置信水平，進(jìn)行多次最優(yōu)潮流計(jì)算后取平均值可得結(jié)果如表5所示。

表5 不同置信水平下的優(yōu)化結(jié)果Tab.5 Optimization results for different confidential levels

由表5可以看出，在相同的置信水平條件下，考慮風(fēng)速相關(guān)性時(shí)常規(guī)機(jī)組發(fā)電成本比獨(dú)立時(shí)大，與上述分析一致。觀察置信水平對結(jié)果的影響，可以看出常規(guī)機(jī)組發(fā)電成本隨置信水平的減小而減小，其原因是允許某些違反約束條件的系統(tǒng)狀態(tài)發(fā)生的概率變大，使最優(yōu)潮流結(jié)果不致過于保守;同理，提高置信水平，系統(tǒng)的發(fā)電成本隨之增大。

5 結(jié)論

本文采用Copula函數(shù)理論構(gòu)建了風(fēng)電場間的相關(guān)性風(fēng)速模型，并生成了相應(yīng)的風(fēng)速樣本空間;將相關(guān)性風(fēng)速應(yīng)用于含風(fēng)電場電力系統(tǒng)概率最優(yōu)潮流模型，以IEEE30節(jié)點(diǎn)測試系統(tǒng)作為算例，通過對計(jì)算結(jié)果的分析，可得出如下結(jié)論:

(1)Copula函數(shù)可以有效地描述風(fēng)速相關(guān)性并生成具有指定相關(guān)系數(shù)的多風(fēng)電場風(fēng)速序列，適用于含風(fēng)電場的電力系統(tǒng)仿真運(yùn)行分析。

(2)風(fēng)速相關(guān)性的大小會影響到系統(tǒng)發(fā)電成本的高低，在制定發(fā)電計(jì)劃或?qū)ο到y(tǒng)進(jìn)行安全性評估時(shí)，應(yīng)考慮風(fēng)速相關(guān)性帶來的影響。

(3)機(jī)會約束規(guī)劃成立的置信水平越高時(shí)，系統(tǒng)發(fā)電成本越高，反之越低，因此在實(shí)際運(yùn)行中，可以根據(jù)經(jīng)驗(yàn)或具體情況調(diào)整置信水平，在保證安全性的前提下使發(fā)電成本最低。

［1］李國清.考慮風(fēng)電綜合利用指標(biāo)約束的輸電規(guī)劃模型［J］.華北電力大學(xué)學(xué)報(bào)，2011，38(2):32－36.

［2］陳金富，陳海焱，段獻(xiàn)忠.含大型風(fēng)電場的電力系統(tǒng)多時(shí)段動態(tài)優(yōu)化潮流 ［J］.中國電機(jī)工程學(xué)報(bào)，2006，26(3):31－35.

［3］肖創(chuàng)英，汪寧渤，陟晶，等.甘肅酒泉風(fēng)電出力特性分析［J］.電力系統(tǒng)自動化，2010，(17):64－67.

［4］楊歡，鄒斌.含相關(guān)性隨機(jī)變量的概率最優(yōu)潮流問題的蒙特卡羅模擬方法 ［J］.電力系統(tǒng)保護(hù)與控制，2012，(19):110－115.

［5］潘雄，周明，孔曉民，等.風(fēng)速相關(guān)性對最優(yōu)潮流的影響 ［J］.電力系統(tǒng)自動化，2013，(6):37－41.

［6］Alicja Lojowska，Dorota Kurowicka，George Papaefthymiou，et al.Stochastic modeling of power demand due to EVs using Copula ［J］.IEEE Trans on Power Systems，2012，27(4):1960－1968.

［7］George Papaefthymiou，Dorota Kurowicka.Using Copulas for modeling stochastic dependence in power system uncertainty analysis ［J］.IEEE Trans on Power Systems，2009，24(1):40 －49.

［8］Cherubini U，Luciano E，Vecchiato W.Copula methods in finance ［M］.England:John Wiley Sons，2004:1－30.

［9］劉寶碇，趙瑞清.隨機(jī)規(guī)劃與模糊規(guī)劃 ［M］.北京:清華大學(xué)出版社，1998.74－94.

［10］Zimmerman R D，Murillo－Sanchez C，Deqing G.Matpower:a MATLAB power system simulation package［EB/OL］.［2011－11－14］.http://www.pserc.cornell.edu/matpower/.