一种适合研究目的使用的牛顿法潮流计算方法 技术领域 本发明涉及一种电力系统牛顿法潮流计算方法， 特别是一种适合研究目的使用的 牛顿法潮流计算方法。
背景技术
潮流计算是电力系统中最基本的一种计算， 也是其它电力系统分析 ( 如静态安全 分析等 ) 的基础。牛顿法潮流计算算法是一种最常用的潮流计算方法， 科研人员经常以牛 顿法潮流计算为基础进行进一步的研究。 实用算法采用稀疏矩阵技术和节点优化编号等高 级技术。 这些技术虽然能大幅度提高潮流计算的速度、 降低内存占用量， 但编程非常麻烦且 不易修改和维护。 为了研究目的， 迫切需要一种易于编程、 修改和调试的牛顿法潮流计算算 法。 牛顿法潮流计算中， 如果节点电压采用极坐标表示为 则为极坐标牛
顿法潮流计算， 其主要方程如下 ：
导纳矩阵为 ：
导纳矩阵可以反映出节点间的连接关系， 如果两个节点不直接连接， 它们的互导 纳为 0。电力网络中每个节点只与几个节点相连， 因此导纳矩阵是一个高度稀疏的矩阵， 它 的大部分元素为 0。
功率偏差方程为 ：
式中， j ∈ i 表示节点 j 与节点 i 直接相连， n 为节点数， m 为 PQ 节点数。 修正方程为 ：
式中，为雅可比矩阵， 其元素计算公式为 ： 当j≠i时
Hij、 Nij、 Mij、 Lij 分别对应雅可比矩阵的 Ji+i-1， Ji+i-1， Ji+i， Ji+i， j+j-1、 j+j、 j+j-1、 j+j。当 j ＝i时
或用下列式子计算 ：式中， Pi、 Qi 分别为节点 i 的计算有功功率和无功功率， 按下式计算 ：潮流计算的修正方程为线性方程， 设线性方程一般形式如下 ：
AX ＝ B (17)
对于线性方程， 可以采用高斯消去法求解， 本算法采用按行消去， 按行回代， 包括 3 个步骤， 分别为 ：
(a) 对系数矩阵的消去和规格化， 公式如下 ：
(b) 对右端常数项的消去和规格化， 公式如下 ：(c) 回代求解， 公式如下 ：极坐标牛顿法潮流计算算法的原理流程图如附图 1 所示， 主要包括以下步骤 ：
步骤 1 ： 原始数据输入和电压初始化 ：
电压初始化一般采用平启动， 即 PV 节点和平衡节点的电压幅值取给定值， PQ 节点 的电压幅值取 1.0 ； 所有电压的相角都取 0.0。这里单位采用标幺值。
步骤 2 ： 形成节点导纳矩阵 Y ；
步骤 3 ： 形成雅可比矩阵 J ：
步骤 4 ： 解方程及修正电压幅值 V、 相角 θ ；
电压修正公式为 ：
步骤 5 ： 节点及支路数据输出。发明内容 为满足研究目的需要， 本发明要提出一种适合研究目的使用的牛顿法潮流计算方 法， 为以牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的牛 顿法潮流计算算法。
为了实现上述目的， 本发明采用极坐标牛顿法潮流计算的基本原理步骤， 在分析 潮流计算各个组成部分的基础上提出了一种通过简单逻辑判断来避免不必要运算以提高 潮流计算计算速度的方法。 本发明的技术方案是在形成雅可比矩阵模块和修正方程高斯消 去模块中通过简单判断来避免不必要的运算。
从式 (3) ～式 (6) 可以看出， 当节点 i 和节点 j 之间没有支路连接时， 它们之间的 互导纳为 0， 不需要计算雅可比矩阵元素， 因此可以通过判断语句避免不必要的计算。同时 雅可比矩阵按 (2n)×(2n) 存储， PV 节点 ΔQ 对应的行和 ΔV 对应的列的元素都为 0， 与平 衡节点有关的行和列的元素也都为 0， 这样虽然增加了内存需求量， 但可以简化节点和方程 行列的对应关系， 大大降低编程难度， 也不会增加计算量。 形成雅可比矩阵模块包括以下步 骤：
步骤 1 ： 设置行号 i ＝ 1。
步骤 2 ： 判断节点 i 是否为平衡节点， 如果是平衡节点转至步骤 15。
步骤 3 ： 设置列号 j ＝ 1。
步骤 4 ： 判断节点 j 是否为平衡节点， 如果是平衡节点转至步骤 13。
步骤 5 ： 判断节点 i 和节点 j 之间的导纳实部 Gij 和虚部 Bij 是否都为 0， 如果都为 0 转至步骤 13。
步骤 6 ： 计算雅可比矩阵元素 Ji+i-1， 如果 j ≠ i， 根据式 (3) 计算， 如果 j ＝ i， j+j-1， 根据式 (11) 计算。
步骤 7 ： 判断节点 j 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 9。
步骤 8 ： 计算雅可比矩阵元素 Ji+i-1，j+j， 如果 j ≠ i， 根据式 (4) 计算， 如果 j ＝ i， 根据式 (12) 计算。
步骤 9 ： 判断节点 i 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 13。
步骤 10 ： 计算雅可比矩阵元素 Ji+i， 如果 j ≠ i， 根据式 (5) 计算， 如果 j ＝ i， j+j-1， 根据式 (13) 计算。
步骤 11 ： 判断节点 j 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 13。
步骤 12 ： 计算雅可比矩阵元素 Ji+i， 如果 j ≠ i， 根据式 (6) 计算， 如果 j ＝ i， 根 j+j， 据式 (14) 计算。
步骤 13 ： 令 j ＝ j+1。
步骤 14 ： 判断 j 是否大于 n， 如果 j 不大于 n， 则返回到步骤 4。
步骤 15 ： 令 i ＝ i+1。
步骤 16 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
从式 (18) 可以看出， 如果则第 k 行对第 i 行的消去就没有必要了， 可以通过判断语句避免不必要的计算 ； 同时通过逻辑判断把雅可比矩阵中主对角元素为 0 的行 ( 主对角元素为 0 表示该行元素都是 0) 跳过去。本发明修正方程高斯消去模块的步骤是 ：
步骤 1 ： 设置当前行号 i ＝ 1。
步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。
步骤 3 ： 设置 k ＝ 1。
步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。
步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。
步骤 6 ： 设置当前列号 j ＝ k+1。
步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。
步骤 8 ： 根据式 (18) 对系数矩阵进行消去运算。
步骤 9 ： 令 j ＝ j+1， 返回到步骤 7。
步骤 10 ： 根据式 (20) 对右端常数项进行消去运算。
步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。
步骤 12 ： 设置当前列号 j ＝ i+1。
步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。
步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。
步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。
步骤 16 ： 根据式 (21) 对右端常数项进行规格化运算。
步骤 17 ： 令 i ＝ i+1。
步骤 18 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
与现有技术相比， 本发明具有以下有益效果 ：
为以牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维 护的牛顿法潮流计算算法。采用本发明对一个 445 节点实际电网进行了计算， 计算时间为 0.484s， 计算速度完全能够满足科研需要。其潮流计算的内存占用量在 10MB 左右， 目前计 算机的内存一般都有几个 GB， 完全能够满足要求。附图说明
本发明共有附图 6 张， 其中 ： 图 1 是极坐标牛顿法潮流计算原理流程框图。 图 2 是本发明雅可比矩阵形成的流程框图。 图 3 是本发明修正方程高斯消去模块流程框图。 图 4 是本发明方法实施例的网络图。 图 5 是输电线路的等值电路。 图 6 是变压器支路的等值电路。具体实施方式
下面结合附图对本发明作进一步地说明。附图 2 和附图 3 中 ε 为给定的小正数， 用来判断浮点数是否为 0， 本发明算例中令 ε ＝ 10-20， 图 4 中编号为母线号。
附图 4 是一个简单电力系统网络， 包括 2 个发电机， 3 个负荷， 2 条变压器支路， 2条 输电线路。共 5 条母线 ( 节点 )， 其中母线 1 为 PV 节点， 此母线的有功功率和电压幅值给 定， 母线 2 为平衡节点， 此母线的电压幅值和相角给定， 其它 3 条母线为 PQ 节点， 这些母线 的有功功率和无功功率给定。支路参数见表 1， 母线数据见表 2， 表中单位为标幺值。
表 1 附图 4 实施例的支路参数
表 2 附图 4 实施例的母线数据
采用本发明方法对附图 4 所示网络进行潮流计算， 步骤如下 ：
步骤 1 ： 原始数据输入和电压初始化。
输入表 1 和表 2 的原始数据。电压初始化采用平启动， 母线 1 和母线 2 电压初值 取给定电压， 都为 1.05， 其它母线电压初值取 1.0。平衡节点的相角取 0.0。
步骤 2 ： 形成节点导纳矩阵 Y。导纳矩阵如下 ：
步骤 3 ： 形成雅可比矩阵 J。第 1 次迭代过程的雅可比矩阵为 ：
步骤 4 ： 解方程及修正电压幅值 V、 相角 θ。 电压修正公式为 ：第 1 次迭代过程中， 解方程得到的修正量及电压幅值 V、 相角 θ 见表 3。算例共迭 代 3 次， 各次迭代结果见表 4。
表 3 附图 4 实施例的第 1 次迭代结果
表 4 附图 4 实施例的各次迭代结果V1/θ1 1.050/8.34° 1.050/7.10° 1.050/7.07° V2/θ2 1.050/0.00° 1.050/0.00° 1.050/0.00° V3/θ3 1.053/-1.82° 1.040/-2.04° 1.040/-2.05° V4/θ4 1.054/-3.95° 1.044/-3.93° 1.044/-3.94° V5/θ5 1.093/4.04° 1.079/3.11° 1.079/3.08°迭代序号 1 2 3步骤 5 ： 节点及支路数据输出。
节点电压和功率计算结果见表 5。支路传输功率计算结果见表 6， 由于支路存在网 损， 支路两端的功率不是互为相反数。表中单位除相角外都为标幺值。
表 5 附图 4 实施例的节点电压和功率计算结果
母线号 1 2 3 4 5
电压幅值 V 1.05000 1.05000 1.03994 1.04416 1.07876电压相角 θ(° ) 7.07065 0.00000 -2.04660 -3.93767 3.08399有功 P 5.00000 2.39012 -1.60000 -3.70000 -2.00000无功 Q 1.75685 2.02699 -0.80000 -1.30000 -1.00000表 6 附图 4 实施例的支路传输功率计算结果 母线号 i 1 2 3 3 母线号 j 5 4 4 5 Pij 5.00000 2.39012 1.32078 -2.92078 Qij 1.75686 2.02699 -0.52292 -0.27708 Pji -5.00000 -2.39012 -1.30988 3.00000 Qji -1.37473 -1.75974 0.45974 0.37473支路号 1 2 3 4
按照图 2 所示流程， 第 1 次迭代时生成雅可比矩阵的步骤如下 ( 形成雅可比矩阵 前， 此矩阵元素全部清零 ) ：
步骤 1 ： 设置行号 i ＝ 1。
形成与第 1 个节点有关的雅可比矩阵第 1 行和第 2 行元素。虽然节点 1 是 PV 节 点， 雅可比矩阵不包含与该节点无功有关的行， 也不包含与该节点电压幅值有关的列。 但为 了计算和编程方便， 仍保留该行和列， 其元素都为 0。
步骤 2 ： 判断节点 i 是否为平衡节点， 如果是平衡节点转至步骤 15。
节点 1 是 PV 节点， 不是平衡节点， 执行步骤 3。
步骤 3 ： 设置列号 j ＝ 1。形成雅可比矩阵的第 1 列和第 2 列元素。
步骤 4 ： 判断节点 j 是否为平衡节点， 如果是平衡节点转至步骤 13。 节点 1 是 PV 节点， 不是平衡节点， 执行步骤 5。 步骤 5 ： 判断节点 i 和节点 j 之间的导纳实部 Gij 和虚部 Bij 是否都为 0， 如果都为0 转至步骤 13。
由于导纳矩阵虚部 B11 ＝ -66.67 不为 0， 执行步骤 6。
步骤 6 ： 计算雅可比矩阵元素 Ji+i-1， 如果 j ≠ i， 根据式 (3) 计算， 如果 j ＝ i， j+j-1，
根据式 (11) 计算。
由于 j ＝ i， 根据式 (11) 计算雅可比矩阵元素 J11 ＝ -66.67。
步骤 7 ： 判断节点 j 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 9。
节点 1 不是 PQ 节点， 转至步骤 9。
步骤 9 ： 判断节点 i 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 13。
节点 1 不是 PQ 节点， 转至步骤 13。
步骤 13 ： 令 j ＝ j+1。j ＝ 1+1 ＝ 2。
步骤 14 ： 判断 j 是否大于 n， 如果 j 不大于 n， 则返回到步骤 4。
n ＝ 5， j ＝ 2 不大于 n， 返回到步骤 4。
步骤 4 ： 判断节点 j 是否为平衡节点， 如果是平衡节点转至步骤 13。
节点 2 是平衡节点， 转至步骤 13。
步骤 13 ： 令 j ＝ j+1。j ＝ 2+1 ＝ 3。 步骤 14 ： 判断 j 是否大于 n， 如果 j 不大于 n， 则返回到步骤 4。
n ＝ 5， j ＝ 3 不大于 n， 返回到步骤 4。
步骤 4 ： 判断节点 j 是否为平衡节点， 如果是平衡节点转至步骤 13。
节点 3 是 PQ 节点， 不是平衡节点， 执行步骤 5。
步骤 5 ： 判断节点 i 和节点 j 之间的导纳实部 Gij 和虚部 Bij 是否都为 0， 如果都为 0 转至步骤 13。
由于导纳矩阵实部 G13 和虚部 B13 都为 0， 转至步骤 13。
步骤 13 ： 令 j ＝ j+1。j ＝ 3+1 ＝ 4。
步骤 14 ： 判断 j 是否大于 n， 如果 j 不大于 n， 则返回到步骤 4。
n ＝ 5， j ＝ 4 不大于 n， 返回到步骤 4。
步骤 4 ： 判断节点 j 是否为平衡节点， 如果是平衡节点转至步骤 13。
节点 4 是 PQ 节点， 不是平衡节点， 执行步骤 5。
步骤 5 ： 判断节点 i 和节点 j 之间的导纳实部 Gij 和虚部 Bij 是否都为 0， 如果都为 0 转至步骤 13。
由于导纳矩阵实部 G14 和虚部 B14 都为 0， 转至步骤 13。
步骤 13 ： 令 j ＝ j+1。j ＝ 4+1 ＝ 5。
步骤 14 ： 判断 j 是否大于 n， 如果 j 不大于 n， 则返回到步骤 4。
n ＝ 5， j ＝ 5 不大于 n， 返回到步骤 4。
步骤 4 ： 判断节点 j 是否为平衡节点， 如果是平衡节点转至步骤 13。
节点 5 是 PQ 节点， 不是平衡节点， 执行步骤 5。
步骤 5 ： 判断节点 i 和节点 j 之间的导纳实部 Gij 和虚部 Bij 是否都为 0， 如果都为 0 转至步骤 13。
由于导纳矩阵虚部 B15 ＝ 63.49 不为 0， 执行步骤 6。
步骤 6 ： 计算雅可比矩阵元素 Ji+i-1， 如果 j ≠ i， 根据式 (3) 计算， 如果 j ＝ i， j+j-1，
根据式 (11) 计算。
由于 j ≠ i， 根据式 (3) 计算雅可比矩阵元素 J19 ＝ 66.67。
步骤 7 ： 判断节点 j 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 9。
节点 5 是 PQ 节点， 执行步骤 8。
步骤 8 ： 计算雅可比矩阵元素 Ji+i-1，j+j， 如果 j ≠ i， 根据式 (4) 计算， 如果 j ＝ i， 根据式 (12) 计算。
由于 j ≠ i， 根据式 (4) 计算雅可比矩阵元素 J1， 10 ＝ 0.000。
步骤 9 ： 判断节点 i 是否为 PQ 节点， 如果不是 PQ 节点转至步骤 13。
节点 1 是 PV 节点， 不是 PQ 节点， 转至步骤 13。
步骤 13 ： 令 j ＝ j+1。j ＝ 5+1 ＝ 6。
步骤 14 ： 判断 j 是否大于 n， 如果 j 不大于 n， 则返回到步骤 4。
n ＝ 5， j ＝ 6 大于 n， 执行步骤 15。
至此与节点 1 有关的雅可比矩阵前两行元素已经计算完毕， 结果见式 (24) 的前两 行。
步骤 15 ： 令 i ＝ i+1。i ＝ 1+1 ＝ 2。 形成与第 2 个节点有关的雅可比矩阵第 3 行和第 4 行元素。
步骤 16 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
n ＝ 5， i ＝ 2 不大于 n， 返回到步骤 2。
步骤 2 ： 判断节点 i 是否为平衡节点， 如果是平衡节点转至步骤 15。
节点 2 是平衡节点， 转至步骤 15。
步骤 15 ： 令 i ＝ i+1。i ＝ 2+1 ＝ 3。
形成与第 3 个节点有关的雅可比矩阵第 5 行和第 6 行元素。
步骤 16 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
n ＝ 5， i ＝ 3 不大于 n， 返回到步骤 2。重复此过程直到 i ＝ 6 大于 n， 结束。形成 的导纳矩阵见式 (24)。
第 1 次迭代时， 得到的潮流计算修正方程如下 ：
按照图 3 所示流程， 第 1 次迭代时高斯消去的步骤如下 ： 步骤 1 ： 设置当前行号 i ＝ 1。对第 1 行进行消去操作。 步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。a11 ＝ -66.67 不为 0， 执行步骤 3。
步骤 3 ： 设置 k ＝ 1。用第 1 行消去。
步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。
k ＝ 1， i ＝ 1， k 不小于 i， 转至步骤 12。
步骤 12 ： 设置当前列号 j ＝ i+1。j ＝ 1+1 ＝ 2。
步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。
n ＝ 10， j ＝ 2 不大于 n， 执行步骤 14。
步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。
结果为 a12 ＝ 0.000/(-66.67) ＝ 0.000。
步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 2+1 ＝ 3， 返回到步骤 13。
步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。
n ＝ 10， j ＝ 3 不大于 n， 执行步骤 14。
步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。
结果为 a13 ＝ 0.000/(-66.67) ＝ 0.000。重复此过程直到 j ＝ 10， 完成系数矩阵 第 1 行规格化运算， 结果除 a11 ＝ -66.67 不变， a19 ＝ -1.000 外， 其它都为 0.000。
步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 10+1 ＝ 11， 返回到步骤 13。
步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。
n ＝ 10， j ＝ 11 大于 n， 转至步骤 16。
步骤 16 ： 根据式 (21) 对右端常数项进行规格化运算。
结果为 b1 ＝ 5.000/(-66.67) ＝ -0.075。
步骤 17 ： 令 i ＝ i+1。i ＝ 1+1 ＝ 2， 对第 2 行进行消去操作。
步骤 18 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
n ＝ 10， i ＝ 2 不大于 n， 返回到步骤 2。
步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。
a22 ＝ 0， 转至步骤 17。
步骤 17 ： 令 i ＝ i+1。i ＝ 2+1 ＝ 3， 对第 3 行进行消去操作。
步骤 18 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
n ＝ 10， i ＝ 3 不大于 n， 返回到步骤 2。
步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。
a33 ＝ 0， 转至步骤 17。
步骤 17 ： 令 i ＝ i+1。i ＝ 3+1 ＝ 4， 对第 4 行进行消去操作。
步骤 18 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
n ＝ 10， i ＝ 4 不大于 n， 返回到步骤 2。
步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。
a44 ＝ 0， 转至步骤 17。
步骤 17 ： 令 i ＝ i+1。i ＝ 4+1 ＝ 5， 对第 5 行进行消去操作。
步骤 18 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。
n ＝ 10， i ＝ 5 不大于 n， 返回到步骤 2。
步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。a55 ＝ -64.24 不为 0， 执行步骤 3。 步骤 3 ： 设置 k ＝ 1。用第 1 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 1， i ＝ 5， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a51 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 1+1 ＝ 2， 用第 2 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 2， i ＝ 5， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a52 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 2+1 ＝ 3， 用第 3 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 3， i ＝ 5， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a53 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 3+1 ＝ 4， 用第 4 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 4， i ＝ 5， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a54 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 4+1 ＝ 5， 用第 5 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 5， i ＝ 5， k 不小于 i， 转至步骤 12。 步骤 12 ： 设置当前列号 j ＝ i+1。j ＝ 5+1 ＝ 6。 步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。 n ＝ 10， j ＝ 6 不大于 n， 执行步骤 14。 步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。 结果为 a56 ＝ -16.62/(-64.24) ＝ 0.259。 步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 6+1 ＝ 7， 返回到步骤 13。 步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。 n ＝ 10， j ＝ 7 不大于 n， 执行步骤 14。 步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。 结果为 a57 ＝ 37.82/(-64.24) ＝ -0.589。 步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 7+1 ＝ 8， 返回到步骤 13。 步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。 n ＝ 10， j ＝ 8 不大于 n， 执行步骤 14。 步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。 结果为 a58 ＝ 9.077/(-64.24) ＝ -0.141。步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 8+1 ＝ 9， 返回到步骤 13。 步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。 n ＝ 10， j ＝ 9 不大于 n， 执行步骤 14。 步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。 结果为 a59 ＝ 26.42/(-64.24) ＝ -0.411。 步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 9+1 ＝ 10， 返回到步骤 13。 步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。 n ＝ 10， j ＝ 10 不大于 n， 执行步骤 14。 步骤 14 ： 根据式 (19) 对系数矩阵进行规格化运算。 结果为 a5， 10 ＝ 7.547/(-64.24) ＝ -0.117。 步骤 15 ： 令 j ＝ j+1， 返回到步骤 13。j ＝ 10+1 ＝ 11， 返回到步骤 13。 步骤 13 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 16。 n ＝ 10， j ＝ 11 大于 n， 转至步骤 16。 步骤 16 ： 根据式 (21) 对右端常数项进行规格化运算。 结果为 b5 ＝ -1.600/(-64.24) ＝ 0.025。 步骤 17 ： 令 i ＝ i+1。i ＝ 5+1 ＝ 6， 对第 6 行进行消去操作。 步骤 18 ： 判断 i 是否大于 n， 如果 i 不大于 n， 则返回到步骤 2 ； 否则结束。 n ＝ 10， i ＝ 6 不大于 n， 返回到步骤 2。 步骤 2 ： 判断矩阵元素 aii 是否为 0， 如果 aii 为 0， 则转至步骤 17。 a66 ＝ -63.98 不为 0， 执行步骤 3。 步骤 3 ： 设置 k ＝ 1。用第 1 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 1， i ＝ 6， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a61 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 1+1 ＝ 2， 用第 2 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 2， i ＝ 6， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a62 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 2+1 ＝ 3， 用第 3 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 3， i ＝ 6， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a63 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 3+1 ＝ 4， 用第 4 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 4， i ＝ 6， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。a64 ＝ 0， 转至步骤 11。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 4+1 ＝ 5， 用第 5 行消去。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 5， i ＝ 6， k 小于 i， 执行步骤 5。 步骤 5 ： 判断矩阵元素 aik 是否为 0， 如果 aik 为 0， 则转至步骤 11。 a65 ＝ 16.62， 执行步骤 6。 步骤 6 ： 设置当前列号 j ＝ k+1。j ＝ 5+1 ＝ 6。 步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。 n ＝ 10， j ＝ 6， j 不大于 n， 执行步骤 8。 步骤 8 ： 根据式 (18) 进行消去运算。 a66 ＝ a66-a65a56 ＝ -63.98-16.62×0.259 ＝ -68.28。 步骤 9 ： 令 j ＝ j+1， 返回到步骤 7。j ＝ 6+1 ＝ 7， 返回到步骤 7 步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。 n ＝ 10， j ＝ 7， j 不大于 n， 执行步骤 8。 步骤 8 ： 根据式 (18) 对系数矩阵进行消去运算。 a67 ＝ a67-a65a57 ＝ -9.077-16.62×(-0.589) ＝ 0.712。 步骤 9 ： 令 j ＝ j+1， 返回到步骤 7。j ＝ 7+1 ＝ 8， 返回到步骤 7 步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。 n ＝ 10， j ＝ 8， j 不大于 n， 执行步骤 8。 步骤 8 ： 根据式 (18) 对系数矩阵进行消去运算。 a68 ＝ a68-a65a58 ＝ 37.82-16.62×(-0.141) ＝ 40.16。 步骤 9 ： 令 j ＝ j+1， 返回到步骤 7。j ＝ 8+1 ＝ 9， 返回到步骤 7 步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。 n ＝ 10， j ＝ 9， j 不大于 n， 执行步骤 8。 步骤 8 ： 根据式 (18) 对系数矩阵进行消去运算。 a69 ＝ a69-a65a59 ＝ -7.547-16.62×(-0.411) ＝ -0.716。 步骤 9 ： 令 j ＝ j+1， 返回到步骤 7。j ＝ 9+1 ＝ 10， 返回到步骤 7 步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。 n ＝ 10， j ＝ 10， j 不大于 n， 执行步骤 8。 步骤 8 ： 根据式 (18) 对系数矩阵进行消去运算。 a6， 10 ＝ a6， 10-a65a5， 10 ＝ 26.42-16.62×(-0.117) ＝ 28.36。 步骤 9 ： 令 j ＝ j+1， 返回到步骤 7。j ＝ 10+1 ＝ 11， 返回到步骤 7 步骤 7 ： 判断 j 是否大于 n， 如果 j 大于 n， 则转至步骤 10。 n ＝ 10， j ＝ 11， j 大于 n， 转至步骤 10。 步骤 10 ： 根据式 (20) 进行对右端常数项进行消去运算。 b6 ＝ b6-a65b5 ＝ -0.670-16.62×0.025 ＝ -1.086。 步骤 11 ： 令 k ＝ k+1， 返回到步骤 4。k ＝ 5+1 ＝ 6， 返回到步骤 4。 步骤 4 ： 判断 k 是否小于 i， 如果 k 不小于 i， 则转至步骤 12。 k ＝ 6， i ＝ 6， k 不小于 i， 转至步骤 12， 对第 6 行进行规格化。重复此过程直到 i＝ 11 大于 n， 结束， 完成对方程的高斯消去工作。得到方程如下 ：
采用本发明和几种对比方法对另一个实施例进行了计算， 该实施例为一个 445 节 点实际电网， 节点数为 445， 支路数为 544。计算时收敛精度为 0.00001。几种潮流计算算法 分别为 ：
方法 1 ： 不采用稀疏矩阵技术。
方法 2 ： 不采用稀疏矩阵技术， 但使用判断避免不必要运算， 为本发明方法。
方法 3 ： 采用稀疏矩阵技术， 但不使用节点优化编号。
方法 4 ： 采用稀疏矩阵技术和节点优化编号。
几种方法的计算时间见表 7， 从表 7 可见仅仅通过简单的逻辑判断就能极大地提 高计算速度， 计算速度与采用稀疏矩阵技术的方法 3 和方法 4 相差不太大， 但编程难度要小 得多， 也利于修改和维护。
表 7 几种牛顿法计算时间比较
潮流计算方法 方法 1 方法 2 方法 3 方法 4计算时间 (s) 5.953 0.484 0.296 0.156不采用稀疏矩阵技术的方法 1 和方法 2 的内存需求量见表 8， 表中矩阵元素为浮点 数据类型。 采用稀疏矩阵技术的方法 3 和方法 4 的内存需求量见表 9， 表中矩阵元素为浮点 数据类型和整型数。 表中导纳矩阵元素数为导纳矩阵元素个数， 导纳矩阵元素为复数， 因此 每个导纳矩阵元素为两个浮点数据。每个双精度浮点数占 8 个字节， 每个整型数占 4 个字 节。
表 8 方法 1 和方法 2 的内存需求量
表 9 方法 3 和方法 4 的内存需求量表 9 中导纳矩阵、 雅可比矩阵和方程系数矩阵采用稀疏矩阵技术存储， 稀疏矩阵 技术一般采用按行存储方式。设一个矩阵行数为 n， 非零元素个数为 m， 需要存储以下非零 元素信息 ：
VA——存储矩阵中非零元素的值， 共 m 个， 为浮点数据类型 ；
JA——存储矩阵中非零元素的列号， 共 m 个， 为整型 ；
IA——存储矩阵中每行第一个非零元素在 VA 中的位置， 共 n 个， 为整型。
内存需求量按下式计算 ：
S ＝ 8m+4m+4n (25)
导纳矩阵元素为复数， 需要两个浮点数单元， 采用稀疏矩阵技术存储时按下式计 算：
S ＝ 16m+4m+4n (26)
由表 8 和表 9 的结果可见， 采用稀疏矩阵技术， 内存占用量明显减少， 同时采用稀 疏矩阵技术和节点优化编号算法的内存占用量更少。 虽然不采用稀疏矩阵技术内存占用量 明显增加， 加上其它数组内存的总内存占用量约为 10MB。考虑到目前计算机一般有几个 GB 的内存， 对于多达 445 节点的大型实际电网来说， 10MB 左右的内存要求并不太多。
本算法可以采用任何一种编程语言和编程环境实现， 如 C 语言、 C++、 FORTRAN、 Delphi 等。开发环境可以采用 Visual C++、 Borland C++Builder、 Visual FORTRAN 等。