基于SOD-IRK的大规模时滞电力系统特征值计算方法

牟倩颖, 叶 华, 刘玉田

(电网智能化调度与控制教育部重点实验室(山东大学), 山东省济南市 250061)

0 引言

随着互联电网规模的不断增大,区域间低频振荡的问题愈加突出[1-2]。基于广域信息的阻尼控制器可以有效抑制区间低频振荡[3-5]。然而,广域通信时滞是影响控制器性能的重要因素。在小干扰稳定性分析[6]和广域阻尼控制器设计[7]中需要考虑广域通信时滞的影响[8-9]。

在频域,时滞被转化为指数项,系统的特征方程成为超越方程。直接求解该方程得到系统的特征值变得异常困难[10]。现有研究通常采用Padé近似有理多项式来逼近时滞环节[11]。然而,在大时滞情况下,Padé近似相位逼近误差较大。近年来,有学者提出了基于谱算子(包括无穷小生成元和解算子)离散化的时滞电力系统特征值计算方法,用于准确求解系统的部分关键特征值。无穷小生成元离散化方法[9,12-14]能够计算给定位移点附近的部分特征值。为了得到系统的全部机电振荡模式,需要进行多次特征值计算。解算子离散化(solution operator discretization,SOD)方法的优势在于通过一次计算便可得到系统最右侧或阻尼比最小的部分特征值。文献[15]提出基于解算子隐式龙格—库塔离散化(SOD with implicit Runge-Kutta,SOD-IRK)的时滞电力系统特征值计算方法,但是该方法不能直接应用于分析大规模时滞电力系统。并且,文献[15]仅考虑了Radau-ⅡA方法求解常微分方程,并未研究其他隐式龙格—库塔(IRK)方法的可用性和性能。

鉴于此,本文提出了一类基于SOD-IRK的大规模时滞电力系统关键特征值计算方法。一方面,采用多种IRK离散化方案离散解算子,得到高度稀疏和结构化的解算子离散化矩阵;另一方面,通过充分利用系统增广状态矩阵和解算子离散化矩阵的稀疏特性,高效地计算大规模时滞电力系统阻尼比最小的部分关键特征值,以避免维数灾问题。

1 问题描述

时滞电力系统的状态方程为[12]:

(1)

(2)

式(1)对应的特征方程为超越方程,即

(3)

式中:λ和v分别为特征值和相应的右特征向量。

(4)

(5)

式中:w∈Rl×1为辅助向量;In为单位矩阵。

为了避免直接求解式(3)和式(4),本文首先采用多种IRK方法对式(1)对应的解算子进行离散化,然后计算解算子离散化矩阵的部分特征值,进而得到时滞电力系统阻尼比最小的部分关键特征值。

2 基于SOD-IRK的特征值计算方法

2.1 基本原理

设X:=C([-τmax,0],Rn×1)是由区间[-τmax,0]到n维实数空间Rn×1映射的连续函数构成的巴拿赫空间,并赋有上确界范数。式(1)所示时滞电力系统的解算子T(h):X→X定义为将X上t时刻初值条件φ映射为t+h时刻系统状态的线性算子,即

(T(h)φ)(t)=Δxh(t)=Δx(t+h)

(6)

式中:h为转移步长,0≤h≤τmax。

在频域,T(h)的特征值μ与时滞电力系统的特征值λ之间存在如下指数关系[16]:

μ=eλhμ∈σ(T(h)){0}

(7)

式中:σ(·)表示谱;表示集合差运算。类似于Cayley变换[17],T(h)将s平面虚轴右侧的特征值λ映射到z平面单位圆之外。

基于SOD-IRK的特征值计算的原理是,首先对T(h)进行IRK离散化,得到其有限维离散化矩阵TNs,这是SOD-IRK方法的核心和关键;其次,计算TNs模值最大的部分特征值μ;最后,由式(7)得到时滞电力系统最右侧的部分关键特征值λ。

T(h)的IRK离散化一般性步骤如下:首先,将时滞区间[-τmax,0]分为长度为h=τmax/N的N个子区间;然后,利用p阶s级IRK的s个横坐标c1,c2,…,cs对每个子区间进行离散化,如图1所示;最后,在集合ΩNs每一个离散点处,根据式(6)推导Δxh与φ之间表达式的系数矩阵,即得TNs。

(8)

图1 离散点集合ΩNsFig.1 Discrete point set ΩNs

本文首先详细说明利用IRK的Radau-ⅡA方法对T(h)进行离散化[18],然后对T(h)的其他IRK离散化方法进行简要总结。

2.2 解算子的Radau-ⅡA离散化方案

1)当t∈[-τmax,-h]时,Δxh(t)为初始状态φ。利用解算子的转移特性,可以得到t=-jh+cqh(q=1,2,…,s;j=2,3,…,N)处Δx(t+h)的估计值为:

(9)

令Δx(0),Δx(1)∈RNsn×1分别为定义在区间[-τmax,0]和[h-τmax,h]上的系统的离散状态变量,其表达式为:

(10)

于是,式(9)写成向量形式:

(11)

2)当t∈[-h,0],t+h-τi<0,i=1,2,…,m时,利用p阶s级IRK的Radau-ⅡA法求解如下常微分方程:

(12)

从而,可推导得到Δx(1)与Δx(0)之间的显式表达式为:

(13)

式中:RNs∈Rsn×sn,ΣNs∈Rsn×Nsn。

(14)

式中:A∈Rs×s为IRK方法Butcher表(A,b,c)的元素[19]。Radau-ⅡA方法的Butcher表具有如下特点:①A可逆;②0

3)综合式(11)和式(13),可得Δx(1)与Δx(0)之间的关系式为:

Δx(1)=TNsΔx(0)

(15)

(16)

(17)

式中:TNs∈RNsn×Nsn为解算子T(h)的IRK离散化矩阵。

2.3 解算子的其他IRK离散化方案

解算子的其他IRK离散化方案与Radau-ⅡA离散化方案的不同之处,可总结为以下两点。

1)当cs≠1时,需要增加离散点cs+1=1以得到每个子区间最右端点处系统状态的估计值。

2)当c1=0时,每个子区间左端点与相邻子区间的右端点重合,即θj-1+c1h=θj+csh或θj-1+c1h=θj+cs+1h,其中j=2,3,…,N。因而,需要去掉这些重复的离散点。

以Lobatto-ⅢA/C离散化方案为例,其系数(A,b,c)有如下特点:①0=c1

(18)

(19)

(20)

类似地,可处理Gauss-Legendre、单对角隐式龙格—库塔法(SDIRK)、Radau-ⅠA和Lobatto-ⅢB等离散化方案(见附录B)。此外,对于p阶IRK方法,特征值计算的局部计算误差为Ο(hp)[19]。

3 SOD-IRK方法的稀疏实现

3.1 旋转—放大预处理

旋转—放大预处理[20],一方面将时滞电力系统阻尼比小于ζ(=sinθ)的部分特征值λ乘以旋转因子e-jθ,从而将它们变换到s平面虚轴右侧,进而利用式(7)将它们映射到z平面单位圆之外;另一方面,将z平面的特征值进行α次乘方,以增加它们之间相对距离。

预处理后,TNs变为TN′s:

(21)

式中:N′=τmax/(αh),其中 表示向上取整。

(22)

(23)

为了保证预处理后时滞为实数,式(23)对时滞做了必要的近似,其对特征值的影响分析可参考文献[20]。

在旋转—放大预处理中,参数α取2或3,即可显著增加T(h)特征值之间的相对距离,改善迭代特征值算法的收敛性。此外,α的取值还需要满足:αh≤τmax。

参数θ的选取依据包括:待求部分关键特征值的分布情况和SOD-IRK方法的应用目的。首先,希望得到的关键特征值应当位于等阻尼比线ζ的右侧。其次,如果希望计算系统较多的关键特征值以设计控制器,则θ应取较大的值,如θ≥5.74°;如果希望计算系统最右侧的少量特征值用以快速、可靠地判断系统的小干扰稳定性,则θ应取较小的值,如1.72°或2.87°。

3.2 稀疏特征值计算

采用隐式重启动Arnoldi算法来计算TN′s模值最大的r个特征值μ′。在计算过程中,计算量最大的操作是迭代生成Krylov向量。下面以Radau-ⅡA离散化方案为例,阐述由第j个Krylov向量qj∈CN′sn×1高效计算第j+1个Krylov向量qj+1的方法。

qj+1=TN′sqj

(24)

由于TN′s的次对角为单位矩阵,qj+1的后(N′-1)sn个元素可通过直接将qj的元素向后位移sn得到。前sn个元素qj+1(1:sn)的计算可以分解为式(25)和式(26)。

z=ΣN′sqj

(25)

(26)

式中:z∈Csn×1为中间向量。

令qj=vec(Q),其中Q∈Cn×N′s,vec(·)表示将矩阵按列重排为列向量的操作。利用Kronecker积的性质[21],式(27)给出了式(25)的高效实现。

(27)

(28)

2)式(28)可改写为以下两步:

(29)

在计算之前,仅对D0做一次稀疏三角分解。于是,式(29)的计算仅是一些稀疏矩阵与向量相乘和两个三角方程的求解。

3.3 牛顿校验

(30)

4 算例分析

4.1 16机68节点系统

16机68节点算例系统的详细数据参见文献[22]。考虑在G2和G5上分别安装广域阻尼控制器[20],反馈时滞和输出时滞分别为τf1=150 ms,τo1=90.5 ms,τf2=70 ms和τo2=37.4 ms。所有测试均在8 GB RAM、3.4 GHz CPU台式计算机上进行。SOD-IRK方法和牛顿校验的收敛精度均为10-6。

首先,验证不进行旋转—放大预处理和牛顿校验,采用各种IRK离散化方案时,TNs逼近T(h)的准确性。各参数分别为h=0.005 s,N=30,α=1和θ=0°。以2级Radau-ⅡA,Radau-ⅠA,Gauss-Legendre,SDIRK,Lobatto-ⅢA,Lobatto-ⅢB和Lobatto-ⅢC离散化方案为例,SOD-IRK方法分别能够准确估计T(h)模值大于0.744 7,0.861 7,0.727 0,0.777 3,0.937 8,0.903 3,0.938 2和0.903 6的特征值μ,分别对应实部大于-59.41,-30.16,-65.45,-50.93,-19.95,-20.66,-12.97和-20.32的特征值λ。部分结果见附录C。

通过对比分析可知,Radau-ⅡA,Radau-ⅠA,Gauss-Legendre和SDIRK离散化方案具有较好的特征值估计精度。此外,尽管各种IRK方法能够准确计算的特征值数目不同,这部分准确特征值都足以分析系统的小干扰稳定性。

接着,分析旋转—放大预处理的影响。在应用旋转—放大预处理(α=2和θ=8.63°)后,各种SOD-IRK方法都具有快速计算系统全部(r=15)机电振荡模式的能力,如图2所示。

图2 预处理后的系统特征值估计值及其准确值λFig.2 Estimated value and exact value λ of eigenvalue after pre-processing

最后,对比分析各种SOD-IRK方法计算系统全部机电振荡模式的效率。由表1可知,所有SOD-IRK方法的计算时间都在10 s以内,且单次迭代时间均小于基于解算子的伪谱离散化(SOD-PS)方法[20](其中切比雪夫多项式阶数为3)。当参数变化时,随着解算子离散化矩阵维数的增大,计算时间随之增加。相比之下,前4种SOD-IRK方法的计算效率略低。

表1 16机68节点系统采用各种SOD方法的效率比较Table 1 Efficiency comparison among various SOD methods in 16-generator 68-bus system

4.2 华北—华中特高压互联电网

在某水平年下华北—华中特高压互联电网上进一步验证各种SOD-IRK方法对大规模时滞电力系统的处理能力。系统含有33 028个节点,2 405台同步发电机,16条高压直流输电线路。考虑在蒙岱海1号和鲁东海7号上分别装设广域阻尼控制器以抑制区域间低频振荡[13]。反馈和输出时滞分别为τf1=120 ms,τo1=100 ms,τf2=100 ms和τo2=80 ms。系统状态变量的维数为n=80 577,代数变量为l=162 718。

首先,利用各种SOD-IRK方法计算系统阻尼小于5%的r=20个特征值。图3给出了2级Radau-ⅡA和Lobatto-ⅢA方法的计算结果。各参数分别设为h=0.006 s,N=20,α=2和θ=2.87°。图中,部分对时滞灵敏度较大的特征值的估计误差较大,需要通过牛顿校验予以消除。此外,系统仍然存在2对负阻尼区间振荡模式和众多失稳和临界失稳局部模式。通过对比分析可知,各种SOD-IRK方法的精度基本相当。

表2给出了不同SOD-IRK方法和SOD-PS方法计算20个特征值的IRA迭代次数和计算时间。通过对比分析可知:①SOD-IRK的单次迭代计算时间远小于SOD-PS;②Lobatto-ⅢA方法的计算效率最高(如表2中红色数据所示)。

图3 华北—华中特高压电网特征值估计值及其准确值λFig.3 Estimated value and exact value λ of eigenvalue in North China-Central China ultra-high voltage power grid

表2 华北—华中特高压电网采用各种SOD方法的效率比较Table 2 Efficiency comparison among various SOD methods in North China-Central ultra-high voltage China power grid

5 结语

本文提出了一类基于SOD-IRK的大规模时滞电力系统特征值的高效计算方法。

2)Radau-ⅡA,Radau-ⅠA,Gauss-Legendre和SDIRK方法具有较好的计算精度,而Lobatto-ⅢA和Radau-ⅡA方法具有较高的计算效率。

本文提出的基于SOD-IRK的大规模时滞电力系统特征值计算方法中,系统所有的状态变量都需要离散化,导致解算子离散化矩阵维数较大,特征值计算时间长。因此,考虑只对与时滞相关的状态变量进行离散化以高效求解时滞电力系统的特征值是后续的研究方向。

本文受到山东大学青年学者未来计划项目(2016WLJH06)支持,特此感谢!

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