极坐标的牛顿法潮流计算通用程序

clear
clc
%极坐标形式的牛顿法潮流计算程序
%B1矩阵为支路参数矩阵，其中1、2：分别为支路首端、末端节点编号；3：支路阻抗；4：接地导纳；5：支路变比；6：支路首端处于K侧为11侧为0；7：单侧对地导纳
%B2矩阵为节点参数矩阵，其中1：该节点的发电机功率；2：该节点负荷功率；3：节点电压初始幅值；4：节点电压初始相位；5：PV节点的电压给定值；6：节点所接的无功补偿设备容量；
%7：节点分类编号：1为平衡节点，2为PQ节点，3为PV节点
n=input（‘请输入节点数:n=‘）;
nl=input（‘请输入支路数:nl=‘）;
isb=input（‘请输入平衡节点编号:isb=‘）;
pr=input（‘请输入误差精度:pr=‘）;
B1=input（‘请输入支路参数矩阵:B1=‘）;
B2=input（‘请输入节点参数矩阵:B2=‘）;
%创建矩阵
Y=zeros（n）;
V=zeros（n1）;
for i=1:nl
    if B1（i6）==0 %支路不含变压器
        p=B1（i1）;q=B1（i2）;
        Y（pq）=Y（pq）-1./B1（i3）;%导纳矩阵的非对角元
        Y（qp）=Y（pq）;
        Y（pp）=Y（pp）+1./B1（i3）+0.5*B1（i4）+B1（i7）;%导纳矩阵的对角元
        Y（qq）=Y（qq）+1./B1（i3）+0.5*B1（i4）+B1（i7）;
    else 
        p=B1（i2）;q=B1（i1）;%支路含有变压器
        Y（pq）=Y（pq）-1./（B1（i3）*B1（i5））;
        Y（qp）=Y（pq）;
        Y（pp）=Y（pp）+1./B1（i3）+B1（i7）; %pi型等效电路
        Y（qq）=Y（qq）+1./（（B1（i5）^2*B1（i3）））+B1（i7）;
    end
end
%分解导纳的实部虚部
G=real（Y）;B=imag（Y）;
for i=1:n
    theta（i1）=B2（i4）;
    V（i1）=B2（i3）;
end
%求出两节点的相位差
Times=0;
%求出PQ节点的个数
m=0;
for i=1:n
    if B2（i7）==2
        m=m+1;
    end
end
M=ones（n-1+m1）;
while max（abs（M））> pr
%形成潮流方程
P=zeros（n1）;
Q=zeros（n1）;
for i=1:n
    for j=1:n
        P（i1）=V（i）*sum（（G（i:）.*cos（theta（i）-theta（:））‘+B（i:）.*sin（theta（i）-theta（:））‘）*V（:））;
        Q（i1）=V（i）*sum（（G（i:）.*sin（theta（i）-theta（:））‘-B（i:）.*cos（theta（i）-theta（:））‘）*V（:））;
    end
end
%求出功率的不平衡量
dP=zero