matlab潮流计算

Line=[1  1  2   0.04    0.25     0.25   

     2  1  3    0.1     0.35    0  

     3  2  3   0.08     0.30    0.25  ];

 %变压器数据  transform=[1-支路编号  2-支路首节点编号  3-支路末节点编号  4-支路电阻(p.u.)  5-支路电抗(p.u.)  6-变压器变比(p.u.)]

transform=[ 1  2   4  0    0.015    1.05

2  3   5  0    0.03     1.05];

%数据预处理

 Nbus=5;                     % 节点数 

 nline=size(Line,1);                   % 线路个数

 ntrans=size(transform,1);             % 变压器个数

 slack=5;       % 平衡节点号

 Npq=4;    % PQ节点的个数 

Ss=0;

%计算节点导纳矩阵

Y=zeros(Nbus);

if nline>=1%判断是否存在线路

    for k=1:nline               %以下处理线路

        t1=Line(k,2); t2=Line(k,3); b2=Line(k,6);%分别取出线路的首端节点编号t1、末端节点编号t2和对地电纳b2

        Yl=1/(Line(k,4)+j*Line(k,5));%计算线路的支路电导Yl

        Y(t1,t1)=Y(t1,t1)+Yl+j*b2;%修正第k条线路首端节点的自导纳

        Y(t1,t2)=Y(t1,t2)-Yl;%修正第k条线路首端节点与末端节点之间的互导纳

        Y(t2,t1)=Y(t2,t1)-Yl;%修正第k条线路末端节点与首端节点之间的互导纳

        Y(t2,t2)=Y(t2,t2)+Yl+j*b2;%修正第k条线路末端节点的自导纳

    end

end

if ntrans>=1%判断是否存在变压器

    for k=1:ntrans   %以下处理变压器

        t1=transform(k,2);t2=transform(k,3);t3=transform(k,6);%分别取出变压器的首端节点编号t1、末端节点编号t2和变比t3

        Yt=1/(transform(k,4)+j*transform(k,5));

        Yt1=Yt/t3;

        Yt2=Yt*(1-t3)/(t3*t3);

        Yt3=Yt*(t3-1)/t3;

        Y(t1,t1)=Y(t1,t1)+Yt1+Yt2;

        Y(t1,t2)=Y(t1,t2)-Yt1;

        Y(t2,t1)=Y(t2,t1)-Yt1;

        Y(t2,t2)=Y(t2,t2)+Yt1+Yt3;

    end

end

G=real(Y);B=imag(Y); %区分节点导纳矩阵的实部和虚部

G

B

%赋初值

delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;

u(1)=1;u(2)=1;u(3)=1;u(4)=1;

p(1)=-0.30;q(1)=-0.18;p(2)=-0.55;q(2)=-0.13;p(3)=0;q(3)=0;

p(4)=0.5;q(4)=1.10;p(5)=0.8;q(5)=0.50;

k=0;precision=1;Npq=4; %Npq分别是网络中的PQ节点数

%[Unbalance]=-[Jacobi][Correction]

while precision>0.00001 %设定误差上限，判断是否继续迭代

    u(5)=1.06; delt(5)=0; %设定平衡节点电压相角与幅值

    k; u; delt;

    for m=1:Npq

        for n=1:Nbus

            pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n))); %由节点电压求得的PQ节点注入有功功率

            qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); %由节点电压求得的PQ节点注入无功功率

        end

        Unbalance(2*m-1)=p(m)-sum(pt); %计算PQ节点有功功率不平衡量

        Unbalance(2*m  )=q(m)-sum(qt); %计算PQ节点无功功率不平衡量

end %[Unbalance]是节点不平衡量矩阵

    for m=1:Npq

        for n=1:Nbus

            h0(n)= u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

            n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

            j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

            l0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

        end

        H(m,m)=sum(h0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));

        N(m,m)=sum(n0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)))-2*u(m)^2*G(m,m);

        J(m,m)=sum(j0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));

        L(m,m)=sum(l0)+u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)))+2*u(m)^2*B(m,m);

        Jacobi(2*m-1,2*m-1)=H(m,m);

        Jacobi(2*m-1,2*m  )=N(m,m);

        Jacobi(2*m  ,2*m-1)=J(m,m);

        Jacobi(2*m  ,2*m  )=L(m,m);

    end %计算m=n情况下的Jacobi矩阵中的子矩阵元素

    for m=1:Npq

        for n=1:Npq

            if m==n

            else

                H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

                J(m,n)= u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

                N(m,n)=-J(m,n);

                L(m,n)= H(m,n);

                Jacobi(2*m-1,2*n-1)=H(m,n);

                Jacobi(2*m-1,2*n  )=N(m,n);

                Jacobi(2*m  ,2*n-1)=J(m,n);

                Jacobi(2*m  ,2*n  )=L(m,n);

            end

        end

    end %计算m≠n情况下的Jacobi矩阵中的子矩阵元素

    Correction=-Jacobi\\(Unbalance');  %计算电压相角和幅值的修正量

    precision = max(abs(Correction)); %取误差最大值

    for m=1:Npq

        delt(m)=delt(m)+Correction(2*m-1); %修正PQ节点电压相角

            u(m)=    u(m)+Correction(2*m  ); %修正PQ节点电压幅值

    end

    k=k+1; %迭代轮数+1

end

k, u, delt, Jacobi, precision

for m=1:Nbus

    U(m)=u(m)*(cos(delt(m))+j*sin(delt(m))); %采用直角坐标系表示电压

    I(m)=Y(Nbus,m)*U(m);                       %计算注入平衡节点的电流

end

Sslack=U(Nbus)*sum(conj(I)) %计算注入平衡节点的功率

for m=1:Nbus

    for n=1:Nbus

        S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n)); %计算线路功率

        Ploss(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

    end

end

S %显示线路功率矩阵，S(m,n)表示假设功率由节点m流向节点n

  %若数值为+则说明实际功率流向与假设方向相同，若数值为-则说明实际功率流向与假设方向相反