CN109447839B - 非零元素的随机生成及形成电力系统雅可比矩阵的方法 - Google Patents

Abstract

非零元素的随机生成及形成电力系统雅可比矩阵的方法，涉及电力系统分析计算领域，主要包括以下步骤：建立按随机顺序分区存储Y阵上三角非零元素的Y(n，d′)数组；根据线路支路数据i、j、r、x、k，分别计算节点的自导纳和互导纳，再根据列号j与m的大小关系累加S_1i或S_2i并按随机顺序存放在PQ区或PV区；根据S_1i、S_2i的数值将二维Y(n，d′)数组中的空置单元剔除，以一维紧凑格式写入Y(n，d′)的数据文件以备调用；将Y(n，d′)数据文件的数据写入Y(n，d)数组；根据Y(n，d)数组的参数按对称性随机分段计算P_i、Q_i；根据Y阵元素与J_ij、J_ji子阵非零元素位置的对应关系以及J_ij、J_ji子阵元素本身和各自的数值对应关系，按(两行+两列)/次方式对称形成J阵。本方法各个环节的速度均优于传统法。

Description

非零元素的随机生成及形成电力系统雅可比矩阵的方法

技术领域

本发明属于电力系统分析计算领域，涉及一种非零元素的随机生成及形成电力系统雅可比矩阵的方法。

背景技术

电力系统计算中广泛应用的节点导纳矩阵Y和雅可比矩阵J都是极度稀疏矩阵，其中Y阵对称而J阵不对称。但如果将J阵元素结构用J_ij子阵表示，则除平衡节点外，Y阵元素结构与J阵中J_ij子阵结构完全相同。此外，虽然J阵不对称，但J阵中J_ij子阵与J_ji子阵的非零位置几乎对称，而这个特点使得Y阵元素与J阵中的J_ij子阵及J_ji子阵元素关系紧密，然而传统方法形成J阵时未能利用上述关系，从而导致J阵形成时间过长。此外，传统法形成J阵的过程中也未利用J_ij子阵与J_ji子阵各个元素之间数值上的对应关系简化J阵对角元素和非对角元素的计算，因此其计算过程繁琐，并有大量的重复计算，从而导致对J阵元素的计算效率低下。

设系统节点数为n，PQ节点数为m。传统方法中不考虑Y元素稀疏性时的Y(n，2n)存贮方式形式简单直观，方便对Y阵数据的处理，但其大量零元素的存贮需要较大的存贮单元，导致读写Y(n，2n)数据文件的时间较长。Y(n，2n)数组中元素结构与J阵中元素结构相似，利用Y(n，2n)数组可以方便地形成J阵，但由于未利用Y阵元素与J_ij子阵及J_ji子阵元素的关系以及J_ij子阵与J_ji子阵元素之间的关系，其形成J阵效率不高。

传统方法中考虑Y阵元素稀疏性的坐标法、顺序法、链表法等虽贮存单元大大减少，但由于其对角元素与非零的非对角元素分开存放，不但存贮结构复杂，也不利于数据的检索、修改、应用，且存贮方式与Y阵、J阵的元素结构没有明确的对应关系，更没有对Y阵元素进行分区存储，因此无法体现Y阵元素与J_ij子阵及J_ji子阵元素的关系以及J_ij子阵与J_ji子阵元素之间的关系，从而也不能利用上述关系，导致其数据文件的读写速度无法达到最佳，且形成J阵的速度也不理想。此外，传统方法中无论考虑稀疏性与否，其节点有功功率P_i和无功功率Q_i的计算效率也都不高。

文献①(一种基于对称稀疏矩阵技术快速求取极坐标牛顿-拉夫逊法潮流的方法：中国，[ZL 201510770899.X])提出根据Y(n，d)存贮方式利用对称稀疏矩阵技术快速形成J阵并消元计算潮流，虽给出了Y阵上三角非零元素的存贮方式，但并未给出具体的计算形成方法，且其形成Y阵的过程中j₁＜j₂---＜j₆的要求需应用大量循环和判断语句，同样导致计算效率不高；且其形成Y阵过程中未将PQ节点与PV节点的Y阵数据分区存储，导致后续利用形成J阵时需要大量判断语句，同样降低计算效率；且在形成J阵时虽利用了J_ij、J_ji子阵非零元素位置的对应关系，但也未利用极坐标牛顿法中各元素数值上的对应关系，同样大大降低了J阵对角元素和非对角元素的计算效率。

发明内容

为克服现有技术的不足，本发明方法提出一种非零元素的随机生成及形成电力系统雅可比矩阵的方法。本申请包括以下步骤：

步骤1：定义按随机顺序分区存储Y阵上三角非零元素的Y(n，d′)数组；

(1)对Y(n，d′)数组按随机顺序仅存储Y阵上三角的非零元素，并将其分为非零计数组、对角元素组、非对角元素组，其中非对角元素组再分为PQ区和PV区；

(2)非零计数组位于第1～2列，分别存放每行对角元以右非零的非对角元素总数S_1i以及对角元以右与PQ节点对应的非零的非对角元素数S_2i；

形成Y阵时S_1i、S_2i均由程序自动累加产生，S_1i可保证准确快速地读写Y阵参数、使数据文件的存储单元最小、数据读写效率最高，并可在后续形成J阵或潮流计算过程大大加快P_i、Q_i的计算；S_2i的作用是可根据J阵中PQ节点与PV节点元素结构的不同以及对应的J_ij子阵及J_ji子阵元素的对应关系不同时，更高效地计算J阵元素，并省去大量判断语句。

(3)对角元素组位于第3～5列，存放对角元行号i和自导纳g_ii、b_ii；

(4)非对角元素组位于第6-80列，再按PQ区和PV区均三列一小组按随机顺序存放与相应对角元素连接的上三角非零的非对角元素的列号j及参数；当j≤m，将与PQ节点对应的列号和互导纳j、g_ij、b_ij存放在第6-50列的PQ区；当j＞m，将与PV和平衡节点对应的列号和互导纳j′、g′_ij、b′_ij存放在第51-80列的PV区；

1)设系统中1～m为PQ节点，m+1～n-1为PV节点，第n节点数为平衡节点。根据系统中上三角的静态非零元素数最大一般不超过15个的特点，再考虑PQ和PV节点之间非零的非对角元素的分布比例，可假设PQ节点非对角元素区域排列在Y(n，d′)数组的第6-50列，PV节点非对角元素区域排列在Y(n，d′)数组的第51-80列，因此d′＝80。

2)非对角元素组中不要求j₁＜j₂＜--＜j_S2i、j′₁＜j′₂＜---＜j′_S3i，即所有的非零元素可以随机存放，但要求i＜j、i＜j′、.j＜j′，即要求按PQ区和PV区仅存放上三角元素。

步骤2：根据线路支路数据i、j、r、x、k，分别计算i、j节点的自导纳和互导纳，再根据列号.j与m的大小关系累加S_1i或S_2i并按随机顺序相继存放在PQ区或PV区；

(1)将初次计算的i、j节点的行号和自导纳g′_ii、b′_ii、g′_jj、b′_jj分别存放在第i行和第j行的对角元组，并将后续计算的i、j节点自导纳分别与g′_ii、b′_ii或g′_jj、b′_jj累加；

(2)若i＜j，则为上三角元素；若列号j≤m，用S_2i统计PQ区的非零元素数，即将第i行S_2i加1，将其列号j与互导纳g_ij、b_ij存放在第i行起始列号为T＝3*(S_2i+1)的非对角元素组的PQ区，即分别存放j、g_ij、b_ij在第T列～第T+2列；若列号j＞m，用S_3i统计PV区的非零元素数，即将第i行的S_3i加1，将其列号j作为j′与互导纳g′_ij、b′_ij存放在第i行起始列号为T′＝48+3*S_3i的非对角元素组的PV区，即分别存放.j′、g′_ij、b′_ij在第T′列～第T′+2列；

(3)若i＞j，则为下三角元素，将其行列号对调；若对调后的列号i≤m，则将第j行的S_2j加1，将其列号i与互导纳g_ji、b_ji存放在第j行起始列号为T＝3*(S_2j+1)的非对角元素组的PQ区，即分别存放i、g_ji、b_ji在第T列～第T+2列；若对调后的列号i＞m，将第j行的S_3j加1，将其列号i作为i′与互导纳g′_ji、b′_ji存放在第j行起始列号为T′＝48+3*S_3i的非对角元素组的PV区，即分别存放i′、g′_ji、b′_ji在第T′列～第T′+2列；

(4)依此循环直至建立仅含Y阵上三角非零元素的Y(n，d′)数组，其中S_1i＝S_2i+S_3i。

由于支路数据的随机读取，因此无论在Y(n，d′)数组中非对角元素组的PQ区还是PV区，其存放的列号及参数也是随机顺序。这种存放方式可省去循环语句中的大量判断，从而大大加快Y阵的形成速度，而Y阵元素的随机存放不但不影响后续J阵的形成，同样可以利用J_ij和J_ji子阵非零元素的对应关系快速形成J阵。

Y(n，d′)数组的二维结构如表1。

表1 Y(n，d′)数组元素存放结构

注意：在Y(n，d′)数组PQ区各行第50列和PV区各行第80列以左的部分单元中几乎均无数据，这些空置单元似乎浪费了Y(n，d′)数组的容量，但实际上将Y((n，d′)数组写入Y(n，d′)数据文件时由于S_1i、S_2i的控制作用，这些空置单元不但不会影响Y(n，d′)数组的形成速度，而且使得所有写入Y(n，d′)数据文件中的Y(n，d′)数组的空置单元也将全部被剔除，因而这些空置单元无论对Y(n，d′)数组还是对Y(n，d′)的数据文件几乎均无影响。

步骤3：根据S_1i、S_2i的数值将二维Y(n，d′)数组中的空置单元剔除，以一维紧凑格式写入Y(n，d′)的数据文件以备调用；

步骤4：将Y(n，d′)数据文件的数据写入Y(n，d)数组；

(1)将一维紧凑格式的Y(n，d′)数据文件的数据写入二维Y(n，d)数组中，每行对角元以右非零的非对角元素总数S_1i以及对角元以右与PQ节点对应的非零的非对角元素数S_2i与Y(n，d′)数组完全相同，再假设与S_1i对应的系统最大非零元素数为S_1max，则d＝3S_1max+5；

(2)Y(n，d′)数组中的最大存储单元为n×d′＝80n，但由于其PQ区和PV区的结构特点，每行的实际列数均为50+3S_3i，因此Y(n，d′)的实际存储单元为n×∑(50+3S_3i)，其PQ区中仍有不少空置单元。Y(n，d)数组的最大存储单元为n×d＝n×(3S_1max+5)≤50n，远小于Y(n，d′)数组中80n的最大存储单元；每行的实际列数为3S_1i+5，因此Y(n，d)的实际存储单元为n×∑(3S_1i+5)，其PQ区无任何空置单元，也远小于Y(n，d′)数组中的实际存储单元。

(3)由于Y(n，d)数组中的PQ区无空置单元，其各行的PV区数据均紧接着相应的PQ区数据，因此Y(n，d)数组中各行PQ区和PV区的分界不在同一列上，但由于S_1i、S_2i控制，完全不影响对Y(n，d)数组数据的应用；

步骤5：根据Y(n，d)数组的参数按对称性随机分段计算P_i、Q_i；

传统方法中极坐标P_i、Q_i的计算式如下：

式(1)中P_i、Q_i的计算需应用Y阵第i行的一整行元素一次性完成计算，由于本发明方法中是随机存贮Y阵对角元及上三角的非零元素，因此必须利用Y阵元素的对称性随机分段计算P_i、Q_i，为此可将式(1)改写成式(2)。

式(2)将P_i、Q_i的计算分解成的对n个节点功率ΔP_i1～ΔP_in、ΔQ_i1～ΔQ_in增量的计算，其中ΔP_ij、ΔQ_ij分别表示P_i、Q_i求和项中的第j项，而ΔP′_ij＝V_j(g_ijcosδ_ij+b_ijsinδ_ij)、ΔQ′_ij＝V_j(g_ijsinδ_ij-b_ijcosδ_ij)。同样可按此计算ΔP_ii、ΔP_ik、ΔQ_ii、ΔQ_ik等。

假设Y阵第i行只有j、i、k、m四列元素非零，且j＜i和i＜k＜m，即第j列位于下三角，而第i、k、m列均位于上三角，此时可将式(2)再分解并简化成式(3)。

根据式(3)计算第i行的P_i时，其ΔP_ij增量需根据Y阵的下三角元素进行计算，而ΔP_ij只能在前面计算第.j行的ΔP_ji时，同时根据对称性计算第i行的ΔP_ij。因此计算第i行的P_i时，只需累加前面计算的ΔP_ij，而根据Y(n，d)数组直接计算上三角部分的功率增量ΔP_ii、ΔP_ik、ΔP_im，然后两者叠加则可得到P_i。同理，在根据Y(n，d)数组计算第i行上三角部分的功率增量ΔP_ii、ΔP_ik、ΔP_im时，同样也要根据对称性计算第k、m行的功率增量ΔP_ki、ΔP_mi，以便在计算第k、m行的P_i时直接进行累加。Q_i计算过程同上。

步骤6：根据Y阵元素与J_ij、J_ji子阵非零元素位置的对应关系以及J_ij、J_ji子阵元素本身和各自的数值对应关系，用Y(n，d)数组按(两行+两列)/次方式对称形成J阵；

(1)利用Y阵元素与J_ij、J_ji子阵非零元素静态位置的对应关系快速形成J阵；

将极坐标牛顿法修正方程式(4)按虚线分为A、B、C、D区，相关子阵标于上。

Y阵中如b_ij≠0可得y_ij≠0；根据对称性可得b_ji≠0和y_ji≠0；再根据Y阵元素结构与J_ij、J_ji子阵非零元素位置的对应关系，可推出子阵J_ij≠0、J_ji≠0；继而可推出J_ij、J_ji子阵中的元素也均不非零。式(4)中J阵中元素的分区计算如下：

1)A区的J_1m与J_m1中，b_1m≠0→J_1m≠0→J_m1≠0→(H_1m、N_1m、M_1m、L_1m)≠0和(H_m1、N_m1、M_m1、L_m1)≠0，即根据Y阵一个元素b_1m可计算J阵八个元素，为“据1计8”；

2)B区的J_1，m+1与C区的J_m+1，1中，b_1，m+1≠0→J_1，m+1≠0→J_m+1，1≠0→(H_1，m+1、M_1，m+1)≠0和(H_m+1，1、N_m+1，1)≠0，即根据Y阵一个元素b_1，m+1可计算J阵四个元素，为“据1计4”；

3)D区的J_m+1，n-1与J_n-1，m+1中，b_m+1，n-1≠0→J_m+1，n-1≠0→J_n-1，m+1≠0→H_n+1，n-1≠0和H_n-1，m+1≠0，即根据Y阵一个元素b_m+1，n-1可计算J阵二个元素，为“据1计2”；

由于Y(n，d)数组中元素的随机存放，因此可根据其上三角的一个非零元素，以J阵中的四个对角元为起点，按(两行+两列)/次对称计算方式，分别计算J阵中的A、B和C、D四个区的8个或4个或2个元素，从而以随机顺序快速形成J阵。

(2)利用J_ij、J_ji子阵元素本身和各自数值对应关系快速完成非对角元素和对角元素的计算；

1)利用J_ij、J_ji子阵非对角元素本身的数值对应关系计算，如J_1m与J_m1中，利用M_1m＝-N_1m、L_1m＝H_1m；

2)利用J_ij、J_ji子阵非对角元素各自的数值对应关系计算；

表面上J_ij与J_ji子阵中元素的数值不等，如H_ij与H_ji、N_ij与M_ji、M_ij与M_ji、L_ji与L_ji，因此传统方法中H_ij与H_ji、N_ij与N_ji、M_ij与M_ji、L_ij与L_ij元素均分开计算，但将各元素计算式重写如下后可以发现以下对应关系。

可以看出：虽然有H_ij≠H_ji，N_ij≠N_ji等，但将其各自分解后可得H_ij＝c₁+d₁和H_ji＝-c₁+d₁、N_ij＝c₂-d₂和H_ji＝c₂+d₂、M_ij＝-(c₂-d₂)＝-N_ij和M_ji＝-(c₂+d₂)＝-N_ji、L_ij＝c₁-d₁＝H_ij和L_ji＝-c₁-d₁＝H_ji。因此在计算H_ij时可同时计算L_ij、H_ji、L_ji，比各自分别计算效率更高。

3)利用J_ii子阵对角元素各自的数值对应关系计算；

传统方法中J阵对角元的计算未利用P_i、Q_i的计算，其计算式中∑的重复计算大大影响计算效率，且所有的对角元素H_ii、N_ii、M_ii、L_ii均分开计算。但将对角元素的计算公式重写如下后也可以发现以下对应关系。

上述等式表明，所有对角元素H_ii、N_ii、M_ii、L_ii的计算均可分为按非对角元素计算和用P_i、Q_i进行修正两步完成，而在未用P_i、Q_i修正时，其以非对角元素计算的部分有与非对角元素同样的特点，即H′_ii＝L′_ii、M′_ii＝-N′_ii。因此在计算非对角元素时可同时根据该特点计算其对角元素的相应部分，然后用P_i或Q_i修正，以此简化并加快对角元素的计算。

本发明方法的有益效果为：

1.采用上三角非零元素Y(n，d)的新型随机存放方式，并分别根据不同的非零元素计数可大大减少Y阵的存贮单元和形成Y阵时的判断语句、提高Y阵的形成速度及Y阵数据文件的读写速度，便于数据的计算、检索、修改等应用；其节点数据按PQ与PV分区存放，可省去形成J阵时的大量判断语句、大大加快形成J阵速度。

2.根据PQ与PV分区、利用Y阵元素与J_ij与J_ji子阵非零元素位置对应的关系，形成J阵时按“据1计8”、“据1计4”、“据1计2”规则分区计算，可按(两行+两列)/次方式对称快速形成J阵。

3.利用J_ii子阵、J_ij与J_ji子阵元素本身和各自数值的对应关系快速对称完成非对角元素和对角元素的计算；

4.利用Y阵元素的对称性随机分段计算P_i、Q_i的新方法，进一步加快P_i、Q_i和J阵对角元的计算。

附图说明

图1为传统方法不考虑元素稀疏性及对称性形成Y(n，2n)数组流程图

图2为本发明方法考虑元素稀疏性及对称性形成Y(n，d′)数组流程图

图3为传统方法用Y(n，2n)数组形成J阵流程图。

图4为本发明方法考虑元素稀疏性及对称性形成Y(n，d′)数组并用Y(n，d)数组形成J阵流程图。

具体实施方式

本申请方法将通过以下实施举例作进一步说明。

实施例1对IEEE-30～-118系统，分别比较用传统法形成Y(n，2n)数组和本发明方法形成Y(n，d′)数组并读写Y(n，2n)、Y(n，d′)数据文件、用Y(n，2n)数组和本发明方法的Y(n，d)数组形成J阵的时间，比较结果如表2所示。

表2传统法和本方法形成并读写Y阵数据文件、形成J阵时间的比较

t_c.w、t_n.w：传统法、本方法形成及写Y阵数据文件时间。

t_n.w/t_c.w：本方法与传统法形成及写Y阵数据文件时间百分比。

t_c.r、t_n.r：传统法、本方法读Y阵数据文件时间。

t_n.r/t_c.r：本方法与传统法读Y阵数据文件时间百分比。

t_c.fj、t_n.fj：传统法、本方法形成J阵时间。

t_n.fj/t_c.fj：本方法与传统法形成J阵时间百分比。

根据表2可以看出：

1.本方法形成及读写Y阵数据文件、形成J阵时间均大大优于传统法。

以IEEE-118系统为例，本方法形成及写Y阵数据文件时间为传统法的15.98％；读Y阵数据文件时间为传统法的11.11％；不考虑读数据文件时，形成J阵时间为传统法的8.52％。

2.电力系统节点数越大，本方法读写数据文件和形成J阵的速度越快。

本方法可以采用任何一种编程语言和编程环境实现。本发明中采用C++编程语言，开发环境是Visual C++，使用电脑型号为联想天启M4500。

Claims (1)

1.非零元素的随机生成及形成电力系统雅可比矩阵的方法，其特征是，包括以下步骤：

步骤1：定义按随机顺序分区存储Y阵上三角非零元素的Y(n，d′)数组；

(1)对Y(n，d′)数组按随机顺序仅存储Y阵上三角的非零元素，并将其分为非零计数组、对角元素组、非对角元素组，其中非对角元素组再分为PQ区和PV区；

(2)非零计数组位于第1～2列，分别存放每行对角元以右非零的非对角元素总数S_1i以及对角元以右与PQ节点对应的非零的非对角元素数S_2i；

(3)对角元素组位于第3～5列，存放对角元行号i和自导纳g_ii、b_ii；

(4)非对角元素组位于第6-80列，再按PQ区和PV区均三列一小组按随机顺序存放与相应对角元素连接的上三角非零的非对角元素的列号j及参数；当j≤m，将与PQ节点对应的列号和互导纳j、g_ij、b_ij存放在第6-50列的PQ区；当j＞m，将与PV和平衡节点对应的列号和互导纳j′、g′_ij、b′_ij存放在第51-80列的PV区；

步骤2：根据线路支路数据i、j、r、x、k，分别计算i、j节点的自导纳和互导纳，再根据列号j与m的大小关系累加S_1i或S_2i并按随机顺序相继存放在PQ区或PV区；

(1)将初次计算的i、j节点的行号和自导纳g′_ii、b′_ii、g′_jj、b′_jj分别存放在第i行和第j行的对角元组，并将后续计算的i、j节点自导纳分别与g′_ii、b′_ii或g′_jj、b′_jj累加；

(2)若i＜j，则为上三角元素；若列号j≤m，则将第i行S_2i加1，将其列号j与互导纳g_ij、b_ij存放在第i行起始列号为T＝3*(S_2i+1)的非对角元素组的PQ区；若列号j＞m，则将第i行的S_3i加1，将其列号j作为j′与互导纳g′_ij、b′_ij存放在第i行起始列号为T′＝48+3*S_3i的非对角元素组的PV区；

(3)若i＞j，则为下三角元素，将其行列号对调；若对调后的列号i≤m，则将第j行的S_2j加1，将其列号i与互导纳g_ji、b_ji存放在第j行起始列号为T＝3*(S_2j+1)的非对角元素组的PQ区；若对调后的列号i＞m，则将第j行的S_3j加1，将其列号i作为i′与互导纳g′_ji、b′_ji存放在第j行起始列号为T′＝48+3*S_3i的非对角元素组的PV区；

(4)依此循环直至仅含Y阵上三角非零元素的Y(n，d′)数组建立，其中S_1i＝S_2i+S_3i；

步骤3：根据S_1i、S_2i的数值将二维Y(n，d′)数组中的空置单元剔除，以一维紧凑格式写入Y(n，d′)的数据文件以备调用；

步骤4：将Y(n，d′)数据文件的数据写入Y(n，d)数组；

将一维紧凑格式的Y(n，d′)数据文件的数据写入二维Y(n，d)数组中，每行对角元以右非零的非对角元素总数S_1i以及对角元以右与PQ节点对应的非零的非对角元素数S_2i与Y(n，d′)数组完全相同；

步骤5：根据Y(n，d)数组的参数按对称性随机分段计算P_i、Q_i；

将传统方法中极坐标P_i、Q_i的计算分解成的对n个节点功率增量的计算，其中ΔP′_ij＝V_j(g_ijcosδ_ij+b_ijsinδ_ij)、ΔQ′_ij＝V_j(g_ijsinδ_ij-b_ijcosδ_ij)，假设Y阵第i行只有j、i、k、m四列元素非零，且j＜i和i＜k＜m，即第j列位于下三角，而第i、k、m列均位于上三角，此时可将P_i、Q_i的计算简化成下式；

上式中ΔP_ij增量需根据Y阵的下三角元素计算，且只能在前面计算第j行的ΔP_ji时，同时根据对称性计算第i行的ΔP_ij；因此直接计算上三角部分的功率增量ΔP_ii、ΔP_ik、ΔP_im，然后再累加前面计算的ΔP_ij即可得到P_i；

步骤6：根据Y阵元素与J_ij、J_ji子阵非零元素位置的对应关系以及J_ij、J_ji子阵元素本身和各自的数值对应关系，用Y(n，d)数组按(两行+两列)/次方式对称形成J阵；

(1)利用Y阵元素与J_ij、J_ji子阵非零元素的对应关系计算J阵元素；

1)对A区的J_1m与J_m1，根据b_1m≠0可计算H_1m、N_1m、M_1m、L_1m和H_m1、N_m1、M_m1、L_m1八个元素；

2)对B区的J_1，m+1与C区的J_m+1，1，根据b_1，m+1≠0可计算H_1，m+1、M_1，m+1和H_m+1，1、N_m+1，1四个元素；

3)对D区的J_m+1，n-1与J_n-1，m+1，根据b_m+1，n-1≠0可计算H_m+1，n-1和H_n-1，m+1二个元素；

(2)利用J_ij、J_ji子阵元素本身和各自数值对应关系计算J阵元素；

1)利用J_ij、J_ji子阵非对角元素本身M_1m＝-N_1m、L_1m＝H_1m的对应关系计算；

2)利用J_ij、J_ji子阵非对角元素H_ij＝c₁+d₁和H_ji＝-c₁+d₁、N_ij＝c₂-d₂和N_ji＝c₂+d₂、M_ij＝-(c₂-d₂)＝-N_ij和M_ji＝-(c₂+d₂)＝-N_ji、L_ij＝c₁-d₁＝H_ij和L_ji＝-c₁-d₁＝H_ji的对应关系计算；

3)利用J_ii子阵对角元素未用P_i、Q_i修正时，H′_ii＝L′_ii、M′_ii＝-N′_ii的对应关系计算。

Images