牛拉法计算C程序

int N; int Type; double e; double f; double Pd; double Qd; double Ps; double Qs; double Bc; };

struct BranchType {

int Nbr; int Nl; int Nr; double R; double X; double Bn; double Kt; };

int n; int nPQ; int nPV; int nPH; int nbr; int ng; int Mark=0; double **G; double **B; double *dS;

double *mid1,*mid2; double *Us; double error=1;

double iteration=0.000001; double **Jacob;

double **invJac; double *dfe;

struct NodeType *Node; struct BranchType *Branch;

void main() {

void LoadData(); void FormY(); void DeltaS(); void FormJacob(); void InvJac(); void UpdateU(); void CalculatePQ(); int kk;

LoadData(); FormY();

printf(\"误差精度 iteration=%lf\\n\ kk=0; DeltaS();

while(error>iteration&&kk<50) {

FormJacob(); UpdateU();

DeltaS(); kk++; }

printf(\"迭代次数为%4d\\n\ CalculatePQ();

printf(\"误差值 error=%e\\n\ printf(\"\\n\"); printf(\"\\n\"); printf(\" {注:\\nN为节点号，Tp为节点类型（其中1为PQ节点，2为PV节点，3为平衡节点），\\nAmp为电压大小，Dlta为相位角，Pd为节点输出的有功功率，Qd为节点输出的无功功率，\\nPs为发电机输出有功功率，Qs为发电机输出无功功率，Bc为并联电容的电抗值}\\n\"); }

void LoadData() {

int i,j;

int tN,tType;

double te,tf,tPd,tQd,tPs,tQs,tBc; FILE *fp;

char filename[50]={\" \

printf(\"请输入数据文件名(提示：数据所在TXT文件)：\"); scanf(\"%s\

if((fp=fopen(filename,\"r\"))==NULL) {

printf(\"cannot open the file:data.txt\\n\"); return; }

fscanf(fp,\"%d\

printf(\"节点个数为：%d\\n\

Node=(struct NodeType *)malloc(sizeof(struct NodeType)*n);

printf(\"调整前的节点参数为：\\n\"); for(i=0;ifscanf(fp,\"%d%d%lf%lf%lf%lf%lf%lf%lf\i].f,&Node[i].Pd,&Node[i].Qd,&Node[i].Ps,&Node[i].Qs,&Node[i].Bc);

for(i=0;iif(Node[i].Type==1) nPQ++;

else if(Node[i].Type==2) nPV++;

else if(Node[i].Type==3) nPH++; }

printf(\"PQ节点个数：%d\\n\ printf(\"PV节点个数：%d\\n\ printf(\"平衡节点个数：%d\\n\ for(j=0;jif(Node[i].Type>Node[i+1].Type) {

tN=Node[i].N;Node[i].N=Node[i+1].N;Node[i+1].N=tN;

tType=Node[i].Type;Node[i].Type=Node[i+1].Type;Node[i+1].Type=tType; te=Node[i].e;Node[i].e=Node[i+1].e;Node[i+1].e=te; tf=Node[i].f;Node[i].f=Node[i+1].f;Node[i+1].f=tf;

tPd=Node[i].Pd;Node[i].Pd=Node[i+1].Pd;Node[i+1].Pd=tPd; tQd=Node[i].Qd;Node[i].Qd=Node[i+1].Qd;Node[i+1].Qd=tQd; tPs=Node[i].Ps;Node[i].Ps=Node[i+1].Ps,Node[i+1].Ps=tPs;

tQs=Node[i].Qs;Node[i].Qs=Node[i+1].Qs;Node[i+1].Qs=tQs; tBc=Node[i].Bc;Node[i].Bc=Node[i+1].Bc;Node[i+1].Bc=tBc; } }

Us=(double *)malloc(sizeof(double)*(n-1)); for(i=0;iprintf(\"支路个数为：%d\\n\

Branch=(struct BranchType *)malloc(sizeof(struct BranchType)*nbr);

for(i=0;ifscanf(fp,\"%d%d%d%lf%lf%lf%lf\h[i].R,&Branch[i].X,&Branch[i].Bn,&Branch[i].Kt); for(i=0;iMark=0;

for(j=0;jif(Branch[i].Nl==Node[j].N&&Mark==0) {

Branch[i].Nl=j+1; Mark=1; } } }

for(i=0;iMark=0;

for(j=0;jif(Branch[i].Nr==Node[j].N&&Mark==0) {

Branch[i].Nr=j+1; Mark=1; } } }

fclose(fp); }

void FormY() {

int i,j; double Z2;

G=(double **)malloc(sizeof(double *)*n); B=(double **)malloc(sizeof(double *)*n); for(i=0;iG[i]=(double *)malloc(sizeof(double)*n); B[i]=(double *)malloc(sizeof(double)*n); }

for(i=0;iG[i][j]=0; B[i][j]=0; }

for(i=0;iZ2=Branch[i].R*Branch[i].R+Branch[i].X*Branch[i].X; if(Branch[i].Kt==0) {

G[Branch[i].Nl-1][Branch[i].Nr-1]-=Branch[i].R/Z2; B[Branch[i].Nl-1][Branch[i].Nr-1]+=Branch[i].X/Z2;

G[Branch[i].Nr-1][Branch[i].Nl-1]=G[Branch[i].Nl-1][Branch[i].Nr-1]; B[Branch[i].Nr-1][Branch[i].Nl-1]=B[Branch[i].Nl-1][Branch[i].Nr-1]; } else {

G[Branch[i].Nl-1][Branch[i].Nr-1]-=Branch[i].R/Z2/Branch[i].Kt; B[Branch[i].Nl-1][Branch[i].Nr-1]+=Branch[i].X/Z2/Branch[i].Kt;

G[Branch[i].Nr-1][Branch[i].Nl-1]=G[Branch[i].Nl-1][Branch[i].Nr-1]; B[Branch[i].Nr-1][Branch[i].Nl-1]=B[Branch[i].Nl-1][Branch[i].Nr-1]; } }

for(i=0;iZ2=Branch[j].R*Branch[j].R+Branch[j].X*Branch[j].X;

if(Branch[j].Kt==0&&(Branch[j].Nl-1==i||Branch[j].Nr-1==i))

{

G[i][i]=G[i][i]+Branch[j].R/Z2; B[i][i]=B[i][i]-Branch[j].X/Z2; }

else if(Branch[j].Kt!=0&&(Branch[j].Nl-1==i||Branch[j].Nr-1==i)) {

G[i][i]=G[i][i]+Branch[j].R/Z2/Branch[j].Kt; B[i][i]=B[i][i]-Branch[j].X/Z2/Branch[j].Kt; } }

for(i=0;iif(Branch[i].Kt==0) {

B[Branch[i].Nl-1][Branch[i].Nl-1]+=Branch[i].Bn; B[Branch[i].Nr-1][Branch[i].Nr-1]+=Branch[i].Bn; } else {

B[Branch[i].Nl-1][Branch[i].Nl-1]-=(1-Branch[i].Kt)/Branch[i].Kt/Branch[i].Kt/Branch[i].X;

B[Branch[i].Nr-1][Branch[i].Nr-1]-=(Branch[i].Kt-1)/Branch[i].Kt/Branch[i].X; } }

for(i=0;iB[i][i]=B[i][i]+Node[i].Bc; }

void DeltaS() {

int i,j;

mid1=(double *)malloc(sizeof(double)*n); mid2=(double *)malloc(sizeof(double)*n); dS=(double *)malloc(sizeof(double)*2*(n-1));

for(i=0;imid1[i]=0; mid2[i]=0;

for(j=0;j{

mid1[i]=mid1[i]+G[i][j]*Node[j].e-B[i][j]*Node[j].f; mid2[i]=mid2[i]+G[i][j]*Node[j].f+B[i][j]*Node[j].e; }

dS[2*i]=Node[i].Ps-Node[i].Pd-(Node[i].e*mid1[i]+Node[i].f*mid2[i]); if(idS[2*i+1]=Node[i].Qs-Node[i].Qd-(Node[i].f*mid1[i]-Node[i].e*mid2[i]); else

dS[2*i+1]=Us[i]*Us[i]-(Node[i].e*Node[i].e+Node[i].f*Node[i].f); }

error=0;

for(i=0;i<2*(n-1);i++) {

if(dS[i]<0&&error<-dS[i]) error=-dS[i];

else if(dS[i]>0&&errorvoid FormJacob() {

int i,j;

Jacob=(double **)malloc(sizeof(double *)*2*(n-1)); for(i=0;i<2*(n-1);i++)

Jacob[i]=(double *)malloc(sizeof(double)*2*(n-1));

for(i=0;i<2*(n-1);i++) for(j=0;j<2*(n-1);j++) Jacob[i][j]=0; for(j=0;jfor(i=0;iif(i!=j) {

Jacob[2*i][2*j]=B[i][j]*Node[i].e-G[i][j]*Node[i].f; Jacob[2*i][2*j+1]=-G[i][j]*Node[i].e-B[i][j]*Node[i].f; } else {

Jacob[2*i][2*i]=B[i][i]*Node[i].e-G[i][i]*Node[i].f-mid2[i]; Jacob[2*i][2*i+1]=-G[i][i]*Node[i].e-B[i][i]*Node[i].f-mid1[i]; } }

for(i=0;iif(i!=j) {

Jacob[2*i+1][2*j]=G[i][j]*Node[i].e+B[i][j]*Node[i].f; Jacob[2*i+1][2*j+1]=B[i][j]*Node[i].e-G[i][j]*Node[i].f; } else {

Jacob[2*i+1][2*i]=G[i][i]*Node[i].e+B[i][i]*Node[i].f-mid1[i]; Jacob[2*i+1][2*i+1]=B[i][i]*Node[i].e-G[i][i]*Node[i].f+mid2[i]; } }

for(i=nPQ;iif(i==j) {

Jacob[2*i+1][2*i]=-2*Node[i].f; Jacob[2*i+1][2*i+1]=-2*Node[i].e; } } } }

void InvJac() {

int i,j,k; double temp;

invJac=(double **)malloc(sizeof(double *)*2*(n-1)); for(i=0;i<2*(n-1);i++)

invJac[i]=(double *)malloc(sizeof(double)*2*(n-1));

for(i=0;i<2*(n-1);i++) for(j=0;j<2*(n-1);j++) {

if(i!=j)

invJac[i][j]=0;

else

invJac[i][j]=1; }

for(i=0;i<2*(n-1);i++) {

for(j=0;j<2*(n-1);j++) {

if(i!=j) {

temp=Jacob[j][i]/Jacob[i][i]; for(k=0;k<2*(n-1);k++) {

Jacob[j][k]-=Jacob[i][k]*temp; invJac[j][k]-=invJac[i][k]*temp; } } } }

for(i=0;i<2*(n-1);i++) if(Jacob[i][i]!=1) {

temp=Jacob[i][i];

for(j=0;j<2*(n-1);j++)

invJac[i][j]=invJac[i][j]/temp; } }

void UpdateU() {

void InvJac(); int i,j;

dfe=(double *)malloc(sizeof(double)2*(n-1)); InvJac();

for(i=0;i<2*(n-1);i++) {

dfe[i]=0;

for(j=0;j<2*(n-1);j++)

dfe[i]-=invJac[i][j]*dS[j]; }

for(i=0;iNode[i].e+=dfe[2*i+1]; Node[i].f+=dfe[2*i];

} }

void CalculatePQ() {

int i,j;

int tN,tType;

double te,tf,tPd,tQd,tPs,tQs,tBc;

mid1[n-1]=0; mid2[n-1]=0; for(j=0;jmid1[n-1]=mid1[n-1]+G[n-1][j]*Node[j].e-B[n-1][j]*Node[j].f; mid2[n-1]=mid2[n-1]+G[n-1][j]*Node[j].f+B[n-1][j]*Node[j].e; }

Node[n-1].Ps=Node[n-1].e*mid1[n-1]; Node[n-1].Qs=-Node[n-1].e*mid2[n-1];

for(i=nPQ;iNode[i].Qs=Node[i].f*mid1[i]-Node[i].e*mid2[i];

for(j=0;jif(Node[i].N>Node[i+1].N) {

tN=Node[i].N;Node[i].N=Node[i+1].N;Node[i+1].N=tN;

tType=Node[i].Type;Node[i].Type=Node[i+1].Type;Node[i+1].Type=tType; te=Node[i].e;Node[i].e=Node[i+1].e;Node[i+1].e=te; tf=Node[i].f;Node[i].f=Node[i+1].f;Node[i+1].f=tf;

tPd=Node[i].Pd;Node[i].Pd=Node[i+1].Pd;Node[i+1].Pd=tPd; tQd=Node[i].Qd;Node[i].Qd=Node[i+1].Qd;Node[i+1].Qd=tQd; tPs=Node[i].Ps;Node[i].Ps=Node[i+1].Ps,Node[i+1].Ps=tPs; tQs=Node[i].Qs;Node[i].Qs=Node[i+1].Qs;Node[i+1].Qs=tQs; tBc=Node[i].Bc;Node[i].Bc=Node[i+1].Bc;Node[i+1].Bc=tBc; } }

for(i=0;ite=sqrt(Node[i].e*Node[i].e+Node[i].f*Node[i].f); tf=atan(Node[i].f/Node[i].e)*180/PI;

Node[i].e=te; Node[i].f=tf; }

printf(\"最终结果为：\\n\");

printf(\" N Tp Amp Dlta Pd Qd Ps Qs Bc\\n\");

for(i=0;iprintf(\"%4d%4d%9.4lf%9.4lf%9.4lf%9.4lf%9.4lf%9.4lf%9.4lf\\n\pe,Node[i].e,Node[i].f,Node[i].Pd,Node[i].Qd,Node[i].Ps,Node[i].Qs,Node[i].Bc); }