CN105354422B - 一种基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法潮流的方法 - Google Patents

Abstract

一种基于对称稀疏矩阵技术快速求取极坐标牛顿‑拉夫逊法潮流的方法，属于电力系统分析计算领域。包括：快速读入仅含非零元素的导纳矩阵Y数据文件；根据Y阵和雅可比矩阵J结构相似的特点快速形成J阵；利用J阵元素的对称稀疏性对J阵快速进行消元和回代求取潮流。本发明无论在读取数据文件、形成J阵、对J阵进行消元和回代计算等方面的计算速度均远快于传统方法中不考虑或考虑元素稀疏性的计算速度。用本发明对各IEEE系统进行验算完全可行。如对IEEE‑118系统，读取数据文件的时间减少约82％，形成J阵和对J阵进行消元和回代的时间减少约66～90％，整个潮流计算的时间减少约80％～85％。且系统节点数越多，本发明的优势越大。

Description

一种基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法 潮流的方法

技术领域

本发明属于电力系统分析计算领域。

背景技术

在大型电力系统的潮流计算中，经常需根据Y阵形成J阵并对其进行消元和回代，而Y、J阵中有大量的零元素且元素结构非常相似。如不考虑利用元素的稀疏性及元素结构的特点，则在J阵的形成过程中需读取较大的Y阵数据文件并进行大量无效J阵元素的计算。这不但使J阵形成时间过长，也会使J阵在消元和回代过程中由于大量零元素的无效计算而大大降低计算效率。如利用了元素的稀疏性，但未利用其元素结构的各个特点，计算效率仍然低下。

传统的极坐标牛顿-拉夫逊法潮流在J阵的形成、J阵的消元和回代计算过程中主要有以下几点不足：

(1)形成J阵所用的Y阵数组形式不合适。

形成J阵所用的Y阵数组为Y(n,2n)，大量零元素的存贮不但需要较大的存贮单元，且数据文件的读取速度也不高。即使考虑Y阵元素的稀疏性采用按坐标存贮、或按顺序存贮、或按链表存贮可省去不少贮存单元，但Y、J阵元素之间没有明显的对应关系，不但数据文件的读取时间较长，也不便于直接计算J阵元素。且由于读取Y阵数据文件的时间远远高于用Y阵形成J阵、对J阵进行消元和回代的时间，因此快速形成J阵和潮流计算的关键在于Y阵元素的存贮和读取方式。

(2)计算I_pi、I_qi或ΔP_i、ΔQ_i过程的不合适，一般存在对大量零元素的无效计算。

(3)形成J阵的方式不合适。

Y、J阵元素的稀疏性及结构极为相似。如果直接根据Y(n,2n)形成J阵，则由于大量零元素的计算导致形成J阵的计算效率极低；如果判断Y(n,2n)中非零元素的方式不合适，如同时判断Y(n,2n)中非零元素的实部和虚部，则相当大量的判断语句仍然导致形成J阵的计算效率极低，即使仅判断Y(n,2n)中非零元素的虚部，同样由于大量的判断语句仍然导致形成J阵的计算效率较低；如果不能按两行/次、或(两行+两列)/次计算J阵元素；或未利用Y_ij≠0时，有J_ij≠0和J_ji≠0的关系等，则形成J阵的计算效率仍然非常低。

(4)消元过程中未利用J阵元素的对称性和稀疏性或利用方式不合适。

如未利用J阵元素的稀疏性；未利用子阵J_ij与J_ji的对称性；未利用构成元素个数不同的子阵J_ij中H_ij≠0时，有N_ij≠0、M_ij≠0、L_ij≠0以及H_ji≠0、N_ji≠0、M_ji≠0、L_ji≠0或M_ij≠0以及H_ji≠0、N_ji≠0或H_ji≠0的关系；未利用对角元以右和以下非零元素位置对称的关系；未利用2种不同非零元素的交叉点来确定相应的计算元素；未利用J阵元素的结构特点按1次判断可对两行/次非零元素规格化、对两列/次非零元素消元；未利用元素的对称性仅计算对角元以右的元素；未利用多次迭代过程中J阵非零元素位置不变的特性避免反复进行非零判断等，都会导致计算效率十分低下或非常低下。

(5)回代过程中未利用消元过程中所记忆的上三角元素的稀疏性，同样导致计算效率低下。

发明内容

为了克服上述现有技术的不足，本发明在快速存贮及读取电力系统稀疏矩阵数据以及快速形成电力系统潮流计算中J阵的基础上，提出一种基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法潮流的方法。

在建立J阵之前，先建立无非零元素的Y阵数据文件(详见中国专利申请2015103648275)。将Y阵的节点分为主节点和子节点。假设系统中各节点的最大连接支路数为l_max，定义Y阵数据文件的数组为Y(n,d)，其中d＝3p+2，p＝l_max+1。分为3组，第1组1列，存贮该行主节点的行号i；第2组1列，存贮与行号对应的节点数S_i，为主节点及与该主节点连接的子节点数之和，S_i值由程序自动累加以保证快速读取对应主节点和子节点的参数，以免对数组Y(n,d)中各行多余存贮单元的读取；第3组共d-2列，按递增顺序存贮主节点和与该主节点相连接的所有非零元素子节点的列号j及参数的实部g_ij和虚部b_ij，分别位于第3～d列，存贮方式如下。

注：第三行数字1～d表示数组Y(n,d)中的实际列号，而同一列号下的所有各行的列号如j_p的数值是不同的，且并非每行都有j_p及其相应参数，即并非每行都能存满相应的列号和参数。

Y(n,d)数组的存贮方式去掉了所有导纳矩阵中实部与虚部均为零的零元素，使所有的非零元素按列号顺序存放以便后续应用。在建立Y(n,d)数组的数据文件基础上，再利用Y阵与J阵元素结构相似、非零元素对应、以及J_ij子阵同组元素之间的关系等特点，在形成J阵时可以免去对大量Y阵元素的非零判断以及对J阵元素的无效计算。再应用本发明提出的基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法潮流的方法，在消元和回代过程中充分利用J阵元素的对称性和稀疏性，可大大加快牛顿法潮流的整体计算速度。

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

步骤1：建立无非零元素的Y阵数据文件；

步骤2：打开数据文件，读取仅含非零元素的Y阵数据文件到Y(n,d)数组；

按Y(n,d)数组方式形成Y阵和建立Y阵数据文件，以便在新的潮流计算等程序中可直接打开并读取Y(n,d)数组的数据文件。与打开并读取Y(n,2n)数组的数据文件相比，可节省大量的存贮单元和读写时间。Y(n,d)数组中的节点数S_i能进一步提高对Y(n,d)数组元素的读取效率。

步骤3：根据Y(n,d)数组中的每行元素，计算J阵中的非零元素；

(1)用Y(n,d)数组可直接计算各个节点的I_pi、I_qi或ΔP_i、ΔQ_i，省去对所有零元素的判断或无效计算。

(2)假设系统的节点数为n，PQ节点数为m，m+1及其后的节点均为PV节点，第n个节点是平衡节点。J阵元素排列和对应的修正方程式如下：

(3)根据Y(n,d)数组中的每行元素，J阵元素的计算过程分为以下四步。

1)在J阵第1～2m行中第1～2m列(J阵虚线左上角)部分的形成过程中，仅计算奇数行的非零元素H_ij、N_ij，利用H_ij＝L_ij，N_ij＝-M_ij的关系得到相应偶数行的非零元素，但对对角元由于H_ij≠L_ij，N_ij≠-M_ij，需对其做补充或单独计算；

2)在第1～2m行中第2m列以右部分(右上角)的形成过程中，两行/次计算相应行的非零元素H_ij，M_ij；

3)在第2m行以下的第1～2m列(左下角)部分的形成过程中，一行/次计算各行的非零元素H_ij，N_ij；

4)在第2m行以下的第2m列以右部分(右下角)的形成过程中，一行/次计算各行的非零元素H_ij。

由于Y(n,d)数组的结构完全反映网络结构的稀疏性，且Y阵与J阵元素结构几乎相同、且非零元素位置相对应，因此在J阵的形成过程中不用对Y阵元素进行非零判断，从而大大提高计算速度。这种方式还可避免仅存贮上三角非零元素的方式中利用对称性确定下三角非零元素的不便和过多的下标转换等。

(4)Y(n,d)数组中主节点的行号可便于数据的检索和检验，节点数S_i可保证对Y(n,d)数组数据的读取效率。

步骤4：利用对称稀疏矩阵技术对J阵进行消元和回代求取潮流；

(1)先定义对含规格化的第k列元素消元时计算中所用的各个元素的如下。

以对角元为界，对角元以右的元素为交叉元素，对角元以下的元素为消元元素，非零的消元元素所在行和非零的交叉元素所在列相交点上的元素为需计算的计算元素。

(2)对J阵第1列元素消元，其非零元素的判断、规格化和消元过程如下。

1)判断J阵中第1行2个对角元H₁₁、N₁₁以右第1～2m列中所有奇数列非零的交叉元素H_1j，如H_1j≠0，可得N_1j≠0及第2行的M_1j≠0、L_1j≠0；根据第1、2行中第1～2m列非零的交叉元素H_1j、N_1j、M_1j、L_1j，按对称性可得对角元以下第1、2列中第1～2m行非零的消元元素H_j1、N_j1、M_j1、L_j1。该部分仅针对J阵元素的左上角部分。

2)判断第1行第2m列以右各列非零的交叉元素H_1j，如H_1j≠0，可得第2行的M_1j≠0；根据第1、2行中第2m列以右非零的交叉元素H_1j、M_1j，按对称性可得对角元以下第1、2列中第2m行以下非零的消元元素H_j1、N_j1。该部分仅针对J阵元素的右上角和左下角部分。

3)规格化第1行不为零的H_1j、N_1j元素；

4)计算第1列非零的消元元素H_j1、M_j1所在行和第1行非零的交叉元素H_1j、N_1j所在列相交的所有元素，完成第1列消元。

(3)对第2列元素消元过程如下。

1)规格化第1行非零的交叉元素判断中得到的第2行非零的交叉元素M_1j、L_1j；

2)计算根据(1)中得到的第2列非零的消元元素N_j1、L_j1所在行(对第1列消元过程中N_j1、L_j1的值可能变化，但其非零性不变)和第2行非零的交叉元素M_1j、L_1j所在列相交的所有元素(对第1列消元过程中M_1j、L_1j的值可能变化，但其非零性不变)，完成第2列消元。

(4)同理，依次循环。

对J阵元素的判断可分以下三部分：

1)判断第1～2m行2个对角元以右的第1～2m列奇数行中奇数列非零的交叉元素，可同时得到第1～2m列奇数行和偶数行及奇数列和偶数列所有非零的交叉元素，如H_ij≠0，可得到N_ij≠0、M_ij≠0、L_ij≠0；按对称性可得对角元以下第1～2m行的第1～2m列中非零的消元元素H_ji、N_ji、M_ji、L_ji。由于只涉及到J阵左上角部分的元素相互对称，且判断1个非零元素可决定8个非零元素，因此可称其为“判1定8”。

2)判断第1～2m行第2m列以右各列奇数行非零的交叉元素H_ij，如H_ij≠0，可得偶数行的M_ij≠0；根据第1～2m行中第2m列以右非零的交叉元素H_ij、M_ij，按对称性可得对角元以下、及第2m行以下的第1～2m列中非零的消元元素H_ji、N_ji。由于只涉及到J阵右上角和左下角部分的元素相互对称，且判断1个非零元素可决定4个非零元素，因此可称其为“判1定4”。

3)判断第2m行以下的第2m列以右各列非零的交叉元素H_ij，如H_ij≠0，按对称性可得对角元以下相应行中非零的消元元素H_ji。由于只涉及到J阵右下角部分的元素相互对称，且判断1个非零元素可决定2个非零元素，因此可称其为“判1定2”。这部分仅针对PV节点部分，全部为H阵元素。实际计算中当系统PV节点数在10％左右时，对这部分可以不做非零判断，而直接进行消元计算。

(5)记住第1次上三角非零元素的坐标，就可利用J阵消元过程中非零元素坐标不变的特性，直接完成后续的多次前代和回代计算。

(6)第1次所记录的上三角非零元素的坐标还可用于回代过程求解Δδ_i、ΔV_i/V_i，进一步提高潮流计算的速度。

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

如果不满足收敛条件，则利用第一次迭代过程中记录的上三角非零元素的坐标继续进行后续的消元和回代计算；如果满足收敛条件，则执行步骤6。

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

步骤4中，对J阵用J_ij子阵进行分析更能简化计算过程。

1)形成J阵的过程中如考虑Y阵元素的稀疏性，则当B_ij＝0，一般可得Y_ij＝0和H_ij＝0，从而可得J_ij＝0和J_ji＝0；若B_ij≠0，则Y_ij≠0，因此可得J_ij≠0和J_ij≠0，即H_ij元素的非零性可由B_ij元素的非零性来决定。由于H_ij元素的特性，在任何情况下均可用H_ij元素的非零性来确定J_ij子阵的非零性。而由于J_ij子阵中H_ij、N_ij、L_ij、M_ij元素的非零性与Y阵元素G_ij、B_ij的非零性一一对应，因此每次消元迭代过程中所形成的上三角元素中非零元素的坐标始终保持不变。

2)由于J_ii子阵非零，所以只需对对角元J_ii子阵以右的J_ij子阵进行非零判断，即只需对J_ij子阵中的H_ij元素进行非零判断。如果H_ij≠0，则J_ij≠0和对角元J_ii子阵以下的J_ji≠0，记录非零的H_ij、N_ij元素的坐标或仅H_ij元素的坐标以便后续应用。由于判断H_ij是否为零，便可确定J_ij与J_ji子阵的非零性，这样在J阵左上三角的第1～2m行和第1～2m列的判断中便减少了7/8的判断语句；在右上三角的第1～2m行和第2m列以右以及左下三角的第2m行以下和第1～2m列的判断中减少了3/4的判断语句；在右下三角的第2m行以下和第2m列以右的判断中减少了1/2的判断语句，即分别为“判1定8”、“判1定4”和“判1定2”。

3)仅计算非零的消元元素所在行与非零的交叉元素所在列相交点上的计算元素，从而大大减少计算元素的计算量，极大地提高了计算效率。而在之后每次对J阵的消元和回代时可直接利用第1次消元过程所记录的上三角非零元素的坐标，又大大减少了后续迭代过程中大量不必要的判断语句。

本发明方法无论在读取数据文件、形成J阵、对J阵进行消元和回代计算等方面的计算速度均远快于传统方法中不考虑或考虑元素稀疏性的计算速度。用本发明方法对各IEEE系统进行验算完全可行。如对IEEE-118系统，读取数据文件的时间减少约82％，形成J阵和对J阵进行消元和回代的时间减少约66～90％，整个潮流计算的时间减少约80％～85％。且系统节点数越多，本发明方法的优势越大。

附图说明

图1为本发明方法进行潮流计算的流程图。

图2为本发明方法进行极坐标牛顿-拉夫逊法潮流计算的流程图。

具体实施方式

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

实施例1。分别比较传统方法中不考虑元素稀疏性和按列判断Y、J阵的非零元素以及本发明方法对IEEE-30、-57、-118节点系统读取Y阵数据文件的时间、形成J阵及消元的时间以及完成潮流计算的总时间，比较结果如表1所示。

表1对IEEE系统读取Y阵数据文件、形成J阵及消元和回代及潮流计算时间的比较

T₁、T₂：传统方法不考虑和考虑元素稀疏性时读取Y(n,2n)数据文件的平均时间。

T₃：本发明方法读取Y(n,d)数据文件的平均时间。

T′₁：传统方法中不考虑元素稀疏性直接形成J阵和对J阵消元和回代计算的累加平均时间。包括：不判断Y(n,2n)的非零元素直接形成J阵、不判断J阵的非零元素直接对J阵元素进行消元和回代。

T′₂：传统方法中部分考虑元素稀疏性时形成J阵和对J阵消元和回代计算的累加平均时间。包括：多次仅判断Y(n,2n)的虚部元素形成J阵、多次按列判断J阵中非零的消元元素进行消元和回代。

T′₃：本发明方法形成J阵和对J阵消元和回代计算的累加平均时间。包括：用Y(n,d)数组直接形成J阵、采用“判1定8”法、“判1定4”法、“判1定2”法、仅计算非零元素交叉点上的元素、仅计算对角元及其以右的非零元素等技巧对J阵元素进行消元和回代。

T”₁：传统方法不考虑元素稀疏性时计算潮流的总平均时间。

T”₂：传统方法部分考虑元素稀疏性时计算潮流的总平均时间。

T”₃：本发明方法计算潮流的总平均时间。

以IEEE-118节点系统为例，根据表1可以看出：

(1)本发明方法读取Y阵数据文件的时间为传统方法的18.85％。

(2)本发明方法形成J阵及对J阵消元和回代的累加时间为传统方法中不考虑元素稀疏性方法的9.26％，为传统方法中部分考虑元素稀疏性方法的33.85％。

(3)本发明方法潮流计算的时间为传统方法中不考虑元素稀疏性方法的15.56％，为传统方法中部分考虑元素稀疏性方法的20.75％。

(4)读取Y阵数据文件的时间基本约占潮流计算的80％，而形成J阵和对J阵消元和回代计算的时间约占潮流计算的20％，说明快速形成J阵和潮流计算的关键在于Y阵数据文件的读取。

因此，可得出以下结论：

(1)在读取Y阵数据文件、形成J阵及对J阵消元和回代、求取潮流等计算过程中，本发明方法的计算速度大大优于传统方法中不考虑或部分考虑稀疏性的情况。

(2)电力系统节点数越多，本发明方法的优势越大。

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

Claims (2)

1.一种基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法潮流的方法，其特征包括以下步骤：

步骤1：建立仅含非零元素的Y阵数据文件；

步骤2：打开数据文件，读取仅含非零元素的Y阵数据文件到Y(n,d)数组；其中，n为系统的节点数，d为Y(n,d)数组中的列数，d＝3p+2，p＝l_max+1，l_max为系统中各节点的最大连接支路数；

步骤3：根据Y(n,d)数组中的每行元素，分四步计算J阵中的非零元素；

1)对J阵第1～2m行中第1～2m列，仅计算奇数行的非零元素H_ij、N_ij，利用H_ij＝L_ij，N_ij＝-M_ij的关系得到相应偶数行的非零元素；

2)对J阵第1～2m行中第2m列以右部分，两行/次计算相应行的非零元素H_ij，M_ij；

3)对J阵第2m行以下的第1～2m列，一行/次计算各行的非零元素H_ij，N_ij；

4)对J阵第2m行以下的第2m列以右部分，一行/次计算各行的非零元素H_ij；

步骤4：利用对称稀疏矩阵技术对J阵进行消元和回代求取潮流：

(1)设系统的PQ节点数为m，对J阵非零元素的判断可分为以下三步：

1)判断第1～2m行中第1～2m列2个对角元以右的奇数行和奇数列1个非零的交叉元素，可同时得到相应的奇数行和偶数行及奇数列和偶数列4个非零的交叉元素，按对称性可得2个对角元以下4个非零的消元元素；

2)判断第1～2m行中第2m列以右各列奇数行1个非零的交叉元素，可得相应的奇数行和偶数行2个非零的交叉元素，按对称性可得相应对角元以下2个非零的消元元素；

3)判断第2m行以下的对角元以右1个非零的交叉元素，按对称性可得相应对角元以下1个非零的消元元素；

(2)对J阵元素的消元可分为以下二步：

1)对第1～2m列元素消元时，判断对角元以右第1～2m列中奇数行的奇数列或第2m列以右奇数行中各列的1个非零元素，可同时确定8个或4个相应的非零元素，分步完成对奇数行的规格化和对奇数列的消元、对偶数行的规格化和对偶数列的消元；

2)对第2m列以右各列消元时，判断每行对角元以右各列的1个非零元素，分别完成对每行对角元以右非零元素的规格化以及对对角元以下相应元素的消元；

在上述消元过程中，记住第1次判断的上三角奇数行非零元素的坐标，利用J阵消元过程中非零元素坐标不变的特性，直接完成后续的多次前代和回代计算；

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

如果不满足收敛条件，则利用第一次迭代过程中记录的上三角非零元素的坐标继续进行后续的消元和回代计算；如果满足收敛条件，则执行步骤6；

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

2.根据权利要求1所述的基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法潮流的方法，其特征是所述的步骤4中，对J阵用J_ij子阵进行分析计算：

1)形成J阵的过程中考虑Y阵元素的稀疏性，则当B_ij＝0，可得Y_ij＝0和H_ij＝0，从而可得J_ij＝0和J_ji＝0；若B_ij≠0，则Y_ij≠0，可得J_ij≠0和J_ji≠0；由于J_ij子阵中H_ij、N_ij、L_ij、M_ij元素的非零性与Y阵元素G_ij、B_ij的非零性一一对应，每次消元迭代过程中所形成的上三角元素中非零元素的坐标始终保持不变；

2)由于J_ii子阵非零，只对对角元J_ii子阵以右的J_ij子阵进行非零判断；如果H_ij≠0，则J_ij≠0和对角元J_ii子阵以下的J_ji≠0，记录非零的H_ij、N_ij元素的坐标或仅H_ij元素的坐标以便后续应用；

3)在之后每次对J阵的消元和回代时可直接利用第1次消元过程所记录的上三角非零元素的坐标。