暂态算法

暂态算法（精选7篇）

暂态算法 篇1

0 引言

作为电力系统仿真的重要组成部分,电磁暂态仿真具有现象刻画准确、应用广泛、数值稳定性好等特点,并与机电暂态仿真共同构成了电力系统暂态仿真的基础,其应用涵盖了电力系统规划、设计、运行及科学研究的各个方面,是了解电力系统暂态复杂行为的必要工具。与机电暂态仿真不同,电磁暂态仿真在精确的电路层面上对系统元件进行建模、分析,并计算得到各种暂态响应的时域波形。这使得电磁暂态仿真从模型、算法到计算结果都有别于机电暂态仿真。电磁暂态仿真最初用于电力系统过电压计算、绝缘配合、次同步谐振、谐波分析、保护及控制装置建模、FACTS与HVDC等方面的研究,其基本理论与方法由Dommel于20世纪60年代末提出[1,2]。近年来,电磁暂态仿真也被广泛用于包括大型风力发电和分布式发电在内的各种新型电能生产方式的研究中[3,4]。

针对不同类型的应用,电磁暂态仿真可分为离线仿真工具和实时仿真器。离线仿真工具包括各种常见的软件包,如EMTP-RV,ATP,EMTDC,MicroTran等,它们可安装在普通PC机或工作站上,面向对时间没有严格要求的各种情况下的仿真计算,尽管这些软件都采用了高效的数值算法,但通常来说,仿真计算时间要远多于所研究暂态现象的持续时间。对于实时仿真器,除了软件技术外,还需要相关硬件装置的配合,以保证仿真时刻与外部时钟的精确同步,由此可以为各种电力系统保护与控制装置提供高度模拟现场实际的测试环境。考虑到经济性和硬件条件的限制,目前实时仿真器并不能完全取代离线的仿真工具。作为实时仿真器的代表,实时数字仿真器(RTDS)被广泛应用于工业和学术界。在国内,中国电力科学研究院和殷图公司也都开发了具有电磁暂态仿真功能的实时仿真器,并在实际系统中得到了应用。

近年来,一些新技术、新方法的出现极大地提高了电磁暂态仿真的精度和计算速度,扩展了电磁暂态仿真的应用领域。本文将首先回顾电磁暂态仿真的理论基础,进而分析电磁暂态仿真面临的问题与挑战,结合这些问题介绍了为提高仿真精度和速度所取得的研究进展,最后对今后的研究工作进行了展望。应该说明的是,精确模拟各种暂态过程离不开准确的数据、恰当的模型以及高效的算法,本文将主要介绍电磁暂态仿真算法。

1 问题与挑战

1.1 理论基础

电力系统电磁暂态仿真本质上可归结为对动力学系统时域响应的求取,包括系统本身的数学模型和与之相适应的数值算法。对电力系统而言,其数学模型包括2类:一类是由系统的网络拓扑结构决定的约束方程,即KCL和KVL方程;另一类则是由系统中各元件自身特性决定的伏安关系方程。其中,第1类约束方程是代数方程,第2类方程则可能是代数方程、微分方程或非线性方程。以图1所示的电感为例,其基本伏安关系方程为式(1)给出的微分方程,在正弦交流稳态电路中式(1)退化为如式(2)所示的相量形式的代数方程:

$\begin{array}{l}v{}_k-v{}_m=L\frac{\text{d}i{}_{km}}{\text{d}t}(1)\\ \dot{V}{}_k-\dot{V}{}_m=\text{j}\omega L\dot{{\rm I}}{}_{km}(2)\end{array}$

式(1)为电磁暂态仿真所采用,而式(2)为机电暂态仿真所采用。

当元件的特性方程具有式(2)的代数形式约束时,电网模型可以用节点方程表示为:

$\mathit{Y}\mathit{U}=\mathit{{\rm I}}(3)$

此外,系统中还存在描述发电机及励磁、调速系统动态特性的微分方程。这样,整个电力系统的数学模型可表示为一组代数—微分(DAE)方程组:

$\{\begin{array}{l}\dot{\mathit{x}}=\mathit{f}(\mathit{x},\mathit{y})\\ 0=\mathit{g}(\mathit{x},\mathit{y})\end{array}(4)$

式(4)是机电暂态仿真的基础,在采用了具体的数值积分方法后,它可以联立求解整个差分后的DAE方程组,或采用交替求解的算法分别求解代数方程组和微分方程组[5]。

与式(2)不同,式(1)并不能直接通过联立节点方程形成式(3)形式的网络方程。一种方法是先采用数值积分方法对式(1)进行差分化,得到代数形式的差分方程。以式(1)为例,应用梯形积分法后得到:

$i{}_{km}(t)=\frac{\text{Δ}t}{2L}(v{}_k(t)-v{}_m(t))+{\rm I}{}_{\text{Η}\text{i}\text{s}\text{t}}(t-\text{Δ}t)(5)$

式中:${\rm I}{}_{\text{Η}\text{i}\text{s}\text{t}}(t-\text{Δ}t)=i{}_{km}(t-\text{Δ}t)+\frac{\text{Δ}t}{2L}(v{}_k(t-\text{Δ}t)-v{}_m(t-\text{Δ}t))$。

式(5)可以认为是一个值为Δt/(2L)的电导与历史项电流源并联的诺顿等效电路形式,如图1所示。由此再通过节点方程联立式(5)的差分方程,形成电磁暂态仿真的基本方程:

$\mathit{G}\mathit{u}=\mathit{i}(6)$

此时式(6)的节点方程中已经包含了具体的数值积分方法,从而将系统的数学模型与数值计算方法融合在一起,这也是式(6)与式(4)的主要区别。

另一种更一般的方法则是形成标准形式的状态方程:

$\{\begin{array}{l}\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}\mathit{u}\\ \mathit{y}=\mathit{C}\mathit{x}+\mathit{D}\mathit{u}\end{array}(7)$

与式(4)一样,用状态方程描述的系统模型同样是独立于数值算法的。具有标准形式的状态方程可使用各种成熟的数值计算程序进行求解,唯一的问题是通常并不能直接得到它。与节点方程相比,状态方程的形成要复杂和困难得多,文献[6]提出了一种由节点关联矩阵和支路数据自动生成状态方程的方法,它便于程序实现,但在速度和效率上没有优势。

式(6)和式(7)构成了电磁暂态仿真的2类基本方法,即状态方程法和节点分析法。以状态方程描述的电力系统模型可采用各种显式或隐式的数值积分方法求解,但无论哪种方法,保证算法的数值稳定性都是至关重要的问题。变步长算法通过误差控制调整仿真步长,使显式方法的数值稳定性也能够得到保证。而对于节点分析法,由于式(6)中已经考虑了具体的数值积分方法,因此,对积分方法本身也提出了要求:首先,算法应该是形式简单的,这样才能使历史项电流源的表达式不至于十分复杂;其次,算法应具有良好的数值稳定性;最后,算法应具有较高的精度。对电力系统这样的中度刚性系统而言,单步具有二阶精度且A稳定的梯形法无疑是最佳选择。考虑到梯形法在稳定性方面的优势,可不必再通过变步长算法进行误差控制,同时,采用定步长可以使式(6)中的G矩阵保持恒定,对于一些元件模型,如分布参数线路中的历史项电流源的计算也更为简便。这样,基于定步长的梯形法得到的节点方程(6)构成了电磁暂态EMTP类程序的基础。状态方程法虽然比较灵活,但节点分析法在计算速度上有更大优势,特别是对于大规模电力系统而言,状态方程法的仿真时间往往是无法接受的,因此,很多电力系统电磁暂态仿真的专门软件都采用了节点分析法的基本框架,它们统称为EMTP类程序,本文也将着眼于EMTP类程序的相关算法。

1.2 存在的问题

电磁暂态仿真理论与方法发展至今已有半个世纪的时间,限于当时的软硬件水平,为了保证计算速度和效率,算法中对一些问题做了必要的简化和处理。随着现代电力技术的不断发展,过去曾采用的一些方法已不能满足当前的要求,电磁暂态仿真面临新的问题与挑战,这包括:电力系统中不断涌现的各种新设备对电磁暂态仿真工具的建模能力提出了更高的要求,在这些元件模型中,有些已经比较成熟,可以纳入到程序的内置模型中,如风力发电系统,但更多的还需要不断发展和完善,提供强大的用户自定义建模功能就显得尤为必要;现代电力技术的发展使得电力电子装置与控制系统越来越多地出现在各种应用场合,精确计及开关动作与控制系统的动态过程对仿真结果起着至关重要的影响;近来,电磁暂态仿真也被越来越多地用于低压配电网的研究[7],随着计算规模的增加,如何在不牺牲计算精度的前提下提高计算速度是一个问题;电磁暂态仿真起始于网络的稳态运行点,但获得需要的稳态解通常并不容易[8],特别是当系统含有大量的电力电子和控制设备时,整个系统的初始化过程就更为复杂和困难,一种折中的办法是采用零状态初始化,这本身就相当于一次暂态过程的计算。文献[9]讨论了EMTP类程序初始化的3种方法,并认为基于潮流的初始化方法更为有效。文献[10]介绍了一种提高控制系统初始化效率的方法,但更为一般的自动初始化方法尚未实现。

2 研究进展

2.1 扩展的网络方程

由于节点方程的一些限制,Ho等人在20世纪70年代提出了改进节点方程[11],能很好地处理理想电压源和各种受控源模型,扩展了网络的建模能力,并被SPICE类程序所采用[12]。考虑到电力系统实际以及式(6)的对称性,直到最近一种增广形式的改进节点方程才被应用到EMTP-RV中[13,14],方程如下:

$\begin{array}{l}\text{节}\text{点}\text{方}\text{程}\\ \text{电}\text{压}\text{源}\text{支}\text{路}\\ \text{受}\text{控}\text{源}\text{支}\text{路}\\ \text{开}\text{关}\text{支}\text{路}\end{array}\left[\begin{array}{cccc}\mathit{Y}{}_{\text{n}}& \mathit{V}{}_{\text{c}}& \mathit{D}{}_{\text{c}}& \mathit{S}{}_{\text{c}}\\ \mathit{V}{}_{\text{r}}& \mathit{V}{}_{\text{d}}& \mathit{D}{}_{\text{V}\text{D}}& \mathit{S}{}_{\text{V}\text{S}}\\ \mathit{D}{}_{\text{r}}& \mathit{D}{}_{\text{D}\text{V}}& \mathit{D}{}_{\text{d}}& \mathit{S}{}_{\text{D}\text{S}}\\ \mathit{S}{}_{\text{r}}& \mathit{S}{}_{\text{S}\text{V}}& \mathit{S}{}_{\text{S}\text{D}}& \mathit{S}{}_{\text{d}}\end{array}\right]\left[\begin{array}{l}\mathit{v}{}_{\text{n}}\\ \mathit{i}{}_{\text{V}}\\ \mathit{i}{}_{\text{D}}\\ \mathit{i}{}_{\text{S}}\end{array}\right]=\left[\begin{array}{l}\mathit{i}{}_{\text{n}}\\ \mathit{v}{}_{\text{b}}\\ \mathit{d}{}_{\text{b}}\\ \mathit{s}{}_{\text{b}}\end{array}\right](8)$

式(8)消除了节点方程在处理不接地理想电压源和理想变压器等模型时的困难,同时可以方便地表示理想开关模型并保证矩阵维数恒定。式(8)可用于稳态和暂态计算,采用直接法即可求解,但与式(6)相比,式(8)已不再是对称矩阵。

2.2 改进计算精度的方法

同早期的电力系统相比,现代电力系统使用了越来越多的电力变流器,这些电力电子装置在进行电磁暂态仿真时表现为大量的开关模型。由于EMTP类程序采用定步长梯形法,算法只能在步长的整数倍时刻改变开关状态,除了梯形法引起的数值振荡外,开关动作时间上的延迟还会导致许多其他问题[15]。因此,为了保证仿真的准确性,必须要求定步长算法能够精确考虑开关的动作时刻。一种方法是在出现开关动作的步长内改用更小的步长积分到开关动作时刻,缺点是需要重新计算各元件的等效电导并增加了一次矩阵求解。更为有效且被广泛使用的方法是采用图2所示的线性插值,此时不需要重新积分就能“还原”到开关动作时刻状态改变以前系统中各变量的值。线性插值算法简单、方便、有效,它假设在相邻的2次开关动作之间的系统特性可以用线性关系拟合,这样的假设在小步长的情况下是合适的。插值算法可以消除非特征谐波,但并不能解决梯形法产生的数值振荡问题。文献[16]通过将插值与临界阻尼调整(CDA)法[17]结合,提出了用于电力电子变换器仿真的“反向追踪”技术。尽管如此,当存在同步开关时,必须在同步开关动作时刻对系统重新进行初始化才能得到正确的仿真结果[18]。文献[18,19,20]分别提出了不同的“重初始化”方法,以解决自然换向和强制换向问题。

控制系统建模与仿真是电磁暂态仿真的重要组成部分,最初是通过子程序TACS(transient analysis of control systems)的方式提供给EMTP的[2]。TACS的开发最早是为了模拟HVDC换流站控制系统,但它随后被应用到包括发电机励磁调速系统建模在内的诸多方面。EMTP中控制系统与电气网络分开求解,使得它们之间存在Δt的时延,如图3所示,这个时延称为外部时延。一些算例表明,外部时延可能会导致数值失稳或不正确的结果。文献[21,22]分别提出了不同的方法,以解决外部时延产生的数值问题。除了外部时延,控制系统内部非线性环节还会产生内部时延,这些时延会影响控制系统内部求解时的数值稳定性。文献[23]提出用牛顿法消除内部时延以获得控制系统的联立解,该方法已被EMTP-RV所采用。文献[24]提出了对电气和控制系统采用联立求解的思想,对此,文献[25]采用一种基于电路的方法加以实现。

包括同步电机和异步电机在内的电机模型是电磁暂态仿真中最为复杂的模型之一,其复杂性不仅体现为高维数和非线性,也包括采用不同于常规元件的坐标系。经过派克变换的电机模型具有恒定的电感系数矩阵,这是它得到广泛应用的主要原因,但却存在电机模型与网络其他部分的接口问题。以异步电机为例,ATP是采用补偿法求解的,此时多个电机之间必须通过输电线模型分隔开,以保证每个电机单独求解[2];而在EMTDC中,异步电机表示为诺顿电流源的形式,它是由上一步的机端电压计算得到的[26],因此存在Δt的时延;MicroTran不仅预报机械量还要预报电气量(磁链)。尽管有文献认为dq轴表示的电机模型在不对称运行条件下会产生不正确的结果[27,28],但文献[29]经过比较后认为采用小步长仿真不对称算例时不会产生问题,但当步长较大时则可能会遇到收敛问题。在相域对电机进行建模的研究最早可以追溯到文献[30]的工作,相域电机模型在数值稳定性方面的优点使其在实时仿真等领域又重新引起人们的重视[31,32],此外,相域模型还更适于电机的饱和特性以及电机内部故障的研究。文献[33]提出的VBR(voltage-behind-reactance)同步电机模型保留定子侧abc相坐标形式,而转子侧采用dq坐标,这样使得电机可以非常方便地与电网接口,文献[34]将这一思想扩展到异步电机模型,算例表明VBR电机模型具有很好的数值稳定性。

2.3 提高计算速度的方法

除了采用高效的稀疏技术外,对网络进行分块是另一种提高网络方程求解效率的思想。由此,基于改进节点方程,文献[35]提出了多区戴维南等值(MATE——multi-area Thevenin equivalent)的方法用于网络方程的计算。其基本思想是:将整个电力系统网络通过一些支路分割成多个子网络,首先分别计算各子网络不含连接支路时的节点电压,然后得到连接支路上的电流,最后叠加连接支路电流的响应得到最终解,整个计算过程与采用戴维南形式的补偿法是一致的。MATE方法的提出为大规模网络的求解提供了极大的便利,一种显而易见的想法是采用并行策略提高计算速度,文献[36]通过节点分裂实现网络分割,进而提出了一种电磁暂态仿真并行算法。除了实现并行计算,MATE还可用于含有大量开关元件的网络求解,通过将网络中含有开关的部分与其他部分分割形成子网络,开关状态的改变仅仅影响其所在的子网络,而不用对其他不含开关的子网络重新进行因子分解,这与文献[26]介绍的通过节点编号实现网络矩阵的局部重新因子分解的思想很相近。此外,MATE还为在不同子系统分别采用不同算法、不同步长进行求解提供了可能。文献[37]在MATE的基础上将子网络继续进行分块提出了多层MATE的方法并用于控制系统和非线性元件的仿真。

传统的电力系统仿真通常采用统一的算法和单一的步长对整个系统进行求解。当系统中同时存在快动态过程和慢动态过程,并且快动态过程迅速衰减时,可采用变步长算法,但是当快动态扰动持续不断时,变步长算法则是低效的。为此,一种思想是将快动态过程和慢动态过程分别进行求解,这种多速率的概念最早见于文献[38],并在近年来得到重视,文献[39]最早将多速率的概念用于电力系统动态仿真。基于多速率仿真的思想,文献[40,41]提出了Latency技术用于电磁暂态层面的仿真。Latency技术基于MATE通过在不同子区域使用不同的仿真步长实现快动态过程与慢动态过程的分区联立求解,在慢动态系统的一个步长内,它通过内插对快动态系统等值,算例表明Latency技术能显著地减少程序的计算量并获得较高精度,文献[42]对Latency技术的精度和数值稳定性进行了详尽的分析。

为了详细考虑FACTS与HVDC对电力系统暂态稳定性的影响,混合仿真[43]技术对这些设备进行了详细的电磁暂态建模,而网络的其余部分仍采用机电暂态模型。由于电磁暂态与机电暂态仿真分别采用不同的坐标系与仿真步长,接口处理是混合仿真的关键问题,文献[44,45]介绍了这方面的工作。电磁暂态与机电暂态仿真在模型、算法之间的差异导致了混合仿真接口的复杂性,同时还会产生频率偏移、基频分量抽取误差等问题[46]。基于动态相量的元件模型简化了混合仿真接口的复杂性,与传统相量相比,动态相量是建立在时变傅里叶分解的基础上,突破了准稳态假定的限制,具有更大的频率带宽,可用于包括电力电子设备在内的电力系统建模[47,48,49]。动态相量法能有目的地选择占主导优势的频率进行相域内的分析,既能够分析电磁暂态过程又能够分析机电暂态过程,可以提供比传统相量法更精确的模型,同时获得比详细时域模型更快的仿真速度。

提高计算速度的各种算法通常还需要硬件装置的配合才能最大限度地发挥其优势,如MATE为仿真的空间并行提供了基础[50],而Latency技术则为算法的时间并行提供了可能,混合仿真可以利用现有的各种机电暂态并行仿真策略。单纯的并行计算可采用各种高性能机群计算机实现,若实现实时仿真则需要各种软硬件平台的配合,特别是要考虑与物理装置的接口问题,目前包括数字信号处理器(DSP)、精简指令集计算机(RISC)、复杂指令集计算机(CISC)和大规模集成电路(VLSI)在内的各种硬件技术在实时仿真中均有应用[51,52],但基于PC机群的实时仿真器可能会由于性能和价格上的优势成为未来的发展方向[53]。

2.4 其他方法

电磁暂态仿真可以得到各种暂态响应的时域波形,而应用动态相量法可以得到这些时域信号的包络线,它们之间通过傅里叶变换联系起来。文献[54]提出了一种基于信号处理的支路友模形式,通过改变移频参数可以分别得到响应信号的瞬时波形及其包络线。与动态相量法一样,其支路友模具有复数形式等效电导和历史项电流源,除了仿真步长,移频参数也被包含在支路友模的计算中,这样的思想也体现在文献[55]中,通过调整阻尼法中的权重系数实现在整个网络的不同部分应用不同数值算法的目的。在此基础上,文献[56]在支路友模中考虑了步长的变化,在保证等效电导恒定的情况下,以算法的精度换取仿真步长在一定范围内的变化,这对于插值算法的再同步过程十分有效。

在频域中,文献[57]总结回顾了基于z变换的方法在电磁暂态仿真中的应用,主要是用于输电线路的建模,对其他元件的建模还需要进一步研究;此外,文献[58,59]通过对频率响应的近似提出了不同的网络等值和化简方法。根匹配法属于传递函数数字仿真的一类方法,通过在z平面上一一对应地确定s域传递函数的零、极点位置,得到输入、输出关系的差分公式,只要原系统是稳定的,不论步长取多大,都能保证模型也是稳定的。文献[60]将根匹配法用于电磁暂态仿真得到具有指数形式的支路友模,该方法具有精度高、数值稳定性好等优点,但由于不是所有元件模型都可以得到s域的传递函数,因此,根匹配法无法取代各种数值积分方法。

3 结语

电力系统电磁暂态仿真理论与方法从建立至今已得到了极大发展,各种新技术、新方法的出现改进了原有算法的不足,提高了仿真精度和计算速度,使电磁暂态仿真有了更强的生命力和更加广泛的应用前景。本文首先回顾了电磁暂态仿真的理论基础,它决定了电磁暂态仿真采用什么样的模型和数值算法,随后从几个方面介绍了电磁暂态仿真算法研究取得的进展,这些工作有的已经比较成熟,开始得到应用;有的还停留在理论研究范畴,需要进一步的发展和完善。

摘要：回顾了基于节点分析的电磁暂态仿真算法的理论基础,在此基础上分析了算法在应用中遇到的问题和挑战,并分别从模型的扩展、计算精度的改进、速度的提高等几个方面介绍了电磁暂态仿真算法所取得的研究进展。

关键词：电磁暂态,仿真,节点分析,插值,梯形法

暂态算法 篇2

关键词：中性点非直接接地,相频特性,暂态零序电流,FTU

0引言

我国配电网多采用中性点非直接接地方式, 单相接地故障快速定位一直是供电企业的难点之一。本文从理论上对中性点非直接接地电网的零序网络进行相频特性分析, 根据暂态零序电流的特征, 提出一种利用相邻FTU检测的暂态零序电流做出快速故障定位的方法。这种方法能够准确定位出中性点非直接接地电网中故障区域, 减少了故障时过渡电阻对故障定位造成的影响, 具有较高的灵敏度。

1 零序网络的相频特性

1.1 单条无故障线路的零序阻抗相频特性

设单条10kV架空线路每千米的零序电阻为R0, 零序电感为L0, 零序电容为C0, 其线路长度为M, 角频率为ω, 则根据典型分布参数模型, 这条线路的输入阻抗Z0如下:

根据式 (1) , 该条10kV架空线路的零序阻抗相频特性如图1所示。

分析图1可以发现, 单条10kV架空线路零序阻抗的相频特性是从+90°~-90°呈周期变换的曲线。随着频率加大, 该架空线路的零序阻抗在容性频带以及感性频带之间相互转换。第一次频带转换的频率在2kHz左右。

1.2 中性点不接地电网的零序阻抗相频特性

假设在多条10kV架空线路构成的系统中, 有一条线路发生了单相接地故障。对于没有发生故障的线路, 其零序阻抗的相频特性只和本条线路相关, 即与图1所示的相频特性相同。而对于故障线路, 其检测到的阻抗与其余所有线路构成的网络有关, 由于该网络内的线路同时存在自身的串联谐振以及相互之间的并联谐振, 零序阻抗具有不规则的感性或容性。即故障线路的相频特性与其他正常线路相关, 故障线路发生串联谐振时, 其频率等于其余线路发生串联谐振的最小值。设有4条线路, 1~3条为正常线路, 第4条线路为故障线路, 则这4条线路的相频特性如图2所示。

由图2可见, 在ω′之前, 所有线路的阻抗都为容性。

2 暂态零序电流的特点

在图2中, 0~ω′频段内各线路的阻抗都为容性, 用电容等效, 这个频段内暂态零序电流为容性电流。对于故障线路, 其出口点监测到的零序电流, 为其余正常线路对地零序电容电流之和, 电流方向为由线路流向母线。在线路发生接地故障时, 在故障点相当于出现了一个虚拟电源, 从该点流出的零序电流向线路上游和下游流去, 故障区段两端检测到的暂态零序电流方向相反, 在幅值上上游电流更大。

因此0~ω′频段内的暂态零序电流具有以下特点: (1) 故障线路的暂态零序电流幅值大于或等于正常线路; (2) 故障线路的暂态零序电流由线路流向母线, 而正常线路的暂态零序电流由母线流向线路; (3) 故障线路的正常段两端电流幅值基本相等, 且电流流向相同, 而故障线路的故障区段两端电流方向相反, 且故障区段流向上游的电流幅值会远大于流向下游的电流幅值。

3 基于暂态零序电流幅值的故障快速定位方法

3.1 建立馈线自动化系统

建立典型的环网馈线自动化系统, 配置FTU、通信网络及FA控制主站。其工作方式为:在线路发生接地故障后, FTU能够检测到过电流, 利用通信网络将故障信息发送给控制主站后, 控制主站对该信息进行分析, 找出故障区段。根据馈线自动化系统给出的结论, 操作分段开关将故障区段隔离, 恢复其余正常线路的供电。

3.2 基于暂态零序电流幅值的故障快速定位算法

设第k条线路在0~ω′频段内暂态零序电流分量为i0k, 其第i个数据有i0ki, 则暂态零序电流的有效值I0k计算公式如下:

式中, n为数据的总个数。

设线路上游的暂态零序电流为I01, 线路下游的暂态零序电流为I02, 则上下游暂态零序电流的比值η如下:

从上文可知, 对于线路正常区段, η约等于1;而对于故障区段, η远大于1。所以, 利用线路各区段两端的FTU检测到的上下游暂态零序电流的比值η, 即可判断故障区段。

当η<ζ (ζ可整定) 时, 则判断该区段两端暂态零序电流幅值相差不大, 为正常区段;当η>ζ (ζ可整定) 时, 则判断该区段两端暂态零序电流幅值相差过大, 为故障区段。

3.3 算法特点

采用上述算法, 由于FTU发送的只是暂态零序电流幅值, 而不是全部暂态零序电流数据, 故传输数据量小, 对通讯网络的压力较小。同时, 这种算法在采集暂态零序电流幅值时不需要时间同步, 即时间不同步也不影响对故障区段的判断。从相频分析可以看出, 暂态零序电流幅值远大于工频分量, 保证了故障发生后监测结果的准确性。

4结语

线路的零序阻抗相频特性呈周期性变化, 第一次线路零序阻抗从容性频带转换为感性频带时的频率远大于工频, 在0~ω′频段内暂态零序电流为容性电流。而线路正常区段和故障区段两端的电流具有明显的不同, 采用基于暂态零序电流幅值的故障快速定位算法, 其传输的数据较少, 对故障的定位具有较高的灵敏性和可靠性。

参考文献

[1]杜刚, 刘迅, 苏高峰.基于FTU和“S”信号注入法的配电网接地故障定位技术的研究[J].电力系统保护与控制, 2010, 38 (12) :73-76.

[2]张彩友.单相接地故障指示器技术现状分析[J].电网技术, 2007, 31 (S2) :280-283.

暂态算法 篇3

电力市场环境下,仅考虑静态安全约束的最优潮流(optimal power flow,OPF)使得系统运行点更容易接近稳定极限,已无法满足现代电力系统的要求,考虑暂态稳定约束的OPF(transient stability constrained OPF,TSOPF)成为新的研究方向。

OPF是一个大规模的非线性优化问题,内点法及其改进算法由于良好的计算效率和收敛特性在该问题求解中得到了广泛应用[1,2,3,4]。考虑暂态稳定约束后,增加了大量的微分代数方程约束,是一类特殊的无穷维优化问题[5]。目前关于TSOPF的研究集中于暂态稳定约束条件的处理。文献[6]将发电机暂态微分代数方程差分后作为等式约束加入到OPF;文献[7]则考虑了多个预想故障,并引入既约导纳矩阵以消去网络方程;文献[8]采用大步长对预想故障集进行故障预选和同调识别,对于小步长精确计算仅考虑“起作用”故障;文献[9,10]采用约束转换技术,通过欧氏空间变换将函数空间的TSOPF转换为相同规模的静态优化问题;文献[11]将该问题分解为最优潮流和最优控制2个子问题,通过交替求解这2个子问题得到原问题的解;文献[12]将暂态能量裕度作为附加约束引入到OPF模型中;文献[13]利用暂态稳定裕度指标及其灵敏度组成的不等式代替暂态稳定约束。逐次线性规划[6]、内点法[7]以及智能优化算法[14]均被成功应用于该问题的求解。

常规方式下,基于时域仿真差分思想的TSOPF将发电机转子运动方程离散化,并作为OPF问题的等式约束。由于任何差分方法都具有截断误差,将其作为等式约束过于苛刻。内点法求解过程中,这种处理方式使得修正方程组的阶数随着系统规模和故障个数急剧上升,导致计算耗时长、内存消耗大甚至不可解的“维数灾”。针对该问题,本文将描述系统暂态过程的差分方程转换为不等式约束,其上下限与所采用的差分方法精度相关,从而大大降低了内点法计算过程中修正方程组的阶数。对3个不同规模的系统进行测试,结果显示与常规方式相比,该方法计算时间明显缩短,占用内存减少,能够求解更大规模系统的TSOPF问题。

1 TSOPF问题描述

考虑多个预想故障的TSOPF问题可表示为:

$\{\begin{array}{l}\underset{\mathit{x}}{{\displaystyle \min }}f(\mathit{x})\\ \text{s}.\text{t}.\mathit{h}(\mathit{x})=0\\ \underset{¯}{\mathit{g}}\le \mathit{g}(\mathit{x})\le \overline{\mathit{g}}\end{array}(1)$

式中:f(x)为该问题的目标函数,为不失一般性,本文采用发电费用最小[7]为目标函数;h(x)为等式约束向量,包括潮流方程、发电机初值约束和暂态微分代数方程离散化后得到的差分方程;g(x)为不等式约束向量,包括发电机有功和无功出力、变压器变比、节点电压和线路有功潮流上下限等静态安全约束以及暂态稳定约束;$\overline{\mathit{g}}$和$\underset{¯}{\mathit{g}}$分别为不等式约束上下限。详细的问题模型可参阅文献[7,8]。

TSOPF问题的等式约束包含大量暂态差分方程。常规方式下,若发电机采用经典模型,负荷用恒定阻抗表示,由隐式梯形法可得如下差分方程[8]:

$\{\begin{array}{l}\text{Δ}\delta {}_i^t(k)=\delta {}_i^t(k)-\delta {}_i^{t-1}(k)-\frac{\omega {}_0\text{Δ}t}{2}(\omega {}_i^t(k)+\\ \omega {}_i^{t-1}(k))=0\\ \text{Δ}\omega {}_i^t(k)=\left(\frac{D{}_i\text{Δ}t}{2{\rm M}{}_i}+1\right)\omega {}_i^t(k)+\left(\frac{D{}_i\text{Δ}t}{2{\rm M}{}_i}-1\right)\cdot \\ \omega {}_i^{t-1}(k)-\frac{\text{Δ}t}{2{\rm M}{}_i}(2{\rm P}{}_{\text{G}i}-{\rm P}{}_{\text{e}i}^t(k)-\\ {\rm P}{}_{\text{e}i}^{t-1}(k))=0\end{array}(2)$

式中:i∈nG,k∈nC,t∈nt,nG,nC,nt分别为发电机、预想故障和积分时段数目;ω0为同步角速度有名值;δti(k),ωti(k)分别为发电机i第k个故障第t时段的功角和角速度相对同步转速的偏差;Δt为积分步长;Di,Mi,PGi分别为发电机i的阻尼系数、惯性常数和有功出力;Ptei(k)为发电机i的电磁功率,可采用既约导纳矩阵表示以消去网络方程[7]。

2 降阶内点算法

基于内点法的TSOPF计算其主要任务是求解对称非正定线性方程组:

$\left[\begin{array}{cc}{\mathit{{\rm H}}}^{\prime }& \nabla {}_{\mathit{x}}\mathit{h}(\mathit{x})\\ \nabla {}_{\mathit{x}}^{\text{Τ}}\mathit{h}(\mathit{x})& 0\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathit{x}\\ \text{Δ}\mathit{y}\end{array}\right]=\left[\begin{array}{c}\mathit{L}{}_{{\mathit{x}}^{\prime }}\\ -\mathit{L}{}_{\mathit{y}}\end{array}\right](3)$

式中:x为控制和状态变量;y为等式约束乘子;H′为障碍矩阵和目标函数、等式及不等式约束对应海森矩阵的线性组合[7];Lx′和Ly为库恩—塔克(KKT)系统的右端残差向量,其定义可参看文献[15]。

式(3)的求解占用了绝大部分计算时间[7],其阶数将直接影响到内点法的计算耗时。由式(3)可知,系数矩阵的阶数取决于等式约束和变量的数目,与不等式约束个数无关。如果设法减少等式约束个数,则可降低式(3)的阶数。常规方式下,差分TSOPF等式约束中绝大部分为暂态差分方程。事实上,由于任何数值算法都存在误差,将这些差分方程作为严格等式约束并无必要。对于初值问题:

${\mathit{v}}^{\prime }=\text{ϕ}(t,\mathit{v})t\ge t{}_0,\mathit{v}(t{}_0)=\mathit{v}{}_0(4)$

采用隐式梯形法差分化,第j步到第j+1步产生的截断误差为[16]:

v(tj+1)-$\mathit{v}(t{}_j)+\frac{\text{Δ}t}{2}(\text{ϕ}(t{}_j,\mathit{v}(t{}_j))+\text{ϕ}(t{}_{j+1},$

v(tj+1)))$=-\frac{(\text{Δ}t){}^3}{12}{\mathit{v}}^{‴}(t{}_j)+{\rm O}((\text{Δ}t){}^4)(5)$

即隐式梯形法是二阶精度的差分方法,存在3阶局部截断误差。实际上,任何一种差分方法都存在误差,将差分方程作为等式约束是非常苛刻的条件,直接导致式(3)阶数增加,带来计算时间和内存消耗过大等问题。在大规模系统、多预想故障的TSOPF求解中,这些问题更加突出。若将这些差分方程转换为不等式,则可以大大减少等式约束个数,即式(2)转换为:

$\{\begin{array}{l}-\beta (\text{Δ}t){}^n\le \text{Δ}\delta {}_i^t(k)\le \beta (\text{Δ}t){}^n\\ -\beta (\text{Δ}t){}^n\le \text{Δ}\omega {}_i^t(k)\le \beta (\text{Δ}t){}^n\end{array}(6)$

式中:Δt为积分步长;n为所用差分方法的精度阶数。

差分不等式约束的上下限设置为±β(Δt)n,其中β为约束系数,β∈[0,1]。在内点法求解过程中,将1个差分方程由等式约束转换为上下限不等式约束,可减少1个等式约束及其乘子(y),同时增加2个松弛变量(l,u)以及2个不等式约束乘子(z,w)。由于松弛变量和乘子通过直接回代求解,计算耗时很少,故这些变量的增加对优化过程的总耗时影响甚小,测试结果也将证明这一点。

如图1所示,假设t0时刻系统中某处发生故障,tcr时刻故障被切除,tmax时刻暂态仿真过程结束。以tb表示暂态过程中等式约束和不等式约束的分隔时刻,即t0～tb时段的暂态差分方程视为等式约束,tb～tmax时段的暂态差分方程视为不等式约束。以隐式梯形法为例,说明本文采用的差分方程不等式化降阶过程。

步骤1:发电机采用经典模型,负荷视为恒定阻抗模型,采用隐式梯形法将发电机转子运动方程差分化,得到差分方程组。

步骤2:以tb为界限,将差分方程组分别转换为等式约束集S1和不等式约束集S2:

$\begin{array}{l}S{}_1=\{\begin{array}{l}\text{Δ}\delta {}_i^t(k)=0\\ \text{Δ}\omega {}_i^t(k)=0\end{array}t{}_0\le t\le t{}_{\text{b}}(7)\\ S{}_2=\{\begin{array}{l}-\beta (\text{Δ}t){}^n\le \text{Δ}\delta {}_i^t(k)\le \beta (\text{Δ}t){}^n\\ -\beta (\text{Δ}t){}^n\le \text{Δ}\omega {}_i^t(k)\le \beta (\text{Δ}t){}^n\end{array}t{}_{\text{b}}<t\le t{}_{\max }(8)\end{array}$

式中:tb∈[t0,tmax]。

如考虑到故障时段对整个暂态过程的重要影响,可令tb=tcr。若取tb=t0,即将所有差分方程处理为不等式约束;若取tb=tmax,即常规方式。对于隐式梯形法,精度阶数n为2。

步骤3:将S1和S2分别作为TSOPF问题的等式约束和不等式约束,采用内点法求解降阶后的TSOPF问题。

由于式(5)中v(tj)的计算比较困难,故差分不等式上下限中引入约束系数β。TSOPF问题中,各积分时段的δt, ωt都与前一时刻的δt-1,ωt-1有关,因此β的取值关系到降阶内点算法的精度。

3 算例分析

3.1 测试算例

本文采用预测—校正内点法对3个算例进行了计算。表1给出了3个测试算例的参数,其中nB,nL分别为系统节点和线路数目,t0=0 s。对于所有算例,均采用第2节提出的处理方式,分别取tb=tcr和tb=t0进行计算,并与常规处理方式所得结果进行对比。本文中也就β和tb的取值对计算耗时和目标函数值的影响进行了分析。程序采用各发电机功角相对于惯性中心(center of inertia,COI)的偏差不超过±100°作为暂态稳定判据。为简化计算,假定变压器变比恒定。计算中最大迭代次数设定为50。采用隐式梯形法进行暂态稳定分析时,积分步长可取0.01 s～ 0.02 s,甚至更长[15],本文采用0.02 s的积分步长,式(8)中精度阶数n为2。

采用32位MATLAB R2009a编程求解TSOPF问题,并用MATLAB集成的常微分方程解法器ODE45对程序求得的摇摆曲线进行验证。仿真环境为Intel Core 2 Duo E8400, 4 GB内存。程序已充分考虑矩阵稀疏性,且计算过程中已自动对稀疏矩阵进行优化编号以减少新增非零注入元。

3.2 测试结果

为叙述方便,称tb=tmax的常规方式为TSO1,tb=tcr和tb=t0的降阶处理方式分别为TSO2和TSO3。表2给出了上述3个算例在各处理方式下TSOPF问题的规模和计算结果。

表2中,re为等式约束数,Nd为线性方程组(3)的阶数,Pr为阶数百分比,Ni为内点法迭代求解次数,Tc为求解式(3)的总耗时,Ta为各次迭代求解式(3)的平均耗时,Tb为各次迭代回代求解及更新松弛变量和乘子的平均耗时,Vobj为目标函数值,Dobj′为当前目标函数值Vobj相对于常规方式下目标函数值VTSO1obj偏差的百分比,即

$D{}_{\text{o}\text{b}\text{j}}´=\frac{|V{}_{\text{o}\text{b}\text{j}}-V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}1}|}{V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}1}}\times 100％(9)$

附录A图A1～图A3给出了上述3个算例在各处理方式下所得的摇摆曲线,可知,在给定的预想故障下,时域仿真得到的摇摆曲线与TSOPF计算所得摇摆曲线基本一致,说明该降阶方法所得最优潮流结果满足预先设定的暂态稳定约束条件。表2中数据Dobj′说明降阶后所得目标函数值与常规方式的偏差在可接受的范围内。

对于所有算例,比较表2中各处理方式下的测试结果Nd,Pr,Tc和Ta可知,降阶后各算例中修正方程组的阶数降低至接近原来的一半,迭代求解该线性方程组的总耗时和平均耗时均大大减少,说明该降阶内点算法对所考虑的算例均具有明显的省时效果。其中对于规模较小的系统(IEEE39),由于常规方式的计算耗时本来就不多,该算法的省时效果并不明显;对于较大规模系统(CASE162),该算法较常规方式节省的时间在10倍左右;对于更大规模的系统(CASE300),常规方式下由于内存限制该问题已不可求解,而经过差分不等式降阶处理后也能够在较短时间内给出优化结果。这些都充分说明了该降阶内点算法在计算耗时和内存耗用量方面的优势。

表2的测试结果显示,每次迭代过程中回代求解松弛变量和乘子的时间Tb相对于求解线性方程组(3)的时间Ta非常少。对于所有算例,TSO2和TSO3对应的Tb与TSO1相比增加很少,说明降阶后由于松弛变量和不等式约束乘子增加带来的计算耗时增长有限,从而印证了第2节的结论。

3.3 β取值分析

图2给出了各算例取tb=tcr时,β取不同值对计算耗时和目标函数值的影响。其中,T1为β取不同值时各算例中求解式(3)的总耗时;T2为表2中各算例TSO2方式下的总耗时;Dobj为当前目标函数值Vobj相对于各算例TSO2方式下目标函数值VTSO2obj的偏差百分比,其定义为:

$D{}_{\text{o}\text{b}\text{j}}=\frac{|V{}_{\text{o}\text{b}\text{j}}-V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}2}|}{V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}2}}\times 100％(10)$

从图2(a)可知,β取值小使得总的计算时间增加(图中β=10-5时,程序达到迭代次数上限,各算例均存在收敛性问题)。图2(b)说明β过大将导致暂态稳定约束过于松弛,目标值偏差较大,优化结果不理想。结合图2(a)和图2(b)可知,对于所考虑的算例,采用隐式梯形法差分化后,β取0.001～0.020可在计算耗时和暂态仿真及目标函数值精度之间取得较好的折中。

表3给出了在所建议的β取值范围内,各算例通过ODE45进行暂态仿真所得发电机的最大相对功角δmax和目标值偏差百分比Dobj′(如式(9)所示)。

从表3结果看出,差分不等式降阶处理所得的暂态仿真和目标函数值精度均在可接受的范围内,说明这种将差分方程部分或全部作为不等式约束的方法保证了足够的计算结果精度,且节省了大量计算时间,因此是可行且有效的。

3.4 tb取值分析

图3给出了不同tb时各算例计算总耗时的变化曲线,其中β按表2所示各算例TSO2方式下取恒定值,T1为tb取不同值时各算例中每次迭代求解式(3)的平均耗时,T2为表2中各算例对应TSO2方式下的平均耗时。在所考虑的算例中,各故障的切除时间tcr均为0.2 s。

从图3看出,随着tb的增大,由于修正方程组阶数不断增大,计算耗时总体上是不断增加的,对于大规模系统则更加明显。这说明对于同一个算例,修正方程组的阶数是影响计算耗时的一个重要因素,从而进一步证明了本文所提出的降阶内点算法的合理性。当tb≤tcr时,计算耗时不会明显增加。注意对于算例CASE300,tb>0.2 s后该问题的求解由于内存不足无法完成。

4 结语

本文提出了一种多预想故障TSOPF问题的降阶内点算法。该方法考虑了差分方法的截断误差,将暂态差分方程转换为不等式,其上下限与所用差分方法的精度相关,并采用内点法求解降阶后的TSOPF问题。对3个不同规模系统的测试结果证明该降阶算法能有效减少计算时间,且测试结果显示出以下特点:

1)降阶内点算法中,线性方程组阶数比常规方式明显减小,其阶数接近常规方式的50%。对于不同规模的系统均可节省大量计算时间,甚至能够求解常规方式下不能求解的问题。

2)测试证明,β取值过小可能导致总的计算时间增加;β取值过大则导致优化结果不理想。若采用隐式梯形法差分化,对于所有算例β,取0.001～0.020时能够取得满意的结果。

3)测试证明,tb取值越大,修正方程组阶数越高,其求解耗时总体上增加,对于大规模系统尤其如此。对于本文的算例,tb≤tcr是一个较合理的取值。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

暂态算法 篇4

中低压电网主要采用的是小电流接地系统,有中性点不接地、经消弧线圈接地和经高阻抗接地3种方式,我国6～66 kV电力系统大多数采用的是中性点不接地或经消弧线圈接地运行方式。小电流接地系统的故障主要是单相接地短路故障,发生单相接地后,故障相对地电压降低,非故障两相的相电压升高,但线电压却依然对称,因而不影响对用户的连续供电,系统可运行1～2 h,这也是小电流接地系统的最大优点。但是若发生单相接地故障时电网长期运行,因非故障的两相对地电压升高,可能引起绝缘的薄弱环节被击穿,进一步扩大为相间短路故障,对电力系统造成重大损失。故很有必要对小电流接地系统的单相接地短路故障的定位进行深入研究,以便快速准确地排除故障。

近年来,国内外学者提出了许多行之有效的配电网故障选线、定位方法,如零序电流法、行波测距法、脉冲信号注入法等[1,2,3,4,5],但这些方法也存在不少缺点。在谐振接地系统中,由于消弧线圈的补偿作用,故障线路零序电流的变化不明显,零序电流法效果不佳;行波测距法的故障点反射波头很难确定;脉冲信号注入法需要特定的装置注入信号,成本投入较高。

本文在对小电流接地系统单相故障暂态零模分量特征分析的基础上,提出了一种利用暂态零模电流的内积、有效值算法进行故障选线,利用相邻检测点的暂态零模电流相关性算法进行故障定位的方法。MATLAB仿真结果表明,该方法不受电压初相角、接地电阻、中性点接地方式的制约,具有较高的准确性和可靠性。

1 小电流接地系统故障暂态零模分量特征分析

小电流接地系统[6]中一旦某条出线发生接地故障,系统中就会产生暂态零模电压和零模电流。图1以工程中经常使用的中性点非直接接地系统为例,画出了发生频率最高的单相接地故障暂态零模电流分布情况。iof为零模虚拟电流源,C0S、C01、C02分别为发电机、线路1和线路2的对地电容,C1、C2分别为线路3上区间M-N、P-Q的对地电容。

小电流接地系统发生单相接地故障时,相当于在故障点处附加一个虚拟电源,在此电源的作用下,M-N区间暂态零模电流的方向为N流向M,P-Q区间暂态零模电流的方向为P流向Q,因此M-N与P-Q区间的零模电流初始极性相反,不具有相似性。M、N、P、Q处的零模电流满足:

式中:ic1为M-N区间的零模对地电容电流;ic2为P-Q区间的零模对地电容电流。

由于ic1相对i0M很小,可以忽略不计,由(1)式知M、N处的暂态零模电流近似相等,即i0M≈i0N;同理,P、Q处的暂态零模电流也近似相等,即i0P≈i0Q,两波形具有较高的相似性;因此,故障点F同侧两相邻检测点(M与N,P与Q)的暂态零模电流幅值基本相同,极性也相同,波形相似。

由式(1)、(3)可知,故障线路出线口处的零模电流为其它健全线路与故障线路对地零模电容电流之和;非故障线路和故障线路故障区段的零模电流为该线路的对地零模电容电流,因此故障线路出线口处的零模电流幅值明显大于非故障线路,且极性相反。又通常情况下故障发生在负荷侧,故障点F上游的线路长度大于故障点F下游的长度,其对地电容和线路电感也大于故障点下游,相应的频率前者小于后者。

故障线路上M、N、P、Q处零模电压具有相同的极性,且幅值相差不大,故零模功率的极性和幅值主要由零模电流决定,故障线路上零模功率分别是由故障点流向母线端与负荷端,非故障线路上则是由母线流向线路。在实际应用中,可以通过馈线终端装置(FTU)检测线路的零模分量,FTU通常安装在配电室或馈线上。

2 线路的相模变换

对于三相电力系统,各相之间存在着复杂的电磁耦合关系,采用Karenbuaer[7]变换对三相行波进行解耦,将相量分解为相互独立的模量,则:

式中:ia、ib、ic为三相系统中的电流值;i0、i1、i2为解耦后的模电流,其中i0即为三相线路的零模电流分量。

由(4)式可知,当系统正常运行时,ia、ib、ic为幅值相等,相位互差120°的交流电,零模电流为零;当系统发生故障时,零模电流不再为零。因此可将暂态零模电流作为故障定位启动信号。

3 基于暂态零模电流的故障选线方法

根据小电流接地系统单相故障的暂态零模分量分布特征,故障线路的零模电流幅值明显大于非故障线路,且极性相反。因此可通过线路零模电流的幅值和极性来确定故障线路,为了确保检测的准确性与灵敏度,设计选线算法如下:

式中:I0k为暂态零模电流的有效值;i0k(t)为第k条线路的暂态零模电流瞬时值。

式中:i0m(t)、i0k(t)分别为第m、k条线路的暂态零模电流瞬时值,选第m条线路为参考线路,依次与其它线路进行暂态零模电流瞬时值的内积运算,φkm>0表示参考线路与第k条线路极性相同,φkm<0表示极性相反。

若参考线路和其它所有线路都极性相反且参考线路的I0最大,则参考线路为故障线路;若参考线路仅和某一条线路极性相反且该线路的I0最大,则该线路为故障线路;若参考线路和其它所有线路极性都相同,则为母线故障。

4 基于暂态零模电流的故障定位方法

4.1 相关性理论

相关函数是描述各量相位、幅值关联性的一个重要的数字特征,可以用相关函数来判断2个函数的相似程度[8,9,10]。

设x(t)、y(t)是两个能量有限的信号,α是常数,τ是位移时间。

式(7)中必定存在一个最优的常数α,使得2个信号在均方误差最小的情况下获得最佳逼近,通过η2的平均值D度量两信号的相似度,则:

对D求关于α的一阶导数,令其一阶导数为零,求得D为最小值时对应的α值,此时有:

其中相关系数:

由以上相关性原理可见,Dmin随着ρxy的增大而减小,Dmin越小表明两个信号的相似度越高,对应的两波形越相似,当ρxy=1时Dmin=0,说明x(t)与y(t+τ)完全相似,对应波形也完全相同。ρxy=-1时Dmin=0,说明x(t)与y(t+τ)极性相反。因此相关系数ρxy可以准确地反映两个信号的相似程度,ρxy越大则两信号相似度越高,波形越相似。

为了便于计算机求解,令T=0,对相关系数ρxy进行离散化处理有:

4.2 故障定位原理

根据相关性原理,将式(11)中的x(n)、y(n)信号用检测到的暂态零模电流代替i01(n)、i02(n),

式中:i01(n)、i02(n)分别为相邻两检测点的暂态零模电流。

发生小电流接地系统单相故障时,非故障区段上两检测点检测到的暂态零模电流信号相似度很高,对应的波形也相似,相关系数|ρxy|>θ(θ为阀值)且接近于1;而检测点检测到的两个信号位于故障点两侧时,两信号的相似度很低,波形也相差较大,相关系数|ρxy|<θ且接近于0。在故障定位时,由母线侧向线路末端依次检测。图2为线路3发生接地故障的电力系统,故障定位时首先检测母线与M处的暂态零模电流并计算出相关系数ρxy,若|ρxy|<θ则故障发生在母线与M区段,若|ρxy|>θ则母线与M之间区段无故障,继续检测下一检测点N处的暂态零模电流,依次类推,直到|ρxy|<θ,检测到故障区段为止。故障定位流程如图3所示。

在求取零模电流相关系数ρxy时,检测到的相邻点暂态零模电流必须是同步的,因此要对采集信号装置进行同步化处理[11],确保算法的准确性。

5 仿真分析

采用MATLAB中的电力系统模块库工具箱,建立工程中最常见的中性点非直接接地系统单相接地故障模型,如图2所示。线路用分布参数模型表示,3条线路长度分别为5 km、11 km、20 km,故障发生在第3条线路L3的距母线10.5 km处,在0.04 s时发生短路故障,持续时间为0.04 s,RF为故障接地电阻,M、N、P、Q为线路3上的4个检测点。变压器为Δ/Y连接,一、二次侧电压分别为110 kV、10 kV,仿真参数设置:开始时间为0 s,终止时间为0.2 s,选择算法为Ode23tb。

图2中,K断开时为中性点不接地运行方式,K闭合时为中性点经消弧线圈接地运行方式,采用广泛应用的过补偿方式,补偿度为7%。采样率为1 MHz。

图4为3条线路在母线出线口处的暂态零模电流波形,可以明显地观察到,系统正常运行时线路的暂态零模电流为0;在0.04～0.08 s时线路出现暂态零模电流,系统存在接地故障,且母线出线口处故障线路的暂态零模电流幅值远大于非故障线路,而极性却相反。

图5为故障点F两侧的检测点N、P处的暂态零模电流波形,可以看出,故障点两侧的暂态零模电流极性相反,且非故障区段(故障点F的上游)N处的暂态零模电流幅值明显大于故障区段(故障点F的下游)p处的暂态零模电流。

电压初相角为0°,故障线路为L3,不同接地电阻时,中性点不接地和经消弧线圈接地方式下,以线路1作为参考线路,与其它线路暂态零模电流瞬时值的内积运算φkm和各线路的暂态零模电流有效值Iok的值如表1所示。

电压初相角为0°,故障区段为N-P,不同接地电阻时,中性点不接地和经消弧线圈接地方式下的相邻检测点暂态零序电流相关系数ρxy的值如表2所示。

接地电阻为5Ω,故障区段为N-P,不同电压初相角时,中性点不接地和经消弧线圈接地方式下的相邻检测点暂态零序电流相关系数ρxy的值如表3所示。

从表1—表3可以看出,在不同中性点运行方式、不同电压初相角、不同接地电阻下,都可利用本文提出的暂态零模电流的内积、有效值运算和相邻检测点的相关系数算法进行小电流接地系统的故障选线、定位,且不需要传输所有的暂态零模电流数据,具有信息传输量小,成本低,不受过渡电阻、接地方式制约,可靠性高等优点。

6 结语

本文提出了一种基于暂态零模电流,结合最大有效值、内积、相关性理论的故障选线、定位方法。该方法不受电压初相角、接地电阻、中性点运行方式的制约,仅需检测暂态零模电流,不需检测零模电压,信息传输量小,成本投入小;综合使用暂态零模电流内积的正负和最大有效值实现故障选线,将相邻检测点的暂态零模电流相关系数与阀值比较实现故障定位,具有较高的准确性和可靠性。

参考文献

[1]王玉梅,季涛.煤矿35 kV配电线路行波故障测距技术[J].电力自动化设备,2006,26(10):42-44.

[2]胥杰,张永健,曹旭,等.基于小波变换提取暂态直流分量故障选相原理的研究[J].华东电力,2011,39(2):262- 265.

[3]王显平,黄益华,唐昆明.基于多重小波变换的输电线行波故障测距[J].中国农村水利水电,2009,(5):112-115.

[4]王新超,张玉海.零序电流比例增量法在小电流接地故??障选线中的应用[J].电力系统保护与控制,2011,39(20): 125-130.

[5]王凤,康怡.基于脉冲信号注入法的小电流接地选线技术[J].电网技术,2008,32(15):90-93.

[6]张保会.电力系统继电保护[M].北京:中国电力出版社, 2004.

[7]季涛,孙波,苑倩倩.中性点非有效接地系统单相接地故障暂态零模电流获取方法研究[J].电力系统保护与控制, 2011,39(17):12-16.

[8]马士聪,徐丙垠,高厚磊,等.检测暂态零模电流相关性的小电流接地故障定位方法[J].电力系统自动化,2008, 32(7):48-52.

[9]王显平,黄益华,唐昆明.基于多重小波变换的输电线行波故障测距[J].中国农村水利水电,2009,(5):112-115.

[10]姚李孝,赵化时,柯丽芳,等.基于小波相关性的配电网单相接地故障测距[J].电力自动化设备,2010,30(1): 71-75.

暂态算法 篇5

针对中国能源资源分布与负荷地区不匹配,长距离、大容量、跨区跨流域水火互济的特高压同步联网送电已成为现代电网发展的趋势。随着大量高压直流输电和灵活交流输电系统的投运以及间歇性新能源发电的接入,电网安全稳定运行风险和控制难度日益增大,实现电网安全稳定的在线分析与控制已成为降低电网安全稳定运行风险的必然选择。

文献[1]提出了基于扩展等面积准则( EEAC) 的“在线预计算、实时匹配”在线紧急控制框架。文献[2]在实际电力系统中应用该类框架,定性评估首摆稳定性,并搜索切机策略。文献[3]详细介绍了分布式并行处理平台和区域稳定控制装置一体化的在线预决策暂态稳定控制系统。文献[4]介绍了预防控制与紧急控制一体化的在线预决策系统的技术框架。

近年来基于在线紧急控制框架[1,5,6,7,8,9,10]研发的在线紧急控制系统已应用到国内外多个电网,但在线策略优化算法方面仍存在2个不足: 一是只考虑了暂态功角稳定约束的优化; 二是只考虑了发电机切除或快关的优化。此外,暂态功角稳定问题的解决, 并不能确保电网暂态电压、频率安全稳定。单纯切除发电机,一方面可能难以确保电网暂态安全稳定, 另一方面可能控制代价过大。

暂态安全稳定紧急控制决策是在保持系统安全稳定的前提下,在含离散多类控制的候选措施空间中以控制代价最小为目标的优化问题。暂态安全稳定约束涉及暂态功角稳定、暂态电压和频率安全稳定。以控制代价为寻优目标进行穷尽式搜索策略的方法,显然不适合在线控制策略计算。为此,大量文献针对紧急控制策略搜索方法[11,12,13,14]、控制模型、择优目标选 择[15]、多类控制 措施间优 化进行研 究[16,17],但满足在线要求的暂态安全稳定紧急控制,在策略算法设计、搜索方向选择、多类控制措施协调搜索、提高在线控制策略计算速度等方面,仍值得探讨。

本文在文献[1,3-4]的基础上,提出直流调制、 切机、切负荷等不同交直流控制措施的控制性能指标和两阶段暂态安全稳定在线紧急控制策略计算方法,以协调优化多类控制约束和多类控制措施。在集群计算系统规模足够大的条件下,所提在线紧急控制策略计算方法最多只需要3次暂态安全稳定时域仿真即可计算出优化的暂态安全稳定在线控制策略,可满足在线计算要求。

1暂态功角稳定紧急控制性能指标

1.1直流系统功率紧急调制措施

根据EEAC理论[18],增加主导群负荷或余下群发电机驱动功率,有利于暂态功角稳定性。跨联于主导群和余下群间的直流系统,等同于主导群负荷和余下群发电机驱动功率,或主导群发电机驱动功率和余下群负荷。通过计算主导群中所有发电机母线与直流系统送、受端交流侧母线间的综合电气距离来确定直流系统功率紧急调制措施的控制性能指标。

采用式( 1) 和式( 2) 分别计算主导群中所有发电机母线与第i个直流系统送、受端交流侧母线之间的综合电气距离Zin. i和Zout. i。

式中: aj为第j个主导群发电机的暂态功角稳定参与因子; zin. i. j和zout. i. j分别为第j个主导群发电机与第i个直流系统送、受端交流侧母线间的电气距离, 采用二端网络输入阻抗评估节点间电气距离[19]; S为主导群中发电机台数。

若| Zin. i|小于| Zout. i| ,则表明第i个直流系统送端整流侧较受端逆变侧更靠近主导机群,第i个直流系统应紧急提升功率,反之,则紧急速降功率。若第i个直流系统确定为紧急提升功率,则将| Zin. i|的倒数作为其暂态功角稳定控制性能指标。若被确定为紧急速降功率,取| Zout. i|的倒数作为其暂态功角稳定控制性能指标。

1.2切机控制措施

对筛选出含有主导群中发电机的切机措施,考虑切机的控制代价及其暂态功角稳定参与因子[20], 采用式( 3) 计算主导发电机群中可选切除发电机中第j1个发电机的控制性能指标

式中:分别为第j1个含有主导群中发电机的切机措施的控制代价和该发电机的暂态功角稳定参与因子。

切机控制代价反映了切除发电机所付出的经济代价,可根据发电机损失的电量和重启动代价等因素人为设定。

1.3解列小电源送出电网控制措施

对于筛选出含有主导群中发电机的解列小电源送出电网措施,考虑切机的控制代价及其暂态功角稳定参与因子,采用式( 4) 来计算第k个解列小电源送出电网的暂态功角稳定控制性能指标Ta. k。

式中:分别为第k个含有主导群发电机的解列小电源送出电网措施中的第j2个发电机的控制代价和暂态功角稳定参与因子; Ng. k为第k个含有主导群中发电机的解列小电源送出电网措施中发电机的总数。

2暂态电压安全稳定紧急控制性能指标

2.1切负荷控制措施

考虑切负荷措施对暂态电压安全稳定薄弱节点的无功电压灵敏度和暂态电压安全稳定薄弱节点的稳定裕度,采用式( 5) 计算切除第j3个负荷对暂态电压安全稳定薄弱节点的综合影响因子采用式( 6) 和式( 7) 分别计算切除单个负荷和ml个负荷的控制性能指标

式中: Bv为暂态电压安全稳定薄弱节点数;为第i3个暂态电压安全稳定薄弱节点的暂态电压安全稳定裕度;为当前运行状态下切除第j3个负荷所连接节点的无功功率对第i3个暂态电压薄弱节点电压的无功电压灵敏度;为切除第j3个负荷的控制代价,反映了切除负荷的经济代价,可根据负荷量及其重要性等因素人为设定。

2.2解列小受端电网控制措施

考虑解列小受端电网对暂态电压安全稳定薄弱节点的无功电压灵敏度和暂态电压安全稳定薄弱节点的稳定裕度,采用式( 8) 和式( 9) 计算解列第k1个解列小受端电网措施对暂态电压安全稳定薄弱节点的综合影响因子和控制性能指标

式中:为第k1个解列小受端电网措施中小受端电网与外网连接的节点数;为第k1个解列小受端电网措施中小受端电网通过第j4个功率受入节点受入的有功功率;为当前运行状态下第k1个解列小受端电网措施中小受端电网第j4个功率受入节点的无功功率对第i3个暂态电压安全稳定薄弱节点电压的无功电压灵敏度;为第k1个解列小受端电网措施中小受端电网受入的总有功功率;为第k1个解列小受端电网措施的控制代价。

3暂态频率安全稳定紧急控制性能指标

3.1切负荷控制措施

考虑切负荷对暂态频率跌落安全薄弱节点的电气距离和暂态频率跌落安全薄弱节点的稳定裕度, 采用式( 10) 计算切除第j5个负荷对暂态频率跌落安全薄弱节点的综合影响因子采用式( 11) 和式( 12) 分别计算切除第j5个负荷和集中切除ml个负荷的暂态频率跌落安全控制性能指标

式中: Bfd为暂态频率跌落安全薄弱节点数;为第i4个暂态频率跌落安全薄弱节点的暂态频率跌落安全裕度;为切除第j5个负荷所连接的节点与第i4个暂态频率跌落安全薄弱节点间的电气距离; | zlmax|为切除第j5个负荷所连接的节点与暂态频率跌落安全薄弱节点间电气距离的模最大值;为切除第j5个负荷的控制代价。

3.2解列小受端电网控制措施

考虑解列小受端电网措施对暂态频率跌落安全薄弱节点的电气距离和暂态频率跌落安全薄弱节点的稳定裕度,采用式( 13) 计算控制措施的综合影响因子采用式( 14) 计算解列小受端电网措施的暂态频率跌落安全控制性能指标

式中:为第k2个解列小受端电网措施中小受端电网与外部电网连接的节点数;为第k2个解列小受端电网措施中小受端电网通过第j6个功率受入节点的有功功率注入;为受入总有功功率为第k2个解列小受端电网措施中小受端电网第j6个功率受入节点与第i4个暂态频率跌落安全薄弱节点间的电气距离; | zlmax' |为第k2个解列小受端电网措施中小受端电网第j6个功率受入节点与暂态频率跌落安全薄弱节点间电气距离的模最大值;为第k2个解列小受端电网措施的控制代价。

3.3切机控制措施

考虑切机对暂态频率安全薄弱节点的电气距离和暂态频率安全薄弱节点的稳定裕度,采用式( 15) 和式( 16) 分别计算切机对暂态频率安全薄弱节点的综合影响因子和暂态频率安全控制性能指标

式中: Bfu为暂态频率安全薄弱节点数;为第i5个暂态频率安全薄弱节点的暂态频率上升安全裕度;为第j7个切机措施中发电机节点与第i5个暂态频率安全薄弱节点之间的电气距离; | zgmax|为第j7个切机措施中发电机节点与暂态频率安全薄弱节点间电气距离的模最大值;为第j7个切机措施的控制代价。

3.4解列小电源送出电网控制措施

考虑解列小电源送出电网措施对暂态频率上升安全薄弱节点的电气距离和暂态频率上升安全薄弱节点的稳定裕度,采用式( 17) 和式( 18) 分别计算控制措施对暂态频率安全薄弱节点的综合影响因子和暂态频率安全控制性能指标

式中:为解列第k3个小电源送出电网措施中小电源送出电网与外部电网连接的节点数;为第k3个解列小电源送出电网中第j8个功率并网节点与第i5个暂态频率安全薄弱节点间的电气距离; | zgmax' |为第k3个解列小电源送出电网中第j8个功率并网节点与暂态频率安全薄弱节点间电气距离的模最大值;为第k3个解列小电源送出电网通过第j8个并网节点送出的有功功率;为送出的总有功功率;为第k3个解列小电源送出电网措施的切除第j9个发电机控制代价。

4暂态安全稳定紧急控制在线策略集群计算方法

本文将暂态安全稳定在线紧急控制策略优化分为2个阶段: 第1阶段计算暂态功角稳定的紧急控制策略; 第2阶段在已确定的暂态功角稳定的紧急控制策略基础上进行暂态电压和频率安全稳定的紧急控制策略计算。

4.1暂态功角稳定紧急控制策略集群搜索

4.1.1直流系统功率紧急升降调制及其挡位确定

首先,按式( 1) 和式( 2) 计算第i个直流系统送、受端交流侧母线与主导群所有发电机母线间的综合电气距离模值| Zin. i|和| Zout. i| ,并确定其功率紧急调制方向。

然后,筛选出控制性能指标与其最大值之比大于设定值 λ0的直流系统功率紧急调制措施,再依据设定的直流功率调制收敛精度 ΔPDC,将紧急功率提升的直流系统紧急提升空间 ΔPup和紧急功率回降的直流系统紧急回降空间 ΔPdown分别均分为M和N份,使得提升和回降搜索空间分别转换为由( 0,1, …,M) 和( - N,- N + 1,…,- 1,0) 的调制挡位构成。其中M = int ( ΔPup/ ΔPDC+ 0. 5 ) ,N =int( ΔPdown/ ΔPDC+ 0. 5) 。

4.1.2切机、解列小电源送出电网措施确定

首先,从备选的切机、解列小电源送出电网措施中筛选出含有主导群发电机的切机、解列小电源送出电网措施。

然后,计算筛选出的切机、解列小电源送出电网措施控制性能指标按式( 19) 和式( 20) 过滤,将满足过滤条件的第j1个切机措施、解列第k个解列小电源送出电网措施作为待组合的措施。

式中: Gmax为切机控制性能指标的最大值; Ta. max为解列小电源送出电网控制性能指标的最大值; λ1为设定值。

4.1.3直流系统功率紧急调制、切机和解列送出小电源措施的组合及评估

按式( 21) 计算直流功率紧急升降调制挡位与切机、解列小电源送出电网措施的枚举组合数Fd. g:

式中:为第ju个直流系统功率紧急提升的最高挡位;为第jd个直流系统功率回降的最低挡位; m为直流功率提升调制措施数; n为直流功率回降调制措施数; h为切机和解列小电源送出电网措施数。

为了确保各个组合措施执行后不会导致稳态频率安全问题,则需要考虑匹配相应的切负荷量,按式( 22) 确定各组合切负荷匹配量 ΔPL. m:

式中: f',K,β 分别为电网当前运行状态下考虑故障清除及组合措施实施后的电网稳态频率、功率频率静态特性系数、网损系数; ΔPG为组合措施中所有切机及解列小电源送出电网的送出有功功率之和; ΔPG0和 ΔPL0分别为因故障清除退出的发电机有功功率之和、负荷有功功率之和; f0为电网当前频率; fd为故障后满足电网安全要求的稳态频率下限。

以振荡中心为割集将电网划分为两部分,将与余下群发电机连接在一起的切负荷措施、解列小受端电网措施作为用于匹配的备选措施,且以其暂态电压安全稳定紧急控制性能指标由大到小的顺序选择用于匹配的负荷控制措施。

将匹配的切负荷措施加入控制措施组合( { 直流调制,切机,解列,匹配切负荷} ) 作为一个计算任务,先按控制代价由小到大形成调度队列,再对含有直流调制措施的计算任务,按直流调制的控制性能指标与调制量的乘积之和除以调制量之和的比值由大到小进行排序,最后对比值相同的计算任务,按调制量之和由小到大进行排序,提交并行集群系统进行暂态安全稳定量化评估。若优先调度的控制措施组合的暂态功角稳定裕度大于0,则将相应的控制措施组合作为暂态功角稳定的在线紧急控制策略, 提前中止后续计算任务。若所有计算任务均暂态功角失稳,则认为搜索不到暂态功角稳定的在线紧急控制策略。

4.2暂态电压和暂态频率安全稳定策略集群搜索

4. 2. 1暂态电压安全稳定措施确定

按可投切电容电抗器所连母线对暂态电压安全稳定薄弱节点的无功电压灵敏度,与投切电容电抗器措施的无功电压灵敏度最大值之比大于设定值 λ2,筛选确定待组合可投电容器、可切电抗器措施。

按切负荷、解列小受端电网的控制性能指标,与所有切负荷、解列小受端电网措施的控制性能指标最大值之比大于设定值 λ3,筛选确定待组合的切负荷措施、解列小受端电网措施。

4.2.2暂态频率跌落安全稳定措施确定

按备选的切负荷、解列小受端电网的控制性能指标,与备选切负荷措施、解列小受端电网控制性能指标最大值之比大于设定值 λ4,筛选确定待组合的暂态频率跌落切负荷、解列小受端电网安全措施。

4.2.3暂态频率上升安全稳定措施确定

按备选的切机、解列小电源送出电网的控制性能指标,与备选切机、解列小电源送出电网的控制性能指标最大值之比大于设定值 λ5,筛选确定待组合的暂态频率上升切机、解列小电源送出电网措施。

4. 3暂态电压和暂态频率安全稳定措施的组合及评估

按式( 23) 计算待组合的可投切电容电抗器、切负荷、解列小受端电网及切机、解列小电源送出电网的枚举控制措施组合数Fc. g. l。

式中: Bb为可投切电容电抗器措施所连节点数; Lj为第jc个节点可投电容器或可切电抗器组数; lv为暂态电压安全稳定切负荷和解列小受端电网措施数; lf为暂态频率跌落安全切负荷和解列小受端电网措施数; hf为暂态频率上升安全切机和解列小电源送出电网措施数。

同样为了确保各个组合措施执行后不会导致稳态频率安全问题,也需要考虑匹配相应的切负荷量, 按式( 24) 确定各组合的匹配切负荷量 ΔPL. m':

式中: f″,K',β'分别为电网当前运行状态下考虑故障清除、暂态功角稳定策略及组合措施实施后的电网稳态频率、功率频率静态特性系数和网损系数; ΔPG'和 ΔPL分别为包括暂态功角稳定策略在内的控制措施组合中所有切机及解列小电源送出电网的送出有功功率之和、切负荷及解列小受端电网的受入有功功率之和。

将匹配的切负荷措施加入控制措施组合( { 可投切电容电抗器,切机,解列,匹配切负荷} ) 作为一个计算任务,先按控制代价由小到大形成调度队列, 再对含有投切电容电抗器措施的计算任务,按投切电容电抗器对暂态电压安全稳定薄弱节点的无功电压灵敏度与投切量的乘积之和除以投切量之和的比值由大到小进行排序,最后对比值相同的计算任务, 按投切量之和由小到大进行排序,提交并行集群系统进行暂态安全稳定量化评估。

若优先调度的控制措施组合可确保电网暂态安全稳定,则将该控制措施组合作为该故障的暂态安全稳定在线紧急控制策略,提前中止后续计算任务。 若所有计算任务均不满足暂态安全稳定,则认为搜索不到该故障的暂态安全稳定在线紧急控制策略。

5策略计算仿真分析

以单环全接线“三华”电网为例,研究方式中复奉直流双极输出功率为2 ×2 500 MW,锦苏直流双极输出功率为2 ×2 800 MW,其他直流均双极满载运行。川渝断面外送功率2 260 MW,华中特高压电网送华北电网1 777 MW,华北特高压电网送华东电网2 828 MW。

暂态电压跌落安全的电压门槛值均取0. 8( 标幺值) ,持续时间设为1 s,暂态频率安全偏移值取 ± 0. 3 Hz,持续时间设为0. 5 s。直流系统功率紧急升降调制、切机和切负荷动作时间,均按检测到故障后延时0. 2 s考虑。复奉、锦苏直流系统按双极满载运行考虑最大提升空间,其他直流按1. 1倍的短时过载能力考虑。

对于洪沟—板桥N - 2故障系统暂态失稳的紧急控制策略搜索如下: 选择复奉、锦苏、龙政、宜华四回直流系统作为直流系统功率紧急升降调制候选措施,经直流控制性能指标计算及筛选,确定复奉、锦苏直流系统功率提升,最高提升依次为4挡和2挡。

按发电机性能指标参与因子筛选后确定切机措施为向家坝等5台可切发电机。

暂态功角稳定控制策略搜索的直流系统功率紧急调制、切机措施组合数为480( 3 × 5 × 25) ,搜索的最优策略为: 复奉直流系统提升2 × 700 MW,锦苏直流系统提升2 × 200 MW。

对于仍存在暂态电压和频率安全稳定的问题,经控制性能指标计算及筛选,确定木里、水洛、遂宁、色尔站作为候选切负荷措施,可能组合数为81( 34) ,进一步的控制策略为: 木里站切负荷200 MW。

灵绍直流双极闭锁后电网暂态失稳的紧急控制策略搜索的结果,同样验证了算法的有效性。详细的算例策略搜索及在线紧急控制策略集群计算流程见附录A。

6结语

本文基于暂态安全稳定量化分析理论,提出了直流系统功率紧急调制、切机、切负荷和无功补偿装置投切以及局部电网解列的控制性能指标计算方法。基于控制性能指标计算技术和集群计算技术, 提出了适合电力系统暂态安全稳定在线紧急控制策略集群计算方法。仿真实例验证了所提控制性能指标计算和在线策略计算方法的有效性。基于本文所开发的适用于特高压交直流多直流系统第二道安全防线在线控制策略的快速计算模块,已应用于宁夏协调控制系统。

暂态算法 篇6

电力负荷监测方法大致可以分为侵入式和非侵入式2类[1]。非侵入式负荷监测(non-intrusive load monitoring,NILM)只需要在电力供给的入口处安装非侵入式的监测设备,就可以实现对整个系统内部负荷的在线监测。相较于传统的侵入式负荷监测方法,NILM不中断设备供电,易于被用户所接受;不需要安装大量的检测设备,节省了购买、安装和维护这些硬件设备所需要的投资和时间,是未来电力负荷监测的重要发展方向之一[2,3,4]。需求侧管理是智能电网最重要的组成部分之一[5],对电力公司而言,NILM有助于电力公司了解电力用户负荷的构成,加强负荷侧管理,通过引导用户合理消费、合理安排负荷的使用时间等措施达到调节峰谷差、降低网损等目的;有助于改善电力负荷的预测精度,为电力系统仿真分析、系统规划提供更准确的数据和模型。对电力用户而言,通过NILM系统,用户能够深入了解企业或者家庭各种用电设备的能源消耗及构成等各种信息,从而为采取有针对性的节能措施提供指导,进而减少不必要的能源开销。

为了提高NILM系统在负荷监测与辨识方面的效率和准确性,研究人员提出了很多的监测方法[2,3,6,7,8]。负荷投切等过程中产生的暂态特征相比用电设备的稳态特征,能为负荷辨识提供更多的有用信息,因此利用负荷暂态特征来辨识负荷及其运行状态是NILM发展的一个主流方向[4,6,9,10]。理论上,依据电气设备的暂态特征,NILM能够实现对整个变电站、建筑物内部负荷集群的有效分解与准确分析,可在一定程度上克服仅利用负荷稳态特征进行辨识的局限性。基于暂态特征进行负荷辨识,准确掌握电气设备的投切时刻是其关键,为此,本文提出了一种有效的暂态事件检测算法———基于滑动窗的双边累积和(cumulative sum,CUSUM)暂态事件自动检测算法,它能根据相关信号(如电压、电流、有功、无功等)准确检测到电气设备投切等引起的系统暂态过程,具有算法简单、抗干扰能力强的特点。

1 非参数化的CUSUM变点检测算法

把负荷投入或者切除等暂态过程定义为一个暂态事件,负荷投切等引起的暂态事件检测可以归结为变点检测问题。所谓变点是指在一个序列或者过程中,在某个未知时刻,序列或者过程的某些统计特性发生了变化[11,12,13]。变点检测就是对变点是否发生以及发生时刻等作出估计。电力负荷的投切等暂态过程往往导致采样序列的某些统计特性发生突变,常见的检测算法为CUSUM变点检测方法。

CUSUM算法是统计过程中常用的算法,它最初是由Page在1954年提出的,随后许多学者对此算法进行了深入的研究。CUSUM算法的理论基础是序贯分析原理中的序贯概率比检验(sequential probability ratio test,SPRT)理论。CUSUM设计思想是对样本数据信息加以累积,将过程的小偏移累积起来,达到放大的效果,从而提高检测过程中对小偏移的灵敏度。由于使用方便、判断准则简单、易于操作等原因,CUSUM算法在工业质量控制、自动故障监测、经济、金融等方面应用广泛。它通过研究信号均值、方差等信息来判断系统是否发生了变化,可分为参数化算法和非参数化算法2类[11,12]。在NILM的应用场合,由于电力系统中负荷是不断变化的,通常很难获得负荷的功率、电压、电流等随机变量的概率分布情况,故非参数化CUSUM算法更为合适。其主要思想是:当监测量的CUSUM明显比正常平稳运行条件下的平均水平高的时候,就意味着系统发生了变化,由此可判断系统出现了暂态过程。由于负荷的投切等通常伴随着功率水平、电压与电流有效值等信号平均值的变化,此时就可以用非参数化的CUSUM变点检测算法确定负荷的投切时刻,其原理如下。

考察随机时间序列X={x(k)},k=1,2,…。x(k)=μ+AkI(k<τ)+(D+Bk)I(k≥τ),其中,Ak和Bk为随机变量,且E(Ak)=E(Bk)=0,D为平均值的变化量,τ为变点发生的时间,I(·)为指示函数,μ为X的均值。定义非参数化CUSUM算法中统计函数g为:g0=0gk=max(0,gk-1+sk{)(1)式中:sk=xk-μ0-β,μ0为变点发生前的平均值,通常假定是已知或可估计的;β为外界引入的噪声,当xk的变化量大于β时才认为负荷水平有可能发生突变。

记检测延迟时间为d。当gk>0时,采样序列发生了变化,但若尚未达到阈值h,则令检测延迟时间d=d+1,否则令d=0。当gk>h时,一个暂态事件被检测出来,采样序列突变发生的时刻可以倒推估计为τ^=k-d。

由于电力系统负荷水平有可能增加也有可能减少,本文采用双边CUSUM算法,其具体算法如下。定义)(2)(3)

假设稳态时的采样序列均值为μ0,下面简要说明双边CUSUM的原理。

1)当序列处于稳态时,CUSUM统计量gk+和gk-是一个在0附近随机波动的变量,即以0为均值的随机变量。

2)若序列发生突变,出现正向偏移,偏移后的序列均值上升为μ1>μ0+β,那么,这个向上的、正的偏移就会不断累积到统计量gk+中,gk+不断增大;反之,如果序列发生负向偏移,偏移后的序列均值下降为μ1<μ0-β,那么向下的、负的偏移就会不断累积到统计量gk-中,gk-不断增大。

3)当gk+>0,但尚未达到阈值h时,令检测延迟时间dk+=dk++1,否则令dk+=0。当gk+>h时,一个均值增加的暂态事件被检测出来,采样序列突变的发生时刻τ可以倒推估计为τ^=k-dk+。类似地,当gk->0,但尚未达到阈值h时,令检测延迟时间dk-=dk-+1,否则令dk-=0。当gk->h时,一个均值减少的暂态事件被检测出来,采样序列突变的发生时刻τ可以倒推估计为τ^=k-dk-。

上述双边CUSUM算法的优点是算法简单、易于实现,但该算法的应用前提是已知采样序列均值μ0。对于负荷不断变化的实际电力系统,均值μ0也是不断变化的,为此本文通过引入滑动窗口来对CUSUM算法进行改进。

2 基于滑动窗的双边CUSUM变点检测

实际电力系统的负荷总在不断变化,故式(2)、式(4)中的μ0是不断变化的,这对直接采用双边CUSUM算法准确确定暂态事件造成了障碍。为此,本文提出了基于滑动窗的双边CUSUM算法,对此加以改进,该算法演进过程如图1所示。

在WM窗口,算法计算采样序列均值μ0;在WD窗口,算法判断系统是否发生了暂态事件。当没有检测到暂态事件时,WM和WD窗口向右滑动到新的采样点。当检测到一个暂态事件时,新的WM和WD窗口滑动到算法检测到变点的时刻。为了防止检测过程中某些暂态事件被漏检,当前WM窗口的末端与上一WD窗口减去检测延迟d后对齐。这是因为,当暂态事件发生在WD检测窗口的尾部时,gk+和gk-可能没有足够的时间上升到阈值h,此时平均值窗口WM的末端被特意调整放置在暂态事件可能的发生时刻之前,这样可以在新的窗口进一步验证是否真正有暂态事件发生,从而降低了漏检的可能性。

双边CUSUM算法中的2个参数h和β的物理意义不是很明确,它们可以用其他更具有物理意义的参数来描述,如暂态事件检测器能够检测到的最小突变量Δmin和检测到此突变量的最大允许延迟Tmax。如图2所示,针对幅度为Δmin的阶跃突变,要求算法在时间Tmax以内检测到。

由此可得,阴影部分面积h=(Δmin-β)Tmax,而在离散采样系统中:h=(4)式中:Tmax=NmaxTs;Nmax为允许的最大延迟采样点数;Ts为采样间隔;β与测量的噪声有关。

由于实际电力系统中,负荷消耗的有功功率一般是缓慢上升而不是阶跃变化的,上述阴影部分并不是一个矩形,此外噪声水平等因素也比较难确定,特引入系数λ1和λ2来修正阴影部分面积,式(4)可以进一步改写为:(5)式中:λ1≤1;λ2≥1。

具体的滑动窗双边CUSUM算法的程序流程如图3所示。

根据对Δmin和Tmax以及抗干扰能力的要求,可以对h等参数进行初始化,然后根据不断读入的数据块来判断系统是否发生了扰动。其中,数据块的采样点数为N,WM和WD窗口中采样点数分别为Nm和Nd,算法检测到暂态事件的时刻为kevent,暂态事件发生的起始时刻的估计值为kstart。Nm的大小表示平均值计算窗口样本个数的多少,它的取值要能尽量代表正常状态下的特征。Nd表示检测窗口样本个数的多少,Nd值不能太大,否则会由于检测窗口过长而引入较长的检测延迟。

3 算例仿真研究

本节通过算例仿真来验证所述暂态事件检测算法的有效性。在PSCAD/EMTDC中构建了一个改造过的IEEE 30节点测试系统,其基本参数见文献[14],网络拓扑结构如图4所示。在节点14处多种负荷通过一个33kV/13.8kV的变压器挂接在33kV母线14下,它简单模拟一个电力工业用户的内部系统。用户负荷包括容量为1.0 MVA的绕线式感应电动机M1(有软启动控制器),容量为0.3 MVA的鼠笼式感应电动机M2(无软启动控制器),多个感性(RL)负荷。由于负荷投切等暂态过程中,通常会伴随着功率水平的变化,可以根据功率水平的变化来判断是否发生了暂态事件。下文通过检测T1点的聚合功率来检测负荷的投切过程。NILM通过测量节点14相关的电气量可以进一步辨识工业用户内部的负荷构成及运行情况等有用信息,如哪些电气设备是耗电大户,哪些电气设备的运行管理存在问题等,以此来制订相应的节能改造措施。此外也为电力公司监测大型工业用户内部负荷种类及运行状态信息提供了便利。

3.1 算例1

额定电压为13.8kV、额定功率为1 MVA、功率因数为0.8的RL负荷启动过程中测量点的有功功率如图5所示,采样频率为1 000 Hz。RL负荷启动前测量点处负荷水平约为1MW,RL负荷启动时暂态功率分量会叠加在原有功率之上。RL负荷在k=420采样点处(约0.42s时)投入,测量点处的有功功率大约发生0.8 MW的跃变。

在MATLAB中实现上文所述的滑动窗双边CUSUM算法。在PSCAD/EMTDC中获得相关仿真数据,然后导入MATLAB中进行检验。为验证算法的抗干扰能力,在PSCAD/EMTDC中获得的相关仿真数据上叠加了标准差为2%的白噪声。本算例中,假定只检测功率大于0.8 MW的负荷波动,即检测算法的分辨率为0.8MW,波动幅度小于0.8 MW的负荷变化都不认为是真正的负荷投切引起的突变。在其他应用场合,综合考虑负荷水平及仪器的测量精度等信息,按照类似的原则可以选择恰当的分辨率。由上文分析可知,可以选择Δmin=0.8 MW,β=0.02。取λ1=0.8,λ2=2.0,Nmax=100,Nd=100,Nm=50,此时可得h=60.8。算法仿真结果如图6所示。在滑动窗1~3数据段内,由于无暂态事件发生,g+,g-,d+,d-均等于0或者很接近于0。在滑动窗4数据段尾段时,RL负荷投入,负荷水平发生变动,如图6所示,g+和d+开始累积,不断增长。但是,由于突变点发生在滑动窗检测窗口的尾部,虽然在此窗口中确实发生了暂态事件,但由于没有足够的时间使g+累积到阈值h,在此窗口内不能判断是否发生了暂态事件。如图6(a)所示,根据滑动窗4的d+,可以计算滑动窗口4内暂态事件可能发生的起始时刻为:k=450-27=423。为了准确检测到暂态事件发生的起始时刻,把滑动窗5检测窗口的起始点修正到k=423点处。在滑动窗5内,如图6(b)所示,g+不断累积,在k=kevent时g+累积到h,可以确定系统发生了暂态事件,此时暂态事件的起始时刻可以倒推为kstart=kevent-d+=523-98=425,非常接近发生突变的真实时刻。若采用普通的CUSUM算法,对于上述暂态事件发生在检测窗口尾部的情况将产生较大误差。算例中,由于负荷功率水平的跃变量等于能够检测到的最小突变量Δmin,故直到接近滑动窗口5的尾部g+才累积到h。

3.2 算例2

本算例校验负荷功率水平发生-1 MW跃变,即大于Δmin时的检测效果。图7所示为某个RL负荷切除时(k=350,t=0.35s)测量点的有功功率。

此算例中,CUSUM算法基本参数同算例1。在滑动窗1~3数据段内,由于无暂态事件发生,g+,g-,d+,d-均等于0或很接近于0。在滑动窗4数据段,RL负荷切除,负荷水平发生变动,如图8所示,g-和d-开始累积,不断增长,在k=kevent时,g-累积到h,可确定系统发生了暂态事件,此时暂态事件起始时刻可倒推为kstart=kevent-d-=437-82=355,同样非常接近突变事件的真实发生时刻。

类似地,在感应电动机负荷投入或切除等过程中,测量点的有功功率水平也会发生变化,此时根据g+和d-也能判断系统是否发生突变以及突变发生的初始时刻。根据采集到的负荷暂态特征信息可以进一步辨识负荷类型、参数及运行状态等信息。

4 结语

对NILM系统而言,负荷投切等过程中产生的暂态特征信息包含了丰富的有用信息。NILM能够利用此暂态信息将检测对象的负荷构成精确地分解到各类用电设备,进而根据电气设备的负荷特性形成整个检测对象的负荷特性,以提高后者的准确性。现在NILM已经与智能电表技术密切结合,来增强后者的检测功能。为了能充分地利用此暂态信息,本文提出了一种基于滑动窗的双边CUSUM暂态事件自动检测算法,此算法简单、抗干扰能力强。数值仿真表明它能够自动检测、辨识负荷投切等引起的暂态过程,从而为NILM后台高级应用程序打下基础,具有一定的工程实用价值。

摘要：非侵入式负荷监测(NILM)是未来电力负荷监测的重要发展方向之一。不同类型电力负荷在投切过程中,通常会表现出独特的暂态特征。据此,NILM能够克服利用负荷稳态特征信息进行负荷辨识的局限性,实现对整个变电站、建筑物内部负荷集群的分解与分析。提出了一种信号预处理算法——基于滑动窗的双边累积和(CUSUM)暂态事件自动检测算法,它能根据相关信号准确检测到负荷投切等引起的暂态过程、发生时刻等重要信息,并触发相关程序把暂态信息记录下来,然后送给后台高级应用程序作进一步处理。算例仿真表明了所述算法的有效性。

暂态算法 篇7

随着智能电网技术的发展,新的发电、用电和储能设备逐渐接入,新的控制策略也越来越多地应用于电网。这些设备和策略一方面改善了电网的运行性能,另一方面也产生了新的故障形式和相应的动态过程。为此,有必要研究并发展准确、快速的电磁暂态仿真工具,提升智能电网设计和分析能力[1-2]。

为满足不断增长的仿真计算需求,需利用先进的计算机技术不断改善电磁暂态仿真工具[3-4]。近年来,高性能并行计算技术有了长足的进步,特别是图像处理器(GPU)在通用科学计算领域的成功应用[5-7],为实现大规模系统的高效仿真提供了新的解决思路。GPU是计算机显示卡的核心处理单元,其开发初衷是为了实现高速的图像渲染。在通用计算领域,GPU最适合求解数据密集、计算量大且并行度高的计算问题,而电磁暂态仿真正属于此范畴。 利用GPU强大的计算能力有望大幅度提升电磁暂态仿真的效率。目前,已有研究通过调用CULA[8]等商业软件包并行求解节点电压方程,实现了最简单的电磁暂态GPU仿真[9-11],但其并未提出适用于GPU的电磁暂态并行仿真算法。

现有电磁暂态并行仿真算法采用“分解—协调” 策略实现系统级并行[12-16],称为粗粒度并行算法。 由于粗粒度并行算法只能在系统层面分解计算量, 可并发运行的线程和计算核心数目十分有限[17],因此,其并不适合在GPU中运行。GPU要求将求解过程分解为大量可并行的简单运算,实现运算级并行。因此,实现适用于GPU电磁暂态仿真算法的难点在于,如何将带有开关过程和复杂控制的电磁暂态仿真求解分解为大量可并行的简单运算,以增加可并发运行的线程和计算核心数目,实现算法的细粒度并行化。针对GPU高性能计算环境的特性,本文在文献[18]提出的改进电磁暂态程序(EMTP)算法基础上,设计并实现了适用于GPU的电磁暂态细粒度并行算法。 文中介绍了改进EMTP算法的总体流程,指出其适用于并行求解, 并设计了适用于GPU的运算级并行算法,包括基于单指令多数据流(SIMD)和基于共享内存的运算级并行。针对改进EMTP算法中各求解步骤的特点,采用两种运算级并行算法,实现了改进EMTP算法的细粒度并行,并通过算例验证了此细粒度并行仿真算法的正确性和高效性。

1改进EMTP算法流程

EMTP算法是电磁暂态仿真的常用方法,具有原理简单、仿真结果准确、效率较高等特点,但开关动作和控制系统求解会造成仿真流程碎片化,不适用于GPU并行求解。文献[18]提出了改进EMTP算法(主要流程见图1),其将EMTP算法中复杂的开关处理逻辑转化为开关周期的迭代求解。在每个开关周期内,求解流程是一致的,这便有利于采用运算级并行策略,实现仿真流程的细粒度并行求解。

2适用于GPU的细粒度运算级并行策略

本节设计了两种细粒度运算级并行策略:一是脱离仿真对象的物理背景,将主要计算转化为SIMD的向量计算模式[19],该转化过程称为“向量化”;二是根据仿真对象的物理属性分解计算,使用共享内存处理各线程之间的耦合关系。

2.1基于SIMD的运算级并行策略

当一组GPU线程访问连续存储于设备内存中的数据并执行相同操作时,其具有较高的计算效率,如图2所示。

从图2可见,SIMD的本质是向量运算。实现SIMD的一个方法是将计算转化为如下向量形式:

式中:x,y,z∈Rn;*表示任意向量运算符。

简单且常见的向量运算有向量加法和向量Hadamard积,即

向量Hadamard积的定义如下:设x,y,z∈Rn,若z=xy,则有zi=xiyi(i=1,2,…,n)。其中xi,yi,zi分别为x,y,z的对应元素。

由于向量加法和Hadamard积运算中,各元素的操作相互独立,故可在GPU上启动大量线程并行求解,其中每个线程计算向量的一个元素。

2.2基于共享内存的运算级并行策略

对于无法转化为向量运算的求解步骤,可根据仿真对象的物理属性分解计算。例如,在执行与元件有关的操作时,可并发与元件个数等量的线程,其中每个线程负责完成一个元件的计算过程。但是该计算过程中可能存在耦合,导致各线程无法完全独立执行,主要有以下两种情况。

1)不同元件的计算过程存在先后顺序,如元件A的输出是元件B的输入,这种情况在控制系统求解中十分普遍。此时,必须在元件A求解完成后再求解元件B,否则将导致计算结果错误。

2)两个元件需对同一块内存执行写操作,例如在形成节点注入电流向量时,需要将连接在同一个节点上不同电气元件的诺顿等值电流加到等值电流向量的同一个元素上。此时,需要避免两个线程同时写同一块内存,否则可能导致累加结果不正确。

对于以上两种情况,均可使用GPU共享内存实现线程间通信和原子操作来正确求解[19],故称这类细粒度并行方法为基于共享内存的运算级并行。

3适用于GPU的电磁暂态细粒度并行算法

改进EMTP算法各环节的计算特征,以及据此可采用的细粒度并行求解策略总结如表1所示。

3.1计算电气元件诺顿等值电流

由于不同端口数目元件的诺顿等值方程维数不同,因此在求解前,需要对电气元件按端口数分组。

对于一般的电气元件,设端口数为l,则其诺顿等值电流向量通常可以表示为上一时步端口电压和电流的线性组合[20],即

式中:ine,ib,ub∈Rl,分别为诺顿等值电流向量、支路电流向量和端口电压向量;p,q∈Rl×l,为系数矩阵。

设系统中存在s个端口数为l的电气元件,令

式中:Ib,Ub,Ine∈Rsl;P,Q,Yne∈Rsl×l;下标1,2,…,s为元件编号。

设P,Q的第k列分别为P(k),Q(k),则有:

这便实现了电气元件诺顿等值电流计算的分组SIMD向量化,可在GPU上并行求解。

3.2形成节点注入电流向量

在形成节点注入电流向量时,需要将每个电气元件的诺顿等值电流值加到注入电流向量的对应元素上。因此在求解时,只需并发与电气元件数目相同的线程,每个线程对应一个电气元件,线程之间通过共享内存通信。其具体实现为:在共享内存中建立一个数组,用于存储节点注入电流向量;每个线程将其对应元件的诺顿等值电流值加到相应的向量元素上,使用原子操作以避免不同线程同时对同一块内存进行操作。

3.3求解节点电压方程

节点电压方程为线性方程组,由于每个开关周期内系统导纳矩阵不变,因而可采取先求逆矩阵、每时步计算矩阵向量乘法的求解策略。其中,后者可以利用SIMD并行求解以提高效率。具体来说,分别采用以下两个并行数学库调用GPU进行求解: 利用CULA对导纳矩阵求逆[8];利用CUBLAS程序库实现矩阵向量乘法[21]。

3.4计算电气元件内部变量

对于端口数为l的电气元件,其支路电流计算方程为[20]:

式中:yne∈Rl×l,为诺顿等值导纳矩阵。

在按照端口数分组之后,类似3.1节,端口电流的求解也可统一为:

式中:Y(k)ne为Yne的第k列。

需要说明的是,部分复杂元件可能需要求解支路电流以外的内部变量,例如发电机需要求解机械方程以获得转子角和角速度,这些计算亦可通过类似方式实现向量化,在此不再赘述。

3.5求解控制系统

基于共享内存的控制系统求解流程见图3。

图中:Nctrl为控制元件总数;Ncou为检测控制元件是否可以求解的计数器。在每次求解时,循环检测每个控制元件的输入节点值是否已知,若已知便求解此控制元件,并将其输出节点标记为已知状态。当所有元件均求解后,便完成了控制系统求解。在求解时,需并发Nctrl个线程,每个线程求解一个控制元件。不同线程间通过共享内存通信,以满足不同控制元件间的依赖关系。具体实现如下。

1)在共享内存中建立一个长度等于控制节点数的一维标志位数组,用于标记节点的状态:若节点已求解,则相应位置标记为“1”;否则标记为“0”。

2)每个线程在求解控制元件前,首先检测该元件所有输入节点的标志位是否都为“1”:若是则可读取输入节点信号,执行求解,并将输出节点的标志位设为“1”;否则继续检测。

3.6开关动作时刻预测和校正

根据文献[18],采用开关动作时刻的线性预测和校正方法:

式中:α为迭代阻尼系数;G(·)为系统求解占空比的等效多项式;τ为开关动作时刻向量;下标f表示预测值;上标j表示第j次迭代值;m表示第m个开关周期。

式(9)和式(10)所示运算均为向量运算,因此可直接采用SIMD并行求解策略。

4算例测试与分析

4.1算例及仿真环境设置

采用三相脉宽调制(PWM)AC-DC变流器系统对细粒度并行仿真算法的正确性和效率进行测试, 其电气拓扑如图4所示。参数如下:直流电容为32 000μF,滤波电阻为0.001 Ω,滤波电感为0.6mH,载波频率为5kHz,电源频率为50Hz, 电源幅值为0.4kV,电源初相角为0rad,直流负载为2Ω。控制系统采用VQ控制模式。仿真所用硬件系统的主要参数如下:CPU型号为Intel?Xeon? E5645,核心数为6个,主频为2.4 GHz,内存为24GB;GPU型号为nVIDIA?TeslaTMC2070, 核心数为448(分为14个流处理器),主频为1.15 GHz,内存为6GB。

4.2正确性测试

设置仿真积分步长为10μs。测试变流器系统在以下两种情景中的动态响应:①交流电压跌落故障,电源幅值于0.05 s时从0.4 kV跌落到0.15kV,并在0.1s时恢复;②直流侧短路故障,变流器直流侧于0.05s时发生电阻性短路故障,短路电阻为0.05Ω,并在0.1s时恢复。

分别采用3种方式对测试系统进行仿真:①利用PSCAD/EMTDC作为标准参照系统;②在CPU上采用改进EMTP串行仿真算法[18](以下简称CPU串行程序);③在GPU上采用细粒度并行仿真算法(以下简称GPU并行程序)。

选取交流A相电流和直流电压,在两种情景下,对比3种方式的仿真波形,得到结果分别如图5和图6所示。可见,GPU并行程序可准确模拟PWM变流器在交流电压跌落以及直流侧短路时交流电流与直流电压的变化过程,其仿真结果与CPU串行程序和PSCAD/EMTDC高度吻合。对于其他仿真结果,三者也有类似的吻合关系,在此不再赘述。

4.3效率测试

在图4所示测试系统的基础上,通过增加并联PWM变流器(及其控制系统)的数目以构造不同规模的仿真算例,测试交流侧电压跌落故障下系统规模不同时并行算法的仿真效率。由此得到CPU串行程序和GPU并行程序的单时步仿真平均耗时如表2和图7所示。

可见,当PWM变流器数量为15个时,GPU并行程序的效率已超过CPU串行程序。随着系统规模增加,GPU并行程序相对于CPU串行程序的加速比增大。当PWM变流器数量为90个时,得到最大加速比为6.74。可以预计,在不超过GPU线程分配数目情况下,增加系统规模将获得更大加速比。

进一步,分别测试电气系统和控制系统的仿真耗时,得到如表3、图8和图9所示测试结果。

分析以上结果,可以得到如下结论。

1)由于PWM变流器算例的电气系统规模较小,故并行计算加速效果并不明显。但在测试中已可以看出,并行算法的求解时间仅仅随电气系统规模增加而稍有增加,在变流器数量为45个时,已经低于快速增加的串行程序仿真耗时。

2)由于不同变流器的控制系统具有很好的并行度,随着控制系统规模增大,细粒度并行仿真耗时基本不变,从而获得了较为明显的加速效果。

由表3可知,系统规模较大时,并行程序求解电气与控制系统耗时不到总体耗时的一半。进一步测试可知,开关动作时刻预测与校正环节,特别是其中的系统导纳矩阵求逆环节占用了将近60%的求解时间,其显著影响并行仿真算法的效率。由于在3.3节所示细粒度并行仿真算法中,并未对导纳矩阵求逆做任何细粒度处理,只是调用数学库CULA并行求解,因此,在进一步研究中需考虑如何提高这一步骤的求解速度,如对CULA进行优化配置,或换用更加高效、专门性的数学库。此外,也可以尝试在此环节中仅对导纳矩阵做简单分解,将部分计算转移到求解节点电压方程环节中,以达到总体计算量最小的目的。

5结语

结合GPU并行计算平台的特点,本文在改进EMTP算法的基础上设计并实现了适用于GPU的电磁暂态细粒度并行算法。仿真测试表明,此细粒度并行算法能够在保持仿真正确性的同时,提高仿真效率。本文对使用GPU进行带有开关过程和复杂控制的大规模电力系统快速电磁暂态仿真进行了有益探索,并指出算法的改进方向为求解节点电压方程中导纳矩阵的细粒度处理。

摘要：智能电网技术的发展需要快速电磁暂态程序(EMTP),而日益广泛应用的图像处理器(GPU)为电磁暂态仿真提供了高效的仿真环境和平台。文中首先提出了细粒度并行算法的运算级并行策略,即基于单指令多数据流(SIMD)的运算级并行策略和基于共享内存的运算级并行策略。随后,设计了应用这两种并行策略的改进电磁暂态细粒度并行算法。三相脉宽调制(PWM)变流器仿真测试表明,适用于GPU的细粒度并行算法能够在保证仿真正确性的同时,显著提高仿真效率,从而验证了基于GPU的细粒度并行仿真算法适用于带有开关过程和复杂控制的大规模电力系统快速电磁暂态仿真应用的可行性。