matlab实现电力系统潮流分析
n=input('请输入节点数:n=');
nl=input('请输入支路数:nl=');
isb=input('请输入平衡节点号:isb=');
pr=input('请输入误差精度:pr=');
B1=input('请输入由支路参数形成的矩阵:B1=');%变压器侧为1，否则为0
B2=input('请输入各节点参数形成的矩阵:B2=');
X=input('请输入由节点号及其对地阻抗形成的矩阵:X=');
Y=zeros(n);U=zeros(1,n);cta=zeros(1,n);V=zeros(1,n);O=zeros(1,n);S1=zeros(nl);
S=zeros(1,n);C=zeros(1,n);D=zeros(1,n);
for i=1:n
if X(i,2)~=0
p=X(i,1);
Y(p,p)=X(i,2);
end
end
for i=1:nl
if B1(i,6)==0
p=B1(i,1);q=B1(i,2);
else
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
end %求导纳矩阵
disp(Y);
G=real(Y);B=imag(Y);
for i=1:n
cta(i)=angle(B2(i,3));
U(i)=abs(B2(i,3));
end
for i=1:n
S(i)=B2(i,1)-B2(i,2);
end
P=real(S);Q=imag(S);
ICT1=0;IT2=1;ci=0;
while IT2~=0
IT2=0;t1=1;t2=1;ci=ci+1;
if ci>500
disp('迭代超过500次，跳出');
break;
end
for i=1:n
if i~=isb
C(i)=0;
D(i)=0;
for j1=1:n
C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
end
DP(t1)=P(i)-C(i);
t1=t1+1;
if B2(i,6)==2
DQ(t2)=Q(i)-D(i);
t2=t2+1;
end
end
end
t1=t1-1;t2=t2-1;
DPQ=; %求DP,DQ
for i=1:t1+t2
if abs(DPQ(i))>pr
IT2=IT2+1;
end
end
H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2);
a1=0;
for i=1:n
if i~=isb
a1=a1+1;
b1=0;
for j1=1:n
if j1~=isb
b1=b1+1;
if j1~=i
H(a1,b1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
elseif j1==i
H(a1,b1)=U(i)^2*B(i,j1)+D(i);
end
end
end
end
end
a1=0;
for i=1:n
if i~=isb
a1=a1+1;
b1=0;
for j1=1:n
if j1~=isb&&B2(j1,6)==2
b1=b1+1;
if j1~=i
N(a1,b1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
elseif j1==i
N(a1,b1)=0-U(i)^2*G(i,j1)-C(i);
end
end
end
end
end
a1=0;
fori=1:n
if i~=isb&&B2(i,6)==2
a1=a1+1;
b1=0;
for j1=1:n
if j1~=isb
b1=b1+1;
if j1~=i
K(a1,b1)= U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
elseif j1==i
K(a1,b1)=U(i)^2*G(i,j1)-C(i);
end
end
end
end
end
a1=0;
for i=1:n
if i~=isb&&B2(i,6)==2
a1=a1+1;
b1=0;
for j1=1:n
if j1~=isb&&B2(j1,6)==2
b1=b1+1;
if j1~=i
L(a1,b1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
elseif j1==i
L(a1,b1)=U(i)^2*B(i,j1)-D(i);
end
end
end
end
end
J=;%求雅可比矩阵
modify=-J\DPQ;
t4=0;
for i=1:n
if B2(i,6)~=1
t4=t4+1;
cta(1,i)=cta(1,i)+modify(t4,1);
end
end
for i=1:n
if B2(i,6)==2
t4=t4+1;
U(1,i)=U(1,i)+U(1,i)*modify(t4,1);
end
end
ICT1=ICT1+1;
end %修正原值
for i=1:n
UU(i)=U(i)*cos(cta(i))+1i*U(i)*sin(cta(i));
end
disp('--------------------------------------------------------------------------------');
disp('各节点电压U为（节点从小到大排列）:');
disp(UU);
disp('--------------------------------------------------------------------------------');
disp('各节点电压相角为（节点从小到大排列）:');
disp(180*angle(UU)/pi);
disp('--------------------------------------------------------------------------------');
disp('按公式计算全部线路功率，结果如下:');
for i=1:nl
if B1(i,6)==0
p=B1(i,1);q=B1(i,2);
else
p=B1(i,2);q=B1(i,1);
end
Si(p,q)=UU(p)*(conj(UU(p))*conj(B1(i,4)./2)+(conj(UU(p)*B1(i,5))-conj(UU(q)))*conj(1./(B1(i,3)*B1(i,5))));%各条支路首端功率Si
f=;
disp(f);
end
for i=1:nl
if B1(i,6)==0
p=B1(i,1);q=B1(i,2);
else p=B1(i,2);q=B1(i,1);
end
Sj(q,p)=UU(q)*(conj(UU(q))*conj(B1(i,4)./2)+(conj(UU(q)./B1(i,5))-conj(UU(p)))*conj(1./(B1(i,3)*B1(i,5))));%各条支路末端功率Sj
f=;
disp(f);
end
disp('--------------------------------------------------------------------------------');
disp('各条支路的功率损耗DS为(顺序同您输入B1时一样):');
for i=1:nl
if B1(i,6)==0
p=B1(i,1);q=B1(i,2);
else p=B1(i,2);q=B1(i,1);
end
DS(i)=Si(p,q)+Sj(q,p);%各条支路功率损耗DS
disp(DS(i));
end
Sp=0;
for i=1:n
Sp=Sp+UU(isb)*conj(Y(isb,i))*conj(UU(i));
end
disp('平衡节点的功率:');
disp(Sp);
页:
[1]