Neural Process Lett (2016) 44:857–886

DOI 10.1007/s11063-016-9500-7

A Novel Adaptive NARMA-L2 Controller Based on

Online Support Vector Regression for Nonlinear Systems

Kemal Uçak1 · Gülay Öke Günel 1

Published online: 3 February 2016

© Springer Science+Business Media New York 2016

Abstract In this study, a novel nonlinear autoregressive moving average (NARMA)-L2 con-
troller based on online support vector regression (SVR) is proposed. The main idea is to obtain
a SVR based NARMA-L2 model of a nonlinear single input single output system (SISO)
by decomposing a single SVR which estimates the nonlinear autoregressive with exogenous
inputs (NARX) model of the system. Consequently, using the obtained SVR-NARMA-L2
submodels, a NARMA-L2 controller is designed. The performance of the proposed SVR
based NARMA-L2 controller has been evaluated by simulations carried out on a bioreactor
system, and the results show that the SVR based NARMA-L2 model and controller attain
good modelling and control performances. Robustness of the controller in the case of system
parameter uncertainty and measurement noise have also been examined.

Keywords Adaptive control · NARMA-L2 controller · NARMA-L2 model · Online support

vector regression · SVR-NARMA-L2 controller

1 Introduction

The way to survival for living organisms is adaptation to the physical world. Similarly, con-
tinuation of the process in control systems depends on the adaptation skills of the controllers
to internal and external factors affecting the system. In nonlinear systems, determination of
the control input to cope with the nonlinear dynamics of the system is a challenging task.
Therefore, representation of complex systems via simple models is crucial to design effective
controllers. An ingenious model can be the key to obtain a high performance control method.
Nonlinear autoregressive moving average (NARMA)-L2 controller, introduced by Narendra

B Kemal Uçak
kemal.ucak@itu.edu.tr
Gülay Öke Günel
gulay.oke@itu.edu.tr
1 Department of Control and Automation Engineering, Faculty of Electrical-Electronics Engineering,
Istanbul Technical University, Maslak, 34469 Istanbul, Turkey

123
858 K. Uçak, G. Öke Günel

and Mukhopadhyay [1], is one of the most effective artificial neural network (ANN) controller
architectures for nonlinear systems. The central idea of this type of controller is to extricate
the control signal from nonlinear dynamics of the system model via Taylor expansion [2–4].
Thus, a nonlinear model is achieved with separate control signal term which is eluded from
nonlinear inner dynamics of the system model. The controller is simply a rearrangement of
neural network system model, which is trained offline in batch form [4,5]. An important
advantage of NARMA-L2 controller is that it does not require the design of extra network
topology or controller training as in model reference adaptive control (MRAC). NARMA-L2
controller has been successfully utilized to control various nonlinear systems in technical
literature. Majstorovic et al. [2] designed a NARMA-L2 controller for a two tank system
which is derived via offline NARMA model of the system. Pedro et al. showed a realistic
implementation of NARMA-L2 slip controller for an anti-lock braking system [3]. Hagan et
al. compared the performance of various controller structures based on ANNs on a continu-
ously stirred tank reactor (CSTR), a single link robot arm and a magnetic levitation system
[4,6]. They have controlled a catalytic CSTR via NARMA-L2 controller. Since NARMA-L2
controller produces an oscillatory control signal, Pukrittayakamee et al. [5] proposed to use a
term which consists of scaled reference and system output in order to alleviate oscillation and
chattering in control signal. Wahyudi et al. [7] smoothed NARMA-L2 controller to control
a single link robot arm since nonlinearities and parametric uncertainties are exceptionally
apparent in a single link manipulator. As underactuated systems are difficult to control and
can be asymptotically stabilized via continuous static feedback law, Akbarimajd and Kia [8]
utilized NARMA-L2 controller for a two-degree of freedom (DOF) underactuated planar
manipulator. Vesselenyi et al. [9] have achieved to control the position of a pneumatic actu-
ator via NARMA-L2 controller. The common feature of the works in literature is utilization
of offline training phase to obtain NARMA-L2 model of the system [1–9].
The controller performance in model based adaptive methods is directly influenced by
modeling inaccuracies. Nonlinear systems can be successfully modeled via various intelligent
modeling techniques such as ANN [10,11], adaptive neuro fuzzy inference system (ANFIS)
[10–12] and support vector regression (SVR) [13–15]. Since the network topologies which
are trained with backpropagation algorithm (ANN and ANFIS) may get stuck at local minima
resulting from their non-convex objective function, it is likely that they obtain system model
only locally. In order to reduce modelling inaccuracies and improve controller performance,
SVR based identification and control methods have been popularly applied in recent years due
to their non-linear prediction and generalization competency. Since the objective function in
SVR is convex, the gradient effects which occur in ANN and ANFIS vanish in SVR and global
extremum is ensured. Therefore, SVR-based controller structures [16,17] have recently been
used as leading model based adaptive controllers instead of ANN based approach, since they
yield a unique solution and possess powerful generalization ability [16,18].
Various controller structures based on SVR have been proposed to control nonlinear
systems in literature. These structures can be grouped under four main headings: adaptive
PID controllers, inverse controllers, internal model controllers (IMCs) and model predictive
controllers (MPCs).
Wanfeng et al. [19] utilized online least square SVR based on sliding window in order
to model Jacobian of the system which is required to obtain update rules for PID parame-
ters via gradient descent. Since the model based on sliding window forgets the transient
dynamics of the system in steady state, Zhao et al. [20] proposed to alter the size of the
sliding window to improve training performance of the model and also the controller per-
formance which depends on its Jacobian approximation competency. The kernel functions
used in SVR to provide nonlinearity have design parameters which have direct efficacy on

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 859

regression performance. Ucak and Oke [16] considered to tune bandwidth parameter of the
kernel function simultaneously with PID controller parameters via gradient descent method.
Yuan et al. [21] proposed a composite controller structure which consists of feedforward and
feedback parts. The feedforward part is derived using system Jacobian which is estimated by
SVR system model. Conventional PID controller is employed as feedback controller in order
to interfuse robustness to the controller structure against disturbances and estimated errors,
and robustness of the control structure has also been examined. Iplikci [17] offered a support
vector machine (SVM)-based PID controller in which the controller parameters are adjusted
via Levenberg–Marquardt algorithm by approximating K-step ahead system Jacobian. The
controller is formed of five components: classical incremental PID controller, ε-SVR non-
linear autoregressive with exogenous inputs (NARX) model of the system, line search block,
control signal correction block and controller parameter tuner. The Jacobian matrix utilized
in parameter tuner block is coalescence of correction block which includes K-step ahead
system Jacobian and a gradient vector which involves first order derivative of control signal
with respect to controller parameters. Control signal correction block produces a term which
is added to the computed control signal since the adjusted controller parameters may not be
adequate to compel the system output to track desired reference. The control signal correc-
tion term is derived via Taylor approximation of the control signal which includes correction
term. Golden section method has been employed in line search block to determine optimum
learning rate for the correction term. The main characteristics of these works is that SVR is
utilized to model both the system and the Jacobian. Unlike previously mentioned studies,
Takao et al. [22] used a support vector classifier (SVC) as a switching structure in a decision
tree framework to assign suitable controller parameters according to the operating point of
the system. The whole operating region is splitted into small uncertainty ranges via a priori
information. Robust PID controllers for each small range are constructed. Then, a single
SVC is designed for each small uncertainty range to determine whether the current operating
conditions of the system belong to that range or not. By combining all SVCs, a main decision
tree which covers the whole operating range of the system is constituted. Based on the final
decision, the appropriate controller for the current state of the system is deployed.
The main aim in control theory is to design a controller which identifies the inverse
dynamics of the system to be controlled as closely as possible. For this purpose, inverse
controllers which mimic inverse dynamics of the system can be designed using data sampled
approaches. The main stalemate of inverse controllers is that there may be no unique inverse
controller, moreover, there is no guarantee for the existence of the system inverse. Liu et al.
[23] employed a SVR to identify inverse dynamics of a nonlinear system and designed an
inverse controller based on this model. The control signal error which is also the approxi-
mation error of SVR is compensated via a PID controller. Wang et al. [24] used an online
SVR to design an inverse controller for a nonlinear system to be controlled. Yuan et al. [25]
utilized two SVR structures as the inverse controller and forward model of the system. The
parameters of controller and system model are optimized online via backpropagation algo-
rithm. The stability and convergence analysis of the controller and model are given in order
to guarantee convergence and fast learning.
Another model based control technique for nonlinear systems is nonlinear internal model
control (NIMC). NIMC is frequently utilized as a nonlinear control method owing to its dis-
turbance rejection capabilities [26] and robustness properties [27]. The drawbacks of NIMC
are that its implementation is convenient, it is restricted only for open-loop stable systems
[28] and it is based on the assumption that the system is reversible. Since the determination
of system model is the most significant stage of the design [26], NIMC performance depends
on accurate modelling. For this purpose, powerful generalization and modeling capacities

123
860 K. Uçak, G. Öke Günel

of SVR have been combined with IMC in [26,27,29,30]. Zhao et al. [26] employed LS-
SVR to avoid complex inverse controller. Sun and Song [29] suggested an adaptive IMC
for nonlinear systems by combining LS-SVR with sequential minimal optimization (SMO)
based pruning algorithm. Wang and Yuan [30] derived an approximate inverse control law
via Taylor expansion of NARMA model where SVR is used to approximate system Jacobian
and they also analyzed the stability of the proposed control law.
Model predictive control (MPC) can successfully be used to handle systems with large
time delay, non-minimum phase properties or unstable dynamics. MPC provides very robust
schemes [31,32] to control highly nonlinear dynamical systems compared to PID controllers
and can cope with structural changes [32–36]. All kinds of MPC are based on the same tactics:
using the future predictions of the system model, a set of future control signals is obtained
by solving, at each sampling instant, a finite-horizon open-loop optimal control problem and
then the first element of the control sequence is applied to the system [31,32]. In MPC, it is
aimed to update the set of control signals so that the system will follow the reference signal
closely by minimizing an objective function which includes future approximated tracking
errors and estimated control signals. For this purpose, a direction vector should be determined
to update the control signal vector. This direction vector is commonly composed of first and
second order derivatives of system output with respect to control inputs which are estimated
via intelligent model of the system. Iplikci [31,32], Zhiying and Xianfang [37] and Shin
et al. [38] deployed SVR based system models in MPC in order to overcome modelling
inaccuracies of highly non-linear systems. The performance of the MPC may devolve if the
model of the controlled system cannot be accurately computed. Iplikci [32] and Zhiying and
Xianfang [37] proposed to utilize online SVR system model in MPC framework in order to
cope with modelling inaccuracies resulting from disturbances or varying dynamics. Shin et
al. [38] proposed to update the parameters of the previously offline trained SVR model via
gradient descent algorithm to prevent deterioration in SVR system model.
In this paper, the strong modelling capability of SVR is combined with the functionality
of the NARMA-L2 controller structure to propose a novel online SVR based NARMA-L2
controller for nonlinear single input single output (SISO) dynamical systems. The NARMA
controller in general consists of two nonlinear submodels which are independent from current
control signal applied to the system. If ANN is used to estimate these submodels, a subnetwork
should be designed for each of them. The parameters of these subnetworks in NARMA based
on ANN can be easily adjusted without knowing the exact outputs of the submodels since
update rules for each parameter of the subnetworks can be derived by using backpropagation
algorithm and chain rule. However, estimation techniques based on ANN have a major
drawback, they are likely to get stuck at local minima and the models attained by ANN
will be valid only locally. This is the main rationale behind our choosing SVR to estimate
NARMA-L2 submodels. The primary advantage of SVR compared to backpropagation based
identification methods is that global extremum is assured, hence the system model prevalent
in all regions is obtained precisely. Nevertheless, if NARMA submodels are estimated by
SVR instead of ANNs, there is a major difficulty. To obtain the optimum mapping function
in SVR the input–output training pairs must be available, however outputs of NARMA
submodels are not known a priori, so estimating NARMA model with SVR is a hard task. In
this paper, we overcome this problem by deriving the NARX model of the system first, and
then obtaining NARMA submodels via NARX model of the system. Thus, it is achieved to
design NARMA-L2 controller using NARMA-L2 model of the system.
The main contribution of the paper is to convert an online SVR model of the system to
an online NARMA-L2 controller, without necessarily any prior information on submodel
outputs. The parameters which are utilized in conversion are optimized via Levenberg–

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 861

Marquardt algorithm by considering K-step ahead future behaviour of the system. The
performance of the proposed SVR controller has been evaluated by simulations carried out
on a bioreactor benchmark system. Robustness of the proposed controller has been examined
by adding measurement noise and parametric uncertainty to the system. The performance
of the controller has been compared with a SVM-based PID controller proposed by Iplikci
[17]. The results indicate that the online SVR-NARMA-L2 controller together with online
SVR model attain good modelling and control performances.
The paper is organized as follows: Sect. 2 describes the basic principles of NARMA-L2
model and controller. Online ε-SVR is summarized in Sect. 3. In Sect. 4, the proposed control
architecture is explained and the optimization problem to derive SVR submodels is given.
In Sect. 5, simulation results and performance analysis of the controller are given together
with an assessment of real time applicability. Also, the proposed method is compared with a
SVM-based PID controller. The paper ends with a brief conclusion in Sect. 6.

2 NARMA-L2 Model and Controller

NARMA-L2 model represents the dynamics of a nonlinear system exactly in a neighborhood

of the equilibrium state [1]. The output of the NARMA-L2 model is as follows:

ŷNARMAn+d = fˆn + ĝn u n , (1)

where fˆn = Fn (xn ) and ĝn = G n (xn ) are two nonlinear functions to be approximated and
xn = [u n−1 , . . . , u n−n u , yn , . . . , yn−n y +1 ]T stands for the current input feature vector of
the model where n u and n y emblematise the number of the past control inputs and system
outputs included in the feature vector, d indicates the relative degree. Fn and G n denote the
submodels computed by some numerical or intelligent estimator. The control signal u n is
successfully extricated from the nonlinearities in the model [2,6,39]. The main superiority
of this form is that it allows us to solve for the required control input which forces the system
output to track the desired reference signal [2,39]. As can be seen from (1), the next control
signal u n is separated from nonlinear dynamics. Thus, the control signal forcing the system
to follow the desired behaviour can be obtained via fˆn and ĝn sub-models. In ANN based
NARMA controllers, design of controller consists of two steps. The first step is to identify
the dynamics of the system to be controlled by training Fn (·) and G n (·) submodels. Second
step is to design NARMA controller and approximate u n using Fn (·) and G n (·) trained at
the previous step. The structure of the NARMA-L2 model for d = 1 is illustrated in Fig. 1.
The parameters of submodels Fn (·) and G n (·) are optimized to minimize the error between
system and model outputs. After identifying the dynamics of the system as in model given
in Fig. 1, NARMA-L2 controller can be designed by substituting reference signal (rn+1 ) in

T + yˆ NARMA n +1
yn D
xn
Fn (.) ∑
+
L

un
T
D
L
Gn (.) ×
Fig. 1 NARMA model

123
862 K. Uçak, G. Öke Günel

Online NARMA
Controller

rn +1 un yn +1
+ ∑ System
−
Fn (.) Gn (.)

xn

T
D
L

T
D
L

Fig. 2 NARMA controller

place of system output (yn+1 ) since the system output is compelled to follow the reference
signal. The control signal produced by the NARMA-L2 controller is obtained as follows:

rn+1 − fˆn
un ∼
= , (2)
ĝn

where fˆn = Fn (xn ) and ĝn = G n (xn ) are estimated outputs of submodels [6,40,41]. The
block diagram of NARMA-L2 controller is illustrated in Fig. 2. The advantage of NARMA-
L2 controller is that it is simple to implement and it requires only two submodels of the
system to be controlled without a separate controller [1,8,40]. However, the control input
computed by NARMA-L2 controllers is typically oscillatory and this is a major drawback
of this method [6]. Another weakness is the relatively complicated structure of NARMA-L2
compared to NARX model since two separate subblocks are used [40]. Also, if the inverse
model is locally unstable or near stability margin the design of NARMA controller is not
possible [40].

3 Online ε-Support Vector Regression

In regression, the aim is to approximate a function from a set of sample data. SVR, first
asserted by Vapnik et al. [42–44], is one of the most effective regression techniques amongst
data sampled modeling methods since it ensures global extremum. Let us consider a training
data set:
T = {xi , yi }i=1
N
xi ∈ X ⊆ R n , yi ∈ R, (3)

where xi is the ith input data and yi is the corresponding output value, N is the size of the
training data and n is the dimension of the input space. The regression function ( ŷi ) which
represents the relationship among training data set in (3) can be characterized via SVR model
in (4):
ŷi = w (xi ) + b, i = 1, 2, . . . , N , (4)

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 863

where w is the weight vector of the regressor in feature space F, (xi ) is the image of input
data in F and b is the bias term [31].
In ε-SVR, an optimal regression surface which represents the given training data as accu-
rately as possible is obtained. The optimal regressor is searched within a predefined ε-tube
characterized with Vapnik’s ε-insensitive loss function. The primal form of the optimization
problem is described as follows [31,45]:

1  N

min  P = w2 + C ξi + ξi , (5)
w, b, ξ, ξ 2
i=1
with the following constraints constituted via Vapnik’s ε-insensitive loss function
yi − w (xi ) − b ≤ ε + ξi ,
w (xi ) + b − yi ≤ ε + ξi , (6)
ξi , ξi ≥ 0 i = 1, 2, . . . , N ,
where P indicates primal objective function, ε is the upper value of tolerable error, ξ ’s and
ξ  ’s are the slack variables representing the deviation from ε tube [31,45]. Since the objective
function in (5) is non-convex with respect to primal variables (w, b, ξ, ξ  ) and the training
algorithm may get stuck at local minima, it is difficult to determine the optimal regressor para-
meters. Therefore, a Lagrangian function which enables to convert the optimization problem
to a new form where global minima is guaranteed can be derived via Lagrangian multiplier
method. The Lagrangian, which is the combined version of primal objective function and its
constraints, can be comprised as follows [31,45]:

1  N
  N
LP = w2 + C ξi + ξi + βi (yi − w (xi ) − b − ε − ξi )
2
i=1 i=1
(7)

N
  
N
 
+ βi w (xi ) + b − yi − ε − ξi + −ηi ξi − ηi ξi ,
i=1 i=1

by inserting non-negative Lagrange multipliers β, β  , η and η . First order optimality con-

ditions for Lagrangian (L P ) are given as

∂Lp  N
= 0 −→ w − βi  (xi ) = 0, (8)
∂w
i=1

∂Lp 
N
 
= 0 −→ βi − βi = 0, (9)
∂b
i=1
∂Lp
= 0 −→ C − βi − ηi = 0, i = 1, 2, . . . , N , (10)
∂ξi
∂Lp
= 0 −→ C − βi − ηi = 0, i = 1, 2, . . . , N . (11)
∂ξi
Lagrangian (L P ) function has to satisfy the above optimality conditions at the solution. If the
optimality conditions in (8–11) are substituted in (7), dual form of the ε-SVR optimization
problem which is a quadratic problem is attained as follows:

1   
N N N N
   
D= βi − βi β j − β j K i j + ε βi + βi − yi βi − βi , (12)
2
i=1 j=1 i=1 i=1

123
864 K. Uçak, G. Öke Günel

subject to
0 ≤ βi ≤ C,
0 ≤ βi ≤ C,
(13)

N
 
βi − βi = 0, i = 1, 2, . . . , N ,
i=1

where K i j = (xi )T (xi ) is the kernel value which represents the similarity among the xi
and xj training samples. The optimal values of the dual variables βi , βi in the constructed
dual optimization problem formulated in (12–13) are obtained via quadratic programming
(QP) solvers. The support vectors are the training data with Lagrange multipliers that have
λi = βi − βi = 0 [17,31,46]. Thus, regression function in (4) can be rewritten with respect
to the support vectors and the corresponding Lagrange multipliers as in (14):

ŷ(x) = λi K (xi , x) + b, λi = βi − βi . (14)
i∈SV

The optimization problem in (12–13) is suitable for offline batch learning algorithms. In
online learning, since every new added data changes the structure of the regression problem
consistently, it is required to derive incremental tuning rules for model parameters by ensuring
the Karush–Kuhn–Tucker (KKT) conditions. For this purpose, a new Lagrangian function
which is derived from dual optimization problem (12–13) can be formulated as

1     
N N N N
  
LD = βi − βi β j − β j K i j + ε βi + βi − yi βi − βi
2
i=1 j=1 i=1 i=1
(15)

N
  
N
  N
 
− δi βi + δi βi − u i (C − βi ) + u i C − βi − z βi + βi .
i=1 i=1 i=1

Thus, KKT optimality conditions for L D are obtained in (16) via first order necessary opti-
mality conditions [47–51]:

∂LD  N
   N N
= β j − β j K i j + ε − yi − δi + u i + z = 0,
∂βi
j=1 i=1 i=1

∂LD 
N
  
N 
N
=− β j − β j K i j + ε + yi − δi + u i − z = 0,
∂βi
j=1 i=1 i=1
 
() () () () () ()
δi ≥ 0, ui ≥ 0, δi βi = 0, ui C − βi = 0. (16)

According to KKT conditions, at most one of the βi and βi should be nonzero and both are
nonnegative [48]. KKT conditions allow the reformulation of SVM for regression by dividing
the whole training data set T into the three main subsets: error support vectors (E), support
vectors (S) and remaining samples (R) which are classified depending on their Lagrange and
margin values in (17) [47].
E = {i | |λi | = C, |h (xi )| ≥ ε} ,
S = {i | 0 < |λi | < C, |h (xi )| = ε} , (17)
R = {i | |λi | = 0, |h (xi )| ≤ ε} .

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 865

(a) (b)

Fig. 3 E, S and R subsets before (a) and after (b) training

Using the function to be approximated and the SVR, error margin can be defined as follows:


N
h (xi ) = f (xi ) − yi = λ j K i j + b − yi . (18)
j=1

Convergence conditions for online SVR algorithm are given in (19) depending on KKT
conditions (16) and error margin function (18)
h (xi ) ≥ ε, λi = −C,
h (xi ) = ε, −C < λi < 0,
−ε ≤ h (xi ) ≤ ε, λi = 0, (19)
h (xi ) = −ε, 0 < λi < C,
h (xi ) ≤ −ε, λi = C.
The geometrical interpretation of incremental learning and these three subsets are visualized
in Fig. 3. In the incremental online SVR algorithm, when a new training data xc is received,
its corresponding λc value is initially set to zero, which is later updated to Δλc [31]. Then,
the largest possible change Δλc is calculated while at the same time keeping the system at
the equilibrium state with respect to the new KKT conditions [31]. Depending on Δλc , the
bias and the Lagrange multipliers of the samples in support set (S) are adjusted. Detailed
information related to recursive algorithm can be achieved via [47,48,52].

4 NARMA Controller Based on Online SVR

4.1 Identification of the NARMA Model by SVR

The first step in the proposed control method is to obtain a NARMA model of the system.
Initially, from the available input–output data, a NARX model is easily computed. In order to
design a NARMA-L2 controller, the dynamics of the NARX model should be decomposed
into a NARMA model. Therefore two separate SVRs are designed as depicted in Fig. 4.
SVRNARX calculates the NARX model of the system from the input–output dataset of the
system, then SVRNARMA-L2 decomposes this model to a NARMA model in order to design
a convenient SVRNARMA-L2 controller. We have used online SVR for modelling so that the

123
866 K. Uçak, G. Öke Günel

(a)
T
yn D yˆ NARX n +1
L
SVRNARX
T
un D
L

(b)
SVRNARMA-L2
T + yˆ NARMA n +1
yn D
xn
Fn (.) ∑
+
L

un
T
D
L
Gn (.) ×

Fig. 4 Decomposition of SVRNARX model (a) to SVRNARMA-L2 model (b)

model and the controller can be adapted in accordance with the changing dynamics of the
system. Since the parameters of the model are updated in online manner, the regression
functions of submodels F(·) and G(·) may alter as the number of discrete time steps n
proceeds, therefore F(·) and G(·) are symbolized as Fn (·) and G n (·) as given in Fig. 4. The
greatest difficulty of SVR based modelling is the necessity of input–output data pairs. In
ANN based modeling, the weights related to Fn (·) and G n (·) subnetworks can be tuned
simply using backpropagation algorithm and chain rule, even if the output of Fn (·) and
G n (·) are not known exactly. In SVR, input–output data pairs must be available to find the
optimal regression function. Since the outputs of Fn (·) and G n (·) are not known exactly, SVR
submodels for Fn (·) and G n (·) cannot be designed separately. Therefore, the focus in this
study is to drive model parameters of SVRNARMA-L2 from the previously obtained SVRNARX
model of a system. To achieve this, firstly a main SVRNARX model is designed as in Fig. 4a,
we call this as the “main model”, then submodels given in Fig. 4b are decomposed depending
on the main model. In SVRNARMA-L2 controller, the aim is to find u n . However as can been
seen in (2), at time step n the model parameters of Fn (·) and G n (·) are not known. Therefore,
the outputs of these submodels can be approximated depending on the model parameters
obtained at the previous step (n − 1). Hence, the output of the submodels at prediction phase
can be expressed as follows:
fˆn− ∼
= Fn−1 (xn ) ,
(20)
ĝ − ∼
n = G n−1 (xn ) ,

where the subscript stands for the time index of the state vector and the superscript is used to
indicate whether the system model derived in time step n or n − 1 is used. More specifically,
“−” in (20) means that the models computed in the previous time step n − 1, are evaluated
with the current state vector xn to obtain fˆn− and ĝn− . Similarly, “+” in fˆn+ and ĝn+ indicates

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 867

that the model obtained in current step (n) has been utilized together with the current state
vector xn as follows:
fˆn+ ∼
= Fn (xn ) ,
(21)
ĝ + ∼= G n (xn ) .
n
In order to compute the control signal (u n ), the parameters of Fn (·) and G n (·) submodels
must be obtained. Therefore, the following decomposition step is crucial to derive these
parameters.

Modeling and Decomposition Step at Time Step n–1: Consider the SVRNARX model of
the system depicted in Fig. 4a. Regression function for this model is given by:

N
ŷNARXn = βi K (xi , xc ) + bβ . (22)
i=1

Since this model is not directly usable in feedback linearization, the main model (SVRNARX )
should be divided to submodels as in Fig. 4b. Let us assume that the corresponding regression
function for Fn−1 (xn−1 ) and G n−1 (xn−1 ) are as follows:

N
+ ∼
fˆn−1 = Fn−1 (xn−1 ) = αi K (xi , xn−1 ) + b f ,
i=1
(23)

N
+
ĝn−1 ∼
= G n−1 (xn−1 ) = θi K (xi , xn−1 ) + bg .
i=1

Thus, output of the SVRNARMA-L2 model is obtained as:

ŷNARMAn = Fn−1 (xn−1 ) + G n−1 (xn−1 ) u n−1
 N 
N 
= αi K (xi , xn−1 ) + b f + θi K (xi , xn−1 ) + bg u n−1
i=1 i=1 (24)

N

= αi + θi u n−1 K (xi , xn−1 ) + b f + bg u n−1 .

i=1

If SVRNARMA-L2 and SVRNARX models are matched in (25), a relation between submodels
of SVRNARMA-L2 and SVRNARX can be attained as in (26).
ŷNARXn ∼
= ŷNARMAn ,

N 
N

βi K (xi , xn−1 ) + bβ = αi + θi u n−1 K (xi , xn−1 ) + b f + bg u n−1 , (25)

i=1 i=1
βi = αi + θi u n−1 ,
bβ = b f + bg u n−1 . (26)
The following assumption has been proposed to utilize this relation effectively and to approx-
imate the parameters of the submodels by this approach.
Assumption: Let us assume that the following relations exist between submodels.
αi = μ1 (·)θi ,
(27)
b f = μ2 (·)bg .

123
868 K. Uçak, G. Öke Günel

Thus, the following equality is obtained:

Fn−1 (xn−1 ) = μ1 (·)G n−1 (xn−1 ) + [μ2 (·) − μ1 (·)] bg . (28)

As can be seen from (28), the derived SVRNARMA-L2 model depends only on G n−1 (xn−1 ) for
the given values of μ1 (·) and μ2 (·). Using (26) and (27), the bias and the Lagrange multipliers
of the SVRNARMA-L2 submodels Fn−1 (xn−1 ) and G n−1 (xn−1 ) are obtained with respect to
parameters of SVRNARX model as follows:
βi
θi = , αi = μ1 (·)θi ,
μ1 (·) + u n−1
(29)
bβ
bg = , b f = μ2 (·)bg .
μ2 (·) + u n−1
Thus, outputs of the internal dynamics (Fn−1 (xn−1 ), G n−1 (xn−1 )) required for con-
troller design can be calculated using main SVRNARX model. The output of the derived
SVRNARMA-L2 model is:

ŷNARMAn = Fn−1 (xn−1 ) + G n−1 (xn−1 ) u n−1

= fˆn−1
+ +
+ ĝn−1 u n−1 , (30)

T
xn−1 = u n−2 , . . . , u n−n u −1 , yn−1 , . . . , yn−n y .

This assumption is crucial in building the relationship between NARX and NARMA models
and hence making it possible to obtain SVRNARMA-L2 model from the main SVRNARX model.

4.2 Controller Design

Prediction Step, Calculation of u n at Time Step n: In this section, the control signal
applied to the system is calculated via submodels Fn−1 (·) and G n−1 (·) which are approxi-
mated using the trained SVRNARX model at previous step (n − 1). The corresponding output
of the submodels against current state vector (xn ) are represented as fˆn− and ĝn− as given
in (31).

T
xn = u n−1 , . . . , u n−n u , yn , . . . , yn−n y +1 ,
fˆn− ∼
= Fn−1 (xn ) , (31)
− ∼
ĝ = G n−1 (xn ) .
n

Thus, the control signal (u n ) produced by SVRNARMA-L2 controller can be calculated as:

rn+1 − fˆn−
un ∼
= , (32)
ĝn−
by means of the submodels obtained at previous step (Fn−1 , G n−1 ). Then, u n is applied
to the system and yn+1 is computed. So, input–output data pair can be arranged for next
training phase of main SVRNARX model. Accordingly, SVRNARX model of the system can
be utilized to derive a SVRNARMA-L2 controller law given in (32) without explicitly knowing
the outputs of the submodels F(·) and G(·) separately. The proposed controller is illus-
trated in Fig. 5 where yn+1 is the output of the system, rn+1 denotes reference signal and
u n is the control signal produced by the controller. It is anticipated that the tracking per-
formance of the control system against varying dynamics is going to improve by using the

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 869

SVR NARMA-L2 Model

yˆ n +1
T +
D
xn
Fn (.) ∑
+
L

T
D
L
Gn (.) ×
−
SVR NARMA-L2 Controller ∑
rn +1 + yn +1
un
∑ System
+
−
Fn (.) Gn (.)

xn

T
D
L

T
D
L

Fig. 5 SVRNARMA-L2 controller based on online SVR

online SVR models. The adaptive structure of the submodels has an impact on the controller
performance.

4.3 Adaptive Predictive SVRNARMA-L2 Controller

The tracking performance of the proposed controller depends on the parameters (μ1 (·), μ2 (·))
utilized to compute the outputs of submodels. We have employed a predictive structure
which takes into account the K-step ahead future behaviour of the system to adjust μ1 (·) and
μ2 (·) parameters effectively. The optimization technique proposed by Iplikci for SVM-based
PID controller [17] has been utilized to optimize μ1 (·) and μ2 (·) parameters. The objective
function given in (33) is optimized via Levenberg–Marquardt algorithm to tune μ1 (·) and
μ2 (·).
1 2
K
1
2
E(μ) = ên+k + λ u n − u n−1 , (33)
2 2
k=1

where ên+k = rn+k − ŷn+k . The parameters can be adjusted using Levenberg–Marquardt
algorithm as follows:

   
μ1new μ1old
−1 T
= + JT J + ηI J e, (34)
μ2new μ2old

123
870 K. Uçak, G. Öke Günel

Fig. 6 Adaptive SVRNARMA-L2 controller based on online SVR

where ⎡ ⎤ ⎡ ∂ ŷn+1 ∂u
∂en+1 ∂ ŷn+1 ∂en+1 ∂ ŷn+1 ŷn+1 ∂u n ⎤
∂ ŷn+1 ∂μ1 ∂ ŷn+1 ∂μ2 − n
− ∂∂u
⎢ ⎥ ⎢ ∂u n ∂μ1 n ∂μ2
⎥
⎢ .. .. ⎥ ⎢ .. .. ⎥
⎢ . . ⎥ ⎢ . . ⎥
⎢
J = ⎢ ∂en+K ∂ ŷn+K ⎥
∂en+K ∂ ŷn+K ⎥ = ⎢ ∂ ŷn+K ∂u n
⎢ ∂ ŷn+K ∂u n ⎥⎥
⎢ ∂ ŷn+K ∂μ1 ⎥ − −
⎦ ⎣ √∂u n ∂μ1 √∂u n ∂μ2 ⎦
∂ ŷn+K ∂μ
⎣ √ √
2

∂ λΔu n ∂ λΔu n ∂u n ∂u n
λ ∂μ λ ∂μ
∂μ1 ∂μ2 1 2
(35)
⎡ ⎤ ⎡ ⎤
∂ ŷn+1 ∂ ŷn+1
− ∂u n − ∂u n
⎢ . ⎥ ⎢ ⎥
⎢ . ⎥  ⎢ .. ⎥  
⎢ . ⎥ ∂u n ∂u n ⎢ . ⎥ bg −gn −bg
= ⎢ ∂ ŷ ⎥ ∂μ1 ∂μ2 = ⎢ ∂ ŷ ⎥ gn = J m J c,
⎢− n+K ⎥ ⎢− n+K ⎥ gn
⎣ ∂u n ⎦ ⎣ ∂u n ⎦
√ √
λ λ
√
and e = [ên+1 , . . . , ên+K λΔu n ]T . The adaptation mechanism for the controller
is illustrated in Fig. 6. If the current state vector and jth support vector are cn =
[u n−1 , . . . , u n−n u , yn , . . . , yn−n y +1 ]T and x j = [x j1 , . . . , x jn u , x jn u +1 , . . . , x jn u +n y ]T ,
respectively, K-step ahead future behaviour of the system can be approximated as
   
ŷn+k = f n ĉn+k + gn ĉn+k u n ,
 T
ĉn+k = u n , . . . , u n u n−1 , . . . , u n+k−n u ŷn+k , . . . , ŷn+1 yn , . . . , yn+k−n y , (36)
         
k n u −k k−1 n y +1−k

where
   
    −dn+k, j
f n ĉn+k = α j K dn+k, j + b f = α j ex p + bf,
2σ 2
j∈SV j∈SV
    (37)
    −dn+k, j
gn ĉn+k = θ j K dn+k, j + bg = θ j ex p + bg .
2σ 2
j∈SV j∈SV

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 871

dn+k, j is defined as the Euclidean distance between the state vector at time step (n + k)th
and jth support vector x j as:

 T    
dn+k, j = cn+k − x j cn+k − x j = DUn+k + DYn+k = An+k , (38)

where

 u)
min(k,n
 2
DUn+k = u n − x j,i
i=1

nu
 2
+ u n+min(k,n u )−i − x j,i δ (n u − k) ,
i = min(k,n u )+1
min(k−1,n y )
(39)
  2
DYn+k = ŷn+k−i − x j,n u +i
i=1
ny
  2  
+ yn+min(k−1,n y )+1−i − x j,n u +i δ ny + 1 − k ,
i = min(k−1,n y )+1

and δ(·) is unit step function. The system Jacobian is given as

∂ ŷn+k ∂ f n (cn+k ) ∂dn+k, j ∂gn (cn+k ) ∂dn+k, j

= + u n + gn (cn+k ) , (40)
∂u n ∂dn+k, j ∂u n ∂dn+k, j ∂u n

where  
∂ f n (cn+k ) −1  −dn+k, j
= α j ex p ,
∂dn+k, j 2σ 2 2σ 2
j∈SV
  (41)
∂gn (cn+k ) −1  −dn+k, j
= θ j ex p ,
∂dn+k, j 2σ 2 2σ 2
j∈SV

and  
∂dn+k, j ∂dn+k, j ∂ An+k ∂ DUn+k ∂ An+k ∂ DYn+k
= + ,
∂u n ∂ An+k ∂ DUn+k ∂u n ∂ DYn+k ∂u n
∂dn+k, j 1 1
= √ = ,
∂ An+k 2 An+k 2dn+k, j
∂ An+k ∂ An+k
= = 1,
∂ DUn+k ∂ DYn+k (42)
∂ DUn+k  u)
min(k,n
 
= 2 u n − x j,i ,
∂u n
i=1
min(k−1,n y )
∂ DYn+k    ∂ ŷn+k−i
= 2 ŷn+k−i − x j,n u +i δ(k − i).
∂u n ∂u n
i=1

123
872 K. Uçak, G. Öke Günel

4.4 Adaptive Predictive SVRNARMA-L2 Controller with Adaptive Filter

As mentioned in Sect. 2, NARMA-L2 controller generally produces a control signal with

more oscillation and chattering compared to other controllers such as MRAC, MPC, etc.
[6]. Since the oscillatory control input contains high frequency components, the unmod-
eled high frequency dynamics of the system can be excited. In order to reduce chattering
and oscillation, the control signal can be filtered [6]. The general form of the filter is
given as:
dnu −m
u fn m=0 qm n z
H (z) = = dde , (43)
u cn 1 + k=1 vkn z −k

where vk , k ∈ {1, . . . , dde } and qm , m ∈ {0, 1, . . . , dnu } are “feed-backward” and “feed-
forward” coefficients of the filter, and dde and dnu indicate the degree of the “feed-backward”
(denominator) and “feed-forward” (numerator) parts of the filter, respectively. The adaptation
mechanism for the controller is illustrated in Fig. 7. The filtered control signal is computed
as

dde 
dnu
u fn = − vkn u f n−k + qm n u cn−m . (44)
k=1 m=0

The tracking performance of the proposed controller depends on the parameters (μ j (·), j ∈
{1, 2}) used to compute output of submodels as well as the filter parameters (qm , vk ). A
predictive structure which takes into account the K-step ahead future behaviour of system
has been proposed to adjust both μ j (·), j ∈ {1, 2}, vk , k ∈ {1, . . . , dde } and qm , m ∈
{0, 1, . . . , dnu } parameters effectively. The objective function given in (33) is optimized
via Levenberg–Marquardt algorithm to tune parameters. The adjustable parameters can be
adapted using Levenberg–Marquardt algorithm [49] as follows:

Fig. 7 Adaptive SVRNARMA-L2 controller based on online SVR with adaptive filter

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 873

⎡ ⎤ ⎡ ⎤
μ1new μ1old
⎢ μ2new ⎥ ⎢ μ ⎥
⎢ ⎥ ⎢ 2old ⎥
⎢ q0new ⎥ ⎢ q ⎥
⎢ ⎥ ⎢ 0old ⎥
⎢ .. ⎥ ⎢ .. ⎥⎢
⎢ . ⎥ ⎢ . ⎥ ⎥ T
−1 T
⎢ ⎥=⎢ + J J + ηI
⎢qd ⎥ ⎢qdnu ⎥ ⎥
J e, (45)
⎢ nu new ⎥
⎢ old
⎥
⎢ v1 ⎥ v1old ⎥
⎢ new ⎥ ⎢
⎢ . ⎥ ⎢ . ⎥
⎣ .. ⎦ ⎢ ⎣ .. ⎦
⎥
vddenew vddeold

where
⎡ ∂en+1 ∂en+1 ∂en+1 ∂en+1 ∂en+1 ∂en+1
⎤
∂μ1 ∂μ2 ∂q0 ··· ∂qdnu ∂v1 ··· ∂vdde
⎢ ⎥
⎢ .. .. .. .. .. .. .. .. ⎥
⎢ . . . . . . . . ⎥
J =⎢
⎢ ∂en+K ∂en+K ∂en+K ∂en+K ∂en+K ∂en+K
⎥
⎥
⎢ ∂μ ∂μ ∂q ··· ∂q ∂v ··· ∂v ⎥
⎣ √ 1 √ 2 √ 0 √ dnu √ 1 √ dde ⎦
∂ λΔu f n ∂ λΔu f n ∂ λΔu fn ∂ λΔu f n ∂ λΔu f n ∂ λΔu f n
∂μ1 ∂μ2 ∂q0 ··· ∂qdnu ∂v1 ··· ∂vdde
⎡ ⎤
− ∂∂μ ŷn+1
− ∂∂μŷn+1 ∂ ŷn+1
∂q0 · · · ∂∂qŷn+1 − ∂ ∂vŷn+1
· · · − ∂∂vŷn+1
⎢ ..
1
..
2
..
dnu
.. ..
1
..
dde
⎥
⎢ .. .. ⎥
⎢ . . . . . . . . ⎥
=⎢ ⎥
⎢ − ∂ ŷn+K − ∂ ŷn+K ∂ ŷn+K · · · ∂ ŷn+K ∂ ŷn+K
− ∂v1 · · · − ∂vd ∂ ŷn+K ⎥
⎣ ∂μ1 ∂μ2 ∂q0 ∂qdnu ⎦
√ ∂Δu fn √ ∂Δu fn √ √ √ √ de
λ ∂μ1 λ ∂μ2 λu n n · · · λu n n−dnu λu fn−1 · · · λu f n−dde
⎡ ⎤
− ∂∂uŷn+1
⎢ fn ⎥
⎢ .. ⎥  
⎢ . ⎥ ∂u fn ∂u nn ∂u fn ∂u nn
=⎢ ⎥ ∂u n ∂μ1 ∂u n ∂μ2 u n n · · · u n n−dnu −u fn−1 · · · −u fn−dde
⎢− ∂ ŷn+K ⎥ n n
⎣ ∂u fn ⎦
√
λ
⎡ ⎤
∂ ŷn+1
− ∂u f
⎢ n ⎥
⎢ .. ⎥  
⎢ . ⎥ b −g −b
=⎢ ⎥ q0n ggn n q0n gng u n n · · · u n−dnu −u f n−1 · · · −u f n−dde = J m J c ,
⎢− ∂ ŷn+K ⎥
⎣ ∂u fn ⎦
√
λ
√ (46)
and e = [ên+1 , . . . , ên+K λΔu f n ]T . The Jacobian matrix J m is derived via SVRNARMA-L2
model of the system as given in (40–42). An adaptive first order filter is adequate to suc-
cessfully filter high frequency components and expeditiously force the system to track the
reference signal for the proposed controller mechanism. Therefore, a first order adaptive low
pass filter with an adjustable parameter is deployed to filter high frequency components of
the control input as given by:
u fn q 0n
H (z) = = , (47)
u cn 1 + (1 − q0n )z −1

where 1 − q0 , q0 are “feed-backward” and “feed-forward” coefficients of the filter, respec-

tively. The adjustable parameters can be adapted using Levenberg–Marquardt algorithm [49]
as follows:

123
874 K. Uçak, G. Öke Günel

⎡ ⎤ ⎡ ⎤
μ1new μ1old

⎣μ2new ⎦ = ⎣μ2old ⎦ + J T J + η I −1 J T e, (48)

q0new q0old
where ⎡ ⎤
∂en+1 ∂ ŷn+1 ∂en+1 ∂ ŷn+1 ∂en+1 ∂ ŷn+1 ∂u f n
∂ ŷn+1 μ1 ∂ ŷn+1 μ2 ∂ ŷn+1 ∂u f n ∂q0
⎢ .. .. .. ⎥
⎢ ⎥
⎢ . . . ⎥
J = ⎢ ∂e ∂u f n ⎥
⎢ n+K ∂ ŷn+K ∂en+K ∂ ŷn+K ∂en+K ∂ ŷn+K ⎥
⎣ ∂ ŷn+K
√ μ1 ∂ ŷn+K
√ μ2 ∂ ŷn+K√ ∂u f n ∂q0 ⎦
∂ λΔu f n ∂ λΔu f n ∂ λΔu f n
∂μ1 ∂μ2 ∂q0
⎡ ∂ ŷn+1 ∂u f ∂u n n ∂ ŷn+1 ∂u f n ∂u n n ∂ ŷn+1
⎤
− ∂u f ∂u nn μ1 − ∂u f n ∂u n n μ2 − ∂u f n (u n n − u f n−1 )
⎢ n
..
n
.. .. ⎥
⎢ ⎥
⎢ . . . ⎥
= ⎢ ∂ ŷ ⎥
⎢− n+K ∂u fn ∂u nn − ∂ ŷn+K ∂u fn ∂u nn − ∂ ŷn+K (u n − u f )⎥
⎣ ∂u fn ∂u nn μ1 ∂u f n ∂u n n μ2 ∂u f n n n−1 ⎦
√ ∂u fn ∂u n √ ∂u fn ∂u n √
λ ∂u n μ1n λ ∂u n μ2n λ(u n n − u f n−1 )
⎡ ⎤ n n (49)
− ∂ ŷn+1
⎢ ∂u fn ⎥
⎢ .. ⎥  
⎢ ⎥ ∂u ∂u ∂u ∂u 
= ⎢ . ⎥ ∂u nfn ∂μn1n ∂u nfn ∂μn2n u n n − u fn−1
⎢− ∂ ŷn+K ⎥ n n
⎣ ∂u fn ⎦
√
λ
⎡ ⎤
∂ ŷn+1
− ∂u f
⎢ n ⎥
⎢ .. ⎥  
⎢ . ⎥ b −g −b 
=⎢ ⎥ q0n ggn n q0n gng u n n − u fn−1 = J m J c .
⎢− ∂ ŷn+K ⎥
⎣ ∂u fn ⎦
√
λ

4.5 Control Procedure

The control procedure for proposed SVRNARMA-L2 controller is as follows:

Step 1: Initialize Controller, Filter and Model parameters.
– Controller parameters: μ j , j ∈ {1, 2}
– Filter parameters: qm , m ∈ {0, 1, . . . , dnu }, vk , k ∈ {1, . . . , dde }
– SVRNARX model parameters: β = bβ = 0, C = 1000, ε = 10−3
Step 2: Train SVRNARX model using (u fn−1 , yn ) data pair
– Set index to n.
– Constitute xn−1 = [u n−2 , . . . , u n−n u −1 yn−1 , . . . , yn−n y ]
– Predict ŷn via SVRNARX model
– Calculate em n = yn − ŷn
If |em n | > ε
Train system model where em n = yn − ŷn
else
Continue with system model parameters obtained at previous step
Step 3: Convert SVRNARX model to SVRNARMA-L2 model
– Compute parameters of Fn−1 (α, b f ) and G n−1 (θ , bg ) via parameters of SVRNARX
(β, bβ ) as in (29)

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 875

– Constitute xn = [u n−1 , . . . , u n−n u yn , . . . , yn−n y +1 ]

– Compute fˆn− and ĝn− via Fn−1 (·) and G n−1 (·)
Step 4: Calculate the control input produced by SVRNARMA-L2 controller (u cn )
– Calculate the control signal (u cn ) using the predictions of the submodels and control law
in (32)
Step 5: Filter the high frequency components of control signal (u fn )
– Compute the filtered control signal via (44)
Step 6: Apply the filtered control signal u f n to the system to compute yn+1
– Training data pair of SVRNARX model for the next step is obtained (u fn , yn+1 )
Step 7: Calculate the Jacobian
– Apply u fn K-times to SVRNARX model in order to constitute Jacobian matrix in (46 or
49)
Step 8: Learning Step for Controller and Filter
– Update Controller parameters (μ j , j ∈ {1, 2}) and filter parameters (qm , m ∈
{0, 1, . . . , dnu }, vk , k ∈ {1, . . . , dde }) using (45 or 48)
Step 9: Increment n = n + 1 and go to step 2.

5 Simulation Results

5.1 The Bioreactor System

The performance of the proposed SVRNARMA-L2 controller is examined on a bioreactor

system. Bioreactor is a tank containing water and cells (e.g., yeast or bacteria) which consume
nutrients (substrate) and produce product (both desired and undesired) and more cells [53].
Since the uncontrolled equations of this system are highly nonlinear and exhibit limit cycles
[53], it is a very difficult system to control and it has frequently been employed as a nonlinear
benchmark problem to examine the performances of proposed adaptive controllers [17,31,
54,55]. The differential equations describing the cell concentration [c1 (t)] and amount of
nutrients [c2 (t)] are given as follows:
c2 (t)
ċ1 (t) = −c1 (t)u(t) + c1 (t) (1 − c2 (t)) e γ (t) ,
c2 (t) 1 + β(t) (50)
ċ2 (t) = −c2 (t)u(t) + c1 (t) (1 − c2 (t)) e γ (t) ,
1 + β(t) − c2 (t)
where c1 (t) is the controlled output of the system (y(t) = c1 (t)), u(t) is the flow rate as
the control signal, γ (t) is nutrient inhibition parameter, β(t) is grow rate parameter [17,
31,53,55]. In the simulations, the max–min limitations for the magnitude of the control
signal are given as u min = 0 and u max = 2; and duration of control signal is choosen as
τmin = τmax = 0.5 s. Since we have assumed that the dynamics of the system is not known,
online SVRNARX has been utilized to identify the unknown dynamics using the input–output
data pairs. Then an equivalent SVRNARMA-L2 model is obtained via proposed method given
in (25–29, 34). Thus, SVRNARMA-L2 submodels can be computed to design the control law
given in (32). The input feature vectors for SVRNARX and SVRNARMA-L2 models are chosen
as Mc = [u n−1 , . . . , u n−n u , yn , . . . , yn−n y +1 ]T where n u = 3 and n y = 3 indicate the
number of the past inputs and outputs, respectively.

123
876 K. Uçak, G. Öke Günel

5.2 SVR Design Parameters

Identification methods based on ε-SVR possess several design parameters which have direct
influence on modeling performance. ε which is the most significant parameter in identification
is the maximum permissible modeling error in training phase of ε-SVR. ε-SVR cannot
accurately capture the dynamics of the system when ε is chosen too large and the controller
may induce unacceptable control performance. If it is chosen too small, identification can
result in too many support vectors than required to model successfully. Therefore ε, the
maximum tolerable error for system model, must be chosen so that an acceptable control
quality [17] is obtained while the number of support vector for SVRNARX model is kept at
a reasonable level. Accordingly, maximum tolerable error is assigned as ε = 10−3 . Also,
the maximum number of the support vector is limited to a predefined value to decrease
computational load in our simulations. As can be seen in (13), C is the boundary of the
parameter search box in which Lagrange multipliers are searched. Since the control algorithm
gets slower as the number of support vectors increases, C is fixed at a large value (1000) in
order to reduce the number of the support vectors as low as possible. The mapping of the input
data to feature space and the nonlinearity in SVR are provided via kernel functions, hence
selection of kernel parameters is also crucial in regression performance. In our simulations,
−||x−y||
we utilized an exponential kernel function with K (x, y) = e 2σ 2 . The value of σ effects
the separability of the data in feature space. Since the kernel matrix stores similarities of the
input data, very small σ maps similar data in input space to distant locations in feature space,
the number of the support vectors for system model increases to capture the system dynamics.
On the other hand, large σ values maps different data in input space to close locations in
feature space and results in loss of the nonlinearity of the kernel that leads to inaccurate
identification of system dynamics. Therefore, the similarities of the data in input space must
be reflected in the feature space as correctly as possible. We have also observed that very
small σ values cause fluctuations in control signal and system output since the number of
the support vectors increases to capture the system dynamics. On the other hand, when σ
is chosen too large, the nonlinearity of the kernel is withered and the control signal is too
slow, so the system cannot be forced to track reference signal accurately. In simulations for
bioreactor system SVRNARX and SVRNARMA-L2 have been designed with σ = 0.75. The
closed-loop performance of the system is observed for different values of prediction horizon
(K ) and penalty parameter (λ) using the performance criteria in (51) and the optimal values
are determined with the grid search algorithm.
∞
1 
2
2
Jcomp = rn+1 − yn+1 + λ u n − u n−1 . (51)
2
n=1

As can be seen from Fig. 8, as λ increases, the fluctuation of control signal is restricted,
therefore, the tracking performance deteriorates. The performance is positively affected as
K increases. Thus, the values of λ and K for bioreactor are determined as K = 9, λ =
0.005, Jcompmin = 0.013235 via the minimum of (51) given in Fig. 8.

5.3 Simulation Results

The tracking performance of controller for the case when no noise is applied to the system
and all parameters are fully known is depicted in Fig. 9. The controller and filter parameters
are illustrated in Fig. 10. In order to assess the performance of the controller with respect
to measurement noise, 30 dB measurement noise is added to the output of the system. The

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 877

0.4

0.3
Jcomp
0.2

0.1

0.05
0.04
8
0.03
6
0.02 4
0.01 2
λ K

Fig. 8 Tracking performance surface with respect to penalty term (λ) and prediction horizon (K ) for bioreactor

(a)
0.15
r(t)
y(t)
0.1

0.05

0 100 200 300 400 500

(b)
0.6
u(t)
0.5
0.4
0.3
0.2
0.1
0 100 200 300 400 500

(c)
12
10
8
6
4
2 #svm(t)
0
0 100 200 300 400 500
Time[sec]

Fig. 9 System output (a), control signal (b), and number of support vectors (c) with adaptive filter for the
case with no noise or parametric uncertainty

123
878 K. Uçak, G. Öke Günel

(a)
0
μ1(t)
-0.02

-0.04

-0.06

-0.08
0 100 200 300 400 500

(b)
0
-0.05
-0.1
-0.15
-0.2
μ2(t)
-0.25
0 100 200 300 400 500

(c)
0.2
0.15
0.1
0.05 q0(t)
0
0 100 200 300 400 500
Time[sec]

Fig. 10 Adaptive controller (a, b) and filter parameters (c) for the case with no noise or parametric uncertainty

tracking performance is illustrated in Fig. 11. Adaptation of the controller and filter parameters
are given in Fig. 12. The robustness of the controller is inspected in terms of parametric
uncertainty, γ (t) is considered as the time-varying system parameter, where its nominal
value is γnom (t) = 0.48. The value of γ (t) is allowed to slowly alter in neighborhood of
its nominal value as γ (t) = 0.48 + 0.06 sin(0.012πt). It can be interpreted from Fig. 13
that the controller can successfully manage to reject the uncertainty in system parameters.
Controller and filter parameters are illustrated in Fig. 14. The maximum number of the support
vectors are limited to 10, 25 and 35 for noiseless, noisy, and parametric uncertainty cases,
respectively.

5.4 Comparison of the Results with SVM-Based PID Controller

A performance comparison of SVRNARMA-L2 controller proposed in this paper and the SVM-
based PID controller implemented by Iplikci [17] has been carried out for cases with no noise,
with measurement noise added to the system and with parametric uncertainty. The tuning
mechanism described in [17] adjusts the parameters of a PID controller via Levenberg–
Marquardt algorithm by using the approximated K-step ahead future system behaviour. The
controller consists of five components: classical incremental PID controller, SVRNARX model
of the system, line search block, control signal correction block and controller parameter
tuner. The Jacobian matrix utilized in parameter tuner block is coalescence of correction

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 879

(a)
0.15
r(t)
y(t)
0.1

0.05

0 100 200 300 400 500

(b)
0.6
u(t)
0.5
0.4
0.3
0.2
0.1
0 100 200 300 400 500

(c)
30

20

10
#svm(t)

0
0 100 200 300 400 500
Time[sec]

Fig. 11 System output (a), control signal (b), and number of support vectors (c) with adaptive filter for the
case with measurement noise

block which includes K-step ahead system Jacobian and a gradient vector which involves first
order derivative of control signal with respect to controller parameters. K-step ahead future
Jacobian of the system is approximated via SVRNARX model of the system. Control signal
correction block produces a correction term for the control signal since the adjusted controller
parameters may not be adequate to force the system output to track the desired reference. The
control signal correction term is derived via Taylor approximation of the control signal term
which includes correction term. Golden section method has been employed in line search
block to determine optimum learning rate for correction term. The simulation results for
the SVRNARMA-L2 controller proposed in this paper are compared with results obtained for
the structure given in [17]. Same prediction horizon (K ) and penalty term (λ) have been
employed in both SVRNARMA-L2 controller and SVM-based PID controller for comparison.
The tracking performances of the controllers for cases with no noise, with measurement noise
and with parametric uncertainty are illustrated in Figs. 15, 16 and 17 for bioreactor system,
respectively. Comparisons of the performance index given by (51) for tracking performance
of both controllers are illustrated in Fig. 18 which shows that SVRNARMA-L2 controller has
better tracking performance than SVM-based PID controller for all cases except the case
with parametric uncertainty. The applicability of the proposed controller in real time is as
significant as its tracking performance. For this purpose, computation times of each operation

123
880 K. Uçak, G. Öke Günel

(a)
0
μ1(t)
-0.01

-0.02

-0.03

0 100 200 300 400 500

(b)
0
-0.02
-0.04
-0.06
-0.08 μ2(t)
-0.1
0 100 200 300 400 500

(c)
0.15

0.1

0.05
q0(t)

0
0 100 200 300 400 500
Time[sec]

Fig. 12 Adaptive controller (a, b) and filter parameters (c) for the case with measurement noise

in the control algorithm have been registered for each case during every sampling period and
the average response times of the operations have been listed in Table 1 for all conditions.
In our simulations, a PC with 2.2 GHz core i7 CPU and 8 GB RAM has been operated for
the implementation of the control algorithm without optimizing codes. As can be seen from
the Table 1, the average response times of the controller accrete for measurement noise and
parametric uncertainty cases since the controller struggles to tolerate the uncertainty in the
system. In simulations, the sampling time is chosen as 100 ms and the total average response
times of the proposed controller are less than 50 ms as given in Table 1. Consequently, it
can be expressed that SVRNARMA-L2 controller can be deployed in real-time applications
conveniently. Moreover, total response times can be minimised and the controller can be
easily employed by implementing the algorithm on a lower level programming environment
such as C/C++ utilizing effective hardwares (e.g., FPGA) in real time application.

6 Conclusion

In this paper, a novel SVRNARMA-L2 controller is proposed where online SVR is used to
estimate SVRNARMA-L2 submodels. Firstly, a SVRNARX model is constructed to approximate
the system, then the parameters of the SVRNARMA-L2 model are derived via SVRNARX model
of the system. Consequently, SVRNARMA-L2 controller is designed using the SVRNARMA-L2

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 881

(a)
0.15
r(t)
y(t)
0.1

0.05

0 100 200 300 400 500

(b)
0.8
u(t)
0.6

0.4

0.2

0 100 200 300 400 500

(c)
0.55
γ(t)

0.5

0.45

0 100 200 300 400 500

(d)
40

30

20

10 #svm(t)

0
0 100 200 300 400 500
Time[sec]

Fig. 13 System output (a), control signal (b), uncertain system parameter (c), and number of support vectors
(d) with adaptive filter for the case with parametric uncertainty

submodels. The main contribution of the paper is that it justifies the use of an online SVRNARX
model of the system directly to derive online SVRNARMA-L2 controller as opposed to existing
works in technical literature where SVRNARX models are generally utilized for adjusting
the parameters of the traditional controllers by approximating system Jacobians. Another
novelty of this work is to use SVR to estimate system model in NARMA-L2 controller
design instead of ANN. In technical literature, NARMA-L2 controllers have previously
been designed using ANNs to approximate system model. However, ANNs are trained with
backpropagation algorithm and have the drawback that they may get stuck at local minima
which will lead to models valid only locally. The main strength of SVR in modelling is that
it minimizes a convex goal function and achieves global minimum. The performance of the

123
882 K. Uçak, G. Öke Günel

(a)
0
μ1(t)
-0.02

-0.04

-0.06

0 100 200 300 400 500

(b)
0
-0.05
-0.1
-0.15
-0.2
μ2(t)
-0.25
0 100 200 300 400 500

(c)
0.25
0.2
0.15
0.1
0.05 q0(t)

0
0 100 200 300 400 500
Time[sec]

Fig. 14 Adaptive controller (a, b) and filter parameters (c) for the case with parametric uncertainty

SVR NARMA-L2
(a)
0.1 r(t)
0.08 y(t)

0.06
0.04
0.02

0 100 200 300 400 500

(b) SVM-based PID

0.1 r(t)
0.08 y(t)

0.06
0.04
0.02

0 100 200 300 400 500

Time[sec]

Fig. 15 Tracking performance of SVRNARMA-L2 controller (a) and SVM-based PID controller (b) for the
case with no noise or parametric uncertainty

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 883

SVRNARMA-L2
(a)
0.1 r(t)
0.08 y(t)

0.06
0.04
0.02

0 100 200 300 400 500

(b) SVM-based PID

0.1 r(t)
0.08 y(t)

0.06
0.04
0.02

0 100 200 300 400 500

Time[sec]

Fig. 16 Tracking performance of SVRNARMA-L2 controller (a) and SVM-based PID controller (b) for the
case with measurement noise

SVRNARMA-L2
(a)
0.1 r(t)
0.08 y(t)

0.06
0.04
0.02

0 100 200 300 400 500

(b) SVM-based PID

0.1 r(t)
0.08 y(t)

0.06
0.04
0.02

0 100 200 300 400 500

(c)
0.55
γ(t)

0.5

0.45

0 100 200 300 400 500

Time[sec]

Fig. 17 Tracking performance of SVRNARMA-L2 controller (a) and SVM-based PID controller (b) for the
case with parametric uncertainty (c)

123
884 K. Uçak, G. Öke Günel

0.025
SVRNARMA-L2
SVM-based PID
0.02

0.015

Jcomp

0.01

0.005

0
Bioreactor Bioreactor Bioreactor
Noiseless Noisy Disturbance

Fig. 18 Tracking performance comparison of the controller with respect to the defined performance index
(51)

Table 1 Computation times (in ms)

Bioreactor
Operations Noiseless Noisy Disturbance

SVRNARX model training 6.8487 24.42 18.8102

Conversion to SVRNARMA-L2 model 0.93768 0.95917 0.95705
K-step Jacobian 6.9377 8.8229 16.0968
Parameter tuning (LM) 0.029813 0.031787 0.037375
Miscellaneous tasks 2.9456 2.9348 3.0298
Total time 17.6995 37.1686 38.9312

proposed controller is evaluated by simulations done on a bioreactor benchmark system.

Additionally, the success of the controller is compared with a SVM-based PID controller.
The robustness of the controller against system parameter uncertainty and measurement noise
have also been examined. The results indicate that the proposed controller is quite successful
in attaining low tracking error, suppressing measurement noise and parametric uncertainties.
In future works, it is planned to extend the derivation of SVRNARMA-L2 model via SVRNARX
to develop new SVR type adaptive controller design methods.

References
1. Narendra KS, Mukhopadhyay S (1997) Adaptive control using neural networks and approximate models.
IEEE Trans Neural Netw 8(3):475–485. doi:10.1109/72.572089
2. Majstorovic M, Nikolic I, Radovic J, Kvascev G (2008) Neural network control approach for a two-tank
system. In: 9th Symposium on neural network applications in electrical engineering (NEUREL 2008),
Belgrade, Serbia
3. Pedro JO, Nyandoro OTC, John S (2009) Neural network based feedback linearisation slip control of an
anti-lock braking system. In: Asian control conference (ASCC 2009), Hong Kong, China
4. De Jesus O, Pukrittayakamee A, Hagan MT (2001) A comparison of neural network control algorithms.
In: International joint conference on neural networks (IJCNN’01), Washington, DC

123
A Novel Adaptive SVRNARMA-L2 Controller for Nonlinear Systems 885

5. Pukrittayakamee A, De Jesus O, Hagan MT (2002) Smoothing the control action for NARMA-L2 con-
trollers. In: 45th Midwest symposium on circuits and systems (MWSCAS 2002), Tulsa, OK
6. Hagan MT, Demuth HB, De Jesus O (2002) An introduction to the use of neural networks in control
systems. Int J Robust Nonlinear Control 12(11):959–985. doi:10.1002/rnc.727
7. Wahyudi M, Mokri SS, Shafie AA (2008) Real time implementation of NARMA L2 feedback linearization
and smoothed NARMA L2 controls of a single link manipulator. In: International conference on computer
and communication engineering, Kuala Lumpur, Malaysia
8. Akbarimajd A, Kia S (2010) NARMA-L2 controller for 2-DoF underactuated planar manipulator. In:
International conference on control, automation, robotics and vision (ICARCV 2010), Singapore
9. Vesselenyi T, Dzitac S, Dzitac I, Manolescu MJ (2007) Fuzzy and neural controllers for a pneumatic
actuator. Int J Comput Commun Control 2(4):375–387
10. Efe MO, Kaynak O (2000) A comparative study of soft-computing methodologies in identification of
robotic manipulators. Robot Auton Syst 30(3):221–230
11. Efe MO, Kaynak O (1999) A comparative study of neural network structures in identification of nonlinear
systems. Mechatronics 9(3):287–300. doi:10.1016/S0957-4158(98)00047-6
12. Denai MA, Palis F, Zeghbib A (2004) ANFIS based modelling and control of non-linear systems: a
tutorial. In: IEEE international conference on systems, man and cybernetics
13. Gretton A, Doucet A, Herbrich R, Rayner PJW, Scholkopf B (2001) Support vector regression for black-
box system identification. In: IEEE workshop on statistical signal processing (SSP 2001), Singapore
14. Rong HN, Zhang GX, Zhang CF (2005) Application of support vector machines to nonlinear system iden-
tification. In: International symposium on autonomous decentralized systems (ISADS 2005), Chengdu,
China
15. Suykens JAK (2001) Nonlinear modelling and support vector machines. In: IEEE instrumentation and
measurement technology conference (IMTC 2001), Budapest, Hungary
16. Ucak K, Oke G (2011) Adaptive PID controller based on online LSSVR with kernel tuning. In: Inter-
national symposium on innovations in intelligent systems and applications (INISTA 2011), Istanbul,
Turkey
17. Iplikci S (2010) A comparative study on a novel model-based PID tuning and control mechanism for
nonlinear systems. Int J Robust Nonlinear Control 20(13):1483–1501. doi:10.1002/rnc.1524
18. Vapnik VN (1999) An overview of statistical learning theory. IEEE Trans Neural Netw 10(5):988–999.
doi:10.1109/72.788640
19. Wanfeng S, Shengdun Z, Yajing S (2008) Adaptive PID controller based on online LSSVM identification.
In: IEEE/ASME international conference on advanced intelligent mechatronics (AIM 2008), Xian, China
20. Zhao J, Li P, Wang Xs (2009) Intelligent PID controller design with adaptive criterion adjustment via
least squares support vector machine. In: 21st Chinese control and decision conference (CCDC 2009),
Guilin, China
21. Yuan XF, Wang YN, Wu LH (2008a) Composite feedforward–feedback controller for generator excitation
system. Nonlinear Dyn 54(4):355–364. doi:10.1007/s11071-008-9334-6
22. Takao K, Yamamoto T, Hinamoto T (2006) A design of PID controllers with a switching structure by
a support vector machine. In: 2006 IEEE international joint conference on neural network (IJCNN),
Vancouver, Canada
23. Liu XJ, Yi JQ, Zhao DB (2005) Adaptive inverse control system based on least squares support vector
machines. In: 2nd international symposium on neural networks (ISNN 2005), Chongqing, China
24. Wang H, Pi DY, Sun YX (2007) Online SVM regression algorithm-based adaptive inverse control. Neu-
rocomputing 70(4–6):952–959. doi:10.1016/j.neucom.2006.10.021
25. Yuan XF, Wang YN, Wu LH (2008b) Adaptive inverse control of excitation system with actuator uncer-
tainty. Neural Process Lett 27(2):125–136. doi:10.1007/s11063-007-9064-7
26. Zhao ZC, Liu ZY, Xia ZM, Zhang JG (2012) Internal model control based on LS-SVM for a class of
nonlinear process. In: International conference on solid state devices and materials science (SSDMS),
Macao, China
27. Zhong WM, Pi DY, Sun YX, Xu C, Chu SZ (2006) SVM based internal model control for nonlinear
systems. In: 3rd International symposium on neural networks (ISNN 2006), Chengdu, China
28. Datta A, Ochoa J (1996) Adaptive internal model control: design and stability analysis. Automatica
32(2):261–266. doi:10.1016/0005-1098(96)85557-9
29. Sun CY, Song JY (2007) An adaptive internal model control based on LS-SVM. In: International sym-
posium on neural networks (ISNN 2007), Nanjing, China
30. Wang YN, Yuan XF (2008) SVM approximate-based internal model control strategy. Acta Autom Sin
34(2):172–179. doi:10.3724/SP.J.1004.2008.00172
31. Iplikci S (2006a) Online trained support vector machines-based generalized predictive control of non-
linear systems. Int J Adapt Control Signal Process 20(10):599–621. doi:10.1002/acs.919

123
886 K. Uçak, G. Öke Günel

32. Iplikci S (2006b) Support vector machines-based generalized predictive control. Int J Robust Nonlinear
Control 16(17):843–862. doi:10.1002/rnc.1094
33. Camacho EF (1993) Constrained generalized predictive control. IEEE Trans Autom Control 38(2):327–
332. doi:10.1109/9.250485
34. Clarke DW, Mohtadi C (1989) Properties of generalized predictive control. Automatica 25(6):859–875.
doi:10.1016/0005-1098(89)90053-8
35. Camacho EF, Bordons C (1999) Model predictive control. Springer, New York
36. Clarke DW, Mohtadi C, Tuffs PS (1987) Generalized predictive control—Part I. The basic algorithm.
Automatica 23(2):137–148. doi:10.1016/0005-1098(87)90087-2
37. Zhiying D, Xianfang W (2008) Nonlinear generalized predictive control based on online SVR. In: 2nd
International symposium on intelligent information technology application, Shanghai, China
38. Shin J, Kim HJ, Park S, Kim Y (2010) Model predictive flight control using adaptive support vector
regression. Neurocomputing 73(4–6):1031–1037. doi:10.1016/j.neucom.2009.10.002
39. Abu-Rub H, Awwad A (2009) Artificial neural networks and fuzzy logic based control of AC motors. In:
IEEE international electric machines and drives conference (IEMDC 2009), Miami, FL
40. Norgaard M, Ravn O, Poulsen NK, Hansen LK (2000) Neural networks for modelling and control of
dynamic systems. Springer, London
41. Ng GW (1997) Application of neural networks to adaptive control of nonlinear systems. In: UMIST
control systems centre series, 4 (Book 4), Research Studies Press Ltd, Taunton, UK
42. Cortes C, Vapnik V (1995) Support-vector networks. Mach Learn 20(3):273–297. doi:10.1023/A:
1022627411411
43. Drucker H, Burges CJC, Kaufman L, Smola A, Vapnik V (1997) Support vector regression machines. In:
Annual conference on neural information processing systems (NIPS), Denver, CO
44. Vapnik V, Golowich SE, Smola A (1997) Support vector method for function approximation, regres-
sion estimation, and signal processing. In: Annual conference on neural information processing systems
(NIPS), Denver, CO
45. Smola AJ, Schölkopf B (2004) A tutorial on support vector regression. Stat Comput 14(3):199–222.
doi:10.1023/B:STCO.0000035301.49549.88
46. Cristianini N, Shawe-Taylor J (2000) An introduction to support vector machines and other kernel-based
learning methods. Cambridge University Press, Cambridge
47. Martin M (2002) On-line support vector machine regression. In: 13th European conference on machine
learning (ECML 2002), Helsinki, Finland
48. Ma J, Theiler J, Perkins S (2003) Accurate online support vector regression. Neural Comput 15(11):2683–
2703
49. Luenberger DG, Ye Y (2008) Linear and nonlinear programming, 3rd edn. Springer, New York
50. Griva I, Nash SG, Sofer A (2009) Linear and nonlinear optimization, 2nd edn. SIAM, Philadelphia
51. Nocedal J, Wright SJ (1999) Numerical optimization. Springer, New York
52. Cauwenberghs G, Poggio T (2001) Incremental and decremental support vector machine learning. In:
Annual neural information processing systems conference (NIPS), Denver, CO
53. Ungar LH (1990) Neural networks for control. In: Miller WT III, Sutton RS, Werbos PJ (eds) A bioreactor
benchmark for adaptive network based process control. MIT Press, Cambridge, pp 387–402
54. Efe MO, Abadoglu E, Kaynak O (1999) A novel analysis and design of a neural network assisted nonlin-
ear controller for a bioreactor. Int J Robust Nonlinear Control 9(11):799–815. doi:10.1002/(SICI)1099-
1239(199909)9:<11799::AID-RNC441>3.0.CO;2-U
55. Efe MO (2007) Discrete time fuzzy sliding mode control of a biochemical process. In: 9th WSEAS
international conference on automatic control, modeling and simulation (ACMOS’07), Istanbul, Turkey

123