-0.015 0
2
4
0
6
8
10
-0.04 0
2
4 t/s
6
8
10
t/s
-0.02 0 0.05 2 4 t/s dθdt(rad/s) 6 8 10
dθdt(rad/s)
0
16
-0.05 0
2
4 t/s
6
8
10
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
图形化显示计算结果
分别绘制在两幅图中
Y/m
1000
1500
500
1000 1000
0 0 500 1000 1500
1500
2000
2000
2500 3000
40
X/m
哈尔滨工业大学控制与仿真中心
研究方法
采用数值积分方法求解微分方程初值问题
欧拉法 梯形法 龙格库塔法
30
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
仿真步骤
编写描述运动方程的m函数 编写数值积分方法求解微分方程初值问题的m文件 初始化微分方程的初值，设定积分步长和时间 调用数值积分程序，在程序中调用m函数求解方程 存储计算结果 图形化显示计算结果
0.005
θ(rad)
0
-0.005
24
-0.01 0
5
10
15
20
哈尔滨工业大学控制与仿真中心 25 30
基于m文件的建模与仿真
问题回顾 单摆系统建模与仿真 拦截系统建模与仿真
25
哈尔滨工业大学控制与仿真中心
拦截系统建模与仿真
系统模型
运动方程
dV P cos α − X − mg sin θ m = dt dθ mV = P sin α + L − mg cos θ dt 2 X = ρV C (α , Ma ) S d 2 2 L = ρV C (α ) S y 2 P =0
27
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
系统模型
相对运动方程
dr = dt VT cosηT − V cosη dq 1 (V sin η − V sin η ) = T T dt r q = σ +η q σ T + ηT = y − yM q = arctan( T ) xT − xM
采用数值积分方法求解微分方程初值问题
欧拉法 梯形法 龙格库塔法
7
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
步骤
编写描述运动方程的m函数 编写数值积分方法求解微分方程初值问题的m文件 初始化微分方程的初值，设定积分步长和时间 调用数值积分程序，在程序中调用m函数求解方程 存储计算结果 图形化显示计算结果
h ( K1 + 2 K 2 + 2 K 3 + K 4 ) 6
dx(t ) = f (t , x(t )) dt
h h K 2 = f (t (k ) + , x(k ) + K1 ) 2 2 h h K 3 = f (t (k ) + , x(k ) + K 2 ) 2 2 K 4 = f (t (k ) + h, x(k ) + hK 3 )
、
单摆系统建模与仿真
图形化显示计算结果
线型控制
0.04
θ(rad) dθdt(rad/s)
0.02
0
-0.02
20
-0.04 0
2
4 t/s
6
8
10
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
图形化显示计算结果
绘制在同一幅图中，分栏
0.02
θ/rad
θ(rad)
0
-0.02 0 0.05
基于MATLAB的系统建模与仿真分析
(基于m文件的建模与仿真)
哈尔滨工业大学控制与仿真中心
4.3 基于m文件的建模与仿真
4.3.1 问题回顾 4.3.2 单摆系统建模与仿真 4.3.3 拦截系统建模与仿真
2
哈尔滨工业大学控制与仿真中心
问题回顾
系统建模与仿真问题
问题不存在解析解，或解析解不直观
1 1 − m sin x
2
dx
6
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
考虑阻力
速度 (rad/s) 比例系数 0.001 0.1 0.005 0.11 0.008 0.14 0.01 0.2 0.02 0.5
假定已知初始角度和角速度，怎么求解？
1 x= x1 (0) x10 x 2 kL g x2 (0) = x20 x x x sin x = − 2 2 2 1 m L
L=1.0 L=0.2 L=0.5 L=5.0 L=10.0
23
5 t/s
10
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
基于编写的程序对单摆系统进行分析
分析摆长和初始条件不变情况下，不同质量时的周期变化
质量(kg) 0.2 0.5 1.0 5.0 10.0
更改程序中的m值
0.01 m=1.0 m=0.2 m=0.5 m=5.0 m=10.0
求解步骤 1.建立m函数描述问题 2.编写数值积分算法，对问题进行求解 3.绘制结果图形
3
哈尔滨工业大学控制与仿真中心
基于m文件的建模与仿真
问题回顾 单摆系统建模与仿真 拦截系统建模与仿真
4
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
系统模型
运动方程
1 x= x1 (0) x10 x 2 kL g x2 (0) = x20 x x x sin x = − 2 2 2 1 m L
8
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
描述运动方程的m函数
pendulumequ.m
function dy=pendulumequ(t,x,k,L,m,g) • dy(1) = x(2); • dy(2) = k*L/m*x(2)*abs(x(2)) - g/L*sin(x(1)); • dy = dy'; end
2
4 t/s
6
8
10
dθdt(rad/s)
dθdt(rad/s)
0
-0.05 0
2
4 t/s
6
8
10
21
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
基于编写的程序对单摆系统进行分析
考虑气动阻力系数
速度 (rad/s) 比例系数 0.001 0.1 0.005 0.11 0.008 0.14 0.01 0.2 0.02 0.5
31
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
程序流程
32
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
描述运动方程的程序代码
33
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
编写龙格库塔法积分求解初值问题
龙格库塔法公式
x(k + 1) = x(k ) +
K1 = f (t (k ), x(k ))
34
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
初始化微分方程的初值，设定积分步长和时间
35
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
初始化微分方程的初值，设定积分步长和时间
36
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
初始化微分方程的初值，设定积分步长和时间
37
哈尔滨工业大学控制与仿真中心
哈尔滨工业大学控制与仿真中心
、
拦截系统建模与仿真
系统模型
运动方程
dxT = V= VT cos θT Tx dt dy T = = V VT sin θT Ty dt dV T =0 dt dθ g cos θT aT − + T = dt V VT T
end
12
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
存储计算结果
角度和角速度构成 的矩阵 时间
13
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
最终的代码
14
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
存储计算结果到数据文件
save result
加载数据文件到工作空间
需要编写插值函数（Matlab函数interp1）
22
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
基于编写的程序对单摆系统进行分析
分析质量和初始条件不变情况下，不同摆长时的周期变化
摆长(m) 0.2 0.5 1.0 5.0 10.0
0.01 更改程序中的L值 0.005
θ(rad)
0 -0.005 -0.01 0
变量定义
摆球质量为m 摆杆长度为L 摆角为θ 阻力系数k
x2 =
dθ dt
x1 = θ
5
哈尔滨工业大学控制与仿真中心
、
单摆系统建模与仿真
系统模型
不考虑阻力系数，小角度时
= θ A cos(
角度大时
T =4
g t +ϕ) L
L K (sin 2 θ 0 ) g