李頂根，劉剛，汪榮青

（1.華中科技大學(xué) 能源與動(dòng)力工程學(xué)院，湖北，武漢 430074；2.浙江機(jī)電職業(yè)技術(shù)學(xué)院，浙江，杭州 310053）

汽油機(jī)常采用油膜動(dòng)態(tài)方程設(shè)計(jì)空燃比補(bǔ)償器，這種基于氣道油膜辨識(shí)的空燃比控制方法，在發(fā)動(dòng)機(jī)穩(wěn)態(tài)運(yùn)行時(shí)，較為適用，而在發(fā)動(dòng)機(jī)瞬態(tài)時(shí)，辨識(shí)誤差較大[1]。美國FTP75測試循環(huán)[2]表明，發(fā)動(dòng)機(jī)80%的排放物均在發(fā)動(dòng)機(jī)冷啟動(dòng)過程的瞬態(tài)運(yùn)行工況下產(chǎn)生?？梢娝矐B(tài)工況下的空燃比精確控制（考慮油膜動(dòng)態(tài)效應(yīng)）是改善發(fā)動(dòng)機(jī)缸內(nèi)空燃比波動(dòng)，從而降低污染物排放的有效手段。為此，本文擬采用子空間辨識(shí)N4SID算法，辨識(shí)多影響參數(shù)下缸內(nèi)空燃比動(dòng)態(tài)波動(dòng)模型，分析多參數(shù)對(duì)辨識(shí)出的空燃比動(dòng)態(tài)模型的影響，分析模型的抗擾性能，為空燃比的控制建模打下基礎(chǔ)。

1 汽油機(jī)油膜動(dòng)態(tài)方程辨識(shí)及空燃比控制建模

汽油機(jī)采用進(jìn)氣道噴射方式，燃油從噴油器噴出后，一部分以蒸汽形式存在于氣道中，另一部分直接附著在壁面上，形成附壁油膜。油膜的存在對(duì)發(fā)動(dòng)機(jī)實(shí)際燃燒空燃比有很大影響，特別是發(fā)動(dòng)機(jī)冷啟動(dòng)時(shí)。研究發(fā)現(xiàn)，美國聯(lián)邦FTP75測試循環(huán)冷啟動(dòng)過程中所排放出的HC及CO占整個(gè)循環(huán)排放物的70%～80%，而很大部分原因歸結(jié)于冷啟動(dòng)過程中空燃比波動(dòng)過大[2]。

因此，對(duì)汽油機(jī)進(jìn)氣道油膜動(dòng)態(tài)模型建模分析的研究引起了國內(nèi)外學(xué)者的廣泛興趣。最有效的油膜物理模型建模研究始于1981年Aquino所提出的X-x模型，他將由噴油器噴出的燃油分成以X分?jǐn)?shù)沉積在壁面上的油膜以及以（1-X）分?jǐn)?shù)懸浮在進(jìn)氣道中的燃油蒸汽及微小液珠。附著在壁面上的油膜的蒸發(fā)時(shí)間常數(shù)是x[3]。之后于1992年，Hendricks E.和Vesterholm 提出了雙時(shí)間常數(shù)模型，該模型額外考慮了燃油蒸汽和液珠由進(jìn)氣道進(jìn)入氣缸的輸運(yùn)時(shí)間，其在某些工況下的模擬結(jié)果較優(yōu)于X-x模型[4]。之后很多國外學(xué)者也對(duì)進(jìn)氣道油膜進(jìn)行過深入的研究，在上述兩種模型的基礎(chǔ)上，提出了不同的改進(jìn)模型[2，5-6]。國內(nèi)有人采用最小二乘法、擴(kuò)展卡爾曼濾波等方法對(duì)油膜X- x模型進(jìn)行辨識(shí)，并在此基礎(chǔ)上設(shè)計(jì)了油膜補(bǔ)償器，對(duì)各種工況下的噴油量進(jìn)行補(bǔ)償達(dá)到精確控制空燃比的目的[7-9]。

上述傳統(tǒng)的辨識(shí)方法，均基于油膜物理模型的推導(dǎo)及建立，過程繁雜，不利于考慮發(fā)動(dòng)機(jī)多運(yùn)行參數(shù)對(duì)缸內(nèi)空燃比的影響。子空間辨識(shí)方法直接從輸入輸出數(shù)據(jù)中提取狀態(tài)空間參數(shù)，可以綜合考慮多參數(shù)對(duì)空燃比模型的影響，同時(shí)在辨識(shí)時(shí)間上要較傳統(tǒng)辨識(shí)方法少，這主要是因?yàn)樽涌臻g方法沒有循環(huán)迭代的過程[10]。

2 研究對(duì)象及模型簡述

試驗(yàn)系統(tǒng)采用某型汽油發(fā)動(dòng)機(jī)為研究對(duì)象，將實(shí)物系統(tǒng)簡化，在GT-Power軟件中搭建相應(yīng)的發(fā)動(dòng)機(jī)模型。發(fā)動(dòng)機(jī)基本參數(shù)見表1。

表1 發(fā)動(dòng)機(jī)的主要參數(shù)

模型主體是GT-Power軟件中自帶的具有高度保真性的算例模型，本文是在該模型的基礎(chǔ)上修改相關(guān)參數(shù)，并進(jìn)行研究。

假設(shè)汽油機(jī)缸內(nèi)空燃比動(dòng)態(tài)模型的狀態(tài)空間方程如下：

且做如下假設(shè)（即系統(tǒng)噪聲與測量噪聲不相關(guān)）：

定義輸入數(shù)據(jù)Hankel矩陣如下：

輸出yk的Hankel矩陣Y02i-1定義同上。

3 子空間N4SID算法及所辨識(shí)的模型

子空間辨識(shí)算法能夠根據(jù)給定輸入輸出數(shù)據(jù)的Hankel矩陣，通過空間正交投影或斜向投影的數(shù)學(xué)方法，直接估計(jì)出系統(tǒng)的狀態(tài)序列。從而將狀態(tài)空間的求解問題轉(zhuǎn)化為最小二乘問題，再求解出狀態(tài)空間矩陣A、B、C、D。它大大簡化了系統(tǒng)的建模過程，去除了繁瑣的物理模型推導(dǎo)過程。

假設(shè)式（1）所描述的系統(tǒng)是可觀的，則可設(shè)計(jì)Kalman濾波器對(duì)其狀態(tài)進(jìn)行估計(jì)。

K為Kalman增益，令

系統(tǒng)方程變?yōu)?/p>

將yk的表達(dá)式代入到狀態(tài)方程中，可得

可推出

在應(yīng)用系統(tǒng)輸入輸出Hankel矩陣進(jìn)行子空間辨識(shí)時(shí)，需作如下假設(shè)進(jìn)行限制[11]。

（1）矩陣（A-KC）的特征值嚴(yán)格分布在單位圓內(nèi)。

（2）系統(tǒng)矩陣（A，C）是可觀的，（A，[B，K]）是可控的。

（3）誤差序列ek是穩(wěn)定的零均值、白噪聲序列，且其二階矩滿足

式中：dij為Kronecker符號(hào)。

（4）對(duì)于開環(huán)數(shù)據(jù)，輸入序列uk與誤差序列ek不相關(guān)。

（5）輸入信號(hào)是準(zhǔn)穩(wěn)態(tài)，且是持續(xù)激勵(lì)的。

針對(duì)上述假設(shè)，進(jìn)行發(fā)動(dòng)機(jī)開環(huán)試驗(yàn)設(shè)計(jì)，選擇某型汽油發(fā)動(dòng)機(jī)在GT-Power軟件中建模，試驗(yàn)轉(zhuǎn)速變化范圍為1 000～6 000 r/min，每100 r/min進(jìn)行數(shù)據(jù)采樣；噴油量變化范圍為3.5～5.5 g/s；進(jìn)氣溫度范圍為280～340 K，此三者為模型辨識(shí)的輸入量。缸內(nèi)空燃比為輸出量。辨識(shí)數(shù)據(jù)分布如圖1所示。

應(yīng)用上述試驗(yàn)數(shù)據(jù)進(jìn)行系統(tǒng)辨識(shí)，數(shù)據(jù)樣本為1 000組。首先，對(duì)式（4）的輸出方程進(jìn)行擴(kuò)展和迭代，得到新的輸出方程為

對(duì)式（5）進(jìn)行迭代，可得Xk的表達(dá)式為

由（9）式可知，僅僅應(yīng)用已知的輸入輸出Hankel矩陣便可估計(jì)系統(tǒng)可觀測矩陣C，其估計(jì)值為

計(jì)算出C矩陣的估計(jì)值后，通過下式計(jì)算狀態(tài)矩陣的估計(jì)值。

從而式（1）所表示的系統(tǒng)模型辨識(shí)問題，轉(zhuǎn)化為一個(gè)最小二乘求解問題[10]。

利用計(jì)算出的系統(tǒng)矩陣的估計(jì)值，可計(jì)算系統(tǒng)殘差的協(xié)方差矩陣，協(xié)方差矩陣的求解仍然轉(zhuǎn)化為最小二乘的求解問題[10]。

根據(jù)上述算法，編寫Matlab計(jì)算程序，程序計(jì)算出的系統(tǒng)階數(shù)圖如圖2所示。

圖2中縱坐標(biāo)為系統(tǒng)奇異值的log函數(shù)值分布，橫坐標(biāo)為程序所設(shè)定的階數(shù)范圍。由圖2可知，當(dāng)系統(tǒng)階數(shù)n大于3后，其所對(duì)應(yīng)的奇異值的log函數(shù)值越來越小，故程序計(jì)算出3為系統(tǒng)的最佳階數(shù)。

確定模型階數(shù)后，即可利用發(fā)動(dòng)機(jī)開環(huán)試驗(yàn)數(shù)據(jù)對(duì)系統(tǒng)進(jìn)行辨識(shí)，辨識(shí)出的系統(tǒng)矩陣為

由于試驗(yàn)設(shè)計(jì)為開環(huán)狀態(tài)下的試驗(yàn)，故直接饋入矩陣D為0，所以程序中未予計(jì)算。

4 結(jié)果分析

針對(duì)計(jì)算出的系統(tǒng)模型，使用未參與辨識(shí)的200組數(shù)據(jù)進(jìn)行驗(yàn)證，驗(yàn)證數(shù)據(jù)的分布如圖3所示。

將驗(yàn)證數(shù)據(jù)代入所辨識(shí)出的系統(tǒng)模型，模型計(jì)算出的數(shù)據(jù)與試驗(yàn)數(shù)據(jù)的比較如圖4所示。

圖4上半部分為預(yù)測出的空燃比與試驗(yàn)值的絕對(duì)誤差分布圖；下半部分為相對(duì)誤差分布圖。由圖4可知，模型預(yù)測的空燃比絕對(duì)誤差分布在-0.2647～0.4475之間，相對(duì)誤差分布在-2.25%～3.59%之間，模型計(jì)算的結(jié)果與試驗(yàn)值吻合程度高。說明子空間辨識(shí)算法N4SID在辨識(shí)空燃比模型方面準(zhǔn)確、有效。由于轉(zhuǎn)速信號(hào)的周期性采樣，計(jì)算結(jié)果也大致呈周期性分布。

驗(yàn)證數(shù)據(jù)的輸入值都假定是準(zhǔn)確和穩(wěn)定的。但在發(fā)動(dòng)機(jī)實(shí)際運(yùn)行過程中，對(duì)數(shù)據(jù)進(jìn)行采集時(shí)均會(huì)夾雜外界的干擾信號(hào)，絕對(duì)準(zhǔn)確無擾動(dòng)的數(shù)據(jù)在試驗(yàn)過程中不可能產(chǎn)生。為此，為驗(yàn)證所辨識(shí)出的模型抗擾動(dòng)的能力，在驗(yàn)證數(shù)據(jù)中分別單獨(dú)對(duì)輸入信號(hào)添加零均值，方差為信號(hào)平均值±1%范圍內(nèi)分布的白噪聲信號(hào)，以此檢驗(yàn)?zāi)Ｐ偷目箶_動(dòng)能力。

添加擾動(dòng)信號(hào)后，空燃比驗(yàn)證誤差分布圖如圖5和圖6所示。

由圖5、圖6可知，轉(zhuǎn)速加入白噪聲信號(hào)后，輸出數(shù)據(jù)的絕對(duì)誤差與相對(duì)誤差幾乎與不添加白噪聲信號(hào)時(shí)的相同，可見辨識(shí)出的模型對(duì)轉(zhuǎn)速的變化有較大的抗擾能力。噴油量信號(hào)加入白噪聲后，空燃比的絕對(duì)誤差與相對(duì)誤差波動(dòng)幅度小，模型對(duì)噴油量的干擾信號(hào)具備抗擾能力。當(dāng)對(duì)進(jìn)氣量添加白噪聲后，模型預(yù)測出的空燃比絕對(duì)誤差、相對(duì)誤差均產(chǎn)生大幅度的波動(dòng)，相對(duì)誤差限躍升至-14%～10%，說明辨識(shí)出的模型對(duì)進(jìn)氣溫度的波動(dòng)比較敏感。分析其原因，可能是用于辨識(shí)的數(shù)據(jù)中進(jìn)氣溫度變化頻率最小，每250個(gè)采樣點(diǎn)才產(chǎn)生進(jìn)氣溫度的階躍變化，由此導(dǎo)致辨識(shí)出的模型對(duì)進(jìn)氣溫度變化敏感。

將辨識(shí)所得的模型應(yīng)用于不同發(fā)動(dòng)機(jī)工況，此處任意選取4個(gè)工況點(diǎn)進(jìn)行對(duì)比，將進(jìn)氣溫度和噴油量分別定為283 K、4.1 g/s，轉(zhuǎn)速分別為1 372 r/min、4 859 r/min、5 128 r/min、5 732 r/min。其所預(yù)測出的發(fā)動(dòng)機(jī)空燃比與發(fā)動(dòng)機(jī)仿真試驗(yàn)所得的空燃比見表2。

表2 模型預(yù)測空燃比與發(fā)動(dòng)機(jī)仿真試驗(yàn)空燃比對(duì)比

由表2可知，在所選取的工況點(diǎn)下，經(jīng)辨識(shí)得出的模型預(yù)測出的空燃比值與發(fā)動(dòng)機(jī)仿真試驗(yàn)所得的空燃比值符合度好，預(yù)測的相對(duì)誤差在±5%以內(nèi)。

經(jīng)以上分析可知，對(duì)于多輸入下的發(fā)動(dòng)機(jī)空燃比動(dòng)態(tài)模型建模，傳統(tǒng)的基于物理分析建模的方法較為復(fù)雜，且不易于考慮多參數(shù)的影響。而應(yīng)用子空間方法針對(duì)試驗(yàn)采集的大量樣本數(shù)據(jù)，進(jìn)行快速高效辨識(shí)，可以解決物理建模繁瑣低效的缺陷，且辨識(shí)的精度高，抗輸入干擾的能力強(qiáng)，是發(fā)動(dòng)機(jī)瞬態(tài)空燃比控制建模的有效方法。

5 結(jié)論

（1）推導(dǎo)了子空間辨識(shí)N4SID算法的主要過程，基于該算法的限制條件在GT-Power軟件中設(shè)計(jì)了發(fā)動(dòng)機(jī)開環(huán)狀態(tài)試驗(yàn)，采集了1 200組各工況下的數(shù)據(jù)，其中1 000組數(shù)據(jù)用于模型的辨識(shí)，200組數(shù)據(jù)用于模型的驗(yàn)證。

（2）辨識(shí)程序辨識(shí)出的3階狀態(tài)空間模型具有最優(yōu)的精度，模型計(jì)算出的空燃比絕對(duì)誤差限為-0.2647～0.4475、相對(duì)誤差限為 -2.25%～ 3.59%。

（3）向驗(yàn)證數(shù)據(jù)輸入信號(hào)添加零均值，方差為各信號(hào)平均值1%的白噪聲后發(fā)現(xiàn)，辨識(shí)出的模型對(duì)轉(zhuǎn)速、噴油量的波動(dòng)抵抗力強(qiáng)。而對(duì)進(jìn)氣溫度變化敏感，究其原因可能在于采樣數(shù)據(jù)中進(jìn)氣溫度變化頻率低，造成模型對(duì)進(jìn)氣溫度的敏感性增強(qiáng)。

（4）模型預(yù)測值與GT-Power仿真試驗(yàn)值的相對(duì)誤差在±5%以內(nèi)，辨識(shí)出模型的精度符合空燃比控制建模需求。

[1]BRUNO B A，SANTAVICCA D A，ZELLO J V. LIF Characterization of Intake Valve Fuel Films during Cold Start in a PFI Engine[C].SAE Paper 2002-01-2751.

[2]NEYACHENKO I I. Method of A/F Control during SI Engine Cold Start[C].SAE Paper 981171.

[3]AQUINO C F. Transient A/F Control Characteristics of the 5 Liter Central Fuel Injected Engine[C]. SAE Paper 810494.

[4]HENDRICKS E，VESTERHOLM T，KAIDANTZIS P，et al. Nonlinear Transient Fuel Film Compensation (NTFC)[C].SAE Paper 930767.

[5]Xu Hongming. Control of A/F Ratio During Engine Transients[C]. SAE Paper 1999-01-1484.

[6]ARSIE I， PIANESE C， RIZZOA G，et al. An Adaptive Estimator of Fuel Film Dynamics in the Intake Port of a Spark Ignition Engine[J]. Control Engineering Practice， 2003，11(3)：303-309.

[7]鐘祥麟. 基于油膜模型的多點(diǎn)噴射汽油機(jī)瞬態(tài)工況控制研究[D].長春：吉林大學(xué)， 2007.Zhong Xianglin.Study on Control of MPI Gasoline Engine in Transient Condition Based on Fuel Film Model[D].Changchun：Jilin University，2007.（in Chinese）

[8]李頂根，舒詠強(qiáng). 汽油機(jī)進(jìn)氣道油膜模型參數(shù)辨識(shí)算法的研究[J].內(nèi)燃機(jī)學(xué)報(bào)，2009，5(4)：15-22.Li Dinggen，Shu Yongqiang. Research of Parameter Identification Algorithm of Intake Port Fuel Film Model for Gasoline Engine[J].Transactions of CSICE，2009，5(4)：15-22.（in Chinese）

[9]姚棟偉，吳鋒，朱一峰，等. 基于 EKF 算法的進(jìn)氣道燃油狀態(tài)觀測器研究[J]. 內(nèi)燃機(jī)學(xué)報(bào)，2010，4(28)：374-379.Yao Dongwei，Wu Feng，Zhu Yifeng，et al.Study on Observer of Fuel Dynamics in Intake Port of SI Engines Based on EKF Algorithm[J].Transactions of CSICE，2010，4(28)：374-379.（in Chinese）

[10]Van OVERSCHEE P， de MOOR B. Subspace Identification for Linear Systems：Theory-Implementation-Application[M].Dordrecht：Kluwer Academic Publishers，1996.

[11]QIN S J. An Overview of Subspace Identification[J].Computers and Chemical Engineering .2006(30)：1502-1513.