Closed-Loop Identification of Multivariable Systems by Optimization

ARTICLE pubs.acs.org/IECR

Closed-Loop Identification of Multivariable Systems by Optimization Method C. Rajapandiyan and M. Chidambaram* Department of Chemical Engineering, Indian Institute of Technology, Madras, Chennai-600036, India ABSTRACT: An optimization method is presented for the closed-loop identification of first-order-plus-time-delay (FOPTD) transfer function models of multivariable systems using step responses. A standard least-squares optimization method is used to obtain the parameters of the FOPTD models by matching the closed-loop step responses of the model with those of the actual process. A simple method is proposed to obtain the initial guess values for the transfer function model parameters from the process main and interaction responses. The effects of measurement noise and controller settings on the identified model parameters were also studied. This method was applied to stable FOPTD and higher-order transfer function models of multivariable systems. The proposed method considerably reduces the computational time (by about a factor of 15) for the optimization when compared with the genetic algorithm method reported by Viswanathan et al. (Ind. Eng. Chem. Res. 2001, 40, 28182826).

1. INTRODUCTION Most chemical processes are multi-inputmulti-output (MIMO) systems. The main characteristic of multivariable systems is the presence of interactions among the loops. The control of such processes is complicated. Multivariable system transfer function matrix models are required to design decentralized/centralized/ decoupled proportionalintegral (PI)/proportionalintegral derivative (PID) controllers.1 Whereas many studies have been reported for the identification of single-input single-output (SISO) systems, fewer studies have been reported for the identification of multivariable systems. The interactions among the loops makes MIMO identification difficult. There are two methods available for the identification of systems: the open-loop identification method and the closed-loop identification. For the open-loop test, an excitation is introduced in each input variable (manipulated variable) of the process, one at a time, to obtain the output responses. Using these output responses, the transfer function matrix is identified. Because there is no interaction in this test, the identification of multivariable systems is easier. Open-loop identification is simple, but the test is sensitive to disturbances and is not applicable for unstable systems. In the closed-loop identification method, some form of excitation (step/pulse of known magnitude) is introduced in each of the set points, one at a time, and the output responses of the actual process are obtained. The transfer function matrix model is obtained using the main and interaction responses. The effects of interactions cannot be neglected in closed-loop identification. The closed-loop test reduces the effects of disturbances on the estimated model parameters. The model parameters obtained by the closed-loop identification method fit the actual process better than those obtained by open-loop identification.25 Further, the controller designed from the model obtained by closed-loop identification exhibits better performance than the controller designed from the open-loop identification.25 The closed-loop identification method is classified in three types:2 direct method, indirect method, and joint inputoutput method. r 2011 American Chemical Society

In the direct method, the process input (additional measurement) and output data are used, and the identification is carried out similarly to open-loop process identification. The works reported by Sung et al.6 and Sung and Lee7 for SISO systems fall under the direct method. Here, the open-loop model is identified by a simple linear least-squares method. The open-loop model, G(s), is defined by a numerator polynomial (mth order in s) and a denominator polynomial (nth order in s). Then, a suitable model reduction method is applied to obtain an first-order-plus-time-delay (FOPTD) or second-order-plus-time-delay (SOPTD) model. Sung and Lee8 extended the direct method for SISO systems to multi-input and single-output (MISO) systems as well. The model parameters for the deterministic part are identified by a nonlinear least-squares (LevenbergMarquardt) method. Sung and Lee and Sung et al.8,9 suggested that initial guess values of the parameters can be obtained, if possible, by physical insight. Otherwise, a continuous-time autoregressive with exogenous input (ARX) model can be identified first, and then a model reduction method can be used to obtain the guess values. This method introduces overparameterization. Methods for identifying the lower-order model directly are preferable to the above approach. In the indirect identification method, by using the known reference signal to the process set point and the measured output response and using the controller settings, the process model parameters are obtained. Viswanathan and Rangaiah10 presented a method for closed-loop identification of SOPTD models using optimization for single-inputsingle-output (SISO) systems. Input excitation was carried out by a step or relay method. Pramod and Chidambaram11,12 applied an optimization method for the identification of stable and unstable process transfer function models and provided a simple method to obtain the initial guesses for the model parameters. Initial guesses obtained Received: June 28, 2011 Accepted: December 20, 2011 Revised: December 14, 2011 Published: December 20, 2011 1324

dx.doi.org/10.1021/ie2013778 | Ind. Eng. Chem. Res. 2012, 51, 1324–1336

Industrial & Engineering Chemistry Research

ARTICLE

by the proposed method considerably reduced the computational time. Sree and Chidambaram13 extended the above optimization method for the identification of unstable FOPTD model with a zero using the closed-loop step response data. Ramachandiran et al.14 proposed a curve fitting method from step responses for process identification. Both open-loop and closed-loop identification methods were studied. Mclntosh and Mahalec15 presented a method for obtaining the steady-state gains of multivariable systems by introducing a step change in one set point at a time. From the output data, steadystate gains are calculated by matrix operations. The routine operating data of a plant under closed-loop operation can also be used for this method. A frequency response technique for SISO systems16 was extended to MIMO systems by Melo and Friedly.17 They presented a method for the closed-loop identification of multivariable systems in the form of a nonparametric frequency response model. A step change in the set point is introduced in each loop sequentially, and the output transient response is obtained. The Fourier transform is used to obtain the frequency response from the measured time response. This method does not give the transfer function matrix. Semino and Scali18 used a sequential relay-based technique to excite each loop one at a time, with all other loops under closed-loop conditions. Their method is an extension of the ATV (autotune variation) identification method. Ham and Kim19 applied an optimization method for process identification of multivariable system from the pulse response; it can also be used for tuning of PID controllers. A pulse set-point change takes less time than a relay feedback test. Mei et al.20 proposed a decentralized closed-loop identification method for a two-inputtwo-output (TITO) process. The coupled closed-loop

two-inputtwo-output system was decoupled equivalently into four individual single-open-loop systems. The multivariable identification process is converted into a multiple single-open-loop identification process. This method requires process input and output data. Mei and Li21 proposed a method for the identification of multivariable integrating processes based on their previous work.22 The integrating process is an open-loop unstable process. The frequency response matrix can be identified by indirect methods from the closed-loop data. The structural information of the integrating process such as the number of integrating factors and their positions in the transfer function matrix can be identified by calculating the value for integrating indexes. Yagang23 used a signal decomposition and analysis method for identifying a multivariable system from the closed-loop test in the frequency domain. The process frequency response matrix can be formulated over important frequency ranges and solved by linear least-squares method to obtain the model parameters. This method requires process input and output response data. In summary, the majority of existing techniques are based on frequency-domain approaches. Comparatively, it should be observed that time-domain identification methods are easier to formulate and conceptually simpler. However, the optimization methods used for time-domain multivariable identification take considerable time. Viswanathan et al.24 reported that the computational time for the identification of an FOPTD TITO system using a genetic algorithm ranges from 1 h to 3 h and 14 min (for different controller settings). The computational time is higher. They presented the upper and lower bounds for a genetic algorithm based on a trial basis for each simulation example. For the limiting values of process gains, they used a method reported by Papastathopoulou and Luyben25 that also requires the process input data to obtain the steady-state gains. Hence, there is a need to propose a method for obtaining reasonable initial guess values for model parameters. For any optimization method, the selection of initial guess values has a significant effect on the computational time and convergence. In the present work, a simple and generalized method is proposed for obtaining reasonable initial guess values for the FOPTD transfer function model parameters. A method to obtain the upper and lower bounds to be used in the optimization routine is also presented. The present work is an extension of the method proposed by Pramod and Chidambaram11,12 and Padmasree and Chidambaram13 for

Figure 1. Decentralized control system of a TITO process.

Figure 2. Simplified diagrams for TITO systems for (a) eq 8 and (b) eq 10. 1325

dx.doi.org/10.1021/ie2013778 |Ind. Eng. Chem. Res. 2012, 51, 1324–1336


ARTICLE

Figure 3. Closed-loop responses of the process and the identified FOPTD transfer function model for example 1 using the PID controller settings.19 (Model and process response curves overlap.) Responses obtained for step changes in (left) yr1 and (right) yr2.

Table 1. Model Identification Results with and without Noise Using the PID Controller19 model ij

Kpij

τij

θij

IAE

ISE

time (NI)a

Model Identification Results without Noise 11

guess converged

21

guess converged

12

guess

5.0505

6.2500

1.0000

12.8009

16.7031

1.0002

5.2391

8.7500

7.0000

6.6004

10.9020

7.0001

7.2450

6.2500

3.0000

converged

18.9062

21.0104

3.0002

22

guess converged

16.1290 19.4026

7.5000 14.4027

3.0000 3.0001

11

guess

6.41 104

1.17 108

1.60 103

7.24 108

9.10 104

1.16 108

1.60 103

7.05 108

0.1017

6.24 104

0.1810

1.90 103

152 (25)

Model Identification Results with Noise converged 21

guess converged

a

5.0505

5.0000

1.0000

12.5038

16.0582

1.0447

5.6799

7.5000

7.0000

6.4410

10.3176

6.9856

8.6169

6.2500

3.0000

0.0969

3.62 104

converged guess

18.2553 16.1290

20.1191 6.2500

2.9997 3.0000

0.2110

2.40 103

converged

19.1683

14.0428

2.9041

12

guess

22

68 (21)

Time in seconds; number of iterations (NI) in parentheses.

the identification of SISO systems to MIMO systems. The proposed method is expected to give rapid and guaranteed convergence. The standard lsqnonlin in Matlab routine is used to solve the optimization problem. This method is applied to FOPTD and higher-order transfer function models of multivariable systems.

decentralized controller matrix with compatible dimensions, respectively, expressed as 2

g11 6 6 g21 GðSÞ ¼ 6 6l 4 gn1

2. PROPOSED METHOD Consider an n-inputn-output multivariable system. G(S) and GC(S) are the process transfer function matrix and the 1326

g12 ::: g22 3 3 3 ⋱ gn1 3 3 3

3 g1n 7 g2n 7 7 l 7 5 gnn

ð1Þ



ARTICLE

Table 2. Effects of Changing Controller Settings on the Number of Iterations and Computational Time for Example 1 ref (parameters: Kc11, τi11, τd11; Kc22, τi22, τd22) 29

computational NIa

time (s) 94

IAEb

Table 3. Effects of Initial Guess Values for the Model Parameters on the Number of Iterations and Computational Time for Example 1 Using the PID Controller Settings19

ISEb

0.0229 6.51 10

guess variation 6

NIa

IAEb

ISEb

computational time (s) 7

+10%

25

0.004

1.88 10

152

10%

25

0.0047

1.62 107

154

Vu et al. (0.7, 17.27; 0.13, 16.08)

16

Loh et al.30

21

132

0.0046 3.35 107

a

28

257

0.0066 4.40 107

closed-loop response is faster than the open-loop response. In MIMO systems, because of interactions among the loops, the closed-loop response is usually assumed to be slower than the open-loop response. This assumption is made only to obtain the initial guess values of the time constant of the open-loop system for the optimization problem. Hence, the initial guess value of the time constant for the model11,12 is taken as (ts/8). The initial guess values for the main-loop process gains (Kp11 and Kp22) are obtained by the following method. For the lowerorder model systems, the value26 of kpkc,max is taken as 4, and the design value of kpkc,des is 2. Hence, the guess value for kp is taken as 2/kc. Similarly, for the higher-order model systems, the value26 of kpkc,max is 1, and the design value for kpkc,des is 0.5. Hence, the guess value for kp is taken as 0.5/kc. The higher-order model and lower-order model are distinguished by checking the initial dynamics of the main responses of the actual process. If the initial slope of the main response is sluggish (zero slope), then we consider the system as higher-order model. Similarly, if the initial slope of the main response has is steep (nonzero slope), then we consider the system as a lower-order model. The initial guess value of the cross loop-process gains (Kp21 and Kp12) are obtained by the following method. The Laplace transforms of the interaction responses y21(S*) and y12(S*) can be expressed as

(0.868, 3.246; 0.0868, 10.4) Viswanathan et al.24 (0.2, 4.0, 3.0; 0.03, 2.0, 3.0) a

NI, number of iterations. b Sum of main and interaction response values.

2

gc11 6 60 GC ðSÞ ¼ 6 6l 4 0

0::: gc22 3 3 3 ⋱ 03 3 3

3 0 7 0 7 7 l 7 5 gcnn

ð2Þ

The controller parameters can be chosen arbitrarily for multivariable systems, such that the closed-loop system is stable with reasonable responses. Consider a decentralized TITO multivariable system as shown in Figure 1. The process transfer function models are identified by FOPTD model gmij ðSÞ ¼

Kmpij eθmij s ðτmij s þ 1Þ

i ¼ 1, 2; j ¼ 1, 2

ð3Þ

In the present work, a step change of known magnitude is introduced into the set point yr1, with all remaining set points unchanged and all other loops kept under closed-loop operation. From the prescribed step change in set point yr1, we obtain the main response y11 and the interaction response y21. Similarly, the same magnitude of step change is introduced in set point yr2, and we obtain the main response y22 and interaction response y12. The response matrix of the TITO system can be expressed as " # y11 y12 Y ðtÞ ¼ ð4Þ y21 y22 The first column in the response matrix (eq 4) is the responses (main and interaction) obtained for the step change in set point yr1 in the first loop. The second column is the responses (main and interaction) obtained for the step change in set point yr2 in second loop. From these step responses, the initial guess values of the model parameters are obtained by the following method. As stated earlier, in any optimization method, the selection of initial guess values plays a vital role for which a simple method is proposed. The initial guess values of time delay for each transfer function model can be taken as the same as the corresponding closed-loop time delay values. The initial guess value for time constant is obtained from the settling time (ts) of the main and interaction responses. From the settling time (to reach 98% of the final steady-state value and remain within the limit), the closed-loop time constant (τc) can be assumed as ts/4. The time constant of the open-loop system is assumed to be less than the closed-loop time constant. In principle, for SISO systems, the

NI, number of iterations. b Sum of main and interaction response values

y21 ðsÞ ¼ y12 ðsÞ ¼

Z t 0

Z t 0

ð5Þ

ð6Þ

y21 ðtÞets dt y12 ðtÞets dt

To evaluate eqs 5 and 6, we need to assume the value of s*. In the present work, the value of s* is considered as 8/ts. Here, ts is the settling time of the closed-loop response. For t g ts, the value of e8 is close to zero. Then, the integral value does not change further. So, to reduce the computation time of the integral, the value of s* is taken as 8/ts. The upper bound for t is ts. The numerical values for y21(S*) and y12(S*) are obtained by substituting the value of S* into eqs 5 and 6. The closed-loop transfer function for the cross-loops is assumed by considering the interaction from the main loop taken as a disturbance on that of the interacting loop. Some signals are ignored for simplicity, as shown in Figure 2, to avoid gp12 (see Figure 2a) in deriving eq 7 and to avoid gp21 (see Figure 2b) in deriving eq 9. We obtain y21 ðsÞ gc11 ðsÞ g21 ðsÞ ¼ yr1 ðsÞ ½1 þ gc11 ðsÞ g11 ðsÞ½1 þ gc22 ðsÞ g22 ðsÞ ð7Þ Substitute the value for yr1(S*) as 1/s* and also substitute the controller settings and the guess values of the model parameters 1327



ARTICLE

Figure 4. Closed-loop responses of the process and the identified FOPTD transfer function model for example 1 (with noise) using the PID controller settings.19 Responses obtained for step changes in (left) yr1 and (right) yr2.

(time delays, time constants, and main-loop process gains) into eq 7. We obtain Kp21 eθ21 s ð1=sÞKc11 ½1 þ ð1=τi11 sÞ þ τd11 s ðτ21 s þ 1Þ )( ) y21 ðsÞ ¼ ( θ22 s Kp11 eθ11 s Kp e 22 1 þ Kc11 ½1 þ ð1=τi11 sÞ þ τd11 s 1 þ Kc22 ½1 þ ð1=τi22 sÞ þ τd22 s ðτ11 s þ 1Þ ðτ22 s þ 1Þ Similarly, the initial guess value for Kp12 is obtained as y12 ðsÞ gc22 ðsÞ g12 ðsÞ ¼ yr2 ðsÞ ½1 þ gc22 ðsÞ g22 ðsÞ½1 þ gc11 ðsÞ g11 ðsÞ ð9Þ Substitute the value for yr2(S*) as 1/s* and also substitute the controller settings and the guess values of the model parameters (time delays, time constants, and main-loop process gains) in eq 9. We obtain

Similarly, the same magnitude of step change is introduced in set point ymr2, and the main response ym22 and interaction response ym12 are obtained. The closed-loop responses are compared with the actual process main and interaction responses. To obtain the optimized model parameters, the optimization problem is formulated to select the model parameters to minimize the sum of squared errors between the model and the actual process responses. minimize f ¼

ðKmpij , τmij , θmij Þ

Kp12 eθ12 s ð1=sÞKc22 ½1 þ ð1=τi22 sÞ þ τd22 s ðτ 12 s þ 1Þ ) y12 ðsÞ ¼ ( Kp22 eθ22 s 1 þ Kc22 ½1 þ ð1=τi22 s Þ þ τd22 s ðτ22 s þ 1Þ ( ) Kp eθ11 s 1 þ Kc11 ½1 þ ð1=τi11 sÞ þ τd11 s 11 ðτ11 s þ 1Þ

ð8Þ

2

2

n

½ymij ðtn Þ yij ðtn Þ2 Δt ∑ ∑ ∑ i¼1 j¼1 1

i ¼ 1, 2; j ¼ 1, 2 ð11Þ

ð10Þ The guess values for Kp21 and Kp12 are obtained from eqs 8 and 10, respectively. After obtaining all of the guess values of the model, the same magnitude of step change is introduced to the model set point ymr1 with all other remaining set points kept unchanged and all other loops under closed-loop conditions. From the prescribed step change in set point ymr1, the main response ym11 and interaction response ym21 are obtained.

TITO systems have two main responses and two interaction responses. In the objective function, n is the number of data points in the response of the process and the model. First, we obtain the error between the responses of each transfer function model with the associated process response for various times (to form an error vector for each model). Further, we take the square of each of the errors and multiply by a fixed sampling time, Δt. The sum of the elements of each vector gives an integral of squared error (ISE) value. There are four ISE values. The sum of these four ISE values is considered as a scalar objective function (refer to eq 11). The model parameters are 1328



ARTICLE

Table 4. FOPTD Model Identification Results with and without Noise for Example 2 model ij

τij

Kpij

θij

IAE

ISE

time (NI)a

FOPTD Model Identification Results without Noise 11 21

guess

1.2500

1.2500

0.1000

converged

0.5114

0.3554

0.2801

guess

2.2314

1.2500

0.1000

converged 12 22

0.9906

0.2998

0.2247

guess

2.2314

1.2500

0.1000

converged

0.9843

0.3283

0.1754

guess

5.0000

1.2500

0.1000

converged

2.4518

0.6369

0.4000

guess

1.2500

1.0000

0.1000

converged

0.5074

0.3354

0.3000

guess

2.2793

1.0000

0.1000

0.0177

1.18 104

0.0515

7.46 104

0.0228

9.39 104

0.0415

3.92 104

69 (22)

FOPTD Model Identification Results with Noise 11 21

converged 12 22 a

0.9813

0.3316

0.1917

guess

2.4484

1.2500

0.1000

converged

0.9793

0.2371

0.2611

guess converged

5.0000 2.5160

1.0000 0.6608

0.1000 0.4000

0.0492

3.64 104

0.1222

2.00 103

0.0348

2.22 104

0.0771

6.84 104

30 (23)


Figure 5. Closed-loop responses of the process and the identified FOPTD model for example 2. (Model and process response curves overlap.) Responses obtained for step changes in (left) yr1 and (right) yr2.

selected by minimizing this scalar function. In eq 11, the sampling time, Δt, need not be included because it is a constant value. In this present work, the optimization problems are solved by using the routine lsqnonlin in Matlab, and the routine implements the trust-region-reflective algorithm. Although the optimization problem is formulated and solved using the standard routine (lsqnonlin) in Matlab, the convergence and computational time fully depend on the initial guesses. The main focus of the present work is to propose the above simple method to obtain the initial guess values of the model parameters. The computational time is dramatically reduced by the proposed guess values.

In the lsqnonlin routine, the values used for the convergence parameters are TolX = 1.0 103, TolFun = 1.0 106, and number of iterations = 400.The closed-loop multivariable system is simulated using Simulink. The fourth-order RungeKutta method with fixed step size is used for performing the integration of the governing differential equations. All of the computational work for the identification of multivariable systems was performed on a Core 2 Quad (Intel/3.00 GHz) personal computer. To evaluate the proposed method, a known transfer function model matrix was considered with suitable decentralized PI/PID controllers. The closed-loop main and interaction responses of the actual process were obtained by the earlier discussed method. 1329



ARTICLE

From the responses, the model parameters were obtained by the proposed method, and the closed-loop main and interaction responses were compared with the actual system with the same decentralized PI/PID controller settings.

3. SIMULATION EXAMPLES In this section, three simulation examples (2 2 FOPTD system, 2 2 higher-order system, and 3 3 FOPTD system) are considered to show the effectiveness of the proposed method. 3.1. Example 1. Consider the Wood and Berry27 (WB) distillation process. The transfer function matrix of the WB column is given by 2 3 12:8es 18:9e3s 6 16:7s þ 1 21s þ 1 7 7 ð12Þ Gp ðsÞ ¼ 6 4 6:6e7s 19:4e3s 5 10:9s þ 1 14:4s þ 1 The decentralized PID controller settings used for the identification are19Kc11 = 0.396, τi11 = 5.926, τd11 = 1.070; Kc22 = 0.124, τi22 = 7.092, τd22 = 1.849. Here, the FOPTD multivariable process is taken to show the effectiveness of this method by direct comparison of the converged parameters of the model with those of the actual process. If we take a higher-order system, we cannot compare the model parameters as such; rather, a comparison of the process and the Table 5. Effects of Initial Guess Values for the Model Parameters on the Number of Iterations and Computational Time for Example 2

a

guess variation

NIa

IAEb

ISEb

computational time (s)

+10%

13

0.1226

0.0014

48

10%

24

0.2089

0.0028

86


model can be made only by response matching. The closed-loop main and interaction responses were obtained by the previously mentioned method using Simulink, and they are shown in Figure 3. The FOPTD model was assumed, and the initial guess values for the model parameters were obtained by the earliermentioned methods (refer to Table 1 for the guess values). Because the responses showed a finite nonzero initial slope, kpkc,des = 2 was used for calculating the guess values of process gains kp11 and kp22. The optimization routine lsqnonlin was used in Matlab, and the sample time of 0.01 s was considered. For all of the simulation examples, in the optimization routine, we fixed the lower-bound value for the time constant as 1/20th of the guess value and the upper bound as five times the initial guess value. Similarly, for time delays, the lower bound was specified as 1/4th of the guess value, and the upper bound was specified as four times the initial guess value. If the initial guess value of the process gain was calculated as negative, then we fixed the upper bound as 1/20th of the guess value and the lower bound as 20 times the guess value. Similarly, for positive guess values of the process gain, we fixed the lower bound as 1/20th of the guess value and the upper bound as 20 times the guess value. The upper and lower bound values used in the optimization routine were found to give quick convergence and also to reduce the computational time. The initial guess values, final converged parameters, numbers of iterations, computational times, integrals of absolute error (IAE), and integrals of squared error (ISE) are listed in Table 1. The converged model parameters match well with the actual FOPTD system. The closed-loop main responses and interaction responses using the identified model parameters with the same controller settings are also shown in Figure 3. Figure 3 shows a good match between the model and the actual process. The computational time for optimization wais considerably reduced from 1 h and 29 min24 to 4 min and 17 s (maximum computational time) because of the proposed method for obtaining the initial guess values of the model parameters. The identification

Figure 6. Closed-loop responses of the corrupted process and the identified FOPTD transfer function model for example 2 (with noise). Responses obtained for step changes in (left) yr1 and (right) yr2. 1330



ARTICLE

Figure 7. Action and interaction for example 3 for step changes in the set points. In each panel, at steady state, the sum of the three responses equals 1.

method was also carried out using the closed-loop main and interaction responses obtained by different PI/PID settings. The converged model parameters show good agreement with the actual process parameters. The numbers of iterations, computational times, and IAE and ISE values are given in Table 2. The effects of initial guess values on the number of iterations, computational time, and IAE and ISE values were also studied by perturbing the guess values by +10% and separately by 10%. The present method shows a good agreement with the actual process parameters. The numbers of iterations, computational

times, and IAE and ISE values are listed in Table 3 for the perturbed guess values. The effects of measurement noise on the guess values and model parameters were also considered by adding random noise (standard deviation = 0.05) to the actual main and interaction responses and using a sample time of 0.1 s. The noisy output was used for feedback control action and also for model identification. To obtain the initial guess values, a smooth curve was first drawn to estimate the initial guess values only. Because of the integral action in the controller, the main response reached 1, and the interaction reached 0 at steady state. 1331



ARTICLE

Table 6. FOPTD Model Identification Results for Example 3 Using the PI Controller33 model ij

Kpij

11 21

guess

1.2739

converged guess converged 31

τij

θij

10.0000

2.6000

0.6598

6.6958

2.6007

4.3027

12.5000

6.5000

1.1093

3.2468

6.5012

guess

22.6872

10.0000

9.2000

converged

IAE

ISE

1.75 103

6.36 108

6.85 103

1.44 106

5.01 102

6.58 105

34.6606

8.1438

9.1991

12

guess

1.7594

12.5000

3.5000

1.06 103

2.12 108

22

converged guess

0.6096 6.4516

8.6365 10.0000

3.5021 3.0000

3.55 103

1.76 107

converged

2.2988

4.9947

3.0009

32

guess

28.3333

12.5000

9.4000

8.62 103

1.80 106

6.09 105

1.55 1010

converged

46.1734

10.8923

9.4012

guess

0.0030

10.0000

1.0000

converged

0.0049

9.0634

0.9996

23

guess

0.0064

8.7500

1.2000

1.50 104

5.95 1010

33

converged guess

0.0100 0.3279

7.0942 6.2500

1.2007 0.8000

1.11 103

2.81 108

converged

0.8500

6.5995

0.8300

13

Figure 8. Simplified diagram for 3 3 multivariable system for calculating guess values for kp31 and kp32.

Hence, a horizontal was line drawn for main response at the value of 1, and a horizontal line was drawn at the value of 0 for the interaction response in the steady-state portion only. In the initial dynamics part, a smooth curve was passed through the midpoint of the response. Alternatively, one can use the least-squares method proposed by Savitzky and Golay,28 which filters the noise and gives a smooth curve. From that smooth curve, the initial guess values were calculated by the proposed method. The identified model parameters were closer to those obtained without measurement noise. The closed-loop main responses and interaction responses using the identified model with the same controller settings are shown in the Figure 4, along with the corrupted process responses. The initial guess values, final converged values, numbers of iterations, computational times, and IAE and ISE values are also listed in Table 1. The proposed method is robust for the controller settings and for the measurement noise. 3.2. Example 2. Consider the Niederlinski31 multivariable transfer function model considered by Viswanathan et al.24 2

0:5 2 2 6 6 ð0:1s þ 1Þ ð0:2s þ 1Þ Gp ðsÞ ¼ 6 1:0 4 ð0:1s þ 1Þð0:2s þ 1Þ2

3 1:0 7 ð0:1s þ 1Þð0:2s þ 1Þ2 7 7 2:4 5 2 ð0:1s þ 1Þð0:2s þ 1Þ ð0:5s þ 1Þ

ð13Þ

Here, the higher-order TITO process transfer function model was approximated by FOPTD models. Decentralized PI controllers were used for this identification process. The decentralized PI controllers parameters were24 Kc11 = 0.4, τi11 = 0.4; Kc22 = 0.1, τi22 = 0.3. The closed-loop responses of the actual process were obtained by the previously mentioned methods using Simulink. The initial guess values were obtained by the proposed method (refer to Table 4 for the guess values). Because the response showed a higher-order model behavior (initial zero slope), kpkc,des = 0.5 was considered for calculating the guess value for the process gains. The initial guess values, final converged parameters, numbers of iterations, computational times, and IAE and ISE values are listed in Table 4. The closed-loop main and interaction responses using the identified FOPTD models with the same controller settings are shown in the Figure 5, along with the actual process responses. This figure shows a good match between the FOPTD model and the actual process. The computational time was significantly reduced to 1 min and 26 s (maximum computational time) when compared to 1 h and 20 min as reported by Viswanathan et al.24 Because the dynamics of the higher-order system of example 2 is fast (settling time is about 10 min), the computational time for the optimization problem is less than that required for example 1 (settling time is about 70 min for example 1). The effects of the initial guess values on the number of iterations, computational time, and IAE and ISE values were also studied by changing the guess values by +10% and separately by 10%. The present method showed good agreement with the actual process parameters for these perturbed guess values. The numbers of iterations, computational times, and IAE and ISE values are listed in Table 5 for these perturbed guess values. The effects of measurement noise were also evaluated on the identified FOPTD model parameters by adding a random noise (standard deviation = 0.05) to the actual main and interaction responses and considering a sample time of 0.1 s. The noisy output was used for the feedback control action and for the 1332



ARTICLE

Figure 9. Closed-loop responses of the process and the identified FOPTD transfer function model for example 3 using a PI controller.33 (Solid line, model; dashed line, process. Model and process response curves overlap.)

Table 7. Effects of Changing Controller Settings on the Number of Iterations and Computational Time for Example 3 ref (parameters: Kc11, τi11, τd11; Kc22, τi22, τd22; Kc33, τi33, τd33) Ham and Kim

19

NIa 50

IAEb 0.6222

ISEb 0.0096

Table 8. Effects of Initial Guess Values of the Model Parameters on the Number of Iterations and Computational Time for Example 3 Using the PID Controller Settings33 variation

30.53

(2.794, 4.230, 1.102; a

0.524, 5.635, 1.110; 3.835, 6.006, 1.069) Chien et al.34

42

0.056

4.08 105

33.23

(1.08, 4.25; 0.2333, 3.32; 2.78, 5.24) a NI, number of iterations. b Sum of main and interaction response values.

model identification. To obtain the initial guess values for the FOPTD models, a smooth curve was first drawn. As discussed earlier, the initial guess values for the model parameters were obtained from the smooth curve. The identified model parameters were closer to those obtained without measurement noise. The closed-loop main responses and interaction responses using the identified FOPTD model with the same controller settings are shown in the Figure 6, along with the actual corrupted process main and interaction responses. The initial guess values, final converged values, numbers of iterations, computational times, and IAE and ISE values for FOPTD model identification with noise are also listed in Table 4. 3.3. Example 3. Consider the distillation column process transfer function matrix reported by Ogunnaike and Ray32 and

computational

guess

computational time (min)

NIa

IAE

ISE

time (min) 4

10%

25

0.0939

4.32 10

20

+10%

27

0.0476

2.37 105

21.31


Ham and Kim19 2 0:66e2:6s 6 6 6:7s þ 1 6 6 1:11e6:5s G¼6 6 3:25s þ 1 6 4 34:68e9:2s 8:15s þ 1

0:61e3:5s 8:64s þ 1 2:3e3s 5:0s þ 1 46:2e9:4s 10:9s þ 1

3 0:0049e1:02 7 9:06s þ 1 7 7 0:01e1:2s 7 7 7:09s þ 1 7 7 5 0:85e0:88s 6:60s þ 1 ð14Þ

The system is interactive, as indicated by its steady-state relative gain array matrix 2 3 2:0451 0:7554 0:2896 6 7 0:1938 7 Λ¼6 ð15Þ 4 0:6730 1:8667 5 0:3721 0:1113 1:4834 Decentralized PI/PID controllers were used for the MIMO identification. The PI controller parameters33 were Kc11 = 1.57, 1333



ARTICLE

Table 9. Guess Values Obtained by the arx Method and the Final Converged Model Parameters (Obtained by lsqnonlin) for Example 1 Kpij

model

θij

IAE

ISE

time (NI)a 37 (11)

11

guess

8.9424

9.0750

0.9900

0.4470

0.0394

21

converged guess

17.0097 12.7846

19.6764 36.8000

0.4689 6.9900

0.4070

0.0617

converged

7.8409

10.8091

7.8595 0.1606

0.0042

0.2570

0.0138

12 22 a

τij

7.3675

12.0110

3.0000

converged

29.3309

29.2621

2.5819

guess

28.1234

33.7490

3.0000

converged

23.8196

19.1251

2.5529

guess


τi11 = 5.96; Kc22 = 0.31, τi22 = 4.81; Kc33 = 6.1, τi33 = 9.60. The step responses of the closed-loop system are shown in Figure 7a for a step change in the set point of y1r, in Figure 7b for a step change in y2r, and in Figure 7c for a step change in y3r. The interactions are significant. The FOPTD model was assumed, and the initial guess values for the model parameters were obtained by using step responses (refer to Table 6 for the guess values). The guess values for the main-loop (diagonal) process gains kp11, kp22, and kp33 were calculated based on the corresponding controller gains as discussed in example 1. Because the responses of the actual process showed a lower-order model behavior (initial zero slope), the kpkc,des value was considered as 2 for calculating the guess value of the process gains (g11, g22, and g33). The nondiagonal process gains (interactions) were calculated using equations similar to eqs 8 and 10. To simplify the calculations, only two loops were considered at a time, and the interaction coming from the remaining loop was considered as negligible. For example, the guess value of the process gain (nondiagonal) kp21 was estimated by considering the interaction coming from main loop only (similar to Figure 2a). The interaction to and from the third loop were assumed to be neglected. This procedure gave explicitly the guess value of the required process gain. Similarly, to obtain the guess values of kp12, an approach similar to Figure 2b was considered. To obtain the guess values for the gains kp31 and kp32, the corresponding Figure 8a,b was considered. Similarly, the guess values for the other gains kp13 and kp23 were calculated. The signs of the guess values of the steady-state gains are important. With these guess values, the convergence of the parameters of the model was achieved, as shown in Table 6. The responses obtained from g13 and g23 are of very low order of magnitude compared with other responses. Hence, the responses obtained from all of the models were normalized between 0 and 1. These normalized data were used for the optimization routine. For simulating the model, the sample time was considered as 0.01 s. The initial guess values, final converged parameters, numbers of iterations, and IAE and ISE values are presented in Table 6. The converged model parameters match well with the actual FOPTD model parameters. The computational time for the optimization was 23 min and 20 s, and the number of iterations was 29. The closed-loop main and interaction responses using the identified model parameters with the same controller settings are also shown in Figure 9. Figure 9 shows a good match between the model and the actual process. For example 3, because the number of parameters to be identified was larger (27 parameters) and the dynamics of the closed-loop system requires a longer settling time (refer to Figure 9), the

computational time for the optimization problem was higher (about 23 min) than those for examples 1 and 2. The identification method was also carried out using the closed-loop main and interaction responses obtained with different PI/PID settings. The converged model parameters showed good agreement with the actual process parameters. The numbers of iterations, computational times, and IAE and ISE values are listed in Table 7.The effects of initial guess values on the number of iteration, computational time, and IAE and ISE values were also studied by perturbing the guess values by +10% and separately by 10%. The model parameters obtained by the proposed method showed good agreement with the actual process parameters for these perturbed guess values. The numbers of iterations, computational times, and IAE and ISE values are listed in Table 8. The initial guess values for the model parameters can also be obtained by assuming the discrete time ARX model35 and solving the linear regression equations to get the parameter estimations. The guess values for time delays were identified from the impulse response plot (by using the Matlab command cra). After the estimate for the time delays had been obtained, the standard routine arx (with na = 3, nb = [2,2]) in Matlab was used to get the arx model. For solving the linear regression equations, the MIMO model was first converted into a MISO (multi-input and single-output) model. The identified higher-order ARX model was reduced to an FOPTD model. The guess values obtained by the above method were used in the nonlinear optimization (lsqnonlin) to get the optimal model parameters. For example 1, the initial guess values, final model parameters, numbers of iterations, computational times, and IAE and ISE values are listed in Table 9. The guess values obtained by the arx method were of the same order of magnitude as the true parameters. However, with these guess values, satisfactory results could not be obtained. The method of getting the initial guess values discussed in section 2 is simpler to use and also gives satisfactory results. The proposed work considers only the identification of FOPTD model for 2 2 and 3 3 systems. The identification of transfer function models with zero and underdamped SOPTD transfer function models of multivariable systems needs further study.

4. CONCLUSIONS Closed-loop identification of FOPTD models of multivariable systems based on an optimization method using step response was presented. Simple methods for obtaining initial guess values 1334


Industrial & Engineering Chemistry Research for the model parameters from the main and interaction responses were proposed for getting quick and guaranteed convergence. By this proposed method, the closed-loop main and interaction responses of the identified model showed good matching with the actual process. The identified model parameters were closer to the actual process parameters. The effects of measurement noise and controller settings on the identification method were also studied. The proposed method was found to be robust to measurement noise and controller settings. The proposed method was applied to 2 2 FOPTD, 2 2 higher-order, and 3 3 FOPTD transfer function matrix models. The computational time for the optimization method is considerably reduced from 1 h and 29 to 4 min and 17 s for 2 2 FOPTD model systems and 1 h and 20 min15 to 1 min and 26 s for 2 2 higher-order systems. The computational times for 3 3 FOPTD model system wer 33 min and 14 s.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ NOMENCLATURE G = process transfer function matrix g = process transfer function GC = controller transfer function matrix gc = controller transfer function gm = model transfer function IAE = integral of absolute error ISE = integral of squared error Kc = controller gain Kmp = model steady-state gain Kp = process steady-state gain t = time tc = closed-loop time constant of the model ts = settling time Y = closed-loop process response matrix y = closed-loop process response ym = closed-loop model response ymr = model set point yr = process set point θ = process time delay θm = model time delay τ = process time constant τd = derivative time τi = integral time τm = model time constant ’ REFERENCES (1) Skogestad, S.; Postlethwaite, I. Multivariable Feedback Control: Analysis and Design, 2nd ed.; John Wiley & Sons: New York, 2005. (2) Gustavsson, I.; Ljung, L.; Soderstrom, T. Identification of Processes in Closed Loop—Identifiability and Accuracy Aspects. Automatica 1977, 13, 59–75. (3) Landau, I. D. Identification in Closed Loop: A Powerful Design Tool (Better Design Models, Simpler Controllers). Control Eng. Pract. 2001, 9, 51–65. (4) Hjalmarsson, H.; Gevers, M.; Bruyne, F. D. For Model-Based Control Design, Closed-Loop Identification Gives Better Performance. Automatica 1996, 32 (12), 1659–1673. (5) Langer, J.; Landau, I. D. Improvement of Robust Digital Control by Identification in the Closed Loop. Application to a 360° Flexible Arm. Control Eng. Pract. 1996, 4 (12), 1637–1646.

ARTICLE

(6) Sung, S. W.; Lee, I. B.; Lee, B. K. On-line Process Identification and Automatic Tuning Method for PID Controllers. Chem. Eng. Sci. 1998, 53 (No 10), 1847–1859. (7) Sung, S. W.; Lee, I. B. On-line Process Identification and PID Controller Auto Tuning. Korean J. Chem. Eng. 1999, 16 (No 1), 45–55. (8) Sung, S. W.; Lee, I. B. Prediction Error Identification Method for Continuous-Time Processes with Time Delay. Ind. Eng. Chem. Res. 2001, 40, 5743–5751. (9) Sung, S. W.; Lee, J.; Lee, I. B. Process Identification and PID Control; John Wiley & Sons (Asia) Pvt Ltd: Singapore, 2009. (10) Viswanathan, P. K.; Rangaiah, G. P. Process Identification from Closed-Loop Response Using Optimization Methods. Trans. Inst. Chem. Eng. A 2000, 78, 528–541. (11) Pramod, S.; Chidambaram, M. Closed Loop Identification of Transfer Function Model for Unstable Bioreactors for Tuning PID Controllers. Bioprocess Biosyst. Eng. 2000, 22, 185–188. (12) Pramod, S.; Chidambaram, M. Closed Loop Identification by Optimization Method for Tuning PID Controllers. Indian Chem. Eng. A 2001, 43 (No.2), 90–94. (13) Sree, R. P.; Chidambaram, M. Identification of Unstable Transfer Model with a Zero by Optimization Method, J. Indian Inst. Sci. 2002, 82, 219–225. (14) Ramachandran, R.; Lakshimanarayanan, S.; Rangaiah, G. P. Process Identification Using Open-Loop and Closed-Loop Step Responses. J. Inst. Eng. Singapore 2005, 45 (6), 1–13. (15) Mclntosh, A. R.; Mahalec, V. Calculation of Steady-State Gains for Multivariable Systems from Closed-Loop Steady-State Data. J. Process Control 1991, 1, 178–186. (16) Rajakumar, A.; Krishnaswamy, P. R. Time to Frequency Domain Conversion of Step Response Data. Ind. Eng. Chem. Process Des. Dev. 1975, 14 (3), 250–256. (17) Melo, D. L.; Friedly, J. C. On-Line, Closed-Loop Identification of Multivariable Systems. Ind. Eng. Chem. Res. 1992, 31, 274–281. (18) Semino, D.; Scali, C. Improved Identification and Autotuning of PI Controllers for MIMO Processes by Relay Techniques. J. Process Control 1998, 8 (3), 219–227. (19) Ham, T. W.; Kim, Y. H. Process Identification and PID Controller Tuning in Multivariable Systems. J. Chem. Eng. Jpn. 1998, 31 (No.6), 941–949. (20) 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. (21) Mei, H.; Li, S. Decentralized Identification for Multivariable Integrating Processes with Time Delays from Closed-Loop Step Tests. ISA Trans. 2007, 46, 189–198. (22) Li, S. Y.; Cai, W. J.; Mei, H.; Xiong, Q. Robust Decentralized Parameter Identification for Two-Input Two-Output Process from Closed-Loop Step Responses. Control Eng. Pract. 2005, 13, 519–531. (23) Yagang, W. Closed-Loop Multivariable Process Identification in the Frequency Domain. In Proceedings of the 27th Chinese Control Conference; Cheng, D. Z. and Li, C., Eds.; IEEE Press: Piscataway, NJ, 2008; pp 440445. (24) Viswanathan, P. K.; Toh, W. K.; Rangaiah, G. P. Closed-Loop Identification of TITO Processes Using Time-Domain Curve Fitting and Genetic Algorithms. Ind. Eng. Chem. Res. 2001, 40, 2818–2826. (25) Papastathopoulou, H. S.; Luyben, W. L. A New Method for the Derivation of Steady State Gains for Multivariable Process. Ind. Eng. Chem. Res. 1990, 29, 366–369. (26) Chidambaram, M. Applied Process Control; Allied Publishers: New Delhi, India, 1998; p 7. (27) Wood, R. K.; Berry, M. W. Terminal Composition Control of a Binary Distillation Column. Chem. Eng. Sci. 1973, 28, 1707–1717. (28) Savitzky, A.; Golay, M. J. E Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Anal. Chem. 1964, 36 (8), 1627–1639. (29) Vu, T. N. L.; Hong, S.; Lee, M. Analytical Design of Robust Multi-Loop PI Controller for Multivariable Processes. In ICROS-SICE 1335



ARTICLE

International Joint Conference; IEEE Press: Piscataway, NJ, 2009; pp 29612966. (30) Loh, P. A.; Hang, C. C.; Quek, K. C.; Vasnani, U. V. Autotuning of Multiloop ProportionalIntegral Controllers Using Relay Feedback. Ind. Eng. Chem. Res. 1993, 32, 1102–1107. (31) Niederlinski, A. A Heuristic Approach to the Design of Linear Multivariable Interacting Control Systems. Automatica 1971, 7, 691–701. (32) Ogunnaike, B. A.; Lemaire, J. P.; Morari, M.; Ray, W. H. Advanced Multivariable Control of a Pilot Plant Distillation Column. AIChE J. 1983, 29 (4), 632–640. (33) Vu, T. N. L.; Lee, M. Multi-Loop PI Controller Design Based on the Direct Synthesis for Interacting Multi-Time Delay Processes. ISA Trans. 2010, 49, 79–86. (34) Chien, I. L.; Huang, H. P.; Yang, J. C. A Simple Multiloop Tuning method for PID Controllers with No Proportional Kick. Ind. Eng. Chem. Res. 1999, 38, 1456–1468. (35) Lyzell, C. Initialization Methods for System Identification. Licentiate Thesis no. 1426, Link€oping University, Link€oping, Sweden, 2009.

1336


Closed-Loop Identification of Multivariable Systems by Optimization

Recommend Documents