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

Abstract

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

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的关系；未利用对角元以右和以下非零元素位置对称的关系；未利用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阵中奇数行的非零元素，而利用H_ij＝L_ij，N_ij＝-M_ij和R_ij＝S_ij＝0的关系得到相应偶数行的非零元素。由于Y(n,d)数组的结构完全反映网络结构的稀疏性，且Y阵与J阵元素结构几乎相同、且非零元素位置相对应，因此在J阵的形成过程中不用对Y阵元素进行非零判断，从而大大提高计算速度。还可避免传统方法中仅存贮上三角非零元素的方式中利用对称性确定下三角非零元素的不便和过多的下标转换等。

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

(5)初步计算J阵中所有的对角元素和非零的非对角元素，再根据计算出的I_pi、I_qi修正所有的对角元素，形成完整的J阵。

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

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

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

(2)对第1列进行消元：先判断J阵中第1行2个对角元H₁₁、N₁₁以右所有奇数列非零的交叉元素H_1j，如H_1j≠0，可得N_1j、M_1j、L_1j均不为零；仅规格化第1行不为零的H_1j、N_1j元素；根据H_1j、N_1j、M_1j、L_1j，按对称性可得对角元以下非零的消元元素H_j1、N_j1、M_j1、L_j1；计算第1列非零的消元元素H_j1、N_j1和M_j1、L_j1所在行和第1行非零的交叉元素H_1j、N_1j所在列相交的所有元素，完成第1列消元。

(3)对第2列进行消元：根据(1)中对第1行对角元以右非零的交叉元素的判断，规格化第2行对应的不为零的M_1j、L_1j元素；而在(1)中也得到了第2列对角元以下非零的消元元素N_j1、L_j1(对第1列消元过程中N_j1、L_j1的值可能变化，但其非零性不会变化)；计算第2列非零的消元元素N_j1、L_j1所在行和第2行非零的交叉元素M_1j、L_1j所在列相交的所有元素(同样对第1列消元过程中M_1j、L_1j的值可能变化，但其非零性不会变化)，完成第2列消元。

(4)同理，依次循环。

只要判断奇数行2个对角元以右奇数列非零的交叉元素，可同时得到奇数行和偶数行及奇数列和偶数列所有非零的交叉元素，即判断1个非零元素可确定4个非零元素；利用对称性，可得到其相应对角元以下奇数列和偶数列及奇数行和偶数行4个非零的消元元素，因此可称此方法为“判1定8”法；对奇数列消元时用奇数列和奇数行的非零元素交互点上的元素进行计算，对偶数列消元时用偶数列和偶数行的非零元素交互点上的元素进行计算。这样就可用对1行对角元以右奇数列元素的判断，分步完成对2行对角元以右非零元素的规格化以及对2列对角元以下相应元素的消元。

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

(6)第1次所记录的上三角奇数行非零元素的坐标还可用于回代过程求解Δe_i、Δf_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是否为零，便可确定J_ij与J_ji子阵的非零性，这样便减少了7/8的判断语句，即“判1定8”。并且只需计算非零的消元元素所在行与非零的交叉元素所在列相交点上的计算元素，从而大大减少计算元素的计算量，极大地提高了计算效率。而在之后每次对J阵的消元和回代时可直接利用第一次消元过程所记录的上三角奇数行非零元素的坐标，又大大减少了后续迭代过程中大量不必要的判断语句。如潮流迭代6次收敛，则本发明方法中对非零元素的判断只需进行1次。

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

附图说明

图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”法、仅计算非零元素交叉点上的元素、仅计算对角元及其以右的非零元素等技巧对J阵元素进行消元和回代。

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

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

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

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

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

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

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

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

因此，可得出以下结论：

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

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

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

Claims (2)

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

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

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

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

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

(1)判断奇数行2个对角元以右奇数列1个非零的交叉元素H_ij，同时得到奇数行和偶数行与奇数列和偶数列4个非零的交叉元素H_ij、N_ij、M_ij、L_ij；

(2)利用子阵J_ij与子阵J_ji的对称性，得到其相应的2个对角元以下奇数列和偶数列与奇数行和偶数行4个非零的消元元素H_ji、N_ji、M_ji、L_ji；

(3)对记录的奇数行的非零元素规格化，根据对称性得到奇数列的非零元素，再对其消元；

(4)根据记录的奇数行的非零元素对偶数行的非零元素规格化，根据对称性得到偶数数列的非零元素，再对其消元；

(5)对奇数列消元时用对奇数列和奇数行的非零元素交互点上的元素进行计算，对偶数列消元时用偶数列和偶数行的非零元素交互点上的元素进行计算；

(6)在上述消元过程中，记录第1次判断的上三角奇数行非零元素H_ij的坐标，利用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_ij≠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元素的坐标以便后续应用；

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