第一篇：電力系統通用潮流計算C語言程序
#include
//節點號
類型
負荷有功
負荷無功母線數據（類型 1＝PV節點，2＝PQ節點，3＝平衡節點）struct BUS { int busno;int type;float Pd;float Qd;};
//發電機數據 節點號
有功發電
電壓幅值 struct Generator { int busno;float Pg;float Vg;};
//支路信息 節點I 節點J R X B/2 k struct Line { int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2// void fun1(double YG[][50],double YB[][50],double e[],double f[],int type[],int N,double W[],double P[],double Q[],double V[]){ double dP=0,dQ=0,dV=0;int i,j;for(i=0;i double A=0,B=0; for(j=0;j A+=YG[i][j]*e[j]-YB[i][j]*f[j]; B+=YG[i][j]*f[j]+YB[i][j]*e[j]; } dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i]; dP=P[i]-e[i]*A-f[i]*B; W[2*i]=dP; dQ=Q[i]-f[i]*A+e[i]*B; if(type[i]==1) W[2*i+1]=dQ; else W[2*i+1]=dV;} } //Jacobi矩陣// void Jacobi(double YG[][50],double YB[][50],double e[50],double f[50],int type[50],int N ,double Ja[100][101]){ int i,j; for(i=0;i for(j=0;j if(i!=j){ if(type[i]==1){ Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]); Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i]; Ja[2*i+1][2*j]=Ja[2*i][2*j+1]; Ja[2*i+1][2*j+1]=-Ja[2*i][2*j]; } else { Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i]; Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i]; Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0; } } else { double a[50]={0},b[50]={0}; for(int k=0;k a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]); b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]); Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i]; Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i]; if(type[i]==1){ Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i]; Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i]; } else { Ja[2*i+1][2*j]=-2*e[i]; Ja[2*i+1][2*j+1]=-2*f[i]; } } } } } } //高斯消元法解方程組函數// void gauss(double a[][101],int n) { int i,j,k;double c; for(k=0;k c=a[k][k]; for(j=k;j<=n;j++)a[k][j]/=c; for(i=k+1;i c=a[i][k]; for(j=k;j<=n;j++)a[i][j]-=c*a[k][j]; } } a[n-1][n]/=a[n-1][n-1]; for(k=n-2;k>=0;k--) for(j=k+1;j void main(){ ifstream fin;int N=0,GS=0,LD=0,ZLs=0;//節點數 發電機數 負荷數 BUS *B; Generator *G;Line *L; //從文本中讀入原始數據到數組中// fin.open(“C:data.txt”);if(!fin) { cout<<“輸入數據文件不存在！”< getchar(); } int m1[50]={0},m2[50]={0};float m3[50],m4[50],m5[50],m6[50]; int i,j,l; 支路數// for(i=0;;i++){ fin>>m1[i]; if(m1[i]==0)break; fin>>m2[i]>>m3[i]>>m4[i]; N++;} B =new BUS[N];for(i=0;i B[i].busno=m1[i]; B[i].type=m2[i]; B[i].Pd=m3[i]; B[i].Qd=m4[i];} for(i=0;;i++){ fin>>m1[i]; if(m1[i]==0)break; fin>>m4[i]>>m3[i]; GS++;} G =new Generator[GS];for(i=0;i G[i].busno=m1[i]; G[i].Pg=m4[i]; G[i].Vg=m3[i];} for(i=0;;i++){ fin>>m1[i]; if(m1[i]==0)break; fin>>m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i]; ZLs++;} L =new Line[ZLs];for(i=0;i L[i].busi=m1[i]; L[i].busj=m2[i]; L[i].R=m3[i]; L[i].X=m4[i]; } L[i].B=m5[i];L[i].k=m6[i]; LD=N-GS;fin.close(); //節點導納矩陣形成// double YB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i][j]=1;} } for(l=0;l } for(i=0;i for(j=i;j K[i][j]=K[j][i];K[j][i]=1; } for(j=0;j if(i!=j){ YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]); YB[i][i]=YB[i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]); } } } //修正后// for(l=0;l // 求A=e+f// double e[50]={0},f[50]={0};double C[100]={0},D[100]={0};for(i=0;i C[2*i]=1;} else C[2*i]=V[i]; } double W[100]={0},Ja[100][101]={0}; //調用Jacobi函數和高斯函數// for(int t=1;t<10;t++){ for(i=0;i<2*N-2;i++){ e[i]=C[2*i]; f[i]=C[2*i+1]; } fun1(YG,YB,e,f,type,N,W,P,Q,V); double it=fabs(W[0]); for(i=1;i<2*N-2;i++){ if(it it=fabs(W[i]);j=i; } } //中間迭代過程// cout< cout< if(it<0.00001)break; Jacobi(YG,YB,e,f,type,N,Ja); for(i=0;i<2*N-2;i++){ Ja[i][2*N-2]=W[i]; } //高斯消元法解方程// gauss(Ja,2*N-2); for(i=0;i<2*N-2;i++){ D[i]=-Ja[i][2*(N-1)]; C[i]+=D[i]; } } //平衡節點// for(i=0;i double a=0,b=0; for(int j=0;j a+=(YG[i][j]*e[j]-YB[i][j]*f[j]); b+=(YB[i][j]*e[j]+YG[i][j]*f[j]); } P[i]=e[i]*a+f[i]*b; Q[i]=f[i]*a-e[i]*b; } //支路// double PZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[100][101]={0};double x1=0,x2=0,y1=0,y2=0,I2=0;for(int k=0;k i=L[k].busi-1;j=L[k].busj-1;x1=e[i]/L[k].k-e[j];y1=f[i]/L[k].k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZL[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k].R;qx[i][j]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k].B);x1=e[j]*L[k].k-e[i];y1=f[j]*L[k].k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZL[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k].R;qx[j][i]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k].B); } //全網數據// int high=1,low=1; double PG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0]*e[0]+f[0]*f[0]);for(k=0;k Vh=sqrt(e[k]*e[k]+f[k]*f[k]);high=k+1;} if(sqrt(e[k]*e[k]+f[k]*f[k]) Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;} } //輸出數據到文件databak.txt// ofstream fout; fout.open(“C:databak.txt”); fout<<“節點”< fout< for(i=0;i fout< } for(j=0;j i=G[j].busno-1; fout< fout<<“支路 ”< fout< for(k=0;k i=L[k].busi-1;j=L[k].busj-1; fout< fout<<“全網數據”< fout< fout< 電力系統潮流計算c語言程序，兩行，大家可以看看，仔細研究，然后在這個基礎上修改。謝謝 #include “stdafx.h” #include #include“Complex.h” #include“wanjing.h” #include“gauss.h” using namespace std; int _tmain(int argc, _TCHAR* argv[]){ int i; //i作為整個程序的循環變量 int N=Bus::ScanfBusNo();//輸入節點個數 int L=Line::ScanflineNo();//輸入支路個數 if((L&&N)==0){return 0;} //如果找不到兩個文件中的任意一個，退出 Line *line=new Line[L];//動態分配支路結構體 Line::ScanfLineData(line);//輸入支路參數 Line::PrintfLineData(line,L);//輸出支路參數 Bus *bus=new Bus[N];//動態分配結點結構體 for(int i=0;i bus[i].Sdelta.real=0; bus[i].Sdelta.image=0;} Bus::ScanfBusData(bus);//輸入節點參數 Bus::PrintfBusData(bus,N);//輸出結點參數 Complex **X;X=new Complex *[N];for(i=0;i Bus::JisuanNodeDnz(X,line,bus,L,N);//計算節點導納矩陣 Bus::PrintfNodeDnz(X,N);//輸出節點導納矩陣 int NN=(N-1)*2;double **JacAug;JacAug=new double *[NN];for(i=0;i double *x;x=new double[NN];int count=1; LOOP: Bus::JisuanNodeI(X,bus,N);//計算節點注入電流 Bus::JisuanNodeScal(X,bus,N);//計算節點功率 Bus::JisuanNodeScal(X,bus,N);//計算節點功率 Bus::JisuanNodeSdelta(bus,N);//計算節點功率差值 Bus::PrintfNodeScal(X,bus,N);//輸出節點功率差值 int icon=wehcon1(bus,N);//whether converbence看迭代是否結束 if(icon==1){ cout<<“icon=”< Bus::JisuanJacAug(JacAug,X,bus,N);//計算雅可比增廣矩陣 // Bus::PrintfJacAug(JacAug,N); gauss::gauss_slove(JacAug,x,NN);//解方程組求出電壓差值 Bus::ReviseNodeV(bus,x,N);//修正節點電壓 // Bus::PrintfNodeV(bus,N); count++; goto LOOP;} else { for(i=0;i { int statemp,endtemp; Complex aa,bb,cc,dd,B; B.real=0; B.image=-line[i].B; statemp=line[i].start; endtemp=line[i].end; aa=Complex::productComplex(Complex::getconj(bus[statemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].stoe=Complex::productComplex(bus[statemp-1].V,dd); aa=Complex::productComplex(Complex::getconj(bus[endtemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].etos=Complex::productComplex(bus[endtemp-1].V,dd); } cout<<“icon=”< Bus::JisuanNodeScal(X,bus,N);//計算節點功率 for(i=0;i { bus[i].Scal.real = bus[i].Scal.real + bus[i].Load.real;//發電機功率=注入功率+負荷功率 bus[i].Scal.image= bus[i].Scal.image+ bus[i].Load.image; bus[i].V=Complex::Rec2Polar(bus[i].V); } cout<<“====節點電壓===============發電機發出功率======”< for(i=0;i { cout<<“節點”<<(i+1)<<'t'; Complex::PrintfComplex(bus[i].V); coutt(bus[i].Scal.real); coutt(bus[i].Scal.image); cout< } cout<<“======線路傳輸功率==========”< for(i=0;i { int statemp,endtemp; statemp=line[i].start; endtemp=line[i].end; cout< Complex::PrintfComplex(Complex::ComDivRea(line[i].stoe,0.01)); Complex::PrintfComplex(Complex::ComDivRea(line[i].etos,0.01)); cout< } } return 0;} #include “stdafx.h” #include #include“Complex.h” #include“wanjing.h” #include“gauss.h” using namespace std; int _tmain(int argc, _TCHAR* argv[]){ int i; //i作為整個程序的循環變量 int N=Bus::ScanfBusNo();//輸入節點個數 int L=Line::ScanflineNo();//輸入支路個數 if((L&&N)==0){return 0;} //如果找不到兩個文件中的任意一個，退出 Line *line=new Line[L];//動態分配支路結構體 Line::ScanfLineData(line);//輸入支路參數 Line::PrintfLineData(line,L);//輸出支路參數 Bus *bus=new Bus[N];//動態分配結點結構體 for(int i=0;i bus[i].Sdelta.real=0; bus[i].Sdelta.image=0;} Bus::ScanfBusData(bus);//輸入節點參數 Bus::PrintfBusData(bus,N);//輸出結點參數 Complex **X;X=new Complex *[N];for(i=0;i Bus::JisuanNodeDnz(X,line,bus,L,N);//計算節點導納矩陣 Bus::PrintfNodeDnz(X,N);//輸出節點導納矩陣 int NN=(N-1)*2;double **JacAug;JacAug=new double *[NN];for(i=0;i double *x;x=new double[NN];int count=1; LOOP: Bus::JisuanNodeI(X,bus,N);//計算節點注入電流 Bus::JisuanNodeScal(X,bus,N);//計算節點功率 Bus::JisuanNodeScal(X,bus,N);//計算節點功率 Bus::JisuanNodeSdelta(bus,N);//計算節點功率差值 Bus::PrintfNodeScal(X,bus,N);//輸出節點功率差值 int icon=wehcon1(bus,N);//whether converbence看迭代是否結束 if(icon==1){ cout<<“icon=”< Bus::JisuanJacAug(JacAug,X,bus,N);//計算雅可比增廣矩陣 // Bus::PrintfJacAug(JacAug,N); gauss::gauss_slove(JacAug,x,NN);//解方程組求出電壓差值 Bus::ReviseNodeV(bus,x,N);//修正節點電壓 // Bus::PrintfNodeV(bus,N); count++; goto LOOP;} else { for(i=0;i { int statemp,endtemp; Complex aa,bb,cc,dd,B; B.real=0; B.image=-line[i].B; statemp=line[i].start; endtemp=line[i].end; aa=Complex::productComplex(Complex::getconj(bus[statemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].stoe=Complex::productComplex(bus[statemp-1].V,dd); aa=Complex::productComplex(Complex::getconj(bus[endtemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].etos=Complex::productComplex(bus[endtemp-1].V,dd); } cout<<“icon=”< Bus::JisuanNodeScal(X,bus,N);//計算節點功率 for(i=0;i { bus[i].Scal.real = bus[i].Scal.real + bus[i].Load.real;//發電機功率=注入功率+負荷功率 bus[i].Scal.image= bus[i].Scal.image+ bus[i].Load.image; bus[i].V=Complex::Rec2Polar(bus[i].V); } cout<<“====節點電壓===============發電機發出功率======”< for(i=0;i { cout<<“節點”<<(i+1)<<'t'; Complex::PrintfComplex(bus[i].V); coutt(bus[i].Scal.real); coutt(bus[i].Scal.image); cout< } cout<<“======線路傳輸功率==========”< for(i=0;i { int statemp,endtemp; statemp=line[i].start; endtemp=line[i].end; cout< Complex::PrintfComplex(Complex::ComDivRea(line[i].stoe,0.01)); Complex::PrintfComplex(Complex::ComDivRea(line[i].etos,0.01)); cout< } } return 0;} #include class Complex//定義復數類 { public: double real;double image;int RecPolar;//0表示直角坐標，1表示極坐標 static Complex CaddC(Complex c1,Complex c2);//求兩個復數和 static Complex subComplex(Complex c1,Complex c2);//求兩個復數差 static Complex productComplex(Complex c1,Complex c2);//求兩個復數積 static Complex divideComplex(Complex c1,Complex c2);//求兩個復數商 static Complex ComDivRea(Complex c1,double r2);//除數 static Complex getconj(Complex c1);//求一個復數共軛 static Complex getinverse(Complex c1);//取倒數 static double getComplexReal(Complex c1);//求一個復數實部 static double getCompleximage(Complex c1);//求一個復數虛部 static void PrintfComplex(Complex c1);//顯示一個復數 static void PrintfmultiComplex(Complex C,int N);//顯示多個復數 static void zeroComplex(Complex c1);//將復數復零 static Complex Rec2Polar(Complex c1);//取極坐標 Complex(){ RecPolar=0;} }; Complex Complex::Rec2Polar(Complex c1)//極坐標表示 { Complex Node;Node.real=sqrt(c1.real*c1.real+c1.image*c1.image);Node.image=atan2(c1.image,c1.real)*180/3.1415926;Node.RecPolar=1;return Node;} Complex Complex::CaddC(Complex c1,Complex c2)//復數加法 { Complex Node; Node.real=c1.real+c2.real; Node.image=c1.image+c2.image; return Node;} Complex Complex::subComplex(Complex c1,Complex c2)//復數減法 { Complex Node; Node.real=c1.real-c2.real; Node.image=c1.image-c2.image; return Node;} Complex Complex::productComplex(Complex c1,Complex c2)//復數乘法 { Complex Node; Node.real=c1.real*c2.real-c1.image*c2.image; Node.image=c1.image*c2.real+c2.image*c1.real; return Node;} Complex Complex::divideComplex(Complex c1,Complex c2)//復數除法 { Complex Node; Node.real=(c1.real*c2.real+c1.image*c2.image)/(pow(c2.real,2)+pow(c2.image,2));Node.image=(c1.image*c2.real-c1.real*c2.image)/(pow(c2.real,2)+pow(c2.image,2));return Node;} Complex Complex::ComDivRea(Complex c1,double r1)//復數除數 { Complex Node;Node.real=c1.real/(r1);Node.image=c1.image/(r1);return Node;} Complex Complex::getconj(Complex c1)//取共軛 { Complex Node; Node.real=c1.real;Node.image=-c1.image; return Node;} Complex Complex::getinverse(Complex c1)//取倒數 { Complex Node;Node.real=1;Node.image=0;Node=(Complex::divideComplex(Node,c1));return Node;} double Complex::getComplexReal(Complex c1)//取實部 { return c1.real;} double Complex::getCompleximage(Complex c1)//取虛部 { return c1.image;} void Complex::PrintfComplex(Complex c1)//按直角坐標輸出 { if(c1.RecPolar==0){ cout.precision(6); cout.width(8); cout.setf(ios::right); cout< ”; cout.precision(6); cout.width(8); cout.setf(ios::left); cout< ”;} else { cout< Complex::zeroComplex(Complex c1)//清零 { c1.real=0;c1.image=0;} class gauss { public: static void gauss_slove(double **a,double *x,int NN);static void gauss_output();}; void gauss::gauss_slove(double **a,double *x,int NN){ int n,i,j,k,*pivrow,**pivcol;double eps,pivot,sum,aik,al; n=NN;pivrow=new int[n];pivcol=new int *[n]; for(i=0;i pivot= fabs(a[k][k]); pivrow[k]=k;//行 pivcol[k][0]=k;pivcol[k][1]=k;//列n*2矩陣 for(i=k;i { for(j=k;j { if(pivot { pivot=fabs(a[i][j]); pivrow[k]=i;//行 pivcol[k][1]=j;//列 } } } if(pivot { cout<<“error”< getchar(); exit(0); } if(pivrow[k]!=k)//行變換 { for(j=k;j<(n+1);j++) { al=a[pivrow[k]][j]; a[pivrow[k]][j]=a[k][j]; a[k][j]=al; } } if(pivcol[k][1]!=k)//列變換 { for(i=0;i { al=a[i][pivcol[k][1]]; a[i][pivcol[k][1]]=a[i][k]; a[i][k]=al; } } if(k!=(n-1))//將矩陣化為上三角形 式 { for(i=(k+1);i { aik=a[i][k]; for(j=k;j<(n+1);j++) { a[i][j]-=aik*a[k][j]/a[k][k]; } } } } x[n-1]=a[n-1][n]/a[n-1][n-1];//解方程 for(i=(n-2);i>=0;i--){ sum=0; for(j=(i+1);j { sum +=a[i][j]*x[j];0.182709 0.016894-0.0310701 -0.0402051 0.156702 -0.0355909-0.0668055 -0.00703229-0.0886481 -0.0129814-0.0390805 -0.0135062-0.1023 -0.0460568 -0.0342827 -0.00382402-0.102896 -0.0184062 } x[i]=(a[i][n]-sum)/a[i][i];} for(k=(n-2);k>=0;k--){ al=x[pivcol[k][1]]; x[pivcol[k][1]]=x[pivcol[k][0]]; x[pivcol[k][0]]=al;} cout<<“節點電壓修正量”< cout< } ====節點功率計算值==== 0.935261 -0.159048 0.573909 0.0789973-0.00289889 -0.00796623-0.0791247 -0.0168362-0.436255 -0.0580392 0.0359139 -0.0106592-0.229118 -0.0885419-0.136179 -0.148207 0.0446243 0.0111298-0.0223764 -0.00695775-0.0237482 -0.198318 -5.24266e-015 -0.0354071 -0.0925078 -1.05629e-015 -0.0391348 0.014529 0.00158644 -0.0258771 -0.109514 icon=1進行第2次迭代 節點電壓修正量 =================-0.00164889-0.000540034-0.00261067-0.00532027-0.00235315-0.00600971-0.00189677-0.00643874-0.0023631-0.00650659-0.00170949-0.0074907-0.00164545-0.00485415-0.00493977-0.0119042-0.00331285-0.0175611-0.00207908 -0.00347744-0.0869347-9.48909e-015-0.0110778-0.0538236-7.53784e-016-0.0168097 7.049e-005-0.00146487-0.00458276 0.00251645 -0.00336375-0.00530645-0.0147816-0.000326161-0.00640487-0.00251701-0.0169829-0.00175286-0.0174333-0.0239063 -0.0119192-0.076014 -0.0160104-0.441997 -0.0750285 0.000250012 3.72542e-005-0.228052 -0.108844-0.100078 -0.105634 0.000410707 0.000378067-0.057497 -0.0195879 0.200039 0.0582563-0.00307326-0.0163809-0.00232773-0.0175806 8.74293e-005-0.0192018 0.000558996-0.0197776-0.000247851-0.0193784-0.00115346-0.0185848-0.00127275-0.0186244-0.00010108-0.0188966 0.000553585-0.0200901-3.76315e-005-0.0208303 0.00308341-0.0219386-0.00195916-0.0205356-0.00184757-0.0076401 0.00197593-0.0245534 0.00434657-0.027534 ====節點功率計算值==== 0.98623 -0.134163 0.583136 0.166278-0.111173 0.199792 -0.0621041 -0.0821379 -0.0350785 -0.0902383 -0.0320461 -0.0951562 -0.0220362 -0.175458 4.72557e-015 -0.0320661 -0.0871134 -7.03489e-017 -0.0350769 0.000273455 1.51804e-005 -0.0240417 -0.10604 icon=1進行第3次迭代 節點電壓修正量 =================-2.67079e-005-2.30128e-006-2.20543e-005-6.00686e-005-2.33043e-005-6.85601e-005-3.22294e-005-2.61107e-005-2.80198e-005-6.6167e-005-2.34528e-005 -0.0739846 0.0227868-0.0158709-0.0248173-0.0179447-0.0578368-0.00890719-0.0337091-0.00693706-0.111601 1.21429e-014-0.0159145-0.0667319 9.24355e-016-0.0228592 7.10354e-005-6.6188e-006-0.00889343-0.0184098 -5.66132e-005-4.4646e-005-1.74668e-005-4.50947e-005-0.000181763-3.81763e-006-0.000286581-6.68993e-005-1.28441e-005-5.17172e-005-0.000223284-4.54717e-005-2.47586e-005 4.32335e-007-0.000258494 1.82635e-005-0.000272051-6.95195e-006-0.000251969 1.11318e-005-0.000279418 5.74737e-005-0.000307368 6.86998e-005-0.000320274 5.38112e-005-0.00031447 3.59531e-005-0.00030494 3.37607e-005-0.000307449 5.26532e-005-0.000310721 6.92761e-005-0.000350373 5.60942e-005-0.00040977 0.000123641-0.000440259 1.36149e-005-0.000426973-1.70227e-005-9.37794e-005 0.000113675-0.000544011 0.000176034-0.000636202 ====節點功率計算值==== 0.986878 -0.133979 0.583 0.167193-0.024 -0.012-0.076 -0.016-0.442 -0.0748606 1.43501e-008 1.07366e-008-0.228 -0.109 -0.0999999 -0.104049 4.51318e-008 8.98835e-008-0.0579999 -0.0199999 0.2 0.0591018-0.112 -0.0749997 0.2 0.0242519-0.062 -0.016-0.082 -0.025-0.035 -0.018 -0.0900001 -0.058-0.032 -0.00899997-0.095 -0.0339999-0.022 -0.00699998-0.175 -0.112 -6.07156e-015 -1.19217e-014-0.032 -0.016-0.087 -0.0669999 7.03078e-017 -9.23979e-016-0.035 -0.0229999 1.09492e-007 4.45699e-008 1.54958e-009 -2.01531e-010-0.024 -0.00899994-0.106 -0.0189996 icon=0,迭代結束。 ====節點電壓===============發電機發出功率====== 節點1 1.05 0。 98.6878-13.3979 節點2 1.045 -1.846。 29.4193 節點3 1.02384-3.83352。 0 節 點25 1.01216-9.68486。 0 0 0 節點4 1.01637-4.55698。 0 節 點26 0.994393 -10.1089。 0 0 0 節點5 1.01 -6.48617。 節 點27 1.02012-9.42025。 0 11.5139 0 節點6 1.01332-5.38073。 0 節 點28 1.00992-5.86244。 0 0 0 節點7 1.00489-6.38368。 0 節 點29 1.00022-10.6579。 0 0 節點8 19.5951 節點9 0 節點10 0 節點11 5.91018 節點12 0 節點13 2.42519 節點14 0 節點15 0 節點16 0 節點17 0 節點18 0 節點19 0 節點20 0 節點21 0 節點22 0 節點23 0 節點24 0 1.01 -5.62974。 1.03905-6.78143。 1.03595-8.69362。 -4.5962。 1.04711-7.80323。 1.05 -6.34392。 1.03242-8.7401。 1.02788-8.86784。 1.03458-8.45044。 1.03051-8.83678。 1.01845-9.5141。 1.01604-9.70326。 1.02022-9.50938。 1.0237-9.17478。 1.02432-9.17024。 1.01802-9.36719。 1.01339-9.68362。 0 20 節 點30 0.988705 -11.5464。 0 0 0 ====== 線路傳輸功率========== 2to1 -57.7373 5.41674i 58.3454 0 -15.1827i 3to1 -39.659 -7.75964i 40.3424 1.78481i 4to2 -30.87 -9.74186i 31.4153 0 3.58352i 4to3 -37.0772 -7.78596i 37.259 6.55964i 5to2 -44.3717 -9.78456i 45.2968 0 4.84242i 6to2 -38.4766 -8.22625i 39.3252 0 2.87667i 6to4 -34.946 1.92384i 35.0885 0 -3.28202i 7to5 -0.16304 -6.41767i 0.171702 0 2.2985i 7to6 -22.637 -4.48233i 22.7745 0 1.44238i 8to6 -11.8939 -5.48098i 11.913 0 3.70557i 6to9 12.3737 -12.3826i -12.3737 0 13.0033i 6to10 10.9107 -3.80907i -10.9107 0 4.53223i 11to9 5.91018i 0 -5.08963i 10to9 -32.652 -2.3712i 32.652 0 3.46974i 4to12 23.5411 -11.5375i -23.5411 0 13.2407i 13to12 2.42519i 1.05 -1.90978i 1.66484i 14to12 -7.9019 -2.06732i 7.97894 30to29 -3.6702 -0.542564i 3.70398 2.22749i 0.606393i 15to12 -18.254 -5.74885i 18.4835 28to8 -1.89152 -3.79982i 1.89395 6.20089i-4.9239i 16to12-7.53872 -2.90237i 7.59633 28to6 -14.7868 -2.82565i 14.8234 3.02352i 0.294601i 15to14-1.69544 -0.461488i 1.70189 請按任意鍵繼續...0.467323i 17to16-4.03014 1.10238i 18to15-6.08074 1.46028i 19to18-2.87549 0.478389i 20to19 6.6418-2.93222i 20to10 -8.8418 3.85077i 17to10-4.96987 4.76656i 21to10-16.1562 9.42843i 22to10-7.87782 4.21401i 22to21 1.34443-2.01837i 23to15-5.59369 2.25006i 24to22-6.48186 2.08163i 24to23-2.38596 0.579814i 25to24-0.167617 0.281364i 26to25 -3.5 2.3674i 27to25 3.39433-2.08638i 28to27 16.1446 3.13006i 29to27-6.10398 1.67047i 30to27-6.92979-1.07089i-1.37839i-0.467767i 2.96679i-3.66679i-4.72911i-9.18162i-4.10132i 2.01969i-2.17981i-2.00141i-0.56401i -0.28102i-2.29999i 2.11848i-2.10093i-1.50639i -1.3574i 4.03872 6.12096 2.88074 -6.62452 8.9242 4.98423 16.2709 7.93248 -1.34378 5.62846 6.53339 2.39369 0.167814 3.54513 -3.37751 -16.1446 6.19083 7.09313 高等電力系統分析 IEEE30節點潮流程序 班級：電研114班 姓名：王大偉 學號：2201100151 ////////////////////////////////////////////////////////////////////// // PQ分解法潮流 // //文件輸入格式：節點總數n(包括聯絡節點),支路數zls //節點數(發電機和負荷)nb,接地電抗數mdk,迭代精度eps //考慮負荷靜特性標志kk2(0考慮),平衡節點號，優化標志(0不優化) //最大迭代次數it1,支路左右節點號izl[],jzl[],支路電阻zr[],電抗zx[] // //支路容納zyk[],節點號nob[]及標志nobt[](0-PQ,-1-PV) //發電機和負荷有功、無功pg[],qg[],pl[],ql[] //電壓v0[](pv節點輸入實際值，PQ節點任輸入一值) //電抗節點號idk[],電抗值dkk[] ////////////////////////////////////////////////////////////////////// #include “math.h” #include “stdio.h” #define NS 2000 //最大節點數 #define NS2 NS * 2 #define NS4 #define ZS #define ZS2 1000 3000 ZS * 2 //NS4、NS必須大于2*zls。//最大支路數 // // // // // // // #define DKS 200 //最大電抗器數 #define N2 ZS * 4 #define N3 ZS * 8 + NS * 4 FILE *fp1, *fp2;char inname[12], outname[12];// fp1輸入數據文件指針 fp2輸出文件指針 // inname[]輸入數據文件名 outname[]輸出數據文件名 int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;// 節點總數n(包括聯絡節點) 支路數(回路數)zls 節點數nb(發電機和負荷) // 接地電抗數mdk 精度eps 平衡節點號mpj // 節點優化(標志)bnsopton(=0節點不優化，!=0節點優化) // 最大迭代次數it1 最低電壓或最大功率誤差節點號dsd // 負荷靜特性標志(=0考慮負荷靜特性)// 支路數(雙回線算一條支路)int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];// izl[],jzl[],idk[]:分別存放左、右節點號和電抗器節點號。// yds[]存放各行非零非對角元素的個數。 // ydz[i]是第 i 行第一個非零非對角元素的首地址，// 即在所有非零非對角元素中的次序號 // iy[]存放列足碼。 int nnew[NS4], old[NS], nob[NS], nobt[NS];// nnew[],old[]存放的是新、舊節點號。// nnew[i]中為i對應的新號 // nob[]存放的是節點號。nobt[]存放的是節點類型, 0: pq節點,-1: pv節點。double eps, dsm, vmin, dph, dqh, af[3];// eps迭代收斂精度，dsm最大功率誤差 // vmin:系統最低電壓值。dph,dqh:系統有、無功損耗。 // af[0]和af[1]分別是負荷有功功率、無功功率靜態特性系數。double v00; // v00: 系統平均電壓 ci,cj分別作為節點i,j的電壓相角的臨時存儲單元。double zr[ZS], zx[ZS], zyk[ZS], dkk[DKS], gii[NS], bii[NS], yg[ZS2], yb[ZS2];double pg[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[]是電壓相角。double w[NS2], kg[3], b[NS2];int newsort[NS4];// newsort[i]存放i對應的老號 void initial();void pqflow();void out();void dataio();void bnsopt();void zlsort(int* nnew);void printo();void printy();void y2();void ya0();void yzb();void jdgl(int kq0);void bbhl(int kq0);void calc();int iabs(int a);void branch_output();void newval(double* aa);void printc();void iswap();void swap();void printf2(double* aa, double* bb, int n);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 find(int k, int a[], int* z);void yzb(int t, int* iu, double* u, double* di, int* nfd);int isgn(int a, int b);void yy1();void y3();void newtoold(); int main(void){ initial();//初始化 pqflow();//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待搜索數組，參數3返回k在a[]中的次序號// int i;for(i = 1;i <= n;i ++) if(iabs(k)== a[i]) { *z = i; return 1; } return 0;} void oldtonew(){ //**** 本函數將輸入數據中的節點號變成從1開始的連續節點號 ****// int i, 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] = %dn”, 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] = %dn”, 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; } } } void initial(){ //**** 本函數進行初始化工作 ****// int i, 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);// sort the r,x and b 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。 // exchange the node before and after sort for(i = 1;i <= n;i ++) yds[i] = 0; // the immediate 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[] is type of node if(nobt[i] ==-1) v[i] = v0[i]; // nob[] is serials numbe //nobt[] =-1: pv節點，v0[]存放的是最后一個節點數據，//對于pv節點，即為該點應維持的電壓值。 //nobt[] = 0: pq節點，v0[]存放的是最后一個節點數據，//對于pq節點，即為系統平均電壓值。printo();//輸出af[]、v00和節點排序后的支路、節點和 //接地電抗數據(僅僅查看中間結果)ya0();//獲得yds[]、ydz[]、列足碼iy[]。(P407)} void printo(){ //**** 輸出af[]、v00和節點排序后的支路、節點和接地電抗數據 ****// int i;fprintf(fp2, “n ******AF AND V0 ******n”);fprintf(fp2, “n %7.3f%7.3f%7.3fn”, af[0], af[1], v00);printc('-', 78);fprintf(fp2, “nn *******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, “nn*******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,“nn******DKK******n”);for(i = 1;i <= mdk;i ++) { fprintf(fp2, “n”); fprintf(fp2, “%8d%8d%7.4f”, idk[i], old[idk[i]], dkk[i]); } } void dataio() { //**** 系統數據初始化 ****// int i;af[0] = 0.6; af[1] = 2.0;//af[0]和af[1]分別是負荷有功功率、無功功率靜態特性系數。 v00 = 1.000;//系統平均電壓 printf(“nplease input the name of data filen”);scanf(“%s”, inname);fp1 = fopen(inname, “r”);printf(“nplease output the name of data filen”);scanf(“%s”, outname);fp2 = fopen(outname, “w”);fscanf(fp1, “%d %d %d %d”, &n, &zls, &nb, &mdk); // the number of node ,branches, node fscanf(fp1, “%lf %d %d %d %d”, &eps, &kk2, &mpj，&bnsopton, &it1); //precision, swing node,sort the node,iteration numbers 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);} void pqflow() { //**** PQ分解法計算潮流,程序框圖見P164圖3-16(從第 7 步起)****// int kq0, iu1[N2], nfd1[NS], iu2[N2], nfd2[NS];int i, t; double u1[N2], u2[N2], di1[NS], di2[NS];yy1(); yzb(0, iu1, u1, di1, nfd1);//form the B matrix of P-0 iteration y2();yzb(1, iu2, u2, di2, nfd2);//form the B matrix of Q-V iteration t = 0; kq0 = 0;kg[0] = kg[1] = 1;do { jdgl(kq0);// calculating the power bbhl(kq0);// find out the maxi if(kq0 == 0) printf(“P: %dt%dt%fn”, t, dsd, dsm);else printf(“Q: %dt%dt%fn”, 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];} } 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);} void out() { //**** 本函數輸出節點和支路數據 ****// zlsort(old); // recover the data if sorted // newtoold();node_output();// node data branch_output(); //branch data printc('-', 78);printc('*', 78);fprintf(fp2, “n”);} void newval(double* aa){ //**** 本函數將舊號換成新號 ****// int i, 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];} void yzb(int t, int* iu, double* u, double* di, int* nfd){ //**** 本函數求因子表 ****// //參數1為標志(t=0 求B',t=1求B'')// //參數2因子表上三角矩陣非零非對角元素的列足碼 //參數3因子表上三角矩陣非零非對角元素的數值 //參數4因子表上三角矩陣對角元素 //參數5因子表上三角各行非零元素個數 int i, j, k, i1, i2; int jj, jj1, jj2, im, x, fd[NS];double ai, b[NS];nfd[1] = 1;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(im = 1;im <= i1;im ++){ jj1 = nfd[im]; jj2 = nfd[im + 1]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, “%10.5f%8i”, u[i], iu[i]);} fprintf(fp2, “n********DI********”);printf1(di, n);} void printf1(double* aa, int n){ //**** 本函數輸出aa[i],i=1-n ****// int i;for(i = 1;i <= n;i ++) { if(i % 5 == 1) fprintf(fp2, “n”); fprintf(fp2, “%9.5f”, aa[i]); } fprintf(fp2, “nn”);} void calc(int* iu, double* u, double* di, int* nfd, double* b){ //**** 本函數利用因子表解線形方程組。(P417圖F1-9) double bi;int i, 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]1; for(j = i2;j >= i1;j--) { k = iu[j]; bi = bi1;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]); } } } void bnsopt(){ iswap(&jzl[i], &jzl[ip]);swap(&zr[i], &zr[ip]);swap(&zx[i], &zx[ip]);swap(&zyk[i], &zyk[ip]);//**** 節點優化 ****// int ii1, ii2, zls2, nomax;int i, j, l, k1, k; int temp;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 <= n1;//總支路數(雙回線算一條支路)ydz[1] = 1;for(i = 1;i <= ngij;yb[ll] = yb[ll]gij / yk;yb[ll] = yb[ll]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]bi;//非對角元素 j 的功率 } } void bbhl(int kq0) { //**** 本函數計算各節點的功率誤差，求最大功率誤差dsm ****// //**** 和常數項b[i]。(程序框圖見P423)****// int i;double vi, vj;double pl0, 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 *((vipl0; // pv節點 if(i == mpj && kq0 == 0) pg[i] = w[i]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)不參與 // 求最大功率誤差和常數項的運算。 } } node_output() // { //**** 輸出節點數據和最小電壓幅值、相角(角度)及其節點號 ****// //****(程序框圖見p426 F1-16)****// double vi, ci;int i, j, oldnumber;printc('+', 72);fprintf(fp2, “n%5s%8s%10s%11s%11s%11s%11sn”, “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);} void branch_output() { //**** 本函數輸出支路數據。(程序框圖見p428 F1-17) int ii, jj;double r, x, yk, zf, vi, vj, ci, cj; int i, j, l;double de, 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%12sn ”, “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 = eifj; zf = r * r + x * x; fii =(de * r + df * x)/ zf; fir =(df * rfir * ei; pji =-fii * ejvi * vi * yk; qji = qji1;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]); } } for(i = 1;i <= nb;i ++){ k = nob[i]; nob[i] = newsort[k]; printf(“nob[%d] = %dn”, i, nob[i]);} for(i = 1;i <= nb-1;i ++){ for(j = i+1;j <= nb;j ++) { if(nob[i] > nob[j]) { k = nob[i]; nob[i] = nob[j]; nob[j] = k; } } } printf(“nob[%d] = %dn”, i, nob[i]);} void yy1() { //**** 本函數形成節點導納陣(不包括接地支路)****// int j1;double r, x, yk, zf, gij, bij;int i, j, i1, l, ll; for(i = 1;i <= n;i ++) { gii[i] = 0.0; bii[i] = 0.0; }// 導納陣對角元(與節點一一對應)先清零。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;bij =-1/x;yg[ll] = yg[ll]bij;gii[i] = gii[i] + gij;bii[i] = bii[i] + bij;gii[j] = gii[j] + gij;bii[j] = bii[j] + bij; if((i!= iabs(izl[l + 1]))||(j!= iabs(jzl[l + 1]))) ll++;} // 打印導納矩陣。對角元實部為gii，虛部為bii，// 非零非對角元實部為yb[]，虛部為yb[]，列足碼為iy[]。 fprintf(fp2, “*******GII(1),BII(1)********n”);printf2(gii,bii,n);} void y3() { //**** 本函數形成節點導納陣，追加接地支路 ****// int j1;double r, x, yk, zf, gij, bij;int i, j, i1, l, ll, kk = 0; for(i = 1;i <= mdk;i ++) { j = idk[i]; bii[j] = bii[j]gij / yk; yb[ll] = yb[ll]-bij / yk; } if((i!= iabs(izl[l + 1]))||(j!= iabs(jzl[l + 1]))) { ll++; kk = 0; } else kk = 1;} // 打印導納矩陣。對角元實部為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]); } } #include void main(){ cout<<“請輸入末端負荷：（先有功Pd再無功Qd）n”;double Pd;double Qd;cin>>Pd>>Qd;cout<<“請輸入末端大概的電壓值Ud：n”;double Ud;cin>>Ud;double data[300][8];//分別為0前面有功，1前面無功，2壓降，3中間有功，4中間無功,5變比,6后面有功，7后面無功 int ii;ii=0;data[ii][0]=0.0;data[ii][1]=0.0;data[ii][2]=Ud;data[ii][3]=Pd;data[ii][4]=Qd;data[ii][5]=1.0;data[ii][6]=0.0;data[ii][7]=0.0;ii++;cout<<“請按從末端到首段的順序輸入線路的參數（分為--節點、輸電線路、變壓器）：n”; while(1){ double a; double b; double c; double d; double e; double f; int x; cout<<“請輸入數值，0代表退出，1代表節點，2代表輸電線路，3代表變壓器：n”; cin>>x; if(x!=0&&x!=1&&x!=2&&x!=3) { cout<<“選擇有誤，請重新輸入！n”; cout<<“請輸入數值，0代表退出，1代表節點，2代表輸電線路，3代表變壓器：n”; cin>>x; } if(x==0) break; if(x==1) { cout<<“請輸入節點參數值n（有功功率P、無功功率Q）：n”; cin>>a>>b; Node N[100]; int i=0; N[i].Pn=a; N[i].Qn=b; data[ii][0]=0.0; data[ii][1]=0.0; data[ii][2]=0.0; data[ii][3]=N[i].Pn; data[ii][4]=N[i].Qn; data[ii][5]=1.0; data[ii][6]=0.0; data[ii][7]=0.0; i++; ii++; } if(x==2) { cout<<“請輸入輸電線路參數值n（線路阻抗R、線路感抗X、電導B）：n”; cin>>a>>b>>c; Transmission_line Tm[100]; int j=0; Tm[j].Rl=a; Tm[j].Xl=b; Tm[j].Bl=c; data[ii][0]=0.0; data[ii][1]=Tm[j].GetQf(data[ii-1][2]); data[ii][6]=0.0; data[ii][7]=Tm[j].GetQb(data[ii-1][2]); data[ii][2]=Tm[j].GetUl(data[ii-1][3]+data[ii-1][1]+data[ii][7],data[ii-1][4]+data[ii-1][0]+data[ii][6],data[ii-1][2]); //參數分別為前節點的有功、無功、電壓 data[ii][3]=Tm[j].GetPl(data[ii-1][3]+data[ii-1][1]+data[ii][7],data[ii-1][4]+data[ii-1][0]+data[ii][6],data[ii-1][2]); data[ii][4]=Tm[j].GetQl(data[ii-1][3]+data[ii-1][1]+data[ii][7],data[ii-1][4]+data[ii-1][0]+data[ii][6],data[ii-1][2]); data[ii][5]=1.0; j++; ii++; } if(x==3) { cout<<“請輸入變壓器銘牌值n（容量Sn、短路損耗Pk、短路電壓百分比Uk%、空載損耗P0、空載電壓百分比I0%、變比k）：n”; cin>>a>>b>>c>>d>>e>>f; Transformer Tf[100]; int k=0; Tf[k].SN=a; Tf[k].Pk=b; Tf[k].Uk=c; Tf[k].P0=d; Tf[k].I0=e; Tf[k].Kk=f; data[ii][0]=Tf[k].GetP0(); data[ii][1]=Tf[k].GetQ0(); data[ii][6]=0.0; data[ii][7]=0.0; data[ii][2]=Tf[k].GetUt(data[ii-1][3]+data[ii-1][1]+data[ii][7],data[ii-1][4]+data[ii-1][0]+data[ii][6],data[ii-1][2]);//參數要改 data[ii][3]=Tf[k].GetPt(data[ii-1][3]+data[ii-1][1]+data[ii][7],data[ii-1][4]+data[ii-1][0]+data[ii][6],data[ii-1][2]); data[ii][4]=Tf[k].GetQt(data[ii-1][3]+data[ii-1][1]+data[ii][7],data[ii-1][4]+data[ii-1][0]+data[ii][6],data[ii-1][2]); data[ii][5]=Tf[k].Kk; k++; ii++; } } cout<<“輸入首端電壓值U1：n”;double U1;cin>>U1;data[ii][0]=0.0;data[ii][1]=0.0; } data[ii][2]=U1;data[ii][3]=data[ii-1][3]+data[ii-1][1]+data[ii][7];data[ii][4]=data[ii-1][4]+data[ii-1][0]+data[ii][6];data[ii][5]=1.0;data[ii][6]=0.0;data[ii][7]=0.0;int jj;jj=ii;if(fabs(data[ii][2]-data[ii-1][2])<10.0)cout<<“末端電壓值Ud=”< 南 京 理 工 大 學 《電力系統穩態分析》 課程報告 姓名 XX 學 號： 5*** 自動化學院 電氣工程 基于牛頓-拉夫遜法的潮流計算例題編程報學院(系): 專 業: 題 目: 任課教師 碩士導師 告 楊偉 XX 2015年6月10號 基于牛頓-拉夫遜法的潮流計算例題編程報告 摘要：電力系統潮流計算的目的在于：確定電力系統的運行方式、檢查系統中各元件是否過壓或者過載、為電力系統繼電保護的整定提供依據、為電力系統的穩定計算提供初值、為電力系統規劃和經濟運行提供分析的基礎。潮流計算的計算機算法包含高斯—賽德爾迭代法、牛頓-拉夫遜法和P—Q分解法等，其中牛拉法計算原理較簡單、計算過程也不復雜，而且由于人們引入泰勒級數和非線性代數方程等在算法里從而進一步提高了算法的收斂性和計算速度。同時基于MATLAB的計算機算法以雙精度類型進行數據的存儲和運算, 數據精確度高，能進行潮流計算中的各種矩陣運算，使得傳統潮流計算方法更加優化。 一 研究內容 通過一道例題來認真分析牛頓-拉夫遜法的原理和方法（采用極坐標形式的牛拉法），同時掌握潮流計算計算機算法的相關知識，能看懂并初步使用MATLAB軟件進行編程，培養自己電力系統潮流計算機算法編程能力。 例題如下：用牛頓-拉夫遜法計算下圖所示系統的潮流分布，其中系統中5為平衡節點，節點5電壓保持U=1.05為定值，其他四個節點分別為PQ節點，給定的注入功率如圖所示。計算精度要求各節點電壓修正量不大于10-6。 二 牛頓-拉夫遜法潮流計算 1 基本原理 牛頓法是取近似解x(k)之后，在這個基礎上，找到比x(k)更接近的方程的根，一步步地迭代，找到盡可能接近方程根的近似根。牛頓迭代法其最大優點是在方程f(x)=0的單根附近時誤差將呈平方減少，而且該法還可以用來求方程的重根、復根。電力系統潮流計算，一般來說，各個母線所供負荷的功率是已知的，各個節點的電壓是未知的（平衡節點外）可以根據網絡結構形成節點導納矩陣，然后由節點導納矩陣列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計算的問題就轉化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對功率平衡方程求偏導，得出對應的雅可比矩陣，給未知節點賦電壓初值，將初值帶入功率平衡方程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節點電壓不平衡量（未知的）構成了誤差方程，解誤差方程，得到節點電壓不平衡量，節點電壓加上節點電壓不平衡量構成節點電壓新的初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。2 基本步驟和設計流程圖 形成了雅克比矩陣并建立了修正方程式，運用牛頓-拉夫遜法計算潮流的核心問題已經解決，已有可能列出基本計算步驟并編制流程圖。由課本總結基本步驟如下： 1)形成節點導納矩陣Y； 2)設各節點電壓的初值，如果是直角坐標的話設電壓的實部e和虛部f；如果是極坐標的話則設電壓的幅值U和相角a； 3)將各個節點電壓的初值代入公式求修正方程中的不平衡量以及修正方程的系數矩陣的雅克比矩陣； 4)解修正方程式，求各節點電壓的變化量，即修正量； 5)計算各個節點電壓的新值，即修正后的值； 6)利用新值從第（3）步開始進入下一次迭代，直至達到精度退出循環； 7)計算平衡節點的功率和線路功率，輸出最后計算結果； ① 公式推導 ② 流程圖 三 matlab編程代碼 clear; % 如圖所示1，2，3，4為PQ節點,5為平衡節點 y=0; % 輸入原始數據，求節點導納矩陣 y(1,2)=1/(0.07+0.21j); y(4,5)=0;y(1,3)=1/(0.06+0.18j); y(1,4)=1/(0.05+0.10j); y(1,5)=1/(0.04+0.12j); y(2,3)=1/(0.05+0.10j); y(2,5)=1/(0.08+0.24j); y(3,4)=1/(0.06+0.18j); for i=1:5 for j=i:5 y(j,i)=y(i,j); end end Y=0; % 求節點導納矩陣中互導納 for i=1:5 for j=1:5 if i~=j Y(i,j)=-y(i,j); end end end % 求節點導納矩陣中自導納 for i=1:5 Y(i,i)=sum(y(i,:)); end Y % Y為導納矩陣 G=real(Y); B=imag(Y);% 輸入原始節點的給定注入功率 S(1)=0.3+0.3j; S(2)=-0.5-0.15j; S(3)=-0.6-0.25j; S(4)=-0.7-0.2j; S(5)=0; P=real(S); Q=imag(S); % 賦初值，U為節點電壓的幅值，a為節點電壓的相位角 U=ones(1,5); U(5)=1.05; a=zeros(1,5); x1=ones(8,1); x2=ones(8,1); k=0; while max(x2)>1e-6 for i=1:4 for j=1:4 H(i,j)=0; N(i,j)=0; M(i,j)=0; L(i,j)=0; oP(i)=0; oQ(i)=0; end end % 求有功、無功功率不平衡量 for i=1:4 for j=1:5 oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j))); end oP(i)=oP(i)+P(i); oQ(i)=oQ(i)+Q(i); end x2=[oP,oQ]'; % x2為不平衡量列向量 % 求雅克比矩陣 % 當i~=j時，求H,N,M,L for i=1:4 for j=1:4 if i~=j H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j))); N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); L(i,j)=H(i,j); M(i,j)=-N(i,j); end end end % 當i=j時，求H,N,M,L for i=1:4 for j=1:5 if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j))) end end N(i,i)=N(i,i)-2*(U(i))^2*G(i,i); L(i,i)=L(i,i)+2*(U(i))^2*B(i,i); end J=[H,N;M,L] % J為雅克比矩陣 x1=-((inv(J))*x2); % x1為所求△x的列向量 % 求節點電壓新值，準備下一次迭代 for i=1:4 oa(i)=x1(i); oU(i)=x1(i+4)*U(i); end for i=1:4 a(i)=a(i)+oa(i); U(i)=U(i)+oU(i); end k=k+1; end k,U,a % 求節點注入功率 i=5; for j=1:5 P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i); Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i); end S(5)=P(5)+Q(5)*sqrt(-1); S % 求節點注入電流 I=Y*U' 四 運行結果 節點導納矩陣 經過五次迭代后的雅克比矩陣 迭代次數以及節點電壓的幅值和相角（弧度數） 節點注入功率和電流 五 結果分析 在這次學習和實際操作過程里：首先，對電力系統分析中潮流計算的部分特別是潮流計算的計算機算法中的牛頓-拉夫遜法進行深入的研讀，弄明白了其原理、計算過程、公式推導以及設計流程。牛頓-拉夫遜法是求解非線性方程的迭代過程，其計算公式為?F?J?X，式中J為所求函數的雅可比矩陣；?X為需要求的修正值；?F為不平衡的列向量。利用x(*)=x(k+1)+?X(k+1)進行多次迭代，通過迭代判據得到所需要的精度值即準確值x(*)。六 結論 通過這個任務，自己在matlab編程，潮流計算，word文檔的編輯功能等方面均有提高，但也暴漏出一些問題：理論知識儲備不足，對matlab的性能和特點還不能有一個全面的把握，對word軟件也不是很熟練，相信通過以后的學習能彌補這些不足，達到一個新的層次。第二篇：電力系統潮流計算程序
第三篇：電力系統潮流計算的C語言實現
第四篇：電力系統潮流計算程序
第五篇：電力系統潮流計算
點擊咨詢 