潮流计算程序牛顿拉夫逊法

function ChenDandan
%% 前处理
clearclcclose all
format compact
P=[0 -0.5 0.2 NaN]‘;
Q=[0 -0.3 NaN NaN]‘;
Us=[NaN NaN 1.05 NaN]‘;%PV节点的电压给定值
U=[1.05 1.05 1.05 1.05]‘;%电压初值
y0=[0 -0.68182i 0.013023i 0.017303i
    0.60434i 0.05i 0 0
    0.013023i 0 0 0.025864i
    0.017303i 0 0.025864i 0];
Y=[2.9049-11.503i 5.3182i -1.66+3.1619i -1.245+2.3714i
    5.3182i -4.6638i 0 0
    -1.66+3.1619i 0 2.49-4.704i -0.83+1.5809i
    -1.245+2.3714i 0 -0.83+1.5809i 2.075-3.9092i];
N=100;%最大迭代次数
varepsilon=1e-10;%收敛精度
%% 牛顿—拉夫逊法计算各节点电压
G=real（Y）;B=imag（Y）;%求G、B
e=real（U）;f=imag（U）;%求e、f
index_PQ=find（~isnan（P）&~isnan（Q））;%PQ节点编号
index_PV=find（~isnan（P）&~isnan（Us））;%PV节点编号
index_P=find（~isnan（P））;%PV节点及PV节点的编号
for k=1:N
    [DeltaPDeltaQDeltaW]=CAL1（PQUsGBefindex_PQindex_PVindex_P）;
    Delta=max（abs（[DeltaP;DeltaQ]））;%最大绝对误差
    if Delta<=varepsilonbreakend%若达到收敛精度条件，则退出循环
    Jacobi=CAL2（GBefindex_PQindex_PVindex_P）;%求雅可比矩阵
    DeltaV=-Jacobi\DeltaW;%求△e、△f
    %修正各节点电压
    e（index_P）=DeltaV（2:2:end）+e（index_P）;
    f（index_P）=DeltaV（1:2:end）+f（index_P）;
end
U=e+f*1j;%电压
%% 后处理
fprintf（‘电压大小：\n‘）
disp（abs（U））

fprintf（‘电压相角（弧度）：\n‘）
disp（angle（U））

n=length（U）;

S_zr=zeros（n1）;
for k=1:n
    S_zr（k）=U（k）*sum（conj（Y（k:））.*conj（U.‘））;
end
fprintf（‘\n节点注入功率为：\n‘）
disp（S_zr）

S=zeros（n）;
for iRow=1:size（S1）
    for iColumn=1:size（S2）
        S（iRowiColumn）=abs（U（iRow））^2*conj（y0（iRowiColumn））+...
            U（iRow）*（conj（U（iRow））-conj（U（iColumn）））*conj（-Y（iRowiColumn））;
    end
end
fprintf（‘\n线路上流动的功率为：\n‘）
disp（S）

Delta=zeros（n）;
for iRow=1:size（S1）
    for iColumn=1:size（S2）
        if iRow~=iColumn
            Delta（iRowiColumn）=S（iRowiColumn）+S（iColumniRow）;