Proceedings of the 17th World Congress

The International Federation of Automatic Control

Seoul, Korea, July 6-11, 2008

Grey-Box Control Oriented Emissions Models

Markus Hirsch* Daniel Alberer** Luigi del Re**

*LCM- Linz Center of Mechatronics, Linz, Austria

(Tel: 0043-732-2468-1842; email: markus.hirsch@lcm.at)
**Institute for Design and Control of Mechatronical Systems,
Johannes Kepler University, Linz, Austria (email: {daniel.alberer, luigi.delre}@jku.at)

Abstract: Further improvements of emission control will require reliable estimation of emissions in real
time. While many progresses are being done in terms of physical sensors, there is a wide agreement that
virtual sensors and more in general real time emission models will play a central role in the next steps.
While there is a deep understanding of the physics of the regulated pollutants, most general emission
models tend to be too complex and poorly parametrized to be used on-line, while most data based models
tend to be either insufficiently precise or of limited scope. To avoid this problem, this paper proposes a
combined approach in which static maps are identified numerically, but the effect of dominant factors, like
cylinder-head temperature and air path dynamics, is included on the basis of physical assumptions.
Differently from most models developed for sensors, this approach is based on pure engine control unit
(ECU) data, i.e. can be used for the computation of optimal control laws. As the paper shows, this strategy
is able to provide not only real time estimation of NOx as a function of the ECU outputs, but also of
particulate matter (PM).

theoretical models are very helpful in this case, their ability to

1. INTRODUCTION
simulate exhaust emissions accurately is limited.
The high mobility requirements have well known negative
Therefore, for on-line use, usually identified models are used
side effects, in particular in terms of environmental damage
instead. The variety of these models is huge, as the numbers
and air pollution by exhaust gas, as carbon dioxide (CO2),
of publications in the last years show. In principle they can be
carbon monoxide (CO), residual hydrocarbons (HC),
divided in two different groups: on the one hand models
particulate matter (PM) and nitrogen oxides (NOx, essentially
(Ouenou-Gamo et al. 1998; Schilling et al. 2006) whose
a mixture of NO and NO2). Accordingly, vehicles with
structure is defined by basic chemical and physical equations
combustion engines may be put on the market only after
(see: Heywood 1988; Warnatz et al. 1999) and only
passing standardized emission tests and require special after-
parameters are adapted and identified respectively to
treatment systems. In the case of Diesel engines, CO and HC
represent the emissions of a specific engine. In general, these
limits can be achieved rather easily with oxidation catalysts,
models can represent global behaviors well, but are often
while NOx as well as PM limits are much more demanding.
overstrained to simulate local, engine dependent, behaviors.
Although these limits can be met by Diesel particulate filters
Another issue for restriction at real engine applications can be
(DPF) and selective catalytic reduction (SCR), a
the use of characteristic values like the heat release rate, often
minimization of the engine row emissions is important to
not available, as input variables available in the ECU.
optimize or even avoid these exhaust aftertreatment devices.
On the other hand, so called black-box models assume no
The increasing number of possible control actions in the
specific model structure. Due the highly nonlinear behavior
combustion process of modern Diesel engines allows a much
of an engine and the huge interaction of the different inputs,
better optimization of the combustion, but complicates the
standard linear identification algorithms like FIR, ARX or
locating process for the optimal set points. The challenge
OE are even locally overstrained. For the identification of
here is finding models that represent this highly nonlinear
nonlinear systems many different possibilities exist (Ljung
process in a way suitable for optimization. Models of
1999, Nelles 2001), each of them having advantages and
combustion processes have been developed for a long time.
disadvantages. Very common algorithms to describe engine
As shown in (Patterson et al. 1994; Li et al. 1995; Wang et
emissions are artificial neural networks (ANN) (Hafner et al.
al. 1999; Jung et al. 2001), these models usually define
2000; Desantes et al. 2002; Ayeb et al. 2005; Galindo et al.
different interacting zones in terms of chemical compositions
2005). Models applying genetic algorithms (del Re et al.
and operation conditions in the combustion chamber. These
2005; Alberer et al. 2005) are also used and brought good
models are often included in finite-element programs and can
results for the standard operating mode (under closed loop
be very useful in finding optimal chamber geometrics by
action of the ECU).
numerical analysis of the combustion process. Though

978-1-1234-7890-2/08/$20.00 © 2008 IFAC 8514 10.3182/20080706-5-KR-1001.3426

17th IFAC World Congress (IFAC'08)
Seoul, Korea, July 6-11, 2008

This restriction – to closed loop operation under the action of Table 1. Input Variables
the ECU – is critical, as it strongly reduces the dimension of
the model space. Indeed, an ECU tends to produce fixed Input varible Unit Description
relationships between control variables (for instance, for a qMI mg/Inj fuel mass of main injection
single speed, temperature and demanded torque, there is one phiMI Deg crank angle position of main
single combination of several control variables, like rail injection before top dead centre
pressure, pre-injection time and amount, main injection time qPI mg/Inj fuel mass of pilot injection
and amount, etc. which have been optimized for the specific tPI µs time of pilot injection before main
operation). While this generally increases the precision of the injection
model derived under closed loop operation for the closed pRail bar fuel pressure in common rail
loop condition, it makes it mostly worthless for the general N rpm engine speed
operation and therefore for optimization, for which a MAP mbar manifold air pressure
sufficiently general model is needed. MAF kg/h manifold fresh air mass flow
This work is concerned in deriving such an open loop model Toil °C engine oil temperature
for a Diesel engine. As the engine cannot be operated without
an ECU, and some combinations of control inputs could even fuel injection
lead to an immediate damage of the engine, a real open loop
30 mg
model identification is not possible, but the problem can be 2200rpm 2700rpm
1300rpm 1800rpm
reduced using tools both at the measurement and at the 25mg 25mg 25mg 25mg
interpretation level: the model is designed using a
combination of a comprehensive steady state model (obtained 950rpm 1300rpm 1800rpm 2200rpm 2700rpm
modifying the set points of the controller so to cover the 15mg 15mg 15mg 15mg 15mg

whole range of combinations of control inputs) and statistical

950rpm 1300rpm 1800rpm 2200rpm 2700rpm
tools are used to reduce the correlation present in data and 5mg 5mg 5mg 5mg 5mg
thus improve the problem condition. 0 mg
engine speed

Additionally the model uses first-principle dynamic 830 rpm 3000 rpm

extensions to account for the effects of the temperature and of

transport phenomena in the air path. The resulting model is Fig. 1. Switching between local models defines the output of
able to represent with a sufficient precision both the NOx and the global model
PM values over the whole engine operating range using only
values available in the ECU. M

2. MODEL DESIGN
y glob = ∑ w ( N , qMI ) ⋅ y
r r
r =1
(1)
M
2.1 General issues
∑ w ( N , qMI ) = 1
r
r =1
The design of a control oriented emission model involves
addressing several aspects, in particular the strong (d) Compared to the uncertain and slow dynamics (delay
nonlinearity of the engine map, the strong correlation time and low pass filter) of the emission measurement
between many measurements and the essentially open devices, the fast injection dynamics can be neglected.
problem of the right formula structure of an NOx and Dynamics of the air path - which are much slower than
especially of a PM estimator. These problems are solved in those of the injection system - could be excluded as the
our case in the following way: output values of this system manifold air flow (MAF)
(a) a data set “as rich as possible” under open loop and boost pressure (MAP), both measured by the
conditions is determined experimentally standard ECU, were taken as inputs for the emission
model. As the air mass flow sensor is not mounted
(b) variable selection methods are used to determine the directly on the manifold but in the upstream of the turbo
input channels with the largest independent impact on charger, a second order filter was included here to
emissions (resulting in table 1) approximate the time lag (based on the results of
(Vierlinger 2005)). The resulting transfer function from
(c) the whole engine map is subdivided in local regions (as the measured value MAF to the real manifold air mass
shown in fig. 1) and the outputs of the global model flow MAF* can be written in the frequency domain as
yglobal (emission values for NOx and PM) are computed follows:
switching between outputs of local models yr, using a
linear combination in overlapping regions (1). 1
MAF * = MAF (2)
T s + 2T ξ s + 1
2 2

Here T defines the time constant of the system ξ the

damping and s the complex Laplace parameter.

8515
17th IFAC World Congress (IFAC'08)
Seoul, Korea, July 6-11, 2008

2.2 Local models normalized standard error of each regressor SE (θˆi ) was
The local models describing the global model are defined by calculated according to (7). For robustness, these values had
in parameters linear regression functions. Eq. (3) shows the to be smaller then a defined level which was decreased
basic structure of these functions: the logarithmic engine iteratively until validation error on a not in the identification
output (NOx and opacity respectively) is defined by a mean data included test data set raised significantly. Otherwise, the
value and a deviation depending on the standard deviation input channel was rejected and a new identification was done.
σ uu multiplied by a sum of input functions fi ( ∆u )
Σ kk
weighted by regressors θi. The logarithm is used to achieve SE (θˆk ) =
constant relative prediction errors over a larger output range θˆk
(error homoscedasticity) and therefore accurate results for −1
lower emissions levels too. Both the output and the input Σ = σˆ
2
⎡⎣φ T φ ⎤⎦ (7)
variables have been standardized by subtraction of the mean
values ( y , u ) and division by the square root of their 1 n

variance ( σ yy , σ uu ).
σˆ =
2

n− p
∑(y k
− yˆ k )
2

k =1

ln( yˆ ) = y + σ yy ⋅ ∆y 2.4 First principle extensions

p

∆y = ∑θ i
⋅ f i ( ∆u ) (3) High combustion chamber wall temperatures have an
increasing influence on NOx by the Zeldovich-formation
i =1
(Warnatz et al. 1999) and a decreasing one on soot. These
u−u
∆u = temperatures, which are normally not known, can be
σ uu approximated by the engine oil temperature Toil, also a
variable of a standard ECU. Actually, this value is a dynamic
∆u represents the vector of standardized input variables and p output of the engine system and therefore not settable but
the number of regressors. Several approaches for the describable only by quite complex definitions of the engine
generation of possible new input variables have been tested. which are depending on much more parameters than those
In conclusion, a polynomial attempt of second order is used here for the emission model. The very low excitation of
sufficient and therefore (3) can be rewritten in matrix Toil at the local identification areas causes huge standard
notation as follows: errors SE in the regressors of the channels containing this
value. Due to the selection criterion described before these
channels were rejected and temperature dependence would
∆y = θ ⋅ ∆u
T
(4) disappear. There could be many approaches to tackle this
problem, but a simple and efficient possibility to include
[ ]
T
∆u= ∆u1 " ∆u p ∆u1 ∆u1 ∆u1∆u 2 " ∆u p ∆u p effects of temperature proven to be multiplying the output of
the nonlinear model by a scaling factor representing
A criterion J to assess the quality of a given estimated model temperature differences between the main temperature of the
is given by the sum of the square residuals: identification data and current oil temperature. The weighting
function for this temperature dependence was – on the basis
n
of physical considerations (Warnatz et al. 1999) - an
J = ∑ ( yk − yˆ k )
2
(5)
k =1
exponential function. The growth factor which is positive for
NOx and negative for PM, is defined by separate warm-up
2.3 Variable selections models tests under constant speed and load tests (Fig. 3).

The estimation problem with the cost function defined by

eq. (5) can be solved by a standard least square method, i.e.

θˆ = (φ φ ) φ Y
T −1 T
(6)

where Y represents the data vector of the target variable and

φ the data matrix of the input vectors (3), albeit without the
convergence properties of the standard ARX due to the
correlation of the inputs. This problem can be reduced either
by well known regularization methods, or even better by
choosing the most representative and therefore least Fig. 2. Structure of the local models
correlated description basis. The procedure follows the basic
ideas of orthogonal forward selection (Henning 2004), and
essentially consists in looking for the sensitivity of the
estimation result to each candidate regressor. Therefore, the

8516
17th IFAC World Congress (IFAC'08)
Seoul, Korea, July 6-11, 2008

550
1500 measured signal
measured points
validation
fitted exponential function identification

500

1000

450
NOx [ppm]

NOx [ppm]
400
500

350

0
0 100 200 300 400 500 600
300
40 50 60 70 80 90 100 110 120 130 time [s]
Temperature [deg C]

Fig. 4. Local model for NOx for the working point 950 rpm
Fig. 3. Measured influence of oil temperature on emissions and 15 mg/cycle of main injection
under constant speed and load.

4. EXPERIMENTAL FACILITIES AND 35 measured signal

validation
IDENTIFICATION 30
identification

The engine used to conduct the experiments was a production 25

2 litre 4 cylinder BMW Diesel engine on an AVL dynamical

test bench. The engine parameters qMI, phiMI, qPI, tPI, 20
Opacity [%]

pRail, positions of the exhaust gas recirculation valve (EGR), 15

and that of the variable geometry turbocharger (VGT),

defining MAF and MAP, have been specified directly and 10

independently. Therefore the ECU was used to reach the 5

operating point, but was partly out of authority during
measurements. Also engine speed N, controlled by the torque 0

of the test bench dynamometer, was changed separately and

2800 2900 3000 3100 3200 3300 3400 3500 3600 3700
uncorrelated to the other inputs. Emission analyses were done time [s]
by the HORIBA-MEXA 7000 for NOx and the opacimeter
AVL439 as an indicator for PM respectively. Fig. 5. Local model for opacity (representing PM) for the
For each of the local models, measurements were done. Here, working point 2200 rpm and 15 mg/cycle of main injection.
emissions response on uncorrelated steps of the input
variables mentioned above were measured. As already
mentioned delay time varies heavily for different operating 6. RESULTS
points on the global working range (Fig. 1). In this case the
description with local models again shows an advantage, Validation of the global model was done by using the
namely adaptation of the delay time. After correction of this, identified global model as a virtual sensor during the NEDC
the static polynomial model (4) and the parameters of the standard driving cycle. Though no data of this cycle has been
transfer function (2) have been identified iteratively to included in the identification process at all, results of good
minimize the cost function J. Furthermore, elimination of quality could be reached as Figures 6 and 7 show. Here is to
single input channels of ∆u was done iteratively by mention that delay times of the measured signals have been
decreasing the limit of maximal allowed normalized standard compensated in order to enable a fair comparison.
error until residuals on a defined data interval, only used for
validation, began to increase significantly. Fig. 4 and 5 show
identification and validation for NOx and opacity at local
levels.

8517
17th IFAC World Congress (IFAC'08)
Seoul, Korea, July 6-11, 2008

⎛ n
⎞
800
measured ⎜ ∑(y k
− yˆ )k
2

⎟
= max ⎜ 1 − ,0⎟
k =1
simulated FIT
700
⎜⎜ ⎟⎟
n

∑(y − y)
2

(8)
600 ⎝ k =1
k
⎠
1 n
500
MSE = ∑(y − yˆ )
2
NOx [ppm]

k k

n k =1

400

Though especially the model for NOx shows good results,

300
difficulties can be seen in the high peaks where the virtual
200
sensor displays even higher values. These and also other
inaccuracies can be explained easily by the fact that here the
100 local models have to be extrapolated. Fast and high fuel
reduction rates from high levels of injections where high
720 740 760 780 800 820 840 860 880 900 values of MAF and especially MAP occur, cause fast changes
time [s] to models of lower mean injection. MAF and especially MAP
values are not decreasing that fast which means that the local
Fig. 6. Comparison between measurements and simulation model of low load has to be extrapolated in a wide range as in
result for NOx (selection of the NEDC) identification data neither points containing such high values
were included nor were possible. Polynomial functions of
low order show much better extrapolation abilities than for
measured example ANN-models but here even these are overstrained.
35 simulated To handle these effects more accurate, suitable regions of
transient identification cycles could be included into the
30
identification process of the local models.
25
These effects are apparently influencing the opacity model
Opacity [%]

20 much more. Though the validation at local models is quite

good, global validation seems to have room for improvement.
15 It is to mention here that especially soot model identification
was very difficult at low speed levels with higher injection.
10
Especially in regions where rail pressure is limited by low
5 engine speed, many effects are caused by local flame
instabilities which are hardly describable.
0
720 740 760 780 800 820 840 860 880 900 General speaking, identification of models at low load was
time [s] more difficult than doing this at higher load areas. This can
be explained by the fact that variation of inputs was hold
Fig. 7. Comparison between measurements and simulation constant for all models and therefore for instance a
result for opacity representing PM (selection of the NEDC) 3 mg/cycle fuel mass change at a mean injection level of
5 mg/cycle has relatively much more influence than doing
In table 2 two characteristic values for identification this at a mean level of 25 mg/cycle.
evaluation the FIT-value and the mean square error (MSE)
(8), are listed for global validation and local validation (mean 6. CONCLUSIONS
values).
The main message of this work is that our method is able to
Table 2. Model Errors provide a sufficient good model. The global structure of
switching between local models shows advantages as
mean mean global global different dynamics as well as sensitivities can be handled
Output local FIT local MSE FIT MSE easily. In combination with the variable selection criterion for
NOx 88% 34 ppm 84% 35ppm the local models, this structure shows a wide range of validity
Opac 75% 3.2 % opc 52% 5.7 % opc as well as high accuracy. Though improvements especially
for the PM model are possible, their necessity has to be
checked by using these for control optimization. Finally, it
could be shown that accurate engine emission models
containing only ECU values are possible and therefore
predestined as basis for control applications.

8518
17th IFAC World Congress (IFAC'08)
Seoul, Korea, July 6-11, 2008

7. ACKNOWLEDGEMENT Schilling, A., A. Amstutz, C. H. Onder and L. Guzzella

(2006) A Real-Time Model for the Prediction of the
This work was supported by the Linz Center of Competence
NOx Emissions in DI Diesel Engines, Proceedings of the
in Mechatronics (LCM) under grant LCM-001
2006 IEEE International Conference on Control
REFERENCES Applications.
Vierlinger, H. (2005) Verbesserung eines Frischluft-
Alberer, D., L. del Re, S. Winkler, P. Langthaler (2005)
massenmodells mit Hilfe von Umgebungs-größen an
Virtual Sensor Design of Particulate and Nitric Oxide in einem modernen Dieselmotor, Diploma Thesis –
a DI Diesel Engine, ICE 2005, paper Nr: 2005-24-063. Johannes Kepler University, Linz, Austria, Institute for
Ayeb, M., H. J. Theuerkauf and Thomas Winsel (2005) SI
Design and Control of Mechatronical Systems.
Engine Emissions Model Based on Dynamic Neural Wang, G., G. Li, Y. Liu, L. Chen, X. Zhang and J. Lu (1999),
Networks and D-Optimality, SAE 2005, Paper Nr. 2005- A Developed Model for Emission Prediction, SAE 1999,
01-0019 . Paper Nr. 1999-01-0233.
Del Re, L., P. Langthaler, C. Furmüller, S. Winkler, M. Warnatz, J., U. Maas and R.W. Dibble (1999), Combustion
Affenzeller (2005) NOx Virtual Sensor Based on
second Edition Springer-Verlag, Berlin Heidelberg.
Structure Identification and Global Optimization, SAE
2005, Paper Nr. 2005-01-0050.
Desantes, J. M., J. J. López, J. M. García and L. Hernández
(2002) Application of Neural Networks for Prediction
and Optimization of Exhaust Emissions in a H.D. Diesel
Engine, SAE 2002, Paper Nr. 2002-01-1144.
Galindo, J., J. M. Luján, J. R. Serrano, L. Hernández (2005),
Combustion simulation of turbocharger HSDI Diesel
engines during transient operation using neural networks,
Applied Thermal Engineering 25, 877-898.
Hafner, M., M. Schüller, O. Nelles, R. Isermann (2000) Fast
neural networks for diesel engine control design, Control
Engineering Practice 8, 1211-1221.
Hennig C. (2004). Modellwahl und Variablenselektion in der
Statistik [Online]. Available: http://www.math.uni-
hamburg.de/home/hennig/lehre/mskript1.pdf (September
10, 2007
Heywood, J. B. (1988). Internal Combustion Engine
Fundamentals, McGrwa-Hill series in mechanical
engineering, Singapore.
Jung, D. and D. Assanis (2001). Multi-Zone DI Diesel Spray
Combustion Model for Cycle Simulation Studies of
Engine Performance and Emissions, SAE 2001, Paper
Nr. 2001-01-1246 .
Li, X. and J. S. Wallace (1995). A Phenomenological Model
for Soot Formation and Oxidation in Direct-Injection
Diesel Engines, SAE 1995, Paper Nr. 952428.
Ljung, L. (1999) System Identification 2nd edition Prentice
Hall PTR, United States of America.
Nelles, O. (2001) Nonlinear System Identification, Springer-
Verlag, Berlin Heidelberg New York.
Ouenou-Gamo, S., M. Ouladsine and A. Rachid (1998)
Measurement and prediction of diesel engine exhaust
emissions, ISA Transactions 37, 135-140.
Patterson, M. A., S.-C. Kong, G. J. Hampson and R. D. Reitz
(1994). Modeling the Effects of Fuel Injection
Characteristics on Diesel Engine Soot and NOx
Emissions, SAE 1994, Paper Nr. 940523.

8519