pfc.m

潮流计算的程序

% 潮流计算程序(PFC)
clear
clc
tic
format long
choice=0;
clc
    fprintf('欢迎使用潮流计算程序！\n');
    fprintf('请选择算例：\n\n');
    fprintf('    1.IEEE14节点算例\n');
    fprintf('    2.IEEE30节点算例\n');
    fprintf('    3.IEEE57节点算例\n');
    fprintf('    4.IEEE118节点算例\n');
    fprintf('    0.退出\n');
    choice = input('\nEnter your choice : ');
    clc
switch choice
    case 1
        load('D:\Program Files\MATLAB\R2008a\work\IEEE14data.mat')
        load('D:\Program Files\MATLAB\R2008a\work\state14.mat')
    case 2
        load('C:\MATLAB6p5p1\work\IEEE30data.mat')
        load('C:\MATLAB6p5p1\work\state30.mat')
    case 3
        load('C:\MATLAB6p5p1\work\IEEE57data.mat')
        load('C:\MATLAB6p5p1\work\state57.mat')
    case 4
        load('C:\MATLAB6p5p1\work\IEEE118data.mat')
        load('C:\MATLAB6p5p1\work\state118.mat')
    case 0
        clc
        break
    otherwise
        clc
        fprintf('Sorry, Invalid Choice !\n');
        break
end
% 形成节点导纳矩阵Y，w=sqrt(-1)
Y=zeros(node_number,node_number);
yground=zeros(node_number,1);
w=sqrt(-1);
% 形成互导纳
for m=1:branch_number
    i=branch(m,2);
    j=branch(m,3);
    type=branch(m,7);
    if type==1 % 一般支路
        g(i,j)=branch(m,4)/(branch(m,4)^2+branch(m,5)^2);
        g(j,i)=g(i,j);
        b(i,j)=-branch(m,5)/(branch(m,4)^2+branch(m,5)^2);
        b(j,i)=b(i,j);
        Y(i,j)=-g(i,j)-w*b(i,j);
        Y(j,i)=Y(i,j);
        yground(i)=yground(i)+w*branch(m,6);
        yground(j)=yground(j)+w*branch(m,6);
        yc(i,j)=branch(m,6);
        yc(j,i)=yc(i,j);
    end
    if type==0 % 变压器支路
        g(i,j)=0;
        g(j,i)=g(i,j);
        b(i,j)=-1/branch(m,5);
        b(j,i)=b(i,j);
        Y(i,j)=-w*b(i,j)/branch(m,8);
        Y(j,i)=Y(i,j);
        yground(i)=yground(i)+w*b(i,j)*(1-branch(m,8))/branch(m,8)^2;
        yground(j)=yground(j)+w*b(i,j)*(branch(m,8)-1)/branch(m,8);
        K(i,j)=branch(m,8);
        K(j,i)=K(i,j);
    end
end
% 形成自导纳（不含对地支路）
for m=1:node_number
    if m==1
        for n=2:node_number
            Y(m,m)=Y(m,m)-Y(m,n);
        end
    elseif m==node_number
        for n=1:node_number-1
            Y(m,m)=Y(m,m)-Y(m,n);
        end
    else
        for n=1:m-1
            Y(m,m)=Y(m,m)-Y(m,n);
        end
        for n=m+1:node_number
            Y(m,m)=Y(m,m)-Y(m,n);
        end
    end
end
% 形成自导纳（含对地支路）
for m=1:node_number
    Y(m,m)=Y(m,m)+yground(m);
end
G=real(Y);
B=imag(Y);
P=zeros(node_number,node_number);
Q=zeros(node_number,node_number);
for i=1:node_number
    state(2*i)=state(2*i)*pi/180;
end
% 节点注入功率
for i=1:node_number
    for j=1:node_number
        P(i,i)=P(i,i)+state(2*i-1)*state(2*j-1)*(G(i,j)*cos(state(2*i)-state(2*j))+B(i,j)*sin(state(2*i)-state(2*j)));
        Q(i,i)=Q(i,i)+state(2*i-1)*state(2*j-1)*(G(i,j)*sin(state(2*i)-state(2*j))-B(i,j)*cos(state(2*i)-state(2*j)));
    end
end
% 计算支路潮流
for m=1:branch_number
    i=branch(m,2);
    j=branch(m,3);
    t=branch(m,7);
    if t==1 % 一般支路
        P(i,j)=state(2*i-1)^2*g(i,j)-state(2*i-1)*state(2*j-1)*(g(i,j)*cos(state(2*i)-state(2*j))+b(i,j)*sin(state(2*i)-state(2*j)));
        P(j,i)=state(2*j-1)^2*g(i,j)+state(2*i-1)*state(2*j-1)*(-g(i,j)*cos(state(2*i)-state(2*j))+b(i,j)*sin(state(2*i)-state(2*j)));
        Q(i,j)=-state(2*i-1)^2*(b(i,j)+yc(i,j))-state(2*i-1)*state(2*j-1)*(g(i,j)*sin(state(2*i)-state(2*j))-b(i,j)*cos(state(2*i)-state(2*j)));
        Q(j,i)=-state(2*j-1)^2*(b(i,j)+yc(i,j))+state(2*i-1)*state(2*j-1)*(g(i,j)*sin(state(2*i)-state(2*j))+b(i,j)*cos(state(2*i)-state(2*j)));
    else % 变压器支路
        P(i,j)=-state(2*i-1)*state(2*j-1)*b(i,j)*sin(state(2*i)-state(2*j))/K(i,j);
        P(j,i)=state(2*i-1)*state(2*j-1)*b(i,j)*sin(state(2*i)-state(2*j))/K(i,j);
        Q(i,j)=-state(2*i-1)^2*b(i,j)/(K(i,j)^2)+state(2*i-1)*state(2*j-1)*b(i,j)*cos(state(2*i)-state(2*j))/K(i,j);
        Q(j,i)=-state(2*j-1)^2*b(i,j)+state(2*i-1)*state(2*j-1)*b(i,j)*cos(state(2*i)-state(2*j))/K(i,j);
    end
end
toc;
format short
t=toc
P;
Q;
PP=sparse(P)
QQ=sparse(Q)