计算机潮流计算资料

方程数： n 1 m n m 1 2(n 1)
i
未知量： ei , f i , i ( PQ PV ) , 2(n 1)
第四讲 电力系统潮流计算
（2）极坐标下的数学方程
将 Vi Vi e
j i
和 Yij Gij jBij 代入 Pi jQi Vi Yij V j
n个节点，m个PQ节点，n m 1个PV节点
编号
1 m
PQ
m 1 n 1
PV
n
平衡点
I YV
I i Yij V j
j 1 n
Vi I i Vi Y ij V j
j 1

n


Pi jQi Vi Yij V j
x(k )
x
讨论收敛区域和收敛条件。又称切线法。
第四讲 电力系统潮流计算
2、设初始点 xo , f ( xo ) 0
f ( xo x ) 0 df f ( xo ) dx f ( xo ) x df dx 1 d2 f x 2 dx 2 x 0
x0
x 2
( 0) x1 0 ( 0) x 2
(1) ( 0) ( 0) x1 x1 x1 (1) ( 0) ( 0) x2 x2 x2
第四讲 电力系统潮流计算
3、多维非线性方程组的迭代公式
基于同样的思想，我们可以得到 n维非线性方程 —牛顿 拉夫逊迭代公式
j 1
j i j i ( G jB ) V e ij ij i j 1 n

n
*

由于本人参加我们电气学院的电气小课堂，主讲的是计算机算法计算潮流这章，所以潜心玩了一个星期，下面整理给大家分享下。

本人一个星期以来的汗水，弄清楚了计算机算法计算潮流的基础，如果有什么不懂的可以发信息到邮箱：zenghao616@接下来开始弄潮流的优化问题，吼吼！电力系统的潮流计算的计算机算法：以MATLAB为环境这里理论不做过多介绍，推荐一本专门讲解电力系统分析的计算机算法的书籍---------《电力系统分析的计算机算法》—邱晓燕、刘天琪编著。

这里以这本书上的例题【2-1】说明计算机算法计算的过程，分别是牛顿拉弗逊算法的直角坐标和极坐标算法、P-Q分解算法。

主要是简单的网络的潮流计算，其实简单网络计算和大型网络计算并无本质区别，代码里面只需要修改循环迭代的N即可，这里旨在弄清计算机算法计算潮流的本质。

代码均有详细的注释.其中简单的高斯赛德尔迭代法是以我们的电稳教材为例子讲，其实都差不多，只要把导纳矩阵Y给你，节点的编号和分类给你，就可以进行计算了，不必要找到原始的电气接线图。

理论不多说，直接上代码：简单的高斯赛德尔迭代法：这里我们只是迭代算出各个节点的电压值，支路功率并没有计算。

S_ij=P_ij+Q_ij=V_i(V_i* - V_j*) * y_ij*可以计算出各个线路的功率在显示最终电压幅角的时候注意在MATLAB里面默认的是弧度的形式，需要转化成角度显示。

clear;clc;%电稳书Page 102 例题3-5%计算网络的潮流分布 --- 高斯-赛德尔算法%其中节点1是平衡节点%节点2、3是PV节点，其余是PQ节点% 如果节点有对地导纳支路%需将对地导纳支路算到自导纳里面%------------------------------------------------%%输入原始数据，每条支路的导纳数值，包括自导和互导纳；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);%由于电路网络的互易性，导纳矩阵为对称的矩阵for i=1:1:5for j=1:1:5y(j,i)=y(i,j);endend%节点导纳矩阵的形成Y=zeros(5,5);%求互导纳for i=1:1:5for j=1:1:5if i~=jY(i,j)=-y(i,j);endendend%求自导纳for i=1:1:5%这句话是说将y矩阵的第i行的所有元素相加，得到自导纳的值 Y(i,i)=sum(y(i,:));end%上面求得的自导纳不包含该节点的对地导纳数值，需要加上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)+0.0246*1i;%导纳矩阵的实部和虚部G = real(Y);B = imag(Y);Qc2=0;Qc3=0;%原始节点功率%这里电源功率为正，负荷功率为负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;%节点功率的P QP = real(S);Q = imag(S);%下面是两个PV节点的无功初始值Q(2) = 0;Q(3) = 0;U=ones(5,1); %1列5行的‘1’矩阵%节点电压初始值U(1)=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是迭代次数for i=2:1:5for j=1:1:iif 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节点计算%下面是把求出的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;Q(3)=imag( U(3)*( conj(Sum_YU0) + conj(Sum_YU1) + 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节点收敛判断Sum_YU0 = 0;Sum_YU1 = 0;endend%节点注入无功，流入为正，流出为负Qc2=Q(2)+0.121-1.045^2 * 0.067;Qc3=Q(3)+0.19-1.01^2 * 0.022;%电压幅值和相角angle_U=angle(U)*180/pi;U=abs(U);S_Line=zeros(5,5);%计算平衡节点功率S_BalanceNode=0;for j=1:1:5S_BalanceNode = S_BalanceNode + U(1) * conj(Y(1,j)*U(j));end%下面由上面算出的电压值求线路的功率%这里计算出来的线路功率的有功、无功%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_Line(3,j)+U(3)*conj(0.022*1i);% end% end%end计算网络的潮流分布 ---- Newton算法(直角坐标)clear;clc;%电稳书Page 102 例题3-5%计算网络的潮流分布 ---- Newton算法(直角坐标)%其中节点1是平衡节点%节点2、3是PV节点，其余是PQ节点% 如果节点有对地导纳支路%需将对地导纳支路算到自导纳里面%------------------------------------------------%%输入原始数据，每条支路的导纳数值，包括自导和互导纳；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);%由于电路网络的互易性，导纳矩阵为对称的矩阵for i=1:1:5for j=1:1:5y(j,i)=y(i,j);endend%节点导纳矩阵的形成Y=zeros(5,5);%求互导纳for i=1:1:5for j=1:1:5if i~=jY(i,j)=-y(i,j);endendend%求自导纳for i=1:1:5%这句话是说将y矩阵的第i行的所有元素相加，得到自导纳的值Y(i,i)=sum(y(i,:));end%上面求得的自导纳不包含该节点的对地导纳数值，需要加上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)+0.0246*1i;%导纳矩阵的实部和虚部G = real(Y);B = imag(Y);%节点2、3需补偿的无功Qc2=0;Qc3=0;%原始节点功率%这里电源功率为正，负荷功率为负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;%节点功率的P QP = real(S);Q = imag(S);%下面是两个PV节点的无功初始值Q(2) = 0;Q(3) = 0;%给点电压初始值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%下面开始计算delta_P/delta_Q/delta_Ufor i=2:1:5for j=1:1:5Sum_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) );enddelta_P(i)=P(i)-e(i)*Sum_GB1-f(i)*Sum_GB2;if i~=2 && i~=3 %不为节点2,3则计算无功delta_Q(i)=Q(i)-f(i)*Sum_GB1+e(i)*Sum_GB2;endif i==2 || i==3 %这里计算delta_U的值，始终为零delta_U(i)=U(i)^2-( e(i)^2 + f(i)^2 );endSum_GB1=0;Sum_GB2=0;end%___________________________________%%下面计算雅克比矩阵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) );endfor jj=2:1:5j=jj-1;if ii~=2 && ii~=3 %PQ节点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)=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)=(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));endelse%PV节点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,2*j-1)=0;J(2*i,2*j)=0;endendendSum_GB1=0;Sum_GB2=0;end%在求解修正方程之前建议把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)];%下面求解修正方程;注意矩阵运算时候的左除和右除的区别delta_ef=-J\delta_PQV;%下面修正各个节点的电压for i=2:1:5e(i)=e(i)+delta_ef(2*(i-1)-1);f(i)=f(i)+delta_ef(2*(i-1));end%到这里第一轮迭代完成end%电压幅值和相角U=e+f*1i;angle_U=angle(U)*180/pi;%节点注入无功，流入为正，流出为负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.045^2 * 0.067;Qc3=Q(3)+0.19-1.01^2 * 0.022;U=abs(U);disp(['Iteration times : ' num2str(cont)]);%显示最终的迭代次数牛顿算法求解潮流 (极坐标):clear;clc;%牛顿算法求解潮流 (极坐标)%计算网络的潮流分布%其中节点5是平衡节点%节点1、2、3是PQ节点，节点4是PV节点% 如果节点有对地导纳支路%需将对地导纳支路算到自导纳里面%------------------------------------------------%%输入原始数据，每条支路的导纳数值，包括自导和互导纳；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.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];%导纳矩阵的实部和虚部G = real(Y);B = imag(Y);%给点电压初始值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))+U_abs(i)*sin(angle_U(i))*1i;%end%原始节点功率%这里电源功率为正，负荷功率为负%下面给点PQ PV节点功率值S=[-0.22-0.14*1i,-0.18-0.09*1i,-0.27-0.13*1i,0.35,0];%节点功率的P QP = real(S);Q = imag(S);%下面是PV节点的无功初始值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;Sum_GB2=0;cont=0;%最外层循环，cont代表迭代的次数，这里可以用约束条件来代替%for cont=1:1:4while max(delta_PQ)>1e-6,%下面计算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_U(j)) - B(i,j)*cos(angle_U(i)-angle_U(j)) );enddelta_P(i)=P(i)-U(i)*Sum_GB1;if i~=4 %不为节点四则计算无功delta_Q(i)=Q(i)-U(i)*Sum_GB2;endSum_GB1=0;Sum_GB2=0;end%_______________________________________________________%%下面计算雅克比矩阵J=zeros(7,7);for ii=1:1:4for jj=1:1:4if ii ~= 4 %PQ节点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,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)=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节点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(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全部放在一个矩阵delta_PQ=[delta_P(1);delta_Q(1);delta_P(2);delta_Q(2);delta_P(3); delta_Q(3);delta_P(4)];%下面求解修正方程;注意矩阵运算时候的左除和右除的区别J=J(1:7,1:7);delta_ef=-J\delta_PQ;%下面修正各个节点的电压for i=1:1:4if i~=4U(i)=U(i)+delta_ef(2*i)*U(i);endangle_U(i)=angle_U(i)+delta_ef(2*i-1);end%到这里第一轮迭代完成end%下面显示出满足条件后的迭代的次数disp(['Iteration times : ' num2str(cont)]);%下面计算平衡节点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/pi;End计算计算法P-Q算法计算潮流：这个算法是由牛顿算法的极坐标形式简化而来。

计算电流及功率
for n=1:nbus Sc = conj(conj(V(n))*(Y(n,n)*V(n))); DP(n) = P(n) - real(Sc); DQ(n) = Q(n) - imag(Sc); end

返回
Qi Lii Vi Vi V j (Gij sin ij Bij cos ij ) 2Vi 2 Bii Vi 2 Bii Qi ji i j i
直流潮流与牛拉法的结合
1.通过直流潮流计算各节点的相角。 2.将结果作为牛拉法初值进行计算。
J ij
Qi Vi V j (Gij cos ij Bij sin ij ) j
Qi Lij V j Vi V j (Gij sin ij Bij cos ij ) Vj
雅克比矩阵的形成
当j＝i 时
H ii
Pi Vi V j (Gij sin ij Bij cos ij ) Vi 2 Bii Qi ji i j i

数据输入：
将各节点的初始值及线路数据输入excel文 件中，并通过xlsread这一函数输入到MATLAB 中。
形成导纳矩阵
1.形成一个零矩阵Y 2.检测线路参数前两列，将Y对应元素加上相应 参数。
形成导纳矩阵


1.非对角线元素：
for n=1:nbr Y(n1(n),n2(n))= Y(n1(n),n2(n))-Y1(n)/(a(k)^2); Y(n2(n),n1(n))= Y(n1(n),n2(n)); end
总体思路：
1.直接使用“牛顿——拉夫逊”法，通过图像 显示每次迭代后的误差并记录程序运行时间。 2.直接使用“PQ分解法”，通过图像显示每次 迭代后的误差并记录程序运行时间。 3.首先使用直流潮流计算，计算各节点电压的 相角，并将计算结果作为 “牛顿——拉夫逊” 法的初值，通过图像显示每次迭代后的误差并 记录程序运行时间。

y34 4 y440
7
将接在同一节点上的接 地导纳进行并联，得
I1
y12 y10 y13 y23 y20
I2
y10 = y120 + y130
⎫ ⎪ y20 = y210 + y230 ⎪ ⎬ y30 = y310 + y320 + y340 ⎪ ⎪ y40 = y440 + y430 ⎭
3 y30 y34 4 y40 I4
& ⎤ L Y1i L Y1 j L Y1n ⎤ ⎡U 1 ⎢ ⎥ ⎥ O M O M O M ⎥⎢ M ⎥ & ⎥ L Yii L Yij L Yin ⎥ ⎢U i ⎥⎢ ⎥ O M O M O M ⎥⎢ M ⎥ & ⎥ L Y ji L Y jj L Y jn ⎥ ⎢U ⎥⎢ j⎥ O M O M O M ⎥⎢ M ⎥ ⎢U ⎥ & L Yni L Ynj L Ynn ⎥ ⎥ ⎦⎢ ⎣ n⎥ ⎦
&⎤ L Z1i L Z1 j L Z1n ⎤ ⎡ I 1 ⎢ ⎥ O M O M O M ⎥ ⎥⎢ M ⎥ &⎥ L Z ii L Z ij L Z in ⎥ ⎢ I i ⎥⎢ ⎥ O M O M O M ⎥⎢ M ⎥ &⎥ L Z ji L Z jj L Z jn ⎥ ⎢ I ⎥⎢ j⎥ O M O M O M ⎥⎢ M ⎥ ⎢I &⎥ L Z ni L Z nj L Z nn ⎥ ⎥ ⎢ ⎦⎣ n⎥ ⎦
22
节点阻抗矩阵
Z =Y
−1
节点阻抗矩阵
Z ii
Z ij
对角元素,称为节点i的自阻抗 非对角元素称为节点i和节点j之间的互阻抗
⎡ Z11 ⎢ M ⎢ ⎢ Z i1 ⎢ ⎢ M ⎢ Z j1 ⎢ ⎢ M ⎢Z ⎣ n1

一、潮流计算的计算机方法对于复杂网络的潮流计算，一般必须借助电子计算机进行。

其计算步骤是：建立电力网络的数学模型，确定计算方法、制定框图和编制程序。

本章重点介绍前两部分，并着重阐述在电力系统潮流实际计算中常用的、基本的方法。

1，电力网络的数学模型电力网络的数学模型指的是将网络有关参数相变量及其相互关系归纳起来所组成的．可以反映网络性能的数学方程式组。

也就是对电力系统的运行状态、变量和网络参数之间相互关系的—种数学描述。

电力网络的数学模型有节点电压方程和回路电流方程等，前者在电力系统潮流计算中广泛采用。

节点电压方程又分为以节点导纳矩阵表示的节点电压方程和以节点阻抗矩阵表示的节点电压方程。

（1）节点导纳矩阵在电路理论课中。

已讲过了用节点导纳矩阵表示的节点电压方程：对于n个节点的网络其展开为：上式中，I是节点注入电流的列向量。

在电力系统计算中，节点注入电流可理解为节点电源电流与负荷电流之和，并规定电源向网络节点的注人电流为正。

那么，只有负荷节点的注入电流为负，而仅起联络作用的联络节点的注入电流为零。

U是节点电压的列向量。

网络中有接地支路时，通常以大地作参考点，节点电压就是各节点的对地电压。

并规定地节点的编号为0。

y是一个n×n阶节点导纳矩阵，其阶数n就等于网络中除参考节点外的节点数。

物理意义：节点i单位电压，其余节点接地，此时各节点向网络注入的电流就是节点i 的自导纳和其余节点的与节点i之间的互导纳。

特点：对称矩阵，稀疏矩阵，对角占优(2) 节点阻抗矩阵对导纳阵求逆，得：其中称为节点阻抗矩阵，是节点导纳矩阵的逆阵。

物理意义：节点i注入单位电流，其余节点不注入电流，此时各节点的电压就是节点i 的自阻抗和其余节点的与节点i之间的互阻抗。

特点：满阵，对称，对角占优2，功率方程、变量和节点分类（1）功率方程已知的是节点的注入功率，因此，需要重新列写方程： **==B B B B B U S I U Y其展开式为： i i i nj j ij U jQ P U Y ~1-=∑= 所以：∑=**=+nj jij i i i U Y U jQ P 1 展开写成极坐标方程的形式：)cos sin ()sin cos (11ij ij ij ij n j j i i ij ij ij ij n j j i i B G U U Q B G U U P δδδδ-=+=∑∑==所以节点的功率方程为：)cos sin ()sin cos (11ij ij ij ij n j j i di Gi i ij ij ij ij nj j i di Gi i B G U U Q Q Q B G U U P P P δδδδ---=∆+--=∆∑∑==（2） 变量分类负荷消耗的有功、无功功率取决于用户，因而是无法控制的，故称为不可控变量或扰动变量。

第四章潮流计算的计算机算法第一节概述潮流计算是电力系统最基本、最常用的计算。

根据系统给定的运行条件、网络接线及元件参数，通过潮流计算可以确定各母线的电压（幅值及相角），各元件中流过的功率、整个系统的功率损耗等。

潮流计算是实现电力系统安全经济发供电的必要手段和重要工作环节。

因此潮流计算在电力系统的规划设计、生产运行、调度管理及科学研究中都有着广泛的应用。

电力系统潮流计算分为离线潮流计算和在线潮流计算。

前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。

本章主要讨论离线潮流计算问题，它的基本算法同样适用于在线潮流计算。

潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种。

自从五十年代计算机应用于电力系统以来，当时求解潮流的方法是以节点导纳矩阵为基础的逐次代入法（导纳法），后来为解决导纳法的收敛性较差的问题，出现了以阻抗矩阵为基础的逐次代入法（阻抗法）。

到六十年代，针对阻抗法占用计算机内存大的问题又出现了分块阻抗法及牛顿－拉夫逊（Newton－Raphson）法。

Newton —Raphson法是数学上解非线形方程式的有效方法，有较好的收敛性。

将N－R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧，使N－R法在收敛性、占用内存、计算速度方面的优点都超过了阻抗法，成为六十年代末期以后普遍采用的方法。

同时国内外广泛研究了诸如非线形规划法、直流法、交流法等各种不同的潮流计算方法。

七十年代以来，又涌现出了更新的潮流计算方法。

其中有1974年由B、Stott、O、Alsac 提出的快速分解法以及1978年由岩本伸一等提出的保留非线性的高129速潮流计算法。

其中快速分解法（Fast decoupled load flow）从1975年开始已在国内使用，并习惯称之为PQ分解法。

由于PQ分解法在计算速度上大大超过N－R法，不但能应用于离线潮流计算，而且也能应用于在线潮流计算。

潮流计算的主要方法
最近几年，随着计算机仿真技术和复杂系统全面发展，潮流计算也受到越来越多的重视。

潮流计算是研究不同电力网络的物理特性和操作规律的一项重要工作。

针对潮流计算的主要方法，总结如下：
一、基于动力学的方法
1. 碰撞模型：根据动力学方法，计算电力系统的运行稳定性。

基于动力学的碰撞模型能够快速而精确地预测两个潮流的变化情况。

2. 时变快速收敛：在碰撞模型的基础上，为快速求解电力系统潮流，提出了时变快速收敛算法。

可以更快地获得潮流解。

二、基于牛顿迭代法的方法
1.牛顿迭代潮流计算方法：根据牛顿迭代法，采用迭代算法，求解电力系统潮流运行状态。

2. 功率流计算方法：计算机基于牛顿迭代法，快速求解节点电能的功率流公式。

可以有效的缩短潮流计算的时间，提高计算效率。

三、基于模糊聚类算法的方法
1. 基于模糊聚类的潮流计算方法：采用模糊聚类算法，对潮流计算进行多维度分析，可以得出最优的潮流结果。

2. 基于模糊划分的多目标模糊控制：根据模糊聚类理论，对潮流算法进行最佳控制，以满足电力网不同优化目标。

四、基于期望最大化的方法
1、基于粒子群优化的潮流计算方法：采用粒子群优化算法，将电力网潮流计算定义为多目标最优化问题，以期望最大化来求解潮流值，提高计算效率。

2、基于遗传算法的潮流计算方法：遗传算法利用进化过程来搜索全局最优解，使用遗传变异原则来改变候选解，以期望最大化来求解潮流计算问题。

高等电力系统分析（潮流计算的计算机算法）PQ分解法潮流计算（IEEE14）目录一、MATLAB源程序二、对支路参数（B1）、节点参数（B2）的说明三、带入数据，运行结果非对角元 对角元MATLAB 源程序 clearclose alln=in put(' 请输入节点数:n=');n 仁in put(' 请输入支路数：n 1=') )isb=in put(' 请输入平衡节点号 :isb=')pr=i nput(' 请输入误差精度:pr=');B1=i nput(' 请输入支路参数：B1 =');B2=i nput(' 请输入节点参数:B2=');n2=in put(' 请输入PQ 节点个数： :n2=');Y=zeros (n); for i=1: n1P=B1(i,1)；q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4)*1j); % Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; %Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; enddisp(' 导纳矩阵Y=');disp(Y)% --------------------------------------------% -------------- 下面是求P,Q,V,0矩阵 -----------------------V=zeros(1, n);O=zeros(1, n);P=zeros(1, n);Q=zeros(1, n); G=real(Y);B=imag(Y); for i=1: nP(i)=B2(i,3);Q(i)=B2(i,4);V(i)=B2(i,5);O(i)=B2(i,6);endB3=B(1: n-1,1: n-1); % 不含平衡节点，由节点导纳虚部构成B4=B(1:n2,1:n2); % 所有 PQ 节点% --------------------------------------------% -------------- 下面是求△ P, △ Q 矩阵 ---------------------DX=0;ICT=1;Mp=1;Mq=1;while ICT~=0m1=1;m2=1;for i=1:n所有节点数 节点数 求 DP,DQ 对V 矩阵求逆 △ P/V △ P/V/B3 △角=-△ P/V/V/B3△ V=- △ Q/V/Bif i~=isbC(i)=O;D(i)=O;for j1=1:n C(i)=C(i)+V(i)*V(j1)*(G(i,j1)*cos(O(i)-O(j1))+B(i,j1)*si n( O(i)-O(j1)))；D(i)=D(i)+V(i)*V(j1)*(G(i,j1)*si n(O(i)-O(j1))-B(i,j1)*cos(O(i)-O( j1)))；endDP(m1)=P(i)-C(i); m1=m1+1;if B2(i,2)==1 DQ(m2)=Q(i)-D(i); m2=m2+1;endendendm1=m1-1;% m2=m2-1; %PQ DPQ=[DP';DQ']; % V 仁 V(:,1:m1);V2=diag(V1);V3=i nv(V2);% H=V3*DP'; %K=-i nv(B3)*H; %- deltO=V3*K;% max1=max(abs(DP));for i=1:m1if max1<prMp=0;elseO(i)=O(i)+deltO(i)'; Mq=1;endendV4=V(:,1:m2);V5=diag(V4);V6=i nv(V5);L=V6*DQ';N=-i nv(B4)*L; deltV=N; %max2=max(abs(DQ));for i=1:m2if max2<prMq=0;elseif B2(i,2)==1;V(i)=V(i)+deltV(i)：Mp=1;endendendif Mp==0&&Mq==0ICT=0;elseICT=1;endDX=DX+1;end% ------------------------------------------% --------------- 迭代结束，开始输出结果disp(' --------------------------------------------');disp(' 迭代次数为：');disp(DX);for i=1: nE(i)=V(i)*cos(O(i))+1j*V(i)*si n(O(i));o(i)= 180*a ngle(E(i))/pi;enddisp(' --------------------------------------------');disp(' 修正后各节点电压标么值为(节点号从小到大排列) :');disp(V);disp(' --------------------------------------------');disp(' 修正后各节点电压相角为(节点号从小到大排列) :');disp(o);% ----------- 计算各个节点的功率disp(' --------------------------------------------');disp(' 各节点的功率为：')；for p=1: nC(p)=0;for q=1: nC(p)=C(p)+conj(Y(p,q)*conj(E(q)));endS(p)=E(p)*C(p);enddisp(S);% ----------- 计算各支路的功率 -------------------------for i=1: n1 p=B1(i,1);q=B1(i,2);Si(p,q)=E(p)*(conj(E(p))*conj(Y(p,p)-Y(p,q))+(conj(E(p))-conj(E(q )))*conj(Y(p,q))); disp(' -------------------------------------------- ');disp(' 各条支路的首端功率为:');disp(Si(p,q));Si(q,p)=E(q)*(conj(E(q))*conj(Y(q,q)-Y(p,q))+(conj(E(q))-conj(E(p )))*conj(Y(p,q))); disp(' -------------------------------------------- ');disp(' 各条支路的末端功率为:');disp(Si(q,p));DS(i)=Si(p,q)+Si(q,p);disp(' -------------------------------------------- ');disp(' 各条支路的功率损耗为:');disp(DS(i));end% --------------- 计算平衡节点功率-----------------------Sp=0;for i=1: nSp=Sp+V( n)*conj(Y( n, i))*conj(V(i));enddisp(' -------------------------------------------- ');disp(' 平衡节点功率为：')；disp(Sp);二、对支路参数（B1）、节点参数（B2）的说明1•节点数：142•支路数：203•支路矩阵B1的各支路参数：起点编号，终点编号，电阻，电抗，电导，电纳[1 2 0.01335 0.04211 0 0 ；1 3 0 0.20912 0 0 J1 4 0 0.55618 0 0 J1 10 0.05811 0.17632 0 0.0341 11 0.06701 0.17103 0 0.01282 10 0.05695 0.17388 0 0.0346 2 12 0 0.25202 0 0 J2 14 0.05403 0.22304 0 0.04923 4 0 0.11001 0 0 J3 13 0 0.17615 0 0 J4 5 0.03181 0.0845 0 0 ；4 9 0.12711 0.27038 0 0 ；5 6 0.08205 0.19207 0 0 ；；6 12 0.09498 0.1989 0 07 8 0.22092 0.19988 0 0 ;[1 1 -0.4782 1 -0.0763 10 04 1 -0.2955 1 -0.096 1 -0.0357 1 -0.0618 1 -0.1359 1 -0.14910 2 0.18311 2 -0.94212 2 -0.11213 214 00.039 1 0;-0.016 1 0;1 0;-0.166 1 0;-0.058 1 0;-0.018 1 0;-0.016 1 0;-0.058 1 0;-0.05 1 0;0 1.045 0;0 1.01 0;0.047 1.70;0 0.174 1.9 0;0 0 1.06 0;]7 12 0.12291 0.25581 0 08 9 0.17093 0.34802 0 0 ;8 12 0.06615 0.13027 0 0 ;10 11 0.04699 0.19797 0 0.043810 14 0.01938 0.05917 0 0.05284.节点参数矩阵B2的各节点参数：（对应的每一列为）节点编号，类型，注入有功，注入无功，电压幅值，电压相位其中节点类型：仁PQ节点，2=PV节点，0二平衡节点三、带入数据，运行结果>> clearclose alln=input('请输入节点数：n=');n1=input('请输入支路数：n仁');isb=i nput('请输入平衡节点号：isb='); pr=input('请输入误差精度:pr=');B1= input('请输入支路参数:B仁')；B2=input('请输入节点参数:B2=');n2=input('请输入PQ节点个数:n2=');Y=zeros (n);for i=1: n1p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4)*1j); %非对角元Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; % 对角元Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j;enddisp('导纳矩阵Y=');disp(Y)% --------------------------------------------% -------------- 下面是求PQ,V,0矩阵-------------V=zeros(1, n);O=zeros(1, n);P=zeros(1, n);Q=zeros(1, n); G=real(Y);B=imag(Y); for i=1: nP(i)=B2(i,3);Q(i)=B2(i,4);V(i)=B2(i,5);O(i)=B2(i,6);endB3=B(1: n-1,1: n-1); %不含平衡节点，由节点导纳虚部构成B4=B(1:n2,1:n2); %所有PQ 节点% --------------------------------------------endendend m1=m1-1;m2=m2-1; DPQ=[DP';DQ']; V1=V(:,1:m1); V2=diag(V1);V3=i nv(V2);H=V3*DP';K=-i nv(B3)*H; deltO=V3*K;% 所有节点数%PQ节点数%求DPQQ%对V矩阵求逆%A P/V%-A P/V/B3%A 角=-A% -------------- 下面是求△ P,A Q矩阵------------DX=0;ICT=1;Mp=1;Mq=1;while ICT~=0m1=1;m2=1;for i=1: nif i~=isbC(i)=0;D(i)=0;for j1=1: nC(i)=C(i)+V(i)*V(j1)*(G(i,j1)*cos(O(i)-O(j1))+B(i,j1)*si n(O(i)-O(j1)));D(i)=D(i)+V(i)*V(j1)*(G(i,j1)*si n(O(i)-O(j1))-B(i,j1)*cos(O(i)-O(j1)));endDP(m1)=P(i)-C(i); m1=m1+1;if B2(i,2)==1DQ(m2)=Q(i)-D(i); m2=m2+1;max1=max(abs(DP));for i=1:m1if max1<prMp=0;elseO(i)=O(i)+deltO(i)';Mq=1;endendV4=V(:,1:m2);V5=diag(V4); V6=i nv(V5);L=V6*DQ';N=-i nv(B4)*L;deltV=N; %△V=-A Q/V/Bmax2=max(abs(DQ));for i=1:m2if max2<prMq=0;elseif B2(i,2)==1;V(i)=V(i)+deltV(i)： Mp=1;endendendif Mp==0&&M q==0ICT=0;elseICT=1;endDX=DX+1;end% ------------------------------------------% --------------- 迭代结束，开始输出结果disp(' --------------------------------------------');disp('迭代次数为：');disp(DX);for i=1: nE(i)=V(i)*cos(O(i))+1j*V(i)*si n(0 (i));o(i)= 180*a ngle(E(i))/pi;enddisp(' --------------------------------------------');disp('修正后各节点电压标么值为(节点号从小到大排列) :');disp(V);disp(' --------------------------------------------');disp('修正后各节点电压相角为(节点号从小到大排列) :');disp(o);% ----------- 计算各个节点的功率-----------------disp(' --------------------------------------------');disp('各节点的功率为:');for p=1: nC(p)=0;for q=1: nC(p)=C(p)+conj(Y(p,q)*conj(E(q)));endS(p)=E(p)*C(p);end disp(S);% ----------- 计算各支路的功率------------------for i=1: n1 p=B1(i,1);q=B1(i,2);Si(p,q)=E(p)*(conj(E(p))*conj(Y(p,p)-Y(p,q))+(conj(E(p))-conj(E(q)))*conj(Y(p,q)));disp(' -------------------------------------------- ');disp('各条支路的首端功率为:');disp(Si(p,q));Si(q,p)=E(q)*(conj(E(q))*conj(Y(q,q)-Y(p,q))+(conj(E(q))-conj(E(p)))*conj(Y(p,q))); disp(' -------------------------------------------- ');disp('各条支路的末端功率为:');disp(Si(q,p));DS(i)=Si(p,q)+Si(q,p);disp(' -------------------------------------------- ');disp('各条支路的功率损耗为:');disp(DS(i));end% --------------- 计算平衡节点功率------------Sp=0;for i=1: nSp=Sp+V( n)*conj(Y( n,i))*conj(V(i));enddisp(' -------------------------------------------- ');disp('平衡节点功率为:');disp(Sp);请输入节点数：n=14请输入支路数：n 1=20请输入平衡节点号:isb=14请输入误差精度：pr=0.00001请输入支路参数:B1=[1 2 0.01335 0.04211 0 0 ;1 3 0 0.20912 0 0 J1 4 0 0.55618 0 0 J1 10 0.05811 0.17632 0 0.0341 11 0.06701 0.17103 0 0.01282 10 0.05695 0.17388 0 0.03462 12 0 0.25202 0 0 J2 14 0.05403 0.22304 0 0.049210.5130 -38.2963i -6.8410 +21.5786i 0.0000 + 4.7819i0.0000 + 1.7980i0.0000 +0.0000i-6.8410 +21.5786i 9.5680 -34.8916i0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 +0.0000i0.0000 + 4.7819i 0.0000 + 0.0000i0.0000 -19.5490i 0.0000 + 9.0901i0.0000 +0.0000i0.0000 + 1.7980i0.0000 + 0.0000i0.0000 + 9.0901i5.3261 -24.2825i-3.9020 +10.3654i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-3.9020 +10.3654i 5.7829-14.7683i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i 0.0000 + 0.0000i-1.8809 +4.4029i0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 +O.OOOOi3 4 0 0.11001 00 ; 3 13 0 0.17615 00 J 4 5 0.03181 0.08450 ; 4 9 0.12711 0.27038 0 0 ; 5 60.08205 0.19207 00 ; 6 12 0.09498 0.1989 0 0 ; 7 8 0.22092 0.19988 0 0 ; 7 12 0.12291 0.25581 00 ; 8 9 0.17093 0.34802 0 0 ; 8 12 0.06615 0.13027 0 0 ;10 11 0.04699 0.19797 00.0438 ; 10 14 0.01938 0.05917 0 0.0528;请输入节点参数：B2=[1 1 -0.478 0.039 2 1 -0.076 -0.016 10;3 1 0 0 1 0;4 1 -0.295 -0.166 1 0; 5 1 -0.09 -0.058 1 0; 6 1 -0.035 -0.018 1 0;7 1 -0.061 -0.016 1 0;8 1 -0.135 -0.058 1 0;9 1 -0.149 -0.05 10; 10 2 0.183 01.0450;11 2 -0.9420 1.01 0;12 2 -0.112 0.047 1.7 0;13 20 0.1741.9 0;14 0 0 0 1.06 0;]请输入PQ 节点个数:n2=9导纳矩阵 Y=1Columns 1 through 5；0.0000 +O.OOOOiO.OOOOi0.0000 + 0.0000i 0.0000i-1.6860 + 5.1158i 0.0000i-1.9860 + 5.0688i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.7011 + 5.1939i0.0000 + 0.0000i0.0000 + 3.9679i0.0000 + 0.0000i-1.0259 + 4.2350i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 5.6770i0.0000 + 0.0000i0.0000 + 0.0000i-1.4240 + 3.0291i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000iColu mns 6 through 100.0000 + 0.0000i 5.115 8i0.0000 + 0.0000i 5.1939i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i-1.8809 + 4.4029i 0.0000i3.8359 - 8.4970i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i -30.1895i0.0000 + 0.0000i 4.7819i-1.9550 + 4.0941i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i4.0150 -5.4279i-2.4890 + 2.2520i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.5260 + 3.1760i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-2.4890 + 2.2520i6.7249 -10.6697i-1.1370 + 2.3150i0.0000 + 0.0000i0.0000 + 0.0000i-3.0989 + 6.1028i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.4240 + 3.0291i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.1370 + 2.3150i2.5610 - 5.3440i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 +0.0000 +0.0000 +0.0000 +0.0000 +0.0000 +0.0000 +-1.6860 +-1.7011 +0.0000 +0.0000 +0.0000 +0.0000 +0.0000 +0.0000 +0.0000 +9.5213-1.1350 +0.0000 +0.0000 +0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -4.9991-1.9860 + 5.0688i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -1.1350 + 4.7819i 3.1210 - 9.7941i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 3.9679i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.9550 + 4.0941i-1.5260 + 3.1760i-3.0989 + 6.1028i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i6.5799 -17.3407i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 5.6770i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 - 5.6770i0.0000 + 0.0000i0.0000 + 0.0000i-1.0259 + 4.2350i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-4.9991 +15.2631i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i6.0250 -19.3961i+15.2631iColu mns 11 through 14迭代次数为:42修正后各节点电压标么值为Colu mns 1 through 10（节点号从小到大排列）1.2128 1.21481.5809 1.04501.5627 1.5364 1.5602 1.6264 1.6792 1.6654Colu mns 11 through 141.0100 1.7000 1.90001.0600修正后各节点电压相角为（节点号从小到大排列）Colu mns 1 through 10-13.1377 -11.9998 -14.4540 -15.0006 -15.3011 -15.7590 -16.4954 -16.3362 -16.0220 -5.2081Colu mns 11 through 14-12.1568 -16.1918 -14.4540-0.0475 --1.5950 -各节点的功率为：Columns 1 through 5 1.4718 - 1.8800i 0.9041 - 1.1423i0.0000 - O.OOOOi -0.1314 + 0.1399i0.0646iColu mns 6 through 10 -0.0300 + 0.0035i 0.0039 - 0.1228i -0.0656 - 0.0724i-0.0352 - 0.1437i1.3905iColu mns 11 through 14 -0.4870 - 1.2396i3.2005 + 3.0414i1.7586 + 3.1847i -2.7057 + 0.3992i各条支路的首端功率为26.1711 +87.9128i各条支路的末端功率为23.5667 +83.4812i各条支路的功率损耗为：4.9738e+01 + 1.7139e+02i各条支路的首端功率为15.2545 +65.3872i各条支路的末端功率为0.2082 +56.8008i各条支路的功率损耗为：1.5463e+01 + 1.2219e+02i各条支路的首端功率为15.3538 +59.6752i各条支路的末端功率为12.6813 +60.6685i各条支路的功率损耗为：2.8035e+01 + 1.2034e+02i各条支路的首端功率为18.4735 +62.4536i各条支路的末端功率为11.6194 +39.6840i各条支路的功率损耗为：3.0093e+01 + 1.0214e+02i各条支路的首端功率为18.0012 +62.4930i各条支路的末端功率为5.5097 +16.2405i各条支路的功率损耗为23.5109 +78.7335i各条支路的首端功率为17.0452 +57.7870i各条支路的末端功率为11.7622 +39.7706i各条支路的功率损耗为28.8075 +97.5576i各条支路的首端功率为13.5220 +59.6677i19.6150 +58.2875i各条支路的功率损耗为：3.3137e+01 + 1.1796e+02i各条支路的首端功率为16.5471 +56.5546i各条支路的末端功率为6.9281 +27.4025i各条支路的功率损耗为23.4752 +83.9571i各条支路的首端功率为-0.2082 +69.5644i各条支路的末端功率为12.7806 +79.1441i各条支路的功率损耗为：1.2572e+01 + 1.4871e+02i各条支路的首端功率为0.0000 +64.5962i各条支路的末端功率为-0.0000 +37.3498i各条支路的功率损耗为：-1.7764e-15 + 1.0195e+02i21.7957 +82.2160i各条支路的末端功率为23.5612 +60.7482i各条支路的功率损耗为：4.5357e+01 + 1.4296e+02i各条支路的首端功率为15.8994 +64.7377i各条支路的末端功率为9.9896 +20.6498i各条支路的功率损耗为25.8891 +85.3875i各条支路的首端功率为18.7607 +47.1609i各条支路的末端功率为15.0088 +33.6107i各条支路的功率损耗为33.7695 +80.7716i各条支路的首端功率为15.4666 +33.8366i各条支路的末端功率为24.5068 +61.3933i39.9734 +95.2299i各条支路的首端功率为18.2997 +21.5842i各条支路的末端功率为25.5956 +35.9107i各条支路的功率损耗为43.8954 +57.4949i各条支路的首端功率为15.7256 +24.3488i各条支路的末端功率为23.3240 +59.2040i各条支路的功率损耗为39.0496 +83.5528i各条支路的首端功率为21.6792 +35.6719i各条支路的末端功率为9.3603 +19.4667i各条支路的功率损耗为31.0395 +55.1385i各条支路的首端功率为27.4696 +46.8495i各条支路的末端功率为27.7461 +67.4149i各条支路的功率损耗为：5.5216e+01 + 1.1426e+02i各条支路的首端功率为10.9761 +38.1226i各条支路的末端功率为4.9835 +14.8560i各条支路的功率损耗为15.9596 +52.9785i各条支路的首端功率为17.4469 +49.3022i各条支路的末端功率为10.7497 +39.1332i各条支路的功率损耗为28.1966 +88.4354i平衡节点功率为：-0.0889 - 0.5670i。

潮流计算机算法熊飞
1、运用计算机计算时，一般要完成一下几个步骤：
（1）建立数学模型。

（2）确定解算方法。

（3）制定程序框图。

（4）编制程序。

（5）上机调试及运算。

2、潮流计算的数学模型
数学模型是指反映电力系统中运行状态参数【如电压、电力、功率等】与网络参数之前的关系，反映网络性能的数学方程式。

3、牛顿-拉夫逊法
（1）牛顿-拉夫逊法不仅在多数情况下没有发散的危险，而且收敛性较强，可以大大节约计算时间，因而得到了广泛的应用。

它最大的特点是初始值的选择要求严格，必须选好恰当的初始值，否则不收敛。

（2）该方法的要点就是把对非线性方程的求解过程转化为反复对相应的线性方程求解的过程，通常称为逐次线性化过程，这是牛顿-拉夫逊法的核心。

（3）牛顿-拉夫逊法潮流计算的求解过程
1、输入原始数据和信息。

2、形成节点导纳矩阵Y。

3、送电压初始值e, f。

4、求不平衡量。

5、计算雅可比矩阵的各个元素。

6、解修正方程。

7、求节点电压新值。

8、判断是否收敛。

9、反复迭代4~7，直至满足第8步的条件。

10、求平衡节点的功率和PV节点的无功功率e及各支路的功率。

附表1：计算机计算潮流程序：%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号clear;n=13;%input('请输入节点数：n=');nl=13;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[];%input('请输入由支路参数形成的矩阵：B1=');B2=[];%input('请输入各节点参数形成的矩阵：B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); %－－－－－－－修改部分－－－－－－－－－－－－ym=0;SB=100;UB=220;%ym=input('您输入的参数是标么值？（若不是则输入一个不为零的数值）'); if ym~=0%SB=input('请输入功率基准值:SB=');%UB=input('请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1矩阵B1=');disp(B1)for i=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2矩阵B2=');disp(B2)end% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于低压侧p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(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; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵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节点电压给定模值endfor i=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%=========================================================== ========P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1:nif i~=isb %非平衡节点C(i)=0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求P',Q'V2=e(i)^2+f(i)^2; %电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素=========if B2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算=================%================= 求取Jacobi矩阵===================for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df=dQ/deX3=X2; % X2=dp/df X3=dQ/deX4=-X1; % X1=dP/de X4=dQ/dfp=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);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素===========DP=P(i)-P1; % PV节点有功误差DV=V(i)^2-V2; % PV节点电压误差for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/dfX5=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);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX5=-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;endendendendend%========= 以上为求雅可比矩阵的各个元素===================== for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%======================================================== ====k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行前的k4行消去for k2=k1:N1 % k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endelsefor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endendend%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵=====for k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); %修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%=============== 计算各输出量=========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);VV=V*UB;disp(VV);disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率S为(节点号从小到大排列)：');disp(S);disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');SS=S*SB;disp(SS);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./( B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);elseSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))-conj(E(q)))*conj(1./( B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);enddisp(Si(p,q));SSi(p,q)=Si(p,q)*SB;ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);%disp(SSi(p,q));disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./( B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);elseSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))-conj(E(p)))*conj(1./( B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p)*SB;ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);%disp(SSj(q,p));disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);%disp(DDS(i));disp('-----------------------------------------------------');endfigure(1);subplot(2,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(2,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;subplot(2,2,3);bar(real(S));ylabel('节点注入有功');grid on;subplot(2,2,4);bar(Siz);ylabel('支路首端无功');grid on;1.冬季最大运行方式潮流计算结果:计算机运行的B1,B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.88+0.545*i 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.88+0.545*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.2+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：节点号电压值相角值支路标号首端功率末端功率支路功率损耗1 242.0000 0 1-2 148.9391+3.327637i -143-27.45476i 5.93913-24.1271i2 231.0000 -3.0341 1-11 89.1956+37.193i -88.19803-61.4687i 0.997519-24.2757i3 227.4154 -4.4383 11-3 88.19803+61.4687i -88-54.5i 0.19803+6.9687i4 226.3304 -4.9004 1-12 77.8364+36.4186i -77.2404-55.7062i 0.596039-19.2876i5 229.8318 -3.9547 12-4 77.2404+55.7062i -77-47.72i 0.24036+7.9862i6 217.7226 -8.2720 1-5 63.5438+12.2204i -61.8773-32.0321i 1.66656-19.8116i7 234.9326 -3.1047 5-6 77.2597+56.3501i -77-47.72i 0.25974+8.6301i8 226.0991 -6.2429 5-7 15.3825-24.3181i 15.5354+0.274108i 0.152961-24.044i9 231.0000 0.5851 7-8 88.20085+61.55009i -88-54.5i 0.20085+7.0501i10 231.0000 3.4950 1-7 40.806-6.05737i -40.0647-22.0555i 0.741264-28.1128i11 236.1931 -1.3348 7-13 -63.6716-39.7687i 64.0965+17.5832i 0.424934-22.1855i12 237.9346 -0.8892 13-9 -64.0965-17.5832i 64.2+21.0187i 0.10348+3.4355i13 238.1868 -2.2028 1-10 79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果图形：2.冬季最小运行方式潮流计算结果：计算机运行的B1B2矩阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.616+0.3817*i 1 0 0 20 0.539+0.3817*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.1+0.06197*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]电压调整措施：变电所1、4变压器变比：+2.5% 水电厂变压器变比：-2.5%5 234.2241 -2.0051 (12--4) 54.0234+42.2706i -53.9-38.17i 0.12342+4.1006i6 226.2159 -4.8577 (1--5) 34.1669+4.22856 -33.6427-28.276i 0.524165-24.0474i7 235.1128 -0.4737 (5--6) 54.0179+37.3169i -53.9-33.4i 0.11789+3.9169i8 223.9941 -2.4603 (5--7) -20.3752-9.04094i 20.5371-15.4558i 0.161947-24.4968i9 231.0000 3.4429 (1--7) 11.0013+1.45186i -10.8258-31.4513i 0.175442-29.9995i10 231.0000 3.9321 (7--8) 53.9768+36.0957i -53.9-33.4i 0.076797+2.6957i11 238.4187 -0.9561 (7--13) -63.6881+10.8115i 64.0874-32.7763i 0.39932-21.9648i12 239.1742 -0.6185 (13--9) -64.0874+32.7763i 64.2-29.0386i 0.11258+3.7377i13 234.9468 0.6942 (1--10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机运行结果的图形：3.夏季最大运行方式计算机计算结果：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.776+0.481*i 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.776+0.481*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：10 231.0000 3.4950 (7,8) 77.7585+53.6634i -77.6-48.1i0.1585+5.5634i11 237.0938 -1.1858 (7,13) -113.4984+9.476406i 114.6926-28.38885i 1.19422-18.9124i12 239.1742 -0.6185 (13,9) -114.6926+28.38885 115-18.18369i0.307384+10.2052i13 233.3912 3.2303 (1,10) -79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果如图：4.夏季最小运行方式：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.543+0.336*i 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.543+0.336*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]10 231.0000 3.9321 (7,8) 54.3735+36.2509i -54.3-33.67i 0.073528+2.5809i11 239.0099 -0.8518 (7,13) -113.4428+24.39707i 114.6754-43.79301i 1.23255-19.3959i12 239.8726 -0.5689 (13,9) -114.6754+43.79301i 115-33.01613i 0.324605+10.7769i13 235.9565 4.4510 (1,10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机计算结果如图：5.夏季故障运行状态:调压及无功补偿措施如下:变电所3的变压器变比为-2.5%,无功补偿容量为20Mvar。

潮流计算研究报告一、潮流计算概述电力工业是国民经济的基础产业,为保证电力系统安全、稳定及可靠的运行,首先需计算电力系统网络潮流。

潮流计算指根据给定的电力系统运行条件及系统接线情况确定整个电力系统各部分的运行状态:各母线电压、各元件流过的功率及系统功率损耗等等。

在电力系统规划、设计及运行方式的研究中,都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性及经济性,所以潮流计算是研究电力系统的一种很重要和很基础的计算。

在数学上潮流计算问题是一组多元非线性方程式求解问题,其解法离不开迭代。

二、潮流计算算法框图三、节点分类本次设计采用直角坐标牛顿拉夫逊法，用C语言编制潮流计算程序。

牛顿法拉夫逊：把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，通常称为逐次线性化过程。

节点分类：1）平衡节点已知节点的电压幅值和相角，求有功和无功，一般设在调频厂（一个）。

2）PQ节点已知节点的有功、无功，求电压幅值和相角，一般设在变电所母线，按给定频率发电的电厂母线（多个）。

3）PV节点已知节点的有功和电压幅值，求无功和电压相角，一般设在有无功补偿的变电所母线，无功可调的发电厂母线（少个）。

四、基本公式1、节点注入电流2、节点功率3、不平衡量4、雅可比矩阵（1）i ≠ j时（2）i=j 时5、线路功率及损耗五、IEEE30拓扑结构IEEE30节点结构六、算法输入原始数据1、支路参数（支路数L=41）始端末端回线数线路类型线路电阻线路电抗线路电纳(B/2)/ 变比2 1 1 0 0.0192 0.0575 0.0528 03 1 1 0 0.0452 0.1852 0.0408 04 2 1 0 0.057 0.1737 0.0368 04 3 1 0 0.0132 0.0379 0.0084 05 2 1 0 0.0472 0.1983 0.0418 06 2 1 0 0.0581 0.1763 0.0374 06 4 1 0 0.0119 0.0414 0.009 07 5 1 0 0.046 0.116 0.0204 07 6 1 0 0.0267 0.082 0.017 08 6 1 0 0.012 0.042 0.009 06 9 1 1 0 0.208 0 0.9786 10 1 1 0 0.556 0 0.96911 9 1 0 0 0.208 0 010 9 1 0 0 0.11 0 04 12 1 1 0 0.256 0 0.93213 12 1 0 0 0.14 0 014 12 1 0 0.1231 0.2559 0 015 12 1 0 0.0662 0.1304 0 016 12 1 0 0.0945 0.1987 0 015 14 1 0 0.221 0.1997 0 017 16 1 0 0.0524 0.1923 0 018 15 1 0 0.1073 0.2185 0 019 18 1 0 0.0639 0.1292 0 020 19 1 0 0.034 0.068 0 020 10 1 0 0.0936 0.209 0 017 10 1 0 0.0324 0.0845 0 021 10 1 0 0.0348 0.0749 0 022 10 1 0 0.0727 0.1499 0 022 21 1 0 0.0116 0.0236 0 023 15 1 0 0.1 0.202 0 024 22 1 0 0.115 0.179 0 024 23 1 0 0.132 0.27 0 025 24 1 0 0.1885 0.3292 0 026 25 1 0 0.2544 0.38 0 027 25 1 0 0.1093 0.2087 0 028 27 1 1 0 0.396 0 0.96829 27 1 0 0.2198 0.4153 0 030 27 1 0 0.3202 0.6027 0 030 29 1 0 0.2399 0.4533 0 028 8 1 0 0.0636 0.2 0.0428 028 6 1 0 0.0169 0.0599 0.013 02、节点参数（节点数N=30）编号节点类型电压实部电压虚部负荷有功负荷无功发电有功发电无功PV节点的V 并联电容1 1 1.05 0 0 0 0 0 1.05 02 3 1.045 0 0.217 0.127 0.8 0 1.045 03 2 1 0 0.024 0.012 0 0 1 04 2 1 0 0.076 0.016 0 0 1 05 3 1.01 0 0.942 0.190 0.5 0 1.01 06 2 1 0 0 0 0 0 1 07 2 1 0 0.228 0.109 0 0 1 08 3 1.01 0 0.30 0.30 0.2 0 1.01 09 2 1 0 0 0 0 0 1 010 2 1 0 0.058 0.02 0 0 1 0.1911 3 1.05 0 0 0 0.2 0 1.05 012 2 1 0 0.112 0.075 0 0 1 013 3 1.05 0 0 0 0.2 0 1.05 014 2 1 0 0.062 0.016 0 0 1 015 2 1 0 0.082 0.025 0 0 1 016 2 1 0 0.035 0.018 0 0 1 017 2 1 0 0.09 0.058 0 0 1 018 2 1 0 0.032 0.009 0 0 1 019 2 1 0 0.095 0.034 0 0 1 020 2 1 0 0.022 0.007 0 0 1 021 2 1 0 0.175 0.112 0 0 1 022 2 1 0 0 0 0 0 1 023 2 1 0 0.032 0.016 0 0 1 024 2 1 0 0.087 0.067 0 0 1 0.04325 2 1 0 0 0 0 0 1 026 2 1 0 0.035 0.023 0 0 1 027 2 1 0 0 0 0 0 1 028 2 1 0 0 0 0 0 1 029 2 1 0 0.024 0.009 0 0 1 0302 1 0 0.106 0.019 0 0 1 0七、潮流程序#include<iostream.h>#include<fstream.h>#include <iomanip.h> //使用setw函数#include<math.h>#include "stdlib.h"#include<stdio.h>#include "math.h"#include <iostream>#include <stdio.h>#include <tchar.h>#include <fstream>#include <cmath>#include <process.h>#include <math.h>#include <malloc.h>using namespace std;void main(){FILE *fp;fp=fopen("result.txt","w");ifstream fin;fin.open("node information.txt");int node_num;fin>>node_num;int n = node_num*(node_num-1);int Jocabi_rank=2*(node_num-1);double * * Jocabi=new double * [Jocabi_rank];int *II=new int[n];int *JJ=new int[n];double *WR=new double[n];double *WX=new double[n];double *WXC=new double[n];double *WNT=new double[n];for(int i=0;i<n;i++){II[i]=0;JJ[i]=0;WR[i]=0;WX[i]=0;WXC[i]=0;WNT[i]=0;}double *GF=new double[n];for(int GF_i=0;GF_i<=n;GF_i++)GF[GF_i]=0;double *BF=new double[n];for(int BF_i=0;BF_i<=n;BF_i++)BF[BF_i]=0;double * * GF1=new double * [node_num+1];for (int GF1_i=0;GF1_i<node_num;GF1_i++){GF1[GF1_i]=new double[node_num+1];for (int GF1_j=0;GF1_j<node_num;GF1_j++){GF1[GF1_i][GF1_j]=0;}} double * * BF1=new double * [node_num+1];for (int BF1_i=0;BF1_i<=node_num;BF1_i++) {BF1[BF1_i]=new double[node_num+1];for (int BF1_j=0;BF1_j<node_num;BF1_j++) {BF1[BF1_i][BF1_j]=0;}} double *GD=new double[node_num];for(int GD_i=0;GD_i<=node_num;GD_i++){GD[GD_i]=0;}double *BD=new double[node_num];for(int BD_i=0;BD_i<=node_num;BD_i++){BD[BD_i]=0;}double * BC=new double[node_num];for(int BC_i=0;BC_i<=node_num;BC_i++){BC[BC_i]=0;}for(i=0;i<n;i++)fin>>II[i]>>JJ[i]>>WR[i]>>WX[i]>>WXC[i]>>WNT[i];//节点导纳矩阵int p,q;double K;for(i=0;i<n;i++){p=0;q=0;K=0.0;if(WNT[i]==0){GF[i]=WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]);BF[i]=(-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]);}if(WNT[i]>0){K=WNT[i];p=II[i];GF[i]=(WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]))/K;BF[i]=((-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]))/K;GD[p-1]=GD[p-1]+((K-1)/K)*(WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]));BD[p-1]=BD[p-1]+((K-1)/K)*(-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]);}if(WNT[i]<0){K=fabs(WNT[i]);q=JJ[i];GF[i]=(WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]))/K;BF[i]=((-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]))/K;GD[q-1]=GD[q-1]+((1-K)/(K*K))*WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]);BD[q-1]=BD[q-1]+((1-K)/(K*K))*(-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]);}}for(int I=1;I<=node_num;I++){BC[I-1]=0.0;for(int i=0;i<n;i++){int J=0;if(II[i]==I){if(WXC[i]<1000)BC[I-1]=BC[I-1]+WXC[i];J=JJ[i];GF1[I-1][J-1]=-GF[i];BF1[I-1][J-1]=-BF[i];GF1[I-1][I-1]=GF1[I-1][I-1]+GF[i];BF1[I-1][I-1]=BF1[I-1][I-1]+BF[i];}}GF1[I-1][I-1]=GF1[I-1][I-1]+GD[I-1];BF1[I-1][I-1]=BF1[I-1][I-1]+BD[I-1]+BC[I-1];} //输出节点导纳矩阵fprintf(fp,"节点导纳矩阵是： \n");for(int m=0;m<node_num;m++){for(int n=0;n<node_num;n++){fprintf(fp,"%9.6f",GF1[m][n]);fprintf(fp,"+j");fprintf(fp,"%9.6f",BF1[m][n]);fprintf(fp," ");}fprintf(fp,"\n");}fprintf(fp,"\n");fprintf(fp,"\n");fprintf(fp,"\n");//$$$$$$$$$$$$$$$$$输出节点导纳矩阵完毕$$$$$$$$$$$$$$$$$$$$$//$$$$$$$$$$$$$$$$$$$求节点初始差值$$$$$$$$$$$$$$$$$$$$$$$$$ double *wv=new double [node_num];double *we=new double [node_num];double *wf=new double [node_num];double *wp_g=new double [node_num];double *wq_g=new double [node_num];double *wp_l=new double [node_num];double *wq_l=new double [node_num];double *wn=new double [node_num];for(int p_i=0;p_i<node_num;p_i++){wv[p_i]=0;we[p_i]=0;wf[p_i]=0;wp_g[p_i]=0;wq_g[p_i]=0;wp_l[p_i]=0;wq_l[p_i]=0;wn[p_i]=0;}double * dp=new double [node_num];for (int dp_i=0;dp_i<node_num;dp_i++)dp[dp_i]=0;double * dq=new double [node_num];for (int dq_i=0;dq_i<node_num;dq_i++)dq[dq_i]=0;double * dvv=new double [node_num];for (int dvv_i=0;dvv_i<node_num;dvv_i++)dvv[dvv_i]=0;double * sm=new double [node_num];double * sn=new double [node_num];int n_m=0 ;ifstream ffin;ffin.open("node power.txt");for(int pi=0;pi<node_num;pi++){ffin>>wp_g[pi]>>wq_g[pi]>>wp_l[pi]>>wq_l[pi]>>we[pi]>>wf[pi]>> wn[pi];}for(int wv_i=0;wv_i<node_num;wv_i++){wv[wv_i]=we[wv_i];}double * w_p=new double [node_num];double * w_q=new double [node_num];for (int w_i=0;w_i<node_num;w_i++){w_p[w_i]=0;w_q[w_i]=0;}//开始迭代，次数为iterative_frequencyfor (intiterative_frequency=0;iterative_frequency<12;iterative_frequency++) {fprintf(fp,"第[");fprintf(fp,"%2d",iterative_frequency);fprintf(fp,"]次迭代结果如下：");fprintf(fp,"\n");n_m=0;for(int Jocabi_i=0;Jocabi_i<Jocabi_rank;Jocabi_i++){Jocabi[Jocabi_i]=new double[Jocabi_rank];for(int Jocabi_j=0;Jocabi_j<Jocabi_rank;Jocabi_j++)Jocabi[Jocabi_i][Jocabi_j]=0;}for(i=0;i<node_num;i++){sm[i]=0;sn[i]=0;}for(i=0;i<node_num;i++){if(wn[i]==1){for(int j=0;j<node_num;j++){sm[i]=sm[i]+(GF1[i][j]*we[j]-BF1[i][j]*wf[j]);sn[i]=sn[i]+(GF1[i][j]*wf[j]+BF1[i][j]*we[j]);}dp[i]=wp_l[i]-we[i]*sm[i]-wf[i]*sn[i];dq[i]=wq_l[i]-wf[i]*sm[i]+we[i]*sn[i];n_m++;}if(wn[i]==2){for(int j=0;j<node_num;j++){sm[i]=sm[i]+(GF1[i][j]*we[j]-BF1[i][j]*wf[j]);sn[i]=sn[i]+(GF1[i][j]*wf[j]+BF1[i][j]*we[j]);}dp[i]=wp_g[i]-we[i]*sm[i]-wf[i]*sn[i];dvv[i]=wv[i]*wv[i]-(we[i]*we[i]+wf[i]*wf[i]);w_p[i]=we[i]*sm[i]+wf[i]*sn[i];w_q[i]=wf[i]*sm[i]-we[i]*sn[i];fprintf(fp,"第[" );fprintf(fp,"%2d",i+1 );fprintf(fp, "]个节点是PV节点，且其节点功率为：");fprintf(fp,"%9.6f",w_p[i]);fprintf(fp,"+j" );fprintf(fp,"%9.6f",w_q[i]);fprintf(fp,"\n" );} if(wn[i]==0){for(int j=0;j<node_num;j++){sm[i]=sm[i]+(we[j]*GF1[i][j]-wf[j]*BF1[i][j]);sn[i]=sn[i]+(wf[j]*GF1[i][j]+we[j]*BF1[i][j]);}w_p[i]=we[i]*sm[i]+wf[i]*sn[i];;w_q[i]=wf[i]*sm[i]-we[i]*sn[i];fprintf(fp,"第[" );fprintf(fp,"%2d",i+1 );fprintf(fp, "]个节点是平衡节点，其节点功率为：");fprintf(fp,"%9.6f",w_p[i]);fprintf(fp,"+j" );fprintf(fp,"%9.6f",w_q[i]);fprintf(fp,"\n" );} //$$$$$$$$$$$$$$$$$$$$求Jocabi矩阵$$$$$$$$$$$$$$$$$$$$$for(int j=0;j<node_num;j++){if(j==i&&wn[i]==1){Jocabi[2*i][2*j]=-sm[i]-GF1[i][i]*we[i]-BF1[i][i]*wf[i]; Jocabi[2*i][2*j+1]=-sn[i]+BF1[i][i]*we[i]-GF1[i][i]*wf[i]; Jocabi[2*i+1][2*j]=sn[i]+BF1[i][i]*we[i]-GF1[i][i]*wf[i]; Jocabi[2*i+1][2*j+1]=-sm[i]+GF1[i][i]*we[i]+BF1[i][i]*wf[i];} if(j==i&&wn[i]==2){Jocabi[2*i][2*j]=-sm[i]-GF1[i][i]*we[i]-BF1[i][i]*wf[i]; Jocabi[2*i][2*j+1]=-sn[i]+BF1[i][i]*we[i]-GF1[i][i]*wf[i]; Jocabi[2*i+1][2*j]=(-2)*we[i];Jocabi[2*i+1][2*j+1]=(-2)*wf[i];} if(j!=i&&wn[i]==1){Jocabi[2*i][2*j]=(-1)*(GF1[i][j]*we[i]+BF1[i][j]*wf[i]);Jocabi[2*i][2*j+1]=BF1[i][j]*we[j]-GF1[i][j]*wf[i];Jocabi[2*i+1][2*j]=BF1[i][j]*we[j]-GF1[i][j]*wf[i];Jocabi[2*i+1][2*j+1]=GF1[i][j]*we[i]+BF1[i][j]*wf[i];} if(j!=i&&wn[i]==2){Jocabi[2*i][2*j]=(-1)*(GF1[i][j]*we[i]+BF1[i][j]*wf[i]);Jocabi[2*i][2*j+1]=BF1[i][j]*we[j]-GF1[i][j]*wf[i];Jocabi[2*i+1][2*j]=0;Jocabi[2*i+1][2*j+1]=0;}}} //$$$$$$$$$$$$$$$$$$$高斯法解方程$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ double * Dalta_w=new double[Jocabi_rank];for(i=0;i<node_num;i++){if(wn[i]==1){Dalta_w[2*i]=dp[i];Dalta_w[2*i+1]=dq[i];}if(wn[i]==2){Dalta_w[2*i]=dp[i];Dalta_w[2*i+1]=dvv[i];}}double * M=new double[Jocabi_rank];double * X=new double[Jocabi_rank];for(i=0;i<Jocabi_rank;i++){X[i]=0;M[i]=0;}double t;int * js=new int [Jocabi_rank];for(i=0;i<Jocabi_rank;i++)js[i]=0;//交换Jocabi矩阵的行，使每列的对角线元素为所有行中最大的元素double * temp_swap_vector=new double[Jocabi_rank+1];for(inttemp_swap_vector_i=0;temp_swap_vector_i<Jocabi_rank;temp_swap_vector_ i++)temp_swap_vector[temp_swap_vector_i]=0;double temp_variable_for_compare=0;int temp_line=0;for(int diagonal_i=0;diagonal_i<Jocabi_rank;diagonal_i++){temp_variable_for_compare=Jocabi[diagonal_i][diagonal_i];for(intJocabi_line=diagonal_i;Jocabi_line<Jocabi_rank;Jocabi_line++) {if(temp_variable_for_compare<fabs(Jocabi[Jocabi_line][diagonal _i])){temp_variable_for_compare=fabs(Jocabi[Jocabi_line][diagonal_i] );temp_line=Jocabi_line;for(int Jocabi_j_j=0;Jocabi_j_j<Jocabi_rank;Jocabi_j_j++) {temp_swap_vector[Jocabi_j_j]=Jocabi[diagonal_i][Jocabi_j_j];Jocabi[diagonal_i][Jocabi_j_j]=Jocabi[temp_line][Jocabi_j_j];Jocabi[temp_line][Jocabi_j_j]=temp_swap_vector[Jocabi_j_j];}t=Dalta_w[diagonal_i];Dalta_w[diagonal_i]=Dalta_w[temp_line];Dalta_w[temp_line]=t;}}}//$$$$$$$$$$$$$$交换Jocabi矩阵的行完毕$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ double * * temp=new double * [Jocabi_rank];for(i=0;i<Jocabi_rank;i++){temp[i]=new double[Jocabi_rank];for(int j=0;j<Jocabi_rank;j++)temp[i][j]=0.0;}for(int k=0;k<Jocabi_rank;k++){if(Jocabi[k][k]==0) cout<<"fail"<<endl;else{temp[k][k]=Jocabi[k][k];for(int j=0;j<Jocabi_rank;j++)Jocabi[k][j]=Jocabi[k][j]/temp[k][k];Dalta_w[k]=Dalta_w[k]/temp[k][k];for(int ii=k+1;ii<Jocabi_rank;ii++){M[ii]=0.0;M[ii]=Jocabi[ii][k];for(int jj=k;jj<Jocabi_rank;jj++)Jocabi[ii][jj]=Jocabi[ii][jj]-M[ii]*Jocabi[k][jj];Dalta_w[ii]=Dalta_w[ii]-M[ii]*Dalta_w[k];}}} X[Jocabi_rank-1]=Dalta_w[Jocabi_rank-1]/(-Jocabi[Jocabi_rank-1][J ocabi_rank-1]);for(i=Jocabi_rank-2;i>=0;i--){t=0;for(k=i+1;k<Jocabi_rank;k++){t=t+(-Jocabi[i][k])*X[k];}X[i]=(Dalta_w[i]-t)/(-Jocabi[i][i]);}for(i=0;i<node_num;i++){if(wn[i]==1){we[i]=we[i]+X[2*i];wf[i]=wf[i]+X[2*i+1];fprintf(fp,"第[");fprintf(fp,"%2d",i+1);fprintf(fp,"]个节点是PQ节点其节点电压实部修正量是：");fprintf(fp,"%9.6f",X[2*i]);fprintf(fp," 其节点电压虚部修正量是：");fprintf(fp,"%9.6f",X[2*i+1]);fprintf(fp,"\n");} if(wn[i]==2){we[i]=we[i]+X[2*i];wf[i]=wf[i]+X[2*i+1];fprintf(fp,"第[");fprintf(fp,"%2d",i+1);fprintf(fp,"]个节点是PV节点其节点电压实部修正量是：");fprintf(fp,"%9.6f",X[2*i]);fprintf(fp," 其节点电压虚部修正量是：");fprintf(fp,"%9.6f",X[2*i+1]);fprintf(fp,"\n");}} fprintf(fp,"\n");fprintf(fp,"\n");//$$$$$$$$$$$$$$$$判断是否收敛并求最后的差值$$$$$$$$$$$$$$$$$$$$$$ double ch=0;double * * s_p=new double * [node_num+1];for(int s_p_i=1;s_p_i<=node_num+1;s_p_i++){s_p[s_p_i]=new double[node_num+1];for(int s_p_j=1;s_p_j<=node_num+1;s_p_j++)s_p[s_p_i][s_p_j]=0;}double * * s_q=new double * [node_num+1];for(int s_q_i=1;s_q_i<=node_num+1;s_q_i++){s_q[s_q_i]=new double[node_num+1];for(int s_q_j=1;s_q_j<=node_num+1;s_q_j++)s_q[s_q_i][s_q_j]=0;}for(int ci=1;ci<=node_num-1;ci++){if((fabs(Dalta_w[ci]))>(fabs(ch))){ch=fabs(Dalta_w[ci]);}}fprintf(fp,"判断是否迭代的最大值是:");fprintf(fp,"%9.6f",ch);fprintf(fp,"\n");fprintf(fp,"\n");if(ch<1e-005){double *UU=new double [node_num];for(int k=0;k<node_num;k++)UU[k]=0.0;for(k=0;k<node_num;k++)UU[k]=we[k]*we[k]+wf[k]*wf[k];double *G_ground=new double [n];double *B_ground=new double [n];for(int i=0;i<n;i++){G_ground[i]=0;B_ground[i]=0;}cout<<"满足迭代要求，计算完成，迭代次数是"<<iterative_frequency<<endl;fprintf(fp,"||||||||||||||||||||||||||||||||||||||||||||||||||||| ||||||||\n");fprintf(fp,"\n");fprintf(fp,"满足迭代要求，计算完成，迭代次数是");fprintf(fp,"%2d",iterative_frequency);fprintf(fp,"\n");for(i=0;i<node_num;i++){fprintf(fp,"节点");fprintf(fp,"[");fprintf(fp,"%2d",i+1);fprintf(fp,"]的电压是：");fprintf(fp,"%9.6f",we[i]);if(wf[i]>0){fprintf(fp," + j");fprintf(fp,"%9.6f",wf[i]);}if(wf[i]<0){fprintf(fp," - j");fprintf(fp,"%9.6f",fabs(wf[i]));}fprintf(fp," 幅值是：");fprintf(fp,"%9.6f",sqrt(we[i]*we[i]+wf[i]*wf[i]));fprintf(fp," 相角是：");fprintf(fp,"%9.6f",atan(wf[i]/we[i])*180/3.1415926);fprintf(fp,"度");fprintf(fp,"\n");}fprintf(fp,"线路功率是:\n");for(int j=0;j<n;j++){double K=0.0;int p=0;int q=0;p=II[j];q=JJ[j];if(p==0&&q==0)break;GF[j]=WR[j]/(WR[j]*WR[j]+WX[j]*WX[j]);BF[j]=(-WX[j])/(WR[j]*WR[j]+WX[j]*WX[j]);if(WNT[j]<0){K=fabs(WNT[j]);G_ground[j]=((1-K)/(K*K))*GF[j];B_ground[j]=((1-K)/(K*K))*BF[j];}if(WNT[j]>0){K=WNT[j];G_ground[j]=((K-1)/K)*GF[j];B_ground[j]=((K-1)/K)*BF[j];}if(WXC[j]<1000)B_ground[j]=B_ground[j]+WXC[j];s_p[p][q]=UU[p-1]*G_ground[j]+UU[p-1]*GF1[p-1][q-1]-(we[p-1]*we[q -1]+wf[p-1]*wf[q-1])*GF1[p-1][q-1]+(we[p-1]*wf[q-1]-we[q-1]*wf[p-1])*BF1[p-1][q-1];s_q[p][q]=-UU[p-1]*B_ground[j]-UU[p-1]*BF1[p-1][q-1]+(we[p-1]*we[ q-1]+wf[q-1]*wf[q-1])*BF1[p-1][q-1]+(we[p-1]*wf[q-1]-we[q-1]*wf[p-1])*GF1[p-1][q-1];fprintf(fp,"S[");fprintf(fp,"%2d",p);fprintf(fp,"][");fprintf(fp,"%2d",q);fprintf(fp,"]=");fprintf(fp,"%9.6f",s_p[p][q]);if(s_q[p][q]>0){fprintf(fp," + j");fprintf(fp,"%9.6f",s_q[p][q]);}if(s_q[p][q]<0){fprintf(fp," - j");fprintf(fp,"%9.6f",fabs(s_q[p][q]));}fprintf(fp,"\n");}break;}}//迭代结束八、IEEE30节点的仿真结果1、节点电压编号电压实部电压虚部电压幅值电压相角节点1 1.05 0i 1.05 0度节点2 1.04445 -0.0337697i 1.045 -1.85186度节点3 1.01843 -0.067472i 1.02066 -3.79037度节点4 1.01033 -0.0798977i 1.01348 -4.52159度节点5 1.0035 -0.114413i 1.01 -6.50445度节点6 1.00662 -0.0943613i 1.01104 -5.35528度节点7 0.996459 -0.111039i 1.00263 -6.35844度节点8 1.00511 -0.0993071i 1.01 -5.64266度节点9 1.03035 -0.122151i 1.03756 -6.76104度节点10 1.02229 -0.156054i 1.03413 -8.67927度节点11 1.04666 -0.0837098i 1.05 -4.57268度节点12 1.03608 -0.141734i 1.04573 -7.78959度节点13 1.0436 -0.115738i 1.05 -6.32836度节点14 1.01898 -0.156444i 1.03092 -8.72848度节点15 1.01403 -0.157951i 1.02626 -8.85357度节点16 1.02183 -0.151553i 1.033 -8.43636度节点17 1.01658 -0.157803i 1.02876 -8.82354度节点18 1.00279 -0.167848i 1.01674 -9.50214度节点19 0.999803 -0.170757i 1.01428 -9.69206度节点20 1.00449 -0.16805i 1.01845 -9.49755度节点21 1.00872 -0.162656i 1.02175 -9.16004度节点22 1.00931 -0.162655i 1.02234 -9.15476度节点23 1.00254 -0.16504i 1.01604 -9.34827度节点24 0.996609 -0.169594i 1.01094 -9.65755度节点25 0.995082 -0.169371i 1.00939 -9.65966度节点26 0.976247 -0.173651i 0.991571 -10.0861度节点27 1.00354 -0.166036i 1.01718 -9.3945度节点28 1.00168 -0.102188i 1.00688 -5.82495度节点29 0.980068 -0.184114i 0.997212 -10.6395度节点30 0.96576 -0.197073i 0.985662 -11.5334度2、线路传输功率首端—末端首末功率末首功率2-1 -57.9109+(8.36395)i 58.5228+(-12.325)i3-1 -39.536+(-7.37948)i 40.2261+(5.83313)i4-2 -30.8909+(-9.48018)i 31.4524+(7.29204)i4-3 -36.9557+(-6.53076)i 37.136+(6.17948)i5-2 -44.4784+(-7.6178)i 45.4077+(7.10769)i6-2 -38.4904+(-7.58429)i 39.3508+(6.2414)i6-4 -34.7811+(4.04816)i 34.9241+(-4.02213)i 7-5 -0.260435+(-7.29197)i 0.278436+(5.2715)i7-6 -22.5396+(-3.60803)i 22.6765+(2.30528)i8-6 -11.9244+(0.486521)i 11.9413+(-1.34664)i 6-9 12.3726+(-12.7417)i -12.3726+(13.3835)i6-10 10.9034+(-3.88366)i -10.9034+(4.61235)i 11-9 20+(6.66081)i -20+(-5.82246)i10-9 -32.6509+(-2.67651)i 32.6509+(3.78043)i4-12 23.6005+(-12.0947)i -23.6005+(13.8474)i 13-12 20+(3.45425)i -20+(-2.93116)i14-12 -7.91606+(-2.10553)i 7.99377+(2.26709)i15-12 -18.28+(-5.90486)i 18.5119+(6.36176)i16-12 -7.55825+(-2.98863)i 7.61675+(3.11163)i15-14 -1.7094+(-0.499524)i 1.71605+(0.505537)i 17-16 -4.04947+(-1.1564)i 4.05825+(1.18863)i18-15 -6.09617+(-1.40631)i 6.13679+(1.48904)i19-18 -2.89082+(-0.495509)i 2.89617+(0.506313)i 20-19 6.6264+(2.93893)i -6.60918+(-2.90448)i 20-10 -8.82641+(-3.63893)i 8.90866+(3.82259)i17-10 -4.95054+(-4.64359)i 4.96464+(4.68038)i21-10 -16.1704+(-9.32613)i 16.2865+(9.57613)i22-10 -7.88783+(-4.19568)i 7.94335+(4.31016)i22-21 1.33022+(1.87505)i -1.32963+(-1.87386)i 23-15 -5.6167+(-2.34288)i 5.65258+(2.41535)i24-22 -6.50437+(-2.23775)i 6.55761+(2.32062)i24-23 -2.40853+(-0.726163)i 2.4167+(0.742882)i25-24 -0.212582+(-0.351348)i 0.212894+(0.351893)i 26-25 -3.5+(-2.29999)i 3.54539+(2.36778)i27-25 3.34908+(2.04752)i -3.3328+(-2.01643)i 28-27 16.1025+(-2.11682)i -16.1025+(3.14712)i 29-27 -6.10422+(-1.50686)i 6.1916+(1.67196)i30-27 -6.92977+(-1.35732)i 7.09411+(1.66665)i30-29 -3.67023+(-0.542645)i 3.70422+(0.60687)i28-8 -1.92156+(-3.12492)i 1.92445+(-1.21858)i28-6 -14.7132+(-3.43133)i 14.7506+(2.24037)i 3、发电机注入的功率编号功率实部功率虚部节点1 98.7489 -6.4919i 节点2 80 41.7051i 节点5 50 16.6537i 节点8 20 29.2679i 节点11 20 6.66081i 节点13 20 3.45425i。

潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。

2.在给定的电力网络上画出等值电路图。

3.运用计算机进行潮流计算。

4.编写设计说明书。

一、设计原理1．牛顿-拉夫逊原理牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。

牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。

电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。

为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。

牛顿—拉夫逊迭代法的一般步骤：（1）形成各节点导纳矩阵Y。

（2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。

（3）计算各个节点的功率不平衡量。

（4）根据收敛条件判断是否满足，若不满足则向下进行。

（5）计算雅可比矩阵中的各元素。

（6）修正方程式个节点电压（7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。

（8）计算平衡节点输出功率和各线路功率2．网络节点的优化１）静态地按最少出线支路数编号这种方法由称为静态优化法。