大规模电力系统小干扰稳定性的特征值分析系统的制作方法
【专利摘要】一种电力系统仿真及分析【技术领域】的大规模电力系统小干扰稳定性的特征值分析系统，包括：用于实现BPA数据文件的无缝读写操作的BPA数据接口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵计算模块、潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算引擎模块、部分特征值计算引擎模块和振荡模式提取和分析模块，采用TTQRE、多种计算方案的IRAM及JDM，实现对实际大规模电网进行全方位、无遗漏的小干扰稳定性特征值仿真分析。
【专利说明】大规模电力系统小干扰稳定性的特征值分析系统
【技术领域】
[0001] 本发明涉及的是一种电力系统仿真及分析【技术领域】的系统，具体是一种采用三种 不同特征值算法的大规模电力系统小干扰稳定性的特征值分析系统。
【背景技术】
[0002] 电力系统小干扰稳定性指的是当电网经历微小扰动后，继续保持同步运行的能 力。电力系统小干扰稳定性一般采用李雅普诺夫第一法作为判定标准。李雅普诺夫第一法 指出，如果系统线性化后的状态矩阵没有出现零或者正实部的特征值，那么就可以判定当 前系统是小干扰稳定的。因此，在电力系统中，特征值计算一直以来都是实现低频功率振荡 模式的识别、各类稳定控制器的装配选址和参数优化、运行参数对控制参数的灵敏度分析、 在线检测振荡数据的模态信息提取等功能的重要前提和基础保障。
[0003] 由于电网规模的不断扩大，自上个世纪80年代开始，电力行业的研究人员就已投 入大量的精力到特征值计算方法和小干扰稳定性分析系统的研究与开发中。至今为止，世 界范围内已出现了诸多用于大规模电力系统小干扰稳定性的特征值分析系统，例如：美国 太平洋瓦斯与电力公司开发的EISEMAN、美国电科院联合加拿大安大略水电局共同开发的 SSSP、加拿大动力技术实验室开发的SSAT、德国西门子公司开发的NEVA、巴西电力科学研 究中心开发的PacDyn、中国电力科学研究院开发的PSASP和PSD - SSAP等。
[0004] 按照待求特征值的数量进行归类，电力系统小干扰稳定性的特征值分析方法分为 全部特征值分析法和部分特征值分析法。虽然几乎所有的大规模电力系统小干扰稳定性 分析系统均包含全部特征值分析和部分特征值分析两个方面，但仍然存在以下的问题和不 足：
[0005] 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,Demme1 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 = 23°GB)寻址的64位。因此，基于双重步位移隐式QR算法和32 位个人计算机的全部特征值分析法，显然已经不能满足当前大规模电力系统小干扰稳定性 分析的计算需求；
[0006] 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 Η· 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， ffu Y Κ. 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. ] 〇
[0007] 虽然IRAM已被广泛集成于各电力系统小干扰稳定性分析系统中，但其算法本身 具备的功能却没能得到最大程度的发挥，例如：PSD - SSAP和PSASP只提供了基于平移-逆 变换的扫频IRAM，而非常适合于计算关键振荡模式的Cayley变换IRAM，在两种方法中均没 有体现。扫频IRAM虽然对于极特征值有着线性的收敛速度，但其搜索行为存在着随机性， 在参数配置不合理的情况下，容易出现"漏根"，即漏掉一些特别关键的振荡模式。当精确 求解修正方程时，JDM具有渐进二阶的收敛速度，但更重要的是，即使不使用谱变换，JDM仍 然能够按照某种指定的规则（实部最大、阻尼比最小），顺序收敛出影响电力系统小干扰稳 定性的那些关键振荡模式。
[0008] 然而，不像IRAM具有开源的数学库ARPACK，由于JDM仅有Matlab下的示例程序， 因此，JDM-直常见于学术研究，而没有集成到任何的大规模电力系统小干扰稳定性特征值 分析系统中。
[0009] "三华"（华中-华东-华北）特高压交/直流电网的互联，一定会出现新增的低 频功率振荡模式，而如今处于建设初期的"强直弱交"电网，极有可能使得其中的某些模式 将要或者已经呈现出不稳定的状态，因此，采用三种数值原理完全不同的特征值算法来分 析大规模电力系统的小干扰稳定性，不仅具有理论价值，更具有现实意义。
【发明内容】
[0010] 本发明针对现有技术存在的上述不足，提出一种大规模电力系统小干扰稳 定性的特征值分析系统，根据PSD - BPA格式的电网数据文件（*. dat和*. swi)，采用 TTQRE(Two - tone small - bulge multishift QR algorithm with aggressive early deflation，带有积极早期收缩策略的两步小块多重步位移QR算法），快速、准确的识别出 影响系统动态稳定的全部机电振荡模式；采用具有多种计算方案的IRAM，提取和分析影响 系统动态稳定的部分关键机电振荡模式；采用具有多种计算方案的JDM，提取和分析影响 系统动态稳定的部分关键机电振荡模式，通过匹配和比较三种算法得到的机电振荡模式， 最终实现对实际大规模电网进行全方位、无遗漏的小干扰动态稳定的振荡模式识别和模态 分析。
[0011] 本发明是通过以下技术方案实现的：
[0012] 本发明涉及一种大规模电力系统小干扰稳定性的特征值分析系统，包括：BPA数 据接口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵计 算模块、潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算 引擎模块、部分特征值计算引擎模块和振荡模式提取和分析模块，其中：BPA数据接口模 块接收电网数据文件，并分别与潮流计算模块和系统建模及其线性化模块相连并提供 API (Application Program Interface,应用程序接口）；稀疏矩阵计算模块分别与潮流计 算模块、系统建模及其线性化模块、状态矩阵计算模块和部分特征值计算引擎模块相连并 提供API ;潮流计算模块根据BPA数据接口模块以及稀疏矩阵计算模块提供的API，实现大 规模电力系统的潮流计算并输出潮流计算结果至系统建模及其线性化模块；系统建模及其 线性化模块根据BPA数据接口模块提供的API和潮流计算模块提供的潮流计算结果，实现 大规模电力系统的线性化建模，即生成增广的系统状态矩阵并分别输出至状态矩阵计算模 块和部分特征值计算引擎模块；状态矩阵计算模块根据稀疏矩阵计算模块提供的API和 系统建模及其线性化模块提供的增广的系统状态矩阵，获得系统状态矩阵并输出至全部特 征值计算引擎模块；全部特征值计算引擎模块根据TTQRE从系统状态矩阵中获得系统状 态矩阵的全部特征值以及部分左/右特征向量并输出至振荡模式提取和分析模块；部分 特征值计算引擎模块根据稀疏矩阵计算模块提供的API，以RCI (Reverse Communication Interface，逆通信接口）的方式执行IRAM或者JDM的数值迭代过程，获得系统状态矩阵的 部分特征值和部分左/右特征向量并输出至振荡模式提取和分析模块；振荡模式提取和分 析模块根据用户指定合并来自全部特征值计算引擎模块或者部分特征值计算引擎模块的 系统状态矩阵的特征值和左/右特征向量，实现机电振荡模式的识别和模态分析，并将所 得的模式信息和模态分析结果以Excel表格的方式输出给用户。 …/ /
[0013] 所述的增广的系统状态矩阵Jaug满足/ / 其中小为动态元件线性 \_J ( JD 化模型拼接而成的分块对角矩阵，1为动态元件与静态元件的联系矩阵，1为静态元件与 动态元件的联系矩阵，JD为修正的系统导纳矩阵。
[0014] 所述的系统状态矩阵S是指：S
[0015] 所述的全部特征值计算引擎模块内包括TTQRE单元，该单元与状态矩阵计算模块 相连，并根据系统状态矩阵获得系统状态矩阵的全部特征值和部分左/右特征向量后，输 出至振荡模式提取和分析模块。
[0016] 所述的部分特征值计算引擎模块包括：IRAM单元和JDM单元，其中：IRAM单元与 系统建模及其线性化模块相连，根据用户的选择来确定采用平移-逆变换结合模最大Ritz 值选择策略还是Cayley变换结合模最大Ritz值选择策略的计算方案，并根据增广的系统 状态矩阵获得系统状态矩阵的部分特征值和部分左/右特征向量后，输出至振荡模式提取 和分析模块；JDM单元与系统建模及其线性化模块相连，根据用户的选择来确定采用平移-逆变换结合模最大Ritz值选择策略、Cayley变换结合模最大Ritz值选择策略、原始矩阵结 合实部最大Ritz值选择策略还是原始矩阵结合阻尼比最小Ritz值选择策略的计算方案， 并根据增广的系统状态矩阵获得系统状态矩阵的部分特征值和部分左/右特征向量后，输 出至振荡模式提取和分析模块。
[0017] 所述的振荡模式识别和模态分析是指：根据每一对左/右特征向量获得相应特征 值的相关因子，根据相关因子获得该特征值的机电回路相关比，若机电回路相关比大于1 且该特征值为复数，则判定该特征值为机电振荡模式，相应的右特征向量为振荡模态；相关 因子模值的大小反应了电力系统动态元件与振荡模式之间的耦合程度，在耦合程度最大的 机组上加装电力系统稳定器，是抑制由该振荡模式诱发的电力系统低频功率振荡的一贯做 法；振荡模态的幅值决定了该振荡模式下机组的摇摆程度，角度决定了该振荡模式下机组 与机组之间的相对摇摆关系，振荡模态是判断振荡模式对机组的影响程度、判断多个机组 是否可以进行动态等值的重要指标。 技术效果
[0018] 与现有技术手段相比，本发明的有益效果包括：
[0019] 1.采用了能够代表当前技术水平的TTQRE，来实现大规模电力系统小干扰稳定性 的全部特征值分析。该算法在特征值迭代求解环节，采用链式小块追赶技术，消除了多重步 位移QR块追赶的位移模糊现象，解决了算法不收敛的问题；采用积极早期收缩策略，减少 了 QR迭代的次数，提高了算法的计算精度；采用矩阵-矩阵的高级运算，提高了 CPU缓存的 利用率，节省了算法的计算时间。从工程实用的角度出发，论证了 QR算法对于大型电力系 统小干扰稳定性全部特征值分析的可行性；
[0020] 2.将IRAM的谱变换统一描述为算子形式，实现了具有平移-逆变换结合模最大 Ritz值选择策略、Cayley变换结合模最大Ritz值选择策略两种计算方案的IRAM。Cayley 变换IRAM的最大优点在于，能够有效识别电力系统中阻尼比小于等于某个指定值的一些 关键振荡模式，相比平移-逆变换IRAM，Cayley变换IRAM的搜索行为基本不存在随机性， 因此，即便是在参数配置不完全合理的情况下，Cayley变换IRAM仍可以识别出绝大部分 的关键振荡模式，很大程度上降低了用IRAM分析大规模电力系统的小干扰稳定性时，出现 "漏根"的可能性；
[0021] 3.在国内外首次以RCI的方式实现了 JDM。RCI不仅完美诠释了特征值算法与电 力系统应用之间的分界点，RCI还使得迭代投影法具备了一致的调用规则。相比Krylov子 空间下的迭代投影法，非Krylov子空间下的JDM不仅具有更快的收敛速度，而且在不使用 谱变换的情况下，JDM仍可以按照某种指定的Ritz值选择策略（阻尼比最小、实部最大）， 依次收敛出想要的特征值。实现了具有4种计算方案的JDM，并将其与IRAM -起集成到了 大规模电力系统小干扰稳定性特征值分析系统中，对JDM的工程实用化起到了推动性的作 用。
【专利附图】
【附图说明】
[0022] 图1为本发明提供的技术方案实施图；
[0023] 图2a和图2b分别为IRAM和JDM的数值迭代过程示意图的左半部分和右半部分；
[0024] 图3为平移-逆变换扫频IRAM的应用效果图；
[0025] 图4为Cayley变换结合模最大Ritz值选择策略IRAM的应用效果图；
[0026] 图5为原始矩阵结合阻尼比最小Ritz值选择策略JDM的应用效果图；
[0027] 图6为原始矩阵结合阻尼比最小Ritz值选择策略JDM的收敛过程图；
[0028] 图7为HD8241系统中某一区间振荡模式的模态分析图。
【具体实施方式】
[0029] 下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行 实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施 例。 实施例1
[0030] 如图1所示，本实施例包括：用于实现BPA数据文件的无缝读写操作的BPA数据接 口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵计算模块、 潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算引擎模块、 部分特征值计算引擎模块和振荡模式提取和分析模块，其中：
[0031] 所述的BPA数据接口模块：由于国内的电网企业偏向于使用BPA的文件格式来 存放实际电网数据，为此，基于面向对象编程思想，采用C++开发了具有动态链接库形式的 BPA数据接口模块。该模块不仅具有BPA数据文件的无缝读写功能，还可以为其他模块提供 API。该模块的【具体实施方式】包括以下3个主要步骤：
[0032] <1>BPA数据文件的层次结构分析，用来确定电网的静态/动态元件参数在BPA数 据文件中的存放位置和存放顺序，实现数据卡片和控制卡片的有效分离；
[0033] 〈2>类结构框架的设计，按照面向对象编程的要求，设计一套与BPA程序用户手册 具有相同层次关系的类结构框架，方便接口模块的维护与扩展；
[0034] <3>BPA数据接口的生成，采用C++对所设计的类结构框架进行编程实现，生成的 数据接口不仅可以提供API，更可以提供独立的应用程序。
[0035] 所述的稀疏矩阵计算模块：任何大规模电力系统的仿真分析与计算，都离不开稀 疏技术。该模块采用了佛罗里达州立大学Timothy A. Davis教授提供的SuiteSparse，来实 现与稀疏矩阵相关的所有运算，以动态链接库的形式为其他模块提供API。该模块主要由5 部分构成：
[0036] 〈1>数据输入/输出，实现稀疏矩阵从硬盘到内存之间的双向数据传输，该部分既 支持2进制格式保存的稀疏矩阵，同时又满足32/64位寻址的需求；
[0037] 〈2>基本矩阵运算，实现稀疏矩阵的一些常规计算，例如：转置、重排、矩阵-向量 乘法、矩阵-矩阵加法、矩阵-矩阵乘法等；
[0038] 〈3>稀疏LU分解，实现稀疏矩阵的符号分析和LU数值分解；
[0039] 〈4>稀疏前代/回代，实现稀疏矩阵的快速前代/回代；
[0040] 〈5>扩展矩阵运算，实现一些扩展的稀疏矩阵运算，例如：计算稀疏矩阵关于某个 子块的Schur补、截断的稀疏矩阵-向量乘法等。
[0041] 现有技术中SuiteSparse是与SuperLU、PARDISO、TAUCS等齐名的稀疏矩阵解算 器，更重要的是，在SuiteSparse源码的基础上，能够开发出专用于电力系统小干扰稳定性 特征值分析的稀疏矩阵计算功能，例如，稀疏矩阵计算模块的第〈4>和第〈5>部分所示的功 能。这些扩展稀疏矩阵计算功能的实际应用效果，将在后面给出。
[0042] 所述的潮流计算模块：利用BPA数据接口模块以及稀疏矩阵计算模块提供的API， 实现大规模电力系统的潮流计算，为系统建模及其线性化模块提供潮流计算结果。
[0043] 所述的系统建模及其线性化模块根据BPA数据接口模块提供的API和潮流计算 模块提供的潮流计算结果，并结合文献[励刚，陈陈.应用设计模式开发小干扰稳定性分 析方法[刀.中国电机工程学报，2002,22(1) :12 - 16.]提出的统一元件联接建模方法， 来对小干扰稳定性分析进行系统线性化建模，同时生成增广的系统状态矩阵Jaug满足： 人^ ^]其中：JA为动态元件线性化模型拼接而成的分块对角矩阵，JB为动态元件 与静态元件的联系矩阵，J。为静态元件与动态元件的联系矩阵，JD为修正的系统导纳矩阵。
[0044] 所述的状态矩阵计算模块：利用稀疏矩阵计算模块提供的API和系统建模 及其线性化模块提供的稀疏矩阵J A、JB、J。和JD，实现系统状态矩阵S的稀疏求解： X m:由于S是稠密矩阵，而^、^是稀疏矩阵，根据稠密矩阵的计算方法 来求解S，一方面计算时间会很长，另一方面内存开销也会很大，于是，该模块采用以下步骤 求解S :
[0045] 〈1>对JD进行符号分析以及稀疏LU分解：PDJ DQD = LDUD ;其中：PD和QD分别为实现 列选主元和近似最小度排序的重排矩阵。
[0046] 〈2>用JA来初始化稠密存储的S :S = full (JA)
[0047] 〈3>利用稀疏矩阵计算模块提供的扩展矩阵运算功能，计算稀疏矩阵Jaug关于子 块J D的Schur补，=S - (Λ仏)?^/砧（):由于电力系统小干扰稳定性分析中的JB和 J。具有非常特殊的稀疏结构，而使得第〈3>步公式的右端项并不需要执行全部的数值计算， 因此，稀疏矩阵计算模块才会在SuiteSparse的源码基础上，实现了"计算稀疏矩阵关于某 个子块的Schur补"的扩展矩阵运算功能。该功能的实际应用效果将在后面给出。
[0048] 所述的全部特征值计算引擎模块：采用TTQRE，实现对状态矩阵S的全部特征值、 部分左/右特征向量的求解，包括以下6个主要步骤：
[0049] 〈1>矩阵平衡，采用矩阵重排分离出已经暴露于S对角线上的部分实特征值，采用 循环缩放减小平衡后矩阵A的2范数：a= (rors(ro);其中：p为重排矩阵，d为对角矩 阵。
[0050] <2>Hessenberg 约化，米用一系列 Householder 反射，将 A 约化为上 Hessenberg 形 式H :H = QTAQ ;其中：Q为一系列Householder反射累积而成的正交矩阵。
[0051] 〈3>链式小块追赶的多重步位移隐式QR迭代，采用一系列Householder反射和 Givens旋转执行链式QR块追赶，采用收缩窗收敛多个特征值，将Η相似变换为拟上三角矩 阵T，S的全部特征值计算结果即为Τ的所有对角元：T = QTAQ ;其中：Q为进一步累积而成 的正交矩阵，包含一系列Householder反射和Givens旋转。
[0052] 〈4>标准Schur型部分特征向量的求解，确定T中想要的特征值Λ后，计算T关于 \VHT = AVH Λ的左特征向量V和右特征向量υ:^ V · *
[0053] 〈5>状态矩阵部分特征向量的求解根据步骤〈1>得到的P和D、步骤〈3>得到的Q、 \l· Π )〇1 步骤〈4>得到的V和U，回推计算S关于Λ的左特征向量V和右特征向量=
[0054] 〈6>将计算得到的全部特征值和部分左/右特征向量，以2进制的形式输入到文本 文件中。
[0055] 现有技术中能够完成上述计算步骤的数值计算库非常多，例如：Intel的MKL、AMD 的 ACML、Netlib 的参考 BLAS/LAPACK、SourceForge 的 ATLAS/LAPACK、Matlab 等。为使得 TTQRE算法能够在所部署的计算机上发挥出最佳的计算性能，本实施例采用中国科学院软 件所张先轶团队提供的(^61^1^3 ￥0.2.8，以及他丨1让提供的参考1^^4〇(￥3.5.0。该模 块的实际应用效果将在上述实施例中给出。
[0056] OpenBLAS是一套针对于不同CPU内核架构的优化BLAS子程序库。它是张先轶团 队在後藤和茂开发的GotoBLAS2的基础上，逐渐开发与完善的多线程并行BLAS子程序库。 对于不同的CPU微架构，例如：Penryn、Athlon、SandyBridge等，OpenBLAS的计算速度比 MKL、ACML、ATLAS、Matlab等要快，并且矩阵规模越大，OpenBLAS的加速效果越明显。
[0057] LAPACK是在BLAS子程序库的基础上，采用Fortran77开发的一套专用于稠密矩 阵数值计算的线性代数方法。LAPACK是早期的线性代数方法UNPACK以及特征值方法 EISPACK的集合体，涵盖了 UNPACK中的线性方程组求解、线性最小二乘逼近以及EISPACK 中的特征值计算等功能。另外，3. 1. 0及其以后版本的LAPACK，封装的就是TTQRE。
[0058] 所述的部分特征值计算引擎模块：以RCI函数回调的方式，执行IRAM或者JDM的 数值迭代过程根据稀疏矩阵计算模块提供的API，实现与稀疏矩阵相关的所有运算。该模块 包含以下11个主要步骤：
[0059] 〈1>读入系统建模及其线性化模块生成的JA、JB、J。和J D ;
[0060] 〈2>拼接Jaug并进行符号分析，保留符号分析的重排矩阵Qaug ;
[0061] 〈3>对JD进行符号分析以及LU数值分解；
[0062] 〈4>从*. cfg文件中读入部分特征值计算引擎的配置信息，包括：所选算法、谱变 换类型、Ritz值选择策略、搜索的特征值个数、最小/最大搜索子空间的维数、收敛精度、最 大迭代次数、平移点σ以及反平移点μ ;
[0063] 〈5>利用平移点σ来修正Jaug的部分对角元，并对修正后的Jaug - σ I进行LU数 值分解：Paug(Jaug-〇 I)Qaug = Lauguaug ;其中：Paug和Qaug分别为实现列选主元和近似最小度排 序的重排矩阵。
[0064] 〈6>若所选算法为IRAM，转入步骤〈7>，否则所选算法为JDM，转入步骤〈8> ;
[0065] 〈7>调用IRAM的RCI函数，执行IRAM的数值迭代，转入步骤〈9> ;
[0066] 〈8>调用JDM的RCI函数，执行JDM的数值迭代，转入步骤〈9> ;
[0067] 〈9>判断数值迭代是否结束，若是，则转入步骤〈11>，若否，则转入步骤〈10> ;
[0068] 〈10>根据谱变换算子Op的类型根据前述步骤中保存的中间变量，实现RCI函数回 调过程中两类矩阵-向量乘法的隐式稀疏求解，求解完毕后回到步骤〈6> ; 第-.......炎：Op(S)xx = j
[0069] 所述的两类矩阵-向量乘法是指：j Λ ;其中：S为系 -灸：L〇p(s)-併_| xx=j 统的状态矩阵，X为IRAM或JDM提供的RCI回调向量，Θ为JDM提供的RCI回调标量，y为 提供给IRAM或JDM的待求向量。
[0070] 〈11>若算法收敛得到了一部分特征值及其右特征向量，则采用两次反幂法来求 解该部分特征值对应的左特征向量，并将所有的计算结果以2进制的形式输入到文本文件 中。
[0071] 如图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形式，就可以直接添加到部分特征值计算引擎模块中，因此，本实施例具有非常强的 可扩展性。
[0072] 所述步骤〈10>中，两类矩阵-向量乘法的隐式稀疏求解，是通过调用稀疏矩阵计 算模块中"截断的稀疏矩阵-向量乘法"这一扩展矩阵运算功能来实现的。
[0073] 以Cayley变换为例，即Op⑶=(S - μ I) X (S - σ I) ' IRAM和JDM均采用以下 的步骤来实现第1类矩阵-向量乘法的隐式稀疏求解。
[0074] <a>用X来初始化与Jaug同维的列向量z = [X ;0];
[0075] <b> 对 z 执行稀疏前代 / 回代：z = QaugX (Uaug\ (Laug\ (PaugX z)))
[0076] <c>记n为S的阶数，取出z的前n行，生成列向量y = z(l:n);
[0077] <d> 调用 2 级 BLAS 子程序 gemv，计算 y = x+(〇 -μ)Χγ。
[0078] 以平移-逆变换为例，即Op (S) = (S - 〇 I) JDM采用以下的步骤来实现第2类 矩阵-向量乘法的隐式稀疏求解。
[0079] <e> 计算临时向量 v :v = QDX (UD\(LD\(PDX (_JcXx))))
[0080] <f> 计算临时向量 u : η = -^σ~'）x + JB^ρ θ
[0081] <g>利用临时平移点α = (1+ θ χ σ )/ θ修正Jaug的部分对角元，并对修正后的 Jaug - α I 进行 LU 数值分解：K人= Κ、，
[0082] <h>用u来初始化与Jaug同维的列向量ζ = [u ;0];
[0083] <i> 对 ζ 执行稀疏前代 / 回代：f = &·? '(l Xf)))
[0084] <j>记n为S的阶数，取出z的前n行，生成列向量y = z(l:n)。
[0085] 所述的振荡模式提取和分析模块：读取全部特征值计算引擎模块或者部分特征值 计算引擎模块的计算结果，实现机电振荡模式的识别和模态分析功能，最终将计算所得的 模式和模态信息，以Excel表格的形式输出。
[0086] 为了进一步说明本实施例系统的正确性和有效性，分别对3个IEEE标准系统和3 个实际系统的小干扰稳定性进行了仿真分析与计算。其中，3个实际系统分别取自华东电网 05、09和15年的夏高规划数据。这6个测试系统的详细信息如表1所示。
[0087] 表1测试系统描述 jji式系统I节点个数I支路个数I发电机个数I状态变量个数I代数变量个数 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 |l!657 丨500 丨3000 |l6482
[0088] 由于测试系统的数据均存储为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的串行计算模式。
[0089] 正确性验证包含以下2个方面：
[0090] 1)分别采用本实施例的TTQRE和PSD -SSAP的QR法，计算6个测试系统的全部特 征值和左/右特征向量，将识别得到的机电振荡模式进行两两匹配后，比较计算结果的最 大绝对偏差，结果如表2所示；
[0091] 2)分别采用本实施例的平移-逆变换IRAM和PSD - SSAP的平移-逆变换IRAM， 计算6个测试系统的部分特征值和左/右特征向量，将识别得到的机电振荡模式进行两两 匹配后，比较计算结果的绝对偏差，结果如表3所示。扫频IRAM采用PSD - SSAP默认的参 数配置。
[0092] 表2全部特征值分析的正确性验证
【权利要求】
1. 一种大规模电力系统小干扰稳定性的特征值分析系统，其特征在于，包括：包括： BPA数据接口模块、用于实现小干扰稳定性特征值分析中的稀疏矩阵相关处理的稀疏矩阵 计算模块、潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块、全部特征值计算 引擎模块、部分特征值计算引擎模块和振荡模式提取和分析模块，其中：BPA数据接口模块 接收电网数据文件，并分别与潮流计算模块和系统建模及其线性化模块相连并提供API ; 稀疏矩阵计算模块分别与潮流计算模块、系统建模及其线性化模块、状态矩阵计算模块和 部分特征值计算引擎模块相连并提供API ;潮流计算模块根据BPA数据接口模块以及稀疏 矩阵计算模块提供的API，实现大规模电力系统的潮流计算并输出潮流计算结果至系统建 模及其线性化模块；系统建模及其线性化模块根据BPA数据接口模块提供的API和潮流计 算模块提供的潮流计算结果，实现大规模电力系统的线性化建模，即生成增广的系统状态 矩阵并分别输出至状态矩阵计算模块和部分特征值计算引擎模块；状态矩阵计算模块根据 稀疏矩阵计算模块提供的API和系统建模及其线性化模块提供的增广的系统状态矩阵，获 得系统状态矩阵并输出至全部特征值计算引擎模块；全部特征值计算引擎模块根据TTQRE 从系统状态矩阵中获得系统状态矩阵的全部特征值以及部分左/右特征向量并输出至振 荡模式提取和分析模块；部分特征值计算引擎模块根据稀疏矩阵计算模块提供的API，以 RCI的方式执行IRAM或者JDM的数值迭代过程，获得系统状态矩阵的部分特征值和部分左 /右特征向量并输出至振荡模式提取和分析模块；振荡模式提取和分析模块根据用户指定 合并来自全部特征值计算引擎模块或者部分特征值计算引擎模块的系统状态矩阵的特征 值和左/右特征向量，实现机电的振荡模式识别和模态分析，并将所得的模式信息和模态 分析结果以Excel表格的方式输出给用户。
2. 根据权利要求1所述的系统，其特征是，所述的增广的系统状态矩阵Jaug满足： /产 - J /其中：JA为动态元件线性化模型拼接而成的分块对角矩阵，JB为动态元件 _JC t 与静态元件的联系矩阵，J。为静态元件与动态元件的联系矩阵，JD为修正的系统导纳矩阵。
3.根据权利要求1所述的系统，其特征是，所述的系统状态矩阵S是指： S - J, 〇
4.根据权利要求1所述的系统，其特征是，所述的全部特征值计算引擎模块内包括 TTQRE单元，该单元与状态矩阵计算模块相连，并根据系统状态矩阵获得系统状态矩阵的全 部特征值和部分左/右特征向量后，输出至振荡模式提取和分析模块。
5.根据权利要求1所述的系统，其特征是，所述的部分特征值计算引擎模块包括：IRAM 单元和JDM单元，其中：IRAM单元与系统建模及其线性化模块相连，根据用户的选择来确定 采用平移-逆变换结合模最大Ritz值选择策略还是Cayley变换结合模最大Ritz值选择策 略的计算方案，并根据增广的系统状态矩阵获得系统状态矩阵的部分特征值和部分左/右 特征向量后，输出至振荡模式提取和分析模块；JDM单元与系统建模及其线性化模块相连， 根据用户的选择来确定采用平移-逆变换结合模最大Ritz值选择策略、Cayley变换结合 模最大Ritz值选择策略、原始矩阵结合实部最大Ritz值选择策略还是原始矩阵结合阻尼 比最小Ritz值选择策略的计算方案，并根据增广的系统状态矩阵获得系统状态矩阵的部 分特征值和部分左/右特征向量后，输出至振荡模式提取和分析模块。
6.根据权利要求1所述的系统，其特征是，所述的振荡模式识别和模态分析是指：根据 每一对左/右特征向量获得相应特征值的相关因子，根据相关因子获得该特征值的机电回 路相关比，若机电回路相关比大于1且该特征值为复数，则判定该特征值为机电振荡模式， 相应的右特征向量为振荡模态；相关因子模值的大小反应了电力系统动态兀件与振荡模式 之间的耦合程度，在耦合程度最大的机组上加装电力系统稳定器，是抑制由该振荡模式诱 发的电力系统低频功率振荡的一贯做法；振荡模态的幅值决定了该振荡模式下机组的摇摆 程度，角度决定了该振荡模式下机组与机组之间的相对摇摆关系，振荡模态是判断振荡模 式对机组的影响程度、判断多个机组是否可以进行动态等值的重要指标。
【文档编号】G06F19/00GK104091092SQ201410366590
【公开日】2014年10月8日 申请日期:2014年7月29日 优先权日:2014年7月29日
【发明者】赵文恺, 严正, 张逸飞, 曹路, 李建华, 周挺辉, 范翔 申请人:上海交通大学, 华东电网有限公司