CN101814738A - 基于启发式能量函数的电力系统无功充裕度评估方法 - Google Patents

Abstract

本发明涉及一种基于启发式能量函数的无功充裕度评估方法。本发明针对现有技术无功充裕度评估非全局性和精度低的缺点，公开了一种基于启发式能量函数的电力系统无功充裕度评估方法，以系统全局的视角高精度的、相对快速的对电力系统无功充裕度极限评估。本发明的技术方案是，基于启发式能量函数的电力系统无功充裕度评估方法，包括以下步骤：a、引入优化乘子测量系统中各节点的临界电压值；b、建立系统启发式能量函数模型，计算系统中各节点的临界能量；c、计算系统的全局临界能量；d、根据各负荷节点的临界能量得到负荷节点的无功充裕度极限。本发明主要用于电力系统的无功充裕度评估，为系统无功补偿提供依据。

Description

基于启发式能量函数的电力系统无功充裕度评估方法

技术领域

本发明涉及电力系统稳定性技术领域，特别涉及电力系统无功功率补偿技术，确切的说，是一种基于启发式能量函数的无功充裕度评估方法。

背景技术

近年来，由于非线性负荷和感性负荷的大量使用，使得电力系统中无功功率不平衡现象越发严重，导致电压失稳甚至崩溃，影响电网的安全运行。电容器、静止无功补偿器(SVC)和新型无功发生器(ASVG)是电力系统中无功功率补偿的主要装置。要实现电力系统经济有效的无功功率补偿，必须首先对电力系统无功充裕度作出正确的评估，这是实现无功功率补偿的关键技术。对于多母线复杂的电力系统，无功功率是一个动态变化的参数，如何快速准确地评估无功充裕度一直是电力系统稳定分析与控制领域的研究热点。电力市场的发展，要求对系统无功备用进行合理定价，进一步促进了无功充裕度的评估研究。总之，系统无功充裕度的评估，无论是维护系统安全，还是适应电力市场的发展，都具有重要意义。

目前在电力系统无功充裕度评估技术方面，P-Q或Q-V平面法以及负荷裕度法是目前应用比较多的方法。P-Q或Q-V平面法主要是将有功P与无功Q之间的关系或者无功Q与电压V之间的关系在平面上描绘出来，然后根据给定的有功P或者电压U的值来确定当前系统某点的无功水平，这种方法不能从全局角度获取无功充裕度。负荷裕度法主要是从系统的整体出力来看每个负荷节点能承担的最大无功负荷量，主要采用非线性规划法求解，其模型中难于包括多种无功源，计算速度很慢，所得测算结果也不够精确。

发明内容

本发明是要解决的技术问题，就是针对现有技术无功充裕度评估非全局性和精度低的缺点，提供一种基于启发式能量函数的电力系统无功充裕度评估方法，以系统全局的视角高精度的、相对快速的对电力系统无功充裕度极限评估。

本发明解决所述技术问题，采用的技术方案是，基于启发式能量函数的电力系统无功充裕度评估方法，包括以下步骤：

a、引入优化乘子测量系统中各节点的临界电压x^cr；

b、建立系统启发式能量函数模型，计算系统中各节点的临界能量V_cri；

c、计算系统的全局临界能量V_global；

d、根据各负荷节点的临界能量得到负荷节点的无功充裕度极限。

具体的，所述临界电压值由下式得到：

其中：x^cr为临界电压；x₀为潮流对解中的初始稳态解；x_l为潮流对解中的低电压解。

进一步的，步骤b中具体包括：

b1、根据线路功率输送关系，列出系统中各节点有功功率和无功功率平衡等式如下：

$f_i(\mathrm{\δ},U)=P_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j[G_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)+B_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)]$

$g_i(\mathrm{\δ},U)={\left(U_i\right)}^{-1}[Q_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)-B_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j))];$

b2、根据能量函数表达式与功率平衡等式的关系，得到系统各节点的临界能量值，其表达式为：

$V_{\mathrm{cri}}=\∫[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}.$

$+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}$

进一步的，步骤c中根据各节点能量的迭加，得到系统全局临界能量，其表达式为：

$V_{\mathrm{global}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_{\mathrm{cri}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{G}_{\mathrm{ij}}\·\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}.$

$+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Q}_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j).U_i{|}_{U_i^s}^{U_i^u}$

特别的，采用并联电容器补偿时，系统全局临界能量表达式为：

$V_{\mathrm{cglobal}}=V_{\mathrm{global}}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_i}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2];$

步骤b2中，各节点临界能量值为：

$V_{\mathrm{crci}}=\∫[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}.$

$+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_i}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2]$

将系统各节点临界能量迭加得到的全系统临界能量，可用于判断系统全局的稳定性，b2中对应临界能量的无功充裕度值即为无功裕度极限。

本发明的有益效果是，能够在保证系统电压稳定和统一的能量框架下，以能量的形式综合给出包含多种无功源的无功充裕度，为无功补偿提供科学依据。并且，借助能量函数本身的计算优势，可以提高计算速度，有助于在正常运行状态下随时获取系统各节点的无功水平，提高电网的安全运行水平。

附图说明

图1是本发明基于启发式能量函数的无功充裕度评估流程图；

图2是IEEE-30系统结构图；

图3是无功补偿前后系统所有节点稳态电压幅值变化曲线；

图4是无功补偿前后系统所有节点临界电压幅值变化曲线；

图5是无功补偿前后系统所有节点稳态电压相角变化曲线；

图6是无功补偿前后系统所有节点临界电压相角变化曲线；

图7是无功补偿前后系统所有节点能量值的变化曲线；

图8(a)、(b)、(c)分别是无功补偿前后节点4、10、24需耗能量值随无功负荷比例k增长的变化曲线；

图9是无功补偿前后系统各负荷节点的无功裕度变化比较。

图3、图4中纵坐标表示电压的标幺值，横坐标为节点标号；图5、图6的纵坐标代表电压相角，单位：弧度，横坐标为节点标号；图7中纵坐标表示能量的标幺值，横坐标为节点标号、图8(a)、(b)、(c)中纵坐标表示能量的标幺值，横坐标为k值；图9中纵坐标表示无功裕度的标幺值，横坐标为节点标号。

具体实施方式

下面结合附图及具体实施方式，详细描述本发明的技术方案。

本发明的技术方案，将多种无功源纳入能量函数框架下，借助能量函数的快速计算优势来完成无功充裕度的评估，其实现方法如下：

步骤(1)：据改进的保留非线性潮流算法求取系统各负荷节点的临界电压值，该步骤实现过程如下：

步骤(1.1)：先给出非线性潮流一般表达式：

f_s＝f(x)                            1.1

式1.1是电力系统代数方程，其中：f_s＝[f_s1，f_s2，...，f_sn]^T，f_s1，L，f_sn为系统负荷节点的有功和无功稳态值；x＝[x₁，x₂，...x_n]^T为状态变量(这里代表电压向量)，f(x)＝[f₁(x)，f₂(x)，...，f_n(x)]^T，为系统有功功率和无功功率表达式。

步骤(1.2)：潮流方程表达成矩阵形式：

$f_s=f\left(x\right)=A\·{\left[\begin{array}{c}{x}_{1}x\\ x_{2}x\\ M\\ x_{n}x\end{array}\right]}^T---1.2$

式中A为系数矩阵。

步骤(1.3)：对(1.2)泰勒展开为：

f_s＝f(x)＝f(x₀)+J·Δx+f(Δx)                    1.3

式中：x₀为x的初始稳态解；J为潮流雅可比矩阵；Δx＝[x-x₀]＝[Δx₁，Δx₂，...，Δx_n]^T为误差修正向量。

步骤(1.4)：1.3式是一个没有截断误差的精确展开式，但是直接解这个等式有困难。本发明采用一个优化乘子μ调整Δx的大小，以加快计算速度。在1.3式中分别对第2项乘以优化乘子μ，第3项乘以μ²，得到

F＝‖f_s-f(x₀)-μJΔx-μ²f(Δx)‖²                    1.4

步骤(1.5)：1.4式对μ求导，得

$\∂F/\∂\mathrm{\μ}=g_3{\mathrm{\μ}}^3+g_2{\mathrm{\μ}}^2+g_1\mathrm{\μ}+g_0=0---1.5$

其中：g₀＝J·Δx·f(x₀)；g₁＝J^T·J·Δx^T·Δx+2f(x₀)·f(Δx)；g₂＝3J·Δx·f(Δx)；g₃＝2f^T·(Δx)·f(Δx)。(优化乘子μ，是一个数学因子，其作用是加快计算速度)

步骤(1.6)：求解1.5式中的μ，取满足1.5式F的值最小的μ为最优值μ^*。由于1.5式得到的是3次多项式，因此，μ^*有3个根μ₁、μ₂、μ₃，设μ₁＜μ₂＜μ₃，在对应一特定运行条件下μ₁、μ₃使得1.4式F为最小值，μ₂使得F为最大值。μ₁调整误差向量到接近并尽可能等于x₀，而μ₃调整误差向量到接近或等于x_l。

步骤(1.7)：根据下式计算低电压解x₁：

x_l＝x₀+μ₃·Δx

初始稳态解x₀的大小应该满足式下式：

ε_m＝max|f_si-f_i|＜ε，i＝1，2，L，2(n-1)

式中：ε_m为系统稳态有功功率值和无功功率值与当前有功值和无功值的最大失配误差；ε为允许的有功功率和无功功率失配误差；n为母线数目。

步骤(1.8)：计算临界电压。由于x₀和x_l与临界点对称，因此临界电压可由下式计算得到：

$x^{\mathrm{cr}}=\frac{x_0+x_l}{2}$

步骤(2)：构建系统的启发式能量函数模型，具体步骤如下：

步骤(2.1)：根据线路功率输送关系，列出系统中各节点有功功率和无功功率平衡等式：

$f_i(\mathrm{\δ},U)=P_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j[G_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)+B_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)]---\mathrm{2.1.1}$

$g_i(\mathrm{\δ},U)={\left(U_i\right)}^{-1}[Q_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)-B_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j))]---\mathrm{2.1.2}$

式中：P_Li是节点i的注入有功功率；Q_Li是节点i的注入无功功率；G_ij是节点i、j之间的电导；B_ij是节点i、j之间的电纳；U_i、U_j是节点i、j的电压幅值；δ_i、δ_j是节点i、j的电压相角。

步骤(2.2)：计算系统中各节点i的临界能量函数表达式。将2.1.1式、2.1.2式代入下式：

$V(X^u,X^s)={\∫}_{X^s}^{X^u}[f(\mathrm{\δ},U),g(\mathrm{\δ},U)]\·\left[\begin{array}{c}\mathrm{d\δ}\\ \mathrm{dU}\end{array}\right]$

其中：X^u代表电压幅值和相角的临界值；X^s代表电压幅值和相角的稳态值。

可得系统中节点i的临界能量V_cri的表达式为：

$V_{\mathrm{cri}}=\∫[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}$

$-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}---2.2$

如果考虑并联电容器补偿，则上式修改为：

$V_{\mathrm{crci}}=\∫[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}$

$+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_{\mathrm{ci}}}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2]---2.2c$

其中，B_ci为节点i投入的电容器电纳；

步骤(2.3)：构建系统的全局临界能量V_global表达式。由于系统的全局临界能量是各节点临界能量的迭加，因此根据2.2式，可得系统全局临界能量表达式为：

$V_{\mathrm{global}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_{\mathrm{cri}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{G}_{\mathrm{ij}}\·\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}$

$+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Q}_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j).U_i{|}_{U_i^s}^{U_i^u}---2.3$

如果考虑并联电容补偿，则此时系统的全局临界能量，记为V_cglobal，其表达：

$V_{\mathrm{cglobal}}=V_{\mathrm{global}}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_i}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2]---2.3c$

步骤(3)：推导参数变化下启发式能量函数表征的无功裕度表达式；具体步骤如下：

步骤(3.1)：考虑发电机无功约束下启发式能量函数表征的无功裕度表达式推导：

步骤(2.2)中的能量函数模型涉及的控制参数为无功负荷和发电机注入各节点的无功功率，按参数将其分解为f(U_i，δ_ij)、f(Q_Gi，U_j)、f(B_cj，U_j)、f(Q_Lj，U_j)几个部分：

$f(U_i,{\mathrm{\δ}}_{\mathrm{ij}})=U_i{U}_j^u{B}_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}}^u-U_i{U}_j^s{B}_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}}^s---\mathrm{3.1.1}$

$f(Q_{\mathrm{Gi}},U_j)=U_i{U}_j{G}_{\mathrm{ij}}\·\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}=-{\∫}_{U_j^s}^{U_j^u}\frac{Q_{\mathrm{Gi}}}{U_j}d{U}_j---\mathrm{3.1.2}$

$f(B_{\mathrm{cj}},U_j)=-B_{\mathrm{cj}}/2[{\left(U_j^u\right)}^2-{\left(U_j^s\right)}^2]---\mathrm{3.1.3}$

$f(Q_{\mathrm{Lj}},U_j)=Q_{\mathrm{Lj}}\ln U_j^u/U_j^s---\mathrm{3.1.4}$

其中：Q_Gi为发电机节点i的无功出力；G_Lj为节点j的无功负荷。

步骤(3.2)：计算系统节点j的临界能量V_crj表达式。将3.1.1、3.1.2、3.1.3、3.1.4式迭加，同时加上有功消耗的能量P_Lj(δ_ij ^u-δ_ij ^s)，有：

$V_{\mathrm{crj}}=P_{\mathrm{Lj}}({\mathrm{\δ}}_{\mathrm{ij}}^u-{\mathrm{\δ}}_{\mathrm{ij}}^s)+f(U_j,{\mathrm{\δ}}_{\mathrm{ij}})+f(Q_{\mathrm{Gi}},U_j)+f(B_{\mathrm{cj}},U_j)+f(Q_{\mathrm{Lj}},U_j)---3.2$

步骤(3.2.1)：计算单台发电机无功出力对系统势能函数的总贡献。将3.2式对Q_Gi求偏导得：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Gi}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Gi}},U_j)}\frac{\∂f(Q_{\mathrm{Gi}},U)}{\∂Q_{\mathrm{Gi}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^u}\frac{\∂X^u}{\∂Q_{\mathrm{Gi}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^s}\frac{\∂X^s}{\∂Q_{\mathrm{Gi}}}$

因在临界点处有

所以有：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Gi}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Gi}},U_j)}\frac{\∂f(Q_{\mathrm{Gi}},U_j)}{\∂Q_{\mathrm{Gi}}}=-\ln \frac{U_j^u}{U_j^s}$

单台发电机无功出力对系统势能函数的总贡献为：

$f^{\′}(Q_{\mathrm{Gi}},U_j)=-{\∫}_{Q_{\mathrm{Gi}\_0}}^{Q_{\mathrm{Gi}\_\max }}\ln \frac{U_j^u}{U_j^s}d{Q}_{\mathrm{Gi}}$

所有发电机无功出力对系统势能函数的总贡献为单台发电机对系统势能函数贡献的迭加，表达为：

$f^{\′}(Q_{\mathrm{Gi}},U_j)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\∫}_{Q_{\mathrm{Gi}\_0}}^{Q_{\mathrm{Gi}\_\max }}\ln \frac{U_j^u}{U_j^s}d{Q}_{\mathrm{Gi}}$

临界能量值表达式修正为：

$V_{\mathrm{crj}}^{\′}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_{\mathrm{Lj}}({{\mathrm{\δ}}_{\mathrm{ij}}}^u-{{\mathrm{\δ}}_{\mathrm{ij}}}^s)+f(U_j,{\mathrm{\δ}}_{\mathrm{ij}})+f^{\′}(Q_{\mathrm{Gi}},U_j)+f(B_{\mathrm{cj}},U_j)+f(Q_{\mathrm{Lj}},U_j)$

式中：Q_{Gi_0}为系统稳态时的输出；Q_{Gi_max}为发电机最大无功输出。在多机系统中，f′(Q_Gi，U_j)是所有发电机无功出力对系统势能函数贡献的总和。

步骤(3.2.2)：计算负荷的无功需求增长对系统势能函数的总贡献：

在考虑发电机无功约束条件下，当节点j无功负荷Q_Lj从初始值Q_{Lj_0}增加到使能量函数达到临界值时的Q_{Lj_max}时，有：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Lj}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Lj}},U_j)}\frac{\∂f(Q_{\mathrm{Lj}},U)}{\∂Q_{\mathrm{Lj}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^u}\frac{\∂X^u}{\∂Q_{\mathrm{Lj}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^s}\frac{\∂X^s}{\∂Q_{\mathrm{Lj}}}$

因在临界点处

所以有：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Lj}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Lj}},U_j)}\frac{\∂f(Q_{\mathrm{Lj}},U_j)}{\∂Q_{\mathrm{Lj}}}=\ln \frac{U_j^u}{U_j^s}$

对上式进行积分，得到无功负荷需要消耗的最大势能表达式为：

$\left|f^{\′}(Q_{\mathrm{Lj}},U_j)\right|=-{\∫}_{Q_{\mathrm{Lj}\_0}}^{Q_{\mathrm{Lj}\_\max }}\ln \frac{U_j^u}{U_j^s}{\mathrm{dQ}}_{\mathrm{Lj}}$

无功负荷需求需要消耗的最大能量等于系统的临界能量值V′_cr，因此可得节点j的无功负荷极限为：

$Q_{\mathrm{Lj}\_\max }=\frac{V_{\mathrm{crj}}^{\′}}{-{\ln U}_j^u/U_j^s}+Q_{\mathrm{Lj}\_0}$

本发明方法将启发式能量函数与改进的保留非线性潮流算法结合，用于构建电力系统无功充裕度评估指标，进行无功充裕性与节点电压水平之间的关系分析。

研究发现，能量函数具有易构建量化指标的优势，在对符合具有能量表征特性的参数进行量化时，能够通过能量函数的优势将多种参数的影响综合于同一框架下以能量的形式完成量化，给出具体的量化值。保留非线性潮流算法是为了改进牛顿法在处理病态条件时的缺陷提高收敛性能而提出的，为了进一步提高收敛速度，本发明在已有方法的基础上引入了优化乘子，目的是快速求解。

图1是本发明基于启发式能量函数的无功充裕度评估流程图。本发明根据改进的PNFA(保留非线性潮流算法)算法求取各负荷节点的临界电压值。计算临界电压值的速度和精度将对能量函数模型的计算速度和精度具有决定性的影响，因此计算临界电压是实现启发式能量函数分析无功裕度的关键一环。传统的PNFA算法较之牛顿法等，计算速度有很大提高，本发明在传统PNFA算法展开式的基础上，引入优化乘子μ，然后求取临界电压值。

已知非线性潮流一般表达式：

f_s＝f(x)

其中：f_s＝[f_s1，f_s2，...，f_sn]^T，f_s1，L，f_sn为给定值；x＝[x₁，x₂，...x_n]^T为状态变量，这里代表电压向量。

则潮流方程可表达成矩阵形式：

$f_s=f\left(x\right)=A\·{\left[\begin{array}{c}{x}_{1}x\\ x_{2}x\\ M\\ x_{n}x\end{array}\right]}^T$

式中A为系数矩阵。对上式泰勒展开为：

f_s＝f(x)＝f(x₀)+J·Δx+f(Δx)

式中：x₀为X的初始值(即潮流对解中的初始稳态解)；J为潮流雅可比矩阵；Δx＝[x-x₀]＝[Δx₁，Δx₂，...，Δx_n]^T为误差修正向量。

采用一个优化乘子μ调整Δx的大小，以加快计算速度。在泰勒展开式中分别对第2项乘以优化乘子μ，第3项乘以μ²，得到：

F＝‖f_s-f(x₀)-μJΔx-μ²f(Δx)‖²

上式对μ求导，得：

$\∂F/\∂\mathrm{\μ}=g_3{\mathrm{\μ}}^3+g_2{\mathrm{\μ}}^2+g_1\mathrm{\μ}+g_0=0$

其中：g₀＝J·Δx·f(x₀)；g₁＝J^T·J·Δx^T·Δx+2f(x₀)·f(Δx)；g₂＝3J·Δx·f(Δx)；g₃＝2f^T·(Δx)·f(Δx)。

求解式中的μ，使F的值最小的μ为最优值μ^*。μ^*有3个根：μ₁、μ₂、μ₃，设μ₁＜μ₂＜μ₃，在对应一特定运行条件下μ₁、μ₃使得F为最小值，μ₂使得F为最大值。μ₁调整误差向量到接近并尽可能等于x₀，而μ₃调整误差向量到接近或等于x_l。然后根据下式计算x_l：

x_l＝x₀+μ₃·Δx

x₀的大小需满足下式：

ε_m＝max|f_si-f_i|＜ε，i＝1，2，L，2(n-1)

式中：ε_m为最大失配误差；ε为潮流误差；n为母线数目。由于x₀和x_l与临界点对称，因此电压临界点为：

$x^{\mathrm{cr}}=\frac{x_0+x_l}{2}$

第二步是本发明的中心环节，目的是形成各负荷节点能量函数表达式：

$f_i(\mathrm{\δ},U)=P_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j[G_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)+B_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)]$

$g_i(\mathrm{\δ},U)={\left(U_i\right)}^{-1}[Q_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)-B_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j))]$

式中：P_Li是节点i的注入有功功率；Q_Li是节点i的注入无功；G_ij是节点i、j之间的电导；B_ij是节点i、j之间的电纳；δ_i、δ_j是节点i、j的电压相角。

将以上两式代入启发式能量函数表达式：

$V(X^u,X^s)={\∫}_{X^s}^{X^u}[f(\mathrm{\δ},U),g(\mathrm{\δ},U)]\·\left[\begin{array}{c}\mathrm{d\δ}\\ \mathrm{dU}\end{array}\right]$

可得系统中第i节点的临界能量V_cri的表达式为：

$V_{\mathrm{cri}}=\∫[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}$

$+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}$

如果考虑并联电容器补偿，则上式修改为：

$V_{\mathrm{cri}}=\∫[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}d{\mathrm{\δ}}_i\\ d{U}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}$

$+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_{\mathrm{ci}}}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2]$

第三步利用第二步得到的系统中节点临界能量函数表达式形成全系统总临界能量函数表达式。系统的全局临界能量函数是系统所有节点临界能量函数的迭加，因此可得系统全局的临界能量为：

$V_{\mathrm{global}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_{\mathrm{cri}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{G}_{\mathrm{ij}}\·\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}$

$+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Q}_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j).U_i{|}_{U_i^s}^{U_i^u}$

如果考虑并联电容补偿，则系统的全局临界能量函数表达为：

$V_{\mathrm{cglobal}}=V_{\mathrm{global}}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_{\mathrm{ci}}}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2]$

第四步是形成参数变化下启发式能量函数表征的无功裕度表达式

形成该式总共需要两步：

①是求取无功补偿前、后临界能量值表达式，首先将能量模型进行分解，在上述的能量函数模型中，涉及的控制参数为无功负荷和发电机注入各节点的无功，按参数将能量模型分解为f(U_i，δ_ij)、f(Q_Gi，U_j)、f(B_cj，U_j)、f(Q_Lj，U_j)几个部分：

$f(U_i,{\mathrm{\δ}}_{\mathrm{ij}})=U_i{U}_j^u{B}_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}}^u-U_i{U}_j^s{B}_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}}^s$

$f(Q_{\mathrm{Gi}},U_j)=U_i{U}_j{G}_{\mathrm{ij}}\·\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}=-{\∫}_{U_j^s}^{U_j^u}\frac{Q_{\mathrm{Gi}}}{U_j}d{U}_j$

$f(B_{\mathrm{cj}},U_j)=-B_{\mathrm{cj}}/2[{\left(U_j^u\right)}^2-{\left(U_j^s\right)}^2]$

$f(Q_{\mathrm{Lj}},U_j)=Q_{\mathrm{Lj}}\ln U_j^u/U_j^s$

将上式进行叠加，同时加上有功消耗的能量，就得j节点的临界能量表达式如下：

$V_{\mathrm{crj}}=P_{\mathrm{Lj}}({\mathrm{\δ}}_{\mathrm{ij}}^u-{\mathrm{\δ}}_{\mathrm{ij}}^s)+f(U_j,{\mathrm{\δ}}_{\mathrm{ij}})+f(Q_{\mathrm{Gi}},U_j)+f(B_{\mathrm{cj}},U_j)+f(Q_{\mathrm{Lj}},U_j)$

接着计算单台发电机无功出力对系统势能函数的总贡献。上式对Q_Gi求偏导得：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Gi}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Gi}},U_j)}\frac{\∂f(Q_{\mathrm{Gi}},U)}{\∂Q_{\mathrm{Gi}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^u}\frac{\∂X^u}{\∂Q_{\mathrm{Gi}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^s}\frac{\∂X^s}{\∂Q_{\mathrm{Gi}}}$

因在临界点处有

所以有：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Gi}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Gi}},U_j)}\frac{\∂f(Q_{\mathrm{Gi}},U_j)}{\∂Q_{\mathrm{Gi}}}=-\ln \frac{U_j^u}{U_j^s}$

所以单台发电机无功出力对系统势能函数的总贡献为：

$f^{\′}(Q_{\mathrm{Gi}},U_j)=-{\∫}_{Q_{\mathrm{Gi}\_0}}^{Q_{\mathrm{Gi}\_\max }}\ln \frac{U_j^u}{U_j^s}d{Q}_{\mathrm{Gi}}$

节点j的临界能量值表达式修正为：

V′_crj＝P_Lj(δ_ij ^u-δ_ij ^s)+f(U_j，δ_ij)+f′(Q_Gi，U_j)+f(B_cj，U_j)+f(Q_Lj，U_j)

式中：Q_{Gi_0}为系统稳态时的输出；Q_{Gi_max}为发电机最大无功输出。在多机系统中，f′(Q_Gi，U_j)是所有发电机无功出力对系统势能函数贡献的总和。式中f(B_cj，U_j)为考虑并联电容器补偿后增加的函数项，补偿前的表达式中没有此项。

②是计算负荷的无功需求增长对系统势能函数的总贡献：

在考虑发电机无功约束条件下，当节点j的无功负荷Q_Lj以增长因子k(无功负荷增长量与初始无功负荷大小的比值)从初始值Q_{Lj_0}增加到使能量函数达到临界值时的Q_{Lj_max}时，有：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Lj}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Lj}},U_j)}\frac{\∂f(Q_{\mathrm{Lj}},U)}{\∂Q_{\mathrm{Lj}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^u}\frac{\∂X^u}{\∂Q_{\mathrm{Lj}}}+\frac{\∂V_{\mathrm{crj}}}{\∂X^s}\frac{\∂X^s}{\∂Q_{\mathrm{Lj}}}$

因在临界点处

所以有：

$\frac{{\mathrm{dV}}_{\mathrm{crj}}}{{\mathrm{dQ}}_{\mathrm{Lj}}}=\frac{\∂V_{\mathrm{crj}}}{\∂f(Q_{\mathrm{Lj}},U_j)}\frac{\∂f(Q_{\mathrm{Lj}},U_j)}{\∂Q_{\mathrm{Lj}}}=\ln \frac{U_j^u}{U_j^s}$

对上式进行积分，得到无功负荷需要消耗的最大能量表达式为：

$\left|f^{\′}(Q_{\mathrm{Lj}},U_j)\right|=-{\∫}_{Q_{\mathrm{Lj}\_0}}^{Q_{\mathrm{Lj}\_\max }}\ln \frac{U_j^u}{U_j^s}{\mathrm{dQ}}_{\mathrm{Lj}}$

无功负荷需要消耗的最大能量等于系统的临界能量值V′_cr，因此可得节点j的无功负荷极限为：

$Q_{\mathrm{Lj}\_\max }=\frac{V_{\mathrm{crj}}^{\′}}{-{\ln U}_j^u/U_j^s}+Q_{\mathrm{Lj}\_0}$

上式即是参数变化下启发式能量函数表征的无功裕度表达式，分别代入相应的V′_crj，即可得到补偿前和补偿后的无功裕度表达式，在系统中，上式为全部的负荷节点无功增加到最大所消耗的最大能量总和。Q_{Lj_max}就是需要求取的最大无功裕度，也是系统各负荷节点能承受的最大无功负荷极限。

图2是IEEE-30母线系统结构图，IEEE-30母线系统共包括6台发电机和30个节点，6台发电机分别安装在标号为1、2、5、8、11、13节点，其中在标号为10、24的节点并联电容器进行无功补偿。标号为3、4、7、10、12、14、15、16、17、18、19、20、21、23、24、26、29、30的节点为负荷节点。

对该IEEE-30母线系统的能量函数模型进行仿真，结果如图3～9所示。

图3、图4分别给出了并联电容器补偿前后系统所有节点稳态电压和临界电压的变化曲线。当在节点10、24处并联电容进行无功补偿后，各节点的稳态电压和临界电压都比补偿前有所提升，这与在2母线系统中得出的结论一致。还可观察到，无论是稳态电压还是临界电压，距离补偿点越远的地方变化越不明显，而距离补偿地点相对越近的地方，电压提升相对多一些。

图5、6是无功补偿前后系统所有节点电压的稳态和临界电压相角的变化情况。稳态和临界电压相角在补偿前后变化较多的节点均集中在标号为10～30的节点之间，说明无功补偿将同时影响补偿地点附近的节点电压和相角。并且，节点10、24在补偿后的稳态电压相角绝对值都比补偿前有所增大，但节点10、24的临界电压相角在补偿前后并没表现出与稳态电压相角一致的变化，节点10在补偿后临界相角绝对值稍微增大，而节点24的相角在补偿后增加不明显，说明无功补偿装置的投入，对临界相角的增加具有抑制作用。

图7是无功补偿前后系统所有节点能量值的变化。将得到的补偿前后稳态电压幅值、稳态电压相角和临界电压幅值、临界电压相角代入启发式能量函数模型，得到能量函数模型的数值仿真结果。从图7虚实曲线的对比发现，并不是所有节点在无功补偿后临界能量值都升高，结合补偿前后各节点临界电压的变化曲线发现，临界电压升高较多的节点的能量值有下降趋势，这是应该注意的问题。

从图8可看出，某节点的无功负荷越大，需要消耗的能量值越大，系统距离临界能量点的距离越小，节点10、24处的能量函数在补偿后增长趋势较补偿前平缓很多，因此到达临界能量值需要的时间变长，说明无功补偿之后，系统稳定性加强。负荷节点4是远离补偿的所有节点之一，从图8a可见，补偿前后节点4的需耗能量值随无功负荷比例k的增长其变化趋势几乎一样，说明补偿与否对该节点的稳定性基本没太大影响，即距离补偿点很远的节点较少受益于无功补偿；从图8b可以看出，补偿后节点10的无功需耗能量值随无功负荷比例k的增长变化趋势比补偿前变化趋缓，说明补偿后该节点的稳定性增强；从图8c可以看出，补偿后节点24的无功需耗能量值随无功负荷比例k的增长变化趋势比补偿前变化趋缓很多，说明补偿后该节点的稳定性增强很多，即无功补偿只对本补偿节点稳定性影响最大。

图9是无功补偿前后各负荷节点的无功裕度变化，从图中虚、实线的对比发现，无功补偿后，除节点3以外，所有节点的无功裕度都得到不同程度的提高。并且，无功补偿后的系统总能量值比补偿前提高0.8左右，说明局部的无功补偿对整个系统稳定性有益。

尽管上面对本发明说明性的具体实施方式进行了描述，但应当清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (6)

1.基于启发式能量函数的电力系统无功充裕度评估方法，包括以下步骤：

a、引入优化乘子测量系统中各节点的临界电压x^cr；

b、建立系统启发式能量函数模型，计算系统中各节点的临界能量V_cri；

c、计算系统全局临界能量V_global；

d、根据各负荷节点的临界能量得到负荷节点的无功充裕度极限。

2.根据权利要求1所述的基于启发式能量函数的电力系统无功充裕度评估方法，其特征在于，所述临界电压值由下式得到：

其中：x^cr为临界电压；x₀为潮流对解中的初始稳态解；x_l为潮流对解中的低电压解。

3.根据权利要求1或2所述的基于启发式能量函数的电力系统无功充裕度评估方法，其特征在于，步骤b中具体包括：

b1、根据线路功率输送关系，列出系统中各节点有功功率和无功功率平衡等式如下：

$f_i(\mathrm{\δ},U)=P_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j[G_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)+B_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)]$

$g_i(\mathrm{\δ},U)={\left(U_i\right)}^{-1}[Q_{\mathrm{Li}}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)-B_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j))];$

b2、根据能量函数表达式与功率平衡等式的关系，得到系统各节点临界能量，其表达式为：

$V_{\mathrm{cri}}=\∫[f_i,g_i]\·\left[\begin{array}{c}{\mathrm{d\δ}}_i\\ {\mathrm{dU}}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}{\mathrm{d\δ}}_i\\ {\mathrm{dU}}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}.$

$+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){\·U}_i{|}_{U_i^s}^{U_i^u}$

4.根据权利要求3所述的基于启发式能量函数的电力系统无功充裕度评估方法，其特征在于，步骤c中根据各节点能量的迭加，得到系统全局临界能量V_global，其表达式为：

$V_{\mathrm{global}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_{\mathrm{cri}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{G}_{\mathrm{ij}}\·\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_i{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}.$

$+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Q}_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j).U_i{|}_{U_i^s}^{U_i^u}$

5.根据权利要求4所述的基于启发式能量函数的电力系统无功充裕度评估方法，其特征在于，采用并联电容器补偿时，系统全局临界能量，记为V_cglobal，其表达式为：

$V_{\mathrm{cglobal}}=V_{\mathrm{global}}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_i}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2].$

6.根据权利要求3所述的基于启发式能量函数的电力系统无功充裕度评估方法，其特征在于，采用并联电容器补偿时，步骤b2中，各节点临界能量值表达式为：

$V_{\mathrm{crci}}=\∫[f_i,g_i]\·\left[\begin{array}{c}{\mathrm{d\δ}}_i\\ {\mathrm{dU}}_i\end{array}\right]={\∫}_{({\mathrm{\δ}}_i^s,U_i^s)}^{({\mathrm{\δ}}_i^u,U_i^u)}[f_i,g_i]\·\left[\begin{array}{c}{\mathrm{d\δ}}_i\\ {\mathrm{dU}}_i\end{array}\right]$

$=P_{\mathrm{Li}}({\mathrm{\δ}}_i^u-{\mathrm{\δ}}_i^s)-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}+Q_{\mathrm{Li}}\ln \frac{U_i^u}{U_i^s}-U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j){|}_{{\mathrm{\δ}}_i^s}^{{\mathrm{\δ}}_i^u}.$

$+\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{B}_{\mathrm{ij}}\cos ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)U_i{|}_{U_i^s}^{U_i^u}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j{G}_{\mathrm{ij}}\sin ({\mathrm{\δ}}_i-{\mathrm{\δ}}_j)\·U_i{|}_{U_i^s}^{U_i^u}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{B_i}{2}[{\left(U_i^u\right)}^2-{\left(U_i^s\right)}^2]$

Images