Version of Record: [Link]

com/science/article/pii/S0952197615002365
Manuscript_bb5c4a041868f0935dc3add29053fe49

An ELM Based Predictive Control Method for HCCI Engines

Vijay Manikandan Janakiramana,∗, XuanLong Nguyenb , Dennis Assanisc

a Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI, USA.
b Department of Statistics, University of Michigan, Ann Arbor, MI, USA.
c Stony Brook University, NY, USA.

Abstract

We formulate and develop a control method for homogeneous charge compression ignition (HCCI) engines using
model predictive control (MPC) and models learned from operational data. An HCCI engine is a highly efficient but
complex combustion system that operates with a high fuel efficiency and reduced emissions compared to the present
technology. HCCI control is a nonlinear, multi-input multi-output problem with state and actuator constraints which
makes controller design a challenging task. In this paper, we propose an MPC approach where the constraints are
elegantly included in the control problem along with optimality in control. We develop the engine models using
experimental data so that the complexity and time involved in the modeling process can be reduced. An extreme
learning machine (ELM) is used to capture the engine dynamic behavior and is used by the MPC controller to evaluate
control actions. We also used a simplified quadratic programming making use of the convexity of the MPC problem
so that the algorithm can be implemented on the engine control unit that is limited in memory. The working and
effectiveness of the proposed MPC methodology has been analyzed in simulation using a nonlinear HCCI engine
model. The controller tracks several reference signals taking into account the constraints defined by HCCI states,
actuators and operational limits.
Keywords: Model Predictive Control, Extreme Learning Machines, Engine Control, HCCI Engine, System
Identification, Neural Networks.

1. Introduction

In recent years, the requirements on automotive performance, emissions and cost have become increasingly strin-
gent. As a consequence, the auto industry has been continuously introducing advanced technology to meet these
demands (EPA, 2012). Invariably, such systems introduce additional complexity and associated challenges in design
and operation. Homogeneous charge compression ignition (HCCI) is an advanced combustion technology integrat-
ing several of the recent automotive advancements including variable valve timing, exhaust gas recirculation (EGR),
intake charge boosting etc (EPA, 2012). Although HCCI engines are very efficient combustion systems, they also

∗ Corresponding author. Postal Address: UARC @ NASA Ames Research Center, Moffett Field, CA, USA. Phone: 001-734-358-6633. Email:

vijai@[Link].

Preprint submitted to Elsevier October 17, 2015

© 2015 published by Elsevier. This manuscript is made available under the Elsevier user license
[Link]
pose several challenges when it comes to practically using it for automobiles. Some of them include the absence of a
direct control over combustion, having a narrow operating region and high sensitivity to ambient conditions. Control
of HCCI combustion is a challenging problem and a model based approach is typically employed (Bengtsson et al.,
2006a; Ravi et al., 2009), where a control-oriented model developed using first principles is used. Although physics
based models are the only choice before and during the engine design stage for making design decisions, such models
usually require expert knowledge and take a long time to develop. However when a matured prototype is available,
controller development needs to be accelerated. To this end, we propose to use data-based predictive models that can
be developed fairly quickly without much expert knowledge.
The HCCI engine is equipped with sensors such as manifold temperature and pressure sensors, in-cylinder pressure
transducers that can provide valuable data about the engine’s behavior. We can draw insights about the general
combustion behavior of the engine including efficiency, emissions and combustion stability using the in-cylinder
pressure measurements. We have shown in our prior work that in-cylinder pressure data can be used to develop
control oriented models of the engine’s state variables and operating envelope using neural networks (Janakiraman
et al., 2013) and support vector machines (Janakiraman et al., 2014, 2015). In this paper, we show how to use these
models to develop an optimal controller for the HCCI engine.
For the HCCI modeling problem, Extreme Learning Machines (ELM) (Huang et al., 2006) were selected for
their fast operation and their extraordinary approximation capabilities to fit nonlinear systems. ELM techniques have
been applied successfully to engine fault detection (Yin et al., 2014; Wong et al., 2015b) and modeling diesel and
biodiesel engines using data (Wong et al., 2013, 2015a; Mozaffari and Azad, 2014) showing its suitability in modeling
combustion engines. Furthermore, ELM models have been successfully applied for automotive control purposes
where control algorithms using ELM have been developed for engine cold start emission control (Azad et al., 2015),
air-ratio control (Wong et al., 2014) and autonomous cruise control of electric vehicles (Mozaffari et al., 2015). The
major advantage of ELM (and data based modeling in general) is that it approximates the real system and makes
no (or minimal) simplifying assumptions about the underlying phenomena. The dynamics of sensors, actuators and
other complex processes which are usually overlooked/hard to model using first principles, can be captured using
the identification method. In addition, for a system like the combustion engine, prototype hardware is typically
available and extensive experimental data can be collected making the data based approach more suitable. Data
based identification is less common for HCCI engines owing to its nonlinear, highly sensitive and unstable behavior.
However, we have shown in our prior work that artificial neural networks (Janakiraman et al., 2013) and support vector
machines (Janakiraman et al., 2014) are indeed suited for accurately predicting the dynamics of HCCI. Artificial neural
networks (based on backpropagation) have issues with local minima while support vector machines result in memory
intensive models making them unsuitable for real engine implementations (Janakiraman et al., 2014). ELMs on the
other hand, solve a convex problem resulting in a global optimal solution and are simple enough to be implemented
onboard the engine control unit (ECU), which makes it a good candidate for the problem considered in this work.
As our modeling approach is a black-box type, we propose to use a model based predictive control (MPC) approach
2
for the HCCI control problem because MPC can efficiently handle black-box type models for evaluating control
candidates (Maciej, 2009). In addition, MPC is well suited for handling operation related and hardware related
constraints in an elegant manner. Even for simple linear systems, MPC has been shown to outperform traditional
control such as PID (Maciejowski, 2002) and LQ control (Franko et al., 2011; Bengtsson et al., 2006b). In general
terms, MPC makes use of a mathematical model of the system and solves an optimization problem with the given
constraints to achieve an optimal control solution (Maciejowski, 2002). In this paper, an ELM engine model is used
to make predictions which are then used by the MPC module to make control decisions for tracking a given reference
command. Model predictive control has been quite popular in the automotive domain for both spark ignited (Re et al.,
2010) and compression ignited engines (Ortner et al., 2006; Ortner and Del Re, 2007). However, the application of
MPC to HCCI engines is at an early stage involving physics based models (Widd et al., 2013; Li and Wang, 2013) and
simple identification models (Bengtsson et al., 2006a). This article advances the MPC application to HCCI engines
by demonstrating it on a data based non-linear model that could capture more complex behavior of the engine for a
wide range of operations. Further, by using concepts from convex optimization, the MPC is converted into a quadratic
programming problem and solved using minimum overhead suitable for implementing on-board the engine ECU.
Thus the main contribution of this article is an integrated control framework that applies existing ideas from ELM,
predictive control and convex optimization to a complex combustion engine. This includes optimal model selection
(tuning model hyper-parameters), training and validation of the model to the HCCI engine data, linearizing the ELM
model considering the structure of ELMs, and a simplified quadratic program update developed for the HCCI system
suitable for real-time application.
As there are several components involved, a road-map of our methodology is shown in Figure 1. Real experimental
data from a prototype HCCI engine is taken offline and used to train ELM models. A supervised batch learning
algorithm is used to train the ELM models and the model is validated for generalization performance using real engine
test data. The ELM models act as a proxy to the actual engine where candidate control trajectories are evaluated by
the MPC controller. The evaluation is done by solving an optimization problem that has a well defined cost function
with constraints. The MPC controller is required to provide optimal control solution at every sampling interval of the
engine. We make use of model linearization and the convexity of the MPC problem to develop a simple quadratic
programming subproblem that can be solved fairly quickly, well within the sampling periods of the engine. For
every sampling interval and for a given reference command, an optimal control trajectory is designed by our control
algorithm and is fed to the engine. After application of the first control action, the engine outputs are sampled and
are fed to the MPC algorithm to calculate the optimal control for the next time step. Each of the above are detailed in
the rest of the paper which is organized as follows. First, we describe the HCCI engine experimental setup, design of
experiments and data collection in Section 2. A system identification procedure using ELM is described in Section 3
where predictive models are developed and validated using real HCCI experimental data. A model predictive control
algorithm using linearized ELM models is derived in Section 4 where linearization is done using the structure of ELM
to avoid numerical gradient calculations. Finally, the ELM based MPC approach is demonstrated in simulations in
3
 ELM Modeling
ELM learning
HCCI Input Data

Model prediction

HCCI Output Data

ELM Model
Offline
learning

Offline trained
nonlinear model
MPC Controller
Linearized
ELM
Linearization

Reference
Cost
command Fast
function
Quadratic
Program Constraints

HCCI
Engine

Figure 1: An integrated control framework for HCCI combustion control showing the use of ELM based engine models and model predictive
controller using a fast quadratic programming solver.

Section 5. The simulation case studies help understand the working and effectiveness of the proposed method.

2. HCCI Engine System And Experimentation

In this section, the HCCI engine basics, experimental setup and data collection strategy are described1 . The
engine is a 4 stroke, inline 4 cylinder gasoline engine that has been modified to enable HCCI operation, including an
increase in compression ratio, having a direct fuel injection, intake boosting and variable valve timing. A schematic
of the experimental setup and instrumentation is shown in Fig. 22 . The sensors of the engine include a fast in-cylinder
pressure transducer, thermocouples at the intake and exhaust manifold. The engine is operated with natural aspiration;
i.e., no or negligible difference between intake and exhaust manifold pressures.
HCCI operation is achieved using a recompression strategy. The engine is a discrete event system where one com-
bustion cycle (720 deg crank angle) represents an event. A snapshot of the cylinder pressure trace of one combustion
cycle along with valve events and fuel injection events are shown in Fig. 33 . During every combustion cycle, three
distinct control actions can be defined - the crank angle at intake valve opening (IVO), crank angle at exhaust valve

1 Some content in this section has been adopted from (Janakiraman et al., 2013, 2015) and included in this article for completeness.
2 This Figure appears in Janakiraman et al. (2015). We hope to get permissions from the journal upon acceptance of this article.
3 This Figure appears in Janakiraman et al. (2013). We hope to get permissions from the journal upon acceptance of this article.

4
 Air in

9 7 3
VVT 15
8
Turbocharger
6
13 11

Throttle Intake Manifold Exhaust Manifold Exhaust out

2 5
1
12 10 4

Coolant 0𝐶 8. Fuel injector 𝑚𝑚𝑚𝑚

Piston 1.
thermocouple
Exhaust manifold 9. Cam phaser −
0
2. 𝐶
thermocouple
Crank shaft 10, Cylinder 1 Runner 0𝐶 ,
𝑚
3.
Exhaust manifold
𝑚𝑚𝑚𝑚 11. T and P 𝑚𝑚𝑚𝑚
pressure transducer
14 Intake manifold 0
Post-turbine 12. 𝐶
4. 0
𝐶 thermocouple
thermocouple
Intake manifold
13. 𝑚𝑚𝑚𝑚
5. Lambda sensor − pressure transducer

Air mass flow 14. RPM sensor 𝑅𝑅𝑅

6. 𝑘𝑘/ℎ
sensor
In-cylinder pressure
15. 𝑏𝑏𝑏
7. Spark plug − transducer

Figure 2: A schematic of the HCCI engine showing relevant sensors and instrumentation.

closing (EVC) and crank angle at start of fuel injection (SOI). The valve events are measured in degrees after exhaust
top dead center (deg eTDC) while SOI is measured in degrees after combustion top dead center (deg cTDC). It should
be noted that the fuel mass (FM in mg/cyc) is injected in a region between EVC and IVO defined as the negative valve
overlap (NVO). The fuel is injected early (during NVO of the previous cycle) and given sufficient time to mix with
air forming a homogeneous mixture. A large fraction of the exhaust gas from the previous combustion cycle (EGR) is
retained to elevate the temperature and hence the reaction rates of the fuel and air mixture. The variable valve timing
capability of the engine enables trapping suitable quantities of exhaust gas in the cylinder. The fuel-air mixture is
given sufficient time to mix homogeneously and auto-ignited during the compression stroke. Since the fuel is mixed
homogeneously throughout the cylinder, combustion of the fuel-air mixture happens almost instantaneously, releasing
heat and pushing the piston down performing work. More details on setup and operation of the HCCI engine can be
found in (Janakiraman et al., 2013, 2015).

5
 cTDC eTDC cTDC eTDC
35
Cycle Cycle Cycle
30 (k−1) (k) (k+1)
25

In−cylinder Pressure
20 SOI

(bar)
15 Data
Record
10
EVC IVO
5

−5
−180 0 180 360 540 720 900 1080
Crank Angle (deg)

Figure 3: The pressure trace of a HCCI combustion cycle showing cycle definition, actuator ranges of intake valve opening (IVO), exhaust valve
closing (EVC), start of injection (SOI). The instant when the data is recorded in each cycle is also shown. Each “hump” in pressure corresponds
to the piston moving up and compressing the gas. Each combustion cycle has a large compression “hump” and a smaller re-compression “hump”
during which the fuel is injected.

2.1. HCCI Inputs and Outputs

Before we dive into the modeling and control aspect of our work, let us introduce some inputs (control variables)
and outputs (response variables) of the engine. HCCI combustion is dominated by chemical kinetics and is limited by
stability constraints such as misfire and knock (Kalghatgi and Head, 2006; Shahbakhti and Koch, 2008). To achieve
a stable combustion, a right proportion of fuel, air and EGR is required for a given thermal condition of the engine.
Further, the combustion must occur at the right instant during the expansion stroke (combustion phasing) for high
efficiency and low emissions. By varying the valve events, the quantity of EGR can be varied while FM and SOI
can control the ignition delay that affects combustion phasing. Thus, the engine operation can be controlled using
quantities such as FM, IVO, EVC and SOI. Other important physical variables that influence the performance of
HCCI combustion include intake manifold temperature T in , intake manifold pressure Pin , mass flow rate of air at
intake ṁin , exhaust gas temperature T ex , exhaust manifold pressure Pex , coolant temperature T c , fuel to air ratio (FA)
etc.
Controlling the engine translates to achieving some optimal performance measured using the following metrics.
The indicated mean effective pressure (IMEP) represents the work output from the engine and is one of the variables to
be controlled. For instance, the driver’s power demand can be interpreted as a change in desired IMEP command given
to the engine controller. The second and most important performance metric is the combustion phasing indicated by
the crank angle at 50% mass fraction burned (CA50). The CA50 is a key quantity that influences engine performance,
efficiency, emissions and stability (Bengtsson et al., 2004; Olsson et al., 2001). Further, additional quantities that give
an indication of the combustion stability such as maximum pressure (Pmax ) and maximum pressure rise rate (Rmax ) in a
combustion cycle are considered performance metrics to be monitored. The equivalent air-fuel ratio (EAFR) indicates
if the mixture in the cylinder is fuel rich or fuel lean is another important parameter to be monitored.
For the HCCI control discussed in this article, the control knobs such as FM, EVC and SOI are considered inputs
while performance variables such as IMEP, CA50, Pmax , Rmax , engine torque and EAFR are considered outputs. The

6
 Fuel Mass (mg/cyc)
11 4

IMEP (bar)
10
3
9

8 2
0 500 1000 1500 2000 0 500 1000 1500 2000

CA50 (deg bCTDC)

EVC (deg aCTDC)
−90 0

−100
−10
−110

−120 −20
0 500 1000 1500 2000 0 500 1000 1500 2000

Rise Rate (bar/deg)

400 10
SOI (deg aETDC)

Max Pressure
350
5
300

250 0
0 500 1000 1500 2000 0 500 1000 1500 2000
Max pressure (bar)

50 1.4

EAFR (−)
40
1.2
30

20 1
0 500 1000 1500 2000 0 500 1000 1500 2000
Engine Cycles Engine Cycles

Figure 4: A subset of closed loop experimental data showing the A-PRBS inputs and the measured engine outputs.

IMEP, CA50, Pmax and Rmax are calculated from the high speed in-cylinder pressure measurements.

2.2. Experiment Design

One of our objectives of this work is to model the HCCI engine using data. The HCCI engine is a dynamic
system where the data is correlated as a sequence. Thus, dynamic data that represents the engine behavior needs to
be recorded and used for model development. This is referred to as system identification Nelles (2001). In order to
obtain sufficiently rich data, the engine needs to be excited by an input signal that has all frequencies and amplitudes
(Nowak and Van Veen, 1993) so that all modes of engine operation can be observed in the data. Such a signal can
be called a persistently exciting signal. Persistent excitation cannot be guaranteed for nonlinear systems (Nowak and
Van Veen, 1993) but in practice, an amplitude modulated pseudo-random binary sequence (A-PRBS) works well for
this problem (Nelles, 2001). A set of transient experiments is conducted at a constant speed of 1800 RPM (rotations
per minute) and naturally aspirated conditions using a feedback controller developed for HCCI using first principles
(Jade et al., 2012). An A-PRBS sequence of reference IMEP and CA50 are given as inputs to the feedback controller
while the controller calculates the commands of FM, EVC and SOI to be given to the engine. The experiments are
conducted and data recorded using specialized engine rapid prototyping hardware. The data is sampled using the AVL
Indiset acquisition system where in-cylinder pressure is sensed every crank angle while IMEP, CA50, Pmax and Rmax
are determined on a per-combustion cycle basis. A subset of the HCCI input and output sensor data can be shown
using Fig. 4. The data is processed, scaled and converted to the required format for ELM training.

7
3. HCCI Engine Model Development Using Extreme Learning Machines

3.1. Learning HCCI Dynamics using Recurrent Models

For the HCCI engine system, both the inputs and the outputs of the engine are available as sensor measurements
and naturally, supervised learning methods can be used. The HCCI engine is a nonlinear dynamic system and sensor
measurements represent discrete time sequences. We consider a nonlinear auto regressive model with exogenous input
(NARX) (Nelles, 2001) as follows

y(k) = fNARX [u(k − 1), .., u(k − nu ), y(k − 1), .., y(k − ny )], (1)

where u(k) ∈ Rud and y(k) ∈ Ryd represent the inputs and outputs of the system respectively, k represents the discrete
time index, fNARX (.) represents the nonlinear function mapping specified by the model, nu , ny represent the number of
past input and output samples required (order of the system) while ud and yd represent the dimension of inputs and
outputs respectively. Let
x(k) = [u(k − 1), .., u(k − nu ), y(k − 1), .., y(k − ny )]T (2)

represent the augmented input vector obtained by appending the input and output measurements from the system. The
engine measurement sequence can be converted to the form of training data

{(x1 , y1 ), ..., (xN , yN )} ∈ X, Y

(3)

where N denotes the number of training samples, X denotes the space of the input features (Here X = Rud nu +yd ny and
Y = Ryd ). The above conversion of system measurements to training data can be used to train an ELM model in a
series-parallel architecture, where the past inputs and outputs of the system (concatenated in x) is used for one-step
ahead predictions (OSAP); i.e., given a set of measurements until time index k, the model predicts the output at time
k + 1 (see equation (4) and Fig. 5a). A parallel architecture on the other hand can be used to perform multiple step
ahead predictions (MSAP) by feeding back the predictions of the OSAP model in a recursive manner (see equation
(5) and Fig. 5b). For further reading on series-parallel and parallel architectures, the reader is referred to (Narendra
and Parthasarathy, 1990).

ŷ(k + 1) = fˆNARX [u(k), .., u(k − nu + 1), y(k), .., y(k − ny + 1)] (4)

ŷ(k + n pred ) = fˆNARX [u(k + n pred − 1), .., u(k − nu + n pred ), ŷ(k + n pred − 1), .., ŷ(k − ny + n pred )]. (5)

The OSAP model is used for training because we can use existing simple training algorithms and once we train an
accurate OSAP model, we can convert it to an MSAP model recursively in a straightforward manner by feeding back
the model predictions to predict successively. The MSAP model can be used for making long term predictions useful
for predictive control discussed in Section 5.

8
 (𝑥0 , 𝑦0 ) (𝑥1 , 𝑦1 ) (𝑥2 , 𝑦2 ) (𝑥𝑁𝑝𝑝𝑝𝑝 , 𝑦𝑁𝑝𝑝𝑝𝑝 )

Model Model Model Model

𝑦0 𝑦�1 𝑦�2 𝑦�𝑁𝑝𝑝𝑝𝑝

(a) Series-Parallel architecture.

𝑥0 𝑥1 𝑥2 𝑥𝑁𝑝𝑝𝑝𝑝

Model Model Model Model

𝑦0 𝑦�1 𝑦�2 𝑦�𝑁𝑝𝑝𝑝𝑝

(b) Parallel Architecture.

Figure 5: Illustration showing the model architectures for HCCI engine identification.

3.2. Extreme Learning Machines

Extreme Learning Machines is an emerging learning paradigm for multi-class classification and regression prob-
lems (Huang et al., 2006, 2012) and has outperformed some state of the art algorithms such as backpropagation neural
nets, support vector machines etc. The highlight of ELM is that the training speed is extremely fast with superior gen-
eralization capabilities. The key enabler for ELM’s training speed is the random assignment of input layer parameters
which do not require adaptation to the data. In such a setup, the output layer parameters can be determined analytically
using linear least squares (See Fig. 6). Some of the attractive features of ELM include the universal approximation
capability, better generalization, the convex optimization problem of ELM resulting in the smallest training error with-
out getting trapped in local minima, availability of a closed form solution eliminating iterative training (Huang et al.,
2006).
ELM training involves solving the following optimization problem
n o
min kHW − Yk2 + λkWk2 (6)
W

φ = H T = ψ(WrT x(k) + br ) ∈ Rnh ×1 (7)

where λ represents the regularization coefficient, Y represents the vector of outputs, ψ represents the hidden layer
activation function (typically sigmoidal or radial basis functions) and Wr , W represents the input and output layer
parameters respectively. Here, nh represents the number of hidden neurons of the ELM model and H represents the
hidden layer output matrix.
9
 Input Neurons

Hidden Neurons

Output Neurons

x D
ŷ
D
ϕ

D
W
Linear
Wᵣ Regression

Random
Projection

Figure 6: Illustration showing an extreme learning machine model.

A prominent feature of ELM is that the nonlinear optimization problem is reduced to a linear parameter estimation
problem. This reduction is made possible by the random assignment of the input layer parameters. The matrix Wr
consists of randomly assigned elements that maps the input vector to a high dimensional feature space while br ∈ Rnh is
a bias component assigned in a random manner similar to Wr . The number of hidden neurons determine the dimension
of the transformed feature space. The random parameters can be assigned based on any continuous random distribution
(Huang et al., 2012) and remains fixed during the training process. The intuition behind the random parametrization
of ELM can be explained as follows. By assigning random weights as described above, along with an activation
function, a regressor φ with several functional combinations (with different order terms) of the input feature vector x
is obtained. For a complex problem, more hidden neurons are required which can be thought of as introducing more
complex feature mappings. When the dimension of the regressor φ is sufficiently high with sufficient inclusion of
higher order (nonlinear) features, it may be mapped to the outputs linearly. Hence the training reduces to a single step
calculation given by equation (8). The ELM decision hypothesis can be expressed as in equation (9)
 −1
W ∗ = H T H + λI H T Y. (8)

f (x) = W T [ψ(WrT x + br )]. (9)

3.3. Model Development and Validation

The engine experimental data is split into training and testing sets. The training set consists of about 16000
cycles of data while the testing set consists of about 7000 cycles. In order to validate multi-step ahead prediction
performance, a separate MSAP-testing set containing about 2400 cycles of data is used. Firstly, the model hyper-
parameters are tuned using a cross-validation approach (Janakiraman et al., 2013, 2014), where a small subset of the
training set is used to train the model with several combinations of hyper-parameters4 and testing its performance

4 We performed a grid-search with hidden neurons taking values of [5,10,20,40,100,1000] and regularization parameter taking values [1E-4,1E-
2,1E-1,0,1,10,100].

10
on the rest of the unseen training set. After tuning the model hyper-parameters using cross-validation, an optimal
model structure is determined. Our optimal HCCI engine model consists of 20 hidden neurons with a regularization
coefficient (λ) of 0.001 and a system order (n0 ) identified (automatically from data) to be 1. Although we trained
with about 16000 data samples, our cross-validation results indicate that a simple model of 20 hidden neurons were
able to capture the engine dynamics. It should be noted that the data is recorded at a constant speed and a narrow
engine operation which may be a reason behind a simpler ELM model. The model is retrained with the optimal
hyper-parameters and evaluated on an unseen test data set. The testing and MSAP-testing data sets are never seen by
the models during training.
The models developed in this section are intended for use in a predictive control framework where real time
predictions are required in the engine ECU. In order to evaluate the prediction capability of the models, a root mean
squared error (RMSE) given by equation (10) is used. The models are trained using a series-parallel architecture
(Narendra and Parthasarathy, 1990) and one-step ahead prediction performance is evaluated on the unseen testing
data set. Finally, the model predictions are evaluated for n pred steps ahead in time using a separate input sequence.
This gives a good measure of the long term predictions as well as the generalization capability of the models to unseen
situations. v
u
t N yd
1 XX
RMS E = (yi − ŷij )2 . (10)
N i=1 j=1 j

The test RMSE for one-step ahead prediction is observed to be 0.1085 while for a 600-step ahead prediction, the test
RMSE is observed to be 0.1092. The predictions are summarized in Fig. 7 for an unseen MSAP-testing data set. The
RMSE values give a good indication of how accurately the model can predict the test data. The MSAP predictionsand
the MSAP RMSE values prove that the model can be used for predicting multiple steps ahead in time and can be used
as a predictive model in the MPC framework that is discussed in the next section.

4. Predictive Control formulation using Extreme Learning Machines

In this section, the MPC problem is formulated based on general predictive control principles (Maciejowski, 2002;
Maciej, 2009) for a class of nonlinear systems using ELM models. The formulation will be applied to the HCCI engine
problem in the later sections. The idea is to linearize the nonlinear ELM models at every instant of the engine operation
and use the linearized model in the MPC framework to evaluate control candidates. From the control candidates, the
MPC picks the one that optimizes a certain objective function. In this section, we first formulate the MPC problem
more generally and then apply to the HCCI engine data in Section 5.
Consider a general class of nonlinear discrete time system,

z(k + 1) = f (z(k), u(k))

y(k) = g(z(k)), (11)

11
 3.4

Pmax (bar) CA50 (deg bCTDC) IMEP (bar)

3.2
3
2.8
2.6
0 100 200 300 400 500 600
−2
−4
−6
−8
−10
−12
−14
0 100 200 300 400 500 600
40

30
(bar/deg CA)

0 100 200 300 400 500 600

6
4
2
0
max

0 100 200 300 400 500 600

Torque (Nm) R

30

20
0 100 200 300 400 500 600
1.4
EAFR

1.2

0 100 200 300 400 500 600

Combustion Cycles

Figure 7: Multi-step prediction performance of the ELM HCCI engine model. The model is predicting the shown 7 engine variables for 600 cycles,
given the values at the first cycle and a 600 cycle history of the inputs such as FM, EVC and SOI. The prediction is shown in red while the actual
data is shown in blue. The inputs are in the form of A-PRBS signals and is not shown because it is common for both prediction and actual data.

where z ∈ Rn , u ∈ Rm , y ∈ R p , k the time index, n, m and p represent the number of states, inputs and outputs
respectively. f (.) and g(.) are continuously differentiable and globally Lipschitz (Ross, 1980) nonlinear functions5 .
The models can be linearized around any operating point (z0 , u0 ) using Taylor’s series expansion as follows

z(k + 1) = f (z0 (k), u0 (k)) + A(z(k) − z0 (k))

+B(u(k) − u0 (k)) + z,ho

y(k) = g(z0 (k)) + C(z(k) − z0 (k)) + y,ho

∂f
" #
A = (z(k), u(k))
∂z z0 (k),u0 (k)
∂f
" #
B = (z(k), u(k))
∂u z0 (k),u0 (k)
∂g
" #
C = (z(k), u(k))
∂z z0 (k),u0 (k)

where A, B and C represent the partial derivatives of the ELM models in the Taylor’s expansion, z,ho and y,ho represent

5 These functions are modeled using ELM as in the previous section

12
higher order terms. The linearized model can be further written as

z(k + 1) = Az(k) + Bu(k) + d1 (k)

y(k) = Cz(k) + d2 (k) (12)

d1 = f (z0 (k), u0 (k)) − (Az0 (k) + Bu0 (k)) + z,ho

d2 = g(z0 (k), u0 (k)) − Cz0 (k) + y,ho .

The following subsections describe the calculation of the system matrices A, B and C followed by the MPC optimiza-
tion problem formulation.

4.1. Calculation of System Matrices

The matrices A, B and C in equation (12) can be determined by exploiting the structure of the ELM model as
follows. Let the augmented input vector (as in equation (2)) to ELM be given by x(k) = [u(k), z(k)]T ∈ Rn+m . The
matrices A, B and C can be determined by calculating the Jacobian of f (.) and g(.) with respect to the augmented input
vector x(k). The ELM model structure can be expressed as

ẑ(k + 1) = W T [ψ(WrT x(k) + br )] (13)

where ψ represents the hidden layer activation function and Wr ∈ Rn+m×nh , W ∈ Rnh ×n represents the input and
output layer parameters respectively. Here, nh represents the number of hidden neurons of the ELM model, φ(k) =
ψ(WrT x(k) + br ) ∈ Rnh ×1 represents the hidden layer output matrix as termed in literature (see Fig. 6). As mentioned
earlier, the matrix Wr consists of randomly assigned elements that maps the input vector to a high dimensional feature
space while br ∈ Rnh is a bias component assigned in a random manner similar to Wr . A sigmoidal activation function
is used that has been identified to work well.
Note that the ELM model in (13) is defined for inputs and outputs normalized to lie between [−1, +1]. Expressing
the model in (13) along with the normalization and de-normalization terms,
z − z  n
max min
o
ẑ(k + 1) = zmin + 1 + W T φ(k)
2
1
φ(k) =   
x(k)−xmin
  . (14)
T
− W 2 −1 +br
1 + e r xmax −xmin
∂f
Then, Jacobian matrix ∂x can be derived (ignoring the time index k) as
∂ f  zmax − zmin  T ∂φ
= W (15)
∂x 2 ∂x
where      
x−xmin
T
− Wr,1 2 xmax −x −1 +br,1
∂φi 2Wr ( j, i)e min
=  !2 ,
∂x j
    
x−xmin
T 2
− Wr,1 −1 +br,1
(xmax, j − xmin, j ) 1 + e xmax −x
min

13
 ∂φ
Wr,i represents the ith column of Wr and br,i represents the ith element of br , ∂x ∈ Rnh ×n+m . Since the augmented vector
∂f
is defined as x(k) = [u(k), z(k)]T , the matrices A and B can be extracted from the Jacobian ∂x as

∂f
" #
= [Bn×m |An×n ]. (16)
∂x n×n+m

A similar jacobian calculation can be done for g(.) and also for other type of activation functions ψ(.). The above
calculations are algebraic and can be efficiently performed online.

4.2. MPC Optimization Problem

The goal of MPC is to force the system output y(k) track a given reference r(k) and also penalize any large
excursion in the input signal u(k). This is achieved by solving the following optimization problem at every time
instant k (naturally corresponds to every combustion cycle of the engine).
Ny
X N
X u −1

min J(k) = kr(k + j|k) − y(k + j|k)k2Q1 + k∆u(k + j|k)k2Q2 (17)

j=1 j=0


umin ≤ u(k + j|k) ≤ umax









subjected to ∆umin ≤ ∆u(k + j|k) ≤ ∆umax (18)






ymin ≤ y(k + j|k) ≤ ymax




where r(k + j|k), y(k + j|k), ∆u(k + j|k) represents the reference, system output and control increment respectively. The
argument (k + j|k) indicates the signal from time index k until k + j being used for solving the optimization problem
at time index k. Here ∆u(k + j|k) = u(k + j|k) − u(k − 1 + j|k), Ny and Nu are prediction horizon (Ny ≥ 1) and control
horizon (0 < Nu ≤ Ny ) respectively, k.kQ1 and k.kQ2 denote weighted Euclidean norm weighted by matrices Q1 and Q2
respectively. The minimum and maximum constraints for u, ∆u and y are given as in equation (18).
By defining the following vectors,

Y(k) = [y(k + 1|k), y(k + 2|k), .., y(k + Ny |k)]T

∆U(k) = [∆u(k|k), ∆u(k + 2|k), .., ∆u(k + Nu − 1|k)]T

R(k) = [r(k + 1|k), r(k + 2|k), .., r(k + Ny |k)]T

and calculating y(k + j|k) recursively using equation (12), the vector Y(k) can be expressed as

Y(k) = Z(k)z(k) + U(k)∆U(k) + V(k)u(k − 1) + D1 (k)d1 (k) + D2 (k)d2 (k) (19)

14
where  

 CB . . 0 

CB + CAB . . .
 
 
 
. . . .
 
 
U =  CB + .. + CANu −1 B
 
. . CB 


 CB + .. + CA u B . . CB + CAB
 N


 
. . . .
 
 

CB + .. + CANy −1 B . . CB + .. + CANy −Nu B


   

 CB 

 CA



CB + CAB 
   
  CA2 
   
V =  .
 , Z = 
  .


   
. .
   
   

CB + CANy −1 B
  
CANy
 
   

 C   I p 
  
 C + CA 

I
 
 p 
  
D1 =   , D2 =  . .

.

 
 

   
.  . 
 
   
  
C + CANy −1 Ip


Q1 ∈ RNy p×Ny p , Q2 ∈ RNu m×Nu m , Y ∈ RNy p×1 ,

∆U ∈ RNu m×1 , R ∈ RNy p×1 , U ∈ RNu p×Nu m ,
V ∈ RNy p×m , Z ∈ RNy p×n , D1 ∈ RNy p×n , D2 ∈ RNy p×p .
The optimization problem in equation (17) can now be expressed in vector form as
n o
min ∆Y T (k)Q1 ∆Y(k) + ∆U T (k)Q2 ∆U(k) (20)
∆U(k)

where

∆Y(k) = R(k) − Y(k)

= R(k) − Z(k)z(k) − U(k)∆U(k)

−V(k)u(k − 1) − D1 (k)d1 (k)

−D2 (k)d2 (k). (21)

Expressing as a quadratic programming problem, the above formulation reduces to

( )
1
min ∆U (k)W1 (k)∆U(k) + W2 (k)∆U(k)
T T
(22)
∆U(k) 2

subjected to E(k) ≤ F(k) (23)

15
where

W1 (k) = 2[U T (k)Q1 U(k) + Q2 ]

W2 (k) = −2U T (k)Q1 [R(k) − Z(k)z(k)

−V(k)u(k − 1) − D1 (k)d1 (k) − D2 (k)d2 (k)]

E(k) = [INu m − INu m H − H U(k) − U(k)]T

 

 ∆Umax 

 
 −∆Umin 
 
Umax − Uk−1
 
 
 

 −{Umin − Uk−1 } 

 
=
 
F(k)  
 Ymax − Z(k)z(k) − V(k)u(k − 1) 
 
 
 −D1 (k)d1 (k) − D2 (k)d2 (k) 

 
 
 

 −{Ymin − Z(k)z(k) − V(k)u(k − 1) 
 
 
−D1 (k)d1 (k) − D2 (k)d2 (k)}
 
 Im 0 . . 0 
 
 Im Im . . 0 
 
 
H =  .
 . . . . 

 . . . . . 
 
 
Im Im . . Im
 

Umin = [umin umin ..umin ]T

Umax = [umax umax ..umax ]T

∆Umin = [∆umin ∆umin ..∆umin ]T

∆Umax = [∆umax ∆umax ..∆umax ]T

Uk−1 = [u(k − 1) u(k − 1)..u(k − 1)]T

Ymin = [ymin ymin ..ymin ]T

Ymax = [ymax ymax ..ymax ]T .

Im , I p and INu m represents identity matrices in Rm×m , R p×p and RNu m×Nu m respectively. The above problem can be
solved using efficient quadratic programming solvers (Maciejowski, 2002). The first value of the optimal control
increment ∆U(k) is applied at the present time instant k and the optimization is solved again for the next time index.
This is done to handle any model inaccuracies and external disturbances.
MPC technology is well developed for linear systems with computationally efficient solvers. However, application
of MPC to a nonlinear system such as the HCCI engine involves solving several quadratic programming subproblem
16
sequentially at every time instant. To give an idea of the engine sampling period, an engine operating at 2500 RPM has
about 48 milliseconds to solve the quadratic program problem. The non-availability of efficient solvers for sequential
quadratic programming problem in real-time and the inability to guarantee a global optimal solution for such problems
constitute the major bottleneck for nonlinear MPC techniques. However, if the system is mildly nonlinear, a linearized
model (such as the one proposed in this paper) can be used to obtain suboptimal control laws in real-time using
quadratic programming. Also, when the quadratic programming (QP) subproblem (22) is strictly convex, a more
efficient and fast QP algorithm can be used (Effati and Ranjbar, 2011). The fast QP algorithm is derived as follows.

4.3. Fast QP-algorithm

Restating the MPC quadratic programming problem from (22),

( )
1
min ∆U T (k)W1 (k)∆U(k) + W2T (k)∆U(k) (24)
∆U(k) 2

subjected to E(k)∆U(k) − F(k) ≤ 0 (25)

where W1 (k) is a symmetric positive definite matrix and the objective function is strictly convex. The lagrangian dual
problem can be formulated as
1 
max inf ∆U T (k)W1 (k)∆U(k) + W2T (k)∆U(k) + λL (E(k)∆U(k) − F(k)) (26)
λL 2

where λL ≥ 0. Taking the derivative with respect to ∆U,

W1 (k)∆U(k) + W2 (k) + E T (k)λL = 0. (27)

The dual problem can be expressed as

1 
max ∆U T (k)W1 (k)∆U(k) + W2T (k)∆U(k) + λL (E(k)∆U(k) − F(k)) (28)
λL 2

subjected to
W1 (k)∆U(k) + W2 (k) + E T (k)λL = 0

y ≥ 0.

Since W1 (k) is positive definite, W1 (k)−1 exists and equation (27) can be solved as follows

∆U ∗ = −W1 (k)−1 (W2 (k) + E T (k)λL ). (29)

Following a direct substitution of (29) in (28), the problem becomes

( )
1 T 1 T
max λL Λ1 λL + λL Λ2 − W2 (k)W1 (k)W2 (k)
T −1
(30)
λL 2 2

subjected to y ≥ 0 (31)

17
where Λ1 = −E(k)W1−1 (k)E T (k) and Λ2 = −F(k) − E(k)W1−1 (k)W2 (k). The problem in (30) can be solved using a
simple gradient ascent method as follows.

λLk+1 = λLk + λ step (Λ1 λL + Λ2 ). (32)

In order to satisfy the constraint in (31), the following modification is made

λLk+1 = max(λLk + λ step (Λ1 λL + Λ2 ), 0) (33)

where λ step defines the step size. The solution of (33) gives the optimal λ∗L which can be substituted in (29) to get the
optimal MPC output ∆U ∗ . The strictly convex property of the MPC problem is made use of in solving the quadratic
programming subproblem efficiently.

5. Application to Control of a Gasoline Homogeneous Charge Compression Ignition Engine

The goal of this section is to design a real-time control algorithm that can control the HCCI engine to track a
given reference command in IMEP and CA50. The IMEP reference represents the driver’s load demand while CA50
reference is calculated a priori for achieving a sweet spot between high efficiency, low emissions and stability in
combustion. The procedure for calculating reference IMEP and CA50 is outside the scope of this article and the
reader is referred to (Killingsworth et al., 2009; Popovic et al., 2006). In this paper, it is assumed that the reference
commands are available and the goal is to determine the optimal control trajectories of FM, EVC and SOI that tracks
the reference commands.
Let us recap the overall control framework shown in Fig.1. An offline trained nonlinear model of the HCCI
engine, developed using ELM in section 3, is used as a proxy for the HCCI engine. At every operating point, the
model is linearized and the system matrices are used to formulate a quadratic program as described in section 4. The
engine hardware limitations as well as operating constraints are included in the optimization problem. The reference
command includes optimal set points of IMEP and CA50 that is to be tracked by the engine. At every time instant, the
solver receives the present state information (z(k)) from sensor measurements and attempts several control trajectories
of FM, EVC and SOI. The optimal trajectory that minimizes the tracking error and control excursion is given to the
engine as input.
As mentioned previously, the control knobs for the HCCI engine such as FM, EVC and SOI vary the amount of
EGR trapped in the cylinder. The EGR influences the temperature and concentration of the mixture, affecting IMEP
and CA50 of the combustion event. The physical states of the system can include a combination of temperature,
concentration, pressure of the gas mixture before combustion but a direct measurement of in-cylinder states is not
feasible for production engines. In-cylinder pressure is the only available measurement from which input-output
models were built in section 3. As a consequence, for the control problem considered here, the states, inputs and

18
outputs of the system are given by

z(k) = [I MEP(k) CA50(k) Pmax (k)

Rmax (k) T b (k) λ(k)]T ∈ R6

u(k) = [F M(k) EVC(k) S OI(k)]T ∈ R3

y(k) = [I MEP(k) CA50(k)]T ∈ R2 .

The dynamic equations of the model can be expressed as in (11). This system represents a nonlinear multi-input
multi-output problem. The step references are designed to vary between 2.6 and 3.2 bar IMEP and between -4 and -10
deg CA50 defined in degrees before combustion TDC. This covers most of the naturally aspirating HCCI operating
range at 1800 RPM. The MPC constraints are defined as follows.

umin = [19 − 121 272]T

umax = [25 − 100 375]T

∆umin = [−6 − 22 − 103]T

∆umax = [6 22 103]T

ymin = [2.1 − 14]T

ymax = [3.55 − 2]T .

It should be noted that the objective function and constraints can be modified in a straightforward manner to include
more performance metrics such as emissions, combustion roughness, knock levels etc.
The HCCI and CA50 measurements from the engine are extremely noisy and has different characteristics. CA50 is
characterized by a highly random process involving noise in the in-cylinder pressure measurements, variability in gas
mixing, heat transfer, difference in cylinder-to-cylinder design and cycle-to-cycle variations (Thring, 1989; Kalghatgi
and Head, 2006; Bengtsson et al., 2004). As a consequence, the variability in CA50 is high. The IMEP is a more
averaged effect of the pressure characteristics in the cylinder and thus its variability is lesser. The noise distribution
of IMEP and CA50 at certain representative operating conditions are summarized in Fig. 8. It can be seen that the
noise is almost Gaussian with zero mean. The top rows represent noise in CA50 signals at different IMEP conditions
while the bottom rows represent noise in IMEP measurements at different CA50 conditions. The MPC framework
considered in this work involves linearization of the models at every operating point defined by sensor measurements.
In order to simulate the controller’s robustness to the high amplitude noise seen in HCCI engine measurements, white
noise signals of zero mean and variances 0.0012 and 1.76 are injected at the outputs of the IMEP and CA50 model
respectively and the models are linearized around the noisy measurements experienced in real engine operation.

5.1. Results And Discussion

In this section, the results of the ELM based MPC applied to the HCCI problem is discussed. A simulation is
conducted with a Matlab real-time target bypassing communication with the engine ECU. The ELM model (obtained
19
 IMEP =3.26 IMEP =3 IMEP =2.63 IMEP =2.88
150 40 15 40

30 30
100 10
20 20
50 5
10 10

0 0 0 0
−10 0 10 −5 0 5 −5 0 5 −5 0 5
CA50 =−8.93 CA50 =−7.45 CA50 =−7.87 CA50 =−8.55
10 20 15 40

15 30
10
5 10 20
5
5 10

0 0 0 0
−0.1 0 0.1 −0.2 0 0.2 −0.2 0 0.2 −0.2 0 0.2

Figure 8: Noise characteristics of CA50 (top rows) and IMEP (bottom rows) as observed in the HCCI engine experimental data for different
operating conditions.

from Section 3.3) is observed to compactly represent the HCCI engine with about 320 parameters (Wr ∈ Rn+m×nh , br ∈
Rnh , W ∈ Rnh ×n ). Further, the fast quadratic programming approach enabled online implementation suitable for the
engine’s high sampling rates. At every simulation step, the inputs and outputs of the engine model are available
as sensor measurements to the MPC system. The HCCI model is linearized around the sampled states and is used
in the online MPC algorithm to determine the optimal control change from its set-point (∆U ∗ ); i.e., the optimal
fuel mass, EVC and SOI to be given to the engine to track the reference commands. After performing several trial
and error simulations, the prediction and control horizons take values Ny = 3 and Nu = 3, tuned for minimum
tracking error. Similarly, the gain matrices are tuned to be Q1 = diag([500 1 500 1 500 1]T ) and Q2 =
diag([20 1 1 20 1 1 20 1 1]T ). The odd diagonal values of Q1 correspond to reducing IMEP tracking
error while the even diagonal values correspond to reducing CA50 tracking error. The relative difference in values
between IMEP and CA50 indicates the relatively high noise level of CA50 and a smaller weight trying to track the
noisy CA50 less aggressively. Q2 is tuned to avoid saturation of the input signals (taken as constraints in the MPC
formulation). The gain corresponding to FM is high indicating a large excursion in FM is penalized so that fuel mass
is slowly varied. This is done to ensure that both FM and SOI from reaching their actuator limits.
To evaluate the MPC performance, a reference command involving a series of random steps of IMEP and CA50
is considered. The step input to the controller can be thought of as shifting the engine from one set point to another
and gives a good idea of the control performance in terms of response speed, transient and steady state performances.
It can be observed from Fig. 9a that the controller is able to track the reference command. The time taken by the
20
 System State
Tracking performance Reference
IMEP

3 24
(bar)

Fuel Mass
2.8

(mg/cyc)
22
0 500 1000 1500 2000
Tracking performance
(deg bTDC)

−5 20
CA50

−10 0 500 1000 1500 2000

0 500 1000 1500 2000

Exhaust Valve Closing

40 −100
max
(bar)

35

(deg aCTDC)
−105
P

30
0 500 1000 1500 2000 −110

4 −115
(bar/deg)
Rmax

3 −120
2
0 500 1000 1500 2000
0 500 1000 1500 2000

30 350
Torque
(Nm)

Start of Injection
25

(deg aETDC)
0 500 1000 1500 2000
345
EAFR

1.3
(−)

1.25
1.2
340
0 500 1000 1500 2000 0 500 1000 1500 2000
Combustion Cycles Combustion Cycles

(a) State trajectories of the HCCI engine model (with noise) (b) Control trajectories of the HCCI engine model (with
using MPC control (case 1: constraint on Rmax not enforced noise) using MPC control (case 1 constraint on Rmax not en-
in the MPC formulation). forced in the MPC formulation). The upper and lower limits
of each actuator is shown in dotted red.

System State
Tracking performance Reference
IMEP

3 24
(bar)

Fuel Mass
(mg/cyc)

2.8
22
0 500 1000 1500 2000
Tracking performance
(deg bTDC)

−5 20
CA50

−10 0 500 1000 1500 2000

0 500 1000 1500 2000
Exhaust Valve Closing

−100
max
(bar)

35
(deg aCTDC)

−105
P

30
0 500 1000 1500 2000 −110

4 −115
(bar/deg)
Rmax

3 −120
2
0 500 1000 1500 2000
0 500 1000 1500 2000

350
Torque

28
(Nm)

26
Start of Injection

345
(deg aETDC)

24
0 500 1000 1500 2000 340

335
EAFR

1.3
(−)

1.25
1.2 330
1.15
0 500 1000 1500 2000 0 500 1000 1500 2000
Combustion Cycles Combustion Cycles

(c) State trajectories of the HCCI engine model (with noise) (d) Control trajectories of the HCCI engine model (with
using MPC control (case 2 constraint on Rmax enforced in the noise) using MPC control (case 2 constraint on Rmax enforced
MPC formulation). in the MPC formulation). The upper and lower limits of each
actuator is shown in dotted red.

Figure 9: Evaluation of MPC controller using HCCI engine simulator.

21
controller to switch the engine from one set point to the next is about 5 cycles for a small step, while it takes about 20
cycles for a larger step. By tuning the MPC gains further, the controller response can be shaped to tradeoff between
quickness in response to overshoot and minimum oscillations around the reference. Care must be taken while tuning
the gains so that the controller does not saturate the actuators. Also from Fig. 9a, the transients of the other states
of the engine including Pmax , Rmax , engine brake torque and fuel richness in the mixture can be observed. These
quantities are obtained from the HCCI engine model for the optimal MPC control trajectories. The control inputs
(optimal trajectories) of MPC such as FM, EVC and SOI can be seen in Fig. 9b. It can be observed that none of the
actuators violate the specified limits.
In the above simulation, the constraints are specified only for the actuators and output variables such as IMEP
and CA50. However in real situations, typically the engine operation is constrained between stability limits such
as misfire, knock and ringing. In order to analyze MPC performance with state constraints and stability limits, an
additional constraint is defined on Rmax . In Fig. 9a, the maximum rate of pressure rise (Rmax ) exceeds 3.5 bar/deg
CA(shown by red dotted line). By including this as a state constraint in the MPC formulation, the second simulation
is performed. The MPC performance is summarized in Fig. 9c and Fig. 9d. It can be observed that the MPC control
is modified so that the engine Rmax is less than 3.5 bar/deg CA. A direct comparison of the MPC control trajectories
in Fig. 9b and Fig. 9d shows that when the constraint on Rmax is active, MPC reduces the SOI command to operate
the engine with less pressure rise rates for the given IMEP and CA50 reference commands. This summarizes the
capability of the proposed approach to compute an optimal solution for the multi-input multi-output control problem.
It has to be noted that, as the constraint is active, it becomes difficult for MPC to achieve close tracking as the available
freedom is restricted by the active constraint. By tuning the gains further, a better tracking may be achieved.
Finally, in order to simulate the engine demands in a vehicle application, a slowly varying sinusoidal reference
command is given both for IMEP and CA50. This is to simulate the engine power demand while accelerating and
slowing down repeatedly. The control performance is summarized in Fig. 10a and Fig. 10b. It can be observed
that the MPC tracks the given reference to a good degree of accuracy with constraint satisfaction. Hence, the above
simulation cases demonstrate the working of ELM based MPC, its manageable computational complexity and its
overall potential for real time HCCI engine control.

6. Conclusions and Future Work

HCCI engine control is a nonlinear, multi-input multi-output problem with state and actuator constraints which
calls for advanced control designs. We formulated and developed an integrated model based control framework and
used machine learning to model the complex combustion behavior. ELM models are shown to compactly represent
the nonlinear dynamics of HCCI and can be used for predicting several steps ahead of time. Model linearization
is facilitated by analytically calculating the gradients using the structure of ELM models. The MPC optimization is
formulated as a convex problem for which a fast quadratic programming method is used. Simulation studies have been

22
 System State
Tracking performance Reference
3.2
IMEP

24
(bar)

Fuel Mass
(mg/cyc)
2.8
2.6 22
0 500 1000 1500 2000
Tracking performance
(deg bTDC)

−5 20
CA50

−10 0 500 1000 1500 2000

0 500 1000 1500 2000

Exhaust Valve Closing

−100
Pmax
(bar)

35

(deg aCTDC)
−105
30
0 500 1000 1500 2000 −110
4 −115
(bar/deg)
Rmax

3
−120
2
0 500 1000 1500 2000
0 500 1000 1500 2000

30
Torque

355
(Nm)

Start of Injection
25

(deg aETDC)
350
0 500 1000 1500 2000

1.3 345
EAFR
(−)

1.2 340
0 500 1000 1500 2000 0 500 1000 1500 2000
Combustion Cycles Combustion Cycles

(a) State trajectories of the HCCI engine model (with noise) using (b) Control trajectories of the HCCI engine model (with noise) us-
MPC control (case 3: sinusoidal reference commands). ing MPC control (case 3: sinusoidal reference commands). The
upper and lower limits of each actuator is shown in dotted red.

Figure 10: Evaluation of MPC controller using HCCI engine simulator.

conducted to demonstrate the working of ELM based MPC for the considered HCCI engine problem. For future, we
focus on extending the method to develop controller for the experimental engine and reformulating MPC to include
vehicle level objectives such as fuel economy and emissions.

Acknowledgment

This material is based upon work supported by the Department of Energy and performed as a part of the ACCESS
project consortium (Robert Bosch LLC, AVL Inc., Emitec Inc.) under the direction of PI Hakan Yilmaz, Robert
Bosch, LLC. X. Nguyen is supported in part by NSF Grants CCF-1115769 and ACI-1047871.

Disclaimer

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither
the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or
implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any infor-
mation, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.
Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or

23
otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States
Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or
reflect those of the United States Government or any agency thereof.

References

Azad, N.L., Mozaffari, A., Hedrick, J.K., 2015. A hybrid switching predictive controller based on bi-level kernel-based {ELM} and online trajectory
builder for automotive coldstart emissions reduction. Neurocomputing, –.
Bengtsson, J., Strandh, P., Johansson, R., Tunestal, P., Johansson, B., 2004. Closed-loop combustion control of homogeneous charge compression
ignition HCCI engines dynamics. Int. J. Adaptive Control and Signal Processing 18, 167–179.
Bengtsson, J., Strandh, P., Johansson, R., Tunestal, P., Johansson, B., Oct 2006a. Model predictive control of homogeneous charge compression
ignition (HCCI) engine dynamics. In: Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications,
2006 IEEE International Symposium on Intelligent Control, 2006 IEEE, pp. 1675–1680.
Bengtsson, J., Strandh, P., Johansson, R., Tunestl, P., Johansson, B., 2006b. Hybrid control of homogeneous charge compression ignition (HCCI)
engine dynamics. International Journal of Control 79 (5), 422–448.
Effati, S., Ranjbar, M., 2011. A novel recurrent nonlinear neural network for solving quadratic programming problems. Applied Mathematical
Modelling 35 (4), 1688 – 1695.
EPA, 2012. Light-duty automotive technology, carbon dioxide emissions, and fuel economy trends: 1975 through 2011. Executive Summary,
United States Environmental Protection Agency.
Franko, S., Koç, I.M., Özsoy, C., Sari, N., 2011. Mpc and lqr type controller design and comparison for an unmanned helicopter. In: Proceedings
of the 2011 Summer Computer Simulation Conference, SCSC ’11, Society for Modeling & Simulation International, Vista, CA, pp. 138–144.
Huang, G.-B., Zhou, H., Ding, X., Zhang, R., 2012. Extreme learning machine for regression and multiclass classification. IEEE Transactions on
Systems, Man, and Cybernetics, Part B 42 (2), 513–529.
Huang, G.-B., Zhu, Q.-Y., Siew, C.-K., 2006. Extreme learning machine: Theory and applications. Neurocomputing 70, 489–501.
Jade, S., Hellström, E., Jiang, L., Stefanopoulou, A.G., 2012. Fuel governor augmented control of recompression HCCI combustion during large
load transients. In: Proc. American Control Conference, pp. 2072–2077.
Janakiraman, V.M., Nguyen, X., Assanis, D., 2013. Nonlinear identification of a gasoline HCCI engine using neural networks coupled with principal
component analysis. Applied Soft Computing 13 (5), 2375 – 2389.
Janakiraman, V.M., Nguyen, X., Sterniak, J., Assanis, D., 2014. A system identification framework for modeling complex combustion dynamics
using support vector machines. In: Informatics in Control, Automation and Robotics. vol. 283 of Lecture Notes in Electrical Engineering.
Springer International Publishing, pp. 297–313.
Janakiraman, V.M., Nguyen, X., Sterniak, J., Assanis, D., Jan 2015. Identification of the dynamic operating envelope of HCCI engines using class
imbalance learning. Neural Networks and Learning Systems, IEEE Transactions on 26 (1), 98–112.
Kalghatgi, G.T., Head, R.A., 2006. Combustion limits and efficiency in a homogeneous charge compression ignition engine. Int. J. Engine Res 7,
215–236.
Killingsworth, N.J., Aceves, S.M., Flowers, D.L., Espinosa-Loza, F., Krstic, M., 2009. HCCI engine combustion-timing control: Optimizing gains
and fuel consumption via extremum seeking. Control Systems Technology, IEEE Transactions on 17 (6), 1350–1361.
Li, D., Wang, D., 2013. Multivariable control of residual-affected hcci engines based on model predictive control. In: SAE Technical Paper 2013-
01-2511.
Maciej, L., 2009. Neural networks in model predictive control. In: Intelligent Systems for Knowledge Management. vol. 252 of Studies in Compu-
tational Intelligence. Springer Berlin Heidelberg, pp. 31–63.
Maciejowski, J.M., 2002. Predictive Control: With Constraint. Pearson Education. Prentice Hall.

24
Mozaffari, A., Azad, N.L., 2014. Optimally pruned extreme learning machine with ensemble of regularization techniques and negative correlation
penalty applied to automotive engine coldstart hydrocarbon emission identification. Neurocomputing 131, 143 – 156.
Mozaffari, A., Vajedi, M., Azad, N.L., 2015. A robust safety-oriented autonomous cruise control scheme for electric vehicles based on model
predictive control and online sequential extreme learning machine with a hyper-level fault tolerance-based supervisor. Neurocomputing 151,
Part 2, 845 – 856.
Narendra, K.S., Parthasarathy, K., Mar. 1990. Identification and control of dynamical systems using neural networks 1 (1), 4–27.
Nelles, O., 2001. Nonlinear System Identification: From Classical Approaches to Neural Networks and Fuzzy Models. Engineering Online Library.
Springer.
Nowak, R., Van Veen, B., april 1993. Nonlinear system identification with pseudorandom multilevel excitation sequences. In: Acoustics, Speech,
and Signal Processing, 1993. ICASSP-93., 1993 IEEE International Conference on, vol. 4, pp. 456 –459 vol.4.
Olsson, J.-O., Tunestal, P., Johansson, B., 2001. Closed-loop control of an HCCI engine. In: SAE Technical Paper 2001-01-1031.
Ortner, P., Del Re, L., May 2007. Predictive control of a diesel engine air path. Control Systems Technology, IEEE Transactions on 15 (3), 449–456.
Ortner, P., Langthaler, P., Garcia Ortiz, J.V., del Re, L., Oct 2006. MPC for a diesel engine air path using an explicit approach for constraint
systems. In: Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications, 2006 IEEE International
Symposium on Intelligent Control, 2006 IEEE, pp. 2760–2765.
Popovic, D., Jankovic, M., Magner, S., Teel, A., 2006. Extremum seeking methods for optimization of variable cam timing engine operation.
Control Systems Technology, IEEE Transactions on 14 (3), 398–407.
Ravi, N., Roelle, M.J., Liao, H.H., Jungkunz, A.F., Chang, C.F., Park, S., Gerdes, J.C., 2009. Model-based control of HCCI engines using exhaust
recompression. In: IEEE Transactions on Control Systems Technology.
Re, L., Allgöwer, F., Glielmo, L., Guardiola, C., Kolmanovsky, I., 2010. Automotive Model Predictive Control: Models, Methods and Applications.
Lecture Notes in Control and Information Sciences. Springer.
Ross, K., 1980. Elementary Analysis: The Theory of Calculus. Undergraduate Texts in Mathematics. Springer.
Shahbakhti, M., Koch, C.R., 2008. Characterizing the cyclic variability of ignition timing in a homogenous charge compression ignition engine
fueled with n-heptane/iso-octane blend fuels. Int. J. Engine Res 9, 361–397.
Thring, R., 1989. Homogeneous-charge compression-ignition engines. SAE paper 892068.
Widd, A., Liao, H.-H., Gerdes, J.C., Tunestål, P., Johansson, R., 2013. Hybrid model predictive control of exhaust recompression hcci. Asian
Journal of Control.
Wong, K.I., Wong, P.K., Cheung, C.S., Vong, C.M., 2013. Modelling of diesel engine performance using advanced machine learning methods
under scarce and exponential data set. Applied Soft Computing 13 (11), 4428 – 4441.
Wong, P., Wong, H., Vong, C., Xie, Z., Huang, S., 2014. Model predictive engine air-ratio control using online sequential extreme learning machine.
Neural Computing and Applications, 1–14.
URL [Link]
Wong, P.K., Wong, K.I., Vong, C.M., Cheung, C.S., 2015a. Modeling and optimization of biodiesel engine performance using kernel-based extreme
learning machine and cuckoo search. Renewable Energy 74, 640 – 647.
Wong, P.K., Zhong, J., Yang, Z., Vong, C.M., 2015b. Sparse bayesian extreme learning committee machine for engine simultaneous fault diagnosis.
Neurocomputing, –.
Yin, G., Zhang, Y.-T., Li, Z.-N., Ren, G.-Q., Fan, H.-B., 2014. Online fault diagnosis method based on incremental support vector data description
and extreme learning machine with incremental output structure. Neurocomputing 128, 224 – 231.

25