CN101699448B - 一种电力系统暂态稳定的分布式仿真方法 - Google Patents

Abstract

一种电力系统暂态稳定的分布式仿真方法属于电力系统暂态稳定分析技术领域，其特征在于，使用基于节点撕裂方法的分区方法使边界协调方程的变量维数降至基于联络线的分区方法的一半，以提高仿真效率；使用边界摄动法得到较为准确的初始Jacobi逆矩阵，以减少迭代次数；综合采用预估下一时步边界节点电压值、多时步同时协调和动态更新Jacobi逆矩阵的方法来加速计算进程；使用基于逆Broyden拟牛顿法的协调算法大幅降低边界协调方程的求解次数，从而提高了仿真的效率。该技术方案数据接口简单，数据传送量小，简化了系统模型，具有接近牛顿法的超线性收敛速度，在广域通信网络中也可以实现超实时的仿真，更加便于在实际电力系统中应用。

Description

一种电力系统暂态稳定的分布式仿真方法

技术领域

一种电力系统暂态稳定的分布式仿真方法属于电力系统分布式仿真技术领域，更具体地说，是涉及到一种电力系统暂态稳定分布式仿真的算法。

背景技术

我国区域电网互联正在不断加强，为了保持互联大系统安全高效运行，有必要对全网进行一体化和无简化的仿真分析。分布式暂稳仿真可有效整合各类异构计算资源，并在保持区域调度中心计算独立性的同时获得全网一致收敛的仿真结果，是实现互联电网一体化仿真分析的有效途径。

目前，关于分布式暂稳仿真的研究主要集中在以下几个方面：

1、分解协调算法和模型的研究；

2、通信方式和通信效率的研究；

3、提高仿真效率的加速方法的研究；

4、分布式暂稳仿真评价指标的研究；

5、分布式暂稳仿真鲁棒性的研究。

总的来说，已有学者初步探讨出几条可行的分布式暂稳仿真技术路线，并采用一些有效的加速方法来调高仿真效率。然而，目前存在的最大的困难在于仿真过程中的通信时间过长。这是因为调度中心间分布式暂稳仿真需要在广域电力通信网络中进行，该网络带宽有限(一般为2Mbps)，且通信延时较高(一般为几十到一百多毫秒)，而分布式暂稳仿真所需通信次数往往较多，这必然导致仿真效率偏低。因此，在保证算法正确性的基础，研究如何减少仿真所需通信次数，是分布式暂稳仿真算法研究的难点所在。

发明内容

本发明的特征在于，在保留数据接口简单、通信量小等优点的基础上，大幅提高了协调计算过程的收敛速度，使其达到超线性，因此仿真所需通信次数较少，即使在高延时广域通信网络中亦可获得较高仿真效率。同时，本算法含有基于节点撕裂方法的分区方法和基于边界摄动的Jacobi逆矩阵生成方法。节点撕裂分区方法可以使协调方程的变量维数降为基于联络线的分区方法的一半，从而提高仿真效率；而边界摄动法可以获得较为准确的初始Jacobi逆矩阵，从而减少了迭代次数。本发明所述的分布式暂稳仿真算法依次含有以下步骤：

步骤(1)初始化

步骤(1.1)采用节点撕裂方式进行网络分区。以图1所示的两分区电力系统为例，S₁、S₂为两个区域电网(以下称子分区，其计算服务简称分区服务)，S_B是由边界节点构成的上级电网(以下称协调侧，其计算服务简称协调服务)，S_1In、S_2In为两个子分区的内部节点，B为边界节点集合，B₁、B₂、B₃是同一边界节点在分区1、分区2和协调侧的表示，是子分区的边界节点在协调侧对应的虚拟节点，i表示子分区序号。基于这种分区方式，边界节点满足：

$\left\{\begin{array}{c}{u}_{\mathrm{Bi}}\∠{\mathrm{\θ}}_{\mathrm{Bi}}=u_{\tilde{B}i}\∠{\mathrm{\θ}}_{\tilde{B}i}\\ P_{\mathrm{Bi}}+{\mathrm{jQ}}_{\mathrm{Bi}}+P_{\tilde{B}i}+{\mathrm{jQ}}_{\tilde{B}i}=0\end{array}\right.,i=\mathrm{1,2}---\left(1\right)$

式(1)中，u、θ为边界节点的电压幅值和相角，P、Q为边界节点的注入有功和无功。

步骤(1.2)设定仿真始末时间和仿真步长。以t_s表示仿真起始时刻，t_e表示仿真结束时刻，h表示仿真步长。t_s通常取0s，t_e根据仿真需要选择，一般取5s~80s，h通常取0.01s~0.02s。

步骤(1.3)设定同时协调的多时步数目。以m表示同时协调的多时步数目(以下简称多时步数)，m根据仿真需要灵活选取。多时步数是对常规逐步仿真模式的重大改进，能够一次性的完成m个时步的仿真，从而提高仿真效率。

步骤(1.4)设置网络故障信息。根据仿真的需要选择网络故障，用于仿真的故障需具备典型性。

步骤(2)设置当前的仿真时刻，并生成对应的仿真时段信息。

时段是与多时步数对应的概念，常规的逐步仿真模式下，每个时段的长度为一个时步h。而在多时步模式下，每个时段含有m个时步。用参数k表示当前仿真的时段序号，第一个时段k取1。对于第k个时段，其始末时间表示为[t_k，t_k+m·h]。

步骤(3)采用边界摄动法获得初始Jacobi逆矩阵。

电力系统暂稳仿真本质上是求解一组微分代数方程：

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

其中，x是状态量向量，y是代数量向量，f是微分方程组，g为代数方程组。

在分布式暂稳仿真过程中，各子分区独立进行本分区的暂稳仿真，协调侧则负责各分区边界信息的匹配。如果每一时步，各子分区内部动态网络潮流达到平衡，且边界节点满足式(3)：

$\left\{\begin{array}{c}{u}_{B1}\∠{\mathrm{\θ}}_{B1}=u_{B2}\∠{\mathrm{\θ}}_{B2}\\ P_{B1}+{\mathrm{jQ}}_{B1}+P_{B2}+{\mathrm{jQ}}_{B2}=0\end{array}\right.---\left(3\right)$

则可判定全网动态潮流一致收敛。

从协调侧看，分布式暂稳仿真过程中，边界节点的电压和注入功率差可表示为隐函数方程即边界协调方程：

$F\left(U_k\right)=\left[\begin{array}{c}\mathrm{\Δ}{P}_k\\ \mathrm{\Δ}{Q}_k\end{array}\right]=\left[\begin{array}{c}{P}_{B1}+P_{B2}\\ Q_{B1}+Q_{B2}\end{array}\right]=0---\left(4\right)$

式(4)中U、ΔP、ΔQ分别表示边界节点电压向量、注入有功功率差、注入无功功率差，下标k表示当前仿真时段。对式(4)进行迭代求解即可完成第k个时段的仿真计算。

初始Jacobi逆矩阵对式(4)的求解效率影响重大，因此需要采用有效的技术手段获得。边界摄动法能够保持分布式仿真通信接口简单，传送数量量少的优点，并能获得足够精确的初始Jacobi逆矩阵。其实施流程见步骤(3.1)~(3.4)：

步骤(3.1)将非线性方程F(U)＝b在U₀做线性化，得到JU₀＝b。其中，U₀表示电压初值，J是U＝U₀处的Jacobi矩阵，b是F(U)的函数值。

步骤(3.2)将U₀的每一维分别加上一个微小量e＝10^-4，得到n个向量[U₀₁+e，U₀₂，…，U_0n]^T，[U₀₁，U₀₂+e，…，U_0n]^T，…，[U₀₁，U₀₂，…，U_0n+e]^T。将这些向量代入F(U)得到b₁，b₂，…，b_n，即：

步骤(3.3)将JI₀＝b式(5)逐列相减，得到：

e·J·I＝[b₁-b b₂-b…b_n-b]    (6)

其中，I是与J同阶的单位矩阵。从式(6)可解出Jacobi逆矩阵J¹，下文用A来表示J¹。

式(5)是边界摄动法的基本原理，在实际应用时只要向需要摄动的分区传递摄动信息即可，边界节点电压无变化的分区无需传递；发往同一分区的多组边界量信息彼此无耦合关系，因此可以一次性打包发送以减少通信开销。摄动向量只含边界节点信息，数据量很小，因此边界摄动环节对分布式暂稳仿真耗时影响甚微。

步骤(4)预估边界节点电压值。

当开始一个新时段的第一次仿真计算时，如果有3组以上的历史数据(当前仿真时步之前的边界节点电压值)，则采用式(7)所示的二次外推方法获得边界节点电压预估值，可以使仿真效率大为提高；如果历史数据不足则采用上一时步电压值作为本时步的预估值。

$U_{c+i}=\left[\begin{array}{ccc}{t}^2& t& 1\end{array}\right]\·\left[\begin{array}{ccc}0.5& -1& 0.5\\ 0.5& -2& 1.5\\ 0& 0& 1.0\end{array}\right]\·\left[\begin{array}{c}{U}_{c-2}\\ U_{c-1}\\ U_c\end{array}\right],t=\mathrm{1,2},...,m---\left(7\right)$

其中，U_c-2，U_c-1，U_c分别是第c-2，c-1，c时步的边界节点电压值，t表示预估的时步序号。在多时步情况下，分别预估当前时刻往后顺延的第1~m个时步的边界节点电压值，并将其组合成列向量，即为多时步情况下的边界节点电压预估值。

步骤(5)采用逆Broyden拟牛顿法计算边界节点电压的修正量，并更新Jacobi逆矩阵。

步骤(5.1)已知A_k，协调服务向分区服务发送边界节点电压U_k，并从分区服务获得边界节点注入功率，经计算得到本时段本次协调计算的偏差向量

F(U_k)＝[ΔP，ΔQ]^T；

步骤(5.2)利用U′_k＝U_k-A_kF(U_k)计算U′_k，U′_k是修正后的边界节点电压值；

步骤(5.3)协调侧向分区发送U′_k并获得分区反馈的F(U′_k)，从而得到s_k＝U′_k-U_k，y_k＝F(U′_k)-F(U_k)；

步骤(5.4)利用 $A_k^{\′}=A_k+(s_k-A_k{y}_k)\frac{{\left(s_k\right)}^T{A}_k}{{\left(s_k\right)}^T{A}_k{y}_k}$ 更新A′_k，并令新的U_k＝U′_k，A_k＝A′_k；

步骤(6)判断本时段仿真是否结束。

判断是否满足收敛条件F(U_k)＜10^-4，如果收敛则转步骤(7)，如果不收敛则返回步骤(5)，并令k＝k+1；

步骤(7)判断全部仿真是否结束。

判断是否满足仿真结束条件，如果t_i≥t_e则仿真结束，输出相关信息；反之，返回步骤(2)，更新时间信息t_i＝t_i+(m+1)×h，并令k＝1。

附图说明：

图1节点撕裂分区方法

图2分解协调流程图

图3东北-华北-华中联网系统中协调算法的收敛速度

图4高速局域网络中算法的加速比

图5广域通信网络中(40ms延时)算法的加速比

具体实施方式：

以东北-华北-华中联网系统(简称NNN系统)为例说明专利说明书中叙述的算法。该系统的相关参数为：

表1  测试系统参数

仿真平台的先关参数为：

表2  仿真平台参数

基于表1和表2的测试系统和仿真平台进行分布式暂稳仿真，具体实施步骤如下：

步骤(1)初始化

步骤(1.1)采用节点撕裂方式进行网络分区。按照地理分布情况将NNN系统分为三个分区，分别是东北电网、华北电网和华中电网。

步骤(1.2)设定仿真开始时刻时间t_s为0s，结束时刻t_e为5s和仿真步长h为0.01s。

步骤(1.3)设定同时协调的多时步数目m为20。

步骤(1.4)设置网络故障信息为华中分区“湘民丰500”母线三相接地，故障持续时间0s到0.09s。

步骤(2)设置当前的仿真时刻t_k(t₁＝0)，并生成对应的仿真时段信息[t_k，t_k+0.2]。

步骤(3)采用边界摄动法获得初始Jacobi逆矩阵。

步骤(3.1)初始电压U₀＝[0.00275 0.05733 0.98044 0.99247]^T，b＝F(U₀)＝[538145×10^-5 0.11037 2.45441×10^-5 0.13446]^T。

步骤(3.2)将U₀的每一维分别加上一个微小量e＝10⁴，得到4个向量：

$U_1=\left[\begin{array}{c}0.0028\\ 0.0573\\ 0.9804\\ 0.9925\end{array}\right],$ $U_2\left[\begin{array}{c}0.0027\\ 0.0574\\ 0.9804\\ 0.9925\end{array}\right],$ $U_3\left[\begin{array}{c}0.0027\\ 0.0573\\ 0.9805\\ 0.9925\end{array}\right],$ $U_4=\left[\begin{array}{c}0.0027\\ 0.0573\\ 0.9804\\ 0.9926\end{array}\right].$

将这些向量代入F(U)得到：

$b_1=\left[\begin{array}{c}3.6000\×{10}^{-4}\\ 1.1037\×{10}^{-1}\\ -3.5252\×{10}^{-6}\\ 1.3446\×{10}^{-1}\end{array}\right],$ $b_2=\left[\begin{array}{c}-1.5172\×{10}^{-5}\\ 1.1686\×{10}^{-1}\\ 1.5834\×{10}^{-5}\\ 1.3394\×{10}^{-1}\end{array}\right],$ $b_3=\left[\begin{array}{c}1.0247\×{10}^{-2}\\ 1.1066\×{10}^{-1}\\ 1.1179\×{10}^{-1}\\ 1.3328\×{10}^{-1}\end{array}\right],$ $b_4=\left[\begin{array}{c}2.0576\×{10}^{-4}\\ 1.1968\×{10}^{-1}\\ -1.1786\×{10}^{-3}\\ 2.4797\×{10}^{-1}\end{array}\right].$

步骤(3.3)应用式(5)、(6)可以解出jacobi逆矩阵：

$A=\left[\begin{array}{cccc}8.8833\×{10}^{-3}& 8.7517\×{10}^{-5}& 7.9404\×{10}^{-4}& 3.7203\×{10}^{-5}\\ 9.1151\×{10}^{-5}& 8.7829\×{10}^{-3}& 2.7257\×{10}^{-5}& 6.9681\×{10}^{-4}\\ -8.1070\×{10}^{-4}& -3.7751\×{10}^{-5}& 8.7010\×{10}^{-3}& 8.6528\×{10}^{-5}\\ -2.7953\×{10}^{-5}& -7.2014\×{10}^{-4}& 8.8933\×{10}^{-5}& 8.6875\×{10}^{-3}\end{array}\right].$

步骤(4)预估边界节点电压值。

采用式(7)所示的二次外推方法获得边界节点电压预估值。第一个时步无历史数据可用，因此将初值作为预估值。另外，故障期间系统波动剧烈，因此采用单时步积分，多时步积分在故障清除后启动。

步骤(5)协调服务向分区服务发送边界节点电压U₁＝U₀＝[0.00275 0.057330.98044 0.99247]^T，分区服务计算出界节点注入功率，并将之返还给协调服务。协调服务对节点注入功率求代数和，由式(4)得到偏差向量F(U₁)＝[538145×10^-5 0.11037 2.45441×10^-5 0.13446]^T；采用逆Broyden拟牛顿法计算边界节点电压的修正量ΔU＝[6.7605×10^-6 8.7260×10^-4 1.5223×10^-51.2450×10^-3]^T，并更新Jacobi逆矩阵：

$A=\left[\begin{array}{cccc}8.8836\×{10}^{-6}& 8.7583\×{10}^{-5}& 7.9405\×{10}^{-4}& 3.7251\×{10}^{-5}\\ 9.1113\×{10}^{-5}& 8.7945\×{10}^{-3}& 2.7433\×{10}^{-5}& 7.0616\×{10}^{-4}\\ -8.1072\×{10}^{-4}& -3.7600\×{10}^{-5}& 8.7011\×{10}^{-3}& 8.6637\×{10}^{-5}\\ -2.8001\×{10}^{-5}& -7.1340\×{10}^{-4}& 8.9147\×{10}^{-5}& 8.7087\×{10}^{-3}\end{array}\right].$

步骤(6)判断本时段仿真是否结束。

因为F(U₁)＞10^-4，本时段尚未收敛，因此返回步骤(5)，并令k＝k+1，直到本时段收敛；

步骤(7)判断全部仿真是否结束。

不断增加时段号，直到所有时段的仿真全部完成。

执行步骤(1)~步骤(8)即完成多时步数为20的一次暂态稳定仿真。为了比较多时步对收敛速度、仿真效率的影响，特将m由1到20逐步增加，其边界方程求解次数、高速局域网的加速比情况和40ms延时的广域互联网络的加速比情况见图3、图4和图5。其中加速比定义为待仿真的动态过程的时间(本算例为5s)与实际仿真耗时的比值。

Claims (1)

1.一种电力系统暂态稳定的分布式仿真方法，其特征在于，所述方法是在计算机中依次按以下步骤进行的：

步骤(1).初始化

步骤(1.1).采用节点撕裂方法进行电力网络的分区；

设定：S₁、S₂为两个区域电网，简称子分区，所提供的计算服务称为分区服务，

SB是由边界节点构成的上级电网，简称协调侧，所提供的计算服务简称协调服务，

S_1In、S_2In为所述的两个子分区的内部节点，

B为边界节点集合，B₁、B₂、B₃是同一边界节点在分区1、分区2和协调侧的表示， 

是子分区的边界节点在协调侧对应的虚拟节点，

i表示子分区序号，

所述的边界节点满足以下方程：

u、θ分别是所述边界节点的电压幅值和相角，P、Q分别是所述边界节点的注入有功和无功；

步骤(1.2).设定仿真始末时间和仿真步长

t_s表示仿真起始时刻，t_s＝0s，

t_e表示仿真结束时刻，t_e＝5s～80s，

h表示单时步的仿真步长，h＝0.01s～0.02s，在多时步情况下每个仿真时段包含多个h； 

步骤(1.3).设定同时协调的多时步数目

m表示同时协调的多时步数目，简称多时步数，采用多时步同时协调能够一次性完成m个时步的仿真；

步骤(1.4).根据仿真的需要选择网络故障，设定网络故障信息；

步骤(2).设定当前的仿真时刻，并生成对应的仿真时段信息

仿真时段信息表示每个时段内含有的m个时步，当前仿真的时段的序号用参数k表示，k＝1表示第一个时段，时段的始末时刻表示为[t_k，t_k+m·h]；

步骤(3).采用边界摄动法获得初始Jacobi逆矩阵，其步骤如下：

步骤(3.1).当前仿真时段k的边界协调方程为：

其中：U、ΔP、ΔQ分别表示边界节点的电压向量、注入有功功率差、注入无功功率差；

步骤(3.2).在电压向量的初值为U₀时，得到线性方程为：

JU₀＝b

其中U₀表示电压初值，J是U＝U₀处的Jacobi矩阵，b是F(U)的函数值；

步骤(3.3).通过电力系统网络向边界节点电压有变化的分区传送摄动信息，把电压向量U₀的每一维分别逐次加上一个微小量e＝10^-4，得到几个向量，表示如下：

[U₀₁+e，U₀₂，…，U_0n]^T，[U₀₁，U₀₂+e，…，U_0n]^T，…，[U₀₁，U₀₂，…，U_0n+e]^T，并代入F(U)，得到： 

n为电压向量U的维数；

步骤(3.4).将JU₀＝b与步骤(3.3)的方程逐列相减，得到：

e·J·I＝[b₁-b b₂-b...b_n-b]

I是与J同阶的单位矩阵，

从上式可解出Jacobi逆矩阵J^-1，用A表示；

步骤(4).预估边界节点电压值

当开始一个新时段的第一次仿真计算时，若：

存在三组以上的历史数据，指的是当前仿真时步之前的边界节点电压值，则按下式所示的二次外推以获得边界节点电压预估值：

U_c-1，U_c-1，U_c分别是第c-2，c-1，c时步的边界节点电压值，

t表示预估的时步序号；

步骤(5).采用逆Broyden拟牛顿法计算边界节点电压的修正量，并更新Jacobi逆矩阵

步骤(5.1).协调侧向各分区发送边界节点电压信息，并接收各分区返回的边界节点注入功率信息

协调服务向各分区服务发送边界节点电压U_k，各分区服务计算边界节点注入功率，并返回协调服务，所述协调服务按下式计算： 

得到：偏差向量F(U_k)＝[ΔP，ΔQ]^T；

步骤(5.2).按下式计算修正后的边界节点电压值U′_k：

U′_k＝U_k-A_kF(U_k)

步骤(5.3).协调侧向分区发送U′_k并获得分区反馈的F(U′_k)，从而得到s_k＝U′_k-U_k，y_k＝F(U′_k)-F(U_k)；

步骤(5.4).利用下式更新A‘_k，并令新的U_k＝U′_k，A_k＝A‘_k：

步骤(6).判断本时段仿真是否结束

判断是否满足收敛条件F(U_k)＜10^-4，如果收敛则转步骤(7)，如果不收敛则返回步骤(5)，并令k＝k+1；

步骤(7).判断所有时段仿真是否结束

判断是否满足仿真结束条件，如果t_i≥t_e则仿真结束，输出相关信息；反之，返回步骤(2)，更新时间信息t_i＝t_i+(m+1)×h，并令k＝1。 

Images