基于改進NSGA-Ⅱ的汽油機標定優(yōu)化研究

康順，李自勝，肖曉萍，陳杰，張楷，尚偉

(1.西南科技大學制造科學與工程學院，四川 綿陽 621010；2.西南科技大學工程技術中心，四川 綿陽 621010；3.綿陽新晨動力機械有限公司，四川 綿陽 621000；4.西南交通大學機械工程學院，四川 成都 610031)

為使發(fā)動機在經(jīng)濟性、動力性及排放性上均處于最佳狀態(tài)，需要各控制參數(shù)在不同工況下做出相應的改變。發(fā)動機標定的目的是找到合適的控制參數(shù)組合，最大限度地提高發(fā)動機性能與降低污染物排放[1]。發(fā)動機電控技術的發(fā)展極大地提高了發(fā)動機的經(jīng)濟性和排放性，但也導致了發(fā)動機控制參數(shù)的增加及參數(shù)間耦合性增強，發(fā)動機的標定工作也變得更加困難。傳統(tǒng)的人工標定所需周期長，試驗量大，已不能滿足現(xiàn)代發(fā)動機的標定要求。為了提高標定效率，基于模型的標定方法被提出并成為近年來國內(nèi)外學者研究的課題[2-3]。

基于模型的標定方法需要依據(jù)試驗設計[4]獲得有限數(shù)據(jù)，建立預測模型，進而優(yōu)化發(fā)動機全工況下控制參數(shù)，工況點優(yōu)化后，將不同工況點的參數(shù)組合在一起，同發(fā)動機轉(zhuǎn)速和負載一起寫入校準圖(MAP)中?？刂茀?shù)的優(yōu)化可分為局部優(yōu)化和全局優(yōu)化。局部優(yōu)化是通過預測模型對單個工況點進行優(yōu)化，獨立搜索控制參數(shù)的最佳組合[5]，這是目前運用最廣泛的方法。全局優(yōu)化則全面考慮所有運行工況點，通過施加大量約束條件，并綜合考慮不同的量級折中[6]，對控制參數(shù)進行整體優(yōu)化。全局優(yōu)化是最全面的方法，但由于其局限性，并未得到廣泛應用。

近年來，國內(nèi)外學者對發(fā)動機預測模型與優(yōu)化方法進行了研究并取得一系列的成果。Francesco等[7]、Martínez等[8]使用少量數(shù)據(jù)建立人工神經(jīng)網(wǎng)絡(ANN)模型并完成了發(fā)動機標定工作，但所建立的模型存在過擬合的風險。柴志剛等[9]基于克里格法(Kriging)生成模型預測發(fā)動機性能、排放等。張文[10]基于高斯過程建立油耗預測模型，并將優(yōu)化問題轉(zhuǎn)化為約束規(guī)劃問題，實現(xiàn)了VVT參數(shù)的標定，但該研究僅針對油耗進行優(yōu)化，并未考慮其他指標。Federico等[11]應用遺傳算法(Genetic Algorithms，GA)快速有效地獲得了更好的標定，但在局部優(yōu)化后，仍需要人工決策生成校準圖。Wong等[12]使用不同機器學習算法，包括極限學習機(Extreme Learning Machine，ELM)、最小二乘支持向量機(Least-Squares Support Vector Machine，LS-SVM)及徑向基神經(jīng)網(wǎng)絡(RBFNN)創(chuàng)建數(shù)據(jù)驅(qū)動模型，并分別采用模擬退火(Simulated Annealing，SA)與粒子群算法(Particle Swarm Optimization，PSO)優(yōu)化生物柴油發(fā)動機，試驗表明：ELM在模型精度與訓練時間上優(yōu)于LS-SVM和RBFNN；PSO在適應度方面優(yōu)于SA，但在局部優(yōu)化中，其目標函數(shù)是各性能指標的加權(quán)和，各目標加權(quán)值的分配帶有較大的主觀性。張文[13]采用NSGA-Ⅱ?qū)蝹€工況點進行局部優(yōu)化，并采用人工插值法完成了VVT參數(shù)的標定，但求得的Pareto解分布并不均勻。

相較于人工標定方法，基于模型的標定方法優(yōu)勢明顯，但也存在一些缺陷，如模型過擬合，優(yōu)化目標單一及局部優(yōu)化后需要人工決策生成校準圖等。針對該問題，提出一種局部優(yōu)化與校準圖平滑步驟相結(jié)合的方法，建立基于高斯過程回歸的發(fā)動機性能預測模型，改進NSGA-Ⅱ算法并用于發(fā)動機控制參數(shù)優(yōu)化，自動生成平滑校準圖獲得標定結(jié)果。基于某3缸汽油機的臺架數(shù)據(jù)，利用GT-Suite仿真環(huán)境驗證了該方法的有效性。

1 發(fā)動機性能模型與工況點優(yōu)化

1.1 發(fā)動機性能預測模型

1.1.1 高斯過程回歸

高斯過程回歸[14](Gaussian Process Regression，GPR)是由貝葉斯框架與統(tǒng)計學習方法相結(jié)合的非參數(shù)模型，它使用極大似然估計進行訓練，并且其性質(zhì)完全由其輸入的均值函數(shù)m和協(xié)方差函數(shù)k決定。對于樣本集合M={(xi,yi)|i=1,2,…n}(其中xi∈RD為輸入向量，yi∈R為輸出變量)，構(gòu)建回歸模型：

y=f(x)+ε。

(1)

(2)

以及觀測值y和預測值f*的聯(lián)合先驗分布：

(3)

式中：K(X,X)為訓練集的協(xié)方差矩陣；K(x*,X)為預測點x*與訓練集X之間的協(xié)方差矩陣；K(X,x*)=K(x*,X)T；K(x*,x*)為預測點自身的協(xié)方差；In為n維單位矩陣。

根據(jù)多元高斯聯(lián)合分布可計算出預測值f*的后驗分布：

(4)

即

(5)

核函數(shù)是GPR的核心，本研究選用具有平滑特性的平方指數(shù)函數(shù)[15](Squared Exponential，SE)作為核函數(shù)，其具體形式如下：

(6)

式中：l為方差尺度；σf為信號方差。參數(shù)集合θ={l,σf,σn}即為超參數(shù)，其中σn為噪聲方差，一般采用極大似然法求得。首先建立訓練樣本條件概率，得負對數(shù)似然函數(shù)L(θ)=-logp(y|X,θ)，并對超參數(shù)求偏導；然后采用共軛梯度法對偏導數(shù)進行最小化得到超參數(shù)得最優(yōu)解。負對數(shù)似然函數(shù)及其對超參數(shù)的偏導數(shù)具體形式如下：

(7)

式中：C=Kn+σnIn;α=(K+σnIn)-1y。

1.1.2 模型變量與數(shù)據(jù)處理

將燃油消耗率(BSFC)、排放顆粒數(shù)量(PN)、NOx、CO及HC排放作為模型的輸出，將表征工況的發(fā)動機轉(zhuǎn)速、扭矩與控制參數(shù)進氣VVT、排氣VVT及點火角共同作為模型的輸入。

在高斯過程回歸模型建立過程中，輸入?yún)?shù)與輸出數(shù)值間差異較大，因此使用Z-score標準化方式對數(shù)據(jù)進行歸一化處理，如式(8)。

(8)

1.1.3 模型建立

高斯過程回歸預測建模步驟如下。

1)對試驗數(shù)據(jù)進行歸一化處理，并將數(shù)據(jù)集劃分為訓練集和測試集，其中訓練集用于訓練模型，測試集用于檢驗模型最終性能。

2)對GPR模型核函數(shù)的超參數(shù)進行優(yōu)化，一般采用極大似然估計方法確定。首先確定訓練樣本的負對數(shù)似然函數(shù)，然后采用共軛梯度法對似然函數(shù)偏導數(shù)進行最小化以得到超參數(shù)最優(yōu)解，具體操作見第1.1.1節(jié)。

3)按照如圖1所示模型結(jié)構(gòu)，利用高斯過程回歸模型對測試樣本進行預測。需要注意的是，并非只需要建一個模型，而是需要建5個模型，分別用于燃油消耗率(BSFC)、CO等5個性能及排放參數(shù)的預測。

圖1 GPR模型結(jié)構(gòu)

1.1.4 評價指標

為了驗證模型的預測能力和準確性，引入了3個評價指標[16]。

1)決定系數(shù)(R2)：

(9)

式中：yi,model為預測值；yi,truth為真實值(試驗值)；ymean為真實值的平均值。R2越大，模型精度越高。

2)相對平均絕對誤差(RAAE)：

(10)

RAAE的值越小，模型精度越高。

3)相對絕對最大誤差(RMAE)：

(11)

RMAE表示局部最大誤差，RMAE的值越小，模型精度越高。

1.2 基于約束違反的改進NSGA-Ⅱ算法

發(fā)動機標定優(yōu)化是多目標優(yōu)化問題，需要在不同的目標之間進行折中。本研究對NSGA-Ⅱ進行改進并應用于發(fā)動機標定優(yōu)化中。

1.2.1 基于高斯過程回歸的NSGA-Ⅱ算法

基于精英策略的快速非支配排序遺傳算法[17](Non-Dominated Sorting Genetic Algorithm-Ⅱ，NSGA-Ⅱ)是一種多目標優(yōu)化算法，NSGA-Ⅱ在運行時需要多次評估目標值，自變量會在可行解范圍內(nèi)任意取值，因此需要事先對每個目標函數(shù)與約束條件建立預測模型，以模型代替試驗臺架，在算法運行過程中提供目標值。本研究優(yōu)化目標是針對單個工況點，使其在其他排放指標(PN、CO及CH排放)的約束條件下盡可能地減少BSFC和NOx排放。將工況參數(shù)與需要標定的控制參數(shù)轉(zhuǎn)化為適合NSGA-Ⅱ的數(shù)據(jù)輸入，這里采用實數(shù)編碼，編碼長度為5。如圖2，每個解都稱為染色體(個體)，每個控制參數(shù)與工況參數(shù)是一個基因，其中控制參數(shù)為等位基因，可進行交叉變異操作，工況參數(shù)為固定基因，不進行任何操作。

圖2 以基因表示工況參數(shù)與控制參數(shù)

NSGA-Ⅱ算法主要環(huán)節(jié)如下。

1)初始化種群：設置工況參數(shù)不變，控制參數(shù)在可行范圍內(nèi)生成初始種群，個體如圖2所示，種群大小為200。

2)父代種群經(jīng)過選擇、交叉和變異生成子代種群，并將兩個種群合并為一個完整種群。

3)計算種群中的個體目標值及約束值，對個體進行非支配排序和擁擠度計算。其中個體目標值和約束值是通過GPR模型預測得到。

4)選擇Pareto等級更高且擁擠度更大的個體組成下一代父代種群。

5)重復步驟2至步驟4，直至到達最大迭代數(shù)。最大迭代數(shù)設置為300。

通過在步驟2交叉操作中改進交叉算子及在步驟3中加入約束違反來對NSGA-Ⅱ進行改進。

1.2.2 改進NSGA-Ⅱ算法

針對NSGA-Ⅱ在發(fā)動機標定優(yōu)化問題上解的分布性和收斂性較差問題，提出一種改進的NSGA-Ⅱ算法，在非支配排序中加入約束違反并引入混合交叉算子。

1)交叉算子的改進

NSGA-Ⅱ算法采用SBX算子，其搜索性能較差。針對這一缺陷，結(jié)合正態(tài)分布交叉(NDX)算子[18]，提出一種混合交叉算子(NDX-SBX)，在迭代初期NDX算子占較大權(quán)重，從而跳出局部最優(yōu)解，在算法后期SBX算子占較大權(quán)重，能在小范圍內(nèi)搜索。(NDX-SBX)公式如式(12)所示。

(12)

2)添加約束違反

NSGA-Ⅱ在判斷解的支配關系時，通常需要比較適應度與擁擠度。在解決多約束問題時，可采用約束違反來影響解的質(zhì)量，從而保證解的多樣性。在判斷解X1優(yōu)于解X2時，需滿足以下條件：

(1)解X1為可行解，解X2為不可行解；

(2)解X1與解X2同為不可行解，但解X1約束違反程度小于解X2；

(3)解X1與解X2同為可行解，但解X1支配解X2。

1.3 校準圖自動平滑

工況點標定優(yōu)化后，需要人工決策生成校準圖。針對這一缺陷，提出一種平滑校準圖自動生成的方法。這里主要針對WLTC循環(huán)工況[19]，使用基于時間權(quán)重的參考點和關鍵工況點(Key Point，KP)的標定優(yōu)化結(jié)果，生成多張校準圖，然后評估這些校準圖的性能與平滑度，自動選取一些優(yōu)秀的校準圖。

1.3.1 基于時間權(quán)重的工況點計算

使用k-近鄰聚類算法計算每個KP的時間權(quán)重，如圖3所示，圖中聚類了20個KP，每一工況點都與最近的KP相關聯(lián)，每個KP用相應的大符號表示，與之相關聯(lián)的工況點用相同顏色的小符號表示。使用3個權(quán)重最大的KP生成控制參數(shù)平面，分別是30 N·m@1 500 r/min，60 N·m@2 000 r/min，90 N·m@2 250 r/min。通過該平面以及剩余的17個KPs生成校準圖。

圖3 WLTC運行工況聚類點

1.3.2 校準圖創(chuàng)建

單個工況點的標定優(yōu)化結(jié)果是一組Pareto解，一組Pareto解中包含多個解。在3個權(quán)重最大的KP中，分別選取一個解，并使用三點法生成不同控制參數(shù)的平面。不同控制參數(shù)的平面即關于發(fā)動機轉(zhuǎn)速和扭矩的函數(shù)，如式(13)：

zi=a×s+b×Tor+c。

(13)

式中：zi為第i個控制參數(shù)的值；s為發(fā)動機轉(zhuǎn)速；Tor為扭矩；a,b,c為常量。點火角參數(shù)示例見圖4。

圖4 點火角平面

平面確定之后，使用該平面與剩余的17個KP的每一個解，如式(14)，選擇與該平面映射點最接近的解，使得曲面變形最小。

ξz,i=min|pi-pi,ref|。

(14)

式中：pi為第i個工況點的解；pi,ref為第i個工況點映射到校準平面上的點；ξz,i為第i個工況點的解到映射點的距離。對點火角參數(shù)的校準見圖5。

圖5 點火角校準圖

但發(fā)動機標定不只標定點火角參數(shù)，需要標定許多控制參數(shù)，如果點火角校準圖足夠平滑，但其他校準圖卻不平滑，則此次標定是失敗的。所以，需要考慮所有控制參數(shù)，如式(15)所示：

(15)

式中：wi表示第i個控制參數(shù)的權(quán)重；nor(ξz,i)表示歸一化偏差距離。由于每個控制參數(shù)的偏差不一樣，所以需要進行歸一化處理。對控制參數(shù)權(quán)重做出以下設置：進氣VVT角、排氣VVT角及基礎點火角的權(quán)重分別為0.3，0.3及0.4。

1.3.3 平滑度評價指標

將校準圖視作規(guī)則曲面，使用高斯曲率評估校準圖的平滑度。高斯曲率[20]為曲面主曲率k1,k2的乘積，對于曲面上任意一點(x,y)，都有其高斯曲率：

(16)

式中，高斯曲率使用一階和二階微分計算曲面的平滑度。校準圖的總體平滑度用校準圖中網(wǎng)格點的高斯曲率之和來表示，如式(17)：

(17)

式中：Kxy為校準圖上單個網(wǎng)格點的高斯曲率。

2 試驗與分析

為驗證所提方法的有效性，對所建預測模型的預測結(jié)果進行分析，將改進的NSGA-Ⅱ算法與原NSGA-Ⅱ算法進行比較，并將改進的NSGA-Ⅱ算法與平滑校準圖自動生成方法結(jié)合，生成標定結(jié)果，比較該結(jié)果與人工標定結(jié)果的性能與平滑度。試驗數(shù)據(jù)來源于某直列3缸汽油機，發(fā)動機的基本參數(shù)見表1。發(fā)動機所用預測模型為所建立的高斯過程回歸模型。模型使用600組訓練數(shù)據(jù)進行訓練，120組測試數(shù)據(jù)進行測試。

表1 發(fā)動機技術參數(shù)

2.1 預測模型及預測結(jié)果分析

表2列出用于建立性能及排放預測模型的部分原始數(shù)據(jù)，模型輸入為轉(zhuǎn)速、扭矩、進氣VVT角、排氣VVT角和點火角，輸出分別為燃油消耗率、NOx排放、排放顆粒數(shù)量(PN)、CO排放及HC排放。表3列出各GPR模型訓練完成后的最優(yōu)超參數(shù)。圖6示出測試數(shù)據(jù)在各模型的擬合效果，展示了所有預測點的發(fā)動機燃油消耗率、PN排放、NOx排放、CO排放及HC排放試驗結(jié)果與預測結(jié)果的比較。從圖6中預測點的分布可以看出，性能與排放參數(shù)的預測值與試驗值具有很好的相關性。表4示出預測結(jié)果的多目標評價。從表4可以看出，GPR模型對5個參數(shù)的預測R2均大于0.96，表明這些模型具有良好的預測能力并能夠用于后續(xù)多目標優(yōu)化。

表3 訓練完成后各模型超參數(shù)

圖6 發(fā)動機性能與排放預測結(jié)果與試驗結(jié)果比較

表4 模型多指標評價

2.2 改進算法收斂性和穩(wěn)定性分析

為驗證改進算法的有效性，采用兩個評價指標來評估算法性能。

1)收斂性評價指標γ，即解集X中每個點到參考集X*的平均最小距離，γ越小，則算法的收斂性越好。

(18)

式中：X為算法獲得的解集；X*為參考集；dis(x,y)為解集與參考集中點的距離。

2)多樣性評價指標Δ，表示解集分布在前沿面上的均勻程度。

(19)

采用公開測試函數(shù)(ZDT1,ZDT2,ZDT3,ZDT4)進行測試，其中種群大小設置為100，最大迭代次數(shù)設置為250，其交叉算子γ1和變異算子γ2分別設置為2和5。

對每一基準函數(shù)測試20次，對每一指標取平均值(見表5)。兩種算法在ZDT1與ZDT4測試函數(shù)對比分別如圖7與圖8所示。從表5可以看出，ENSGA-Ⅱ在保證收斂性的同時能夠提高分布性，但在基準函數(shù)ZDT3中表現(xiàn)欠佳。總體來看，ENSGA-Ⅱ收斂性和分布性都有所提高。

表5 兩種算法基準函數(shù)測試結(jié)果對比

圖7 兩種算法對ZDT1函數(shù)測試結(jié)果對比

圖8 兩種算法對ZDT4函數(shù)測試結(jié)果對比

2.3 發(fā)動機標定優(yōu)化實驗與結(jié)果分析

在NSGA-Ⅱ算法中，非支配排序和擁擠度計算都是基于目標函數(shù)的結(jié)果，因此需將GPR模型嵌入到其中，將部分GPR模型作為目標函數(shù)，部分GPR模型作為約束條件，進行多目標優(yōu)化。優(yōu)化目標是針對單個工況點，使其在其他排放指標的約束條件下盡可能地減少燃油消耗量和NOx。目標函數(shù)定義如式(20)：

(20)

式中：wnwe為進氣VVT參數(shù)；wnwa為排氣VVT參數(shù)；spark為點火角參數(shù)。定義約束條件如式(21)：

(21)

以90 N·m@2 250 r/min工況點為例，闡述試驗步驟：

1)初始化種群。固定工況參數(shù)轉(zhuǎn)速和扭矩不變，隨機生成如[2 250,90,40,25,20]在可行解范圍內(nèi)個體，個體中5個參數(shù)分別代表轉(zhuǎn)速、扭矩、進氣VVT、排氣VVT及點火角。其中種群數(shù)量為200。

2)選擇、交叉和變異。選擇策略采用二元錦標賽機制；交叉策略采用第1.2.2節(jié)中改進的交叉策略，交叉算子γ1=2，變異算子γ2=5。

3)計算種群中的個體目標值及約束值，然后對個體進行非支配排序和擁擠度計算。其目標函數(shù)即為BSFC模型(基于高斯過程回歸建立BSFC模型)和NOx模型，約束條件為CO模型、HC模型及PN模型。

4)選擇Pareto等級更高且擁擠度更大的個體組成下一代父代種群。

5)重復步驟2至步驟4，直至到達最大迭代數(shù)。設置最大迭代數(shù)為300。

GPR模型聯(lián)合改進的NSGA-Ⅱ優(yōu)化算法流程見圖9。

圖9 GPR聯(lián)合ENSGA-Ⅱ算法流程圖

在CPU主頻為2.9 GHz、內(nèi)存為8 GB的計算機，Pycharm2022平臺上用改進NSGA-Ⅱ進行發(fā)動機標定優(yōu)化，設置種群數(shù)量為200，最大迭代次數(shù)為300，設置交叉算子γ1=2，變異算子γ2=5，優(yōu)化一個工況點大概需要30 min。由于工況點數(shù)較多，這里僅展示個別工況點優(yōu)化結(jié)果，圖10與圖11分別示出90 N·m@2 250 r/min與60 N·m@2 000 r/min工況點標定優(yōu)化結(jié)果。從圖10與圖11中可以看出，Pareto解的分布性較好且生成了多組可行解(發(fā)動機控制參數(shù))，表明ENSGA-Ⅱ算法在發(fā)動機標定優(yōu)化中具有較好的效果。

圖10 90 N·m@2 250 r/min工況點標定優(yōu)化結(jié)果

圖11 60 N·m@2 000 r/min工況點標定優(yōu)化結(jié)果

2.4 平滑度校準結(jié)果與比較

將關鍵工況點標定優(yōu)化結(jié)果與平滑校準圖自動生成的方法結(jié)合，標定結(jié)果見圖12至圖14，圖12至圖14分別為進氣VVT角，排氣VVT角以及點火角校準圖。

圖12 進氣VVT角校準圖

圖13 排氣VVT角校準圖

將平滑度與性能相對均衡的校準圖與人工校準圖相對比，分析校準圖平滑度及其性能。圖15示出進氣VVT，排氣VVT以及基礎點火角校準圖與人工標定的校準圖平滑度之間的比較。

圖15 平滑度比較

從圖15可以看出，算法校準的校準圖比人工校準的校準圖更加平滑。為了評估校準圖在發(fā)動機上的性能，利用建模軟件GT-Suite建立車輛模型。傳動系模塊中的駕駛員模型通過PID控制根據(jù)期望和實際車速計算所需扭矩，然后調(diào)整制動踏板開度，以匹配模擬行駛車速，其中測試循環(huán)設置為WLTC。將人工標定的校準圖與算法標定的校準圖數(shù)據(jù)分別寫入發(fā)動機MAP中。測試表明，在WLTC循環(huán)測試工況下，算法標定的校準圖與人工標定的校準圖相比，燃油消耗率降低1%，氮氧化物NOx排放降低了5%。

3 結(jié)束語

建立并使用GPR模型對燃油消耗率、NOx等性能和排放參數(shù)進行預測，并將GPR模型作為改進NSGA-Ⅱ算法的目標函數(shù)，進而求解發(fā)動機局部優(yōu)化問題。將校準圖平滑步驟與發(fā)動機局部優(yōu)化結(jié)合在一起，生成多張校準圖。通過采集某3缸汽油機的臺架試驗數(shù)據(jù)，并使用所提出的方法進行標定，試驗結(jié)果表明，所建立的GPR模型R2均大于0.96，局部優(yōu)化后的Pareto解更加均勻，生成的校準圖在WLTC循環(huán)工況下測試，經(jīng)濟性和排放性比人工標定的更好。該方法可為減少發(fā)動機試驗工作量提供參考。