Decoupling internal model control for non-square processes based on equivalent transfer function

Abstract

A novel design method of decoupling internal model control is proposed for non-square processes with multiple time delays that are often encountered in complicated industrial processes. The method can obtain a realizable decoupling controller of non-square processes with more inputs than outputs by inserting some compensated terms, which are derived analytically. Meanwhile, based on the relative normalized gain array, an equivalent transfer function matrix is introduced to approximate the pseudo-inverse of the process transfer function matrix, which makes the design of decoupling internal model control simple and easy to calculate. Filters are added to the control structure to improve the robustness. Simulation results have proved the effectiveness and reliability of the proposed method.

Keywords

Decoupling equivalent transfer function (ETF)internal model control (IMC)non-square processes relative normalized gain array (RNGA)

Introduction

In a complex industrial production process, it is always necessary to use multiple controlled variables and multiple manipulated variables to describe a controlled system, and we often encounter a class of systems in which the number of input variables and output variables are usually not equal, which are known as non-square systems. Although non-square processes in industry are very common, the literature on controller design of non-square system is still sparse. A traditional approach towards the control of non-square processes is to square up or square down the systems to square systems through adding or deleting appropriate inputs (manipulated variables) or outputs (controlled variables). Thereafter, most well developed control strategies for square systems (where the number of inputs and outputs are equal) can be applied directly.

However, adding inputs or outputs to be measured can be costly, while removing inputs or outputs will lead to a poor performance because of the information that is neglected and may even make the system unstable through the possible introduction of right half plane (RHP) poles. Therefore, a promising approach to the above problems is to design a controller for the original non-square system (Morari et al., 1985). Treiber and Hoffman (1986) have showed that the performance achieved by using the full non-square model of a vacuum distillation unit is superior to the result obtained by employing a square subsystem, and the non-square system is illustrated to possess a better control character and stronger robustness compared with any of the possible squared-down systems designed by an independent robust controller (Loh and Chiu, 1997).

Internal model control (IMC) has been proved an effective method to resolve time-delay and decoupling problems (Garcia and Morari, 1982, 1985). On the basis of the IMC principle, the design of a decoupling controller for a square system is generally based on the inverse matrix of original system (Liu et al., 2006). However, considering non-square processes, IMC cannot be applied directly, as it cannot obtain the traditional sense of inversion. In addition, it is also difficult to select the non-minimum phase part of a transfer matrix with multi-delay for optimal performance. To overcome these questions, a modified IMC method for non-square processes with multiple time delays and right-half-plant zeros is proposed (Chen et al., 2008). In order to improve dynamic quality and robustness of non-square systems with first-order time delay, the IMC method is used to design Smith delay compensation to compensate for shortcomings caused by static decoupling, and also overcome the impact of the model error on the system performance by model approximation (Chen et al., 2011). Because IMC is a method of realizing the controller design, relying heavily on the inverse or generalized inverse of process model, it is important to research a simple and easy method to calculate the inverse or generalized inverse.

In this paper, a novel design method of decoupling IMC is proposed for non-square processes with different time delays. The method can realize decoupling of non-square processes by inserting some compensated terms, which are derived analytically. Meanwhile, an equivalent transfer function (ETF) matrix was used to approximate the pseudo-inverse of the process transfer function matrix, which makes the design of decoupling IMC simple and easy to calculate. The simulation results show the effectiveness of the proposed method. However, it is worth mentioning that the proposed method is specific to plants with more inputs than outputs, as IMC is based on the assumption of perfect control, but for systems with more outputs than inputs, perfect control is in the least-square sense. Obviously, the proposed method cannot be applied to control it perfectly; how to achieve the control of that should be further studied.

The design of decoupling internal model control

The decoupling IMC structure

A block diagram of IMC control structure is shown in Figure 1, where $Y_{p} (s)$ and $Y_{m} (s)$ are the outputs of the process and its model, respectively; $U (s)$ is the input of the process; $R (s)$ represents the set-point; $D (s)$ is the external disturbance; $G_{p} (s)$ and $G_{m} (s)$ represent the transfer functions of the process and the process model, respectively, with $n$ inputs and $m$ outputs $(m \leq n)$ given by the equation:

G_{p} (s) = G_{m} (s) = [\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 n} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 n} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{m 1} (s) & g_{m 2} (s) & \dots & g_{mn} (s) \end{matrix}]

(1)

where $g_{ij} (s) = g_{ij 0} (s) e^{- θ_{ij} s}$ is the transfer function from jth input to ith output, $g_{ij 0} (s)$ is strictly proper and stable transfer function, and $θ_{ij} \geq 0$ is the corresponding time delay of process transmission. $G_{IMC} (s)$ is the transfer function matrix of the decoupling IMC with $m$ inputs and $n$ outputs. In order to improve the robustness, a low-pass filter $G_{f} (s)$ is added in the control structure. The key to decoupling the IMC structure is the design of $G_{IMC} (s)$ , which considers the decoupling and the controller together.

The design of internal model controller

The design procedure of the controllers $G_{IMC} (s)$ and $G_{f} (s)$ will be detailed in this section.

From Figure 1, the closed-loop transfer matrix $H (s)$ between $Y_{p} (s)$ and $R (s)$ is derived as:

H (s) = G_{P} (s) G_{IMC} (s) {[I + (G_{P} (s) - G_{m} (s)) G_{f} (s) G_{IMC} (s)]}^{- 1}

(2)

Assuming $G_{P} (s) = G_{m} (s)$ , $H (s)$ is reduced to:

H (s) = G_{P} (s) G_{IMC} (s) = G_{m} (s) G_{IMC} (s)

(3)

According to the traditional IMC theory, the system is decoupled if and only if $H (s)$ is decoupled (diagonal and non-singular), which is formed as:

H (s) = diag {h_{ii}}; h_{ii} \neq 0; i = 1, 2, \dots m

(4)

Then, from (3), the decoupling internal model controller can be written as:

G_{IMC} (s) = G_{m}^{- 1} (s) H (s)

It is obvious that $G_{IMC} (s)$ can be obtained if $G_{m}^{- 1} (s)$ and $H (s)$ are determined. The design procedure of $G_{IMC} (s)$ can be summarized in the following steps.

First, determine $G_{m}^{- 1} (s)$ . For non-square systems, the traditional sense of the inverse does not exist due to the process model being a non-square matrix. Hence the Moore–Penrose pseudo-inverse is introduced, written as $G_{m}^{*} (s)$ . For matrix $G_{m} (s)$ , Moore–Penrose pseudo-inverse is calculated as (Zhang, 2007):

G_{m}^{*} (s) = G_{m}^{H} (s) {[G_{m} (s) G_{m}^{H} (s)]}^{- 1}

(6)

where $G_{m}^{H} (s)$ is the Hermitian matrix of $G_{m} (s)$ .

Next, choose the appropriate $H (s)$ . Consider the process model $G_{m} (s)$ , which can be factored as:

G_{m} (s) = G_{m +} (s) G_{m -} (s)

(7)

where $G_{m -} (s)$ is stable and minimum phase with no predictors, while $G_{m +} (s)$ contains all the dead time and RHP zeros of $G_{m} (s)$ . The traditional internal model controller is designed as the inverse of the model cascaded with a filter $F (s)$ , which used to compensate for the model mismatch and ignore error of non-minimum, formed as:

G_{IMC} (s) = G_{m -}^{- 1} (s) F (s)

(8)

Then (3) becomes to

H (s) = G_{m +} (s) F (s)

(9)

Here, $F (s) = diag {\frac{1}{{(λ_{i} s + 1)}^{n_{i}}}, i = 1, 2, \dots, m}$ is an $m$ -dimension filter matrix, which guarantees that the controller is rational. The filter order $n_{i}$ is selected to be sufficiently large to guarantee the internal model controller proper, while $λ_{i}$ is an adjustable parameter that determines the speed of response. $G_{m +} (s)$ is the part of non-minimum phase of $G_{m} (s)$ , which contains all the time delays (their inversion would require a prediction) in the stable process and RHP zeros (their inversion would yield an unstable controller), and satisfies $| G_{m +} (j ω) | = I$ . For single-input and single-output systems, it is easy to factor $G_{p} (s)$ , and $G_{m +} (s)$ is obtained, but for multivariable systems with different delays, $G_{m +} (s)$ is very difficult to obtain. However, its general form has been derived by Chen et al. (2008), which is a diagonal matrix, and the ith element of the diagonal matrix corresponds to the ith column of $G_{m}^{*} (s)$ . The form of $G_{m +} (s)$ is:

G_{m +} (s) = diag {e^{- θ_{Ki} s} Π_{1}^{D_{i}} {(\frac{- s + z_{p}}{s + z_{p}^{*}})}^{W_{i}}, i = 1, 2, \dots, m}

(10)

where $θ_{Ki}$ is the maximum predictive term in the ith column of $G_{m}^{*} (s)$ ; $z_{p} {p = 1, 2, \dots D_{i}}$ is the RHP pole in corresponding column, and $z_{p}^{*}$ is the conjugate complex of $z_{p}$ ; $W_{i}$ is the number of the same pole while $D_{i}$ is the number of different poles existing in the ith column of $G_{m}^{*} (s)$ and always meets the condition of $D_{i} < m$ . The essence of (10) is a compensator, to compensate for the predictive term and RHP poles introduced by $G_{m}^{*} (s)$ . From the above, the non-square decoupling internal model controller is therefore calculated by

G_{IMC} = G_{m}^{*} (s) G_{m +} (s) F (s)

(11)

The design of feedback filter

In industry, model inaccuracies widely exist, and thus it is necessary to include filters in the feedback signals, in order to make the system maintain the original performance. When uncertainty exists in the control system, the system robust stability can be guaranteed by increasing the value of control parameter $λ_{i}$ . However, the corresponding response speeds of the system output will be decreased. Hence, in order to improve the robustness while obtaining a fast response speed, a first-order filter $G_{f} (s)$ (Normey-Rico et al., 1997) is recommended for addition in the original IMC structure, as shown in Figure 1.

For the $m \times n$ non-square process $G_{p} (s)$ , the first-order filter $G_{f} (s)$ can be designed as:

G_{f} (s) = diag {\frac{1}{T_{f} s + 1}, i = 1, 2, \dots, m}

(11)

where $T_{f} = θ_{Fi} / α$ is time constant; $α$ is a constant and $θ_{Fi}$ denotes the time delay term of the ith control loop. In terms of (10) and (11), $θ_{Fi}$ satisfies the equation $θ_{Fi} = θ_{Ki}$ .

Using the cut-frequency of the filter as a third of the frequency of these oscillations, $T_{f}$ can be chosen as $T_{f} = θ_{Fi} / 2$ . This choice gives a good solution for dead-time errors up to 30%. If necessary, in some particular cases, the desired performance and robustness can be obtained by tuned $α$ .

Internal model controller design based on equivalent transfer function

The internal model controller designed by the above procedure usually has a complicated expression, especially when the system dimension is high. In this section, the ETF method is employed to find a simple description for $G_{m}^{*} (s)$ that is easy to calculate. Besides providing a very simple approach for obtaining an internal model controller, it is also a good approximation to the above theoretical controller.

For multivariate systems, ETF refers to open equivalent transfer function of a loop in the case of other loops being closed. Then how to obtain the ETF model is an important issue. Relative normalized gain array (RNGA) was provided by Xiong et al. (2007), and using the array, the change of steady state gain, time constant and delay time can be obtained when other loops are closed, and then the approximate model of ETF can also be obtained based on original transfer function (Shen et al., 2010). In this paper, it is expanded to non-square processes.

Equivalent transfer function for non-square processes

For non-square processes, by assuming that the process is under perfect control, relative gain array (RGA) and RNGA is extended to non-square multivariable systems (Chang and Yu, 1990). The non-square relative gain (NRG) is defined as the ratio of the open-loop gain to the closed-loop gain. Then the open-loop gain is easy to obtain by the transfer function, while the closed-loop gain is not that straightforward. Under the assumption of ‘perfect control’, which is to say the designed controller makes all outputs $y_{k} (s), k \neq i$ , other than the output $y_{i} (s)$ , which are kept at set points, a simple computational algorithm for non-square RGA (NRGA) can be obtained.

Similar to the square RGA, the NRG and the NRGA are defined as

λ_{ij}^{N} = \frac{k_{ij}}{{\hat{k}}_{ij}}, i = 1, 2, \dots, m; j = 1, 2, \dots, n

(12)

and

Λ^{N} = [\begin{matrix} λ_{11}^{N} & λ_{12}^{N} & \dots & λ_{1 n}^{N} \\ λ_{21}^{N} & λ_{22}^{N} & \dots & λ_{2 n}^{N} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ λ_{m 1}^{N} & λ_{m 2}^{N} & \dots & λ_{mn}^{N} \end{matrix}] = K^{N} \otimes {(K^{N})}^{* T}

(13)

respectively. In form (12), $k_{ij}$ is the open-loop gain when all loops are open, ${\hat{k}}_{ij}$ is the closed-loop gain when all loops other than the loop i–j are closed and under perfect control. $K^{N} = {[k_{ij}]}_{m \times n}$ represents the open-loop gain array, and ${(K^{N})}^{* T}$ is the corresponding transposed matrix of pseudo-inverse; the operator $\otimes$ is the Hadamard product.

To describe the dynamic properties, the normalized gain $k_{N,_{ij}}^{N}$ is defined as

k_{N, ij}^{N} ≜ \frac{k_{ij}}{σ_{ij}}, i = 1, 2, \dots, m, j = 1, 2, \dots, n

(14)

In the above expressions, $σ_{ij}$ is the average residence time (ART) of loops i–j, defined as $σ_{ij} = \frac{\int_{0}^{\infty} [y (\infty) - y (t)] dt}{y (\infty)}$ , $y (t)$ is open step response function and $y (\infty)$ is open steady state value of step response.

The non-square RNGA (NRNGA) can be written as

\begin{matrix} Φ^{N} & = [\begin{matrix} ϕ_{11} & ϕ_{12} & \dots & ϕ_{1 n} \\ ϕ_{21} & ϕ_{22} & \dots & ϕ_{2 n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ϕ_{m 1} & ϕ_{m 2} & \dots & ϕ_{mn} \end{matrix}] \\ = [\begin{matrix} \frac{k_{N, 11}^{N}}{{\hat{k}}_{N, 11}^{N}} & \frac{k_{N, 12}^{N}}{{\hat{k}}_{N, 12}^{N}} & \dots & \frac{k_{N, 1 n}^{N}}{{\hat{k}}_{N, 1 n}^{N}} \\ \frac{k_{N, 21}^{N}}{{\hat{k}}_{N, 21}^{N}} & \frac{k_{N, 22}^{N}}{{\hat{k}}_{N, 22}^{N}} & \dots & \frac{k_{N, 2 n}^{N}}{{\hat{k}}_{N, 2 n}^{N}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{k_{N, m 1}^{N}}{{\hat{k}}_{N, m 1}^{N}} & \frac{k_{N, m 2}^{N}}{{\hat{k}}_{N, m 2}^{N}} & \dots & \frac{k_{N, mn}^{N}}{{\hat{k}}_{N, mn}^{N}} \end{matrix}] = K_{N}^{N} \otimes {(K_{N}^{N})}^{* T} \end{matrix}

(15)

Furthermore, the relative average residence time $γ_{ij}$ , which is defined as the ratio of loop i–j average residence time between when other loops are closed and when other loops are open, and the non-square relative average residence time array (NRARTA) can be obtained by

Γ^{N} = [\begin{matrix} γ_{11} & γ_{12} & \dots & γ_{1 n} \\ γ_{21} & γ_{22} & \dots & γ_{2 n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ γ_{m 1} & γ_{m 1} & \dots & γ_{mn} \end{matrix}]

(16)

where

γ_{ij} = \frac{{\hat{σ}}_{ij}}{σ_{ij}}, i = 1, 2, \dots, m, j = 1, 2, \dots, n

(17)

From (13)–(15), NRARTA is calculated by

Γ^{N} = Φ^{N} ⊙ Λ^{N}

(18)

where ⊙ indicates element-by-element division.

Previous investigations have shown that almost all industry processes can be simplified to a first-order-plus-time-delay (FOPTD) model. For simplicity, consider the FOPTD model for a transfer function

g_{ij} (s) = \frac{k_{ij}}{τ_{ij} s + 1} e^{- θ_{ij} s}, i = 1, 2, \dots, m, j = 1, 2, \dots, n

(19)

The average residence time is formed as $σ_{ij} = τ_{ij} + θ_{ij}$ , and the relative average residence times $γ_{ij}$ represents the dynamic changes of the transfer function $g_{ij} (s)$ when other loops are closed. Using the definition of NRARTA, it is possible to write

\begin{matrix} {\hat{σ}}_{ij} = γ_{ij} \times σ_{ij} = γ_{ij} \times τ_{ij} + γ_{ij} \times θ_{ij}, \\ i = 1, 2, \dots, m, j = 1, 2, \dots, n \end{matrix}

(20)

Then the changes of transfer function elements when other loops closed can be uniquely determined. An ETF element when other loops are closed can be approximated by a transfer function element and has the same form as the open-loop transfer function element, where the steady state gain, the time constant and the time delay are scaled by NRGA, NRNGA and NRARTA, respectively, i.e.

{\hat{g}}_{ij} (s) = \frac{{\hat{k}}_{ij}}{{\hat{τ}}_{ij} s + 1} e^{- {\hat{θ}}_{ij} s} = \frac{k_{ij}}{λ_{ij}} \frac{1}{γ_{ij} τ_{ij} s + 1} e^{- γ_{ij} θ_{ij} s}

(21)

where ${\hat{g}}_{ij} (s)$ is ETF element. When ${\hat{g}}_{ij} (s)$ are calculated for all of the input/output combinations of a multivariable process, the results are in an array form, and the ETF array $\hat{G} (s)$ is obtained.

The relationship between equivalent transfer function array and system transfer function matrices

For square systems, Xiong et al. (2007) established the relationship between $G^{- 1} (s)$ and the ETF matrix, where

G^{- 1} (s) = {\hat{G}}^{T} (s)

(22)

For a non-square system, the non-square dynamic RGA (NDRGA) is defined as

\begin{matrix} Λ^{N} (s) = G (s) \otimes \hat{G} (s) & = [\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 n} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 n} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{m 1} (s) & g_{m 2} (s) & \dots & g_{mn} (s) \end{matrix}] \\ \otimes [\begin{matrix} \frac{1}{{\hat{g}}_{11} (s)} & \frac{1}{{\hat{g}}_{12} (s)} & \dots & \frac{1}{{\hat{g}}_{1 n} (s)} \\ \frac{1}{{\hat{g}}_{21} (s)} & \frac{1}{{\hat{g}}_{22} (s)} & \dots & \frac{1}{{\hat{g}}_{2 n} (s)} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{1}{{\hat{g}}_{m 1} (s)} & \frac{1}{{\hat{g}}_{m 2} (s)} & \dots & \frac{1}{{\hat{g}}_{mn} (s)} \end{matrix}] \end{matrix}

(23)

Under the assumption of perfect control, the NDRGA can also be calculated by

\begin{matrix} Λ^{N} (s) = G (s) \otimes G^{*} (s) & = [\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 n} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 n} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{m 1} (s) & g_{m 2} (s) & \dots & g_{mn} (s) \end{matrix}] \\ \otimes {[\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 n} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 n} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{m 1} (s) & g_{m 2} (s) & \dots & g_{mn} (s) \end{matrix}]}^{* T} \end{matrix}

(24)

Comparing (24) and (25):

[\begin{matrix} \frac{1}{{\hat{g}}_{11} (s)} & \frac{1}{{\hat{g}}_{12} (s)} & \dots & \frac{1}{{\hat{g}}_{1 n} (s)} \\ \frac{1}{{\hat{g}}_{21} (s)} & \frac{1}{{\hat{g}}_{22} (s)} & \dots & \frac{1}{{\hat{g}}_{2 n} (s)} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{1}{{\hat{g}}_{m 1} (s)} & \frac{1}{{\hat{g}}_{m 2} (s)} & \dots & \frac{1}{{\hat{g}}_{mn} (s)} \end{matrix}] = {[\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 n} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 n} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{m 1} (s) & g_{m 2} (s) & \dots & g_{mn} (s) \end{matrix}]}^{* T}

(25)

Taking the transpose of both sides of (25), we obtain

G^{*} (s) = {\hat{G}}^{T} (s)

(26)

From the above, substituting (26) into (11), the decoupling internal model controller is calculated by

G_{IMC} (s) = {\hat{G}}_{m}^{T} (s) G_{m +} (s) F (s)

(27)

By using ETF, the design of $G_{IMC} (s)$ will be not only simple and easy to calculate, but also have high-dimensional capacity.

Furthermore, carrying on the analysis on the impact of ETF model approximation, when adding decoupling controller, the closed-loop transfer matrix $\bar{H} (s)$ is expressed as

\bar{H} (s) = {[{\bar{h}}_{ij} (s)]}_{m \times m} = G_{m} (s) {\hat{G}}_{m}^{T} (s) G_{m +} (s) F (s) = G_{m +} (s) F (s)

(28)

Assume the best loops pairing are i–i. First considering the diagonal elements, to describe them clearly, we obtain the further expansion of (28)

\begin{matrix} {\bar{h}}_{ii} (s) = g_{mi +} (s) f_{i} (s) \\ [\frac{k_{i 1}}{{\hat{k}}_{i 1}} \frac{{\hat{τ}}_{i 1} s + 1}{τ_{i 1} s + 1} e^{- (θ_{i 1} - {\hat{θ}}_{i 1}) s} + \dots + \frac{k_{im}}{{\hat{k}}_{im}} \frac{{\hat{τ}}_{im} s + 1}{τ_{im} s + 1} e^{- (θ_{im} - {\hat{θ}}_{im}) s}] \\ = {\bar{h}}_{ii}^{1} + \dots + {\bar{h}}_{ii}^{m} \end{matrix}

(29)

where $λ_{ij}^{N} = \frac{k_{ij}}{{\hat{k}}_{ij}}$ may be negative, then the control system will appear to have an inverse response, which leads to poor dynamic performance. Therefore, in order to overcome this defect, the following constraint condition must be added to correct the ETF

θ_{i j^{*}} - {\hat{θ}}_{i j^{*}} \leq θ_{ij} - {\hat{θ}}_{ij}, j^{*} \neq j

(30)

where $θ_{i j^{*}}$ represents the time delay of the best loops i–j^*, which are determined by NRGA, and $λ_{i j^{*}}^{N}$ is positive and close to $1$ . The constraint ensures the response of the best loops i–j^* before other loops on the output $y_{i} (s)$ , and makes the system obtain satisfactory performance and stability.

Then, considering off-diagonal elements, it can be expressed that

{\bar{h}}_{ij} (s) = g_{mi +} (s) f_{i} (s) [\sum_{l = 1}^{n} \frac{k_{il}}{{\hat{k}}_{jl}} \frac{{\hat{τ}}_{il} s + 1}{τ_{jl} s + 1} e^{- (θ_{il} - {\hat{θ}}_{jl}) s}]

(31)

It has been proved that $\sum_{l = 1}^{n} \frac{k_{il}}{{\hat{k}}_{jl}} = 0$ (Shen et al., 2010), so in order to improve the decoupling, the following constraint is expected to be established

θ_{i 1} - {\hat{θ}}_{j 1} = θ_{i 2} - {\hat{θ}}_{j 2} = \dots = θ_{in} - {\hat{θ}}_{jn}

(32)

Taken together, for the best control effect, both conditions of (30) and (32) are expected to be established. However, if that is impossible, it needs to meet the constraint (32) first of all. Then update the time delay ${\hat{θ}}_{ij}$ of the ETF array according the above constraints.

Robust analysis and performance index

As a matter of fact, modelling errors and system uncertainties are inevitable and cannot be ignored when modelling a real plant, and the present of uncertainties usually causes instability and poor performance of control systems, thereby leading to very complex dynamical behaviours. It is necessary to analyse the robust stability of the proposed control system, and in this subsection, we analyse the system robustness.

As for the standard assumptions in robustness analysis, the nominal IMC system is stable. Suppose that the control system with the additive uncertainty is denoted by

G_{p} (s) = G_{m} (s) + Δ (s)

(33)

where $Δ (s) = diag {Δ_{i} (s), i = 1, 2, \dots, m}$ is uncertain weight, and it is assumed to be stable and satisfies $‖ Δ (j ω) ‖ \leq 1, \forall ω$ . Figure 2 describes the control system with the additive uncertainty.

For convenience and for analysing robust stability, the general $M - Δ$ structure is achieved in Figure 2. For a nominal system, the transfer function $M (s)$ is expressed as:

M (s) = - G_{IMC} (s) G_{f} (s) Δ (s)

(34)

According to the robust stability theorem (Morari and Zafiriou, 1989), the proposed control system is robustly stable for $Δ (s)$ if and only if

\bar{σ} (M) = \bar{σ} (G_{IMC} G_{f} (s) Δ) < 1, \forall ω

(35)

where $\bar{σ} (•)$ denotes the maximum singular value.

If $Δ (s)$ is known, Equation (35) can be employed to evaluate the robust stability of the proposed control system. In other words, if the magnitude plots of the left side of (35) is less than 1 $\forall ω$ , the system is robustly stable.

Simulation results

Two plant models are considered to demonstrate the effectiveness of the proposed design methods.

Example 1

Consider the shell control problem with three inputs and two outputs (Ogunnaike and Ray, 1979). The transfer function matrix is given as

G_{p} (s) = [\begin{matrix} \begin{matrix} \frac{4.05 e^{- 81 s}}{50 s + 1} & \frac{1.77 e^{- 84 s}}{60 s + 1} \end{matrix} & \begin{matrix} \frac{5.88 e^{- 81 s}}{50 s + 1} \end{matrix} \\ \begin{matrix} \frac{5.39 e^{- 54 s}}{50 s + 1} & \frac{5.72 e^{- 42 s}}{60 s + 1} \end{matrix} & \begin{matrix} \frac{6.9 e^{- 45 s}}{40 s + 1} \end{matrix} \end{matrix}]

Step 1. Calculating the ETF array ${\hat{G}}_{m} (s)$ .

The open-loop gain, time constant, and time delay of the system are given as

\begin{matrix} K^{N} = [\begin{matrix} 4.05 & 1.77 & 5.88 \\ 5.39 & 5.72 & 6.9 \end{matrix}], T = [\begin{matrix} 50 & 60 & 50 \\ 50 & 60 & 40 \end{matrix}], \\ Θ = [\begin{matrix} 81 & 84 & 81 \\ 54 & 42 & 45 \end{matrix}] \end{matrix}

A simple calculation gives

Λ^{N} = K^{N} \otimes {(K^{N})}^{* T} = [\begin{matrix} 0.3203 & - 0.5946 & 1.2744 \\ - 0.0170 & 1.5733 & - 0.5563 \end{matrix}]

Φ^{N} = K_{N}^{N} \otimes {(K_{N}^{N})}^{* T} = [\begin{matrix} 0.6662 & 0.9585 & 0.9585 \\ - 0.3174 & 1.6153 & - 0.2978 \end{matrix}]

Γ^{N} = Φ^{N} ⊙ Λ^{N} = [\begin{matrix} - 4.8992 & 1.0377 & 4.8615 \\ 5.7741 & - 0.6536 & - 4.1205 \end{matrix}]

which results in the following ETF parameters

\begin{matrix} \hat{K} = K^{N} ⊙ Λ^{N} = [\begin{matrix} 12.6459 & - 2.9766 & 4.6140 \\ - 316.9046 & 3.6356 & - 12.4033 \end{matrix}] \\ \hat{Θ} = Γ^{N} \otimes Θ = [\begin{matrix} 81 & 81 & 84 \\ 45 & 45 & 42 \end{matrix}] \end{matrix}

\hat{T} = Γ^{N} \otimes T = [\begin{matrix} 104.0130 & 63.0397 & 37.6080 \\ 933.1931 & 61.6007 & 21.4159 \end{matrix}]

Substituting these ETF parameters into the original transfer function matrix, the ETF matrix is obtained as

{\hat{G}}_{m}^{T} (s) = [\begin{matrix} \frac{104.0130 s + 1}{12.6459} e^{81 s} & \frac{933.1931 s + 1}{- 316.9046} e^{45 s} \\ \frac{63.0397 s + 1}{- 2.9766} e^{81 s} & \frac{61.6007 s + 1}{3.6356} e^{45 s} \\ \frac{37.6080 s + 1}{4.6140} e^{84 s} & \frac{21.4159 s + 1}{- 12.4033} e^{42 s} \end{matrix}]

Step 2. Deigning decoupling internal model controller $G_{IMC} (s)$ .

According to Equation (10), the form of $G_{m +} (s)$ is selected as

G_{m +} (s) = diag {e^{- 84 s}, e^{- 45 s}}

In view of the robustness and response time, filters are selected as

F (s) = diag {\frac{1}{75 s + 1}, \frac{1}{35 s + 1}}

Thus, the final decoupling controller is obtained as

G_{IMC} (s) = [\begin{matrix} \frac{104.0130 s + 1}{12.6459 (75 s + 1)} e^{- 3 s} & \frac{933.1931 s + 1}{- 316.9046 (35 s + 1)} \\ \frac{63.0397 s + 1}{- 2.9766 (75 s + 1)} e^{- 3 s} & \frac{61.6007 s + 1}{3.6356 (35 s + 1)} \\ \frac{37.6080 s + 1}{4.6140 (75 s + 1)} & \frac{21.4159 s + 1}{- 12.4033 (35 s + 1)} e^{- 3 s} \end{matrix}]

Step 3. Designing of filter $G_{F} (s)$ . The feedback loop filters are considered

G_{F} (s) = diag {\frac{1}{42 s + 1}, \frac{1}{27 s + 1}}

At a perfect matching of a process plant and model, the set-point responses of the system for a step change in $r_{1}$ and $r_{2}$ are shown in Figures 3 and 4, where the solid line represents the method in this paper, the dashed line, dot dash line and red solid line represent the methods proposed by Sarma and Chidambaram (2005), that by Guo (Jin et al., 2010) and the squared method, respectively, while the disturbance responses are shown in Figures 5 and 6, where distractions introduced the system for 0.5 and −0.5 step change at 500s in $d_{1}$ and $d_{2}$ , respectively.

It can be seen from the simulation results that the method proposed in this paper makes two outputs nearly decoupled and driven to their steady-state values. Compared with those given by the other three methods, the proposed method shows smaller overshoot, faster tracking features and disturbance rejecting performance.

In order to verify the robustness of the proposed method, we increase the gains and the dead times of the process elements by 20%, and, conversely, the time constants are decreased by 20% and the process becomes:

{\tilde{G}}_{p} (s) = [\begin{matrix} \frac{4.86 e^{- 97.2 s}}{40 s + 1} & \frac{2.124 e^{- 100.8 s}}{48 s + 1} & \frac{5.88 e^{- 97.2 s}}{40 s + 1} \\ \frac{6.468 e^{- 64.8 s}}{40 s + 1} & \frac{6.86 e^{- 50.4 s}}{48 s + 1} & \frac{8.28 e^{- 54 s}}{32 s + 1} \end{matrix}]

Under the condition of the set-point and the controller parameters being constant, the set-point responses are shown in Figures 7 and 8. From the results, we conclude that the proposed method still offers good tracking features and strong robustness in the case of the given perturbation.

To evaluate the decoupling result and dynamic performance comprehensively, the integral of squared error index (ISE; Lin et al., 1999), which considers the sum of the ISE values in response and interaction, is chosen as the performance index. To make a performance comparison among these four methods, the matching ISE values corresponding to the set-point responses and the mismatching with perturbations are shown in Table 1. From the ISE values, it is observed that the proposed method gives slower ISE values and shows no apparent variation during process model mismatches. This is a clear indication that the method proposed in this paper demonstrates good robustness.

Table 1.

Integral of squared error index (ISE) values when process plant matches and mismatches with process model.

Method		ISE_y1–r1	ISE_y2–r1	ISE_y1–r2	ISE_y2–r2	Sum of ISE value
Match with process model
	Proposed method	46.5900	0.1914	0.0228	22.3826	69.1868
	Sama’s method	109.1022	82.3699	44.4929	46.8515	282.8165
	Guo’s method	50.6014	36.1234	0.1758	44.6217	131.5223
	Squared method	72.3218	41.6950	3.3412	45.9312	163.2892
Mismatch with perturbation of +20% in gains, the dead times and −20% in time constants
	Proposed method	47.8865	0.3971	0.0549	23.3894	71.7279
	Sama’s method	308.3789	509.3730	143.8811	237.3390	1198.972
	Guo’s method	62.9583	48.2591	0.3786	46.5188	158.1148
	Squared method	186.3200	310.7543	15.26394	162.6485	674.9867

Judging from the simulated curves of these figures, we can conclude that the proposed controller works well. Furthermore, the tracking performance is better and the speed of response is faster than the other methods when there is a model mismatch.

Example 2

Let us consider another example of a hot oil fractionator with four inputs and two outputs (Ogunnaike and Ray, 1994)

G_{p} (s) = [\begin{matrix} \begin{matrix} \frac{4.05 e^{- 81 s}}{50 s + 1} & \frac{1.77 e^{- 84 s}}{60 s + 1} \end{matrix} & \begin{matrix} \frac{5.88 e^{- 81 s}}{50 s + 1} & \frac{1.44 e^{- 81 s}}{40 s + 1} \end{matrix} \\ \begin{matrix} \frac{4.38 e^{- 60 s}}{33 s + 1} & \frac{4.42 e^{- 66 s}}{44 s + 1} \end{matrix} & \begin{matrix} \frac{7.2}{19 s + 1} & \frac{1.26}{32 s + 1} \end{matrix} \end{matrix}]

Step 1. Calculate the ETF array ${\hat{G}}_{m} (s)$ . Similarly to Example 1, the ETF matrix is given as

\begin{matrix} {\hat{G}}_{m}^{T} (s) = {[\begin{matrix} \begin{matrix} \frac{60.5909 s + 1}{4.6764} e^{60 s} & \frac{- 8.8095 s + 1}{- 2.1881} e^{63 s} \end{matrix} & \begin{matrix} \frac{- 18.8415 s + 1}{7.9532} & \frac{21.6829 s + 1}{7.0746} \end{matrix} \\ \begin{matrix} \frac{12.7674 s + 1}{- 8.7532} e^{60 s} & \frac{- 1.0061 s + 1}{2.5399} e^{66 s} \end{matrix} & \begin{matrix} \frac{- 197.3751 s + 1}{- 58.8664} & \frac{10.1295 s + 1}{- 10.7225} \end{matrix} \end{matrix}]}^{T} \end{matrix}

Step 2. Design the decoupling internal model controller $G_{IMC} (s)$ . In view of the robustness and response time, filters are selected as

F (s) = diag {\frac{1}{145 s + 1}, \frac{1}{95 s + 1}}

Then, the final decoupling controller is obtained as

G_{IMC} (s) = [\begin{matrix} \begin{matrix} \frac{60.5909 s + 1}{4.6764 (145 s + 1)} e^{- 3 s} \\ \frac{- 8.8095 s + 1}{- 2.1881 (145 s + 1)} \end{matrix} & \begin{matrix} \frac{12.7674 s + 1}{- 8.7532 (95 s + 1)} e^{- 6 s} \\ \frac{- 1.0061 s + 1}{2.5399 (95 s + 1)} \end{matrix} \\ \begin{matrix} \frac{- 18.8415 s + 1}{7.9532 (145 s + 1)} e^{- 63 s} \\ \frac{21.6829 s + 1}{7.0746 (145 s + 1)} e^{- 63 s} \end{matrix} & \begin{matrix} \frac{- 197.3751 s + 1}{- 58.8664 (95 s + 1)} e^{- 66 s} \\ \frac{10.1295 s + 1}{- 10.7225 (95 s + 1)} e^{- 66 s} \end{matrix} \end{matrix}]

Step 3. Design the filter $G_{F} (s)$ . Feedback loop filters are considered

G_{F} (s) = diag {\frac{1}{42 s + 1}, \frac{1}{32 s + 1}}

To show the improvement of our method, the centralized PI controller $G_{PID} (s)$ is also designed for the plant in the example. By taking the pseudo-inverse of the steady-state gain matrix, Davison’s method (Sarma and Chidambaram, 2003) is used to design the controller $G_{PID} (s)$ . The controller form of $G_{PID} (s)$ is given by:

G_{PID} (s) = K_{C} + \frac{K_{I}}{s}

where $K_{C} = δ [G_{p}^{*} (s = 0)]$ and $K_{I} = ε [G_{p}^{*} (s = 0)]$ are the proportional and integral gain matrices, respectively; $[G_{p}^{*} (s = 0)]$ is called the rough tuning matrix, and $δ$ and $ε$ are the fine tuning parameters, which generally range from 0 to 1.

The steady-state gain matrix is obtained as

G_{p} (s = 0) = [\begin{matrix} 4.05 & 1.77 & 5.88 & 1.44 \\ 4.38 & 4.42 & 7.2 & 1.26 \end{matrix}]

The corresponding pseudo-inverse is obtained as:

G_{p}^{*} (s = 0) = {[\begin{matrix} 0.2138 & - 0.457 & 0.1257 & 0.1414 \\ - 0.1142 & 0.3937 & - 0.017 & - 0.0933 \end{matrix}]}^{T}

Thus, the controller transfer function matrix can be written as

\begin{matrix} G_{PID} (s) = K_{C} + \frac{K_{I}}{s} \\ = [\begin{matrix} \begin{matrix} 0.2138 δ + \frac{0.2138 ε}{s} \\ - 0.457 δ - \frac{0.457 ε}{s} \end{matrix} & \begin{matrix} - 0.1142 δ - \frac{0.1142 ε}{s} \\ 0.3937 δ + \frac{0.3937 ε}{s} \end{matrix} \\ \begin{matrix} 0.1257 δ + \frac{0.1257 ε}{s} \\ 0.1414 δ + \frac{0.1414 ε}{s} \end{matrix} & \begin{matrix} - 0.017 δ - \frac{0.017 ε}{s} \\ - 0.0933 δ - \frac{0.0933 ε}{s} \end{matrix} \end{matrix}] \end{matrix}

To obtain good closed-loop performances, the tuning parameters are selected as $δ = 0.37$ and $ε = 0.01$ . At a perfect match of a process plant and model, the set-point responses of the system for a two-step change in $r_{1}$ and $r_{2}$ are shown in Figures 9 and 10, while the disturbance responses are shown in Figures 11 and 12, where distractions are introduced the system for a 0.5 and −0.5 step change at 1000 s in loop 1 and loop 2, respectively.

G_{p} (s) = [\begin{matrix} \begin{matrix} \frac{4.86 e^{- 97.2 s}}{40 s + 1} & \frac{2.124 e^{- 100.8 s}}{48 s + 1} \end{matrix} & \begin{matrix} \frac{7.056 e^{- 97.2 s}}{40 s + 1} & \frac{1.728 e^{- 97.2 s}}{32 s + 1} \end{matrix} \\ \begin{matrix} \frac{5.256 e^{- 72 s}}{26.4 s + 1} & \frac{5.304 e^{- 79.2 s}}{35.2 s + 1} \end{matrix} & \begin{matrix} \frac{8.64}{15.2 s + 1} & \frac{1.512}{25.6 s + 1} \end{matrix} \end{matrix}]

When the model is mismatched, the set-point responses are shown in Figure 13. It can be seen from the simulation results that the method proposed in this paper has a smaller overshoot, faster tracking features, less interaction, better anti-interference performance, and is driven to steady-state values under the given perturbation.

From the simulation curves of these figures, we can conclude that the proposed controller has better tracking performance and faster response speed than the PID control method and squared method whether the model is matched or not.

As can be seen from the above two examples, the proposed IMC based on ETF not only has a good performance on tracking given value and overcoming disturbance, but also has good decoupling and control quality both in model match and mismatch cases. Moreover, the designed controller is easy to compute and realize. In brief, the designed controller works well, and can achieve a certain control level.

Conclusion

A new method for decoupling IMC, focusing on a class of non-square stable processes with multiple time delays, has been presented in this paper. The method in this paper proceeded from the desired diagonal transfer function matrix of the system, based on the stability and realizability, to design a compensator, and then obtained the final decoupling internal model controller based on the pseudo-inverse. Furthermore, the calculation of the pseudo-inverse based on ETF matrix was very simple, straightforward and easy to understand, and the ETF element has the same form as the open-loop transfer function element, so the decoupling internal model controller has a simple form. The results of simulation have illustrated that our approach gives better control performance and robust performance than other methods. However, there are several limitations of our present study, which should be taken into consideration. as possible limitations of the study, it should be noted that all our results were obtained under the assumption of perfect control, which barely exists in practical engineering. Without question, for the design of multi-variable control systems, there often exists conservation by using the proposed method. Another limitation is that we did not extend this method to non-square systems with fewer inputs than outputs, while the experiments were directed only toward another kind of non-square system with more inputs than outputs, probably because it is difficult to achieve good results in cases where the number of control variables are fewer than output variables. An additional research issue that should be tackled is how to design IMC controllers for non-square systems with fewer inputs than outputs, which we are currently planning.

Footnotes

Acknowledgements

The authors are grateful to the anonymous reviewers for their valuable recommendations.

Funding

This work was partially supported by the financial support of the National Science Foundation of China (grant number 61273132) and the National Grand Fundamental Research 973 Program of China (grant number 2007CB714300) and partly by the Automation Institute Beijing University of Chemical Technology.

References

Chang

(1990) The relative gain for non-square multivariable systems. Chemical Engineering Science 45(5): 1309–1323.

Chen

(2011) A new control method for MIMO first order time delay non-square systems. Journal of Process Control 21(4): 538–546.

Chen

Sun

. (2008) Modified internal model control and its application in non-square processes. Control and Decision 23(5): 581–584.

Garcia

Morari

(1982) Internal model control. 1. A unifying review and some new result. Industrial & Engineering Chemistry Process Design and Development 21(2): 308–323.

Garcia

Morari

(1985) Internal model control. 2. Design procedure for multivariable systems. Industrial & Engineering Chemistry Process Design and Development 24(2): 472–484.

Jin

Guo

Liu

. (2010) Decoupling internal model control for non-square process with time delays. International Conference on Measuring Technology and Mechatronics Automation 1: 898–901.

Lin

Jang

Shieh

. (1999) Generalized multivariable dynamic artificial neural network modeling for chemical processes. Industrial & Engineering Chemistry Research 38: 4700–4711.

Liu

Zhang

(2006) Analytical design of decoupling internal model control (IMC) scheme for two-input–two-output (TITO) processes with time delays. Industrial & Engineering Chemistry Research 45(9): 3149–3160.

Loh

Chiu

(1997) Robust decentralized control of non-square systems. Chemical Engineering Communications 158(1): 157–180.

10.

Morari

Grimm

Oglesby

. (1985) Design of resilient processing plants-VII: Design of energy measurement systems for unstable plants-new insights. Chemical Engineering Science 40(2): 187–198.

11.

Morari

Zafiriou

(1989) Robust Process Control. Englewood Cliffs, NJ: Prentice Hall.

12.

Normey-Rico

Bordon

Camacho

(1997) Improving the robustness of dead-time compensating PI controllers. Control Engineering Practice 5(6): 801–810.

13.

Ogunnaike

Ray

(1979) Multivariable controller design for linear systems having multiple time delays. AIChE Journal 25(6): 1043–1057.

14.

Ogunnaike

Ray

(1994) Process Dynamics, Modeling and Control. New York: Oxford University Press.

15.

Sarma

KLN

Chidambaram

(2003) Design of centralized controllers for nonsquare systems. In: Proceedings of the National Conference on Instrumentation and Control, pp. 5–6.

16.

Sarma

KLN

Chidambaram

(2005) Centralized PI/PID controllers for non-square systems with RHP zeros. Indian Institute Science 85(4): 201–214.

17.

Shen

Cai

(2010) Normalized decoupling control for high-dimensional MIMO processes for application in room temperature control HVAC systems. Control Engineering Practice 18(6): 652–664.

18.

Treiber

Hoffman

(1986) Multivariable constraint control using a frequency domain design approach. In: Proceeding of the 3rd Conference on Chemistry Process Control, pp. 185–230.

19.

Xiong

Cai

(2007) Equivalent transfer function method for PI/PID controller design of MIMO processes. Journal of Process Control 17(8): 665–673.

20.

Zhang

(2007) Pseudo-inverse multivariate/matrix-variate distributions. Journal of Multivariate Analysis 98(8): 1684–1692.