牛拉法潮流计算.docx
- 文档编号:2144348
- 上传时间:2022-10-27
- 格式:DOCX
- 页数:16
- 大小:19.51KB
牛拉法潮流计算.docx
《牛拉法潮流计算.docx》由会员分享,可在线阅读,更多相关《牛拉法潮流计算.docx(16页珍藏版)》请在冰豆网上搜索。
牛拉法潮流计算
自动化07-1班段佳
functionnl;
%------------------------------------------------------------------------
%===================================================================
%======================牛顿——拉夫逊法==============================
%===========================潮流计算=================================
%===================================================================
%-----------------------------------------------------------------------
%%%---------------使用说明部分---------------------------
display('%%本程序的功能是用牛顿——拉夫逊法进行潮流计算');
display('%%本程序要求用户按照一定的格式将电力系统的参数制成excel表格,系统运行时将从excel中加载这些参数,随后后即可进行潮流计算');
display('%%为了方便运算,用户再给系统节点进行编号时,请按照先PQ节点,再PV节点,最后平衡节点的顺序从小到大编号');
display('%%电力系统潮流计算excel格式——支路参数:
1、支路首端号;2、末端号;3、支路阻抗;4、支路对地电纳;5、支路的变比K:
1;6、支路首端处于K侧为1,1侧为0');
display('%%电力系统潮流计算excel格式——节点参数:
1、节点号;2、电压大小;3、相位角;4、发电机有功;5、发电机无功;6、负载有功;7、负载无功;8、节点类型');
%===================================================================
%==============================数据准备==============================
%===================================================================
%%---------------------电力系统数据加载部分-----------------------------------------------
clear
x=0;
Branch=0;%支路参数
Note=0;%节点参数
[filename,pathname]=uigetfile('*.xls','pleasechoosetheexcelfilewithyourpowersystemparameters');%从外部excel导入电力系统潮流计算相关参数
try
iffilename~=0
x=xlsread([pathname,filename],'sheet1','A3:
F3');
Branch=xlsread([pathname,filename],'sheet1','A5:
G10');%读支路参数
Note=xlsread([pathname,filename],'sheet1','A15:
H19');%读节点参数
end
catch
%进行出错处理
errmsg=lasterr;
errordlg(errmsg,'SaveasError');
rethrow(lasterror);
end
%%---------------------支路参数初始化部分-----------------------------------------------
SB=100;UB=220;n=1;m=1;pr=0.0001;
SB=x(5);%功率基准值
UB=x(6);%电压基准值
n=x
(1);%节点数
nl=x
(2);%支路数
m=x(3);%PQ节点的个数
pr=x(4);;%误差精度
B1(:
1)=Branch(:
1);%1、支路首端号
B1(:
2)=Branch(:
2);%2、末端号
B1(:
3)=Branch(:
3)+Branch(:
4)*i;%3、支路阻抗
B1(:
4)=Branch(:
5)*i;%4、支路对地电纳
B1(:
5)=Branch(:
6);%5、支路的变比K:
1;
B1(:
6)=Branch(:
7);%6、支路首端处于K侧为1,1侧为0'
%%
%%---------------------节点参数初始化部分--------------------------------------------------
U=ones(n,1);
a=zeros(n,1);
Ps=zeros(n,1);
Qs=zeros(n,1);
P=zeros(n,1);
Q=zeros(n,1);
detp=zeros(n-1,1);
detq=zeros(m,1);
deta=zeros(n-1,1);
detu=zeros(m,1);
k=0;%迭代次数
U=Note(:
2);%各节点电压初始值(标幺值)
a=Note(:
3);%各节点电压相位初始值(弧度)
Gp=Note(:
4);%各节点发电机有功功率初始值(标幺值)
Gq=Note(:
5);%各节点发电机无功功率初始值(标幺值)
Lp=Note(:
6);%各节点负载有功功率初始值(标幺值)
Lq=Note(:
7);%各节点负载无功功率初始值(标幺值)
type=Note(:
8);%节点类型,PQ节点=1,PV节点=2,平衡节点=3
forh=1:
n
Ps(h)=Gp(h)-Lp(h);%各节点注入的有功功率
Qs(h)=Gq(h)-Lq(h);%各节点注入的无功功率
end
%%%---------------------导纳矩阵计算部分-----------------------------------------------------
Y=zeros(n);
forh=1:
nl%支路数
ifB1(h,6)==0%左节点处于低压侧(6、支路首端处于K侧为1,1侧为0)
p=B1(h,1);q=B1(h,2);%1、支路首端号;2、末端号;
Y(p,q)=Y(p,q)-1./B1(h,3);%非对角元3、支路阻抗;4、支路对地电纳;5、支路的变比;
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1./B1(h,3)+B1(h,4);
Y(q,q)=Y(q,q)+1./B1(h,3)+B1(h,4);
else
p=B1(h,1);q=B1(h,2);%1、支路首端号;2、末端号;
Y(p,q)=Y(p,q)-1./(B1(h,3)*B1(h,5));%非对角元3、支路阻抗;4、支路对地电纳;5、支路的变比;
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1./B1(h,3)+B1(h,4);
Y(q,q)=Y(q,q)+1./(B1(h,3)*B1(h,5)^2)+B1(h,4);
end
end
%导纳矩阵显示
disp('导纳矩阵Y=');
disp(Y)
%%%-------------OK,至此潮流计算所需的数据已经准备好了----------------
%===================================================================
%==============================潮流计算==============================
%===================================================================
%u(i)=e(i)+jf(i);Y(ij)=G(ij)+jB(ij);
G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部
%============================计算失配功率初始值detp\detq==========================
forh=1:
n-1
s=0;
forj=1:
n
s=s+U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));
end
P(h)=U(h)*s;
end
forh=1:
n-1
s=0;
forj=1:
n
s=s+U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));
end
Q(h)=U(h)*s;
end
forh=1:
n-1
detp(h)=Ps(h)-P(h);
end
forh=1:
m
detq(h)=Qs(h)-Q(h);
end
%============================不满足精度要求则进入循环==========================
while(max(abs(detp))>=pr|max(abs(detq))>=pr)%%不满足精度要求则循环
%=================================求取Jacobi矩阵===============================
H=zeros(n-1,n-1);
N=zeros(n-1,m);
K=zeros(m,n-1);
L=zeros(m,m);
forh=1:
n-1
forj=1:
n-1
ifh==j
H(h,j)=U(h)^2*B(h,j)+Q(h);
else
H(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));
end
end
end
forh=1:
n-1
forj=1:
m
ifh==j
N(h,j)=-U(h)^2*G(h,j)-P(h);
else
N(h,j)=-U(h)*U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));
end
end
end
forh=1:
m
forj=1:
n-1
ifh==j
K(h,j)=U(h)^2*G(h,j)-P(h);
else
K(h,j)=U(h)*U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));
end
end
end
forh=1:
m
forj=1:
m
ifh==j
L(h,j)=U(h)^2*B(h,j)-Q(h);
else
L(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));
end
end
end
%========================解修正方程,得到修正量detu,deta=====
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 牛拉法 潮流 计算
![提示](https://static.bdocx.com/images/bang_tan.gif)