%本程序的功能是用牛頓——拉夫遜法進行潮流計算
%?B1矩陣：1、支路首端號；???2、末端號；???3、支路阻抗；???4、線路對地電納?（或變壓器導納）；
%		??5、支路的變比；???6、支路首端處于K側為1，1側為0；
%??????????7、線路/變壓器標識（0/1）變壓器參數當支路首端處于K側標識為1時歸算至末端側，0歸算至首端側
%?B2矩陣：1、該節點發電機功率；??????2、該節點負荷功率；
%?????????3、PQ節點電壓初始值，或平衡節點及PV節點電壓的給定值
%?????????4、節點所接無功補償并聯電容（感）的電納
%?????????5、節點分類標號：1為平衡節點（應為1號節點）；2為PQ節點；3為PV節點；
clear;
isb=1;??????%input（‘請輸入平衡母線節點號：isb=‘）;
pr=0.00001;????%input（‘請輸入誤差精度：pr=‘）;
%---------------------------------------------------
n=6;%input（‘請輸入節點數：n=‘）;
nl=3;%input（‘請輸入支路數：nl=‘）;
B1=[1?2?01i?0?]
B2=[0?????0????????1.20??????0???1;
???0???2.1+1i??????1.00??????0???2;
???0???1.8+0.40i???1.00??????0???2;
???0???1.6+0.8i????1.00??????0???2;
???0???3.7+1.3i????1.00??????0???2;
???0???-5+0i???????1.20??????0???3]?
%input（‘請輸入各節點參數形成的矩陣：?B2=‘）;
%X=[1?0;2?0;3?0;4?0;5?0;6?0]
%-------------------------------------------------------------
%n=4;%input（‘請輸入節點數：n=‘）;nl=4;%input（‘請輸入支路數：nl=‘）;
%B1=[1?2?4+16i?0?1?0?0;1?3?4+16i?0?1?0?0;2?3?2+8i?0?1?0?0;2?4?1.49+48.02i?0?11/110?0?1]?%input（‘請輸入由支路參數形成的矩陣：?B1=‘）;
%B2=[0?0?115?0?1;0?0?110?0?2;0?20+4i?110?0?2;0?10+6i?10?0?2]?%input（‘請輸入各節點參數形成的矩陣：?B2=‘）;
%-------------------------------------------------------------
Y=zeros（n）;e=zeros（1n）;f=zeros（1n）;V=zeros（1n）;sida=zeros（1n）;S1=zeros（nl）;
%?%?%-----------求導納矩陣------------------------
%for?i=1:n
???%?if?X（i2）~=0;
??????%??p=X（i1）;
???????%?Y（pp）=1/X（i2）;
????%end
%end
for?i=1:nl								%從1到n1（總支路數）
	if?B1（i7）==1????????%-----------如果是變壓器支路--------
????	if?B1（i6）==0						%左節點（首端）處于1側
???????		p=B1（i1）;q=B1（i2）;
????	else????????????????????????????????%左節點（首端）處于K側
???????		p=B1（i2）;q=B1（i1）;
????	end
????	Y（pq）=Y（pq）-1./（B1（i3）*B1（i5））;	%非對角元
????	Y（qp）=Y（pq）;??????????????????????%非對角元
???		Y（qq）=Y（qq）+1./（B1（i3）*B1（i5）^2）;	%對角元K側
????	Y（pp）=Y（pp）+1./B1（i3）+B1（i4）;		%對角元1側+勵磁導納?
	else????????%------------否則為線路支路--------------------
		p=B1（i1）;q=B1（i2）;
		Y（pq）=Y（pq）-1./B1（i3）;	%非對角元
????	Y（qp）=Y（pq）;??????????????%非對角元
???		Y（qq）=Y（qq）+1./B1（i3）+B1（i4）./2.0000;	%對角元j側+線路電納的一半
????	Y（pp）=Y（pp）+1./B1（i3）+B1（i4）./2.0000;	%對角元i側+線路電納的一半
	end
end
disp（‘導納矩陣?Y=‘）;disp（Y）;
%-----------給定各節點初始電壓及給定各節點注入功率--------------------------
G=real（Y）;B=imag（Y）;			%分解出導納陣的實部和虛部??
for?i=1:n						%給定各節點初始電壓的實部和虛部????
????e（i）=real（B2（i3））;
????f（i）=imag（B2（i3））;
????V（i）=abs（B2（i3））;				%PV、平衡節點及PQ節點電壓模值?
end
for?i=1:n						%給定各節點注入功率????
????S（i）=B2（i1）-B2（i2）;?		%i節點注入功率SG-SL??
????B（ii）=B（ii）+B2（i4）;		%i節點無功補償量（電納值）
end
%==================用牛頓-拉夫遜法迭代求解非線性代數方程（功率方程）=======================
P=real（S）;Q=imag（S）;????????????%分解出各節點注入的有功和無功功率
ICT1=0;IT2=1;N0=2*n;N1=N0+1;a=0;?%迭代次數ICT1、a；不滿足收斂要求的節點數IT2
while?IT2~=0???????????????????%?N0=2*n?雅可比矩陣的階數；N1=N0+1擴展列
??????IT2=0;a=a+1;
??????JZ=[‘Jacobi矩陣第（‘num2str（a）‘）次消去運算‘];JZ1=[‘Jacobi矩陣第（‘num2str（a）‘）次回代運算‘];JZ0=[‘功率方程第（‘num2str（a）‘）次差值：‘];
????%----------------求取各個節點的功率及功率偏差及PV節點的電