張發(fā)平，張書暢，武鍇，張云賀，閻艷

(1.北京理工大學(xué) 機械與車輛學(xué)院，北京 100081;2.北京電子工程總體研究所，北京 100085)

0 引言

坦克作為現(xiàn)代戰(zhàn)爭的主要路基武器載體，行進間射擊精度和行駛平順性是其重要的設(shè)計技術(shù)指標(biāo)［1-3］，這些指標(biāo)要通過對坦克動力學(xué)模型的參數(shù)優(yōu)化來保證?；谔摂M樣機技術(shù)的多體動力學(xué)建模仿真與參數(shù)優(yōu)化是目前研究坦克等高機動履帶車輛動力學(xué)性能的重要方法［4］。在進行參數(shù)優(yōu)化時，需要大量重復(fù)調(diào)用仿真模型以獲取響應(yīng)計算目標(biāo)值。然而由于該類車輛總體復(fù)雜，動力學(xué)仿真模型往往是高度密集計算、求解周期長、難以在有限計算資源內(nèi)大量求解，成為制約參數(shù)優(yōu)化效率和精度的瓶頸之一［5-6］。

使用代理模型技術(shù)構(gòu)建設(shè)計變量和目標(biāo)函數(shù)的近似響應(yīng)面來提高優(yōu)化效率，是現(xiàn)有研究的主流方向之一。如借助近似模型對履帶車輛多體動力學(xué)模型的參數(shù)進行修正來解決模型仿真與試驗結(jié)果的匹配問題［7-8］，基于全局敏感性分析采用代理模型對復(fù)雜非線性機械系統(tǒng)進行參數(shù)優(yōu)化［9］等。然而，這些方法需要大量采樣且優(yōu)化精度不高，因此如何進一步提高效率和精度仍需研究。

基于代理模型的全局優(yōu)化(SGO) 可以有效降低“昂貴模型”評估的次數(shù)，適用于包含高性能仿真模型的設(shè)計優(yōu)化過程［10-11］。傳統(tǒng)的SGO 方法需大量增加樣本來降低代理模型的全局誤差，很容易產(chǎn)生維度災(zāi)難問題。為了解決這一問題，近年來衍生出基于代理模型的全局抽樣優(yōu)化(SGSO) 方法，在SGO方法的基礎(chǔ)上增加代理模型的靜態(tài)和動態(tài)采樣方法更新策略，在優(yōu)化過程中選擇樣本點來更新代理模型。

更新抽樣方法采用基于預(yù)測最優(yōu)解的發(fā)掘方法［12］，將預(yù)測最優(yōu)解作為新增樣本點，通過提高代理模型的精度，來有效避免優(yōu)化收斂得到非最優(yōu)解。基于設(shè)計空間填充的采樣策略，通過引入規(guī)則，將基于預(yù)測最優(yōu)的發(fā)掘法和基于誤差的發(fā)掘法融合。方法包括由Jones［13］提出的基于改善概率的全局優(yōu)化(PGO) 算法和Jones 等［14］提出的基于期望改善度的全局優(yōu)化(EGO) 算法，以及期望改善度(EI) 方法［15］。這些方法雖有較強的全局探索能力，但局部發(fā)掘性能較差，當(dāng)各設(shè)計變量維數(shù)較高時收斂速度會大幅降低。

另外一類提高效率和精度的研究是基于設(shè)計空間縮減的策略方法。主要包括信賴域(TR) 方法和分割平面(CP) 方法。TR 法［16］提出TR 縮放指標(biāo)，在每次迭代中根據(jù)縮放指標(biāo)的大小判斷設(shè)計空間是該收縮還是該擴張。Cheng 等［17］將TR 思想和模式采樣(MPS) 方法結(jié)合，顯著提高了MPS 的高維優(yōu)化效率。TR 法求解局部最優(yōu)解的效率較高，但在全局探索方面的能力稍顯薄弱。CP 法由Wang 等［18］提出，其思想是通過設(shè)計空間中的超平面分割響應(yīng)面，保留代理模型預(yù)測值較小的部分，并在該區(qū)域使用DOE 方法采集新的樣本點，CP 的邊界由最優(yōu)化方法確定。CP 法是單向設(shè)計空間縮減法，無法增大尋優(yōu)子空間，因此很可能陷入局部最優(yōu)。

針對以上問題，本文以坦克動力學(xué)模型參數(shù)優(yōu)化為研究對象，提出一種基于代理模型進化的參數(shù)優(yōu)化方法，旨在確保優(yōu)化精度的前提下大幅減少仿真模型調(diào)用次數(shù)，提高優(yōu)化效率。首先，進行坦克機械系統(tǒng)拓?fù)潢P(guān)系及約束分析，建立坦克行進間多體動力學(xué)仿真模型;其次，以垂向加權(quán)加速度均方根與俯仰角加速度均方根組合構(gòu)建目標(biāo)函數(shù)，并選出對目標(biāo)函數(shù)較敏感的設(shè)計變量;隨后，提出基于代理模型進化的參數(shù)優(yōu)化方法，在重點區(qū)域持續(xù)更新代理模型;最后，使用該方法進行模型參數(shù)優(yōu)化，得到高精度代理模型和最佳參數(shù)組合。

1 坦克行進間多體動力學(xué)建模

1.1 坦克機械系統(tǒng)拓?fù)浣Y(jié)構(gòu)及約束建模

現(xiàn)代履帶坦克部件眾多，組成結(jié)構(gòu)復(fù)雜，性能指標(biāo)多。根據(jù)經(jīng)驗，在此選擇影響性能的技術(shù)指標(biāo):車輛行駛的平順性為主要成分［19-20］。忽略動力和傳動系統(tǒng)，某型號坦克主要由上裝系統(tǒng)和底盤行動系統(tǒng)兩部分組成。上裝系統(tǒng)通過上、下座圈與底盤系統(tǒng)連接。所構(gòu)建的被動懸掛坦克行進間多體動力學(xué)模型共包含619 個部件、660 個約束。部件和約束明細(xì)見表1 和表2。根據(jù)部件之間的連接關(guān)系以及部件之間的約束狀況，簡化后的坦克機械系統(tǒng)拓?fù)浣Y(jié)構(gòu)圖如圖1 所示。這些零部間的約束共包含普通接觸約束24 個，間隙碰撞接觸約束12 個，接觸摩擦184 個，剛度阻尼力矩18 個。其多體動力學(xué)模型的構(gòu)建可通過理論方法描述約束力與狀態(tài)變量間的關(guān)系。具體可用多體動力學(xué)仿真軟件RecurDyn 進行建模仿真，通過添加直線彈簧阻尼器和扭轉(zhuǎn)彈簧阻尼器實現(xiàn)各約束，其中直線彈簧阻尼器和扭轉(zhuǎn)彈簧阻尼器的非線性剛度和阻尼特性數(shù)據(jù)設(shè)置為樣條曲線。

圖1 坦克多體動力學(xué)模型拓?fù)浣Y(jié)構(gòu)圖Fig.1 Topological structure of the multi-body dynamics model

表1 坦克多體動力學(xué)模型部件明細(xì)Table 1 Details of tank multi-body dynamics model cemponents

表2 坦克多體動力學(xué)模型約束明細(xì)Table 2 Constraint details of tank multi-body dynamics model

1.2 路面譜的三維不平度建模及包絡(luò)效應(yīng)處理

1.2.1 路面譜的三維不平度建模

坦克等履帶車輛行駛過程中路面激勵影響車輛的性能。根據(jù)國家標(biāo)準(zhǔn)GB7031—1986 車輛振動輸入路面平度表示方法規(guī)定的路面不平度公式量化如下:

式中:Gq(n) 和Gq(n0) 分別為路面不平度功率譜密度和標(biāo)準(zhǔn)功率譜密度;n 為路面不平度空間頻率;n0為其對應(yīng)的標(biāo)準(zhǔn)空間頻率;w 為路面功率譜密度的頻率指數(shù)。一般采用諧波疊加法［10］建立如下三級路譜。左右相干路面不平度表達式如下:

式中:ql(x) 和qr(x) 分別為左、右路面不平度表達式;ni為將空間頻率nmin＜n ＜nmax平均劃分為N 個區(qū)間中第i 個區(qū)間的中心頻率;αl為［0，1］區(qū)間內(nèi)的隨機數(shù)，表示左側(cè)相位系數(shù);αr為右側(cè)相位系數(shù)，決定了左右兩側(cè)的相干性，其計算由相位角擬合，見文獻［21］;d 為負(fù)重輪輪距;αn為［0，1］區(qū)間內(nèi)的隨機數(shù);Ai表示諧波振動幅值，

Gq(ni) 為功率譜密度，可以由國家標(biāo)準(zhǔn)GB7031 進行計算得到，Δn 為空間頻率間隔。

GB7031 根據(jù)功率譜密度將路面劃分為A～H 8 個等級，H 級的平度系數(shù)的幾何平均值最大，為262 144，代表等級中最嚴(yán)苛的路面。因此，本文建立D、F、H 三個等級路譜，分別代表不同路面不平度的等級進行仿真，如圖2 所示。

圖2 D、F、H 級路面左右相干三維不平度模型Fig.2 Left-right coherent three-dimensional roughness model of class D，F(xiàn) and H pavements

1.2.2 履帶板包絡(luò)效應(yīng)的路面譜濾波建模

履帶裝甲車輛在路面行駛時，因履帶板對地面具有包絡(luò)效應(yīng)，作用在負(fù)重輪上的高程激勵發(fā)生變化。可以將履帶對地面的濾波作用等效成1 階低通濾波器對路面譜的作用，標(biāo)準(zhǔn)形式為

式中:G(ω) 為濾波器幅頻特性，ω 為路面激勵頻率;G0為零頻增益;ωc為濾波器截止角頻率;G(n) 為履帶濾波函數(shù);S0為靜態(tài)靈敏度;nc為空間截止頻率。履帶對地面的濾波函數(shù)與自身結(jié)構(gòu)和車速有關(guān)。通過相應(yīng)等級的路面不平度和對應(yīng)車速下濾波函數(shù)對位移譜進行濾波，得到時域信號作為去履帶后的路面激勵輸入，對履帶車輛懸掛系統(tǒng)做動態(tài)特性分析。

1.3 坦克整車多體動力學(xué)仿真模型

基于機械系統(tǒng)部件拓?fù)潢P(guān)系及約束模型，采用多體動力學(xué)仿真軟件RecurDyn 建立坦克行進間多體動力學(xué)模型。所構(gòu)建的整車RecurDyn 模型如圖3 所示。

圖3 坦克整車多體動力學(xué)RecurDyn 模型Fig.3 Tank multi-body dynamics RecurDyn model

2 懸掛系統(tǒng)參數(shù)優(yōu)化目標(biāo)函數(shù)及變量確定

2.1 懸掛系統(tǒng)參數(shù)優(yōu)化目標(biāo)函數(shù)

坦克等履帶車輛行駛過程中，駕駛艙地板位置的垂向加速度和坦克車體的俯仰角加速度是重要的評價指標(biāo)［20，22-23］，因此將坦克在不同路面上行駛過程中車體駕駛艙位置的垂向加速度均方根aya和俯仰角加速度均方根值azw加權(quán)，作為車輛性能綜合指標(biāo)，共同組成如下優(yōu)化目標(biāo)函數(shù):

式中:wya和wzw分別為兩者的權(quán)重，用來平衡兩目標(biāo)項的數(shù)量級。該目標(biāo)函數(shù)沒有具體物理意義，只是保證兩者之和的最小化。目標(biāo)函數(shù)的取值越小，車輛行進間的綜合性能越好。

垂向加速度均方根aya和俯仰角加速度均方根值azw的計算根據(jù)國際標(biāo)準(zhǔn)ISO 2631-1:1997 人體承受全身振動評價進行，此標(biāo)準(zhǔn)認(rèn)為人體對0.5～80 Hz區(qū)間的不同頻率振動的敏感程度不同。因此采用加權(quán)加速度均方根值對舒適性進行量化，通過快速傅里葉變換(FFT) 將仿真得到的駕駛艙位置垂向加速度時域譜變換為功率譜

式中:Ga(f) 為垂向加速度頻率譜;aya(t) 為垂向加速度時域譜。在此基礎(chǔ)上，通過ISO2631-1 規(guī)定的0.5～80 Hz 區(qū)間頻率加權(quán)函數(shù)對頻譜加權(quán)，計算均方根值如下:

式中:F 為總頻率尺度;ω(f) 為頻率分段加權(quán)函數(shù):

聯(lián)合式(7)～式(8)，可以得到垂向加權(quán)加速度功率譜密度。

坦克車體的俯仰角加速度時域均方根值量化計算方法如下:

式中:azw(t) 為車體俯仰角加速度時域譜;T 為總時間歷程。

工程上這些參數(shù)均通過坦克車輛主動懸掛系統(tǒng)的設(shè)計參數(shù)或控制參數(shù)來控制，因此選擇懸掛系統(tǒng)參數(shù)作為優(yōu)化變量來進行優(yōu)化。

2.2 懸掛系統(tǒng)參數(shù)優(yōu)化設(shè)計變量

坦克等履帶車輛的參數(shù)優(yōu)化效率受設(shè)計變量的維度和設(shè)計空間的影響，高維度的設(shè)計變量和大尺寸的設(shè)計空間常常導(dǎo)致計算成本驟升。因此為了提高優(yōu)化效率，根據(jù)車輛動力學(xué)相關(guān)理論，選擇6 個對目標(biāo)函數(shù)敏感的懸掛系統(tǒng)參數(shù)作為設(shè)計變量來控制設(shè)計變量維度，分別為:左右兩側(cè)1、6 號扭桿彈簧剛度系數(shù)，左右兩側(cè)2、5 號扭桿彈簧剛度系數(shù)，左右兩側(cè)3、4 號扭桿彈簧剛度系數(shù)，左右兩側(cè)1、6 號扭桿彈簧阻尼系數(shù)，左右兩側(cè)2、5 號扭桿彈簧阻尼系數(shù);左右兩側(cè)3、4 號扭桿彈簧阻尼系數(shù)。另一方面，通過基于代理模型進化的參數(shù)優(yōu)化方法在參數(shù)優(yōu)化過程種持續(xù)縮小設(shè)計空間尺寸，節(jié)約計算資源，提高優(yōu)化效率。6 個設(shè)計變量的位置示意圖如圖4 所示。圖4 中，a1～a6分別為相應(yīng)車輪間的水平距離，mb為車體質(zhì)量，zb為車體的垂向位移，zw1～zw6分別為相應(yīng)車輪相的垂向位移，zt1～zt6分別為相應(yīng)車輪位置處路面的激勵位移。

圖4 優(yōu)化變量位置示意圖Fig.4 Schematic diagram of optimized variable location

王欽龍等［8］在履帶車輛模型修正工作中通過極差法確定了對目標(biāo)函數(shù)較敏感的設(shè)計變量分別是車體俯仰轉(zhuǎn)動慣量Izz、懸掛系統(tǒng)剛度k 和懸掛系統(tǒng)阻尼c。車體轉(zhuǎn)動慣量Izz受整車各部件質(zhì)量、質(zhì)心等多因素影響，不具備優(yōu)化價值。

2.3 懸掛系統(tǒng)參數(shù)優(yōu)化聯(lián)合仿真過程

基于代理模型進化優(yōu)化算法的履帶車輛參數(shù)優(yōu)化過程由算法程序、控制程序和坦克行進間多體動力學(xué)模型聯(lián)合仿真環(huán)節(jié)所組成。其中，算法程序為采用MATLAB 軟件自研的代理模型進化優(yōu)化算法，控制程序依托Simulink 由RecurDyn 聯(lián)合仿真模塊和其他基礎(chǔ)組件構(gòu)成，實現(xiàn)坦克行進間多體動力學(xué)模型的參數(shù)傳入和響應(yīng)采集功能，以動力學(xué)仿真軟件RecurDyn 為基礎(chǔ)平臺，根據(jù)第1 節(jié)提供的機械系統(tǒng)拓?fù)浼凹s束信息實現(xiàn)坦克行進間多體動力學(xué)模型建模，以及3 種三維路面譜建模和行進動力學(xué)仿真。具體結(jié)構(gòu)如圖5 所示。

3 初始代理模型構(gòu)建

坦克行進間多體動力學(xué)仿真模型部件及約束眾多，單次仿真求解效率極低(本文仿真配置計算機下一次仿真約20 min) 。為解決該問題，本節(jié)將構(gòu)建初始代理模型表征設(shè)計參數(shù)和目標(biāo)函數(shù)之間的響應(yīng)關(guān)系，為基于代理模型進化的參數(shù)優(yōu)化方法確定優(yōu)化搜索方向。

初始訓(xùn)練樣本集的空間均布性決定了初始代理模型能否有效地表征設(shè)計變量和目標(biāo)值在整個設(shè)計空間的響應(yīng)面關(guān)系。針對2.2 節(jié)給出的坦克優(yōu)化設(shè)計變量及對應(yīng)的可行空間，使用基于隨機進化算法的最優(yōu)拉丁立方(ESEA-OLHD) 法進行試驗設(shè)計，獲得6D(D 為設(shè)計變量維度) 個樣本，通過2.3 節(jié)聯(lián)合仿真計算得到響應(yīng)，生成初始訓(xùn)練集。使用擴展徑向基函數(shù)(RBF) 代理模型［24］構(gòu)建6 維設(shè)計變量與目標(biāo)之間的初始響應(yīng)面。擴展RBF 代理模型詳細(xì)過程見文獻［11］。

在D、F、H 三個等級路面條件下，ESEA-OLHD試驗設(shè)計法采樣獲得訓(xùn)練樣本集，據(jù)此構(gòu)建三級路面6 維設(shè)計變量和目標(biāo)函數(shù)間的初始代理模型，如圖6 所示為以圖形方式表示的初始代理模型，因篇幅原因只列出D 級路面。圖6 中，K1為設(shè)計變量中的剛度，C1為設(shè)計變量中的阻尼。初始代理模型將在優(yōu)化過程中持續(xù)更新，不斷提高擬合精度。

4 基于代理模型進化的參數(shù)優(yōu)化方法

為了解決傳統(tǒng)工程設(shè)計優(yōu)化問題大量調(diào)用仿真模型導(dǎo)致的優(yōu)化效率低問題，本文提出設(shè)計空間的三層分解方案、基于代理模型進化的參數(shù)優(yōu)化方法，以及在可能最優(yōu)解的重點區(qū)域的代理模型動態(tài)更新，從而提高優(yōu)化精度，降低樣本使用量。

4.1 多層設(shè)計空間縮減策略

在工程設(shè)計優(yōu)化初始階段，設(shè)計變量范圍設(shè)定往往較大，以便能覆蓋所有可能的設(shè)計結(jié)果，但帶來的后果是大設(shè)計范圍導(dǎo)致空間搜索成本過高。本節(jié)將設(shè)計變量所構(gòu)成的空間劃分為三層空間，分別為初始設(shè)計空間(OA)、局部過渡子空間(MA) 和重要子空間(IA) 。各子空間關(guān)系示意圖如圖7 所示。

圖7 某次迭代過程中三層子空間示意圖Fig.7 Schematic diagram of three-layer subspaces in an iterative process

4.1.1 初始設(shè)計空間(OA)

OA 是優(yōu)化變量的初始范圍，也是首次迭代的基準(zhǔn)空間，覆蓋變量所有可能取值。在OA 空間中通過極值法和高斯誤差法的動態(tài)采樣策略，可以改善代理模型在整個初始設(shè)計空間的精度，通過基于ADAM 梯度下降的局部尋優(yōu)法和松弛判別式尋找可能被遺漏的極小值，為算法提供跳出局部最優(yōu)解的能力。

4.1.2 局部過渡子空間(MA)

MA 是過渡子空間，可由圖8 所示的二次聚類法把OA 劃分為多個子MA 空間。每一個MA 空間包含優(yōu)化問題的一個或多個局部極小值。在各MA空間中，采用ADAM 梯度下降法尋找MA 空間中的局部極小值并更新代理模型，確定本次迭代的預(yù)測最優(yōu)解。

圖8 基于二次聚類確定過渡區(qū)間流程圖Fig.8 Flow chart of determining transition intervals based on quadratic clustering

使用確定MA 子空間的步驟如下:

1) 首次聚類種群驗證。通過基準(zhǔn)空間樣本點數(shù)量確定聚類數(shù)量區(qū)間，以區(qū)間中的每一個取值為先驗值，執(zhí)行模糊C 均值(FCM) 聚類，計算每一次聚類結(jié)果的評價指標(biāo)，衡量聚類種群的有效性。將評價指標(biāo)最優(yōu)的聚類數(shù)量值作為最佳聚類數(shù)量，并保存最佳聚類中心坐標(biāo)。

2) 二次聚類劃分種群。以最佳聚類數(shù)量和最佳聚類中心坐標(biāo)為先驗條件，對樣本集進行基于核函數(shù)的模糊C 均值(KFCM) 聚類，得到cbest個簇，每個簇都包含局部極小值。選擇平均響應(yīng)較低的前50%(根據(jù)經(jīng)驗取值) 簇作為各MA 子空間，根據(jù)各簇包含樣本各維坐標(biāo)的最大、最小值確定子空間范圍。文獻［25 -26］提供了FCM 和KFCM 的詳細(xì)執(zhí)行過程。

4.1.3 重要子空間(IA)

IA 是最可能包含全局最優(yōu)解的區(qū)域，是設(shè)計空間的進一步收縮。通過ADAM 梯度下降法在IA 子空間中開展優(yōu)化搜索工作，使用基于極值法和高斯誤差聯(lián)合采樣策略在IA 子空間中重點采樣，在可能存在全局最優(yōu)解的鄰域深入探索并持續(xù)精化代理模型。確定IA 的具體操作步驟如下:

1) 通過二次聚類法確定若干過渡空間MA;

2) 對所有MA 空間中的局部最小值排序，得到最佳預(yù)測值作為IA 子空間的中心;

3) 以IA 空間隸屬的MA 子空間邊界和最佳預(yù)測值的各維度坐標(biāo)為依據(jù)，計算IA 子空間各維度邊界。MA 數(shù)量大于1 或等于1 兩種情況下IA 區(qū)間邊界IA_range 的計算方法如下:

式中:prebest為當(dāng)前最優(yōu)預(yù)測值所在點的各維度坐標(biāo)向量;MAurapnge和MAlroawnge分別為MA 空間各維度上界和下界;size(MA) 為MA 空間的數(shù)量。通過在IA 子空間中進行多種形式，更多數(shù)量的采樣，提高代理模型在該區(qū)域的精度，進而促進全局最優(yōu)解的捕獲，且IA 區(qū)間的范圍相對于MA 區(qū)間得到了進一步縮減。多層空間縮減的算法偽代碼見附錄1。

4) 根據(jù)搜索效率和精度的綜合考慮，每3 次迭代后進行IA 空間確定，防止過早收斂。

4.2 多起點局部尋優(yōu)及代理模型動態(tài)更新

4.2.1 基于ADAM 梯度下降法的局部尋優(yōu)

優(yōu)化迭代過程中，需要在MA 和IA 子空間進行局部尋優(yōu)，獲得局部極小值、確定新增樣本點，更新代理模型。使用用于深度學(xué)習(xí)的ADAM 梯度下降法［27］，從多個起點進行搜索來捕獲空間中的多個局部極小值。

在MA 空間和IA 空間執(zhí)行基于ADAM 的局部尋優(yōu)算法，流程見附錄2。需將尋優(yōu)搜索限制在MA或IA 空間范圍中，在搜索過程中接近上下邊界的變量維度將停留在邊界，其他維度依然保持向負(fù)梯度方向搜索，直至到達預(yù)設(shè)最大迭代次數(shù)stepnum或滿足收斂條件‖deltaX‖2＜ε，ε 為一個足夠小的值，一般取ε×10-8。

4.2.2 基于動態(tài)采樣策略的代理模型更新

新增樣本點是優(yōu)化迭代進程中代理模型更新的重要來源，本節(jié)提出基于極值法和高斯誤差法聯(lián)合的動態(tài)采樣策略。

4.2.2.1 基于極值法的新增樣本點

該方法將獲得子空間中可能是全局最優(yōu)點的候選點，執(zhí)行步驟如下:

步驟1確定各空間(如MA) 包含的觀測樣本點集以及邊界。

步驟2在空間邊界內(nèi)從多起點進行ADAM 梯度下降搜索，得到各子空間內(nèi)所有局部極小值，將其按預(yù)測響應(yīng)從小到大排序后取各子空間前p 個局部極小值，p 的取值如下:

式中:s 為本次迭代生成的子空間數(shù)量。需要說明的是，p 和s 的取值是考慮到效率精度的平衡。

步驟3剔除新增樣本點集中的兩類冗余:1) 從多起始點搜索得到的局部極小值可能收斂到同一位置;2) 得到的新增樣本點可能與原有樣本集中的點重合。對此，計算新增樣本點集中各點歐式距離，剔除距離相近的新增樣本點;逐一計算新增樣本集和原始觀測樣本集的歐氏距離，剔除冗余點。歐式距離評價公式如下:

式中:x1、x2為兩個樣本點;‖＊‖2為向量L2 范數(shù);LHSrange(1，:) 和LHSrange(2，:) 分別為初始設(shè)計空間各維度的上界和下界。

4.2.2.2 基于高斯誤差法加點

該方法將獲得該次迭代關(guān)注區(qū)域內(nèi)代理模型誤差較大位置的樣本點。高斯誤差的計算包含三部分:1) 利用空間(如IA) 樣本集訓(xùn)練該區(qū)域的高斯過程回歸模型(GRF)，得出IA 空間各處高斯過程預(yù)測值IA_GRF_pre 和方差σ;2) 計算IA 子空間各處GRF 和RBF 預(yù)測響應(yīng)偏差Δpre=|IA_GRF_pre-IA_RBF_pre|;3) 分別將高斯過程方差和預(yù)測響應(yīng)偏差進行標(biāo)準(zhǔn)化:

式中:σnor為標(biāo)準(zhǔn)化方差;Δprenor為預(yù)測響應(yīng)偏差。以上兩個參數(shù)加權(quán)，得到總誤差值SE=0.2Δprenor+0.8σnor，選擇總誤差值最小的3 個點作為新增樣本集。

高斯過程方差和高斯過程回歸模型的理論內(nèi)容見文獻［28］。優(yōu)化迭代中動態(tài)采樣策略在MA 和IA 子空間采樣流程見附錄1 的多層空間縮減策略算法和附錄2 的基于ADAM 梯度下降的局部尋優(yōu)法算法。

4.3 基于動態(tài)采樣的局部最優(yōu)解跳出

在滿足判別標(biāo)準(zhǔn)下在OA 空間中使用動態(tài)采樣法更新代理模型避免算法陷入局部最優(yōu)解，偽代碼如附錄3，分為以下4 個步驟:

步驟1執(zhí)行判斷:若某次迭代新增樣本點為0或觀測最優(yōu)值連續(xù)3 次迭代中沒有發(fā)生變化，則判定當(dāng)前迭代過程是否陷入局部最優(yōu)。

步驟2執(zhí)行動態(tài)采樣:在OA 中執(zhí)行高斯誤差采樣和極值采樣，得到預(yù)測響應(yīng)誤差較大和局部最優(yōu)位置的新增樣本點。

步驟3更新基準(zhǔn)空間:若OA 中最小的局部極小值OAbest與當(dāng)前全局最優(yōu)解prebest(iter) (iter 為迭代次數(shù)) 滿足松弛判別式(11)，則更新全局最優(yōu)解，合并基準(zhǔn)空間OAbest的坐標(biāo)作為下次迭代基準(zhǔn)空間邊界。

4.4 優(yōu)化流程

基于代理模型進化的參數(shù)優(yōu)化方法流程如圖9所示。流程包含4 個部分:1) 循環(huán)迭代前，初始化參數(shù)并訓(xùn)練初始代理模型;2) 進入迭代后，使用多層空間縮減策略確定三類子空間，并從三類子空間中交替局部尋優(yōu)并獲取新增樣本點;3) 執(zhí)行跳出局部最優(yōu)解的措施，確定新增樣本集;4) 更新基準(zhǔn)空間和代理模型，計算目標(biāo)值并判斷是否終止。

圖9 坦克參數(shù)優(yōu)化總體技術(shù)流程Fig.9 Overall technical process of tank parameter optimization

5 基于進化優(yōu)化算法的懸掛系統(tǒng)參數(shù)優(yōu)化

5.1 坦克懸掛系統(tǒng)設(shè)計參數(shù)優(yōu)化過程

分別針對D 級、F 級和H 級3 種道路環(huán)境進行基于代理模型進化的參數(shù)優(yōu)化，以盡可能少的樣本得到單路面情況下的最佳設(shè)計變量組合以及單路面情況下設(shè)計變量和目標(biāo)值的高精度響應(yīng)面。其中，D 級路面優(yōu)化問題的數(shù)學(xué)表達式如下:

在本文優(yōu)化問題中，多層設(shè)計空間縮減策略、基于代理模型進化的參數(shù)優(yōu)化方法各子方法參數(shù)以及RecurDyn 行進間多體模型仿真參數(shù)設(shè)置如表1 所示。

表1 優(yōu)化過程中算法及仿真模型參數(shù)設(shè)置Table 1 Algorithm and simulation model parameter setting

5.2 3 種路面條件的優(yōu)化結(jié)果分析

圖10 為D 級路面10 次優(yōu)化中最小全局最優(yōu)解的收斂曲線。由圖10 可知，每3 次迭代目標(biāo)值有較大幅度下降，第8～13 次迭代中陷入局部最優(yōu)，算法采取相應(yīng)策略跳出局部最優(yōu)，于14 次迭代中得到更小的目標(biāo)值。表2 所示為10 次優(yōu)化中坦克在D 級、F 級和H 級路面下優(yōu)化前后懸掛系統(tǒng)參數(shù)優(yōu)化的目標(biāo)值，變化率分別為［27.2%，34.8%］、［23.3%，26.6%］、［19.6%，21.3%］，樣本使用量分別為［91，108］、［91，110］、［93，109］。表3 為優(yōu)化前后的參數(shù)值。表2 表明，該方法的優(yōu)化參數(shù)可使得3 種路面下行進間車輛性能綜合指標(biāo)平均提升分別為32.4%、24.5%和20.4%，且可綜合將傳統(tǒng)方法平均樣本使用量的約500～1 000 次降低到80～150 次，在同等算力下節(jié)省時間100 h 以上(以本文所用計算機配置單次仿真周期20 min 計算) 。

表2 基于代理模型進化優(yōu)化算法的10 次單路面優(yōu)化(懸掛系統(tǒng)參數(shù)優(yōu)化目標(biāo)函數(shù)值)Table 2 10 times single road optimization based on agent model evolutionary optimization algorithm(objective function value of suspension parameter optimization)

表3 3 種路面下的最優(yōu)參數(shù)組合Table 3 Optimal parameter combination for each grade of pavement

圖10 D 級路面譜優(yōu)化收斂曲線Fig.10 Convergence curve of spectrum optimization of grade D pavement

使用H 級路面下的最優(yōu)參數(shù)組合和初始參數(shù)組合進行仿真，得到坦克行進間加速度功率譜對比曲線如圖11 所示，從中可以看出優(yōu)化后0～50 Hz振動頻率區(qū)間內(nèi)功率譜密度均值和極值都有了顯著降低。

圖11 H 級路面下優(yōu)化前后功率譜密度對比Fig.11 Comparison of power spectral densities before and after optimization for H grades of pavement

6 結(jié)論

本文以履帶車輛多體動力學(xué)模型參數(shù)優(yōu)化為研究對象，提出了基于代理模型進化算法的參數(shù)優(yōu)化方法，以提高復(fù)雜多體動力學(xué)系統(tǒng)優(yōu)化的效率，并用懸掛系統(tǒng)的參數(shù)優(yōu)化進行了驗證。得出主要結(jié)論如下:

1) 提出了復(fù)雜結(jié)構(gòu)履帶車輛動力學(xué)參數(shù)優(yōu)化的代理模型進化優(yōu)化方法，可實現(xiàn)異構(gòu)模型間信息無縫傳遞，以及代理模型進化和參數(shù)優(yōu)化的有效融合，提高了參數(shù)優(yōu)化的精度和優(yōu)化過程的效率。

2) 在3 種路面下對坦克動力學(xué)數(shù)值仿真的結(jié)果表明，本文提出的基于代理模型進化的參數(shù)優(yōu)化方法與傳統(tǒng)優(yōu)化方法相比，優(yōu)化精度提高21.3%，降低樣本使用量31.5%。同時，該方法的優(yōu)化參數(shù)可使得3 種路面下行進間車輛性能綜合指標(biāo)平均提升分別為32.4%、24.5%和20.4%。

3) 本文優(yōu)化方法本質(zhì)上是一種單目標(biāo)優(yōu)化法，在求解多目標(biāo)優(yōu)化問題時，需要通過權(quán)重分配系數(shù)建出單目標(biāo)代理模型后再求解。下一步研究中將引入多目標(biāo)概率期望改善等方法，使算法具有本質(zhì)多目標(biāo)優(yōu)化能力，并結(jié)合具體的工程實際驗證。

附錄1

多層空間縮減策略算法

Begin

初始化:初始采樣個數(shù)m 和維度D，其中樣本個數(shù)m 推薦的數(shù)量為二維問題［20，40］，當(dāng)維度2 ＜n ＜10 時，推薦采樣數(shù)量為［6n，8n］;n ＞10 時推薦［50，70］;RBF 代理模型及訓(xùn)練迭代次數(shù)為80，初始各個維度可行設(shè)計空間為LHSrange;

得到昂貴初始樣本集:

Xolh 優(yōu)化拉丁超立方采樣在LHSrange 中得到樣本點矩陣mxD

Yexp 通過復(fù)雜模型計算初始樣本集Xolh 的目標(biāo)值，即響應(yīng)值，矩陣大小為m×1

Sexp 將初始樣本點坐標(biāo)和其對應(yīng)的響應(yīng)值合并為初始昂貴樣本集矩陣，大小為m×(D+1)劃分訓(xùn)練集和測試集:

Sexp_train 和Sexp_test 通過留一交叉驗證法將初始昂貴樣本集劃分為訓(xùn)練集和驗證集，用來更高精度地訓(xùn)練代理模型

訓(xùn)練代理模型:

RBF(Sexp_train，Sexp_test) 訓(xùn)練RBF 神經(jīng)網(wǎng)絡(luò)作為代理模型，并通過交叉驗證提高代理模型精度

save NET_RBF_iter1.mat-mat net 保存該次迭代過程中訓(xùn)練好的代理模型

End

空間縮減:

Input: Sexp，NET_RBF_iter1.mat，iter

初始化: FCM 和KFCM 最大迭代次數(shù)為100，終止條件為1 ×10-5，KFCM 使用高斯核函數(shù)，參數(shù)sigma 設(shè)置為150;

MA——首先進行二次聚類，劃分種群過后進行一定篩選工作，得到s 個MA

Two-stage clustering:

Cbest——最佳聚類個數(shù)

Para_miu_best——種群驗證后得到的與最佳聚類個數(shù)對應(yīng)的聚類中心

Label_KFCM——得到種群劃分矩陣，存儲類表標(biāo)識，數(shù)據(jù)點坐標(biāo)以及觀測響應(yīng)

The determination of MA——將二階段聚類得到的種群通過響應(yīng)值的大小進行排序，根據(jù)最小響應(yīng)值從小到大的順序?qū)⑺蟹N群排序，取前20%的s 個種群作為s 個過渡區(qū)間。如果種群數(shù)量小于5，則只選出1 個種群作為MA。

MA_plt——得到每個過渡區(qū)間包含的樣本點

MA_range——得到每個過渡區(qū)間對應(yīng)的可行空間范圍

IA——在確定MA 的基礎(chǔ)上，通過分層尋優(yōu)法判斷哪個空間是IA 空間

The determination of MA——確定s 個過渡區(qū)間以及響應(yīng)的樣本點和空間范圍

分層尋優(yōu)法——通過分層尋優(yōu)法尋找局部極小值

Local_min_MA——得到每個過渡區(qū)間內(nèi)的局部極小值

The determination of IA

IA_plt——得到重要空間包含的所有樣本點

IA_best——得到重要空間的最佳預(yù)測值

ImpAre_range——得到重要空間的區(qū)間范圍

Output: MA_plt and MA_range/IA_plt，IA_best and ImpAre_range

End

附錄2

基于ADAM 梯度下降的局部尋優(yōu)法算法

初始化:學(xué)習(xí)率，最大迭代步長，初始參數(shù)，衰減速率

Input:

MA_plt and MA_range/ IA_plt and ImpAre_range——MA

空間或IA 空間包含的樣本點以及對應(yīng)的空間范圍，假設(shè)某一空間中包含樣本點S 個

NET_RBF_iter.mat——本次迭代訓(xùn)練的RBF 模型

Whileiter_opt ＜step_num——設(shè)置最大循環(huán)次數(shù)為step_num

grad——計算數(shù)值梯度

Square gradient——累積平方梯度

更新參數(shù)

deltaX=-grad* 學(xué)習(xí)率——計算該次迭代的樣本點各維坐標(biāo)增量

Xiter_opt+1=xiter_opt +deltaX() ——更新尋優(yōu)樣本點的坐標(biāo)

If ‖deltaX‖2＜1 ×10-8——所有維度增量都小于1 -e8，終止循環(huán)

Break;

End if

End while

Output: Local_min——得到大小為S×(D+1) 的局部極小值矩陣

End

附錄3

優(yōu)化流程及跳出局部最優(yōu)策略算法總體迭代流程

While iter ＜itermax——進入循環(huán)體，itermax 為最大迭代次數(shù)

Iter=Iter++——迭代次數(shù)+1

If mod (iter，3)==0——如果當(dāng)前迭代次數(shù)iter 是3 的整數(shù)倍

IA_plt，IA_best，BS_range——執(zhí)行多層空間縮減算法，得到重要空間IA，下次迭代基準(zhǔn)空間以及包含的觀測樣本點

Sample_best——執(zhí)行基于極值的加點方法，得到包含可能全局最優(yōu)點的新增樣本點

Sample_gua——執(zhí)行高斯誤差加點，得到IA 空間誤差較大區(qū)域的新增樣本點

Sample_new——合并Sample_best 和Sample_gua 后，執(zhí)行樣本點篩選方法，剔除內(nèi)外部距離過近的樣本，得到最終新增樣本集

Else

MA_plt，BS_range——執(zhí)行多層空間縮減算法，得到過渡空間MA

Sample_best——執(zhí)行基于極值的加點方法

Sample_new——執(zhí)行樣本點篩選方法

判斷并執(zhí)行未知區(qū)域探索部分

If iter ＞5——5 次迭代后，計算當(dāng)前最優(yōu)預(yù)測值和前5 次最優(yōu)預(yù)測值的差CRT

CRT=Pre_best(iter) -Pre_best(iter-5)

Else

CRT=1 ×1010

End if

If size(Sample_new)=0 or CRT ＜1 ×10-6——如果本次迭代過生成新增樣本點個數(shù)為0 或CRT 小于閾值

執(zhí)行稀疏采樣方法在OA 中采樣

Sample_OA——OA 空間中執(zhí)行高斯誤差采樣，生成新增樣本點

OA_best——OA 空間中部分觀測點為起始點使用尋優(yōu)方法確定OA 最優(yōu)預(yù)測值

Sample_new=［Sample_new;Sample_OA］——合并新增樣本點

If OA_best-pre_best(iter) ＜pre_best(iter)/pre_best(1)

Pre_best(iter)=OA_best——更新最優(yōu)預(yù)測值

BS_range——擴大并更新下次迭代基準(zhǔn)空間

Sample_new=［Sample_new;Sample_OA;OA_best］——合并新增樣本點

End if

End if

Sample_new_rdc——對新增樣本點集執(zhí)行去冗余操作，生成無冗余樣本矩陣

Ysample_new——使用復(fù)雜模型計算新增樣本點的觀測值

Sexp——［Sexp;［Sample_new_rdc，Ysample_new］］將新增樣本坐標(biāo)和響應(yīng)合并到整體樣本集中

NET_RBF_iter+1.mat 更新RBF 代理模型

End While

Out put:pre_best，RBF——最終輸出最優(yōu)解以及最終RBF 代理模型