计算导纳矩阵.docx
《计算导纳矩阵.docx》由会员分享，可在线阅读，更多相关《计算导纳矩阵.docx（13页珍藏版）》请在冰豆网上搜索。
计算导纳矩阵
//****************计算导纳矩阵*******************
G[1][1]=z12r/（z12r*z12r+z12m*z12m）+k*k*z13r/（z13r*z13r+z13m*z13m）+z14r/（z14r*z14r+z14m*z14m）;
B[1][1]=-z12m/（z12r*z12r+z12m*z12m）-k*k*z13m/（z13r*z13r+z13m*z13m）-z14m/（z14r*z14r+
z14m*z14m）+y140+y120;
G[2][2]=z12r/（z12r*z12r+z12m*z12m）+z24r/（z24r*z24r+z24m*z24m）;
B[2][2]=-z12m/（z12r*z12r+z12m*z12m）-z24m/（z24r*z24r+z24m*z24m）+y240+y120;
G[3][3]=z13r/（z13r*z13r+z13m*z13m）;
B[3][3]=-z13m/（z13r*z13r+z13m*z13m）;
G[4][4]=z14r/（z14r*z14r+z14m*z14m）+z24r/（z24r*z24r+z24m*z24m;
B[4][4]=-z14m/（z14r*z14r+z14m*z14m）-z24m/（z24r*z24r+z24m*z24m）+y240+y140;
G[1][2]=G[2][1]=-z12r/（z12r*z12r+z12m*z12m）;
B[1][2]=B[2][1]=z12m/（z12r*z12r+z12m*z12m）;
G[1][3]=G[3][1]=-k*z13r/（z13r*z13r+z13m*z13m）;
B[1][3]=B[3][1]=k*z13m/（z13r*z13r+z13m*z13m）;
G[1][4]=G[4][1]=-z14r/（z14r*z14r+z14m*z14m）;
B[1][4]=B[4][1]=z14m/（z14r*z14r+z14m*z14m）;
G[2][3]=G[3][2]=0.0;
B[2][3]=B[3][2]=0.0;
G[2][4]=G[4][2]=-z24r/（z24r*z24r+z24m*z24m）;
B[2][4]=B[4][2]=z24m/（z24r*z24r+z24m*z24m）;
G[3][4]=G[4][3]=0.0;
B[3][4]=B[4][3]=0.0;
for（i=1;i<5;i++）
{for（j=1;j<5;j++）
{printf（"%f+%fj",G[i][j],B[i][j]）;
printf（" "）;
}
printf（"\n"）;//形成节点导纳矩阵
//计算各节点不平衡量
loop1:
printf（"迭代次数k1=%d\n",k1）;
for （i=1;i<5;i++）
{float a=0,b=0;
for（j=1;j<5;j++）
{a+=G[i][j]*e[j]-B[i][j]*f[j];
b+=G[i][j]*f[j]+B[i][j]*e[j];
}
P[i]=Ps[i]-（e[i]*a+f[i]*b）;//计算有功功率的增量
Q[i]=Qs[i]-（f[i]*a-e[i]*b）;//计算无功功率的增量
V32=V3s*V3s-e[3]*e[3];
}
printf（"有功功率增量P[1]=%f",P[1]）;
printf（" ,"）;
printf（"有功功率增量P[2]=%f",P[2]）;
printf（" ,"）;
printf（"有功功率增量P[3]=%f",P[3]）;
printf（"无功功率增量Q[1]=%f",Q[1]）;
printf（" ,"）;
printf（"无功功率增量Q[2]=%f",Q[2]）;
printf（" ,"）;
printf（"电压增量V32=%f",V32）;
printf（"\n"）;
//****形成雅克比矩阵**********************
for（j=1;j<4;j++）
{if（1==j）
{float c=0,d=0;
int m;
for（m=1;m<5;m++）
{c+=G[1][m]*e[m]-B[1][m]*f[m];
d+=G[1][m]*f[m]+B[1][m]*e[m];
}
J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1];
J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1];
}
else
{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1];
J[1*N][j*N]=G[1][j]*e[1]-B[1][j]*f[1];
J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1];
}
}
//********计算修正方程*************
for（i=1;i{L[i][i]=1;
}
L[i][1]=J[i][1]/U[1][1];
}
for（n=2;n{
for（j=n;j{
sigma1=0;
for（s=0;s<=n-1;s++）
sigma1+=L[n][s]*U[s][j];
U[n][j]=J[n][j]-sigma1;
}
for（i=n;i{
sigma2=0;
for（s=0;s<=n-1;s++）
sigma2+=L[i][s]*U[s][n];
L[i][n]=（J[i][n]-sigma2）/U[n][n];
}
}
b[1]=P[1];
b[2]=Q[1];
b[3]=P[2];
b[4]=Q[2];
b[5]=P[3];
b[6]=V32;
for（i=1;i{
sigma1=0;
for（n=1;n<=i-1;n++）
for（i=1;i{
U[1][i]=J[1][i];
sigma1+=L[i][n]*y[n];
y[i]=b[i]-sigma1;
}
for（i=M-1;i>=1;i--）
{
sigma2=0;
for（n=i+1;nsigma2+=U[i][n]*x[n];
x[i]=（y[i]-sigma2）/U[i][i];
}
xe[1]=-x[1];xe[2]=-x[3];xe[3]=-x[5];
xf[1]=-x[2];xf[2]=-x[4];xf[3]=-x[6];
printf（"节点电压:
\n"）;
for（i=1;i<4;i++）
{e[i]+=xe[i];
f[i]+=xf[i];
}
for（i=1;i<4;i++）
{printf（"e[%d]=",i）;
printf（"%f",e[i]）;
printf（" ,"）;
}
for（i=1;i<4;i++）
{printf（"f[%d]=",i）;
printf（"%f",f[i]）;
printf（" ,"）;
}
%给定节点电压初值
e=[111.11.05];
f=[0000];
%给定PQ节点有用和无功功率，以及PV节点有功和电压
P1=-0.3;
Q1=-0.18;
P2=-0.55;
Q2=-0.13;
P3=0.5;
V3=1.1;
%手动输入节点导纳矩阵
G=[1.042093,-0.588235,0,-0.453858;
-0.588235,1.069005,0,-0.480769;
0,0,0,0;
-0.453858,-0.480769,0,0.934627
];
B=[-8.242876,2.352941,3.666667,1.891074;
2.352941,-4.727377,0,2.403846;
3.666667,0,-3.333333,0;
1.891074,2.403846,0,-4.261590
];
maxP=0;
maxQ=0;
maxV=0;
I=[0,0;
0,0;
0,0;
0,0
];
forv=1:
15%迭代次数
forn=1:
4
I（1,1）=I（1,1）+G（1,n）*e（n）-B（1,n）*f（n）;
I（1,2）=I（1,2）+G（1,n）*f（n）+B（1,n）*e（n）;
end
forn=1:
4
I（2,1）=I（2,1）+G（2,n）*e（n）-B（2,n）*f（n）;
I（2,2）=I（2,2）+G（2,n）*f（n）+B（2,n）*e（n）;
end
forn=1:
4
I（3,1）=I（3,1）+G（3,n）*e（n）-B（3,n）*f（n）;
I（3,2）=I（3,2）+G（3,n）*f（n）+B（3,n）*e（n）;
end
forn=1:
4
I（4,1）=I（4,1）+e（n）*e（n）+f（n）*f（n）;
I（4,2）=0;
end
H=[];
N=[];
M=[];
L=[];
R=[];
S=[];
J=[];
%求不平衡量
P1=-0.30-e
（1）*I（1,1）-f
（1）*I（1,2）;
Q1=-0.18-f
（1）*I（1,1）+e
（1）*I（1,2）;
P2=-0.55-e
（2）*I（2,1）-f
（2）*I（2,2）;
Q2=-0.13-f
（2）*I（2,1）+e
（2）*I（2,2）;
P3=0.50-e（3）*I（3,1）-f（3）*I（3,2）;
V3=1.10^2-I（4,1）;
%分块计算雅各比矩阵元素
form=1:
3
forn=1:
3
if（m==n）
H（m,m）=B（m,m）*e（m）-G（m,m）*f（m）-I（m,2）;
N（m,m）=-G（m,m）*e（m）-B（m,m）*f（m）-I（m,1）;
M（m,m）=G（m,m）*e（m）+B（m,m）*f（m）-I（m,1）;
L（m,m）=B（m,m）*e（m）-G（m,m）*f（m）+I（m,2）;
R（m,m）=-2*f（m）;
S（m,m）=-2*e（m）;
else
H（m,n）=B（m,n）*e（m）-G（m,n）*f（m）;
N（m,n）=-G（m,n）*e（m）-B（m,n）*f（m）;
M（m,n）=-N（m,n）;
L（m,n）=H（m,n）;
R（m,m）=0;
S（m,m）=0;
end
end
end
%确定雅克比矩阵各元素
J=[H（
升级会员 