ch12 數(shù)字仿真實驗

1、實驗十二 數(shù)字仿真實驗38612.1數(shù)字仿真的基本概念38612.1.1 所謂數(shù)字仿真38612.1.2 數(shù)字仿真的分類38712.1.3 數(shù)字仿真的原理38712.2 數(shù)字仿真的方法38812.2.1 解決計算機模擬的一般步驟：38812.2.2蒙特卡羅法38812.2.3離散型隨機變量的模擬39212.2.4連續(xù)型隨機變量的模擬40012.3 系統(tǒng)仿真軟件SIMULINK40412.3.1 SIMULINK簡介40412.3.2 SIMULINK建模40512.3.3 SIMULINK仿真407實驗任務(wù)408411實驗十二 數(shù)字仿真實驗學(xué)習(xí)目的了解計算機模擬仿真在解決實際問題時的重要性，理

2、解使用計算機仿真的原理及方法。掌握計算機模擬仿真的三大類方法：蒙特卡羅法、連續(xù)系統(tǒng)模擬和離散事件系統(tǒng)模擬。會利用MATLAB進行模擬仿真實驗。實驗指導(dǎo)在實際問題中，我們通常會面對一些帶隨機因素的復(fù)雜系統(tǒng)，用分析方法建模常常需要作許多簡化假設(shè)，這樣進行處理過后的模型與我們面臨的實際問題可能相差很遠，以致求解得到答案根本無法應(yīng)用，這時，計算機模擬幾乎成為唯一的選擇。12.1數(shù)字仿真的基本概念我們將一些具有特定的功能、相互之間有一定的內(nèi)在聯(lián)系的對象所組成的總體稱為一個系統(tǒng)。模擬就是利用物理的、數(shù)學(xué)的模型類比、模仿現(xiàn)實系統(tǒng)及其演變過程，以尋求變化規(guī)律的一種方法。模擬的基本思想是建立一個實驗的模型，這個

3、模型包含所研究系統(tǒng)的主要特點，這樣做的目的就是通過對這個實驗?zāi)Ｐ偷倪\行，獲得所要研究系統(tǒng)的必要信息。12.1.1 所謂數(shù)字仿真數(shù)字仿真是利用計算機對一個實際系統(tǒng)的結(jié)構(gòu)和行為進行動態(tài)演示，以評價或預(yù)測該系統(tǒng)的行為效果。它是解決較復(fù)雜的實際問題的一條有效途徑。數(shù)字仿真的一般過程是先對研究對象進行數(shù)學(xué)描述，然后構(gòu)建仿真算法，最后在計算機中編程并運行實現(xiàn)。對比于物理模擬通常花費較大、周期較長，且在物理模型上改變系統(tǒng)結(jié)構(gòu)和系數(shù)都較困難的諸多缺陷，數(shù)字模擬花費較少、周期較短，且易修改、可重用，有更強的系統(tǒng)適應(yīng)能力。但是計算機模擬也有缺陷，比如系統(tǒng)數(shù)學(xué)模型不易建立、程序調(diào)試復(fù)雜等。數(shù)字仿真可以用于研制產(chǎn)品或

4、設(shè)計系統(tǒng)的全過程中，包括方案論證、技術(shù)指標(biāo)確定、設(shè)計分析、生產(chǎn)制造、試驗測試、維護訓(xùn)練、故障處理等各個階段。具體來說，下列情形適用于數(shù)字仿真：（1）在一個實際系統(tǒng)還沒有建立起來之前，要對系統(tǒng)的行為或結(jié)果進行分析研究時，數(shù)字仿真是一種行之有效的方法。（2）在有些真實系統(tǒng)上做實驗會影響系統(tǒng)的正常運行，這時進行計算機模擬就是為了避免給實際系統(tǒng)帶來不必要的損失。如在生產(chǎn)中任意改變工藝參數(shù)可能會導(dǎo)致廢品，在經(jīng)濟活動中隨意將一個決策付諸行動可能會引起經(jīng)濟混亂。（3）當(dāng)人是系統(tǒng)的一部分時，他的行為往往會影響實驗的效果，這時運用系統(tǒng)進行仿真研究，可以排除人的主觀因素的影響。（4）在實際系統(tǒng)上做實驗時，由于系統(tǒng)

5、誤差和偶然誤差的存在，我們很難正確地判斷實驗結(jié)果的好壞，這時運用計算機模擬，就可以保證每次操作的條件相同，排除誤差。12.1.2 數(shù)字仿真的分類計算機模擬分為動態(tài)模擬和靜態(tài)模擬，數(shù)值分析中的蒙特卡羅法就是典型的靜態(tài)模擬，動態(tài)模擬可以分為連續(xù)系統(tǒng)模擬和離散事件系統(tǒng)模擬。在連續(xù)系統(tǒng)模擬和離散事件系統(tǒng)模擬中時間具有重要作用，我們都是要觀察系統(tǒng)在時間過程中的變化。連續(xù)系統(tǒng)模擬研究系統(tǒng)的狀態(tài)隨時間連續(xù)變化的情況，在解決實際問題時，一般要建立微分方程模型，先確定系統(tǒng)的連續(xù)狀態(tài)變化量，然后將它在時間上進行適當(dāng)?shù)碾x散量化處理，并由此模擬系統(tǒng)的運行狀態(tài)。而離散事件系統(tǒng)模擬討論的是系統(tǒng)狀態(tài)只在一些離散時間點上，由

6、于隨機時間的推進而發(fā)生變化，其問題解決模型一般用流程圖或網(wǎng)絡(luò)來表示。連續(xù)系統(tǒng)模擬中雖然有時也考慮一些隨機作用，但經(jīng)常還是把它當(dāng)作確定性問題去考慮，而在離散事件系統(tǒng)模擬中事件的出現(xiàn)和系統(tǒng)狀態(tài)變量幾乎總都是隨機的。12.1.3 數(shù)字仿真的原理 現(xiàn)實世界充滿不確定性，我們所研究的現(xiàn)實對象往往難以擺脫隨機因素的影響。概率論是用數(shù)學(xué)的思想和方法處理和研究隨機現(xiàn)象的一個有效的工具。但有時它還難以用來處理復(fù)雜系統(tǒng)中的隨機性。數(shù)字仿真是使用計算機研究和處理復(fù)雜系統(tǒng)的隨機性一條重要途徑。對隨機現(xiàn)象進行模擬，實質(zhì)上要給出隨機變量的模擬。隨機數(shù)的產(chǎn)生方法主要有以下幾種：物理方法：一是放射性物質(zhì)隨機蛻變；二是電子管回

7、路的熱噪聲，其噪聲電壓的大小表示不同的隨機數(shù)。查隨機數(shù)表（Rand Table）（1955年由美國蘭德公司編制，有隨機數(shù)100萬個。）隨機數(shù)表中的數(shù)字具有均勻的隨機性，沒有周期性。按照一定算法在計算機內(nèi)產(chǎn)生偽隨機數(shù)：由于第i+1個隨機數(shù)是由第i個按一定公式推算出來的，故并非真正的隨機數(shù)。如表12-1所示，利用MATLAB中命令可生成各種分布下產(chǎn)生隨機數(shù)。表12-1 MATLAB中各種分布下隨機數(shù)生成命令常見的分布函數(shù)MATLAB語句常見的分布函數(shù)MATLAB語句均勻分布均勻分布指數(shù)分布柏松分布標(biāo)準(zhǔn)正態(tài)分布正態(tài)分布二項分布在實際應(yīng)用中用哪種隨機數(shù)生成法，要針對具體的系統(tǒng)做出與之相適應(yīng)的選擇。例如

8、：（1）當(dāng)研究對象視為大量相互獨立的隨機變量之和，且其中每一種變量對總和的影響都很小時，可以認(rèn)為該對象服從正態(tài)分布。（2）考試分?jǐn)?shù)的偏差、射擊命中點與目標(biāo)的偏差、人的身高、體重等，都可近似看成服從正態(tài)分布。（3）排隊服務(wù)系統(tǒng)中顧客到達率為常數(shù)時的到達間隔、故障率為常數(shù)時零件的壽命都服從指數(shù)分布。指數(shù)分布在排隊論、可靠性分析中有廣泛應(yīng)用。（4）涉及到排隊系統(tǒng)、產(chǎn)品檢驗時可用到泊松分布。12.2 數(shù)字仿真的方法 12.2.1 解決計算機模擬的一般步驟：進行數(shù)字仿真一般要進過四個步驟：系統(tǒng)分析，模型構(gòu)造，運行與改進和輸出結(jié)果。圖12-1中表示系統(tǒng)分析，表示模型構(gòu)造，表示運行與改進，表示輸出結(jié)果。第一

9、步：系統(tǒng)分析，明確目標(biāo)。就是要明確問題和提出總體方案。首先要把被仿真系統(tǒng)的內(nèi)容表達清楚，弄清仿真的目的，然后選擇描述這些目標(biāo)的主要環(huán)節(jié)和狀態(tài)變量，明確定義所研究問題的范圍、邊界和初始條件，并充分估計初始條件對系統(tǒng)主要性能的影響。第二步：模型構(gòu)造、采集數(shù)據(jù)。包括建立模型、收集數(shù)據(jù)、編寫程序、程序驗證和模型確認(rèn)等。建立模型就是選擇合適的仿真方法，如時間步長法、事件表法等，確定系統(tǒng)的初始狀態(tài)， 設(shè)計整個系統(tǒng)的仿真流程。最后選擇合適的通用語言或仿真語言編寫、調(diào)試程序。第三步：模型的運行與改進。首先確定一些具體的運行方案，如初始條件、參數(shù)、步長、重復(fù)次數(shù)等，然后輸入數(shù)據(jù)，運行程序，直到符合實際系統(tǒng)的要求

10、及精度為止。第四步：模型輸出、統(tǒng)計分析。包括提供文件的清單，記錄重要的中間結(jié)果，輸出格式要有利于用戶了解整個仿真過程，分析和使用仿真結(jié)果。12.2.2蒙特卡羅法蒙特卡羅法（Random simulation）又稱隨機抽樣技巧或統(tǒng)計試驗方法，它利用隨機數(shù)進行統(tǒng)計試驗，以求得的統(tǒng)計特征值（如均值、概率等）作為待解決問題的數(shù)值解。蒙特卡羅法的基本思想原理是首先建立一個概率模型，使所求問題的解正好是該模型的參數(shù)或其他有關(guān)的特征量。然后通過模擬、統(tǒng)計，即多次隨機抽樣實驗，統(tǒng)計出某事件發(fā)生的百分比。只要實驗次數(shù)很大，該百分比便近似于事件發(fā)生的概率。一、蒙特卡羅法的原理根據(jù)車比雪夫定理，設(shè),是相獨立的隨機變

11、量序列，它們服從相同的分布，且有有限的數(shù)學(xué)期望和方差 ，則,的算術(shù)平均值當(dāng)時按概率1收斂于，即對于任意有：由中心極限定理得到：圖12-2 針在平行線間的位置 即當(dāng)n很大時，近似服從標(biāo)準(zhǔn)正態(tài)分布。二、蒙特卡羅法的實例1、浦豐(Buffon)投針問題例12.1如圖12-2所示，若干條相互平行線間的距離為，針的長度為，將針投出后其與平行線相交的概率為,求出其概率。模型建立 通過計算可知針與平行線相交的概率為將針投次，若有次與平行線相交，則可用作為的近似值。由此可計算出的近似值為要使結(jié)果更加準(zhǔn)確就需使盡量大一些。計算機模擬 如圖12-2所示，設(shè)針投到地面上的位置可以用一組參數(shù)（x,）來描述，x為針中心

12、的坐標(biāo)，為針與平行線的夾角。任意投針，就是意味著x與都是任意取的，但x的范圍限于0，a，夾角的范圍限于0，。假設(shè)是服從區(qū)間0，a上的均勻分布，是服從區(qū)間0，上的均勻分布。在此情況下，針與平行線相交的數(shù)學(xué)條件是按照某種抽樣法，產(chǎn)生隨機變量的可能值，例如進行次抽樣，得到樣本值，其中，統(tǒng)計出滿足不等式的次數(shù)，然后利用MATLAB進行計算機編程模擬。%浦豐(Buffon)投針問題程序a=2; % 2a為平行線間距離l=1; % 2l為針長N=100000; % N為試驗次數(shù)x=a*rand(1,N); % 針中心縱坐標(biāo)theta=pi*rand(1,N); % 針傾斜角n=sum(x<=sin(

13、theta); % n為相交次數(shù)P=n/N %相交的概率C=2*l/a*N/n % 圓周率的近似值當(dāng)取，時，時，模擬得到：概率，圓周率, 這與真實值已經(jīng)很接近了。2、定積分計算我們知道在近似計算中Simpson方法是一種常用的數(shù)值積分方法，但對于不規(guī)則的被積函數(shù)及多維情形這種方法往往無能為力。而蒙特卡洛方法計算這種積分就較為方便了，雖然其精度較差，但在某些場合往往有其出現(xiàn)的必要性?？紤]計算積分，其中為計算上述積分，我們考慮取矩形區(qū)域，設(shè)為在上均勻分布的隨機變量，它有概率密度則落入的概率為即.那么我們只要能估計出即可給出的近似值。S（HitMiss圖12-3 隨機投點法如圖12-3所示，為估計，

14、產(chǎn)生個相互獨立的在上均勻分布的隨機數(shù) 。用表示滿足的數(shù)目，稱為“Hit”的數(shù)目，根據(jù)大數(shù)定律，可以用估計，從而得到的近似值例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運行結(jié)果為647.7242，并顯示圖12-4.圖12-4隨機投點法計算計算積分的另一種方法是將其表示成某個隨機變量的數(shù)學(xué)期望，我們只要將積分改寫為這里假定為任意概率密度函數(shù)且滿足，這樣就可以將積分表示為如下的數(shù)學(xué)期望一般情況下，我們可以取X為上均勻分布的隨機變量，即于是有為了估計積分，只要產(chǎn)生N個上均勻分布的隨機數(shù)，并用樣本均值來估計積分即可。例12.3 試用蒙特卡洛模擬中的平均值法計算積分計算機模擬 編寫MATLAB程序如下：% 平均值法計算積分su

16、m=0;N=50000; %模擬次數(shù)x=5*rand(1,N);y=20.*x.3;I=mean(y)運行結(jié)果為621.4690，與真值的相對誤差為0.56%.12.2.3離散型隨機變量的模擬離散系統(tǒng)是指系統(tǒng)狀態(tài)只在有限的時間點或可數(shù)的時間點上有隨機事件發(fā)生的系統(tǒng)，在這些點上發(fā)生的離散事件改變了系統(tǒng)的狀態(tài)，并假設(shè)系統(tǒng)狀態(tài)的變化是在該時間點上瞬時完成的，即離散事件是瞬時出現(xiàn)的。一、離散型隨機變量模擬的原理設(shè)離散隨機變量的可能值為，相應(yīng)的概率為，設(shè)， ，將作為分點，把區(qū)間分為一系列小區(qū)間。設(shè)是上均勻隨機變量，則有：=，因此可以用隨機變量落在內(nèi)的情況來模擬離散的隨機變量的取值情況。二、離散型隨機變量

17、模擬的實例1、系統(tǒng)可靠性分析例12.4某生產(chǎn)電子產(chǎn)品的企業(yè)，要對某型號的產(chǎn)品平均無故障運行時間做出估計。該產(chǎn)品由A、B、C三個部件串聯(lián)而成。因此，當(dāng)這三個部件中任何一個部件發(fā)生故障而失效時，則該電子產(chǎn)品也即告失效。如果根據(jù)該產(chǎn)品投入運行后再對其無故障運行時間做出估計，則費用較高?，F(xiàn)在企業(yè)已經(jīng)得到每一種部件的有關(guān)運行試驗記錄資料，其中包括用來確定部件失效時間的概率分布。表12-2 部件失效概率分布和隨機數(shù)取值表A部件B部件C部件失效時間概率隨機數(shù)失效時間概率隨機數(shù)失效時間概率隨機數(shù)40.101020.0500560.202050.2113030.1061570.3215060.3316040.2

18、163580.25517570.2618050.3366590.15769080.15819560.256690100.1910090.05960070.19100模型建立 設(shè)為A、B、C三個部件運行時間失效的概率。為A、B、C三個部件運行超過時間的概率。因為該產(chǎn)品由A、B、C三個部件串聯(lián)而成，故該產(chǎn)品運行超過時間的概率為該產(chǎn)品運行時間失效的概率為，其中。計算機模擬 利用MATLAB進行計算，編程如下：% 產(chǎn)品失效問題% PA,PB,PC為產(chǎn)品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;%計算產(chǎn)品運行時間超過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%計算產(chǎn)品運行t時間失效的概率P=1,F(1:9)-F(1:10)t=1:10;Et=sum(t.*P) %計算產(chǎn)品運行失效的平均時間經(jīng)計算，該產(chǎn)品運行失效的概率分布如下表：表12-3 產(chǎn)品運行失效的概率分布失效時間02345678910概率00.05000.10000.26500.34000.21300.0320000該產(chǎn)

20、品運行失效平均時間為此題也可用數(shù)字模擬方法。計算機模擬 隨機生成一系列服從0, 100均勻分布的整數(shù)向量,查表12-2中隨機數(shù)分段范圍，確定部件A、B、C失效的時間。例如，表12-4列出了10次模擬實驗結(jié)果。表12-4 產(chǎn)品運行失效模擬實驗結(jié)果產(chǎn)品失效時間仿真表序號ABC產(chǎn)品失效時間隨機數(shù)失效時間隨機數(shù)失效時間隨機數(shù)失效時間1336244528425067268596313519479944828204869455969177286630588620657245957126580243852174915541599105103865155885經(jīng)計算，10次模擬得到產(chǎn)品平均失效時間為4.8小時

21、。然后，利用MATLAB進行計算機編程模擬。程序如下：function PFail()% 產(chǎn)品失效問題計算機模擬N=100; %模擬次數(shù)ZN=ceil(100*rand(N,3); %生成一系列隨機整數(shù)% 計算部件失效時間tA=PAFail(ZN(:,1);tB=PBFail(ZN(:,2);tC=PCFail(ZN(:,3);% 計算產(chǎn)品失效時間t=min(tA;tB;tC,1);fprintf('產(chǎn)品失效平均時間為%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次，得到產(chǎn)品運行失效平均時間為4.70. 與真值4.6620的相對誤差為0.82%.2、存貯系統(tǒng)的模擬例12.5 有某種貨物的存貯系統(tǒng)，市場對這種貨物的需求量和訂貨提前期都是隨機的，它們

24、的概率分布如下：現(xiàn)在考慮訂貨、存貯、缺貨損失三項費用：訂貨費用每次25元，訂貨量每次20單位，訂貨點為15單位。（即存貨低于15單位時訂貨，但已訂貨未到前不再訂）存貯費每件每周10元，缺貨損失費每件每周500元。對于缺貨，貨到后不補，設(shè)開始時存貨為20單位。試模擬計算14周的運行情況：并求訂貨費用、存貯費用、缺貨費用以及周平均費用。模型建立 設(shè)訂貨費用每次元，存貯費每件每周元，缺貨損失費每件每周元。記訂貨量每次單位，訂貨點為單位，周末的存貨量為單位。按照存貯策略的要求，當(dāng)周末的存貨量時，訂貨量為；當(dāng)且未定貨時，訂貨量為，當(dāng)且訂貨未到時，訂貨量為.單位設(shè)第k周末的存貨量為單位，第k周末的訂貨到貨

25、量為單位，第k周末的訂貨提前期為周，開始時， ，。當(dāng)時，；當(dāng)時，若則更新，并且，否則，。設(shè)第k周的需求量為單位，則.假設(shè)周初交付本周需求量，從而第k周的訂貨費用為，當(dāng)時，第k周的存貯費用為，當(dāng)時，第k周的缺貨損失費用為，故第k周的費用為.計算機模擬 設(shè)開始時存貨為20單位。隨機生成一系列服從0, 1均勻分布的整數(shù)向量,分別查看表12-5和表12-6中隨機數(shù)所對應(yīng)累積概率范圍，確定各周需求量和訂貨提前期。然后根據(jù)所建模型計算各周的訂貨費用、存貯費用、缺貨費用以及周平均費用。編程如下：function store()% 存貯系統(tǒng)的模擬N=10; %模擬次數(shù)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 % 子函數(shù)計算需求量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 % 子函數(shù)計算訂貨提前期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運行后的結(jié)果見表12-7. 周平均費用為237.50元。表12-7 存貯系統(tǒng)模擬過程周次需求量提前期到貨量存儲量訂貨費用存貯費用缺貨費用周費用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、排隊系統(tǒng)在離散型隨機變量模擬所解決的問題中，排隊系統(tǒng)是一類很常見又非常重要的問題，且這類問題廣泛存在于銀行、醫(yī)院、超市等現(xiàn)實生活中的領(lǐng)域。由一個或多個服務(wù)臺和一些等待服務(wù)的顧客組成的離散時間系統(tǒng)稱為排隊系統(tǒng)。這里列舉一個簡單的排隊系統(tǒng)問題，并從中了解離散型隨機變量模擬的思想

31、方法。例12.6一個理發(fā)店有兩位服務(wù)員和,顧客隨機地到達該理發(fā)店，每分鐘有一個顧客到達和沒有顧客到達的概率均是，其中的顧客理發(fā)僅用分鐘,另外的顧客用分鐘。試對前分鐘的情況進行仿真，模擬出完成服務(wù)的個數(shù)及顧客平均等待時 間。模型建立 本題屬于多服務(wù)臺隨即服務(wù)模型。假設(shè)客源是無窮的，且隊長不受限制，顧客接受先來后到的服務(wù)規(guī)則。設(shè)顧客總的等待時間為，第個顧客的到達時刻為，第個顧客開始接受服務(wù)的時刻為，第個顧客服務(wù)結(jié)束的時刻為，第個顧客與第個顧客之間的時間間隔為，第個顧客的服務(wù)時間為。這樣便容易得到如下關(guān)系：計算機模擬 上述排隊系統(tǒng)的仿真算法步驟如下：（1）將時間區(qū)間0,60等分為60份，即每個時間間

32、隔為1分鐘。（2）應(yīng)用MATLAB隨機數(shù)生產(chǎn)函數(shù)，按照等概率原則，模擬出前分鐘中顧客到達情況： 、；（3）再按照的顧客理發(fā)僅用分鐘,另外的顧客用分鐘的顧客接受服務(wù)時間原則，模擬出顧客 、所對應(yīng)的接受服務(wù)的時間、；（4）確定每個顧客接受服務(wù)的時間。當(dāng)系統(tǒng)剛開始運行時，第一、第二個顧客可以直接接受服務(wù)，而對于以后的顧客，其接受服務(wù)的時刻的判斷取決于三者中最大者。即：；（5）計算顧客的累計等待時間。顧客總的等待時間為，得到：；（6）判斷第個顧客開始接受服務(wù)的時刻為是否超出分鐘。若是，則總得接受服務(wù)的人數(shù)為：,顧客平均等待時間為：；若否，重復(fù)步驟（3）。編寫進行數(shù)字模擬的MATLAB程序如下：% 排隊

33、系統(tǒng)模擬clear allN=60; %模擬階段數(shù)num=0; %完成服務(wù)的人數(shù) %模擬每階段顧客到達情況，當(dāng)r1=1時，有顧客到達，否則，無顧客到達r1=binornd(1,0.5,N,1); c=find(r1); %確定顧客到達的時間m=length(c); %計算0, N內(nèi)到達顧客數(shù)%模擬每位顧客服務(wù)情況，當(dāng)時，服務(wù)時間為5分鐘，否則，服務(wù)時間為8分鐘r2=binornd(1,0.6,m,1); y=8*ones(m,1);y(r2=1)=5;%初始化前兩位顧客服務(wù)結(jié)束時刻b=zeros(m,1); %顧客服務(wù)時間e=zeros(m,1); %顧客服務(wù)結(jié)束時刻w=zeros(m,1);

34、 %顧客等待服務(wù)時間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('序號 到達時間 服務(wù)開始時間 服務(wù)結(jié)束時間 等待服務(wù)時間');for k=1:num fprintf('%2d %6d %

35、6d %6d %6dn',k,c(k),b(k),e(k),w(k)endfprintf('%d分鐘內(nèi)完成服務(wù)的人數(shù)為%dn',N,num);fprintf('顧客平均等待時間為%2.1f分鐘n',mean(w); 運行后的結(jié)果見表12-8，60分鐘內(nèi)完成服務(wù)的人數(shù)為12，顧客平均等待時間為4.4分鐘。表12-8 排隊系統(tǒng)模擬過程序號到達時間服務(wù)開始時間服務(wù)結(jié)束時間等待服務(wù)時間133702471233101217241217225513222696162631107233136882636411092841461310304651161131515520

36、123255592312.2.4連續(xù)型隨機變量的模擬連續(xù)系統(tǒng)模擬研究系統(tǒng)的狀態(tài)隨時間連續(xù)變化的情況。連續(xù)系統(tǒng)雖然研究的是確定性問題，但在實際處理中，是將連續(xù)狀態(tài)變量在時間上進行離散化處理，并由此模擬系統(tǒng)的運行狀態(tài)，即其狀態(tài)由于隨機時間的推進而發(fā)生變化。因此，連續(xù)系統(tǒng)的計算機模擬只能是近似的，只要這種近似達到一定的精度，一般就可以滿足要求。一、連續(xù)型隨機變量模擬的原理處理連續(xù)型隨機變量模擬問題有多種方法，其中反函數(shù)法是最常用的。其方法核心是通過求概率分布的反函數(shù)產(chǎn)生隨機數(shù)。眾所周知隨機變量的概率分布函數(shù)為定義在 區(qū)間的單調(diào)遞增函數(shù)，設(shè)為區(qū)間的均勻隨機變量, 令,只要求出反函數(shù),即為具有概率分布函

37、數(shù)的隨機機數(shù)。由概率論的理論可以證明和相同的概率分布。若的概率密度為，由,是區(qū)間上均勻分布的隨機變量，如果給定區(qū)間上均勻分布的隨機數(shù), 則具有給定分布的隨機數(shù)可由方程中解出：二、連續(xù)型隨機變量模擬的實例應(yīng)用1、鹽水濃度問題例12.7 在水池中有2000m3水，其中含鹽2kg?，F(xiàn)以每分鐘6m3的速度向池中注入含鹽率為0.5kg/m3的鹽水，同時又以每分鐘4m3的速度從水池流出攪拌均勻的鹽水。求出池中含鹽量隨時間的變化規(guī)律，并算出含鹽率達到0.2kg/m3時所用的時間。問題分析 初始時刻水池中含鹽率為0.001kg/m3，注入鹽水的濃度為0.5kg/m3，故水池中鹽水的濃度逐漸升高。模型構(gòu)建與求解

38、 設(shè)t分鐘時刻水池中含有kg的鹽。考察t, t+t時間內(nèi)水池中含鹽量的變化，有：兩邊除以t，且使t0，則可得微分方程模型：在MATLAB命令窗口中輸入如下指令：x=dsolve('Dx=3-2*x/(1000+t)','x(0)=2')運行后輸出結(jié)果：x = 1000+t-998000000/(1000+t)2結(jié)果表明，微分方程的解為最后，由 得容易算出含鹽率達到0.2kg/m3時所用的時間為t184.84（分鐘）計算機仿真 設(shè)初始時刻時間為，時間步長為1分鐘，設(shè)第k分鐘池水中含鹽量為，池水體積為，含鹽率為，假設(shè)在時間區(qū)間k,k+1)池水中的含鹽率不變，則每10

39、分鐘輸出一次體積、含鹽量、含鹽率，直到含鹽率達到0.2 kg/m3，結(jié)束程序，流程圖如圖12-5所示。圖12-5 鹽水濃度隨時間變化模擬流程圖根據(jù)流程圖12-5，編寫MATLAB程序如下：clc, clear all;time=0; %時間初始化volume=2000; %狀態(tài)初始化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);運行后輸出結(jié)果見表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中結(jié)果表明，第185分鐘時，池水中含鹽量剛好超過0.2 kg/m3，這與精確解184.84分鐘十分的接近

42、。2、追擊性問題例12.8 正方形的個頂點各有人。在某一時刻，人同時出發(fā)以勻速，按順時針方向追逐下一人，如果他們始終保持對準(zhǔn)目標(biāo)，則最終按螺旋狀曲線于中心點。試求出這種情況下每個人的行進軌跡。模型建立 建立平面直角坐標(biāo)系，設(shè),，取時間間隔為，計算每一點在各個時刻的坐標(biāo)。設(shè)某點在時刻的坐標(biāo)為，則在時刻的坐為其中，。計算機模擬 設(shè)置初始化條件為，取足夠小的,當(dāng)時算法結(jié)束。對每一個點連接它在各時刻的位置,即得所求運動軌跡。然后根據(jù)上述建立的模型，用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)通過仿真動畫,追擊的過程可以一目了然地表現(xiàn)出來。12.3 系統(tǒng)仿真軟件SIMULINK12.3.1 SIMULINK簡介SIMULINK是MATLAB中一個用來對動態(tài)系統(tǒng)進行建模、仿真和分析的軟件包。SIMULINK 中的“SIMU”一詞表示可用于計算機仿真，而“LINK”一詞表示它能進行系統(tǒng)連接，即把一系列模塊連接起來，構(gòu)成復(fù)雜的系統(tǒng)模型。它支持線性和非線性系統(tǒng)，連續(xù)和離散模型等。SIMULINK可以使用顯示模塊可以在運行仿真時觀察到仿真結(jié)果，還可以改變參數(shù)并立即就可以觀察到

46、它的變化。SIMULINK提供有可自定義的定制模塊庫、模型編輯器和模型求解器。模塊庫的界面為 SIMULINK模塊庫瀏覽器，使用命令simulink或工具條中的快捷鈕可以進入模塊庫瀏覽器。如圖12-7所示，在該瀏覽器中，模塊庫以文件夾形式組織，由SIMULINK通用庫和專業(yè)庫組成，而模塊單元以文件形式保存。圖12-7 Simulink模塊庫瀏覽器如圖12-7所示，SIMILINK模塊通用庫按功能進行分為：Continuous（連續(xù)模塊），discontinuous（非連續(xù)模塊），Discrete（離散模塊），Logic and Bit Operations（邏輯和位操作），Lookup Tab

47、les（查找表模塊），Math Operations（數(shù)學(xué)模塊），Model Verification（模型檢驗?zāi)K），Model-Wide Utilities（模型擴充模塊），Ports&Subsystems（端口和子系統(tǒng)模塊），Signal Attributes（信號屬性模塊），Signal Routing（信號線路模塊），Sinks（接收器模塊），Sources（輸入源模塊）和User Defined Function（用戶自定義函數(shù)模塊）等。在SIMILINK中，各功能模塊的參數(shù)描述都可以由用戶通過該模塊的模塊屬性對話框進行操作給出或修改。在“File”下點擊“NewModel

48、”或單擊模塊庫瀏覽器的新建圖標(biāo)就可以進入模型空間組建模型，如圖12-8所示。在模型編輯器中，可以直接從模塊庫中將模塊用鼠標(biāo)點擊和拖拉過來就可以拼接成所需要的模型，它可以以M文件的形式保存下來。圖12-8 模型編輯器SIMILINK的模型窗口的主菜單有文件、編輯、查看、仿真、格式設(shè)定、工具與幫助七項菜單選項。模型求解器位于模型瀏覽器的菜單“Simulation”中,Simulation(仿真)菜單包括：l start (Ctrl +T)：啟動或暫停仿真。l stop：停止仿真。l Configuration Parameters (Ctrl +E)：設(shè)置仿真參數(shù)l Normal：常規(guī)標(biāo)準(zhǔn)仿真。l

49、 Accelerator：加速仿真。12.3.2 SIMULINK建模一個典型的SIMULINK模型包括如下三種要素： 信號源模塊 被模擬的系統(tǒng)模塊 輸出顯示模塊 系統(tǒng)模塊作為中心模塊是SIMULINK仿真建模所要解決的主要部分；信號源為系統(tǒng)的輸入，它包括常數(shù)信號源函數(shù)信號發(fā)生器（如正弦和階躍函數(shù)波等）和用戶自己在Matlab中創(chuàng)建的自定義信號或Matlab工作間中三種。輸出模塊主要在Sinks庫中。SIMULINK模型并不一定要包含全部的三種要素，在實際應(yīng)用中通?？梢匀鄙倨渲械囊粋€或兩個。例如，若要模擬一個系統(tǒng)偏離平衡位置后的恢復(fù)行為，就可以建立一個沒有輸入而只有系統(tǒng)模塊加一個顯示模塊的模型

50、。在某種情況下，也可以建立一個只有源模塊和顯示模塊的系統(tǒng)。若需要一個由幾個函數(shù)復(fù)合的特殊信號，則可以使用源模塊生成信號并將其送入Matlab工作間或文件中。利用SIMULINK建立物理系統(tǒng)和數(shù)學(xué)系統(tǒng)的仿真模型，關(guān)鍵是對SIMULINK提供的功能模塊進行操作，即用適當(dāng)?shù)姆绞綄⒏鞣N模塊連接在一起。利用SIMULINK建模與仿真的步驟為： 啟動SIMULINK，打開SIMULINK模塊庫； 打開空白模型窗口，用鼠標(biāo)拖拉模塊構(gòu)成模型； 設(shè)置仿真參數(shù)，進行仿真； 輸出仿真結(jié)果，以觀察仿真的結(jié)果。下面以數(shù)學(xué)方程的建模為例，說明SIMULINK的建模過程。例12.9 構(gòu)建攝氏向華氏溫度的變換器分析：借助SI

51、MULINK模塊化框圖實現(xiàn)TF=9/5(TC)+32。 功能框圖如圖12-9所示：圖12-9 方程TF=9/5(TC)+32的功能框圖解：步驟1：模塊定位與選取如表所示。功能需求模塊庫模塊單元攝氏溫度（信號）sourceRamp(斜坡發(fā)生器)信號增益mathgain常量元素mathconstant求和運算mathsum信號顯示sinkScope(示波器)步驟2：按功能框圖配置模塊元件。打開模塊庫瀏覽器，進入模型空間組建模型。直接從模塊庫中將模塊拖拉過來拼接成所需要的模型，如圖12-10所示。圖12-10 模型所需模塊元件圖步驟3: 添加信號線。 如圖12-11所示，共增加4條線路。線路1：信號發(fā)生器信號增益，線路2：信號增益信號求和單元，線路3：常量單元信號求和單元，線路4：信號