基于靈敏度分析的區間不確定性穩健設計

許煥衛 李沐峰 王 鑫 胡 聰 張遂川

電子科技大學機械與電氣工程學院，成都，611731

0 引言

在工程實際中，復雜技術裝備在設計和優化時，目標函數和約束條件通常都不是線性的。更為困難的是，一些信息通常情況下是匱乏的、不確定和不精確的，如溫度、應力、零件形狀尺寸、操作方式和運行軌跡等的變化，以及建立數學模型時由于認知所限帶來的誤差等。這些不定因素的存在往往導致復雜技術裝備的性能對不確定因素更加敏感，性能波動幾率大大增加，進而影響產品的質量。另外，在設計和優化的過程中往往涉及多種因素，如果全部考慮，不僅會使問題復雜化，還會浪費大量資源。由此，有必要篩選出對產品性能影響較大的因素著重考慮，同時適當忽略影響較小因素的不確定性。

穩健設計理論因具有抗干擾的屬性，現已被廣泛地應用在各領域中[1-3]。穩健優化主要包括兩大類：第一類是概率穩健優化，該類方法需知道輸入參數的概率分布，概率穩健優化被應用到可靠性優化[4]、協同優化[5]、分析目標級聯(ATC)策略中[6]，然而，這些穩健設計優化僅適用于有連續目標函數和約束函數的單目標優化問題，且往往是一種特定的情形，如文獻[6-9]；第二類是以區間分析方法為主的非概率穩健優化[10-13]，因其只需獲得輸入參數的取值區間，而無需確切的分布，故區間分析有很強的適用性。

靈敏度分析[14](sensitivity analysis，SA)的最終目標是識別輸入的可變性對輸出變異性的貢獻。阮文斌等[15]將基于方差和基于失效概率的全局靈敏度分析方法用到復合材料結構中，成功分析了隨機輸入變量對復合材料模型輸出響應量(最大位移和臨界強度比方差)的貢獻大小以及對兩個響應失效概率的影響；張義民等[16]提出了單自由度振動系統的可靠性靈敏度分析方法，放松了對隨機參數的分布類型和激勵類型的限制；CANNAV[17]在總結現有靈敏度分析方法的基礎上提出了一種基于靈敏度的模型選擇準則，并將該準則成功應用到火山源模型的量化擬合中；HAMEL等[18]提出了一種基于多目標的區間不確定下靈敏度分析的改進設計方法，可以在確保可行性的條件下對隨機變量進行輕微調整；邱志平等[19]利用區間擴張理論及其性質定義了結構區間的相對和絕對靈敏度，可在不求導的條件下通過區間運算求得變量靈敏度值，降低了靈敏度求解難度。

由上述可以看出，靈敏度分析被廣泛地應用到各領域中，但是傳統靈敏度分析的精確度依賴于數據的完整程度。本文為解決產品優化設計過程中信息少、涉及因素多等難點，利用穩健設計及靈敏度分析和區間不確定理論的優越性，提出基于靈敏度分析的區間不確定性穩健設計方法。該方法的優勢在于：不僅能在設計之初篩選出對產品性能影響較大的因素，而且還可以選擇性地忽略影響較小因素的不確定性，進而取得性能函數穩健的效果。

1 靈敏度分析

靈敏度分析[20-21]能夠反映某一或某幾個參數變化或外界噪聲對系統輸出的影響程度。

1.1 全局靈敏度分析——Sobol’法

全局靈敏度分析就是綜合考慮輸入變量對性能函數的整體影響，進而掌握非單調、非線性、非疊加的性能函數整體特性。在全局靈敏度分析方法中，Sobol’法[22]是目前常用的方法。設空間單元體為Ω(k)，輸入為k維，表示為

Ω(k)={x|0≤xi≤1;i=1,2,…,k}

(1)

Sobol’法的關鍵步驟在于將性能函數f(x)準確地轉化為

(2)

式(2)中共有2k個子項，采用多重積分進行分解。式(2)中，f0是常數，其余各項對其所包含的每一個因素的積分為0，即

(3)

式(3)的各個子項彼此正交，即

(4)

式(2)中的分解唯一，且各階子項均可通過多重積分運算得到，即

(5)

式中，x-i為去除xi之后的其他變量；x-(ij)為去除xi和xj之后的其他變量。

同樣可類比得到其他高階子項。從而模型f(x)的總方差為

(6)

則式(2)中的各階偏方差可表示為

(7)

i1,i2,…,is=1,2，…,k

對式(2)在整個Ω(k)域先平方后積分，再結合式(4)可得

(8)

那么s階的靈敏度為

(9)

此處Si是因素xi的一階靈敏度系數，反映了因素xi對性能函數f(x)的主要影響程度，Sij(i≠j)為二階靈敏度系數，用來表述兩參數交叉對性能函數的共同影響程度。依此類推，S1,2,…,k反映了k個因素之間的彼此交叉影響。

總效應指數是指因素xi對性能函數f(x)整體影響程度，常用來評價單個參數對性能函數的全部影響。

由式(8)可知

(10)

1.2 Sobol’ 法計算簡化公式

Sobol’法概念雖然簡單，但求解卻十分困難，因此可由蒙特卡羅積分法進行化簡，則式(6)、式(7)中的f0、D及Di，可簡化為

(11)

x(-i)m={x1m,x2m,…,x(i-1)m,x(i+1)m,…,xkm}

(12)

式中，n為蒙特卡羅估計的采樣數；xm為Ω(k)空間的采樣點；上標(1)、(2)代表x的兩個n×k維采樣數據。

2 不確定性分析

工程實際中各種不確定性因素的變化將導致復雜技術裝備的性能不穩定，甚至會導致重大的損失。傳統方法通過提高制造精度等手段來保證產品質量，但因代價過大往往不可取。不確定性優化方法主要有隨機規劃[23]和模糊規劃[24]兩類，然而，在實際應用中，由于信息不完整使得這兩類方法有較大的局限性。區間數優化因其所需信息少，操作便捷，故與工程實際有優良的匹配性，成為一種新的不確定性工程優化方法。

2.1 區間不確定性分析

根據區間分析[25]，區間模型可以定義為

AI=[AL,AR]={x|AL≤x≤AR,x∈R}

(13)

式中，上標L、R 、I分別表示區間下界、區間上界、區間,當AL=AR時，區間退化為一實數a。

區間還可定義如下：

AI=〈Ac,Aw〉={x|Ac-Aw≤x≤Ac+Aw}

(14)

其中，Ac和Aw分別為區間AI的中點和半徑，它們的關系如圖1所示。

圖1 區間的幾何描述Fig.1 The geometric description of interval

定義γ(AI)為區間AI的不確定性水平，表達式為[26]

(15)

2.2 區間可能度及區間可能度轉換模型

因區間表示的是一個范圍，因此需要進行相應轉換以判斷大小及優劣。文獻[27]給出了一種區間可能度關系的計算式：

P(AI≤BI)=

(16)

在區間BI退化為一實數b的情況下，有如下新的可能度關系：

(17)

同理,當區間AI退化為一實數a時，有如下可能度關系：

(18)

圖2、圖3、圖4分別為式(16)、式(17)、式(18)所對應的幾何描述。

圖2 區間AI和區間BI所有可能的位置關系Fig.2 All possible position relationships of interval AI and BI

圖3 區間AI和實數b可能的位置關系Fig.3 The possible position relation of interval AI and real b

圖4 實數a和區間BI和可能的位置關系Fig.4 The possible position relation of real a and interval BI

兩個區間的位置關系確定后，基于區間可能度的確定性優化模型為[27]

(19)

(20)

上述穩健模型可以保證在絕對滿足fI(X)≤VI的情況下，找到不確定因素對目標函數影響最小、性能最穩健、可靠的優化設計解。

3 區間不確定性穩健設計

當前的工程問題日趨復雜，通常需考慮眾多因素，且已知的設計信息較少。這種情況下，傳統優化方法因考慮太多因素而導致計算繁雜且效率不高。本文通過整合全局靈敏度分析、不確定分析、穩健設計方法，構造了基于靈敏度分析的區間不確定性穩健設計框架。具體步驟如下：

(1)根據設計信息建立傳統優化數學模型。

(2)對目標函數f(X)和重要的約束條件g(X)運用Sobol’理論進行靈敏度分析，求出第i個設計變量xi對應的目標函數和重要約束條件的靈敏度值Si。

(3)若Si≤ε(ε值根據實際工程需要而定)，則對應的設計變量xi可設置為常數a，否則對其進行區間分析，確定其取值區間AI，令xi=AI。

(4)利用區間分析改進常規優化模型，利用式(20)建立區間可能度穩健設計模型，從而達到魯棒效果，采用雙層嵌套理論[27]求解所建立的區間可能度模型。

(5)判斷得到的穩健解，如滿足目標函數輸出區間與理想值之間的最大誤差δ<0.1的條件，則程序結束，否則返回步驟(3)重新分析不靈敏項。

該框架可以有效解決靈敏度分析中對隨機變量信息完整度的高要求問題，又可以在區間不確定分析中，降低影響較大因素因波動太大造成的影響。該優化框架可以避開影響較小因素的干擾，重點考慮影響較大的因素，有效降低設計變量維度，因此能夠大大降低模型復雜度，有效縮短計算時間。另外，因為該模型最終輸出為一功能區間，因此還可以在乏信息、乏數據的產品設計初期當做預測模型來預測產品的性能，為進一步地優化提供參考。圖5為基于靈敏度分析的區間不確定性穩健設計框架及求解流程圖。

圖5 流程圖Fig.5 The flow chart

4 工程實例

該工程實例來源于文獻[28]，設計車輛的某型轉向機構，示意圖見圖6。各設計變量相互獨立且服從正態分布。設計變量為軸距L、主銷中心距B、轉向梯形臂長l、梯形底角θ。不確定參數為運動副間隙r1～r4。運動副間隙r1～r4在相應區間上相互獨立且正態分布。相關設計信息如表1所示。

圖6 考慮間隙的汽車轉向機構示意圖Fig.6 Steering Trapezoidal Structure Considering Motion Pair Clearance表1 設計信息Tab.1 Design information

設計變量L(m)B(m)l(m)θ(°)標準差σ0.020.010.010.05均值μ3.652.30.2567不確定參數r1(mm)r2(mm)r3(mm)r4(mm)變化區間(0.5,1.5)(0.2,0.8)(0.5,1.5)(0.2,0.8)

根據相關設計要求,各個設計變量的邊界約束條件為：3.5 m≤L≤3.7 m，2.2 m≤B≤2.5 m， 0.125 m≤l≤0.5 m；根據設計信息可知，在整個轉向過程中轉向機構的運動誤差不能超過3°，并將其設為極限函數，要求其可靠度R≥0.99。當內側轉向輪的轉角為β(通常小于20°)時，外側車輪偏轉角

α0=f0(β)=arccot(cotβ+B/L)

(21)

根據圖6所示的轉向機構，實際的外向車輪偏轉角

(22)

將各相鄰桿件間的配合間隙值代入式(22)，就可以得到實際的外向車輪偏轉角。

4.1 目標函數確定

轉向機構的轉向精度是衡量轉向機構性能好壞的一個重要指標，本文將其作為目標函數。在內側車輪轉角β從最小角度0°轉到最大角度20°的過程中，外側車輪的理想轉向角與實際轉向角之間應盡可能保持一致，即

(23)

式中，α0i為理論值。

在滿足設計要求的前提下，式(23)可以保證實際值的波動在盡可能小的同時最接近目標值，也就是保證目標函數的方差和均值同時盡可能小。

4.2 重要約束條件

當外側車輪的轉角誤差過大時，輪胎的磨損將加劇，為此有

g1(x,r)=|αi-α0i|≤3°

(24)

式中，x為設計變量；r為不確定性參數，體現在相應的桿件中。

在轉向機構中，橫拉桿和轉向梯形臂之間的夾角是一個變化的值，該夾角的最小銳角被稱為最小傳動角。最小傳動角過大將會導致轉向機構的各個桿件之間出現“死點”， 最小傳動角過小又會導致各個桿件受力過大，進而對桿件的強度要求也會相應提高，因此最小傳動角的計算公式如下：

g2(x,r)=

(25)

其中，γmax為最大極限轉角。計算得到的最小傳動角范圍為30°～50°。

轉向梯形角約束如下：

(26)

除上述3個約束外，還有設計變量及不確定性變量的邊界約束。

建立如下一般設計模型：

s.t.g2(x,r)-3°≤0

30°≤g2(x,r)≤50°

0°≤g3(x,r)≤5°

2.2 m≤B≤2.5 m 3.5 m≤L≤3.7 m

0.125 m≤l≤0.5 m

r1∈(0.000 5,0.001 5)mr2∈(0.000 2,0.000 8)m

r3∈(0.000 5,0.001 5)mr4∈(0.000 2,0.000 8)m

接著對目標函數f(X)、重要約束條件g2進行靈敏度分析(g1可包含在目標函數中，g3涉及的參數過少不必進行靈敏度分析),結果如表2所示。

表2 靈敏度分析Tab.2 Sensitivity analysis

綜合分析表2數據可以看出，當ε=0.05時，對目標函數而言，因SL+∑Sri=0.03，小于0.05，故在目標函數中均可忽略其不確定性而設為常值；同理在約束函數g1中因SB+SL+∑Sri=0.025 6，小于0.05，則設計變量B、L及不確定因素r1～r4的不確定性均可忽略而設為常值。表中某些參數靈敏度值為0是因為其值太小而被忽略。

利用區間分析對常規優化設計模型進行改進，可得如下基于靈敏度分析的區間不確定性穩健設計模型：

s.t.P(g1(x,r)≤3°)≥0.99

P(g2(x,r)≥[30°,50°])≥0.9

P(g3(x,r)=[0°,5°])≥0.9

2.2 m≤B≤2.5 m 3.5 m≤L≤3.7 m

0.125 m≤l≤0.5 m

r1∈(0.000 5,0.001 5)mr2∈(0.000 2,0.000 8)m

r3∈(0.000 5,0.001 5)mr4∈(0.000 2,0.000 8)m

為檢驗該模型可行性，設計兩套方案與文獻[28]所提方案及初始解做對比：方案一目標函數和約束條件中L、r2、r4為定值，即設L=3.65 m，r2=r4=0.000 5 m；方案二目標函數中L及不確定因素r1～r4均為定值，這里設L=3.65 m，r1=r3=0.001 m，r2=r4=0.000 5m。采用雙層嵌套理論及粒子群算法對本文實例取自于文獻[28]設計模型進行求解，得到的優化結果如表3所示。

4.3 結果分析

為驗證本文所提方法，將誤差考慮在內，即將兩方案中非靈敏項變差的影響考慮在內時，得到內側角β變化時(0°～20°)相對應的理論外側轉角α0、穩健設計轉角α及本文所求轉角區間，結果如表4所示。由表4可知，內側轉向角β從0°到20°變化的過程中，本文兩種方案結果均在理想值附近。比較這兩方案的最大誤差，δ1max=0.070、δ2max=0.074，均小于設定值0.1，因此模型的解是合理的；不確定性水平γ1max=0.004、γ2max=0.006，說明在去除了不靈敏項的變差影響后，運用本文所提模型依然能夠得到高質量解。對比兩方案的求解效率，方案一由于不考慮3個變量的不確定性，求解時間平均為6.37 s，而方案二去除了5個參數的不確定性，平均求解時間為5.25 s，縮短了13.34%，說明此模型可以在有效減小計算量、簡化模型復雜度的前提下提高效率；從解的質量上看，所忽略變量的總靈敏度值小于規定值時，解的精度和穩健性并不會由于舍棄這些參數而發生質的變化。由此在兼顧計算效率和計算精度的前提下，可以綜合考慮工程需要與模型復雜度等因素選擇去除不靈敏變量的數量。此工程實例驗證了本文所提方法的正確性，且作為一個預測模型是完全滿足實際工程要求的。

表3 優化結果Tab.3 Optimization results

表4 結果對比Tab.4 Result comparison

5 結論

(1)利用全局靈敏度分析——Sobol’法降低模型復雜度，利用區間可能度理論量化不確定因素，進而建立了基于靈敏度分析的區間不確定性穩健設計方法。

(2)將本文所提方法應用到汽車轉向機構的優化設計中，通過兩種優化結果的對比分析可知，所提方法在模型精確度、復雜度、穩健性、求解效率等方面都具有一定的優勢。由于最終給出的結果是一個功能輸出區間，故在乏信息、乏數據的條件下，該工程穩健設計預測模型完全可以滿足工程實際，為設計人員后續的進一步優化設計提供了參考。

(3)本文并沒有嚴格定義可以舍棄其不確定性的靈敏度范圍，也沒有嚴格從數學角度推導當參數靈敏度很小時舍棄其不確定性的正確性，這些需要設計人員根據實際工程優化的具體需求而定。