c语言调用原始数据提示数据打开失败请高手帮忙
/* FUNCTION : POWER FLOW *//* WRITTEN BY : HUANG&YANG&TONG */
#include <stdio.h>
#include <math.h>
/*** 宏定义 ***/
#define eP 0.00001
#define eQ 0.00001
#define Y_(i,j) (*(*(Y+i)-i+j))
#define Yij (*(Yi+j))
#define Yji (*(*(Y+j)-j+i))
#define Pji Yji.G*cos(tmp)+Yji.B*sin(tmp)
#define Pij Yij.G*cos(tmp)+Yij.B*sin(tmp)
#define Qji Yji.G*sin(tmp)-Yji.B*cos(tmp)
#define Qij Yij.G*sin(tmp)-Yij.B*cos(tmp)
/*** 结构体变量定义 ***/
struct Nodetype /* 节点功率 */
{
float P,Q;
};
struct Linetype /* 线路类型 */
{
float G,B,B0,k;
};
/*** 子函数声明 ***/
void in_node(); /* 读节点功率 */
void in_line(); /* 读线路参数 */
void B_Form(); /* 生成BP、BQ矩阵 */
void factor(); /* 求因子表 */
void solve(float **B,float *X,int N); /* 解方程组 */
void PrtNode(); /*打印节点参数*/
void ErrorMsg(int Flag); /* 错误提示信息 */
/*** 全局变量声明 ***/
int Node; /*节点数*/
int *Num; /*保存原始节点序号*/
int NP,NQ=0; /*PV+PQ、PQ节点数*/
struct Nodetype *No; /*节点数据*/
struct Linetype **Y; /*线路参数*/
float **BP,**BQ; /*有功、无功简化雅克比矩阵B*/
float *V; /*节点电压有效值*/
float *Dlta; /*节点电压相角值*/
unsigned int count=0; /*迭代计数*/
int i,j,k,m; /*通用下标值*/
float tmp; /*临时数据暂存*/
char *Type; /*节点类型*/
FILE *in,*out; /*输入、输出文件指针*/
/****************** 主函数 ******************/
/**↓****↓****↓** 主函数 **↓****↓****↓**/
int main(void)
{
int kp=1,kq=1; /*P、Q精度判断：1-不满足，0-满足精度要求*/
float *dP,*dPi,*dQ,*dQi; /*ΔP、ΔQ*/
float Dltai;
struct Linetype *Yi;
struct Nodetype *Noi;
float tP,tQ;
if((in=fopen("Data.txt","r"))==NULL) ErrorMsg(1);
if((out=fopen("out.txt","w"))==NULL) ErrorMsg(2);
in_node(); /*读取节点参数并统计PQ、PV节点数*/
in_line(); /*读取线路参数并形成Y矩阵*/
B_Form(); /*形成B(BP&BQ)矩阵*/
factor(); /*求B因子式(仍保存于BP&BQ)*/
for(i=0;i<NQ;i++) *(V+i)=1; /*对PQ节点电压V赋初值*/
dP=(float *)malloc(sizeof(float)*NP); /*ΔP*/
dQ=dP; /*ΔQ*//*ΔP、ΔQ不同时存在，故而可共用空间*/
loop: /****开始迭代****/
if(kp==0&&kq==0) goto loopEnd;
count++; /*迭代次数加一*/
if(count==65535) ErrorMsg(99); /*不收敛，退出*/
kp=0;
for(i=0;i<NP;i++)
{
dPi=dP+i;
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
*dPi=0;
for(j=0;j<Node;j++)
{
tmp=Dltai-*(Dlta+j); /*tmp即δij*/
if(i>j) *dPi+=*(V+j)*(Pji);
else *dPi+=*(V+j)*(Pij);
} /*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
*dPi*=*(V+i);
*dPi=(*(No+i)).P-*dPi; /*求得ΔPi*/
if(fabs(*dPi)>0x8fffffff) ErrorMsg(99); /*不收敛，退出*/
if(fabs(*dPi)>eP) kp=1; /*有不满足精度的ΔP即令kp=1*/
*dPi/=*(V+i); /*求得常数项ΔPi/Vi*/
}
if(kp==0) goto loopQ;
solve(BP,dP,NP);
for(i=0;i<NP;i++) *(Dlta+i)+=(*(dP+i)/(*(V+i))); /*修正相角δ+=Δδ*/
loopQ:
if(kp==0&&kq==0) goto loopEnd;
kq=0;
for(i=0;i<NQ;i++)
{
dQi=dQ+i;
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
*dQi=0;
for(j=0;j<Node;j++)
{
tmp=Dltai-*(Dlta+j); /*tmp即δij*/
if(i>j) *dQi+=*(V+j)*(Qji);
else *dQi+=*(V+j)*(Qij);
} /*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
*dQi*=*(V+i);
*dQi=(*(No+i)).Q-*dQi; /*求得ΔQi*/
if(fabs(*dQi)>0x8fffffff) ErrorMsg(99); /*不收敛，退出*/
if(fabs(*dQi)>eQ) kq=1; /*有不满足精度的ΔQ即令kq=1*/
*dQi/=*(V+i); /*求得常数项ΔQi/Vi*/
}
if(kq==0) goto loop;
solve(BQ,dQ,NQ);
for(i=0;i<NQ;i++) *(V+i)+=*(dQ+i); /*修正PQ节点电压V+=ΔV*/
goto loop; /*无功迭代，则必定需要下一轮回迭代判断*/
loopEnd: /****迭代结束****/
free(dP); /*释放内存空间*/
/****计算PV节点和平衡节点的无功功率Q****/
for(i=NQ;i<Node;i++)
{
Noi=No+i;
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
for(j=0;j<Node;j++)
{
tmp=Dltai-*(Dlta+j); /*tmp即δij*/
if(i>j) (*Noi).Q+=*(V+j)*(Qji);
else (*Noi).Q+=*(V+j)*(Qij);
} /*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
(*Noi).Q*=*(V+i);
}
/****计算平衡节点的有功功率P****/
i=NP;
Noi=No+i;
Dltai=*(Dlta+i);
for(j=0;j<Node;j++)
{
tmp=Dltai-*(Dlta+j); /*tmp即δij*/
(*Noi).P+=*(V+j)*(Pji);
} /*注意到Y矩阵为上三角矩阵，i>j时要交换下标*/
(*Noi).P*=*(V+i);
/****输出最终结果****/
fprintf(out,"\n\n【 潮流计算结果(节点) 】 ( 迭代次数k=%3d )\n",count-1);
PrtNode();
/****计算全部线路功率****/
fprintf(out,"\n\n【 潮流计算结果(线路) 】\n");
fprintf(out," 线路 P Q\n");
for(k=0;k<Node;k++)
{
i=*(Num+k);
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
Noi=No+i;
for(m=0;m<Node;m++)
{
j=*(Num+m);
if(j==i) continue;
tmp=Dltai-*(Dlta+j); /*tmp即δij*/
if(j<i)
{
if(Yji.B==0) continue; /*若Bij=0，则节点i、j无直接联系*/
tP=*(V+j)*(Pji);
tP=*(V+i)*Yji.G-tP;
tP*=*(V+i);
tQ=-*(V+j)*(Qji);
tQ-=*(V+i)*(Yji.B-Yji.B0/Yji.k);
tQ*=*(V+i);
}
else
{
if(Yij.B==0) continue; /*若Bij=0，则节点i、j无直接联系*/
tP=*(V+j)*(Pij);
tP=*(V+i)*Yij.G-tP;
tP*=*(V+i);
tQ=-*(V+j)*(Qij);
tQ-=*(V+i)*(Yij.B-Yij.B0);
tQ*=*(V+i);
}
fprintf(out,"S[%d,%d]= (%10.6f,%10.6f)\n",k+1,m+1,-tP,-tQ);
}
}
fclose(out);
system("cmd /c start out.txt");
return(0);
}
/**↑****↑****↑** 主函数 **↑****↑****↑**/
/****************** 主函数 ******************/
/**************** 子函数：读节点数据 ****************/
void in_node()
{
struct Nodetype *Noi; /*临时变量*/
fscanf(in,"%d %d",&Node,&k);/*读取节点数Node*/
NP=Node-1; /*PV+PQ节点数，即非平衡节点数目*/
Num=(int *)malloc(sizeof(int)*Node); /*开Node个空间，每节点一个*/
V=(float *)malloc(sizeof(float)*Node); /*电压*/
Dlta=(float *)malloc(sizeof(float)*Node); /*电压相角*/
No=(struct Nodetype *)malloc(sizeof(struct Nodetype)*Node);/*节点功率*/
j=1;
while(k!=0) /*若k=0，表明节点数据读取完毕*/
{
switch(k)
{
case 1:k=NQ;NQ++;break; /*NQ统计PQ节点个数*/
case 2:k=NP-j;j++;break; /*从NP-1个空间倒着保存PV节点*/
case 3:k=NP;break; /*平衡节点*/
default:ErrorMsg(3);
}
Noi=No+k;
fscanf(in,"%d %f %f %f %f",&i,&(*Noi).P,&(*Noi).Q,V+k,Dlta+k);
i--; /*节点编号减一，以和程序表达方式兼容*/
*(Num+i)=k; /*第i个Num元素中存放i节点在No中的下标*/
fscanf(in,"%d",&k); /*读取节点类型*/
}
if(NQ+j!=Node) ErrorMsg(4); /*检验节点数据个数是否够Node个*/
fprintf(out,"【 节点参数表 】\n");
PrtNode();
fprintf(out," 总节点：%d\n PQ节点：%d\n PV节点：%d\n",Node,NQ,NP-NQ);
}
/************ 子函数：读线路数据，并形成节点导纳矩阵Y ************/
void in_line()
{
struct Linetype *Yi;
float R,X,k,B;
m=sizeof(struct Linetype);
Y=(struct Linetype **)malloc(m*Node); /*先开Node行,每一个节点一行*/
for(i=0;i<Node;i++) /*再在第i行上面开辟Node-i列*/
{ /*即以上三角存储 Y 矩阵*/
*(Y+i)=(struct Linetype *)malloc(m*(Node-i));
Yi=*(Y+i)-i;
for(j=i;j<Node;j++) {Yij.G=Yij.B=Yij.B0=Yij.k=0;} /*初始化*/
}
while(feof(in)==0) /*文件指针到文件末*/
{
fscanf(in,"%d %d %f %f %f %f",&i,&j,&R,&X,&k,&B);
i--;j--;
i=*(Num+i); /*转换节点号为该节点在程序中的储存编号*/
j=*(Num+j);
(*(*(Y+i))).B+=B; /* 将对地充电导纳累加到自导纳 */
(*(*(Y+j))).B+=B;
if(k!=1.0)
{
X*=k;R=0;
tmp=(1-k)/X;
(*(*(Y+i))).B+=tmp; /*将变压器的对地充电容纳累加到自导纳*/
(*(*(Y+j))).B+=-(tmp/k);
B=tmp;
k=-k;
}
if(i>j) {tmp=i;i=j;j=tmp;k=1/k;B*=k;}
Yi=*(Y+i)-i; /* 以Yi代替*(Y+i)-i，简化表达式并避免重复计算 */
Yij.B0=B; /*保存ij0、ji0对地充电电容到Bij0*/
Yij.k=k; /*且有B0ji=B0ij/k*/
tmp=R*R+X*X;
R/=tmp;
X/=tmp;
Yij.G=-R; /* 生成互导纳 */
Yij.B=X;
(*(*(Y+i))).G+=R; /* 将线路互导纳累加到自导纳 */
(*(*(Y+i))).B+=-X;
(*(*(Y+j))).G+=R;
(*(*(Y+j))).B+=-X;
}
fclose(in);
fprintf(out,"\n【 节点导纳矩阵 Y 】\n");
for(k=0;k<Node;k++)
{
i=k;
i=*(Num+i); /*查取第i节点在程序中存储序号*/
for(j=0;j<k;j++) fprintf(out,"\t\t\t");
for(m=k;m<Node;m++)
{
j=*(Num+m); /*查取第j节点在程序中存储序号*/
if(i<j) fprintf(out,"(%10.6f,%10.6f) ",Y_(i,j).G,Y_(i,j).B);
else fprintf(out,"(%10.6f,%10.6f) ",Y_(j,i).G,Y_(j,i).B);
}
fprintf(out,"\n");
}
}
/**************** 子函数：生成BP、BQ矩阵 ****************/
void B_Form()
{
float *BPi,*BQi;
struct Linetype *Yi;
int size=sizeof(float);
BP=(float **)malloc(size*NP); /*以上三角存储*/
for(i=0;i<NP;i++) *(BP+i)=(float *)malloc(size*(NP-i));
for(i=0;i<NP;i++)
{
BPi=*(BP+i)-i; /* 以BPi代替*(BP+i)-i，避免重复计算 */
Yi=*(Y+i)-i;
for(j=i;j<NP;j++) *(BPi+j)=Yij.B; /*(BPi+j)即相当于BP[i][j]*/
}
BQ=BP; /* BP包含BQ，BP左上角的NQ*NQ子阵即BQ */
}
/**************** 子函数：求因子表 ****************/
void factor()
{
float *BPi,*BPk,*BQi;
for(i=0;i<NP;i++)
{
BPi=*(BP+i)-i;
for(k=0;k<i;k++)
{
BPk=*(BP+k)-k;
tmp=(*(BPk+i))/(*(BPk+k));
for(j=i;j<NP;j++) (*(BPi+j))-=tmp*(*(BPk+j));
}
*(BPi+i)=1/(*(BPi+i));
for(j=i+1;j<NP;j++) *(BPi+j)*=*(BPi+i);
}
}
/**************** 子函数：解方程组 ****************/
void solve(float **B,float *X,int N)
{
float *Bi,*Xi;
for(i=0;i<N;i++) *(X+i)=-*(X+i);/*对常数项取负*/
/****对常数列进行前代****/
for(i=0;i<N;i++)
{
Bi=*(B+i)-i;
Xi=X+i;
for(j=i+1;j<N;j++) *(X+j)-=*(Bi+j)**Xi;
*Xi*=*(Bi+i);
}
/****回代以求解方程组****/
for(i=N-1;i>=0;i--)
{
Bi=*(B+i)-i;
Xi=X+i;
for(j=N-1;j>i;j--) *Xi-=*(Bi+j)**(X+j);
}
}
/****打印节点参数****/
void PrtNode()
{
struct Nodetype *Noi;
fprintf(out,"节点 类型 P Q V δ\n");
for(i=0;i<Node;)
{
j=*(Num+i); /*查取第i节点在程序中存储序号*/
Noi=No+j;
if(j<NQ) Type="PQ"; else Type="PV";
if(j==NP) Type="BS";
fprintf(out,"%3d %s %10.6f %10.6f %10.6f %10.6f\n",
++i,Type,(*Noi).P,(*Noi).Q,*(V+j),*(Dlta+j)/0.017453);
}
}
/**************** 子函数：错误信息 ****************/
void ErrorMsg(int Flag)
{
switch(Flag)
{
case 1:printf("\n\tError(1): Failed to Open File \"Data.txt\"!");break;
case 2:printf("\n\tError(2): Failed to Creat File \"out.txt\"!");break;
case 3:printf("\n\tError(3): Node Data Error, Please Check!");break;
case 4:printf("\n\tError(4): Lack Node Data, Please Check!");break;
case 99:printf("k=%d\n\tError(99): It's Emanative!",count);break;
}
getch();
fclose(out);
exit(Flag);
}
