基于自適應近似模型的GTS模型低風阻尾板優化

胡興軍 劉一塵 李金成 蘭巍 張揚輝 王靖宇?

(1．吉林大學 汽車仿真與控制國家重點實驗室，吉林 長春 130012；2．吉林大學 汽車工程學院，吉林 長春 130012)

隨著能源問題的不斷凸顯，商用車因其能耗較高、受空氣阻力影響大的特點，逐漸成為氣動減阻領域的研究熱點。Hakansson[1]通過改變重型卡車貨箱外形及增加各種氣動附加裝置，發現合理布置的氣動附加裝置可以達到與流線型貨箱同等的氣動減阻效果，而并不壓縮貨箱容積，從而說明了氣動附加裝置的重要性。其中，尾板由于形狀簡單、容易布置，且減阻潛能大而獲得了廣泛關注。Salati等[2]通過對某汽車列車梯形尾板及圓錐形氣動附加裝置進行研究和優化，發現當梯形尾板傾角為13°時，汽車列車氣動阻力系數隨著梯形尾板的板長增加呈現出先增加后降低的趨勢；Ellis等[3]通過在尾板中間增加副尾板，使得某汽車列車在0°～5°偏航角的平均氣動阻力系數降低最高達26%；Hassaan等[4]基于無車輪GTS模型增加尾板結構，并以尾板長度和傾角作為變量探究最優阻力方案減阻機理，發現尾板越長氣動阻力系數越小，傾角為15°時氣動阻力系數最小，并指出C柱渦引起的上洗氣流變化是控制尾部上洗渦、下洗渦分布的關鍵因素；國內的胡興軍等[5]通過對比5種不同的尾板方案指出，推遲氣流分離的尾板方案相比產生分離渦的尾板方案可以獲得更大的氣動減阻效果。

另一方面，氣動減阻中的試驗設計(DOE)和優化也是研究重點。張凱歌[6]通過使用靜態SVR(支持向量回歸)近似模型，對重型卡車尾板進行優化，取得了氣動減阻8.2%的成果；Doyle等[7]應用遺傳算法(GA)直接驅動CFD程序，從而實現對二維卡車尾板幾何的自動尋優，但其執行的GA算法相當于每次產生300個樣本點，對于計算量巨大的三維流場仿真很難應用?？傮w來講，先構建靜態近似模型，再利用某種尋優算法對近似模型進行尋優獲得最優點是優化設計的常見思路。但是，靜態近似模型主要依靠較多的樣本點以實現全局擬合精度，存在數據集獲取時間長、優化效率低、樣本量計算困難的問題。因此，部分學者轉而研究追求最優點擬合精度的自適應近似模型。Jeyasingham等[8]利用ExaADAPT構建自適應采樣法，結合RSM(響應面方法)近似模型對某皮卡車型的30個參數進行了氣動減阻優化研究，發現該自適應方法可以加速RSM近似模型的收斂過程，并獲得預測精度更高的RSM近似模型；Huang等[9]通過構建自適應SVR近似模型對T型管液壓成型加載路徑進行優化，自適應SVR近似模型只迭代27次便收斂得到Pareto前沿。由此說明自適應近似模型具有較大的優化效率。

本文使用粒子群算法(PSO)完成最小二乘支持向量回歸(LSSVR)模型的調參過程，構建自適應PSO-LSSVR近似模型，并采用全局和局部自適應近似模型的方式進一步提高最優解區域的擬合精度，盡量避免陷入局部最優解，同時基于Branin函數為測試函數進行驗證；最后將自適應PSO-LSSVR近似模型用于GTS模型低風阻尾板的優化設計上。

1 自適應PSO-LSSVR近似模型的建立

1.1 LSSVR算法

LSSVR算法是將SVR中求解二次規劃的問題轉換為求解線性方程組的問題，由此可以顯著加快收斂速度[10]，算法過程如下。

對于訓練集S=(xi,yi),i=1,2,…,n，其中xi∈Rm，作為訓練集的m維輸入；yi∈R，作為一維輸出。在SVR中，將原樣本空間中的非線性函數估計問題轉化為高維特征空間的線性函數估計問題的表達式如下：

(1)

式中：w=(w1,w2,…，wn)T,為權值系數向量；φ(x)=(φ1(x),…,φn(x))T，為輸入映射函數；b為實常量。利用結構風險最小化原理求解，將問題轉換為一個等式約束優化問題：

(2)

s.t.yi=wTφ(xi)+b+ei，i=1,2，…,n。

式中，γ為正則化系數，ei為誤差項，可在后式中約去。求解上述優化問題，使用Lagrange等式并根據KKT條件得：

(3)

式中，En為n階單位矩陣，1n=(1,1,…,1)T，y=(y1,y2,…,yn)T，α=(α1,α2,…,αn)T,αi是拉格朗日乘子，Ω=k(xi,xj),k(xi,xj)=φT(xi)φ(yj)，i、j=1,2,…,n。求解α和b后，LSSVR 算法可以表述為

(4)

為了避免“維數災難”問題，LSSVR引入了核函數[11]。常用的核函數包括線性核函數、多項式核函數、RBF核函數及Sigmoid核函數等。RBF核函數具有普適性強的特點，可以適用于任意分布的樣本集[12]，因此本文構建的LSSVR模型采用RBF核函數。RBF核函數表達式如下：

(5)

式中，σ為核函數的寬度。本文的LSSVR模型使用MATLAB LSSVM Toolbox V1.6構建，在該工具箱中，LSSVR模型的σ2和Γ(懲罰因子)決定擬合和泛化能力。已有的文獻采用交叉驗證法[12]、網格搜索法[13]、核校準[14]等方法進行調參，本文針對上述兩個參數，構建二維坐標(σ2，Γ)，采用PSO進行自適應調參。

1.2 PSO-LSSVR近似模型構建方法

PSO算法是一種全局優化的進化算法[15]，因其具有相對較高的求解效率和較快的收斂速度，被廣泛用于模型調參、函數優化等領域。

PSO算法的原理是采用一種僅有速度和位置兩個屬性的粒子來模擬鳥類覓食的過程，其更新速度和位置的公式如下：

vi+1=ωvi+c1(pbest,i-xi)+c2(gbest,i-xi)

(6)

xi+1=xi+vi+1

(7)

式中：i為迭代次數；vi為粒子第i次迭代的速度向量；xi是粒子第i次迭代的位置向量；pbest,i是粒子歷史最優解的位置向量；gbest,i是粒子群當前最優解的位置向量；ω是慣性系數，用于調整全局尋優與局部尋優能力，一般取0至1之間；c1和c2被分別稱為每個粒子的個體學習因子和社會學習因子，通常取c1=c2=2。

均方根誤差(RMSE)因其量綱與數據集相同，更加便于解釋誤差，因此本文使用驗證集的RMSE作為PSO算法的適應度函數。RMSE的計算表達式如下：

(8)

圖1 PSO-LSSVR近似模型構建方法

1.3 自適應PSO-LSSVR近似模型構建方法

在自適應近似模型構建中，未達到收斂標準之前，每次迭代產生的解將會被重新放入數據集中，這些解被稱為校正點。校正點的選取對模型擬合效果有很大影響，目前較為常見的思路是響應最優策略[9]，即將每次迭代的最優解作為校正點，但Jones等[16]指出響應最優策略容易使得近似模型陷入局部最優解。

另一方面，在利用PSO算法進行LSSVR模型的調參時，將會根據驗證集數據RMSE選擇最佳超參數。如果驗證集是每次迭代后的數據集最優解時，理論上可以獲得在數據集最優解處擬合精度最高的近似模型，從而加快尋優速度；但如果數據集的最優解不是全局最優解，則更容易陷入局部最優。

綜上，為盡量避免自適應PSO-LSSVR近似模型陷入局部最優解，本文提出一種自適應近似模型構建方法。借鑒逐步縮減求解域[17]的方法，將近似模型分為全局自適應近似模型和局部自適應近似模型兩部分，其中全局自適應近似模型的主要目的是較為精確地獲得全局最優解的可能位置和數值，為局部自適應近似模型提供邊界條件；局部自適應近似模型利用全局自適應近似模型迭代的所有校正點，在全局最優解可能的位置處構建新的自適應近似模型，從而增加全局最優解附近的擬合精度。這樣，由于不直接將全局自適應近似模型的最終迭代解作為全局最優解，也就減小了陷入局部最優解的可能，而構建局部自適應近似模型的樣本點主要來自以往的校正點，從而加強了樣本點的利用率。

同時，為了增強自適應近似模型在最優點附近的搜索能力，自適應近似模型的驗證集將從上次迭代預測的最優解實際值和數據集最優解中隨機選擇，為了保證收斂結果有意義，同一驗證集選擇方案必須至少保持兩次迭代步。全局和局部校正點的收斂標準如下[18]:

(9)

(10)

1.4 優化問題求解策略

自適應PSO-LSSVR近似模型針對單目標優化問題的求解策略如圖2所示。

圖2 自適應PSO-LSSVR近似模型構建方法

1.5 自適應PSO-LSSVR近似模型尋優算例

1.5.1 測試函數

為驗證自適應PSO-LSSVR近似模型的有效性，使用Branin函數[19]作為單目標優化問題測試函數，并使用自適應PSO-LSSVR近似模型進行全局尋優。

Branin函數表達式如下：

(11)

其中，x1∈[-5,10]，x2∈[0,15]，目標函數為minfBF(x1,x2)。該函數理論最優解約為0.397 887，具有3個最優點。

初始全局數據集選取20組樣本點，為兼顧訓練和驗證效果，初始驗證集為數據集中函數值最小的兩組樣本點。MIGA算法個體數為6 000。

1.5.2 Branin函數尋優過程

對初始數據集構建自適應PSO-LSSVR近似模型。取ε1=5%，ε2=1%，N=20。迭代13次后，全局自適應PSO-LSSVR近似模型收斂，進入局部近似模型階段，迭代3次后，局部自適應PSO-LSSVR近似模型收斂。

自適應PSO-LSSVR近似模型在Branin函數上的校正點迭代過程如圖3所示。

圖3 Branin函數自適應尋優過程

1.5.3 與其他常用近似模型尋優的對比

為進一步說明自適應PSO-LSSVR近似模型相對于其他常用近似模型的優勢，在Branin函數的基礎上，分別使用與自適應PSO-LSSVR近似模型收斂后相等的樣本點數量，采用iSight軟件構建靜態二階近似模型(2階RSM)、三階RSM近似模型(3階RSM)、靜態Kriging近似模型和靜態RBF神經網絡近似模型(靜態RBF NN)，并使用相同參數的MIGA算法進行全局尋優，從而對比得到的最優值及其精度差異。各近似模型在Branin函數上的尋優結果如表1所示。

表1 各近似模型在Branin函數上的尋優結果1)Table 1 Optimization results in Branin function with different approximation models

綜上所述，自適應PSO-LSSVR近似模型在同等樣本量情況下可以獲得更優秀的尋優能力及預測效果，證明其較適合由于樣本量獲得較為困難而產生的小樣本優化問題。本文將應用自適應PSO-LSSVR近似模型對GTS模型低風阻尾板進行優化設計。

2 GTS模型低風阻尾板優化

2.1 GTS模型的計算設置及對標過程

本文采用GTS模型作為車輛研究對象。GTS模型是Sandia National Laboratories Ground Transportation System的基礎車輛模型[20]，車輛外形是一種具有貨箱的半掛汽車列車，車輛長19.805 m、寬2.590 m、高4.113 m。Gutierrez等[20- 21]對1:8 GTS模型進行了風洞試驗，試驗結果顯示，非側風工況時，車輛尾部流場呈現出幾乎對稱分布的流場結構。由于本文僅計算非側風工況穩態情況，為減少計算量，可認為流場左右對稱分布，而只計算半車身。

計算域方案參照SAE J2966-2017標準[22]，為減少計算域寬度對流場的影響，在推薦值基礎上增加至15倍半車寬，計算域全長220.0 m(前3倍車長、后7倍車長)、寬19.6 m(15倍半車寬)、高24.0 m(6倍車高)，阻塞比為1.06%。設置3層加密區，按照1:2增長率逐層加密，并對車尾分離區進行了再次加密。多面體網格具有生成與計算效率高、對梯度及局部流動分布預測性好的特點[23]，因此本文使用多面體網格進行劃分，網格劃分如圖4所示。

圖4 GTS模型網格劃分

數值計算使用STAR-CCM+軟件。采用RANS穩態計算方式和SSTk-ω湍流模型進行數值計算,整車邊界層厚度為2.16 mm，劃分為15層，使得全車絕大部分滿足y+≤1的要求，而僅在駕駛室側部邊緣處，由于氣流加速使得此處y+～4，為縮小可能帶來的計算誤差，使用全y+壁面處理功能[1]。殘差收斂標準為10-4，氣動阻力系數CD收斂后取后300步CD均值作為最終CD計算值。各邊界條件如表2所示。

表2 GTS模型仿真邊界條件Table 2 Boundary condition of simulation of GTS model

Gutierrez等[20- 21]的試驗雷諾數(以車寬作為特征長度)為1.6×106，滿足商用車縮比模型風洞試驗雷諾數應大于0.5×106的要求[24]，可認為進入雷諾數自準區，氣動阻力系數不再改變，這與Storms等[25]的試驗結果變化趨勢一致,因此可用于氣動阻力系數對標過程。網格無關性結果與試驗的對標結果如表3所示。

表3 氣動阻力系數試驗對標結果與網格無關性Table 3 Benchmark of CD and grid independency

1)ΔCD指CD的相對變化量。

由表3可以發現，當半車身體網格數為857萬時，氣動阻力系數幾乎不變，且與試驗值的誤差小于5%，而進一步增加體網格數使得計算量增長巨大，因此使用857萬網格方案進行后續計算。

2.2 GTS模型低風阻尾板的自適應優化過程

對GTS模型增加尾板，并通過自適應PSO-LSSVR近似模型研究最優尾板設計。初始方案尾板形狀及幾何參數如圖5所示。

圖5 GTS模型尾板設計變量及初始值

使用相同的仿真策略，得到增加尾板后GTS模型的氣動阻力系數為0.312 3，相比無尾板GTS模型下降15.14%。確定上尾板傾角θ1、下尾板傾角θ2、側尾板傾角θ3和尾板長度L作為設計變量。由于改型前后車輛正投影面積不變，故可將氣動阻力系數作為目標函數。

根據后擾流器的設計準則[6]，尾板的長度應滿足下式：

(12)

式中，A為商用車的正投影面積，γ′為比例系數，通常取0.28～0.57。計算可得尾板的長度范圍為885～1 805 mm。然而在一定情況下，尾板過長不利于側風穩定，因此選擇長度約束條件為500～1 000 mm，尾板各角度的約束條件均為5°～20°。整個優化問題描述可以表達為

minf(θ1,θ2,θ3,L)

(13)

s.t.θ1,θ2,θ3∈[5,20],

L∈[500,1000]。

使用最優拉丁超立方抽樣方法對樣本空間進行抽樣。初始數據集為22組樣本點，由于數據集因變量CD極差只有0.06左右，與各自變量極差相差過大，因此對數據集進行歸一化以提升LSSVR算法訓練效率[26]。歸一化算法如下：

(14)

其中，X=(θ1,θ2,θ3,L,CD)。

對初始數據集構建自適應PSO-LSSVR近似模型。為更加貼近實際優化需求，使用反歸一化的氣動阻力系數控制收斂過程，并取ε1=1%，ε2=0.5%，N=20。初始化近似模型驗證集為2組樣本點。

迭代6次后全局自適應PSO-LSSVR近似模型收斂，最優結果范圍縮小至θ1，norm∈[0.191,0.572]，θ2，norm∈[0.030,0.296]，θ3，norm∈[0.020,0.273]，Lnorm∈[0.939,1.000]，下標norm表示歸一化后。迭代3次后局部自適應PSO-LSSVR近似模型收斂。迭代過程如圖6所示。

圖6 氣動阻力系數自適應尋優過程

最優氣動阻力系數為CD,norm=-0.018 4,位于θ1,norm=0.523 3、θ2,norm=0.278 4、θ3,norm=0.216 7、Lnorm=0.976 6處。反歸一化結果為CD=0.282 5，位于θ1=12.85°、θ2=9.18°、θ3=8.25°、L=988.31 mm處。計算得到該方案的氣動阻力系數為0.283 0，自適應近似模型預測誤差僅為0.18%，證明了構建的自適應PSO-LSSVR近似模型可以用于GTS模型低氣動阻力尾板的優化。

2.3 最優低風阻尾板減阻機理

相對于初始尾板工況，最優尾板工況氣動阻力系數下降達9.38%。由圖7、圖8壓力系數分布可知，相比于初始方案，優化后方案在尾板附近流場壓力出現回升，由圖9可知，優化后車尾表面壓力系數明顯提高；由圖10可以發現，優化后方案的總壓為零等值面縮小，說明分離區減小，從而降低了壓差阻力。觀察圖7、圖8的速度矢量分布，優化后的尾板使得C柱下洗氣流強度減弱，使得上洗氣流得到充分發展，上洗渦(B渦)更靠近車尾，下洗渦(A渦)和上洗渦(B渦)呈現出對稱分布的特點，駐點上移，使得尾部分離區面積減小，側板附近的分離渦(C渦)尺度減小。優化后方案的對稱面尾流場分布，與Hassaan等[4]對無車輪GTS模型尾板研究中相似板長最低阻力工況的尾流場分布類似。圖11是使用λ2渦識別方法得到的尾流分離區渦分布圖，可以看出，由于上洗氣流強度的減弱，優化后流向渦尺度明顯減小，且位置不再緊貼地面。

圖7 對稱面速度矢量及壓力系數分布

圖8 Z=2.3 m處的速度矢量及壓力系數分布

圖9 車尾表面壓力系數分布

圖10 總壓為零等值面

3 結論

(1)通過合理選擇初始數據集樣本量和驗證集，基于PSO-LSSVR機器學習算法的自適應近似模型在數據集樣本量一定的情況下，于 Branin函數上的尋優能力和擬合精度超過了常用的4種近似模型，說明自適應PSO-LSSVR近似模型優化方法具有較大的快速尋優潛力。

圖11 λ2=-5等值面

(2)自適應PSO-LSSVR近似模型可以通過對最優解構建局部近似模型的方法,增強模型在最優解處的擬合精度。

(3)將自適應PSO-LSSVR近似模型優化方法應用在GTS模型低風阻尾板優化上，僅使用31組樣本便得出了最優方案，最優方案相比初始方案減阻達9.38%，且局部近似模型擬合誤差僅為0.18%。由此說明自適應PSO-LSSVR近似模型優化方法在保證尋優精度的情況下，尋優效率提升明顯，適用于單樣本點獲取耗時較多的小樣本單目標優化問題。