管湘源，儲(chǔ)江偉，高偉健，趙小婷

(東北林業(yè)大學(xué) 交通學(xué)院，哈爾濱 150040)

0 引言

發(fā)動(dòng)機(jī)是車輛的動(dòng)力系統(tǒng)，是一個(gè)復(fù)雜的非線性、時(shí)變性及滯后系統(tǒng)[1]，使用真實(shí)發(fā)動(dòng)機(jī)來設(shè)計(jì)控制器是一項(xiàng)費(fèi)時(shí)費(fèi)力的任務(wù)。因此，許多研究人員通過建立發(fā)動(dòng)機(jī)數(shù)學(xué)模型來設(shè)計(jì)、測(cè)試和優(yōu)化其控制器。汽車行業(yè)已經(jīng)轉(zhuǎn)向基于模型的校準(zhǔn)來尋求解決方案[2]。當(dāng)前主流發(fā)動(dòng)機(jī)建模方法有兩種：一種是解析法，即依據(jù)發(fā)動(dòng)機(jī)的機(jī)理和相應(yīng)參數(shù)關(guān)系而建立仿真數(shù)學(xué)模型；第二種是實(shí)驗(yàn)法，即通過發(fā)動(dòng)機(jī)臺(tái)架試驗(yàn)數(shù)據(jù)建立仿真模型[3-5]。

通過解析法建立數(shù)學(xué)模型，通常需要根據(jù)數(shù)據(jù)擬合(一般為曲面擬合法、多元插值法和綜合法)來確定相關(guān)系數(shù)(即曲線或曲面的待定系數(shù))，進(jìn)而得到最終的數(shù)學(xué)擬合公式。利用Matlab擬合工具箱進(jìn)行數(shù)據(jù)擬合，具有求解便捷和精度較高等特點(diǎn)，因而在機(jī)械、物理、電力、水利、醫(yī)藥和能源等多種學(xué)科領(lǐng)域均有廣泛應(yīng)用。在機(jī)械領(lǐng)域應(yīng)用方面，鄭嘯洲等[6]通過Matlab擬合工具箱對(duì)某型發(fā)動(dòng)機(jī)試驗(yàn)數(shù)據(jù)進(jìn)行擬合處理，建立發(fā)動(dòng)機(jī)的數(shù)學(xué)模型。在物理領(lǐng)域應(yīng)用方面，陳文芳等[7]和魏奶萍[8]通過Matlab曲線擬合工具箱處理PN結(jié)的伏安特性數(shù)據(jù),擬合得到相應(yīng)的PN結(jié)伏安特性函數(shù)表達(dá)式；郭悅韶[9]通過Matlab曲線擬合工具箱處理霍爾效應(yīng)實(shí)驗(yàn)測(cè)量數(shù)據(jù)，擬合得到磁場分布曲線；范雋宏等[10]通過Matlab曲線擬合工具箱對(duì)軟磁材料的單位損耗數(shù)據(jù)的曲線擬合分析，得出用于變壓器電磁計(jì)算的擬合方程；岳鵬等[11]以測(cè)量玻璃熱膨脹系數(shù)和折射率溫度系數(shù)實(shí)驗(yàn)、太陽能電池基本特性測(cè)量實(shí)驗(yàn)為例，通過Matlab曲線擬合工具箱處理物理實(shí)驗(yàn)數(shù)據(jù)，進(jìn)行擬合實(shí)驗(yàn)曲線。在電力領(lǐng)域應(yīng)用方面，李太興等[12-13]通過Matlab曲線擬合工具箱進(jìn)行火電廠運(yùn)行的機(jī)組負(fù)荷優(yōu)化分配；高艷平等[14]通過Matlab曲線擬合工具箱分析某發(fā)電廠的地基沉降實(shí)測(cè)資料，并擬合出多種不同的地基沉降預(yù)測(cè)模型。在水利領(lǐng)域應(yīng)用方面，王杰等[15]通過Matlab曲線擬合工具箱求解分析抽水試驗(yàn)資料，并確定承壓含水層的導(dǎo)水系數(shù)和導(dǎo)壓系數(shù)。在醫(yī)藥領(lǐng)域應(yīng)用方面，曹俊涵等[16]通過Matlab曲線擬合工具箱處理藥物溶出數(shù)據(jù)，并計(jì)算藥物溶出度Weibull分布參數(shù)。在能源領(lǐng)域應(yīng)用方面，高鵬[17]通過Matlab曲線擬合工具箱對(duì)浮選尾礦灰分?jǐn)?shù)據(jù)與圖像灰度數(shù)據(jù)進(jìn)行曲線擬合；何定武[18]通過Matlab曲線擬合工具箱分析風(fēng)機(jī)特性曲線并計(jì)算風(fēng)機(jī)實(shí)際工況下節(jié)電率。

在上述應(yīng)用中，經(jīng)過驗(yàn)證，由Matlab曲線擬合工具箱得到的擬合結(jié)果具有較高精度，因此本研究也選用Matlab曲線擬合工具箱作為曲線擬合工具，以某款發(fā)動(dòng)機(jī)的外特性和萬有特性曲線為例，通過Matlab曲線擬合工具箱cftool對(duì)提取的數(shù)據(jù)按照不同的擬合方法進(jìn)行公式擬合，并根據(jù)擬合優(yōu)度檢驗(yàn)結(jié)果選擇出最佳擬合公式，此公式即為最終確定的外特性曲線的數(shù)學(xué)模型。

1 發(fā)動(dòng)機(jī)數(shù)值建模原理

以發(fā)動(dòng)機(jī)臺(tái)架試驗(yàn)數(shù)據(jù)的離散點(diǎn)或萬有特性圖數(shù)字化提取的坐標(biāo)點(diǎn)數(shù)據(jù)作為樣點(diǎn)，選擇一種合適的函數(shù)逼近方式，將它重構(gòu)成確定的數(shù)學(xué)模型，這是建模的關(guān)鍵。

發(fā)動(dòng)機(jī)數(shù)學(xué)建模的原理是基于函數(shù)逼近，即函數(shù)的近似表示問題。函數(shù)逼近的方法多種多樣，具體分類如圖1所示。

圖1 函數(shù)逼近方法分類Fig.1 Function approximation method classification

由圖1可知，函數(shù)逼近的常規(guī)方法為擬合法、插值法和綜合法。

1.1 曲面(線)擬合法

曲面(線)擬合法是選定近似函數(shù)H(x)后，不要求H(x)通過已知樣本點(diǎn)，只要求在某種意義下H(x)在這些樣點(diǎn)的總偏值最小。采用多元線性回歸方法進(jìn)行曲面擬合，構(gòu)造發(fā)動(dòng)機(jī)穩(wěn)態(tài)輸出轉(zhuǎn)矩、油門開度以及發(fā)動(dòng)機(jī)轉(zhuǎn)速的關(guān)系。發(fā)動(dòng)機(jī)的負(fù)荷特性數(shù)學(xué)模型構(gòu)建為[19]：

(1)

式中：Me為發(fā)動(dòng)機(jī)輸出轉(zhuǎn)矩，Nm；α為發(fā)動(dòng)機(jī)的油門開度，%；ne為發(fā)動(dòng)機(jī)的轉(zhuǎn)速，r/min；B為系數(shù)矩陣。

通過最小二乘原理來確定系數(shù)矩陣B。建立的回歸方程為[20]：

(2)

式中：b,e為系數(shù)矩陣B的余數(shù)。

寫成矩陣形式如下：

Z=G·B+E。

(3)

G的列數(shù)k與多項(xiàng)式階數(shù)s的關(guān)系為：k=(s+1)(s+2)/2。其擬合值和測(cè)量值的逼近程度可以用擬合度來評(píng)價(jià)。

1.2 多元插值法

插值法是給出函數(shù)f(x)的一些樣點(diǎn)值，選定某些便于計(jì)算的函數(shù)，要求函數(shù)通過已知樣點(diǎn)，由此確定函數(shù)H(x)為f(x)近似值。

采用多元插值法，構(gòu)造發(fā)動(dòng)機(jī)穩(wěn)態(tài)輸出轉(zhuǎn)矩、油門開度以及發(fā)動(dòng)機(jī)轉(zhuǎn)速的關(guān)系。發(fā)動(dòng)機(jī)的負(fù)荷特性數(shù)學(xué)模型構(gòu)建為:

Me(t)=f(ne(t),αe(t))。

(4)

式中：Me(t)為發(fā)動(dòng)機(jī)輸出轉(zhuǎn)矩，Nm；αe為發(fā)動(dòng)機(jī)的油門開度，%；ne為發(fā)動(dòng)機(jī)的轉(zhuǎn)速，rad/s。

要計(jì)算的點(diǎn)ne、αe落在4個(gè)網(wǎng)格點(diǎn)中間，即

(5)

定義：

(6)

則：

(7)

其中，cij是插值系數(shù)，插值問題的核心就是求出插值系數(shù)，插值算法不同所求出的插值系數(shù)各異。

Matlab中的曲線擬合等相應(yīng)的工具箱或cubic、linear、neast、spline、interpn、interpft、ndgrid等函數(shù)可以用來實(shí)現(xiàn)插值的功能。下面對(duì)比分析幾種常用的插值方法的特點(diǎn)及用途，如圖2所示。

圖2 幾種常用的插值方法Fig. 2 Several interpolation methods

1.3 綜合法

綜合法是融合多項(xiàng)式擬合與一元函數(shù)插值的方法。它首先根據(jù)最小二乘原理，擬合得到每個(gè)轉(zhuǎn)速下的T=f(α)；然后對(duì)于某一工況點(diǎn)A，根據(jù)轉(zhuǎn)速判斷最臨近于工況點(diǎn)A的3條負(fù)載特性曲線，再用該曲線的多項(xiàng)式計(jì)算出來A點(diǎn)α對(duì)應(yīng)的轉(zhuǎn)矩Ti-1、Ti、Ti+1的值，最后用插值法求得T(α)[23]。

2 萬有特性圖數(shù)字化方法

傳統(tǒng)數(shù)據(jù)提取方法是進(jìn)行手工測(cè)量和估算，此種方法不但誤差大，而且效率低。因此，有必要采取數(shù)字化方法提取萬有特性曲線圖的數(shù)據(jù)。選用圖像數(shù)字化軟件Get Data作為數(shù)據(jù)提取的工具，其提取數(shù)據(jù)流程如圖3所示。

提取數(shù)據(jù)前，需對(duì)圖像進(jìn)行預(yù)處理，包括圖像的放大、截取和擺正等清晰化過程。良好的預(yù)處理可以在一定范圍內(nèi)降低視覺誤差和分辨率誤差，從而使獲得的數(shù)據(jù)更為精確。

文中選取某款汽油機(jī)，根據(jù)廠方提供的萬有特性曲線(圖4)和外特性曲線(圖5)的圖像數(shù)據(jù)，對(duì)圖上的數(shù)據(jù)進(jìn)行采集，提取為數(shù)據(jù)表形式，為后期建立數(shù)學(xué)模型奠定基礎(chǔ)。

圖3 Get Data提取數(shù)據(jù)流程圖Fig.3 Get Data extraction data flow chart

將萬有特性曲線(圖4)和外特性曲線(圖5)分別導(dǎo)入Get Data軟件，經(jīng)過設(shè)置曲線顏色、設(shè)置背景色、設(shè)置基準(zhǔn)點(diǎn)與取值范圍的幾個(gè)前期過程后，選擇曲線相應(yīng)坐標(biāo)點(diǎn)，進(jìn)行對(duì)應(yīng)數(shù)據(jù)的提取。萬有特性曲線的數(shù)據(jù)提取，如圖6所示。外特性曲線的數(shù)據(jù)提取，如圖7所示。

圖4 發(fā)動(dòng)機(jī)萬有特性曲線Fig.4 Engine universal characteristic curve

圖5 發(fā)動(dòng)機(jī)外特性曲線Fig.5 Engine external characteristic curve

圖6 萬有特性圖坐標(biāo)點(diǎn)數(shù)據(jù)提取Fig. 6 Universal feature map coordinate point data extraction

圖7 外特性圖坐標(biāo)點(diǎn)數(shù)據(jù)提取Fig.7 External characteristic map coordinate point data extraction

提取的所有坐標(biāo)點(diǎn)的數(shù)據(jù)將顯示在工作區(qū)間內(nèi)，可將數(shù)據(jù)導(dǎo)出為表格形式，見表1。

表1 導(dǎo)出的部分?jǐn)?shù)據(jù)

除此之外，也可以選用Origin Pro軟件進(jìn)行數(shù)據(jù)提取，其提取具體流程如圖8所示。

3 發(fā)動(dòng)機(jī)數(shù)學(xué)建模

3.1 外特性模型

外特性建模采用曲線擬合的方法。首先初步確定外特性曲線的函數(shù)表達(dá)式(即確定近似函數(shù))，進(jìn)而選擇相應(yīng)的擬合方式得到最優(yōu)擬合系數(shù)，最終得到完整的外特性數(shù)學(xué)模型。

圖8 Origin Pro提取數(shù)據(jù)流程圖Fig.8 Origin Pro extraction data flow

3.1.1 函數(shù)形式確定

發(fā)動(dòng)機(jī)外特性是研究汽車動(dòng)力性的基礎(chǔ),它是在發(fā)動(dòng)機(jī)節(jié)氣門全開的情況下,發(fā)動(dòng)機(jī)轉(zhuǎn)矩Me(或功率Pe)、燃油消耗率be與發(fā)動(dòng)機(jī)曲軸轉(zhuǎn)速ne之間的關(guān)系,可以看作是發(fā)動(dòng)機(jī)轉(zhuǎn)速的一元函數(shù)[24]。

分別用x表示ne，y表示Me、Pe或be，n為擬合多項(xiàng)式的最高階數(shù)。

假定其函數(shù)形式為：

y=p1+p2x+p3x2+…+pnxn。

(8)

3.1.2 擬合方式選擇

本文以Matlab的曲線擬合工具箱cftool作為擬合工具。cftool中包含的擬合方法有：Exponential指數(shù)擬合、Fourier傅里葉擬合、Gaussian高斯法、Interpolant內(nèi)插法、Polynomial多項(xiàng)式(1～9階)、Rational有理擬合、Power冪指數(shù)擬合、Smoothing spline平滑樣條擬合、Sum of Sin Functions正弦曲線擬合、Weibull兩個(gè)參數(shù)的Weibull分布，并支持Custom Equations用戶自定義。

本文選用多項(xiàng)式擬合方法，并分析擬合效果，在后續(xù)擬合優(yōu)度檢驗(yàn)一節(jié)中給出詳細(xì)說明。

發(fā)動(dòng)機(jī)外特性模型按以下流程進(jìn)行擬合：將發(fā)動(dòng)機(jī)轉(zhuǎn)矩Me(或功率Pe)、燃油消耗率be和發(fā)動(dòng)機(jī)轉(zhuǎn)速ne的數(shù)據(jù)分別導(dǎo)入到cftool工具箱中，發(fā)動(dòng)機(jī)轉(zhuǎn)速ne設(shè)置為x，燃油消耗率be、發(fā)動(dòng)機(jī)功率Pe和發(fā)動(dòng)機(jī)轉(zhuǎn)矩Me分別依次設(shè)置為y。隨后設(shè)置權(quán)重，默認(rèn)值為1。完成初始數(shù)據(jù)設(shè)定后，選擇擬合方法。擬合運(yùn)算執(zhí)行后，其結(jié)果在Results中顯示，如圖9所示。

圖9 擬合結(jié)果Fig.9 The results of fitting

在擬合過程中，選擇“Center and scale”(標(biāo)準(zhǔn)化)，可減少出現(xiàn)病態(tài)矩陣的情況。采用擬合算法時(shí)應(yīng)盡量選擇合適的最高冪次[25]，但并非次數(shù)越高越好。因?yàn)楫?dāng)次數(shù)過高，容易出現(xiàn)龍格(Runge)現(xiàn)象，導(dǎo)致在試驗(yàn)點(diǎn)以外的數(shù)據(jù)上偏差可能很大。

3.1.3 數(shù)學(xué)模型確定

通過多項(xiàng)式擬合得到外特性的最優(yōu)擬合系數(shù)，最終得到外特性數(shù)學(xué)模型，見表2。

表2 外特性曲線的最佳擬合公式Tab.2 Best fit formula of external characteristic curve

3.2 萬有特性模型(油耗模型)

發(fā)動(dòng)機(jī)的燃油消耗率be是發(fā)動(dòng)機(jī)轉(zhuǎn)速和有效轉(zhuǎn)矩(ne,Me)的函數(shù)，設(shè)發(fā)動(dòng)機(jī)萬有特性的函數(shù)類型為多項(xiàng)式[26-29]，其回歸模型的矩陣形式為[30-31]：

be=[1,ne,Me,ne2,neMe,Me2,…,neMel-1,Mel]×

(9)

式中：a0、a2、…、ak-1為模型中的待定系數(shù)；e0、e2、…、ek-1為隨機(jī)誤差；l為多項(xiàng)式的最高冪次；k為多項(xiàng)式的項(xiàng)數(shù)，k=(l+1)(l+2)/2[29-30]。

發(fā)動(dòng)機(jī)萬有特性曲線的數(shù)學(xué)模型在cftool工具箱中的建立流程類似于外特性曲線的建立流程，在此不做贅述。

4 擬合優(yōu)度檢驗(yàn)

為評(píng)價(jià)擬合的性能，需對(duì)擬合結(jié)果進(jìn)行擬合優(yōu)度檢驗(yàn)。具體的檢驗(yàn)方法為：

(1)查看圖的擬合結(jié)果“Examine the Graphical Fit Results”。選擇Axes Limits，放大x軸坐標(biāo)，顯示曲線的變化趨勢(shì)，若擬合曲線的整體趨勢(shì)與預(yù)期趨勢(shì)一致，則表明擬合效果較好。

(2)檢查數(shù)值結(jié)果“Evaluate the Numerical Fit Results”。擬合結(jié)果指標(biāo)“Goodness of fit”中包含4個(gè)指標(biāo)，分別為：SSE(the sum of squares due to error)誤差平方和、R-square復(fù)相關(guān)系數(shù)、AdjustedR-square調(diào)整自由度的復(fù)相關(guān)系數(shù)、RMSE(Root mean squared error)均方根誤差。

當(dāng)SSE值和RMSE值越小，R-square值越接近于1時(shí)，說明擬合效果越好。

(3)若結(jié)果中顯示“Fit computation did not converge”或“Equation is badly conditioned”，則說明擬合效果較差。可通過改變最高階次數(shù)、調(diào)整擬合算法或選擇“居中和縮放”來加以改善擬合效果[10]。

(4) 對(duì)于多項(xiàng)式擬合，多項(xiàng)式階數(shù)越高，若擬合系數(shù)的置信區(qū)間接近0，說明可能為過度擬合。

本研究以外特性曲線中燃油消耗率be與發(fā)動(dòng)機(jī)轉(zhuǎn)速ne的擬合為例，根據(jù)之前所確定的外特性函數(shù)形式，選擇多項(xiàng)式擬合方法，分別作1～9階的多項(xiàng)式擬合，每階分別選擇OFF(常規(guī)最小二乘法)、LAR(最小絕對(duì)殘差法)、Bisquare(加權(quán)最小二乘法，默認(rèn)算法)3種回歸模式，進(jìn)行擬合效果對(duì)比，如圖10—圖12所示。圖10—圖12中分別標(biāo)出各擬合結(jié)果指標(biāo)的極大值和極小值。當(dāng)SSE和RMSE值越小，R-square值越接近于1時(shí)，說明擬合效果越好。

圖10顯示RMSE在9階多項(xiàng)式擬合的OFF模式(常規(guī)最小二乘法)下取得極小值，圖11顯示SSE在9階多項(xiàng)式擬合的OFF模式(常規(guī)最小二乘法)下取得極小值，圖12顯示在9階多項(xiàng)式擬合的OFF模式(常規(guī)最小二乘法)下R-square和AdjustedR-square值均趨近與1。

綜上所述，在9階多項(xiàng)式擬合的OFF模式(常規(guī)最小二乘法)時(shí)，擬合效果最好。由此，得到最終擬合公式為：

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階多項(xiàng)式擬合RMSE值Fig.10 1—9 order polynomial fitting RMSE value

.

將最佳擬合(9階多項(xiàng)式擬合的OFF模式)和最差擬合(此處選取1階多項(xiàng)式擬合的LAR)的擬合效果圖作對(duì)比，如圖13所示，兩者的殘差圖作對(duì)比，如圖14所示。

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

圖14 最佳擬合和最差擬合的殘差圖對(duì)比Fig. 14 Comparison of residual plots for best & worst fit

由圖13和圖14反映出，若擬合曲線與擬合多項(xiàng)式的基底高度相關(guān)(整體趨勢(shì)一致)，且殘差絕對(duì)值較小，則表明擬合效果較好。

5 結(jié)論

利用Matlab曲線擬合工具箱cftool進(jìn)行曲線擬合更為簡單便捷。這種方式特別適用于自編程序或聯(lián)合仿真。

初步擬定的函數(shù)關(guān)系式對(duì)后續(xù)擬合方法的選擇至關(guān)重要。另外，在擬合過程中，根據(jù)初步擬定的函數(shù)關(guān)系式，可以多加嘗試不同的擬合方法，對(duì)比其SSE(誤差平方和)、R-square(復(fù)相關(guān)系數(shù))、Adjusted R-square(調(diào)整自由度的復(fù)相關(guān)系數(shù))、RMSE(均方根誤差)，據(jù)此選擇最佳的擬合方法，得到更精確的擬合公式，提高后續(xù)仿真模擬的精度，增強(qiáng)健壯性和可靠性。