Ind. Eng. Chem. Res. 2004, 43, 6125-6135
6125
Robust Identification of Continuous Parametric Models Based on Multiple Sinusoidal Testing under Slow or Periodic Disturbances Shyh-Hong Hwang* and Han-Chern Ling Department of Chemical Engineering, National Cheng Kung University, Tainan, Taiwan 70101, R.O.C.
Shing-Jia Shiu Refining & Manufacturing Research Center, Chinese Petroleum Corporation, Chia-Yi, Taiwan 60036, R.O.C.
This paper addresses the identification of models using multiple sinusoidal forcing subject to practical difficulties such as unknown initial states, offsets, slow and periodic disturbances, noise, and unknown model structures. A linear regression equation is derived by integrating the system equation excited by a single sinusoid and extended to the case of concurrent multiple sinusoids. The regression equation can be used to estimate the model parameters including time delay in a least-squares fashion. A simple scheme based on the condition number of the matrix formed by the regression vector is presented to infer the best model order. Two-stage sinusoidal testing, represented as two sets of sinusoids applied sequentially, is then developed to generate response data that are informative enough to deal with the aforementioned identification difficulties. This requires only the application of the regression equation for concurrent multiple sinusoids modified in a clever and sequential manner. 1. Introduction The most direct way of obtaining a continuous dynamic model of a process is to find the parameters that fit the experimental step response. Smith1 proposed graphical techniques to identify a first- or second-order transfer function with time delay from a step response. These techniques utilize limited data points selected for fitting and encounter difficulties in distinguishing between second- and third- or high-order systems. Wang and Zhang2 proposed a least-squares method to identify a dynamic model using the entire step-response curve. The method is based on a linear regression equation derived for the least-squares estimation of the composite parameters from which the model parameters including time delay can be recovered. To diminish the effect of measurement noise, the method incorporates the instrumental variable approach in its estimation algorithms.2-4 The next level of dynamic testing is direct sinusoidal testing, which is a useful way of obtaining precise frequency-response data of a plant.5,6 The input of a plant is varied sinusoidally at a specified frequency. After all transients have died out and a steady oscillation in the output has become established, one point on the frequency-response curve is identified from the steady oscillation data. The complete frequency-response curve can then be found experimentally by varying the frequency of the sinusoidal input many times over the range of interest. This procedure can be very time-consuming, especially when applied to a slow process. Pulse testing is considered to be more efficient than direct sinusoidal testing in that it can yield the complete frequency-response curve in a single experimental test.6 * To whom correspondence should be addressed. Fax: 8866-2344496. Tel.: 886-6-2086969 ext. 62661. E-mail: shhwang@ mail.ncku.edu.tw.
A pulse of suitable shape and width is introduced in an input variable of the process. The recorded input and output pulses are employed to calculate frequencyresponse data in the range of interest through simplifications of one-sided Fourier transforms.5 A transfer function model can then be obtained by fitting frequencyresponse data using a graphical or an iterative leastsquares technique.6,7 Ham and Kim8 presented a method to acquire a second-order plus delay transfer function based on a rectangular pulse response. However, the method utilizes only the latter portion of the response after the input pulse vanishes and is essentially a stepresponse approach. The feasibility of the above continuous identification methods is doubtful for the requirements that the plant test should be started from a consistent steady state, that no disturbances are present during the test, and that the model structure (order and delay) is given a priori.9 Considering that, in practice, the initial states of the system are possibly neither known nor steady, some least-squares approaches treat these states as additional parameters to be estimated.10-13 However, as the number of parameters to be estimated increases, these approaches become sensitive to noise and require that the input signal be more complicated than step and sinusoidal waves to ensure a unique solution to the parameter estimation problem. Sagara and Zhao14 established a multiple-integral filter to produce filtered signals for model identification without involving unknown initial states. Nevertheless, their method is rather sensitive to noise and model structure mismatch. Many continuous identification methods that utilize time-consuming corrective action have been developed to remedy the effect of an offset arising in the course of a plant test.15-17 However, in addition to offsets, processes are also subject to disturbances that vary periodically or slow disturbances such as trends and drift.6,18,19 As an example, a periodic disturbance in cooling water temperature can often be closely tied to
10.1021/ie030706c CCC: $27.50 © 2004 American Chemical Society Published on Web 08/07/2004
6126 Ind. Eng. Chem. Res., Vol. 43, No. 19, 2004
diurnal fluctuations in ambient conditions. Although identification under diverse disturbances is often treated in discrete time, it is seldom considered in continuous time.19,20 This paper is intended to develop methods that identify a continuous parametric model with time delay using multiple sinusoidal testing in the face of practical difficulties. In other words, the methods should allow the test to be started from any states and reject the effects of offsets and slow or periodic disturbances. If the model order is not presumed, the methods should infer the best order from the measured output. Moreover, these methods should be robust with respect to noise and model structure mismatch.
sinusoidal input is applied at time zero, the following integrals hold for t g d
(-1)iA Iiu(t - d) ) cos[ω(t - d) + φ + (i - 1)π/2] + ωi i (-1)j-1A cos[φ + (j - 1)π/2] (t - d)i-j , (i - j)! j)1 ωj i ) 1, 2, ..., n (6)
∑
n
A linear nth-order plus delay system can be expressed as
y(n)(t) + an-1y(n-1)(t) + ‚‚‚ + a1y(1)(t) + a0y(t) )
an-jIjy(t) + ξ2A sin[ω(t - d) + φ] + ∑ j)1 ξ1 A cos[ω(t - d) + φ] +
b1u(1)(t - d) + b0u(t - d) + l(t) (1) where y(t) and u(t) are, respectively, the output and the input measurements; n is the system order; and d denotes the time delay. The term l(t) is used to account for unknown disturbances arising during the identification experiment. To avoid the use of various time derivatives, u(i)(t) and y(i)(t), we define a multiple integral of the signal x(t) as
∑ j)0
- d)jηn-1-j j!
(7)
∑ j)0
[
]
(-1)jA cos(φ + jπ/2)
i
ηi )
bi-j
ωj+1
,
i ) 0, 1, ..., n - 1 (8a)
and ξ1 and ξ2 can be expressed as
ξ1 )
int[(n+1)/2](-1)jb
∫0t∫0τ ‚‚‚ ∫0τ x(τ1) dτ1 dτ2 ‚‚‚ dτj, j
n-1(t
where ηi is the following function of A, ω, and φ
bn-1u(n-1)(t - d) + bn-2u(n-2)(t - d) + ‚‚‚ +
Ijx(t) )
]
Substituting the above expression into eq 4 gives
y(t) ) -
2. Single Sinusoidal Testing
[
∑ j)1
n+1-2j
ω2j-1
(8b)
2
j ) 0, 1, ..., n (2)
Suppose that the system is initially at a known steady state, that the input signal remains constant for at least a period of d before the identification test is started, and that no disturbances occur during the test. Thus, expressing y(t) and u(t) as deviation variables leads to the following zero conditions
y(i)(0) ) 0,
i ) 0, 1, ..., n - 1
u(t) ) 0,
-d e t < 0
l(t) ) 0
ξ2 )
(3c)
∑ j)1
n-2j
ω2j
(8c)
In the above equations, int(z) rounds the real number z to the nearest integer toward zero. Expanding the terms containing d in eq 7 results in the following linear regression equation
y(t) ) φ(t)Tθ,
(3a) (3b)
int(n/2)(-1)jb
tgd
(9)
where
φ(t) ) [-I1y(t) -I2y(t) ‚‚‚ -Iny(t) 1 t ‚‚‚ tn-1 A sin(ωt + φ) A cos(ωt + φ)]T
Operating on eq 1 with the multiple integrals in eq 2 gives rise to
h1 θ h 2 ‚‚‚ θ h n+2 ]T θ ) [an-1 an-2 ‚‚‚ a0 θ
y(t) ) -an-1I1y(t) - an-2I2y(t) - ‚‚‚ - a0Iny(t) + bn-1I1u(t - d)+ bn-2I2u(t - d) + ‚‚‚ + b0Inu(t - d) (4)
Invoking eq 9 for sampling instants t ) ti (gd), i ) 1, 2, ..., N (.2n + 2) and applying the ordinary least-squares yields the estimation of θ as N
Equation 4 represents the underlying identification model obtained via multiple integrals and under the zero conditions. Consider a single sinusoidal signal given by
u(t) )
{
A sin (ωt + φ), t g 0 0, t
