电力系统计算机辅助分析试验报告.docx
《电力系统计算机辅助分析试验报告.docx》由会员分享，可在线阅读，更多相关《电力系统计算机辅助分析试验报告.docx（49页珍藏版）》请在冰点文库上搜索。
电力系统计算机辅助分析试验报告
附录：
学生实验报告表头格式
电力工程学院学生实验报告
实验课程名称：
《电力系统计算机辅助分析》
开课实验室：
计算中心4082014年12月11日
年级、专业、班
学号
姓名
成绩
实验项目名称
电力系统计算机辅助分析
指导教师
教师评语
教师签名：
年月日
实验一MATLAB软件的基本操作及程序的基本结构
一、程序流程、源码及仿真结果
启动MATLAB软件，识别出常用的四个窗口――命令窗口（CommandWindow）、历史命令窗口（CommandHistory）、工作空间浏览器（WorkspaceBrowser）。
熟悉Matlab工作环境。
2、上机求解下面的电路问题，脚本文件已经给出。
已知某个ＲＣ电路的端电压的表达式为：
区间时，试绘制电压的波形。
提示：
在MATLAB的M文件编辑器中键入以下命令语句，并保存为ex.m：
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
Clear,clc,close
t=0:
0.5:
10;
u=6*exp（-2*t）;
plot（t,u）;title（‘RC电路的电压响应曲线’）;xlabel（‘时间/s’）;ylabel（‘电压/v’）
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5、对于如下电路图，已知
V，
，
，
，
，自己编写程序，画出
和
的波形（四个周期）。
提示：
写出
的向量表达式用于
的计算（写出串联段和并联段的阻抗，利用分压公式），写出
和
的幅值和相角就能写出波形表达式，设定波形的长度（四个周波）和采样频率（两点间隔）后作图。
注意画图时使用标题，横纵坐标及网格线等命令。
3.
%%%%%%%%%%%%%%%%%
clear,clc,close
mUsm=200;%的幅值
aUsm=（pi/180）*（+90）;%的相角
Usm=mUsm*exp（j*aUsm）;%用于计算
w=100*pi;
r1=22;l=0.05;
r2=33;c=0.0001;
z1=r1+j*w*l;
z2=r2/（j*w*c*r2+1）;
Ucm=Usm*z2/（z1+z2）;
mUcm=abs（Ucm）%求幅值
aUcm=angle（Ucm）%求相角
T=4*2*pi/w;%用w表示，T表示四个周波
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,us,t,uc）%注意区分两条曲线的线型和颜色
gridon
title（''）
xlabel（'','fontsize',16）%fontsize’,16表示横坐标的字体大小为16
ylabel（''）
legend（'us','uc'）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
一、讨论
1、MATLAB软件是什么？
有什么特点？
答：
Matlab是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。
用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程（M文件）后再一起运行。
特点：
1.编程效率高，2.用户使用方便，3.扩充能力强，交互性好，4.移植性和开放性很好，5.语句简单，内涵丰富，6.高效方便的矩阵和数组运算，7.方便的绘图功能。
3.如何获得命令或函数的帮助文档？
答：
MATLAB的各个函数，不管是内建函数、M文件函数、还是MEX文件函数等，一般它们都有M文件的使用帮助和函数功能说明，各个工具箱通常情况下也具有一个与工具箱名相间的M文件用来说明工具箱的构成内容等，在MATLAB命令窗口中，可以通过指令来获取这些纯文本的帮助信息。
通常能够起到帮助作用，获取帮助信息的指令有help、lookfor、which、doc、get、type等。
实验二 电力系统计算中常用的数值算法及电力网络的数学模型
1、用因子表法求解电路方程中的节点电压；
A=[1/10+1/15-1/10-1/15-1
-1/101/10+1/50+1/30-1/500
-1/15-1/501/15+1/50+1/300
1000];
B=[00010];
[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时，用因子表进行求解，体会该算法带来的益处；
1电压源改为2V
②电压源改为4V
③电压源改为6V
④电压源改为8V
3、上机验证课本例17.6－17.7的计算结果；
①例17.6
②例17.7
4、设计系统变更时修改导纳矩阵的函数（参考课本P.323），并调试通过，给出导纳矩阵修正函数的调用规则及其主要的实现代码。
%n=input（'请输入节点数:
n='）;
%nl=input（'请输入支路数:
nl='）;
%B=input（'请输入由支路参数形成的矩阵:
B='）;
%X=input（'请输入由节点号及其对地阻抗形成的矩阵:
X='）;
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.20];
%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）;
仿真结果：
X=
10
20
30
40
导纳矩阵Y=:
0.9346-4.2949i-0.4808+2.4038i-0.4539+1.8911i0
-0.4808+2.4038i1.0690-4.7568i-0.5882+2.3529i0
-0.4539+1.8911i-0.5882+2.3529i1.0421-7.5773i0+3.3333i
000+3.3333i0-3.3333i
请输入追加支路参数形成的矩阵:
B1=[340.3i01.20];
新的导纳矩阵Y1=
0.9346-4.2949i-0.4808+2.4038i-0.4539+1.8911i0
-0.4808+2.4038i1.0690-4.7568i-0.5882+2.3529i0
-0.4539+1.8911i-0.5882+2.3529i1.0421-9.0440i0+4.0000i
000+4.0000i0-3.3333i
思考题：
1、高斯消去法与因子表法有何异同？
答：
因子表法是以高斯消去法为基础求解线性方程组的一种方法，当线性方程组的系数矩阵不变，而常数项发生变化时，采用因子表法可避免系数矩阵的重复运算。
因子表法首先要建立因子表，然后用因子表对常数项进行消去、回代，从而求出线性方程组的解。
2、节点导纳阵与节点阻抗阵之间有何关系，二者之间如何转换？
电力系统的网络方程式一般都用节点方程式表示，节点方程有导纳型和阻抗型两种。
节点导纳矩阵具有对称、稀疏、可以根据网络接线图直接写出、易于修改的特点；节点阻抗矩阵是节点导纳矩阵的逆阵，它具有对称、满阵的特点，它不能由网络接线图直接写出，通常由支路追加法形成。
实验三 短路电流计算
合并课本第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（Iazl（i））;
end
disp（'各支路B相电流Ibzl分别为（顺序同你输入B时一样）:
'）;
fori=1:
nl
disp（Ibzl（i））;
end
disp（'各支路C相电流Iczl分别为（顺序同你输入B时一样）:
'）;
fori=1:
nl
disp（
升级会员 