CN104376158B - 一种用于矩阵指数的暂态仿真多时间尺度输出方法 - Google Patents

Abstract

一种用于矩阵指数的暂态仿真多时间尺度输出方法，首先确定待仿真电力系统的结构与参数，得到微分方程形式的暂态模型，初始化并启动仿真程序；在每一个仿真步长内采用矩阵指数方法形成状态矩阵，并确定本步长可用的连续平方环节个数；计算矩阵指数函数的Padé逼近矩阵，并利用这一矩阵的连续平方运算和矩阵向量乘法运算得到离散时间节点之间多个中间节点上的状态向量的近似值，仿真推进一个步长；依此迭代进行，直到仿真结束。本发明实现了矩阵指数积分方法的仿真计算步长与结果输出步长之间的解耦，使矩阵指数积分方法良好的数值精度和刚性处理能力能得到充分的发挥，解决了矩阵指数积分方法在多时间尺度特性的电力系统暂态仿真领域的适用性问题。

Description

一种用于矩阵指数的暂态仿真多时间尺度输出方法

技术领域

本发明涉及一种暂态仿真的稠密输出方法。特别是涉及一种用于矩阵指数的暂态仿真多时间尺度输出方法。

背景技术

电力系统是由发电、输电、配电、用电的大量设备及其相关辅助系统组成的复杂系统，涉及电磁暂态、机械暂态、电化学过程、热力学过程等多种动态过程，有显著的多时间尺度特征。

在传统电力系统数字仿真领域，针对电力系统动态过程的不同时间尺度分别发展出了电磁暂态仿真和机电暂态仿真两种建模仿真方法，二者从数学模型到仿真算法都有不同的特征。其中，电磁暂态仿真主要反映系统中电场与磁场相互耦合影响产生的电气量变化过程，采用详细的动态建模和微秒级仿真步长，得到从工频到几十kHz频谱范围内的三相电压电流瞬时值波形；机电暂态仿真则主要关注故障、切机切负荷等大扰动发生后，发电机转子功角摇摆等机械变化过程，对电力网络及与之相连的机组绕组采用交流稳态建模，模型较电磁暂态仿真进行了简化，通常采用毫秒级仿真步长，能够捕捉工频以下的系统动态过程。

近年来，我国电力系统实际运行出现的一些仿真场景，为传统电力系统数字仿真提出了新的问题。输电系统方面，随着特高压技术和直流输电技术的发展与应用，交直流互联电网规模逐渐扩大；配电系统方面，可再生电源和电力电子装置的大量接入，也使得传统的无源配网成为有源配网。这些运行工况的复杂暂态特性，客观上需要使用建模更详细的电磁暂态仿真程序进行研究分析，然而，较大的系统规模和较长的仿真时间无疑给电磁暂态仿真提出了新的挑战，需要结合问题特性从仿真算法层面提出针对性的改进。

电力系统电磁暂态仿真本质上可归结为对动力学系统时域响应的求取，它包括系统本身的数学模型和与之相适应的数值算法。

当前，电力系统电磁暂态仿真基本框架可分为两类，包括节点分析法(NodalAnalysis)和状态变量分析法(State-Variable Analysis)。基于节点分析框架的电磁暂态仿真方法可概括为先采用某种数值积分方法(通常为梯形积分法)将系统中动态元件的特性方程差分化，得到等效的计算电导与历史项电流源并联形式的诺顿等效电路，对其求解即可得到系统中各节点电压的瞬时值。节点分析法广泛应用于EMTP、PSCAD/EMTDC等专业的电力系统电磁暂态仿真程序中，工程上也称基于节点分析框架的电磁暂态仿真工具为EMTP类程序。节点分析法的主要优势体现在程序实现难度和仿真计算效率方面，但由于节点电导方程本身已将数值积分方法与系统模型融为一体，导致EMTP类程序在求解算法选择方面缺乏灵活性与开放性，同时节点电导方程也不能给出系统本身的特征信息。

状态变量分析法属于一般性建模方法(general purpose modeling)，不仅适于电路与电力系统仿真，同样也适于其它工程领域的动力学系统的建模与仿真。Matlab/SimPowerSystems软件是状态变量分析框架下暂态仿真程序的典型代表。与节点分析框架相比，状态方程在模型的计算求解方面具有高度的开放性和灵活性，可方便地选择与问题相适应的数值积分方法，同时能够提供关于系统各种特征的丰富信息(如系统的特征值)，进而能够从全局角度了解系统的动态特性，为各种快速、准确、高效的仿真算法的开发与测试工作提供了便利条件。

应用状态变量分析的基础是形成状态方程描述的电力系统暂态仿真模型。改进节点法Modified Nodal Analysis(MNA)通过KCL、KVL等约束关系以及元件伏安特性构造得到MNA模型，再经过一定的正规化处理(regularization)即可得到状态方程模型；也可采用一般支路类方法，如Automated State Model Generator(ASMG)方法直接构造得到。基于这些方法得到的电力系统模型能够很容易的与本发明所提出的算法进行接口。

在数值算法方面，传统数值积分方法可分为显式和隐式两类，不同积分方法所具有的数值稳定性和数值精度各不相同。一般来说，隐式方法处理仿真模型中刚性特征的能力较强。电力系统由于动态过程时间尺度差异较大，系统模型表现出一定刚性，这使得主流电磁暂态软件EMTP类程序采用隐式方法以保证数值稳定性。从计算开销方面来看，隐式方法在每一时步内需求解线性方程组，极大限制了其在大规模系统的应用能力。与之相对的，传统显式方法无需迭代，在每一时步内的运算量较小，但其有限的数值稳定域使得仿真步长受到约束，综合来看对刚性系统的仿真性能不佳。

矩阵指数积分方法(Exponential Integrator)是近年来从应用数学领域发端的一种数值积分方法。它使用矩阵指数算子e^hA精确描述动态系统的线性变化规律，可以准确求解形如

x'(t)＝Ax(t)+Bu(t)

的线性动力学系统，并具有计算效率高、刚性处理能力强等特点。矩阵指数积分方法已经在诸如应用物理、化学工程等领域得到一定应用。然而，对于隐式梯形法主导的电力系统仿真领域，矩阵指数积分方法的单步计算量较大，在相同步长条件下，尽管矩阵指数积分方法能提供更为精确的仿真结果，但仿真软件使用者更为关心的计算时长问题反而更为恶化；如果按照得到与梯形法相同误差水平的仿真结果为准则设置仿真步长，计算时间问题虽然得到解决，但此时仿真步长较大，往往无法提供足够的采样频率，可能遗漏系统动态过程中的快动态细节。这一两难的处境极大限制了矩阵指数积分方法在电力系统数字仿真领域的应用，是矩阵指数暂态仿真急需解决的问题。

发明内容

本发明所要解决的技术问题是，提供一种用于矩阵指数的暂态仿真多时间尺度输出方法，能够在每一个时步内，获得具有良好精度的依用户设置的若干稠密输出点，使得仿真可以在采用较大仿真步长计算的同时，得到较小步长的输出结果，从而覆盖待研究电力系统在多个时间尺度上的动态过程。

本发明所采用的技术方案是：一种用于矩阵指数的暂态仿真多时间尺度输出方法，包括如下步骤：

1)根据仿真算例输入文件，确定待仿真电力系统的结构与参数，得到微分方程形式的暂态模型x'＝f(x)，x为状态变量；

2)设置仿真结束时间T，设置仿真计算时刻t_n＝0，进行仿真初始化计算；

3)记t_n时刻的仿真步长为h，状态向量为x_n，设定需要进行稠密输出的平方环节的个数N_S，N_S是正整数；

4)采用矩阵指数方法形成时步[t_n,t_n+h]所需的状态矩阵A_n；

5)确定[t_n,t_n+h]时步内，矩阵指数计算所用的连续平方环节个数s，若s<N_S，临时将s的值赋值给N_S，否则N_S维持不变，记本时步的输出为X，X是2^Ns个列向量组成的矩阵；

6)计算Padé逼近矩阵F＝r_km(h/2^sA_n)，r_km是矩阵指数函数的Padé逼近函数；

7)进行(s-N_S)次F的连续平方，并将结果矩阵赋值给F，将F和x_n的矩阵向量乘积赋值给X的第一列，作为t_n+h/2^Ns时刻状态向量的值，如果N_S＝1，跳至步骤11)，否则设定循环变量i＝1，并进入下一步骤；

8)F＝F²，将F和x_n的矩阵向量乘积赋值给X的第2ⁱ列，作为t_n+h/2^Ns-i时刻状态向量的值；

9)记X的1至2ⁱ-1列组成的矩阵为X_i，将F和X_i的矩阵乘积赋值给X的第2ⁱ+1至2ⁱ⁺¹-1列，作为{t_n+(2ⁱ+1)/2^Nsh,t_n+(2ⁱ+2)/2^Nsh,…,t_n+(2ⁱ⁺¹-1)/2^Nsh}这2ⁱ-1个时刻的状态向量的值；

10)更新i＝i+1，如果i＝N_S，结束本时步的仿真计算和输出，将X的值写入输出文件后进入步骤11)，否则返回步骤8)；

11)更新t_n＝t_n+h，仿真向前推进一个步长，如果t_n小于仿真结束时间T，返回步骤3)继续进行仿真计算，否则仿真结束。

步骤4)所述的采用矩阵指数方法形成时步[t_n,t_n+h]所需的状态矩阵A_n的具体形式是：

其中，是非线性函数f(x)在x_n取值时的Jacobian矩阵。

步骤5)所述的确定[t_n,t_n+h]时步内矩阵指数计算所用的连续平方环节个数s的具体方法为：

(1)计算仿真步长h与状态矩阵A_n的乘积矩阵的列范数，记为Nm；

(2)如果Nm≤θ，其中θ是通过数值分析得到的阈值，则本时步没有能够进行稠密输出的平方环节，取s＝0，否则进入下一步骤；

(3)选取小数M和整数E，使得0.5≤M<1，M×2^E＝Nm/θ；

(4)如果M＝0.5，取s＝E-1，否则取s＝E。

本发明的一种用于矩阵指数的暂态仿真多时间尺度输出方法，在一定程度上实现了矩阵指数积分方法的仿真计算步长与结果输出步长之间的解耦，充分利用了Scaling&Squaring方法求解矩阵指数函数计算过程中的中间结果，得到精度不差于原始输出的稠密输出点。本发明能够在每一个时步内，获得具有良好精度的依用户设置的若干稠密输出点，使得仿真可以在采用较大仿真步长计算的同时，得到较小步长的输出结果，从而覆盖待研究电力系统在多个时间尺度上的动态过程。在本发明的工作方式下，大的计算步长需求和小的输出步长需求之间的矛盾得以解决，矩阵指数积分方法良好的数值精度和刚性处理能力得到了充分的发挥，解决了矩阵指数积分方法在多时间尺度特性的电力系统暂态仿真领域的适用性问题。

附图说明

图1是本发明的一种用于矩阵指数的暂态仿真多时间尺度输出方法的流程图；

图2是本发明实施例中同步发电机组的结构示意图；

图3是图2的等值电路图；

图4是同步发电机组定子q轴上的电流i_qs曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，黑色虚线表示传统暂态仿真方法得到的仿真曲线；

图5是同步发电机组定子d轴上的电流i_ds曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，黑色虚线表示传统暂态仿真方法得到的仿真曲线；

图6是同步发电机组励磁绕组上的电流i_fd曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，黑色虚线表示传统暂态仿真方法得到的仿真曲线；

图7是同步发电机组定子q轴上的电压v_qs曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，传统暂态仿真方法得到的仿真曲线用黑色虚线表示；

图8是同步发电机组定子d轴上的电压v_ds曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，黑色虚线表示传统暂态仿真方法得到的仿真曲线；

图9是同步发电机组转子转速ω_r曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，黑色虚线表示传统暂态仿真方法得到的仿真曲线；

图10是同步发电机组功角δ曲线图；

图中，黑色实线表示本发明所述方法得到的仿真曲线，黑色虚线表示传统暂态仿真方法得到的仿真曲线。

具体实施方式

下面结合实施例和附图对本发明的一种用于矩阵指数的暂态仿真多时间尺度输出方法做出详细说明。

本发明的一种用于矩阵指数的暂态仿真多时间尺度输出方法，算法流程图如图1所示。采用同步发电机出口三相短路这一工况作为实施例进行说明，本发明包括如下步骤：

1)根据仿真算例输入文件，确定待仿真电力系统的结构与参数，得到微分方程形式的暂态模型x'＝f(x)，x为状态变量；

对于本实施例，同步发电机组的结构与等效电路如图2、图3所示，其参数为：

2)设置仿真结束时间T，设置仿真计算时刻t_n＝0，进行仿真初始化计算；

对于本实施例，设定仿真时长T＝3.0s，仿真从潮流解启动，状态初值为x_n＝

[-5.065205167639663e+03

1.555227402985741e+04

-3.798689647800140e+03

2.439057797934730e+04

1.704945533934731e+04

11.780972450961723e+00

-1.256941276453054e+00]

3)记t_n时刻的仿真步长为h，状态向量为x_n，设定需要进行稠密输出的平方环节的个数N_S，N_S是正整数；

对于本实施例，采用定步长h＝0.5/60≈8.3333ms，即每个工频周波计算两步，属于机电暂态仿真的步长范围；设置N_S＝5，即每步采用本发明所述算法计算5个稠密输出，每周波总计输出10个仿真结果，可以捕捉到各绕组上电流的工频分量；

4)采用矩阵指数方法形成时步[t_n,t_n+h]所需的状态矩阵A_n，所述A_n的具体形式是：

其中，是非线性函数f(x)取值为x_n时的Jacobian矩阵；

5)确定[t_n,t_n+h]时步内，矩阵指数计算所用的连续平方环节个数s，若s<N_S，临时将s的值赋值给N_S，否则N_S维持不变，记本时步的输出为X，X是2^Ns个列向量组成的矩阵；

确定平方环节个数s的具体方法为：

(1)计算仿真步长h与状态矩阵A_n的乘积矩阵的列范数，记为Nm；

(2)如果Nm≤θ，其中θ是通过数值分析得到的阈值，则本时步没有能够进行稠密输出的平方环节，取s＝0，否则进入下一步骤；

(3)选取小数M和整数E，使得0.5≤M<1，M×2^E＝Nm/θ；

(4)如果M＝0.5，取s＝E-1，否则取s＝E；

6)计算Padé逼近矩阵F＝r_km(h/2^sA_n)，r_km是矩阵指数函数的Padé逼近函数；

7)进行(s-N_S)次F的连续平方，并将结果矩阵赋值给F，将F和x_n的矩阵向量乘积赋值给X的第一列，作为t_n+h/2^Ns时刻状态向量的值，如果N_S＝1，跳至步骤11)，否则设定循环变量i＝1，并进入下一步骤；

8)F＝F²，将F和x_n的矩阵向量乘积赋值给X的第2ⁱ列，作为t_n+h/2^Ns-i时刻状态向量的值；

9)记X的1至2ⁱ-1列组成的矩阵为X_i，将F和X_i的矩阵乘积赋值给X的第2ⁱ+1至2ⁱ⁺¹-1列，作为{t_n+(2ⁱ+1)/2^Nsh,t_n+(2ⁱ+2)/2^Nsh,…,t_n+(2ⁱ⁺¹-1)/2^Nsh}这2ⁱ-1个时刻的状态向量的值；

10)更新i＝i+1，如果i＝N_S，结束本时步的仿真计算和输出，将X的值写入输出文件后进入步骤11)，否则返回步骤8)；

11)更新t_n＝t_n+h，仿真向前推进一个步长，如果t_n小于仿真结束时间T，返回步骤3)继续进行仿真计算，否则仿真结束。

本发明的一种用于矩阵指数的暂态仿真多时间尺度输出方法，执行仿真计算的计算机硬件环境为Intel Core2Q84002.66GHz CPU，内存容量2GB；软件环境为Windows 7操作系统。

将本发明所提出的一种用于矩阵指数的暂态仿真多时间尺度输出方法所得的仿真结果(dense)，与传统机电暂态仿真方法所得的仿真结果(traditional)进行比较。

图4至图6比较了三个绕组上电流的仿真波形。图中可以看出，传统暂态仿真方法只能得到电流波形中的低频慢动态，而实际的电流波形则包含有指数衰减的工频分量动态，这一时间尺度的动态过程可以被本发明所述算法计算得到。

图7、图8比较了定子绕组dq轴上电压的仿真波形。图中可以看出，对于不包含有快动态细节的动态过程，本发明所述算法与传统暂态仿真算法可以得到很好的一致性。

图9、图10比较了同步电机转子的转速和功角的仿真波形。图中本发明所述算法与传统暂态仿真算法得到了很好的一致性。特别的，对于图9，故障发生初期绕组电流上的工频量造成转速小幅度的波动这一现象也通过本发明所述算法反映了出来。这充分验证了本发明所述算法在捕获不同时间尺度动态细节上的显著优势。

以上算例测试结果证明，本发明提出的一种用于矩阵指数的暂态仿真多时间尺度输出方法具有良好的可行性与适用性，为解决具有多时间尺度特性的电力系统暂态仿真提供了一种很好的解决思路。

Claims (2)

1.一种用于矩阵指数的暂态仿真多时间尺度输出方法，其特征在于，包括如下步骤：

1)根据仿真算例输入文件，确定待仿真电力系统的结构与参数，得到微分方程形式的暂态模型x'＝f(x)，x为状态变量；

2)设置仿真结束时间T，设置仿真计算时刻t_n＝0，进行仿真初始化计算；

3)记t_n时刻的仿真步长为h，状态向量为x_n，设定需要进行稠密输出的平方环节的个数N_S，N_S是正整数；

4)采用矩阵指数方法形成时步[t_n,t_n+h]所需的状态矩阵A_n；

5)确定[t_n,t_n+h]时步内，矩阵指数计算所用的连续平方环节个数s，若s<N_S，临时将s的值赋值给N_S，否则N_S维持不变，记本时步的输出为X，X是2^Ns个列向量组成的矩阵；

所述的确定[t_n,t_n+h]时步内矩阵指数计算所用的连续平方环节个数s的具体方法为：

(1)计算仿真步长h与状态矩阵A_n的乘积矩阵的列范数，记为Nm；

(2)如果Nm≤θ，其中θ是通过数值分析得到的阈值，则本时步没有能够进行稠密输出的平方环节，取s＝0，否则进入下一步骤；

(3)选取小数M和整数E，使得0.5≤M<1，M×2^E＝Nm/θ；

(4)如果M＝0.5，取s＝E-1，否则取s＝E；

6)计算Padé逼近矩阵F＝r_km(h/2^sA_n)，r_km是矩阵指数函数的Padé逼近函数；

7)进行(s-N_S)次F的连续平方，并将结果矩阵赋值给F，将F和x_n的矩阵向量乘积赋值给X的第一列，作为t_n+h/2^Ns时刻状态向量的值，如果N_S＝1，跳至步骤11)，否则设定循环变量i＝1，并进入下一步骤；

8)F＝F²，将F和x_n的矩阵向量乘积赋值给X的第2ⁱ列，作为t_n+h/2^Ns-i时刻状态向量的值；

9)记X的1至2ⁱ-1列组成的矩阵为X_i，将F和X_i的矩阵乘积赋值给X的第2ⁱ+1至2ⁱ⁺¹-1列，作为{t_n+(2ⁱ+1)/2^Nsh,t_n+(2ⁱ+2)/2^Nsh,…,t_n+(2ⁱ⁺¹-1)/2^Nsh}这2ⁱ-1个时刻的状态向量的值；

10)更新i＝i+1，如果i＝N_S，结束本时步的仿真计算和输出，将X的值写入输出文件后进入步骤11)，否则返回步骤8)；

11)更新t_n＝t_n+h，仿真向前推进一个步长，如果t_n小于仿真结束时间T，返回步骤3)继续进行仿真计算，否则仿真结束。

2.权利要求1所述的一种用于矩阵指数的暂态仿真多时间尺度输出方法，其特征在于，步骤4)所述的采用矩阵指数方法形成时步[t_n,t_n+h]所需的状态矩阵A_n的具体形式是：

$A_n=\left[\begin{array}{cc}\frac{\&part;f}{\&part;x}\left(x_n\right)& f\left(x_n\right)-\frac{\&part;f}{\&part;x}\left(x_n\right)x_n\\ 0& 0\end{array}\right]$

其中，是非线性函数f(x)在x_n取值时的Jacobian矩阵。