【DOC】-潮流计算程序至雅克比矩阵
潮流计算程序至雅克比矩阵
%本程序的功能是用牛顿——拉夫逊法进行潮流计算
% X1矩阵:1、支路首端号;2、末端号;3、支路阻抗;4、支路对地电纳
% 5、支路的变比;6、支路首端处于K侧为1，1侧为0
% X2矩阵:1、该节点发电机功率;2、该节点负荷功率;3、节点电压初始值
% 4、PV节点电压V的给定值;5、初相角给定值 6、节点所接的无功补偿设备的容量
% 7、节点分类标号:1为PQ节点;2为PV节点;3为平衡节点;
clear;
n=6; %input('请输入节点数:n=')
nl=7; %input('请输入支路数:nl=')
n2=5; %input('请输入已知条件中含有有功功率的节点数:n2=')
n3=4; %input('请输入已知条件中含有无功功率的节点数:n3=')
isb=6; %input('请输入平衡母线节点号:isb=')
pr=0.00001; %input('请输入误差精度:pr=')
X1=[1 2 0.000+0.300i 0 1.025 1;
1 4 0.097+0.407i 0 1 0;
1 6 0.123+0.518i 0 1 0;
2 5 0.282+0.640i 0 1 0;
3 5 0.723+1.050i 0 1 0;
3 4 0.000+0.133i 0 1.100 0;
4 6 0.080+0.370i 0 1 0]
%input('请输入由支路参数形成的矩阵:X1=')
X2=[0 50.0+5.0i 1.000 0 0 0 1;
0 30.0+18.0i 1.000 0 0 0 1;
0 55.0+13.0i 1.000 0 0 0 1;
0 0 1.000 0 0 0 1;
50.1 0 1.000 1.100 0 0 2;
0 0 1.000 1.050 0 0 3]
%input('请输入各节点参数形成的矩阵:X2=')
X3=[1 0;
2 0;
3 0;
4 0;
5 0;
6 0]
%input('请输入由节点号及其对地阻抗形成的矩阵:X3=')
Y=zeros(n);u=zeros(1,n);v=zeros(1,n);delt=zeros(1,n);JJ=zeros(n2
+n3);%设置各个参数矩阵的形式
for i=1:n
if X3(i,2)~=0;
a=X3(i,1);
Y(a,a)=1./X3(i,2);
end
end
for i=1:nl
if X1(i,6)==0
a=X1(i,1);b=X1(i,2);
else
a=X1(i,2);b=X1(i,1);
end
Y(a,b)=Y(a,b)-1./(X1(i,3)*X1(i,5));
Y(b,a)=Y(a,b);
Y(b,b)=Y(b,b)+1./(X1(i,3)*X1(i,5)^2)+X1(i,4)./2;
Y(a,a)=Y(a,a)+1./X1(i,3)+X1(i,4)./2;
end %求导纳矩阵
disp('导纳矩阵 Y=')
disp(Y)
G=real(Y);B=imag(Y); %将导纳矩阵的实部和虚部整理出来
for i=1:n
u(i)=X2(i,3);
v(i)=X2(i,4);
delt(i)=X2(i,5);
end %根据X2矩阵得出各个节点电压的初始给定值和相角 for i=1:n
S(i)=X2(i,1)-X2(i,2);
B(i,i)=B(i,i)+X2(i,6);
end %计算各个节点的发电机功率与负荷功率的差值，以及检查各个节点是否有无功补偿装置
P=real(S); Q=imag(S); %将各个节点的发电机功率与负荷功率的差值的实部和虚部整理出来
disp(P);
disp(Q);
k=0; %设置迭代次数初始值
pr=1; %给定一个误差精度
while pr>0.00001
for a=1:n2
for b=1:n
pt(a)=u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-
delt(b)));
end
pp(a)=P(a)-sum(pt);
end
for a=1:n3
for b=1:n
qt(a)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-
delt(b)));
end
qq(a)=Q(a)-sum(qt);
end %求得各个节点的有功功率和无功功率的误差
disp(qq);
disp(pp);
for a=1:n2
for b=1:n
hh(b)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-
delt(b)));
end
w1(a)=sum(hh);
for a=1:n3
for b=1:n
nn(b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));
kk(b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));
ll(b)=-u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-
B(a,b)*cos(delt(a)-delt(b)));
end
w2(a)=sum(nn);
w3(a)=sum(kk);
w4(a)=sum(ll);
end %求式(3-33)中的包括j=i在内所计算的结果
for a=1:n2
for b=1:n2
if a==b
H(a,b)=w1(a)-u(a)^2*(G(a,a)*sin(delt(a)-delt(a))-
B(a,a)*cos(delt(a)-delt(a)));
else
H(a,b)=-u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-
B(a,b)*cos(delt(a)-delt(b)));
end
end
end %计算H单元的元素
for a=1:n2
for b=1:n3
if a==b
N(a,b)=w2(a)+2*u(a)^2*G(a,a)+u(a)^2*(G(a,a)*cos(delt(a)-
delt(a))+B(a,a)*sin(delt(a)-delt(a))); else
N(a,b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));
end
end
end %计算N单元的元素
for a=1:n3
for b=1:n2
if a==b
K(a,b)=w3(a)+u(a)^2*(G(a,a)*cos(delt(a)-delt(a))-
B(a,a)*sin(delt(a)-delt(a)));
else
K(a,b)=u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));
end
end
end %计算K单元的元素
for a=1:n3
for b=1:n3
if a==b
L(a,b)=w4(a)+2*u(a)^2*G(a,a)+u(a)^2*(G(a,a)*sin(delt(a)-delt(a))+B(a,a)*cos(delt(a)-delt(a))); else
L(a,b)=-u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-
B(a,b)*cos(delt(a)-delt(b)));
end
end
end %计算L单元的元素
JJ=[H(1,1) H(1,2) H(1,3) H(1,4) H(1,5) N(1,1) N(1,2) N(1,3) N(1,4);
H(2,1) H(2,2) H(2,3) H(2,4) H(2,5) N(2,1) N(2,2) N(2,3) N(2,4);
H(3,1) H(3,2) H(3,3) H(3,4) H(3,5) N(3,1) N(3,2) N(3,3) N(3,4);
H(4,1) H(4,2) H(4,3) H(4,4) H(4,5) N(4,1) N(4,2) N(4,3) N(4,4);
H(5,1) H(5,2) H(5,3) H(5,4) H(5,5) N(5,1) N(5,2) N(5,3) N(5,4);
K(1,1) K(1,2) K(1,3) K(1,4) K(1,5) L(1,1) L(1,2) L(1,3) L(1,4);
K(2,1) K(2,2) K(2,3) K(2,4) K(2,5) L(2,1) L(2,2) L(2,3) L(2,4);
K(3,1) K(3,2) K(3,3) K(3,4) K(3,5) L(3,1) L(3,2) L(3,3) L(3,4);
K(4,1) K(4,2) K(4,3) K(4,4) K(4,5) L(4,1) L(4,2) L(4,3) L(4,4)] %排列雅可比矩阵的位置
disp(JJ)
for a=1:n2
DP(a,1)=pp(a);
end
for a=1:4
DP(a+n2,1)=qq(a);
end %放置有功功率和无功功率误差的位置
uu=-inv(JJ)*DP; %inv表示JJ的逆矩阵
pr=max(abs(uu)); %检验是否满足要求 for a=1:n2
delt(a)=delt(a)+uu(a);
end
for a=1:n3
u(a)=u(a)+uu(a+n2);
end %对方程式进行修正 k=k+1;
k-1,delt',u'
end
因篇幅问题不能全部显示，请点此查看更多更全内容