电力系统稳态分析牛顿拉夫逊法

1、0引言潮流是配电网络分析的基础，用于电网调度、运行分析、操作模拟和设计 规划，同时也是电压优化和网络接线变化所要参考的内容。潮流计算通过数值 仿真的方法把电力系统的详细运行情况呈现给工作人员，从而便于研究系统在 给定条件下的稳态运行特点。随着市场经济的发展，经济利益是企业十分看重 的，而线损却是现阶段阻碍企业提高效益的一大因素。及时、准确的潮流计算 结果，可以给出配电网的潮流分布、理论线损及其在网络中的分布，从而为配 电网的安全经济运行提供参考。从数学的角度来看，牛顿-拉夫逊法能有效进行非线性代数方程组的计算且具有二次收敛的特点，具有收敛快、精度高的特 点，在输电网中得到广泛应用。随着现代计算

2、机技术的发展，利用编程和相关 软件，可以更好、更快地实现配电网功能，本文就是结合牛顿-拉夫逊法的基本原理，利用C+程序进行潮流计算，计算结果表明该方法具有良好的收敛性、 可靠性及正确性。1牛顿-拉夫逊法基本介绍1.1潮流方程对于N个节点的电力网络（地作为参考节点不包括在内），如果网络结构 和元件参数已知，则网络方程可表示为：Y&1&（ 1-1）式中，Y为N*N阶节点导纳矩阵；V为N*1维节点电压列向量；I&为N*1维节 点注入电流列向量。如果不计网络元件的非线性，也不考虑移相变压器，则Y为对称矩阵。电力系统计算中，给定的运行变量是节点注入功率，而不是节点注入电 流，这两者之间有如下关系：E&

3、S&（ 1-2）式中，S&为节点的注入复功率，是N*1维列矢量；S&为S&的共轭；E diag V&是由节点电压的共轭组成的N*N阶对角线矩阵。由(1-1 )和(1-2),可得：s&EY&上式就是潮流方程的复数形式，是N维的非线性复数代数方程组。将其展开，有:pjQiV& YjV&j=1,2 ,.,N(1-3)j i式中，j i表示所有和i相连的节点j ,包括j i o将节点电压用极坐标表示，即令V& V i ,代入式(1-3)中则有：P jQiVi Gj ijBij V jV V Gj ijBijcos ijj sin ij故有：P V j iV G cos ijBjsin ji=1,2,N

4、(1-4)Q V V G sin ijBjcos ijj i式(1-4)是用极坐标表示的潮流方程。而节点功率误差：P PSP Vi Vj(Gj cos j Bj sin J (1-5)j iQi QiSP VVj(Gij cos ij Bq sin ij)(1-6)j i式中：PiSP，QiSP为节点i给定的有功功率及无功功率。1.2牛顿-拉夫逊法基本原理牛拉法的一般描述牛拉法是把非线性方程式的求解过程变成反复对相应的线性方程式的求解 过程，即非线性问题通过线性化逐步近似，这就是牛拉法的核心。下面以非线 性方程式的求解过程来进行说明。设电力网络的节点功率方程一般形式如下:SP/、y y x （

5、 1-7）式中，ySP为节点注入功率给定值；y为ySP对应的物理量和节点电压之间的函数表达式；x为节点电压。写成功率偏差的形式：SPf x y y x 0（ 1-8）应用牛拉法求解如下。在给定的初值x 0处将式（1-8）作一阶泰勒展开:定义J耳为潮流方程的雅克比矩阵,xJ。为J在x 0处的值，则有:xJ01f x用x修正x 0就得到x的新值。如果用表示迭代次数，写成一般的表达式,有:k 1 f x k（1-9）对于潮流收敛的情况,应比xk更接近于解点。收敛条件为:maxf i x k由简单迭代法收敛性分析的结论知，越接近解点，牛顿-拉夫逊法收敛越快，它具有二阶收敛速度。由图1.1可以直观地了解

6、牛拉法的步骤:图1.1牛顿-拉夫逊法的几何解释极坐标的牛顿-拉夫逊法在极坐标中，f x有如下的形式:,V,VpSPQSP,V,V(1-10)共2n-r个方程，状态变量为:T vtn V1 V2 L Vn r共2n-r个待求量。r个PV节点的电压幅值给定,不需求解。潮流雅克比矩阵的维数是(2n-r) *(2n-r)，结构如下:npV7_QV7 n r1，为使雅克比上式右侧的对电压幅值的偏导数项中的电压幅值的阶数减少了 矩阵的各部分子矩阵具有一致的形式，在实际计算中，常将该项乘以电压幅T值，并选取 V / VV,/ V,V2/ V2LVnr/ Vnr作为待求的修正量，则雅克比矩阵可写成:将式（1-

7、10）和（1-11）代入式它修正x直到max.将式（1-11）用下式表示:其中每个字块的计算公式如下:Hh对角元素:非对角元素:NiiMiHijNjMjLj二 VQAVT(1-11)(1-9)为止。的修正方程即可求得x的修正量x，用Bii卫VVQjQ-V.V j2牛顿法潮流计算步骤2.1程序流程图vnM nVMM MiiVLMV Hj Vj, HjVNjVj,NjVMjV，MjVLjV IBiBijQ疋P才P疋Qcos ijGj cosNjHjij(1-12)sinjBj sinij(1-13)在了解了牛拉法的原理之后，明确程序编写思路，如图2.1、2.2所示。其中图2.1中的“计算电压幅值和

8、角度”步骤较多，单独用图2.2表示出来。图2.1牛顿法计算潮流的程序框图图22电压幅值和角度求解步骤框图当不符合收敛的条件“ amontk1 ”时，即认为计算不收敛。具体程序见附录。2.2计算步骤下面讨论的是极坐标形式的牛顿法求解过程，大致分为以下几个步骤： 形成节点导纳矩阵； 给各节点电压设初值（V（0）, i（0）； 根据式（1-12）、（ 1-13）生成雅克比矩阵（H、N、M、L）; 将节点电压初值代入式（1-5）、式（1-6），求出修正方程式的常数项向量P, Qi ； 求解修正方程，得到电压幅值和角度； 判断是否收敛，若收敛，计算平衡节点和线路功率； 输出结果，并结束。3算例3.1系统

9、模型本文以图3.1所示电力网络为例，调用基于牛顿-拉夫逊法的C+程序其中节点4设为平衡节点，电压标幺值为1.05,计算误差为0.0000013.2输入与输出将图3.1所示模型的相关数据放在data.dat文件中$ Efibcsgft /isual C+ + - d牛就拉夫逊搖日上曰文他B 輛吒)童看GO捶入工程0 會 IMffi 京口血 m(H匕 Qy |叮|虚错|恤|JIdl丄1料1|1.05 &. 90999101 1 2 9.1 0.4 氛 91528? 1 4 9.12 0.5 0.017203 2 4 0.08 0.4 0.0141301130 0.3 fl.909a01 1 0 0

10、 0.3 0.1B2 2 0 0.55 U.133 3 0,5 0 0 Q013 1-19 00图3.2输入节点和支路数据对各个数字含义的解释如下:网络模型有四个节点，四条支路，编号见图3.1。第一个零下面三行数为支路参数，分别表示三条支路的起始和终止节点编号，后面的为电阻、电抗和 电纳，电导均为0,例如：1 2 0.1 0.4 0.01528第二个零下面的为变压器支 路，各数字意义同支路参数。接下去三行均为节点参数，分别表示注入有功功 率和无功功率。调用text.cpp文件，得到运行结果，见图3.3和图3.4。pSirk * C/- U5er&Ad minjt rdtorDesktopg片、

11、牛和拉夫逊法De b jgtext.exe计昇精 = le-B06导纳矩阵为.1-04209-8.24368j -0.S8823S*2-35294jB+3.6&703j-G.4S38Se*1.89107J-0.588235+2.35294J1.069-4.72738j 0*0J 0.48079+2.40385jB+3-66703J00j0-3.33333j Q+0J-0.453858+1.89107J-0-480769 *2-4B385j 0+0j0.934627-4.26159j第1迭代的节点电压幅值和角度呦度）0.?9345?0.?763111.11.05-0.00881675-0-107

12、S130.1151360第2迭代的节点电压幅值和角度恋度0,9848510.96S0441.1105-0.0070932-0.1104820.1190830第迭代的节点电压幅值和角度t弧度S.9846230.9647721.1丄.05-0.0867344-0-1125270.117552 Q第4迭代的节点电压幅值和角度那度）fl.9846190.96476?1.11.Q5-0.00871534-0-1125?0.11751图3.3运行结果1t C ：MJse rs Ad nn i n i stratoDe sttop返片 牛顿垃天辺法 ibLigltwxt 申 f3,98461?0.?&476

13、71,11.05-0.00871534-0Al25790-117510平衡节点功率为：9.3788G*jfl.264895线路上的功率为：S=9-2459S1+J-0.S=-B.5 +J-0-029001S=-0434&S63+j-0_136187S2i)=-0.2431&+j0-0110535& =B-J0S=0*J0S=-0.312?49+j-014036&-&.50.097016soa-e*j0S=0*ja531?=B*jaI& =0.0482113 +J0.10-164S(42=0-319671+j0_160255S43=0+j0ls=0+jaPress on# key to cont

14、inue图3.4运行结果23.3结果分析将上述仿真结果整理为表格 3.1、3.2，其中“ + ”表示节点i输出功率给节 点j，“-”表示节点j输出功率给i(纵向为i，横向为j)。表3.1节点有功功率输入与输出节点号一 123410+0.245981-0.5-0.0465632-0.2431600-0.31294930.500040.04821430.31967100表3.2节点无功功率输入与输出节点号123410-0.014708-0.029001-0.13618720.011050500-0.1403630.09701600040.104640.16025500根据表格计算:节点 1 有功功

15、率：0+0.245981-0.5-0.046563=-0.300582无功功率：0-0.014708-0.029001-0.136187=-0.179896节点2有功功率：-0.24316+0+0-0.312949=-0.556109无功功率：0.0110505+0+0-0.14036=-0.1293095节点3有功功率：0.5+0+0+0=0.5无功功率：0.097016+0+0+0=0.097016节点4有功功率：0.0482143+0.319671+0+0=0.3678853无功功率:0.10464+0.160255+0+0=0.264895根据已知条件，两个PQ节点的注入有功、无功分别

16、为:P仁0.3，Q1= 0.18 ；P2=0.55，Q2=0.130.3005820.30.3100%0.194%潮流计算误差：QQ0.180.179896100%0.0578%0.18巳 巳0.5561090.55100% 1.11%0.55QQ20.130.1293095100%0.531%0.13P可见，误差均在允许范围内线路损耗:S120.245981j 0.014708( 0.24316j 0.0110505)0.002821j 0.0036575S130.5 j0.029001(0.5 j 0.097016)j 0.068015S140.046563j 0.136187(0.048

17、2143j 0.10464)0.0016513j 0.031547S240.312949j 0.14036(0.319671j 0.160255)0.006722j 0.0198953.4结论通过上面的分析与计算，验证了程序的正确性。由于编写过程的不足，线路损耗没能直接计算出来，而是需要手算，比较遗憾。程序在运行过程中，需 要区分三种不同的节点，这由子程序保证实现。相比于快速分解法，牛拉法程 序较为复杂，但更精确一点，潮流误差较小。4总结本文基于牛顿-拉夫逊潮流算法的基本原理，利用 C+编程计算了一个4节 点简单电力网络的潮流，并验证了运行成果，误差在允许范围之内。因为牛拉 法计算过程中要不断

18、生成新的雅各布矩阵，所以相对来说占用内存较多，但收 敛速度快，这在程序运行过程中可以体现出来。本文程序并不是特别实用，因 为真正的电力网络不可能只有几个节点，而且各种电力设备的情况也会复杂很 多，因此程序会变得非常大，占用极大内存。但是，我还是通过这次练习，进 一步巩固了书本上的理论知识，了解了实际操作的过程步骤。最后，感谢杨伟 老师的悉心指导！参考文献1 朱红，赵琦，王庆宝 C+程序设计教程M.北京：清华大学出版社，2010.2 张伯明，陈寿孙，严正.高等电力网络分析M.北京：清华大学出版社，2007.3 吴明波.牛顿-拉夫逊法在潮流计算中的应用J.内蒙古科技与经济，2011,21：111-

19、112，115.4 明日科技.VisualC+从入门到精通M.北京：清华大学出版社，2012.5 顾洁，陈章潮，徐禧.一种新的配电网潮流算法-改进牛顿拉夫逊法J.华东电力，2000，5:10-12.6 朱文强.牛顿-拉夫逊法在配电网中的应用J.水利科技，2004，3:55-56,58.7 刘明波，谢敏，赵维新.大电网最优潮流计算M.北京：科学出版社，2010.附录/*/#includeiostream.h#includemath.h#include#include#include#define nodeNumber readParameter6using namespace std;/clas

20、s powerFlowCalculation / 导纳矩阵的计算类private:int readDataAmount; / 存放元件参数的读取个数 int balanceNodeindex;/ 平衡结点号double balanceNodeVoltage; / 平衡节点电压幅值 double balanceNodePowerP; / 平衡节点功率 double balanceNodePowerQ;double calculationAccuracy; / 计算精度 double readParameter7; / 宏定义 double *conductance;double *suscept

21、ance ;double *admittanceAmplitude;double *admittanceAngle;double *lineData;int PVNodeNumber; /PV 节点数double *jacobiMatrix; / 雅克比矩阵double *PQData; / 矩阵，存放 PQ 节点的数据double *PVData; / 矩阵，存放 PV 节点的数据double *voltageAmplitude; / 电压幅值double *voltageAngle; / 电压角度double *constantVector; / 常数向量double *lineConsu

22、meG;double *lineConsumeB;ifstream instream;public:/下标的转换int converIndex(int i,int j)int serial;serial=(i-1)*nodeNumber+j;return serial;/线路的计算void countLineBranch(double para,double*G,double*B)double GIJ,BIJ;int nii,njj,nij,nji;GIJ=para3/(para3*para3+para4*para4);BIJ=-para4/(para3*para3+para4*para4);

23、 nij=converIndex(para1,para2);nji=converIndex(para2,para1);nii=converIndex(para1,para1);njj=converIndex(para2,para2);Gnji=-GIJ;Gnij=-GIJ;Bnji=-BIJ;Bnij=-BIJ;Gnii+=GIJ;Gnjj+=GIJ;Bnii+= (BIJ+para5);Bnjj+= (BIJ+para5);/变压器的计算，规定 para1 为理想变压器侧， para2 为归算阻抗侧 void countTranformer(double para,double*G ,dou

24、ble*B)double GIJ,BIJ,b_i,b_j;int nii,njj,nij,nji;GIJ=para3/(para3*para3+para4*para4);BIJ=-para4/(para3*para3+para4*para4); b_i=BIJ*(1-para5)/para5/para5;b_j=BIJ*(para5-1)/para5;nij=converIndex(para1,para2);nji=converIndex(para2,para1); nii=converIndex(para1,para1);njj=converIndex(para2,para2);Gnij=-

25、GIJ/para5;Bnij=-BIJ/para5;Gnji=-GIJ/para5;Bnji=-BIJ/para5;Gnii+= GIJ/para5;Bnii+= (BIJ/para5+b_i);Gnjj+= GIJ/para5;Bnjj+= (BIJ/para5+b_j); /对地支路的计算 void countGroundBranch(double para,double*G ,double*B) int nii; nii=converIndex(para1,para1); Gnii+=para2;Bnii+=para3; /导纳矩阵的显示 void display()coutendl 导

26、纳矩阵为： endl;int nij;for(int i=1;i=nodeNumber;i+)for(int j=1;j=nodeNumber;j+)nij=converIndex(i,j); if(susceptance nij0) coutconductancenijsusceptance nijj ;elsecoutconductancenij+susceptance nijj ;coutendl;coutendl;/显示任意秩矩阵void displayMatrix(double *matr,int ranki,int rankj)int nij;for(int i=1;i=ranki

27、;i+)for(int j=1;j=rankj;j+)nij=(i-1)*rankj+j;coutmatrnij ; coutendl;/coutendl;/构造函数 , 打开文件 powerFlowCalculation(string input) instream.open(input(); assert(instream.is_open();/析构函数，关闭文件powerFlowCalculation()instream.close();/初始化，根据节点数开辟数组void initi() nodeNumber=readParameter0; / 读取节点数PVNodeNumber=0;

28、conductance=new doublenodeNumber*nodeNumber; susceptance =new doublenodeNumber*nodeNumber; admittanceAmplitude=new doublenodeNumber*nodeNumber; admittanceAngle=new doublenodeNumber*nodeNumber;PQData=new double(nodeNumber-1)*6; / 开辟 PQData PVData=new double(nodeNumber-1)*5; / 开辟 PVData voltageAmplitu

29、de=new doublenodeNumber; voltageAngle=new doublenodeNumber;jacobiMatrix=new doublenodeNumber*nodeNumber*4;constantVector=new doublenodeNumber*2; lineData=new doublenodeNumber*(nodeNumber-1)/2*5; lineConsumeG=new doublenodeNumber*nodeNumber; lineConsumeB=new doublenodeNumber*nodeNumber; for(int a=1;a

30、=nodeNumber*(nodeNumber-1)/2*5;a+) lineDataa=0.0;for(int i=1;i=nodeNumber;i+) / 电压向量的初始化voltageAmplitudei=1.0;voltageAnglei=0.0;for(int j=1;j=(nodeNumber*nodeNumber);j+) / 导纳初始化conductancej=0.0;susceptance j=0.0;admittanceAmplitudej=0.0;admittanceAnglej=0.0;lineConsumeGj=0.0;lineConsumeBj=0.0;for(in

31、t k=0;k=nodeNumber*nodeNumber*4;k+) / 雅克比矩阵初始化jacobiMatrixk=0.0;for(int n=1;n=nodeNumber*2;n+) / 常数向量初始化constantVectorn=0.0;/计算导纳矩阵主程序void countAdmittance()int lineDataSerial=1;for(int facility=0;facilityreadParameter0;if (readParameter00.5) break; / 序号为 0 时退出循环 for(int i=1;ireadParameteri;if(facili

32、ty=1)lineDatalineDataSerial=readParameteri; lineDataSerial+;if(facility=2)lineDatalineDataSerial=readParameteri;if(i=5) lineDatalineDataSerial=0.0; lineDataSerial-; lineDatalineDataSerial*=readParameteri; lineDataSerial-;lineDatalineDataSerial*=readParameteri; lineDataSerial+=2;/end for if(i=5)lineD

33、ataSerial+;/ 不同元件的计算/end for if(facility=2) /end for loop iswitch(facility) case 0:initi();balanceNodeindex=(int)readParameter1;balanceNodeVoltage=readParameter2; calculationAccuracy=readParameter3; cout 节点数 :nodeNumberendl 平衡节点 :readParameter1endl 平衡节点电压 :readParameter2endl 计算精度 :readParameter3endl

34、; break;case 1: countLineBranch(readParameter,conductance,susceptance ); break;case 2:countTranformer(readParameter,conductance,susceptance ); break;case 3:countGroundBranch(readParameter,conductance,susceptance ); break;/end switch/end for loop/end for (facility)/end countAdmittance() double vector

35、Angle(double a,double b)/*雅克比矩阵*/if(fabs(a)0.0000001 & fabs(b)0.0000001) return 0.0;double alpha;alpha=atan(fabs(b)/fabs(a);if(a=0)alpha=3.1415926-alpha;else if(a=0 & b=0 & b=0) alpha=-alpha;return alpha; /将直角坐标的导纳矩阵转换为极坐标形式 void converAdmittanceMatrix()int n;for(int i=1;i=nodeNumber;i+)for(int j=1;

36、j=nodeNumber;j+) n=converIndex(i,j);admittanceAmplituden=sqrt(conductancen*conductancen+susceptancen*susceptan cen);admittanceAnglen=vectorAngle(conductancen,susceptancen);/读取PQ节点和PV节点的数据void readNodePowerData()for(int j=1;jPQDataj;double dustbin; / 将 PQData 与 PVData 之间的 0 去掉instreamdustbin;for(int

37、i=1;iPVDatai; if(i-1)%5=0) if(PVDatai=0) break;PVNodeNumber+=1;/end for loop read PVData/电压幅值和角度赋初值void vectorAssignment() /PQ 节点的值已经在初始化中赋值了，这里给平衡节点和 PV 节点赋值int n;voltageAmplitudebalanceNodeindex=balanceNodeV oltage; voltageAnglebalanceNodeindex=0.0;for(int i=1;i=PVNodeNumber;i+)n=PVDataconverIndex

38、(i,2); voltageAmplituden=PVDataconverIndex(i,3); voltageAnglen=0.0;/生成 L 矩阵double creatL(int i,int j)double sum=0.0;double alpha=0.0;double beta=0.0;if(i=j) alpha=admittanceAngleconverIndex(i,i); for(int k=1;k=nodeNumber;k+) beta=-admittanceAngleconverIndex(i,k); sum+=admittanceAmplitudeconverIndex(

39、i,k)*voltageAmplitudek*sin(beta);/end for loop k sum*=voltageAmplitudei;sum=voltageAmplitudei*voltageAmplitudei*admittanceAmplitudeconverIndex(i,i)*sin(alpha );sum=-sum;elsebeta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*vol

40、tageAmplitudej*sin(beta); return sum;/生成N矩阵double creatN(i nt i,i nt j)double sum=0.0;double alpha=0.0;double beta=0.0;if(i=j)alpha=admittanceAngleconverIndex(i,i);for(int k=1;k=nodeNumber;k+)beta=-admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*cos(beta);/

41、end for loop k sum*=voltageAmplitudei;sum+=voltageAmplitudei*voltageAmplitudei*admittanceAmplitudeconverIndex(i,i) *cos(alpha);sum=-sum;elsebeta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=- voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*cos(beta);return

42、 sum;/生成 J 矩阵double creatJ(int i,int j)double sum=0.0;double beta;if(i=j)for(int k=1;k=nodeNumber;k+)if(k=i) continue; beta=voltageAnglei-voltageAnglek- admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*cos(beta); /end for loop ksum*=voltageAmplitudei; sum=-s

43、um;else beta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej* cos(beta);/end for ifreturn sum;/生成H矩阵double creatH(i nt i,i nt j)double sum=0.0; double beta;if(i=j) for(int k=1;k=nodeNumber;k+)if(k=i) continue; b

44、eta=voltageAnglei-voltageAnglek- admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*sin(beta);sum*=voltageAmplitudei;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=-1*voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplit

45、udej*sin(beta);/end for ifreturn sum;/生成雅各比矩阵void creatJacobiMatrix()int serial;int serialH;int serialN;int serialJ;int serialL;for(int i=1;i=2*(nodeNumber-1);i=i+2)for(int j=1;j=2*(nodeNumber-1);j=j+2)serialH=(i-1)*nodeNumber*2+j;serialN=(i-1)*nodeNumber*2+j+1;serialJ=(i-1)*nodeNumber*2+j+nodeNumbe

46、r*2;serialL=(i-1)*nodeNumber*2+j+nodeNumber*2+1;jacobiMatrixserialH=creatH(i+1)/2,(j+1)/2);jacobiMatrixserialN=creatN(i+1)/2,(j+1)/2);jacobiMatrixserialJ=creatJ(i+1)/2,(j+1)/2);jacobiMatrixserialL=creatL(i+1)/2,(j+1)/2);/end for loop j/end for loop ifor(int n=(nodeNumber-PVNodeNumber)*2;n=nodeNumber

47、*2;n=n+2)for(int m=1;m=nodeNumber*2;m+)if(n=m) serial=(m-1)*nodeNumber*2+n; jacobiMatrixserial=1.0;elseserial=(m-1)*nodeNumber*2+n; jacobiMatrixserial=0.0; serial=(n-1)*nodeNumber*2+m; jacobiMatrixserial=0.0;/end for if/end for loop m/end for loop n serial=nodeNumber*nodeNumber*4-nodeNumber*2-1; jacobiMatrixserial=1.0;/ rank 为 6 的下标转换int converSerial(int i,int j,int rank=6)int serial;serial=(i-1)*rank+j;return serial;/生成常数向量void creatConstantVector()int i; /i 为 P、 Q 的下标double sum=0.0;double beta=0.0;for(int n=1;n=(nodeNumber-1)*2;n+) /n 为向量的下标i=(n+1)/2; /i 为节点号sum=0.0