隐式改进欧拉法在暂态稳定分析中的应用

蒋德珑 尹淑萍 王克文 代旭光

(郑州大学电气工程学院，河南 郑州 450001)

0 引言

经济和技术的进步使电力系统得到了迅速发展，但也对电力系统稳定性提出了更高的要求［1］。多年来，大量文献研究了电力系统的专业化和自动化，在电力系统暂态稳定自动化分析中更是取得了长足的进步［2-11］。如采用计算出扰动后各发电机转子间相对角度随时间的变化曲线(即摇摆曲线)来判断系统的稳定性［6-9］;对暂态稳定程序进行模块化处理，使程序具有较好的扩展性和灵活性［10-11］。

相对显式解法，隐式解法有其明显的优越性。在总结前人研究的基础上，本文将经典的改进欧拉法与隐式解法相结合，用于电力系统的暂态稳定计算，以有效地改善计算的精度和速度，提高暂态稳定自动化分析的水平，为电力系统的安全稳定运行提供更为有力的保障。

1 牛顿法潮流计算

牛顿法是数学中解决非线性方程式的典型方法，它通过对非线性方程逐次线性化求解方程的根，并能保持较好的收敛性。潮流计算是研究电力系统暂态稳定的基础，它根据给定的发电运行方式和系统接线方式来确定整个电力系统各部分的运行状态，包括各母线的电压、各元件流过的功率和系统的功率损耗等，为暂态稳定分析提供初始运行状态。当前，牛顿法被广泛应用于电力系统潮流计算，其计算的核心问题是修正方程式的建立和求解［12］。本文采用经典的牛顿法潮流计算方法，其求解过程大致分为以下几步:首先给定节点导纳矩阵YB和节点电压初值e(0)、f(0)，代入功率偏差量的表达式中，求解 ΔP(0)、ΔQ(0)、(ΔV2)(0)，再求出修正方程式的系数矩阵(雅可比矩阵)各元素;接着通过解修正方程式得到修正量Δe(0)、Δf(0)，然后对节点电压 e(1)=e(0)-Δe(0)、f(1)=f(0)-Δf(0)进行修正，进而利用 e(1)、f(1)求得 ΔP(1)、ΔQ(1)、(ΔV2)(1);最后判断结果是否收敛，如收敛，则求出各支路潮流并输出结果，否则以e(1)、f(1)为初值重新迭代。

2 稳定分析数学模型的确定

整个电力系统的模型，在数学上可以统一描述成如下一般形式的微分-代数方程组:式中:x为微分方程组中描述系统动态特性的状态变量;y为代数方程组中系统的运行参量。本文重点关注发电机微分方程的求解。

2.1 发电机数学模型

考虑到暂态稳定的计算精度和系统的复杂性，本文采用忽略定子电磁暂态但考虑转子阻尼绕组作用的五阶发电机模型(E'q、E″q、E'd、δ、ω)，即考虑 f绕组、D绕组、Q绕组的电磁暂态以及转子运动的机电暂态。

转子电压平衡方程为:

转子运动方程为:

式中:T'd为励磁绕组本身的时间常数;T″d为励磁绕组短路、定子绕组开路时d轴阻尼绕组的时间常数;T″q为定子绕组开路时q轴阻尼绕组的时间常数;TJ为发电机组的惯性时间常数;E'q为发电机q轴的暂态电动势;E″q为发电机q轴的次暂态电动势;E″d为发电机d轴的次暂态电动势;δ为功角;ω为转子角频率。将式(2)与式(3)联立，即可得到同步发电机的五阶数学模型。

2.2 隐式改进欧拉法

对于一般的非线性微分方程式:

确定时间步长 Δt=h，由(tn，xn)点推算(tn+1，xn+1)点时，改进欧拉法的递推公式即为:

改进欧拉法的核心思想是取某时段始点导数值与终点导数值的平均值作为该时刻折线段的斜率，进而得到更加精确的计算结果。而用f(xn，tn)和f(xn+1，tn+1)代替积分曲线［tn，tn+1］区间始点和终点的斜率，即可将式(5)进一步变换为:

式(6)中只有xn+1为未知数，因而利用求解代数方程式的方法就可以求得xn+1，这就得到了隐式改进欧拉法。

传统的改进欧拉法是利用已知量，根据递推公式逐次递推计算出相应时段终点的函数值，因此是一种显式解法。隐式改进欧拉法不是给出递推公式，而是首先把微分方程化为差分方程，然后利用求解差分方程的方法确定函数值，其特点就是把微分方程的求解问题转换成一系列代数方程的求解过程。因此，隐式解法有两个明显的优点:一是当微分方程和代数方程需要联立求解时，利用隐式解法便于消除交接误差;二是隐式解法可以采取较大的步长。

2.3 差分方程

对于电力系统暂态计算涉及到的发电机、电动机和控制器等数学模型中的微分方程，皆可以将其化为差分方程［13-14］。下面以同步发电机的数学模型为例，利用隐式改进欧拉法，将五阶发电机模型中的微分方程化为差分方程。

由此得到的用于描述各发电机、电动机、控制器等差分方程组，可与用于描述电力网络的代数方程组联立求解，发挥了牛顿法所具有的数值稳定、无交接误差和收敛速度快等特点。

3 初值计算及运行参数的求解

根据上述潮流计算法，可以计算出电力系统遭受扰动以前的正常运行状态，得到电力网络的运行参数，并由此求出电力系统中有关元件的状态参数初始值。这些都是进行电力系统暂态稳定分析的先决条件。

3.1 初值的计算

在进行暂态稳定分析之前，必须求出暂态过程计算所需要的初值条件，即首先要确定求解微分方程所需要的初值。

在简化的暂态稳定计算程序中，初值的计算包括求系统扰动瞬间发电机的暂态电势、转子角度、原动机的机械功率和等值导纳等。这些参数在电力系统扰动的瞬间是不会突变的［10］。因此，可以由扰动前的正常运行状态计算得到。

3.2 与网络方程联立求解

在电力系统中，电力网络将系统中看起来相互独立的所有暂态元件联系在一起，在暂态过程中的任一时刻，各暂态元件注入网络的电流不但由其本身的特性决定，且整个电力网络必须满足基尔霍夫定律。其中，前者由各暂态元件自身的代数方程描述，后者反映在电力网络方程中。因此，为了求解网络方程，需要列出各暂态元件自身的代数方程，并对其进行处理，从而与网络方程联立求解。将发电机作为电流源，其节点注入电流的表达式为:

对于负荷模型，只需对导纳矩阵的相应对角块进行变化，即可得到新的网络方程。

4 算例分析

利用上述方法编写暂态稳定计算程序，对典型的9节点三机系统进行仿真计算，三机系统示意图如图1所示。

图1 三机系统示意图Fig.1 Three-machine system

图1中，系统有3台发电机、3个负荷和9条支路，频率为60 Hz，支路数据和发电机参数见文献［14］。

假设系统在0 s受到扰动，在线路5-7靠近母线7处发生三相接地短路，则故障在5个周波(约0.08333 s)后由断开线路5-7而被消除。

本文采用牛顿法潮流计算得出正常运行情况下的系统潮流，如表1所示。以其作为初始运行条件，由隐式改进欧拉法暂态稳定分析程序可计算各发电机功角在故障过程中的变化情况。

表1 潮流计算结果Tab.1 Results of load flow calculation

在暂态稳定计算程序中，各发电机用五阶模型模拟，各负荷用恒定阻抗模拟，电力网络用导纳矩阵描述，微分方程用隐式改进欧拉法求解，网络方程用直接法求解。程序计算了从短路故障开始后，系统在2 s内的暂态过程，得出计及发电机凸极效应时各发电机的功角和相对摇摆角，其中数值积分采用0.001 s的步长。

系统故障过程中各发电机的功角数据如表2所示，计算时间步长取 0.01 s。表 2 中，1-1、2-2、3-3分别表示发电机①、②、③的功角。限于篇幅，表2中只截取了以0.2 s为间隔的数据。

表2 功角数据Tab.2 Data of power angle

系统故障过程中各发电机功角的相对角度如表3所示，即发电机相对摇摆角，计算时间步长取0.01 s，其中，下标2-1表示发电机②相对于发电机①的角度，3-1表示发电机③相对于发电机①的角度。限于篇幅，表3中只截取了以0.4 s为间隔的数据，“*”表示极值时刻。

表3 相对摇摆角数据Tab.3 Data of relative swing angle

由表2和表3可知，虽然发电机的功角不断增大，但其相对摇摆角的变化并不大，第一摆时的最大摇摆角度 为 δ2-1=126.47259°(t=0.59s)、δ3-1=94.32671°(t=0.64 s)，第二摆时的最大摇摆角度为δ2-1=126.74058°(t=1.88 s)、δ3-1=95.87463°(t=1.92 s)，比第一摆时的角度稍大，2 s内各发电机间的最大相对摇摆角均小于180°。

为了便于分析，利用Matlab7.1对2 s内电力系统发生三相短路故障的功角暂态过程进行仿真，发电机功角和相对摇摆角仿真结果如图2和图3所示。图2表示系统故障过程中各发电机的功角变化曲线，反映出3台发电机的功角逐渐增大，说明有功出力逐渐增大。图3表示系统故障过程中各发电机的功角的相对变化，反映出在计及发电机凸极效应时，各发电机间的最大相对摇摆角均小于180°。通过与图2的比较可以发现，虽然3台发电机的功角逐渐增大，但其相对摇摆角的变化并没有超过180°，说明系统在发生该故障时仍能保持稳定，即系统是暂态稳定的，与表2和表3中数据反映的情况一致。

算例结果证明，该方法能够较好地应用于电力系统暂态稳定的判定过程，同时用隐式解法对微分方程的处理，也更有利于在暂态计算过程中对各变量进行修正，使发电机功角的计算结果更加接近真实值。

5 结束语

本文采用隐式改进欧拉法对电力系统暂态计算中各元件的数学模型进行处理，通过联立求解的方法进行暂态稳定计算，最后对电力系统多机运行方式下发电机功角的暂态过程进行了分析。

与传统的改进欧拉法相比，该方法具有积分精度高和无交接误差的优点，可以满足长期仿真计算对数值精度的要求。同时，本文介绍的方法保持了电力系统内在的模块性，在程序设计上具有良好的可扩展性与灵活性。由于选择了较为理想的系统，忽略了励磁系统和调相机对电力系统暂态的影响，采用线性化模型得到的功角摇摆曲线还存在一定的误差，因此，下一步可以采用保留高阶项的模型，以求进一步改善计算精度。

［1］倪以信，陈寿孙，张宝霖.动态电力系统的理论和分析［M］.北京:清华大学出版社，2002.

［2］孙晖，赵菁.建立专业电力自动化监控系统的探讨［J］.自动化仪表，2003，24(6):15-18.

［3］唐雯，罗安，唐杰，等.电力监控系统图形绘制及数据处理一体化平台［J］.自动化仪表，2007，28(8):56-58.

［4］张海生，范征宇.S7-300 PLC和VC++.NET在电力监控系统中的应用［J］.自动化仪表，2006，27(8):50-52.

［5］程若发，冯士芬，江晓舟.基于ARM的同步发电机故障录波系统的设计［J］.自动化仪表，2007，28(9):11-14.

［6］西安交通大学，清华大学.电力系统计算［M］.北京:水利电力出版社，1978.

［7］Al-Taee M A，Al-Azzawi F J，Al-Taee A A，et al.Real-time assessment of power system transient stability using rate of change of kinetic energy method［C］//IEEE Proceedings of Generation Transmission and Distribution，2001，148(6):505-510.

［8］吴颖，赵岩，蒋传文，等.风电场穿透功率极限计算方法及发展［J］.自动化仪表，2008，29(11):7-11.

［9］赵艳雷，童建忠.改进欧拉法在电力系统暂态分析中的应用和软件设计［J］.微计算机信息，2004，20(5):98-99.

［10］Kasztenny B，Kezunovic M.A method for linking different modeling techniques for accurate and efficient simulation［J］.IEEE Transactions on Power Systems，2000，15(1):65-72.

［11］Suyono H，Nor K M，Yusof S.Component based development for transient stability power system simulation software［C］//Proceedings of Power Engineering Conference，PECon，2003:16-21.

［12］陈珩.电力系统稳态分析［M］.2版.北京:中国电力出版社，1995.

［13］武东亚，王克文，马洁.基于概率潮流的多运行方式下发电机功角的暂态过程分析［J］.电气应用，2009，28(6):64-67.

［14］王锡凡，方万良，杜正春.现代电力系统分析［M］.北京:科学出版社，2003.