- Email: [email protected]
Closed-loop identification of systems using hybrid Box–Jenkins structure and its application to PID tuning
Closed-loop identification of systems using hybrid Box-Jenkins structure and its application to PID tuning Quanshan Li, Dazi Li, Liulin Cao PII: DOI: Reference:
S1004-9541(15)00323-7 doi: 10.1016/j.cjche.2015.08.026 CJCHE 381
To appear in: Received date: Revised date: Accepted date:
25 May 2015 5 June 2015 25 June 2015
Please cite this article as: Quanshan Li, Dazi Li, Liulin Cao, Closed-loop identification of systems using hybrid Box-Jenkins structure and its application to PID tuning, (2015), doi: 10.1016/j.cjche.2015.08.026
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Received 25 May 2015 Received in revised form 5 June 2015 Accepted 25 June 2015
PT
基于混合 Box-Jenkins 混合模型的闭环系统辨识及其在 PID 整定中的应用
CR I
Closed-loop identification of systems using hybrid Box-Jenkins structure and its application to PID tuning*1 Quanshan LI (李全善)1, 2, Dazi LI (李大字)1, **, Liulin CAO (曹柳林)1 1
Institute of Automation, Beijing University of Chemical Technology, Beijing 100029, China;
2
US
Beijing Century Robust Technology Co. Ltd., Beijing 100020, China
Abstract The paper describes a closed-loop system identification procedure for hybrid continuous-time
MA N
Box-Jenkins models and demonstrates how it can be used for IMC based PID controller tuning. An instrumental variable algorithm is used to identify hybrid continuous-time transfer function models of the Box-Jenkins form from discrete-time prefiltered data, where the process model is a continuous-time transfer function, while the noise is represented as a discrete-time ARMA process. A novel penalized
D
maximum-likelihood approach is used for estimating the discrete-time ARMA process and a circulatory
PT E
noise elimination identification method is employed to estimate process model. The input-output data of a process are affected by additive circulatory noise in a closed-loop. The noise-free input-output data of the process are obtained using the proposed method by removing these circulatory noise components.
AC CE
The process model can be achieved by using instrumental variable estimation method with prefiltered noise-free input-output data. The performance of the proposed hybrid parameter estimation scheme is evaluated by Monte Carlo simulation analysis. Simulation results illustrate the efficacy of the proposed procedure. The methodology has been successfully applied in tuning of IMC based flow controller and a practical application demonstrates the applicability of the algorithm. Key Words hybrid Box-Jenkins models, ARMA models, instrumental variable, closed-loop identification, PID tuning
1 . Introduction Process models play an important role in design and implementations of several modern or traditional control theories. The model of a process to be controlled is essential in controller parameters tuning, model predictive control (MPC) and other model based control approaches [1-3]. The modelling and identification of a system have
*
Supported by the National Natural Science Foundation of China (60974031, 60704011， 61174128).
ACCEPTED MANUSCRIPT been a subject of attention to many researchers. It is well known that most physical systems are continuous-time (CT) in nature, and can conveniently be described by CT models. These models are attractive because some physical interpretation may be
PT
attached to their components, more appealing to engineers and system operators to
CR I
understand their behaviours.
Some identification algorithms are only optimal in the case of additive white noise. However, in practical situations the additive noise may not be nice white. It is likely that the
US
noise will be a coloured noise process. In this situation, it is possible to use a hybrid modelling approach, where an Auto-Regressive and Moving Average Model (ARMA) for the noise part
MA N
is estimated and used in the prefiltering operation, while the plant model would be in continuous time, as suggested in [4].
There are several approaches for ARMA process identification. In the Gaussian case, the
D
prediction error method (PEM) is an efficient estimator for ARMA process and it has been
PT E
proved that optimal instrument variable (IV) method has the same asymptotic accuracy as the PEM [5, 6]. Gaussian maximum likelihood estimation (MLE) for causal and invertible ARMA time series models has been demonstrated theoretically to be consistent and
AC CE
asymptotically normal [7].
Process model identification is the procedure of developing models from experimental
data. A good choice for model identification is applying input-output data from open-loop tests. In this situation, the noise sequence is not correlated to the input sequence. Model estimates can be directly developed without problems of inconsistency and bias [8]. For many industrial processes, production restrictions and safety are often strong reasons for not allowing identification experiments in open-loop. In such cases, experimental data can only be obtained under closed-loop conditions [9, 10]. The main difficulty in closed-loop identification is due to the correlation between the disturbances and the control signal induced by the loop [11]. Thus, most of the structures that are effective for open-loop identification result in models that are biased and not consistent in parameters. Several alternatives are available to cope with this problem [12, 13]. Although these methods have
** To whom correspondence should be addressed. E-mail: [email protected]
ACCEPTED MANUSCRIPT given effective identification results, seldom methods make fully use of sampling data of a loop for closed-loop identification. CT model identification has some difficulties because of the presence of the derivative
PT
operators associated with the input and output signals. Preprocessing input and output signals method has been proved in such a way that the undesirable derivative operations are
CR I
favourably realized [14]. Linear filter methods, due to their effectiveness and simplicity, have been widely applicated in the field of CT model identification over a considerable time
US
period. The use of state variable filters (SVF), Poisson Moment Functional filters (PMF) and refined instrumental variable approaches for continuous time model (RIVC) filters are well
MA N
known. The structures of these filters are given blow [15].
SVF: PMF: s s n
RIVC:
1 A( s)
(1)
and are filter coefficients, s is Laplace operator and A(s) is the
D
where
n 1
PT E
denominator polynomial of the identifying process model. It is confirmed that RIVC filter method is a rather simpler approach, which can not only avoid the differentiation of noisy signals but also improve statistical efficiency and
AC CE
yield lower variance estimates simultaneously [5, 16]. In this paper, RIVC filters are employed for data prefiltering. An improved
maximum-likelihood method with a penalized part is proposed for a better estimate and the facility of implementation from a practical standpoint. An improved closed-loop identification approach with IV estimate is used to identify hybrid CT Box-Jenkins (B-J) model. Convergence and consistency analysis of this class of methods can be found in [17, 18]. The proposed method is a circulatory noise elimination identification method with closed-loop sampled data. The identification data involve the setpoint input of the loop, the process input signal and the process output signal. All of them are used for closed-loop identification. The method gets noise-free input and output data of the process by removing circulatory noise component and identifies the process model using IV method with these data.
ACCEPTED MANUSCRIPT 2 . Preliminaries 2.1 Hybrid CT B-J model identification
There are some difficulties in estimation of a CT process with coloured noise,
PT
because of the theoretical and practical problems associated with the estimation of purely
CR I
stochastic continuous-time AR or ARMA models. A hybrid B-J form has been used to identify this process successfully, where the process model is a CT transfer function, while the noise is represented as a discrete-time ARMA process [4, 16].
US
For simplicity of presentation, the formulation and solution of the CT estimation
MA N
problem will be restricted to the case of a linear, single-input and single-output system. A time-delay on the system input is not considered for simplicity here but is easy to accommodate. Consider a continuous time model described by Laplace transfer function
D
of the form
PT E
G( s)
n n 1 B( s ) b0 s b b1s b bnb na A( s ) s a1s n a 1 an a
(2)
where ai , i 1, , na and bi , i 0, , nb are the coefficients, nb and n a are the
AC CE
orders of the numerator and denominator, respectively ( na nb ), and s is Laplace operator. Such that, the input u(t ) and the noise-free output response x (t ) is described by
b0 p nb b1 p nb 1 bnb B( p ) x ( t ) G ( p )u ( t ) u(t ) u(t ) A( p) p na a1 p na 1 a na
(3)
where G( p ) , A( p ) and B( p ) are operators corresponding to G(s ) , A(s ) and , B(s ) , respectively. A( p ) and B( p ) are assumed to be coprime. The additive coloured noise (t ) has rational spectral density, and then a suitable parametric representation is the ARMA model in the form of Eq. (4).
1 c1 z 1 cnc z nc D( z 1 ) (t ) w(t ) w(t ) C ( z 1 ) 1 d1 z 1 d nd z nd
(4)
ACCEPTED MANUSCRIPT where w(t ) is a white noise process with zero mean and variance 2 , and z 1 denotes the backward-shift operator z 1 (k ) (k 1) .
B( p ) D( z 1 ) u (t ) w(t ) A( p) C ( z 1 )
(5)
CR I
y (t ) G( p )u(t ) (t )
PT
The hybrid B-J models can be described as
In most practical situations, the input and output signals u(t ) and y (t ) will be
US
sampled in discrete time. In the case of uniform sampling with a constant interval Ts ,
MA N
these sampled signals will be denoted by u(tk ) and y (tk ) . The output observation equation then takes the form
y (t k ) x (t k ) (t k )
B( p ) D( z 1 ) u(t k ) w(t k ) k 1,, N A( p) C ( z 1 )
(6)
D
2.2 Instrumental variable estimation
PT E
The instrumental variable parameter estimation, with some advantages of unbiased, consistent and asymptotically normal, has been frequently used in applications such as
AC CE
system identification, time series analysis and econometrics [19]. An instrumental variable optimization procedure can also be employed to iteratively adjust the parameters that characterize the unknown process model in the hybrid B-J model, until they converge to an optimal solution. From Eq. (6), it can be concluded that
w(tk ) where Eq. (7) can be further rewritten as
C ( z 1 ) B( p ) y ( tk ) u ( tk ) 1 D( z ) A( p )
(7)
ACCEPTED MANUSCRIPT
C ( z 1 ) A( p) B( p) y (tk ) u (tk ) 1 D( z ) A( p) A( p) n n 1 b0 p nb b1 p nb 1 bnb C ( z 1 ) p a a1 p a ana y ( t ) u (tk ) k 1 D( z ) A( p) A( p)
C ( z 1 ) C ( z 1 ) na p y ( t ) a y (tk ) k n a D( z 1 ) A( p) D( z 1 ) A( p) C ( z 1 ) C ( z 1 ) nb p u ( t ) b u (tk ) k nb D( z 1 ) A( p) D( z 1 ) A( s)
CR I
b0
(8)
PT
w(tk )
yf( na ) (tk ) ana yf (tk ) b0u (fnb ) (tk ) bnb uf (tk ) (i )
(i )
US
where yf (tk ) and uf (tk ) are sampled values of the CT prefiltered variables obtained in the following manner:
C ( z 1 ) pi D( z 1 ) A( p )
(9)
MA N
f (i )
(10)
uf(i ) (t k ) f (i)u(tk )
(11)
D
yf(i ) (t k ) f (i) y(tk )
PT E
The associated, linear-in-the-parameter estimation model then takes the form
yf( na ) (tk ) a1 yf( na 1) (tk ) b0uf( nb ) (tk )
ana yf(0) (tk )
bnb uf(0) (tk ) w(tk )
(12)
AC CE
The estimation model can be written in the form
yf( na ) (tk ) f (tk )θ w(tk )
(13)
where
f (tk ) [ yf( na 1) (tk )
yf(0) (tk ) uf( nb ) (tk )
θ [a1
bnb ]T
ana b0
uf(0) (tk )]
(14)
The IV variable xˆ (t ) is generated from the following “auxiliary model”
x (t )
B( p) u (t ) A( p)
(15)
An associated IV vector ˆf (tk ) is defined as follows:
ˆf (tk ) [ xf( n 1) (tk ) a
xf(0) (tk ) uf( nb ) (tk )
uf(0) (tk )]
(16)
x (fi ) (t k ) f (i ) x(tk )
(17)
The IV optimization problem can now be stated in the form
ACCEPTED MANUSCRIPT N θˆ j ( N ) ˆ T (tk ) (tk ) k 1
1 N
ˆ k 1
T
(tk ) yf( na ) (tk )
(18)
iteration
based
on
the
appropriately
prefiltered
D N u(tk ); y(tk )k 1 .
input
and
output
data
CR I
N
3. A penalized maximum-likelihood method
PT
where the θˆ j ( N ) is the IV estimation of the system model parameter vector at the jth
US
The instrumental variable algorithms are not statistically efficient (minimum variance) in the situation that the identification data involved coloured noise and
MA N
prefiltered [4, 16]. The noise estimation is necessary to improve the statistical efficiency of the parameter estimation. On the basis of research on the ARMA identification algorithms, an efficient maximum-likelihood method with a penalized component is
D
proposed to obtain an accurate estimate of ARMA noise model.
PT E
Consider an ARMA( nc , nd ) process
1 d1 z 1 d nd z nd D( z 1 ) w ( k ) w(k ) C ( z 1 ) 1 c1 z 1 cnc z nc
AC CE
(k )
(19)
The corresponding difference equation is
(k ) cn (k nc ) w(k ) d n w(k nd ) c
d
(20)
The following assumptions are made [20, 21]: A1. C (z ) and D(z ) have no common roots. In other words, ( nc , nd ) are the
minimal orders of the ARMA model in the form of Eq. (19). It is also assumed that the orders ( nc , nd ) are known. A2. The ARMA representation in Eq. (19) is causal and invertible, which means C ( z ) 0 for all z 1 and D( z ) 0 for all z 1 . Then, the Maximum Likelihood estimator βˆ can be obtained [22] by
ACCEPTED MANUSCRIPT 1 βˆ arg min ( L( β )) arg min (ln( N β β where βˆ [c1
(21)
dnb ]T is maximum likelihood estimator, ˆ(k ) is
cna d0
(k ) and rk 1 is the mean squared error of prediction
PT
the best one-step predictor of
( (k ) ˆ(k )) 2 )) rk 1 k 1 N
CR I
which is denoted as
rk 2 E ( (k 1) ˆ(k 1)) 2
(22)
US
A drastic simplification in the calculation of one-step predictors ˆ( k 1)
MA N
and rk can be achieved by the following Eqs. [23, 24]
(23)
PT E
D
ˆ( k 1) 0 k 0 k ˆ( k 1) ( ( k 1 j ) ˆ( k 1 j )), 1 k m k, j j 1 ˆ ( k 1) c1 ( k ) cnc ( k 1 n c ) nd k , j ( ( k 1 j ) ˆ( k 1 j )) k m j 1
The coefficients k , j and mean squared errors rk are achieved recursively from Eq.
AC CE
(24) with K (i , j ) defined in Eq. (25). r0 K (1,1) h 1 1 k ,k h rh K (k 1, h 1) h ,h j k ,k j r j h 0,1,, k 1 j 0 k 1 rk K (k 1, k 1) k2,k j r j j 0
(24)
The autocovariances K (i, j ) are then found from 1 i, j m (i j ), nc (i j ) c ( h | i j |) min( i, j ) m max( i, j ) 2m h h 1 K (i , j ) nd d d min( i, j ) m h h |i j | h 0 0 otherwise
where the convention d h 0 for h nd is adopted.
(25)
ACCEPTED MANUSCRIPT The autocovariance function
(k ) based on the difference equations is obtained by
multiplying each side of Eq. (20) by
2 (t k ) and taking expectations, namely:
PT
(k ) c1 (k 1) c nc (k n c ) d j j k 0 k max( n c , n d 1) k j nd (k ) c1 (k 1) c nc (k n c ) 0 k max( n c , n d 1)
(26)
The coefficients j is determined by the relation ( z 1 ) z j D( z 1 ) / C ( z 1 ) j
CR I
j 0
j
(defining d 0 1 ， d j 0 for j nd
US
and is obtained by equating coefficients of z and c j 0 for j nc ).
MA N
j ch j h d j 0 h j j ch j h 0 0 h nc
0 j max( nc , nd 1) j max( nc , nd 1)
(27)
D
The estimator βˆ of β can be obtained using the penalized maximum-likelihood
PT E
(PML) method by minimizing the criterion function defined as follows: 2 N 1 N ( (k ) ˆ(k )) 2 1 ˆ(k )) ( 1 ) ln 1 ( ( k ) (28) 2 N k 1 rk 1 2 N k 1 J 1 ( β ) (1 ) J 2 ( β )
AC CE
J ( β)
where (0.9,1) is a factor given by 1
1 0.9 * Itr , and Itr and Itrmax are the Itrmax
current number and the maximum number of iterations, respectively. Apply a Taylor expansion at β * to Eq. (28), it can be derived T
J ( β*) J ( β ) J ( β*) ( β β*) 0( β ) β
(29)
Differentiating J ( β ) partially with respect to β , the following equations are satisfied:
J ( β ) J ( β*) 2 J ( β*) ( β β*) 0 2 β β β
(30)
1
2 J ( β*) J ( β*) β β * 2 β β
(31)
ACCEPTED MANUSCRIPT Lemma 1: As n , for k 1,, nc nd ,
PT
1 n ( (k ) ˆ(k )) 2 rk 1 P 0 n k 1 rk21 β Proof: It follows from the argument of Borckwell and Davis (10.8.47) [23] that for
CR I
1 k nc nd , with βˆ β0 , for any 0 and βˆ β0 for all large n , it can be formulated that
US
rk 1 ( ) k K1C1 ˆ
(32)
MA N
where K1 (0, ) is positive constant and 0 C1 1 . Therefore, from Eq. (22) and the Cauchy-Schwarz inequality
(33)
D
1 n ( (k ) ˆ(k )) 2 rk 1 1 P E KK1C1k 0 2 n rk 1 β ˆ n ββ k 1
AC CE
PT E
where K 0 is a constant. Thus,
1 n ( (k ) ˆ(k )) 2 rk 1 P 0 n k 1 rk21 β
(34)
It can be derived from Eq. (28) and lemma 1 with ignoring of high-order partial
derivatives that
J 1 ( β ) 1 β 2N
1 N
[( ( k ) ˆ( k )) 2 / rk 1 ] β k 1 N
(35)
( ( k ) ˆ( k )) ( ( k ) ˆ( k )) rk 1 β k 1 N
N
J 2 ( β ) β
( (k ) ˆ(k ))
N 1 ( (k ) ˆ(k )) k 1
2 J1 ( β) 1 β 2 N
N
( (k ) ˆ(k )) β k 1 N
k 1
2
1 ( (k ) ˆ(k )) ( (k ) ˆ(k )) T r β β k 1 k 1
(36)
N
(37)
ACCEPTED MANUSCRIPT
2 J 2 ( β) β 2
2 n 1 ( ( k ) ˆ( k )) k 1 2 2 n N 1 ( ( k ) ˆ( k )) k 1 T
PT
N ( ( k ) ˆ( k )) N ( ( k ) ˆ( k )) β β k 1 k 1
(38)
CR I
Such that the following iterative algorithm is achieved
1
(39)
US
2 J ( β*) J ( β*) β j β j 1 2 β β
MA N
where β j is the estimator vector at the jth iteration.
The penalized maximum-likelihood (PML) method is able to obtain a more accurate estimation of the stationary ARMA noise model. It involves an iterative relaxation approach to optimization. The main steps in the algorithm are as follows:
D
Step 1) Estimate a high order AR model with the golden section method.
PT E
ˆ (k ) associated with the high order AR Step 2) Generate the residual error series w model and apply conventional least squares manner to generate initial
AC CE
1 ˆ ( z 1 ) polynomials. estimates of the Cˆ ( z ) and D
Step 3) Compute the best one-step predictors ˆ( k ) of
(k ) and rk 1 based on
Eqs. (23) and (24).
1 ˆ ( z 1 ) in Step 4) Obtain the j th estimations of the polynomials Cˆ j ( z ) and D j
an ARMA noise model according to the recursion formula given in Eq. (39). Step 5) Calculate the j th cost function J j , if J j is less than a predefined threshold, or there is no noticeable change for the cost function then stop, else turn the algorithm to Step 3.
ACCEPTED MANUSCRIPT 4. Closed-loop CT model identification 4.1 Closed-loop hybrid models
Consider a stable, linear, SISO feedback control system consisting of a PID
PT
controller Gc ( s) , a process Gm ( s) and a noise process Gn ( s) , as shown in Fig. 1 [25,
MA N
US
CR I
26].
Figure 1 Closed-loop configuration
In Fig. 1, r (s ) is the setpoint input to the closed-loop system. u(s ) is the output of
D
the controller and is affected by additive circulatory noise. The noisy measured
PT E
output y (s ) is the sum of process output and the additive ARMA noise (s ) . w(s ) is the zero mean, normally distributed white noise source to the ARMA noise model. The
AC CE
closed-loop system is assumed to be given by the following relations:
y ( s) Gm ( s)u ( s) ( s) u ( s) r ( s)Gc ( s) y( s)Gc ( s) ( s) G ( s) w( s) n
(40)
The process is denoted by n n 1 B( s) b0 s b b1s b bnb G m ( s) A( s) s na a1s na 1 ana
(41)
where ai , i 1,, na and bi , i 0,, nb are the coefficients, nb and na are the orders of the numerator and denominator, respectively ( na nb ), and s is Laplace operator. The transfer function in time-domain corresponding to the transfer function Gm ( s) is described by
Gm ( p)
n n 1 B( p) b0 p b b1 p b A( p) p na a1 p na 1
bnb ana
(42)
ACCEPTED MANUSCRIPT The controller and the noise model are represented in the following forms: n E ( p) ene p e e1 p 1 n F ( p) fn f p f f0 f0
1 D( z 1 ) 1 d1 z Gn ( z ) C ( z 1 ) 1 c1 z 1
d nd z nd
(44)
cnc z nc
CR I
1
(43)
PT
Gc ( p)
It is assumed that all of the CT signals r(t ) , u(t ) and y (t ) are sampled at a
US
uniform sampling interval T , and the coloured disturbance (t ) is affected the closed loop. The setpoint input signal r (t ) is also assumed to be uncorrelated with the noise
MA N
disturbance (t ) . The CT model identification problem can be described as finding the process model parameters from finite data sets r (tk )k 1 , u(tk )k 1 and y (tk )k 1 . The N
N
N
hybrid model is then formulated as
D
y (tk ) Gm ( p, θ )u (tk ) Gn ( z 1 , β ) w(tk ) ana
b0
bnb ]T
β [c1
cnc
d1
d nd ]T
PT E
θ [a1
(45)
AC CE
where θ is process model parameter vector, β is ARMA model parameter vector and
w(tk ) is a zero-mean, normally distributed, DT white noise sequence. With these notations and assumptions, the closed-loop system can be described as
GmGc 1 y (tk ) 1 G G r (tk ) 1 G G (t k ) m c m c G G c c u (t ) r (t ) (t k ) k 1 GmGc k 1 GmGc
(46)
The noise from model can be derived from Eq. (46) as
GmGc x(tk ) 1 G G r (tk ) m c v(t ) Gc r (t ) k 1 GmGc k
(47)
ACCEPTED MANUSCRIPT where v(tk ) is an underlying noise-free input to the system, and x(tk ) is the corresponding output from the system which is not contaminated by noise. They can be
PT
used as sources on instrumental variables 4.2 A novel two-stage closed-loop identification
CR I
In order to obtain better estimations of the process models from sampled closed-loop system data, a circulatory noise elimination identification method (ECIV) is
US
presented based on the analysis of the structure of Fig. 1. Referring to a traditional closed-loop system with coloured noise, it can be seen that the process input generated by
MA N
the controller is affected by additive circulatory noise. The output from the process is also affected by additive circulatory noise. By removing these circulatory noise components, the noise-free input and output of the process can be obtained. The process
D
model can be achieved by using IV method with prefiltered noise-free input-output data.
PT E
The novel identification scheme is carried out by two-stage estimation procedures. The initial values of the process model are obtained in the first stage. In the second stage, the input v (tk ) and output yc (tk ) are first calculated, and both of them do not include
AC CE
the circulatory noise component. Subsequently, the process model between v (tk ) and yc (tk ) is estimated. Consider a closed-loop system, which has the structure as shown in Fig. 1, can be
represented by the following model:
y (tk ) Gmu (tk ) (t k ) u (tk ) Gc r (tk ) Gc y (t k )
(48)
The following expression can be derived:
GmGc GmGc y (tk ) 1 G G r (tk ) 1 G G (t k ) (t k ) m c m c G G c c u (t ) r (t ) (t k ) k 1 GmGc k 1 GmGc
(49)
Define the circulatory noise input vn (tk ) and circulatory noise output x n (tk ) as
ACCEPTED MANUSCRIPT Gc v n (tk ) 1 G G (t k ) m c x (t ) GmGc (t ) k n k 1 GmGc
PT
(50)
From Eqs. (47), (49) and (50) it is not difficult to get
CR I
y ( t k ) x ( t k ) xn ( t k ) ( t k ) u(tk ) v(tk ) vn (tk )
(51)
US
Substituting u(tk ) in Eq. (51) into the following equation
y(tk ) Gmu(tk ) (t k )
(52)
MA N
The process model without circulatory noise can be obtained:
yc (tk ) Gm ( p)v(tk ) (t k )
(53)
The proposed ECIV algorithm based on Eqs. (52) and (53) is an instrumental
D
variable estimation for identifying hybrid B-J models from closed-loop DT prefiltered
PT E
data. The first stage of the scheme based on Eq. (52) can be described in steps as follows:
ˆ ( s) and Bˆ ( s) Step 1) Initialization: generate initial estimates of the A 1 1
AC CE
polynomials using some estimation methods such as LSSVF method [27].
ˆ ( p) Step 2) Set j 2 , iteratively estimate the process model ploynominals A j and Bˆ j ( p ) .
Step 3) Calculate x(t )
ˆ ( p) by x(t ) Gmu(t ) using A j 1
and Bˆ j 1 ( p ) ,
respectively. 1
1
Step 4) Obtain estimates of the polynomials C j ( z ) and D j ( z ) in a DT ARMA( nc , nd ) noise model based on the estimated noise sequence
(tk ) y(tk ) x(tk ) , using a standard ARMA model estimation algorithm.
ACCEPTED MANUSCRIPT (i )
(i )
Step 5) Generate both the prefiltered input uf (t k ) , output yf (t k ) and (i )
instrumental variable xf (t k ) by equations (10), (11) and (17).
PT
ˆ ( s ) and Bˆ ( s ) using the Step 6) Obtain updated estimates of the polynomials A j j IV algorithm based on these prefiltered data by equation (18).
CR I
Step 7) Set j j 1 and calculate the variance of the residuals. If the variance does not have noticeable change, then stop the procedure. Otherwise, go to
US
Step 3).
MA N
The second stage can be carried out after the first stage is finished and the initial parameter vector θˆ has been obtained. The second stage consists of the following steps: Step 1) Generate estimates of x (tk ) and v (tk ) by Eq. (47) using parameter
D
vector θˆ j 1 .
PT E
Step 2) use Eq. (50) to calculate xn (tk ) and vn (tk ) . Step 3) Compute yc (tk ) according to yc (tk ) y(tk ) xn (tk ) .
AC CE
Step 4) Generate the prefiltered variables ycfc(i ) (t k ) , vfc(i ) (t k ) and xfc(i ) (t k ) from yc (tk ) , v (tk ) and x (tk ) by CT filter f c (i) p / A( p) i
i na .
Step 5) Obtain an optimal estimate of the noise model parameter vector βˆ j based on the estimated noise sequence
(tk ) y(tk ) x(tk ) using proposed
PML algorithm. Step 6) Use the estimated noise model parameters in βˆ j to define the DT 1
1
filter fd 1/ Gn C ( z ) / D( z ) . Then, compute ycf(i ) (t k ) , vf(i ) (t k ) and
xf(i ) (t k ) from ycf(i ) (t k ) , vf(i ) (t k ) and xf(i ) (t k ) by f d .
ACCEPTED MANUSCRIPT Step 7) Based on these prefiltered data, generate an updated estimate θˆ j of the process model parameter vector by Eq. (18). The prefiltered process vector
f (tk ) [ ycf(na 1) (tk )
ycf(0) (tk ) vf(nb ) (tk )
vf(0) (tk )] (54)
ˆf (tk ) [ xf(na 1) (tk )
CR I
PT
f (tk ) and instrumental variable vector ˆf (tk ) are described as
vf(0) (tk )] (55)
xf(0) (tk ) vf(nb ) (tk )
US
Step 8) Set j j 1 and calculate the variance of the residuals. If the variance does not have noticeable change, then, stop the procedure. Otherwise, go to
MA N
Step 1).
5. Simulation studies
An example is given to illustrate the effectiveness of the proposed algorithm. The
D
simulation model is a closed-loop transfer function model with ARMA process. The
PT E
hybrid B-J model estimation is evaluated by comparative Monte Carlo Simulation (MCS) analysis based on 100 random realizations involving the following ARMA
AC CE
model:
Gn ( z 1 )
1 0.9 z 1 1 0.8z 1
(56)
G c ( s)
1.6s 0.8 2s
(57)
The controller is given by
The process model is described as
G m ( s)
0.5s 1 s 0.2s 4 2
(58)
In this example, the white noise variance is 1 and the setpoint input r (t ) is a 2
binary ±1, piece-wise constant excitation signal. The system response can be calculated exactly at the sampling instances via appropriate Zero-Order Hold (ZOH) discretization of the closed-loop system. The sampling interval is chosen as t 0.01 s and the number of data points N 3000 . The initialized values of IV routine are obtained from
ACCEPTED MANUSCRIPT the LSSVF method. The cut-off frequency of the SVF filter is set to 10 rad.sec-1. The variance of the additive noise
(tk ) is adjusted to obtain a signal-to-noise ratio of 10
PT
dB. The comparisons of the MCS estimation results are shown in Table 1 and Table 2.
CR I
In order to compare the statistical performance of the different approaches, the computed mean βˆ and the standard error SE of the ensemble about the mean estimates are
US
presented, as well as the mean-absolute error ( MAE ) described as
oj ˆ j 100 oj
MA N
MAE
where ˆ j is the j th element of the estimated parameter vector, while true value of the parameter.
(59)
oj denotes the
D
Three estimate results are presented in Tables 1 and 2. Some are obtained from the
PT E
proposed PML algorithm, while others achieved from more conventional maximum likelihood estimation using the PEM algorithm in the Matlab System Identification
AC CE
Toolbox and the third results are attained by IVARMA or CLRIVC in CONSTID [5]. Table 1 Comparisons of MCS estimation results of ARMA Method
c1
d1
True Values
-0.8
0.9
PML
-0.798
0.900
SE
0.013
0.009
MAE
0.228
PEM
-0.799
0.898
SE
0.013
0.010
MAE
0.307
IVARMA
-0.797
0.898
SE
0.019
0.062
MAE
0.597
Form all of these results, it can be seen that the PML algorithm has a somewhat better estimates than those of PEM method and IVARMA approach. Table 2 Comparisons of MCS estimation of three different methods parameter
a1
a2
b0
b1
ACCEPTED MANUSCRIPT True value
0.2 0.20 1 0.02 9
ECIV
SE
4
-0.5
4.000
-0.501
0.022
0.036
MAE
1 1.00 0 0.02 3
0.4462 -0.504
0.043
0.027
MAE
4.5459 0.20 2 0.03 5
PEM
4.041
-0.504
0.233
0.043
0.96 8 0.18 0
US
SE
0.99 6 0.05 0
PT
SE
4.005
CR I
CLRIVC
0.20 3 0.02 1
MAE
18.6884
From Table 2, it is shown that all of these methods deliver the same good results.
MA N
However, as expected, the ECIV method improves the statistical efficiency and also provides consistent estimates of the B-J model parameters.
6. Engineering Applications
D
Because PID controllers are widely used, tuning methods for them have been
PT E
actively researched for the past several decades. Among all methods, the internal model control (IMC) based tuning methodology has gained wide spread acceptance by
AC CE
industrial users. This is because it has only one tuning parameter, which can obviously indicate the relationship between tuning parameters, closed-loop response and robustness to guarantee control performance [28, 29]. It is known that model-based PID tuning schemes, such as IMC-PID, are virtually
dependent on the identification models. Therefore, it is needed to achieve accurate models from industrial process data. This section is devoted to the research of a practical case to show the performance of the proposed method in process industries. From an understanding of the physical process such as some flow control system, it is known that time delay of the process can be negligible and the process can be described as a stable linear time-invariant system. In order to obtain an accurate and simple model for tuning a PI controller, both the orders of CT Model and ARMA model are assumed to be 2 and 1, respectively. The hybrid B-J model is described as follows:
ACCEPTED MANUSCRIPT y (t )
b0 p b1 1 d1 z 1 u ( t ) w(t ) p 2 a1 p a2 1 c1 z 1
(60)
In the following the closed-loop data sets collected from an ethylene unit are chosen
PT
for illustration. The input-output signals are sampled uniformly with a sampling interval
CR I
0.2 minute and the observation time is 30 minutes. Fig. 2 depicts the input-output data
D
MA N
US
sets of the loop 400FIC0808. It is seen that their measurements are somewhat noisy.
PT E
Figure 2 Process input and output data of loop 400FIC0808. Upper diagram: setpoint and output signal; lower diagram: controller output signal.
With these input-output data obtained from the DCS system, the ECIV
AC CE
identification algorithm is carried out. The identification results of the loop are demonstrated in Table 3. Fig. 3 depicts the output curves and the estimated result curves. In the figure, SP are setpoint values, PV are observation values and OBJ are identification results. From Fig. 3, it can be seen that the proposed method gives a better estimate for the flow control system. Hereafter, with these estimated parameters, it is not a difficult task to tune PI controller based on IMC rule. Table 3 Closed-loop identification results Loop
a1
a2
b0
b1
c1
d1
400FIC0808
3.36
5.11
3.78
5.80
0.81
0.24
Figure 3 Identification result curves of close loop 400FIC0808
ACCEPTED MANUSCRIPT Based on the identification process models, the IMC controller Gimc has the following structure:
(61)
PT
s 2 a1s a2 Gimc ( s) (b0 s b1 )( s 1)
CR I
where a1 , a2 , b0 and b1 are parameters of process model, is IMC constant and s is Laplace operator.
US
The PID controller Gc can be obtained from DCS system and its formula is described as
MA N
Gc kc (1
1 Td s 1 )( ) Ti s 1 0.1Td s
(62)
where kc is controller gain and Ti is integral time constant, and Td is derivative time
D
constant.
PT E
Utilizing an IMC based PID Controller tuning strategy to closed-loop 400FIC0808, the IMC and PID parameter are achieved and listed in Table 4. Fig. 4 are step responses of PID tuning results. In Fig. 4, SP are setpoint values, PVS is the
AC CE
closed-loop step response with old PID parameters, and PVN is step response with new PID parameters.
Table 4 IMC-PID tuning results of loop 400FIC0808
Loop
kc
Ti
Td
400FIC0808
0.50
0.44
0.25
0.00
Figure 4 PID tuning result curves of loop 400FIC0808
As illustrated in Fig. 4, the controller with new PID parameters has a better performance than the old ones. The new tuning PID parameters are set in the controllers
ACCEPTED MANUSCRIPT and the loop performs quite satisfactorily. Fig. 5 demonstrates the control performance of
MA N
US
CR I
PT
loop 400FIC0808 before and after PID tuning.
Figure 5 Comparisons of performing result for loop 400FIC0808 (top panel: before PID tuning; lower panel: after PID tuning; SP: setpoit values; PV: measured values)
It can be seen from Fig. 5 that the loop 400FIC0808 control effectiveness increases
D
obviously through the observation data. The performing results show that the tuning
PT E
parameters are suitable for the controller to control the process object. The result also indicates the accuracy of the process model identification. Meanwhile, the performance further proves the effectiveness and practical application ability of the closed-loop
AC CE
Box-Jenkins model identification algorithm.
7. Conclusions
This paper has developed a closed-loop Box-Jenkins models identification method.
From simulation results obtained in the ARMA case and the open-loop case, it is clear that the proposed PML method, as expected, performs somewhat better than alternative PEM and IVARMA approaches. The Monte Carlo Simulation studies of closed-loop case demonstrate the effectiveness of the ECIV algorithm. Finally, a real-world application has demonstrated the effectiveness and practical application ability of the proposed method.
REFERENCE 1 Raghunath Bajarangbali, Somanath Majhi, Saurabh Pandey, “Identification of FOPDT and SOPDT process dynamics using closed loop test”, ISA Trans., 53(4): 1223-1231 (2014). 2 CONG Erding, HU Minghui, TU Shandong , SHAO Huihe, “A new optimal control system design for chemical processes”, Chinese J. Chem. Eng., 21(12): 1341-1346 (2013).
ACCEPTED MANUSCRIPT 3 WANG Hong, XIE Lei, SONG Zhihuan , “A review for model plant mismatch measures in process monitoring”, Chinese J. Chem. Eng., 20(6), 1039-1046 (2012). 4 Young Peter J W, Garnier Hugues, Gilson Marion, “An optimal instrumental variable approach for identifying hybrid continuous-time Box-Jenkins models”, In: Proceedings of the 14th IFAC Symposium on System Identification( SYSID06), Newcastle, 225-230 (2006).
PT
5 Young Peter C, “An instrumental variable approach to ARMA model identification and estimation”, In: Proceedings of the 14th IFAC Symposium on System Identification( SYSID06), Newcastle,
CR I
410-415 (2006).
6 Stoica Petre, Soderstrom Torsten, “Friedlander Benjamin. Optimal instrumental variable estimates of the AR parameters of an ARMA process”, IEEE Trans. Automat. Contr., 30(11): 1066-1074 (1985).
US
7 Yao Qiwei, Brockwell Peter J., “Gaussian maximum likelihood estimation for ARMA models I Time Series”, J. Time Ser. Anal., 27(6): 857-875 (2006).
MA N
8 Lemma D., Ramasamy M., “Closed-loop Identification of Systems with Uncertain Time Delays using ARX-OBF Structures”, J. Process Control, 21(6): 1148-1154 (2011). 9 ZHANG Cong, YANG Fan, YE Hao, “Informative property of the data set in a single-input single-output (SISO) closed-loop system with a switching controller”, Chinese J. Chem. Eng., 20(6) 1128-1135 (2012).
10 Miroslav R. Mataušek, Tomislav B. Šekara, “A fast closed-loop process dynamics
D
characterization”, ISA Trans., 53(2): 489-496 (2014).
PT E
11 Tóth Roland, Laurain Vincent, Gilson Marion,“Garnier Hugues. Instrumental variable scheme for closed-loop LPV model identification”, Automatica, 48(9): 2314-2320 (2012). 12 Qibing Jin, Zhu Wang, Ruigeng Yang, Jing Wang, “An effective direct closed loop identification method for linear multivariable systems with colored noise”, J. Process Control, 24(5): 485-492
AC CE
(2014).
13 Mathieu Pouliquen, Olivier Gehan, Eric Pigeon, “Bounded-error identification for closed-loop systems”, Automatica, 50 (7): 1884-1890 (2014).
14 Young Peter C, “An instrumental variable method for real-time identification of a noisy process”, Automatica, 6(2): 271-287 (1970).
15 Unbehauen H., Rao G. P., “Continuous-time approaches to system identification—A survey”, Automatica, 26(1): 23-35 (1990).
16 Young Peter C, Garnier Hugues, “Identification and estimation of continuous-time rainfall-flow models”, In: Proceedings of the 14th IFAC Symposium on System Identification(SYSID06), Newcastle, 1276-1281 (2006). 17 Liu X, Wang J, Zheng W, “Convergence analysis of refined instrumental variable method for continuous-time system identification”, IET Control Theory Appl., 5(7): 868-877 (2011). 18 Gilson Marion, Garnier Hugues, Young Peter C, “Van den Hof Paul. Instrumental variable methods for closed-loop continuous-time model identification”, In: Identification of continuous-time models from sampled data, Springer, London, 133-160 (2008). 19 Stoica Petre, Jansson Magnus, “Estimating optimal weights for instrumental variable methods”, Digit. Signal Process., 11(3): 252-268 (2001). 20 Chen Han-Fu, Yang Jun-Mei, “Strongly consistent coefficient estimate for errors-in-variables models”, Automatica, 41(6): 1025-1033 (2005). 21 Song Qi-jiang, Chen Han-Fu, “Identification of errors-in-variables systems with ARMA
ACCEPTED MANUSCRIPT observation noises”, Syst. Control Lett., 57(5): 420-424 (2008). 22 Yao Qiwei, Brockwell Peter J, “Gaussian maximum likelihood estimation for ARMA models I Time Series”, J. Time Ser. Anal., 27(6): 857-875 (2006). 23 Brockwell Peter J., “Estimation for ARMA Models”, In: Time series: theory and methods, Springer, New York, 238-329 (1991).
PT
24 Brockwell Peter J., Davis Richard A., “Modeling and Forecasting with ARMA Processes”, In: Introduction to time series and forecasting, Taylor & Francis, New York, 137-178 (2002).
CR I
25 Young Peter C, “Identification of Transfer Function Models in Closed-Loop”, In: Recursive Estimation and Time-Series Analysis, Springer, Berlin, 271-287 (2012).
26 Young Peter C, Garnier Hugues, Gilson Marion, “Simple refined IV methods of closed-loop
(SYSID09), Saint-Malo, 1151-1156 (2009).
US
system identification”, In: Proceedings of the 15th IFAC Symposium on System Identification 27 Welsh James, Goodwin Graham C, Garnier Hugues, “A simple method for bias reduction in time
MA N
domain least squares parameter estimation”, In: Proceedings of the 3rd IFAC Symposium on System, Structure and Control(SSSC07), Brazil , 590-595 (2007). 28 Naik K Anil, Srikanth P, Negi Pankaj, “IMC Tuned PID Governor Controller for Hydro Power Plant with Water Hammer Effect”, Procedia Technol., 4: 845-853 (2012). 29 JIN Qibing, LIU Qie, WANG Qi, TIAN Yuqi, WANG Yuanfei, “PID Controller Design Based on the Time Domain Information of Robust IMC Controller Using Maximum Sensitivity”, Chinese J.
AC CE
PT E
D
Chem. Eng., 21(5): 529-536 (2013).
S1004-9541(15)00323-7 doi: 10.1016/j.cjche.2015.08.026 CJCHE 381
To appear in: Received date: Revised date: Accepted date:
25 May 2015 5 June 2015 25 June 2015
Please cite this article as: Quanshan Li, Dazi Li, Liulin Cao, Closed-loop identification of systems using hybrid Box-Jenkins structure and its application to PID tuning, (2015), doi: 10.1016/j.cjche.2015.08.026
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Received 25 May 2015 Received in revised form 5 June 2015 Accepted 25 June 2015
PT
基于混合 Box-Jenkins 混合模型的闭环系统辨识及其在 PID 整定中的应用
CR I
Closed-loop identification of systems using hybrid Box-Jenkins structure and its application to PID tuning*1 Quanshan LI (李全善)1, 2, Dazi LI (李大字)1, **, Liulin CAO (曹柳林)1 1
Institute of Automation, Beijing University of Chemical Technology, Beijing 100029, China;
2
US
Beijing Century Robust Technology Co. Ltd., Beijing 100020, China
Abstract The paper describes a closed-loop system identification procedure for hybrid continuous-time
MA N
Box-Jenkins models and demonstrates how it can be used for IMC based PID controller tuning. An instrumental variable algorithm is used to identify hybrid continuous-time transfer function models of the Box-Jenkins form from discrete-time prefiltered data, where the process model is a continuous-time transfer function, while the noise is represented as a discrete-time ARMA process. A novel penalized
D
maximum-likelihood approach is used for estimating the discrete-time ARMA process and a circulatory
PT E
noise elimination identification method is employed to estimate process model. The input-output data of a process are affected by additive circulatory noise in a closed-loop. The noise-free input-output data of the process are obtained using the proposed method by removing these circulatory noise components.
AC CE
The process model can be achieved by using instrumental variable estimation method with prefiltered noise-free input-output data. The performance of the proposed hybrid parameter estimation scheme is evaluated by Monte Carlo simulation analysis. Simulation results illustrate the efficacy of the proposed procedure. The methodology has been successfully applied in tuning of IMC based flow controller and a practical application demonstrates the applicability of the algorithm. Key Words hybrid Box-Jenkins models, ARMA models, instrumental variable, closed-loop identification, PID tuning
1 . Introduction Process models play an important role in design and implementations of several modern or traditional control theories. The model of a process to be controlled is essential in controller parameters tuning, model predictive control (MPC) and other model based control approaches [1-3]. The modelling and identification of a system have
*
Supported by the National Natural Science Foundation of China (60974031, 60704011， 61174128).
ACCEPTED MANUSCRIPT been a subject of attention to many researchers. It is well known that most physical systems are continuous-time (CT) in nature, and can conveniently be described by CT models. These models are attractive because some physical interpretation may be
PT
attached to their components, more appealing to engineers and system operators to
CR I
understand their behaviours.
Some identification algorithms are only optimal in the case of additive white noise. However, in practical situations the additive noise may not be nice white. It is likely that the
US
noise will be a coloured noise process. In this situation, it is possible to use a hybrid modelling approach, where an Auto-Regressive and Moving Average Model (ARMA) for the noise part
MA N
is estimated and used in the prefiltering operation, while the plant model would be in continuous time, as suggested in [4].
There are several approaches for ARMA process identification. In the Gaussian case, the
D
prediction error method (PEM) is an efficient estimator for ARMA process and it has been
PT E
proved that optimal instrument variable (IV) method has the same asymptotic accuracy as the PEM [5, 6]. Gaussian maximum likelihood estimation (MLE) for causal and invertible ARMA time series models has been demonstrated theoretically to be consistent and
AC CE
asymptotically normal [7].
Process model identification is the procedure of developing models from experimental
data. A good choice for model identification is applying input-output data from open-loop tests. In this situation, the noise sequence is not correlated to the input sequence. Model estimates can be directly developed without problems of inconsistency and bias [8]. For many industrial processes, production restrictions and safety are often strong reasons for not allowing identification experiments in open-loop. In such cases, experimental data can only be obtained under closed-loop conditions [9, 10]. The main difficulty in closed-loop identification is due to the correlation between the disturbances and the control signal induced by the loop [11]. Thus, most of the structures that are effective for open-loop identification result in models that are biased and not consistent in parameters. Several alternatives are available to cope with this problem [12, 13]. Although these methods have
** To whom correspondence should be addressed. E-mail: [email protected]
ACCEPTED MANUSCRIPT given effective identification results, seldom methods make fully use of sampling data of a loop for closed-loop identification. CT model identification has some difficulties because of the presence of the derivative
PT
operators associated with the input and output signals. Preprocessing input and output signals method has been proved in such a way that the undesirable derivative operations are
CR I
favourably realized [14]. Linear filter methods, due to their effectiveness and simplicity, have been widely applicated in the field of CT model identification over a considerable time
US
period. The use of state variable filters (SVF), Poisson Moment Functional filters (PMF) and refined instrumental variable approaches for continuous time model (RIVC) filters are well
MA N
known. The structures of these filters are given blow [15].
SVF: PMF: s s n
RIVC:
1 A( s)
(1)
and are filter coefficients, s is Laplace operator and A(s) is the
D
where
n 1
PT E
denominator polynomial of the identifying process model. It is confirmed that RIVC filter method is a rather simpler approach, which can not only avoid the differentiation of noisy signals but also improve statistical efficiency and
AC CE
yield lower variance estimates simultaneously [5, 16]. In this paper, RIVC filters are employed for data prefiltering. An improved
maximum-likelihood method with a penalized part is proposed for a better estimate and the facility of implementation from a practical standpoint. An improved closed-loop identification approach with IV estimate is used to identify hybrid CT Box-Jenkins (B-J) model. Convergence and consistency analysis of this class of methods can be found in [17, 18]. The proposed method is a circulatory noise elimination identification method with closed-loop sampled data. The identification data involve the setpoint input of the loop, the process input signal and the process output signal. All of them are used for closed-loop identification. The method gets noise-free input and output data of the process by removing circulatory noise component and identifies the process model using IV method with these data.
ACCEPTED MANUSCRIPT 2 . Preliminaries 2.1 Hybrid CT B-J model identification
There are some difficulties in estimation of a CT process with coloured noise,
PT
because of the theoretical and practical problems associated with the estimation of purely
CR I
stochastic continuous-time AR or ARMA models. A hybrid B-J form has been used to identify this process successfully, where the process model is a CT transfer function, while the noise is represented as a discrete-time ARMA process [4, 16].
US
For simplicity of presentation, the formulation and solution of the CT estimation
MA N
problem will be restricted to the case of a linear, single-input and single-output system. A time-delay on the system input is not considered for simplicity here but is easy to accommodate. Consider a continuous time model described by Laplace transfer function
D
of the form
PT E
G( s)
n n 1 B( s ) b0 s b b1s b bnb na A( s ) s a1s n a 1 an a
(2)
where ai , i 1, , na and bi , i 0, , nb are the coefficients, nb and n a are the
AC CE
orders of the numerator and denominator, respectively ( na nb ), and s is Laplace operator. Such that, the input u(t ) and the noise-free output response x (t ) is described by
b0 p nb b1 p nb 1 bnb B( p ) x ( t ) G ( p )u ( t ) u(t ) u(t ) A( p) p na a1 p na 1 a na
(3)
where G( p ) , A( p ) and B( p ) are operators corresponding to G(s ) , A(s ) and , B(s ) , respectively. A( p ) and B( p ) are assumed to be coprime. The additive coloured noise (t ) has rational spectral density, and then a suitable parametric representation is the ARMA model in the form of Eq. (4).
1 c1 z 1 cnc z nc D( z 1 ) (t ) w(t ) w(t ) C ( z 1 ) 1 d1 z 1 d nd z nd
(4)
ACCEPTED MANUSCRIPT where w(t ) is a white noise process with zero mean and variance 2 , and z 1 denotes the backward-shift operator z 1 (k ) (k 1) .
B( p ) D( z 1 ) u (t ) w(t ) A( p) C ( z 1 )
(5)
CR I
y (t ) G( p )u(t ) (t )
PT
The hybrid B-J models can be described as
In most practical situations, the input and output signals u(t ) and y (t ) will be
US
sampled in discrete time. In the case of uniform sampling with a constant interval Ts ,
MA N
these sampled signals will be denoted by u(tk ) and y (tk ) . The output observation equation then takes the form
y (t k ) x (t k ) (t k )
B( p ) D( z 1 ) u(t k ) w(t k ) k 1,, N A( p) C ( z 1 )
(6)
D
2.2 Instrumental variable estimation
PT E
The instrumental variable parameter estimation, with some advantages of unbiased, consistent and asymptotically normal, has been frequently used in applications such as
AC CE
system identification, time series analysis and econometrics [19]. An instrumental variable optimization procedure can also be employed to iteratively adjust the parameters that characterize the unknown process model in the hybrid B-J model, until they converge to an optimal solution. From Eq. (6), it can be concluded that
w(tk ) where Eq. (7) can be further rewritten as
C ( z 1 ) B( p ) y ( tk ) u ( tk ) 1 D( z ) A( p )
(7)
ACCEPTED MANUSCRIPT
C ( z 1 ) A( p) B( p) y (tk ) u (tk ) 1 D( z ) A( p) A( p) n n 1 b0 p nb b1 p nb 1 bnb C ( z 1 ) p a a1 p a ana y ( t ) u (tk ) k 1 D( z ) A( p) A( p)
C ( z 1 ) C ( z 1 ) na p y ( t ) a y (tk ) k n a D( z 1 ) A( p) D( z 1 ) A( p) C ( z 1 ) C ( z 1 ) nb p u ( t ) b u (tk ) k nb D( z 1 ) A( p) D( z 1 ) A( s)
CR I
b0
(8)
PT
w(tk )
yf( na ) (tk ) ana yf (tk ) b0u (fnb ) (tk ) bnb uf (tk ) (i )
(i )
US
where yf (tk ) and uf (tk ) are sampled values of the CT prefiltered variables obtained in the following manner:
C ( z 1 ) pi D( z 1 ) A( p )
(9)
MA N
f (i )
(10)
uf(i ) (t k ) f (i)u(tk )
(11)
D
yf(i ) (t k ) f (i) y(tk )
PT E
The associated, linear-in-the-parameter estimation model then takes the form
yf( na ) (tk ) a1 yf( na 1) (tk ) b0uf( nb ) (tk )
ana yf(0) (tk )
bnb uf(0) (tk ) w(tk )
(12)
AC CE
The estimation model can be written in the form
yf( na ) (tk ) f (tk )θ w(tk )
(13)
where
f (tk ) [ yf( na 1) (tk )
yf(0) (tk ) uf( nb ) (tk )
θ [a1
bnb ]T
ana b0
uf(0) (tk )]
(14)
The IV variable xˆ (t ) is generated from the following “auxiliary model”
x (t )
B( p) u (t ) A( p)
(15)
An associated IV vector ˆf (tk ) is defined as follows:
ˆf (tk ) [ xf( n 1) (tk ) a
xf(0) (tk ) uf( nb ) (tk )
uf(0) (tk )]
(16)
x (fi ) (t k ) f (i ) x(tk )
(17)
The IV optimization problem can now be stated in the form
ACCEPTED MANUSCRIPT N θˆ j ( N ) ˆ T (tk ) (tk ) k 1
1 N
ˆ k 1
T
(tk ) yf( na ) (tk )
(18)
iteration
based
on
the
appropriately
prefiltered
D N u(tk ); y(tk )k 1 .
input
and
output
data
CR I
N
3. A penalized maximum-likelihood method
PT
where the θˆ j ( N ) is the IV estimation of the system model parameter vector at the jth
US
The instrumental variable algorithms are not statistically efficient (minimum variance) in the situation that the identification data involved coloured noise and
MA N
prefiltered [4, 16]. The noise estimation is necessary to improve the statistical efficiency of the parameter estimation. On the basis of research on the ARMA identification algorithms, an efficient maximum-likelihood method with a penalized component is
D
proposed to obtain an accurate estimate of ARMA noise model.
PT E
Consider an ARMA( nc , nd ) process
1 d1 z 1 d nd z nd D( z 1 ) w ( k ) w(k ) C ( z 1 ) 1 c1 z 1 cnc z nc
AC CE
(k )
(19)
The corresponding difference equation is
(k ) cn (k nc ) w(k ) d n w(k nd ) c
d
(20)
The following assumptions are made [20, 21]: A1. C (z ) and D(z ) have no common roots. In other words, ( nc , nd ) are the
minimal orders of the ARMA model in the form of Eq. (19). It is also assumed that the orders ( nc , nd ) are known. A2. The ARMA representation in Eq. (19) is causal and invertible, which means C ( z ) 0 for all z 1 and D( z ) 0 for all z 1 . Then, the Maximum Likelihood estimator βˆ can be obtained [22] by
ACCEPTED MANUSCRIPT 1 βˆ arg min ( L( β )) arg min (ln( N β β where βˆ [c1
(21)
dnb ]T is maximum likelihood estimator, ˆ(k ) is
cna d0
(k ) and rk 1 is the mean squared error of prediction
PT
the best one-step predictor of
( (k ) ˆ(k )) 2 )) rk 1 k 1 N
CR I
which is denoted as
rk 2 E ( (k 1) ˆ(k 1)) 2
(22)
US
A drastic simplification in the calculation of one-step predictors ˆ( k 1)
MA N
and rk can be achieved by the following Eqs. [23, 24]
(23)
PT E
D
ˆ( k 1) 0 k 0 k ˆ( k 1) ( ( k 1 j ) ˆ( k 1 j )), 1 k m k, j j 1 ˆ ( k 1) c1 ( k ) cnc ( k 1 n c ) nd k , j ( ( k 1 j ) ˆ( k 1 j )) k m j 1
The coefficients k , j and mean squared errors rk are achieved recursively from Eq.
AC CE
(24) with K (i , j ) defined in Eq. (25). r0 K (1,1) h 1 1 k ,k h rh K (k 1, h 1) h ,h j k ,k j r j h 0,1,, k 1 j 0 k 1 rk K (k 1, k 1) k2,k j r j j 0
(24)
The autocovariances K (i, j ) are then found from 1 i, j m (i j ), nc (i j ) c ( h | i j |) min( i, j ) m max( i, j ) 2m h h 1 K (i , j ) nd d d min( i, j ) m h h |i j | h 0 0 otherwise
where the convention d h 0 for h nd is adopted.
(25)
ACCEPTED MANUSCRIPT The autocovariance function
(k ) based on the difference equations is obtained by
multiplying each side of Eq. (20) by
2 (t k ) and taking expectations, namely:
PT
(k ) c1 (k 1) c nc (k n c ) d j j k 0 k max( n c , n d 1) k j nd (k ) c1 (k 1) c nc (k n c ) 0 k max( n c , n d 1)
(26)
The coefficients j is determined by the relation ( z 1 ) z j D( z 1 ) / C ( z 1 ) j
CR I
j 0
j
(defining d 0 1 ， d j 0 for j nd
US
and is obtained by equating coefficients of z and c j 0 for j nc ).
MA N
j ch j h d j 0 h j j ch j h 0 0 h nc
0 j max( nc , nd 1) j max( nc , nd 1)
(27)
D
The estimator βˆ of β can be obtained using the penalized maximum-likelihood
PT E
(PML) method by minimizing the criterion function defined as follows: 2 N 1 N ( (k ) ˆ(k )) 2 1 ˆ(k )) ( 1 ) ln 1 ( ( k ) (28) 2 N k 1 rk 1 2 N k 1 J 1 ( β ) (1 ) J 2 ( β )
AC CE
J ( β)
where (0.9,1) is a factor given by 1
1 0.9 * Itr , and Itr and Itrmax are the Itrmax
current number and the maximum number of iterations, respectively. Apply a Taylor expansion at β * to Eq. (28), it can be derived T
J ( β*) J ( β ) J ( β*) ( β β*) 0( β ) β
(29)
Differentiating J ( β ) partially with respect to β , the following equations are satisfied:
J ( β ) J ( β*) 2 J ( β*) ( β β*) 0 2 β β β
(30)
1
2 J ( β*) J ( β*) β β * 2 β β
(31)
ACCEPTED MANUSCRIPT Lemma 1: As n , for k 1,, nc nd ,
PT
1 n ( (k ) ˆ(k )) 2 rk 1 P 0 n k 1 rk21 β Proof: It follows from the argument of Borckwell and Davis (10.8.47) [23] that for
CR I
1 k nc nd , with βˆ β0 , for any 0 and βˆ β0 for all large n , it can be formulated that
US
rk 1 ( ) k K1C1 ˆ
(32)
MA N
where K1 (0, ) is positive constant and 0 C1 1 . Therefore, from Eq. (22) and the Cauchy-Schwarz inequality
(33)
D
1 n ( (k ) ˆ(k )) 2 rk 1 1 P E KK1C1k 0 2 n rk 1 β ˆ n ββ k 1
AC CE
PT E
where K 0 is a constant. Thus,
1 n ( (k ) ˆ(k )) 2 rk 1 P 0 n k 1 rk21 β
(34)
It can be derived from Eq. (28) and lemma 1 with ignoring of high-order partial
derivatives that
J 1 ( β ) 1 β 2N
1 N
[( ( k ) ˆ( k )) 2 / rk 1 ] β k 1 N
(35)
( ( k ) ˆ( k )) ( ( k ) ˆ( k )) rk 1 β k 1 N
N
J 2 ( β ) β
( (k ) ˆ(k ))
N 1 ( (k ) ˆ(k )) k 1
2 J1 ( β) 1 β 2 N
N
( (k ) ˆ(k )) β k 1 N
k 1
2
1 ( (k ) ˆ(k )) ( (k ) ˆ(k )) T r β β k 1 k 1
(36)
N
(37)
ACCEPTED MANUSCRIPT
2 J 2 ( β) β 2
2 n 1 ( ( k ) ˆ( k )) k 1 2 2 n N 1 ( ( k ) ˆ( k )) k 1 T
PT
N ( ( k ) ˆ( k )) N ( ( k ) ˆ( k )) β β k 1 k 1
(38)
CR I
Such that the following iterative algorithm is achieved
1
(39)
US
2 J ( β*) J ( β*) β j β j 1 2 β β
MA N
where β j is the estimator vector at the jth iteration.
The penalized maximum-likelihood (PML) method is able to obtain a more accurate estimation of the stationary ARMA noise model. It involves an iterative relaxation approach to optimization. The main steps in the algorithm are as follows:
D
Step 1) Estimate a high order AR model with the golden section method.
PT E
ˆ (k ) associated with the high order AR Step 2) Generate the residual error series w model and apply conventional least squares manner to generate initial
AC CE
1 ˆ ( z 1 ) polynomials. estimates of the Cˆ ( z ) and D
Step 3) Compute the best one-step predictors ˆ( k ) of
(k ) and rk 1 based on
Eqs. (23) and (24).
1 ˆ ( z 1 ) in Step 4) Obtain the j th estimations of the polynomials Cˆ j ( z ) and D j
an ARMA noise model according to the recursion formula given in Eq. (39). Step 5) Calculate the j th cost function J j , if J j is less than a predefined threshold, or there is no noticeable change for the cost function then stop, else turn the algorithm to Step 3.
ACCEPTED MANUSCRIPT 4. Closed-loop CT model identification 4.1 Closed-loop hybrid models
Consider a stable, linear, SISO feedback control system consisting of a PID
PT
controller Gc ( s) , a process Gm ( s) and a noise process Gn ( s) , as shown in Fig. 1 [25,
MA N
US
CR I
26].
Figure 1 Closed-loop configuration
In Fig. 1, r (s ) is the setpoint input to the closed-loop system. u(s ) is the output of
D
the controller and is affected by additive circulatory noise. The noisy measured
PT E
output y (s ) is the sum of process output and the additive ARMA noise (s ) . w(s ) is the zero mean, normally distributed white noise source to the ARMA noise model. The
AC CE
closed-loop system is assumed to be given by the following relations:
y ( s) Gm ( s)u ( s) ( s) u ( s) r ( s)Gc ( s) y( s)Gc ( s) ( s) G ( s) w( s) n
(40)
The process is denoted by n n 1 B( s) b0 s b b1s b bnb G m ( s) A( s) s na a1s na 1 ana
(41)
where ai , i 1,, na and bi , i 0,, nb are the coefficients, nb and na are the orders of the numerator and denominator, respectively ( na nb ), and s is Laplace operator. The transfer function in time-domain corresponding to the transfer function Gm ( s) is described by
Gm ( p)
n n 1 B( p) b0 p b b1 p b A( p) p na a1 p na 1
bnb ana
(42)
ACCEPTED MANUSCRIPT The controller and the noise model are represented in the following forms: n E ( p) ene p e e1 p 1 n F ( p) fn f p f f0 f0
1 D( z 1 ) 1 d1 z Gn ( z ) C ( z 1 ) 1 c1 z 1
d nd z nd
(44)
cnc z nc
CR I
1
(43)
PT
Gc ( p)
It is assumed that all of the CT signals r(t ) , u(t ) and y (t ) are sampled at a
US
uniform sampling interval T , and the coloured disturbance (t ) is affected the closed loop. The setpoint input signal r (t ) is also assumed to be uncorrelated with the noise
MA N
disturbance (t ) . The CT model identification problem can be described as finding the process model parameters from finite data sets r (tk )k 1 , u(tk )k 1 and y (tk )k 1 . The N
N
N
hybrid model is then formulated as
D
y (tk ) Gm ( p, θ )u (tk ) Gn ( z 1 , β ) w(tk ) ana
b0
bnb ]T
β [c1
cnc
d1
d nd ]T
PT E
θ [a1
(45)
AC CE
where θ is process model parameter vector, β is ARMA model parameter vector and
w(tk ) is a zero-mean, normally distributed, DT white noise sequence. With these notations and assumptions, the closed-loop system can be described as
GmGc 1 y (tk ) 1 G G r (tk ) 1 G G (t k ) m c m c G G c c u (t ) r (t ) (t k ) k 1 GmGc k 1 GmGc
(46)
The noise from model can be derived from Eq. (46) as
GmGc x(tk ) 1 G G r (tk ) m c v(t ) Gc r (t ) k 1 GmGc k
(47)
ACCEPTED MANUSCRIPT where v(tk ) is an underlying noise-free input to the system, and x(tk ) is the corresponding output from the system which is not contaminated by noise. They can be
PT
used as sources on instrumental variables 4.2 A novel two-stage closed-loop identification
CR I
In order to obtain better estimations of the process models from sampled closed-loop system data, a circulatory noise elimination identification method (ECIV) is
US
presented based on the analysis of the structure of Fig. 1. Referring to a traditional closed-loop system with coloured noise, it can be seen that the process input generated by
MA N
the controller is affected by additive circulatory noise. The output from the process is also affected by additive circulatory noise. By removing these circulatory noise components, the noise-free input and output of the process can be obtained. The process
D
model can be achieved by using IV method with prefiltered noise-free input-output data.
PT E
The novel identification scheme is carried out by two-stage estimation procedures. The initial values of the process model are obtained in the first stage. In the second stage, the input v (tk ) and output yc (tk ) are first calculated, and both of them do not include
AC CE
the circulatory noise component. Subsequently, the process model between v (tk ) and yc (tk ) is estimated. Consider a closed-loop system, which has the structure as shown in Fig. 1, can be
represented by the following model:
y (tk ) Gmu (tk ) (t k ) u (tk ) Gc r (tk ) Gc y (t k )
(48)
The following expression can be derived:
GmGc GmGc y (tk ) 1 G G r (tk ) 1 G G (t k ) (t k ) m c m c G G c c u (t ) r (t ) (t k ) k 1 GmGc k 1 GmGc
(49)
Define the circulatory noise input vn (tk ) and circulatory noise output x n (tk ) as
ACCEPTED MANUSCRIPT Gc v n (tk ) 1 G G (t k ) m c x (t ) GmGc (t ) k n k 1 GmGc
PT
(50)
From Eqs. (47), (49) and (50) it is not difficult to get
CR I
y ( t k ) x ( t k ) xn ( t k ) ( t k ) u(tk ) v(tk ) vn (tk )
(51)
US
Substituting u(tk ) in Eq. (51) into the following equation
y(tk ) Gmu(tk ) (t k )
(52)
MA N
The process model without circulatory noise can be obtained:
yc (tk ) Gm ( p)v(tk ) (t k )
(53)
The proposed ECIV algorithm based on Eqs. (52) and (53) is an instrumental
D
variable estimation for identifying hybrid B-J models from closed-loop DT prefiltered
PT E
data. The first stage of the scheme based on Eq. (52) can be described in steps as follows:
ˆ ( s) and Bˆ ( s) Step 1) Initialization: generate initial estimates of the A 1 1
AC CE
polynomials using some estimation methods such as LSSVF method [27].
ˆ ( p) Step 2) Set j 2 , iteratively estimate the process model ploynominals A j and Bˆ j ( p ) .
Step 3) Calculate x(t )
ˆ ( p) by x(t ) Gmu(t ) using A j 1
and Bˆ j 1 ( p ) ,
respectively. 1
1
Step 4) Obtain estimates of the polynomials C j ( z ) and D j ( z ) in a DT ARMA( nc , nd ) noise model based on the estimated noise sequence
(tk ) y(tk ) x(tk ) , using a standard ARMA model estimation algorithm.
ACCEPTED MANUSCRIPT (i )
(i )
Step 5) Generate both the prefiltered input uf (t k ) , output yf (t k ) and (i )
instrumental variable xf (t k ) by equations (10), (11) and (17).
PT
ˆ ( s ) and Bˆ ( s ) using the Step 6) Obtain updated estimates of the polynomials A j j IV algorithm based on these prefiltered data by equation (18).
CR I
Step 7) Set j j 1 and calculate the variance of the residuals. If the variance does not have noticeable change, then stop the procedure. Otherwise, go to
US
Step 3).
MA N
The second stage can be carried out after the first stage is finished and the initial parameter vector θˆ has been obtained. The second stage consists of the following steps: Step 1) Generate estimates of x (tk ) and v (tk ) by Eq. (47) using parameter
D
vector θˆ j 1 .
PT E
Step 2) use Eq. (50) to calculate xn (tk ) and vn (tk ) . Step 3) Compute yc (tk ) according to yc (tk ) y(tk ) xn (tk ) .
AC CE
Step 4) Generate the prefiltered variables ycfc(i ) (t k ) , vfc(i ) (t k ) and xfc(i ) (t k ) from yc (tk ) , v (tk ) and x (tk ) by CT filter f c (i) p / A( p) i
i na .
Step 5) Obtain an optimal estimate of the noise model parameter vector βˆ j based on the estimated noise sequence
(tk ) y(tk ) x(tk ) using proposed
PML algorithm. Step 6) Use the estimated noise model parameters in βˆ j to define the DT 1
1
filter fd 1/ Gn C ( z ) / D( z ) . Then, compute ycf(i ) (t k ) , vf(i ) (t k ) and
xf(i ) (t k ) from ycf(i ) (t k ) , vf(i ) (t k ) and xf(i ) (t k ) by f d .
ACCEPTED MANUSCRIPT Step 7) Based on these prefiltered data, generate an updated estimate θˆ j of the process model parameter vector by Eq. (18). The prefiltered process vector
f (tk ) [ ycf(na 1) (tk )
ycf(0) (tk ) vf(nb ) (tk )
vf(0) (tk )] (54)
ˆf (tk ) [ xf(na 1) (tk )
CR I
PT
f (tk ) and instrumental variable vector ˆf (tk ) are described as
vf(0) (tk )] (55)
xf(0) (tk ) vf(nb ) (tk )
US
Step 8) Set j j 1 and calculate the variance of the residuals. If the variance does not have noticeable change, then, stop the procedure. Otherwise, go to
MA N
Step 1).
5. Simulation studies
An example is given to illustrate the effectiveness of the proposed algorithm. The
D
simulation model is a closed-loop transfer function model with ARMA process. The
PT E
hybrid B-J model estimation is evaluated by comparative Monte Carlo Simulation (MCS) analysis based on 100 random realizations involving the following ARMA
AC CE
model:
Gn ( z 1 )
1 0.9 z 1 1 0.8z 1
(56)
G c ( s)
1.6s 0.8 2s
(57)
The controller is given by
The process model is described as
G m ( s)
0.5s 1 s 0.2s 4 2
(58)
In this example, the white noise variance is 1 and the setpoint input r (t ) is a 2
binary ±1, piece-wise constant excitation signal. The system response can be calculated exactly at the sampling instances via appropriate Zero-Order Hold (ZOH) discretization of the closed-loop system. The sampling interval is chosen as t 0.01 s and the number of data points N 3000 . The initialized values of IV routine are obtained from
ACCEPTED MANUSCRIPT the LSSVF method. The cut-off frequency of the SVF filter is set to 10 rad.sec-1. The variance of the additive noise
(tk ) is adjusted to obtain a signal-to-noise ratio of 10
PT
dB. The comparisons of the MCS estimation results are shown in Table 1 and Table 2.
CR I
In order to compare the statistical performance of the different approaches, the computed mean βˆ and the standard error SE of the ensemble about the mean estimates are
US
presented, as well as the mean-absolute error ( MAE ) described as
oj ˆ j 100 oj
MA N
MAE
where ˆ j is the j th element of the estimated parameter vector, while true value of the parameter.
(59)
oj denotes the
D
Three estimate results are presented in Tables 1 and 2. Some are obtained from the
PT E
proposed PML algorithm, while others achieved from more conventional maximum likelihood estimation using the PEM algorithm in the Matlab System Identification
AC CE
Toolbox and the third results are attained by IVARMA or CLRIVC in CONSTID [5]. Table 1 Comparisons of MCS estimation results of ARMA Method
c1
d1
True Values
-0.8
0.9
PML
-0.798
0.900
SE
0.013
0.009
MAE
0.228
PEM
-0.799
0.898
SE
0.013
0.010
MAE
0.307
IVARMA
-0.797
0.898
SE
0.019
0.062
MAE
0.597
Form all of these results, it can be seen that the PML algorithm has a somewhat better estimates than those of PEM method and IVARMA approach. Table 2 Comparisons of MCS estimation of three different methods parameter
a1
a2
b0
b1
ACCEPTED MANUSCRIPT True value
0.2 0.20 1 0.02 9
ECIV
SE
4
-0.5
4.000
-0.501
0.022
0.036
MAE
1 1.00 0 0.02 3
0.4462 -0.504
0.043
0.027
MAE
4.5459 0.20 2 0.03 5
PEM
4.041
-0.504
0.233
0.043
0.96 8 0.18 0
US
SE
0.99 6 0.05 0
PT
SE
4.005
CR I
CLRIVC
0.20 3 0.02 1
MAE
18.6884
From Table 2, it is shown that all of these methods deliver the same good results.
MA N
However, as expected, the ECIV method improves the statistical efficiency and also provides consistent estimates of the B-J model parameters.
6. Engineering Applications
D
Because PID controllers are widely used, tuning methods for them have been
PT E
actively researched for the past several decades. Among all methods, the internal model control (IMC) based tuning methodology has gained wide spread acceptance by
AC CE
industrial users. This is because it has only one tuning parameter, which can obviously indicate the relationship between tuning parameters, closed-loop response and robustness to guarantee control performance [28, 29]. It is known that model-based PID tuning schemes, such as IMC-PID, are virtually
dependent on the identification models. Therefore, it is needed to achieve accurate models from industrial process data. This section is devoted to the research of a practical case to show the performance of the proposed method in process industries. From an understanding of the physical process such as some flow control system, it is known that time delay of the process can be negligible and the process can be described as a stable linear time-invariant system. In order to obtain an accurate and simple model for tuning a PI controller, both the orders of CT Model and ARMA model are assumed to be 2 and 1, respectively. The hybrid B-J model is described as follows:
ACCEPTED MANUSCRIPT y (t )
b0 p b1 1 d1 z 1 u ( t ) w(t ) p 2 a1 p a2 1 c1 z 1
(60)
In the following the closed-loop data sets collected from an ethylene unit are chosen
PT
for illustration. The input-output signals are sampled uniformly with a sampling interval
CR I
0.2 minute and the observation time is 30 minutes. Fig. 2 depicts the input-output data
D
MA N
US
sets of the loop 400FIC0808. It is seen that their measurements are somewhat noisy.
PT E
Figure 2 Process input and output data of loop 400FIC0808. Upper diagram: setpoint and output signal; lower diagram: controller output signal.
With these input-output data obtained from the DCS system, the ECIV
AC CE
identification algorithm is carried out. The identification results of the loop are demonstrated in Table 3. Fig. 3 depicts the output curves and the estimated result curves. In the figure, SP are setpoint values, PV are observation values and OBJ are identification results. From Fig. 3, it can be seen that the proposed method gives a better estimate for the flow control system. Hereafter, with these estimated parameters, it is not a difficult task to tune PI controller based on IMC rule. Table 3 Closed-loop identification results Loop
a1
a2
b0
b1
c1
d1
400FIC0808
3.36
5.11
3.78
5.80
0.81
0.24
Figure 3 Identification result curves of close loop 400FIC0808
ACCEPTED MANUSCRIPT Based on the identification process models, the IMC controller Gimc has the following structure:
(61)
PT
s 2 a1s a2 Gimc ( s) (b0 s b1 )( s 1)
CR I
where a1 , a2 , b0 and b1 are parameters of process model, is IMC constant and s is Laplace operator.
US
The PID controller Gc can be obtained from DCS system and its formula is described as
MA N
Gc kc (1
1 Td s 1 )( ) Ti s 1 0.1Td s
(62)
where kc is controller gain and Ti is integral time constant, and Td is derivative time
D
constant.
PT E
Utilizing an IMC based PID Controller tuning strategy to closed-loop 400FIC0808, the IMC and PID parameter are achieved and listed in Table 4. Fig. 4 are step responses of PID tuning results. In Fig. 4, SP are setpoint values, PVS is the
AC CE
closed-loop step response with old PID parameters, and PVN is step response with new PID parameters.
Table 4 IMC-PID tuning results of loop 400FIC0808
Loop
kc
Ti
Td
400FIC0808
0.50
0.44
0.25
0.00
Figure 4 PID tuning result curves of loop 400FIC0808
As illustrated in Fig. 4, the controller with new PID parameters has a better performance than the old ones. The new tuning PID parameters are set in the controllers
ACCEPTED MANUSCRIPT and the loop performs quite satisfactorily. Fig. 5 demonstrates the control performance of
MA N
US
CR I
PT
loop 400FIC0808 before and after PID tuning.
Figure 5 Comparisons of performing result for loop 400FIC0808 (top panel: before PID tuning; lower panel: after PID tuning; SP: setpoit values; PV: measured values)
It can be seen from Fig. 5 that the loop 400FIC0808 control effectiveness increases
D
obviously through the observation data. The performing results show that the tuning
PT E
parameters are suitable for the controller to control the process object. The result also indicates the accuracy of the process model identification. Meanwhile, the performance further proves the effectiveness and practical application ability of the closed-loop
AC CE
Box-Jenkins model identification algorithm.
7. Conclusions
This paper has developed a closed-loop Box-Jenkins models identification method.
From simulation results obtained in the ARMA case and the open-loop case, it is clear that the proposed PML method, as expected, performs somewhat better than alternative PEM and IVARMA approaches. The Monte Carlo Simulation studies of closed-loop case demonstrate the effectiveness of the ECIV algorithm. Finally, a real-world application has demonstrated the effectiveness and practical application ability of the proposed method.
REFERENCE 1 Raghunath Bajarangbali, Somanath Majhi, Saurabh Pandey, “Identification of FOPDT and SOPDT process dynamics using closed loop test”, ISA Trans., 53(4): 1223-1231 (2014). 2 CONG Erding, HU Minghui, TU Shandong , SHAO Huihe, “A new optimal control system design for chemical processes”, Chinese J. Chem. Eng., 21(12): 1341-1346 (2013).
ACCEPTED MANUSCRIPT 3 WANG Hong, XIE Lei, SONG Zhihuan , “A review for model plant mismatch measures in process monitoring”, Chinese J. Chem. Eng., 20(6), 1039-1046 (2012). 4 Young Peter J W, Garnier Hugues, Gilson Marion, “An optimal instrumental variable approach for identifying hybrid continuous-time Box-Jenkins models”, In: Proceedings of the 14th IFAC Symposium on System Identification( SYSID06), Newcastle, 225-230 (2006).
PT
5 Young Peter C, “An instrumental variable approach to ARMA model identification and estimation”, In: Proceedings of the 14th IFAC Symposium on System Identification( SYSID06), Newcastle,
CR I
410-415 (2006).
6 Stoica Petre, Soderstrom Torsten, “Friedlander Benjamin. Optimal instrumental variable estimates of the AR parameters of an ARMA process”, IEEE Trans. Automat. Contr., 30(11): 1066-1074 (1985).
US
7 Yao Qiwei, Brockwell Peter J., “Gaussian maximum likelihood estimation for ARMA models I Time Series”, J. Time Ser. Anal., 27(6): 857-875 (2006).
MA N
8 Lemma D., Ramasamy M., “Closed-loop Identification of Systems with Uncertain Time Delays using ARX-OBF Structures”, J. Process Control, 21(6): 1148-1154 (2011). 9 ZHANG Cong, YANG Fan, YE Hao, “Informative property of the data set in a single-input single-output (SISO) closed-loop system with a switching controller”, Chinese J. Chem. Eng., 20(6) 1128-1135 (2012).
10 Miroslav R. Mataušek, Tomislav B. Šekara, “A fast closed-loop process dynamics
D
characterization”, ISA Trans., 53(2): 489-496 (2014).
PT E
11 Tóth Roland, Laurain Vincent, Gilson Marion,“Garnier Hugues. Instrumental variable scheme for closed-loop LPV model identification”, Automatica, 48(9): 2314-2320 (2012). 12 Qibing Jin, Zhu Wang, Ruigeng Yang, Jing Wang, “An effective direct closed loop identification method for linear multivariable systems with colored noise”, J. Process Control, 24(5): 485-492
AC CE
(2014).
13 Mathieu Pouliquen, Olivier Gehan, Eric Pigeon, “Bounded-error identification for closed-loop systems”, Automatica, 50 (7): 1884-1890 (2014).
14 Young Peter C, “An instrumental variable method for real-time identification of a noisy process”, Automatica, 6(2): 271-287 (1970).
15 Unbehauen H., Rao G. P., “Continuous-time approaches to system identification—A survey”, Automatica, 26(1): 23-35 (1990).
16 Young Peter C, Garnier Hugues, “Identification and estimation of continuous-time rainfall-flow models”, In: Proceedings of the 14th IFAC Symposium on System Identification(SYSID06), Newcastle, 1276-1281 (2006). 17 Liu X, Wang J, Zheng W, “Convergence analysis of refined instrumental variable method for continuous-time system identification”, IET Control Theory Appl., 5(7): 868-877 (2011). 18 Gilson Marion, Garnier Hugues, Young Peter C, “Van den Hof Paul. Instrumental variable methods for closed-loop continuous-time model identification”, In: Identification of continuous-time models from sampled data, Springer, London, 133-160 (2008). 19 Stoica Petre, Jansson Magnus, “Estimating optimal weights for instrumental variable methods”, Digit. Signal Process., 11(3): 252-268 (2001). 20 Chen Han-Fu, Yang Jun-Mei, “Strongly consistent coefficient estimate for errors-in-variables models”, Automatica, 41(6): 1025-1033 (2005). 21 Song Qi-jiang, Chen Han-Fu, “Identification of errors-in-variables systems with ARMA
ACCEPTED MANUSCRIPT observation noises”, Syst. Control Lett., 57(5): 420-424 (2008). 22 Yao Qiwei, Brockwell Peter J, “Gaussian maximum likelihood estimation for ARMA models I Time Series”, J. Time Ser. Anal., 27(6): 857-875 (2006). 23 Brockwell Peter J., “Estimation for ARMA Models”, In: Time series: theory and methods, Springer, New York, 238-329 (1991).
PT
24 Brockwell Peter J., Davis Richard A., “Modeling and Forecasting with ARMA Processes”, In: Introduction to time series and forecasting, Taylor & Francis, New York, 137-178 (2002).
CR I
25 Young Peter C, “Identification of Transfer Function Models in Closed-Loop”, In: Recursive Estimation and Time-Series Analysis, Springer, Berlin, 271-287 (2012).
26 Young Peter C, Garnier Hugues, Gilson Marion, “Simple refined IV methods of closed-loop
(SYSID09), Saint-Malo, 1151-1156 (2009).
US
system identification”, In: Proceedings of the 15th IFAC Symposium on System Identification 27 Welsh James, Goodwin Graham C, Garnier Hugues, “A simple method for bias reduction in time
MA N
domain least squares parameter estimation”, In: Proceedings of the 3rd IFAC Symposium on System, Structure and Control(SSSC07), Brazil , 590-595 (2007). 28 Naik K Anil, Srikanth P, Negi Pankaj, “IMC Tuned PID Governor Controller for Hydro Power Plant with Water Hammer Effect”, Procedia Technol., 4: 845-853 (2012). 29 JIN Qibing, LIU Qie, WANG Qi, TIAN Yuqi, WANG Yuanfei, “PID Controller Design Based on the Time Domain Information of Robust IMC Controller Using Maximum Sensitivity”, Chinese J.
AC CE
PT E
D
Chem. Eng., 21(5): 529-536 (2013).