一种适合研究目的使用的直角坐标牛顿法潮流计算方法

1.本发明涉及一种电力系统牛顿法潮流计算方法，特别是一种适合研究目的使用的直角坐标牛顿法潮流计算方法。


背景技术：

2.电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据电力系统给定的运行条件和网络结构确定整个电力系统的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。由于具有收敛可靠、计算速度较快及内存需求适中的优点，牛顿法成为当前潮流计算的主流方法。牛顿法分为极坐标和直角坐标两种形式，两种形式的牛顿法潮流计算都在电力系统中得到了广泛的应用。
3.在直角坐标牛顿法潮流计算中，节点i的电压采用直角坐标表示为：
[0004][0005]
如图1所示，常规的直角坐标牛顿法潮流计算方法，主要包括以下步骤：
[0006]
a、输入原始数据和初始化电压。
[0007]
根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点注入有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为pq节点；节点注入有功功率和电压幅值已知、节点注入无功功率和电压相角未知的节点称为pv节点；节点电压幅值和电压相角已知、节点注入有功功率和无功功率未知的节点称为平衡节点。
[0008]
电压初始化采用平启动，即pv节点和平衡节点的电压实部取电压幅值的给定值，pq节点的电压实部取1.0；所有节点的电压虚部都取0.0。这里单位采用标幺值。
[0009]
b、形成节点导纳矩阵。
[0010]
c、设置迭代计数t＝0。
[0011]
d、计算节点功率及节点功率不平衡量和节点电压不平衡量，求最大不平衡量δw
max
。
[0012]
节点功率计算公式为：
[0013][0014]
式中，pi、qi分别为节点i的有功功率和无功功率；ei、fi分别为节点i的电压实部和虚部；ek和fk分别为节点k的电压实部和虚部；g
ik
、b
ik
分别为节点导纳矩阵元素y
ik
的实部和虚部；n为节点数。
[0015]
设第1～m号节点为pq节点，第m 1～n-1号节点为pv节点，第n号节点为平衡节点，节点功率不平衡量和节点电压不平衡量计算公式为：
[0016][0017]
式中，δpi、δqi分别为节点i的有功功率不平衡量和无功功率不平衡量；p
is
、q
is
分别为节点i给定的注入有功功率和注入无功功率；δu
i2
为节点i的电压不平衡量；u
is
为节点i给定的电压幅值；ui为节点i的电压幅值；m为pq节点数。
[0018]
平衡节点不参与迭代计算，不需要计算节点功率不平衡量和节点电压不平衡量。
[0019]
求各节点功率不平衡量和节点电压不平衡量中绝对值最大的值，称为最大不平衡量。
[0020]
e、判断最大不平衡量绝对值|δw
max
|是否小于收敛精度ε；如果小于收敛精度ε，执行步骤i；否则，执行步骤f。
[0021]
f、形成雅可比矩阵j。
[0022]
潮流计算的基本方程是非线性方程，通常采用逐次线性化方法迭代求解。
[0023]
线性化得到的修正方程(3)，用来求电压实部和虚部的修正量。
[0024][0025]
式中，j为(2n-2)
×
(2n-2)阶的雅可比矩阵，h、n、m、l、r、k分别为雅可比矩阵的六个分块子矩阵，其维数分别为(n-1)
×
(n-1)阶、(n-1)
×
(n-1)阶、m
×
(n-1)阶、m
×
(n-1)阶、(n-m-1)
×
(n-1)阶、(n-m-1)
×
(n-1)阶；δp＝[δp1,
…
,δp
n-1
]
t
为节点有功功率不平衡量列向量，上标t表示转置；δq＝[δq1,
…
,δqm]
t
为节点无功功率不平衡量列向量；为节点电压不平衡量列向量；δe＝[δe1,
…
,δe
n-1
]
t
为节点电压实部修正量列向量；δf＝[δf1,
…
,δf
n-1
]
t
为节点电压虚部修正量列向量。
[0026]
雅可比矩阵j的元素(i≠j时)计算公式如下：
[0027][0028][0029][0030]
[0031][0032][0033]
式中，h
ij
为节点i的有功功率不平衡量对节点j的电压实部的偏导数；n
ij
为节点i的有功功率不平衡量对节点j的电压虚部的偏导数；m
ij
为节点i的无功功率不平衡量对节点j的电压实部的偏导数；l
ij
为节点i的无功功率不平衡量对节点j的电压虚部的偏导数；r
ij
为节点i的电压不平衡量对节点j的电压实部的偏导数；k
ij
为节点i的电压不平衡量对节点j的电压虚部的偏导数。
[0034]
雅可比矩阵j的元素(i＝j时)计算公式如下：
[0035][0036][0037][0038][0039][0040][0041]
式中，h
ii
为节点i的有功功率不平衡量对节点i的电压实部的偏导数；n
ii
为节点i的有功功率不平衡量对节点i的电压虚部的偏导数；m
ii
为节点i的无功功率不平衡量对节点i的电压实部的偏导数；l
ii
为节点i的无功功率不平衡量对节点i的电压虚部的偏导数；r
ii
为节点i的电压不平衡量对节点i的电压实部的偏导数；k
ii
为节点i的电压不平衡量对节点i的电压虚部的偏导数。
[0042]
g、解修正方程及修正节点电压实部e和虚部f。
[0043]
解修正方程(3)，得到节点电压实部和虚部修正量。
[0044]
节点电压修正公式为：
[0045][0046]
式中，上标(t)表示第t次迭代；δei和δfi分别为节点i的电压实部修正量和电压虚部修正量。
[0047]
h、令t＝t 1，返回步骤d进行下一次迭代。
[0048]
i、输出节点及支路数据。
[0049]
式(3)所示的线性方程中，pq节点有相关的δp和δq方程，pv节点有相关的δp和δu2方程，而平衡节点相关的δp、δq和δu2方程都不存在。这样pq节点和pv节点对应两个方程，而平衡节点则一个都没有；而且平衡节点不需要求解δe和δf，其余每个节点则都需要求解δe和δf。这样节点和方程的对应关系需要根据节点类型确定，类型变化了，对应关系也要发生变化，给编程和调试带来了很大麻烦。虽然可以规定前m个节点为pq节点，第m 1～n-1号节点为pv节点，第n号节点为平衡节点。但这样规定，输入数据时有可能需要根据这一规定对节点重新编号；另外节点类型发生变化时节点和方程的对应关系还需要重新调整。
[0050]
为此，中国专利cn201010509556.5提出一种适合研究目的使用的牛顿法潮流计算方法，为以极坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的牛顿法潮流计算算法。该专利潮流计算的雅可比矩阵按(2n)
×
(2n)存储，pv节点δq对应的行和δu对应的列的元素都为0，与平衡节点有关的行和列的元素也都为0，这样虽然增加了内存需求量，但可以简化节点和方程系数矩阵行列的对应关系，大大降低编程难度，也不会增加计算量。求解潮流计算的修正方程时通过判断，跳过雅可比矩阵中主对角元素为0的行(主对角元素为0表示该行元素都是0，没有对应的方程)，不进行处理。
[0051]
中国专利cn201010509556.5的方法虽然简化了节点和方程系数矩阵行列的对应关系，大大降低了编程难度，但在求解潮流计算修正方程时，需要特殊处理，无法使用编程语言提供的性能优异的线性方程求解函数。


技术实现要素：

[0052]
为解决现有技术存在的上述问题，本发明要提出了一种适合研究目的使用的直角坐标牛顿法潮流计算方法，为以直角坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的直角坐标牛顿法潮流计算算法。
[0053]
为了实现上述目的，本发明对雅可比矩阵进行改进如下：雅可比矩阵仍按(2n)
×
(2n)存储，平衡节点有关的行和列的元素都为0，但对角元不清零，保留原来计算值。这样就可以使用编程语言提供的性能优异的线性方程求解函数求解潮流计算的修正方程了。
[0054]
本发明的技术方案如下：一种适合研究目的使用的直角坐标牛顿法潮流计算方法，包括以下步骤：
[0055]
a、输入原始数据和初始化电压。
[0056]
b、形成节点导纳矩阵。
[0057]
c、设置迭代计数t＝0。
[0058]
d、计算节点功率，计算各节点的不平衡量δw，求最大不平衡量δw
max
。
[0059]
按式(1)计算节点功率；
[0060]
不平衡量δw为2n维列向量，由有功功率不平衡量、无功功率不平衡量、电压不平衡量组成。
[0061]
如果节点i是pq节点，按式(17)计算不平衡量：
[0062][0063]
式中，δpi、δqi分别为节点i的有功功率不平衡量和无功功率不平衡量；p
is
、q
is
分别为节点i给定的注入有功功率和注入无功功率。
[0064]
如果节点i是pv节点，按式(18)计算不平衡量：
[0065][0066]
式中，δu
i2
为节点i的电压不平衡量；u
is
为节点i给定的电压幅值。
[0067]
如果节点i是平衡节点，按式(19)计算不平衡量：
[0068][0069]
求各节点有功功率不平衡量、无功功率不平衡量和电压不平衡量中绝对值最大的值，作为最大不平衡量δw
max
。
[0070]
e、判断最大不平衡量绝对值|δw
max
|是否小于收敛精度ε；如果小于收敛精度ε，执行步骤i；否则，执行步骤f。
[0071]
f、形成雅可比矩阵j。
[0072]
形成雅可比矩阵j的步骤如下：
[0073]
f1、不考虑节点类型，按式(20)-式(27)计算雅可比矩阵元素，形成(2n)
×
(2n)的雅可比矩阵j；
[0074]
雅可比矩阵j的元素(i≠j时)计算公式如下：
[0075][0076][0077][0078][0079]
式中，j
i,j
为雅可比矩阵第i行第j列的元素；为节点i的有功功率不平衡量对节点j的电压实部的偏导数；为节点i的有功功率不平衡量对节点j的电压虚部的偏导
数；为节点i的无功功率不平衡量对节点j的电压实部的偏导数；为节点i的无功功率不平衡量对节点j的电压虚部的偏导数。
[0080]
雅可比矩阵j的元素(i＝j时)计算公式如下：
[0081][0082][0083][0084][0085]
式中，为节点i的有功功率不平衡量对节点i的电压实部的偏导数；为节点i的有功功率不平衡量对节点i的电压虚部的偏导数；为节点i的无功功率不平衡量对节点i的电压实部的偏导数；为节点i的无功功率不平衡量对节点i的电压虚部的偏导数。
[0086]
f2、设置节点计数i＝1；
[0087]
f3、判断节点i是否为pv节点，如果不是pv节点，转到步骤f8；
[0088]
f4、设置计数j＝1；
[0089]
f5、按式(28)-式(31)计算雅可比矩阵元素j
n i,j
，j
n i,n j
；
[0090]
雅可比矩阵j的元素(i≠j时)计算公式如下：
[0091][0092][0093]
式中，为节点i的电压不平衡量对节点j的电压实部的偏导数；为节点i的电压不平衡量对节点j的电压虚部的偏导数。
[0094]
雅可比矩阵j的元素(i＝j时)计算公式如下：
[0095][0096][0097]
式中，为节点i的电压不平衡量对节点i的电压实部的偏导数；为节点i的电压不平衡量对节点i的电压虚部的偏导数。
[0098]
f6、令j＝j 1；
[0099]
f7、判断j是否大于n，如果j大于n，转到步骤f8；否则返回到步骤f5；
[0100]
f8、令i＝i 1；
[0101]
f9、判断i是否大于n，如果i大于n，转到步骤f10；否则返回到步骤f3；
[0102]
f10、设置节点计数i＝1；
[0103]
f11、判断节点i是否为平衡节点，如果不是平衡节点，转到步骤f21；
[0104]
f12、设置计数k＝1；
[0105]
f13、判断k＝i是否成立，如果成立，转到步骤f16；
[0106]
f14、令j
i,k
＝0；
[0107]
f15、令j
k,i
＝0；
[0108]
f16、判断k＝n i是否成立，如果成立，转到步骤f19；
[0109]
f17、令j
n i,k
＝0；
[0110]
f18、令j
k,n i
＝0；
[0111]
f19、令k＝k 1；
[0112]
f20、判断k是否大于2n，如果k大于2n，转到步骤f21；否则返回到步骤f13；
[0113]
f21、令i＝i 1；
[0114]
f22、判断i是否大于n，如果i大于n，转到步骤g；否则返回到步骤f11。g、解修正方程(32)及修正节点电压实部e和虚部f。
[0115][0116]
式中，j为(2n)
×
(2n)阶的雅可比矩阵；δe＝[δe1,
…
,δen]
t
为节点电压实部修正量列向量，上标t表示转置；δf＝[δf1,
…
,δfn]
t
为节点电压虚部修正量列向量；δw为2n维不平衡量列向量。
[0117]
节点电压修正公式为：
[0118][0119]
式中，上标(t)表示第t次迭代；δei和δfi分别为节点i的电压实部修正量和电压虚部修正量。
[0120]
h、令t＝t 1，返回步骤d进行下一次迭代。
[0121]
i、输出节点及支路数据。
[0122]
与现有技术相比，本发明具有以下有益效果：
[0123]
1、本发明潮流计算的雅可比矩阵按(2n)
×
(2n)存储，可以简化节点和方程系数矩阵行列的对应关系，大大降低了编程难度，也不会增加计算量。
[0124]
2、对于平衡节点不存在的方程，雅可比矩阵对应行中除对角元外的其他元素都为0，用来表示不存在对应的方程。由于对角元不为0，可以使用编程语言提供的线性方程求解函数求解潮流计算的修正方程，可以大大降低编程难度。编程语言提供的线性方程求解函数经过了优化，计算速度很快，稳定性也非常高，可以大大提高潮流计算的计算速度和程序稳定性。
[0125]
3、本发明潮流计算的节点编号灵活，编号与节点类型无关，不要求pq节点编号在前，平衡节点为最后一个节点，不需要根据程序设计要求改变系统节点的原有编号。
附图说明
[0126]
本发明共有附图2张。其中：
[0127]
图1是常规直角坐标牛顿法潮流计算的流程图。
[0128]
图2是本发明雅可比矩阵形成的流程框图。
具体实施方式
[0129]
下面结合附图对本发明进行进一步地说明。根据图1所示的直角坐标牛顿法潮流计算的流程图和图2所示的本发明雅可比矩阵形成的流程框图，对一个实际大型电网进行了潮流计算。该实际大型电网有445个节点，含有大量的小阻抗支路，为了使得常规潮流计算方法也能够计算，把这些小阻抗支路改为正常支路。潮流计算的收敛精度为0.00001。
[0130]
作为对比，采用以下3种方法对该实际大型电网进行了潮流计算：
[0131]
常规方法：常规的直角坐标牛顿法潮流计算方法；
[0132]
专利方法：专利cn201010509556.5的方法；
[0133]
对比方法：采用专利cn201010509556.5的方法，但其中潮流计算的修正方程求解采用matlab的线性方程求解函数。
[0134]
各种潮流计算方法都采用matlab的m文件编程，不同潮流计算方法的计算时间见表1。
[0135]
表1不同潮流计算方法的计算时间(s)
[0136]
方法常规方法专利方法对比方法本发明计算时间23.10330321.780145系数矩阵奇异，无法求解0.316549
[0137]
由表1可见，对于修改后的445节点实际电力系统算例，常规直角坐标牛顿法潮流计算方法和专利方法的计算时间接近，本发明方法由于解方程时采用matlab编程语言提供的线性方程求解函数计算，计算速度要快得多。对比方法的雅可比矩阵中的有些对角元为0，造成修正方程的系数矩阵奇异，无法用编程语言提供的线性方程求解函数求解。
[0138]
本发明可以采用任何一种编程语言和编程环境实现，如c语言、c 、fortran、delphi、matlab的m文件编程语言等。开发环境可以采用visual c 、borland c builder、visual fortran、matlab等。
[0139]
本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。