武汉轻工大学毕业设计（论文）设计(论文)题目：基于潮流计算的稀疏技术研究
姓名曹乐
学号 100408521
院（系）电气与电子工程学院
专业电气工程及其自动化
指导老师仰彩霞
2014 年 6 月 2 日
目录
摘要 ............................................................................................................................... I ABSTRCT .................................................................................................................... I I 1 绪论 (1)
1.1 课题研究目的和意义 (1)
1.2 稀疏技术在电力系统潮流计算中的研究 (1)
1.3 本文完成的主要工作和章节安排 (2)
2 电力系统潮流计算 (4)
2.1 电力系统潮流计算的数学模型 (4)
2.1.1 节点电压方程 (4)
2.1.2 功率方程 (5)
2.2 电力系统潮流计算的方法 (5)
2.2.1 高斯-赛德尔潮流计算 (6)
2.2.2 牛顿-拉夫逊潮流计算 (8)
3 稀疏存储技术 (15)
3.1 静态数组结构存储方法 (15)
3.1.1 散居格式 (15)
3.1.2 按行（列）存储格式 (15)
3.1.3 三角检索存储格式 (16)
3.2动态数组结构存储方法 (17)
3.2.1 带行指针数组的单链表表示法 (17)
3.2.2 十字链表表示法 (18)
3.2.3 二叉单元链表数组 (19)
4 用于牛顿-拉夫逊潮流计算的稀疏技术 (21)
4.1 节点分块雅可比矩阵 (21)
4.2 二层链表的应用 (22)
4.3 二层链表的改进 (23)
4.4 实验与结果 (24)
总结 (25)
致谢 (29)
参考文献 (30)
附录 (30)
摘要
潮流计算是电力分析最基础的电气运算之一，其效率一直是研究人员关心的问题。目前，稀疏技术已经广泛应用于工程计算、科学分析等诸多领域。因此,结合潮流计算方法的特点进一步提高稀疏技术的应用效率，对电力系统各种分析计算都具有一定意义。创建一种特殊的数据存储结构，即二层链表，该结构可以实现导纳矩阵与节点分块雅可比矩阵的同结构存储，以及导纳矩阵各元素到雅可比矩阵的定位、查询功能，提高了雅可比矩阵的形成与修正效率。由于导纳矩阵和雅克比矩阵具有对称性，可进一步将二层链表存储单元数量减半，形成下三角二层链表结构，能更好的节省计算机的存储空间；同时，减少了待修正雅可比矩阵与功率不平衡中共同元素的运算量，提高潮流计算效率。
关键词：潮流计算，稀疏技术，二层链表
ABSTRCT
Power flow calculation is one of most basic eletrical operation applied in the electrical system analysis,its efficiency has been a matter of concern to researchers.Currently,the sparse technique has been widely used in the areas of engineering calculations,scientific analysis.Therefore,to improve the efficiency of the sparse technique by combing the power flow calculation features could enhance the efficiency of the electrical operations.We proposed a special kind of two-layer linked list to make the bus admittance matrix and the bus block jacobian matrix stored,and the function of location searching between the two matrix implemented under the same structure.And the efficiency of the foemation and modification of jacobian matrix has benn largely improved.The research then used the symmetry of the symmetric matrix during the factorization,halved the amount of the layer two linked list storage cell,and formed the lower triangular and two-layer linked list.So the storage space could be better used.At the mean time,through exacting the elements from the being revised jacobian matrix and the unbalanced power to reduce the operation so as to enhance the efficiency.
Keywords:Power Flow Calculation,Sparse Technique,Linked List
1 绪论
1.1 课题研究目的和意义
潮流计算是电力系统分析中最根本的电气运算。潮流计算主要是用来研究电力系统规划和网络运行过程中提出的各种问题。对于规划中的电力系统，通过潮流计算可以有效的检验所提出的规划方案是否能满足运行的要求；对于运行中的电力网络，通过潮流计算可以立刻确定各个负荷的变化和部分网络结构改变对电力系统运行稳定性的影响、系统中各条母线电压是否在允许的范围内、以及系统中各种元件(线路电缆、变压器等)是否会出现过负荷工作。
为了提高电力系统潮流计算的速度和效率，提出了基于潮流计算的稀疏技术研究。稀疏数据是一类零元素所占比例很大，非零元素所占比例较小的数据。一般而言，数据的有效运算是非零元素之间的运算，但稀疏数据中大量的零元素将会给计算带来大量的冗余度，降低计算机算法的速度和效率。在电力系统分析与计算中存在大量具有稀疏性的矩阵，如潮流计算中的雅可比矩阵。因此，充分利用电力系统中数据的稀疏性，运用相关的稀疏算法，可以有效的提高电力系统潮流计算的速度和效率。
1.2 稀疏技术在电力系统潮流计算中的研究
稀疏技术已被广泛的应用于电力系统潮流计算中，受到了广大电力专家、学者的高度重视并且取得了大量的研究成果。
文献[1]就曾提出过基于二维链表的稀疏存储方法，该文献的作者朱凌志、安宁成功地将该方法应用于潮流计算中。通过对传统二维链表的存储结构的改进、预先在稀疏矩阵中增加相关的冗余元素、优化潮流计算电压和功率方程的结构，实现了对潮流计算方程的改进。文献[2]通过研究数组存储方式和链表存储方式在牛顿-拉夫逊潮流计算中的效率问题、以及内存占用情况，证实链表存储方式在内存空间占用、计算效率上有明显的优势。文献[3]提出了一种基于十字链表动态生成的稀疏技术，这种十字链表能同时存储电力网络中的拓扑结构信息和参数信息，能提高电力网络中潮流计算的效率，同时还能在运行状态变化及故障时对电力网络结构进行有关的修改。为了提高牛顿-拉夫逊潮流计算效率，文献[4]引入了基于消去树理论中的符号因子分解技术以及改进的LU数值分解算法。在符号因子分解中引入消去树理论是对LU分解可能出现的填充元位置进行定位，避免在之后的数值计算时频繁的插入非零元素。文献[5]提出来一种静态存储方法VEPD，该方法有利于存储具有稀疏性质的矩阵。作者在文章指出，采
用VEPD法存储矩阵所占用的内存空间较小，且元素的查找、修改、删除等操作方便快捷。
另外，稀疏技术对配电网潮流计算效率的提高也有一定的意义。为了减小大型配电网潮流计算的冗余度，文献[6]分析和研究了配电网的拓扑结构特点，作者提出了利用树状链表和递归搜索方法得出反映配电网树状链表关系的数据结构。该数据结构具有存储过程不受输入数据影响、直观反映电力网络结构、避免存取大规模矩阵等优点，从而有效降低配电网潮流计算时间。文献[7]通过理论研究，推导出基于逆流编号法的配电网牛顿-拉夫逊法潮流过程与等值递推法相一致的结论，为实现配电网潮流计算方法多样性与准确性提供了理论依据。文献[8]提出了一种在辐射网络的导纳矩阵消去分解过程中不出现注入元的节点优化编号方法，即逆流编号法。逆流编号法指出：任意网络节点的编号一定要大于其顺流节点的编号或者任意节点的编号必须小于其逆流节点的编号。作者经研究证明该方法能有效降低辐射电网节点导纳矩阵的运算量，节约计算机内存空间。在文献[9]中，作者系统地分析和对比各种配电系统节点编号方案，对采用不同潮流计算方法（节点电压法、支路电流法）所对应的节点编号方案进行了概括、分析，此文献对配电网潮流计算有一定的指导意义。
1.3 本文完成的主要工作和章节安排
本文是以潮流计算中的节点导纳矩阵与节点分块雅可比矩阵作为研究对象，利用稀疏存储技术，实现对上述两个矩阵同结构存储。本文的理论部分，详细地介绍潮流计算的数学模型、计算方法以及稀疏存储技术；在潮流计算程序部分，作者将以牛顿-拉夫逊潮流计算法编写程序，并进行Matlab程序仿真。
本文的章节具体安排如下：
第一章绪论
第二章电力系统潮流计算
第三章稀疏技术
第四章用于牛顿-拉夫逊潮流计算的稀疏技术
第五章结论
注：
1.第一章写出本文的研究目的和意义，并对目前稀疏技术在电力系统潮流计算中的研究现状加以介绍；
2.第二章写出了电力系统潮流计算的原理，介绍了两种潮流计算的方法，其中作者详细的介绍了牛顿-拉夫逊潮流计算法；
3.第三章作者介绍了稀疏技术内容，着重对各类稀疏存储技术进行了详细的讲解与说明；
4.第四章作者将稀疏技术运用到电力系统的潮流计算中，旨在提高计算的速度和效率，并且给出Matlab编写的算法实例；
5.第五章总结。
2 电力系统潮流计算
所谓电力系统潮流计算，就是在已知电网接线方式、参数及运行条件的情况下，计算出电力系统稳态运行时各母线电压、各支路电流、功率及网络损耗。
2.1 电力系统潮流计算的数学模型
电力网络的数学模型是由网络的有关参数和变量所组成的可反映网络性能的数学方程式组，如节点电压方程和回路电流方程。电力系统潮流计算中多使用节点电压方程，原因是节点电压方程所描述的运行状态直观、易于处理运行条件的变化，且方程个数少于回路电流方程。
2.1.1 节点电压方程
由电路原理可知，对于具有n 个独立节点的电力系统，可写出n 个节点电压的矩阵方程为
?
???????
???????????????
???????????
????
??????????????????????????=????????
??????????????
????
???n nn n n n n n n n U U U U Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y
I I I I 321321333323122322211131211321 (2.1)
简记为
YU I =
(2.2)
式中 T
n I I I I I ??
?
??
????=?
?
?
?
3
2
1
——节点注入电流的转置列向量，
其阶数为1?n ； T
n U U U U U ??
?
??????=?
???321——节点电压的转置列向量，其阶数为1?n ；
????
?
???
?????
??????????????????????????=nn n n n n n n Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y Y 3
21333
323122322211131211——节点导纳矩阵，其阶数为n n ?。
2.1.2 功率方程
电力网络的节点电压方程YU I =就是潮流计算时的数学模型。如果已知各节点的电流I ，直接求解线性的节点电压方程YU I =即可。但是，工程中已知的条件既不是节点电压U ，也不是节点电流I ，而是各节点的注入功
率 。将 代入上述节点电压方程，并写成其展开形式为
∑=?*
*
=+n j j ij i i i U Y U jQ P 1
)( ),,2,1(n i ???= (2.3)
式（2.3）为一非线性方程组，通常称为功率方程。i P 、i Q 分别为节点i 向网络注入的有功功率和无功功率。
功率方程随着节点电压向量表示形式不同有两种不同的形式。 ①直角坐标形式的功率方程
若节点电压向量以直角坐标形式即i i i jf e U +=?
表示，导纳用ij ij ij jB G Y +=形式表示，则将这两个关系代入式（2.3）后展开，并将功率的实部和虚部分别列写，即可得到分离的有功、无功功率方程为
)
()()
()(1
1
11j ij j ij n
j i j ij j ij n
j i j ij j ij n
j i j ij j ij n
j i i e B f G e f B e G f Q e B f G f f B e G e P +--=-+-=∑∑∑∑==== (2,4)
②极坐标形式的功率方程
若节点电压以极坐标形式即i
j i i e U U δ=?
或i i i i i jU U U δδsin cos +=?
表示，并将其导纳的复数表示式一起代入（2.3）的功率方程，经整理后也可分别写成有功、无功的方程式为
)
cos sin ()
sin cos (1
1ij ij ij ij n
j j i i ij ij ij ij n
j j i i B G U U Q B G U U P δδδδ-=+=∑∑== (2.5)
式中 j i ij δδδ-=——i 节点电压与j 节点电压的相角差。
2.2 电力系统潮流计算的方法
潮流计算的功率方程是一组非线性代数方程，电力系统潮流计算就是求解这组非线性方程组。求解非线性方程组在数学上有两种途径，即采用迭代法直接求解非线性方程组或将非线性方程组线性化后迭代求解。直接求解非线性方程组有
T
n S S S S S ??
????=???? 321**?=i i i U S I
雅可比迭代法、高斯-赛德尔法等；将非线性方程组线性化后迭代求解的方法主要有牛顿-拉夫逊法等。
早期，潮流计算以高斯-赛德尔迭代法较为普遍。高斯-赛德尔迭代法是数学上求解线性或非线性方程组的一种常用的迭代方法。在潮流计算中，高斯-赛德尔迭代法可分导纳矩阵迭代法和阻抗矩阵迭代法两种，前者是以节点导纳矩阵为基础建立的赛德尔迭代格式;后者是以节点阻扰矩阵为基础建立的赛德尔迭代格式。
目前，牛顿-拉夫逊法和快速解耦法是两种应用最为普遍的方法，也是其他潮流计算方法的基础。牛顿-拉夫逊法的原理是将非线性方程组转化为线性方程组，通过迭代求解线性方程组得到解，牛顿-拉夫逊法具有较好的收敛性。快速解耦法是由牛顿-拉夫逊法极坐标形式简化而来的，它以节点有功功率偏差作为电压相角的修正量，以节点无功功率偏差作为电压幅值的修正量，有功功率和无功功率迭代是分开进行的。快速解耦法的修正方程系数矩阵的维数较低且为常数，因而计算速度较快，占用的内存更少，被广泛使用。
潮流计算的研究广泛而活跃，为改进牛顿-拉夫逊法的某些不足，采用了包括泰勒级数高阶项的更精确模型——保留非线性潮流计算[]11,10，一定程度上提高了潮流计算的收敛性能。另外，专家、学者以及研究人员针对一些特殊的潮流计算开发了对应的方法，例如直流潮流算法[]12、随机潮流算法[]13、三相潮流算法[]14、谐波潮流算法[]15以及用于交直流混联的潮流算法[]16。随着人工智能理论的发展，人工神经网络[]17、遗传算法[]17和模糊算法[]19等也逐渐被引用到潮流计算中。但是目前为止，无论新的或是改进的算法，仍不能取代牛顿-拉夫逊法和快速解耦法的地位。
2.2.1 高斯-赛德尔潮流计算
用高斯-赛德尔法计算电力系统潮流，首先需将功率方程式（2.3）改写，然后经行迭代计算。
假设待计算的系统有n 个独立节点，其中1个平衡节点，m 个PQ 节点，
)1(+-m n 个PV 节点，则由功率方程式（2.3）可解得?
i U 为
????
??????--=∑≠=?*?
n
i j j j ij i i i ii i U Y U jQ P Y U 1
1 (2.6)
再将式（2.6）改写为高斯-赛德尔迭代法的格式，即
??
????---=∑∑+=?+?-=*+?n
i j k j ij k j i j ij k i i i ii k i
U Y U Y U jQ P Y U 1)()1(11)()
1(1 (2.7)
则可对上式迭代进行潮流计算。
若有一网络，其节点1为平衡节点，其余1-n 个节点都为PQ 节点。由于平衡节点不参与迭代，则用高斯-赛德尔迭代法潮流计算的公式为
??
????----=??
????----=??
????----=??
????----=????*+?????*+?????*+?????*+?)(4)(33)(2211)()
1()(4)(343)(242141)
(44
444
)
1(4)(1)(434)(232131)(33
333)
1(3
)(1)(441)(331121)(22
222)
1(21111k n n k n k n n k n n
n nn k n
k n n k k k k k n n k k k k k n n k k k k U Y U Y U Y U Y U jQ P Y U U Y U Y U Y U Y U jQ P Y U U Y U Y U Y U Y U jQ P Y U U Y U Y U Y U Y U jQ P Y U （2.8）
式中
?
1U ——平衡节点1的电压； k ——迭代次数。
对于PQ 节点，由于其功率是给定的，故只要给出节点电压的初始值 ， 即可按式（2.7）进行迭代计算。
对于PV 节点，因其无功功率是未知量，只能在迭代开始时给定初始值)0(i Q ，此后的迭代值必须在逐次的迭代过程中计算得出。因此，对PV 节点进行每一次迭代计算，必须要完成一下几步计算：
①修正节点电压。按式（2.7）计算所得到的电压并不一定刚好等于PV 节点所要求的电压幅值is U 。为满足此条件，可只保留节点电压的相位角k i δ，而将其幅值用给定值is U 代替进行修正，即
)()
(k i is k i U U δ∠=?
(2.9)
②计算节点的无功功率。第k 次迭代的节点无功功率计算式为
??
????+=???
???=∑∑-==**+**?*?)(Im Im 11)()1()()()()
(i j n i j k ij ij k j ij k i k i k i k i
U Y U Y U I U Q
(2.10)
③用式（2.10）计算所得的无功功率)(k i Q ，应该满足下面的约束条件，即
max )(min i k i i Q Q Q ≤≤
如果min )(i k i Q Q <则应令min )(i k i Q Q =；如果max )(i k i Q Q >，则应令max )(i k i Q Q =。 上述情况表明无功功率已达到极限值，不能再用调节无功功率来保持PV 节点的电压值。此时由于无功功率已限定为极限值，因而PV 节点就转变为PQ 节点，故计算应按PQ 节点进行。
)
0(i U ?
完成上述三步计算后，才可应用式（2.7）计算节点电压的新值。对于平衡节点，由于电压的幅值和相位角为给定值，故无需迭代计算。
每次迭代完成后，应根据给定的任意小数ε，用下面的迭代收敛判据检验，即
ε≤-=??+??
max
)
()
1(max k i k i
i U U U
如果该式满足时，停止迭代计算，)1(+k i U 或)(k i U 即为所求得的节点电压。 求出节点电压后，即可计算各条线路的潮流、平衡节点功率以及各支路的功率损耗。由图2.1可推出线路中的功率计算公式如下：
[]
[]
ji
ji ji i j j j j ij j ji ij ij ij j i i i i ij
i ij jQ P y U U y U U I U S jQ P y U U y U U I U S +=-+==+=-+==*****?
*?
*
****?
*
?
)()(00
(2.11)
图2.1 线路功率计算示意图
由于式（2.8）中的1=S 为平衡节点，故平衡节点的功率可以S i =代入式（2.8）求出。即
∑=**?
=+=n
j j sj s s s s U Y U jQ P S 1
(2.12)
各线路上损耗的功率为
ij ij ji ij ij Q j P S S S ?+?=+=?
(2.13)
2.2.2 牛顿-拉夫逊潮流计算
牛顿-拉夫逊潮流计算的核心问题是雅可比修正方程的建立和求解。为说明这一修正方程式的建立过程，首先对网络中各类节点的编号作如下的约定：
（1）网络中共有n 个节点，编号为n ,,3,2,1???，其中1个平衡节点，编号为s 。 （2）网络中有m 个PQ 节点，编号为m ,,3,2,1???，其中包括编号为s 的平衡节点。 （3）网络中有)1(+-m n 个PV 节点，编号为n m m ,,2,1???++。 直角坐标形式的牛顿-拉夫逊潮流计算
式2.1-2.3给出了节点的注入有功功率、无功功率与电压幅值方程。式中，i P 、
i Q 与i U 分别为节点i 的注入有功功率、无功功率与电压幅值；i e 与i f 分别为节点i 电压的实部与虚部，并且i i i jf e U +=；ij G 与ij B 分别为支路j i -的电导与电纳，ii G 与ii B 为节点i 的自电导与自电纳。
[]
[]
2
221
1)()()()(i i i i n
j j ij j ij i j ij j
ij
i
i n
j j ij j ij i j ij j
ij
i
U f e Q e B f G e f B e
G f P e B f G f f B e
G e =+=+--=++-∑∑==
(2.13)
对于PQ 节点，其功率误差方程可表示为
[][
]
∑∑==++--=-=?++--=-=?n j j ij j ij i j ij j ij i i i is i n
j j ij j ij i j ij j ij i is i is i e B f G e f B e G f Q Q Q Q e B f G f f B e G e P P P P 1
1)
()()
()( (2.14)
式中
is P 、is Q ——节点i 的给定有功功率、无功功率。
对于PV 节点，其有功功率、电压误差方程可表示为
[
])
()
()(222
2221
i i is i is i n
j j ij j ij i j ij j ij i is i is i f e U U U U e B f G f f B e G e P P P P --=-=?++--=-=?∑=
(2.15)
由于平衡节点只有一个，电压不参加迭代，其电压为
s s s jf e U +=?
(2.16)
上述的功率误差方程和电压误差方程构成的方程组即为牛顿-拉夫逊法潮流计算所要求解的非线性方程组。非线性方程组中的待求量为各节点的电压的实部
i e 和虚部i f 。除了一个平衡节点的电压已知外，共有)1(2-n 个未知数待求。在具有n 个独立节点的电力系统中，PQ 节点有m 个，PV 节点有)1(+-m n 个，则有功功率误差方程i P ?共有)1(-n 个，无功功率误差方程i Q ?共有m 个，电压误 差方程2i U ?共有)1(+-m n 个所以上述非线性方程组中共有方程)1(2-n 个，与方程中的未知个数相同，方程组有唯一的非零解。
把各个节点的电压变量用初始值与修正值的形式表示为
)
0()
0()
1()0()0()1(i
i
i
i i i f f f e e e ?-=?-=
将此关系式带到式（2.14）、式（2.15）中，在初值)1(i e 、)1(i f 附近的)0(i e ?、)
0(i f ?范围内将其展开为泰勒级数，略去含有)0(i e ?、)0(i f ?的二次以上各项就可得到修正方程。省略掉迭代次数的符号，可写出其矩阵形式为
?
??
???
?
?
?
??????????
???????????????????????????????
?????????
????????
??
?
??
???????
???????
????????????????
???????????????
??????????????????????????????????????????????????????????????????
??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????=????????????????
???????????????????????????--++----+-+---------+-+------+-+++++++++-+-+++++++++--++--++--++--++--++11111112
11
2
11
2
11
2
12
12
11
2
1
1
2
1
111
11
11
1111
11
1
12
112
112
112
12
1
2
112
1
1
2
1
111
11
11
1111
111111111111111111111111111111111111111111
1
21121
111n n m m m m n n n n m n m n m
n m
n n n
n n n n m n m n m
n m
n n n n m n m m m m m m m m m m m n m n m m m m m m
m m
m m m n m n m m m m m m m m m m m n m n m m m m m m m m m m m
n n m m m m n n m m m m n n m m m m f e f e f e f e f U e U f U e U f U e U f U e U f P e P f P e P f P e P f P e P f U e U f U e U f U e U f U e U f P e P f P e P f P e P f P e P
f Q e Q f Q e Q f Q e Q f Q e Q f P e P f P e P f P e P f P e
P f Q e Q f Q e Q f Q e Q f Q e Q f P e P e P e P f P e P f P e
P U P U P Q P Q P
(2.17)
式中雅可比矩阵的各个元素分别为：
非对角线元素)(j i ≠时：
)()
()()(2
2
=???==???=+=??=
-=??=-=??=
+-=??=j
i ij j
i ij i ij i ij j i
ij i ij i ij j
i ij i ij i ij j i ij i ij i ij j i ij f U S e U R f B e G f Q L f G e B e Q J f G e B f P N f B e G e P H (2.18)
对角元素)(j i =时：
i
j
i ij i
j
i ij i ij i ij n j j ij j ij j
i
ij i ij i ij n j j ij j ij j
i ij i ii i ii n j j ij j ij j i ij i ii i ii n
j j ij j ij j i ij f f U S e e U R f B e G f B e G f Q L f G e B e B f G e Q J f G e B e B f G f P N f B e G f B e G e P H 22)()()
()()()()()(2
2
1
111-=???=-=???=++--=??=
-++=??=-++-=??=
+---=??=∑∑∑∑==== (2.19)
分析上述表达式可以看到，雅可比矩阵具有以下特点：
1）各元素是节点电压的函数，每迭代一次各节点电压都要发生改变，因而各元素也随之发生变化；
2）非对称矩阵；
3）互导纳0=ij Y 时，与之对应的非对角线元素亦为零，此外因非对角线元素
，故雅可比矩阵是稀疏矩阵。 牛顿-拉夫逊潮流计算的主要流程可简述如下： a. 形成导纳矩阵；
b. 设置各节点电压初始值)0(i e 、)0(i f
c. 将初始值代入式（2.14）、（2.15），计算各节点功率及电压的偏移量)0(i P ?、
)0(i Q ?、2)0(i U ?；
d. 运用式（2.18）、式（2.19）计算雅可比矩阵中的各个元素；
e. 解修正方程式（2.17），求出各个节点电压的修正量)0(i e ?、)0(i e ?
f. 求新的电压初始值，其公式为
)
0()
0()
1()0()0()1(i
i
i
i i i f f f e e e ?-=?-=
g. 用新的初始值代入式（2.14）、式（2.15），计算新的各个节点功率及电压偏移)1(i P ?、)1(i Q ?和2)1(i U ?；
h. 检查计算是否已经收敛，由式（2.14）可知，如电压趋近与真实解时，其功率偏移量趋向零。因而其收敛判据为
ε
(max )(max max
)(,)
(k k k i Q P x f
其中ε为预先给定的小数。若不收敛返回到第（3）步重新迭代，若收敛进行下一步。
i. 计算各线路中的功率分布及平衡节点功率，最后打印出来计算的结果。各线路中的功率平衡节点功率可安式（2.12）及式（2.13）计算。
牛顿-拉夫逊潮流计算直角坐标式的流程图如图2.3所示。
022=???=???j
i j i f U e U
图2.3 牛顿-拉夫逊潮流计算直角坐标式的流程图
开始
形成节点导纳矩阵
输入原始数据 设节点电压(0)(0)i i e f ，
≠置迭代次数0k = 置节点号i=1 按式（2-18），（2-19）计算雅克比按式（2-14）计算PQ 节点的()
k i P ?，()k i Q ?，PV 节点的()k i P ?，k i U ?求解修正方程式，得()k i e ?，
()k i f ? 雅克比矩阵是否已全部形
计算平衡节点及PV 节点功求()max ||k e ?，()max ||k f ? 迭代次数 i=i+1
()()max max ||,||k k e f ε??≤潮流计算完成
计算各节点电压的新值：
(1)()()k k k i e e e +=+? k k k
i f f f +=+?
极坐标形式的牛顿-拉夫逊潮流计算
节点电压不仅可以用直角坐标形式表示，也可以用极坐标形式表示。因此，牛顿-拉夫逊潮流计算时的修正方程还有另一种表达形式。为建立极坐标形式的修正方程式，可仍令ij ij ij jB G Y +=，但令)sin (cos i i i j i j U e U U i δδδ+==?
（i δ为第
i 个节点的电压相角，ij δ为第i 个节点与第j 个节点的电压相角之差：
j i ij δδδ-=）。由此，可以将式（2.14）改写为
∑∑==--=-=?+-=-=?n
j ij ij ij ij j i is i is i n
j ij ij ij ij j i is i is i B G U U Q Q Q Q B G U U P P P P 1
1)
cos sin ()
sin cos (δδδδ (2.20)
对于具有n 个独立节点，其中有)1(+-m n 个PV 节点的网络，式（2.20）组成的方程中共有1-+m n 个方程式。采用极坐标表示时，方程组个数较采用直角坐标表示少)1(+-m n 个。因对PV 节点，采用极坐标表示时，待求的只有电压的相角i δ和注入无功功率i Q ，而采用直角坐标表示时，待求的有电压的实数部分i e 、虚数部分i f 和注入无功功率I Q 。所以极坐标形式的牛顿-拉夫逊潮流计算的未知变量少了)1(+-m n 个，方程数也相应少了)1(+-m n 个。这样可建立修正方程式的矩阵形式为： ??
?
???????????-=????????-U U L K N H
Q P 1δ (2.21)
式中
?????
?
????????????=?-121n P P P P
?????
?
????????????=?m Q Q Q Q 21
????
?
?
????????????=?-121n δδδδ ???
?
?
?????????????=?m U U U U 21
?
?????
?????
?=m U U U U 21 H 是)1()1(-?-n n 阶方阵；N 是m n ?-)1(阶矩阵；K 是)1(-?n m 阶矩阵；L 是
m m ?阶方阵。各矩阵的元素分别为：
j
j
i
ij j i
ij j j
i
ij j i
ij U U Q L Q J U U P N P H )(
)(???=???=
???=???=
δδ
(2.22)
将式（2.20）代入（2.22）求偏导数可得雅可比矩阵中各元素的表达式如下： 非对角线元素)(j i ≠：
)
cos sin ()(
)sin cos ()
sin cos ()(
)cos sin (ij ij ij ij j i j j
i
ij ij ij ij ij j i j
i
ij ij ij ij ij j i j j
i
ij ij ij ij ij j i j
i
ij B G U U U U Q L B G U U Q J B G U U U U P N B G U U P H δδδδδδδδδδ--=???=+=???=
+-=???=--=???=
(2.23)
对角线元素)(j i =：
∑∑∑∑====--=-=+-=-=+--=--=-+=+=n
j ij ij ij ij j i ii i
i ii i
ii n
j ij ij ij ij j i ii i i ii i ii n
j ij ij ij ij j i ii i
i ii i
ii n
j ij ij ij ij j i ii i
i ii i
ii B G U U B U Q B U L B G U U G U P G U J B G U U G U P G U N B G U U B U Q B U L 1
221221
221
22)
cos sin ()
sin cos ()
sin cos ()
cos sin (δδδδδδδδ (2.24)
3 稀疏存储技术
在科学计算与工程分析中，常常需要处理一类稀疏性的数据，这类数据所包含的零元素比例较大，而非零元素所占比例小，我们将此类数据称之为稀疏数据。
所谓稀疏技术，是指通过设计算法和编制程序，尽可能避免储存稀疏数据中的零元素以及对这些元素进行计算。
稀疏存储技术的核心内容是只存储与非零元素相关的信息，从而节省储存空间并提高计算机的运算速度。具体的操作方式有很多，如果按照存储原理可以分为两大类：一是静态数组结构存储方法，二是动态链式结构存储方法。
3.1 静态数组结构存储方法
此类方法采用静态数组存储稀疏矩阵非零元素的数值和位置信息，常用的方法有散居格式、按行（列）存储格式、三角检索存储格式等。
3.1.1 散居格式
在散居格式中，对n m ?阶稀疏矩阵A ，其非零元素共有τ个，令ij a 是A 中第i 行第j 列非零元素。可以定义三个数组，按下面的存储格式存储矩阵A 中非零元素的信息：
V A ——存储A 中非零元素ij a 的值，共τ个 JA ——存储A 中非零元素ij a 的行指标i ，共τ个 IA ——存储A 中非零元素ij a 的列指标j ，共τ个
因此，散居格式存储稀疏矩阵时，共需要3τ个存储单元。在散居格式中，矩阵A 中非零元在数组中的位置可任意排列，修改灵活。但是，其存储顺序无一定规律，检索起来不方便，最差的情况下可能性要在整个数组中查找一遍。
3.1.2 按行（列）存储格式
按行（列）存储格式的原理是按行(列)顺序依次存储n n ?阶矩阵A 中的非零元，同一行(列)元素依次排在一起。以按行为例，其存储格式是：
V A ——按行存储矩阵A 中非零元ij a 的值，共m 个； JA ——按行存储矩阵A 中非零元ij a 的列号，共m 个； IA ——记录A 中每行第一个非零元素在VA 中的位置，共n 个。
由按行（列）存储格式的特点，要查找第i 行的非零元素,在V A 中取出从
IA(i)=k 到1－1)＋IA(i =k 共IA(i)－1)＋IA(i 个非零元就是A 中第i 行的全部非零
元,非零元的值是V A(k)，列号为JA(k)。要查找第i 行第j 列的元素ij a 在V A 中的位置，则对k 从IA(i)到1－1)＋IA(i ，判断列号JA(k)是否等于j ，如果相等，那么V A(k)就是要查找的非零元ij a
3.1.3 三角检索存储格式
用三角检索存储格式来存储非零元素的方法，特别适用于稀疏矩阵的三角分解。以按行、列存储n n ?阶矩阵A 的上下三角部分非零元素为例来说明三角检索存储格式的原理。 其存储格式是：
U ——存储矩阵A 的上三角部分的非零元素的值，按行依次存储 JU ——存储矩阵A 的上三角部分的非零元素的列号
IU ——存储矩阵A 中上三角部分每行第一个非零元素在U 中的位置(首地址) L ——存储矩阵A 中下三角部分的非零元素的值，按列依次存储 IL ——存储矩阵A 的下三角部分的非零元素的行号
JL ——存储矩阵A 中下三角部分每列第一个非零元素在L 中的位置(首地址) D ——存储矩阵A 的对角元素的值，其检索下标不需要存储
三角检索存储格式示例：
????
???
???=4443
42332322211412110
00000a a a a a a a a a a A
则矩阵A 的存储可表示为：
[]
[][][]
[]
[][]
4433
12114342
21
2314
12
43214424431342a a a a D JL IL a a a L IU JU a a a U =======
IU （3）为4，表明矩阵A 上三角部分第3行的第1个非零元如果有的话应在U 的第4个位置，而U 表中第4个位置没有非零元素，为了检索方便，IU(3)仍应赋值4。有了IU 表即可知道A 的上三角部分第i 行的非零元的数目。
如果要查找矩阵A 的上三角第i 行所有非零元素，只要取出从[]IU(i)U 到
[]1-1)IU(i U +共IU(i)-1)IU(i +个元素即可，JU(k)指出了该元素的列号，U(k)是该非零元素的值。对于按列存储的格式进行查找的情况类同。