系統(tǒng)仿真的MATLAB實(shí)現(xiàn)

1、第七章 系統(tǒng)仿真的MATLAB實(shí)現(xiàn)由于計(jì)算機(jī)技術(shù)的高速發(fā)展，我們可以借助計(jì)算機(jī)完成系統(tǒng)的數(shù)字仿真。綜前所述，數(shù)字仿真實(shí)質(zhì)上是根據(jù)被研究的真實(shí)系統(tǒng)的模型，利用計(jì)算機(jī)進(jìn)行實(shí)驗(yàn)研究的一種方法。仿真的主要過(guò)程是：建立模型、仿真運(yùn)行和分析研究仿真結(jié)果。仿真運(yùn)行就是借助一定的算法，獲得系統(tǒng)的有關(guān)信息。MATLAB是一種面向科學(xué)與工程計(jì)算的高級(jí)語(yǔ)言，它集科學(xué)計(jì)算、自動(dòng)控制、信號(hào)處理、神經(jīng)網(wǎng)絡(luò)和圖像處理等學(xué)科的處理功能于一體，具有極高的編程效率。MATLAB是一個(gè)高度集成的系統(tǒng)，MATLAB提供的Simulink是一個(gè)用來(lái)對(duì)動(dòng)態(tài)系統(tǒng)進(jìn)行建模、仿真和分析的軟件包，它支持線(xiàn)性和非線(xiàn)性系統(tǒng)，能夠在連續(xù)時(shí)間域、離散時(shí)

2、間域或者兩者的混合時(shí)間域里進(jìn)行建模，它同樣支持具有多種采樣速率的系統(tǒng)。在過(guò)去幾年里，Simulink已經(jīng)成為數(shù)學(xué)和工業(yè)應(yīng)用中對(duì)動(dòng)態(tài)系統(tǒng)進(jìn)行建模時(shí)使用得最為廣泛的軟件包。 MATLAB仿真有兩種途徑：（1）MATLAB可以在SIMULINK窗口上進(jìn)行面向系統(tǒng)結(jié)構(gòu)方框圖的系統(tǒng)仿真；（2）用戶(hù)可以在MATLAB的COMMAND窗口下，用運(yùn)行m文件，調(diào)用指令和各種用于系統(tǒng)仿真的函數(shù)，進(jìn)行系統(tǒng)仿真。這兩種方式可解決任意復(fù)雜系統(tǒng)的動(dòng)態(tài)仿真問(wèn)題，前者編輯靈活，而后者直觀性強(qiáng)，實(shí)現(xiàn)可視化編輯。 下面介紹在MATLAB上實(shí)現(xiàn)幾類(lèi)基本仿真。7.1 計(jì)算機(jī)仿真的步驟在學(xué)習(xí)計(jì)算機(jī)仿真以前，讓我們先總結(jié)一下計(jì)算機(jī)仿真的

3、步驟。計(jì)算機(jī)仿真，概括地說(shuō)是一個(gè)“建模實(shí)驗(yàn)分析”的過(guò)程，即仿真不單純是對(duì)模型的實(shí)驗(yàn)，還包括從建模到實(shí)驗(yàn)再到分析的全過(guò)程。因此進(jìn)行一次完整的計(jì)算機(jī)仿真應(yīng)包括以下步驟：（1）列舉并列項(xiàng)目每一項(xiàng)研究都應(yīng)從說(shuō)明問(wèn)題開(kāi)始，問(wèn)題由決策者提供或由熟悉問(wèn)題的分析者 提供。（2）設(shè)置目標(biāo)及完整的項(xiàng)目計(jì)劃目標(biāo)表示仿真要回答的問(wèn)題、系統(tǒng)方案的說(shuō)明。項(xiàng)目計(jì)劃包括人數(shù)、研究費(fèi)用以及每一階段工作所需時(shí)間。（3）建立模型和收集數(shù)據(jù)模型和實(shí)際系統(tǒng)沒(méi)有必要一一對(duì)應(yīng)，模型只需描述實(shí)際系統(tǒng)的本質(zhì)或者描述系統(tǒng)中所研究部分的本質(zhì)。因此，最好從簡(jiǎn)單的模型開(kāi)始，然后進(jìn)一步建立更復(fù)雜的模型。 （4）編制程序和驗(yàn)證利用數(shù)學(xué)公式、邏輯公式和算法

4、等來(lái)表示實(shí)際系統(tǒng)的內(nèi)部狀態(tài)和輸入/輸出的關(guān)系。建模者必須決定是采用通用語(yǔ)言如MATLAB、FORTRAN、C還是專(zhuān)用仿真語(yǔ)言來(lái)編制程序。在本教材中，我們選擇的是MATLAB和其動(dòng)態(tài)仿真工具Simulink。（5）確認(rèn)確認(rèn)指確定模型是否精確地代表實(shí)際系統(tǒng)。它不是一次完成，而是比較模型和實(shí)際系統(tǒng)特性的差異，不斷對(duì)模型進(jìn)行校正的迭代過(guò)程。（6）實(shí)驗(yàn)設(shè)計(jì) 確定仿真的方案、初始化周期的長(zhǎng)度、仿真運(yùn)行的長(zhǎng)度以及每次運(yùn)行的重復(fù) 次數(shù)。 （7）生產(chǎn)性運(yùn)行和分析 通常用于估計(jì)被仿真系統(tǒng)設(shè)計(jì)的性能量度。利用理論定性分析、經(jīng)驗(yàn)定性分析或系統(tǒng)歷史數(shù)據(jù)定量分析來(lái)檢驗(yàn)?zāi)Ｐ偷恼_性，利用靈敏度分析等手段來(lái)檢驗(yàn)?zāi)Ｐ偷姆€(wěn)定性。

5、 （8）文件清單和報(bào)表結(jié)果 （9）實(shí)現(xiàn)圖7.1是計(jì)算機(jī)仿真的程序圖。 圖7.1 計(jì)算機(jī)仿真程序流圖7.2 基于數(shù)值積分法的連續(xù)系統(tǒng)仿真MATLAB的工具箱提供了各種數(shù)值積分方法函數(shù)，這些函數(shù)是ODE23、ODE45、ODE113和ODE15s。這些函數(shù)均是m文件，還有一個(gè)函數(shù)是ode1.C，是直接用C語(yǔ)言編寫(xiě)的。 函數(shù)ode23( )是用Runge-Kutta法求解微分方程。它是一種采用三階積分算法、二階誤差估計(jì)、變積分步長(zhǎng)的低階積分算法，調(diào)用格式為T(mén), Y = ode23 ( 'F', TSPAN, YO, OPTIONS )其中，F(xiàn)為系統(tǒng)模型文件名，模型為y' =

6、f( t, y )形式；TSPAN = To TFINAL 為積分計(jì)算時(shí)間，初值為T(mén)o，終值為T(mén)FINAL；YO為系統(tǒng)輸出初始值；OPTIONS選項(xiàng)積分計(jì)算相對(duì)允差 'RelTol' 和絕對(duì)允差 'AbsTol'，當(dāng)缺省時(shí)， Reltol1e-3, AbsTol1e-6T為計(jì)算點(diǎn)時(shí)間向量，Y為微分方程的解。函數(shù)ode45( )也是用Runge-Kutta法求解微分方程，它是變步長(zhǎng)的一種中等階次積分算法，調(diào)用格式為T(mén), Y = ode45 ( 'F' , TSPAN, YO, OPTIONS )各項(xiàng)含義同上。函數(shù)ode113( )是變階的Adams

7、-Bashforth-Moulton，用變階方法解微分方程，采用多步法，調(diào)用格式為T(mén), Y = ode113 ( 'F', TSPAN, YO, OPTIONS )各項(xiàng)含義同上。函數(shù)odel5s( )采用改進(jìn)的Gear法解微分方程，調(diào)用格式為 T, Y = odel5s ( 'F', TSPAN, YO, OPTIONS )各項(xiàng)含義同上。MATLAB還提供了解微分方程函數(shù)ode23s。函數(shù)ode1.C是用Euler法求解微分方程。在SIMULINK中，系統(tǒng)仿真有變步距和固定步距兩種工作方式。在變步距仿真中，可選用ode45、ode23、ode113、ose15s

8、和ode23s。在固定步距仿真中，可選用材ode5、ode4、ode3、ode2和ode1。ode5是5階Runge-Kutta法，ode4是4階Runge-Kutta法，odel是Euler法?！纠?.1】求微分方程。用MATLAB編寫(xiě)程序：建立一個(gè)m函數(shù)文件dfun.mfuncion ydhfn ( t，x ) ysqrt ( 3*x )5 ; 在MATLAB COMMAND窗口下運(yùn)行 % MATLAB PROGRAM 7-1 t, x ode23 ( 'dfun'，0 l0，1) ; Plot ( t，x ) ; 基于數(shù)值積分法的連續(xù)系統(tǒng)的數(shù)字仿真對(duì)于一個(gè)n階微分方程表示

9、的連續(xù)系統(tǒng)，也可以用具有n個(gè)狀態(tài)變量的狀態(tài)空間表達(dá)式來(lái)描述：狀態(tài)方程實(shí)際上是n個(gè)一階微分方程組成的方程組。如果應(yīng)用四階Runge-Kutta法進(jìn)行仿真計(jì)算。則該式可以改寫(xiě)為式中，K1，K2，K3，K4均為n維列向量。系統(tǒng)輸出為ykCXkDU( tk )類(lèi)似地，可將各種數(shù)值積分方法用于連續(xù)系統(tǒng)的仿真中?！纠?.2】用MATLAB編寫(xiě)圖7.2所示液壓控制系統(tǒng)的仿真程序。圖7.2 系統(tǒng)方框圖用MATLAB編寫(xiě)仿真程序，采用四階Runge-Kutta法：% MATLAB PROGRAM 7-2% * Input system data * %調(diào)入數(shù)據(jù)文件hynat;% Input system fun

10、ction; %調(diào)入系統(tǒng)模型ypfun1 = ' valve 'ypfun = ' hysys '% Initialization %初始化yref = 5000; % Referent value of system outputx0 = 0 0; % Initial value of servo valvey0 = 0 0 0; % Initial value of hydraulic cylinderu0 = 0; t0 = 0; % Start time of simulation tfinal = 1; % End time of simulation

11、tsamp = 0.01; % Sample periodh = 0.001; % Simulation step sizeY1imit = 512;Ilimit = 40;max_epoch = fix ( tfinal / h )1;t = t0;u1 = u0;x = x0'y = y0;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;

12、uout ( i ) = 1;yd ( i ) = Kf * y ( l );yout ( i, : ) = y'% The main loop for i = 1 : max_epoch% Compute output of valvesvl = feval ( ypfun1, t, u, x, av, bv );sv2 = fcval (ypfun1, t+h/2, u, x+h*sv1/2, av, hv ); %模塊valve仿真sv3 = feval (ypfun1, t+h/2, u, x+h*sv2/2, av, hv );sv4 = feval (ypfun1, t+h

13、, u, x+h*sv3, av, hv);x = x+h*(svl+2*sv2+2*sv3+sv4)/6;vo = cv*x; % Compute output of cylinderxl = feval (ypfun, t, vo, y, a, b);s2 = feval (ypfun, t+h/2, vo, y+h*s1/2, a, h); %模塊hysya 仿真s3 = feval (ypfun, t+h/2, vo, y+h*s2/2, a, b);s4 = feval (ypfun, t+h u, vo, y+h*s3, a, b);y = y+h* (sl+2*s2+2*s3+s

14、4)/6;i=i+lt=t+h;tout(i) =t;uout(i) =u;yd(i)=Kf*y(l);yout(i, :) = y'% Discrete control process %離散控制量計(jì)算if abs (round (t/tsamp)-t/tsamp)<le9 ye = yrefy(l)*Kf;ul=Kd*ye;% Saturation blockif ul>Ylimitxl=ylimit ;else if ul<ylimitxl =Ylimit;else xl = ul;endend% D/A conversionx2 = Kda*xl;% Ampl

15、ifier gain x3=Ka*x2;% V/I conversion u4=Kvi*x3;% limit of currentif u4>Ilimitx4=Ilimit;else if u4<Ilimitx4=Ilimit;else x4 =u4;end endu = Kq*x4;end % for discrete section end % for main loop% save data to file %存儲(chǔ)仿真數(shù)據(jù)hout = tout uout yout;save hout.dat hout -ascii;plot (tout, yd, ' y');

16、gird;m 函數(shù)文件 valve .mfunction xd = valve (t, u, x, av, bv) xd = av*x+hv*u;m 函數(shù)文件hysys. Mfunction yd=hysys(t, u, y, a, b) yd = a* y+h*u 上例的仿真結(jié)果如圖7.3所示。圖7.3 例7.2系統(tǒng)仿真結(jié)果7.3 基于離散相似法的連續(xù)系統(tǒng)仿真 前面討論的連續(xù)系統(tǒng)仿真采用的方法是數(shù)值積分法。對(duì)于連續(xù)系統(tǒng)仿真還有一種方法是離散相似法。所謂離散相似法是首先將連續(xù)系統(tǒng)模型離散化，得到等價(jià)的或相似的離散化的模型，然后對(duì)相似的離散模型進(jìn)行仿真計(jì)算。根據(jù)這一原理，首先應(yīng)將連續(xù)時(shí)間系統(tǒng)模型

17、轉(zhuǎn)換為等價(jià)的離散時(shí)間系統(tǒng)模型。連續(xù)系統(tǒng)離散化處理是通過(guò)采樣保持器來(lái)實(shí)現(xiàn)的。從一個(gè)連續(xù)時(shí)間系統(tǒng)轉(zhuǎn)換為相應(yīng)的離散時(shí)間模型，兩者的等價(jià)特性取決于采樣保持器。 MATLAB提供了將連續(xù)系統(tǒng)模型轉(zhuǎn)換為離散時(shí)間系統(tǒng)模型的函數(shù)C2D，調(diào)用格式為sysd = c2d (sys, Ts)其中，sys為線(xiàn)性連續(xù)時(shí)間系統(tǒng)；Ts為采樣時(shí)間；sysd為等價(jià)的離散時(shí)間系統(tǒng)。當(dāng)采樣保持器為零階保持器時(shí)； sysd = c2d (sys, Ts, method)其中，method為離散化方法，MATLAB提供以下幾種離散化方法，可以選用： 'zoh '為零階保持器 'foh'為一階保持器 &#

18、39;tustion'為雙線(xiàn)性變換法， 'prewarp'為改進(jìn)的雙線(xiàn)性變換法 'matched'使連續(xù)和離散系統(tǒng)具有匹配的DC增益【例7.3】連續(xù)系統(tǒng)傳遞函數(shù)，采樣一階采樣保持器，采樣周期為，求其離散化系統(tǒng)模型，并比較離散前后系統(tǒng)階躍響應(yīng)。用MATLAB編寫(xiě)程序： % MATLAB PROGRAM 7-3sys = tf ( l -1 , 14 5 , 'td' , 0.35 );sysd = c2d ( sysc, 0.l, 'for' )step ( sysc, sysd );運(yùn)行結(jié)果：Transfer functi

19、on: Sampling time: 0.1圖7.4 例7.3離散前后系統(tǒng)階躍響應(yīng)比較7.4 離散事件系統(tǒng)仿真的MATLAB實(shí)現(xiàn)應(yīng)用實(shí)例MATLAB提供了兩個(gè)基本的函數(shù)用于產(chǎn)生隨機(jī)數(shù)。它們是rand和randn。1rand它的作用是產(chǎn)生(0, 1)間均勻分布的隨機(jī)數(shù)或向量。其調(diào)用形式有Y = rand ( n ) 返回一個(gè)n維的U(0,l)隨機(jī)數(shù)方陣，n必須是一個(gè)標(biāo)量Y = rand (m, n)或Y = rand (m n) 返回一個(gè)m×n維的U(0, l)隨機(jī)數(shù)矩陣，n必須是一個(gè)標(biāo)量 Y = rand (size(A) ) 返回和矩陣A相同維數(shù)的矩陣 rand 返回單個(gè)的隨機(jī)數(shù)

20、s = rand ( 'state' ) 返回均勻分布隨機(jī)數(shù)發(fā)生器的當(dāng)前狀態(tài)矢量，其維數(shù)是35。這種形式一般不常用。MATLAB 5以上的版本使用了多種隨機(jī)數(shù)發(fā)生器，它可以產(chǎn)生 2-53，253 內(nèi)的所有浮點(diǎn)數(shù)。它的序列周期為21492?！纠?.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函數(shù)可以起到聲明變量的作用a = 0 00 0>&

21、gt;rand (size(a)ans = 0.1996 0.5385 0.3031 0.9201>>(3-1)*rand(2, 3)+l %產(chǎn)生1,3區(qū)間上的均勻分布隨機(jī)數(shù)ans = 2.0506 1.0689 2.5374 1.6137 2.4307 1.11902randn其作用是產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)或隨機(jī)矢量，其調(diào)用形式和rand類(lèi)似?！纠?.5】>>randn(1,2)ans =-0.4326 -1.6656>>randn(1, 2) / 2+1 %產(chǎn)生N(1, 0.25)的正態(tài)分布。ans = l.0627 1.14383randperm其作

22、用是產(chǎn)生一個(gè)隨機(jī)置換，調(diào)用形式p = randperm(n) 產(chǎn)生1:n的一個(gè)隨機(jī)置換?！纠?.6】>>p = randperm(4)p = 2 3 4 1以rand和randn為基礎(chǔ)可以產(chǎn)生給定分布的隨機(jī)數(shù)，例如泊松分布和二項(xiàng)分布。M文件函數(shù)posong演示了如何產(chǎn)生泊松分布的隨機(jī)數(shù)。function p = posong(a, n, m)w = ones(n, m)*(-a);th = exp(w);r = ones(n, m)dv = r > th;p = zeros(n, m);while nnz(dv ) > 0r = r.*(rand(n, m).*dv);

23、dv = r > th;p = p+dv;end在這里我們對(duì)它的MATLAB的實(shí)現(xiàn)，簡(jiǎn)要地說(shuō)明一下。參數(shù)a指泊松分布的參數(shù)，n, m則指返回的隨機(jī)矢量的維數(shù)。這個(gè)程序中有幾個(gè)地方需要注意：（1）MATLAB支持矢量的邏輯運(yùn)算，運(yùn)算結(jié)果是一個(gè)元素為0或1的矢量。如“r > th”這條語(yǔ)句。（2）函數(shù)nnz的作用是統(tǒng)計(jì)一個(gè)矢量非0元素的個(gè)數(shù)。在這里的作用是判斷是不是r中所有的元素小于exp(-a)。（3）運(yùn)算符“.*” 表示數(shù)組元素對(duì)元素的乘法，如 a1, a2, ，an.*b1, b2, , bn = al*bl, a2*b2, , an*bn總之這個(gè)程序比較充分的利用了MATLAB

24、基于向量的特點(diǎn)，產(chǎn)生二項(xiàng)分布隨機(jī)數(shù)M文件函數(shù)，請(qǐng)大家自己嘗試。其實(shí)，像上面的這些工作可以由MATLAB完成，MATLAB的統(tǒng)計(jì)工具箱提供了一個(gè)各種分布隨機(jī)數(shù)的發(fā)生器函數(shù)。安裝MATLAB的統(tǒng)計(jì)工具箱，就可以使用這些函數(shù)。表7.1列出了這些函數(shù)。你可以通過(guò)help命令在MATLAB界面上來(lái)獲得這些函數(shù)的詳盡使用方法。下面作為示例來(lái)看看random的用法。4Random它的作用是產(chǎn)生一個(gè)由name參數(shù)指定的分布的隨機(jī)數(shù)。其調(diào)用形式為：R = random (name, A, M, N ), name指定分布的形式，A說(shuō)明該分布的參數(shù)，M，N規(guī)定產(chǎn)生的隨機(jī)數(shù)的矩陣維數(shù)，M為行數(shù)，N為列數(shù)，它們的缺省

25、值都是1（MATLAB所有函數(shù)涉及維數(shù)的參數(shù)，順序都是行在前，列在后）。由于name指定的分布可能不止由一個(gè)參數(shù)來(lái)表示，所以調(diào)用形式要做相應(yīng)的變化。R = random(name, A, B, M, N)，指定的分布含有兩個(gè)參數(shù)。R = random(name, A, B, C, M, N)，指定的分布含有三個(gè)參數(shù)。例如，泊松分布只含有一個(gè)參數(shù)，就要選用第一種形式來(lái)調(diào)用。>>random ( 'poiss', 4, 2, 2)ans =6 46 1表7.1 MATLAB統(tǒng)計(jì)工具箱里的隨機(jī)數(shù)發(fā)生器報(bào)童問(wèn)題仿真報(bào)童問(wèn)題是一個(gè)古典的概率統(tǒng)計(jì)分析問(wèn)題，雖然問(wèn)題本身并不復(fù)雜，

26、但作為一個(gè)演示實(shí)例，可以反映出離散事件系統(tǒng)計(jì)算機(jī)仿真的很多特征。 1報(bào)童問(wèn)題 一報(bào)童從報(bào)刊發(fā)行中心定報(bào)后零售，每賣(mài)一份報(bào)紙可賺錢(qián)a元，若定報(bào)后賣(mài)不出去，則可再退回發(fā)行處，此時(shí)每退一份報(bào)要賠錢(qián)b元。雖然每天賣(mài)出報(bào)的份數(shù)是隨機(jī)的，但報(bào)童可根據(jù)以往賣(mài)報(bào)情況的統(tǒng)計(jì)來(lái)獲得每天賣(mài)k份的概率p(k)，試求報(bào)童每天期望受益達(dá)到最大的定報(bào)量z'。 2數(shù)學(xué)摸型 設(shè)報(bào)童每天訂報(bào)z份，而報(bào)紙每天賣(mài)出y份，我們假設(shè)y的分布為考慮到報(bào)童每天的損失有如下兩種情形。（1）供過(guò)于求。因退貨造成的平均損失為：（2）供不應(yīng)求。因缺貨造成的平均損失為所以，每天的期望損失費(fèi)（也可以從總收益的角度來(lái)考慮）為現(xiàn)在我們的目標(biāo)是求出使

27、得每天期望損失最小的定報(bào)量，換言之，就是使報(bào)童的每天期望總收益達(dá)到最大。寫(xiě)成一個(gè)目標(biāo)函數(shù)的形式約束條件如z的取值范圍，要受到報(bào)童的資本多少的影響。只有在特殊的概率分布情況下，我們才可以推導(dǎo)出C(z)的解析形式，并通過(guò)求極值的方法來(lái)求解。但在實(shí)際的應(yīng)用中，這樣的思路往往是行不通的。但是可以通過(guò)枚舉所有可能的訂報(bào)量，求出對(duì)應(yīng)的平均損失，進(jìn)行比較求出滿(mǎn)足條件的z，這里搜索域通常是有限的。以上就是一個(gè)比較簡(jiǎn)單的計(jì)算機(jī)仿真方法。3報(bào)童問(wèn)題的計(jì)算機(jī)仿真對(duì)于給定每一訂報(bào)量Z值，利用離散隨機(jī)變量采樣算法產(chǎn)生給定分布的隨機(jī)數(shù)R，用來(lái)表示報(bào)童當(dāng)天賣(mài)出的報(bào)紙數(shù)，從而可以計(jì)算出一天的損失以及一個(gè)階段的平均損失。這里比

28、較關(guān)鍵的一點(diǎn)是如何產(chǎn)生服從給定分布的隨機(jī)變量，這個(gè)內(nèi)容在本教材中不再詳盡介紹。而且在實(shí)際的應(yīng)用中，分布并非總是給定了的，需要我們收集數(shù)據(jù)，并從中辨識(shí)分布，進(jìn)行參數(shù)估計(jì)。圖7.5是根據(jù)上述思路設(shè)計(jì)出的仿真流程框圖。其中各變量含義如下：Tm一輪試驗(yàn)的預(yù)定模擬天數(shù) T一輪實(shí)驗(yàn)的仿真天數(shù)累計(jì)值Z訂報(bào)量 Z最優(yōu)訂報(bào)量G定報(bào)量Z之上界 S1損失值之累計(jì)值S最小損失值這里，a和b是這個(gè)問(wèn)題的兩個(gè)參數(shù)。圖7.5 報(bào)童問(wèn)題計(jì)算機(jī)仿真流程框圖4計(jì)算機(jī)仿真根據(jù)圖7.5所示框圖，我們不難寫(xiě)出報(bào)童問(wèn)題的仿真程序，其MATLAB源碼如下。function superz, supers = baotong ( tm, g,

29、 a, b)z = 1;supers =1000;while z<g% r=poisondis (5, 1, tm);r = round (2*randn (1, tm)+5);%產(chǎn)生均勻分布隨機(jī)數(shù)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;if supers>=aver_ssupers = aver_s;superz = z;end;z = z +1;end;上面的代碼，是在隨機(jī)變量均勻分布的假設(shè)下編寫(xiě)的，對(duì)于其他的分布，大家只要用相應(yīng)的隨機(jī)數(shù)發(fā)生器進(jìn)行替換即可

30、。 5仿真結(jié)果和輸出數(shù)據(jù)分析 在分布為均勻分布，參數(shù)a=0.2，b=0.4的條件下，在MATLAB命令窗口運(yùn)行仿真程序就可以得出仿真的結(jié)果。例如 >> z, s = baotong (5, 10, 0.2, 0.4) z = 4s = 0.2400此外，我們還可以改變參數(shù)a、b的值，觀察相應(yīng)最優(yōu)值的變化，用MATLAB可以很方便的畫(huà)出C-a，C-b的變化曲線(xiàn)。6報(bào)童問(wèn)題模擬系統(tǒng)的推廣和應(yīng)用報(bào)童仿真問(wèn)題的改進(jìn)推廣和應(yīng)用可以有以下幾個(gè)方面：（1）求每天的賣(mài)出報(bào)數(shù)服從任意分布的情況下，使報(bào)童收益最大意義下的最優(yōu)訂報(bào)量Z；（2）對(duì)報(bào)紙總發(fā)行量進(jìn)行測(cè)算；（3）適當(dāng)修改仿真系統(tǒng)，可將其用于其它

31、系統(tǒng)，例如企業(yè)的訂貨和庫(kù)存策略研究。7.5 SIMULINK動(dòng)態(tài)仿真 SIMULINK是MATLAB重要軟件包，用于對(duì)動(dòng)態(tài)系統(tǒng)仿真，它適用于連續(xù)系統(tǒng)和離散系統(tǒng)，也適用線(xiàn)性系統(tǒng)和非線(xiàn)性系統(tǒng)。它采用系統(tǒng)模塊直觀地描述系統(tǒng)典型環(huán)節(jié)，因此可十分方便地建立系統(tǒng)模型而不需要花較多時(shí)間編程。正由于這些特點(diǎn)，SIMULINK廣泛流行，被認(rèn)為是最受歡迎的仿真軟件。本節(jié)對(duì)SIMULINK和它的MATLAB實(shí)現(xiàn)，作一個(gè)簡(jiǎn)單介紹，為大家今后運(yùn)用SIMULINK打下基礎(chǔ)。SIMULINK實(shí)際上是面向結(jié)構(gòu)的系統(tǒng)仿真軟件。利用它進(jìn)行系統(tǒng)仿真的步驟 如下：（1） 啟動(dòng)SIMULINK，進(jìn)入SIMULINK窗口；（2） 在SI

32、MULINK窗口下，借助SIMULINK模塊庫(kù)，創(chuàng)建系統(tǒng)框圖模型并調(diào)整模塊參數(shù)；（3） 設(shè)置仿真參數(shù)后，啟動(dòng)仿真；（4） 輸出仿真結(jié)果。用戶(hù)首先進(jìn)入MATLAB COMMAND窗口，鍵入Simulink，立即彈出Simulink模塊庫(kù)窗口，如圖7.6所示。Library：Simulink窗口是Simulink模塊庫(kù)窗口，共有7個(gè)標(biāo)準(zhǔn)模塊子庫(kù)。一個(gè)MATLAB演示程序圖標(biāo)。這七個(gè)標(biāo)準(zhǔn)模塊子庫(kù)是Source（信號(hào)源）、Sink（顯示輸出）、Discrete（離散）、Linear（線(xiàn)性）、Nonlinear（非線(xiàn)性）、Connections（連接）、Blocksets & Toolboxes

33、。每個(gè)模塊子庫(kù)又包含許多常用的模塊，供用戶(hù)建模時(shí)使用，SIMULINK提供的庫(kù)模塊清單見(jiàn)附錄。圖7.6 Simulink模塊庫(kù)窗口 系統(tǒng)框圖模型建立的過(guò)程如下： 1建立模型窗口 建立新的模型窗口常有四種方法： （1）在MATLAB COMMAND窗口下，鍵入Simulink，彈出Simulink模塊庫(kù)窗口同時(shí)，也彈出一個(gè)Untitled窗口，該窗口為未取名的模型窗口。用戶(hù)可在該窗口下建立新的系統(tǒng)框圖模型。（2）在Simulink窗口下，用鼠標(biāo)選取菜單 File 中 New 子菜單的 Mode 后，會(huì)彈出一個(gè)untitled窗口，如圖7.7所示。該窗口供用戶(hù)建立系統(tǒng)框圖模型。圖7.7 模型窗口

34、（3）若模型文件已存在，Simulink窗口下，選擇菜單 File 中 Open 命令，輸入文件名，即打開(kāi)一個(gè)已存在的方框圖模型。（4）若模型文件已存在，在MATLAB COMMAND窗口下，選擇 File中Run Script命令，輸入文件名，也可打開(kāi)一個(gè)已存在的方框圖模型。（5）若模型文件已存在，在MATLAB COMMAND窗口下，鍵入文件名，也可打開(kāi)一個(gè)已存在的模型。2選取模塊或模塊組 在建立框圖模型過(guò)程中，需進(jìn)行如拷貝、刪除模塊等操作。必須首先選擇模塊或模塊組，具體操作如下： （1）在模型或模塊庫(kù)的窗口內(nèi)，找出所需模塊圖標(biāo)，用鼠標(biāo)左鍵單擊。圖標(biāo)四角出現(xiàn)黑圓點(diǎn)，表示該模塊已被選中。 （

35、2）在模型或模塊庫(kù)窗口內(nèi)，用鼠標(biāo)左鍵在窗口矩形邊界兩個(gè)對(duì)角單擊一下，即生成一個(gè)邊界框?qū)⑺鑾讉€(gè)框塊圖標(biāo)包圍，松開(kāi)鼠標(biāo)，則邊界框內(nèi)模型和連接線(xiàn)出現(xiàn)黑圓點(diǎn)，表示這些模型（包括在連接線(xiàn)）均被選中。 用同樣方法可以選取一個(gè)系統(tǒng)框圖模型的全部模塊。 3模塊拷貝及刪除用戶(hù)在建立自己模型的時(shí)候，常常需要從Simulink模塊庫(kù)、其他模塊庫(kù)或其他模型窗口中復(fù)制所需的模塊并移動(dòng)至自己的模型窗口內(nèi)。有兩種操作方法：鼠標(biāo)拖動(dòng)方法和菜單命令法。 鼠標(biāo)拖動(dòng)法如下：（1） 打開(kāi)模塊庫(kù)窗口或模型窗口。（2） 將鼠標(biāo)移至要拷貝的模塊圖標(biāo)上，按下鼠標(biāo)左鍵并保持。（3） 移動(dòng)鼠標(biāo)將模塊圖標(biāo)拖至目標(biāo)模型窗口一定位置。（4） 松開(kāi)鼠

36、標(biāo)左鍵，模塊圖標(biāo)保留在目標(biāo)模型窗口內(nèi)，模塊拷貝完成。 菜單命令復(fù)制方法如下：（1） 打開(kāi)模塊庫(kù)窗口或模型窗口，按上述方法選取要復(fù)制的模塊。（2） 從Edit菜單中選取Copy命令。（3） 用鼠標(biāo)單擊目標(biāo)模型適當(dāng)位置。（4）從Edit菜單中選取Paste命令，要拷貝的模塊出現(xiàn)在目標(biāo)模型窗口內(nèi)，模塊復(fù)制完成。 模塊的刪除有兩種方法：（1）選取要?jiǎng)h除的模型，從Edit菜單中選取Clear或Cut命令，用Cut命令刪除的模塊允許使用Paste粘貼在另一個(gè)地方。（2）選取要?jiǎng)h除的模塊，并按Del鍵。 SIMULINK允許模塊更名，圖標(biāo)大小改變、模塊圖標(biāo)移動(dòng)、模塊圖標(biāo)旋轉(zhuǎn)等操作。 4模塊參數(shù)設(shè)置 用鼠標(biāo)雙

37、擊待設(shè)置參數(shù)的模塊圖標(biāo)。打開(kāi)模塊對(duì)話(huà)框，按對(duì)話(huà)框欄目中提供的信息，輸入或改變模塊參數(shù)。 按Close模塊參數(shù)設(shè)置或修改完成。 各模塊對(duì)話(huà)框需設(shè)置參數(shù)不同。大家打開(kāi)對(duì)話(huà)框后，一目了然。 5模塊連接線(xiàn) 模塊之間的連接線(xiàn)是信號(hào)線(xiàn)，每根連接線(xiàn)都表示標(biāo)量或向量信號(hào)的傳輸，連接線(xiàn)的箭頭表示信號(hào)流向。連接線(xiàn)把一個(gè)模塊的輸出端口和另一個(gè)模塊的輸入端口連接起來(lái)，也可以利用分支線(xiàn)把一個(gè)模塊的輸出端口和幾個(gè)模塊的輸入端口連接 起來(lái)。 模塊間連接線(xiàn)生成方法如下：（1） 將鼠標(biāo)置于某模塊的輸出端口上，立即呈現(xiàn)出一個(gè)十字形光標(biāo)。（2） 拖動(dòng)十字光標(biāo)至另一模塊的輸入端口，在兩模塊之間生成一條帶箭頭的連線(xiàn)，模塊連線(xiàn)完成。 分

38、支線(xiàn)的生成方法如下： 用鼠標(biāo)鍵單擊被選取好信號(hào)線(xiàn)，按鼠標(biāo)右鍵對(duì)準(zhǔn)該信號(hào)線(xiàn)后顯示十字形光標(biāo)，拖動(dòng)它會(huì)生成信號(hào)線(xiàn)的一條分支線(xiàn)。按此法，可生成多條分支線(xiàn)。 6模型文件取名及保存 一旦把模型窗口上各模型連接起來(lái)，一個(gè)系統(tǒng)方框圖模型建立工作就已完成。 選擇模型窗口File菜單中Save as命令，彈出對(duì)話(huà)框，填入模型文件名。如hexam，按保存(S)，模型以文件名hexam.mdl保存在用戶(hù)的子目錄內(nèi)。按上面方法建立一個(gè)簡(jiǎn)單系統(tǒng)框圖模型，如圖7.8所示，文件名為hexam。圖7.8 簡(jiǎn)單系統(tǒng)框圖模型 系統(tǒng)仿真運(yùn)行常有兩種方法進(jìn)行，一種是利用SIMULINK菜單命令，二是在MABLAB COMMAND窗口

39、下輸入有關(guān)命令。1SIMULINK模型窗口下的仿真運(yùn)行 在SIMULINK模型窗口下進(jìn)行仿真操作簡(jiǎn)單、直觀，不必記憶命令的語(yǔ)法規(guī)則，人機(jī)交互方式選擇或修改仿真參數(shù)，模型參數(shù)等。具體操作如下：（1） 打開(kāi)系統(tǒng)模型窗口。（2）從菜單Simulation中選取Parameters，彈出Simulation Parameters對(duì)話(huà)框，如圖7.9所示。圖7.9 仿真參數(shù)對(duì)話(huà)框Simulation Parameters對(duì)話(huà)框用三頁(yè)管理仿真參數(shù)。 Solver頁(yè)，設(shè)置或修改參數(shù)有：仿真時(shí)間：開(kāi)始時(shí)間，終止時(shí)間仿真方法類(lèi)型： 變步距（積分算法一）：discrete ode45, ode23, odel13,

40、 odel5s, ode23s 定步距（積分算法二）：ode5，ode4, ode3, ode2, ode1, discrete 初始步距 最大步距 Workspace I/O頁(yè)設(shè)置或修改參數(shù)有： 從工作空間調(diào)用輸入 把數(shù)據(jù)存入工作空間（時(shí)間、狀態(tài)、輸出） 狀態(tài)值（初始值調(diào)入，初始值存放） 存儲(chǔ)選擇 Diagnostics頁(yè)設(shè)置或選擇項(xiàng)目：仿真過(guò)程警告信息（3）應(yīng)用仿真參數(shù)在通過(guò)Simulation Parameters對(duì)話(huà)框設(shè)置或修改好仿真參數(shù)后，準(zhǔn)備應(yīng)用到用戶(hù)模型仿真中，按對(duì)話(huà)中Apply鍵，再按對(duì)話(huà)框中Close鍵，關(guān)閉對(duì)話(huà)框。（4）仿真開(kāi)始和停止從Simulation菜單中選取Star

41、t命令。仿真立即開(kāi)始運(yùn)行，這時(shí)Start變?yōu)镾top。若要停止仿真，從Simulation菜單中選擇Stop。仿真運(yùn)行立即停止。若要使仿真運(yùn)行暫停，可從Simulation菜單中選擇Pause，這時(shí)Pause變?yōu)镃ontinue。若要使仿真繼續(xù)運(yùn)行，選擇Continue。2MATLAB的COMMAND窗口下的仿真運(yùn)行 在MATLAB COMMAND窗口下運(yùn)行SIMULINK模型仿真也可進(jìn)行，方法是在COMMAND窗口下調(diào)用函數(shù)SIM。調(diào)用格式為t, x, y = sim ( 'model', timespan, option, ut, p1, pn)其中，t為返回的仿真時(shí)間向量

42、；y為返回的輸入矩陣，從輸出模塊輸出的數(shù)據(jù)；model為框圖模型文件名或MEX-文件或M-文件S函數(shù)；timespan = tStart，tFinal；Options為仿真參數(shù)選項(xiàng)，由函數(shù)Simset( )產(chǎn)生；ut為加在第一級(jí)輸入模塊的外部輸入，如一個(gè)函數(shù)u(t)等；p1, pn為S-函數(shù)參數(shù)，如Sys = Sfun(t，x，u，flay，p1,p2, pn)僅用于M-文件和MEX-文件系統(tǒng)。例如，圖7.8所示系統(tǒng)，在MATLAB COMMAND窗口下仿真t, y = Sim ( 'hexm' );plot ( t, y )比較在SIMILINK模型窗口小的仿真結(jié)果，可見(jiàn)兩種

43、方法仿真結(jié)果相同。由于在MATLAB的COMMAND窗口下運(yùn)行SIMULINK模型不如直接在SIMULINK模型窗口下仿真方便、直觀，因此不推薦使用。利用Mum模塊和Scope模塊可顯示多個(gè)變量的運(yùn)動(dòng)軌跡，如圖7.10所示。圖7.10 顯示多個(gè)變量軌跡SIMULINK可通過(guò)下列三種方法之一繪制輸出軌跡曲線(xiàn)。1利用Scope模塊在仿真期間，Scope模塊動(dòng)態(tài)地顯示輸出軌跡的一個(gè)簡(jiǎn)單例子如圖7.10所示。Scope模塊還具有放大、縮小功能以便更清晰顯示感興趣的區(qū)域，Scope模塊還可將數(shù)據(jù)存儲(chǔ)在工作空間（Workspace）內(nèi)。利用XY圖形模型，可顯示兩個(gè)變量之間變化軌跡。2利用變量返回值利用Co

44、nnection模塊子庫(kù)的Out模塊，建立如圖7.11的模型。仿真運(yùn)行后，輸出量yout保存在工作空間內(nèi)，可供用戶(hù)調(diào)用，繪制輸出軌跡。圖7.11 Out模塊應(yīng)用也可從Simulation菜單調(diào)用Simulation Parameters對(duì)話(huà)框的Workspace I/O頁(yè)面，定義時(shí)間變量為tout，輸出變量為yout。仿真結(jié)束后，在MATLAB COMMAND窗口下鍵入plot ( tout, yout )觀察輸出軌跡。3利用To Workspace模塊建立如圖7.12所示的SIMULINK模型，模型中采用兩個(gè)To Workspace模塊，一個(gè)用于記錄時(shí)間變量t，一個(gè)用來(lái)記錄變量y。變量名t和y通過(guò)打開(kāi)To Workspace模塊對(duì)話(huà)