第一篇：电力系统潮流计算程序
#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=”< 电力系统潮流计算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 南 京 理 工 大 学 《电力系统稳态分析》 课程报告 姓名 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' 四 运行结果 节点导纳矩阵 经过五次迭代后的雅克比矩阵 迭代次数以及节点电压的幅值和相角（弧度数） 节点注入功率和电流 五 结果分析 在这次学习和实际操作过程里：首先，对电力系统分析中潮流计算的部分特别是潮流计算的计算机算法中的牛顿-拉夫逊法进行深入的研读，弄明白了其原理、计算过程、公式推导以及设计流程。牛顿-拉夫逊法是求解非线性方程的迭代过程，其计算公式为FJX，式中J为所求函数的雅可比矩阵；X为需要求的修正值；F为不平衡的列向量。利用x(*)=x(k+1)+X(k+1)进行多次迭代，通过迭代判据得到所需要的精度值即准确值x(*)。六 结论 通过这个任务，自己在matlab编程，潮流计算，word文档的编辑功能等方面均有提高，但也暴漏出一些问题：理论知识储备不足，对matlab的性能和特点还不能有一个全面的把握，对word软件也不是很熟练，相信通过以后的学习能弥补这些不足，达到一个新的层次。 #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< 电力系统潮流计算程序设计 姓名：韦应顺 学号：2011021052 电力工程学院 牛顿—拉夫逊潮流计算方法具有能够将非线性方程线性化的特点，而使用MATLAB语言是由于MATLAB语言的数学逻辑强，易编译。 【】【】1.MATLAB程序12 Function tisco %这是一个电力系统潮流计算的程序 n=input（‘n请输入节点数：n=’）； m=input(‘请输入支路数：m=’);ph=input(‘n请输入平衡母线的节点号：ph=’)； B1=input（‘n请输入支路信号：B1=’）;%它以矩阵形式存贮支路的情况，每行存贮一条支路 %第一列存贮支路的一个端点 %第二列存贮支路的另一个端点 %第三列存贮支路阻抗 %第四列存贮支路的对地导纳 %第五列存贮变压器的变比，注意支路为1 %第六列存贮支路的序号 B2=input(‘n请输入节点信息：B2=’)； %第一列为电源侧的功率 %第二列为负荷侧的功率 %第三列为该点的电压值 %第四列为该点的类型：1为PQ，2为PV节点，3为平衡节点 A=input（‘n请输入节点号及对地阻抗：A=’）； ip=input（‘n请输入修正值：ip=’）； %ip为修正值);Y=zeros（n）； Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）； Y(p,q)=Y(p,q);f=zeros（1，n）； no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2； End for i=1：n G=real(Y);if A（i2）=0 B=imag(Y);p=A(i1)； Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while max(abs(De))>ip&max(abs(Dfi)>ip t=t+1; e=e+De； f=f+Df; [C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)； J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)； [De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v); jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3)); b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2; c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%该子程序是用来求取Df for i=1:n If i=ph C(i)=0; D(i)=0; For j=i:n C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If B2(i4)=2 p=2*i-1; Df(p)=P(i)-P1;p=p+1;else p=2*i-1; Df(p)=P(i)-P1;p=p+1; Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no=[];end function [De,Df]=hxf(J,Df,ph,n,no)%该子函数是为求取De Df DX=JDf;DX1=DX; x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0; For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else DX=[DX1,0,0];End k=0; [x,y]=size(DX);For i=1:2:x K=k+1; Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i); X2=G(i,j)*e(i)+B(i,j)*f(i);%该子程序是用来求取jacci矩阵 for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End 2实例验证 【例题】设有一系统网络结线见图1，各支路阻抗和各节点功率均已以标幺值标示于图1中，其中节点2连接的是发电厂，设节点1电压保持U1＝1.06定值，试计算其中的潮流分布，请输入节点数：n=5 请输入支路数：m=7 请输入平衡母线的节点号：ph=l 请输入支路信息： BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 请输入节点信息： B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 请输入节点号及对地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 请输入修正值：ip=0．000 0l 参考文献 [1]陈珩．电力系统稳定分析[M]．北京：中国电力出版社，2002：139—187． [2]郑阿奇．MATLAB实用教程[M]．北京：电子工业出版社，2005：1-243． [3] 束洪春，孙士云，等．云电送粤交商流混联系统全过 程动态电压研究[J】．中国电力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn锄ic Voltage stabil时of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4． [4] 朱新立，汤涌，等．大电网安全分析的全过程动态仿 真技术[J】．电网技术，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn锄ic simulation for the stabilhy a眦lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28． [5]Roytelm锄I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．发电机励磁系统数学模型及参 数对电网动态稳定性分析结果影响的研究[J】．继电 器，2007，35(21)：22-27． SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g锄e翰to璐baScd ∞de诅iled excitation system models柚d parameterS t0 power铲id dyn锄ic stabil时【J】．Relay，2007，35(2 1)：22-27． [7] 方思立，朱方．快速励磁系统对系统稳定的影响[J】．中 国电机工程学报，1986，6(1)：20．28． FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞ excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28． [8] 刘取．电力系统稳定性及发电机励磁控制[M】．北京： 中国电力出版社，2007． LIU Qu．Power system S诅bility锄d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro诅ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996． [10] 陈利芳，陈天禄．浅谈自并励励磁系统在大容量机组 中的应用【J】．继电器，2007，35(1)：8l培4． CHEN Li-f抽岛CHEN Tian—lIL Application of 辩l仁exci组tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84． [11] 方思立，刘增煌，孟庆和．大型汽轮发电机自并励励 磁系统的应用条件【J]．中国电力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu锄g，MENG Qin争hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63． [12]梁小冰，黄方能．利用EMTDC进行长持续时间过程 的仿真研究【J】．电网技术，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan争眦ng．How to cany out simulalion of long dul‘ation processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57． [13]王卉，陈楷，彭哲，等．数字仿真技术在电力系统中 的应用及常用的几种数字仿真工具【J】．继电器，2004，32(21)：7l一75． wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眦hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75． [14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI时studies【s】．第二篇：电力系统潮流计算程序
第三篇：电力系统潮流计算
第四篇：电力系统通用潮流计算C语言程序
第五篇：电力系统潮流计算程序设计
![下载电力系统潮流计算程序[推荐5篇]word格式文档](http://static.xiexiebang.com/skin/default/images/icon_word.png){kind=icon}
点击咨询 