5.3 MATLAB 常微分建模与SIMULINK仿真

Solution of van der Pol Equation, µ = 1000 2 1.5 1 0.5 solution y 1 0 -0.5 -1 -1.5 -2 -2.5Biblioteka 0 济南大学控制科学与工程学院
500
1000 1500 计算机仿真技术 time t
2000
2500
3000
9
【例5.12】van der Pol Equation ：参数µ定义为系统的输入量
济南大学控制科学与工程学院 计算机仿真技术
%写成列向量
5
dy1 = y2 dt dy2 = − y1 + 1 − y12 y2 dt
(
)
主程序，用ode45算法求解： tspan=[0 20]； %仿真时间，行向量 y0=[2;0]; %初始值，列向量 [t,y] = ode45(@vdp1,tspan,y0) %得到时间t与变量y，这里，options为默认值 画图表示模拟结果： %分别画出y1与y2随时间变化的动态曲线，注意t与y的行数或列数要相同 plot(t,y(:,1),‘-’,t,y(:,2),‘--’) title('Solution of van der Pol Equation, \mu = 1'); xlabel('time t'); ylabel('solution y'); legend(‘y_1’,‘y_2’); %在图上表示y1与y2标志
求解刚性问题的ODE solver：
• • • • ode15s: 使用Gear方法，也叫向后微分公式，为多步法，如果已知系统 为刚性问题或ode45无法解决问题，可优先采用ode15s。 ode23s：基于 改进Rosenbrock法，一步法 ，当精度要求 不高时，比 ode15s运算更快。当ode15s无法解决问题，可采用ode23s。 ode23t ode23tb
函数程序： function dydt = vdp1000(t,y) dydt1 = y(2); dydt2 =1000* ( 1 - y(1)^2 ) * y(2) - y(1); dydt=[dydt1；dydt2]; 主程序：用ode45算法求解 tspan=[0 3000]; y0=[2;0]; tic [t,y] = ode45(@vdp1000,tspan,y0); toc
济南大学控制科学与工程学院
& x1 = x2 & x2 = −α sin x1 − βx2
Solution of ODE system 2 x1 x2 1.5
1 solution x
0.5
0
0 0.5 计算机仿真技术
-0.5
1
1.5
2 time t
2.5
3
3.5 11
4
5.4 Simulink仿真
+ -
K s ( s + 1)(s + 10)
用SIMULINK模拟 SIMULINK
济南大学控制科学与工程学院
计算机仿真技术
13
[例5.7] 如下闭环系统在不同K值的阶跃响应？K=5 K=10 K=20
+ -
K s ( s + 1)(s + 10)
用SIMULINK模拟 SIMULINK
济南大学控制科学与工程学院
济南大学控制科学与工程学院
计算机仿真技术
10
【 例 5.13】 给 出 如 下 微 分 方 程 模 型 ， α=10 ， β=5 ， 用 MATLAB模拟该方程，初始值为[0 2]，模拟时间t=[0 20]
函数程序： function dx = nl(t,x) dx1 = x(2); dx2 = -10*sin(x(1))-5*x(2); dx=[dx1;dx2]; 主程序： tspan=[0 20]; x0=[0;2]; [t,x] = ode45(@nl,tspan,x0); plot(t,x(:,1),‘-’,t,x(:,2),‘--’) title('Solution of ODE system'); xlabel('time t'); ylabel('solution x'); legend(‘x_1’,‘x_2’);
计算机仿真技术 3
济南大学控制科学与工程学院
ODE求解语法： [t,y] = solver(odefun,tspan,y0,options) • odefun: 要运行的微分方程系统，一般可由函数定义： dydt = odefun(t,y)，t 和 y 为函数的输入量，分别代表时间和系 统变量，dy 则表示对变量的微分值。 tspan：时间向量，表示微分方程求解的时间轴 y0：变量的初始值 options：微分方程求解参数，如精度，算法等要求，一般可不 修改，直接采用默认值。 输出项为[t, y]，为得到的微分方程时间变量与状态变量，可用 于结果输出与绘图。 当函数中需要引用一些参数或输入时，可定义为yp0，相应的， 微分方程系统函数也要引用： dydt = odefun(t, y，yp0)
济南大学控制科学与工程学院
计算机仿真技术
6
Solution of van der Pol Equation, µ = 1 3 y1 2 y2
1 solution y
0
-1
-2
-3
0
2
4
6
8
10 time t
12
14
16
18
20
济南大学控制科学与工程学院
计算机仿真技术
7
【例5.11】van der Pol Equation （stiff 刚性）：刚性方程中，存在 着系统变量在很短时间内变化很大的现象，使用非刚性系统的算 法无法解决或运算时间过长，这时，采用刚性系统的算法。 参数µ由1增加到1000
• Simulink模块库 • 建立一个简单模型：
• 信号源、传递函数模型、示波器图形输出，进行连 线 • 模块命名、颜色设定、旋转… • 设置仿真参数：solve页面，选择仿真时间，选择算 法，设定算法的参数。 • 设置workspace
济南大学控制科学与工程学院
计算机仿真技术
12
[例5.7] 如下闭环系统在不同K值的阶跃响应？K=5 K=10 K=20
dy1 = y2 dt dy2 = − y1 + µ 1 − y12 y2 dt
(
)
首先，写出该系统函数，保存在M-file文件中： function dydt = vdp1(t,y) dydt1 = y(2); dydt2 = ( 1 - y(1)^2 ) * y(2) - y(1); dydt=[dydt1;dydt2];
函数程序： function dydt = vdp_mu(t,y,mu) dydt1 = y(2); dydt2 =mu* ( 1 - y(1)^2 ) * y(2) - y(1); dydt=[dydt1;dydt2]; 主程序： tspan=[0 3000]; y0=[2;0]; mu=100; options=[]; %参数选择为默认 [t,y] = ode15s(@vdp,tspan,y0,[],mu); plot(t,y(:,1)) title('Solution of van der Pol Equation, \mu = 100'); xlabel('time t');ylabel('solution y');
计算机仿真技术
济南大学控制科学与工程学院 授课教师：李实
cse_lis@ 1教1007室
5.3 常微分方程建模仿真
• 5.3.1 ODE解函数 • 5.3.2 常微分方程仿真举例
济南大学控制科学与工程学院
计算机仿真技术
2
5.3.1 ODE解函数
求解非刚性问题的ODE solver：
• • • ode45: 基于显式龙格-库塔法4-5阶公式，为一步法，一般情况下，在大 多数非刚性问题求解中首先采用ode45作为求解办法。 ode23: 基于显式龙格-库塔法2-3阶公式，比ode45低阶，在精度要求不 高或涉及到部分刚性问题时，比ode45运算更快速。 ode113: 使用Adams-Bashforth-Moulton，即线性多步法。需要精度较高 时，比ode45运算速度快，但属于多步法，需要多个之前时间数值。
计算机仿真技术 4
• • • •
[t,y] = solver(odefun,tspan,y0,yp0, options) •
济南大学控制科学与工程学院
5.3.2 常微分方程仿真举例
【例5.10】van der Pol Equation （Non-stiff 非刚性）： 变量在t=0点的初始值为：y1=2，y2=0，参数µ=1 用MATLAB仿真该系统，时间为0到20
计算机仿真技术
14
【例5.12】van der Pol Equation ：参数µ定义为系统的输入量
在SIMULINK中模拟
济南大学控制科学与工程学院
计算机仿真技术
15
济南大学控制科学与工程学院 计算机仿真技术 8
用ode15s算法求解： tic [t,y] = ode15s(@vdp1000,tspan,y0) toc plot(t,y(:,1)) title('Solution of van der Pol Equation, \mu = 1000'); xlabel('time t');ylabel('solution y_1');