CN104679937A - 一种适于隐式投影算法的误差估计及参数自适应调节方法 - Google Patents

Abstract

一种适于隐式投影算法的误差估计及参数自适应调节方法，充分利用了隐式投影算法数值稳定域与其算法参数M基本不相关的数值稳定特性，针对具有刚性特征的有源配电系统稳定性仿真模型，同时考虑仿真过程中可能出现的故障与开关操作等事件，分别根据事前误差估计方法和事后误差估计方法对隐式投影算法每个投影积分步的计算误差进行事前估计和事后估计，并通过算法参数的自适应变化实现对投影算法仿真误差的控制，较固定参数的隐式投影算法具有明显优势，尤其适于在大扰动下具有显著多时间尺度特征的大规模有源配电系统稳定性仿真问题的快速、准确求解，大大提高了投影积分算法的适用性，为高效、可靠的有源配电系统仿真程序的开发奠定了坚实的基础。

Description

一种适于隐式投影算法的误差估计及参数自适应调节方法

技术领域

本发明涉及一种有源配电系统稳定性仿真方法。特别是涉及一种适用于基于隐式投影积分算法的含大量分布式电源的有源配电系统稳定性仿真应用的适于隐式投影算法的误差估计及参数自适应调节方法。

背景技术

近年来，随着分布式发电(distributed generation,DG)、储能及需求侧响应等技术在配电层面的广泛应用，传统的单向无源配电系统正逐渐过渡为由于分布式电源接入而具有双向潮流的有源配电系统。有源配电系统是具备组合控制各种分布式能源(DER，如DG、可控负荷、储能等)能力的复杂配电系统，其中涉及静态直流环节、交流电机、电力电子装置及其控制系统等多个部分，且各环节动态特性的时间尺度跨度极大，使得有源配电系统具有较为显著的多时间尺度特征，对配电系统的规划与运行提出了新的挑战。DG的大规模接入会导致配电网的动态响应特性发生变化进而影响整个电力系统的动态特性，特别是受大扰动时的动态特性。在系统层面，相关问题的分析与研究往往无法直接在实际系统上进行试验，因此必须采用有效的数字仿真工具作为研究手段，一种可靠、高效的数字仿真工具对有源配电系统动态特性的分析具有重要作用。

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

基于详细建模的电磁暂态仿真是目前有源配电系统暂态分析中较常用的仿真方法，适合观察并分析系统各元件详细的动态特性变化。然而，随着DG接入的逐渐增多，仿真系统规模逐渐增大，其计算负担也会随之大幅增加，计算效率受到严重影响，进而限制了电磁暂态仿真对含大规模分布式电源接入的有源配电系统运行分析与研究的应用。基于准稳态建模的有源配电系统稳定性仿真着重于分析系统运行时其工频电气量在系统扰动下(开关操作、故障、分布式电源及负荷波动等)的动态响应特性，仿真计算效率明显高于电磁暂态仿真，广泛应用于大规模有源配电系统的动态分析和相关运行策略的验证等方面。

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

有源配电系统全系统稳定性仿真模型可通过一个高维微分-代数方程组(differential-algebraic equation,DAE)描述，如式(1)所示。

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

式中，为微分方程，为代数方程，为系统状态变量，代表电机转子转速、电力电子器件控制系统和负荷动态参数等，为代数变量，表征母线电压幅值与相角。有源配电系统数学模型的求解一般通过特定的数值算法来实现。有源配电系统的多时间尺度特征在数学上体现为“刚性”特征，因此，有源配电系统稳定性仿真在数学上可以归结为求解一个刚性微分-代数方程组的初值问题，其对所采用的数值算法的精度和数值稳定性要求更高。

隐式投影积分方法是针对具有刚性特征有源配电系统提出的一种固定参数、固定步长的稳定性仿真方法。隐式投影方法基于交替求解思想对式(1)所示的DAE方程组进行求解，并采用隐式投影积分算法对其中的常微分方程(ordinary differential equation,ODE)进行求解，其投影积分步的积分步骤为：首先以步长h进行若干步的小步长积分计算，而后根据小步长的计算结果，基于隐式预测-校正方法，以步长Mh进行一个大步长投影积分计算。其中，小步长积分计算过程称为内部积分器，须选取具有二阶以上精度的显式数值积分方法；大步长投影积分过程称为外部积分器。

隐式投影积分算法是一种二阶精度算法，其数值稳定性几乎不受算法参数—外部积分器积分步长倍数M值的影响，即随着M的变化，隐式投影算法的数值稳定域基本不发生变化，同时，在某些情况下，隐式投影算法较传统隐式梯形算法具有更好的算法性能。然而，由于隐式投影积分算法采用固定步长、固定参数进行积分计算，如果选取较大的M值进行积分计算，虽然投影算法的计算效率会得到显著提升，但当仿真过程中系统出现诸如短路故障、开关开/断、负荷投/切等操作或事件引起的快动态过程时，系统各变量单位时间的波动较大，外部积分器的大步长积分计算可能会产生较大的仿真误差乃至无法准确反映出此时系统的快动态特性；而若一味地采用较小的M值，在系统处于稳定运行或慢动态过程时，各变量单位时间的波动相对较小，会产生诸多不必要的小步长计算过程，继而严重影响仿真计算效率。

可见，提出一种基于误差估计的隐式投影算法参数自适应调节方法，实现含大规模分布式电源接入的有源配电系统快速、准确的仿真计算与动态分析，对隐式投影算法的进一步发展与应用具有重要意义。

发明内容

本发明所要解决的技术问题是，提供一种能够实现仿真过程中隐式投影算法参数随估计误差的自适应调节的适于隐式投影算法的误差估计及参数自适应调节方法。

本发明所采用的技术方案是：一种适于隐式投影算法的误差估计及参数自适应调节方法，包括如下步骤：

1)读取电力系统拓扑连接关系、元件参数和仿真计算参数，其中仿真计算参数包括仿真终止时间T，仿真步长h，隐式投影算法内部积分器的积分步数k，隐式投影算法外部积分器步长倍数的初始值M_ini和最大值M_max，最大允许误差并设置仿真故障及操作事件；

2)设置当前隐式投影算法内部积分器的积分步数s＝1，s为正整数；

3)设置仿真时间t＝t+h，采用隐式投影算法内部积分器计算一个步长得到电力系统该时刻的状态变量x_n+s和代数变量y_n+s，并设置s＝s+1；

4)根据步骤1)设置的仿真故障及操作事件，判断电力系统当前时刻是否存在故障或操作，若存在，即故障或操作事件的发生时间T_event＝t，返回步骤2)，否则进行下一步骤；

5)判断隐式投影算法内部积分器的积分步数s是否大于k+1，是则进入下一步骤，若s≤k+1，则返回步骤3)；

6)设当前隐式投影算法外部积分器以M_n-1h作为积分步长，其中，M_n-1为上一步的隐式投影算法外部积分器步长倍数，利用事前误差估计方法对隐式投影算法内部积分器积分k+1步以及隐式投影算法外部积分器以步长M_n-1h积分一步的局部截断误差进行估计，得到事前估计误差

7)根据估计误差，估计当前隐式投影算法外部积分器步长倍数M_n；

8)根据步骤1)设置的故障及操作事件，判断t～t+M_nh时间内是否发生故障或操作，若发生，即t<T_event<t+M_nh，则进入步骤10)，否则进行下一步骤；

9)设置隐式投影算法外部积分器积分步长H＝M_nh，设置仿真时间t＝t+H，利用隐式投影算法外部积分器得到系统该时刻的状态变量和代数变量然后转入步骤11)；

10)设置隐式投影算法外部积分器外部积分步长H＝T_event-t，设置仿真时间t＝T_event，利用隐式投影算法外部积分器得到故障或操作发生前系统的状态变量和代数变量

11)通过事后误差估计方法对隐式投影算法外部积分器的事后估计误差进行估计；

12)将与进行比较，若则拒绝本次隐式投影算法外部积分器的积分计算，设置仿真时间t＝t-H，返回步骤7)，否则进入下一步骤；

13)判断仿真时间t是否达到仿真终止时间T，若t＝T，则仿真结束，否则返回步骤2)，依据步骤2)至13)反复进行直至仿真结束。

步骤1)所述的隐式投影算法内部积分器采用基于显式改进欧拉法的微分代数方程交替求解方法，积分步数k要求大于等于3，隐式投影算法外部积分器采用基于隐式预测-校正方法的微分代数方程交替求解方法进行计算，隐式投影算法外部积分器步长倍数的初始值M_ini值为1。

步骤6)所述的事前误差估计方法是：

设x_n为准确值，经隐式投影算法内部积分器积分k+1步及隐式投影算法外部积分器以步长M_n-1h积分一步后的局部截断误差估计值，通过下式得到

${\hat{T}}_{n+k+1+M_{n-1}}^{<\mathrm{a\; priori}>}\&ap;\left|\right|[\frac{1}{6}k(k+1)(k+2)+\frac{{(k+1)}^2}{2}{M}_{n-1}+\frac{k+1}{2}{M}_{n-1}^2+\frac{1}{4}{M}_{n-1}^3]h^3{x}_n^{\left(3\right)}\left|\right|$

式中，状态变量的三阶微分项通过拉格朗日插值方法得到，||x||为欧几里得距离，分别如下式所示

$x_n^{\left(3\right)}\&ap;-\frac{1}{h^3}(x_n-{3x}_{n+1}+{3x}_{n+2}-x_{n+3})$

$\left|\right|x\left|\right|={\left(x^{T}x\right)}^{\frac{1}{2}}$

当对第一步隐式投影算法外部积分器的事前误差进行估计时，上一步的隐式投影算法外部积分器步长倍数M_n-1，即M₀取M_ini。

步骤7)所述的当前隐式投影算法外部积分器步长倍数M_n利用估计误差，并由下式估计得到

$M_n={\left|\frac{{\hat{T}}_{\mathrm{tol}}}{\hat{T}}\right|}^{\frac{1}{3}}{M}_{\mathrm{old}}$

式中，为事前估计误差或事后估计误差M_old为事前误差估计时采用的M_n-1或经事后估计误差判断后被拒绝的隐式投影算法外部积分器的外部积分步长倍数M_{n_rejected}。

步骤11)所述的事后误差估计方法是：

隐式投影算法外部积分器采用隐式预测-校正方法对微分方程组进行积分求解，经过隐式投影算法外部积分器的预测步计算得到的系统状态变量记为经过隐式投影算法外部积分器的校正步计算得到的状态变量记为则隐式投影算法外部积分器的事后估计误差由下式得到

${\hat{T}}_{n+k+1+M_n}^{<\mathrm{a\; posteriori}>}\&ap;\left|\right|x_{n+k+1+M_n}^{<1>}-x_{n+k+1+M_n}^{<0>}\left|\right|.$

本发明的一种适于隐式投影算法的误差估计及参数自适应调节方法，充分利用了隐式投影算法数值稳定域与其算法参数M基本不相关的数值稳定特性，分别根据事前误差估计方法和事后误差估计方法对投影算法每个投影积分步的计算误差进行事前估计和事后估计，并通过算法参数的自适应变化实现对投影算法仿真误差的控制，较固定参数的隐式投影算法具有明显优势，尤其适于在大扰动下具有显著多时间尺度特征的大规模有源配电系统稳定性仿真问题的快速、准确求解，大大提高了投影积分算法的适用性，为高效、可靠的有源配电系统仿真程序的开发奠定了坚实的基础。

附图说明

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

图2是不同算法参数下隐式投影算法和显式投影算法的数值稳定域；

图3是图2中A的局部放大示意图；

图4是IEEE123节点有源配电系统算例结构图；

图5是61节点处的光伏电池并网电压仿真结果及局部放大图；

图6是自适应参数隐式投影算法M值和事前估计误差变化曲线；

图7是自适应参数隐式投影算法与固定参数隐式投影算法精度比较(对数坐标系)；

图8是低压有源配电系统算例结构图；

图中1：第一蓄电池；2：第一光伏电池；3：第二光伏电池；4：第二蓄电池；M1：中压母线；S1：开关；L1～L19：低压母线；Load1～Load7：负荷；

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

图10是第二光伏电池有功功率输出仿真结果及局部放大图；

图11是自适应参数隐式投影算法M值和事前估计误差变化曲线。

具体实施方式

下面结合实施例和附图对本发明的一种适于隐式投影算法的误差估计及参数自适应调节方法做出详细说明。

本发明的一种适于隐式投影算法的误差估计及参数自适应调节方法，属于事前-事后误差估计与控制方法。隐式投影积分算法是针对刚性有源配电系统稳定性仿真问题而提出的一种具有二阶精度的数值积分算法，其数值稳定性几乎不受其算法参数——外部积分器积分步长倍数M值的影响，且在某些特定情况下，隐式投影算法较传统隐式梯形算法具有更好的算法性能。然而，由于隐式投影积分算法采用固定步长、固定参数进行积分计算，如果选取较大的M值进行积分计算，虽然投影算法的计算效率会得到显著提升，但当仿真过程中系统出现诸如短路故障、开关开/断、负荷投/切等操作或事件引起的快动态过程时，系统各变量单位时间的波动较大，外部积分器的大步长积分计算可能会产生较大的仿真误差乃至无法准确反映出此时系统的快动态特性；而若一味地采用较小的M值，在系统处于稳定运行或慢动态过程时，各变量单位时间的波动相对较小，会产生诸多不必要的小步长计算过程，继而严重影响仿真计算效率。因此，一种符合隐式投影算法特点的误差估计和参数自适应调节方法对隐式投影算法的进一步发展与应用具有重要意义。

本发明提出的一种适于隐式投影算法的误差估计及参数自适应调节方法，充分利用了隐式投影算法数值稳定域与其算法参数M基本不相关的数值稳定特性，分别根据事前误差估计方法和事后误差估计方法对投影算法每个投影积分步的计算误差进行事前估计和事后估计，并通过算法参数的自适应变化实现对投影算法仿真误差的控制，较固定参数的隐式投影算法具有明显优势，尤其适于在大扰动下具有显著多时间尺度特征的大规模有源配电系统动态仿真问题的快速、准确求解，大大提高了投影积分算法的适用性，为高效、可靠的有源配电系统仿真程序的开发奠定了坚实的基础。

本发明基于事前误差估计和事后误差估计方法，对隐式投影算法进行误差估计进而实现算法参数的自适应调节，其基本思想为：首先隐式投影算法内部积分器计算完成后，在进入外部积分器之前，利用上一步投影积分步的外部积分器步长倍数M_n-1，根据事前误差估计方法对本次投影积分步的局部截断误差进行事前估计；然后，基于事前估计误差预测当前外部积分器的积分步长倍数M_n，并进行外部积分计算；最后，基于外部积分计算结果，利用事后误差估计方法对M_n的预测值进行验证，若满足误差允许条件，则进入下一个投影积分步计算，若不满足，则拒绝本次外部积分步的计算，基于事后估计误差估计算法参数M_n并重新进行外部积分计算，依此反复进行直至满足误差允许条件。

如附图1所示，本发明的一种适于隐式投影算法的误差估计及参数自适应调节方法，包括如下步骤：

1)读取电力系统拓扑连接关系、元件参数和仿真计算参数，其中仿真计算参数包括仿真终止时间T，仿真步长h，隐式投影算法内部积分器的积分步数k，隐式投影算法外部积分器步长倍数的初始值M_ini和最大值M_max，最大允许误差并设置仿真故障及操作事件；

所述的隐式投影算法内部积分器是采用申请号为CN201410425443.5所公开的隐式投影算法内部积分器，采用基于显式改进欧拉法的微分代数方程交替求解方法，积分步数k要求大于等于3，隐式投影算法外部积分器采用基于隐式预测-校正方法的微分代数方程交替求解方法进行计算，隐式投影算法外部积分器步长倍数的初始值M_ini值为1。

所述的隐式投影算法外部积分器是采用申请号为CN201410425443.5所公开的隐式投影算法外部积分器。

2)设置当前隐式投影算法内部积分器的积分步数s＝1，s为正整数；

3)设置仿真时间t＝t+h，采用隐式投影算法内部积分器计算一个步长得到电力系统该时刻的状态变量x_n+s和代数变量y_n+s，并设置s＝s+1；

4)根据步骤1)设置的仿真故障及操作事件，判断电力系统当前时刻是否存在故障或操作，若存在，即故障或操作事件的发生时间T_event＝t，返回步骤2)，否则进行下一步骤；

5)判断隐式投影算法内部积分器的积分步数s是否大于k+1，是则进入下一步骤，若s≤k+1，则返回步骤3)；

6)设当前隐式投影算法外部积分器以M_n-1h作为积分步长，其中，M_n-1为上一步的隐式投影算法外部积分器步长倍数，利用事前误差估计方法对隐式投影算法内部积分器积分k+1步以及隐式投影算法外部积分器以步长M_n-1h积分一步的局部截断误差进行估计，得到事前估计误差

所述的事前误差估计方法是：

设x_n为准确值，经隐式投影算法内部积分器积分k+1步及隐式投影算法外部积分器以步长M_n-1h积分一步后的局部截断误差估计值，通过下式得到

${\hat{T}}_{n+k+1+M_{n-1}}^{<\mathrm{apriori}>}\&ap;\left|\right|[\frac{1}{6}k(k+1)(k+2)+\frac{{(k+1)}^2}{2}{M}_{n-1}+\frac{k+1}{2}{M}_{n-1}^2+\frac{1}{4}{M}_{n-1}^3]h^3{x}_n^{\left(3\right)}\left|\right|$

式中，状态变量的三阶微分项通过拉格朗日(Lagrange)插值方法得到，||x||为欧几里得距离(Euclidean Length)，分别如下式所示

$x_n^{\left(3\right)}\&ap;-\frac{1}{h^3}(x_n-{3x}_{n+1}+{3x}_{n+2}-x_{n+3})$

$\left|\right|x\left|\right|={\left(x^{T}x\right)}^{\frac{1}{2}}$

当对第一步隐式投影算法外部积分器的事前误差进行估计时，上一步的隐式投影算法外部积分器步长倍数M_n-1，即M₀取M_ini。

7)根据估计误差，估计当前隐式投影算法外部积分器步长倍数M_n；

所述的当前隐式投影算法外部积分器步长倍数M_n利用估计误差，并由下式估计得到

$M_n={\left|\frac{{\hat{T}}_{\mathrm{tol}}}{\hat{T}}\right|}^{\frac{1}{3}}{M}_{\mathrm{old}}$

式中，为事前估计误差或事后估计误差M_old为事前误差估计时采用的M_n-1或经事后估计误差判断后被拒绝的隐式投影算法外部积分器的外部积分步长倍数M_{n_rejected}。

8)根据步骤1)设置的故障及操作事件，判断t～t+M_nh时间内是否发生故障或操作，若发生，即则t<T_event<t+M_nh，进入步骤10)，否则进行下一步骤；

9)设置隐式投影算法外部积分器积分步长H＝M_nh，设置仿真时间t＝t+H，利用隐式投影算法外部积分器得到电力系统该时刻的状态变量和代数变量然后进入步骤11)；

10)设置隐式投影算法外部积分器外部积分步长H＝T_event-t，设置仿真时间t＝T_event，利用隐式投影算法外部积分器得到故障或操作发生前系统的状态变量和代数变量

11)通过事后误差估计方法对隐式投影算法外部积分器的事后估计误差进行估计；

所述的事后误差估计方法是：

隐式投影算法外部积分器采用隐式预测-校正方法对微分方程组进行积分求解，经过隐式投影算法外部积分器的预测步计算得到的系统状态变量记为经过隐式投影算法外部积分器的校正步计算得到的状态变量记为则隐式投影算法外部积分器的事后估计误差由下式得到：

${\hat{T}}_{n+k+1+M_n}^{<\mathrm{aposteriori}>}\&ap;\left|\right|x_{n+k+1+M_n}^{<1>}-x_{n+k+1+M_n}^{<0>}\left|\right|.$

12)将与进行比较，若则拒绝本次隐式投影算法外部积分器的积分计算，设置仿真时间t＝t-H，返回步骤7)，否则进入下一步骤；

13)判断仿真时间t是否达到仿真终止时间T，若t＝T，则仿真结束，否则返回步骤2)，依据步骤2)至13)反复进行直至仿真结束。

下面给出具体实例：

本实例中隐式投影算法内部积分器取显式改进欧拉法时，根据绝对稳定条件，分别得到不同参数下隐式投影算法在hλ平面中的数值稳定域，如附图2和附图3所示。可以看出，隐式投影积分算法的数值稳定域与改进欧拉法的稳定域基本一致，且不随外部积分器步长倍数M的改变而发生明显变化，因此可以认为隐式投影算法的数值稳定域与其算法参数M值基本不相关。

本实例在C++编程环境中实现了本发明提出的一种适于隐式投影算法的误差估计及参数自适应调节方法，通过将自适应参数隐式投影算法与固定参数隐式投影算法的仿真结果与计算性能进行比较以验证了本发明方法的正确性和有效性，并以小步长显式4阶龙格-库塔法(以下简称RK4算法)为基准，测试了本发明方法的误差控制效率。执行仿真测试的硬件平台为Intel(R)Core(TM)i5-3470CPU3.20GHz，4GB RAM的PC机；软件环境为32位Windows7操作系统。

首先，本实例采用IEEE 123节点测试算例(如附图4所示)对本发明方法进行验证。IEEE123节点算例描述了一个结构复杂的辐射状配电网络，共有123个节点，电压等级为4.16kV，其内部考虑了多种形式的负荷，并在节点150处与外部网络相连。本实例在附图4中虚线框内的节点处共接入50个容量为30kWp，有功功率输出为20.4kW的光伏发电系统。采用自适应参数隐式投影算法对测试算例进行稳定性仿真计算，设置仿真时间为5s，仿真步长为0.5ms，2.5s时刻系统1号母线发生三相短路故障，0.1s后故障清除。

将本发明的一种适于隐式投影算法的误差估计及参数自适应调节方法应用于隐式投影算法中进行仿真计算，并将自适应参数隐式投影算法与固定参数隐式投影算法的仿真结果进行比较，算法参数设置如表1所示。

61号节点母线电压仿真结果如附图5所示，实线、虚线和点划线为DIgSILENT，自适应参数隐式投影算法以及固定参数隐式投影算法的仿真计算结果。由图可知，自适应参数隐式投影算法的仿真结果与DIgSILENT仿真结果基本一致，且在系统发生故障后的暂态运行过程中，自适应参数算法与固定参数算法相比可以更准确地反映系统的动态运行特性。

表1 隐式投影算法参数设置

自适应参数投影算法M_n值和事前估计误差在仿真过程中的变化情况如附图6所示，可以看出，在仿真开始时刻到2.5s期间，系统处于稳态运行状态，估计误差几乎为零，M_n值迅速增大至最大值M_max；2.5s～2.6s系统发生短路故障而进入暂态运行过程，估计误差急剧增大，M_n值迅速减小至1以满足误差允许条件；故障清除后，系统经过短暂的暂态运行后逐渐恢复稳定，估计误差在故障恢复期间存在小幅振荡并逐渐减小，因此M_n值也经过小幅摆动缓慢上升，并在系统恢复稳态运行时达到最大值M_max。

为验证本发明方法的仿真计算精度，以步长为0.1ms的RK4算法为基准，在对数坐标系中比较自适应参数隐式投影算法和固定参数隐式投影算法相对于RK4算法的绝对误差，如附图7所示。图中，实线为固定参数隐式投影算法的仿真误差，虚线为自适应参数隐式投影算法的仿真误差，在系统稳态运行过程中，自适应参数投影算法的M_n值取值较大，仿真误差大于固定参数投影算法，但此时系统处于稳态运行，这种误差差异基本不会影响系统动态运行特性的分析；系统发生故障而进入暂态运行过程时，自适应参数投影算法的仿真误差小于固定参数算法，可以更准确地反映系统的动态运行特性，满足有源配电系统动态仿真分析的实际需求。

本实例以步长取0.5ms的RK4算法为基准，选取不同的仿真步长和算法参数，比较自适应参数隐式投影算法和固定参数隐式投影算法的仿真计算效率，结果如表2所示。

表2 算法性能比较

由表2可知，对于固定参数隐式投影算法，与传统RK4算法相比，其计算效率得到了一定程度的提升。然而，在系统故障运行的过程中，若M值取值较大，固定参数投影算法由于误差快速累积而最终数值不收敛；而对于自适应参数隐式投影算法，系统故障期间在满足仿真精度的同时，在系统稳态运行过程中M_n达到最大值运行，大大缩短了仿真时间。另外，表2括号中代表了因事后估计误差不满足误差允许条件而进行重算的次数，分别在附图6中“×”标记的位置进行重算。可以看出，外部积分器的重算主要集中在故障恢复M_n逐渐上升的过程中，这是由于在故障恢复过程中M_n值上升过快，使得仿真计算误差超过误差允许值而引起。因此，事后误差估计方法对M_n值的事前估计具有校正作用。综上所述，本发明提出的一种适于隐式投影算法的误差估计及参数自适应调节方法可以快速准确地模拟系统发生短路故障的运行特性，其有效性得到验证。

另外，为验证本发明方法对具有变结构特征的有源配电系统稳定性仿真问题的适应性，本实例基于一个含分布式电源的低压有源配电系统算例(如附图8所示)对本发明方法进行测试。低压有源配电系统算例电压等级为400V，主馈线通过0.4/10kV变压器接至中压母线M1处，变压器采用常用的DYn11联结方式，低压侧设有无功补偿电容，主馈线节点间距为50m，采用三相对称线路与负荷。另外，算例中接入了多种类型的分布式电源，包括：具备最大功率跟踪控制的光伏发电系统和蓄电池储能系统，各分布式电源控制方式、接入容量及有功功率输出如表3所示。设置仿真时间为9s，仿真步长为0.5ms。2.0s时刻低压有源配电系统开关S1断开，系统由并网运行模式切换至孤岛运行模式；4.7s时刻S1开关闭合，系统由孤岛运行模式切换至并网运行模式。

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

将可变参数的隐式投影积分算法与固定参数的隐式投影算法的仿真结果进行比较，算法参数与表1相同，L16母线电压及第二蓄电池有功功率输出的仿真结果及局部放大图如附图9和附图10所示，与前文结果相同，自适应参数隐式投影算法仿真结果与DIgSILENT的仿真结果基本一致。仿真过程中自适应参数投影算法事前估计误差与M_n值的变化情况如附图11所示，图中，实线表示M_n，虚线表示事前估计误差。可以看出，在系统稳态运行过程中，事前估计误差接近于零，M_n逐渐增大到最大值M_max；2s时刻开关S1断开，系统结构发生改变，各物理量发生剧烈变化，估计误差迅速增大，M_n迅速下降至1，并在系统暂态运行过程中一直以较小值运行；4s时刻开关S1闭合，在闭合瞬间，估计误差由于系统结构发生变化而增大，M_n也随之迅速降为最小值，随着系统逐渐恢复到稳态运行，估计误差减小为较小值，M_n逐渐恢复到最大值M_max直至仿真结束。

以步长取0.5ms的RK4算法为基准，选取不同的仿真步长和算法参数，分别比较自适应参数隐式投影算法和固定参数隐式投影算法的计算效率，比较结果如表4所示。

表4 算法性能比较

与表2结论相同，固定参数投影算法虽然可以一定程度上提高计算效率，但随着M值的增大，在系统孤岛运行过程中算法无法维持收敛，且k值的增大依然无法改善固定参数算法的数值收敛情况。而自适应参数隐式投影算法可以在M_max取30时依旧保持数值收敛，进一步节省了计算用时。

综上所述，本发明提出的一种适于隐式投影算法的误差估计及参数自适应调节方法，充分利用了隐式投影算法数值稳定域与其算法参数基本不相关的数值稳定特性，基于误差估计方法实现了隐式投影算法的自适应调节。自适应参数隐式投影算法较固定参数隐式投影算法在速度和数值收敛性方面具有明显优势，尤其适于在大扰动下具有显著多时间尺度特征的大规模有源配电系统动态仿真问题的快速、准确求解，大大提高了投影积分算法的适用性，为高效、可靠的有源配电系统仿真程序的开发奠定了坚实的基础。

Claims (5)

1.一种适于隐式投影算法的误差估计及参数自适应调节方法，其特征在于，包括如下步骤：

1)读取电力系统拓扑连接关系、元件参数和仿真计算参数，其中仿真计算参数包括仿真终止时间T，仿真步长h，隐式投影算法内部积分器的积分步数k，隐式投影算法外部积分器步长倍数的初始值M_ini和最大值M_max，最大允许误差并设置仿真故障及操作事件；

2)设置当前隐式投影算法内部积分器的积分步数s＝1，s为正整数；

3)设置仿真时间t＝t+h，采用隐式投影算法内部积分器计算一个步长得到电力系统该时刻的状态变量x_n+s和代数变量y_n+s，并设置s＝s+1；

4)根据步骤1)设置的仿真故障及操作事件，判断电力系统当前时刻是否存在故障或操作，若存在，即故障或操作事件的发生时间T_event＝t，返回步骤2)，否则进行下一步骤；

5)判断隐式投影算法内部积分器的积分步数s是否大于k+1，是则进入下一步骤，若s≤k+1，则返回步骤3)；

6)设当前隐式投影算法外部积分器以M_n-1h作为积分步长，其中，M_n-1为上一步的隐式投影算法外部积分器步长倍数，利用事前误差估计方法对隐式投影算法内部积分器积分k+1步以及隐式投影算法外部积分器以步长M_n-1h积分一步的局部截断误差进行估计，得到事前估计误差

7)根据估计误差，估计当前隐式投影算法外部积分器步长倍数M_n；

8)根据步骤1)设置的故障及操作事件，判断t～t+M_nh时间内是否发生故障或操作，若发生，即t<T_event<t+M_nh，则进入步骤10)，否则进行下一步骤；

9)设置隐式投影算法外部积分器积分步长H＝M_nh，设置仿真时间t＝t+H，利用隐式投影算法外部积分器得到系统该时刻的状态变量和代数变量然后转入步骤11)；

10)设置隐式投影算法外部积分器外部积分步长H＝T_event-t，设置仿真时间t＝T_event，利用隐式投影算法外部积分器得到故障或操作发生前系统的状态变量和代数变量

11)通过事后误差估计方法对隐式投影算法外部积分器的事后估计误差进行估计；

12)将与进行比较，若 ${\hat{T}}_{n+k+1+M_n}^{<\mathrm{aposteriori}>}>1.1{\hat{T}}_{\mathrm{tol}},$ 则拒绝本次隐式投影算法外部积分器的积分计算，设置仿真时间t＝t-H，返回步骤7)，否则进入下一步骤；

13)判断仿真时间t是否达到仿真终止时间T，若t＝T，则仿真结束，否则返回步骤2)，依据步骤2)至13)反复进行直至仿真结束。

2.根据权利要求1所述的一种适于隐式投影算法的误差估计及参数自适应调节方法，其特征在于，步骤1)所述的隐式投影算法内部积分器采用基于显式改进欧拉法的微分代数方程交替求解方法，积分步数k要求大于等于3，隐式投影算法外部积分器采用基于隐式预测-校正方法的微分代数方程交替求解方法进行计算，隐式投影算法外部积分器步长倍数的初始值M_ini值为1。

3.根据权利要求1所述的一种适于隐式投影算法的误差估计及参数自适应调节方法，其特征在于，步骤6)所述的事前误差估计方法是：

设x_n为准确值，经隐式投影算法内部积分器积分k+1步及隐式投影算法外部积分器以步长M_n-1h积分一步后的局部截断误差估计值，通过下式得到

${\hat{T}}_{n+k+2+M_{n-1}}^{<\mathrm{apriori}>}\left|\right|[\frac{1}{6}k(k+1)(k+2)+\frac{{(k+1)}^2}{2}{M}_{n-1}+\frac{k+1}{2}{M}_{n-1}^2+\frac{1}{4}{M}_{n-1}^3]h^3{x}_n^{\left(3\right)}\left|\right|$

式中，状态变量的三阶微分项通过拉格朗日插值方法得到，||x||为欧几里得距离，分别如下式所示

$x_n^{\left(3\right)}\&ap;-\frac{1}{h^3}(x_n-{3x}_{n+1}+{3x}_{n+2}-x_{n+3})$

$\left|\right|x\left|\right|={\left(x^{T}x\right)}^{\frac{1}{2}}$

当对第一步隐式投影算法外部积分器的事前误差进行估计时，上一步的隐式投影算法外部积分器步长倍数M_n-1，即M₀取M_ini。

4.根据权利要求1所述的一种适于隐式投影算法的误差估计及参数自适应调节方法，其特征在于，步骤7)所述的当前隐式投影算法外部积分器步长倍数M_n利用估计误差，并由下式估计得到

$M_n={\left|\frac{{\hat{T}}_{\mathrm{tol}}}{\hat{T}}\right|}^{\frac{1}{3}}{M}_{\mathrm{old}}$

式中，为事前估计误差或事后估计误差M_old为事前误差估计时采用的M_n-1或经事后估计误差判断后被拒绝的隐式投影算法外部积分器的外部积分步长倍数M_{n_rejected}。

5.根据权利要求1所述的一种适于隐式投影算法的误差估计及参数自适应调节方法，其特征在于，步骤11)所述的事后误差估计方法是：

隐式投影算法外部积分器采用隐式预测-校正方法对微分方程组进行积分求解，经过隐式投影算法外部积分器的预测步计算得到的系统状态变量记为经过隐式投影算法外部积分器的校正步计算得到的状态变量记为则隐式投影算法外部积分器的事后估计误差由下式得到

${\hat{T}}_{n+k+1+M_n}^{<\mathrm{aposteriori}>}\&ap;\left|\right|x_{n+k+1+M_n}^{<1>}-x_{n+k+1+M_n}^{<0>}\left|\right|.$