本发明涉及一种电力系统牛顿法潮流计算方法，特别是一种适合研究目的使用的极坐标牛顿法潮流计算方法。
背景技术：
：电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据给定的运行条件和网络结构确定整个网络的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。极坐标牛顿法潮流计算方法是一种最常用的潮流计算方法，科研人员经常以极坐标牛顿法潮流计算为基础进行进一步地研究。实用的商业软件采用稀疏矩阵技术和节点优化编号等高级技术。这些技术虽然能大幅度提高潮流计算的速度、降低内存占用量，但编程非常麻烦且难以修改和维护，不易增加新的功能，因而不适合科研人员用于研究目的使用。Matlab软件以矩阵为最基本的数据单位，可以方便地处理各种矩阵和向量运算，也可以很方便自然地处理复数类型，其指令表达式与数学中常用的形式很接近，还有大量常见实用的函数，给编程带来很大便利。Matlab软件简单易用、代码短小易操作，易于编程和调试，计算功能强大，同时还具有非常强大的可视化图形处理和交互式功能，为科学研究以及工程应用提供了一种高效的编程工具，目前已经成为许多科学领域的基本工具和首选平台，在各种科学和工程计算领域得到了广泛的应用。为了适应越来越多的科研人员需要在Matlab平台上以极坐标牛顿法潮流计算为基础进行进一步地研究的需求，迫切需要一种基于Matlab软件的易于编程、修改和调试的极坐标牛顿法潮流计算方法。根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。牛顿法潮流计算分为两类：牛顿法潮流计算中节点电压采用极坐标表示时的计算方法，称为极坐标牛顿法潮流计算方法；牛顿法潮流计算中节点电压采用直角坐标表示时的计算方法，称为直角坐标牛顿法潮流计算方法。极坐标牛顿法潮流计算主要方程如下：节点导纳矩阵为：式中，Yik为节点导纳矩阵元素，当下标i≠k时，为节点i和节点k的互导纳，当下标i＝k时，为节点i的自导纳；n为节点数。节点功率方程为：Pi=UiΣk=1nUk(Gikcosθik+Biksinθik)i=1,...,nQi=UiΣk=1nUk(Giksinθik-Bikcosθik)i=1,...,n---(2)]]>式中，Pi、Qi分别为节点i的节点有功功率和无功功率；Ui、Uk分别为节点i和节点k的节点电压幅值；θik＝θi-θk，θi和θk分别为节点i和节点k的节点电压相角；Gik、Bik分别为节点导纳矩阵元素Yik的实部和虚部。功率偏差方程为：ΔPi=PiS-Pi=PiS-UiΣk=1nUk(Gikcosθik+Biksinθik)i=1,...,n-1ΔQi=QiS-Qi=QiS-UiΣk=1nUk(Giksinθik-Bikcosθik)i=1,...,m---(3)]]>式中，ΔPi、ΔQi分别为节点i的节点有功功率偏差和无功功率偏差；PiS、QiS分别为节点i给定的节点注入有功功率和注入无功功率；m为PQ节点数。修正方程组为：ΔPΔQ=JΔθΔU/U=HNMLΔθΔU/U---(4)]]>式中，J为雅可比矩阵，H、N、M、L为雅可比矩阵的分块子矩阵。雅可比矩阵各元素计算公式为：当j≠i时Hij=∂ΔPi∂θj=-UiUj(Gijsinθij-Bijcosθij)---(5)]]>Nij=∂ΔPi∂UjUj=-UiUj(Gijcosθij+Bijsinθij)---(6)]]>Mij=∂ΔQi∂θj=UiUj(Gijcosθij+Bijsinθij)---(7)]]>Lij=∂ΔQij∂UjUj=-UiUj(Gijsinθij-Bijcosθij)---(8)]]>当j＝i时Hii=∂ΔPi∂θi=UiΣk=1k≠inUk(Giksinθik-Bikcosθik)---(9)]]>Nii=∂ΔPi∂UiUi=-UiΣk=1k≠inUk(Gikcosθik+Biksinθik)-2Ui2Gii---(10)]]>Mii=∂ΔQi∂θi=-UiΣk=1k≠inUk(Gikcosθik+Biksinθik)---(11)]]>Lii=∂ΔQi∂UiUi=-UiΣk=1k≠inUk(Giksinθik-Bikcosθik)+2Ui2Bii---(12)]]>或用下列公式计算：Hii=Ui2Bii+Qi---(13)]]>Nii=-Ui2Gii-Pi---(14)]]>Mii=Ui2Gii-Pi---(15)]]>Lii=Ui2Bii-Qi---(16)]]>式中，Pi、Qi分别为节点i的有功功率和无功功率，按式(2)计算。如图1-2所示，现有极坐标牛顿法潮流计算方法，主要包括以下步骤：A、原始数据输入和电压初始化；原始数据包括线路和变压器支路数据、节点注入有功功率和无功功率、节点电压幅值、节点无功补偿数据，以及收敛精度、最大迭代次数。电压初始化采用平启动，即PV节点和平衡节点的节点电压幅值取给定值，PQ节点的节点电压幅值取1.0；所有节点电压的相角都取0.0。这里单位采用标幺值。B、形成节点导纳矩阵；根据输入的线路和变压器支路数据形成如式(1)所示的节点导纳矩阵。C、形成雅可比矩阵；按式(5)～式(16)计算雅可比矩阵的各元素。D、计算节点功率及功率偏差；按式(2)计算节点功率，按式(3)计算节点功率偏差。E、解方程及修正节点电压幅值U和相角θ；解修正方程组(4)，求出电压幅值修正量列向量ΔU及电压相角修正量列向量Δθ。电压修正公式为：Ui(t+1)=Ui(t)-ΔUi(t),i=1,...,m---(17)]]>θi(t+1)=θi(t)-Δθi(t),i=1,...,n-1---(18)]]>式中，上标(t)表示第t次迭代的值；ΔUi和Δθi分别为节点i的节点电压幅值修正量和节点电压相角修正量。F、判断功率最大不平衡量|ΔP|max和|ΔQ|max是否都小于收敛精度ε；如果都小于收敛精度ε，进行步骤G，否则返回步骤C进行下一次迭代；G、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。直接采用上述原理实现的极坐标牛顿法潮流计算软件计算速度较慢，商业使用的极坐标牛顿法潮流计算软件采用稀疏矩阵技术和节点优化编号技术，比较复杂，不适合科研人员以此为基础进一步进行科学研究。因此，中国专利ZL201010509556.5提出了一种适合研究目的使用的牛顿法潮流计算方法，为以极坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的牛顿法潮流计算方法，其特点如下：(1)不采用稀疏矩阵技术和节点优化编号，大大降低了算法的实现难度；(2)通过简单逻辑判断来避免不必要运算，提高了潮流计算的计算速度。中国专利ZL201010509556.5所提出方法为以极坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供了一个易于修改和维护的极坐标牛顿法潮流计算方法。该方法采用C语言等编译型编程语言实现时速度很快，但使用Matlab这类解释型编程语言实现时计算速度则很慢，同时该专利也没有充分利用Matlab擅长矩阵运算和复数运算的特点。因此需要一个充分利用Matlab的特点且计算快速的潮流计算方法供在Matlab平台上进行科学研究的科研人员使用。技术实现要素：为解决现有技术存在的上述问题，本发明要提出一种基于Matlab的极坐标牛顿法潮流计算方法，可以充分利用Matlab特有的擅长矩阵运算和复数运算的特点，同时又有较快计算速度的潮流计算方法。为了实现上述目的，本发明的技术方案如下：一种基于Matlab的极坐标牛顿法潮流计算方法，采用矩阵运算和复数运算。包括以下步骤：A、原始数据输入和电压初始化；B、形成节点导纳矩阵；C、形成雅可比矩阵及计算节点功率；Matlab擅长矩阵运算和复数运算，因此采用Matlab编程，应该推导出基于矩阵运算和复数运算的雅可比矩阵计算方法。雅可比矩阵元素与节点类型有关，常规形成雅可比矩阵时要判断节点类型，根据节点类型决定哪些节点需要形成雅可比矩阵元素。这样处理对于通过循环来实现的算法，容易处理，但不适合通过矩阵整体运算形成雅可比矩阵的方法。因此，本发明形成雅可比矩阵时，不判断节点类型，对所有节点都形成雅可比矩阵元素，之后再去掉不需要的行和列。下面推导采用矩阵运算计算雅可比矩阵元素和节点功率的公式。对式(5)～式(8)的分析，可得Mij＝-Nij(19)Lij＝Hij(20)因此，先求Hij和Nij，求出Hij和Nij后，自然就能得到Mij和Lij。推导由矩阵运算计算雅可比矩阵的公式之前，先看看雅可比矩阵各元素如何用复数或相量表示。从式(5)和式(6)可见，两式中都包含UiUj和θij项，说明雅可比矩阵各元素中应该包含乘以的共轭的项，即包含项；观察式(5)和式(6)括号内的表达式，可知雅可比矩阵各元素中还应该包含Yij的共轭项。因此Hij、Nij、Mij、Lij可以由下式生成：Jij0=U·iU^jY^ij---(21)]]>式中，上标(^)表示复数的共轭；为节点电压列向量。式(21)可以看成由下式与节点导纳矩阵对应位置元素相乘得到的矩阵J0的第i行第j列的元素：式(22)则可由节点电压列向量与其共轭转置相乘得到，因此雅可比矩阵分块子矩阵可由下式得到：J0=(U·U·H).*Y^---(23)]]>式中，J0为雅可比初始计算矩阵；上标H表示矩阵的共轭转置；.*表示两矩阵对应行列的元素相乘。由J0得到初始雅可比矩阵分块子矩阵为：H0＝-Im(J0)(24)N0＝-Re(J0)(25)M0＝Re(J0)(26)L0＝-Im(J0)(27)式中，H0、N0、M0、L0为初始雅可比矩阵的分块子矩阵；Re表示取矩阵元素的实部；Im表示取矩阵元素的虚部。由式(24)～式(27)得到的初始雅可比矩阵分块子矩阵的非对角元已经是雅可比矩阵元素了，对角元还需要修正。由式(13)～式(16)，并考虑到sinθii＝0和cosθii＝1，得：Hii＝-UiUi(Giisinθii-Biicosθii)+Qi(28)Nii＝-UiUi(Giicosθii+Biisinθii)-Pi(29)Mii＝UiUi(Giicosθii+Biisinθii)-Pi(30)Lii＝-UiUi(Giisinθii-Biicosθii)-Qi(31)式(28)～式(31)中等式右侧的第1项就是H0、N0、M0、L0的对角元，因此只需对得到的H0、N0、M0、L0用Pi和Qi修正，就能得到雅可比矩阵分块子矩阵对角元。观察式(2)和式(6)，以及N0，可得Pi=-Σk=1nNik0---(32)]]>式中，为矩阵N0第i行第k列的元素。观察式(2)和式(5)，以及H0，可得Qi=-Σk=1nHik0---(33)]]>由式(24)、式(25)、式(32)、式(33)得到节点复功率为：S~i=Pi+jQi=Σk=1nJik0---(34)]]>式中，为节点i的复功率。由各节点复功率组成的节点复功率列向量可以用Matlab的一个矩阵求和函数实现：S~=sum(J0,2)---(35)]]>式中，为节点复功率列向量；sum为Matlab的矩阵求和函数；2表示对矩阵每一行的元素求和。用节点复功率对雅可比矩阵分块子矩阵对角元进行修正如下：Hii=Hii0+Im(S~i),i=1,...,n---(36)]]>Nii=Nii0-Re(S~i),i=1,...,n---(37)]]>Mii=Mii0-Re(S~i),i=1,...,n---(38)]]>Lii=Lii0-Im(S~i),i=1,...,n---(39)]]>形成雅可比矩阵及计算节点功率，包括以下步骤：C1、计算雅可比初始计算矩阵J0；C2、计算节点复功率；C3、由J0计算初始的雅可比矩阵分块子矩阵；C4、用节点复功率对雅可比矩阵分块子矩阵对角元进行修正；C5、由雅可比矩阵分块子矩阵形成雅可比矩阵；J=HNML---(40)]]>C6、对雅可比矩阵进行调整，去掉PV节点无功功率偏差对应的行以及平衡节点有功功率偏差和无功功率偏差对应的行；去掉PV节点电压幅值修正量对应的列以及平衡节点电压幅值修正量和电压相角修正量对应的列，结束。D、计算节点功率偏差；式(3)计算节点功率偏差公式写成矩阵运算的形成为：ΔP=PS-Re(S~)ΔQ=QS-Im(S~)---(41)]]>式中，ΔP、ΔQ分别为节点有功功率偏差列向量和无功功率偏差列向量；PS、QS分别为节点给定的注入有功功率列向量和注入无功功率列向量。计算得到的节点功率偏差列向量ΔP和ΔQ中去掉PV节点无功功率偏差及平衡节点有功功率偏差和无功功率偏差。E、解方程及修正电压幅值U和相角θ；直接调用Matlab软件的解线性方程组算法解修正方程组(4)，求出电压幅值修正量列向量ΔU及电压相角修正量列向量Δθ。对电压进行修正的式(17)和式(18)改写成矩阵形式为：U(t+1)＝U(t)-ΔU(t)(42)θ(t+1)＝θ(t)-Δθ(t)(43)式中，上标(t)表示第t次迭代的值；ΔU和Δθ分别为节点电压幅值修正量列向量和节点电压相角修正量列向量。F、判断功率最大不平衡量|ΔP|max和|ΔQ|max是否都小于收敛精度ε；如果都小于收敛精度ε，进行步骤G，否则返回步骤C进行下一次迭代。G、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。与现有技术相比，本发明具有以下有益效果：1、本发明提出的方法在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。2、本发明提出的方法采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能。3、由于Matlab对矩阵运算实行优化，采用矩阵运算要比按矩阵元素循环运算编程快得多，同时直接调用Matlab的方程求解算法，也大大提高了计算速度。实践证明，本发明的方法既方便了科研人员对程序进行编写、修改和调试，同时计算速度也基本接近了在C语言平台上实现的速度，为科研人员的科研工作提供了一个优秀的分析工具。附图说明本发明共有附图4张。其中：图1是现有极坐标牛顿法潮流计算的流程图。图2是现有极坐标牛顿法形成雅可比矩阵的流程图。图3是本发明极坐标牛顿法潮流计算的流程图。图4是本发明形成雅可比矩阵和计算节点功率的流程图。具体实施方式下面结合附图对本发明进行进一步地说明，按照图3-4所示流程对一个修改后的445节点实际系统算例进行了计算。445节点实际大型电力系统有445个节点，544条支路，含有大量的小阻抗支路。为了对各种方法进行比较，把这些小阻抗支路改为正常阻抗支路以满足各种方法的要求。采用本发明和几种对比方法对445节点实际系统算例进行了计算，计算时收敛精度为0.00001。几种潮流计算算法分别为：方法1：中国专利ZL201010509556.5方法，采用Matlab语言实现。方法2：中国专利ZL201010509556.5方法，采用Matlab语言实现，但解方程直接调用Matlab的解方程算法。方法3：中国专利ZL201010509556.5方法，采用Matlab语言实现，但解方程直接调用Matlab的解方程算法，同时对矩阵变量预定义维数。方法4：本发明方法。几种方法的计算时间见表1，计算时间不包括数据读入和输出及支路功率计算的时间。表1几种极坐标牛顿法潮流计算计算时间比较潮流计算方法计算时间(s)方法118.303方法23.095方法30.927方法40.519从表1可见，直接采用Matlab实现中国专利ZL201010509556.5方法，计算时间很长；采用Matlab实现专利ZL201010509556.5方法时，如果直接调用Matlab的解方程方法，计算速度就能大大提高，Matlab的解方程方法比较成熟稳定，有利于算法稳定性，Matlab的解方程方法的调用也非常简单，简化了编程，使程序更加清晰；如果同时对矩阵变量预定义维数，避免程序在执行过程中不断扩充矩阵的大小，可以大大提高计算速度。本发明的计算结果表明形成雅可比矩阵时采用矩阵运算技术可以进一步较大幅度地提高计算速度，程序也进一步简化。本发明可以在任何版本的MATLAB编程语言实现，但建议使用较新版本的MATLAB语言。本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。当前第1页1 2 3