基于代理模型的多體動力學(xué)模型修正技術(shù)

方劍光，高云凱，徐成民，張玉婷

(1．同濟大學(xué)汽車學(xué)院，上海 201804;2．悉尼大學(xué)航空、機械和機電一體化工程學(xué)院，悉尼，NSW 2006)

前言

在汽車耐久性虛擬開發(fā)流程中，獲得車身或其他部件正確的邊界載荷，是準(zhǔn)確預(yù)測其疲勞壽命的先決條件。通常，這些邊界載荷從多體動力學(xué)仿真中提取［1－4］，因此建立高精度能反映實際動力學(xué)特性的模型對于后續(xù)的疲勞壽命預(yù)測尤為重要。

然而，建立多體模型過程中經(jīng)常遇到的困難是模型中各個部件動力學(xué)特性表征參數(shù)的模擬，例如懸架系統(tǒng)阻尼參數(shù)等。為了獲取這些參數(shù)，就要求進(jìn)行室內(nèi)臺架或室外道路試驗。而實際開發(fā)過程中很有可能會受到實驗成本和開發(fā)周期等限制，這些參數(shù)無法獲取，只能粗略估計，這就將造成模型的動力學(xué)性能較差。在這種情況下，對這些估計的參數(shù)進(jìn)行簡單的人工調(diào)整，往往很難達(dá)到理想的效果:首先，對每個參數(shù)進(jìn)行調(diào)整無疑會降低效率;而且，由于非線性因素的存在，盲目地調(diào)整并不能得到正確的結(jié)果;此外，隨著參數(shù)數(shù)量的增加，參數(shù)之間的交互作用也將增加人工調(diào)整的難度。

因此，本文中引入?yún)?shù)識別方法以實現(xiàn)多體模型的修正。文獻(xiàn)［5］中將響應(yīng)面模型用于5自由度動力學(xué)系統(tǒng)的參數(shù)識別，對于線性系統(tǒng)取得了成功，而對于非線性系統(tǒng)還有待深入。文獻(xiàn)［6］中將參數(shù)識別定義為逆問題，在時域內(nèi)通過全局和子結(jié)構(gòu)法識別了彈簧和阻尼器的非線性參數(shù)。文獻(xiàn)［7］中將正交分解與混合適應(yīng)性群體智能算法相結(jié)合，根據(jù)測量的加速度信號在時域內(nèi)識別結(jié)構(gòu)的剛度和阻尼參數(shù)。

本文中將代理模型技術(shù)、優(yōu)化技術(shù)與參數(shù)識別技術(shù)有機結(jié)合，以某試車場路面為例，高效地反算汽車多體動力學(xué)表征參數(shù)，實現(xiàn)多體模型的修正，旨在最大限度地降低汽車的開發(fā)成本和縮短開發(fā)周期。

1 模型修正方法

模型修正過程首先要求建立一個目標(biāo)函數(shù)，即模型誤差評價函數(shù)，它反映了一個物理系統(tǒng)的數(shù)值模型和試驗?zāi)Ｐ椭g在輸出響應(yīng)上的差異。一個最小二乘問題的目標(biāo)函數(shù)可以定義為

式中:N代表實驗曲線上取樣點的個數(shù)，見圖1;x為多體模型中待修正參數(shù)向量;Wi為第i個取樣點的權(quán)系數(shù);Ti為第i個取樣點處的實驗測量值;fi(x)為第i個取樣點處近似模型的響應(yīng)值;si為第i個取樣點的殘差比例系數(shù);ei(x)為第i個取樣點的殘差值。本文中通過仿真和實驗的時域曲線之間的差異來評價模型誤差，N=512，Wi和si等于1。定義模型誤差函數(shù)之后，模型修正問題就轉(zhuǎn)化為通過調(diào)整模型參數(shù)使誤差最小化的優(yōu)化問題。

本文中將結(jié)合實驗設(shè)計、代理模型技術(shù)和二次規(guī)劃算法，實現(xiàn)用于耐久性分析的動力學(xué)模型的修正，流程如圖2所示。

2 多體動力學(xué)建模及參數(shù)化

用MSC．ADAMS軟件建立車輛的多體動力學(xué)模型，如圖3所示。其中，前懸架為麥弗遜懸架，由減振器總成、擺臂、橫向穩(wěn)定桿和轉(zhuǎn)向系統(tǒng)組成，后懸架為非獨立懸架，由鋼板彈簧、減振器、限位塊和車橋組成。車身采用質(zhì)量點建模，與懸架通過襯套連接。

在多體動力學(xué)模型中，非線性因素的存在無疑大大提高了模型的復(fù)雜程度，使模型修正更為棘手。考慮到減振器總成參數(shù)對整車動力學(xué)性能的重要性，本文中將懸架系統(tǒng)前減振器的阻尼特性和緩沖塊剛度特性作為調(diào)整參數(shù)，考慮到左右對稱性，以底盤左前懸架處與車身之間的相對位移信號來衡量多體模型的精確性，圖4給出了已有同類車型減振器的相應(yīng)特性。

由阻尼特性曲線可以看出，載荷與速度是多段線性關(guān)系:緩沖塊剛度曲線可分為3段，當(dāng)位移較小(即緩沖塊起作用)時載荷為0，之后載荷線性增加，當(dāng)緩沖塊受擠壓到一定程度時載荷成指數(shù)增加。因此，阻尼和剛度特性曲線可表述為

3 篩選參數(shù)

從眾多影響模型精度的因素中篩選出最主要的幾個因素，它們對參數(shù)識別的效率非常重要。在工業(yè)界，正交實驗設(shè)計方法［8］已被廣泛用作從眾多影響因素中快速篩選出幾個重要因素的一種快速有效的方法。采用三水平的L27正交表分析了表1中11個因素中對位移響應(yīng)的影響。圖5為各因素不同水平下誤差評價函數(shù)的平均值。挑出影響最顯著的5個變量 (v3、k1、k2、k3和 k5)作為模型修正參數(shù)，優(yōu)化問題的變量為x=［v3k1k2k3k5］T。

表1 正交試驗設(shè)計參數(shù)信息

4 代理模型

在工程優(yōu)化中，直接將仿真模型與優(yōu)化算法耦合會導(dǎo)致優(yōu)化過程的低效，因為優(yōu)化迭代中仿真分析需要消耗大量的計算成本。此外，優(yōu)化過程中多體動力學(xué)仿真中存在積分誤差等帶來的模型不穩(wěn)定性，使得單次仿真失效而導(dǎo)致整個優(yōu)化過程的終止。根據(jù)以上原因，采用基于代理模型的優(yōu)化策略。代理模型可視為代替仿真模型，近似表達(dá)設(shè)計變量和響應(yīng)關(guān)系的傳遞函數(shù)。將其代替仿真模型進(jìn)行優(yōu)化，可有效提高優(yōu)化效率。

(1)響應(yīng)面模型 響應(yīng)面模型是最典型的代理模型，最早在1951年提出［9］。在工程上廣泛應(yīng)用的二階響應(yīng)面模型的數(shù)學(xué)表達(dá)式為

式中:y^為響應(yīng)的估計值;xi為第i個輸入變量;a0、ai、aii和aij為待定系數(shù)，由最小二乘法求得。

(2)徑向基模型 徑向基函數(shù)(radial basis function，RBF)是在1971年首次提出的一種用于離散多元數(shù)據(jù)的擬合方法［10］，其表達(dá)式為

式中:μ為多項式模型;‖x－xi‖是點x與xi的歐氏距離;φ是基函數(shù);λi是基函數(shù)在取樣點xi處的未知加權(quán)系數(shù)。因此，徑向基函數(shù)模型實際上是由基于歐氏距離的徑向?qū)ΨQ基函數(shù)和多項式函數(shù)的線性組合。

關(guān)于基函數(shù)φ的形式，高斯函數(shù)和多重二次函數(shù)具有較好的整體特性［11］，其基函數(shù)形式如表2所示。c為自由形狀參數(shù)。

表2 基函數(shù)

(3)Kriging模型 Kriging模型由全局近似模型與局部偏差結(jié)合而成［12］:式中:y(x)為待擬合的響應(yīng)函數(shù);f(x)通常是多項式響應(yīng)面近似模型，一般被看作常數(shù)，提供設(shè)計空間的全局近似模型;Z(x)是期望為零，方差為σ2和協(xié)方差為非零的隨機偏差，提供局部偏差，使得Kriging模型可在ns個樣本點做插值。Z(x)的協(xié)方差矩陣為

Cov［Z(xi)，Z(xj)］ = σ2R［R(xi，xj)］ (5)式中:R為相關(guān)矩陣，是對角線上全為1的(ns×ns)對稱矩陣;R(xi，xj)為ns中任意兩個樣本點xi和xj的相關(guān)函數(shù)，一般采用高斯相關(guān)函數(shù):

式中:θk為用于擬合近似模型的未知相關(guān)系數(shù);xi

k和xj

k分別為樣本點xi，xj第k個元素。

響應(yīng)y在未知位置x點的估計值y^為

其中:

式中:y為ns個樣本點的響應(yīng)列向量;f為長度為ns的單位列向量;rT(x)為ns個樣本點 x1，…，xn{

}

s與未知位置x所組成的相關(guān)向量:

其中:

詳細(xì)過程參考文獻(xiàn)［13］。

代理模型僅僅是對復(fù)雜仿真分析模型的簡化和近似，因此，在求出代理模型的待估參數(shù)后，須對代理模型進(jìn)行統(tǒng)計檢驗，評估代理模型對真實響應(yīng)的逼近程度。在精度不高的代理模型基礎(chǔ)上獲得的最優(yōu)結(jié)果不但不能指導(dǎo)設(shè)計，反而會誤導(dǎo)設(shè)計，因此，對代理模型的誤差進(jìn)行評估至關(guān)重要。常用的誤差估計準(zhǔn)則有:決定系數(shù) R2、相對平均絕對誤差(RAAE)和相對最大絕對誤差(RMAE)。

式中:yi、y^和y－分別為第i個檢驗點響應(yīng)量的實測值、響應(yīng)量的預(yù)測值和響應(yīng)量實測值的平均值;q為檢驗點的個數(shù)。

從上述的定義式可以看出，R2的取值在［0，1］區(qū)間內(nèi)，R2越接近于1，表明回歸擬合的效果越好。此外，R2和RAAE為全局指標(biāo)，而RMAE是一個局部指標(biāo)，其重要性不如R2和RAAE。但是，即使R2較大，RAAE較小，大的RMAE也能表明在設(shè)計空間的某個區(qū)域的精度較差。在實際應(yīng)用中，經(jīng)常根據(jù)R2的值來判斷代理模型的擬合精度。

建立高精度的近似模型，很大程度上取決于對設(shè)計空間的采樣技術(shù)。用合理的實驗設(shè)計方法均勻分布樣本點，可以有效地體現(xiàn)設(shè)計空間的特征，保證近似模型的精度。文中采用文獻(xiàn)［14］中提出的優(yōu)化拉丁方(optimal latin hypercube，OLH)對空間采樣，生成隨機均勻分布的50個樣本點。

采用上述的3種代理模型，對樣本點進(jìn)行擬合，評價各種模型的精度，如表3所示。由表可見，徑向基模型在整體和局部的精度都很差;而響應(yīng)面和Kriging模型表現(xiàn)出更好的性能，最終選取Kriging模型用于后期的參數(shù)修正。

表3 代理模型的精度評價

5 模型修正結(jié)果

基于建立的車身與底盤間相對位移的Kriging模型，本文中采用序列二次規(guī)劃優(yōu)化算法使模型誤差E(x)最小化。E(x)的優(yōu)化歷程如圖6所示，可見經(jīng)過103步迭代，E(x)由122收斂于25.67。將優(yōu)化得到的參數(shù)代入式(2)和式(3)，得到減振器阻尼和緩沖塊剛度特性。圖7給出了修正前后阻尼和剛度特性。根據(jù)修正結(jié)果，更新多體動力學(xué)模型。圖8為車身與底盤左前懸架處相對位移修正仿真值和實驗值的對比。可以看出，經(jīng)參數(shù)識別后的修正模型仿真結(jié)果與試驗非常接近。

6 結(jié)論

(1)使用正交試驗設(shè)計，從眾多參數(shù)中篩選出對模型準(zhǔn)確性影響顯著的參數(shù)，可以提高模型修正的效率。

(2)采用優(yōu)化拉丁方實驗設(shè)計對設(shè)計空間均勻采樣，并通過3種類型代理模型的對比研究發(fā)現(xiàn)，Kriging模型具有較高的精度。

(3)模型修正后的仿真結(jié)果逼近實驗值，可為后期耐久性的準(zhǔn)確預(yù)測奠定基礎(chǔ)。

(4)利用模型修正技術(shù)獲得汽車懸架非線性動力學(xué)特性，可大大降低開發(fā)成本，縮短開發(fā)周期。該技術(shù)可以推廣至其他相關(guān)領(lǐng)域。

［1］ 高云凱，李翠，崔玲，等．燃料電池大客車車身疲勞壽命仿真分析［J］．汽車工程，2010，32(1):7－12．

［2］ 孟瑾，朱平，胡志剛．基于多體動力學(xué)和有限元法的車身結(jié)構(gòu)疲勞壽命預(yù)測［J］．中國公路學(xué)報，2010，23(4):113－120．

［3］ 趙婷婷，李長波，王軍杰，等．基于有限元法的某微型貨車車身疲勞壽命分析［J］．汽車工程，2011，33(5):428－432．

［4］ Kim H S，Yim H J，Kim CB．Computational Durability Prediction of Body Structure in Propotype Vehicles［J］．International Journal of Automotive Technology，2000，23(1)．

［5］ Rutherford A C，Inman D J，Park G，et al．Use of Response Surface Metamodels for Identification of Stiffness and Damping Coefficients in a Simple Dynamic System［J］．Shock and Vibration，2005(12):317－331．

［6］ Kumar R K，Shankar K．Parametric Identification of Structures with Nonlinearities Using Global and Substructure Approaches in the Time Domain［J］．Advances in Structural Engineering，2009，12(2):195－210．

［7］ Rama Mohan Rao A，Lakshmi K．Parameter Estimation Technique Combining Pod with Hybrid Adaptive Swarm Intelligence Algorithm［C］．Proceedings of International Conference on Advances in Civil Engineering ACE-2011，21－23 October 2011:57－63．

［8］ Mohan L V，Shunmugam M S．An Orthogonal Array Based Optimization Algorithm for Computer-aided Measurement of Worm Surface［J］．International Journal Advanced Manufacture Technology，2006，30(5 －6):434 －443．

［9］ Box G E P，Wilson K B．On the Experimental Attainment of Optimum Conditions［J］．Journal of the Royal Statistical Society B，1951，13(1):1 －45．

［10］ Hardy R L．Multiquadric Equations of Topography and Other Irregular Surfaces［J］．Journal of Geophysical Research，1971，76(8):1905－1915．

［11］ Mcdonald DB，Grantham W J，Tabor W L，et al．Response Surface Model Development for Global/local Optimization Using Radial Basis Functions［C］．Presented at the 8th AIAA Symposium on Multidisciplinary Analysis and Optimization，Long Beach，CA，2000．

［12］ Simpson T W，Mistree F．Kriging Models for Global Approximation in Simulation-based Multidisciplinary Design Optimization［J］．Aiaa Journal，2001，39(12):2233 －2241．

［13］ Cressie N．Geostatistics［J］．The American Statistician，1989，43(4):197－202．

［14］ Jin R，Chen W，Sudjianto A．An Efficient Algorithm for Constructing Optimal Design of Computer Experiments［J］．Journal of Statistical Planning and Inference，2005，134(1):268－287．