电力系统分析课程设计潮流计算

《电力系统分析课程设计潮流计算》由会员分享，可在线阅读，更多相关《电力系统分析课程设计潮流计算（33页珍藏版）》请在装配图网上搜索。

1、目录摘要11.任务及题目要求22.计算原理32.1牛顿拉夫逊法简介32.2牛顿拉夫逊法的几何意义73计算步骤74.结果分析9小结11参考文献12附录：源程序13本科生课程设计成绩评定表33摘要电力系统的出现，使高效，无污染，使用方便，易于调控的电能得到广泛应用，推动了社会生产各个领域的变化，开创了电力时代，发生率第二次技术革命。电力系统的规模和技术水准已经成为一个国家经济发展水平的标志之一。电力系统稳态分析包括潮流计算和静态安全分析。电力系统潮流计算是电力系统最基本的计算，也是最重要的计算。所谓潮流计算，就是已知电网的接线方式与参数及运行条件，计算电力系统稳态运行各母线电压、个支路电流与功率及

2、网损。对于正在运行的电力系统，通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提供原始数据。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在

3、电力系统运行状态的实时监控中，则采用在线潮流计算。关键词：电力系统 潮流计算 牛顿-拉夫逊法1.任务及题目要求对如下所诉系统编程进行潮流计算：节点数：3 支路数：3 计算精度：0.00010支路1： 0.0300+j0.0900 12支路2： 0.0200+j0.0900 23支路3： 0.0300+j0.0900 31节点1：PQ节点，S（1）=-0.5000-j0.2000节点2：PQ节点，S（2）=-0.6000-j0.2500节点3：平衡节点，U（3）=1.00000.0000经分析可知题目所给的系统为三节点组成的环形回路，且均为线路没有变压器。有两个为PQ节点，这类节点的有功功率P和

4、无功功率Q是给定的，节点电压是待求量。节点3为平衡节点，在潮流分布算出来以前，网络中的功率损失是未知的，因此，至少要有一个节点的有功功率不能给定，这个节点承担了系统的有功功率平衡，故称之为平衡节点。对系统进行潮流计算，要求各线路的功率分布及功率损耗，未知节点的电压等，并进行编程运行得到结果。2.计算原理潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如各母线上的电压（幅值及相角），网络中的功率分布及功率损耗等。常用的计算潮流分析的方法有牛顿拉夫逊法，PQ分解法等。本次设计采用牛顿拉夫逊法进行计算。牛顿拉夫逊法(简称牛顿法)在数学上是求解非线性代数方程式

5、的有效方法，其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求解的过程，即通常所称的逐次线性化过程。2.1牛顿拉夫逊法简介牛顿拉夫逊法（NewtonRaphson法）是求解非线性方程代数方程组的有效迭代计算方法。在牛顿拉夫逊法的每一次迭代过程中，对非线性方程通过线性化处理逐步近似。下面以单变量加以说明。设有单变量非线性方程 (2-1)求解此方程时。先给出解的近似值它与真解的误差为，则将满足方程，即 (2-2)将(2-2)式左边的函数在附近展成泰勒级数，于是便得 （2-3）式中,分别为函数在处的一阶导数，.，n阶导数。如果差值很小，2-3式右端的二次及以上阶次的各项均可略去。于是

6、，2-3便简化为 0 (2-4)这是对于变量的修正量的现行方程式，亦称修正方程式。解此方程可得修正量 (2-5)用所求的去修正近似解，变得 (2-6)由于2-6是略去高次项的简化式，因此所解出的修正量也只是近似值。修正后的近似解同真解仍然有误差。但是，这样的迭代计算可以反复进行下去，迭代计算的通式是 (2-7)迭代过程的收敛判据为 (2-8)或 (2-9)式中，为预先给定的小正数。 这种解法的几何意义可以从图1得到说明。函数yf(x)为图中的曲线。f(x)0的解相当于曲线与x轴的交点。如果第k次迭代中得到，则过点作一切线，此切线同x轴的交点便确定了下一个近似值。由此可见，牛顿拉夫逊法实质上就是

7、切线法，是一种逐步线性化的方法。 应用牛顿法求解多变量非线性方程组2-1时，假定已给出各变量的初值，. ，令，. 分别为各变量的修正量，使其满足方程 (2-10) 将上式中的n个多元函数在初始值附近分别展成泰勒级数,并略去含有，二次及以上阶次的各项，便得. (2-11)方程式2-11也可以写成矩阵形式 (2-12) 方程式2-11是对于修正量， 的线性方程组,称为牛顿法的修正方程式.利用高斯消去法或三角分解法可以解出修正量，。然后对初始近似值进行修正 (i=1,2,.,n) (2-13)如此反复迭代，在进行k1次迭代时，从求解修正方程式 (2-14)得到修正量，并对各变量进行修正 (i=1,2

8、,n) (2-15)也可以缩写为 (2-16)和 (2-17) 式中的X和分别是由n个变量和修正量组成的n维列向量；F(X)是由n个多元函数组成的n维列项量；J是n阶方阵，称为雅可比矩阵，它的第i、j个元素是第n个函数对第j个变量的偏导数；上角标(k)表示阵的每一个元素都在点处取值。迭代过程一直到满足收敛判据 (2-18) 或 (2-19)为止。和为预先给定的小正数。 将牛顿拉夫逊法用于潮流计算，要求将潮流方程写成形如方程式2-1的形式。由于节点电压可以采用不同的坐标系表示，牛顿拉夫逊法潮流计算也将相应的采用不同的计算公式。2.2牛顿拉夫逊法的几何意义牛顿-拉夫逊法几何意义：图中为假设的初始值

9、它和真值之间的差值为：，然后求得：，其与真值之间的差值为：。如此类推，求得逼近真值。图1牛顿法的几何解释3计算步骤用牛顿拉夫逊法计算电力系统潮流的基本步骤：求解过程大致可以分为以下步骤：(1)形成节点导纳矩阵；(2)将各节点电压设初值 ；(3)将节点初值代入相关求式，求出修正方程式的常数项向量；(4)将节点电压初值代入求式，求出雅可比矩阵元素；(5)求解修正方程，求修正向量；(6)求取节点电压的新值；(7)检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步；(8)计算支路功率分布，节点无功功率和平衡节点功率。图2中示出牛顿拉夫逊法计算潮流的程序框

10、图。启动输入原始数据形成节点导纳阵给定电压初值e（0）、f（0）k=0根据公式计算用公式计算雅克比矩阵各元素高斯法解修正方程，求修正节点电压是计算平衡节点功率及全部线路功率输出k=k+1否图2计算流程图4.结果分析经过用C语言编程运算得到第一次迭代结果为： 节点导纳矩阵为: 6.66667+j-20.00000-3.33333+j10.00000-3.33333+j10.00000-3.33333+j10.00000 5.68627+j-20.58824-2.35294+j10.58824-3.33333+j10.00000-2.35294+j10.58824 5.68627+j-20.588

11、24 Gauss_Seidel迭代结果：第1次迭代PQ节点e及f的值：0.983500 -0.019500第1次迭代PQ节点e及f的值：0.973610 -0.033997第2次迭代PQ节点e及f的值：0.969642 -0.036485第2次迭代PQ节点e及f的值：0.965875 -0.042673 Jacobi 矩阵:19.31452 6.80237-9.57480-3.59699-7.5856018.98468 3.59699-9.57480-9.51651-3.6463119.92097 5.76424 3.64631-9.51651-6.9773519.36506 不平衡值：-0.

12、114256-0.054375-0.002280-0.007412输出 df,de 节点为 1 的电压不平衡值： df=-0.00508 de=-0.00267 节点为 2 的电压不平衡值： df=-0.00006 de=-0.00037-输出迭代过程中的电压值 节点为 1 的电压修正新值： f=-0.04156 e= 0.96698 节点为 2 的电压修正新值： f=-0.04273 e= 0.96551 还不能达到题目中的精度要求，继续进行迭代，经6次迭代后得到结果。 各节点电压为: U1=0.963186 + j-0.042323 U2=0.962333 + j-0.045732 U3=

13、1.000000 + j0.000000 平衡节点功率为: S3=1.118788+j0.518288 输电线路功率如下: * 输电线路1-2的功率: 0.035453 + j-0.004289 * 输电线路2-1的功率: -0.035412 + j0.004412 * 输电线路2-3的功率: -0.564587 + j-0.254058 * 输电线路3-2的功率: 0.572847 + j0.291225 * 输电线路3-1的功率: 0.545942 + j0.227063 * 输电线路1-3的功率: -0.535453 + j-0.195598 输电线路上损耗的功率为: * 输电线路1-2

14、上的功率损耗: 0.000041 + j 0.000123 * 输电线路2-3上的功率损耗: 0.008259 + j 0.037167 * 输电线路3-1上的功率损耗: 0.010488 + j 0.031465 网络总损耗: 0.018789 + j 0.068755结果满足了题目中所需的精度要求，并且计算出了潮流计算要求的各节点电压，网络中的功率分布及功率损耗。小结在这次课程设计中，我强烈感觉到自己在很多方面的不足，对别人的依赖性比较强。我想我会在以后的学习中不断去发现自己的不足，并一一改正，希望在以后的学习中不要犯同样的错误。在这次课程设计中，我们尽量按照老师的要求做，但在具体的操作过

15、程中，还是出现了很多的问题。搞完这个课程设计让我感觉电力系统分析是一门很有用的课程。因为我对它的学到的知识比较少。在很多时候我很多东西都不了解。并且走了很多的弯路。而且我感觉自己的知识不够连贯。好些时候都出现了卡壳的情况。这次课程设计后，我一定要重新对电力系统分析这门课程做进一步的了解。对在此过程中遗留下的问题做好好的研究。争取早点对电力系统分析这门课程有个全方位的了解。为在以后的毕业课程设计中多些方案。也为我子以后走上工作岗位，提升自己的专业技能，打下扎实的基础。还有就是在十几天的课程设计中，使我养成了很好的学习习惯，和对学习知识的严谨的态度，同时也养成了积极查阅相关资料的好习惯，好习惯的养

16、成是来之不易的，我相信在以后的学习和工作中，我将继续保持这些良好的习惯，并积极努力的学习。让自己更上一层楼。也感谢班上同学给了我很多宝贵的意见和参考，使我获益很多。参考文献：1 何仰赞 温增银 电力系统分析 华中科技大学出版社 20062 何仰赞 温增银 电力系统分析题解 华中科技大学出版社 20063 王锡凡 现代电力系统分析 科学出版社 20044 于永源 杨绮雯 电力系统分析（第二版）中国电力出版社 20075 薛定宇 控制系统计算机辅助设计 清华大学出版社 1996附录：源程序#include #include #define M 40 /*矩阵阶数*/#define N 10 /*最

17、大迭代次数*/ int i,j,k,a,b,l,z; /* 循环变量 */ int n, /* 节点数 */ m, /* 支路数 */ dd, /*对地支路数*/ pq, /* PQ节点数 */ pv, /* PV节点数 */ byq; /* 变压器数 */ float eps, /* 精度 */ bb, /* 变比 */ max; /*指标标志*/ float yyM; /*中间变量*/ static float GMM,BMM,B1MM,B2MM,B3MM; struct jiedian /* 定义节点结构体 */ int num,s; /* num为节点号，s为节点类型*/ float

18、p,q,u,f,v; float dp,dq,df,du; jiedianM; struct zhilu /* 定义支路结构体 */ int num; int p1,p2; /*支路的两个节点*/ float r,x; /*支路的电阻与电抗*/ zhiluM; FILE *fp1,*fp2; /*文件指针*/ void data1() /* 读取数据函数 */ int h,numb,a,b; a=1; fp1=fopen(d:input.txt,r); if(fp1=NULL) printf( can not open file !n); exit(0); fscanf(fp1,%d,%d,%

19、d,%d,%d,%d,%fn,&n,&m,&dd,&pq,&pv,&byq,&eps); /*输入节点数,支路数,对地支路数,PQ节点数，PV节点数，变压器数和精度*/ j=1;k=pq+1; for(i=1;i=n;i+) /*输入节点类型的输入功率和节电电压初值*/ fscanf(fp1,%d,%d,&numb,&h);if(h=1) /*类型h=1是PQ节点*/ fscanf(fp1,%f,%f,%f,%fn,&jiedianj.p,&jiedianj.q,&jiedianj.u,&jiedianj.f); /*输入PQ节点的数据*/ jiedianj.num=numb; jiedian

20、j.s=h; j+; if(h=2) /*类型h=2是pv节点*/ fscanf(fp1,%f,%f,%f,%fn,&jiediank.p,&jiediank.v,&jiediank.u,&jiediank.f); /*输入PV节点的数据*/ jiediank.num=numb; jiediank.s=h; k+; if(h=3) /*类型h=3是平衡节点*/ fscanf(fp1,%f,%fn,&jiediann.u,&jiediann.f); /*输入平衡节点的数据*/ jiediann.num=numb; jiediann.s=h; for(i=1;i=m;i+) /*输入支路阻抗*/

21、fscanf(fp1,%d,%d,&zhilui.num,&h);if(h=0) /*输入非变压器支路阻抗*/fscanf(fp1,%d,%d,%f,%fn, &zhilui.p1,&zhilui.p2,&zhilui.r,&zhilui.x); if(h=1) /*输入变压器支路阻抗*/ fscanf(fp1,%f,%d,%d,%f,%fn,&bb,&zhilui.p1,&zhilui.p2,&zhilui.r, &zhilui.x); zhilui.r=zhilui.r*bb; zhilui.x=zhilui.x*bb; b=m+a; zhilub.num=b; zhilub.p1=zhi

22、lui.p1; zhilub.p2=0; zhilub.r=zhilui.r*bb/(1-bb); zhilub.x=zhilui.x*bb/(1-bb); zhilub+1.num=b zhilub+1.p1=zhilui.p2; zhilub+1.p2=0; zhilub+1.r=zhilui.r/(bb-1); zhilub+1.x=zhilui.x/(bb-1); a=a+2; fclose(fp1); if(fp2=fopen(d:output.txt,w)=NULL) printf( can not open file!n); exit(0); fprintf(fp2,nn * 原

23、始数据 *n); fprintf(fp2, =n) fprintf(fp2, 节点数:%2d 支路数:%2d 对地支路数:%2d 变压器数:%2d n PQ节点数:%2d PV节点数:%2d 精度:%fn,n,m+2*byq,dd,byq,pq,pv,eps); fprintf(fp2, -n); for(i=1;i=pq;i+) fprintf(fp2, PQ节点 节点%2d P%d=%+f Q%d=%+fn, jiediani.num,jiediani.num,jiediani.p,jiediani.num,jiediani.q); for(i=pq+1;i=pq+pv;i+) fprin

24、tf(fp2, PV节点 节点%2d P%d=%+f V%d=%+fn, jiediani.num,jiediani.num,jiediani.p,jiediani.num,jiediani.v); fprintf(fp2, 平衡节点 节点%2d u%d=%+f f%d=%+fn, jiediann.num,jiediann.num,jiediann.u,jiediann.num,jiediann.f); fprintf(fp2, -n); for(i=1;i=m+2*byq;i+) fprintf(fp2, 支路%2d 相关节点:%2d,%2d R=%+f X=%+fn, i,zhilui.

25、p1,zhilui.p2,zhilui.r,zhilui.x); fprintf(fp2, =n); void form_y() /* 形成节点导纳矩阵函数 */ float S; for(i=1;i=m+2*byq;i+) /*节点导纳矩阵的主对角线上的导纳*/ for(j=0;j=n;j+) if(zhilui.p1=j)|(zhilui.p2=j) S=zhilui.r*zhilui.r+zhilui.x*zhilui.x; if(S=0) continue; Gjj+=zhilui.r/S; Bjj+=-zhilui.x/S; for(i=1;i=m+2*byq;i+) /*对地导纳*

26、/ j=zhilui.p1; k=zhilui.p2; if(j=0)|(k=0) S=zhilui.r*zhilui.r+zhilui.x*zhilui.x; if(S=0) G0i=B0i=Gi0=Bi0=0;continue; if(j=0) G0k+=zhilui.r/S; B0k+=-zhilui.x/S; Gk0=G0k; Bk0=B0k; if(k=0) G0j+=zhilui.r/S; B0j+=-zhilui.x/S; Gj0=G0j; Bj0=B0j; for(k=1;k=m-dd;k+) /*节点导纳矩阵非主对角线上的导纳*/ i=zhiluk.p1; j=zhiluk.

27、p2;S=zhiluk.r*zhiluk.r+zhiluk.x*zhiluk.x;if(S=0) continue;Gij+=-zhiluk.r/S;Bij+=zhiluk.x/S;Gji=Gij;Bji=Bij; for(i=1;i=(pq+pv);i+) /*形成B矩阵*/ for(j=1;j=(pq+pv);j+) B1ij=Bjiediani.numjiedianj.num; for(i=1;i=pq;i+) /*形成B矩阵*/ for(j=1;j=pq;j+) B2ij=Bjiediani.numjiedianj.num; /*输出节点导纳矩阵、系数矩阵B、系数矩阵B*/ fprin

28、tf(fp2,nn * 计算结果 *n); fprintf(fp2,=n) fprintf(fp2,n 节点导纳矩阵为:); for(i=1;i=n;i+) fprintf(fp2,n ); for(j=1;j=n;j+) fprintf(fp2,%+8.6f+j%+8.6f ,Gij,Bij); fprintf(fp2,n -n); fprintf(fp2,n 节点导纳矩阵实部G为:); for(i=1;i=n;i+) fprintf(fp2,n ); for(j=1;j=n;j+) fprintf(fp2,%+8.6f ,Gij); fprintf(fp2,n -n); fprintf(f

29、p2,n 节点导纳矩阵虚部B为:); for(i=1;i=n;i+) fprintf(fp2,n ); for(j=1;j=n;j+) fprintf(fp2,%+8.6f ,Bij); fprintf(fp2,n -n); fprintf(fp2,n 系数矩阵B为:); for(i=1;i=pv+pq;i+) fprintf(fp2,n ); for(j=1;j=pv+pq;j+) fprintf(fp2,%+8.6f ,B1ij); fprintf(fp2,n -n); fprintf(fp2,n 系数矩阵B为:); for(i=1;i=pq;i+) fprintf(fp2,n ); fo

30、r(j=1;j=pq;j+) fprintf(fp2,%+8.6f ,B2ij); fprintf(fp2,n =n); void gaoshi(int A,float aaMM,float bbM) /*高斯计算函数*/int i,j,z; int k=0; float a; for(i=k+1;iA;i+) for(z=0;(z+i)A;z+) aai+zk=aai+zk/aakk; /*校验过程*/ bbi+z=bbi+z-aai+zk*bbk; for(j=k+1;j=0;i-) a=0; for(j=i+1;jA;j+) a=aaij*bbj+a; bbi=(bbi-a)/aaii;

31、 float solve(float max) /*迭代函数*/ float aaM,bbM;/*求节点有功功率不平衡量dp*/ for(i=1;i=(pq+pv);i+) /*求dp*/ jiediani.dp=jiediani.p; for(j=1;j=n;j+) jiediani.dp=jiediani.dp-jiediani.u*jiedianj.u*(Gjiediani.numjiedianj.num*cos(jiediani.f-jiedianj.f)+Bjiediani.numjiedianj.num*sin(jiediani.f-jiedianj.f); for(i=1;i=(

32、pq+pv);i+) /*定义各节点的dP/U*/ if(jiediani.s!=3) yyi=fabs(jiediani.dp/jiediani.u); max=yy1; /*求dP/U的最大值*/ for(i=1;imax)max=yyi; for(i=1;i=(pq+pv);i+) fprintf(fp2, 有功功率不平衡量dP%d=%+fn,jiediani.num,jiediani.dp); fprintf(fp2,n -n); /*求各节点的相位角f，df*/ a=1; /*装换为高斯函数的标准形式*/ for(i=1;i=n;i+) if(jiediani.s=3)a=a+1;c

33、ontinue; aai-a=jiediani.dp/jiediani.u; b=1; for(j=1;j=n;j+) if(jiedianj.s=3)b=b+1;continue; B3i-aj-b=0.0-B1i+1-aj+1-b*jiediani.u; gaoshi(pv+pq),B3,aa); /*代入高斯函数求解方程*/ a=1; for(i=1;i=n;i+) if(jiediani.s=3)a+;continue; jiediani.df=aai-a; jiediani.f=jiediani.f+jiediani.df; for(i=1;i=(pq+pv);i+) fprintf

34、(fp2, df%d=%f 电压的相位角f%d=%+fn,jiediani.num,jiediani.df,jiediani.num,jiediani.f); fprintf(fp2,n -n); /*求节点的无功功率不平衡量dq*/ for(i=1;i=n;i+) /*求dq*/ jiediani.dq=jiediani.q; for(j=1;j=n;j+) jiediani.dq=jiediani.dq-jiediani.u*jiedianj.u*(Gjiediani.numjiedianj.num*sin(jiediani.f-jiedianj.f)-Bjiediani.numjiedi

35、anj.num*cos(jiediani.f-jiedianj.f); for(i=1;i=pq;i+) /*定义各节点的dQ/U*/ yyi=fabs(jiediani.dq/jiediani.u); for(i=1;imax)max=yyi; for(i=1;i=pq;i+) fprintf(fp2, 无功功率不平衡量dQ%d=%+fn,jiediani.num,jiediani.dq); fprintf(fp2,n -n); /*计算各节点电压*/ for(i=1;i=n;i+) /*装换为高斯函数的标准形式*/ for(j=1;j=n;j+) B3i-1j-1=0-B2ij; a=1;

36、 for(i=1;i=n;i+) if(jiediani.s!=1)a=a+1;continue; bbi-a=jiediani.dq/jiediani.u; gaoshi(pq,B3,bb); /*代入高斯函数求解方程*/ a=1; for(i=1;i=n;i+) if(jiediani.s=3)a+;continue; jiediani.du=bbi-a; jiediani.u=jiediani.u+jiediani.du; for(i=1;i=n;i+) if(jiediani.s=1) fprintf(fp2, dU%d=%f 电压的大小U%d=%+fn,jiediani.num,ji

37、ediani.du,jiediani.num,jiediani.u); fprintf(fp2,n =n); return(max); /*返回循环标志值*/float mul_Re(x1,y1,x2,y2) /*求实部*/float x1,x2,y1,y2;float x; x=x1*x2-y1*y2; return(x); float mul_Im(x1,y1,x2,y2) /*求虚部*/float x1,x2,y1,y2; float y; y=x1*y2+x2*y1; return(y); void data2() /* 潮流计算结果 */ float sp=0,sq=0,mo1,mo

38、2,mo3,mo4,mo5,mo6,ff,ttM; int i1,j1,i11=0,j11=0; static float PMM,QMM; fprintf(fp2,nn 各节点电压为:); for(i=1;i=n;i+) ttM=jiediani.u;jiediani.u=(jiediani.u)*(cos(jiediani.f); /*由极坐标转化为直角坐标*/ jiediani.f=(ttM)*(sin(jiediani.f); fprintf(fp2,n U%d=%+8.6f + j %+8.6f, jiediani.num,jiediani.u,jiediani.f); /*输出直角

39、坐标节点电压*/ fprintf(fp2,n =n); fprintf(fp2,nn 平衡节点功率为:); for(i=1;i=n;i+) i1=jiediann.num; j1=jiediani.num; sp+=mul_Re(Gi1j1,-Bi1j1,jiediani.u,-jiediani.f); sq+=mul_Im(Gi1j1,-Bi1j1,jiediani.u,-jiediani.f); jiediann.p=mul_Re(jiediann.u,jiediann.f,sp,sq); jiediann.q=mul_Im(jiediann.u,jiediann.f,sp,sq); fp

40、rintf(fp2,n S%d=%+8.6f+j %+8.6f,jiediann.num,jiediann.p,jiediann.q); fprintf(fp2,n =n); fprintf(fp2,nn 线路功率如下:n); for(k=1;k=(m-dd);k+) i1=zhiluk.p1; j1=zhiluk.p2; for(l=1;l=n;l+) if(jiedianl.num=i1)i=l; if(jiedianl.num=j1)j=l; sp=mul_Re(jiediani.u-jiedianj.u),(-jiediani.f+jiedianj.f),-Gi1j1,Bi1j1);

41、sq=mul_Im(jiediani.u-jiedianj.u),(-jiediani.f+jiedianj.f),-Gi1j1,Bi1j1); mo1=jiediani1.u*jiediani1.u+jiediani1.f*jiediani1.f; Pi1j1=mul_Re(jiediani.u,jiediani.f,sp,sq)+mo1*Gi10; Qi1j1=mul_Im(jiediani.u,jiediani.f,sp,sq)-mo1*Bj10; sp=mul_Re(jiedianj.u-jiediani.u),(-jiedianj.f+jiediani.f),-Gj1i1,Bj1i1); sq=mul_Im(jiedianj.u-jiediani.u),(-jiedianj.f+jiediani.f),-Gj1i1,Bj1i1); mo2=jiedianj1.u*jiedianj1.u+jiedianj1.f*jiedianj1.f; Pj1i1=mul_Re(jiedianj.u,jiedianj.f,sp,sq)+mo2*Gj10; Qj1i1=mul_Im(jiedianj.u,jiedian