計算機模擬MATLAB實現課件

《計算機模擬MATLAB實現課件》由會員分享，可在線閱讀，更多相關《計算機模擬MATLAB實現課件（36頁珍藏版）》請在裝配圖網上搜索。

1、實驗目的實驗目的實驗內容實驗內容學習計算機模擬的基本過程與方法。學習計算機模擬的基本過程與方法。1 1、模擬的概念。、模擬的概念。4 4、實驗作業、實驗作業。3、計算機模擬實例。、計算機模擬實例。2、產生隨機數的計算機命令。、產生隨機數的計算機命令。連續系統模擬實例:追逐問題追逐問題離散系統模擬實例:排隊問題排隊問題用蒙特卡洛法解非線性規劃問題用蒙特卡洛法解非線性規劃問題返回計算機模擬實例計算機模擬實例模擬的概念模擬的概念 模擬就是利用物理的、數學的模型來類比、模仿現實系統及其演變過程，以尋求過程規律的一種方法。模擬的基本思想是建立一個試驗模型，這個模型包含所研究系統的主要特點通過對這個實驗模

2、型的運行，獲得所要研究系統的必要信息模擬的方法模擬的方法1、物理模擬：對實際系統及其過程用功能相似的實物系統去模仿。例如，軍事演習、船艇實驗、沙盤作業等。物理模擬通?；ㄙM較大、周期較長，且在物理模型上改變系統結構和系數都較困難。而且，許多系統無法進行物理模擬，如社會經濟系統、生態系統等。在實際問題中，面對一些帶隨機因素的復雜系統，用分析方法建模常常需要作許多簡化假設，與面臨的實際問題可能相差甚遠，以致解答根本無法應用。這時，計算機模擬幾乎成為唯一的選擇。在一定的假設條件下，運用數學運算模擬系統的運行，稱為數學模擬?，F代的數學模擬都是在計算機上進行的，稱為計算機模擬。2、數學模擬 計算機模擬可以

3、反復進行，改變系統的結構和系數都比較容易。蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一種應用隨機數來進行計算機模擬的方法此方法對研究的系統進行隨機觀察抽樣，通過對樣本值的觀察統計，求得所研究系統的某些參數例例1 1在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進行干擾和破壞為躲避我方打擊，敵方對其陣地進行了偽裝并經常變換射擊地點 經過長期觀察發現，我方指揮所對敵方目標的指示有50是準確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人 現在希望能用某種方式把我方將要對敵人實施的20次打擊結

4、果顯現出來，確定有效射擊的比率及毀傷敵方火炮的平均值。分析分析：這是一個概率問題，可以通過理論計算得到相應的概率和期望值.但這樣只能給出作戰行動的最終靜態結果，而顯示不出作戰行動的動態過程.為了能顯示我方20次射擊的過程，現采用模擬的方式。需要模擬出以下兩件事：1.問題分析問題分析2 2 當指示正確時，我方火力單位的射擊結果情況當指示正確時，我方火力單位的射擊結果情況1 1 觀察所對目標的指示正確與否觀察所對目標的指示正確與否模擬試驗有兩種結果，每一種結果出現的概率都是1/2 因此，可用投擲一枚硬幣的方式予以確定可用投擲一枚硬幣的方式予以確定，當硬幣出現正面時為指示正確，反之為不正確 模擬試驗

5、有三種結果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6)這時可用投擲骰子的方法來確定可用投擲骰子的方法來確定：如果出現的是、三個點：則認為沒能擊中敵人；如果出現的是、點：則認為毀傷敵人一門火炮；若出現的是點：則認為毀傷敵人兩門火炮2.符號假設符號假設i：要模擬的打擊次數；k1：沒擊中敵人火炮的射擊總數；k2：擊中敵人一門火炮的射擊總數；k3：擊中敵人兩門火炮的射擊總數E：有效射擊比率；E1：20次射擊平均每次毀傷敵人的火炮數3.模擬框圖模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數?k1=k1+1k2=k

6、2+1k3=k3+1k1=k1+1i20?E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20停止硬幣正面?YNNY1，2，34，564.模擬結果模擬結果消滅敵人火炮數消滅敵人火炮數試驗試驗序號序號投硬幣投硬幣結果結果指示指示正確正確指示指示不正確不正確擲骰子擲骰子結果結果正正反正正反正正反反消滅敵人火炮數消滅敵人火炮數試驗試驗序號序號投硬幣投硬幣結果結果指示指示正確正確指示指示不正確不正確擲骰子擲骰子結果結果正反正反正正正正反正從以上模擬結果可計算出：E=7/20=0.35 20322041201301E=0.55.理論計算理論計算設：觀察所對目標指示正確確觀察所對

7、目標指示不正10jA0：射中敵方火炮的事件；A1：射中敵方一門火炮的事件；A2：射中敵方兩門火炮的事件則由全概率公式：E=P(A0)=P(j=0)P(A0j=0)+P(j=1)P(A0j=1)=25.02121021P(A1)=P(j=0)P(A1j=0)+P(j=1)P(A1j=1)=613121021P(A2)=P(j=0)P(A2j=0)+P(j=1)P(A2j=1)=1216121021E1=33.012126116.結果比較結果比較 理論計算和模擬結果的比較分類項目無效射擊有效射擊平均值模擬理論 返回 雖然模擬結果與理論計算不完全一致，但它卻能更加真實地表雖然模擬結果與理論計算不完全

8、一致，但它卻能更加真實地表達實際戰斗動態過程達實際戰斗動態過程 用蒙特卡洛方法進行計算機模擬的步驟用蒙特卡洛方法進行計算機模擬的步驟：1 設計一個邏輯框圖，即模擬模型這個框圖要正確反映系統各部分運行時的邏輯關系。2 模擬隨機現象可通過具有各種概率分布的模擬隨機數來模擬隨機現象產生模擬隨機數的計算機命令產生模擬隨機數的計算機命令 在Matlab軟件中，可以直接產生滿足各種分布的隨機數，命令如下：2產生mn階，均勻分布的隨機數矩陣：rand(m,n)rand(m,n)產生一個，均勻分布的隨機數：randrand1產生mn階a，b均勻分布U（a，b）的隨機數矩陣：unifrnd(a,b,m,n)un

9、ifrnd(a,b,m,n)產生一個a，b均勻分布的隨機數：unifrnd(a,b)unifrnd(a,b)當只知道一個隨機變量取值在（a，b）內，但不知道（也沒理由假設）它在何處取值的概率大，在何處取值的概率小，就只好用U（a，b）來模擬它。例 1的計算機模擬3.產生 mn 階均值為，方差為的正態分布的隨機數矩陣：n no or rm mr rn nd d (,m,n)產生一個均值為，方差為的正態分布的隨機數：normrnd(,)To Matlab（rnd）當研究對象視為大量相互獨立的隨機變量之和，且其中每一種變量對總和的影響都很小時，可以認為該對象服從正態分布。機械加工得到的零件尺寸的偏差

10、、射擊命中點與目標的偏差、各種測量誤差、人的身高、體重等，都可近似看成服從正態分布。4產生 mn 階期望值為的指數分布的隨機數矩陣：e ex xp pr rn nd d(,m,n)若連續型隨機變量X的概率密度函數為 其中 0為常數，則稱X服從參數為 的指數分布指數分布。000)(xxexft指數分布的期望值為 1排隊服務系統中顧客到達率為常數時的到達間隔、故障率為常數時零件的壽命都服從指數分布。指數分布在排隊論、可靠性分析中有廣泛應用。注意：注意：Matlab中，產生參數為 的指數分布的命令為exprnd()1例例 顧客到達某商店的間隔時間服從參數為顧客到達某商店的間隔時間服從參數為0.10.

11、1的指數分布的指數分布 指數分布的均值為指數分布的均值為1/0.1=101/0.1=10。指兩個顧客到達商店的平均間隔時間是指兩個顧客到達商店的平均間隔時間是1010個單位時間個單位時間.即平均即平均1010個個單位時間到達單位時間到達1 1個顧客個顧客.顧客到達的間隔時間可用顧客到達的間隔時間可用exprnd(10)exprnd(10)模擬。模擬。設離散型隨機變量X的所有可能取值為0,1,2,且取各個值的概率為其中 0為常數，則稱X服從參數為 的帕松分布帕松分布。,2,1,0,!)(kkekXPk5產生 mn 階參數為的帕松分布的隨機數矩陣：poissrnd poissrnd(,m,n)帕松

12、分布在排隊系統、產品檢驗、天文、物理等領域有廣泛應用。帕松分布的期望值為如相繼兩個事件出現的間隔時間服從參數為 的指數分布，則在單位時間間隔內事件出現的次數服從參數為 的泊松分布即單位時間內該事件出現k次的概率為：,2,1,0,!)(kkekXPk反之亦然。指數分布與帕松分布的關系：(1)(1)指兩個顧客到達商店的平均間隔時間是指兩個顧客到達商店的平均間隔時間是1010個單位時間個單位時間.即平均即平均1010個單位時間到達個單位時間到達1 1個顧客個顧客.(2)(2)指一個單位時間內平均到達指一個單位時間內平均到達0.10.1個顧客個顧客例例 (1)(1)顧客到達某商店的間隔時間服從參數為顧

13、客到達某商店的間隔時間服從參數為0.10.1的指數分布的指數分布 (2)(2)該商店在單位時間內到達的顧客數服從參數為該商店在單位時間內到達的顧客數服從參數為0.10.1的帕松分布的帕松分布 返回例例2 2敵坦克分隊對我方陣地實施突襲，其到達規律服從泊松分布，平均每分鐘到達輛（1）模擬敵坦克在分鐘內到達目標區的數量，以及在第、分鐘內各到達幾輛坦克（2）模擬在3分鐘內每輛敵坦克的到達時刻。（1）用poissrnd(4)進行模擬。To Matlab（poiss）（2）坦克到達的間隔時間應服從參數為4的負指數分布，用exprnd（1/4）模擬。To Matlab（time）連續系統模擬實例:追逐問題

14、追逐問題 狀態隨時間連續變化的系統稱為連續系統連續系統。對連續系統的計算機模擬只能是近似的，只要這種近似達到一定的精度，也就可以滿足要求。例例 追逐問題追逐問題:如圖,正方形ABCD的四個頂點各有一人.在某一時刻,四人同時出發以勻速v=1米/秒按順時針方向追逐下一人,如果他們始終保持對準目標,則最終按螺旋狀曲線于中心點O.試求出這種情況下每個人的行進軌跡.OBCDA1.建立平面直角坐標系:A(x1,y1),B(x2,y2),C(x3,y3),D(x4,y4).2.取時間間隔為t,計算每一點在各個時刻的坐標.設某點在 t 時刻的坐標為:),(iiyx則在tt時刻的坐標為:)sin,cos(tvy

15、tvxii其中 dxxii1cos dyyii1sin 2121)()(iiiiyyxxd3.取足夠小的,d時結束算法.4.對每一個點，連接它在各時刻的位置,即得所求運動軌跡.求解過程求解過程:To Matlab(chase)返回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 onendd=20;while(d0.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;

16、y(i)=y(i)+v*dt*(y(i+1)-y(i)/d;plot(x(i),y(i),.),hold on end end計算程序計算程序:To Matlab(chase)返回離散系統模擬實例離散系統模擬實例:排隊問題排隊問題 排隊論排隊論主要研究隨機服務系統的工作過程。在排隊系統中，服務對象的到達時間和服務時間都是隨機的。排隊論通過對每個個別的隨機服務現象的統計研究，找出反映這些隨機現象平均特性的規律，從而為設計新的服務系統和改進現有服務系統的工作提供依據。對于排隊服務系統,顧客常常注意排隊的人是否太多,等候的時間是否長,而服務員則關心他空閑的時間是否太短.于是人們常用排隊的長度、等待的

17、時間及服務利用率等指標來衡量系統的性能.1 系統的假設：系統的假設：（1）顧客源是無窮的；顧客源是無窮的；（2）排隊的長度沒有限制；排隊的長度沒有限制；（3）到達系統的顧客按先后順序依次進入服務，到達系統的顧客按先后順序依次進入服務，即即“先到先服務先到先服務”。單服務員的排隊模型單服務員的排隊模型：在某商店有一個售貨員，顧客陸續來到，在某商店有一個售貨員，顧客陸續來到，售貨員逐個地接待顧客當到來的顧客較多時，一部分顧客便須排隊售貨員逐個地接待顧客當到來的顧客較多時，一部分顧客便須排隊等待，被接待后的顧客便離開商店設：等待，被接待后的顧客便離開商店設：1 1顧客到來間隔時間服從參數為顧客到來間

18、隔時間服從參數為0.10.1的指數分布的指數分布對顧客的服務時間服從，上的均勻分布對顧客的服務時間服從，上的均勻分布排隊按先到先服務規則，隊長無限制排隊按先到先服務規則，隊長無限制 假定一個工作日為假定一個工作日為8 8小時，時間以分鐘為單位。小時，時間以分鐘為單位。11模擬一個工作日內完成服務的個數及顧客平均等待時間模擬一個工作日內完成服務的個數及顧客平均等待時間t t22模擬模擬100100個工作日，求出平均每日完成服務的個數及每日顧客的平均個工作日，求出平均每日完成服務的個數及每日顧客的平均等待時間。等待時間。2 2 符號說明符號說明 w w：總等待時間；：總等待時間；c ci i：第：

19、第i i個顧客的到達時刻；個顧客的到達時刻；b bi i：第：第i i個顧客開始服務時刻；個顧客開始服務時刻；e ei i：第：第i i個顧客服務結束時刻個顧客服務結束時刻 x xi i：第第i-1i-1個顧客與第個顧客與第i i個顧客之間到達的間隔時間個顧客之間到達的間隔時間 y yi i：對第對第i i個顧客的服務時間個顧客的服務時間c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci=ci-1+xiei=bi+yibi=max(ci,ei-1)t3 模擬框圖模擬框圖初始化：令i=1，ei-1=0，w=0產生間隔時間隨機數xi參數為0.1的指數分布ci=xi,bi=xi 產生服務時

20、間隨機數yi4,15的均勻分布ei=bi+yi累計等待時間：w=w+bi-ci準備下一次服務：i=i+1產生間隔時間隨機數xi參數為0.1的指數分布ci=ci-1+xi 確定開始服務時間：bi=max(ci,ei-1)bi480?YNi=i-1,t=w/i輸出結果：完成服務個數：m=i 平均等待時間：t停止1模擬一日To Matlab(simu1)2模擬100日To Matlab(simu2)返回用蒙特卡洛法解非線性規劃問題用蒙特卡洛法解非線性規劃問題對于非線性規劃問題非線性規劃問題：min f(X)XnE s.t.0)(Xgi i=1,2,m jjjbxa j=1,2,n用蒙特卡洛法求解的基

21、本思想求解的基本思想是：在估計的區域 (x1,x2,xn)|xj,jjba,j=1,2,n內隨機取若干實驗點，然后從試驗點中找出可行點，再從可行點中選擇最小點基本假設基本假設 試驗點的第j個分量xj服從aj,bj內的均勻分布符號假設符號假設 P:試驗點總數；MAXP：最大試驗點總數；K:可行點總數；MAXK：最大可行點數;X*：迭代產生的最優點；Q：迭代產生的最小值 f(X*)，其初始值為計算機所能表示的最大數求解過程求解過程 先產生一個隨機數作為初始試驗點,以后則將上一個試驗點的第j個分量隨機產生，其它分量不變而產生一新的試驗點這樣，每產生一個新試驗點只需一個新的隨機數分量當KMAXK或PM

22、AXP時停止迭代框框 圖圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數xj=aj+R(bj-aj)j=1,2,nj=0j=j+1,p=p+1PMAXP?YNxj=aj+R(bj-aj)gi(X)0?i=1,2nYNjMAXK?YN輸出X,Q,停止YN 例例 max 21212221382xxxxxxz s.t 10321 xx 01x 02x 在Matlab軟件包中編程，共需三個文件:randlp.m,mylp.m,lpconst.m.主程序為randlp.m.%mylp.m%mylp.mfunction z=mylp(x)%目標函數z=2*x(1)2+x(2)2-x(1)*x(

23、2)-8*x(1)-3*x(2);%轉化為求最小值問題%lpconst.mlpconst.mfunction lpc=lpconst(x)%約束條件if 3*x(1)+x(2)-10=-0.5%約束條件的誤差為5.0 lpc=1;else lpc=0;end%randlp.m%randlp.m function sol,r1,r2=randlp(a,b,n)%隨機模擬解非線性規劃debug=1;a=0;%試驗點下界b=10;%試驗點上界n=1000;%試驗點個數r1=unifrnd(a,b,n,1);%n1階的a,b均勻分布隨機數矩陣r2=unifrnd(a,b,n,1);sol=r1(1)r

24、2(1);z0=inf;for i=1:n x1=r1(i);x2=r2(i);lpc=lpconst(x1 x2)；if lpc=1 z=mylp(x1 x2)；if zz0 z0=z；sol=x1 x2；end endendTo Matlab(randlp)返回實驗作業實驗作業2.某報童以每份 0.03 元的價格買進報紙,以 0.05 元的價格出售.根據長期統計,報紙每天的銷售量及百分率為 銷售量200210220230240250百分率0.100.200.400.150.100.05已知當天銷售不出去的報紙，將以每份 0.02 元的價格退還報社.試用模擬方法確定報童每天買進報紙數量,使報

25、童的平均總收入為最大?1、編一個福利彩票電腦選號的程序。4.某設備上安裝有四只型號規格完全相同的電子管，已知電子管壽命為1000-2000小時之間的均勻分布。當電子管損壞時有兩種維修方案，一是每次更換損壞的那一只；二是當其中一只損壞時四只同時更換。已知更換時間為換一只時需1小時，4只同時換為2小時。更換時機器因停止運轉每小時的損失為20元，又每只電子管價格10元，試用模擬方法決定哪一個方案經濟合理？5.導彈追蹤問題：導彈追蹤問題：設位于坐標原點的甲艦向位于x軸上點A(1,0)處的乙艦發射導彈，導彈頭始終對準乙艦.如果乙艦以最大的速度(是常數)沿平行于y軸的直線行駛，導彈的速度是5，模擬導彈運行

26、的軌跡.又乙艦行駛多遠時，導彈將它擊中？初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i20?E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20停止硬幣正面?YNNY1，2，34，56投擲硬幣的計算機模擬投擲硬幣的計算機模擬1、產生服從U（0，1）的隨機數R12、將區間0，1兩等分：若 ，則對應硬幣正面 若 ，則對應硬幣反面5.001R15.01R擲骰子的計算機模擬擲骰子的計算機模擬1、產生服從U（0，1）的隨機數R22、將區間0，1六等份：若 ，則對應骰子點數為1 若 ，則對應骰子點數為2 若 ，則對應骰子點數為3 若 ，則對應骰子點數為4 若 ，則對應骰子點數為5 若 ，則對應骰子點數為66102 R62612 R63622 R65642 R64632 R1652 R初始化:i=0,k1=0,k2=0,k3=0i=i+1R2=?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i20?E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20停止R1=0.5YNNYR25/6To Matlab(liti1)