(完整word版)潮流計(jì)算方法

1、由于本人參加我們電氣學(xué)院的電氣小課堂， 主講的是計(jì)算機(jī)算法計(jì)算潮流這章， 所以潛 心玩了一個(gè)星期，下面整理給大家分享下。本人一個(gè)星期以來的汗水， 弄清楚了計(jì)算機(jī)算法計(jì)算潮流的基礎(chǔ)， 如果有什么不懂的可 以發(fā)信息到郵箱： zenghao616 接下來開始弄潮流的優(yōu)化問題，吼吼！電力系統(tǒng)的潮流計(jì)算的計(jì)算機(jī)算法：以 MATLAB 為環(huán)境 這里理論不做過多介紹，推薦一本專門講解電力系統(tǒng)分析的計(jì)算機(jī)算法的書籍 電力系統(tǒng)分析的計(jì)算機(jī)算法 邱曉燕、劉天琪編著。這里以這本書上的例題【 2-1 】說明計(jì)算機(jī)算法計(jì)算的過程，分別是牛頓拉弗遜算法的 直角坐標(biāo)和極坐標(biāo)算法、 P-Q 分解算法。主要是簡(jiǎn)單的網(wǎng)絡(luò)的潮流計(jì)

2、算，其實(shí)簡(jiǎn)單網(wǎng)絡(luò)計(jì)算 和大型網(wǎng)絡(luò)計(jì)算并無(wú)本質(zhì)區(qū)別，代碼里面只需要修改循環(huán)迭代的 N 即可，這里旨在弄清計(jì) 算機(jī)算法計(jì)算潮流的本質(zhì)。代碼均有詳細(xì)的注釋 .其中簡(jiǎn)單的高斯賽德爾迭代法是以我們的電穩(wěn)教材為例子講， 其實(shí)都差不多， 只要把導(dǎo) 納矩陣 Y 給你，節(jié)點(diǎn)的編號(hào)和分類給你，就可以進(jìn)行計(jì)算了，不必要找到原始的電氣接線 圖。理論不多說，直接上代碼：簡(jiǎn)單的高斯賽德爾迭代法： 這里我們只是迭代算出各個(gè)節(jié)點(diǎn)的電壓值，支路功率并沒有計(jì)算。S_ij=P_ij+Q_ij=V_i(V_i* - V_j*) * y_ij*可以計(jì)算出各個(gè)線路的功率 在顯示最終電壓幅角的時(shí)候注意在 MATLAB 里面默認(rèn)的是弧度的形

3、式，需要轉(zhuǎn)化成角 度顯示。clear;clc;%電穩(wěn)書 Page 102 例題 3-5%計(jì)算網(wǎng)絡(luò)的潮流分布 -高斯- 賽德爾算法%其中節(jié)點(diǎn) 1 是平衡節(jié)點(diǎn)%節(jié)點(diǎn) 2、3是PV節(jié)點(diǎn)，其余是 PQ節(jié)點(diǎn)% 如果節(jié)點(diǎn)有對(duì)地導(dǎo)納支路%需將對(duì)地導(dǎo)納支路算到自導(dǎo)納里面%輸入原始數(shù)據(jù)，每條支路的導(dǎo)納數(shù)值，包括自導(dǎo)和互導(dǎo)納；y=zeros(5,5);y(1,2)=1/(0.0194+0.0592*1i);y(1,5)=1/(0.054+0.223*1i);y(2,3)=1/(0.04699+0.198*1i);y(2,4)=1/(0.0581+0.1763*1i);%由于電路網(wǎng)絡(luò)的互易性，導(dǎo)納矩陣為對(duì)稱的矩陣f

4、or i=1:1:5for j=1:1:5y(j,i)=y(i,j);endend%節(jié)點(diǎn)導(dǎo)納矩陣的形成Y=zeros(5,5);%求互導(dǎo)納for i=1:1:5for j=1:1:5if i=jY(i,j)=-y(i,j);endendend%求自導(dǎo)納for i=1:1:5%這句話是說將 y 矩陣的第 i 行的所有元素相加，得到自導(dǎo)納的值 Y(i,i)=sum(y(i,:);end %上面求得的自導(dǎo)納不包含該節(jié)點(diǎn)的對(duì)地導(dǎo)納數(shù)值，需要加上 Y(2,2)=Y(2,2)+0.067*1i;Y(3,3)=Y(3,3)+0.022*1i;Y(4,4)=Y(4,4)+0.0187*1i;Y(5,5)=Y(

5、5,5)+0.0246*1i;%導(dǎo)納矩陣的實(shí)部和虛部G = real(Y);B = imag(Y);Qc2=0;Qc3=0;%原始節(jié)點(diǎn)功率 %這里電源功率為正，負(fù)荷功率為負(fù)S(1)=0;S(2)=-0.217-0.121*1i+Qc2*1i;S(3)=-0.749-0.19*1i+Qc3*1i;S(4)=-0.658+0.039*1i;S(5)=-0.076-0.016*1i;%節(jié)點(diǎn)功率的 P QP = real(S);Q = imag(S);%下面是兩個(gè) PV節(jié)點(diǎn)的無(wú)功初始值Q(2) = 0;Q(3) = 0;U=ones(5,1);%1列5行的 1 '矩陣%節(jié)點(diǎn)電壓初始值U(1)=

6、1.06;U(2)=1.045;U(3)=1.01;U_reg=U;Sum_YU0=0; %中間變量Sum_YU1=0; %中間變量for cont=1:1:6%這里的 cont 是迭代次數(shù)for i=2:1:5for j=1:1:i if i=jSum_YU0 = Sum_YU0 + Y(i,j)*U_reg(j);endendfor j=i+1:1:5Sum_YU1 = Sum_YU1 + Y(i,j)*U(j);endU(i)=( (P(i)-Q(i)*1i ) / conj(U(i) - Sum_YU0 - Sum_YU1 ) / Y(i,i); U_reg(i)=U(i);%PV節(jié)點(diǎn)

7、計(jì)算%下面是把求出的 U2、 U3 只保留其相位，幅值不變 if i=2angle_U2 = angle(U(2); U(2)=1.045*cos(angle_U2)+1.045*sin(angle_U2)*1i;Q(2)=imag( U(2)*( conj(Sum_YU0) + conj(Sum_YU1) conj(Y(2,2)*U(2) ) );endif i=3angle_U3 = angle(U(3); U(3)=1.01*cos(angle_U3)+1.01*sin(angle_U3)*1i;+ conj(Sum_YU1)Q(3)=imag( U(3)*( conj(Sum_YU0)

8、 conj(Y(3,3)*U(3) ) );end% 下面做越界檢查%if Q(4)>Q_Max% Q(4) = Q_Max; %end%if Q(4)<Q_Min% Q(4) = Q_Min;%end%下面可以做 PV節(jié)點(diǎn)收斂判斷Sum_YU0 = 0;Sum_YU1 = 0;endend%節(jié)點(diǎn)注入無(wú)功，流入為正，流出為負(fù)Qc2=Q(2)+0.121-1.0452 * 0.067;Qc3=Q(3)+0.19-1.012 * 0.022;%電壓幅值和相角angle_U=angle(U)*180/pi;U=abs(U);S_Line=zeros(5,5);%計(jì)算平衡節(jié)點(diǎn)功率S_Bal

9、anceNode=0;for j=1:1:5S_BalanceNode = S_BalanceNode + U(1) * conj(Y(1,j)*U(j); end %下面由上面算出的電壓值求線路的功率%這里計(jì)算出來的線路功率的有功、無(wú)功%for i=1:1:5% for j=i:1:5% if i=j% S_Line(i,j)=U(i)*( conj(U(i)-conj(U(j) ) * conj(y(i,j);% end% if i=2% %S_Line(2,j)=S_Line(2,j)+U(2)*conj(0.067*1i);% end% if i=3% %S_Line(3,j)=S_L

10、ine(3,j)+U(3)*conj(0.022*1i);% end% end %end計(jì)算網(wǎng)絡(luò)的潮流分布 Newton算法( 直角坐標(biāo) )clear;clc;%電穩(wěn)書 Page 102 例題 3-5%計(jì)算網(wǎng)絡(luò)的潮流分布 Newton 算法( 直角坐標(biāo) ) %其中節(jié)點(diǎn) 1 是平衡節(jié)點(diǎn)%節(jié)點(diǎn) 2、3是PV節(jié)點(diǎn)，其余是 PQ節(jié)點(diǎn)% 如果節(jié)點(diǎn)有對(duì)地導(dǎo)納支路 %需將對(duì)地導(dǎo)納支路算到自導(dǎo)納里面%輸入原始數(shù)據(jù)，每條支路的導(dǎo)納數(shù)值，包括自導(dǎo)和互導(dǎo)納； y=zeros(5,5);y(1,2)=1/(0.0194+0.0592*1i); y(1,5)=1/(0.054+0.223*1i);y(2,3)=1/(0

11、.04699+0.198*1i);y(2,4)=1/(0.0581+0.1763*1i);%由于電路網(wǎng)絡(luò)的互易性，導(dǎo)納矩陣為對(duì)稱的矩陣for i=1:1:5 for j=1:1:5y(j,i)=y(i,j);endend%節(jié)點(diǎn)導(dǎo)納矩陣的形成Y=zeros(5,5);%求互導(dǎo)納for i=1:1:5for j=1:1:5 if i=jY(i,j)=-y(i,j);endendend%求自導(dǎo)納for i=1:1:5%這句話是說將 y 矩陣的第 i 行的所有元素相加，得到自導(dǎo)納的值 Y(i,i)=sum(y(i,:);end%上面求得的自導(dǎo)納不包含該節(jié)點(diǎn)的對(duì)地導(dǎo)納數(shù)值，需要加上Y(2,2)=Y(2,

12、2)+0.067*1i;Y(3,3)=Y(3,3)+0.022*1i;Y(4,4)=Y(4,4)+0.0187*1i;Y(5,5)=Y(5,5)+0.0246*1i;%導(dǎo)納矩陣的實(shí)部和虛部G = real(Y);B = imag(Y);%節(jié)點(diǎn) 2、3需補(bǔ)償?shù)臒o(wú)功Qc2=0;Qc3=0;%原始節(jié)點(diǎn)功率%這里電源功率為正，負(fù)荷功率為負(fù)S(1)=0;S(2)=-0.217-0.121*1i+Qc2*1i;S(3)=-0.749-0.19*1i+Qc3*1i;S(4)=-0.658+0.039*1i;S(5)=-0.076-0.016*1i;%節(jié)點(diǎn)功率的 P QP = real(S);Q = imag

13、(S);%下面是兩個(gè) PV節(jié)點(diǎn)的無(wú)功初始值Q(2) = 0;Q(3) = 0;%給點(diǎn)電壓初始值e=1.06,1.045,1.01,1,1;f=0,0,0,0,0;U=e+f*1i;delta_U=zeros(1,5);delta_P=zeros(1,5);delta_Q=zeros(1,5);delta_PQV=ones(8,1);Sum_GB1=0;Sum_GB2=0;cont=0;while max(delta_PQV > 1e-6),cont=cont+1;%for cont=1:1:3%下面開始計(jì)算 delta_P/delta_Q/delta_Ufor i=2:1:5for j=

14、1:1:5 Sum_GB1=Sum_GB1 + ( G(i,j)*e(j) - B(i,j)*f(j) ); Sum_GB2=Sum_GB2 + ( G(i,j)*f(j) + B(i,j)*e(j) );end delta_P(i)=P(i)-e(i)*Sum_GB1-f(i)*Sum_GB2;if i=2 && i=3%不為節(jié)點(diǎn) 2,3 則計(jì)算無(wú)功delta_Q(i)=Q(i)-f(i)*Sum_GB1+e(i)*Sum_GB2;endif i=2 | i=3%這里計(jì)算 delta_U 的值，始終為零delta_U(i)=U(i)2-( e(i)2 + f(i)2 );en

15、dSum_GB1=0;Sum_GB2=0;end %下面計(jì)算雅克比矩陣 J=zeros(8,8);for ii=2:1:5i=ii-1;for j=1:1:5Sum_GB1=Sum_GB1 + ( G(ii,j)*e(j) - B(ii,j)*f(j) ); Sum_GB2=Sum_GB2 + ( G(ii,j)*f(j) + B(ii,j)*e(j) );end for jj=2:1:5 j=jj-1;if ii=2 && ii=3%PQ節(jié)點(diǎn)if ii=jjJ(2*i-1,2*i-1)=-Sum_GB1-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);J(2*i-

16、1,2*i)=-Sum_GB2+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);J(2*i,2*i-1)=Sum_GB2+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);J(2*i,2*i)=-Sum_GB1+G(ii,ii)*e(ii)+B(ii,ii)*f(ii); elseJ(2*i-1,2*j-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii); J(2*i-1,2*j)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii); J(2*i,2*j-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii); J(2*i,2*j)=(

17、G(ii,jj)*e(ii)+B(ii,jj)*f(ii); end else%PV節(jié)點(diǎn)if ii=jjJ(2*i-1,2*i-1)=-Sum_GB1-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);J(2*i-1,2*i)=-Sum_GB2+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);J(2*i,2*i-1)=-2*e(ii);J(2*i,2*i)=-2*f(ii);elseJ(2*i-1,2*j-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii);J(2*i-1,2*j)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);J(2*i,

18、2*j-1)=0;J(2*i,2*j)=0;endendendSum_GB1=0;Sum_GB2=0;end全部放在一個(gè)矩陣%在求解修正方程之前建議把 delta_P 和 delta_Q,delta_U delta_PQV=delta_P(2);delta_U(2);delta_P(3);delta_U(3);delta_P(4) ;delta_Q(4);delta_P(5);delta_Q(5);%下面求解修正方程 ; 注意矩陣運(yùn)算時(shí)候的左除和右除的區(qū)別 delta_ef=-Jdelta_PQV;%下面修正各個(gè)節(jié)點(diǎn)的電壓for i=2:1:5e(i)=e(i)+delta_ef(2*(i-1

19、)-1);f(i)=f(i)+delta_ef(2*(i-1);end%到這里第一輪迭代完成end%電壓幅值和相角U=e+f*1i;angle_U=angle(U)*180/pi;%節(jié)點(diǎn)注入無(wú)功，流入為正，流出為負(fù)Sum_YU=0;for i=2:1:3for j=1:1:5Sum_YU = Sum_YU + Y(i,j)*U(j);endQ(i)=imag( U(i)*conj( Sum_YU ) );Sum_YU=0;endQc2=Q(2)+0.121-1.0452 * 0.067;Qc3=Q(3)+0.19-1.012 * 0.022;U=abs(U);num2str(cont);dis

20、p( 'Iteration times :%顯示最終的迭代次數(shù)牛頓算法求解潮流 ( 極坐標(biāo) ):clear;clc;%牛頓算法求解潮流 ( 極坐標(biāo) ) %計(jì)算網(wǎng)絡(luò)的潮流分布%其中節(jié)點(diǎn) 5 是平衡節(jié)點(diǎn)%節(jié)點(diǎn) 1、2、3是PQ節(jié)點(diǎn)，節(jié)點(diǎn) 4是PV節(jié)點(diǎn)% 如果節(jié)點(diǎn)有對(duì)地導(dǎo)納支路 %需將對(duì)地導(dǎo)納支路算到自導(dǎo)納里面%輸入原始數(shù)據(jù)，每條支路的導(dǎo)納數(shù)值，包括自導(dǎo)和互導(dǎo)納； Y=0.8381-3.7899*1i,-0.4044+1.6203*1i,0,0,-0.4337+2.2586*1i;-0.4044+1.6203*1i,0.7769-3.3970*1i,-0.3726+1.8557*1i,0,0

21、;0,-0.3726+1.8557*1i,1.1428-7.0210*1i,-0.5224+4.1792*1i,-0.2739+1. 2670*1i; .0,0,-0.5224+4.1792*1i,0.5499-4.3591*1i,0; . -0.4337+2.2586*1i,0,-0.2739+1.2670*1i,0,0.7077-3.4437*1i;%導(dǎo)納矩陣的實(shí)部和虛部G = real(Y);B = imag(Y);%給點(diǎn)電壓初始值U = 1,1,1,1,1.05; angle_U=0,0,0,0,0;%for i=1:1:5% U(i)=U_abs(i)*cos(angle_U(i)+

22、U_abs(i)*sin(angle_U(i)*1i; %end%原始節(jié)點(diǎn)功率 %這里電源功率為正，負(fù)荷功率為負(fù) %下面給點(diǎn) PQ PV 節(jié)點(diǎn)功率值 S=-0.22-0.14*1i,-0.18-0.09*1i,-0.27-0.13*1i,0.35,0;%節(jié)點(diǎn)功率的 P QP = real(S);Q = imag(S);%下面是 PV 節(jié)點(diǎn)的無(wú)功初始值Q(4) = 0;delta_P=zeros(1,5);delta_Q=zeros(1,5);%delta_angleU=zeros(1,4);%delta_absU=zeros(1,4);delta_PQ=ones(8,1);Sum_GB1=0;

23、Sum_GB2=0;cont=0;%最外層循環(huán)， cont 代表迭代的次數(shù)，這里可以用約束條件來代替%for cont=1:1:4while max(delta_PQ)>1e-6,%下面計(jì)算 delta_P/delta_Q/delta_Ucont=cont+1;for i=1:1:4for j=1:1:5Sum_GB1=Sum_GB1 + U(j)*(G(i,j)*cos(angle_U(i)-angle_U(j) +B(i,j)*sin(angle_U(i)-angle_U(j) );Sum_GB2=Sum_GB2 + U(j)*(G(i,j)*sin(angle_U(i)-angle

24、_U(j) -B(i,j)*cos(angle_U(i)-angle_U(j) );enddelta_P(i)=P(i)-U(i)*Sum_GB1;if i=4%不為節(jié)點(diǎn)四則計(jì)算無(wú)功delta_Q(i)=Q(i)-U(i)*Sum_GB2; endSum_GB1=0;Sum_GB2=0;end%下面計(jì)算雅克比矩陣J=zeros(7,7);for ii=1:1:4for jj=1:1:4if ii = 4%PQ節(jié)點(diǎn)if ii=jjJ(2*ii-1,2*ii-1)=U(ii)2*B(ii,ii)+Q(ii);J(2*ii-1,2*ii)=-U(ii)2*G(ii,ii)-P(ii);J(2*ii,

25、2*ii-1)=U(ii)2*G(ii,ii)-P(ii);J(2*ii,2*ii)=U(ii)2*B(ii,ii)-Q(ii);elseJ(2*ii-1,2*jj-1)=-U(ii)*U(jj)*( G(ii,jj)*sin(angle_U(ii)-angle_U( jj) - B(ii,jj)*cos(angle_U(ii)-angle_U(jj) );J(2*ii-1,2*jj)=-U(ii)*U(jj)*( G(ii,jj)*cos(angle_U(ii)-angle_U(jj)+ B(ii,jj)*sin(angle_U(ii)-angle_U(jj) );J(2*ii,2*jj-1

26、)=U(ii)*U(jj)*( G(ii,jj)*cos(angle_U(ii)-angle_U(jj)+ B(ii,jj)*sin(angle_U(ii)-angle_U(jj) );J(2*ii,2*jj)=-U(ii)*U(jj)*( G(ii,jj)*sin(angle_U(ii)-angle_U(jj) B(ii,jj)*cos(angle_U(ii)-angle_U(jj) );endelse%PV節(jié)點(diǎn)if ii=jjJ(2*ii-1,2*ii-1)=U(ii)2*B(ii,ii)+Q(ii); J(2*ii-1,2*ii)=-U(ii)2*G(ii,ii)-P(ii);elseJ

27、(2*ii-1,2*jj-1)=-U(ii)*U(jj)*( G(ii,jj)*sin(angle_U(ii)-angle_U( jj) - B(ii,jj)*cos(angle_U(ii)-angle_U(jj) );J(2*ii-1,2*jj)=-U(ii)*U(jj)*( G(ii,jj)*cos(angle_U(ii)-angle_U(jj)+ B(ii,jj)*sin(angle_U(ii)-angle_U(jj) );endendendend%在求解修正方程之前建議把 delta_ef 和 delta_ef 全部放在一個(gè)矩陣 delta_PQ=delta_P(1);delta_Q(

28、1);delta_P(2);delta_Q(2);delta_P(3); delta_Q(3);delta_P(4);%下面求解修正方程 ; 注意矩陣運(yùn)算時(shí)候的左除和右除的區(qū)別 J=J(1:7,1:7);delta_ef=-Jdelta_PQ;%下面修正各個(gè)節(jié)點(diǎn)的電壓for i=1:1:4 if i=4U(i)=U(i)+delta_ef(2*i)*U(i);endangle_U(i)=angle_U(i)+delta_ef(2*i-1);end%到這里第一輪迭代完成end%下面顯示出滿足條件后的迭代的次數(shù)disp( 'Iteration times : 'num2str(co

29、nt);%下面計(jì)算平衡節(jié)點(diǎn) 5的功率 PQfor j=1:1:5Sum_GB1=Sum_GB1+ U(j)*(G(5,j)*cos(angle_U(5)-angle_U(j)+B(5,j)*sin(angle_U(5)-angle_U(j) );Sum_GB2=Sum_GB2+ U(j)*(G(5,j)*sin(angle_U(5)-angle_U(j)-B(5,j)*cos(angle_U(5)-angle_U(j) );endP(5)=U(5)*Sum_GB1;Q(5)=U(5)*Sum_GB2;%下面將相角用角度表示for i=1:1:5angle_U(i)=angle_U(i)*180

30、/pi;End計(jì)算計(jì)算法 P-Q 算法計(jì)算潮流： 這個(gè)算法是由牛頓算法的極坐標(biāo)形式簡(jiǎn)化而來。 clear;clc;%牛頓算法求解潮流 %計(jì)算網(wǎng)絡(luò)的潮流分布 %其中節(jié)點(diǎn) 5 是平衡節(jié)點(diǎn)%節(jié)點(diǎn) 1、2、3是PQ節(jié)點(diǎn)，節(jié)點(diǎn) 4是PV節(jié)點(diǎn)% 如果節(jié)點(diǎn)有對(duì)地導(dǎo)納支路 %需將對(duì)地導(dǎo)納支路算到自導(dǎo)納里面%輸入原始數(shù)據(jù)，每條支路的導(dǎo)納數(shù)值，包括自導(dǎo)和互導(dǎo)納；Y=0.8381-3.7899*1i,-0.4044+1.6203*1i,0,0,-0.4337+2.2586*1i; -0.4044+1.6203*1i,0.7769-3.3970*1i,-0.3726+1.8557*1i,0,0;0,-0.3726+1

31、.8557*1i,1.1428-7.0210*1i,-0.5224+4.1792*1i,-0.2739+1. 2670*1i; .0,0,-0.5224+4.1792*1i,0.5499-4.3591*1i,0; . -0.4337+2.2586*1i,0,-0.2739+1.2670*1i,0,0.7077-3.4437*1i;%導(dǎo)納矩陣的實(shí)部和虛部G = real(Y);B = imag(Y);%給點(diǎn)電壓初始值U = 1,1,1,1,1.05; angle_U=0,0,0,0,0;%原始節(jié)點(diǎn)功率%這里電源功率為正，負(fù)荷功率為負(fù)%下面給點(diǎn) PQ PV 節(jié)點(diǎn)功率值 S=-0.22-0.14*1

32、i,-0.18-0.09*1i,-0.27-0.13*1i,0.35,0;%節(jié)點(diǎn)功率的 P QP = real(S);Q = imag(S);%下面是 PV 節(jié)點(diǎn)的無(wú)功初始值Q(4) = 0;%下面計(jì)算出無(wú)功和有功迭代的系數(shù)矩陣B_P=B(1:4,1:4);%有功迭代系數(shù)矩陣 n-1 即除開平衡節(jié)點(diǎn)以外的所有節(jié)點(diǎn)B_Q=B(1:3,1:3);%無(wú)功迭代系數(shù)矩陣 m 個(gè) 即PQ節(jié)點(diǎn)的個(gè)數(shù)%下面是相關(guān)變量的定義 delta_P=ones(1,5); delta_Q=ones(1,5);delta_PQ=ones(8,1);delta_angle=ones(4,1);Sum_GB1=0;Sum_GB

33、2=0;cont=0;%最外層循環(huán)， cont 代表迭代的次數(shù)，這里可以用約束條件來代替%for cont=1:1:5while max(delta_angle) > 1e-6,%下面計(jì)算 delta_P/delta_Q/delta_Ucont=cont+1;for i=1:1:4for j=1:1:5Sum_GB1=Sum_GB1 + U(j)*(G(i,j)*cos(angle_U(i)-angle_U(j) +B(i,j)*sin(angle_U(i)-angle_U(j) );Sum_GB2=Sum_GB2 + U(j)*(G(i,j)*sin(angle_U(i)-angle_

34、U(j) -B(i,j)*cos(angle_U(i)-angle_U(j) );enddelta_P(i)=P(i)-U(i)*Sum_GB1;if i=4%不為 PV 節(jié)點(diǎn)四則計(jì)算無(wú)功delta_Q(i)=Q(i)-U(i)*Sum_GB2;endSum_GB1=0;Sum_GB2=0;end%下面計(jì)算 delta_P(i)/U(i) delta_Q(i)/U(i) delta_QV=zeros(3,1);delta_PV=zeros(4,1);for i=1:1:4if i=4 delta_QV(i)=delta_Q(i)/U(i);enddelta_PV(i)=delta_P(i)/U(i);end%下面計(jì)算 del