Ind. Eng. Chem. Res. 2002, 41, 4295-4302
4295
System Identification Method for Hammerstein Processes Su Whan Sung* Department of Chemical & Biomolecular Engineering and Center for Ultramicrochemical Process Systems, Korea Advanced Institute of Science and Technology, 373-1 Guseong-dong, Yuseong-gu, Daejeon 305-701, Korea
A new strategy is proposed to identify Hammerstein-type nonlinear processes composed of a nonlinear static function and a linear dynamic subsystem. It completely separates the identification problem of the linear dynamic subsystem from that of the nonlinear static function using a special test signal. Then, quite useful advantages over previous approaches are guaranteed: we can use existing well-established linear system identification methods and their asymptotic properties of parameter estimates as well as excellent techniques for minimal parametrization to identify the linear dynamic subsystem of the Hammerstein process. Also, we can identify the nonlinear static function of the Hammerstein process analytically without any iterative optimization. 1. Introduction It is a well-known fact that linear system identification approaches such as prediction error, instrumental variablemethods,1-3 andsubspaceidentificationmethods4-7 have limitations in describing nonlinear processes. Many researchers have exerted much of their efforts in developing nonlinear system identification methods to overcome the limitations of the linear approaches. These nonlinear approaches can be categorized into two groups according to the model type. One group adopts general types of nonlinear models such as nonlinear autoregressive, moving average, and exogenous input (NAMAX) models, neural networks, Volterra series+, or other various orthogonal series to describe nonlinear dynamics. The other group uses particular types of nonlinear models composed of linear dynamic subsystem and memoryless nonlinear functions such as Wiener, Hammerstein, Hammerstein-Wiener, and so forth.8 If we do not have to consider practical issues, certainly, the system identification methods for the general types are better than those for the particular types in approximating complex nonlinear dynamics. However, the particular types of nonlinear models are definitely preferred if the process is close to one of the particular types. Also, they have advantages over the general approaches in practical issues such as computation time, initial parameter guess, test signal generation, overparametrization problems, and so forth. It has been proven that the particular types can describe the nonlinear dynamics of many chemical, biological, and electrical processes with much less computation cost and simpler plant tests than those for the general types of nonlinear models.9-13 In this study, we will focus on the modeling of the Hammerstein process of Figure 1 among the particular types of the nonlinear processes, where the nonlinear static function precedes the linear dynamic subsystem. Numerous nonlinear system identification methods have been proposed to identify the Hammerstein process. They parametrized the static nonlinear function * To whom all correspondence should be addressed. E-mail:
[email protected]. Telephone: +82-42-869-3992. Fax: +82-42-869-3910.
Figure 1. Hammerstein-type nonlinear process.
in various manners such as using polynomials,9,14,15 a multilayered feedforward neural network,16 or piecewise linear mapping.17 All of them had no choice to use an iterative optimization procedure because they could not separate the identification problem of the linear part from that of the nonlinear part of the Hammerstein model, requiring expensive computation cost and suffering from parameter convergence and initial setting problems compared to analytical approaches. Several approaches have been proposed to identify complex nonlinear static functions without iterative optimization. For examples, Pottmann et al.18 used Kolmogorov-Gabor polynomials to describe highly nonlinear dynamics. An optimal two-stage identification algorithm was proposed, which extracts the model parameters using singular value decomposition (SVD) after estimating an oversized adjustable parameter vector.19 However, their approaches are unnecessarily complicated and/or many redundant adjustable parameters should be estimated because they cannot separate the identification problem of the linear dynamic subsystem from that of the nonlinear static function. Some investigators10,13 used an autotune approach, which utilizes several relay feedback tests with both different relay heights and different dynamic elements inserted in the control loop to determine the Hammerstein model. Even though the autotune approach can identify the nonlinear model efficiently with simple algebraic calculations and simple plant tests, only several points of the nonlinear static function can be estimated because the available information from the relay test is restricted. In this research, I propose a new system identification method for the Hammerstein process. It completely separates the two identification problems of the nonlinear static function and the linear dynamic subsystem. Thus, we can directly utilize existing well-established linear system identification methods1-7 and the riches of theoretical analysis for asymptotic properties of
10.1021/ie0109206 CCC: $22.00 © 2002 American Chemical Society Published on Web 07/25/2002
4296
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002
estimates as well as excellent techniques for minimal realization to identify the linear dynamic subsystem of the Hammerstein process. Furthermore, we can analytically identify the nonlinear static function of the Hammerstein process without any iterative optimization.
In this section, the nonlinear continuous-time Hammerstein process will be explained. And, a special test signal will be designed to completely separate the two identification problems of the linear dynamic subsystem and the nonlinear static function. Finally, we will derive the proposed identification method. 2.1. Hammerstein-type Nonlinear Process. The following Hammerstein process is considered.
v(t) ) p1u(t) + p2u2(t) + ‚‚‚ + pmum(t) for u(t) e umax v(t) ) pm+1 for u(t) > umax
(1)
dx(t) ) Ax(t) + Bv(t - Td) + Ke(t) dt
(2)
y(t) ) Cx(t) + e(t)
(3)
where u(t) and y(t) denote the process input and the process output, respectively. x(t) is the n-dimensional state. e(t) is a piecewise white noise with a holding time of ∆t and E[e2(t)] ) σ2. Td denotes the input time delay. The system matrixes A, B, and K have the following respective forms:
[ ] 0 0 1 l 0 0
0 0 0 l 0 0
... ... ... l ... ...
0 0 0 l 0 1
-an -an-1 -an-2 l -a2 -a1
b1v(n-1)(t - Td) + b2v(n-2)(t - Td) + ... + bn-1v(1)(t - Td) + bnv(t - Td) + k1e(n-1)(t - ∆t) + k2e(n-2)(t - ∆t) + ... + kn-1e(1)(t) + kne(t) (8) y(t) ) z(t) + e(t)
2. Proposed Identification Method
0 1 0 A) l 0 0
z(n)(t) + a1z(n-1)(t) + ... + an-1z(1)(t) + anz(t) )
(4)
B ) [bn bn-1 bn-2 ... b2 b1 ]T
(5)
K ) [kn kn-1 kn-2 ... k2 k1 ]T
(6)
C ) [0 0 0 ... 0 1 ]
(7)
Here, eq 1 is the nonlinear static function and the system of eqs 2 and 3 is the linear dynamic subsystem of the Hamerstein process. Note that we parametrize the nonlinear static function using two conditions of u(t) e umax and u(t) > umax in eq 1. The second condition of u(t) > umax is to represent “saturation”. Even though only “saturation” is exemplified for simplicity, if necessary, we can add more conditions, as many as we want to describe a complex nonlinear function. Refer to two papers2,20 for detailed information on the parametrization for the linear stochastic process. The problem in this research is to estimate the system matrixes A, B, K, Td, and pi, i ) 1, 2, ..., m + 1, from the discretesampled-output data, the given process input, and umax. Remarks. (1) For readers who are familiar with time series expressions, eqs 2 and 3 can be rewritten equivalently like the following continuous-time autoregressive, moving average, exogenous input (ARMAX) model.
(9)
where y(i)(t) denotes the i-th derivative of the continuoustime signal y(t). (2) We choose the polynomial form of eq 1 as the nonlinear static function of the Hammerstein process because it allows us to identify the nonlinear static function without any iterative procedures (which will be discussed later) and it is popular for the nonlinear fitting. Even though we exemplify the polynomials in this research, any other types of parametrization (for example, any orthogonal series) are possible if v(t) is linear with respect to the model parameters of pi, i ) 1, 2, ..., m + 1. (3) We can describe input multiplicity without any problems because the polynomials of eq 1 can have multiple zeros up to m. Sometimes, it is not enough to fit a complex nonlinear static function using only one polynomial. This limitation can be overcome simply by splitting the fitting region with enough conditions. (4) We can use various discrete-time parametrizations such as discrete-time ARMAX time-series or discretetime state space processes instead of the continuoustime state space parametrization of eqs 2 and 3. Choosing a different one among various options of continuous-time, discrete-time, time-series, and state space does not affect the context of the proposed identification method. 2.2. Process Activation. In this research, I suggest a special test signal to activate the Hammerstein process as well as dramatically simplify the identification procedure. The signal satisfies the following conditions: (1) It must contain a binary (two-step) signal of which one step value must be zero and the other step value must be nonzero. (2) It must contain a multistep signal. (3) It must persistently excite the Hammerstein process. As an example, consider the test signal of Figure 2 composed of a two-step (binary) signal from t ) 0 to t ) 50 and a multistep signal from t ) 50 to t ) 100. It satisfies the first and second conditions. We use a random-binary-sequence (RBS) signal of 0 and 0.5 for the two-step signal and a uniformly distributed random noise between - 0.5 and 0.5 for the multistep signal to satisfy the third condition. The above-mentioned condition is a key to separate the identification problem of the linear dynamic subsystem from that of the nonlinear static function completely. To understand this, consider the following argument: we have one degree of freedom to multiply pi, i ) 1, 2,..., m + 1, by any nonzero value (denoted by R) because we can multiply B in eq 2 by the reciprocal (1/R) of the nonzero value not to change the map from the input (u(t)) to the output (y(t)). Thus, if the process input (u(t)) is a two-step signal satisfying the first condition, we can assume u(t) ) v(t) because we have one degree of freedom to change B to compensate the assumption as shown in Figure 3. Now, it is clear that the two-step signal activates the linear dynamic subsystem without activating the nonlinearity of the nonlinear static function. On the other hand, the multistep
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002 4297
mental variable methods,1-3 and subspace system identification methods.4-7 Remarks. (1) All asymptotic properties of parameter estimates and all techniques to determine the model order for linear systems are valid for the proposed identification method because the nonlinear Hammerstein process actually is linear for the two-step signal. This allows us to enjoy elegant and powerful tools for the linear system identification methods. For example, if we use the continuous-time prediction error identification method2 to identify the deterministic part in eq 2 through minimizing the mean square error norm of the prediction error, the estimates are unbiased and the expected error covariance matrix is as follows:
E[(θˆ - θ)(θˆ - θ)T] )
Figure 2. Process activation using the proposed test signal: (a) activated process output; (b) input signal.
signal is to activate the nonlinearity of the nonlinear static function of the Hammerstein process. 2.3. Identification of the Linear Dynamic Subsystem. The core idea to separate completely the identification problem of the linear dynamic subsystem from that of the nonlinear static function is that we can assume u(t) ) v(t) without loss of generality if we use a two-step signal whose one-step value is zero with the other step value being any nonzero value, as in Figure 2. Then, with the assumption of v(t) ) u(t), we can identify straightforwardly the linear dynamic subsystem from u(t) of the two-step signal input and y(t) of the measured process output using existing linear system identification methods such as prediction error, instru-
( ∑( )( ) )
1 1
N
N Ni)1
∂yˆ (ti,θ) ∂yˆ (ti,θ) ∂θ
∂θ
T -1
σo2
where θˆ ) [aˆ 1aˆ 2...aˆ nbˆ 1bˆ 2...bˆ nT ˆ d]T and θ are the estimate vector and actual system parameters, respectively. σo is the variance of the output noise. ∂yˆ (ti,θ)/∂θ denotes the partial derivative of the model output with respect to the actual system parameter vector of θ, and N is the total data number. And, we can use various techniques to determine the model order such as by considering physical insight, examining the spectral analysis estimate, testing ranks in covariance matrixes, examining the information matrix, minimizing Akaike’s information theoretic criterion (AIC), and checking whiteness of residuals. For details, refer to ref 1. (2) Previous approaches do not allow us to utilize useful asymptotic properties developed for linear systems and should estimate an unnecessarily oversized parameter vector because they cannot separate the identification problem of the linear dynamic subsystem from that of the nonlinear static function. (3) Before we estimate the model parameters, we can check if the nonlinear process is Hammerstein-type. Assume that we obtain the process output (named y1(t)) activated by the first binary signal and, next, the second process output (named y2(t)) activated by the second binary signal of a different magnitude. For the Hammerstein process, the first process output is proportional to the second process output (that is, y1(t) ) Ry2(t), where R is a constant).8 Therefore, we can determine if the nonlinear process is Hammerstein-type by checking how much the magnitude of y1(t) - Ry2(t) is close to zero after calculating the R value minimizing the magnitude. There is also another available paper21 to identify the model structure, which uses relay feedback tests. 2.4. Identification of the Nonlinear Function. Once we estimate the linear dynamic subsystem, we can
Figure 3. Actual nonlinear static function (a) and an assumed nonlinear static function without loss of generality (b).
4298
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002
identify the nonlinear static function by minimizing the following cost function from the whole input-output data sets in Figure 2.
min pˆ i,i)1,2,...,m+1
[
V(pˆ 1,pˆ 2,...,pˆ m+1) )
0.5 N
∑(y(ti) - yˆ (ti))2
N i)1
]
(10)
Then, we can solve the differential equations 17-19 and calculate the first derivative of the cost function of eqs 15 and 16. The second derivative of eq 10 is
∂2V(P ˆ) ∂P ˆ2
)-
1
∂2yˆ (ti)
N
∑
Ni)1
(y(ti) - yˆ (ti))
subject to 2
∑
u(t) e umax vˆ (t) ) pˆ m+1 for u(t) > umax
(11)
dxˆ (t) ˆ (y(t) - yˆ (t)) (12) )A ˆ xˆ (t) + B ˆ vˆ (t - T ˆ d) + K dt yˆ (t) ) Cxˆ (t)
(13)
t0 ) 0 < t1 < ... < tN-1 < tN
(14)
where y(t) and yˆ (t) denote the measured process output and the predicted model output, respectively. A ˆ, B ˆ, K ˆ, and T ˆ d represent the estimates of the linear dynamic subsystem obtained in the previous section. The matrix A ˆ - K ˆ C is assumed stable. Equations 12 and 13 represent the optimal one-step-ahead predictor for the linear deterministic-stochastic process of eqs 2 and 3. The optimization problem of eq 10 can be solved analytically because v(t) is linear with respect to the adjustable parameters, pˆ i, i ) 1, 2, ..., m + 1, and the parameters do not affect the initial value of xˆ (t) and u(t). To prove this and derive the analytical solution, consider the following arguments: The derivatives of the cost function with respect to the adjustable parameters can be calculated as follows: From eqs 10 and 15-19 are derived
∂P ˆ
)-
[
1
∂yˆ (ti)
N
∑
Ni)1
(y(ti) - yˆ (ti))
∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ∂yˆ (t) ... ) ∂P ˆ ∂pˆ 1 ∂pˆ 2 ∂pˆ m+1
∂P ˆ
]
(15)
T
(16)
where P ˆ ) [pˆ 1, pˆ 2, ..., pˆ m+1]T. From eqs 11-13 and 17-19 are obtained
( )
[ ][ ]
Ni)1 ∂P ˆ
m
vˆ (t) ) pˆ 1u(t) + pˆ 2u (t) + ... + pˆ mu (t) for
∂V(P ˆ)
+ ∂P ˆ2 1 N ∂yˆ (ti) ∂yˆ (ti)
∂yˆ (t) ∂xˆ (t) )C , i ) 1, 2, ..., m + 1 ∂pˆ i ∂pˆ i
(17)
∂vˆ (t - T ˆ d) ∂xˆ (t) ∂yˆ (t) d ∂xˆ (t) )A ˆ +B ˆ -K ˆ , dt ∂pˆ i ∂pˆ i ∂pˆ i ∂pˆ i i ) 1, 2, ..., m + 1 (18) ∂vˆ (t) ∂vˆ (t) ) ui(t), i ) 1, 2, ..., m and ) 0 for ∂pˆ i ∂pˆ m+1 u(t) e umax ∂vˆ (t) ∂vˆ (t) ) 0, i ) 1, 2, ..., m and ) 1 for ∂pˆ i ∂pˆ m+1 u(t) > umax (19) In eq 18, the initial value of ∂xˆ (t)/∂pˆ i|t)0 is zero because pˆ i cannot change the initial value of the state vector.
∂P ˆ
T
(20)
Here, we realize from eq 19 that ∂2v(t)/∂pˆ i ∂pˆ j, i, j ) 1, 2, ..., m + 1, are zero because the process input is independent of the adjustable parameters. Also, the second derivatives of the state vector are zero at t ) 0 (i.e. ∂2xˆ (t)/∂pˆ i ∂pˆ j|t)0 ) 0, i, j ) 1, 2, ..., m + 1) because the initial value of xˆ (t) is always constant. Then, we reach an important conclusion that ∂2xˆ (t)/∂pˆ i ∂pˆ j, i, j ) 1, 2, ..., m + 1, are always zero or, equivalently, ∂2yˆ (t)/ ∂P ˆ 2 are zero. So, we obtain the following equations from eq 20.
ˆ) ∂2V(P ∂P ˆ2
)
1
N
∑
[ ][ ] ∂yˆ (ti) ∂yˆ (ti)
Ni)1 ∂P ˆ
ˆ) ∂jV(P ∂P ˆi
T
∂P ˆ
) 0, j > 2
(21)
(22)
From eq 22, it is clear that the cost function is quadratic. Then, we can obtain the following linear equation.
∂V(P ˆ ) ∂V(P ˆ) ˆ) ∂2V(P ) |Pˆ )0 + | ˆ )0P ˆ 2 P ∂P ˆ ∂P ˆ ∂P ˆ
(23)
Note, the linear equation is not an approximation but an exact one because the cost function is quadratic. Now, we finally obtain the following optimal solution by setting ∂V(P ˆ )/∂P ˆ ) 0 in eq 23.
[ | ][ | ]
P ˆ* ) -
-1
∂2V(P ˆ) ∂P ˆ
2
P ˆ )0
∂V(P ˆ) ∂P ˆ
P ˆ )0
(24)
In summary, we can easily calculate the first and the second derivative in eq 24 from eqs 15-19 and 21 by solving the differential equation 18. Then we obtain the optimal solution of eq 24 for the nonlinear static function of the Hammerstein model. Remarks. (1) We obtain the analytic solution of eq 24 for the nonlinear static function without any iterative procedure. So, there are no problems related to convergence, computation time, and initial settings. Previous iterative approaches to identify the nonlinear static function of the Hammerstein process cannot be free from these problems. (2) Because we identify the nonlinear static function after we estimate the linear dynamic subsystem, it is difficult to analyze the parameter convergence of the nonlinear part due to the effects of the modeling error of the linear dynamic subsystem. As a restricted case, if we can assume the identified linear dynamic subsystem is exact, the asymptotic properties of the proposed nonlinear identification method are similar to those of the discrete-time or continuous-time prediction error identification methods,1,2 like the following:
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002 4299
E[(P ˆ - P)(P ˆ - P)T] )
( ∑( )( ) )
1 1
N
N Ni)1
∂yˆ (ti,P) ∂yˆ (ti,P) ∂P
T -1
σ2
∂P
where P ˆ ) [pˆ 1,pˆ 2,...,pˆ m+1]T is the estimate vector and P denotes actual system parameters. ∂yˆ (ti,P)/∂P denotes the partial derivative of the model output with respect to the actual system parameter vector of P and N is the total data number. For details, refer to the corresponding papers.1,2 3. Case Study Two processes are simulated to test the proposed identification method: one is a simple Hammerstein process, and the other is a more complex Hammerstein process that has input saturation and input multiplicity. Also, the performance of the proposed identification method is experimentally confirmed for a thermal microsystem, and the use of the identified Hammerstein model for controller design is exemplified. Example 1. Consider the following simple Hammerstein-type nonlinear process that has a third-order plus time delay as the linear dynamic subsystem and an exponential function as the nonlinear static function:
v(t) ) (1 - exp(-5u(t))
(25)
y(s) exp(-0.3s) ) v(s) s4 + 4s3 + 6s2 + 4s + 1
(26)
where u(t), v(t), and y(t) denote the process input, the output of the nonlinear static function, and the process output, respectively. Figure 2 shows the activated process output by the RBS signal followed by the uniformly distributed random noise. Here, the minimum switching time is 1.0 for both signals. Note that the dynamic behavior of the process output in the RBS signal region is totally different from that of the process output in the uniformly distributed random noise region due to the strong nonlinearity. A third-order plus time delay model is used to approximate the linear dynamic subsystem of the Hammerstein process. Note that we can choose any favorite one among numerous well-established linear system identification methods to estimate the linear dynamic subsystem model. In this study, I obtained the following model from the input-output data sets from t ) 0 to t ) 50 using the continuous-time prediction error identification method2 that provides the optimal parameters in minimizing the prediction error.
1.837 exp(-0.572s) yˆ (s) (27) ) 3 vˆ (s) 2.590s + 4.834s2 + 3.733s + 1.000 Next, I estimated the following nonlinear static function using the analytical solution of eq 24 from the inputoutput data sets from t ) 0 to t ) 100 as follows.
vˆ (t) ) 2.733u(t) - 6.684u2(t) + 10.811u3(t) 15.855u4(t) + 18.026u5(t) - 7.644u6(t) (28) To check the accuracy of the identified model of eqs 27 and 28, rewrite the actual process of eqs 25 and 26 as follows:
y(s) 1.837 exp(-0.3s) ) 4 v(s) s + 4s3 + 6s2 + 4s + 1
(29)
Figure 4. Identification result for the linear dynamic subsystem in example 1.
Figure 5. Identification result for the nonlinear static function in example 1.
v(t) )
(1 - exp(-5u(t))) 1.837
(30)
As shown in eqs 27 and 28 and 29 and 30, there are structural plant/model mismatches. Nevertheless, the estimated linear dynamic subsystem is excellent, as shown in Figure 4 because the two-step signal does not activate the nonlinearity of the nonlinear static function and the large time delay of the third-order plus time delay model of eq 27 can compensate effectively the order difference between the process of eq 29 and the model of eq 27. And, the actual nonlinear static function of eq 30 and that of the estimated polynomial of eq 28 are compared in Figure 5. The identification result for the nonlinear static function is almost exact for the whole activated region. Previous approaches using relay feedback can identify only several points of the nonlinear static function. The other previous approaches should solve iterative nonlinear optimization problems or should estimate unnecessarily many adjustable parameters to identify the nonlinear static function. In contrast with them, the proposed nonlinear identification method can identify the whole activated nonlinear region without any iterative procedures. Also, it provides the user an option to adopt his/her favorite one
4300
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002
Figure 7. Identification result for the linear dynamic subsystem in example 2.
Figure 6. Process activation using the proposed test signal in example 2: (a) activated process output; (b) input signal.
among existing linear system identification methods and well-established techniques to parametrize the model minimally for the linear dynamic subsystem identification. Example 2. The second example is the following Hammerstein-type nonlinear process, of which the static nonlinear function shows input multiplicity and saturation. The linear dynamic subsystem is the same as that for example 1, but there is measurement noise. 2
3
v(t) ) 5.3u(t) - 2.8u(t) + 0.43u(t)
for u(t) e 4.5
v(t) ) 6.33375 for u(t) > 4.5 y(s) )
exp(-0.3s) 4
v(s) + e(s)
s + 4s3 + 6s2 + 4s + 1
Figure 8. Identification result for the nonlinear static function in example 2.
yˆ (s) 1.005 exp(-0.577s) (33) ) vˆ (s) 2.714s3 + 4.892s2 + 3.819s + 1.000 Next, I estimated the following nonlinear static function using the analytical solution of eq 24 from the inputoutput data sets from t ) 0 to t ) 100.
vˆ (t) ) 5.257u(t) - 2.783u2(t) + 0.429u3(t) for u(t) g 4.5
(31)
v(t) ) 6.328 for u(t) > 4.5
(32)
The estimated linear dynamic subsystem of eq 33 and the nonlinear static function of eq 34 are almost exact for the whole activated region, as shown in Figures 7 and 8, confirming that the proposed method can effectively incorporate input multiplicity, input saturation, and measurement noises; also, we can approximate a more complex nonlinearity by further splitting the fitting region with more conditions. Example 3. The proposed identification method is applied to the thermal microsystem of Sung et al. (2002).22 It is composed of a silicon-based microreactor and hardware/software for automation. The siliconbased microreactor is integrated with a thin-film plati-
where e(t) is a white noise of 0.05 variance. The process is activated by the RBS signal (composed of 0.0 and 4.030 620 1) followed by the uniformly distributed random noise (between 0.0 and 7.0), as shown in Figure 6. Here, the minimum switching time is 1.0 for both signals. I obtained the following linear model from the inputoutput data sets from t ) 0 to t ) 50 using the continuous-time prediction error identification method.2
(34)
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002 4301
Figure 9. Process activation using the proposed test signal for the thermal microsystem: (a) activated process output; (b) input signal.
Figure 10. Identification results for the thermal microsystem: (a) nonlinear model; (b) linear model.
num sensor and a heater by semiconductor fabrication technologies. The thermal system is very important for DNA polymerase chain reaction (PCR), which requires a rapid temperature control to shorten the total running time and a precise temperature control for high efficiency. The thermal microsystem adjusts the voltage of the platinum heater (control output) to control the temperature of the microreactor (process output). It measures the reactor temperature by reading the resistance of the platinum sensor. For detailed information on the system, refer to the corresponding paper.22 Figure 9 shows the activated process data using the proposed special test signal. I estimated the linear dynamic subsystem of the Hammerstein model using the continuous-time prediction error identification method2 from the input-output data between t ) 0 and t ) 50 and the nonlinear static function using the analytical solution of eq 24 from the input-output data between t ) 0 and t ) 130 as follows.
The Hammerstein model of eqs 35 and 36 well describes the nonlinear dynamic behavior of the thermal microsystem while the linear model identified from the input-output data between t ) 0 and t ) 50 shows poor model performance, as shown in Figure 10. Figure 11 illustrates how to use the identified Hammerstein model for controller design, where, ys is the set point and PI denotes the linear proportional-integral controller. The control strategy of Figure 11a linearizes the nonlinearity of eq 35 so that the dotted-line-box is approximately equivalent to the linear system of eq 36. Then, the linear PI controller controls the linear system of eq 36 as shown in Figure 11b. So, we can achieve the same control performance as expected from the linear control system. Figure 12 compares the control results of the linear PI controller and the nonlinear PI controller of Figure 11a, where the tuning parameters are tuned for the operating region between 45 and 50 °C and the sampling time is 0.055 s. As expected, the linearizing nonlinear controller shows much better performances than the linear controller. The control performances of the nonlinear PI controller are almost the same for all operating regions while the linear PI controller shows acceptable performance for only the operating region between 45 and 50 °C.
vˆ (t) ) ˆf(u(t)) ) 0.8335(u(t) - 1.5) + 0.2354(u(t) - 1.5)2 (35) yˆ (s) 43.7587s + 8.8236 ) vˆ (s) s2 + 3.0252s + 0.2653
(36)
4302
Ind. Eng. Chem. Res., Vol. 41, No. 17, 2002
and techniques for minimal realization. Also, we can identify the nonlinear static function analytically without any iterative nonlinear optimization. Acknowledgment This work was supported by Korea Research Foundation Grant KRF-2001-003-E00291. Literature Cited
Figure 11. Nonlinear PI controller (a) and equivalent linear control system (b).
Figure 12. Control results for the thermal microsystem: (a) nonlinear PI controller; (b) linear PI controller.
4. Conclusions A new strategy to identify Hammerstein-type nonlinear processes has been proposed. It completely separates the two identification problems of the linear dynamic subsystem and the nonlinear static function using a special test signal. The proposed method has remarkable advantages over previous approaches: we can use any existing well-established linear system identification methods for the linear dynamic subsystem identification. This means we can utilize elegant theories for asymptotic properties of parameter estimates
(1) Ljung L. System identification: Theory for the user; PrenticeHall: Englewood Cliffs, NJ, 1987. (2) Sung, S. W.; Lee, I. Prediction Error Identification Method for Continuous-time Processes with Time Delay. Ind. Eng. Chem. Res. 2001, 40, 5743. (3) So¨derstro¨m, T.; Stoica, P. System identification; PrenticeHall: Englewood Cliffs, NJ, 1989. (4) Larimore, W. E. Canonical Variate Analysis in Identification, Filtering, and Adaptive Control. Proceedings of the 29th IEEE Conference On Decision and Control, 1990, Honolulu, HI; 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) Haber, R.; Unbehauen, H. Structure Identification of Nonlinear Dynamic Systems - A Survey on Input/Output Approches. Automatica 1990, 26, 651. (9) Eskinat, E.; Johnson, S. H.; Luyben, W. L. Use of Hammerstein Models in Identification of Nonlinear Systems. AIChE J. 1991, 37, 255. (10) Luyben, W. L.; Eskinat, E. Nonlinear Auto-tune Identification. Int. J. Control 1994, 59, 595. (11) Huang, H.; Lee, M.; Tang, Y. Identification of Wiener Model using Relay Feedback Test. J. Chem. Eng. Jpn. 1998, 31, 604. (12) Maksym, G.; Kearney, R. E.; Bates, H. T. Nonparametric Block-Structured Modeling of Lung Tissue Strip Mechanics. Ann. Biomed. Eng. 1998, 26, 242. (13) Balestrino, A.; Landi, A.; Ould-Zmirli, M.; Sani, L. Automatic Nonlinear Auto-Tuning Method for Hammerstein Modeling of Electrical Drives. IEEE Trans. Ind. Electron. 2001, 48, 645. (14) Zhu, Y. Identification of Hammerstein Models for Control using ASYM. Int. J. Control. 2000, 73, 1692. (15) Westwick, D.; Kearney, R. Separable Least Squares Identification of Nonlinear Hammerstein Models: Application to Stretch Reflex Dynamics. Ann. Biomed. Eng. 2001, 29, 707. (16) Al-duwaish, H.; Karim, M. N. A New Method for the Identification of Hammerstein Model. Automatica 1997, 33, 1871. (17) Vo¨ro¨s, J. Modeling and Parameter Identification of Systems With Multisegment Piecewise-Linear Characteristics. IEEE Trans. Autom. Control 2002, 47, 184. (18) Pottmann, M.; Unbehauen, H.; Seborg, D. E. Application of a General Multi-Model Approach for Identification of Highly Nonlinear ProcessessA Case Study. Int. J. Control 1993, 57, 97. (19) Bai, E. An Optimal Two-Step Identification Algorithm for Hammerstein-Wiener Nonlinear Systems. Automatica 1998, 34, 333. (20) Pintelon, R.; Schoukens, J.; Rolain, Y. Box-Jenkins Contiuous-time Modeling. Automatica 2000, 36, 983. (21) Huang, H. P.; Lee, M. W.; Tsai, C. Y. Structure Identification for Block-oriented Nonlinear Models Using Relay Feedback Tests. J. Chem. Eng. Jpn. 2001, 34, 748. (22) Sung, S. W.; Park, S.; Yoon, D. S.; Lim, G. Modeling and Control of Thermal Microsystems. J. Micromech. Microeng., submitted.
Received for review November 14, 2001 Revised manuscript received June 6, 2002 Accepted June 24, 2002 IE0109206
