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

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 2 O) [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 Paretooptimal 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 (J K −1 kg −1 , W m −1 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. The first monolith in the TWC TWC2 The second monolith in the TWC

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 oxygenstorage 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 lowdimensional 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.
An illustration of the assumed TWC geometry is shown in Figs. 1a and 1b. 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 , . . . , L N (where L n = L) to be different for each slice. We also define the lengths L 1,2 , L 2,3 , . . . , 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. 1b.
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, (a) TWC body. Engine exhaust travels from left to right. Figure derived from [1].
(b) Detail of TWC channels. we model the difference in temperature between the radial center and radial periphery of each slice as as ∆ T . Note that ∆ T is not axially resolved, i.e. we assume that ∆ T is identical for all axial slices, i.e. the periphery temperature of slice 1 is T 1 + ∆ T , slice 2 is T 2 + ∆ 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 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.

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 profilê 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.T (0, r) = 0, r = [0, R]. Note thatT 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.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 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 (whereT (t, r) is scaled to the range [0, 1], and R = 1) are shown in Fig. 2.
In this paper we interpolate the radial TWC temperature profile by using precomputed solutions to the above flatplate problem. More specifically, we assume that the radial temperature profile in slice n is given by A 1T (t , r) + A 2 for a given time instant t . Here, A 1 and A 2 are selected so thatT (t , 0) = T n andT (t , R) = T n + ∆ T (i.e.T is scaled and offset to match the known radial center and periphery temperatures), and t 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 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ṁ 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 aṡ

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 We will assume all reactions are balanced, i.e. for a TWC with 100% conversion efficiency all CO, H 2 , C 3 H 6 , and C 3 H 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)] 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 and m correspond to the number of carbon and hydrogen atoms respectively 1 in each molecule of the fuel. With the RON95 E10 fuel studied in this paper we used the supplierspecified value of m /2n = 0.258. Typically, nitrogen oxides (NO and NO 2 ) and hydrocarbon (C 3 H 6 and C 3 H 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 3 H 6 to C 3 H 8 .
We model the reaction rate k s n,m of an emissions species s in any given cell n, m using an Arrhenius expression of form where R is the ideal gas constant, T n,m is the temperature of cell n, m, E s a is the activation energy of emission species s, and A s is the apparent pre-exponential factor for species s.
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 s n,m . 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 H2 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 H2 a = E CO a and A H2 = A 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 constanttemperature approximation we can explicitly solve Eq. (6) as 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 where V slice,n is the gas volume of slice n and ν 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 where OFA is the open frontal area of the TWC, defined by In this paper we extend [1] to use physically meaningful SI units for all parameters. Hereṁ exh is the exhaust massflow (kg s −1 ), P is the absolute pressure (Pa) in the TWC, which is typically close to the ambient pressure, and R specific is the specific gas constant (J K −1 kg −1 ) for tailpipe ratio of N 2 , CO 2 , O 2 , and H 2 O 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 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 s) [7, p. 64] we hypothesize that moderately-varying dynamic operation with transitions on the order of 0.5 -1 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ṁ s,out n,m =ṁ s,in n,m e −k s n,m tr,n , and, by the conservation of mass, the converted massflow is triviallyṁ For convenience, we also define the conversion efficiency of the entire TWC for a given emission species as 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 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 By Eqs. (15a), (15b), (16a) and (16b) the total temperaturedependent heat of reaction generated in each cell is thus

C. Temperature dynamics
We model the temperature dynamics using a heat balance ODE. Introducing the relative length-weighting matrix we can define the heat-balance ODEs as where P ctr = P ax − P rad + P con,ctr + P exo,ctr (20a) P per = P ax + P rad + P con,per + P exo,per − P amb . (20b) Here, m TWC is the mass of the TWC, c p is its specific heat, P ctr and P per are N ×1 vectors corresponding to the total power developed in the radial center and periphery respectively, and dTctr dt and dTper dt are N × 1 vectors representing the temperature derivative in the radial center and periphery. With these terms, we can construct the total state vector ODE as where mean([x 1 , x 2 , . . . , x n ]) corresponds to the arithmetic mean of the elements in x, i.e. 1/n n i=1 x i . Note that we can interpret N mean(W L ( dTper dt − dTctr dt )) 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.
1) Axial conduction: The axial heat conduction power P 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 where q is the heat flux (W/m 2 ), T 1 and T 2 are two known temperatures (K), and k is the thermal conductivity of the material (W m −1 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 where k 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. 1a. 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 ax (1 − OFA)R 2 π. 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.
2) Radial conduction: The radial heat conduction is modeled in a manner similar to the axial heat conduction. The radial flux is modeled as i.e. a temperature difference of ∆ 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 where 1 is the ones vector of size N × 1.

3) Convection:
The convection heat powers P con,ctr and P 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 ≤ L/5 ∀n.
This gives the convection powers as . . .
whereṁ exh is the exhaust massflow (kg s −1 ), c p,exh is the constant-pressure specific heat of the exhaust gases (J kg −1 K −1 ), and T exh is the temperature of the exhaust gas fed into the TWC (K). 4) Exothermic power: The exothermic power terms P exo,ctr and P 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 Note that the term m−1 M −1 varies from 0 to 1 as m varies from 1 to M .

5) Ambient losses:
The heat losses to the ambient environment are modeled as conductive, with a total flux of Here, k amb is the effective thermal conductivity of the insulating material and t 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 D. Two-monolith structure 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 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 iṡ i.e. we use the same weighting scheme previously defined in Eq. (2a).

III. EXPERIMENTAL SETUP
As in our previous study in [1], the experimental setup consisted of a production Volvo Cars two liter in-line fourcylinder 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 mm type-K with a grounded hot junction 2 , 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.
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 L min −1 STP, which was the maximum flowrate 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.

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 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°C), followed by a heated conditioning unit (190°C) 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 nondispersive 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   Fig. 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. studying transient emission concentration changes moderately well using an instrument rack primarily intended for steadystate 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.

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.
1) Steady-state analysis: The goal of this test was to identify the steady-state engine-out emissions and the associated steadystate 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.
2) TWC cold-start characterization: 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  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 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°C, • 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 lightoff 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.

A. Steady-state
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 T T 0 , T T R/3 , T T 2R/3 , and T T 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. T T 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 i.e. we let f 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 . 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. 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;

B. Cold-start
• The radial temperature distribution is significant throughout both TWC's, with a temperature difference between the radial center and periphery of up to 100°C near lightoff. • 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  TS1AA  TS1AB   TS2BA  TS2BB  TS2BC  TS2BD  TS2BF  TS2BG  TS2BJ   TS3AA  TS3AB   TS4AA  TS4AB   TS5BA  TS5BB  TS5BD  TS5BF  TS5BG  TS5BI TS6AA TS6AB Fig. 11. Representative cold-start temperature evolution. Here shown for load point 1 in Table V (  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.

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 train = [1,3,5,7,9] as a training set, and i valid = [2,4,6,8,10] as a validation set. All model tuning was done solely using i train , allowing us to later study the model's accuracy by studying results of applying the model to the validation set's operating points.

D. Reaction rate parameters
Here, we consider tuning the per-species reaction-rate parameters A s and E s a , 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 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) meas,TWC1 , where k = 0, 1, 2, . . . indicates the time-sample of the state vector, sampled at a rate of 10 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 channels 3 . We started the tuning process by first determining an initial guess, where we selected E s a 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 5.63 · 10 9 using the MADSPositiveBasis2N polling method and with the UseCompletePoll flag set. Here,ṁ(k) s,out 3,1,TWC1,i 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,ṁ(k) s meas,TWC1out,i 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 where t wash and w are the washcoat thickness and loading respectively for each TWC. 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.
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 ∈ [0, 30] for NO x and at t ∈ [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.

E. Temperature dynamics parameters
With the reaction rate parameters determined we consider the temperature dynamics parameters k ax , k ra , k 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 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 Fig. 13. Emissions profiles for 1000 rpm, 2 bar BMEP load point. See Fig. 11 for the associated measured temperature evolution and Fig. 14a for the associated state evolution.

Tuned Parameter
Value Unit then applying a 1-norm penalty to the deviation between the simulated and measured states, i.e.
where T 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 s, T (0) sim is initialized as T (0) sim = T (0) meas , T meas is defined by Eq. (36) for TWC1, and for TWC2 defined as where As in the reaction rate parameters, we used the patternsearch method to determine the optimal parameters. We supplied the initial guess for k ax , k ra , k 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 · 10 −3 m, L 2 = 102 · 10 −3 m, and L 3 = 20 · 10 −3 m. Table VIII lists the parameter values found after tuning, as well as the known (i.e. assigned) fixed model parameters. An illustration of a representative temperature evolution is shown in Fig. 14a and Fig. 14b. 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 (a) Measured and simulated temperature evolution for TWC1 with 3 axial slices.
(b) Measured and simulated temperature evolution for TWC2. 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 discretetime 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.

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 ∆ s i (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. Fig. 15. Increasing the number of resolved slices (Fig. 15a) and/or heating the TWC more quickly (Fig. 15b) reduces the modeling error caused by the limited ability to represent a delay.
corresponds to the i'th TWC for emission species s gives 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 engineout 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 ±20% to ±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.

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 coldstart temperature and emission trajectories and compare the results for different weightings of fuel efficiency and tailpipe emissions.

A. Problem formulation
We introduce the optimal control problem as subject to Here, x is a state vector corresponding to both TWC's (i.e. [T TWC1 ; T 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), BSFC(k) is the mean BSFC associated with operating point u(k), Λ is a 3 × 1 tuning parameter that balances the relative weight given to fuel-efficient operation and minimizing emissions (where smaller Λ prioritizes the BSFC and larger Λ 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 + Λ T · b, as this is equivalent to minimizing a while limiting b ≤ 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 Λ rather than λ to avoid confusion with conventional notation where λ is used to denote the air-fuel ratio.
Note that Eq. (44a) is formulated as an undiscounted infinitehorizon 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 i.e. we resolve the first axial slice in TWC1 with fairly high detail, while the remaining slices and ∆ T is more coarsely resolved.

B. Optimal results
We have solved Eq. (44) for a range of different normalized weights Λ n , defined element-wise as  Table X. We use the normalized Λ n for ease of reference as Λ 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 Λ n = 1 as a lower bound of relevant values to consider.
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 Λ n and two initial conditions. As expected, with increasing Λ 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 Λ 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 s, i.e. a relatively short heating interval [3].
These sections can largely be seen for other values of Λ n , with shorter times allocated to the initial-and intermediary 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 Λ 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 openloop heating strategy for some Λ n .

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 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 selected to give the same average BSFC as the case for Λ 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 ∈ [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 g/kWh, in comparison to the suboptimal controller's 262 g/kWh. This difference corresponds to a 33% reduction relative to the minimum BSFC of 250 g/kWh, indicating the potential for fuel savings by using a closed-loop cold-start strategy.

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. Furthermore, we can reduce the resolution of T 1 for temperatures significantly below and above light-off by resolving T 1 with 100°C increments for temperatures below 100°C and above 400°C, i.e.  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 · 9500 ≈ 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 online 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 wellcaptured 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 coldstart 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.