牛拉法潮流计算文档格式.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=====
升级会员 