CN103646152B - 一种基于矩阵指数的电力系统电磁暂态仿真方法 - Google Patents

Abstract

一种基于矩阵指数的电力系统电磁暂态仿真方法，首先在状态分析框架下建立待研究系统的电磁暂态仿真模型：设定仿真步长Δt及相关仿真参数后，启动仿真程序；在每一个仿真步长内，利用当前时刻激励源g(t)、各阶导数信息、状态矩阵A组成增广的状态矩阵采用Scaling and Squaring算法求解矩阵指数与状态变量向量x(t)与p维单位向量组成的增广状态向量进行矩阵向量乘法，结果作为当前时刻的状态变量，由输出方程得到输出向量y(t)，写入输出文件，仿真推进一个步长；并在下一个仿真步长内，首先需判断是否有短路、开路和开关动作。本发明能够实现结构复杂且高度不对称的、模型具有刚性特性的电力系统电磁暂态建模仿真。

Description

一种基于矩阵指数的电力系统电磁暂态仿真方法

技术领域

本发明涉及一种用于电力系统电磁暂态仿真建模方法。特别是涉及一种基于矩阵指数的电力系统电磁暂态仿真方法

背景技术

电力系统电磁暂态仿真主要反映系统中电场与磁场的相互影响产生的电气量的变化过程，可得到几十kHz到工频范围内的三相电压电流瞬时值波形。为了准确获取系统动态过程中的高频特性，电磁暂态的仿真步长通常是微秒级。随着直流输电技术的发展与广泛应用，交直流系统互联给电磁暂态仿真提出了新的挑战。近年来我国电力系统实际运行中出现的一些新的仿真场景，如直流输电送端孤岛运行的电磁暂态仿真，系统规模大，仿真时间长，使用传统电磁仿真算法会消耗较多计算资源和计算时间，需要结合模型特性从算法层面提出针对性的改进。

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

当前，电力系统电磁暂态仿真基本框架可分为两类，包括节点分析法（nodalanalysis）和状态变量分析法（state-variable analysis）。基于节点分析框架的电磁暂态仿真方法可概括为先采用某种数值积分方法（通常为梯形积分法）将系统中动态元件的特性方程差分化，得到等效的计算电导与历史项电流源并联形式的诺顿等效电路，此时联立整个电气系统的元件特性方程形成节点电导矩阵，如式（1）所示，对其求解即可得到系统中各节点电压的瞬时值。

Gu=i （1）

式（1）所示的节点电导矩阵为线性方程组，可使用各种成熟的线性稀疏矩阵算法库进行求解。节点分析法广泛应用于EMTP、PSCAD/EMTDC等专业的电力系统电磁暂态仿真程序中，工程上也称基于节点分析框架的电磁暂态仿真工具为EMTP类程序。节点分析法的主要优势体现在程序实现难度和仿真计算效率方面，但由于式（1）的节点电导方程本身已将数值积分方法与系统模型融为一体，导致EMTP类程序在求解算法选择方面缺乏灵活性与开放性，同时式（1）已不能给出系统本身的特征信息。

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

应用状态变量分析的基础是形成式（2）所示标准形式的状态-输出方程，此时系统中的电源作为输入u。

$\begin{array}{c}\stackrel{\·}{x}=A^{\′}x+B^{\′}u\\ y=C^{\′}x+D^{\′}u\end{array}---\left(2\right)$

在电力系统仿真领域，式（2）可以由改进节点法Modified Nodal Analysis（MNA）通过KCL、KVL等约束关系以及元件伏安特性进行构造而得到MNA模型，再经过一定的正规化处理（regularization）转化而来。MNA模型是形如式（3）的状态-输出方程。

$\begin{array}{c}C\stackrel{\·}{x}+\mathrm{Gx}=\mathrm{Bu}\\ y-\mathrm{Lx}\end{array}---\left(3\right)$

也可以采用一般支路等方法，如Automated State Model Generator（ASMG）方法直接构造得到。基于这些方法得到的电力系统模型能够很容易的与本发明所采用的状态变量分析框架下的电磁暂态仿真程序进行接口。

在数值算法方面，传统数值积分方法可分为显式和隐式两类，不同积分方法所具有的数值稳定性和数值精度各不相同。一般来说，隐式方法处理仿真模型中刚性特征的能力较强。电力系统由于动态过程时间尺度差异较大，系统模型表现出一定刚性，这使得主流电磁暂态软件EMTP类程序采用隐式方法以保证数值稳定性。从计算开销方面来看，隐式方法在每一时步内需求解线性方程组，极大限制了其仿真在大规模系统的能力。与之相对的，传统显式方法无需迭代，在每一时步内的运算量较小，但其有限的数值稳定域使得仿真步长受到约束，综合来看对刚性系统的仿真性能不佳。对于现代电力系统来说，系统中既存在微秒级的电力电子开关动态过程，又存在同步机组的励磁、调速等秒级的机电动态过程，时间尺度差异极大，系统刚性特征十分显著。充分利用状态方程框架在数值算法选择方面的灵活性，结合电力系统电磁暂态仿真的应用场景与特殊需求，发展合适的数值积分方法，是提高电力系统电磁暂态仿真的计算性能和应用前景的关键。

随着电力电子设备如FACTS、HVDC，以及分布式电源并网系统等在现代电力系统内的广泛应用，电磁暂态和机电暂态之间的界限逐渐模糊。由于电磁暂态仿真采用了系统级的详细模型，能准确反映电力电子开关动作，能与现场运行的实际情况较好吻合，其在实际工程应用中显得尤为重要。在一些特定的仿真场景，如直流输电系统对次同步谐振的影响，或由电力电子装置驱动的大惯性机组的运行情况等，通常要求进行数十秒的电磁暂态仿真模拟。这些场景客观上要求电磁暂态仿真采取尽可能大的仿真步长，以缩短仿真计算时间。电力电子设备在电力系统中的广泛应用，使得电力系统动态过程的多时间尺度特性更为明显，仿真模型刚性更强。

发明内容

本发明所要解决的技术问题是，提供一种适合处理电磁暂态模型刚性特性，保证数值稳定性，同时具有较高精度，并可增大仿真步长的数值积分的基于矩阵指数的电力系统电磁暂态仿真方法。

本发明所采用的技术方案是：一种基于矩阵指数的电力系统电磁暂态仿真方法，包括如下步骤：

1）在状态分析框架下建立待研究电力系统的电磁暂态仿真模型，形式为：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+g\left(t\right)$

y(t)=Lx(t)

其中，A是状态矩阵，g(t)是激励源，L是输出矩阵，x(t)是状态变量，y(t)是输出变量；

2）启动仿真程序，设定仿真时间T_end，仿真步长Δt，多项式近似系数p，设定当前时刻为仿真起始时刻t_n=t₀，依照仿真需要设置仿真初值x_n=x₀，其中的下标_n表示t_n时刻的相关变量；

3）利用当前时刻激励源及当前时刻激励源各阶导数信息，组成增广的状态矩阵

4）采用Scaling and Squaring算法求解计算矩阵指数

5）状态变量向量与p维单位向量 $e_p=\left[\begin{array}{c}0\\ .\\ .\\ .\\ 1\end{array}\right]$ 组成增广状态向量，与步骤4）得到的矩阵指数进行矩阵向量乘法，得到 $e^{\mathrm{\Δt}\hat{A}}\left[\begin{array}{c}{x}_n\\ e_p\end{array}\right],$ 取结果向量的前n行赋值给x_n；

6）由y_n=Lx_n得到t_n时刻的输出向量并写入输出数据文件，更新当前时刻为下一时刻t_n=t_n+Δt，仿真向前推进一个步长；

7）比较当前时刻t_n与仿真时间T_end，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则继续到步骤8）；

8）判断当前时刻是否有短路和开路系统故障，或开关动作导致系统模型发生改变的事件，若有，重新建立系统仿真模型，更新状态矩阵A、输出矩阵L；若没有，则跳过更新；

9）返回到步骤3）继续进行下一时刻计算，依此迭代进行，直到仿真结束。

步骤3）所述的利用当前时刻激励源及当前时刻激励源各阶导数信息，组成增广的状态矩阵具体方法如下：

计算当前时刻的激励源g(t_n)，及当前时刻激励源g(t_n)的1至p-1阶导数，如果g(t)的解析形式未知，则用数值微分近似t_n时刻g(t)的各阶导数，依照W=[g^(p-1)(t₀),…,g’’(t₀),g’(t₀),g(t₀)]的格式将这些向量组成矩阵W，按照 $\hat{A}=\left[\begin{array}{cc}A& W\\ 0& J\end{array}\right]\∈C^{(n+p)\×(n+p)}$ 的格式更新矩阵，其中矩阵J的格式是固定的 $J=\left[\begin{array}{cc}0& I_{p-1}\\ 0& 0\end{array}\right]\∈C^{p\×p},$ 其中I是单位矩阵，n是系统维数，n为大于0的整数。

步骤4）所述的采用Scaling and Squaring算法求解计算矩阵指数包括如下步骤：

（1）读入算法参数θ，b₁，b₂，…b₁₃，其中θ是scaling因子选择阈值，b_i是Padé逼近系数；

（2）求解scaling因数s，使得s是满足||A/2^s||₁<θ的最小正整数，其中||·||₁是指矩阵的1范数；

（3）对矩阵A进行scaling运算，即A=A/2^s；

（4）求解矩阵e^A的[13/13]阶Pade有理近似r_A，其过程为：

矩阵运算A₂=A²，A₄=A₂A₂，A₆=A₂A₄；

矩阵运算U=A[A₆(b₁₃A₆+b₁₁A₄+b₉A₂)+b₇A₆+b₅A₄+b₃A₂+b₁I]；

矩阵运算V=A₆(b₁₂A₆+b₁₀A₄+b₈A₂)+b₆A₆+b₄A₄+b₂A₂+b₀I；

由矩阵线性方程(-U+V)r_A=U+V求解r_A；

（5）通过对r_A连续平方的方式，计算得到，将这一结果作为scaling操作前原矩阵A的矩阵指数e^A的近似。

本发明的一种基于矩阵指数的电力系统电磁暂态仿真方法，在状态分析框架下，基于矩阵指数运算实现了常微分方程初值问题的求解，能够实现结构复杂且高度不对称的、模型具有刚性特性的电力系统电磁暂态建模仿真，并在大的仿真步长下保证仿真的数值稳定性和数值精度。本发明的方法不仅适于传统电力系统电磁暂态仿真，更由于其具有的良好的刚性处理能力、大步长仿真能力，特别适用于多时间尺度特征显著的现代电力系统的电磁暂态仿真，所以本发明的方法具有很好的适用性。

附图说明

图1是IEEE123节点基准电网算例结构图；

图2是IEEE123节点基准电网算例特征值分布图；

图3是本发明的一种基于矩阵指数的电力系统电磁暂态仿真方法的流程图；

图4是本发明中求解矩阵指数e^A所采用的Scaling and Squaring算法流程图；

图5是IEEE123节点算例中配电变压器上流过的C相电流；

图6是IEEE123节点算例中配电变压器上流过的C相电流在区间[0.0332，0.03355]上的放大图，图中可清楚看出本发明提出的矩阵指数方法与隐式梯形法之间的数值精度比较；

图7是IEEE123节点算例中配电变压器上流过的C相电流在区间[0.01664，0.01666]上的放大图，图中可清楚看出本发明提出的矩阵指数方法与4阶显式龙格库塔法之间的数值精度比较；

图8是IEEE123节点算例中配电变压器上流过的C相电流，采用不同数值积分方法时的绝对误差比较。

具体实施方式

在电力系统电磁暂态仿真中，一次系统的主要设备如输电线路、电力变压器等元件可由其相应的等效电路表示，可以采用MNA或ASMG等自动建模方法（本发明采用ASMG方法）建立形如式（4）的电力系统电磁暂态模型。

$\begin{array}{c}C\stackrel{\&CenterDot;}{x}\left(t\right)+\mathrm{Gx}\left(t\right)=\mathrm{Bu}\left(t\right)\\ y\left(t\right)=\mathrm{Lx}\left(t\right)\end{array}---\left(4\right)$

其中，向量x(t)包含节点电压和一部分支路电流，向量u(t)包含网络中的独立电压源和独立电流源，矩阵C是电感和电容矩阵，G是电导矩阵，B是独立电源的输入矩阵。在采用MNA时，矩阵C可能是奇异的，需要采用额外的正规化方法进行处理。由此即可得到A=-C^{- 1}G，g(t)=C^-1Bu(t)。此时系统模型成为式（5）形式。

$\begin{array}{c}\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+g\left(t\right)\\ y\left(t\right)=\mathrm{Lx}\left(t\right)\end{array}---\left(5\right)$

电磁仿真程序中包含的控制系统，如同步电机的励磁、调速、PSS，通常采用控制框图或传递函数等形式建模。它们可以很容易地转化为形如式（5）的状态方程模型。这样就可以由描述一次系统元件连接关系与参数、控制系统结构和参数的输入数据文件，得到可用于状态变量框架下的电磁暂态仿真系统模型。

当系统发生短路、断线故障等导致系统模型发生改变的事件时，应由MNA或ASMG方法重新形成系统的状态模型。

一般来说，根据经典线性系统理论，形如（5）所示的一般性动态系统从给定初始状态（6）开始的运动轨迹，可以使用矩阵指数算子和卷积运算解析地表示为式（7）形式的Volterra积分方程。

x(t₀)=x₀ （6）

$x\left(t\right)=e^{(t-t_0)A}{x}_0+\underset{t_0}{\overset{t}{\&Integral;}}{e}^{(t-\mathrm{\&tau;})A}g\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;}---\left(7\right)$

其中，原微分方程组中线性部分，即矩阵A的作用体现在矩阵指数e^(t-t0)A一项中，即通常所称的状态转移矩阵。矩阵A的实质含义是状态量之间的相互耦合影响关系，包含了电力系统模型的主要信息。对这一部分的准确求解确保了数值积分算法的数值精度。激励项的作用则体现在卷积积分中，这一项表示状态量受电压源、电流源的直接影响。出于计算复杂度考虑，通过引入族函数定义对卷积项的求解进行转化。族函数的成员定义为：

族函数成员之间满足以下递归关系：

基函数为

电力系统仿真中的电压源、电流源通常满足一定的光滑性条件，可以采用p-1次多项式近似激励源g(t)，则式（7）重新组织为至这p+1个族函数的加权求和。本发明中选取多项式为p-1阶截断泰勒展开式，即式（10）。将其带入到式（7），记h=t-t₀为第n时步仿真步长，x_n和x_n+1分别是离散时间点t_n和t_n+1上真值x(t_n)和x(t_n+1)的数值近似，则结合族函数定义即可得到式（11）。

$g\left(t\right)\&ap;\underset{k=0}{\overset{p-1}{\mathrm{\&Sigma;}}}\frac{g^{\left(k\right)}\left(t_0\right)}{k!}{(t-t_0)}^k---\left(10\right)$

此时从数学角度来看，已将常微分方程的初值问题求解转化为矩阵指数函数的求解。通过调用成熟、高效的矩阵指数数值算法，（本发明采用Scaling and Squaring方法），每一时步求解p+1个族函数与向量的乘积，即可对该系统进行电磁暂态仿真计算。

为了减少计算量，本发明进一步设计了如下算法：

令W=[g^(p-1)(t₀),…,g’’(t₀),g’(t₀),g(t₀)]∈C^n×p，并由分块形式组成矩阵

$\hat{A}=\left[\begin{array}{cc}A& W\\ 0& J\end{array}\right]\&Element;C^{(n+p)\&times;(n+p)}---\left(12\right)$

$J=\left[\begin{array}{cc}0& I_{p-1}\\ 0& 0\end{array}\right]\&Element;C^{p\&times;p}---\left(13\right)$

那么式（14）所得结果与式（11）等价。利用这一方法，可将式（11）中的p+1个n维族函数求解合并为一个维数为(n+p)的族函数求解。

$x_{n+1}=\left[\begin{array}{cc}{I}_n& 0\end{array}\right]e^{h_n\hat{A}}\left[\begin{array}{c}{x}_n\\ e_p\end{array}\right]---\left(14\right)$

本发明方法通过矩阵指数算子的引入使积分方法具有A稳定性，克服了系统刚性对数值稳定性的限制；同时，通过对系统的线性动态过程进行精确求解，避免了传统数值积分方法所遇到的局部截断误差问题，使算法具有很高的数值精度，这使得仿真可以采用大步长进行计算。本发明所提出算法的性质可以充分满足当前电磁暂态仿真的新需求，且算法的形式可以与MNA方法、ASMG方法等电力系统自动建模方法方便地进行接口，为进一步进行电磁暂态仿真算法的开发奠定了良好的基础。

如图3所示，本发明的一种基于矩阵指数的电力系统电磁暂态仿真方法，具体包括如下步骤：

1）在状态分析框架下建立待研究电力系统的电磁暂态仿真模型，形式为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+g\left(t\right)$

y(t)=Lx(t)

其中，A是状态矩阵，g(t)是激励源，L是输出矩阵，上标^T是转置运算，x(t)是状态变量，y(t)是输出变量；

2）启动仿真程序，设定仿真时间T_end，仿真步长Δt，多项式近似系数p，设定当前时刻为仿真起始时刻t_n=t₀，依照仿真需要设置仿真初值x_n=x₀，其中的下标_n表示t_n时刻的相关变量；

3）利用当前时刻激励源及当前时刻激励源各阶导数信息，组成增广的状态矩阵具体是指计算当前时刻的激励源g(t_n)，及当前时刻激励源g(t_n)的1至p-1阶导数，如果g(t)的解析形式未知，则用数值微分近似t_n时刻g(t)的各阶导数，依照W=[g^(p-1)(t₀),…,g’’(t₀),g’(t₀),g(t₀)]的格式将这些向量组成矩阵W，按照 $\hat{A}=\left[\begin{array}{cc}A& W\\ 0& J\end{array}\right]\&Element;C^{(n+p)\&times;(n+p)}$ 的格式更新矩阵，其中矩阵J的格式是固定的 $J=\left[\begin{array}{cc}0& I_{p-1}\\ 0& 0\end{array}\right]\&Element;C^{p\&times;p},$ 其中I是单位矩阵，n是系统维数，n为大于0的整数。

4）采用Scaling and Squaring算法求解计算矩阵指数如图4所示，包括如下步骤：

（1）读入算法参数θ，b₁，b₂，…b₁₃，其中θ是scaling因子选择阈值，b_i是Padé逼近系数；

（2）求解scaling因数s，使得s是满足||A/2^s||₁<θ的最小正整数，其中||·||₁是指矩阵的1范数；

（3）对矩阵A进行scaling运算，即A=A/2^s；

（4）求解矩阵e^A的[13/13]阶Pade有理近似r_A，其过程为：

矩阵运算A₂=A²，A₄=A₂A₂，A₆=A₂A₄；

矩阵运算U=A[A₆(b₁₃A₆+b₁₁A₄+b₉A₂)+b₇A₆+b₅A₄+b₃A₂+b₁I]；

矩阵运算V=A₆(b₁₂A₆+b₁₀A₄+b₈A₂)+b₆A₆+b₄A₄+b₂A₂+b₀I；

由矩阵线性方程(-U+V)r_A=U+V求解r_A；

（5）通过对r_A连续平方的方式，计算得到将这一结果作为scaling操作前原矩阵A的矩阵指数e^A的近似。

5）状态变量向量与p维单位向量 $e_p=\left[\begin{array}{c}0\\ .\\ .\\ .\\ 1\end{array}\right]$ 组成增广状态向量，与步骤4）得到的矩阵指数进行矩阵向量乘法，得到 $e^{\mathrm{\&Delta;t}\hat{A}}\left[\begin{array}{c}{x}_n\\ e_p\end{array}\right],$ 取结果向量的前n行赋值给x_n；

6）由y_n=Lx_n得到t_n时刻的输出向量并写入输出数据文件，更新当前时刻为下一时刻t_n=t_n+Δt，仿真向前推进一个步长；

7）比较当前时刻t_n与仿真时间T_end，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则继续到步骤8）；

8）判断当前时刻是否有短路和开路系统故障，或开关动作导致系统模型发生改变的事件，若有，重新建立系统仿真模型，更新状态矩阵A、输出矩阵L；若没有，则跳过更新；

9）返回到步骤3）继续进行下一时刻计算，依此迭代进行，直到仿真结束。

下面结合实例，对本发明的一种基于矩阵指数的电力系统电磁暂态仿真方法做出进一步的说明。

以IEEE123节点算例作为实施例，其系统结构如图1所示。该网络电压等级为4.16kV，在节点150处通过变电站与115kV大电网相连，内部设有单相、三相负荷及线路，是一个复杂不对称的辐射状供电网络。通过特征值分析可以发现，系统状态矩阵中含有大量远离实轴的特征值，如图2所示。这一特征在仿真曲线中即表现为高频振荡的动态过程，其频率最高可达1MHz以上，属于典型刚性的高度振荡（highly oscillatory）系统。对于这一刚性系统，传统数值积分方法进行求解会遇到显著的数值精度和稳定性问题。

第一步，在状态分析框架下建立待研究电力系统的电磁暂态仿真模型，其形式为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+g\left(t\right)$

y(t)=Lx(t)

对于本实例，采用ASMG进行系统建模，状态向量x(t)包含支路电流和部分节点电压，激励源g(t)是115kV母线的等效电压源；

第二步，启动仿真程序，设定仿真时间T_end，仿真步长Δt，多项式近似系数p。设定当前时刻为仿真起始时刻t_n=t₀，依照仿真需要设置仿真初值x_n=x₀。本实例从零状态启动，设定仿真时间0.06s，仿真步长5us，多项式近似阶数为1，即分段线性近似；

第三步，计算当前时刻的激励源g(t_n)，及其1至p-1阶导数。如果g(t)的解析形式未知，则用数值微分近似t_n时刻g(t)的各阶导数。依照W=[g^(p-1)(t₀),…,g’’(t₀),g’(t₀),g(t₀)]的格式将这些向量组成矩阵W。按照 $\hat{A}=\left[\begin{array}{cc}A& W\\ 0& J\end{array}\right]\&Element;C^{(n+p)\&times;(n+p)}$ 的格式更新矩阵，其中矩阵J的格式是固定的 $J=\left[\begin{array}{cc}0& I_{p-1}\\ 0& 0\end{array}\right]\&Element;C^{p\&times;p}.$ 本实例中激励源形式已知，即 $g\left(t\right)=C^{-1}B\left[\begin{array}{c}{V}_g\cos \left(\mathrm{\&omega;t}\right)\\ V_g\cos (\mathrm{\&omega;t}-2/3\mathrm{\&pi;})\\ V_g\cos (\mathrm{\&omega;t}+2/3\mathrm{\&pi;})\end{array}\right],$ $g^{\&prime;}\left(t\right)=C^{-1}B\left[\begin{array}{c}-\mathrm{\&omega;}{V}_g\sin \left(\mathrm{\&omega;t}\right)\\ -\mathrm{\&omega;}{V}_g\sin (\mathrm{\&omega;t}-2/3\mathrm{\&pi;})\\ -\mathrm{\&omega;}{V}_g\sin (\mathrm{\&omega;t}+2/3\mathrm{\&pi;})\end{array}\right],$ $\hat{A}=\left[\begin{array}{ccc}A& g\left(t\right)& g^{\&prime;}\left(t\right)\\ 0& 0& 1\\ 0& 0& 0\end{array}\right];$

第四步，采用Scaling and Squaring算法求解计算矩阵指数具体实现方式如下：

1)读入算法参数θ，b₁，b₂，…b₁₃，其中θ是scaling因子选择阈值，b_i是Padé逼近系数；

2)求解scaling因数s，使得s是满足||A/2^s||₁<θ的最小正整数；

3)对矩阵A进行scaling运算，即A=A/2^s;

4)求解矩阵e^A的[13/13]阶Pade有理近似r_A，其过程为：

矩阵运算A₂=A²，A₄=A₂A₂，A₆=A₂A₄；

矩阵运算U=A[A₆(b₁₃A₆+b₁₁A₄+b₉A₂)+b₇A₆+b₅A₄+b₃A₂+b₁I]；

矩阵运算V=A₆(b₁₂A₆+b₁₀A₄+b₈A₂)+b₆A₆+b₄A₄+b₂A₂+b₀I；

由矩阵线性方程(-U+V)r_A=U+V求解r_A；

5)通过连续平方的方式，计算得到将这一结果作为scaling操作前原矩阵A的矩阵指数近似；

第五步，状态变量向量与p维单位向量 $e_p=\left[\begin{array}{c}0\\ .\\ .\\ .\\ 1\end{array}\right]$ 组成增广状态向量，与第四步得到的矩阵指数进行矩阵向量乘法，得到 $e^{\mathrm{\&Delta;t}\hat{A}}\left[\begin{array}{c}{x}_n\\ e_p\end{array}\right],$ 取结果向量的前n行赋值给x_n；

第六步，由y_n=Lx_n得到t_n时刻的输出向量并写入输出数据文件。更新当前时刻为下一时刻t_n=t_n+Δt，仿真向前推进一个步长；

第七步，比较当前时刻t_n与仿真时间T_end，判断是否已经抵达仿真结束时刻。若已经达到，则仿真结束；若未达到，则继续到第八步；

第八步，判断当前时刻是否有短路和开路等系统故障，或开关动作等导致系统模型发生改变的事件，若有，则重新建立系统仿真模型，更新状态矩阵A、输出矩阵L；若没有，则跳过更新；

第九步，回到第三步继续进行下一时刻计算。依此迭代进行，直到仿真结束。

执行仿真计算的计算机硬件环境为Intel Core2Q84002.66GHz CPU，内存容量2GB；软件环境为Windows7操作系统。

分别将本发明所提出的矩阵指数方法（Matrix Exponential）与典型隐式方法——梯形法（Trap）以及典型显式方法——4阶龙格库塔法（RK4）进行比较，并采用有误差控制的变步长Adams方法作为仿真的准确值。对于矩阵指数方法和隐式梯形法，仿真步长设置为5us；龙格库塔法由于数值稳定性限制，需要采用更小的步长才能抑制误差的数值振荡，仿真步长采用0.1us，是其他仿真方法的1/50。

附图5至图8比较了本发明提出的矩阵指数方法和其它几种数值的仿真结果，其中图5为总体图，图6和图7为放大图，图8为绝对误差曲线。

从放大的图6中可以看出，基于矩阵指数的电磁仿真方法所得仿真结果与小步长的RK4方法、以及作为准确值的变步长Adams法较为接近，而梯形法的仿真结果则与其他方法有显著差异，曲线中存在不合理的尖峰，不能真实反映系统较高频率上的振荡特性。

从进一步放大的图7中可以看出，本发明提出的基于矩阵指数的电磁仿真方法所得仿真结果在时轴的每一个离散点，即图7中的黑色圆圈，相比较于小步长的龙格库塔法，与准确值更靠近。

图8给出了不同仿真方法所得曲线与准确值之间绝对误差的比较，图中的纵坐标是10为底的对数坐标。本发明提出的基于矩阵指数的电磁仿真方法绝对误差在1e-3量级，而传统的隐式梯形法、小步长的显式4阶龙格库塔法绝对误差分别在1e-1、1e-2.5量级。充分验证了本发明提出的仿真方法在数值稳定性和数值精度上的突出优势。

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

Claims (3)

1.一种基于矩阵指数的电力系统电磁暂态仿真方法，其特征在于，包括如下步骤：

1）在状态分析框架下建立待研究电力系统的电磁暂态仿真模型，形式为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+g\left(t\right)$

y(t)=Lx(t)

其中，A是状态矩阵，g(t)是激励源，L是输出矩阵，x(t)是状态变量，y(t)是输出变量；

2）启动仿真程序，设定仿真时间T_end，仿真步长Δt，多项式近似系数p，设定当前时刻为仿真起始时刻t_n=t₀，依照仿真需要设置仿真初值x_n=x₀，其中的下标_n表示t_n时刻的相关变量；

3）利用当前时刻激励源及当前时刻激励源各阶导数信息，组成增广的状态矩阵

4）采用Scaling and Squaring算法求解计算矩阵指数

5）状态变量向量与p维单位向量 $e_p=\left[\begin{array}{c}0\\ .\\ .\\ .\\ 1\end{array}\right]$ 组成增广状态向量，与步骤4）得到的矩阵指数进行矩阵向量乘法，得到 $e^{\mathrm{\&Delta;t}\hat{A}}\left[\begin{array}{c}{x}_n\\ e_p\end{array}\right],$ 取结果向量的前n行赋值给x_n；

6）由y_n=Lx_n得到t_n时刻的输出向量并写入输出数据文件，更新当前时刻为下一时刻t_n=t_n+Δt，仿真向前推进一个步长；

7）比较当前时刻t_n与仿真时间T_end，判断是否已经抵达仿真结束时刻，若已经达到，则仿真结束；若未达到，则继续到步骤8）；

8）判断当前时刻是否有短路和开路系统故障，或开关动作导致系统模型发生改变的事件，若有，重新建立系统仿真模型，更新状态矩阵A、输出矩阵L；若没有，则跳过更新；

9）返回到步骤3）继续进行下一时刻计算，依此迭代进行，直到仿真结束。

2.根据权利要求1所述的一种基于矩阵指数的电力系统电磁暂态仿真方法，其特征在于，步骤3）所述的利用当前时刻激励源及当前时刻激励源各阶导数信息，组成增广的状态矩阵具体方法如下：

计算当前时刻的激励源g(t_n)，及当前时刻激励源g(t_n)的1至p-1阶导数，如果g(t)的解析形式未知，则用数值微分近似t_n时刻g(t)的各阶导数，依照W=[g^(p-1)(t₀),…,g’’(t₀),g’(t₀),g(t₀)]的格式将这些向量组成矩阵W，按照 $\hat{A}=\left[\begin{array}{cc}A& W\\ 0& J\end{array}\right]\&Element;C^{(n+p)\&times;(n+p)}$ 的格式更新矩阵，其中矩阵J的格式是固定的 $J=\left[\begin{array}{cc}0& I_{p-1}\\ 0& 0\end{array}\right]\&Element;C^{p\&times;p},$ 其中I是单位矩阵，n是系统维数，n为大于0的整数。

3.根据权利要求1所述的一种基于矩阵指数的电力系统电磁暂态仿真方法，其特征在于，步骤4）所述的采用Scaling and Squaring算法求解计算矩阵指数包括如下步骤：

（1）读入算法参数θ，b₁，b₂，…b₁₃，其中θ是scaling因子选择阈值，b_i是Padé逼近系数；

（2）求解scaling因数s，使得s是满足||A/2^s||₁<θ的最小正整数，其中||·||₁是指矩阵的1范数；

（3）对矩阵A进行scaling运算，即A=A/2^s；

（4）求解矩阵e^A的[13/13]阶Pade有理近似r_A，其过程为：

矩阵运算A₂=A²，A₄=A₂A₂，A₆=A₂A₄；

矩阵运算U=A[A₆(b₁₃A₆+b₁₁A₄+b₉A₂)+b₇A₆+b₅A₄+b₃A₂+b₁I]；

矩阵运算V=A₆(b₁₂A₆+b₁₀A₄+b₈A₂)+b₆A₆+b₄A₄+b₂A₂+b₀I；

由矩阵线性方程(-U+V)r_A=U+V求解r_A；

（5）通过对r_A连续平方的方式，计算得到，将这一结果作为scaling操作前原矩阵A的矩阵指数e^A的近似。