- Email: [email protected]
Continuous-time prediction error identification method for delayed processes
IFAC
Copyright © IFAC Dynamics and Control of Process Systems, Jejudo Island, Korea, 200 I
D
0
[>
Publications www.elsevier.comllocate/ifac
CONTINUOUS-TIME PREDICTION ERROR IDENTIFICATION METHOD FOR DELAYED PROCESSES Su Whan Sung l " and In-Beum Lee 2 JDepartment
of Chemical Engineering, Korean Institute of Science and Technology 373-1 Kusong-dong, Yusong-gu, Taejon, 305-701, South Korea (*E-mail:[email protected]) 2Department of Chemical Engineering, Pohang University of Science and Technology San 31 Hyoja-dong, Pohang, 790-784, South Korea
Abstract: We propose a continuous-time prediction error method to identify combined deterministicstochastic continuous-time delayed processes. It minimizes the prediction error using the Newton optimization method with analytical derivatives of the objective function with respect to the adjustable parameters that contains the time delay. The proposed method can incorporate directly both small and large sampling time as well as irregular sampling time. Unbiased estimates are guaranteed like discretetime prediction error method. Copyright © 20011FAC Key words: time delay, continuous-time, prediction error identification method approaches under certain circumstances. So, many continuous-time system identification methods have been proposed like using block pulse functions, Laguerre polynomials, Poisson moment functionals. These functions produce initial value problem. Recently, many other methods of no initial value problem have been developed: linear integral transform or macro-difference expressions (Sagara and Zhao (1989» , complex variable transform (Johansson (1994» , delta operator (Middleton and Goodwin (1990» and ball-shape weighting integral transform (Sung et al. (1998». Several papers have been presented to analyze the of numerical derivatives and effects remove/compensate the bias in the estimation problem (Soderstrom and Mossberg (2000» . Recently, Sung et al. (2000) justified the usefulness of continuous-time identification methods against discrete-time ones in control point of view and proposed a continuous-time subspace system identification method.
1. INTRODUCTION
The main focus of this paper is to develop a prediction error identification method that can identify first, time delay second, system parameters in continuous-time domain third, incorporate a small/large and irregular sampling time. In practice, the time delay is much common in chemical industries due to transportation, measurement delay, and long transmission lines in pneumatic systems. Also, the time delay has been used widely to approximate the complex and high order dynamic system together with a low order process model like the first or second order one. But, most combined deterministic-stochastic system identification methods cannot treat the time delay in a systematic manner. On the other hand, most researchers in system identification field have concentrated their efforts on developing discrete-time identification methods like prediction error, instrumental variable or subspace identification methods (Ljung (1987), Soderstrom and Stoica (1989), Larimore (1990), Verhaegen (1994), Van Overschee and De Moor (1994». But, we had better use continuous-time
The above-mentioned previous methods have several disadvantages: Previous discrete-time approaches have problems that an irregular sampling period cannot be incorporated without major modification since the adjustable parameters are not constant for the irregular sampling period.
This work is supported by the BK21 Program. 481
Also, a small sampling period and/or a continuoustime input cannot be treated with acceptable accuracy. On the other hand, previous continuoustime approaches using integral transforms or delta operator cannot manipulate a large sampling period.
o -a. o -a._, o - a._2
000 o 0
A=
Furthermore, the time delay cannot be identified in a systematic manner.
o
o
o o
0 0
0 0
b•.,
-a 2
0
-a, b•.2
b•.m
b._,., b._,.
2
In this research, we propose a continuous-time prediction error method to identify combined
(3)
B = b._2., b._2.2
deterministic-stochastic continuous-time delayed processes. Compared with previous discrete-time identification methods, the proposed method does not suffer from the small sampling period or irregular sampling problem. Also, it can be applied
Here, (1) and (2) can be rewritten equivalently like the following continuous-time Auto-Regressive, Moving Average, eXogenous input (ARMAX) process.
to the large sampling period that cannot be incorporated by the previous continuous-time approach. Unbiased estimates are guaranteed from the proposed prediction error method like discrete-
z (· )et) + alz (·-I )et) +
+ a._lz(l) et) + a.z{t)
= b,.,ui·-'\t - d , ) +
time prediction error methods. Also, it can identify the time delay term as well as the other system
+ b•.,u, (t - d , ) +
+b'.mU~·-') (t -d m )+
+k,e (· -' ){t-f.t)+
parameters simultaneously.
yet)
+b•.mu m(t -d .. )
+k. e{t-At)
=z{t) + e(t)
(4)
DELAYED
where, z{n)(t) denotes the n-th derivative of the
PROCESSES In this paper, we consider the following multi-input, single-output (MISO) process.
continuous-time signal z (t) . The noise e(t) is assumed as a discrete-time Gaussian noise with a
2.
CONTINUOUS-TIME
dx(t) = Ax(t)+Buv(t)+Ke(t-At) dt y(t)
=Cx(t) + e(t)
(1)
e(t)=e{t;), I; ~t
(2)
CONTINUOUS-TIME 3. PROPOSED PREDICTION ERROR METHOD Like discrete-time prediction error methods (Ljung (1987», we estimate the system parameters minimizing the prediction error. In this research,
where uv{t) and y (t) denote the m-dimensional delayed process input and the scalar process output, respectively. x(t) is the n-dimensional state. The
first, we identify the deterministic part ( A, B, D) by solving the output error minimization problem and second, the stochastic part (K) is identified.
problem in this research is to estimate the system matrices A,B , K and D=[d,
d2
sampling time At as follows.
d .. f
in
the uD{t) vector from the measured input-output 3.1 IDENTIFICATION OF DETERMINISTIC PART We estimate the deterministic part of the combined
data sets which are the following respective forms: uv(t)=[u,(t-d,) K
= [k.
C=[O
k. _,
0 0
k. _2
u .. (t-d .. W
deterministic-stochastic model by solving the following output error minimization.
k 2 kIf
o 1]
IJljll[V(A,B,D) = O.S i(y(t;)- Y
subject to
482
tN
i=1
dX(t) = Ax(t) + Bu D (t) dt
(7)
y(t) = Cx(t)
(8) (9)
where, y(t) and y(t) denote the process output and the model output, respectively. To solve this optimization problem, we use the following Newton method, which repeats (10) until the parameters converge within tolerance. 8 d (J) =
=
,n
= 1,2,
,m
,k
9 ( . _ I) _ [C!2V(9 d U d } P C!9~ () d
, k = 1,2,
[a a l
1»]-1[C!V(8C!9Aj - I»]
l •1
b
2•1
n. b b 1
l •m
2•m
(19)
wh..-,. G, =[00 b
(18)
,m
or
(10)
d
an b b
2
/ = 1,2,
10
ha> 1 fo, th, ,-
n.
m
(11)
k+ I-th component and 0 for others. u/(t-d/)
where, j denotes the iteration number and the rate constant is 0 < pSI. By choosing a large rate constant we can increase the convergence speed. But, if the parameters fluctuate or diverge the rate constant should be reduced.
denotes
The derivatives of the objective function with respect to the adjustable parameters can be calculated as follows: From (6), (12)-(16) are derived.
We need to consider the following equation to calculate (19).
the
/-th
delayed
uD(t)=[ul(t-dl )
=[b
Bk
n •k
process
input
of
um(t-dm)f
bn _ l.t
b 2•k
is
blkf
the
k-th
column vector of B.
C!Uk(t -d k ) = C!d k
duk(t -d k ) dt
(20)
Then, (19) can be rewritten like the following .
ay(t)
dy(t)
!!.-(ax(t) dt dd k
dbll
dbn •1
,k = 1,2,
"\(}2
U
da k
da k
dY.
(14)
,n
+~
'
k
=1,2,
,n /
=1,2,
(21)
,m
The second derivative of (6) is 2 2 a V(9 d ) =_~ f(y(t)_y'(t»a y(t;) (t -t_)
From (8) and (7), (14)-(16) and (17)-(19) are obtained. k = 1,2,
ax(t) _ B dU k (t - d k) k dt k
So, with a numerical derivative, we can approximately calculate dU k(t - d k) / dt and can solve (21).
(13)
dY(t) = C dX(t)
l= A ad
d
£...
I
I
.=1
d
"\92
U
f [ay(t; )][ay(t; )]T (t; _ t;_I) C!9 C!9
tN ;=1
,m (15)
tN
I
•
I
d
(22)
d
Here, the second summation becomes more dominant than the first summation when the solution is close to the true. Then, we can neglect the first term of the right hand side in (22).
db k ./ (16)
483
In summary, we calculate the prediction error y(t) - .9(t) for given Od (j -I) from (7) and (8) and the first and second derivatives of the objective function with respect to the adjustable parameters from (12)-(21) and (23). Then, we estimate the
updated parameters 8Aj) from (10). Repeat this procedure until converge.
Remarks: I . The time delay and the other system parameters can be estimated simultaneously through minimizing the prediction (output) error. 2. In many cases, we should remove spikes, outliers or no meaningful measurements due to electrical problems or some actions (replacement, cleaning etc.) for sensor maintenance. The proposed approach can incorporate directly this irregular sampling period due to measurement problems. Also, it can manipulate the small sampling time that cannot be treated efficiently by a discrete-time previous prediction error method for a constant covariance noise with respect to the sampling time. 3. Like discrete-time prediction error method, even though the actual output error is a colored noise that uncorrelated with the process input and the process output, the proposed strategy gives unbiased estimates. Consider the following actual process noises. dx(t) dt
= A.x(t) + Bu(t),
y(t)
3.2 IDENTIFICATION OF STOCHASTlC PART In a similar way, we estimate the stochastic part of the combined deterministic-stochastic model by solving the prediction error minimization for given deterministic part of A and B.
= Cx(t) + V n (t)
where vn (t) is the actual output noise. Then the expected objective function is like this.
- - =-L.(y(t; O.5f V(A,B) )- Y-(t ; ))2( t; -tH ) tN
where y d (t) denotes the actual deterministic part (Cx(t» of the process output. Here, note that the expected second term is constant and the expected third term is zero since the noises are uncorrelated with the process input. Therefore, minimizing the expected objective function of (6) is the same with minimizing the first term. So, the optimal estimates are unbiased regardless of the actual output noises when the process satisfies the moderate condition that a unique solution minimizing the first term is the actual parameters. 4. The continuous-time model identification method estimates the matrix A directly. On the other hand, we can say that the discrete-time model identification method calculates the matrix A indirectly after the step of estimating exp(-AL\t) (sampling period dependent matrix). The difference between the direct (continuoustime) identification and the indirect (discrete-time) in totally different identification results identification performances. When the sampling period is very small exp( - ALlt) converges to identity matrix. So, small uncertainties in estimating exp( - AM) make the reconstructed matrix A unstable or severely erroneous. Contrarily, the same uncertainty for the direct estimation of A might be negligible. For details, refer to Sung et al. (2000) . So, the continuous-time identification method is recommended for the small sampling period case.
;=}
mjn[V(K)
= 0.5
i
t N ~I
K
{e(t;)y (t; -
ti-J)]
(24)
subject to dX(t) dt
= Ax(t) + Bu(t) + Ke(t;_I)'
y(t) = Cx(t)
If y(t;) e(t; ) = 0
484
ti-J 5,
t < t; (25) (26)
is measured e(t;) = y(t;) - .9«(), else (27)
a y(t) / ao;
In the identification of the deterministic part, we
uncorrelated with
can choose any arbitrary t;. But, we restrict t; as a multiple of a small sampling time f:J (that is,
is close to the true. Then, we approximate the
when the solution
second derivative as follows .
= il1t,i = 0,1,2, »
in the stochastic part identification for mathematical simplicity. Using (27), we can incorporate unmeasured process outputs. We chose the most probable value as zero since we assume the noise e(t;) is a Gaussian noise. t;
2
(34)
In summary, we calculate the prediction error yet) - Ht) for given Bs U -1) from (25)-(26) and
To solve this optimization problem, we use the Newton method, which repeats (28) until the parameters converge within tolerance.
o ( ')=0 , ]
( ' -1)-
, ]
the derivatives of the objective function with respect to the adjustable parameters from (30)-(33) and (34). Then, we estimate the updated parameters
1 [a v(o,UP ao »]-I[aV(O,(J-I))](28) ao 2
,2
B,U)
,
converge.
(29)
where, j
denotes the iteration number and the
rate constant is 0 < p
~
Remarks:
1. The first and second
The initial values of the adjustable parameters can be determined by previous approaches. But, in this research, we recommend zeros as initial settings for simplicity because we believe that the proposed non linear optimization does not suffer much from local optima problem for MISO process.
derivatives of the objective function with respect to the adjustable parameters can be calculated as follows : From (24), (30)-(32) are obtained.
aV(oJ = _~ feet,) oy(t,) (t; aos t N .=1 aos
(30)
-tj-j)
ay~t)lT
4. SIMULATION RESULTS The following process is simulated to confirm the identification performance of the proposed method and to compare it with an existing discrete-time prediction error method (OE in Matlab).
(31)
ak. From (26) and (25), (32) and (33) are obtained.
ay~t) = C ax~t) , 1=1,2, ,n
akl
(32)
akl
If y(t;_I) is measured
1=1,2,
!!...(ax(.t))= A "aX(.t), dt
akl
ak,
where, G,
=[0 0
,nel se
1=1,2,
t'_1 ~ t < t.
10 I 0
r
, n (33)
has 1 for the n-/+ I-th
component and 0 for others. The
error
yet) - Y(t)
is
l
2
dt
dt
d z(t) d z(t) dz(t) (t ) -u - (t - 0 .5) + 3+ 3--+z l 2
(35)
yet) = z(t) + e(t)
(36)
dt
The process input is a Random-Binary-Sequence (RBS) with the minimum switching time of 0.5. The variance of e(tJ is 0.0001. The proposed and discrete-time prediction error method (oe function) in Matlab identified the third order model. We assumed that the actual time delay information is available when apply the discrete-time approach. But, we identified the time delay as well as other system parameters when apply the proposed method. We did 100 experiments for 0.01 sampling period and compared the worst cases in Figure l. The discrete-time approach with the exactly known
!!.-.(dX(t)]=A,dX(t) )_K"cdHt;-I) " , + G"( le tj-J " dt ak l dk l dk, ,t;_ I~t
from (28). Repeat this procedure until
approximately
485
time models", IEEE Transactions on Signal Processing 42887-897. Larimore, WE . (1990) "Canonical variate analysis in identification, filtering, and adaptive control" In Proc. 29'" IEEE Con! On Decision and Control, Honolulu, Hawaii, 596-604. Ljung L. (1987) System identification: Theory for the user. Englewood Cliffs, NJ: Prentice-Hall. Middleton, R. H. and Goodwin, G. C. (1990) Digital Control and Estimation Englewood Cliffs, NJ: Prentice-Hall. Sagara, S., and Zhao, Z. (1989) "Recursive identification of transfer function matrix in continuous systems via linear integral filter", Int. J. Control 50, 457-477. Soderstrom, T. and Stoica, P. (1989) System identification. Englewood Cliffs, NJ: Prentice-Hall. Soderstrom, T. and Mossberg, M. (2000) "Performance evaluation of methods for identifying continuous-time autoregressive processes", Automatica 36, 53-59. Sung, S.W, Lee, I. and Lee, 1. (1998) "New Process Identification Method for Automatic Design of PID Controllers", Automatica 34, 513520. Sung, S.W, Lee, S.Y., Kwak, HJ. and Lee, I. (2000) "Continuous-time Subspace System Identification Method", submitted to Ind. Eng. Chem. Res. Van Overschee, P. and De Moor, B. (1994) "N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems", Automatica 30, 75-93. Verhaegen, M. (1994) "Identification of the deterministic part of MIMO state space model given in innovations form from input-output data", Automatica 30, 1877-1883.
time delay shows very poor identification results for the small sampling time. Contrarily, the proposed method gives almost same plot with the real process. On the other hand, both the discretetime approach and the proposed method show very good identification results for the large sampling period. We confirmed from extensive simulations that the discrete-time approaches give acceptable accuracy for usual (large) sampling time but unacceptable accuracy for a small sampling time. The proposed approach gives good results for both sampling times.
--
Figure 1. Identified results (solid - plant, dotted proposed, dashed - discrete-time approaches) 5. CONCLUSIONS We proposed a continuous-time prediction error identification method to identify combined deterministic-stochastic continuous-time delayed processes. It has advantages compared with previous continuous-time and discrete-time identification methods since it can incorporate a large and a small sampling period test data as well as an irregular sampling period and the time delay in a systematic manner. Also, it is more attractive than any other previous works minimizing equation errors because more meaningful criterion for system identification should be based on the output error minimization. 6. REFERENCES Johansson, R. (1994) "Identification of continuous-
486
Copyright © IFAC Dynamics and Control of Process Systems, Jejudo Island, Korea, 200 I
D
0
[>
Publications www.elsevier.comllocate/ifac
CONTINUOUS-TIME PREDICTION ERROR IDENTIFICATION METHOD FOR DELAYED PROCESSES Su Whan Sung l " and In-Beum Lee 2 JDepartment
of Chemical Engineering, Korean Institute of Science and Technology 373-1 Kusong-dong, Yusong-gu, Taejon, 305-701, South Korea (*E-mail:[email protected]) 2Department of Chemical Engineering, Pohang University of Science and Technology San 31 Hyoja-dong, Pohang, 790-784, South Korea
Abstract: We propose a continuous-time prediction error method to identify combined deterministicstochastic continuous-time delayed processes. It minimizes the prediction error using the Newton optimization method with analytical derivatives of the objective function with respect to the adjustable parameters that contains the time delay. The proposed method can incorporate directly both small and large sampling time as well as irregular sampling time. Unbiased estimates are guaranteed like discretetime prediction error method. Copyright © 20011FAC Key words: time delay, continuous-time, prediction error identification method approaches under certain circumstances. So, many continuous-time system identification methods have been proposed like using block pulse functions, Laguerre polynomials, Poisson moment functionals. These functions produce initial value problem. Recently, many other methods of no initial value problem have been developed: linear integral transform or macro-difference expressions (Sagara and Zhao (1989» , complex variable transform (Johansson (1994» , delta operator (Middleton and Goodwin (1990» and ball-shape weighting integral transform (Sung et al. (1998». Several papers have been presented to analyze the of numerical derivatives and effects remove/compensate the bias in the estimation problem (Soderstrom and Mossberg (2000» . Recently, Sung et al. (2000) justified the usefulness of continuous-time identification methods against discrete-time ones in control point of view and proposed a continuous-time subspace system identification method.
1. INTRODUCTION
The main focus of this paper is to develop a prediction error identification method that can identify first, time delay second, system parameters in continuous-time domain third, incorporate a small/large and irregular sampling time. In practice, the time delay is much common in chemical industries due to transportation, measurement delay, and long transmission lines in pneumatic systems. Also, the time delay has been used widely to approximate the complex and high order dynamic system together with a low order process model like the first or second order one. But, most combined deterministic-stochastic system identification methods cannot treat the time delay in a systematic manner. On the other hand, most researchers in system identification field have concentrated their efforts on developing discrete-time identification methods like prediction error, instrumental variable or subspace identification methods (Ljung (1987), Soderstrom and Stoica (1989), Larimore (1990), Verhaegen (1994), Van Overschee and De Moor (1994». But, we had better use continuous-time
The above-mentioned previous methods have several disadvantages: Previous discrete-time approaches have problems that an irregular sampling period cannot be incorporated without major modification since the adjustable parameters are not constant for the irregular sampling period.
This work is supported by the BK21 Program. 481
Also, a small sampling period and/or a continuoustime input cannot be treated with acceptable accuracy. On the other hand, previous continuoustime approaches using integral transforms or delta operator cannot manipulate a large sampling period.
o -a. o -a._, o - a._2
000 o 0
A=
Furthermore, the time delay cannot be identified in a systematic manner.
o
o
o o
0 0
0 0
b•.,
-a 2
0
-a, b•.2
b•.m
b._,., b._,.
2
In this research, we propose a continuous-time prediction error method to identify combined
(3)
B = b._2., b._2.2
deterministic-stochastic continuous-time delayed processes. Compared with previous discrete-time identification methods, the proposed method does not suffer from the small sampling period or irregular sampling problem. Also, it can be applied
Here, (1) and (2) can be rewritten equivalently like the following continuous-time Auto-Regressive, Moving Average, eXogenous input (ARMAX) process.
to the large sampling period that cannot be incorporated by the previous continuous-time approach. Unbiased estimates are guaranteed from the proposed prediction error method like discrete-
z (· )et) + alz (·-I )et) +
+ a._lz(l) et) + a.z{t)
= b,.,ui·-'\t - d , ) +
time prediction error methods. Also, it can identify the time delay term as well as the other system
+ b•.,u, (t - d , ) +
+b'.mU~·-') (t -d m )+
+k,e (· -' ){t-f.t)+
parameters simultaneously.
yet)
+b•.mu m(t -d .. )
+k. e{t-At)
=z{t) + e(t)
(4)
DELAYED
where, z{n)(t) denotes the n-th derivative of the
PROCESSES In this paper, we consider the following multi-input, single-output (MISO) process.
continuous-time signal z (t) . The noise e(t) is assumed as a discrete-time Gaussian noise with a
2.
CONTINUOUS-TIME
dx(t) = Ax(t)+Buv(t)+Ke(t-At) dt y(t)
=Cx(t) + e(t)
(1)
e(t)=e{t;), I; ~t
(2)
CONTINUOUS-TIME 3. PROPOSED PREDICTION ERROR METHOD Like discrete-time prediction error methods (Ljung (1987», we estimate the system parameters minimizing the prediction error. In this research,
where uv{t) and y (t) denote the m-dimensional delayed process input and the scalar process output, respectively. x(t) is the n-dimensional state. The
first, we identify the deterministic part ( A, B, D) by solving the output error minimization problem and second, the stochastic part (K) is identified.
problem in this research is to estimate the system matrices A,B , K and D=[d,
d2
sampling time At as follows.
d .. f
in
the uD{t) vector from the measured input-output 3.1 IDENTIFICATION OF DETERMINISTIC PART We estimate the deterministic part of the combined
data sets which are the following respective forms: uv(t)=[u,(t-d,) K
= [k.
C=[O
k. _,
0 0
k. _2
u .. (t-d .. W
deterministic-stochastic model by solving the following output error minimization.
k 2 kIf
o 1]
IJljll[V(A,B,D) = O.S i(y(t;)- Y
subject to
482
tN
i=1
dX(t) = Ax(t) + Bu D (t) dt
(7)
y(t) = Cx(t)
(8) (9)
where, y(t) and y(t) denote the process output and the model output, respectively. To solve this optimization problem, we use the following Newton method, which repeats (10) until the parameters converge within tolerance. 8 d (J) =
=
,n
= 1,2,
,m
,k
9 ( . _ I) _ [C!2V(9 d U d } P C!9~ () d
, k = 1,2,
[a a l
1»]-1[C!V(8C!9Aj - I»]
l •1
b
2•1
n. b b 1
l •m
2•m
(19)
wh..-,. G, =[00 b
(18)
,m
or
(10)
d
an b b
2
/ = 1,2,
10
ha> 1 fo, th, ,-
n.
m
(11)
k+ I-th component and 0 for others. u/(t-d/)
where, j denotes the iteration number and the rate constant is 0 < pSI. By choosing a large rate constant we can increase the convergence speed. But, if the parameters fluctuate or diverge the rate constant should be reduced.
denotes
The derivatives of the objective function with respect to the adjustable parameters can be calculated as follows: From (6), (12)-(16) are derived.
We need to consider the following equation to calculate (19).
the
/-th
delayed
uD(t)=[ul(t-dl )
=[b
Bk
n •k
process
input
of
um(t-dm)f
bn _ l.t
b 2•k
is
blkf
the
k-th
column vector of B.
C!Uk(t -d k ) = C!d k
duk(t -d k ) dt
(20)
Then, (19) can be rewritten like the following .
ay(t)
dy(t)
!!.-(ax(t) dt dd k
dbll
dbn •1
,k = 1,2,
"\(}2
U
da k
da k
dY.
(14)
,n
+~
'
k
=1,2,
,n /
=1,2,
(21)
,m
The second derivative of (6) is 2 2 a V(9 d ) =_~ f(y(t)_y'(t»a y(t;) (t -t_)
From (8) and (7), (14)-(16) and (17)-(19) are obtained. k = 1,2,
ax(t) _ B dU k (t - d k) k dt k
So, with a numerical derivative, we can approximately calculate dU k(t - d k) / dt and can solve (21).
(13)
dY(t) = C dX(t)
l= A ad
d
£...
I
I
.=1
d
"\92
U
f [ay(t; )][ay(t; )]T (t; _ t;_I) C!9 C!9
tN ;=1
,m (15)
tN
I
•
I
d
(22)
d
Here, the second summation becomes more dominant than the first summation when the solution is close to the true. Then, we can neglect the first term of the right hand side in (22).
db k ./ (16)
483
In summary, we calculate the prediction error y(t) - .9(t) for given Od (j -I) from (7) and (8) and the first and second derivatives of the objective function with respect to the adjustable parameters from (12)-(21) and (23). Then, we estimate the
updated parameters 8Aj) from (10). Repeat this procedure until converge.
Remarks: I . The time delay and the other system parameters can be estimated simultaneously through minimizing the prediction (output) error. 2. In many cases, we should remove spikes, outliers or no meaningful measurements due to electrical problems or some actions (replacement, cleaning etc.) for sensor maintenance. The proposed approach can incorporate directly this irregular sampling period due to measurement problems. Also, it can manipulate the small sampling time that cannot be treated efficiently by a discrete-time previous prediction error method for a constant covariance noise with respect to the sampling time. 3. Like discrete-time prediction error method, even though the actual output error is a colored noise that uncorrelated with the process input and the process output, the proposed strategy gives unbiased estimates. Consider the following actual process noises. dx(t) dt
= A.x(t) + Bu(t),
y(t)
3.2 IDENTIFICATION OF STOCHASTlC PART In a similar way, we estimate the stochastic part of the combined deterministic-stochastic model by solving the prediction error minimization for given deterministic part of A and B.
= Cx(t) + V n (t)
where vn (t) is the actual output noise. Then the expected objective function is like this.
- - =-L.(y(t; O.5f V(A,B) )- Y-(t ; ))2( t; -tH ) tN
where y d (t) denotes the actual deterministic part (Cx(t» of the process output. Here, note that the expected second term is constant and the expected third term is zero since the noises are uncorrelated with the process input. Therefore, minimizing the expected objective function of (6) is the same with minimizing the first term. So, the optimal estimates are unbiased regardless of the actual output noises when the process satisfies the moderate condition that a unique solution minimizing the first term is the actual parameters. 4. The continuous-time model identification method estimates the matrix A directly. On the other hand, we can say that the discrete-time model identification method calculates the matrix A indirectly after the step of estimating exp(-AL\t) (sampling period dependent matrix). The difference between the direct (continuoustime) identification and the indirect (discrete-time) in totally different identification results identification performances. When the sampling period is very small exp( - ALlt) converges to identity matrix. So, small uncertainties in estimating exp( - AM) make the reconstructed matrix A unstable or severely erroneous. Contrarily, the same uncertainty for the direct estimation of A might be negligible. For details, refer to Sung et al. (2000) . So, the continuous-time identification method is recommended for the small sampling period case.
;=}
mjn[V(K)
= 0.5
i
t N ~I
K
{e(t;)y (t; -
ti-J)]
(24)
subject to dX(t) dt
= Ax(t) + Bu(t) + Ke(t;_I)'
y(t) = Cx(t)
If y(t;) e(t; ) = 0
484
ti-J 5,
t < t; (25) (26)
is measured e(t;) = y(t;) - .9«(), else (27)
a y(t) / ao;
In the identification of the deterministic part, we
uncorrelated with
can choose any arbitrary t;. But, we restrict t; as a multiple of a small sampling time f:J (that is,
is close to the true. Then, we approximate the
when the solution
second derivative as follows .
= il1t,i = 0,1,2, »
in the stochastic part identification for mathematical simplicity. Using (27), we can incorporate unmeasured process outputs. We chose the most probable value as zero since we assume the noise e(t;) is a Gaussian noise. t;
2
(34)
In summary, we calculate the prediction error yet) - Ht) for given Bs U -1) from (25)-(26) and
To solve this optimization problem, we use the Newton method, which repeats (28) until the parameters converge within tolerance.
o ( ')=0 , ]
( ' -1)-
, ]
the derivatives of the objective function with respect to the adjustable parameters from (30)-(33) and (34). Then, we estimate the updated parameters
1 [a v(o,UP ao »]-I[aV(O,(J-I))](28) ao 2
,2
B,U)
,
converge.
(29)
where, j
denotes the iteration number and the
rate constant is 0 < p
~
Remarks:
1. The first and second
The initial values of the adjustable parameters can be determined by previous approaches. But, in this research, we recommend zeros as initial settings for simplicity because we believe that the proposed non linear optimization does not suffer much from local optima problem for MISO process.
derivatives of the objective function with respect to the adjustable parameters can be calculated as follows : From (24), (30)-(32) are obtained.
aV(oJ = _~ feet,) oy(t,) (t; aos t N .=1 aos
(30)
-tj-j)
ay~t)lT
4. SIMULATION RESULTS The following process is simulated to confirm the identification performance of the proposed method and to compare it with an existing discrete-time prediction error method (OE in Matlab).
(31)
ak. From (26) and (25), (32) and (33) are obtained.
ay~t) = C ax~t) , 1=1,2, ,n
akl
(32)
akl
If y(t;_I) is measured
1=1,2,
!!...(ax(.t))= A "aX(.t), dt
akl
ak,
where, G,
=[0 0
,nel se
1=1,2,
t'_1 ~ t < t.
10 I 0
r
, n (33)
has 1 for the n-/+ I-th
component and 0 for others. The
error
yet) - Y(t)
is
l
2
dt
dt
d z(t) d z(t) dz(t) (t ) -u - (t - 0 .5) + 3+ 3--+z l 2
(35)
yet) = z(t) + e(t)
(36)
dt
The process input is a Random-Binary-Sequence (RBS) with the minimum switching time of 0.5. The variance of e(tJ is 0.0001. The proposed and discrete-time prediction error method (oe function) in Matlab identified the third order model. We assumed that the actual time delay information is available when apply the discrete-time approach. But, we identified the time delay as well as other system parameters when apply the proposed method. We did 100 experiments for 0.01 sampling period and compared the worst cases in Figure l. The discrete-time approach with the exactly known
!!.-.(dX(t)]=A,dX(t) )_K"cdHt;-I) " , + G"( le tj-J " dt ak l dk l dk, ,t;_ I~t
from (28). Repeat this procedure until
approximately
485
time models", IEEE Transactions on Signal Processing 42887-897. Larimore, WE . (1990) "Canonical variate analysis in identification, filtering, and adaptive control" In Proc. 29'" IEEE Con! On Decision and Control, Honolulu, Hawaii, 596-604. Ljung L. (1987) System identification: Theory for the user. Englewood Cliffs, NJ: Prentice-Hall. Middleton, R. H. and Goodwin, G. C. (1990) Digital Control and Estimation Englewood Cliffs, NJ: Prentice-Hall. Sagara, S., and Zhao, Z. (1989) "Recursive identification of transfer function matrix in continuous systems via linear integral filter", Int. J. Control 50, 457-477. Soderstrom, T. and Stoica, P. (1989) System identification. Englewood Cliffs, NJ: Prentice-Hall. Soderstrom, T. and Mossberg, M. (2000) "Performance evaluation of methods for identifying continuous-time autoregressive processes", Automatica 36, 53-59. Sung, S.W, Lee, I. and Lee, 1. (1998) "New Process Identification Method for Automatic Design of PID Controllers", Automatica 34, 513520. Sung, S.W, Lee, S.Y., Kwak, HJ. and Lee, I. (2000) "Continuous-time Subspace System Identification Method", submitted to Ind. Eng. Chem. Res. Van Overschee, P. and De Moor, B. (1994) "N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems", Automatica 30, 75-93. Verhaegen, M. (1994) "Identification of the deterministic part of MIMO state space model given in innovations form from input-output data", Automatica 30, 1877-1883.
time delay shows very poor identification results for the small sampling time. Contrarily, the proposed method gives almost same plot with the real process. On the other hand, both the discretetime approach and the proposed method show very good identification results for the large sampling period. We confirmed from extensive simulations that the discrete-time approaches give acceptable accuracy for usual (large) sampling time but unacceptable accuracy for a small sampling time. The proposed approach gives good results for both sampling times.
--
Figure 1. Identified results (solid - plant, dotted proposed, dashed - discrete-time approaches) 5. CONCLUSIONS We proposed a continuous-time prediction error identification method to identify combined deterministic-stochastic continuous-time delayed processes. It has advantages compared with previous continuous-time and discrete-time identification methods since it can incorporate a large and a small sampling period test data as well as an irregular sampling period and the time delay in a systematic manner. Also, it is more attractive than any other previous works minimizing equation errors because more meaningful criterion for system identification should be based on the output error minimization. 6. REFERENCES Johansson, R. (1994) "Identification of continuous-
486