Ind. Eng. Chem. Res. 2001, 40, 5743-5751
5743
Prediction Error Identification Method for Continuous-Time Processes with Time Delay Su Whan Sung* Department of Chemical Engineering and Center for Ultramicrochemical Process Systems, Korea Advanced Institute of Science and Technology, 373-1 Kusong-dong, Yusong-Gu, Taejon 305-701, Korea
In-Beum Lee Department of Chemical Engineering, Pohang University of Science and Technology, San 31 Hyoja Dong, Pohang 790-784, Korea
We propose a continuous-time prediction error identification method to identify combined deterministic-stochastic continuous-time processes with time delay. It minimizes the prediction error using the Levenberg-Marquardt optimization method with exact derivatives of the objective function with respect to the adjustable parameters that include the time delay. Compared with previous discrete-time identification methods, the proposed method does not suffer from a small sampling time problem. Also, while previous continuous-time approaches using transforms cannot treat a large sampling time, the proposed method can incorporate directly both small and large sampling times as well as irregular sampling time. It can determine the time delay systematically; meanwhile, previous methods use ad hoc approaches. We derive the error covariance matrix and justify the small sampling problem of discrete-time identification methods. 1. Introduction In practice, the time delay is very common in chemical industries because of transportation. Also, time delay plays an important role in describing high-order dynamics with a low-order model like the first (second) order plus time delay model. However, most prediction error methods to identify combined deterministicstochastic processes cannot treat the time delay systematically. The goal of this research is to identify combined deterministic-stochastic continuous-time processes with time delay through minimization of the prediction error. Until now, most researchers in the system identification field have concentrated their efforts on developing discrete-time identification methods such as prediction error, instrumental variable, or subspace identification methods.1-6 We guess several reasons for the incomparable popularity of discrete-time identification over continuous-time identification. First, the mathematics for discrete time is more easily and better established. The stochastic theory of discrete-time approaches is very simple and clear compared with that of continuous-time approaches. We need to calculate time derivatives in continuous-time approaches. This makes continuoustime approaches complicated. Second, estimating the physical parameters of the continuous-time physical model is not important if the model is to design controllers. For these reasons, discrete-time approaches have been recognized as more efficient than continuoustime approaches. Certainly, the above-mentioned arguments against continuous-time approaches are convincible in many cases. However, we should use continuous-time approaches under certain circumstances. First, when we need to estimate physical parameters, continuous-time * To whom all correspondence should be addressed. Tel: +82-42-869-3992. Fax: +82-42-869-3910. E-mail: swsung@ mail.kaist.ac.kr.
approaches are definitely preferred. Second, when the process input is a physical signal like temperature, pressure, level, etc., we should use continuous-time approaches because the physical signal is always a continuous-time signal rather than a discrete-time signal from a zero-order holder. Third, when the sampling time is small, continuous-time approaches are recommended. Someone might say “it is also possible to use discrete-time approaches by reducing the sampling time for the continuous-time input case”. However, this is not recommendable because discrete-time approaches suffer from numerically ill-conditioning for a small sampling time and the frequency modeling errors increase inversely proportional to the sampling time for a certain type of noise (which will be discussed later).7,8 Many continuous-time system identification methods have been proposed such as using Walsh functions,9 block pulse functions,10 Legendre or Laguerre polynomials,11,12 multiple integrators,13,14 and Poisson moment functionals.15 These functions produce an initial value problem. Recently, many other methods free from an initial value problem have been developed: linear integral transform or macrodifference expressions,16,17 digital filtering approaches,18 a complex variable transform,19,20 a δ operator,8,21 and a bell-shaped weighting integral transform.22 Several papers have been presented to analyze the effects of numerical derivatives and to remove/compensate the bias in the estimation problem.23-25 Several papers among them have justified the necessity of the continuous-time system identification from the viewpoint of estimating the physical parameters of the physical continuous-time linear system. Recently, the usefulness of continuous-time identification methods against discrete-time ones from the control point of view was justified.7 The above-mentioned previous methods have several disadvantages: Previous discrete-time approaches have problems that an irregular sampling time cannot be incorporated without major modification because the
10.1021/ie0100636 CCC: $20.00 © 2001 American Chemical Society Published on Web 11/02/2001
5744
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001
adjustable parameters are not constant for the irregular sampling time. Also, a small sampling time and/or a continuous-time input cannot be treated with acceptable accuracy (which will be exemplified later). On the other hand, previous continuous-time approaches using integral transforms or a δ operator cannot manipulate a large sampling time because they use numerical integrations or numerical derivatives of the process output to realize the transforms. The most recent version of the system identification MATLAB toolbox can identify the parametrized continuous-time state-space model. An oelmdf function of the CONTSID MATLAB toolbox can identify the continuous-time output error model. The functions of the two packages for the continuous-time process can overcome the small sampling time problem of the previous discrete-time approaches and the large sampling time of the previous continuous-time approaches. However,, the user should specify the time delay with ad hoc methods of inspecting the initial impulse responses or exhaustive search. Also, the continuous-time identification functions of the system identification MATLAB toolbox use numerical differentiation to calculate the derivatives of the cost function, and the continuous-time modeling is for a regular sampling time. The oelmdf function of the CONTSID toolbox can incorporate the irregular sampling time, and it uses exact differential information in the optimization step. However, it still cannot identify the time delay, and it is only for the deterministic model. In this research, we propose a continuous-time prediction error method to identify combined deterministicstochastic continuous-time processes with time delay. It identifies the system parameters by minimizing the prediction error using the Levenberg-Marquardt optimization method with exact derivatives of the objective function with respect to the adjustable parameters. Compared with previous discrete-time identification methods, the proposed method does not suffer from the small sampling time or irregular sampling problem. Also, it can be applied to the large sampling time that cannot be incorporated by the previous continuous-time approach using transforms. Also, it can identify the time delay term as well as the other system parameters simultaneously in a systematic way. We derive error covariance matrices for the proposed identification method and justify the small sampling time problem of discrete-time identification methods. This paper is organized as follows: Section 2 defines the problem in this research. In section 3, we develop the proposed continuous-time prediction error method and remark on the properties of the proposed method. Section 4 gives numerical examples to confirm the identification performance of the proposed method and to compare it with the existing discrete-time prediction error method. Section 5 summarizes the main results. 2. Continuous-Time Processes with Time Delay In this paper, we consider the following multi-input, single-output (MISO) process.
dx(t) ) Ax(t) + Budelay(t) + Ke(t) dt y(t) ) Cx(t) + e(t)
(1) (2)
where udelay(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 problem in this research is to estimate the system matrices A, B, K, and D ) [d1 d2 ... dm]T in the udelay(t) vector from the discrete-sampled-output data and the given process input.
udelay(t) )
[
]
[u1(t-d1) u2(t-d2) ... um-1(t-dm-1) um(t-dm)]T (3) 0 1 0 A) l 0 0
[
0 0 1 l 0 0
bn,1 bn-1,1 b B ) n-2,1 l b2,1 b1,1
0 0 0 l 0 0
‚‚‚ ‚‚‚ ‚‚‚ l ‚‚‚ ‚‚‚
bn,2 bn-1,2 bn-2,2 l b2,2 b1,2
0 0 0 l 0 1
-an -an-1 -an-2 l -a2 -a1
‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚
bn,m bn-1,m bn-2,m l b2,m b1,m
K ) [kn kn-1 kn-2 ... k2 k1]T
(4)
]
C ) [0 0 0 ... 0 1]
(5)
(6) (7)
Here, (1) and (2) can be rewritten equivalently like the following continuous-time autoregressive, moving average, exogenous input (ARMAX) process.
z(n)(t) + a1z(n-1)(t) + ... + an-1z(1)(t) + anz(t) ) b1,1u1(n-1)(t-d1) + ‚‚‚ + bn,1u1(t-d1) + ... + b1,mum(n-1)(t-dm) + ... + bn,mum(t-dm) + k1e(n-1)(t) + ... + kne(t) y(t) ) z(t) + e(t) (8) where z(n)(t) denotes the nth derivative of the continuous-time signal z(t). The noise e(t) is assumed as a piecewise constant Gaussian noise with a time length ∆t as follows.
e(t) ) e(ti), ti e t < ti + ∆t E{e2(ti)} ) σ2
(9)
Remarks. 1. The chosen process of (1) and (2) in this research is a hybrid model structure combining the continuoustime deterministic part and the stochastic part represented by a piecewise continuous-time white noise. A similar model is treated in a recent paper.26 They identify the continuous-time model for a regular sampling time in the frequency domain. It does not consider the estimation problem of the time delay. 2. (1) and (2) can be rewritten as follows:
x(k) ) A∆tx(k-1) + B∆tudelay(k - 1) + K∆te(k-1) (10) y(k) ) Cx(k) + e(k) where A∆t ) exp(A∆t), B∆t )
∫∆t 0 exp(Aτ)B
(11) dτ, and K∆t )
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5745
∫∆t 0 exp(Aτ)K dτ and k is sampling number for the sampling time of ∆t. As shown in (10) and (11), the stochastic model chosen in this research is actually the standard discrete-time stochastic model. 3. The proposed method is to identify the parameters (A, B, and K) of the differential equation model of (1) and (2), while discrete-time identification methods are to identify the parameters (A∆t, B∆t, and K∆t) of the difference equation model of (10) and (11). This difference results in different performances in terms of the sampling time (which will be discussed later). 4. Some parameters among the adjustable parameters of (3)-(6) can be fixed when we apply the proposed identification method. For example, if we know there are no zeros, then we can fix all elements of B as zeros except the first-row vector, or if there is definitely no time delay, we can fix the time delay of the above model as zero. 3. Proposed Continuous-Time Prediction Error Method Like discrete-time prediction error methods,1 we estimate the system matrices by minimizing the prediction error. In this research, first, we identify the deterministic part (A, B, and D) by solving the output error minimization problem and, second, the stochastic part (K) is identified. This two-step approach is surely suboptimal compared with a simultaneous identification. However, we find the usefulness of the two-step approach when considering practical aspects: (1) In industrial chemical plants, the time length of the intentional perturbation for the identification should be as short as possible. That is, the data length is not enough to identify both the deterministic and stochastic parts. So, it is desirable that we use the intentionally perturbed data to identify only the deterministic part and additionally use much longer historical data sets to identify the stochastic part. (2) A simultaneous approach would identify A - KC of the closed-loop observer rather than A (or K) itself, and then the accuracy of A (or K) itself can be poor or even unstable when there are big structural plant/model mismatches or nonlinearities. So, we develop the proposed method on the basis of the two-step approach. However, if the simultaneous identification should be done for some reason, it can be developed easily by combining the two separated steps. 3.1. Identification of the Deterministic Part. We estimate the deterministic part of the combined deterministic-stochastic model from the discrete-sampled data y(ti), i ) 1, ..., N, and known u(t) by solving the following output error minimization.
[
min V(A ˆ ,B ˆ ,D ˆ) )
A ˆ ,B ˆ ,D ˆ
0.5 N
∑(y(ti) - yˆ (ti))2
N i)1
]
(12)
model output, respectively. To solve this optimization problem, we use the following Levenberg-Marquardt method like Ljung,1 which repeats (16) until the parameters converge within tolerance.
θˆ d(j) ) θˆ d(j-1) -
[
∂2V(θˆ d(j-1)) ∂θˆ d2
][ -1
+ RI
]
∂V(θˆ d(j-1)) ∂θˆ d
(16)
θˆ d ) [aˆ 1aˆ 2 ... aˆ nbˆ 1,1bˆ 2,1 ... bˆ n,1 ... bˆ 1,mbˆ 2,m ... bˆ n,m, dˆ 1, dˆ 2, ..., dˆ m]T (17)
where j denotes the iteration number and R is a small positive value that can be updated every iteration to compromise between the robustness and the convergence rate. The derivatives of the objective function with respect to the adjustable parameters can be calculated as follows: From (12), (18) is derived.
∂V(θˆ d) ∂θˆ d
)-
1
N
∑(y(ti) - yˆ (ti)) Ni)1
∂yˆ (ti) ∂θˆ d
∂yˆ (t) ) ∂θˆ d
[
(18)
]
∂yˆ (t) ∂yˆ (t) ∂yˆ (t) T ∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ... ... ... ... ... ... ∂aˆ 1 ∂aˆ n ∂bˆ 1,1 ∂bˆ n,1 ∂bˆ 1,m ∂bˆ n,m ∂dˆ 1 ∂dˆ m (19)
From (13) and (14), (20)-(23) are obtained.
∂yˆ (t) ∂xˆ (t) )C ∂θˆ d ∂θˆ d
( )
∂xˆ (t) d ∂xˆ (t) )A ˆ - Gkyˆ (t), k ) 1, 2, ..., n dt ∂aˆ k ∂aˆ k
(20)
(21)
( )
∂xˆ (t) d ∂xˆ (t) )A ˆ + Gkul(t-dl), k ) 1, 2, ..., n, dt ∂bˆ k,l ∂bˆ k,l l ) 1, 2, ..., m (22)
( )
∂uk(t-dk) ∂xˆ (t) d ∂xˆ (t) )A ˆ + Bk , k ) 1, 2, ..., m dt ∂dˆ k ∂dˆ k ∂dk (23) where
subject to
dxˆ (t) )A ˆ xˆ (t) + B ˆ udelay(t) dt
(13)
yˆ (t) ) Cxˆ (t)
(14)
t0 ) 0 < t1 < ... < tN-1 < tN
(15)
where y(t) and yˆ (t) denote the process output and the
has 1 for the (n - k + 1)th component and 0 for others. ul(t-dl) denotes the lth input of udelay(t) ) [u1(t-d1) ... um(t-dm)]T. Bk ) [bn,k bn-1,k ... b2,k b1,k]T is the kth column vector of B. The initial values of ∂xˆ (t)/∂θˆ d|t)0 are zero because the parameter of θˆ d does not affect the
5746
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001
initial values of the states. We need to note the following equation to calculate (23).
∂uk(t-dk) duk(t-dk) )∂dk dt
(24)
where vn(t) is the actual output noise. Then the expected objective function is like this.
V(A ˆ ,B ˆ) ) )
Then, (23) can be rewritten like the following.
0.5 N tN
0.5 N tN
( )
duk(t-dk) ∂xˆ (t) d ∂xˆ (t) )A ˆ - Bk , k ) 1, 2, ..., m dt ∂dˆ k ∂dˆ k dt (25)
)
∂θˆ d2
)-
1
N
∑
Ni)1
(y(ti) - yˆ (ti))
∂2yˆ (ti)
tN
1
∑
0.5 N
[ ][ ] ∂yˆ (ti) ∂yˆ (ti) ∂θˆ d
T
(26)
When the solution is close to true, we can neglect the first term of the right-hand side in (26).
∂θˆ d2
≈
1
N
∑
[ ][ ] ∂yˆ (ti) ∂yˆ (ti)
Ni)1 ∂θˆ d
tN
(vn(ti))2(ti - ti-1) + ∑ i)1
∂θˆ d
N
1 tN
Ni)1 ∂θˆ d
∂2V(θˆ d)
(yd(ti) - yˆ (ti))2(ti - ti-1) + ∑ i)1
+
∂θˆ d2 N
(yd(ti) + vn(ti) - yˆ (ti))2(ti - ti-1) ∑ i)1
0.5 N
With a numerical derivative, we can calculate duk(tdk)/dt from the given process input and can solve (25). The second derivative of (12) is
∂2V(θd)
(y(ti) - yˆ (ti))2(ti - ti-1) ∑ i)1
T
(27)
In summary, we calculate xˆ (ti) and ∂xˆ (t)/∂θˆ |t)ti for given θˆ d(j-1) by selecting them at every ti while continuously solving the ordinary differential equations of (13), (21), (22), and (25) simultaneously. Next, it is straightforward to calculate the updated parameters θˆ d(j) from (16). Repeat this procedure until the parameters converge. Remarks. 1. The time delay and the other system parameters can be estimated simultaneously through minimization of the prediction (output) error. Also, we can fix some parameters of (13) as mentioned in section 2. In this case, we can simply delete corresponding equations from the above-derived equations to incorporate the modification. 2. In many cases, we should remove spikes, outliers, or nonmeaningful measurements due to electrical problems or some actions (replacement, cleaning, etc.) for sensor maintenance. The proposed approach can incorporate directly this irregular sampling time due to measurement problems. Look at the objective function of (12). The discrete measurement data can be sampled at arbitrary (irregular) instances. Also, it can manipulate the small sampling time that cannot be treated efficiently by previous discrete-time methods (which will be discussed later). 3. Like discrete-time prediction error methods, even though the actual output error is a colored noise that is uncorrelated with the process input, the proposed strategy gives unbiased estimates. Consider the following actual process noise.
vn(ti)(yd(ti) - yˆ (ti))(ti - ti-1) ∑ i)1
where yd(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 because the noise is uncorrelated with the process input. Therefore, minimizing the expected objective function of (12) is the same as minimizing the first term. So, the optimal estimates are unbiased regardless of the actual output noise when the process satisfies the moderate condition that a unique solution minimizing the first term is the actual parameters. 4. The initial values of the adjustable parameters can be guessed from physical insight or determined roughly by previous approaches. For example, after identifying the continuous-time high-order ARX model using previous identification methods, we can obtain a low order plus time delay model with a model reduction procedure.27 The estimated model can be used as the initial parameter for the proposed method. 5. There are various techniques to determine the model structure such as considering physical insight, examining the spectral analysis estimate, testing ranks in covariance matrices, examining the information matrix, minimizing Akaike’s information theoretic criterion (AIC), and checking the whiteness of residuals. For details, refer to work by Ljung.1 These techniques can prevent overparametrization. 3.2. Expected Error Covariance Matrix. Consider the following Taylor series around the real parameter of θd with the assumption that N is sufficiently large.1
∂V(θˆ d,N) ∂V(θd,N) ∂2V(θd) ) + (θˆ d - θd) ∂θˆ d ∂θd ∂θ 2
(28)
d
Because the proposed method estimates θˆ d satisfying ∂V(θˆ d,N)/∂θˆ d ) 0, the following is obtained from (28).
(θˆ d - θd)(θˆ d - θd)T )
( )( ∂2V(θd)
-1
∂θd2
)(
)( ) T
∂V(θd,N) ∂V(θd,N) ∂θd ∂θd
∂2V(θd) ∂θd2
-T
(29)
where
dx(t) ) Ax(t) + Bu(t) dt
∂V(θd,N)
y(t) ) Cx(t) + vn(t)
∂θd
)-
1
N
∂yˆ (ti,θd)
∑vo(ti,θd) Ni)1
∂θd
(30)
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5747
∂2V(θd)
)
1
N
∑
Ni)1
∂θd2
[ ][ ] ∂yˆ (ti,θd) ∂yˆ (tj,θd) ∂θd
T
(31)
∂θd
vo(ti,θd) ) y(ti) - yˆ (ti,θd)
(32)
(29) is the error covariance matrix because the proposed method gives unbiased estimates. Here, the expected covariance matrix of (30) is as follows.
E
[(
)]
)(
∂V(θd,N) ∂V(θd,N) ∂θd
1 1
N N
∑∑
N Ni)1j)1
T
) ∂θd ∂yˆ (ti,θd) ∂yˆ (tj,θd)
( )( ) ∂θd
T
∂θd
E[vo(ti,θd) vo(tj,θd)] (33)
It should be noted that NE[(∂V(θd,N)/∂θd)(∂V(θd,N)/∂θd)T] becomes a constant matrix for a sufficiently large N, and let us call the matrix Q as follows:
1
N N
∑∑
Ni)1j)1
( )( ) ∂yˆ (ti,θd) ∂yˆ (tj,θd) ∂θd
∂θd
E[vo(ti,θd) vo(tj,θd)] ) Q (34)
( ) ( ) -1
∂2V(θd)
Q
-T
∂θd2
(35)
If the output noise vo(ti) is a white noise of constant variance σ2, (34) and (35) can be simplified as follows: N
∑
Ni)1
( )( ) ∂yˆ (ti,θd) ∂yˆ (ti,θd) ∂θd
σ2 ) Q
(36)
∂θd
E[(θˆ d - θd)(θˆ d - θd) ] )
( ∑( )( ) ) N
N Ni)1
∂yˆ (ti,θd) ∂yˆ (ti,θd) ∂θd
2
σ (37)
∂θd
E[(θˆ ∆t - θ∆t)(θˆ ∆t - θ∆t)T] ) 1 1 N ∂yˆ (ti,θ∆t) ∂yˆ (ti,θ∆t) N Ni)1
∂θ∆t
)(
∂θ∆t
∂θ∆t
))
σ2 (39)
Here, it should be noted that there is 1/∆t2 that means the expected modeling error of discrete-time identification methods may increase as the sampling time decreases. For example, if we fix the identification time length as a constant value of tf and use a regular sampling time then, N ) tf/∆t. Then, (37) for the continuous-time case becomes
E[(θˆ d - θd)(θˆ d - θd)T] )
( ∑( )( ) )
∆t 1
∂yˆ (ti,θd) ∂yˆ (ti,θd)
N
tf Ni)1
∂θd
T -1
σ2 (40)
∂θd
∂θ∆t
))
E[(θˆ d - θd)(θˆ d - θd)T] ) 1 1 N ∂yˆ ∆t(ti,θ∆t) ∂yˆ ∆t(ti,θ∆t)
( ∑(
tf∆t Ni)1
T -1
Up to here, we estimated the expected error covariance matrix of the proposed continuous-time prediction error method. In a similar way, we can derive the following expected error covariance matrix of discrete-time identification methods1 if we assume that the process output is measured with a regular sampling time of ∆t and the output error is a white noise:
( ∑(
)(
T -1
and (39) for the discrete-time case becomes
T
T
1 1
( ∑(
∆t N Ni)1
T
2 1 ∂ V(θd) E[(θˆ d - θd)(θˆ d - θd) ] ) N ∂θ 2 d
1
E[(θˆ d - θd)(θˆ d - θd)T] ) 1 1 N ∂yˆ ∆t(ti,θ∆t) ∂yˆ (ti,θ∆t) 2
Then, the expected error covariance matrix of (29) becomes T
are A∆t ) exp(A∆t) and B∆t ) ∫∆t 0 exp(Aτ) dτ B, while those identified by continuous-time methods are A and B. That is, the parameters themselves are totally different between the two approaches. This results in totally different behaviors for a small sampling time.7 If the sampling time is small, the modeling error of ˆ ∆t ) discrete-time identification methods, i.e., A∆t - A ˆ ∆t ) ∫∆t exp(Aτ) dτ B exp(A∆t) - exp(A ˆ ∆t) and B∆t - B 0 exp(A ˆ τ) dτ B ˆ from (10) and (11) can be ap- ∫∆t 0 ˆ ∆t ) (A - A ˆ )∆t and B∆t - B ˆ ∆t ) proximated to A∆t - A (B - B ˆ )∆t. This means that the system parameters of discrete-time models are those of the continuous-time model multiplied by the sampling time ∆t, that is, θ∆t - θˆ ∆t ) (θd - θˆ d)∆t if the sampling time is small enough and the same state-space realization method is used. Then, (38) can be rewritten as follows:
T -1
σ2 (38)
where θ∆t and θˆ ∆t are the real parameters of discretetime processes and the estimates of discrete-time identification methods. From (37) and (38), we can recognize that both give the same decreasing rate of 1/N. So, we can say that, from the accuracy point of view, the two approaches of continuous-time and discrete-time methods do not show remarkable differences. However, the system parameters identified by discrete-time methods
∂θ∆t
)(
∂θ∆t
))
T -1
σ2 (41)
Now, from the comparison of (40) and (41), we see a remarkable difference between the continuous-time case and the discrete-time case: as the sampling time is decreased, the modeling error of the continuous-time case decreases; contrarily, that of the discrete case increases. This is the reason why discrete-time identification methods suffer from the small sampling time problem and why continuous-time identification methods give better accuracy as the sampling time is decreased. In summary, 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 estimating exp(A∆t) (sampling time dependent matrix). The difference between the direct (continuous-time) identification and the indirect (discretetime) identification results in totally different identification performances even though the parameter convergence rates with respect to the data number are the same as those shown in (37) and (38). As an example, we derived the expected error covariance matrix for a special case of a constant variance noise. Clearly, the discrete-time identification method cannot manipulate the small sampling time case as shown in (41). Also,
5748
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001
we can recognize from (37) and (38) that the relative merit of the continuous-time approach is still valid even if the variance of the noise depends on the sampling time. 3.3. Identification of the Stochastic 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 parts of A ˆ and B ˆ.
[
K ˆ
]
0.5 N
∑{eˆ (ti)}2 N i)1
min V(K ˆ) )
(42)
( )
∂yˆ (ti-1) ∂xˆ (t) d ∂xˆ (t) )A ˆ + Gleˆ (ti-1) - K ˆ , dt ∂kˆ l ∂kˆ l ∂kˆ l ti-1 e t < ti, l ) 1, 2, ..., n else
( )
∂xˆ (t) d ∂xˆ (t) )A ˆ , ti-1 e t < ti, l ) 1, 2, ..., n (51) dt ∂kˆ l ∂kˆ l where
subject to
dxˆ (t) )A ˆ xˆ (t) + B ˆ udelay(t) + K ˆ eˆ (ti-1), ti-1 e t < ti dt (43) yˆ (t) ) Cxˆ (t)
(44)
has 1 for the (n - l + 1)th component and 0 for others. The error y(t) - yˆ (t) is approximately uncorrelated with ∂2yˆ (t)/∂θˆ s2 when the solution is close to true. Then, we approximate the second derivative as follows.
If y(ti) is measured
∂2V(θˆ s)
eˆ (ti) ) y(ti) - yˆ (ti), else eˆ (ti) ) 0
(45)
In the identification of the deterministic part, we can choose any arbitrary ti. However, we restrict ti as a multiple of the basic time length ∆t (that is, ti ) i∆t, i ) 1, 2, ...) in the stochastic part identification for mathematical simplicity. Using (45), we can incorporate unmeasured process outputs. We chose the most probable value of the unmeasured data as zero because the noise e(ti) is assumed to be a white noise. To solve this optimization problem, we use the Levenberg-Marquardt method like Ljung,1 which repeats (46) until the parameters converge within tolerance.
θˆ s(j) ) θˆ s(j-1) -
[
∂2V(θˆ s(j-1)) 2
∂θˆ s
][ -1
+ RI
]
∂V(θˆ s(j-1)) ∂θˆ s
θˆ s ) [kˆ 1 kˆ 2 ... kˆ n]T
(46) (47)
where j denotes the iteration number and R is a small positive value and can be updated every iteration to compromise between the robustness and the convergence rate. The first and second derivatives of the objective function with respect to the adjustable parameters can be calculated as follows: From (42), (48) is derived.
∂V(θˆ s) ∂θˆ s
)-
1
N
∑eˆ (ti)
Ni)1
[
∂yˆ (ti)
(48)
∂θˆ s
]
∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ) ... ∂θˆ s ∂kˆ 1 ∂kˆ 2 ∂kˆ n
T
If y(ti-1) is measured
≈
N
∑
[ ][ ] ∂yˆ (ti) ∂yˆ (ti)
Ni)1 ∂θˆ s
∂θˆ s
T
(52)
In summary, we calculate xˆ (ti) and ∂xˆ (t)/∂θˆ |t)ti for given θˆ s(j-1) by selecting them at every ti while continuously solving the ordinary differential equations of (43) and (51) simultaneously. Next, it is straightforward to calculate the updated parameters θˆ s(j) from (46). Repeat this procedure until the parameters converge. Remarks. 1. Because we choose the two-step approach that identifies the deterministic part and the stochastic part separately, it is difficult to analyze the parameter convergence of the stochastic part because of the effects of the modeling error of the deterministic part. However, if we assume the deterministic part is exact, the same expected error covariance matrix with the deterministic case of (35) and (37) can be derived through the same procedure of the deterministic case. 2. The initial settings are chosen zeros in this research. This is equivalent to choosing the output error model as the initial model. We guess that there are nearly no big problems such as a local minimum and numerical instability due to the initial settings because it identifies only K vector. We confirmed the argument through extensive simulations (refer to Table 2 for an example). 4. Simulation Results The following continuous-time process with time delay 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).
(49) d3z(t) dt
From (43)-(45), (50) and (51) are obtained.
∂yˆ (t) ∂xˆ (t) )C , l ) 1, 2, ..., n ∂kˆ l ∂kˆ l
∂θˆ s2
1
(50)
3
+3
d2z(t) dt
2
+3
dz(t) + z(t) ) u(t-0.4) + e(t) dt (53)
y(t) ) z(t) + e(t)
(54)
The process is activated from t ) 0 and 30 by a random binary sequence (RBS) with a minimum switching time
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5749 Table 1. Typical Convergence Pattern for the Deterministic Identification estimate for iteration residual no. a1 ) 30 a2 ) 3.0 a3 ) 1.0 b1 )1.0 d ) 0.4 variance 0 1 2 3 4 5 6 7 8 9 10
6.0000 7.8227 7.0728 5.6455 4.1115 2.6553 2.9177 2.9737 2.9752 2.9752 2.9752
7.0000 4.8396 5.8892 6.1189 4.0829 2.6571 2.9116 2.9757 2.9776 2.9776 2.9776
2.0000 2.5193 2.2089 2.1749 1.4408 0.8558 0.9617 0.9872 0.9880 0.9880 0.9880
1.0000 2.1421 2.0698 2.1609 1.4158 0.8652 0.9648 0.9896 0.9904 0.9904 0.9904
0.2000 0.8179 0.4255 0.5690 0.5094 0.4012 0.3898 0.3961 0.3964 0.3967 0.3971
1.3056 0.4642 0.0398 0.0124 0.0129 0.0136 0.0101 0.0101 0.0101 0.0101 0.0101
of 1.0 and a magnitude of (5. The variance of e(ti) is 0.01. The proposed and discrete-time prediction error method (oe function) in MATLAB identified the deterministic part of the third order plus time delay model for different sampling times of 0.001 and 0.2. When we apply the discrete-time method, we assume that the time delay is exactly known; meanwhile, the proposed method estimates the time delay additionally. Table 1 shows a typical convergence pattern of the proposed deterministic identification method for the 0.001 sampling time case. R is updated according to Marquardt’s compromise28 with the initial R value of 0.1. We chose poor initial settings on purpose even though we can chose good initial settings in a systematic way as mentioned in section 3.1. Considering the confirmed robustness of the Levenberg-Marquardt method and extensive simulation results for various initial settings, we conclude that the nonlinear optimization method gives us acceptable robustness. Figures 1 and 2 show the worst result for the deterministic part identification among 200 Monte Carlo simulations for 0.001 and 0.2 sampling times, respec-
tively. The discrete-time approach with the small sampling time of 0.001 shows poor identification results due to the small sampling time problem as shown in Figure 1. The discrete-time approach for the large sampling time of 0.2 shows better results than the case of the small sampling time. We also reached the same conclusion with the armax function in MATLAB that identifies a discrete-time ARMAX input model. Contrarily, the proposed method gives almost the same plot as the real process for both cases of small and large sampling times. Table 2 and Figure 3 show the identified result for the stochastic part of the proposed method with the given deterministic part of Figure 2. In this case, we collected historical data sets from t ) 0 and 2000 of the closed-loop system controlled by a simple proportional controller of u(t) ) -y(t). The estimated model in the Laplace domain is as follows:
y(s) )
(
exp(-0.416s)
)
u(s) + s + 3.358s2 + 3.316s + 1.117 -0.022s2 - 0.067s + 0.999 1+ 3 e(s) (55) s + 3.358s2 + 3.316s + 1.117 3
(
)
Table 2 shows a typical convergence pattern of the proposed stochastic identification method. R is updated according to Marquardt’s compromise28 with the initial R value of 0.1. We recognize that the zero initial setting for the nonlinear optimization gives us acceptable robustness. From Figure 3, we can realize that our twostep approach gives acceptable accuracy. Those simulations of Figures 1-3 exemplify that the proposed strategy shows good results for the small and large sampling times. Also, it can incorporate the input time delay efficiently. Previous continuous-time approaches using integral or δ transforms cannot be applied to the large sampling time because the numerical integration
Figure 1. Deterministic identification results for 0.001 sampling time: proposed, dotted line; discrete-time PEM, dashed line; process, solid line.
5750
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001
Figure 2. Deterministic identification results for 0.2 sampling time: proposed, dotted line; discrete-time PEM, dashed line; continuous process, solid line; exact discretized process, dash-dotted line.
Figure 3. Stochastic identification results for 0.2 sampling time: proposed, dotted line; process, solid line.
or derivative of the process output is impossible for the large sampling time. All previous prediction error identification methods do not consider the time delay in a systematic way. 5. Conclusions We proposed a continuous-time prediction error identification method for continuous-time processes with
time delay. It has advantages compared with previous continuous-time and discrete-time identification methods because it can incorporate small, large, and irregular sampling times as well as the input time delay. We justified why discrete-time approaches show an inferior estimation performance to continuous-time approaches by comparing the expected error covariance matrices.
Ind. Eng. Chem. Res., Vol. 40, No. 24, 2001 5751 Table 2. Typical Convergence Pattern for the Stochastic Identification estimate for
iteration no.
k1 ) 0.0
k2 ) 0.0
k3 ) 1.0
residual variance
0 1 2 3 4 5 6 7 8 9 10
0.0000 -0.1713 -0.0253 -0.0231 -0.0221 -0.0221 -0.0221 -0.0221 -0.0221 -0.0221 -0.0221
0.0000 0.4322 -0.0611 -0.0617 -0.0669 -0.0667 -0.0668 -0.0667 -0.0667 -0.0667 -0.0667
0.0000 0.7160 1.0426 0.9936 1.0002 0.9993 0.9994 0.9994 0.9994 0.9994 0.9994
0.0105 0.0102 0.0101 0.0101 0.0101 0.0101 0.0101 0.0101 0.0101 0.0101 0.0101
Acknowledgment This work is supported by the BK21 Project and Center for Ultramicrochemical Process Systems sponsored by KOSEF. Literature Cited (1) Ljung, L. System identification: Theory for the user; PrenticeHall: Englewood Cliffs, NJ, 1987. (2) So¨derstro¨m, T.; Stoica, P. System identification; PrenticeHall: Englewood Cliffs, NJ, 1989. (3) Forssell, U.; Ljung, L. Closed-loop identification revisited. Automatica 1999, 35, 1215. (4) Larimore, W. E. Canonical variate analysis in identification, filtering, and adaptive control. Proceedings of the 29th IEEE Conference on Decision and Control, Honolulu, HI, 1990; p 596. (5) Verhaegen, M. Identification of the deterministic part of MIMO state space model given in innovations form from inputoutput data. Automatica 1994, 30, 1877. (6) Van Overschee, P.; De Moor, B. N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica 1994, 30, 75. (7) Sung, S. W.; Lee, S. Y.; Kwak, H. J.; Lee, I. Continuoustime Subspace System Identification Method. Ind. Eng. Chem. Res. 2001, 40, 2886. (8) Middleton, R. H.; Goodwin, G. C. Digital Control and Estimation; Prentice-Hall: Englewood Cliffs, NJ, 1990. (9) Rao, G. P.; Palanisamy, K. R. Improved algorithm for parameter identification in continuous systems via Walsh functions. Proc. Inst. Electr. Eng., Part D 1983, 130, 9. (10) Cheng, B.; Hsu, N. Analysis and parameter estimation of bilinear systems via block-pulse functions. Int. J. Control 1982, 36, 53. (11) Hwang, C.; Guo, T. Transfer-function matrix identification
in MIMO systems via shifted Legendre polynomials. Int. J. Control 1984, 39, 807. (12) Chou, C. T.; Verhaegen, M.; Johansson, R. Continuoustime Identification of SISO Systems using Laguerre Functions. IEEE Trans. Signal Process. 1999, 47, 349. (13) Mathew, A. V.; Fairman, F. W. Transfer function matrix identification, IEEE Trans. Circuits Syst. 1974, CAS-21, 584. (14) Whitfield, A. H.; Messali, N. Integral-equation approach to system identification. Int. J. Control 1987, 45, 1431. (15) Garnier, H.; Gilson, M.; Zheng, W. X. A bias-eliminated least-squares method for continuous-time model identification of closed-loop systems. Int. J. Control 2000, 73, 38. (16) Sagara, S.; Zhao, Z. Recursive identification of transfer function matrix in continuous systems via linear integral filter. Int. J. Control 1989, 50, 457. (17) Eitelberg, E. Continuous-time system representation with exact macro-difference expressions. Int. J. Control 1988, 47, 1207. (18) Sagara, S.; Yang, Z.; Wada, K. Recursive identification algorithms for continuous systems using an adaptive procedure. Int. J. Control 1991, 53, 391. (19) Johansson, R. Identification of continuous-time models. IEEE Trans. Signal Process. 1994, 42, 887. (20) Johansson, R.; Verhaegen, M.; Chou, C. T. Stochastic theory of continuous-time state-space identification. IEEE Trans. Signal Process. 1999, 47, 41. (21) Kuznetsov, A. G.; Bowyer, R. O.; Clark, D. W. Estimation of multiple order models in the delta-domain. Int. J. Control 1999, 71, 629. (22) Sung, S. W.; Lee, I.; Lee, J. New Process Identification Method for Automatic Design of PID Controllers. Automatica 1998, 34, 513. (23) So¨derstro¨m, T.; Fan, H.; Carlsson, B.; Stefano, B. Least squares parameter estimation of continuous-time ARX models from discrete-time data. IEEE Trans. Autom. Control 1997, 42, 659-672. (24) So¨derstro¨m, T.; Mossberg, M. Performance evaluation of methods for identifying continuous-time autoregressive processes. Automatica 2000, 36, 53. (25) Fan, H.; So¨derstro¨m, T.; Mossberg, M.; Carlsson, B.; Zou, Y. Estimation of continuous-time AR process parameters from discrete-time data. IEEE Trans. Signal Process. 1999, 47, 1232. (26) Pintelon, R.; Schoukens, J.; Rolain, Y. Box-Jenkins contiuous-time modeling. Automatica 2000, 36, 983. (27) Sung, S. W.; Lee, I. Limitations and Countermeasures of PID Controllers. Ind. Eng. Chem. Res. 1996, 35, 2596. (28) Reklaitis, G. V.; Ravindran, A.; Ragsdell, K. M. Engineering Optimization; John Wiley & Sons: New York, 1983.
Received for review January 22, 2001 Revised manuscript received June 4, 2001 Accepted September 5, 2001 IE0100636
