CN103700036A - 一种适于电力***多时间尺度的暂态稳定性投影积分方法 - Google Patents

Abstract

本发明的一种适于电力***多时间尺度的暂态稳定性投影积分方法，采用显式交替求解算法对电力***数学模型进行求解，关于电力***数学模型中微分方程的求解则采用投影积分算法：首先进行若干小步长的积分运算，对应***的快动态过程；而后根据小步长积分计算结果，进行一个投影步长的计算，与***的慢动态过程对应。其中，小步长积分计算过程称为内部积分器，采用数值稳定性较好的显式四阶龙格-库塔法以提高算法的稳定性；大步长积分计算过程称为投影积分器。本发明不仅适用于传统电力***暂态稳定性仿真，也适用于智能配电网动态仿真问题，在满足仿真数值稳定性和数值精度要求的同时，提升了仿真计算速度，为高效、可靠的电力***仿真程序的开发奠定了良好的基础。

Description

一种适于电力***多时间尺度的暂态稳定性投影积分方法

技术领域

本发明涉及一种电力***暂态稳定性仿真方法。特别是涉及一种适于电力***多时间尺度的暂态稳定性投影积分方法。

背景技术

电力***是由发电、变电、输电、配电、用电等设备和相应的辅助***组成的复杂大***。在电力***层面，相关问题的分析与研究往往无法直接在实际***上进行试验，因此必须采用有效的数字仿真作为重要的研究手段。

在电力***时域仿真的研究中，针对***动态过程的不同时间尺度分别发展出电磁暂态仿真、机电暂态仿真和中长期动态仿真三种电力***数字仿真方法，三者从元件数学模型到仿真计算方法都具有不同的特征。电力***电磁暂态仿真侧重于***中电场与磁场相互影响产生的电压电流的变化过程；机电暂态仿***要研究电力***在大扰动下（如故障、切机、切负荷、重合闸操作等情况）的动态行为和保持同步稳定运行的能力，即暂态稳定性，所关注的时间范围通常为几秒至几十秒，因而也称为暂态稳定性仿真；中长期动态过程仿真是电力***受到扰动后较长过程的动态仿真，即通常的电力***长过程动态稳定计算。

电力***暂态稳定性仿真除关注传统电力***的暂态稳定运行能力外，还侧重于分析含各种分布式电源及储能装置的智能配电***运行时其工频电气量在***扰动下（开关操作、故障、分布式电源及负荷波动等）的动态响应特性，此时也可称为智能配电***动态仿真。智能配电***接入了种类繁多的分布式电源和大量的电力电子装置，包括旋转电机和各种静态直流型分布式电源，具有明显的多时间尺度特征。电力***暂态稳定性仿真的多时间尺度问题在数学上可以归结为刚性问题，其对所采用的数值算法的精度和稳定性要求更高。

电力***暂态稳定性仿真本质上可归结为对动力学***时域响应的求取，分为数学建模和模型求解两部分。暂态稳定性仿真首先根据元件间拓扑关系将电力***各元件模型构成全***模型，形成一组联立的微分-代数方程组，然后以稳态工况或潮流解为初值，求解扰动下的数值解，即逐步求得***状态量和代数量随时间的变化曲线。

电力***数字仿真建模的过程则是根据***仿真关注的时间尺度范围，由物理原型抽象出数学模型。电力***暂态稳定性仿真中的数学模型包括两部分：描述设备动态特征的微分方程和描述设备之间电气联系的代数方程。

动态设备之间的电气连接关系在运行过程中可能改变，如负荷的投切、机组的启停、线路开断和重合闸等操作，若计及继电保护装置，还应包含大量连续和（或）离散的逻辑时变参数。一般可将电力***数学模型通过一个高维非线性且连续自治的微分-代数方程组来描述，如式(1)所示。

$\left\{\begin{array}{c}\stackrel{\·}{x}=f(x,y)\\ 0=g(x,y)\end{array}\right.---\left(1\right)$

式中，

为微分方程，为代数方程，

为***状态变量，代表电机转子转速、电力电子器件控制***和负荷动态参数等，

为代数变量，表征母线电压幅值与相角。数学模型的求解一般通过特定的数值算法和相应的仿真程序来实现。因此，电力***暂态稳定性仿真在数学上可以归结为求解一个微分-代数方程组的初值问题。

暂态稳定性仿真算法按照对于式(1)中微分方程和代数方程解算形式的不同可以分为交替求解法和联立求解法两大类。交替求解法首先采用特定的数值积分算法，根据初始化计算结果求解微分方程，得到本时步状态变量的值，然后将其代入到代数方程中求解，得到该时步代数变量的值，最后再将代数变量代入微分方程进行下一时步状态变量求解，以此类推实现微分-代数方程组的交替求解；联立求解法则是将微分方程差分化之后，和代数方程联立成一个完整的代数方程组，同时求解状态变量和代数变量。

对于式(1)中的微分方程，除少数可得到解析解以外，大多数只能采用数值解法进行求解。微分方程数值解求解有两种基本方法：一种是把近似解表示成有限个独立函数之和；另一种是差分法，该方法常用于电力***数字仿真。差分法是寻求在一系列离散点上x的近似值的一类方法的统称，两相邻离散点之间的距离为步长，记为h。根据状态变量x(t)在t_n+1时刻的数值解x(t_n+1)能否由t_n时刻的值直接求出，差分法又可分为单步法(one step method)和线性多步法(linear multistep method)，其中，单步法在电力***暂态稳定性仿真中应用广泛。

根据求解过程的不同，单步法可分为显式积分方法和隐式积分方法，显式积分方法可根据当前时刻状态变量直接计算下一时刻状态变量，而隐式积分方法则需要对含有当前时刻和下一时刻状态变量的方程进行求解才能求得下一时刻状态变量。常见的显式积分方法包括欧拉法、改进欧拉法和龙格-库塔法，而隐式积分方法主要有后向欧拉法和隐式梯形法。电力***动态过程时间尺度差异较大，具有较明显的刚性特征。显式积分方法每一时步内的运算量较小，但由于其数值稳定性较差，即在某一时步计算中产生的误差在以后逐步积分过程中不断累积而可能导致计算结果数值不收敛，因此针对刚性问题求解往往只能采取较小的仿真步长，仿真速度受到极大限制。隐式积分方法虽然数值稳定性较好，可以在刚性问题的求解过程中保证数值稳定性，但其每一时步都需要迭代求解方程组，相比显式积分算法其计算与编程工作复杂，极大地限制了其在大规模电力***暂态稳定性仿真中的应用。

可见，开发一种计算简单、数值稳定性好、适用于刚性问题的电力***暂态稳定性仿真方法十分重要。

发明内容

本发明所要解决的技术问题是，提供一种不仅适用于传统电力***暂态稳定性仿真，也适用于智能配电网动态仿真问题的适于电力***多时间尺度的暂态稳定性投影积分方法。

本发明所采用的技术方案是：一种适于电力***多时间尺度的暂态稳定性投影积分方法，包括如下步骤：

1）输入电力***潮流数据、动态元件数据以及故障及操作设置信息，设置仿真参数；

2）根据***算例参数和动态数据形成全***节点导纳矩阵，其中包括***节点导纳矩阵和并入动态元件的节点导纳矩阵，并重新进行潮流校核；

3）对***算例的动态元件进行仿真初始化计算；

4）设置仿真计算时间t=0；

5）设置当前投影积分算法内部积分器的积分步数s=1；

6）采用投影积分算法内部积分器对电力***模型积分一个步长，步长为h，然后设置仿真计算时间t=t+h，内部积分器的积分步数s=s+1，内部积分器采用显式交替求解方法对微分-代数方程组进行求解，对其中的微分方程求解则采用显式四阶龙格-库塔法；

7）根据步骤1）输入的故障及操作设置信息判断***是否发生故障或操作，若发生，则返回步骤5），否则进入下一步骤；

8）判断仿真计算时间t是否大于设置时间T，若计算时间t大于设置时间T，则仿真结束，否则进入下一步骤；

9）判断投影积分算法内部积分步数s是否大于设置步数k+1，若不大于，则返回步骤6），否则进入下一步骤；

10）根据步骤1）输入的故障及操作设置信息判断未来Mh时间内是否存在故障或操作，若存在，则返回步骤5），否则进入下一步骤；

11）判断未来Mh时间内仿真计算时间t是否会超过设置时间T，若超过，则返回步骤5），否则进入下一步骤；

12）根据内部积分器得到的x(t_n+k)和x(t_n+k+1)的值进行一步外部积分步，步长为内部积分器积分步长的M倍，利用式x(t_n+k+1+M)=(M+1)x(t_n+k+1)-Mx(t_n+k)得到x(t_n+k+1+M)，然后求解代数方程g(x(t_n+k+1+M),y(t_n+k+1+M))=0得到y(t_n+k+1+M)，设置t=t+Mh；

13）判断此时仿真计算时间t是否达到设置的仿真时间T，若达到，则仿真结束，否则返回步骤5），依此反复进行直至仿真结束。

步骤1）中所述的仿真参数，包括仿真时间T、仿真步长h和投影积分算法参数，其中，所述的投影积分算法参数包括内部积分器的积分步数k和外部积分器步长相对内部积分器步长的倍数M。

所述的内部积分器算法是：

（1）采用显式四阶龙格-库塔法积分k步，步长为h，时间从t_n到t_n+k，由x(t_n)得到x(t_n+k)，递推公式如下；

$x\left(t_{n+1}\right)=x\left(t_n\right)+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)$

其中

k₁=hf(x(t_n),t_n)

$k_2=\mathrm{hf}(x\left(t_n\right)+\frac{k_1}{2},t_n+\frac{h}{2})$

$k_3=\mathrm{hf}(x\left(t_n\right)+\frac{k_2}{2},t_n+\frac{h}{2})$

k₄=hf(x(t_n)+k₃,t_n+h)；

（2）再次采用显式四阶龙格-库塔法积分一步，步长仍为h，得到x(t_n+k+1)。

所述的外部积分器算法是：

根据内部积分器得到的x(t_n+k)和x(t_n+k+1)的值，利用下式进行一步外部积分步，步长为内部积分步长的M倍，时间从t_n+k+1到t_n+k+1+M，得到x(t_n+k+1+M)：

x(t_n+k+1+M)=(M+1)x(t_n+k+1)-Mx(t_n+k)。

本发明的一种适于电力***多时间尺度的暂态稳定性投影积分方法，考虑了电力***的多时间尺度特征，采用显式交替求解方法对描述电力***模型的微分-代数方程组进行交替求解，对其中的微分方程利用投影积分算法进行差分求解。本发明的方法不仅适用于传统电力***暂态稳定性仿真，也适用于智能配电网动态仿真问题，具有一般性，在满足仿真数值稳定性和数值精度要求的同时，提升了仿真计算速度，为高效、可靠的电力***仿真程序的开发奠定了良好的基础。本发明的方法计算简单、数值稳定性好、适用于具有刚性问题的电力***暂态稳定性仿真。

附图说明

图1是本发明方法的整体流程图；

图2是低压配电***算例结构图；

图中1：燃料电池；2：第二蓄电池；3：第二光伏电池；4：第一光光伏电池；

5：第一蓄电池；M1：中压母线；L1～L19：低压母线；Load1～Load7：负荷；

图3是L16母线电压仿真结果及局部放大图；

图4是L17母线电压仿真结果及局部放大图；

图5是L17母线频率仿真结果及局部放大图；

图6是L17母线电压投影积分算法仿真结果较DIgSILENT相对误差。

具体实施方式

下面结合实施例和附图对本发明的一种适于电力***多时间尺度的暂态稳定性投影积分方法做出详细说明。

本发明的一种适于电力***多时间尺度的暂态稳定性投影积分方法，属于显式积分方法，考虑了电力***的多时间尺度特征，采用显式交替求解方法对电力***模型进行交替求解，对其中的微分方程利用投影积分算法进行差分求解。不仅适用于传统电力***暂态稳定性仿真，也同样适用于含各种分布式电源及储能装置的智能配电***动态仿真。智能配电***接入了种类繁多的分布式电源和大量的电力电子装置，包括旋转电机和各种静态直流型分布式电源，具有明显的多时间尺度特征。因此，电力***暂态稳定性仿真的多时间尺度问题在数学上可以归结为刚性问题，其对所采用的数值算法的精度和稳定性要求更高。本发明的方法具有一般性，不仅适用于传统电力***暂态稳定性仿真，也适用于含分布式电源及储能的智能配电网动态仿真。本发明在满足仿真数值稳定性和数值精度要求的同时，提升了仿真计算速度，为高效、可靠的电力***仿真程序的开发奠定了良好的基础。

本发明采用显式交替求解算法对电力***数学模型进行求解，关于电力***数学模型中微分方程的求解则采用投影积分算法：首先进行若干小步长的积分运算，对应***的快动态过程；而后根据小步长积分计算结果，进行一个投影步长的计算，与***的慢动态过程对应。其中，小步长积分计算过程称为内部积分器，采用数值稳定性较好的显式四阶龙格-库塔法（explicit four-order Runge-Kutta method）以提高算法的稳定性；大步长积分计算过程称为外部积分器，通过大步长投影计算以提升算法效率。

如图1所示，本发明的一种适于电力***多时间尺度的暂态稳定性投影积分方法，包括如下步骤：

1）输入电力***潮流数据、动态元件数据以及故障及操作设置信息，设置仿真参数；

所述的仿真参数，包括仿真时间T、仿真步长h和投影积分算法参数，其中，所述的投影积分算法参数包括内部积分器的积分步数k和外部积分器步长相对内部积分器步长的倍数M。

所述的内部积分器算法是：

（1）采用显式四阶龙格-库塔法（explicit four-order Runge-Kutta method）积分k步，步长为h，时间从t_n到t_n+k，由x(t_n)得到x(t_n+k)，递推公式如下；

$x\left(t_{n+1}\right)=x\left(t_n\right)+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)$

其中

k₁=hf(x(t_n),t_n)

$k_2=\mathrm{hf}(x\left(t_n\right)+\frac{k_1}{2},t_n+\frac{h}{2})$

$k_3=\mathrm{hf}(x\left(t_n\right)+\frac{k_2}{2},t_n+\frac{h}{2})$

k₄=hf(x(t_n)+k₃,t_n+h)

（2）再次采用显式四阶龙格-库塔法积分一步，步长仍为h，得到x(t_n+k+1)。

所述的外部积分器算法是：

根据内部积分器得到的x(t_n+k)和x(t_n+k+1)的值，利用下式进行一步外部积分步，步长为内部积分步长的M倍，时间从t_n+k+1到t_n+k+1+M，得到x(t_n+k+1+M)：

x(t_n+k+1+M)=(M+1)x(t_n+k+1)-Mx(t_n+k)

2）根据***算例参数和动态数据形成全***节点导纳矩阵，其中包括***节点导纳矩阵和并入动态元件的节点导纳矩阵，并重新进行潮流校核；

3）对***算例的动态元件进行仿真初始化计算；

4）设置仿真计算时间t=0；

5）设置当前投影积分算法内部积分器的积分步数s=1；

6）采用投影积分算法内部积分器对电力***模型积分一个步长，步长为h，然后设置仿真计算时间t=t+h，内部积分器的积分步数s=s+1，内部积分器采用显式交替求解方法对微分-代数方程组进行求解，对其中的微分方程求解则采用显式四阶龙格-库塔法；

所述的电力***模型是下式所示的一个高维非线性且连续自治的微分-代数方程组：

$\left\{\begin{array}{c}\stackrel{\·}{x}=f(x,y)\\ 0=g(x,y)\end{array}\right.$

式中，

为微分方程，用于描述电力***设备动态特征；

为代数方程，用于描述电力***设备之间的电气联系；

为***状态变量，代表同步电机转子转速、电力电子器件控制***和负荷动态参数等；

为代数变量，表征母线电压幅值与相角。

对电力***模型显式差分化后得到下式：

$\left\{\begin{array}{c}x\left(t_{n+1}\right)=x\left(t_n\right)+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)\\ g(x\left(t_{n+1}\right),y\left(t_{n+1}\right))=0\end{array}\right.$

所述的显式交替求解方法对显式差分化后的电力***模型进行求解包括如下步骤：

（1）根据t_n时刻电力***的状态变量x(t_n)和代数变量y(t_n)，计算向量k₁=hf(x(t_n),y(t_n))；

（2）计算向量

然后求解代数方程g(x₁,y₁)=0得出y₁，计算向量k₂=hf(x₁,y₁)；

（3）计算向量

然后求解代数方程g(x₂,y₂)=0得出y₂，计算向量k₃=hf(x₂,y₂)；

（4）计算向量

然后求解代数方程g(x₃,y₃)=0得出y₃，计算向量k₄=hf(x₃,y₃)；

（5）计算t_n+1时刻

然后求解代数方程g(x(t_n+1),y(t_n+1))=0得出y(t_n+1)。

7）根据步骤1）输入的故障及操作设置信息判断***是否发生故障或操作，若发生，则返回步骤5），否则进入下一步骤；

8）判断仿真计算时间t是否大于设置时间T，若计算时间t大于设置时间T，则仿真结束，否则进入下一步骤；

9）判断投影积分算法内部积分步数s是否大于设置步数k+1，若不大于，则返回步骤6），否则进入下一步骤；

10）根据步骤1）输入的故障及操作设置信息判断未来Mh时间内是否存在故障或操作，若存在，则返回步骤5），否则进入下一步骤；

11）判断未来Mh时间内仿真计算时间t是否会超过设置时间T，若超过，则返回步骤5），否则进入下一步骤；

12）根据内部积分器得到的x(t_n+k)和x(t_n+k+1)的值进行一步外部积分步，步长为内部积分器积分步长的M倍，利用式x(t_n+k+1+M)=(M+1)x(t_n+k+1)-Mx(t_n+k)得到x(t_n+k+1+M)，然后求解代数方程g(x(t_n+k+1+M),y(t_n+k+1+M))=0得到y(t_n+k+1+M)，设置t=t+Mh；

13）判断此时仿真计算时间t是否达到设置的仿真时间T，若达到，则仿真结束，否则返回步骤5），依此反复进行直至仿真结束。

下面给出具体实例：

数字仿真和电网计算程序（DIgSILENT PowerFactory）是德国DIgSLENTGmbH公司开发的一款商业电力***仿真软件。本实例以C++编程语言环境为基础，实现了面向电力***暂态稳定性仿真的投影积分仿真方法，通过含分布式电源的低压配电***算例（附图2）对本发明的方法进行测试验证，并与商业软件DIgSILENT PowerFactory进行比较。

低压配电***算例电压等级为400V，主馈线通过0.4/10kV变压器接至中压母线M1处，变压器采用常用的DYn11联结方式，低压侧设有电容无功补偿***，主馈线节点间距为50m。考虑到配电***的实际运行特点，测试算例中既有三相对称和不对称负荷，也有单相负荷。另外，算例中接入了多种类型的分布式电源，包括：具备最大功率跟踪控制的光伏发电***，燃料电池发电***和蓄电池储能***，各分布式电源控制方式、接入容量及输出功率如表1所示。

表1分布式电源控制方式、接入容量及输出功率

采用本发明提出的投影积分算法对测试***算例进行动态仿真计算，设置仿真时间为5s，仿真步长为0.3ms，2.0s时刻L1母线发生三相短路故障，2.1s时刻故障清除。关于算法参数的选取本发明以k=6，M=2，k=8，M=2和k=8，M=3为例对该算法进行测试，在具体实现时，算法参数可根据实际应用情况，在满足数值稳定性和数值精度的条件下进行任意取值，本发明的实施对此不做限制。执行仿真计算的硬件平台为Intel(R)Core(TM)i5-3470CPU3.20GHz，4GB RAM的PC机；软件环境为32位Windows7操作***。

将本发明提出的一种适于电力***多时间尺度的暂态稳定性投影积分方法与步长取0.1ms时DIgSILENT PowerFactory的仿真结果进行比较，仿真结果如附图3至附图5所示，附图6为不同参数投影算法较DIgSILENT仿真结果的相对误差。比较分析结果表明投影积分算法具有较好的精度，其仿真结果与DIgSILENT吻合，而显式交替求解方法由于存在一定的交接误差，导致计算结果存在一定的误差，且随着k值的减小或M值的增大，投影积分算法与DIgSILENT仿真结果的偏差逐渐增大，但最大误差维持在0.15%以内，显示出该算法具有良好的数值精度。

本发明提出的一种适于电力***多时间尺度的暂态稳定性投影积分方法可以实现对其内部积分算法，即显式四阶龙格-库塔法仿真速度的提升，将选取不同参数的投影积分算法的仿真计算时间与同步长显式四阶龙格-库塔法以及商业软件DIgSILENT定步长仿真进行比较，如表2所示。

表2算法性能比较

表2给出了选取不同参数时投影积分算法与同步长显式四阶龙格-库塔法以及商业软件DIgSILENT定步长仿真计算时间的对比。可以看出，投影积分算法的仿真计算时间明显小于传统显式四阶龙格-库塔法和DIgSILENT，且随着k值的减小或M值的增大，投影算法仿真速度的提升越明显，相比显式四阶龙格-库塔法其最大加速比可达47.9%。由此可见，本发明提出的一种适于电力***多时间尺度的暂态稳定性投影积分方法可以实现仿真计算效率的提升。

Claims (4)

1.一种适于电力***多时间尺度的暂态稳定性投影积分方法，其特征在于，包括如下步骤：

1）输入电力***潮流数据、动态元件数据以及故障及操作设置信息，设置仿真参数；

2）根据***算例参数和动态数据形成全***节点导纳矩阵，其中包括***节点导纳矩阵和并入动态元件的节点导纳矩阵，并重新进行潮流校核；

3）对***算例的动态元件进行仿真初始化计算；

4）设置仿真计算时间t=0；

5）设置当前投影积分算法内部积分器的积分步数s=1；

6）采用投影积分算法内部积分器对电力***模型积分一个步长，步长为h，然后设置仿真计算时间t=t+h，内部积分器的积分步数s=s+1，内部积分器采用显式交替求解方法对微分-代数方程组进行求解，对其中的微分方程求解则采用显式四阶龙格-库塔法；

7）根据步骤1）输入的故障及操作设置信息判断***是否发生故障或操作，若发生，则返回步骤5），否则进入下一步骤；

8）判断仿真计算时间t是否大于设置时间T，若计算时间t大于设置时间T，则仿真结束，否则进入下一步骤；

9）判断投影积分算法内部积分步数s是否大于设置步数k+1，若不大于，则返回步骤6），否则进入下一步骤；

10）根据步骤1）输入的故障及操作设置信息判断未来Mh时间内是否存在故障或操作，若存在，则返回步骤5），否则进入下一步骤；

11）判断未来Mh时间内仿真计算时间t是否会超过设置时间T，若超过，则返回步骤5），否则进入下一步骤；

12）根据内部积分器得到的x(t_n+k)和x(t_n+k+1)的值进行一步外部积分步，步长为内部积分器积分步长的M倍，利用式x(t_n+k+1+M)=(M+1)x(t_n+k+1)-Mx(t_n+k)得到x(t_n+k+1+M)，然后求解代数方程g(x(t_n+k+1+M),y(t_n+k+1+M))=0得到y(t_n+k+1+M)，设置t=t+Mh；

13）判断此时仿真计算时间t是否达到设置的仿真时间T，若达到，则仿真结束，否则返回步骤5），依此反复进行直至仿真结束。

2.根据权利要求1所述的一种适于电力***多时间尺度的暂态稳定性投影积分方法，其特征在于，步骤1）中所述的仿真参数，包括仿真时间T、仿真步长h和投影积分算法参数，其中，所述的投影积分算法参数包括内部积分器的积分步数k和外部积分器步长相对内部积分器步长的倍数M。

3.根据权利要求2所述的一种适于电力***多时间尺度的暂态稳定性投影积分方法，其特征在于，所述的内部积分器算法是：

（1）采用显式四阶龙格-库塔法积分k步，步长为h，时间从t_n到t_n+k，由x(t_n)得到x(t_n+k)，递推公式如下；

$x\left(t_{n+1}\right)=x\left(t_n\right)+\frac{1}{6}(k_1+{2k}_2+{2k}_3+k_4)$

其中

k₁=hf(x(t_n),t_n)

$k_2=\mathrm{hf}(x\left(t_n\right)+\frac{k_1}{2},t_n+\frac{h}{2})$

$k_3=\mathrm{hf}(x\left(t_n\right)+\frac{k_2}{2},t_n+\frac{h}{2})$

k₄=hf(x(t_n)+k₃,t_n+h)；

（2）再次采用显式四阶龙格-库塔法积分一步，步长仍为h，得到x(t_n+k+1)。

4.根据权利要求2所述的一种适于电力***多时间尺度的暂态稳定性投影积分方法，其特征在于，所述的外部积分器算法是：

根据内部积分器得到的x(t_n+k)和x(t_n+k+1)的值，利用下式进行一步外部积分步，步长为内部积分步长的M倍，时间从t_n+k+1到t_n+k+1+M，得到x(t_n+k+1+M)：

x(t_n+k+1+M)=(M+1)x(t_n+k+1)-Mx(t_n+k)。

Images