HUNGARIAN JOURNAL 
OF INDUSTRIAL CHEMISTRY 

VESZPREM 
Vol. 29. pp. 129- 134 (2001) 

IDENTIFICATION OF NONLINEAR SYSTEMS USING GAUSSIAN MIXTURE 
OF LOCAL MODELS 

J. ABONYI, T. CHOVAN and F. SZEIFBRT 

(Department of Process Engineering, University ofVeszprem, P.O. Box 158, H-8201, HUNGARY) 

Received: October 8, 2001 

Identification of operating regime based models of nonlinear dynamic systems is addressed. The operating regimes and 
the parameters of the local linear models are identified directly and simultaneously based on the Expectation 
Maximization (EM) identification of Gaussian Mixture Model (GMM). The proposed technique is demonstrated by 
means of the identification of a neutralization reaction in a continuously stirred tank reactor. 

Keywords: Operating regime based model, expectation maximization, Takagi-Sugeno fuzzy model, nonlinear system, 
neutralization reaction 

Introduction 

The problem of a successful model based control 
application arises from difficulties in system modeling 
[1, 2]. This difficulty stems from lack of knowledge or 
understanding of the process to be controlled [3]. While 
it may not be possible to find process information that is 
universally applicable, .it would certainly be worthwhile 
to examine what types of process-knowledge would be 
most relevant for specific operating points of the 
process. This type of local understanding, in fact, will 
be a key to identifying reliable local models with a 
limited amount of data. The model that has a range of 
validity less than the operating regime of the process is 
called local model, as opposed to a global model that is 
valid in the full range of operation. Global modeling is a 
complicated task because of the need to describe the 
interactions between a large number of phenomena that 
appear globally. Local modeling, on the other hand, 
may be considerably simpler, because locally there may 
be a smaller number of phenomena that are relevant, 
and their interactions are simpler [4]. The modeling 
framework that is based on combining a number of local 
models. where each local model has a predefined 
operating region in which the local model is valid is 
called operating regime based model [5], where the 
local models are combined into a global model using an 
interpolation technique as it is illustrated in Fig. I. 

The main advantage of this framework is its 
transparency. Both the concept of operating regimes and 
the model structure are easy to understand. This is 

important, since the model structure can be interpreted 
in terms of operating regimes, but also quantitatively in 
terms of individual local models. 

The operating regime of the local models can be also 
represented by fuzzy sets [6]. This representation is 
appealing, since many systems change behaviors 
smoothly as a function of the operating point, and the 
soft transition between the regimes introduced by the 
fuzzy set representation captures this feature in an 
elegant fashion. 

Fuzzy modeling and identification proved to be 
effective tools for the approximation of uncertain 
nonlinear systems because of the ability to combine 
expert knowledge and measured data. Fuzzy models use 
if-then rules to describe the process through a collection 
of locally valid relationships. The antecedents (if-parts) 
of the rules divide the input space into several fuzzy 
subspaces, while the consequents (then-parts) describe 
the local behavior of the system in these fuzzy 
subspaces [7]. In this paper the local models are linear. 

The contribution of this paper is two-fold: 

• A new method for the identification of operating 
regime based models is proposed based on EM 
identification of Gaussian Mixtures Model. 

• Method to transform the obtained mt1del imo 
Takagi-Sugeno fuzzy model is presented. 

The paper is organized as follows. Section 2 pre~ents 
the structure of the operating regime based model along 
with the methods for its transformation into a fuzzy 
model. In Section 3, the identification algorithm of the 



130 

y(k) 

where the l/J;(xk) function describes the operating 

regime of the i -th local linear model defined by the 

ei =[a; ,bi r parameter vector. 
The operating regime of the local models can also be 

represented by fuzzy sets [6]. Hence, the entire global 
model Eq.(3) can be conveniently represented by 
Takagi-Sugeno fuzzy rules: 

R;: If xk is A;(xk) then Yk =a; xk +hi' i = l, ... ,c (4) 

u(k) where A; (xk) represents the multi variable membership 

function that describes the fuzzy set A; while a; and 

Fig. I Example for an operating regime based model. The 
operating region defined by the current input u(k) and output 

y(k) of the system is decomposed into four regimes. · 

model is proposed. An application example - the 
identification of a pH process is given in Section 3. 
Conclusions are given in Section 4. 

Operating Regime based Modeling of Dynamical 
System 

Nonlinear dynamic systems are often represented in the 
Nonlinear AutoRegressive with eXogenous input 
(NARX) model form, which establishes a nonlinear 
relationship between the past inputs and outputs and the 
predicted output: 

y(k +1) = f(y(k) •.•. ,y(k -ny),u(~-nd), .•. ,u(k -n")) (1) 

Here, n.,. and n
11 

denote the maximum lags considered 
for the output, and input terms, respectively, nd <nil is 
the discrete dead time, and f represents the mapping of 
the NARX model. 

The aim of this paper is to develop an algorithm for 
the identification of a transparent and easily 
interpretable model 

(2) 

based on some available training pattern y1 

and x4 = [tu ~ .... .xd Y t where k = 1~ ... ~ N denotes the 
index: of the k .. th data that can be used for 
identification. When Eq.(2) is used to represent a 
NARX modet the training pattern is formed to have a 
model that gives one-step-ahead prediction: 

Yt = y(k +ll, 
xt = [v(k) •... ~y(k -n~),u(k -nJ)• ...• u(k -n,.,).f. 

The operating regime based model of the system 
given by Eq.(2J is formulated as: 

r 

5,., = L9i(x,) (al xh +hi) (3) 
l=l 

bi are the parameters of the local linear model. 

Usually, the antecedent proposition "xk is A;(xk)" 

is expressed as a logical combination of simple 
propositions with univariate fuzzy sets defined for the 
individual components of xk , often in a conjunction 

form: 

R
1
• : If x 1 k is A1• 1 (x1 k) and ... and xn k is A; n (xn k) , , , , ' ' (5) 

then Yk =a; xk +h; 

In this case, the degree of fulfillment of a rule is 
calculated as the product of the degree of fulfillment of 
the fuzzy sets in the rule 

/3;(xk) = A;(xk) = Ii: Ai,j(xj,k) (6) 
j=l' 

The rules of the fuzzy model are aggregated using 
the fuzzy mean formula 

c 

L wi /3;(xk)(a;xk +b) 
Yk=~~~~---c----------- (7) 

LW; f3;(Xk) 
i=l 

where w; = [0,1] is the weight of the rule that represents 
the desired impact of the rule. The value of W; is often 

chosen by the designer of the fuzzy system based on his 
or her belief in the goodness and accuracy of the i -th 
rule. When such knowledge is not available W; is set as 

W; = 1, \;fj. 

As Eq.(7) and Eq.(3) show, fuzzy models are 
identical to operating regime based models as the 
operating region of the local linear models are defined 
by the normalized rule fulfillments: 

W; {3i(xk) 
r (8) 

l:wi /3;(XA) 
i=l 

To represent the Ai.i(xu) fuzzy set, in this paper 

Gaussian membership function is used 

(9) 



where v,. . represents the center and a,~ . the. variance of 
,J ,) 

the Gaussian function. 
The use of Gaussian membership functions allows 

the compact formulation of Eq.(6): 

,B;(x,) ~ A;(x,) ~exp( -~(x,- v;l (F;'T' (x,- vj)) (10) 

where v~ = [ v1,j' ... , v n) represents the center of the i -

th multivariate Gaussian and (F:Xf1 stands for the 
inverse of a diagonal matrix that contains the variances: 

[a~; 0 
XX 0 

2 

Fi = : 
(j2,i 

0 0 

1] (11) 
This compact formulation of the operating regime 

based model suggests that the model is not only 
equivalent to a TS fuzzy model but it is functionally 
identical to Generalized Radial Basis Function Network 
(GBFN) [8]. 

When there is a correlation among input variables of 

the model, F;xx is not a diagonal matrix as it was shown 

in Eq.(ll). In this case the decomposition of Ai(xk) to 

~./xj,k) fuzzy sets by 

A ( II
n ( ) IIn ( 1 (xj,k -zvi,jl) (12) 

i xk) = i=l ~.i xi,k = i=l exp 2 ai,j 

is not possible directly. The proposed method to solve 
this problem is based on the eigenvector projection [9] 
or transformed input-domain approach [10]. The 
approach is based on the calculation of the eigenvalues 

Ali and the eigenvectors t~ , of the ~xx matrix, where 
j = l, ... ,n. Using the eigenvalues and the eigenvectors 

the following fuzzy model that has no correlation in its 
transformed input space can be obtained: 

R; : If (t~)T xk is A;,1 and ... and (t~)r xk is A;,n (13
) 

then yk =a; xk + b1 

where the Gaussian membership functions are defined 
as 

h - _ ( i )T - _ ( i T x d -2 _ ( 11 )2 were x ... - t. x .. , v.,- t.) v. an 0';
1
·- A

1
· 

}." J " l,J J l , 

denote the transformed input variable, the cluster center 
and variance, respectively. 

The aim of the remaining part of the paper is to 
propose a new identification technique for the 
identification of the model presented above. 

131 

EM algorithm for Identification of Mixture of 
Gaussians Model 

The Expectation Maximization (EM) algorithm is an 
iterative algorithm for the computation of maximum 
likelihood parameter estimates when the observations 
can be viewed as incomplete data. The EM algorithm is 
widely used for parameter estimation of the mixture of 
models, in particular the mixture of Gaussians model 
[11]. 

The basics of EM are the following. Suppose we 
know the observed values of a random variable z and 
we wish to model the density of z using a model 
parameterized by rJ . 

Each observation consists of n + 1 measured 
variables, grouped into an n +!-dimensional column 

vector z k = [zl,k' ... 'Zn+l,k r . A set of N observations is 
denoted by Z and represented as a matrix: 

[ z,,, Zr,z z,,N l 
z = Zz,l Zz,z Zz,N (15) 

Zn:I,l Zn+l,2 Zn:!.N 

As the identification is performed on the available 
identification data, Z is divided into a regression data 
matrix X and a regression vector y 

(16) 

(17) 

In the pattern recognition terminology, the columns 
of Z called patterns or objects, the rows are called the 
features or attributes, and Z is called the pattern matrix. 

EM obtains parameter estimates fj which maximize 

the likelihood L(rJ) = p(z~) of the data. The EM 
assumes that this estimation is intractable and the values 
of a missing or hidden random variable h would make 
the problem more tractable. Let p(z,h~) denote the 
joint probability of z and h parameterized by 11 • It is 
assumed that z and h are such that maximizing the 
complete data likelihood Lr(rJ) = p(z,hjr]) is more 

tractable than maximizing L(q). However, the values of 

h are not known. The EM algorithm tackles this 
problem by iteratively generating a probability over ~he 
values h and estimating the parameters whtch 
maximize the expected value of L< VJ) with respect to 

b. 
The mixture of Gaussians model represents the 

p(z~) probability density function that b e\panded in 

a sum over the c clusters 



132 

c c 

p(z~) = L p(z,1Ji) = LP<zl1]j)p(1Ji) 
Next we compute the remaining parameters of the 

(18) cluster, the mean 
i=l i=l 

where 1J is the set of the parameters 1J = {1J;Ii = l, ... ,c} 
of the model and p(1J;) denotes the unconditioned 

cluster probabilities normalized to satisfy 

:L;=l p(1J;);:::: 1. 

The p(zj1Ji) distribution generated by the i -th 

cluster is represented by Gaussian like 

p(z~-) 1 exp(-.!.(z-v.l (F.f1 (z-v.)) (19) 
I ( 27C ff .JiFJ 2 l I l 

where 1J; represents the parameters of the i -th cluster, 

1]; ={v;,Fzli=l, ... ,c}. As each data point zk is 
generated by one and only one of the component 
Gaussians, the hidden random variable hk for the EM 

algorithm it is the label of the component Gaussian to 
which z k belongs. 

Based on this assumption, the EM algorithm for 
maximum likelihood parameter estimation is the 
following. The algorithm starts with an initial guess 

11'0' of the parameters and repeatedly applies the 
following two steps to generate successively better 
parameter estimates: 

Initialization: Initialize the means v; to randomly 
picked data points from Z and the covariance 
matrices Fj to unit matrices.' Set p(1Ji) = 1/ c for all 
i. 

Expectation (E) step: In the E-step we assume the 
current cluster parameter to be correct and evaluate 
the posterior probabilities that relate each data point in 
the conditional probability p (rJtjz). These posterior 

probabilities can be interpreted as the probability that 
a particular piece of data was generated by a particular 
cluster. By using Bayes Theorem, 

(. ·lz>= p(zJ1J;) p(1J;) p Y/, ' ( ) p z 

{21) 

Maximization (M) step: In the M-step we assume the 
current data distribution to be correct and find the 
parameters of the clusters that maximize the likelihood 
of the data. 

According to this~ the unconditional probabilities are 
calculated as 

f f 
P (1JiJz) 

v. = zp(z~-)dz= z--p(z)dz 
I I p(1Ji) 

(23) 

and in a similar way the cluster weighted covariance 
matrices 

N 

L(zk -v)(zk -v;l p(1Jijzk) 
k=l (24) 

In the above it has been shown that the Gaussian 
Mixture Model models the p(z) = p(x, y) joint density 

of the response variable yk and the regressors xk as a 

mixture of c multivariate n + 1 dimensional Gaussian 
functions. The conditional density p(yjx) is also a 

mixture of Gaussians model and the regression E[yjx] 

is 

I J I 
J y p(y,x)dy J y p(y,x)dy 

E[yx]= y p(yx)dy 
p(x) J p(x,y)dy (25) 

= i. [xr 1]8; ]p(~1J;) p(1J;) i. p (1J;Ix) [xr 1]8;] 
i=l p(x) i=l 

where ei denotes the parameter vector of the local 
model and p(1J;jx) denotes the probability that the i -th 

Gaussian component generated the regression vector x : 

p (Tl;) exp(-.!(x- v':)T (F:"T1 (x -·v':) l 
(21r )'2 .Jiif4 2 ' ' ' J (26) 

:t p(ij;) exp(-.!(x-v':l (F.""'f1 (x-v':)l 
i=l (21r )'2 .Jiif4 2 ' l ' J 

where the ~xx is obtained by the partitioning of the F; 
covariance matrix 

F -[· Ft FtY] 
i- ~]X F/Y 

(27) 

The optimal ei parameter vector of the local models 
can be obtained as: 

a1 = {FJxx t F;xy 
b; =v( -afv: 

or in a more compact form: 

(J = L.!"~lF.t.t \-1 ,Y -F'J'X{Fxx \-1 xJ 
t l.i'i \ i } ~ li i \ i J V; 

(28) 

(29) 



pH 

Fig.2 The considered continuous stirred tank reactor. 

526 

524 

522 

520 

j518 

516 

514 

512 
6.5 7.5 8.5 9 9.5 10 10.5 11 

pH(k) 

Fig.3 Orientation of the identified clusters. 

This is identical to the weighted, also called local, 
parameter estimation approach that does not estimate all 
parameters simultaneously, because the parameters of 
the local models are estimated separately using a set of 
local estimation criteria 

where xe denotes the extended regression matrix 
obtained by adding a unitary column to X, Xe = [X 1], 
and <Pi denote a diagonal matrix having membership 
degrees in its diagonal elements. 

[

f.li,l 0 

((). = 0 J.li,2 
l : : 

0 0 

(31) 

The weighted least-squares estimate of the 
consequent rule parameters is given by 

(32) 

Note that the resulted Gaussian Mixture of Local 
Models defined by Eq.(25) is identical to the operating 
regime based model given by Eq.(3) when the 

133 

6.5 

"o 200 400 600 800 ~000 
Time [time-step=0.2 min] 

Fig.4 Free-run test of the model on the validation data. 

weighting function is chosen as ¢i(x) = p(1J;jx). 
Furthermore, the model is also identical to TS fuzzy 
models, when the membership functions are identified 
as it was shown in Section 2. 

Application to the identification of a pH process 

The identification of the pH (the concentration of 
hydrogen ions) in a continuous stirred tank reactor 
(CSTR) is a well-known benchmark problem. The 
CSTR depicted in Fig.2. has two input streams: sodium 
hydroxide and acetic acid. The dynamic model for the 
pH in the tank is given in the Appendix. Po, collection 
of the data, a sampling interval of 0.2 min was used. 

As the process can be modeled as a first !Jrder 
dynamic system [12], the fuzzy model consisu, uf the 
following rules: 

Ri : If (pH(k), FNaOH (k)] is A (33) 
then pH (k + 1) = a1,1 ·pH (k) + a1•2 • F NuOil {k) + b, 

By using transformed input variables: 

Ri: If ~:.1 • pH(k)+t{,2 ·FNaou(k)j is As.t and 
&~.1 • pH(k)+t~.2 • FNuOil(k)] is A1.z (34} 

then pH(k + 1) = ai,t' pH(k)+a,,2 • F,.a011 (k) + b, 

Four local models were identified. The number of 
local models were determined by cross-validation. The 

operating region (the orientation of Ft matrices) of 
these models are shown in Fig.3. 

The identified model was tested with a validation 
data set. The model was used for one-step-ahead and 
simulation {infinite-step~ahead) prediction of the pH. 
The later experiments is depicted in Fig.4. 

The results were compared with the performance of 
models obtained by using Fuzzy Model ldentl tkation 
Toolbox [9]. As Table 1 shows, the propo'\ed method 
shows superior performance over this adv~mced tool 
developed for identification of nonline;1r dynamic 
systems. 



134 

Table 1 Comparison of Model Performances (mean squares of 
prediction errors). 

Method One-ste Simulation 
FMID [9] 
Proposed 

0.0241 
0.009 

Conclusions 

0.2835 
0.0956 

A new algorithm for the identification of nonlinear 
systems is proposed that is based on the Expectation 
Maximization identification of Gaussian Mixtures 
Model. A method to extract Takagi-Sugeno fuzzy 
models from Gaussian Mixtures Model is presented. 
The resulted fuzzy models are based on the transformed 
input-domain approach, which allows the effective 
partition of the input space and enables the 
interpretability of the model. The performance of the 
proposed modeling technique was demonstrated in the 
identification of the a pH process. 

Acknowledgement 

The financial support of the Hungarian Ministry of 
Culture and Education (FKFP-0023/2000, FKFP-
0073/200 1) and the Hungarian Science Foundation 
(T023157} is greatly acknowledged. Janos Abonyi is 
grateful for the financial support of the Janos Bolyai 
Research Fellowship of the Hungarian Academy of 
Science. Tibor Chovan is supported by the Hans Pape 
foundation. 

Model of the pH process 

A dynamic model of the pH in a tank can be obtained by 

considering the material balances on [Na+] and the 

total acetate [HAC+ AC-] and assuming that acid-base 
equilibrium and electroneutrality relationships hold 
[12]. 

F e.~tc • [HAClin - ( F HAc + FNuOH) ·[HAC+ AC-] = 
=Vd[HAC+AC-] 

dt 
Sodium ion balance: 

F,.,itOH ·[NaOH},,-(FnM: + fl..,uoH )·[Na+]= V d[~;+J 
HAC equilibrium: 

[AC"}·[H .. ] 
----=K 

{HAC] " 
Water equilibrium: 

[H .. l·[OH-1= K)! 
Electroneutra1ity: 

{Na.,.}+[H""J=[OH-}+[AC-1 

The pH can be calculated from the previous equations 
as 

Table 2 Parameters used in the simulations. 

Parameter 

FNaOH 

[NaOH!n 
[HACJn 

[Na+] 

Description 

Volume of the tank 

Flow rate of acetic acid 

Flow rate of NaOH 

Inlet concentration of NaOH 

Nominal 
Value 

1000 [1] 

81 [1/min] 

515 [1/min] 

0.05 [molll] 

Inlet concentration of acetic 0.32 [mol/1] 
acid 

Initial concentration of 0.0432 [molll] 
sodium in the CSTR 

r -] Initial concentration of acetate 0 0432 [mol!l] LHAC + AC in the CSTR . 
Ka Acid equilibrium constant 1.753 ·10-5 

Water equilibrium constant 

[H+]3 + [H+] 2 (Ka + [Na+]) + [H+]([Na+]Ka -

-[HAC+AC-]Ka -K.,.)-KwKa =0 

pH== log[H+] 

The parameters used in our simulations are taken 
from [12] and are given in Table 2. 

REFERENCES 

1. LEITCH R.: Comput. Control Eng. J., 1992, July, 
153-163 

2. COIT B. J., DURHAM R. G. and SULLIVAN G. R.: 
Comput. Chern. Eng., 1989, 13,973-984 

3. SJOBERG J., ZHANG Q., LJUNG L., BENVENISTE A., 
DEYWN B., GWRENNEC P-Y., HJALMARSSON H. 
and JUDITSKY A.: Automatica, 1995,31, 1691-1724 

4. MURRAY-SMITH R.: A Local Model Network 
Approach to Nonlinear Modelling, PhD Thesis, 
University of Strathclyde, Computer Science 
Department, 1994 

5. MURRAY-SMITH R. and JOHANSEN T. A. (Ed): 
Multiple Model Approaches to Nonlinear Modeling 
and Control, Taylor & Francis, London, UK, 1997 

6. BABUSKA R. and VERBRUGGEN H. B.: Fuzzy Set 
Methods for Local Modeling and Identification, in 
Multiple Model Approaches to Nonlinear Modeling 
and Control (Ed. MURRAY-SMITH R., JOHANSEN, T. 
A.), Taylor & Francis. London, UK, 75-100, 1997 

7. TAKAGI T. and SUGENO M.: IEEE Trans. Syst. Man 
Cybem.,l985, 15(1), 116-132 

8. HUNT K.J., HAAS R. and MURRAY-SMITH R.: IEEE 
Trans. Neural Netw., 1996, 7(3), 776-781 

9. BABUSKA R.: Fuzzy Modeling for Control, Kluwer 
Academic Publishers, Boston, MA, 1998 

10. KIM E., PARK M., KIM S. and PARK M.: IEEE 
Trans. Fuzzy Syst., 1998, 6, 596-604 

11. SCHONER B.: Probabilistic Cheracterization and 
Synthesis of Complex Driven Systems, Ph.D. 
Thesis~ Massachusetts Institute of Technology, 
2000 

12. BRAT N. V. and MCAVOY T. J.: Comput. Chern. 
Eng., 1992, 16,271-281 


	Page 130 
	Page 131 
	Page 132 
	Page 133 
	Page 134 
	Page 135