管湘源,儲江偉,高偉健,趙小婷
(東北林業大學 交通學院,哈爾濱 150040)
發動機是車輛的動力系統,是一個復雜的非線性、時變性及滯后系統[1],使用真實發動機來設計控制器是一項費時費力的任務。因此,許多研究人員通過建立發動機數學模型來設計、測試和優化其控制器。汽車行業已經轉向基于模型的校準來尋求解決方案[2]。當前主流發動機建模方法有兩種:一種是解析法,即依據發動機的機理和相應參數關系而建立仿真數學模型;第二種是實驗法,即通過發動機臺架試驗數據建立仿真模型[3-5]。
通過解析法建立數學模型,通常需要根據數據擬合(一般為曲面擬合法、多元插值法和綜合法)來確定相關系數(即曲線或曲面的待定系數),進而得到最終的數學擬合公式。利用Matlab擬合工具箱進行數據擬合,具有求解便捷和精度較高等特點,因而在機械、物理、電力、水利、醫藥和能源等多種學科領域均有廣泛應用。在機械領域應用方面,鄭嘯洲等[6]通過Matlab擬合工具箱對某型發動機試驗數據進行擬合處理,建立發動機的數學模型。在物理領域應用方面,陳文芳等[7]和魏奶萍[8]通過Matlab曲線擬合工具箱處理PN結的伏安特性數據,擬合得到相應的PN結伏安特性函數表達式;郭悅韶[9]通過Matlab曲線擬合工具箱處理霍爾效應實驗測量數據,擬合得到磁場分布曲線;范雋宏等[10]通過Matlab曲線擬合工具箱對軟磁材料的單位損耗數據的曲線擬合分析,得出用于變壓器電磁計算的擬合方程;岳鵬等[11]以測量玻璃熱膨脹系數和折射率溫度系數實驗、太陽能電池基本特性測量實驗為例,通過Matlab曲線擬合工具箱處理物理實驗數據,進行擬合實驗曲線。在電力領域應用方面,李太興等[12-13]通過Matlab曲線擬合工具箱進行火電廠運行的機組負荷優化分配;高艷平等[14]通過Matlab曲線擬合工具箱分析某發電廠的地基沉降實測資料,并擬合出多種不同的地基沉降預測模型。在水利領域應用方面,王杰等[15]通過Matlab曲線擬合工具箱求解分析抽水試驗資料,并確定承壓含水層的導水系數和導壓系數。在醫藥領域應用方面,曹俊涵等[16]通過Matlab曲線擬合工具箱處理藥物溶出數據,并計算藥物溶出度Weibull分布參數。在能源領域應用方面,高鵬[17]通過Matlab曲線擬合工具箱對浮選尾礦灰分數據與圖像灰度數據進行曲線擬合;何定武[18]通過Matlab曲線擬合工具箱分析風機特性曲線并計算風機實際工況下節電率。
在上述應用中,經過驗證,由Matlab曲線擬合工具箱得到的擬合結果具有較高精度,因此本研究也選用Matlab曲線擬合工具箱作為曲線擬合工具,以某款發動機的外特性和萬有特性曲線為例,通過Matlab曲線擬合工具箱cftool對提取的數據按照不同的擬合方法進行公式擬合,并根據擬合優度檢驗結果選擇出最佳擬合公式,此公式即為最終確定的外特性曲線的數學模型。
以發動機臺架試驗數據的離散點或萬有特性圖數字化提取的坐標點數據作為樣點,選擇一種合適的函數逼近方式,將它重構成確定的數學模型,這是建模的關鍵。
發動機數學建模的原理是基于函數逼近,即函數的近似表示問題。函數逼近的方法多種多樣,具體分類如圖1所示。

圖1 函數逼近方法分類Fig.1 Function approximation method classification
由圖1可知,函數逼近的常規方法為擬合法、插值法和綜合法。
曲面(線)擬合法是選定近似函數H(x)后,不要求H(x)通過已知樣本點,只要求在某種意義下H(x)在這些樣點的總偏值最小。采用多元線性回歸方法進行曲面擬合,構造發動機穩態輸出轉矩、油門開度以及發動機轉速的關系。發動機的負荷特性數學模型構建為[19]:
(1)
式中:Me為發動機輸出轉矩,Nm;α為發動機的油門開度,%;ne為發動機的轉速,r/min;B為系數矩陣。
通過最小二乘原理來確定系數矩陣B。建立的回歸方程為[20]:
(2)
式中:b,e為系數矩陣B的余數。
寫成矩陣形式如下:
Z=G·B+E。
(3)
G的列數k與多項式階數s的關系為:k=(s+1)(s+2)/2。其擬合值和測量值的逼近程度可以用擬合度來評價。

插值法是給出函數f(x)的一些樣點值,選定某些便于計算的函數,要求函數通過已知樣點,由此確定函數H(x)為f(x)近似值。
采用多元插值法,構造發動機穩態輸出轉矩、油門開度以及發動機轉速的關系。發動機的負荷特性數學模型構建為:
Me(t)=f(ne(t),αe(t))。
(4)
式中:Me(t)為發動機輸出轉矩,Nm;αe為發動機的油門開度,%;ne為發動機的轉速,rad/s。
要計算的點ne、αe落在4個網格點中間,即
(5)
定義:
(6)
則:
(7)
其中,cij是插值系數,插值問題的核心就是求出插值系數,插值算法不同所求出的插值系數各異。
Matlab中的曲線擬合等相應的工具箱或cubic、linear、neast、spline、interpn、interpft、ndgrid等函數可以用來實現插值的功能。下面對比分析幾種常用的插值方法的特點及用途,如圖2所示。

圖2 幾種常用的插值方法Fig. 2 Several interpolation methods
綜合法是融合多項式擬合與一元函數插值的方法。它首先根據最小二乘原理,擬合得到每個轉速下的T=f(α);然后對于某一工況點A,根據轉速判斷最臨近于工況點A的3條負載特性曲線,再用該曲線的多項式計算出來A點α對應的轉矩Ti-1、Ti、Ti+1的值,最后用插值法求得T(α)[23]。
傳統數據提取方法是進行手工測量和估算,此種方法不但誤差大,而且效率低。因此,有必要采取數字化方法提取萬有特性曲線圖的數據。選用圖像數字化軟件Get Data作為數據提取的工具,其提取數據流程如圖3所示。
提取數據前,需對圖像進行預處理,包括圖像的放大、截取和擺正等清晰化過程。良好的預處理可以在一定范圍內降低視覺誤差和分辨率誤差,從而使獲得的數據更為精確。
文中選取某款汽油機,根據廠方提供的萬有特性曲線(圖4)和外特性曲線(圖5)的圖像數據,對圖上的數據進行采集,提取為數據表形式,為后期建立數學模型奠定基礎。

圖3 Get Data提取數據流程圖Fig.3 Get Data extraction data flow chart
將萬有特性曲線(圖4)和外特性曲線(圖5)分別導入Get Data軟件,經過設置曲線顏色、設置背景色、設置基準點與取值范圍的幾個前期過程后,選擇曲線相應坐標點,進行對應數據的提取。萬有特性曲線的數據提取,如圖6所示。外特性曲線的數據提取,如圖7所示。

圖4 發動機萬有特性曲線Fig.4 Engine universal characteristic curve

圖5 發動機外特性曲線Fig.5 Engine external characteristic curve

圖6 萬有特性圖坐標點數據提取Fig. 6 Universal feature map coordinate point data extraction

圖7 外特性圖坐標點數據提取Fig.7 External characteristic map coordinate point data extraction
提取的所有坐標點的數據將顯示在工作區間內,可將數據導出為表格形式,見表1。

表1 導出的部分數據
除此之外,也可以選用Origin Pro軟件進行數據提取,其提取具體流程如圖8所示。
外特性建模采用曲線擬合的方法。首先初步確定外特性曲線的函數表達式(即確定近似函數),進而選擇相應的擬合方式得到最優擬合系數,最終得到完整的外特性數學模型。

圖8 Origin Pro提取數據流程圖Fig.8 Origin Pro extraction data flow
3.1.1 函數形式確定
發動機外特性是研究汽車動力性的基礎,它是在發動機節氣門全開的情況下,發動機轉矩Me(或功率Pe)、燃油消耗率be與發動機曲軸轉速ne之間的關系,可以看作是發動機轉速的一元函數[24]。
分別用x表示ne,y表示Me、Pe或be,n為擬合多項式的最高階數。
假定其函數形式為:
y=p1+p2x+p3x2+…+pnxn。
(8)
3.1.2 擬合方式選擇
本文以Matlab的曲線擬合工具箱cftool作為擬合工具。cftool中包含的擬合方法有:Exponential指數擬合、Fourier傅里葉擬合、Gaussian高斯法、Interpolant內插法、Polynomial多項式(1~9階)、Rational有理擬合、Power冪指數擬合、Smoothing spline平滑樣條擬合、Sum of Sin Functions正弦曲線擬合、Weibull兩個參數的Weibull分布,并支持Custom Equations用戶自定義。
本文選用多項式擬合方法,并分析擬合效果,在后續擬合優度檢驗一節中給出詳細說明。
發動機外特性模型按以下流程進行擬合:將發動機轉矩Me(或功率Pe)、燃油消耗率be和發動機轉速ne的數據分別導入到cftool工具箱中,發動機轉速ne設置為x,燃油消耗率be、發動機功率Pe和發動機轉矩Me分別依次設置為y。隨后設置權重,默認值為1。完成初始數據設定后,選擇擬合方法。擬合運算執行后,其結果在Results中顯示,如圖9所示。

圖9 擬合結果Fig.9 The results of fitting
在擬合過程中,選擇“Center and scale”(標準化),可減少出現病態矩陣的情況。采用擬合算法時應盡量選擇合適的最高冪次[25],但并非次數越高越好。因為當次數過高,容易出現龍格(Runge)現象,導致在試驗點以外的數據上偏差可能很大。
3.1.3 數學模型確定
通過多項式擬合得到外特性的最優擬合系數,最終得到外特性數學模型,見表2。

表2 外特性曲線的最佳擬合公式Tab.2 Best fit formula of external characteristic curve
發動機的燃油消耗率be是發動機轉速和有效轉矩(ne,Me)的函數,設發動機萬有特性的函數類型為多項式[26-29],其回歸模型的矩陣形式為[30-31]:
be=[1,ne,Me,ne2,neMe,Me2,…,neMel-1,Mel]×

(9)
式中:a0、a2、…、ak-1為模型中的待定系數;e0、e2、…、ek-1為隨機誤差;l為多項式的最高冪次;k為多項式的項數,k=(l+1)(l+2)/2[29-30]。
發動機萬有特性曲線的數學模型在cftool工具箱中的建立流程類似于外特性曲線的建立流程,在此不做贅述。
為評價擬合的性能,需對擬合結果進行擬合優度檢驗。具體的檢驗方法為:
(1)查看圖的擬合結果“Examine the Graphical Fit Results”。選擇Axes Limits,放大x軸坐標,顯示曲線的變化趨勢,若擬合曲線的整體趨勢與預期趨勢一致,則表明擬合效果較好。
(2)檢查數值結果“Evaluate the Numerical Fit Results”。擬合結果指標“Goodness of fit”中包含4個指標,分別為:SSE(the sum of squares due to error)誤差平方和、R-square復相關系數、AdjustedR-square調整自由度的復相關系數、RMSE(Root mean squared error)均方根誤差。
當SSE值和RMSE值越小,R-square值越接近于1時,說明擬合效果越好。
(3)若結果中顯示“Fit computation did not converge”或“Equation is badly conditioned”,則說明擬合效果較差??赏ㄟ^改變最高階次數、調整擬合算法或選擇“居中和縮放”來加以改善擬合效果[10]。
(4) 對于多項式擬合,多項式階數越高,若擬合系數的置信區間接近0,說明可能為過度擬合。
本研究以外特性曲線中燃油消耗率be與發動機轉速ne的擬合為例,根據之前所確定的外特性函數形式,選擇多項式擬合方法,分別作1~9階的多項式擬合,每階分別選擇OFF(常規最小二乘法)、LAR(最小絕對殘差法)、Bisquare(加權最小二乘法,默認算法)3種回歸模式,進行擬合效果對比,如圖10—圖12所示。圖10—圖12中分別標出各擬合結果指標的極大值和極小值。當SSE和RMSE值越小,R-square值越接近于1時,說明擬合效果越好。
圖10顯示RMSE在9階多項式擬合的OFF模式(常規最小二乘法)下取得極小值,圖11顯示SSE在9階多項式擬合的OFF模式(常規最小二乘法)下取得極小值,圖12顯示在9階多項式擬合的OFF模式(常規最小二乘法)下R-square和AdjustedR-square值均趨近與1。
綜上所述,在9階多項式擬合的OFF模式(常規最小二乘法)時,擬合效果最好。由此,得到最終擬合公式為:
be=1.128×ne9+0.7566×ne8-7.806×ne7-
4.719×ne6+20.37×ne5+7.679×ne4-
27.07×ne3+15.97×ne2+17.85×ne+
200.1。
(10)

圖10 1~9階多項式擬合RMSE值Fig.10 1—9 order polynomial fitting RMSE value
.
將最佳擬合(9階多項式擬合的OFF模式)和最差擬合(此處選取1階多項式擬合的LAR)的擬合效果圖作對比,如圖13所示,兩者的殘差圖作對比,如圖14所示。

圖13 最佳擬合和最差擬合的擬合效果對比Fig.13 Comparison of fitting effects in the best & worst fit

圖14 最佳擬合和最差擬合的殘差圖對比Fig. 14 Comparison of residual plots for best & worst fit
由圖13和圖14反映出,若擬合曲線與擬合多項式的基底高度相關(整體趨勢一致),且殘差絕對值較小,則表明擬合效果較好。
利用Matlab曲線擬合工具箱cftool進行曲線擬合更為簡單便捷。這種方式特別適用于自編程序或聯合仿真。
初步擬定的函數關系式對后續擬合方法的選擇至關重要。另外,在擬合過程中,根據初步擬定的函數關系式,可以多加嘗試不同的擬合方法,對比其SSE(誤差平方和)、R-square(復相關系數)、Adjusted R-square(調整自由度的復相關系數)、RMSE(均方根誤差),據此選擇最佳的擬合方法,得到更精確的擬合公式,提高后續仿真模擬的精度,增強健壯性和可靠性。