CN104052063A - 一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法 - Google Patents

Abstract

本发明公开了一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，舍弃了根据一元二次方程实根解的存在条件以及扩展雅克比矩阵的行列式建立电压稳定性局部指标的方法，通过求解系统在鞍结分岔点处的电压和临界负荷水平因子，建立计及ZIP静态负荷影响的电压稳定性局部指标，它准确地反映了负荷电压特性对系统电压稳定性的影响；对于求解临界负荷水平因子时遇到的非线性方程组，考虑到局部牛顿法的收敛性对初值的选取敏感以及收敛速度等问题，采用全局收敛且自适应最优迭代步长的NLEQ-RES法求解非线性方程组，改善了求解非线性方程组时数值计算的收敛性及其收敛速度。

Description

一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法

技术领域

本发明涉及一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法。

背景技术

电力系统电压失稳可能导致大规模停电事故，造成巨大经济损失，严重影响社会生活，利用节点电压指标实现对电网电压稳定性的快速监测对确保电力系统的安全稳定运行意义重大。现有的节点电压指标法一般都基于简单传输系统，假设负荷为恒功率负荷模型，以负荷侧母线电压为变量建立方程，直接或间接利用方程实根解的存在条件建立指标，其稳定临界点是系统最大功率传输点。但已有众多研究指出，当负荷为非恒功率模型时，系统电压失稳点不是系统最大功率传输点而是鞍结分岔点，实际上由于电力电子开关元件的大规模使用，恒功率负荷模型已不能真实的反映系统负荷特性。

为了在节点指标法中计及负荷电压特性的影响，有文献提出了改进的L指标L₁指标，以及在阻抗模比值指标的基础上，将ZIP负荷的恒阻抗负荷和恒电流负荷直接折算到构造的单电源功率传输等值系统的等值阻抗和等值电源中，计算VSLB指标，或者基于L_p指标和L_q指标的思想提出了一种计及负荷电压特性影响的电压稳定性指标。但它们都采用或者基于方程解的存在条件建立稳定性指标的理论思路，指标值不能准确地计及负荷电压特性对系统电压稳定性的影响。另外，有文献利用计入了静态负荷模型的节点扩展雅克比矩阵行列式是否等于零作为电压稳定的判据，但逼近分岔点时，由于除最小特征值外，雅克比矩阵的其它特征值也可能较大，使得扩展雅克比矩阵行列式远大于零，因此它在实际应用时有误差大的缺陷。

发明内容

本发明的目的是克服现有技术的不足，舍弃仅根据方程实根解的存在条件和雅克比矩阵行列式建立指标的方法，通过引入简单传输系统达到鞍结分岔点时的特征方程，给出了一种计及静态负荷特性的功率偏差指标，准确地反映了负荷静态特性对系统电压稳定性的影响；对于求解计及静态负荷特性的节点电压稳定指标遇到的非线性方程组问题，本发明将NLEQ-RES算法应用于求解提出的电压稳定指标，它克服了用常规牛顿法求解上述非线性方程组初值选取困难以及收敛速度慢的问题。

本发明实现其发明目的所采用的技术方案是：一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，包括以下步骤：

A、根据含负荷节点简单传输系统达到鞍结分岔点条件建立非线性函数F₁(X)、F₂(X)、F₃(X)

基于已建立的含负荷节点简单传输系统的模型及参数，根据该系统达到鞍结分岔点的特征方程建立三个非线性函数F₁(X)、F₂(X)、F₃(X)，各函数的变量X分别为向量[U_r,k_p,max]的转置[U_r,k_p,max]^T、[U_r,k_q,max]的转置[U_r,k_q,max]^T和[U_r,k_s,max]的转置[U_r,k_s,max]^T，U_r为被监测节点电压幅值，k_p,max、k_q,max和k_s,max分别为被监测节点有功功率、无功功率、视在功率负荷因子最大值；

B、利用NLEQ-RES算法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0得到k_p,max、k_q,max和k_s,max

根据A中建立的非线性函数F₁(X)、F₂(X)、F₃(X)，采用具有全局收敛且能自适应最优迭代步长的基于收敛准则和自适应置信区间策略的带残差全局牛顿法即NLEQ-RES法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0，得到被监测节点有功功率、无功功率、视在功率负荷因子最大值k_p,max、k_q,max和k_s,max，每次迭代求解时均自适应调整迭代步长；

C、计算被监测节点计及静态负荷特性的电压稳定指标

根据B中计算得到的负荷因子最大值k_p,max、k_q,max和k_s,max，得到被监测节点的负荷裕度k_p,margin、k_q,margin和k_s,margin，计算反映有功功率、无功功率以及视在功率裕度的指标I_VSIP、I_VSIQ和I_VSIS，基于此构造并计算被监测节点计及静态负荷特性的电压稳定指标I_VSI；

D、判断

根据C步计算出的被监测节点的电压稳定指标，判断被监测节点是否为薄弱节点；如果是薄弱节点，则进行电压稳定预警，否则不作预警；判断结束后，返回步骤A继续执行下一时刻的电压稳定监测。由于电力系统会随着时间的变化而变化，在本专利中的表现即是含负荷节点简单传输系统参数会随时间而变化，步骤A-D给出了当前时刻下电压稳定的监测方法，当执行完步骤D的判断之后，当前时刻的电压稳定状态已经明了，当前时刻下的算法结束。然而电力系统会随时间而变化，因此，当下一时刻时，算法会重新执行步骤A，继续监测下一时刻的系统电压稳定状态。

所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，所述A步根据含负荷节点简单传输系统达到鞍结分岔点条件建立非线性函数F₁(X)、F₂(X)、F₃(X)的具体步骤为：

已建立的含负荷的被监测节点r处的简单传输系统参数如下：U_s∠δ为电源端电压，U_r∠0为负荷侧电压，Z＝R+jX为线路阻抗，Y_L＝Y_p+jY_q为ZIP静态负荷中恒阻抗负荷分量的导纳值，I_L＝I_p+jI_q为ZIP静态负荷中恒电流负荷分量的电流值，S_con＝P_con+jQ_con为ZIP静态负荷中恒功率负荷分量的功率值；

被监测节点负荷因子最大值为有功功率负荷因子最大值k_p,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₁(X)为：

$F_1\left(X\right)=\left(\begin{array}{c}{f}_1\left(X\right)\\ g_1\left(X\right)\end{array}\right)---\left(1\right)$

其中：方程组的变量X＝[U_r,k_p,max]^T，U_r为被监测节点电压幅值，k_p,max为被监测节点有功功率负荷因子最大值；

$f_1\left(X\right)=U_r^4+b_1{U}_r^2+c_1---\left(2\right)$

$\begin{array}{c}{g}_1\left(X\right)=4U_r^3+{2b}_1{U}_r^2+4({\mathrm{RY}}_p{k}_{p,\max }+{\mathrm{XY}}_q)U_r^3+2(R^2+X^2)\·\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)k_{p,\max }^2(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)(I_q+{2Y}_q{U}_r)]\end{array}---\left(3\right)$

b₁、c₁的表达式为式(4)、(5)，

b₁＝2R(P_con+I_pU_r+Y_pU_r ²)k_p,max+2X(Q_con+I_qU_r+Y_qU_r ²)-U_s ²    (4)

$c_1=(R^2+X^2)[{(P_{\mathrm{con}}+I_p{U}_r+Y_p{U_r}^2)}^2{k}_{p,\max }^2+{(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U_r}^2)}^2]---\left(5\right)$

被监测节点负荷因子最大值为无功功率负荷因子最大值k_q,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₂(X)为：

$F_2\left(X\right)=\left(\begin{array}{c}{f}_2\left(X\right)\\ g_2\left(X\right)\end{array}\right)---\left(6\right)$

其中：方程组的变量X＝[U_r,k_q,max]^T，U_r为被监测节点电压幅值，k_q,max为被监测节点有功功率负荷因子最大值；

$f_2\left(X\right)=U_r^4+b_2{U}_r^2+c_2---\left(7\right)$

$\begin{array}{c}{g}_2\left(X\right)={4U}_r^3+{2b}_2{U}_r^2+4({\mathrm{RY}}_p+{\mathrm{XY}}_q{k}_{q,\max })U_r^3+2(R^2+X^2)\·\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)k_{q,\max }^2(I_q+{2Y}_q{U}_r)]\end{array}---\left(8\right)$

b₂、c₂的表达式为式(9)、(10)，

b₂＝2R(P_con+I_pU_r+Y_pU_r ²)+2X(Q_con+I_qU_r+Y_qU_r ²)k_q,max-U_s ²    (9)

c₂＝(R²+X²)[(P_con+I_pU_r+Y_pU_r ²)²+(Q_con+I_qU_r+Y_qU_r ²)²k_q,max ²]    (10)

被监测节点负荷因子最大值为视在功率负荷因子最大值k_s,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₃(X)为：

$F_3\left(X\right)=\left(\begin{array}{c}{f}_3\left(X\right)\\ g_3\left(X\right)\end{array}\right)---\left(11\right)$

其中：方程组的变量X＝[U_r,k_s,max]^T，U_r为被监测节点电压幅值，k_s,max为被监测节点有功功率负荷因子最大值；

$f_3\left(X\right)=U_r^4+b_3{U}_r^2+c_3---\left(12\right)$

$\begin{array}{c}{g}_3\left(X\right)={4U}_r^3+{2b}_3{U}_r^2+{4k}_{s,\max }({\mathrm{RY}}_p+{\mathrm{XY}}_q)U_r^3+{2k}_{s,\max }^2(R^2+X^2)\·\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)(I_q+{2Y}_q{U}_r)]\end{array}---\left(13\right)$

b₃、c₃的表达式为式(14)、(15)，

b₂＝2k_s,max[R(P_con+I_pU_r+Y_pU_r ²)+X(Q_con+I_qU_r+Y_qU_r ²)]-U_s ²    (14)

$c_3=(R^2+X^2)k_{s,\max }^2[{(P_{\mathrm{con}}+I_p{U}_r+Y_p{U_r}^2)}^2+{(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U_r}^2)}^2]---\left(15\right)$

所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，所述B步利用NLEQ-RES算法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0得到k_p,max、k_q,max和k_s,max的步骤为：

B1、得到用于求解被监测节点有功功率负荷因子最大值k_p,max的非线性方程组F(X)＝0

令F(X)＝F₁(X)、f(X)＝f₁(X)、g(X)＝g₁(X)，构造非线性方程组：

$F\left(X\right)=\left(\begin{array}{c}f\left(X\right)\\ g\left(X\right)\end{array}\right)=0---\left(16\right)$

其中：方程组的变量X＝[U_r,k_p,max]^T，U_r为被监测节点电压幅值，k_p,max为被监测节点有功功率负荷因子最大值；

B2、判断第k(k＝0,1,…)次牛顿迭代是否满足收敛条件

如果||F(X^(k))||≤ε，执行B10；其中：ε为任意满足0<ε≤10^-5的常数；

否则，执行B3；

B3、计算第k(k＝0,1,…)次牛顿迭代时F(X^(k))＝0的雅克比矩阵F'(X^(k))

$F^{\&prime;}\left(X^{\left(k\right)}\right)=\left[\begin{array}{cc}{\frac{\&PartialD;f}{{\&PartialD;U}_r}|}_k& {\frac{\&PartialD;f}{{\&PartialD;k}_{p,\max }}|}_k\\ {\frac{\&PartialD;g}{{\&PartialD;U}_r}|}_k& {\frac{\&PartialD;g}{{\&PartialD;k}_{p,\max }}|}_k\end{array}\right]---\left(17\right)$

$\begin{array}{c}{\frac{\&PartialD;f}{{\&PartialD;U}_r}|}_k=4{\left(U_r^{\left(k\right)}\right)}^3+2b{\left(U_r^{\left(k\right)}\right)}^2+4({\mathrm{RY}}_p{k}_{p,\max }^{\left(k\right)}+{\mathrm{XY}}_q){\left(U_r^{\left(k\right)}\right)}^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(k_{p,\max }^{\left(k\right)}\right)}^2(I_p+{2Y}_p{U}_r^{\left(k\right)})\\ +(Q_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_q{\left(U_r^{\left(k\right)}\right)}^2)(I_q+{2Y}_q{U}_r^{\left(k\right)})]\end{array}---\left(18\right)$

$\begin{array}{c}{\frac{\&PartialD;f}{{\&PartialD;k}_{p,\max }}|}_k=2R(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(U_r^{\left(k\right)}\right)}^2\\ +2(R^2+X^2){(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2)}^2{k}_{p,\max }^{\left(k\right)}\end{array}---\left(19\right)$

$\begin{array}{c}{\frac{\&PartialD;g}{{\&PartialD;U}_r}|}_k=12{\left(U_r^{\left(k\right)}\right)}^2+4b{U}_r^{\left(k\right)}+2[2R(I_p+{2Y}_p{U_r}^{\left(k\right)})k_{p,\max }^{\left(k\right)}+2X(I_q+{2Y}_q{U_r}^{\left(k\right)})]{\left(U_r^{\left(k\right)}\right)}^2\\ +12({\mathrm{RY}}_p{k}_{p,\max }^{\left(k\right)}+{\mathrm{XY}}_q){\left(U_r^{\left(k\right)}\right)}^2+2(R^2+X^2)[{(I_p+{2Y}_p{U}_r^{\left(k\right)})}^2{\left(k_{p,\max }^{\left(k\right)}\right)}^2\\ +{2Y}_p(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(k_{p,\max }^{\left(k\right)}\right)}^2+{(I_q+{2Y}_q{U}_r^{\left(k\right)})}^2\\ +{2Y}_q(Q_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_q{\left(U_r^{\left(k\right)}\right)}^2)]\end{array}---\left(20\right)$

$\begin{array}{c}{\frac{\&PartialD;g}{{\&PartialD;k}_{p,\max }}|}_k=4R(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(U_r^{\left(k\right)}\right)}^2+4R{Y}_p{\left(U_r^{\left(k\right)}\right)}^3\\ +4(R^2+X^2){(P_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2)}^2{k}_{p,\max }^{\left(k\right)}\end{array}---\left(21\right)$

B4、计算第k(k＝0,1,…)次牛顿迭代时的修正量ΔX^(k)

ΔX^(k)＝-[F'(X^(k))]^-1F(X^(k))    (22)

B5、计算第k(k＝0,1,…)次牛顿迭代时的阻尼因子估计值[λ_k]

$\left[{\mathrm{\&lambda;}}_k\right]=\left\{\begin{array}{cc}1& k=0\\ \min (\mathrm{1,1}/[h_k\left]\right)& k>0\end{array}\right.---\left(23\right)$

$\begin{array}{cc}\left[h_k\right]=\frac{\left|\right|F\left(X^{\left(k\right)}\right)\left|\right|}{\left|\right|F\left(X^{(k-1)}\right)\left|\right|}\&CenterDot;\left[h_{k-1}^{\&prime;}\right]& k>0---\left(24\right)\end{array}$

B6、根据阻尼因子判断方程组是否奇异

如果[λ_k]≤λ_min，则表明沿着牛顿路径牛顿解逼近于雅克比矩阵奇异点，结束计算，执行B10；其中：λ_min为任意满足0<λ_min≤10^-5的常数；

否则，执行B7；

B7、计算第k(k＝0,1,…)次牛顿迭代时的方程的解X^(k+1)和F(X^(k))

X^(k+1)＝X^(k)+[λ_k]ΔX^(k)    (25)

将X^(k+1)代入函数，计算得到F(X^(k))；

B8、计算监测量[h_k′]

$\left[h_k^{\&prime;}\right]=\frac{2\left|\right|F\left(X^{(k+1)}\right)-(1-[{\mathrm{\&lambda;}}_k\left]\right)F\left(X^{\left(k\right)}\right)\left|\right|}{{\left(\right[{\mathrm{\&lambda;}}_k\left]\right)}^2\left|\right|F\left(X^{\left(k\right)}\right)\left|\right|}---\left(26\right)$

B9、判断第k次牛顿迭代的迭代值X^(k)

如果||F(X^k+1)||≥(1-0.25·[λ_k])||F(X^k)||，

[λ_k′]＝min(0.5·[λ_k],1/[h_k′])    (27)

令[λ_k]＝[λ_k′]，返回B6；

否则，返回B2；

B10、得到被监测节点有功功率负荷因子最大值k_p,max

由B2～B9得到非线性方程组F(X^(k))＝0的解为：X＝[U_r,k_p,max]^T＝X^(k)；

由此得到k_p,max；

B11、利用NLEQ-RES算法求解被监测节点负荷因子最大值k_q,max和k_s,max

令B1中F(X)＝F₂(X)、f(X)＝f₂(X)、g(X)＝g₂(X)，构造非线性方程组F(X)＝0，得到关于k_q,max的方程组式(16)，然后利用方法B2～B10计算得到k_q,max；

令B1中F(X)＝F₃(X)、f(X)＝f₃(X)、g(X)＝g₃(X)，构造非线性方程组F(X)＝0，得到关于k_s,max的方程组式(16)，然后利用步骤B2～B10计算得到k_s,max。

所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，所述C步计算被监测节点计及静态负荷特性的电压稳定指标的步骤为：

C1、计算被监测节点的负荷裕度k_p,margin、k_q,margin和k_s,margin

$\left\{\begin{array}{c}{k}_{p,m\arg \mathrm{in}}=k_{p,\max }-1\\ k_{q,m\arg \mathrm{in}}=k_{q,\max }-1\\ k_{s,m\arg \mathrm{in}}=k_{s,\max }-1\end{array}\right.---\left(28\right)$

C2、计算反映有功功率、无功功率以及视在功率裕度的指标I_VSIP、I_VSIQ和I_VSIS

$\left\{\begin{array}{c}{I}_{\mathrm{VSIP}}=\left|\frac{k_{p,m\arg \mathrm{in}}}{k_{p,\max }}\right|\\ I_{\mathrm{VSIQ}}=\left|\frac{k_{q,m\arg \mathrm{in}}}{k_{q,\max }}\right|\\ I_{\mathrm{VSIS}}=\left|\frac{k_{s,m\arg \mathrm{in}}}{k_{s,\max }}\right|\end{array}\right.---\left(29\right)$

由I_VSIP、I_VSIQ和I_VSIS可以直接得到各负荷裕度

C3、构造电压稳定指标I_VSI

$I_{\mathrm{VSI}}=\min \left(\right|\frac{k_{p,m\arg \mathrm{in}}}{k_{p,\max }}|,|\frac{k_{q,m\arg \mathrm{in}}}{k_{q,\max }}|,|\frac{k_{s,m\arg \mathrm{in}}}{k_{s,\max }}\left|\right)---\left(30\right)$

所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，所述D步判断的具体方法为：

当I_VSI>0时，系统工作在分岔曲线上半支，系统电压稳定；I_VSI<0时，系统失去电压稳定，进行电压稳定预警；I_VSI＝0时，系统处于临界情况，进行电压稳定预警。

本发明的技术效果在于，舍弃了根据一元二次方程实根解的存在条件以及扩展雅克比矩阵的行列式建立电压稳定性局部指标的方法，通过求解系统在鞍结分岔点处的电压和临界负荷水平因子，建立计及ZIP静态负荷影响的电压稳定性局部指标，它准确地反映了负荷电压特性对系统电压稳定性的影响；对于求解临界负荷水平因子时遇到的非线性方程组，考虑到局部牛顿法的收敛性对初值的选取敏感以及收敛速度等问题，采用全局收敛且自适应最优迭代步长的NLEQ-RES法求解非线性方程组，改善了求解非线性方程组时数值计算的收敛性及其收敛速度。

下面结合附图和具体的实施方式，对本发明作进一步详细的说明。

附图说明

图1为本发明所用到的已建立的含负荷的被监测节点r处的简单传输系统。

图2为在IEEE50机测试系统中母线51含不同负荷成份时母线51的电压稳定指标轨迹及计及了负荷静态特性影响的V-Q雅可比矩阵最小特征值轨迹。

图3为在改进IEEE50机测试系统中母线120含不同负荷成份时母线120的电压稳定指标轨迹及计及了负荷静态特性影响的V-Q雅可比矩阵最小特征值轨迹。

图2(a)、图2(b)、图3(a)和图3(b)中，横轴为功率，单位为标幺值(pu)，纵轴为指标值。由“——”构成的曲线，为采用本发明实施例方法得到的电压稳定指标轨迹；由“-----”构成的曲线，为1减L₁指标轨迹；由构成的曲线，为最小特征值轨迹；由构成的曲线，为1减VSLB指标轨迹。

具体实施方式

本方法实施例包括以下步骤：

A、根据已建立的含负荷节点简单传输系统达到鞍结分岔点条件建立非线性函数F₁(X)、F₂(X)、F₃(X)

图1示出已建立的含负荷的被监测节点r处的简单传输系统，该系统可以利用专利ZL200910164280.9的方法基于广域量测系统的量测数据辨识得到。其中，U_s∠δ为电源端电压，U_r∠0为负荷侧电压，Z_L＝R+jX为线路阻抗，负荷为ZIP静态负荷，其中阻抗负荷分量的导纳值为Y_L＝Y_p+jY_q，电流负荷分量的电流值为I_L＝I_p+jI_q，恒功率负荷分量为S_con＝P_con+jQ_con。本例中对简单传输系统参数采用带遗忘因子的平方根最小二乘估计进行辨识，辨识时遗忘因子取值等于0.15，跳数n＝2。

被监测节点负荷因子最大值为有功功率负荷因子最大值k_p,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₁(X)为：

$F_1\left(X\right)=\left(\begin{array}{c}{f}_1\left(X\right)\\ g_1\left(X\right)\end{array}\right)---\left(1\right)$

其中：方程组的变量X＝[U_r,k_p,max]^T，U_r为被监测节点电压幅值，k_p,max为被监测节点有功功率负荷因子最大值；

$f_1\left(X\right)=U_r^4+b_1{U}_r^2+c_1---\left(2\right)$

$\begin{array}{c}{g}_1\left(X\right)=4U_r^3+{2b}_1{U}_r^2+4({\mathrm{RY}}_p{k}_{p,\max }+{\mathrm{XY}}_q)U_r^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)k_{p,\max }^2(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)(I_q+{2Y}_q{U}_r)]\end{array}---\left(3\right)$

b₁、c₁的表达式为式(4)、(5)，

b₁＝2R(P_con+I_pU_r+Y_pU_r ²)k_p,max+2X(Q_con+I_qU_r+Y_qU_r ²)-U_s ²    (4)

$c_1=(R^2+X^2)[{(P_{\mathrm{con}}+I_p{U}_r+Y_p{U_r}^2)}^2{k}_{p,\max }^2+{(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U_r}^2)}^2]---\left(5\right)$

被监测节点负荷因子最大值为无功功率负荷因子最大值k_q,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₂(X)为：

$F_2\left(X\right)=\left(\begin{array}{c}{f}_2\left(X\right)\\ g_2\left(X\right)\end{array}\right)---\left(6\right)$

其中：方程组的变量X＝[U_r,k_q,max]^T，U_r为被监测节点电压幅值，k_q,max为被监测节点有功功率负荷因子最大值；

$f_2\left(X\right)=U_r^4+b_2{U}_r^2+c_2---\left(7\right)$

$\begin{array}{c}{g}_2\left(X\right)={4U}_r^3+{2b}_2{U}_r^2+4({\mathrm{RY}}_p+{\mathrm{XY}}_q{k}_{q,\max })U_r^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)k_{q,\max }^2(I_q+{2Y}_q{U}_r)]\end{array}---\left(8\right)$

b₂、c₂的表达式为式(9)、(10)，

b₂＝2R(P_con+I_pU_r+Y_pU_r ²)+2X(Q_con+I_qU_r+Y_qU_r ²)k_q,max-U_s ²    (9)

c₂＝(R²+X²)[(P_con+I_pU_r+Y_pU_r ²)²+(Q_con+I_qU_r+Y_qU_r ²)²k_q,max ²]    (10)

被监测节点负荷因子最大值为视在功率负荷因子最大值k_s,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₃(X)为：

$F_3\left(X\right)=\left(\begin{array}{c}{f}_3\left(X\right)\\ g_3\left(X\right)\end{array}\right)---\left(11\right)$

其中：方程组的变量X＝[U_r,k_s,max]^T，U_r为被监测节点电压幅值，k_s,max为被监测节点有功功率负荷因子最大值；

$f_3\left(X\right)=U_r^4+b_3{U}_r^2+c_3---\left(12\right)$

$\begin{array}{c}{g}_3\left(X\right)={4U}_r^3+{2b}_3{U}_r^2+{4k}_{s,\max }({\mathrm{RY}}_p+{\mathrm{XY}}_q)U_r^3+{2k}_{s,\max }^2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)(I_q+{2Y}_q{U}_r)]\end{array}---\left(13\right)$

b₃、c₃的表达式为式(14)、(15)，

b₂＝2k_s,max[R(P_con+I_pU_r+Y_pU_r ²)+X(Q_con+I_qU_r+Y_qU_r ²)]-U_s ²    (14)

$c_3=(R^2+X^2)k_{s,\max }^2[{(P_{\mathrm{con}}+I_p{U}_r+Y_p{U_r}^2)}^2+{(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U_r}^2)}^2]---\left(15\right)$

B、利用NLEQ-RES算法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0得到k_p,max、k_q,max和k_s,max

B1、得到用于求解被监测节点有功功率负荷因子最大值k_p,max的非线性方程组F(X)＝0

令F(X)＝F₁(X)、f(X)＝f₁(X)、g(X)＝g₁(X)，构造非线性方程组：

$F\left(X\right)=\left(\begin{array}{c}f\left(X\right)\\ g\left(X\right)\end{array}\right)=0---\left(16\right)$

其中：方程组的变量X＝[U_r,k_p,max]^T，U_r为被监测节点电压幅值，k_p,max为被监测节点有功功率负荷因子最大值；

B2、判断第k(k＝0,1,…)次牛顿迭代是否满足收敛条件

如果||F(X^(k))||≤ε，执行B10；其中：ε为任意满足0<ε≤10^-5的常数；

否则，执行B3；

B3、计算第k(k＝0,1,…)次牛顿迭代时F(X^(k))＝0的雅克比矩阵F'(X^(k))

$F^{\&prime;}\left(X^{\left(k\right)}\right)=\left[\begin{array}{cc}{\frac{\&PartialD;f}{{\&PartialD;U}_r}|}_k& {\frac{\&PartialD;f}{{\&PartialD;k}_{p,\max }}|}_k\\ {\frac{\&PartialD;g}{{\&PartialD;U}_r}|}_k& {\frac{\&PartialD;g}{{\&PartialD;k}_{p,\max }}|}_k\end{array}\right]---\left(17\right)$

$\begin{array}{c}{\frac{\&PartialD;f}{{\&PartialD;U}_r}|}_k=4{\left(U_r^{\left(k\right)}\right)}^3+2b{\left(U_r^{\left(k\right)}\right)}^2+4({\mathrm{RY}}_p{k}_{p,\max }^{\left(k\right)}+{\mathrm{XY}}_q){\left(U_r^{\left(k\right)}\right)}^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(k_{p,\max }^{\left(k\right)}\right)}^2(I_p+{2Y}_p{U}_r^{\left(k\right)})\\ +(Q_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_q{\left(U_r^{\left(k\right)}\right)}^2)(I_q+{2Y}_q{U}_r^{\left(k\right)})]\end{array}---\left(18\right)$

$\begin{array}{c}{\frac{\&PartialD;f}{{\&PartialD;k}_{p,\max }}|}_k=2R(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(U_r^{\left(k\right)}\right)}^2\\ +2(R^2+X^2){(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2)}^2{k}_{p,\max }^{\left(k\right)}\end{array}---\left(19\right)$

$\begin{array}{c}{\frac{\&PartialD;g}{{\&PartialD;U}_r}|}_k=12{\left(U_r^{\left(k\right)}\right)}^2+4b{U}_r^{\left(k\right)}+2[2R(I_p+{2Y}_p{U_r}^{\left(k\right)})k_{p,\max }^{\left(k\right)}+2X(I_q+{2Y}_q{U_r}^{\left(k\right)})]{\left(U_r^{\left(k\right)}\right)}^2\\ +12({\mathrm{RY}}_p{k}_{p,\max }^{\left(k\right)}+{\mathrm{XY}}_q){\left(U_r^{\left(k\right)}\right)}^2+2(R^2+X^2)[{(I_p+{2Y}_p{U}_r^{\left(k\right)})}^2{\left(k_{p,\max }^{\left(k\right)}\right)}^2\\ +{2Y}_p(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(k_{p,\max }^{\left(k\right)}\right)}^2+{(I_q+{2Y}_q{U}_r^{\left(k\right)})}^2\\ +{2Y}_q(Q_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_q{\left(U_r^{\left(k\right)}\right)}^2)]\end{array}---\left(20\right)$

$\begin{array}{c}{\frac{\&PartialD;g}{{\&PartialD;k}_{p,\max }}|}_k=4R(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(U_r^{\left(k\right)}\right)}^2+4R{Y}_p{\left(U_r^{\left(k\right)}\right)}^3\\ +4(R^2+X^2){(P_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2)}^2{k}_{p,\max }^{\left(k\right)}\end{array}---\left(21\right)$

B4、计算第k(k＝0,1,…)次牛顿迭代时的修正量ΔX^(k)

ΔX^(k)＝-[F'(X^(k))]^-1F(X^(k))    (22)

B5、计算第k(k＝0,1,…)次牛顿迭代时的阻尼因子估计值[λ_k]

$\left[{\mathrm{\&lambda;}}_k\right]=\left\{\begin{array}{cc}1& k=0\\ \min (\mathrm{1,1}/[h_k\left]\right)& k>0\end{array}\right.---\left(23\right)$

$\begin{array}{cc}\left[h_k\right]=\frac{\left|\right|F\left(X^{\left(k\right)}\right)\left|\right|}{\left|\right|F\left(X^{(k-1)}\right)\left|\right|}\&CenterDot;\left[h_{k-1}^{\&prime;}\right]& k>0---\left(24\right)\end{array}$

B6、根据阻尼因子判断方程组是否奇异

如果[λ_k]≤λ_min，则表明沿着牛顿路径牛顿解逼近于雅克比矩阵奇异点，结束计算，执行B10；其中：λ_min为任意满足0<λ_min≤10^-5的常数；

否则，执行B7；

B7、计算第k(k＝0,1,…)次牛顿迭代时的方程的解X^(k+1)和F(X^(k))

X^(k+1)＝X^(k)+[λ_k]ΔX^(k)    (25)

将X^(k+1)代入函数，计算得到F(X^(k))；

B8、计算监测量[h_k′]

$\left[h_k^{\&prime;}\right]=\frac{2\left|\right|F\left(X^{(k+1)}\right)-(1-[{\mathrm{\&lambda;}}_k\left]\right)F\left(X^{\left(k\right)}\right)\left|\right|}{{\left(\right[{\mathrm{\&lambda;}}_k\left]\right)}^2\left|\right|F\left(X^{\left(k\right)}\right)\left|\right|}---\left(26\right)$

B9、判断第k次牛顿迭代的迭代值X^(k)

如果||F(X^k+1)||≥(1-0.25·[λ_k])||F(X^k)||，

[λ_k′]＝min(0.5·[λ_k],1/[h_k′])    (27)

令[λ_k]＝[λ_k′]，返回B6；

否则，返回B2；

B10、得到被监测节点有功功率负荷因子最大值k_p,max

由B2～B9得到非线性方程组F(X^(k))＝0的解为：X＝[U_r,k_p,max]^T＝X^(k)；

由此得到k_p,max；

B11、利用NLEQ-RES算法求解被监测节点负荷因子最大值k_q,max和k_s,max

令B1中F(X)＝F₂(X)、f(X)＝f₂(X)、g(X)＝g₂(X)，构造非线性方程组F(X)＝0，得到关于k_q,max的方程组式(16)，然后利用方法B2～B10计算得到k_q,max；

令B1中F(X)＝F₃(X)、f(X)＝f₃(X)、g(X)＝g₃(X)，构造非线性方程组F(X)＝0，得到关于k_s,max的方程组式(16)，然后利用方法B2～B10计算得到k_s,max。

C、计算被监测节点计及静态负荷特性的电压稳定指标

C1、计算被监测节点的负荷裕度k_p,margin、k_q,margin和k_s,margin

$\left\{\begin{array}{c}{k}_{p,m\arg \mathrm{in}}=k_{p,\max }-1\\ k_{q,m\arg \mathrm{in}}=k_{q,\max }-1\\ k_{s,m\arg \mathrm{in}}=k_{s,\max }-1\end{array}\right.---\left(28\right)$

C2、计算反映有功功率、无功功率以及视在功率裕度的指标I_VSIP、I_VSIQ和I_VSIS

$\left\{\begin{array}{c}{I}_{\mathrm{VSIP}}=\left|\frac{k_{p,m\arg \mathrm{in}}}{k_{p,\max }}\right|\\ I_{\mathrm{VSIQ}}=\left|\frac{k_{q,m\arg \mathrm{in}}}{k_{q,\max }}\right|\\ I_{\mathrm{VSIS}}=\left|\frac{k_{s,m\arg \mathrm{in}}}{k_{s,\max }}\right|\end{array}\right.---\left(29\right)$

由I_VSIP、I_VSIQ和I_VSIS可以直接得到各负荷裕度。

C3、构造电压稳定指标I_VSI

$I_{\mathrm{VSI}}=\min \left(\right|\frac{k_{p,m\arg \mathrm{in}}}{k_{p,\max }}|,|\frac{k_{q,m\arg \mathrm{in}}}{k_{q,\max }}|,|\frac{k_{s,m\arg \mathrm{in}}}{k_{s,\max }}\left|\right)---\left(30\right)$

D、判断

当I_VSI>0时，系统工作在分岔曲线上半支，系统电压稳定；I_VSI<0时，系统失去电压稳定，进行电压稳定预警；I_VSI＝0时，系统处于临界情况，进行电压稳定预警。

仿真实验一：

采用IEEE50机测试系统对本例的方法进行仿真监测，IEEE50机测试系统包括50台发电机，145各节点，平衡节点为节点145，系统功率基准为100MVA。仿真时，处理节点阻抗负荷为网络中对地阻抗支路，系统中恒功率负荷为节点负荷；潮流计算工具为PSS/E30.0，利用其非解耦牛顿拉夫逊法(FNSL)模块计算节点负荷变化时的潮流。简单传输系统参数采用带遗忘因子的平方根最小二乘估计进行辨识得到，辨识时遗忘因子取值等于0.15，跳数n＝2；求解k_p,max、k_q,max和k_s,max时，λ_min＝10^-5和ε＝10^-5。

在IEEE50机测试系统中，母线51为负荷母线。图2给出了母线51上负荷改变时母线51的电压稳定指标轨迹及计及了负荷静态特性影响的V-Q雅可比矩阵最小特征值轨迹，图中横坐标是负荷阻抗值，纵坐标为稳定指标值，最小特征值轨迹为计及了负荷静态特性影响的V-Q雅可比矩阵最小特征值轨迹；图3中，符号PQ、Z和I分别代表恒功率、恒阻抗和恒电流成份，例如10％Z、10％I和80％PQ表示改变负荷时恒阻抗、恒电流和恒功率成份分别占总负荷的10％、10％和80％。(注：图2和图3中，L₁指标被处理为1减L₁指标，VSLB指标被处理为1减VSLB指标，当系统达到稳定临界点时其理论值等于0，这有助于与提出的电压稳定指标比较)。

图2(a)中母线51上负荷为恒功率负荷，图2(b)中母线51上负荷由80％恒功率成份和20％恒阻抗成份组成；如图2(a)和图2(b)所示，当负荷为恒功率负荷或除恒功率负荷外还含有恒阻抗负荷时，VSLB指标和提出的电压稳定指标都能正确的指示分岔点，而L₁指标可能有大的误差。图2(c)中母线51上负荷由50％恒功率成份和50％恒电流成份组成，图2(d)中母线51上负荷由80％恒功率成份、10％恒阻抗成份和10％恒电流成份组成；如图2(c)和图2(d)所示，提出的电压稳定指标都能正确的指示分岔点，但如果负荷中含有恒电流成份较大时，则VSLB指标可能会有大的误差。

表1为在IEEE50机测试系统中母线51含不同负荷模型时鞍结分岔点附近母线51的电压稳定指标比较，进一步比较了改变母线51负荷逼近分岔点附近时母线51的几种典型的计及了负荷静态特性的电压稳定指标值的精度，λ_min是计及了负荷静态特性影响的扩展雅克比矩阵最小特征值，符号PQ、Z和I的含义与图2一致。如表1所示，基于扩展雅克比矩阵行列式的指标值和L₁指标值有大的偏差，VSLB指标精度受负荷恒电流成份比例的影响较大，而提出的电压稳定指标都能正确的反映分岔点。

表1

仿真实验二：

在IEEE50机测试系统中，母线120为负荷母线。图3给出了母线120上负荷改变时母线120的电压稳定指标轨迹以及计及了负荷静态特性影响的V-Q雅可比矩阵最小特征值轨迹，图中横坐标、纵坐标、各符号以及曲线类型定义都与图3相同。

图3(a)中母线120上负荷为恒功率负荷，图3(b)中母线120上负荷由80％恒功率成份和20％恒阻抗成份组成；如图3(a)和图3(b)所示，当负荷为恒功率负荷或除恒功率负荷外还含有恒阻抗负荷时，VSLB指标和提出的电压稳定指标都能正确的指示分岔点，而L₁指标受负荷恒阻抗成份的影响出现了大的误差。图3(c)中母线51上负荷由50％恒功率成份和50％恒电流成份组成，图3(d)中母线51上负荷由60％恒功率成份、20％恒阻抗成份和20％恒电流成份组成；如图3(c)和图3(d)所示，如果负荷中含有恒电流成份较大时，VSLB指标可能会有大的误差，且L₁指标受负荷恒电流成份的影响也出现了大的误差，但提出的电压稳定指标都正确的指示了分岔点。

而对于基于扩展雅克比矩阵行列式的指标，由于受雅克比矩阵中其他特征值的影响，在逼近分岔点时，它也有较大误差。例如，在母线120上负荷由60％恒功率成份、20％恒阻抗成份和20％恒电流成份组成的情况下，当系统扩展雅克比矩阵的最小特征值λ_min等于0.0079285时，对于母线120节点，而对应基于扩展雅克比矩阵行列式的指标值等于2.1245，提出的电压稳定指标值等于0.000001。

Claims (5)

1.一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，其特征在于，包括以下步骤：

A、根据含负荷节点简单传输系统达到鞍结分岔点条件建立非线性函数F₁(X)、F₂(X)、F₃(X)

基于已建立的含负荷节点简单传输系统的模型及参数，根据该系统达到鞍结分岔点的特征方程建立三个非线性函数F₁(X)、F₂(X)、F₃(X)，各函数的变量X分别为向量[U_r,k_p,max]的转置[U_r,k_p,max]^T、[U_r,k_q,max]的转置[U_r,k_q,max]^T和[U_r,k_s,max]的转置[U_r,k_s,max]^T，U_r为被监测节点电压幅值，k_p,max、k_q,max和k_s,max分别为被监测节点有功功率、无功功率、视在功率负荷因子最大值；

B、利用NLEQ-RES算法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0得到k_p,max、k_q,max和k_s,max

根据A中建立的非线性函数F₁(X)、F₂(X)、F₃(X)，采用具有全局收敛且能自适应最优迭代步长的基于收敛准则和自适应置信区间策略的带残差全局牛顿法即NLEQ-RES法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0，得到被监测节点有功功率、无功功率、视在功率负荷因子最大值k_p,max、k_q,max和k_s,max，每次迭代求解时均自适应调整迭代步长；

C、计算被监测节点计及静态负荷特性的电压稳定指标

根据B中计算得到的负荷因子最大值k_p,max、k_q,max和k_s,max，得到被监测节点的负荷裕度k_p,margin、k_q,margin和k_s,margin，计算反映有功功率、无功功率以及视在功率裕度的指标I_VSIP、I_VSIQ和I_VSIS，基于此构造并计算被监测节点计及静态负荷特性的电压稳定指标I_VSI；

D、判断

根据C步计算出的被监测节点的电压稳定指标，判断被监测节点是否为薄弱节点；如果是薄弱节点，则进行电压稳定预警，否则不作预警；判断结束后，返回步骤A继续执行下一时刻的电压稳定监测。

2.如权利要求1所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，其特征在于：所述A步根据含负荷节点简单传输系统达到鞍结分岔点条件建立非线性函数F₁(X)、F₂(X)、F₃(X)的具体步骤为：

已建立的含负荷的被监测节点r处的简单传输系统参数如下：U_s∠δ为电源端电压，U_r∠0为负荷侧电压，Z＝R+jX为线路阻抗，Y_L＝Y_p+jY_q为ZIP静态负荷中恒阻抗负荷分量的导纳值，I_L＝I_p+jI_q为ZIP静态负荷中恒电流负荷分量的电流值，S_con＝P_con+jQ_con为ZIP静态负荷中恒功率负荷分量的功率值；

被监测节点负荷因子最大值为有功功率负荷因子最大值k_p,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₁(X)为：

$F_1\left(X\right)=\left(\begin{array}{c}{f}_1\left(X\right)\\ g_1\left(X\right)\end{array}\right)---\left(1\right)$

其中：方程组的变量X＝[U_r,k_p,max]^T，U_r为被监测节点电压幅值，k_p,max为被监测节点有功功率负荷因子最大值；

$f_1\left(X\right)=U_r^4+b_1{U}_r^2+c_1---\left(2\right)$

$\begin{array}{c}{g}_1\left(X\right)=4U_r^3+{2b}_1{U}_r^2+4({\mathrm{RY}}_p{k}_{p,\max }+{\mathrm{XY}}_q)U_r^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)k_{p,\max }^2(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)(I_q+{2Y}_q{U}_r)]\end{array}---\left(3\right)$

b₁、c₁的表达式为式(4)、(5)，

b₁＝2R(P_con+I_pU_r+Y_pU_r ²)k_p,max+2X(Q_con+I_qU_r+Y_qU_r ²)-U_s ²    (4)

$c_1=(R^2+X^2)[{(P_{\mathrm{con}}+I_p{U}_r+Y_p{U_r}^2)}^2{k}_{p,\max }^2+{(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U_r}^2)}^2]---\left(5\right)$

被监测节点负荷因子最大值为无功功率负荷因子最大值k_q,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₂(X)为：

$F_2\left(X\right)=\left(\begin{array}{c}{f}_2\left(X\right)\\ g_2\left(X\right)\end{array}\right)---\left(6\right)$

其中：方程组的变量X＝[U_r,k_q,max]^T，U_r为被监测节点电压幅值，k_q,max为被监测节点有功功率负荷因子最大值；

$f_2\left(X\right)=U_r^4+b_2{U}_r^2+c_2---\left(7\right)$

$\begin{array}{c}{g}_2\left(X\right)={4U}_r^3+{2b}_2{U}_r^2+4({\mathrm{RY}}_p+{\mathrm{XY}}_q{k}_{q,\max })U_r^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)k_{q,\max }^2(I_q+{2Y}_q{U}_r)]\end{array}---\left(8\right)$

b₂、c₂的表达式为式(9)、(10)，

b₂＝2R(P_con+I_pU_r+Y_pU_r ²)+2X(Q_con+I_qU_r+Y_qU_r ²)k_q,max-U_s ²    (9)

c₂＝(R²+X²)[(P_con+I_pU_r+Y_pU_r ²)²+(Q_con+I_qU_r+Y_qU_r ²)²k_q,max ²]    (10)

被监测节点负荷因子最大值为视在功率负荷因子最大值k_s,max时，根据含负荷节点简单传输系统达到鞍结分岔点条件建立的非线性函数F₃(X)为：

$F_3\left(X\right)=\left(\begin{array}{c}{f}_3\left(X\right)\\ g_3\left(X\right)\end{array}\right)---\left(11\right)$

其中：方程组的变量X＝[U_r,k_s,max]^T，U_r为被监测节点电压幅值，k_s,max为被监测节点有功功率负荷因子最大值；

$f_3\left(X\right)=U_r^4+b_3{U}_r^2+c_3---\left(12\right)$

$\begin{array}{c}{g}_3\left(X\right)={4U}_r^3+{2b}_3{U}_r^2+{4k}_{s,\max }({\mathrm{RY}}_p+{\mathrm{XY}}_q)U_r^3+{2k}_{s,\max }^2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r+Y_p{U}_r^2)(I_p+{2Y}_p{U}_r)+(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U}_r^2)(I_q+{2Y}_q{U}_r)]\end{array}---\left(13\right)$

b₃、c₃的表达式为式(14)、(15)，

b₂＝2k_s,max[R(P_con+I_pU_r+Y_pU_r ²)+X(Q_con+I_qU_r+Y_qU_r ²)]-U_s ²    (14)

$c_3=(R^2+X^2)k_{s,\max }^2[{(P_{\mathrm{con}}+I_p{U}_r+Y_p{U_r}^2)}^2+{(Q_{\mathrm{con}}+I_q{U}_r+Y_q{U_r}^2)}^2]---\left(15\right)$

3.如权利要求1所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，其特征在于：所述B步利用NLEQ-RES算法求解方程组F₁(X)＝0、F₂(X)＝0、F₃(X)＝0得到k_p,max、k_q,max和k_s,max的步骤为：

B1、得到用于求解被监测节点有功功率负荷因子最大值k_p,max的非线性方程组F(X)＝0

令F(X)＝F₁(X)、f(X)＝f₁(X)、g(X)＝g₁(X)，构造非线性方程组：

$F\left(X\right)=\left(\begin{array}{c}f\left(X\right)\\ g\left(X\right)\end{array}\right)=0---\left(16\right)$

其中：方程组的变量X＝[U_r,k_p,max]^T，U_r为被监测节点电压幅值，k_p,max为被监测节点有功功率负荷因子最大值；

B2、判断第k(k＝0,1,…)次牛顿迭代是否满足收敛条件

如果||F(X^(k))||≤ε，执行B10；其中：ε为任意满足0<ε≤10^-5的常数；

否则，执行B3；

B3、计算第k(k＝0,1,…)次牛顿迭代时F(X^(k))＝0的雅克比矩阵F'(X^(k))

$F^{\&prime;}\left(X^{\left(k\right)}\right)=\left[\begin{array}{cc}{\frac{\&PartialD;f}{{\&PartialD;U}_r}|}_k& {\frac{\&PartialD;f}{{\&PartialD;k}_{p,\max }}|}_k\\ {\frac{\&PartialD;g}{{\&PartialD;U}_r}|}_k& {\frac{\&PartialD;g}{{\&PartialD;k}_{p,\max }}|}_k\end{array}\right]---\left(17\right)$

$\begin{array}{c}{\frac{\&PartialD;f}{{\&PartialD;U}_r}|}_k=4{\left(U_r^{\left(k\right)}\right)}^3+2b{\left(U_r^{\left(k\right)}\right)}^2+4({\mathrm{RY}}_p{k}_{p,\max }^{\left(k\right)}+{\mathrm{XY}}_q){\left(U_r^{\left(k\right)}\right)}^3+2(R^2+X^2)\&CenterDot;\\ [(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(k_{p,\max }^{\left(k\right)}\right)}^2(I_p+{2Y}_p{U}_r^{\left(k\right)})\\ +(Q_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_q{\left(U_r^{\left(k\right)}\right)}^2)(I_q+{2Y}_q{U}_r^{\left(k\right)})]\end{array}---\left(18\right)$

$\begin{array}{c}{\frac{\&PartialD;f}{{\&PartialD;k}_{p,\max }}|}_k=2R(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(U_r^{\left(k\right)}\right)}^2\\ +2(R^2+X^2){(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2)}^2{k}_{p,\max }^{\left(k\right)}\end{array}---\left(19\right)$

$\begin{array}{c}{\frac{\&PartialD;g}{{\&PartialD;U}_r}|}_k=12{\left(U_r^{\left(k\right)}\right)}^2+4b{U}_r^{\left(k\right)}+2[2R(I_p+{2Y}_p{U_r}^{\left(k\right)})k_{p,\max }^{\left(k\right)}+2X(I_q+{2Y}_q{U_r}^{\left(k\right)})]{\left(U_r^{\left(k\right)}\right)}^2\\ +12({\mathrm{RY}}_p{k}_{p,\max }^{\left(k\right)}+{\mathrm{XY}}_q){\left(U_r^{\left(k\right)}\right)}^2+2(R^2+X^2)[{(I_p+{2Y}_p{U}_r^{\left(k\right)})}^2{\left(k_{p,\max }^{\left(k\right)}\right)}^2\\ +{2Y}_p(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(k_{p,\max }^{\left(k\right)}\right)}^2+{(I_q+{2Y}_q{U}_r^{\left(k\right)})}^2\\ +{2Y}_q(Q_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_q{\left(U_r^{\left(k\right)}\right)}^2)]\end{array}---\left(20\right)$

$\begin{array}{c}{\frac{\&PartialD;g}{{\&PartialD;k}_{p,\max }}|}_k=4R(P_{\mathrm{con}}+I_p{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2){\left(U_r^{\left(k\right)}\right)}^2+4R{Y}_p{\left(U_r^{\left(k\right)}\right)}^3\\ +4(R^2+X^2){(P_{\mathrm{con}}+I_q{U}_r^{\left(k\right)}+Y_p{\left(U_r^{\left(k\right)}\right)}^2)}^2{k}_{p,\max }^{\left(k\right)}\end{array}---\left(21\right)$

B4、计算第k(k＝0,1,…)次牛顿迭代时的修正量ΔX^(k)

ΔX^(k)＝-[F'(X^(k))]^-1F(X^(k))    (22)

B5、计算第k(k＝0,1,…)次牛顿迭代时的阻尼因子估计值[λ_k]

$\left[{\mathrm{\&lambda;}}_k\right]=\left\{\begin{array}{cc}1& k=0\\ \min (\mathrm{1,1}/[h_k\left]\right)& k>0\end{array}\right.---\left(23\right)$

$\begin{array}{cc}\left[h_k\right]=\frac{\left|\right|F\left(X^{\left(k\right)}\right)\left|\right|}{\left|\right|F\left(X^{(k-1)}\right)\left|\right|}\&CenterDot;\left[h_{k-1}^{\&prime;}\right]& k>0---\left(24\right)\end{array}$

B6、根据阻尼因子判断方程组是否奇异

如果[λ_k]≤λ_min，则表明沿着牛顿路径牛顿解逼近于雅克比矩阵奇异点，结束计算，执行B10；其中：λ_min为任意满足0<λ_min≤10^-5的常数；

否则，执行B7；

B7、计算第k(k＝0,1,…)次牛顿迭代时的方程的解X^(k+1)和F(X^(k))

X^(k+1)＝X^(k)+[λ_k]ΔX^(k)    (25)

将X^(k+1)代入函数，计算得到F(X^(k))；

B8、计算监测量[h_k′]

$\left[h_k^{\&prime;}\right]=\frac{2\left|\right|F\left(X^{(k+1)}\right)-(1-[{\mathrm{\&lambda;}}_k\left]\right)F\left(X^{\left(k\right)}\right)\left|\right|}{{\left(\right[{\mathrm{\&lambda;}}_k\left]\right)}^2\left|\right|F\left(X^{\left(k\right)}\right)\left|\right|}---\left(26\right)$

B9、判断第k次牛顿迭代的迭代值X^(k)

如果||F(X^k+1)||≥(1-0.25·[λ_k])||F(X^k)||，

[λ_k′]＝min(0.5·[λ_k],1/[h_k′])    (27)

令[λ_k]＝[λ_k′]，返回B6；

否则，返回B2；

B10、得到被监测节点有功功率负荷因子最大值k_p,max

由B2～B9得到非线性方程组F(X^(k))＝0的解为：X＝[U_r,k_p,max]^T＝X^(k)；

由此得到k_p,max；

B11、利用NLEQ-RES算法求解被监测节点负荷因子最大值k_q,max和k_s,max

令B1中F(X)＝F₂(X)、f(X)＝f₂(X)、g(X)＝g₂(X)，构造非线性方程组F(X)＝0，得到关于k_q,max的方程组式(16)，然后利用方法B2～B10计算得到k_q,max；

令B1中F(X)＝F₃(X)、f(X)＝f₃(X)、g(X)＝g₃(X)，构造非线性方程组F(X)＝0，得到关于k_s,max的方程组式(16)，然后利用步骤B2～B10计算得到k_s,max。

4.如权利要求1所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，其特征在于：所述C步计算被监测节点计及静态负荷特性的电压稳定指标的步骤为：

C1、计算被监测节点的负荷裕度k_p,margin、k_q,margin和k_s,margin

$\left\{\begin{array}{c}{k}_{p,m\arg \mathrm{in}}=k_{p,\max }-1\\ k_{q,m\arg \mathrm{in}}=k_{q,\max }-1\\ k_{s,m\arg \mathrm{in}}=k_{s,\max }-1\end{array}\right.---\left(28\right)$

C2、计算反映有功功率、无功功率以及视在功率裕度的指标I_VSIP、I_VSIQ和I_VSIS

$\left\{\begin{array}{c}{I}_{\mathrm{VSIP}}=\left|\frac{k_{p,m\arg \mathrm{in}}}{k_{p,\max }}\right|\\ I_{\mathrm{VSIQ}}=\left|\frac{k_{q,m\arg \mathrm{in}}}{k_{q,\max }}\right|\\ I_{\mathrm{VSIS}}=\left|\frac{k_{s,m\arg \mathrm{in}}}{k_{s,\max }}\right|\end{array}\right.---\left(29\right)$

由I_VSIP、I_VSIQ和I_VSIS可以直接得到各负荷裕度

C3、构造电压稳定指标I_VSI

$I_{\mathrm{VSI}}=\min \left(\right|\frac{k_{p,m\arg \mathrm{in}}}{k_{p,\max }}|,|\frac{k_{q,m\arg \mathrm{in}}}{k_{q,\max }}|,|\frac{k_{s,m\arg \mathrm{in}}}{k_{s,\max }}\left|\right)---\left(30\right)$

5.如权利要求1所述的一种基于计及静态负荷特性的节点电压稳定指标的电网电压稳定监测方法，其特征在于：所述D步判断的具体方法为：

当I_VSI>0时，系统工作在分岔曲线上半支，系统电压稳定；I_VSI<0时，系统失去电压稳定，进行电压稳定预警；I_VSI＝0时，系统处于临界情况，进行电压稳定预警。