电力系统分析课程设计牛顿拉夫逊潮流计算.docx
《电力系统分析课程设计牛顿拉夫逊潮流计算.docx》由会员分享,可在线阅读,更多相关《电力系统分析课程设计牛顿拉夫逊潮流计算.docx(21页珍藏版)》请在冰豆网上搜索。
![电力系统分析课程设计牛顿拉夫逊潮流计算.docx](https://file1.bdocx.com/fileroot1/2022-10/9/25b6c50b-7de9-4bbb-b710-02d4370aca3f/25b6c50b-7de9-4bbb-b710-02d4370aca3f1.gif)
电力系统分析课程设计牛顿拉夫逊潮流计算
1设计任务及初步分析
1.1设计任务
条件:
节点数:
3支路数:
3计算精度:
0.00010
支路1:
0.0300+j0.0900
1┠—————□—————┨2
支路2:
0.0200+j0.0900
2┠—————□—————┨3
支路3:
0.0300+j0.0900
3┠—————□—————┨1
节点1:
PQ节点,S
(1)=-0.5000-j0.2000
节点2:
PQ节点,S
(2)=-0.6000-j0.2500
节点3:
平衡节点,U(3)=1.0000∠0.0000
要求:
编写程序计算潮流
1.2初步分析
潮流计算在数学上可归结为求解非线性方程组,其数学模型简写如下:
2牛顿-拉夫逊法简介
2.1概述
牛顿-拉夫逊法是目前求解非线性方程最好的一种方法。
这种方法的特点就是把对非线性方程的求解过程变成反复对相应的线性方程求解的过程,通常称为逐次线性化过程,就是牛顿-拉夫逊法的核心。
牛顿-拉夫逊法的基本原理是在解的某一邻域内的某一初始点出发,沿着该点的一阶偏导数——雅可比矩阵,朝减小方程的残差的方向前进一步,在新的点上再计算残差和雅可矩阵继续前进,重复这一过程直到残差达到收敛标准,即得到了非线性方程组的解。
因为越靠近解,偏导数的方向越准,收敛速度也越快,所以牛顿法具有二阶收敛特性。
而所谓“某一邻域”是指雅可比方向均指向解的范围,否则可能走向非线性函数的其它极值点,一般来说潮流由平电压即各母线电压(相角为0,幅值为1)启动即在此邻域内。
2.2一般概念
对于非线性代数方程组
即
(2-1)
在待求量
的某一个初始计算值
附件,将上式展开泰勒级数并略去二阶及以上的高阶项,得到如下的线性化的方程组
(2-2)
上式称之为牛顿法的修正方程式。
由此可以求得第一次迭代的修正量
(2-3)
将
和
相加,得到变量的第一次改进值
。
接着再从
出发,重复上述计算过程。
因此从一定的初值
出发,应用牛顿法求解的迭代格式为
(2-4)
(2-5)
上两式中:
是函数
对于变量
的一阶偏导数矩阵,即雅可比矩阵
;
为迭代次数。
由式(2-4)和式子(2-5)可见,牛顿法的核心便是反复形成求解修正方程式。
牛顿法当初始估计值
和方程的精确解足够接近时,收敛速度非常快,具有平方收敛特性。
2.3潮流计算的修正方程
运用牛顿-拉夫逊法计算潮流分布时,首先要找出描述电力系统的非线性方程。
这里仍从节点电压方程入手,设电力系统导纳矩阵已知,则系统中某节点(
节点)电压方程为
从而得
进而有
(2-6)
式(2-6)中,左边第一项为给定的节点注入功率,第二项为由节点电压求得的节点注入功率。
他们二者之差就是节点功率的不平衡量。
现在有待解决的问题就是各节点功率的不平衡量都趋近于零时,各节点电压应具有的价值。
由此可见,如将式(2-6)作为牛顿-拉夫逊中的非线性函数
,其中节点电压就相当于变量
。
建立了这种对应关系,就可列出修正方程式,并迭代求解。
但由于节点电压可有两种表示方式——以直角做表或者极坐标表示,因而列出的迭代方程相应地也有两种,下面分别讨论。
2.4直角坐标表示的修正方程
节点电压以直角坐标表示时,令
、
,且将导纳矩阵中元素表示为
,则式(2-7)改变为
(2-7)
再将实部和虚部分开,可得
(2-8)
这就是直角坐标下的功率方程。
可见,一个节点列出了有功和无功两个方程。
对于
节点(
),给定量为节点注入功率,记为
、
,则由式(2-8)可得功率的不平衡量,作为非线性方程
(2-9)
式中
、
——分别表示第
节点的有功功率的不平衡量和无功功率的不平衡量。
对于
节点(
),给定量为节点注入有功功率及电压数值,记为
、
,因此,可以利用有功功率的不平衡量和电压的不平衡量表示出非线性方程,即有
(2-10)
式中
为电压的不平衡量。
对于平衡节点(
),因为电压数值及相位角给定,所以
也确定,不需要参加迭代求节点电压。
因此,对于
个节点的系统只能列出
个方程,其中有功功率方程
个,无功功率方程
个,电压方程
个。
将式(2-9)、式(2-10)非线性方程联立,称为
个节点系统的非线性方程组,且按泰勒级数在
、
(
)展开,并略去高次项,得到以矩阵形式表示的修正方程如下。
(2-11)
上式中雅可比矩阵的各个元素则分别为
将(2-11)写成缩写形式
(2-12)
对雅可比矩阵各元素可做如下讨论:
当
时,对于特定的
,只有该特定点的
和
是变量,于是雅可比矩阵中各非对角元素表示为
当
时,雅可比矩阵中各对角元素的表示式为
由上述表达式可知,直角坐标的雅可比矩阵有以下特点:
1)雅可比矩阵是
阶方阵,由于
、
等等,所以它是一个不对称的方阵。
2)雅可比矩阵中诸元素是节点电压的函数,在迭代过程中随电压的变化而不断地改变。
3)雅可比矩阵的非对角元素与节点导纳矩阵
中对应的非对角元素有关,当
中的
为零时,雅可比矩阵中相应的
、
、
、
也都为零,因此,雅可比矩阵也是一个稀疏矩阵。
3程序设计
3.1程序流程图
图3-1程序流程图
3.2源程序
n=3
n1=3
isb=3
pr=0.0001
B1=[120.03+0.09i010;130.03+0.09i010;230.02+0.09i010]
B2=[0-0.5-0.2i1002;0-0.6-0.25i1002;001101]
X=[10;20;30]
%X=input('节点号和对地参数:
X=');
Y=zeros(n);
Times=1;%置迭代次数为初始值
%创建节点导纳矩阵
fori=1:
n1
ifB1(i,6)==0%不含变压器的支路
p=B1(i,1);
q=B1(i,2);
Y(p,q)=Y(p,q)-1/B1(i,3);
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);
Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);
else%含有变压器的支路
p=B1(i,1);
q=B1(i,2);
Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1/B1(i,3);
Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));
end
end
Y
OrgS=zeros(2*n-2,1);
DetaS=zeros(2*n-2,1);%将OrgS、DetaS初始化
%创建OrgS,用于存储初始功率参数
h=0;
j=0;
fori=1:
n%对PQ节点的处理
ifi~=isb&B2(i,6)==2
h=h+1;
forj=1:
n
OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
end
end
end
fori=1:
n%对PV节点的处理,注意这时不可再将h初始化为0
ifi~=isb&B2(i,6)==3
h=h+1;
forj=1:
n
OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
end
end
end
OrgS
%创建PVU用于存储PV节点的初始电压
PVU=zeros(n-h-1,1);
t=0;
fori=1:
n
ifB2(i,6)==3
t=t+1;
PVU(t,1)=B2(i,3);
end
end
PVU
%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量
h=0;
fori=1:
n%对PQ节点的处理
ifi~=isb&B2(i,6)==2
h=h+1;
DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);
end
end
t=0;
fori=1:
n%对PV节点的处理,注意这时不可再将h初始化为0
ifi~=isb&B2(i,6)==3
h=h+1;
t=t+1;
DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;
end
end
DetaS
%创建I,用于存储节点电流参数
i=zeros(n-1,1);
h=0;
fori=1:
n
ifi~=isb
h=h+1;
I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));
end
end
I
%创建Jacbi(雅可比矩阵)
Jacbi=zeros(2*n-2);
h=0;
k=0;
fori=1:
n%对PQ节点的处理
ifB2(i,6)==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==j%对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag