3节电力系统状态估计(WLS算法)
三相电压不平衡度评估的算法原理
（10） 其中n= 1 ,2 ,3 , ..., N-1。则第n 次谐波的幅值Cn为
，当n取1时就可以得到基波的幅值。 但是这里存在一个计算量的问题，也就是实现算法的程序执行时间 问题。考虑x(n)是长度为N的复数序列的一般情况，对某一个k值，直接 计算X(k)值需要N次复数乘法，(N-1)次复数加法。因此，对所有N个k 值，共需 次复数乘法，以及N(N-1)次复数加法运算。当N>>1时，N(n1)≈N2。由上述可见，N点DFT的乘法和加法运算次数均与 成正比。当N 较大时，运算量相当可观。所以，必须减少其运算量，才能使DFT在工 程计算中得到应用。于是J.W.Cooley和J.W.Tukey于1965年根据DFT导出 了快速傅立叶变换算法（FFT）。迄今为止，快速傅立叶变换的发展方 向主要有两个：一个是针对N等于2的整数次幂的算法，如基2算法、基4 算法和分裂基算法等；另一个是N不等于2的整数次幂的算法，它是以 Winograd为代表的一类算法。因为FFT是DFT的一种快速算法，所以FFT 的运算结果必然满足DFT的基本性质。它使用一些算法上的技巧大大减 少了DFT的运算量，使得计算机计算FFT时的速度更快。
但是对称分量法包含电压矢量计算，需要测量三相电压的大小和相 位，这就提高了对仪器性能的要求。于是，有些标准就推荐了一些基于 线电压方均根值的三相电压不平衡度估算法。
2.4.2基于线电压的其它不平衡度计算方法 线电压的特点就是没有零序分量，不必考虑中性点位移。基于线电
电力系统稳态分析课件
a
(2 ln 1 2 ln 3 r
Dab Dbc Dca
r
2
)ia
107
(2 ln
3
Dab Dbc Dca r
r
2
)ia
107
(2 ln
Dm r
r
2
)ia
107
式中 Dm 3 DabDbcDca 称为三相导线的互几何均距,简称
几何均距。
a相导线每米的正序等值电感为：
0.779r
在近似计算中，可以取架空线路的电抗为 0.40/ km
（3）分裂导线三相架空线路的电抗
分裂导线采用了改变导线周围的磁场分布，等效地增 加了导线半径，从而减少了导线电抗。
分裂导线一般是将每相导线分裂为若干根，布置在正 多角形的顶点上，实际应用中分裂数不超过４根。
d
二分裂
2
d12
d 23
ln
1 r
ib
ln
1 Dab
ic
ln
1 Dca
107
(2) a
r
2
ia
2ia
ln
1 r
ib
ln
1 Dca
ic
ln
1 Dbc
107
(3) a
r
2
ia
2ia
ln
1 r
ib ln
1 Dbc
r 导线材料的相对导磁系数，对铜、铝，r 1
f 交流电频率（Hz）
Dm 几何均距（mm或cm），Dm 3 Dab Dbc Dca
将f=50Hz,μr=1代入上式
配电网的状态估计综述与展望
电力系统状态估计算法综述与展望许勇 12031140 电气工程摘要:电力系统状态估计是能量管理系统(EMS)的重要组成部分，状态估计的核心部分是状态估计算法，简要介绍了电力系统状态估计的基本概念和数学模型，阐述了近年来研究较多的几种电力系统状态估计算法，其中包括:最小二乘算法、快速分解算法、正交变换算法、量测变换状态估计算法以及分区协调算法等，并对各个算法的优缺点加以评述，同时，介绍了几种目前新出现的状态估计方法，即新息图状态估计、谐波状态估计、配电网状态估计等。
最后，对状态估计方法中值得研究的方面进行了展望。
关键词：电力系统状态估计能量管理系统算法0 引言电力系统的运行方式和网络结构的日趋复杂，对现代化调度系统提出了必须准确、快速、全面地掌握电力系统实际运行方式和运行状态的要求。
以计算机为基础的现代能量管理系统(EMS)的出现，是电力系统自动化理论与技术上的一次飞跃，实现了从传统的经验型调度到现代化分析型调度的迈进。
电力系统状态估计是EMS的重要组成部分。
EMS的各种高级应用如电压稳定性分析、暂态稳定性分析和安全约束调度等都要依赖状态估计所提供的实时可靠数据[1-4]。
随着数学、控制理论的不断进步，新技术的不断涌现，电力系统状态估计算法也有了较快的发展，其估计精度和实时性也有了实质性的提高，进一步确保了电力系统的安全性和经济性。
状态估计也称为滤波，它是利用实时量测系统的冗余度来提高数据精度，自动排除随机干扰引起的错误信息的手段，估计或预报系统的运行状态(或轨迹)的一种方法[5]。
本文简要介绍了电力系统状态估计的基本概念和数学模型，阐述了包括:最小二乘算法、快速分解算法、正交变换算法、量测变换状态估计算法以及分区协调算法等近年来研究较多的状态估计算法，并对各个算法的优缺点加以评述。
同时，介绍了几种目前新出现的状态估计方法，即新息图状态估计、谐波状态估计、配电网状态估计等。
最后，对状态估计算法的研究进行展望。
内标法的原理
内标法的原理内标法是一种基于功率频率特性的电力系统状态估计方法。
它通过测量电网节点的电压和电流数据，利用功率频率特性的线性关系来推断系统的状态，并用估计值替代未知量，从而实现电力系统状态的在线估计。
内标法的原理是建立在以下几个基本假设的基础上：1. 电力系统是合理模型：内标法假设电力系统可以用定常线性模型描述，在短时间内不会发生剧烈的变化。
2. 电力系统是功率频率特性模型：内标法利用电力系统功率频率特性的线性关系进行状态估计。
功率频率特性告诉我们，当电力系统处于平衡状态时，电压和电流的幅值和相位之间存在一定的线性关系。
3. 电力系统无扰动：内标法假设电力系统没有外界扰动，即电力系统的状态估计仅基于电力系统自身的测量数据。
具体来说，内标法的计算过程如下：1. 测量数据采集：根据测量点的位置，采集电力系统各个节点的电压和电流数据。
2. 准备正常状态基准值：选取一组正常状态的测量数据作为基准值，用于计算功率频率特性模型的系数。
3. 计算高频状态：根据功率频率特性模型，计算出高频状态下测量数据的估计值，其中包括功率、功率相位和电流幅值等。
4. 计算低频状态：根据高频状态的估计值和测量数据，利用线性方程组的解法（如高斯消元法或最小二乘法），计算出低频状态下的估计值，如电压相位和电压幅值等。
5. 更新状态估计：根据低频状态的估计值，更新电力系统各个节点的状态估计值。
内标法的优点是能够在电力系统实时运行过程中进行状态估计，并且不需要额外的测量设备。
然而，内标法也有一些限制，包括对电力系统模型的要求高、对测量数据的准确性要求高等。
因此，在应用内标法时需要仔细选择和处理测量数据，以提高状态估计的准确性和可靠性。
状态估计和潮流计算
状态估计和潮流计算
状态估计
状态估计vs潮流计算
状态估计与潮流计算的已知量的差别
状态估计与潮流的计算⽅法的差别
状态估计与潮流计算的计算结果的差别
调度员潮流的验收标准
（1）能进⾏给定（历史、当前或预想）断⾯的潮流计算；
（2）能模拟发电机出⼒、负荷功率和变压器分接头位置的调节，模拟各种断路器、隔离开关的变位操作，各种模拟操作应能多重组合；
（3）能够模拟发电机的启停和负荷、电容器、电抗器、变压器的投切；
（4）能模拟系统的解列、并⽹、合环、解环；
（5）发电机、负荷开断或调节时应考虑功率缺额或功率过剩在其它可调机组中的分配；
（6）可以计算⽹络解列后若⼲个独⽴⽹络的潮流；
（7）能改变各种控制参数（收敛精度和迭代限次等）进⾏潮流计算；
（8）能进⾏越限报警，可⽤列表或图形显⽰越限信息；
（9）能保存和调⽤潮流计算结果；
（10）运⾏⽅式的改变应能在⼀次接线图上直接进⾏；
（11）潮流计算结果应能在画⾯上直观明了的显⽰和打印输出。
第七章 EMS能量管理系统
解非线性方程7－4－15 求△x(l管理系统
4。P—Q分解法状态估计 P—Q分解法状态估计沿用稳态潮流计算 中P—Q分解法的思路： 高压系统有功功率主要与各结点电压向 量的角度有关； 无功功率主要受各节点电压幅值影响。 将雅可比矩阵进行简化。
第七章 EMS能量管理系统
第七章 EMS能量管理系统
四、自动发电控制 （Automatic Generation Control， 简称AGC） 发电出力与负荷平衡； 保持系统频率为额定值； 使净区域联络线潮流与计划相等； 最小化区域运行成本。
第七章 EMS能量管理系统
（五）不良数据的辨识 除了不良数据点的残差呈现出超过检测 阈值外，还有一些正常测点的残差也超过 阈值，这种现象称为残差污染。 在多个不良数据情况下，由于相互作用 可能导致部分或全部不良数据测点上的残 差近于正常残差现象，这称为残差淹没。
第七章 EMS能量管理系统
（五）不良数据的辨识
三、实时调度与事故预想
（一）静态安全分析 (Security Analysis) 静态安全分析功能对多种给定运行方 式(状态)进行预想事故分析，对会引起 线路过负荷、电压越限和发电机功率越 限等对电网安全运行构成威胁的故障进 行警示，从而对整个电网的安全水平进 行评估，找出系统运行的薄弱环节。
第七章 EMS能量管理系统
第七章 EMS能量管理系统
（三）最小二乘估计法 对于某个状态x的真值进行测量时， 一般都使仪表的测量值z与x之间是函数 关系，即z＝h（x）。任何仪表都有误 差，设此误差为v，那么测量值z与x之间 的关系为： z＝h（x）＋v
第七章 EMS能量管理系统
（三）最小二乘估计法 最小二乘估计状态变量的估计值，对 应的估计值与测量值z之间的误差平方最 小。
状态估计的原理和作用
状态估计的原理和作用1. 原理状态估计（State Estimation）是指通过系统模型及测量数据，利用数学和统计的方法来确定系统的状态变量的一种方法。
它通常用于系统控制和监测中，能够帮助我们实时获得系统的精确状态信息。
状态估计的基本原理可以分为以下几个步骤：1.1 系统建模首先，需要对系统进行数学建模，将系统的动态行为用数学方程描述。
常见的系统模型有线性方程、非线性方程、概率模型等。
1.2 状态方程系统的状态方程描述了状态变量如何随时间变化的关系。
通常采用微分方程或差分方程来表示。
1.3 观测方程系统的观测方程描述了观测变量与状态变量之间的关系。
观测方程通常是状态方程的线性组合，但也可以是非线性方程。
1.4 测量数据通过传感器等设备，获取系统的测量数据。
测量数据可以是离散的样本数据，也可以是连续的时间序列数据。
1.5 估计方法基于系统模型和测量数据，利用数学和统计推断方法，推导出系统的状态估计方法。
常见的状态估计方法有最大似然估计、卡尔曼滤波、粒子滤波等。
1.6 状态估计根据估计方法，将测量数据代入系统模型，计算出系统的状态变量的估计值。
估计值可以是离散的时间序列，也可以是连续的曲线。
2. 作用状态估计在实际应用中起着重要的作用，具体包括：2.1 系统监测状态估计可以实时准确地监测系统的状态信息，帮助我们了解系统的运行情况。
例如，在航空航天领域，状态估计可以用于检测飞行器的姿态、速度等状态变量，以确保飞行器的稳定和安全。
2.2 系统控制状态估计可以提供准确的状态信息，用于系统控制。
通过与控制算法结合，可以实现对系统的准确控制，提高系统的性能。
例如，在自动驾驶领域，状态估计可以用于估计车辆的位置和速度，从而实现智能驾驶。
2.3 故障诊断状态估计可以用于故障诊断，帮助我们快速准确地判断系统是否发生故障，并找出故障原因。
通过与故障诊断算法结合，可以实现对系统故障的自动检测和诊断。
例如，在工业生产中，状态估计可以用于监测设备的运行状态，及时发现故障。
状态估计方法
状态估计方法概述状态估计是指通过观察和测量数据，推断系统当前状态的一种方法。
它在许多领域中都有广泛的应用，包括工程、金融、医学等。
基于数学模型的状态估计方法一种常见的状态估计方法是基于数学模型的方法。
这种方法通过使用已知的系统模型和测量数据，利用数学运算来估计出系统的状态。
卡尔曼滤波器卡尔曼滤波器是一种常用的基于数学模型的状态估计方法。
它通过使用系统的动态模型和传感器的测量数据，计算出系统当前状态的最优估计。
卡尔曼滤波器能够处理含有噪声和不确定性的测量数据，提供较为准确的状态估计结果。
扩展卡尔曼滤波器扩展卡尔曼滤波器是卡尔曼滤波器的一种扩展版本。
它适用于非线性系统和非线性观测模型。
扩展卡尔曼滤波器通过对系统和观测模型进行线性化，以及使用卡尔曼滤波器相似的算法来进行状态估计。
基于数据驱动的状态估计方法除了基于数学模型的方法，还有一种基于数据驱动的状态估计方法。
这种方法不依赖于系统的数学模型，而是通过大量的数据样本来研究系统的状态和变化规律，从而进行状态估计。
神经网络方法神经网络方法是一种常用的基于数据驱动的状态估计方法。
通过构建和训练神经网络模型，可以利用大量的输入数据来估计系统的状态。
神经网络方法能够处理非线性系统和高维数据，并具有较强的适应性和泛化能力。
支持向量机方法支持向量机方法是另一种常用的基于数据驱动的状态估计方法。
它通过构建并训练支持向量机模型，利用输入数据来进行状态估计。
支持向量机方法适用于分类和回归问题，并能够处理非线性系统和高维数据。
总结状态估计方法是推断系统当前状态的重要工具。
基于数学模型的方法可以利用系统模型和测量数据进行状态估计，包括卡尔曼滤波器和扩展卡尔曼滤波器。
基于数据驱动的方法通过大量数据样本来学习和推断系统状态，包括神经网络方法和支持向量机方法。
在实际应用中，可以根据具体问题的特点选择合适的状态估计方法来获得准确的结果。
电力系统稳态分析计算
标称电压220kV、长度 的单回输电线路， 例题 标称电压 、长度200km的单回输电线路，已知末 的单回输电线路 ~ 端负荷 SL =120 + j10 MVA，电压 2=215kV，线路参数为： ，电压U ，线路参数为： r1＝0.108Ω/km，x1＝0.427Ω/km，b1＝2.665×10－6S/km， Ω ， Ω ， × ， 试求线路始端的电压和功率。 试求线路始端的电压和功率。
阻抗中的功率损耗
2 2 P2 + Q2 1202 + 2.322 1202 + 2.322 ~ P2 + Q2 ∆S = 2 2 R + j 2 2 X = ×21.6 + j ×85.4 2 2 U2 U2 215 215
= 6.73+ j26.61 (M A V )
~ ~ ~ S1 = S2 + ∆S =120 − j2.32 + 6.73+ j26.61 (MVA)
线路始端的电压
U1 = (U2 + ∆U2 ) + jδU2 = (U2 + = 215 +
•
P R + Q2 X P X −Q2 R 2 )+ j 2 U2 U2
120×21.6 + (-2.32)×85.4 120×85.4 − (−2.32) ×21.6 +j 215 215 = 226.13+ j47.9 (kV )
• ∗ ~ S2 =U2 I 2
(MVA)
2 P2 + Q2 S2 I2 = = 2 U2 U2
阻抗中的功率损耗
2 2 P2 + Q2 P2 + Q2 ~ 2 2 ∆S = ∆P + j∆Q = I2 R + jI2 X = 2 2 R + j 2 2 X (MVA) U2 U2
基于TWR的WLS和KF融合室内定位方法
基于TWR的WLS和KF融合室内定位方法目录一、内容概要 (3)1.1 研究背景与意义 (3)1.2 国内外研究现状 (4)1.3 论文结构安排 (6)二、相关理论和技术概述 (7)2.1 定位技术综述 (8)2.1.1 常见定位技术分类 (10)2.1.2 室内定位技术特点 (11)2.2 TWR技术原理 (12)2.3 WLS算法介绍 (13)2.3.1 WLS基本概念 (13)2.3.2 WLS在定位中的应用 (15)2.4 卡尔曼滤波器(KF)原理 (16)2.4.1 KF基本概念 (17)2.4.2 KF在动态系统中的作用 (18)三、基于TWR的WLS和KF融合定位模型设计 (19)3.1 融合模型需求分析 (20)3.2 WLS与KF融合策略 (21)3.3 融合模型数学描述 (22)3.4 模型实现关键技术 (23)3.4.1 数据预处理 (24)3.4.2 参数估计 (25)3.4.3 模型校正 (26)四、实验设计与分析 (27)4.1 实验环境搭建 (29)4.2 实验方案设计 (30)4.3 实验数据采集 (31)4.4 实验结果分析 (32)4.4.1 定位精度评估 (33)4.4.2 算法性能对比 (34)4.4.3 结果讨论 (36)五、系统实现 (37)5.1 系统架构设计 (38)5.2 关键模块开发 (39)5.2.1 数据处理模块 (41)5.2.2 定位计算模块 (42)5.2.3 用户界面模块 (43)5.3 系统测试与验证 (44)5.3.1 测试方案 (46)5.3.2 测试结果 (47)5.3.3 系统优化 (49)六、结论与展望 (50)6.1 研究总结 (51)6.2 存在的问题及改进方向 (52)6.3 未来工作展望 (53)一、内容概要本文档提出了基于传播时差融合的室内定位方法，该方法首先通过TWR技术获取来自信号源的时延差信息，利用这些差值计算接收节点到信号源的相对距离，以构建位置估计模型。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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。