# Cold-Start Modeling and On-Line Optimal Control of the Three-Way Catalyst

Jonathan Lock, Kristoffer Clasén, Jonas Sjöblom, Tomas McKelvey

###### Abstract

We present a three-way catalyst (TWC) cold-start model, calibrate the model
based on experimental data from multiple operating points, and use the model
to generate a Pareto-optimal cold-start controller suitable for implementation
in standard engine control unit hardware. The TWC model is an extension of a
previously presented physics-based model that predicts carbon monoxide,
hydrocarbon, and nitrogen oxides tailpipe emissions. The model axially and
radially resolves the temperatures in the monolith using very few state
variables, thus allowing for use with control-policy based optimal control
methods. In this paper we extend the model to allow for variable axial
discretization lengths, include the heat of reaction from hydrogen gas
generated from the combustion engine, and reformulate the model parameters to
be expressed in conventional units. We experimentally measured the temperature
and emission evolution for cold-starts with ten different engine load points,
which was subsequently used to tune the model parameters (e.g. chemical
reaction rates, specific heats, and thermal resistances). The simulated
cumulative tailpipe emission modeling error was found to be typically -20% to
+80% of the measured emissions. We have constructed and simulated the
performance of a Pareto-optimal controller using this model that balances fuel
efficiency and the cumulative emissions of each individual species. A
benchmark of the optimal controller with a conventional cold-start strategy
shows the potential for reducing the cold-start emissions.

##  I Introduction

The Three-Way Catalyst (TWC) is used in nearly all conventional vehicles with
spark-ignited (SI) engines to reduce the level of harmful emissions generated
by the combustion engine that would otherwise exit the tailpipe. The toxic
emissions generated by the combustion engine, broadly categorized as nitrogen
oxides (NO_x), carbon monoxide (CO), and residual hydrocarbons (THC), are in
the TWC converted to primarily form non-toxic nitrogen gas (N_2), carbon
dioxide (CO_2), and water (H_2O) [2, 3]. Modern TWCs are very effective at
removing emissions, with conversion efficiencies of over 95% being commonplace
[3] and in some cases significantly higher, as was found in the experimental
results in this paper. However, the TWC must be sufficiently hot to function,
which is in ordinary operation maintained by virtue of the hot exhaust gases
passing through it from the combustion engine. However, when a vehicle is
cold-started (i.e. started after the TWC has had sufficient time to cool to
the ambient temperature) the tailpipe emissions are much larger until the TWC
is heated to its ordinary operation temperature, an interval typically taking
on the order of 40–100 seconds[3]. These cold-start emissions are very
significant, and for many regulatory test procedures are responsible for
60–80% of the emissions generated from an entire test (which are for reference
on the order of 30 minutes).

Several methods for reducing cold-start emissions have been studied from a
multiple perspectives. These range from methods of constructing the TWC that
reduces the cold-start time [4], methods for preheating the TWC before
starting the combustion engine [5], and control schemes that focus on
controlling the combustion engine’s operation to limit the emissions generated
during the cold-start [3, 6, 7, 8, 9].

In this paper our goal is to develop a model-based optimal TWC cold-start
controller that can be feasibly be implemented in a standard engine control
unit (ECU). We will consider a conventional SI engine and TWC, where the
engine’s load point can be freely controlled during the cold-start, making the
controller suitable for e.g. hybrid vehicles. More specifically, we view the
TWC cold-start problem as determining the optimal engine speed, load, and
spark timing to apply over time while balancing the conflicting goals of
maximizing fuel efficiency and minimizing the cumulative emissions. Generating
a Pareto-optimal controller will therefore both require a dynamic thermal
model of the TWC as well as a suitable optimal control design method. This
paper can naturally be divided into two parts, one where we develop a TWC
cold-start model, and one where we evaluate the performance of an optimal
controller generated using said model.

The model presented in this paper extends on a TWC model previously developed
by the authors [1]. In the previous work, we derived a thermal model of the
TWC with very few state variables suited for fast off-line simulation or on-
line control systems. The model resolved both axial and radial temperature
variations from a first-principles perspective. In this paper we extend the
previous model by allowing the axial discretization to vary along the TWC’s
length, model the heat generation caused by oxidation of hydrogen in the
exhaust gas, and reformulate the tuning parameters (heat capacity, thermal
conductivity, etc) to be expressed in well-known units
($\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$,
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$, etc).
Furthermore, we expand on our previous work by here considering a TWC
consisting of two separate monoliths and use separate training and validation
datasets for tuning and evaluation.

In this paper we will first briefly discuss categories of existing models and
their strengths, weaknesses, and relevant applications. Following this we will
introduce our extended model and the experimental setup used to calibrate the
model. Finally, we will study the experimental results and evaluate the
simulated performance of the optimal control scheme using the calibrated
model. A listing of all abbreviations is shown in Table I, while all model
parameters and their units are presented in Table II.

TABLE I: Table of used abbreviations. CAbTDC | Crank Angle before Top Dead Center
---|---
CO | Carbon Monoxide
ECU | Engine Control Unit
SI | Spark Ignited
NO_x | Nitrogen Oxides
SA | Spark Angle
THC | Total Hydrocarbon
TWC | Three Way Catalyst
TWC1 | The first monolith in the TWC
TWC2 | The second monolith in the TWC
TABLE II: TWC parameters. Parameter | Unit |  Description
---|---|---
$L$ | $\mathrm{m}$ |  Total TWC length
$R$ | $\mathrm{m}$ |  TWC radius
$L_{n}$ | $\mathrm{m}$ |  Length of axial slice $n$
$L_{n-1,n}$ | $\mathrm{m}$ |  Center-to-center distance between axial slice $n-1$ and $n$
$t_{w}$ | $\mathrm{m}$ |  TWC wall thickness
$l_{c}$ | $\mathrm{m}$ |  TWC channel width
$\mathrm{OFA}$ | - |  TWC open frontal area
$m_{\mathrm{TWC}}$ | $\mathrm{kg}$ |  Monolith mass
$c_{p}$ | $\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$ |  Net monolith specific heat
$\dot{m}_{\mathrm{exh}}$ | $\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$ |  Exhaust massflow
$c_{p,\mathrm{exh}}$ | $\mathrm{J}\text{\,}{\mathrm{K}}^{-1}$ |  Exhaust gas specific heat
$T_{\mathrm{amb}}$ | $\mathrm{\SIUnitSymbolCelsius}$ |  Ambient temperature
$T_{\mathrm{exh}}$ | $\mathrm{\SIUnitSymbolCelsius}$ |  Temperature of gas feeding TWC
$k_{\mathrm{ax}}$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ |  Effective axial thermal conductivity
$k_{\mathrm{rad}}$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ |  Effective radial thermal conductivity
$k_{\mathrm{amb}}$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ |  Effective thermal conductivity to ambient
$t_{\mathrm{amb}}$ | $\mathrm{m}$ |  Thickness of insulation to ambient
$N$ | - |  Number of axial TWC segments
$M$ | - |  Number of resolved radial TWC channels
$T$ | $\mathrm{\SIUnitSymbolCelsius}$ |  Temperature state vector
$T_{n}$ | $\mathrm{\SIUnitSymbolCelsius}$ |  Central temperature in axial slice $n$
$\Delta_{T}$ | $\mathrm{\SIUnitSymbolCelsius}$ |  Center/periphery temperature difference
$k_{n,m}^{s}$ | $\mathrm{mol}\text{\,}{\mathrm{s}}^{-1}$ |  Reaction rate, species $s$ in cell $n,m$
$A^{s}$ | - |  Pre-exponential factor, species $s$
$E_{a}^{s}$ | $\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}$ |  Activation energy of emission species $s$
$y_{n,m}^{s}$ | - |  Mole fraction, species $s$ in cell $n,m$
$\dot{m}_{n,m}^{\mathrm{s,conv}}$ | $\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$ |  Mass conversion, species $s$, cell $n,m$
$\dot{m}_{n,m}^{\mathrm{s,in}}$ | $\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$ |  Incoming mass, species $s$, cell $n,m$
$\dot{m}_{n,m}^{\mathrm{s,conv}}$ | $\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$ |  Outgoing mass, species $s$, cell $n,m$
$\dot{m}^{s}_{\mathrm{tp}}$ | $\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$ |  Tailpipe massflow, species $s$
$t_{r,n}$ | $\mathrm{s}$ |  Gas residence time in slice $n$
$p_{\mathrm{TWC}}$ | $\mathrm{Pa}$ |  Absolute pressure in TWC
$P_{\mathrm{ctr}}$ | $\mathrm{W}$ |  Heating power, radial center
$P_{\mathrm{per}}$ | $\mathrm{W}$ |  Heating power, radial periphery
$P_{\mathrm{ax}}$ | $\mathrm{W}$ |  Heating power, axial conduction
$P_{\mathrm{rad}}$ | $\mathrm{W}$ |  Heating power, radial conduction
$P_{\mathrm{con,ctr}}$ | $\mathrm{W}$ |  Heating power, convection, center
$P_{\mathrm{con,per}}$ | $\mathrm{W}$ |  Heating power, convection, periphery
$P_{n,m}$ | $\mathrm{W}$ |  Exothermic power generated in cell $n,n$
$P_{\mathrm{exo,ctr}}$ | $\mathrm{W}$ |  Weighted central exothermic power
$P_{\mathrm{exo,per}}$ | $\mathrm{W}$ |  Weighted peripheral exothermic power
$P_{\mathrm{amb}}$ | $\mathrm{W}$ |  Heating power, loss to ambient
$q_{\mathrm{ax}}$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ |  Axial heat flux
$q_{\mathrm{rad}}$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ |  Radial heat flux
$q_{\mathrm{amb}}$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ |  Heat flux to ambient

###  I-A Literature survey

The topic of modeling the TWC for cold-start purposes has been considered from
a wide range of perspectives. Some authors derive fully three-dimensional or
two-dimensional models that capture many of the fundamentally complex chemical
kinetics, transport dynamics, and temperature dynamics [10, 11, 12]. Though
accurate, models of this caliber are very computationally demanding and
primarily suited for in-depth analysis.

One commonly used method of reducing the computational demand is to exploit
the characteristic structure of modern TWC’s, which consist of a large number
of parallel channels. Assuming each channel’s construction and composition is
identical it is sufficient to study the behavior of a single channel and
afterwards scale the result by the number of channels in the TWC. This
approximation is typically referred to as the single channel approximation or
single channel model, and results in models with significantly reduced
computational demands.

Several authors[13, 14, 15, 16, 17, 18, 19] model the TWC with the single
channel approximation, and depending on the number of modeled chemical
phenomenon and the level of spatial resolution some models can approach
realtime simulation speeds on powerful PC’s (i.e. where generating one
second’s worth of data requires one second of computation time). Beyond
modeling the temperature dependence of the TWC, many of these models also
include terms to capture a TWC phenomenon where oxygen is absorbed and
released in the TWC. The stored oxygen greatly influences the TWC’s capacity
to convert emissions, where CO and THC are more effectively oxidized when the
stored level of oxygen is large, while NO_x is more effectively reduced when
the level is low.

These models allow for simulating either the whole TWC or a representative
channel to a fairly high degree of accuracy. However, they are primarily of
use for analyzing the performance of the TWC in an off-line manner, for
instance in a TWC design process. In this paper we instead focus on on-line
TWC control, and in particular consider the cold-start problem. This requires
a model that is significantly simpler computationally, both as optimal control
methods place specific demands on the complexity and structure of the system
models, and as the ECU has a very limited computational capacity. There are
several classes of control-oriented models, ranging from models that
approximate the spatially varying temperature distribution as a scalar
temperature [20, 6, 7, 8], to more complex models that axially resolve the TWC
temperature and/or include oxygen-storage terms [21, 22, 23]. The simpler
models are fairly well suited for direct use with on-line optimal control
methods, while most of the more complex models are used primarily as a
starting point for creating a suboptimal controller.

The model presented in this paper (an extension of [1]) has been constructed
for the specific purpose of subsequently being implemented for optimal control
in conventional ECUs. More specifically, it is well-suited to control-policy
based optimal control methods, where the optimal control signal (e.g. engine
speed, load, spark angle, and so on) is precomputed in an offline phase and
stored in a table for a discrete set of TWC temperatures. A subsequent
realtime controller can ultimately determine the optimal control signal by
consulting the table of stored temperatures and associated optimal control
signals [24, 25]. Controllers of this class are very powerful, as they allow
for nearly arbitrarily nonlinear model dynamics, costs, and constraints, but
are limited in that their memory demand scales exponentially with the number
of state variables and require the states variables to either be measured or
estimated. In an effort to limit the number of state variables we have chosen
to not dynamically model the stored oxygen in the TWC as this is not as
significant as the temperature dynamics during a cold-start [20, 1].

##  II TWC Model

We will in this section introduce the TWC model. This model is based on and
extends a model previously presented by the authors [1]. The model is extended
by allowing the size of the axially discretized slices to vary over the length
of the catalyst, reformulated so all parameters are based on easily determined
physical parameters, and we consider the case where two separate TWC monoliths
are placed in series.

The model can naturally be divided into three distinct subsections; one
modeling the chemical kinetics, one modeling the temperature dynamics, and one
that interpolates the low-dimensional state variables to a higher-dimensional
temperature distribution using a physics-based method. We will initially
consider a single TWC monolith, and later return to the case where two are
placed in series.

(a) TWC body. Engine exhaust travels from left to right. Figure derived from
[1].

(b) Detail of TWC channels.

Figure 1: Illustration of geometry and measurement definitions for a single
TWC monolith.

An illustration of the assumed TWC geometry is shown in Figs. 1(a) and 1(b).
Specifically, we assume that the TWC is cylindrical with radius $R$ and length
$L$. We make the modeling choice of dividing the TWC into $N$ different axial
slices, and extend the previously presented model [1] that assumed equally-
sized slices by allowing the associated lengths $L_{1},L_{2},\dots,L_{N}$
(where $\sum L_{n}=L$) to be different for each slice. We also define the
lengths $L_{1,2},L_{2,3},\dots,L_{N-1,N}$ as the axial distances between the
midpoints of neighboring slices. Finally, we assume the TWC has a monolithic
structure with square, axially traversing channels of wall thickness $t_{w}$
and channel length $l_{c}$, as illustrated in Fig. 1(b).

Importantly, we make the approximation of axially discretizing the TWC
temperature. More specifically, we model the TWC temperature in the radial
center of slice $n$ as $T_{n}$. This implies that the temperature in the
radially central channel is modeled as $N$ segments of constant temperature.
Furthermore, we model the difference in temperature between the radial center
and radial periphery of each slice as as $\Delta_{T}$. Note that $\Delta_{T}$
is not axially resolved, i.e. we assume that $\Delta_{T}$ is identical for all
axial slices, i.e. the periphery temperature of slice 1 is $T_{1}+\Delta_{T}$,
slice 2 is $T_{2}+\Delta_{T}$, and so on. Finally, as is described in more
detail in Section II-A, we use an interpolation scheme to approximate the
temperature at $M$ different radial locations ranging from the radial center
to the periphery. Ultimately, a single TWC monolith is at any given instance
in time characterized by the state variable vector

$\displaystyle T$ $\displaystyle=\begin{bmatrix}T_{1}\\\ T_{2}\\\ \dots\\\
T_{N}\\\ \Delta_{T}\end{bmatrix}\,.$ (1)

In the following subsections we will detail the individual parts of the model.
First, we define how the state variable is used to generate a dense
representation of the TWC temperature that is axially and radially resolved.
This is followed by the chemical kinetics model that determines the conversion
efficiency of the TWC as well as the heating power generated by the exothermic
reactions. Finally, we introduce the thermal model that generates the state
variable derivative and define the interface between the two separate TWC
monoliths.

###  II-A Radial temperature interpolation

The model presented in this work uses the single-channel approximation for
fundamental TWC material properties, while resolving the radial temperature
profile by simulating several parallel channels corresponding to different
radial positions with different associated temperatures. This allows for
capturing the experimentally observed behavior where the periphery of the TWC
is significantly colder than the radially central sections (as we will be
discussed in Section IV-B).

More specifically, we model the radial temperature profile $\hat{T}(t,r)$, at
time $t$ and radius $r$, as a solution to the transient heat equation in a
flat circular disc with radius $R$ and an initial temperature of zero, i.e.
$\hat{T}(0,r)=0,r=[0,R]$. Note that $\hat{T}$ has no relation to the state
variable $T$ or the axial slice temperatures $T_{n}$ despite the similar
notation. Furthermore, we assume a Dirichlet boundary condition, i.e.
$\hat{T}(t,R)=0$, and assume that the plate develops a constant homogeneous
power. This power is intended to be analogous to the power delivered to a
slice in the TWC by convection, axial conduction, and exothermic heat
generation. Solving the time-evolution of $\hat{T}(t,r)$ is a textbook problem
(e.g. [26, p. 148]) with a solution that can be expressed as a Fourier-Bessel
series. Solving this numerically over time and radii can also be done easily
(for instance with MATLAB’s pdepe function), generating the radially-varying
time evolution of the plate’s temperature. Normalized solutions (where
$\hat{T}(t,r)$ is scaled to the range $[0,1]$, and $R=1$) are shown in Fig. 2.

Figure 2: Representative solutions to the transient heat equation in a flat
disc with a homogeneous power term and a Dirichlet boundary condition
$\hat{T}(R)=1$. The displayed temperatures and radii have been normalized to
the range $[0,1]$. Figure reused from [1].  Figure 3: Fully-resolved TWC, here
shown for $N=3$ and $M=4$. Figure reused from [1].

In this paper we interpolate the radial TWC temperature profile by using
precomputed solutions to the above flat-plate problem. More specifically, we
assume that the radial temperature profile in slice $n$ is given by
$A_{1}\hat{T}(t^{\prime},r)+A_{2}$ for a given time instant $t^{\prime}$.
Here, $A_{1}$ and $A_{2}$ are selected so that $\hat{T}(t^{\prime},0)=T_{n}$
and $\hat{T}(t^{\prime},R)=T_{n}+\Delta_{T}$ (i.e. $\hat{T}$ is scaled and
offset to match the known radial center and periphery temperatures), and
$t^{\prime}$ is selected to give a radial temperature profile that matches the
experimentally measured temperature profile for a given engine operating point
(more on this in Section IV-A). In summary, we interpolate the radial
temperature distribution using the known radially central and peripheral
temperatures and the instantaneous engine operating point. Letting $M$ denote
the number of independent single-channel models we wish to resolve, this
interpolation scheme allows us to convert the $N+1$ state variables to a
representation with $M$ single-channel models of $N$ segments. This is
illustrated in Fig. 3, where cells $([1\dots M],1)$ are fed with the incoming
gas, cell $(m,n+1)$ is fed with the output of cell $(m,n)$, and the output
from cells $([1\dots M],N)$ are combined to form the total exhaust from the
TWC.

As the TWC is assumed to be circular, regions near the periphery have a larger
associated area than regions near the axial center. This is taken into account
in our model by approximating the massflow in the physical TWC as equal at all
locations, and scaling the proportion of gas passing through each channel to
match. Letting $\dot{m}_{\mathrm{exh}}$ be the total exhaust massflow from the
engine and scaling by the relative area of an annular ring with a major radius
of $m/M$ and a minor radius of $(m-1)/M$ gives the massflow into a given cell
as

$\displaystyle\dot{m}_{m,1}$ $\displaystyle=\frac{\pi(m^{2}-(m-1)^{2})}{\pi
M^{2}}\dot{m}_{\mathrm{exh}}$
$\displaystyle=\frac{(2m-1)}{M^{2}}\dot{m}_{\mathrm{exh}}$ (2a)
$\displaystyle\dot{m}_{m,n+1}$ $\displaystyle=\dot{m}_{m,n}\,.$ (2b)

###  II-B Chemical kinetics

The total range of chemical reactions occurring in the TWC are very complex
and involve a wide range of compounds. However, there are fewer that
contribute to the legislated emissions or significant heat generation. We will
therefore limit our scope to net reactions (i.e. without considering
intermediary steps). This is done both for simplicity, and as a detailed
approach would require the addition of numerous state variables that track the
concentration of the emission species and their intermediaries in the TWC.
[19, 27, 1] give the most significant reactions (apart from Eq. 3b) as

$\displaystyle\ce{2CO+O_{2}}$ $\displaystyle\ce{->2CO_{2}}$ (3a)
$\displaystyle\ce{2H_{2}+O{2}}$ $\displaystyle\ce{->2H_{2}O}$ (3b)
$\displaystyle\ce{2NO+2CO}$ $\displaystyle\ce{->N_{2}+2CO_{2}}$ (3c)
$\displaystyle\ce{2NO_{2}}$ $\displaystyle\ce{->N_{2}+2O_{2}}$ (3d)
$\displaystyle\ce{2C_{3}H_{6}+9O_{2}}$ $\displaystyle\ce{->6CO_{2}+6H_{2}O}$
(3e) $\displaystyle\ce{C_{3}H_{8}+5O_{2}}$
$\displaystyle\ce{->3CO_{2}+4H_{2}O}\,.$ (3f)

We will assume all reactions are balanced, i.e. for a TWC with 100% conversion
efficiency all CO, H_2, C_3H_6, and C_3H_8 emitted from the engine are fully
oxidized and all NO and NO_2 are fully reduced.

In this paper we include the heat generated by the oxidation of hydrogen gas,
i.e. (3b), which is generated in the engine by the water-gas shift reaction.
As H_2 is not typically experimentally measured we will instead estimate its
mole-fraction from the measured mole-fractions of CO and CO_2, as given by [2,
eq. (4.68)]

$\displaystyle y_{\ce{H_{2}}}$
$\displaystyle=\frac{y_{\ce{H_{2}O}}y_{\ce{CO}}}{Ky_{\ce{CO_{2}}}}$ (4a) where
$\displaystyle y_{\ce{H_{2}O}}$
$\displaystyle=\frac{m^{\prime}}{2n^{\prime}}\frac{y_{\ce{CO}}+y_{\ce{CO_{2}}}}{1+y_{\ce{CO}}/(Ky_{\ce{CO_{2}}})+(m^{\prime}/2n^{\prime})(y_{\ce{CO}}+y_{\ce{CO_{2}}})}\,.$
(4b)

Here, $y$ corresponds to the mole-fraction of each corresponding compound, $K$
is a constant value set to 3.8 [2, eq (4.63)], and $n^{\prime}$ and
$m^{\prime}$ correspond to the number of carbon and hydrogen atoms
respectively111[2] indicates in (4.68) that $n^{\prime}$ and $m^{\prime}$
correspond to hydrogen and carbon respectively. Consulting the previous
derivations instead indicates that $n^{\prime}$ and $m^{\prime}$ correspond to
carbon and hydrogen. in each molecule of the fuel. With the RON95 E10 fuel
studied in this paper we used the supplier-specified value of
$m^{\prime}/2n^{\prime}=0.258$.

Typically, nitrogen oxides (NO and NO_2) and hydrocarbon (C_3H_6 and C_3H_8)
emissions are lumped together and denoted as NO_x and THC respectively [2, pp.
572–597]. We will in this paper assume a constant ratio of 99:1 for NO to NO_2
as indicated by [2, p. 578], and by [28, 29] a constant ratio of 3:1 for
C_3H_6 to C_3H_8.

We model the reaction rate $k_{n,m}^{s}$ of an emissions species $s$ in any
given cell $n,m$ using an Arrhenius expression of form

$\displaystyle k_{n,m}^{s}=A^{s}e^{\frac{-E_{a}^{s}}{RT_{n,m}}}\,,$ (5)

where $R$ is the ideal gas constant, $T_{n,m}$ is the temperature of cell
$n,m$, $E_{a}^{s}$ is the activation energy of emission species $s$, and
$A^{s}$ is the apparent pre-exponential factor for species $s$. Letting
$y_{n,m}^{s}$ indicate the mole fraction of emission species $s$ in cell
$n,m$, we model the evolution of the mole fraction as

$\displaystyle\frac{\mathrm{d}y_{n,m}^{s}}{\mathrm{d}t}$
$\displaystyle=-k_{n,m}^{s}y_{n,m}^{s}\,.$ (6)

Note that we do not include an inhibition factor in Eq. 6 in order to limit
the complexity of the model and the model tuning process. However, including
an inhibition term (e.g. as in [18]) is viable and would not interfere with
the optimal control method that will be described in Section V. Furthermore,
note that Eq. 6 does not include an oxygen concentration term. As the engine
is operated stoichiometrically the O_2 concentration is fairly constant,
implying that it can be lumped into $k_{n,m}^{s}$. This is beneficial as we
avoid the need to explicitly measure or model the O_2 concentration.

Though Eq. 4 allows for generating an estimate of the hydrogen gas
concentration for a given CO and CO_2 concentration (quantities which are
easily measured with conventional emissions-measurement equipment), as
$y^{\ce{H_{2}}}$ is not typically measured it is difficult to tune the
associated reaction rate parameters. By [18, 29] we have chosen to instead
model the reaction rate of Eq. 3b as identical to that of Eq. 3a, i.e. we
assume $E_{a}^{\ce{H_{2}}}=E_{a}^{\ce{CO}}$ and $A^{\ce{H_{2}}}=A^{\ce{CO}}$.
Finally, as H_2 is not typically viewed as a problematic emission species we
will in the remainder of this paper only consider Eq. 3b from the perspective
of determining the heat of reaction, in contrast to CO, THC, and NO_x
emissions which both contribute with their associated heat of reaction and
whose tailpipe emissions are important to track.

As described in more detail in [1], the gas residence time in each axial slice
is short enough for the monolith temperature to be close-to constant. Using
this constant-temperature approximation we can explicitly solve Eq. 6 as

$\displaystyle y_{n,m}^{s}(t_{r,n})=y_{n,m}^{s}(0)e^{-k_{n,m}^{s}t_{r,n}}$ (7)

for a residence time in slice $n$ of $t_{r,n}$. We will approximate the
residence time by assuming a plug-flow reactor model (i.e. assuming there is
no axial dispersion), giving

$\displaystyle t_{r,n}$ $\displaystyle=\frac{V_{\mathrm{slice},n}}{\nu}$ (8)

where $V_{\mathrm{slice},n}$ is the gas volume of slice $n$ and $\nu$ is the
volumetric flow-rate of the exhaust gases. Using the geometry of the TWC as
defined in Fig. 1, and using the ideal gas law we can approximate Eq. 8 as

$\displaystyle t_{r,n}$ $\displaystyle=\frac{\mathrm{OFA}\cdot L_{n}\pi
R^{2}}{\dot{m}_{\mathrm{exh}}R_{\mathrm{specific}}T_{n,m}p_{\mathrm{TWC}}^{-1}}\,,$
(9)

where $\mathrm{OFA}$ is the open frontal area of the TWC, defined by

$\displaystyle\mathrm{OFA}$ $\displaystyle=(l_{c}-t_{w})^{2}l_{c}^{-2}\,.$
(10)

In this paper we extend [1] to use physically meaningful SI units for all
parameters. Here $\dot{m}_{\mathrm{exh}}$ is the exhaust massflow
($\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$), $P$ is the absolute pressure
($\mathrm{Pa}$) in the TWC, which is typically close to the ambient pressure,
and $R_{\mathrm{specific}}$ is the specific gas constant
($\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$) for
tailpipe ratio of N_2, CO_2, O_2, and H_2O that was experimentally measured
for a hot TWC. This specific gas constant is used as it is easily determined
and the remaining gases only marginally contribute to $R_{\mathrm{specific}}$.

Ultimately, Eqs. 5 to 10 give a simple physics-based model of the most
significant reactions that occur in the TWC that takes temperature, gas
composition, and residence time into account.

Note that we have implicitly assumed that the incoming gas composition is
time-invariant, as this significantly reduces the number of required state
variables. (Explicitly modeling a time-varying incoming gas concentration can
require an additional $3N$ state variables, one for each emission species
concentration in each slice.) This implies that the model is suited for quasi-
static combustion engine operation, where the engine-out emission species and
massflow varies slowly with respect to the residence time in the TWC.
Fortunately, as the residence time in the entire TWC is fairly short (on the
order of 0.05 – 0.1 $\mathrm{s}$) [7, p. 64] we hypothesize that moderately-
varying dynamic operation with transitions on the order of 0.5 – 1
$\mathrm{s}$ will show accuracy similar to that of constant engine operation.

By Eq. 7, we can compute the massflow emitted from cell $n,m$ as

$\displaystyle\dot{m}_{n,m}^{s,\mathrm{out}}$
$\displaystyle=\dot{m}_{n,m}^{s,\mathrm{in}}e^{-k_{n,m}^{s}t_{r,n}}\,,$ (11)

and, by the conservation of mass, the converted massflow is trivially

$\displaystyle\dot{m}_{n,m}^{s,\mathrm{conv}}$
$\displaystyle=\dot{m}_{n,m}^{s,\mathrm{in}}-\dot{m}_{n,m}^{s,\mathrm{out}}\,.$
(12)

This lets us model the tailpipe emissions of emission species $s$ as the sum
of the outputs from each individual cell in the last axial segment, i.e.

$\displaystyle\dot{m}_{\mathrm{tp}}^{s}$
$\displaystyle=\sum_{m=1}^{M}\dot{m}_{m,N}^{s}\,.$ (13)

For convenience, we also define the conversion efficiency of the entire TWC
for a given emission species as

$\displaystyle\eta^{s}$
$\displaystyle=1-\frac{\dot{m}_{\mathrm{tp}}^{s}}{\dot{m}^{s}_{\mathrm{exh}}}\,,$
(14)

which we can view as the proportion of emissions converted in the TWC.

We generate an estimate of the exothermic reaction power generated by the
above reactions by computing the (temperature-dependent) heat of reaction for
each mole of reactant species as

$\displaystyle\mathrm{d}H_{\ce{CO}}$
$\displaystyle=H_{0,\ce{CO_{2}}}-H_{0,\ce{CO}}-1/2H_{0,\ce{O_{2}}}$ (15a)
$\displaystyle\mathrm{d}H_{\ce{H_{2}}}$
$\displaystyle=H_{0,\ce{H_{2}O}}-H_{0,\ce{H_{2}}}-1/2H_{0,\ce{O_{2}}}$ (15b)
$\displaystyle\mathrm{d}H_{\ce{NO}}$
$\displaystyle=1/2H_{0,\ce{N_{2}}}+H_{0,\ce{CO_{2}}}-H_{0,\ce{NO}}-\mathrm{d}H_{\ce{CO}}$
(15c) $\displaystyle\mathrm{d}H_{\ce{NO}_{2}}$
$\displaystyle=1/2H_{0,\ce{N_{2}}}+H_{0,\ce{O_{2}}}-H_{0,\ce{NO_{2}}}$ (15d)
$\displaystyle\mathrm{d}H_{\ce{C_{3}H_{6}}}$
$\displaystyle=3H_{0,\ce{CO_{2}}}+3H_{0,\ce{H_{2}O}}-H_{0,\ce{C3H6}}-9/2H_{0,\ce{O_{2}}}$
(15e) $\displaystyle\mathrm{d}H_{\ce{C_{3}H_{8}}}$
$\displaystyle=3H_{0,\ce{CO_{2}}}+4H_{0,\ce{H_{2}O}}-H_{0,\ce{C3H8}}-5H_{0,\ce{O_{2}}}\,.$
(15f)

For brevity, we have not explicitly stated the temperature dependence of the
above terms but include their temperature dependence in the numerical model.
We use the Shomate equation and reference constants given by the NIST
(available at https://webbook.nist.gov) to compute the numerical values of the
above terms. Using the previous concentration ratios for the lumped terms
gives the effective reaction power

$\displaystyle\mathrm{d}H_{\ce{NO_{x}}}$
$\displaystyle=(99\mathrm{d}H_{\ce{NO}}+\mathrm{d}H_{\ce{NO_{2}}})/100$ (16a)
$\displaystyle\mathrm{d}H_{\ce{THC}}$
$\displaystyle=(3\mathrm{d}H_{\ce{C_{3}H_{6}}}+\mathrm{d}H_{\ce{C_{3}H_{8}}})/4\,.$
(16b)

By Eqs. 15a, 15b, 16a and 16b the total temperature-dependent heat of reaction
generated in each cell is thus

$\displaystyle P_{n,m}$
$\displaystyle=\dot{m}_{n,m}^{\ce{CO},\text{conv}}\cdot\mathrm{d}H_{\ce{CO}}+\dot{m}_{n,m}^{\ce{H_{2}},\text{conv}}\cdot\mathrm{d}H_{\ce{H_{2}}}$
$\displaystyle+\dot{m}_{n,m}^{\ce{NO_{x}},\text{conv}}\cdot\mathrm{d}H_{\ce{NO_{x}}}+\dot{m}_{n,m}^{\ce{THC},\text{conv}}\cdot\mathrm{d}H_{\ce{THC}}\,.$
(17)

###  II-C Temperature dynamics

We model the temperature dynamics using a heat balance ODE. Introducing the
relative length-weighting matrix

$\displaystyle\mathbf{W}_{L}$
$\displaystyle=\begin{bmatrix}L_{1}/L&0&\dots&0\\\ 0&L_{2}/L&\dots&0\\\
\vdots&&\ddots&\vdots\\\ 0&0&\dots&L_{N}/L\\\ \end{bmatrix},$ (18)

we can define the heat-balance ODEs as

$\displaystyle
m_{\mathrm{TWC}}\mathbf{W}_{L}c_{p}\frac{\mathrm{d}T_{\mathrm{ctr}}}{\mathrm{d}t}$
$\displaystyle=P_{\mathrm{ctr}}$ (19a) $\displaystyle
m_{\mathrm{TWC}}\mathbf{W}_{L}c_{p}\frac{\mathrm{d}T_{\mathrm{per}}}{\mathrm{d}t}$
$\displaystyle=P_{\mathrm{per}}\,,$ (19b)

where

$\displaystyle P_{\mathrm{ctr}}$
$\displaystyle=P_{\mathrm{ax}}-P_{\mathrm{rad}}+P_{\mathrm{con,ctr}}+P_{\mathrm{exo,ctr}}$
(20a) $\displaystyle P_{\mathrm{per}}$
$\displaystyle=P_{\mathrm{ax}}+P_{\mathrm{rad}}+P_{\mathrm{con,per}}+P_{\mathrm{exo,per}}-P_{\mathrm{amb}}\,.$
(20b)

Here, $m_{\mathrm{TWC}}$ is the mass of the TWC, $c_{p}$ is its specific heat,
$P_{\mathrm{ctr}}$ and $P_{\mathrm{per}}$ are $N\times 1$ vectors
corresponding to the total power developed in the radial center and periphery
respectively, and $\frac{\mathrm{d}T_{\mathrm{ctr}}}{\mathrm{d}t}$ and
$\frac{\mathrm{d}T_{\mathrm{per}}}{\mathrm{d}t}$ are $N\times 1$ vectors
representing the temperature derivative in the radial center and periphery.
With these terms, we can construct the total state vector ODE as

$\displaystyle\frac{\mathrm{d}T}{\mathrm{d}t}$
$\displaystyle=\begin{bmatrix}\frac{\mathrm{d}T_{\mathrm{ctr}}}{\mathrm{d}t}\\\
N\mathrm{mean}(\mathbf{W}_{L}(\frac{\mathrm{d}T_{\mathrm{per}}}{\mathrm{d}t}-\frac{\mathrm{d}T_{\mathrm{ctr}}}{\mathrm{d}t}))\end{bmatrix}\,,$
(21)

where $\mathrm{mean}([x_{1},x_{2},\dots,x_{n}])$ corresponds to the arithmetic
mean of the elements in $x$, i.e. $1/n\sum_{i=1}^{n}x_{i}$. Note that we can
interpret
$N\mathrm{mean}(\mathbf{W}_{L}(\frac{\mathrm{d}T_{\mathrm{per}}}{\mathrm{d}t}-\frac{\mathrm{d}T_{\mathrm{ctr}}}{\mathrm{d}t}))$
as corresponding to the average difference between radially central and
peripheral powers, weighted by the relative length of each slice.

The power terms in the right hand side of Eqs. 20a and 20b are separated into
axial, radial, convection, exothermic, and ambient loss terms respectively,
which we will define below. This extends on our previous work [1], which used
lumped-element parameters without an explicit power-balance formulation.

####  II-C1 Axial conduction

Figure 4: Heat flux between two materials of known temperature.

The axial heat conduction power $P_{\mathrm{ax}}$ is modeled by Fourier’s heat
law [26]. Using conventional notation, the heat flux between two materials of
constant temperature separated by a material of thickness $l$ is in general

$\displaystyle q$ $\displaystyle=k\frac{T_{2}-T_{1}}{l},$ (22)

where $q$ is the heat flux ($\mathrm{W}\mathrm{/}\mathrm{m}^{2}$), $T_{1}$ and
$T_{2}$ are two known temperatures ($\mathrm{K}$), and $k$ is the thermal
conductivity of the material
($\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$), as
illustrated in Fig. 4.

In this paper, we extend our previous model [1] by modeling the axial flux
between successive axial slices as

$\displaystyle q_{\mathrm{ax}}$
$\displaystyle=k_{\mathrm{ax}}\begin{bmatrix}\dfrac{T_{2}-T_{1}}{L_{1,2}}\\\\[8.61108pt]
\dfrac{T_{3}-T_{2}}{L_{2,3}}\\\ \vdots\\\
\dfrac{T_{N}-T_{N-1}}{L_{N-1,N}}\end{bmatrix}$ (23)

where $k_{\mathrm{ax}}$ is the axial thermal conductivity and $L_{n,n+1}$ is
the distance between the center of axial slice $n$ and $n+1$, as illustrated
in Fig. 1(a). Note that for $N$ axial slices we thus have $N-1$ axial fluxes
between slices.

We model the power associated with each flux term by scaling by the surface
area of the solid mass of the TWC, i.e.
$q_{\mathrm{ax}}(1-\mathrm{OFA})R^{2}\pi$. We can then model the total
developed power in each axial slice due to conduction as the difference in
incoming and outgoing power fluxes, i.e.

$\displaystyle P_{\mathrm{ax}}=\begin{bmatrix}0\\\
q_{\mathrm{ax}}(1-\mathrm{OFA})R^{2}\pi\end{bmatrix}-\begin{bmatrix}q_{\mathrm{ax}}(1-\mathrm{OFA})R^{2}\pi\\\
0\end{bmatrix}.$ (24)

####  II-C2 Radial conduction

The radial heat conduction is modeled in a manner similar to the axial heat
conduction. The radial flux is modeled as

$\displaystyle q_{\mathrm{rad}}$
$\displaystyle=k_{\mathrm{rad}}\dfrac{\Delta_{T}}{R/2}\,,$ (25)

i.e. a temperature difference of $\Delta_{T}$ and separation of $R/2$.
Approximating the surface area conducting heat as that of a cylinder with half
the radius of the TWC and length equal to the TWC’s length gives a developed
radial conduction power of

$\displaystyle P_{\mathrm{rad}}$ $\displaystyle=q_{\mathrm{rad}}\pi
RL\mathbf{1}\,,$ (26)

where $\mathbf{1}$ is the ones vector of size $N\times 1$.

####  II-C3 Convection

The convection heat powers $P_{\mathrm{con,ctr}}$ and $P_{\mathrm{con,per}}$
are modeled under the assumption that each cell is sufficiently long and
narrow for the gas temperature to reach the cell temperature, i.e. the gas
travels slowly enough to reach thermal equilibrium with the TWC walls. This
was considered in our previous model [1], where we found that five axial
slices was a suitable upper limit. As we now also allow for the slice lengths
to vary, it is thus reasonable to require that $L_{n}\leq L/5\,\forall n$.

This gives the convection powers as

$\displaystyle P_{\mathrm{con,ctr}}$
$\displaystyle=\dot{m}_{\mathrm{exh}}c_{p,\mathrm{exh}}\begin{bmatrix}T_{\mathrm{exh}}-T_{1}\\\
T_{1}-T_{2}\\\ \vdots\\\ T_{N-1}-T_{N}\end{bmatrix}$ (27) $\displaystyle
P_{\mathrm{con,per}}$
$\displaystyle=\dot{m}_{\mathrm{exh}}c_{p,\mathrm{exh}}\begin{bmatrix}T_{\mathrm{exh}}-(T_{1}+\Delta_{T})\\\
(T_{1}+\Delta_{T})-(T_{2}+\Delta_{T})\\\ \vdots\\\
(T_{N-1}+\Delta_{T})-(T_{N}+\Delta_{T})\end{bmatrix}$ (28)

where $\dot{m}_{\mathrm{exh}}$ is the exhaust massflow
($\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$), $c_{p,\mathrm{exh}}$ is the
constant-pressure specific heat of the exhaust gases
($\mathrm{J}\text{\,}{\mathrm{kg}}^{-1}\text{\,}{\mathrm{K}}^{-1}$), and
$T_{\mathrm{exh}}$ is the temperature of the exhaust gas fed into the TWC
($\mathrm{K}$).

####  II-C4 Exothermic power

The exothermic power terms $P_{\mathrm{exo,ctr}}$ and $P_{\mathrm{exo,per}}$
are modeled by weighting the densely-resolved single-channel exothermic power
into an effective central and peripheral powers. Here we use a linear
weighting scheme as a first approximation, given as

$\displaystyle P_{\mathrm{exotherm,ctr}}$
$\displaystyle=\sum_{m=1}^{M}P_{n,m}(1-\frac{m-1}{M-1})$ (29) $\displaystyle
P_{\mathrm{exotherm,per}}$
$\displaystyle=\sum_{m=1}^{M}P_{n,m}\frac{m-1}{M-1}\,.$ (30)

Note that the term $\frac{m-1}{M-1}$ varies from 0 to 1 as $m$ varies from 1
to $M$.

####  II-C5 Ambient losses

The heat losses to the ambient environment are modeled as conductive, with a
total flux of

$\displaystyle q_{\mathrm{amb}}$
$\displaystyle=k_{\mathrm{amb}}\begin{bmatrix}(T_{1}+\Delta_{T}-T_{\mathrm{amb}})t_{\mathrm{amb}}^{-1}\\\
(T_{2}+\Delta_{T}-T_{\mathrm{amb}})t_{\mathrm{amb}}^{-1}\\\ \vdots\\\
(T_{N}+\Delta_{T}-T_{\mathrm{amb}})t_{\mathrm{amb}}^{-1}\,.\\\ \end{bmatrix}$
(31)

Here, $k_{\mathrm{amb}}$ is the effective thermal conductivity of the
insulating material and $t_{\mathrm{amb}}$ is its associated thickness.
Modeling the exposed surface area as a cylinder with radius and length equal
to the whole TWC (and thus neglecting heat loss through the circular ends of
the cylinder) gives the power loss to the ambient environment as

$\displaystyle P_{\mathrm{amb}}$ $\displaystyle=q_{\mathrm{amb}}2\pi RL\,.$
(32)

###  II-D Two-monolith structure

Figure 5: The series-connected TWC is modeled as two independent TWC’s with
the gas exiting the first being mixed and fed into the second.

In this paper we extend the model previously presented in [1] by studying a
TWC with series-coupled monoliths. More specifically, we consider two
physically separated TWC’s, where the gas leaving the first is assumed to be
completely mixed and then fed into the second, as illustrated in Fig. 5. We
model this by assuming two completely independent sets of TWC parameters (as
listed in Table II), and let the exhaust gas first travel from the engine
through the first TWC. As a first approximation, the gas leaving the first TWC
is assumed to be perfectly mixed (both with respect to temperature and
emission species concentrations) and then fed through the second TWC, which
then finally exits to the tailpipe.

Ideal mixing implies that the temperature of the gas feeding the second TWC is
given as

$\displaystyle T_{\mathrm{exh,TWC2}}$
$\displaystyle=\sum_{m=1}^{M_{\mathrm{TWC1}}}\frac{m^{2}-(m-1)^{2}}{M_{\mathrm{TWC1}}^{2}}T_{N,m,\mathrm{TWC1}}\,,$
(33)

where i.e. the gas leaving the first TWC (TWC1) is combined and scaled by the
its relative flow rate. Trivially, we also have that the emission species
concentration entering TWC2 is

$\displaystyle\dot{m}_{1,m,\mathrm{TWC2}}^{s,\mathrm{in}}$
$\displaystyle=\frac{m^{2}-(m-1)^{2}}{M_{\mathrm{TWC2}}^{2}}\dot{m}_{\mathrm{tp,TWC1}}^{s}\,,$
(34)

i.e. we use the same weighting scheme previously defined in Eq. 2a.

##  III Experimental setup

TABLE III: Engine parameters.  Engine type | VEA Gen I, VEP4 MP
---|---
Number of cylinders | Four, in-line
Displaced volume | 1969 cc
Bore/Stroke | 82 mm/93.2 mm
Compression ratio | 10.8:1
Valve train | DOHC, 16 valves
Intake camshaft | Variable 0-48°CA advance
Exhaust camshaft | Variable 0-30°CA retard
Ignition system | DCI, standard J-gap spark plugs
Fuel system/Injection pressure | DI/200 bar
Fuel | Gasoline RON95 E10
Start of injection | 308-340 CAbTDCf
Boosting system | Turbocharger
Rated power/Rated torque | 187 kW/350 Nm
Stoichiometric air/fuel ratio | 14.01:1

As in our previous study in [1], the experimental setup consisted of a
production Volvo Cars two liter in-line four-cylinder direct injected spark
ignited turbocharged engine rated for 187 kW and 350 Nm, as listed in Table
III. The engine was connected to an electrical dynamometer that regulated the
engine speed and measured the generated torque. A prototyping ECU was used to
sample and change engine parameters. The TWC was close-coupled to the
turbocharger outlet.

The TWC was instrumented with 28 thermocouples (14 in each monolith) and three
exhaust gas sampling locations. The thermocouples, 0.5 $\mathrm{mm}$ type-K
with a grounded hot junction222Manufacturer: RS PRO, model number: 847-1110,
were inserted into a TWC channel and held by friction. A close-up of the
instrumented TWC is shown in Fig. 6, which also shows the direction of gas
flow through the two monoliths. A more detailed drawing of the TWC
construction and the thermocouple locations is shown in Fig. 7.

Several thermocouples failed during the experimental campaign. We believe this
to be due to the combination of fairly sharply bending the thermocouples in
order to reach the required monolith channels and a high level of vibration in
an initial TWC mounting fixture. Fortunately, the most critical sensors (in
slices 1, 3, 4, and 6) were fully functional and only sensors in slices 2 and
5 were damaged. We excluded data from the damaged sensors in our analysis.

Figure 6: The instrumented TWC, with exposed ceramic (left) and metallic
(right) sections. The exhaust gas flow is shown with the highlighted arrow.
The cover to the left is open here for illustrative purposes and tightly
connected to the main body during operation.  Figure 7: Detailed cross-section
of the TWC structure and thermocouple locations. A total of six axial
positions were measured (three per monolith) as indicated. Each axial position
was sampled in one of two configurations, A and B, as indicated. Each specific
sensor is referenced as TSXYZ, where X is the slice number (1-6), Y is the
slice type (A/B), and Z is the sensor position (A/B or A-J).

An auxiliary air feed was added to the exhaust manifold, which allowed for
flushing the entire exhaust subsystem with room-temperature air. By running
the engine in fuel-cut mode (i.e. disabling fuel injection and motoring the
engine with the dynamometer) and injecting auxiliary air into the exhaust
manifold the exhaust aftertreatment system could be cooled to under 100°C in
approximately 5 minutes. The auxiliary airflow was set to 1000
$\mathrm{L}\text{\,}{\mathrm{min}}^{-1}$ STP, which was the maximum flow-rate
supported by the mass flow controller. The auxiliary airflow was completely
disabled during normal operation (i.e. when fuel injection was enabled). A
photograph of the experimental setup is shown in Fig. 8, where the engine is
visible and the TWC is highlighted. A schematic representation of the
experimental set-up and gas flows is shown in Fig. 9, which also highlights
the auxiliary air feed, exhaust gas flow, and gas sampling locations.

The emission sampling points after TWC1 and after TWC2 measured the emissions
exiting the radially central channel of each respective TWC. As both TWCs at
times displayed a large radial temperature differential, this implies that the
average emissions leaving each TWC can be significantly different from the
emissions measured at the radial center.

Figure 8: The experimental setup, with the TWC housing and heat shield
highlighted. The turbocharger is just visible to the right of the TWC. The
exhaust from the TWC is fed down through the visible ducting.  Figure 9:
Schematic representation of the experimental setup, with exhaust gas passing
through the turbocharger, through monolith 1 (M1), monolith 2 (M2), and
finally exiting to the tailpipe. Exhaust gases were sampled at three
locations; directly after the turbocharger, between the monoliths, and after
monolith 2.

###  III-A Data Acquisition

Emissions signals from instruments, fuel consumption, and dynamometer readings
were sampled with a National Instruments DAQ and an associated LabVIEW
program. Engine temperatures, pressures, and the air-fuel ratio was sampled
using acquisition units over a CAN ETAS module. All thermocouples were of type
K. Fuel massflow was measured with a Coriolis meter. All parameters were
sampled at a 10 $\mathrm{Hz}$ rate.

Exhaust gases were sampled from three different locations (as illustrated in
Fig. 9). All sampled gases were extracted with a heated hose (180
$\mathrm{\SIUnitSymbolCelsius}$), followed by a heated conditioning unit (190
$\mathrm{\SIUnitSymbolCelsius}$) with a heated filter and pump. Emissions
concentrations were measured with separate instruments. THC emissions were
measured using a flame ionization detector, NO_x using a chemiluminescence
analyzer, and CO using a non-dispersive infrared detector. The propagation
delay and axial dispersion in hoses and instruments was identified by
recording the measured engine-out emissions during the transition from fuel-
cut operation to normal operation (as will be described in Section III-B).
With this data we compensated for the propagation delay and applied a first-
order high-pass filter to mitigate some of the axial dispersion. This
compensation was applied to the remaining two sampling locations, allowing for
studying transient emission concentration changes moderately well using an
instrument rack primarily intended for steady-state analysis.

As our experimental set-up only allowed for measuring the emissions at one
location at any given time it was crucial for the engine-out emissions to be
consistent between different runs. Due to this we chose to run the combustion
engine in stationary operation, with the goal of maximizing the exhaust gas
composition repeatability. We hypothesize that using hardware that measures
the emission species at every sample point simultaneously would allow for non-
stationary engine operating during cold-start tests.

###  III-B Measurement procedure

The emission measurement equipment was calibrated before measurements using
calibration gases and the engine was heated to its working temperature by
operating it at a moderate load until the coolant reached its working
temperature. The engine was kept warm during the entire test procedure,
implying that the cold-starts studied in this paper refer to the case where
the TWC is initially cold while the engine is at operating temperature.
Furthermore, the TWC was instrumented with heated lambda sensors, and the
engine operated with the conventional closed-loop lambda control scheme during
TWC cold-start tests.

####  III-B1 Steady-state analysis

TABLE IV: Tested steady-state engine operating points. The listed speeds, loads, and spark angles (SA) are setpoint values.  Speed [RPM] | BMEP [bar] |  SA [CAbTDC]
---|---|---
1000 | 2 |  [12, 14, 16, 18, 20, 22, 24]
1000 | 5 |  [6, 8, 10, 12, 14, 16, 18]
1000 | 8 |  [-2, 2, 4]
1500 | 2 |  [12, 14, 16, 18, 20, 21, 22, 24]
1500 | 5 |  [6, 8, 10, 12, 14, 16, 18]
1500 | 8 |  [-2, 2, 4, 6, 8, 10]
2000 | 2 |  [16, 18, 20, 22, 24, 26, 28]
2000 | 5 |  [8, 10, 12, 14, 16, 18, 20]
2000 | 8 |  [2, 4, 6, 8, 10, 12]
2000 | 10 |  [-4, 4]
2000 | 12 |  [-2, 2, 4]
2000 | 14 |  [-2, 2]
2500 | 2 |  [14, 16, 18, 20, 22, 24, 26]
2500 | 5 |  [12, 14, 16, 18, 20, 22, 24]
2500 | 8 |  [8, 10, 12, 14, 16]
2500 | 13 |  [4, 6]
3000 | 8 |  [8, 10, 12, 14, 16, 18]

The goal of this test was to identify the steady-state engine-out emissions
and the associated steady-state radial temperature distribution in the TWC.
This was performed by statically running the engine at a given speed and BMEP
and sweeping the spark angle from the default value and retarding it to the
edge of combustion stability. Table IV lists the tested speeds, BMEPs, and
spark angles tested.

####  III-B2 TWC cold-start characterization

TABLE V: TWC cold-start load points. SA set to ECU default value. Index | Speed [rpm] | BMEP [bar] | SA [CAbTDC]
---|---|---|---
1 | 1000 | 2 | 24
2 | 1000 | 5 | 18
3 | 1500 | 5 | 18
4 | 1500 | 2 | 24
5 | 2000 | 2 | 28
6 | 2000 | 5 | 20
7 | 3000 | 8 | 18
8 | 1000 | 8 | 4
9 | 1500 | 8 | 10
10 | 2000 | 8 | 12

The goal of this test was to characterize the cold-start parameters of the two
TWC’s. The combustion engine was kept at a warm and constant temperature
throughout these tests, i.e. we evaluated the behavior of a cold TWC and warm
engine. We performed this experiment by

  * •

disabling fuel injection (i.e. motoring the engine with the dynamometer) and
opening the auxiliary air valve until all the TWC thermocouples reported a
temperature of under 100 $\mathrm{\SIUnitSymbolCelsius}$,

  * •

first closing the auxiliary air valve, and then immediately enabling ordinary
fuel injection until the TWC reached near-equilibrium temperature and
emissions.

This procedure was repeated for each emission sample point for each of the
load points listed in Table V. These test points were chosen so that some
generated a heating profile that gave a long time to light-off (primarily load
points 1–5), while others reached light-off more quickly (load points 6–10).
The low-load points gave a longer data stream and reduced the relative error
due to axial dispersion in the emission sampling lines. The remaining load
points (points 6–10) were more representative of a conventional heating
strategy, where light-off is reached more quickly. Furthermore, the load
points were characterized by BMEP rather than IMEP due to limitations in the
measurement equipment. The engine was kept at a warm and constant temperature
to ensure that the exhaust gas composition was not influenced by changes in
friction from one load point to the next. It is plausible that regulating for
a given IMEP would give an exhaust gas profile that is less sensitive to
engine temperature.

##  IV Experimental Results

###  IV-A Steady-state

TABLE VI: TWC1 steady-state data for representative load points (speed, BMEP, and spark advance). Figures shown with measured values. Speed | BMEP | SA | $\dot{m}_{\mathrm{exh}}$ | BSFC |  CO |  NO_x |  THC | $TT_{0}$ | $TT_{R/3}$ | $TT_{2R/3}$ | $TT_{R}$
---|---|---|---|---|---|---|---|---|---|---|---
[RPM] | [bar] | [CAbTDC] | [g/s] | [g/kWh] | [ppm] | [ppm] | [ppm] | [$\mathrm{\SIUnitSymbolCelsius}$] | [$\mathrm{\SIUnitSymbolCelsius}$] | [$\mathrm{\SIUnitSymbolCelsius}$] | [$\mathrm{\SIUnitSymbolCelsius}$]
1010 | 6.9 | 10 | 13 | 257 | 3580 | 2500 | 427 | 568 | 567 | 566 | 537
3000 | 7.9 | 18 | 41 | 259 | 7620 | 2550 | 284 | 839 | 839 | 832 | 813
1000 | 1.8 | 12 | 5 | 410 | 6090 | 391 | 548 | 620 | 615 | 609 | 569
1000 | 2.0 | 24 | 5 | 384 | 6650 | 1150 | 684 | 498 | 497 | 491 | 442
1000 | 4.7 | 6 | 9 | 299 | 5950 | 1200 | 459 | 525 | 526 | 519 | 472
1000 | 5.0 | 18 | 9 | 280 | 6550 | 2070 | 533 | 552 | 552 | 545 | 503
1570 | 4.8 | 6 | 16 | 313 | 7720 | 533 | 345 | 726 | 725 | 719 | 684
1570 | 5.1 | 19 | 15 | 279 | 9070 | 1350 | 459 | 668 | 667 | 661 | 627
2000 | 5.0 | 8 | 21 | 316 | 6730 | 279 | 325 | 715 | 716 | 710 | 674
2000 | 5.2 | 21 | 19 | 274 | 9150 | 686 | 442 | 696 | 697 | 693 | 664
3000 | 8.1 | 8 | 46 | 277 | 7060 | 1610 | 195 | 893 | 893 | 884 | 864
3010 | 8.2 | 19 | 43 | 258 | 7560 | 2590 | 328 | 803 | 805 | 797 | 773

The steady-state experimental results were used to generate a table of the
mean equilibrium TWC temperatures, engine-out emissions, exhaust massflow, and
engine BSFC for each of the load points listed in Table IV. Representative
data is shown in Table VI for TWC1. The parameters $TT_{0}$, $TT_{R/3}$,
$TT_{2R/3}$, and $TT_{R}$ correspond to the mean thermocouple temperature for
thermocouples in slice 2 at radius $r=0$, $r=R/3$, $r=2R/3$, $r=R$
respectively (i.e. $TT_{R/3}$ is the mean of TS2BB, TS2BG, TS2BJ).
Importantly, we also also use this experimental data to determine the radial
interpolation profile outlined in Section II-A, i.e. for each load point we
generate an associated interpolated radial temperature profile. More
specifically, for each load point we assume an interpolation function of form

$\displaystyle f_{\mathrm{interp}}$ $\displaystyle=\hat{T}(t^{\prime},r)$
(35a) where $\displaystyle t^{\prime}$
$\displaystyle=\underset{t}{\mathrm{argmin}}\sum_{n=0}^{3}|\hat{T}(t,nR/3)-TT_{nR/3}|\,,$
(35b)

i.e. we let $f_{\mathrm{interp}}$ be the optimal solution in the one-norm
sense that minimizes the deviation between the measured temperatures and the
solution to the heat equation over all time $t^{\prime}$. The one-norm is
consistently used in this paper in an effort to reduce the effect of outliers.
Measured temperatures and the associated interpolation function is shown Fig.
10 for two representative load points.

(a) 1000 RPM, 8 bar, SA 10 CAbTDC

(b) 3000 RPM, 8 bar, SA 18 CAbTDC

Figure 10: Measured normalized radial temperature distribution and least-
squares interpolation for $M=100$ at two representative load points. Figure
11: Representative cold-start temperature evolution. Here shown for load point
1 in Table V (1000 RPM, 2 bar BMEP). Thermocouple locations listed in Fig. 7.
Damaged thermocouples (TS2BE, TS2BH, TS2BI, TS5BC, TS5BE, TS5BH, TS5BJ)
excluded from plots.

###  IV-B Cold-start

Figure 11 shows a representative cold-start temperature evolution, here for
1000 RPM and 2 bar BMEP. We show this specific load point as it gives the
longest system dynamics. We can draw several useful conclusions from this
test;

  * •

The radial temperature distribution is significant throughout both TWC’s, with
a temperature difference between the radial center and periphery of up to 100
$\mathrm{\SIUnitSymbolCelsius}$ near light-off.

  * •

The first TWC shows no major azimuth temperature variation.

  * •

The second TWC does show variations along the azimuth, with the hottest
regions nearer the bottom section of the cross-section (see Fig. 12). We
hypothesize that this is due to increased massflow near the lower sections, as
the sharp bend in the TWC housing causes an uneven pressure distribution
across the inlet to the second TWC.

  * •

TWC2 shows less pronounced axial temperature variations when compared to TWC1.
This could plausibly be due to the length of TWC2, which is only half of TWC1.

These results are consistent for the other load points, which display similar
results.

Based on these results we have chosen to model TWC1 as consisting of three
axial slices, while TWC2 is modeled with a single axial slice. This gives a
total of 3+1 state variables for the TWC1 and 1+1 state variables for TWC2,
i.e. a total of six state variables. Though additional slices would have the
benefit of improving accuracy, the dynamic-programming based optimal control
method we use in Section V is suited for no more than 4–6 state variables. We
have chosen to allocate more slices to the first TWC as it displays the most
significant axial temperature variations.

Figure 12: Projection of TS5Bx thermocouples onto the plane. Approximate gas
flow shown by gray arrows. The thermocouple closest to the bottom (TS5BF)
reaches the highest temperature in Fig. 3, while the thermocouple closest to
the top (TS5BD) reaches the lowest.

###  IV-C Model tuning

The primary goal of the experimental work is to generate measurement data that
is used to tune the TWC model. The tuning process was divided into two
distinct sections where we first tuned the reaction rate parameters, and
afterwards tuned the temperature dynamics parameters. The problem was divided
into two sections to reduce the number of degrees of freedom in each
optimization step. Furthermore, as the temperature dynamics depends on the
exothermic power, it is prudent to first determine the reaction rate
parameters.

Of the 10 cold-start simulations listed in Table V, we designated operating
points $i_{\mathrm{train}}=[1,3,5,7,9]$ as a training set, and
$i_{\mathrm{valid}}=[2,4,6,8,10]$ as a validation set. All model tuning was
done solely using $i_{\mathrm{train}}$, allowing us to later study the model’s
accuracy by studying results of applying the model to the validation set’s
operating points.

###  IV-D Reaction rate parameters

Here, we consider tuning the per-species reaction-rate parameters $A^{s}$ and
$E_{a}^{s}$, giving a total of six parameters to tune per TWC. With the
experimental setup as described in Section III, the gas composition entering
and leaving TWC1’s radially central channel is well-measured. However, the gas
composition entering TWC2 is not as well characterized, as the gas composition
leaving TWC1 is inhomogeneous (due to the large radial temperature gradient)
and partially mixed before entering TWC2. Due to this, we have chosen to first
tune the reaction rate parameters for TWC1, and afterwards make use of the
identical precious metal composition of TWC1 and TWC2 (which differ only in
their washcoat thickness and loading). This allows us to estimate TWC1’s
reaction rate parameters using experimental data and then compute the
equivalent parameters for TWC2.

With respect to TWC1, for a given set of reaction-rate parameters we used the
(measured) temperature evolution of each axial slice to simulate the outgoing
emission concentration. More specifically, we let the measured state evolution
for each operating point be

$\displaystyle T_{\mathrm{meas,TWC1}}$
$\displaystyle=[T_{1},T_{2},T_{3},\Delta_{T}]$ (36a) where $\displaystyle
T_{1}$ $\displaystyle=\mathrm{TS1AA}$ (36b) $\displaystyle T_{2}$
$\displaystyle=\mathrm{TS2BA}$ (36c) $\displaystyle T_{3}$
$\displaystyle=\mathrm{TS3AA}$ (36d) $\displaystyle\Delta_{T}$
$\displaystyle=\mathrm{mean}([\mathrm{TS1AB}-\mathrm{TS1AA},\mathrm{TS2BD}-\mathrm{TS2BA},$
$\displaystyle\mathrm{TS3AB}-\mathrm{TS3AA}])\,,$ (36e)

i.e. the measured radially central temperatures and the mean difference
between the radial center and periphery respectively. This gives a state
vector time-evolution $T(k)_{\mathrm{meas,TWC1}}$, where $k=0,1,2,\dots$
indicates the time-sample of the state vector, sampled at a rate of 10
$\mathrm{Hz}$. In this and later stages we simulated 100 radial channels. We
chose to simulate a relatively large number of channels as the kinetics
submodel was implemented in a semi-parallel manner that did not require
significantly longer to evaluate than for instance 10 channels333This was in
part due to the nature of our model implementation in MATLAB. As the chemical
kinetics submodel consists of only basic arithmetic operations we could
implement the model in a manner that effectively makes use of MATLAB’s
numerically efficient matrix operations. This gave a model implementation
whose execution time was dominated by the the number of function calls, i.e.
the sample rate and simulation time, rather than the number of radial
channels..

We started the tuning process by first determining an initial guess, where we
selected $E_{a}^{s}$ as determined by [18] and $A^{s}$ was selected to give a
light-off temperature near the experimentally measured behavior. We then used
MATLAB’s patternsearch utility (a zero’th order gradient-descent optimization
method) to minimize the 1-norm penalty

$J^{*}=\min_{E_{a,\mathrm{TWC1}}^{s},A_{\mathrm{TWC1}}^{s}}\sum_{i\in
i_{\mathrm{train}}}\sum_{k=1}^{K}|\dot{m}(k)_{3,1,\mathrm{TWC1},i}^{s,\mathrm{out}}\\\
-\dot{m}(k)_{\mathrm{meas,\mathrm{TWC1out},i}}^{s}|$ (37)

using the MADSPositiveBasis2N polling method and with the UseCompletePoll flag
set. Here, $\dot{m}(k)_{3,1,\mathrm{TWC1},i}^{s,\mathrm{out}}$ is the
simulated concentration of emission species $s$ emitted by the third (i.e.
last) axial slice at the radial center at sample $k$ and operating point $i$,
$\dot{m}(k)_{\mathrm{meas,\mathrm{TWC1out},i}}^{s}$ is the measured emissions
at the same time instance and operating point, and $i$ iterates over the
operating points in the training dataset.

With the kinetics parameters for TWC1 determined, we estimate the parameters
for TWC2 by assuming an identical activation energy and with the pre-
exponential term scaled by the amount of catalytically active material. This
gives the estimates

$\displaystyle E_{a,\mathrm{TWC2}}^{s}$
$\displaystyle=E_{a,\mathrm{TWC1}}^{s}$ (38a) $\displaystyle
A_{\mathrm{TWC2}}^{s}$ $\displaystyle=\frac{t_{\mathrm{wash,TWC2}}\cdot
w_{\mathrm{TWC2}}}{t_{\mathrm{wash,TWC1}}\cdot
w_{\mathrm{TWC1}}}A_{\mathrm{TWC1}}^{s}\,,$ (38b)

where $t_{\mathrm{wash}}$ and $w$ are the washcoat thickness and loading
respectively for each TWC.

TABLE VII: Tuned reaction rate parameters. Parameter | Value | Unit
---|---|---
$E_{a,\mathrm{TWC1,2}}^{\ce{CO}}$ | $84.0\cdot 10^{3}$ | $\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}$
$E_{a,\mathrm{TWC1,2}}^{\ce{NO_{x}}}$ | $82.1\cdot 10^{3}$ | $\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}$
$E_{a,\mathrm{TWC1,2}}^{\ce{THC}}$ | $51.0\cdot 10^{3}$ | $\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}$
$A_{\mathrm{TWC1}}^{\ce{CO}}$ | $59.5\cdot 10^{9}$ | -
$A_{\mathrm{TWC1}}^{\ce{NO_{x}}}$ | $27.6\cdot 10^{6}$ | -
$A_{\mathrm{TWC1}}^{\ce{THC}}$ | $11.2\cdot 10^{9}$ | -
$A_{\mathrm{TWC2}}^{\ce{CO}}$ | $23.8\cdot 10^{9}$ | -
$A_{\mathrm{TWC2}}^{\ce{NO_{x}}}$ | $16.1\cdot 10^{6}$ | -
$A_{\mathrm{TWC2}}^{\ce{THC}}$ | $5.63\cdot 10^{9}$ | -

Table VII lists the identified parameters for both TWCs, which are on the same
order of magnitude as literature suggests [18, 30]. Note that the
patternsearch method is similar to gradient-descent methods, and as the
problem is non-convex is therefore not guaranteed to return a globally optimal
solution.

Figure 13: Emissions profiles for 1000 rpm, 2 bar BMEP load point. See Fig. 11
for the associated measured temperature evolution and Fig. 14(a) for the
associated state evolution.

Time-resolved plots of the measured and simulated emissions profiles and the
simulated conversion efficiencies are shown in Fig. 13 for the lowest-load
operating point. The simulated outgoing emissions are shown for the simulated
radially central channel, which corresponds to the location of the measured
emissions. The most significant deviations are seen at $t\in[0,30]$ for NO_x
and at $t\in[0,50]$ for THC. We hypothesize that the former is due to
adsorption and the latter due to the poor transient response of the
measurement equipment.

###  IV-E Temperature dynamics parameters

With the reaction rate parameters determined we consider the temperature
dynamics parameters $k_{\mathrm{ax}}$, $k_{\mathrm{ra}}$, $k_{\mathrm{amb}}$,
$c_{p}$, and $L_{n}$ separately for each TWC. Note that as TWC2 is modeled
with only a single axial slice, there is no modeled axial conduction and
$k_{\mathrm{ax}}$ and $L_{n}$ therefore have no meaning. We tuned each TWC to
the measured data by assigning the initial state to the measured temperature
at the start of the cold-start test and then applying a 1-norm penalty to the
deviation between the simulated and measured states, i.e.

$\displaystyle J^{*}$
$\displaystyle=\min_{k_{\mathrm{ax}},k_{\mathrm{ra}},k_{\mathrm{amb}},c_{p},L_{n}}\sum_{i\in
i_{\mathrm{train}}}\sum_{k=1}^{K}|T(k)_{\mathrm{meas},i}-T(k)_{\mathrm{sim},i}|$
(39)

where $T_{\mathrm{sim}}$ is the simulated state evolution generated by solving
Eq. 21 using an explicit fourth-order Runge-Kutta solver with a fixed time-
step of 0.1 $\mathrm{s}$, $T(0)_{\mathrm{sim}}$ is initialized as
$T(0)_{\mathrm{sim}}=T(0)_{\mathrm{meas}}$, $T_{\mathrm{meas}}$ is defined by
Eq. 36 for TWC1, and for TWC2 defined as

$\displaystyle T_{\mathrm{meas}}$ $\displaystyle=[T_{1},\Delta_{T}]$ (40)
where $\displaystyle T_{1}$
$\displaystyle=\mathrm{mean}([\mathrm{TS4AA},\mathrm{TS5BA},\mathrm{TS6AA}])$
(41) $\displaystyle\Delta_{T}$
$\displaystyle=\mathrm{mean}([\mathrm{TS4AB}-\mathrm{TS4AA},$
$\displaystyle\mathrm{TS5BD}-\mathrm{TS5BA},\mathrm{TS6AB}-\mathrm{TS6AA}])\,.$
(42)

As in the reaction rate parameters, we used the patternsearch method to
determine the optimal parameters. We supplied the initial guess for
$k_{\mathrm{ax}}$, $k_{\mathrm{ra}}$, $k_{\mathrm{amb}}$, and $c_{p}$ by
setting them to the values specified by the TWC manufacturer, and $L_{n}$ to
geometrically ideal values, where we assume the thermocouples are placed in
the center of each slice. Referencing Fig. 7 gives the initial guess
$L_{1}=20\cdot 10^{-3}$ $\mathrm{m}$, $L_{2}=102\cdot 10^{-3}$ $\mathrm{m}$,
and $L_{3}=20\cdot 10^{-3}$ $\mathrm{m}$. Table VIII lists the parameter
values found after tuning, as well as the known (i.e. assigned) fixed model
parameters.

TABLE VIII: Tuned and fixed temperature dynamics parameters. Tuned Parameter | Value | Unit
---|---|---
$k_{\mathrm{ax,TWC1}}$ | 319 | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
$k_{\mathrm{ra,TWC1}}$ | 46.6 | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
$k_{\mathrm{amb,TWC1}}$ | 0.421 | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
$c_{p,\mathrm{TWC1}}$ | 2318 | $\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$
$L_{1,\mathrm{TWC1}}$ | $32.1\cdot 10^{-3}$ | $\mathrm{m}$
$L_{2,\mathrm{TWC1}}$ | $48.2\cdot 10^{-3}$ | $\mathrm{m}$
$L_{3,\mathrm{TWC1}}$ | $41.8\cdot 10^{-3}$ | $\mathrm{m}$
$k_{\mathrm{ra,TWC2}}$ | 4.53 | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
$k_{\mathrm{amb,TWC2}}$ | 0.602 | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
$c_{p,\mathrm{TWC2}}$ | 2360 | $\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$
Fixed Parameter | Value | Unit
$\mathrm{OFA}_{\mathrm{TWC1}}$ | 0.935 | -
$\mathrm{OFA}_{\mathrm{TWC2}}$ | 0.846 | -
$m_{\mathrm{TWC1}}$ | 0.418 | $\mathrm{kg}$
$m_{\mathrm{TWC2}}$ | 0.248 | $\mathrm{kg}$
$c_{p,\mathrm{exh}}$ | 1050 | $\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$
$t_{\mathrm{amb,TWC1}}$ | $10\cdot 10^{-3}$ | $\mathrm{m}$
$t_{\mathrm{amb,TWC2}}$ | $10\cdot 10^{-3}$ | $\mathrm{m}$
$T_{\mathrm{amb}}$ | 25 | $\mathrm{\SIUnitSymbolCelsius}$

An illustration of a representative temperature evolution is shown in Fig.
14(a) and Fig. 14(b). Though the first slice of TWC1 and TWC2 capture the
measured temperature evolution well, the second and third slices of TWC1 do
not capture the characteristic delay shown in the measured data. We
hypothesize that this is independent of the chosen tuning parameters and an
inherent limitation of our modeling assumption of a small number of axial
slices. More specifically, as a general discrete-time delay of $n$ samples
(trivially) requires storing the values of the $n$ samples, our shown model
can only represent a true delay of three samples, i.e. an insignificant 0.3
seconds, before the last axial segment starts displaying a positive
temperature derivative. This can be alleviated somewhat either by increasing
the number of axial segments or by considering cold-starts with a less
prominent delay, as is shown in Fig. 15.

(a) Measured and simulated temperature evolution for TWC1 with 3 axial slices.

(b) Measured and simulated temperature evolution for TWC2.

Figure 14: Temperature evolution for 1000 rpm, 2 bar BMEP operating point.

(a) Measured and simulated temperature evolution at the 1000 rpm, 2 bar BMEP
operating point, TWC1 model modified to use 6 axial slices. Other model
parameters unchanged.

(b) Measured and simulated temperature evolution for TWC1 at the 1500 rpm, 8
bar BMEP load point.

Figure 15: Increasing the number of resolved slices (Fig. 15(a)) and/or
heating the TWC more quickly (Fig. 15(b)) reduces the modeling error caused by
the limited ability to represent a delay.

###  IV-F Cumulative emissions accuracy

With the model tuned, we will now turn to quantitatively evaluating the TWC
model’s accuracy. Here, we consider the relative difference between the
cumulative measured and simulated emissions (i.e. cold-start “bag emissions”)
for each emission species and TWC. Using the notation where $\Delta_{i}^{s}$
corresponds to the $i$’th TWC for emission species $s$ gives

$\displaystyle\Delta_{1}^{s}$
$\displaystyle=\frac{\sum_{k=0}^{K}\dot{m}(k)^{s,\mathrm{out}}_{3,1,\mathrm{TWC1}}}{\sum_{k=0}^{K}\dot{m}(k)^{s}_{\mathrm{meas,TWC1out}}}-1$
(43a) $\displaystyle\Delta_{2}^{s}$
$\displaystyle=\frac{\sum_{k=0}^{K}\dot{m}(k)^{s,\mathrm{out}}_{1,1,\mathrm{TWC2}}}{\sum_{k=0}^{K}\dot{m}(k)^{s}_{\mathrm{meas,TWC2out}}}-1\,.$
(43b)

An illustration of the simulated and measured cold-start emissions is shown in
Fig. 16 for the 1000 rpm, 2 bar BMEP load point, along with the associated
cumulative simulation error. The figure indicates that one significant
contribution to the cumulative error is due to inaccuracies in the measurement
equipment (which is primarily designed for analyzing stationary operation).
This is most clearly seen during the first 20 seconds of operation for THC,
where the measured emissions are significantly larger than than the engine-out
emissions. It is possible that the NO_x emissions are also incorrectly
measured, as the measured emissions are only half of the engine-out emissions
after 3-4 seconds (while light-off occurs after approximately 60 seconds at
this load point). However, it is also plausible that the low NO_x emissions
are correctly measured and this anomaly is instead due to unmodeled adsorption
in the TWC.

A table listing the relative cumulative tailpipe error is shown in Table IX
for each load point. The validation dataset (the lower half) displays an
accuracy comparable to the training dataset (the upper half), indicating that
the model is not overfitted. Furthermore, though there is a significant degree
of variability between the measured and predicted cumulative emissions this is
only somewhat worse than a significantly more complex model [18] which
displays a typical cumulative error on the order of $\pm 20\%$ to $\pm 50\%$.
Furthermore, we hypothesize that the most significant outliers (e.g. the 1500
rpm, 2 bar load point) are to some extent due to process variability and/or
measurement error. Consulting the time evolution for this load point (Fig. 17,
here shown for CO emissions) indicates that this can be a factor, as the
majority of the modeling error arises after 10 seconds when the engine-out
emissions significantly increase but a similar increase is not seen in the
measured emissions.

Figure 16: Measured and simulated emissions for each emission species at a 1000 rpm, 2 bar BMEP load point. TABLE IX: Relative model accuracy for training (upper half) and validation (lower half) load points. Speed | BMEP | $\Delta_{2}^{\ce{CO}}$ | $\Delta_{2}^{\ce{NO_{x}}}$ | $\Delta_{2}^{\ce{THC}}$
---|---|---|---|---
[RPM] | [bar] | [-] | [-] | [-]
997 | 1.99 | -11.6% | -11.6% | -39.6%
1500 | 4.84 | +106.5% | -20.8% | +66.0%
2000 | 2.07 | +44.8% | +34.2% | -7.8%
3000 | 8.17 | +39.5% | +38.9% | -18.2%
1500 | 8.09 | +72.7% | +56.4% | -10.7%
998 | 4.93 | +2.0% | -19.5% | -23.4%
1500 | 2.07 | +38.6% | +41.7% | -10.8%
2010 | 4.93 | +73.3% | +75.8% | +13.1%
999 | 7.96 | +81.2% | +65.8% | -13.8%
2000 | 8.08 | +63.0% | +48.6% | -26.7%
Figure 17: Measured and simulated emissions for CO emissions at the 1500 rpm,
5 bar BMEP load point. Note the large increase in engine emissions after 10
seconds, giving rise to an increase in simulated emissions without an
associated increase in measured emissions.

##  V Optimal Control

We will here illustrate optimal cold-start control as one application of the
presented model. Specifically, we generate an optimal state-feedback
controller suitable for on-line operation that balances the combustion
engine’s fuel efficiency and tailpipe emissions. We model the combustion
engine exhaust using a static mean-value engine model, and allow the
controller to freely choose the engine’s speed, BMEP, and spark angle [31]. We
will then use the controller to generate simulated cold-start temperature and
emission trajectories and compare the results for different weightings of fuel
efficiency and tailpipe emissions.

###  V-A Problem formulation

We introduce the optimal control problem as

$\displaystyle J^{*}$
$\displaystyle=\min_{u}\lim_{K\to\infty}\sum_{k=0}^{K}\mathrm{BSFC}(k)+\Lambda^{T}\cdot\begin{bmatrix}\dot{m}_{\mathrm{tp}}^{\mathrm{\ce{CO}}}(k)\\\
\dot{m}_{\mathrm{tp}}^{\mathrm{\ce{NO_{x}}}}(k)\\\
\dot{m}_{\mathrm{tp}}^{\mathrm{\ce{THC}}}(k)\end{bmatrix}$ (44a) subject to
$\displaystyle x(k+1)$ $\displaystyle=f_{d}(x(k),u(k))$ (44b) $\displaystyle
g(x(k),u(k))$ $\displaystyle\leq 0\,.$ (44c)

Here, $x$ is a state vector corresponding to both TWC’s (i.e.
$[T_{\mathrm{TWC1}};T_{\mathrm{TWC2}}]$), $u$ is a discrete control variable
corresponding to the requested operating point of the combustion engine (i.e.
an integer value that indexes the operating points in Table IV),
$\mathrm{BSFC}(k)$ is the mean BSFC associated with operating point $u(k)$,
$\Lambda$ is a $3\times 1$ tuning parameter that balances the relative weight
given to fuel-efficient operation and minimizing emissions (where smaller
$\Lambda$ prioritizes the BSFC and larger $\Lambda$ prioritizes the level of
emissions), $f_{d}$ is the system dynamics given by solving Eq. 21 for a given
sample time for each TWC, and $g$ is a constraint function that bounds $u$ to
the integer values that index the tested operating points and bounds $x$ to
safe TWC temperatures.

The cost function Eq. 44a is specifically formulated to be of form
$a+\Lambda^{T}\cdot b$, as this is equivalent to minimizing $a$ while limiting
$b\leq B$, i.e. minimizing the average BSFC while limiting the vector of
cumulative emissions to a given level. The same structure is also commonly
seen in Equivalent Consumption Minimization Strategy (ECMS) controllers [32,
33] for the equivalent purpose balancing fuel consumption and electric energy
consumption. We here notationally use $\Lambda$ rather than $\lambda$ to avoid
confusion with conventional notation where $\lambda$ is used to denote the
air-fuel ratio.

Note that Eq. 44a is formulated as an undiscounted infinite-horizon problem,
as this penalizes the BSFC and emissions without requiring time to heat the
TWC to be explicitly specific or known beforehand. Furthermore, by permitting
the engine to operate at any of the points in Table IV we also allow the
engine power to freely vary. The hybrid vehicle cold-start problem is one
example of an application that is well-suited to this cost formulation, as the
electric machine can typically either supply or consume the difference between
the combustion engine power and traction power.

We have solved Eq. 44 using a method developed by the authors [34] based on
approximate dynamic programming and similar to policy iteration methods. The
method, _Undiscounted Control Policy generation by Approximate Dynamic
Programming_ (UCPADP) extends on existing approximate dynamic programming
policy iteration methods by allowing for undiscounted problem formulations,
i.e. infinite-horizon problems where the cost function does not decay with
increasing $k$. In principle, we can solve Eq. 44 without using the above
method by setting $K$ to a sufficiently large value and using a conventional
ADP method [24, 35] to generate a solution. However, it is difficult to
manually determine a sufficiently but not excessively large value $K$.
Conveniently, the UCPADP method also returns a sufficient horizon, which for
the specific TWC and cost formulation studied here was found to be 145
seconds.

One major benefit with UCPADP and other policy iteration methods is that the
optimal control signal can be represented as a control law, i.e. the optimal
control signal can be simply tabulated by the state values. This implies that
a controller can be implemented by simply looking up the optimal control for
the current state. However, this does require knowledge of the current system
state, either by direct measurement or by a state observer that estimates the
system state. Furthermore, note that the selection of the engine’s operating
point is not formulated as a dynamic problem, which in principle implies that
there is no cost associated with rapidly changing the engine’s operating
point. Though our solutions did not exhibit very impractical operating point
changes, we have in our presented results applied a 5-second rolling average
filter to the engine’s target speed, BMEP, and spark angle. Though this makes
for solutions that are not optimal with respect to Eq. 44, we believe that
this results in more suitable engine behavior with somewhat damped transients.

Solving Eq. 44 using an ADP methods first requires the state and control
variables to be discretized. As $u$ is inherently discrete (indexing the
operating points in Table IV) we only need to discretize the states $x$. A
denser discretization will give a solution closer to the true optimal
solution, but at cost of increased memory and computational demand. We have
chosen to discretize the states in TWC1 as

$\displaystyle T_{1}$ $\displaystyle=[0,25,50,75,\dots,900]$ (45a)
$\displaystyle T_{2}$ $\displaystyle=[0,100,200,\dots,900]$ (45b)
$\displaystyle T_{3}$ $\displaystyle=[0,100,200,\dots,900]$ (45c)
$\displaystyle\Delta_{T}$ $\displaystyle=[-200,100]$ (45d) and for TWC2 as
$\displaystyle T_{1}$ $\displaystyle=[0,100,200,\dots,900]$ (45e)
$\displaystyle\Delta_{T}$ $\displaystyle=[-200,100]\,,$ (45f)

i.e. we resolve the first axial slice in TWC1 with fairly high detail, while
the remaining slices and $\Delta_{T}$ is more coarsely resolved.

###  V-B Optimal results

We have solved Eq. 44 for a range of different normalized weights
$\Lambda_{n}$, defined element-wise as

$\displaystyle\Lambda_{n}^{s}$
$\displaystyle=\frac{\Lambda^{s}}{\min\mathrm{BSFC}/\min\dot{m}_{\mathrm{exh}}^{s}}$
(46)

with results listed in Table X. We use the normalized $\Lambda_{n}$ for ease
of reference as $\Lambda_{n}=[1,1,1]$ in some sense equally weighs the fuel
consumption and engine-out emissions. As tailpipe emissions approach zero when
the TWC heats up we can thus view $\Lambda_{n}=1$ as a lower bound of relevant
values to consider.

TABLE X: Performance of optimal controller for varying $\Lambda_{n}$ during a
cold-start (25 $\mathrm{\SIUnitSymbolCelsius}$, no radial distribution) and
half-warm start (200 $\mathrm{\SIUnitSymbolCelsius}$ some radial
distribution). Green-colored cells indicate cases where $\Lambda_{n}$
penalizes only one emission species, while other are ignored (red), and can be
compared with the equally-penalized case (blue).

(a) $T_{1-3,\mathrm{TWC1}}=T_{1,\mathrm{TWC2}}=25$
$\mathrm{\SIUnitSymbolCelsius}$,
$\Delta_{T,\mathrm{TWC1}}=\Delta_{T,\mathrm{TWC2}}=0$
$\mathrm{\SIUnitSymbolCelsius}$.

$\Lambda_{n}^{\ce{CO}}$ | $\Lambda_{n}^{\ce{THC}}$ | $\Lambda_{n}^{\ce{NO_{x}}}$ |  CO |  THC |  NO_x | BSFC | $m_{\mathrm{fuel}}$
---|---|---|---|---|---|---|---
[-] | [-] | [-] | [$\mathrm{mg}$] | [$\mathrm{mg}$] | [$\mathrm{mg}$] | [$\mathrm{g}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$] | [$\mathrm{g}$]
0 | 0 | 0 | 521 | 19.5 | 262 | 250 | 165
1 | 1 | 1 | 521 | 19.5 | 262 | 250 | 165
10 | 10 | 10 | 423 | 18.7 | 153 | 252 | 155
0 | 0 | 10 | 423 | 18.6 | 152 | 252 | 155
0 | 10 | 0 | 521 | 19.5 | 262 | 250 | 165
10 | 0 | 0 | 388 | 21.2 | 239 | 250 | 162
$10^{2}$ | $10^{2}$ | $10^{2}$ | 526 | 35.0 | 41 | 262 | 130
0 | 0 | $10^{2}$ | 687 | 42.6 | 38 | 262 | 126
0 | $10^{2}$ | 0 | 494 | 19.7 | 242 | 250 | 162
$10^{2}$ | 0 | 0 | 340 | 24.6 | 229 | 251 | 157
$10^{3}$ | $10^{3}$ | $10^{3}$ | 529 | 32.7 | 29 | 271 | 117
0 | 0 | $10^{3}$ | 530 | 32.9 | 29 | 270 | 116
0 | $10^{3}$ | 0 | 432 | 11.5 | 151 | 254 | 158
$10^{3}$ | 0 | 0 | 338 | 33.4 | 245 | 253 | 143
$10^{4}$ | $10^{4}$ | $10^{4}$ | 526 | 37.0 | 26 | 298 | 90
0 | 0 | $10^{4}$ | 668 | 85.1 | 15 | 356 | 60
0 | $10^{4}$ | 0 | 433 | 11.4 | 145 | 255 | 155
$10^{4}$ | 0 | 0 | 340 | 33.0 | 232 | 255 | 140

(b) $T_{1-3,\mathrm{TWC1}}=T_{1,\mathrm{TWC2}}=200$
$\mathrm{\SIUnitSymbolCelsius}$,
$\Delta_{T,\mathrm{TWC1}}=\Delta_{T,\mathrm{TWC2}}=-50$
$\mathrm{\SIUnitSymbolCelsius}$.

$\Lambda_{n}^{\ce{CO}}$ | $\Lambda_{n}^{\ce{THC}}$ | $\Lambda_{n}^{\ce{NO_{x}}}$ |  CO |  THC |  NO_x | BSFC | $m_{\mathrm{fuel}}$
---|---|---|---|---|---|---|---
[-] | [-] | [-] | [$\mathrm{mg}$] | [$\mathrm{mg}$] | [$\mathrm{mg}$] | [$\mathrm{g}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$] | [$\mathrm{g}$]
0 | 0 | 0 | 265 | 11.9 | 153 | 250 | 165
1 | 1 | 1 | 265 | 11.9 | 153 | 250 | 165
10 | 10 | 10 | 188 | 14.2 | 100 | 251 | 155
0 | 0 | 10 | 187 | 14.1 | 99 | 251 | 155
0 | 10 | 0 | 266 | 11.9 | 153 | 250 | 165
10 | 0 | 0 | 207 | 14.3 | 129 | 250 | 162
$10^{2}$ | $10^{2}$ | $10^{2}$ | 219 | 19.5 | 20 | 258 | 138
0 | 0 | $10^{2}$ | 271 | 22.8 | 21 | 257 | 137
0 | $10^{2}$ | 0 | 219 | 12.1 | 126 | 250 | 161
$10^{2}$ | 0 | 0 | 167 | 23.0 | 120 | 251 | 149
$10^{3}$ | $10^{3}$ | $10^{3}$ | 212 | 17.4 | 16 | 267 | 125
0 | 0 | $10^{3}$ | 235 | 19.4 | 13 | 264 | 126
0 | $10^{3}$ | 0 | 201 | 9.0 | 80 | 253 | 157
$10^{3}$ | 0 | 0 | 131 | 32.0 | 122 | 252 | 132
$10^{4}$ | $10^{4}$ | $10^{4}$ | 203 | 26.9 | 12 | 303 | 86
0 | 0 | $10^{4}$ | 208 | 35.3 | 6 | 322 | 75
0 | $10^{4}$ | 0 | 212 | 9.9 | 76 | 255 | 152
$10^{4}$ | 0 | 0 | 168 | 36.0 | 81 | 265 | 117

We have simulated the performance of the optimal controller and list the
cumulative emissions, fuel efficiency, and consumed fuel in Table X for
several different $\Lambda_{n}$ and two initial conditions. As expected, with
increasing $\Lambda_{n}$ the sum of penalized emissions decrease, while the
mean BSFC increases. This data also indicates that the potential for reducing
NO_x emissions is significantly larger than CO and THC emissions, as shown in
the last 3 rows where NO_x emissions are reduced by 94% compared to the
unpenalized case, while CO and THC emissions are reduced by 35% and 41%
respectively. Furthermore, there seems to be some degree of conflict with
respect to the individual emissions, as penalizing one species tends to
increase the production of others. This indicates that the solutions shown in
Table X are Pareto optimal, i.e. a given emission species mass cannot be
reduced without either increasing another species’ or the BSFC.

An illustration of the controller’s time-evolution is shown in Fig. 18 for
$\Lambda_{n}=[10^{2},10^{2},10^{2}]$. The cold-start trajectory can be divided
into three sections;

_$t <5$: Initial heating phase._ The engine-out species massflow is kept low
(reducing tailpipe emissions) and the BSFC is not prioritized. At the end of
this phase the first axial slice is hot enough to convert emissions at low
mass-flows.

_$5 <t<25$: Intermediary phase._ With increasing conversion efficiency the
engine-out emissions are gradually allowed to increase, allowing for the BSFC
to be increasingly prioritized. At the end of this phase the first two axial
slices are hot enough to convert emissions at the massflow associated with the
minimum-BSFC operating point.

_$t >25$: Sufficiently-heated phase._ Here the TWC is sufficiently hot for
operation at the minimum-BSFC operating point, which the engine is statically
operated at while the TWC converts virtually all emissions. Note that the
entire TWC is well above light-off after 40 $\mathrm{s}$, i.e. a relatively
short heating interval [3].

These sections can largely be seen for other values of $\Lambda_{n}$, with
shorter times allocated to the initial- and intermediary heating phases as
$\Lambda_{n}$ decreases and longer times with larger $\Lambda_{n}$.

Figure 18: Simulated trajectories for $\Lambda_{n}=[10^{2},10^{2},10^{2}]$. SA
indicates the spark angle, MF the engine-out emissions (solid) and tailpipe
emissions (dashed), and $\eta$ the net conversion efficiency. Figure 19:
Simulated trajectories for $\Lambda_{n}=[10^{2},10^{2},10^{2}]$ with an
initial condition $T_{1-3,\mathrm{TWC1}}=T_{1,\mathrm{TWC2}}=200$
$\mathrm{\SIUnitSymbolCelsius}$,
$\Delta_{T,\mathrm{TWC1}}=\Delta_{T,\mathrm{TWC2}}=-50$
$\mathrm{\SIUnitSymbolCelsius}$.

Though an open-loop control scheme could be easily implemented, i.e. by
“playing back” the speed, BMEP, and spark angle trajectory shown in Fig. 18
without any temperature feedback, this type of controller is potentially
sensitive to system variations. This includes both the initial temperature of
the TWC (where a hotter initial condition will reach light-off more quickly)
as well as variations in the exhaust gas temperature due to the fuel’s
composition, combustion variability, and so on. We found that the optimal
control trajectory for half-warm starts is for some values of $\Lambda_{n}$
nearly identical to a time-shifted version of Fig. 18, as exemplified in Fig.
19, while others displayed significant differences. This indicates the
potential for implementing a quasi-optimal open-loop heating strategy for some
$\Lambda_{n}$.

###  V-C Comparison to suboptimal control

We have compared the optimal controller with a traditional suboptimal heating
strategy. The suboptimal controller was defined such that combustion engine is
run at a constant operating point for $t^{\prime}$ seconds, and then switches
to the minimum-BSFC operating point. The initial operating point was chosen to
be the same as the point chosen by the optimal controller at $t=0$, and
$t^{\prime}$ selected to give the same average BSFC as the case for
$\Lambda_{n}=[10^{2},10^{2},10^{2}]$. A comparison of the optimal and
suboptimal controllers is listed in Table XI, and the time-evolution is shown
in Fig. 20. The CO and THC emissions are virtually identical, but the NO_x
emissions are reduced by 35% in the optimal controller. We can see the source
of this in Fig. 20, where there is significant NO_x slip at $t\in[9,13]$ when
the engine transitions from the heating phase to the minimum-BSFC operation
phase. Though the heating phase could be extended, this would be at cost of
reduced average BSFC.

We can also compare the optimal and suboptimal controllers for the half-warm
start case. It we consider the same optimal and suboptimal controllers, Table
X indicates that the optimal controller attains a mean BSFC of 258
$\mathrm{g}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$, in comparison to the
suboptimal controller’s 262
$\mathrm{g}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$. This difference
corresponds to a 33% reduction relative to the minimum BSFC of 250
$\mathrm{g}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$, indicating the
potential for fuel savings by using a closed-loop cold-start strategy.

Figure 20: Simulated trajectory for the suboptimal controller. Note the significant NO_x slip at $t=10$ $\mathrm{s}$, which is not present in Fig. 18. TABLE XI: Comparison of the optimal (here shown for $\Lambda_{n}=[10^{2}10^{2}10^{2}]$) and sub-optimal controllers and the relative reduction in emissions for the optimal controller.  Cont |  CO |  THC |  NO_x | BSFC | $m_{\mathrm{fuel}}$
---|---|---|---|---|---
[-] | [$\mathrm{mg}$] | [$\mathrm{mg}$] | [$\mathrm{mg}$] | [$\mathrm{g}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$] | [$\mathrm{g}$]
Optimal | 526 | 35.0 | 41.0 | 262 | 130
Suboptimal | 525 | 34.9 | 62.9 | 262 | 143
Difference | -0.2% | -0.3% | 34.8% | - | -

###  V-D Memory footprint

Though manageable in a PC, the memory demand associated with the
discretization in Eq. 45 (with 148,000 permutations) can be problematic in an
ECU that has a wide range of other tasks to perform. However, we can apply a
simple space reduction scheme to significantly reduce the used memory. Rather
than store the full discretization, we can reduce the number of stored
elements by letting successive axial slices only be resolved for temperatures
equal to or below the preceding slices (rounded up to the nearest discretized
value), i.e. storing a table of form similar to that in Table XII.

TABLE XII: Representative table of stored states. Here sorted from first to last column in ascending order. TWC1 | TWC2
---|---
$T_{1}$ | $T_{2}$ | $T_{3}$ | $\Delta_{T}$ | $T_{1}$ | $\Delta_{T}$
0 | 0 | 0 | -200 | 0 | -200
0 | 0 | 0 | -200 | 0 | 100
0 | 0 | 0 | -200 | 100 | -200
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
250 | 300 | 300 | 100 | 200 | -200
250 | 300 | 300 | 100 | 200 | 100
250 | 300 | 300 | 100 | 300 | -200
250 | 300 | 300 | 100 | 300 | 100
275 | 0 | 0 | -200 | 0 | -200
275 | 0 | 0 | -200 | 0 | 100
275 | 0 | 0 | -200 | 100 | -200
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
900 | 900 | 900 | 100 | 900 | 100

Furthermore, we can reduce the resolution of $T_{1}$ for temperatures
significantly below and above light-off by resolving $T_{1}$ with 100
$\mathrm{\SIUnitSymbolCelsius}$ increments for temperatures below 100
$\mathrm{\SIUnitSymbolCelsius}$ and above 400 $\mathrm{\SIUnitSymbolCelsius}$,
i.e.

$\displaystyle T_{1}$
$\displaystyle=[0,100,125,150,\dots,350,375,400,500,\dots,900]\,.$ (47)

Reducing the range of considered values in this manner reduces the number of
stored states to 9500 permutations. Each state permutation is associated with
an optimal engine speed, BMEP, and spark angle. Assuming 4 bits of information
are allocated for each parameter (allowing resolving the speed, BMEP and spark
angle to 16 different and independent values) gives a total storage
requirement of 12 bits per state permutation, for a total non-volatile memory
requirement of $12/8\cdot 9500\approx 13.9$ KiB, which is feasible with
existing ECU hardware. More sophisticated compression schemes have the
potential to further reduce the required memory, for instance by using
decision trees to avoid the need to exhaustively storing every state
permutation in regions where the optimal control is constant.

##  VI Conclusions

In this paper we have extended a physics-based TWC model previously presented
by the authors [1] suited for on-line optimal control. The previously
presented model resolves both axial and radial temperature variations while
limiting the number of state variables, allowing for use with optimal control
methods that construct an optimal control policy (e.g. nonlinear state-
feedback and explicit MPC). In this paper we extended the model to support
varying axial discretization lengths, use tuning parameters expressed in well-
known SI units, model heat generation by the oxidation of hydrogen, consider a
TWC consisting of two separate monoliths of different construction, and use a
more rigorous evaluation method with separate tuning and validation datasets.
Finally, we have used the model to generate a near-optimal controller [34]
that can easily be implemented in existing ECU hardware, requiring no more
than 13.9 KiB (14250 bytes) of nonvolatile memory and at virtually no
computational cost (as the optimal control is given by a simple linear
interpolation operation and a linear rolling average filter). The specific
construction of the cost function allows for systematically trading off fuel
consumption and each individual emission species, giving the ability to tune
the cold-start controller to minimize fuel consumption while individually
limiting the specific level CO, THC, and NO_x emissions.

Our experimental study, though limited by the measurement equipment, shows the
potential for use both for off-line simulation as well as for generating a
near-optimal cold-start controller. Though we experimentally studied the case
of a warm engine and a cold TWC for improved experimental repeatability, we
hypothesize that the controller can be extended to the cold engine case by a
suitable update to the combustion engine exhaust model. Though the measured
predictive accuracy is fairly low (with cumulative cold-start emissions
typically estimated at -20% to +80% of the measured emissions), it is likely
that the experimental setup significantly contributes to this error. The
second monolith is situated after a sharp bend, giving a temperature
distribution that is not particularly well-captured with an axi-radial model.
However, as the majority of the emissions can be converted in the first
monolith for low to moderate load-points this inaccuracy might not be of great
importance. It may therefore be a prudent design decision to solely model and
control the first monolith dynamics in an effort to further reduce the memory
requirements of the controller.

We have simulated the performance of the Pareto-optimal cold-start controller
for several different relative weightings of fuel efficiency (BSFC) and
cumulative emissions for each emissions species, i.e. different points on the
Pareto front. For one representative weighting the optimal controller gives
NO_x emissions that are 35% lower than a traditional cold-start controller
with otherwise identical BSFC and CO and THC emissions. This indicates that an
optimal controller that is generated using the presented model has the
potential to reduce the cold-start emissions, as well as allowing for
systematically adjusting the trade-off between each emission species and fuel
consumption. Furthermore, for some regions on the Pareto-front the optimal
controllers display similar speed, load, and spark angle trajectories for
varying initial TWC temperatures (up to a shift in time). It is therefore
plausible that a close-to optimal controller could be implemented with only a
single temperature sensor by “playing back” a section of the temporally-
resolved optimal control trajectory based on the measured temperature.

Relevant future work includes performing an experimental study using
measurement equipment more suited to transient conditions and that is capable
of measuring emissions at two locations simultaneously. Furthermore, it would
be prudent to experimentally validate the performance of the presented
controller, which in turn requires a method for measuring or estimating the
temperatures in the TWC. Additionally, we hypothesize that characterizing the
engine emissions for different average air-fuel ratios near stoichiometry
could allow for the optimal controller to gain an additional degree of freedom
in balancing the ratio of CO and THC emissions to the NO_x emissions during a
cold-start.

## Funding

This work was performed within the Combustion Engine Research Center at
Chalmers (CERC) with financial support from the Swedish Energy Agency.

## References

  * [1] Jonathan Lock, Kristoffer Clasén, Jonas Sjöblom and Tomas McKelvey  “A Control-Oriented Spatially Resolved Thermal Model of the Three-Way-Catalyst”  In _SAE International_ , 2021  URL: https://www.sae.org/publications/technical-papers/content/2021-01-0597/
  * [2] John B. Heywood  “Internal Combustion Engine Fundamentals”  McGraw Hill, 1988
  * [3] Raffael Hedinger, Philipp Elbert and Christopher Onder  “Optimal Cold-Start Control of a Gasoline Engine”  In _Energies_ 10.10, 2017, pp. 1548  DOI: 10.3390/en10101548
  * [4] Karthik Ramanathan, David H. West and Vemuri Balakotaiah  “Optimal design of catalytic converters for minimizing cold-start emissions”  In _Catalysis Today_ 98.3, 2004, pp. 357–373  DOI: 10.1016/j.cattod.2004.08.003
  * [5] Jianbing Gao et al.  “Review of thermal management of catalytic converters to decrease engine emissions during cold start and warm up”  In _Applied Thermal Engineering_ 147, 2019, pp. 177–187  DOI: 10.1016/j.applthermaleng.2018.10.037
  * [6] Markus Schori, Thomas Boehme, Torsten Jeinsch and Matthias Schultalbers  “Optimal catalytic converter heating in hybrid vehicles”  In _SAE Technical Paper_ , 2014
  * [7] R. Pannag  “Coldstart Modeling and Optimal Control Design for Automotive SI engine”, 2009
  * [8] Giovanni Fiengo, Luigi Glielmo, Stefania Santini and Gabriele Serra  “Control of the exhaust gas emissions during the warm-up process of a TWC-equipped SI engine”  In _IFAC Proceedings Volumes_ 35.1, 2002, pp. 301–306  DOI: 10.3182/20020721-6-ES-1901.01518
  * [9] Farzad Keynejad and Chris Manzie  “Suboptimal Cold Start Strategies for Spark Ignition Engines”  In _IEEE Transactions on Control Systems Technology_ 21.4, 2013, pp. 1295–1308  DOI: 10.1109/TCST.2012.2203821
  * [10] Joachim Braun et al.  “Three-dimensional simulation of the transient behavior of a three-way catalytic converter”, 2002
  * [11] David KS Chen, Edward J Bissett, Se H Oh and David L Van Ostrom  “A three-dimensional model for the analysis of transient thermal and conversion characteristics of monolithic catalytic converters”  In _SAE transactions_ , 1988, pp. 177–189
  * [12] Kyriacos Zygourakis  “Transient operation of monolith catalytic converters: a two-dimensional reactor model and the effects of radially nonuniform flow distributions”  In _Chemical Engineering Science_ 44.9, 1989, pp. 2075–2086
  * [13] GN Pontikakis, GS Konstantas and AM Stamatelos  “Three-way catalytic converter modeling as a modern engineering design tool”  In _J. Eng. Gas Turbines Power_ 126.4, 2004, pp. 906–923
  * [14] Naoki Baba, Katsuyuki Ohsawa and Shigeki Sugiura  “Numerical approach for improving the conversion characteristics of exhaust catalysts under warming-up condition”  In _SAE transactions_ , 1996, pp. 2064–2079
  * [15] Se H Oh and James C Cavendish  “Transients of monolithic catalytic converters. Response to step changes in feedstream temperature as related to controlling automobile emissions”  In _Industrial & Engineering Chemistry Product Research and Development_ 21.1, 1982, pp. 29–37
  * [16] Marcelo Real, Raffael Hedinger, Benjamín Pla and Christopher Onder  “Modelling three-way catalytic converter oriented to engine cold-start conditions”  In _International Journal of Engine Research_ 22.2, 2021, pp. 640–651  DOI: 10.1177/1468087419853145
  * [17] Xieyang Yan et al.  “Modeling Three-Way Catalyst Converters During Cold Starts And Potential Improvements”  In _SAE International_ , 2019, pp. 2019–01–2326  DOI: 10.4271/2019-01-2326
  * [18] Karthik Ramanathan and Chander Shekhar Sharma  “Kinetic Parameters Estimation for Three Way Catalyst Modeling”  In _Industrial & Engineering Chemistry Research_ 50.17, 2011, pp. 9960–9979  DOI: 10.1021/ie200726j
  * [19] E.P. Brandt, Yanying Wang and J.W. Grizzle  “Dynamic modeling of a three-way catalyst for SI engine exhaust emission control”  In _IEEE Transactions on Control Systems Technology_ 8.5, 2000, pp. 767–776  DOI: 10.1109/87.865850
  * [20] Byron T. Shaw, Gerald D. Fischer and J. Hedrick  “A SIMPLIFIED COLDSTART CATALYST THERMAL MODEL TO REDUCE HYDROCARBON EMISSIONS”  In _IFAC Proceedings Volumes_ 35.1, 2002, pp. 307–312  DOI: 10.3182/20020721-6-ES-1901.01519
  * [21] Pierre Michel  “Optimizing fuel consumption and pollutant emissions of gasoline-HEV with catalytic converter”  In _Control Engineering Practice_ , 2017, pp. 8
  * [22] Nasser L. Azad, Pannag R. Sanketi and J. Hedrick  “Determining Model Accuracy Requirements for Automotive Engine Coldstart Hydrocarbon Emissions Control”  In _Journal of Dynamic Systems, Measurement, and Control_ 134.5, 2012, pp. 051002  DOI: 10.1115/1.4006217
  * [23] Zhaoxuan Zhu, Shawn Midlam-Mohler and Marcello Canova  “Development of physics-based three-way catalytic converter model for real-time distributed temperature prediction using proper orthogonal decomposition and collocation”  In _International Journal of Engine Research_ , 2019, pp. 146808741987612  DOI: 10.1177/1468087419876127
  * [24] Dimitri P. Bertsekas  “Dynamic Programming and Optimal Control”  Athena Scientific, 2017
  * [25] Arthur Earl Bryson  “Applied optimal control: optimization, estimation and control”  CRC Press, 1975
  * [26] J.D. Logan, F.W. Gehring and P.R. Halmos  “Applied Partial Differential Equations”, Springer Undergraduate Mathematics Series  Springer, 1998
  * [27] K.. Pattas et al.  “Transient Modeling of 3-Way Catalytic Converters”  In _International Congress & Exposition_, 1994  DOI: 10.4271/940934
  * [28] Steffen Tischer et al.  “Three-Way-Catalyst Modeling - A Comparison of 1D and 2D Simulations”  In _SAE World Congress & Exhibition_, 2007, pp. 2007–01–1071  DOI: 10.4271/2007-01-1071
  * [29] Irman Svraka and Österdahl Wetterhag Linus  “Model Based Catalyst Control”, 2019
  * [30] Sung Bong Kang et al.  “Detailed reaction kinetics for double-layered Pd/Rh bimetallic TWC monolith catalyst”  In _Chemical Engineering Journal_ 241, 2014, pp. 273–287  DOI: 10.1016/j.cej.2013.12.039
  * [31] S.. Chan and J. Zhu  “The Significance of High Value of Ignition Retard Control on the Catalyst Lightoff”  In _SAE International_ , 1996, pp. 962077  DOI: 10.4271/962077
  * [32] Cristian Musardo, Giorgio Rizzoni, Yann Guezennec and Benedetto Staccia  “A-ECMS: An Adaptive Algorithm for Hybrid Electric Vehicle Energy Management”  In _European Journal of Control_ 11.4, 2005, pp. 509–524  DOI: https://doi.org/10.3166/ejc.11.509-524
  * [33] G. Paganelli et al.  “Equivalent consumption minimization strategy for parallel hybrid powertrains”  In _Vehicular Technology Conference. IEEE 55th Vehicular Technology Conference. VTC Spring 2002 (Cat. No.02CH37367)_ 4, 2002, pp. 2076–2081 vol.4  DOI: 10.1109/VTC.2002.1002989
  * [34] Jonathan Lock and Tomas McKelvey  “Undiscounted Control Policy Generation for Continuous-Valued Optimal Control by Approximate Dynamic Programming”, 2021  arXiv:2104.11093 [math.OC]
  * [35] Dimitri P. Bertsekas  “Dynamic Programming and Optimal Control”  Athena Scientific, 2012