???Ui?ei?jfi ???Yij?Gij?jBij将上式代入Pi?jQi?Ui*?Yj?1nijUj的右端，展开并分出实部和虚部，便得
?nn?Pi?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)??j?1j?1 （1） ?nn?Qi?fi?(Gijej?Bijfj)?ei?(Gijfj?Bijej)?j?1j?1?按照上节的分类，PQ节点的有功功率和无功功率是给定的，第i个节点的给定功率设为Pis和
Qis。假定系统中的第1，2，··· ，m节点为PQ节点，对其中每一个节点可列方程：
nn??Pi?Pis?Pi?Pis?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)?0??j?1j?1 ?nn??Qi?Qis?Qi?Qis?fi?(Gijej?Bijfj)?ei?(Gijfj?Bijej)?0?j?1j?1? (i=1，2，…，m) (2) PV节点的有功功率和节点电压幅值是给定的。假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一节点可以列方程：
nn???Pi?Pis?Pi?Pis?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)?0 ?j?1j?1??U2?U2?U2?U2?(e2?f2)?0isiisii?i (i=m+1，m+2，···，n-1） (3)
第n号节点为平衡点，其电压Un?en?jfn是给定的，故不参加迭代。
式（2)和（3）总共包含了2(n-1)个方程，待求的变量有e1,f1,e2,f2,?,en?1,fn?1也是2（n-1）个。对其近似解的修正量，可以由方程式(4)来确定：
?W??J?U (4） 式中，
???P1???e1???Q???f?1???1????????????P?em???m???Qm???f??W??? ?U??m?
??Pm?1???em?1???U2???f?m?1???m?1?????????P???e?n?1???n?1?2????Un?1????fn?1??14
??P??P??P??P??P??P??P???P?11111111......??e?f1?em?fm?em?1?fm?1?en?1?fn?1?1????Q??Q??Q??Q??Q??Q??Q??Q11111111??......??e1?f1?em?fm?em?1?fm?1?en?1?fn?1????????????????P??P??P??P??P??P??P??Pmmmmmmmm??......??e1?f1?em?fm?em?1?fm?1?en?1?fn?1????Qm??Qm??Qm??Qm??Qm??Qm??Qm??Qm?......???e?f?e?f?e?f?e?f1mmm?1m?1n?1n?1?J????P1??Pm?1??Pm?1??Pm?1??Pm?1??Pm?1??Pm?1??Pm?1? ?m?1......??e?f1?em?fm?em?1?fm?1?en?1?fn?1?1??22222222??U??U??U??U??U??U??U??Um?1??m?1m?1m?1m?1m?1m?1m?1......??e?f1?em?fm?em?fm?1?en?1?fn?1?1??????????????P??P??P??P??P??P??P??Pn?1n?1n?1n?1n?1n?1n?1n?1??......??e1?f1?em?fm?em?1?fm?1?en?1?fn?1???22222222???Un?1??Un?1...??Un?1??Un?1??Un?1??Un?1...??Un?1??Un?1???f1?em?fm?em?fm?1?en?1?fn?1???e1?上述方程中雅可比矩阵的各元素，可以对式（2)和（3）求偏导数获得。当j=i时，对角元素是
n???Pi??e???(Gijej?Bijfj)?Gijej?Biifij?1?in???Pi??(Gf?Be)?Be?Gf?ijjijjiiiiii??fj?1i?n???Qi?(Gf?Be)?Be?Gf?ijjijjiiiiii??eij?1?n??Qi????(Gijej?Bijfj)?Giiei?Biifi??fij?1 （5） 2???Ui???2ei??ei???Ui2??2fi???fi当j?i时，矩阵中非对角元素是
???P??Qii????(Gijei?Bijfi)??fj??ej???P??Q?ii??Bijei?Gijfi （6） ??f?ej?j???U2??U2ii???0?fj???ej由以上表达式不难看出，雅可比矩阵有以下特点：
（1）雅可比矩阵中的诸元素都是节点电压的函数，因此在迭代过程中，它们将随着各节点电压的
变化而不断地改变。 （2）矩阵是不对称的。
由式（6）可以看出，当导纳矩阵中的非对角元素Yij为零时，雅可比矩阵中相对应的元素也是零，即矩阵是非常稀疏的。因此，修正方程的求解同样可以应用稀疏矩阵的求解技巧。正是由于这
15
一点才使牛顿-拉夫逊法获得广泛的应用。
用牛顿法计算潮流时，步骤如下：
(0)(0)（1）给出各节点电压初始值ei、fi。,
（2）将以上电压初始值代入式（2）和（3），求出修正方程式的常数项向量?Pi(0)、?Qi(0)、
?U2i素。
(0)。
（3）将电压初始值代入式（5）和（6），求出修正方程式中系数矩阵（雅可比矩阵）的各元（4）解修正方程式（4），求出修正量?ei（5）修正各节点电压
(0)，?fi(0)。
ei（6）将ei、fi(1)(1)?ei(0)??ei，fi(1)?fi(0)??fi(0)
(1)(0)(1)代入式（2）和（3），求?Pi，?Qi，?U2i(1)(1)。
（7）校验是否收敛，即
如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果；如果不收敛，转回第（2）步进行下一次迭代计算，直到收敛为止。
2. 程序代码及说明 n=input('请输入节点数:n='); nl=input('请输入支路数:nl=');
isb=input('请输入平衡母线节点号:isb='); pr=input('请输入误差精度:pr=');
B1=input('请输入由各支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2=');
X=input('请输入由节点号和接地支路参数形成的矩阵:X=');
Y=zeros(n); e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);S1=zeros(nl); for i=1:n %节点数 if X(i,2)~=0; p=X(i,1);
Y(p,p)=1./X(i,2); %接地支路, “./”点除代表矩阵对应元素相除 end end
for i=1:nl %支路数
if B1(i,6)==0 %折算到哪一侧的标志,0非标准变比在q侧,1非标准变比在在p侧 p=B1(i,1);q=B1(i,2); % B1(i,1), B1(i,2)为支路编号 else
p=B1(i,2);q=B1(i,1); end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元, 节点q的自导纳, B1(i,4)为支路对地容抗
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元,节点p的自导纳 end
16
%输出导纳矩阵 disp('导纳矩阵Y='); disp(Y);
G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部
for i=1:n %给定各节点初始电压的实部和虚部 e(i)=real(B2(i,3)); f(i)=imag(B2(i,3));
V(i)=B2(i,4); %PV节点电压给定模值 end
for i=1:n %给定各节点注入功率 S(i)=B2(i,1)-B2(i,2); B(i,i)=B(i,i)+B2(i,5); end
P=real(S);Q=imag(S);
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; while IT2~=0 IT2=0;a=a+1; for i=1:n if i~=isb C(i)=0; D(i)=0; for j1=1:n
C(i)= C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); D(i)= D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); end
P1=C(i)*e(i)+f(i)*D(i); Q1=f(i)*C(i)-D(i)*e(i); V2=e(i)^2+f(i)^2; if B2(i,6)~=3 DP=P(i)-P1; DQ=Q(i)-Q1; for j1=1:n
if j1~=isb&j1~=i
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X3=X2; X4=-X1;
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ; m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; elseif j1==i&j1~=isb
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
17
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; end end else
DP=P(i)-P1; DV=V(i)^2-V2; for j1=1:n
if j1~=isb&j1~=i
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X5=0; X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2; elseif j1==i&j1~=isb
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X5=-2*e(i); X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; J(m,q)=X2; end end end end end
%求雅可比矩阵 for k=3:N0
k1=k+1;N1=N; for k2=k1:N1
J(k,k2)=J(k,k2)./J(k,k); end
J(k,k)=1; if k~=3; k4=k-1; for k3=3:k4 for k2=k1:N1
J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2); end
J(k3,k)=0;
18