熊芬芬,陳江濤,任成坤,張立,李澤賢
1 北京理工大學 宇航學院,北京 100081
2 中國空氣動力研究與發展中心,四川 綿陽 621000
復雜工程產品(例如飛行器、船舶、汽車等)高水平、高效率的開發設計對國民經濟及國防事業的發展有著舉足輕重的作用。為了縮短設計周期、降低開發成本、滿足產品不斷提升的性能需求和更新換代需要,20 世紀60 年代中期以來,有限元分析(FEA)和計算流體力學(CFD)等計算機仿真技術和優化方法被廣泛應用于復雜系統的設計。模型確認、不確定性分析和優化設計是復雜工程產品基于仿真的開發設計中涉及的兩大關鍵內容。前者主要在于保證仿真模型的高保真度,從而為可靠的產品性能分析和設計提供保障;后者旨在通過設計提高產品性能在各種不確定性下的可靠性和穩健性。對于二者而言,需要解決的首要問題皆為不確定性傳播(uncertainty propagation, UP),這也一直是工程優化領域重要的理論課題之一。然而,隨著工程系統設計的復雜化、多學科化,仿真分析模型計算規模和計算量增長,響應函數高維、強非線性,給不確定性傳播帶來“維數災難”、精度低、可靠性差等諸多難題。
作為一種高效的概率不確定性傳播和量化理論,混沌多項式(polynomial chaos, PC)方法[1-2]由于其堅實的數學基礎和良好性能,近些年在學術界和工業界受到廣泛關注。混沌多項式方法實質上相當于將隨機變量表示為一組正交多項式的加權和,構建一個隨機代理模型,不確定性傳播就直接在這個代理模型上進行。它能夠對具有任意分布類型的隨機變量實現較為精確的近似,且理論上當條件滿足,即可獲得指數收斂速度。而且,一旦混沌多項式模型構建完成,輸出響應的統計矩、失效概率以及概率密度函數都能非常方便地得到。相比于傳統的蒙特卡羅仿真(Monte Carlo simulation, MCS)方法,混沌多項式方法在保證精度的前提下,可大幅降低計算量;相比于傳統的一階[3]、二階[4-5]可靠性分析方法,對非線性函數具有更高的不確定性傳播精度,且無需函數的導數信息,應用起來更加靈活。目前,混沌多項式方法在機械[6]、土木[7-8]、材料[9]、電力電子[10-11]、汽車[12]、航空航天[13-19]、船舶[20-23]、控制[14,24-25]等領域都得到了廣泛研究和應用,在航空航天和船舶領域應用尤為廣泛。
由于構建混沌多項式模型的計算量隨著不確定性輸入維數的增加呈指數增長,高維下面臨嚴重的“維數災難”問題,這也是目前阻礙混沌多項式方法在工程問題中廣泛應用的最大障礙之一。因此,如何解決或緩解混沌多項式的“維數災難”難題,降低計算量,一直是學術界的研究熱點。另一方面,混沌多項式方法基于概率理論,僅能處理隨機不確定性,然而實際問題中存在大量由于知識或數據不足而導致的認知不確定性,例如在CFD 建模中由于認知不夠,湍流模型的建立存在模型不確定性和相關參數選取的不確定性。因此,如何充分利用混沌多項式的優勢,將其擴展為能處理認知不確定性,是關注的另一熱點。本文將首先給出不確定性傳播的基本概念及其在模型確認和不確定性優化中的作用,然后系統地介紹混沌多項式方法的構建原理及其各類變種,接著針對當前混沌多項式應用中面臨的“維數災難”、計算量大、無法直接處理認知不確定性等問題介紹相應的解決策略,最后對混沌多項式的研究進行展望。
不確定性傳播也稱不確定性分析,是研究各種不確定性對產品系統性能(泛指系統輸出)的影響規律的方法。簡而言之,不確定性傳播就是在給定輸入的不確定性信息下,如何估算輸出響應的不確定性信息。不確定性主要分為兩大類:隨機不確定性和認知不確定性,前者無法控制或減少,主要基于概率理論進行研究,而后者主要由認知或數據不足導致,可以減少。目前有多種研究認知不確定性的方法,例如模糊[26]、區間[27]、證據[28-29]理論等。由于概率方法發展成熟,具有堅實的數學基礎,很多情況下結合貝葉斯定理也可以處理認知不確定性,因此目前概率方法應用最為廣泛。本文介紹的混沌多項式方法建立在隨機概率空間,主要用于處理隨機不確定性。圖1展示了概率不確定性傳播的基本概念,從數學上具體描述為:在隨機輸入x=[x1, ···,xd]存在不確定性的情況下(此時x1, ···,xd的不確定性可以用其概率密度函數、累積分布函數、或均值和方差描述),計算輸出響應y的不確定性信息,包括均值、方差、失效概率、概率密度函數(PDF)等。

圖1 概率不確定性傳播示意圖Fig. 1 Probabilistic uncertainty propagation
對于基于仿真的復雜產品性能分析和優化設計,要想保證精度,構建高可信度的仿真模型是關鍵。然而,由于物理過程的復雜性及人們的認知偏差,物理建模與數值模擬始終存在不確定性。以CFD 計算為例,其存在幾何模型、模型假設、模型參數(例如湍流模型封閉系數)、迭代方法等各類不確定性,從而嚴重影響數值模擬結果的可信度, 使得決策需承擔很大風險。以飛行器結構有限元分析數值模擬為例,傳統的確定性有限元分析方法并未考慮有限元建模與試驗測量中普遍存在的不確定性,在實際應用中僅以變量的均值對問題進行分析與描述。或者以圍繞均值附近的“安全系數”這類概念粗略估計問題的隨機性,基于工程師個人的工程經驗得到“安全”或“不安全”這類結論。若以此作為飛行器結構設計的關鍵性決策依據,顯然是不盡合理的。若不確定性無法得到有效量化,使用與真實結果存在較大差異的數值模擬進行設計,極有可能導致真實系統達不到預期的性能要求,引入潛在風險。對于數值模擬模型的確認,由于存在大量模型不確定性和參數不確定性,所以必須首先分析不確定性對數值模擬結果的影響,實現有效的不確定性傳播,對仿真結果進行不確定性量化,評估其可信度,進而進行模型修正或重選,確保仿真結果的高可信度。事實上,模型修正通常涉及不確定性傳播的反問題求解,不確定性傳播也是其中的核心之一。不確定性量化分析在模型確認總體技術路線中占據重要地位,該問題的研究是近年來國內外研究人員關注的熱點。事實上,美國機械工程師協會在其計算固體力學驗證和確認(V&V)項目指南中,已將不確定性量化分析列為模型確認中必不可少的關鍵研究內容[30]。圖2 以CFD 模型確認為例,展示了不確定性傳播在其中的作用,在模型確認過程中要反復調用不確定性傳播模塊進行不確定性量化,這也是其中最耗時的部分。實際上由于不確定性高維、數值模擬較為耗時,“維數災難”是不確定性傳播面臨的最大難題,因此提高其計算效率尤為重要。

圖2 不確定性傳播在模型確認中的作用Fig. 2 UP in model validation
由于混沌多項式的高效性,其在數值模擬模型確認中應用非常廣泛,尤其對于CFD 數值模擬等通常較為耗時的問題,混沌多項式方法能大幅降低計算量,因此其在CFD 模型確認的不確定性量化中應用尤為廣泛。Meldi 等[31]將混沌多項式方法用于空間演化混合層流大渦模擬(LES)的不確定性量化,并分析了不同模擬參數的敏感性。王瑞利等[32]將非嵌入式混沌多項式方法應用于爆轟模型,對平面爆轟問題和散心爆轟問題進行了不確定性分析;趙輝等[33]發展了非嵌入式混沌多項式方法,研究了湍流模型系數的不確定性對RAE2822 翼型跨聲速繞流模擬的影響和材料物性參數的不確定性對燒蝕熱響應預測的影響,為工程中多變量不確定性量化問題提供了有效的解決方案;劉全等[34]將混沌多項式方法應用于拉式流體模型,對Sod 激波管問題進行了不確定度量化;Enderle 等[35]將非嵌入式混沌多項式方法應用于湍流噴霧燃燒模擬,通過對所用的雷諾平均Navier-Stokes 模型進行不確定性分析實現了對變量降維,進而實現了模型優化。Schaefer 等[36-37]將混沌多項式方法用于跨聲速壁面束縛流湍流模型封閉系數的不確定性量化和敏感度分析,提出了一種基于混沌多項式方法的用于工業級規模氣動分析的不確定性量化框架,考慮了湍流模型封閉系數、來流工況、網格收斂誤差等不確定性。
確定性條件下的優化技術已經成功運用到諸多工程設計問題中,但是設計條件的變化,例如載荷、材料特性和操作環境的變化,往往使得工程系統存在很大的不確定性。以往由于數學處理和計算速度等方面的原因,通常將這些不確性量作為確定量處理,導致產品的設計性能對這些不確定性因素非常敏感,穩健性低,或在不確定性下產品性能無法滿足約束,設計失效,極有可能導致災難性的后果。因此,產生了不確定性下的優化設計。根據設計理念的不同,可將其分為穩健優 化 設計(robust design optimization, RDO)[38-39](見圖3)和基于可靠性的優化設計(reliabilitybased design optimization, RBDO)[40-41](見圖4),也有將二者結合的,這就是所謂的穩健可靠性優化。通過在設計過程中考慮設計變量、設計參數、設計決策和系統分析模型等不確定性因素的影響,來降低不確定性對系統性能的影響。如圖3 所示,穩健優化所得最優解,雖然相對確定性優化損失了一定的產品性能,但是卻對不確定性的敏感程度更低,具有更強的穩健性。如圖4 所示,基于可靠性的優化所得最優解離約束邊界相對較遠,雖然一定程度上可能會導致設計的保守,但是能保證設計始終落在可行性區域內,系統是安全的;而確定性最優解由于設計變量或參數存在不確定性,最優解很大可能違反約束。

圖3 確定性最優和穩健性最優Fig. 3 Deterministic and robust optimums

圖4 確定性最優和可靠性最優Fig. 4 Deterministic and reliability-based optimums
圖5 展示了不確定性優化的流程圖。由圖可知,不確定性傳播是其中必不可少的模塊,用于計算目標和約束函數的均值、標準差或失效概率。不確定性傳播的精度和效率幾乎決定了整個優化的精度和效率[42]。如上所述,對于復雜工程系統設計,仿真分析模型通常較為耗時、非線性程度高且高維,因此計算量大是其不確定性傳播面臨的主要難題。

圖5 不確定性傳播在不確定性優化中的地位Fig. 5 UP in design optimization under uncertainty
由于混沌多項式的高效性,其在不確定性優化設計中也得到了大量應用。蔡宇桐等[19]針對壓氣機葉片制造中的加工誤差,應用非嵌入式混沌多項式進行不確定性量化,得到了增強葉片氣動性能穩健性的優化方向;胡晚亭等[8]將混沌多項式方法應用于Winkler 地基沉降計算,得到參數變化對地基沉降可靠性的影響,進而明確了設計方向;李冬琴等[22]利用混沌多項式法分析多維隨機不確定性因素對船舶優化方案的影響,完成了船舶多學科穩健設計優化研究,有效減少和避免了船舶設計優化方案失效的可能性。Lie 等[43]將混沌多項式方法與禁忌搜索算法結合起來對天線陣列反問題進行穩健優化,相較于拉丁超立方抽樣方法極大提高了計算效率并且保證了穩健性能的計算精度;Mandur 等[44]在模型參數不確定性傳播的穩健優化問題中使用混沌多項式方法,展現出了比蒙特卡羅法更強大的計算精度和計算效率優勢;Huang 等[45]針對探測器火星大氣進入的軌跡穩健優化問題,采用混沌多項式方法獲得認知不確定條件下軌跡性能的近似解析函數,建立了一個雙環嵌套的穩健優化模型。
混沌多項式方法分為嵌入式(intrusive)和非嵌入式(non-intrusive)[46]。嵌入式混沌多項式方法主要應用于動力學系統的不確定性傳播,其主要思想是將不確定性源進行混沌多項式模型表達,然后將其代入到動力學微分方程,利用Galerkin投影,將隨機微分方程轉換為一組更高維的以混沌多項式模型系數為狀態量的確定性微分方程組,通過求解該微分方程組得到混沌多項式系數,從而完成不確定性傳播[47]。由于嵌入式混沌多項式方法需要進入系統方程內部,對方程進行擴維修改,不可避免地會引入數值誤差。另一方面,當前的系統分析模型通常都經過了反復標定驗證和確認,對其進行模型內部調整可能性非常小。此外,隨著仿真軟件的大量應用,復雜工程產品設計分析基本都涉及黑箱型響應函數。因此,嵌入式混沌多項式方法應用非常有限,非嵌入式混沌多項式方法將響應函數看作一個黑箱,僅關注輸入和輸出的映射關系,構建輸出響應的混沌多項式代理模型,從而實現不確定性傳播。由于簡單易于操作,相比嵌入式混沌多項式,其研究應用更加廣泛。本文主要介紹非嵌入式混沌多項式方法。
以響應函數y=g(x) (x=[x1,···,xi,···,xd])為例,混沌多項式理論下輸出響應表示為p階截斷的混沌多項式模型。

式中:bi為 待求的混沌多項式系數; Φi為正交多項式;正交多項式的總項數為P+1=(d+p)!/d!p!;ξ=[ξ1,···,ξi,···,ξd], 為標準隨機向量,其中 ξi與原隨機變量xi存 在一定的轉換關系,具體與xi和 ξi的分布形式有關。
一旦混沌多項式系數bi求出,就可直接在混沌多項式模型上運行MCS,得到輸出響應y的隨機概率特性,包括均值 μy、 方差 σy、概率分布等。在式(1) 的基礎上,也可得到低階統計矩關于混沌多項式系數的如下解析表達:

式中,E 表示求期望,是一個運算符。
對于混沌多項式方法,主要涉及兩方面內容:正交多項式 Φi的構建和混沌多項式系數bi的求取,其中后者涉及響應函數的調用,是計算量消耗的主要來源。
混沌多項式方法最早由Wiener[48]提出,以埃爾米特(Hermite)正交多項式作為基函數,被稱為Wiener 混沌多項式。根據Cameron-Martin 理論,它對正態分布型輸入具有指數收斂速度,但是對于其它的分布類型,則收斂速度明顯降低。這是因為Hermite 正交多項式的權函數剛好與正態分布的概率密度函數具有相同形式。利用上述特點,Xiu 等[2]根據Askey 方案中各隨機分布類型概率密度函數與正交多項式權函數一一對應的關系,針對不同隨機輸入類型,采用相應類型的一元正交多項式作為基函數,對其進行直接張量積操作(tensor product)得到由多元正交多項式Φi(ξ)(i=0,1,···,P)構建的混沌多項式模型,將混沌多項式理論擴展到了廣義混沌多項式(general PC,gPC)方法,使其對Askey 方案中所包含的正態、均勻、指數、貝塔(Beta)、伽馬(Gamma)等多種隨機分布類型均具有指數收斂速度。當輸入變量不滿足Askey 方案中所列的5 種分布時,則通常將其首先轉換為標準正態分布,然后選取Hermite正交多項式作為基函數。gPC 方法極大地拓展了混沌多項式理論的應用范圍,使其迅速被大量應用。
對于一般的常見分布類型,gPC 方法能夠得到較好的不確定性傳播結果。但是在實際工程應用中,存在諸多不屬于甚至與Askey 方案相差較遠的分布類型,例如質量、擴散系數、剛度系數等物理參數均服從對數正態分布。此時,通常需要利用不確定性變換將其轉化為Askey 方案中的分布類型,即:將原隨機變量xi(i=1,···,d)轉換為Askey 方案中所示的某類標準隨機變量 ξi,而變換自然會引入一定的誤差,降低不確定性分析的精度。針對隨機分布類型的復雜性和多樣性特點,不少研究提出針對任意的隨機分布類型,自行構建正交多項式。Witteveen 等[49]提出了基于Gram-Schmidt 正交分解的混沌多項式方法;Zhang 等[50]和Xu 等[51]提出了基于斯蒂爾吉斯(Stieltjes)過程的混沌多項式方法。這些研究分別利用Gram-Schmidt 正交分解和Stieltjes 過程,構建相對于各維隨機輸入分布的最優一元正交多項式基函數。在此基礎上,類似于gPC 通過張量積操作,構建p階 截 斷的 多 元正 交 多項 式 Φi(x)(i=0,1,···,P)。相比于gPC,這些混沌多項式方法適用范圍更加廣泛,可應對任意隨機輸入分布類型,無需利用變換,提高了收斂速度。
上述所有混沌多項式方法都建立在已知隨機輸入完整的概率分布函數基礎之上,而在實際工程應用中隨機參數的信息可能以各種形式存在,如離散的原始數據樣本,尤其對于復雜系統,往往由于價格高昂、耗時太長等,難以得到其完整的概率密度函數。此時,上述基于各維隨機輸入完整概率分布函數的混沌多項式方法不再適用。為此,Oladyshkin 等[52]提出了一種數據驅動混沌多項式(data-driven PC)方法,該方法可應對各種類型的隨機分布以及離散的原始數據(無需概率分布函數),并具有良好的收斂性和精度。Wang 等[53]在此基礎上進一步將投影法(Galerkin projection)引入到數據驅動混沌多項式中計算混沌多項式系數,提出了相應的高斯節點和權值的計算方法,進一步提高了其靈活性。數據驅動混沌多項式方法的主要思想為,根據各維隨機輸入變量的離散數據或概率分布函數,計算其統計矩,利用正交多項式的正交性,推導矩匹配(moment-matching)方程,進而通過匹配隨機輸入變量一定階次的統計矩,完成最優一元正交多項式的構建。需要注意的是,當離散數據不足時,一元正交多項式的構建會引入一定的誤差。
考慮到在實際工程應用中,常有隨機輸入變量相關的情況。例如,結構中的材料屬性與疲勞屬性[54]、飛行器氣動噪聲與外表面分布的隨機載荷[55]等。針對上述問題,多數研究學者采用變換方法,如正交變換[56]、Rosenblatt 變換[57]和Nataf變換[58]等,將相關隨機輸入變量轉換為相互獨立的標準正態隨機變量。但是這些變換方法均屬于非線性變換,會導致變換后的響應函數呈現強非線性特征,特別是當相關隨機輸入變量服從復雜的非正態概率分布或響應函數非線性較強時,在不確定性分析中會引入非常大的計算誤差。此外,這些轉換方法均要依賴隨機變量完整的分布函數,當相關隨機輸入變量以離散數據形式存在時,上述變換方法顯然均不再適用。為此,在數據驅動混沌多項式方法基礎上,Lin 等[59]提出了一種能夠直接處理輸入變量相關性的混沌多項式方法,其能避免變換帶來的誤差,且能應對離散數據形式的相關輸入,擴展了混沌多項式方法的適用范圍。Paulson 等[60]針對嵌入式混沌多項式方法,利用Gram-Schmidt 正交變換,提出了一種考慮輸入變量相關性的方法。隨后,Wang 等[61]將該方法應用到了隨機潮流。這類考慮相關性的數據驅動混沌多項式方法依然采用矩匹配方法建立一元正交多項式基,只不過在求解矩匹配方程的時候引入了相關隨機輸入變量混合矩的概念,進而考慮變量的相關性。
這類自行構建正交多項式的方法顯然提高了混沌多項式方法的適用范圍,具有更加廣闊的應用前景,但是構建正交多項式基的過程難免引入近似或數值等誤差。尤其是針對隨機變量為復雜分布(如雙峰或多峰),數據驅動混沌多項式方法可能需要匹配很高階次的統計矩,矩匹配方程的求解極易出現奇異,導致正交多項式基的構建精度難以保證。gPC 方法由于實現方便,穩健性相對更高,且理論上增加混沌多項式階次到一定程度即可應對任意形式的輸入分布,因此,目前依然是應用最為廣泛的混沌多項式方法。
混沌多項式模型的構建首先必須給定混沌多項式模型的階數。事實上,混沌多項式模型階數p對不確定性傳播的計算量和精度具有重要影響。常用的做法是進行階數的收斂性測試,先從較低混沌多項式模型階次出發進行不確定性傳播,增加階次繼續進行不確定性傳播,同時考慮先前樣本的重復利用,直到相鄰兩次的結果(如輸出響應y的均值和方差)變化不大,則認為當前混沌多項式階數滿足要求。但是對于某些特殊問題,非線性程度高,分布函數甚至呈現非規則雙峰的情況,則需要高達p=15 的階數才能獲得滿意的精度[62]。對于工程不確定性分析,往往涉及非線性黑箱型響應函數,對于混沌多項式模型階次的確定并非易事。Hu 和Youn[63]利用留一法(leaveone-out, Loo)評估混沌多項式的精度,構建混沌多項式階次的循環,不斷增加混沌多項式階次,直至滿足精度要求為止,從而決定最終的混沌多項式階數。目前這方面的研究并不多,主要原因是對于一般的實際工程系統,通常系統輸出響應主要受各維輸入變量以及其低階交叉項的影響,高階交叉項的影響較小[64],因此,混沌多項式模型階數p=2 或3 就能滿足工程實際需求,出現高階混沌多項式模型的情況較少。
混沌多項式系數求取是混沌多項式進行不確定性傳播的關鍵,目前主要有投影法和回歸法兩類方法。
2.4.1 投影法
利用Galerkin 投影方法,將式(1)兩邊同時依次投影到各正交多項式 Φj(ξ)上,得

根據內積的定義,并利用正交多項式的正交性,整理式(4),得


2.4.2 回歸法

Xiong 等[62]在上述SRSM 基礎上,進一步提出了加權隨機響應面方法(Weighted SRSM),通過引入樣本權值的概率,考慮其在概率空間的分布特性,在相同樣本的情況下,提高了不確定性傳播的精度。與確定性下的響應面相似,回歸中樣本的性能很大程度上決定了響應面的好壞,因此決定了SRSM 不確定性傳播的精度。關于抽樣方法,Hosder 等[69]從精度和收斂性方面,對隨機抽樣、拉丁超立方設計(Latin hypercube design, LHD)和Hemmasy 抽樣方法進行了綜合的分析比較,并推薦用兩倍于未知混沌多項式系數個數的樣本(即N=2(P+1))可以得到比較滿意的結果。Isukapalli[67]提出在Wiener 混沌多項式方法中運用Hermite 積分節點,也就是Hermite 正交多項式的根,作為回歸樣本來求解系數。在該方法的基礎上,可從Askey 方案中所列正交多項式出發,分別求解p+1 階正交多項式的根,對這些根(一維空間)進行直接張量積操作,得到多維空間的全因子設計樣本,然后在這些樣本點上進行最小二次回歸,得到混沌多項式系數。這種抽樣策略使得樣本點大多數集中在概率空間中的高頻率區域,類似于重要性抽樣的原理,因此提高了混沌多項式系數的估算精度。由于這種抽樣方法對一維樣本采用直接張量積,得到的多維樣本的組合數為(p+1)d,呈指數增長。出于計算量的考慮,目前仍局限于啟發式地選取其中部分樣本進行回歸,通常選取概率空間出現頻率較大的樣本點。也有研究提出利用單項求容積法則(monomial cubature rules, MCR)來產生回歸樣本,構建混沌多項式模型,這就是所謂的PC-MCR 方法。該方法運用MCR 產生樣本點,由于樣本點數目少,可以全部用來估算混沌多項式系數,不僅保證了不確定性傳播的精度,而且與前面提到的高斯積分點方法相比大大減少了所需的樣本數。
回歸法中由于涉及矩陣求逆等運算,針對高維問題可能會較為繁瑣耗時,因此回歸法通常適用于低維問題。隨機抽樣、拉丁超立方設計、高斯積分點、MCR 抽樣方法在混沌多項式系數計算方面具有相似的精度,MCR 方法雖然計算效率方面最為可觀,但由于在目前的科學計算軟件中并未見到成熟算法,而且可選取的樣本形式有多種,具體選何種較為主觀,因此其應用相對較少。
混沌多項式方法具有較高的精度和效率,但是其計算量通常隨著隨機輸入的維數呈指數增長,高維下存在“維數災難”問題,這也是目前混沌多項式在實際應用中面臨的最大難題。關于該問題,目前產生了諸多方法去緩解“維數災難”。
為了解決維數災難,最常用的方法是減少全階混沌多項式模型中的正交多項式項數。其中,雙曲線截斷(hyperbolic truncation)策略[70]是最常用的一種,其主要思想為:一般的實際工程系統的輸出響應主要受各維輸入變量以及其低階交叉項的影響,受高階交叉項的影響較小[64],從而可以直接去除某些多項式項,達到降低計算量的目的。
傳統的全階混沌多項式模型從一元正交多項式基出發,利用直接張量積構建多元正交多項式,當混沌多項式模型階數為p時,需要滿足

式中:ai為 多元正交多項式 Φi中第i維變量對應的一元正交多項式的階數;p為混沌多項式模型的階數,也就是式(1)中的正交多項式Φi(i=0,1,···,P)的最高階次。
在直接張量積下,式(1)中正交多項式的總項數為P+1=(d+p)!/d!p!,可見總項數隨著d的增加呈指數增長。為了降低計算量,考慮以下雙曲線截斷策略,去除部分高階交叉多項式項,以達到減少混沌多項式系數個數、降低計算量的目的。

式中:q為自定義稀疏因子;I為滿足不等式的a=(a1,···,ad)的集合。
以p={3,4,5,6}和q={1,0.75,0.5}為例,在二維空間中上述截斷方法產生的階數組合如圖6 所示,圖中藍色的圓圈表示滿足式(10)中不等式要求的a=(a1,···,ad)的集合。顯然,通過引入雙曲線截斷策略,確實去除了一部分高階交叉項,對降低計算量具有重要作用。而且,相同的階數p下,q越小,去除的高階交叉項越多,計算量降低越明顯。當q=1(圖6中第一行)時,式(10)退化為常規的全階混沌多項式方法。

圖6 不同稀疏因子下一元正交多項式的階數組合Fig. 6 Order combinations of orthogonal polynomial with different sparse factors
除了這類直接去除高階交叉項多項式項的截斷策略,Hampton 和Doostan[71]提出了各向異性多項式階數(anisotropic polynomial order)的概念,在構建混沌多項式模型的過程中,認為 Φi中對應于各維變量的一元多項式的最高階次并非相同,通過引入誤差指數,自適應地搜尋最優的各維正交多項式階數a=(a1,···,ad),使之滿足

可見,當pi=p(i=1,···,d)時,該方法退化為全階混沌多項式模型。實際中,輸出響應關于各維隨機變量的非線性是不同的,因此上述考慮各向異性的多項式階數確定方法能夠在計算量一定的情況下,將更多資源配置于非線性較強的維度上,從而提高精度。
多項式基截斷策略能夠在一定程度上降低計算量,但是在進行截斷的時候稀疏因子的選取較為主觀。如果選取的值過小,雖然可以很大程度降低計算量,但是由于截除了過多的高階交叉項,精度的損失較為嚴重。因此,該方法適合于高階交叉項影響較小的問題。
稀疏重構的主要思想是在全階混沌多項式模型的基礎上,通過去除對輸出響應影響不大的正交多項式項 Φi(i=0,1,···,P),減少混沌多項式系數的個數,從而降低計算量,是當前應對混沌多項式“維數災難”難題最有效的途徑之一。稀疏重構下構建的稀疏混沌多項式模型可表示為

其中,混沌多項式系數的個數為P1+1, 顯然P1+1 ?P+1。
對于稀疏混沌多項式方法,關鍵在于如何去發掘非重要的正交多項式項。Blatman 和Sudret[72]提出了一種自適應算法來自動檢測重要的正交多項式項,認為那些能夠最大程度降低預測誤差的正交多項式項為重要的項并予以保留,從而減少混沌多項式系數的個數。隨后,他們又提出利用最小角回歸(least angle regression, LAR)方法去發掘非重要的正交多項式項,進而將其從全階混沌多項式模型中去除,以降低計算混沌多項式系數的計算量[73]。王豐剛[74]對基于LAR 的稀疏混沌多項式方法展開了研究,通過諸多算例測試發現,隨著維數和混沌多項式階數的增長,全階混沌多項式所需的樣本數量明顯增加,導致混沌多項式求解系數的回歸矩陣出現“病態”,難以求解,然而稀疏混沌多項式則幾乎不受此影響,同時還能保持較高的計算精度。Hu 和Youn[63]針對工程可靠性分析和設計,提出在全階混沌多項式模型的正交多項式集合中,通過引入誤差指數,自適應地篩選誤差指數最大的雙變元正交多項式項來構建稀疏混沌多項式模型。陳光宋等[75]采用最小絕對收縮和選擇算子(LASSO)回歸自動選擇混沌多項式的重要項及其展開系數,進一步由混沌多項式系數解析獲得全局靈敏度系數。Cheng等[76]提出了一種可觀測響應保持同倫下的自同態調制(D-MORPH)算法,用于構建稀疏混沌多項式模型。Diaz 等[77]、Tsilifis 等[78]和陳江濤等[79]分別研究了利用壓縮感知(compressed sensing)方法構建稀疏混沌多項式模型,壓縮感知方法是圖像和信號處理領域興起的新方法,能夠高效地重構稀疏信號,需要的采樣點數目小于自由度的個數。若混沌多項式展開是稀疏的,可通過求解以下優化問題得到:

式中,0 范數表示b中非零元素的個數,式中所涉及的變量符號參見式(7)。考慮到式(13) 為NPhard 問題,求解很難,且實際應用中考慮到測量噪聲的情況,不要求 ψb=Y精確滿足,則式(13)一般變為

在求解過程中,需要指定截斷誤差ε。除了以上通過發掘全階混沌多項式模型中非重要的正交多項式來構建稀疏混沌多項式的方法,也有研究將混沌多項式表示為高維模型表達(high dimensional model representation, HDMR)的架構,認為各變量單獨作用和雙變量共同作用的HDMR函數就已經可以較為精確地描述系統的輸出,忽略交叉高階項,從而構建所謂的稀疏混沌多項式模型,以達到降低混沌多項式模型中正交多項式項數的目的[80]。
LAR、壓縮感知等稀疏重構方法在構建稀疏模型的過程中,需要通過迭代發掘重要的正交多項式項,尤其對于非線性較強且維數較高(d>10)的問題,收斂過程可能非常慢,且收斂過程受所選取的樣本點和問題的非線性程度影響非常大。
3.1 和3.2 節旨在通過減小全階混沌多項式模型中正交多項式的項數來降低計算量。稀疏網格數值積分則通過在投影法計算混沌多項式系數中采取Smoyak 算法,生成積分節點,相對于全因子節點大為減小積分節點的個數,從而緩解“維數災難”難題。關于這方面的研究目前非常多。Winokur[81]針對混沌多項式研究了一種自適應稀疏網格數值積分方法,大大降低了計算量,并將其成功應用于2004 年9 月伊萬颶風穿過墨西哥灣時的海洋環流模型中。Wu 等[82]將稀疏網格數值積分用于混沌多項式,求解了不確定性下的翼型氣動優化。Xiong 等[83]將基于稀疏網格數值積分的混沌多項式方法應用于火箭彈穩健優化。稀疏網格數值積分確實能在一定程度上降低中、低維(d<10)不確定性傳播問題的混沌多項式系數的計算量,根據維數可以預估函數調用次數,而且具有較強的穩健性,是當前緩解維數災難較為可靠的方法,但對于高維問題其應對能力非常有限,計算量依然非常大。
在基于仿真的工程設計中,分析模型(例如CFD和FEA)往往具有高度非線性、計算耗時的特點。采用混沌多項式方法直接基于高精度仿真模型進行不確定性傳播,也同樣面臨計算量大的問題。對于工程系統設計,由于學科分工愈來愈細,分析方法和仿真建模手段逐漸多元化。隨著設計進程的推進,往往伴隨有多種不同精度和計算量的分析模型產生。例如,某小型飛機涉及多個固定翼和螺旋槳的氣動耦合分析,可利用多種多可信度氣動仿真工具實現,包括簡單低階模型(二維渦模型和葉素理論螺旋槳模型)、中精度模型(渦格法和Euler CFD,其中螺旋槳模型均采用激勵盤)、高精度模型(基于雷諾平均Navier-Stokes方程的 CFD)。此外,還有地面試驗數據和飛行試驗數據。通常認為模型的精度越高,耗費代價越大,其可生成的樣本數量也越少。為降低計算量,充分利用多個分析模型,產生了多可信度建模(multi-fidelity modeling)的思想,也稱多模型融合(model fusion),在大量廉價的低精度樣本基礎上,利用少量高精度樣本點為導引或修正,通過融合不同精度和計算量的樣本數據,建立多可信度代理模型(multi-fidelity metamodel),能夠在保證代理模型精度的同時,盡可能地降低計算量。在該方面具有代表性的研究有Kennedy 與O'Hagan提出的基于高斯隨機過程及上述差值(或比值)模型的多層級co-kriging 方法[84]等。
由于多可信度建模理論在降低計算量方面效果顯著,有學者基于混沌多項式方法開展了多可信度不確定性傳播的研究,通過以少量高精度樣本為引導,融合大量低精度樣本數據并建立不同精度模型之間的修正混沌多項式模型,來提高低精度混沌多項式模型的可信度,從而達到降低計算量的目的。Ng 和Eldred[85]提出了基于混沌多項式和稀疏網格數值積分的多可信度不確定性傳播方法;Palar 等[86]提出了基于最小二次回歸的多可信度混沌多項式方法。上述多可信度混沌多項式方法采用加法修正的多可信度建模策略,即在低精度混沌多項式模型上進行加法項混沌多項式修正,雖有效解決了特定的問題,但依然存在不足。例如:Ng 和Eldred 提出的方法需要采集稀疏網格數值積分點處的樣本,其數量和位置都不是任意的;Matteo[87]提出的方法因利用線性回歸法計算混沌多項式系數,雖然在采樣策略上相對靈活,但需要高/低精度樣本點嵌套。為此,Berchier 提出用低精度混沌多項式模型來預測非嵌套樣本點處的響應值,但這必然會引入一定的預測誤差,降低了多可信度混沌多項式模型的精度。Yan 和Zhou[88]針對貝葉斯推理反問題,提出了一種基于自適應抽樣的加法修正多可信度混沌多項式方法。Cheng 等[89]提出了一種基于高斯過程回歸的多層級多可信度稀疏混沌多項式方法。Wang等[90]基于高斯隨機過程,將多層級co-kriging 方法從確定性多可信度建模領域擴展到不確定性量化,構建一種了多可信度混沌多項式方法,研究表明,相較于常用的基于加法修正的多可信度混沌多項式方法精度大幅提高,且該方法能夠應對模型的精度水平為非層次型的情況。除了以上兩類基于加/乘法修正和高斯隨機過程的多可信度混沌多項式方法,還有研究提出利用輸出空間映射(output space mapping)技術實現多可信度混沌多項式模型構建[16],其主要思想為建立低精度混沌多項式模型到高精度混沌多項式模型的映射關系,映射關系最常見的為線性映射。
由于多可信度建模在降低計算量方面具有巨大潛能,而且切合目前產品設計過程中存在多種類型分析模型或數據的現狀,多可信度混沌多項式研究較為活躍。目前,基于加法修正的多可信度混沌多項式方法由于形式和實施簡單,應用最多。基于高斯隨機過程的多可信度混沌多項式方法精度和靈活性高,但由于利用高斯隨機過程,建模中參數估計存在計算繁瑣和穩健性不足的缺點,這也是目前基于高斯隨機過程的建模方法存在的共性問題。整體而言,多可信度混沌多項式方法在降低計算量方面具有巨大潛力,但是與確定性領域的多可信度建模方法相似,需要多可信度分析模型預測趨勢一致,存在泛化能力不足的問題,而且如何配置各個多可信度分析模型的樣本個數以保證不確定性傳播的精度,目前也缺乏科學系統的方法。
穩健優化在工程設計得到廣泛應用,混沌多項式由于精度高被大量用于穩健優化。基于梯度的尋優算法由于效率高在穩健優化中應用很多,而在基于梯度尋優的穩健優化中,需要計算目標和約束函數相對于設計變量的局部敏度信息,也稱設計靈敏度(design sensitivity)。此外,該設計靈敏度還能為設計決策者提供參考,從而在實際中對產品質量波動進行管控。常規的做法是直接利用差分法計算這些設計靈敏度,但這必然會消耗一定的計算量,尤其當設計變量高維時,計算量較大。
混沌多項式用于局部靈敏度分析的研究不多。Ren 等[91]針對基于數值積分計算混沌多項式系數的穩健優化,推導了設計靈敏度,無需調用任何額外的響應函數。該方法被用于基于混沌多項式的翼型氣動穩健優化,相比于直接依賴差分法的穩健優化,可降低高達30%的計算量。考慮到系數的求取基于Galerkin 投影和高斯數值積分,見式(5),其中分子的期望值 E[yΦi(ξ)]計算都將表示為高斯節點上的函數響應值的加權和,最終設計靈敏度的計算歸結到期望值關于節點的偏導?E[yΦi(ξ)]/?xj.ij,其中xj.ij為高斯積分節點。分別針對FFNI 和SGNI,針對不同的節點分布情況,結合拉格朗日插值近似技術,利用節點上的函數值,在不調用任何響應函數的情況下,推導了半解析形式的設計靈敏度。
全局靈敏度分析(global sensitivity analysis, GSA)能夠量化各不確定性因素對系統響應不確定性的影響程度,進而可適當忽略那些影響小的不確定性因素,在不影響精度的前提下降低不確定性傳播的計算量,目前已成為提高不確定性傳播效率的有效途徑之一。全局靈敏度分析方法大體可以分為基于回歸的方法和基于方差的方法[92],其中基于Sobol'靈敏度指數的方差分析法因其簡單有效的特點得到了廣泛應用。
混沌多項式方法將隨機輸出響應表示為一組正交多項式的加權組合,基于混沌多項式模型可解析得到響應方差,從而非常方便地計算Sobol'靈敏度指數。最終的靈敏度指數是關于混沌多項式系數的解析表達,相當于混沌多項式系數計算的副產品。該方法相比于MCS 計算量大為降低,因此基于混沌多項式的Sobol'方差分析法成為一種較為常用的全局靈敏度分析方法。
基于混沌多項式方法的Sobol'靈敏度指數表示為

式中,DPC為總方差,可通過式(3) 計算得到。由式(15) 可知,一旦混沌多項式模型構建好,就可在混沌多項式模型基礎上進行全局靈敏度分析,得到關于混沌多項式系數的靈敏度指數的表達,無需任何額外的函數調用。
Sudret[92]將廣義混沌多項式應用到全局靈敏度分析中,得到了關于混沌多項式系數的解析式的Sobol'靈敏度指數。Palar 等[93]基于加法修正多可信度混沌多項式方法,進行了基于Sobol'靈敏度指數的全局靈敏度分析。Cheng 等[94]基于支持向量機回歸提出了用于全局靈敏度分析的自適應混沌多項式方法。王晗等[95]提出了基于稀疏多項式混沌展開的孤島微電網全局靈敏度分析方法,并將其用于準確、快速地辨識影響系統運行狀態的關鍵輸入隨機變量。卜令澤[96]對基于混沌多項式的全局靈敏度分析方法展開深入研究,為大型復雜結構的靈敏度與可靠度分析提供了新思路。為識別影響自動裝填機構剛度的核心關鍵參數,孫佳等[97]采用基于混沌多項式展開的全局靈敏度分析方法,從32 個自動裝填機構參數中提取出了6 個影響剛度的核心關鍵參數。王娟[98]開展了基于混沌多項式的Sobol'全局靈敏度分析,考慮了輸入變量線性相關性的情況。也有研究首先在全階混沌多項式模型上進行雙曲線截斷,構建稀疏混沌多項式模型,在此基礎上進行全局靈敏度分析,進而發掘Sobol'靈敏度指數較小的非重要變量,在不確定性傳播中不予以考慮,從而實現降維,然后在降維后的變量空間構建階次更高的稀疏混沌多項式,達到了提高精度并降低計算量的目的[85]。
當考慮隨機(aleatory)和認知(epistemic)混合不確定性時,常規的不確定性傳播過程變為一個典型的雙層循環,外層考慮認知不確定性,內層實質為認知不確定性固定于某值情況下的隨機不確定性傳播。由于混沌多項式方法的高效性,近年來出現了不少將混沌多項式方法用于處理隨機和認知混合不確定性傳播的研究工作。其基本思想為:外層利用概率盒(P-box)、證據、模糊、區間以及似然估計等方法考慮認知不確定性,內層進行基于混沌多項式的隨機不確定性傳播,從而避免大量調用耗時的真實響應函數,降低計算量。



圖7 基于概率盒和混沌多項式的混合不確定性傳播Fig. 7 Mixed UP with PC and P-box
證 據 理 論(evidence theory)是 由Dempster 率先提出的,后經Shafer 系統完善,故又稱為Dempster-Shafer 理論[102],它是對經典概率理論的一種擴展,使用概率邊界反映所有可能結果集合冪集的信任度。證據理論通過識別框架、基本可信度分配(basic probability assignment, BPA)、可信度函數Bel()和 似真度函數Pl()這3 個基本概念構成了不確定性建模框架。讀者可參考文獻[29, 94, 103]獲取關于證據理論的具體介紹。利用證據理論做混合不確定性傳播時,基于混沌多項式求取輸出響應的累積信度函數(cumulative belief function,CBF)和累積似真度函數(cumulative plausibility function, CPF)來表征識別框架內系統響應的不確定特性,如圖8 所示。在得到輸出響應y的CBF 和CPF 之后,可進一步建立適合不確定性優化設計的類似于均值和方差的目標準則[104]。

圖8 證據理論下輸出響應y 的CBF 和CPF 曲線Fig. 8 CBF and CPF of output response with Dempster-Shafer theory
模糊數(fuzzy number)是常規實數的一般化,其含義是它不引用一個值,而是引用一組可能的值,其中每個可能的值都有自己的權重(范圍0~1),稱為隸屬函數。在進行混合不確定性傳播時,可以將模糊數看作一種靈活的P-box 的形式,即對每個 αi-cut 的區間進行基于混沌多項式的不確定性傳播[105-106],最終得到系統輸出響應y的隸屬度函數。圖9 為大致過程的示意圖,其中縱軸表示隸屬度函數p(x)值 , α ∈[0,1]是 α-cut 水平。

圖9 基于混沌多項式和模糊理論的混合不確定性傳播Fig. 9 Mixed uncertainty propagation with PC and fuzzy theory
區間模型一般定義如下:

式中,上標I,L,U 分別表示區間、區間下界和區間上界。基于區間數的混合不確定性傳播可轉換為在區間變量的范圍內求解優化問題[107-108]。對于任意的區間變量所在區間上的某個值,固定認知不確定性變量于該值,構建隨機變量的混沌多項式模型,從而可方便地得到輸出響應的均值、方差或失效概率。進一步通過尋優,得到輸出響應均值、方差或失效概率的最大值和最小值。
不同于以上4 種方法,似然方法(likelihood)[109]最終依然采用概率理論來進行不確定性傳播。似然估計允許同時處理點數據和區間數據,它將以稀疏點數據和區間數據存在的認知不確定性輸入變量表示為概率形式的不確定性,通過最大化似然函數,得到認知不確定性變量的概率密度函數(PDF),該PDF 相當于是在當前已知數據下,該變量的平均PDF。概率表征使得該方法可非常方便地應用于各類基于概率的不確定性分析及優化設計理論和方法,而基于概率的方法具有嚴格的理論基礎,發展較為成熟,因此應用起來非常方便。
混沌多項式方法作為一種較為成熟高效的不確定性傳播方法,近幾年得到了長足發展,已經在各個領域得到了廣泛應用。在穩健優化設計領域,由于需要在優化的每個迭代點構建混沌多項式模型,對于高維且非線性較強的優化問題,往往涉及大量迭代,計算量顯著上升,而數值模擬模型確認不涉及該迭代,相比于傳統的不確定性傳播方法,混沌多項式的優勢凸顯,因此混沌多項式方法在數值模擬模型確認方面最具應用潛力。相比其它混沌多項式的變種,廣義混沌多項式方法由于具有相對較強的穩健性,針對不同的輸入分類和函數形式通常能得到較為滿意的結果,是目前應用最廣泛的混沌多項式方法。但其大部分應用基本停留于較簡單的工程問題,問題涉及的規模及維數都不太高(d<10),高維問題鮮有報道。在應對“維數災難”方面,稀疏網格方法穩健性好、泛化能力相對較強,但僅適用于中、低維問題(d<10);稀疏重構策略計算量降低最為顯著,是最具應用潛力的途徑之一,但在適用范圍和泛化能力上還需改進。多可信度混沌多項式是降低計算量的有效途徑之一,但是對于高、低精度模型樣本數目的確定還缺乏較為科學有效的方法,目前基本停留在不斷試湊的模式,且泛化能力較差。此外,實際問題大都涉及黑箱型響應函數,圍繞混沌多項式開發較為通用化的程序軟件,根據實際計算資源預算自動確定有效的混沌多項式階次和樣本點個數,以方便實際工程應用,建立混沌多項式在各類典型問題應用的演示實例,這也是需要進一步探索的問題。實際工程中可能涉及高達上百維的不確定性傳播問題,目前的混沌多項式方法依然會面臨較為突出的“維數災難”問題,這方面還需挖掘新的理論和方法,例如嘗試采用小樣本深度學習技術構建輸入輸出的深度神經網絡模型,取代原系統響應分析,以降低計算量。