昆明理工大学 电力工程 学院学生实验报告.docx
《昆明理工大学 电力工程 学院学生实验报告.docx》由会员分享，可在线阅读，更多相关《昆明理工大学 电力工程 学院学生实验报告.docx（55页珍藏版）》请在冰豆网上搜索。
昆明理工大学电力工程学院学生实验报告
昆明理工大学电力工程学院学生实验报告
实验课程名称：
《电力系统计算机辅助分析》
开课实验室：
计算机中心4082014年12月9日
年级、专业、班
电自114
学号
姓名
成绩
实验项目名称
电力系统计算机辅助分析实验
指导教师
谢静
教师签名：
年月日
实验一 MATLAB软件的基本操作及程序的基本结构
实验内容：
1.验证实验
2.对于如下电路图，已知
V，
，
，
，
，自己编写程序，画出
和
的波形（四个周期）
程序：
clear
mUsm=200;
aUsm=（pi/180）*（+90）;
Usm=mUsm*exp（j*aUsm）;
R1=22,L=0.05,R2=33,C=0.0001;
w=100*pi;
Z1=R1+（w*0.05*j）;
Z2=（1/（1/R2+j*w*C））;
Ucm=Z2/（Z1+Z2）*Usm;
mUcm=abs（Ucm）;
aUcm=angle（Ucm）;
T=4*2*pi/w;
N=100;
dt=T/N;
t=0:
dt:
T;
fork=1:
101
us（k）=mUsm*cos（w*t（k）+aUsm）;
uc（k）=mUcm*cos（w*t（k）+aUcm）;
end
plot（t,uc,'r--'）;
holdon
plot（t,us）;
title（'电压uc，电压us'）;
xlabel（'时间/s','fontsize',16）;
ylabel（'电压'）;
legend（'uc','us'）;
gridon
波形
思考题：
1．matlab软件是什么？
有什么特点？
答：
MATLAB是美国MathWorks公司自20世纪80年代中期推出的数学软件，具有优秀的数值计算能力和卓越的数据可视化能力，可以提供与矩阵有关的强大的数据处理和图形显示功能。
特点有：
1.编程效率高，2.用户使用方便，3.扩充能力强，交互性好，4.移植性和开放性很好，5.语句简单，内涵丰富，6.高效方便的矩阵和数组运算，7.方便的绘图功能。
3.如何获得命令或函数的帮助文档？
答：
MATLAB的各个函数，不管是内建函数、M文件函数、还是MEX文件函数等，一般它们都有M文件的使用帮助和函数功能说明，各个工具箱通常情况下也具有一个与工具箱名相间的M文件用来说明工具箱的构成内容等，在MATLAB命令窗口中，可以通过指令来获取这些纯文本的帮助信息。
通常能够起到帮助作用，获取帮助信息的指令有help、lookfor、which、doc、get、type等。
实验二 电力系统计算中常用的数值算法及电力网络的数学模型
1.用因子表法求解电路方程中的节点电压；
程序：
clear
clc
close
%U1=10
%-U1/10+U2*（1/10+1/30+1/50）-U3/50=0
%-U1/15-U2/50+U3*（1/15+1/50+1/30）=0
U1=2;
A=[100;-1/101/10+1/30+1/50-1/50;-1/15-1/501/15+1/50+1/30]
B=[U1;0;0]
[n,m]=size（A）;
fori=1:
n
A（i,i）=1/A（i,i）;
forj=i+1:
n
A（i,j）=A（i,j）*A（i,i）;
end
fork=i+1:
n
forj=i+1:
n
A（k,j）=A（k,j）-A（k,i）*A（i,j）;
end
end
end
disp（'矩阵A的因子表为：
'）;
disp（A）
fori=1:
n
B（i）=B（i）*A（i,i）;
forj=i+1:
n
B（j）=B（j）-A（j,i）*B（i）;
end
end
fori=n-1:
-1:
1
forj=i+1:
-1:
2
B（j-1）=B（j-1）-A（j-1,i+1）*B（i+1）;
end
end
disp（'在因子表的基础上求解线性方程组的解为：
x='）;
disp（B）
2、将图2.1中的电压源分别改为2V、4V、6V、8V时，用因子表进行求解，体会该算法带来的益处；
3上机验证课本例17.6－17.7的计算结果
①例17.6
clear
clc
close
n=5;nl=5;
B=[120.03i01.050;
230.08+0.3i0.5i10;
340.015i01.051;
250.1+0.35i010;
350.04+0.25i0.5i10];
X=[10;20;30;40;50];
Y=zeros（n）;
fori=1:
n
ifX（i,2）~=0;
p=X（i,1）;Y（p,p）=1./X（i,2）;
end
end
fori=1:
nl
ifB（i,6）==0
p=B（i,1）;q=B（i,2）;
else
p=B（i,2）;q=B（i,1）;
end
Y（p,q）=Y（p,q）-1./（B（i,3）*B（i,5））;
Y（q,p）=Y（p,q）;
Y（q,q）=Y（q,q）+1./（B（i,3）*B（i,5）^2）+B（i,4）./2;
Y（p,p）=Y（p,p）+1./B（i,3）+B（i,4）./2;
end
disp（'导纳矩阵Y=:
'）;disp（Y）
②例17.7
n=3;
nl=6;
B=[012i010;024i010;122i010;
0320i010;238i010;135i010];
m=0;Z=zeros（n）;
fork1=1:
nl
p=B（k1,1）;q=B（k1,2）;
ifB（k1,6）==0
l=1./B（k1,5）;
elsel=B（k1,5）;
end
ifp==0
ifq>m%追加接地树支
Z（q,q）=B（k1,3）;m=m+1;
else%追加接地连支
fork=1:
m,
Z（k,m+1）=-Z（k,q）;Z（m+1,k）=-Z（q,k）;
end
Z（m+1,m+1）=Z（q,q）+B（k1,3）;
forl1=1:
m
fork=1:
m
Z（l1,k）=Z（l1,k）-Z（l1,m+1）*Z（m+1,k）./Z（m+1,m+1）;
end
Z（l1,m+1）=0;
end
fork=1:
m+1
Z（m+1,k）=0;
end
end
elseifq>m%追加不接地树支
fork=1:
m
Z（k,q）=Z（k,p）*l;Z（q,k）=Z（p,k）*l;
end
Z（q,q）=l^2*Z（p,p）+l^2*B（k1,3）;
m=m+1;
else
fork=1:
m%追加不接地连支
Z（k,m+1）=l*Z（k,p）-Z（k,q）;
Z（m+1,k）=l*Z（p,k）-Z（q,k）;
end
Z（m+1,m+1）=l^2*Z（p,p）+Z（q,q）
-2*l*Z（p,q）+l^2*B（k1,3）;
forl1=1:
m
fork=1:
m
Z（l1,k）=Z（l1,k）-Z（l1,m+1）*Z（m+1,k）./Z（m+1,m+1）;
end
Z（l1,m+1）=0;
end
fork=1:
m+1
Z（m+1,k）=0;
end
end
end
end
disp（'阻抗矩阵Z='）;
disp（Z）
4.设计系统变更时修改导纳矩阵的函数（参考课本P.323），并调试通过，给出导纳矩阵修正函数的调用规则及其主要的实现代码。
n=4;nl=4;B=[120.08+0.4i010;230.1+0.4i010;
340.3i010;130.12+0.5i010];
X=[10;20;30;40]
Y=zeros（n）;
fori=1:
n
ifX（i,2）~=0;
p=X（i,1）;
Y（p,p）=1./X（i,2）;
end
end
fori=1:
nl
ifB（i,6）==0
p=B（i,1）;q=B（i,2）;
else
p=B（i,2）;q=B（i,1）;
end
Y（p,q）=Y（p,q）-1./（B（i,3）*B（i,5））;
Y（q,p）=Y（p,q）;
Y（q,q）=Y（q,q）+1./（B（i,3）*B（i,5）^2）+B（i,4）./2;
Y（p,p）=Y（p,p）+1./B（i,3）+B（i,4）./2;
end
disp（'导纳矩阵Y=:
'）;
disp（Y）
B1=input（'请输入追加支路参数形成的矩阵:
B1='）;
%B1=[340.3i01.10];
%B1=[120.3i010];
%B1=[140.3i010];
%B1=[350.3i010];
z=B1（3）;
%ifB1（6）==0
%p=B1
（1）;q=B1
（2）;
%else
%p=B1
（2）;q=B1
（1）;
%end
p=B1
（1）;q=B1
（2）;
ifB1（6）==1
B1（5）=1/（B1（5））;
end
ifB1
（2）>n
Y（p,p）=Y（p,p）+1/z;Y（p,q）=-1/z;Y（q,p）=-1/z;Y（q,q）=（1/z）+B1（4）;
end
ifY（p,q）==0
Y（p,p）=Y（p,p）+1/z;Y（q,q）=Y（q,q）+1/z;Y（p,q）=-1/z;Y（q,p）=-1/z;
end
ifY（p,q）~=0&&B1
（2）<=n
fori=1:
nl
ifB（i,1）==B1
（1）&&B（i,2）==B1
（2）
ifB（i,5）==B1（5）
Y（p,p）=Y（p,p）+1/z;Y（q,q）=Y（q,q）+1/z;
Y（p,q）=Y（p,q）-1/z;Y（q,p）=Y（q,p）-1/z;
else
K1=B1（5）;K=B（i,5）;
Y（p,p）=Y（p,p）+（K1^2-K^2）*（1/z）;Y（q,q）=Y（q,q）;Y（p,q）=Y（p,q）-（K1-K）*（1/z）;Y（q,p）=Y（q,p）-（K1-K）*（1/z）;
end
end
end
end
ifB1
（2）<=n
n=n;
Y1=zeros（n）;
fori=1:
n
forj=1:
n
Y1（i,j）=Y（i,j）;
end
end
end
ifB1
（2）>n
n=n+1;
Y1=zeros（n）;
fori=1:
n
forj=1:
n
Y1（i,j）=Y（i,j）;
end
end
end
disp（'新的导纳矩阵Y1='）;
disp（Y1）;
思考题：
1.高斯消去法与因子表法有何异同？
答；因子表法是高斯消去法的另一种变化形式其不同之处在于（以按行消去过程为例）：
因子表是对系数矩阵和对常数项的消去及规格化分开写，用下三角及对角元素可对常数项进行消去运算，并利用上三角元素则可进行回代运算。
高斯消去是通过增广矩阵回代过程即可求出方程组的全部解。
2．网络矩阵的修改函数对电力系统计算有何作用
答：
不需要每次去修改因子表，减少了计算量，使计算更加方便
实验三 短路电流计算
1.合并课本第18章例程。
即要求实现计算指定节点各种故障条件下，各支路的故障电流分布情况。
其中对于实用算法，可只提供计算电抗。
提示：
尽量使用先前实验中已封装好的函数。
2.用新程序验证第18章的各例题。
（注意短路电流曲线绘制时，至少包含5~7个点）
W=3;
switchW
case1
clc
clearall
NF=1
fori1=1:
NF
clear;
n1=3
n2=3
n0=3
nl1=5
nl2=5
nl0=5
Lf=1%单相接地短路为1横向故障为0
f=3
If=0
zf=0
zg=0
B1=[010.15i010;020.075i010;120.1i010;130.1i010;230.1i010]
B2=[010.15i010;020.075i010;120.1i010;130.1i010;230.1i010]
B0=[010.05i010;020.025i010;120.2i010;130.2i010;230.2i010]
A1=[111;1-1./2-i*sqrt（3）/2-1./2+i*sqrt（3）/2;1-1./2+i*sqrt（3）/2-1./2-i*sqrt（3）/2]
Vcs=[1;1;1]
Z1=zeros（n1）;Z2=zeros（n2）;Z0=zeros（n0）;Y1=zeros（n1）;Y2=zeros（n2）;Y0=zeros（n0）;V1=zeros（n1,1）;V2=zeros（n2,1）;
V0=zeros（n0,1）;I1=zeros（nl1,1）;I2=zeros（nl2,1）;I0=zeros（nl0,1）;
form1=1:
3
m=0;
ifm1==1
nl=nl1;B=B1;Z=Z1;n=n1;
elseifm1==2
nl=nl2;B=B2;Z=Z2;n=n2;
elseifm1==3
nl=nl0;B=B0;Z=Z0;n=n0;
end
m=0;Z=zeros（n）;V=zeros（n）;I=zeros（nl）;
fork1=1:
nl
p=B（k1,1）;q=B（k1,2）;
ifB（k1,6）==0
k=1./B（k1,5）;
elsek=B（k1,5）;
end
ifp==0
ifq>m%追加接地树支
Z（q,q）=B（k1,3）;m=m+1;
else%追加接地连支
fori=1:
m
Z（i,m+1）=-Z（i,q）;Z（m+1,i）=-Z（q,i）;
end
Z（m+1,m+1）=Z（q,q）+B（k1,3）;
fori=1:
m
forj=1:
m
Z（i,j）=Z（i,j）-Z（i,m+1）*Z（m+1,j）./Z（m+1,m+1）;
end
Z（i,m+1）=0
end
fori=1:
m+1
Z（m+1,i）=0
end
end
elseifq>m%追加不接地树支
fori=1:
m
Z（i,q）=Z（i,p）*k;Z（q,i）=Z（p,i）*k;
end
Z（q,q）=k^2*Z（p,p）+k^2*B（k1,3）;m=m+1;
else
fori=1:
m%追加不接地连支
Z（i,m+1）=k*Z（i,p）-Z（i,q）;Z（m+1,i）=k*Z（p,i）-Z（q,i）;
end
Z（m+1,m+1）=k^2*Z（p,p）+Z（q,q）-2*k*Z（p,q）+k^2*B（k1,3）;
fori=1:
m
forj=1:
m
Z（i,j）=Z（i,j）-Z（i,m+1）*Z（m+1,j）./Z（m+1,m+1）;
end
Z（i,m+1）=0;
end
fori=1:
m+1
Z（m+1,i）=0;
end
end
end
end
ifm1==1
Z1=Z;
elseifm1==2
Z2=Z;
elseifm1==3
Z0=Z;
end
end%以上程序求出正,负,零序对应的阻抗阵
ifIf==0%求出横向故障的阻抗阵
Z1（f,f）=Z1（f,f）;Z2（f,f）=Z2（f,f）;Z0（f,f）=Z0（f,f）;
else%求出纵向故障的阻抗阵
Z1（f,f）=Z1（f,f）+Z1（1,1）-2*Z1（f,1）;Z2（f,f）=Z2（f,f）+Z2（1,1）-2*Z2（f,1）;Z0（f,f）=Z0（f,f）+Z0（1,1）-2*Z0（f,1）;
end%-----------------故障点自阻抗--------------------------------
ifLf==1%求出正负零序短路节点的电流标幺值
Z4=Z2（f,f）+Z0（f,f）+3*zf;k2=1;;k0=1;I1（f）=Vcs（f）./（Z1（f,f）+Z4）;I2（f）=k2*I1（f）;I0（f）=k0*I1（f）;a1=I1（f）;a2=I2（f）;a0=I0（f）;
elseifLf==2%两相短路
Z4=Z2（f,f）+2*zf;k2=-1;k0=0;
I1（f）=Vcs（f）./（Z1（f,f）+Z4）;I2（f）=k2*I1（f）;
I0（f）=k0*I1（f）;a1=I1（f）;
a2=I2（f）;a0=I0（f）;
elseifLf==3%两相接地短路
Z4=zf+（Z2（f,f）+zf）*（Z0（f,f）+zf+3*zg）./（Z2（f,f）+Z0（f,f）+2*zf+3*zg）;k2=-（Z0（f,f）+zf+3*zg）./（Z2（f,f）+Z0（f,f）+2*zf+3*zg）;
k0=-（Z2（f,f）+zf）./（Z2（f,f）+Z0（f,f）+2*zf+3*zg）;I1（f）=Vcs（f）./（Z1（f,f）+Z4）;I2（f）=k2*I1（f）;I0（f）=k0*I1（f）;a1=I1（f）;
a2=I2（f）;a0=I0（f）;
elseifLf==4%单相断线
Z4=zf+（Z2（f,f）+zf）*（Z0（f,f）+zf）./（Z2（f,f）+Z0（f,f）+2*zf）;k2=-（Z0（f,f）+zf）./（Z2（f,f）+Z0（f,f）+2*zf）;
k0=-（Z2（f,f）+zf）./（Z2（f,f）+Z0（f,f）+2*zf）;I1（f）=Vcs（f）./（Z1（f,f）+Z4）;I2（f）=k2*I1（f）;I0（f）=k0*I1（f）;a1=I1（f）;a2=I2（f）;a0=I0（f）;
elseifLf==5%两相断线
Z4=Z2（f,f）+Z0（f,f）+3*zf;k2=1;k0=1;I1（f）=Vcs（f）./（Z1（f,f）+Z4）;I2（f）=k2*I1（f）;I0（f）=k0*I1（f）;a1=I1（f）;a2=I2（f）;a0=I0（f）;
end
disp（'故障点处的正序电流标幺值为:
'）;
disp（a1）;
disp（'故障点处的负序电流标幺值为:
'）;
disp（a2）;
disp（'故障点处的零序电流标幺值为:
'）;
disp（a0）;
fori=1:
n1%求出各节点的电压标幺值
V1（i）=Vcs（i）-I1（f）*Z1（i,f）;%正序
end
fori=1:
n2
V2（i）=-I2（f）*Z2（i,f）;%负序
end
fori=1:
n0
V0（i）=-I0（f）*Z0（i,f）;%零序
end
disp（'各节点的正序电压标幺值（节点号从小到大排）:
'）;
disp（V1）;
disp（'各节点的负序电压标幺值（节点号从小到大排）:
'）;
disp（V2）;
disp（'各节点的零序电压标幺值（节点号从小到大排）:
'）;
disp（V0）;
fori=1:
3%求各支路的各序电流标幺值
ifi==1
B=B1;V=V1;nl=nl1;
elseifi==2
B=B2;V=V2;nl=nl2;
elseifi==3
B=B0;V=V0;nl=nl0;
end
forj=1:
nl
ifB（j,6）==0
k=B（j,5）;
elsek=1./B（j,5）;
end
p=B（j,1）;q=B（j,2）;
ifp==0
e=0;b=B（j,3）;
I（j）=（e-V（q）./k）./b;%接地支路电流
else
I（j）=（V（p）-V（q）./k）./B（j,3）;%不接地支路电流
end
end
ifi==1
I1=I;
elseifi==2
I2=I;
elseifi==3
I0=I;
end
end
disp（'各支路的正序电流标幺值为（顺序同您输入B时一样）:
'）;
fori=1:
nl
disp（I1（i））;
end
disp（'各支路的负序电流标幺值为（顺序同您输入B时一样）:
'）;
fori=1:
nl
disp（I2（i））;
end
disp（'各支路的零序电流标幺值为（顺序同您输入B时一样）:
'）;
fori=1:
nl
disp（I0（i））;
end
fori=1:
3%求故障点A,B,C三相的总电流,各个节点A,B,C三相各相的总电压,各支路A,B,C三相各相的总电流
ifi==1%故障点电流
s0=a0;s1=a1;s2=a2;
elseifi==2%各序支路电压
s0=V0;s1=V1;s2=V2;
else%各序支路电流
s0=I0;s1=I1;s2=I2;
end
forj=1:
3%不对称系数
CA=s0*A1（j,1）+s1*A1（j,2）+s2*A1（j,3）;
ifj==1
A=CA;
elseifj==2
B=CA;
elseifj==3
C=CA;
end
end
ifi==1
Iadl=A;Ibdl=B;Icdl=C;
elseifi==2
Vajy=A;Vbjy=B;Vcjy=C;
else
Iazl=A;Ibzl=B;Iczl=C;
end
end
disp（'故障点处A相电流Iadl为:
'）;disp（Iadl）;disp（'故障点处B相电流Ibdl为:
'）;disp（Ibdl）;
disp（'故障点处C相电流Icdl为:
'）;disp（Icdl）;disp（'各节点处A相电压Vajy分别为（节点号从小到大排）:
'）;disp（Vajy）;
disp（'各节点处B相电压Vbjy分别为（节点号从小到大排）:
'）disp（Vbjy）;
disp（'各节点处C相电压Vcjy分别为（节点号从小到大排）:
'）;disp（Vcjy）;
disp（'各支路A相电流Iazl分别为（顺序同你输入B时一样）:
'）;
fori=1:
nl
disp（I
升级会员 