matlab系统仿真设计Word文档格式.docx
《matlab系统仿真设计Word文档格式.docx》由会员分享，可在线阅读，更多相关《matlab系统仿真设计Word文档格式.docx（26页珍藏版）》请在冰豆网上搜索。
图7.1计算机仿真程序流图
7.2基于数值积分法的连续系统仿真
7.2.1数值积分法的MATLAB实现
MATLAB的工具箱提供了各种数值积分方法函数，这些函数是ODE23、ODE45、ODE113和ODE15s。
这些函数均是m文件，还有一个函数是ode1.C，是直接用C语言编写的。
函数ode23（）是用Runge-Kutta法求解微分方程。
它是一种采用三阶积分算法、二阶误差估计、变积分步长的低阶积分算法，调用格式为
[T,Y]=ode23（'
F'
TSPAN,YO,OPTIONS）
其中，F为系统模型文件名，模型为y'
=f（t,y）形式；
TSPAN=[ToTFINAL]为积分计算时间，初值为To，终值为TFINAL；
YO为系统输出初始值；
OPTIONS选项积分计算相对允差'
RelTol'
和绝对允差'
AbsTol'
，当缺省时，
Reltol＝1e-3,AbsTol＝1e-6
T为计算点时间向量，Y为微分方程的解。
函数ode45（）也是用Runge-Kutta法求解微分方程，它是变步长的一种中等阶次积分算法，调用格式为
[T,Y]=ode45（'
TSPAN,YO,OPTIONS）
各项含义同上。
函数ode113（）是变阶的Adams-Bashforth-Moulton，用变阶方法解微分方程，采用多步法，调用格式为
[T,Y]=ode113（'
函数odel5s（）采用改进的Gear法解微分方程，调用格式为
[T,Y]=odel5s（'
MATLAB还提供了解微分方程函数ode23s。
函数ode1.C是用Euler法求解微分方程。
在SIMULINK中，系统仿真有变步距和固定步距两种工作方式。
在变步距仿真中，可选用ode45、ode23、ode113、ose15s和ode23s。
在固定步距仿真中，可选用材ode5、ode4、ode3、ode2和ode1。
ode5是5阶Runge-Kutta法，ode4是4阶Runge-Kutta法，odel是Euler法。
【例7.1】求微分方程
。
用MATLAB编写程序：
建立一个m函数文件dfun.m
funciony＝dhfn（t，x）
y＝sqrt（3*x）＋5;
在MATLABCOMMAND窗口下运行
%MATLABPROGRAM7-1
[t,x]＝ode23（'
dfun'
，[0l0]，1）;
Plot（t，x）;
7.2.2基于数值积分法的连续系统的数字仿真
对于一个n阶微分方程表示的连续系统，也可以用具有n个状态变量的状态空间表达式来描述：
状态方程
实际上是n个一阶微分方程组成的方程组。
如果应用四阶Runge-Kutta法进行仿真计算。
则该式可以改写为
式中，K1，K2，K3，K4均为n维列向量。
系统输出为
yk＝CXk＋DU（tk）
类似地，可将各种数值积分方法用于连续系统的仿真中。
【例7.2】用MATLAB编写图7.2所示液压控制系统的仿真程序。
图7.2系统方框图
用MATLAB编写仿真程序，采用四阶Runge-Kutta法：
%MATLABPROGRAM7-2
%******Inputsystemdata*****%调入数据文件
hynat;
%Inputsystemfunction;
%调入系统模型
ypfun1='
valve'
;
ypfun='
hysys'
%Initialization%初始化
yref=5000;
%Referentvalueofsystemoutput
x0=[00];
%Initialvalueofservovalve
y0=[000];
%Initialvalueofhydrauliccylinder
u0=0;
t0=0;
%Starttimeofsimulation
tfinal=1;
%Endtimeofsimulation
tsamp=0.01;
%Sampleperiod
h=0.001;
%Simulationstepsize
Y1imit=512;
Ilimit=40;
max_epoch=fix（tfinal/h）－1;
t=t0;
u1=u0;
x=x0'
y=y0;
tout=zeros（max_epoch,1）;
uout=zeros（max_epoch,1）;
yd=zeros（max_epoch,1）;
yout=zeros（max_epoch,legth（y））;
i=1;
tout（i）=t;
uout（i）=1;
yd（i）=Kf*y（l）;
yout（i,:
）=y'
%Themainloop
fori=1:
max_epoch
%Computeoutputofvalve
svl=feval（ypfun1,t,u,x,av,bv）;
sv2=fcval（ypfun1,t+h/2,u,x+h*sv1/2,av,hv）;
%模块valve仿真
sv3=feval（ypfun1,t+h/2,u,x+h*sv2/2,av,hv）;
sv4=feval（ypfun1,t+h,u,x+h*sv3,av,hv）;
x=x+h*（svl+2*sv2+2*sv3+sv4）/6;
vo=cv*x;
%Computeoutputofcylinder
xl=feval（ypfun,t,vo,y,a,b）;
s2=feval（ypfun,t+h/2,vo,y+h*s1/2,a,h）;
%模块hysya仿真
s3=feval（ypfun,t+h/2,vo,y+h*s2/2,a,b）;
s4=feval（ypfun,t+hu,vo,y+h*s3,a,b）;
y=y+h*（sl+2*s2+2*s3+s4）/6;
i=i+l
t=t+h;
tout（i）=t;
uout（i）=u;
yd（i）=Kf*y（l）;
yout（i,:
）=y'
%Discretecontrolprocess%离散控制量计算
ifabs（round（t/tsamp）-t/tsamp）<
le－9
ye=yref－y（l）*Kf;
ul=Kd*ye;
%Saturationblock
iful>
Ylimit
xl=ylimit;
elseiful<
－ylimit
xl=－Ylimit;
elsexl=ul;
end
%D/Aconversion
x2=Kda*xl;
%Amplifiergain
x3=Ka*x2;
%V/Iconversion
u4=Kvi*x3;
%limitofcurrent
ifu4>
Ilimit
x4=Ilimit;
elseifu4<
－Ilimit
x4=－Ilimit;
elsex4=u4;
end
u=Kq*x4;
end%fordiscretesection
end%formainloop
%savedatatofile%存储仿真数据
hout=[toutuoutyout];
savehout.dathout-ascii;
plot（tout,yd,'
y'
）;
gird;
m函数文件
valve.m
functionxd=valve（t,u,x,av,bv）
xd=av*x+hv*u;
m函数文件
hysys.M
functionyd=hysys（t,u,y,a,b）
yd=a*y+h*u
上例的仿真结果如图7.3所示。
图7.3例7.2系统仿真结果
7.3基于离散相似法的连续系统仿真
前面讨论的连续系统仿真采用的方法是数值积分法。
对于连续系统仿真还有一种方法是离散相似法。
所谓离散相似法是首先将连续系统模型离散化，得到等价的或相似的离散化的模型，然后对相似的离散模型进行仿真计算。
根据这一原理，首先应将连续时间系统模型转换为等价的离散时间系统模型。
连续系统离散化处理是通过采样保持器来实现的。
从一个连续时间系统转换为相应的离散时间模型，两者的等价特性取决于采样保持器。
MATLAB提供了将连续系统模型转换为离散时间系统模型的函数C2D，调用格式为
sysd=c2d（sys,Ts）
其中，sys为线性连续时间系统；
Ts为采样时间；
sysd为等价的离散时间系统。
当采样保持器为零阶保持器时；
sysd=c2d（sys,Ts,method）
其中，method为离散化方法，MATLAB提供以下几种离散化方法，可以选用：
①'
zoh'
为零阶保持器
②'
foh'
为一阶保持器
③'
tustion'
为双线性变换法，
④'
prewarp'
为改进的双线性变换法
⑤'
matched'
使连续和离散系统具有匹配的DC增益
【例7.3】连续系统传递函数
，采样一阶采样保持器，采样周期为
，求其离散化系统模型，并比较离散前后系统阶跃响应。
%MATLABPROGRAM7-3
sys=tf（[l-1],[145],'
td'
0.35）;
sysd=c2d（sysc,0.l,'
for'
）
step（sysc,sysd）;
运行结果：
Transferfunction:
Samplingtime:
0.1
图7.4例7.3离散前后系统阶跃响应比较
7.4离散事件系统仿真的MATLAB实现应用实例
7.4.1随机变量的产生
MATLAB提供了两个基本的函数用于产生随机数。
它们是rand和randn。
1．rand
它的作用是产生（0,1）间均匀分布的随机数或向量。
其调用形式有
Y=rand（n）返回一个n维的U（0,l）随机数方阵，n必须是一个标量
Y=rand（m,n）或Y=rand（[mn]）返回一个m×
n维的U（0,l）随机数矩阵，n必须是一个标量
Y=rand（size（A））返回和矩阵A相同维数的矩阵
rand返回单个的随机数
s=rand（'
state'
）返回均匀分布随机数发生器的当前状态矢量，其维数是35。
这种形式一般不常用。
MATLAB5以上的版本使用了多种随机数发生器，它可以产生[2-53，253]内的所有浮点数。
它的序列周期为21492。
【例7.4】
>
rand（3,4）
ans=
0.95280.59820.83680.3759
0.70410.84070.51870.8986
0.95390.44280.02220.4290
a=zero（2,2）%zeros函数可以起到声明变量的作用
a=
00
10
rand（size（a））
0.19960.5385
0.30310.9201
（3-1）*rand（2,3）+l%产生[1,3]区间上的均匀分布随机数
2.05061.06892.5374
1.61372.43071.1190
2．randn
其作用是产生标准正态分布的随机数或随机矢量，其调用形式和rand类似。
【例7.5】
randn（1,2）
-0.4326-1.6656
randn（1,2）/2+1%产生N（1,0.25）的正态分布。
l.06271.1438
3．randperm
其作用是产生一个随机置换，调用形式
p=randperm（n）产生1:
n的一个随机置换。
【例7.6】
p=randperm（4）
p=
2341
以rand和randn为基础可以产生给定分布的随机数，例如泊松分布和二项分布。
M文件函数posong演示了如何产生泊松分布的随机数。
functionp=posong（a,n,m）
w=ones（n,m）*（-a）;
th=exp（w）;
r=ones（n,m）
dv=r>
th;
p=zeros（n,m）;
whilennz（dv）>
0
r=r.*（rand（n,m）.*dv）;
p=p+dv;
在这里我们对它的MATLAB的实现，简要地说明一下。
参数a指泊松分布的参数，n,m则指返回的随机矢量的维数。
这个程序中有几个地方需要注意：
（1）MATLAB支持矢量的逻辑运算，运算结果是一个元素为0或1的矢量。
如“r>
th”这条语句。
（2）函数nnz的作用是统计一个矢量非0元素的个数。
在这里的作用是判断是不是r中所有的元素小于exp（-a）。
（3）运算符“.*”表示数组元素对元素的乘法，如[a1,a2,…，an].*[b1,b2,…,bn]=[al*bl,a2*b2,…,an*bn]
总之这个程序比较充分的利用了MATLAB基于向量的特点，产生二项分布随机数M文件函数，请大家自己尝试。
其实，像上面的这些工作可以由MATLAB完成，MATLAB的统计工具箱提供了一个各种分布随机数的发生器函数。
安装MATLAB的统计工具箱，就可以使用这些函数。
表7.1列出了这些函数。
你可以通过help命令在MATLAB界面上来获得这些函数的详尽使用方法。
下面作为示例来看看random的用法。
4．Random
它的作用是产生一个由name参数指定的分布的随机数。
其调用形式为：
R=random（name,A,M,N）,
name指定分布的形式，A说明该分布的参数，M，N规定产生的随机数的矩阵维数，M为行数，N为列数，它们的缺省值都是1（MATLAB所有函数涉及维数的参数，顺序都是行在前，列在后）。
由于name指定的分布可能不止由一个参数来表示，所以调用形式要做相应的变化。
R=random（name,A,B,M,N），指定的分布含有两个参数。
R=random（name,A,B,C,M,N），指定的分布含有三个参数。
例如，泊松分布只含有一个参数，就要选用第一种形式来调用。
random（'
poiss'
4,2,2）
64
61
表7.1MATLAB统计工具箱里的随机数发生器
7.4.2离散事件系统仿真的一个实例——报童问题仿真
报童问题是一个古典的概率统计分析问题，虽然问题本身并不复杂，但作为一个演示实例，可以反映出离散事件系统计算机仿真的很多特征。
1．报童问题
一报童从报刊发行中心定报后零售，每卖一份报纸可赚钱a元，若定报后卖不出去，则可再退回发行处，此时每退一份报要赔钱b元。
虽然每天卖出报的份数是随机的，但报童可根据以往卖报情况的统计来获得每天卖k份的概率p（k），试求报童每天期望受益达到最大的定报量z'
2．数学摸型
设报童每天订报z份，而报纸每天卖出y份，我们假设y的分布为
考虑到报童每天的损失有如下两种情形。
（1）供过于求。
因退货造成的平均损失为：
（2）供不应求。
因缺货造成的平均损失为
所以，每天的期望损失费（也可以从总收益的角度来考虑）为
现在我们的目标是求出使得每天期望损失最小的定报量，换言之，就是使报童的每天期望总收益达到最大。
写成一个目标函数的形式
约束条件如z的取值范围，要受到报童的资本多少的影响。
只有在特殊的概率分布情况下，我们才可以推导出C（z）的解析形式，并通过求极值的方法来求解。
但在实际的应用中，这样的思路往往是行不通的。
但是可以通过枚举所有可能的订报量，求出对应的平均损失，进行比较求出满足条件的z，这里搜索域通常是有限的。
以上就是一个比较简单的计算机仿真方法。
3．报童问题的计算机仿真
对于给定每一订报量Z值，利用离散随机变量采样算法产生给定分布的随机数R，用来表示报童当天卖出的报纸数，从而可以计算出一天的损失以及一个阶段的平均损失。
这里比较关键的一点是如何产生服从给定分布的随机变量，这个内容在本教材中不再详尽介绍。
而且在实际的应用中，分布并非总是给定了的，需要我们收集数据，并从中辨识分布，进行参数估计。
图7.5是根据上述思路设计出的仿真流程框图。
其中各变量含义如下：
Tm——一轮试验的预定模拟天数T——一轮实验的仿真天数累计值
Z——订报量Z′——最优订报量
G——定报量Z之上界S1——损失值之累计值
S——最小损失值
这里，a和b是这个问题的两个参数。
图7.5报童问题计算机仿真流程框图
4．计算机仿真
根据图7.5所示框图，我们不难写出报童问题的仿真程序，其MATLAB源码如下。
function[superz,supers]=baotong（tm,g,a,b）
z=1;
supers=1000;
whilez<
g
%r=poisondis（5,1,tm）;
r=round（2*randn（1,tm）+5）;
%产生均匀分布随机数
t=1;
s=0;
dv=z>
r;
s=sum（（（z-r）*b）.*dv）;
s=s+sum（（（r-z）*a）.*（l-dv））;
aver_s=s/tm;
ifsupers>
=aver_s
supers=aver_s;
superz=z;
end;
z=z+1;
上面的代码，是在随机变量均匀分布的假设下编写的，对于其他的分布，大家只要用相应的随机数发生器进行替换即可。
5．仿真结果和输出数据分析
在分布为均匀分布，参数a=0.2，b=0.4的条件下，在MATLAB命令窗口运行仿真程序就可以得出仿真的结果。
例如
[z,s]=baotong（5,10,0.2,0.4）
z=
4
s=
0.2400
此外，我们还可以改变参数a、b的值，观察相应最优值的变化，用MATLAB可以很方便的画出C-a，C-b的变化曲线。
6．报童问题模拟系统的推广和应用
报童仿真问题的改进推广和应用可以有以下几个方面：
（1）求每天的卖出报数服从任意分布的情况下，使报童收益最大意义下的最优订报量Z′；
（2）对报纸总发行量进行测算；
（3）适当修改仿真系统，可将其用于其它系统，例如企业的订货和库存策略研究。
7.5SIMULINK动态仿真
SIMULINK是MATLAB重要软件包，用于对动态系统仿真，它适用于连续系统和离散系统，也适用线性系统和非线性系统。
它采用系统模块直观地描述系统典型环节，因此可十分方便地建立系统模型而不需要花较多时间编程。
正由于这些特点，SIMULINK广泛流行，被认为是最受欢迎的仿真软件。
本节对SIMULINK和它的MATLAB实现，作一个简单介绍，为大家今后运用SIMULINK打下基础。
SIMULINK实际上是面向结构的系统仿真软件。
利用它进行系统仿真的步骤如下：
（1）启动SIMULINK，进入SIMULINK窗口；
（2）在SIMULINK窗口下，借助SIMULINK模块库，创建系统框图模型并调整模块参数；
（3）设置仿真参数后，启
升级会员 