程序的主要代码：
n=input('请输入节点数n=');
na=input('请输入支路数na=');
isb=input('请输入平衡节点母线号isb=');
jd=input('请输入误差精度jd=');
B1=input('请输入由支路参数形成的矩阵B1=');
B2=input('请输入由节点参数形成的矩阵B2=');
L=input('请输入由节点号及其对地阻抗形成的矩阵L=');
nb=input('请输入P-Q节点数nb=');
Y=zeros(n);Z=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);
O=zeros(1,n);
for i=1:na
if B1(i,6)==0
a=B1(i,1);b=B1(i,2);
else a=B1(i,2);b=B1(i,1);
end
Y(a,b)=Y(a,b)-1./(B1(i,3)*B1(i,5));
Z(a,b)=Z(a,b)-1./(B1(i,3));
Y(b,a)=Y(a,b);
Z(b,a)=Z(a,b);
Y(b,b)=Y(b,b)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
Z(b,b)=Z(b,b)+1./(B1(i,3));
Y(a,a)=Y(a,a)+1./(B1(i,3))+B1(i,4)./2;
Z(a,a)=Z(a,a)+1./(B1(i,3));
end
G=real(Y);B=imag(Z);CI=imag(Y);
for i=1:n
S(i)=B2(i,1)-B2(i,2);
CI(i,i)=CI(i,i)+B2(i,5);
end
P=real(S);Q=imag(S);
for i=1:n
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
end
for i=1:n
if B2(i,6)==2
V(i)=sqrt(e(i)^2+f(i)^2);
O(i)=atan(f(i)./e(i));
end
end
for i=2:n
if i==n
B(i,i)=1./B(i,i);
else
IT1=i+1;
for j1=IT1:n
B(i,j1)=B(i,j1)./B(i,i);
end
B(i,i)=1./B(i,i);
for k=i+1:n
for j1=i+1:n
B(k,j1)=B(k,j1)-B(k,i)*B(i,j1);
end
end
end
end
a=0;b=0;
for i=1:n
if B2(i,6)==2
a=a+1;k=0;
for j1=1:n
if B2(j1,6)==2
k=k+1;
A(a,k)=CI(i,j1);
end
end
end
end
for i=1:nb
if i==na
A(i,i)=1./A(i,i);
else
k=i+1;
for j1=k:nb
A(i,j1)=A(i,j1)./A(i,i);
end
A(i,i)=1./A(i,i);
for k=i+1:nb
for j1=i+1:nb
A(k,j1)=A(k,j1)-A(k,i)*A(i,j1);
end
end
end
end
NT2=1;NT1=0;kp=1;kq=1;K=1;NCT=0;NT3=1;
while NT2~=0|NT3~=0
NT2=0;NT3=0;
for i=1:n
if i~=isb
C(i)=0;
for k=1:n
C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+CI(i,k)*sin(O(i)-O(k)));
end
CP1(i)=P(i)-V(i)*C(i);
CP(i)=CP1(i)./V(i);
NCT=abs(CP1(i));
if NCT>=jd
NT2=NT2+1;
end
end
end
Np(k)=NT2;
if NT2~=0
for i=2:n
CP(i)=B(i,i)*CP(i);
if i~=n
IT1=i+1;
for k=IT1:n
CP(k)=CP(k)-B(k,i)*CP(i);
end
else
for LZ=3:i
L=i+3-LZ;
NC4=L-1;
for MZ=2:NC4
I=NC4+2-MZ;
CP(I)=CP(I)-B(I,L)*CP(L);
end
end
end
end
for i=2:n
O(i)=O(i)-CP(i);
end
kq=1;L=0;
for i=1:n
if B2(i,6)==2
C(i)=0;L=L+1;
for k=1:n
C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-CI(i,k)*cos(O(i)-O(k)));
end
DQ1(i)=Q(i)-V(i)*C(i);
DQ(L)=DQ1(i)./V(i);
NCT=abs(DQ1(i));
if NCT>=jd
NT3=NT3+1;
end
end
end
else
kp=0;
if kq~=0;
L=0;
for i=1:n
if B2(i,6)==2
C(i)=0;L=L+1;
for k=1:n
C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-CI(i,k)*cos(O(i)-O(k)));
end
DQ1(i)=Q(i)-V(i)*C(i);
DQ(L)=DQ1(i)./V(i);
NCT=abs(DQ1(i));
end
end
end
end
Nq(K)=NT3;
if NT3~=0
L=0;
for i=1:nb
DQ(i)=A(i,i)*DQ(i);
if i==nb
for LZ=2:i
L=i+2-LZ;
NC4=L-1;
for MZ=1:NC4
I=NC4+1-MZ;
DQ(I)=DQ(I)-A(I,L)*DQ(L);
end
end
else
IT1=i+1;
for k=IT1:nb
DQ(k)=DQ(k)-A(k,i)*DQ(i);
end
end
end
L=0;
for i=1:n
if B2(i,6)==2
L=L+1;
V(i)=V(i)-DQ(L);
end
end
kp=1;
K=K+1;
else
kq=0;
if kp~=0
K=K+1;
end
end
for i=1:n
Dp(K-1,i)=V(i);
end
end
disp('迭代次数');
disp(K);
disp('每次没有达到精度要求的有功功率个数为');
disp(Np);
disp('每次没有达到精度要求的无功功率个数为');
disp(Nq);
for k=1:n
E(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;
O(k)=O(k)*180./pi;
end
disp('各节点的实际电压标幺值E(节点号从小到大排列):');
disp(E);
disp('各节点的电压大小V为(节点号从小到大排列):');
disp(V);
disp('各节点的电压相角O为(节点号从小到大排列):');
disp(O);
for a=1:n
C(a)=0;
for b=1:n
C(a)=C(a)+conj(Y(a,b))*conj(E(b));
end
S(a)=E(a)*C(a);
end
disp('各节点的功率S为(节点号从小到大排列):');
disp(S);
disp('各条支路的首端功率Si为(顺序同您输入B1时一样):');
for i=1:na
if B1(i,6)==0
a=B1(i,1);b=B1(i,2);
else
a=B1(i,2);b=B1(i,1);
end
Si(a,b)=E(a)*(conj(E(a))*conj(B1(i,4)./2)+(conj(E(a)*B1(i,5))-conj(E(b)))*conj(1./(B1(i,3)*B1(i,5))));
disp(Si(a,b));
end
disp('各条支路的末端功率Sj为(顺序同您输入B1时一样):');
for i=1:na
if B1(i,6)==0
a=B1(i,1);b=B1(i,2);
else
a=B1(i,2);b=B1(i,1);
end
Sj(b,a)=E(b)*(conj(E(b))*conj(B1(i,4)./2)+(conj(E(b)./B1(i,5))-conj(E(a)))*conj(1./(B1(i,3)*B1(i,5))));
disp(Sj(b,a));
end
disp('各条支路的功率损耗DS为(顺序同您输入B1时一样)::');
for i=1:na
if B1(i,6)==0
a=B1(i,1);b=B1(i,2);
else
a=B1(i,2);b=B1(i,1);
end
DS(i)=Si(a,b)+Sj(b,a);
disp(DS(i));
end
for i=1:K
Cs(i)=i;
for j=1:n
Dp(K,j)=Dp(K-1,j);
end
end
disp('以下是每次迭代后各节点的电压值（如图所示）');
plot(Cs,Dp),xlabel('迭代次数'),ylabel('电压'),title('电压迭代次数曲线');下载本文
