項(xiàng)涌涌 潘柏松 羅路平 謝少軍

浙江工業(yè)大學(xué)機(jī)械工程學(xué)院，杭州，310014

0 引言

由于存在物理試驗(yàn)成本和客觀條件等限制，計(jì)算機(jī)仿真建模和分析已逐漸應(yīng)用于設(shè)計(jì)及評估航空航天、機(jī)械等領(lǐng)域的復(fù)雜系統(tǒng)。高精度的數(shù)值模型是開展有效仿真分析的基礎(chǔ)， 滿足使用要求的數(shù)值模型才可用于仿真分析。工程實(shí)際中，普遍存在概率與區(qū)間混合不確定性情況[1-2]，有些模型參數(shù)可獲得充足樣本數(shù)據(jù)來構(gòu)建精確的概率分布；而有些模型參數(shù)試驗(yàn)成本高及難度大，難 以獲得足夠樣本數(shù)據(jù)來精確描述，只能估計(jì)其變化區(qū)間。雖然多數(shù)情況下估計(jì)區(qū)間的不確定性水平僅略大于真實(shí)區(qū)間的不確定性水平，但概率與區(qū)間混合時(shí)會(huì)導(dǎo)致模型輸出響應(yīng)的不確定性水平過高，從而使得模型可信度大幅降低。

目前模型確認(rèn)方法廣泛應(yīng)用于評價(jià)和提高建?？尚哦龋钤缬擅绹茉床繎?yīng)用于戰(zhàn)略武器管理的可靠性評估和決策[3-4]。由于工程領(lǐng)域?qū)Ψ抡娣治隹尚哦纫蟮奶岣?，模型確認(rèn)方法正大量應(yīng)用于計(jì)算力學(xué)、汽車等工程領(lǐng)域[5-6]。隨著模型確認(rèn)技術(shù)的應(yīng)用推廣，人們逐漸認(rèn)識(shí)到僅僅實(shí)施評估模型精確度的模型確認(rèn)過程并不能滿足工程實(shí)際需求，需對未達(dá)到要求的模型進(jìn)行模型修正[7]。

模型確認(rèn)方法可分為基于貝葉斯理論的貝葉斯方法[8]和基于統(tǒng)計(jì)及推斷方法的非貝葉斯方法[9-10]兩大類。KENNEDY等[8]提出了基于貝葉斯理論的模型修正方法；LIU等[11]針對熱傳導(dǎo)問題的模型修正，采用貝葉斯方法來提高模型精度；HASSELMAN等[12]針對汽車碰撞模型提出了一種基于模型響應(yīng)值均方差結(jié)果的修正方法。

近年來，為確保仿真工作具有實(shí)際應(yīng)用價(jià)值和意義，對模型確認(rèn)問題的研究已受到國內(nèi)學(xué)者們的重視。張保強(qiáng)等[13]將貝葉斯方法用于圣地亞熱傳導(dǎo)問題的模型修正中，并建立相應(yīng)的貝葉斯框架；何成等[14]分析了基于徑向基模型的不確定性結(jié)構(gòu)動(dòng)力學(xué)模型非概率型區(qū)間修正方法；肖釗等[15]針對試驗(yàn)數(shù)據(jù)較少的情況，利用區(qū)間方法量化不確定性，并提出了對應(yīng)模型確認(rèn)方法；方圣恩等[16]提出了針對結(jié)構(gòu)參數(shù)不確定性的隨機(jī)模型優(yōu)化方法；姜東等[17]基于區(qū)間分析提出了不確定性結(jié)構(gòu)動(dòng)力學(xué)模型修正方法。綜上可知，人們已開展了一些模型確認(rèn)問題的工作，并取得了一定成果，但針對概率與區(qū)間混合不確定性情況的研究并不多，且總體上就模型確認(rèn)問題并未達(dá)成統(tǒng)一共識(shí)。

由于存在不確定性傳播的問題，在概率與區(qū)間混合不確定性情況下，模型輸出響應(yīng)的數(shù)值結(jié)果為非精確概率。本文通過數(shù)值結(jié)果非精確概率累積分布函數(shù)(cumulative distribution function,CDF)與試驗(yàn)結(jié)果經(jīng)驗(yàn)累積分布函數(shù)(empirical cumulative distribution function，ECDF)之間的面積極大值求解模型確認(rèn)準(zhǔn)則參數(shù)，并對模型確認(rèn)結(jié)果進(jìn)行評估，對于不滿足確認(rèn)評估標(biāo)準(zhǔn)的模型，以區(qū)間變量表征的不確定性參數(shù)為對象實(shí)施修正，最后將模型確認(rèn)方法應(yīng)用于圣地亞熱傳導(dǎo)問題的算例中，并驗(yàn)證了所提方法的有效性。

1 參數(shù)不確定性模型確認(rèn)過程

實(shí)際工程模型中普遍存在隨機(jī)不確定性參數(shù)、認(rèn)知不確定性參數(shù)、確定性參數(shù)混合情況，本文針對區(qū)間變量表征認(rèn)知不確定性參數(shù)，概率變量表征隨機(jī)不確定性參數(shù)及確定性參數(shù)的混合情況進(jìn)行了研究分析，得到模型函數(shù)可表示為

Y=g(X,Z,M)

(1)

式中，X為區(qū)間變量矢量；Z為概率變量矢量；M為確定性變量矢量；Y為模型輸出響應(yīng)(隨機(jī)變量)。

受不確定性混合影響，模型中只要存在參數(shù)不確定性水平略大于真實(shí)水平的情況，就會(huì)使模型輸出響應(yīng)的不確定性水平過高，從而需進(jìn)行模型確認(rèn)以保證響應(yīng)準(zhǔn)確度和建模可信度。

模型確認(rèn)的目的是量化數(shù)值模型與真實(shí)物理系統(tǒng)的差異程度，并基于差異程度確定模型預(yù)測精度。在模型確認(rèn)過程中，對于未達(dá)到模型確認(rèn)要求的數(shù)值模型，需通過模型修正以減小模型誤差，從而使數(shù)值模型滿足預(yù)期要求。在參數(shù)不確定性模型中對模型參數(shù)進(jìn)行修正，概率與區(qū)間混合不確定性情況的修正參數(shù)為區(qū)間變量表征的不確定性參數(shù)。

真實(shí)物理系統(tǒng)的量化表征是以輸出響應(yīng)的試驗(yàn)數(shù)據(jù)為基礎(chǔ)，響應(yīng)試驗(yàn)誤差會(huì)直接影響確認(rèn)結(jié)果的可信度，因此本文將響應(yīng)試驗(yàn)誤差視為數(shù)值模型的隨機(jī)不確定性參數(shù)來開展模型確認(rèn)過程。模型確認(rèn)實(shí)施步驟包括模型確認(rèn)試驗(yàn)，不確定量化，應(yīng)用確認(rèn)準(zhǔn)則進(jìn)行模型確認(rèn)評估，對不滿足要求的模型進(jìn)行參數(shù)修正，以及模型應(yīng)用等，見圖1。

圖1 參數(shù)不確定性模型確認(rèn)過程Fig.1 Validation process of parameteruncertainty model

模型確認(rèn)試驗(yàn)和確認(rèn)準(zhǔn)則在模型確認(rèn)過程中至關(guān)重要，模型確認(rèn)試驗(yàn)的設(shè)計(jì)、執(zhí)行、分析是以定量確定數(shù)值模型及其載體(仿真軟件)在多大程度上能夠代表實(shí)際系統(tǒng)為目的[6]。不同工程實(shí)際問題應(yīng)根據(jù)確認(rèn)的需求對系統(tǒng)、子系統(tǒng)、部件和單元的各個(gè)層面分別進(jìn)行試驗(yàn)測試。模型確認(rèn)準(zhǔn)則主要可分為圖像對比準(zhǔn)則、頻率準(zhǔn)則、經(jīng)典假設(shè)檢驗(yàn)法、貝葉斯因子、面積度量準(zhǔn)則等[18]。在模型確認(rèn)過程中，需根據(jù)實(shí)際問題應(yīng)用合理的確認(rèn)準(zhǔn)則以獲得可信的確認(rèn)結(jié)果。

2 概率與區(qū)間混合不確定性下參數(shù)不確定性模型確認(rèn)理論

2.1 區(qū)間變量

區(qū)間變量的數(shù)學(xué)定義為[19]

XI=[XL,XU]=[x|XL≤x≤XU,x∈R]

(2)

式中，x為實(shí)數(shù)集變量；XL、XU分別為區(qū)間下限和區(qū)間上限。

區(qū)間變量也可利用區(qū)間中點(diǎn)和區(qū)間半徑來表示，即

XI=[XC-XR,XC+XR]

(3)

式中，XC為區(qū)間中點(diǎn)；XR為區(qū)間半徑。

應(yīng)用區(qū)間變量來量化不確定性，區(qū)間可靠度是評價(jià)量化結(jié)果的指標(biāo)之一，區(qū)間可靠度H可用區(qū)間置信系數(shù)來定義[20]，即

H=PXI{XL≤x≤XU}∝XR

(4)

其中，PXI為區(qū)間變量XI包含參數(shù)所有值的概率。H與區(qū)間半徑XR成正比(即區(qū)間半徑XR越大，區(qū)間變量的置信水平越高，區(qū)間估計(jì)越可靠)。

2.2 模型輸出響應(yīng)數(shù)值結(jié)果

不確定性傳播中概率與區(qū)間混合不確定性情況下模型的輸出響應(yīng)數(shù)值結(jié)果為非精確概率，本文采用概率盒(probability box , P-box)來表征[21]。如圖2所示，概率盒將上下兩條累積分布函數(shù)(CDF)作為邊界，數(shù)學(xué)定義為[22]

(5)

圖2 概率盒示意圖Fig.2 Schematic diagram of probability box

2.3 模型確認(rèn)與評估

面積度量準(zhǔn)則[23]是一種有效的模型定量確認(rèn)準(zhǔn)則，以仿真值CDF與試驗(yàn)值ECDF的面積差值作為準(zhǔn)則值。為確保概率與區(qū)間混合不確定性情況下模型確認(rèn)的有效性，本文考慮確認(rèn)準(zhǔn)則值的量綱一性，以面積度量準(zhǔn)則為基礎(chǔ)，根據(jù)數(shù)值結(jié)果概率盒上下邊界與試驗(yàn)結(jié)果ECDF之間的面積極大值定義了模型確認(rèn)準(zhǔn)則，具體表達(dá)式如下：

(6)

式中，VP為模型確認(rèn)準(zhǔn)則參數(shù)；max{·}為求解最大集合元素函數(shù);B為輸出響應(yīng)試驗(yàn)數(shù)據(jù)的中點(diǎn);Se(y)為響應(yīng)試驗(yàn)結(jié)果ECDF。

圖3 概率盒邊界與試驗(yàn)結(jié)果ECDF的面積差值圖Fig.3 Different area between probability boxboundary and the ECDF of test results

模型確認(rèn)準(zhǔn)則可量化試驗(yàn)結(jié)果與仿真結(jié)果的差異，求解得到的VP值越大，表明數(shù)值模型和真實(shí)物理對象的吻合程度越低，為此提出模型確認(rèn)評估標(biāo)準(zhǔn)，即

VP≤δ

(7)

其中，δ為確認(rèn)標(biāo)準(zhǔn)值，根據(jù)模型實(shí)際應(yīng)用條件取某一常數(shù)值。當(dāng)數(shù)值模型有較高要求時(shí)，需設(shè)定較小的δ值，反之則可設(shè)定較大的δ值。

2.4 模型修正

當(dāng)確認(rèn)準(zhǔn)則值大于標(biāo)準(zhǔn)值δ時(shí)，需對模型參數(shù)進(jìn)行修正。模型參數(shù)修正實(shí)質(zhì)為優(yōu)化待修正參數(shù)，本文以模型確認(rèn)評估標(biāo)準(zhǔn)和模型參數(shù)初始值 為約束條件，以待修正區(qū)間不確定性參數(shù)總置信水平最大為目標(biāo)實(shí)施參數(shù)修正，得到優(yōu)化模型如下：

(8)

式中，H(X)為待修正區(qū)間不確定性參數(shù)總置信水平函數(shù)；Θ為模型參數(shù)初始值域。

3 算例(圣地亞熱傳導(dǎo)問題)

圖4為一維瞬態(tài)圣地亞熱傳導(dǎo)模型示意圖[25]，數(shù)學(xué)模型可表示為

(9)

x=(q,L,t,γ)β=(k,cV)

式中，θi為初始溫度，這里為25 ℃；q為熱流密度；L為板厚；k為材料熱導(dǎo)率；cV為體積熱容；γ為表面距離；t為時(shí)間；x為受控輸入?yún)?shù)矢量；β為待修正參數(shù)矢量。

圖4 熱傳導(dǎo)模型示意圖Fig.4 Schematic diagram of thermal challengeproblem model

圣地亞熱傳導(dǎo)模型的預(yù)測目標(biāo)為受控輸入?yún)?shù)矢量x=(3 500 W/m2,0.019 m,1 000 s,0 m)時(shí)，溫度大于900 ℃的概率小于0.01，即

P{θ(x,β)>900 ℃}<0.01

(10)

圣地亞熱傳導(dǎo)問題中試驗(yàn)誤差無法避免，本文通過考慮模型輸出響應(yīng)溫度的試驗(yàn)誤差來開展模型確認(rèn)，則熱傳導(dǎo)問題數(shù)值模型可表示為

θ(x,β,εe)=θ(x,β)+εe

(11)

其中，受控輸入?yún)?shù)x為確定性變量矢量；待修正參數(shù)β為認(rèn)知不確定性變量矢量，應(yīng)用區(qū)間變量表征；響應(yīng)試驗(yàn)誤差εe為隨機(jī)不確定性變量，應(yīng)用概率分布表征。

模型確認(rèn)熱傳導(dǎo)問題試驗(yàn)分為材料特性(material characterization，MC)試驗(yàn)，整體確認(rèn)(ensemble validation，EN)試驗(yàn)，認(rèn)證(accreditation validation，AC)試驗(yàn)，所有試驗(yàn)根據(jù)數(shù)據(jù)量的大小分為低水平、中水平和高水平等級。由于高水平試驗(yàn)數(shù)據(jù)量與其他兩個(gè)水平的數(shù)據(jù)量差距較小，為確保模型確認(rèn)的有效性，故本文選用高水平試驗(yàn)數(shù)據(jù)進(jìn)行熱傳導(dǎo)問題模型確認(rèn)，試驗(yàn)數(shù)據(jù)詳見文獻(xiàn)[25]。

溫度對待修正參數(shù)有較大影響，本文考慮待修正參數(shù)的溫度依賴性求解函數(shù)方程，分析試驗(yàn)數(shù)據(jù)可知，體積熱容cV為獨(dú)立變量，則體積熱容區(qū)間變量的表達(dá)式如下：

(12)

式中，N為試驗(yàn)數(shù)據(jù)數(shù)量。

材料熱導(dǎo)率k對溫度存在依賴性，建立關(guān)于溫度θ的一維線性回歸模型，表達(dá)式如下：

k=a+bθ+ω

(13)

根據(jù)MC試驗(yàn)數(shù)據(jù)求得參數(shù)a為0.051 2，b為2.25×10-5，誤差項(xiàng)ω～N(0,0.005 32)，根據(jù)正態(tài)分布求解誤差項(xiàng)區(qū)間變量形式(即ωI=[μ-3σ,μ+3σ])，得到含區(qū)間變量的方程如下：

k=0.051 2+2.25×10-5θ+[μ-3σ，μ+3σ]

(14)

式中，μ為正態(tài)分布均值；σ為正態(tài)分布標(biāo)準(zhǔn)差。

根據(jù)文獻(xiàn)[7]可得試驗(yàn)誤差εe的分布為εe～N(0,5.252)。

熱傳導(dǎo)問題模型確認(rèn)中可用函數(shù)表示式來表示預(yù)測溫度與各參數(shù)之間的關(guān)系：θ=g(x,β,εe)，將區(qū)間變量(式(12)、式(14))及各概率分布代入上述函數(shù)表達(dá)式，可得熱傳導(dǎo)問題模型確認(rèn)方程如下：

(15)

熱傳導(dǎo)問題模型確認(rèn)方程是以θ為變量的非線性方程，可進(jìn)行優(yōu)化求解。本算例的模型確認(rèn)標(biāo)準(zhǔn)值設(shè)定為δ=0.3，應(yīng)用所提確認(rèn)準(zhǔn)則進(jìn)行模型確認(rèn)，熱傳導(dǎo)問題模型確認(rèn)評估的表達(dá)式為

(16)

圣地亞熱傳導(dǎo)問題設(shè)計(jì)進(jìn)行了不同的EN試驗(yàn)和AC試驗(yàn)，EN試驗(yàn)的高水平數(shù)據(jù)為在0～1 000 s范圍內(nèi)，以間隔100 s選取得到10個(gè)時(shí)間點(diǎn)的4組試驗(yàn)數(shù)據(jù)，每組試驗(yàn)重復(fù)4次；AC試驗(yàn)高水平數(shù)據(jù)為在0～1 000 s范圍內(nèi)，以間隔50 s選取得到20個(gè)時(shí)間點(diǎn)的3組試驗(yàn)數(shù)據(jù)，每組試驗(yàn)重復(fù)2次。表1和表2分別列出了熱傳導(dǎo)問題各EN試驗(yàn)參數(shù)和AC試驗(yàn)參數(shù)。

表1 整體確認(rèn)試驗(yàn)參數(shù)

表2 認(rèn)證試驗(yàn)參數(shù)

本文綜合應(yīng)用EN試驗(yàn)和AC試驗(yàn)的高水平數(shù)據(jù)進(jìn)行模型確認(rèn)，根據(jù)表1和表2中的具體參數(shù)值，將不確定性參數(shù)表征值代入式(15)，運(yùn)用雙層采樣法得到θ的概率盒模型，最后根據(jù)式(16)分別求得各時(shí)刻下的VP(t)。

圖5a為EN試驗(yàn)條件下未修正確認(rèn)準(zhǔn)則的參數(shù)曲線圖，可以看出，所有試驗(yàn)序號各時(shí)刻下的確認(rèn)準(zhǔn)則參數(shù)值均小于0.3，滿足確認(rèn)評估要求。圖5b為AC試驗(yàn)條件下未修正確認(rèn)準(zhǔn)則的參數(shù)曲線圖，可以看出，在0～1 000 s范圍內(nèi)AC1試驗(yàn)的確認(rèn)準(zhǔn)則參數(shù)值滿足要求；而AC2試驗(yàn)的確認(rèn)準(zhǔn)則參數(shù)值有部分不滿足要求，且t=200 s為確認(rèn)準(zhǔn)則參數(shù)曲線極值點(diǎn)(即VP(t=200 s)為0.36>δ)；AC3試驗(yàn)中t=450 s為曲線極值點(diǎn)(即VP(t=450 s)為0.7>δ)，且大部分AC3試驗(yàn)確認(rèn)準(zhǔn)則的參數(shù)值不滿足要求。

(a)整體確認(rèn)試驗(yàn)

(b)認(rèn)證試驗(yàn)圖5 未修正模型確認(rèn)準(zhǔn)則參數(shù)值Fig.5 Values of uncorrected model validationcriteria parameters

在熱傳導(dǎo)問題中，AC2試驗(yàn)和AC3試驗(yàn)的確認(rèn)準(zhǔn)則參數(shù)值不滿足確認(rèn)評估要求，應(yīng)進(jìn)行模型參數(shù)修正。由于AC3試驗(yàn)條件下確認(rèn)準(zhǔn)則參數(shù)值不滿足要求的比例遠(yuǎn)大于AC2試驗(yàn)條件下的比例，且AC3試驗(yàn)條件下準(zhǔn)則參數(shù)值曲線的極值點(diǎn)VP(t=450 s)=0.7大于AC2試驗(yàn)條件下準(zhǔn)則參數(shù)值曲線極值點(diǎn)VP(t=200 s)=0.36，故依據(jù)AC3試驗(yàn)進(jìn)行模型參數(shù)修正。AC3試驗(yàn)的準(zhǔn)則參數(shù)值曲線為開口向下的拋物線，在t=450 s時(shí)達(dá)到最高點(diǎn)，所以修正后滿足VP(t=450 s)≤δ的條件，就可使所有曲線均滿足要求?；贏C3試驗(yàn)條件下t=450 s的確認(rèn)結(jié)果修正模型參數(shù)，得到優(yōu)化模型如下：

(17)

目標(biāo)函數(shù)經(jīng)簡化后，運(yùn)用遺傳算法求解優(yōu)化模型，表3為模型參數(shù)修正結(jié)果，表4為修正前后參數(shù)值對比。

基于修正后參數(shù)值進(jìn)行模型確認(rèn)，圖6a和圖6b分別為修正后EN試驗(yàn)和AC試驗(yàn)條件下的模型確認(rèn)準(zhǔn)則參數(shù)值曲線。由圖6可知，修正后全部確認(rèn)試驗(yàn)所有時(shí)刻點(diǎn)的模型確認(rèn)準(zhǔn)則參數(shù)值均小于0.3，數(shù)值模型滿足使用要求。

表3 模型參數(shù)修正結(jié)果

表4 參數(shù)修正前后結(jié)果對比

(a)整體確認(rèn)試驗(yàn)

(b)認(rèn)證試驗(yàn)圖6 已修正模型確認(rèn)準(zhǔn)則參數(shù)值Fig.6 Values of model validation criteria parametersafter correction

將修正后的模型參數(shù)代入式(10)求解，采用蒙特卡羅法抽樣(抽樣次數(shù)為1×106次)可得到失效概率。圖7為修正后溫度θ的累積分布函數(shù)曲線，預(yù)測結(jié)果見表5。由表5可知，模型未修正時(shí)，溫度超過900 ℃的概率為0.078 1，模型修正后，溫度超過900 ℃的概率為0.016 2，預(yù)測精度得到明顯提高。

圖7 修正后溫度的累積分布曲線Fig.7 The CDF curve of temperature after correction

模型預(yù)測未修正修正后P{θ(x,β)>900 ℃}0.078 10.016 2

熱傳導(dǎo)問題中失效概率的目標(biāo)要求不能大于0.01，采用本文方法得到的模型預(yù)測概率結(jié)果也不滿足要求，這與同類文獻(xiàn)[11,26-28]的結(jié)論一致。采用不同方法的模型修正對比見表6，其中文獻(xiàn)[11,27]基于貝葉斯原理實(shí)施模型修正。

表6 不同方法模型修正對比

由表6可知，三種方法的修正后模型預(yù)測概率十分相近，但文獻(xiàn)[11,27]所用試驗(yàn)數(shù)據(jù)量是本文方法所用試驗(yàn)數(shù)據(jù)量的3倍以上。基于貝葉斯的模型修正方法以增加試驗(yàn)數(shù)據(jù)量為基礎(chǔ)，且模型參數(shù)先驗(yàn)分布的確定存在較大主觀性，因此，在工程實(shí)際試驗(yàn)難度較大、試驗(yàn)成本較高的情況下，本文方法具有明顯的優(yōu)勢。

4 結(jié)論

(1)在概率與區(qū)間混合不確定性情況下，本文的模型確認(rèn)準(zhǔn)則依據(jù)數(shù)值結(jié)果概率盒邊界與試驗(yàn)結(jié)果經(jīng)驗(yàn)累積分布函數(shù)之間的面積極大值進(jìn)行定義，可獲得有效可信的模型確認(rèn)結(jié)果值，為后續(xù)確認(rèn)評估及修正奠定了基礎(chǔ)。

(2)以待修正區(qū)間不確定性參數(shù)總置信水平最大為修正目標(biāo)，以模型確認(rèn)評估標(biāo)準(zhǔn)和模型參數(shù)初始值為約束條件，目標(biāo)函數(shù)簡化后應(yīng)用遺傳算法修正模型參數(shù)，獲得的修正結(jié)果具有可信性。

(3)將本文方法應(yīng)用于圣地亞熱傳導(dǎo)問題，經(jīng)模型確認(rèn)后熱傳導(dǎo)模型的預(yù)測精度有明顯提高，模型可信度和響應(yīng)準(zhǔn)確度有較大提升，獲得了與同類文獻(xiàn)相同的預(yù)測結(jié)果。修正后模型預(yù)測概率與基于貝葉斯的模型修正方法得到的修正后預(yù)測概率的結(jié)果十分相近，但本文方法所用試驗(yàn)數(shù)據(jù)量有明顯減小。在試驗(yàn)難度較大、試驗(yàn)成本較高的情況下，本文方法比基于貝葉斯的模型修正方法更具優(yōu)勢，具有較好的應(yīng)用前景。