牛顿拉夫逊潮流计算[整理版]

牛顿拉夫逊潮流计算[整理版] 牛顿拉夫逊潮流计算

//整个程序为:

//牛拉法解潮流程序

#include

#include

#define N 4 //节点数

#define n_PQ 2 //PQ节点数

#define n_PV 1 //PV节点数

#define n_br 5 //串联支路数

void main()

{

void disp_matrix(float *disp_p,int disp_m,int disp_n); //矩阵显示函数

float Us[2*N]={1.0,0,1.0,0,1.05,0,1.05,0}; //电压初值

float Ps[N]={0,-0.5,0.2}; //有功初值

float Qs[N]={0,-0.3}; //无功初值

float G[N][N],B[N][N]; //各几点电导电纳

struct

//阻抗参数

{

int nl; //左节点

int nr; //右节点

float R; //串联电阻值

float X; //串联电抗值

float Bl; //左节点并联电导

float Br; //右节点并联电纳

}ydata[n_br]={

{1,2,0,0.1880,-0.6815,0.6040},

{1,3,0.1302,0.2479,0.0129,0.0129},

{1,4,0.1736,0.3306,0.0172,0.0172},

{3,4,0.2603,0.4959,0.0259,0.0259},

{2,2,0,0.05,0,0}

};

float Z2; //Z^2=R^2+X^2 各串

联阻抗值的平方

float e[N],f[N],dfe[2*(N-1)]; //e,f存储电压的x轴分量和y轴分量,dfe 存储电压修正值

float mid1[N],mid2[N],dS[2*(N-1)]; //mid1、mid2存储计算雅克比行列式对角线元素的中间值,dS存储

PQU的不平衡量

float Jacob[2*(N-1)][2*(N-1)],inv_J[2*(N-1)][2*(N-1)];

//雅克比行列式

float dPQU=1.0; //PQU不平衡量最大值

int kk=0; //迭代次数

int i,j,k;

float t;

float Pij[n_br]; //存储线路i->j的有功

float Qij[n_br]; //存储线路i->j的无功

float Pji[n_br]; //存储线路j->i的有功

float Qji[n_br]; //存储线路j->i的无功

float dPij[n_br]; //存储线路i->j的有功损耗

float dQij[n_br]; //存储线路i->j的无功损耗

float AA,BB,CC,DD; //存储线路潮流计算时的中间

值

//形成导纳矩阵

--------------------------------------------------------------------------------------------------

for(i=0;i

for(j=0;j

{

G[i][j]=0;

B[i][j]=0;

}

for(i=0;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-1][ydata[i].nr-1]=(-ydata[i].R)/Z2; B[ydata[i].nl-1][ydata[i].nr-1]=ydata[i].X/Z2;

G[ydata[i].nr-1][ydata[i].nl-1]=(-ydata[i].R)/Z2; B[ydata[i].nr-1][ydata[i].nl-1]=ydata[i].X/Z2;

//对角元素

G[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].R/Z2;

G[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].R/Z2;

B[ydata[i].nl-1][ydata[i].nl-1]+=(-ydata[i].X/Z2); B[ydata[i].nr-1][ydata[i].nr-1]+=(-ydata[i].X/Z2); //并联导纳等效导纳值

B[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].Bl;

B[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].Br;

}

else

{

G[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].R;

B[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].X;

}

}

printf("G=\n");

disp_matrix(*G,N,N);

printf("B=\n");

disp_matrix(*B,N,N);

//分离e,f