潮流分布的计算

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

本人一个星期以来的汗水，弄清楚了计算机算法计算潮流的基础，如果有什么不懂的可以发信息到邮箱：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算法计算潮流：这个算法是由牛顿算法的极坐标形式简化而来。

牛顿拉夫逊法计算潮流步骤牛顿拉夫逊法（Newton-Raphson Method）是一种常用于计算潮流的数值求解方法。

它是基于潮流计算的功率流方程的非线性特性而设计的，通过迭代求解来逼近潮流计算的稳态解。

下面将介绍牛顿拉夫逊法计算潮流的基本步骤。

首先，我们需要明确潮流计算的目标，即确定电力系统中各节点的电压相角和幅值。

这些节点是电力系统中的发电机、负荷和交流输电线路的连接点。

通过潮流计算，我们可以得到各节点的电压相角和幅值，从而分析系统的功率分布、电压稳定性等运行特性。

接下来，我们需要建立电力系统的潮流计算模型。

这个模型中，我们需要考虑发电机的注入功率、负荷的吸收功率、线路的传输损耗等因素。

通过利用功率流方程，我们可以将这些因素表示为电压、功率和导纳之间的方程。

然后，我们需要进行初始化操作。

在进行牛顿拉夫逊法迭代计算之前，我们需要对电力系统的各节点进行初始电压值的设定。

这些初始值可以根据经验或者历史数据来得到，但需要满足物理约束条件，如一致性、电压幅值在合理范围内等。

接下来，我们进入迭代计算的过程。

首先，我们需要对系统的节点进行编号，然后选择某一节点作为基准节点，其他节点相对于基准节点的电压相角进行计算。

然后，我们根据节点注入功率和导纳矩阵的关系，得到节点注入电流。

接着，我们根据节点注入电流和电压相角的关系，计算各节点的电压相角和幅值的改变量。

在计算改变量后，我们需要对节点电压进行更新。

更新后，我们判断系统是否达到收敛条件。

如果满足收敛条件，则停止迭代，得到最终的潮流计算结果；如果不满足收敛条件，则继续进行下一轮迭代计算。

最后，我们对潮流计算结果进行分析和验证。

通过比较计算得到的结果和实际运行数据进行对比，我们可以评估潮流计算的准确性。

同时，我们还可以通过故障分析、电压稳定性评估等手段对电力系统进行优化和改进。

总而言之，牛顿拉夫逊法是一种常用的求解潮流计算问题的方法。

它通过迭代求解潮流计算的功率流方程，逼近潮流计算的稳态解。

~ S LDc

j
B2 2
U
2 N
S~b

S~LDb

j
B1 2
U
2 N

j
B2 2
U
2 N
由此将问题转化为：已知
U A ,
j
B1 2
U
2 N
,
S~b ,
S~c
的潮流计算。
~
A SA
~ S1
S~1
S~1
b
~ S2
S~2
S~2
c
U A
Z1
Z2
a.反推功率：

j
B1 2
UHale Waihona Puke 2 NS~bS~c
~ S1
①
S~1
S~2
I1
I1 Z
B j
S~Y 1
2
S~2 ②
I2
B j
2
~ S2
U 2
S~Y 2
求导纳中的功 率损耗S~Y1，S~Y 2；
末端：S~Y 2

U 2
(
j
B 2
U 2 )


j
B 2

U
2 2
首端：S~Y 1

U 1

(
j
B 2
U1 )
jB
~ S LD

30
j15MVA
2
~ SY 2
已知 r1 0.27 / km, x1 0.423 / km
b1 2.69 106 s / km, l 150km, 双回线路
解：R 1 0.27150 20.25 X 1 0.423150 31.725

电力系统潮流计算电力系统潮流计算是电力系统运行分析中的重要环节。

它通过对电力系统中各节点的电压、相角以及功率等参数进行计算和分析，从而得出电力系统的稳态运行状态。

本文将从潮流计算的基本原理、计算方法、应用及其发展等方面进行阐述。

一、潮流计算的基本原理电力系统潮流计算的基本原理是基于潮流方程建立的。

潮流方程是一组非线性的方程，描述了电力系统中各节点的电压、相角以及功率之间的关系。

潮流计算的目的就是求解这组非线性方程，以确定电力系统的电压幅值、相角及有功、无功功率的分布情况。

二、潮流计算的基本方法潮流计算的基本方法主要有直接法、迭代法以及牛顿-拉夫逊法。

直接法是通过直接求解潮流方程得到电力系统的潮流状况，但对于大规模复杂的电力系统来说，直接法计算复杂度高。

迭代法是通过对电力系统的节点逐个进行迭代计算，直到满足预设的收敛条件。

牛顿-拉夫逊法是一种较为高效的迭代法，它通过近似潮流方程的雅可比矩阵，实现了计算的高效和稳定。

三、潮流计算的应用潮流计算在电力系统运行与规划中起着重要作用。

首先，潮流计算可以用于电力系统的稳态分析，确定电力系统在各种工况下的电压、相角等参数，以判断电力系统是否存在潮流拥挤、电压失调等问题。

其次，潮流计算还可以用于电力系统的优化调度，通过调整电力系统的发电机出力、负荷组织等参数，以改善电力系统的经济性和可靠性。

此外，潮流计算还可以用于电力系统规划，通过对电力系统进行潮流计算，可以为新建电源、输电线路以及变电站等设备的规划和选择提供科学依据。

四、潮流计算的发展随着电力系统的规模不断扩大和复杂度的提高，潮流计算技术也得到了迅速的发展。

传统的潮流计算方法在计算效率和计算精度上存在一定的局限性。

因此，近年来研究者提出了基于改进的迭代方法、高精度的求解算法以及并行计算等技术，以提高潮流计算的速度和准确性。

此外，随着可再生能源的不断融入电力系统，潮流计算还需要考虑多种能源的互联互通问题，这对潮流计算提出了新的挑战，需要进一步的研究和改进。