tisco.m

电力系统中牛拉直角坐标方法的潮流程序

function tisco
% 这是一个电力系统潮流计算的程序
n=input('\n 请输入节点数:n=') ;
m=input('\n 请输入支路数: m =') ;
ph=input('\n 请输入平衡母线的节点号:ph=') ;
B1=input('\n 请输入支路信息: B1=') ;
% 它以矩阵形式存贮支路的情况, 每行存贮一条支路
% 第一列存贮支路的一个端点
% 第二列存贮支路的另一个端点
% 第三列存贮支路的阻抗
% 第四列存贮支路的对地导纳
% 第五列存贮变压器的变比
% 第六列存贮支路的序号
B2=input('\n 请输入节点信息: B2=') ;
% 第一列为电源侧的功率
% 第二列为负荷侧的功率
% 第三列为该点的电压值
% 第四列为该点的类型: 1 为PQ 节点,2 为PV 节点,3 为平衡节点
A =input('\n 请输入节点号及对地阻抗: A=') ;
ip=input('\n 请输入修正值: ip=') ;
% ip 为修正值
Y=zeros(n) ;
e=zeros(1,n) ;
f=zeros(1,n) ;
no=2*ph-1;
for i=1:n
if A(i,2)~=0
   p=A(i,1);
   Y(p,p)=1./A(i,2);
end
end
for i=1:m
p=B1(i,1) ;
q=B1(i,2) ;
Y(p,p)=Y(p,p)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
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,4)./2;
end
G=real(Y) ;
B=imag(Y) ;
for i=1:n
e(i)=real(B2(i,3)) ;
f(i)=imag(B2(i,3)) ;
S(i)=B2(i,1)-B2( i,2) ;
V(i)=B2(i,3);
end
P=real(S) ;
Q=imag(S) ;
disp('===========================================导纳矩阵Y为=============================================');
disp(Y);
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no) ;
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no) ;
[De,Df]=hxf(J,DF,ph,n,no);
t=0;
while max(abs(De))>ip&max(abs(Df))>ip
    t=t+1;
    e=e+De;
    f=f+Df;
    if t>15
        disp('不收敛');
        break
    end
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no) ;
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no) ;
[De,Df]=hxf(J,DF,ph,n,no);
end
v=e'+f'*j;
for i=1:n
hh(i)=conj(Y(ph,i)*v(i));
end
S(ph)=sum(hh)*v(ph) ;
B2(ph,1)=S(ph) ;
V=abs(v) ;
jd=angle(v)*180/pi;
result1=[A(:,1),real(v),imag(v),V,jd,real(S.'),imag(S.'),real(B2( :,1)),imag(B2( :,1)),real(B2(:,2)),imag(B2(:,2))];
for i=1:m
a(i)=conj((v(B1(i,1))/B1(i,5)-v(B1(i,2)))/B1(i,3));
b(i)=v(B1(i,1))*a(i)-j*B1(i,4)*v(B1(i,1))^2/2;
c(i)=-v(B1(i,2))*a(i)-j*B1(i,4)*v(B1(i,2))^2/2;
end
result2=[B1(:,6),B1(:,1),B1(:,2),real(b.') ,imag(b.') ,real(c.'),imag(c.'),real(b.'+c.'),imag(b.'+c.')];
disp('===========================================雅克比矩阵===============================================');
disp(J);
disp('===================================各节点电压为（按节点序号排列）=====================================');
disp(v);
disp('==================================各节点电压幅值为（按节点序号排列）===================================');
disp(result1(:,4));
disp('==================================各节点电压相角为（按节点序号排列）===================================');
disp(result1(:,5));
disp('========================================各支路线路功率损耗数据========================================');
disp(result2);

function [C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)
% 该子程序是用来求取D F
for i=1:n
if i~=ph
   C(i) =0;
   D(i) =0;
for j=1:n
C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j) ;
D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j) ;
end
P1=C(i)*e(i)+D(i)*f(i) ;
Q1=C(i)*f(i)-D(i)*e(i) ;
V1=e(i)^2+f(i)^2;
if B2(i,4) ==2
p=2*i-1;
DF(p)=P(i)-P1;
p=p+1;
DF(p)=V(i)^2-V1^2;
else
p=2*i-1;
DF(p)=P(i)-P1;
p=p+1;
DF(p)=Q(i)-Q1;
end
end
end
DF=DF';
if ph~=n
DF(no,:)=[ ] ;
DF(no,:)=[ ] ;
end

function J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)
% 该子程序是用来求取jacci矩阵
for i=1:n
switch B2(i,4)
case 3
continue
case 1
for j=1:n
if j~=i&j~=ph
    X1=G(i,j)*f(i)-B(i,j)*e(i) ;
    X2=G(i,j)*e(i)+B(i,j)*f(i) ;
    X3=-X2;
    X4=X1;
    p=2*i-1;
    q=2*j-1;
    J(p,q) =X1;
    m=p+1;
    J(m,q)=X3;
    q=q+1;
    J(p,q)=X2;
    J(m,q)=X4;
else if j==i&j~=ph
        X1=D(i)+G(i,i)*f(i)-B(i,i)*e(i) ;
        X2=C(i)+G(i,i)*e(i)+B(i,i)*f(i) ;
        X3=C(i)-G(i,i)*e(i)-B(i,i)*f(i) ;
        X4=-D(i)+G(i,i)*f(i)-B(i,i)*e(i) ;
        p=2*i-1;
        q=2*j-1;
        J(p,q) =X1;
        m=p+1;
        J(m,q)=X3;
        q=q+1;
        J(p,q)=X2;
        J(m,q)=X4;
end
end
end
case 2
  for j=1:n
      if j~=i& j~=ph
         X1=G(i,j)*f(i)-B(i,j)*e(i) ;
         X2=G(i,j)*e(i)+B(i,j)*f(i) ;
         X3=0;
         X4=0;
         p=2*i-1;
         q=2*j-1;
         J(p,q)=X1;
         m=p+1;
         J(m,q)=X3;
         q=q+1;
         J(p,q)=X2;
         J(m,q)=X4;
else if j==i& j~=ph
        X1=D(i)+G(i,i)*f(i)-B(i,i)*e(i) ;
        X2=C(i)+G(i,i)*e(i)+B(i,i)*f(i) ;
        X3=0;
        X4=0;
        p=2*i-1;
        q=2*j-1;
        J(p,q)=X1;
        m=p+1;
        J(m,q)=X3;
        q=q+1;
        J(p,q)=X2;
        J(m,q)=X4;
end
end
end
end
end
if ph~=n
J( no,:) =[ ] ;
J( no,:) =[ ] ;
J( :,no) =[ ] ;
J( :,no) =[ ] ;
end

function [De,Df]=hxf(J,DF,ph,n,no)
% 该子函数是为求取De Df
DX =J\DF;
DX1=DX ;
x1=length(DX1);
if ph~=n
  DX(no)=0;
  DX(no+1)=0;
for i=(no+2)(x1+2)
    DX(i)=DX1(i-2) ;
end
else
    DX =[ DX1;0;0] ;
end
k=0;
[x,y]=size(DX);
for i=1:2:x
k=k+1;
Df(k)=DX(i) ;
De(k)=DX(i+1) ;
end