chaoliu.asv

电力系统潮流计算程序

function y
clear;
clc;
%nl=5;%网络支路数
%n=5;%网络节点数
pr=0.00001;%迭代精度
V=[1.05+0i 1 1.05;1+0i 2 0;1+0i 2 0;1+0i 2 0;1.05+0i 3 1.05];
%——V中数据依次为——各节点电压初始值、节点分类号('1'为平衡节点，'2'为PQ节点，'3'为PV节点)、节点给定功率
%节点电压、节点分类号、节点所接的无功补偿设备容量
b=[2 1 0.03i 0 1.05;3 2 0.08+0.3i 0.5i 1;4 2 0.1+0.35i 0 1;4 3 0.04+0.25i 0.5i 1;3 5 0.015i 0 1.05]
%——b中数据依次为——支路始端号、支路末端号、支路阻抗、线路电容、支路变比
%b=[1 5 0.15i 0 1;5 6 0.0025+0.025i 0.0438i 1;6 7 0.001+0.01i 0.0175i 1;6 2 0.15i 0 1;7 8 0.011+0.11i 0.148i 1;8 9 0.001+0.01i 0.0175i 1;9 10 0.0025+0.025i 0.0438i 1;10 3 0.15i 0 1;10 4 0.15i 0 1]
%支路始端号、支路末端号、支路阻抗、线路电容、支路变比
S=[0+0i;-3.7-1.3i;-2-1i;-1.6-0.8i;5];
%各节点的注入功率
%S=[0+0i;-3.7-1.3i;-7-5i;-1.6-0.8i;5];%各节点的注入功率
nl=size(b,1);
n=size(V,1);
disp('支路数：')
disp(nl)
disp('节点数：')
disp(n)
w1=zeros(2*n-2,1);  
P=real(S);Q=imag(S);e=zeros(1,n);f=zeros(1,n);E=zeros(1,n);
Y=zeros(n);
for i=1:nl%导纳矩阵生成
    p=b(i,1);q=b(i,2);
    Y(p,q)=Y(p,q)-1./(b(i,3)*b(i,5));
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./b(i,3)+b(i,4)./2;
    Y(p,p)=Y(p,p)+1./(b(i,3)*b(i,5)^2)+b(i,4)./2;
end
disp('导纳矩阵Y：')
disp(Y)
U=zeros(1,n);
G=real(Y);B=imag(Y);
for i=1:n
e(i)=real(V(i,1));f(i)=imag(V(i,1));
U(i)=V(i,3);B(i,i)=B(i,i)+V(i,3);
end
T=0;co=0;d=0;
while T==0
    A=0;co=co+1;
for i=2:n %生成雅可比矩阵和功率修正量
    for j=2:n
        x=0;x1=0;
        if V(i,2)==2
            for r=1:n
                x=x+(e(i)*(G(i,r)*e(r)-B(i,r)*f(r))+f(i)*(G(i,r)*f(r)+B(i,r)*e(r)));
                x1=x1+(f(i)*(G(i,r)*e(r)-B(i,r)*f(r))-e(i)*(G(i,r)*f(r)+B(i,r)*e(r)));
            end
            w(2*i-1)=P(i)-x;
            w(2*i)=Q(i)-x1;
        else if V(i,2)==3
                for r=1:n
                    x=x+(e(i)*(G(i,r)*e(r)-B(i,r)*f(r))+f(i)*(G(i,r)*f(r)+B(i,r)*e(r)));
                end
                w(2*i-1)=P(i)-x;
                w(2*i)=U(i)^2-(e(i)^2+f(i)^2);
            end
        end                        
        h=0;h1=0;
        if V(i,2)==2
            if i==j
                for r=1:n
                    if r==i
                        continue
                    end
                    h=h+(G(i,r)*f(r)+B(i,r)*e(r));
                    h1=h1+(G(i,r)*e(r)-B(i,r)*f(r));
                end
                J(2*i-1,2*j-1)=2*G(i,i)*f(i)+h;
                J(2*i-1,2*j)=2*G(i,i)*e(i)+h1;
                J(2*i,2*j-1)=-2*B(i,i)*f(i)+h1;
                J(2*i,2*j)=-2*B(i,i)*e(i)-h;
            else
                J(2*i-1,2*j-1)=-B(i,j)*e(i)+G(i,j)*f(i);
                J(2*i-1,2*j)=G(i,j)*e(i)+B(i,j)*f(i);
                J(2*i,2*j-1)=-G(i,j)*e(i)-B(i,j)*f(i);
                J(2*i,2*j)=-B(i,j)*e(i)+G(i,j)*f(i);
            end 
        else if V(i,2)==3
                if i==j
                    for r=1:n
                        if r==i
                           continue
                        end
                    h=h+(G(i,r)*f(r)+B(i,r)*e(r));
                    h1=h1+(G(i,r)*e(r)-B(i,r)*f(r));
                    end
                    J(2*i-1,2*j-1)=2*G(i,i)*f(i)+h;
                    J(2*i-1,2*j)=2*G(i,i)*e(i)+h1;
                    J(2*i,2*j-1)=2*f(i);
                    J(2*i,2*j)=2*e(i);
                else
                    J(2*i-1,2*j-1)=-B(i,j)*e(i)+G(i,j)*f(i);
                    J(2*i-1,2*j)=G(i,j)*e(i)+B(i,j)*f(i);
                    J(2*i,2*j-1)=0;
                    J(2*i,2*j)=0;
                end
            end
        end
    end
end
%disp(J)
%disp(w)
for i=3:2*n%高斯消去法求电压修正量
    for j=3:2*n
        J1(i-2,j-2)=J(i,j);
    end
end
for i=3:2*n
    w1(i-2)=w(i);
end
u=zeros(2*n-2,1);
N=2*n-2;
for k=1:N
    m=0;
    for i=k+1:N        
        m=J1(i,k)./J1(k,k);
        w1(i)=w1(i)-m*w1(k);         
            for j=k+1:N
              J1(i,j)=J1(i,j)-m*J1(k,j);
            end
    end
end
u(N)=w1(N)./J1(N,N);
for i=N-1:-1:1
    c=0;
    for k=i+1:N       
        c=c+J1(i,k)*u(k);
        u(i)=(w1(i)-c)./J1(i,i); 
    end
end
%disp(u)
        for i=1:2*n-2
            Jd=abs(u(i));
            if Jd>pr
                A=A+1;
            end
        end
        bm(co)=A;
        if A==0
            T=1;
        else
            for i=1:n-1
                f(i+1)=f(i+1)+u(2*i-1);
                e(i+1)=e(i+1)+u(2*i);
            end
            d=d+1;
        end
    end
disp('迭代次数=')
disp(d)
disp('每次不满足个数=')
disp(bm)
for i=1:n
    V1(i)=sqrt(e(i)^2+f(i)^2);
    O(i)=atan(f(i)./e(i))*180./pi;
end
disp('各节点电压大小=')
disp(V1)
disp('各节点电压相角=')
disp(O)
E=complex(e,f);
disp('节点电压=')
disp(E)
 for i=1:n%各节点功率
    o1=0;
    for j=1:n
        o1=o1+conj(Y(i,j))*conj(E(j));
    end
    S1(i)=E(i)*o1;
end
disp('各节点功率=')
disp(S1)
for i=1:nl%各支路首末端功率
    p=b(i,1);q=b(i,2);
    if b(i,5)==1
        S2(p,q)=E(p)*(conj(E(p))*conj(b(i,4)./2)+(conj(E(p))-conj(E(q)))*conj(1./b(i,3)));
        S2(q,p)=E(q)*(conj(E(q))*conj(b(i,4)./2)+(conj(E(q))-conj(E(p)))*conj(1./b(i,3)));
    else
        S2(q,p)=-E(q)*conj(((E(p)./b(i,5))-E(q))./b(i,3));
        S2(p,q)=E(p)*conj((E(p)./b(i,5)-E(q))*(1./(b(i,5)*b(i,3))));
    end
end
disp('各支路首端功率=')
for i=1:nl
    p=b(i,1);q=b(i,2);
    disp(S2(p,q))
end
disp('各支路末端功率=')
for i=1:nl
    p=b(i,1);q=b(i,2);
    disp(S2(q,p))
end
disp('各支路功率损耗=')
for i=1:nl%网络总损耗
    p=b(i,1);q=b(i,2);
    DS(p,q)=S2(p,q)+S2(q,p);
    disp(DS(p,q))
end
disp('网络总损耗=')
D=0;
for i=1:n
    D=D+S1(i);
end
disp(D)