一种电力系统全过程动态仿真的数值积分方法

摘要

本发明提出了一种电力系统全过程动态仿真的数值积分方法，能够使固定步长的梯形积分和变步长的Gear积分两种方法在全过程动态仿真中有机结合。该方法在电力系统全过程仿真的机电暂态过程中采用固定步长的梯形积分法，动态元件的微分方程和电力网络的代数方程进行迭代求解；在中长期动态过程中采用变步长的Gear法，动态元件的微分方程和电力网络的代数方程联立求解；固定步长和变步长两种方法在仿真中依据一定的策略自动切换。该方法结合了原有两种数值积分法的优点，在保证数值计算稳定性和计算精度的前提下，克服了现有的全过程动态仿真程序在仿真电力系统机电暂态过程时由于步长过小而导致的计算效率偏低的缺陷。

说明书

技术领域

本发明涉及一种电力系统的模拟与计算，特别涉及一种电力系统全过程动态仿真的数值积分方法。

背景技术

电力系统全过程动态仿真是指将电力系统机电暂态、中期和长期过程有机地统一起来进行数字仿真，是电力系统规划设计、调度运行和分析研究的主要工具。电力系统全过程动态仿真能够描述了电力系统受到扰动之后整个连续的动态过程。电力系统从机电暂态过程到中长期的动态过程是连续的，并不是截然分开的：暂态过程对中长期过程有影响，中长期过程对后续新的暂态过程也有作用，它们之间往往是密不可分的。

全过程动态仿真主要涉及有发电机及其励磁系统、动力系统、动态负荷(如感应电动机和同步电动机)、直流输电系统、电力电子元件(如FACTS)等众多动态元件和输电网络、静态负荷等组成的非线性动态系统。全过程动态仿真模拟电力系统动态元件及其控制系统的过渡过程，属于时域仿真。仿真中混合着快速和慢速动态过程，是典型的刚性非线性动态系统。描述这一非线性动态系统的是一组高阶的微分方程组和非线性代数方程组。时域仿真计算就是采用适当的数值积分方法求解微分代数方程组的初值问题。数值积分方法是整个全过程动态仿真技术的基础和核心。

现有全过程动态仿真程序的数值积分方法采用Gear类变步长方法，例如：瑞典的SIMPOW程序和法国的EUROSTAG程序。这种方法的优点是暂态及中长期动态过程的模型和数值积分方法在计算中不需要切换，在中长期动态过程中可以使用较大的步长进行仿真。但在应用实践中，发现这种方法的主要缺点和问题为：仿真计算速度慢，特别是机电暂态过程中计算步长过小，导致仿真时间很长。这种根源于数值积分法的缺点和问题，降低了全过程动态仿真程序的实用价值，不利于程序的大规模推广和工程应用。

而现有的机电暂态稳定仿真中大多采用固定步长的梯形积分法，例如：国内使用很广泛的PSD-BPA暂态稳定分析程序和PSASP程序。这种方法的优点是编程简单，机电暂态过程中计算步长比Gear类法大，计算速度快，效率高。但是它的缺点和问题为：在中长期动态过程中步长与机电暂态过程基本相同，仿真速度较慢。

发明内容

本发明的目的在于，提供一种电力系统全过程动态仿真的数值积分方法。

本发明在电力系统全过程动态仿真的数值积分方法中，使固定步长的梯形积分和变步长的Gear积分两种方法能够在全过程动态仿真中有机结合：(1)在机电暂态过程中采用固定步长梯形积分法，其中，动态元件和动态系统的微分方程和电力网络的代数方程交替迭代求解。(2)在中长期动态过程中采用变步长的Gear类方法，动态元件的微分方程和电力网络的代数方程联立求解。该方法能够在保证数值稳定性和仿真精度的前提下，大大提高仿真计算速度。

电力网络的代数方程以下简称网络方程。

其中，包括以下步骤：

步骤101：开始仿真，设置仿真时间t＝0；

步骤102：判断当前积分方法的类型；

步骤103：如果是固定步长方法，则设置迭代次数k1为0，开始使用固定步长的梯形积分方法进行计算；

步骤110：如果是变步长的Gear方法，则进行预测，并设置迭代次数k2为0，开始使用变步长的Gear积分方法进行计算；

步骤120：判断积分是否成功；

步骤121：判断有无故障或操作发生；

步骤122：如果有故障或操作发生，则计算故障t⁺时刻的积分初值；

步骤123：进行积分方法选择；

步骤124：判断仿真结束时间是否到；

步骤125：如果仿真结束时间未到，则仿真时间增加一个步长，重复上述步骤102—124；

在步骤103—110之间，还包括如下步骤：

步骤104：使用梯形积分法计算动态元件的微分方程；

步骤105：求网络方程中的注入电流；

步骤106：求解网络方程，得到母线电压；

步骤107：判断是否收敛；

步骤108：判断是否达到最大迭代次数；

步骤109：如果没有达到最大迭代次数限制，则迭代次数k1加1；

在步骤110—120之间，还包括如下步骤：

步骤111：解整个系统的联立方程组，进行校正计算；

步骤112：判断是否收敛；

步骤113：截断误差计算；

步骤114：判断截断误差是否小于容许值；

步骤115：变阶变步长计算，求取积分的步长和阶数；

步骤116：判断是否达到最大迭代次数；

步骤117：迭代次数k2加1；

步骤118：判断步长h是否小于设置的最小步长；

步骤119：变阶变步长计算，求取积分的步长和阶数。

本发明的有益效果是：在国内外全过程动态仿真领域首次提出固定步长的梯形积分法和变步长的Gear法相结合的混合数值积分方法。在仿真中可以根据系统动态过程的特点自动选择合适的积分方法，从而在保证数值稳定性和仿真精度的前提下，大大缩短了仿真时间，提高了程序的计算效率，克服了现有的全过程动态仿真程序在仿真电力系统机电暂态过程时由于步长过小而导致的计算效率偏低的缺陷。

附图说明

图1为本发明的一种电力系统全过程动态仿真的数值积分方法的流程图。

具体实施方式

以下，参考附图1详细描述本发明的一种电力系统全过程动态仿真的数值积分方法。

电力系统复杂的长时间动态过程中混合着快速动态过程和慢速动态过程。电力系统受到扰动后的暂态过程一般在扰动后10秒左右结束。然而在严重故障或连锁反应故障的冲击下，系统发电和负荷之间的有功功率或无功功率可能出现长期持续偏移的不平衡状态而引起潮流、电压和频率等电气量和原动机系统变量的长期变化过程，并最终导致系统失去稳定。在电力系统中长期过程中，一般不直接导致损失负荷，而是激发一系列系统解列形成孤岛、失去电源，并最终由于自动低压低频减载或由于整个孤岛系统的崩溃而损失负荷。

电力系统全过程动态仿真用于研究在受到干扰之后系统较长时间的机电过渡过程。研究的时间从几秒到数十分钟甚至若干个小时，时间跨度大。不仅包含了系统暂态稳定过程，而且与暂态稳定相比，全过程动态仿真变得更复杂。它的主要特点是：①模型的种类众多，微分代数方程的阶数高。除覆盖电力系统暂态稳定、中期和长期过程所需要的动态元件模型外，还要考虑更多的自动装置模型。例如，暂态稳定不予考虑的锅炉及其控制、水力系统模型、自动发电控制、变压器分节头的自动调整等。②仿真时间跨度大，可从几秒到几十分钟。必须采用具有自动变阶变步长功能的数值积分算法，步长范围可从小于1毫秒到几十秒。在系统快变阶段，采用小步长；在慢变阶段，采用大步长，以缩短计算时间。③控制系统中的最大最小时间常数相差大，有的相差近千倍，使系统的刚性比加大。如：锅炉汽包时间常数Cd典型值为125秒，而快速励磁系统的时间常数典型值为0.01秒，电力电子装置时间常数可能更小。而一般的暂态稳定计算中，不考虑中长期过程动态元件的过程，系统的刚性不明显。④变步长算法和微分方程刚性度的加强要求联立求解大型方程组，这是由于通常用于暂态稳定计算的微分方程和代数方程交替求解算法不能使用较大的步长计算，否则，交接误差会越来越大，导致数值计算不稳定。

本发明的电力系统全过程动态仿真的数值积分方法主要由三部分组成：固定步长的梯形积分法、变步长的Gear法和两种方法的自动切换策略。单独的固定步长梯形积分法和单独的变步长的Gear法(小于等于2阶)的数值稳定域都包含复平面的左半部分，满足电力系统全过程动态仿真的要求，且目前分别都已在成熟的电力系统仿真程序中得到应用。本发明的特点在于在仿真中可以使两种积分方法得到有机结合，且能够扬长避短，根据系统动态过程的特点自动选择合适的积分方法，从而在保证数值稳定性和仿真精度的前提下，大大缩短了仿真时间，提高了程序的计算效率。克服了以往使用单独的定步长梯形积分法不能适应中长期动态仿真的需要；而单独的变步长Gear法在电力系统的暂态稳定阶段积分步长过小，仿真速度过慢的缺陷。

积分方法切换策略的依据机电暂态过程和中长期动态过程的不同特点：①机电暂态过程属于系统的快变阶段，电网最低母线电压一般在0.8pu以下，求解变量值变化剧烈。此时，如果采用Gear法，则步长会较小，一般小于0.01sec.，甚至常常低于0.001sec.，求解较慢；如果采用固定步长的梯形积分法，则虽然迭代次数会达到2～3次以上，但是固定步长积分法使用微分代数方程组交替求解，不需要联立求解大型方程组，所以速度相对较快。②中长期动态过程属于系统的慢变阶段，电网最低母线电压一般在0.8pu以上，求解变量值变化缓慢。此时，如果采用Gear法，则步长会较大，一般大于0.05sec.，甚至达数秒以上；如果采用固定步长的梯形积分法，则虽然迭代次数一般为2～3次，但由于步长一般固定在0.01sec.左右，仿真速度不如大步长的Gear法快。由①和②知，根据母线电压、步长和固定步长积分法的迭代次数可以判断系统处于机电暂态或中长期动态过程。此外，为保证可靠的判断系统动态过程，一种算法要持续计算一定步数。

电力系统全过程动态仿真中需要求解的是初始值已知的微分代数方程：

式1为微分方程，表示电力系统动态元件特性，是系统的状态方程，其中，y由y₁和y₂两部分变量组成，y₁为n₁个应变量的状态向量(微分变量)，由全过程动态计算中应考虑的各种动态元件决定；y₂为n₂个应变量的代数向量(代数变量)；t为自变量(时间)。

式2为代数方程，表示电力系统静态元件特性，主要是系统的网络方程，变量的含义同式1。

式3表示系统变量y的初值。

仿真目的是在y₁、y₂和t分别等于y₁₀、y₂₀和t₀的初始条件下使用数值积分方法求解作为t的函数的y₁和y₂。

定义N为n₁和n₂之和，即微分代数方程组的阶数。电力系统全过程动态仿真对象通常是大规模交直流电力系统，其规模可从几千到上万个母线和支路、数千台发电机及其控制系统、若干条到数十条直流输电线路。因此，式1和式2中的方程阶数可以到数万阶以上。

对于变步长的Gear积分法，其计算步骤主要包括预测、校正、截断误差计算、自动变阶变步长控制四个步。Gear法起步时使用1阶，由于只有2阶及以下的Gear法的稳定域能够覆盖复平面的左半平面，所以使用的最大阶数为2阶。这与常见的Gear法大部分是相同的，其中有所不同的是，在校正计算中，对微分和代数变量采用了EUROSTAG程序的处理方法：即微分变量采用ADAMS法，代数变量采用BDF法。ADAMS和BDF法的使用体现在方法的系数向量L_ADAMS和L_BDF上。

此外，Gear法中需要用到Nordsieck向量z，z由变量y的一阶和二阶导数y′(或)、y＂(或)组成，用以代替前两步的y或f(y，t)数据。

$z_n={[y_n,{\mathrm{hy}}_n^{\prime },\frac{h^2}{2!}{y}_n^{\prime \prime }]}^T$  (式4)

式4中：n为上一步的计算步数；h为步长；z_n表示上一步的Nordsieck向量。

对于固定步长的梯形积分方法，为说明方便，电力系统的微分代数方程组(包括动态元件及其控制系统、静态负荷、直流输电系统、FACTS元件、网络方程等)表示为如下形式：

式5和式6中：X为系统状态向量(微分变量)；为母线电压向量；I为母线电流向量；Y_N为根据电力网络及其参数形成的导纳矩阵，是网络方程的重要组成部分，一般为高阶稀疏矩阵。需要说明的是：式5～式6与式1～式2表达的方程式内容和含义完全相同的，只是表达形式不同。

根据梯形积分公式，对于时刻t＝t_n，X＝X_n，V＝V_n，求X在时刻t_n+1＝t_n+Δt的解X_n+1，可用式7的差分方程求解：

$X_{n+1}=X_n+\frac{\mathrm{\Delta t}}{2}[f(X_n,V_n)+f(X_{n+1},V_n)]$  (式7)

V在时刻t_n+1的解V_n+1，可由式6得：

I(X_n+1，V_n+1)＝Y_NV_n+1  (式8)

由式7和式8可见，X_n+1和V_n+1是一组隐式非线性方程的解，必须迭代求解。

附图1是根据本发明实施例的适合电力系统全过程动态仿真的数值积分方法的流程图。

如附图1所示，根据发明实施例的适合电力系统全过程动态仿真的数值积分方法包括如下步骤：

步骤101：开始积分，设置积分法为固定步长的梯形积分法；设置仿真时间t为0；积分步数n为0。

以下说明中的t_n表示是第n个时步，本步骤为t₀。

步骤102：判断积分方法为固定步长的梯形积分法和变步长的Gear法中的哪一个。

如果是固定步长的梯形积分，则执行步骤103；否则执行步骤110。

步骤103：设置固定步长的梯形积分法的迭代次数k1为0，并取前一时间步的值作为初值。

本步骤开始到步骤109结束为固定步长的梯形积分法在一个时步上的求解，该方法的说明使用的电力系统方程组形式为式5～式6。

本步骤中取前一时间步t_n的值作为初值，即

X_n+1⁽⁰⁾＝X_n，V_n+1⁽⁰⁾＝V_n

步骤104：使用梯形积分法计算动态元件的微分方程。

该步骤的目的是求解微分方程式5中的X。根据梯形积分公式，对于时刻t＝t_n，X＝X_n，V＝V_n，求X在时刻t_n+1＝t_n+Δt的第k1+1步迭代的解X_n+1^(k1+1)，可用式9的差分方程求解：

${X_{n+1}}^{(k1+1)}=X_n+\frac{\mathrm{\Delta t}}{2}[f(X_n,V_n)+f({X_{n+1}}^{\left(k1\right)},{V_n}^{\left(k1\right)})]$  (式9)

求得系统状态向量X_n+1^(k1+1)。

步骤105：求解网络方程的注入电流I^(k1+1)(X_n+1^(k1+1)，V_n+1^(k1))。

存在注入电流的元件主要有：发电机、负荷、无功补偿装置、直流输电等。本步骤利用X_n+1^(k1+1)和V_n+1^(k1)求出其向网络方程的注入电流I^(k1+1)(X_n+1^(k1+1)，V_n+1^(k1))。

步骤106：求解网络方程，见式10，得到新的母线电压向量V_n+1^(k1+1)。

I^(k1+1)(X_n+1^(k1+1)，V_n+1^(k1))＝Y_NV_n+1^(k1+1)  (式10)

求出系统电压向量V_n+1^(k1+1)。

式6中的导纳矩阵Y_N除了因故障或操作而改变外，在相邻两次故障或操作之间保持不变。因此为了节省时间和内存，每当因故障或操作而对导纳矩阵进行修改以后，应随即对它进行稀疏三角分解，以反复应用于解式10的消去、回代过程，直到下一次故障或操作为止。

步骤107：检查两次迭代的最大电流误差值，判断迭代计算是否收敛。

最大电流误差值的计算方法是：①计算相邻两次迭代之间电流差向量I_ERR，即I_ERR＝I^(k1+1)-I^(k1)，求出向量I_ERR中绝对值最大的元素值，该值就是最大电流误差值。

如果最大电流误差在收敛允许值内，则迭代计算收敛，完成本次交替计算，置积分成功标志为1，跳到步骤120。

如果最大电流误差大于收敛允许值，则以X_n+1^(k1+1)和V_n+1^(k1+1)为新的初值，执行步骤108，继续迭代。

步骤108：判断是否达到最大迭代次数，并进行相关处理

如果达到最大迭代次数，则表示无法求出此时刻的微分代数方程组的解，置积分成功标志为0，执行步骤120。

如果没有达到最大迭代次数，则执行步骤109。

步骤109：迭代次数k1加1，并跳到步骤104。

步骤110：设置变步长的Gear法的迭代次数k2为0，并对计算变量进行预测。

本步骤开始到步骤119进行变步长的Gear法在一个时步上的求解，该方法的说明使用的电力系统方程组形式为式1～式2。。

本步骤中预测下一步(即n+1步)的Nordsieck向量z_n+1，见式11：

$z_{n+1}^{\left(0\right)}={\mathrm{Pz}}_n$  (式11)

式11中P是Pascal上三角矩阵，它的第(i，j)元素是：当j≥i时为$\left(\begin{array}{c}j\\ i\end{array}\right)=\frac{j!}{(j-i)!i!},j\ge i\ge 0,$当j<i时为零。

预测后的Nordsieck向量z_n+1，即微分变量和代数变量也可写为：$z_{1,n+1}={[y_{1,n+1}^{\left(0\right)},h{\stackrel{·}{y}}_{1,n+1}^{\left(0\right)},\frac{h^2}{2!}{\stackrel{··}{y}}_{1,n+1}^{\left(0\right)}]}^T$和$z_{2,n+1}={[y_{2,n+1}^{\left(0\right)},h{\stackrel{·}{y}}_{2,n+1}^{\left(0\right)},\frac{h^2}{2!}{\stackrel{··}{y}}_{2,n+1}^{\left(0\right)}]}^T$

步骤111：校正计算

在变步长Gear法的迭代过程中，电力系统全过程动态仿真中的微分代数方程组(式1和式2)变换为如下的校正计算方程组(参见式12和式13)。方程组的未知数(求解量)为：Δy_1，n+1和Δy_2，n+1，这个方程组通常采用牛顿法进行求解：

如果牛顿法迭代求解收敛，则求得Δy_1，m+1和Δy_2，m+1，置收敛标志为1，并使用如下校正公式进行计算：

$\left(\begin{array}{c}{z}_{1,n+1}=z_{1,n+1}^{\left(0\right)}+L_{\mathrm{ADMAS}}{\mathrm{\Delta y}}_{1,n+1}\\ z_{2,n+1}=z_{2,n+1}^{\left(0\right)}+L_{\mathrm{BDF}}{\mathrm{\Delta y}}_{2,n+1}\end{array}\right)$  (式14)

对于一阶积分方法，式14可写为如下的式15和式16：

$\left(\begin{array}{c}{y}_{1,n+1}=y_{1,n+1}^{\left(0\right)}+L_{0,\mathrm{ADMAS}}{\mathrm{\Delta y}}_{1,n+1}\\ {h\stackrel{·}{y}}_{1,n+1}={h\stackrel{·}{y}}_{1,n+1}^{\left(0\right)}+L_{1,\mathrm{ADMAS}}{\mathrm{\Delta y}}_{1,n+1}\end{array}\right)$  (式15)

$\left(\begin{array}{c}{y}_{2,n+1}=y_{2,n+1}^{\left(0\right)}+L_{0,\mathrm{BDF}}{\mathrm{\Delta y}}_{2,n+1}\\ {h\stackrel{·}{y}}_{2,n+1}={h\stackrel{·}{y}}_{2,n+1}^{\left(0\right)}+L_{1,\mathrm{BDF}}{\mathrm{\Delta y}}_{2,n+1}\end{array}\right)$  (式16)

对于二阶积分方法，式14可写为如下的式17和式18：

$\left(\begin{array}{c}{y}_{1,n+1}=y_{1,n+1}^{\left(0\right)}+L_{0,\mathrm{ADMAS}}{\mathrm{\Delta y}}_{1,n+1}\\ {h\stackrel{·}{y}}_{1,n+1}={h\stackrel{·}{y}}_{1,n+1}^{\left(0\right)}+L_{1,\mathrm{ADMAS}}{\mathrm{\Delta y}}_{1,n+1}\\ \frac{h^2}{2!}{\stackrel{··}{y}}_{1,n+1}=\frac{h^2}{2!}{\stackrel{··}{y}}_{1,n+1}^{\left(0\right)}+L_{2,\mathrm{ADMAS}}{\mathrm{\Delta y}}_{1,n+1}\end{array}\right)$  (式17)

$\left(\begin{array}{c}{y}_{2,n+1}=y_{2,n+1}^{\left(0\right)}+L_{0,\mathrm{BDF}}{\mathrm{\Delta y}}_{2,n+1}\\ {h\stackrel{·}{y}}_{2,n+1}={h\stackrel{·}{y}}_{2,n+1}^{\left(0\right)}+L_{1,\mathrm{BDF}}{\mathrm{\Delta y}}_{2,n+1}\\ \frac{h^2}{2!}{\stackrel{··}{y}}_{2,n+1}=\frac{h^2}{2!}{\stackrel{··}{y}}_{2,n+1}^{\left(0\right)}+L_{2,\mathrm{BDF}}{\mathrm{\Delta y}}_{2,n+1}\end{array}\right)$  (式18)

如果牛顿法迭代求解不收敛，则置不收敛标志为0。

对于一阶积分方法，向量L_ADAMS和L_BDF的系数相同为：

$\left(\begin{array}{c}{L}_{0,\mathrm{ADMAS}}=1\\ L_{1,\mathrm{ADMAS}}=1\end{array}\right)$和$\left(\begin{array}{c}{L}_{0,\mathrm{BDF}}=1\\ L_{1,\mathrm{BDF}}=1\end{array}\right)$

对于二阶积分方法，向量L_ADAMS和L_BDF的系数分别为：$\left(\begin{array}{c}{L}_{0,\mathrm{ADMAS}}=\frac{1}{2}\\ L_{1,\mathrm{ADMAS}}=1\\ L_{2,\mathrm{ADMAS}}=\frac{1}{2}\end{array}\right)$和$\left(\begin{array}{c}{L}_{0,\mathrm{BDF}}=\frac{2}{3}\\ L_{1,\mathrm{BDF}}=1\\ L_{2,\mathrm{BDF}}=\frac{1}{3}\end{array}\right)$

步骤112：根据步骤111中的收敛标志，判断校正计算是否收敛

如果收敛，执行步骤113；否则，执行步骤116。

步骤113：截断误差计算

在积分过程中，为使计算保证在规定的精度内，同时又尽量减少计算量，需要不断自动改变阶次和步长。为了达到这一目的，需要进行截断误差计算。

K阶方法第n+1步的截断误差，见式19：

  (式19)

式19中：C_K+1与阶有关的常数；是第n+1步计算所得的Z_n+1与前一步的Z_n这两个向量的最后一个分量相减的值。

仿真中常使用相对误差ε，即要求每步的相对误差小于事先规定的一个值ε₀。ε的计算见式20。

$ϵ=\sqrt{\underset{j=1}{\overset{N}{\Sigma }}{\left[\frac{C_{K+1}K!▿Z_{\left(j\right)n+1}^K}{y_{\left(j\right)}}\right]}^2}$  (式20)

式11中，y_(j)为积分到这一步时的值，若y_(j)＝0，则应取y_(j)＝1；K为阶次；N为变量y的个数；对于K阶积分方法，C_K+1的值为

步骤114：根据相对误差ε，判断截断误差是否小于精度容许值ε₀

ε在步骤113中计算得到，ε₀是预先设定的精度要求(容许值)。

如果ε≤ε₀，则认为这一步积分成功，置积分成功标志为1，转入步骤115，进行变阶变步长计算。

如果ε>ε₀，则认为这一步积分不成功，置积分成功标志为0，转入步骤118进行相应处理。

步骤115：变阶变步长计算

当前步长为h_n+1。对于下一步计算，设所用的步长变化系数为R，则下一步计算步长为Rh_n+1。由于截断误差是估计值，而非准确计算值，因此，R值应乘以一个安全系数，在原阶不变的情况下此系数取在变阶变步长的情况下，考虑到误差估计更间接，且变阶要增加额外的工作量，多费时间，此系数取得更加严格，降阶情况下取升阶情况下取

本步骤中，先计算阶次不变时的R_K、降阶时的R_K-1和升阶时的R_K+1值：

$R_K=\frac{1}{1.2}\sqrt[K+1]{\frac{{ϵ}_0}{\underset{j=1}{\overset{N}{\Sigma }}\left[\frac{C_{K+1}K!▿Z_{\left(j\right)n+1}^K}{y_{\left(j\right)}}\right]}}$  (式21)

$R_{K-1}=\frac{1}{1.3}\sqrt[K+1]{\frac{{ϵ}_0}{\underset{j=1}{\overset{N}{\Sigma }}\left[\frac{C_{K+1}K!▿Z_{\left(j\right)n+1}^K}{y_{\left(j\right)}}\right]}}$  (式22)

式22中，是Z_(j)n+1向量中的最后一个分量。

$R_{K+1}=\frac{1}{1.4}\sqrt[K+2]{\frac{{ϵ}_0}{\underset{j=1}{\overset{N}{\Sigma }}\left[\frac{C_{K+1}K!{▿}^2{Z}_{\left(j\right)n+1}^K}{y_{\left(j\right)}}\right]}}$  (式23)

在积分过程中，R_K、R_K-1和R_K+1都进行计算，选取其中最大的R，以确定下一步采用合理的阶和步长。如果R_K最大，则不变阶只变步长；如果R_K-1最大，则降阶变步长；如果R_K+1最大，则升阶变步长。

变阶变步长后，Nordsieck向量Z_n的变换为新值新值的计算方法见如下的式24～式27。

不变阶只变步长时，有：

  (式24)

降阶变步长时，有

${\hat{Z}}_n={[y_n,{\alpha }_{k-1}{\mathrm{hy}}_n^{\prime },···,\frac{{\alpha }_{k-1}^{k-1}{h}^{k-1}}{(k-1)!}{y}_n^{(k-1)}]}^T$  (式25)

升阶变步长时，有

${\hat{Z}}_n={[y_n,{\alpha }_{k+1}{\mathrm{hy}}_n^{\prime },···,\frac{{\alpha }_{k+1}^{k+1}{h}^{k+1}}{(k+1)!}{y}_n^{(k+1)}]}^T$  (式26)

式26中，是未知的，需用向后差分求得

$y_n^{(k+1)}\approx \frac{y_n^{\left(k\right)}-y_{n-1}^{\left(k\right)}}{h}$  (式27)

式24～式27中，k是阶次。

仿真中，每步都改变步长和阶不一定好，故做了下面的考虑：

1)如果这一步积分失败，则估计R；

2)在上一次改变阶或步长以后的k+1步估计R(使阶或步长的改变不频繁)；

3)如果在上次估计R时步长不放大，则在其10步后估计R(以减少过于频繁的附加计算)。

步骤116：判断迭代次数k2是否达到限定次数。

限定次数通常设为5。如果k2超过限定次数，则说明此式刻t_n+1无解，置积分成功标志为0，跳到步骤120。如果k2没有超过限定次数，则执行步骤117。

步骤117：迭代次数k2加1。

步骤118：判断是否达到最小步长。

最小步长可设为0.00001sec.，如果当前步长已经达到最小步长值，则说明此式刻t_n+1无解，置积分成功标志为0，跳到步骤120。否则，执行步骤119。

步骤119：本步阶段误差大于容许值后的变阶变步长。

该步骤的变阶变步长的方法同步骤115中完全一样。本步骤结束后，跳到步骤110，重新进行该时步的计算。

步骤120：根据本步积分成功标志，判断是否继续积分计算。

如果积分成功标志为1，继续计算，则执行步骤121；否则，如果积分成功标志为0，则结束仿真。

步骤121：判断此时有无故障或操作发生。

如果有，则执行步骤122；否则，执行步骤123。

步骤122：求取故障或操作后瞬间t⁺时刻的积分初值。

故障或操作发生后，首先修改导纳阵，然后求取故障或操作后瞬间t⁺时刻的母线电压及其它变量的积分初值。下面针对通常考虑的对称故障和操作为例说明导纳阵的修改方法：

对于对称短路或操作，包括三相短路和元件的三相开断、串联电容的强行补偿等。除了三相短路外，其它的操作显然都可以处理为网络中某些接地支路或不接地支路的参数发生改变，从而很容易修改导纳阵元素。对于三相短路，如果计及电弧影响，可以在短路点接入相应的阻抗，否则可以接入一个足够大的导纳，使短路点的对地阻抗接近于零。然后，根据短路点的位置，修改导纳阵中的对应元素。

求取故障或操作后瞬间t⁺时刻的母线电压，是一个迭代过程，迭代中包含两步，①计算注入电流。②根据式6求得故障后母线电压。第②步得到新的故障后母线电压，返回第①步继续计算注入电流，再求取新的故障后母线电压，直到相邻两次迭代中的母线电压误差小于允许值为止。

求得故障或操作后瞬间t⁺时刻的母线电压后，由于状态(微分)变量的值在故障或操作前后保持不变，则其它变量的积分初值很容易求得。

步骤123：下一步仿真时积分方法的选择。

积分方法的选择，即根据切换策略选取固定步长的梯形积分和变步长的Gear积分法中的一种。由于固定步长的梯形积分法和变步长的数值积分法都是自启动的，且两者单独使用时计算结果完全相同，所以，计算中两种方法可以根据系统动态过程进行自动平滑切换。在机电暂态过程中使用固定步长积分法；在中长期动态过程中使用变步长数值积分。

(1)切换到固定步长积分法的判断条件

①系统发生故障或操作之后导致网络结构改变或进行了切机切负荷等动态元件的投入或切除，例如：输电线路发生短路故障而跳开、发电机切机等。

②在采用变步长的数值积分法时，步长小于0.01秒，且持续30步。

(2)切换到变步长积分法的判断条件

①采用固定步长的积分法的迭代次数只有两次，且持续40步。

②电网最低母线电压大于0.85pu，迭代次数只有两次，且持续20步。

步骤124：判断仿真结束时间是否已到。

如果仿真结束时间是否已到，则结束仿真；否则，执行步骤125。

步骤125：将仿真时间增加一个时间步长h，即n＝n+1，t＝t+h，继续计算下一个时步。

重复步骤102至步骤123，直到仿真时间结束。

以上是为了使本领域普通技术人员理解本发明，而对本发明进行的详细描述，但可以想到，在不脱离本发明的权利要求所涵盖的范围内还可以做出其它的变化和修改，这些变化和修改均在本发明的保护范围内。