基于Matlab语言的电力系统最小二乘法状态估计毕业论文.doc

重庆邮电大学本科毕业设计（论文）基于Matlab语言的电力系统最小二乘法状态估计毕业论文目 录前 言1第一章 电力系统状态估计概述2第一节 电力系统状态估计的发展历史2第二节 电力系统状态估计的主要内容3第三节 状态估计的发展方向6第四节 论文主要内容及章节安排7第二章 算法基础8第一节 数据结构8一、三角形表8二、链接表9第二节 排序算法的实现10第三节 稀疏矩阵求解的实现10第四节 节点编号优化11第三章 状态估计13第一节 状态估计的数学模型13一、状态估计的量测方程13二、量测误差方差矩阵16三、状态估计准则17第二节 状态估计的基本方法18一、无约束加权最小二乘法18二、正交变换法18三、混合法19四、带约束的加权最小二乘法19五、Hachtel矩阵法20第三节 算法比较20第四章 最小二乘法状态估计22第一节 算法计算流程22第二节 IEEE4节点算例分析24第三节 IEEE30节点测试系统仿真30一、初始参数30二、仿真实验31三、结果分析33结 论34致 谢35附录37- 36 -前 言随着电力系统的迅速发展，电力系统的结构和运行方式日趋复杂，电力系统调度中心的自动化水平也需要逐步由低级向高级发展。现代化的调度系统要求能迅速、准确而全面地掌握电力系统的实际运行状态，预测和分析系统的运行趋势，对运行中发生的各种问题提出对策，并要提供下一步的决策。从而保证电力系统运行的安全性和经济性。在现代的调度系统中，计算机已经成为最重要的一环。计算机的高级自动化功能主要体现在它所具备的程序的功能。高级在线应用程序的特点是要对大量实时数据进行处理与分析，以确定电力系统的安全与经济运行状况，因此保证电力系统实时数据的质量是进一步提高计算机在线应用水平的关键。为了建立可靠而完整的实时数据库，通常有两种途径：从硬件的途径可以增加量测设备和运动设备，并提高其精度、速度与可靠性；从软件的途径，可以采用现代的状态估计技术，对数据进行实时处理。但是对测量与运动设备提出过高的要求会导致技术和经济上付出过大的代价。如果在具备一定水平的硬件基础上，采用状态估计技术则能充分发挥已有硬件设备的潜力，提高数据的精度，补充测点和量测良母的不足，排除偶然的错误信息和数据，提高整个数据系统的质量与可靠性。状态估计也被称为滤波，他是利用实时测量系统的冗余度来提高数据精度，自动排除数据干扰所引起的错误信息，估计或预测系统的运行状态。状态估计作为近代计算机实时数据处理的手段，首先应用于宇宙飞船、卫星、导弹、潜艇和飞船的跟踪、导航和控制中。电力系统状态估计的研究也是由卡尔曼滤波开始的，但是根据电力系统的特点，即状态估计的主要处理对象是某一时间段上的高维空间问题，而且对测量误差的统计知识又不够清楚，因此目前很多电力系统实际采用的状态估计算法是最小二乘法。第一章 电力系统状态估计概述第一节 电力系统状态估计的发展历史在电力工业发展初期，发电厂都建在用户附近，电厂规模较小，电力系统也是简单而孤立的。运行人员在发电机、开关设备等电力元件的近旁直接监视设备状态并进行手工操作，例如人工操作开关、调节发电机的出力和电压等。这种工作方式的效果与运行人员的素质和精神状态有关，往往不能及时而正确地进行调节和控制。特别是在发生事故时，往往来不及对事故的发生和发展做出反应而使事故扩大。随着工农业生产和人民生活用电的增长，电力系统内的发电设备及其出力不断增加，供电范围也不断扩大。在这种情况下，设备现场人工就地监视和操作已不能满足电力系统运行的需要了。为了保证电力系统安全运行和向用户供应合格电能，出现了单一功能的自动装置。这些装置有故障自动切除装置(即继电保护装置，自动切除出现故障的发电机、变压器和输电线路等设备)、自动操作和调节装置(如断路器自动操作、发电机自动调压和自动调速装置等)和远距离信息自动传输装置(即远动装置)。为了提高电力系统供电的可靠性和运行的经济性，逐步地将孤立的电力系统连接起来发展成了跨地区的电力系统。由于电力系统中每座发电厂和变电站的运行值班人员只知道本厂(站)的运行情况，对系统内其它厂(站)的运行情况及电力系统的运行结构不清楚，所以在跨地区的电力系统形成之后，就必须建立一个机构对电力系统的运行进行统一管理和指挥，合理调度电力系统中各发电厂的出力并及时综合处理影响整个电力系统正常运行的事故和异常情况。这个机构就是电力系统调动所，也称电力系统调度中心。随着电力系统的迅速发展，电力系统的结构和运行方式日趋复杂，电力系统调度中心的自动化水平也需要逐步由低级向高级发展。现代化的调度系统要求能迅速、准确而全面地掌握电力系统的实际运行状态，预测和分析系统的运行趋势，对运行中发生的各种问题提出对策，并决定下一步的决策。从而保证电力系统运行的安全性和经济性1。但是，电力系统遥测设备经常受随机误差、仪表误差等误差之患，因此用这样粗糙的系统行为信息来判断系统状态，显然是不能满足要求的。对系统状态的估计是控制的必要条件，因此要改变系统状态，首先要知道它处于什么状态。然而，已被广泛应用于飞机和宇航系统的数据分析和估计理论，直到六十年代末七十年代初才开始应用于电力系统的在线数据处理。1968年丰田淳一作出了用卡尔曼滤波方法做负荷预报和水库来水预报的文章，它已经属于状态估计在电力系统中应用的研究。然而状态估计在电力系统中被广泛研究和实际应用，却是针对实时潮流问题进行的。按照目前习惯的说法，“电力系统状态估计”一词的含义就是指实时潮流的状态估计2。1969年美国麻省理工学院的许怀丕(F.C.Schweppe)等人提出了基本加权最小二乘法的状态估计，其特点是收敛性能好，估计质量高。然而由于这种算法的计算量和使用内存比较大，难以用于大型电力系统的实时计算。之后，H.P.Horisberger等人吸取潮流计算经验而建立的快速分解状态估计算法，兼顾了计算速度、收敛性、使用内存和对各种类型测量量的适应性等方面的优点，可以看成是基本加权最小二乘法状态估计的实用形式。接着，美国电力公司(American Electric Power)的道帕兹恩(J. F. Dopazo)等人提出了测量变换估计算法，它也属于最小二乘法的总体算法，其特点是仅用支路潮流测量值，计算速度快、使用内存少和程序简单，虽然难以处理结点注入型测量量，但并不妨碍其实用性，在1975年就投入了实际运行。在同一时期，美国邦那维尔电力系统的拉森(R. E. Larson)等人提出了卡尔曼型的逐次估计算法，但由于电力系统状态量的维数较高，不得不采用对角化的状态估计误差协方差矩阵，因此这样虽然有节省内存和提高计算速度的优点，却降低了收敛性能和估计质量而妨碍了实用性。其后在美国的其它电力公司以及挪威、瑞典、日本、法国、英国、澳大利亚、意大利和前苏联等国相继开展这方面的研究工作。最早应用状态估计程序的是挪威水利电力局(Tokle)所属的较小的电网和美国电力公司(AEP)所属的较大的电网；至70年代末80年代初，世界上约有十几个电网在正常运行中使用了状态估计程序。状态估计在电力系统中所得到的效果己被肯定，新设计的电力系统调度中心都应包含这一新的功能。自70年代末开始，我国北京、广东和华东等电力系统先后与有关科研机构和高等院校合作开展了状态估计课题的研究工作。80年代初北京电力系统进行了状态估计的实时试验。第二节 电力系统状态估计的主要内容电力系统的各种遥测、遥信信息是通过远动装置转送到调度中心的，由于远动装置的误差及在传送过程中各个环节所造成的误差，使这些数据存在不同程度的误差和不可靠性。此外，由于量测装置在数量上或者种类上的限制，往往不可能得到完整的、足够的电力系统分析所需要的数据。为解决上述问题，除了不断改善量测与传输系统外，还可以采用数学处理的方法来提高量测数据的可靠性和完整性。因此，电力系统状态估计就是为适应这一需要而提出来的。从掌握电力系统运行情况的要求来看，总是希望能由足够多的测量信息通过远动装置送到调度中心，但从经济性与可靠性来看，只能要求将某些必不可少的信息送到调度中心，通常称足够表征电力系统特征所需要最少数目的变量为电力系统的状态变量。电力系统状态估计就是要求能在测量量有误差的情况下，通过计算以得到可靠的并且为数最少的状态变量值3。为了满足状态估计计算的上述需要，对电力系统的量的测量在数量上要求有一定的裕度。通常将全系统中独立量测量的数目与状态量数目之比，称为冗余度。只有具有足够冗余度的量测条件，才可能通过计算机以状态估计算法来提高实时信息的可靠性与完整性，建立实时数据库。由于电力系统远动装置的工作情况是会经常变化的，当远动信息量严重不足时，状态估计无法工作。因此，在状态估计之前需要进行可观测性检验。如果系统中某些部分被判为是不可观测的，无法通过状态估计建立实时数据库，则应把它从状态估计的计算中退出来，或者用增加人工设置的虚拟量测量或者称为伪量测数据来使它变成可观测的4。协同状态估计进行工作的是不良数据的检测与辨识，如果有误差很大的，一般没有随机性的数据，就应该将它剔除，并重新进行状态估计，最终建立起完整的电力系统实时数据库。由于电力系统状态估计必须在几分钟内完成，因此它通常可以跟踪节点负荷的变化规律，在必要时可用来提供补充的量测量。因此，状态估计的计算结果也可以用于负荷预测。电力系统状态估计的整个功能流程框图，如图1-1所示。图1-1 电力系统状态估计功能流程框图由此可见，作为状态估计的核心部分-状态估计计算，可以根据量测系统量测量的时域界定，将状态估计算法划分为动态和静态两种：动态状态估计算法考虑的是不同时刻下的量测量之间的联系与影响，静态状态估计计算则仅对同一时刻端面下的量测量进行估计分析，从而确定系统的状态变量。由于受到实际系统的运行限制，如数学模型的维数很大、通道传送量少、传送速度慢以及测点时间难于同步等原因，使得动态状态估计目前仍处于理论研究阶段，未真正投入实际使用。本文以下所述状态估计，如无特别说明，均指静态估计。电力系统状态估计的基本步骤如下所示，一般包括：模型假设、状态估计5、检测6、和辩识。(1) 模型假设：是指在给出网络接线状态和网络参数的条件下，确定量测函数方程和量测误差方阵的过程。(2) 状态估计：是计算状态估计值的过程，即是使残差的加权内积达到最小的状态值。(3) 检测：即检查量测值中是否存在不良数据或网络接线状态中是否存在错误信息的过程。(4) 辩识：是确定具体不良数据或网络接线错误的过程。电力系统量测误差来源大体分为以下几种：(1) PT/CT的误差；(2) A/D转换引起的误差；(3) 数据传输时受到干扰引起的误差；(4) 采集数据的不同时性引起的数据误差；(5) 运行中三相不平衡及功率因数的变化引起的误差；数据的不齐全性主要是由于RTU分散不均或太少引起的状态估计主要完成以下功能：(1) 根据网络方程的最佳估计标准(一般为最小二乘准则)，对生数据进行处理,以得到最接近于系统真实状态的最佳估计值，提高数据精度；(2) 对生数据进行不良数据的检测和辨识，剔除或修正不良数据；(3) 推算出完整而精确的电力系统的各种电气量；(4) 根据遥测量估计电网的实际开关状态,纠正偶然出现的错误的开关状态信息，以保证数据库中电网接线方式的正确性；(5) 可以应用状态估计算法以利有的数据预测未来的趋势和可能出现的状态。这些预测的数据丰富了数据库的内容，为安全分析与运行计划等程序提供必要的计算条件；(6) 如果把某些可疑或未知的参数作为状态量处理时，也可以用状态估计的方法估出这些参数的值；(7) 通过状态估计程序的离线模拟试验，可以确定电力系统合理的数据采集与传送系统，即确定合适的测点数量及其合理分布；综上所述：电力系统的状态估计程序输入的是低精度、不完整、存在不良数据的产生，而输出的则是精度高、完整可靠的熟数据。第三节 状态估计的发展方向状态估计是当代电力系统能量管理系统(Energy Management System，EMS)的重要组成部分，尤其在电力市场环境中发挥更重要的作用7。状态估计问题的提出激发了学者的研究兴趣，他们以数学、控制理论和其它新理论为指导，根据当代的计算机软件和硬件条件，结合电力系统的特点，在理论方面进行了大量研究。同时，以状态估计软件使用为目标，针对实际工程面临的问题，探索和总结出许多可行的宝贵经验。状态估计的理论研究促进工程应用，而状态估计软件的工程应用也推动了状态估计理论的研究和发展。然而，状态估计领域仍然存在不少问题未得到妥善解决，随着电力系统规模的不断扩大，电力工业管理体制向市场化迈进，电力系统监控规模不断扩大和各种新理论、新技术的不断涌现，无论从理论方面还是从实际应用需求方面，状态估计领域仍有许多问题需要深入研究8。状态估计领域在以下方面有重要的研究价值：基于GPS相位测量角技术的实时状态估计问题9；面向大系统，开发计算速度快、数值稳定性好的算法，缩短状态估计的执行周期问题；多种类型和多个相关坏数据的检测和识别问题，各类坏数据的特征抽取问题10；测量误差相关情况下的状态估计问题；抗差估计理论应用于电力系统状态估计的进一步研究问题11；新理论应用于电力系统状态估计的理论探讨和实用化的可行性研究问题12。第四节 论文主要内容及章节安排本文主要以MATLAB为工具来实现电力系统状态估计方面的研究与分析。全文共分四章，结构安排如下：第一章介绍了课题研究背景、主要内容以及状态估计未来的发展方向。第二章介绍了状态估计的算法基础，主要包括数据结构、稀疏矩阵的求解以及节点编号的优化。第三章阐述了状态估计的数学模型，介绍了几种常用的状态估计算法，并比较了几种算法的优劣性。第四章运用最小二乘法，结合IEEE4节点对最小二乘法状态估计的流程进行了分析，并将最小二乘法运用到IEEE30节点中验证其有效性。第二章 算法基础电力网络方程指的是将网络的有关参数和变量及其相互关系归纳起来所组成的、可反映网络性能的数学方程式。常用的是节点导纳矩阵的节点电压方程其中是节点导纳矩阵。它对角线上的元素是自导纳(与该节点相连接的所有支路导纳的和)，非对角线上的元素为互导纳(两个节点之间支路的导纳的负值之和)。电力系统的计算分析最终都归结为迭代求解一个线性代数方程组。由于电力系统网络各节点往往只与其中几个节点有支路相连,而与其他大部分节点无关联,因而导致导纳矩阵中存在很多零元素,使矩阵变得很稀疏。因此在实际应用中应采用稀疏矩阵技术。稀疏矩阵技术是求解线性代数方程组的一种极其重要的手段。稀疏矩阵技术有以下几个显著的特点13：节省存储单元。由于只存储矩阵的非零元素,对于稀疏性很大的矩阵可以节省大量的存储单元。尤其是在进行电力系统的潮流分析和状态估计运算中。在运算过程中由于只处理非零元素，避免了像或这样的无效运算，因而节省了运算时间。在求解过程中通过对方程组未知变量消去次序的适当选择，使矩阵在变换过程中一直保持稀疏性。第一节 数据结构要实现稀疏矩阵技术,首先就要采用一个好的存储稀疏矩阵的数据结构，既要节省存储单元，又要便于对非零元素的检索和处理。存储稀疏矩阵的数据结构一般采用了两种：一种是静态存储，另一种是动态存储。静态存储主要考虑便于访问而不考虑存储单元的变化，通常在稀疏矩阵的行列排序已确定，对方程组求解时使用；动态存储的最大特点是便于非零元素的插入和删除。通常在稀疏矩阵的行列排序中使用。下面简述一下常用的两种数据结构。一、三角表三角表同时使用了行指针/列标和列指针/行标两种存储原理。它适合于LU分解，是一种实用的静态存储格式。由于LU分解过程对于上三角阵U来说是按行进行的，对于下三角阵L来说是按列进行的，因此U阵中的非零元素采用行指针/列标格式存储，L阵中的非零元素采用列指针/行标格式存储。对于给定的矩阵，所有的非零元素在存储时按下述次序编号：首先从左上角开始到右下角为止顺序给阶数个对角元素编号。然后，按行列交替次序地给非零非对角元素编号。通常整个三角形表由四个一维数组组成：VAL(K)：存放第K号元素的数值；ROCO(K)：存放第K号元素的行列号(若为对角元素)，列号(若为U阵元素)或行号(若为L阵元素)；URP(K)：存放U阵第K行的行指针；LCP(K)：存放L阵第K列的列指针。二、链表链表是一种重要的动态存储格式,其特点是矩阵中的非零元素不需要按任何特定的次序来存放，它除了采用行列号和行列指针来存储每个非零元素在矩阵中的位置信息外，还设置了专门的指针来指示该元素的相邻元素的位置，使每行和每列的非零元素分别连接成链。数据连接成行链和列链后，使得对非零元素的寻找，插入和删除十分方便。文中使用的链表是双链表。一般的双链表由九个一维数组组成：VAL(K)：存放第K号元素的数值；ROW(K)：存放第K号元素的行号；COL(K)：存放第K号元素的列号；.UP(K)：存放第K号元素的上邻元素编号,若无上邻元素则置零；DOWN(K)：存放第K号元素的下邻元素编号,若无下邻元素则置零；LEFT(K)：存放第K号元素的左邻元素编号,若无左邻元素则置零；RIGHT(K)：存放第K号元素的右邻元素编号,若无右邻元素则置零；RP(K)：存放矩阵第K行行链的行指针；CP(K)：存放矩阵第K列列链的列指针。第二节 排序算法的实现在进行LU分解运算时有可能使矩阵中原来的零元素变为非零元素。我们称这样的非零元素为填入(fillin)。矩阵的稀疏性可能使求解线性代数方程组的运算量大大减少，但是填入的出现却使矩阵的稀疏性削弱，不但增加了矩阵占用的存储单元，而且增加了后续LU分解的运算量。因此我们总希望填入数越来越少。矩阵的行列排序可以减少LU分解过程中的填入数。在建立稀疏矩阵类库时，采用了四种排序方法：（1）田尼一沃克(Tinner-Walker)算法该算法对LU分解过程中的每一步主元的选择原则是：所选择的主元使该步分解产生的填入数最少。若同时有几个候选主元都满足上述要求时，则可以选取其中的一个作为该步的主元。（2）马柯维兹算法该算法对LU分解的每一步主元的选择原则是:所选择的主元使该步分解所需的长运算(LOP)数最少。若同时有几个候选主元都满足上述要求时，则可以选其中的一个作为该步的主元。（3）全主元消去法该算法对LU分解的每一步主元的选择原则是:在执行到第K步时，要把矩阵中包含第K行和第K列右下角部分的所有元素作为后选主元进行搜索，找出绝对值最大者作为该步的主元。（4）列主元消去法该算法对LU分解的每一步主元的选择原则是:在执行到第K步时，要把第K列中包含第K个元素以下的所有元素作为后选主元进行搜索，找出绝对值最大者作为该步的主元。第三节 稀疏矩阵求解的实现排序过程结束后，网络的节点己经按照稀疏矩阵的要求重新编号，这时便可以建立网络方程并进行求解。方程的求解包括LU分解和前消后代两个部分，其总的要求是对非零元素进行运算。由于在排序过程中所有的填入已经确定并插入，在求解方程时矩阵非零元素的分布结构不再变化，因此宜采用静态存储格式的数据结构。在本文中使用的是三角表。稀疏矩阵求解过程分为以下几步：（1）符号LU分解：LU分解的主要工作量是修正余下的子矩阵的元素。当矩阵结构给定时，每一步分解所需要修正的元素是一定的，所以在数值求解之前应该先进行一次模拟分解,确定每步分解所需修正元素的编号，预先将它们存入数组LUP中，以后每次分解相同结构的矩阵时，就可以根据数组LUP中存储的编号，直接从数组VAL中取出所需修正的元素的数值进行运算，这样可以避免多次重复的数据检索。这种模拟分解建立数组LUP的过程就称为符号LU分解。（2）数值LU分解：这一步才实现真正的理论上的LU分解。LU分解的每一步都包含两类修正运算。在第K步时，第一类修正就是：将U阵中第K行的所有非零元素除以主元值VA)。第二类修正就是：对主行主列右下角子矩阵的所有非零元素的修正。数值LU分解过程只需对非零元素的数值进行运算，不需要任何对非零元素的检索操作。（3）前消后代：由于许多网络问题需要对结构相同但元素值不同的向量进行多次重复求解方程组的运算，因此与LU分解类似，对于前消后代过程也引入了符号前消和符号后代的概念。所以本文中的前消后代过程包括了四个部分：确定并插入在前消过程中右端向量产生的填入；符号前消，建立前消需修正元素的编号表；数值前消；数值后代。右端向量填入的确定和插入，在三角形表的数据结构中，右端向量的非零元素紧接着矩阵的最后一个非零元素,按自上而下的次序编号并存放在数组VAL中，同时在数组ROCO中存放它们的行号。在进行前消后代过程前应先确定一下右端向量的填入。第四节 节点编号优化节点编号优化就是找到一种网络节点重新编号方案，使按此构成的节点导纳矩阵或雅可比矩阵分解产生的注入元大大减少，尽量保持原有的稀疏度14。节点编号优化的方法有如下三种：（1）静态优化法-按静态联结支路数的多少编号这种方法最简单。它是先统计好网络中各节点联结的支路数后，按联结支路数的多少，由少到多，顺序编号。当有节点联结支路数相同时，可以按任意次序对这些结进行编号。这种编号方法的依据是，在导纳矩阵中联结支路数最少的节点所对应的行中非零元素也最少，因此在分解过程中产生注入元素的可能性也比较小。（2）半动态优化法-按动态联结支路数的多少编号这种方法最常用。运用这种方法时，先只编一个联结支路数最少的节点号，并立即将其消去；再编消去第一个节点后联结支路数最少的节点号，再立即将其消去，依此类推。所以要这样是由于消去某节点后，可能出现新增支路而使余下节点联结的支路数发生变化，不宜一次将所有节点号都编好。（3）动态优化法-按动态增加支路数的多少编号这种方法不常用。运用这种方法时，不首先进行节点编号，而是首先寻找消去后出现的新支路数最少的节点,并为其编号，且立即将其消去;然后再寻找第二个消去后出现的新支路数最少的节点，为其编号，也立即将其消去，依此类推。这样可以保证逐个消去节点时出现的新支路数(即注入元数)最少。显然，同一网络按这三种方法所编节点号往往不相同。最严格的方法是动态优化法，但是其计算量比半动态优化法大得多。第三章 状态估计电力系统状态估计15一般包括网络拓扑分析、可观测性检验、估计计算和不良数据的处理四个基本步骤16。其中状态估计算法是状态估计程序的核心部分。本章首先介绍了在给出网络接线和网络参数的条件下确定量测函数方程和量测误差方差阵的过程，接着介绍了估计算法中最常用的方法类型：最小二乘算法17, 18。具体分析了几种常见最小二乘算法的特点，指出了他们的优缺点及其适应范围。第一节 状态估计的数学模型一、状态估计的量测方程电力系统的运行状态可以用节点电压模值、电压相角、线路有功与无功潮流、节点有功与无功注入等物理量来表示。状态估计的目的就是应用经量测量得到的上述物理量通过估计计算来求得能表征系统运行状态的状态变量。电力系统静态运行的状态变量，通常取节点电压模值与电压相角。当有一个平衡节点时，N个节点的电力系统状态变量维数为n=2N1。如果假定电气接线与参数都已知，根据状态变量不难求取各个支路的有功、无功潮流及所有节点的注入量测。在估计中，状态变量需要借助量测方程式，即联系状态向量与量测向量之间函数关系来间接求得。在考虑有量测噪声式，它们之间的关系可以写成： 式中：为维的量测量向量；为量测函数向量。其中： 为量测噪声向量，其表达式为： 很容易写出状态变量x与支路潮流的非线性函数表达式，称为节点电压量测方程式；也可以写出节点注入量测功率与支路潮流的非线性函数表达式，称之为功率量测方程式。表3-1列出五种基本的量测方式。第一种量测其维数为2N1，显然没有任何冗余度，这在状态估计中是不实际的。第五种量测方式具有最高的维数和冗余度，但是所需要的投资太高，也是不现实的。因此，实际电力系统量测方式是第一种到第四种的组合。表 3-1 五种基本量测方程测量方式z的分量方程式h(x)z的维数（1）除平衡节点外所有节点的注入功率式（3-4）式（3-5）2N1（2）除了（1）的量测外再加上所有的节点的电压模值式（3-4）式（3-5）式（3-14）3N1（3）每条支路两侧的有功、无功潮流式（3-6）式（3-7）4M（4）除了（3）的量测外再加上所有的节点的电压模值式（3-6）式（3-7）式（3-8）4M+N（5）完全的量测系统式（3-4）式（3-15）4(M+ N)1注：为节点数；为支路数。表3-1中的各种方程式，当用图3-1中所标的量并以直角坐标形式表示时，节点注入功率方程式为： 由节点到节点的支路潮流为： 上四式中：、和分别为节点电压的不合虚部；、和为图2-1所示的p形路线原件模型中的参数；和为导纳矩阵元素。 图 3-1 p形线路元件模型图图 3-2 形线路元件导纳模型图图 3-3 变压器等值电路图交流电力系统中的潮流方程也可以用极坐标形式表示。线路的等值电路如图 3-2所示。节点注入量测功率量测的极坐标形式为： 上式中： ，约定、分别为节点、的电压幅值；和为导纳矩阵元素。支路-上节点侧线路潮流的极坐标表示形式为： 如变压器的等值电路如图 3-3 所示。变压器支路侧潮流方程的极坐标形式为： 式中：为变压器非标准变比。为标准侧，变比为1；为非标准变比侧，变比为；为变压器标准测的电纳，有，其中为变压器标准侧电抗。、和的关系如下： 二、量测误差方差矩阵用量测量来估计系统的状态存在若干不确定或者不精确的因素，概括起来有以下内容：（1）数学模型不完善。测量数学模型中通常往往包含有工程性的近似处理。除此以外，还可能存在模型中所采用参数不精确的问题，还有当网络结构变化时，所采用的结构模型不能及时更新。上述问题中属于参数不精确的，通常可用参数估计方法来解决；属于网络结构错误的，则采用网络接线错误的检测与辨识来解决。（2）测量系统的系统误差。这是由于仪表不精确，通道不完善所引起的。它的特点是误差恒为正或负且没有随机性。一般这类数据属于不良数据。清除这类误差的方法，主要是依靠提高测量系统的精确性与可靠性，也可以用软件方法来检测与辨识，然后找出不良数据，并通过增加量测系统的冗余度来补救，但这仅是一种辅助手段。（3）随机误差。这是量测系统中不可避免的。其特点就是小误差比大误差出现的概率大，正负误差出现的概率相等，即概率密度曲线对称于零值或误差的数学期望为零。在状态估计式（3-1）和（3-3）中的误差向量就是指的这种误差。测量的随机误差或噪声向量v是均值为零的高斯白噪声，由于不同时间的测量之间是不相关的，而且在一般情况下，不同测量的误差之间也是不相关的。误差的概率密度或者协方差很难由测量或计算来确定，因此在实际应用中常用测量设备的误差来确定。记每个测量误差的方差为。测量误差的方差阵，可以写成每个测量误差方差的对角阵。 各个量测值不可避免地带有随机误差，量测值与被量测的物理量的真值之间总是有差异的。即使被测量的物理量没有变化，重复测量得到的量测值也是不会完全相同的。如果根据理想的量测方程,由量测的量测值来求取系统的状态量，并假定量测方程是线性的。这样由量测量来求解状态量就是解线性方程的问题。一般量测量的维数大于状态量的维数，即方程数大于未知量数。由于量测量的误差，线性方程组存在矛盾方程而无解。但这样的系统仍然是可观测的，虽然不能直接解方程组，但可以用拟合的办法根据带误差的量测量求出系统状态在某种估计准则意义下的最优估计值。所谓优化总是对一定的目标函数来讲的。对于给定的目标函数，当状态量的估计值为最优时，目标函数取极值。最小方差估计、极大验后估计和极大似然估计，这三种估计方法都是统计学的估计方法，虽然有较好的估计质量，但是都要求事先掌握较多的随机矢量的统计特性。这些要求在电力系统状态估计的实际计算中是不容易做到的，因此也是难于实现的。三、状态估计准则状态估计准则是指求解状态变量的原则，电力系统采用的估计准则大多是极大似然估计，即求解状态变量使量测值被观测到的可能性最大。根据给定不同的目标函数，可以得到不同的估计准则。目前有以下估计准则：加权最小二乘准则(WLS)、非二次准则(non-quadratic)、加权最小绝对值(WLAV)、LMS (least median of squares)和LTS (least trimmed squares)。加权最小二乘准则假设量测量严格服从正态分布；LMS、LTS假设量测量服从拉普拉斯分布；非二次准则假设量测量服从Huber；分布。目前应用最多的是WLS估计准则。它的优点是模型简单，计算量小，对理想正态分布的量测量，估计具有最优一致且无偏等优良统计特性。在量测方程式中，由于量测误差是一个随机变量，所以量测量也是一个随机变量。当给定量测矢量以后，加权最小二乘方(WLS)状态估计的目标函数： 状态估计就是求解满足(3-17)式的状态向量。第二节 状态估计的基本方法一、无约束加权最小二乘法要使目标函数达到最小值，加权最小二乘法(3-17)式所给出的状态估计量必须满足极值条件，即 其中 为雅可比矩阵非线性方程(3-18)式的解可以通过将其线性化再迭代求解线性化的方程求得。 得： 求解该式首先形成增益矩阵并对其进行三角分解求得因子表。二、正交变换法线性加权最小二乘方式在每次迭代中的目标函数为： 简化得： 设为一正交矩阵，即，为单位阵，使得 其中为一上三角矩阵，则(3-23)式可进一步化为 当 时，取得最小值。三、混合法由(3-23)式可得到增益矩阵 这样方程(3-19)式可变为 混合法就是用正交因子迭代求解方程(3-27)。四、带约束的加权最小二乘法实际系统中往往存在T型接线。T型接线点处的虚拟量测量Z，其量测方程为 虚拟量测量的量测误差为零，其权重在理论上应为无穷大，其中i为T型接线点的节点号。为了避免因权重值相差悬殊而导致方程数值不稳定等式约束。根据线性规划理论，引入拉格朗日乘子向量，则加权最小二乘法状态估计的拉格朗日函数为 五、Hachtel矩阵法定义残差方程 Hachtel方法是把残差方程式作为等式约束增广到带约束的加权最小二乘法中，采用一组更大和更稀疏的方程组来描述被求解问题，以提高数值稳定性。则带约束加权最小二乘方状态估计的朗格朗日函数变为： 第三节 算法比较本节对上一节所述的状态估计的五种算法从以下几方面进行比较：(1) 数值稳定性。(2) 计算效率。(3) 实现复杂性。1. 正交变换法和无约束加权最小二乘法方程系数矩阵的条件数是方程病态条件的量度，由于无约束加权最小二乘法中增益矩阵的条件数是正交变换法中H阵条件数(大于l)的平方，所以正交变换法的数值稳定性比无约束加权最小二乘法好。在(3-23)式中H阵的QR因子化比无约束加权最小二乘法中的增益矩阵G的因子化需要更多的计算量。此外，正交变换法不能利用有功和无功解耦有效地实现，而无约束加权最小二乘法的快速解耦算法非常有效，并在能量管理系统中得以广泛应用。解耦算法的主要优点是只需在进入循环迭代前因子化一次。正交矩阵Q是非稀疏的，并且它的维数会很高，对于量测冗余度较大的系统尤其如此。保留Q阵将需要很大的存贮量，如果在每次迭代时都形成QR因子，就失去解耦算法在计算上的优越性，因此，正交变换法的解耦算法是不可行的。2. 混合法和正交变换法混合法和正交变换法的差别在于两个方程(3-25)和(3-27)的求解过程不同。后一个方程的右端向量分为两部分，一部分为遥测量，一部分为虚拟量测量，设虚拟量测量的权重值是遥测量权重值的倍，当很大时，第二部分占优势，第一部分将会在迭代过程中被当作舍入误差去掉，在方程接近收敛时更是如此，这对方程的数值稳定性很不利，所以混合法的数值稳定性比正交变换法差。此外，混合法不需存贮Q阵，而且能够很容易地用快速解耦算法实现。3. 带约束的最小二乘法和无约束加权最小二乘法由于把虚拟量测量作为等式约束处理，带约束的最小二乘法避免了因权重值相差悬殊而产生病态的可能性，从数值稳定性来看，带约束的最小二乘法比无约束加权最小二乘法好。带约束的最小二乘法的系数矩阵是不定的，这就必须使用较为复杂的因子化方法。4. Hachtel方法和带约束加权最小二乘法Hachtel方法和带约束加权最小二乘法方法的系数矩阵都是对称的，但都是不定的，都需要比较复杂的优化和排序方法。Hachtel方法系数矩阵的维数比带约束加权最小二乘法方法高，但这对于系数矩阵而言并不重要,所以两种方法的计算量相当。5. Hachtel方法和混合法Hachtel方法和混合法是两种在数值稳定性和计算效率方面较好的状态估计方法，他们都可以用快速解藕算法有效地实现。文献通过实验对这两种方法进行了比较,比较结果表明Hachtel方法的数值稳定比混合法要好，而且Hachtel方法的计算效率也较混合法高。通过综合比较，程序使用了加权最小二乘法。第四章 最小二乘法状态估计加权最小二乘法状态估计的基本算法类似于电力系统潮流计算中的牛顿法，可以直接用于较小的电力系统，而适用于大型的电力系统快速分解状态估计算法是由此演化而来的。因为基本算法的估计质量高，可以作为各种状态估计算法的比较基准。第一节 算法计算流程基本加权最小二乘法状态估计程序框图示于图4-1中，其运行步骤如下所述。图4-1 基本加权最小二乘法状态估计程序框图步骤（1）：程序初始化当在线运行的状态估计程序初始启动时，以及在接线状态出现变化时应进行程序初始化。程序初始化的内容有：对状态量赋初值；节点次序优化；形成节点导纳矩阵和分配内存等。在线跟踪估计是在网络接线不出现变化的条件下，将前一次采样的状态估计值作为下一次采样状态估计的初值，这样程序无需进行初始化，直接进入步骤（2）。步骤（2）：输入测量数据进行一次量测采样，由缓冲区读入数据。步骤（3）：回复迭代计数器：。步骤（4）：由现有的状态量x计算各测量量的计算值和雅可比矩阵。这里的状态变量可以选用极坐标或者直角坐标来表示，根据测量量的类型和地点自动选取相对应的公式进行计算。由z和计算出残差r和目标函数值，并由雅可比矩阵计算信息矩阵和自由矢量。步骤（5）：解线性方程组，求状态修正量，并选取其中绝对值最大者，记为。线性方程组的系数矩阵是对称矩阵，可以选用平方根分解法或LDL分解法来求解，国外使用平方根分解法较多。对于较大的电力系统应该采用稀疏矩阵的程序技巧，以提高计算的速度。步骤（6）：收敛检查达到收敛标准以内结束计算，转出口；否则转到步骤（7）继续迭代。步骤（7）：修正状态量运用公式，并将迭代计数器加1：，返回步骤（4）继续迭代。为了避免无休止的迭代，对迭代次数应该加以限制。以上的步骤（3）步骤（7）组成了一个子过程，既可以独立进行状态估计，也可以作为某些不良数据辨识的子程序。正确的状态估计结果将通过一个输出程序送入到数据库中，根据需要可以送出各节点的电压幅值和相角，各节点的注入功率和各支路的潮流计算等数据。第二节 IEEE4节点算例分析用四节点的简化系统作为算例来说明电力系统的运行模拟。网络示意图如图4-2所示，相关参数见表4-1。图4-2 四节点试验系统阻抗图系统由4个节点和4条支路组成，1、4号节点是电厂，2、3号节点是负荷，13号支路是输电线，4号支路是变压器。表4-1 支路数据表支路号起点i终点j电阻R电抗X对地电纳变比K1120.522.660.00014-2133.058.370.00056-3230.412.150.00044-4340.007.50-1.05下面按照图4-1的步骤详细介绍IEEE4节点的基本加权最小二乘法状态估计的计算过程。在初始化阶段，由测量量模拟得到的各测量量的标准差、权重和一次采样的随机测量值列入表4-2中，节点编号的结果列在表4-3中。表4-2 测量数据（电压单位为V，功率单位为MW）测点号测量项目测量值测量误差真值标准差权重1P118.71-2.1420.861.300.582Q134.05-1.7435.791.400.503P1241.791.4540.331.430.484Q1237.96-0.3938.351.420.495P13-19.100.36-19.471.290.596Q13-2.240.31-2.561.180.717P2117.55-2.0119.571.290.598Q21-10.740.37-11.121.240.649P32132.572.14130.422.030.2410Q2342.161.2641.351.440.4811P2-49.490.50-50.001.500.4412Q2-41.93-1.93-40.001.430.4813P14-193.226.77-199.992.490.1614Q14-71.08-0.85-70.221.630.3715V1112.150.65111.501.110.8016V2110.99-0.71111.711.110.80表4-3 节点内外编号对照表节点内部编号1234节点外部编号4231在下面的计算过程中，除了输入和输出部分外，一律使用节点内部编号。按照内部编号，我们可以写出系统的节点导纳矩阵如表4-4所示。表4-4 节点导纳矩阵节点12341-j0.13330j0.12690200.1563-j0.8103-0.0855+j0.4487-0.0707+j0.36213j0.1269-0.0855+j0.44870.1240-j0.6742-0.0384+j0.105440-0.0707+j0.3621-0.0384+j0.10540.1092-j0.4668将参考节点的测量电压作为初始电压，个节点由同一电压开始计算，选取初始状态变量的电压幅值为112.15kV，选取初始状态变量的电压相角为0rad。首先读入一组测量数据，如表4-2所示；置迭代计数器的初值k=1；在第一次迭代过程中，需要根据前面章节提到的公式计算出雅可比矩阵H、测量量的计算值h。根据测量量的类型和地点信息，可以在第三章中找到相对应的计算公式，便可以由现有的状态值和导纳矩阵计算出测量量i所对应的一行雅可比矩阵元素Hi和测量量i的计算值，其中i=1,2,16，便可以求出整个雅可比矩阵H和测量量的计算值h，分别列于表4-5和表4-6中。表4-5 第一次迭代过程中的雅可比矩阵H状态量测点1234567100-7.93-4555.10-4.31-1326.7312.24200-40.61890.45-11.82483.4652.28300-7.93-4555.10007.93400-40.61890.450040.5850000-4.31-1326.764.3160000-11.82483.4611.70700004.311326.76-4.318000011.70-483.46-11.82900-9.59-5645.679.595645.6701000-50.381076.615