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

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

牛顿拉夫逊潮流计算

//整个程序为:

//牛拉法解潮流程序

#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;ifor(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

for(i=0;i{

e[i]=Us[2*i];

f[i]=Us[2*i+1];

}

//主程序

===========================================================

========================================== while(dPQU>0.00001)

{

//计算功率不平衡量

for(i=0;i{

mid1[i]=0;

mid2[i]=0;

for(j=0;j{

mid1[i]=mid1[i]+G[i][j]*e[j]-B[i][j]*f[j];

mid2[i]=mid2[i]+G[i][j]*f[j]+B[i][j]*e[j];

}

dS[2*i]=Ps[i]-(e[i]*mid1[i]+f[i]*mid2[i]);

if(idS[2*i+1]=Qs[i]-(f[i]*mid1[i]-e[i]*mid2[i]);

else

dS[2*i+1]=Us[2*i]*Us[2*i]-(e[i]*e[i]+f[i]*f[i]);

}

dPQU=0;

for(i=0;i<2*(N-1);i++)

{

if(dS[i]<0&&dPQU<-dS[i])

dPQU=-dS[i];

else if(dS[i]>0&&dPQUdPQU=dS[i];

}

if(dPQU>0.00001)

{

kk++;

//形成雅克比行列式

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

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

for(i=0;i<2*(N-1);i++)

for(j=0;j<2*(N-1);j++)

Jacob[i][j]=0;

for(j=0;j{

//求H,N

for(i=0;i{

if(i!=j)

{

Jacob[2*i][2*j]=B[i][j]*e[i]-G[i][j]*f[i];

Jacob[2*i][2*j+1]=-G[i][j]*e[i]-B[i][j]*f[i];

}

else

{

Jacob[2*i][2*i]=B[i][i]*e[i]-G[i][i]*f[i]-mid2[i];

Jacob[2*i][2*i+1]=-G[i][j]*e[i]-B[i][j]*f[i]-mid1[i];

}

}

//求J,L

for(i=0;i{

if(i!=j)

{

Jacob[2*i+1][2*j]=G[i][j]*e[i]+B[i][j]*f[i];

Jacob[2*i+1][2*j+1]=B[i][j]*e[i]-G[i][j]*f[i];

}

else

{

Jacob[2*i+1][2*i]=G[i][j]*e[i]+B[i][j]*f[i]-mid1[i];

Jacob[2*i+1][2*i+1]=B[i][j]*e[i]-G[i][j]*f[i]+mid2[i];

}

}

//求R,S

for(i=n_PQ;i{

if(i==j)

{

Jacob[2*i+1][2*i]=-2*f[i];

Jacob[2*i+1][2*i+1]=-2*e[i];

}

}

}

//雅克比行列式求逆

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

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

//初始化inv_J[N][N]

for(i=0;i<2*(N-1);i++)

for(j=0;j<2*(N-1);j++)

{

if(i!=j)

inv_J[i][j]=0;

else

inv_J[i][j]=1;

}

//将原矩阵化简为对角阵

for(i=0;i<2*(N-1);i++)

{

for(j=0;j<2*(N-1);j++)

{

if(i!=j)

{

t=Jacob[j][i]/Jacob[i][i];

for(k=0;k<2*(N-1);k++)

{

Jacob[j][k]-=Jacob[i][k]*t;

inv_J[j][k]-=inv_J[i][k]*t;

}

}

}

}

//原矩阵各对角元素化为1，画出逆矩阵

for(i=0;i<2*(N-1);i++)

if(Jacob[i][i]!=1)

{

t=Jacob[i][i];

for(j=0;j<2*(N-1);j++)

inv_J[i][j]=inv_J[i][j]/t;

}

//求电压修正值

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

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

for(i=0;i<2*(N-1);i++)

{

dfe[i]=0;

for(j=0;j<2*(N-1);j++)

dfe[i]-=inv_J[i][j]*dS[j];

}

for(i=0;i{

e[i]+=dfe[2*i+1];

f[i]+=dfe[2*i];

}

}

else

break;

}

//循环结束

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

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

//求平衡节点功率

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

mid1[N-1]=0;

mid2[N-1]=0;

for(j=0;j{

mid1[N-1]=mid1[N-1]+G[N-1][j]*e[j]-B[N-1][j]*f[j];

mid2[N-1]=mid2[N-1]+G[N-1][j]*f[j]+B[N-1][j]*e[j];

}

Ps[N-1]=e[N-1]*mid1[N-1]+f[N-1]*mid2[N-1];

Qs[N-1]=f[N-1]*mid1[N-1]-e[N-1]*mid2[N-1];

for(i=n_PQ;iQs[i]=f[i]*mid1[i]-e[i]*mid2[i];

//---------------------------------------------------------------------------------------------------------------------------

// 显示输出结果

printf(\"kk=%d\\n\

printf(\"P=\");

for(i=0;iprintf(\"%9.4f\

printf(\"\\nQ=\");

for(i=0;iprintf(\"%9.4f\

printf(\"\\ne=\");

for(i=0;iprintf(\"%9.4f\

printf(\"\\nf=\");

for(i=0;iprintf(\"%9.4f\

printf(\"\\n\");

//求线路上的潮流

//计算S[i][j]

for(i=0;i{

if(ydata[i].nl!=ydata[i].nr)

{

Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);

AA=-f[ydata[i].nl-1]*ydata[i].Bl+(e[ydata[i].nl-1]-e[ydata[i].nr-1])*ydata[i].R/Z2+(f[ydata[i].nl-1]-f[ydata[i].nr-1])*ydata[i].X/Z2;

BB=-e[ydata[i].nl-1]*ydata[i].Bl-(f[ydata[i].nl-1]-f[ydata[i].nr-1])*ydata[i].R/Z2+(e[ydata[i].nl-1]-e[ydata[i].nr-1])*ydata[i].X/Z2;

Pij[i]=e[ydata[i].nl-1]*AA-f[ydata[i].nl-1]*BB;

Qij[i]=e[ydata[i].nl-1]*BB+f[ydata[i].nl-1]*AA;

printf(\"S[%d][%d]=%9.4f+j%9.4f\\n\

dPij[i]=Pij[i]+Pji[i];

dQij[i]=Qij[i]+Qji[i];

}

}

printf(\"\\n\");

//计算S[j][i]

for(i=0;i{

if(ydata[i].nl!=ydata[i].nr)

{

Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);

CC=-f[ydata[i].nr-1]*ydata[i].Br+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].R/Z2+(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].X/Z2;

DD=-e[ydata[i].nr-1]*ydata[i].Br-(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].R/Z2+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].X/Z2;

Pji[i]=e[ydata[i].nr-1]*CC-f[ydata[i].nr-1]*DD;

Qji[i]=e[ydata[i].nr-1]*DD+f[ydata[i].nr-1]*CC;

printf(\"S[%d][%d]=%9.4f+j%9.4f\\n\

}

}

printf(\"\\n\");

//计算dS[i][j]

for(i=0;i{

if(ydata[i].nl!=ydata[i].nr)

{

dPij[i]=Pij[i]+Pji[i];

dQij[i]=Qij[i]+Qji[i];

printf(\"dS[%d][%d]=%9.4f+j%9.4f\\n\

dPij[i],dQij[i]);

}

}

printf(\"\\n\");

}//主程序结束

//矩阵显示函数

===========================================================

============================================= disp_matrix(float *disp_p,int disp_m,int disp_n)

void

{

int i,j;

for(i=0;i{

for(j=0;jprintf(\"%9.4f\

printf(\"\\n\");

}

printf(\"\\n\");

}