牛拉法（极坐标形式）

%本程序的功能是用牛顿-拉夫逊法进行潮流计算
%潮流结果：U——节点电压
%         I——平衡节点处电流
%         Sslack——平衡节点处注入功率
%         S——显示线路功率矩阵，S（mn）表示假设功率由节点m流向节点n若数值为+则说明实际功率流向与假设方向相同，若数值为-则说明实际功率流向与假设方向相反
%         Ploss——线路网损
%潮流开始
%输入线路数据  
%Line=[1-支路编号  2-线路首端节点号 3-线路末端节点号  4-支路电阻 5-支路电抗  6-支路电纳（注意：此处取的是B/2）]
%Line=input（‘请输入由支路参数形成的矩阵:Line=‘）;
load Line.txt  %读入线路数据   
%输入变压器数据  
%transform=[1-支路编号  2-支路首节点编号  3-支路末节点编号  4-支路电阻（p.u.）  5-支路电抗（p.u.）  6-变压器变比（p.u.）]
%transform=input（‘请输入变压器支路数据‘）;
load transform.txt  %读入变压器数据   
%输入节点数据
%node=[1-节点号 2-节点类型（PQ节点：1PV节点：2平衡节点：3） 3-电压幅值初值 4-电压相角初值 5-节点负荷有功 6-节点负荷无功]
%node=input（‘请输入节点数据‘）;
load node.txt  %读入节点数据   

%数据预处理
Nbus=input（‘请输入节点数：n=‘）;
slack=input（‘请输入平衡节点号：slack=‘）;
Npq=input（‘请输入PQ节点个数：Npq=‘）;
Npv=input（‘请输入PV节点个数：Npv=‘）;
pre=input（‘请输入误差精度：pre=‘）;
nline=size（Line1）;                   % 线路个数
ntrans=size（transform1）;             % 变压器个数

Ss=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算节点导纳矩阵
Y=zeros（Nbus）;
if nline>=1%判断是否存在线路
    for k=1:nline               %以下处理线路
        t1=Line（k2）; t2=Line（k3）; b2=Line（k6）;%分别取出线路的首端节点编号t1、末端节点编号t2和对地电纳b2
        Yl=1/（Line（k4）+j*Line（k5））;%计算线路的支路电导Yl
        Y（t1t1）=Y（t1t1）+Yl+j*b2;%修正第k条线路首端节点的自导纳
        Y（t1t2）=Y（t1t2）-Yl;%修正第k条线路首端节点与末端节点之间的互导纳
        Y（t2t1）=Y（t2t1）-Yl;%修正第k条线路末端节点与首端节点之间的互导纳
        Y（t2t2）=Y（t2t2）+Yl+j*b2;%修正第k条线路末端节点的自导纳
    end
end
if ntrans>=1%判断是否存在变压器
    for k=1:ntrans   %以下处理变压器
        t1=transform（k2）;t2=transform（k3）;t3=transform（k6）;%分别取出变压器的首端节点编号t1、末端节点编号t2和变比t3
        Yt=1/（transform（k4）+j*transform（k5））;
        Yt1=Yt/t3;
        Yt2=Yt*（1-t3）/（t3*t3）;
        Yt3=Yt*（t3-1）/t3;
        Y（t1t1）=Y（t1t1）+Yt1+Yt2;
        Y（t1t2）=Y（t1t2）-Yt1;
        Y（t2t1）=Y（t2t1）-Yt1;
        Y（t2t2）=Y（t2t2）+Yt1+Yt3;
    end
end
G=real（Y）;B=imag（Y）; %区分节点导纳矩阵的实部和虚部

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%设各节点电压的初值
u=node（:3）;
delt=node（:4）;
p=node（:5）;
q=node（:6）;
k=0;precision=1;%设定迭代次数初值和误差初值
%[Unbalance]=-[Jacobi][Correction]

while precision>pre %设定误差上限，判断是否继续迭代
    k; u; delt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算PQ节点和PV节点的功率不平衡量
    indpq=find（node（:2）==1）;
    for m=1:Npq
        for n=1:Nbus
            pt（n）=u（indpq（m））*u（n）*（G（indpq（m）n）*cos（delt（indpq（m））-delt（n））+B（indpq（m）n）*sin（delt（indpq（m））-delt（n）））; %由节点电压求得的PQ节点注入有功功率
            qt（n）=u（indpq（m））*u（n）*（G（indpq（m）n）*sin（delt（indpq（m））-delt（n））-B（indpq（m）n）*cos（delt（indpq（m））-delt（n）））; %由节点电压求得的PQ节点注入无功功率
        end
        Unbalance（2*m-1）=p（indpq（m））-sum（pt）; %计算PQ节点有功功率不平衡量
        Unbalance（2*m  ）=q（indpq（m））-sum（qt）; %计算PQ节点无功功率不平衡量
end %[Unbalance]是节点不平衡量矩阵
    indpv=find（node（:2）==2）;
    for m=1:Npv
        for n=1:Nbus
         pt（n）=u（indpv（m））*u（n）*（G（indpv（m）n）*cos（delt（indpv（m））-delt（n））+B（indpv（m）n）*sin（delt（indpv（m））-delt（n）））; %由