IEEE30节点系统潮流计算

admin · 发表于 2023-11-25 14:17:29

// 30节点系统潮流计算程序.cpp : Defines the entry point for the console application.
#include "stdafx.h"
//*************************************************************************************************************//
//IEEE-30节点系统潮流计算程序
#include<stdio.h>
#include<math.h>
#define N 30 //节点数
#define n_PQ 24 //PQ节点数
#define n_PV 5 //PV节点数
#define n_br 41 //串联支路数
#define n_C 2 //装有并联电容节点数
//#define n_L 1
#define PI 3.1415926
double K1,K2,K3,K4; //变压器变比
double Yc1,Yc2;
//double Xmcr;
void main()
{
//void disp_matrix(double *disp_p,int disp_m,int disp_n); //矩阵显示函数
//节点电压赋初值，PV节点电压幅值已知，相角置0；平衡节点电压幅值和相角均已知；PQ节点电压幅值设1，相角设0.
// 1 2 3 4 5 6 7 8 9 10
/*double Us[2*N]={1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,
1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,
1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0338,0,1.0058,0,1.0230,0,1.0913,0,1.0883,0,1.0500,0
}; */
/*double Us[2*N]={1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,
1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,
1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0500,0,1.0500,0,1.0500,0,1.0500,0,1.0500,0,1.0500,0
};*/
double Us[2*N]={1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,
1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,
1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0000,0,1.0500,0
};
//PV节点发电机有功赋初值（除平衡节点外）
// 1 2 3 4 5 6 7 8 9 10
double Ps[N]={0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
0.0000,0.0000,0.0000,0.0000,0.5756,0.2456,0.3500,0.1793,0.1691
};
//发电机无功初值置0
double Qs[N]={0};
//各节点有功负荷赋值
// 1 2 3 4 5 6 7 8 9 10
double PL[N]={0.024,0.076,0.000,0.228,0.000,0.058,0.112,0.062,0.082,0.035,
0.090,0.032,0.095,0.022,0.175,0.000,0.032,0.087,0.000,0.035,
0.000,0.000,0.024,0.106,0.217,0.942,0.300,0.000,0.000
};
//各节点无功负荷赋值
// 1 2 3 4 5 6 7 8 9 10
double QL[N]={0.012,0.016,0.000,0.109,0.000,0.020,0.075,0.016,0.025,0.018,
0.058,0.009,0.034,0.007,0.112,0.000,0.016,0.067,0.000,0.023,
0.000,0.000,0.009,0.019,0.127,0.190,0.300,0.000,0.000
};
//double K1=1.0155;K2=0.9629;K3=1.0129;K4=0.9581; //变压器当前变比
//double Yc1=0.19;Yc2=0.04;
double K1=1;K2=1;K3=1;K4=1; //变压器当前变比
double Yc1=0;Yc2=0;
//Xmcr=8.0;
//double Y[N][N],G[N][N],B[N][N],B1[N-1][N-1],B2[n_PQ][n_PQ],invB1[N-1][N-1],invB2[n_PQ][n_PQ]; //阻抗参数
double G[N][N],B[N][N],B1[N-1][N-1],B2[n_PQ][n_PQ],invB1[N-1][N-1],invB2[n_PQ][n_PQ]; //阻抗参数
//结构体定义线路参数
struct
{
int nl; //左节点
int nr; //右节点
double R; //串联电阻值
double X; //串联电抗值
double Bl; //左节点充电电纳
double Br; //右节点充电电纳
}ydata[n_br+n_C]={ // 左右 R X B/2 B/2
/*1*/ {29,24,0.0192,0.0575,0.0264,0.0264}, //一般支路参数
/*2*/ {29,0,0.0452,0.1852,0.0204,0.0204},
/*3*/ {24,1,0.0570,0.1737,0.0184,0.0184},
/*4*/ {0,1,0.0132,0.0379,0.0042,0.0042},
/*5*/ {24,25,0.0472,0.1983,0.0209,0.0209},
/*6*/ {24,2,0.0581,0.1763,0.0187,0.0187},
/*7*/ {1,2,0.0119,0.0414,0.0045,0.0045},
/*8*/ {25,3,0.0460,0.1160,0.0102,0.0102},
/*9*/ {2,3,0.0267,0.0820,0.0085,0.0085},
/*10*/ {2,26,0.0120,0.0420,0.0045,0.0045},
/*11*/ {4,27,0.0000,0.2080,0.0000,0.0000},
/*12*/ {4,5,0.0000,0.1100,0.0000,0.0000},
/*13*/ {6,28,0.0000,0.1400,0.0000,0.0000},
/*14*/ {6,7,0.1231,0.2559,0.0000,0.0000},
/*15*/ {6,8,0.0662,0.1304,0.0000,0.0000},
/*16*/ {6,9,0.0945,0.1987,0.0000,0.0000},
/*17*/ {7,8,0.2210,0.1997,0.0000,0.0000},
/*18*/ {9,10,0.0824,0.1932,0.0000,0.0000},
/*19*/ {8,11,0.1070,0.2185,0.0000,0.0000},
/*20*/ {11,12,0.0639,0.1292,0.0000,0.0000},
/*21*/ {12,13,0.0340,0.0680,0.0000,0.0000},
/*22*/ {5,13,0.0936,0.2090,0.0000,0.0000},
/*23*/ {5,10,0.0324,0.0845,0.0000,0.0000},
/*24*/ {5,14,0.0348,0.0749,0.0000,0.0000},
/*25*/ {5,15,0.0727,0.1499,0.0000,0.0000},
/*26*/ {14,15,0.0116,0.0236,0.0000,0.0000},
/*27*/ {8,16,0.1000,0.2020,0.0000,0.0000},
/*28*/ {15,17,0.1150,0.1790,0.0000,0.0000},
/*29*/ {16,17,0.1320,0.2700,0.0000,0.0000},
/*30*/ {17,18,0.1885,0.3292,0.0000,0.0000},
/*31*/ {18,19,0.2554,0.3800,0.0000,0.0000},
/*32*/ {18,20,0.1093,0.2087,0.0000,0.0000},
/*33*/ {20,22,0.2198,0.4153,0.0000,0.0000},
/*34*/ {20,23,0.3202,0.6027,0.0000,0.0000},
/*35*/ {22,23,0.2399,0.4533,0.0000,0.0000},
/*36*/ {26,21,0.0636,0.2000,0.0214,0.0214},
/*37*/ {2,21,0.0169,0.0599,0.0065,0.0065},
/*38*/ {2,4,0.0000,K1*0.2080,(K1-1)/(K1*0.2080),(1-K1)/(K1*K1*0.2080)}, //变压器支路参数
/*39*/ {2,5,0.0000,K2*0.5560,(K2-1)/(K2*0.5560),(1-K2)/(K2*K2*0.5560)},
/*40*/ {1,6,0.0000,K3*0.2560,(K3-1)/(K3*0.2560),(1-K3)/(K3*K3*0.2560)},
/*41*/ {21,20,0.0000,K4*0.3960,(K4-1)/(K4*0.3960),(1-K4)/(K4*K4*0.3960)},
{5,5,0.0000,Yc1,0.0000,0.0000}, //5号节点有并联电容，电纳标幺值为0.19
{17,17,0.0000,Yc2,0.0000,0.0000}, //17号节点有并联电容，电纳标幺值为0.02
};
double Z2; //Z^2=R^2+X^2 各串联阻抗值的平方
double U_amp[N],U_ang[N]; //电压幅值和相角
double dU_amp[N],dU_ang[N]; //电压幅值和相角修正值
double mid1[N],mid2[N],dP[N-1],dQ[n_PQ]; //mid1、mid2存储计算雅克比行列式对角线元素的中间值,dS存储PQU的不平衡量
double MIN=0.00001; //PQU不平衡量最大值
double dPQU; //PQU不平衡量
double Ploss;
//double S[2*n_br]; //线路功率
double P[n_br][n_br];
double Q[n_br][n_br];
double Ps_total,Qs_total; //发电机总有功、总无功
double PL_total,QL_total; //负荷总有功、总无功
//double Pc; //电容器的功率
int kk=0; //迭代次数
int kkmax=60; //最大迭代次数
int i,j,k;
double t;
//形成导纳矩阵--------------------------------------------------------------------------------------------------
for(i=0;i<N;i++)
for(j=0;j<N;j++) //各节点导纳赋初值为
{
G[i][j]=0;
B[i][j]=0;
}
for(i=0;i<n_br+n_C;i++)
{
if(ydata[i].nl!=ydata[i].nr) //左节点号与右节点号不同
{
Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X); //阻抗的平方
//串联阻抗等效导纳值
//非对角元素
G[ydata[i].nl][ydata[i].nr]=(-ydata[i].R)/Z2;
B[ydata[i].nl][ydata[i].nr]=ydata[i].X/Z2;
G[ydata[i].nr][ydata[i].nl]=(-ydata[i].R)/Z2;
B[ydata[i].nr][ydata[i].nl]=ydata[i].X/Z2;
//对角元素
G[ydata[i].nl][ydata[i].nl]+=ydata[i].R/Z2;
G[ydata[i].nr][ydata[i].nr]+=ydata[i].R/Z2;
B[ydata[i].nl][ydata[i].nl]+=(-ydata[i].X/Z2);
B[ydata[i].nr][ydata[i].nr]+=(-ydata[i].X/Z2);
//节点自导纳需加上充电导纳值
B[ydata[i].nl][ydata[i].nl]+=ydata[i].Bl;
B[ydata[i].nr][ydata[i].nr]+=ydata[i].Br;
}
else //左节点号=右节点号，即节点有并联阻抗的情况
{
B[ydata[i].nl][ydata[i].nl]+=ydata[i].X;
}
}
//输出节点导纳矩阵
//printf("Y=\n");
//for(i=0;i<N;i++)
//{
// for(j=0;j<N;j++)
// {
// printf("%.6f+j%.6f ",G[i][j],B[i][j]);
// }
// printf("\n");
//}
////printf("\n");
//printf("G=\n"); //矩阵G
//disp_matrix(*G,N,N);
//printf("B=\n"); //矩阵B
//disp_matrix(*B,N,N);
//形成B',B''矩阵
for(i=0;i<N-1;i++)
for(j=0;j<N-1;j++)
{
B1[i][j]=B[i][j];
if(i<n_PQ&&j<n_PQ)
B2[i][j]=B[i][j];
}
/*printf("B1=\n"); //输出B1矩阵
disp_matrix(*B1,N-1,N-1);
printf("B2=\n"); //输出B2矩阵
disp_matrix(*B2,n_PQ,n_PQ);*/
//B'求逆-----------------------------------------------------------------------------------
for(i=0;i<N-1;i++)
for(j=0;j<N-1;j++)
{
if(i!=j)
invB1[i][j]=0; //B'的逆矩阵，非对角线元素初值置0
else
invB1[i][j]=1; //B'的逆矩阵，对角线元素初值置1
}
for(i=0;i<N-1;i++)
{
for(j=0;j<N-1;j++)
{
if(i!=j)
{
t=B1[j][i]/B1[i][i];
for(k=0;k<N-1;k++)
{
B1[j][k]-=B1[i][k]*t;
invB1[j][k]-=invB1[i][k]*t;
}
}
}
}
for(i=0;i<N-1;i++)
if(B1[i][i]!=1)
{
t=B1[i][i];
for(j=0;j<N-1;j++)
invB1[i][j]=invB1[i][j]/t;
}
//printf("invB1=\n");
//disp_matrix(*invB1,N-1,N-1);
//B''求逆
for(i=0;i<n_PQ;i++)
for(j=0;j<n_PQ;j++)
{
if(i!=j)
invB2[i][j]=0;
else
invB2[i][j]=1;
}
for(i=0;i<n_PQ;i++)
{
for(j=0;j<n_PQ;j++)
{
if(i!=j)
{
t=B2[j][i]/B2[i][i];
for(k=0;k<n_PQ;k++)
{
B2[j][k]-=B2[i][k]*t;
invB2[j][k]-=invB2[i][k]*t;
}
}
}
}
for(i=0;i<n_PQ;i++)
if(B2[i][i]!=1)
{
t=B2[i][i];
for(j=0;j<n_PQ;j++)
invB2[i][j]=invB2[i][j]/t;
}
/*printf("invB2=\n");
disp_matrix(*invB2,n_PQ,n_PQ);*/
//分离并输出U,δ
for(i=0;i<N;i++)
{
U_amp[i]=Us[2*i]; //电压幅值
U_ang[i]=Us[2*i+1]; //电压相角
}
/* //printf("U_amp=");
for(i=0;i<N;i++)
printf("U[%d]=%.4f ",i,U_amp[i]);
printf("\n");
//printf("\nU_ang=");
for(i=0;i<N;i++)
printf("δ[%d]=%.2f ",i,U_ang[i]);
printf("\n\n");
*/
//主程序=====================================================================================================
for(kk=0;kk<kkmax;kk++)
{
for(i=0;i<N-1;i++)
{
mid1[i]=0; //中间变量,初值置
mid2[i]=0;
for(j=0;j<N;j++) //
{
mid1[i]+=U_amp[j]*(G[i][j]*cos(U_ang[i]-U_ang[j])+B[i][j]*sin(U_ang[i]-U_ang[j])); //计算Pi
mid2[i]+=U_amp[j]*(G[i][j]*sin(U_ang[i]-U_ang[j])-B[i][j]*cos(U_ang[i]-U_ang[j])); //计算Qi
}
dP[i]=Ps[i]-PL[i]-U_amp[i]*mid1[i]; //计算△Pi
if(i<n_PQ)
dQ[i]=Qs[i]-QL[i]-U_amp[i]*mid2[i];
}
//printf("kk=%d\n",kk);
//printf("dP="); //第k次迭代dP
//for(i=0;i<N-1;i++)
// printf("%9.6f ",dP[i]);
//printf("\ndQ=");
//for(i=0;i<n_PQ;i++)
// printf("%9.6f ",dQ[i]);
//求不平衡量dP、dQ的最大值
dPQU=0;
for(i=0;i<N-1;i++)
{
if(fabs(dP[i])>dPQU)
dPQU=fabs(dP[i]);
}
for(i=0;i<n_PQ;i++)
{
if(fabs(dQ[i])>dPQU)
dPQU=fabs(dQ[i]);
}
//printf("\ndPQU=%9.6f",dPQU);
//printf("\n");
if(dPQU<MIN)
break;
else
{
//求P,δ修正值
for(i=0;i<N-1;i++)
{
dP[i]=dP[i]/U_amp[i]; //△Pi/Vi
}
for(i=0;i<N-1;i++)
dU_ang[i]=0; //△δ初值置0
for(i=0;i<N-1;i++)
{
for(j=0;j<N-1;j++)
dU_ang[i]+=-invB1[i][j]*dP[j]; //求Vi*△δi
dU_ang[i]=dU_ang[i]/U_amp[i]; //求△δ
}
for(i=0;i<N-1;i++)
U_ang[i]+=dU_ang[i]; //修正δi
//Q,U修正值
for(i=0;i<n_PQ;i++)
dQ[i]=dQ[i]/U_amp[i]; //△Qi/Vi
for(i=0;i<n_PQ;i++)
dU_amp[i]=0; //dU初值置0
for(i=0;i<n_PQ;i++)
for(j=0;j<n_PQ;j++)
dU_amp[i]+=-invB2[i][j]*dQ[j];
for(i=0;i<n_PQ;i++)
U_amp[i]+=dU_amp[i]; //修正Ui
}
}
//循环结束---------------------------------------------------------------------------------------------------------------
//求平衡节点功率---------------------------------------------------------------------------------------------------
mid1[N-1]=0; //平衡节点编号为第N-1号
mid2[N-1]=0;
for(j=0;j<N;j++)
{
mid1[N-1]+=U_amp[j]*(G[N-1][j]*cos(U_ang[N-1]-U_ang[j])+B[N-1][j]*sin(U_ang[N-1]-U_ang[j]));
mid2[N-1]+=U_amp[j]*(G[N-1][j]*sin(U_ang[N-1]-U_ang[j])-B[N-1][j]*cos(U_ang[N-1]-U_ang[j]));
}
Ps[N-1]=U_amp[N-1]*mid1[N-1];
Qs[N-1]=U_amp[N-1]*mid2[N-1];
//PV节点的无功出力-------------------------------------------------------------------------------------------------
for(i=n_PQ;i<N-1;i++)
Qs[i]=U_amp[i]*mid2[i]+QL[i]; //发电机的等效注入无功
//总发电机有功、无功计算和总负荷有功、无功计算-------------------------------------------------------------------------------------------------
Ps_total=0;Qs_total=0; //发电机总有功、总无功初值置0
PL_total=0;QL_total=0; //负荷总有功、总无功初值置0
for(i=0;i<N;i++)
{
Ps_total+=Ps[i];
Qs_total+=Qs[i];
PL_total+=PL[i];
QL_total+=QL[i];
}
//全部线路功率-------------------------------------------------------------------------------------------------
for(i=0;i<n_br;i++)
{
int a,b; //定义两个中间变量，表示首末节点
a=ydata[i].nl;
b=ydata[i].nr;
double angle=PI/2-atan(ydata[i].X/ydata[i].R);
double Z=sqrt(ydata[i].R*ydata[i].R+ydata[i].X*ydata[i].X);
//printf("%d,%d,%10.6f,%10.6f,%10.6f,%10.6f\n",a,b,U_amp[a],U_amp[b],Z,angle/PI*180);
P[a][b]=U_amp[a]*U_amp[a]*sin(angle)/Z+U_amp[a]*U_amp[b]*sin(U_ang[a]-U_ang[b]-angle)/Z;
P[b][a]=U_amp[b]*U_amp[b]*sin(angle)/Z+U_amp[a]*U_amp[b]*sin(U_ang[b]-U_ang[a]-angle)/Z;
Q[a][b]=U_amp[a]*U_amp[a]*cos(angle)/Z-U_amp[a]*U_amp[b]*cos(U_ang[a]-U_ang[b]-angle)/Z-U_amp[a]*U_amp[a]*ydata[i].Bl;
Q[b][a]=U_amp[b]*U_amp[b]*cos(angle)/Z-U_amp[a]*U_amp[b]*cos(U_ang[b]-U_ang[a]-angle)/Z-U_amp[b]*U_amp[b]*ydata[i].Br;
/*printf("S[%d][%d]=%10.6f+j%.6f\n",a,b,P[a][b],Q[a][b]);
printf("S[%d][%d]=%10.6f+j%.6f\n",b,a,P[b][a],Q[b][a]);*/
}
//并联电容器输出功率-------------------------------------------------------------------------------------------------
//Pc=-U_amp[6]*U_amp[6]/ydata[n_br+n_C-1].X; //输出功率用负号表示
//printf("Pc=%.4f\n",Pc);
//有功损耗--------------------------------------------------------------------------------------------------------------
Ploss=0;
for(i=0;i<N;i++)
{
Ploss+=Ps[i]-PL[i];
}
//电压平均值--------------------------------------------------------------------------------------------------------------
double U_average=0; //所有节点的平均电压
double U_average2=0; //PQ节点的平均电压
for(i=0;i<N;i++)
{
U_average+=U_amp[i];
}
U_average=U_average/N;
for(i=0;i<n_PQ;i++)
{
U_average2+=U_amp[i];
}
U_average2=U_average2/n_PQ;
//---------------------------------------------------------------------------------------------------------------------------
// 输出结果
printf("K1=%.4f, K2=%.4f, K3=%.4f, K4=%.4f, Yc1=%.2f, Yc2=%.2f",K1,K2,K3,K4,Yc1,Yc2); //输出各变压器变比、电容器电纳值
printf("\n\nkk=%d",kk); //最大迭代次数
printf("\n电压幅值="); //电压幅值
for(i=0;i<N;i++)
printf("%10.4f",U_amp[i]);
printf("\n电压相角="); //电压相角
for(i=0;i<N;i++)
printf("%10.4f",U_ang[i]/PI*180); //弧度转换成°
//输出发电机总有功、无功和负荷的总有功、无功
printf("\n\n发电机总有功=%.4f, 发电机总无功=%.4f, 负荷总有功=%.4f，负荷总无功=%.4f",Ps_total,Qs_total,PL_total,QL_total);
printf("\n\n各发电机的注入功率");
printf("\nPs="); //各发电机注入有功功率
for(i=n_PQ;i<N;i++)
printf("%10.4f",Ps[i]);
printf("\nQs="); //各发电机注入无功功率
for(i=n_PQ;i<N;i++)
printf("%10.4f",Qs[i]);
printf("\n\n网损=%.4f, 所有节点平均电压=%.4f, PQ节点平均电压=%.4f",Ploss,U_average,U_average2);
printf("\n\n");
//getchar();
}
//========================================================================================================
/*void disp_matrix(double *A,int m,int n) //矩阵显示函数
{
int i,j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
printf("%10.6f ",*(A+i*m+j));
printf("\n");
}
printf("\n");
}*/

复制代码

		自动登录	找回密码
密码			立即注册

IEEE30节点系统潮流计算

相关帖子