[发明专利]一种三相四线低压配电网的状态估计方法

权利要求书

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之间支路导纳的相反数；

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

Hij(1:8,1:8)=Gijdt-BijdtBijdtGijdt8×8,d∈B1,t∈B1Hij(9:11,1:8)=[0]3×8---(2)]]>

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三相节点电压的不平衡量，计算公式分别为：

Δzi=(U^idn)2-((eid-ein)2+(fid-fin)2),d∈BP---(7)]]>

式中：为低压配电网分表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的对角阵元素，计算公式为：

wi*=Ri-1e-rNi2/2δ2---(9)]]>

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

rNi=|ri|/|(ΔR)i|,|(ΔR)i|>γ|ri|/γ,|(ΔR)i|≤γ---(10)]]>

式中：残差方差的检测门槛值γ通过试探法或经验值选取；为量测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电压实部的偏导数，且∂fsa∂esa=0,∂fsb∂esb=3,∂fsc∂esc=-3;]]>

Hsj1(1:3,1:8)=2Gsjdt(Gsjdtejt-Bsjdtfjt)+2Bsjdt(Gsjdtfjt+Bsjdtejt)-2Bsjdt(Gsjdtejt-Bsjdtfjt)+2Gsjdt(Gsjdtfjt+Bsjdtejt)3×8THsj2(1:3,1:8)=[0]3×8Hsj3(1,1:8)=Σd∈BP{Gsjdt(esd-esn)+Bsjdt(fsd-fsn)}Σd∈BP{-Bsjdt(esd-esn)+Gsjdt(fsd-fsn)}1×8THsj4(1,1:8)=Σd∈BP{-Bsjdt(esd-esn)+Gsjdt(fsd-fsn)}Σd∈BP{-Gsjdt(esd-esn)-Bsjdt(fsd-fsn)}1×8T---(12)]]>

式中：雅可比子矩阵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)形成其雅克比矩阵元素，计算公式为：

His(1:8,1:3)=[Gisdp-Bisdp∂fsp∂esp;Bisdp+Gisdp∂fsp∂esp]8×3,d∈B1,p∈BPHis(9:11,1:3)=[0]3×3---(13)]]>

式中：和分别为节点导纳矩阵元素的实部和虚部；B_P、B₁、p的意义同公式(11)；d表示电气节点组合B₁中的任一节点；

ΔHii(1:3,1:3)=-ΔHii(5:7,5:7)=[μid]3×3ΔHii(5:7,8)=-ΔHii(1:3,4)=[μid]3×1ΔHii(8,5:7)=-ΔHii(4,1:3)=[μid]1×3ΔHii(4,4)=-ΔHii(8,8)=Σd∈BPμidΔHii(1:3,5:7)=ΔHii(5:7,1:3)=[-ηid]3×3ΔHii(1:3,8)=ΔHii(5:7,4)=[ηid]3×1ΔHii(4,5:7)=ΔHii(8,1:3)=[ηid]1×3ΔHii(4,8)=ΔHii(8,4)=-Σd∈BPηidΔHii(9:11,1:3)=[2(eid-ein)]3×3ΔHii(9:11,4)=[-2(eid-ein)]3×1ΔHii(9:11,5:7)=[2(fid-fin)]3×3ΔHii(9:11,8)=[-2(fid-fin)]3×1---(14)]]>

式中：ΔH_ii为雅可比子矩阵H_ii的修正项子阵，i为非平衡端点；B_P、d、及的意义同公式(5)；和均为修正项子阵ΔH_ii的矩阵元素，其计算公式分别为：

ηid=Q^idn((eid-ein)2-(fid-fin)2)-2P^idn(eid-ein)(fid-fin)((eid-ein)2+(fid-fin)2)2μid=P^idn((eid-ein)2-(fid-fin)2)-2Q^idn(eid-ein)(fid-fin)((eid-ein)2+(fid-fin)2)2,d∈BP---(15)]]>

A=HTWHB=HTWΔz---(16)]]>

式中：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)步，进行重新迭代计算。