An LPV Adaptive Observer for Updating a Map Applied to an MAF Sensor in a Diesel Engine

Abstract

In this paper, a new method for mass air flow (MAF) sensor error compensation and an online updating error map (or lookup table) due to installation and aging in a diesel engine is developed. Since the MAF sensor error is dependent on the engine operating point, the error model is represented as a two-dimensional (2D) map with two inputs, fuel mass injection quantity and engine speed. Meanwhile, the 2D map representing the MAF sensor error is described as a piecewise bilinear interpolation model, which can be written as a dot product between the regression vector and parameter vector using a membership function. With the combination of the 2D map regression model and the diesel engine air path system, an LPV adaptive observer with low computational load is designed to estimate states and parameters jointly. The convergence of the proposed algorithm is proven under the conditions of persistent excitation and given inequalities. The observer is validated against the simulation data from engine software enDYNA provided by Tesis. The results demonstrate that the operating point-dependent error of the MAF sensor can be approximated acceptably by the 2D map from the proposed method.

1. Introduction

Accurate exhaust gas recirculation (EGR) rate control and air-fuel ratio (AFR) control are important technologies to satisfy the increasingly stringent emission regulations, which are dependent on the precise calculation of the EGR rate and AFR [1,2]. The accuracy of the EGR rate and AFR can be improved by a mass air flow (MAF) sensor, in which a sensor element is heated to a fixed temperature, and the difference in temperature attributed to heat transfer from the heating element to the air flow is a measure of the air mass flow [3,4,5]. However, there are many different local flow fields within the inlet piping due to the three-dimensional turbulence flow, leading to measurement biases in the MAF sensor installed between the air filter and the intake manifold. In addition, the MAF sensor is also subjected to aging phenomena owing to the accumulation of dust on the sensing element, which causes the deterioration of the measurement accuracy [6,7]. These errors will bring about an inaccurate EGR rate and AFR and have adverse impacts on the emission performance of diesel engine.

It is difficult to accurately establish an analytical model for the MAF sensor error. In view of the relatively low computational load, maps (or lookup tables) have been widely used to characterize systems where the functional relationship is unavailable or too complex to represent analytically [8]. Therefore, the relative error of the MAF sensor is described as a one-dimensional (1D) map taking compressor mass air flow as input [2]. In order to track MAF sensor aging, the extended Kalman filter (EKF) for updating maps is presented in [9,10,11], in which the 1D map is represented as a piecewise linear interpolation model and the map parameters are considered as parameter states. Due to the piecewise linear interpolation model having the characteristic of partition calculation and due to the the map input being able to enter only one input interval of the 1D map at any time, then only two parameter states participating in linear interpolation are observable and the other not. Therefore, the error covariance matrix elements of EKF corresponding to the locally unobservable parameter states will increase linearly. Although the solution is to limit this growth in [9,10,11], the convergence of EKF with a confined covariance matrix cannot be guaranteed. In addition, the measurement error of the MAF sensor depends on the engine operating point, which is usually defined as fuel mass injection quantity and engine speed. The 1D map representing MAF sensor error ignores the engine speed, reducing the accuracy when the diesel engine is run over a wide speed range.

The adaptive observer with the advantage of simple convergence conditions is an alternative method for updating maps. Recursive algorithms designed for joint estimation of states and parameters in state space systems are usually known as adaptive observers, and some early works with adaptive observers to jointly estimate states and parameters in multi-input-multi-output linear time varying systems can be found in [12,13]. In order to estimate sensor faults, adaptive observers for linear time varying systems with unknown parameters in output equations have been studied [14,15]. However, the existing adaptive observers cannot directly update maps.

In this paper, an adaptive observer is developed to update the map, in which the MAF sensor error is described as a two-dimensional (2D) map taking the operating point as the input to improve the model accuracy comparing the 1D map. Then, two problems are studied. First, in order to expediently analyze and design the parameter estimation method, the input-output relationship of the MAF sensor error 2D map is expressed as a dot product between the regression vector and the unknown parameter vector. Second, based on the linear parameter varying (LPV) system of the diesel engine with EGR and variable geometry turbocharger (VGT), a 2D map estimation method with a simple structure and low computational load is designed to facilitate the algorithm implementation.

This paper is organized as follows. In Section 2, the 2D map is expressed as the dot product between the regression vector and the unknown parameter vector, and the estimation problem for a class of LPV systems with an unknown parameter vector is given. In Section 3, the LPV adaptive observer is proposed, as well as the convergence analysis. Simulation results from enDYNA are presented in Section 4, and the conclusions are summarized in Section 5.

2. Problem Formulation

2.1. A Diesel Engine Air Path LPV Model

Figure 1 shows the model structure of a diesel engine with EGR and VGT, and the model can be expressed as [16]:

$$\begin{array}{cc}{\dot{p}}_{im}\hfill & =\frac{R_a{T}_{im}}{V_{im}}\left(W_c+W_{egr}-W_{ei}\right)\hfill \\ {\dot{p}}_{em}\hfill & =\frac{R_e{T}_{em}}{V_{em}}\left(W_f+W_{ei}-W_t-W_{egr}\right)\hfill \\ {\dot{\omega }}_t\hfill & =\frac{P_t{\eta }_m-P_c}{J_t{\omega }_t}\hfill \end{array}$$

(1)

where $W_c$ is the compressor mass air flow, $W_{egr}$ is the EGR mass flow, $W_{ei}$ is the cylinder mass flow, $W_f$ is the fuel rate injected to cylinder, $W_t$ is the turbine mass flow, $P_t$ is the turbine power, $P_c$ is the compressor power, ${\eta }_m$ is the turbocharger mechanical efficiency, $p_{im}$ is the intake manifold pressure, $p_{em}$ is the exhaust manifold pressure and ${\omega }_t$ is the turbine speed.

Figure 1. Schematic of the diesel engine model with exhaust gas recirculation (EGR) and variable geometry turbocharger (VGT).

Figure 1. Schematic of the diesel engine model with exhaust gas recirculation (EGR) and variable geometry turbocharger (VGT).

Meanwhile, $W_c$, $W_{egr}$, $W_{ei}$, $W_t$, $W_f$, $P_c$ and $P_t{\eta }_m$ in Equation (1) can be obtained as follows:

$$\begin{array}{c}{\displaystyle W_c=\frac{p_{amb}\pi R_c^3·{\Phi }_c\left(p_{im},{\omega }_t\right)}{R_a{T}_{amb}}{\omega }_t,W_{egr}=\frac{A_{egr}\left(u_{egr}\right)·{\Psi }_{egr}\left(p_{im},p_{em}\right)}{\sqrt{T_{em}{R}_e}}{p}_{em}}\hfill \\ {\displaystyle W_{ei}=\frac{{\eta }_{vol}\left(p_{im},n_e\right)·n_e{V}_d}{120R_a{T}_{im}}{p}_{im},W_t=\frac{A_{vgtmax}·f_{{\Pi }_t}\left(p_{em}\right)·f_{vgt}\left(u_{vgt}\right)}{\sqrt{T_{em}{R}_e}}{p}_{em}}\hfill \\ {\displaystyle W_f=\frac{{10}^{-6}}{120}{n}_{cyl}{n}_e{u}_{\delta },P_c=\frac{W_c\left(p_{im},{\omega }_t\right)·c_{pa}{T}_{amb}·\left({\Pi }_c^{1-1/\phantom{1{\gamma }_a}{\gamma }_a}-1\right)}{{\eta }_c\left(p_{im},{\omega }_t\right)}}\hfill \\ {\displaystyle P_t{\eta }_m={\eta }_{tm}\left(p_{em},{\omega }_t\right)·W_t\left(p_{em},u_{vgt}\right)·c_{pe}{T}_{em}·\left(1-{\Pi }_t^{1-1/\phantom{1{\gamma }_e}{\gamma }_e}\right)}\hfill \end{array}$$

(2)

However, it is difficult to estimate the measurement error of the MAF sensor based on the complicated nonlinear model Equation (1). In order to simply present the state space equation and the error estimation, define variables:

$$\begin{array}{c}{\rho }_1=\frac{p_{amb}\pi R_c^3·{\Phi }_c\left(p_{im},{\omega }_t\right)}{R_a{T}_{amb}},{\rho }_2=\frac{A_{egr}\left(u_{egr}\right)·{\Psi }_{egr}\left(p_{im},p_{em}\right)}{\sqrt{T_{em}{R}_e}}\hfill \\ {\rho }_3=\frac{{\eta }_{vol}\left(p_{im},n_e\right)·n_e{V}_d}{120R_a{T}_{im}},{\rho }_4=\frac{A_{vgtmax}·f_{{\Pi }_t}\left(p_{em}\right)·f_{vgt}\left(u_{vgt}\right)}{\sqrt{T_{em}{R}_e}}\hfill \\ {\rho }_5=\frac{c_{pa}{p}_{amb}\pi R_c^3·{\Phi }_c\left(p_{im},{\omega }_t\right)·\left({\Pi }_c^{1-1/\phantom{1{\gamma }_a}{\gamma }_a}-1\right)}{J_t{\omega }_t{R}_a·{\eta }_c\left(p_{im},{\omega }_t\right)}\hfill \end{array}$$

(3)

According to Equation (3), the variables ${\rho }_i(i=1,2,3,4,5)$ are available in real-time since $p_{im}$, $p_{em}$, ${\omega }_t$, $u_{egr}$, $u_{vgt}$, $u_{\delta }$ and $n_e$ can be measured or estimated online. Therefore, the nonlinear model Equation (1) can be cast into an LPV system:

$$\begin{array}{c}\hfill \dot{x}=A\left(\rho \right)x+E\end{array}$$

(4)

where:

$$\begin{array}{c}{\displaystyle \rho =\left(\begin{array}{ccccc}{\rho }_1& {\rho }_2& {\rho }_3& {\rho }_4& {\rho }_5\end{array}\right)}\hfill \\ x=\left(\begin{array}{c}{p}_{im}\\ p_{em}\\ {\omega }_t\end{array}\right),A\left(\rho \right)=\left(\begin{array}{ccc}-a_1{\rho }_3& a_1{\rho }_2& a_1{\rho }_1\\ a_2{\rho }_3& -a_2{\rho }_2-a_2{\rho }_4& 0\\ 0& 0& -{\rho }_5\end{array}\right)\hfill \\ E=\left(\begin{array}{c}0\\ a_2{W}_f\\ P_t{\eta }_m/\phantom{P_t{\eta }_m{J}_t{\omega }_t}{J}_t{\omega }_t\end{array}\right),a_1=\frac{R_a{T}_{im}}{V_{im}},a_2=\frac{R_e{T}_{em}}{V_{em}}\hfill \end{array}$$

(5)

In order to determine the bounds on the parameter vector ρ, a simulation study is performed using a 1.9 L four-cylinder common rail turbo diesel engine of enDYNA provided by Tesis [17,18]. The bounds of the parameter vector ρ are found using the simulation data from enDYNA over the European Transient Cycle (ETC), Federal Test Procedure 75 (FTP75) and New European Drive Cycle (NEDC) [19,20,21]. Then, the results are listed in Table 1. It follows that each parameter ${\rho }_i$ from parameter vector ρ is bounded by a minimum and maximum value ${\underline{\rho }}_i$ and ${\overline{\rho }}_i$.

Table 1. Bounds on the parameter vector ρ under three conditions. ETC, European Transient Cycle; FTP75, Federal Test Procedure 75; NEDC, New European Drive Cycle.

Parameter	ETC	FTP75	NEDC
${\rho }_1$	$\left[19,150\right]\times {10}^{-7}$	$\left[146,110\right]\times {10}^{-7}$	$\left[49,130\right]\times {10}^{-7}$
${\rho }_2$	$\left[-1.65,1.6\right]\times {10}^{-7}$	$\left[-0.56,1.15\right]\times {10}^{-7}$	$\left[-0.003,1.0\right]\times {10}^{-7}$
${\rho }_3$	$\left[1.09,6.4\right]\times {10}^{-7}$	$\left[1.18,6\right]\times {10}^{-7}$	$\left[1.1,5.4\right]\times {10}^{-7}$
${\rho }_4$	$\left[0.61,6.01\right]\times {10}^{-7}$	$\left[0.75,5.4\right]\times {10}^{-7}$	$\left[0.62,4.7\right]\times {10}^{-7}$
${\rho }_5$	$\left[0,3.15\right]$	$\left[0,2.78\right]$	$\left[0,2.67\right]$

Table 1. Bounds on the parameter vector ρ under three conditions. ETC, European Transient Cycle; FTP75, Federal Test Procedure 75; NEDC, New European Drive Cycle.

Parameter	ETC	FTP75	NEDC
${\rho }_1$	$\left[19,150\right]\times {10}^{-7}$	$\left[146,110\right]\times {10}^{-7}$	$\left[49,130\right]\times {10}^{-7}$
${\rho }_2$	$\left[-1.65,1.6\right]\times {10}^{-7}$	$\left[-0.56,1.15\right]\times {10}^{-7}$	$\left[-0.003,1.0\right]\times {10}^{-7}$
${\rho }_3$	$\left[1.09,6.4\right]\times {10}^{-7}$	$\left[1.18,6\right]\times {10}^{-7}$	$\left[1.1,5.4\right]\times {10}^{-7}$
${\rho }_4$	$\left[0.61,6.01\right]\times {10}^{-7}$	$\left[0.75,5.4\right]\times {10}^{-7}$	$\left[0.62,4.7\right]\times {10}^{-7}$
${\rho }_5$	$\left[0,3.15\right]$	$\left[0,2.78\right]$	$\left[0,2.67\right]$

2.2. 2D Map Description for the MAF Sensor Error

The intake manifold pressure $p_{im}$, turbine speed ${\omega }_t$ and compressor mass air flow $W_c$ are the outputs of interest to analyze the MAF sensor error, which is:

$$\begin{array}{c}\hfill y=\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)=\left(\begin{array}{c}{p}_{im}\\ {\omega }_t\\ W_c\end{array}\right)\end{array}$$

(6)

where $y_1={\left(\begin{array}{cc}{p}_{im}& {\omega }_t\end{array}\right)}^{\mathrm{T}},y_2=W_c$. Due to the existence of MAF sensor error, the output Equation (6) becomes:

$$\begin{array}{c}\hfill y_m=\left(\begin{array}{c}{y}_{1m}\\ y_{2m}\end{array}\right)=\left(\begin{array}{c}{y}_1\\ y_2+\Delta W_c\end{array}\right)\end{array}$$

(7)

where $y_m$ is the measured value from sensors. ${\Delta W}_c$ is the measurement error of the MAF sensor, which depends on the engine operating point (fuel mass injection quantity $u_{\delta }$ and engine speed $n_e$), i.e., $\Delta W_c\left(u_{\delta }.n_e\right)$. Since it is difficult to accurately build an analytical model for $\Delta W_c\left(u_{\delta }.n_e\right)$, a 2D map is adopted in this paper to describe $\Delta W_c\left(u_{\delta },n_e\right)$. Therefore, define the partition of the 2D map input $\upsilon =\left(u_{\delta },n_e\right)$ as:

$$\begin{array}{c}\hfill a=u_{\delta }^1<u_{\delta }^2<\dots <u_{\delta }^{p_1}=b\\ c=n_e^1<n_e^2<\dots <n_e^{p_2}=d\end{array}$$

(8)

where $a,b\in \mathbb{R}$ are the minimum and maximum values of $u_{\delta }$ and $p_1$ is the number of the grid points in $\left[a,b\right]$. $c,d\in \mathbb{R}$ are the minimum and maximum values of $n_e$, and $p_2$ is the number of the grid points in $\left[c,d\right]$.

Assume that the measurement error of the input grid points $\left(u_{\delta }^i,n_e^j\right)$ is ${\theta }^{i,j}$, i.e.,

$$\begin{array}{c}{\theta }^{i,j}=\Delta W_c\left(u_{\delta }^i,n_e^j\right)\\ i=1,2,\cdots p_1;j=1,2,\cdots p_2\end{array}$$

(9)

Then, for $\forall \upsilon \in \left[u_{\delta }^i,u_{\delta }^{i+1}\right]\times \left[n_e^j,n_e^{j+1}\right]$, $\forall i\in \left[1,2,\cdots ,p_1-1\right]$ and $\forall j\in \left[1,2,\cdots ,p_2-1\right]$, we can hold the $n_e$ value fixed and apply one dimensional (1D) linear interpolation in the $u_{\delta }$ direction. Using the Lagrange form, the result is:

$$\begin{array}{c}{q}_{u_{\delta }}^j\left({\theta }^{i,j}\right)=\frac{u_{\delta }^{i+1}-u_{\delta }}{u_{\delta }^{i+1}-u_{\delta }^i}{\theta }^{i,j}+\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i}{\theta }^{i+1,j}\\ u_{\delta }\in \left[u_{\delta }^i,u_{\delta }^{i+1}\right],i=1,2,\cdots p_1-1\end{array}$$

(10)

Equation (10) can then be used to linearly interpolate along the $n_e$ dimension to yield the piecewise bilinear interpolation model of the measurement error $\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)$ as:

$$\begin{array}{c}\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)=\frac{n_e^{j+1}-n_e}{n_e^{j+1}-n_e^j}{q}_{u_{\delta }}^j\left({\theta }^{i,j}\right)+\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}{q}_{u_{\delta }}^{j+1}\left({\theta }^{i,j}\right)\\ \upsilon \in \left[u_{\delta }^i,u_{\delta }^{i+1}\right]\times \left[n_e^j,n_e^{j+1}\right],i=1,\cdots p_1-1;j=1,\cdots p_2-1\end{array}$$

(11)

For the undefined region $\upsilon \in \left(\mathbb{R}\times \mathbb{R}\right)\backslash \left(\left[a,b\right]\times \left[c,d\right]\right)$, we extend Equation (11) to the final result:

$$\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)=\left\{\begin{array}{cc}{q}_{u_{\delta }}^1\left({\theta }^{i,j}\right)& ,u\in \mathbb{R}\times {\mathbb{R}}_{n_e}^0\\ {\displaystyle \frac{n_e^{j+1}-n_e}{n_e^{j+1}-n_e^j}{q}_{u_{\delta }}^j\left({\theta }^{i,j}\right)+\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}{q}_{u_{\delta }}^{j+1}\left({\theta }^{i,j}\right)}& ,u\in \mathbb{R}\times {\mathbb{R}}_{n_e}^j\\ q_{u_{\delta }}^{p_2}\left({\theta }^{i,j}\right)& ,u\in \mathbb{R}\times {\mathbb{R}}_{n_e}^{p_2}\end{array}\right.$$

(12)

where:

$$\begin{array}{c}{q}_{u_{\delta }}^j\left({\theta }^{i,j}\right)=\left\{\begin{array}{cc}{\theta }^{1,j}& ,u_{\delta }\in {\mathbb{R}}_{u_{\delta }}^0\\ {\displaystyle \frac{u_{\delta }^{i+1}-u_{\delta }}{u_{\delta }^{i+1}-u_{\delta }^i}{\theta }^{i,j}+\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i}{\theta }^{i+1,j}}& ,u_{\delta }\in {\mathbb{R}}_{u_{\delta }}^i\\ {\theta }^{p_1,j}& ,u_{\delta }\in {\mathbb{R}}_{u_{\delta }}^{p_1}\end{array}\right.\\ i=1,2,\cdots p_1-1;j=1,2,\cdots p_2-1\end{array}$$

(13)

and:

$${\mathbb{R}}_{u_{\delta }}^k=\left\{\begin{array}{cc}\left(-\infty ,u_{\delta }^1\right],\hfill & k=0\hfill \\ \left(u_{\delta }^k,u_{\delta }^{k+1}\right],\hfill & k=1,\dots ,p_1-1\hfill \\ \left(u_{\delta }^{p_1},+\infty \right),\hfill & k=p_1\hfill \end{array}\right.$$

(14)

$${\mathbb{R}}_{n_e}^l=\left\{\begin{array}{cc}\left(-\infty ,n_e^1\right],\hfill & l=0\hfill \\ \left(n_e^l,n_e^{l+1}\right],\hfill & l=1,\dots ,p_2-1\hfill \\ \left(n_e^{p_2},+\infty \right),\hfill & l=p_2\hfill \end{array}\right.$$

(15)

For the purposes of estimating unknown parameter ${\theta }^{i,j}$ in $\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)$ expediently, Equation (12) in vector-vector form is needed. According to the input interval Equations (14) and (15), we define membership function as:

$$\begin{array}{c}{\delta }_{u_{\delta }}^k=\left\{\begin{array}{cc}1,& u_{\delta }\in {\mathbb{R}}_{u_{\delta }}^k\\ 0,& other\end{array}\right.\\ k=0,1,\cdots ,p_1\end{array}$$

(16)

and:

$$\begin{array}{c}{\delta }_{n_e}^l=\left\{\begin{array}{cc}1,& n_e\in {\mathbb{R}}_{n_e}^l\\ 0,& other\end{array}\right.\\ l=\mathrm{0,}1,\dots ,p_2\end{array}$$

(17)

Using membership function Equations (16) and (17), Equation (12) becomes:

$$\begin{array}{cc}\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)\hfill & ={\delta }_{n_e}^0{q}_{u_{\delta }}^1\left({\theta }^{i,j}\right)+\sum _{j=1}^{p_2-1}{\delta }_{n_e}^j\left(\frac{n_e^{j+1}-n_e}{n_e^{j+1}-n_e^j}{q}_{u_{\delta }}^j\left({\theta }^{i,j}\right)+\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}{q}_{u_{\delta }}^{j+1}\left({\theta }^{i,j}\right)\right)+{\delta }_{n_e}^{p_2}{q}_{u_{\delta }}^{p_2}\left({\theta }^{i,j}\right)\hfill \\ & ={\Psi }_{n_e}·q_{u_{\delta }}\left({\theta }^{i,j}\right)\hfill \end{array}$$

(18)

where:

$$\begin{array}{c}\hfill {\Psi }_{n_e}={\left(\begin{array}{c}{\displaystyle {\delta }_{n_e}^0+{\delta }_{n_e}^1\frac{n_e^2-n_e}{n_e^2-n_e^1}}\\ {\displaystyle {\delta }_{n_e}^1\frac{n_e-n_e^1}{n_e^2-n_e^1}+{\delta }_{n_e}^2\frac{n_e^3-n_e}{n_e^3-n_e^2}}\\ ⋮\\ {\displaystyle {\delta }_{n_e}^{p_1-2}\frac{n_e-n_e^{p_2-2}}{n_e^{p_2-1}-n_e^{p_2-2}}+{\delta }_{n_e}^{p_1-1}\frac{n_e^{p_2}-n_e}{n_e^{p_2}-n_e^{p_2-1}}}\\ {\displaystyle {\delta }_{n_e}^{p_1-1}\frac{n_e-n_e^{p_2-1}}{n_e^{p_2}-n_e^{p_2-1}}+{\delta }_{n_e}^{p_1}}\end{array}\right)}^{\mathrm{T}},q_{u_{\delta }}\left({\theta }^{i,j}\right)=\left(\begin{array}{c}{q}_{u_{\delta }}^1\left({\theta }^{i,j}\right)\\ q_{u_{\delta }}^2\left({\theta }^{i,j}\right)\\ ⋮\\ q_{u_{\delta }}^{p_2-1}\left({\theta }^{i,j}\right)\\ q_{u_{\delta }}^{p_2}\left({\theta }^{i,j}\right)\end{array}\right)\end{array}$$

(19)

and:

$$\begin{array}{cc}{q}_{u_{\delta }}^j\left({\theta }^{i,j}\right)\hfill & ={\delta }_{u_{\delta }}^0{\theta }^{1,j}+\sum _{i=1}^{p_1-1}{\delta }_{u_{\delta }}^i\left(\frac{u_{\delta }^{i+1}-u_{\delta }}{u_{\delta }^{i+1}-u_{\delta }^i}{\theta }^{i,j}+\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i}{\theta }^{i+1,j}\right)+{\delta }_{u_{\delta }}^{p_1}{\theta }^{p_1,j}\hfill \\ & ={\Psi }_{u_{\delta }}·{\theta }_{u_{\delta }}^j\hfill \end{array}$$

(20)

where:

$$\begin{array}{c}\hfill {\Psi }_{u_{\delta }}={\left(\begin{array}{c}{\displaystyle {\delta }_{u_{\delta }}^0+{\delta }_{u_{\delta }}^1\frac{u_{\delta }^2-u_{\delta }}{u_{\delta }^2-u_{\delta }^1}}\\ {\displaystyle {\delta }_{u_{\delta }}^1\frac{u_{\delta }-u_{\delta }^1}{u_{\delta }^2-u_{\delta }^1}+{\delta }_{u_{\delta }}^2\frac{u_{\delta }^3-u_{\delta }}{u_{\delta }^3-u_{\delta }^2}}\\ ⋮\\ {\displaystyle {\delta }_{u_{\delta }}^{p_1-2}\frac{u_{\delta }-u_{\delta }^{p_1-2}}{u_{\delta }^{p_1-1}-u_{\delta }^{p_1-2}}+{\delta }_{u_{\delta }}^{p_1-1}\frac{u_{\delta }^{p_1}-u_{\delta }}{u_{\delta }^{p_1}-u_{\delta }^{p_1-1}}}\\ {\displaystyle {\delta }_{u_{\delta }}^{p_1-1}\frac{u_{\delta }-u_{{}_{\delta }}^{p_1-1}}{u_{\delta }^{p_1}-u_{\delta }^{p_1-1}}+{\delta }_{u_{\delta }}^{p_1}}\end{array}\right)}^{\mathrm{T}},{\theta }_{u_{{}_{\delta }}}^j=\left(\begin{array}{c}{\theta }^{1,j}\\ {\theta }^{2,j}\\ ⋮\\ {\theta }^{p_1-1,j}\\ {\theta }^{p_1,j}\end{array}\right)\end{array}$$

(21)

and ${\Psi }_{n_e}\in {\mathbb{R}}^{1\times p_1},{\Psi }_{n_e}\in {\mathbb{R}}^{1\times p_2},{\theta }_{u_{{}_{\delta }}}^j\in {\mathbb{R}}^{p_1\times 1},q_{u_{\delta }}\left({\theta }^{i,j}\right)\in {\mathbb{R}}^{p_2\times 1}$.

Now, following Equations (18)–(21), $\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)$ can be written as a dot product between regression vector $\Psi \left(\upsilon \right)$ and unknown parameter vector θ as follows:

$$\begin{array}{c}\hfill \Delta W_{c,T}\left(\theta ,\upsilon \right)=\Psi \left(\upsilon \right)·\theta ,\forall \upsilon \in \mathbb{R}\times \mathbb{R}\end{array}$$

(22)

where:

$$\begin{array}{c}\hfill \Psi \left(\upsilon \right)={\Psi }_{n_e}·\left(\begin{array}{cccc}{\Psi }_{u_{\delta }}& 0& \cdots & 0\\ 0& {\Psi }_{u_{\delta }}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\Psi }_{u_{\delta }}\end{array}\right),\theta =\left(\begin{array}{c}{\theta }_{u_{\delta }}^1\\ {\theta }_{u_{\delta }}^2\\ ⋮\\ {\theta }_{u_{\delta }}^{p_2}\end{array}\right)\end{array}$$

(23)

and $\Psi \left(\upsilon \right)\in {\mathbb{R}}^{1\times p},\theta \in {\mathbb{R}}^{p\times 1},p=p_1·p_2$.

With the combination of Equations (4), (7) and (22), the diesel engine air path LPV model can be described by the following state space equation:

$$\begin{array}{cc}\hfill \dot{x}& =A\left(\rho \right)x+E\hfill \\ \hfill y_m& =C\left(\rho \right)x+G\Psi \left(\upsilon \right)\theta \hfill \end{array}$$

(24)

where:

$$\begin{array}{c}\hfill C\left(\rho \right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 0& {\rho }_1\end{array}\right),G=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\end{array}$$

(25)

Equation (24) indicates that the estimation of the MAF sensor error $\Delta W_c\left(u_{\delta },n_e\right)$ becomes joint estimation of state x and parameter θ for LPV system Equation (24).

3. Adaptive Observer Design

The observer to estimate state x and parameter θ jointly for the LPV system Equation (24) is given:

$$\begin{array}{cc}\dot{\widehat{x}}\hfill & =A\left(\rho \right)\widehat{x}+E+L\left(y_m-C\left(\rho \right)\widehat{x}-G\Psi \left(\upsilon \right)\widehat{\theta }\right)\hfill \\ \dot{\widehat{\theta }}\hfill & =\Gamma \Psi {\left(\upsilon \right)}^{\mathrm{T}}\left(y_{m2}-C_2\left(\rho \right)\widehat{x}-\Psi \left(\upsilon \right)\widehat{\theta }\right)\hfill \end{array}$$

(26)

where $C_2\left(\rho \right)=\left(\begin{array}{ccc}0& 0& {\rho }_1\end{array}\right)$, $\widehat{x}\in {\mathbb{R}}^{3\times 1}$ is the state estimate, $\widehat{\theta }\in {\mathbb{R}}^{p\times 1}$ is the parameter estimate, gain $\Gamma \in {\mathbb{R}}^{p\times p}$ is the positive definite diagonal matrix and $L\in {\mathbb{R}}^{3\times 3}$ is the feedback gain matrix.

The asymptotical stability of the proposed algorithm Equation (26) is analyzed in the following theorem.

Theorem 1. 

If the following Conditions (1) and (2) hold, then LPV adaptive observer Equation (26) is asymptotically stable, i.e., for any initial conditions $x\left(0\right),\widehat{x}\left(0\right),\widehat{\theta }\left(0\right)$ and parameter vector θ, the errors $\widehat{x}-x$ and $\widehat{\theta }-\theta$ tend to zero asymptotically when $t\to \infty$.

(1) There exist matrices L, $P=P^{\mathrm{T}}>0$, $Q=Q^{\mathrm{T}}>0$ and constant ${\epsilon }_1,{\epsilon }_2>0$, such that the following set of linear matrix inequalities (LMIs) is feasible for $\forall {\rho }_i\in \left[{\underline{\rho }}_i,{\overline{\rho }}_i\right]$, $i=1,2,3,4,5$:

$$M\left(\rho \right)=\left[\begin{array}{ccc}{A}_{cl}{\left(\rho \right)}^{\mathrm{T}}P+P{A}_{cl}\left(\rho \right)+Q& PLG& C_2^{\mathrm{T}}\left(\rho \right)\\ G^{\mathrm{T}}{L}^{\mathrm{T}}{P}^{\mathrm{T}}& -{\epsilon }_{1}I& 0\\ C_2\left(\rho \right)& 0& -{\epsilon }_{2}I\end{array}\right]<0$$

(27)

$$\begin{array}{c}\hfill 2-{\epsilon }_1-{\epsilon }_2>0\end{array}$$

(28)

where $A_{cl}\left(\rho \right)=A\left(\rho \right)-LC\left(\rho \right)$ .

(2) There exists map input υ, such that the regression vector $\Psi \left(\upsilon \right)$ is persistently exciting, i.e., $\exists {\delta }_1,{\delta }_2>0;\exists T>0;\forall t\ge 0$:

$$\begin{array}{c}\hfill {\delta }_1{I}_p\le {\int }_t^{t+T}\Psi {\left(\upsilon \left(\tau \right)\right)}^{\mathrm{T}}\Psi \left(\upsilon \left(\tau \right)\right)d\tau \le {\delta }_2{I}_p\end{array}$$

(29)

Proof. 

Set the estimation error $\tilde{x}=\widehat{x}-x,\tilde{\theta }=\widehat{\theta }-\theta$. Notice that $\dot{\theta }=0$; the error dynamic system between Equations (24) and (26) is:

$$\begin{array}{cc}\hfill \dot{\tilde{x}}& =\left(A\left(\rho \right)-LC\left(\rho \right)\right)\tilde{x}-LG\Psi \left(\upsilon \right)\tilde{\theta }\hfill \\ \hfill \dot{\tilde{\theta }}& =-\Gamma \Psi {\left(\upsilon \right)}^{\mathrm{T}}\left(C_2\left(\rho \right)\tilde{x}+\Psi \left(\upsilon \right)\tilde{\theta }\right)\hfill \end{array}$$

(30)

A valid Lyapunov function candidate is considered as $V={\eta }^{\mathrm{T}}P\eta +{\tilde{\theta }}^{\mathrm{T}}{\Gamma }^{-1}\tilde{\theta }$. For $\forall \tilde{x}\ne 0$, the derivative of V along with the error dynamic system Equation (30) is:

$$\begin{array}{cc}\hfill \dot{V}& =2{\tilde{x}}^{\mathrm{T}}P\dot{\tilde{x}}+2{\tilde{\theta }}^{\mathrm{T}}{\Gamma }^{-1}\dot{\tilde{\theta }}\hfill \\ & =2{\tilde{x}}^{\mathrm{T}}P{A}_{cl}\left(\rho \right)\tilde{x}-2{\tilde{\theta }}^{\mathrm{T}}\Psi {\left(\upsilon \right)}^{\mathrm{T}}\Psi \left(\upsilon \right)\tilde{\theta }-2{\tilde{x}}^{\mathrm{T}}PLG\Psi \left(\upsilon \right)\tilde{\theta }-2{\tilde{\theta }}^{\mathrm{T}}\Psi {\left(\upsilon \right)}^{\mathrm{T}}{C}_2\left(\rho \right)\tilde{x}\hfill \end{array}$$

(31)

There exist ${\epsilon }_1,{\epsilon }_2>0$, such that the following inequalities hold:

$$\begin{array}{c}-2{\tilde{x}}^{\mathrm{T}}PLG\Psi \left(\upsilon \right)\tilde{\theta }⩽{∥\frac{1}{\sqrt{{\epsilon }_1}}{\tilde{x}}^{\mathrm{T}}PLG∥}^2+{∥\sqrt{{\epsilon }_1}\Psi \left(\upsilon \right)\tilde{\theta }∥}^2\\ -2{\tilde{\theta }}^{\mathrm{T}}\Psi {\left(\upsilon \right)}^{\mathrm{T}}{C}_2\left(\rho \right)\tilde{x}⩽{∥\frac{1}{\sqrt{{\epsilon }_2}}{C}_2\left(\rho \right)\tilde{x}∥}^2+{∥\sqrt{{\epsilon }_2}{\tilde{\theta }}^{\mathrm{T}}\Psi {\left(\upsilon \right)}^{\mathrm{T}}∥}^2\end{array}$$

According to Condition (1) and Equation (31), the following inequality holds:

$$\begin{array}{cc}\hfill \dot{V}& \le {\tilde{x}}^{\mathrm{T}}\left(A_{cl}{\left(\rho \right)}^{\mathrm{T}}P+P{A}_{cl}\left(\rho \right)+\frac{1}{{\epsilon }_1}PLG{G}^{\mathrm{T}}{L}^{\mathrm{T}}P+\frac{1}{{\epsilon }_2}{C}_2^{\mathrm{T}}\left(\rho \right)C_2\left(\rho \right)\right)\tilde{x}\hfill \\ & -\left(2-{\epsilon }_1-{\epsilon }_2\right){\tilde{\theta }}^{\mathrm{T}}\Psi {\left(\upsilon \right)}^{\mathrm{T}}\Psi \left(\upsilon \right)\tilde{\theta }\hfill \\ & <-{\tilde{x}}^{\mathrm{T}}Q\tilde{x}<0\hfill \end{array}$$

That is $\dot{V}<-\omega \left(t\right)<0$, for $\forall {\rho }_i\in \left[{\underline{\rho }}_i,{\overline{\rho }}_i\right],i=1,2,3,4,5$, where $\omega \left(t\right)={\tilde{x}}^{\mathrm{T}}Q\tilde{x}$. Based on the Lyapunov stability theory, we know that the equilibrium $\tilde{x}=0$ and $\tilde{\theta }=0$ are stable. Now, integrating $\dot{V}<-\omega \left(t\right)$ from zero to t yields:

$$\begin{array}{c}\hfill V\left(t\right)+{\int }_0^t\omega \left(\tau \right)d\tau <V\left(0\right)\end{array}$$

(32)

and this means that ${\int }_0^t\omega \left(\tau \right)d\tau <V\left(0\right)$ since $V>0$. Therefore, we have $\underset{t\to \infty }{lim}{\int }_0^t\omega \left(\tau \right)d\tau \le V\left(0\right)$, and this implies that $\underset{t\to \infty }{lim}{\int }_0^t\omega \left(\tau \right)d\tau$ exists and is finite. By Barbalat’s Lemma [22], we know that $\underset{t\to \infty }{lim}\omega \left(t\right)=0$, and this leads to $\underset{t\to \infty }{lim}\tilde{x}\left(t\right)=0$.

Under Condition (2), the vector $\Psi \left(\upsilon \right)$ is persistently exciting, that is we have $\underset{t\to \infty }{lim}\tilde{\theta }\left(t\right)=0$ [22]. ☐

Remark 1. 

With the concept of multi-convexity [23], the solution of the infinite LMI Equation (27) can be reduced to be a solution of the finite LMIs for the vertex set, that is:

$$\begin{array}{c}M\left(w\right)<0,\forall w\in \mathcal{V}=\left\{\left.\left(w_1,\cdots ,w_5\right)\right|w_i\in \left\{{\underline{\rho }}_i,{\overline{\rho }}_i\right\},i=1,2,3,4,5\right\}\end{array}$$

(33)

Therefore, feedback gain L can be obtained by the solution of inequality Equations (28) and (33).

Remark 2. 

With the membership function ${\delta }_{u_{\delta }}^k,{\delta }_{n_e}^l$ in Equations (16) and (17), we know that ${\delta }_{u_{\delta }}^k$ = ${\delta }_{n_e}^l$ = 1 when $\upsilon \in {\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$ and ${\delta }_{u_{\delta }}^k={\delta }_{n_e}^l=0$ when $\upsilon \notin {\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$. Therefore, the regression vector $\Psi \left(\upsilon \right)$ is a sparse vector.

According to the partition of the map input $\upsilon =\left(u_{\delta },n_e\right)$ defined in Equation (8) and the piecewise bilinear interpolation model Equation (12), the input υ (engine operating point) moves in only one region ${\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$ at any time, and only the parameters ${\widehat{\theta }}^{i,j}$ corresponding to the region ${\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$ can participate in the interpolation. That is, for $\forall \upsilon \in {\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$:

Case 1: $\left(k,l\right)\in \left\{0,p_1\right\}\times \left\{0,p_2\right\}$. Only one parameter ${\widehat{\theta }}^{i,j},\left(i,j\right)\in \left\{1,p_1\right\}\times \left\{1,p_2\right\}$ takes part in the interpolation, i.e., $\Delta W_{c,T}\left({\widehat{\theta }}^{i,j},\upsilon \right)={\widehat{\theta }}^{i,j}$.

Case 2: $\left(k,l\right)\in \left\{1,2,\cdots p_1-1\right\}\times \left\{0,p_2\right\}$. Two parameters ${\widehat{\theta }}^{i,j},{\widehat{\theta }}^{i+1,j},\left(i,j\right)\in \left\{1,2,\cdots p_1-1\right\}\times \left\{1,p_2\right\}$ take part in the interpolation, i.e.,

$$\Delta W_{c,T}\left({\widehat{\theta }}^{i,j},\upsilon \right)=\frac{u_{\delta }^{i+1}-u_{\delta }}{u_{\delta }^{i+1}-u_{\delta }^i}{\widehat{\theta }}^{i,j}+\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i}{\widehat{\theta }}^{i+1,j}$$

Case 3: $\left(k,l\right)\in \left\{0,p_1\right\}\times \left\{1,2,\cdots p_2-1\right\}$. Two parameters ${\widehat{\theta }}^{i,j},{\widehat{\theta }}^{i,j+1},\left(i,j\right)\in \left\{1,p_1\right\}\times \left\{1,2,\cdots p_2-1\right\}$ take part in the interpolation, i.e.,

$$\Delta W_{c,T}\left({\theta }^{i,j},\upsilon \right)=\frac{n_e^{j+1}-n_e}{n_e^{j+1}-n_e^j}{\widehat{\theta }}^{i,j}+\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}{\widehat{\theta }}^{i,j+1}$$

Case 4: $\left(k,l\right)\in \left\{1,2,\cdots p_1-1\right\}\times \left\{1,2,\cdots p_2-1\right\}$. Four parameters ${\widehat{\theta }}^{i,j}$, ${\widehat{\theta }}^{i+1,j}$, ${\widehat{\theta }}^{i,j+1}$, ${\widehat{\theta }}^{i+1,j+1}$, $\left(i,j\right)\in \left\{1,2,\cdots p_1-1\right\}\times \left\{1,2,\cdots p_2-1\right\}$ take part in the interpolation, i.e.,

$$\begin{array}{cc}\hfill \Delta W_{c,T}\left({\widehat{\theta }}^{i,j},\upsilon \right)& =\frac{n_e^{j+1}-n_e}{n_e^{j+1}-n_e^j}\frac{u_{\delta }^{i+1}-u_{\delta }}{u_{\delta }^{i+1}-u_{\delta }^i}{\widehat{\theta }}^{i,j}+\frac{n_e^{j+1}-n_e}{n_e^{j+1}-n_e^j}\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i}{\widehat{\theta }}^{i+1,j}\hfill \\ & +\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}\frac{u_{\delta }^{i+1}-u_{\delta }}{u_{\delta }^{i+1}-u_{\delta }^i}{\widehat{\theta }}^{i,j+1}+\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i}{\widehat{\theta }}^{i+1,j+1}\hfill \end{array}$$

In order to expediently discuss the convergence of the parameter estimate ${\widehat{\theta }}^{i,j}$ corresponding to different regions ${\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$, a local regression vector ${\Psi }_l\left(\upsilon \right)$ is defined based on the above four classifications of the region partition as follow:

$${\Psi }_l\left(\upsilon \right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if} \upsilon \in \left({\mathbb{R}}_{u_{\delta }}^0\cup {\mathbb{R}}_{u_{\delta }}^{p_1}\right)\times \left({\mathbb{R}}_{n_e}^0\cup {\mathbb{R}}_{n_e}^{p_2}\right)\hfill \\ \left(\begin{array}{cc}1-{\eta }_1& {\eta }_1\end{array}\right)\hfill & \mathrm{if} \upsilon \in {\mathbb{R}}_{u_{\delta }}^i\times \left({\mathbb{R}}_{n_e}^0\cup {\mathbb{R}}_{n_e}^{p_2}\right)\hfill \\ \left(\begin{array}{cc}1-{\eta }_2& {\eta }_2\end{array}\right)\hfill & \mathrm{if} \upsilon \in \left({\mathbb{R}}_{u_{\delta }}^0\cup {\mathbb{R}}_{u_{\delta }}^{p_1}\right)\times {\mathbb{R}}_{n_e}^j\hfill \\ {\left(\begin{array}{c}\left(1-{\eta }_1\right)\left(1-{\eta }_2\right)\\ {\eta }_1\left(1-{\eta }_2\right)\\ \left(1-{\eta }_1\right){\eta }_2\\ {\eta }_1{\eta }_2\end{array}\right)}^{\mathrm{T}}\hfill & \mathrm{if} \upsilon \in {\mathbb{R}}_{u_{\delta }}^i\times {\mathbb{R}}_{n_e}^j\hfill \end{array}\right.$$

(34)

where:

$$\begin{array}{c}{\eta }_1=\frac{u_{\delta }-u_{\delta }^i}{u_{\delta }^{i+1}-u_{\delta }^i},{\eta }_2=\frac{n_e-n_e^j}{n_e^{j+1}-n_e^j}\\ i=1,2,\cdots p_1-1;j=1,2,\cdots p_2-1\end{array}$$

(35)

When $\upsilon \in {\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$, regression vector $\Psi \left(\upsilon \right)$ in Equation (26) can be replaced by local regression vector ${\Psi }_l\left(\upsilon \right)$; then, observer Equation (26) can be replaced by:

$$\begin{array}{cc}\hfill \dot{\widehat{x}}& =A\left(\rho \right)\widehat{x}+E+L\left(y_m-C\left(\rho \right)\widehat{x}-G{\Psi }_l\left(\upsilon \right){\widehat{\theta }}_l^{i,j}\right)\hfill \\ \hfill {\dot{\widehat{\theta }}}_l^{i,j}& ={\Gamma }_l{\Psi }_l{\left(\upsilon \right)}^{\mathrm{T}}\left(y_{m2}-C_2\left(\rho \right)\widehat{x}-{\Psi }_l\left(\upsilon \right){\widehat{\theta }}_l^{i,j}\right)\hfill \end{array}$$

(36)

where ${\widehat{\theta }}_l^{i,j}$ is the local parameter estimate of appropriate size and ${\Gamma }_l$ is a local positive definite diagonal matrix of appropriate size.

According to Theorem 1, the local parameter estimate ${\widehat{\theta }}_l^{i,j}$ is convergent if local regression vector ${\Psi }_l\left(\upsilon \right)$ is persistently exciting. Meanwhile, the parameter estimate $\widehat{\theta }$ is also convergent if the trajectory of the map input υ passes through all of the interpolation regions ${\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{u_2}^l$.

There are heavy matrices calculated in real time for the covariance matrix equation of EKF in [9,10,11], preventing it from being implemented in commercial electronic control units (ECUs) for map adaptation. Nevertheless, the computational burden of the proposed observer Equation (26) without the additional matrix equation is lower. Moreover, the number of parameter estimates $\widehat{\theta }$ updated in Equation (26) is no more than four at any time; then, the computational load can be further reduced by stopping estimating ${\widehat{\theta }}^{i,j}$ corresponding to $\upsilon \notin {\mathbb{R}}_{u_{\delta }}^k\times {\mathbb{R}}_{n_e}^l$.

Remark 3. 

For the area S where the trajectory of the map input υ does not move, the parameters ${\widehat{\theta }}^{i,j}$ corresponding to the interpolation region belonging to S cannot be estimated by observer Equation (26). In order to get the map parameters corresponding to S, an extrapolation model can be taken as follows:

$$\begin{array}{c}\hfill \Delta W_{c,e}\left(u_{\delta },n_e\right)=a_2{u}_{\delta }^2+a_1{u}_{\delta }+b_2{n}_e^2+b_1{n}_e+c_2{u}_{\delta }{n}_e+c_1\end{array}$$

(37)

where $a_2,a_1,b_2,b_1,c_2,c_1$ are polynomial parameters. Based on the data from the estimated map parameters, extrapolation model Equation (37) can be fitted by polynomial fitting approach, and then map parameters corresponding to S can be obtained.

4. Simulation Results

In this section, the simulation study of 2D map estimation is presented in the environment of a 1.9 L four-cylinder common rail turbo diesel engine of enDYNA, in which the ETC and FTP75 are used as test conditions, respectively. The observer architecture is illustrated in Figure 2, where $\Delta W_c\left(u_{\delta },n_e\right)$ is the additive reference error as the true measurement error from enDYNA.

Figure 2. Schematic diagram of the adaptive observer.

Figure 2. Schematic diagram of the adaptive observer.

Bounds on the parameter vector ρ are presented in Table 1. When the inequality Equations (28) and (33) are solved with ${\epsilon }_1=0.25$ and ${\epsilon }_2=0.11$, the gain matrix L can be given by:

$$\begin{array}{c}\hfill L=\left(\begin{array}{ccc}\text{1635.83}& -\text{124.95}& \text{6.52}\times {\text{10}}^{-7}\\ -\text{16,063.4}& \text{3135.43}& \text{1.67}\times {\text{10}}^{-7}\\ \text{4.2}& \text{29.7}& \text{4.55}\times {\text{10}}^{-20}\end{array}\right)\end{array}$$

(38)

The initial values of observer Equation (26) used in the simulation are $\widehat{x}\left(0\right)={\left[\begin{array}{ccc}9.8\times {10}^5& 9.8\times {10}^5& 0\end{array}\right]}^{\mathrm{T}},\widehat{\theta }\left(0\right)=0$, and the parameter gain is $\Gamma =200I$. Here, the reference error $\Delta W_c\left(u_{\delta },n_e\right)$ assumed as MAF sensor measurement error is depicted in Figure 3, which is superimposed on the signal $W_c$ in enDYNA as the measured value $y_{m2}$ in the simulation.

Figure 3. Reference error $\Delta W_c\left(u_{\delta },n_e\right)$ used as the mass air flow (MAF) sensor error in enDYNA.

Figure 3. Reference error $\Delta W_c\left(u_{\delta },n_e\right)$ used as the mass air flow (MAF) sensor error in enDYNA.

4.1. 2D Map Estimation under ETC

There are three parts of the ETC representing three different driving conditions, including urban, rural and motorway driving. Due to the engine speed range from ETC Part 1 covering the other two parts, ETC Part 1 is employed as the test condition in this section. Accordingly, the fuel mass injection quantity $u_{\delta }$ and engine speed $n_e$ from ETC Part 1 are plotted in Figure 4a, and the trajectory of the operating point $\upsilon =\left(u_{\delta },n_e\right)$ is depicted in Figure 4b, in which the trajectory does not move in area $S=\left[\left(0,56\right)\times \left(0,800\right)\right]\cup \left[\left(40,56\right)\times \left(0,2000\right)\right]$. According to the range $u_{\delta }\in \left[0,56\right]$ and $n_e\in \left[0,3100\right]$ from Figure 4a, an average partition can be respectively given as [0:4:56] and [0:250:3000].

Figure 4. Evolution of operating point $\upsilon =\left(u_{\delta },n_e\right)$ during ETC part one. (a) Evolution of $u_{\delta }$ and $n_e$. (b) Trajectory of operating point $\upsilon =\left(u_{\delta },n_e\right)$.

Figure 4. Evolution of operating point $\upsilon =\left(u_{\delta },n_e\right)$ during ETC part one. (a) Evolution of $u_{\delta }$ and $n_e$. (b) Trajectory of operating point $\upsilon =\left(u_{\delta },n_e\right)$.

The estimation results of the MAF sensor error using observer Equation (26) under ETC Part 1 are shown in Figure 5a, in which the map parameters have been estimated, except area S. According to Remark 3, the map parameters corresponding to area S can be obtained from the extrapolation model Equation (37). Based on the estimated map parameters from Figure 5a, the polynomial parameters in Equation (37) can be fitted as follows:

$$\begin{array}{c}{a}_2=9.12\times 1{\text{0}}^{-6},a_1=-4.39\times 1{\text{0}}^{-4},b_2=1.30\times 1{\text{0}}^{-9}\hfill \\ b_1=-5.74\times 1{\text{0}}^{-6},c_2=-5.37\times 1{\text{0}}^{-8},c_1=0.0012887\hfill \end{array}$$

(39)

The parameters corresponding to area S obtained from Equation (37) are presented in Figure 5b, which can roughly reflect the trend of the map.

In order to evaluate the accuracy of the estimated 2D map shown in Figure 5a, the comparison between the reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map during the ETC segment is presented in Figure 6a. Accordingly, the true mass air flow $y_2$, measured mass air flow $y_{m2}$ and map compensation are shown in Figure 6b. The mean relative error between reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and estimated 2D map is $10.41\%$, which demonstrates that the measured output $y_{m2}$ of the MAF sensor after map correction can approximate the true value of $W_c$ acceptably.

Figure 5. The estimated 2D map for the MAF sensor error. (a) Estimation results of the 2D map under ETC Part 1; (b) Extrapolation results based on the estimated 2D map.

Figure 5. The estimated 2D map for the MAF sensor error. (a) Estimation results of the 2D map under ETC Part 1; (b) Extrapolation results based on the estimated 2D map.

Figure 6. 2D map compensation during the ETC segment. (a) Simulation results between reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map; (b) Simulation results of true mass air flow $y_2$, measured mass air flow $y_{m2}$ and the map compensation.

Figure 6. 2D map compensation during the ETC segment. (a) Simulation results between reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map; (b) Simulation results of true mass air flow $y_2$, measured mass air flow $y_{m2}$ and the map compensation.

4.2. 2D Map Estimation under FTP75

In order to verify the effectiveness of the proposed method under different conditions, the cold start transient phase of the FTP75 is used in this section. Accordingly, $u_{\delta }$ and $n_e$ are plotted in Figure 7a, and the trajectory of υ is depicted in Figure 7b, in which the trajectory does not move in area $S=\left[\left(0,56\right)\times \left(0,650\right)\right]\cup \left[\left(25,37\right)\times \left(0,2000\right)\right]\cup \left[\left(37,56\right)\times \left(0,3500\right)\right]$.

Figure 7. Evolution of operating point $\upsilon =\left(u_{\delta },n_e\right)$ FTP75 cold start transient phase. (a) Evolution of $u_{\delta }$ and $n_e$; (b) Trajectory of operating point $\upsilon =\left(u_{\delta },n_e\right)$.

Figure 7. Evolution of operating point $\upsilon =\left(u_{\delta },n_e\right)$ FTP75 cold start transient phase. (a) Evolution of $u_{\delta }$ and $n_e$; (b) Trajectory of operating point $\upsilon =\left(u_{\delta },n_e\right)$.

Figure 8. The estimated 2D map for the MAF sensor error. (a) Estimation results of the 2D map under FTP75 cold start transient phase; (b) Extrapolation results based on the estimated 2D map.

Figure 8. The estimated 2D map for the MAF sensor error. (a) Estimation results of the 2D map under FTP75 cold start transient phase; (b) Extrapolation results based on the estimated 2D map.

The estimation results of the MAF sensor error under the cold start transient phase of the FTP75 are shown in Figure 8a, and the polynomial parameters in Equation (37) are fitted as follows:

$$\begin{array}{c}{a}_2=2.10\times 1{\text{0}}^{-5},a_1=-8.62\times 1{\text{0}}^{-4},b_2=1.46\times 1{\text{0}}^{-9}\hfill \\ b_1=-6.40\times 1{\text{0}}^{-6},c_2=-1.84\times 1{\text{0}}^{-8},c_1=0.0032316\hfill \end{array}$$

(40)

The map added the parameters corresponding to area S are shown in Figure 8b, which can also roughly reflect the trend of the map. Under the FTP75 segment, the comparison between the reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map from Figure 8a is shown in Figure 9a. Accordingly, the MAF sensor measured value $y_{m2}$ using map compensation is presented in Figure 9b. The mean relative error between reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map is $5.28\%$, demonstrating that the measured output $y_{m2}$ after map correction can approximate the true value of $W_c$ acceptably.

Figure 9. 2D map compensation during the FTP75 segment. (a) Simulation results between reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map; (b) Simulation results of true mass air flow $y_2$, measured mass air flow $y_{m2}$ and the map compensation.

Figure 9. 2D map compensation during the FTP75 segment. (a) Simulation results between reference error $\Delta W_c\left(u_{\delta },n_e\right)$ and the estimated 2D map; (b) Simulation results of true mass air flow $y_2$, measured mass air flow $y_{m2}$ and the map compensation.

5. Conclusions

A method for updating and storing sensor bias from different operating points is developed and investigated. This method achieves simultaneous estimation of model states and map parameters and applies to updating the MAF sensor error 2D map in the engine. The map in the form of a vector-vector dot product is given to conveniently analyze and design the parameter estimation method. An LPV adaptive observer to estimate map parameters is designed, which has the advantage of a simple structure and low computational load. Under ETC Part 1 and the cold start transient phase of the FTP75, the effectiveness of the presented algorithm is verified and validated in the engine software enDYNA. The results demonstrate that the proposed method can estimate the MAF sensor error acceptably.