IEEE30节点系统潮流计算
// 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={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={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={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={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={0};
//各节点有功负荷赋值
// 1 2 3 4 5 6 7 8 9 10
double PL={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={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,G,B,B1,B2,invB1,invB2; //阻抗参数
double G,B,B1,B2,invB1,invB2; //阻抗参数
//结构体定义线路参数
struct
{
int nl; //左节点
int nr; //右节点
double R; //串联电阻值
double X; //串联电抗值
double Bl; //左节点充电电纳
double Br; //右节点充电电纳
}ydata={ // 左 右 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,U_ang; //电压幅值和相角
double dU_amp,dU_ang; //电压幅值和相角修正值
double mid1,mid2,dP,dQ; //mid1、mid2存储计算雅克比行列式对角线元素的中间值,dS存储PQU的不平衡量
double MIN=0.00001; //PQU不平衡量最大值
double dPQU; //PQU不平衡量
double Ploss;
//double S; //线路功率
double P;
double Q;
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=0;
B=0;
}
for(i=0;i<n_br+n_C;i++)
{
if(ydata.nl!=ydata.nr) //左节点号与右节点号不同
{
Z2=(ydata.R)*(ydata.R)+(ydata.X)*(ydata.X); //阻抗的平方
//串联阻抗等效导纳值
//非对角元素
G.nl].nr]=(-ydata.R)/Z2;
B.nl].nr]=ydata.X/Z2;
G.nr].nl]=(-ydata.R)/Z2;
B.nr].nl]=ydata.X/Z2;
//对角元素
G.nl].nl]+=ydata.R/Z2;
G.nr].nr]+=ydata.R/Z2;
B.nl].nl]+=(-ydata.X/Z2);
B.nr].nr]+=(-ydata.X/Z2);
//节点自导纳需加上充电导纳值
B.nl].nl]+=ydata.Bl;
B.nr].nr]+=ydata.Br;
}
else //左节点号=右节点号，即节点有并联阻抗的情况
{
B.nl].nl]+=ydata.X;
}
}
//输出节点导纳矩阵
//printf("Y=\n");
//for(i=0;i<N;i++)
//{
// for(j=0;j<N;j++)
// {
// printf("%.6f+j%.6f ",G,B);
// }
// 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=B;
if(i<n_PQ&&j<n_PQ)
B2=B;
}
/*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=0; //B'的逆矩阵，非对角线元素初值置0
else
invB1=1; //B'的逆矩阵，对角线元素初值置1
}
for(i=0;i<N-1;i++)
{
for(j=0;j<N-1;j++)
{
if(i!=j)
{
t=B1/B1;
for(k=0;k<N-1;k++)
{
B1-=B1*t;
invB1-=invB1*t;
}
}
}
}
for(i=0;i<N-1;i++)
if(B1!=1)
{
t=B1;
for(j=0;j<N-1;j++)
invB1=invB1/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=0;
else
invB2=1;
}
for(i=0;i<n_PQ;i++)
{
for(j=0;j<n_PQ;j++)
{
if(i!=j)
{
t=B2/B2;
for(k=0;k<n_PQ;k++)
{
B2-=B2*t;
invB2-=invB2*t;
}
}
}
}
for(i=0;i<n_PQ;i++)
if(B2!=1)
{
t=B2;
for(j=0;j<n_PQ;j++)
invB2=invB2/t;
}
/*printf("invB2=\n");
disp_matrix(*invB2,n_PQ,n_PQ);*/
//分离并输出U,δ
for(i=0;i<N;i++)
{
U_amp=Us; //电压幅值
U_ang=Us; //电压相角
}
/* //printf("U_amp=");
for(i=0;i<N;i++)
printf("U[%d]=%.4f",i,U_amp);
printf("\n");
//printf("\nU_ang=");
for(i=0;i<N;i++)
printf("δ[%d]=%.2f",i,U_ang);
printf("\n\n");
*/
//主程序=====================================================================================================
for(kk=0;kk<kkmax;kk++)
{
for(i=0;i<N-1;i++)
{
mid1=0; //中间变量,初值置
mid2=0;
for(j=0;j<N;j++) //
{
mid1+=U_amp*(G*cos(U_ang-U_ang)+B*sin(U_ang-U_ang)); //计算Pi
mid2+=U_amp*(G*sin(U_ang-U_ang)-B*cos(U_ang-U_ang)); //计算Qi
}
dP=Ps-PL-U_amp*mid1; //计算△Pi
if(i<n_PQ)
dQ=Qs-QL-U_amp*mid2;
}
//printf("kk=%d\n",kk);
//printf("dP="); //第k次迭代dP
//for(i=0;i<N-1;i++)
// printf("%9.6f",dP);
//printf("\ndQ=");
//for(i=0;i<n_PQ;i++)
// printf("%9.6f",dQ);
//求不平衡量dP、dQ的最大值
dPQU=0;
for(i=0;i<N-1;i++)
{
if(fabs(dP)>dPQU)
dPQU=fabs(dP);
}
for(i=0;i<n_PQ;i++)
{
if(fabs(dQ)>dPQU)
dPQU=fabs(dQ);
}
//printf("\ndPQU=%9.6f",dPQU);
//printf("\n");
if(dPQU<MIN)
break;
else
{
//求P,δ修正值
for(i=0;i<N-1;i++)
{
dP=dP/U_amp; //△Pi/Vi
}
for(i=0;i<N-1;i++)
dU_ang=0; //△δ初值置0
for(i=0;i<N-1;i++)
{
for(j=0;j<N-1;j++)
dU_ang+=-invB1*dP; //求Vi*△δi
dU_ang=dU_ang/U_amp; //求△δ
}
for(i=0;i<N-1;i++)
U_ang+=dU_ang; //修正δi
//Q,U修正值
for(i=0;i<n_PQ;i++)
dQ=dQ/U_amp; //△Qi/Vi
for(i=0;i<n_PQ;i++)
dU_amp=0; //dU初值置0
for(i=0;i<n_PQ;i++)
for(j=0;j<n_PQ;j++)
dU_amp+=-invB2*dQ;
for(i=0;i<n_PQ;i++)
U_amp+=dU_amp; //修正Ui
}
}
//循环结束---------------------------------------------------------------------------------------------------------------
//求平衡节点功率---------------------------------------------------------------------------------------------------
mid1=0; //平衡节点编号为第N-1号
mid2=0;
for(j=0;j<N;j++)
{
mid1+=U_amp*(G*cos(U_ang-U_ang)+B*sin(U_ang-U_ang));
mid2+=U_amp*(G*sin(U_ang-U_ang)-B*cos(U_ang-U_ang));
}
Ps=U_amp*mid1;
Qs=U_amp*mid2;
//PV节点的无功出力-------------------------------------------------------------------------------------------------
for(i=n_PQ;i<N-1;i++)
Qs=U_amp*mid2+QL; //发电机的等效注入无功
//总发电机有功、无功计算和总负荷有功、无功计算-------------------------------------------------------------------------------------------------
Ps_total=0;Qs_total=0; //发电机总有功、总无功初值置0
PL_total=0;QL_total=0; //负荷总有功、总无功初值置0
for(i=0;i<N;i++)
{
Ps_total+=Ps;
Qs_total+=Qs;
PL_total+=PL;
QL_total+=QL;
}
//全部线路功率-------------------------------------------------------------------------------------------------
for(i=0;i<n_br;i++)
{
int a,b; //定义两个中间变量，表示首末节点
a=ydata.nl;
b=ydata.nr;
double angle=PI/2-atan(ydata.X/ydata.R);
double Z=sqrt(ydata.R*ydata.R+ydata.X*ydata.X);
//printf("%d,%d,%10.6f,%10.6f,%10.6f,%10.6f\n",a,b,U_amp,U_amp,Z,angle/PI*180);
P=U_amp*U_amp*sin(angle)/Z+U_amp*U_amp*sin(U_ang-U_ang-angle)/Z;
P=U_amp*U_amp*sin(angle)/Z+U_amp*U_amp*sin(U_ang-U_ang-angle)/Z;
Q=U_amp*U_amp*cos(angle)/Z-U_amp*U_amp*cos(U_ang-U_ang-angle)/Z-U_amp*U_amp*ydata.Bl;
Q=U_amp*U_amp*cos(angle)/Z-U_amp*U_amp*cos(U_ang-U_ang-angle)/Z-U_amp*U_amp*ydata.Br;
/*printf("S[%d][%d]=%10.6f+j%.6f\n",a,b,P,Q);
printf("S[%d][%d]=%10.6f+j%.6f\n",b,a,P,Q);*/
}
//并联电容器输出功率-------------------------------------------------------------------------------------------------
//Pc=-U_amp*U_amp/ydata.X; //输出功率用负号表示
//printf("Pc=%.4f\n",Pc);
//有功损耗--------------------------------------------------------------------------------------------------------------
Ploss=0;
for(i=0;i<N;i++)
{
Ploss+=Ps-PL;
}
//电压平均值--------------------------------------------------------------------------------------------------------------
double U_average=0; //所有节点的平均电压
double U_average2=0; //PQ节点的平均电压
for(i=0;i<N;i++)
{
U_average+=U_amp;
}
U_average=U_average/N;
for(i=0;i<n_PQ;i++)
{
U_average2+=U_amp;
}
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);
printf("\n电压相角="); //电压相角
for(i=0;i<N;i++)
printf("%10.4f",U_ang/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);
printf("\nQs="); //各发电机注入无功功率
for(i=n_PQ;i<N;i++)
printf("%10.4f",Qs);
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");
}*/
页:
[1]