CN105703359B - 一种对称稀疏因子表法在直角坐标pq分解法潮流计算中的应用 - Google Patents

Abstract

一种对称稀疏因子表法在直角坐标PQ分解法潮流计算中的应用，属于电力系统分析计算领域。本发明从给定结构的数据文件中读取数据，大大提高了数据读取速度和I_pi、I_qi或ΔP_i、ΔQ_i的计算速度；利用对称稀疏技术快速形成因子表；记录其上三角非零元素的列号及其个数；将对称稀疏技术用于后续的快速前代及回代计算；在有功、无功迭代中引入塞德尔方式加快潮流计算速度。如对IEEE‑118系统，本发明与传统方法相比，读取数据文件、形成因子表、有功及无功迭代的时间及潮流计算总时间分别为后者的6.12％、6.15％、11.86％和6.16％。且系统节点数越多，本发明的优势越大。

Description

一种对称稀疏因子表法在直角坐标PQ分解法潮流计算中的 应用

技术领域

本发明属于电力系统分析计算领域，涉及直角坐标PQ分解法潮流计算。

背景技术

牛顿-拉夫逊法(牛顿法)是电力系统潮流计算中最常用的方法，按计算过程可分直角坐标和极坐标牛顿法，简化后可分别派生出直角坐标和极坐标PQ分解法。从计算过程上来看，极坐标牛顿法比直角坐标牛顿法修正方程式的方程数要少，雅可比矩阵元素的种类和数量也少，有时潮流迭代次数也少一次，因此极坐标比直角坐标牛顿法应用似乎更为广泛。同理，极坐标PQ分解法得到了广泛应用，而几乎无文献介绍直角坐标PQ分解法。实际上，由于直角坐标潮流计算中不涉及三角函数，因此在考虑稀疏性和直角坐标潮流计算系数矩阵结构特点的情况下其存贮单元数、潮流迭代速度等指标并不亚于极坐标法潮流计算。

电力系统中PQ分解法潮流计算中涉及到多个对称稀疏矩阵应用。若不考虑稀疏性，大量零元素的存贮、读取及计算过程中会导致存贮空间的极大浪费、读写数据文件及计算过程时间过长。因此，结合直角坐标PQ分解法系数矩阵的特点、考虑稀疏技术对元素的存贮方式及形成因子表的过程不但可大量节省存贮单元，还可大大减少数据文件的存贮、读取及计算过程的时间。

传统直角坐标PQ分解法潮流计算需3个数组，导纳矩阵Y用于计算节点电流(I_pi、I_qi)和节点功率(ΔP_i、ΔQ_i)，系数矩阵B′、B″用于求取节点电压的虚部增量Δf_i和实部增量Δe_i。传统方法中对3个数组数据的存贮、读取、应用等存在以下不足：

(1)B′、B″阵元素的获取方式不合理。

B′、B″阵元素的构成对潮流收敛性或收敛速度影响较大，因此其构成应各不相同。如果B′、B″阵元素直接取自Y阵元素的虚部，或如果仅B′、B″阵元素相同、而不同于Y阵元素的虚部，都可能影响PQ分解法潮流的迭代次数或收敛性。

(2)Y、B′、B″阵元素的存贮和读取方式不合理。

设系统节点数为n，PQ节点数为m，传统直角坐标PQ分解法中对应的数组分别为Y(n,2n)、B′(n-1,n-1)、B″(n-1,n-1)，这种存贮方式会造成存贮空间的极大浪费和读写数据文件及计算时间过长，且用Y(n,2n)数组计算I_pi、I_qi或ΔP_i、ΔQ_i效率极为低下。如果用坐标存贮、顺序存贮、链表存贮等方法，尽管可节省存贮单元，但由于对角元与非对角元分开存贮使得数据的读取过程繁琐，不利于其计算及处理。此外，坐标存贮、顺序存贮、链表存贮等方法比较适用于实数元素的存贮，对复数元素的存贮则略为复杂。此外，在对Y、B′、B″阵元素的存贮和读取时未利用3个矩阵中绝大部分非零元素的位置相同的特点。

(3)Y、B′、B″阵数据文件的存贮个数较多、读写时间较长。

如不考虑稀疏性按传统方法分别存贮Y、B′、B″阵的3个数据文件，则在潮流计算程序中要分别打开3个数据文件并读取数据；如考虑稀疏性按链表存贮方法分别对Y、B′、B″阵进行存贮，则存贮文件的个数可达到9个，不利于PQ分解法的实时计算。

(4)Y、B′、B″阵非零元素位置之间的关系未充分利用。

形成Y、B′阵的过程中，在1～(n-1)阶范围内，其静态的非零元素的位置和个数完全一样；而在对Y、B′阵的消元过程中，在1～(n-1)阶范围内，其动态的非零元素的位置和个数也完全一样。在第1～m行，B′、B″阵各行非零元素的位置和个数在静态和动态情况下完全相同；而在第(m+1)～(n-1)行，B″阵的中除对角元为-2外，其余元素均为零，与B′阵第(m+1)～(n-1)行中在静态和动态情况下非零元素的位置和个数不同，且在消元过程中B″阵除对角元由于规格化发生变化外，其余元素均不变化。传统方法中由于未利用该特性而导致存贮单元的浪费和因子表形成过程中一些不必要的计算。

(5)对B′、B″阵形成因子阵的方式不合适。

1)形成因子阵的方法不合适。因P_i-Δf_i和Q_i-Δe_i迭代均为常系数线性方程，可用因子表法或LR、LDU、CU三角分解法等求解，一般多用LDU三角分解法。但LDU三角分解法需求解L、D、U三个因子阵，而因子表法计算过程简单，且比LDU三角分解法的计算速度更快。

2)形成因子表时稀疏技术应用不足。B′、B″阵和Y阵元素结构几乎相同，且极度稀疏，若不使用稀疏技术，则计算效率极为低下。三角分解法中稀疏技术应用并不多见，而因子表法中稀疏技术的应用极为广泛，但并不完全到位。

3)对B′、B″阵因子表分开形成。在PQ节点的范围内，B′、B″阵的元素结构完全相同，如果不能根据对B′阵形成因子表的过程快速形成B″阵的因子表，则形成B″阵因子表时要浪费时间。

(6)用因子表进行前代和回代过程中稀疏技术应用不到位。

形成因子表后，有时传统方法会利用因子表上三角的非零元素快速回代求取Δf_i、Δe_i，但各行非零元素的个数和列号分开存贮未能使其计算方式达到最优，而对ΔP_i/e_i、ΔQ_i/e_i的前代计算则未利用稀疏技术。

(7)有功、无功迭代过程中未引入塞德尔迭代方式。

引入塞德尔迭代方式不但对PQ分解法的计算速度有较大影响，有时还会影响其收敛性，而传统PQ分解法中较少引入塞德尔迭代方式。

由于综上所述原因，传统直角坐标PQ分解法的计算速度远未达到最优，从而几乎不被应用。

发明内容

为了克服上述现有技术的不足，本发明提供了一种对称稀疏因子表法在直角坐标PQ分解法潮流计算中的应用。

本发明读取的数据文件为仅由非零元素构成的三个虚拟数组Y(n,3d₁)、B′(n-1,d₂)、B″(m,d₃)所组成的A(n,d)数据文件，加速数据文件的读取和I_pi、I_qi或ΔP_i、ΔQ_i的计算；利用对称稀疏技术加速因子表的形成；记录B′(n-1,n-1)数组因子表B′^(n-2)的形成过程中上三角非零元素的个数和列号，并将其用于加速B″(m,n-1)数组因子表B″^(m-1)的形成；将对称稀疏技术继续用于因子表B′^(n-2)、B″^(m-1)中加速后续的前代及回代计算；在迭代过程中引入塞德尔迭代方式，加速有功、无功迭代。

本发明是通过以下技术方案实现的。

本发明所述的一种对称稀疏因子表法在直角坐标PQ分解法潮流计算中的应用，包括以下步骤：

步骤1：定义数组Y(n,3d₁+1)、B′(n-1,n-1)、B″(m,n-1)；

步骤2：从A(n,d)数据文件中将虚拟数组Y(n,3d₁)、B′(n-1,d₂)及B″(m,d₂)的数据分别读入至Y(n,3d₁+1)、B′(n-1,n-1)、B″(m,n-1)数组；

(1)在应用PQ分解法潮流程序之前先建立无非零元素的A(n,d)数据文件，将虚拟数组Y(n,3d₁)、B′(n-1,d₂)、B″(m,d₂)对应的数据全部存于A(n,d)数组中，其中d＝3d₁+2d₂+3，Y(n,3d₁)对应的列数为3d₁，B′(n-1,d₂)、B″(m,d₂)对应的列数均为d₂，“+3”分别为行号列和S_1max、S_2max所在的计数列，即表1中的第1～3列；且根据数组结构有d₁＝S_1max，d₂＝S_2max，S_1max、S_2max为数组各行中最大的非零元素节点之和。

(2)传统方式中，Y阵的阶数为n，而B′、B″阵阶数均为(n-1)。在1～(n-1)阶范围内，Y、B′阵的阶数相同时其非零元素的个数和列号完全一样；在第1～m行范围内，B′、B″阵各行非零元素的个数和列号完全相同，而在第(m+1)～(n-1)行，B″阵除对角元为-2外，其余元素均为零，与B′阵不同。因此，只要给定Y(n,3d₁)各行中对应的父节点和非零元素子节点之和为S_1max，B′(n-1,d₂)和B″(m,d₂)对应的为S_2max。所以Y(n,3d₁)中只需存贮主节点及非零元素子节点的列号j和参数，而B′(n-1,d₂)、B″(m,d₂)中只需存贮所有对应的非零元素参数，而非零元素的列号j可根据S_2max由Y(n,3d₁)中得到。

(3)利用S_i1快速将A(n,d)中Y(n,3d₁)的数据读入Y(n,3d₁+1)数组，“+1”表示将S_i1读入到Y(n,3d₁+1)数组的第1列；利用S_i2快速将B′(n-1,d₂)、B″(m,d₂)的数据分别读入B′(n-1,n-1)、B″(m,n-1)数组，但不用将S_i2读入到B′(n-1,n-1)、B″(m,n-1)数组，且B″阵第(m+1)～(n-1)行的元素除对角元为-2外，其余元素均为零，因此也无需读入。

(4)A(n,d)数组分为“行号组”、“节点数组”、“Y阵组”、“B′阵组”、“B″阵组”，其元素存贮方式如表1所示。

表1A(n,d)数组元素的存贮方式

注：表1中除含最大非零元素节点对应的行外，并非所有的存贮单元都有数据，而S_i1、S_i2可保证实际存贮单元数远远小于最大存贮单元数，进一步提高数据的存贮和读写效率。

(5)从A(n,d)数据文件中利用S_i1值和S_i2值快速将虚拟数组Y(n,3d₁)、B′(n-1,d₂)、B″(m,d₃)对应的数据分别读入Y(n,3d₁+1)、B′(n-1,n-1)及B″(m,n-1)数组，利用S_i1值和Y(n,3d₁+1)可快速计算I_pi、I_qi或ΔP_i、ΔQ_i，用B′(n-1,n-1)、B″(m,n-1)求解Δf_i、Δe_i。

步骤3：根据对称稀疏技术快速形成B′(n-1,n-1)、B″(m,n-1)的因子表B′^(n-2)、B″^(m-1)，并记录形成B′^(n-2)因子表中上三角每行非零元素的个数和列号，而B′^(n-2)因子表中第1～m行非零元素的个数及列号可直接用于B″^(m-1)的形成过程中；

(1)用对称稀疏技术快速形成B′(n-1,n-1)的因子表B′^(n-2)，自动统计记录因子表中每行上三角非零元素u′_ij的个数S′_i和列号，以便在后续计算中根据S′_i和相应的列号直接取用B′^(n-2)的非零元素反复对P_i-Δf_i方程进行快速前代和回代，计算要点如下。

1)快速确定非零元素和计算元素。根据对角元以右非零元素的值和位置按对称性获得对角元以下非零元素的值和位置，仅计算各非零元素行与列交互点上的元素。

在“逐行规格化、按列消元”过程中，B′^(n-2)的每行元素在规格化前，其对角元以右的非零元素与对角元以下的非零元素数值相等、位置对称；在规格化后，对角元以右的非零元素与对角元以下的非零元素数值不等、但位置仍然对称。因此可根据对角元以右的非零元素的值和位置按对称性快速获得对角元以下的非零元素的值和位置，并利用相应的非零元素所在行与所在列的交互点来确定需计算的元素，大大提高对非零元素的判断和计算效率。

2)省去对下三角非零元素的计算。根据1)，对对角元以右的元素进行规格化前，先将其赋给相应位置的对角元以下的素。因此在形成因子表的过程中只需计算对角元及其以右的、非零元素行列交互点上的元素，可减少约50％非零的非对角元素的计算。

3)自动统计并记录上三角非零元素个数和列号。自动统计记录B′^(n-2)中每行上三角非零元素u′_ij的个数S′_i和列角标，并存入相应数组，以便在后续计算中根据S′_i和相应的列号直接取用B′^(n-2)的非零元素反复快速进行前代和回代，免去该过程对零元素的无效计算。

(2)根据对B′^(n-2)中记录的S′_i和列号快速形成B″(m,n-1)的因子表B″^(m-1)，以便在后续计算中根据S′_i和相应的列号直接取用B″^(m-1)的非零元素反复对Q_i-Δe_i方程进行快速前代和回代，计算要点如下；

1)B′、B″阵的实际阶数均为(n-1)。但在形成B′、B″阵的过程中，第1～m行B′、B″阵各行静态非零元素b′_ij、b″_ij的个数S_i2和列号完全相同，而在第(m+1)～(n-1)行，B″阵除对角元为-2外，其余元素均为零，与B′阵各行静态非零元素的个数S_i2和列号完全不同。

2)在形成B′、B″阵因子表的过程中，第1～m行B′、B″阵因子表各行动态非零元素u′_ij、u″_ij的个数S′_i和列号完全相同；而在第(m+1)～(n-1)行，B″阵因子表除对角元外，其余元素仍全部为零，与B′阵因子表各行动态非零元素的个数S′_i和列号完全不同。

3)在形成因子表的消元过程中，对B″阵的前m行的m列，在对对角元以右的元素进行规格化前，可根据对角元以右非零元素的值和位置获得对角元以下非零元素的值和位置，省去对下三角非零元素的计算；对B″阵的前m行的(m+1)～(n-1)列，仅需计算其非零元素交互点上的元素；对B″阵(m+1)～(n-1)行的所有元素，由于只有对角元-2，其它元素均为零，且在消元过程中，除对角元-2在规格化时变化外，所有零元素均不变化，因此无需在B″阵的因子表中考虑(m+1)～(n-1)行非对角元的动态元素，即对B″阵的B″(m,n-1)数组只要形成含PQ节点部分的因子表B″^(m-1)。

4)由于在形成B″阵因子表B″^(m-1)的过程中，其前m行非零元素u″_ij的个数和列号与B′阵因子表B′^(n-2)中非零元素u′_ij的个数S′_i和列号完全一样，因此无需对其再进行非零判断和自动统计记录，而在后续计算中可根据对应B″^(m-1)的S′_i和相应的列号直接取用B″^(m-1)的非零元素反复对Q_i-Δe_i方程进行快速前代和回代。

步骤4：根据对称稀疏技术利用因子表B′^(n-2)中下三角非零元素l′_ji和S′_i快速前代计算ΔP_i/e_i，利用其上三角非零元素u′_ij和S′_i快速回代求取Δf_i。利用因子表B″^(m-1)中下三角非零元素l″_ji和S′_i快速前代计算ΔQ_i/e_i，利用其上三角非零元素u″_ij和S′_i快速回代求取Δe_i，并在有功、无功迭代过程中引入塞德尔方式；

(1)通过因子表形成过程中记录的B′^(n-2)中上三角非零元素u′_ij的位置，按对称性得到下三角非零元素l′_ji的位置，用l′_ji元素和S′_i快速前代计算ΔP_i/e_i；用u′_ij元素和S′_i快速回代求取Δf_i；通过B″^(m-1)中上三角非零元素u″_ij的位置，按对称性得到下三角非零元素l″_ji的位置，用l″_ji元素和S′_i快速前代计算ΔQ_i/e_i；用u″_ij元素和S′_i快速回代求取Δe_i，大大加快前代、回代过程。

(2)引入塞德尔迭代方式后，PQ分解法的迭代过程为：

即先计算有功，然后迭代有功求出f_i的修正量和新值，再用f_i的新值计算无功，继而迭代无功求出e_i的修正量和新值，再用e_i的新值计算有功，依次循环加快收敛。

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

如果不满足收敛条件，则利用本次迭代得到的Δf_i、Δe_i继续进行下一次的前代和回代计算；如果满足收敛条件，则执行步骤6。

步骤6：计算平衡节点的功率及支路功率并输出计算结果。

本发明方法与传统的PQ分解法相比，具有以下优点。

(1)仅读取A(n,d)数组1个数据文件，而A(n,d)数组中仅存贮相应非零元素的数据，存贮单元数及数据文件的读取时间大大减少。

(2)可直接根据Y(n,3d₁+1)数组高效计算I_pi、I_qi或ΔP_i、ΔQ_i。

(3)形成B″阵和B″阵的因子表时考虑其PV节点在形成和消元过程的特点，对(n-1)阶的B″阵仅需形成B″(m,n-1)数组，而对B″(m,n-1)数组只要形成含PQ节点部分的B″^(m-1)因子表。

(4)根据对称稀疏性快速获得B′、B″阵的因子表，并仅记录B′阵因子表中上三角非零元素的个数S′_i及列号，以便后续计算中直接取用B′^(n-2)、B″^(m-1)阵的非零元素反复快速前代和回代。

(5)利用下三角非零元素l′_ji、l″_ji和S′_i可快速前代计算ΔP_i/e_i、ΔQ_i/e_i。

(6)利用上三角非零元素u′_ij、u″_ij和S′_i可快速回代求取Δf_i、Δe_i。

(7)有功、无功迭代过程中引入塞德尔迭代方式，加快潮流收敛。

附图说明

图1为传统直角坐标PQ分解法潮流的计算流程图。

图2为本发明方法用于直角坐标PQ分解法潮流的计算流程图。

具体实施方式

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

实施例1。分别用传统的直角坐标PQ分解法和本发明法用于直角坐标PQ分解法对IEEE-30、-57、-118节点系统进行潮流计算，比较其读取数据文件、形成因子表、有功无功迭代和潮流计算(总)的平均计算时间。计算结果如表2所示。

表2传统方法和本发明对IEEE系统进行PQ分解潮流计算过程时间的比较

t_r.c、t_f.c、t_i.c、t_p.c：传统PQ分解法中不考虑稀疏性时，读取数据文件、形成因子表、有功无功迭代的平均计算时间及潮流计算(总)时间；

t_f.s.c、t_p.s.c：传统方法PQ分解法中仅判断对角元以下的非零元素形成因子表的平均计算时间及潮流计算(总)时间(除形成因子表外，其他同传统方法中不考虑元素稀疏性)；

t_r.new、t_f.new、t_i.new、t_p.new：本发明中读取数据文件时间、形成因子表、有功无功迭代的平均计算时间及PQ分解法潮流计算(总)时间；

t_r.new/t_r.c、t_f.new/t_f.c、t_i.new/t_i.c、t_p.new/t_p.c：本发明与传统方法不考虑稀疏性时读取数据文件时间、形成因子表时间、有功无功迭代时间、潮流计算(总)时间的百分比。

t_tfs.c/t_x.c、t_p.s.c/t_p.c：传统方法中仅判断对角元以右的非零元素形成因子表的时间和PQ分解法潮流计算(总)时间与传统方法中不考虑任何稀疏性时相应计算结果的百分比。

以IEEE-118节点系统为例，计算结果比较分析如下：

(1)读取数据文件数及时间的比较

传统因子表法分别将Y、B′、B″阵元素分别存放在Y(n,2n)、B′(n-1,n-1)和B″(n-1,n-1)数组中，需读取3个数据文件，而本发明方法存放在A(n,d)数组中，仅需读取1个数据文件。本发明读取数据文件的时间仅为传统方法的6.12％。

(2)形成因子表时间的比较

本发明形成因子表的时间仅为传统方法的6.15％；如果在传统方法中判断对角元以下的非零元素形成因子表，则其为传统方法中不考虑稀疏性时的7.84％。

(3)本发明方法有功、无功迭代时间仅为传统方法的11.86％。

(4)潮流计算(总)时间的比较

本发明潮流计算(总)的时间仅为传统方法的6.16％。传统方法中判断对角元以下的非零元素形成因子表与传统方法中不考虑稀疏性时形成因子表的速度相比虽然大大提高，但前者潮流计算的时间与后者相比几乎没有变化。

上述计算结果表明，本发明提出的方法与传统的直角坐标PQ分解法相比，无论在读取数据文件、形成因子表、有功无功迭代还是潮流计算方面，计算速度均大大加快，说明本发明方法的计算效率极高。且电力系统节点数越多，本发明方法优势越大。

本发明可以采用任何一种编程语言和编程环境实现，这里采用Visual C++语言，开发环境是Visual C++，运行平台是Intel(R)Core i7-4790CPU@3.60GHZ，内存8.00GB。

Claims (1)

1.一种对称稀疏因子表法在直角坐标PQ分解法电力系统潮流计算中的应用，其特征包括以下步骤：

步骤1：定义数组Y(n,3d₁+1)、B′(n-1,n-1)、B″(m,n-1)；其中，n为系统节点数，m为PQ节点数；

步骤2：从A(n,d)数据文件中将虚拟数组Y(n,3d₁)、B′(n-1,d₂)及B″(m,d₂)的数据分别读入至Y(n,3d₁+1)、B′(n-1,n-1)、B″(m,n-1)数组，其中Y(n,3d₁+1)用于快速计算I_pi、I_qi或ΔP_i、ΔQ_i，B′(n-1,n-1)、B″(m,n-1)数组用于求解Δf_i、Δe_i；其中，d₁＝S_1max，d₂＝S_2max，S_1max、S_2max为数组各行中最大的非零元素节点之和，d＝3d₁+2d₂+3；

步骤3：根据对称稀疏技术快速形成B′(n-1,n-1)、B″(m,n-1)的因子表B′^(n-2)、B″^(m-1)，并记录形成B′^(n-2)因子表中每行上三角非零元素的个数和列号，而B′^(n-2)因子表中第1～m行非零元素的个数及列号可直接用于B″^(m-1)的形成过程中；

(1)用对称稀疏技术快速形成B′(n-1,n-1)的因子表B′^(n-2)，自动统计记录因子表中上三角每行非零元素u′_ij的个数S′_i和列号，在后续计算中根据S′_i和相应的列号直接取用B′^(n-2)的非零元素反复对P_i-Δf_i方程进行快速前代和回代；

(2)根据对B′^(n-2)中记录的S′_i和列号用对称稀疏技术快速形成B″(m,n-1)的因子表B″^(m-1)，在后续计算中根据S′_i和相应的列号直接取用B″^(m-1)的非零元素反复对Q_i-Δe_i方程进行快速前代和回代；其中，在形成B′、B″阵和形成B′、B″阵因子表的过程中，第1～m行B′、B″阵静态和动态非零元素的个数S_i2和S′_i以及列号完全相同，而在第m+1～n-1行，B″阵除对角元外，其余元素均为零，与B′阵静态和动态非零元素的个数S_i2和S′_i以及列号完全不同；

步骤4：根据对称稀疏技术利用因子表B′^(n-2)中下三角非零元素l′_ji和S′_i快速前代计算ΔP_i/e_i，利用其上三角非零元素u′_ij和S′_i快速回代求取Δf_i；利用因子表B″^(m-1)中下三角非零元素l″_ji和S′_i快速前代计算ΔQ_i/e_i，利用其上三角非零元素u″_ij和S′_i快速回代求取Δe_i，并在有功、无功迭代过程中引入塞德尔方式；

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

如果不满足收敛条件，则利用本次迭代得到的Δf_i、Δe_i继续进行下一次的前代和回代计算；如果满足收敛条件，则执行步骤6；

步骤6：计算平衡节点的功率及支路功率并输出计算结果。

Images