基于对称直角坐标的快速极坐标牛顿-拉夫逊潮流方法与流程

本发明属于电力系统分析计算领域，具体涉及一种基于对称直角坐标的快速极坐标牛顿-拉夫逊潮流方法。

背景技术：

牛顿-拉夫逊法(牛顿法)是电力系统最常用的潮流计算方法，根据坐标体系和修正方程的不同，可分为直角坐标和极坐标牛顿法。两者计算原理类似，迭代次数和潮流计算速度可能有所不同。设系统节点数为n、pq节点数为m，直角坐标牛顿法与极坐标牛顿法的主要问题和特点如下：

1.修正方程式个数不同。极坐标牛顿法个数为(n-1+m)，直角坐标牛顿法修正方程式个数为2(n-1)，极坐标牛顿法比直角坐标牛顿法少pv节点所对应的δvi²的方程个数(n-1-m)，还少求对应的(n-1-m)个δvi²的值。一般当pv节点数较多时，极坐标牛顿法具有较大的速度优势；当pv节点数较少时，直角坐标牛顿法具有一定速度优势。

2.雅可比矩阵j中元素的分类和数量不同。极坐标牛顿法j阵中的元素分h、n、m、l四类，直角坐标牛顿法j阵中元素分h、n、m、l、r、s六类。如果计算rij、sij元素，直角坐标牛顿法j阵元素的计算个数为4(n-1)²，极坐标牛顿法j阵元素的计算个数为(n-1+m)²；如果不计算rij、sij元素，直角坐标牛顿法j阵元素计算量仍比极坐标牛顿法多2(n-1-m)个。

3.潮流计算的迭代次数可能不同。一般极坐标牛顿法迭代次数与直角坐标牛顿法迭代次数相同，但在有些系统潮流计算中，极坐标牛顿法比直角坐标牛顿法少1次。

4.极坐标和直角坐标牛顿法的j阵元素不对称，因此j阵元素的计算效率较低。

5.j阵元素和节点电流ipi、iqi或节点功率δpi、δqi的计算方式不同。极坐标牛顿法中均含大量三角函数计算，对计算速度影响较大，而直角坐标牛顿法无此问题。

6.潮流计算过程中角度与弧度的转换。极坐标牛顿法求解修正方程式得到弧度，有时要进行弧度与角度的转换，而直角坐标牛顿法无此问题。

上述第1～3条使得极坐标牛顿法比直角坐标牛顿法更有速度优势；第4条对两者均有影响；第5条又使得极坐标牛顿法的速度不如直角坐标牛顿法；第6条的影响可忽略不计。因此，如果对极坐标牛顿法实现j阵元素的对称计算并减少三角函数计算，则可使极坐标牛顿法的计算速度大大超过直角坐标牛顿法。

纯数学计算中，由于计算过程严格按照给定的数学公式进行计算，所以很难实现计算过程的拆分计算或综合计算，也就很难提高计算速度。因此在极坐标牛顿法中对j阵元素计算中仅考虑同一子阵jij中元素hij、nij、mij、lij的关系，而不考虑不同子阵jij和jji中元素的关系，所以很难实现j阵元素的对称计算。且由于极坐标牛顿法的修正方程式和j阵元素数量均远少于直角标牛顿法，因此减少三角函数计算的问题几乎被忽略不计了。

电力系统工程计算中，可根据电力系统的特点和实际工程计算过程需要对计算过程拆分计算或综合计算，以提高计算速度。如果将j阵中各个元素的对应关系扩展到各个子阵的对应关系，则可拆分j阵元素计算，实现j阵元素的部分对称计算；如果对j阵元素等计算部分进行三角变换，则可大大减少三角函数计算；如果将取倒后的对角元素作为规格化因子，则可大大除法计算；如果用四角规则对修正方程式进行消元计算，则可将高斯消元法计算原理简单化，简化计算过程、便于编程。

技术实现要素：

为解决上述第4～5条的问题，提高极坐标牛顿法的计算效率和计算速度，本发明提出一种基于对称直角坐标的快速极坐标牛顿-拉夫逊潮流方法。

将极坐标和直角坐标形式的节点电压分别表示如下：

其之间关系为：ei＝vicosδi，fi＝visinδi

本发明是通过以下技术方案实现的，主要步骤如下：

步骤1：打开数据文件，读取y阵及系统数据到相关数组；

步骤2：三角变换极坐标牛顿法潮流中δpi、δqi计算式并计算；

(1)极坐标牛顿法潮流中δpi、δqi计算式如下；

(2)三角变换后δpi、δqi计算式如下；

步骤3：拆分极坐标牛顿法j阵元素的计算过程，建立j阵元素的部分对称关系并对其进行三角变换，按“二行+二列”方式计算j阵元素；

设系统节点数为n，pq节点数为m，m+1及其后节点均为pv节点，第n节点为平衡节点，极坐标牛顿法修正方程式中j阵元素排列如下：

(1)拆分j阵元素的计算过程，建立j阵元素的部分对称关系以简化元素计算；

1)拆分计算非对角子阵jij和jji中非对角元素(j≠i)，建立部分对称关系。

传统方法形成j阵时最多仅利用一个jij子阵中元素(hij、nij、mij、lij)关系，如hij＝lij和nij＝-mij，仅略微简化元素计算。但由于hij≠hji、nij≠nji、mij≠mji、lij≠lji，因此无法利用jij子阵中元素和对应子阵jji中元素(hji、nji、mji、lji)的关系，从而必须分别求取jij和jji子阵中元素，形成j阵速度极不理想。但如果对对应子阵jij和jji中元素进行拆分计算，则能建立jij和jji子阵中元素的对应关系，加快j阵形成速度。各jij和jji子阵元素拆分计算如下。

上述有三角函数的均是传统计算式，无法显示元素的部分对称性，形成j阵时只能分别计算jij和jji子阵中的元素。用变量c、d表示拆分后的元素计算式，不但可分别显示jij和jji子阵中hij和hji、nij和nji、mij和mji、lij和lji元素的部分对称性，还能显示jij和jji子阵中不同元素之间的关系如mij＝-nij和mji＝-nji、lij＝hij和lji＝hji。因此，计算出c1、d1，就可直接得到hij、hji、lij、lji四个元素；而计算出c2、d2，就可直接得到nij、nji、mij、mji四个元素。

上述计算表明，尽管数学上对应子阵jij和jji中元素hij≠hji、nij≠nji、mij≠mji、lij≠lji关系不对称，但在拆分计算后对应子阵jij和jji元素有部分对称关系。因此利用jij和jji子阵中元素的部分对称关系，可大大简化j阵下三角元素计算，提高j阵形成速度。

2)拆分计算对角子阵jii中对角元素(j＝i)，建立部分对称关系。

传统方法中对角子阵jii中元素关系为hii≠lii和nii≠-mii，因此必须分别计算hii、nii、mii、lii，从而影响j阵形成速度。如果对对角子阵jii中元素进行拆分计算，建立类似非对角子阵jij中元素对应关系，同样可加快j阵形成速度。各元素拆分计算如下。

上述式(1)是传统计算式，由于生成j阵之前需要求取功率不平衡量δpi、δqi因而已经求解了pi、qi，若用式(1)计算各个元素则含对pi、qi的重复计算。式(2)的简化计算式则可免去对pi、qi的重复计算。但由于hii≠lii和nii≠-mii，因此式(1)、式(2)均无法利用元素的部分对称性，须分别计算jii子阵中元素。式(3)是是拆分后的计算式，其中h′ii、n′ii、m′ii、l′ii是按非对角元方式计算的对角元素，有l′ii＝h′ii，m′ii＝-n′ii，即在不考虑pi、qj修正的情况下，jii子阵中元素也具有部分对称关系。

因此在计算非对角元素时同时利用同组元素关系计算相应的对角元素部分，所有元素计算完成后再用pi、qi修正对角元素，同样可加快j阵形成速度。

(2)三角变换j阵元素计算式，以减少三角函数计算；

1)非对角元素(j≠i)的三角变换：

简化上述计算式后有：

hij＝c1+d1、hji＝-c1+d1；nij＝c2-d2、nji＝c2+d2；

mij＝-c2+d2＝-nij、mji＝-c2-d2＝-nji；lij＝c1+d1＝hij、lji＝-c1+d1＝hji

其中各c、d具体公式如下：

c1＝gij(eifj-fiej)，d1＝bij(eiej+fifj)；c2＝-gij(eiej+fifj)，d2＝bij(fiej-eifj)

2)对角元素(j＝i)的三角变换：

由于对对角元素可先按非对角元方式计算出h′ii、n′ii、m′ii、l′ii，再按前述三角变换的pi、qi修正，因此只需对h′ii、n′ii、m′ii、l′ii进行三角如下变换：

h′ii＝gii(eifi-fiei)+bii(eiei+fifi)＝c1+d1

n′ii＝-gii(eiei+fifi)-bii(fiei-eifi)＝c2-d2

m′ii＝gii(eiei+fifi)+bii(fiei-eifi)＝-c2+d2＝-n′ii

l′ii＝gii(eifi-fiei)+bii(eiei+fifi)＝c1+d1＝h′ij

上述三角变换可大大省略极坐标j阵元素中三角函数计算，且利用了j阵对角元素之间的对应关系，因此可大大提高极坐标牛顿法中j阵的形成速度。

(3)按“二行+二列”方式计算j阵元素以提高j阵形成速度；

传统方法中一般是按每次一行或每次二行方式计算j阵元素，未应用y阵元素和j阵元素的对称性。本发明中由于按对称法计算j阵元素，因此只需要读入y阵对角元素和上三角元素，大大减少了读入y阵数据时间；同时可按“二行+二列”方式对称形成j阵以提高j阵形成速度。

步骤4：将取倒的对角元素作为规格化因子，用四角规则完成修正方程式的消元计算并回代求取δvi、δδi；

(1)对j阵元素规格化前先将对角元素取倒作为规格化因子，大大减少除法计算；

传统高斯消元法是将未取倒的对角元素作为规格化因子，因子表法是形成因子表后再将对角元素取倒，程序中有大量的除法计算。取倒后的对角元素作为规格化因子可使计算速度提高约1％。

(2)用四角规则完成对修正方程式的消元计算而无需消元计算公式；

设对第k列的元素规格化消元前后的简化矩阵及各元素定义如下：

对角元素(参考元素，规格化因子)；

规格化前交叉元素(规格化后^(k))，对角元素同行以右；

消元元素，对角元素同列以下；

计算元素前值(新值^(k))，与消元元素同行、交叉元素同列交互点上。

将对角元素取倒作为规格化因子，乘以以右所有交叉元素得完成规格化，再对以下所有消元元素消元得其以右各个计算元素前值新值

可将高斯消元计算公式用文字表述为：对角元素为参考元素，计算元素新值＝其前值-消元元素*规格化后交叉元素。据此文字表述可直接写出各个计算元素的新值如下：

由于参与计算的四个相关元素均在矩形的四个角上，故称四角规则。因此根据四角规则，即根据元素在矩阵中位置可直接写出相应消元计算式而无需消元计算公式，更无需考虑元素上下标。四角规则特别便于对消元计算过程的理解和编程。

(3)求出电压幅值和相角的新值vi^(k+1)＝vi^(k)+δvi^(k)、δi^(k+1)＝δi^(k)+δδi^(k)，再通过三角变换得到电压的实部和虚部ei^(k+1)、fi^(k+1)，再计算δpi、δqi。

步骤5：判断是否满足收敛条件；

如果δpi、δqi不满足收敛条件，则跳转到步骤3，继续形成新的j阵，再进行消元和回代计算；如果满足收敛条件，则执行步骤6。

步骤6：结束迭代并输出结果。

本发明的有益效果为：本方法根据电力系统潮流计算的特性和工程计算的特点拆分极坐标牛顿法j阵元素的计算过程，建立j阵元素的部分对称关系以大大减少j阵元素的计算；对j阵元素和δpi、δqi的计算式进行三角变换以大大减少三角函数的计算；按“二行+二列”方式计算j阵元素以快速取用y阵数据和快速生成j阵；将取倒后的对角元素作为规格化因子以大大减少规格化中除法计算；用四角规则完成对修正方程式的消元和回代计算以简化计算过程、便于编程。本发明方法大大提高了j阵的形成速度和极坐标牛顿法的潮流计算速度，简化了极坐标牛顿法的计算过程。

附图说明

图1为极坐标牛顿法潮流计算框图。

图2为本发明方法潮流计算流程图。

具体实施方式

本发明将通过以下实施例作进一步说明。

实施例1。分别用传统极坐标牛顿法(方法1)、j阵元素部分对称的极坐标牛顿法(发明方法1)、直角坐标解法的极坐标牛顿法(发明方法2)和本发明方法(发明方法3)计算ieee-14、-30、-57、-118节点系统潮流，其潮流计算的平均时间、生成整个j阵的平均时间和迭代次数等计算结果及比较如表1所示。

以ieee-118节点系统为例，根据表1计算结果可以看出：

(1)j阵元素部分对称的极坐标牛顿法、直角坐标解法的极坐标牛顿法和本发明方法生成整个j阵的平均时间分别为传统极坐标牛顿法的56.19％、18.85％、10.66％，说明采用发明方法1、发明方法2均会大大提高j阵形成的速度，而综合发明方法3更具优势，可提高约90％，且随着系统规模的加大而增加。

(2)传统极坐标牛顿法、j阵元素部分对称的极坐标牛顿法、直角坐标解法的极坐标牛顿法和本发明方法的步骤1～3生成整个j阵的平均时间分别为其潮流计算平均时间的14.35％、8.61％、3.90％、2.21％，同样说明发明方法1、发明方法2、发明方法3对提高j阵形成速度的作用，也同样随着系统规模的加大而增加。

(3)j阵元素部分对称的极坐标牛顿法、直角坐标解法的极坐标牛顿法和本发明方法潮流计算的平均时间分别为传统极坐标牛顿法的93.71％、69.34％、69.29％，说明采用发明方法1、发明方法2、发明方法3均会有效提高潮流计算速度。但由于潮流计算过程中，生成j阵时间远远小于求解线性方程时间，因此本发明方法提高潮流计算速度效果远不如提高j阵形成速度效果，但仍可将潮流计算速度提高约30％。

(4)本发明方法中还采用了将取倒的对角元素作为规格化因子的方法，该方法对生成j阵速度没有影响，而对潮流计算速度可提高约1％，但由于影响不大而省去了该计算结果的比较。此外，本发明方法中应用的四角规则只是将潮流计算中修正方程式的计算过程和计数原理简单化，特别适于编程计算。

表1各种方法对ieee系统潮流计算时间、生成j阵时间和迭代次数比较

tpf：潮流计算的平均时间；t∑j：生成整个j阵的平均时间；ins：潮流计算迭代次数。

本发明方法根据电力系统潮流计算的特性和工程计算的特点拆分极坐标牛顿法j阵元素的计算过程，建立j阵元素的部分对称关系以大大减少j阵元素的计算；对j阵元素和δpi、δqi的计算式进行三角变换以大大减少三角函数的计算；按“二行+二列”方式计算j阵元素以快速取用y阵数据和快速生成j阵；将取倒后的对角元素作为规格化因子以大大减少规格化中除法计算；用四角规则完成对修正方程式的消元和回代计算以简化计算过程、便于编程。本发明方法大大提高了j阵的形成速度和极坐标牛顿法的潮流计算速度，简化了极坐标牛顿法的计算过程。

本发明方法可以采用任何一种编程语言和编程环境实现，这里采用c++编程语言，开发环境是visualc++。