Identication of multivariable delay processes in presence of nonzero initial conditions and disturbances.

Liu, Min ; Wang, Qing-Guo ; Hang, Chang Chieh 等

INTRODUCTION

Identification and control of single variable processes have been well studied (Astrom and Hagglund, 1995; Ljung, 1999). In the past decades, many identification methods have been proposed (Gawthrop, 1984; Unbehauen and Rao, 1987; Wang and Gawthrop, 2001; Garnier et al., 2003). However, most industrial processes are multivariable in nature (Ogunnaike and Ray, 1994). To achieve performance requirement, modern advanced controllers based on process models are implemented (Ikonen and Najim, 2002) and identifications of multivariable processes are in great demand (Cott, 1995; Zhu, 1998). In the context of continuous process identification, many methods have been proposed for multivariable case, for example, Whitfield and Messali (1987), Wang and Zhang (2001a) and Garnier et al. (2007). An important issue with continuous process identification is time delay. Its estimation needs special attention.

Typical test signals used in process identifications include step, pulse, relay, pseudo-random binary sequence (PRBS) and sinusoidal functions. Among them, step, pulse and relay experiments are more popular for their simplicity. In Wang and Zhang (2001a), relay tests are applied. The frequency responses from the inputs to the outputs are obtained by applying the FFT. The process step response is constructed by using the inverse FFT to each process channel. Integral identification methods are then used to recover all the process model parameters including time delay. Their method is very robust in the face of noise. However, their identification tests and those used in Wang et al. (2003) require zero initial conditions and no significant disturbance. For easy applications, these assumptions should be removed. Recently, based on novel inte gration techniques, robust identification methods have been proposed for single variable time delay processes in the presence of nonzero initial conditions and dynamic disturbances (Hwang and Lai, 2004; Wang et al., 2005; Ahmed et al., 2006; Wang et al., 2006; Liu et al., 2007). In Hwang and Lai (2004) and Wang et al. (2005), identifications from pulse tests were proposed. In Wang et al. (2006), relay-tests-based identification is presented. In Ahmed et al. (2006), an iterative method based on a linear filter is proposed. An improved general method was developed in Liu et al. (2007). Extending these SISO identification methods to MIMO cases is of great interest and value.

In this paper, an integral identification method is presented for multivariable processes with multiple time delays. It adopts the integral technique and can work under nonzero initial con ditions and dynamic disturbances. The unknown multiple time delays can be identified without iteration. The effectiveness of the proposed method is demonstrated through simulation and real-time implementation. This paper is organized as follows. In the second section the identification method is developed for two-input and two-output (TITO) time delay processes. Simulation examples are given in the third section. The proposed method is extended to the general cases in the fourth section. In the fifth section the proposed method is applied to a physical thermal control system. Conclusions are drawn in the sixth section.

TITO PROCESSES

To introduce our method with simplicity and clarity, let us consider a TITO continuous-time delay process first,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where [Y.sub.1](s) and [Y.sub.2](s) are the Laplace transforms of two outputs, [y.sub.1](t) and [y.sub.2](t), U1(s) and [U.sub.2](s) are the Laplace transforms of two inputs, [u.sub.1](t) and [u.sub.2](t), and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] and j =1, 2. The given TITO process may be decomposed into 2 two input and single-output sub-processes, which can be described as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Let the common denominator of [G.sub.i1] and [G.sub.i2] [[beta].sub.i](s). We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

The equivalent differential equations are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (1)

where [w.sub.i](t) accounts for the unknown disturbances and biases. Our task is to identify [a.sub.i,k], [b.sub.ij,k] and [d.sub.ij] from some tests on the process. During the identification test, two separate sets of piecewise step signals are applied on two inputs at t = 0, respectively. The test signals under consideration are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where 1(t) is the unit step, [K.sub.1] [greater than or equal to]1 and [t.sub.1,k], k =1,...,K1 are the switching time instants of [u.sub.1](t), and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where [K.sub.2] [greater than or equal to]1 and [t.sub.2,k], k =1,...,[K.sub.2] are the switching time instants of [u.sub.2](t). Such forms of [u.sub.i], i =1, 2, cover many types of test signals such as steps, rectangular pulses, rectangular doublet pulses, PRBS signals and the relay feedback output.

To eliminate those derivatives in (1), we introduce a multiple integration operator,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (2)

Integrating (1) with (2) ni times yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (3)

Its left-hand side is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (4)

where the last term corresponds to the initial conditions of the output. In the right-hand side, it follows that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Suppose that there holds

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (5)

where qi is an integer. Equation (5) covers a wide range of disturbances (Hwang and Lai, 2004) with its simplest as the static disturbance for which [w.sub.i](t)= [c.sub.l](t), [P.sub.ni] [w.sub.i](t)= [ct.sup.ni]/[n.sub.i]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] and [q.sub.i]=[n.sub.i]

Equation (3) is then cast into the following regression linear in a new parameterization:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where [[gamma].sub.i](t) = [[gamma].sub.i](t),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (6)

The first ni elements in ?i are the model parameter ai,k:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (7)

[[theta].sub.i,k], k = [n.sub.i] + 1,...,[2n.sub.i] + 1 are functions of [d.sub.i1] and [b.sub.i1,k],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (8)

are functions of [d.sub.i2] and [b.sub.i2,k],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (9)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

account for the collective effects of the initial conditions and the disturbances. Note that all the elements in ?i(t) should be mutually independent over the real number field to enable identifiability of the parameter vector, [[theta].sub.i]. This is not the case if [t.sub.1,k] = [t.sub.2,k] for all k, for which

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

p = 0,...,[n.sub.i], become dependent of each other. This should be avoided by the identification test design.

One invokes (6) for t = [t.sub.0],...,[t.sub.N], to get

[[PSI].sub.i][[theta].sub.i] = [[gamma].sub.i] (10)

where [[PSI].sub.i] = [[theta].sub.i]([t.sub.0]),...,[theta]([t.sub.N])]T and GAMMA]i = [[gamma].sub.i]([t.sub.0]),...,[[gamma]([t.sub.N])].sup.T]. The ordinary least-squares method can be applied to find the solution

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

In the presence of noise in the measurement of the process output, the instrumental variable (IV) method is adopted to guarantee the identification consistency. For our case, the instrumental variable [Z.sub.i](t) is chosen as

It should be pointed out that for a selected t, the value of some elements of [[theta].sub.i] depend on [d.sub.i1], [d.sub.i2], which are to be identified and unknown. It is possible to estimate a range of [d.sub.i1] and [d.sub.i2] (Hwang and Lai, 2004). In many engineering applications, one can have simple reliable and probably conservative estimation of the range of time delay from knowledge of the process. For example, the range of transportation delay due to a long pipe can be easily estimated based on the pipe length and fluid speed range. Besides, one may start with a rough estimated delay range and use the proposed method to find [d.sub.i1] and [d.sub.i2], estimates of [d.sub.i1] and [d.sub.i2]. Then with [d.sub.-i1] and [d.sub.i2], one retunes the ranges of time delays and apply the proposed method again to achieve a better estimation. Let [d.sub.i1] and [d.sub.i2] be in the ranges of [[d.sub.i1],[d.sub.i1]] and [[d.sub.i2],[d.sub.i2]], respectively. Define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where Tend is the ending time of the identification test. Then, t should be taken in the set of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

to apply (10). There is no need to solve the estimation equation for each delay within the estimated range. Once the estimate ranges of time delays are given, time delays can be obtained by solving some polynomial equations without iteration. Then, all other parameters than delays are determined accordingly.

Once [[theta].sub.i] is estimated by applying the least-squares method or IV method, the model parameters can be recovered. From (8) for k = [2n.sub.i] + [1-m.sub.i1],..., [2.sub.ni] + 1, [b.sub.i1,k], k = 0,...,[m.sub.i1] can be expressed as the functions of [d.sub.i1] and [[theta].dub.i,k],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Substitute [b.sub.i1,k], k = 0,...,[m.sub.i1] into (8) for k = [2.sub.ni] + [1-m.sub.i1], and we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Equation (12) is solved to get [d.sub.i1] and [b.sub.i1],k, k = 0,...,[m.sub.i1] are then obtained from (11). Similarly, we can find [d.sub.i2] from the following algebraic equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[b.sub.i2k], k = 0,...,[m.sub.i2], are then calculated as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

The proposed method will lead to [m.sub.ij] + 1 estimates for [d.sub.ij], just like Wang and Zhang (2001b) and Hwang and Lai (2004). By inspecting the lag between the input and output signals, the selection can simply be made. The selection can also be made by virtue of the con sistency between various sets of [b.sub.ij,k] and [d.sub.ij] and those ignored relations (Hwang and Lai, 2004).

SIMULATION STUDIES

Example 1. Consider the well-known Wood-Berry binary distillation column plant:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (13)

The equivalent differential equations are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (14)

Case A

Assume that [w.sub.1](t) = 1(t) and [w.sub.2](t) = 0.51(t) and the identification test starts from nonzero initial conditions: [y.sub.1](0) = -1, [y.sub.1](1)(0) = 1, [y.sub.2](0) = 0.5 and [y.sub.2] (1)(0) = 2. The test12 signals, [u.sub.1](t) and [u.sub.2](t), are both pulse signals, [u.sub.1](t) = 1(t) - 1(t - 60), and [u.sub.2](t) = 1(t) - 1(t - 30).

The process inputs and outputs are shown in Figure 1 and the sampling interval is 0.02. Suppose that 0 [less than or equal to][d.sub.11] 2, 0 [less than or equal to][d.sub.12] [less than or equal to]6. It leads to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (14)

[T.sub.1] and [T.sub.2] have some elements in common and these elements are included in T . In other word, the elements in T are members of both [T.sub.1] and [T.sub.2]. This can be seen clearly in Figure 2. Choose t = [t.sub.0],...,[t.sub.N] in T, [n.sub.1] = 2, [m.sub.11] = [m.sub.12] = 1 and [q.sub.1] = 2. The proposed method leads to two estimates for [d.sub.11]: one is -39.05 and the other is 1.02. The time delay must be positive so that we choose [d.sub.11] = 1.02. The proposed method also leads to two estimates for [d.sub.12]: -29.11 and 3.02. For the same reason, we choose [d.sub.12] = 3.02. The first sub-process is then obtained as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Suppose that 0 less than or equal to] [d.sub.21] [less than or equal to]14, 0 [less than or equal to][d.sub.22] [less than or equal to]6. The proposed method with [n.sub.2] = 2, [m.sub.21] = [m.sub.22] = 1 and [q.sub.2] = 2 leads to the second subprocess as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[FIGURE 1 OMITTED]

The identification error, ERR = {ER[R.sub.ij]}, is measured by the worst case error,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (15)

where [G.sub.ij]([j[omega].sub.k]) and [G.sub.ij]([j[omega].sub.i]) are the estimated frequency response and the actual ones. The Nyqusit curve for a phase ranging from 0 to -[pi]is considered, because this part is the most significant for control design. For this example, the identification error is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

In real applications, numerical integration is employed to calculate the multiple integration of the output and this introduces errors. Better identification results can be obtain by sampling the process response with a small sampling interval. If the sampling interval is 0.2, the proposed method leads to the identification error as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

In this case, the identification result is still acceptable. If the sampling interval is chosen as 1 and the identification error is obtained as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

The identification error is very large. From these simulations, one can find that a small sampling interval leads to good identification results. Generally, chemical processes have a slow response. With the development of computer technologies, the sampling interval can be set very small and enough data can be obtained easily for use in process identification.

Case B

This is the same as Case A except that process outputs are subject to changing disturbances, where [w.sub.1](t) and [w.sub.2](t) are simulated by letting 1(t) pass through the transfer functions of

1/15 + 1s and -3/20s + 1, respectively. The proposed method, with [n.sub.1] = [n.sub.2] = 2, [m.sub.11] = [m.sub.12] = [m.sub.21] = [m.sub.22] = 1 and [q.sub.1] = [q.sub.2] = 3, leads to

[FIGURE 2 OMITTED]

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

The identification error is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Case C

This is the same as Case B except that a white noise is added to corrupt the outputs. The noise-to-signal ratio defined by

NSR = mean (abs (noise))/mean (abs (signal))

denoted [N.sub.1] and

NSR = variance (noise)/variance(signal)

(denoted [N.sub.2]) are used to represent a noise level. Let the outputs be corrupted with noise of [N.sub.1] = 15%, 25% and 40% or N2 = 3%, 7% and 18%, respectively. Suppose that the estimated ranges of time delays are 0.5 [less than or equal to][d.sub.11] [less than or equal to]1.5, 2 [less than or equal to][d.sub.12] [less than or equal to]4, 6 [less than or equal to][d.sub.21] [less than or equal to]9 and 2 [less than or equal to][d.sub.22] [less than or equal to]4. The identified parameters are expressed as the mean and standard deviation of each estimate from 20 noisy simulations and shown in Table 1.

In case of noise, we may also start with rough estimated delay ranges given in Case A and use the proposed method to find [d.sub.ij], estimates of [d.sub.ij]. Then with [d.sub.ij], we retune the ranges of time delays and apply the proposed method again to achieve a better estimation. For example, in case of [N.sub.1] = 15%, one identification test is applied. The proposed method, with 0 [less than or equal to][d.sub.11]] [less than or equal to]2,

0 [less than or equal to][d.sub.12] [less than or equal to]6, 0 [less than or equal to][d.sub.21] [less than or equal to] 14 and 0 [less than or equal to][d.sub.22] [less than or equal to]6, leads to [d.sub.11] = 1.08, [d.sub.12] = 2.88, [d.sub.21] = 7.23 and [d.sub.22] = 3.05, with the identification error of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

We then retune the ranges of the time delays as above and the proposed method leads to a smaller identification error

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Example 2. Consider a TITO system,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

A closed-loop relay feedback is applied in this example. The relay feedback system is shown in Figure 3. The relay unit is described as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (16)

where e(t) and u(t) are the relay input and output, respectively. The relay experiment is applied at t = 0 with u + = 1, u_ = -1, [epsilon]+ = 0.8 and [epsilon] = -0.8 under zero initial conditions and nonzero static disturbances of [w.sub.1] = [w.sub.2] =0.51(t). The process inputs and outputs are shown in Figure 4 and the sampling interval is 0.02. Suppose that 2 [less than or equal to] [d.sub.11] [less than or equal to]3, 2 [less than or equal to] [d.sub.12] [less than or equal to]3,2 [less than or equal to] [d.sub.21] [less than or equal to]3 and 3 [less than or equal to] [d.sub.22] [less than or equal to]4. The proposed method, with [n.sub.1] = [n.sub.2] = 2, [m.sub.11] = [m.sub.12] = [m.sub.21] = [m.sub.22] = 1 and [q.sub.1] = [q.sub.2] = 2, leads to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

with the identification error as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

GENERAL MIMO PROCESSES

The TITO identification method is now extended to a general MIMO process. Consider a process with l inputs and m outputs, Y(s) = G(s)U(s),

where Y(s) = [[Y.sub.1](s) ... [Y.sub.i](s) ... [Y.sub.l](s)]T is the output vector, U(s) = [[U.sub.1](s) ... [U.sub.j](s) ... [U.sub.m](s)]T is the input vector, and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

the process transfer function matrix. The given MIMO process may be decomposed into l sub-processes, which can be described as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Let the common denominator of all [G.sub.ij],j =1,...,m be [[beta].sub.i](s). We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

The equivalent differential equations are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (17)

The inputs under considerations are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where [t.sub.j,k] is the kth switch instant of [u.sub.j](t). Integrating (17) with (2) [n.sub.i] times yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (18)

The left-hand side is (4) again. For the right-hand side, it follows that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII] (19)

Equation (18) can be rearranged as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

where [[gamma].sub.i](t) = [[gamma].sub.i](t),

Note that the first [n.sub.i] elements of [[theta].sub.i] are the same as (7). [[theta].sub.i,k], k = j([n.sub.i] + 1) + 1,...,j([n.sub.i] + 1) + [n.sub.i], and j = 1,...,m, are combinations of the model parameters [b.sub.ij,k], k = 0,...,mij and [d.sub.ij], and are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[[theta].sub.i,k], k =(m + 1)([n.sub.i] + 1), ..., (m + 1)([n.sub.i] + 1) + [q.sub.i] account for the effects of the aforementioned nonzero conditions and the disturbances.

Suppose that [d.sub.ij], j = 1,...,m are in the ranges of [[d.sub.ij], [d.sub.ij]]. Define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Then, t should be taken in the set of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

One invokes (19) for t in T with t = [t.sub.0],[t.sub.1],...,[t.sub.N], and they give [[PSI].sub.i][[THETA].sub.i] = [[GAMMA].sub.i] (21)

where [[PSI].sub.i] = [[[theta].sub.i]([t.sub.0]) ,..., [[theta].sub.i]([t.sub.N])]T and [[GAMMA].sub.i] = [[gamma].sub.i] ([t.sub.0]),..., [[[gamma].sub.i]([t.sub.N])].sup.T]. The ordinary least-square method can be applied to find the solution; in the presence of noise in the measurement of the process output, the instrumental variable (IV) method is adopted to guarantee the identification consistency. Once ?i is estimated by applying the least-squares method or IV method, the model parameters can be recovered. We can recover [d.sub.ij] from [[theta].sub.i,k], k = j([n.sub.i] + 1),...,j([n.sub.i] + 1) + [n.sub.i], using the following algebraic equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

i = 1.....l and j = 1,....m

Once [d.sub.ij] are obtained, the parameter [b.sub.ij,k] are then calculated as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

Example 3. Consider a system in Vasnani (1995)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

Suppose that [w.sub.1](t) = 100 1(t), [w.sub.2](t) = 20 1(t) and [w.sub.3](t) = 100 1(t) and the identification test starts from nonzero initial conditions: [y.sub.1](0) = [y.sub.2](0) = [y.sub.3](0) = 1, [y.sub.1] (1)(0) = [y.sub.2](1)(0) = [y.sub.3] (1)(0) = 0.5 and [y.sub.1] (2)(0) = [y.sub.2] (2)(0) = -0.2. The process inputs and outputs are shown in Figure 5 and the sampling interval is 0.02. Let 0 < [d.sub.11] < 7, 0 < [d.sub.12] < 7, 1 < [d.sub.13] < 6, 0 < [d.sub.21] < 7, 0 < [d.sub.22] < 5, 0 < [d.sub.23] < 6, 2 < [d.sub.31] < 7, 1 < [d.sub.32] <7 and 0 < [d.sub.33] < 7. The proposed method with [n.sub.i] = 3 [m.sub.ij] = 2 and [q.sub.i] = 3, where i = 1, 2, 3 and j = 1, 2, 3, leads to the following MIMO transfer function matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

with the identification error as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

REAL-TIME TESTING

The proposed method is also applied to a temperature chamber system in our lab. The ex periment set-up consists of two parts: a thermal chamber set (which is made by National Instruments Corp. and shown in Figure 6) and a personal computer with data acquisition cards and LabVIEW software. The system has two inputs: one is to control a 12V light with 20W Halogen bulb, the other is to control a 12V fan. The system output is the temperature of the temperature chamber. Extra transport delays are simulated by using LabVIEW soft ware. An identification test is applied at t = 0. The process inputs and the output are given in Figure 7 and the sampling interval is 0.1 second. [u.sub.1](t) in Figure 7 is used to control the fan speed, and [u.sub.2](t) is used to control the light intensity. First, we estimate the range of time delays roughly: 0 [less than or equal to][d.sub.11] [less than or equal to] 0.8 and 0 [less than or equal to][d.sub.12] [less than or equal to]0.8. Applying the proposed method with [n.sub.1] = 2, [m.sub.11] = [m.sub.12] = 1 and [q.sub.1] = 2, the estimated time delays are obtained as [d.sub.11] = 0.555 and [d.sub.12] = 0.354. Based on these estimated time delays, we can detune the ranges of time delay more accurately: 0.3 [less than or equal to][d.sub.11] [less than or equal to]0.7

and 0.2 [less than or equal to][d.sub.12] [less than or equal to]0.6. Applying the proposed identification method again and one obtains the model as follows,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE TO ASCII]

[FIGURE 7 OMITTED]

If the disturbance is static, the initial conditions can then be estimated (Hwang and Lai, 2004). The estimated response of the model and the real one are shown in Figure 7 for comparison. After the identification test, another cross-evaluation test is applied and the actual response is recorded. The model obtained in the above is used to predict the process response. The actual response and the predicted one are then compared in Figure 8. The effectiveness of the proposed method is obvious.

CONCLUSIONS

Most industrial processes are multivariable in nature and have time delays. Implementation of modern advanced controllers such as the model predictive control explicitly makes use of process models. Thus, the identification of multivariable processes with time delay is in great demand. In this paper, an integral identification method has been presented for multivariable processes with multiple time delays. It adopts the integral technique and can work under nonzero initial conditions and dynamic disturbances. The unknown multiple time delays can be obtained with out iteration. The permissible identification tests include all popular tests used in applications. Only a reasonable amount of computations is required. It is shown through simulation and real time implementation that the proposed method can yield accurate and robust identification results.

[FIGURE 8 OMITTED]

ACKNOWLEDGEMENTS

This work was sponsored by the Ministry of Education's AcRF Tier 1 funding, R-263-000-306-112, Singapore.

REFERENCES

Ahmed, S., B. Huang and S. L. Shah, "Parameter and Delay Estimation of Continuous-Time Models Using a Linear Filter," J. Process Control 16, 323-331 (2006).

Astrom, K. J. and T. Hagglund, "PID Controller: Theory, Design, and Tuning," Instrument Society of America, Research Triangle Park, North Carolina (1995).

Cott, B. J., "Introduction to the Process Identification Workshop at the 1992 Canadian Chem ical Engineering Conference," J. Process Control 5(2), 67-69 (1995).

Garnier, H., M. Gilson and A. Richard, "Continuous-Time Model Identification from Sam pled Data, Implementation Issues and Performance Evaluation," Int. J. Control 76(13), 1337-1357 (2003).

Garnier, H., M. Gilson, P. C. Young and E. Huselstein, "An Optimal IV Technique for Identifying Continuous-Time Transfer Function Model of Multiple Input Systems," Control Eng. Practice 46(15), 471-486 (2007).

Gawthrop, P. J., "Parametric Identification of Transient Signals," IMA J. Math. Control Inform. 1, 117-128 (1984).

Hwang, S. H. and S. T. Lai, "Use of Two-Stage Least-Square Algorithms for Identification of Continuous Systems with Time Delay Based on Pulse Responses," Automatica 40, 1561-1568 (2004).

Ikonen, E. and K. Najim, "Advanced Process Identification and Control," Marcel Dekker, New York, U.S.A. (2002).

Liu, M., Q. G. Wang, B. Huang and C. C. Hang, "Improved Identification of Continuous-Time Delay Processes from Piecewise Step Tests," J. Process Control 17, 51-57 (2007).

Ljung, L., "System Idenfication: Theory for the User," 2nd ed., Prentice Hall PTR, Upper Saddle River, NJ, U.S.A. (1999)

Ogunnaike, B. A. and K. H. Ray, "Process Dynamic, Modeling and Control," Oxford University Press, New York, U.S.A. (1994).

Unbehauen, H. and G. P. Rao. "Identification of Continuous Systems," North-Holland, Ams terdam (1987).

Vasnani, U. V., "Towards Relay Feedback Auto-Tuning of Multi-Loop Systems," PhD Thesis, Department of ECE, National University of Singapore, Singapore (1995).

Wang, L. and P. J. Gawthrop, "On the Estimation of Continuous Time Transfer Functions," Int. J. Control 74, 889-904 (2001).

Wang, Q. G., T. H. Lee and C. Lin, "Relay Feedback: Analysis, Identification and Control," Springer-Verlag, London, U.K. (2003)

Wang, Q.-G. and Y. Zhang, "A Novel FFT-Based Robust Multivariable Process Identification," Ind. Eng. Chem. Res. 40, 2485-2494 (2001a).

Wang, Q.-G. and Y. Zhang, "Robust Identification of Continuous Systems with Dead-Time from Step Responses," Automatica 37, 377-390 (2001b).

Wang, Q.-G., M. Liu and C. C. Hang, "Simplified Identification of Time Delay Systems with Non-Zero Initial Conditions from Pulse Tests," Ind. Eng. Chem. Re s. 44(19), 7591-7595 (2005).

Wang, Q.-G., M. Liu, C. C. Hang and W. Tang, "Robust Process Identification from Relay Tests in the Presence of Nonzero Initial Conditions and Disturbance," Ind. Eng. Chem. Re s. 45(12), 4063-4070 (2006).

Whitfield, A. H. and N. Messali, "Integral-Equation Approach to System Identification," Int. J. Control 45, 1431-1445 (1987).

Zhu, Y. C., "Multivariable Process Identification For MPC: The Asymptotic Method and Its Applications," J. Process Control 8(2), 101-115 (1998).

Manuscript received January 4, 2007; revised manuscript received March 15, 2007; accepted for publication April 8, 2007.

Min Liu [1], Qing-Guo Wang [1]*, Chang Chieh Hang [1] and Wei Tang [2]

[1.] Department of Electrical and Computer Engineering, National University of Singapore, Singapore 119260

[2.] Faculty of Paper-Making Engineering, Shaanxi University of Science & Technology, Xianyang 712081, China

* Author to whom correspondence may be addressed. E-mail address: elewqg@nus.edu.sg

Table 1. Estimated model parameters of Example 1

 [N.sub.1] = 15% [N.sub.1] = 25%
 ([N.sub.2] = 3%) ([N.sub.2] = 7%)

[a.sub.1,1] 0.1101 [+ or -] 0.0076 0.1116 [+ or -] 0.0135
[a.sub.1,0] 0.0029 [+ or -] 0.0007 0.0031 [+ or -] 0.0008
[a.sub.11,1] 0.7746 [+ or -] 0.0235 0.7695 [+ or -] 0.0249
[a.sub.11,0] 0.0389 [+ or -] 0.0064 0.0393 [+ or -] 0.0106
[d.sub.11] 1.0232 [+ or -] 0.0801 1.0523 [+ or -] 0.1489
[b.sub.12,1] -0.9117 [+ or -] 0.0342 -0.9045 [+ or -] 0.0334
[b.sub.12,0] -0.0549 [+ or -] 0.0076 -0.0573 [+ or -] 0.0086
[d.sub.12] 3.0499 [+ or -] 0.1254 3.0421 [+ or -] 0.1440
[a.sub.2,1] 0.1554 [+ or -] 0.0097 0.1581 [+ or -] 0.0143
[a.sub.2,0] 0.0060 [+ or -] 0.0005 0.0061 [+ or -] 0.0009
[b.sub.21,1] 0.6066 [+ or -] 0.0345 0.6127 [+ or -] 0.0445
[b.sub.21,0] 0.0403 [+ or -] 0.0056 0.0413 [+ or -] 0.0078
[d.sub.21] 6.9337 [+ or -] 0.2134 6.9397 [+ or -] 0.2045
[b.sub.22,1] -1.3642 [+ or -] 0.0375 -1.3663 [+ or -] 0.0486
[b.sub.22,0] -0.1130 [+ or -] 0.0096 -0.1156 [+ or -] 0.0122
[d.sub.22] 3.0548 [+ or -] 0.0841 3.0678 [+ or -] 0.1103

 [N.sub.1] = 40%
 ([N.sub.2] = 18%)

[a.sub.1,1] 0.1126 [+ or -] 0.0225
[a.sub.1,0] 0.0032 [+ or -] 0.0013
[a.sub.11,1] 0.7284 [+ or -] 0.1765
[a.sub.11,0] 0.0395 [+ or -] 0.0187
[d.sub.11] 0.9909 [+ or -] 0.3228
[b.sub.12,1] -0.9046 [+ or -] 0.0555
[b.sub.12,0] -0.0579 [+ or -] 0.0143
[d.sub.12] 3.0561 [+ or -] 0.2381
[a.sub.2,1] 0.1607 [+ or -] 0.0242
[a.sub.2,0] 0.0063 [+ or -] 0.0015
[b.sub.21,1] 0.6126 [+ or -] 0.0753
[b.sub.21,0] 0.0429 [+ or -] 0.0133
[d.sub.21] 6.9233 [+ or -] 0.3372
[b.sub.22,1] -1.3620 [+ or -] 0.0835
[b.sub.22,0] -0.1180 [+ or -] 0.0206
[d.sub.22] 3.0796 [+ or -] 0.1857