📄 p-q潮流计算.cpp
字号:
#include <cmath>
#include <iostream>
#define N 10
//#define n 5 /*备用*/
using namespace std;
int main()
{
int n; //统计节点个数
int i,j,k; //循环变量
int r; //PV节点个数
//int r=0;
int m; //PQ节点个数
double R[N][N];/*R[5][5]={{0,0.02,0.08,0,0},{0.02,0,0.06,0.06,0.04},{0.08,0.06,0,0.01,0},
{0,0.06,0.01,0,0.08},{0,0.04,0,0.08,0}};*/ //节点间的电阻
double X[N][N];/*X[5][5]={{0,0.06,0.24,0,0},{0.06,0,0.18,0.18,0.12},{0.24,0.18,0,0.03,0},
{0,0.18,0.03,0,0.24},{0,0.12,0,0.24,0}};*/ //节点间的电抗
double g1[N][N]={0}; //节点间的电导
double b1[N][N]={0}; //节点间的电纳
double G[N][N]={0};/*={{6.25,-5.0,-1.25,0,0},{-5.0,10.834,-1.667,-1.667,-2.5},
{-1.25,-1.667,12.917,-10.0,0},{0,-1.667,-10.0,12.917,-1.25},
{0,-2.5,0,-1.25,3.75}};
导纳矩阵电导*/
double B[N][N]={0};/*={{-18.75,15.0,3.75,0,0,},{15.0,-32.5,5.0,5.0,7.5},
{3.75,5.0,-38.75,30.0,0},{0,5.0,30.0,-38.75,3.75},
{0,7.50,0,3.75,-11.25}};
导纳矩阵电纳*/
double P[N]={0};//={0.2,-0.45,-0.4,-0.6}; //节点有功初值，平衡节点除外
double Q[N]={0};//={0.2,-0.15,-0.05,-0.1}; //节点无功初值，平衡节点除外
//double P[5]={0,0.2,-0.45,-0.4,-0.6}; //节点有功初值，平衡节点除外
//double Q[5]={0,0.2,-0.15,-0.05,-0.1}; //节点无功初值，平衡节点除外
double B1[N-1][N-1]; //系数矩阵1
double B2[N-1][N-1]; //系数矩阵2
double A1[N-1][N-1]; //系数矩阵1的逆矩阵
double A2[N-1][N-1]; //系数矩阵2的逆矩阵
double a[N][2*N]; //求逆矩阵时用的矩阵
double w=0,s=0,l,m1; //求逆矩阵时用的变量
double U[N]; //节点电压幅值
double rad[N]; //节点电压幅角
//double U[N]={1.06,1,1,1,1};
//double rad[N]={0,0,0,0,0};
double cp[N]; //有功不平衡量
double cq[N]; //无功不平衡量
cout << "请输入网络节点数：";
cin >> n ;
cout << "请输入PV节点个数(r<n):" << "r=" ;
cin >> r;
m=n-r-1;
cout << "PQ节点个数(n-r-1)为：" << "m=" << m << endl;
cout << "请输入节点间电阻：\n"; /*************************************/
for( i=0 ; i<n ; i++ ) /* 0 0.02 0.08 0 0 */
{ /* 0.02 0 0.06 0.06 0.04 */
for( j=0 ; j<n ; j++ ) /* 0.08 0.06 0 0.01 0 */
{ /* 0 0.06 0.01 0 0.08 */
cin >> R[i][j]; /* 0 0.04 0 0.08 0 */
} /*************************************/
}
cout << "\n请输入节点间电抗：\n"; /*************************************/
for( i=0 ; i<n ; i++ ) /* 0 0.06 0.24 0 0 */
{ /* 0.06 0 0.18 0.18 0.12 */
for( j=0 ; j<n ; j++ ) /* 0.24 0.18 0 0.03 0 */
{ /* 0 0.18 0.03 0 0.24 */
cin >> X[i][j]; /* 0 0.12 0 0.24 0 */
} /*************************************/
}
/*********************************************************************************************
接着计算节点间的导纳值
*********************************************************************************************/
for( i=0 ; i<n ; i++ )
{
for( j=0 ; j<n ; j++ )
{
if( (R[i][j]!=0)||(X[i][j]!=0) )
{
g1[i][j]=R[i][j]/(R[i][j]*R[i][j]+X[i][j]*X[i][j]);
b1[i][j]=-X[i][j]/(R[i][j]*R[i][j]+X[i][j]*X[i][j]);
}
else
{
g1[i][j]=0;
b1[i][j]=0;
}
}
}
/*********************************************************************************************
下面形成导纳矩阵
*********************************************************************************************/
for( i=0 ; i<n ; i++ )
{
for( j=0 ; j<n ; j++ )
{
if(j!=i) G[i][j]=-g1[i][j];
else for( k=0 ; k<n ; k++ )
{
G[i][j]+=g1[i][k];
}
}
}
for( i=0 ; i<n ; i++ )
{
for( j=0 ; j<n ; j++ )
{
if(j!=i) B[i][j]=-b1[i][j];
else for( k=0 ; k<n ; k++ )
{
B[i][j]+=b1[i][k];
}
}
}
/*for( i=0 ; i<n ; i++ )
{
cout << endl;
for( j=0 ; j<n ; j++ )
{
cout << B[i][j] << " ";
}
}*/
/*********************************************************************************************
下面形成系数矩阵B1，B2
*********************************************************************************************/
for( i=0 ; i<n-1 ; i++ )
{
for( j=0 ; j<n-1 ; j++ )
{
B1[i][j]=B[i+1][j+1];
}
}
for( i=0 ; i<n-r-1 ; i++ )
{
for( j=0 ; j<n-r-1 ; j++ )
{
B2[i][j]=B[i+1][j+1];
}
}
/*for( i=0 ; i<n-1 ; i++ ) /*输出系数矩阵B1
{
cout << endl;
for( j=0 ; j<n-1 ; j++ )
{
cout << B1[i][j] << " " ;
}
}
for( i=0 ; i<n-r-1 ; i++ ) /*输出系数矩阵B2
{
cout << endl;
for( j=0 ; j<n-r-1 ; j++ )
{
cout << B2[i][j] << " ";
}
}*/
/*********************************************************************************************
下面形成系数矩阵B1逆矩阵A1
*********************************************************************************************/
int n1=n-1;
for (i=0;i<n1;i++)
{
for (j=0;j<n1;j++)
{
a[i][j]=B1[i][j];
}
}
for (i=0;i<n1;i++)
{
for (j=n1;j<2*n1;j++)
{
if((i+n1)==j)
a[i][j]=1;
else
a[i][j]=0;
}
}
for (i=0;i<n1;i++)
{
w=a[i][i];
for(j=i;j<2*n1;j++)
{
a[i][j]=a[i][j]/w;
}
for(k=i+1;k<n1;k++)
{
s=a[k][i];
for(j=i;j<2*n1;j++)
{
a[k][j]=a[k][j]-s*a[i][j];
}
}
}
for(i=n1-1;i>=0;i--)
{
l=a[i][i];
for(j=i;j<2*n1;j++)
{
a[i][j]=a[i][j]/l;
}
for(k=i-1;k>=0;k--)
{
m1=a[k][i];
for(j=i;j<2*n1;j++)
{
a[k][j]=a[k][j]-m1*a[i][j];
}
}
}
for(i=0;i<n1;i++)
{
for(j=n1;j<2*n1;j++)
{
A1[i][j-n1]=a[i][j] ;
}
}
/*for(i=0;i<n1;i++)
{
cout << endl;
for(j=0; j<n1 ; j++ )
{
cout << A1[i][j] << " ";
}
}*/
/*********************************************************************************************
下面形成系数矩阵B2的逆矩阵A2
*********************************************************************************************/
int n2=n-r-1;
for (i=0;i<n2;i++)
{
for (j=0;j<n2;j++)
{
a[i][j]=B1[i][j];
}
}
for (i=0;i<n2;i++)
{
for (j=n1;j<2*n2;j++)
{
if((i+n2)==j)
a[i][j]=1;
else
a[i][j]=0;
}
}
for (i=0;i<n2;i++)
{
w=a[i][i];
for(j=i;j<2*n2;j++)
{
a[i][j]=a[i][j]/w;
}
for(k=i+1;k<n2;k++)
{
s=a[k][i];
for(j=i;j<2*n2;j++)
{
a[k][j]=a[k][j]-s*a[i][j];
}
}
}
for(i=n2-1;i>=0;i--)
{
l=a[i][i];
for(j=i;j<2*n2;j++)
{
a[i][j]=a[i][j]/l;
}
for(k=i-1;k>=0;k--)
{
m1=a[k][i];
for(j=i;j<2*n2;j++)
{
a[k][j]=a[k][j]-m1*a[i][j];
}
}
}
for(i=0;i<n2;i++)
{
for(j=n1;j<2*n2;j++)
{
A2[i][j-n2]=a[i][j] ;
}
}
/*for(i=0;i<n2;i++)
{
cout << endl;
for(j=0; j<n2 ; j++ )
{
cout << A2[i][j] << " ";
}
}*/
/*********************************************************************************************
设置节点电压初值：U[i] ， rad[i]
*********************************************************************************************/
cout << "请输入平衡节点电压幅值:" ;
cin >> U[0];
cout << "请输入其余" << m << "个节点电压幅值初值：" << endl;/********************************/
for( i=1 ; i<=m ; i++ ) /* 1.06 1 1 1 1 */
{
cin >> U[i];
}
cout << "请输入" << n << "个节点电压相位角初值:" << endl;
for( i=0 ; i<n ; i++ ) /* 0 0 0 0 0 */
{
cin >> rad[i];
}
/*********************************************************************************************
输入注入功率
*********************************************************************************************/
cout << "请输入" << n-1 << "个节点注入有功初值：(平衡节点除外）\n";
for( i=1 ; i<n ; i++ )
{
cin >> P[i]; // 0.20 -0.45 -0.40 -0.60
}
cout << "请输入" << n-1 << "个节点注入无功初值：(平衡节点除外）\n";
for( i=1 ; i<n ; i++)
{
cin >> Q[i]; // 0.20 -0.15 -0.05 -0.10
}
/*********************************************************************************************
开始迭代
*********************************************************************************************/
for( k=0 ; k<=6 ; k++ )
{
/**********************************************************************************
计算有功不平衡量cp[i]
***********************************************************************************/
for( i=1 ; i<n ; i++ )
{
double s1=0;
for( j=0 ; j<n ; j++ )
{
s1+=U[j]*(G[i][j]*cos(rad[i]-rad[j])+B[i][j]*sin(rad[i]-rad[j]));
}
cp[i]=P[i]-U[i]*s1;
//cout << cp[i] << " " ;
}
//cout << endl;
/********************************************************************************
计算电压相位角的变量c_rad[i]和相位角的新值rad[i]
*******************************************************************************/
{
double a1[N][N]; //系数矩阵
double b[N]; //常数
double x[N];
double c_rad[N];
for( i=0 ; i<n-1 ; i++ )
{
for( j=0 ; j<n-1 ; j++ )
a1[i][j]=A1[i][j];
}
for( i=0 ; i<n-1 ; i++ )
{
b[i]=cp[i+1]/U[i+1]; //cout << b[i] << " ";
}
for( i=0 ; i<n-1 ; i++ )
{
double sum=0;
for( j=0 ; j<n-1 ; j++ )
{
sum+=a1[i][j]*b[j];
}
x[i]=sum;
}
for( i=0; i<n-1 ; i++ )
{
c_rad[i]=-x[i]/U[i+1];
}
for( i=0 ; i<n-1 ; i++ )
{
rad[i+1]=rad[i+1]+c_rad[i]; //cout << rad[i+1] << " ";
}
}
//cout << endl;
/**********************************************************************************
计算无功不平衡量cq[i]
***********************************************************************************/
for( i=1 ; i<n ; i++ )
{
double s1=0;
for( j=0 ; j<n ; j++ )
{
s1+=U[j]*(G[i][j]*sin(rad[i]-rad[j])-B[i][j]*cos(rad[i]-rad[j]));
}
cq[i]=Q[i]-U[i]*s1;
//cout << cq[i] << " " ;
}
//cout << endl;
/********************************************************************************
计算电压大小的变量CU[i]和电压的新值U[i]
*******************************************************************************/
{
double a1[N][N]; //系数矩阵
double b[N]; //常数
double x[N];
double CU[N];
for( i=0 ; i<m ; i++ )
{
for( j=0 ; j<m ; j++ )
a1[i][j]=A2[i][j];
}
for( i=0 ; i<m ; i++ )
{
b[i]=cq[i+1]/U[i+1]; //cout << b[i] << " ";
}
for( i=0 ; i<n-1 ; i++ )
{
double sum=0;
for( j=0 ; j<n-1 ; j++ )
{
sum+=a1[i][j]*b[j];
}
x[i]=sum;
}
for( i=0; i<n-1 ; i++ )
{
CU[i]=-x[i]/U[i+1]; //cout << CU[i] << " ";
}
for( i=0 ; i<n-1 ; i++ )
{
U[i+1]=U[i+1]+CU[i]; //cout << U[i+1] << " ";
}
}
/*********************************************************************************************
输出电压大小和相位角大小
*********************************************************************************************/
cout << endl;
for( i=1 ; i<n ; i++ )
{
cout << rad[i] << " " << U[i] << " ";
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -