第46卷第21期电力系统保护与控制Vol.46 No.21 2018年11月1日Power System Protection and Control Nov. 1, 2018 DOI: 10.7667/PSPC171534
基于Copula函数及Ronblatt变换的含相关性
概率潮流计算
王 涛，王 淳，李成豪
(南昌大学电气与自动化系，江西 南昌 330031)
摘要：随着风电渗透率的增加，在考虑风电出力随机性和间歇性对配电网影响的同时，风电场间风速的相关性应予以考虑。将Ronblatt变换与Copula函数相结合以处理风电场间复杂的非线性关系，提出一种基于Copula函数及Ronblatt变换的概率潮流半不变量算法。所提方法能准确捕捉到输入变量间的非线性相关关系，且具备半不变量法计算效率高的优点。以IEEE33节点网络接入风电场对所提方法进行测试。结果证明了所提方法的有效性、准确性和实用性，所提方法计算结果较考虑线性相关性算法的计算结果更接近实际运行情况。
关键词：概率潮流；半不变量法；Copula函数；Ronblatt变换；相关性
Probabilistic load flow calculation bad on Copula function and Ronblatt transformation
considering correlation among input variables
WANG Tao, WANG Chun, LI Chenghao
(Department of Electrical and Automatic Engineering, Nanchang University, Nanchang 330031, China)
Abstract: With the increa of wind power penetration, the correlation of wind speed among wind farms should be considered while considering the randomness and intermittency of wind power.In this paper, Ronblatt transformation and Copula function are combined to deal with the complex nonlinear relationship among wind farms, and a probabilistic load flow cumulant method bad on Copula function and Ronblatt transformation is propod. The method can capture the nonlinear correlation among input variables accurately and has the advantage of high efficiency of cumulant method.
The IEEE33-bus system including wind farms is ud for testing the propod method. The results verify the validity, accuracy and practicability of the propod method, and the calculation results of the propod method are clor to the actual operating condition than that of the method which considers linear correlation.
This work is supported by National Natural Science Foundation of China (No.51467012).
Key words: probabilistic load flow; cumulant method; Copula function; Ronblatt transformation; correlation
0 引言
实际电网中存在诸多不确定因素[1-2]，如分布式电源出力的随机性和间歇性以及负荷的波动。概率潮流计算是分析不确定性对电力系统影响的重要方法，目前概率潮流的计算方法主要有蒙特卡罗模拟法[3]、半不变量解析法[4]和点估计法[5]。蒙特卡罗模拟法能获得很好的精度，但计算耗时，常用于对其他计算方法的校验。半不变量解析法将输入变量与输出变量的关系线性化，求解状态变量的概率分布。
基金项目：国家自然科学基金项目资助(51467012) 点估计法采用近似公式求取状态变量的统计特征。
在电网运行中，还存在大量相关性因素，如邻近的分布式电源之间往往具有相关性，相关性因素对系统的运行产生重要的影响[6]。因此不仅需要考虑分布式电源出力的随机性和间歇性，还需考虑其相关性。文献[7]提出基于Nataf变换的蒙特卡罗模拟法，分析了分布式电源和负荷相关性对电网可靠性的影响。文献[8]提出一种基于Cholesky分解的半不变量算法，分析了不同相关系数对系统运行特性的影响。
文献[9]采用奇异值分解代替Cholesky分解，提出一种考虑相关性的扩展准蒙特卡罗随机潮流算法。以上变换均是基于Pearson线性系数，无
王 涛，等 基于Copula 函数及Ronblatt 变换的含相关性概率潮流计算 - 19 -
法捕捉到变量间的非线性关系。近年来许多学者将广泛应用于金融和统计分析中的Copula 函数用于分布式电源间相关性的研究[10-14]，Copula 函数不仅能描述变量间的线性相关性关系，还能捕捉到变量间的非线性相关性。基于Copula 函数描述的变量包含相关性，已有文献多采用蒙特卡罗模拟法对其进行概率潮流计算，计算时效低。研究一种既能准确描述变量间的相关性，且求解快速的概率潮流计算方法非常必要。基于此，本文提出一种基于Copula 函数及Ronblatt 变换的含相关性半不变量概率潮流算法。
1 风电场风速概率分布
建立风速的概率模型主要有参数法和非参数法：参数法基于假定的概率分布函数，通过估计其参数，进而得到概率分布；非参数估计法不需要分布的先验知识和任何概率形式的假设，是一种从数据样本本身出发，研究数据分布特征的方法。不同地区的风速通常呈现不同的概率分布特性，参数估计模型只对某一地区的风速有较高拟合度。非参数核密度估计法已引起学者的关注并在电力系统领域得到一定的应用[15-17]，风速的非参数核密度估计模型较参数模型具有更高的准确度和通用性[18]。
本文采用非参数核密度估计法建立风速的分布模型，其具体求法参见文献[19]。窗宽和核函数是
核密度估计的重要参数，本文选取高斯核函数，即
21
())2K x x =- (1)
窗宽H 的求取参考文献[20]。
0.21.06H T σ-≈ (2)
式中：σ为样本的标准差；T 为样本容量。
2 基于Copula 函数含相关性的风速联合分
布模型
2.1 Copula 理论
假设n 个随机变量X 1、X 2、 、X n 的边缘分布函数分别为F 1(x 1)、F 2(x 2)、 、F n (x n )，其联合分布函数为H (X 1、X 2、 、X n )，Copula 函数是描述边缘分布和联合分布的连接函数。即有： 121122()(()()())n n n H X X X C F x F x F x = 、、、、、、 (3)
式中，C 为Copula 连接函数。
Copula 函数构造联合模型对各变量的边缘分布没有特定的要求，且相比其他构造联合模型的方式具有更高精度。
Copula 函数的种类繁多，考虑其计算精度和计算难度，当前应用最为广泛的Copula 函数共两大类：Archimedes Copula 函数(包含Frank 、Gumbel 和Clayton Copula 函数等)和椭圆Copula 函数(包含正态和t-Copula 函数等)。其中Archimedes Copula 函数在收敛率和收敛速度上普遍优于椭圆Copula 函数，因此本文选择表1的三类Archimedes Copula 函数进行风速联合模型的构造。
表1 阿基米德Copula 函数 Table 1 Archimedes Copula functions
函数名称 Copula 函数 参数
Frank
1
1
(e 1)
1
ln(1)(e
1)
n
u i n C θ
θ
θ=---=-+
-∏
0θ≠
Gumbel 1/1
exp(((ln )))n
i i C u θθ==--∑
1θ≥ Clayton
1/1
周星驰金典台词
(1(1))n
i i C u θθ--==+-∑
1,0θθ≥-≠
2.2 最优Copula 函数的选择
在选择最优Copula 函数前，需要基于初始样本确定每个备选Copula 函数的参数θ，对于θ 的求取方法，有极大似然估计、分布估计、半参数估计和非参数估计等方法，本文采取应用较为广泛的极大似然估计法。
对于Copula 函数的求取，本文采用基于核密度估计的最短距离法，用式(4)所示的基于乘积核的核密度(,)F x y ∧
作为真实分布函数的估计值，按式(5)计算每个备选Copula 函数和(,)F x y ∧
的距离d k ，选择d k 最小的Copula 函数为最优Copula 。
11(,)T t t t X Y x x y y F x y T h h ΦΦ∧=⎛⎫⎛⎫--=⋅ ⎪ ⎪
⎝⎭⎝⎭∑ (4) 12
妈妈梦见自己孩子丢了21()T k k t d F C ∧=⎧⎫⎪⎪=-⎨⎬
⎪⎪⎭⎩
∑
(5) 式中，C k 为第k 个备选Copula 函数。
3 基于Ronblatt 的风速独立化转换
第2.2节建立的风速联合分布模型包含相关性，而半不变量法要求输入变量相互独立，因此需将含相关性的风速独立化，常用的方法是Rackwitz- Fiessler 逆变换、Nataf 逆变换和Ronblatt 逆变换。Rackwitz-Fiesslerb 逆变换假设变换过程中变量间相关性不变，或者忽略相关性的变化，将导致较大误差；Nataf 逆变换利用变量的边缘概率密度及相关系数矩阵进行变换，相关系数矩阵无法捕捉到变量间
- 20 - 电力系统保护与控制
的非线性关系，当很多样本分布类型或者变量间的联合分布函数不服从高斯分布时，该方法可能存在较大误差。在已知变量联合分布的前提下，Ronblatt 逆变换采用变量的条件概率密度和边缘概率密度进行变换，变换过程中条件概率密度可准确反映变量间的相关性，因此不受变量分布类型以及相关性是否线性的影响，是一种精确的变换方法[21]。本文采用Ronblatt 逆变换将考虑相关性的风速独立化。
若一组随机变量T 12(,,,)d X X X = X 具有联合概率密度函数12(,,,)d f x x x ，则该联合概率密度函数可以表示为一系列的条件概率密度函数的乘积，即
121121121|21|121(,,,)()(|) (|,,,)d d d X X X X X X X d d f x x x f x f x x f x x x x --=
(6)
式中，121|121(|,,,)i d X X X X i i f x x x x -- ，(2,,)i d = 为随机变量X i 的条件概率密度函数。
对于二维变量，由式(6)可得变量X 2在X 1下的条件概率密度如式(7)。
2111212112212|121121212112122(,)
电大专业
(|)=
()((),())()()
大学生性行为
=()
((),())()
X X X X X X X X X X X f x x f x x f x c F x F x f x f x f x c F x F x f x =⋅⋅
(7)
式中：c 12为变量X 1、X 2的Copula 密度函数；1X F 、2X F 和1X f 、2X f 分别为随机变量X 1、X 2的累积分布
函数和概率密度函数。
因此随机变量X 2的条件累积分布函数为 2121122|21|212121222
不鸣则已
(|)(|)d =
((),())()d x X X X X x
X X X F x x f x x x c F x F x f x x -∞-∞
=⋅⎰⎰
(8)
为方便表述，此处引入h 方程表示二维条件分布，如式(9)所示。
聚美优品创始人|(,)(|)i j i j i j X X i j h r r F x x =
(9)
式中：()i i X i r F x =；()j j X j r F x =。
Ronblatt 变换直接将一组相关非正态变量T 12(,,,)d X X X = X 转化成独立的标准正态变量T 12(,,,)d U U U = U ，根据等概率边缘变换原则有
12112112|21|121()()()(|)
()(|,,,)d d X X X d X X X X d d u F x u F x x u F
x x x x ΦΦΦ-=⎧⎪
=⎪⎨
⎪
⎪=⎩
(10) 则可得独立标准正态变量T 12(,,,)d U U U = U ，表示为
12112111112|211
|121()(|)(|,,,)d d X X X d X X X X d d u F x u F x x u F x x x x ΦΦΦ-----⎧⎡⎤
=⎣⎦⎪
⎪⎡⎤=⎪⎣⎦⎨
⎪
⎪⎡⎤=⎪⎣⎦⎩
(11) 式中，Φ为标准正态分布变量样本的累积分布函数。
以上是Ronblatt 变换正过程，其逆过程可求得独立的变量。在通过Copula 函数求得变量边缘分布和条件分布的基础上，由Ronblatt 逆变换生成输入变量X 1、X 2的独立样本，其基本步骤如下。
1) 通过蒙特卡罗产生服从独立标准正态分布变量U 1、U 2的n 个样本。
2) 根据式(10)，有
1111(())X x F Φu -= (12) 将U 1的n 个样本代入式(12)求得变量X 1的n 个样本。
3) 令111()X r F x =，222()X r F x =，再结合式(9)和式(10)，可得
121212(,())r h r u Φ-= (13)
将U 2和X 1的n 个样本代入式(13)求得变量r 2
的n 个样本。
4) 将变量r 2的n 个样本代入式(14)即可求得变量X 2的n 个样本。
2122()X x F r -= (14)
对独立化后的风速样本，按式(15)所示的风电
场风速与有功出力模型求取风电场有功出力，再根据恒定的功率因数求得无功出力。
ci 12ci r
r r co co 0 () 0 v v k v k v v v P v P v v v v v <⎧⎪
+<<⎪⎨
<<⎪⎪>⎩
= (15) 式中：P r 为风机额定出力；v ci 、v r 和v co 分别为切入风速、额定风速和切出风速；1r r ci /()k P v v =-；21ci k k v =-。
4 基于Copula 函数及Ronblatt 变换的半
不变量法概率潮流计算
基于Copula 函数及Ronblatt 变换的半不变量法的核心在于根据输入变量的初始样本构建Copula 联合分布模型并求得条件分布函数及产生满足相关条件的独立样本，再以简化为线性的交流
王涛，等基于Copula函数及Ronblatt变换的含相关性概率潮流计算- 21 -
潮流方程为基础，利用半不变量的齐次性和可加性，再根据线性关系求取输出变量的半不变量，最后通过Gram-Charlier级数展开的方法计算其概率分布。计算步骤如下：
1) 根据风速的历史测量数据采用非参数核密度计算每个风电场风速的分布函数。
2) 根据风速间的相互关系，按2.2节的方法选择合适的Copula函数。
3) 根据所求Copula函数采用蒙特卡罗抽样技术产生含相关性的风速样本。
4) 利用逆Ronblatt变换将含相关性的风速样本转换成独立风速样本。
5) 对独立的风速样本利用概率潮流半不变量法计算输出变量的概率分布。
5 算例分析
5.1 算例简介
算例选取IEEE33 节点配电网络，负荷波动按正态分布，其原确定性数据作为期望μ，负荷的标准差σ取期望的10%。在网络的末端节点18、33分别接入额定功率为0.5 MW的风机。取根节点电压标幺值为1 p.u.。风速取自加拿大魁北克和三河两市实测数据[22]，取恒定功率因数为0.98，风机的额定风速为15 m/s，切入风速为3 m/s，切出风速为25 m/s。
仿真平台采用Intel四核CPU，频率为3.60 GHz，内存为16.0 GB。根据算例参数及本文算法流程，用Matlab R2014a编制相应程序，进行考虑风电场相关性的概率潮流计算，并以基于Copula函数的蒙特卡罗仿真法所得结果作为参照，比较本文所提算法的计算精度。
5.2所提方法性能评估
为验证所提算法的有效性和准确性，本文进行了大量的蒙特卡罗模拟仿真。仿真发现对该算例抽样10 000次、20 000次及50 000次蒙特卡罗模拟的结果非常接近，抽样10 000次相对于抽样20 000次蒙特卡罗模拟输出变量期望的最大变化率为0.120%，相对于抽样50 000次计算结果的最大变化率为0.411%。考虑计算时间及计算精度，以10 000次蒙特卡罗模拟计算结果作为参考值，比较所提方法的有效性和准确性。
由表2中统计误差可知，所提算法和蒙特卡罗法计算结果基本一致，说明所提算法在输出变量统计数字特征方面具有较高的准确度。
为了更好地验证本文所提方法在输出变量概率分布特性方面的计算精度，此处采用文献[23]定义的方差和的根均值(ARMS)来描述两种方法的计算误差，其定义如下。
表2 节点电压幅值和支路功率误差统计
Table 2 Analysis results of bus voltage magnitude and branch flow
电压幅值/p.u. 有功功率/MW 无功功率/Mvar 类型期望
(10-2)
标准差
驾校广告宣传语
(10-2)
期望
(10-2)
标准差
(10-2)
期望
(10-2)
标准差
(10-2)
平均误差0.007 0.089 0.006 0.138 0.152 0.026
最大误差0.020 0.262 0.027 0.408 0.938 0.127
RMS
A N
γ=(16)
式中：
i
开放社会C
M和
i
R
C分别为用蒙特卡罗法和本文所提算
法求得输出变量在统计点i处的累积概率值；γ表示输出变量的类型，包括节点电压幅值U、支路有功功
率P和无功功率Q；
RMS
Aγ为误差指标；N为取点总数。
ARMS计算结果见表3。输出变量的ARMS指标平均值的最大值为0.000 638，ARMS指标最大值的最
大值为0.000 905，说明所提算法能够准确地求得各输出变量的概率分布特性。以离风机接入点较近的节点16、支路16-17和离风机接入点较远的节点20、支路20-21为例，其ARMS指标分别为0.000 853、0.000 25和0.000 887、0.000 166，图1为累积分布曲线。虽然节点16以及支路16-17离风机接入点较近，受风电随机性和相关性影响较为严重，所提方法仍然能够准确求得其概率分布，进一步证明了所提算法具有较高的精度。
表3 输出变量ARMA统计结果
Table 3 Result of output variable ARMA
输出变量平均值/10-2最大值/10-2
电压幅值/p.u.0.063 80.090 5
有功功率/p.u.0.032 80.087 2
无功功率/p.u.0.019 90.055 5
图1 输出变量累积分布曲线
Fig. 1 Cumulative distribution curve of output variables
- 22 - 电力系统保护与控制
为了进一步验证所提方法的优势，此处将所提方法与传统半不变量法、基于Nataf 的半不变量法进行比较，计算系统节点电压幅值、支路有功的期望、标准差如图2所示，节点电压幅值标准差和支路有功的标准差如图3所示。图中ca1采用传统半不变量法，不考虑风速相关性；ca2采用Nataf 变换的半不变量法，考虑风速为线性相关性；ca3采用本文方法，考虑风速为非线性相关性。
图2 节点电压幅值、支路有功的期望
Fig. 2 Expectation of voltage amplitude and active power
图3 节点电压幅值、支路有功的标准差
Fig. 3 Standard deviation of voltage amplitude and active power
从图2和图3可以看到：三种情形下节点电压幅值、支路功率的期望变化较小，但标准差有明显的区别。以标准差的平均值为指标，分析相关性对
系统运行状态的影响[20]。ca2 较ca1系统电压幅值标准差的平均值增加4.22%，支路功率标准差
的平均值增加10.21%；ca3较ca1系统电压幅值标准差的平均值增加16.98%，支路功率标准差的平均值增加15.77%。说明风速相关性对系统基准运行点的影响很小，但会增加系统的波动幅度，使网络运行安全性降低，这与文献[24]所得结论吻合，证明了所提方法的合理性，考虑非线性相关性后波动幅度进一步增大，说明以线性关系不能准确描述风速相关性对系统的影响，所提方法能更精确地捕捉到这种影响，所提算法在输入变量相关性的描述上优于其他解析法。从表4可以看出所提算法在计算速度上有明显优势。
表4 不同算法计算时间比较
Table 4 Calculation time of different algorithms
计算 方法 基于Copula 函数 的蒙特卡罗法
传统半不 变量法 基于Nataf 的半 不变量法 所提 方法 计算时间/s
62.13
0.85
0.93
0.89
以上算例表明所提方法兼具Copula 函数和半不变量法的优点，在充分考虑风电场非线性相关性的同时，获得半不变量法高效的计算速度。 6 结论
本文将Copula 函数与Ronblatt 变换相结合，并引入到半不变量概率潮流计算中。提出一种可准确处理输入变量间相关性的概率潮流半不变量法。Copula 函数能捕捉到风速间的线性相关性和非线性相关性；Ronblatt 变换使其适用于要求输入变量相互独立的概率潮流半不变量法，较其他变换方法具有过程简单、变换精确的特点。
以IEEE33节点接入风机的配电网络进行测试，结论如下：
1) 算例分析验证了所提方法的有效性、准确性和实用性。
2) 所提方法兼具copula 函数与半不变量解析法的优点，在准确捕捉输入变量相关性的同时，具有较高的计算效率。
3) 风电场间的相关性会增加系统的波动，所提方法较其他考虑线性相关性的方法能更准确地捕捉到这种波动。 参考文献
[1] 张玉, 莫寒, 张烈平. 基于模糊支持向量机的光伏发
电量预测[J]. 热力发电, 2017, 46(1): 116-120.
ZHANG Yu, MO Han, ZHANG Lieping. Photovoltaic power prediction bad on fuzzy support vector machine
