CN107181260A - 基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 - Google Patents
基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 Download PDFInfo
- Publication number
- CN107181260A CN107181260A CN201710557642.5A CN201710557642A CN107181260A CN 107181260 A CN107181260 A CN 107181260A CN 201710557642 A CN201710557642 A CN 201710557642A CN 107181260 A CN107181260 A CN 107181260A
- Authority
- CN
- China
- Prior art keywords
- mrow
- node
- mover
- matrix
- mtd
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
- H02J3/04—Circuit arrangements for ac mains or ac distribution networks for connecting networks of the same frequency but supplied from different sources
- H02J3/06—Controlling transfer of power between connected networks; Controlling sharing of load between connected networks
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
- H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
Landscapes
- Engineering & Computer Science (AREA)
- Power Engineering (AREA)
- Supply And Distribution Of Alternating Current (AREA)
Abstract
本发明公开了一种基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法，采用矩阵运算和复数运算并使用Matlab提供的稀疏矩阵技术。本发明在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。本发明主要部分都采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能。本发明采用Matlab的稀疏矩阵技术并使用矩阵运算，大幅度地提高了计算速度，同时Matlab的稀疏矩阵使用非常方便，可以像全矩阵一样用行列号直接使用稀疏矩阵的元素，为科研人员的科研工作提供了一个优秀的分析工具。
Description
技术领域
本发明涉及一种电力系统牛顿法潮流计算方法，特别是一种适合研究目的使用的直角坐标牛顿法潮流计算方法。
背景技术
电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据给定的运行条件和网络结构确定整个网络的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。直角坐标牛顿法潮流计算方法是一种最常用的潮流计算方法，科研人员经常以直角坐标牛顿法潮流计算为基础进行进一步地研究。实用的商业软件采用C语言等高级编程语言编写且使用稀疏矩阵技术和节点优化编号等高级技术。这些技术虽然能大幅度提高潮流计算的速度、降低内存占用量，但编程非常麻烦且难以修改和维护，不易增加新的功能，因而不适合科研人员用于研究目的使用。
Matlab软件以矩阵为最基本的数据单位，可以方便地处理各种矩阵和向量运算，也可以很方便自然地处理复数类型，其指令表达式与数学中常用的形式很接近，还有大量常见实用的函数，给编程带来很大便利。Matlab软件简单易用、代码短小易操作，易于编程和调试，计算功能强大，同时还具有非常强大的可视化图形处理和交互式功能，为科学研究以及工程应用提供了一种高效的编程工具，目前已经成为许多科学领域的基本工具和首选平台，在各种科学和工程计算领域得到了广泛的应用。为了适应越来越多的科研人员需要在Matlab平台上以直角坐标牛顿法潮流计算为基础进行进一步地研究的需求，迫切需要一种基于Matlab软件的易于编程、修改和调试的直角坐标牛顿法潮流计算方法。
根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。
牛顿法潮流计算分为两类：牛顿法潮流计算中节点电压采用极坐标表示时的计算方法，称为极坐标牛顿法潮流计算方法；牛顿法潮流计算中节点电压采用直角坐标表示时的计算方法，称为直角坐标牛顿法潮流计算方法。直角坐标牛顿法潮流计算主要方程如下：
节点导纳矩阵为：
式中，Yik为节点导纳矩阵元素，当下标i≠k时，为节点i和节点k的互导纳，当下标i＝k时，为节点i的自导纳；n为节点数。
节点功率方程为：
式中，Pi、Qi分别为节点i的节点有功功率和无功功率；ei、ek分别为节点i和节点k的节点电压相量实部；fi和fk分别为节点i和节点k的节点电压相量虚部；Gik、Bik分别为节点导纳矩阵元素Yik的实部和虚部。
功率不平衡量及电压平方不平衡量方程为：
式中，ΔPi、ΔQi分别为节点i的节点有功功率不平衡量和无功功率不平衡量；ΔUi 2为节点i的节点电压平方不平衡量；Pis、Qis分别为节点i给定的节点注入有功功率和注入无功功率；Uis为节点i给定的节点电压幅值；m为PQ节点数。
功率不平衡量及电压平方不平衡量方程也可以表示为：
式中，ai、bi分别为节点i的节点电流相量的实部和虚部，为：
修正方程组为：
式中，J为雅可比矩阵，H、N、M、L、R、K为雅可比矩阵的分块子矩阵。
雅可比矩阵各元素计算公式为：
当j≠i时
Hij＝-Gijei-Bijfi (7)
Nij＝Bijei-Gijfi (8)
Mij＝Bijei-Gijfi (9)
Lij＝Gijei+Bijfi (10)
Rij＝0 (11)
Kij＝0 (12)
当j＝i时
Hii＝-Giiei-Biifi-ai (13)
Nii＝Biiei-Giifi-bi (14)
Mii＝Biiei-Giifi+bi (15)
Lii＝Giiei+Biifi-ai (16)
Rii＝-2ei (17)
Kii＝-2fi (18)
如图1-2所示，现有直角坐标牛顿法潮流计算方法，主要包括以下步骤：
A、原始数据输入和电压初始化；
原始数据包括线路和变压器支路数据、节点注入有功功率和无功功率、节点电压幅值、节点无功补偿数据，以及收敛精度、最大迭代次数。
电压初始化采用平启动，即PV节点和平衡节点的节点电压实部取给定值，PQ节点的节点电压实部取1.0；所有节点电压的虚部都取0.0。这里单位采用标幺值。
B、形成节点导纳矩阵；
根据输入的线路和变压器支路数据形成如式(1)所示的节点导纳矩阵。
C、形成雅可比矩阵；
按式(5)、式(7)-(18)计算雅可比矩阵的各元素。
D、计算节点功率及功率不平衡量和电压平方不平衡量；
按式(2)计算节点功率，按式(3)计算节点功率不平衡量和节点电压平方不平衡量。
E、解方程及修正节点电压实部e和虚部f；
解修正方程组(6)，求出电压实部修正量列向量Δe及电压虚部修正量列向量Δf。
电压修正公式为：
fi (t+1)＝fi (t)-Δfi (t)i＝1，…，n-1 (20)
式中，上标(t)表示第t次迭代的值；Δei和Δfi分别为节点i的节点电压实部修正量和节点电压虚部修正量。
F、判断最大不平衡量|ΔP|max、|ΔQ|max和|ΔU2|max是否都小于收敛精度ε；如果都小于收敛精度ε，转步骤G，否则返回步骤C进行下一次迭代；
G、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。
直接采用上述原理实现的直角坐标牛顿法潮流计算软件计算速度较慢，商业使用的直角坐标牛顿法潮流计算软件采用稀疏矩阵技术和节点优化编号技术，比较复杂，不适合科研人员以此为基础进一步进行科学研究。因此，中国专利ZL201610864281.4提出了一种基于Matlab的直角坐标牛顿法潮流计算方法，可以充分利用Matlab特有的擅长矩阵运算和复数运算的特点，设计出了简洁又有较快计算速度的潮流计算方法，为以直角坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的牛顿法潮流计算方法，其特点如下：
(1)在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析；
(2)采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能；
(3)采用矩阵运算并直接调用Matlab的方程求解算法，大大提高了计算速度。
中国专利ZL201610864281.4提出的一种基于Matlab的直角坐标牛顿法潮流计算方法，为从事电力系统研究的科研人员提供了一种基于Matlab平台的易于修改和维护且计算较为快速的直角坐标牛顿法潮流计算方法。该方法采用Matlab实现，充分利用Matlab擅长矩阵运算和复数运算的特点，但未使用稀疏矩阵技术，计算速度相对较慢，仍有待进一步提高计算速度。
发明内容
为解决现有技术存在的上述问题，本发明要提出一种基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法，充分利用Matlab特有的擅长矩阵运算和复数运算的特点，并且运用Matlab提供的稀疏矩阵技术，实现提高潮流计算的计算速度的目的。
为了实现上述目的，本发明的技术方案如下：基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法，采用矩阵运算和复数运算并使用Matlab提供的稀疏矩阵技术，包括以下步骤：
A、原始数据输入和电压初始化；
B、记录相关节点类型的节点号；
直角坐标牛顿法潮流计算的修正方程组方程个数及变量个数与电力系统的节点类型有关，ΔP方程组中没有平衡节点有功功率不平衡量对应的方程，ΔQ方程组中仅有PQ节点无功功率不平衡量对应的方程，ΔU2方程组中仅有PV节点电压平方不平衡量对应的方程；变量则不包含平衡节点的相角变量和电压幅值变量。
为了提高计算速度，形成雅可比矩阵及方程右端向量时先不考虑节点类型，形成雅可比矩阵及方程右端向量后解修正方程时，再去掉无关的行和列。为此，设置3个数组记录有关节点类型的节点号，其中数组bt1记录PV节点的节点号，数组bt2记录PQ节点和PV节点的节点号，数组bt记录雅可比矩阵及方程右端向量需要的行列号。
记录相关节点类型的节点号使用Matlab的find函数实现：
bt1＝find(bus_type＝＝PV) (21)
bt2＝find(bus_type～＝Vθ) (22)
式中，bus_type为节点类型列向量；～＝为不等于关系运算符；＝＝为等于关系运算符；Vθ为平衡节点类型；PV为PV节点类型。
形成数组bt2后，把数组bt2的所有元素都加上节点数n后，再添加到原数组bt2后形成数组bt，用来记录雅可比矩阵及方程右端向量需要的行列号：
bt＝[bt2 bt2+n] (23)
C、形成节点导纳矩阵，并转化为稀疏矩阵Y；
D、形成雅可比矩阵及计算节点功率；
采用Matlab编程，推导出基于稀疏矩阵技术的雅可比矩阵的矩阵运算和复数运算计算方法。
雅可比矩阵元素与节点类型有关，常规形成雅可比矩阵时要判断节点类型，根据节点类型决定哪些节点需要形成雅可比矩阵元素。这样处理对于通过循环来实现的算法，容易处理，但不适合通过矩阵整体运算形成雅可比矩阵的方法。因此，本发明形成雅可比矩阵时，不判断节点类型，对所有节点都形成雅可比矩阵元素，之后再去掉不需要的行和列。
下面推导采用矩阵运算计算雅可比矩阵元素和节点功率的公式。
对式(7)～式(10)的分析，得：
Mij＝Nij (24)
Lij＝-Hij (25)
因此，先求Hij和Nij，求出Hij和Nij后，自然就能得到Mij和Lij。
推导由矩阵运算计算雅可比矩阵的公式之前，先看看雅可比矩阵各元素如何用复数或相量表示。
根据式(7)和式(8)，雅可比矩阵各元素为Yij与的共轭的乘积，即因此Hij、Nij、Mij、Lij由下式生成：
式中，上标(^)表示复数的共轭。
由式(26)元素组成的矩阵为：
式(27)的矩阵元素为电压相量的共轭和导纳矩阵元素的乘积，矩阵由下式得到：
式中，J0为雅可比初始计算矩阵；为节点电压共轭值列向量形成的稀疏对角矩阵；Y为稀疏导纳矩阵。
由J0得到初始雅可比矩阵分块子矩阵为：
H0＝-Re(J0) (29)
N0＝Im(J0) (30)
M0＝Im(J0) (31)
L0＝Re(J0) (32)
式中，H0、N0、M0、L0为初始雅可比矩阵的分块子矩阵；Re表示取矩阵元素的实部；Im表示取矩阵元素的虚部。
由式(29)～式(32)得到的初始雅可比矩阵分块子矩阵的非对角元素已经是雅可比矩阵元素了，对角元素还需要修正。
式(13)～式(16)中等式右侧的第1和第2项就是H0、N0、M0、L0的对角元素，因此只需对得到的H0、N0、M0、L0用ai和bi修正，即得到雅可比矩阵分块子矩阵对角元素。
修正雅可比矩阵分块矩阵对角元素要用到节点电流相量的实部ai和虚部bi，节点电流相量为：
式中，为节点电流相量列向量；Y为稀疏导纳矩阵；为节点电压列向量。
复功率列向量为节点电压列向量与节点电流相量的共轭值列向量对应行的元素相乘如下：
式中，为节点复功率列向量；为节点电压相量列向量；为节点电流相量的共轭值列向量；.*表示两向量对应行的元素相乘。
由于ai、bi分别为节点i的节点电流相量的实部和虚部，因此
式中，为节点i的节点电流相量。
用节点电流相量对雅可比矩阵分块子矩阵对角元素进行修正如下：
式(37)～式(40)写成矩阵运算形式为
式中，为节点电流相量形成的稀疏对角矩阵。
式(29)～式(32)代入到式(41)～式(44)，得
电压平方不平衡量对电压偏导的非对角元素为0，对角元素为：
Rii＝-2ei i＝m+1,…,n-1 (49)
Kii＝-2fi i＝m+1,…,n-1 (50)
式(49)和式(50)写成矩阵形式如下：
式中，为节点电压相量形成的稀疏对角矩阵。
用R、K中PV节点对应的行替换M、L中PV节点对应行：
M(bt1)＝R(bt1) (53)
L(bt1)＝K(bt1) (54)
形成雅可比矩阵及计算节点功率，包括以下步骤：
D1、计算雅可比初始计算矩阵J0；
D2、计算节点电流相量列向量
D3、计算节点复功率列向量
D4、由J0和计算雅可比矩阵分块子矩阵H、N、M和L；
D5、计算雅可比矩阵分块子矩阵R、K；
D6、用R、K修正M、L；
D7、由雅可比矩阵分块子矩阵H、N、M和L形成雅可比矩阵；
E、计算节点功率不平衡量和节点电压平方不平衡量；
将式(4)计算节点功率不平衡量及电压平方不平衡量方程写成矩阵运算的形成为：
式中，ΔP、ΔQ分别为节点有功功率不平衡量列向量和无功功率不平衡量列向量；ΔU2为节点电压平方不平衡量列向量；Ps、Qs分别为节点给定的注入有功功率列向量和注入无功功率列向量；Us为节点电压给定值列向量；为节点电压列向量。
计算最大有功功率不平衡量ΔPmax、最大无功功率不平衡量ΔQmax和最大电压平方不平衡量
用ΔU2中PV节点对应的行替换ΔQ中PV节点对应行：
ΔQ(bt1)＝ΔU2(bt1) (57)
F、解方程及修正电压实部e和虚部f；
由步骤D得到雅可比矩阵J和步骤E得到节点功率偏差向量ΔP和ΔQ，构造成潮流计算的修正方程如下：
直接调用Matlab软件的解线性方程组算法解修正方程组(58)，求出节点电压实部修正量列向量Δe及节点电压虚部修正量列向量Δf。调用Matlab软件的解线性方程组算法时，用数组bt去掉雅可比矩阵中不需要的行和列及不平衡量中不需要的行。
对电压进行修正的式(19)和式(20)改写成矩阵形式为：
e(t+1)＝e(t)-Δe(t) (59)
f(t+1)＝f(t)-Δf(t) (60)
式中，上标(t)表示第t次迭代的值；Δe和Δf分别为节点电压实部修正量列向量和节点电压虚部修正量列向量；e和f分别为电压实部列向量和电压虚部列向量。
计算节点电压实部和电压虚部后，按下式形成节点电压相量
G、判断最大不平衡量|ΔP|max、|ΔQ|max和|ΔU2|max是否都小于收敛精度ε；如果都小于收敛精度ε，转步骤H，否则返回步骤D进行下一次迭代。
H、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。
与现有技术相比，本发明具有以下有益效果：
1、本发明提出的方法在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。
2、本发明提出的方法主要部分都采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能；使用矩阵运算也大大提高了计算速度。
3、本发明采用Matlab的稀疏矩阵技术，较大幅度地提高了计算速度，同时Matlab的稀疏矩阵使用非常方便，可以像全矩阵一样用行列号直接使用稀疏矩阵的元素，也不需要设计稀疏存储结构。
4、发明根据稀疏矩阵运算特点，修改了雅可比矩阵、复功率和计算不平衡量的计算公式，使得程序简化，速度提高。实践证明，本发明的方法既方便了科研人员对程序进行编写、修改和调试，同时计算速度也基本接近了在C语言平台上实现的速度，为科研人员的科研工作提供了一个优秀的分析工具。
附图说明
本发明共有附图4张。其中：
图1是现有直角坐标牛顿法潮流计算的流程图。
图2是现有直角坐标牛顿法形成雅可比矩阵的流程图。
图3是本发明直角坐标牛顿法潮流计算的流程图。
图4是本发明形成雅可比矩阵和计算节点功率的流程图。
具体实施方式
下面结合附图对本发明进行进一步地说明，按照图3-4所示流程对一个修改后的445节点实际系统算例进行了计算。
445节点实际大型电力系统有445个节点，544条支路，含有大量的小阻抗支路。为了对各种方法进行比较，把这些小阻抗支路改为正常阻抗支路以满足各种方法的要求。
采用本发明和已有专利方法对445节点实际系统算例进行了计算，计算时收敛精度为0.00001。两种潮流计算方法分别为：
方法1：中国专利ZL201610864281.4方法，未采用稀疏矩阵技术；
方法2：本发明方法。
两种方法的潮流计算的计算时间见表1，计算时间不包括数据读入和输出及支路功率计算的时间。
表1 两种直角坐标牛顿法潮流计算计算时间比较
| 潮流计算方法 | 计算时间(s) |
| 方法1 | 0.562 |
| 方法2 | 0.061 |
从表1可见，中国专利ZL201610864281.4方法没有采用稀疏矩阵技术，计算时间较长；本发明采用稀疏矩阵技术并对某些计算公式根据稀疏矩阵的特点进行改进可以明显提高计算速度，计算时间仅为现有专利方法的1/9；同时采用Matlab的稀疏矩阵技术比较简单方便。
本发明方法可以在任何版本的MATLAB编程语言实现，但建议使用较新版本的MATLAB语言。
本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。
Claims (1)
1.基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法，包括以下步骤：
A、原始数据输入和电压初始化；
电压初始化采用平启动，即PV节点和平衡节点的节点电压实部取给定值，PQ节点的节点电压实部取1.0；所有节点电压的虚部都取0.0；这里单位采用标幺值；
B、记录相关节点类型的节点号；
直角坐标牛顿法潮流计算的修正方程组方程个数及变量个数与电力系统的节点类型有关，ΔP方程组中没有平衡节点有功功率不平衡量对应的方程，ΔQ方程组中仅有PQ节点无功功率不平衡量对应的方程，ΔU2方程组中仅有PV节点电压平方不平衡量对应的方程；变量则不包含平衡节点的相角变量和电压幅值变量；
设置3个数组记录有关节点类型的节点号，其中数组bt1记录PV节点的节点号，数组bt2记录PQ节点和PV节点的节点号，数组bt记录雅可比矩阵及方程右端向量需要的行列号；
记录相关节点类型的节点号使用Matlab的find函数实现：
bt1＝find(bus_type＝＝PV) (1)
bt2＝find(bus_type～＝Vθ) (2)
式中，bus_type为节点类型列向量；～＝为不等于关系运算符；＝＝为等于关系运算符；Vθ为平衡节点类型；PV为PV节点类型；
形成数组bt2后，把数组bt2的所有元素都加上节点数n后，再添加到原数组bt2后形成数组bt，用来记录雅可比矩阵及方程右端向量需要的行列号：
bt＝[bt2bt2+n] (3)
C、形成节点导纳矩阵，并转化为稀疏矩阵Y；
其特征在于：
D、形成雅可比矩阵及计算节点功率；
D1、计算雅可比初始计算矩阵J0；
<mrow>
<msup>
<mi>J</mi>
<mn>0</mn>
</msup>
<mo>=</mo>
<msub>
<mover>
<mi>U</mi>
<mo>^</mo>
</mover>
<mi>D</mi>
</msub>
<mi>Y</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，J0为雅可比初始计算矩阵；为节点电压共轭值列向量形成的稀疏对角矩阵；Y为稀疏导纳矩阵；
D2、计算节点电流相量列向量
<mrow>
<mover>
<mi>I</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>=</mo>
<mi>Y</mi>
<mover>
<mi>U</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，为节点电流相量列向量；为节点电压相量列向量；
D3、计算节点复功率列向量
复功率列向量为节点电压列向量与节点电流相量的共轭值列向量对应行的元素相乘如下：
<mrow>
<mover>
<mi>S</mi>
<mo>~</mo>
</mover>
<mo>=</mo>
<mover>
<mi>U</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>.</mo>
<mo>*</mo>
<mover>
<mi>I</mi>
<mo>^</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，为节点电流相量的共轭值列向量；.*表示两向量对应行元素相乘；
D4、由J0和计算雅可比矩阵分块子矩阵H、N、M和L；
由J0和计算雅可比矩阵分块子矩阵为：
<mrow>
<mi>H</mi>
<mo>=</mo>
<mo>-</mo>
<mi>Re</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>J</mi>
<mn>0</mn>
</msup>
<mo>+</mo>
<msub>
<mover>
<mi>I</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>D</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>N</mi>
<mo>=</mo>
<mi>Im</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>J</mi>
<mn>0</mn>
</msup>
<mo>-</mo>
<msub>
<mover>
<mi>I</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>D</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>M</mi>
<mo>=</mo>
<mi>Im</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>J</mi>
<mn>0</mn>
</msup>
<mo>+</mo>
<msub>
<mover>
<mi>I</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>D</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>9</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>L</mi>
<mo>=</mo>
<mi>Re</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>J</mi>
<mn>0</mn>
</msup>
<mo>-</mo>
<msub>
<mover>
<mi>I</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>D</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>10</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，H、N、M、L为雅可比矩阵的分块子矩阵；为节点电流相量形成的稀疏对角矩阵；Re表示取矩阵元素的实部；Im表示取矩阵元素的虚部；
D5、计算雅可比矩阵分块子矩阵R、K：
<mrow>
<mi>R</mi>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<mi>Re</mi>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>U</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>D</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>K</mi>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<mi>Im</mi>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>U</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>D</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>12</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，为节点电压相量形成的稀疏对角矩阵；
D6、用R、K修正M、L；
用R、K中PV节点对应的行替换M、L中PV节点对应行：
M(bt1)＝R(bt1) (13)
L(bt1)＝K(bt1) (14)
D7、由雅可比矩阵分块子矩阵H、N、M和L形成雅可比矩阵J；
<mrow>
<mi>J</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>H</mi>
</mtd>
<mtd>
<mi>N</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>M</mi>
</mtd>
<mtd>
<mi>L</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
</mrow>
E、计算节点功率不平衡量和节点电压平方不平衡量；
按下式计算节点功率不平衡量和节点电压平方不平衡量：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>&Delta;</mi>
<mi>P</mi>
<mo>=</mo>
<msub>
<mi>P</mi>
<mi>s</mi>
</msub>
<mo>-</mo>
<mi>Re</mi>
<mrow>
<mo>(</mo>
<mover>
<mi>S</mi>
<mo>~</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&Delta;</mi>
<mi>Q</mi>
<mo>=</mo>
<msub>
<mi>Q</mi>
<mi>s</mi>
</msub>
<mo>-</mo>
<mi>Im</mi>
<mrow>
<mo>(</mo>
<mover>
<mi>S</mi>
<mo>~</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msup>
<mi>&Delta;U</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<msubsup>
<mi>U</mi>
<mi>s</mi>
<mn>2</mn>
</msubsup>
<mo>-</mo>
<mo>|</mo>
<mover>
<mi>U</mi>
<mo>&CenterDot;</mo>
</mover>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>16</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，ΔP、ΔQ分别为节点有功功率不平衡量列向量和无功功率不平衡量列向量；ΔU2为节点电压平方不平衡量列向量；Ps、Qs分别为节点给定的注入有功功率列向量和注入无功功率列向量；Us为节点电压给定值列向量；
计算最大有功功率不平衡量ΔPmax、最大无功功率不平衡量ΔQmax和最大电压平方不平衡量
用ΔU2中PV节点对应的行替换ΔQ中PV节点对应行：
ΔQ(bt1)＝ΔU2(bt1) (17)
F、解方程及修正电压实部e和虚部f；
由步骤D得到雅可比矩阵J和步骤E得到节点功率不平衡量向量ΔP和ΔQ，构造成潮流计算的修正方程如下：
<mrow>
<mi>J</mi>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&Delta;</mi>
<mi>e</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>&Delta;</mi>
<mi>f</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&Delta;</mi>
<mi>P</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>&Delta;</mi>
<mi>Q</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>18</mn>
<mo>)</mo>
</mrow>
</mrow>
2
直接调用Matlab软件的解线性方程组算法解修正方程组(18)，求出电压实部修正量向量Δe及电压虚部修正量向量Δf；调用Matlab软件的解线性方程组算法时，用数组bt去掉雅可比矩阵中不需要的行和列及不平衡量中不需要的行；
按下式对节点电压实部和虚部进行修正：
e(t+1)＝e(t)-Δe(t) (19)
f(t+1)＝f(t)-Δf(t) (20)
式中，上标(t)表示第t次迭代的值；e和f分别为电压实部列向量和电压虚部列向量；
计算节点电压实部和电压虚部后，按下式形成节点电压相量
<mrow>
<mover>
<mi>U</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>=</mo>
<mi>e</mi>
<mo>+</mo>
<mi>j</mi>
<mi>f</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>21</mn>
<mo>)</mo>
</mrow>
</mrow>
G、判断最大不平衡量|ΔP|max、|ΔQ|max和|ΔU2|max是否都小于收敛精度ε；如果都小于收敛精度ε，进行步骤H，否则返回步骤D进行下一次迭代；
H、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710557642.5A CN107181260B (zh) | 2017-07-10 | 2017-07-10 | 基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710557642.5A CN107181260B (zh) | 2017-07-10 | 2017-07-10 | 基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107181260A true CN107181260A (zh) | 2017-09-19 |
| CN107181260B CN107181260B (zh) | 2019-10-29 |
Family
ID=59845079
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710557642.5A Expired - Fee Related CN107181260B (zh) | 2017-07-10 | 2017-07-10 | 基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107181260B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107482636A (zh) * | 2017-09-22 | 2017-12-15 | 大连海事大学 | 一种电力系统潮流计算的支路功率矩阵计算方法 |
| CN108233382A (zh) * | 2017-11-29 | 2018-06-29 | 广西大学 | 一种提取直角坐标潮流方程雅可比矩阵的方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101662148A (zh) * | 2009-09-25 | 2010-03-03 | 大连海事大学 | 一种直角坐标牛顿法潮流计算的电压初值设置方法 |
| CN104037764A (zh) * | 2014-07-03 | 2014-09-10 | 大连海事大学 | 一种雅可比矩阵改变的直角坐标牛顿法潮流计算方法 |
| CN106356859A (zh) * | 2016-09-29 | 2017-01-25 | 大连海事大学 | 一种基于Matlab的直角坐标牛顿法潮流计算方法 |
| CN106602570A (zh) * | 2017-01-25 | 2017-04-26 | 大连海事大学 | 一种基于Matlab的快速分解法潮流计算方法 |
-
2017
- 2017-07-10 CN CN201710557642.5A patent/CN107181260B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101662148A (zh) * | 2009-09-25 | 2010-03-03 | 大连海事大学 | 一种直角坐标牛顿法潮流计算的电压初值设置方法 |
| CN104037764A (zh) * | 2014-07-03 | 2014-09-10 | 大连海事大学 | 一种雅可比矩阵改变的直角坐标牛顿法潮流计算方法 |
| CN106356859A (zh) * | 2016-09-29 | 2017-01-25 | 大连海事大学 | 一种基于Matlab的直角坐标牛顿法潮流计算方法 |
| CN106602570A (zh) * | 2017-01-25 | 2017-04-26 | 大连海事大学 | 一种基于Matlab的快速分解法潮流计算方法 |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107482636A (zh) * | 2017-09-22 | 2017-12-15 | 大连海事大学 | 一种电力系统潮流计算的支路功率矩阵计算方法 |
| CN107482636B (zh) * | 2017-09-22 | 2020-03-13 | 大连海事大学 | 一种电力系统潮流计算的支路功率矩阵计算方法 |
| CN108233382A (zh) * | 2017-11-29 | 2018-06-29 | 广西大学 | 一种提取直角坐标潮流方程雅可比矩阵的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107181260B (zh) | 2019-10-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106356859B (zh) | 一种基于Matlab的直角坐标牛顿法潮流计算方法 | |
| Ghatak et al. | A fast and efficient load flow technique for unbalanced distribution system | |
| CN102013680B (zh) | 一种电力系统快速分解法潮流计算方法 | |
| CN106602570B (zh) | 一种基于Matlab的快速分解法潮流计算方法 | |
| CN106229988B (zh) | 一种基于Matlab的极坐标牛顿法潮流计算方法 | |
| CN107196306B (zh) | 基于Matlab稀疏矩阵的快速分解法潮流计算方法 | |
| CN103018534B (zh) | 确定谐波电压的方法及系统 | |
| CN107194131A (zh) | 基于Matlab稀疏矩阵的极坐标牛顿法潮流计算方法 | |
| CN103810646B (zh) | 一种基于改进投影积分算法的有源配电系统动态仿真方法 | |
| CN104899396B (zh) | 一种修正系数矩阵的快速分解法潮流计算方法 | |
| CN106786493A (zh) | 一种多馈入直流相互作用因子的实用计算方法 | |
| CN106709243B (zh) | 含小阻抗支路电网的补偿法极坐标牛顿法潮流计算方法 | |
| CN101976838A (zh) | 一种适合研究目的使用的牛顿法潮流计算方法 | |
| CN107181260A (zh) | 基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 | |
| CN106532712B (zh) | 含小阻抗支路电网的补偿法直角坐标牛顿法潮流计算方法 | |
| Zhao et al. | Holomorphic embedding power flow for AC/DC hybrid power systems using Bauer's Eta algorithm | |
| CN104979840B (zh) | 一种主动配电网三相无功优化方法 | |
| CN106410811B (zh) | 首次迭代小阻抗支路端点改变雅可比矩阵的潮流计算方法 | |
| CN110417022A (zh) | 矩阵乘法运算提取雅可比元素的配电网三相潮流计算方法 | |
| CN107658880B (zh) | 基于关联矩阵运算的快速分解法系数矩阵计算方法 | |
| CN105703363A (zh) | 一种基于线电压的不接地配电网三相潮流计算方法 | |
| CN110336287B (zh) | 一种基于雅可比元素提取的配电系统三相潮流计算方法 | |
| CN110336288B (zh) | 基于矩阵运算提取雅可比元素的配电网三相潮流计算方法 | |
| CN107944682B (zh) | 基于Matlab矩阵运算的潮流计算导纳矩阵计算方法 | |
| CN107834562B (zh) | 基于Matlab矩阵运算的快速分解法系数矩阵计算法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20191029 Termination date: 20200710 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |