基于多级高阶辛Runge-Kutta方法的暂态稳定性并行计算方法.pdf

第３９卷第１１期２０１１年６月１日电力系统保护与控制ＰｏｗｅｒＳｙｓｔｅｍＰｒｏｔｅｃｔｉｏｎａｎｄＣｏｎｔｒｏｌＶｏ１．３９Ｎ０．１１Ｊｕｎ．１．２０１１基于多级高阶辛Ｒｕｎｇｅ－Ｋｕｔｔａ方法的暂态稳定性并行计算方法汪芳宗，何一帆（三峡大学电气与新能源学院，湖北宜昌４４３００２）摘要：将级２阶的辛Ｒｕｎｎｇｅ－Ｋｕｔｔａ方法用于电力系统暂态稳定性计算，利用矩阵分裂技巧以及矩阵求逆运算的松弛方法，导出了一种新的暂态稳定性并行计算方法，具有较好的时间并行特性和超线性收敛性。利用ＩＥＥＥ１４５节点系统，对导出的并行算法进行了仿真测试和评估。仿真测试结果表明，所提出的并行算法具有很好的收敛性，有效地解决了时间并行度与收敛性之间的矛盾，可以获得较高的加速比和很好的并行计算效率。关键词：暂态稳定性；辛几何方法；并行算法；矩阵分裂；松弛牛顿法Ｐａｒａｌｌｅｌｃｏｍｐｕｔａｔｉｏｎｏｆｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙｂａｓｅｄｏｎｍｕｌｔｉ・－ｌｅｖｅｌｈｉｇｈ・－ｏｒｄｅｒｓｙｍｐｌｅｃｆｉｃＲｕｎｇｅ・－ＫｕｔｔａｍｅｔｈｏｄＷＡＮＧＦａｎｇ－ｚｏｎｇ，ＨＥＹｉ・ｆａｎ（ＥｌｅｃｔｒｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ＆ＲｅｎｅｗａｂｌｅＥｎｅｒｇｙＳｃｈｏｏｌ，ＣｈｉｎａＴｈｒｅｅＧｏｒｇｅｓＵｎｉｖｅｒｓｉｔｙ，Ｙｉｃｈａｎｇ４４３００２，Ｃｈｉｎａ）—Ａｂｓｔｒａｃｔ：Ｉｎｔｈｉｓｐａｐｅｒ，ｔｈｅ－ｓｔａｇｅ２ｓ－ｏｒｄｅｒｓｙｍｐｌｅｃｔｉｃＲｕｎｇｅＫｕｔｔａｍ￣ｈｏｄｉｓａｄｏｐｔｅｄｆｏｒｎｕｍｅｒｉｃａｌｃａｌｃｕｌａｔｉｏｎｏｆｐｏｗｅｒｓｙｓｔｅｍｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙ．ＢｙａｎａｒｔｆｕｌｓｐｌｉｔｔｉｎｇｏｆＪａｃｏｂｉａｎｍａｔｒｉｘａｎｄｕｓｉｎｇｒｅｌａｘａｔｉｏｎｔｅｃｈｎｉｑｕｅｏｆｍａｔｒｉｘｉｎｖｅｒｓｅ，ａｎｅｗｐａｒａｌｌｅｌａｌｇｏｒｉｔｈｍ——ｆｏｒｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙｃｏｍｐｕｔａｔｉｏｎｈａｓｂｅｅｎｄｅｒｉｖｅｄ．Ｔｈｅｐｒｏｐｏｓｅｄａｌｇｏｒｉｔｈｍｉｓｏｆｃｏｍｐｌｅｔｅｐａｒａｌｌｅｌ－ｉｎｔｉｍｅ，ａｎｄｈａｓｓｕｐｅｒｌｉｎｅａｒｃｏｎｖｅｒｇｅｎｃｅｒａｔｅ．Ｆｏｒｔｅｓｔ，ＩＥＥＥ１４５－ｂｕｓｐｏｗｅｒｓｙｓｔｅｍｉｓｕｓｅｄ，ａｎｄｔｈｒｏｕｇｈｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｔｈｅｐｒｏｐｏｓｅｄａｌｇｏｒｉｔｈｍｈａｓｂｅｅｎ—ｃｏｍｐａｒｅｄｗｉｔｈｔｈｅｃｏｎｖｅｎｔｉｏｎ￣ｐａｒａｌｌｅｌｉｎ－ｔｉｍｅＮｅｗｔｏｎａｐｐｒｏａｃｈ．Ｔｈｅｔｅｓｔｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅｐｒｏｐｏｓｅｄａｌｇｏｒｉｔｈｍｈａｓｇｏｏｄｃｏｎｖｅｒｇｅｎｃｅｒａｔｅ，ｃａｎｒｅｓｏｌｖｅｔｈｅｃｏｎｔｒａｄｉｃｔｉｏｎｂｅｔｗｅｅｎｔｈｅｐａ—ｒａｌｌｅｌ－ｉｎｔｉｍｅｄｅｇｒｅｅａｎｄｔｈｅｃｏｎｖｅｒｇｅｎｃｅｒａｔｅ，ａｎｄｃａｎｏｂｔａｉｎｈｉｇｈｓｐｅｅｄｕｐａｎｄｐａｒａｌｌｅｌｃａｌｃｕｌａｔｉｏｎｅｆｆｉｃｉｅｎｃｙ．Ｋｅｙｗｏｒｄｓ：ｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙ；ｓｙｍｐｌｅｃｔｉｃｇｅｏｍｅｔｒｙａｌｇｏｒｉｔｈｍ；ｐａｒａｌｌｅｌａｌｇｏｒｉｔｈｍ；ｍａｔｒｉｘｓｐｌｉｔｔｉｎｇ；ｒｅｌａｘｅｄＮｅｗｔｏｎｍｅｔｈｏｄ中图分类号：ＴＭ７１２文献标识码：Ａ—文章编号：１６７４．３４１５（２００９）２０００２２．０５０引言并行计算技术是实现大规模电力系统实时和超实时计算这一目标的有效途径Ｌ１Ｊ。有关电力系统暂态稳定性并行计算方法，大致可以分为三类：空间并行、时间并行以及时间一空间并行。空间并行类方法又可以分为两类：一类是基于网络分割或矩阵分块、分裂思想的粗粒度空间并行方法，例如波形松弛法；另一类是基于矢量计算的细粒度空间并行方法。粗粒度的空间并行方法适用于多ＣＰＵ结构，而细粒度的空间并行方法只能采用矢量或向量处理器。国内有关空间并行计算方法的研究，绝大多数是采用粗粒度的空间并行。时间并行算法是通过同时求解多个积分步长来实现并行求解。很易理解，时间并行是一种粗粒度的并行。将时间并行与空间并行相结合，就是时间一空间并行计算。空间并行算法从理论上讲比较简单；时间并行计算方法比较复杂，因为在电力系统暂态稳定性计算中，待求变量在每一时刻的状态强烈依赖于前一时刻的状态，亦即不同时刻的状态变量彼此强烈相关。因此，时间并行算法研究工作中所存在的主要问题，是如何有效解决算法的时间并行度与收敛性之间的矛盾。近年来，研究人员已提出了不少新的数值积分方法，其中一类是著名的辛几何算法（以下简称辛算法）［４－６］。辛算法是由我国已故著名学者冯康先生及其研究小组，针对传统的Ｒｕｎｇｅ．Ｋｕｔｔａ方法不能保持Ｈａｍｉｌｔｏｎｉａｎ系统的辛几何结构以及具有人为耗散等缺点而提出的。这一新方法的提出，为Ｈａｍｉｌｔｏｎｉａｎ系统的分析计算，同时也为微分方程数值计算方法的研究开辟了一个全新的领域。本文将—多级、高阶辛ＲｕｎｇｅＫｕｔｔａ方法（以下简称辛ＲＫ方法）用于电力系统暂态稳定性的计算，利用多级、高阶辛ＲＫ方法所具有的内在时间并行特性，导出汪芳宗，等基于多级高阶辛Ｒｕｎｇｅ・Ｋｕａａ方法的暂态稳定性并行计算方法－２３－了一种新的暂态稳定性并行计算方法。１辛Ｒｕｎｇｅ＝Ｋｕｔｔａ算法自冯康先生在上世纪８０年代后期首创Ｈａｍｉｌｔｏｎｉａｎ系统的辛算法后】，辛算法的研究已取得了十分丰富的成果，研究人员已提出了多种辛差分格式即辛算法【６】口限于主题，本文只介绍辛ＲＫ方法－。对给定的常微分方程初值问题＝ｆ（ｔ，），ｘ（ｔ。）＝。（１）—其级的ＲｕｎｇｅＫｕｔｔａ方法的一般形式为∑∈ｌＹ＝＋厂（ｆ（１，）｛ｊ＝ｌ（２）ｌＸｎ＋∑一＋ｂ，厂（）Ｌｊ＝ｌ∈式中：ｈ为积分步长；Ｊ，，（（１，））为真解（ｒ）在，』＝ｔ＋ｃｈ时的逼近，即Ｙ（ｆ＋ｃ，）。该ＲＫ方法也可用Ｂｕｔｃｈｅｒ表喘表示，这里∈…∈Ａ＝（），Ｃ：（ＣＩＣ２，，ｃｓ）Ｒ…’∈ｂ＝（６Ｉ，６２，，ｂ）Ｒ其中∑≠∑∈ｂ，＝１，ｂ０；ａ＝Ｃｉｆ（１，）定义Ｂ＝ｄｇａｇ（ｂ，）．Ｍ＝ＢＡ＋ＡＢ－ｂｂ‘∈』ｈｉ＝ｔｉ－－ｔ，＿ｌ，ｆ（１，Ｓ＋１）‘１ｔｉ＝＋ｃ『，＝ｔｎ。＝＋因此，１个Ｓ级的全隐辛ＲＫ方法具有内在和“”真实的Ｓ个时间并行度，而且相当于用传统的积分方法同时积分Ｓ＋１步，如图ｌ所示。为叙述方“”便，将后一种特性称为等效时间并行度，即Ｓ级“”的全隐辛ＲＫ方法具有Ｓ＋１个等效时间并行度。ｌ蹦．¨．［二ｔｔｌｔ２…ｔｔ图１Ｓ级辛ＲＫ方法的内在时间并行特性Ｆｉｇ．１ＩｌｌｕｓｔｒａｔｉｏｎｏｆｐａｒａｌｌｅｌｉｓｍｆｏｒｓｙｍｐｌｅｃｔｉｃＲＫｍｅｔｈｏｄ若用隐式梯形积分方法同样以变步长同时在Ｓ＋１个时间点上进行求解，则这种并行算法仍然只具有２阶计算精度，因为隐式梯形积分方法本身只是２阶算法。但是，上述Ｓ级的全隐辛ＲＫ方法则具有２阶的精度。依据局部截断误差，若用隐式梯形积分方法以固定步长，同时在Ｓ＋１个时间点上进行求解，在不考虑误差累积的前提条件下，这种并行求解的总体截断误差为Ｄ（ｈ３）；若采用Ｓ级２阶的全隐辛ＲＫ方法以步长ｈ积分１步，则其截”断误差为Ｏ（ｈ）。因此，在保持２种计算方法计算精度大致相近的条件下，即Ｏ（ｈ）Ｄ（｜ｉ２），则有ｈ（）。举一个简单的例子：取＝０．０５ｓｅｃ．，Ｓ＝３，则ｈ０．２７７＞４ｈ也就是说，一个３级的辛ＲＫ方法积分ｌ步至少可以等同于用隐式梯形积分法同时或连续积分４步。这就是多级高阶辛ＲＫ方法所具有的内在时间并行特性。多级相当于多步；每增加１级就相当于增加了１个时间并行度。洙古要耋暑曼经证明ｎ伽：当且仅当＝０时，上２暂态稳定性并行计算方法述方法是辛算法。‘。日心’。工刀¨¨开具体的辛ＲＫ方法有２种结构形式：当为下三角阵时，方法称为对角隐式辛ＲＫ方法；当为满阵时，方法称为全隐辛ＲＫ方法。对角隐式辛ＲＫ方法的级、阶不可能很高。迄今为止，已推导出的对”角隐式辛ＲＫ方法最高级阶为４级４阶。基于并行计算的目的，本文主要研究Ｓ级２阶的全隐辛ＲＫ方法。研究人员已经证明：级２阶的全隐辛ＲＫ方法是Ａ一稳定的【ｌ引。该方法的具体参数参见附录Ａ。多级高阶辛ＲＫ方法除具有辛算法本身的优点外【ｌ引，还具有内在的时间并行特性。如图１所示，１个Ｓ级的全隐辛ＲＫ方法，其一步的计算需要在Ｓ个时间点上同时进行求解，这相当于用变步长同时积分Ｓ＋１步。的表达式如下：２．１算法推导电力系统暂态稳定性数值计算可用下述模型表示：ｊ＝（）（３）～ｌ０＝ｇ（ｘ，＇，）∈’式中：为状态变量，ｘＲ；ｌ，为辅助变量，一∈般是网络节点电压，Ｒ（行为网络节点数）。定义∈∈ｙ，＝ｘ（ｔ＋Ｃ，ｈ）Ｒ。，ｆ（１，）∈∈ｚ＝ｖ（ｔ＋ｃ，ｈ）Ｒ。，ｆ（１，Ｓ）＝……（Ｊ，，Ｊ＇，，Ｊ，ＴＴ，Ｖｚ（ｚ，ｚ：Ｔ，，ｚ）…Ｆ（ｙ，ｚ）＝（ｆ（Ｙ。，），ｆ（ｙ，Ｚ２），，ｆ（Ｙ，ｚ））则用级的全隐辛ＲＫ方法对方程式（３）进行电力系统保护与控制求解的基本公式为：Ｘ＝ｅ０＋ｈ（Ａ０１）Ｆ（ｙ，ｚ）（４）∈０＝ｇ（，），（１，）（５）＝＋ｈ（ｂ０Ｉ）Ｖ（ｙ，ｚ）（６）上述方程中，ｅ是所有分量均为１的ｓ维矢量，Ｊ为ｍｘｍ维单位矩阵，０表示矩阵的Ｋｒｏｎｅｃｋｅｒ积（也称直积或张量积）。很显然，式（６）的计算首∈先需要同时在Ｓ个时间点上求解出Ｙｉ，ｚ，，ｆ（１，Ｓ）。将方程式（４）变换成（Ａ＠』）一Ｘ＝（＠Ｊ）。。０）＋ｈＦ（ｙ，ｚ）定义ｔ２＝Ａ～＝（ｑ）＝∑∈∈（ｇ）ｔａｘｉ，ｆ（１，）…伍＝（、戤，，，则上式可以进一步写成（ＱＩ）Ｘ＝＋ｈＦ（ｙ，ｚ）（７）定义仃＝（１，）ｆ＝ｌ＝ｚｆ’—∈：Ｑ＠Ｉｈ・ｄｉａｇ（ｏ－ｆ），ｆｆ，ｌ，Ｓ）（８）＝∈ｄｉａｇ（ｑｆ），＝ｄｉａｇ（ｒｆ），Ｄ４＝ｄｉａｇ（￣ｆ），ｆ（１，Ｓ）则联立方程（５）和（７）并利用牛顿法进行求解可以导出：［ＤＩＤ４ＬＡｙ—１ｌＪｌ（，ｚ）Ｊ……ＡＸ＝（，ｊ，，）ＴＡＶ＝（，，，）上式中，‘…ｚ）＝（，，，）∑∈，：＝，＋（Ｊ，，，ｚ，）一，＝（ｇＪ，），ｆ（１，），…∈ｚ）＝（Ｔ，ｓ：Ｔ，，￡），ｓ＝－ｇ（ｙ，ｚ），ｉ（１，Ｓ）定义∈Ｊｉ＝仃ｆ一ｆｆｆ，Ｐ＝ｄｉａｇ（Ｌ），ｆ（１，Ｓ）则从方程（９）可以推出：△（ＱＩ一Ｐ）＝Ｒ（ｙ，ｚ）＋ｈＤ（Ｊ，，ｚ）若下列条件成立ｐ［ｈ（ＱＪ）Ｐ］＝ｈ・［（ｏＩ）Ｐ］＜１（１０）式中］表示矩阵的谱半径，则利用矩阵求逆运算的松弛方法【３】可以导出：—（ＱＩｈｅ）～［Ｉ＋ｈ（Ａ０Ｊ）Ｐ］（Ｊ）在此情况下，最终可以导出：∑∈＋（口ｆ（１，）（１１）ｊ＝ｌ＝∑∈（），，＝，，ｆ（１，）（１２）ｊ＝ｌ△△∈：，：一，＝。ｇ（ｙ，ｚ），ｆ（１，Ｓ）（１３），＝一∈［ｇ（Ｊ，，ｚ）＋ｆ］，ｉ（１，）（１４）至此，本文导出了暂态稳定性的并行计算方法。在上述计算公式中，盯，，，ｆ，，相当于或就是传统逐步串行计算时的４个雅可比矩阵，它们均为稀疏矩阵；是网络节点导纳矩阵，通常情况下可以保持为常数矩阵。在上述算法的推导过程中，将方程式（４）转换为方程式（７）以及方程式（８）中所采用的矩阵分裂方法，是算法推导的２个关键步骤。从上述计算公式可以看出：在每个时问点上的求解过程中，只需对，进行稀疏三角分解。因此，该方法在每个时间点上进行求解时每次迭代的计算量不会大于传统的单步串行计算时每步迭代的计算量。上述算法具有很好的并行特性。不同时间点上“”的求解过程不存在相互等待的问题，只是在几个步骤上需要共享各自独立计算出的几个中间向量（，等的计算），而这可以通过共享存储方式或数据通信方式来实现。因此，Ｓ个时间点上的求解完全可以同步实现。很易理解，上述算法易于按时间和空间组织并行计算：不同时间点上的计算采用不同的处理器；同一时间点上的计算任务适合于矢量计算，亦即细粒度的空间并行计算。基于这一思路，可以采用ＧＰＵ或经重构后的ＦＰＧＡ来具体实现算法的并行装配。２．２算法收敛性分析数值积分方法的好坏，主要看算法的数值稳定性、计算精度和计算速度。如上所述，级２阶的全隐辛ＲＫ方法是稳定的。本文算法的收敛性，主要取决于方程式（１０）：若该方程成立，则本文导出的算法具有超线性收敛性，因为松弛牛顿法是超线性收敛的。为此，可以推导出：—汪芳宗，等基于多级高阶辛ＲｕｎｇｅＫｕｔｔａ方法的暂态稳定性并行计算方法－２５・≤≤（，）Ｐ］Ｉ刮・Ｉｌｅｌｌ＝ｌｌＡＩＩ・ｍａｘ（Ｉ］Ｊ￣１］）表１传统时间并行方法收敛性测试结果ｍａｘ（１ｌ小ＩＩ．１）—————Ｔａｂ．１Ｃｏｎｖｅｒｇｅｎｃｅｒａｔｅｏｆｃｌａｓｓｉｃｐａｒａ＿ｌｌｅｌ－ｉｎ－ｔｉｍｅｍｅｔｈｏｄ上式中，表示范数。对上式取列和范数（用．表示），并定义。＝，ｍａＸ（（ｌｌ＋ｌ・Ｉ。）＝≤则当ｈ，时，方程式（１ｏ）成立。，凡｜为常数矩阵（见附录Ａ），对不同级阶的辛ＲＫ方法，有（５１Ｏ．８７３１４。对具体的电力系统，是很容易计算出来的。因此，利用下式可以确定暂态稳定性并行计算的最大可用步长，即ｈ１．１４５３／７＇（１５）换言之，当步长满足上述条件时，本文所导出的算法是超线性收敛的。３仿真测试结果及分析并行算法的好坏，主要是看算法的并行度、加速比及并行效率。对时间并行计算而言，影响算法并行性能的主要因素是算法的收敛性：一般情况下，时间并行度愈高，算法的收敛性就会下降，因而很难获得高的加速比和并行效率【ｌ引。为评估本文所导出的算法的并行性能，利用ＩＥＥＥ１４５节点系统，对本文算法进行了仿真测试。该系统含５０台发电机。仿真测试中，为简单起见，发电机均采用经典模型；故障设定为在７号母线处发生三相短路，经０．１Ｓ后切除；暂态过程计算时程设定为１．５Ｓ。依据公式（１５）可以算出：利用本文方法进行计算时最大可用步长可达０．３４ｓ。也就是说，当步长ｈ０．３４时，本文算法将具备超线性收敛性。３．１算法收敛性测试为测试本文算法的收敛性，将本文方法与传统的时间并行牛顿计算方法进行了比较。所谓传统的“”时间并行牛顿计算方法，就是利用隐式梯形积分方法，以ｈ，＝０．０５ｓ的步长，一次性同时积分７７步（ｒ／即为时间并行度），而且使用严格的牛顿法进行“”整体迭代求解。时间并行牛顿计算方法具有最好的收敛性，但难以完全并行化Ｌ３，１引。表１是传统时间并行方法的收敛性测试结果。其中，收敛精度设定为￡＝１０（即牛顿残差向量范数小于该值时停止迭代），ｋ为牛顿求解过程中所需的最大迭代次数。“”参照上述情况，按等效时间并行度考虑，对本文方法的收敛性进行了测试。表２是在保持本“文算法的计算精度不低于传统的时间并行牛顿计”算方法的前提下，对本文算法的收敛性测试结果。表２本文方法收敛性测试结果Ｔａｂ．２Ｃｏｎｖｅｒｇｅｎｃｅｒａｔｅｏｆｔｈｅｐｒｏｐｏｓｅｄｐａｒａｌｌｅｌｍｅｔｈｏｄ４对比表ｌ和表２可以看出：在保持本文算法的“”计算精度不低于传统的时间并行牛顿计算方法的前提下，本文算法的收敛性与时间并行（严格）牛顿计算方法的收敛性大致相当。因此，本文算法很好地解决了时间并行度与收敛性之间的矛盾。３．２算法并行性能仿真测试表３是本文算法在采用不同级数即时间并行度时相对于经典的串行计算方法所获得的加速比和并行效率。该结果是在串行计算环境下进行模拟并行计算的测试结果。如２．１节所述，本文算法在Ｓ个时间点上的求解基本上是完全解耦即相互独立的，因此，在模拟并行计算中，可以测出每个时间点上求解过程所耗机时即ＣＰＵ时间，这个ＣＰＵ时间加上少量需集中处理的计算任务所耗机时，就是或相当于（真实）并行计算时间。经典的串行计算方法是采用常用的隐式梯形积分方法，并以ｈ＝０．０５ｓ的步长逐步进行积分求解，其总共所耗机时即串行ＣＰＵ时间为１２８８ｍｓ。表３本文算法并行性能仿真测试结果Ｔａｂ．３Ｓｉｍｕｌａｔｅｄｓｐｅｅｄｕｐｓａｎｄｅｆｆｉｃｉｅｎｃｉｅｓｏｆｔｈｅｐｒｏｐｏｓｅｄｐａｒａｌｌｅｌｍｅｔｈｏｄ步长ｈ３ｈＮ４ｈＮ５ｈＮ６ｈＮ（模拟）并行计算时间／ｍｓ加速比并行效率１４２．１１１８．４９８．９７８．７上述测试结果只是对本文算法的时间并行性能进行仿真测试的结果，算法的空间并行性能无法通过仿真测试来获得。对比相关或同类并行算法的性能，从表３可以看出：单以时间并行而言，本文算法具有很好的并行性能，可以达到较高的时间并行加速比并获得很高的并行效率。需要说明的是：通常情况下，并行算法的并行效率一般不会超过１００％。本文算法在Ｓ＝２，３时获得了超过常规的并行效率，这主要是得益于以下两点：．２６．电力系统保护与控制Ｓ级的算法具有Ｓ＋１个等效时间并行度；经典的串行计算方法尽管每步只需２次迭代，但在求解和迭代中需要对雅可比矩阵进行三角分解，而本文算法在每个时间点上的迭代求解中无需对雅可比矩阵进行三角分解（．是网络节点导纳矩阵，只需一次性三角分解）。４结论本文将多级高阶辛ＲＫ方法用于暂态稳定性计算，导出了一种新的暂态稳定性并行计算方法。理论分析及仿真测试结果表明，该算法具有以下优点：１）该算法具有内在的时间并行特性和超线性收敛性。２）多级高阶辛ＲＫ方法具有更高的计算精度。在保持相同的计算精度的前提下，本文算法的收敛性与时间并行（严格）牛顿计算方法的收敛性大致相当。因此，本文算法很好地解决了时间并行度与收敛性之间的矛盾，可以获得较高的时间并行加速比和很好的并行效率。有关本文算法的具体并行实现及验证结果将另文发表。附录Ａ２级４阶辛ＲＫ方法的Ｂｕｔｃｈｅｒ表为：０．２ｌｌ３２０．７８８６８０．２５００００．５３８６８０．５－０．０３８６８０．２５０００Ｏ．５３级６阶方法的Ｂｕｔｃｈｅｒ表为：０．１１２７０ｌ０．１３８８９－０．０３５９８０．００９７９—０．５１０．３００２６０．２２２２２０．０２２４９０．８８７３０｝０．２６７９９０．４８０４２０．１３８８９４级８阶方法的Ｂｕｔｃｈｅｒ表为：０．０６９４３０．３３００１０．６６９９９０．９３０５７０．０８６９６０．１８８１２Ｏ．１６７１９０．１７７４８０．１７３９３－０．０２６６００．１６３０４０．３５３９５０－３１３４５０．３２６０７０．Ｏ１２６３—０．０２７８８０．１６３０４０．３５２６８０．３２６０７５级１０阶方法的Ｂｕｔｃｈｅｒ表为：０．０４６９ｌ０．２３０７７０．５０．７６９２３０．０５９２３Ｏ．１２８１５０．１ｌ３７８Ｏ．１２ｌ２３０．１１６８８０．１１８４６－０．０１９５７Ｏ．１１９６６０．２６００００．２２９０００．２４４９１０．２３９３１０．０ｌｌ２５—０．０２４５９０．１４２２２０．３０９０４０．２７３ｌ９０．２８４４４－０．００３５６０．００６７４—０．０１４ｌ９０．０８６９６０．１７３９３－ｏ．００５５９０．０１０３２—０．０２０６９０．１】９６６Ｏ．ＯＯ１５９—０．００２７７０．００４６８—０．００９６９０．２３９３１０．１１８４６参考文献［１］赵建成，谢小荣，穆钢．基于ＷＡＭＳ的电力系统暂态稳定的快速预测［Ｊ］．继电器，２００５，３３（７）：１．５．—ＺＨＡＯＪｉａｎｇｃｈｅｎｇ，ＸＩＥＸｉａｏ－ｒｕｎｇ，ＭＵＧａｎｇ．—ＷＡＭＳｂａｓｅｄｒｅａｌ－ｔｉｍｅｐｒｅｄｉｃｔｉｏｎｏｆｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙ—ｆｏｒｕｌｔｉｍａｃｈｉｎｅｐｏｗｅｒｓｙｓｔｅｍｓ［Ｊ］．Ｒｅｌａｙ，２００５，３３（７）：１．５．［２］石恒初．基于Ｐｃ机群的电力系统暂态稳定评估［Ｊ］．电力系统保护与控制，２００９，３７（１０）：５－９，１４．—ＳＨＩＨｅｎｇｃｈｕ．Ｐｅｒｓｏｎａｌｃｏｍｐｕｔｅｒｃｌｕｓｔｅｒｂａｓｅｄｐｏｗｅｒｓｙｓｔｅｍｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙａｓｓｅｓｓｍｅｎｔ［Ｊ］．ＰｏｗｅｒＳｙｓｔｅｍＰｒｏｔｅｃｔｉｏｎａｎｄＣｏｎｔｒｏｌ，２００９，３７（１０）：５－９，１４．［３］汪芳宗．基于高度并行松驰牛顿方法的暂态稳定性实时分析计算的并行算法［Ｊ］．中国电机工程学报，１９９９，１９（１１、：１４－１８．ＷＡＮＧＦａｎｇ－ｚｏｎｇ．ＰａｒａｌｌｅｌａｌｇｏｒｉｔｈｍｏｆｈｉｇｈｌｙｐａｒａｌｌｅｌｒｅｌａｘｅｄＮｅｗｔｏｎｍｅｔｈｏｄｆｏｒｒｅａ１．ｔｉｍｅｓｉｍｕｌａｔｉｏｎｏｆｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉ州Ｊ］．ＰｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＣＳＥＥ，１９９９，１９（１１）：１４・１８．［４］ＦｅｎｇＫ．Ｉ）ｉｆｆｅｒｅｎｃｅｓｃｈｅｍｅｓｆｏｒＨａｍｉｌｔｏｎｉａｎｆｏｒｍａｌｉｓｍａｎｄｓｙｍｐｌｅｃｔｉｃｇｅｏｍｅｔｒｙ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌ—Ｍａｔｈｅｍａｔｉｃｓ，１９８６，４（３）：２７９２８９．［５］ＦｅｎｇＫ，ＱｉｎＭＺ．ＨａｍｉｌｔｏｎｉａｎａｌｇｏｒｉｔｈｍｓｆｏｒＨａｍｉｌｔｏｎｉａｎｄｙｎａｍｉｃａｌｓｙｓｔｅｍｓ［Ｊ］．ＰｒｏｇｒｅｓｓｉｎＮａｔｕｒａｌ—Ｓｃｉｅｎｃｅ，１９９１，１（２）：１０５１１６．［６］冯康，秦孟兆．哈密尔顿系统的辛几何算法［Ｍ】．杭州：浙江科学技术出版社，２００３．ＦＥＮＧＫａｎｇ，ＱＩＮＭｅｎｇ－ｚｈａｏ．Ｓｙｍｐｌｅｃｔｉｃｇｅｏｍｅｔｒｉｃａｌｇｏｒｉｔｈｍｆｏｒｈａｍｉｌｔｏｎｉａｎｓｙｓｔｅｍｓ［Ｍ］．Ｈａｎｇｚｈｏｕ：ＺｈｅｊｉａｎｇＳｃｉｅｎｃｅ＆ＴｅｃｈｎｏｌｏｇｙＰｒｅｓｓ．２００３．［７］肖爱国，李寿佛．辛Ｒｕｎｇｅ．Ｋｕｔｔａ方法的特征与构造［Ｊ］Ｊ．高等学校计算数学学报，１９９５，１７（３）：２１３．２２２．—ＸＩＡＯＡｉ－ｇｕｏ，ＬＩＳｈｏｕｆｕ．Ｃｈａｒａｃｔｅｒｉｚａｔｉｏｎａｎｄｃｏｎｓｔｒｕｃｔｉｏｎｏｆｓｙｍｐｌｅｃｔｉｃｒｕｎｇｅ－ｋｕｔｔａｍｅｔｈｏｄｓ［Ｊ］．ＮｕｍｅｒｉｃａｌＭａｔｈｅｍａｔｉｃｓ：ＡＪｏｕｒｎａｌｏｆＣｈｉｎｅｓｅ—Ｕｎｉｖｅｒｓｉｔｉｅｓ，１９９５，１７（３）：２１３２２２．［８］ＳｕｎＧ．ＡｓｉｍｐｌｅｗａｙｃｏｎｓｔｒｕｃｔｉｎｇｓｙｍｐｌｅｃｔｉｃＲｕｎｇｅ－Ｋｕｔｔａｍｅｔｈｏｄｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＭａｔｈｅｍａｔｉｃｓ，２０００，１８（１）：６１－６８．［９］ＳｕｎＧ．ＣｏｎｓｔｒｕｃｔｉｏｎｏｆｈｉｇｈｏｒｄｅｒｓｙｍｐｌｅｃｔｉｃＲｕｎｇｅ－ＫｕＲａｍｅｔｈｏｄｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＭａｔｈｅｍａｔｉｃｓ，１９９３，１１（３）：２５０－２６０．［１０］Ｓａｎｚ－ＳｅｍａＪＭ．Ｒａｎｇｅ－ＫｕＲａｓｃｈｅｍｅｓｆｏｒＨａｍｉｌｔｏｎｉａｎｓｙｓｔｅｍｓ［Ｊ］．ＢＩＴＮｕｍｅｒｉｃａｌＭａｔｈｅｍａｔｉｃｓ，１９８８，２８（４）：８７７．８８３．［１１］蒋长锦．四级四阶对角隐式辛Ｒｕｎｇｅ．ＫｕＲａ方法参数计算［Ｊ］．数值计算与计算机应用，２００２，２５（３）：１６１－１６６．ＪＩＡＮＧＣｈａｎｇ－ｊｉｎ．Ｏｎｃｏｍｐｕｔｅｏｆｐａｒａｍｅｔｅｒｓｆｏｒ４－ｓｔａｇｅ—４－ｏｒｄｅｒｄｉａｇｏｎａｌｌｙｉｍｐｌｉｃｉｔｓｙｍｐｌｅｃｔｉｃＲｕｎｇｅＫｕｔｔａｍｅｔｈｏｄｓ［Ｊ］．ＪｏｕｒｎａｌｏｎＮｕｍｅｒｉｃａｌＭｅｔｈｏｄｓａｎｄＣｏｍｐｕｔｅｒＡｐｐｌｉｃａｔｉｏｎｓ，２００２，２５（３）：１６１－１６６．［１２］李寿佛．刚性微分方程算法理论【Ｍ】．长沙：湖南科学技术出版社，１９９７．—ＬＩＳｈｏｕｆｕ．Ｔｈｅｏｒｙｏｆｃｏｍｐｕｔａｔｉｏｎａｌｍｅｔｈｏｄｓｆｏｒｓｔｉｆｆｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎｓ［Ｍ］．Ｃｈａｎｇｓｈａ：ＨｕｎａｎＳｃｉｅｎｃｅ＆ＴｅｃｈｎｏｌｏｇｙＰｒｅｓｓ，１９９７．［１３］汪芳宗，何一帆．电力系统暂态稳定性数值计算的几种新方法及其比较［Ｊ］．电力系统保护与控制，２００９，３７（２３）：１５．１９．（下转第３２页ｃｏｎｔｉｎｕｅｄｏｎｐａｇｅ３２）－３２一电力系统保护与控制（上接第２１页ｃｏｎｔｉｎｕｅｄｆｒｏｍｐａｇｅ２１）［１４］赵晋泉，江晓东，张伯明．一种用于静态稳定分析的故障参数化连续潮流模型［Ｊ］．电力系统自动化，２００４，２８（１４）：４５－４９．———ＺＨＡＯＪｉｎｑｕａｎ，ＪＩＡＮＧＸｉａｏｄｏｎｇ，ＺＨＡＮＧＢｏｍｉｎｇ．Ａｎｅｗｃｏｎｔｉｎｇｅｎｃｙｐａｒａｍｅｔｅｒｉｚａｔｉｏｎｃｏｎｔｉｎｕａｔｉｏｎｐｏｗｅｒｆｌｏｗｍｏｄｅｌｆｏｒｓｔｅａｄｙｓｔａｂｉｌｉｔｙａｎａｌｙｓｉｓ［Ｊ］．ＡｕｔｏｍａｔｉｏｎｏｆＥｌｅｃｔｒｉｃＰｏｗｅｒＳｙｓｔｅｍｓ，２００４，２８（１４）：４５－４９．［１５］ＤｏｂｓｏｎＩ，ＬｕＬＭ．Ｖｏｌｔａｇｅｃｏｌｌａｐｓｅｐｒｅｃｉｐｉｔａｔｅｄｂｙｔｈｅｉｍｍｅｄｉａｔｅｃｈａｎｇｅｉｎｓｔａｂｉｌｉｔｙｗｈｅｎｇｅｎｅｒａｔｏｒｒｅａｃｔｉｖｅｐｏｗｅｒｌｉｍｉｔｓａｒｅｅｎｃｏｕｎｔｅｒｅｄ［Ｊ］．ＩＥＥＥＴｒａｎｓｏｎＣｉｒｃｕｉｔｓ—ａｎｄＳｙｓｔｅｍｓＩ：ＦｕｎｄａｍｅｎｔａｌＴｈｅｏｒｙａｎｄＡｐｐｌｉｃａｔｉｏｎｓ．１９９２，３９（９）：７６２－７６６．［１６］赵晋泉，江晓东，张伯明．一种静态电压稳定临界点的识别和计算方法［Ｊ］．电力系统自动化，２００４，２８（２３）：２８。３２．ＺＨＡＯＪｉｎ－ｑｕａｎ，ＪＩＡＮＧＸｉａｏ－ｄｏｎｇ，ＺＨＡＮＧＢｏ－ｍｉｎｇ．Ａｎａｐｐｒｏａｃｈｆｏｒｉｄｅｎｔｉｆｉｃａｔｉｏｎａｎｄｃｏｍｐｕｔａｔｉｏｎｏｆｓｔａｔｉｃｖｏｌｔａｇｅｓｔａｂｉｌｉｔｙｃｒｉｔｉｃａｌｐｏｉｎｔ［Ｊ］．ＡｕｔｏｍａｔｉｏｎｏｆＥｌｅｃｔｒｉｃＰｏｗｅｒＳｙｓｔｅｍｓ，２００４，２８（２３）：２８－３２．［１７］ＣｈｉａｎｇＨＤ，ＦｌｕｅｃｋＡＪ，ＳｈａｈＫＳ，ｅｔａ１．ＣＰＦＬＯＷ：Ａ—ｐｒａｃｔｉｃａｌｔｏｏｌｆｏｒｔｒａｃｉｎｇｐｏｗｅｒｓｙｓｔｅｍｓｔｅａｄｙｓｔａｔｅｓｔａｔｉｏｎａｒｙｂｅｈａｖｉｏｒｄｕｅｔｏｌｏａｄａｎｄｇｅｎｅｒａｔｉｏｎｖａｒｉａｔｉｏｎｓ［Ｊ］．ＩＥＥＥＴｒａｎｓｏｎＰｏｗｅｒＳｙｓｔｅｍｓ，１９９５，１０（２）６２３．６３４．收稿日期：２０１０－０６－０７；修回Ｅｔ期：２０１０－１１－２４作者简介：赵晋泉（１９７２－），男，博士，教授，博士生导师，主要从事电力系统优化运行、静态稳定评估与控制和电力市场定价理论方面的研究工作。Ｅ－ｍａｉｌ：ｊｑｚｈａｏｚ＠ｔｏｍ．ｃｏｍ（上接第２６页ｃｏｎｔｉｎｕｅｄｆｒｏｍｐａｇｅ２６）—ＷＡＮＧＦａｎｇｚｏｎｇ，ＨＥＹｉ－ｆａｎ．Ｓｅｖｅｒａｌｎｅｗｎｕｎａｅｒｉｃａｌｍｅｔｈｏｄｓａｎｄｔｈｅｉｒｃｏｍｐａｒａｔｉｖｅｓｔｕｄｉｅｓｆｏｒｐｏｗｅｒｓｙｓｔｅｍｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙａｎａｌｙｓｉｓ［Ｊ］．ＰｏｗｅｒＳｙｓｔｅｍＰｒｏｔｅｃｔｉｏｎａｎｄＣｏｎｔｒｏｌ，２００９，３７（２３）：１５－１９．［１４］汪芳宗，杨力森，高学军，等．基于松弛牛顿法的互联电网并行分布式潮流计算方法［Ｊ］．继电器，２００７，３５（１６）：１８－２２．—ＷＡＮＧＦａｎｇ－ｚｏｎｇ，ＹＡＮＧＬｉｓｅｎ，ＧＡＯＸｕｅ－ｊｕｎ，ｅｔａ１．ＰａｒａｌｌｅｌａｎｄｄｉｓｔｒｉｂｕｔｅｄｌｏａｄｆｌｏｗｓｏｌｕｔｉｏｎｏｆｉｎｔｅｒｃｏｎｎｅｃｔｅｄｐｏｗｅｒｓｙｓｔｅｍｓｕｓｉｎｇｒｅｌａｘｅｄＮｅｗｔｏｎｍｅｔｈｏｄ［Ｊ］．Ｒｅｌａｙ，２００７，３５（１６）：１８－２２．［１５］汪芳宗．基于高度并行松驰牛顿方法的暂态稳定性实时分析计算方法的并行装配［Ｊ］．中国电机工程学报，１９９９，１９（１１）：１８－２１．ＷＡＮＧＦａｎｇ－ｚｏｎｇ．ＰａｒａｌｌｅｌｉｍｐｌｅｍｅｎｔａｔｉｏｎｏｆｈｉｇｈｌｙｐａｒａｌｌｅｌｒｅｌａｘｅｄＮｅｗｔｏｎｍｅｔｈｏｄｆｏｒｒｅａｌ－ｔｉｍｅｓｉｍｕｌａｔｉｏｎｏｆｔｒａｎｓｉｅｎｔｓｔａｂｉｌｉｔｙ［Ｊ］．ＰｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＣＳＥＥ，１９９９，１９（１１、：１８－２１．收稿日期：２０１０＂０６－１８；修回日期：２０１０－０８－０２作者简介：汪芳宗（１９６６－），男，博士，教授，主要从事电力系统—自动化、电工新技术及应用的研究工作；Ｅｍａｉｌ：ｆｚｗａｎｇ＠ｃｔｇｕ．ｅｄｕ．ｃｎ何一帆（１９８４－），男，硕士，主要从事电力系统分析计算与控制的研究工作。