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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%設(shè)各節(jié)點(diǎn)電壓的初值
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%計(jì)算PQ節(jié)點(diǎn)和PV節(jié)點(diǎn)的功率不平衡量
????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é)點(diǎn)電壓求得的PQ節(jié)點(diǎn)注入有功功率
????????????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é)點(diǎn)電壓求得的PQ節(jié)點(diǎn)注入無(wú)功功率
????????end
????????Unbalance（2*m-1）=p（indpq（m））-sum（pt）;?%計(jì)算PQ節(jié)點(diǎn)有功功率不平衡量
????????Unbalance（2*m??）=q（indpq（m））-sum（qt）;?%計(jì)算PQ節(jié)點(diǎn)無(wú)功功率不平衡量
end?%[Unbalance]是節(jié)點(diǎn)不平衡量矩陣
????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）））;?%由