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節(jié)點(diǎn)的電壓給定值
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;%最大迭代次數(shù)
varepsilon=1e-10;%收斂精度
%%?牛頓—拉夫遜法計(jì)算各節(jié)點(diǎn)電壓
G=real（Y）;B=imag（Y）;%求G、B
e=real（U）;f=imag（U）;%求e、f
index_PQ=find（~isnan（P）&~isnan（Q））;%PQ節(jié)點(diǎn)編號
index_PV=find（~isnan（P）&~isnan（Us））;%PV節(jié)點(diǎn)編號
index_P=find（~isnan（P））;%PV節(jié)點(diǎn)及PV節(jié)點(diǎn)的編號
for?k=1:N
????[DeltaPDeltaQDeltaW]=CAL1（PQUsGBefindex_PQindex_PVindex_P）;
????Delta=max（abs（[DeltaP;DeltaQ]））;%最大絕對誤差
????if?Delta<=varepsilonbreakend%若達(dá)到收斂精度條件，則退出循環(huán)
????Jacobi=CAL2（GBefindex_PQindex_PVindex_P）;%求雅可比矩陣
????DeltaV=-Jacobi\DeltaW;%求△e、△f
????%修正各節(jié)點(diǎn)電壓
????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節(jié)點(diǎ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）;
?