2014數(shù)學(xué)建模之計算機仿真

1、數(shù)學(xué)建模之,計算機仿真,天目學(xué)院基礎(chǔ)部 帥昌浩 tel:18067679183,China Undergraduate Mathematical Contest in Modeling,概述,計算機科學(xué)技術(shù)的迅猛發(fā)展，給許多學(xué)科帶來了巨大的影響．計算機不但使問題的求解變得更加方便、快捷和精確，而且使得解決實際問題的領(lǐng)域更加廣泛．計算機適合于解決那些規(guī)模大、難以解析化以及不確定的數(shù)學(xué)模型．例如對于一些帶隨機因素的復(fù)

2、雜系統(tǒng)，用分析方法建模常常需要作許多簡化假設(shè)，與面臨的實際問題可能相差甚遠，以致解答根本無法應(yīng)用，這時仿真幾乎成為人們的唯一選擇．在歷屆的美國和中國大學(xué)生的數(shù)學(xué)建模競賽（MCM）中，學(xué)生們經(jīng)常用到計算機仿真方法去求解、檢驗等．計算機仿真(computer simulation)是建模過程中較為重要的一類方法．,計算機仿真的基本概念,計算機仿真，是根據(jù)已知的信息和知識（如數(shù)學(xué)、物理規(guī)律），利用計算機模擬現(xiàn)實情況和系統(tǒng)的演變過程，發(fā)現(xiàn)新的知

3、識或規(guī)律，從而解決問題的一種方法,計算機仿真的基本概念,獨立于理論研究與實驗研究的認識世界的第三中方法,計算機仿真的基本概念,計算機仿真的特點,代價小、時間短、可重復(fù)、參數(shù)設(shè)置靈活,計算機仿真可以解決以下5類問題:,(1) 難以用數(shù)學(xué)公式表示的系統(tǒng),或者沒有建立和求解的有效方法.(2)雖然可以用解析的方法解決問題,但數(shù)學(xué)的分析與計算過于復(fù)雜,這時計算機仿真可能提供簡單可行的求解方法.(3)希望能在較短的時間內(nèi)觀察到系統(tǒng)發(fā)展的全過程

4、,以估計某些參數(shù)對系統(tǒng)行為的影響.(4)難以在實際環(huán)境中進行試驗和觀察時,計算機仿真是唯一可行的方法,如太空飛行的研究.(5)需要對系統(tǒng)或過程進行長期運行比較,從大量方案中尋找最優(yōu)方案.,計算機仿真的分類,計算機仿真在計算機中運行實現(xiàn),不怕破壞,易修改,可重用,安全經(jīng)濟,不受外界條件和場地空間的限制.仿真分為靜態(tài)仿真(static simulation)和動態(tài)仿真(dynamic simulation).數(shù)值積分中的蒙特卡洛方法（

5、統(tǒng)計模擬方法）是典型的靜態(tài)仿真．動態(tài)仿真又分為連續(xù)系統(tǒng)仿真和離散系統(tǒng)仿真．連續(xù)系統(tǒng)是指狀態(tài)變量隨著時間連續(xù)變化的系統(tǒng)，例如傳染病的檢測與預(yù)報系統(tǒng).離散系統(tǒng)是指系統(tǒng)狀態(tài)變量只在有限的時間點或可數(shù)的時間點上發(fā)生變化的系統(tǒng),例如排隊系統(tǒng).,仿真系統(tǒng)，必須設(shè)置一個仿真時鐘(simulate clock)，它能將時間從一個時刻向另一個時刻進行推進，并且能隨時反映系統(tǒng)時間的當(dāng)前值．其中，模擬時間推進方式有兩種:時間步長法(均勻間隔時間推進法,連續(xù)系

6、統(tǒng)常用)和事件步長法(下次事件推進法,離散系統(tǒng)常用) ．,主要內(nèi)容,一: 準(zhǔn)備知識:隨機數(shù)的產(chǎn)生二:隨機變量的模擬三:連續(xù)系統(tǒng)的模擬-時間步長法四: 離散系統(tǒng)的模擬-事件步長法●五：蒙特卡洛方法,一: 準(zhǔn)備知識:隨機數(shù)的產(chǎn)生,由于仿真研究的實際系統(tǒng)要受到多種隨機因素的作用和影響,在仿真過程中必須處理大量的隨機因素.要解決此問題的前提是確定隨機變量的類型和選擇合適的隨機數(shù)產(chǎn)生的方法.對隨機現(xiàn)象進行模擬,實質(zhì)是要給出隨機變量的模擬

7、,也就是說要利用計算機隨機產(chǎn)生一系列數(shù)值,使它們服從一定的概率分布,稱這些數(shù)值為隨機數(shù).最基本,最常用的是(0,1)區(qū)間內(nèi)均勻分布的隨機數(shù).其他分布的隨機數(shù)均可利用它來產(chǎn)生.,1:產(chǎn)生模擬隨機數(shù)的計算機命令,在MATLAB中,可以直接產(chǎn)生滿足各種分布的隨機數(shù),命令如下:常見的分布函數(shù) MATLAB語句 均勻分布U[0,1] R=rand(m,n) 均勻分布U[a,b]

8、 R=unifrnd(a,b,m,n) 指數(shù)分布E(λ) R = exprnd(λ, m , n)正態(tài)分布N(mu,sigma) R=normrnd(mu,sigma,m,n) 標(biāo)準(zhǔn)正態(tài)分布N(0,1) R=randn(m,n) 二項分布B(n,p) R=binornd(n,p,m,n) 泊松分布 P(λ ) R=poi

9、ssrnd (λ, m , n)以上語句均產(chǎn)生m× n 的矩陣．,2:案例分析,例1: unifrnd(2,3)unifrnd(1,32,1,4)normrnd(1,2) normrnd(1,2,2,3)rand(2,3)randn(2,3),2:案例分析,ans =2.8132ans =1.3057 5.3056 7.2857 7.1604ans = 0.2527ans =

10、 2.7429 0.0219 2.7759 2.2756 0.0992 -0.9560ans = 0.6038 0.1988 0.7468 0.2722 0.0153 0.4451ans = -0.0945 -1.3089 -0.2440 -0.2141 0.8248 -0.1778>>,2:案例分析,例2:敵空戰(zhàn)部

11、隊對我方港口進行空襲,其到達規(guī)律服從泊松分布,平均每分鐘到達4架飛機.(1)模擬敵機在3分鐘內(nèi)到達目標(biāo)區(qū)域的數(shù)量,以及在第1,2,3分鐘內(nèi)各到達幾架飛機;(2)模擬在3分鐘內(nèi)每架飛機的到達時刻.分析:(1) n1=poissrnd(4), n2=poissrnd(4), n3=poissrnd(4),n=n1+n2+n3 (2) 由排隊論知識,敵機到達規(guī)律服從泊松分布等價于敵機到達港口的間隔時間服從參數(shù)為1/4的指數(shù)分布,故

12、可由指數(shù)分布模擬每架飛機的到達時刻.,2:案例分析,cleart=0;j=0; %到達的飛機數(shù) while t<3 j=j+1 t=t+exprnd(1/4)end,二:隨機變量的模擬,利用均勻分布的隨機數(shù)可以產(chǎn)生具有任意分布的隨機變量的樣本,從而可以對隨機變量的取值情況進行模擬.1 連續(xù)型隨機變量的模擬具有給定分布的連續(xù)型隨機變量可以利用在區(qū)間(0

13、,1)上均勻分布的隨機數(shù)來模擬,最常用的方法是逆變換法.結(jié)論:若隨機變量Y有連續(xù)的分布函數(shù)F(y), 則Z與Y有相同的分布.,1 連續(xù)型隨機變量的模擬,若已知Y的概率密度為如果給定區(qū)間(0,1)上均勻分布的隨機數(shù) ,則具有給定分布Y的隨機數(shù) 可由方程 解出.例 :模擬服從參數(shù)為 的指數(shù)分布時,由

14、 可得,2 離散型隨機變量的模擬,設(shè)隨機變量X的分布律為:,有相同的發(fā)生的概率.因此我們可以用隨機變量R落在小區(qū)間內(nèi)的情況來模擬離散的隨機變量X的取值情況.,2 離散型隨機變量的模擬,例 3 :隨機變量 表示每分鐘到達銀行柜臺的顧客數(shù).X的分布列見下表,試模擬10分鐘內(nèi)顧客到達柜臺的情況. 表1

15、10分鐘內(nèi)顧客到達柜臺的情況 Xk 0 1 2 pk 0.4 0.3 0.3分析:因為每分鐘到達柜臺的人數(shù)是隨機的,所以可用計算機隨機生成一組(0,1)的數(shù)據(jù),由X的概率分布情況,可認為隨機數(shù)在(0,0.4)范圍內(nèi)時沒有顧客光顧,在[0.4,0.7)時,有一個顧客光顧,在[0.7,1)時,有兩個顧客光顧. 從而有MATLAB程序:,2 離散型隨機變量的模擬,r

16、=rand(1,10);for i=1:10; if r(i)<0.4 n(i)=0; elseif 0.4<=r(i)&r(i)<0.7 n(i)=1; else n(i)=2; end;

17、 end r n,三:連續(xù)系統(tǒng)的模擬-時間步長法,對連續(xù)系統(tǒng)的計算機模擬是近似地獲取系統(tǒng)狀態(tài)在一些離散時刻點上的數(shù)值．在一定假設(shè)條件下，利用數(shù)學(xué)運算模擬系統(tǒng)的運行過程．連續(xù)系統(tǒng)模型一般是微分方程，它在數(shù)值模擬中最基本的算法是數(shù)值積分算法．例如有一系統(tǒng)可用微分方程來描述: 已知輸出量y的初始條件 ，現(xiàn)在要求出輸出量y隨時間變化的過程y(

18、t)。,1 理論介紹,最直觀的想法是：首先將時間離散化，令 ，稱為第k步的計算步距（一般是等間距的），然后按以下算法計算狀態(tài)變量在各時刻 上的近似值: 其中初始點 按照這種作法即可求出整個的曲線．這種最簡單的數(shù)值積分算法稱為歐拉法．除此之外，還有其他一些算法． 􀀢,1 理

19、論介紹,因此，連續(xù)系統(tǒng)模擬方法是：首先確定系統(tǒng)的連續(xù)狀態(tài)變量，然后將它在時間上進行離散化處理，并由此模擬系統(tǒng)的運行狀態(tài)．模擬過程分為許多相等的時間間隔,時間步長的長度可以根據(jù)實際問題分別取為秒,分,小時,天等.仿真時鐘按時間步長等距推進,每次推進都要掃描系統(tǒng)中所有活動,按預(yù)定計劃和目標(biāo)進行分析,計算,記錄系統(tǒng)狀態(tài)的變化,直到滿足某個終止條件為止.,,例:追逐問題1．問題提出假設(shè)正方形ABCD的四個頂點處各站一人．在某一時刻，四人同時

20、以勻速v沿順時針方向追逐下一個人，并且在任意時刻他們始終保持追逐的方向是對準(zhǔn)追逐目標(biāo)，例如，A追逐B(yǎng)，任意時刻A始終向著B追．可以證明四人的運動軌跡將按螺旋曲線狀匯合于中心O． 怎樣證明呢？有兩種證明方法．一是分別求出四人的運動軌跡曲線解析式，求證四條曲線在某時刻相交于一點．另一方法則是用計算機模擬將四人的運動軌跡直觀地表示在圖形上.,2 應(yīng)用舉例,2．建立模型及模擬方法,模擬步驟： 1）建立平面直角坐標(biāo)系． 2）以時間間隔Δt進行

21、采樣，在每一時t計算每個人在下一時t+Δt時的坐標(biāo)． 3）不妨設(shè)甲的追逐對象是乙，在時間t時，甲,2．建立模型及模擬方法,4)選取足夠小的 ，模擬到 時為止5）連接四人在各時刻的位置，就得到所求的軌跡．根據(jù)以上模擬步驟，編出MATLAB程序(simu2.m)如下：,%取v=1,t=12,A,B,C,D點的坐標(biāo)分另為（0，10），（10，10），（10，0），(0, 0)v=1;dt=0.05;

22、d=20;x=[0 0 0 10 10 10 10 0];x(9)=x(1);x(10)=x(2);holdaxis('equal')axis([0 10 0 10]);,3 MATLAB實現(xiàn),for k=1:2:7plot(x(k),x(k+1),'.' )endwhile(d>0.1)for i=1:2:7d=sqrt((x(i)-x(i+1))^2+(x(i+1)-x(i

23、+3))^2);x(i)=x(i)+v*dt*(x(i+2)-x(i))/d;x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1))/d;plot(x(i),x(i+1),'.')endx(9)= x(1);x(10)= x(2);endhold,3 MATLAB實現(xiàn),3 MATLAB實現(xiàn),例9.6 水池含鹽量問題,某水池有2000m3水,其中含鹽2kg,以6m3/min的速率向水池內(nèi)注入含鹽

24、為0.5kg/m3的鹽水,同時又以4m3/min的速率從水池流出攪拌均勻的鹽水.試用計算機仿真該水池內(nèi)鹽水的變化過程,并每隔10min計算水池中水的體積,含鹽量,含鹽率.欲使池中鹽水含鹽率達到0.2kg/m3,需經(jīng)過多長時間?分析:這是一個連續(xù)系統(tǒng),首先要將系統(tǒng)離散化,在一些離散點上進行考察,這些離散點的間隔就是時間步長.可取步長為1min,即隔1min考察一次系統(tǒng)的狀態(tài),并相應(yīng)地記錄和分析.在注入和流出的作用下,池中水的體積與含

25、鹽量,含鹽率均隨時間變化,初始時刻含鹽率為0.001kg/m3,以后每分鐘注入含鹽率為0.5kg/m3的水6m3,流出混合均勻的鹽水為4m3,當(dāng)池中水的含鹽率達到0.2kg/m3時,仿真過程結(jié)束.,例9.6 水池含鹽量問題,記T時刻的體積為w(m3)，水的含鹽量為s(kg)，水的含鹽率為r=s/w(kg/m3)，每隔1min池水的動態(tài)變化過程如下：每分鐘水的體積增加6-4=2(m3)；每分鐘向池內(nèi)注入鹽6×0.5=3(kg)；

26、每分鐘向池外流出鹽4×r(kg)；每分鐘池內(nèi)增加鹽3-4×r(kg).本例還可以用微分方程建立數(shù)學(xué)模型，并求出它的解析解，這個解析解就是問題的精確解，有興趣的讀者可以按照這個思路求出該問題的精確解，考察相應(yīng)時刻精確解與仿真解的差異，還可以進一步調(diào)整仿真過程的時間步長，通過與精確解的比較來研究時間步長的大小對仿真度的影響。,MATLAB實現(xiàn),clearh=1; %時間步長為1s0=2; %初始

27、含鹽2kgw0=2000; %初始水池有水2000m3r0=s0/w0; %初始濃度s(1)=s0+0.5*6*h-4*h*r0; %一分鐘后的含鹽量w(1)=w0+2*h; %一分鐘后水池中的鹽水體積r(1)=s(1)/w(1); %一分鐘后的濃度 t(1)=h;y(1)=(2000000+3000000*h+3000*h^2+h^3)/(1000+h)

28、^2;for i=2:200 t(i)=i*h; s(i)=s(i-1)+0.5*6*h-4*h*r(i-1); %第i步后的含鹽量 w(i)=w(i-1)+2*h; %第i步后的鹽水體積 r(i)=s(i)/w(i); %第i步后的鹽水濃度 y(i)=(2000000+3000000*t(i)+3000*t(i)^2+

29、t(i)^3)/(1000+t(i))^2; m=floor(i/10);,MATLAB實現(xiàn),if i/10-m0.2 %若第i步后的鹽水濃度大于0.2 t02=i*h; r02=r(i); break endend[t02,r02][10*tm',sm',rm'] % '表示逆subplot(1,2,

30、1),plot(t,s,'blue');hold onsubplot(1,2,2),plot(t,y,'red');,四: 離散系統(tǒng)的模擬-事件步長法,離散系統(tǒng)(discrete system)是指系統(tǒng)狀態(tài)只在有限的時間點或可數(shù)的時間點上有隨機事件驅(qū)動的系統(tǒng)．例如排隊系統(tǒng)（queue system），顯然狀態(tài)量的變化只是在離散的隨機時間點上發(fā)生．假設(shè)離散系統(tǒng)狀態(tài)的變化是在一個時間點上瞬間完成的． 常

31、用的是事件步長法（下次事件推進法）．其過程是：置模擬時鐘的初值為0，跳到第一個事件發(fā)生的時刻，計算系統(tǒng)的狀態(tài)，產(chǎn)生未來事件并加入到隊列中去，跳到下一事件，計算系統(tǒng)狀態(tài)，……，重復(fù)這一過程直到滿足某個終止條件為止．,例9.7 收款臺前的排隊過程,假設(shè):(1) 顧客到達收款臺是隨機的,平均時間間隔為0.5min,即間隔時間服從lamda=2的指數(shù)分布.(2) 對不同的顧客收款和裝袋的時間服從正態(tài)分布N(1,1/3).試模擬20位顧客到

32、收款臺前的排隊情況,我們關(guān)心的問題是每個顧客的平均等待時間,隊長及服務(wù)員的工作效率,例9.7 收款臺前的排隊過程,分析:單服務(wù)臺結(jié)構(gòu)的排隊系統(tǒng)有兩類原發(fā)事件即到來和離去,顧客到來的后繼事件是顧客接受服務(wù),顧客離去的后繼事件是服務(wù)臺尋找服務(wù),這四類事件各自的子程序框圖如圖1所示.假設(shè):t(i)為第i位顧客到達時刻;t2(i)為第i位顧客受到服務(wù)的時間(隨機變量);T(i)為第i位顧客離去時刻.將第i位顧客到達作為件事發(fā)生:t(i+1)

33、-t(i)=r(i)(隨機變量);平衡關(guān)系:當(dāng)t(i+1) ≥ T(i)時, T(i+1)= t(i+1) + t2(i+1);否則, T(i+1)= T(i) + t2(i+1).,圖1 到來事件子程序,,,,系統(tǒng)人數(shù)+1,產(chǎn)生下一個顧客到來時刻,調(diào)用接受服務(wù)事件的程序,,圖 2 離去事件子程序,,系統(tǒng)人數(shù)-1,,,置服務(wù)臺”閑”,已服務(wù)人數(shù)+1,調(diào)用尋找服務(wù)事件子程序,,圖3 :接受服務(wù)事件子程序,,服務(wù)臺空,,,置服務(wù)臺”忙”,

34、產(chǎn)生服務(wù)結(jié)束時刻,登記到事件表,,,,排隊人數(shù)+1,否,是,,,圖4 :尋找服務(wù)事件子程序,,排隊人數(shù)≥1,,,置服務(wù)臺”忙”,產(chǎn)生服務(wù)結(jié)束時刻,排隊人數(shù)-1,,,,排隊人數(shù)+1,否,是,,,五：蒙特卡洛方法,在用傳統(tǒng)方法難以解決的問題中，有很大一部分可以用概率模型進行描述．由于這類模型含有不確定的隨機因素，分析起來通常比確定性的模型困難．有的模型難以作定量分析，得不到解析的結(jié)果，或者是雖有解析結(jié)果，但計算代價太大以至不能使用．在這種情

35、況下，可以考慮采用Monte Carlo方法。Monte Carlo方法是計算機模擬的基礎(chǔ)，它的名字來源于世界著名的賭城——摩納哥的蒙特卡洛，其歷史起源于1777年法國科學(xué)家蒲豐提出的一種計算圓周π的方法——隨機投針法，即著名的蒲豐投針問題。,五：蒙特卡洛方法,Monte Carlo方法的基本思想是首先建立一個概率模型，使所求問題的解正好是該模型的參數(shù)或其他有關(guān)的特征量．然后通過模擬一統(tǒng)計試驗，即多次隨機抽樣試驗（確定m和n），統(tǒng)計出某

36、事件發(fā)生的百分比．只要試驗次數(shù)很大，該百分比便近似于事件發(fā)生的概率．這實際上就是概率的統(tǒng)計定義．利用建立的概率模型，求出要估計的參數(shù)．蒙特卡洛方法屬于試驗數(shù)學(xué)的一個分支．,五：蒙特卡洛方法,蒙特卡洛方法適用范圍很廣泛，它既能求解確定性的問題，也能求解隨機性的問題以及科學(xué)研究中的理論問題．例如利用蒙特卡洛方法可以近似地計算定積分，即產(chǎn)生數(shù)值積分問題．,蒙特卡洛法求圓周率,clearn=50000X=rand(n,1);Y=rand(

37、n,1); k=0;for i=1:n; if X(i)^2+Y(i)^2<=1 k=k+1; endend4*k/n,普豐投針求圓周率的M程序,function pi_value=pinpi(k,a,l) %k投針次數(shù)%a線距%l針長（必須小于等于d）%pi_value返回pi值m=0; if l>a; fprintf('error：針長必須小于等于%d\

38、n',a); fprintf('請重新調(diào)用函數(shù)pinpi(k,d,l) \n'); pi_value=0;else for i=1:k if a*rand(1)<=l*sin(pi*rand(1)) m=m+1; end end p=m/k;pi_value=2*l/(a*p); fprintf('投針法求得pi=%d\n',pi_valu

39、e);end,,format long g>> pinpi(100000,4,3)投針法求得pi=3.143797e+000ans = 3.14379728795087,任意曲邊梯形面積的近似計算,一個古老的問題：用一堆石頭測量一個水塘的面積．應(yīng)該怎樣做呢？測量方法如下：假定水塘位于一塊面積已知的矩形農(nóng)田之中．如圖8．2所示．隨機地向這塊農(nóng)田扔石頭使得它們都落在農(nóng)田內(nèi)．被扔到農(nóng)田中的石頭可能濺

40、上了水，也可能沒有濺上水，估計被“濺上水的”石頭量占總的石頭量的百分比．試想如何利用這估計的百分比去近似計算該水塘面積？,任意曲邊梯形面積的近似計算,任意曲邊梯形面積的近似計算,結(jié)合圖8．2中的圖形(1)分析，只要已知各種參數(shù)及函數(shù)（a，b，H，f(x)），有以下兩種方法可近似計算水塘面積． 1．隨機投點法 1）賦初值：試驗次數(shù)n=0，成功次數(shù)m=0；規(guī)定投點試驗的總次數(shù)N； 2）隨機選擇m個數(shù)對xi,yi,1<i<m

41、，其中ya<xi<b, 0<yi<H，置 n＝n＋l； 3）判斷n≤N，若是，轉(zhuǎn)4，否則停止計算；4）判斷條件 (表示一塊濺水的石頭)是否成立，若成立則置m=m+1，轉(zhuǎn)2，否則轉(zhuǎn)2；5）計算水塘面積的近似值S=H×(b-a)×m/N．,任意曲邊梯形面積的近似計算,2．平均值估計法1）產(chǎn)生[a,b]區(qū)間的均勻隨機數(shù)x,i=1,2…,N􀀢

42、;2) 計算f(xi）， i=1,2…,N3）計算該方法的特點是估計函數(shù)f(x)在[a,b]上的平均值，面積近似等于該平均值乘以(b-a)．,例9 庫存問題,在物資的供應(yīng)過程中,由于到貨和銷售不可能做到同步同量,故總要保持一定的庫存儲備.如果庫存過多,就會造成積壓浪費以及保管費的上升;如果庫存過少,會造成缺貨.如何選擇庫存和訂貨策略,就是一個需要研究的問題.庫存問題有多種類型,一般都比較復(fù)雜,下面討論一種簡單的情形. 某電動車行

43、的倉庫管理人員采取一種簡單的訂貨策略,當(dāng)庫存降低到P輛電動車時就向廠家,例9 庫存問題,訂貨,每次訂貨Q輛,如果某一天的需求量超過了庫存量,商店就有銷售損失和信譽損失,但如果庫存量過多,就會導(dǎo)致資金積壓和保管費增加.若現(xiàn)在有如下表所示的兩種庫存策略,試比較選擇一種策略以使總費用最少. 表 訂貨方案 方案 重新訂貨點P/輛 重新訂貨量Q/輛 方案1 125

44、 150 方案2 150 250,例9 庫存問題,這個問題的已知條件是:(1) 從發(fā)出訂貨到收到貨物需隔3天.(2) 每輛電動車保管費為0.50元/天,每輛電動車的缺貨損失為1.60元/天,每次的訂貨費為75元.(3)每天電動車需求量是0到99之間均勻分布的隨機數(shù).(4)原始庫存為110輛,假設(shè)第一天沒有發(fā)出訂貨.,例9 庫存問題,分析: 這一問題用解析法討論比較麻煩,但用計算

45、機按天仿真?zhèn)}庫貨物的變動情況卻很方便.我們以30天為例,依次對這兩種方案進行仿真,最后比較各方案的總費用,從而就可以作出決策.計算機仿真時的工作流程是早上到貨,全天銷售,晚上訂貨,輸入一些常數(shù)和初始數(shù)據(jù)后,以一天為時間步長進行仿真.首先檢查這一天是否為預(yù)訂到貨日期,如果是,則原有庫存量加Q,并把預(yù)定到貨量清為0;如果不是,則庫存量不變.接著用計算機中的隨機函數(shù)仿真隨機需求量,若庫存量大于需求量,則新的庫存量減去需求量;反之,則新庫存量

46、變?yōu)?,并且要在總費用上加缺貨損失,然后檢查實際,例9 庫存問題,庫存量加上預(yù)定到貨量是否小于重新訂貨點P,如果是,則需要重新訂貨,這時就加一次訂貨費.如此重復(fù)運行30天,即可得所需費用總值.由此比較這兩種方案的總費用,即得最好方案.,,cleardays=30;P=[125,150];Q=[150,250];cost=[0,0];arrivalinterval=2;storagefee=0.5;lossfee=1.6;bo

47、okfee=75;storage0=110;booknumber=0;arrivedate=0;nr=rand(days,1);for i=1:2 storage(1)=storage0; n=round(99*nr(1)); sale=n; remain=storage(1)-n; if remain<=P(i); booknumber=Q(i

48、); arrivedate=4; orderfee=bookfee; else orderfee=0; end storage(1)=remain; cost(i)=cost(i)+remain*storagefee+orderfee; for j=2:days dh=j;

49、 if abs(dh-arrivedate)<0.01 storage(j)=storage(j-1)+booknumber; booknumber=0;,,arrivedate=j; else storage(j)=storage(j-1);

50、 end n=round(99*nr(j)); if storage(j)>=n sale=n; remain=storage(j)-n; shortagenumber=0; else

51、 sale=storage(j); remain=0; shortagenumber=n-storage(j); end storage(j)=remain; if remain+booknumber<=P(i);

52、 booknumber=Q(i); arrivedate=dh+arrivalinterval; orderfee=bookfee; else orderfee=0; end cost(i)=cost(i)+remain*sto

53、ragefee+shortagenumber*lossfee+orderfee; end; mincost=min(cost); end cost mincost,例9.10 趕火車過程仿真,一列火車從A站經(jīng)B站開往C站,某人每天趕往B站乘這趟火車.已知火車從A站到B站運行時間為均值30min,標(biāo)準(zhǔn)差為2min的正態(tài)隨機變量.火車大約在下午1點離開A站,離開時刻的頻率分布見表1.這

54、個人到達B站時的頻率分布表見表2.用計算機仿真火車開出,火車到達B站,這個人到達B站的情況,并給出能趕上火車的仿真結(jié)果. 表1 離開時刻的頻率分布出發(fā)時刻T 1:00 1:05 1:10頻率 0.7 0.2 0.1,例9.10 趕火車過程仿真,表2 人到達B站時的頻率分布 到達時刻T 1:28 1:30 1:32 1:34 頻

55、率: 0.3 0.4 0.2 0.1分析:引入以下變量:T1為火車從A站開出的時刻;T2為火車從A站運行到B站所需要的時間;T3為此人到達B站的時刻.則有表3. 表3 T1, T2的頻率分布T1,/min 0 5 10 T2/min 28 30 32 34 p 0.7 0.2 0.1 0.3 0.4 0.2 0.1,例9.10 趕火車過程

56、仿真,顯然,這位旅客要趕上火車的條件是T3< T1+ T2,可以通過計算機模擬出這三個時間,再檢驗是否滿足T3< T1+ T2.如果滿足,即能夠趕上火車,若不然,則不能趕上火車.,火車運行時間的仿真程序,x=randn(10000,1);for i=1:10000 y(i)=30+2*x(i);end,開車時間的仿真程序,s1=0; s2=0; s3=0; x=rand(10000,1);

57、 for i=1:10000 if x(i)0.9 s3=s3+1; end end[s1/10000, 1-s1/10000-s3/10000,s3/10000],人到達時刻的仿真程序,s1=0; s2=0; s3=0;s4=0; x=rand(10000,1); for i=1:10000 if x(i)<

58、;0.3 s1=s1+1; elseif x(i)<0.7 s2=s2+1; else if x(i)<0.9 s3=s3+1; else s4=s4+1; end

59、 endend[s1/10000, s2/10000,s3/10000,s4/10000],趕上火車的仿真程序,s=0;x1=rand(10000,1);x2=rand(10000,1);x3=randn(10000,1);for i=1:10000if x1(i)<0.7T1=0;elseif x1(i)<0.9T1=5; else T1=10;

60、 endT2=30+2*x3(i);if x2(i)<0.3 T3=28;elseif x2(i)<0.7 T3=30; elseif x2(i)<0.9 T3=32; else T3=34; endif T3<T1+T2