CN103532137A - 一种三相四线低压配电网的状态估计方法 - Google Patents

Abstract

一种三相四线低压配电网状态估计方法，属于配电网的三相状态估计技术领域。本发明方法利用计算机，通过程序，首先输入采集的任一时间断面的低压配抄表电力数据和网络结构及参数信息并初始化，然后计算节点导纳矩阵和雅可比矩阵常数项部分，接着综合量测方程和T端点零注入等式约束，以电压的实虚部为状态变量，计算不平衡量、指数型权函数及雅可比矩阵，最后更新状态变量，并收敛性判断，实现三相四线低压配电网状态估计。本发明能独立计算节点有功与无功功率和抑制坏数据对估计结果的不良影响，具有强抗差性和良好的收敛性及数值稳定性等特点。本发明方法广泛应用于求解三相四线低压配电网的状态估计。

Description

一种三相四线低压配电网的状态估计方法

技术领域

本发明属于配电网的三相状态估计技术领域，具体涉及一种三相四线低压配电网状态估计方法。

背景技术

随着经济的发展和人民生活水平的提高，低压配电网电力用户对供电质量提出了越来越高的要求，因此加强对低压配电网的运行监视与控制显得尤为重要。状态估计是能量管理系统的重要组成部分，其结果直接影响电网调度的智能化分析与决策，并为配电管理系统提供可靠而完整的系统运行状态信息，是配电系统分析的数据之源，对配电网安全可靠、优质运行发挥重大作用，有重要的研究价值。

对三相四线制低压配电网，长期以来，由于缺少用户电力的实测数据信息，无法进行低压配电网的潮流计算及线损分析，低压配电网状态估计也存在数据不足而实现困难等问题。近年来，智能电网中用户智能电表和配变终端信息采集系统在380/220V的低压配电网中大量推广应用，从而能够同时采集三相电力与电量信息，为低压配电网的用电管理提供充裕的数据支撑，并为低压配电网的状态估计及其基础上的用电监察提供了可能。但至今未见三相四线低压配电网状态估计的研究报道。

现有的配电网状态估计方法，如2011年第31卷第19期《中国电机工程学报》中“支持大规模电流量测的配网抗差状态估计方法”一文，公开的方法是以支路首端功率和支路电流幅值的平方作为待求量，采用指数型状态估计目标函数的抗差估计模型。该方法综合利用了电压幅值量测、功率量测和电流量测，能有效处理辐射电网和弱环网的三相状态估计问题。但该方法的主要缺点是：①要求功率量测成对出现且P、Q量测的权系数相同，影响状态估计结果精度；②等效复电流量测在迭代过程中变化较大，收敛性差，导致状态估计的计算效率低；③主要适用于以支路电流幅值量测为主的三相三线中、高压配网。在三相四线低压配电网状态估计中，由于线制、量测量、量测方程和约束方程等与现有的三相三线接地配电网状态估计有较大差异，因此方法只能对辐射电网和弱环网的三相三线制配电网进行状态评估，不能对三相四线低压配电网进行状态估计。

发明内容

本发明的目的是针对现有配电网状态估计研究方法的不足，提供了一种三相四线低压配电网的状态估计方法，通过充分利用量测值和合理分配权值，该方法具有强抗差性和良好的收敛性、收敛速度及在线应用前景等特点，且不要求功率量测成对出现以及P、Q量测的权系数相同。本方法充分利用智能电表的电力量测数据，能有效求解三相四线低压配电网的状态估计问题，同时有效抑制坏数据和量测误差对估计结果的影响，改善了三相四线低压配电网状态估计的收敛性能。

实现本发明目的之技术方案是：一种三相四线低压配电网状态估计方法，利用计算机，通过程序，首先输入数据终端采集到的任一时间断面的低压配抄表电力数据和网络结构及参数信息并初始化，然后计算节点导纳矩阵和雅可比矩阵常数项部分，接着综合考虑相对中性点的注入电流量测方程、中性点的虚拟量测方程、相对中性点的电压幅值量测方程和T端点的零注入等式约束方程，以相对地电压和中性点电压的实虚部为状态变量，应用指数型权函数，计算相应的不平衡量、指数型权函数及雅可比矩阵，最后更新状态变量，进行收敛性判断，来实现三相四线低压配电网的状态估计。所述方法的具体步骤如下：

(1)输入基础数据及初始化

1)输入基础数据

首先输入任一时间断面的低压配总表及各用户分表的瞬时抄表电力数据和网络结构及参数信息。即输入低压配总表在任一时间断面下的瞬时抄表总有功功率和总无功功率以及三相电压幅值和三相电流幅值，各用户分表在同一时间断面下的瞬时抄表各相有功功率和各相无功功率数据以及三相电压幅值，网络结构及参数信息为网络线路电阻、电抗和电纳参数、线路的额定电压、功率基准。

2)参数初始化

设置n阶单位矩阵R^-1（即n阶单位矩阵R^-1的对角线元素全为1，非对角线元素全为0），n为状态估计中实际的量测变量个数；初始化最大迭代次数Tmax为40～60；收敛精度ε为10^-3～10^-5及残差方差的检测门槛值γ=5，并设置迭代次数time=1。

(2)形成节点导纳矩阵和雅可比矩阵常数项部分

第(1)步完成后，计算节点导纳矩阵，计算公式为：

式中：B₁＝{a,b,c,n}，表示端点包含中性点n的电气节点组合，其中a、b、c分别表示abc三相电气节点；x、d、t表示电气节点组合B₁中的任一节点；φ_i为与端点i直接相联但不包括端点i的端点集合；j为端点集合φ_i中的任一端点；

(x∈B₁&x≠d)表示端点i中节点d和x之间的并联支路导纳；

表示端点i中节点d与端点j中节点t直接相联的支路导纳；表示端点i中节点d对地电容得到的对地并联导纳；自导纳表示与端点i节点d直接相联的所有支路导纳和对地电容导纳

之和，互导纳

(t≠d)表示端点i中节点d和t之间支路导纳

的相反数与端点i且与i直接相联端点节点d和t之间所有支路导纳之和，互导纳

(k≠i)则表示端点i和k中节点d和t之间支路导纳的相反数。

得到节点导纳矩阵后，计算雅可比矩阵常数项子矩阵，计算公式为：

$\left\{\begin{array}{c}{H}_{\mathrm{ij}}(1:\mathrm{8,1}:8)={\left[\begin{array}{cc}{G}_{\mathrm{ij}}^{\mathrm{dt}}& -B_{\mathrm{ij}}^{\mathrm{dt}}\\ B_{\mathrm{ij}}^{\mathrm{dt}}& G_{\mathrm{ij}}^{\mathrm{dt}}\end{array}\right]}_{8\×8},d\∈B_1,t\∈B_1\\ H_{\mathrm{ij}}(9:\mathrm{11,1}:8)={\left[0\right]}_{3\×8}\end{array}\right.---\left(2\right)$

C_ij(1:8,1:8)＝H_ij(1:8,1:8)    (3)

式中：i和j为端点编号；H_ij和C_ij分别为雅可比矩阵H和C的常数项子矩阵；

和

分别为节点导纳矩阵元素

的实部和虚部；[0]表示全零矩阵，下标3×8表示该矩阵为3行8列；d、t及B₁的意义同公式(1)。

(3)计算不平衡量和指数型权函数及雅可比矩阵

第(2)步完成后，计算三相四线低压配电网中量测量的不平衡量和指数型权函数及雅可比矩阵，具体步骤如下：

1)计算不平衡量

基于公式(4)计算三相四线低压配电网总表端点s的电流不平衡量和有功功率不平衡量及无功功率不平衡量，计算公式为：

式中：B_P＝{a,b,c}，表示端点不包含中性点n的电气节点组合，其中a、b、c分别表示abc三相电气节点；d表示电气节点组合B_P中的任一节点；B₁、t的意义同公式(1)；s表示三相四线低压配电网总表的端点；

为包括总表端点s且与其直接相联端点的集合；k为端点集合

中的任一端点；Δz_s为基于公式(4)计算端点s的电流不平衡量和功率不平衡量；(d∈B_P)、

和

分别为低压配电网总表s相对中性点n的三相电流幅值、三相总注入有功和三相总注入无功的瞬时抄表量测量；

和分别为端点k中节点t的电压

的实部和虚部；

和

分别为端点s中节点d的电压

的实部和虚部；和

分别为端点s的中性点n的电压的实部和虚部；

和

分别为节点导纳矩阵元素的实部和虚部。

基于公式(5)计算三相四线低压配电网用户分表端点i中a、b、c三相节点注入有功功率及注入无功功率的不平衡量，计算公式为：

式中：

和

分别为低压配电网用户分表i中节点d相对中性点n的注入有功和注入无功的瞬时抄表量测量。i表示三相四线低压配电网中任意某个用户分表端点；

和

(d∈B_p)分别为端点i中节点d的电压

的实部和虚部；

和

分别为端点i的中性点n电压的实部和虚部；Δz_i为基于公式(5)计算端点i的节点注入功率的不平衡量；

和分别为节点导纳矩阵元素

的实部和虚部；

为包括端点i且与其直接相联端点的集合；k为端点集合

中的任一端点；

B_P、d、B₁及t的意义同公式(4)。

基于公式(6)计算三相四线低压配电网用户分表端点i的中性点注入有功功率及注入无功功率的不平衡量，基于公式(7)计算三相四线低压配电网用户分表端点i中a、b、c三相节点电压的不平衡量，计算公式分别为：

$\mathrm{\Δ}{z}_i={\left({\hat{U}}_i^{\mathrm{dn}}\right)}^2-({(e_i^d-e_i^n)}^2+{(f_i^d-f_i^n)}^2),d\∈B_P---\left(7\right)$

式中：

为低压配电网分表i相对中性点n的电压幅值的瞬时抄表量测量；

和

分别为节点导纳矩阵元素的实部和虚部；Δz_i为用户分表端点i计算得到的功率和电压不平衡量；

k、B_P、d、B₁及t的意义同公式(5)。

基于公式(8)计算零注入电流端点m中a、b、c三相及中性点注入电流的不平衡量Δc_m，其计算公式为：

式中：

和

分别为节点导纳矩阵元素

的实部和虚部；

为包括T端点m且与其直接相联端点的集合；k为端点集合

中的任一端点；的意义同公式(4)；B₁、t、d的意义同公式(1)。

2)计算指数型权函数

基于公式(9)-公式(10)计算指数型权函数W的对角阵元素，计算公式为：

$w_i^{*}=R_i^{-1}{e}^{{-r}_{\mathrm{Ni}}^2/{2\mathrm{\δ}}^2}---\left(9\right)$

式中：w_i ^*为量测i的权函数；

为量测i的固定权重；δ为概率密度函数的标准差，其初值由δ²＝r^TPr/(m-n)计算得到，r为残差向量，r^T为残差向量的转置，P为先验矩阵

m为量测量的个数，n为状态变量的个数。r_Ni为量测i的标准化残差，计算公式为

$r_{\mathrm{Ni}}=\left\{\begin{array}{cc}\left|r_i\right|/\sqrt{\left|{\left(\mathrm{\ΔR}\right)}_i\right|},& \left|{\left(\mathrm{\ΔR}\right)}_i\right|>\mathrm{\γ}\\ \left|r_i\right|/\sqrt{\mathrm{\γ}},& \left|{\left(\mathrm{\ΔR}\right)}_i\right|\≤\mathrm{\γ}\end{array}\right.---\left(10\right)$

式中：残差方差的检测门槛值γ通过试探法或经验值选取；

为量测i对应的残差；R为测量误差方差的对角阵，对角线元素为测量误差的方差；Δ＝Ι-H(H^TWH)^-1H^TW；Ι为n阶单位阵；W为指数型权函数对角阵，其对角元素等于权函数，即W_ii＝w_i ^*；为雅可比矩阵，H^T为雅克比矩阵的转置；(H^TWH)^-1表示矩阵H^TWH相乘后的逆矩阵。

3)形成雅可比矩阵

基于公式(11)-公式(15)形成雅可比矩阵H和C，根据公式(16)计算A和B。

对于平衡节点，基于公式(11)和(12)形成其雅克比矩阵元素，计算公式为：

式中：雅可比子矩阵H_ss＝[H_ss1(1:3,1:3);H_ss2(1:3,1:3);H_ss3(1,1:3);H_ss4(1,1:3)；]和

分别为节点导纳矩阵元素

和

的实部，

和

分别为节点导纳矩阵元素和的虚部；

和

分别为端点f中节点t电压的实部和虚部，f为端点集合

中的任一端点；

和

分别为端点s中的节点p电压的实部和虚部；

B_P、B₁、d、t、s的意义同公式(4)；p表示电气节点组合B_P中的任一节点；

为平衡端点s中节点p电压虚部对平衡端点s中节点p电压实部的偏导数，且 $\frac{{\∂f}_s^a}{{\∂e}_s^a}=0,\frac{{\∂f}_s^b}{{\∂e}_s^b}=\sqrt{3},$ $\frac{{\∂f}_s^c}{{\∂e}_s^c}=-\sqrt{3}.$

$\left\{\begin{array}{c}{H}_{\mathrm{sj}1}(1:\mathrm{3,1}:8)={\left[\begin{array}{c}2G_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t-B_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t)+2B_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t+B_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t)\\ -2B_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t-B_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t)+2G_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t+B_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t)\end{array}\right]}_{3\×8}^T\\ H_{\mathrm{sj}2}(1:\mathrm{3,1}:8)={\left[0\right]}_{3\×8}\\ H_{\mathrm{sj}3}(\mathrm{1,1}:8)={\left[\begin{array}{c}\underset{d\∈B_P}{\mathrm{\Σ}}\{G_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)+B_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\\ \underset{d\∈B_P}{\mathrm{\Σ}}\{-B_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)+G_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\end{array}\right]}_{1\×8}^T\\ H_{\mathrm{sj}4}(\mathrm{1,1}:8)={\left[\begin{array}{c}\underset{d\∈B_P}{\mathrm{\Σ}}\{-B_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)+G_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\\ \underset{d\∈B_P}{\mathrm{\Σ}}\{-G_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)-B_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\end{array}\right]}_{1\×8}^T\end{array}\right.---\left(12\right)$

式中：雅可比子矩阵H_sj＝[H_sj1(1:3,1:8);H_sj2(1:3,1:8);H_sj3(1,1:8);H_sj4(1,1:8)]；s为平衡端点；j为非平衡端点；

和

分别为节点导纳矩阵元素

的实部和虚部；

和

分别为端点j中节点t电压的实部和虚部；B_P、

d及t的意义同公式(4)。

对于非平衡节点，基于公式(13)～(15)形成其雅克比矩阵元素，计算公式为：

$\left\{\begin{array}{c}{H}_{\mathrm{is}}(1:\mathrm{8,1}:3)={[G_{\mathrm{is}}^{\mathrm{dp}}-B_{\mathrm{is}}^{\mathrm{dp}}\frac{{\∂f}_s^p}{{\∂e}_s^p};B_{\mathrm{is}}^{\mathrm{dp}}+G_{\mathrm{is}}^{\mathrm{dp}}\frac{{\∂f}_s^p}{{\∂e}_s^p}]}_{8\×3},d\∈B_1,p\∈B_P\\ H_{\mathrm{is}}(9:11,1:3)={\left[0\right]}_{3\×3}\end{array}\right.---\left(13\right)$

式中：

和

分别为节点导纳矩阵元素

的实部和虚部；

B_P、B₁、p的意义同公式(11)；d表示电气节点组合B₁中的任一节点。

$\left\{\begin{array}{c}{\mathrm{\ΔH}}_{\mathrm{ii}}(1:\mathrm{3,1}:3)=-\mathrm{\Δ}{H}_{\mathrm{ii}}(5:\mathrm{7,5}:7)={\left[{\mathrm{\μ}}_i^d\right]}_{3\×3}\\ \mathrm{\Δ}{H}_{\mathrm{ii}}(5:\mathrm{7,8})=-\mathrm{\Δ}{H}_{\mathrm{ii}}(1:\mathrm{3,4})={\left[{\mathrm{\μ}}_i^d\right]}_{3\×1}\\ \mathrm{\Δ}{H}_{\mathrm{ii}}(\mathrm{8,5}:7)=-\mathrm{\Δ}{H}_{\mathrm{ii}}(\mathrm{4,1}:3)={\left[{\mathrm{\μ}}_i^d\right]}_{1\×3}\\ {\mathrm{\ΔH}}_{\mathrm{ii}}\left(\mathrm{4,4}\right)=-{\mathrm{\ΔH}}_{\mathrm{ii}}\left(\mathrm{8,8}\right)=\underset{d\∈B_P}{\mathrm{\Σ}}{\mathrm{\μ}}_i^d\\ {\mathrm{\ΔH}}_{\mathrm{ii}}(1:3,5:7)=\mathrm{\Δ}{H}_{\mathrm{ii}}(5:7,1:3)={[-{\mathrm{\η}}_i^d]}_{3\×3}\\ \mathrm{\Δ}{H}_{\mathrm{ii}}(1:3,8)=\mathrm{\Δ}{H}_{\mathrm{ii}}(5:7,4)={\left[{\mathrm{\η}}_i^d\right]}_{3\×1}\\ \mathrm{\Δ}{H}_{\mathrm{ii}}(4,5:7)=\mathrm{\Δ}{H}_{\mathrm{ii}}(8,1:3)={\left[{\mathrm{\η}}_i^d\right]}_{1\×3}\\ {\mathrm{\ΔH}}_{\mathrm{ii}}(4,8)={\mathrm{\ΔH}}_{\mathrm{ii}}(8,4)=-\underset{d\∈B_P}{\mathrm{\Σ}}{\mathrm{\η}}_i^d\\ {\mathrm{\ΔH}}_{\mathrm{ii}}(9:11,1:3)={\left[2(e_i^d-e_i^n)\right]}_{3\×3}\\ {\mathrm{\ΔH}}_{\mathrm{ii}}(9:\mathrm{11,4})={[-2(e_i^d-e_i^n)]}_{3\×1}\\ {\mathrm{\ΔH}}_{\mathrm{ii}}(9:\mathrm{11,5}:7)={\left[2(f_i^d-f_i^n)\right]}_{3\×3}\\ {\mathrm{\ΔH}}_{\mathrm{ii}}(9:\mathrm{11,8})={[-2(f_i^d-f_i^n)]}_{3\×1}\end{array}\right.---\left(14\right)$

式中：ΔH_ii为雅可比子矩阵H_ii的修正项子阵，i为非平衡端点；B_P、d、

及

的意义同公式(5)；

和

均为修正项子阵ΔH_ii的矩阵元素，其计算公式分别为：

$\left\{\begin{array}{c}{\mathrm{\η}}_i^d=\frac{{\hat{Q}}_i^{\mathrm{dn}}({(e_i^d-e_i^n)}^2-{(f_i^d-f_i^n)}^2)-2{\hat{P}}_i^{\mathrm{dn}}(e_i^d-e_i^n)(f_i^d-f_i^n)}{{({(e_i^d-e_i^n)}^2+{(f_i^d-f_i^n)}^2)}^2}\\ {\mathrm{\μ}}_i^d=\frac{{\hat{P}}_i^{\mathrm{dn}}({(e_i^d-e_i^n)}^2-{(f_i^d-f_i^n)}^2)-2{\hat{Q}}_i^{\mathrm{dn}}(e_i^d-e_i^n)(f_i^d-f_i^n)}{{({(e_i^d-e_i^n)}^2+{(f_i^d-f_i^n)}^2)}^2}\end{array}\right.,d\∈B_P---\left(15\right)$

$\left\{\begin{array}{c}A=H^T\mathrm{WH}\\ B=H^T\mathrm{W\Δz}\end{array}\right.---\left(16\right)$

式中：B_p、d、

及

的意义同公式(5)；W、H和H^T的意义同公式(10)，

表示迭代值为x₀时的不平衡量。

(4)状态变量更新和收敛性判断

1)状态变量更新

第(3)步完成后，根据公式(17)计算状态变量的修正量Δx^(time)，然后更新状态变量，得到状态变量新值，即：x^(time+1)=x^(time)+Δx^(time)，time=time+1，计算公式为

Δx＝A^-1B-A^-1C^T((CA^-1C^T)^-1(CA^-1B-Δc))    (17)

式中，time为计算迭代次数；A^-1为矩阵A的逆矩阵；C^T为矩阵C的转置；Δc＝-c(x₀)。

2)收敛性判断

当状态变量的修正量Δx^(time)满足max(|Δx^(time)|)<ε，则结束迭代计算，输出结果；当max(|Δx^(time)|)≥ε且迭代次数time≥Tmax，则停止迭代，输出“不收敛!”。

当max(|Δx^(time)|)≥ε且迭代次数time<Tmax，转入进行如下步骤：若(δ^(time))²≥0.01，令(δ^(time+1))²＝0.001(δ^(time))²；若(δ^(time))²＜0.01，则令(δ^(time+1))²＝(δ^(time))²。并使迭代次数time增加1，返回第(3)步，进行重新迭代计算。

本发明采用上述技术方案后，主要有以下效果：

①与现有技术的权系数赋值方式比较，本发明采用权重随残差变化的指数型权函数，抗差性好，能有效抑制坏数据对估计结果的不良影响，且能够极大改善坏数据与量测误差共存时的收敛性和数值稳定性。

②在指数型权函数计算时引入阈值判据，采用试探法选取残差方差检测门槛值的最佳值。采用该小残差权重校正策略能改善状态估计方法的收敛性和收敛精度结果，具有强抗差性和良好的收敛性。

③本发明方法通过充分利用量测值，能独立计算节点有功和无功功率，能够避免功率量测要求成对出现且权系数相同的问题。

本发明方法可广泛有效应用于求解三相四线低压配电网的状态估计问题，对三相四线低压配电网的用电监测与管理具有良好的实用价值和在线应用前景。

附图说明

图1为本发明方法的程序流程框图；

图2为应用本发明方法进行三相四线低压配电网状态估计的两端点低压配电网等值电路；

图3为IEEE-13节点修正系统的网络图。

图中，i和j为端点编号；

(d∈B₁)分别是端点i、j中电气节点d对中性点的注入电流相量和对地的电压相量；B₁＝{a,b,c,n}，表示端点包含中性点的电气节点组合；

是端点i和j间节点d对中性点的支路电流相量；

和

(d∈B₁,d≠x)为端点i和j间线路支路的自阻抗和互阻抗；

和

为端点i和j中同一端点节点d和x之间的并联支路导纳；

和为端点i和j中节点d对地电容得到的对地并联导纳。

具体实施方式

下面结合具体实施方式，进一步说明本发明。

实施例

如图1-3所示，一种基于三相四线低压配电网状态估计方法的具体步骤如下：

(1)输入基础数据及初始化

1）输入基础数据

首先输入任一时间断面的低压配总表及各用户分表的瞬时抄表电力数据和网络结构及参数信息。即输入低压配总表在任一时间断面下的瞬时抄表总有功功率和总无功功率以及三相电压幅值和三相电流幅值，各用户分表在同一时间断面下的瞬时抄表各相有功功率和各相无功功率数据以及三相电压幅值，网络结构及参数信息为网络线路电阻、电抗和电纳参数、线路的额定电压、功率基准。

2）参数初始化

设置n阶单位矩阵R^-1（即n阶单位矩阵R^-1的对角线元素全为1，非对角线元素全为0），n为状态估计中实际的量测变量个数129；初始化最大迭代次数Tmax、收敛精度ε及残差方差的检测门槛值γ。根据1991年第6卷第3期《IEEE Transactions on Power Systems》中“Radialdistribution test feeders”一文关于IEEE-13节点系统的标准数据，输入网络结构、负荷数据及相关参数信息，并设置端点2的a相节点和端点13的c相节点为不良数据点，且正常量测均带有一定的量测误差，该误差符合正态分布。输入最大迭代次数为：Tmax=50；收敛精度为：ε=10^-4；残差方差的检测门槛值为：γ=5；令time=1。

(2)形成节点导纳矩阵和雅可比矩阵常数项部分

第(1)步完成后，分别计算该网络的节点导纳矩阵和雅可比矩阵常数项部分，计算公式为技术方案中的公式(1)-公式(3)。

根据IEEE-13节点修正系统的网络结构和参数信息，按照技术方案中的公式(1)，计算得到该网络的节点导纳矩阵Y。

$Y={\left[\begin{array}{ccccc}0.8221-2.9896i& -0.2096+0.4785i& ...& 0& 0\\ -0.2096+0.4785i& 0.8221-2.9896i& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& ...& 3.29-11.96i& -0.8396+1.9169i\\ 0& 0& ...& -0.8396+1.9169i& 3.29-11.96i\end{array}\right]}_{52\&times;52}$

根据IEEE-13节点修正系统的网络结构和节点导纳矩阵Y，按照技术方案中的公式(2)和公式(3)计算得到雅可比矩阵的常数项J0。

${J0=\left[\begin{array}{ccccccc}0.8221& -1.0384& 0.6192& ...& 0& 0& 0\\ -0.2096& 6.0002& 0.6192& ...& 0& 0& 0\\ -0.2096& -1.0384& -4.3561& ...& 0& 0& 0\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& 3.29& -0.8396& -0.8396\\ 0& 0& 0& ...& -0.8396& 3.29& -0.8396\\ 0& 0& 0& ...& -0.8396& -0.8396& 3.29\end{array}\right]}_{104\&times;104}$

(3)计算不平衡量、指数型权函数和雅可比矩阵

第(2)步完成后，基于量测方程、等式约束方程、指数型权函数的定义及电流注入型牛顿法计算相应的不平衡量、指数型权函数和雅可比矩阵，计算公式为技术方案中的公式(4)-公式(16)。

1)计算不平衡量

根据量测方程、等式约束方程、指数型权函数的定义及电流注入型牛顿法，按照技术方案中的公式(4)-公式(8)，计算得到量测方程与等式约束方程对应的不平衡量Δz和Δc。

以第1次迭代计算的结果举例，按技术方案中的公式(4)-公式(8)，计算得到不平衡量Δz和Δc为：

$\mathrm{\&Delta;z}={\left[\begin{array}{c}-2.1798\\ 0.4958\\ -0.6979\\ ...\\ 0.0017\\ 0.0021\\ 0.0015\end{array}\right]}_{129\&times;1};\mathrm{\&Delta;c}={\left[\begin{array}{c}6.8962E-10\\ 1.0709E-10\\ -6.5779E-10\\ 1.0914E-11\\ 0\\ 1.5029E-10\\ -1.8485E-10\\ 3.5652E-10\end{array}\right]}_{8\&times;1}$

2)计算指数型权函数

根据指数型权函数W的定义和小残差权重校正策略，按照技术方案中的公式(9)-公式(10)计算指数型权函数W和量测的标准化残差r_N。

以第1次迭代计算的结果举例，按技术方案中的公式(9)-公式(10)，计算得到指数型权函数W和量测的标准化残差r_N为：

$W={\left[\begin{array}{ccccccc}0.9999& 0& 0& ...& 0& 0& 0\\ 0& 1& 0& ...& 0& 0& 0\\ 0& 0& 1& ...& 0& 0& 0\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& 1& 0& 0\\ 0& 0& 0& ...& 0& 1& 0\\ 0& 0& 0& ...& 0& 0& 1\end{array}\right]}_{129\&times;129};r_N={\left[\begin{array}{c}0.9748\\ 0.2217\\ 0.3121\\ ...\\ 0.0008\\ 0.0010\\ 0.0007\end{array}\right]}_{129\&times;1}^T$

3)形成雅可比矩阵

根据量测方程、等式约束方程及电流注入型牛顿法，按照技术方案中的公式(11)-公式(15)计算雅可比矩阵H和C，并按照技术方案中的公式(16)计算A和B。

以第1次迭代计算的结果举例，按技术方案中的公式(11)-公式(15)，计算得到雅可比矩阵H和C为：

$H={\left[\begin{array}{ccccc}4.6179E+04& -2.6220E+03& ...& 0& 0\\ 1.3110E+03& -9.2357E+04& ...& 0& 0\\ 1.3110E+03& -2.6220E+03& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& ...& 0& 6.9125E+02\\ 0& 0& ...& 0& 4.0898E+03\\ 0& 0& ...& 3.4999E+03& -3.4999E+03\end{array}\right]}_{129\&times;99}$

$C={\left[\begin{array}{ccccc}...& 0.8396& -11.96& 1.9169& ...\\ ...& 0.8396& 1.9169& -11.96& ...\\ ...& 0.8396& 1.9169& 1.9169& ...\\ ...& -3.29& 1.9169& 1.9169& ...\\ ...& -1.9169& -3.29& 0.8396& ...\\ ...& -1.9169& 0.8396& -3.29& ...\\ ...& -1.9169& 0.8396& 0.8396& ...\\ ...& 11.96& 0.8396& 0.8396& ...\end{array}\right]}_{8\&times;99}$

以第1次迭代计算的结果举例，按技术方案中的公式(16)计算得到A和B为：

$A={\left[\begin{array}{ccccc}2.2416E+09& -4.1023E+08& ...& 0& 0\\ -4.1023E+08& 8.9634E+09& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& ...& 1.2250E+07& -1.2249E+07\\ 0& 0& ...& -1.2249E+07& 2.9454E+07\end{array}\right]}_{99\&times;99}$

$B={\left[\begin{array}{c}-1.0094E+05\\ -3.8244E+04\\ ...\\ 1.4398E+03\\ -1.4406E+03\end{array}\right]}_{99\&times;1}$

(4)状态变量更新和收敛性判断

1）状态变量更新

第(3)步完成后，按照技术方案中的公式(17)计算状态变量的修正量Δx^(time)，然后更新状态变量，得到状态变量新值，即：x^(time+1)=x^(time)+Δx^(time)，time=time+1。

以第1次迭代计算的结果举例，按技术方案中的公式(17)，计算得到状态变量的修正量Δx⁽¹⁾为：

$\mathrm{\&Delta;}{x}^{\left(1\right)}={\left[\begin{array}{ccccccc}...& -0.1423& 0.0106& 0.6434& -0.1243& -0.39& ...\end{array}\right]}_{1\&times;99}^T$

2）收敛性判断

若状态变量的修正量max(|Δx^(time)|)<ε，则停止迭代，输出结果；若max(|Δx^(time)|)>ε且迭代次数time>Tmax，则停止迭代，输出“不收敛!”。否则，转入进行如下步骤：若(δ^(time))²≥0.01，令(δ^(time+1))²＝0.001(δ^(time))²；若(δ^(time))²＜0.01，则令(δ^(time+1))²＝(δ^(time))²。并使迭代次数time加1，返回第(3)步，重新进行迭代。

以第1次迭代计算的结果举例，此时，ε=10^-4，Tmax=50，max(|Δx⁽¹⁾|)=9.6159>ε，time=1<Tmax。根据收敛性判断，进行如下步骤：按照技术方案中的δ²＝r^TPr/(m-n)，计算得到概率密度函数的标准差的初值为：(δ⁽¹⁾)²=9.4421E+03。此时，(δ⁽¹⁾)²>0.01，更新后，(δ⁽²⁾)²=0.001(δ⁽¹⁾)²=9.4421，time=time+1=2。返回第(3)步，重新进行迭代。

根据前面的步骤，迭代6次后满足收敛条件，此时max(|Δx⁽⁶⁾|)=1.5867E-06<ε，(δ⁽⁶⁾)²=0.0094，估计结果如下表1所示。

表1  三相负荷有功功率的量测估计结果

其中，粗体字标识的数据为量测坏数据。此时，可通过量测值与估计值之差(即为残差)由大至小排序及数值大小来辨识不良数据点，并通过本发明方法在一定程度上修正量测数据，减轻不良数据及量测误差对估计结果的影响。

结合IEEE-13节点修正系统算例及表1中的量测估计结果，通过残差排序和数值大小可以辨识出端点2的a相节点和端点13的c相节点为不良数据点，这与具体实施方式(1)输入基础数据中设定的量测坏数据点(或不良数据点)情况完全一致。

实验效果

以附图3所示的IEEE-13节点修正系统为对象，设计以下仿真算例验证本发明方法的有效性。

为简化及构造IEEE-13节点修正系统，对原IEEE-13节点系统进行如下基本修改：以大地和平衡端点的a相节点为零电位参考点，仅有电源端点(配变低压侧)的中性点直接接地，无端点接地阻抗，将支路型号均设置为501，且相间距设置相同，为该支路a-c间的相间距。

对IEEE-13节点修正系统，构造以下算例：1)系统存在不良数据，且正常量测均带有一定的量测误差。其中，真值通过三相四线配电网潮流计算获得。取原正常量测值的相反数为不良数据，量测误差构造参照1985年水利电力出版社出版的《电力系统状态估计》一书。设置端点2的a相节点和端点13的c相节点为不良数据点。除T端点5和平衡端点1外，其余端点均配置相对中性点的注入有功、无功功率及电压幅值量测。平衡端点1配置相对中性点的三相总注入有功、总注入无功、三相注入电流幅值及电压幅值量测。测点固定权重设置为单位阵，残差方差的检测门槛值γ取值为5，收敛精度为10^-4。2)基于算例1)，构造如下算例：2a)EFLS估计方法中无阈值判据(即不增设残差方差的检测门槛值)，量测数据和其它量的设置均同算例1)；2b)除EFLS估计方法中增设残差方差的检测门槛值γ为4外，其它条件同算例2a)；2c)除EFLS估计方法中残差方差的检测门槛值γ设置为3外，其它条件同算例2a)。

为验证本发明方法的特点，根据算例中的仿真条件分别计算量测值与真值之差、量测值与估计值之差以及估计值与真值之差用以分析。以量测值与估计值之差为例说明该计算方法，设第i个量测值为m_i，对应的估计值为n_i，则量测值与估计值之差定义为(|m_i-n_i|)²。单相和为差值按相分别相加，三相和为差值三相求和。

(1)坏数据抑制的抗差效果

算例1)情况下三相功率的量测估计结果如下表2所示，其中，粗体字标识数据为坏数据，表中数据均为有名值。在收敛精度为10^-4且确定量测误差及不良数据情况下，算例1)迭代6次收敛，收敛可靠，速度快。

表2  算例1)下三相功率的量测估计结果

由上述计算结果可见，本发明方法采用权重随残差变化的指数型权函数，抗差性好，可有效抑制坏数据对估计结果的不良影响，极大改善坏数据与量测误差共存时的收敛性和数值稳定性。同时，零注入等式约束的考虑可使得零注入功率端点(T端点)严格满足系统真实的潮流状态，与实际情况相符。

(2)小残差权重校正策略的收敛性改善效果

在收敛精度为10^-4且确定量测误差及不良数据情况下，算例2a)状态估计不收敛，算例2b)迭代8次收敛，算例2c)迭代10次收敛。算例2b)-2c)收敛情况下三相功率估计值与真值之差的计算结果如下表3所示，表中数据均为有名值。

表3  算例2b)-2c)下三相功率估计值与真值之差结果

由上述计算结果可见，阈值判据的引入及残差方差检测门槛值γ的适当取值可极大改善状态估计方法的收敛性，提高估计结果精度，即：小残差权重校正策略能够较好的改善状态估计方法的收敛性和收敛精度结果。本发明方法中，γ的取值采用试探法获得。

Claims (1)

1.一种三相四线低压配电网状态估计方法，利用计算机，通过程序，实现三相四线低压配电网的状态估计，其特征在于所述方法的具体步骤如下：

(1)输入基础数据及初始化

1)输入基础数据

首先输入任一时间断面的低压配总表及各用户分表的瞬时抄表电力数据和网络结构及参数信息，即输入低压配总表在任一时间断面下的瞬时抄表总有功功率和总无功功率以及三相电压幅值和三相电流幅值，各用户分表在同一时间断面下的瞬时抄表各相有功功率和各相无功功率数据以及三相电压幅值，网络结构及参数信息为网络线路电阻、电抗和电纳参数、线路的额定电压、功率基准；

2)参数初始化

设置n阶单位矩阵R^-1，n为状态估计中实际的量测变量个数；初始化最大迭代次数Tmax为40～60；收敛精度ε为10^-3～10^-5及残差方差的检测门槛值γ=5，并设置迭代次数time=1；

(2)形成节点导纳矩阵和雅可比矩阵常数项部分

第(1)步完成后，计算节点导纳矩阵，计算公式为：

式中：B₁＝{a,b,c,n}，表示端点包含中性点n的电气节点组合，其中a、b、c分别表示abc三相电气节点；x、d、t表示电气节点组合B₁中的任一节点；φ_i为与端点i直接相联但不包括端点i的端点集合；j为端点集合φ_i中的任一端点；(x∈B₁&x≠d)表示端点i中节点d和x之间的并联支路导纳；

表示端点i中节点d与端点j中节点t直接相联的支路导纳；表示端点i中节点d对地电容得到的对地并联导纳；自导纳

表示与端点i节点d直接相联的所有支路导纳和对地电容导纳

之和，互导纳

(t≠d)表示端点i中节点d和t之间支路导纳

的相反数与端点i且与i直接相联端点节点d和t之间所有支路导纳之和，互导纳

(k≠i)则表示端点i和k中节点d和t之间支路导纳

的相反数；

得到节点导纳矩阵后，计算雅可比矩阵常数项子矩阵，计算公式为：

$\left\{\begin{array}{c}{H}_{\mathrm{ij}}(1:\mathrm{8,1}:8)={\left[\begin{array}{cc}{G}_{\mathrm{ij}}^{\mathrm{dt}}& -B_{\mathrm{ij}}^{\mathrm{dt}}\\ B_{\mathrm{ij}}^{\mathrm{dt}}& G_{\mathrm{ij}}^{\mathrm{dt}}\end{array}\right]}_{8\&times;8},d\&Element;B_1,t\&Element;B_1\\ H_{\mathrm{ij}}(9:\mathrm{11,1}:8)={\left[0\right]}_{3\&times;8}\end{array}\right.---\left(2\right)$

C_ij(1:8,1:8)＝H_ij(1:8,1:8)    (3)

式中：i和j为端点编号；H_ij和C_ij分别为雅可比矩阵H和C的常数项子矩阵；

和

分别为节点导纳矩阵元素

的实部和虚部；[0]表示全零矩阵，下标3×8表示该矩阵为3行8列；d、t及B₁的意义同公式(1)；

(3)计算不平衡量和指数型权函数及雅可比矩阵

第(2)步完成后，计算三相四线低压配电网中量测量的不平衡量和指数型权函数及雅可比矩阵，具体步骤如下：

1)计算不平衡量

基于公式(4)计算三相四线低压配电网总表端点s的电流不平衡量和有功功率不平衡量及无功功率不平衡量，计算公式为：

式中：B_P＝{a,b,c}，表示端点不包含中性点n的电气节点组合，其中a、b、c分别表示abc三相电气节点；d表示电气节点组合B_P中的任一节点；B₁、t的意义同公式(1)；s表示三相四线低压配电网总表的端点；

为包括总表端点s且与其直接相联端点的集合；k为端点集合

中的任一端点；Δz_s为基于公式(4)计算端点s的电流不平衡量和功率不平衡量；(d∈B_P)、

和

分别为低压配电网总表s相对中性点n的三相电流幅值、三相总注入有功和三相总注入无功的瞬时抄表量测量；

和

分别为端点k中节点t的电压

的实部和虚部；和分别为端点s中节点d的电压

的实部和虚部；

和

分别为端点s的中性点n的电压的实部和虚部；和

分别为节点导纳矩阵元素

的实部和虚部；

基于公式(5)计算三相四线低压配电网用户分表端点i中a、b、c三相节点注入有功功率及注入无功功率的不平衡量，计算公式为：

式中：

和

分别为低压配电网用户分表i中节点d相对中性点n的注入有功和注入无功的瞬时抄表量测量；i表示三相四线低压配电网中任意某个用户分表端点；

和

(d∈B_p)分别为端点i中节点d的电压

的实部和虚部；

和

分别为端点i的中性点n电压

的实部和虚部；Δz_i为基于公式(5)计算端点i的节点注入功率的不平衡量；和分别为节点导纳矩阵元素

的实部和虚部；

为包括端点i且与其直接相联端点的集合；k为端点集合

中的任一端点；

B_P、d、B₁及t的意义同公式(4)；

基于公式(6)计算三相四线低压配电网用户分表端点i的中性点注入有功功率及注入无功功率的不平衡量，基于公式(7)计算三相四线低压配电网用户分表端点i中a、b、c三相节点电压的不平衡量，计算公式分别为：

$\mathrm{\&Delta;}{z}_i={\left({\hat{U}}_i^{\mathrm{dn}}\right)}^2-({(e_i^d-e_i^n)}^2+{(f_i^d-f_i^n)}^2),d\&Element;B_P---\left(7\right)$

式中：

为低压配电网分表i相对中性点n的电压幅值的瞬时抄表量测量；

和分别为节点导纳矩阵元素

的实部和虚部；Δz_i为用户分表端点i计算得到的功率和电压不平衡量； k、B_P、d、B₁及t的意义同公式(5)；

基于公式(8)计算零注入电流端点m中a、b、c三相及中性点注入电流的不平衡量Δc_m，其计算公式为：

式中：和分别为节点导纳矩阵元素的实部和虚部；

为包括T端点m且与其直接相联端点的集合；k为端点集合

中的任一端点；

的意义同公式(4)；B₁、t、d的意义同公式(1)；

2)计算指数型权函数

基于公式(9)-公式(10)计算指数型权函数W的对角阵元素，计算公式为：

$w_i^{*}=R_i^{-1}{e}^{{-r}_{\mathrm{Ni}}^2/{2\mathrm{\&delta;}}^2}---\left(9\right)$

式中：w_i ^*为量测i的权函数；

为量测i的固定权重；δ为概率密度函数的标准差，其初值由δ²＝r^TPr/(m-n)计算得到，r为残差向量，r^T为残差向量的转置，P为先验矩阵

m为量测量的个数，n为状态变量的个数；r_Ni为量测i的标准化残差，计算公式为

$r_{\mathrm{Ni}}=\left\{\begin{array}{cc}\left|r_i\right|/\sqrt{\left|{\left(\mathrm{\&Delta;R}\right)}_i\right|},& \left|{\left(\mathrm{\&Delta;R}\right)}_i\right|>\mathrm{\&gamma;}\\ \left|r_i\right|/\sqrt{\mathrm{\&gamma;}},& \left|{\left(\mathrm{\&Delta;R}\right)}_i\right|\&le;\mathrm{\&gamma;}\end{array}\right.---\left(10\right)$

式中：残差方差的检测门槛值γ通过试探法或经验值选取；

为量测i对应的残差；R为测量误差方差的对角阵，对角线元素为测量误差的方差；Δ＝Ι-H(H^TWH)^-1H^TW；Ι为n阶单位阵；W为指数型权函数对角阵，其对角元素等于权函数，即W_ii＝w_i ^*；

为雅可比矩阵，H^T为雅克比矩阵的转置；(H^TWH)^-1表示矩阵H^TWH相乘后的逆矩阵；

3)形成雅可比矩阵

基于公式(11)-公式(15)形成雅可比矩阵H和C，根据公式(16)计算A和B；

对于平衡节点，基于公式(11)和(12)形成其雅克比矩阵元素，计算公式为：

式中：雅可比子矩阵H_ss＝[H_ss1(1:3,1:3);H_ss2(1:3,1:3);H_ss3(1,1:3);H_ss4(1,1:3)]；

和分别为节点导纳矩阵元素

和的实部，

和

分别为节点导纳矩阵元素

和

的虚部；

和

分别为端点f中节点t电压

的实部和虚部，f为端点集合

中的任一端点；和

分别为端点s中的节点p电压

的实部和虚部；

B_P、B₁、d、t、s的意义同公式(4)；p表示电气节点组合B_P中的任一节点；

为平衡端点s中节点p电压虚部对平衡端点s中节点p电压实部的偏导数，且 $\frac{{\&PartialD;f}_s^a}{{\&PartialD;e}_s^a}=0,\frac{{\&PartialD;f}_s^b}{{\&PartialD;e}_s^b}=\sqrt{3},\frac{{\&PartialD;f}_s^c}{{\&PartialD;e}_s^c}=-\sqrt{3};$

$\left\{\begin{array}{c}{H}_{\mathrm{sj}1}(1:\mathrm{3,1}:8)={\left[\begin{array}{c}2G_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t-B_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t)+2B_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t+B_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t)\\ -2B_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t-B_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t)+2G_{\mathrm{sj}}^{\mathrm{dt}}(G_{\mathrm{sj}}^{\mathrm{dt}}{f}_j^t+B_{\mathrm{sj}}^{\mathrm{dt}}{e}_j^t)\end{array}\right]}_{3\&times;8}^T\\ H_{\mathrm{sj}2}(1:\mathrm{3,1}:8)={\left[0\right]}_{3\&times;8}\\ H_{\mathrm{sj}3}(\mathrm{1,1}:8)={\left[\begin{array}{c}\underset{d\&Element;B_P}{\mathrm{\&Sigma;}}\{G_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)+B_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\\ \underset{d\&Element;B_P}{\mathrm{\&Sigma;}}\{-B_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)+G_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\end{array}\right]}_{1\&times;8}^T\\ H_{\mathrm{sj}4}(\mathrm{1,1}:8)={\left[\begin{array}{c}\underset{d\&Element;B_P}{\mathrm{\&Sigma;}}\{-B_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)+G_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\\ \underset{d\&Element;B_P}{\mathrm{\&Sigma;}}\{-G_{\mathrm{sj}}^{\mathrm{dt}}(e_s^d-e_s^n)-B_{\mathrm{sj}}^{\mathrm{dt}}(f_s^d-f_s^n)\}\end{array}\right]}_{1\&times;8}^T\end{array}\right.---\left(12\right)$

式中：雅可比子矩阵H_sj＝[H_sj1(1:3,1:8);H_sj2(1:3,1:8);H_sj3(1,1:8);H_sj4(1,1:8)]；s为平衡端点；j为非平衡端点；和

分别为节点导纳矩阵元素

的实部和虚部；

和

分别为端点j中节点t电压

的实部和虚部；B_P、

d及t的意义同公式(4)；

对于非平衡节点，基于公式(13)～(15)形成其雅克比矩阵元素，计算公式为：

$\left\{\begin{array}{c}{H}_{\mathrm{is}}(1:\mathrm{8,1}:3)={[G_{\mathrm{is}}^{\mathrm{dp}}-B_{\mathrm{is}}^{\mathrm{dp}}\frac{{\&PartialD;f}_s^p}{{\&PartialD;e}_s^p};B_{\mathrm{is}}^{\mathrm{dp}}+G_{\mathrm{is}}^{\mathrm{dp}}\frac{{\&PartialD;f}_s^p}{{\&PartialD;e}_s^p}]}_{8\&times;3},d\&Element;B_1,p\&Element;B_P\\ H_{\mathrm{is}}(9:11,1:3)={\left[0\right]}_{3\&times;3}\end{array}\right.---\left(13\right)$

式中：

和

分别为节点导纳矩阵元素的实部和虚部；

B_P、B₁、p的意义同公式(11)；d表示电气节点组合B₁中的任一节点；

$\left\{\begin{array}{c}{\mathrm{\&Delta;H}}_{\mathrm{ii}}(1:\mathrm{3,1}:3)=-\mathrm{\&Delta;}{H}_{\mathrm{ii}}(5:\mathrm{7,5}:7)={\left[{\mathrm{\&mu;}}_i^d\right]}_{3\&times;3}\\ \mathrm{\&Delta;}{H}_{\mathrm{ii}}(5:\mathrm{7,8})=-\mathrm{\&Delta;}{H}_{\mathrm{ii}}(1:\mathrm{3,4})={\left[{\mathrm{\&mu;}}_i^d\right]}_{3\&times;1}\\ \mathrm{\&Delta;}{H}_{\mathrm{ii}}(\mathrm{8,5}:7)=-\mathrm{\&Delta;}{H}_{\mathrm{ii}}(\mathrm{4,1}:3)={\left[{\mathrm{\&mu;}}_i^d\right]}_{1\&times;3}\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}\left(\mathrm{4,4}\right)=-{\mathrm{\&Delta;H}}_{\mathrm{ii}}\left(\mathrm{8,8}\right)=\underset{d\&Element;B_P}{\mathrm{\&Sigma;}}{\mathrm{\&mu;}}_i^d\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}(1:3,5:7)=\mathrm{\&Delta;}{H}_{\mathrm{ii}}(5:7,1:3)={[-{\mathrm{\&eta;}}_i^d]}_{3\&times;3}\\ \mathrm{\&Delta;}{H}_{\mathrm{ii}}(1:3,8)=\mathrm{\&Delta;}{H}_{\mathrm{ii}}(5:7,4)={\left[{\mathrm{\&eta;}}_i^d\right]}_{3\&times;1}\\ \mathrm{\&Delta;}{H}_{\mathrm{ii}}(4,5:7)=\mathrm{\&Delta;}{H}_{\mathrm{ii}}(8,1:3)={\left[{\mathrm{\&eta;}}_i^d\right]}_{1\&times;3}\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}(4,8)={\mathrm{\&Delta;H}}_{\mathrm{ii}}(8,4)=-\underset{d\&Element;B_P}{\mathrm{\&Sigma;}}{\mathrm{\&eta;}}_i^d\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}(9:11,1:3)={\left[2(e_i^d-e_i^n)\right]}_{3\&times;3}\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}(9:\mathrm{11,4})={[-2(e_i^d-e_i^n)]}_{3\&times;1}\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}(9:\mathrm{11,5}:7)={\left[2(f_i^d-f_i^n)\right]}_{3\&times;3}\\ {\mathrm{\&Delta;H}}_{\mathrm{ii}}(9:\mathrm{11,8})={[-2(f_i^d-f_i^n)]}_{3\&times;1}\end{array}\right.---\left(14\right)$

式中：ΔH_ii为雅可比子矩阵H_ii的修正项子阵，i为非平衡端点；B_P、d、及

的意义同公式(5)；

和

均为修正项子阵ΔH_ii的矩阵元素，其计算公式分别为：

$\left\{\begin{array}{c}{\mathrm{\&eta;}}_i^d=\frac{{\hat{Q}}_i^{\mathrm{dn}}({(e_i^d-e_i^n)}^2-{(f_i^d-f_i^n)}^2)-2{\hat{P}}_i^{\mathrm{dn}}(e_i^d-e_i^n)(f_i^d-f_i^n)}{{({(e_i^d-e_i^n)}^2+{(f_i^d-f_i^n)}^2)}^2}\\ {\mathrm{\&mu;}}_i^d=\frac{{\hat{P}}_i^{\mathrm{dn}}({(e_i^d-e_i^n)}^2-{(f_i^d-f_i^n)}^2)-2{\hat{Q}}_i^{\mathrm{dn}}(e_i^d-e_i^n)(f_i^d-f_i^n)}{{({(e_i^d-e_i^n)}^2+{(f_i^d-f_i^n)}^2)}^2}\end{array}\right.,d\&Element;B_P---\left(15\right)$

$\left\{\begin{array}{c}A=H^T\mathrm{WH}\\ B=H^T\mathrm{W\&Delta;z}\end{array}\right.---\left(16\right)$

式中：B_p、d、

及

的意义同公式(5)；W、H和H^T的意义同公式(10)，

表示迭代值为x₀时的不平衡量；

(4)状态变量更新和收敛性判断

1)状态变量更新

第(3)步完成后，根据公式(17)计算状态变量的修正量Δx^(time)，然后更新状态变量，得到状态变量新值，即：x^(time+1)=x^(time)+Δx^(time)，time=time+1，计算公式为

Δx＝A^-1B-A^-1C^T((CA^-1C^T)^-1(CA^-1B-Δc))    (17)

式中，time为计算迭代次数；A^-1为矩阵A的逆矩阵；C^T为矩阵C的转置；Δc＝-c(x₀)；

2)收敛性判断

当状态变量的修正量Δx^(time)满足max(|Δx^(time)|)<ε，则结束迭代计算，输出结果；当max(|Δx^(time)|)≥ε且迭代次数time≥Tmax，则停止迭代，输出“不收敛!”；

当max(|Δx^(time)|)≥ε且迭代次数time<Tmax，转入进行如下步骤：若(δ^(time))²≥0.01，令(δ^(tim+e1))²＝0.001(δ^(tim)e)²；若(δ^(time))²＜0.01，则令(δ^(time+1))²＝(δ^(time))²；并使迭代次数time增加1，返回第(3)步，进行重新迭代计算。

Images