Monte-Carlo 方法在尺寸鏈方程組計算中的應用

蔡偉妹，仲其伐

(南京理工大學 機械工程學院，江蘇 南京 210094)

0 引言

引信作為一種微型化和精密化的軍工產品，技術和尺寸精度要求都比較高，因此尺寸鏈計算是引信產品設計中的一個重要環節，是引信設計的特點之一。引信產品的尺寸鏈計算主要為零件設計、產品裝配、制造工藝、發火性能以及自毀性能等的計算和校核提供數據，通過理論計算及時發現產品設計的缺陷，改善產品的可生產性、可裝配性以及產品使用的安全性與可靠性。

引信產品設計一般采用極大值極小值法計算尺寸鏈。文獻［1］從時間發生的概率大小和生產經濟性以及可生產性方面出發，提出用Monte-Carlo 方法解尺寸鏈的問題，將Monte-Carlo 方法應用到一般尺寸鏈的計算中;文獻［2］系統的介紹了求解由方程組組成的尺寸鏈問題的方法;文獻［3］通過Monte-Carlo 法求解多環尺寸鏈估算裝配的成功率問題;文獻［4］和文獻［5］則應用了Monte-Carlo 方法研究球轉子旋轉偏心情況，針對相關原則的同軸度誤差的處理問題，通過對計算條件的假設與抽樣檢驗，得出了精度比較高的計算仿真結果。文獻［7］解決了用常規方法對尺寸鏈組成環的尺寸分布服從非正態分布而無法求解的問題。

本文所要做的工作，是在文獻［1］～［6］所做工作的基礎上，將Monte-Carlo 方法具體應用到由方程組構成的尺寸鏈計算的實際問題中，進行分析求解，得出所求尺寸鏈封閉環的尺寸分布情況，進而求出封閉環的極值。

1 尺寸鏈計算方法

一般來說，進行尺寸鏈計算的主要目的是求解極值。所使用的計算方法一般有三類:極大值極小值法、概率統計法和Monte-Carlo 法，另外還有高等數學中的函數微分求極值法。

1.1 極大值極小值法

極大值極小值法是尺寸鏈計算中應用得較為廣泛的一種計算方法。應用該方法進行計算的基本原則［7］是:

1)計算封閉環最大尺寸時，所有增環都取最大極限尺寸，所有減環都取最小極限尺寸。

2)計算封閉環最小尺寸時，所有增環都取最小極限尺寸，所有減環都取最大極限尺寸。

3)計算封閉環基本尺寸時，所有組成環都取基本尺寸。

4)封閉環的最大最小尺寸的差值即為封閉環的公差。

極大值極小值法考慮了極端不利的情況，一般情況下極端條件出現的可能性很小，把這種極端情況考慮進去所計算出來的結果使得尺寸公差范圍小，生產的難度變大。一般來說，極大值極小值法適用于精度要求比較高，使用環境比較極端多變并要求產品裝配的100%互換的產品設計上。從引信設計角度和使用環境來說，引信的尺寸鏈計算一般采用極大值極小值法。

1.2 概率統計法

概率統計法從概率原理出發，考慮尺寸鏈中各組成環的實際分布情況，能客觀的反映加工精度和尺寸分布的本質。概率統計法的計算方法是，將各組成環的尺寸分布看作某一特定形態的分布，然后根據這一形態分布進行計算［7-8］。

一般情況下，當各組成環的分布規律按正態分布時，其封閉環尺寸也必然符合正態分布規律;當各組成環分布不按正態分布，在組成環數不太少，各組成環變化范圍大小相差不大的情況下，封閉環的尺寸仍趨于正態分布;當尺寸鏈組成環數較少而各組成環分布又偏離正態分布較大時，封閉環也將偏離正態分布。

一般而言，應用概率統計法計算引信零部件的尺寸鏈時，必須明確知道各組成環的尺寸分布形態。但是要明確知道各組成環的分布律是很困難的，需要長時間的統計資料積累，因此，概率統計法只適用于尺寸分布律已知的情況下的尺寸鏈計算。

1.3 微分求導法

在尺寸鏈計算中，由函數式表示的尺寸鏈可利用微分法來求解，該方法的本質是通過微分來求解函數的極值。在使用微分法進行尺寸鏈計算之前，必須建立封閉環與各組成環之間的函數關系。

設封閉環與各組成環的函數關系為:

式中:f—封閉環;

n—組成環個數。

當各個組成環xi之間相互獨立時，對函數全微分，得:

假設x1，…，xi－1，xi+1，…xn為已知常數，xi為未知變量，求出xi在其取值范圍內的增減性，從而判斷其取何值時使得f 取極值。

當各組成環相關時，對式(3)進行求偏微分，即:

假設x1，…，xi－1，xi+1，…，xn為已知常數，xi為未知變量，求出xi在其取值范圍內的增減性，從而判斷其取何值時使得f 取極值。

由以上可知，在封閉環與各組成環之間的函數關系比較明確的情況下，可通過微分求導法求出封閉環的極值，從原理上看，該方法是通用的。但是，對于由方程組組成的尺寸鏈，對方程組進行求導不僅繁瑣，并且各組成環的增減性不容易判斷，因此該方法就存在著一定的局限性。

1.4 Monte-Carlo 方法

a)Monte-Carlo 方法的特點

Monte-Carlo 方法被稱為隨機模擬方法(random simulation)，也被稱為隨機抽樣(random sampling)方法。該方法主要用于解決確定性的數學問題和隨機性問題［9-10］，是一種獨具風格的數值計算方法。它的理論基礎來源于概率的大數定理和伯努利定理，其優點以及與其他方法的不同點可歸納如下:

1)Monte-Carlo 方法及其程序結構簡單，只需要產生符合要求的隨機數，通過重復抽樣，求得平均值即可。

2)收斂的概率性和收斂速度與問題的維數無關，Monte-Carlo 方法可適用于多維問題的求解。Monte-Carlo方法的收斂是概率意義下的收斂，其收斂速度比一般數值方法的收斂速度要慢得多。

3)Monte-Carlo 方法通用性強，適用范圍廣，在求解問題時受條件限制的影響小。

b)計算機偽隨機數的產生和抽樣方法

在應用Monte-Carlo 方法模擬某問題的求解過程時，需要產生各種概率分布的隨機變量。服從［0，1］分布的隨機變量是最簡單、最基本并且是最重要的隨機變量也成為隨機數。其他分布的隨機變量的抽樣是通過隨機數來實現的。

引信零部件的尺寸分布可看作服從以下三類分布［5］:

1)均勻分布。服從［a，b］分布的偽隨機數可由式(4)生成:

2)正態分布。當影響因素具有確定的公稱值和公差范圍時，便可認為該影響因素服從正態分布。生成該類分布的表達式為:

其中:tN01，j=)，i=1，2，…，n。

3)瑞利分布。服從瑞利分布的隨機數可由式(6)生成:

隨機抽樣的方法一般有直接抽樣法、舍選抽樣法、復合抽樣法、復合舍選抽樣法、近似抽樣法和變換抽樣法六種。針對不同的影響因素分別進行隨機抽樣，通過計算便可獲得相應的尺寸分布結果。

需要補充說明的是，在生成具有正態分布特性的隨機數后，需要對不在公差范圍內的數值進行剔除，以保證計算結果的正確性。

2 應用Monte-Carlo 方法計算尺寸鏈的實例

現在以兩個例子具體說明Monte-Carlo 方法在尺寸鏈方程組中的實際應用。

2.1 求某引信體上端退刀槽斜面的寬度

a)問題的引出

某引信體及其相關尺寸如圖1 所示。在該引信體零件中，需要計算引信體上端退刀槽斜面的寬度，為檢查風帽與引信體安裝正確性提供數據。該斜面的寬度是通過一組方程組來表示的，無法通過簡單的極大值極小值法來進行計算，也無法通過微分求導法和概率統計法得出結果。

圖1 引信體及相關尺寸

為此，首先通過數學上的幾何關系，以引信體的中心軸線為y 軸來建立一個直角坐標系，并分別將以上各個相關尺寸用a，b，c，d，e，f，g，h，α 來表示，引信體坐標系建立如圖2 所示。

圖2 引信體零件坐標系

在該坐標系中，圖中標注的A，B，C，D 四點的坐標分別表示為(－a/2，e－d)，(－c/2，h)，(－b/2，e－df+g)，(x，y)。根據A，B 兩點坐標與半徑r1，r2之間的數學關系以及C，D 兩點的關系，通過建立關系式，得到如下方程組:

方程組中的l 即為所要求解的退刀槽斜面的寬度。

通過上述所列方程組可以發現，一個看似很簡單的問題，已經轉化為一系列的解釋式，并且該方程組無解析解。在此，將Monte-Carlo 方法應用到該問題的求解過程中。

b)計算的基本假設

在應用Monte-Carlo 方法進行尺寸鏈計算前，需要對引信體零件相關尺寸作以下假設:

1)不考慮有關零件表面形狀誤差影響;

2)不考慮有關零件同軸度誤差的影響;

3)假設零件的相關尺寸均服從正態分布，取其公差帶中心為散布中心，公差帶寬度取6σ，這相當于取工藝能力系數為1(引信生產工藝能力系數一般為1.2～1.5);

4)所有尺寸均以抽樣105次作為整體，以保證所得結果與實際相符合。

c)Monte-Carlo 法求解

通過上述假設，將Monte-Carlo 方法與MATLAB 編程相結合對該問題進行求解，程序流程圖如圖3 所示。

圖3 程序流程圖

在模擬的過程中，先隨機生成1.2×105個隨機數，將不在尺寸公差范圍內的隨機數剔除，然后再從符合要求的隨機數中隨機抽取105個隨機數代表整體，得出如圖4 所示頻數直方圖。

圖4 退刀槽斜面寬度尺寸頻數分布直方圖

注:正態分布抽樣的置信度為99.73%。

d)計算結果分析

通過上面的抽樣分析可以知道:

1)在置信度為99.73%的情況下，引信體上端退刀槽斜面寬度的最大值為1.238 mm，最小值為1.01 mm，平均值為1.12 mm，其分布曲線基本上呈正態分布。

2)在方程組沒有解析解的情況下，Monte-Carlo 方法給出了較為精確的計算結果，并且可以看出，在隨機抽樣得出的結果中，該尺寸值的取值是隨機的，并且大量分布在1.05 mm～1.2 mm 之間。

2.2 風帽收口前的軸向間隙

a)問題分析

某引信安裝風帽后的裝配簡圖如圖5 所示，在本文中只給出引信體與風帽兩零件裝配的示意圖。

圖5 某引信安裝風帽后的裝配簡圖

在求解該問題時，可應用極大值極小值法;但通過計算發現，極大值極小值法計算結果的最小值為負值，與實際不符合。究其原因，是尺寸鏈計算中連續用了兩個通過其他尺寸鏈計算出來的引用尺寸。根據誤差傳遞可知，計算的結果是不準確的。下面根據2.1 條中所用方法來建立幾何關系方程，坐標系的建立與2.1 條相同，如圖6所示。

圖6 裝配圖坐標系

方程的建立過程在此不再贅述。在該問題中，要解決的是求出圖中B，C 兩點的坐標，其中已知B，C 兩點的坐標分別表示為(－i/2，y2)，(x1，e－d－f)。

將兩點的坐標分別代入方程 (x－a0)2+(y－b0)2=，即可求得兩點的坐標值，問題的解即為:

b)基本假設

根據引信體與風帽的裝配情況，現作以下假設:

1)不考慮有關零件表面形狀誤差影響;

2)不考慮有關零件同軸度誤差的影響;

3)假設風帽軸線與引信體軸線重合;

4)假設零件的相關尺寸均服從正態分布，取其公差帶中心為散布中心，公差帶寬度取6σ;

5)所有尺寸均以抽樣105次作為整體，以保證所得結果與實際相符合。

c)Monte-Carlo 法求解

應用Monte-Carlo 方法求解得到的結果如圖7 所示的頻數直方圖。

圖7 風帽收口前軸向間隙尺寸頻數分布直方圖

d)計算結果分析

通過抽樣分析可知:

1)風帽收口前軸向間隙的最大值為0.989 6 mm，最小值為0.435 5 mm，平均值為0.684 2 mm，其分布曲線基本上呈正態分布。

2)在極大值極小值法計算中，最小值為負值，并且最大值與最小值的差值比較大，而該結果中的最小值都大于0.4 mm，因此是與實際結果是相符合的。

3 結論

通過上述計算可知:

1)Monte-Carlo 方法通過隨機抽樣模擬得出的值，可準確的反映尺寸鏈封閉環的尺寸分布情況，并得出在尺寸鏈封閉環的極值。

2)Monte-Carlo 方法可解決尺寸鏈計算中遇到的復雜的尺寸鏈計算問題，通過大量的模擬運算，得出問題的統計特征值。

3)Monte-Carlo 方法結合MATLAB 程序設計可計算組成環除正態分布外的其他分布規律的尺寸鏈，程序簡單，節省時間。

［1］王克沖.用蒙特卡羅法解尺寸鏈問題［J］.華東工學院學報，1985，(01):53-59.

［2］魏傳禮，殷仁龍，劉之生.尺寸鏈的計算及數值解法［J］.華東工學院學報，1989，(01):52-59.

［3］王磊，丁紅宇，楊東拜，等.基于Monte Carlo 方法求解多環尺寸鏈的分析研究［J］.制造業信息化，2002，(01):28-30.

［4］聞泉，王雨時，陳會光.引信機構旋轉偏心的蒙特卡羅模擬［J］.探測與控制學報，2007，29(3):72-75，79.WEN Quan，WANG Yu-shi，CHEN Hui-guang.Simulation for Spin Eccentricity of Devicein Fuze with Monte Carlo Method［J］.Journal of Detection ＆Control，2007，29(3):72-75，79.

［5］陳會光.引信球轉子機構研究［D］.南京:南京理工大學，2007.

［6］陳松濤.蒙特卡洛法在尺寸鏈求解中的應用研究［C］.江蘇省計量測試學術論文集，2005:374-379.

［7］引信設計手冊編寫組.引信設計手冊［M］.北京:國防工業出版社，1987.

［8］李柱主編.互換性與測量技術基礎［M］.北京:計量出版社，1985.

［9］徐鐘濟.蒙特卡羅方法［M］.上海:上海科學技術出版社，1985.

［10］朱本仁編.蒙特卡羅方法引論［M］.濟南:山東大學出版社，1986.