%本程序的功能是用牛頓-拉夫遜法進行潮流計算
%潮流結果：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）））;?%由