牛拉法潮流计算文档格式.docx
《牛拉法潮流计算文档格式.docx》由会员分享,可在线阅读,更多相关《牛拉法潮流计算文档格式.docx(16页珍藏版)》请在冰豆网上搜索。
1;
6、支路首端处于K侧为1,1侧为0'
%%电力系统潮流计算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],'
A5:
G10'
%读支路参数
Note=xlsread([pathname,filename],'
A15:
H19'
%读节点参数
end
catch
%进行出错处理
errmsg=lasterr;
errordlg(errmsg,'
SaveasError'
rethrow(lasterror);
%%---------------------支路参数初始化部分-----------------------------------------------
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、支路首端号
2)=Branch(:
2);
%2、末端号
3)=Branch(:
3)+Branch(:
4)*i;
%3、支路阻抗
4)=Branch(:
5)*i;
%4、支路对地电纳
5)=Branch(:
6);
%5、支路的变比K:
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(:
%各节点电压初始值(标幺值)
a=Note(:
3);
%各节点电压相位初始值(弧度)
Gp=Note(:
4);
%各节点发电机有功功率初始值(标幺值)
Gq=Note(:
5);
%各节点发电机无功功率初始值(标幺值)
Lp=Note(:
%各节点负载有功功率初始值(标幺值)
Lq=Note(:
%各节点负载无功功率初始值(标幺值)
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);
nl%支路数
ifB1(h,6)==0%左节点处于低压侧(6、支路首端处于K侧为1,1侧为0)
p=B1(h,1);
q=B1(h,2);
%1、支路首端号;
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
Y(p,q)=Y(p,q)-1./(B1(h,3)*B1(h,5));
%非对角元3、支路阻抗;
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==========================
n-1
s=0;
forj=1:
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;
s=s+U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));
Q(h)=U(h)*s;
detp(h)=Ps(h)-P(h);
m
detq(h)=Qs(h)-Q(h);
%============================不满足精度要求则进入循环==========================
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:
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)));
N(h,j)=-U(h)^2*G(h,j)-P(h);
N(h,j)=-U(h)*U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));
K(h,j)=U(h)^2*G(h,j)-P(h);
K(h,j)=U(h)*U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));
L(h,j)=U(h)^2*B(h,j)-Q(h);
L(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));
%========================解修正方程,得到修正量detu,deta=====