**1. Introduction**

Batch and continuous systems are of multivariable in nature. A multivariable system is one in which one input not only affects its own outputs but also one or more other outputs in the plant. Multivariable processes are difficult to control due to the presence of the interactions. Increase in complexity and interactions between inputs and outputs yield degraded process behavior. Such processes are found in process industries as they arise from the design of plants that are subject to rigid product quality specifications, are more energy efficient, have more material integration, and have better environmental performance. Most of the unit operations in process industry require control over product rate and quality by adjusting one/more inputs to the process; thus making multivariable systems. For example, chemical reactors, distillation column, heat exchanger, fermenters are typical multivariable processes in industry. In case of chemical reactor, the output variables are product composition and temperature of reaction mass. The input variables are reactant or feed flow rate and energy added to the system by heating/ cooling through jackets. Product composition can be controlled by manipulating feed rate whereas rate of reaction (thereby temperature) can be controlled by changing addition/ removal rate of energy. But, while controlling product composition, temperature is affected; similarly, while controlling temperature of reaction mass, the composition gets affected, thus, exhibiting interactions between input and output variables. Distillation is widely used for separating components from mixture in refineries. Composition of top and bottom products are controlled by adjusting energy input to the column. A common scheme is to use reflux flow to control top product composition whilst heat input is used to control bottom product composition. However, changes in reflux also affect bottom product composition and component fractions in the top product stream are also affected by changes in heat input. Hence, loop interactions occur in composition control of distillation column. Thus, unless proper precautions are taken in terms of control system design, loop interactions can cause performance degradation and instability. Control system design needs availability of linear models for the multivariable system.

The basic and minimum process model for multivariable system is considered here as 2x2 system. The outputs of the loops are given by

<sup>\*</sup> Corresponding Author

$$\begin{aligned} y\_1 &= \mu\_1 G\_{P11} + \mu\_2 G\_{P12} \\ y\_2 &= \mu\_2 G\_{P21} + \mu\_2 G\_{P22} \end{aligned} \tag{1}$$

Identification and Control of Multivariable Systems – Role of Relay Feedback 107

Moreover, with weak interactions and with large dimensional systems they induce to go for more criteria for selection of pairs. Morari resiliency index (MRI) is also used to select in/out

of the plant transfer function matrix G(iw). The set of manipulated variables that gives the largest minimum singular value over the frequency range of interest is the best. The MRI is a measure of the inherent ability of the process (control structure) to handle disturbances, model plant mismatches, changes in operating conditions, etc. The larger the value of MRI, the more resilient the control structure. Dynamic Relative Gain Array (DRGA) is defined to extend the RGA notion to non-zero frequencies. The RGA provides only limited knowledge about when to use multivariable controllers and gives no indication of how to choose multivariable controller structures. A somewhat different approach for investigating channel interaction was therefore employed by Conley and Salgado (2000) and Salgado and Conley (2004) when considering observability and controllability gramians in so called Participation Matrices (PM). In a similar approach Wittenmark and Salgado (2002) introduced the Hankel Interaction Index Array (HIIA). These gramian based interaction measures seem to overcome most of the disadvantages of the RGA. One key property of these is that the whole frequency range is taken into account in one single measure. Furthermore, these measures seem to give appropriate suggestions for controller structures selection. The use of the system H2 norm as a base for an interaction measure has been proposed by Birk and Medvedev (2003) as an alternative to the HIIA. But, dynamic simulation is a powerful tool to be used to test the viability of a control scheme during various process disturbances. Controllers for MIMO systems can be of either multiloop (controllers are designed only for diagonal elements of process TF) or multivariable (controllers are designed for all the elements of the MIMO TF). Multiloop control scheme has an edge over multivariable as the former can work even if a single loop fails. In presence of interactions between input/output, the process need to be decoupled and then multiloop controllers can be designed. When interaction effects produce a significant deterioration in control system performance, decoupling control should be considered. One of the most powerful and simplest ways of reducing or eliminating interaction is by altering manipulated and / or controlled variables. Improvement of closedloop performance needs proper tuning of controller parameters that requires process model structure and estimation of respective parameters. There are many methods to select input/output pairs or to design control structures, design control strategy (either PID or IMC or predictive or heuristics etc.) and tuning of controller parameters in literature. But because of hazy pictures on above selections, till today, it is difficult to choose correct pairs, carryout interaction analysis and choose tuning rules. Thus the aim of this chapter is to bring out a clear picture of identifying process parameters and designing controller for MIMO systems. The rest of the chapter is carried out as follows: section 2 discusses identification methods of multivariable systems. Interaction analysis is explained in section 3. Control structure selection and determination of input/output pairs are given in section 4. Tuning of controllers is presented in section 5. Stability analysis for multivariable systems is provided in section 6. At

Most of the chemical and bio-chemical processes are multivariable in nature, having more than one input and outputs. Estimation of process parameters is a key element in

where is eigenvalue. The MRI is the minimum singular value *(g)* 

pairs. *MRI G*

 *<sup>P</sup>*( ) *j*

the end, conclusion is drawn.

**2. System identification** 

where yi are system outputs and ui are the system inputs, G is system transfer functions. Eqn (1) can be expressed as *<sup>P</sup> y G u* where

$$y = \begin{bmatrix} y\_1 & y\_2 \end{bmatrix}^T; \quad u = \begin{bmatrix} u\_1 & u\_2 \end{bmatrix}^T; \quad G\_P = \begin{bmatrix} G\_{P11} & G\_{P12} \\ G\_{P21} & G\_{P22} \end{bmatrix} \tag{2}$$

In order to achieve desired quality, specified output characteristics at the cost of spending optimum inputs one needs to design a controller and run the plant under closed loop so that optimal production of product under safe operation. The first thing we need is to select input-output pairs, i.e., which output should be controlled by which input? This needs knowledge in control structure selection or interaction analysis. In the next section, a brief state of art on interaction analysis is presented.

Relative gain array (RGA) (Bristol 1966) is the most discussed method for analyzing interactions and it is based on steady state gain information of MIMO processes. Control loops should have input-output pairs which give positive relative gains that have values which are as close as unity as possible. It is dependent on process models, independent of scaling of inputs and outputs and can include all ways of pairing in a single matrix. Niederlinski index (NI) is a useful tool to analyse interactions and stability of the control loop pairings determined using process gain matrix. NI is found by the following formula,

$$\text{NI} = \frac{|\mathbf{G}\_P|}{\prod\_{i=1}^n g\_{ii}} \quad |\_{\text{SS}} \text{ where each element of } \mathbf{G}\_P \text{ is rational and is open} \\ \text{top stable. The values of } \prod\_{i=1}^n g\_{ii}$$

NI need to be positive. A negative value of NI will imply that the system is un-stable. Ni is used to check if the system (more than 2x2) is unstable or not. NI will detect instability introduced by closing the other control loops. Generally, NI is not used for systems with time delays. Any loop pairing is unacceptable if it leads to a control system conguration for which the NI is negative. But both RGA & NI do not provide dynamic information on the process transients. They do not give information on change in in/op pairing for instances when there is a sudden load disturbance. Singular value decomposition (SVD) is a useful tool to determine whether a system will be prone to control loop interactions resulting in sensitivity problems that rises from model mismatch in process gains. SVD considers directional changes in the disturbances. SVD is applied to steady state gain matrix that is decomposed into product of three matrices,

*<sup>T</sup> S UV* where U is matrix of normalized eigen vectors of *<sup>T</sup> GGP* , is diagonal matrix of eigenvalues and V is matrix of normalized eigenvectors of *<sup>T</sup> G GP P* The condition number (CN) is defined as ratio between maximum and minimum eigenvalues. Generally if the CN < 50 then the system is not prone to sensitivity problems (a small error in process gain will not cause a large error in the controller's reactions). The greater the CN value, the harder it is for the system in question to be decoupled. An ideal system would have a CN number of one, where each control variable controls a single distinct output variable. CN value tells us how easy it is to decouple a system. Though SVD has good geometric interpretation in terms of selection of measurement and pairing of variables, SVD depends on input-output scaling.

1 1 11 2 12 2 2 21 2 22 *P P P P*

where yi are system outputs and ui are the system inputs, G is system transfer functions.

 [ ] ; [ ] ; = *T T P P*

In order to achieve desired quality, specified output characteristics at the cost of spending optimum inputs one needs to design a controller and run the plant under closed loop so that optimal production of product under safe operation. The first thing we need is to select input-output pairs, i.e., which output should be controlled by which input? This needs knowledge in control structure selection or interaction analysis. In the next section, a brief

Relative gain array (RGA) (Bristol 1966) is the most discussed method for analyzing interactions and it is based on steady state gain information of MIMO processes. Control loops should have input-output pairs which give positive relative gains that have values which are as close as unity as possible. It is dependent on process models, independent of scaling of inputs and outputs and can include all ways of pairing in a single matrix. Niederlinski index (NI) is a useful tool to analyse interactions and stability of the control loop pairings determined using process gain matrix. NI is found by the following formula,

NI need to be positive. A negative value of NI will imply that the system is un-stable. Ni is used to check if the system (more than 2x2) is unstable or not. NI will detect instability introduced by closing the other control loops. Generally, NI is not used for systems with time delays. Any loop pairing is unacceptable if it leads to a control system conguration for which the NI is negative. But both RGA & NI do not provide dynamic information on the process transients. They do not give information on change in in/op pairing for instances when there is a sudden load disturbance. Singular value decomposition (SVD) is a useful tool to determine whether a system will be prone to control loop interactions resulting in sensitivity problems that rises from model mismatch in process gains. SVD considers directional changes in the disturbances. SVD is applied to steady state gain matrix that is

*y y y u uu G G G* 

*P*

where each element of GP is rational and is openloop stable. The values of

where U is matrix of normalized eigen vectors of *<sup>T</sup> GGP* , is diagonal matrix of

eigenvalues and V is matrix of normalized eigenvectors of *<sup>T</sup> G GP P* The condition number (CN) is defined as ratio between maximum and minimum eigenvalues. Generally if the CN < 50 then the system is not prone to sensitivity problems (a small error in process gain will not cause a large error in the controller's reactions). The greater the CN value, the harder it is for the system in question to be decoupled. An ideal system would have a CN number of one, where each control variable controls a single distinct output variable. CN value tells us how easy it is to decouple a system. Though SVD has good geometric interpretation in terms of selection of measurement and pairing of variables, SVD depends on input-output scaling.

(1)

21 22

(2)

*P P G G*

*y uG uG y uG uG* 

 11 12 1 2 1 2

Eqn (1) can be expressed as *<sup>P</sup> y G u* where

state of art on interaction analysis is presented.

decomposed into product of three matrices,

1

*g* 

*G*

*NI*

 

*<sup>T</sup> S UV* 

 | *<sup>P</sup> n SS ii i*

Moreover, with weak interactions and with large dimensional systems they induce to go for more criteria for selection of pairs. Morari resiliency index (MRI) is also used to select in/out pairs. *MRI G <sup>P</sup>*( ) *j* where is eigenvalue. The MRI is the minimum singular value *(g)*  of the plant transfer function matrix G(iw). The set of manipulated variables that gives the largest minimum singular value over the frequency range of interest is the best. The MRI is a measure of the inherent ability of the process (control structure) to handle disturbances, model plant mismatches, changes in operating conditions, etc. The larger the value of MRI, the more resilient the control structure. Dynamic Relative Gain Array (DRGA) is defined to extend the RGA notion to non-zero frequencies. The RGA provides only limited knowledge about when to use multivariable controllers and gives no indication of how to choose multivariable controller structures. A somewhat different approach for investigating channel interaction was therefore employed by Conley and Salgado (2000) and Salgado and Conley (2004) when considering observability and controllability gramians in so called Participation Matrices (PM). In a similar approach Wittenmark and Salgado (2002) introduced the Hankel Interaction Index Array (HIIA). These gramian based interaction measures seem to overcome most of the disadvantages of the RGA. One key property of these is that the whole frequency range is taken into account in one single measure. Furthermore, these measures seem to give appropriate suggestions for controller structures selection. The use of the system H2 norm as a base for an interaction measure has been proposed by Birk and Medvedev (2003) as an alternative to the HIIA. But, dynamic simulation is a powerful tool to be used to test the viability of a control scheme during various process disturbances. Controllers for MIMO systems can be of either multiloop (controllers are designed only for diagonal elements of process TF) or multivariable (controllers are designed for all the elements of the MIMO TF). Multiloop control scheme has an edge over multivariable as the former can work even if a single loop fails. In presence of interactions between input/output, the process need to be decoupled and then multiloop controllers can be designed. When interaction effects produce a significant deterioration in control system performance, decoupling control should be considered. One of the most powerful and simplest ways of reducing or eliminating interaction is by altering manipulated and / or controlled variables. Improvement of closedloop performance needs proper tuning of controller parameters that requires process model structure and estimation of respective parameters. There are many methods to select input/output pairs or to design control structures, design control strategy (either PID or IMC or predictive or heuristics etc.) and tuning of controller parameters in literature. But because of hazy pictures on above selections, till today, it is difficult to choose correct pairs, carryout interaction analysis and choose tuning rules. Thus the aim of this chapter is to bring out a clear picture of identifying process parameters and designing controller for MIMO systems. The rest of the chapter is carried out as follows: section 2 discusses identification methods of multivariable systems. Interaction analysis is explained in section 3. Control structure selection and determination of input/output pairs are given in section 4. Tuning of controllers is presented in section 5. Stability analysis for multivariable systems is provided in section 6. At the end, conclusion is drawn.

## **2. System identification**

Most of the chemical and bio-chemical processes are multivariable in nature, having more than one input and outputs. Estimation of process parameters is a key element in

Identification and Control of Multivariable Systems – Role of Relay Feedback 109

The discrete transfer function has three parameters that need to be identified: dead time (D)

where *<sup>n</sup> y* is the predicted value of the current output of the process. For a FOPDT process,

*y* 

 

0 0 1 1 ................ ..........

*y u y u*

*yN uN*

model as time goes on with a reduced computational cost.

1 1

<sup>1</sup> *T T*

In the state space form the relationship between the input, noise and output signals are written as a system of first-order differential or difference equations using auxillary state vectors. Transfer function in laplace domain is converted to state space form using a

The beginning of the 1990s witnesses the birth of a new type of linear system identification algorithms, called subspace method. Subspace identification methods are indeed attractive since a state-space realization can be directly estimated from input/output data without nonlinear optimization. Furthermore, these techniques are characterized by the use of robust numerical tools such as RQ factorization and the singular values decomposition (SVD). Interesting from numerical point of view, the batch subspace model identification (SMI) algorithms are not usable for online implementation because of the SVD computational complexity. Indeed, in many online identification scenarios, it is important to update the

Linear subspace identification methods are concerned with systems and models of the form

is state matrix and y is outputs.

 *<sup>y</sup>*

<sup>1122</sup> <sup>1122</sup> .... .... *nn n nb n nb n n na n na y bu bu b u ay ay a y* (2.3)

*y k ay k bu k* 1 1 1 1 (2.4)

and <sup>1</sup>

1 *a b*

(2.6)

*k k kk* <sup>1</sup> *x Ax Bu w* (2.7)

*e* (2.5)

) contained in b1, and a1.

contained in nk, and other two parameters of the model (kp and

The discrete output can be represented in the following form:

equation (2.3) can be written as

where

where 

**2.2 State-space model** 

sampling period of 0.1s

**2.2.1 Subspace method** 

which can be written in matrix form as

The parameters a1 and b1 are calculated using

is the parameter vector

multivariable controller design. Thus, as better performance is achieved by model based tuning algorithms, estimation of model structures are necessary from either open-loop or closed-loop data. This is due to the fact that tuning rules are based on model structures & parameters. As their exist advantages and disadvantages in both of these identification strategies, for example, open-loop responses may show unstable behavior with certain inputs, whereas, closed-loop strategy needs more excitation to yield observable response. Here we use mostly used methods of identification for multivariable systems. Least square method (Tungnait 1998) is an old but reliable technique that was in use to estimate multivariable parameters of open-loop systems. But, MIMO systems with interactions may not yield satisfactory transfer function estimates with these techniques. Overschee and Moor (1994) proposed subspace method of identification that mostly applies to identification of multivariable state space models. This method involves more computational time. Practical industrial plants are easy to identify in closed-loop using relay feedback method (Astrom and Hagguland 1984) and Yu (1999) explains advances in autotuning using sequential identification. System identification is the method of estimating parameters from system's input/output data using numerical techniques:

#### **2.1 Transfer function identification**

Model structures and parameters of transfer function are constructed from observed plant input output data. Transfer function models are developed using three schemes: (a) Least square (b) subspace and (c) sequential identification method. These approximations made out through each of the methods carry errors that propagate to controller tuning and in turn deteriorates the overall performance.

#### **2.1.1 Least-squares method**

Least-squares method, used to reduce the mean square error, is very simple and more numerically stable and can be used to identify the unknown parameters of the 2x2 MIMO transfer function model from the input (u) and output (y) data. Though any type of forcing function (step, pulses or a sequence of positive and negative pulses) can be used, a very popular sequence of inputs, "Pseudo-random binary sequence" (PRBS) is made use of in the present work.

Let us consider a process with continuous transfer function

$$\frac{y(s)}{u(s)} = G(s) = \frac{K\_p e^{-Ds}}{\tau s + 1} \tag{2.1}$$

The pulse transfer function of this process with a zero-order hold is

$$\frac{\partial \ln(z)}{\partial \ln(z)} = \text{HG}(z) = \frac{K\_p \left(1 - b\right)}{z^{nk} \left(z - b\right)} = \frac{K\_p \left(1 - b\right) z^{-1}}{z^{nk} \left(1 - b z^{-1}\right)} = \frac{z^{-nk} \left(b\_0 + b\_1 z^{-1}\right)}{1 + a\_1 z^{-1}}\tag{2.2}$$

where *<sup>k</sup> <sup>s</sup> n D <sup>T</sup>* ; Ts=sampling period;

$$b = e^{-\frac{T}{\lambda}\sqrt{\varepsilon}}; b\_o = 0; b\_1 = k\_p \left(1 - b\right); a\_1 = -b\_1$$

The discrete transfer function has three parameters that need to be identified: dead time (D) contained in nk, and other two parameters of the model (kp and ) contained in b1, and a1. The discrete output can be represented in the following form:

$$\overline{y\_n} = b\_1 \mu\_{n-1} + b\_2 \mu\_{n-2} + \dots + b\_{nb} \mu\_{n-nb} - a\_1 y\_{n-1} - a\_2 y\_{n-2} - \dots - a\_{na} y\_{n-na} \tag{2.3}$$

where *<sup>n</sup> y* is the predicted value of the current output of the process. For a FOPDT process, equation (2.3) can be written as

$$y(k) + a\_1 y(k-1) = b\_1 \mu(k-1) \tag{2.4}$$

which can be written in matrix form as

$$y = \phi \theta + e\tag{2.5}$$

where

108 Introduction to PID Controllers – Theory, Tuning and Application to Frontier Areas

multivariable controller design. Thus, as better performance is achieved by model based tuning algorithms, estimation of model structures are necessary from either open-loop or closed-loop data. This is due to the fact that tuning rules are based on model structures & parameters. As their exist advantages and disadvantages in both of these identification strategies, for example, open-loop responses may show unstable behavior with certain inputs, whereas, closed-loop strategy needs more excitation to yield observable response. Here we use mostly used methods of identification for multivariable systems. Least square method (Tungnait 1998) is an old but reliable technique that was in use to estimate multivariable parameters of open-loop systems. But, MIMO systems with interactions may not yield satisfactory transfer function estimates with these techniques. Overschee and Moor (1994) proposed subspace method of identification that mostly applies to identification of multivariable state space models. This method involves more computational time. Practical industrial plants are easy to identify in closed-loop using relay feedback method (Astrom and Hagguland 1984) and Yu (1999) explains advances in autotuning using sequential identification. System identification is the method of estimating parameters from system's

Model structures and parameters of transfer function are constructed from observed plant input output data. Transfer function models are developed using three schemes: (a) Least square (b) subspace and (c) sequential identification method. These approximations made out through each of the methods carry errors that propagate to controller tuning and in turn

Least-squares method, used to reduce the mean square error, is very simple and more numerically stable and can be used to identify the unknown parameters of the 2x2 MIMO transfer function model from the input (u) and output (y) data. Though any type of forcing function (step, pulses or a sequence of positive and negative pulses) can be used, a very popular sequence of inputs, "Pseudo-random binary sequence" (PRBS) is made use of in the

<sup>1</sup>

 

; 0; 1 ; 1 1

*o p b e b b k ba b*

*y s K ep G s us s*

1 1

*u z z zb z bz a z*

*y z K b K bz z b bz*

*p p nk nk* *Ds*

1 1

(2.2)

(2.1)

 <sup>1</sup> <sup>1</sup> 0 1

1

1 1

*nk*

 

input/output data using numerical techniques:

Let us consider a process with continuous transfer function

The pulse transfer function of this process with a zero-order hold is

*Ts*

*HG z*

*<sup>T</sup>* ; Ts=sampling period;

**2.1 Transfer function identification** 

deteriorates the overall performance.

**2.1.1 Least-squares method** 

present work.

where *<sup>k</sup> <sup>s</sup> n D*

$$\phi = \begin{bmatrix} -y(0) & \mu(0) \\ -y(1) & \mu(1) \\ \cdots & \cdots & \cdots \\ -y(N-1) & \mu(N-1) \end{bmatrix} \quad \text{and} \quad \theta = \begin{bmatrix} a\_1 \\ b\_1 \end{bmatrix}$$

The parameters a1 and b1 are calculated using

$$\boldsymbol{\Theta} = \left(\boldsymbol{\phi}^T \boldsymbol{\phi}\right)^{-1} \boldsymbol{\phi}^T \boldsymbol{y} \tag{2.6}$$

where is the parameter vector is state matrix and y is outputs.

#### **2.2 State-space model**

In the state space form the relationship between the input, noise and output signals are written as a system of first-order differential or difference equations using auxillary state vectors. Transfer function in laplace domain is converted to state space form using a sampling period of 0.1s

#### **2.2.1 Subspace method**

The beginning of the 1990s witnesses the birth of a new type of linear system identification algorithms, called subspace method. Subspace identification methods are indeed attractive since a state-space realization can be directly estimated from input/output data without nonlinear optimization. Furthermore, these techniques are characterized by the use of robust numerical tools such as RQ factorization and the singular values decomposition (SVD). Interesting from numerical point of view, the batch subspace model identification (SMI) algorithms are not usable for online implementation because of the SVD computational complexity. Indeed, in many online identification scenarios, it is important to update the model as time goes on with a reduced computational cost.

Linear subspace identification methods are concerned with systems and models of the form

$$\mathbf{x}\_{k+1} = A\mathbf{x}\_k + B\mathbf{u}\_k + w\_k \tag{2.7}$$

$$y\_k = \mathbb{C}\mathfrak{x}\_k + Du\_k + v\_k \tag{2.8}$$

Identification and Control of Multivariable Systems – Role of Relay Feedback 111

All the above identification methods involve more computations and many offline methods. These difficulties can be avoided easily by using another method of estimation technique,

Based on the concept of sequential auto tuning (Shen & Yu, 1994) method each controller is designed in sequence. Let's consider a 2-by-2 MIMO system with a known pairing *y*1 1 *u* and *y*2 2 *u* under decentralized PI control (Figure 1). Initially, an ideal / biased relay is placed between 1 *y* and *u*<sup>1</sup> , while loop 2 is on manual (Figure 2a). Following the relayfeedback test, a controller can be designed from the ultimate gain and ultimate frequency. The next step is to perform relay-feedback test between 2 *y* and *u*2 while loop 1 is on automatic (Figure 2b). A controller can also be designed for loop 2 following the relayfeedback test. Once the controller on the loop 2 is put on automatic, another relay-feedback experiment is performed between 1 *y* and *u*<sup>1</sup> , (Figure 2c). Generally, a new set of tuning constants is found for the controller in loop 1. This procedure is repeated until the controller parameters converge. Typically, the controller parameters converge in 3 - 4 relay-feedback

**G11**

**G21**

**G12**

**G22**

**u1** 

**Y1** 

**Y2** 

**Y1**

**Y2** 

**G21**

**G12**

**G22**

**u2** 

**Relay**

**Controll**

**u2** 

**R1 G11**

**Relay**

**u1**

namely, relay feedback method as explained below:

**-**

**2.3 Sequential identification** 

tests for 2 x 2 systems.

**R1** 

**(a)**

**(b)**

**R2** 

with

$$E\begin{bmatrix} \begin{pmatrix} w\_p \\ v\_p \end{pmatrix} \begin{pmatrix} w\_q \end{pmatrix}^T & v\_q^{\top} \end{pmatrix} \begin{bmatrix} \begin{pmatrix} Q & S \\ S^T & R \end{pmatrix} \delta\_{pq} \end{bmatrix} \tag{2.9}$$

The vectors *mx*<sup>1</sup> *u R <sup>k</sup>* and *lx*<sup>1</sup> *<sup>k</sup> y R* are the measurements at time instant k of, respectively, the m inputs and l outputs of the process. The vector xk is the state vector of the process at discrete time instant k, *lx*<sup>1</sup> *<sup>k</sup> v R* and *nx*<sup>1</sup> *w R <sup>k</sup>* are unobserved vector signals, vk is called the measurement noise and wk is called the process noise. It is assumed that they are zero mean, stationary white noise vector sequences and uncorrelated with the inputs uk. *nxn A R* is the system matrix, *nxm B R* is the input matrix, *lxn C R* is the output matrix while *lxm D R* is the direct feed-through matrix. The matrices *nxn Q R* , *nxl S R* and *lxl R R* are the covariance matrices of the noise sequences wk and vk.

In subspace identification it is typically assumed that the number of available data points goes to infinity, and that the data is ergodic. The main problem of identification is arranged as follows:

Given a large number of measurements of the input uk and the output yk generated by the unknown system described by equations (2.7)-(2.9). The task is to determine the order n of the unknown system, the system matrices A, B, C, D up to within a similarity transformation and an estimate of the matrices Q, S and R.

Subspace identification algorithms always consist of two steps:

Step 1: Make a weighted projection of certain subspace generated from the data, to find an

estimate of the extended observability matrix, *<sup>i</sup>* and/or an estimate *Xi* of the state sequence *Xi* of the unknown system

Step 2: Retrieve the system matrices (A, B, C, D and Q, S, R) and from either this extended observability matrix ( *<sup>i</sup>* ) or the estimated states.

Fig. 1. Flow chart of subspace algorithm.

All the above identification methods involve more computations and many offline methods. These difficulties can be avoided easily by using another method of estimation technique, namely, relay feedback method as explained below:

### **2.3 Sequential identification**

110 Introduction to PID Controllers – Theory, Tuning and Application to Frontier Areas

*<sup>w</sup> Q S E wv v S R*

 

the m inputs and l outputs of the process. The vector xk is the state vector of the process at

the measurement noise and wk is called the process noise. It is assumed that they are zero mean, stationary white noise vector sequences and uncorrelated with the inputs uk. *nxn A R* is the system matrix, *nxm B R* is the input matrix, *lxn C R* is the output matrix while *lxm D R* is the direct feed-through matrix. The matrices *nxn Q R* , *nxl S R* and

In subspace identification it is typically assumed that the number of available data points goes to infinity, and that the data is ergodic. The main problem of identification is arranged

Given a large number of measurements of the input uk and the output yk generated by the unknown system described by equations (2.7)-(2.9). The task is to determine the order n of the unknown system, the system matrices A, B, C, D up to within a similarity

Step 1: Make a weighted projection of certain subspace generated from the data, to find an

Step 2: Retrieve the system matrices (A, B, C, D and Q, S, R) and from either this extended

*<sup>i</sup>* and/or an estimate *Xi*

*q q T pq*

with

as follows:

*<sup>p</sup> T T*

The vectors *mx*<sup>1</sup> *u R <sup>k</sup>* and *lx*<sup>1</sup>

discrete time instant k, *lx*<sup>1</sup>

*p*

*lxl R R* are the covariance matrices of the noise sequences wk and vk.

*<sup>i</sup>* ) or the estimated states.

transformation and an estimate of the matrices Q, S and R. Subspace identification algorithms always consist of two steps:

estimate of the extended observability matrix,

Fig. 1. Flow chart of subspace algorithm.

sequence *Xi* of the unknown system

observability matrix (

*k k kk y Cx Du v* (2.8)

(2.9)

of the state

*<sup>k</sup> y R* are the measurements at time instant k of, respectively,

*<sup>k</sup> v R* and *nx*<sup>1</sup> *w R <sup>k</sup>* are unobserved vector signals, vk is called

Based on the concept of sequential auto tuning (Shen & Yu, 1994) method each controller is designed in sequence. Let's consider a 2-by-2 MIMO system with a known pairing *y*1 1 *u* and *y*2 2 *u* under decentralized PI control (Figure 1). Initially, an ideal / biased relay is placed between 1 *y* and *u*<sup>1</sup> , while loop 2 is on manual (Figure 2a). Following the relayfeedback test, a controller can be designed from the ultimate gain and ultimate frequency. The next step is to perform relay-feedback test between 2 *y* and *u*2 while loop 1 is on automatic (Figure 2b). A controller can also be designed for loop 2 following the relayfeedback test. Once the controller on the loop 2 is put on automatic, another relay-feedback experiment is performed between 1 *y* and *u*<sup>1</sup> , (Figure 2c). Generally, a new set of tuning constants is found for the controller in loop 1. This procedure is repeated until the controller parameters converge. Typically, the controller parameters converge in 3 - 4 relay-feedback tests for 2 x 2 systems.

Identification and Control of Multivariable Systems – Role of Relay Feedback 113

From the identified step response models of Gp12,CL(s) and Gp22,CL(s), we can obtain their frequency response data and, by fitting them, we can get approximate low order models. Time domain modeling is obtained using equations (2.15) and (2.16) for 2x2 and 3x3 MIMO

22 2

*t t <sup>c</sup> <sup>n</sup> p p i n i k k yk e <sup>e</sup> y t*

1 1

<sup>21</sup> <sup>22</sup>

2 2 1 1

1 1

111

33

Wood and Berry (1973) (WB) reported a column for methanol-water separation with transfer

12.8 18.9 16.7 1 21 1

*e e x L s s x V e e s s*

6.6 19.4 10.9 1 14.4 1

The compositions of top (xD) and bottom (xB) products expressed in wt% of methanol are controlled variables. The reflux (L) and the reboiler (V) steam flow rates are the manipulated inputs are expressed in lb/min. time constants are in minutes. Feed flow rate is disturbance. Here the input variables are liquid (L) and vapour (V) flow rates (where as feed (F) flow rate is the load); outputs are distillate (xD) and bottom (xB) compositions. This plant given by

On applying least square algorithms to individual transfer function elements of an unknown 2x2 MIMO process (WB column) the estimated transfer function is obtained as shown in Table 1.The output (y) and input data (to original WB plant transfer function) are

On applying subspace algorithms to an unknown 2x2 MIMO process (WB column) the

Step 1: From the transfer function matrix State space representation matrices are calculated. Step 2: A, B, C and D matrices are simulate to get output data for a random input signal.

*s s <sup>B</sup>*

7 3

 

*s s*

*u*

22 2

*t t <sup>c</sup> n i p p <sup>i</sup>*

3 23 33

<sup>222</sup> <sup>1</sup>

*ttt*

21 2 22

21 2 22

*k k yk e <sup>e</sup>*

12 3

1

*e*

*a e a e a e*

*i*

<sup>33</sup>

*<sup>t</sup> <sup>c</sup> <sup>i</sup> <sup>p</sup> <sup>i</sup>*

*e*

*D*

Eq.(2.15) is considered as actual or real plant-model in present work.

used to form matrix. The parameters a1 and b1 were calculated using Eq.(2.6).

3 33 <sup>3</sup> <sup>2</sup> <sup>2</sup> <sup>1</sup>

 <sup>21</sup> <sup>22</sup> 21 22

21 21

*u u*

<sup>2</sup> 2 2

 

*e e*

3 23 33

3

*uuu*

*ppp*

222

*eee*

 

> 

> >

 

> 

(2.17)

(2.15)

(2.16)

*u u*

<sup>2</sup> 2 2 2 2 11 1

 

 

*e e*

process with FOPDT models using relay feedback test as:

*i*

 

**2.4 Process dynamics of example under study** 

33 3

*k k*

function as given below

following steps are followed

Fig. 2. Sequential method of tuning for 2x2 multivariable system. Steps are: (a) followed by (b) and followed by (c).

In order to proceed with sequential identification, it is necessary to derive closed-loop transfer functions for the above mentioned schemes. The following notations will be used for 2-by- 2 MIMO system:

$$\mathbf{g}\_{p}(\mathbf{s}) = \begin{pmatrix} \mathcal{S}\_{p11} & \mathcal{S}\_{p12} \\ \mathcal{S}\_{p21} & \mathcal{S}\_{p22} \end{pmatrix} \mathbf{G}\_{c}(\mathbf{s}) = \begin{pmatrix} \mathcal{G}\_{c1} & 0 \\ 0 & \mathcal{G}\_{c2} \end{pmatrix} \quad y(\mathbf{s}) = \begin{pmatrix} y\_{1} \\ y\_{2} \end{pmatrix} \text{ and } \boldsymbol{\mu}(\mathbf{s}) = \begin{pmatrix} \boldsymbol{\mu}\_{1} \\ \boldsymbol{\mu}\_{2} \end{pmatrix} \tag{2.10}$$

Thus, when perturbation is introduced in the second input u2, transfer functions for the input u2(s) are

$$y\_1 = \mathcal{G}\_{p12, \text{CL}}\left(s\right)u\_2\left(s\right) = -\frac{\mathcal{G}\_{p12}\left(s\right)}{1 + \mathcal{G}\_{p11}\left(s\right)\mathcal{G}\_{c1}\left(s\right)}u\_2\left(s\right) \tag{2.11}$$

$$\mathbf{g}\_{p2}(s) = \mathbf{G}\_{p22,CL}(s)\boldsymbol{\mu}\_2(s) = \left[\mathbf{g}\_{p22}(s) - \frac{\mathbf{g}\_{p21}(s)\mathbf{G}\_{c1}(s)\mathbf{g}\_{p12}(s)}{1 + \mathbf{g}\_{p11}(s)\mathbf{G}\_{c1}(s)}\right]\boldsymbol{\mu}\_2(s) \tag{2.12}$$

By applying the above identification method to the 2nd loop (by collecting output y2 for the change in input u1), we can obtain models for Gp12,CL(s) and Gp22,CL(s). Then, we have

$$\mathcal{G}\_2 = \mathcal{G}\_{p21,CL}(s)u\_1(s) = -\frac{\mathcal{G}\_{p21}(s)}{1 + \mathcal{G}\_{p22}(s)\mathcal{G}\_{c2}(s)}u\_1(s) \tag{2.13}$$

$$\mathbf{g}\_{1}(s) = \mathbf{G}\_{p11,CL}(s)\boldsymbol{u}\_{1}(s) = \left[\mathbf{g}\_{p11}(s) - \frac{\mathbf{g}\_{p21}(s)\mathbf{G}\_{c2}(s)\mathbf{g}\_{p12}(s)}{1 + \mathbf{g}\_{p22}(s)\mathbf{G}\_{c2}(s)}\right]\boldsymbol{u}\_{1}(s) \tag{2.14}$$

**u1**

Fig. 2. Sequential method of tuning for 2x2 multivariable system. Steps are: (a) followed by

**u2**

In order to proceed with sequential identification, it is necessary to derive closed-loop transfer functions for the above mentioned schemes. The following notations will be used

> 2 0

<sup>1</sup>

**G11**

**G21**

**G12**

**G22**

2 *y y s <sup>y</sup>* 

<sup>12</sup>

21 1 12

*g sG s* (2.11)

11 1 1 *p cp*

*p c g sG sg s*

*g sG s*

<sup>21</sup>

21 2 12

*g sG s* (2.13)

22 2 1 *p cp*

*p c g sG sg s*

*g sG s*

and <sup>1</sup>

**Y1**

**Y2** 

*u s*

2 *u*

(2.10)

(2.12)

(2.14)

*u* 

*c*

*G* 

Thus, when perturbation is introduced in the second input u2, transfer functions for the

11 1 1 *p*

*p c g s*

22 2 1 *p*

*p c g s*

1 12, 2 2

2 22, 2 22 2

*y s G su s g s u s*

By applying the above identification method to the 2nd loop (by collecting output y2 for the change in input u1), we can obtain models for Gp12,CL(s) and Gp22,CL(s). Then, we have

2 21, 1 1

1 11, 1 11 1

*y s G su s g s u s*

*y G su s u s*

*y G su s u s*

<sup>1</sup>

**Contro**

**Relay** 

*c*

*p CL*

*p CL p*

*p CL*

*p CL p*

*G s*

0 *c*

*G*

(b) and followed by (c).

for 2-by- 2 MIMO system:

*p*

input u2(s) are

11 12

**R1** 

**(c)**

**R2**

*g s g g*

21 22 *p p*

*p p*

 

*g g*

From the identified step response models of Gp12,CL(s) and Gp22,CL(s), we can obtain their frequency response data and, by fitting them, we can get approximate low order models. Time domain modeling is obtained using equations (2.15) and (2.16) for 2x2 and 3x3 MIMO process with FOPDT models using relay feedback test as:

$$y\_n = k\_{21} \left( 1 - e^{-\frac{t}{\tau\_{21}}} \left( \frac{2}{1 + e^{-\frac{p\_u}{2\tau\_{21}}}} \right) \right) - \frac{k\_{22} k\_{v2}}{\tau\_{i2}} \left( 1 + \left( \tau\_{i2} - \tau\_{22} \right) e^{-\frac{t}{\tau\_{22}}} \left( \frac{2}{1 + e^{-\frac{p\_u}{2\tau\_{22}}}} \right) \right) y\_n \left( t - 1 \right) \tag{2.15}$$

$$\begin{split} \boldsymbol{y}\_{n} &= k\_{21} \left( 1 - e^{-\frac{\boldsymbol{t}}{\tau\_{21}}} \left( \frac{2}{1 + e^{-\frac{\boldsymbol{p}\_{s}}{2\tau\_{21}}}} \right) \right) - \frac{k\_{22} k\_{c2}}{\tau\_{12}} \left( 1 + (\tau\_{i2} - \tau\_{22}) e^{-\frac{\boldsymbol{t}}{\tau\_{22}}} \left( \frac{2}{1 + e^{-\frac{\boldsymbol{p}\_{s}}{2\tau\_{21}}}} \right) \right) \\ &- \left[ \left[ 1 - \left( a\_{1} e^{-\frac{\boldsymbol{\mathcal{H}}\_{s}}{\tau\_{13}}} \left( \frac{2}{1 + e^{-\frac{\boldsymbol{p}\_{s}}{2\tau\_{13}}}} \right) + a\_{2} e^{-\frac{\boldsymbol{\mathcal{H}}\_{s}}{\tau\_{23}}} \left( \frac{2}{1 + e^{-\frac{\boldsymbol{p}\_{s}}{2\tau\_{21}}}} \right) \right) \right] - \left( \frac{2}{1 + e^{-\frac{\boldsymbol{p}\_{s}}{2\tau\_{21}}}} \right) \right] \\ &- \left[ \frac{k\_{33} k\_{c3}}{\tau\_{13}} \left[ 1 + (\tau\_{i3} - \tau\_{33}) e^{-\frac{\boldsymbol{\mathcal{H}}\_{s}}{\tau\_{33}}} \left( \frac{2}{1 + e^{-\frac{\boldsymbol{p}\_{s}}{2\tau\_{33}}}} \right) \right] \end{split} \tag{2.16}$$

#### **2.4 Process dynamics of example under study**

Wood and Berry (1973) (WB) reported a column for methanol-water separation with transfer function as given below

$$
\begin{bmatrix} \mathbf{x}\_D \\ \mathbf{x}\_B \end{bmatrix} = \begin{bmatrix} \frac{12.8e^{-s}}{16.7s+1} & \frac{-18.9e^{-3s}}{21s+1} \\ \frac{6.6e^{-7s}}{10.9s+1} & \frac{-19.4e^{-3s}}{14.4s+1} \end{bmatrix} \begin{bmatrix} L \\ V \end{bmatrix} \tag{2.17}
$$

The compositions of top (xD) and bottom (xB) products expressed in wt% of methanol are controlled variables. The reflux (L) and the reboiler (V) steam flow rates are the manipulated inputs are expressed in lb/min. time constants are in minutes. Feed flow rate is disturbance. Here the input variables are liquid (L) and vapour (V) flow rates (where as feed (F) flow rate is the load); outputs are distillate (xD) and bottom (xB) compositions. This plant given by Eq.(2.15) is considered as actual or real plant-model in present work.

On applying least square algorithms to individual transfer function elements of an unknown 2x2 MIMO process (WB column) the estimated transfer function is obtained as shown in Table 1.The output (y) and input data (to original WB plant transfer function) are used to form matrix. The parameters a1 and b1 were calculated using Eq.(2.6).

On applying subspace algorithms to an unknown 2x2 MIMO process (WB column) the following steps are followed

Step 1: From the transfer function matrix State space representation matrices are calculated.

Step 2: A, B, C and D matrices are simulate to get output data for a random input signal.

Identification and Control of Multivariable Systems – Role of Relay Feedback 115

The identified models and actual plant model are compared (Table-2.1). It is found that

**LEAST SQUARE SUBSPACE SEQUENTIAL** 

− − − − ⎡ ⎤ <sup>−</sup> ⎢ ⎥ + + − ⎣ ⎦ + +

12.799 18.9 16.7541 1 21.054 1 6.600 19.398 10.9505 1 14.448 1

Table 2.1. Actual and estimated multivariable transfer functions using different methods

After identifying the model structures and estimating process parameters of the models,

MIMO systems came into use in chemical industries as the processes were redesigned to improve efficiency. Multivariable control involves the objective of maintaining several controlled variables at independent set points. Interaction between inputs and output cause a manipulated variable to affect more than one controlled variable. The various control schemes studied here are the decentralized, centralized and decoupled systems. In decentralized structure, diagonal controllers are used. Hence they result in systems having n controllers. The centralized control systems have n x n controllers. In decoupled systems the process interactions are decoupled before they can actually reach and affect

Centralized control scheme is a full multivariable controller where the controller matrix is not a diagonal one. The decentralized control scheme is preferred over the centralized control scheme mainly because the control system has only n controlling n output variables, and the operator can easily understand the control loops. However, the design methods of such decentralized controllers require first pairing of input-output variables, and tuning of controllers requires trial and error steps. The centralized control system requires n x n controllers for controlling n output variables using n manipulated variables. But if we are calculating the control action using a computer, then this problem of requiring n x n controllers does not exist. The advantage of the centralized controller is easy to tune even with the knowledge of the steady state gain matrix alone, multivariable PI controllers can be

For the centralized structure, Internal model control-proportional integral tuning is adopted, based on studies on the studies and recommendations of Reddy et al (1997) on the design of centralized PI controllers for a Multi-stage flash desalination plant using Davison,

The IMC-PID tuning relations are used in tuning the controller. For a first order system of

, the PI controller settings are as follows:

7 3

*s s s s e e s s e e s s*

3

11, 2

*CL*

*CL*

22, 2

3

− ⎛ ⎞ <sup>+</sup> <sup>=</sup> ⎜ ⎟ + + ⎝ ⎠ <sup>+</sup>

*s*

− + ⎛ ⎞ <sup>=</sup> ⎜ ⎟ + + ⎝ ⎠ <sup>+</sup>

*s*

−

*e s <sup>G</sup> s s <sup>s</sup>*

*e s <sup>G</sup> s s <sup>s</sup>*

6.4 60 1 42.25 11.7 1 44 1 9.655 0.4774 1 1453 58.15 1 2.741 1

**COLUMN ESTIMATED TRANSFER FUNCTION** 

3

7 3 12.692 18.89 16.41 1 21.46 1 6.4 19.6 10.95 1 14.481 1

*s s s s e e s s e e s s*

next work is to select a suitable control strategy for the process.

− − − − ⎡ ⎤ <sup>−</sup> ⎢ ⎥ + + − ⎣ ⎦ + +

subspace identification method gives better result/

**ACTUAL WB** 

7 3 12.8 18.9 16.7 1 21 1 6.6 19.4 10.9 1 14.4 1

the processes.

easily designed.

the form <sup>1</sup>

*Ds <sup>p</sup> k e s*

Maciejowski and Tanttu-Lieslehto methods.

**3.1 Centralized structure** 

*s s s s e e s s e e s s*

− − − − ⎡ ⎤ <sup>−</sup> ⎢ ⎥ + + − ⎣ ⎦ + +

3

**3. Different control strategies** 

Step 3: From the output and input data Henkel matrix are formed and LQ decomposition method is used to spilt the matrix

Step 4: Then Singular value decomposition method is used to estimate A, B, C and D matrices.

Step 5: From estimated matrices the transfer function were found.

Fig. 3. Comparison of responses between actual (solid) and identified (Sequential identification, dashed line) models of WB column

Mostly, the purpose of identification of transfer functions is to design controller for the system in order to achieve desired performance. Three methods of identifications (two in openloop mode and the other in closed-loop mode) are used to identify the two-input-twooutput process, WB column. Least square and subspace methods have been used to identify the process in openloop and sequential identification technique is used to estimate the process in closedloop.

The identified models and actual plant model are compared (Table-2.1). It is found that subspace identification method gives better result/


Table 2.1. Actual and estimated multivariable transfer functions using different methods

After identifying the model structures and estimating process parameters of the models, next work is to select a suitable control strategy for the process.
