完整word版潮流计算方法

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

本人一个星期以来的汗水，弄清楚了计算机算法计算潮流的基础，如果有什么不懂的可以发信息到邮箱：zenghao616@https://www.360docs.net/doc/248360754.html,

接下来开始弄潮流的优化问题，吼吼！

电力系统的潮流计算的计算机算法：以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:5

for j=1:1:5

y(j,i)=y(i,j);

end

end

%节点导纳矩阵的形成

Y=zeros(5,5);

%求互导纳

for i=1:1:5

for j=1:1:5

if i~=j

Y(i,j)=-y(i,j);

end

end

end

%求自导纳

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 Q

P = 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:5

for j=1:1:i

if i~=j

Sum_YU0 = Sum_YU0 + Y(i,j)*U_reg(j);

end

end

for j=i+1:1:5

Sum_YU1 = Sum_YU1 + Y(i,j)*U(j);

end

U(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==2

angle_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)) ) );end

i==3 if angle_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)

conj(Y(3,3)*U(3)) ) );end下面做越界检查%

%if Q(4)>Q_Max

% Q(4) = Q_Max;%end

%if Q(4)

% Q(4) = Q_Min;

%end

%下面可以做PV节点收敛判断

Sum_YU0 = 0;

Sum_YU1 = 0;

end

end

%节点注入无功，流入为正，流出为负

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:5

S_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:5

for j=1:1:5

y(j,i)=y(i,j);

end

end

%节点导纳矩阵的形成

Y=zeros(5,5);

%求互导纳

for i=1:1:5

for j=1:1:5

if i~=j

Y(i,j)=-y(i,j);

end

end

end

%求自导纳

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 Q

P = 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_U

i=2:1:5for

for j=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 %不为节点2,3则计算无功

delta_Q(i)=Q(i)-f(i)*Sum_GB1+e(i)*Sum_GB2;

end

if i==2 || i==3 %这里计算delta_U的值，始终为零

delta_U(i)=U(i)^2-( e(i)^2 + f(i)^2 );

end

Sum_GB1=0;Sum_GB2=0;

end

%___________________________________%

%下面计算雅克比矩阵

J=zeros(8,8);

for ii=2:1:5

i=ii-1;

for j=1:1:5

Sum_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节点

if ii==jj

J(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);

else

J(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));

end

else%PV节点

if ii==jj

J(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);

else

J(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;

end

end

end

Sum_GB1=0;Sum_GB2=0;

end