3机9节点系统的临态稳定计算程序设计和算法原理
大家好，今天给大家介绍一下。matlab3机9节点系统的临态稳定计算程序设计与算法原理(论文 程序源码)。
- 3机9节点系统的临态稳定计算程序设计和算法原理
-
- 1、项目简介
- 2、难度指数
- 3.运行环境:
- 4.项目详解:
- 5.下载源码：
1、项目简介
- 本程序是基于隐形梯形积分法的临态稳定分析程序IEEE 3机9节点模型，在pt三相短路故障发生在7号节点，并发生在ct始终模拟故障切除(断开5号节点和7号节点之间的线路)，最后绘制2号和1号发电机的功率角差(δ21)与时间的关系曲线。发电机模型采用三阶模型，忽略阻尼，电枢电阻为零；励磁系统只有测量和放大环节，忽略反馈和电压调节器；恒定阻抗负载。程序是在Matlab R2009b调试和运行。
2、难度指数
项目难度：中难度 适用场景：毕业设计及相关领域的应用研究
3.运行环境:
环境配置： 本项目使用MATLAB版本为Matlab R2009b 注：可适配绝大部分matlab版本 项目技术： 隐形梯形积分法 临时稳定分析 发电机模型 三阶模型 潮流计算等
4.项目详解:
提示：以下是项目的详细介绍，项目源代码和项目设计文件下载地址见文末。
4.1模型推导及相关公式 4.1.1发电机模型 发电机采用三阶模型，其隐形梯形积分法推导如下： 发电机转子运动方程 
式中，=w-1，TJ和t单位为s，wB=314rad/s，其它变量为标准值。根据梯形积分法则，得到tn~tn 1步差分方程为
其中 由发电机转子绕组暂态方程 根据梯形积分法的差异化，得到 式中 发电机定子电压方程 将其转化为xy坐标下为 4.1.2 励磁系统模型 如图1所示，励磁系统传输函数框图：
励磁系统的基本方程是通过图1传递函数
根据梯形积分法，编写相应的差分方程并整理出来
联立式13、10和19，得到系统的差异方程(略下标n 1)
使用前面的公式可以获得公式中的参数。
4.1.3 负荷模型 负载模型被视为恒定阻抗并入导纳矩阵。负载等值并联导纳为
4.1.4 消除网络模型 由于只有发电机端节点注入电流不为零，其他节点注入电流为零，除发电机外的其他节点可以消除网络。消除过程如下：
下标n代表发电机节点，r代表其他节点。
到目前为止，网络差异方程已经完成，即式20和25。N机系统有6个N也有6个未知量N使用牛顿迭代法可以解决个方程。
4.2程序流程图 如下图所示：
故障在t=pt时刻加入，并在t=ct如果系统稳定，始终清除故障t=t 1继续迭代牛顿法，直到系统不稳定，最后绘制发电机功率角差与时间的关系曲线。
4.3变量说明 LN：支路参数。格式如下：
线路状态为0或1。0代表支路为线路，1代表支路为变压器。 GEN：发电机参数。格式如下：
发电机节点类型为0或1。0代表节点PQ节点，1代表节点PV节点。 LOAD：负载参数。格式如下：
ROTOR：发电机暂态参数如下：
EXC：励磁系统参数如下：
sp：常规参数。格式如下：
nsp：非常参数。格式如下：
其它变量定义如下：
4.4输出结果 系统故障在t=pt随时添加，故障是7号节点的三相短路ct始终切断故障，切断5号和7号节点之间的支路。Matlab模拟结果如下： 图是当pt=0.05s，ct=0.167s2号和3号发电机和1号发电机之间的功角差曲线图。从图中可以看出，0~2.4s在此期间，任何功率角差不超过180，系统暂时稳定；但2.4s将来失去稳定。 图是当pt=0.05s，ct=0.168s2号和3号发电机和1号发电机之间的功率角差曲线图。如图所示，系统在2号s内部失去了稳定。
进一步增加故障清除时间，得到下图所示的图形pt=0.05s，ct=0.4s。可见系统在清除故障前已失去稳定性。 可以得出以下结论：当系统发生故障时，随着故障切除时间的延迟，系统的暂态稳定时间范围越来越小。因此，当系统发生故障时，尽快检测和清除故障是确保系统暂态稳定的主要条件。如果采取一些稳定措施，系统将更有可能恢复稳态。
4.5部分程序源代码
function x=nagauss2(a,b) %用途：列主元Gauss线性方程组消除法解ax=b %格式：x=nagauss2(a,b,flag) a为系数矩阵，b为右端列向量,x为解向量 n=length(b);a=[a,b]; for k=1:(n-1) %选主元 [ap,p
]
=
max
(
abs
(
a
(k
:n
,k
)
)
)
; p
=p
+k
-
1
;
if p
>k t
=
a
(k
,
:
)
;
a
(k
,
:
)
=
a
(p
,
:
)
;
a
(p
,
:
)
=t
; end
%消元
a
(
(k
+
1
)
:n
,
(k
+
1
)
:
(n
+
1
)
)
=
a
(
(k
+
1
)
:n
,
(k
+
1
)
:
(n
+
1
)
)
-
a
(
(k
+
1
)
:n
,k
)
/
a
(k
,k
)
*
a
(k
,
(k
+
1
)
:
(n
+
1
)
)
;
a
(
(k
+
1
)
:n
,k
)
=
zeros
(n
-k
,
1
)
; end
%回代 x
=
zeros
(n
,
1
)
;
x
(n
)
=
a
(n
,n
+
1
)
/
a
(n
,n
)
;
for k
=n
-
1
:
-
1
:
1
x
(k
,
:
)
=
(
a
(k
,n
+
1
)
-
a
(k
,
(k
+
1
)
:n
)
*
x
(
(k
+
1
)
:n
)
)
/
a
(k
,k
)
; end function
[nsp
]
=
parametersolve
(sp
,nsp
,EXC
,ROTOR
,U
)
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%求解非定常参数子程序
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%sp
:
[adita
,Pm0
;aw
;aq1
;aEf_1
;aq2
;aEf
;aUTR
;Uref
;Eq0
]
;
%
%
%nsp
:
[IG0
;Eq0_1
;dita0
;Pe0
;UTR0
;Ef0
;Id0
;dw0
;VR
;bdita_1
;bEf_1
;bq
;bUTR
,bEf
]
; aw
=
sp
(
3
,
:
)
;aEf
=
sp
(
7
,
:
)
;Pm0
=
sp
(
2
,
:
)
;aUTR
=
sp
(
8
,
:
)
; Uref
=
sp
(
9
,
:
)
;aq1
=
sp
(
4
,
:
)
;Id
=
nsp
(
7
,
:
)
;Pe
=
nsp
(
4
,
:
)
; UTR
=
nsp
(
5
,
:
)
;dw
=
nsp
(
8
,
:
)
;Ef
=
nsp
(
6
,
:
)
;Eq_1
=
nsp
(
2
,
:
)
; Ef0
=
sp
(
10
,
:
)
;VR
=
nsp
(
9
,
:
)
;dita
=
nsp
(
3
,
:
)
;adita
=
sp
(
1
,
:
)
;
for i
=
1
:
3
bw
(i
)
=
aw
(i
)
*
(
Pm0
(i
)
-
Pe
(i
)
)
+
dw
(i
)
;
bdita
(i
)
=
adita
(i
)
*
dw
(i
)
+
dita
(i
)
;
bdita_1
(i
)
=
(
bw
(i
)
+
aw
(i
)
*
Pm0
(i
)
)
*
adita
(i
)
+
bdita
(i
)
;
bUTR
(i
)
=
aUTR
(i
)
*
(
-
2
*
UTR
(i
)
+
abs
(
U
(i
)
)
)
+
UTR
(i
)
;
bEf
(i
)
=
aEf
(i
)
*
(
2
*
EXC
(i
,
3
)
*
Uref
(i
)
-
2
*
VR
(i
)
-
EXC
(i
,
3
)
*
UTR
(i
)
)
+
VR
(i
)
;
bEf_1
(i
)
=
bEf
(i
)
-
aEf
(i
)
*
EXC
(i
,
3
)
*
bUTR
(i
)
+
Ef0
(i
)
;
bq
(i
)
=
aq1
(i
)
*
(
Ef
(i
)
-
2
*
Eq_1
(i
)
-
(
ROTOR
(i
,
3
)
-
ROTOR
(i
,
4
)
)
*
Id
(i
)
)
+
Eq_1
(i
)
; end
nsp
(
10
,
:
)
=bdita_1
;
nsp
(
11
,
:
)
=bEf_1
;
nsp
(
12
,
:
)
=bq
;
nsp
(
13
,
:
)
=bUTR
;
nsp
(
14
,
:
)
=bEf
; function
[J
,F
]
=
Jform
(sp
,U
,nsp
,ROTOR
,Y
)
%
%
%形成雅可比矩阵子程序
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%sp
:
[adita
,Pm0
;aw
;aq1
;aEf_1
;aq2
;aEf
;aUTR
;Uref
;Eq0
]
;
%
%nsp
:
[IG0
;Eq0_1
;dita0
;Pe0
;UTR0
;Ef0
;Id0
;dw0
;VR
;bdita_1
;bEf_1
;bq
;bUTR
]
; I
=
nsp
(
1
,
:
)
;dita
=
nsp
(
3
,
:
)
;aw
=
sp
(
3
,
:
)
;aq1
=
sp
(
4
,
:
)
;adita
=
sp
(
1
,
:
)
;bq
=
nsp
(
12
,
:
)
; aq2
=
sp
(
6
,
:
)
;aEf_1
=
sp
(
5
,
:
)
;Eq_1
=
nsp
(
2
,
:
)
;bdita_1
=
nsp
(
10
,
:
)
;bEf_1
=
nsp
(
11
,
:
)
; Ix
=
real
(I
)
;Iy
=
imag
(I
)
;Xd_1
=
ROTOR
(
:
,
4
)
;Xq
=
ROTOR
(
:
,
8
)
; Ux
=
real
(U
)
;Uy
=
imag
(U
)
;unit
=
1
;a
=
[
0
,
0
,
0
]
;b
=
[
0
,
0
,
0
]
;J
=
zeros
(
18
,
18
)
;F
=
zeros
(
18
,
1
)
;
for i
=
1
:
3
J
(
6
*
(i
-
1
)
+
1
,
6
*
(i
-
1
)
+
1
)
=unit
;
J
(
6
*
(i
-
1
)
+
1
,
6
*
(i
-
1
)
+
3
)
=
adita
(i
)
*
aw
(i
)
*
Ix
(i
)
;
J
(
6
*
(i
-
1
)
+
1
,
6
*
(i
-
1
)
+
4
)
=
adita
(i
)
*
aw
(i
)
*
Iy
(i
)
;
J
(
6
*
(i
-
1
)
+
1
,
6
*
(i
-
1
)
+
5
)
=
adita
(i
)
*
aw
(i
)
*
Ux
(i
)
;
J
(
6
*
(i
-
1
)
+
1
,
6
*
(i
-
1
)
+
6
)
=
adita
(i
)
*
aw
(i
)
*
Uy
(i
)
;
J
(
6
*
(i
-
1
)
+
2
,
6
*
(i
-
1
)
+
1
)
=
-
aq2
(i
)
*
(
Ix
(i
)
*
cos
(
dita
(i
)
)
+
Iy
(i
)
*
sin
(
dita
(i
)
)
)
;
J
(
6
*
(i
-
1
)
+
2
,
6
*
(i
-
1
)
+
2
)
=unit
;
J
(
6
*
(i
-
1
)
+
2
,
6
*
(i
-
1
)
+
3
)
=
-
aq1
(i
)
*
aEf_1
(i
)
*
Ux
(i
)
/
abs
(
U
(i
)
)
;
J
(
6
*
(i
-
1
)
+
2
,
6
*
(i
-
1
)
+
4
)
=
-
aq1
(i
)
*
aEf_1
(i
)
*
Uy
(i
)
/
abs
(
U
(i
)
)