[发明专利]一种计及一次调频不确定性的概率潮流分析方法

摘要

本发明公开提出了一种计及一次调频不确定性的概率潮流分析方法，构造计及一次调频特性的动态潮流模型，应用点估计方法处理模型中的不确定问题，采用牛拉法求解确定性潮流问题，从而求出扰动不确定情况下，计及系统一次调频特性及其系数不确定性时的系统频率、节点电压等状态变量的统计特征值。在计及一次调频特性的动态潮流模型设计中，考虑发电机和负荷单位调节功率必然存在的不确定性，因此可以得到更接近实际情况的状态变量，为相关的分析计算提供更准确的数据。该方法可以更为有效准确地描述系统的运行状态，为电力系统进一步的分析计算提供基础。同时，点估计法具有速度快，低阶准确率高等优点，为本方法的实现提供了有效途径。

主权项

一种计及一次调频不确定性的概率潮流分析方法，其特征在于，包括以下步骤：1)输入确定参数及随机变量的分布信息；输入常规潮流计算所需要的基础信息、随机变量的分布信息和随机变量的数量；所述输入常规潮流计算所需要的基础信息，包括：网络拓扑数据，支路阻抗，母线并联电导电纳，节点类型，节点负荷功率，PV节点上发电机出力数据，PQ节点上发电机出力数据，发电机的额定容量及功率因素，发电机的调差系数标幺值Rg和系统负荷扰动功率预测值；所述随机变量的分布信息，随机变量包括：节点负荷扰动预测误差ΔPL_pre、负荷的单位调节功率kL和发电机的单位调节功率kG；所述节点负荷扰动预测误差ΔPL_pre服从正态分布，确定其正态分布的期望和方差；所述发电机的单位调节功率kG服从二项分布，确定发电机单位调节功率标幺值的两点取值KGg＝1/Rg和KGg＝0，及其对应的概率；所述负荷的单位调节功率kL服从均匀分布，确定其标幺值kLg的分布区间[a',b']；所述随机变量的数量为输入系统中随机变量的数量，总个数m，其中随机变量ΔPL_pre的个数为mp,随机变量kL的个数为mkL，随机变量kG的个数为mkG；2)计算基准状态下的潮流分布；结合步骤1获得的所述常规潮流基础信息，建立如下常规潮流模型；PGi-PLi0=UiΣj=1j=nUj(Gijcosδij+Bijsinδij)QGi-QLi0=UiΣj=1j=nUj(Gijsinδij-Bijcosδij)---(1)]]>公式1是基于极坐标的常规潮流方程，其中PGi,QGi分别表示节点i所接发电机的有功出力和无功出力；分别表示基准状态下，节点i上负荷吸收的有功功率和无功功率；Ui,Uj分别表示节点i的电压幅值和节点j的电压幅值，δij＝θi‑θj表示节点i、j间电压的相角差，其中θi和θj分别表示节点i和节点j电压的相角；Gij,Bij分别表示节点i、j间的电导和电纳，若节点i、j间无支路则对应的导纳为零；运用牛拉法求解上述方程，得到基准状态下发电机有功出力无功出力节点电压幅值和相角3)构建计及系统一次调频特性的动态潮流模型；结合步骤1获得的所述常规潮流基础信息以及步骤2求解得到的所述基准状态变量，构建负荷扰动时，计及负荷和发电机一次调频作用的动态潮流模型如下：PGi-PLi=UiΣj=1j=nUj(Gijcosδij+Bijsinδij)QGi0-QLi0=UiΣj=1j=nUj(Gijsinδij-Bijcosδij)PGi=PGi0-kGiΔfPLi=PLi′+kLiΔf=(PLi0+ΔPLi)+kLiΔf---(2)]]>公式2中，系统发生扰动后，计及一次调频作用时节点i上发电机的实际有功出力用PGi表示，负荷吸收的实际有功功率用PLi表示；kGi表示节点i上发电机的单位调节功率,kLi表示负荷的单位调节功率；Δf＝f‑f0为系统的频率偏移，其中f0和f分别表示基准状态的系统频率(50Hz)和扰动后的系统频率；P′Li表示扰动发生瞬间节点i上的有功，表示扰动前基准状态节点i上的有功负荷,ΔPLi表示基准状态节点i上发生的实际负荷扰动功率；4)根据分布信息，计算随机变量所需的统计特征值；对概率密度为fX(x)的变量X,求取下列统计特征值：期望μX：μX=∫-∞+∞xfX(x)dx---(3)]]>标准差σX：σX=(∫-∞+∞(x-μX)2fX(x)dx)---(4)]]>X的j阶中心距Mj(X)：Mj(X)=∫-∞+∞(x-μX)jfX(x)dx---(5)]]>定义λX_j表示随机变量X的j阶中心距与标准差σX的j次方的比值：λX_j=Mj(X)σXj---(6)]]>其中，λX_1＝0，λX_2＝1，λX_3和λX_4分别表示随机变量X的偏度系数和峰度系数；结合公式(3)‑(6)，并根据步骤1获得的所述随机变量分布信息，计算所述随机变量的期望、均方差、偏度系数和峰度系数；计算过程如下：I)负荷扰动预测误差；实际的扰动功率ΔPL由负荷扰动的预测值PL_pre及预测误差ΔPL_pre两部分组成，计算表达式如公式7：ΔPL＝PL_pre+ΔPL_pre     (7)用变量X1表示预测误差ΔPL_pre，各统计特征值求取如下：利用公式8计算扰动预测误差ΔPL_pre的偏度系数λX1_3：λX1_3=M3(X1)σX13=0σX13=0---(8)]]>利用公式9计算扰动预测误差ΔPL_pre的峰度系数λX1_4：λX1_4=M4(X1)σX14=3σX14σX14=3---(9)]]>II)发电机的单位调节功率；机组的单位调节功率kG服从二项分布，p表示参与一次调频的概率，对应的单位调节功率的值为1/R，其中R是机组的调差系数；q＝1‑p表示不参与一次调频的概率，对应的单位调节功率的值为0；将步骤1中输入的调差系数标幺值Rg转换为其有名值R，转换公式如下：R=RgfNPGN---(10)]]>其中fN＝50Hz，为系统额定频率；PGN是机组的额定有功功率，通过公式11求取：PGN＝SGNcosα     (11)其中SGN机组的额定容量，α为功率因素；用变量X2表示发电机的单位调节功率kG，其各统计特征值求取如下：利用公式12计算发电机的单位调节功率kG的期望μX2：μX2=p1R+q×0=pR---(12)]]>利用公式13计算发电机的单位调节效率kG的均方差σX2：σX2=p(1R-pR)2+q(0-pR)2=pqR---(13)]]>利用公式14计算发电机的单位调节功率kG的偏度系数λX2_3：λX2_3=M3(X2)σX23=p(1/R-p/R)3+q(0-p/R)3(pq/R)3=q2-p2pq---(14)]]>利用公式15计算发电机的单位调节功率kG的峰度系数λX2_4：λX2_4=M4(X2)σX24=p(1/R-p/R)4+q(0-p/R)4(pq/R)4=q3+p3pq---(15)]]>III)负荷的单位调节功率；负荷的单位调节功率标幺值kLg在[a',b']区间内均匀分布；与步骤II一样，需要先将标幺值区间[a',b']转换为有名值区间[a,b]，按照公式16转换：kL=kL·PLNfN---(16)]]>其中PLN表示额定有功；fN＝50Hz，为系统额定频率；用变量X3表示负荷的单位调节功率kL，其各统计特征值求取如下：利用公式17计算负荷的单位调节功率kL的期望μX3：μX3=b-a2---(17)]]>利用公式18计算负荷的单位调节功率kL的均方差σX3：σX3=∫ab(x-μX3)21b-adx=b-a23---(18)]]>利用公式19计算负荷的单位调节功率kL的偏度系数λX3_3：λX3_3=M3(X3)σX33=∫ab(x-μX3)31b-adxσX33=0---(19)]]>利用公式20计算负荷的单位调节功率kL的峰度系数λX3_4：λX3_4=M4(X3)σX34=∫ab(x-μX3)41b-adxσX34=95---(20)]]>5)根据2m+1点估计法构造随机变量的样本点，确定样本值及其权重系数；首先求取各样本值对应的权重系数；根据点估计原理，对系统中的所有随机变量zi(zi∈{z1,z2,…zm})分别按照公式21构造三个样本值zi,1,zi,2,zi,3；zi,j＝μi+ξi,jσi i＝1,…m；j＝1,2,3     (21)其中ξi,j表示第i个随机变量的第j个位置系数，通过公式22求取；ξi,1=λi,32+λi,4-34λi,32ξi,2=λi,32-λi,4-34λi,32ξi,3=0---(22)]]>其中λi3和λi4分别表示步骤4中求取的第i个随机变量的偏度系数和峰度系数；基于公式22得到的参数，利用公式23求取各样本值对应的权重系数wi,1、wi,2和wi,3；wi,1=1ξi12-ξi1ξi2wi,2=1ξi12-ξi1ξi2wi,3=1m-1λi,4-λi,32---(23)]]>其中m表示系统中随机变量的个数；然后，基于得到的负荷扰动预测误差的三个样本点，按照步骤4中的公式7，计算得到负荷扰动的样本；最后构造系统样本点；构造随机变量样本值zi,j(i＝1,…m；j＝1,2,3)对应的系统样本点时，随机变量zi＝zi,j，其余变量取均值，即样本点Z(i,j)＝(μ1,…,zi,j,…μm)；其中Z(1,3)＝Z(2,3)＝…＝Z(m,3)＝(μ1,…,μi,…μm)，只需进行一次计算，故对有m个随机变量的系统，只需完成2m+1次确定性潮流计算即可；6)参数初始化；完成步骤5后，确定潮流计算所需的初始参数；设置最大迭代次数Tmax，收敛精度ε＝10‑5～10‑3；根据步骤2得到的基准状态数据，确定发电机有功出力初值、无功出力初值分别为和节点i电压幅值和相角初值分别为和系统频率初值fN＝50Hz；利用公式24计算节点i的注入无功Qi_input0=QGi0-QLi0---(24)]]>7)采用牛拉法计算各样本对应的系统状态变量值；完成步骤6后，基于步骤5得到的由负荷扰动、发电机单位调节功率和负荷单位调节功率组成的系统样本和步骤3构造的动态潮流模型，分别计算3m个样本点对应的确定性潮流，获得该样本点对应的系统频率f，节点电压幅值U和相角θ，发电机的实际有功出力PG和负荷吸收的实际有功功率PL，计算过程包括；i)设置迭代次数初值，计算初始不平衡量；首先设置迭代次数初值t＝0；利用公式25计算基准状态发生扰动ΔPLi瞬间的初始不平衡量；ΔPi0=Pi_input0-Ui0Σj=1j=nUj0(Gijcosδij0+Bijsinδij0)(i∈[pq;pv;ref])ΔQi0=Ui_input0-Ui0Σj=1j=nUj0(Gijsinδij0-Bijcosδij0)(i∈pq)---(25)]]>公式25中为基准状态节点i和节点j电压的相角差，即；δij0=θi0-θj0---(26)]]>表示扰动后节点i的注入有功，根据公式27求取；Pi_input0=PGi0-P′Li0P′Li0=PLi0+ΔPLi=PLi0+PLi_pre+ΔPLi_pre---(27)]]>其中和分别表示基准状态下节点i上发电机的有功出力和有功负荷；表示扰动瞬间节点i的有功负荷；PLi_pre和ΔPLi_pre分别表示扰动预测值和扰动预测误差；此时系统的不平衡量为ΔF；ii)计算雅克比矩阵；基于步骤3提出的动态潮流模型，计算雅克比矩阵；J=HNMKLTWOV---(29)]]>其中H为npq+npv阶的方阵；Hij=∂ΔPi∂θj=UitUjt(Gijsinδijt-Bijcosδijt)(i≠j)-UitΣk=1k≠inUkt(Giksinδikt-Bikcosδikt)(i=j)(i∈[pq;pv];j∈[pq;pv])---(30)]]>N为(npq+npv)×npq阶的矩阵；Nij=∂ΔPi∂UjUj=UitUjt(Gijcosδijt+Bijsinδijt)(i≠j)UitΣk=1nUkt(Gikcosδikt+Biksinδikt)+(Uit)2Gii(i=j)(i∈[pq;pv];j∈pq)---(31)]]>M为(npq+npv)×1阶的矩阵；Mi=∂ΔPi∂f=(kGi+kLi)(i∈[pq:pv])---(32)]]>K为npq×(npq+npv)阶的矩阵；Kij=∂ΔQi∂θj=-UitUjt(Gijcosδijt+Bijsinδijt)(i≠j)UitΣk=1k≠inUkt(Gikcosδikt+Biksinδikt)(i=j)(i∈pq;j∈[pq;pv])---(33)]]>L为npq×npq阶的矩阵；Lij=∂ΔQi∂UjUj=UitUjt(Gijsinδijt-Bijcosδijt)(i≠j)UitΣk=1nUkt(Giksinδikt-Bikcosδikt)-(Uit)2Bii(i=j)(i∈pq;j∈pq)---(34)]]>T为npq×1阶的矩阵；Ti=∂Qi∂f=0(i∈pq)---(35)]]>W为1×(npq+npv)阶的矩阵；Wj=∂ΔPi∂θj=UitUjt(Gijsinδijt-Bijcosδijt)(i=ref,j∈[pq;pv])---(36)]]>O为1×npq阶的矩阵；Oj=∂ΔPi∂UjUj=UitUjt(Gijcosδijt+Bijsinδijt)(i=ref;j∈pq)---(37)]]>v为1×1阶的矩阵；V=∂ΔPreff=kGi+kLi(i=ref)---(38)]]>其pq,pv分别表示PQ节点序号矩阵，PV节点序号矩阵；ref为平衡节点的序号；n,npv,npq,nref(nref＝1)分别表示系统总的节点个数，PQ节点个数，PV节点个数和平衡节点个数；迭代次数自增，即t＝t+1；iii)解修正方程，计算修正量；基于步骤i和步骤ii，求解修正方程39，得到修正量；ΔX＝‑J‑1ΔF  (39)ΔX为变量的修正量，参见公式40；利用公式41对电压修正量进行变换；ΔUi=ΔUiUiUit-1(i∈pq)---(41)]]>其中Uit‑1表示上次迭代得到的节点i的电压幅值，初次迭代时指给定的电压幅值初值；进一步，结合公式40和41，计算得到的修正量修正原有的变量，得到新值；Xt＝Xt‑1+ΔX   (42)iv)计算不平衡量并进行收敛判断；按照公式25计算新值下的不平衡量ΔF；判断不平衡量ΔF是否满足收敛判据；max{|ΔF|}≤ε   (44)若收敛判据满足则迭代结束转到步骤v；若收敛判据不满足且迭代次数t＜Tmax则跳转至步骤ii，继续迭代直至满足收敛条件；若判据不满足且t＝Tmax，说明迭代不收敛，跳出循环，转至步骤v；v)结果输出；判断此时t是否满足t＜Tmax，若满足则此时得到的结果X即为样本Z(i,j)对应的状态变量值；否则表示潮流迭代不收敛；8)计算变量x的统计特征值；基于步骤7所得的3m个系统样本计算出的潮流结果以及步骤5计算得到的各样本点对应的权重系数，根据公式45‑47计算变量的期望和方差；E[xk]=Σi=1mΣj=13wi,j(X(i,j))k---(45)]]>其中X(1,3)＝X(2,3)＝…＝X(m,3)    (46)μx=E[x]σx=E[x2]-μx2---(47)]]>其中X(i,j)表示样本Z(i,j)计算得到的变量X的值；μx和σx分别表示变量X的期望和均方差；变量X包含节点电压相角，电压幅值和频率，