C语言-潮流计算实现

word文档可自由复制编辑C语言程序设计潮流计算学院自动化学院专业班级学号姓名联系方式本程序潮流计算部分采用牛顿拉夫逊极坐标法进行计算，求解一次多元方程采用高斯列主元分解法进行求解。根据工程实际，在存储文件实时记录产生文件时间。此外本程序特意增加了文件查看功能，方便文件的查看。程序代码：#include<stdio.h>#include<malloc.h>#include<stdlib.h>#include<math.h>#include<string.h>#include<time.h>#include<conio.h>voiddataprepare(void);voidinitial(void);voidYc(void);voidshowy(void);doubledetpqc(void);voidshowdetav(void);voidJrc(void);voidshowdetpq(void);voidGauss(void);voidshowj(void);voidshowsolution(void);charsave2file(void);charoption(void);intchose(void);intlist(charfilename[],inti,chardstn[]);structPQV{charname[5];doublevb;doublep;doubleq;doublev;doubleag;}*pqv;structY{ doubleG; doubleB;}**y;structLZ{ charname1[5]; doublevb1; charname2[5]; doublevb2; doubler; doublex; doubleb;}*lz;structTZ{ charname1[5]; doublevb1; charname2[5]; doublevb2; doublex; doublek1; doublek2;}*tz;charfsource[20];intpqnum=0,pvnum=0,lznum=0,tznum=0;inttemp1=0,temp2=0,temp3=0,temp4=0;intdiedai=0;intnumw=1,numofw=0;double*detpq,*detav,*Pi,*Qi,error,wucha,**J;/*detpq节点功率的误差量，detav修正量，Pi节点的有功功率*/charc,ch,dstname[20],filename[20];intmain(){opt:c=option(); if(c=='1')gotoopt1; elseif(c=='2') { list("1.bin",1,filename); system("modecon:cols=100lines=30&color07"); while(c!='0') { system("cls"); printf("################################计算结果文件查看################################\n"); printf("[8]上一个[2]下一个[5]查看[-]删除文件【0】退出\n\n"); list("1.bin",0,filename); chose(); } gotoopt; } elseif(c=='0')exit(0); elsegotoopt;opt1:dataprepare();init:initial(); Yc(); printf("输入允许误差量："); scanf("%lf",&wucha); printf("选择模式：a.详情；b.简单：");redo:fflush(stdin); c=getchar(); if(c!='a'&&c!='b'&&c!='A'&&c!='B') {printf("输入错误！重新输入：");gotoredo;} if(c!='B'&&c!='b')showy(); error=detpqc(); showdetpq(); while(error>wucha&&error<100) { if(c!='B'&&c!='b') { printf("\t\t是否继续迭代<空格：是R/r:重赋初值E/e:退出潮流计算>？\n"); ifdie: fflush(stdin); ch=getch(); switch(ch) { case'':break; case'R':; case'r':gotoinit;break; case'E':; case'e':gotoopt;break; default:gotoifdie; } } diedai++; Jrc(); if(c!='B'&&c!='b') { printf("\n>>>>>>>>>>>>>>>>>>>第%d次迭代<<<<<<<<<<<<<<<<<<<<<\n",diedai); showj(); } Gauss(); if(c!='B'&&c!='b')showdetav(); for(temp1=0;temp1<pqnum+pvnum;temp1++) /*修正电压*/ pqv[temp1].ag=pqv[temp1].ag-detav[temp1]; for(temp1=0;temp1<pqnum;temp1++) pqv[temp1].v=pqv[temp1].v-detav[temp1+pqnum+pvnum]; error=detpqc(); if(c!='B'&&c!='b') showdetpq(); } if(error>=300||error<0) { printf("\t\t迭代不收敛<R/r:重赋初值E/e:退出潮流计算>？\n"); what: fflush(stdin); ch=getch(); switch(ch) { case'R':; case'r':gotoinit;break; case'E':; case'e':gotoopt;break; default:gotowhat; } } showsolution(); ch=save2file(); if(ch=='1')gotoopt; else gotoinit;}voiddataprepare(void){ FILE*fp; charch,type[4]; pqnum=0,pvnum=0,lznum=0,tznum=0,temp1=0,temp2=0,temp3=0,temp4=0; printf("\n 请输入潮流计算的数据文件\n --");getf:fflush(stdin); gets(fsource); while((fp=fopen(fsource,"r"))==NULL) /*打开文件*/ { printf("\t\t不能打开文件！！\n\t\t重新输入请键入“Y”，其他任意键退出:"); fflush(stdin); ch=getchar(); if(ch=='Y'||ch=='y') { printf("\n\t\t潮流计算的数据文件："); gotogetf; } else exit(0); } system("modecon:cols=130lines=50&color06"); while(!feof(fp)) /*计算各种节点数目*/ { ch=fgetc(fp); if(ch==10) { if(fgets(type,3,fp)) { if(strcmp(type,"pq")==0)pqnum++; if(strcmp(type,"pv")==0)pvnum++; if(strcmp(type,"lz")==0)lznum++; if(strcmp(type,"tz")==0)tznum++; } } } pqv=(structPQV*)malloc((pqnum+pvnum+1)*sizeof(structPQV)); /*根据节点数目开辟各节点储存空间*/ lz=(structLZ*)malloc((lznum)*sizeof(structLZ)); tz=(structTZ*)malloc((tznum)*sizeof(structTZ)); rewind(fp); while(!feof(fp)) /*读取和输入各节点数据*/ { ch=fgetc(fp); if(ch==10) { if(!fgets(type,3,fp))break; if(strcmp(type,"pq")==0) { fscanf(fp,"%s%lf%lf%lf",pqv[temp1].name,&pqv[temp1].vb,&pqv[temp1].p,&pqv[temp1].q); temp1++; } if(strcmp(type,"pv")==0) { fscanf(fp,"%s%lf%lf\t%lf",pqv[temp2+pqnum].name,&pqv[temp2+pqnum].vb,&pqv[temp2+pqnum].p,&pqv[temp2+pqnum].v); temp2++; } if(strcmp(type,"ph")==0) { fscanf(fp,"%s%lf\t\t%lf%lf",pqv[pqnum+pvnum].name,&pqv[pqnum+pvnum].vb,&pqv[pqnum+pvnum].v,&(pqv[pqnum+pvnum].ag)); pqv[pqnum+pvnum].ag=pqv[pqnum+pvnum].ag/180*3.1415926; } if(strcmp(type,"lz")==0) { fscanf(fp,"%s%lf%s%lf%lf%lf%lf",lz[temp3].name1,&lz[temp3].vb1,lz[temp3].name2,&lz[temp3].vb2,&lz[temp3].r,&lz[temp3].x,&lz[temp3].b); temp3++; } if(strcmp(type,"tz")==0) { fscanf(fp,"%s%lf%s%lf%lf%lf%lf",tz[temp4].name1,&tz[temp4].vb1,tz[temp4].name2,&tz[temp4].vb2,&tz[temp4].x,&tz[temp4].k1,&tz[temp4].k2); temp4++; } } } fclose(fp); puts("***************************************************************************************************************************"); puts(" 电力系统潮流计算 "); puts("***************************************************************************************************************************"); printf("从文件中获得的数据如下：\n"); printf("节点类型名称额定电压有功标幺无功标幺电压标幺相角弧度\n"); for(temp1=0;temp1<=pqnum+pvnum;temp1++) {if(temp1<pqnum) printf("PQ%s%9.3lf%9.3lf%9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[temp1].p,pqv[temp1].q); elseif(temp1<pqnum+pvnum) printf("PV%s%9.3lf%9.3lf%9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[temp1].p,pqv[temp1].v); else printf("平衡%s%9.3lf%9.3lf%9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[temp1].v,pqv[temp1].ag); } printf("\n"); printf("线路节点1额定电压节点2额定电压电阻电抗对地导纳\n"); for(temp1=0;temp1<lznum;temp1++) printf("输电线路%s%9.3lf%s%9.3lf%9.3lf%9.3lf%9.5lf\n",lz[temp1].name1,lz[temp1].vb1,lz[temp1].name2,lz[temp1].vb2,lz[temp1].r,lz[temp1].x,lz[temp1].b); printf("\n"); printf("线路节点1额定电压节点2额定电压电抗值一次比二次\n"); for(temp1=0;temp1<tznum;temp1++) printf("变压器%s%9.3lf%s%9.3lf%9.3lf%9.3lf：%9.3lf\n",tz[temp1].name1,tz[temp1].vb1,tz[temp1].name2,tz[temp1].vb2,tz[temp1].x,tz[temp1].k1,tz[temp1].k2);}voidinitial(void){ charc;re_insert: for(temp1=0;temp1<pqnum+pvnum;temp1++) { if(temp1<pqnum) { printf("请输入第%d个节点%s(pq)的电压幅值：",temp1+1,pqv[temp1].name); fflush(stdin); scanf("%lf",&pqv[temp1].v); printf("请输入第%d个节点%s(pq)的电压相角：",temp1+1,pqv[temp1].name); fflush(stdin); scanf("%lf",&pqv[temp1].ag); pqv[temp1].ag=pqv[temp1].ag/180*3.1415926; } if(temp1>=pqnum) { printf("请输入第%d个节点%s(pv)的电压相角：",temp1+1,pqv[temp1].name); fflush(stdin); scanf("%lf",&pqv[temp1].ag); } } printf("\t\t\t按ESC重新赋值，其他任意键继续！\n"); c=getch(); if(c==27)gotore_insert;}voidYc(void){ inti,j,temp,flag=0; y=(structY**)malloc((pqnum+pvnum+1)*sizeof(structY*)); for(temp1=0;temp1<pqnum+pvnum+1;temp1++) y[temp1]=(structY*)malloc((pqnum+pvnum+1)*sizeof(structY)); for(i=0;i<=pqnum+pvnum;i++) { for(j=0;j<=(pqnum+pvnum);j++) { if(i!=j) { flag=0; for(temp=0;temp<lznum;temp++)/*当两节点间为线路时*/ if((strcmp(pqv[i].name,lz[temp].name1)==0&&strcmp(pqv[j].name,lz[temp].name2)==0)||(strcmp(pqv[i].name,lz[temp].name2)==0&&strcmp(pqv[j].name,lz[temp].name1)==0)) { flag++; y[i][j].G=-lz[temp].r/(pow(lz[temp].r,2)+pow(lz[temp].x,2)); y[i][j].B=lz[temp].x/(pow(lz[temp].r,2)+pow(lz[temp].x,2)); } for(temp=0;temp<tznum;temp++)/*当两节点间为变压器*/ if((strcmp(pqv[i].name,tz[temp].name1)==0&&strcmp(pqv[j].name,tz[temp].name2)==0)||(strcmp(pqv[i].name,tz[temp].name2)==0&&strcmp(pqv[j].name,tz[temp].name1)==0)) { flag++; y[i][j].G=0; y[i][j].B=1/tz[temp].x*tz[temp].k1/tz[temp].k2; } if(flag==0){y[i][j].G=0;y[i][j].B=0;} if(flag>1)printf("\n两节点间出现多条支路，本程序暂不支持!!\n"); } } } for(i=0;i<=pqnum+pvnum;i++) { y[i][i].G=0,y[i][i].B=0; for(temp=0;temp<lznum;temp++)/*当与另一节点间为线路时*/ if((strcmp(pqv[i].name,lz[temp].name1)==0)||(strcmp(pqv[i].name,lz[temp].name2)==0)) { y[i][i].G+=lz[temp].r/(pow(lz[temp].r,2)+pow(lz[temp].x,2)); y[i][i].B+=-lz[temp].x/(pow(lz[temp].r,2)+pow(lz[temp].x,2))+lz[temp].b/2; } for(temp=0;temp<tznum;temp++)/*当与另一节点间为变压器*/ { if(strcmp(pqv[i].name,tz[temp].name1)==0) { y[i][i].B+=-(1/tz[temp].x*tz[temp].k1/tz[temp].k2+1/tz[temp].x*(tz[temp].k2-tz[temp].k1)/tz[temp].k2); } if(strcmp(pqv[i].name,tz[temp].name2)==0) { y[i][i].B+=-(1/tz[temp].x*tz[temp].k1/tz[temp].k2+1/tz[temp].x*(tz[temp].k1-tz[temp].k2)*tz[temp].k1/pow(tz[temp].k2,2)); } } }}voidshowy(void){ printf("\n"); printf("导纳矩阵Y:\n"); for(temp1=0;temp1<=pqnum+pvnum;temp1++) { for(temp2=0;temp2<=pqnum+pvnum;temp2++) printf("%9.6lf+j%9.6lf",y[temp1][temp2].G,y[temp1][temp2].B); printf("\n"); } printf("\n");}doubledetpqc(void){ doublemax=0; detpq=(double*)malloc((2*pqnum+pvnum)*sizeof(double));/*为detpq误差量*/ detav=(double*)malloc((2*pqnum+pvnum)*sizeof(double)); /*修正量*/ Pi=(double*)malloc((pqnum+pvnum)*sizeof(double)); /**/ Qi=(double*)malloc((pqnum+pvnum)*sizeof(double)); for(temp1=0;temp1<2*(pqnum+pvnum);temp1++) { if(temp1<pqnum+pvnum) /*求P（pq,pv)误差量pqnum+pvnum个*/ { Pi[temp1]=0; for(temp2=0;temp2<=(pqnum+pvnum);temp2++) Pi[temp1]=Pi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*cos(pqv[temp1].ag-pqv[temp2].ag)+y[temp1][temp2].B*sin(pqv[temp1].ag-pqv[temp2].ag)); Pi[temp1]=pqv[temp1].v*Pi[temp1]; detpq[temp1]=pqv[temp1].p-Pi[temp1]; if(fabs(detpq[temp1])>max)max=fabs(detpq[temp1]); } if(temp1>=pqnum+pvnum) /*求Q（pq）误差量pqnum个*/ { Qi[temp1-pqnum-pvnum]=0; for(temp2=0;temp2<=(pqnum+pvnum);temp2++) Qi[temp1-pqnum-pvnum]=Qi[temp1-pqnum-pvnum]+pqv[temp2].v*(y[temp1-pqnum-pvnum][temp2].G*sin(pqv[temp1-pqnum-pvnum].ag-pqv[temp2].ag)-y[temp1-pqnum-pvnum][temp2].B*cos(pqv[temp1-pqnum-pvnum].ag-pqv[temp2].ag)); Qi[temp1-pqnum-pvnum]=pqv[temp1-(pqnum+pvnum)].v*Qi[temp1-pqnum-pvnum]; detpq[temp1]=pqv[temp1-(pqnum+pvnum)].q-Qi[temp1-pqnum-pvnum]; if(fabs(detpq[temp1])>max&&temp1!=2*pqnum+pvnum)max=fabs(detpq[temp1]); } } returnmax;}voidshowdetav(void){ printf("解得:\n"); for(temp1=0;temp1<2*pqnum+pvnum;temp1++) if(temp1<pqnum) printf("Δδ%d：%lf",temp1,detav[temp1]); else printf("ΔV%d:%lf",temp1-pqnum,detav[temp1]); printf("\n");}voidJrc(void){ inti,j; doubleagij; J=(double**)malloc((2*pqnum+pvnum)*sizeof(double*));/*为雅可比矩阵分配内存空间*/ for(temp1=0;temp1<2*pqnum+pvnum;temp1++) J[temp1]=(double*)malloc((2*pqnum+pvnum)*sizeof(double)); for(i=0;i<pqnum+pvnum;i++) { for(j=0;j<pqnum+pvnum;j++)/*Hij*/ { if(i!=j) { agij=pqv[i].ag-pqv[j].ag; J[i][j]=-pqv[i].v*pqv[j].v*(y[i][j].G*sin(agij)-y[i][j].B*cos(agij)); } else J[i][j]=pqv[i].v*pqv[i].v*y[i][i].B+Qi[i]; } for(j=0;j<pqnum;j++)/*Nij*/ { if(i!=j) { agij=pqv[i].ag-pqv[j].ag; J[i][j+pqnum+pvnum]=-pqv[i].v*(y[i][j].G*cos(agij)+y[i][j].B*sin(agij)); } else J[i][j+pqnum+pvnum]=-pqv[i].v*y[i][i].G-Pi[i]; } } for(i=0;i<pqnum;i++) { for(j=0;j<pqnum+pvnum;j++)/*Kij*/ { if(i!=j) { agij=pqv[i].ag-pqv[j].ag; J[i+pqnum+pvnum][j]=pqv[i].v*pqv[j].v*(y[i][j].G*cos(agij)+y[i][j].B*sin(agij)); } else J[i+pqnum+pvnum][j]=pqv[i].v*pqv[i].v*y[i][i].G-Pi[i]; } for(j=0;j<pqnum;j++)/*Lij*/ { if(i!=j) { agij=pqv[i].ag-pqv[j].ag; J[i+pqnum+pvnum][j+pqnum+pvnum]=-pqv[i].v*(y[i][j].G*sin(agij)-y[i][j].B*cos(agij)); } else J[i+pqnum+pvnum][j+pqnum+pvnum]=pqv[i].v*y[i][i].B-Qi[i]; } }}voidshowj(void){ printf("\n"); printf("雅可比矩阵J:\n"); for(temp1=0;temp1<2*pqnum+pvnum;temp1++) { for(temp2=0;temp2<2*pqnum+pvnum;temp2++) printf("%9.6lf",J[temp1][temp2]); printf("\n"); } printf("\n");}voidGauss(void) { intxiao,big,i,j,n=2*pqnum+pvnum; doublebiggest,sta,E; doubleconj=0,conpq=0; for(xiao=0;xiao<2*pqnum+pvnum-1;xiao++) { big=xiao; biggest=fabs(J[xiao][xiao]); for(i=xiao+1;i<2*pvnum+pqnum;i++) /*找出列主元*/ if(fabs(J[i][xiao])>biggest) { big=i; biggest=fabs(J[i][xiao]); } for(i=0;i<2*pqnum+pvnum;i++) /*交换两行*/ { conj=J[xiao][i]; J[xiao][i]=J[big][i]; J[big][i]=conj; } conpq=detpq[xiao]; detpq[xiao]=detpq[big]; detpq[big]=conpq; for(i=xiao+1;i<n;i++) /*消元*/ { sta=J[i][xiao]; for(j=0;j<n;j++) J[i][j]=J[i][j]-sta*J[xiao][j]/J[xiao][xiao]; detpq[i]=detpq[i]-sta*detpq[xiao]/J[xiao][xiao]; } } detav[n-1]=detpq[n-1]/J[n-1][n-1]; /*回代*/ for(i=n-2;i>=0;i--) {E=0; for(j=i+1;j<2*pqnum+pvnum;j++) E=E+J[i][j]*detav[j]; detav[i]=(detpq[i]-E)/J[i][i]; }}voidshowdetpq(void){ printf("\n"); printf("功率不平衡量:\n"); for(temp1=0;temp1<2*pqnum+pvnum;temp1++) {if(temp1<pqnum+pvnum) printf("ΔP%d:%lf",temp1+1,detpq[temp1]); else printf("ΔQ%d:%lf",temp1-pqnum,detpq[temp1]); } printf("\n最大功率不平衡量为（绝对值）：%lf\n",error); printf("\n");}voidshowsolution(void){ printf("\n>>>>>>>>>>>>>>>>>>>迭代结束<<<<<<<<<<<<<<<<<<<<<\n"); printf("\t\t总共进行%d次迭代:\n",diedai); for(temp1=0;temp1<pqnum+pvnum;temp1++) printf("节点%d<%s>电压：幅值%lf相角%lf\n",temp1+1,pqv[temp1].name,pqv[temp1].v,pqv[temp1].ag); for(temp1=0;temp1<(pqnum+pvnum+1);temp1++) { Pi[temp1]=0; for(temp2=0;temp2<=(pqnum+pvnum);temp2++) { Pi[temp1]=Pi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*cos(pqv[temp1].ag-pqv[temp2].ag)+y[temp1][temp2].B*sin(pqv[temp1].ag-pqv[temp2].ag)); Qi[temp1]=Qi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*sin(pqv[temp1].ag-pqv[temp2].ag)-y[temp1][temp2].B*cos(pqv[temp1].ag-pqv[temp2].ag)); } Pi[temp1]=Pi[temp1]*pqv[temp1].v; Qi[temp1]=Qi[temp1]*pqv[temp1].v; printf("节点%d<%s>功率: %8.6lf+j%8.6lf\n",temp1+1,pqv[temp1].name,Pi[temp1],Qi[temp1]); } printf("\n");}charsave2file(void){ FILE*fp; inti; time_tt; charfilename[20],c; printf("\t\t是否保存文件<Y/y>："); fflush(stdin); c=getchar(); if(c!='Y'&&c!='y')gotonext; printf("输入要储存数据的文件名(默认后缀.data)："); fflush(stdin); scanf("%s",filename); strcat(filename,".data"); if((fp=fopen(filename,"w"))==NULL) { printf("无法执行存储操作！！"); exit(0); } time(&t); fprintf(fp,"\t\t\t潮流计算结果\n保存时间：%s\n对应数据文件：%s\n",ctime(&t),fsource); fprintf(fp,"<标幺值>\n"); for(i=0;i<pqnum;i++) fprintf(fp,"PQ节点%s，电压幅值：%10.6lf相角:%10.6lf度有功功率：%10.6lf无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3Pi[i],Qi[i]); for(i=pqnum;i<pqnum+pvnum;i++) fprintf(fp,"PV节点%s，电压幅值：%10.6lf相角:%10.6lf度有功功率：%10.6lf无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3Pi[i],Qi[i]); fprintf(fp,"平衡节点%s，电压幅值：%10.6lf相角:%10.6lf度有功功率：%10.6lf无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3Pi[i],Qi[i]); fclose(fp); printf("\t\t\t文件保存成功!!!\n");next:puts("\t1.返回主界面2.重新计算此文件数据:");nex:fflush(stdin); c=getch(); if(c!='1'&&c!='2')gotonex; returnc;}charoption(void){ charc; system("modecon:cols=75lines=25&color16");cls:puts("***************************************************************"); puts(" 电力系统潮流计算 "); puts("***************************************************************\n\n"); printf(" 选择您的操作：\n\n\t\t1.潮流计算 2.查看数据文件0.退出程序\n\t\t\t\t"); fflush(stdin); c=getch(); returnc;}intlist(charfilename[],inti,chardstn[]) { FILE*fp; charch,c; intshu=0,j=0,k=0; charname[50]; system("dir/x*.data>1.bin2>nul");if((fp=fopen("1.bin","r"))==NULL) { printf("读取列表有误！"); system("pause"); return0; }while((ch=fgetc(fp))!=EOF) { if(ch=='2'&&c==10) { shu++; if(shu!=i) { fseek(fp,47L,1); fgets(name,50L,fp); } if(shu==i) { fseek(fp,35L,1); fgets(dstn,13,fp); if(dstn[0]=='') { fseek(fp,1L,1); fgets(dstn,13,fp); dstn[strlen(dstn)-1]=0; } return0; } if(i==0) { printf("\t\t"); if(numw==shu)printf("\b\b->"); printf("%3d.%s",shu,name); numofw=shu; } ch=10; } c=ch; }fclose(fp);return0;}intchose(void){ chardel[30]={"del"}; FILE*fp; fflush(stdin); c=getch(); if(c=='2'&&numw<numofw)numw++; if(c=='8'&&numw>1)numw--; if(c=='-') { list("1.bin",numw,filename); strcat(del,filen