1、 自動化車床管理的數(shù)學模型摘 要本文研究的是自動化車床管理問題，該問題屬于離散型隨機事件的優(yōu)化模型，目的是使管理得到最優(yōu)化。首先我們借用maltlab中的lillietest函數(shù)對題目給出的100次刀具故障記錄的數(shù)據(jù)進行了數(shù)據(jù)處理和假設(shè)檢驗(見附錄一)，樣本數(shù)據(jù)與正態(tài)分布函數(shù)擬合得很好，從而接受了數(shù)據(jù)符合正態(tài)分布的假設(shè)，求得刀具壽命的概率密度函數(shù)的期望=600，標準差=196.6296，積分后求得刀具壽命的分布函數(shù)。對于問題(1)，我們建立起離散型隨機事件模型，以合格零件的平均損失期望作為目標函數(shù)，借用概率論與數(shù)理統(tǒng)計的方法列出方程組，并利用matlab以窮舉法(見附錄二)得出最優(yōu)檢查間隔為1

2、8個，最優(yōu)刀具更新間隔為368個，合格零件的平均損失期望為5.17元。對于問題(2)，我們建立單值目標函數(shù)最優(yōu)化模型，以平均合格零件的損失期望作為目標函數(shù)，并由題所給條件列出約束條件表達式。最后借用matlab編程求解(見附錄三)得出最優(yōu)檢查間隔為32個，最優(yōu)刀具更新間隔為320個，合格零件的平均損失期望為7.46元。對于問題(3)，我們采取的優(yōu)化策略是：進行一次檢查，如果是合格品則再進行一次檢查，后一次檢查為不合格品則換刀。在做定量分析時，我們將問題(2)中的目標函數(shù)和方程組在問題(3)的條件上做了相應改變，利用matlab用窮舉法求解(見附錄四)得出優(yōu)檢查間隔為32個，最優(yōu)刀具更新間隔為3

3、20個，合格零件的平均損失期望為6.40元。由結(jié)果可以看出問題(3)的檢查間隔和刀具更新間隔與問題(2)的結(jié)果相同，但合格零件的平均損失期望降低了1.06元。說明問題(3)的檢查方式較問題(2)更優(yōu)。關(guān)鍵詞：離散型隨機事件優(yōu)化模型 概率理論 擬合優(yōu)度 窮舉法1問題重述1.1問題背景我國是一個工業(yè)化大國，其中自動化車床生產(chǎn)在我國工業(yè)生產(chǎn)中扮演著舉足輕重的角色。因此能否對于自動化車床進行高效經(jīng)濟地管理直接關(guān)系到工業(yè)生產(chǎn)是否可以做到“低消耗，高產(chǎn)出”。對于自動化機床管理進行優(yōu)化符合我國“可持續(xù)發(fā)展”的戰(zhàn)略，同時對于環(huán)境資源的節(jié)約保護有著突出貢獻。對于一個工業(yè)化企業(yè)而言，在日趨激烈的市場競爭中，“成本

4、最小化，效率最大化”已經(jīng)成為其至關(guān)重要的生存之道。所以大到國家，小至企業(yè)，對“自動化車床管理”的研究都給予了高度重視。如今，數(shù)學模型分析已經(jīng)成為對該問題進行研究的主要途徑。1.2需要解決的問題一道工序用自動化車床連續(xù)加工某種零件，由于刀具損壞等原因該工序會出現(xiàn)故障，其中刀具損壞故障占95%, 其它故障僅占5%。工序出現(xiàn)故障是完全隨機的, 假定在生產(chǎn)任一零件時出現(xiàn)故障的機會均相同。工作人員通過檢查零件來確定工序是否出現(xiàn)故障?，F(xiàn)積累有100次刀具故障記錄，故障出現(xiàn)時該刀具完成的零件數(shù)如附錄表一。現(xiàn)計劃在刀具加工一定件數(shù)后定期更換新刀具。 已知參數(shù)： (1)故障時產(chǎn)出的零件損失費用 f=200元/件

5、； (2)進行檢查的費用 t=10元/次； (3)發(fā)現(xiàn)故障進行調(diào)節(jié)使恢復正常的平均費用 d=3000元/次(包括刀具費)； (4)未發(fā)現(xiàn)故障時更換一把新刀具的費用 k=1000元/次。 問題一：假定工序故障時產(chǎn)出的零件均為不合格品，正常時產(chǎn)出的零件均為合格品, 試對該工序設(shè)計效益最好的檢查間隔（生產(chǎn)多少零件檢查一次）和刀具更換策略。 問題二：如果該工序正常時產(chǎn)出的零件不全是合格品，有2%為不合格品；而工序故障時產(chǎn)出的零件有40%為合格品，60%為不合格品。工序正常而誤認有故障停機產(chǎn)生的損失費用為1500元/次。對該工序設(shè)計效益最好的檢查間隔和刀具更換策略.問題三:在問題三的情況下, 可否改進檢

6、查方式獲得更高的效益。 2模型的假設(shè)及符號說明2.1模型的假設(shè)假設(shè)1：假設(shè)在生產(chǎn)任一零件時出現(xiàn)故障的機會均相等。假設(shè)2：假設(shè)生產(chǎn)剛啟動時使用的刀具都是新的。假設(shè)3：假設(shè)生產(chǎn)任一零件時所需的時間相同。假設(shè)4：假設(shè)提供的刀具故障記錄數(shù)據(jù)是獨立同分布的.假設(shè)5： 假設(shè)無論刀具損壞故障還是其它故障, 發(fā)生故障并使恢復正常的平均費用均為3000元每次。假設(shè)6： 假設(shè)提供的刀具故障記錄數(shù)據(jù)是獨立同分布的。假設(shè)7： 假設(shè)發(fā)現(xiàn)故障和停機維修的時間可以忽略不計。假設(shè)8： 假設(shè)檢查時不停止生產(chǎn),在檢查出不合格零件時才停止再進行維修。假設(shè)9: 假設(shè)每次檢查只能檢查一個零件。假設(shè)10：假設(shè)5%的其它故障可以忽略不計。

7、2.2符號說明符號說明f故障時產(chǎn)生的零件損失費用200元/件t檢查的費用10元/次d發(fā)現(xiàn)故障進行調(diào)節(jié)使恢復正常的平均費用3000元/次k未發(fā)現(xiàn)故障時更換一把新刀具的費用1000元/次刀具平均壽命樣本方差Tc檢查零件的單位時間間隔T定期換刀的單位時間間隔T(c)以檢測時間間隔為時，系統(tǒng)工序合格零件的單位期望損失Tc*以經(jīng)濟損失最小為目標的最優(yōu)檢查的時間間隔T*以經(jīng)濟損失最小為目標的最優(yōu)的換刀間隔T(c)*在和的情況下，系統(tǒng)工序合格零件的最小單位期望損失f(x)系統(tǒng)的失效概率密度F(x)累積失效概率密度，即壽命分布函數(shù)3問題分析在自動化車床生產(chǎn)流程中，由于刀具損壞等原因會使工序出現(xiàn)故障，工序故障的

8、出現(xiàn)是完全隨機的。工作人員通過檢驗零件來確定是否出現(xiàn)故障，并且決定在刀具加工一定的零件后更換刀具。當發(fā)生故障時要及時維修，如果檢修周期太長，故障不能及時發(fā)現(xiàn)，會給生產(chǎn)帶來損失；檢查周期太短又會增加費用。在理論上我們首先將問題轉(zhuǎn)化為概率模型。通過分析題目所給的100次刀具故障記錄，我們通過繪圖分析假設(shè)刀具的壽命服從正態(tài)分布。再通過假設(shè)檢驗，我們決定接受這一假設(shè)。問題1中我們建立離散型隨機事件模型I。我們選擇一個周期T。目標函數(shù)要求目標函數(shù)取最小值的情況下求解檢查間隔和道具更換策略。U總分為兩種情況：故障發(fā)生在換刀之前與故障發(fā)生在換刀之后。我們分別求解這兩種情況。問題(2)中，可能會產(chǎn)生2中誤判。

9、誤檢：工序正常時，由于檢查到不合格品而停機，會產(chǎn)生一個費用。漏檢：工序不正常工作時由于有40%的合格品，會因為檢測到合格品而不換刀，導致不合格品增加。問題(3)中，由于問題(2)中由于工序正常時產(chǎn)出的零件不全是合格品，有2%為不合格品；而工序故障時產(chǎn)出的零件有40%為合格品，60%為不合格品。這樣會導致誤檢和漏檢，從而增加了損失。我們基于這種情況，建立當檢查到合格品時再檢查一零件，若任然是合格品則判定工序正常，否則，判定工序出現(xiàn)故障。這樣雖然會增加檢查成本。但是，也會大大減少誤檢。從而，可能使得損失減小。我們基于這樣一種思想對模型II進行改進。4數(shù)據(jù)處理及分析4.1刀具故障完成零件個數(shù)的數(shù)據(jù)統(tǒng)

10、計分析 我們將刀具壽命統(tǒng)計數(shù)據(jù)繪圖如下，再通過分布擬合檢驗法可以證明刀具故障數(shù)據(jù)近似服從正態(tài)分布。(編程見附錄二)4.1.2正態(tài)分布擬合度檢驗：函數(shù) lillietestH = lillietest(X) 表示對輸入向量X進行Lilliefors測試，顯著性水平為0.05。H = lillietest(X,alpha) 表示在水平alpha而非5%下施行Lilliefors測試，alpha在0.01和0.2之間。H,P,LSTAT,CV = lillietest(X,alpha) P為接受假設(shè)的概率值，P越接近于0，則越應該拒絕正態(tài)分布的原假設(shè)；LSTAT為測試統(tǒng)計量的值，CV為是否拒絕原假設(shè)

11、的臨界值。說明: H為測試結(jié)果，若H=0，則可以認為X是服從正態(tài)分布的；若H=1，則可以否定X服從正態(tài)分布。運用matlab計算可得(見附錄一) h = 0 p = 0.5000 l = 0.0421cv = 0.0890h=0表示接受正態(tài)分布的假設(shè)；p=0.5000表示服從正態(tài)分布的概率很大；統(tǒng)計量的值l=0.421小于接受假設(shè)的臨界值cv=0.0890，因而接受假設(shè)(測試水平為0.05)通過以上4個指標，可以得出結(jié)論：樣本中所給的數(shù)據(jù)與正態(tài)分布函數(shù)擬合得很好，我們接受這一假設(shè)。4.1.3概率密度函數(shù)求解.繪制正態(tài)分布函數(shù)圖象如下則累計失效概率密度函數(shù)(壽命分布函數(shù))為：5問題一的解答5.1

12、模型I的建立5.1.1確定目標函數(shù) 以合格零件單位期望為目標函數(shù)的目標函數(shù)：5.1.2問題一的求解I如果故障發(fā)生在換刀之后：換刀間隔T尚未出現(xiàn)故障時一次更新所消耗費用S1：(1)檢查費用：檢查費用等于檢查的次數(shù)乘以單次檢查所需的費用，即g1t一次換刀前未出現(xiàn)故障的過程的檢查次數(shù)，等于固定換刀間隔T除以檢查周期Tc所得的整數(shù)部分(2)換刀費用：k(3)不合格零件損失費用：0其中換刀前未出現(xiàn)故障的情況下總的損失費用為U1:換刀間隔T前尚未出現(xiàn)故障發(fā)生這種情況時的更新間隔均為T出現(xiàn)的次數(shù)等于刀具更新的總次數(shù)乘以以T為更新間隔情況下?lián)Q刀前仍未出現(xiàn)故障的概率，即,因此定期換刀前未出現(xiàn)故障的情況下的總損失

13、U1等于這種情況下的刀具更新次數(shù)乘以單位更新過程的損失費用S1即有：所以：如果故障發(fā)生在換刀之前更新過程所消耗的費用S2(1)發(fā)生故障時的維修換刀費用：d(2)故障維修前所有的損失費用：由于故障發(fā)生的隨機性，因此可以發(fā)生在T內(nèi)的任何位置因此這部分的損失費用等于對于周期T內(nèi)任意點x處發(fā)生故障所造成的損失與在x處可能發(fā)生故障的概率的乘積進行積分的結(jié)果，即其中表示在一個換刀周期T內(nèi)任意的x處發(fā)生故障概率。任意位置發(fā)生故障的損失費Wx=檢查費用+零件損失費a檢查費用等于檢查次數(shù)乘以單次檢查費用，即 (g2為發(fā)生故障x前的檢查次數(shù)，等于所得的整數(shù)部分) b零件損失費用等于從發(fā)生故障到維修檢查之間產(chǎn)生的不

14、合格零件數(shù)乘以單個零件的損失費用，即 (H為發(fā)生故障的檢查間隔內(nèi)產(chǎn)生的合格零件數(shù)，即發(fā)生故障前的所有合格零件數(shù)除以檢查間隔所得的余數(shù))所以換刀前已出現(xiàn)故障的情況下的損失總費用U2：定期換刀前出現(xiàn)故障的情況下的總損失U2等于這種情況下的刀具更新次數(shù)乘以單位更新過程的損失費用S2,即：換刀間隔T前就出現(xiàn)故障這時在故障發(fā)生后進行檢查并進行維修換刀，從而完成一個更新過程這種情況下總的發(fā)生次數(shù)等于總的更換次數(shù)乘以系統(tǒng)中發(fā)生這種情況的概率，即 ，F(xiàn)(T)是為以T為更新周期的情況下工序出現(xiàn)故障的概率，即為前面的數(shù)據(jù)處理中的累計失效概率密度函數(shù)，當t=T情況下F(T)的結(jié)果?？倱p失費用U總：合格零件總數(shù)： 系

15、統(tǒng)工序合格零件的單位期望損失：5.1.3綜上所述，我們對問題(1)建立離散型隨機事件模型目標函數(shù)：其中N表示：N次更新過程(每次換刀或維修換刀為一次更新過程)。5.1.4問題(1)的求解結(jié)果及分析運用matlab求解得每隔18個零件就該檢查一次，每隔368個零件換一次刀具。合格零件的平均損失期望結(jié)果分析：我們知道刀具的平均壽命為600.即每平均生產(chǎn)600個零件刀具就不能再使用。368 a=45936262454250958443374881550561245243498264074256570659368092665316448773460842811535938445275525137814

16、74388824538862659775859755649697515628954771609402960885610292837473677358638699634555570844166061062484120447654564339280246687539790581621724531512577496468499544645764558378765666763217715310851; u=sum(a)/100 %求期望值 b=60060060060060060060060060060060060060060060060060060060060060060060060060060060

17、0600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600600附錄二600600600600600600600600600600600600600600600600600600600600600600600600600600600600600; sum(a-b).*(a-b)/(100-1) o=sqrt(sum(a-b).*(a-b)/(100-1) %求標準差h,p,l,cv=lillie

18、test(a) %擬合度測試繪制正態(tài)分布函數(shù)圖象clearx=200:0.1:1000;f=exp(-(x-600).2./(2*(196.62922).*(1/(sqrt(2*pi)*196.6296); plot(x,f,+r)問題一的程序syms x Ta0=100;for T=368:402 for Tc=15:19 g1=ceil(T/Tc); g2=11;%ceil(x/Tc); f=exp(-(x-600).2./(2*(196.62922)*(1/(sqrt(2*pi)*196.6296); F=int(f,x,0,T);F=subs(F); h=rem(F,Tc);k1=in

19、t(10*(g2+1).*f./F,x,0,T)+int(200*(Tc-h).*f./F,x,0,T); f1=(10*g1+1000).*(1-F)+F*(3000+subs(k1,T); a=taylor(x.*f/F,x,0,6); f2=(1-F)*T+F*int(a,x,0,T);f2=subs(f2,T); Ta=f1/f2; Ta=simplify(Ta); if (TaTa0)附錄三 T1=T;Tc1=Tc;Ta0=Ta; else continue; end endendfprintf(最有檢查間隔為：n);T1fprintf(最有換刀間隔為：n);Tc1在7.0以上就可以

20、用taylor命令直接泰勒展開了，taylor(f,x,a,n)命令，使f函數(shù)泰勒展開，其中f為函數(shù)表達式，x為函數(shù)中的變量，在a點展開，n為展開的項數(shù)。要畫圖的話，用taylortool，默認函數(shù)xcosx，可以改。問題二的程序syms x Ta0=100;for T=301:321 for Tc=18:32 x0=ceil(T*rand(); f0=200;d=3000;t=10; g1=ceil(T/Tc); g2=11;%ceil(x/Tc); f=exp(-(x-600).2./(2*(196.62922)*(1/(sqrt(2*pi)*196.6296); %概率密度函數(shù) F=in

21、t(f,x,0,T); %分布函數(shù) F=subs(F); h=rem(F,Tc); x1=ceil(T-h); if x0x1 k1=int(40*(g2-1)+(1-0.4(g1-g2+1)*d+0.6*f0*(2*h+0.4*(1-0.4(g1-g2)/0.6*Tc)+(1-0.4(g1-g2+1)/0.6*t).*f./F,x,0,T)+int(taylor(200/5附錄四0*x.*f./F,x,0,11),x,0,T); f1=(40*g1+1000+T*f0/50)*(1-F)+F*subs(k1,T); a=taylor(0.98*x+0.4*(2*h+0.4*(1-0.4(g1

22、-g2)/0.6*Tc).*f/F,x,0,6); f2=(1-F)*T*0.98+F*int(a,x,0,T); f2=subs(f2,T); else k1=int(40*(g2-1)+(1-0.4(g1-g2+1)*d+0.6*f0*(2*h+0.4*(1-0.4(g1-g2)/0.6*Tc)+(1-0.4(g1-g2+1)/0.6*t).*f./F,x,0,T)+int(taylor(200/50*x.*f./F,x,0,11),x,0,T); f1=(40*g1+1000+T*f0/50)*(1-F)+F*subs(k1,T); a=taylor(0.4*T+0.58*x).*f/F

23、,x,0,6); f2=(1-F)*T*0.98+F*int(a,x,0,T); f2=subs(f2,T);end第三問程序說明：syms x Ta0=100;for T=301:321 for Tc=18:32 x0=ceil(T*rand(); f0=200;d=3000;t=10; g1=ceil(T/Tc); g2=11;%ceil(x/Tc); f=exp(-(x-600).2./(2*(196.62922)*(1/(sqrt(2*pi)*196.6296); %概率密度函數(shù) F=int(f,x,0,T); %分布函數(shù) F=subs(F); h=rem(F,Tc);附錄五 x1=ceil(T-h); if x0x1 k1=i