%本程序的功能是用牛頓-拉夫遜法進行潮流計算
%潮流結(jié)果：U——節(jié)點電壓
%?????????I——平衡節(jié)點處電流
%?????????Sslack——平衡節(jié)點處注入功率
%?????????S——顯示線路功率矩陣，S（mn）表示假設(shè)功率由節(jié)點m流向節(jié)點n若數(shù)值為+則說明實際功率流向與假設(shè)方向相同，若數(shù)值為-則說明實際功率流向與假設(shè)方向相反
%?????????Ploss——線路網(wǎng)損
%潮流開始
%輸入線路數(shù)據(jù)??
%Line=[1-支路編號??2-線路首端節(jié)點號?3-線路末端節(jié)點號??4-支路電阻?5-支路電抗??6-支路電納（注意：此處取的是B/2）]
%Line=input（‘請輸入由支路參數(shù)形成的矩陣:Line=‘）;
load?Line.txt??%讀入線路數(shù)據(jù)???
%輸入變壓器數(shù)據(jù)??
%transform=[1-支路編號??2-支路首節(jié)點編號??3-支路末節(jié)點編號??4-支路電阻（p.u.）??5-支路電抗（p.u.）??6-變壓器變比（p.u.）]
%transform=input（‘請輸入變壓器支路數(shù)據(jù)‘）;
load?transform.txt??%讀入變壓器數(shù)據(jù)???
%輸入節(jié)點數(shù)據(jù)
%node=[1-節(jié)點號?2-節(jié)點類型（PQ節(jié)點：1PV節(jié)點：2平衡節(jié)點：3）?3-電壓幅值初值?4-電壓相角初值?5-節(jié)點負荷有功?6-節(jié)點負荷無功]
%node=input（‘請輸入節(jié)點數(shù)據(jù)‘）;
load?node.txt??%讀入節(jié)點數(shù)據(jù)???

%數(shù)據(jù)預(yù)處理
Nbus=input（‘請輸入節(jié)點數(shù)：n=‘）;
slack=input（‘請輸入平衡節(jié)點號：slack=‘）;
Npq=input（‘請輸入PQ節(jié)點個數(shù)：Npq=‘）;
Npv=input（‘請輸入PV節(jié)點個數(shù)：Npv=‘）;
pre=input（‘請輸入誤差精度：pre=‘）;
nline=size（Line1）;???????????????????%?線路個數(shù)
ntrans=size（transform1）;?????????????%?變壓器個數(shù)

Ss=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%計算節(jié)點導納矩陣
Y=zeros（Nbus）;
if?nline>=1%判斷是否存在線路
????for?k=1:nline???????????????%以下處理線路
????????t1=Line（k2）;?t2=Line（k3）;?b2=Line（k6）;%分別取出線路的首端節(jié)點編號t1、末端節(jié)點編號t2和對地電納b2
????????Yl=1/（Line（k4）+j*Line（k5））;%計算線路的支路電導Yl
????????Y（t1t1）=Y（t1t1）+Yl+j*b2;%修正第k條線路首端節(jié)點的自導納
????????Y（t1t2）=Y（t1t2）-Yl;%修正第k條線路首端節(jié)點與末端節(jié)點之間的互導納
????????Y（t2t1）=Y（t2t1）-Yl;%修正第k條線路末端節(jié)點與首端節(jié)點之間的互導納
????????Y（t2t2）=Y（t2t2）+Yl+j*b2;%修正第k條線路末端節(jié)點的自導納
????end
end
if?ntrans>=1%判斷是否存在變壓器
????for?k=1:ntrans???%以下處理變壓器
????????t1=transform（k2）;t2=transform（k3）;t3=transform（k6）;%分別取出變壓器的首端節(jié)點編號t1、末端節(jié)點編號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）;?%區(qū)分節(jié)點導納矩陣的實部和虛部

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%設(shè)各節(jié)點電壓的初值
u=node（:3）;
delt=node（:4）;
p=node（:5）;
q=node（:6）;
k=0;precision=1;%設(shè)定迭代次數(shù)初值和誤差初值
%[Unbalance]=-[Jacobi][Correction]

while?precision>pre?%設(shè)定誤差上限，判斷是否繼續(xù)迭代
????k;?u;?delt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%計算PQ節(jié)點和PV節(jié)點的功率不平衡量
????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）））;?%由節(jié)點電壓求得的PQ節(jié)點注入有功功率
????????????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）））;?%由節(jié)點電壓求得的PQ節(jié)點注入無功功率
????????end
????????Unbalance（2*m-1）=p（indpq（m））-sum（pt）;?%計算PQ節(jié)點有功功率不平衡量
????????Unbalance（2*m??）=q（indpq（m））-sum（qt）;?%計算PQ節(jié)點無功功率不平衡量
end?%[Unbalance]是節(jié)點不平衡量矩陣
????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）））;?%由