华电-通用潮流上机程序.doc

#include "stdafx.h"
#include <iostream>
#include <fstream>
#include "NEquation.h"
#include "math.h"
#include "config.h"
using namespace std;

void GetData(int& nBus, int& nBranch)//Read the data
{
FILE *fp;
int i;

if((fopen_s(&fp,"data\\data5.txt","r"))!=0)
{
printf("Can not open the file named 'data.txt' \n");
return;
}
fscanf_s(fp,"%d,%d",&nBus,&nBranch,100);
for(i=0;i<nBus;i++)
{
fscanf_s(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&gBus[i].No,&gBus[i].Voltage,&gBus[i].Phase,
&gBus[i].GenP,&gBus[i].GenQ,&gBus[i].LoadP,&gBus[i].LoadQ,&gBus[i].Type,100);
}

for(i=0;i<nBranch;i++)
{
fscanf_s(fp,"%d,%d,%d,%f,%f,%f,%f",&gBranch[i].No,&gBranch[i].No_I,&gBranch[i].No_J,
&gBranch[i].R,&gBranch[i].X,&gBranch[i].B,&gBranch[i].k,100);
}

fclose(fp);
}

void GetYMatrix(int& nBus, int& nBranch)
{
int i,j,bus1,bus2;
float r,x,d,g,b,k;
FILE *fp;

for(i=0;i<nBus;i++)
{
for(j=0;j<nBus;j++)
{
gY_G[i][j]=0;
gY_B[i][j]=0;
}
}

for(i=0; i<nBranch; i++)
{
bus1=gBranch[i].No_I-1;
bus2=gBranch[i].No_J-1;
r=gBranch[i].R;
x=gBranch[i].X;
d=r*r+x*x;
g=r/d;
b=-x/d;
if(gBranch[i].k==0)
{
if(bus1 != bus2)
{
gY_G[bus1][bus1]=gY_G[bus1][bus1]+g;
gY_G[bus2][bus2]=gY_G[bus2][bus2]+g;
gY_G[bus1][bus2]=gY_G[bus1][bus2]-g;
gY_G[bus2][bus1]=gY_G[bus2][bus1]-g;
gY_B[bus1][bus1]=gY_B[bus1][bus1]+b+gBranch[i].B;
gY_B[bus2][bus2]=gY_B[bus2][bus2]+b+gBranch[i].B;
gY_B[bus1][bus2]=gY_B[bus1][bus2]-b;
gY_B[bus2][bus1]=gY_B[bus2][bus1]-b;
}
else
{
break;
}

}
else
{
// add your codes
}
}

// output the Y matrix
if((fopen_s(&fp,"data\\ymatrix1.txt","w"))!=0)
{
printf("Can not open the file named 'ymatrix.txt' \n");
return ;
}

fprintf(fp,"---Y Matrix---\n");
for(i=0;i<nBus;i++)
{
for(j=0;j<nBus;j++)
fprintf_s(fp,"Y(%d,%d)=(%10.5f,%10.5f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);

}
fclose(fp);
}

void SetInitial(int& nBus)
{
int i;

for(i=0;i<nBus;i++)
{
if(gBus[i].Type==2)
{
gf[i]=gBus[i].Voltage*sin(gBus[i].Phase);
ge[i]=gBus[i].Voltage*cos(gBus[i].Phase);
}
else
{
gf[i]=0;
ge[i]=1;
}
}
}

void GetUnbalance(int& nBus)

{ int i,j;
float sum1=0,sum2=0;
FILE *fp;

for(i=0;i<nBus;i++)
{
if (gBus[i].Type!=2)
{
sum1=0;
sum2=0;
for(j=0;j<nBus;j++)
{
sum1=sum1+ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]); //chuqu balance point add your codes
sum2=sum2+gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
}
gDelta_P[i]= gBus[i].GenP-gBus[i].LoadP-sum1;
gDelta_Q[i]= gBus[i].GenQ-gBus[i].LoadQ-sum2;

if (gBus[i].Type==1)
gDelta_Q[i]= gBus[i].Voltage* gBus[i].Voltage-(gf[i]*gf[i]+ge[i]*ge[i]);
}
}

for (i=0;i<nBus;i++)
{
j=2*i;
gDelta_PQ[j]=gDelta_P[i];
gDelta_PQ[j+1]=gDelta_Q[i];
}

}

void GetJaccobi(int& nBus)
{
int i,j,k;
float s1=0,s2=0;
FILE *fp;

for (i=0;i<2*nBus;i=i+2) // add your codes

for(j=0;j<2*nBus;j=j+2)
if (i/2==j/2) // i is or not=j.
{ s1=0;
s2=0;
for (k=0;k<nBus;k++)
{
s1=s1+gY_G[i/2][k]*gf[k]+gY_B[i/2][k]*ge[k];
s2=s2+gY_G[i/2][k]*ge[k]-gY_B[i/2][k]*gf[k];
}
gJaccobi[i][j]=gY_G[i/2][j/2]*gf[j/2]-gY_B[i/2][j/2]*ge[j/2]+s1;
gJaccobi[i][j+1]=gY_G[i/2][j/2]*ge[j/2]+gY_B[i/2][j/2]*gf[j/2]+s2;

if (gBus[i/2].Type==0) //PQ pioint and PV Point difference.
{
gJaccobi[i+1][j]=-(gY_G[i/2][j/2]*ge[j/2]+gY_B[i/2][j/2]*gf[j/2])+s2;
gJaccobi[i+1][j+1]=gY_G[i/2][j/2]*gf[j/2]-gY_B[i/2][j/2]*ge[j/2]-s1;
}

else //pv bus.
{
gJaccobi[i+1][j]=2*gf[i/2];
gJaccobi[i+1][j+1]=2*ge[i/2];
}
}
else
{
gJaccobi[i][j]=gY_G[i/2][j/2]*gf[i/2]-gY_B[i/2][j/2]*ge[i/2];
gJaccobi[i][j+1]=gY_G[i/2][j/2]*ge[i/2]+gY_B[i/2][j/2]*gf[i/2];
if (gBus[i/2].Type==0)
{
gJaccobi[i+1][j]=-gJaccobi[i][j+1];
gJaccobi[i+1][j+1]=gJaccobi[i][j];
}

else
{
gJaccobi[i+1][j]=0;
gJaccobi[i+1][j+1]=0;
}

}

// output Jaccobimatrix.
if((fopen_s(&fp,"data\\Jaccobimatrix.txt","w"))!=0)
{
printf("Can not open the file named 'Jaccobimatrix.txt' \n");
return ;
}

fprintf(fp,"---JaccobiMatrix---\n");

for(i=0;i<2*nBus;i++)
{
fprintf_s(fp,"%d ,",i+1);

for(j=0;j<2*nBus;j++)
{
fprintf_s(fp,"[%10.5f]",gJaccobi[i][j]);
}
fprintf_s(fp,"\n");

}
fclose(fp);

}

void GetRevised(int& nBus)
{
int i,j;
NEquation ob1;
FILE *fp;

ob1.SetSize(2*nBus);

for(i=0;i<2*nBus;i++)
{
for(j=0;j<2*nBus;j++)
ob1.Data(i,j)=gJaccobi[i][j];

ob1.Value(i)=gDelta_PQ[i];
}

ob1.Run();

for(i=0;i<2*nBus;i++)
{
gDelta_fe[i]=ob1.Value(i);
}

for (i=0;i<nBus;i++)
{
gDelta_f[i]=gDelta_fe[2*i];
gDelta_e[i]=gDelta_fe[2*i+1];
}

// output delta_f and delta_e.

if((fopen_s(&fp,"data\\Delta_f and Delta_e matrix.txt","w"))!=0)
{
printf("Can not open the file named 'Delta_f and Delta_e matrix.txt' \n");
return ;
}

fprintf(fp,"---Delta_f and Delta_e Matrix---\n");

for(i=0;i<2*nBus;i++)
{
fprintf(fp," Delta %d=%10.8f \n",i+1,gDelta_fe[i]);

}
fprintf_s(fp,"\n");
fclose(fp);
}

void GetNewValue(int& nBus)
{
int i;
FILE *fp;
for (i=0;i<nBus;i++)
if(gBus[i].Type!=2)
{
gf[i]=gf[i]+gDelta_f[i];
ge[i]=ge[i]+gDelta_e[i]; }
if((fopen_s(&fp,"data\\gf and ge matrix.txt","w"))!=0)
{
printf("Can not open the file named 'gf and ge matrix.txt' \n");
return ;
}

fprintf(fp,"---gf and ge Matrix---\n");

for(i=0;i<nBus;i++)
{
fprintf(fp," gf[%d]=%10.6f, ge[%d]=%10.6f\n",i+1,gf[i],i+1,ge[i]);

}
fprintf_s(fp,"\n");
fclose(fp);

}

int main(int argc, char* argv[])
{
int i,Count_Num,j,n,nBus,nBranch;
float maxValue,s1,s2,k,su1,su2,deltalinep[12][12]={0,0},deltalineq[12][12]={0,0},deltap[12][12]={0,0},deltaq[12][12]={0,0};
FILE *fp;
float g,b,d,r,x;
GetData(nBus,nBranch);
GetYMatrix(nBus,nBranch);
SetInitial(nBus);

for(Count_Num=0;Count_Num<=20;Count_Num++)
{
GetUnbalance(nBus);
GetJaccobi(nBus);
GetRevised(nBus);
GetNewValue(nBus);

maxValue=fabs(gDelta_fe[0]);
for(i=1;i<=2*(nBus-1)-1;i++)
{
if(maxValue<fabs(gDelta_fe[i]))
{
maxValue=fabs(gDelta_fe[i]);
}
}
if(maxValue<Precision)
{
break;
}
}

printf("countnumber %d\n",Count_Num);
for(i=0;i<nBus;i++)
{
printf("bus(%d).voltage=%10.6f,angle=%.6f , e=%.6f,f=%.6f\n",i+1,sqrt(ge[i]*ge[i]+gf[i]*gf[i]),atan(gf[i]/ge[i])/3.1415926*180,ge[i],gf[i]);
}

//output zhuru power.
for (i=0;i<nBus;i++)
{
if(gBus[i].Type!=0)
{
su1=0;
su2=0;
for(j=0;j<nBus;j++)
{
su1=su1+ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
su2=su2+gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
}
gBus[i].GenP=su1; gBus[i].GenQ=su2;
printf("zhurup[%d]=%10.7f zhuruq[%d]=%10.7f \n",i+1,gBus[i].GenP,i+1,gBus[i].GenQ);
}
}

for(n=0; n<nBranch; n++)
{
i=gBranch[n].No_I-1;
j=gBranch[n].No_J-1;
r=gBranch[n].R;
x=gBranch[n].X;
d=r*r+x*x;
g=r/d;
b=-x/d;
if(gBranch[n].k==0)
{
if(i!=j)
{ deltap[i][j]= (ge[i]*(ge[i]-ge[j])+gf[i]*(gf[i]-gf[j]))*g+(ge[i]*gf[j]-ge[j]*gf[i])*b;

deltaq[i][j]=-(ge[i]*ge[i]+gf[i]*gf[i])*gBranch[n].B+(ge[i]*gf[j]-ge[j]*gf[i])*g-(ge[i]*(ge[i]-ge[j])+gf[i]*(gf[i]-gf[j]))*b;

deltap[j][i]=(ge[j]*(ge[j]-ge[i])+gf[j]*(gf[j]-gf[i]))*g+(ge[j]*gf[i]-ge[i]*gf[j])*b;

deltaq[j][i]=-(ge[j]*ge[j]+gf[j]*gf[j])*gBranch[n].B+(ge[j]*gf[i]-ge[i]*gf[j])*g-(ge[j]*(ge[j]-ge[i])+gf[j]*(gf[j]-gf[i]))*b;

}
else break;
}
else
{
k=gBranch[n].k;
deltap[i][j]= (ge[i]*ge[i]+gf[i]*gf[i])*k*(k-1)*g+(ge[i]*(ge[i]-ge[j])+gf[i]*(gf[i]-gf[j]))*g*k+(ge[i]*gf[j]-ge[j]*gf[i])*b*k;

deltaq[i][j]=-(ge[i]*ge[i]+gf[i]*gf[i])*k*(k-1)*b+(ge[i]*gf[j]-ge[j]*gf[i])*g*k-(ge[i]*(ge[i]-ge[j])+gf[i]*(gf[i]-gf[j]))*b*k;

deltap[j][i]=(ge[j]*ge[j]+gf[j]*gf[j])*(1-k)*g+(ge[j]*(ge[j]-ge[i])+gf[j]*(gf[j]-gf[i]))*g*k+(ge[j]*gf[i]-ge[i]*gf[j])*b*k;

deltaq[j][i]=-(ge[j]*ge[j]+gf[j]*gf[j])*(1-k)*b+(ge[j]*gf[i]-ge[i]*gf[j])*g*k-(ge[j]*(ge[j]-ge[i])+gf[j]*(gf[j]-gf[i]))*b*k;
}

}

for (i=0;i<nBus-1;i++)
for(j=i+1;j<nBus;j++)
{
deltalinep[i][j]=deltap[i][j]+deltap[j][i];

deltalineq[i][j]=deltaq[i][j]+deltaq[j][i];

printf("deltalinep[%d][%d]= %10.7f, deltalineq[%d][%d]=%10.7f \n",i+1,j+1,deltalinep[i][j],i+1,j+1,deltalineq[i][j]);
}

s1=0;s2=0;
for(i=0;i<nBus;i++)
for(j=0;j<nBus;j++)
{
s1=s1+deltalinep[i][j];
s2=s2+deltalineq[i][j];
}

printf("totallineP=%10.7f, totallineQ=%10.7f\n",s1,s2);

if((fopen_s(&fp,"data\\bus voltage and zhuru power.txt","w"))!=0)
{
printf("Can not open the file named 'bus voltage and zhuru power.txt' \n");
return 0;
}

fprintf(fp,"---bus voltage and zhuru power---\n");

for(i=0;i<nBus;i++)
{
fprintf_s(fp,"zhurup[%d]=%10.7f, zhuruq[%d]=%10.7f \n",i+1,gBus[i].GenP-gBus[i].LoadP,i+1,gBus[i].GenQ-gBus[i].LoadQ);

}
fprintf_s(fp,"\n");
for(i=0;i<nBus;i++)
{
fprintf_s(fp,"bus(%d).voltage=%10.5f,,angle=%10.5f\n",i+1,sqrt(ge[i]*ge[i]+gf[i]*gf[i]),atan(gf[i]/ge[i]));

}
fprintf_s(fp,"\n");

for(i=0;i<nBus;i++)
for(j=i+1;j<nBus;j++)
fprintf_s(fp,"deltalinep[%d][%d]= %10.5f, deltalineq[%d][%d]=%10.5f \n",i+1,j+1,deltalinep[i][j],i+1,j+1,deltalineq[i][j]);

fprintf_s(fp,"\n");
fclose(fp);

return (0);
}

4,4
1,1.00,0.00,0.00,0.00,0.30,0.18,0
2,1.00,0.00,0.00,0.00,0.55,0.13,0
3,1.10,0.00,0.50,0.00,0.00,0.00,1
4,1.05,0.00,0.00,0.00,0.00,0.00,2
1,1,2,0.10,0.41,0.01528,0.0
2,1,3,0.00,0.30,0.00000,1.1
3,1,4,0.12,0.50,0.01920,0.0
4,2,4,0.08,0.40,0.01413,0.0