谢应昭，卢继平

（重庆大学 输配电装备及系统安全与新技术国家重点实验室，重庆 400044）

0 引言

随着风力发电逐步向大规模、高集中开发的方向发展，电力系统必然会面临大量可再生能源的接入问题。风力发电的规模性接入，对电力系统运行的各方面都会带来明显影响，电压稳定研究作为系统安全稳定运行的重要部分，将面临更多的不确定因素。传统电压稳定分析都是基于确定性模型，忽略了系统中的随机注入功率的不确定性，以系统最危险工作模式为研究对象进行电压稳定性评估，评估结果过于保守，灵活性不强。为了弥补确定性方法的不足，在电力系统不确定问题研究中引入概率分析方法，这些方法可分为以下3类：蒙特卡洛MC（Monte Carlo）法、解析法以及点估计法。其中，蒙特卡洛法可以得到高精度的随机变量统计信息，但由于计算效率低下，一般仅用作其他概率分析方法的验证标准；解析法通过卷积技术或输出变量的累积量求解概率密度函数，这种方法假定输入和输出随机变量成线性关系，不能反映系统运行的实际情况；点估计法根据已知变量的概率分布求解未知变量的各阶矩信息，计算量小、精度高，是一种理想的概率分析方法。上述3类方法中，只有蒙特卡洛法可直接得出随机变量概率分布，后面2类方法还需要结合其他方法得出随机变量的概率分布。

现有文献主要集中于概率潮流计算问题的研究［1-5］，关于电压稳定问题概率分析的文献较少［6-9］：文献［6］将点估计法应用于电压稳定概率分析的计算，但并未计及风电功率的影响；文献［7］研究了考虑分布式电源的静态电压稳定概率问题，研究过程忽略了输入变量的相关性，用级数展开求取概率分布时，存在一定误差；文献［8］采用随机响应面法对电压稳定概率问题进行求解，取得了较为理想的分析结果；文献［9］对含风电系统进行了电压稳定概率分析，文中采用Nataf变换处理相关输入变量，蒙特卡洛方法求解概率问题，求解过程较为复杂。

本文提出了一种基于多项式正态变换PNT（Polynomial Normal Transformation）和最大熵估计的含风电系统电压稳定概率分析方法。利用多项式正态变换方法处理输入变量之间的相关性，结合点估计法将电压稳定概率分析问题转化为确定性负荷裕度非线性优化模型的求解，并计算负荷裕度的统计特征，采用最大熵估计方法估计输出变量概率分布函数。文中还分析了风速参数变化、相关系数矩阵变化、风电接入容量变化以及风电场功率因数变化对负荷裕度统计特征的影响。本文方法解决了点估计法无法处理输入变量相关性和得出输出变量概率分布的不足，保持了点估计法计算简便、精度高的优点。算例分析表明，本文所提方法具有令人满意的计算精度，且计算效率高，变换方法简洁灵活，能够得到全面反映系统电压稳定性的概率信息。

1 多项式正态变换

概率问题中包含不同概率分布的随机变量，有些变量之间还具有相关性，分布形式的不同和变量的相关性都会对分析结果产生影响。因此有必要对输入变量进行预处理，以减小这些因素带来的影响。多项式正态变换可以同时考虑多个非正态相关的随机变量，利用统一的变换步骤将其变换到独立正态空间，变换方法简洁有效，具有良好的通用性。

1.1 单变量正态变换原理

多项式正态变换是利用标准正态分布随机变量z的多项式对非正态分布的随机变量x进行求解。随机变量x的r阶多项式正态变换表达式为：

其中，a0、a1、a2、…、ar为多项式变换系数。

多项式阶数r可以取不同值，得到不同的多项式变换公式，其中最常用的有二阶、三阶和五阶多项式。尽管高阶多项式可利用高阶矩信息提高多项式正态变换的近似度，但高阶多项式对近似度的改善程度还不明确，而计算过程却更加复杂［10］，因此阶数的确定需要结合具体计算要求进行选择。除了阶数的确定，正态变换的另一关键在于变换系数ai的求取。ai的求解算法主要包括积矩法、L-矩法、最小二乘法以及Fisher-Cornish级数展开法，文献［11］对这些算法做了详细的数值分析和对比，可作为求解算法的选择依据。结合已有文献结论［10-11］，本文选择三阶多项式作为正态变换多项式，并采用L-矩法对变换系数进行求解。

随机变量 x 的 q 阶 L-矩 λq（q=1，2，3，4）与多项式变换系数 ai（i=0，1，2，3）满足如式（2）所示函数关系［11］。

利用前4阶L-矩得出变换系数，即可将随机变量x转换为标准正态变量z进行处理。

本文采用Weibull分布对风速进行刻画，可直接利用Weibull分布参数对前4阶L-矩进行计算［12］，如式（3）所示。

其中，τ3=λ3/λ2和τ4=λ4/λ2分别为 L-偏度和 L-峰度；kw和c分别为Weibull分布的形状参数和尺度参数。通过式（3）即可完成服从Weibull分布的随机变量多项式正态变换。

1.2 多变量多项式正态变换

多随机变量概率分布模型的建立需同时考虑各变量的边际分布和变量间的相关性因素。由1.1节可知，各边际分布随机变量均可采用单变量正态变换方法进行变换处理，本节解决的是把原始随机变量的相关系数矩阵变换为对应的标准正态变量的相关系数矩阵，实现多变量的多项式正态变换。

假设 X=［x1，x2，…，xn］T表示由 n 个随机变量 xi组成的随机向量，对应的相关系数矩阵为Z=［z1，z2，…，zn］T表示由 n 个标准正态变量 zi组成的随机向量，对应的相关系数矩阵为则两相关系数矩阵对应元素满足如下关系［11］：

其中，ami和 amj（m=0，1，2，3）分别为 xi和 xj的多项式变换系数分别为xi和xj的标准差和期望。依次求解RX中各元素对应的等式方程，并选择满足不等式条件的解作为ρzizj的值，即可得到标准正态变量相关系数矩阵RZ。

结合1.1节和1.2节，可以得到由独立标准正态随机变量 S=［s1，s2，…，sn］T求取原始随机变量的完整公式：

其中，Am=［am1，am2，…，amn］T（m=0，1，2，3）为由 1.1 节得出的n个随机变量xi对应的多项式变换系数；L为标准正态变量相关系数矩阵RZ经过Cholesky分解后所得的下三角矩阵。

通过多项式正态变换，为点估计法处理含相关性的随机变量提供数据基础。

2 点估计法

点估计法是求解非线性函数Y=h（X）在随机变量 X=［x1，x2，…，xn］T作用下的统计分布的一种有效方法，具有良好的计算精度和较小的计算复杂度。其原理是首先求解一个r×n的采样点集F1（X）和每个采样点对应的权重系数，然后利用式（6）计算函数变量Y的前2r-1矩信息：

其中，（ μ1，μ2，…，xi，k，…，μn）和 ωi，k分别为随机变量采样点和对应的权重系数，xi，k为X的某一分量估计值， μl（l≠i）为其他分量的自身期望；r和 n分别为随机变量xi的估计点数和变量个数。研究表明当r＞3时，采样点集F1（X）和权重系数集ω可能包含复数解，不符合实际情况，因此r的数值通常取为3。

将随机变量xi的一个估计点选为自身期望μi的三点估计法称为2n+1估计方案。该方案的采样点集 F1（X）中有 n 个相同采样点（ μ1，μ2，…，μi，…，μn），只需进行一次计算即可完成对这n个采样点的处理，具有很高的计算效率。2n+1估计方案下的位置系数 ξi，k和对应权重系数 ωi，k计算公式如下［13］：

其中，i=1，2，…，n；λi，3和 λi，4分别为随机变量 xi的偏度系数和峰度系数；σi和μi分别为xi的标准差和期望。

点估计法适用于输入随机变量相互独立的情况，若输入变量之间存在相关性，可利用第1节的方法进行处理。

3 基于最大熵的概率分布估计

任意随机变量x的熵定义为：

其中，H（x）和 f（x）分别为随机变量 x的熵和概率密度函数。

利用最大熵原理POME（Principle Of Maximum Entropy）估计随机变量的概率密度函数的基本思路为：给定随机变量x的相关统计信息，建立约束条件下的最大熵优化模型，在候选概率分布解集Φ［f（x）］中，选择使H（x）最大的 f（x）作为随机变量 x的概率密度函数最优估计。POME数学模型如下［14-15］：

其中，φ0（x）=1、φn（x）（n=1，2，…，N）为 N+1 个已知函数，称为基函数集为给定的N+1个随机变量x的相关统计信息。POME模型的解可用φn（x）的函数进行表达：

其中，λn（n=0，1，…，N）为拉格朗日参数。求解模型式（8）—（10）得出拉格朗日参数，即可获得概率密度函数 f（x）的最优估计。

本文选择幂函数形式 φn（x）=xn（n=0，1，…，N）作为POME模型的基函数集，则POME模型变为：

4 含风电系统的电压稳定概率分析方法

4.1 风电功率模型

风速的概率分布可用双参数Weibull分布进行描述。风电场风速vi的概率密度函数如式（13）所示：

其中，kw和c分别为Weibull分布的形状参数和尺度参数。

风电机组输出功率与风速的函数关系为：

其中，vci、vco和vR分别为切入风速、切出风速和额定风速；PR为风电机组额定功率。

4.2 负荷方向随机模型

电压稳定分析中，通常选择各节点负荷作为负荷增长方向，以保证各节点负荷功率因数不变且按同比例增加。由于系统测量、估计等方面的误差，实际的节点负荷并不是常数，存在一定的随机性，可选择正态分布反映其随机特点。

假定节点i的有功负荷增长方向IPLi满足以基态有功PLi为均值、以σLi为标准差的正态分布，则其概率密度函数如下［8］：

对应的无功负荷增长方向IQLi由式（16）表示为：

其中，QLi为节点i的基态无功功率。

4.3 负荷裕度非线性优化模型

负荷裕度为运行人员提供了系统从当前运行点到电压崩溃点距离的直观量度，是评估系统电压稳定性的最有效指标。选择负荷裕度为目标函数，同时计及系统等式和不等式约束，构建计算最大负荷裕度的非线性优化模型如下［16-19］：

其中，jєi表示节点j与节点i直接相连，包括j=i的情况；SB为系统节点集合；SG为发电机节点集合；SR为无功源集合；SL为支路集合；PGi和QRi分别为节点i的电源发出的有功和无功功率；PWi和QWi分别为节点i的风电输出有功、无功功率；PLi和QLi分别为节点i的有功负荷和无功负荷；ρ为负荷增长系数，即负荷裕度；αij=δi- δj-φij；Ui和 δi分别为节点 i的电压幅值和相角；Yij和φij分别为节点导纳矩阵对应位置的导纳元素幅值和相角；Sij为节点i与节点j之间的支路潮流大小；下划线变量和上划线变量分别对应各自变量的下限值和上限值。式中等式约束由含负荷裕度参数的扩展潮流方程组成，不等式约束由系统静态安全约束组成。

4.4 含风电系统电压稳定概率分析步骤

在负荷裕度优化模型中加入风电输出功率和节点负荷的随机特性，并采用概率方法进行求解，即可实现含风电系统电压稳定的概率分析评估。

本文假定风电场风速和节点负荷的概率分布参数已知，并给出相应的风速相关系数矩阵，则含风电系统的电压稳定概率分析具体步骤如下。

a.确定输入随机变量个数n，在独立标准正态空间S中对n个标准正态变量进行采样，并求解方程（7），形成基本采样点集F1（S）和对应权重系数集ω。

b.将Weibull分布参数和风速相关系数矩阵代入式（2）—（4），求得相应的变换矩阵L和多项式变换系数 Am=［am1，am2，…，amn］T（m=0，1，2，3）。

c.将基本采样点集 F1（S）代入式（5），得到输入随机变量的采样点集F1（X）。

d.将F1（X）代入负荷裕度非线性优化模型进行求解，得出对应的负荷裕度离散点集F2（ρ）。非线性优化模型采用预测-校正原对偶内点法进行求解。

e.将 ω 和 F2（ρ）代入式（6），计算负荷裕度的各阶矩信息。

f.将负荷裕度各阶矩信息代入模型式（12）进行求解，即可得到负荷裕度概率分布的最优估计。

5 算例分析

5.1 系统概况

本文以IEEE 30和IEEE 118节点系统为基础，对含风电场的电力系统电压稳定进行概率分析，程序实现平台为MATLAB R2009a，所用计算机的CPU主频为2.1 GHz，内存为2 GB。风电场按恒功率因数方式运行，各风电场原始相关参数如表1所示，表中接入台数列“/”前后数字分别表示118节点和30节点的标准接入台数，风电场原始相关系数矩阵RW为：

表1 风电场相关参数Table1 Parameters of wind farms

IEEE测试系统数据参见文献［20］，IEEE 30节点系统中，风电场分别接入节点9、25、28，系统总负荷为189.2 MW；IEEE 118节点系统中，风电场分别接入节点 23、39、114、117，系统总负荷为 3668 MW。

负荷增长方式为全网负荷同时增加，各节点负荷按基态功率因数等比例增长。各负荷分量服从以基态负荷为均值、标准差为5%的正态分布。

5.2 基于最大熵的概率分布估计结果

应用基于最大熵的概率分布估计方法对算例系统进行含风电系统的电压稳定裕度概率计算，得出相应的分布曲线，并同时与40000次蒙特卡洛模拟结果进行对比分析。

图1、2和图3、4分别给出了IEEE 30和 118节点系统在最大熵估计和蒙特卡洛模拟2种方法下对应的负荷裕度概率密度分布曲线和累计概率分布曲线，图中负荷裕度为标幺值，后同。从图中可以看出，由最大熵估计方法得出的分布曲线较为合理地拟合了蒙特卡洛仿真结果。对比图1和图3中由最大熵估计方法得出的概率密度曲线可知，图1中的密度曲线发生了畸变，存在“翘尾现象”，而图3中的密度曲线变形较小，仅在其尾部有一定程度的收缩。结合5.1节的算例数据可知，IEEE 30节点系统的风电接入比例大于IEEE 118节点系统，这表明非正态变量比例的增加会对最大熵概率分布估计结果产生影响。观察图2和图4可知，2种方法所得累计概率分布曲线较为接近，较好地反映了负荷裕度的累积分布特性。

表2给出了服从Weibull分布的风电场风速随机变量多项式正态变换系数，相应的变换矩阵L为：

表3给出了基于最大熵的概率分布估计参数，可用对应的估计函数描述负荷裕度的概率密度分布。表4给出了2种方法下不同算例系统的负荷裕度均值（标幺值，后同）、标准差（标幺值，后同）以及2种方法计算结果的相对误差。由表4数据可得，均值的相对误差都小于1%，标准差的相对误差都小于4%，计算精度令人满意，在工程应用的误差要求范围之内。表5给出了2种方法下的计算仿真规模和计算时间。由表5可得，最大熵估计方法的计算时间仅为蒙特卡洛方法的0.14%和0.6%，大幅节省了计算时间且能获得较为满意的计算结果，为电压稳定概率分析方法的工程实际应用实现提供了可能性。

图1 IEEE 30节点系统负荷裕度概率密度分布Fig.1 Probabilistic density distribution of load margin of IEEE 30-bus system

图2 IEEE 30节点系统负荷裕度累计概率分布Fig.2 Cumulative probabilistic density distribution of load margin of IEEE 30-bus system

图3 IEEE 118节点系统负荷裕度概率密度分布Fig.3 Probabilistic density distribution of load margin of IEEE 118-bus system

图4 IEEE 118节点系统负荷裕度累计概率分布Fig.4 Cumulative probabilistic density distribution of load margin of IEEE 118-bus system

5.3 风电参数影响分析

为了研究风速概率密度函数变化和风速相关系数矩阵变化对电压稳定裕度的影响，本节假设在不同风速分布参数和相关系数矩阵组合下对电压稳定问题进行概率求解和结果分析，得出风电参数变化对电压稳定裕度概率分布所带来的影响。

表2 风速随机变量多项式正态变换系数Table 2 Polynomial normal transformation coefficients of wind speed random variables

表3 基于最大熵原理的概率分布估计参数Table 3 Parameters of probability distribution estimation based on maximum entropy

表4 两算例系统计算结果及误差分析Table 4 Calculative results and error analysis for two case systems

表5 两算例系统仿真规模和计算时间Table 5 Simulation scale and calculation time of two case systems

本节以IEEE 118节点系统为研究对象，分别以表1的风电场原始风速参数和原始相关系数矩阵RW为标准，同方向增加或减小风速分布参数数值和相关系数矩阵中的元素数值，分别形成以下3个新风速参数场景和3个新相关系数矩阵场景：风速参数场景取名为“较小风速场景”、“较大风速场景1”和“较大风速场景2”，各场景的具体参数如表6所示；相关系数矩阵场景取名为“零相关系数矩阵”、“较小相关系数矩阵”和“较大相关系数矩阵”，其中“零相关系数矩阵”用4阶单位矩阵I4×4描述，“较小相关系数矩阵”和“较大相关系数矩阵”分别用RWS和RWL描述，各系数矩阵参数如下：

表6 不同场景下的风速分布参数Table 6 Parameters of wind speed distribution for different cases

图5为不同风速参数场景下计算得出的负荷裕度概率密度曲线。观察图中曲线，可以看出风速参数对概率密度分布的影响：随着风速参数的逐渐增大，概率密度曲线朝着扁平化的趋势发展，曲线顶点对应的概率密度数值逐渐减小，负荷裕度数值逐渐增大。这说明风速参数的增加有助于系统负荷裕度的提高，但随着风速参数的增加，负荷裕度随机变量的变化也逐渐变大，风速变化所带来的波动性更加明显。

图5 不同风速下的负荷裕度概率密度曲线Fig.5 Probabilistic density curve of load margin for different wind speeds

表7所示为不同风速参数和风速相关系数矩阵组合下的负荷裕度计算结果，通过对比表中数据，验证了图5所得相关结论的正确性。由上述分析可知，在选择风能资源时，并不能只考虑风资源本身的变化因素，还要考虑所接入系统的电压稳定承受能力等影响因素，这样才能将风能资源和所接系统进行合理的匹配。

表7中还给出风速相关系数矩阵变化时的负荷裕度模型计算结果。观察表中数据可知，相关系数变化会对负荷裕度计算结果产生影响：随着风速相关性增加，有些负荷裕度均值逐渐减小，有些负荷裕度均值先减后增，相关系数的变化引起负荷裕度均值的变化；随着相关系数增大，负荷裕度标准差逐渐增大，表明其波动更加明显。因此在负荷裕度概率分析时，不能忽略风速相关性，必须加以考虑，否则分析结果与实际情况偏差较大，不能正确反映负荷裕度变化的概率特性。本文方法所得结果和蒙特卡洛模拟所得结果的相对误差如表7所示。从表中可以看出，采用本文方法所得的负荷裕度结果与蒙特卡洛模拟所得结果较为接近，结果均值误差较小，标准差误差偏大，但均在可接受范围之内，可认为本文方法所得结果正确有效。综上可知，采用基于多项式正态变换的点估计法能较好地处理风速相关性问题，计算结果反映了负荷裕度变化趋势，相对蒙特卡洛模拟方法，计算时间缩短为原来的1%以下，适合于对电力系统概率问题的分析求解。

5.4 风电场接入容量影响分析

随着风电场的发展，风电容量的接入比例会越来越大，本节就不同的风电接入比例对负荷裕度概率模型进行求解，观察容量变化带来的影响。标准风电接入容量为230MW，考虑不同的标准接入容量倍数，得到不同风电接入比例下的负荷裕度均值和标准差，具体数据如表8所示。可以看出，风电比例越高，负荷裕度的均值和标准差均有所增加，这是因为一方面风电的接入提供了电源支持，使得负荷裕度水平有所提高，但是由于风电固有的随机特性，容量的增加也让波动更加明显，标准差随之增大。因此，如何有效地减小风电波动性带来的影响是大规模利用风电必须考虑的问题，否则大量的风电接入会给系统运行的稳定性造成不利的影响。

表7 不同风电参数组合下的负荷裕度模型计算结果Table 7 Calculative results of load margin for different wind power parameter combinations

表8 不同风电接入容量下的负荷裕度模型计算结果Table 8 Calculative results of load margin for different wind power penetration capacities

5.5 风电场功率因数影响分析

风电场发电控制通常采用恒功率策略，不同的功率因数设定值对负荷裕度有所影响。本节采用最大熵估计法快速求解不同功率因数下的负荷裕度模型，计算结果如表9所示。表中负功率因数表示风电场需要消耗无功。

表9 不同功率因数下的负荷裕度计算结果Table 9 Calculative results of load margin for different power factors

由表9可以看出，随着功率因数的增加，负荷裕度均值逐渐增加，标准差逐渐减小。这是因为随着功率因数的增加，整个风电场的无功消耗逐渐降低，无功输出逐渐增加。由于风电场能够向系统提供无功，因此系统的电压稳定情况得以改善，电压稳定裕度有所增加，能够承受更大的系统变化，在同样的有功波动情况下，能使电压稳定裕度的均值更大，而标准差更小，电压稳定特性更加优良。表中还给出了不同相关系数矩阵条件下，功率因数变化时的负荷裕度计算结果，不同条件下的计算结果变化趋势一致，验证了前述的原因分析。

6 结论

本文提出了一种基于多项式正态变换和最大熵估计的含风电系统电压稳定概率分析方法。该方法结合多项式正态变换、点估计法以及最大熵估计，能够对考虑相关性的含风电系统电压稳定问题进行完整的概率分析，得出全面的概率信息。本文将风电输入功率和节点负荷功率视为随机变量，以蒙特卡洛模拟作为评价标准，验证了本文算法的计算精度，讨论了不同风速参数、相关系数矩阵变化、不同风电接入容量以及风电场功率因数变化对电压稳定裕度的影响。算例分析结果表明，基于多项式正态变换和最大熵估计的电压稳定概率分析方法能够获得正确的负荷裕度概率分布，且与蒙特卡洛模拟相比，具有更高的计算效率，能够满足负荷裕度概率信息及时获取的要求，为运行人员控制决策提供了重要参考。