第六章 控制系统仿真
6.1 状态空间法的系统仿真
6.2 非线性环节的仿真
6.3 离散系统的仿真
6.4 SIMULINK仿真




仿

信号流图
状态空间
单输入单输出系统
( SISO)
主要研究内容
通过系统的数学模型和计算方法,编写程序运算语句,使之能自动求
解各环节变量的动态变化情况,从而得到关于系统输出和所需要的中间各
变量的有关数据、曲线等,以实现对控制系统性能指标的分析与设计。
实现步骤
根据数学模型、要求
的精度和时间,确定
数值计算方法
按算法要求通过分解、
综合、等效变换等方
法转换成适于在计算
机上运行的公式
上机调试并不断改进,
满足系统各项动态性能
指标,并得到理想的仿
真结果
用合适的开发语
言进行算法编程
和实现
多输入多输出系统
( MIMO)
6.1 状态空间法的系统仿真
一、四阶龙格 -库塔( Runge-Kutta)法
1,Runge-Kutta法推导
''ny (,)f x y实际上,可以由对 求导得到
''
'''
22
( ) ( )
2
xy
x x x y x y y y y x y
x x x y x y y y y
y f f f
y f f f f f f f f f f f
f f f f f f f f f
? ? ?
? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ?
高阶微商就更复杂了。
为了提高精度,用 r阶展开计算公式,2' '' ( )
1 2 ! !
r
r
n n n n n
hhy y h y y y
r? ? ? ? ? ?
我们考虑计算不同点上 (,)f x y 的值,由欧拉公式:
用高阶差分代替各级导数这一思想得到各阶 Runge-Kutta法。
定义运算格式为,1 (,)nnk h f x y?
1
1
(,),1,2,
i
i n i n ij j
j
k h f x a h y b k i n
?
?
? ? ? ??
计算不同点上的函数值的线性组合
,,i ij ia b c 为待定参数
由式( 1)用 h幂次重新整理得
1
1
n
n n i i
i
y y c k?
?
?? ?
通过台劳展开得
1 (,)n n n ny y h f x y? ?? ( 1)
' 2 ''1( ) ( ) ( ) ( )
2!n n n ny x h y x h y x h y x? ? ? ? ?
相比较系数,尽可能重合到较高幂次,而求得待定参数。
以 i=2为例,只计算两次函数 (,)f x y 的值。
1 1 1 2 2
1
2 2 2 1 1
(,)
(,)
nn
nn
nn
y y c k c k
k h f x y
k h f x a h y b k
?
? ? ?
?
? ? ?
1 2 2 21,,,c c a b 12,kk (,)nnxy如何求出 四个系数,将 在 点展开
'1 (,) ( )n n nk hf x y hy x??
32 2 2 1 1 2 2 1 1(,) [ (,( ) ) ] ( )n n n n x yk h f x a h y b k h f x y x a h f b k f o h? ? ? ? ? ? ?
,xyff,nnxy为在 处取值

' 2 321
1 1 2 2 2 1
2
( ) ( ) ( ) ( ) ( )n n n x yby y x c c hy x c a h f f k o ha? ? ? ? ? ? ?
如此有:
' 2 ''1( ) ( ) ( ) ( )
2!n n n n
y x h y x h y x h y x? ? ? ?
21
1 2 2 2
2
11,,1
2
bc c c a
a
? ? ? ?
如果取
1 2 2 2 1
1,1
2c c a b? ? ? ?则
,则满足上式有递推公式:
1 1 2
1
21
11
22
(,)
(,)
nn
nn
nn
y y k k
k hf x y
k hf x h y k
?
? ? ?
?
? ? ?
其局部截断误差为 3()oh
如果进而考虑计算两次函数值 (,)f x y,而提高截断误差阶次呢? k2多展开一项:
1
2 2 2 2 3
2 2 2 1 1 2 2 2 1 1 2 1
2
2 2 3 2 42 1 2 1 2 1
22 2
2 2 2
1
[ (,( ) ) ( ) ] ( 2 ) ( )
2!
1
(,( ) ) ( ) ( 2 ( )
2
n n x y x x x y y y
n n x y x x x y y y
k h f x y x a h f b k f a h f a b h k f b k f o h
b b b
h f x y x a h f f f a h f f f f f o h
a a a
? ? ? ? ? ? ?
? ? ? ? ? ? ?
而有
2' 2 2 3 2 4
2 1 2 1 2 11 1 2 2 2 2 2
2
2 2 2
1( ) ( ) ( ) ( ) ( 2 ) ( )
2n n n x y x x x y y y
b b by y x c c h y x a c h f f f c a h f f f f f o h
a a a? ? ? ? ? ? ? ? ? ?但
23' 2 2 4
1( ) ( ) ( ) ( ) ( 2 ) ( )2 ! 3 !n n n x y x x x y y y x y y
hhy x y x h y x f f f f f f f f f f f f o h
? ? ? ? ? ? ? ? ? ? ? ?
注意:
''' ' 2( ) ( ) 2n x y x x x y y y y x y yy x f f f f f f f f f f f f f? ? ? ? ? ? ?
2,x y yf f f f (,)f x y两式中的 难以消掉,因此达不到两次计算函数值
为此必须增加函数值的计算。
来达到 4阶精度,
一般常用是计算 4阶函数值,得到每步截断误差为 5()oh 的四阶 Runge-Kutta法
n 1 1 2 3 4
1
1
2
2
3
43
1
( 2 2 )
6
(,) ;
(,) ;
22
(,) ;
22
(,)
n
nn
nn
nn
nn
y y k k k k
k hf x y
kh
k hf x y
kh
k hf x y
k hf x h y k
?
? ? ? ? ?
?
? ? ?
? ? ?
? ? ?
其递推公式为,
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
???
???
???
?
?
)22(
6
),(
)
2
,
2
(
)
2
,
2
(
),(
43211n
34
23
12
1
kkkk
h
xx
hkxhtfk
k
h
x
h
tfk
k
h
x
h
tfk
xtfk
n
nn
nn
nk
nn
2、根据四阶龙格 -库塔法的递推公式:
已知开环系统的状态方程为
??
?
?
??
Cxy
BuAxx?
采用四阶龙格 -库塔法进行求解和仿真,其求解步骤和方法如下,:
1、由,可知 ;BuAxx ??? BuAxxtf ??),(
2、仿真算法
3,由 时刻的状态为, 得到1?nt 1?nx
11 ?? ? nn Cxy取 不断递推,便可得到所需时刻各点的状态变量 和输出量 。
Nn,,2,1,0 ?? )( ntx )( nty
二、闭环系统的模型建立
典型闭环控制系统的方框图
对 SISO系统 r,u,y,v 均为标量, 由图可知, 得vyru ??
)( vyrBAxx ????
可得到系统的闭环状态方程,BrxABrxB v CAx ?????
b)(?
Cxy ?又由
ode4( )函数实现上述算法,其程序框图如下:
Y
开始
求四阶 R u n g - K u t t a 法各次斜率 k
j
计算
1143211n
),22(
6
x
???
??????
nnn
Cxykkkk
h
x
t =t
f
结束
计算闭环状态空间系数矩阵 A
b
=A - B v C
N
1?
?
nn
xx
ode4( )函数调用格式为:
[t,y]=ode4(A,B,C,D,x0,h,r,v,t0,tf)
其中,{A,B,C,D}为系统的系数矩阵,x0为状态向量初值,
h为仿真步长,r为输入信号的幅值,v为反馈系数,
t0为仿真的起始时间,tf终止时间,
输出值 t为仿真时间,y为输出量。
ode4( )函数
ode45( )函数
MATLAB中的 ode45() 函数可实现四阶 /五阶龙格 -库塔算法,其调用格式为:
[t,y]=ode45('f',tspa,x0)
其中,f为定义的常微分方程函数名,
tspa为起止时间向量,
x0为初始状态向量,
MATLAB仿真程序 l602
例 6-1,已知系统的开环传递函数为:
在零初始条件下,当输入信号的幅值为 1时,试绘制单位负反馈系统的仿真曲线。
22 )125.0(
)12()(
?
??
ss
sksG
MATLAB仿真程序 l601
?令 K=0.1,1,10,绘制单位负反馈系统的仿真曲线。
例 6-2,已知系统的状态方程为:
在零初始条件下, 阶跃信号的输入幅值为 100,试应用状态空间法对系统进行仿真 。
uxx
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
12
4
3
182424
10912
10914
?
?当 K=0.1时,令 v=0.1,1,10,绘制反馈系统的仿真曲线。
6.2 非线性环节的仿真
一、饱和非线性
饱和非线性环节的数学表达式:
11
11
11
r
c r r
r
s u s
u u s u s
s u s
? ? ??
?? ? ? ?
?
? ?
?
uro
-s1
s1
uc
Saturation( )函数
Saturation( )函数调用格式为:
uc=Saturation(ur,s1)
其中,ur为输入向量
s1为非线性环节的饱和值
uc为输出向量
N
开始
ur ≥ s1?
ur>0?
uc=s1 uc=-s1 uc=ur
返回
NY
Y
Saturation( )函数实现上述算法,其程序框图如下:
二、死区非线性
死区非线性环节的数学表达式:
?
?
?
?
?
??
???
???
?
11
11
11
0
susu
sus
susu
u
rr
r
rr
c
O
s1
-s1
uc
ur
DeadZone( )函数
DeadZone( )函数调用格式为:
uc=DeadZone(ur,s1)
其中,ur为输入向量
s1为死区非线性环节的死区值
uc为输出向量
DeadZone( )函数实现上述算法,其程序框图如下:
N
开始
ur ≥s 1?
ur >0?
uc=ur-s1 uc=0
返回
NY
Y
uc=ur+s1
三、间隙非线性
间隙非线性环节的数学表达式:
?
?
?
?
?
?
?
??
??
???
???
?
00
00
00
00
1
1
crcs
crcs
crr
crr
c
uuu
uuu
uusu
uusu
u
??
??
??
??




s1
-s1
uc
urO
backlash( )函数
backlash( )函数调用格式为:
[uc,uss]=backlash(urs,ur,ucs,s1)
其中,ur为输入向量,uc为输出向量
urs,ucs为 ur,uc前一时刻的值
uss为 下次运算保留的输入值
s1为环节的间隙宽度值
Y
Y
N
N
N
开始
ur >urs?
ur-s1≥ ucs?
uc=ur-s1 uc=ur+s1
返回
NY
Y
uc=ucs
ur<urs?
ur+s1≤ ucs?
uss=ur
backlash( )函数实现上述算法,其程序框图如下:
6.3 离散系统的仿真
y(t)e(k)+
-
r(t)
保持器 被控对象
数字
计算机 D/AA/D
数字控制系统方框图
离散系统的数学模型一般用差分方程和离散状态方程来描述。
一, 差分方程法
差分方程描述法系统仿真的步骤如下:
? 根据系统的结构图, 在适当位置加设虚拟采样开关和保持器;
? 将原系统转换成状态空间形式, 并按指定的采样周期, 依照离散化方法, 将系统
离散化, 并得到离散化的状态方程, 即系统的差分方程
? 输入系统初始化参数, 并根据差分方程编写仿真程序 。
?,,,k
D u ( k )C x ( k )y ( k )
H u ( k )G x ( k ))x ( k
210
1
?
?
?
?
??
???
diffstate( )函数实现上述算法,其程序框图如下:
Y
开始

)()()(
)()()1(
kDukCxky
kHukGxkx
??
???
计算结束?
数据输出
结束
输入系统参数 G
,
,H,C,D,
初始值 x 0,y 0,运行参数 T,N
N
1?? kk
diffstate( )函数调用格式为:
[t,xx]=diffstate(G,H,x0,u0,N,T)
其中,{G,H}为差分方程的系统矩阵,
x0为初始状态值,u0输入向量,
N为仿真点数,T为采样周期,t采样时间序列,
xx为状态向量,
diffstate( )函数
MATLAB仿真程序 l603
例 6-3,已知定常离散系统的状态方程为:
给定初始状态为,试求解,并绘制其仿真曲线。
)(11)(116.0 10)1( kukxkx ????????????? ????
????????? 11)0(,1)( xku )(kx
二,Z变换法
+
-
y(t)r(t)
)1(
1
?sss
e Ts??1
例 6-4,已知系统如图所示, 试求零初始条件下系统的开环和闭环离散化状态方
程, 并绘制闭环阶跃响应曲线 。
MATLAB仿真程序 l604
系统的开环离散化状态方程为:
)(0952.0 0048.0)(9048.00 0952.01)1( kukxkx ???????????????
)(0952.0 0048.0)(9048.00952.0 0952.09952.0)1( kukxkx ????????????? ???
系统的闭环离散化状态方程为:
6.4 SIMULINK仿真
一,Simulink的基本操作
? 在 MATLAB的命令窗口直接键入, Simulink”并回车;
? 单击 MATLAB工具条上的 Simulink 图标;
? 在 MATLAB菜单上选 File→New→Model 。
?运行 Simulink 的三种方式:
?常用的标准模块
?模块的操作
? 模块的选取
? 模块的复制、剪切、删除、移动
? 模块的连接
? 模块参数的设置
例 6-5,已知单位负反馈二阶系统的开环传递函数为
试绘制单位阶跃响应的实验结构图 。
sssG 47.4
10)(
2 ??
10
s +4.47s2
Transfer FcnStep Scope
?模块的外形的调整
1,改变模块的大小
2,调整模块的方向
3,给模块加阴影
?模块名的处理
1,模块名的显示与消隐
2,修改模块名
例 6-6:对例 6-5的结构图进行模块处理
?à ?ó ?÷
ê? 2¨?÷μ¥ ?? ?×? ? D? o?
10
s +4.47s2
′? μY oˉ êy
二、系统仿真及参数设置
? 算法设置 ( Solver)
?设置起始时间和终止时间( Simulation time)
?算法设置 (Solver option)
?算法类型设置
?仿真算法设置
?设置输出选项
变步长解法
定步长解法
? 工作空间设置 ( Workspace I/O)
?从工作空间读入数据( Load from workspace )
?保存数据到工作空间 (Save to workspace)
?初始化状态模块
三,simulink仿真分析
sim()函数调用格式为:
[t,x,y]=sim(‘model’,timespan,option,ut)
其中,‘ model’为 Simulink生成的模型文件名;
timespan为仿真时间设置,可指定终止时间和起止时间;
option是用于设置初始条件、步长与容许误差等值,
ut为外部输入信号。
系统默认仿真算法为变步长的四阶 /五阶龙格 -库塔积分算法即 ode45()函

1、带延迟环节 系统 的仿真
例 6-7,已知系统的 Simulink结构图如图所示, 其中:延迟环节的延迟时间为
1秒, 文件名为 exp1.mdl,试绘制系统的仿真曲线 。
1
O u t 1
?ó 3ù ?·? ú
1
s+1
±? ?? ?? ?ó
1
I n 1
MATLAB仿真程序 l607
2,含非线性环节系统的仿真
例 6-8,已知含有饱和非线性环节的系统如图所示, 初始状态为零, 饱和非线性
环节的饱和值 c=0.5,试对系统含有饱和非线性环节前后进行仿真, 并绘制其单位
阶跃响应曲线 。
+
-
y ( t ) u ( t )
c )1(
1
?ss
1
Out1
1
s +s2
Transfer Fcn
1
In1
解,1、绘制系统含有饱和非线性环节前后的 Simulink结构图,如图所示;
含饱和非线性环节系统的 Simulink结构图( exp3.mdl)
2
Out2
1
Out1
1
s +s2
Transfer FcnSaturation
1
In1
原系统的 Simulink结构图( exp2.mdl)
2、运行 MATLAB仿真程序 l608,得到系统含有饱和非线性环节前后的单位阶跃响应曲线
3,离散系统的仿真
例 6-9,系统的 Simulink结构图如图所示, 文件名为 exp4.mdl,试对离散线性系
统进行仿真并求取其频域指标 。
1
O u t 1
0, 7 z+ 0, 0 8
z - 1, 2 z+ 0, 3 72
1
I n 1
解,运行 MATLAB仿真程序 l609,得到离散系统的仿真曲线
综合应用
例 6-10,系统的结构如下图所示, 设输入幅值为 10,间隙非线性的宽度为 1,试对
包含非线性环节前后的系统进行仿真 。
+
-
u(t)
s5.01
1
? s1.01
1
?s
s
101
105
?
?
s
1 y(t)
解,绘制 SIMULINK系统结构图 l610
例 6-11,已知钢铁厂车间加热炉传递函数与温度传感器及其变送器的传递函数
模型分别为:

设定控制所用的 PID调节器传递函数为
试对系统的 PID控制进行分析, 设计与仿真 。
se
sG
80
01 1120
9.9 ?
?? 110
1 0 7.0
02 ?? sG
ss
ssG
c 1 4 55 2 1
5.12 4 09 2 8 6
2
2
?
???
解,绘制 SIMULINK系统结构图 l611
作业
习题 1,4,7