基于半不变量及最大熵的概率谐波潮流算法

曾　江， 蔡东阳， 黄德华

(华南理工大学电力学院， 广东省广州市 510641)

0　引言

随着大量非线性负载的接入，电网中的谐波污染日益严重，对电网中谐波传播与分布的研究工作具有重要意义，谐波潮流计算是研究谐波分布特性的有效手段[1-3]。传统谐波潮流计算旨在求解特定工况下电网各节点基波、谐波电压以及支路潮流的确定值。但实际电网中存在着诸多不确定因素，如负荷变化、系统随机故障、网络结构变化等。此外，新能源发电技术的快速发展也加剧了电力系统的不确定性[4-5]，确定性谐波潮流已不能满足工程应用的需求。因此，考虑系统的不确定性，进行概率谐波潮流计算[6-7]，对电网谐波管理、薄弱环节分析具有重要意义。

目前，概率谐波潮流算法主要有解析法、模拟法和拟合法三种。解析法[8-9]指构建待求随机变量的概率密度解析式进行求解，概率谐波潮流的求解实质上是多个随机变量的叠加，从理论上可以用卷积法进行求解，但卷积涉及大量积分计算，当随机变量数量增多时，运算将变得十分困难，因此不适用于大规模网络的概率谐波潮流计算。模拟法[10-11]主要包含蒙特卡洛法，其思路是根据网络注入量的概率密度函数进行反复抽样计算，得到多组谐波潮流的结果，最后进行统计分析得到节点状态量的概率密度。为了保证计算结果的精度，可能需要进行上万甚至几十万次的运算，对运算时间和内存空间消耗巨大，且难以进行灵敏度分析。但由于蒙特卡洛法具有编程简单、准确度可靠等优点，通常作为对照实验方法使用。拟合法包括点估计法[12]、仿射法[13]、级数拟合法[14]等，该方法不直接求出随机变量的概率密度，而是通过求取随机变量的数字特征，通过数值拟合方法，获得概率密度函数，具有运算量小、精度可靠等特点。其中仿射法常用于改善区间分析的保守性，能够快速计算获得谐波的分布区间，但无法分析区间内的概率分布；点估计法通常需要构造2m或2m+1个场景(m为随机变量个数)进行计算，当随机变量数量增多时，计算量会明显增大，有一定的局限性。

本文提出一种基于半不变量及最大熵的概率谐波潮流算法，根据谐波电流的监测样本数据计算电网谐波电压的概率密度函数，属于拟合法的范畴。在计算多个随机谐波电流共同作用产生的谐波电压概率密度函数时，为了避免卷积运算，引入半不变量，半不变量能够与样本数据的高阶矩互相转换，易于求取，利用半不变量的可加性与齐次性简化概率求解过程，且计算量不会随着随机变量的增多而明显增加。根据最大熵原理能够在给定的约束条件下确定一种主观成分最少、最接近实际情况的概率密度。以待求谐波的高阶矩作为约束条件构建最大熵模型，可以快速准确地拟合出谐波概率密度函数，进而求取谐波95%概率值等特征参数。将本文算法在4节点系统上与解析法比较，以及在IEEE 57节点系统上与蒙特卡洛法比较，验证本文算法的有效性。

1　谐波潮流方程及数学基础

1.1　谐波潮流方程

在谐波潮流计算中，通常将谐波源等值为恒流源模型或诺顿等效模型[15]，本文采用恒流源模型进行等值。假设谐波潮流计算中的谐波次数为h，计算电网各元件谐波参数，生成谐波导纳矩阵Yh。根据谐波注入电流情况和谐波导纳矩阵，可以列出谐波潮流方程：

Ih=YhUh

(1)

式中：Ih为h次谐波注入电流；Yh为h次谐波导纳矩阵；Uh为节点的h次谐波电压。

当谐波源注入电流和网络参数已知时，可运用稀疏技术对式(1)进行求解，得到电网各节点谐波电压。若对谐波导纳矩阵求逆可得谐波阻抗矩阵，相应地，式(1)可改写为：

(2)

式中：Zh为h次谐波阻抗矩阵。

使用计算机求解谐波潮流时，采用大规模电网的计算框架，谐波导纳矩阵可根据电网拓扑结构和基波参数直观地求得，采用稀疏矩阵技术储存数据和求解谐波潮流方程，有利于节省内存和提高计算速度，与基波潮流的成熟算法类似，以上方法适用于各种规模的电网计算。

1.2　半不变量

当考虑谐波源注入电流的随机性时，式(2)则无法满足计算要求。此时将谐波电流视为随机变量进行概率谐波潮流计算，并引入半不变量简化计算过程。

半不变量是随机变量很重要的一种数字特征，又称累积量，半不变量序列可唯一确定随机变量的分布规律[16]。根据定义，半不变量可通过对随机变量的特征函数取自然对数，再进行泰勒展开求取而得[17]，计算过程复杂，在实际运用中，通常将半不变量与高阶原点矩进行互相转换。对随机变量X，其半不变量与高阶原点矩的关系如下：

(3)

式中：gk为k阶半不变量；ak为k阶原点矩。

利用随机变量X的样本数据可直接求得ak，有

(4)

式中：xi为随机变量X的第i个可能取值；pi表示取值为xi时的概率。半不变量具有可加性和齐次性两个重要性质，在随机分析的简化计算中起到关键作用。

可加性：若X(t)为n个互相独立的随机变量之和，即X(t)=X(1)+X(2)+…+X(n)，则X(t)的k阶半不变量为：

(5)

齐次性：若随机变量X(t)与随机变量X(1)呈线性关系X(t)=aX(1)+b，则X(t)的k阶半不变量为：

(6)

可以看出，半不变量与高阶矩之间的转换、随机变量的半不变量之间的转换均可轻易实现，因此半不变量在随机变量分析中具有独特优势，能够大大简化运算过程。

1.3　最大熵原理

本文采用最大熵原理，实现高阶矩对概率密度函数的拟合。最大熵又称最大信息熵，其拟合随机变量概率密度的中心思想是：在满足给定的约束条件下，使得信息熵取最大值所对应的概率密度分布是最接近实际情况的分布。自提出以来，最大信息熵原理已在通信、交通、气象等领域获得成功应用，在电力系统中也越来越受到关注[18-20]。

假设X是一个离散型随机变量，取值为xi时对应的概率为pX(xi)，则X的信息熵为：

(7)

根据最大熵原理的思想，当信息熵取得最大值时，得到的随机变量概率密度函数最接近事实，含主观误差最小。以随机变量X的k阶原点矩ak和归一化条件作为约束，对信息熵建立最大值规划模型：

(8)

(9)

(10)

谐波电流、电压属于连续型随机变量，但谐波监测样本数据是一组离散数列，因此对谐波概率密度函数的拟合也转换为对离散值概率分布的拟合。

2　基于半不变量及最大熵的概率谐波潮流

2.1　谐波潮流方程线性化

省略谐波次数h，将式(2)的谐波潮流方程用相量的形式表示：

(11)

假设电网结构不发生改变，即上式中Zij不变，谐波源谐波电流Ij为节点注入量，节点谐波电压Ui是节点状态量，可以写出方程：

U=f(I)

(12)

假设谐波电流的监测样本中电流幅值数据为{I}，选取幅值期望值作为谐波潮流基准运行点，即I0=E(I)，代入式(11)得到基准运行点节点谐波电压为U0。

将式(12)在基准运行点处进行泰勒展开，并忽略2次及以上的高次项，可得

U=U0+ΔU=f(I0)+f′(I0)ΔI

(13)

式中：ΔI为谐波源谐波电流误差量；ΔU为节点谐波电压误差量；f′(I0)为ΔU对ΔI的灵敏度系数。

将上式中的谐波电流电压写成幅值和相角的形式，灵敏度系数用矩阵表示为：

(14)

式中：R的元素Rij=∂Ui/∂Ij;S的元素Sij=∂Ui/∂θIj;J的元素Jij=∂θUi/∂Ij;T的元素Tij=∂θUi/∂θIj。

上述参数均是基于谐波潮流基准运行点进行的计算，结果为常数，可见式(14)是一组线性方程组。

2.2　随机变量数字特征求取

采用线性化谐波方程时，谐波源谐波电流的幅值和相角与节点谐波电压的幅值和相角可记为基准量和误差量之和，为了统一分析，将这4个随机变量表示为：

X=X0+ΔX

(15)

式中：X为状态量；X0为基准量；ΔX为误差量。

由于X0是一个常数列，由式(3)和式(4)可得k阶原点矩和半不变量为：

(16)

(17)

根据半不变量的可加性，可得X与ΔX的k阶半不变量之间的关系如下：

(18)

根据半不变量的可加性和齐次性，利用式(14)的线性关系，可计算得到谐波电压幅值和相角误差量的半不变量为：

(19)

式中：gΔUk为节点谐波电压幅值误差量的半不变量；为节点谐波电压相角误差量的半不变量;gΔIk为节点谐波电流幅值误差量的半不变量；为节点谐波电流相角误差量的半不变量。

得到谐波电压半不变量后，根据式(3)转化为原点矩，利用最大熵原理进行拟合获得概率分布。

2.3　算法流程

基于半不变量和最大熵的线性化概率谐波潮流算法步骤如下。

1)读取电网基础信息，确定待求谐波潮流谐波次数，列出谐波潮流方程。

2)读取谐波源电流样本数据，确定谐波潮流基准运行点，计算该基准运行点下的谐波潮流。

3)在基准运行点处，将谐波潮流方程线性化，计算谐波电压误差量对电流误差量的灵敏度矩阵。

4)处理谐波源电流样本数据，计算谐波电流幅值误差量的半不变量和相角误差量的半不变量。

5)利用半不变量的可加性和齐次性，根据式(19)计算待求节点谐波电压幅值误差量的半不变量和相角误差量的半不变量。

6)与基准值进行叠加，根据式(18)计算谐波电压幅值半不变量和相角半不变量。

7)计算谐波电压幅值原点矩和相角原点矩。

8)建立最大熵模型，使用合适的数值分析方法求解，拟合出谐波电压的概率分布。

3　算例分析

3.1　4节点系统算例

为了说明本文方法的适用性，对一个简单4节点系统进行分析，从而便于采用卷积算法计算准确的概率分布曲线。该系统包含4个节点、3条线路、2台发电机和1台变压器，网络结构见附录A图A1。基波下的网络参数可查阅文献[21]。

节点1和3均带有非线性负载，两者互相独立。为了使卷积计算简单，只研究幅值的概率密度。假设节点1的5次谐波电流幅值服从正态分布I1～N(0.08,0.022)，节点3的5次谐波电流幅值服从正态分布I3～N(0.06,0.012)。计算网络中各节点5次谐波电压的概率密度。

首先，采用卷积法进行求解，获得准确的概率曲线。由于只研究幅值的概率密度，谐波方程中的复数运算变为模值运算，根据正态分布的卷积特性，可直接计算各节点谐波电压的概率密度函数。

利用本文方法计算节点谐波电压的概率密度，与卷积法进行对比，两种方法得到的概率密度曲线基本一致，具体可参考附录A图A2，证明了本文方法的有效性。当网络规模以及谐波源数量进一步增大时，使用直接法时卷积计算量和复杂程度大大增加，不利于实际应用，此时运用本文的方法将具有明显优势。

3.2　IEEE 57节点系统算例

采用IEEE 57节点系统作为研究对象。该系统包含57个节点、7台发电机、63条输电线路以及17台变压器，拓扑结构如图1所示。节点19，20，26，45，54，55带有非线性负载，产生并向电网中注入5次谐波。电网中发电机、变压器等设备及其他节点产生的谐波忽略不计。

图1IEEE 57节点电力系统
Fig.1IEEE 57-bus power system

以上6个谐波注入点的谐波样本取自广东某电网谐波监测点的监测数据。每个节点的谐波样本包含约16 000个5次谐波电流幅值数据，其概率密度详见附录B图B1。

运用蒙特卡洛法进行104次重复计算，结果与本文方法进行对比。由于该电力系统中节点较多，受篇幅所限，每个节点的计算结果无法一一列出，本文随机选取节点4，10，15，25，38，53进行对比分析。

节点5次谐波电压概率密度计算结果对比如图2所示。从图中可以看出，本文方法与蒙特卡洛法得到的节点谐波电压概率密度曲线基本吻合，从而证明本文方法能够有效评估谐波电压的概率密度。

基于本文方法和蒙特卡洛法评估的谐波电压幅值统计数据如表1所示，两种方法得到的谐波电压期望值和标准差相差较小，具有近似的平均值和离散度。

图2　谐波电压幅值概率密度曲线Fig.2　Probability density curves of harmonic voltage amplitude

表1谐波电压幅值统计数据
Table1Statisticaldataofharmonicvoltageamplitude

为了进一步评估本文方法的准确性，以蒙特卡洛法的结果作为对照组，计算节点谐波电压概率密度的平均误差和最大误差，用来评估概率密度曲线吻合程度，如表2所示，除了节点10以外，其他节点的概率密度评估值误差都很小，与蒙特卡洛法评估结果有较好的吻合程度。

表2　谐波电压概率密度误差分析Table 2　Analysis on probability density error of harmonic voltage

表3展示了两种方法计算而得的节点谐波电压幅值的95%概率值，本文方法得到的结果与蒙特卡洛法的结果基本相同，误差在5%以内。

表3　谐波电压幅值95%概率值Table 3　95% probability values of harmonic voltage amplitude

蒙特卡洛法因为简单准确的优点被广泛应用于随机分析方法的有效性验证中，但其计算量巨大，且随着网络规模的扩大和随机变量的增加，计算量呈指数性增长。对比蒙特卡洛法和本文方法在同一平台上的计算耗时，排除网络拓扑和参数读入、谐波导纳矩阵生成所需的时间，蒙特卡洛法耗时138.28 s，而本文方法耗时1.53 s，在计算速度方面具有明显的优势。

综上所述，本文方法得到的结果具有较高的准确度，能够合理地评估节点谐波电压的概率分布。评估得到的谐波电压期望值、标准差、95%概率值均与蒙特卡洛法的计算结果基本相同，证明本文方法能对节点谐波的分布预期、离散度做出较准确的判断，对分布规律的评估可能会出现局部误差，原因在于谐波方程线性化的线性误差、最大熵拟合的计算误差以及蒙特卡洛法本身的误差。本文方法较其他评估方法的优点主要在于能够避免大量重复模拟运算，运算时间和占用内存空间少，对于电网结构或谐波源的变化，仅需修改相应参数即可重新计算，且计算过程中无需对随机变量的分布规律做任何假设，能够避免主观判断误差，具有客观准确性。

4　结语

由于电力系统中的谐波分布具有明显的随机性，确定性谐波潮流计算无法满足实际应用的需求，本文结合数据半不变量、线性化谐波潮流方程以及最大熵规划模型，提出一种概率谐波计算方法。该方法根据谐波源电流的概率分布，计算并拟合出电网各个节点的谐波电压的概率分布曲线，进而能够为电网薄弱环节和潜在风险的评估、谐波管理等工作提供指导意见，具有重要的实际应用价值。本文方法能够避免复杂的卷积运算和大量的蒙特卡洛模拟计算，具有较高的准确性。如何减少线性化和曲线拟合等环节的误差，仍需进一步研究。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。