发明名称 一种电力系统状态估计方法
摘要 本发明属于电力系统运行与控制技术领域,尤其涉及一种基于插值矩阵的含相角量测装置(Phasor Measurement Unit-PMU)的电力系统状态估计方法。本发明将SCADA量测更新的间隔期间通过PMU插值得到的区域U状态等效为区域U的先验状态信息,并在每次PMU采样时刻都对此先验状态信息进行不断的更新,直到下一时刻的SCADA和PMU量测同时到达时,插值矩阵更新,从而进一步计算出区域U各个节点的实时运行状态,这样就充分兼顾了PMU量测精度高、采样速度快的优点,能够更为准确和实时的追踪系统的运行状态。
申请公布号 CN103972884B 申请公布日期 2016.03.02
申请号 CN201410166675.3 申请日期 2014.04.24
申请人 西南交通大学 发明人 张葛祥;赵俊博
分类号 H02J3/00(2006.01)I;G06F19/00(2011.01)I 主分类号 H02J3/00(2006.01)I
代理机构 成都宏顺专利代理事务所(普通合伙) 51227 代理人 李顺德;王睿
主权项 一种电力系统状态估计方法,其特征在于,包括以下步骤:S1、在电网中配置PMU,同时对电网进行划分,具体如下:S11、在发电机和一级、二级负荷节点处安装PMU;S12、根据S11所述PMU的配置情况将电网划分为PMU可观测区域和PMU不可直接观测但SCADA可观测区域,其中,PMU可观测区域记作区域O,PMU不可直接观测但SCADA可观测区域记作区域U,N<sub>O</sub>为区域O的节点数,N<sub>U</sub>为区域U的节点数,电网中总节点数为N=N<sub>O</sub>+N<sub>U</sub>;S2、读取电网的数据信息,所述数据信息包括:区域O的网络参数、拓扑结构和线路阻抗,区域U的网络参数、拓扑结构和线路阻抗;S3、根据S2所述电网中的数据信息并结合形成电力系统中相应的节点导纳矩阵和支路‑节点关联矩阵;S4、读取电网的SCADA量测配置信息,所述SCADA量测配置信息包括:节点电压幅值量测、节点电流幅值量测、节点功率注入量测和线路潮流量测;S5、计算插值矩阵H<sup>0</sup>,所述插值矩阵H<sup>0</sup>由S3所得的节点导纳矩阵经插值处理得到,具体插值处理方法如下:S51、电力系统通用节点电流方程为[I<sub>bus</sub>]=[Y<sub>bus</sub>]·[E<sub>bus</sub>],其中,I<sub>bus</sub>是节点注入电流向量,Y<sub>bus</sub>是节点导纳矩阵,E<sub>bus</sub>是节点电压向量;S52、根据S12所述电网划分将节点导纳矩阵分为区域O节点自导纳矩阵Y<sub>OO</sub>,区域O节点互导纳矩阵Y<sub>OU</sub>,区域U的自导纳矩阵Y<sub>UU</sub>,连接区域O和区域U的互导纳矩阵Y<sub>UO</sub>;S53、得到节点电流方程<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>I</mi><mi>O</mi></msub></mtd></mtr><mtr><mtd><msub><mi>I</mi><mi>U</mi></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>Y</mi><mrow><mi>O</mi><mi>O</mi></mrow></msub></mtd><mtd><msub><mi>Y</mi><mrow><mi>O</mi><mi>U</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mi>Y</mi><mrow><mi>U</mi><mi>O</mi></mrow></msub></mtd><mtd><msub><mi>Y</mi><mrow><mi>U</mi><mi>U</mi></mrow></msub></mtd></mtr></mtable></mfenced><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>E</mi><mi>O</mi></msub></mtd></mtr><mtr><mtd><msub><mi>E</mi><mi>U</mi></msub></mtd></mtr></mtable></mfenced><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000011.GIF" wi="544" he="169" /></maths>其中,I<sub>O</sub>是区域O中节点注入电流,E<sub>O</sub>是区域O中的节点电压,I<sub>U</sub>是区域U中的节点注入电流,E<sub>U</sub>是区域U中的节点电压;S54、将I<sub>U</sub>等效转化为一个负荷导纳向量Y<sub>U</sub>,即,当区域U中的节点数目为N<sub>U</sub>,那么区域U中各个节点处的电流注入可以表示为<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mo>&lsqb;</mo><msub><mi>I</mi><mi>U</mi></msub><mo>&rsqb;</mo><mo>=</mo><mo>&lsqb;</mo><msup><mrow><mo>(</mo><mfrac><msub><mi>S</mi><mi>i</mi></msub><msub><mi>E</mi><mi>i</mi></msub></mfrac><mo>)</mo></mrow><mo>*</mo></msup><mo>&rsqb;</mo><mo>=</mo><mo>&lsqb;</mo><mrow><mo>(</mo><msub><mi>P</mi><mi>i</mi></msub><mo>-</mo><msub><mi>jQ</mi><mi>i</mi></msub><mo>)</mo></mrow><mo>/</mo><msubsup><mi>E</mi><mi>i</mi><mo>*</mo></msubsup><mo>&rsqb;</mo><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>...</mo><msub><mi>N</mi><mi>U</mi></msub><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000012.GIF" wi="1026" he="188" /></maths>其中,S<sub>i</sub>为节点i处的复功率,P<sub>i</sub>为节点i处的有功功率,Q<sub>i</sub>为节点i处的无功功率,[Y<sub>U</sub>]=[Y<sub>L</sub>]=[I/E]=[(P<sub>i</sub>‑jQ<sub>i</sub>)/|E|<sub>i</sub><sup>2</sup>],i=1,2,...,N<sub>U</sub>,Y<sub>L</sub>为一个N<sub>U</sub>×N<sub>U</sub>的对角矩阵;S55、将S54所得电流注入表达式[I<sub>U</sub>]代入S53所述节点电流方程,得到<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>I</mi><mi>O</mi></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>Y</mi><mrow><mi>O</mi><mi>O</mi></mrow></msub></mtd><mtd><msub><mi>Y</mi><mrow><mi>O</mi><mi>U</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mi>Y</mi><mrow><mi>U</mi><mi>O</mi></mrow></msub></mtd><mtd><msub><mi>Y</mi><mi>T</mi></msub></mtd></mtr></mtable></mfenced><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>E</mi><mi>O</mi></msub></mtd></mtr><mtr><mtd><msub><mi>E</mi><mi>U</mi></msub></mtd></mtr></mtable></mfenced><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000021.GIF" wi="534" he="161" /></maths>然后计算得到Y<sub>UO</sub>·E<sub>O</sub>+Y<sub>T</sub>·E<sub>U</sub>=0,其中,Y<sub>T</sub>=Y<sub>UU</sub>+Y<sub>L</sub>;S56、得到系统的插值矩阵H<sup>0</sup>=‑Y<sub>T</sub><sup>‑1</sup>·Y<sub>UO</sub>,其中,矩阵Y<sub>T</sub><sup>‑1</sup>·Y<sub>UO</sub>是稀疏的,所述矩阵Y<sub>T</sub><sup>‑1</sup>·Y<sub>UO</sub>维数是N<sub>U</sub>×N<sub>O</sub>;S6、根据区域O中基于PMU量测插值得到的先验信息(x<sub>k</sub>)并利用改进加权最小二乘法估计出区域U中各节点的节点状态向量<img file="FDA00008324368000000211.GIF" wi="84" he="85" />并由此计算出各节点处的注入复功率,所述区域U的节点处的注入复功率计算方法如下:S61、电力系统在第k次采样时的状态x<sub>k</sub>可以用z<sub>k</sub>=h(x<sub>k</sub>)+v<sub>k</sub>表示,其中,x是状态向量,z<sub>k</sub>为量测向量,h(·)是m维非线性函数向量,v<sub>k</sub>是服从正态分布的随机白噪声,即v<sub>k</sub>~N(0,R<sub>k</sub>),R<sub>k</sub>是量测误差的方差;S62、考虑到PMU采样速率通常比SCADA快很多,在SCADA量测更新的间隔期间通过PMU插值得到的区域U状态可以等效为区域U的先验信息;S63、对于区域U,以节点电压为状态变量为例,根据S62所述区域U的先验信息有如下的优化目标函数:<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><munder><mrow><mi>m</mi><mi>i</mi><mi>n</mi></mrow><msub><mi>x</mi><mi>k</mi></msub></munder><mi>J</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msup><mrow><mo>&lsqb;</mo><msub><mi>z</mi><mi>k</mi></msub><mo>-</mo><mi>h</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>)</mo></mrow><mo>&rsqb;</mo></mrow><mi>T</mi></msup><msub><mi>W</mi><mi>k</mi></msub><mo>&lsqb;</mo><msub><mi>z</mi><mi>k</mi></msub><mo>-</mo><mi>h</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>)</mo></mrow><mo>&rsqb;</mo><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>)</mo></mrow><msubsup><mi>P</mi><mi>k</mi><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>)</mo></mrow><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000022.GIF" wi="1455" he="131" /></maths>其中,<img file="FDA0000832436800000023.GIF" wi="63" he="75" />是k时刻的先验状态,在两次SCADA量测采样间隔期间,<img file="FDA0000832436800000024.GIF" wi="59" he="77" />可由式<img file="FDA0000832436800000025.GIF" wi="699" he="99" />得到,P<sub>k</sub>是N<sub>U</sub>×N<sub>U</sub>阶的误差协方差矩阵,对上面的目标函数进行优化求解可得:<maths num="0005" id="cmaths0005"><math><![CDATA[<mrow><mo>&lsqb;</mo><msup><mi>Q</mi><mi>T</mi></msup><msub><mi>W</mi><mi>k</mi></msub><mi>Q</mi><mo>+</mo><msubsup><mi>P</mi><mi>k</mi><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><mo>&rsqb;</mo><mi>&Delta;</mi><mi>x</mi><mo>=</mo><msup><mi>Q</mi><mi>T</mi></msup><msub><mi>W</mi><mi>k</mi></msub><mi>&Delta;</mi><mi>z</mi><mo>+</mo><msubsup><mi>P</mi><mi>k</mi><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><mi>&Delta;</mi><mover><mi>x</mi><mo>&OverBar;</mo></mover><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000026.GIF" wi="804" he="95" /></maths>其中<maths num="0006" id="cmaths0006"><math><![CDATA[<mrow><mi>&Delta;</mi><mover><mi>x</mi><mo>&OverBar;</mo></mover><mo>=</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>-</mo><msub><mi>x</mi><mi>k</mi></msub><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000027.GIF" wi="304" he="83" /></maths>W<sub>k</sub>是对角权重矩阵,<maths num="0007" id="cmaths0007"><math><![CDATA[<mrow><msub><mi>W</mi><mi>k</mi></msub><mo>=</mo><msubsup><mi>R</mi><mi>k</mi><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><mo>=</mo><msup><mrow><mo>(</mo><mi>d</mi><mi>i</mi><mi>a</mi><mi>g</mi><mo>&lsqb;</mo><msubsup><mi>&sigma;</mi><mrow><mi>k</mi><mo>,</mo><mn>1</mn></mrow><mn>2</mn></msubsup><mo>,</mo><msubsup><mi>&sigma;</mi><mrow><mi>k</mi><mo>,</mo><mn>2</mn></mrow><mn>2</mn></msubsup><mo>...</mo><msubsup><mi>&sigma;</mi><mrow><mi>k</mi><mo>,</mo><mi>N</mi></mrow><mn>2</mn></msubsup><mo>&rsqb;</mo><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000028.GIF" wi="769" he="114" /></maths><img file="FDA0000832436800000029.GIF" wi="297" he="84" />表示在k时刻各个量测噪声方差,Q是h(x)关于x的雅克比矩阵,所述权重W<sub>k</sub>是由新的自适应量测权重分配函数来计算,<img file="FDA00008324368000000210.GIF" wi="551" he="131" />其中,λ是一个大于0的自然数;S64、得到电力系统状态向量的迭代等式x<sub>k+1</sub>=x<sub>k</sub>+Δx;S65、当Δx<10<sup>‑6</sup>,则算法收敛,输出当前电力系统各个节点的电压和相角值,进而估计出区域U中各节点处的注入复功率;S7、判断电力系统是否出现较大扰动,如果<img file="FDA0000832436800000031.GIF" wi="466" he="117" />则进入S10,如果<img file="FDA0000832436800000032.GIF" wi="464" he="112" />则进入S8;S8、更新插值矩阵,插值矩阵的更新式子为H<sup>0</sup>+ΔH,其中,ΔH为电力系统在运行点发生变化造成的插值矩阵H<sup>0</sup>的偏移量,所述插值矩阵更新方法如下:S81、在给定的参考运行点‘0’处有:<img file="FDA0000832436800000033.GIF" wi="415" he="102" />其中,所述给定的参考运行点‘0’为稳态时电力系统的初始运行状态;S82、电力系统的运行点可以由<maths num="0008" id="cmaths0008"><math><![CDATA[<mrow><mi>H</mi><mo>=</mo><mo>-</mo><msup><mrow><mo>(</mo><msub><mi>Y</mi><mrow><mi>U</mi><mi>U</mi></mrow></msub><mo>+</mo><mi>d</mi><mi>i</mi><mi>a</mi><mi>g</mi><mo>&lsqb;</mo><msup><mrow><mo>(</mo><msub><mi>S</mi><mi>i</mi></msub><mo>)</mo></mrow><mo>*</mo></msup><mo>/</mo><mo>|</mo><msub><mi>E</mi><mi>i</mi></msub><msup><mo>|</mo><mn>2</mn></msup><mo>&rsqb;</mo><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>&CenterDot;</mo><msub><mi>Y</mi><mrow><mi>U</mi><mi>O</mi></mrow></msub><mo>=</mo><mo>-</mo><msubsup><mi>Y</mi><mi>T</mi><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><mo>&CenterDot;</mo><msub><mi>Y</mi><mrow><mi>U</mi><mi>O</mi></mrow></msub></mrow>]]></math><img file="FDA0000832436800000034.GIF" wi="1018" he="119" /></maths>确定,即有Y<sub>T</sub>H=‑Y<sub>UO</sub>,对此等式两边求导可得ΔY<sub>T</sub>·H+Y<sub>T</sub>·ΔH=0,从而可以得到ΔH=‑Y<sub>T</sub><sup>‑1</sup>·ΔY<sub>T</sub>·H由于Y<sub>T</sub>=Y<sub>UU</sub>+diag[(S<sub>i</sub>)<sup>*</sup>/|E<sub>i</sub>|<sup>2</sup>],那么对Y<sub>T</sub>求导可得:<img file="FDA00008324368000000312.GIF" wi="555" he="122" />Y<sub>UU</sub>是常数矩阵,求导后为0;S83、如果初始运行点相对应的复功率为<img file="FDA00008324368000000313.GIF" wi="93" he="80" />定义复功率偏移量为:<img file="FDA0000832436800000035.GIF" wi="362" he="104" />,其中,<img file="FDA0000832436800000036.GIF" wi="70" he="80" />为当前时刻的注入复功率;那么插值矩阵的偏移量可以表示为<maths num="0009" id="cmaths0009"><math><![CDATA[<mrow><mi>&Delta;</mi><mi>H</mi><mo>=</mo><mo>-</mo><msup><mrow><mo>(</mo><msubsup><mi>Y</mi><mi>T</mi><mn>0</mn></msubsup><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>&CenterDot;</mo><mo>&lsqb;</mo><mi>d</mi><mi>i</mi><mi>a</mi><mi>g</mi><mo>&lsqb;</mo><msubsup><mi>&Delta;S</mi><mi>U</mi><mo>*</mo></msubsup><mo>/</mo><mo>|</mo><msub><mi>E</mi><mi>i</mi></msub><msup><mo>|</mo><mn>2</mn></msup><mo>&rsqb;</mo><mo>&rsqb;</mo><mo>&CenterDot;</mo><mi>H</mi><mo>,</mo></mrow>]]></math><img file="FDA0000832436800000037.GIF" wi="810" he="119" /></maths>从而插值矩阵的更新式子为H<sup>0</sup>+ΔH;S84、将S6中已经估计出的U区域各个节点的注入复功率值和电压值带入插值矩阵的更新式子中就得到了更新后的插值矩阵;S9、通过插值矩阵和式子<maths num="0010" id="cmaths0010"><math><![CDATA[<mrow><msub><mi>E</mi><mi>U</mi></msub><mo>=</mo><mrow><mo>(</mo><msup><mi>H</mi><mn>0</mn></msup><mo>+</mo><mi>&Delta;</mi><mi>H</mi><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>E</mi><mi>O</mi></msub><mo>=</mo><mrow><mo>(</mo><msup><mi>H</mi><mn>0</mn></msup><mo>-</mo><msup><mrow><mo>(</mo><msubsup><mi>Y</mi><mi>T</mi><mn>0</mn></msubsup><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>&CenterDot;</mo><mo>&lsqb;</mo><mi>d</mi><mi>i</mi><mi>a</mi><mi>g</mi><mo>&lsqb;</mo><msubsup><mi>S</mi><mi>U</mi><mo>*</mo></msubsup><mo>/</mo><mo>|</mo><msub><mi>E</mi><mi>i</mi></msub><msup><mo>|</mo><mn>2</mn></msup><mo>-</mo><msup><mrow><mo>(</mo><msubsup><mi>S</mi><mi>U</mi><mn>0</mn></msubsup><mo>)</mo></mrow><mo>*</mo></msup><mo>/</mo><mo>|</mo><msubsup><mi>E</mi><mi>U</mi><mn>0</mn></msubsup><msup><mo>|</mo><mn>2</mn></msup><mo>&rsqb;</mo><mo>&rsqb;</mo><mo>&CenterDot;</mo><mi>H</mi><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>E</mi><mi>O</mi></msub></mrow>]]></math><img file="FDA0000832436800000038.GIF" wi="1580" he="140" /></maths>对区域U各节点状态进行计算,得到当前时刻区域U各个节点的电压,其中,E<sub>U</sub>为区域U中的节点电压,E<sub>O</sub>为区域O中的节点电压,<img file="FDA0000832436800000039.GIF" wi="70" he="82" />和<img file="FDA00008324368000000310.GIF" wi="75" he="79" />分别为初始运行点和当前时刻的注入复功率,<img file="FDA00008324368000000311.GIF" wi="76" he="78" />和E<sub>i</sub>分别为初始运行点和当前时刻的电压;S10、判断收敛性,当目前电力系统量测采样次数k<TH<sub>2</sub>时,算法收敛,则输出区域U中各个节点的电压,当k≥TH<sub>2</sub>时,则进入S5,其中,TH<sub>2</sub>为预先设定的系统量测总采样次数。
地址 610031 四川省成都市二环路北一段111号