系统仿真的MATLAB实现.docx
《系统仿真的MATLAB实现.docx》由会员分享，可在线阅读，更多相关《系统仿真的MATLAB实现.docx（33页珍藏版）》请在冰豆网上搜索。
系统仿真的MATLAB实现
第七章系统仿真的MATLAB实现
由于计算机技术的高速发展，我们可以借助计算机完成系统的数字仿真。
综前所述，数字仿真实质上是根据被研究的真实系统的模型，利用计算机进行实验研究的一种方法。
仿真的主要过程是：
建立模型、仿真运行和分析研究仿真结果。
仿真运行就是借助一定的算法，获得系统的有关信息。
MATLAB是一种面向科学与工程计算的高级语言，它集科学计算、自动控制、信号处理、神经网络和图像处理等学科的处理功能于一体，具有极高的编程效率。
MATLAB是一个高度集成的系统，MATLAB提供的Simulink是一个用来对动态系统进行建模、仿真和分析的软件包，它支持线性和非线性系统，能够在连续时间域、离散时间域或者两者的混合时间域里进行建模，它同样支持具有多种采样速率的系统。
在过去几年里，Simulink已经成为数学和工业应用中对动态系统进行建模时使用得最为广泛的软件包。
MATLAB仿真有两种途径：
（1）MATLAB可以在SIMULINK窗口上进行面向系统结构方框图的系统仿真；
（2）用户可以在MATLAB的COMMAND窗口下，用运行m文件，调用指令和各种用于系统仿真的函数，进行系统仿真。
这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。
下面介绍在MATLAB上实现几类基本仿真。
7.1计算机仿真的步骤
在学习计算机仿真以前，让我们先总结一下计算机仿真的步骤。
计算机仿真，概括地说是一个“建模—实验—分析”的过程，即仿真不单纯是对模型的实验，还包括从建模到实验再到分析的全过程。
因此进行一次完整的计算机仿真应包括以下步骤：
（1）列举并列项目每一项研究都应从说明问题开始，问题由决策者提供或由熟悉问题的分析者提供。
（2）设置目标及完整的项目计划目标表示仿真要回答的问题、系统方案的说明。
项目计划包括人数、研究费用以及每一阶段工作所需时间。
（3）建立模型和收集数据模型和实际系统没有必要一一对应，模型只需描述实际系统的本质或者描述系统中所研究部分的本质。
因此，最好从简单的模型开始，然后进一步建立更复杂的模型。
（4）编制程序和验证利用数学公式、逻辑公式和算法等来表示实际系统的内部状态和输入/输出的
关系。
建模者必须决定是采用通用语言如MATLAB、FORTRAN、C还是专用仿真语言来编制程序。
在本教材中，我们选择的是MATLAB和其动态仿真工具Simulink。
（5）确认
确认指确定模型是否精确地代表实际系统。
它不是一次完成，而是比较模型和实际系统特性的差异，不断对模型进行校正的迭代过程。
（6）实验设计确定仿真的方案、初始化周期的长度、仿真运行的长度以及每次运行的重复次数。
（7）生产性运行和分析
通常用于估计被仿真系统设计的性能量度。
利用理论定性分析、经验定性分析或系统历史数据定量分析来检验模型的正确性，利用灵敏度分析等手段来检验模型的稳定性。
8）文件清单和报表结果
9）实现
图7.1是计算机仿真的程序图
图7.1计算机仿真程序流图
7.2基于数值积分法的连续系统仿真
7.2.1数值积分法的MATLAB实现
MATLAB的工具箱提供了各种数值积分方法函数，这些函数是ODE23、
ODE45、ODE113和ODE15S。
这些函数均是m文件，还有一个函数是odel.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为计算点时间向量，丫为微分方程的解。
函数ode45（）也是用Runge-Kutta法求解微分方程，它是变步长的一种中等阶次积分算法，调用格式为
[T,Y]=ode45（'F',TSPAN,YO,OPTIONS）
各项含义同上。
函数ode113（）是变阶的Adams-Bashforth-Moulton，用变阶方法解微分方程，
采用多步法，调用格式为
[T,Y]=ode113（'F',TSPAN,YO,OPTIONS）
各项含义同上。
函数odel5s（）采用改进的Gear法解微分方程，调用格式为
[T,Y]=odel5s（'F',TSPAN,YO,OPTIONS）
各项含义同上。
MATLAB还提供了解微分方程函数ode23s。
函数ode1.C是用Euler法求解微分方程。
在SIMULINK中，系统仿真有变步距和固定步距两种工作方式。
在变步距仿
真中，可选用ode45ode23ode113ose15s和ode23s在固定步距仿真中，可选用材ode5、ode4ode3ode2和ode1。
ode5是5阶Runge-Kutta法，ode4是4阶Runge-Kutta法,odel是Euler法。
【例7.1】求微分方程x==x5,0_t_10,x0=1。
用MATLAB编写程序：
建立一个m函数文件dfun.m
funciony=dhfn（t,x）
y=sqrt（3*x）+5;
在MATLABCOMMAND窗口下运行
%MATLABPROGRAM7-1
[t,x]=ode23（'dfun',[010],1）;
Plot（t,x）;
7.2.2基于数值积分法的连续系统的数字仿真
对于一个n阶微分方程表示的连续系统，也可以用具有n个状态变量的状态空
间表达式来描述：
X二AXBU
丫二CXDU
状态方程X=AXBU实际上是n个一阶微分方程组成的方程组。
如果应用四
阶Runge-Kutta法进行仿真计算。
则该式可以改写为
Xk^Xk-（K!
2K22K3K4）
6
KAXkBUk
心二A（Xk-KJBU（tk7）
22
K3二A（Xk7K2BU（tk-））
22
K4=A（Xkha）BU（tkh）
式中，Ki，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';
%Initializationyref=5000;
x0=[00];
y0=[000];
u0=0;
t0=0;
tfinal=1;
tsamp=0.01;
h=0.001;
Y1imit=512;
Ilimit=40;
%调入系统模型
%初始化
%Referentvalueofsystemoutput
%Initialvalueofservovalve
%Initialvalueofhydrauliccylinder
%Starttimeofsimulation
%Endtimeofsimulation
%Sampleperiod
%Simulationstepsize
max_epoch=fix（tfinal/h）—1;t=t0;
u1=u0;
x=x0';
y=y0;
%模块valve仿真
%模块hysya仿真
%离散控制量计算
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';
%Themainloopfori=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）;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）;
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+lt=t+h;tout（i）=t;uout（i）=u;yd（i）=Kf*y（l）;yout（i,:
）=y';
%Discretecontrolprocessifabs（round（t/tsamp）-t/tsamp）ul=Kd*ye;
%Saturationblockiful>Ylimit
xl=ylimit;elseiful<－ylimitxl=－Ylimit;elsexl=ul;
end
end%D/Aconversionx2=Kda*xl;
%Amplifiergainx3=Ka*x2;
%V/Iconversionu4=Kvi*x3;
%limitofcurrentifu4>Ilimit
x4=Ilimit;
elseifu4<－Ilimitx4=－Ilimit;
elsex4=u4;end
end
u=Kq*x4;
end
end%savedatatofilehout=[toutuoutyout];savehout.dathout-ascii;plot（tout,yd,'y'）;gird;
%fordiscretesection
%formainloop
%存储仿真数据
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提供以下几种离散化方法，可以选用:
1'zoh'为零阶保持器
2’foh'为一阶保持器
3讪。
n'为双线性变换法，"齐=需
4’prewarp为改进的双线性变换法
5'matched使连续和离散系统具有匹配的DC增益
【例7.3】连续系统传递函数H（s）二十1e^35s，采样一阶采样保持器，采样
s2+4s+5
周期为Ts=0.1s，求其离散化系统模型，并比较离散前后系统阶跃响应
用MATLAB编写程序：
%MATLABPROGRAM7-3
sys=tf（[I-1],[145],'td',0.35）;
sysd=c2d（sysc,0.l,'for'）
step（sysc,sysd）;
运行结果：
Transferfunction:
0.011^30.045®2-0.05622-0.009104zW-1.629八5+0.670ZM
Samplingtime:
0.1
015
启<
4175j■-*I=———
如叫051.62432
Tire（vec.J
图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
必须是一个标量
rand
s=rand（'state'）
返回均匀分布随机数发生器的当前状态矢量，其维数是35。
这种形式一般
不常用
-5353MATLAB5以上的版本使用了多种随机数发生器，它可以产生［2-53，253］内的所有浮点数。
它的序列周期为21492。
例7.4】
>>rand（3,4）ans=
0.9528
0.5982
0.8368
0.3759
0.7041
0.8407
0.5187
0.8986
0.9539
0.4428
0.0222
0.4290
>>a=zero（2,2）
%zeros函数可以起到声明变量的作用
a=
00
00>>rand（size（a））
ans=
0.19960.53850.30310.9201
%产生［1,3］区间上的均匀分布随机数
>>（3-1）*rand（2,3）+l
ans=
2.05061.06892.5374
1.61372.43071.1190
2．randn
其作用是产生标准正态分布的随机数或随机矢量，其调用形式和rand类似
例7.5】
>>randn（1,2）ans=
%产生N（1,0.25）的正态分布
-0.4326-1.6656>>randn（1,2）/2+1ans=
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）;
dv=r>th;p=p+dv;
end
在这里我们对它的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）ans=
64
61
表7.1MATLAB统计工具箱里的随机数发生器
功能说團
betfimd
Beta分布
binomd
:
进制分带
Expmd
捋数分布
Fmd
F分布
gaiTun
側马井帝
gcomd
几何分布
hygcmd
却凡何分布
nhinnid
负二冷廿布
ncfmd
非申心F舟布
nctmd
非申心1号布
iwrmrnd
泊松分布
randjum
一揩定的分布
mylmd
瑞剧廿布
tnid
T分布
unidmd
unifmd
均匀分布
WEibrndI唏怕并布
7.4.2离散事件系统仿真的一个实例一一报童问题仿真
报童问题是一个古典的概率统计分析问题，虽然问题本身并不复杂，但作为一
个演示实例，可以反映出离散事件系统计算机仿真的很多特征。
1•报童问题
一报童从报刊发行中心定报后零售，每卖一份报纸可赚钱a元，若定报后卖不出去，则可再退回发行处，此时每退一份报要赔钱b元。
虽然每天卖出报的份数是随机的，但报童可根据以往卖报情况的统计来获得每天卖k份的概率p（k），试求报
童每天期望受益达到最大的定报量z'。
2•数学摸型
设报童每天订报z份，而报纸每天卖出y份，我们假设y的分布为
p（y=k）=Pkk=0,1,2,
考虑到报童每天的损失有如下两种情形。
（1）供过于求。
因退货造成的平均损失为：
QO
G（z-k）pk
k=0
（2）供不应求。
因缺货造成的平均损失为
QO
C2二a'（k-z）pk
k=z1
所以，每天的期望损失费（也可以从总收益的角度来考虑）为
C（z）=CiC2
现在我们的目标是求出使得每天期望损失最小的定报量，换言之，就是使报童
的每天期望总收益达到最大。
写成一个目标函数的形式
z=Cm?
n（z）
约束条件如z的取值范围，要受到报童的资本多少的影响。
只有在特殊的概率分布情况下，我们才可以推导出C（z）的解析形式，并通过求极值的方法来求解。
但在实际的应用中，这样的思路往往是行不通的。
但是可以通过枚举所有可能的订报量，求出对应的平均损失，进行比较求出满足条件的z,这
里搜索域通常是有限的。
以上就是一个比较简单的计算机仿真方法。
3.报童问题的计算机仿真
对于给定每一订报量Z值，利用离散随机变量采样算法产生给定分布的随机数R,用来表示报童当天卖出的报纸数，从而可以计算出一天的损失以及一个阶段的平均损失。
这里比较关键的一点是如何产生服从给定分布的随机变量，这个内容在
本教材中不再详尽介绍。
而且在实际的应用中，分布并非总是给定了的，需要我们
收集数据，并从中辨识分布，进行参数估计。
图7.5是根据上述思路设计出的仿真
流程框图。
其中各变量含义如下:
S最小损失值
这里，a和b是这个问题的两个参数
4．计算机仿真
根据图7.5所示框图，我们不难写出报童问题的仿真程序，其MATLAB源码如下。
function[superz,supers]=baotong（tm,g,a,b）
z=1;
supers=1000;
whilez%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;
end;
上面的代码，是在随机变量均匀分布的假设下编写的，对于其他的分布，大家只要用相应的随机数发生器进行替换即可
5．仿真结果和输出数据分析
在分布为均匀分布，参数a=0.2，b=0.4的条件下，在MATLAB命令窗口运行仿真程序就可以得出
升级会员 