基于代理模型的多體動力學模型修正技術

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

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

前言

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

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

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

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

1 模型修正方法

模型修正過程首先要求建立一個目標函數，即模型誤差評價函數，它反映了一個物理系統的數值模型和試驗模型之間在輸出響應上的差異。一個最小二乘問題的目標函數可以定義為

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

本文中將結合實驗設計、代理模型技術和二次規劃算法，實現用于耐久性分析的動力學模型的修正，流程如圖2所示。

2 多體動力學建模及參數化

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

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

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

3 篩選參數

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

表1 正交試驗設計參數信息

4 代理模型

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

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

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

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

式中:μ為多項式模型;‖x－xi‖是點x與xi的歐氏距離;φ是基函數;λi是基函數在取樣點xi處的未知加權系數。因此，徑向基函數模型實際上是由基于歐氏距離的徑向對稱基函數和多項式函數的線性組合。

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

表2 基函數

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

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

式中:θk為用于擬合近似模型的未知相關系數;xi

k和xj

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

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

其中:

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

}

s與未知位置x所組成的相關向量:

其中:

詳細過程參考文獻［13］。

代理模型僅僅是對復雜仿真分析模型的簡化和近似，因此，在求出代理模型的待估參數后，須對代理模型進行統計檢驗，評估代理模型對真實響應的逼近程度。在精度不高的代理模型基礎上獲得的最優結果不但不能指導設計，反而會誤導設計，因此，對代理模型的誤差進行評估至關重要。常用的誤差估計準則有:決定系數 R2、相對平均絕對誤差(RAAE)和相對最大絕對誤差(RMAE)。

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

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

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

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

表3 代理模型的精度評價

5 模型修正結果

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

6 結論

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

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

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

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

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

［2］ 孟瑾，朱平，胡志剛．基于多體動力學和有限元法的車身結構疲勞壽命預測［J］．中國公路學報，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．