📄 nl_lf.c
字号:
/******************************************************************************/
/* */
/* 电 力 系 潮 流 计 算 程 序 */
/* */
/* 牛 顿 - 拉 夫 逊 法 */
/* */
/* 极 坐 标 形 式 */
/* */
/* 输入文件 input.txt */
/* 校核文件 check.txt */
/* 输出文件 output.txt */
/* */
/******************************************************************************/
/* 预编译 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* 宏定义 */
#define node_max 700
/* 定义全局变量 */
struct data
{
int i;
int j;
double a;
double b;
double c;
};
struct data1
{
int i;
double a;
double b;
};
/* num_node 节点数 */
/* num_line 线路数 */
/* num_tran 变压器数 */
/* num_bran 支路数 */
/* num_gene 发电机数 */
/* num_load 负荷数 */
int num_node,num_line,num_tran,num_bran,num_gene,num_load;
int iter_max=10; // 最大允许迭代次数
int conv,k;
double error_max=1.0e-6; //容许最大迭代误差
double **yr,**yi,*vm,*va;
double *ppq_node,*ppq_load;
double *pq_node[node_max],*pq_load[node_max];
struct data *line,*tran,*gene;
struct data1 *load;
FILE *fpi,*fpc,*fpo;
/* 申明子函数 */
void matrix_adm(); // 形成导纳矩阵
void ntlf(); // 牛顿-拉夫逊潮流迭代
void eliminate(int,double *,double *,int,int *,int *,int *,double *); // 消去运算
void zero(double *,int,int); //数组清零
double **TwoArrayAlloc_D(int,int); //定义两维整型数组
void TwoArrayFree_D(double **); //释放两维整型数组
/* 主函数 */
void main()
{
int i,j,ii,jj;
struct data *p;
struct data1 *p1;
double deg=1.8e2/3.1415927; // 用于弧度化度
double pij,qij,pji,qji; // 支路两端有功和无功
double r,x,b,v1,v2,a1; // 支路参数和两端电压几相位差
printf("\n * * * 开 始 运 行 程 序 * * *\n\n");
/* 打开文件 */
if((fpi=fopen("input.txt","r"))==NULL)
{
printf("Can not open input file!\n");
exit(0);
}
if((fpc=fopen("check.txt","w"))==NULL)
{
printf("Can not open check file!\n");
exit(0);
}
if((fpo=fopen("output.txt","w"))==NULL)
{
printf("Can not open output file!\n");
exit(0);
}
/* 从输入文件中读入原始信息 */
fscanf(fpi,"%d,%d,%d,%d,%d;",&num_node,&num_line,&num_tran,&num_gene,&num_load);
/* 输出原始信息至校核文件 */
fprintf(fpc," * * * 系 统 原 始 信 息 及 数 据 * * *:\n\n");
fprintf(fpc," 1. 节点数: %d\n",num_node);
fprintf(fpc," 2. 线路数: %d\n",num_line);
fprintf(fpc," 3. 变压器数: %d\n",num_tran);
fprintf(fpc," 4. 发电机数: %d\n",num_gene);
fprintf(fpc," 5. 负荷数: %d\n",num_load);
/* 判断节点数是否超过宏定义数 */
if(num_node>node_max)
{
printf(" 注意, 不正常停止 !!!\n\n 需要改变宏定义'node_max'为： %d\n",num_node);
fprintf(fpo," 注意, 不正常停止 !!!\n\n 需要改变宏定义'node_max'为： %d\n",num_node);
exit(0);
}
/* 读入线路数据 */
if(num_line>0)
{
line=(struct data *)calloc(num_line+1,sizeof(struct data)); // 用于存放线路数据
if(line==NULL)
{
printf("No memory in line!\n");
fprintf(fpo,"No memory in line!\n");
exit(0);
}
for(p=line+1;p<=line+num_line;p++)
fscanf(fpi,"%d,%d,%lf,%lf,%lf;",&p->i,&p->j,&p->a,&p->b,&p->c);
// 依次读入线路两端的节点i号,节点j号,电阻、电抗和电纳
fprintf(fpc," 线路数据:\n");
fprintf(fpc," 序号 节点i 节点j 电阻 电抗 电纳\n");
for(i=1;i<=num_line;i++)
fprintf(fpc,"%8d%8d%8d%14.5f%14.5f%14.5f\n",i,line[i].i,line[i].j,line[i].a,line[i].b,line[i].c);
fprintf(fpc,"\n");
}
else
fprintf(fpc," 系统中没有线路!\n\n");
/* 读入变压器数据 */
if(num_tran>0)
{
tran=(struct data *)calloc(num_tran+1,sizeof(struct data)); // 用于放变压器数据
if(tran==NULL)
{
printf("No memory in tran!\n");
fprintf(fpo,"No memory in tran!\n");
exit(0);
}
for(p=tran+1;p<=tran+num_tran;p++)
fscanf(fpi,"%d,%d,%lf,%lf,%lf;",&p->i,&p->j,&p->a,&p->b,&p->c);
// 依次读入变压器支路两端的节点i号,节点j号,电阻、电抗和变比
fprintf(fpc," 变压器数据:\n");
fprintf(fpc," 序号 节点i 节点j 电阻 电抗 变比\n");
for(i=1;i<=num_tran;i++)
fprintf(fpc,"%8d%8d%8d%14.5f%14.5f%14.5f\n",i,tran[i].i,tran[i].j,tran[i].a,tran[i].b,tran[i].c);
fprintf(fpc,"\n");
}
else
fprintf(fpc," 系统中没有变压器!\n\n");
/* 读入发电机数据 */
if(num_gene>0)
{
gene=(struct data *)calloc(num_gene+1,sizeof(struct data)); // 用于放发电机数据
if(gene==NULL)
{
printf("No memory in gene!\n");
fprintf(fpo,"No memory in gene!\n");
exit(0);
}
for(p=gene+1;p<=gene+num_gene;p++)
fscanf(fpi,"%d,%d,%lf,%lf,%lf;",&p->i,&p->j,&p->a,&p->b,&p->c);
// 依次读入发电机端的节点号,类型,有功,无功和电压
fprintf(fpc," 发电机数据:\n");
fprintf(fpc," 序号 节点号 类型 有功 无功 电压\n");
for(i=1;i<=num_gene;i++)
fprintf(fpc,"%8d%8d%8d%14.5f%14.5f%14.5f\n",i,gene[i].i,gene[i].j,gene[i].a,gene[i].b,gene[i].c);
fprintf(fpc,"\n");
}
else
{
fprintf(fpo," 系统中没有发电机 ?!\n\n");
exit(0);
}
/* 读入负荷数据 */
if(num_load>0)
{
load=(struct data1 *)calloc(num_load+1,sizeof(struct data1));// 用于放负荷数据
if(load==NULL)
{
printf("No memory in load!\n");
fprintf(fpo,"No memory in load!\n");
exit(0);
}
for(p1=load+1;p1<=load+num_load;p1++)
fscanf(fpi,"%d,%lf,%lf;",&p1->i,&p1->a,&p1->b);
// 依次读入负荷所在节点号、有功和无功
fprintf(fpc," 负荷数据:\n");
fprintf(fpc," 序号 节点号 有功 无功\n");
for(i=1;i<=num_load;i++)
fprintf(fpc,"%8d%8d%14.5f%14.5f\n",i,load[i].i,load[i].a,load[i].b);
fprintf(fpc,"\n");
}
else
{
fprintf(fpo," 系统中没有负荷 ?!\n\n");
exit(0);
}
/* 关闭输入文件 */
if(fclose(fpi)!=0)
{
printf("Input file cannot be closed!\n");
fprintf(fpo,"Input file cannot be closed!\n");
exit(1);
}
/* 申请数组变量空间 */
yr=TwoArrayAlloc_D(num_node+1,num_node+1); // 用于存放导纳矩阵实部
yi=TwoArrayAlloc_D(num_node+1,num_node+1); // 用于存放导纳矩阵虚部
vm=(double *)calloc(num_node+1,sizeof(double)); // 用于存放节点电压幅值
if(vm==NULL)
{
printf("No memory in vm!\n");
fprintf(fpo,"No memory in vm!\n");
exit(0);
}
va=(double *)calloc(num_node+1,sizeof(double)); // 用于存放节点电压相位
if(va==NULL)
{
printf("No memory in va!\n");
fprintf(fpo,"No memory in va!\n");
exit(0);
}
ppq_node=(double *)calloc((num_node+1)*2,sizeof(double)); // 用于存放节点功率
if(ppq_node==NULL)
{
printf("No memory in ppq_node!\n");
fprintf(fpo,"No memory in ppq_node!\n");
exit(0);
}
for(i=0;i<=num_node;i++)
{
pq_node[i]=&ppq_node[i*2];
}
ppq_load=(double *)calloc((num_node+1)*2,sizeof(double)); // 用于存放负荷功率
if(ppq_load==NULL)
{
printf("No memory in ppq_load!\n");
fprintf(fpo,"No memory in ppq_load!\n");
exit(0);
}
for(i=0;i<=num_node;i++)
{
pq_load[i]=&ppq_load[i*2];
}
/* 形成导纳矩阵 */
matrix_adm();
/* 输出导纳矩阵 */
fprintf(fpc," 导纳矩阵对角线元素:\n");
fprintf(fpc," 节点号 实部 虚部\n");
for(i=1;i<=num_node;i++)
fprintf(fpc,"%8d%14.5f%14.5f\n",i,yr[i][i],yi[i][i]);
fprintf(fpc,"\n");
fprintf(fpc," 导纳矩阵非零非对角线元素:\n");
fprintf(fpc," 节点号 节点号 实部 虚部\n");
for(i=1;i<=num_node-1;i++)
for(j=i+1;j<=num_node;j++)
if((yr[i][j]!=0) || (yi[i][j]!=0))
fprintf(fpc,"%8d%8d%14.5f%14.5f\n",i,j,yr[i][j],yi[i][j]);
fprintf(fpc,"\n");
fprintf(fpc," 总非零元素个数%8d%\n",num_bran);
/* 用牛顿法计算潮流 */
ntlf();
/* 不收敛情况 */
if(conv==1)
{
fprintf(fpo," * * * 注 意 ! 潮 流 不 收 敛 * * *\n");
printf(" * * * 注 意 ! 潮 流 不 收 敛 * * *\n");
exit(0);
}
printf("\n");
fprintf(fpo,"\n");
/* 潮流输出 */
fprintf(fpo," * * * * * * 潮 流 计 算 结 果 * * * * * *\n\n");
fprintf(fpo,"1. 节 点 电 压 及 功 率 :\n\n");
fprintf(fpo," 节点号 电压幅值 电压相位 发电机有功 发电机无功 负荷有功 负荷无功\n");
for(i=1;i<=num_node;i++)
fprintf(fpo,"%8d%13.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",
i,vm[i],va[i]*deg,pq_node[i][0]+pq_load[i][0],pq_node[i][1]+pq_load[i][1],
pq_load[i][0],pq_load[i][1]);
fprintf(fpo,"\n");
fprintf(fpo,"2. 支 路 功 率 分 布 :\n\n");
fprintf(fpo," 节点号 节点号 Pij Qij Pji Qji\n");
for(i=1;i<=num_line;i++) // 线路功率输出
{
ii=line[i].i; jj=line[i].j;
r=line[i].a; x=line[i].b;
b=r*r+x*x; r=r/b; x=-x/b;
b=line[i].c;
v1=vm[ii]; v2=vm[jj];
a1=va[ii]-va[jj];
pij=v1*(v1*r+v2*(-r*cos(a1)-x*sin(a1)));
qij=v1*(-v1*(b+x)+v2*(-r*sin(a1)+x*cos(a1)));
pji=v2*(v2*r+v1*(-r*cos(-a1)-x*sin(-a1)));
qji=v2*(-v2*(b+x)+v1*(-r*sin(-a1)+x*cos(-a1)));
fprintf(fpo,"%8d%8d%12.4f%12.4f%12.4f%12.4f\n",ii,jj,pij,qij,pji,qji);
}
for(i=1;i<=num_tran;i++)
{
ii=tran[i].i; jj=tran[i].j;
r=tran[i].a; x=tran[i].b;
b=r*r+x*x; r=r/b; x=-x/b;
b=tran[i].c;
v1=vm[ii]; v2=vm[jj];
a1=va[ii]-va[jj];
pij=v1*(v1*r+v2*(-r*cos(a1)-x*sin(a1))/b);
qij=v1*(-v1*x+v2*(-r*sin(a1)+x*cos(a1))/b);
pji=v2*(v2*r/b+v1*(-r*cos(-a1)-x*sin(-a1)))/b;
qji=v2*(-v2*x/b+v1*(-r*sin(-a1)+x*cos(-a1)))/b;
fprintf(fpo,"%8d%8d%12.4f%12.4f%12.4f%12.4f\n",ii,jj,pij,qij,pji,qji);
}
/* 关闭输出文件 */
if(fclose(fpo)!=0)
{
printf("Output file cannot be closed!\n");
exit(1);
}
/* 释放申请过的空间 */
free(line);
free(tran);
free(gene);
free(load);
free(ppq_node);
free(ppq_load);
TwoArrayFree_D(yr);
TwoArrayFree_D(yi);
}
/* 形成导纳矩阵forming admittance matrix */
void matrix_adm()
{
int i,j;
double r,x,b;
struct data *p,*end;
end=line+num_line;
for(p=line+1;p<=end;p++) // 形成有关线路部分的导纳矩阵
{
i=p->i;
j=p->j;
r=p->a;
x=p->b;
b=r*r+x*x;
r/=b;
x/=-b;
if(i==j)
{
yr[i][i]+=r;
yi[i][i]+=x;
continue;
}
b=p->c;
yr[i][i]+=r;
yi[i][i]+=x+b;
yr[j][j]+=r;
yi[j][j]+=x+b;
yr[i][j]+=-r;
yi[i][j]+=-x;
yr[j][i]+=-r;
yi[j][i]+=-x;
}
end=tran+num_tran;
for(p=tran+1;p<=end;p++) // 形成有关变压器部分的导纳矩阵
{
i=p->i;
j=p->j;
r=p->a;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -