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

摘要

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

说明书

技术领域

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

背景技术

在大型电力系统的潮流计算中，经常需根据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++编程语言，开 发环境是VisualC++。