一种用于矩阵指数的暂态仿真多时间尺度输出方法
【专利摘要】一种用于矩阵指数的暂态仿真多时间尺度输出方法，首先确定待仿真电力系统的结构与参数，得到微分方程形式的暂态模型，初始化并启动仿真程序；在每一个仿真步长内采用矩阵指数方法形成状态矩阵，并确定本步长可用的连续平方环节个数；计算矩阵指数函数的Padé逼近矩阵，并利用这一矩阵的连续平方运算和矩阵向量乘法运算得到离散时间节点之间多个中间节点上的状态向量的近似值，仿真推进一个步长；依此迭代进行，直到仿真结束。本发明实现了矩阵指数积分方法的仿真计算步长与结果输出步长之间的解耦，使矩阵指数积分方法良好的数值精度和刚性处理能力能得到充分的发挥，解决了矩阵指数积分方法在多时间尺度特性的电力系统暂态仿真领域的适用性问题。
【专利说明】一种用于矩阵指数的暂态仿真多时间尺度输出方法
【技术领域】
[0001] 本发明涉及一种暂态仿真的稠密输出方法。特别是涉及一种用于矩阵指数的暂态 仿真多时间尺度输出方法。
【背景技术】
[0002] 电力系统是由发电、输电、配电、用电的大量设备及其相关辅助系统组成的复杂系 统，涉及电磁暂态、机械暂态、电化学过程、热力学过程等多种动态过程，有显著的多时间尺 度特征。
[0003] 在传统电力系统数字仿真领域，针对电力系统动态过程的不同时间尺度分别发展 出了电磁暂态仿真和机电暂态仿真两种建模仿真方法，二者从数学模型到仿真算法都有不 同的特征。其中，电磁暂态仿真主要反映系统中电场与磁场相互耦合影响产生的电气量变 化过程，采用详细的动态建模和微秒级仿真步长，得到从工频到几十kHz频谱范围内的三 相电压电流瞬时值波形；机电暂态仿真则主要关注故障、切机切负荷等大扰动发生后，发电 机转子功角摇摆等机械变化过程，对电力网络及与之相连的机组绕组采用交流稳态建模， 模型较电磁暂态仿真进行了简化，通常采用毫秒级仿真步长，能够捕捉工频以下的系统动 态过程。
[0004] 近年来，我国电力系统实际运行出现的一些仿真场景，为传统电力系统数字仿真 提出了新的问题。输电系统方面，随着特高压技术和直流输电技术的发展与应用，交直流互 联电网规模逐渐扩大；配电系统方面，可再生电源和电力电子装置的大量接入，也使得传统 的无源配网成为有源配网。这些运行工况的复杂暂态特性，客观上需要使用建模更详细的 电磁暂态仿真程序进行研究分析，然而，较大的系统规模和较长的仿真时间无疑给电磁暂 态仿真提出了新的挑战，需要结合问题特性从仿真算法层面提出针对性的改进。
[0005] 电力系统电磁暂态仿真本质上可归结为对动力学系统时域响应的求取，它包括系 统本身的数学模型和与之相适应的数值算法。
[0006] 当前，电力系统电磁暂态仿真基本框架可分为两类，包括节点分析法（Nodal Analysis)和状态变量分析法（State-Variable Analysis)。基于节点分析框架的电磁暂 态仿真方法可概括为先采用某种数值积分方法（通常为梯形积分法）将系统中动态元件的 特性方程差分化，得到等效的计算电导与历史项电流源并联形式的诺顿等效电路，对其求 解即可得到系统中各节点电压的瞬时值。节点分析法广泛应用于EMTP、PSCAD/EMTDC等专 业的电力系统电磁暂态仿真程序中，工程上也称基于节点分析框架的电磁暂态仿真工具为 EMTP类程序。节点分析法的主要优势体现在程序实现难度和仿真计算效率方面，但由于节 点电导方程本身已将数值积分方法与系统模型融为一体，导致EMTP类程序在求解算法选 择方面缺乏灵活性与开放性，同时节点电导方程也不能给出系统本身的特征信息。
[0007] 状态变量分析法属于一般性建模方法（general purpose modeling)，不仅适于 电路与电力系统仿真，同样也适于其它工程领域的动力学系统的建模与仿真。Matlab/ SimPowerSystems软件是状态变量分析框架下暂态仿真程序的典型代表。与节点分析框架 相比，状态方程在模型的计算求解方面具有高度的开放性和灵活性，可方便地选择与问题 相适应的数值积分方法，同时能够提供关于系统各种特征的丰富信息（如系统的特征值）， 进而能够从全局角度了解系统的动态特性，为各种快速、准确、高效的仿真算法的开发与测 试工作提供了便利条件。
[0008] 应用状态变量分析的基础是形成状态方程描述的电力系统暂态仿真模型。改进节 点法Modified Nodal Analysis(MNA)通过KCUKVL等约束关系以及元件伏安特性构造得到 MNA模型，再经过一定的正规化处理（regularization)即可得到状态方程模型；也可采用 一般支路类方法，如Automated State Model Generator (ASMG)方法直接构造得到。基于 这些方法得到的电力系统模型能够很容易的与本发明所提出的算法进行接口。
[0009] 在数值算法方面，传统数值积分方法可分为显式和隐式两类，不同积分方法所具 有的数值稳定性和数值精度各不相同。一般来说，隐式方法处理仿真模型中刚性特征的能 力较强。电力系统由于动态过程时间尺度差异较大，系统模型表现出一定刚性，这使得主流 电磁暂态软件EMTP类程序采用隐式方法以保证数值稳定性。从计算开销方面来看，隐式方 法在每一时步内需求解线性方程组，极大限制了其在大规模系统的应用能力。与之相对的， 传统显式方法无需迭代，在每一时步内的运算量较小，但其有限的数值稳定域使得仿真步 长受到约束，综合来看对刚性系统的仿真性能不佳。
[0010] 矩阵指数积分方法（Exponential Integrator)是近年来从应用数学领域发端的 一种数值积分方法。它使用矩阵指数算子eM精确描述动态系统的线性变化规律，可以准确 求解形如
[0011] X'（t) = Ax (t)+Bu ⑴
[0012] 的线性动力学系统，并具有计算效率高、刚性处理能力强等特点。矩阵指数积分方 法已经在诸如应用物理、化学工程等领域得到一定应用。然而，对于隐式梯形法主导的电力 系统仿真领域，矩阵指数积分方法的单步计算量较大，在相同步长条件下，尽管矩阵指数积 分方法能提供更为精确的仿真结果，但仿真软件使用者更为关心的计算时长问题反而更为 恶化；如果按照得到与梯形法相同误差水平的仿真结果为准则设置仿真步长，计算时间问 题虽然得到解决，但此时仿真步长较大，往往无法提供足够的采样频率，可能遗漏系统动态 过程中的快动态细节。这一两难的处境极大限制了矩阵指数积分方法在电力系统数字仿真 领域的应用，是矩阵指数暂态仿真急需解决的问题。
【发明内容】
[0013] 本发明所要解决的技术问题是，提供一种用于矩阵指数的暂态仿真多时间尺度输 出方法，能够在每一个时步内，获得具有良好精度的依用户设置的若干稠密输出点，使得仿 真可以在采用较大仿真步长计算的同时，得到较小步长的输出结果，从而覆盖待研究电力 系统在多个时间尺度上的动态过程。
[0014] 本发明所采用的技术方案是：一种用于矩阵指数的暂态仿真多时间尺度输出方 法，包括如下步骤：
[0015] 1)根据仿真算例输入文件，确定待仿真电力系统的结构与参数，得到微分方程形 式的暂态模型x' = f (X)，X为状态变量；
[0016] 2)设置仿真结束时间T，设置仿真计算时刻tn = 0,进行仿真初始化计算；
[0017] 3)记tn时刻的仿真步长为h，状态向量为xn，设定需要进行稠密输出的平方环节 的个数Ns，Ns是正整数；
[0018] 4)采用矩阵指数方法形成时步[tn，tn+h]所需的状态矩阵A n ;
[0019] 5)确定[tn，tn+h]时步内，矩阵指数计算所用的连续平方环节个数s，若s〈N s，临时 将s的值赋值给Ns，否则Ns维持不变，记本时步的输出为X，X是2Ns个列向量组成的矩阵；
[0020] 6)计算Pad6逼近矩阵F = rkm(h/2sAn)，'是矩阵指数函数的Pad6逼近函数；
[0021] 7)进行（S-Ns)次F的连续平方，并将结果矩阵赋值给F，将F和Xn的矩阵向量乘 积赋值给X的第一列，作为t n+h/2Ns时刻状态向量的值，如果Ns = 1，跳至步骤11)，否则设 定循环变量i = 1，并进入下一步骤；
[0022] 8) F = F2,将F和Xn的矩阵向量乘积赋值给X的第21列，作为tn+h/2 Ns_i时亥Ij状态 向量的值；
[0023] 9)记X的1至2^1列组成的矩阵为Xi，将F和Xi的矩阵乘积赋值给X的第2kl 至 2i+1-l 列，作为{^+0+1)/2?, tn+(2i+2)/2Nsh，…，tn+(2i+1-l)/2 Nsh}这 2Ll 个时刻的状 态向量的值；
[0024] 10)更新i = i+1，如果i = Ns，结束本时步的仿真计算和输出，将X的值写入输出 文件后进入步骤11)，否则返回步骤8);
[0025] 11)更新tn = tn+h，仿真向前推进一个步长，如果tn小于仿真结束时间T，返回步 骤3)继续进行仿真计算，否则仿真结束。
[0026] 步骤4)所述的采用矩阵指数方法形成时步[tn，tn+h]所需的状态矩阵A n的具体 形式是：
[0027]
【权利要求】
1. 一种用于矩阵指数的暂态仿真多时间尺度输出方法，其特征在于，包括如下步骤： 1) 根据仿真算例输入文件，确定待仿真电力系统的结构与参数，得到微分方程形式的 暂态模型X' = f (X)，X为状态变量； 2) 设置仿真结束时间T，设置仿真计算时刻tn = 0,进行仿真初始化计算； 3) 记tn时刻的仿真步长为h，状态向量为Xn，设定需要进行稠密输出的平方环节的个 数Ns，Ns是正整数； 4) 采用矩阵指数方法形成时步[tn，tn+h]所需的状态矩阵An; 5) 确定[tn，tn+h]时步内，矩阵指数计算所用的连续平方环节个数s，若s〈Ns，临时将s 的值赋值给Ns，否则Ns维持不变，记本时步的输出为X，X是2Ns个列向量组成的矩阵； 6) 计算Pad6逼近矩阵F = rkm(h/2sAn)，rkm是矩阵指数函数的Pad6逼近函数； 7) 进行（S-Ns)次F的连续平方，并将结果矩阵赋值给F，将？和\的矩阵向量乘积赋 值给X的第一列，作为t n+h/2Ns时刻状态向量的值，如果Ns = 1，跳至步骤11)，否则设定循 环变量i = 1，并进入下一步骤； 8. F = F2,将F和Xn的矩阵向量乘积赋值给X的第21列，作为tn+h/2 Ns4时刻状态向量 的值； 9) 记X的1至2^1列组成的矩阵为Xi，将F和Xi的矩阵乘积赋值给X的第2^1至 2 i+1_l 列，作为{tn+(2^1)/2?, tn+(2^2)/2?,…，tn+(2i+1-l)/2Nsh}这 2^1 个时刻的状态 向量的值； 10) 更新i = i+Ι，如果i =Ns，结束本时步的仿真计算和输出，将X的值写入输出文件 后进入步骤11)，否则返回步骤8); 11) 更新tn = tn+h，仿真向前推进一个步长，如果tn小于仿真结束时间T，返回步骤3) 继续进行仿真计算，否则仿真结束。
2. 权利要求1所述的一种用于矩阵指数的暂态仿真多时间尺度输出方法，其特征在 于，步骤4)所述的采用矩阵指数方法形成时步[t n，tn+h]所需的状态矩阵An的具体形式 是：
其中，是非线性函数f (X)在Xn取值时的Jacobian矩阵。
3. 权利要求1所述的一种用于矩阵指数的暂态仿真多时间尺度输出方法，其特征在 于，步骤5)所述的确定[tn，t n+h]时步内矩阵指数计算所用的连续平方环节个数s的具体 方法为： (1) 计算仿真步长h与状态矩阵An的乘积矩阵的列范数，记为Nm ; (2) 如果NmS Θ，其中Θ是通过数值分析得到的阈值，则本时步没有能够进行稠密输 出的平方环节，取s = 0,否则进入下一步骤； ⑶选取小数M和整数E，使得0·5彡M〈l，MX2E = Nm/θ; (4)如果 M = 0· 5,取 s = E-1，否则取 s = E。
【文档编号】G06F17/50GK104376158SQ201410619247
【公开日】2015年2月25日 申请日期:2014年11月5日 优先权日:2014年11月5日
【发明者】王成山, 富晓鹏, 李鹏, 丁承第, 于浩, 于力, 郭晓斌, 许爱东, 董旭柱, 吴争荣, 田兵, 王建邦, 魏文潇 申请人:天津大学, 中国南方电网有限责任公司电网技术研究中心, 南方电网科学研究院有限责任公司