clear;clc;
a=[%首節點t1	末節點t2	電阻r?????電抗x???容納之半ban??變比k
1		????	2		0.010		0.085		0.088	1;
1				3		0.017		0.092		0.079	1;
2				4		0.032		0.161		0.153	1;
3				6		0.039		0.170		0.179	1;
4				5		0.0085		0.072		0.0745??1;
5		????	6		0.0119		0.1008		0.1045	1;
9				1		0			0.0576		0		1;
7				4		0			0.0625		0		1;
8				6		0			0.0586		0		1;
];
b=[%節點	電壓幅值?電壓相角	?節點有功	節點無功	節點類型?/“1”為PQ節點，“2”為PV節點，“3”為平衡節點
1		1.0	?		0		0			0		1;
2		1.0	?		0		-1.2500		-0.5	1;
3		1.0	?		0		-0.9		-0.3	1;
4		1.0	?		0		0			0		1;
5		1.0	?		0		-1			-0.35	1;
6		1.0??????	0		0			0		1;
7		1.025	?	0		1.6300		0		2;
8		1.025	?	0		0.8500		0		2;
9		1.040	?	0		0			0		3;
];

%?a=load（‘zhilu.txt‘）;
%?b=load（‘jiedian.txt‘）;
N=size（a1）;%支路數
Nbus=size（b1）;%節點數
pq=0;
for?k=1:Nbus?%PQ節點個數
????if?b（k6）==1?
????????pq=pq+1;
????end
end
Y=zeros（Nbus）;%節點導納矩陣（生成一個Nbus*Nbus的零矩陣）
for?k=1:N
????t1=a（k1）;t2=a（k2）;r=a（k3）;x=a（k4）;ban=a（k5）;K=a（k6）;
????Y（t1t1）=Y（t1t1）+1/（r+j*x）+j*ban;?%電力工程P175
????Y（t1t2）=Y（t1t2）-1/（K*（r+j*x））;
????Y（t2t1）=Y（t2t1）-1/（K*（r+j*x））;
????Y（t2t2）=Y（t2t2）+1/（K*K*（r+j*x））+j*ban;
end
G=real（Y）;B=imag（Y）;??????%提取導納矩陣實部虛部
?
precision=1;???????????%誤差精度
t=0?;%存儲迭代次數
while?precision>0.00001???%只要精度一直大于0.00001就循環
????P=zeros（Nbus1）;%存儲所有節點的有功
????Q=zeros（Nbus1）;%存儲所有節點的無功
????for?m=1:Nbus?%求pv、pq和平衡節點的有功
????????for?n=1:Nbus
????????????P（m1）=P（m1）+b（m2）*b（n2）*（G（mn）*cos（b（m3）-b（n3））+B（mn）*sin（b（m3）-b（n3）））;
????????%Pm=Um*E（1:n）Un[Gmn*cos（&m-&n）+Bmn*sin（&m-&n）]
????????end
????end
????for?m=1:Nbus?%求pq、pv和平衡節點的無功
????????for?n=1:Nbus
????????????Q（m1）=Q（m1）+b（m2）*b（n2）*（G（mn）*sin（b（m3）-b（n3））-B（mn）*cos（b（m3）-b（n3）））;
???????end
????end
????deltp=b（1:Nbus-14）-P（1:Nbus-11）;%pq和pv節點有功差（除去第9行平衡節點）
????deltq=b（1:pq5）-Q（1:pq1）;%pq節點無功差（只有pq有Q）
????deltPQ=[deltp;deltq];?%功率不平衡量
????H=zeros（Nbus-1）;%極坐標形式不用計算平衡節點
????for?m=1:Nbus-1?%求H矩陣（雅可比HNJL）
????????for?n=1:Nbus-1
????????????if?m~=n??%不等于，非對角元素，對&求偏導
????????????????H（mn）=-b（m2）*b（n2）*（G（mn）*sin（b（m3）-b（n3））-B（mn）*cos（b（m3）-b（n3）））;
????????????else
????????????????H（mm）=b（m2）*b（m2）*B（mm）+Q（m1）;%陳珩書P84
????????????end
????????end
????end
????N=zeros（Nbus-1pq）;%求N矩陣
????for?m=1:Nbus-1?
????????for?n=1:pq
????????????if?m~=n
????????????????N（mn）=-b（m2）*b（n2）*（G（mn）*cos（b（m3）-b（n3））+B（mn）*sin（b（m3）-b（n3）））;
????????????else
????????????????N（mm）=-b（m2）*b（m2）*G（mm）-P（m1）;
????????????end
????????end
????end
????J=zeros（pqNbus-1）;
????for?m=1:pq?%求J矩陣
????????for?n=1:Nbus-1
????????????if?m~=n
????????????????J（mn）=b（m2）*b（n2）*（G（mn）*cos（b（m3）-b（n3））+B（mn）*sin（b（m3）-b（n3）））;
????????????else
????????????????J（mm）=b（m2）*b（m2）*G（mm）-P（m1）;
????????????end
????????end
????end
????L=zeros（pqpq）;
????for?m=1:pq?%求L矩陣
????????for?n=1:pq
????????????if?m~=n
????????????????L（mn）=-b（m2）*b（n2）*（G（mn）*sin（b（m3）-b（n3））-B（mn）*cos（b（m3）-b（n3）））;
????????????else
?

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????3731??2018-08-06?16:34??3機9節點\motor3bus9.m

?????文件?????418408??2018-08-06?16:48??3機9節點\牛拉法潮流計算程序（附3機9節點結果對比）.docx

?????目錄??????????0??2018-08-06?16:48??3機9節點

-----------?---------??----------?-----??----

???????????????422139????????????????????3