function?[delPQ]?=deltaPQ（Ybus）?????????%計算功率偏差
G=real（Y）;
B=imag（Y）;
[nymy]=size（Y）;
[nbmb]=size（bus）;
nSW=0;????????????????????????%計數：平衡節點個數
nPV=0;????????????????????????%計數：PV節點個數
nPQ=0;????????????????????????%計數：PQ節點個數
for?I=1:nb????????????????????%?nb為總節點數
????type=bus（I6）;
????if?type==3
???????nSW=nSW+1;?????????????%?increment?swing?bus?counter
????elseif?type==2
???????nPV=nPV?+1;????????????%?increment?PV?bus?counter
????else
???????nPQ=nPQ+1;?????????????%?increment?PQ?bus?counter
????end
end
ang=bus（:3）;?????????????????%把bus第3列相角存放到矩陣ang中
V=bus（:2）;???????????????????%把bus第2列相角存放到矩陣V中
A=zeros（nymy）;???????????????%矩陣A存放節點相角差
P=bus（:4）;???????????????????%把bus第4列輸入有功功率存放到矩陣P中
Q=bus（:5）;???????????????????%把bus第5列輸入無功功率存放到矩陣Q中
for?k=1:ny????????????????????%求各節點相角差存放到矩陣A中
????A（k:）=ang（k）-ang（:）;
end
for?j=1:（nPQ+nPV）?????????????%計算PQ、PV節點的有功功率偏差dP存放到矩陣delP中
????delP（j）=P（j）;
????for?k=1:ny
????????delP（j）=delP（j）-V（j）*V（k）*（G（jk）*cos（A（jk））+B（jk）*sin（A（jk）））;
????end
end
for?j=1:nPQ??????????????????%計算PQ節點的無功功率偏差dQ存放到矩陣delQ中
????delQ（j）=Q（j）;
????for?k=1:ny
????????delQ（j）=delQ（j）-V（j）*V（k）*（G（jk）*sin（A（jk））-B（jk）*cos（A（jk）））;
????end
end
delPQ=[delPdelQ]‘;??????????%按dP、dQ節點的順序形成一個單列矩陣存放到矩陣delPQ中

????
????
????
????