基于sod-ps的时滞电力系统稳定性判别方法
【技术领域】
[0001] 本发明涉及一种基于S0D-PS的时滞电力系统稳定性判别方法。
【背景技术】
[0002] 广域测量系统（Wide-AreaMeasurementSystem，WAMS)的出现给大规模互联电力 系统稳定分析与控制的发展带来新的契机。基于WAMS提供的广域信息的互联电网低频振 荡控制，通过引入有效反映区间振荡模式的广域反馈信号，能够获得较好的阻尼控制性能， 其为解决互联电网中的区域间低频振荡问题，进而提高系统的输电能力提供了新的控制手 段，具有良好而又广泛的应用前景。
[0003] 广域信号在由不同通信介质（如光纤、电话线、数字微波、卫星等）组成的WAMS 通信网络中传输和处理时，存在几十到几百毫秒间变化的通信延时。时滞是导致系统控 制律失效、运行状况恶化和系统失稳的一种重要诱因【[l]WuHX，TsakalisKS，HeydtG T.Evaluationoftimedelayeffectstowide-areapowersystemstabilizerdesign. IEEETrans.PowerSyst.，2004, 19(4) : 1935-1941.】。因此，利用广域测量信息进行电力系 统闭环控制时，必须计及时滞的影响。
[0004] 发明【[2]牛新生，叶华，王春义，等.基于Pad6近似的时滞电力系统特征值计 算与稳定性判别方法.201210271783. 8: [P].】利用Pade近似多项式逼近时滞环节，进而计 算系统最右侧的关键特征值，并判断系统的时滞稳定性。
[0005] 发明【[3]叶华，王燕燕，刘玉田·基于EI⑶的大规模时滞电力系统特征值计算 方法· 201510055743. 3.中国，201510055743. 3[P]·】提出了一种基于显示I⑶(Explicit IGD，EIGD)的大规模时滞电力系统特征值计算。利用计算得到的系统最右侧的关键特征 值，可以判断系统在固定时滞下的稳定性。这些时滞稳定性判别方法，均需要通过多次扫描 [0. 1，2. 5]Hz低频振荡频率范围内、靠近虚轴的关键特征值，才能判断系统的时滞稳定性。
【发明内容】
[0006] 本发明的目的就是为了解决上述问题，提供一种基于S0D-PS(伪谱离散化解算子 (Pseudospectraldiscretizationofsolutionoperator,S0D-PS))的时滞电力系统稳定 性判别方法，该方法只需要计算解算子离散化矩阵的模值最大的一个特征值或一对共辄特 征值，就可以判断出系统在固定时滞下的稳定性。具有计算量小，判别准确的特点。
[0007] 为了实现上述目的，本发明采用如下技术方案：
[0008] 基于S0D-PS的电力系统时滞稳定性判别方法，包括如下步骤：
[0009] 步骤（1):建立时滞电力系统模型；依据时滞电力系统模型的特征值与时滞电力 系统模型的解算子特征值之间的关系，将计算时滞电力系统模型的特征值转化成计算解算 子的特征值；从而将判断时滞电力系统稳定性的问题转化为计算解算子的模值最大的特征 值问题；
[0010] 步骤（2):采用伪谱方法对解算子进行离散化，得到解算子的离散化矩阵；
[0011] 步骤（3):采用序贯法或子空间法来计算步骤（2)得到的解算子的离散化矩阵的 模值最大的特征值μ;
[0012] 步骤（4):根据特征值μ的大小来判断时滞电力系统的稳定性。
[0013] 所述步骤（4)的步骤为：
[0014] 若特征值μ的模值大于1，则时滞电力系统处于小干扰不稳定的状态；
[0015] 若特征值μ的模值等于1，则时滞电力系统处于临界稳定的状态；
[0016] 若特征值μ的模值小于1，则时滞电力系统处于渐进稳定的状态。
[0017] 所述步骤（1)的时滞电力系统模型如下：
[0018]
[0019] 式中：为电力系统的状态变量向量，η为系统状态变量总数。t为当 前时刻。〇〈τ'τ2〈···〈τΓ··〈τ 时滞环节的时滞常数，其中最大的时滞为τn。 i,.e：rx"(i=(U, .,,，7W)为系统状态矩阵。Ax(t)为t时刻系统状态变量的增量，ΔΧα-τι) 为时刻系统状态变量的增量，Μ(?)为t时刻系统状态变量导数的增量。Δχ(0)为系 统状态变量的初始值（即初始条件），并简写为…
[0020] 式（1. 1)表示的线性化系统的特征方程为：
[0021]
[0022] 式中：λ为特征值，ν为特征值对应的右特征向量。
[0023] 所述解算子的定义如下：
[0024] 解算子T(h)定义为将初值条件@转移到h(转移步长，0彡h彡τJ时刻之后时 滞电力系统解分段的线性算子。
[0025]
[0026] 其中，s为积分变量，Θ为变量，和分别为〇和h+θ时刻时滞电力系 统的状态。。
[0027] 所述时滞电力系统模型的特征值与解算子特征值之间的关系：
[0028]由谱映射定理可知，解算子T(h)的特征值μ与时滞电力系统的特征值λ之间具 有如下关系：
[0029]
[0030] 式中：〇 (T(h))表示解算子的谱。
[0031] 所述步骤（2)的步骤如下：
[0032] 与解算子T(h)对应的、标准基形式的离散化矩阵ΤΜιΝ表示如下：
[0033]TMN -TM+TM，N*(lNn_UN) *UM，N(1.5)
[0034] 式中：
[0039] 式（1.5)中，M和N为给定正整数，1_为版阶单位阵（对角线元素为1，其他元素 为0)，上标-1表示矩阵求逆运算。
[0040] 式（1.6)中，Q为给定正整数，
，1MX1为元素全为1的Μ维 列向量，1?$为（Q-2)M阶单位阵，1"为11阶单位阵，0 (Μ+1)Χ% (Μ+1)ΧΜ阶零矩阵，
…?为Kronecker积运算。矩阵ΤΜ为高度稀疏矩阵，并 与时滞电力系统状态矩阵沁(/=()，…，"?>无关。
[0041] 式（1.7)中，
为QMXN阶零矩阵，
[0042] 式（1.8)中，
i=0, "'πιΑ的元素完全由 拉格朗日系数决定。
[0043] 式（1.9)中，
，i= 0, "'m，R的元素完全由拉格朗日系 数决定。
[0044] 所述步骤（3)的步骤如下：
[0045] 假设在第k次迭代时，需要计算ΤΜιΝ与向量 的乘积，
具体步骤如下：
[0046]步骤（3-1):将向量vk按列压缩为一矩阵…，i= 1，…，QM+1。相应地，有：vk=vec(Vk)，其中，vec( ·)为将矩阵压缩为列向量的运算。
[0047] 步骤（3-2):计算/ =C/uv *v*eM-1;
[0048]步骤（3-3):计算/=(/#
[0049]步骤（3-4):计算wk=TM*vk+T'M,N*qk。
[0050] 所述步骤（3-2)的步骤如下：
[0051] 将式（1.8)代入，可得：
[0053] 式中：上标T表示矩阵转置。k为第k次迭代；m为时滞的个数。
[0054] 由上述分析可知，要计算？14，首先要计算為尽1，i=0，···，πι，然后再进行求 和，最后再压缩为一个Νη维的列向量。
[0055] 值得注意的是，為哪…,4^+1]可以稀疏实现，从而减小计算量，提高 计算效率。
[0056] 所述步骤（3-3)的步骤如下：
[0057] 将式（1.9)代入后，
没有显示表达。因而，这 里采用迭代算法来计算^= (INn_UN)Vs。求解过程中，涉及矩阵向量乘法运算b= (INn_UN)y，其中ftelR*''1。
[0058] 首先，将向量y按列压缩为一矩阵
只eKBxl,i= 1，…，No
[0059] 进而，可以得到b= (INn_UN)y的稀疏实现步骤如下：
[0061] 值得注意的是，篇)作=扃* ...,抑]可以稀疏实现，从而减小计算量，提高计 算效率。