# Identification of Multimodel LPV Models with Asymmetric Gaussian Weighting Function

Authors: Jie You; Jiangang Lu; Yucai Zhu; Qinmin Yang; Jianhua Zhu; Jiangyin Huang; Youxian Sun
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

This paper is concerned with the identification of linear parameter varying (LPV) systems by utilizing a multimodel structure. To improve the approximation capability of the LPV model, asymmetric Gaussian weighting functions are introduced and compared with commonly used symmetric Gaussian functions. By this mean, locations of operating points can be selected freely. It has been demonstrated through simulations with a high purity distillation column that the identified models provide more satisfactory approximation. Moreover, an experiment is performed on real HVAC (heating, ventilation, and air-conditioning) to further validate the effectiveness of the proposed approach.


---

## Body

## 1. Introduction

In nonlinear system identification, numerous black box modeling structures have been developed in pieces of literature. Among them, nonlinear AR(MA)X and neural-network models are often used by researchers, due to their mature theoretical results. However, the complexity in those structures and calculation cost limit their applicability especially in process control applications. Meanwhile, block-oriented nonlinear models such as Hammerstein models and Wiener models consisting of linear time invariant (LTI) dynamics and static (memoryless) nonlinearities have been studied widely. Although they are comparatively simpler, only the model nonlinearity in static gains is integrated, which cannot ensure the accuracy and efficiency of the model, particularly for nonlinear plants operating over a large number of different operating points.Recently, linear parameter varying (LPV) model identification has attracted great attention from both academia and practitioners [1]. Many significant breakthroughs have taken place in the intervening years to mature the underlying theory. According to [2], the existing LPV identification approaches can be categorized in terms of their LPV model structures. Subspace approaches are studied by Verdult and Verhaegen [3, 4] and Felici et al. [5]; orthonormal basis related functions are used by Tóth et al. [6]; transfer function LPV models are discussed by Bamieh and Giarré [7], Previdi and Lovera [8], Wei [9], Butcher et al. [10], and Laurain et al. [1]. On the other hand, the practicality of the LPV approaches has caught the attention of practitioners outside of academia, and interests in the LPV approach are evidenced by the applications of such methods to aerospace systems including high performance aircraft, missiles, and turbofan engines [11].In studying input-output LPV methods, most available references are parameter interpolation based, which consider the parameters of the transfer function as nonlinear functions of the scheduling variable. Because complex nonlinear functions appear in the denominator of the transfer function, model stability cannot be guaranteed, which may result in numerical problems during model identification [12]. Besides, the input excitation signal for this representation results in too much upset, which can be costly or even unrealistic in practice [13]. To circumvent these difficulties, Zhu and Xu [12] proposed a multimodel LPV model based on blended linear models. The identification method is relatively simple and the stability of this kind of LPV models has been proved in Huang et al. [2].Proper weighting functions are required to combine all the local linear models into a global LPV model in the multimodel LPV structure. Several options are available in the literature, namely, linear weight function [14], polynomial function [12], cubic spline function [15], and Gaussian weight function [13]. However, due to large number of parameters to be estimated and unconstrained shapes for polynomial or cubic spline function, they will easily lead to an ill-conditioned and may not be suitable for identification of complex nonlinear processes [2]. Meanwhile, the linear weight function is not accurate enough to capture the full dynamic behaviors of a nonlinear process. In the realm of multimodel structures [16] and fuzzy sets [17], Gaussian weighting functions have been widely adopted, which have relative small number of parameters and naturally better functions than linear functions. However, a disadvantage of Gaussian weighting functions is that the operating points should have an equal distance with respect to a scheduling variable, while other functions do not have this restriction, which brings large inconveniences and limits their capabilities for usage in actual applications.Therefore, this paper aims at improving the performance in identifying multimodel LPV models by adopting the asymmetric Gaussian weighing function, such that multilocal models can be smoothly interpolated to approximate the global dynamical behaviors of the process. This means that the locations of the operating points can be freely selected. The parameters of the proposed LPV models are estimated by using Levenberg-Marquardt method. The LPV model with two scheduling variables using asymmetric Gaussian weighing is also discussed. Simulation study is employed to demonstrate the efficiency of the proposed LPV model identification scheme. Further, the experiment is conducted on an HVAC system in our lab. It should be noted that there are some references available in the literature [18] about the modeling of the HVAC system. However, the identification with multimodel LPV structure is very rare. The experimental results also show that the asymmetric Gaussian weighing function can obtain improved results for multi model LPV identification in real industrial processes.The rest of the paper is organized as follows. A brief revisit of the LPV model identification is given in Section2. The identification of multimodel LPV approach with asymmetric Gaussian is developed in Section 3. In Sections 4 and 5, simulation example and experiment studies are given to demonstrate the effectiveness of the proposed approach. Section 6 gives the conclusion.


## 2. LPV Model Identification

### 2.1. System Description of LPV Model

Given a multiinput single-output (MISO) LPV system, denote them inputs as {u1(t),…,um(t)} at time t and output as y(t). Denote w(t) as the scheduling variable which is a measured process variable from the process or can be calculated from measurable process variables. Throughout this paper, it is assumed that w(t)∈[wminwmax], where wmin and wmax are the low and high limits of w(t). The data generating the input-output LPV system can be described by the following equation [1]:
(1)y(t)=∑i=1mGi(q,w)ui(t)+v(t),
where
(2)Gi(q,w)=Bi(q,w)Ai(q,w),Bi(q,w)=[b1i(w)q-1+···+bni(w)q-n]q-di,i=1,…,m,Ai(q,w)=1+a1i(w)q-1+⋯+ani(w)q-n,
where Gi(q,w) is a LPV transfer function from the ith input to the output. Bi(q,w) and Ai(q,w) are polynomials of q-1, which denotes unit delay operator. di is the delay from the ith input to the output, v(t) is a stationary stochastic process with zero mean and bounded variance, and n is the order of the model.Moreover, polynomial method is a commonly used parameterization method to represent the LPV model; that is, the parametersbji(w) and aji(w) are replaced by polynomial functions of w(t) as follows:
(3)bji(w)=γ1i,j+γ2i,jw+⋯+γnβi,jwnβ-1,j=1,…,n,aji(w)=γ1i,j+γ2i,jw+⋯+γnαi,jwnα-1,j=1,…,n,
where nα and nβ are the order of the polynomial functions.Therefore, (1)–(3) formulate the input-output LPV model, which is called parameter-interpolation input-output LPV model by Zhao et al. [19].


### 2.2. Revisit of Multimodel LPV Model Identification

Zhu and Xu [12] pointed out that the LPV model in (1)–(3) is not easy to identify due to its complex structure and thus proposed a simpler LPV model structure, which is in a form of blended local linear models. The basic principle of the approach is to first identify several local linear models at fixed operating points. The global model is then obtained by interpolation via certain weighting functions, which is called multimodel LPV model by Huang et al. in [2] or model interpolation based input-output LPV (MI-IO-LPV) by Zhao et al. in [19]. This model also falls into the multimodel structure or operating regimebased model studied by many authors; see, for example, Johansen and Foss [20] and Boukhris et al. [16].To identify the multimodel LPV model, local linear models at each fixed operating point should be determined first. Assume thatl operating points are firstly set as
(4)wmin≤w1<w2<⋯<wl≤wmax.Obviously, the choice ofl is a tradeoff of computing cost and model accuracy. Then, the corresponding local linear models can be represented as
(5)y^k(t)=∑i=1mG^ik(q)ui(t),k=1,2,…,l,
where y^k(t) denotes the output of the linear model at w(t)=wk. The transfer function G^ik(q), i=1,2,…,m of the kth local linear model in (5) can be expressed in the following form:
(6)G^ik(q)=Bik(q)Aik(q)=[b1i,kq-1+⋯+bni,kq-n]q-di,k1+a1i,kq-1+⋯+ani,kq-n.The parameters to be estimated for each local model can be written as(7)θ^ki=[a1i,k⋯ani,kb1i,k⋯bni,kdi,k](n+n+1)×1.Parameters to be estimated for alll local linear models can be denotes as
(8)ΘL=[θ^11⋯θ^1mθ^21⋯θ^2m⋯θ^l1⋯θ^lm].Then, instead of building the complex LPV model in form of (1)–(3), the multimodel global LPV model is adopted to represent the process along the operating trajectory as
(9)y^(t)=∑k=1lηk(w(t))y^k(t)+v(t),
where ηk(w(t)), k=1,2,…,l, are weight functions of corresponding local linear models which are essentially static functions of the scheduling variable w(t). v(t) is the white noise as defined in (1).For multiinput multioutput (MIMO) processes, the procedure can be repeated for each output. As a matter of fact, this representation is more suitable for modeling processes which operate on a few fixed operating points with quick transition between neighboring operating points. Essentially, with the use of local linear models and model interpolating philosophy, the identification method is relatively simple and the stability of this kind of LPV models can be guaranteed readily [2]. Therefore, the multimodel LPV model considerably simplifies the task of the input-output LPV method identification. Several pieces of literatures [2, 13, 19] have verified that the multimodel LPV in (5)–(9) can be a good approximation of the real process along its operating-trajectory, while the commonly used parameter-interpolation LPV models in (1)–(3) fail to obtain an acceptable result in a case study [2]. To combine all the local linear models into a global multimodel LPV model, proper weighting functions are required for interpolation, which post a large impact on the accuracy of the global model.Several common weighting functions are available in the literature, for example, linear weight function [14], polynomial function [12], cubic spline function [15], and Gaussian weight function [13]. By comparison, Gaussian weighting functions can achieve the best performance in terms of good input-output fit and accurate step responses [2]. However, a disadvantage of the Gaussian weighting functions is that the operating points for local linear models should have identical distance with respect to a scheduling variable, which involves large inconvenience in practical usage. To overcome this drawback, developing the multimodel LPV approach with a more suitable weighting function motivates us to conduct this research. Hence, this paper can be considered as an important extension of the work by Huang et al. [2].Before proceeding, identification test methods for multimodel LPV model utilized here will be briefly discussed. Tests are performed at some fixed operating points plus transition tests, which only cover the trajectory of the process operation. Through the whole operating range, small test signals should be added to the inputs or set points to persistently excite the system.


## 2.1. System Description of LPV Model

Given a multiinput single-output (MISO) LPV system, denote them inputs as {u1(t),…,um(t)} at time t and output as y(t). Denote w(t) as the scheduling variable which is a measured process variable from the process or can be calculated from measurable process variables. Throughout this paper, it is assumed that w(t)∈[wminwmax], where wmin and wmax are the low and high limits of w(t). The data generating the input-output LPV system can be described by the following equation [1]:
(1)y(t)=∑i=1mGi(q,w)ui(t)+v(t),
where
(2)Gi(q,w)=Bi(q,w)Ai(q,w),Bi(q,w)=[b1i(w)q-1+···+bni(w)q-n]q-di,i=1,…,m,Ai(q,w)=1+a1i(w)q-1+⋯+ani(w)q-n,
where Gi(q,w) is a LPV transfer function from the ith input to the output. Bi(q,w) and Ai(q,w) are polynomials of q-1, which denotes unit delay operator. di is the delay from the ith input to the output, v(t) is a stationary stochastic process with zero mean and bounded variance, and n is the order of the model.Moreover, polynomial method is a commonly used parameterization method to represent the LPV model; that is, the parametersbji(w) and aji(w) are replaced by polynomial functions of w(t) as follows:
(3)bji(w)=γ1i,j+γ2i,jw+⋯+γnβi,jwnβ-1,j=1,…,n,aji(w)=γ1i,j+γ2i,jw+⋯+γnαi,jwnα-1,j=1,…,n,
where nα and nβ are the order of the polynomial functions.Therefore, (1)–(3) formulate the input-output LPV model, which is called parameter-interpolation input-output LPV model by Zhao et al. [19].


## 2.2. Revisit of Multimodel LPV Model Identification

Zhu and Xu [12] pointed out that the LPV model in (1)–(3) is not easy to identify due to its complex structure and thus proposed a simpler LPV model structure, which is in a form of blended local linear models. The basic principle of the approach is to first identify several local linear models at fixed operating points. The global model is then obtained by interpolation via certain weighting functions, which is called multimodel LPV model by Huang et al. in [2] or model interpolation based input-output LPV (MI-IO-LPV) by Zhao et al. in [19]. This model also falls into the multimodel structure or operating regimebased model studied by many authors; see, for example, Johansen and Foss [20] and Boukhris et al. [16].To identify the multimodel LPV model, local linear models at each fixed operating point should be determined first. Assume thatl operating points are firstly set as
(4)wmin≤w1<w2<⋯<wl≤wmax.Obviously, the choice ofl is a tradeoff of computing cost and model accuracy. Then, the corresponding local linear models can be represented as
(5)y^k(t)=∑i=1mG^ik(q)ui(t),k=1,2,…,l,
where y^k(t) denotes the output of the linear model at w(t)=wk. The transfer function G^ik(q), i=1,2,…,m of the kth local linear model in (5) can be expressed in the following form:
(6)G^ik(q)=Bik(q)Aik(q)=[b1i,kq-1+⋯+bni,kq-n]q-di,k1+a1i,kq-1+⋯+ani,kq-n.The parameters to be estimated for each local model can be written as(7)θ^ki=[a1i,k⋯ani,kb1i,k⋯bni,kdi,k](n+n+1)×1.Parameters to be estimated for alll local linear models can be denotes as
(8)ΘL=[θ^11⋯θ^1mθ^21⋯θ^2m⋯θ^l1⋯θ^lm].Then, instead of building the complex LPV model in form of (1)–(3), the multimodel global LPV model is adopted to represent the process along the operating trajectory as
(9)y^(t)=∑k=1lηk(w(t))y^k(t)+v(t),
where ηk(w(t)), k=1,2,…,l, are weight functions of corresponding local linear models which are essentially static functions of the scheduling variable w(t). v(t) is the white noise as defined in (1).For multiinput multioutput (MIMO) processes, the procedure can be repeated for each output. As a matter of fact, this representation is more suitable for modeling processes which operate on a few fixed operating points with quick transition between neighboring operating points. Essentially, with the use of local linear models and model interpolating philosophy, the identification method is relatively simple and the stability of this kind of LPV models can be guaranteed readily [2]. Therefore, the multimodel LPV model considerably simplifies the task of the input-output LPV method identification. Several pieces of literatures [2, 13, 19] have verified that the multimodel LPV in (5)–(9) can be a good approximation of the real process along its operating-trajectory, while the commonly used parameter-interpolation LPV models in (1)–(3) fail to obtain an acceptable result in a case study [2]. To combine all the local linear models into a global multimodel LPV model, proper weighting functions are required for interpolation, which post a large impact on the accuracy of the global model.Several common weighting functions are available in the literature, for example, linear weight function [14], polynomial function [12], cubic spline function [15], and Gaussian weight function [13]. By comparison, Gaussian weighting functions can achieve the best performance in terms of good input-output fit and accurate step responses [2]. However, a disadvantage of the Gaussian weighting functions is that the operating points for local linear models should have identical distance with respect to a scheduling variable, which involves large inconvenience in practical usage. To overcome this drawback, developing the multimodel LPV approach with a more suitable weighting function motivates us to conduct this research. Hence, this paper can be considered as an important extension of the work by Huang et al. [2].Before proceeding, identification test methods for multimodel LPV model utilized here will be briefly discussed. Tests are performed at some fixed operating points plus transition tests, which only cover the trajectory of the process operation. Through the whole operating range, small test signals should be added to the inputs or set points to persistently excite the system.


## 3. Multimodel LPV Model with Asymmetric Gaussian Weights

A preferable choice of determining the model weights is Gaussian function. The representation of Gaussian weighting functions can be written as(10)ηk(w(t))=α~k(w(t))∑j=1lα~j(w(t)),k=1,2,…,l,
where
(11)α~k(w(t))=exp[-12(w(t)-wkσk)2]
and σk represents the width coefficients of the kth local linear model.


### 3.1. Identification of Multimodel LPV Model with Asymmetric Gaussian Weights

To overcome the drawback of the Gaussian weights, asymmetric Gaussian functions are introduced here. The weights can be written as(12)ηk(w(t))=αk(w(t))∑j=1lαj(w(t)),k=1,2,…,l,
where
(13)αk(w(t))={exp[-12(w(t)-wkσk1)2],w(t)<wk,exp[-12(w(t)-wkσk2)2],w(t)≥wk=exp[-12×(-sign(w(t)-wk)+12)(w(t)-wkσk1)2000000000×(w(t)-wkσk1)2]+exp[-12×(sign(w(t)-wk)+12)(w(t)-wkσk1)20000000000×(w(t)-wkσk2)2]-1,
where σk1 and σk2 called left and right variances, respectively, are the width coefficients of the scheduling variable; sign(·) is the sign function.The asymmetric Gaussian function in (12) are normalized is in the range of zero to one which prevents negative values. Essentially, this maintains the constrained shape and smooth properties of the Gaussian function. Since asymmetric Gaussian function has different left and right variances, the operating points for local linear models can have different distances with respect to a scheduling variable. While for traditional symmetric Gaussian functions, the range in w(t) should be divided with equal intervals.Remark 1.
For a given number of operating points and a given operating range, the selection of the operating points is based on trial-and-errors. By combing the prior knowledge of empirical data, if we can locate operate points to the place where the process exhibits strong nonlinearity and parameter varying properties, a more accuracy model can be obtained. In the meantime, even though the model accuracy can be improved by selecting more operating points, it inevitably means more local linear models and more test costs. Moreover, in many industrial processes, the process cannot run around some “favorable” operating points for test purpose due to safety and economic consideration. Therefore, the selection of operation points with unequal distances provides an effective and convenient way to improve the accuracy of the LPV model.Obviously, only two parameters need to be estimated for each weighting function. The parameter vector to be estimated for alll weighting functions can be defined as
(14)ΘG=[σ11σ21⋯σl1σ12σ22⋯σl2].Then, all the parameters to be estimated can now be written in a compact form:(15)Θ=[ΘGΘL]T.To estimate theΘ in (15), nonlinear optimization algorithm is desired which can minimize the following output error loss function:
(16)f(Θ)=∑t=1N[y(t)-∑k=1lηk(w(t))y^k(t)]2.


### 3.2. Parameter Estimation Scheme

Several nonlinear numerical optimization algorithms are available for this purpose, such as Gauss-Newton algorithm, steepest descent method, expectation-maximization algorithms and Levenberg-Marquardt method [2, 21, 22]. In this work, the Levenberg-Marquardt algorithm is used. The iterative optimization flow is as follows.Step 1.
InitializeΘ as Θ0=[ΘG0ΘL0]T. Denote the convergence criterion ∥Θ^(p+1)-Θ^(p)∥2≤ψ and set the iteration count p=0.Step 2.
Calculate simulated prediction output as(17)y^(t,Θ^(p))=∑k=1lηk(w(t)∣Θ^(p))∑i=1mBik(q)Aik(q)ui(t).Step 3.
Calculate the gradient matrix ofy^(t,Θ^(p))(18)J(t,Θ^(p))=∂∂Θy^(t,Θ^(p))=[∂∂σk1y^(t,Θ^(p))∂∂σk2y^(t,Θ^(p))∂∂ajiy^(t,Θ^(p))∂∂bjiy^(t,Θ^(p))∂∂diy^(t,Θ^(p))](2l+2mn+m)×1,0000000(i=1,2,…,m;j=1,2,…,n;k=1,2,…,l),
where each one of the entries can be readily obtained by using basic calculus.
Thereafter, negative gradient vector is represented as(19)d(p)=-[1N∑t=1N[J(t∣Θ^(p))·J(t∣Θ^(p))T]+μ(p)I]-1·[-1N∑t=1N[J(t∣Θ^(p))·e(t∣Θ^(p))]],
where e(t∣Θ^(p))=y(t)-y^(t,Θ^(p)), I is a unit matrix with proper dimension, and the scalar μ(p) is the Levenberg-Marquardt parameter.Step 4.
Calculateα(p) in minα(p)>0f(Θ^(p)+α(p)×d(p)) and update the parameters by Θ^(p+1)=Θ^(p)+α(p)×d(p), p=p+1.Step 5.
If∥Θ^(p+1)-Θ^(p)∥2≤ψ is satisfied, stop iteration and Θ^(p) is the optimized parameters in (16); otherwise, go to Step 2.It is noticed that in terms of the nonlinear numerical optimization problems, the initial values have a great influence on the estimation results. Hence, it is desirable to have a few appropriate initial guesses in order to secure the optimality of the algorithm. In order to circumvent the problem of local minimum, in the proposed method, the proper initial values of all the parametersΘL in the local linear models can be obtained by utilizing several linear identification methods [23] using the process input-output data collected around each corresponding operating point firstly. Afterwards, genetic algorithm (GA) is employed in attempt to search for the appropriate initial guess for the parameters ΘG for each weighting function associated with relevant local models.Following the procedure in the literature [24, 25], a pool of candidates (chromosomes) is randomly selected first and the fitness of each individual in the population is evaluated through certain fitness calculation function. Afterwards, different genetic operators including selection, reproduction, and mutation are applied to maintain the genetic diversity. The algorithm is implemented iteratively until certain termination criterion is met. Out of all the procedures adopted in GA, fitness measurement of the candidates is of particular importance as it directly influences the solution of the algorithm. For system parameters estimation, the mean square errors (MSE) of the predictions of the identified model against the real process data are calculated based on the following equation:
(20)MSEΘG=∑t=1N(y(t)-y^(t))2N,
where N in (20) represents the number of data points.


### 3.3. Extension to Two Scheduling Variables

Up to now, only one scheduling variable is assumed to explain the idea. If one scheduling variable may not be enough to describe the complex nonlinear dynamics of the industrial process, LPV models with two scheduling variables can be easily introduced by using the similar way in Huang et al. [2].Assume that that there are linear models atl1×l2 operating points:
(21)w1(t):wmin1≤w11<w21<⋯<wl11≤wmax1;w2(t):wmin2≤w12<w22<⋯<wl22≤wmax2.Then, the LPV model with two scheduling variables can be obtained by interpolating the local linear models. Consider(22)y^(t)=∑k=1l1∑h=1l2[∑i=1mG^ik,h(q)ui(t)+v(t)ηk,h(w1(t),w2(t))000000×∑i=1mG^ik,h(q)ui(t)+v(t)],
where v(t) is a Gaussian distributed noise with zero mean and variance o2 and ηk,h(w1(t),w2(t)) are asymmetric Gaussian weighting functions, which can be rewritten as
(23)ηk,h(w1(t),w2(t))=αk,h(w1(t),w2(t))∑r=1l1∑z=1l2αr,z(w1(t),w2(t)),0000000000000000(k=1,2,…,l1;h=1,2,…,l2),αk,h(w1(t),w2(t))={exp[-12((w1(t)-wk1σk,h1,1)2+(w2(t)-wh1σk,h2,1)2)],w1(t)<wk1,w2(t)<wh1,exp[-12((w1(t)-wk1σk,h1,1)2+(w2(t)-wh1σk,h2,2)2)],w1(t)<wk1,w2(t)≥wh1,exp[-12((w1(t)-wk1σk,h1,2)2+(w2(t)-wh1σk,h2,1)2)],w1(t)≥wk1,w2(t)<wh1,exp[-12((w1(t)-wk1σk,h1,2)2+(w2(t)-wh1σk,h2,2)2)],w1(t)≥wk1,w2(t)≥wh1.Note that more operating point tests and transition tests may be needed for multimodel LPV with two scheduling variables.


## 3.1. Identification of Multimodel LPV Model with Asymmetric Gaussian Weights

To overcome the drawback of the Gaussian weights, asymmetric Gaussian functions are introduced here. The weights can be written as(12)ηk(w(t))=αk(w(t))∑j=1lαj(w(t)),k=1,2,…,l,
where
(13)αk(w(t))={exp[-12(w(t)-wkσk1)2],w(t)<wk,exp[-12(w(t)-wkσk2)2],w(t)≥wk=exp[-12×(-sign(w(t)-wk)+12)(w(t)-wkσk1)2000000000×(w(t)-wkσk1)2]+exp[-12×(sign(w(t)-wk)+12)(w(t)-wkσk1)20000000000×(w(t)-wkσk2)2]-1,
where σk1 and σk2 called left and right variances, respectively, are the width coefficients of the scheduling variable; sign(·) is the sign function.The asymmetric Gaussian function in (12) are normalized is in the range of zero to one which prevents negative values. Essentially, this maintains the constrained shape and smooth properties of the Gaussian function. Since asymmetric Gaussian function has different left and right variances, the operating points for local linear models can have different distances with respect to a scheduling variable. While for traditional symmetric Gaussian functions, the range in w(t) should be divided with equal intervals.Remark 1.
For a given number of operating points and a given operating range, the selection of the operating points is based on trial-and-errors. By combing the prior knowledge of empirical data, if we can locate operate points to the place where the process exhibits strong nonlinearity and parameter varying properties, a more accuracy model can be obtained. In the meantime, even though the model accuracy can be improved by selecting more operating points, it inevitably means more local linear models and more test costs. Moreover, in many industrial processes, the process cannot run around some “favorable” operating points for test purpose due to safety and economic consideration. Therefore, the selection of operation points with unequal distances provides an effective and convenient way to improve the accuracy of the LPV model.Obviously, only two parameters need to be estimated for each weighting function. The parameter vector to be estimated for alll weighting functions can be defined as
(14)ΘG=[σ11σ21⋯σl1σ12σ22⋯σl2].Then, all the parameters to be estimated can now be written in a compact form:(15)Θ=[ΘGΘL]T.To estimate theΘ in (15), nonlinear optimization algorithm is desired which can minimize the following output error loss function:
(16)f(Θ)=∑t=1N[y(t)-∑k=1lηk(w(t))y^k(t)]2.


## 3.2. Parameter Estimation Scheme

Several nonlinear numerical optimization algorithms are available for this purpose, such as Gauss-Newton algorithm, steepest descent method, expectation-maximization algorithms and Levenberg-Marquardt method [2, 21, 22]. In this work, the Levenberg-Marquardt algorithm is used. The iterative optimization flow is as follows.Step 1.
InitializeΘ as Θ0=[ΘG0ΘL0]T. Denote the convergence criterion ∥Θ^(p+1)-Θ^(p)∥2≤ψ and set the iteration count p=0.Step 2.
Calculate simulated prediction output as(17)y^(t,Θ^(p))=∑k=1lηk(w(t)∣Θ^(p))∑i=1mBik(q)Aik(q)ui(t).Step 3.
Calculate the gradient matrix ofy^(t,Θ^(p))(18)J(t,Θ^(p))=∂∂Θy^(t,Θ^(p))=[∂∂σk1y^(t,Θ^(p))∂∂σk2y^(t,Θ^(p))∂∂ajiy^(t,Θ^(p))∂∂bjiy^(t,Θ^(p))∂∂diy^(t,Θ^(p))](2l+2mn+m)×1,0000000(i=1,2,…,m;j=1,2,…,n;k=1,2,…,l),
where each one of the entries can be readily obtained by using basic calculus.
Thereafter, negative gradient vector is represented as(19)d(p)=-[1N∑t=1N[J(t∣Θ^(p))·J(t∣Θ^(p))T]+μ(p)I]-1·[-1N∑t=1N[J(t∣Θ^(p))·e(t∣Θ^(p))]],
where e(t∣Θ^(p))=y(t)-y^(t,Θ^(p)), I is a unit matrix with proper dimension, and the scalar μ(p) is the Levenberg-Marquardt parameter.Step 4.
Calculateα(p) in minα(p)>0f(Θ^(p)+α(p)×d(p)) and update the parameters by Θ^(p+1)=Θ^(p)+α(p)×d(p), p=p+1.Step 5.
If∥Θ^(p+1)-Θ^(p)∥2≤ψ is satisfied, stop iteration and Θ^(p) is the optimized parameters in (16); otherwise, go to Step 2.It is noticed that in terms of the nonlinear numerical optimization problems, the initial values have a great influence on the estimation results. Hence, it is desirable to have a few appropriate initial guesses in order to secure the optimality of the algorithm. In order to circumvent the problem of local minimum, in the proposed method, the proper initial values of all the parametersΘL in the local linear models can be obtained by utilizing several linear identification methods [23] using the process input-output data collected around each corresponding operating point firstly. Afterwards, genetic algorithm (GA) is employed in attempt to search for the appropriate initial guess for the parameters ΘG for each weighting function associated with relevant local models.Following the procedure in the literature [24, 25], a pool of candidates (chromosomes) is randomly selected first and the fitness of each individual in the population is evaluated through certain fitness calculation function. Afterwards, different genetic operators including selection, reproduction, and mutation are applied to maintain the genetic diversity. The algorithm is implemented iteratively until certain termination criterion is met. Out of all the procedures adopted in GA, fitness measurement of the candidates is of particular importance as it directly influences the solution of the algorithm. For system parameters estimation, the mean square errors (MSE) of the predictions of the identified model against the real process data are calculated based on the following equation:
(20)MSEΘG=∑t=1N(y(t)-y^(t))2N,
where N in (20) represents the number of data points.


## 3.3. Extension to Two Scheduling Variables

Up to now, only one scheduling variable is assumed to explain the idea. If one scheduling variable may not be enough to describe the complex nonlinear dynamics of the industrial process, LPV models with two scheduling variables can be easily introduced by using the similar way in Huang et al. [2].Assume that that there are linear models atl1×l2 operating points:
(21)w1(t):wmin1≤w11<w21<⋯<wl11≤wmax1;w2(t):wmin2≤w12<w22<⋯<wl22≤wmax2.Then, the LPV model with two scheduling variables can be obtained by interpolating the local linear models. Consider(22)y^(t)=∑k=1l1∑h=1l2[∑i=1mG^ik,h(q)ui(t)+v(t)ηk,h(w1(t),w2(t))000000×∑i=1mG^ik,h(q)ui(t)+v(t)],
where v(t) is a Gaussian distributed noise with zero mean and variance o2 and ηk,h(w1(t),w2(t)) are asymmetric Gaussian weighting functions, which can be rewritten as
(23)ηk,h(w1(t),w2(t))=αk,h(w1(t),w2(t))∑r=1l1∑z=1l2αr,z(w1(t),w2(t)),0000000000000000(k=1,2,…,l1;h=1,2,…,l2),αk,h(w1(t),w2(t))={exp[-12((w1(t)-wk1σk,h1,1)2+(w2(t)-wh1σk,h2,1)2)],w1(t)<wk1,w2(t)<wh1,exp[-12((w1(t)-wk1σk,h1,1)2+(w2(t)-wh1σk,h2,2)2)],w1(t)<wk1,w2(t)≥wh1,exp[-12((w1(t)-wk1σk,h1,2)2+(w2(t)-wh1σk,h2,1)2)],w1(t)≥wk1,w2(t)<wh1,exp[-12((w1(t)-wk1σk,h1,2)2+(w2(t)-wh1σk,h2,2)2)],w1(t)≥wk1,w2(t)≥wh1.Note that more operating point tests and transition tests may be needed for multimodel LPV with two scheduling variables.


## 4. Simulation Example

As one of the most important unit operations in the chemical industry, high purity distillation has been widely investigated and accepted as a highly challenging process for nonlinear identification and control [26, 27]. When the column is operated over a relatively wide operation region, it shows a significant nonlinear behavior. In this section, the performance of the multimodel LPV models with asymmetric Gaussian is verified by identifying the high purity distillation column. The test data are generated by simulating a rigorous physical model of the distillation column, which is operated in the LV-configuration [28]. The inputs or manipulated variables (MVs) are reflux (LT) and boil-up rate (VB); the outputs or controlled variables (CVs) are the top product composition (yD) and bottom product composition (xB). A schematic representation of the high purity distillation column is shown in Figure 1.Figure 1
Diagram of the distillation column underLV configuration.Because the top product composition(yD) is the most important quality index of the process and has the significant influence over the process dynamics, it is chosen as the scheduling variable, which reflects the operating condition of the process. The choice of the numbers of operating points is a tradeoff of computing cost and model accuracy. Without loss of generality, 3 typical operating points with unequal distances are selected as
(24)yD:w(t)=0.95,w(t)=0.976,w(t)=0.99.Thereafter, to identify the global LPV model throughout the whole operating trajectory, more identification tests are preformed along the scheduling variableyD with length of 39,000 min. The sampling time is 1 min. The test signal is chosen as a generalized binary noise (GBN) signal with average switch time of 100 min. A filtered white noise sequence with 1% of variance is added into the output yD and xB to represent measurement noises. Second order output error local linear model is identified at each fixed operating point. Then, the global LPV model with two inputs and two outputs is achieved by utilizing all the training data. To validate the identified model, additional data are generated throughout the operating trajectory that is different from the one for the training data generation. In other words, the feasibility of the identified LPV model is verified by data attained from the process operating on different operating trajectory.Figure2 shows the input excitation dataset of the validation test. As can be seen from Figures 3 and 4, the identification results for the multimodel LPV model using Gaussian weights and asymmetric Gaussian weights are compared with the measured outputs. It has to be noted that when using Gaussian weights the operating points w(t) should be set with equal intervals. Therefore, the three operating points for LPV model with Gaussian weights are predetermined as 0.95, 0.97, and 0.99.Input excitation data.
(a)
(b)Outputs of the simulated and of the LPV model with Gaussian weights.
(a)
(b)Outputs of the simulated and of the LPV model with asymmetric Gaussian weights.
(a)
(b)Table1 compares the output errors of linear models and LPV models using Gaussian weights and asymmetric Gaussian weights. The output errors are defined as
(25)Error%=[(var(y(t)-y^(t))var(y(t)))]×100,
where y(t) is the output of real process and y^(t) is the predicted output of the LPV model.Table 1
Model output error comparison.
Models
y
D
x
B
Linear model
0.5668
%
22.5239
%
LPV model
Gaussian weights
0.2395
%
3.2305
%
LPV model
Asymmetric Gaussian weights
0.1811
%
1.0421
%As pointed out in Huang et al. [2], the simulation output error alone may not be sufficient to assure good model quality for the whole range and can even be misleading in nonlinear system identification. Thence, one good supplemental way to verify model quality is to check its step responses, which is often utilized in industrial applications. Step responds of the models and the actual process at other two operating points: 0.98 and 0.986 are shown in Figures 5 and 6, respectively.Figure 5
Step responds atyD=0.98.Figure 6
Step responds atyD=0.986.One can see that the proposed multimodel LPV with asymmetric Gaussian weights can provide more satisfactory approximation to the dynamics of the distillation column operating on a large number of different operating points.


## 5. Experimental Study

To further illustrate the effectiveness of the proposed approach, nonlinear process identification experiment is designed and conducted on an HVAC (heating, ventilation, and air-conditioning) system in our lab. A photograph of the system with R22 as the refrigerant and a simplified schematic diagram are displayed in Figures7(a) and 7(b), respectively. The system contains two cycles: a refrigeration cycle and a chilled water cycle. The cooling capacity of the system is 3.06 kW with compressor power input of 1.04 kW. A Danfoss plate heat exchanger with 3.06 kW heat exchanging rate is selected as the evaporator. Chilled water driven by a pump is circling in the cycle with a maximum flow rate of 0.7 kg/s.(a) Photograph of the HVAC system. (b) Diagram of the HVAC system.
(a)
(b)Heat transfer occurs in the evaporator between the moving chilled water and refrigerant. As a result, the chilled water is cooled in the evaporator, as shown in Figure7. Obviously, given a fixed flow rate of chilled water, the temperature of it Tc is mainly influenced by the frequency of the compressor Fcomp. Under different frequencies of chilled water pump Fwater (different flow rates of chilled water), the transfer function from the frequency of compressor Fcomp to the temperature of chilled water Tc varies.Owing to significant effect of the frequency of chilled water pumpFwater on process dynamics, it is chosen as the scheduling variable, which is from 15 Hz to 30 Hz. The frequency of the compressor Fcomp is the input variable and the temperature of the chilled water Tc is the output variable. By varying the frequency of chilled water pump at different operating points, the process can be represented using multimodel LPV model structure.Consider three operating points at 15 Hz, 26.5 Hz, and 30 Hz for asymmetric Gaussian weight method, respectively, while the three operating points for Gaussian weight counterpart are picked as 15 Hz, 22.5 Hz, and 30 Hz. Then, three 2nd-order local linear models are identified first at the three preselected operating points. Their step responses are shown in Figure8, which implies that the dynamic behaviors of the process are different at different operating points. As can be seen from Figure 9, the input test signal is chosen as a generalized binary noise (GBN) signal and scheduling variable Fwater varies gradually from 15 Hz to 30 Hz. The output chilled water temperature is measured online at a sampling rate of 1 s. The total identification test lasted for 26579 s. During this, data obtained from the beginning 15,552 s are considered as the model training data and the other 11,027 s is used to collect validation data. The normalized weighting values using asymmetric Gaussian function and Gaussian function are given in Figure 10.Figure 8
Step responds of the three local linear models.Validation input data and variation of HVAC process scheduling variable.
(a)
(b)Normalized weighting functions for the local linear models.
(a)
(b)The comparison results are shown in Figure11. Table 2 compares the output errors of linear models and LPV models using Gaussian weights and asymmetric Gaussian weights. Step responds of the real process and LPV models are compared at other four operating points, 19, 21, 24.5, and 28.5, which are displayed in Figure 12. These results show that the estimated LPV model with asymmetric Gaussian weights is more accurate in capturing the actual nonlinear process dynamics than that with Gaussian weights.Table 2
Model output error comparison.
Models
Output error
Linear model
8.9126
%
LPV model
Gaussian weights
4.5925
%
LPV model
Asymmetric Gaussian weights
2.6249
%Figure 11
The comparison results among LPV models with Gaussian weights, LPV models with asymmetric Gaussian weights, and real process output.Figure 12
Step responds atFwater=19, 21, 24.5, and 28.5 Hz.


## 6. Conclusion

In this paper, the asymmetric Gaussian weighting function is introduced to identify the multimodel LPV model. Due to the flexibility of choosing uneven operating points for local linear models over the scheduling variable, the LPV model can be more accurate and effective. To substantiate the feasibility of the proposed approach, a simulation example of a distillation column and an experiment study with an HVAC system are performed. It has been demonstrated that with asymmetric Gaussian weighing, the accuracy of the multimodel LPV model can be improved considerably.



---

*Source: 840628-2013-12-17.xml*