电力系统计算报告

电力系记录算报告院（系）电气工程及自动化授课教师胡林献姓名张远实学号13S106052

P-Q分解潮流法简述P-Q分解潮流法基本原理、计算过程、计算框图。用C语言编制P-Q分解法潮流程序，并用电科院6机22节点系统加以验证。规定采用稀疏技术、因子表技术和节点优化技术，并考虑负荷静态特性。P-Q分解潮流法基本原理：P-Q分解法即是基于采用极坐标形式表达牛顿法，其依照电力系统实际运营状态线路参数R/X普通很小状况，对求解修正量修正方程系数矩阵加以简化，使其变为常数阵(即所谓等斜率)，且P、Q迭代解耦。这样可减少每次迭代计算时间，提高计算速度，又不影响最后成果，因而是普通选用一种办法。但在低电压配电网中，当线路R/X比值很大时，也许浮现不收敛状况，此时应考虑更换其他办法。计算过程：形成有功迭代和无功迭代简化雅克比矩阵和给定PQ节点初值和各节点电压相角初值作有功迭代，计算，解修正方程，得各节点电压相角修正值。作无功迭代，计算，解修正方程，得各节点电压幅值修正量。返回第三步，继续迭代到满足规定为止。计算框图流程：YY输入原始数据YY输入原始数据形成矩阵及并进行三角分解设PQ节点电压初值，各节点电压相角初值制迭代计数k=0用公式计算不平衡功率，计算置解修正方程求置用公式计算不平衡功率，计算置解修正方程求置计算平衡节点功率及所有线路功率输出YY算例描述用电科院6机22节点算例进行验证，详细参数和网络拓扑图如下所示表1线路参数支路号首末端节点号支路电阻支路电抗对地电纳/217-80.01060.07400.027-90.01470.10400.038-90.00340.01310.049-220.05590.21800.1954512-130.002450.02551.395614-190.00340.02000.0716-190.05780.21800.1887816-200.01630.06620.2353916-210.03740.17800.1641016-180.00330.03330.01119-210.01140.03700.01220-220.02140.08590.30081321-220.01500.06070.2198148-220.05370.19000.16531511-120.00330.03431.08797表2变压器支路数据支路号首末端节点号电阻电抗变比17-10.00.01501.05029-20.00.02171.075322-30.00.01241.100419-40.00.06401.025518-50.00.03751.050617-60.00.03371.000710-90.0-0.0021.000811-100.00.01801.000915-120.00.01801.0001017-130.00.01001.0001115-140.0-0.0021.0001216-170.00.00101.027表3并联电容器数据节点号电纳12-1.3665160.5017611-1.366512-1.3665表4母线功率数据母线号发电机输出有功发电机输出无功负荷有功负荷无功15.96311.73550.000.00263.20.000.0000.0000.0054.33.340.000.006-0.011.00.000.0070.000.000.000.0080.000.002.871.4490.000.003.762.21100.000.000.000.00110.000.000.000.00120.000.000.000.00130.000.000.000.00140.000.000.000.00150.000.000.000.00160.000.005.02.9170.000.000.000.00180.000.004.32.6190.000.000.8640.662200.000.000.7190.474210.000.000.70.5220.000.002.2651.69表5无功可调母线数据母线号电压幅值（标幺值）无功下限值无功上限制11.0-51031.0-5561.0-56表6发电机参数母线号暂态电抗转子惯性时间常数阻尼系数发电机有功出力下限发电机有功出力上限10.0150140.820.00316.520.038230.000.001.26.630.039679.500.001.58.2540.121015.680.0048039.200.001.025.6160.19762.620.000.21.1图1电科院6机22节点系统图计算成果IVCAPLQLPGQG11.000000.000000.000000.000005.963121.7354920.97384-11.3118930.000000.000006.000003.031.00000-27.4590570.000000.000003.100003.1473141.02190-25.1686100.000000.000001.600000.7000051.04392-28.2056560.000000.000004.300003.3400061.00000-37.5665660.000000.00000-0.010000.9166471.02697-5.2471700.000000.000000.000000.0000080.96847-19.788813-2.87000-1.440000.000000.0000090.98081-19.738201-3.76000-2.210000.000000.00000100.97985-19.3043440.000000.000000.000000.00000110.99050-23.1737710.000000.000000.000000.00000120.99359-30.4990340.000000.000000.000000.00000130.98198-35.5280900.000000.000000.000000.00000141.00020-30.6958610.000000.000000.000000.00000151.00103-30.7202820.000000.000000.000000.00000160.99308-27.750976-5.00000-2.900000.000000.00000170.96911-37.5466410.000000.000000.000000.00000180.98360-37.696751-4.30000-2.600000.000000.00000191.00776-31.018334-0.86400-0.662000.000000.00000201.01475-35.674534-0.71900-0.474000.000000.00000211.101546-32.104404-0.70000-0.500000.000000.00000221.05792-29.749725-2.26500-1.690000.000000.00000IJPIJQIJPJIQJI175.963121.73549-5.96312-1.15693296.000003.0-6.00000-2.141953223.100003.14731-3.10000-2.905314191.600000.70000-1.60000-0.513085184.300003.34000-4.30000-2.31986617-0.010000.916640.01000-0.88833783.481200.74380-3.353840.14533792.481910.41314-2.393680.2111289-0.28186-0.839220.284700.850168220.76571-0.74611-0.712140.595619103.63860-0.48476-3.638600.456749220.71038-0.64457-0.668940.3995110113.63860-0.45674-3.638600.7088611123.38055-2.04951-3.59392-1.1859912133.38055-1.10216-3.35200-1.3230612150.21337-0.40993-0.213370.4138313173.352001.32306-3.35200-1.188391415-0.213370.413390.21337-0.4138314190.21337-0.41339-0.212640.417721617-3.34200-2.060233.342002.0767116180.000270.28284-0.00000-0.280141619-0.51158-0.085670.52751-0.231991620-0.59323-0.399260.59959-0.049581621-0.55346-0.142840.56509-0.1326619210.42113-0.33465-0.417880.345192022-1.31859-0.424421.35499-0.075842122-0.84721-0.712540.861080.29604源程序#include"math.h"#include"stdio.h"#defineNS //最大节点数#defineNS2 NS*2#defineNS4 1000 //NS4、NS必要不不大于2*zls。#defineZS 3000 //最大支路数#defineZS2 ZS*2#defineDKS 200 //最大电抗器数#defineN2 ZS*4#defineN3 ZS*8+NS*4FILE*fp1，*fp2;charinname[12]，outname[12];// fp1输入数据文献指针 fp2输出文献指针// inname[]输入数据文献名 outname[]输出数据文献名intn，zls，nb，mdk，mpj，bnsopton，it1，dsd，kk2，nzls;// 节点总数n(涉及联系节点)支路数(回路数)zls节点数nb(发电机和负荷)// 接地电抗数mdk 精度eps平衡节点号mpj// 节点优化(标志)bnsopton(=0节点不优化，!=0节点优化)// 最大迭代次数it1 最低电压或最大功率误差节点号dsd// 负荷静特性标志(=0考虑负荷静特性)// 支路数(双回线算一条支路)intizl[ZS]，jzl[ZS]，idk[DKS]，yds[NS]，ydz[NS]，iy[ZS2];// izl[],jzl[],idk[]:分别存储左、右节点号和电抗器节点号。// yds[]存储各行非零非对角元素个数。// ydz[i]是第i行第一种非零非对角元素首地址，// 即在所有非零非对角元素中顺序号//iy[]存储列足码。intnnew[NS4]，old[NS]，nob[NS]，nobt[NS];// nnew[],old[]存储是新、旧节点号。// nnew[i]中为i相应新号// nob[]存储是节点号。nobt[]存储是节点类型，0：pq节点，-1：pv节点。doubleeps，dsm，vmin，dph，dqh，af[3];// eps迭代收敛精度，dsm最大功率误差// vmin:系统最低电压值。dph,dqh:系统有、无功损耗。// af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。doublev00; //v00：系统平均电压 ci,cj分别作为节点i,j电压相角暂时存储单元。doublezr[ZS]，zx[ZS]，zyk[ZS]，dkk[DKS]，gii[NS]，bii[NS]，yg[ZS2]，yb[ZS2];doublepg[NS]，qg[NS]，pl[NS]，ql[NS]，v0[NS]，v[NS]，va[NS];// 支路电阻zr[] 支路电抗zx[] 输电线路充电容纳zyk[](y0/2)// 接地电抗dkk[] 对角元实部gii[] 对角元虚部// 非对角元实部yg[] 非对角元虚部yb[]// pg[],qg[],pl[],ql[]:发电机，负荷功率实、虚部//v[]是电压幅值,va[]是电压相角。doublew[NS2]，kg[3]，b[NS2];intnewsort[NS4];// newsort[i]存储i相应老号voidinitial();voidpqflow();voidout();voiddataio();voidbnsopt();voidzlsort(int*nnew);voidprinto();voidprinty();voidy2();voidya0();voidyzb();voidjdgl(intkq0);voidbbhl(intkq0);voidcalc();intiabs(inta);voidbranch_output();voidnewval(double*aa);voidprintc(charaa，intn);//voidprintc();voidiswap(int*m，int*n);//voidiswap();//voidswap();voidswap(double*m，double*n);voidprintf2(double*aa，double*bb，intn);voidcalc(int*iu，double*u，double*di，int*nfd，double*b);voidprinti(int*aa，intn);voidprintf1(double*aa，intn);intfind(intk，inta[]，int*z);voidyzb(intt，int*iu，double*u，double*di，int*nfd);intisgn(inta，intb);voidyy1();voidy3();voidnewtoold();voidnode_output()；//intpq_main(void){ initial(); //初始化 pqflow(); //pq潮流计算 out(); //输出节点和支路数据 return1;}intisgn(inta，intb){//**** 本函数功能返回值为a绝对值b符号 ****////参数1提供值，参数2提供符号// if(b<0) if(a>0) a=-a; returna;}intfind(intk，inta[]，int*z){//**** 本函数查找a[]中与否有fabs(k)有则返回0，无则返回1 ****////参数1为待查找量，参数2待搜索数组，参数3返回k在a[]中顺序号// inti; for(i=1；i<=n；i++) if(iabs(k)==a[i]) { *z=i; return1; } return0;}voidoldtonew(){//**** 本函数将输入数据中节点号变成从1开始持续节点号 ****// inti，j，k，ii1，ii2，zls2，k1，k2，k3，k4，ip; zls2=zls+zls; for(i=1；i<=zls2；i++) newsort[i]=0; ii1=0; for(i=1；i<=zls；i++) { k=izl[i]; if(!find(k，newsort，&ii2)) { ii1++; newsort[ii1]=iabs(k); } k=jzl[i]; if(!find(k，newsort，&ii2)) { ii1++; newsort[ii1]=iabs(k); } } for(i=1；i<=ii1-1；i++) { for(j=i+1；j<=ii1；j++) { if(newsort[i]>newsort[j]) { k=newsort[i]; newsort[i]=newsort[j]; newsort[j]=k; } } } for(i=1；i<=zls；i++) { k=izl[i]; if(find(k，newsort，&ii2)) { izl[i]=isgn(ii2，k); } else printf("error!"); k=jzl[i]; if(find(k，newsort，&ii2)) { jzl[i]=isgn(ii2，k); } else printf("error!"); printf("izl[%d]=%d，jzl[%d]=%d\n"，i，izl[i]，i，jzl[i]); } for(i=1；i<=nb；i++) { for(j=1；j<=n；j++) if(nob[i]==newsort[j]) { nob[i]=j; break; } printf("nob[%d]=%d\n"，i，nob[i]); } for(j=1；j<=n；j++) { if(mpj==newsort[j]) { mpj=j; break; } } //电抗器节点号转变 for(j=1；j<=mdk；j++) { for(i=1；i<=n；i++) { if(idk[j]==newsort[i]) { idk[j]=i; break; } } }//hu-1-8 printc('-'，58); fprintf(fp2，"\n\n*******newsort*******\n"); for(i=1；i<=n；i++){ fprintf(fp2，"%8d%8d\n"，i，newsort[i]);} fprintf(fp2，"\n*******newsort*******\n");//hu-1-8}voidinitial(){//**** 本函数进行初始化工作 ****// inti，k1; dataio();//输入原始数据 oldtonew();//转化为新号 if(bnsopton==0)//节点不优化，新节点号即为老节点号。 for(i=1；i<=n；i++) { old[i]=i； nnew[i]=i; } else bnsopt(); //节点优化 mpj=nnew[mpj];//mpj:平衡节点 zlsort(nnew); //sortther,xandb for(i=1；i<=mdk；i++) { k1=idk[i]； idk[i]=nnew[k1]; } for(i=1；i<=n；i++) { v[i]=v00； va[i]=0.0; }//所有节点电压幅值初值都为1.000(v00)，电压相角初值都为0。 //exchangethenodebeforeandaftersort for(i=1；i<=n；i++) yds[i]=0；//theimmediate for(i=1；i<=nb；i++) { k1=nnew[nob[i]]； yds[k1]=nobt[i]； } for(i=1；i<=n；i++) nobt[i]=yds[i]; newval(pg); newval(qg); newval(pl); newval(ql); newval(v0); for(i=1；i<=n；i++)//nobt[]istypeofnode if(nobt[i]==-1) v[i]=v0[i]；//nob[]isserialsnumbe //nobt[]=-1：pv节点，v0[]存储是最后一种节点数据， //对于pv节点，即为该点应维持电压值。 //nobt[]=0：pq节点，v0[]存储是最后一种节点数据， //对于pq节点，即为系统平均电压值。 printo(); //输出af[]、v00和节点排序后支路、节点和 //接地电抗数据(仅仅查看中间成果) ya0();//获得yds[]、ydz[]、列足码iy[]。(P407)}voidprinto(){//**** 输出af[]、v00和节点排序后支路、节点和接地电抗数据 ****// inti; fprintf(fp2，"\n******AFANDV0******\n"); fprintf(fp2，"\n%7.3f%7.3f%7.3f\n"，af[0]，af[1]，v00); printc('-'，78); fprintf(fp2，"\n\n*******ZLB*******\n"); for(i=1；i<=zls；i++){ fprintf(fp2，"\n"); fprintf(fp2，"%8d%8d%8d%8d"，izl[i]，jzl[i]，old[abs(izl[i])]，old[abs(jzl[i])]); fprintf(fp2，"%9.4f%9.4f%9.4f"，zr[i]，zx[i]，zyk[i]);} printc('-'，78); fprintf(fp2，"\n\n*******BUS*******\n"); for(i=1；i<=nb；i++){ fprintf(fp2，"\n"); fprintf(fp2，"%8d%8d%8d"，nob[i],old[nob[i]]，nobt[i]); fprintf(fp2，"%9.4f%9.4f%9.4f%9.4f%9.4f"，pg[i]，qg[i]，pl[i]，ql[i]，v0[i]);} printc('-'，78); fprintf(fp2,"\n\n******DKK******\n"); for(i=1；i<=mdk；i++){ fprintf(fp2，"\n"); fprintf(fp2，"%8d%8d%7.4f"，idk[i]，old[idk[i]]，dkk[i]);}}voiddataio(){//**** 系统数据初始化 ****// inti; af[0]=0.6； af[1]=2.0;//af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。 v00=1.000;//系统平均电压 printf("\npleaseinputthenameofdatafile\n"); scanf("%s"，inname); fp1=fopen(inname，"r"); printf("\npleaseoutputthenameofdatafile\n"); scanf("%s"，outname); fp2=fopen(outname，"w"); fscanf(fp1，"%d%d%d%d"，&n，&zls，&nb，&mdk)； //thenumberofnode,branches，node fscanf(fp1，"%lf%d%d%d%d"，&eps，&kk2，&mpj， &bnsopton，&it1)； //precision，swingnode,sortthenode,iterationnumbers for(i=1；i<=zls；i++){ fscanf(fp1，"%d%d"，&izl[i]，&jzl[i]); fscanf(fp1，"%lf%lf%lf"，&zr[i]，&zx[i]，&zyk[i]);} for(i=1；i<=nb；i++) { fscanf(fp1，"%d%d"，&nob[i]，&nobt[i]); fscanf(fp1，"%lf%lf%lf%lf%lf"，&pg[i]，&qg[i]，&pl[i]， &ql[i]，&v0[i]); } for(i=1；i<=mdk；i++) { fscanf(fp1，"%d%lf"，&idk[i]，&dkk[i])； } fclose(fp1);}voidpqflow() {//**** PQ分解法计算潮流 ****// intkq0，iu1[N2]，nfd1[NS]，iu2[N2]，nfd2[NS]; inti，t； doubleu1[N2]，u2[N2]，di1[NS]，di2[NS]; yy1()； yzb(0，iu1，u1，di1，nfd1)；//formtheBmatrixofP-0iteration y2(); yzb(1，iu2，u2，di2，nfd2)；//formtheBmatrixofQ-Viteration t=0； kq0=0; kg[0]=kg[1]=1; do { jdgl(kq0)；//calculatingthepower bbhl(kq0)；//findoutthemaxi if(kq0==0) printf("P：%d\t%d\t%f\n"，t，dsd，dsm); else printf("Q：%d\t%d\t%f\n"，t，dsd，dsm); if(fabs(dsm)>eps) { kg[kq0]=1; if(kq0==0) calc(iu1，u1，di1，nfd1，b); if(kq0==1) calc(iu2，u2，di2，nfd2，b); for(i=1；i<=n；i++) { if(kq0==0) va[i]=va[i]-b[i]/v00; else v[i]=v[i]-b[i]; } } else kg[kq0]=0; if(kq0==0) kq0=1; else { kq0=0； t++; } if(t>it1) break; }while((fabs(dsm)>eps)||(kg[kq0]!=0)); fprintf(fp2，"\n%s%d"，"times="，t);}voidout(){//**** 本函数输出节点和支路数据 ****// zlsort(old)；//recoverthedataifsorted// newtoold(); node_output()；//nodedata branch_output()；//branchdata printc('-'，78); printc('*'，78); fprintf(fp2，"\n");}voidnewval(double*aa){//**** 本函数将旧号换成新号 ****// inti，k1; for(i=1；i<=n；i++) b[i]=0.0; for(i=1；i<=nb；i++) { k1=nnew[nob[i]]； b[k1]=aa[i]； } for(i=1；i<=n；i++) aa[i]=b[i];}voidyzb(intt，int*iu，double*u，double*di，int*nfd){//**** 本函数求因子表 ****////参数1为标志(t=0求B',t=1求B'')////参数2因子表上三角矩阵非零非对角元素列足码//参数3因子表上三角矩阵非零非对角元素数值//参数4因子表上三角矩阵对角元素//参数5因子表上三角各行非零元素个数 inti，j，k，i1，i2； intjj，jj1，jj2，im，x，fd[NS]; doubleai，b[NS]; nfd[1]=1; for(i=1；i<=NS；i++)//hu-1-9 { u[i]=0.0; }//hu-1-9 for(i=1；i<=n；i++) { //nobt[]存储是节点类型，0：pq节点，-1：pv节点。 if(((t!=1)||(nobt[i]!=-1))&&i!=mpj)//<---| {//| for(j=i+1；j<=n；j++)// | b[j]=0.0；//| b[i]=bii[i]；//| if((kk2==0)&&(t==1)&&(nobt[i]!=-1))//存在(t==1)状况，不多余。 b[i]=b[i]+af[1]*ql[i]/v0[i]/v0[i];//af[1] i1=ydz[i]; i2=ydz[i+1]-1; for(j=i1；j<=i2；j++) { k=iy[j]； b[k]=yb[j]; } b[mpj]=0.0; if(t==1) for(j=1；j<=n；j++) if(nobt[j]==-1) b[j]=0.0; i1=i-1; for(im=1；im<=i1；im++) { jj1=nfd[im]; jj2=nfd[im+1]-1; for(jj=jj1；jj<=jj2；jj++) { if(iu[jj]==i) { ai=u[jj]/di[im]; for(k=jj；k<=jj2；k++) { j=iu[k]; b[j]=b[j]-ai*u[k]; } break; } } } x=nfd[i]; di[i]=1.0/b[i]; ai=di[i]; k=0; i1=i+1; for(j=i1；j<=n；j++) { if(fabs(b[j])>1.0e-15) { u[x]=b[j]*ai; iu[x]=j; k++; x++; } } fd[i]=k; } else { fd[i]=0； di[i]=0.0; } nfd[i+1]=nfd[i]+fd[i]; } fprintf(fp2，"\n********U*********"); for(i=1；i<x；i++) { if(i%3==1) fprintf(fp2，"\n"); fprintf(fp2，"%9.5f%7d"，u[i]，iu[i]); }// fprintf(fp2，"\n***********xxxx**************");// fprintf(fp2，"%8d%10.5f%8d",x,u[38],fd[22]);fprintf(fp2，"\n***********fd**************"); printi(fd，n);fprintf(fp2，"\n***********nfd**************"); printi(nfd，n); fprintf(fp2，"\n********DI********"); printf1(di，n);}voidprintf1(double*aa，intn){//**** 本函数输出aa[i],i=1-n ****// inti; for(i=1；i<=n；i++){ if(i%5==1) fprintf(fp2，"\n"); fprintf(fp2，"%9.5f"，aa[i]);} fprintf(fp2，"\n\n");}voidcalc(int*iu，double*u，double*di，int*nfd，double*b){//**** 本函数运用因子表解线形方程组。(P417图F1-9) ****// doublebi; inti，j，k，i1，i2; for(i=1；i<=n；i++)//前代过程。{ bi=b[i]; i1=nfd[i]; i2=nfd[i+1]; for(j=i1；j<i2；j++) { k=iu[j]； b[k]=b[k]-bi*u[j]； } b[i]=bi*di[i]; } for(i=n；i>=1；i--)//回代过程。{ bi=b[i]; i1=nfd[i]; i2=nfd[i+1]-1; for(j=i2；j>=i1；j--) { k=iu[j]； bi=bi-b[k]*u[j]; } b[i]=bi;}}voidzlsort(int*nnew) //// 。{//**** 本函数进行支路数据排序 ****////小节点号放左边，大节点号放右边////左右皆按从小到大顺序排列// intip，k1，k2，k3，k4; inti，j; for(i=1；i<=zls；i++){ k3=izl[i]; k4=jzl[i]; k1=iabs(k3); k2=iabs(k4); // 原节点号。 izl[i]=isgn(nnew[k1]，k3)；//新节点号。 jzl[i]=isgn(nnew[k2]，k4); k3=izl[i]; k4=jzl[i]; k1=iabs(k3); k2=iabs(k4); if(k1>k2) { izl[i]=k4； jzl[i]=k3; }} for(i=1；i<=zls-1；i++){ ip=i; k1=iabs(izl[i]); k3=iabs(jzl[i]); for(j=i+1；j<=zls；j++) { k2=iabs(izl[j]); k4=iabs(jzl[j]); if(k2<k1||(k2==k1&&k4<k3)) { ip=j； k1=k2； k3=k4; } } if(i!=ip) { iswap(&izl[i]，&izl[ip]); iswap(&jzl[i]，&jzl[ip]); swap(&zr[i]，&zr[ip]); swap(&zx[i]，&zx[ip]); swap(&zyk[i]，&zyk[ip]); }}}voidbnsopt(){//**** 节点优化 ****// intii1，ii2，zls2，nomax; inti，j，l，k1，k； inttemp; zls2=zls+zls; for(i=1；i<=zls2；i++) old[i]=nnew[i]=0;//先清零。由此可知：NS4、NS必要不不大于2*zls。 for(i=1；i<=zls；i++) { old[i]=iabs(izl[i])； old[i+zls]=iabs(jzl[i]); } //变压器节点号由正变负，old[]前zls个为左节点号，后zls个为右节点号。 for(i=1；i<=zls2；i++)// 冒泡法排序。{ k1=i+1; for(j=k1；j<=zls2；j++) if(old[i]>old[j]) iswap(&old[i]，&old[j]); //互换整数old[i]、old[j]。小节点号排在支路左侧。} nomax=old[zls2];//nomax即是最大节点号。Iee30.dat----30 l=1; for(i=1；i<=n；i++){ ii1=old[l]; old[i]=ii1; for(j=l；j<=zls2；j++) { ii2=old[j]; if(ii1!=ii2) { l=j； break; } nnew[i]++; } } for(i=1；i<=n-1；i++){ for(j=i+1；j<=n；j++) if(nnew[j]<=nnew[i]) if((nnew[j]!=nnew[i])||(old[j]<old[i])) { iswap(&old[i]，&old[j]); iswap(&nnew[i]，&nnew[j]); } } for(i=1；i<=nomax；i++) nnew[i]=0; for(i=1；i<=n；i++) { j=old[i]； nnew[j]=i; }//hu-1-8 printc('-'，68); fprintf(fp2，"\n\n*******inewold*******\n"); for(i=1；i<=n；i++){ fprintf(fp2，"%8d%8d%8d\n"，i，nnew[i],old[i]);} fprintf(fp2，"\n*******inewold*******\n");//hu-1-8}voidya0()//(P407){//**** 本函数获得yds[]、ydz[]、列足码iy[] ****// inti，j，l，ll; for(i=1；i<=n；i++) yds[i]=0；//yds[]存储各行非零非对角元素个数。 ll=1; for(l=1；l<=zls；l++){ i=iabs(izl[l]); j=iabs(jzl[l]); if(i==j) continue; iy[ll]=j; if((i!=iabs(izl[l+1]))||(j!=iabs(jzl[l+1]))) { ll++； //ll记录总支路数(双回线算一条支路) yds[i]++; }} nzls=ll-1；//总支路数(双回线算一条支路) ydz[1]=1; for(i=1；i<=n-1；i++) ydz[i+1]=yds[i]+ydz[i]；//由yds[]得ydz[]。//ydz[i]是第i行第一种非零非对角元素首地址，// 即在所有非零非对角元素中顺序号。 fprintf(fp2，"\n*******YDZ********")； printi(ydz，n); fprintf(fp2，"\n"); fprintf(fp2，"\n*******YDS********"); printi(yds，n); fprintf(fp2，"\n\n");}voidprinti(int*aa，intn){//**** 本函数输出aa[1]-aa[n]值 ****// inti; for(i=1；i<=n；i++){ if(i%10==1) fprintf(fp2，"\n"); fprintf(fp2，"%5d"，aa[i]);}}voidy2() {//**** 本函数形成节点导纳阵，一次形成，不分B'和B'' ****// intj1; doubler，x，yk，zf，gij，bij; inti，j，i1，l，ll； for(i=1；i<=n；i++) { gii[i]=0.0； bii[i]=0.0; }//导纳阵对角元(与节点一一相应)先清零。 for(i=1；i<=mdk；i++) { j=idk[i]； bii[j]=-1.0/dkk[i]; }//计算接地支路导纳，只影响导纳阵对角元(自导纳)。 for(i=1；i<=zls；i++) { yg[i]=0.0； yb[i]=0.0; }//导纳阵非零非对角元(与支路一一相应)先清零。 ll=1; for(l=1；l<=zls；l++) { i1=izl[l]; // 支路左节点号。 j1=jzl[l]; // 支路右节点号。 i=iabs(i1); // 变压器支路有一节点号为负值。 j=iabs(j1); if(i==j) continue；//排除左、右节点号相等状况。 r=zr[l]; x=zx[l]; yk=zyk[l]; // zr[],zx[],zyk[]:支路三参数。 zf=r*r+x*x; gij=r/zf; bij=-x/zf; if((i1>0)&&(j1>0))// 不是变压器支路。是普通支路。 { yg[ll]=yg[ll]-gij; yb[ll]=yb[ll]-bij; gii[i]=gii[i]+gij; bii[i]=bii[i]+bij+yk; gii[j]=gii[j]+gij; bii[j]=bii[j]+bij+yk； } else // 变压器支路。 { if(j1<0) { i=iabs(j1)； j=iabs(i1); } // 若非原则变比在右(j)侧，则左、右互换，保证非原则变比在左侧。 gii[j]=gii[j]+gij; bii[j]=bii[j]+bij; // 原则变比侧。 gii[i]=gii[i]+gij/yk/yk； // yk=zyk[]，对变压器支路指非原则变比(设在节点号为负一侧)。 bii[i]=bii[i]+bij/yk/yk； // 非原则变比侧阻抗计算。 yg[ll]=yg[ll]-gij/yk; yb[ll]=yb[ll]-bij/yk; } //◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆ //注意此处解决双回线和导纳阵一不同样 if((i1!=izl[l+1])||(j1!=jzl[l+1])) ll++; }// 打印导纳矩阵。对角元实部为gii，虚部为bii，// 非零非对角元实部为yb[]，虚部为yb[]，列足码为iy[]。 fprintf(fp2，"*******GII,BII********"); printf2(gii,bii,n); fprintf(fp2，"\n*******YYYYY********"); for(i=1；i<=nzls；i++){ if(i%2==1) fprintf(fp2，"\n"); fprintf(fp2，"%10.4f%10.4f%8d"，yg[i]，yb[i]，iy[i]);}}voidprintf2(double*aa，double*bb，intn){//**** 本函数打印aa[i],bb[i]i=1-n ****// inti; for(i=1；i<=n；i++){ if(i%2==1) fprintf(fp2，"\n"); fprintf(fp2，"%5d%10.4f%10.4f"，i，aa[i]，bb[i]);}}voidjdgl(intkq0){//**** 本函数计算节点功率(流程图见p421) ****// doubleai，bi；//ai,bi是暂时工作单元。 doublevi，vj，ci，sn，cs; inti，j，k，i1，i2； for(i=1；i<=n；i++) w[i]=0.0;//w[]存储节点功率,一方面清零。 for(i=1；i<=n；i++){ vi=v[i]；//v[]存储节点电压幅值。 ai=-bii[i]；//bii[]存储导纳阵对角元虚部(gii[]+j*bii[]) if(kq0==0) ai=gii[i]；//gii[]存储导纳阵对角元实部(gii[]+j*bii[]) w[i]=w[i]+vi*vi*ai； //对角元素节点功率p=vi*vi*gii,q=vi*vi*bii if(i<n)//导纳阵最后一行没有非对角元，故i<n而不能i=n. { i1=ydz[i]； //ydz[i]是第i行第一种元素首地址， //既在所有非零非对角元素中顺序号。 i2=ydz[i+1]-1；//即为第i行非零元素个数。 for(k=i1；k<=i2；k++)//对第i行所有非零元素进行操作。 { if(kq0!=0) { ai=-yb[k]； bi=yg[k]; }//yb[]存储导纳阵非对角元虚部(yg[]+j*yb[]) else { ai=yg[k]； bi=yb[k]; }//yg[]存储导纳阵非对角元实部(yg[]+j*yb[]) j=iy[k]; //iy[]存储是列足码。 vj=vi*v[j]；//v[]存储节点电压幅值。 ci=va[i]-va[j]; //va[]存储节点电压相角，ci即为节点i和 //节点j之间相角差(i-j)。 sn=sin(ci)；//sin()近似计算式 cs=cos(ci); //cos()近似计算式 bi=bi*vj*sn; ai=ai*vj*cs; w[i]=w[i]+ai+bi；//非对角元素i功率 w[j]=w[j]+ai-bi；//非对角元素j功率 } }}}voidbbhl(intkq0){//**** 本函数计算各节点功率误差，求最大功率误差dsm ****////**** 和常数项b[i]。(程序框图见P423) ****// inti; doublevi，vj; doublepl0，pg0; dsm=0.0；//dsm即为最大功率误差。 for(i=1；i<=n；i++) { vi=v[i]；//v[]存储节点电压幅值。 vj=v0[i]; //v0[]存储节点初始电压幅值。v0[]存储是最后一种节点数据。 //vi[] //nopt[]=-1：pv节点，对于pv节点，即为该点应维持电压值。 //nopt[]=0：pq节点，对于pq节点，即为系统平均电压值。 //vj此时表达是节点正常电压会定值。 if(kq0==0) { pl0=pl[i]; pg0=pg[i]； }//负荷p,发电机p else{ pl0=ql[i]； pg0=qg[i]； }//负荷q,发电机q if(kk2==0) pl0=pl0*((vi-vj)*af[kq0]/vj+1.0)； //考虑负荷静特性。(公式F1-39,P423) if(nobt[i]==-1&&kq0==1) qg[i]=w[i]-pl0；//pv节点 if(i==mpj&&kq0==0) pg[i]=w[i]-pl0；//平衡节点 b[i]=pg0+pl0-w[i]; //pv节点(nobt[]=-1)和平衡节点(mpj)不参加 //求最大功率误差和常数项运算 if(((kq0!=1)||(nobt[i]!=-1))&&(i!=mpj)) { if(fabs(b[i])>fabs(dsm)) { dsm=b[i]； dsd=i; }//dsm即为最大功率误差,dsd存储其相应节点号。 b[i]=b[i]/vi；//计算修正方程式常数项。 } else b[i]=0.0; //pv节点(nobt[]=-1)和平衡节点(mpj)不参加 //求最大功率误差和常数项运算。 }}voidnode_output()//{//**** 输出节点数据和最小电压幅值、相角(角度)及其节点号 ****////**** (程序框图见p426F1-16) ****// doublevi，ci; inti，j，oldnumber; printc('+'，72); fprintf(fp2，"\n%5s%8s%10s%11s%11s%11s%11s\n"，"I"， "V"，"CA"，"PL"，"QL"，"PG"，"QG"); vmin=v[1]; dsd=1; for(i=1；i<=n；i++) { j=nnew[i]; oldnumber=newsort[i];//转化为相应旧号 ci=va[j]*180.0/3.1416；//弧度转化为角度。 vi=v[j]; if(vi<vmin) { vmin=vi； dsd=j; }//vmin即为最小电压,dsd存储其相应新节点号。 fprintf(fp2，"\n%5d%11.5f%12.6f"，oldnumber，vi，ci); fprintf(fp2，"%11.5f%11.5f%11.5f%11.5f"，pl[j]，ql[j], pg[j]，qg[j]); } printc('-'，72);}voidbranch_output(){//**** 本函数输出支路数据。(程序框图见p428F1-17) ****// intii，jj; doubler，x，yk，zf，vi，vj，ci，cj； inti，j，l; doublede，df，ei，ej，fi，fj，fii，fir，pij，pji，qij，qji; dph=0.0；//记录系统有功网损。 dqh=0.0；//记录系统无功网损。 fprintf(fp2，"\n%5s%5s%10s%12s%12s%12s\n"，"I"，"J"，"PIJ", "QIJ"，"PJI"，"QJI"); for(i=1；i<=mdk；i++){ j=idk[i]; dkk[i]=v[j]*v[j]/dkk[i];} for(l=1；l<=zls；l++){ ii=iabs(izl[l])；//izl[]：支路左节点号。 jj=iabs(jzl[l])；//jzl[]：支路右节点号。 i=nnew[ii]; j=nnew[jj]；//转换为新节点号。 ii=newsort[ii]; jj=newsort[jj];//转化为相应旧号 r=zr[l]; x=zx[l]; yk=zyk[l]; vi=v[i]；//v[]：电压幅值。 ci=va[i]；//va[]：电压相角。 vj=v[j]; cj=va[j]; //支路左、右节点电压值由极坐标转换为直角坐标 ei=vi*cos(ci)； fi=vi*sin(ci)； //ei：支路左节点电压实部，fi：支路左节点电压虚部。 ej=vj*cos(cj)； fj=vj*sin(cj)； //ej：支路右节点电压实部，fj：支路右节点电压虚部。 if((izl[l]<0)||(jzl[l]<0))//变压器支路。 { if(izl[l]<0) { ei=ei/yk； fi=fi/yk; }//yk=zyk[l] else { ej=ej/yk； fj=fj/yk; } yk=0.0; } de=ei-ej; df=fi-fj; zf=r*r+x*x; fii=(de*r+df*x)/zf; fir=(df*r-de*x)/zf; pij=fii*ei+fir*fi; qij=fii*fi-fir*ei; pji=-fii*ej-fir*fj; qji=-fii*fj+fir*ej; qij=qij-vi*vi*yk; qji=qji-vj*vj*yk; dph=dph+pij+pji; dqh=dqh+qij+qji; fprintf(fp2，"\n%5d%5d%12.5f%12.5f%12.5f%12.5f"，ii，jj，pij, qij，pji，qji); } fprintf(fp2，"\n\n%8s%19s%18s%13s"，"DPH"，"DQH"，"VMIN"，"DSD"); fprintf(fp2，"\n%10.5f%18.5f%18.5f%12d"，dph，dqh，vmin，newsort[old[dsd]]); fprintf(fp2，"\n\n%22sTHEENDOFOUTPUT"，"");}voidiswap(int*m，int*n){//**** 本函数互换m和n值 ****// intp; p=*m； *m=*n； *n=p;}voidswap(double*m，double*n){//**** 本函数互换m和n值 ****// doublep; p=*m； *m=*n； *n=p;}intiabs(inta){//**** 本函数求绝对值 ****// a=(a>=0)?a：-a; returna;}voidprintc(charaa，intn){//**** 本函数输出n个aa字符 ****// inti; fprintf(fp2，"\n"); for(i=1；i<=n；i++) fprintf(fp2，"%c"，aa);}voidnewtoold(){ inti，j，k，ii1，ii2，zls2，k1，k2，k3，k4，ip; zls2=zls+zls; ii1=0; for(i=1；i<=zls；i++) { k=izl[i]; ii2=newsort[(int)fabs(k)]; izl[i]=isgn(ii2，izl[i]); k=jzl[i]; ii2=newsort[(int)fabs(k)]; jzl[i]=isgn(ii2，jzl[i]); printf("izl[%d]=%d，jzl[%d]=%d\n"，i，izl[i]，i，jzl[i]); } for(i=1；i<=zls-1；i++){ ip=i; k1=iabs(izl[i]); k3=iabs(jzl[i]); for(j=i+1；j<=zls；j++) { k2=iabs(izl[j]); k4=iabs(jzl[j]); if(k2<k1||(k2==k1&&k4<k3)) { i