电气项目工程及其自动化专业综合实训一{新}

1、电气工程及其自动化专业综合实训一原始资料网络接线如图，各支路阻抗和各节点功率均以标幺值标于图中,其中节点2连接的实际是发额定功率的发电厂。设节点1的电压保持为1.06，试计算图示系统中的潮流分布。（计算精度10-5）二、等值电路图s；4 = nL4-i- jn.m4U-0 厲 5+*|ai5AI 3IICl.umrQDbQlOC+JOOO.og-i-jD.50tJg4jLj .151/=+JO,1电力系统潮流计算二、报告要求:1、写出潮流计算的原理 2、写出功率方程（极坐标或直角坐标）3、写出雅克比矩阵元素的表达式或 PQ分解法中的系数矩阵4、利用MATLABL语言或其他语言编程调试 5、心得

2、或体会 6、附程序代码，打印1、写出潮流计算的原理潮流是指在发电机母线上功率被注入网络，而在变（配）电站的 母线上接入负荷，其间，功率在网络中流动。对于这种流动的功率, 电力生产部门称之为潮流。以电力网络潮流、电压计算为主要内容的 电力网络稳态行为特性计算的目的在于估计对用户电力供应的质量以 及为电力网运行的安全性与经济性评估提供基础数据。配电网潮流计 算是配电网络分析的基础，配电网的网络重构、无功功率优化、状态 估计和故障处理都需要用到配电网潮流数据。电力系统稳态运行应满足以下要求：1）满足系统经济性运行的要求，每一台发电机的输出必须接近于 预先设定值;2）必须确保联络线潮流低于线路热极限和

3、电力系统稳定极限;3）必须保持某些中枢点母线上的电压水平在容许范围内，必要时 用无功功率补偿计划来达到;4）区域电网是互联系统的一部分，必须执行合同规定的输送至邻 网的联络线功率计划;5）用故障前的潮流控制策略使事故扰动效应最小化。通常情况下，输电线路电压在轻载时会较高，重载时会较低，电 压调整是指在负载由轻载到满载变化过程中实时调整线路电压满足运行要求；对于超高压输电线路，线路电压维持在额定电压的士5浓内,实际运行时，通常电压调整约为10%。对于低压输电线路，电压调整 数值为10%包含了变压器本身的电压降落。3.1.1潮流计算的基本物理量潮流计算是电力系统分析中的一种最基本的计算，它的任务是

4、对 给定的运行条件确定系统的运行状态，就是在三相平衡稳态状态下计 算电力系统中每条母线的电压幅值和相角，其中每一设备如传输线和 变压器中的有功和无功潮流，以及各设备的损耗都需要计算出来。潮流计算采用电力系统的单线图，对于任意一条母线i ,需要以下四个变量描述：电压幅值Ui、相角 电网供给母线的有功Pi、无功Qi。若某一电力系统有N个节点，则共 有4N个变量，对于每条母线，这些变量中的两个指定为输入数据，其 它的两个是潮流程序所要计算的未知量。为方便起见，在图 3.1中传 送给母线i的功率可分为发电机发出和负载吸收两部分。也就是Pi = PGi - PLiQi = QGi - QLi图3-1节点

5、的变量每条母线被归分为以下三种母线类型中的某一类: 1）平衡节点，一般一个系统只有一个平衡节点。在潮流分布算出以前，网络中的功率损耗是未知的，因此，至少有一个节点的有功功率P和无功功率Q不能给定。另外必须选定一个节点，制定其电压相角为零，作为其它节点电压相位的参考，这个节点叫基准节点。为了计算方便，常将平衡节点和基准节点设在同一个节点上。为方便起见在本书中把它标号为母线 1。平衡节点是电压参考节点，该母线的是给定值，作为输入。潮流程数据，典型取标幺值 序计算P1和Q1。因为平衡节点的P、Q事先无法确定，为使潮流计 算结果符合实际，常把平衡节点选在有较大调节裕量的发电机节 点，潮流计算结束时若平

6、衡节点的有功功率、无功功率和实际情况 不符，就要调整其他节点的边界条件以使平衡节点的功率满足实际 允许范围。2) PQ节点，Pi和Qi是输入数据。这类节点的有功功率 Pi和 无功功率Qi是给定的，潮流计算程序计算节点电压幅值 Ui和相角。负荷节点和无功功率注入的联络节点都属于这类节点。有些情况下，系统中某些发电厂送出的功率在一定时间内为固定时，该发电厂母线也可以作为PQ节点。在一个典型的潮流程序中绝大多数母线作为 PQ 节点。3）PU节点（电压控制母线），Pi和Ui是输入数据。这类节点的有功功率Pi和节点电压幅值Ui是给定的，潮流程序计算节点的。这类无功功率Qi和电压相角 节点必须具有足够的无

7、功可调容量，用以保持给定的节点电压幅 值。在电力系统中这类节点的数目较少。例如与发电机、并联补偿 电容器或者静止无功系统相连的母线。设备无功功率最大值QGimax 和最小值QGimin都是输入数据。另一个例子是与抽头可调节变压 器相连的母线，用潮流程序计算抽头的位置。注意，当母线i是无发电机相连接的负载母线时，Pi=-PLi为 负值；也就是说，在图3.1中给母线i提供的有功为负值。如果负 荷是感性的，Qi=-QLi为负值。综上所述，若系统中有n个节点，n为网络中除去参考节点外 的节点数，本书中以大地为参考节点，选第 1个节点为平衡节点, 剩下的n-1个节点中有r个PU节点，则有n-r-1个PQ

8、节点。因此, 除了平衡节点外，有n-1个节点的注入有功功率、n-r-1个PQ节点 的注入无功功率和r个PU节点电压幅值为已知量。3.1.2潮流计算的数学模型在稳态潮流计算中，电力系统各元件（参数）等效成一个有源 网络。将发电机和负荷用无阻抗线从网络中抽出，剩下的是由接地和不接地支路组成的无源线性网络，可以用导纳矩阵（Y）或阻抗 矩阵（Z）来描述。采用导纳矩阵时，节点电流和节点电压构成以下方程:(3-1 )其中：丫为nXn阶导纳矩阵，其阶数n为网络中除去参考节点 外的节点数，如果不考虑网络元件的非线性及变压器的相位偏移,Y为对称矩阵,为nxi维节点注入电流列向量，在电力系统计算中，节点注入电流可

9、理解为该节 点电源电流与负荷电流之和，并规定流入节点电流为正。因此仅有负荷的节电电流就为负值，某些仅起联络作用的节点，图3-2中节点n=3,其注入电流为零。nxi维节点电压列向量。网络中有接地支路时，节点电压通常指该 节点的对地电压，以大地作为参考节点，并规定其编号为零。图3-2运用节点电压法时的电力网路等值电路对于第i个节点，展开为如下形式:(3-3)若采用阻抗矩阵可表示为展开为:在潮流计算时一般以节点电压方程进行。节点导纳矩阵与阻抗 矩阵互为逆阵，在短路计算时可直接利用导纳矩阵求逆得到阻抗矩 阵以求得短路点的短路电流。由于实际系统中一般不给出节点电流而是节点功率，因此将式(3-2 )中的节

10、点注入电流用节点注入功率来表示为:(3-4)如果节点电压用极坐标表示，令n个节点电力系统的潮流方程的一般形式是(3-5)(3-6)若米用直角坐标系，节点电压可表示为导纳矩阵元素可以表示为将上述表达式带入式（2-8）的右端，展开并分出实部和虚部,便得(3-7)可见，原来电流电压的线性方程组变换为功率和电压的非线性 方程组，式（3-6）（ 3-7 ）就是潮流计算的基本方程。它是一组共 有n个非线性方程组成的复数方程组，如果把实部和虚部分开便得 到2n个实数方程，由该方程组可解出2n个运行参数。但是每一个 节点都有P、Q U 四个运行变量，共有4n个运行参数，所以要 事先给定其余2n个参数。这就要根

11、据节点的分类，将每个节点的4 个运行参数中的两个作为原始数据，另外两个作为待求量。3.1.3潮流计算的约束条件为了保证电力系统的正常运行，潮流问题中某些变量应满足一 定的约束条件，常用的约束条件有:（1）所有节点电压必须满足从保证电能质量和供电安全的要求来看， 电力系统的所有电气 设备都必须运行在额定电压附近。PU节点的电压幅值必须按上述条 件给定。因此，这一定约束条件主要是对 PQ节点而言。(2) 所有电源节点的有功功率和无功功率必须满足PQ节点的有功功率和无功功率以及 PU节点有功功率，在给定 时就必须满足此条件。因此，对平衡节点的P和Q以及PU节点的Q 应按此条件进行检验。(3) 某些节

12、点之间电压的相位差应满足为了保证系统运行的稳定性，要求某些输电线两端的电压相位差不 超过一定的数值。因此，潮流计算可以归结为求解一组非线性方程 组，并使其解答满足一定的约束条件。如果不能满足，则应修改某 些变量。甚至修改系统的运行方式，重新进行计算。2、写出功率方程（极坐标或直角坐标）3、写出雅克比矩阵元素的表达式或PQ分解法中的系数矩解：（1）、将阻抗转换为导纳矩阵:=0.4 5-1-JO. 15s4=O,4-FjO,Ob5.19-J2.77s2=0,254-j0.25I心-J3.546.70-J17E 斗l,7£-J4.3LL.72-J4.311.33-J3.4?5V55=-C,

13、6+tU(2)、形成导纳矩阵:YT 二6.46-jd 31-5.13+J277-13+口 5400-5.19+J2.7710.a3-j25.D9-1.27+J4.3J -1.27+j4.31 -2+卜-1 27tj3 54-1.72+j4.31P.69-j25.0P -6.70+j 17.240D-1.72+J4.31g.&j25-L.3S+j3.402+j60-1 3E+j34工把-p45(3)、B', B形成以及其逆矩阵。由题可知除1为平衡节点外，其他节点均为PQ节点，系数B', B阶数相同，又因对该等值网络不存在去除与有功功率和电压或无功功率和电压大小关系较小因数

14、的可能性，这两个矩阵B'，B完全相同，它们就由导纳矩阵的虚数部分中除第一行和第二行的各个元素组成。B' =B''25.094.314.3164.3125.0917.2404.3117.24253.45603.459.45由此可见，网络的节点导纳矩阵为奇异矩阵，但他的虚数部分的子矩阵B'和B，则是非记异矩阵，可以求逆矩阵，其逆矩阵:B1B'10.07680.05520.06110.07110.05520.11920.10170.07220.06110.10170.13270.08730.07110.07220.08730.1828、计算各节点有功功

15、率不平衡量VPi。取Ui1.06,1 0;U2(0) U3(0) U4(0) U5(0) 1.0;2(0)（0）3(0)45(0)0 ,按下式计算各节点有功率不平衡:VP(0)j nPiUi(0)Uj(0) (Gij cos ijj 1Bij sin ij)Vp(0) 0.251.0*1.06( 1.59o©2.77*sir0) 1.0*1.0*l0.63Co0 25.9ir01.O*1.O1.72Co0 4.31Sirf) 1.0*10 1.72£o0 4.3*sir01.0*1.*( 2cos06sin0) 0.584相似的可得:VP3（0）0.9Vp4(0)0.4 V

16、P5(0)0.6 ；（5）、计算各节点电压的相位角i由下列矩阵方程式：b' 1(vp(0)/U(0)(U(0)V(0)0.07680.0552 0.06110.07110.58141.0*(0)20.05520.1192 0.10170.07220.9 /1.01.0*(0)30.06110.1017 0.13270.08730.4 /1.0=1.0*(0)40.07110.0722 0.08730.18230.6 /1.01.0*(0)5从而可得：(0)20.0721(1)2(0)20.0721(0)30.1592(1)3(0)30.1592(0)40.1615; (1)4(0)40

17、.1615(0)50.1679(1)5(0)50.1679(6) 、计算各节点无功功率不平衡VQi按下列计算各节点无功功率不平衡:VQi(0)Qjj nUi(0)U.(0)i Jj 1(GjSin jBj cos j)心0)1.0 1.01.0 1.01.0 1.01.0 1.0VQ®0.46527V4®0.777443Q®0.757824;(7) 、计算各节点电压的大小Ui由下列矩阵方程式:(VQ (0) /U (0) ) VU (0)0.25 1.0 1.065.19 Sin( 0.07210) 2.77cos(0.0721 0)(10.63si n(0)25

18、.09cos0)1.72 sin(0.07210 0.15920)4.31 cos(0.072100.15920)1.72 sin(0.072100.1615)4.31 coS 0.072100.16150)2 sin(0.0721 0.16790)6 cos(0.0721 0.16790)0.264267相似可得:0.07680.05520.0611 0.711 7.264267U2(0)0.05520.11920.10仃0.0722 0.465275VU®0.06110.10170.13270.08730.777443VJ4®0.07110.07220.08730.18

19、28 0.757824U5(0)Vt(0)VU(0)V4(0)V5(0)0.6849620.5902290.6604900.756140UU(0)VU(0)0.315038UU2(0)VU(0)0.409771UU(0)V4(0)0.339510UU(0)V5(0)0.2438964、利用MATLABL语言或其他语言编程调试PQ分解法潮流程序/文件输入格式：节点总数n（包括联络节点）,支路数zls/节点数（发电机和负荷）nb,接地电抗数mdk,迭代精度eps/考虑负荷静特性标志kk2（0考虑）,平衡节点号，优化标志（0不优化）/最大迭代次数it1,支路左右节点号izl,jzl,支路电阻zr,/

20、支路容纳zyk,节点号nob及标志nobt（0-PQ-1-PV）/发电机和负荷有功、无功pg,qg,pl,ql/电压v0（pv节点输入实际值，PQ节点任输入一值）/ /电抗zx/电抗节点号idk,电抗值dkk #in elude "math.h"#i nclude "stdio.h"#defi ne NS#defi ne NS2#defi ne NS4#defi ne ZS#defi ne ZS2#defi ne DKS#defi ne N2#defi ne N32000NS * 21000 3000ZS * 2200ZS * 4ZS * 8 + NS

21、* 4/最大节点数/NS4、NS必须大于 2*zls。最大支路数/最大电抗器数FILE *fp1, *fp2; char inn ame12, out name12; / fp1输入数据文件指针fp2 /inname输入数据文件名输出文件指针outname输出数据文件名int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, n zls; / 荷） / / / / /节点总数n（包括联络节点） 支路数（回路数）zls 节点数nb（发电机和负接地电抗数mdk精度eps平衡节点号mpj节点优化（标志）bnsopton（=0节点不优化，！=0节点优化） 最

22、大迭代次数it1 最低电压或最大功率误差节点号dsd 负荷静特性标志（=0考虑负荷静特性）支路数（双回线算一条支路）int izlZS, jzlZS, idkDKS, ydsNS, ydzNS, iyZS2; / / / / /izl,jzl,idk:分别存放左、右节点号和电抗器节点号。yds存放各行非零非对角元素的个数。ydzi是第i行第一个非零非对角元素的首地址， 即在所有非零非对角元素中的次序号iy 存放列足码。int nn ewNS4, oldNS, nobNS, no btNS; / / / 点。double eps, dsm, vmin, dph, dqh, af3;/eps迭代收

23、敛精度，dsm最大功率误差/nn ew,old存放的是新、旧节点号。nnewi中为i对应的新号nob存放的是节点号。nobt存放的是节点类型，0: pq 节点,-1: pv 节vmin:系统最低电压值。dph,dqh:系统有、无功损耗。/af0和af1分别是负荷有功功率、无功功率静态特性系数。double v00;/ v00:系统平均电压 ci,cj分别作为节点i,j7元0的电压相角的临时存储单double zrZS, zxZS, zykZS, dkkDKS, giiNS, biiNS, ygZS2, ybZS2;double p gNS, qgNS, plNS, qlNS, v0NS, vN

24、S, vaNS;支路电阻zr支路电抗zx输电线路充电容纳zyk(y0/2)接地电抗dkk对角元实部gii对角元虚部非对角元实部yg非对角元虚部ybpg,qg, pl,ql：发电机，负荷功率实、虚部v是电压幅值,va是电压相角。/ double wNS2, kg3, bNS2;int n ewsortNS4;/ newsorti 存放i对应的老号void in itial();void p qflow();void out();void dataio();void bnsop t();void zlsort(i nt* nn ew);void prin to();void prin ty();v

25、oid y2();void yaO();void yzb();void jdgl(i nt kq0);void bbhl(i nt kq0);void calc();int iabs(i nt a);void bran ch_out pu t();void n ewval(double* aa);void prin tc();void iswa p();void swa p();void prin tf2(double* aa, double* bb, int n);void calc(i nt* iu, double* u, double* di, i nt* nfd, double* b)

26、;void prin ti(i nt* aa, int n);void prin tf1(double* aa, int n);int fin d(i nt k, int a, int* z);void yzb(i nt t, i nt* iu, double* u, double* di, i nt* n fd);int isg n(i nt a, int b);void yy1();void y3();void n ewtoold();int main (void)/初始化/pq潮流计算/输出节点和支路数据in itial(); P qflow(); out();return 1;int

27、isg n(i nt a, int b)/* 本函数功能返回值为a的绝对值b的符号 */ /参数1提供值，参数2提供符号/if (b < 0)if (a > 0) a = -a;return a;int fin d(i nt k, i nt a, i nt* z)*/*本函数查找a中是否有fabs(k)有则返回0,无则返回1/参数1为待查找量，参数2待搜索数组，参数3返回k在a中的次序号/int i;for (i = 1; i <= n; i +) if(iabs(k) = ai) *z = i; return 1;return 0;void oldt on ew()/*本函

28、数将输入数据中的节点号变成从1开始的连续节点号*/int i, j, k, ii1, ii2, zls2, k1, k2, k3, k4, i p; zls2 = zls + zls;for (i = 1; i <= zls2; i +) n ewsorti = 0;ii1 = 0;for (i = 1; i <= zls; i +) k = izli;if (!fin d(k, newsort, & ii2)ii1 +;n ewsortii1 = iabs(k);k = jzli;if (!fin d(k, newsort, & ii2)ii1 +;n ewsor

29、tii1 = iabs(k);for (i = 1; i <= ii1-1; i +)for (j = i+1; j <= ii1; j +)if (n ewsorti > n ewsortj)k = n ewsorti;n ewsorti = n ewsortj;n ewsortj = k;for (i = 1; i <= zls; i +)k = izli;if (fin d(k, newsort, & ii2)izli = isg n( ii2, k);elseprin tf("error!");k = jzli;if (fin d(k

30、, newsort, & ii2)jzli = isg n( ii2, k);elseprin tf("error!");prin tf("izl%d = %d, jzl%d = %dn", i, izli, i, jzli);for (i = 1; i <= nb; i +)for (j = 1; j <= n; j +)if (n obi = n ewsortj)nobi = j;break;prin tf(" nob%d = %dn", i, no bi);for (j = 1; j <= n; j +

31、)if (mpj = n ewsortj) mpj = j；break;电抗器节点号转变mdk; j +)/for (j = 1; j <=for (i = 1; i <= n; i +) if (idkj = n ewsorti) idkj = i;break;void in itial() /*本函数进行初始化工作*/int i, k1;dataioO;/输入原始数据oldto new();/转化为新号if (bn so pton = 0) /节点不优化，新节点号即为老节点号。for (i = 1; i <= n; i +)oldi = i; nn ewi = i;els

32、ebn so pt();/节点优化mpj = nn ewmpj;/mpj:平衡节点zlsort (nn ew); / sort the r,x and b for (i = 1; i <= mdk; i +)k1 = idki;idki = nn ewk1;for (i = 1; i <= n; i +)vi = vOO;vai = 0.0; /所有节点的电压幅值初值都为1.000(v00)，电压相角初值都为0。/ excha nge the node before and after sortfor (i = 1; i <= n; i +)ydsi = 0;/ the im

33、mediatefor (i = 1; i <= nb; i +)k1 = nnewn obi;ydsk1 = nobti;for (i = 1; i <= n; i +)n obti = ydsi;n ewval( pg);n ewval(qg);n ewval( pl);n ewval(ql);n ewval(v0);for (i = 1; i <= n; i +) / no bt is type of nodeif (n obti = -1)vi = v0i;/ no b is serials nu mbe/nobt = -1: pv 节点，v0存放的是最后一个节点数据，

34、/对于pv节点，即为该点应维持的电压值。/nobt = 0: pq 节点，v0存放的是最后一个节点数据，/对于pq节点，即为系统平均电压值。prin to();/输出af、vOO和节点排序后的支路、节点和/接地电抗数据(仅仅查看中间结果)yaO();/ 获得 yds、ydz、列足码 iy 。( P407 ) void prin to() /* 输出af、vOO和节点排序后的支路、节点和接地电抗数据*/int i;fprin tf(fp2, "n *AF AND VO *n");fprin tf(fp2, "n %7.3f%7.3f%7.3fn", afO

35、, af1, vOO);prin tc('-', 78);fprin tf(fp2, "nn *zLB*n");for (i = 1; i <= zls; i +)fprin tf(fp2, "n");fprin tf(fp2, "%8d%8d%8d%8d", izli, jzli, oldabs(izli), oldabs(jzli);fprin tf(fp2, "%9.4f%9.4f%9.4f", zri, zxi, zyki);prin tc('-', 78);fprin

36、 tf(fp2, "nn*BUS*n");for (i = 1; i <= nb; i +)fprin tf(fp2, "n");fprin tf(fp2, "%8d%8d%8d", no bi,old nobi, nobti);fprin tf(fp2, "%9.4f%9.4f%9.4f%9.4f%9.4f", p gi, qgi, pli, qli, vOi);prin tc('-', 78);fprin tf(fp2,"nn*DKK* n"); for (i = 1;

37、 i <= mdk; i +)fprin tf(fp2, "n");fprin tf(fp2, "%8d%8d%7.4f', idki, oldidki, dkki);void dataio()/* 系统数据初始化 */int i;af0 = 0.6;af1 = 2.0;/af0 和af1分别是负荷有功功率、无功功率静态特性系数。 v00 = 1.000;/ 系统平均电压prin tf("np lease input the n ame of data filen");sca nf("%s", inn ame)

38、;fp1 = fopen (i nn ame, "r");prin tf("np lease out put the n ame of data filen");sca nf("%s", out name);fp2 = fopen(outn ame, "w");fsca nf(fp1, "%d %d %d %d", &n, &zls, &nb, & mdk);/ the nu mber of node ,bra nches, nodefsca nf(fp1, &q

39、uot;%lf %d %d %d %d", &eps, & kk2, & mpj,&bn so pto n, & t1);/p recisi on, swing no de,sort the no de,iterati on nu mbers for (i = 1; i <= zls; i +)fscan f(fp1, "%d %d", & izli, &jzli);fscan f(fp1, "%lf %lf %lf ", &zri, & zxi, &zyki

40、);for (i = 1; i <= nb; i +)fscan f(fp1, "%d %d", &n obi, &n obti);fscan f(fp1, "%lf %lf %lf %lf %lf", &p gi, & qgi, &p li, &qli, &v0i);for (i = 1; i <= mdk; i +)fscan f(fp1, "%d %lf", & idki, & dkki);fclose(fp1); void p qflow()i

41、nt kqO, iu1N2, nfd1NS, iu2N2, nfd2NS;int i, t;double u1N2, u2N2, di1NS, di2NS;yyi();yzb(0, iu1, u1, di1, nfd1); /form the B matrix of P-0 iterationy2();yzb(1, iu2, u2, di2, nfd2); /form the B matrix of Q-V iteration t = 0;kq0 = 0;kg0 = kg1 = 1;dojdgl(kq0); / calculat ing the po werbbhl(kq0); / find

42、out the maxi if (kq0 = 0)printf("P: %dt%dt%fn", t, dsd, dsm);elseprin tf("Q: %dt%dt%fn", t, dsd, dsm);if (fabs(dsm) > eps) kgkqO=1; if (kqO = 0) calc(iu1, u1, di1, n fd1, b);if (kqO = 1)calc(iu2, u2, di2, nfd2, b); for (i = 1; i <= n; i +) if(kq0 = 0 )vai = vai - bi / v00;

43、elsevi = vi - bi;elsekgkq0 = 0;if(kq0 = 0)kq0 = 1;elsekq0 = 0;t +;if(t > it1)break;while(fabs(dsm) > eps) | (kgkq0 != 0); fprin tf(fp2, "n%s%d", "times = ", t);void out() /*本函数输出节点和支路数据*/zlsort(old); / recover the data if sorted /n ewtoold();no de_out pu t(); / node data br

44、an ch_out pu t(); /branch data prin tc('-', 78);prin tc('*', 78);fprin tf(fp2, "n"); void n ewval(double* aa) /* 本函数将旧号换成新号*/int i, k1;for (i = 1; i <= n; i +) bi = 0.0;for (i = 1; i <= nb; i +) k1 = nnewn obi; bk1 = aai;for (i = 1; i <= n; i +) aai = bi;void yzb(i

45、 nt t, i nt* iu, double* u, double* di, i nt* nfd) /* 本函数求因子表 */参数1为标志（t=0 求B',t=1 求B''）/参数2因子表上三角矩阵非零非对角元素的列足码参数3因子表上三角矩阵非零非对角元素的数值 参数4因子表上三角矩阵对角元素参数5因子表上三角各行非零元素个数int i, j, k, i1, i2;int jj, jj1, jj2, im, x, fdNS;double ai, bNS;n fd1 = 1;for (i = 1; i <= n; i +)节点。存在(t = 1)/nobt存放的是

46、节点类型，0: pq 节点,-1: pvif (t != 1) II (n obti != -1) && i != mpj) / <-| / |for (j = i + 1; j <= n; j +)/bj = 0.0;/|bi = biii;/|if (kk2 = 0) && (t = 1) && (no bti != -1)/的情况，不多余。bi = bi + af1 * qli / v0i / v0i;/af111 = ydzi;12 = ydzi + 1 - 1;for (j = i1; j <= i2; j +)k =

47、 iyj;bk = ybj;bmpj = 0.0;if (t = 1) for (j = 1; j <= n; j +)if (n obtj = -1) bj = 0.0;i1 = i - 1;for (im = 1; im <= i1; im +)jj1 = n fdim;jj2 = nfdim + 1 - 1;for (jj = jj1; jj <= jj2; jj +)if(iujj = i)ai = ujj / diim;for(k = jj; k <= jj2; k +) j = iuk;bj = bj - ai * uk; break;(三=7启)壬一启(L

48、 HHg&w(+ 二u Hv 二 L 上)04 二u 一U丄上兰ee SHM(u lu一 Ee *e-qnop)E£d p一0>宀 (=*-CI * 5= 7启)主启 宀 二二n二二n 送g oL&= -g)主启(三=7启)壬一启(L HHe&w(+ 二X HV 二 L 上)04-(=*=* 5= 7启)主启宀兰E +曰EU丄L+二EU宀6.0 hnp6丄二ES-宀->1 丄二E宀宀匸+X_+土J丄合* Bq H 2n(gTeo L A (mq)sq国二一(+u Hv匸H DO4二 + 一 H 匸兰一P上e 兰q、O L hnp兰 Eu H XII

49、fprin tf(fp2, "%9.5f", aai); fprin tf(fp2, "nn"); void calc(i nt* iu, double* u, double* di, i nt* nfd, double* b)II*本函数利用因子表解线形方程组。(P417图F1-9)*ll前代过程。回代过程。double bi;int i, j, k, i1, i2;for (i = 1; i <= n; i +) IIbi = bi;11 = n fdi;12 = n fdi + 1;for (j = i1; j < i2; j +) k

50、 = iuj;bk = bk - bi * uj;bi = bi * dii;for (i = n; i >= 1; i -) IIbi = bi;11 = n fdi;12 = n fdi + 1 - 1;for (j = i2; j >= i1; j -) k = iuj;bi = bi - bk * uj;bi = bi;void zlsort(i nt* nnew)Illi*IIII* 本函数进行支路数据排序II小节点号放左边，大节点号放右边II左右皆按从小到大顺序排列/int ip, k1, k2, k3, k4;int i, j；for (i = 1; i <=

51、zIs; i +) k3 = izli;k4 = jzli;k1 = iabs(k3);k2 = iabs(k4);/原节点号。izli = isg n(nn ewk1, k3); /新节点号。jzli = isg n(nn ewk2, k4);k3 = izli;k4 = jzli;k1 = iabs(k3);k2 = iabs(k4);if (k1 > k2)izli = k4;jzli = k3;for (i = 1; i <= zis - 1; i +) ip = i;k1 = iabs(izli);k3 = iabs(jzli);for (j = i + 1; j <= zls; j +)k2 = iabs(izlj);k4 = iabs(jzlj);if(k2 < k1 | (k2 = k1 && k4 < k3) ip = j;k1 = k2;k3 = k4;if(i != ip)iswa p(&izli, & izli p);iswa p(&jzli, &jzl