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

1、电气工程及其自动化专业综合实训一一. 原始资料网络接线如图，各支路阻抗和各节点功率均以标么值标于图中， 其屮节点2连接的实际是发额定功率的发电厂。设节点1的电压 保持为1.06,试计算图示系统中的潮流分布。（计算精度10，二、等值电路图二、报告要求：1、写出潮流计算的原理2、写出功率方程（极坐标或直角坐标）3、写出雅兑比矩阵元素的表达式或pq分解法屮的系数矩阵4、利用matlable语言或其他语言编程调试5、心得或体会6、附程序代码，打印潮流计算步骤：1、写出潮流计算的原理潮流是指在发电机母线上功率被注入网络，而在变（配）电站的 母线上接入负荷，其间，功率在网络中流动。对于这种流动的功率， 电

2、力生产部门称之为潮流。以电力网络潮流、电压计算为主要内容的 电力网络稳态行为特性计算的目的在于估计对用户电力供应的质量以 及为电力网运行的安全性与经济性评估提供基础数据。配电网潮流计 算是配电网络分析的基础，配电网的网络重构、无功功率优化、状态 估计和故障处理都需要用到配电网潮流数据。电力系统稳态运行应满足以下要求：1）满足系统经济性运行的要求，每一台发电机的输出必须接近于 预先设定值；2）必须确保联络线潮流低于线路热极限和电力系统稳定极限；3）必须保持某些中枢点母线上的电压水平在容许范围内，必要时 用无功功率补偿计划来达到；4）区域电网是互联系统的一部分，必须执行合同规定的输送至邻 网的联络

3、线功率计划；5）用故障前的潮流控制策略使事故扰动效应最小化。通常情况下，输电线路电压在轻载时会较高，重载时会较低，电 压调整是指在负载由轻载到满载变化过程中实时调整线路电压满足运 行要求；对于超高压输电线路，线路电压维持在额定电压的±5%之内， 实际运行吋，通常电压调整约为10%。对于低压输电线路，电压调整 数值为10%，包含了变压器本身的电压降落。3. 1. 1潮流计算的基本物理量潮流计算是电力系统分析中的一种最基本的计算，它的任务是对 给定的运行条件确定系统的运行状态，就是在三相平衡稳态状态下计 算电力系统中每条母线的电压幅值和相角，其中每一设备如传输线和 变压器中的有功和无功潮

4、流，以及各设备的损耗都需要计算出来。潮流计算采用电力系统的单线图，对于任意一条母线i，需要以下 四个变量描述：电压幅值ui、相角电网供给母线的有功pi、无功 qi。若某一电力系统有n个节点，则共有4n个变量，对于每条母线， 这些变量中的两个指定为输入数据，其它的两个是潮流程序所要计算 的未知量。为方便起见，在图3. 1中传送给母线i的功率可分为发电 机发出和负载吸收两部分。也就是pi = pgi - pli每条母线被归分为以下三种母线类型中的某一类:1）平衡节点，一般一个系统只有一个平衡节点。在潮流分布 算出以前，网络中的功率损耗是未知的，因此，至少有一个节点的 有功功率p和无功功率q不能给定

5、。另外必须选定一个节点，制定 其电压相角为零，作为其它节点电压相位的参考，这个节点叫基准 节点。为了计算方便，常将平衡节点和基准节点设在同一个节点上。 为方便起见在本书中把它标号为母线1。平衡节点是电压参考节点， 该母线的7是给定值，作为输入数据，典型取标么值i&q、潮 流程序计算p1和q1。因为平衡节点的p、q事先无法确定，为使潮 流计算结果符合实际，常把平衡节点选在有较大调节裕量的发电机 节点，潮流计算结朿时若平衡节点的有功功率、无功功率和实际情 况不符，就要调整其他节点的边界条件以使平衡节点的功率满足实 际允许范围。2）pq节点，pi和qi是输入数据。这类节点的有功功率pi和 无

6、功功率qi是给定的，潮流计算程序计算节点电压幅值ui和相角 5:。负荷节点和无功功率注入的联络节点都属于这类节点。有些情 况下，系统中某些发电厂送出的功率在一定吋间内为固定时，该发 电厂母线也可以作为pq节点。在一个典型的潮流程序中绝大多数 母线作为pq节点。3）pu节点（电压控制母线），pi和ui是输入数据。这类节点 的有功功率pi和节点电压幅值ui是给定的，潮流程序计算节点的 无功功率qi和电压相角这类节点必须具有足够的无功可调容量，用以保持给定的节点电压幅值。在电力系统中这类节点的数目较少。例如与发电机、并联补偿电容器或者静止无功系统相连的母 线。设备无功功率最大值qgimax和最小值q

7、gimin都是输入数掘。 兒一个例子是与抽头可调节变压器相连的母线，用潮流程序计算抽 头的位置。注意，当母线i是无发电机相连接的负载母线时，pi=-plj为 负值；也就是说，在图3. 1中给母线i提供的有功为负值。如果负 荷是感性的，qi=-qli为负值。综上所述，若系统中有n个节点，n为网络中除去参考节点外 的节点数，本书中以大地为参考节点，选第1个节点为平衡节点， 剩下的n-1个节点中有r个pu节点，则有n-r-1个pq节点。因此， 除了平衡节点外，有n-1个节点的注入有功功率、n-r-1个pq节点 的注入无功功率和r个pu节点电压幅值为已知量。3.1.2潮流计算的数学模型在稳态潮流计算屮

8、，电力系统各元件（参数）等效成一个有源 网络。将发电机和负荷用无阻抗线从网络中抽出，剩不的是由接地 和不接地支路组成的无源线性网络，可以用导纳矩阵（y）或阻抗 矩阵（z）来描述。采用导纳矩阵时，节点电流和节点电压构成以下方程：i = yu （3-1）其中：y为nxn阶导纳矩阵，其阶数n为网络中除去参考节点 外的节点数，如果不考虑网络元件的非线性及变压器的相位偏移，y为对称矩阵，i为nxl维节点注入电流列向量，在电力系统计算中，节点注入电流可理解为该节点电源电流与负荷电流之和，并规 定流入节点电流为正。因此仅有负荷的节电电流就为负值，某些仅 起联络作用的节点，图3-2中节点n=3,其注入电流为零

9、。u为nxl 维节点电压列向量。网络中有接地支路时，节点电压通常指该节点 的对地电压，以人地作为参考节点，并规定其编号为零。图3-2运用节点电压法时的电力网路等值电路 对于第i个节点，展开为如下形式：=(/=1，2，"、门)(3-2)若采用阻抗矩阵可表示为：u=zi展开为：斉(3-3)在潮流计算时一般以节点电压方程进行。节点导纳矩阵与阻抗 矩阵互为逆阵，在短路计算时可直接利用导纳矩阵求逆得到阻抗矩 阵以求得短路点的短路电流。由于实际系统中一般不给出节点电流而是节点功率，因此将式(3-2)中的节点注入电流用节点注入功率来表示为:si = pjq = uj；=ui(3-4)如果节点电压用

10、极坐标表示，令 ui = u,= u,h =q.+jb.j=u,"，乃n个节点电力系统的潮流方程的一般形式是=pi+jq 川知严，)>1，2，.，(3-5)>1j-l或zujco叫-中 +乓sin(4-)_q =uizujcos( -.)(3-6)若采用直角坐标系，节点电压可表示为导纳矩阵元素可以表示为a = q + j马将上述表达式带入式(2-8)的右端，展开并分出实部和虚部,便得二ms（似.j-l(3_?)可见，原来电流电压的线性方程组变换为功率和电压的非线性 方程组，式(3-6) (3-7)就是潮流计算的基本方程。它是一组共 wn个非线性方程组成的复数方程组，如果把

11、实部和虚部分开便得 到2n个实数方程，巾该方程组讨解出2n个运行参数。但是每一个 节点都有p、q、u、四个运行变量，共有4n个运行参数，所以要 事先给定其余2n个参数。这就要根据节点的分类，将每个节点的4 个运行参数中的两个作为原始数据，另外两个作为待求量。3. 1.3潮流计算的约束条件为了保证电力系统的正常运行，潮流问题中某些变量应满足一 定的约束条件，常用的约束条件有：(1) 所有节点电压必须满足脈 o =从保证电能质量和供电安全的要求来看，电力系统的所有电气 设备都必须运行在额定电压附近。pu节点的电压幅值必须按上述条 件给定。因此，这一定约束条件主要是对pq节点而言。(2) 所有电源节

12、点的有功功率和无功功率必须满足max qgi 2gimaxpq节点的有功功率和无功功率以及pu节点有功功率，在给定吋就必须满足此条件。因此，对平衡节点的p和q以及pu节点的q 应按此条件进行检验。(3) 某些节点之间电压的相位差应满足为了保证系统运行的稳定性，要求某些输电线两端的电压相位差不 超过一定的数值。因此，潮流计算可以归结为求解一组非线性方程 组，并使其解答满足一定的约朿条件。如果不能满足，则应修改某 些变量。甚至修改系统的运行方式，重新进行计算。2、写出功率方程(极坐标或直角坐标)3、写出雅克比矩阵元素的表达式或pq分解法中的系数矩(2)、形成导纳矩阵6.46-j6.31-5.19+

13、j2.77-1.27+j3.540-5.19+j2.77-1.27+j3.5410.63-j25.09 -1.72+j4.31 -1.27+j4.319.69-j25.09-1.27+j4.31-6.70+j17.24-2+j600-1.72+j4.31 -6.70+j 17.24 9.8-j25 -1.38+j3.40-2+j60-1.38+j3.43.38-j9.45(3)、b' , b形成以及其逆矩阵。巾题可知除1为平衡节点外，其他节点均为pq节点，系数b'，b" 阶数相同，又因对该等值网络不存在去除与有功功率和电压或无功功 率和电压大小关系较小因数的可能性，这

14、两个矩阵b' , b"完全相 同，它们就由导纳矩阵的虚数部分中除第一行和第二行的各个元素组-25.094.3 14.3 164.3 1-25.0917.2404.3 117.24-253.45603.45-9.45=b"由此可见，网络的节点导纳矩阵为奇异矩阵，但他的虚数部分的子 矩阵b'和b"，则是非记昇矩阵，可以求逆矩阵，其逆矩阵:-0.07 6 8-0.05 5 2- 0.0 6 1 1- 0.07 1 1- 0.0 5 5 2 -0.06 1 1 -0.07 1 1-0.1192 -0.10 17 -0.07 22-0.10 17 -0.13

15、27 - 0.0 8 7 3-0.0 7 2 2 - 0.0 8 7 3 -0.1828（4）、计算各节点有功功率不平衡量p，。取 q =1级4 =q"2(0) =u ="4(0) ="5(o) =i.q4(0)=岑0) =4(0) =4(0) =o，按 下式计算各节点有功率不平衡：叫(°) ="/",-(g" cos4 + % sin)7=1口 乃(。)=0.25-1.0*1.06- (-1.59taso-4 -277)*0)-lo*w*(10.63*aso+25.9咖0) 一1.0*1.0(1.72*cr«o+

16、4.31* 咖0) -1.0*1.0(1.72*6r»o+4.31*50) -l0*1.0*(-2oos0+6sin0) =0.5814相似的可得：凡=-0即凡=-0.4j：朽=0.6;（5）、计算各节点电压的相位角5；:由下列矩阵方程式：b r1 旧户。）/f/o）=（/口汐。）一-0.0768 -o.q552 -0.0611 -0.071 f_ 0.5814"l.o* j2(0)"-0.0552 -0.1192 -0.1017 -0.0722-0.9/1.01.0*么(0)-0.0611 -0.1017 -0.1327 -0.0873-0.4/1.0二1.0*

17、 j4(0)-0.0711 -0.0722 -0.0873 -o.1823_-0.6/1.0lo*(o)_从而可得:"-0.0721"、-0.07210)-0.1592參+-0.1592么(0)-0.16155么么(0)-0.1615-0.16795(1)-0.1679(6)、计算各节点无功功率不平衡ne,。按下列计算各节点无功功率不平衡：j=n 2,(°) = a - s fa(0)c/°) - (a sin 氏-' cos 氏)7=122(0)=0.251.0xl.06xx5.19xsin(-0.07210)-2.77cos(-0.0721

18、-0) 一 1 .oxl.ox (10.63sin(0) - 25.09cos 0) -l.oxl.ox1.72xsin(-0.07210+0.15920)4.31xcos(-0.07210+0.15920) 一 l.oxl.ox1.72xsin(-0.07210+0.1615) 4.31xcos(-0.07210+0.16150) 一 l.oxl.ox-2xsin(-0.0721+0.16790)-6xcos(-0.0721+0.16790) =-0.264267相似可得:q(o) =-o.465275og(o) =-q777wq(o) =-0.757824;(7)、计算各节点电压的大小。巾

19、下列矩阵方程式：-(b')_1(e(o)/"(o) =口/(0)-7.25w7-q4s275-q7774i3_-q757824-q(»s8 -qq552 -q0611 -q711 -qq552 -q1192 -q1017 -qo722 -q0611 -q1017 -q1327 -q0873 -q0711 -qo722 -q0673 -q1828f/0)r-q6s«s-q39qb29xy<0)-qffiow_-q756140_+f/2(0)q315cb8q409771巧)ifq33951oq243s64、利用matlable语言或其他语言编程调试pq分

20、解法潮流程序/文件输入格式：节点总数n （包括联络节点），支路数zl s/节点数（发电机和负荷）nb，接地电抗数mdk，迭代精度ops/考虑负荷静特性标志kk2 （0考虑），肀衡节点号，优化标志（0不优化）/最大这代次数itl,支路左右节点号izl, jzl，支路电附zr，电抗zx/支路容纳 zyk,节点号 nob及标志 nobt （o-pq, -1-pv）/发电机和负荷有功、无功pg, qg, pl, ql 口/电压v0 （pv节点输入实邮值，pq爷点任输入一值）/电抗节点号idk,电抗值dkk/ttdefine ns define ns2 #dcfinc ns4 define zs def

21、ine zs2 ttdefine dks define n2 #dcfinc n3include "math, h" include "stdio. h"2000/最大节点数ns * 21000/ns4、ns 必须大于 2*zls。3000/最大支路数zs * 2200/最大电抗器数zs * 4zs * 8 + ns * 4file *fpl, *fp2;char inname12, outname12;/ fpl输入数据文件指针fp2输出文件指针/ inname输入数掘文件名outname输出数掘文件名int n, zls, nb, mdk, mpj,

22、 bnsopton, itl, dsd, kk2, nzls;/节点总数n（包拈联络节点）支路数（冋路数）zls节点数nb（发电机和负 荷）/接地电抗数mdk精度eps 平衡节点号mpj/节点优化（标志）bnsopton（=0节点不优化，！ =0节点优化）/最大迭代次数itl最低电压或最大功率误差节点号dsd/负萜静特性标志（=0考虑负萜静特性）/支路数（双回线算条支路）int izlzs, jzlzs, idkdks, ydsns, ydzns, iyzs2;/ izl, jzl，idk:分别存放左、右节点号和电抗器节点号。/ yds存放各行非零非对角元素的个数。/ ydzi是第i行第一个非

23、零非对角元素的首地址，/即在所有非零非对角元素屮的次序号/ iy口存放列足码。int nnewns4, oldns, nobns, nobtns;/ nnew, old存放的是新、旧节点号。/ nnew i 中为i对应的新号/ nob存放的是节点号。nobt存放的是节点类型，0: pq节点，_1: pv节 点。double eps, dsm, vmin, dph, dqh, af 3;/ eps这代收敛精度，dsm最大功率误差/ vmin:系统最低电压值。dph，dqh:系统有、无功损耗。/ af0和afl分别是负荷冇功功率、无功功率静态特性系数。double v00;/_/ v00:系统t均

24、电压ci,cj分别作为节点i，j的电压相角的临时存储单yg odouble zrzs, zxzs, zykzs, dkkdks, giins, biins, ygzs2, ybzs2;double pgns, qgns, pins, qlns, vons, vns, vans;/支路电阻zr支路电抗zx输电线路充电容纳zyk (yo/2)/接地电抗dkk对角元实部gii对角元虚部/非对角元实部yg非对角元虚部yb/ pg, qg, pl, ql:发电机，负荷功率实、虚部/ v是电压幅值，va是电压相角。double wns2, kg3, bns2;int newsortns4;/ newsor

25、ti存放i对应的老号void initial 0;void pqflowo ;void out ();void dataioo ;void bnsopt ();void zlsort (int* nnew);void printoo ;void printy ();void y2 ();void yao();void yzb ();void jdgl (int kqo);void bbhl (int kqo);void calc 0;int iabs(int a);void branch_output();void newval(double氺 aa);void printc ();void

26、iswap();void swap();void printf2 (double* aa, double氺 bb， int n);int find(int k， void yzb(int t, int isgn(int a， void yyl0; void y30 ; void ncwtoold()void calc (int氺 iu， double氺 u， double氺 di, int氺 nfd， double氺 b); void printi (int* aa，int n); void printf1(double* aa, int n);int a口， int* iu, int b);

27、int* z);double氺 u， double* di, int* nfd);int main(void)initial () ; /初始化pqflowo ; /pq潮流计算out () ;/输出节点和支路数据return 1;int isgn(int a，int b)/*木函数功能返回值为a的绝对值b的符号 */参数1提供值，参数2提供符号/if (b < 0) if (a > 0) a = -a;return a;int find (int k，int a，int* z)/*本函数查找a中是否有fabs(k)有则返冋0，无则返冋1*/参数1为待查找量，参数2待搜索数组，参数

28、3返回k在a屮的次序号/int i;for (i = 1; i = n; i +)if(iabs(k) = ai)*z = i; return 1;return 0;void oldtonewo/*本函数将输入数据中的节点号变成从1开始的连续节点号int i, j, k, iil, ii2, zls2, kl, k2, k3, k4, ip; zls2 = zls + zls;for (i 二 1; i 二 zls2; i +) ncwsorti = 0;iil = 0;for (i 二 1; i <二 zls; i +)k = izli;if (!find(k, newsort, &a

29、mp;ii2)iil +;newsortiil = iabs(k);k = jzli;if (!find(k, newsort, &ii2)iil +；ncwsortiil = iabs(k);for (i = 1; i 二 iil-1; i +)for (j = i+1; j <= iil; j +)if (newsorti > newsortj)k = newsorti; newsorti = newsortj; newsortj = k;for (i = 1; i <= zls; i +)k 二 izli;if (find(k, newsort, &ii

30、2)izl i = isgn(ii2, k);elseprintf (error!)； k =if (find(k, newsort, &ii2)jzli = isgn(ii2, k);elseprintf (error! )；printf ("izl %d = %d, jzl %d = %dn", i, izli, i, jzli);for (i = 1; i <= nb; i +)for (j = 1; j <= n; j +) if (nobi = newsortj)nobi = j; break;printf ("nob %d = %d

31、n", i, nobi);for (j = 1; j二 n; j +)if (mpj = newsortj)mpj = j； break;/电抗器节点号转变for (j 二 1; j <= mdk; j +)for (i = 1; i <= n; i +)if (idkj = newsorti)idkj = i; break;void initial ()/*木函数进行初始化工作*/int i, kl;dataioo ;/输入原始数据 oldtoncwo ;/转化为新号if (bnsopton = 0)/节点不优化，新节点号即为老节点号。for (i = 1; i <

32、;= n; i +)oldi = i; nnewi = i;elsebnsopt () ;/节点优化mpj : nnewmpj ;/mpj:平衡节点 zlsort (nncw); / sort the r, x and b for (i = 1; i <= mdk; i +)kl 二 idki; idki = nnewkl;for (i = 1; i <= n; i +)vi二 voo; vai = 0. 0; /所宥节点的电压幅值初值都为l.ooo(voo),电压相角初值都为0。 / exchange the node before and after sortfor (i =

33、1; i <= n; i +)ydsi二 0;/ the immediatefor (i = 1; i <= nb; i +)kl = nnewnobi; ydskl = nobti;for (i = 1; i <= n; i +)nobti = ydsi;newval (pg);newval (qg);newval(pi);newval (ql);newval (vo);for (i 二 1; i <= n; i +)/ nobt口 is type of nodeif (nobti=二 一1)vi二 voi ;/ nob is serials numbe/nobt =

34、 -1: pv节点，v0存放的是最。一个节点数掘， /对于pv节点，即为该点应维持的电压值。/nobt = 0: pq节点，v0存放的是最后一个节点数据， /对于pq节点，即为系统平均电压值。 printoo ；/输出af、v00和节点排序£；的支路、节点和/接地电抗数据（仅仅查看中间结果）ya0（） ;/获得 yds 口、ydz、列足码 iy。（ p407 ）void printoo/*输出af、v00和节点排序后的支路、节点和接地电抗数据*/int i;fprintf(fp2, n 料and vo 料*n"); fprintf (fp2, "n %7. 3f%

35、7. 3f%7. 3fn", af0, afl, voo); printc ，78);fprintf (fp2, nn*7丄b*n");for (i = 1; i <= zls; i +)fprintf (fp2, "rt);fprintf(fp2, 8d%8d%8d%8d", izli, jzli, oldabs(izli), oldabs(jzli);fprintf (fp2, 9. 4f%9. 4f%9. 4f，zri, zxi, zyki);printc c，78);fprintf (fp2,；for (i = 1; i <= nb;

36、 i +)fprintf(fp2, "rt);fprintf(fp2, "%8d%8d%8d", nobi,oldnobi，nobti); fprintf (fp2, ，z%9. 4f%9. 4f%9. 4f%9. 4f%9. 4f，z, pgi, qgi, pli,qli，voi);printcc，78);fprintf (fp2, "nn*林林dkk*n"); for (i = 1; i <= mdk; i +)fprintf(fp2, "rt);fprintf(fp2, "%8d%8d%7.4f", i

37、dki, oldidki, dkki);void dataioo/*系统数据初始化*/int i;af0 = 0.6;afl =2.0;/af0和afl分别是负荷有功功率、无功功率静态特性系数。 v00 = 1.000;/系统平均电压printf(nplease input the name of data filen)；scanf (s，inname);fpl = fopen(inname, r)；printf(nplease output the name of data filen)；scanf (s，outname);fp2 = fopen(outname, ”w”)；fscanf (

38、fpl, "%d %d %d %d", &n, &zls, &nb, &mdk);/ the number of node ,branches, nodefscanf (fpl, /z%lf %d %d %d %d", &eps, &kk2, &mpj,&bnsopton, &itl);/precision, swing node,sort the node, iteration numbers for (i = 1; i <= zls; i +)fscanf (fpl, "%

39、d %d", &izli, &jzli);fscanf(fpl, lf %lf %lf ", &zri, &zxi, &zyki);for (i = 1; i <= nb; i +)fscanf(fpl, "%d %d", &nobi, &nobti);fscanf(fpl, lf %lf %lf %lf %lf", &pgi, &qgi, &pli,&qli， &v0i);for (i = 1; i <= mdk; i +)fscanf

40、(fpl, "%d %lf", &idki, &dkki);fclose(fpl);void pqflow()int kqo, iuln2, nfdlns, iu2n2, nfd2ns;int i, t;double uln2, u2n2, dilns, di2ns;yyl 0 ;yzb(0, iul, ul, dil, nfdl); /form the b matrix of p-0 iteration y2 0 ;yzb(1, iu2, u2, di2, nfd2); /form the b matrix of qv iteration t = 0;kq

41、o = 0;kg0 = kgl = 1;dojdgl (kqo); / calculating the power bbhl(kqo); / find out the maxi if (kqo = 0)printf("p: %dt%dt%fn", t, dsd, dsm); elseprintf(q: %dt%dt%fn", t, dsd, dsm);if (fabs (dsm) > eps)kgkqo:l; if (kqo = 0)calc (iul，ul,dil, nfdl,b);if (kqo = 1)calc(iu2,u2,di2, nfd2,b);

42、for (i = 1; ij二n; i +)iif (kqo =0 )vai=vai - bi/ voo;elsevi=vi- bi;elsekgkqo = 0;if(kqo = 0) kqo = 1;elsekqo = 0; t +;if(t > itl) break;while(fabs(dsm) > eps) | | (kgkqo != 0); fprintf(fp2, n%s%d", "times = ", t);void out ()/*本函数输出节点和支路数据*/zlsort (old);/ recover the data if sort

43、ed/ newtooldo ;nodc_output();/ node databranch_output(); /branch data printcc，78);printcc*，78); fprintf (fp2, "n"); void newval(double氺 aa)/*木函数将ih号换成新号*/int i, kl;for (i = 1; i <= n; i +) bi = 0.0;for (i = 1; i <= nb; i +)kl = nnewnobi; bkl = aai;for (i = 1; i <= n; i +) aai = bi

44、;void yzb(int t, int* iu, double* u, double* di, int* nfd)/*木函数求因子表 */参数1为标志（t=0求b，t=l求b ）/参数2因子表上三角矩阵非零非对角元素的列足码 /参数3因子表上三角矩阵非零非对角元素的数值 /参数4因子表上三角矩阵对角元素/参数5因子表上三角各行非零元索个数int i, j, k, il, int jj, jjl, jj2, double ai, bns; nfdl = 1; for (i = 1; i <=i2;im, x, fdns;n; i +)/nobt存放的是节点类型，0: pq节点，-1: p

45、v节点。 if (t != 1) | (nobti != -1) && i != mpj) / <|/ ifor (j = i + 1; j <= n; j +)/ ibj = 0.0;/ ibi = biii;/ iif (kk2 = 0) && (t = 1) && (nobti != -1)/ 存在(t = 1)的情况，不多余。bi = bi + afl * qli / v0i / v0i ;/af1 il = ydzi;i2 = ydzi + 1 - 1;for (j = il; j <= i2; j +)k = iyj

46、; bk = ybj;bmpj = 0. 0;if (t 1)for (j = 1; j <= n; j +)if (nobtj = -1)bj = 0. 0;il 二 i - 1;for (im = 1; im <= i 1; im +)jjl = nfdim;jj2 = nfdim +1 - 1;for (jj = jjl; jj <= jj2; jj +)if(iujj = i)ai 二 ujj / diim; for(k = jj; k 二 jj2; k +) j = iuk;bj = bj - ai * ukbreak;x 二 nfdi; dii二 1.0 / bi

47、; ai = di i; k = 0;il = i + 1;for (j = il; j <= n; j +)if (fabs(bj)1. oc-15)ux = bj * ai; iux二 j; k+; x+;fdi = k;elsefdi = 0; dii二 0.0;nfdi+l = nfdi + fdi;fprintf (fp2,"n氺氺氺氺氺氺氺氺u氺氺氺氺氺氺氺氺氺")；for (i = 1; i <= x; i +)if(i % 3 = 1)fprintf(fp2, "rt);fprintf (fp2, "%10.5f%8i&quo

48、t;, ui, iui);fprintf（fp2,"n氺氺氺氺氺氺氺氺di氺氺氺氺氺氺氺氺"）；printf 1 （di, n）;void printf1(double氺 aa, int n)/*本函数输出aai，t r no 1 f t t1 zzlnii<i=l一n本本本本/i +）if(i % 5 = 1)fprintf (fp2, "rt);fprintf (fp2, "%9. 5f", aai);fprintf(fp2, "nn");void calc(int* iu, double* u, double*

49、di, int* nfd, double* b)/*本函数利用因子表解线形方程组。（p417阁f1-9）*/doublebi;int i,j, k, il,i2;for (if=1; i <=n;i +)/前代过程。ibi=bi;il:nfdi;i2=nfdi +1；for(j = il；j <i2; j+)k = iuj;bk二 bk - bi * uj;bi = bi * di i;for (i = n; i >= 1; i ) / 冋代过程。 bi = bi;11 = nfdi;12 = nfdi + 1 - 1;for (j = i2; j >= il; j )

50、k = iuj;bi = bi - bk * uj;bi = bi;void zlsort (int* nnew) /*木函数进行支路数损排序 */小节点号放左边，大节点号放右边/左右皆按从小到大顺序排列/int ip, kl, k2, k3, k4; int i, j;for (i = 1; i <= zls; i +)k3 二 izl i; k4 = jzl i; kl = iabs (k3);k2 = iabs(k4) ;/ 鹿节点号。izli = isgn(nnewkl, k3) ;/ 新节点号。jzl i = isgn (nnewk2, k4);k3 = izl i; k4 =

51、 jzli; kl = iabs (k3); k2 = iabs (k4); if (kl > k2)izli = k4; jzli = k3;k2 = iabs (izl j);k4 = iabs (jzl j);if(k2 < kl | (k2 = kl && m < k3)ip = j； kl = k2; k3 = k4;if(i != ip)iswap(&izli, &izlip); iswap(&jzli, &jzlip); swap(&zri, &zrip); swap(&zxi, &

52、zxip); swap(&zyki, &zykip);void bnsopt ()/*节点优化 */int iil, ii2, zls2, nomax;int i, j, 1, kl, k;int temp;zls2 = zls + zls;for (i = 1; i = zls2; i +)oldi = nnewi = 0;/先清零。由此可知:ns4、ns必须大于2*zls。 for (i = 1; i <= zls; i +)oldi = iabs(izli); oldi + zls = iabs (jzli);/变压器节点号由正变负，old前zls个为左节点号，后zls个为右节点 号。for (i = 1; i <= zls2; i +)/ 冒泡法排序。kl = i + 1;for (j = kl; j <= zls2; j +) if (oldi > oldj)iswap(&oldi, &oldj);/交换整数oldi、oldjo小节点号排在支路左侧。nomax = old