Matlab实现潮流计算程序你可以通过这个链接引用该篇文章:http://lntianxia.bokee.com/viewdiary.16330339.html
程序代码如下:
111111.
%读入数据
clc
clear
filename='123.txt';
a=textread(filename)
n=a(1,1);
pinghengjd=a(1,2);
phjddianya=a(1,3);
jingdu=a(1,4);
b=zeros(1,9);
j1=0;
[m1,n1]=size(a);
for i1=1:m1
if a(i1,1)==0
j1=j1+1;
b(j1)=i1;
end
end
b;
%矩阵分块
a1=a(b(1)+1:b(2)-b(1)+1,1:n1);
a2=a(b(2)+1:b(3)-1,1:n1);
a3=a(b(3)+1:b(4)-1,1:n1);
a4=a(b(4)+1:b(5)-1,1:n1);
a5=a(b(5)+1:b(6)-1,1:n1);
%设置初值
vcz=1;
dcz=0;
kmax=20;
k1=0;
%求节点导纳矩阵
a11=zeros(4,6);
for i0=1:3
for j0=1:6
a11(i0,j0)=a1(i0,j0);
a11(4,j0)=a2(1,j0);
end
end
a11;
linei=a11(1:4,2);
linej=a11(1:4,3);
liner=a11(1:4,4);
linex=a11(1:4,5);
lineb=a11(1:4,6);
branchi=0;
branchj=0;
branchb=0;
G=zeros(4,4);
B=zeros(4,4);
for k=1:4
i2=linei(k,1);
j2=linej(k,1);
r=liner(k,1);
x=linex(k,1);
b=0;
GIJ=r/(r*r+x*x);
BIJ=-x/(r*r+x*x);
if k>=4 & lineb(k)~=0
k0=lineb(k);
G(i2,j2)=-GIJ/k0;
G(j2,i2)=G(i2,j2);
B(i2,j2)=-BIJ/k0;
B(j2,i2)=B(i2,j2);
G(i2,i2)=G(i2,i2)+GIJ/k0/k0;
B(i2,i2)=B(i2,i2)+BIJ/k0/k0;
else
G(j2,i2)=-GIJ;
G(i2,j2)=G(j2,i2);
B(j2,i2)=-BIJ;
B(i2,j2)=B(j2,i2);
G(i2,i2)=G(i2,i2)+GIJ;
b=lineb(k);
B(i2,i2)=B(i2,i2)+BIJ+b;
end
G(j2,j2)=G(j2,j2)+GIJ;
B(j2,j2)=B(j2,j2)+BIJ+b;
end
G;
B;
B=B.*i;
Yf=G+B
Y=abs(Yf);
alf=angle(Yf);
%赋Jacobian矩阵参数
P=zeros(n,1);
Q=zeros(n,1);
Pd=zeros(1,n);
Qd=zeros(1,n);
dP=zeros(1,n);
dQ=zeros(1,n);
PG=a4(:,3);
PD=a4(:,5);
QG=a4(:,4);
QD=a4(:,6);
i8=a4(:,2);
for j8=1:length(i8)
P(i8(j8))=PG(i8(j8))-PD(i8(j8));
Q(i8(j8))=QG(i8(j8))-QD(i8(j8));
end
delt=zeros(n,1);
V=ones(n,1);
V(3)=1.10;
V(4)=1.05;
ddelt=zeros(n,1);
dV=zeros(n,1);
A=zeros(2*n,2*n);
B=zeros(2*n,1);
Jacobian=Jaco(V,delt,n,Y,alf)
%求取矩阵功率
for j5=1:kmax
disp(['第' int2str(j5) '次计算结果'])
if k>=kmax
break
end
for i10=1:4
Pd(i10)=0;
Qd(i10)=0;
for j10=1:n
Pd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(delt(i10)-delt(j10)-alf(i10,j10));
Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(delt(i10)-delt(j10)-alf(i10,j10));
end
end
for i4=1:3
dP(i4)=P(i4)-Pd(i4);
end
for j4=1:2
dQ(j4)=Q(j4)-Qd(j4);
end
A=Jaco(V,delt,n,Y,alf)
for i14=1:n
B(i14*2-1)=-dP(i14);
B(i14*2)=-dQ(i14);
end
if max(abs(B))>jingdu
X=A\B;
for i16=1:n
ddelt(i16)=X(2*i16-1);
dV(i16)=X(2*i16)*V(i16);
end
V=V+dV
delt=delt+ddelt
else
break
end
disp('----------------')
end
%流氓算法
% for ii=1:2
% V(ii)=V(ii)+dV(ii);
% end
% V
222222.
function A=Jaco(V,delt,n,Y,alf)
%计算Jacobian矩阵
for i7=1:n
Hd1(i7)=0;
Jd1(i7)=0;
for j7=1:n
Hd1(i7)=Hd1(i7)+V(i7)*Y(i7,j7)*V(j7)*sin(delt(i7)-delt(j7)-alf(i7,j7));
Jd1(i7)=Jd1(i7)+V(i7)*Y(i7,j7)*V(j7)*cos(delt(i7)-delt(j7)-alf(i7,j7));
end
end
for i6=1:n
for j6=1:n
if i6~=j6
H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-delt(j6)-alf(i6,j6));
N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));
J(i6,j6)=-N(i6,j6);
L(i6,j6)=H(i6,j6);
else
H(i6,i6)=Hd1(i6)-V(i6)*Y(i6,i6)*V(i6)*sin(delt(i6)-delt(j6)-alf(i6,j6));
J(i6,j6)=-Jd1(i6)+V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));
N(i6,j6)=-Jd1(i6)-V(i6)*Y(i6,i6)*V(i6)*cos(alf(i6,i6));
L(i6,i6)=-Hd1(i6)+V(i6)*Y(i6,i6)*V(i6)*sin(alf(i6,i6));
end
end
end
%修正Jacobian矩阵
for j9=3
for i9=1:n
N(i9,j9)=0;
L(i9,j9)=0;
J(j9,i9)=0;
L(j9,i9)=0;
end
end
L(j9,j9)=1;
for j9=4
for i9=1:n
H(i9,j9)=0;
N(i9,j9)=0;
J(i9,j9)=0;
L(i9,j9)=0;
H(j9,i9)=0;
N(j9,i9)=0;
J(j9,i9)=0;
L(j9,i9)=0;
end
end
H(j9,j9)=1;
L(j9,j9)=1;
%Jaco=[H N;J L];
%Jaco=zeros(2*n,2*n);
for i11=1:n
for j11=1:n
Jaco(2*i11-1,2*j11-1)=H(i11,j11);
Jaco(2*i11-1,2*j11)=N(i11,j11);
Jaco(2*i11,2*j11-1)=J(i11,j11);
Jaco(2*i11,2*j11)=L(i11,j11);
end
end
A=Jaco;
33333.
数据:
4 4 1.05 0.00001
0
1 1 2 0.1 0.40 0.01528
2 1 4 0.12 0.50 0.01920
3 2 4 0.08 0.40 0.01413
0
1 1 3 0 0.3 0.90909091
0
0
1 1 0 0 0.30 0.18
2 2 0 0 0.55 0.13
3 3 0.5 0 0 0
0
1 3 1.10 0 0
0
