Closed-Loop Identification of Dynamic Models for Multivariable

1460

Ind. Eng. Chem. Res. 2011, 50, 1460–1472

Closed-Loop Identification of Dynamic Models for Multivariable Systems with Applications to Monitoring and Redesign of Controllers Jyh-Cheng Jeng* and Yan-You Lin Department of Chemical Engineering and Biotechnology, National Taipei UniVersity of Technology, Taipei 106, Taiwan

Closed-loop identification is usually desirable, but more difficult compared to open-loop identification because of the correlation between the noise and the process inputs. This article presents two new closed-loop identification methods for multivariable systems to solve this problem. Based on one closed-loop test, the complete models, including process and disturbance dynamic models, can be identified simultaneously. The finite impulse responses (FIRs) of the closed-loop systems are first estimated by an efficient recursive algorithm derived from a subspace identification method. Then, the fast Fourier transform (FFT) and inverse FFT (IFFT) techniques are used to construct the frequency responses and FIR models, respectively, of process and disturbance dynamics. These algorithms do not need prior knowledge of the process or assumptions about the model structure. Applications of the proposed method to control system monitoring and controller redesign are thus presented. These applications can be accomplished before the process model has been identified and would be very helpful for improving the control performance. The simulation results demonstrate the effectiveness of the proposed identification method and its application. 1. Introduction Closed-loop identification refers to the identification of process models using data measured under feedback control, and it plays an important role in improving existing underperforming controllers, also enabling the identification of systems that cannot operate in open-loop fashion. However, correlation between the noise in the process measurements and process inputs imposes a fundamental limitation on the application of standard open-loop identification methods with closed-loop data.1 Over the years, closed-loop identification methods have been introduced into the literature, and they fall into three main groups: direct, indirect, and joint input-output approaches.2 In the direct approach, measurements of the process inputs and outputs are used the same way as in the open-loop identification approach, ignoring any feedback. The feedback mechanism does not need to be known a priori. In the indirect approach, some closed-loop systems are identified between a reference signal to the process and/or controller output. Then, the process model is retrieved using knowledge of the operating controller. The feedback mechanism thus needs to be known. In the joint input-output approach, the process output and the process input are jointly regarded as the output of a system driven by the reference signal and noise. The process model is determined from an estimate of this augmented system. The closed-loop structure needs to be known, but the controller can be unknown. Closed-loop identification using prediction error methods (PEMs) was summarized in a survey article of Forssell and Ljung,2 who showed the connections between these approaches and displayed similarities and differences in the asymptotic properties of the resulting estimates. However, special parametrizations of the dynamics and noise models are needed using PEMs, and such parametrizations are not straightforward when prior knowledge of the process is not available. Subspace identification methods (SIMs) have enjoyed tremendous development in the past 15 years and offered an attractive alternative to PEMs because of their simple and general parametrization * To whom correspondence should be addressed. Tel.: 886-227712171 ext. 2540. E-mail: [email protected].

for multivariable systems.3 Moreover, in SIMs, the subspace matrices are directly calculated from the input-output data Hankel matrices without any iteration compared to the iterative or nonlinear optimization schemes used in PEMs. The fundamental problem for closed-loop identification is the correlation between the unmeasured noise and the process input.1 Therefore, the traditional open-loop SIMs are biased under closed-loop conditions, which requires special treatment and/or additional assumptions.3,4 The models identified by PEMs and SIMs are typically in a discrete form. However, some continuous transfer function (CTF) models are preferred in industrial process control. Wang et al.5 proposed a direct identification approach using a closedloop step or relay test, where the frequency responses of the process are first identified and then transformed into a parametric CTF model. Silva et al.6 proposed an indirect closed-loop identification technique based on Laguerre series expansions of the closed-loop step response. Then, the step response of the open-loop process is reconstructed, and a low-order CTF model is identified. A similar method was also proposed by Park et al.7 These methods apply only to single-input/singleoutput (SISO) systems. Recently, based on the ideas in Wang et al.,5 Mei et al.8 presented the decentralized closed-loop parameter identification for multivariable processes from step responses. However, their method is applicable only to square multi-input/multi-output (MIMO) processes under decentralized control and requires sequential step tests, so that the experimental time is long for higher-dimension systems. In addition, most methods reported in the literature focus on the identification of process model, but the disturbance model for the measured disturbance is usually not considered. The disturbance model is very useful for the design of advanced control systems, such as feedforward control. To overcome the difficulties mentioned above, two novel closed-loop identification methods for multivariable systems are proposed in this article. The first one is an indirect approach, and the second one can be classified as a joint input-output approach. The complete models, including process and disturbance dynamic models, can be identified using closed-loop data

10.1021/ie1010146  2011 American Chemical Society Published on Web 10/12/2010

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

from one test. The algorithms do not require any prior parametrization about the models. The finite impulse response (FIR) models of several closed-loop systems are first estimated by a recursive algorithm derived from a SIM. The use of a SIM for indirect closed-loop identification can bypass the special treatment required for traditional direct closed-loop identification. Also, the recursive algorithm makes the estimation of FIR models more efficient. Then, the fast Fourier transform (FFT) technique is used to construct the frequency responses of the process and disturbance dynamics. Finally, the FIR models of the process and disturbance dynamics are calculated by the inverse FFT (IFFT) technique. Furthermore, the information obtained during the identification procedures can be applied to control system monitoring and controller redesign to improve the control performance of existing controllers. Therefore, the proposed methods, including identification, monitoring, and controller retuning, can be combined to develop an intelligent control system. The remainder of this article is organized as follows: Section 2 presents the subspace-based recursive identification of FIR model. In section 3, two identification algorithms using closedloop data are proposed. The applications of the proposed procedures to monitoring and controller redesign are provided in section 4. The effectiveness of the proposed methods is demonstrated through simulation examples in section 5. Finally, conclusions are drawn in section 6. 2. Identification of FIR Model The FIR model (hi, i ) 1, 2, 3, ...) has a number of obvious advantages from the viewpoint of system identification.9 For example, the determination of an FIR model requires less a priori knowledge than do determinations of parametric models, and this model can be identified more satisfactorily in the presence of noise. For open-loop stable systems, the FIR will decay to zero after some i > p. The typical method for identifying the FIR model is the least-squares estimation. However, before proceeding with the estimation, the value of the settling-time parameter p must be chosen in advance. To obtain an accurate result, p should be picked according to the condition hi ≈ 0 for i > p, which might result in some computational complexities, although, usually, the solution can be repeated with progressively increasing p values until a satisfactory fit has been achieved. However, a large p value increases the computational difficulties associated with high-order matrix inversion in the least-squares estimation. To overcome this difficulty, an alternative FIR estimation method based on part of the subspace identification algorithm is proposed. 2.1. Alternative Method for FIR Estimation. A brief description of the part of the subspace identification method that is useful for FIR estimation is first summarized as follows. The innovation form of state-space representation of linear systems is given by xk+1 ) Axk + Buk + Kek yk ) Cxk + Duk + ek

(1)

where xk ∈ Rn, uk ∈ Rnu, yk ∈ Rny, and ek ∈ Rny are the system state, input, output, and white noise (innovations), respectively. The indexes n, nu, and ny denote the number of states, inputs, and outputs, respectively. Matrices A, B, C, and D are system matrices with appropriate dimensions, and K is the steady-state Kalman gain. Based on the innovation form in eq 1, an extended state-space model can be formulated as3,4

Yf ) ΓNXf +

HdNUf

+

1461

HsNEf

(2)

where the subscript f denotes the future horizon. In eq 2, the future input, output, and innovation block-Hankel matrices Uf, Yf, and Ef, respectively, are defined as

Uf } UN|2N-1

Yf } YN|2N-1

Ef } EN|2N-1

( ( (

) ) )

uN uN+1 · · · uN+j-1 uN+j uN+1 uN+2 · · · } ∈ RnuN×j · ·. l l l u2N-1 u2N · · · u2N+j-2

(3)

yN yN+1 · · · yN+j-1 yN+j yN+1 yN+2 · · · ∈ RnyN×j } · ·. l l l y2N-1 y2N · · · y2N+j-2

(4)

eN eN+1 · · · eN+j-1 eN+j eN+1 eN+2 · · · ∈ RnyN×j } · ·. l l l e2N-1 e2N · · · e2N+j-2

(5)

where the subscript k1|k2 indicates that the first column of the matrix starts from time k1 and ends by time k2. The state matrix, Xf, is defined as Xf } XN ) (xN xN+1 · · · xN+j-1 ) ∈ Rn×j

(6)

The extended observability matrix ΓN is given by

( )

C CA ΓN ) ∈ RnyN×n l CAN-1

(7)

The lower-triangular block-Toeplitz matrices HNd and HNs are given by

HdN )

( (

D CB l

0 D l

··· 0 ··· 0 · ·. l

(8)

··· 0 ··· 0 nyN×nyN · ·. l ∈ R

(9)

CAN-2B CAN-3B · · · D

HsN )

I CK l

0 I l

) )

∈ RnyN×nuN

CAN-2K CAN-3K · · ·

I

The Kalman state Xf is unknown, but it can be estimated from past input and output data as4 Xf ) Lp

( )

Yp } LpWp Up

(10)

where the subscript p denotes the past horizon; Lp is the subspace matrix corresponding to past data; and the past input and output block-Hankel matrices Up and Yp, respectively, are defined as

Up } U0|N-1 }

(

u0 u1 l uN-1

)

u1 · · · uj-1 uj u2 · · · ∈ RnuN×j · ·. l l uN · · · uN+j-2

(11)

1462


(

Yp } Y0|N-1 }

y0 y1 l yN-1

)

y1 · · · yj-1 yj y2 · · · ∈ RnyN×j l · ·. l yN · · · yN+j-2

Substituting eq 10 into eq 2, one can write

The lower block of eq 16 can be rewritten as

(12)

Yf ) LwWp + HdNUf + HsNEf

(13)

where Lw ) ΓNLp. Thus, the subspace matrices Lw and HNd can be identified from data Hankel matrices using regression techniques such as least-squares method, provided that the following conditions are satisfied:4 (i) The input uk is uncorrelated with ek. (ii) The input uk is persistently exciting. (iii) The number of measurements is sufficiently large (i.e., j is sufficiently large). With the input-output data, the solution of this least-squares problem is given by (Lw HdN ) ) Yf

( ) Wp Uf

†

) Yf(WpT UfT )

[( )

]

Wp (WpT UfT ) Uf

-1

(14)

The matrix Lw is useful for the estimation of system matrices (A, B, C, D, K), but it is not needed here because our goal is to identify the FIR model, which can be extracted from the subspace matrix HNd . The relationship between FIR sequences, h ) {hi}i)1,2,3,..., and state-space matrices is given as hi ) CAi-1B,

i ) 1, 2, 3, ...

(15)

Notice that each hi is a matrix with dimensions ny × nu. From eq 8, it is found that the first column of matrix HNd contains the first N sequences of the FIR. In other words, the first N sequences of the FIR (i.e., hi, i ) 0, 1, ..., N - 1) can be obtained from the solution of HNd in eq 14. However, as the number of the variables increases, the size of the data Hankel matrices can be prohibitively high, especially for systems with a long settling time. The numerical problem becomes a potential concern.10 Here, it is proposed that the value of N can be chosen small enough that the computational difficulty associated with high-order matrix inversion in eq 14 is avoided. Furthermore, the whole sequence of the FIR model can be computed piece-by-piece using a recursive algorithm as described in the next section. 2.2. Recursive Algorithm for FIR Estimation. After the first N sequences of the FIR (or HNd ) have been obtained, the next N sequences of the FIR can be computed in a recursive manner. The extended version of eq 13, by extending the end point of future-horizon subscript f from (2N - 1) to (3N - 1), can be formulated as

(

) ( )

(

)( ) ( )

Lw YN|2N-1 Hd2N,11 Hd2N,12 UN|2N-1 + ) (2) Wp + d Y2N|3N-1 Lw H2N,21 Hd2N,22 U2N|3N-1 Hs2N,1 E (16) Hs2N,2 N|3N-1

d d d ) H2N,22 ) HNd , H2N,12 ) 0, and where H2N,11

Hd2N,21

(

CAN-1B CANB ) l

CAN-2B CAN-1B l

··· ··· ·

·.

CB CAB l

CA2N-2B CA2N-3B · · · CAN-1B

)

∈ RnyN×nuN

(17)

Y2N|3N-1 - HdNU2N|3N-1 ) Lw(2)Wp + Hd2N,21UN|2N-1 + Hs2N,2EN|3N-1 (18) d As a result, the matrix (L(2) w H2N,21) can be solved by right-handside of eq 14, except that Yf is replaced with (Y2N3N-1 HNd U2N3N-1). Then, the next N sequences of FIR (i.e., hi, i ) N, d . N + 1, ..., 2N - 1) are obtained in the first column of H2N,21 Likewise, the recursive algorithm, for m ) 3, 4, ..., is derived as described below. d ) is solved by the equation (1) The matrix (Lw(m) HmN,m1

[( )

d (Lw(m) HmN,m1 ) ) Y˜f(WpT UfT )

Wp (WpT UfT ) Uf

]

-1

(19)

where Y˜f ) YmN|(m+1)N-1 d d d HmN,m3 · · · HmN,mm )U2N|(m+1)N-1 (HmN,m2

(20)

and HdmN,mi designates the (m, i) block of HdmN. Notice that HdmN,mi d d ) H(m-i+1)N,(m-i+1)1 for i ) 2, 3, ..., (m - 1), and HmN,mm ) HNd . (2) The FIR sequences hi, for i ) (m - 1)N, (m - 1)N + 1, d . ..., mN - 1, are obtained in the first column of HmN,m1 In this way, the FIR sequences can be computed piece-bypiece until all of the nonzero FIR parameters are identified. There are two obvious advantages of using the proposed recursive algorithm for FIR estimation. (i) First, it is possible to avoid the numerical problems associated with inversion of a large matrix. To identify the complete FIR model by the traditional subspace method given in section 2.1, the number of rows, N, should be chosen such that the last impulse response coefficient (last element of the first column of HNd ) approaches zero. For multivariable systems with long settling times, numerical problems associated with the inversion of the large matrix, with dimensions (2nu + ny)N × (2nu + ny)N, in eq 14 occur. These numerical problems can be avoided by choosing a small number of N and then identifying the complete FIR model using the recursive algorithm. Consequently, as shown in eq 19, only one computation of low-order matrix inversion is needed in this procedure of FIR estimation. (ii) Second, it is possible to reduce the required length of the experiments. The required length of the experiments for generating data is related to the number of columns in the data Hankel matrices, j, which plays a role similar to that of the number of observations in regression analysis. Thus, the value of j should be chosen sufficient large to obtain a consistent estimate, and typically, j . max(nuN, nyN), as this reduces noise sensitivity. Because the number N can be chosen small in the recursive algorithm, the number j required for accurate identification is thus smaller than that in the traditional subspace identification method (without recursive calculation). Therefore, the use of the recursive algorithm can significantly reduce the required length of the experiments. The following is an example to illustrate the effectiveness of the recursive algorithm. Consider a 2 × 2 process with transfer function matrix as:

[

-3e-0.5s 2e-s 3s + 1 4s + 1 Gp(s) ) -3e-s e-0.5s 2s + 1 2.5s + 1

]

(21)


1463

Figure 1. Comparison of FIR models identified using the methods without and with recursive algorithm (solid line, actual; dashed line, without recursive algorithm; dotted line, with recursive algorithm).

Figure 2. Block diagram of closed-loop control system.

To identify the FIR model, the process is excited simultaneously from both inputs with a pseudorandom binary sequence (PRBS) shifting between 1 and -1. The sampling time is taken as 0.1, and 900 samples are collected. The methods of SIMbased FIR estimation without and with recursive algorithm, where the number of N is chosen as 150 and 30, respectively, are compared. A total number of 150 FIR parameters are identified using both methods. The complexity of computation is greater in the method without recursive algorithm because an inversion of 900 × 900 matrix is needed, while only an inversion of 180 × 180 matrix is needed in the method with recursive algorithm. The results of identification are compared in Figure 1, where the method with recursive algorithm is shown to outperform the method without recursive algorithm. The reason is that the number of available data (j) is not large enough compared with the number of N in the method without recursive algorithm. This result clearly verifies that the use of recursive algorithm can reduce the required length of the experiments for FIR identification. In general, the proposed method can be applied to identify the FIR models of any stable multivariable systems, provided that the input signal is persistently exciting and uncorrelated with noise. In addition, the number of measurements should be taken sufficiently rich, as this reduces the variance errors of estimated parameters. However, this limitation on the required number of data points can be alleviated, to certain extent, using the recursive algorithm.

and noise signals, respectively; and Gp(s), Gd(s), and Gc(s) are the process, disturbance, and controller transfer function matrices, respectively, with appropriate dimensions. The identification problem is to estimate models for Gp(s) and Gd(s) using closed-loop data. Under closed-loop control, the FIR estimation algorithm presented in the previous section using u and y cannot be directly applied because of the correlation between process input u and noise signal e. Therefore, two new closed-loop identification methods are presented in this section. The first method is an indirect approach, where the controller Gc(s) is assumed known. The second one can be classified as a joint input-output approach, where the information of controller Gc(s) is not needed. 3.1. Algorithm 1. In Figure 2, the process output y is given by y ) GpGc(I + GpGc)-1r + (I + GpGc)-1Gdd } Myr(s)r + Myd(s)d (22) where Myr(s) ) GpGc(I + GpGc)-1 is the closed-loop transfer function relating r to y and Myd(s) ) (I + GpGc)-1Gd is the closed-loop transfer function relating d to y. The first step is to estimate the FIR models of Myr(s) and Myd(s). For this purpose, the closed-loop system is excited by simultaneously applying pseudorandom binary sequence (PRBS) at all the reference signals r and disturbance variables d. Because r and d are uncorrelated with the noise e, the FIR models of Myr(s) and Myd(s) (i.e., hyr and hyd, respectively) can be obtained using the algorithm given in section 2 by taking [r d]T as system input and y as system output. A system’s frequency response is the Fourier transform of its impulse response.11 Thus, with the identified hyr, the frequency response of Myr(s) can be computed using discrete Fourier transform (DFT), that is N-1

Myr(ωk) )

∑h

3. Model Identification Using Closed-Loop Data Consider a multi-input, multi-output control system, as illustrated in Figure 2, where u ∈ Rnu is the process input; d ∈ Rnd is the disturbance; r, y, e ∈ Rny are the reference, output,

-j2πik/N

yr,ie

i)0

,

k ) 0, 1, 2, ...,

N 2

(23)

where ωk ) 2πk/(NTs), N is the length of the FIR model, and Ts is the sampling interval. Likewise, the frequency response of Myd(s), namely, Myd(ωk), can be computed using DFT on hyd. The DFT can be efficiently calculated using the fast Fourier

1464


transform (FFT) algorithm. In addition, higher resolution is obtained in the frequency response by padding the FIR with zeros.11 The process transfer function matrix can be written in terms of closed-loop transfer function as Gp ) (I - Myr)-1MyrGc-1, so that its frequency response can be solved using the frequency responses of Myr and Gc by Gp(ωk) ) [I - Myr(ωk)]-1Myr(ωk) Gc(ωk)-1

(24)

at the frequency values of ωk. On the other hand, the disturbance transfer function matrix can be written as Gd ) (I + GpGc)Myd, so that its frequency response can be solved using the frequency responses of Myd, Gp, and Gc by Gd(ωk) ) [I + Gp(ωk) Gc(ωk)]Myd(ωk)

(25)

Consequently, the FIR of the process transfer function (i.e., hp) can be constructed by calculating the inverse FFT (IFFT) of Gp(ωk), that is N/2

hp,i )

∑

1 G (ω )ej2πik/N, N k)0 p k

i ) 0, 1, 2, ..., N - 1

(26) Likewise, the FIR of the disturbance transfer function (i.e., hd) can be constructed by calculating the IFFT of Gd(ωk). In industrial process control, the controllers are usually proportional-integral (PI) or proportional-integral-derivative (PID) to eliminate offset. When the controller contains an integrator, the use of eqs 24 and 25 would, in general, not be reliable in the low-frequency range because the frequency response of Gc(s) approaches infinity. Therefore, another approach is proposed to get around this problem. 3.2. Algorithm 2. In Figure 2, the process input u is given by u ) Gc(I + GpGc)-1r - Gc(I + GpGc)-1Gdd } Mur(s)r + Mud(s)d (27) where Mur(s) ) Gc(I + GpGc)-1 is the closed-loop transfer function relating r to u and Mud(s) ) -Gc(I + GpGc)-1Gd is the closed-loop transfer function relating d to u. With the closedloop data, the FIR models of Mur(s) and Mud(s) (i.e., hur and hud, respectively) can be obtained using the algorithm given in section 2 by taking [r d]T as system input and u as system output. Then, frequency responses of Mur(s) and Mud(s) can be calculated using FFT on hur and hud, respectively. It seems that the frequency responses of Gp(s) and Gd(s) can be calculated by the following equations Gp(ωk) ) Myr(ωk)[Mur(ωk)]-1

(28)

Gd(ωk) ) -[Mur(ωk)]-1Mud(ωk)

(29)

However, depending on the type of controller Gc, Mur(s) might not be strictly proper, which prevents its frequency responses from being correctly calculated using FFT. It is proposed here to use a stable and strictly proper filter F(s), such as a firstorder lag. That is, if the process input u is first allowed to pass through F(s), then the closed-loop transfer function relating r to the filtered process input, uF, is given by Mur,F(s) ) Gc(I + GpGc)-1F. The FIR model of Mur,F(s) (i.e., hur,F) can be identified using the algorithm given in section 2 by taking r as system input and the filtered process input, uF, as system output.

Because F(s) can be chosen to make Mur,F(s) strictly proper, the frequency response of Mur,F(s) can be correctly calculated using FFT on hur,F. As a result, the frequency response of Gp(s) can be solved as Gp(ωk) ) Myr(ωk)[Mur,F(ωk)]-1F(ωk)

(30)

Also, the frequency response of Gd(s) can be solved as Gd(ωk) ) -[Mur,F(ωk)]-1Mud(ωk) F(ωk)

(31)

Consequently, the FIR models of the process and disturbance transfer function can be constructed by calculating the IFFTs of Gp(ωk) and Gd(ωk), respectively. In this algorithm, knowledge of the controller is not needed, so the numerical problem encountered in algorithm 1 is avoided. A schematic diagram for the proposed identification procedures is as shown in Figure 3. Remark 1. The statistical properties for the subspace-based identification method have been studied in the literature (e.g., Van Overschee and De Moor12), where it has been shown that subspace identification can yield a consistent estimation of the parameters under open-loop conditions. Because the proposed algorithms utilize the subspace identification method in an openloop manner, the property of consistency can be applied. Remark 2. With the identified FIR models, one can directly reconstruct the step response models, Sp and Sd, which are very suitable for designing model predictive controllers (MPC), such as dynamic matrix controller (DMC). If parametric transfer function models are preferred, several algorithms reported in the literature5,6 can be applied to identify parametric models from the step responses. Remark 3. Although the presented methods for identification of disturbance models are theoretically correct, it should be noted that the disturbance models can be identified only in limited circumstances. Practically, it is possible to identify the disturbance models in cases where the measurements of corresponding disturbance variables, which are uncorrelated with the set-point change, are available. For example, the measured disturbances are available for feedforward control. 4. Applications to Monitoring and Retuning of Controllers This section presents applications of the proposed method to control system monitoring and controller retuning. These applications can be accomplished using the information obtained from the intermediate step of the proposed identification procedures. 4.1. Monitoring of Maximum Closed-Loop Log Modulus. The goal of control system monitoring is to provide information that can be used to assess the current status of the existing controller, which can assist process engineers in determining whether a redesign of the controller is necessary. The maximum closed-loop log modulus, Lc,max, provides a good measure of the compromise made between performance and robustness in the frequency domain for control systems and, thus, can indicate the appropriateness of the controller parameters to the changing process dynamics. In the literature, online identification of Lc,max using a relay feedback test for SISO13 and MIMO14,15 control system has been proposed. However, these methods require several relay experiments, where a series of relays has to be placed between the process and controller and the number of relay experiments increases exponentially with the dimension of the systems, for example, three relay tests


1465

Figure 3. Schematic diagram for the proposed identification procedure and its applications.

for a 2 × 2 system and seven relay tests for a 3 × 3 system.15 Therefore, it can be time-consuming to apply these methods in practice because of the tedious procedure involved. Here, a novel method developed from the proposed identification procedure is presented. For multivariable control systems, the maximum closed-loop log modulus is defined as16

( |

Lc,max ) 20 log max ω

-1 + det(I + GpGc) det(I + GpGc)

|)

(32)

From the definition of Myr, it is straightforward to show that I + GpGc ) (I - Myr)-1

(33)

PI/PID controller because this method can bypass the tedious iterative procedures and can give a satisfactory result. Most importantly, the proposed identification approach is very suitable for incorporation with the design method. In Huang et al.,18 the design of multiloop controllers is decomposed into tasks of design for controllers in a number of equivalent and independent single loops. The transfer function that describes the effective transmission in each equivalent loop is considered as the effective open-loop process (EOP) of that loop. For example, the EOPs in a 2 × 2 system (gj1 and gj2) can be expressed as gj1(s) ) g11 - g12g22-1g21C2

gj2(s) ) g22 - g21g11-1g12C1

(35)

By combining eqs 32 and 33, Lc,max can be written as

( |

Lc,max ) 20 log max ω

-1 + det[(I - Myr)-1] det[(I - Myr)-1]

|)

(34)

Because the frequency response of Myr(s) has been identified in the proposed algorithm (see eq 23), Lc,max can be directly calculated using eq 34. Notice that the identification of Lc,max requires no a priori information of the process and controller. A large Lc,max value means that the control system can tolerate a small amount of uncertainties. For control system monitoring, the identified Lc,max from the proposed method can be compared with +2n dB, which is the design criterion used in the multiloop control system by Luyben.17 4.2. Retuning of Multiloop Controllers. Despite the recent advances in multivariable control design methods, the control structure used in the chemical process industries is still dominated by multiloop PI/PID controllers. There is a need to redesign the controller when the identified Lc,max from the proposed method indicates that the current status of the existing control system is not appropriate. In this article, a direct method proposed by Huang et al.18 is adopted to redesign the multiloop

where gij designates the element of process transfer function matrix and Ci designates the single-loop complementary sensitivity function of the ith loop. The difficulty of modeling each EOP, due to the needs of the controller dynamics in other loops, is overcome by a benchmark formulation of Ci as 0.6 (1 + 0.4θis) -θis e θi s Ci ) 0.6 (1 + 0.4θis) -θis e 1+ θi s

(36)

where θi denotes the time delay of the ith diagonal open-loop process. With these formulated EOPs, design of controllers can be carried out directly and independently without referring to the controller dynamics of other loops. Model-based tuning formulas for multiloop PI/PID parameters are then formulated in terms of simple parametric models, or in terms of the ultimate gain and ultimate frequency, of these EOPs. In the work of Huang et al.,18 the process transfer function matrix is assumed to be known to compute the frequency

1466


Figure 4. Closed-loop data for example 1.

responses of EOPs by eq 35. However, this process transfer function matrix might not be available in practice. When the process transfer function matrix is unknown, the proposed identification procedure is readily applicable because the frequency response of process transfer function is identified by eq 24 or eq 30 and the time delay θi can be extract from the identified FIR model of hyr. With these pieces of information, the frequency response of each EOP can be computed using eq 35, and then a simple parametric model of each EOP can be found to have a best fit to its frequency response. Consequently, the multiloop PI/PID controller can be redesigned using the tuning formulas given in Huang et al.18 A schematic diagram for the applications of the proposed identification procedure is shown in Figure 3.

method. In practice, the measurements are generally corrupted by noises. In the context of system identification, noise-to-signal ratio (NSR) is usually defined as5 NSR )

mean[abs(noise)] mean[abs(signal)]

(37)

The identification accuracy is evaluated by the sum of the squared errors (SSE) between the identified and the actual FIR models. 5.1. Example 1: 2 × 2 Process. Consider the closed-loop system of Figure 2 with the following process, disturbance, and controller transfer function matrices

[

] [ ]

-3e-0.5s 2e-s 3s + 1 4s + 1 Gp(s) ) , -3e-s e-0.5s 2s + 1 2.5s + 1

5. Simulation Examples The proposed identification methods are now applied to typical 2 × 2 and 3 × 3 processes, including a nonlinear stirred mixing tank process. Examples 1 and 2 are used to demonstrate the results of identification, whereas examples 3 and 4 are mainly used to demonstrate the applications of the proposed

e-2s 2.8s + 1 Gd(s) ) , 1.5e-s 2.4s + 1 1 1+ 0 2s Gc(s) ) 1 -0.1 1 + 0 s

[

(

)

]

(38)

Table 1. Calculated SSEs of Two Proposed Algorithms under Different Values of NSR for Example 1a 0% NSR Gp11 Gp12 Gp21 Gp22 Gd11 Gd21 a

A1

20% A2

-4

1.65 × 10 3.92 × 10-4 7.05 × 10-5 1.85 × 10-4 5.99 × 10-5 2.96 × 10-5

A1 -7

2.51 × 10 8.88 × 10-7 7.06 × 10-7 1.68 × 10-6 1.09 × 10-7 3.70 × 10-6

A1, algorithm 1; A2, algorithm 2.

0.0015 0.0058 4.76 × 10-4 0.0027 7.56 × 10-4 3.57 × 10-4

50% A2

A1 -5

8.40 × 10 6.08 × 10-4 8.07 × 10-5 7.36 × 10-4 1.22 × 10-5 2.66 × 10-5

0.0026 0.0110 0.0011 0.0072 0.0012 6.55 × 10-4

100% A2 -4

4.00 × 10 0.0034 4.89 × 10-4 0.0044 6.68 × 10-5 1.55 × 10-4

A1

A2

0.0032 0.0196 0.0024 0.0196 0.0013 0.0011

0.0015 0.0132 0.0019 0.0174 2.53 × 10-4 6.14 × 10-4


1467

Figure 5. Identified FIR models of process and disturbance dynamics under 20% NSR for example 1 (solid line, actual; dashed line, identified by algorithm 1; dotted line, identified by algorithm 2).

Two proposed closed-loop identification algorithms are used to estimate models of process and disturbance dynamics under different values of NSR. The system is excited simultaneously from r and d with a PRBS shifting between 1 and -1. The sampling time is taken as 0.1, and the closed-loop data are plotted in Figure 4. The filter F(s) used in algorithm 2 is 1/(5s + 1). The FIR models of closed-loop systems are estimated using the proposed recursive algorithm 10 times, each time for 30 parameters. The calculated SSE of the identified results is given in Table 1, and the identified FIR models of process and disturbance dynamics under 20% NSR are shown in Figure 5. It can be seen that the identification results are satisfactory in the presence of noises. Moreover, the accuracy of algorithm 2 is better than that of algorithm 1, as algorithm 1 cannot identify the steady-state gain very well in this example. The reason is that, as mentioned in section 3, the frequency response of the process cannot be accurately computed by algorithm 1 in the low-frequency range due to the controller containing an integrator. Figure 6 presents comparisons between the actual and the identified frequency responses (Nyquist plots) of Gp(s) and Gd(s), which verifies this conclusion. 5.2. Example 2: 3 × 3 Process. Consider the 3 × 3 process of a pilot-scale distillation column for the separation of ethanol and water as follows19

[

Gp(s) ) -0.61e-3.5s -00049e-s 0.66e-2.6s 6.7s + 1 8.64s + 1 9.06s + 1 -2.36e-3s -0.012e-1.2s 1.11e-6.5s 3.25s + 1 5s + 1 7.09s + 1 0.87(11.61s + 1)e-s -33.68e-9.2s 46.2e-9.4s 8.15s + 1 10.9s + 1 (3.89s + 1)(18.8s + 1)

]

(39)

The process is controlled by a diagonal PI controller with the parameters Kc1 ) 1.2, τI1 ) 6.3; Kc2 ) -0.001, τI2 ) 0.1; and Kc3 ) 0.005, τI3 ) 0.1. Closed-loop data are collected with

a sampling time of 1 time unit. Only algorithm 2 is applied in this example because algorithm 2 can perform quite well when the controller has integral mode. The calculated SSE of identified results is given in Table 2, and the identified process step responses under 20% NSR are shown in Figure 7. It is illustrated that the identified models match the actual ones well. The step responses of Gp13 and Gp23 show a little oscillatory behavior because the magnitude of responses for Gp13 and Gp23 are much smaller than the others, so that the effect of noise becomes significant for these two models. 5.3. Example 3: 2 × 2 Wood and Berry Column. To demonstrate the proposed applications of the identification method, the following process transfer function matrix of Wood and Berry distillation column20 for separating methanol from water is considered

[

12.8e-s -18.9e-3s 16.7s + 1 21s + 1 Gp(s) ) -19.4e-3s 6.6e-7s 10.9s + 1 14.4s + 1

]

(40)

The above process is controlled by decentralized PI controllers with Kc1 ) 0.2, τI1 ) 4.44 and Kc2 ) -0.04, τI1 ) 2.67. To monitor the appropriateness of the current controller, the proposed procedure is applied to identify Lc,max of the current control system using simulated closed-loop data with 20% NSR. The obtained Lc,max value is 9.96 dB, which is very close to the actual value of 10.1. According to this result, it is concluded that the current control system lacks robustness, and a redesign of the controller is necessary. Hence, the proposed procedure is again applied to retune the decentralized PI controllers. Based on the identified process frequency responses, the frequency responses of two EOPs are computed, and the following two parametric models for them are further obtained gj1(s) )

6.90e-1.28s , 7.07s + 1

gj2(s) )

-9.06e-3.93s 3.92s + 1

(41)

1468


Figure 6. Identified frequency responses of process and disturbance dynamics under 20% NSR for example 1 (solid line, actual; dashed line, identified by algorithm 1; dotted line, identified by algorithm 2).

Figure 7. Identified process step responses under 20% NSR for example 2 (solid line, actual; dashed line, identified by algorithm 2).

Table 2. Calculated SSEs of Algorithm 2 under Different Values of NSR for Example 2 NSR

0%

10%

15%

20%

50%

Gp11 Gp12 Gp13 Gp21 Gp22 Gp23 Gp31 Gp32 Gp33

3.18 × 10-5 5.43 × 10-4 4.52 × 10-8 2.54 × 10-4 0.0300 4.28 × 10-7 0.0763 2.2639 0.0025

4.49 × 10-5 0.0013 8.72 × 10-7 5.46 × 10-4 0.0477 2.22 × 10-5 0.1183 5.3982 0.0056

4.42 × 10-5 0.0026 1.95 × 10-6 5.69 × 10-4 0.0855 5.31 × 10-5 0.1480 13.4986 0.0136

5.82 × 10-5 0.0052 3.65 × 10-6 8.38 × 10-4 0.1473 1.00 × 10-4 0.2001 25.6495 0.0229

2.39 × 10-4 0.0238 2.32 × 10-5 0.0051 0.6466 6.34 × 10-4 0.9863 123.699 0.1189


1469

Figure 8. Set-point responses of the original and retuned control system for example 3.

By making use of these parametric models, the PI controller parameters are calculated as Kc1 ) 0.48, τI1 ) 7.07 and Kc2 ) -0.066, τI2 ) 3.92. Figure 8 compares the set-point responses for the original and retuned system. Clearly, the original control system shows a more oscillatory response and the performance is greatly improved after retuning the controllers. This result verifies the effectiveness of the proposed procedures for monitoring and controller redesign. 5.4. Example 4: Nonlinear Stirred Mixing Tank Process. A stirred mixing tank21 as depicted in Figure 9 is considered. In this system, the liquid level, h, and the temperature of the liquid in the tank, T, are controlled by manipulating the hot-stream flow rate, qh, and the cold-stream flow rate, qc. The following nonlinear ordinary differential equations describe the dynamics of the mixing tank Ac FcpAc

dh ) qh + qc - φ√h dt

d(hT) ) Fcp(qhTh + qcTc - φ√hT) dt

The proposed procedure is employed to monitor the robustness of the current control system. Notice that the recursive FIR estimation is not applied in this nonlinear system because the problem of convergence might be encountered. The identified Lc,max value is found to be 13.85 dB. Therefore, it is possible to redesign the controller to make the system more robust. Using the proposed procedure, the frequency responses of two EOPs are computed, and the following two parametric models resulted gj1(s) )

0.59e-0.56s , 0.69s + 1

gj2(s) )

(44)

As a result, the multiloop PI controller parameters are calculated as Kc1 ) 1.25, τI1 ) 0.69 and Kc2 ) -0.25, τI2 ) 0.53. The

(42) (43)

where Ac is the cross-sectional area of the tank, F is the liquid density, cp is the liquid heat capacity, φ is a proportional constant, Th is the hot-stream temperature, and Tc is the coldstream temperature. Here, the hot-stream temperature, Th, is considered as a measured disturbance. The process conditions studied in the simulation are as follows: Th ) 65 °C, Tc ) 15 °C, Ac ) 1 m2, and φ ) 10 m5/2min-1. Multiloop PI controllers are used to control the liquid level and temperature at the set points of 4 m and 50 °C, respectively. The controller parameters are Kc1 ) 0.8, τI1 ) 0.3 in the level loop and Kc2 ) -0.1, τI2 ) 0.15 in the temperature loop. The measurement time delays in the level loop and temperature loop are both 0.5 min. A sampling interval of 0.1 min is used in the simulation, and all measurements are corrupted by 20% noise.

-2.53e-0.50s 0.53s + 1

Figure 9. Diagram of a stirred mixing tank.

1470


Figure 10. Closed-loop responses for a -5 °C step disturbance in hot-stream temperature for example 4.

Figure 11. Identified process and disturbance step responses under 20% NSR for example 4 (solid line, actual; dashed line, identified by algorithm 2).

Lc,max of the retuned system is identified as 5.11 dB. Figure 10 compares the load responses of the original and retuned system for a -5 °C step change in hot stream temperature. Clearly, the original control system shows a more oscillatory response due to a large value of Lc,max, and the retuned system gives a better response. Finally, the FIR models of the process and disturbance dynamics are identified from closed-loop data, and these

identified open-loop step responses are compared with the actual ones in Figure 11. Because a change in the hot-stream temperature has no effect on the liquid level under open-loop conditions, the corresponding FIR model identified approaches zero. The result shows that the proposed identification method also works well for a nonlinear system. Here, the proposed method is compared with a subspace-identification-based method presented by Kadali and Huang10 that identifies the process


1471

Figure 12. Comparison of process step responses identified using the proposed method and the method in Kadali and Huang10 for example 4 (solid line, actual; dashed line, identified by algorithm 2; dotted line, identified by method in Kadali and Huang10).

dynamic matrix from closed-loop data. The same data are used for the identification in both methods, and the results are compared in Figure 12. It is shown that the proposed algorithm leads to more accurate models than the method of Kadali and Huang.10

controller monitoring, and controller retuning, can be combined to develop an intelligent control system.

By making use of these identified FIR (or step response) models of process and disturbance dynamics, a model predictive controller (MPC) can be implemented to achieve better control performance. The responses of the MPC are shown in Figure 10. The MPC parameters are tuned as follows: prediction horizon ) 20, control horizon ) 5, and equal weighting on all inputs and outputs. It can be seen that the MPC outperforms the multiloop PI controllers because the disturbance model is used in MPC for prediction, which verifies the appropriateness of the identified models.

The author thanks the National Science Council (NSC) of Taiwan for supporting this research under Grant NSC-98-2221E-027-018.

6. Conclusions This article presents two novel closed-loop identification methods, an indirect approach and a joint input-output approach, for multivariable systems. An efficient algorithm has been derived from subspace identification methods to first estimate the FIR models of closed-loop systems in a recursive manner. The open-loop models are then retrieved through the frequency domain by FFT and IFFT techniques. Models of both process and disturbance dynamics can be identified from one closed-loop test. Parameterizations of the models and iterative or nonlinear optimization schemes are not needed. Simulation results have demonstrated the effectiveness of the proposed identification methods even if the measurements are corrupted by noise. It is also found that the proposed joint input-output approach (algorithm 2) can lead to more accurate results than the indirect approach (algorithm 1) when the controller has integral mode. Moreover, practical applications of the proposed procedure for monitoring Lc,max and controller redesign are also provided. These applications would be very useful for improving the performance of existing underperforming controllers. Therefore, the proposed methods, including model identification,

Acknowledgment

Literature Cited (1) Ljung, L. System Identification: Theory for the User; Prentice Hall: Upper Saddle River, NJ, 1999. (2) Forssell, U.; Ljung, L. Closed-Loop Identification Revisited. Automatica 1999, 35, 1215–1241. (3) Qin, S. J. An Overview of Subspace Identification. Comput. Chem. Eng. 2006, 30, 1502–1513. (4) Huang, B.; Kadali, R. Dynamic Modeling, PredictiVe Control and Performance Monitoring; Springer: London, 2008. (5) Wang, Q. G.; Zhang, Y.; Guo, X. Robust Closed-Loop Identification with Application to Auto-tuning. J. Process Control 2001, 11, 519–530. (6) Silva, R.; Sbarbaro, D.; Barra, B. A. L. Closed-Loop Process Identification under PI Control: A Time-Domain Approach. Ind. Eng. Chem. Res. 2006, 45, 4671–4678. (7) Park, H. I.; Sung, S. W.; Lee, I. B.; Lee, J. On-line Process Identification Using the Laguerre Series for Automatic Tuning of the Proportional-Integral-Derivative Controller. Ind. Eng. Chem. Res. 1997, 36, 101–111. (8) Mei, H.; Li, S.; Cai, W. J.; Xiong, Q. Decentralized Closed-Loop Parameter Identification for Multivariable Processes from Step Responses. Math. Comput. Simul. 2005, 68, 171–192. (9) Hsia, T. C. System Identification: Least-Squares Methods; D.C. Heath and Company: Lexington, MA, 1977. (10) Kadali, R.; Huang, B. Estimation of the Dynamic Matrix and Noise Model for Model Predictive Control Using Closed-Loop Data. Ind. Eng. Chem. Res. 2002, 41, 842–852. (11) Smith, S. W. The Scientist and Engineer’s Guide to Digital Signal Processing; California Technical Publishing: San Diego, CA, 1997. (12) Van Overschee, P.; De Moor, B. Subspace Identification for Linear Systems: Theory Implementation Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (13) Chiang, R. C.; Yu, C. C. Monitoring Prodcedure for Intelligent Control: On-Line Identification of Maximum Closed-Loop Log Modulus. Ind. Eng. Chem. Res. 1993, 32, 90–99.

1472


(14) Ju, J.; Chiu, M. S. Relay-Based On-Line Monitoring Procedures for 2 × 2 and 3 × 3 Multiloop Control Systems. Ind. Eng. Chem. Res. 1997, 36, 2225–2230. (15) Ju, J.; Chiu, M. S. A Fast Fourier Transform Approach for OnLine Monitoring of the Maximum Closed-Loop Log Modulus. Ind. Eng. Chem. Res. 1998, 37, 1045–1050. (16) Luyben, W. L. Process Modeling, Simulation, and Control for Chemical Engineers; McGraw-Hill: New York, 1990. (17) Luyben, W. L. Simple Method for SISO Controllers in Multivariable Systems. Ind. Eng. Chem. Process Des. DeV. 1986, 25, 654. (18) Huang, H. P.; Jeng, J. C.; Chiang, C. H.; Pan, W. A Direct Method for Multi-loop PI/PID Controller Design. J. Process Control 2003, 13, 769– 786.

(19) Ogunnaike, B. A.; Ray, W. H. Process Dynamics, Modeling, and Control; Oxford University Press: New York, 1994. (20) Wood, R. K.; Berry, M. W. Terminal Composition Control of a Binary Distillation Column. Chem. Eng. Sci. 1973, 28, 1707–1717. (21) Wang, B.; Chiu, M. S. On-Line Monitoring for Multiloop Control Systems Under Load Disturbance. Trans. Inst. Chem. Eng. 2000, 78, 605– 611.

ReceiVed for reView May 4, 2010 ReVised manuscript receiVed September 17, 2010 Accepted September 20, 2010 IE1010146

Closed-Loop Identification of Dynamic Models for Multivariable

Recommend Documents