连续系统的数字仿真1_08解析
第09讲 连续系统的数字仿真
9
ode45：四/五阶龙格－库塔法，它采用的是单步法， 也就是说它在计算y(tn)时，它仅需要最近处理时刻的 结果y(tn-1)。该算法对于大多数系统有效，最常用， 但不适用于刚性（stiff）系统。 discrtet：当Simulink检测到模型中没有连续状态 时所选择的一种算法。
10
定步长(Fixed-step)算法
15
单任务模式(Single
Tasking)
该模式不检查模块间的采样速率转换。该模式 对于建造单任务系统模型非常有用，在此类系统中， 任务同步不是问题。
自动模式(Auto)
当选用此模式时，如果模型中所有模块运行于 同样的采样率下， Simulink 使用单任务模式；如果 模型包含有不同采样速率运行的模块，则使用多任 务模式。
25
(1) 利用示波器模块(Scope)得到输出结果
当利用示波器模块作输出时，它不仅会自动地将 仿真的结果从示波器上实时地显示出来。而且也可同 时把示波器缓冲区存储的数据，送到MATLAB工作 空间指定的变量中保存起来，以便利用绘图命令在 MATLAB命令窗口里绘制出图形。
26
在示波器模块的窗口中，利用快捷按钮“ ”， 可打开如图6-35所示的示波器模块参数(parameters) 对话框。示波器参数对话框中有两个页面，图6-35 （a）为一般参数设置（General），图6-35（b）为 数据存储参数设置（Data history）。
第六讲 连续系统仿真概述祥解
26
今希望用直角坐标输出，故还引入两个定义变量 x y(21)，y)* cos[ y(2)] y(22) y(1)* sin[ y(2)]
根据卫星的发射速度，可以建立起方程组的初始条件： y(1)0=6400km(卫星到地心的距离，即为地球之半径)，
任务失败。
计算机仿真技术
6
Xb＝Vb*t
Yb=100
Xf’=Vf*cosθ
Yf’=Vf*sinθ
tanθ=(100-Yf)/(Xb-Xf) 描述方程－》龙格库塔法求解位置－》求解我机角度
及位置－》判断距离时间－》end
计算机仿真技术
7
二阶系统：汽车轮子悬置系统
y’’(t)+p1y’(t)+y(t)=1.0
Y ( y1 , y2 , yn ) P ( p1, p2 , pk )
U (u1, u2 ,um )
G ( g1 , g2 , gn )
是 n 维状态向量 是 k 维参数向量 是 m 维控制向量 是n 维函数向量 计算机仿真技术
24
Y (0) [ y1 (0), y2 (0), yn (0)] 是给定的初始状态
计算机仿真技术
4
微分方程数学模型
一般形式
dny d n1 y dy d n1u d n2 u a1 n1 a n1 a n y c1 n1 c2 n1 cn u n dt dt dt dt dt
n为系统的阶次，ai为系统的结构参数，ci为输入函数的
16
连续系统数字仿真中离散化最基本的算法是数值积分 f ( y ,u ,t ) 的系统，已知系统变量 算法。对于形如 y 的初始条件 y (t 0 ) y 0 ，现在要求y随时间变化的过 程y(t)。计算过程可以这样考虑：首先求出初始点 y(t 0 ) y0 的 f (t 0 ，y0 ) ，微分方程可以写作：
系统工程之连续系统仿真
b1s n 1 b2 s n 2 bn Y (s) G ( s) n n 1 U ( s ) s a1s an 1s an Y (s) X (s) X (s) U (s) (b1s
n 1
b2 s
n2
1 bn ) n s a1s n 1 an 1s an
Part 1 连续系统模型概述
系统的概念：所谓系统，是由相互联系、相互作 用的若干部分，以一定的结构组成的具有特定功 能的整体。 模型的概念：对系统的特征与变化规律的一种定 量抽象，是人们用以认识事物的一种手段。
•模型是现实系统的一种抽象，是在一定假设条件下对系统 的转化。 •模型中必须包含系统中的主要因素。 •模型中必须反映出主要因素之间的逻辑关系和数学关系。
系统仿真的应用
•对已有系统进行分析时，采用系统仿真时仿真成 为系统分析器。 •对尚未有的系统进行设计时，采用仿真技术考察 其性能是否满足预定要求，这时仿真是系统设计 器。 •在系统运行时，利用仿真模型作为观测器。 •在系统运行前，利用仿真模型作为预测器。 •利用仿真模型作为训练器，训练系统操作人员， 仿真成为训练仿真器。
系统
系统建模
模型
仿真试验
计算机
系统是研究的对象 模型是系统的抽象 仿真是对模型实验
仿真建模
计算机仿真的三个活动 • 系统建模是通过对实际系统的观测和检测，在忽略次要因 素及不可检测变量的基础上，用物理或数学的方法进行描 述，从而获得实际系统的简化近似模型。 • 仿真建模是将系统模型转化为仿真模型的过程，仿真模型 反映了系统模型同计算机之间的关系，它能为计算机所接 受并在其上运行。 • 仿真实验是将系统的仿真模型置于计算机上运行的过程。
系统工程导论 第五章 系统建模与仿真 第五节连续系统仿真与离散系统仿真
每天接待顾客人 次
发生天数
30～39 5
40～49 25
50～59 40
60～69 28
70～79以上 2
每位顾客购货金额 （元）
发生天数
10～19 20～29 30～39 40～49
40
30
15
10
50～59以 上
5
据此可以列出相应的概率分布，今若以随机数01、02、…、98、 99、00来表示上述概率分布
若没有事件发生，则仿真时钟继续等距推进；若有事件发生，则 记录此事件区间，从而可以得到有关事件的时间参数。
若有若干事件同时发生，除了记录该事件的时间参数外，还需事 先规定这种情况下对各类事件处理的优先顺序。
5.5 连续系统仿真与离散系统仿真
2）仿真建模策略
如何建立仿真模型，也就是采用何种方法推进仿真时钟，建立起系 统中各类实体之间的逻辑联系，一般有三种策略，这就是事件调度法、 活动扫描法和进程交互法。
随机数的生成方法大致有如下三种： 1）随机数表(random number table)法 即由人们在事先人为地产生
出一批均匀随机数，并制成表格形式备用。当需要使用它时，直接调 用这张随机数表就可以了。
2）随机数发生器 即在计算机上附加一个能产生随机数的装置，如 附加一个某种放射粒子的发射源装置，由于发射源在单位时间内发射 的粒子数量是随机的，所以计数器记录下来的数值就是随机数了。
① 事件调度法（Event Scheduling）。用事件的观点来分析真实系统， 通过定义事件及每个事件发生时系统状态的变化，按时间顺序确定并执 行每个事件发生时有关的逻辑关系，这就是事件调度法的基本思想。它 直接对事件加以调度。
所有事件连同其发生时间均放在事件表中。模型中有一个时间控制 模块，不断地从事件表中选择具有最早发生时间的事件，推进仿真时钟 到该事件发生的时间，并调用与该事件类型相应的事件处理模块，处理 完后在返回事件控制模块。
连续系统的数字仿真1_08解析
图3-2 义
梯形法的几何意
xk 1
xk
h[ f 2
(tk , xk )
f
(tk1, xk1)]
8
由于上式右边包含未知量xk+1,所以每一步都 必须通过迭代求解，每一步迭代的初值xk+1(0 通常采用欧拉公式来计算，因此梯形法每一步
迭代公式为
xk
(0) 1
xk
hf
(tk , xk )
xk
( 1
一般控制系统的输出动态响应在开始段变化较快， 到最后变化将会很缓慢。这时，计算可以采用变步长 的方法，即在开始阶段步长取得小一些，在最后阶段 取得大一些，这样即可以保证计算的精度，也可以加 快计算的速度。
23
对于一般工程计算，计算精度要求并不 太高，故常用定步长的方法。作为经验 数据，当采用四阶龙格－库塔法作数值 积分计算时，取计算步长
j 1
（3-8）
c1 0, m 1,2,, v
即：
k1 f (tk , xk )
k2 f (tk c2h, xk b21k1h)
k3 f (tk c3h, xk b31k1h b32k2h) 将式（3-7）展 开成h的幂级数并与微分方程式 （3-1）精确解式（3-6）逐项比较,便可求得式
t k 1
x(tk1) x(tk ) f (t, x(t))dt
tk
通常假设离散点t0,t1,…,tn是等距离的，即tk+1tk=h, 称h为计算步长或步距。
5
当t>t0时，x(t)是未知的，因此式（3-2） 右端的积分是求不出的。为了解决这个问题， 把积分间隔取得足够小，使得在tk与tk+1之间 的 f(t,x(t)) 可 以 近 似 看 作 常 数 f(tk,x(tk)) ， 这样便得到用矩形公式积分的近似公式
连续系统仿真概述
第一篇 连续系统仿真篇第一章 连续系统仿真概论按系统模型的特征分类，可以有连续系统仿真及离散事件系统仿真两大类。
本篇讨论连续系统仿真问题。
过程控制系统、调速系统、随动系统等这类系统称作连续系统，它们共同之处是系统状态变化在时间上是连续的，可以用方程式（常微分方程、偏微分方程、差分方程）描述系统模型。
1.1 连续系统模型描述连续系统仿真中的数学模型有很多种，但基本上可分为三类：连续时间模型、离散时间模型及连续－离散混合模型。
本节将对它们的线性定常形式作一介绍，1.2节将介绍几种模型结构变换的方法。
1.1.1 连续时间模型如果一个系统的输入量u(t)，输出量y(t)，系统的内部状态变量x(t)都是时间的连续函数，那么我们可以用连续时间模型来描述它。
系统的连续时间模型通常可以有以下几种表示方式：常微分方程，传递函数，权函数，状态空间描述．本节仅对其一般描述形式作一简要介绍，偏微分方程将在第8章讨论。
1.常微分方程常微分方程可用(1.1)式表示：u c dt ud c dt u d c y a dt dy a dt y d a dt y d n n n n n n n n n nn +++=++++------- 1221111111(1.1) 其中n 为系统的阶次，a i n i (,,,,)=012 为系统的结构参数，),,2,1(n j c j =为输入函数的结构参数，它们均为实常数。
2.传递函数若系统的初始条件为零，即系统在 t=0时已处于一个稳定状态，那么对（1.1）式两边取拉氏变换后可得：)()()()()()()(2211111s U c s U sc s U sc s Y a s sY a s Y s a s Y s n n n n n n n +++=++++----稍加整理，并记：jn j jn n j jjn s a s cs U s Y s G ∑∑=--=-==10)()()( (1.2)(1.2)式称为系统的传递函数。
实验二 面向系统结构图的连续系统数字仿真实验 matlab程序
4.离散相似法编程
时域离散编程求解 程序：
clear num=[8 10]; den=[0.1 1 0 0]; [A,B,C,D]=tf2ss(num,den) sysc=ss(A,B,C,D); T=0.05; sysd=c2d(sysc,T); Ad=sysd.a; Bd=sysd.b;
0 0.5 1
程序：
clear all A = [-1000.25 999.75 0.5;999.75 -1000.25 0.5;0 0 0]; t = 0; x = [1 -1 1]; h = 0.05; M = round(10/h); for k = 1:M tt = t(k) + h; x1 = -exp(-0.5*tt)+exp(-2000*tt)+1; x2 = -exp(-0.5*tt)-exp(-2000*tt)+1; x3 = 0; temp = [x1 x2 x3]; x = [x; temp]; t = [t; tt]; end T = 0.003125; m = 4; t_frog = 0; x_frog = [1 -1 1]; start_time_frog = clock; F = expm(A*T); for k = 0:m tt = 2^(k) * T; temp = F*x_frog(1,:)'; temp = temp'; x_frog = [x_frog; temp]; t_frog = [t_frog; tt]; if k<m F=F*F; end end qT = 2^(m)*T; for k = 1:M-1 tt = t_frog(m+k+1) + qT; temp = F*x_frog(m+k+1,:)'; temp = temp'; x_frog = [x_frog; temp]; t_frog = [t_frog; tt]; end pass_time_frog = etime(clock,start_time_frog); hh = 0.001; end pass_time_rk4=etime(clock,start_time_rk4); t_display = t(2:31,:); x1_display = x(2:31,1); x1_frog_display = x_frog(6:35,1); x1_rk4_display = x_rk4(51:50:1501,1); disp(' 时间 解析解 蛙跳法 RK4法 ') t_rk4 = 0; x_rk4 = [1 -1 1]; start_time_rk4=clock; for k = 1:50*M tt = k*hh; K1 = A*x_rk4(k,:)'; K2 = A*(x_rk4(k,:)'+hh*K1/2); K3 = A*(x_rk4(k,:)'+hh*K2/2); K4 = A*(x_rk4(k,:)'+hh*K3); temp = x_rk4(k,:)'+hh*(K1/2+2*K2+2*K3+K4)/6; temp = temp'; x_rk4 = [x_rk4; temp]; t_rk4 = [t_rk4; tt];
连续系统数字仿真数值积分法
main() { float x11,x21,x31,x10,x20,x30,R; int i; x10=0;x20=0;x30=0;outputy[0][0]=0; ST=10;DT=0.005; LP=ST/DT;VN=1;R=20; for (i=1;i<=LP;i++) { x11=x10-DT*x30+DT*R; x21=DT/0.5*x10+(1-DT/0.5)*x20; x31=DT/0.1*x20+(1-DT/0.1)*x30; outputy[0][i]=x31; x10=x11;x20=x21;x30=x31; } dispcurve(); } 程序为llx56.c
4．2 梯形法 为了提高仿真精度，离散-再现环节采取图4．3 的形式。 (4-11)
X (k 1) X (k )
T [e(k ) e(k 1)] 2
(4-12)
式(4-12)称为梯形公式，其几何解释如图4．4所示。
[例4．1] 已知一多变量系统的结构框图如图4．5所示，请 用梯形公式对此系统进行仿真，并输出y1、y2的仿真结果。
假设线性系统的状态空间描述为
AX BU X Y CX DU
式中：X为n×1维状态向量；U为r×l维输入 向量；A为n×n维系统矩阵；B为n×r维输入 矩阵；Y为m×1维输出向量；C为m×n维输出 矩阵；D为m×r维传递矩阵。
(4-3) (4-4) (4-5) (4-6)
4.3 龙格-库塔(Runge-Kutta)法
在非实时仿真中，有时需要更高的精度。 4．3节和4．4节中将再介绍两种更精确的方 法及其离散-再现环节。
连续系统仿真的方法.
第3章 连续系统仿真的方法3.1 数值积分法连续系统数值积分法，就是利用数值积分方法对广微分方程建立离散化形式的数学模型——差分方程，并求其数值解。
可以想象在数学计算机上构造若干个数字积分器，利用这些数字积分器进行积分运算。
在数字计算机上构造数字积分器的方法就是数值积分法，因而数字机的硬件特点决定了这种积分运算必须是离散和串行的。
把被仿真系统表示成一阶微分方程组或状态方程的形式。
一阶向量微分方程及初值为()(),00t Y Y t Y ⎧⎫⎪⎪⎨⎬⎪⎪⎩⎭Y =F =（3-1）其中，Y 为n 维状态向量，F （t ，Y ）为n 维向量函数。
设方程（3-1）在011,,,,n n t t t t t +=…处的形式上的连续解为()()()()n+1n+1t t n+10t t t =Y t +,(),n Y F t Y dt Y t F t Y dt=+⎰⎰（3-2）设 n =()n Y Y t ，令1n n n Y Y Q +=+（3-3）则有：()1n+1t n Y Y +=也就是说，1(,)n nt n t Q F t Y dt +≈⎰（3-4）如果n Y 准确解()n Y t 为近似值，n Q 是准确积分值的近似值，则式（3-4）就是式（3-2）的近似公式。
换句话说，连续系统的数值解就转化为相邻两个时间点上的数值积分问题。
因此，所谓数值解法，就是寻求初值问题（3-1）的真解在一系列离散点12n t t t <…<…上的近似解12,,,n Y Y Y ……，相邻两个时间离散点的间隔1n n n t t +=-h ，称为计算步距或步长，通常取n =h h 为定值。
可见，数值积分法的主要问题归结为对函数(,)F t y 的数值积分问题，即如何求出该函数定积分的近似解。
为此，首先要把连续变量问题用数值积分方法转化成离散的差分方程的初值问题，然后根据已知的初值条件0y ，逐步地递推计算后续时刻的数值解(1,2,)i y i =…。
第4章连续系统按环节离散化的数字仿真
eA[(k+1)T −τ ]Bu(τ )dτ ∫
kT
令τ=kT+t则上式可得
x[(k +1)T] = eAT x(kT) + ∫ eA(T −τ ) Bu(kT + t)dt
0 T
(4-2)
4
当系统输入u(t)给定时，便可根据式（42）求出系统离散化状态方程的解。由于u(t) 一般为时间t的函数，而且是未知的，故对于 两相邻采样时刻之间的输入u(kT+t)常用以下 两种方法近似处理。 (0 < t < T) 1.令 u(kT +t) ≈ u(kT). 这相当于在系统的输入端加了一个采样开关 和零阶保持器。．
当t=(k+1)T时
uh[(k +1)T] = u(kT) u(kT) − u[(k −1)T] & u[(k +1)T] = T
(4-6)
13
今对典型环节中系数a,b,c,d 的不同情况，求离散 状态变量式输出量的解。 1．当a≠0,b=0(相应有比例、微分和比例微分等环节) 时，由式（4-4）可得
（4-13）
cT cT z[(k +1)] == z(kT) + u(kT) + u[(k −1)T] 2b 2b
同样由式（4-8）和式（4-6）两式可得 d x[(k +1)T] = z[(k +1)T] + u(kT) （4-14） b
17
今将以上三种情况下的典型环节的仿真模 型归纳为一个统一公式
x(kT) = eA(kT ) x(0) + ∫ eA(kT −τ ) B (τ )dτ u
0 kT
x[(k +1)T] = eA[(k+1)T ] x(0) +
第2章 连续系统数字仿真的基本算法
f (t , y)
f (t , y) 在区间 (tk , tk 1 ) 上定积分的近似解
qk ，并且数值方法的共同特点是步进式的（即所谓递推式
先从标量形式开始讨论。考虑一阶微分方程及初值问题
y f (t , y ) y ( t0 ) y 0
(2.5)
解析解为
上一页 下一页 返回
(2.2)
在 t t0 , t1,, tk , tk 1 上的解析解为
y(tk 1 ) y(t0 ) y( t k )
tk 1 t0 tk 1 tk
f ( t , y ) dt f ( t , y ) dt
(2.3)
希望用公式
yk 1 yk qk
返回
tk 1 tk
f (t, y )dt
下一页
2.1.2 欧拉法
欧拉（Euler）法是一种最简单的数值积分算法，对于
y f (t , y ) y ( t0 ) y 0
( 在区间tk , tk 1 )
上求积分，有
tk 1 tk
y(tk 1 ) y(tk )
上一页
下一页
返回
2.1.1 数值积分算法的基本原理
连续系统仿真的数值积分算法就是利用数值积分法将常 微分方程（组）描述的连续系统变换成离散形式的仿真 模型——差分方程（组）。 为了能在计算机上进行求解，首先要把被仿真系统的数 学模型表示为一阶微分方程组或状态空间模型。
上一页
下一页
返回
假设一阶向量微分方程及初值问题为
第2章 连续系统数字仿真的基本算法 2.1 数值积分算法 2.2 数值积分算法的基本分析 2.3 连续系统仿真的离散相似算法 2.4 常用快速数字仿真算法 2.5 实时数字仿真算法 小结
一、连续系统的仿真分析 例： 蹦极跳系统的数学模型
单位延迟模块的参数设置
三、 线性系统仿真分析
例 ：线性离散系统仿真分析 1. 数字滤波器的数学描述 2. 低通数字滤波器的差分方程描述：
y(n) 1.6 y(n 1) 0.7 y(n 2) 0.04u(n) 0.08u(n 1) 0.04u(n 2)
系统的Z变换 Y ( z ) 0.04 0.08z 1 0.04z 2 U ( z) 1 1.6 z 1 0.7 z 2
y(n) 1.6 y(n 1) 0.7 y(n 2) 0.04u(n) 0.08u(n 1) 0.04u(n 2)
此数字滤波器为线性离散系统，使用滤波器形式对其进 行描述如下：
Y ( z ) 0.04 0.08z 1 0.04z 2 U ( z) 1 1.6 z 1 0.7 z 2
通信系统输出信号 原始锯齿波信号
延迟与失真
图5.40 通信系统输出与原始锯齿波信号比较
混合系统设计之二：行驶控制系统 汽车行驶控制系统是应用非常广泛的控制系统之一， 其主要目的是对汽车速度进行合理的控制。系统的工 作原理如下： （1）汽车速度操纵机构的位置发生改变以设置汽车 的速度，这是因为操纵机构的不同位置对应着不同的 速度。 （2）测量汽车的当前速度，并求取它与指定速度的 差值。 （3）由速度差值信号驱动汽车产生相应的牵引力，并 由此牵引力改变汽车的速度直到其速度稳定在指定的 速度为止。 由系统的工作原理来看，汽车行驶控制系统为典型 的反馈控制系统。下面建立此系统的 Simulink 模型并 进行仿真分析。
一、连续系统的仿真分析
例： 蹦极跳系统的数学模型
二、离散系统的仿真分析
例： 人口变化系统的数学模型 这是一个简单的人口变化模型。在此模型中，设某一年的人口数 目为p，其中表示年份n，它与上一年的人口、人口繁殖速率以及新 增资源所能满足的个体数目之间的动力学方程由如下的差分方程所描 述：
面向结构图的连续系统数字仿真
课程设计面向结构图的连续系统数题目字仿真学院计算机科学与信息工程学院专业自动化班级2010级2班学生姓名小指导教师吴诗贤2013 年12 月20 日面向结构图的连续系统数字仿真姓名：陶园班级：10自动化3班学号：2010133330摘要根据自动控制系统中面向结构图的数字仿真的基本思想，探讨了仿真过程中典型环节的规范性、系统的连接矩阵、仿真求解、程序框图问题，并应用到实际的范例当中，并分析了结果总结了相关特点和相关结论。
自动控制系统常常是由许多环节组成的，要应用数字仿真方法对系统进行分析和研究，首先需要求出总的传递函数，再转化为状态空间表达式的形式，然后对其求解。
当改变系统某一环节的参数时，尤其是要改变小闭环中某一环节的参数时，以上整个过程又需要重新计算，这对研究对象参数变化对整个控制系统的影响是十分不便的，为了克服这些缺点，同时大多数从事自动化工作的科技人员更习惯于用结构图的形式来分析和研究控制系统，为此产生了面向结构图的仿真方法。
该方法只需将各个环节的参数及各环节间的连接方式输入计算机，仿真程序就能自动求出闭环系统的状态空间表达式。
本课程设计主要介绍典型环节参数和连接关系构成闭环系统的状态方程的方法，而动态响应的计算，仍采用四阶龙格-库塔法。
这种方法具有便于研究各个环节参数对系统的影响，并可以得到每个环节的动态响应，以及对多输入输出系统的进行仿真的有点。
关键字：结构图；典型环节；连接矩阵；数字仿真；1、设计任务已知某一系统结构如下图所示，编写matlab程序求a分别为2,4,6,8,10,12时输出量y的动态响应。
图12、需求分析及概要设计2.1 需求分析根据上述设计任务我们可以基本明确在我们课程设计当中应该明确以下几个方面：✓熟悉在数字计算机仿真技术中常用的四阶龙格-库塔算法。
✓明确在面向结构图的连续系统数字仿真，典型环节及其系数矩阵确定。
✓明确各连接矩阵的确定。
✓能够熟练运用MATLAB仿真软件。
第3章 连续系统仿真方法
在多数情况下，假设具有同样的积分误差，变步长仿真的时间小于固定步长的仿真时间。
变步长法不适用于实时仿真。
3.2.2.4仿真算法的选择与比较
算法的误差和稳定性
仿真算法主要包括两种误差：截断误差、舍入误差。
截断误差：基于台劳展开公式的数值计算方法都存在截断误差，一般差分公式若局部截断误差为 ，则称有r阶精度，即方法是念。
所谓稳定性是指误差的积累是否受到控制的问题。如果在每步计算过程中，前面积累的舍入误差对实际误差的影响是减弱的，则计算方法是稳定的；反之，则可能由于误差的恶性增长而变得不稳定。如果计算过程发生不稳定情况，计算结果将失去意义。
通常，用试验方程：
来判断积分算法的稳定性。若积分公式为：
•对于低精度问题：可采用欧拉法；
•要提高精度：可减小步长，但增加了计算量和舍入误差。
•因而采用高阶算法，如：4阶R-K法
用数字仿真的方法对微分方程的数值积分是通过某种数值计算方法来实现的。任何一种计算方法只能是原积分的一种近似。
因此，连续系统仿真，从本质上是从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到离散模型来近似原连续模型。
如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。
计算精度和速度是常见的一对矛盾，也是数字仿真重要解决的问题之一。
3.2连续系统仿真算法
3.2.1线性连续系统仿真算法
3.2.1.1线性连续系统数学模型的几种表示方法
⑴高阶微分方程
⑵传递函数
⑶状态方程
⑴、⑵只能给出输入输出的关系，⑶可描述中间的状态。
以上只是对SISO系统。
仿真技术实验程序及思考题解答(仅供参考)
实验一 连续系统的数字仿真一、实验目的1． 熟悉Matlab 中m 文件的编写；2． 掌握龙格－库塔法的基本原理。
二、实验设备计算机、MATLAB 软件三、实验内容假设单变量系统如图所示。
试根据四阶龙格－库塔法，求系统输出y 的动态响应。
1．首先把原系统转化为状态空间表达式：⎪⎩⎪⎨⎧=+=•CXy bu AX X ，根据四阶龙格－库塔公式，可得到： ⎪⎩⎪⎨⎧=++++=+++1143211)22(6k k k k CX y K K K K h X X （1） 其中： ⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=+++=+++=+=)()()2()2()2()2()(3423121h t bu hK X A K h t bu K h X A K h t bu K h X A K t bu AX K k k k k k k k k （2） 根据（1）、（2）式编写仿真程序。
2．在Simulink 环境下重新对上述系统进行仿真，并和1中结果进行比较。
四、实验结果及分析要求给出系统输出响应曲线，并分析计算步长对龙格－库塔法的影响。
计算步长对龙格－库塔法的影响：单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加，不但引起计算量的增大，而且可能导致舍入误差严重积累，因此同积分的数值计算一样，微分方程的解法也有选择步长的问题。
源程序：r=5;numo=[1];deno=[1 4 8 5];numh=1;denh=1;[num,den]=feedback(numo,deno,numh,denh);[A,b,C,d]=tf2ss(num,den);Tf=input('仿真时间 Tf= ');h=input('计算步长 h=');x=[zeros(length(A),1)];y=0;t=0;for i=1:Tf/h;K2=A*(x+h*K1/2)+b*r;K3=A*(x+h*K2/2)+b*r;K4=A*(x+h*K3)+b*r;x=x+h*(K1+2*K2+2*K3+K4)/6;y=[y;C*x];t=[t;t(i)+h];endplot(t,y)Tf=5 h=0.02五、思考题1．试说明四阶龙格－库塔法与计算步长关系，它与欧拉法有何区别。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本章内容
(1) 熟悉在数字计算机仿真技术中常用的几种数值积分法 ，
特别是四阶龙格-库塔法； (2) 典型环节及其系数矩阵的确定； (3) 各连接矩阵的确定； (4) 利用MATLAB在四阶龙格-库塔法的基础上，对以状 态
空间表达式和方框图描述的连续系统进行仿真； (5) 了解以增广矩阵法为基础的连续系统的快速仿真方法1 。
xk 1
xk
hxk
1 2!
h2
xk
hp p!
xk( p)
0(h p31-)5）
可以看出，提高截断误差的阶次，便可提高其精
度，但是由于计算各阶导数相当麻烦，所以直接采
用泰勒级数公式是不适用的，为了解决提高精度问
Hale Waihona Puke 题，龙格和库塔两人先后提出了间接使用泰勒级数
公式的方法，即用函数值f (t,x)的线性组合来代替f
j 1
（3-8）
c1 0, m 1,2,, v
即：
k1 f (tk , xk )
k2 f (tk c2h, xk b21k1h)
k3 f (tk c3h, xk b31k1h b32k2h) 将式（3-7）展 开成h的幂级数并与微分方程式 （3-1）精确解式（3-6）逐项比较,便可求得式
t k 1
x(tk1) x(tk ) f (t, x(t))dt
tk
通常假设离散点t0,t1,…,tn是等距离的，即tk+1tk=h, 称h为计算步长或步距。
5
当t>t0时，x(t)是未知的，因此式（3-2） 右端的积分是求不出的。为了解决这个问题， 把积分间隔取得足够小，使得在tk与tk+1之间 的 f(t,x(t)) 可 以 近 似 看 作 常 数 f(tk,x(tk)) ， 这样便得到用矩形公式积分的近似公式
干个离散点a≤ t0 < t1<… <tn≤b处的近似值
x(t1),x(t2),…,x(tn)。
4
3.1.1 欧拉法
欧拉法又称折线法或矩形法，是最简单也是最 早的一种数值方法
将式（3-1）中的微分方程两边进行积分，得
tk1 dx(t)
t k 1
dt f (t, x(t))dt
tk dt
tk
即
图3-2 义
梯形法的几何意
xk 1
xk
h[ f 2
(tk , xk )
f
(tk1, xk1)]
8
由于上式右边包含未知量xk+1,所以每一步都 必须通过迭代求解，每一步迭代的初值xk+1(0 通常采用欧拉公式来计算，因此梯形法每一步
迭代公式为
xk
(0) 1
xk
hf
(tk , xk )
xk
( 1
或简化为 x(tk1) x(tk ) f (tk , x(tk ))h
这就是欧拉公x式k。1 xk f (tk , xk )h
6
以x(t0)=x0作为初始值，应用欧拉公式， 就可以一步步求出每一时刻tk的xk值，
即 k =0,x1≈x0+f(t0,x0)h
k=1, x2≈x1+f(t1,x1)h x2
┇
x1
k=n-1,
x
(t2 ,x2)
(t1 ,x1)
xn≈xn-1+f(tn-1,xn-1)h 这样式(3-1)的解x(t) 就求出来了。欧拉法的计 算虽然比较简单，但精度 较低。图3-1为欧拉法的 几何解释。
x0
t
t0 t1 t2 图3-1 欧拉法的几何解释
7
3.1.2 梯形法
由上可知欧拉公式中的积分是用矩形面积
xk 1
xk
hfk
h2 2!
(
f
' k
f' xk
fk ) 0(h p1（) 3-6）
为了避免计算式（3-6）中的各阶导数项，可令xk+1由
以下多项式表示。
v
xk1 xk h amkm m1
（3-7）
12
式中 am为待定因子，v为使用f函数值的个数，
km满足下列方程
m1
km f (tk cmh, xk bmjk jh)
hf
(tk , xk )
xk
(1) 1
xk
h[ f 2
(tk , xk )
f
(tk
1
,
xk
(0) 1
)]
通常称这类方法为预估－校正方法。它首
先根据欧拉公式计算出xk+1的预估值
xk+1(0) ，然后再对它进行校正，以得到更
准确的近似值xk+1(1)。
10
3.1.4 龙格－库塔法
根据泰勒级数将式（3-1）在tk+1=tk+h时刻的解 xk+1=x(tk+h) 在tk附近展开，有
（3-7）和式（3-8）中的系数am ,bmj和cm等。
13
现以v=2为例，来说明这些参数的确定方法。
用数字计算机来仿真或模拟一个连续控 制系统的目的就是求解系统的数学模型。 由控制理论知，一个n阶连续系统可以被描 述成由n个积分器组成的模拟结构图。因此 利用数字计算机来进行连续系统的仿真， 从本质上讲就是要在数字计算机上构造出n 个数字积分器，也就是让数字计算机进行n 次数值积分运算。可见，连续系统数字仿 真中的最基本的算法是数值积分算法。
2
3
3.1 数值积分法
连续系统通常把数学模型化为状态空间表达式，为
了对n阶连续系统在数字计算机上仿真及求解，就要 采用数值积分法来求解系统数学模型中的n个一阶微 分方程。
设n阶连续系统由以下n个一阶微分方程组成
dx(t)
dt
f (t, x(t))
（3-1）
x(t0 ) x0
所谓数值积分法，就是要逐个求出区间［a ,b］内若
(t,x)的导数，然后按泰勒公式确定其中的系数, 这
样既能避免计算f (t,x)的导数，又可以提高数值计算
精度，其方法如下。
11
因
xk f (tk , xk ) fk
xk
df
(t, dt
x)
t tk
f f dx
( t
x
) dt
t tk
fk f xk fk
x xk
x xk
故式（3-5）可写成
R1)
xk
h[ 2
f
(tk , xk )
f
(tk
1
,
x(R) k 1
)]
(3-3）
式中 迭代次数R=0,1,2,…
9
3.1.3 预估－校正法
虽然梯形法比欧拉法精确，但是由于每 一步都要进行多次叠代，计算量大，为了 简化计算，有时只对式（3-3）进行一次叠 代就可以了，因此可得
xk
(0) 1
xk
f(tk,xk)h 来近似的。 由图3-2知，用矩形面积
f dx
dt
c
tkabtk+1代替积分，其误差就
是图中阴影部分。为了提高精
a
b
度现用梯形面积tkactk+1来代
替积分，即
tk1
h
tk
t k 1
t
tk
f
(t, x(t))dt
[ 2
f
(tk , xk )
f
(tk1, xk1)]
于是可得梯形法的计算公式为