CN110649624A - 一种电力系统潮流并行计算方法 - Google Patents

Abstract

本发明公开了一种电力系统潮流并行计算方法，该方法包括以下步骤：1)基于牛拉法对电力系统进行潮流计算，获取非线性潮流方程，并依据逐次线性原理得到潮流修正方程组；2)基于BICGSTAB法求解潮流修正方程组，获取由改进PPAT预处理、改进Jacobi预处理组成的两阶段预处理算法，并对潮流修正方程组的系数矩阵即雅可比矩阵进行预处理，以改善BICGSTAB法的收敛性；3)基于CPU‑GPU异构平台运行步骤2)的并行计算过程。本发明不仅计算速度较快，同时鲁棒性较强，可适用于一些传统方法无法收敛的系统；本发明能够实现电力系统潮流的快速求解。

Description

一种电力系统潮流并行计算方法

技术领域

本发明涉及电力系统领域，尤其涉及一种电力系统潮流并行计算方法。

背景技术

潮流计算是电力系统安全稳定分析、故障计算、经济调度等问题的基础。准确、快速的潮流计算结果是电网安全可靠运行的重要保证。随着区域电网互联规模持续扩大、可再生能源大规模并网、电力电子装备大规模应用，电力系统潮流计算的规模和复杂度急剧增加。同时，电网智能化进程不断加深，电网精细化运行要求不断提高，运行人员对电力系统潮流计算精度和效率提出了更高的要求。

牛拉法是目前电力系统潮流计算最为广泛采用的方法，该方法基于逐次线性逼近原理将非线性潮流方程的求解转化为对线性修正方程组的求解^[1]。而修正方程组的求解是牛拉法中最为耗时的部分，约占整个潮流计算耗时的80％^[2]。因此，减少该部分求解耗时是提升电力系统潮流计算效率的关键。

当前，修正方程组的求解主要采用直接法和迭代法。直接法^[3]基于稀疏矩阵分解技术，逐一求解各修正量的解，在此过程中各修正量间具有极强的相关性，不易实现修正方程组求解的并行化^[4]。同时，由于直接法在对雅可比矩阵进行分解时，存在非零元素的填充问题，当系统规模增大时，会导致数据存储量大幅增加，不利于大规模互联电力系统潮流的快速求解。相对直接法，迭代法在求解修正方程组时，通过不断修正修正量的解向量，使修正量的解向量逐步逼近潮流方程实际解向量。该过程主要涉及稀疏矩阵和向量间运算，数据相关性弱，具有良好并行特性，适合于处理高维修正方程组。

迭代法主要分为古典迭代法与Krylov子空间迭代法两类。古典迭代法收敛性较差，目前已较少应用；Krylov子空间迭代法因其优良的收敛性与数值稳定性，已在潮流计算、小干扰分析、暂态稳定性分析等领域得到应用^[5]-[8]。但每种Krylov子空间法均有一定适用范围，采用Krylov子空间法求解潮流修正方程组时，针对雅可比矩阵的不对称不定特征，常用的Krylov子空间法为广义最小残差(generalized minimal residual，GMRES)法和BICGSTAB法。GMRES法可较精确实现修正方程组的求解，且具有良好的收敛性，但在构造Krylov子空间基底时需要长递推式，系统规模增大时将导致迭代次数、数据存储量和计算量急剧增加，不利于高维修正方程组的快速求解。相对GMRSE法，BICGSTAB法以其存储空间稳定性和短递推特性等优势，已在高维修正方程组的求解中得到了应用。但采用BICGSTAB法求解修正方程组时，需对潮流雅可比矩阵进行预处理^[9]，压缩其特征值分布区间，从而提高BICGSTAB法收敛性。而预处理器设计目前仍是制约BICGSTAB法应用于电力系统潮流并行计算的一个瓶颈。

如何快速、准确地实现电力系统潮流计算是当前研究重点。

发明内容

本发明提供了一种电力系统潮流并行计算方法，保证了潮流计算结果的正确性，对于大规模的系统，本发明不仅计算速度较快，同时鲁棒性较强，可适用于一些传统方法无法收敛的系统；与同样基于CPU-GPU异构平台的采用经典ILU(0)预处理器或单一预处理器的潮流算法相比，本发明同样体现出了明显的优势，能够实现电力系统潮流的快速求解，详见下文描述：

一种电力系统潮流并行计算方法，所述方法包括以下步骤：

1)基于牛拉法对电力系统进行潮流计算，获取非线性潮流方程，并依据逐次线性原理得到潮流修正方程组；

2)基于BICGSTAB法求解潮流修正方程组，获取由改进PPAT预处理、改进Jacobi预处理组成的两阶段预处理算法，并对潮流修正方程组的系数矩阵即雅可比矩阵进行预处理，以改善BICGSTAB法的收敛性；

3)基于CPU-GPU异构平台运行步骤2)的并行计算过程。

其中，所述改进PPAT预处理具体为：

借鉴快速解耦潮流算法中有功和无功雅可比矩阵B′、B″组成的分块对角矩阵B可作为雅可比矩阵A的预处理器的特点：

式中，B是A的近似Schur补；设M₀是由A^T的子方阵H^T、L^T组成的分块对角矩阵：

式中，M₀为A的近似Schur补；以M₀作为M_P的初始稀疏模式，完成M_P的构造，即为改进PPAT预处理器；M_P保持了A的类对角占优特性，但其非零元素数量为A的一半；

其中，两阶段预处理的第一阶段具体为：

改进后的PPAT预处理器M_P表示为：

修正方程经改进PPAT预处理后为：

M_PAΔx＝M_Pb

则修正量Δx的系数矩阵A_P为：

其中，M_P1、M_P4为改进PPAT预处理器M_P的对角分块矩阵，A为雅可比矩阵，b为有功无功变化量，A_P1、A_P2、A_P3、A_P4为预处理后新形成的系数矩阵的四个分块矩阵。

其中，所述改进Jacobi预处理具体为：

系数矩阵A_P由四个分块矩阵组成，子矩阵A_P2、A_P3不是方阵但存在满秩子方阵A_P2 ^*、A_P3 ^*，且A_P2 ^*、A_P3 ^*具有类对角占优性；将A_P1、A_P4、A_P2 ^*、A_P3 ^*的主对角线元素取出，并按照其在A_P中的位置，重新组成稀疏矩阵A_J，A_J可描述为四分块矩阵，形如：

式中，S、T^*、U^*、W均为对角方阵，各矩阵对角元素满足：

式中，A_J是高度稀疏的矩阵，其子矩阵S与W为对角矩阵，子矩阵T与U中各行和列的非零元素至多不超过1个；对A_J求逆，得到改进Jacobi预处理器M_J；

其中，两阶段预处理的第二阶段具体为：

修正方程组再经改进Jacobi预处理器后为：

A_PM_Jy＝M_Pb

待求量y的系数矩阵A_PJ为：

A_PJ＝A_PM_J

Δx＝M_Jy

其中，M_J为改进Jacobi预处理。

本发明提供的技术方案的有益效果是：

1、本发明可实现电力系统潮流的准确求解，与传统牛拉法潮流计算结果相比，二者相对误差极小；

2、本发明提出一种改进PPAT预处理与改进Jacobi预处理，并提出一种二者相结合的两阶段预处理方法，使用该预处理方法对雅可比矩阵进行预处理，可有效改善雅可比矩阵特征值分布，改善BICGSTAB法的收敛性和鲁棒性，大幅提高修正方程组的求解效率，进而有效提高潮流计算效率；

3、两阶段预处理的潮流算法与经典的ILU(0)预处理的潮流算法相比，计算速度更快，计算效率更高，鲁棒性更强，且ILU(0)预处理方法的数据相关性强，不利于并行计算，而两阶段预处理方法具有较好的并行特性，适合于借助CPU-GPU异构平台进行加速；

4、本发明所提两阶段预处理的潮流算法与其他常见预处理的潮流算法相比，在规模较大系统中的迭代次数更少，计算耗时更少，且具有更好的适应性，为大规模互联电力系统潮流的快速求解提供了可行思路。

附图说明

图1为一种电力系统潮流并行计算方法的流程图；

图2是本发明算法框架图；

图3为算例误差最大时的PQ节点电压幅值柱形图；

图4为算例计算用时比较柱形图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

实施例1

为了实现电力系统潮流的快速、准确求解，本发明实施例提供了一种电力系统潮流并行计算方法，该方法包括以下步骤：

101：基于牛拉法进行电力系统潮流计算，得到非线性潮流方程，并依据逐次线性原理得到潮流修正方程组；

102：基于BICGSTAB法求解潮流修正方程组，获取由改进PPAT预处理、改进Jacobi预处理组成的两阶段预处理算法，并对潮流修正方程组的系数矩阵即雅可比矩阵进行预处理，以改善BICGSTAB法的收敛性；

103：基于CPU-GPU异构平台实现两阶段预处理BICGSTAB法求解修正方程组的并行化，实现电力系统潮流的快速求解。

综上所述，本发明实施例通过上述步骤101-步骤103，基于CPU-GPU异构平台，将两阶段预处理和BICGSTAB法相结合并行求解修正方程组，提高电力系统潮流计算效率。

实施例2

本发明实施例的总体框架如图2所示，可分为CPU部分和GPU部分，并行加速由GPU部分完成。

其中，GPU部分主要包括：

采用BICGSTAB法求解两阶段预处理后的修正方程组的过程，本发明实施例将该过程称为内层迭代；

其中，CPU部分主要包括：

牛拉法的流程控制，包括修正方程组的形成，预处理器的形成，迭代变量的修正，收敛性判断等过程。本发明实施例将该过程称为外层迭代，其实质为牛顿法潮流计算迭代过程。

下面结合具体的计算公式、附图对实施例1中的方案进行进一步地介绍，详见下文描述：

本发明实施例以下均以线性修正方程组式(1)为例，介绍如何构造预处理器对雅可比矩阵进行预处理。

AΔx＝b (1)

式中，A为雅可比矩阵，Δx为修正量，b为电压幅值与相角组成的列向量。

改进PPAT预处理器基于PPAT预处理器改进得到。

雅可比矩阵A的转置A^T为：

PPAT预处理器旨在构造与A^T具有相同稀疏模式的近似逆。有关PPAT预处理器的构造过程详见文献[10]。由于雅可比矩阵A某些行的非零元素较多，导致所需存储空间较大。为减少数据存储量，本发明对PPAT预处理器稀疏模式进行了适当简化，提出一种改进PPAT预处理器，记为M_P，具体改进如下：

借鉴快速解耦潮流算法中有功和无功雅可比矩阵B′、B″组成的分块对角矩阵B可作为雅可比矩阵A的良好预处理器的特点：

式中，B是A的近似Schur补。设M₀是由A^T的子方阵H^T、L^T组成的分块对角矩阵：

式中，M₀也为A的近似Schur补。以M₀作为M_P的初始稀疏模式，完成M_P的构造，即为改进PPAT预处理器。M_P保持了A的类对角占优特性，但其非零元素数量约为A的一半，减少了数据的存储量和计算量。

改进后的PPAT预处理器M_P可表示为：

修正方程经改进PPAT预处理后为：

M_PAΔx＝M_Pb (6)

则式(6)中修正量Δx的系数矩阵A_P为：

式(1)经过上述改进PPAT预处理后，求解修正方程组转化为求解式(6)所示与式(1)同解的线性方程组。

本发明根据A_P及其子方阵的类对角占优性，对Jacobi预处理器的稀疏模式进行了适当扩充，提出一种改进Jacobi预处理器，记为M_J，具体改进如下：

由式(7)可知，系数矩阵A_P由四个分块矩阵组成，虽然其子矩阵A_P2、A_P3不是方阵，但A_P2、A_P3存在满秩子方阵A_P2 ^*、A_P3 ^*，且A_P2 ^*、A_P3 ^*具有类对角占优性。因此，将A_P1、A_P4、A_P2 ^*、A_P3 ^*的主对角线元素取出，并按照其在A_P中的位置，重新组成稀疏矩阵A_J。A_J可以描述为四分块矩阵，形如：

式中，S、T^*、U^*、W均为对角方阵。各矩阵对角元素满足：

式中，A_J是高度稀疏的矩阵，其子矩阵S与W为对角矩阵，子矩阵T与U中各行和列的非零元素至多不超过1个。对A_J求逆，即得到改进Jacobi预处理器M_J。该运算不会引入任何新的非零元素，故不会引起数据存储量的变化。此外，M_J的构造过程主要涉及高度稀疏矩阵的运算，相对Jacobi预处理器，改进方法增加的计算量极少。

式(6)经改进Jacobi预处理后为：

A_PM_Jy＝M_Pb (10)

则(10)中待求量y的系数矩阵A_PJ为：

A_PJ＝A_PM_J (11)

此时，求解式(11)得到的列向量y并非式(1)所示修正方程组的解Δx，但二者满足：

Δx＝M_Jy (12)

为了简化描述，本发明实施例将上述首先采用改进PPAT预处理器，然后采用改进Jacobi预处理器处理雅可比矩阵的过程称为两阶段预处理，下文统一用两阶段预处理进行描述。

本发明具体步骤如下：

201：读取电力系统的原始节点数据，构建节点导纳矩阵，得到式(1)中的G_ij、B_ij，给定节点电压幅值V和相角δ的初始值；

202：设置迭代次数m＝0，计算精度ε₁，最大迭代次数N₁；

203：将V和δ带入式(1)，判断||f(x)||<ε₁是否成立，若成立，退出牛拉法潮流计算，算法收敛；否则，继续以下步骤；

204：判断m是否达到设定值N₁，如已达到，退出牛拉法潮流计算，算法不收敛；若未达到，令m＝m+1；

205：得到式(1)形式的修正方程组；

其中，步骤205主要包括：

1)极坐标下电力系统节点功率方程为：

式中，V_i、V_j为节点i和j的电压幅值，δ_ij为节点i与j的电压相角差，P_i与Q_i分别为节点i的有功与无功注入量。

2)式(13)可描述为形如F(x)＝0的形式，其中x为状态变量V与δ组成的列向量。将式(13)按照泰勒级数展开，舍弃二阶及以上项得：

式中，

3)式(14)可进一步描述为式(1)所示的形式。

206：采用两阶段预处理BICGSTAB法基于GPU进行修正方程组的求解；

其中，该步骤206包括：

1)若m＝1，构造M_P与M_J，继续步骤2)；否则，跳过此步骤，直接执行步骤2)；

2)使用M_P与M_J对A进行两阶段预处理，将得到的新方程组各部分传入GPU；

3)设置y初值y⁽⁰⁾，允许误差ε₂，求解初始残差r₀＝M_Pb-A_PJy⁽⁰⁾；

4)选择

使得

一般取

并令p₁＝r₀，k＝1；

5)计算

q_k＝r_k-1-α_kA_PJp_k，ω_k＝(q_k,A_PJq_k)/(A_PJq_k,A_PJq_k)；

6)计算y^(k)＝y^(k-1)+α_kp_k+ω_kq_k，r_k＝q_k-ω_kA_PJq_k；

7)若||r_k||₂≤ε₂，迭代收敛，转至步骤9)；否则，继续步骤8)；

8)计算

p_k+1＝r_k+β_k(p_k-ω_kA_PJp_k)，令k＝k+1,转至步骤5)；

9)求解得到Δx^(m)。

207：将GPU计算所得Δx^(m)传回CPU，修正迭代量V和δ，转至步骤203。

综上所述，本发明实施例通过上述步骤201-步骤207以牛拉法为基础，提出一种电力系统潮流并行计算方法，以提升潮流计算速率。

实施例3

下面结合具体的图3、图4、以及表1、表2对实施例1和2中的方案进行可行性验证，详见下文描述：

本实例以不同测试系统进行验证、分析。除特殊说明外，以下算例中牛拉法的计算精度为0.001，求解修正方程组的BICGSTAB法计算精度为0.001。

为检验本发明潮流计算的准确性，分别对比了传统牛拉法与本发明对不同测试系统的潮流计算结果。图3为二者误差最大时的PQ节点电压幅值，表1详细给出了误差最大时两种算法所得PV节点无功功率和节点电压相角的数值及相对误差。显然，由图3可知：对于不同规模的测试系统，本方法计算的PQ节点电压幅值与传统牛拉法基本一致，二者只存在极小误差。进一步，由表1可知：对于不同规模的测试系统，本方法得到的PV节点无功功率和节点电压幅值与传统牛拉法近似相等，二者的最大相对误差均小于0.15％。验证了本发明方法在潮流计算中具有较高的计算精度。

表1误差最大时的无功功率与电压相角

为验证本方法的有效性，将本方法与基于无非零元素填充的ILU(0)预处理的潮流算法进行对比，结果如表2所示。表2中结果为求解修正方程组时采用不同预处理方法的平均迭代次数。

表2使用不同预处理器时潮流算法的迭代次数

注：“/”前后分别为外层牛拉法迭代次数与内层预处理BICGSTAB法求解修正方程组的平均迭代次数

由表2的牛拉法迭代次数对比结果可知：受内部修正方程组求解精度的影响，同一系统迭代次数可能稍有波动，但正常只需迭代3-6次就能收敛。对比预处理BICGSTAB法的平均迭代次数可知：在规模较小系统中，ILU(0)预处理BICGSTAB法迭代次数最少；当计算规模增大到一定程度(如9241节点)后，采用ILU(0)预处理的BICGSTAB法所需迭代次数明显增加，此时改进PPAT预处理和两阶段预处理则能更好地减少BICGSTAB法的迭代次数，而两阶段预处理减少迭代次数的效果最为明显。

值得注意的是：对于单独采用ILU(0)预处理、改进PPAT预处理或改进Jacobi预处理时BICGSTAB法无法收敛的大规模测试系统，如case_ACTIVSg70k，采用两阶段预处理仍可收敛。由此可见：本方法所需迭代次数少，鲁棒性强，为大规模互联电力系统潮流的快速求解提供了可行方案。

为评估本方法计算效率，图4基于CPU平台和CPU-GPU异构平台分别对比了ILU(0)预处理的潮流算法与本方法的计算耗时。对比图4中在同一平台上分别采用ILU(0)预处理和两阶段预处理的潮流算法的计算耗时可知：在规模较小系统中，ILU(0)预处理潮流算法计算时间最短；随着系统规模增大，ILU(0)预处理潮流算法计算耗时急剧增加，而两阶段预处理潮流算法优势逐渐显现，且系统规模越大，优势越显著。比较同一预处理的潮流算法分别在CPU平台和CPU-GPU异构平台上的计算耗时可知：对于规模较小系统，基于CPU-GPU异构平台的并行计算耗时比基于CPU平台的串行计算耗时久，主要原因是GPU并行加速时间难以弥补CPU与GPU间的数据交互时间；当节点规模增加到3000左右时，基于CPU-GPU平台的并行计算优势逐渐显现，且系统规模越大，优势越明显。由图4中基于CPU-GPU异构平台的两阶段预处理潮流算法计算耗时曲线可知：本方法的潮流计算耗时与系统规模近似呈现线性关系，且变化较为平缓，说明本方法具有良好的收敛性及可伸缩性。

为分析不同预处理潮流算法加速情况，以在CPU平台上ILU(0)预处理潮流算法计算耗时为基准值，统计基准值与其它潮流算法计算耗时的比值(记为加速比)，分析不同潮流算法相对于ILU(0)预处理潮流算法的加速情况，结果如表3所示。表3中若CPU平台上ILU(0)预处理潮流算法不收敛，即潮流计算的加速比无基准值时，则仅记录算法收敛信息。

表3不同潮流算法加速比

由表3可知：对于规模较小系统，ILU(0)预处理潮流算法用时最少；随着系统规模扩大，两阶段预处理潮流算法优势逐渐显现。此外，在规模较大系统中，基于CPU-GPU异构平台的ILU(0)预处理潮流算法虽然能一定程度上提高潮流计算效率，但加速效果不佳，甚至与基于CPU平台的两阶段预处理潮流算法相比，仍存在较大差距。

而本方法基于CPU-GPU异构平台的两阶段预处理潮流算法具有最高的计算效率，与基于CPU平台的ILU(0)预处理潮流算法相比，最高加速比可达40.70倍。对ILU(0)预处理潮流算法不收敛的大规模测试系统，本方法依然收敛。由此可见：本方法在求解大规模互联电力系统潮流时具有更高的效率和更强的鲁棒性，能够很好地实现大规模互联电力系统快速潮流计算。

上述分析表明：本方法基于CPU-GPU异构平台，采用两阶段预处理和BICGSTAB法相结合的并行求解方法在较大规模测试系统潮流计算中，可有效提高潮流计算效率，且具有良好的算法稳定性，较传统潮流计算方法的优势更为明显。

参考文献

[1]唐坤杰，董树锋，宋永华.基于不完全LU分解预处理迭代法的电力系统潮流算法[J].中国电机工程学报，2017，37(S1)：55-62.

[2]蔡大用，陈玉荣.用不完全LU分解预处理的不精确潮流计算方法[J].电力系统自动化，2002，26(8)：11-14.

[3]张永杰，孙秦.大型稀疏线性方程组符号LU分解法[J].计算机工程与应用，2007，43(28)：29-30.

[4]刘洋，周家启，谢开贵，等.基于Beowulf集群的大规模电力系统方程并行PCG求解[J].电工技术学报，2006，21(3)：105-111.

[5]Li X，Li F.GPU-based two-step preconditioning for conjugategradient method in power flow[C].IEEE Power&Energy Society General Meeting，Denver，CO，USA 2015.

[6]Li X，Li F，Clark J M.Exploration of multifrontal method with GPU inpower flow computation[C].IEEE Power&Energy Society General Meeting.2013：1-5.

[7]Zhou G，Bo R，Chien L，et al.GPU-accelerated algorithm for on-lineprobabilistic power flow[J].IEEE Transactions on Power Systems，2017，PP(99)：1-1.

[8]Benzi M.Preconditioning techniques for large linear systems：asurvey[J].Journal of Computational Physics，2002，182(2)：418-477.

[9]Yu Z，Huang S，Shi L，et al.GPU-based JFNG method for power systemtransient dynamic simulation[C].International Conference on Power SystemTechnology.IEEE，2014.

[10]迟利华，刘杰，李晓梅.稀疏近似逆预条件子及其并行计算[J].计算机学报，2000，23(3)：255-260.

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种电力系统潮流并行计算方法，其特征在于，所述方法包括以下步骤：

1)基于牛拉法对电力系统进行潮流计算，获取非线性潮流方程，并依据逐次线性原理得到潮流修正方程组；

2)基于BICGSTAB法求解潮流修正方程组，获取由改进PPAT预处理、改进Jacobi预处理组成的两阶段预处理算法，并对潮流修正方程组的系数矩阵即雅可比矩阵进行预处理，以改善BICGSTAB法的收敛性；

3)基于CPU-GPU异构平台运行步骤2)的并行计算过程。

2.根据权利要求1所述的一种电力系统潮流并行计算方法，其特征在于，所述改进PPAT预处理具体为：

借鉴快速解耦潮流算法中有功和无功雅可比矩阵B′、B″组成的分块对角矩阵B可作为雅可比矩阵A的预处理器的特点：

式中，B是A的近似Schur补；设M₀是由A^T的子方阵H^T、L^T组成的分块对角矩阵：

式中，M₀为A的近似Schur补；以M₀作为M_P的初始稀疏模式，完成M_P的构造，即为改进PPAT预处理器；M_P保持了A的类对角占优特性，但其非零元素数量为A的一半；

则两阶段预处理的第一阶段具体为：

改进后的PPAT预处理器M_P表示为：

修正方程经改进PPAT预处理后为：

M_PAΔx＝M_Pb

则修正量Δx的系数矩阵A_P为：

其中，M_P1、M_P4为改进PPAT预处理器M_P的对角分块矩阵，A为雅可比矩阵，b为有功无功变化量，A_P1、A_P2、A_P3、A_P4为预处理后新形成的系数矩阵的四个分块矩阵。

3.根据权利要求2所述的一种电力系统潮流并行计算方法，其特征在于，所述改进Jacobi预处理具体为：

系数矩阵A_P由四个分块矩阵组成，子矩阵A_P2、A_P3不是方阵但存在满秩子方阵A_P2 ^*、A_P3 ^*，且A_P2 ^*、A_P3 ^*具有类对角占优性；将A_P1、A_P4、A_P2 ^*、A_P3 ^*的主对角线元素取出，并按照其在A_P中的位置，重新组成稀疏矩阵A_J，A_J可描述为四分块矩阵，形如：

式中，S、T^*、U^*、W均为对角方阵，各矩阵对角元素满足：

式中，A_J是高度稀疏的矩阵，其子矩阵S与W为对角矩阵，子矩阵T与U中各行和列的非零元素至多不超过1个；对A_J求逆，得到改进Jacobi预处理器M_J；

则两阶段预处理的第二阶段具体为：

修正方程组再经改进Jacobi预处理器后为：

A_PM_Jy＝M_Pb

待求量y的系数矩阵A_PJ为：

A_PJ＝A_PM_J

Δx＝M_Jy

其中，M_J为改进Jacobi预处理。

Images