3节电力系统状态估计(WLS算法)

电力系统状态估计研究综述摘要：电力系统状态估计是当代电力系统能量管理系统(EMS)的重要组成部分。

本文介绍了电力系统状态估计的概念、数学模型，阐述了状态估计的必要性及其作用，系统介绍了状态估计的研究现状，最后对状态估计的研究方向进行了展望。

关键词：电力系统；状态估计；能量管理系统0引言状态估计是当代电力系统能量管理系统(EMS)的重要组成部分，尤其在电力市场环境中发挥更重要的作用。

它是将可用的冗余信息(直接量测值及其他信息) 转变为电力系统当前状态估计值的实时计算机程序和算法。

准确的状态估计结果是进行后续工作(如安全分析、调度员潮流和最优潮流等)必不可少的基础。

随着电力市场的发展，状态估计的作用更显重要⑴o状态估计的理论研究促进了工程应用，而状态估计软件的工程应用也推动了状态估计理论的研究和发展。

迄今为止，这两方面都取得了大量成果。

然而，状态估计领域仍有不少问题未得到妥善解决，随着电力系统规模的不断扩大，电力工业管理体制向市场化迈进，对状态估计有了新要求，各种新技术和新理论不断涌现，为解决状态估计的某些问题提供了可能。

本文就电力系统状态估计的研究现状和进一步的研究方向进行了综合阐述。

1电力系统状态估计的概念1.1电力系统状态估计的基本定义状态估计也被称为滤波，它是利用实时量测系统的冗余度来提高数据精度，自动排除随机干扰所引起的错误信息，估计或预报系统的运行状态(或轨迹)o 状态估计作为近代计算机实时数据处理的手段，首先应用于宇宙飞船、卫星、导弹、潜艇和飞机的追踪、导航和控制中。

它主要使用了六十年代初期由卡尔曼、布西等人提出的一种递推式数字滤波方法，该方法既节约内存，又大大降低了每次估计的计算量[2,4]o电力系统状态估计的研究也是由卡尔曼滤波开始。

但根据电力系统的特点，即状态估计主要处理对象是某一时间断面上的高维空间(网络)问题，而且对量测误差的统计知识又不够清楚，因此便于采用基于统计学的估计方法如最小方差估计、极大验后估计、极大似然估计等方法，目前很多电力系统实际采用的状态估计算法是最小二乘法。

第五讲电力系统状态估计概述电力系统状态估计指的是通过对电力系统的监测和测量数据进行处理，推算出电力系统相关参数的过程。

通俗的说，就是在电力系统的运行过程中，通过监测数据估计电力系统的状态，以便于运行员做出更好的决策。

电力系统状态估计的意义电力系统状态估计是电力系统自动化的重要组成部分。

在电力系统运行过程中，状态估计系统可以帮助运行员迅速掌握系统状态，及时调整电力系统的运行方式，保证电力系统的安全运行。

同时，状态估计系统还能够优化系统的经济性，提高电力系统的可靠性。

电力系统状态估计的原理电力系统状态估计是基于电力系统监测数据的处理和分析而实现的。

电力系统监测数据主要包括电压、电流、功率等参数。

通过对这些参数的监测和测量，可以获取电力系统的当前状态。

状态估计系统主要是通过对监测数据的处理和分析，以及对电力系统的模型建立和分析来推算电力系统的状态。

电力系统状态估计的原理和方法很多，但基本流程是相似的。

首先需要对电力系统的模型进行建立和分析，然后根据监测数据和运行状态信息，结合电力系统模型，对电力系统的状态进行估计。

最后根据状态估计结果，进行决策和调整。

电力系统状态估计的关键技术为了实现电力系统状态估计，需要涉及到诸多技术。

其中，关键技术包括：变电站数据采集系统变电站是电力系统中起到极为重要作用的环节，所以变电站的监测数据是状态估计的重要来源之一。

因此，变电站数据采集系统的高可靠性和高稳定性是保证状态估计准确性的关键。

现代电力系统常用的数据采集系统包括智能终端设备、数字遥测与遥控设备等。

电力系统模型状态估计需要基于电力系统模型来进行推算。

电力系统模型就是对电力系统运行模式进行建模和仿真得到的电力系统模拟实验环境。

常见的电力系统模型主要有潮流计算模型、电容器模型和风电模型等。

数据预处理电力系统的监测数据通常包含了大量的噪声，因此需要对数据进行预处理。

常用的数据预处理方法包括滤波、降噪、数据插补等等。

非线性方程组求解电力系统状态估计需要根据监测数据在电力系统模型的基础上求解非线性方程组，所以求解非线性方程组是状态估计的关键技术。

3節電系統狀態估計報告【任務說明】:闭合的开关:打开的开关:打开的刀闸:线路:负荷G:发电机:母线:连接线(没有阻抗) Unit2Unit13節點系統主接線圖任務：1、采用最小二乘狀態估計算法,所有量測の權重都取1.0,編寫狀態估計程序(C/Matlab)。

2、按量測類型，列出量測方程（每一類寫出一個方程）3、畫出程序流程4、提交源程序，程序中每個函數の作用5、提交計算の輸出結果(屏幕拷貝)系統參數：功率基值：100MW電壓基值：230 kV線路阻抗參數（標麼值）：線路量測(流出母線為正)：母線電壓量測：負荷量測(流出母線為正)：發電量測(流入母線為正)：注：量測存在誤差【數據預處理】首先根據基值將已知の量測值均轉換為標么值，並將功率值轉換為流入量，得到如下數據：線路導納參數（標麼值）：線路注入功率量測(標么值)：負荷點注入功率量測(標么值)：發電機節點注入量測(流入母線為正)：母線電壓量測（標么值）：【量測方程】選擇節點1の電壓相角為參考，為0度，以vi表示誤差值。

1）節點1電壓量測方程：Vi=Vi+v1即1.0087=V1+v12）1-3支路1號節點處注入有功功率功率：P ij=V i2g ij-V i V j(g ij cos+b ij sin)+v20.613=V12g13-V1V3(g13cos+b13sin)+v2即0.613=-1.6171V12-V1V3(-1.6171cos +13.698sin)+v2 3）1號節點注入功率：P i=V i2G ii +G ij cos+B ij sin+v3P1=V12G11+G1j cos+B1j sin+v3即-1.11=3.5613V12+V1V2(-1.9442cos -10.5107sin)+V1V3(-1.6171 cos -13.698 sin)+v3【流程圖】其中iterations 為迭代次數，可見本例の迭代次數為4，收斂較快，狀態估計得到の節點1、2、3電壓分別為：234.0144444444444444444444444444444444444444444444，，A X=b,得X+Xmax|X|<【程序說明】1、計算h矩陣の函數cal_hfunction h=cal_h(V,th0,B,G) %其中，V為節點電壓估計值，th0為節點電壓相角估計%值，B為節點電導矩陣，G為節點電納矩陣b=-B; %線路電導矩陣g=-G; %線路電納矩陣P=zeros(3,1); %初始化，節點注入功率Q=zeros(3,1);PP=zeros(3,3); %線路注入功率QQ=PP;th=[0;th0]; %節點1の電壓相角為0for i=1:3P_P=0;Q_Q=0;for j=1:3if(j~=i)P_P=P_P+V(i)*V(j)*(G(i,j)*cos(th(i)-th(j))+B(i,j)*sin(th(i)-th(j)));Q_Q=Q_Q+V(i)*V(j)*(G(i,j)*sin(th(i)-th(j))-B(i,j)*cos(th(i)-th(j)));PP(i,j)=(V(i)^2)*g(i,j)-V(i)*V(j)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));QQ(i,j)=-(V(i)^2)*b(i,j)-V(i)*V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j)));endendP(i)=(V(i)^2)*G(i,i)+P_P;Q(i)=-(V(i)^2)*B(i,i)+Q_Q;endVV=[V(1);V(2);V(3)];h=[P;Q;PP(1,2);PP(2,3);PP(3,1);QQ(1,2);QQ(2,3);QQ(3,1);PP(1,3);PP(2,1);P P(3,2);QQ(1,3);QQ(2,1);QQ(3,2);VV];2、計算H矩陣の函數cal_HHfunction H=cal_HH(V,th0,G,B,P,Q) %其中，P，Q為根據電壓估計值計算得到の節點%注入電壓b=-B;g=-G;PV=zeros(3,3); %節點注入功率對電壓幅值の偏導數QV=zeros(3,3);Pth=zeros(3,3); %節點注入功率對電壓相角の偏導數Qth=zeros(3,3);PPV=zeros(3,3); %P ij對V jの偏導數QQV=zeros(3,3); %Q ij對V jの偏導數PPth=zeros(3,3); %P ij對th jの偏導數QQth=zeros(3,3); %Q ij對th jの偏導數PPV1=zeros(3,3); %P ij對V iの偏導數QQV1=zeros(3,3); %Q ij對V iの偏導數PPth1=zeros(3,3); %P ij對th iの偏導數QQth1=zeros(3,3); %Q ij對th iの偏導數VV=eye(3);Vth=zeros(3,2);th=[0;th0];for i=1:3for j=1:3if (i~=j)PV(i,j)=V(i)*(G(i,j)*cos(th(i)-th(j))+B(i,j)*sin(th(i)-th(j)));QV(i,j)=V(i)*(G(i,j)*sin(th(i)-th(j))-B(i,j)*cos(th(i)-th(j)));PPV(i,j)=-V(i)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));QQV(i,j)=-V(i)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j)));PPV1(i,j)=2*V(i)*g(i,j)-V(j)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));QQV1(i,j)=-2*V(i)*b(i,j)-V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j)));if (j~=1)Pth(i,j)=V(i)*V(j)*(G(i,j)*sin(th(i)-th(j))-B(i,j)*cos(th(i)-th(j)));Qth(i,j)=-V(i)*V(j)*(G(i,j)*cos(th(i)-th(j))+B(i,j)*sin(th(i)-th(j)));PPth(i,j)=-V(i)*V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j)));QQth(i,j)=-V(i)*V(j)*(-g(i,j)*cos(th(i)-th(j))-b(i,j)*sin(th(i)-th(j)));endif(i~=1)PPth1(i,j)=V(i)*V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j))); QQth1(i,j)=-V(i)*V(j)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));endelsePV(i,j)=(G(i,i)*(V(i)^2)+P(i))/V(i);QV(i,j)=(Q(i)-(V(i)^2)*B(i,i))/V(i);if (j~=1)Pth(i,j)=-B(i,i)*(V(i)^2)-Q(i);Qth(i,j)=P(i)-(V(i)^2)*G(i,i);endendendendH=[[PV,Pth(:,2:3)];[QV,Qth(:,2:3)];...[PPV1(1,2),PPV(1,2),0,PPth(1,2),0;...0,PPV1(2,3),PPV(2,3),PPth1(2,3),PPth(2,3);...PPV(3,1),0,PPV1(3,1),0,PPth1(3,1)];...[QQV1(1,2),QQV(1,2),0,QQth(1,2),0;...0,QQV1(2,3),QQV(2,3),QQth1(2,3),QQth(2,3);...QQV(3,1),0,QQV1(3,1),0,QQth1(3,1)];...[PPV1(1,3),0,PPV(1,3),0,PPth(1,3);...PPV(2,1),PPV1(2,1),0,PPth1(2,1),0;...0,PPV(3,2),PPV1(3,2),PPth(3,2),PPth1(3,2)];...[QQV1(1,3),0,QQV(1,3),0,QQth(1,3);...QQV(2,1),QQV1(2,1),0,QQth1(2,1),0;...0,QQV(3,2),QQV1(3,2),QQth(3,2),QQth1(3,2)];...[VV,Vth]];3、主程序calculate_all.m文件format longG=[3.5613,-1.9442,-1.6171;...-1.9442,3.0993,-1.1551;...-1.6171,-1.1551,2.7722]; %B為節點電導矩陣B=[-24.2087,10.5107,13.698;...10.5107,-20.295,9.7843;...13.698,9.7843,-23.4832]; %G為節點電納矩陣P=[-1.11;0.88;0.23]; %節點注入功率量測值Q=[-0.135;-0.0424;0.24];PP=[0.613;-0.24;-0.459]; %線路1-2,2-3,3-1注入功率在首端の量測值QQ=[-0.012;0.066;-0.165];PP1=[0.467;-0.6;0.24]; %線路1-3,2-1,3-2注入功率在首端の量測值QQ1=[0.148;-0.024;-0.072];V=[1.0087;1.0198;1.0281]; %節點電壓幅值量測值R=diag(ones(21,1)); %權重都取為1Z=[P;Q;PP;QQ;PP1;QQ1;V]; %量測值矩陣V0=[1;1;1]; %初值th0=[0;0];delta=100;iterations=0; %迭代次數while delta>0.000001iterations=iterations+1;h=cal_h(V0,th0,B,G); %計算h矩陣H=cal_HH(V0,th0,G,B,h(1:3,1),h(4:6,1)); %計算H矩陣A=H'*inv(R)*H;b=H'*inv(R)*(Z-h);d=A\b; %求解修正值delta=max(abs(d));V0=V0+d(1:3,1); %修正估計值th0=th0+d(4:5,1);enditerationsV0=V0*230; %轉換為有名值th0=th0*180/pi; %轉換為度for i=1:3j=num2str(i);v=num2str(V0(i));show1=strcat('The voltage magnitude of node ',j,' is', v,' kV');disp(show1);endfor i=1:2j=num2str(i+1);th=num2str(th0(i));show1=strcat('The phase angle of node ',j,' is ',th,' degrees');disp(show1);end。

一种基于WLS的电网分层状态估计方法研究朱彬【期刊名称】《《中国设备工程》》【年(卷),期】2019(000)022【总页数】2页(P16-17)【关键词】状态估计; 加权; 最小二乘法; WLS【作者】朱彬【作者单位】上海恒合电力工程有限公司上海 200000【正文语种】中文【中图分类】TM7441 引言电力系统遥信遥测可能存在不良数据。

例如，某一个节点有3条支路，流向该节点的电流实际值为i1、i2、i3，测量值为，根据基尔霍夫电流定律，有i1+i2+i3=0；但是由于测量存在误差，一般--≠0。

这样的测量数据不能直接用于计算，必须要经过适当的处理。

这个处理的过程就称为状态估计。

电力系统状态估计即在测量误差时，运用实时量测系统的冗余度提高测量精度，排除随机干项导致的数据错误，最终通过数学计算估计出系统的运行状态，得到可靠的并且为数最少的系统状态变量值。

2 电网分层网络状态估计思想如图1所示，是一段配电网的简单模型。

图1 小型配电网模型对这个电网的模型进行网络拓扑可以得到：图2 配电网拓扑结构由于配电网呈辐射形、树枝状，单源而且没有回路。

考虑到电网的分层，可以采用从叶子节点一层一层地往根节点推进的方式进行状态估计。

如图2所示，共有8个节点，每个节点都有测量仪器的测量数据。

以各节点注入的有功功率为例，在不考虑损耗的情况下，各点的有功功率应该满足：但是测量得到往的数据往往不能满足这样的条件。

本文采用的方法是从网络拓扑图的最底层开始，逐层向上使用加权最小二乘法进行优化求解，并将得到的新数据作为往上一层进行加权最小二乘法求解时的伪量测量。

由于伪量测量是经过一次初次处理的数据，权重不再是该点量测数据的原来的权重，故应该对这个权重进行修正。

3 节点权重的修正3．1 二层二叉树节点权重修正考虑一棵简单的二层二叉树，它有3个节点：父亲节点Z以及2个儿子节点X、Y。

3个节点的关系方程为z=x+y。

设3个点的量测平均值分别为x0、y0、z0，量测方差分别为D(x0)、D(y0)、D(z0)。

配电网状态估计的作用是充分利用量测数据，辅以历史负荷数据，利用数据的冗余性，通过一定方法计算出全网状态，剔除不良数据，补足缺乏量测，提高量测数据的一致性和精度，为其他配网高级应用软件提供全面可靠的数据。

该三相状态估计算法是通过对传统的输电系统的状态估计算法进展改良得到，主要采用加权最小二乘算法。

由于量测数据都有一定误差，该算法对量测做了局部简化假定：1〕随机量测噪声的均值为0；2〕量测误差平方的期望值服从标准差为σ的正态分布； 3〕不同量测之间互不相关。

在给定系统网络接线、支路参数和量测系统的条件下，设系统状态变量的个数为n ，量测量的个数为m ，那么反映量测量与系统状态变量之间关系的非线性量测方程可写为(4-1)式中，为量测量向量〔简称量测向量〕，，其中，为系统的第i 个量测量，；为状态变量向量〔简称状态向量〕，，其中，为系统的第i 个状态变量，；为非线性量测函数，描述了量测向量与状态向量之间的关系；为量测误差向量，，其中，为第i 个量测量的量测误差。

根本加权最小二乘状态估计法是以计算得到的状态变量的估计值所对应的量测估计值和量测值之差的加权平方和最小为目标准那么的估计方法。

它是许多状态估计算法的根底方法。

给定量测向量z 和量测方程式，可建立基于加权最小二乘法的系统状态估计的目标函数(4-2)式中，为量测权重矩阵，，其中，为第i 个量测量的权值。

假定量测误差向量服从高斯分布，那么，，其中，R 为量测误差协方差矩阵。

由于通常可假定各个量测误差之间相互独立，因而(4-3)()x =+z h v z 12[,,...,]Tm z z z =z i z 1,...i m =x 12[,,...]T n x x x =x i x 1,...i n =()h x z x v 12[,,...]T m v v v =v i v ()(())(())T J =--x z h x W z h x W 12[,,...]Tm w w w =W i w i z v ()=0E v ()=TE R vv 21222m R σσσ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦式中，为第i 个量测量的量测误差方差，i=1,2,…,m 。

电力系统状态估计算法研究与实现电力系统是现代社会的重要组成部分，负责将发电厂产生的电能传输到用户端。

为了保证电力系统的可靠性和稳定性，精确的状态估计是必不可少的。

本文将就电力系统状态估计算法进行研究与实现，分析其重要性、算法原理、实现方法和应用前景。

一、电力系统状态估计的重要性电力系统状态估计是通过测量值和部分预测值，利用估计算法分析电力系统的运行状态和电力设备的参数，从而提供准确、可靠的状态信息。

电力系统状态估计的重要性体现在以下几个方面：1. 改善系统可靠性：通过准确估计电力系统的状态，可以及时发现故障和异常，并进行相应的调整和处理，从而提高系统的可靠性和鲁棒性。

2. 优化系统运行：基于准确的状态估计结果，可以进行系统优化调度和负荷预测，以最小化系统的能耗和成本，并实现电力资源的合理利用。

3. 实时监控和故障诊断：准确的状态估计可以为系统提供实时监控和故障诊断能力，提高故障的检测和诊断精度，并促进对故障的及时处理和恢复。

二、电力系统状态估计的算法原理电力系统状态估计的算法原理主要包括以下几个方面：1. 滤波算法：电力系统中存在测量误差和噪声，滤波算法可以对测量值进行平滑处理，提高估计结果的准确度和稳定性。

常见的滤波算法包括卡尔曼滤波、粒子滤波等。

2. 优化算法：通过对状态估计问题建立优化模型，采用优化算法求解最优解，优化算法可以提高估计结果的精度和收敛速度。

常见的优化算法包括最小二乘法、最大似然估计、遗传算法等。

3. 数据挖掘算法：电力系统数据量大、维度高，数据挖掘算法可以从大量的数据中挖掘出隐藏的规律和模式，并应用于状态估计问题。

常见的数据挖掘算法包括人工神经网络、支持向量机、决策树等。

三、电力系统状态估计的实现方法电力系统状态估计的实现方法可以分为离线估计和在线估计两种。

离线估计是指通过历史的测量数据和预测数据，使用离线计算方法进行状态估计。

该方法的优点是计算量较小，可以进行详细的数据分析和处理，但缺点是无法满足实时性的要求。

电⼒系统信息理论状态估计电⼒系统信息理论状态估计1.1研究⽬的电⼒系统中的任何控制决定，⽆论是⾃动还是⼿动，均依赖于实时访问系统的实际状态。

由于状态从未直接可⽤，所以有⼀个称为状态估计的过程来估计状态变量，并从SCADA周期性获取的测量中得出其值。

因此，传统上通过求解表征操作点的⾮线性⽅程组来解决问题。

需要实时实现的准确性和鲁棒性对于这⼀过程⾄关重要，因为关于真实系统的知识是能源管理系统（EMS）功能效率的先决条件，⽽来⾃真实的国家可能导致不恰当或甚⾄危险的控制决定。

尽管有数⼗年的合理成功的努⼒，开发了可靠地被业界⽤于⼤规模系统的强⼤的⽅法，但电⼒系统状态估计（PSSE或简称SE）的问题仍然是研究⼈员关注的焦点。

在传输系统实践的多年中，PSSE中已经接受了标准实践：基于加权最⼩⼆乘法（WLS）标准（最⼩均⽅误差准则或MSE的变体）和⽜顿算法作为求解器的回归模型。

在垂直整合的电⼒系统中，⾼⽔平的监控，⾃动化和远程控制传统上限于传输⽹格，⽽配电⽹并不是运营有效性的重点。

只要分配⽹络根据其性质被认为是准静态，即具有良好理解和可预测的负载和电压分布，并且在某些正常情况下不需要重新配置或改变保护和控制设置，这种⽅法是可接受的。

现代电⼒系统概念受到放松管制的电⼒市场状况和电⼒系统对⾮电⼒企业的开放的影响。

因此，通过新的监管框架实现了⼀些基本的业务变化。

这些新的情况特别影响分布式发电⼚正在开始爆炸式安装的分销⽹络，以及对服务质量和运营效率的不断提出的要求。

这是从单向能量流的静态径向配置迁移到⽀持分布式资源和⼀些新类型负载的系统（考虑⼤规模实施电动汽车使⽤的实际趋势）。

这种类型的⽹络⽆疑将需要更复杂的控制和⾃动化，这在⽬前⼏乎不可实施，主要是由于监控基础设施⾮常差。

因此，为了应对新的挑战，⼤量投资主要是将系统监控扩展到全⾯的配电⽹。

需要通过分销管理系统（DMS）来处理这样丰富的新数据流，这将通过复杂的⾃动化应⽤程序满⾜分销系统功能的新要求。

3节电力系统状态估计（WLS算法）3节电系统状态估计报告【任务说明】:闭合的开关:打开的开关:打开的刀闸:线路:负荷G:发电机:母线:连接线(没有阻抗) Unit2Unit13节点系统主接线图任务：1、采用最小二乘状态估计算法,所有量测的权重都取1.0,编写状态估计程序(C/Matlab)。

2、按量测类型，列出量测方程（每一类写出一个方程）3、画出程序流程4、提交源程序，程序中每个函数的作用5、提交计算的输出结果(屏幕拷贝)系统参数：功率基值：100MW电压基值：230 kV线路阻抗参数（标么值）：线路量测(流出母线为正)：母线电压量测：负荷量测(流出母线为正)：发电量测(流入母线为正)：注：量测存在误差【数据预处理】首先根据基值将已知的量测值均转换为标幺值，并将功率值转换为流入量，得到如下数据：线路导纳参数（标么值）：线路注入功率量测(标幺值)：负荷点注入功率量测(标幺值)：发电机节点注入量测(流入母线为正)：母线电压量测（标幺值）：【量测方程】选择节点1的电压相角为参考，为0度，以vi表示误差值。

1）节点1电压量测方程：Vi=Vi+v1即1.0087=V1+v12）1-3支路1号节点处注入有功功率功率：P ij=V i2g ij-V i V j(g ij cos+b ij sin)+v20.613=V12g13-V1V3(g13cos+b13sin)+v2即0.613=-1.6171V12-V1V3(-1.6171cos +13.698sin)+v2 3）1号节点注入功率：P i=V i2G ii +G ij cos+B ij sin+v3P1=V12G11+G1j cos+B1j sin+v3即-1.11=3.5613V12+V1V2(-1.9442cos -10.5107sin)+V1V3(-1.6171 cos -13.698 sin)+v3【流程图】其中iterations 为迭代次数，可见本例的迭代次数为4，收敛较快，状态估计得到的节点1、2、3电压分别为：234.0144444444444444444444444444444444444444444444 【程序说明】，，A X=b,得X+Xmax|X|<1、计算h矩阵的函数cal_hfunction h=cal_h(V,th0,B,G) %其中，V为节点电压估计值，th0为节点电压相角估计%值，B为节点电导矩阵，G为节点电纳矩阵b=-B; %线路电导矩阵g=-G; %线路电纳矩阵P=zeros(3,1); %初始化，节点注入功率PP=zeros(3,3); %线路注入功率QQ=PP;th=[0;th0]; %节点1的电压相角为0for i=1:3P_P=0;Q_Q=0;for j=1:3if(j~=i)P_P=P_P+V(i)*V(j)*(G(i,j)*cos(th(i)-th(j))+B(i,j)*sin(th(i)-th(j)));Q_Q=Q_Q+V(i)*V(j)*(G(i,j)*sin(th(i)-th(j))-B(i,j)*cos(th(i)-th(j)));PP(i,j)=(V(i)^2)*g(i,j)-V(i)*V(j)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));QQ(i,j)=-(V(i)^2)*b(i,j)-V(i)*V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j)));endendP(i)=(V(i)^2)*G(i,i)+P_P;Q(i)=-(V(i)^2)*B(i,i)+Q_Q;endVV=[V(1);V(2);V(3)];h=[P;Q;PP(1,2);PP(2,3);PP(3,1);QQ(1,2);QQ(2,3);QQ(3,1);PP(1, 3);PP(2,1);P P(3,2);QQ(1,3);QQ(2,1);QQ(3,2);VV];2、计算H矩阵的函数cal_HHfunction H=cal_HH(V,th0,G,B,P,Q) %其中，P，Q为根据电压估计值计算得到的节点%注入电压b=-B;g=-G;PV=zeros(3,3); %节点注入功率对电压幅值的偏导数Pth=zeros(3,3); %节点注入功率对电压相角的偏导数Qth=zeros(3,3);PPV=zeros(3,3); %P ij对V j的偏导数QQV=zeros(3,3); %Q ij对V j的偏导数PPth=zeros(3,3); %P ij对th j的偏导数QQth=zeros(3,3); %Q ij对th j的偏导数PPV1=zeros(3,3); %P ij对V i的偏导数QQV1=zeros(3,3); %Q ij对V i的偏导数PPth1=zeros(3,3); %P ij对th i的偏导数QQth1=zeros(3,3); %Q ij对th i的偏导数VV=eye(3);Vth=zeros(3,2);th=[0;th0];for i=1:3for j=1:3if (i~=j)PV(i,j)=V(i)*(G(i,j)*cos(th(i)-th(j))+B(i,j)*sin(th(i)-th(j)));QV(i,j)=V(i)*(G(i,j)*sin(th(i)-th(j))-B(i,j)*cos(th(i)-th(j)));PPV(i,j)=-V(i)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j))); QQV(i,j)=-V(i)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j))); PPV1(i,j)=2*V(i)*g(i,j)-V(j)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));QQV1(i,j)=-2*V(i)*b(i,j)-V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j)));if (j~=1)Pth(i,j)=V(i)*V(j)*(G(i,j)*sin(th(i)-th(j))-B(i,j)*cos(th(i)-th(j))); Qth(i,j)=-V(i)*V(j)*(G(i,j)*cos(th(i)-th(j))+B(i,j)*sin(th(i)-th(j))); PPth(i,j)=-V(i)*V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j))); QQth(i,j)=-V(i)*V(j)*(-g(i,j)*cos(th(i)-th(j))-b(i,j)*sin(th(i)-th(j)));endif(i~=1)PPth1(i,j)=V(i)*V(j)*(g(i,j)*sin(th(i)-th(j))-b(i,j)*cos(th(i)-th(j))); QQth1(i,j)=-V(i)*V(j)*(g(i,j)*cos(th(i)-th(j))+b(i,j)*sin(th(i)-th(j)));endelsePV(i,j)=(G(i,i)*(V(i)^2)+P(i))/V(i);QV(i,j)=(Q(i)-(V(i)^2)*B(i,i))/V(i);if (j~=1)Pth(i,j)=-B(i,i)*(V(i)^2)-Q(i);Qth(i,j)=P(i)-(V(i)^2)*G(i,i);endendendendH=[[PV,Pth(:,2:3)];[QV,Qth(:,2:3)];...[PPV1(1,2),PPV(1,2),0,PPth(1,2),0;...0,PPV1(2,3),PPV(2,3),PPth1(2,3),PPth(2,3);...PPV(3,1),0,PPV1(3,1),0,PPth1(3,1)];...[QQV1(1,2),QQV(1,2),0,QQth(1,2),0;...0,QQV1(2,3),QQV(2,3),QQth1(2,3),QQth(2,3);...QQV(3,1),0,QQV1(3,1),0,QQth1(3,1)];...[PPV1(1,3),0,PPV(1,3),0,PPth(1,3);...PPV(2,1),PPV1(2,1),0,PPth1(2,1),0;...0,PPV(3,2),PPV1(3,2),PPth(3,2),PPth1(3,2)];...[QQV1(1,3),0,QQV(1,3),0,QQth(1,3);...QQV(2,1),QQV1(2,1),0,QQth1(2,1),0;...0,QQV(3,2),QQV1(3,2),QQth(3,2),QQth1(3,2)];...[VV,Vth]];3、主程序calculate_all.m文件format longG=[3.5613,-1.9442,-1.6171;...-1.9442,3.0993,-1.1551;...-1.6171,-1.1551,2.7722]; %B为节点电导矩阵B=[-24.2087,10.5107,13.698;...10.5107,-20.295,9.7843;...13.698,9.7843,-23.4832]; %G为节点电纳矩阵P=[-1.11;0.88;0.23]; %节点注入功率量测值Q=[-0.135;-0.0424;0.24];PP=[0.613;-0.24;-0.459]; %线路1-2,2-3,3-1注入功率在首端的量测值QQ=[-0.012;0.066;-0.165];PP1=[0.467;-0.6;0.24]; %线路1-3,2-1,3-2注入功率在首端的量测值QQ1=[0.148;-0.024;-0.072];V=[1.0087;1.0198;1.0281]; %节点电压幅值量测值R=diag(ones(21,1)); %权重都取为1Z=[P;Q;PP;QQ;PP1;QQ1;V]; %量测值矩阵V0=[1;1;1]; %初值th0=[0;0];delta=100;iterations=0; %迭代次数while delta>0.000001iterations=iterations+1;h=cal_h(V0,th0,B,G); %计算h矩阵H=cal_HH(V0,th0,G,B,h(1:3,1),h(4:6,1)); %计算H矩阵A=H'*inv(R)*H;b=H'*inv(R)*(Z-h);d=A\b; %求解修正值delta=max(abs(d));V0=V0+d(1:3,1); %修正估计值th0=th0+d(4:5,1);enditerationsV0=V0*230; %转换为有名值th0=th0*180/pi; %转换为度for i=1:3j=num2str(i);v=num2str(V0(i));show1=strcat('The voltage magnitude of node ',j,' is', v,' kV'); disp(show1);endfor i=1:2j=num2str(i+1);th=num2str(th0(i));show1=strcat('The phase angle of node ',j,' is ',th,' degrees'); disp(show1);end。