P-CS 不確定性量化模型與其性能數據驅動更新方法1)

何佳琦 賈曉璇 吳偉達 鐘杰華 羅陽軍,2)

* (大連理工大學航空航天學院,遼寧大連 116024)

?(中國運載火箭研究院,北京 100076)

引言

實際工程問題往往包含多源不確定性變量,如荷載、結構自身材料屬性(楊氏模量、磁導率等)、結構服役環境(磁場、溫度場等)、結構加工誤差以及傳感器測試誤差等.這些不確定性信息直接或間接影響著工程結構形式設計[1]、結構性能評估與預測[2]以及在役結構損傷識別[3]等工作的開展與決策.建立合適的數學模型來量化這些不確定性變量是開展考慮上述不確定性的結構分析與設計工作的前提與基礎.

傳統處理不確定性變量的方法通常將其量化為概率統計變量[4-6],但概率模型只適用于先驗信息充足且測試樣本容量足夠大的客觀不確定性.對于無充足先驗知識、僅具有少量數據樣本的客觀或主觀不確定性變量,Ben-Hahn[7-8]在20 世紀90 年代放棄傳統概率理論框架,提出了用凸集模型描述變量不確定性的非概率思想.近年來,大量學者針對非概率模型的建立與結構性能非概率測度描述方法做了一系列的研究工作[9-12].

對于具備有限先驗知識與有限數據樣本的情況,一些學者還結合區間理論提出了非精確概率不確定性量化方法,核心思想是用一個概率的區間模型取代精確概率值來表現由于先驗知識與數據樣本不充足而引入的主觀不確定性,并形成了一系列相關方法論,主要包括DS 證據理論[13-14]、P-Boxes 模型[15-16]與模糊概率模型(fuzzy probabilities)[17-18].其中,DS 證據理論基于信任函數和似然函數給出上限概率與下限概率,為不確定性量化尤其是多源不確定信息融合提供了合理且較為完備的理論框架,已廣泛應用于多傳感器監測信息融合[19]、失效模式分析[20]、結構可靠性評估[21]等領域.P-Boxes 模型本質上是通過求解一系列優化問題尋找到非精確概率不確定性變量的邊界概率密度累積函數,進而用這些邊界概率密度累積函數界定變量的不確定性范圍.目前,已有一些學者利用P-Boxes 模型開展多源不確定性量化分析[22]、結構剩余壽命評估[23],以及考慮不確定性的侵入式結構有限元分析[24]等問題的研究.模糊概率模型是模糊集理論與概率模型的結合,隨機不確定性由概率模型描述,而概率模型的不精確性由模糊集隸屬度函數(membership function)[25-26]描述.目前已有一些學者基于模糊概率模型提出了結構可靠性評估方法[18,25-29]與疲勞強度分析方法[30].考慮到結構在服役期間,不確定性因素會隨時間變化,一些學者針對時變不確定性問題開展了研究.Li 等[31]基于P-Box 模型提出了一種用于量化時變不確定性的非精確隨機過程模型;Ping 等[32]基于正交級數展開、稀疏網格數值積分以及混沌多項式展開等數值方法提出了對于非高斯隨機過程的不確定性傳播算法;考慮到時變不確定性的結構失效概率在一個時間段為一區間值,Wang等[33]提出了一種非概率區間過程模型,Meng 等[34]結合概率模型(量化隨機不確定性變量)與凸集模型(量化功能函數在某一時間段的變化范圍)提出了一種求解時變可靠性指標的迭代算法;為避免因大規模求解功能函數而造成時間成本過高的問題,一些學者結合Kriging 代理模型提出了效率較高的失效概率區間計算方法[35-36].

盡管上述各類模型與算法根據各自的適用范圍已廣泛應用于不同領域各類問題中,然而,多源不確定性量化依舊存在下述兩個問題難以兼顧解決:(1)實際工程問題往往既包含可用概率模型量化的不確定性變量,也包含適用于非概率模型量化的不確定性變量,兩者之間存在復雜層級關系;(2)一些不確定性變量在使用過程中隨時間變化,且在使用過程中難以直接測量,需要根據間接性能測試信息在使用工程中更新不確定性量化模型.這些不確定性變量既包括隨時間變化的客觀不確定性變量如環境載荷(冰、雪、風、波浪等)、結構本身屬性(楊氏模量、屈服強度、阻尼系數等),也包括隨著更多相關信息加入而不確定性逐步減小的主觀不確定性變量.對于問題(1),盡管P-Boxes 模型可融合概率不確定性變量與非概率不確定性變量,但由于該模型求解較為復雜,很難基于P-Boxes 模型進行基于性能測試數據驅動的量化更新計算.對于問題(2),一些學者基于貝葉斯理論[37]與證據合成理論[38]分別針對概率模型[39]與非概率模型[40]提出了更新算法.但考慮到兼顧問題(1)與問題(2),目前尚沒有能統一描述多源不確定性且易于性能數據驅動更新的不確定性量化模型.

本文為兼顧解決上述問題,提出一種P-CS(probability-convex set)不確定性量化模型,用以實現多源、多類型不確定性的統一量化以及性能測試數據驅動模型的更新計算.基于該模型,可將概率模型、非概率模型以及非精確概率模型在同一框架下統一表達,本文將分別用一維、二維與三維P-CS 模型詳細展示P-CS 模型的構建過程并討論其概率與非概率特性.進而,基于貝葉斯理論提出P-CS 模型參數在性能數據測試樣本驅動下的更新算法,并分別采用懸臂梁結構與框架結構算例來展示更新計算過程并驗證算法的適用性.

1 P-CS 不確定性量化模型

1.1 P-CS 不確定性量化模型表達形式

當給定隨機變量的概率模型后,凸集模型Ω用來表征: 在某一概率水平 α 下,一簇概率模型的隨機變量分位數Xα相對于隨機變量的分位數的變化范圍.P-CS 模型可定義如下

1.2 一維P-CS 模型構建方法

由表1 可知,只有針對“已知變量可能的概率分布形式,概率參數不確定但有界”的非精確客觀不確定性變量,P-CS 模型的構建區別于傳統的概率模型與非概率模型.對于這種情況,其概率分布為一族函數,如圖1 所示.

表1 多源不確定性變量的P-CS 量化模型Table 1 P-CS models for multi-source uncertainty

圖1 非精確概率模型示意圖Fig.1 Schematic illustration of imprecise probabilistic model

式(3)中,Xα關于變量 μ與 σ 為線性變化,因此

1.3 多維P-CS 模型構建方法

由此,無量綱橢球模型可表達為

圖2 橢球模型旋轉角度示意圖Fig.2 Schematic illustration of the rotation angle of ellipsoidal model

在實際工程中,對于非精確概率變量很難給出清晰的不確定性變量之間的相關性關系,一般只能憑借工程經驗判斷變量之間是否相關.因此,本文將非概率橢球模型中關于變量間的相關性進行模糊處理.考慮到隨著后續測試樣本的輸入,可逐步降低變量間相關性的模糊程度,使非概率多橢球模型逐漸變得精確.故而,初始多橢球模型可構建為包絡超立方盒模型 Ωcubic的最小體積超橢球模型.在幾何學上,可包絡Ni維區間模型的最小橢球,其主軸方向 θ0沿區間長度最大的維度,其第k維度半軸長與變化范圍區間長度dk=存在下述對應關系式

2 性能測試信息驅動的P-CS 不確定性量化模型更新算法

結構在服役過程中往往會產生新的測試數據,這些新的測試數據的加入使得主觀不確定性減小,不確定變量之間非概率相關性模糊程度降低,表現為P-CS 模型中表征不確定性范圍的多橢球凸集模型體積減小與橢球方向角的改變.實時性能測試數據驅動P-CS 模型更新就是根據實時性能測試數據對P-CS 模型中表征不確定性范圍的多橢球凸集模型參數進行更新計算,本質上就是在計及實時性能測試數據的基礎上尋求當前最可能的參數值組合.為了尋求可能性最大的參數值組合,這里把參數值的可能性用參數“信度”這一概念加以量化.

對標準化參數取值范圍做m個劃分Ai(i=1,2,···,m),Ai滿足如式(13)所列出的條件,將X落入Ai的概率定義為信度,記為Pl(X∈Ai)或簡略為Pl(Ai) .

圖3 二維焦元示意圖Fig.3 Schematic illustration of 2-dimensional focus elements

多個焦元Bi,i=1,2,···,n的并集的信度等于每個焦元信度的和,即

2.1 似然信度定義

其中,fPDF(·) 為標準正態分布概率密度函數,τi表示第i個多維橢球模型參數.

2.2 基于貝葉斯理論更新P-CS 模型

考慮到基于不同多橢球按照優化列式(17)計算得到的性能可能取值范圍會存在重疊現象.而根據式(16)的定義,當一個樣本值落入到重疊部分,不同的焦元會獲得相同的似然信度.不同焦元對應的多橢球模型體積不同,顯然基于體積越大的多橢球模型可獲得范圍更大的性能取值域.因此,在式(16)對似然信度定義的基礎上,對多橢球模型的體積進行了“懲罰”,即定義焦元對應的多橢球模型體積與焦元似然信度成反比例,在P-CS 模型更新過程中,對于給定P-CS 模型 C(m),測試樣本為 ν=ν0的似然信度可進一步定義為

其中,Vm為P-CS 模型 C(m)中多橢球模型體積.

基于貝葉斯方法,可根據性能測試樣本數據更新多橢球模型參數取值劃分的信度.對于第j個參數 τj,其信度更新貝葉斯公式如下

其中,Pl(ν=ν0) 表示參數τj取任意值,測試樣本為ν=ν0的似然信度,根據式(15),信度具有可加性,所以

對無初始信度信息的問題,參數先驗信度可定義為每一個焦元等信度,即

在更新計算過程中,可按下式計算更新后橢球模型參數的參數值,即

2.3 更新算法流程

考慮到焦元數量M=R2N-1,為提高計算效率,本文采取一種根據參數后驗信度分布,逐漸減小參數值取值范圍的多次數更新策略.在讀入樣本數據后,更新算法流程圖如圖4 所示,具體步驟詳述如下.

圖4 更新計算方法流程圖Fig.4 A flow chart of the updating method

步驟1:根據式(18)確定各初始參數取值范圍,并分別將各初始參數取值范圍做R個劃分,建議R≤6,并根據式(23)定義初始先驗信度分布,根據式(17)計算每一個焦元對應的性能函數值的期望變化范圍Ie,進而根據式(21)做第k=1 次更新計算,得出參數后驗分布.

步驟2:若在第k次更新計算后,第j個參數的第i個劃分的后驗信度則在參數取值范圍中將區間去掉,得到取值范圍

步驟3:將取值范圍做R個劃分,并根據式(23)定義初始先驗信度分布,根據式(17)計算每一個焦元對應的性能函數值的期望變化范圍Ie,進而根據式(21)做更新計算,得出參數后驗分布.

步驟4:重復步驟2 與步驟3,直至在第K次更新計算中,每一個參數的劃分區間都足夠小,則終止多次更新計算.可定義多次更新計算終止條件為

考慮到可能會出現某一參數在其取值區間內變化對結構性能影響很小的情況,會使該參數在該取值區間內的每個劃分的后驗信度近似相等,在 1/R附近以較小的幅值波動,即可理解為該參數在該取值區間內可以取任意值.對于上述情況的出現,可終止對該參數的進一步更新計算,可將該取值區間的中心點值作為該參數本次更新計算的結果,并在第k=k+1次更新計算中不再改變該參數的取值.因此,可定義第j個參數的多次更新終止條件

步驟5:在多次更新計算終止后,根據式(24)逐次倒序計算參數取值.若步驟3 中共進行K次更新計算,即

3 數值算例

3.1 P-CS 量化模型構建

3.1.1 一維不確定性變量算例

假定待量化的某不確定性變量描述為X～Nμ ∈[3.5,4.5],σ2∈[0.5,1.5] .該不確定性變量本質上可量化為一族無限個精確概率模型的組合.取該族內8 個概率密度累積函數可繪制成圖5.

圖5 概率密度累積函數族Fig.5 Family of probability density cumulative functions

概率密度累積函數族的邊界分別由四個正態分布概率密度累積函數構成,分別為X～N(4.5,1.5),X～N(4.5,0.5),X～N(3.5,1.5),X～N(3.5,0.5),如圖6 所示.

圖6 概率密度累積函數族的邊界Fig.6 Boundaries for the family of probability density cumulative functions

圖7 等概率水平條件下不確定性變量變化范圍驗證Fig.7 Validation of the variation range of uncertainties under equal probability levels

3.1.2 二維不確定性變量算例

對于給定某一概率水平 α,即可根據式(30)～式(33)分別計算Xα相對于的變化范圍 Ω .考慮到本算例變量間存在相關性,可以用一個最小體積的橢圓模型將在某一概率水平 α下Xα的所有可能取值包絡.選取2 個隨機變量分位數=(4,3) 與=(3,2.6),可得出在該隨機變量分位數處,由等概率原則得到的非概率橢圓凸集模型范圍,如圖8 所示.

圖8 二維 P-CS 模型示意圖Fig.8 Schematic illustration of the 2D P-CS model

為驗證與討論橢圓模型對上述一簇概率隨機變量分位數的包絡性,在本問題的概率密度函數族中選取16 組邊界概率密度函數(μ1,σ1,μ2,σ2分別取最大值與最小值進行排列組合),并分別計算出與隨機變量分位數與等概率水平下的其他16 組概率密度函數對應的分位數.如圖9 所示,橢圓中心點為隨機變量分位數與.對于同等概率水平下的,其16 個點(圖9(a)中每個角點包含4 個位置重復的點)皆位于橢圓模型邊界上;對于同等概率水平下的,其12 個點位于橢圓模型內,4 個角點皆位于橢圓模型邊界上.因此,可以驗證由本文提出的凸集構建方法可構造出包絡等概率水平下一簇隨機變量分位數的最小橢圓模型.3 維及3 維以上多維橢球模型與二維橢圓模型情況完全類似,不再多加討論.

圖9 等概率水平下橢圓模型包絡性驗證Fig.9 Validation of elliptic model envelope under equal probability levels

3.1.3 多源不確定性輸入構建P-CS 模型

圖10 多源不確定性P-CS 模型示意圖Fig.10 Schematic illustration of the multidimensional P-CS model

3.2 懸臂梁結構算例

本算例為一懸臂梁算例,其結構形式如圖11 所示,懸臂梁長度為L=1 m,橫截面慣性矩為I=1×10-4m4,材料彈性模量為E=100 GPa .結構承受集中載荷P1與均布載(荷P2.)集中載荷與均布載荷存在不確定性,P1～,μ1∈[3.5,4.5]kN,σ1=1 kN,P2∈[1.5,2.5] kN/m.假定均布載荷P1與集中載荷P2之間存在相關性,性能測試信息為梁端處轉角位移 β,具體樣本數值統計為表2.由于本文算例為驗證算法合理性,故文中性能測試樣本皆為數值模擬計算產生,非試驗測試值.

表2 懸臂梁結構算例性能測試樣本Table 2 Performance test samples of the example of a cantilevered beam structure10-4/rad

圖11 懸臂梁結構Fig.11 A cantilevered beam structure

根據結構力學原理,可推導出懸臂梁端轉角β與載荷之間的數學表達關系式為

圖12 第1 次更新計算后參數后驗信度分布Fig.12 Posterior credibility distributions after the first updating

根據第一次更新后的信度分布情況,選取下一次更新的參數取值范圍為Pr1∈[1.207,1.311],Pr2∈[1,1.104],θ∈[π/4,π/2] .將上述取值范圍分別做4 個劃分,形成64 個焦元.在第二次更新計算后,參數后驗信度如圖13 所示.

圖13 第2 次更新計算后參數后驗信度分布Fig.13 Posterior credibility distributions after the second updating

根據第二次更新后的信度分布情況,選取下一次更新的參數取值范圍為Pr1∈[1.207,1.259],Pr2∈[1,1.052],θ∈[5π/16,π/2] .將上述取值范圍分別做4 個劃分,形成64 個焦元.在第三次更新計算后,參數后驗信度如圖14 所示.

圖14 第3 次更新計算后參數后驗信度分布Fig.14 Posterior credibility distributions after the third updating

共進行了3 次更新計算后,對于參數 θ,根據式(26),=0,對于參數Pr1與Pr2,根據式(25),ε1=0.03 ＜0.05,滿足了多次更新終止條件.3 次更新計算共形成了192 個焦元,最終信度分布如圖15 所示,最終參數取值根據式(27) 計算為Pr1=1.239,Pr2=1.038,θ=1.205 .為驗證本文焦元選取策略的準確性,對初始參數取值范圍分別做40 個劃分,形成64 000 個焦元,經一次更新計算后根據式(24)計算得到參數取值為Pr1=1.232,Pr2=1.037,θ=0.996 .對于參數Pr1與Pr2,計算誤差分別為0.6%與0.1%.對于參數 θ,由其后驗分布可以看出,當θ ∈[5π/16,π/2]時,其參數真實值落在每一個劃分的信度均等,即參數可以取為該范圍內任何值.由此可見,本文提供的焦元選取策略可以在較少焦元計算量的基礎上獲得較高的計算精度.

圖15 懸臂梁結構算例最終參數后驗信度分布Fig.15 Final posterior credibility distributions of the example of a cantilevered beam structure

在更新過程中,每一個新的樣本點輸入后皆按2.3 節所述的多次更新策略更新參數劃分的后驗信度,進而根據式(27)計算得出當前的P-CS 模型參數.在得到更新后的模型參數后,可根據優化列式(17)計算出當前梁端轉角上限βup與下限βlow,其變化趨勢與均值點處橢圓凸集模型變化如圖16所示.

圖16 懸臂梁結構算例更新過程圖Fig.16 The updating process of the example of a cantilevered beam structure

從圖16 中可以看出,樣本測試值一直處于基于當前P-CS 模型的梁端轉角上限與下限范圍內.前10 個性能測試樣本點變化范圍較小,表示不確定性變量波動范圍較小,P-CS 模型中表征不確定性范圍的橢圓模型逐步變小,同時梁端轉角上限與下限范圍也在縮小.當第11 個性能測試樣本點輸入后,該次采樣點較之前采樣點數值變化很大,這說明不確定性變量從第11 個點開始波動到了另一個范圍,該范圍顯然不在當前P-CS 模型的不確定性范圍內,所以P-CS 中表征不確定性范圍的橢圓模型變大,同時梁端轉角上限與下限范圍也在變大.而后續采樣點較第11 個采樣點變化較小,所以橢圓模型變小同時梁端轉角范圍也在變小.當第21 個性能測試樣本點輸入后,該次采樣點較第11 至第20 采樣點數值變化很大,這說明不確定性變量從第21 個點開始波動到了另一個范圍,而新的變化范圍仍在當前P-CS 模型范圍內,因此P-CS 模型中表征不確定性范圍的橢圓模型繼續變小,同時梁端轉角上限與下限范圍也在繼續變小.

3.3 框架結構算例

本算例為一框架結構,如圖17 所示.該結構由3 根梁構成,兩個梁之間為剛性連接,2 號梁與 3 號梁在3 號節點與4 號節點施加固支約束.1 號梁的長度L1=1.44 m,2 號梁與3 號梁的長度L2=L3=0.96 m.3 根梁的橫截面相同,皆為邊長為 0.1 m 的正方形.結構承受垂直方向的均布載荷P1與水平方向局部載荷P2,P3,其中水平方向的均布載荷P2與P3大小相等,方向相反.不確定性變量包括: 橫向梁(1號梁)材料彈(性模量)E1,E1∈[29,31] GPa(,載荷不)確定性P1～,P23=P2=P3～μ1=μ2∈[36,44] kN,σ1=σ2=2 kN,載荷不確定性變量P1與P23相關.2 號梁與3 號梁的材料彈性模量分別為E2=E3=30 GPa .性能測試信息為角點1 處轉角位移 β,具體樣本數值統計為表3.由于本文算例為驗證算法合理性,故文中性能測試樣本皆為數值模擬計算產生,非試驗測試值.可基于有限元方法建立2 號角點轉角位移 β 與不確定變量之間的數學關系,即β=g(X) .

圖17 框架結構Fig.17 A frame structure

表3 框架結構算例性能測試樣本(10-4/rad)Table 3 Performance test samples of the example of a frame structure (10-4/rad)

圖18 框架結構算例最終參數后驗信度分布Fig.18 Final posterior credibility distributions of the example of a frame structure

圖18 框架結構算例最終參數后驗信度分布(續)Fig.18 Final posterior credibility distributions of the example of a frame structure (continued)

隨著性能測試信息 β*輸入,每一個新的樣本點輸入后皆按3.3 節所述的多次更新策略更新參數劃分的后驗信度,進而根據式(27) 計算得出當前的P-CS 模型參數.在得到更新后的模型參數后,可根據優化列式(17)計算出當前轉角位移的上限 βup與下限βlow,其變化趨勢與均值點處多橢球凸集模型變化如圖19 所示.

圖19 框架結構算例更新過程圖Fig.19 The updating process of the example of a frame structure

從圖19 中可以看出,樣本測試值一直處于基于當前P-CS 模型的轉角位移上限與下限范圍內.前10 個性能測試樣本點變化范圍較小,表示不確定性變量波動范圍較小,P-CS 模型中表征不確定性范圍的橢圓模型逐步變小,同時轉角位移上限與下限范圍也在縮小.當第11 個性能測試樣本點輸入后,該次采樣點較之前采樣點數值變化很大,而后續采樣點較第11 個采樣點變化較小,這說明不確定性變量從第11 個點開始波動到了另一個范圍,因此P-CS模型中表征不確定性范圍的橢圓模型變大,同時轉角位移上限與下限范圍也在變大.

4 結論

本文提出的P-CS 不確定性量化模型可以將概率隨機變量、非概率凸集變量與非精確概率隨機變量統一在一個模型框架下進行表征,從而實現對多源不確定性變量的統一模型量化.進而基于P-CS 模型,提出了一種性能測試信息驅動P-CS 模型參數更新的算法.

P-CS 不確定性量化模型由隨機變量概率模型部分與多橢球非概率模型兩部分組成,既具有概率特性也具有非概率特性.在給定某一概率水平的前提下,根據等概率原則可推導計算出在隨機變量分位數,與其具有相同概率特性的不確定性變量取值范圍.在數值算例1 中,本文詳細闡述了一維、二維與三維P-CS 模型構建過程與其概率、非概率特性.

基于貝葉斯原理,提出了一種性能測試信息驅動P-CS 模型參數更新算法.采用參數信度概念取代傳統貝葉斯理論框架內的條件概率進行更新計算,從而實現根據性能測試信息更新P-CS 模型參數的目的.數值算例中,本文分別提供了承受不確定性均布載荷與集中載荷的懸臂梁結構算例以及考慮材料屬性與載荷不確定性耦合的框架結構算例,上述兩個力學算例驗證了本文提出的更新算法的適用性.

本文提出的P-CS 不確定性量化模型可解決多源、時變不確定性變量的量化與更新問題,并且該更新算法不需要在結構服役過程中對不確定性變量進行直接測量采樣.該量化模型與其更新方法有著很強的實際工程應用價值,有望應用于考慮時變不確定性的結構形面、振動主動控制問題中.