中国石油大学2010届电力系统分析大作业

电力系统分析大作业

牛顿拉夫逊潮流计算

电气1002

1、牛顿-拉夫逊法潮流算法的基本原理

电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部的运行参数，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。

2、节点导纳矩阵的算法

根据定义直接求取节点导纳矩阵时，注意以下几点：

(1) 节点导纳矩阵是方阵，其阶数就等于网络中除去参考节点外的节点数。参考节点一般取大地，编号为零。

(2) 节点导纳矩阵是稀疏矩阵，其各行非零非对角元素就等于与该行相对应节点所连接的不接地支路数。

(3) 节点导纳矩阵的对角元素就等于该节点所连接导纳的总和。因此，与没有接地支路的节点对应的行或列中，对角元素为非对角元素之和的负值。

(4) 节点导纳矩阵的非对角元素等于连接节点i，j支路导纳的负值。因此，在一般情况下，节点导纳矩阵的对角元素往往大于非对角元素的负值。

(5) 节点导纳矩阵一般是对称矩阵，这是网络的互易特性所决定的。从而，一般只要求取这个矩阵的上三角或下三角部分。

(6) 网络中的变压器。

5. 牛顿—拉夫逊迭代法的步骤及算法流程图：

（1）形成各节点导纳矩阵Y。

（2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。（3）计算各个节点的功率不平衡量。

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

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

（6）修正方程式个节点电压

（7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。（8）计算平衡节点输出功率和各线路功率

流程图如下：

6、心得体会

把潮流计算的任务完成后，有一些体会，主要体现在以下三个方面：

（1）、对Matlab这个工具本身有了一些了解，包括它的集成开发环境，调试运行方式，程序的基本语句以及控制方式。其中感受最深刻的是Matlab的数据结构，很多书上介绍说Matlab的数据单元是矩阵，但自己看了以后并没有什么体会，在编写调试程序时，出于C 语言的习惯，我们把节点的初始电压，节点的输出功率都设置成单个数据的形式，但在写功率平衡方程用到For语句时，单个数据就很不方便，于是又改成了矩阵形式，在后续的程序中，也均使用矩阵做为数据单元。通过编写Matlab程序，对“Matlab的基本数据单元是矩阵”这句话有了一定的了解。

（2）、对潮流计算也比以前有了进一步的体会，在学习潮流计算时，虽然依次学习了节点导纳矩阵，功率方程、雅可比矩阵，但不能将它们联系起来，更不知道其中的原委，通过程序的编写，知道了其中的联系，也知道了每个方程、矩阵在计算中的作用。

（3）、在接到这个课外课题时，由于以前没有用过Matlab，刚开始感觉有些茫然，不知道从何做起，但后来通过看书，对Matlab有了一定了解，又觉得和C语言差不多，没什么难的，真正开始写程序时，发现既不是我们想的那么简单，也不是难得无法下手，说它不是那么简单，是因为Matlab和C语言有一定区别，它的数据结构主要是矩阵，另外功率方程里的求和部分也不像WORD里那样，可以用公式编辑器写出，得用循环实现。说它不是难得无法下手，是因为通过我们查阅资料和自己调试，最后完成了潮流计算的程序并得到了收敛的结果。

7、附上程序代码及详细注释：

%本程序是用牛顿-拉夫逊法极坐标进行潮流计算.

%第一列为节点注入参数

%第二列为节点负荷功率参数，第三列为节点电压参数，第六列为节点类型参数

%其中‘1’为平衡节点，‘2’为PQ节点，‘3’为PV节点.

%B1是支路参数矩阵，第一列和第二列是节点编号，B2为节点参数矩阵.

%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号.

%第三列为之路串联阻抗参数.

%第四列为含支路的对地导纳矩阵.

%第五列为含变压器之路的变压器变比.

%第六列表式是否含有变压器的参数，‘1’表示含有变压器，‘0’表示不含有.

%源程序代码如下：

clc%清屏

clear

disp('各输出结果如下：')

n=5;

nl=5;

isb=5;

pr=0.0001;

B1=[1 2 0.025+0.08i 0.07i 1 0;

1 3 0.03+0.1i 0.09 1 0;

2 3 0.02+0.06i 0.05i 1 0;

4 2 0.1905i 0 1.0522 1;

5 3 0.1905i 0 1.0522 1];%变压器侧为1，否则为0

B2=[0 -0.8055-0.5320i 1 0 0 2;

0 -0.18-0.12i 1 0 0 2;

0 0 1 0 0 2;

0 0.5 1.0522 1.0522 0 3;

0 0 1.0522 1.0522 0 1];%题中已知矩阵

disp('——————【已知量：】——————')

disp('【支路参数矩阵】')

disp(B1)

disp('【节点参数矩阵】')

disp(B2)

X=[1 0;2 0;3 0;4 0;5 0];%输入节点序号

Y=zeros(n);U=zeros(1,n);cta=zeros(1,n);V=zeros(1,n);O=zeros(1,n);S1=zeros(nl); for i=1:n

if X(i,2)~=0;

p=X(i,1);

Y(p,p)=X(i,2);

end

end

%————————————求导纳矩阵———————————————

for i=1:nl%从1到n1（总支路数）

if B1(i,6)==0

p=B1(i,1);q=B1(i,2);

else p=B1(i,2);q=B1(i,1);

end

Y(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;%对角元

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;%对角元

end

disp('——————【计算量：】——————')

disp('【导纳矩阵】 Y=');

disp(Y)%输出导纳矩阵

%-----------给定各节点初始电压及给定各节点注入功率-------------

G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部

for i=1:n

cta(i)=angle(B2(i,3));

U(i)=abs(B2(i,3));

V(i)=B2(i,4);%PV、平衡节点及PQ节点电压模值

end

for i=1:n%给定各节点注入功率

S(i)=B2(i,1)-B2(i,2);%i节点注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量(电纳值)

end

%=============用牛顿-拉夫逊法迭代求解非线性代数方程（功率方程）============== P=real(S);Q=imag(S);%分解出各节点注入的有功和无功功率

ICT1=0;IT2=1;%迭代次数ICT1；不满足收敛要求的节点数IT2

while IT2~=0

IT2=0;t1=1;t2=1;

for i=1:n

if i~=isb

C(i)=0;

D(i)=0;

for j1=1:n %第i行共n列(n个节点间互导纳及节点电压相乘即电流)

C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));

D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1))); end

DP(t1)=P(i)-C(i);

t1=t1+1;

if B2(i,6)==2

DQ(t2)=Q(i)-D(i);

t2=t2+1;

end

end

end

t1=t1-1;t2=t2-1;

DPQ=[DP';DQ']; %求DP,DQ

%============== 以下为形成完整的Jacobi矩阵 ===============

for i=1:t1+t2

if abs(DPQ(i))>pr

IT2=IT2+1;

end

end

H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2);

for i=1:t1

for j1=1:t1

if j1~=isb&j1~=i

H(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1))); elseif j1~=isb&j1==i

H(i,j1)=U(i)^2*B(i,j1)+D(i);

end

end

for i=1:t1

for j1=1:t2

if j1~=isb&j1~=i

N(i,j1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1))); elseif j1~=isb&j1==i

N(i,j1)=0-U(i)^2*G(i,j1)-C(i);

end

end

end

for i=1:t2

for j1=1:t1

if j1~=isb&j1~=i

K(i,j1)=

U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));

elseif j1~=isb&j1==i

K(i,j1)=U(i)^2*G(i,j1)-C(i);

end

end

end

for i=1:t2

for j1=1:t2

if j1~=isb&j1~=i

L(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1))); elseif j1~=isb&j1==i

L(i,j1)=U(i)^2*B(i,j1)-D(i);

end

end

end

J=[H,N;K,L];%输出雅可比矩阵

disp('【雅可比矩阵：】');

disp(J);

modify=-J\DPQ;

Dcta=modify([1:t1],:);

t3=U(:,[1:t2]);

DU=diag(t3,0)*modify([t1+1:t1+t2],:);

t4=1;

for i=1:t1

if B2(i,6)~=1

cta(1,i)=cta(1,i)+Dcta(t4,1);

t4=t4+1;

end

t5=1;

for i=1:t2

if B2(i,6)==2

U(1,i)=U(1,i)+DU(t5,1);

t5=t5+1;

end

end

ICT1=ICT1+1;

end %求修正值

for i=1:n

UU(i)=U(i)*cos(cta(i))+1i*U(i)*sin(cta(i));

end

for p=1:n

c(p)=0;

for q=1:n

c(p)=c(p)+conj(Y(p,q))*conj(UU(q));

end

s(p)=UU(p)*c(p);

end

disp('-------------------------------------------------------------------'); disp('【各节点电压U为（节点从小到大排列）:】');

disp(UU);

disp('-------------------------------------------------------------------'); disp('【各节点电压相角为（节点从小到大排列）：】');

disp(180*angle(UU)/pi);

disp('-------------------------------------------------------------------'); disp('【按公式计算全部线路功率，结果如下:】');

for i=1:nl

if B1(i,6)==0

p=B1(i,1);q=B1(i,2);

else p=B1(i,2);q=B1(i,1);

end

Si(p,q)=UU(p)*(conj(UU(p))*conj(B1(i,4)./2)+(conj(UU(p)*B1(i,5))-conj(UU(q)))*c onj(1./(B1(i,3)*B1(i,5))));%各条支路首端功率Si

f=[p,q,Si(p,q)];

disp(f);

end

for i=1:nl

if B1(i,6)==0

p=B1(i,1);q=B1(i,2);

else p=B1(i,2);q=B1(i,1);

end

Sj(q,p)=UU(q)*(conj(UU(q))*conj(B1(i,4)./2)+(conj(UU(q)./B1(i,5))-conj(UU(p)))* conj(1./(B1(i,3)*B1(i,5))));%各条支路末端功率Sj

f=[q,p,Sj(q,p)];

disp(f);

end

disp('------------------------------------------------------------------'); disp('【各条支路的功率损耗DS为(顺序和输入的B1时一样):】');

for i=1:nl

if B1(i,6)==0

p=B1(i,1);q=B1(i,2);

else p=B1(i,2);q=B1(i,1);

end

DS(i)=Si(p,q)+Sj(q,p);%各条支路功率损耗DS

disp(DS(i));

end

Sp=0;%最后求平衡节点功率

for i=1:n

Sp=Sp+UU(isb)*conj(Y(isb,i))*conj(UU(i));

end

disp('【平衡节点的功率：】');

disp(Sp);

8、程序运行结果截图(此为matlab2013运行截图)