ch12 数字仿真实验

1、实验十二 数字仿真实验38612.1数字仿真的基本概念38612.1.1 所谓数字仿真38612.1.2 数字仿真的分类38712.1.3 数字仿真的原理38712.2 数字仿真的方法38812.2.1 解决计算机模拟的一般步骤：38812.2.2蒙特卡罗法38812.2.3离散型随机变量的模拟39212.2.4连续型随机变量的模拟40012.3 系统仿真软件SIMULINK40412.3.1 SIMULINK简介40412.3.2 SIMULINK建模40512.3.3 SIMULINK仿真407实验任务408411实验十二 数字仿真实验学习目的了解计算机模拟仿真在解决实际问题时的重要性，理

2、解使用计算机仿真的原理及方法。掌握计算机模拟仿真的三大类方法：蒙特卡罗法、连续系统模拟和离散事件系统模拟。会利用MATLAB进行模拟仿真实验。实验指导在实际问题中，我们通常会面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，这样进行处理过后的模型与我们面临的实际问题可能相差很远，以致求解得到答案根本无法应用，这时，计算机模拟几乎成为唯一的选择。12.1数字仿真的基本概念我们将一些具有特定的功能、相互之间有一定的内在联系的对象所组成的总体称为一个系统。模拟就是利用物理的、数学的模型类比、模仿现实系统及其演变过程，以寻求变化规律的一种方法。模拟的基本思想是建立一个实验的模型，这个

3、模型包含所研究系统的主要特点，这样做的目的就是通过对这个实验模型的运行，获得所要研究系统的必要信息。12.1.1 所谓数字仿真数字仿真是利用计算机对一个实际系统的结构和行为进行动态演示，以评价或预测该系统的行为效果。它是解决较复杂的实际问题的一条有效途径。数字仿真的一般过程是先对研究对象进行数学描述，然后构建仿真算法，最后在计算机中编程并运行实现。对比于物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难的诸多缺陷，数字模拟花费较少、周期较短，且易修改、可重用，有更强的系统适应能力。但是计算机模拟也有缺陷，比如系统数学模型不易建立、程序调试复杂等。数字仿真可以用于研制产品或

4、设计系统的全过程中，包括方案论证、技术指标确定、设计分析、生产制造、试验测试、维护训练、故障处理等各个阶段。具体来说，下列情形适用于数字仿真：（1）在一个实际系统还没有建立起来之前，要对系统的行为或结果进行分析研究时，数字仿真是一种行之有效的方法。（2）在有些真实系统上做实验会影响系统的正常运行，这时进行计算机模拟就是为了避免给实际系统带来不必要的损失。如在生产中任意改变工艺参数可能会导致废品，在经济活动中随意将一个决策付诸行动可能会引起经济混乱。（3）当人是系统的一部分时，他的行为往往会影响实验的效果，这时运用系统进行仿真研究，可以排除人的主观因素的影响。（4）在实际系统上做实验时，由于系统

5、误差和偶然误差的存在，我们很难正确地判断实验结果的好坏，这时运用计算机模拟，就可以保证每次操作的条件相同，排除误差。12.1.2 数字仿真的分类计算机模拟分为动态模拟和静态模拟，数值分析中的蒙特卡罗法就是典型的静态模拟，动态模拟可以分为连续系统模拟和离散事件系统模拟。在连续系统模拟和离散事件系统模拟中时间具有重要作用，我们都是要观察系统在时间过程中的变化。连续系统模拟研究系统的状态随时间连续变化的情况，在解决实际问题时，一般要建立微分方程模型，先确定系统的连续状态变化量，然后将它在时间上进行适当的离散量化处理，并由此模拟系统的运行状态。而离散事件系统模拟讨论的是系统状态只在一些离散时间点上，由

6、于随机时间的推进而发生变化，其问题解决模型一般用流程图或网络来表示。连续系统模拟中虽然有时也考虑一些随机作用，但经常还是把它当作确定性问题去考虑，而在离散事件系统模拟中事件的出现和系统状态变量几乎总都是随机的。12.1.3 数字仿真的原理 现实世界充满不确定性，我们所研究的现实对象往往难以摆脱随机因素的影响。概率论是用数学的思想和方法处理和研究随机现象的一个有效的工具。但有时它还难以用来处理复杂系统中的随机性。数字仿真是使用计算机研究和处理复杂系统的随机性一条重要途径。对随机现象进行模拟，实质上要给出随机变量的模拟。随机数的产生方法主要有以下几种：物理方法：一是放射性物质随机蜕变；二是电子管回

7、路的热噪声，其噪声电压的大小表示不同的随机数。查随机数表（Rand Table）（1955年由美国兰德公司编制，有随机数100万个。）随机数表中的数字具有均匀的随机性，没有周期性。按照一定算法在计算机内产生伪随机数：由于第i+1个随机数是由第i个按一定公式推算出来的，故并非真正的随机数。如表12-1所示，利用MATLAB中命令可生成各种分布下产生随机数。表12-1 MATLAB中各种分布下随机数生成命令常见的分布函数MATLAB语句常见的分布函数MATLAB语句均匀分布均匀分布指数分布柏松分布标准正态分布正态分布二项分布在实际应用中用哪种随机数生成法，要针对具体的系统做出与之相适应的选择。例如

8、：（1）当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。（2）考试分数的偏差、射击命中点与目标的偏差、人的身高、体重等，都可近似看成服从正态分布。（3）排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。指数分布在排队论、可靠性分析中有广泛应用。（4）涉及到排队系统、产品检验时可用到泊松分布。12.2 数字仿真的方法 12.2.1 解决计算机模拟的一般步骤：进行数字仿真一般要进过四个步骤：系统分析，模型构造，运行与改进和输出结果。图12-1中表示系统分析，表示模型构造，表示运行与改进，表示输出结果。第一

9、步：系统分析，明确目标。就是要明确问题和提出总体方案。首先要把被仿真系统的内容表达清楚，弄清仿真的目的，然后选择描述这些目标的主要环节和状态变量，明确定义所研究问题的范围、边界和初始条件，并充分估计初始条件对系统主要性能的影响。第二步：模型构造、采集数据。包括建立模型、收集数据、编写程序、程序验证和模型确认等。建立模型就是选择合适的仿真方法，如时间步长法、事件表法等，确定系统的初始状态， 设计整个系统的仿真流程。最后选择合适的通用语言或仿真语言编写、调试程序。第三步：模型的运行与改进。首先确定一些具体的运行方案，如初始条件、参数、步长、重复次数等，然后输入数据，运行程序，直到符合实际系统的要求

10、及精度为止。第四步：模型输出、统计分析。包括提供文件的清单，记录重要的中间结果，输出格式要有利于用户了解整个仿真过程，分析和使用仿真结果。12.2.2蒙特卡罗法蒙特卡罗法（Random simulation）又称随机抽样技巧或统计试验方法，它利用随机数进行统计试验，以求得的统计特征值（如均值、概率等）作为待解决问题的数值解。蒙特卡罗法的基本思想原理是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量。然后通过模拟、统计，即多次随机抽样实验，统计出某事件发生的百分比。只要实验次数很大，该百分比便近似于事件发生的概率。一、蒙特卡罗法的原理根据车比雪夫定理，设,是相独立的随机变

11、量序列，它们服从相同的分布，且有有限的数学期望和方差 ，则,的算术平均值当时按概率1收敛于，即对于任意有：由中心极限定理得到：图12-2 针在平行线间的位置 即当n很大时，近似服从标准正态分布。二、蒙特卡罗法的实例1、浦丰(Buffon)投针问题例12.1如图12-2所示，若干条相互平行线间的距离为，针的长度为，将针投出后其与平行线相交的概率为,求出其概率。模型建立 通过计算可知针与平行线相交的概率为将针投次，若有次与平行线相交，则可用作为的近似值。由此可计算出的近似值为要使结果更加准确就需使尽量大一些。计算机模拟 如图12-2所示，设针投到地面上的位置可以用一组参数（x,）来描述，x为针中心

12、的坐标，为针与平行线的夹角。任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。假设是服从区间0，a上的均匀分布，是服从区间0，上的均匀分布。在此情况下，针与平行线相交的数学条件是按照某种抽样法，产生随机变量的可能值，例如进行次抽样，得到样本值，其中，统计出满足不等式的次数，然后利用MATLAB进行计算机编程模拟。%浦丰(Buffon)投针问题程序a=2; % 2a为平行线间距离l=1; % 2l为针长N=100000; % N为试验次数x=a*rand(1,N); % 针中心纵坐标theta=pi*rand(1,N); % 针倾斜角n=sum(x<=sin(

13、theta); % n为相交次数P=n/N %相交的概率C=2*l/a*N/n % 圆周率的近似值当取，时，时，模拟得到：概率，圆周率, 这与真实值已经很接近了。2、定积分计算我们知道在近似计算中Simpson方法是一种常用的数值积分方法，但对于不规则的被积函数及多维情形这种方法往往无能为力。而蒙特卡洛方法计算这种积分就较为方便了，虽然其精度较差，但在某些场合往往有其出现的必要性。考虑计算积分，其中为计算上述积分，我们考虑取矩形区域，设为在上均匀分布的随机变量，它有概率密度则落入的概率为即.那么我们只要能估计出即可给出的近似值。S（HitMiss图12-3 随机投点法如图12-3所示，为估计，

14、产生个相互独立的在上均匀分布的随机数 。用表示满足的数目，称为“Hit”的数目，根据大数定律，可以用估计，从而得到的近似值例12.2 试用蒙特卡洛模拟中的随机投点法计算积分.计算机模拟 编写MATLAB程序如下：% 随机投点法计算积分nt=0; %the total number of the test nofin=0; %the number of point in the area for x=0:0.1:5; for tmp=0:0.1:50 y=500*rand(1); y0=4*x3; nt=nt+1; if(y<=y0) nofin=nofin+1; plot(x,y,

15、9;r'); hold on end end end sofsq=x*y0; %计算矩形的面积p=nofin/nt; s=p*sofsq运行结果为647.7242，并显示图12-4.图12-4随机投点法计算计算积分的另一种方法是将其表示成某个随机变量的数学期望，我们只要将积分改写为这里假定为任意概率密度函数且满足，这样就可以将积分表示为如下的数学期望一般情况下，我们可以取X为上均匀分布的随机变量，即于是有为了估计积分，只要产生N个上均匀分布的随机数，并用样本均值来估计积分即可。例12.3 试用蒙特卡洛模拟中的平均值法计算积分计算机模拟 编写MATLAB程序如下：% 平均值法计算积分su

16、m=0;N=50000; %模拟次数x=5*rand(1,N);y=20.*x.3;I=mean(y)运行结果为621.4690，与真值的相对误差为0.56%.12.2.3离散型随机变量的模拟离散系统是指系统状态只在有限的时间点或可数的时间点上有随机事件发生的系统，在这些点上发生的离散事件改变了系统的状态，并假设系统状态的变化是在该时间点上瞬时完成的，即离散事件是瞬时出现的。一、离散型随机变量模拟的原理设离散随机变量的可能值为，相应的概率为，设， ，将作为分点，把区间分为一系列小区间。设是上均匀随机变量，则有：=，因此可以用随机变量落在内的情况来模拟离散的随机变量的取值情况。二、离散型随机变量

17、模拟的实例1、系统可靠性分析例12.4某生产电子产品的企业，要对某型号的产品平均无故障运行时间做出估计。该产品由A、B、C三个部件串联而成。因此，当这三个部件中任何一个部件发生故障而失效时，则该电子产品也即告失效。如果根据该产品投入运行后再对其无故障运行时间做出估计，则费用较高。现在企业已经得到每一种部件的有关运行试验记录资料，其中包括用来确定部件失效时间的概率分布。表12-2 部件失效概率分布和随机数取值表A部件B部件C部件失效时间概率随机数失效时间概率随机数失效时间概率随机数40.101020.0500560.202050.2113030.1061570.3215060.3316040.2

18、163580.25517570.2618050.3366590.15769080.15819560.256690100.1910090.05960070.19100模型建立 设为A、B、C三个部件运行时间失效的概率。为A、B、C三个部件运行超过时间的概率。因为该产品由A、B、C三个部件串联而成，故该产品运行超过时间的概率为该产品运行时间失效的概率为，其中。计算机模拟 利用MATLAB进行计算，编程如下：% 产品失效问题% PA,PB,PC为产品A,B,C失效概率分布PA=0 0 0 0.1 0.2 0.3 0.2 0.15 0.05 0;PB=0 0.05 0.1 0.2 0.3 0.25 0

19、.1 0 0 0;PC=0 0 0 0 0 0.2 0.3 0.25 0.15 0.1;%计算产品运行时间超过t的概率F=zeros(1,10);for t=1:10 F(t)=(1-sum(PA(1:t)*(1-sum(PB(1:t)*(1-sum(PC(1:t);end%计算产品运行t时间失效的概率P=1,F(1:9)-F(1:10)t=1:10;Et=sum(t.*P) %计算产品运行失效的平均时间经计算，该产品运行失效的概率分布如下表：表12-3 产品运行失效的概率分布失效时间02345678910概率00.05000.10000.26500.34000.21300.0320000该产

20、品运行失效平均时间为此题也可用数字模拟方法。计算机模拟 随机生成一系列服从0, 100均匀分布的整数向量,查表12-2中随机数分段范围，确定部件A、B、C失效的时间。例如，表12-4列出了10次模拟实验结果。表12-4 产品运行失效模拟实验结果产品失效时间仿真表序号ABC产品失效时间随机数失效时间随机数失效时间随机数失效时间1336244528425067268596313519479944828204869455969177286630588620657245957126580243852174915541599105103865155885经计算，10次模拟得到产品平均失效时间为4.8小时

21、。然后，利用MATLAB进行计算机编程模拟。程序如下：function PFail()% 产品失效问题计算机模拟N=100; %模拟次数ZN=ceil(100*rand(N,3); %生成一系列随机整数% 计算部件失效时间tA=PAFail(ZN(:,1);tB=PBFail(ZN(:,2);tC=PCFail(ZN(:,3);% 计算产品失效时间t=min(tA;tB;tC,1);fprintf('产品失效平均时间为%3.2fn',mean(t);end function t=PAFail(n)t(n>=0 & n<=10)=4;t(n>=11 &a

22、mp; n<=30)=5;t(n>=31 & n<=60)=6;t(n>=61 & n<=80)=7;t(n>=81 & n<=95)=8;t(n>=96 & n<=100)=9;end function t=PBFail(n)t(n>=0 & n<=5)=2;t(n>=6 & n<=15)=3;t(n>=16 & n<=35)=4;t(n>=36 & n<=65)=5;t(n>=66 & n<=90)=6;t

23、(n>=91 & n<=100)=7;end function t=PCFail(n)t(n>=0 & n<=20)=6;t(n>=21 & n<=50)=7;t(n>=51 & n<=75)=8;t(n>=76 & n<=90)=9;t(n>=91 & n<=100)=10;end模拟100次，得到产品运行失效平均时间为4.70. 与真值4.6620的相对误差为0.82%.2、存贮系统的模拟例12.5 有某种货物的存贮系统，市场对这种货物的需求量和订货提前期都是随机的，它们

24、的概率分布如下：现在考虑订货、存贮、缺货损失三项费用：订货费用每次25元，订货量每次20单位，订货点为15单位。（即存货低于15单位时订货，但已订货未到前不再订）存贮费每件每周10元，缺货损失费每件每周500元。对于缺货，货到后不补，设开始时存货为20单位。试模拟计算14周的运行情况：并求订货费用、存贮费用、缺货费用以及周平均费用。模型建立 设订货费用每次元，存贮费每件每周元，缺货损失费每件每周元。记订货量每次单位，订货点为单位，周末的存货量为单位。按照存贮策略的要求，当周末的存货量时，订货量为；当且未定货时，订货量为，当且订货未到时，订货量为.单位设第k周末的存货量为单位，第k周末的订货到货

25、量为单位，第k周末的订货提前期为周，开始时， ，。当时，；当时，若则更新，并且，否则，。设第k周的需求量为单位，则.假设周初交付本周需求量，从而第k周的订货费用为，当时，第k周的存贮费用为，当时，第k周的缺货损失费用为，故第k周的费用为.计算机模拟 设开始时存货为20单位。随机生成一系列服从0, 1均匀分布的整数向量,分别查看表12-5和表12-6中随机数所对应累积概率范围，确定各周需求量和订货提前期。然后根据所建模型计算各周的订货费用、存贮费用、缺货费用以及周平均费用。编程如下：function store()% 存贮系统的模拟N=10; %模拟次数c1=25; %订货费用()c2=10;

26、%存贮费(元/(件周)c3=500; %缺货损失费(元/(件周)S=25; %每次订货量（件）s=15; %订货点（件）% 初始化需求量，到货量，存储量x=20; %初始的存货量（件）u=zeros(N,1); %初始的到货量（件）t=0; %初始的订货提前期（周）C=zeros(N,1); % 模拟过程disp('周次 需求量 提前期 到货量 存储量 订货费用 存贮费用 缺货费用 周费用'); for k=1:N d=need(rand(); % 计算各周需求量 x=x-d; Corder=c1*u(k); % 订货费用 Cstore=c2*max(x,0); % 存贮费用

27、Clost=c3*max(-x,0); % 缺货费用 C(k)=Corder+Cstore+Clost; % 周费用 if t>0 t=t-1; else if x<s t=advance(rand(); % 计算订货提前期 u(k+t)=S; else t=0; end end x=x+u(k); fprintf('%3d %d %2d %2d %2d %5.2f %5.2f %5.2f %5.2fn', k, d, t, u(k), x, Corder, Cstore, Clost, C(k);end fprintf('周平均费用为%5.2fn'

28、,sum(C)/N);end % 子函数计算需求量function d=need(r)d=zeros(size(r);d(r>0.02 & r<=0.1)=1;d(r>0.1 & r<=0.32)=2;d(r>0.32 & r<=0.66)=3;d(r>0.66 & r<=0.84)=4;d(r>0.84 & r<=0.93)=5;d(r>0.93)=6;end % 子函数计算订货提前期function t=advance(r)t=zeros(size(r);t(r>0.23 &am

29、p; r<=0.48)=1;t(r>0.48 & r<=0.78)=2;t(r>0.78 & r<=0.9)=3;t(r>0.9 & r<=0.96)=4;t(r>0.96)=5;end运行后的结果见表12-7. 周平均费用为237.50元。表12-7 存贮系统模拟过程周次需求量提前期到货量存储量订货费用存贮费用缺货费用周费用1300170.00170.000.00170.002520120.00120.000.00120.00331090.0090.000.0090.004502529625.0040.000.00665

30、.005400250.00250.000.00250.006300220.00220.000.00220.007200200.00200.000.00200.008200180.00180.000.00180.0095025380.00130.000.00130.0010300350.00350.000.00350.003、排队系统在离散型随机变量模拟所解决的问题中，排队系统是一类很常见又非常重要的问题，且这类问题广泛存在于银行、医院、超市等现实生活中的领域。由一个或多个服务台和一些等待服务的顾客组成的离散时间系统称为排队系统。这里列举一个简单的排队系统问题，并从中了解离散型随机变量模拟的思想

31、方法。例12.6一个理发店有两位服务员和,顾客随机地到达该理发店，每分钟有一个顾客到达和没有顾客到达的概率均是，其中的顾客理发仅用分钟,另外的顾客用分钟。试对前分钟的情况进行仿真，模拟出完成服务的个数及顾客平均等待时 间。模型建立 本题属于多服务台随即服务模型。假设客源是无穷的，且队长不受限制，顾客接受先来后到的服务规则。设顾客总的等待时间为，第个顾客的到达时刻为，第个顾客开始接受服务的时刻为，第个顾客服务结束的时刻为，第个顾客与第个顾客之间的时间间隔为，第个顾客的服务时间为。这样便容易得到如下关系：计算机模拟 上述排队系统的仿真算法步骤如下：（1）将时间区间0,60等分为60份，即每个时间间

32、隔为1分钟。（2）应用MATLAB随机数生产函数，按照等概率原则，模拟出前分钟中顾客到达情况： 、；（3）再按照的顾客理发仅用分钟,另外的顾客用分钟的顾客接受服务时间原则，模拟出顾客 、所对应的接受服务的时间、；（4）确定每个顾客接受服务的时间。当系统刚开始运行时，第一、第二个顾客可以直接接受服务，而对于以后的顾客，其接受服务的时刻的判断取决于三者中最大者。即：；（5）计算顾客的累计等待时间。顾客总的等待时间为，得到：；（6）判断第个顾客开始接受服务的时刻为是否超出分钟。若是，则总得接受服务的人数为：,顾客平均等待时间为：；若否，重复步骤（3）。编写进行数字模拟的MATLAB程序如下：% 排队

33、系统模拟clear allN=60; %模拟阶段数num=0; %完成服务的人数 %模拟每阶段顾客到达情况，当r1=1时，有顾客到达，否则，无顾客到达r1=binornd(1,0.5,N,1); c=find(r1); %确定顾客到达的时间m=length(c); %计算0, N内到达顾客数%模拟每位顾客服务情况，当时，服务时间为5分钟，否则，服务时间为8分钟r2=binornd(1,0.6,m,1); y=8*ones(m,1);y(r2=1)=5;%初始化前两位顾客服务结束时刻b=zeros(m,1); %顾客服务时间e=zeros(m,1); %顾客服务结束时刻w=zeros(m,1);

34、 %顾客等待服务时间b(1)=c(1);e(1)=b(1)+y(1); num=num+1;w(1)=0;b(2)=max(c(2),e(1);e(2)=b(2)+y(2); num=num+1;w(2)=b(2)-c(2);for k=3:m b(k)=max(c(k),e(k-1),e(k-2); e(k)=b(k)+y(k); if e(k)>N break; end num=num+1; w(k)=b(k)-c(k);enddisp('序号 到达时间 服务开始时间 服务结束时间 等待服务时间');for k=1:num fprintf('%2d %6d %

35、6d %6d %6dn',k,c(k),b(k),e(k),w(k)endfprintf('%d分钟内完成服务的人数为%dn',N,num);fprintf('顾客平均等待时间为%2.1f分钟n',mean(w); 运行后的结果见表12-8，60分钟内完成服务的人数为12，顾客平均等待时间为4.4分钟。表12-8 排队系统模拟过程序号到达时间服务开始时间服务结束时间等待服务时间133702471233101217241217225513222696162631107233136882636411092841461310304651161131515520

36、123255592312.2.4连续型随机变量的模拟连续系统模拟研究系统的状态随时间连续变化的情况。连续系统虽然研究的是确定性问题，但在实际处理中，是将连续状态变量在时间上进行离散化处理，并由此模拟系统的运行状态，即其状态由于随机时间的推进而发生变化。因此，连续系统的计算机模拟只能是近似的，只要这种近似达到一定的精度，一般就可以满足要求。一、连续型随机变量模拟的原理处理连续型随机变量模拟问题有多种方法，其中反函数法是最常用的。其方法核心是通过求概率分布的反函数产生随机数。众所周知随机变量的概率分布函数为定义在 区间的单调递增函数，设为区间的均匀随机变量, 令,只要求出反函数,即为具有概率分布函

37、数的随机机数。由概率论的理论可以证明和相同的概率分布。若的概率密度为，由,是区间上均匀分布的随机变量，如果给定区间上均匀分布的随机数, 则具有给定分布的随机数可由方程中解出：二、连续型随机变量模拟的实例应用1、盐水浓度问题例12.7 在水池中有2000m3水，其中含盐2kg。现以每分钟6m3的速度向池中注入含盐率为0.5kg/m3的盐水，同时又以每分钟4m3的速度从水池流出搅拌均匀的盐水。求出池中含盐量随时间的变化规律，并算出含盐率达到0.2kg/m3时所用的时间。问题分析 初始时刻水池中含盐率为0.001kg/m3，注入盐水的浓度为0.5kg/m3，故水池中盐水的浓度逐渐升高。模型构建与求解

38、 设t分钟时刻水池中含有kg的盐。考察t, t+t时间内水池中含盐量的变化，有：两边除以t，且使t0，则可得微分方程模型：在MATLAB命令窗口中输入如下指令：x=dsolve('Dx=3-2*x/(1000+t)','x(0)=2')运行后输出结果：x = 1000+t-998000000/(1000+t)2结果表明，微分方程的解为最后，由 得容易算出含盐率达到0.2kg/m3时所用的时间为t184.84（分钟）计算机仿真 设初始时刻时间为，时间步长为1分钟，设第k分钟池水中含盐量为，池水体积为，含盐率为，假设在时间区间k,k+1)池水中的含盐率不变，则每10

39、分钟输出一次体积、含盐量、含盐率，直到含盐率达到0.2 kg/m3，结束程序，流程图如图12-5所示。图12-5 盐水浓度随时间变化模拟流程图根据流程图12-5，编写MATLAB程序如下：clc, clear all;time=0; %时间初始化volume=2000; %状态初始化salt=2; rate=salt/volume;disp('时间 水的体积 含盐量 含盐率'); %模拟开始while (rate<0.2) time=time+1; volume=volume+2; salt=salt+3-4*rate; rate=salt/volume; if mod(

40、time,10)=0 fprintf('%4d %6d %4.2f %3.2fn',time,volume,salt,rate); endendfprintf('%4d %6d %4.2f %3.2fn',time,volume,salt,rate);运行后输出结果见表12-9.表12-9 池中盐水变化情况时间水的体积含盐量含盐率时间水的体积含盐量含盐率10202031.690.01571102220300.240.135220204060.810.02981202240324.660.144930206089.370.04341302260348.690.15

41、43402080117.400.05641402280372.350.1633502100144.910.06901502300395.660.1720602120171.930.08111602320418.630.1804702140198.480.09271702340441.270.1886802160224.570.10401802360463.580.1964902180250.210.11481852370474.620.20031002200275.430.1252表12-9中结果表明，第185分钟时，池水中含盐量刚好超过0.2 kg/m3，这与精确解184.84分钟十分的接近

42、。2、追击性问题例12.8 正方形的个顶点各有人。在某一时刻，人同时出发以匀速，按顺时针方向追逐下一人，如果他们始终保持对准目标，则最终按螺旋状曲线于中心点。试求出这种情况下每个人的行进轨迹。模型建立 建立平面直角坐标系，设,，取时间间隔为，计算每一点在各个时刻的坐标。设某点在时刻的坐标为，则在时刻的坐为其中，。计算机模拟 设置初始化条件为，取足够小的,当时算法结束。对每一个点连接它在各时刻的位置,即得所求运动轨迹。然后根据上述建立的模型，用MATLAB编程。%追击问题仿真程序v=1;dt=0.05;x=0 0 10 10;y=0 10 10 0; for i=1:4 plot(x(i),y(

43、i),'.'),hold onend d=20;subplot(2,2,4);while(d>0.1) x(5)=x(1);y(5)=y(1); for i=1:4 d=sqrt(x(i+1)-x(i)2+(y(i+1)-y(i)2); x(i)=x(i)+v*dt*(x(i+1)-x(i)/d; y(i)=y(i)+v*dt*(y(i+1)-y(i)/d; plot(x(i),y(i),'.'),hold on endend图12-6 追击问题仿真分解图将整个仿真过程分解(如图12-6)，我们可以直观地看到整个追击过程。为了更直观的看到整个追击过程，利用

44、MATLAB编制的仿真动画程序。%追击问题的仿真动画程序v=1;dt=0.05;x=0 0 10 10;y=0 10 10 0;for i=1:4 plot(x(i),y(i),'.'),hold onend d=20;axis equal,M=moviein(16); while(d>0.1) x(5)=x(1);y(5)=y(1); for i=1:4 d=sqrt(x(i+1)-x(i)2+(y(i+1)-y(i)2); x(i)=x(i)+v*dt*(x(i+1)-x(i)/d; y(i)=y(i)+v*dt*(y(i+1)-y(i)/d; plot(x(i),y

45、(i),'.'), M(:,i)=getframe; endhold onendmovie(M,30)通过仿真动画,追击的过程可以一目了然地表现出来。12.3 系统仿真软件SIMULINK12.3.1 SIMULINK简介SIMULINK是MATLAB中一个用来对动态系统进行建模、仿真和分析的软件包。SIMULINK 中的“SIMU”一词表示可用于计算机仿真，而“LINK”一词表示它能进行系统连接，即把一系列模块连接起来，构成复杂的系统模型。它支持线性和非线性系统，连续和离散模型等。SIMULINK可以使用显示模块可以在运行仿真时观察到仿真结果，还可以改变参数并立即就可以观察到

46、它的变化。SIMULINK提供有可自定义的定制模块库、模型编辑器和模型求解器。模块库的界面为 SIMULINK模块库浏览器，使用命令simulink或工具条中的快捷钮可以进入模块库浏览器。如图12-7所示，在该浏览器中，模块库以文件夹形式组织，由SIMULINK通用库和专业库组成，而模块单元以文件形式保存。图12-7 Simulink模块库浏览器如图12-7所示，SIMILINK模块通用库按功能进行分为：Continuous（连续模块），discontinuous（非连续模块），Discrete（离散模块），Logic and Bit Operations（逻辑和位操作），Lookup Tab

47、les（查找表模块），Math Operations（数学模块），Model Verification（模型检验模块），Model-Wide Utilities（模型扩充模块），Ports&Subsystems（端口和子系统模块），Signal Attributes（信号属性模块），Signal Routing（信号线路模块），Sinks（接收器模块），Sources（输入源模块）和User Defined Function（用户自定义函数模块）等。在SIMILINK中，各功能模块的参数描述都可以由用户通过该模块的模块属性对话框进行操作给出或修改。在“File”下点击“NewModel

48、”或单击模块库浏览器的新建图标就可以进入模型空间组建模型，如图12-8所示。在模型编辑器中，可以直接从模块库中将模块用鼠标点击和拖拉过来就可以拼接成所需要的模型，它可以以M文件的形式保存下来。图12-8 模型编辑器SIMILINK的模型窗口的主菜单有文件、编辑、查看、仿真、格式设定、工具与帮助七项菜单选项。模型求解器位于模型浏览器的菜单“Simulation”中,Simulation(仿真)菜单包括：l start (Ctrl +T)：启动或暂停仿真。l stop：停止仿真。l Configuration Parameters (Ctrl +E)：设置仿真参数l Normal：常规标准仿真。l

49、 Accelerator：加速仿真。12.3.2 SIMULINK建模一个典型的SIMULINK模型包括如下三种要素： 信号源模块 被模拟的系统模块 输出显示模块 系统模块作为中心模块是SIMULINK仿真建模所要解决的主要部分；信号源为系统的输入，它包括常数信号源函数信号发生器（如正弦和阶跃函数波等）和用户自己在Matlab中创建的自定义信号或Matlab工作间中三种。输出模块主要在Sinks库中。SIMULINK模型并不一定要包含全部的三种要素，在实际应用中通常可以缺少其中的一个或两个。例如，若要模拟一个系统偏离平衡位置后的恢复行为，就可以建立一个没有输入而只有系统模块加一个显示模块的模型

50、。在某种情况下，也可以建立一个只有源模块和显示模块的系统。若需要一个由几个函数复合的特殊信号，则可以使用源模块生成信号并将其送入Matlab工作间或文件中。利用SIMULINK建立物理系统和数学系统的仿真模型，关键是对SIMULINK提供的功能模块进行操作，即用适当的方式将各种模块连接在一起。利用SIMULINK建模与仿真的步骤为： 启动SIMULINK，打开SIMULINK模块库； 打开空白模型窗口，用鼠标拖拉模块构成模型； 设置仿真参数，进行仿真； 输出仿真结果，以观察仿真的结果。下面以数学方程的建模为例，说明SIMULINK的建模过程。例12.9 构建摄氏向华氏温度的变换器分析：借助SI

51、MULINK模块化框图实现TF=9/5(TC)+32。 功能框图如图12-9所示：图12-9 方程TF=9/5(TC)+32的功能框图解：步骤1：模块定位与选取如表所示。功能需求模块库模块单元摄氏温度（信号）sourceRamp(斜坡发生器)信号增益mathgain常量元素mathconstant求和运算mathsum信号显示sinkScope(示波器)步骤2：按功能框图配置模块元件。打开模块库浏览器，进入模型空间组建模型。直接从模块库中将模块拖拉过来拼接成所需要的模型，如图12-10所示。图12-10 模型所需模块元件图步骤3: 添加信号线。 如图12-11所示，共增加4条线路。线路1：信号发生器信号增益，线路2：信号增益信号求和单元，线路3：常量单元信号求和单元，线路4：信号