基于通用生成函數的結構可靠性優化策略

周金宇 朱達偉 王志凌

1.金陵科技學院機電工程學院，南京，211169 2.江蘇理工學院機械工程學院，常州，213001

0 引言

不確定性是客觀世界的重要特征，普遍存在于工業產品的設計、制造和使用過程。工程問題中，加工尺寸、外載荷、材料等往往存在不確定性，設計時若不考慮這些不確定性，可能得到不可靠的設計結果[1]。近年來，計算技術飛躍發展，考慮不確定性的優化設計方法逐步應用于航空航天、遠洋深海、高速機車、特種裝備等工程領域。其中，基于概率理論和隨機不確定性的可靠性設計優化(reliability-based design optimization，RBDO)是理論較成熟、應用最廣泛的考慮不確定性的優化設計方法[2]。

根據優化流程的迭代結構，可將RBDO方法分為雙循環法、單循環法和解耦法。雙循環法是解決RBDO問題的最基本方法，具有外層設計變量優化與內層可靠性分析相互耦合的雙層嵌套結構。雙循環法計算成本高、效率低，因此學者相繼提出了單循環法和解耦法。單循環法利用近似等效條件替代可靠性約束，避免優化循環內嵌的可靠性分析循環。為進一步提高優化求解效率，LI等[3]基于可靠設計空間，將RBDO問題完全轉換為確定性優化問題；ZHOU等[4]提出的基于順序逼近的兩階段RBDO方法顯著減小了計算量。對于解耦法，DU等[5]引入偏移向量，將設計變量優化與可靠性分析解耦分離，提出了序列優化與可靠性評估(sequential optimization and reliability assessment, SORA)方法；TORII等[6]提出一種適用于各類可靠性分析方法的通用解耦方法；HAO等[7]提出一種增強步長調整(enhanced step length adjustment, ESLA)迭代算法和基于二階可靠性分析的逐步增強順序優化與可靠性評估方法(SSORA-SORM)。

工程結構RBDO通常面臨的兩大技術瓶頸是復雜功能函數的調用成本高和可靠度指標的求解精度低。一方面，實際工程計算中含復雜功能函數的問題時，往往需借助有限元分析等成本高昂的學科工具。近年來，研究人員借助多項式響應面(RSM)、人工神經網絡(ANN)、支持向量機(SVM)、克里金(Kriging)模型等代理模型，通過試驗樣本的高效訓練建立分析變量與復雜功能函數之間的簡明映射，在一定程度上解決了復雜功能函數調用成本高的難題，因此功能函數調用次數不總是算法評估的核心指標。另一方面，針對工程結構RBDO中常見的變量非正態分布(如指數分布、均勻分布、多峰分布)和功能函數非線性場合，可靠度指標求解的常用方法(FORM、SORM)存在較大誤差，甚至因誤差不可控、迭代不收斂導致方法失效。Monte Carlo模擬(MCS)法是高精度求解方法，但計算成本極高，一般難以勝任工程結構RBDO的可靠性分析。因此，在可控成本下實施可靠度指標的高精度求解，已成為工程結構RBDO亟待攻克的現實難題。

上世紀80年代USHAKOV[8]提出通用生成函數(universal generating function, UGF)概念以來，UGF法被引入工程領域并取得豐碩成果。文獻[9-11]在多狀態系統研究領域中應用并進一步發展該方法，使之逐步成為系統概率分析的重要工具。LISNIANSKI[12]基于變量離散化思想，利用UGF計算連續狀態系統可靠性指標的界后，UGF法逐步擴展到具有連續型隨機變量的結構體系概率分析[13]。

筆者在結構RBDO中引入UGF，針對隨機變量為非正態分布和功能函數非線性的RBDO問題，提出一種基于UGF的可靠性設計優化策略。通過系列響應面建立分析變量與概率指標的動態映射，消除傳統雙循環法的內層循環；利用UGF法替代傳統矩法完成響應面動態更新過程中的可靠性分析。該方法能在可控的計算成本內顯著提高隨機變量/參數非正態、功能函數非線性場合的工程結構RBDO的求解精度。

1 結構可靠性分析的UGF法

1.1 連續型隨機變量的UGF

對于連續型隨機變量s，累積分布函數及概率密度函數分別為FS(s)和fS(s)，將s在其定義域(smin,smax)內近似均勻離散化為m個點，分別記作s1、s2、…、sm。離散點si(i=1,2,…,m)對應的概率為

(1)

式中，δ為離散步長，δ=(smax-smin)/m。

這樣，可根據離散數據集{(si,pi)|i=1，2，…，m}定義連續型隨機變量s的UGF：

(2)

式中，離散值si為隨機變量的第i個狀態值，si=smin+(i-0.5)δ;z為UGF模型中的默認字符，僅用于表示函數結構，無實際含義。

1.2 基于UGF的可靠性分析

對于涉及n維連續型隨機向量S=(S1,S2,…,Sn)的工程結構，可靠性分析可利用式(1)、式(2)獲得S第j個分量的UGF：

(3)

設結構的功能函數為G(S)，G(S)>0時，結構可靠，否則失效。對結構進行可靠性分析時，需對各隨機變量UGF進行復合運算，獲得描述總體性能分布的結構UGF：

(4)

其中，?G為復合算子，由復合算子根據結構物理內涵定義不同變量UGF之間的運算規則。

結構UGF包含總體性能分布信息，可用來計算各類可靠性指標[9]。設結構的性能分布用UGF表示為

(5)

則依據該UGF的概率項進行條件求和，得到結構可靠度

(6)

其中，M為結構UGF的離散狀態組合總數，可考慮對性能值接近的狀態組合進行同類項合并，且通常M≤mn；ψ(·)為條件求和算子；I(·)為示性函數，當其自變量大于0時取1，否則取0。

由式(6)可知，基于UGF的可靠度計算原理是，對各隨機變量離散狀態組合進行條件枚舉。借助UGF復合算子內稟的普適性、遞推性、可分性、互換性等優異特性[9]，該方法適用于任意分布的隨機變量，以及任意形態的功能函數。

1.3 連續型隨機變量的非均勻離散化

UGF法已在多狀態系統分析中得到成功應用，因存在組合爆炸的壁壘，因此在連續狀態結構系統分析中的應用受限。對于連續狀態的結構系統，通過低密度離散化將連續型變量描述為UGF往往達不到預期精度，高密度離散化將在后繼復合運算中消耗較高計算成本，甚至導致組合爆炸，而借助同類項合并的增效技術在較多場合下無效。因此，本文針對結構RBDO問題，對隨機變量進行非均勻離散化，在減少離散狀態組合數的情況下保證可靠度指標的求解精度。

非均勻離散化的基本思路是當給定的離散狀態數m因計算成本原因而限定為較小值時，選擇結構RBDO當前迭代步的最大概然失效點(most probable failure point, MPP)為敏感點SMPP。以該點為中心，依照等比級數對各隨機變量進行敏感點密集、邊緣點稀疏的離散化，使SMPP鄰域內的離散點密集化，以保證極限狀態敏感區的概率分析精度。隨機變量Si在敏感區可實現的最小離散步長為

(7)

式中，simax、simin分別為Si的最小值和最大值；SMPP(i)為SMPP的第i個分量；q為等比級數的公比，通常取1

這里的SMPP為極限狀態超曲面上距離隨機空間的均值點μS最近的點，通過求解下面的數學模型可得敏感點SMPP：

(8)

2 基于UGF的RBDO方法

RBDO求解的代表性數學模型為

(9)

2.1 指標函數的響應面模型

分析RBDO數學模型不難發現，不確定性功能函數Gq(d,X,P)與設計點D密切相關，每個迭代步的可靠性分析均基于當前設計點Dk。因此，Dk和該點對應的當前可靠度指標βk之間存在如下映射關系：

Dk→βk

(10)

確定指標函數的響應面模型時，可選用不含交叉項的一次和二次響應面，這兩種響應面模型結構簡單、運算量小且精度滿足求解要求，各自表達式分別為

(11)

(12)

式中，nd為設計變量的個數；Dk,i為設計點Dk第i個分量；a、bi、ci為待定系數。

初始響應面的構建依賴于初始設計點D0鄰域的N對輸入輸出數據：

B0=[D0,1D0,2…D0,N]
β0=(β0,1,β0,2, …,β0,N)

式中，B0為初始迭代步的響應面輸入數據矩陣；D0,i為D0鄰域的第i個試驗點向量；β0為初始迭代步的響應面輸出數據向量；β0,i為對應輸入數據的第i個可靠度指標值。

對于一次響應面模型和二次響應面模型，N=2nd+1。

首先，確定初始設計點D0與初始迭代步的響應面輸入數據矩陣B0。具體做法是：通過等距抽樣在設計空間中生成若干點樣本，對這些樣本進行約束函數評估，選定位于可行域且目標函數值較小的點為初始設計點D0。以此為中心，通過下列規則確定B0中的試驗點向量：

(13)

其中，α為調節系數，可取0.05～1。θ(t)(t=1,2,…,n)為擴展向量，其分量若與確定性設計變量相對應，則取值為初始設計點到該分量定義域最近邊界的距離；其分量若與不確定性設計變量相對應，則取值為相應隨機變量的標準差。

(14)

式中，βT為許用可靠度指標。

2.2 指標函數響應面的動態更新

圖1 基于UGF的RBDO算法流程

3 算例分析

為比較不同算法的求解精度，計算最優點處可靠度指標的相對誤差：

(15)

式中，βMCS為在不同算法最優點處采用Monte Carlo模擬求得的可靠度指標。

3.1 算例1

本例為通用數值算例[14]，涉及2個確定性設計變量d1、d2和2個隨機參數X1、X2，且功能函數為非線性。RBDO數學模型為

(16)

其中，X1、X2均服從正態分布，其均值μ1、μ2分別為5和3，標準差σ1、σ2分別為1.5和0.9；許用可靠度指標βT=2.32；隨機參數的離散狀態數取30。

根據式(11)構建迭代開始時的響應面和迭代收斂時的響應面，如圖2所示。不同算法的目標函數值迭代過程如圖3所示，優化結果見表1，其中，N為功能函數調用次數；算法1為基于FORM的雙循環法；算法2為SORA法；算法3為本文所提方法，其可靠性分析環節由UGF法完成；算法4為本文所提方法的變形，其可靠性分析環節由MCS法完成。

(a)迭代開始時響應面

圖3 不同算法目標函數迭代過程

由表1可以看出，所提算法與傳統算法均可收斂到最優解。算法3的精度高于傳統方法(算法1和算法2)，算法4的精度更高，但對功能函數的調用過多(5×106次)。因此，與傳統方法相比，算法3在計算成本可控的條件下獲得了更高精度，算法可行。值得一提的是，本例中的隨機變量均為正態變量，所以傳統算法也能有效實施并獲得一定精度。

表1 不同RBDO算法的算例1優化結果

3.2 算例2

本算例涉及設計隨機變量X3和X4，其中，X3服從標準差為0.6的正態分布，X4服從標準差為0.6的均勻分布，兩設計隨機變量的均值μ3和μ4為不確定性設計變量。RBDO數學模型為

(17)

其中，許用可靠度指標βT=2。設計隨機變量的離散狀態數取30，不同算法的優化結果見表2。

由表2數據可以發現，優化求解問題含高度非正態(如均勻分布)的隨機變量時，若采用矩法完成可靠性分析，則當量正態化處理極易導致顯著誤差[15]，因此算法1、算法2均無法收斂到最優解，方法失效。算法3和算法4均可收斂到最優點，且精度較高。兼顧求解效率，算法3最優，解決了傳統算法對高度非正態隨機變量無法求解的普遍難題。

表2 算例2不同RBDO算法優化結果

3.3 算例3

短梁結構如圖4所示。矩形截面的高為h、寬為b，自由端受到雙軸彎矩M1、M2和軸向力F作用，材料的屈服強度為Y。

圖4 短梁結構示意圖

確定性設計變量d=(b,h)，隨機參數P=(M1,M2,F,Y)。各隨機參數相互獨立，其中，M1、M2、Y服從正態分布，F服從均勻分布，各自統計信息見表3。該結構RBDO的數學模型為

表3 隨機參數統計信息

(18)

其中，許用可靠度指標βT=3，隨機參數的離散狀態數取10，不同算法的優化結果見表4。

分析表4數據可知，RBDO問題中的隨機變量增多且存在非正態(如均勻分布)隨機變量時，算法1可以收斂，但在最優解處的可靠度不滿足概率約束的要求，求解精度極低；算法2無法收斂，方法失效；算法3以較高精度收斂到最優解。算法3為抑制多變量UGF復合運算引發的組合爆炸，將隨機變量的離散狀態數減為10，但由于算法采用了非均勻離散化技術，故仍可保證較高的求解精度。

表4 不同RBDO算法的算例3優化結果

3.4 算例4

本算例對單級齒輪傳動減速箱的局部結構進行可靠性優化設計，要求在滿足齒輪強度可靠度R≥0.99的條件下結構總質量最小。已知傳動比為3.2，指定齒寬系數d1、模數d2、主動輪齒數d3、主動軸的最小直徑d4、從動軸的最小直徑d5為確定性設計變量。齒輪傳動的輸入轉矩T服從均值為108.8 N·m、標準差為8.704 N·m的均勻分布；齒輪的彎曲疲勞強度極限σF,lim服從均值380 MPa、標準差30.4 MPa的對數正態分布；齒輪的接觸疲勞強度極限σH,lim服從均值500 MPa、標準差44 MPa的對數正態分布。該優化問題的數學模型為

(19)

其中，確定性設計變量d=(d1,d2,d3,d4,d5)，隨機參數P=(T,σF,lim,σH,lim)，許用可靠度指標βT=2.33，隨機參數的離散狀態數取20，不同算法的優化結果見表5。

表5 不同RBDO算法的算例4優化結果

由表5所示的優化結果可知，由于隨機變量為均勻分布和對數正態分布，且同時存在2個概率約束，因此算法1、算法2均無法得到收斂解而失效，算法3可較好地收斂到最優解。與算法4相比，算法3的求解精度滿足工程需求并具有更高的求解效率。

4 結論

(1)針對通常的結構RBDO，本文方法比傳統方法具有更高的求解精度和可控的計算效率。

(2)RBDO問題涉及非正態隨機變量、高度非線性功能函數時，傳統方法存在求解精度低或無法收斂的劣勢；本文方法的可靠性分析由UGF法完成，不受隨機變量非正態、功能函數高度非線性的影響，具有較高的魯棒性。

(3)RBDO問題存在高維隨機空間時，所提方法可借助非均勻離散化技術，在保持隨機變量離散狀態數不變的條件下，提高了優化求解的精度。