大规模电力系统小干扰稳定性的特征值分析系统

摘要

一种电力系统仿真及分析技术领域的大规模电力系统小干扰稳定性的特征值分析系统，包括：用于实现BPA数据文件的无缝读写操作的BPA数据接口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵计算模块、潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算引擎模块、部分特征值计算引擎模块和振荡模式提取和分析模块，采用TTQRE、多种计算方案的IRAM及JDM，实现对实际大规模电网进行全方位、无遗漏的小干扰稳定性特征值仿真分析。

说明书

技术领域

本发明涉及的是一种电力系统仿真及分析技术领域的系统，具体是一种采用三种不同特征值算法的大规模电力系统小干扰稳定性的特征值分析系统。 

背景技术

电力系统小干扰稳定性指的是当电网经历微小扰动后，继续保持同步运行的能力。电力系统小干扰稳定性一般采用李雅普诺夫第一法作为判定标准。李雅普诺夫第一法指出，如果系统线性化后的状态矩阵没有出现零或者正实部的特征值，那么就可以判定当前系统是小干扰稳定的。因此，在电力系统中，特征值计算一直以来都是实现低频功率振荡模式的识别、各类稳定控制器的装配选址和参数优化、运行参数对控制参数的灵敏度分析、在线检测振荡数据的模态信息提取等功能的重要前提和基础保障。 

由于电网规模的不断扩大，自上个世纪80年代开始，电力行业的研究人员就已投入大量的精力到特征值计算方法和小干扰稳定性分析系统的研究与开发中。至今为止，世界范围内已出现了诸多用于大规模电力系统小干扰稳定性的特征值分析系统，例如：美国太平洋瓦斯与电力公司开发的EISEMAN、美国电科院联合加拿大安大略水电局共同开发的SSSP、加拿大动力技术实验室开发的SSAT、德国西门子公司开发的NEVA、巴西电力科学研究中心开发的PacDyn、中国电力科学研究院开发的PSASP和PSD‐SSAP等。 

按照待求特征值的数量进行归类，电力系统小干扰稳定性的特征值分析方法分为全部特征值分析法和部分特征值分析法。虽然几乎所有的大规模电力系统小干扰稳定性分析系统均包含全部特征值分析和部分特征值分析两个方面，但仍然存在以下的问题和不足： 

1)全部特征值分析法的核心模块仍采用的是上个世纪60年代Kublanovskay和Francis提出的双重步位移隐式QR算法，从而导致国内电力行业人士普遍认为：对于大型电力系统的小干扰稳定性分析而言，QR算法存在内存不足、计算时间很长、计算出的特征值误差很大、算法可能不收敛等问题[王康，金宇清，甘德强，等.电力系统小信号稳定分析与控制综述[J].电力自动化设备，2009(5):10‐19.薛禹胜，郝思鹏，刘俊勇.关于低频振荡分析方法的评述[J].电力系统自动化，2009，33(3):1‐8.中国电力科学研究院，PSASP7.0版小干扰计算用户手册[R]，北京:中国电力科学研究院，2010.中国电力科学研究院，PSD‐SSAP小干扰稳定性分析程序用户手册(2.5.2版)[R]，北京:中国电力科学研究院，2012.]。不可否认，早期的双重步位移隐式QR算法的确存在以上所述的种种问题。然而，随着数值计算方法以及计 算机软硬件技术的不断进步，QR算法早已完成了从双重步位移—大块多重步位移—链式小块多重步位移—带有积极早期收缩策略的两步小块多重步位移的演变[Francis J G F.The QR transformation a unitary analogue to the LR transformation—Part 1[J].The Computer Journal，1961，4(3):265‐271.Francis J G F.The QR transformation—part 2[J].The Computer Journal，1962，4(4):332‐345.Bai Z，Demmel J.On a Block Implementation of Hessenberg Multishift QR Iteration[J].International Journal of High Speed Computing，1989，1(1):97‐112.Braman K，Byers R，Mathias R.The multishift QR algorithm.part I:Maintaining well‐focused shifts and level 3 performance[J].SIAM Journal on Matrix Analysis and Applications，2002，23(4):929‐947.Braman K，Byers R，Mathias R.The multishift QR algorithm.Part II:Aggressive early deflation[J].SIAM Journal on Matrix Analysis and Applications，2002，23(4):948‐973.]，个人计算机也早已从仅能支持4GB寻址的32位，变成了能够支持16EB(1EB＝2³⁰GB)寻址的64位。因此，基于双重步位移隐式QR算法和32位个人计算机的全部特征值分析法，显然已经不能满足当前大规模电力系统小干扰稳定性分析的计算需求； 

2)部分特征值分析法是目前分析大规模电力系统小干扰稳定性的主流方法。部分特征值分析法的核心模块是一类迭代投影法(Iterative projection methods)[Bai Z J，et al.Templates for the solution of algebraic eigenvalue problems:a practical guide[M].Siam，2000.]。在近10年来的大规模电力系统小干扰稳定性分析的特征值算法研究中，出现频率最多的迭代投影法有两个，一个是Krylov子空间下的IRAM(Implicitly Restarted Arnoldi Method，隐式重启动Arnoldi算法)[Kim D J，Moon Y H.Application of the implicit restarted Arnoldi method to the small‐signal stability of power systems[J].Journal of Electrical Engineering&Technology，2007，2(4):428‐433.仲悟之，宋新立，汤涌，等.基于多进程的电力系统频域特征值并行搜索算法[J].电力系统自动化，2010(21):11‐16.]，另一个则是非Krylov子空间下的Jacobi‐Davidson方法(JDM)[杜正春，刘伟，方万良，等.基于Jacobi‐Davidson方法的小干扰稳定性分析中关键特征值计算[J].中国电机工程学报，2005，25(14):19‐24.Tsai S H，Lee C Y，Wu Y K.Efficient calculation of critical eigenvalues in large power systems using the real variant of the Jacobi–Davidson QR method[J].IET generation，transmission&distribution，2010，4(4):467‐478.]。 

虽然IRAM已被广泛集成于各电力系统小干扰稳定性分析系统中，但其算法本身具备的功能却没能得到最大程度的发挥，例如：PSD‐SSAP和PSASP只提供了基于平移‐逆变换的扫频IRAM，而非常适合于计算关键振荡模式的Cayley变换IRAM，在两种方法中均没有体现。扫频IRAM虽然对于极特征值有着线性的收敛速度，但其搜索行为存在着随机性，在参数配置不合 理的情况下，容易出现“漏根”，即漏掉一些特别关键的振荡模式。当精确求解修正方程时，JDM具有渐进二阶的收敛速度，但更重要的是，即使不使用谱变换，JDM仍然能够按照某种指定的规则(实部最大、阻尼比最小)，顺序收敛出影响电力系统小干扰稳定性的那些关键振荡模式。 

然而，不像IRAM具有开源的数学库ARPACK，由于JDM仅有Matlab下的示例程序，因此，JDM一直常见于学术研究，而没有集成到任何的大规模电力系统小干扰稳定性特征值分析系统中。 

“三华”(华中‐华东‐华北)特高压交/直流电网的互联，一定会出现新增的低频功率振荡模式，而如今处于建设初期的“强直弱交”电网，极有可能使得其中的某些模式将要或者已经呈现出不稳定的状态，因此，采用三种数值原理完全不同的特征值算法来分析大规模电力系统的小干扰稳定性，不仅具有理论价值，更具有现实意义。 

发明内容

本发明针对现有技术存在的上述不足，提出一种大规模电力系统小干扰稳定性的特征值分析系统，根据PSD‐BPA格式的电网数据文件(*.dat和*.swi)，采用TTQRE(Two‐tone small‐bulge multishift QR algorithm with aggressive early deflation，带有积极早期收缩策略的两步小块多重步位移QR算法)，快速、准确的识别出影响系统动态稳定的全部机电振荡模式；采用具有多种计算方案的IRAM，提取和分析影响系统动态稳定的部分关键机电振荡模式；采用具有多种计算方案的JDM，提取和分析影响系统动态稳定的部分关键机电振荡模式，通过匹配和比较三种算法得到的机电振荡模式，最终实现对实际大规模电网进行全方位、无遗漏的小干扰动态稳定的振荡模式识别和模态分析。 

本发明是通过以下技术方案实现的： 

本发明涉及一种大规模电力系统小干扰稳定性的特征值分析系统，包括：BPA数据接口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵计算模块、潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算引擎模块、部分特征值计算引擎模块和振荡模式提取和分析模块，其中：BPA数据接口模块接收电网数据文件，并分别与潮流计算模块和系统建模及其线性化模块相连并提供API(Application Program Interface，应用程序接口)；稀疏矩阵计算模块分别与潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块和部分特征值计算引擎模块相连并提供API；潮流计算模块根据BPA数据接口模块以及稀疏矩阵计算模块提供的API，实现大规模电力系统的潮流计算并输出潮流计算结果至系统建模及其线性化模块；系统建模及其线性化模块根据BPA数据接口模块提供的API和潮流计算模块提供的潮流计算结果，实现大规模电力系统的线性化建模，即生成增广的系统状态矩阵并分别输出至状态矩阵计算模块和部分特征值计算引擎模块；状态矩阵计算模块根据稀疏矩阵计 算模块提供的API和系统建模及其线性化模块提供的增广的系统状态矩阵，获得系统状态矩阵并输出至全部特征值计算引擎模块；全部特征值计算引擎模块根据TTQRE从系统状态矩阵中获得系统状态矩阵的全部特征值以及部分左/右特征向量并输出至振荡模式提取和分析模块；部分特征值计算引擎模块根据稀疏矩阵计算模块提供的API，以RCI(Reverse Communication Interface，逆通信接口)的方式执行IRAM或者JDM的数值迭代过程，获得系统状态矩阵的部分特征值和部分左/右特征向量并输出至振荡模式提取和分析模块；振荡模式提取和分析模块根据用户指定合并来自全部特征值计算引擎模块或者部分特征值计算引擎模块的系统状态矩阵的特征值和左/右特征向量，实现机电振荡模式的识别和模态分析，并将所得的模式信息和模态分析结果以Excel表格的方式输出给用户。 

所述的增广的系统状态矩阵J_aug满足：$J_{\mathrm{aug}}=\left(\begin{array}{cc}{J}_A& J_B\\ J_C& J_D\end{array}\right),$其中：J_A为动态元件线性化模型拼接而成的分块对角矩阵，J_B为动态元件与静态元件的联系矩阵，J_C为静态元件与动态元件的联系矩阵，J_D为修正的系统导纳矩阵。 

所述的系统状态矩阵S是指：

所述的全部特征值计算引擎模块内包括TTQRE单元，该单元与状态矩阵计算模块相连，并根据系统状态矩阵获得系统状态矩阵的全部特征值和部分左/右特征向量后，输出至振荡模式提取和分析模块。 

所述的部分特征值计算引擎模块包括：IRAM单元和JDM单元，其中：IRAM单元与系统建模及其线性化模块相连，根据用户的选择来确定采用平移‐逆变换结合模最大Ritz值选择策略还是Cayley变换结合模最大Ritz值选择策略的计算方案，并根据增广的系统状态矩阵获得系统状态矩阵的部分特征值和部分左/右特征向量后，输出至振荡模式提取和分析模块；JDM单元与系统建模及其线性化模块相连，根据用户的选择来确定采用平移‐逆变换结合模最大Ritz值选择策略、Cayley变换结合模最大Ritz值选择策略、原始矩阵结合实部最大Ritz值选择策略还是原始矩阵结合阻尼比最小Ritz值选择策略的计算方案，并根据增广的系统状态矩阵获得系统状态矩阵的部分特征值和部分左/右特征向量后，输出至振荡模式提取和分析模块。 

所述的振荡模式识别和模态分析是指：根据每一对左/右特征向量获得相应特征值的相关因子，根据相关因子获得该特征值的机电回路相关比，若机电回路相关比大于1且该特征值为复数，则判定该特征值为机电振荡模式，相应的右特征向量为振荡模态；相关因子模值的大小反应了电力系统动态元件与振荡模式之间的耦合程度，在耦合程度最大的机组上加装电力系统稳定器，是抑制由该振荡模式诱发的电力系统低频功率振荡的一贯做法；振荡模态的幅值决定了该振荡模式下机组的摇摆程度，角度决定了该振荡模式下机组与机组之间的相对摇摆关系， 振荡模态是判断振荡模式对机组的影响程度、判断多个机组是否可以进行动态等值的重要指标。 

技术效果 

与现有技术手段相比，本发明的有益效果包括： 

1.采用了能够代表当前技术水平的TTQRE，来实现大规模电力系统小干扰稳定性的全部特征值分析。该算法在特征值迭代求解环节，采用链式小块追赶技术，消除了多重步位移QR块追赶的位移模糊现象，解决了算法不收敛的问题；采用积极早期收缩策略，减少了QR迭代的次数，提高了算法的计算精度；采用矩阵‐矩阵的高级运算，提高了CPU缓存的利用率，节省了算法的计算时间。从工程实用的角度出发，论证了QR算法对于大型电力系统小干扰稳定性全部特征值分析的可行性； 

2.将IRAM的谱变换统一描述为算子形式，实现了具有平移‐逆变换结合模最大Ritz值选择策略、Cayley变换结合模最大Ritz值选择策略两种计算方案的IRAM。Cayley变换IRAM的最大优点在于，能够有效识别电力系统中阻尼比小于等于某个指定值的一些关键振荡模式，相比平移‐逆变换IRAM，Cayley变换IRAM的搜索行为基本不存在随机性，因此，即便是在参数配置不完全合理的情况下，Cayley变换IRAM仍可以识别出绝大部分的关键振荡模式，很大程度上降低了用IRAM分析大规模电力系统的小干扰稳定性时，出现“漏根”的可能性； 

3.在国内外首次以RCI的方式实现了JDM。RCI不仅完美诠释了特征值算法与电力系统应用之间的分界点，RCI还使得迭代投影法具备了一致的调用规则。相比Krylov子空间下的迭代投影法，非Krylov子空间下的JDM不仅具有更快的收敛速度，而且在不使用谱变换的情况下，JDM仍可以按照某种指定的Ritz值选择策略(阻尼比最小、实部最大)，依次收敛出想要的特征值。实现了具有4种计算方案的JDM，并将其与IRAM一起集成到了大规模电力系统小干扰稳定性特征值分析系统中，对JDM的工程实用化起到了推动性的作用。 

附图说明

图1为本发明提供的技术方案实施图； 

图2a和图2b分别为IRAM和JDM的数值迭代过程示意图的左半部分和右半部分； 

图3为平移‐逆变换扫频IRAM的应用效果图； 

图4为Cayley变换结合模最大Ritz值选择策略IRAM的应用效果图； 

图5为原始矩阵结合阻尼比最小Ritz值选择策略JDM的应用效果图； 

图6为原始矩阵结合阻尼比最小Ritz值选择策略JDM的收敛过程图； 

图7为HD8241系统中某一区间振荡模式的模态分析图。 

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施， 给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。 

实施例1 

如图1所示，本实施例包括：用于实现BPA数据文件的无缝读写操作的BPA数据接口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵计算模块、潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算引擎模块、部分特征值计算引擎模块和振荡模式提取和分析模块，其中： 

所述的BPA数据接口模块：由于国内的电网企业偏向于使用BPA的文件格式来存放实际电网数据，为此，基于面向对象编程思想，采用C++开发了具有动态链接库形式的BPA数据接口模块。该模块不仅具有BPA数据文件的无缝读写功能，还可以为其他模块提供API。该模块的具体实施方式包括以下3个主要步骤： 

<1>BPA数据文件的层次结构分析，用来确定电网的静态/动态元件参数在BPA数据文件中的存放位置和存放顺序，实现数据卡片和控制卡片的有效分离； 

<2>类结构框架的设计，按照面向对象编程的要求，设计一套与BPA程序用户手册具有相同层次关系的类结构框架，方便接口模块的维护与扩展； 

<3>BPA数据接口的生成，采用C++对所设计的类结构框架进行编程实现，生成的数据接口不仅可以提供API，更可以提供独立的应用程序。 

所述的稀疏矩阵计算模块：任何大规模电力系统的仿真分析与计算，都离不开稀疏技术。该模块采用了佛罗里达州立大学Timothy A.Davis教授提供的SuiteSparse，来实现与稀疏矩阵相关的所有运算，以动态链接库的形式为其他模块提供API。该模块主要由5部分构成： 

<1>数据输入/输出，实现稀疏矩阵从硬盘到内存之间的双向数据传输，该部分既支持2进制格式保存的稀疏矩阵，同时又满足32/64位寻址的需求； 

<2>基本矩阵运算，实现稀疏矩阵的一些常规计算，例如：转置、重排、矩阵‐向量乘法、矩阵‐矩阵加法、矩阵‐矩阵乘法等； 

<3>稀疏LU分解，实现稀疏矩阵的符号分析和LU数值分解； 

<4>稀疏前代/回代，实现稀疏矩阵的快速前代/回代； 

<5>扩展矩阵运算，实现一些扩展的稀疏矩阵运算，例如：计算稀疏矩阵关于某个子块的Schur补、截断的稀疏矩阵‐向量乘法等。 

现有技术中SuiteSparse是与SuperLU、PARDISO、TAUCS等齐名的稀疏矩阵解算器，更重要的是，在SuiteSparse源码的基础上，能够开发出专用于电力系统小干扰稳定性特征值分析的稀疏矩阵计算功能，例如，稀疏矩阵计算模块的第<4>和第<5>部分所示的功能。这些扩展稀疏矩阵计算功能的实际应用效果，将在后面给出。 

所述的潮流计算模块：利用BPA数据接口模块以及稀疏矩阵计算模块提供的API，实现大规模电力系统的潮流计算，为系统建模及其线性化模块提供潮流计算结果。 

所述的系统建模及其线性化模块根据BPA数据接口模块提供的API和潮流计算模块提供的潮流计算结果，并结合文献[励刚，陈陈.应用设计模式开发小干扰稳定性分析方法[J].中国电机工程学报，2002，22(1):12‐16.]提出的统一元件联接建模方法，来对小干扰稳定性分析进行系统线性化建模，同时生成增广的系统状态矩阵J_aug满足：$J_{\mathrm{aug}}=\left(\begin{array}{cc}{J}_A& J_B\\ J_C& J_D\end{array}\right)$其中：J_A为动态元件线性化模型拼接而成的分块对角矩阵，J_B为动态元件与静态元件的联系矩阵，J_C为静态元件与动态元件的联系矩阵，J_D为修正的系统导纳矩阵。 

所述的状态矩阵计算模块：利用稀疏矩阵计算模块提供的API和系统建模及其线性化模块提供的稀疏矩阵J_A、J_B、J_C和J_D，实现系统状态矩阵S的稀疏求解：由于S是稠密矩阵，而J_A、J_B、J_C、J_D是稀疏矩阵，根据稠密矩阵的计算方法来求解S，一方面计算时间会很长，另一方面内存开销也会很大，于是，该模块采用以下步骤求解S： 

<1>对J_D进行符号分析以及稀疏LU分解：P_DJ_DQ_D＝L_DU_D；其中：P_D和Q_D分别为实现列选主元和近似最小度排序的重排矩阵。 

<2>用J_A来初始化稠密存储的S：S＝full(J_A) 

<3>利用稀疏矩阵计算模块提供的扩展矩阵运算功能，计算稀疏矩阵J_aug关于子块J_D的Schur补：由于电力系统小干扰稳定性分析中的J_B和J_C具有非常特殊的稀疏结构，而使得第<3>步公式的右端项并不需要执行全部的数值计算，因此，稀疏矩阵计算模块才会在SuiteSparse的源码基础上，实现了“计算稀疏矩阵关于某个子块的Schur补”的扩展矩阵运算功能。该功能的实际应用效果将在后面给出。 

所述的全部特征值计算引擎模块：采用TTQRE，实现对状态矩阵S的全部特征值、部分左/右特征向量的求解，包括以下6个主要步骤： 

<1>矩阵平衡，采用矩阵重排分离出已经暴露于S对角线上的部分实特征值，采用循环缩放减小平衡后矩阵A的2范数：A＝(PD)^-1S(PD)；其中：P为重排矩阵，D为对角矩阵。 

<2>Hessenberg约化，采用一系列Householder反射，将A约化为上Hessenberg形式H：H＝Q^TAQ；其中：Q为一系列Householder反射累积而成的正交矩阵。 

<3>链式小块追赶的多重步位移隐式QR迭代，采用一系列Householder反射和Givens旋转执行链式QR块追赶，采用收缩窗收敛多个特征值，将H相似变换为拟上三角矩阵T，S 的全部特征值计算结果即为T的所有对角元：T＝Q^TAQ；其中：Q为进一步累积而成的正交矩阵，包含一系列Householder反射和Givens旋转。 

<4>标准Schur型部分特征向量的求解，确定T中想要的特征值Λ后，计算T关于Λ的左特征向量V和右特征向量U：$\left(\begin{array}{c}{V}^{H}T={\mathrm{\Lambda V}}^H\\ \mathrm{TU}=\mathrm{U\Lambda }\end{array}\right);$

<5>状态矩阵部分特征向量的求解根据步骤<1>得到的P和D、步骤<3>得到的Q、步骤<4>得到的V和U，回推计算S关于Λ的左特征向量V和右特征向量U：$\left(\begin{array}{c}U=\mathrm{PDQU}\\ V={\mathrm{PD}}^{-1}\mathrm{QV}\end{array}\right);$

<6>将计算得到的全部特征值和部分左/右特征向量，以2进制的形式输入到文本文件中。 

现有技术中能够完成上述计算步骤的数值计算库非常多，例如：Intel的MKL、AMD的ACML、Netlib的参考BLAS/LAPACK、SourceForge的ATLAS/LAPACK、Matlab等。为使得TTQRE算法能够在所部署的计算机上发挥出最佳的计算性能，本实施例采用中国科学院软件所张先轶团队提供的OpenBLAS v0.2.8，以及Netlib提供的参考LAPACK v3.5.0。该模块的实际应用效果将在上述实施例中给出。 

OpenBLAS是一套针对于不同CPU内核架构的优化BLAS子程序库。它是张先轶团队在後藤和茂开发的GotoBLAS2的基础上，逐渐开发与完善的多线程并行BLAS子程序库。对于不同的CPU微架构，例如：Penryn、Athlon、SandyBridge等，OpenBLAS的计算速度比MKL、ACML、ATLAS、Matlab等要快，并且矩阵规模越大，OpenBLAS的加速效果越明显。 

LAPACK是在BLAS子程序库的基础上，采用Fortran77开发的一套专用于稠密矩阵数值计算的线性代数方法。LAPACK是早期的线性代数方法LINPACK以及特征值方法EISPACK的集合体，涵盖了LINPACK中的线性方程组求解、线性最小二乘逼近以及EISPACK中的特征值计算等功能。另外，3.1.0及其以后版本的LAPACK，封装的就是TTQRE。 

所述的部分特征值计算引擎模块：以RCI函数回调的方式，执行IRAM或者JDM的数值迭代过程根据稀疏矩阵计算模块提供的API，实现与稀疏矩阵相关的所有运算。该模块包含以下11个主要步骤： 

<1>读入系统建模及其线性化模块生成的J_A、J_B、J_C和J_D； 

<2>拼接J_aug并进行符号分析，保留符号分析的重排矩阵Q_aug； 

<3>对J_D进行符号分析以及LU数值分解； 

<4>从*.cfg文件中读入部分特征值计算引擎的配置信息，包括：所选算法、谱变换类型、Ritz值选择策略、搜索的特征值个数、最小/最大搜索子空间的维数、收敛精度、最大迭代 次数、平移点σ以及反平移点μ； 

<5>利用平移点σ来修正J_aug的部分对角元，并对修正后的J_aug–σI进行LU数值分解：P_aug(J_aug-σI)Q_aug＝L_augU_aug；其中：P_aug和Q_aug分别为实现列选主元和近似最小度排序的重排矩阵。 

<6>若所选算法为IRAM，转入步骤<7>，否则所选算法为JDM，转入步骤<8>； 

<7>调用IRAM的RCI函数，执行IRAM的数值迭代，转入步骤<9>； 

<8>调用JDM的RCI函数，执行JDM的数值迭代，转入步骤<9>； 

<9>判断数值迭代是否结束，若是，则转入步骤<11>，若否，则转入步骤<10>； 

<10>根据谱变换算子Op的类型根据前述步骤中保存的中间变量，实现RCI函数回调过程中两类矩阵‐向量乘法的隐式稀疏求解，求解完毕后回到步骤<6>； 

所述的两类矩阵‐向量乘法是指：其中：S为系统的状态矩阵，x为IRAM或JDM提供的RCI回调向量，θ为JDM提供的RCI回调标量，y为提供给IRAM或JDM的待求向量。 

<11>若算法收敛得到了一部分特征值及其右特征向量，则采用两次反幂法来求解该部分特征值对应的左特征向量，并将所有的计算结果以2进制的形式输入到文本文件中。 

如图2a和图2b(IRAM和JDM的数值迭代过程示意图的左半部分和右半部分)给出了上述步骤<6>～步骤<10>的具体实施过程，由图2a和图2b可以看出，RCI函数回调的执行过程是非常复杂的，但正如Jack Dongarra教授在文献[Dongarra J，Eijkhout V，Kalhan A.Reverse communication interface for linear algebra templates for iterative methods[J].UT，CS‐95‐291，May，1995.]的引言中所述：“逆通信是一项技术，通过它我们可以将迭代算法中各种操作的实现细节给隐藏起来。一方面，算法开发者可以无需考虑用户是采用什么样的数据结构来保存迭代求解器所需的矩阵，另一方面，用户可以根据自己想要的方式，来实现迭代求解器所要求的运算”。虽然本实施例仅提供了IRAM和JDM两种特征值算法的具体实施方式，但可以看出，对于其他的迭代投影法，只要稍作修改将其部署为RCI形式，就可以直接添加到部分特征值计算引擎模块中，因此，本实施例具有非常强的可扩展性。 

所述步骤<10>中，两类矩阵‐向量乘法的隐式稀疏求解，是通过调用稀疏矩阵计算模块中“截断的稀疏矩阵‐向量乘法”这一扩展矩阵运算功能来实现的。 

以Cayley变换为例，即Op(S)＝(S–μI)×(S–σI)^‐1，IRAM和JDM均采用以下的步骤来实现第1类矩阵‐向量乘法的隐式稀疏求解。 

<a>用x来初始化与J_aug同维的列向量z＝[x；0]； 

<b>对z执行稀疏前代/回代：z＝Q_aug×(U_aug\(L_aug\(P_aug×z))) 

<c>记n为S的阶数，取出z的前n行，生成列向量y＝z(1:n)； 

<d>调用2级BLAS子程序gemv，计算y＝x+(σ–μ)×y。 

以平移‐逆变换为例，即Op(S)＝(S–σI)^‐1，JDM采用以下的步骤来实现第2类矩阵‐向量乘法的隐式稀疏求解。 

<e>计算临时向量v：v＝Q_D×(U_D\(L_D\(P_D×(-J_C×x)))) 

<f>计算临时向量u：$u=-\frac{(J_A-\mathrm{\sigma I})\times x+J_B\times v}{\theta }$

<g>利用临时平移点α＝(1+θ×σ)/θ修正J_aug的部分对角元，并对修正后的J_aug–αI进行LU数值分解：${\hat{P}}_{\mathrm{aug}}(J_{\mathrm{aug}}-\mathrm{\alpha I})Q_{\mathrm{aug}}={\hat{L}}_{\mathrm{aug}}{\hat{U}}_{\mathrm{aug}}$

<h>用u来初始化与J_aug同维的列向量z＝[u；0]； 

<i>对z执行稀疏前代/回代：$Q_{\mathrm{aug}}\times ({\hat{U}}_{\mathrm{aug}}\backslash ({\hat{L}}_{\mathrm{aug}}\backslash ({\hat{P}}_{\mathrm{aug}}\times z)))$

<j>记n为S的阶数，取出z的前n行，生成列向量y＝z(1:n)。 

所述的振荡模式提取和分析模块：读取全部特征值计算引擎模块或者部分特征值计算引擎模块的计算结果，实现机电振荡模式的识别和模态分析功能，最终将计算所得的模式和模态信息，以Excel表格的形式输出。 

为了进一步说明本实施例系统的正确性和有效性，分别对3个IEEE标准系统和3个实际系统的小干扰稳定性进行了仿真分析与计算。其中，3个实际系统分别取自华东电网05、09和15年的夏高规划数据。这6个测试系统的详细信息如表1所示。 

表1测试系统描述 

测试系统 节点个数 支路个数 发电机个数 状态变量个数 代数变量个数 IEEE118 118 186 34 340 236 IEEE300 300 411 69 690 600 IEEE600 600 823 139 1390 1200 HD4171 4171 5318 375 3750 8342 HD5473 5473 7440 472 4720 10946 HD8241 8241 11657 500 3000 16482

由于测试系统的数据均存储为PSD‐BPA的文件格式，因此，借助PSD‐SSAP v2.5.2作为精度的标准参考，借助Matlab R2012a作为速度的标准参考，验证本实施例提供的小干扰稳定 性特征值分析系统的正确性和有效性。所有测试均在同一个计算环境下进行：Windows XP操作系统、4GB内存、Intel 2.5GHz的双核CPU E5200。由于该CPU为Penryn微架构，因此，在OpenBLAS的makefile.rule文件中，指定TARGET＝PENRYN。由于PSD‐SSAP为串行计算模式，为公平起见，本实施例和Matlab也指定为单核心CPU的串行计算模式。 

正确性验证包含以下2个方面： 

1)分别采用本实施例的TTQRE和PSD‐SSAP的QR法，计算6个测试系统的全部特征值和左/右特征向量，将识别得到的机电振荡模式进行两两匹配后，比较计算结果的最大绝对偏差，结果如表2所示； 

2)分别采用本实施例的平移‐逆变换IRAM和PSD‐SSAP的平移‐逆变换IRAM，计算6个测试系统的部分特征值和左/右特征向量，将识别得到的机电振荡模式进行两两匹配后，比较计算结果的绝对偏差，结果如表3所示。扫频IRAM采用PSD‐SSAP默认的参数配置。 

表2全部特征值分析的正确性验证 

表3基于平移‐逆变换扫频IRAM部分特征值分析的正确性验证 

由表2和表3可以看出，本实施例是正确的。产生振荡模式个数出现偏差的原因在于， 

本实施例和PSD‐SSAP采用了不同的方法来计算相关因子。 

有效性验证包含以下3个方面： 

1)分别采用本实施例的状态矩阵计算模块、PSD‐SSAP、Matlab的稀疏矩阵计算功能，显示计算6个测试系统的状态矩阵，三者的计算时间如表4所示； 

2)分别采用本实施例的TTQRE、PSD‐SSAP的QR法、Matlab的eig函数，计算6个测试系统的全部特征值和左/右特征向量，三者的计算时间如表5所示； 

3)分别采用本实施例的平移‐逆变换IRAM和PSD‐SSAP的平移‐逆变换IRAM，计算6个测试系统的部分特征值和左/右特征向量，两者的计算时间如表6所示。同样地，扫频IRAM采用PSD‐SSAP默认的参数配置。 

表4状态矩阵计算的有效性验证 

表5全部特征值和特征向量计算的有效性验证 

测试系统 本实施例的TTQRE/s PSD‐SSAP的QR法/s Matlab的eig函数/s IEEE118 0.203 0.540 0.408 IEEE300 1.705 4.190 2.029 IEEE600 13.144 31.115 14.371 HD4171 236.007 694.957 253.122 HD5473 450.542 1812.102 482.496 HD8241 128.400 376.078 139.143

表6基于平移‐逆变换扫频IRAM部分特征值分析的有效性验证 

测试系统 本实施例的平移‐逆变换扫频IRAM/s PSD‐SSAP的平移‐逆变换扫频IRAM/s IEEE118 1.058 1.206 IEEE300 1.523 1.618 IEEE600 3.224 3.388 HD4171 63.072 76.291

 

HD5473 72.184 85.682 HD8241 24.825 27.071

表4充分证明了本实施例的稀疏矩阵计算模块中，“计算稀疏矩阵关于某个子块的Schur补”这一扩展矩阵运算功能的有效性。表5充分证明了本实施例的全部特征值计算引擎模块中，TTQRE算法的有效性。表6充分证明了稀疏矩阵计算模块中，“截断的稀疏矩阵‐向量乘法”这一扩展矩阵运算功能的有效性。 

为了考察本实施例中提供的Cayley变换结合模最大Ritz值选择策略的IRAM单元、原始矩阵结合阻尼比最小Ritz值选择策略的JDM单元，对于识别大规模电力系统关键振荡模式的实际应用效果，取上述实施例中的HD8241系统进行仿真测试。关键振荡模式指的是阻尼比不大于5％的特征值。 

图3所示为平移‐逆变换扫频IRAM识别出的机电振荡模式。扫频IRA识别出115个模式，其中72个为阻尼比小于5％的关键模式，而全部模式有499个，其中87个为关键模式，因此，扫频IRA方法出现了“漏根”，特别是在1.2Hz附近。 

图4所示为Cayley变换结合模最大Ritz值选择策略IRAM识别出的机电振荡模式。Cayley变换IRAM一共迭代70次，耗时60.058秒，识别出了全部87个阻尼比小于5％的关键机电振荡模式，因此，没有出现“漏根”。 

图5所示为原始矩阵结合阻尼比最小Ritz值选择策略JDM识别出的机电振荡模式。该计算方案的JDM一共迭代109次，耗时208.007秒，识别出了全部87个阻尼比小于5％的关键机电振荡模式，因此，没有出现“漏根”。 

图6所示为原始矩阵结合阻尼比最小Ritz值选择策略JDM的数值迭代收敛过程。虽然本实施例中JDM的计算效率不及IRAM，但从该图可以看出，JDM能够按照阻尼比由小到大的顺序，依次收敛出关键振荡模式。 

由于本系统识别出的机电振荡模式和模态信息均以Excel表格的形式输出，因此，用户既可以采用Excel自带的功能，又可以采用其他的方式来使用计算结果。 

仍以HD8241系统为例，给出本系统仿真分析结果的两种应用场景： 

1)在Excel中将原始矩阵结合阻尼比最小Ritz值选择策略JDM和平移‐逆变换扫频IRAM得到的机电振荡模式按照阻尼比升序排列后，发现IRAM漏掉了一个阻尼比为0.04303542、频率为1.19297526Hz的关键振荡模式，将该振荡模式的相关因子按照模值降序排列后，发现与该振荡模式最相关的机组是”ZCN_1”和”ZCN_2”，于是，在这两台机组上分别加装电力系统稳定器，可以有效提高HD8241系统的动态稳定性； 

2)在Matlab中读取保存在Excel文件中的JDM计算得到的模态信息，找到阻尼比为0.02234904、频率为0.62502012Hz的关键振荡模式，由于该模式的振荡频率较低，可以判断 其为区域间振荡模式，为了找出该模式下的振荡区域，对该模式的模态信息进行作图结果如图7所示，由图7可以看出在该模式下，“振荡区域1”和“振荡区域2”的机组之间进行相对摇摆，因此，分析结果表明这两个振荡区域内的所有机组可以进行动态等值。