【正文】 %非對(duì)角元 Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)./。 %對(duì)角元j側(cè)+線路電納的一半 Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./。 %對(duì)角元i側(cè)+線路電納的一半 endenddisp(39。導(dǎo)納矩陣 Y=39。)。disp(Y)。%給定各節(jié)點(diǎn)初始電壓及給定各節(jié)點(diǎn)注入功率G=real(Y)。B=imag(Y)。 %分解出導(dǎo)納陣的實(shí)部和虛部 for i=1:n %給定各節(jié)點(diǎn)初始電壓的實(shí)部和虛部 e(i)=real(B2(i,3))。 f(i)=imag(B2(i,3))。 V(i)=abs(B2(i,3))。 %PV、平衡節(jié)點(diǎn)及PQ節(jié)點(diǎn)電壓模值 endfor i=1:n %給定各節(jié)點(diǎn)注入功率 S(i)=B2(i,1)B2(i,2)。 %i節(jié)點(diǎn)注入功率SGSL B(i,i)=B(i,i)+B2(i,4)。 %i節(jié)點(diǎn)無功補(bǔ)償量(電納值)end%==================用牛頓拉夫遜法迭代求解非線性代數(shù)方程（功率方程）=======================P=real(S)。Q=imag(S)。 %分解出各節(jié)點(diǎn)注入的有功和無功功率ICT1=0。IT2=1。N0=2*n。N1=N0+1。a=0。 %迭代次數(shù)ICTa；不滿足收斂要求的節(jié)點(diǎn)數(shù)IT2while IT2~=0 % N0=2*n 雅可比矩陣的階數(shù)；N1=N0+1擴(kuò)展列 IT2=0。a=a+1。 JZ=[39。Jacobi矩陣第(39。,num2str(a),39。)次消去運(yùn)算39。]。JZ1=[39。Jacobi矩陣第(39。,num2str(a),39。)次回代運(yùn)算39。]。JZ0=[39。功率方程第(39。,num2str(a),39。)次差值：39。]。 %求取各個(gè)節(jié)點(diǎn)的功率及功率偏差及PV節(jié)點(diǎn)的電壓偏差 for i=1:n %n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) p=2*i1。m=p+1。C(i)=0。D(i)=0。 for j1=1:n %第i行共n列(n個(gè)節(jié)點(diǎn)間互導(dǎo)納及節(jié)點(diǎn)電壓相乘即電流) C(i)=C(i)+G(i,j1)*e(j1)B(i,j1)*f(j1)。%Σ(Gij*ejBij*fj) D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1)。%Σ(Gij*fj+Bij*ej) end %求i節(jié)點(diǎn)有功和無功功率P39。,Q39。的計(jì)算值 P1=C(i)*e(i)+f(i)*D(i)。%節(jié)點(diǎn)功率P計(jì)算eiΣ(Gij*ejBij*fj)+fiΣ(Gij*fj+Bij*ej) Q1=C(i)*f(i)e(i)*D(i)。%節(jié)點(diǎn)功率Q計(jì)算fiΣ(Gij*ejBij*fj)eiΣ(Gij*fj+Bij*ej) V2=e(i)^2+f(i)^2。 %電壓模平方 %===求取功率差及PV節(jié)點(diǎn)電壓模平方差 ========= if i~=isb %非平衡節(jié)點(diǎn)(PQ或PV節(jié)點(diǎn)) if B2(i,5)~=3 %非PV節(jié)點(diǎn)(只能是PQ節(jié)點(diǎn)) J(m,N1)=P(i)P1。 %PQ節(jié)點(diǎn)有功功率差J(m,N1)擴(kuò)展列△P J(p,N1)=Q(i)Q1。 %PQ節(jié)點(diǎn)無功功率差J(p,N1)擴(kuò)展列△Q else %PV節(jié)點(diǎn)================== J(m,N1)=P(i)P1。 %PV節(jié)點(diǎn)有功功率差J(m,N1)擴(kuò)展列△P J(p,N1)=V(i)^2V2。 %PV節(jié)點(diǎn)電壓模平方差J(p,N1)擴(kuò)展列△U end end %(if i~=isb) 非平衡節(jié)點(diǎn)(PQ或PV節(jié)點(diǎn)) end %(for i=1:n) n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) for m=1:N0 JJN1(m)=J(m,N1)。 end disp(JZ0)。disp(JJN1)。 %判斷功率偏差量及PV節(jié)點(diǎn)的電壓偏差量是否滿足要求 for k=3:N0 %除去平衡節(jié)點(diǎn)2號(hào)以外的所有節(jié)點(diǎn) DET=abs(J(k,N1))。 if DET=pr。 %PQ節(jié)點(diǎn)的功率偏差量及PV節(jié)點(diǎn)的電壓偏差量是否滿足要求 IT2=IT2+1。 %不滿足要求的節(jié)點(diǎn)數(shù)加1 end end ICT2(a)=IT2。 %不滿足要求的節(jié)點(diǎn)數(shù)；a為迭代次數(shù) ICT1=ICT1+1。 %迭代次數(shù) if ICT2(a)==0。 %當(dāng)前不滿足要求的節(jié)點(diǎn)數(shù)為零 break %退出迭代運(yùn)算 end %以上為求取各個(gè)節(jié)點(diǎn)的功率及功率偏差及PV節(jié)點(diǎn)的電壓偏差 %================= 求取Jacobi矩陣形成修正方程 =================== for i=2:n %n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) if i~=isb %非平衡節(jié)點(diǎn)(PQ或PV節(jié)點(diǎn)) if B2(i,5)~=3 %下面是針對(duì)PQ節(jié)點(diǎn)來求取Jacobi矩陣的元素 =========== C(i)=0。D(i)=0。 for j1=1:n %第i行共n列(n個(gè)節(jié)點(diǎn)間互導(dǎo)納及節(jié)點(diǎn)電壓相乘即電流) C(i)=C(i)+G(i,j1)*e(j1)B(i,j1)*f(j1)。%Σ(Gij*ejBij*fj) D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1)。%Σ(Gij*fj+Bij*ej) end for j1=2:n %第i行共n列(2n個(gè)Jacobi矩陣元素dP/de及dP/df或dQ/de及dQ/df) if j1~=isbamp。j1~=i %非平衡節(jié)點(diǎn)amp。非對(duì)角元 X1=G(i,j1)*e(i)B(i,j1)*f(i)。 % X1=dP/de=dQ/df=X4 X2=B(i,j1)*e(i)G(i,j1)*f(i)。 % X2=dP/df=dQ/de=X3 X3=X2。 % X2=dp/df X3=dQ/de X4=X1。 % X1=dP/de X4=dQ/df p=2*i1。q=2*j11。 J(p,q)=X3。m=p+1。 % X3=dQ/de J(p,N)=DQ節(jié)點(diǎn)無功功率差 J(p,N)=DQ。 J(m,q)=X1。q=q+1。 % X1=dP/de J(m,N)=DP節(jié)點(diǎn)有功功率差 J(m,N)=DP。 J(p,q)=X4。J(m,q)=X2。 % X4=dQ/df X2=dp/df elseif j1==iamp。j1~=isb %非平衡節(jié)點(diǎn)amp。對(duì)角元 X1=C(i)G(i,i)*e(i)B(i,i)*f(i)。% dP/de X2=D(i)+B(i,i)*e(i)G(i,i)*f(i)。% dP/df X3=D(i)+B(i,i)*e(i)G(i,i)*f(i)。 % dQ/de X4=C(i)+G(i,i)*e(i)+B(i,i)*f(i)。% dQ/df p=2*i1。q=2*j11。J(p,q)=X3。%擴(kuò)展列△Q J(p,N)=DQ。 m=p+1。 J(m,q)=X1。q=q+1。J(p,q)=X4。%擴(kuò)展列△P J(m,N)=DP。 J(m,q)=X2。 end end else %if B2(i,5)~=3 % 否則（即為PV節(jié)點(diǎn)） %=============== 下面是針對(duì)PV節(jié)點(diǎn)來求取Jacobi矩陣的元素 =========== for j1=1:n if j1~=isbamp。j1~=i %非平衡節(jié)點(diǎn)amp。非對(duì)角元 X1=G(i,j1)*e(i)B(i,j1)*f(i)。 % dP/de X2=B(i,j1)*e(i)G(i,j1)*f(i)。 % dP/df X5=0。X6=0。 p=2*i1。q=2*j11。J(p,q)=X5。 % PV節(jié)點(diǎn)電壓誤差J(p,N)=DV。 m=p+1。 J(m,q)=X1。q=q+1。J(p,q)=X6。 % PV節(jié)點(diǎn)有功誤差J(m,N)=DP。 J(m,q)=X2。 elseif j1==iamp。j1~=isb %非平衡節(jié)點(diǎn)amp。對(duì)角元 X1=C(i)G(i,i)*e(i)B(i,i)*f(i)。% dP/de X2=D(i)+B(i,i)*e(i)G(i,i)*f(i)。% dP/df X5=2*e(i)。 X6=2*f(i)。 p=2*i1。q=2*j11。J(p,q)=X5。 % PV節(jié)點(diǎn)電壓誤差J(p,N)=DV。 m=p+1。 J(m,q)=X1。q=q+1。J(p,q)=X6。 % PV節(jié)點(diǎn)有功誤差J(m,N)=DP。 J(m,q)=X2。 end end end %(if B2(i,5)~=3 else) end %(if i~=isb) end %(for i=1:n)n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) JZ0=[39。形成的第(39。,num2str(a),39。)次Jacobi矩陣：39。]。 disp(JZ0)。disp(J)。 %=============================== 以上為形成完整的Jacobi矩陣 ================================ %====下面用高斯消去法對(duì)由Jacobi矩陣形成的修正方程進(jìn)行求解(按列消去、回代) ========== for k=3:N0 % N0=2*n （從第三行開始，第一、二行是平衡節(jié)點(diǎn)） for k1=k+1:N1 % 從k+1列的Jacobi元素到擴(kuò)展列的△P、△Q 或 △U J(k,k1)=J(k,k1)./J(k,k)。% 用K行K列對(duì)角元素去除K行