A comparative study on existing and new methods to design internal model controllers for non-square systems

Abstract

This paper focuses on a challenging problem in the internal model control (IMC) strategy: the model inversion to design the IMC controller for non-square systems. Several existing approaches for the synthesis of a specific inversion of the identified model will be presented in this paper to deal with the differences between the system’s inputs and outputs numbers. The non-square effective relative gain is firstly presented. It consists of the measurement of interactions between the loops of the system in order to square the system and make it invertible. The equivalent transfer function method is presented as well. It is based on tuning the pseudo-inverse of the process to design the internal model controller. These methods are then compared with a novel proposed model inversion approach based on virtual outputs method. Virtual adding is considered in order to obtain an invertible square transfer matrix to design the internal model controller. This simple yet effective method ensures robust control performance. Its efficiency and availability, as compared with other presented methods, is illustrated through simulations on an overactuated system with three inputs/two outputs.

Keywords

Internal model control model inversion non-square systems non-square effective relative gain equivalent transfer function virtual outputs overactuated system

Introduction

Processes having multiple inputs and multiple outputs (MIMO) occur frequently in the modern complex industrial fields. They can be partitioned into square and non-square systems (Moro and Bonavitat, 1997). Processes with equal number of inputs and outputs are called square systems. Non-square systems are then referred to unequal number of inputs and outputs processes (Ben Omrane et al., 2012; Jin et al., 2014; Qibing et al., 2009; Quan et al., 2012). The control of such systems is challenging due to the existence of interactions among input and output variables.

Several researches have been done in the control field of non-square systems. A simple control method consists of squaring non-square systems by adding or eliminating the appropriate number of the input or output variables. However, this control method may lead to an increase of the control cost when adding variables or reduced control performance when removing inputs or outputs, due to neglected dynamics (Chen et al, 2011; Rao and Baram, 2006; Reeves and Arkum, 1989).

This paper deals with the synthesis of the internal model control (IMC) designed for non-square systems. This control strategy is considered as an excellent vigor to settle the control problems of MIMO systems due to its simple design, strong robustness, and its good tracking ability performance (Ben Omrane et al., 2012; Qibing et al., 2009). The IMC incorporates a model of the system. The Tuning of the IMC controller corresponds to a specific inversion of an appropriate identified model, since direct inversion is rarely achievable. Considering non-square systems makes the IMC controller design much more complicated, as it cannot get the traditional sense of inversion (Fossard, 1972; Dhahri et al., 2016; Qibing et al., 2012).

Some available results in literature consider several model inversion methods to design the IMC controller for non-square systems. The pseudo inverse method can be considered to inverse the transfer matrix of the non-square model. However, the application of this method may be difficult and very complicated.

The non-square effective gain (NERGA) method is as well proposed in Quan et al. (2012) to inverse the transfer matrix of the model. This method is considered as a criterion to choose a square subsystem from a non-square system by adding or removing appropriate inputs or outputs. Once the system is squared, controller synthesis procedures based on the IMC strategy, developed for square systems, are applied directly. In Pavarati and Subbulekshmi (2011) and Luan et al. (2014), the equivalent transfer function (ETF) method were proposed to design the IMC controller by approximating the pseudo-inverse of the non-square process in a simple and easy way. In this paper, a novel inversion approach based on virtual outputs (VO) method is proposed. It consists on attaching lines to the non-square transfer matrix of the model in order to square it. This virtual adding aims to obtain a square transfer matrix and have no influence on the response of the system (Fossard, 1972).

In the present work, the control performances of non-square effective relative gain, ETF and VO methods for inversion the model process to design the internal model controller for non-square systems are evaluated and compared.

This paper is organized as follows. Section 2 describes the internal model controller design for non-square systems and explains the proposed model inversion methods. Section 3 shows simulation results with concluding remarks in section 4.

IMC design

In this section, we present the IMC. On the one hand, we describe its basic principle and its properties. On the other hand, we present the structure of the IMC controller as well as the study of its stability. Finally, solutions for the problem of inversion are presented. They aim to design a controller that is easy to carry out and close as possible to the model inverse (Yaman et al., 2016).

The IMC structure

The block diagram of IMC structure for a multivariable system is exposed in Figure 1 (Chu et al. 2018).

Where: G(s) is the transfer matrix of the process with “m” inputs and “n” outputs; C_IMC(s) is the internal model controller designed to control the process G(s); M(s) is the process.

Considering the case of perfect modeling, so we have G(s) =M(s). From Figure 1, the following equation can be derived for nominal systems

(\begin{matrix} y \\ u \end{matrix}) = (\begin{matrix} G (s) & [I_{m} - G (s) C_{IMC} (s)] G (s) \\ C_{IMC} (s) & - C_{IMC} (s) G (s) \end{matrix}) (\begin{matrix} r \\ d \end{matrix})

(1)

Where I_m is the identify matrix

The closed-loop system shown in figure 1 is internally stable if and only if all four elements of the matrix in (1) are stable.

The internal model controller design

The IMC controller designed for multivariable systems presenting or not multiple time delays and / or with a non-minimum phase shift is presented by Figure 2 (Dhahri et al., 2016).

In order to ensure good control performance and tracking feature, the IMC controller must satisfy the following equation

C_{IMC} (s) M (s) = I_{m}

(2)

From Figure 2, the closed-loop transfer matrix C_IMC(s) between e and u is derived as follows (Dhahri et al., 2016)

C_{IMC} (s) = K_{2} (I_{m} + K_{1} M (s))^{- 1} K_{1}

(3)

where K₁ is the inversion matrix of the form $K_{1} = a I_{m}, a \in R^{+}$

In steady state, the gain matrix K₂ is written in the following form

K_{2} = (I_{m} + K_{1} M (0)) \times (K_{1} M (0))^{- 1}

(4)

where M(0) is the of static gains matrix of the model M(s).

The gain matrix K₂ aims to compensate the static errors of the system, while K₁ is chosen in order to ensure the stability of the controller and to achieve the inverse of the model M(s).

On the basics of the IMC principle, the synthesis of a controller equal to the inverse of the model transfer matrix ensures perfect tracking and good control performances. Nevertheless, the inversion in a direct way is practically impossible, especially for non-square systems. To overcome this problem, several methods of non-square model inversions are proposed and developed in the next sections.

Model inversion using NERGA method

This method is based on the analysis and pairing the loops in order to overcome the singularity of system transfer matrix (Xiong and Jian, 2006). The non-square effective relative gain array “NERGA” is recognized as a spread of relative gain array “RGA”. The RGA has proved its ability to choose the input-output pairing for diagonal control design by illustrating the intensity of coupling within the different process loops and exhibiting the impact degree of each loop by other circuit (Shinde et al., 2015). The NERGA corresponds to loop pairing criterion by employing the steady-state gains and bandwidth of the non-square process transfer function elements.

Consider a non-square system with “n” rows and “m” columns, presented by the following process transfer function matrix G(s)

G (s) = (\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 m} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 m} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{n 1} (s) & g_{n 2} (s) & \dots & g_{nm} (s) \end{matrix})

(5)

The g_ij(s) presents the process transfer function element linking the i^th input and j^th output.

The effective gain matrix EG(s) (Quan et al., 2012) is given by the following expression

EG (s) = [e_{ij} (s)] = G (0) \otimes Ω

(6)

where: G(0) the static gain matrix and Ω the bandwidth matrix are expressed as follows (Quan et al., 2012; Xiong and Jian, 2006)

G (0) = (\begin{matrix} g_{11} (0) & g_{12} (0) & \dots & g_{1 m} (0) \\ g_{21} (0) & g_{22} (0) & \dots & g_{2 m} (0) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{n 1} (0) & g_{n 2} (0) & \dots & g_{nm} (0) \end{matrix})

(7)

and

Ω = (\begin{matrix} ω_{11} & ω_{12} & \dots & ω_{1 m} \\ ω_{21} & ω_{22} & \dots & ω_{2 m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ω_{n 1} & ω_{n 2} & \dots & ω_{nm} \end{matrix})

(8)

where $ω_{ij}$ presents the bandwidths.

We report that the operator ⊗ is the Hadamard product. It consists on applying an element by element multiplication.

The effective relative gain ψ_ij(s) between the output y_i and the input u_j is given by the ratio ψ_ij = e_ij/ ê_ij. The term ê_ij indicates the effective gain between y_i and u_j when other loops are closed (Chang and Yu, 1990; Xiong and Jian, 2006).

Once the gains of all the existing combinations between the various inputs and outputs of the system are calculated, we can obtain the NERGA matrix defined by formula (9)

ψ = (\begin{matrix} ψ_{11} & ψ_{12} & \dots & ψ_{1 m} \\ ψ_{21} & ψ_{22} & \dots & ψ_{2 m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ψ_{n 1} & ψ_{n 2} & \dots & ψ_{nm} \end{matrix}) = EG \otimes {(E G^{+})}^{T}

(9)

where EG and $E G^{+}$ are, respectively, the effective gain matrix and its pseudo inverse.

E G^{+} = {[EG]}^{H} . {[EG . E G^{H}]}^{- 1}

(10)

Throw the non-square effective relative gain matrix, the sum of the elements in each row and column, respectively noted by RS and CS, can be calculated and expressed as follows

\begin{matrix} RS = {[\begin{matrix} \sum_{j = 1}^{n} ψ 1 j & \sum_{j = 1}^{n} ψ 2 j & \dots & \sum_{j = 1}^{n} ψ nj \end{matrix}]}^{T} \\ CS = {[\begin{matrix} \sum_{i = 1}^{n} ψ i 1 & \sum_{i = 1}^{n} ψ i 2 & \dots & \sum_{i = 1}^{n} ψ in \end{matrix}]}^{T} \end{matrix}

(11)

For overactuated systems (m>n), the sum of elements in each row RS is equal to 1, and the sum of elements in each column CS is between 0 and 1. If CS(i) is adjacent or near to null, the i^th input loop admits a small impact on the output loop. Thus, the corresponding loop can be easily overlooked and consequently deleted. In the other case, if the process is under-actuated (m<n), the sum of elements in each column CS is equal to 1, and the sum of elements in each row RS is located in the interval [0 1]. If RS(i) is near to zero, the i^th input loop admits a little impact on the output loop. Consequently, this loop can be neglected and removed from the system (Quan et al., 2012; Schaper et al., 1999; Xiong et al, 2005).

Model inversion using ETF

The equivalent transfer function (ETF) method has been developed to approximate the inverse of the process transfer matrix based on the measurement of dynamic interactions of non-square systems (Jin and Jiang, 2013; Jin et al., 2014).

The majority of the industrial systems are processes of second order plus dead time (SOPDT). The SOPDT transfer function is described below

g_{ij} (s) = \frac{k_{ij}}{a_{ij} s^{2} + b_{ij} s + 1} e^{- θ_{ij} s}

(12)

The ETF is given by equation (13); it takes into account loop interactions

{\hat{g}}_{ij} (s) = \frac{{\hat{k}}_{ij}}{a_{ij} s^{2} + b_{ij} s + 1} e^{- {\hat{θ}}_{ij} s}

(13)

The expression of the pseudo-inverse of non-square process $G_{m}^{*} (s)$ obtained by the ETF method is expressed as follows (Luan et al., 2014, 2015)

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

(14)

where

\hat{G} (s) = (\begin{matrix} \frac{1}{{\hat{g}}_{11} (s)} & \frac{1}{{\hat{g}}_{12} (s)} & \dots & \frac{1}{{\hat{g}}_{1 m} (s)} \\ \frac{1}{{\hat{g}}_{21} (s)} & \frac{1}{{\hat{g}}_{22} (s)} & \dots & \frac{1}{{\hat{g}}_{2 m} (s)} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{1}{{\hat{g}}_{n 1} (s)} & \frac{1}{{\hat{g}}_{n 2} (s)} & \dots & \frac{1}{{\hat{g}}_{nm} (s)} \end{matrix})

(15)

The basic steps required to obtain ETF for a non-square system, are given by the following algorithm.

Algorithm ETF

❖ Step 1: Calculate the open-loop gain array (K^N) and the time delay of the system.

\begin{matrix} K^{N} = {[k_{ij}]}_{n \times m} = {[g_{ij} (j 0)]}_{n \times m} \\ i = 1, 2 \dots, n, j = 1, 2, \dots, m \end{matrix}

(16)

Θ = [\begin{matrix} θ_{11} & θ_{12} & \dots & θ_{1 m} \\ θ_{21} & θ_{22} & \dots & θ_{2 m} \\ \begin{matrix} ⋮ \end{matrix} & ⋮ & ⋱ & ⋮ \\ θ_{n 1} & θ_{n 2} & \dots & θ_{nm} \end{matrix}]

(17)

❖ Step 2: Calculate the non-square relative gain array (NRGA) (Jin et al., 2014; Shen et al., 2010).

The NRGA between output y_i and input variable u_j is defined as the ratio of two gains representing, $k_{ij}$ the open-loop gain when all loops are open and ${\hat{k}}_{ij}$ the closed-loop gain when every loops other than the loop i–j are closed (Liao et al., 2018)

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

(18)

The NRGA in matrix form is defined as follows

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

(19)

where ${(K^{N})}^{* T}$ is the transposed of the pseudo inverse of G(s) and ⊗ is the Hadamard product.

❖ Step 3: Calculate the non-square relative normalized gain array (NRNGA) (Jin et al, 2014; Xiong et al., 2007).

The NRNGA can be calculated as described below

\begin{matrix} Φ^{N} & = [\begin{matrix} ϕ_{11} & ϕ_{12} & \dots & ϕ_{1 m} \\ ϕ_{21} & ϕ_{22} & \dots & ϕ_{2 m} \\ \begin{matrix} ⋮ \end{matrix} & ⋮ & ⋱ & ⋮ \\ ϕ_{n 1} & ϕ_{n 2} & \dots & ϕ_{nm} \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 m}^{N}}{{\hat{k}}_{N, 1 m}^{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 m}^{N}}{{\hat{k}}_{N, 2 m}^{N}} \\ \begin{matrix} ⋮ \end{matrix} & ⋮ & ⋱ & ⋮ \\ \frac{k_{N, n 1}^{N}}{{\hat{k}}_{N, n 1}^{N}} & \frac{k_{N, n 2}^{N}}{{\hat{k}}_{N, n 2}^{N}} & \dots & \frac{k_{N, nm}^{N}}{{\hat{k}}_{N, nm}^{N}} \end{matrix}] \\ = K_{N}^{N} \otimes {(K_{N}^{N})}^{* T} \end{matrix}

(20)

where $k_{N, ij}^{N}$ is the normalized gain describing the dynamic properties of a particular transfer function g_ij(s). It is defined as follows (Jin et al., 2014; Luan et al., 2017; Shen et al., 2010)

k_{N, ij}^{N} \overset{Δ}{=} \frac{k_{ij}}{σ_{ij}} = \frac{k_{ij}}{τ_{ij} + θ_{ij}}, i = 1, 2 \dots, n, j = 1, 2, \dots m

(21)

$σ_{ij}$ is the average residence time (ART) . It indicates the response speed of the controlled variable y_i to the manipulated variable u_j. The ART can be obtained by the following expression

σ_{ij} = \frac{\int_{0}^{\infty} | y (\infty) - y (t) | dt}{y (\infty)}

(22)

where $y (t)$ and $y (\infty)$ are, respectively, the open loop step response and the steady state response value.

The normalized gain matrix $K_{N}^{N}$ (He et al., 2009; Ravitaja et al., 2016) is expressed as follows by extending equation (21)

K_{N}^{N} = k_{ij} ⊙ T_{ar}^{N}

(23)

where $T_{ar}^{N} = {[τ_{ij} + θ_{ij}]}_{n \times m}$ and ⊙ indicates element by element division of the corresponding elements.

❖ Step 4: Calculate the non-square relative average residence time array (NRARTA) (Jin et al., 2014; Shen et al., 2010).

The NRARTA can be calculated by

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

(24)

where

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

(25)

❖ Step 5: Calculate the ETF parameters (Jin et al., 2014; Xiong et al., 2007).

The gain of the ETF matrix is defined as follows

\hat{K} = K^{N} ⊙ Λ^{N}

(26)

The time delays of the ETF matrix are expressed as described below

\hat{Θ} = Γ^{N} \otimes Θ

(27)

The obtained ETF matrix serves to evaluate the internal model controller.

The IMC controller (Jin et al., 2014) is calculated as follows:

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

(28)

where

F (s) = diag {\frac{1}{{(λ_{i} s + 1)}^{n_{i}}}, i = 1, 2 \dots n}

(29)

F(s) is a low pass filter. It ensures the realization of the IMC controller. The parameter $λ_{1}$ is the filter time constant. It determines the speed of response. n_i is the minimum order for the filter to guarantee a strictly proper internal model controller.

$G_{m +} (s)$ is the process model part that contains the time delays and the right half-plane zeros.

The mathematical description of $G_{m +} (s)$ is given by equation (30) (Wang et al., 2002)

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

(30)

where $θ_{k_{i}}$ is the maximum advance term in the i-th column of pseudo inverse $G_{m}^{*} (s)$ , $z_{p}$ is the right half-plane poles in corresponding column, $z_{p}^{*}$ is the conjugate complex of $z_{p}$ , W_i is the multiplicity of the pole and D_i is the number of different poles in the i-th column of $G_{m}^{*} (s)$ .

A filter $G_{f} (s)$ is a well needed in order to improve control performance, as settling time response, and the robustness of the IMC towards model uncertainties (Normey et al., 1997).

The first order filter $G_{f} (s)$ is described below

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

(31)

where $T_{f} = \frac{θ_{k_{i}}}{2}$ , is the filter time constant.

Model inversion using VO method

Over-actuated systems are multivariable systems with more control inputs than outputs. As a result, the transfer matrix of an over-actuated linear system is not square. The model M(s) must be chosen very close to G(s) the process transfer matrix. The controller C_IMC(s) is chosen as an approximate inverse of the model. The idea consists on squaring the model to be inverted and then eliminating excess outputs that will be released to the system.

A new inversion technique of the model M(s) by means of the VO is proposed. It consists on attaching lines to the non-square transfer matrix of the system G(s). This virtual adding aims to obtain a square transfer matrix, which allows the inversion operation. The excess outputs are then eliminated from the system.

To maintain the original performances of the studied system, we have to remove the added variables through the programmation part. The VO are only used to square the model transfer matrix (Foassard, 1972).

Once the VO method is applied, we can now deal with a non-singular matrix. This method is more effective for non-square model inversion as compared with the pseudo inversion that can give several solutions and the controller remains unsuitable.

For overactuated systems, the VO amount is referred to a transfer matrix of dimension ((m-n), m).

The studied overactuated system can be shown through equation (32)

(\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{n} \end{matrix}) = (\begin{matrix} g_{11} (s) & g_{12} (s) & \dots & g_{1 n} (s) \\ g_{21} (s) & g_{22} (s) & \dots & g_{2 m} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{n 1} (s) & g_{m 2} (s) & \dots & g_{nm} (s) \end{matrix}) (\begin{matrix} u_{1} \\ u_{2} \\ ⋮ \\ u_{m} \end{matrix})

(32)

In order to simplify the study and to avoid inversion problems, the new ((m-n), m) transfer matrix that contains the added variables transfer matrix can be chosen as first-order transfer functions or constants satisfying the stability criterion. It is presented as follows

(\begin{matrix} y_{n + 1} \\ y_{n + 2} \\ ⋮ \\ y_{m} \end{matrix}) = (\begin{matrix} g_{n + 1, 1} (s) & g_{n + 1, 2} (s) & \dots & g_{n + 1, m} (s) \\ {g_{n +}}_{2, 1} (s) & {g_{n +}}_{2, 2} (s) & \dots & g_{n + 2, m} (s) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ g_{m 1} (s) & g_{m 2} (s) & \dots & g_{mm} (s) \end{matrix}) (\begin{matrix} u_{1} \\ u_{2} \\ ⋮ \\ u_{m} \end{matrix})

(33)

The proposed IMC structure proposed in Touati et al. (2013) for square multivariable systems is then modified in order to eliminate the excess (m-n) outputs. A block called (elimination of (n-m) outputs) is inserted into the basic IMC structure (cf. Figure 1). The aim of this block is to eliminate the (m-n) excess outputs by means of the usual logical operators, as depicted in Figure 3.

The new IMC structure of multivariable over-actuated system is then shown in Figure 4.

Study of the stability of the controller

The IMC controller C_IMC(s) expressed by equation (3) can be reformulated using the state space representation as follows

{\begin{matrix} x_{c} (t) = A_{c} x (t) + B_{c} e (t) \\ y_{c} (t) = C_{c} x (t) \end{matrix}

(34)

where: $x (t) \in R^{n} and y (t) \in R^{m}$ are, respectively, the state vector and the system output vector. The system (A_c, B_c, C_c) is assumed to be controllable.

The stability analysis can be processed by the knowledge of the state matrix A_c, which depends of the gain matrix K₁.

In this paper, the Lyaponov method is considered to study the controller stability (Maamri et al., 2006).

The stability condition highlighted by Lyapunov stability theory requires that the system (34) is stable if and only if there exists a positive definite matrix $P = P^{T} > 0$ , satisfying the following Lyapunov inequality (Maamri et al., 2006)

P > 0, P A_{c} + P^{T} A_{c} < 0

(35)

This work is restricted to the use of the Linear matrix inequality technique to study the stability of the internal model controller and to calculate the range of the gain K1 that ensures the stability of the open loop system (Adelipour et al., 2014).

The resolution of an LMI problem with Matlab version 2012a is presented according to the following algorithm:

First, we initialize the controller matrix A_c.

Second, we start to describe and initialize the LMI description.This is easily done by the function Setlmis([]) , this function does not have any parameter.

Then we indicate,with the function Lmivar , the decision variables, what are the problem unknowns and what type they are. For our case, the unknowns are the matrix P which is symmetric and the gain matrix K₁

In the next step, we indicate the different terms that make the LMI depending on whether they are constant or that they answer for P. We define our LMI : $P A_{c} + P^{T} A_{c} < 0$ , by function lmiterm() .

The final step, before launching the numerical resolution, consists on indicating that the definition of the problem LMI is finished, and that the function is feasable by the corresponding solver Feasp .

Solver for LMI feasibility problems is L(x) < R(x). This solver minimizes t subject to L(x) < R(x) + t*I. where L(.) and R(.) are affine functions of some structured matrix variables x₁,…, x_n. A simple example is the Lyapunov inequality defined by equation (35).

The best value of time t should be negative for feasibility simulation results. The equation (34) is stable in the sense of Lyapunov.

Simulation results

In order to underline the interest and the performance of the proposed approaches, the air path of a turbocharged diesel engine (Lamara et al., 2012) is considered. This overactuated system is 4 cylinders with 2 liters and 80.9 KW. The three pneumatic air path actuators are wastegate (WG), EGR valve (EGR) and intake throttle (TH). The air path system of the turbocharged Diesel engine is described by the Figure 5 (Lamara et al., 2012).

The dynamics representation of the air path of a turbocharged diesel engine, with three inputs (EGR, TH and WG) and two outputs (Q_Air and P_Col), is given by the following transfer matrix (Lamara et al., 2012)

\begin{matrix} (\begin{matrix} Q_{air} \\ P_{col} \end{matrix}) & = (\begin{matrix} \frac{- 2771.69 e^{- 0.033 s}}{(s + 45) (s + 19)} & \frac{- 1331.31 e^{- 0.026 s}}{(s + 14) (s + 60)} & \frac{188.6 e^{- 0.44 s}}{(s + 17) (s + 7)} \\ \frac{1928.22 e^{- 0.003 s}}{(s + 64) (s + 15.1)} & \frac{- 1397.03 e^{- 0.11 s}}{(s + 42.8) (s + 14.8)} & \frac{100.23 e^{- 0.55 s}}{(s + 20) (s + 1)} \end{matrix}) \\ (\begin{matrix} EGR \\ TH \\ WG \end{matrix}) \end{matrix}

(36)

The eigen values of G(s) $(λ_{1} \dots λ_{12})$ have negative real parts, so the open loop system is stable

\begin{matrix} λ_{1} = - 64, λ_{2} = - 60, λ_{3} = - 45, λ_{4} = - 42.80, \\ λ_{5} = - 20, λ_{6} = - 19, λ_{7} = - 17, λ_{8} = - 15.1, \\ λ_{9} = - 14.79, λ_{10} = - 14, λ_{11} = - 7, λ_{12} = - 1 \end{matrix}

The simulation study is considered for a four-cylinder diesel engine speed of 3000 rpm, according to the parameters described in Table 1.

Table 1.

Diesel engine operational parameters.

Torque (Nm)	$[50 100]$	$[100 150]$
Time (s)	$t_{1} \in [14 22]$	$t_{2} \in [22 30]$
Q_air (g/s)	42	60
P_col (mbar)	1460	1820

Internal model controller design based on NERGA

Starting with the application of the non-square effective relative gain array (NERGA) approach. The (2×3) system is squared down to a (2×2) system as we can delete either an input or an output that has a little effect on our process.

To obtain the effective gain matrix of our system, we should firstly calculate both static gain matrix G(0) and the bandwidth matrix Ω presented by the following equations

G (0) = (\begin{matrix} - 3.2417 & - 1.5849 & 1.5849 \\ 1.9953 & - 2.2055 & 5.0115 \end{matrix})

(37)

Ω = (\begin{matrix} 16.5398 & 13.2954 & 6.1262 \\ 14.3256 & 13.3811 & 0.9952 \end{matrix})

(38)

Using the Hadamard product, the effective gain matrix EG is equal to

EG = (\begin{matrix} - 53.6170 & - 21.0719 & 9.7095 \\ 28.5839 & - 29.5121 & 4.9872 \end{matrix})

(39)

Throw the formula (9); we obtain our NERGA expressed by:

Ψ = (\begin{matrix} 0.7060 & 0.2541 & 0.0398 \\ 0.2875 & 0.6877 & 0.0248 \end{matrix})

(40)

Once the NERGA is evaluated, we can now calculate the sum of the elements in each row and column. $RS = {[\begin{matrix} 1 & 1 \end{matrix}]}^{T}$ and $CS = {[\begin{matrix} 0.9935 & 0.9418 & 0.0646 \end{matrix}]}^{T}$ . From the results shown previously, we notice that the third term in CS is the closest term to zero among the other two terms. Thus, the third input loop of the system can be neglected and consequently removed.

As the third input loop is now deleted, the new square transfer matrix describing the system is expressed as shown below

(\begin{matrix} Q_{air} \\ P_{col} \end{matrix}) = (\begin{matrix} \frac{- 2771.69 e^{- 0.033 s}}{(s + 45) (s + 19)} & \frac{- 1331.31 e^{- 0.026 s}}{(s + 14) (s + 60)} \\ \frac{1928.22 e^{- 0.003 s}}{(s + 64) (s + 15.1)} & \frac{- 1397.03 e^{- 0.11 s}}{(s + 42.8) (s + 14.8)} \end{matrix}) (\begin{matrix} EGR \\ TH \end{matrix})

(41)

Thus, the new non-square effective relative gain array is calculated as

Ψ = (\begin{matrix} 0.7243 & 0.2757 \\ 0.2757 & 0.7243 \end{matrix})

(42)

We can obviously claim that the coupling problem between the different input and output loops has disappeared.

Consider the case of imperfect modeling. So, $G (s) \neq M (s)$ . M(s) is described by the following transfer matrix (43)

M (s) = (\begin{matrix} \frac{- 2000 e^{- 0.03 s}}{(s + 30) (s + 10)} & \frac{- 10001 e^{- 0.02 s}}{(s + 11) (s + 50)} \\ \frac{1500 e^{- 0.004 s}}{(s + 60) (s + 14)} & \frac{- 1400 e^{- 0.05 s}}{(s + 30) (s + 13)} \end{matrix})

(43)

The eigen values of M(s) $(λ_{1} \dots λ_{8})$ have negative real parts, so the model is stable

\begin{matrix} λ_{1} = - 60, λ_{2} = - 50, λ_{3} = λ_{4} = - 30, λ_{5} = - 14, \\ λ_{6} = - 13, λ_{7} = - 11, λ_{8} = - 10 \end{matrix}

The application of the LMI approach allows us to deal with quadratic stability of the internal model controller $C_{IMC} (s)$ that is detailed by equation (3). LMI has been performed in MATLAB environment. By solving the LMI in equation (35), we get the following equation

P = [\begin{matrix} 0.0351 & - 0.0496 & 0.0013 & - 0.0161 & - 0.0046 & 0.0325 & - 0.0006 & 0.0115 \\ - 0.0496 & 0.5149 & 0.0198 & 0.0667 & - 0.0166 & 0.0387 & - 0.0080 & - 0.0134 \\ 0.0013 & 0.0198 & 0.0269 & - 0.0390 & 0.0003 & 0.0189 & 0.0038 & 0.0186 \\ - 0.0161 & 0.0667 & - 0.0390 & 0.1743 & - 0.0141 & - 0.0166 & - 0.0209 & - 0.0502 \\ - 0.0046 & - 0.0166 & 0.0003 & - 0.0141 & 0.0273 & - 0.0394 & - 0.0001 & 0.0124 \\ 0.0325 & 0.0387 & 0.0189 & - 0.0166 & - 0.0394 & 0.2471 & - 0.0186 & 0.0620 \\ - 0.0006 & - 0.0080 & 0.0038 & - 0.0209 & - 0.0001 & - 0.0186 & 0.0558 & - 0.0873 \\ 0.0115 & - 0.0134 & 0.0186 & - 0.0502 & 0.0124 & 0.0620 & - 0.0873 & 0.3637 \end{matrix}]

(44)

The gain interval K₁ ensuring the stability of the controller is given by

0 ≺ K_{1} ≺ 0.20 \times I_{2}

(45)

For the studied system, two chosen gain matrices K₁ are considered: $K_{1} = 0.1 \times I_{2}$ and $K_{1} = 0.001 \times I_{2}$ . Their corresponding gain matrices K₂ are represented, respectively, as follows

\begin{matrix} K_{2} = {(\begin{matrix} - 1.13 & 1.53 \\ - 1.93 & - 2.14 \end{matrix})}_{K_{1} = 0.1} and \\ K_{2} = {(\begin{matrix} - 212.93 & 153.69 \\ - 193.51 & - 313.37 \end{matrix})}_{K_{1} = 0.001} \end{matrix}

(46)

Figures 6 and 7 display, respectively, the system control actions (EGR and TH)) and system outputs (Q_air and P_col ) for $K_{1} = 0.1 \times I_{2}$ and $K_{1} = 0.001 \times I_{2}$ .

It can be seen that the proposed control low ensures an accurate system. The system outputs Q_air and P_col have good tracking features. However, the output Q_air presents a relatively important overshoot due to the negative gains of the model. The control action EGR presents a small overshoot for $K_{1} = 0.1 \times I_{2}$ which can be explained by the effect of the choice of a. Increasing a ( $K_{1} = a \times I_{2}$ ) gives best system performance. However, control action peaks are more significant.

Internal model controller design based on ETF

This section is dedicated to the ETF approach. First of all, we will start with the resolution of the ETF algorithm to design the internal model controller. Then, we will apply this algorithm for the studied system to test its robustness and performance.

The steady-state gain of the system G(s) is given as follows

K^{N} = [\begin{matrix} - 3.2417 & - 1.5849 & 1.5849 \\ 1.9953 & - 2.2055 & 5.0115 \end{matrix}]

(47)

The time delay of the system is given below

Θ = [\begin{matrix} 0.033 & 0.026 & 0.44 \\ 0.003 & 0.11 & 0.55 \end{matrix}]

(48)

The average residence time array is obtained as

T = [\begin{matrix} 0.0749 & 0.0881 & 0.2017 \\ 0.0819 & 0.0909 & 1.05 \end{matrix}]

(49)

Thus, the non-square RGA is

\land = [\begin{matrix} 0.7737 & 0.1351 & 0.0912 \\ 0.1869 & 0.1157 & 0.6974 \end{matrix}]

(50)

The non-square RNGA is calculated as follows

Φ^{N} = [\begin{matrix} 0.5082 & 0.4734 & 0.0185 \\ 0.4912 & 0.4757 & 0.0331 \end{matrix}]

(51)

The non-square relative average residence time array is obtained as

Γ^{N} = [\begin{matrix} 0.6568 & 3.5034 & 0.2026 \\ 2.6279 & 4.1107 & 0.0475 \end{matrix}]

(52)

According to the previous results the ETF parameters are

\hat{K} = [\begin{matrix} - 4.1897 & - 11.7302 & 17.3875 \\ 10.6742 & - 19.0597 & 7.1864 \end{matrix}]

(53)

\hat{Θ} = [\begin{matrix} 0.0217 & 0.0911 & 0.0892 \\ 0.0079 & 0.4522 & 0.0261 \end{matrix}]

(54)

Using the above results and equation (15), the ETF matrix is calculated below

\hat{G} (s) = (\begin{matrix} \frac{(0.0012 s^{2} + 0.0749 s + 1) e^{0.0217 s}}{- 4.1897} & \frac{(0.001 s^{2} + 0.0819 s + 1) e^{0.0079}}{10.6742} \\ \frac{(0.0012 s^{2} + 0.0881 s + 1) e^{0.0911 s}}{- 11.7302} & \frac{(0.0016 s^{2} + 0.0909 s + 1) e^{0.4522}}{- 19.0597} \\ \frac{(0.0084 s^{2} + 0.2017 s + 1) e^{- 0.089}}{17.3875} & \frac{(0.05 s^{2} + 1.05 s + 1) e^{0.0261}}{7.1864} \end{matrix})

(55)

$G_{m +} (s)$ can be obtained using equation (30)

G_{m +} (s) = diag {e^{- 0.0911 s}, e^{- 0.4522 s}}

(56)

The filter time constants can be chosen as: $λ_{1} = 0.2, λ_{2} = 0.5$ and the order of filter n_i = 2.

F (s) = diag {\frac{1}{{(0.2 s + 1)}^{2}}, \frac{1}{{(0.5 s + 1)}^{2}}}

(57)

Hence, the feedback filter is

G_{F} (s) = diag {\frac{1}{0.0165 s + 1}, \frac{1}{0.055 s + 1}}

(58)

Finally, a stable IMC controller can be obtained from equation (28). It is represented by the following transfer matrix

\begin{matrix} C_{IMC} (s) = \\ (\begin{matrix} \frac{(0.0012 s^{2} + 0.0749 s + 1) e^{- 0.0694}}{- 4.1897 {(0.2 s + 1)}^{2}} & \frac{(0.001 s^{2} + 0.0819 s + 1) e^{- 0.4443}}{10.6742 {(0.5 s + 1)}^{2}} \\ \frac{(0.0012 s^{2} + 0.0881 s + 1)}{- 11.7302 {(0.2 s + 1)}^{2}} & \frac{(0.0016 s^{2} + 0.0909 s + 1)}{- 19.0597 {(0.5 s + 1)}^{2}} \\ \frac{(0.0084 s^{2} + 0.2017 s + 1) e^{- 0.0019}}{17.3875 {(0.2 s + 1)}^{2}} & \frac{(0.05 s^{2} + 1.05 s + 1) e^{- 0.4261}}{- 4.1897 {(0.5 s + 1)}^{2}} \end{matrix}) \end{matrix}

(59)

Figure 8 displays the system inputs and outputs evolution. Set-point tracking is established for both outputs. However, important overshoot in Q_air response is noticed. It can lead to serious problems for real time implementation, despite the fact that it is gradually diminished due to the feedback loop. The P_co response presents as well a slower tracking feature as compared with the results obtained for NERGA method.

So, the application of the ETF-IMC method has reduced control performance as compared with the NERGA method.

Internal model controller design based on VO

Dealing now with the new VO method proposed in this paper, and applied to the same system studied previously.

The model M(s) is chosen close to the process G(s). For our case the process G(s) is of dimension (2×3). The VO method considered for this overactuated systems, consists on adding ((3-2), 3)) vector to the original transfer system in order to square it.

The model M(s) is then chosen of dimension (3×3) so that it can be invertible. The elements of the transfer matrix M(s) are the same as G(s) and we will add a third line to the matrix G(s). The elements of this line are transfer functions with fast responses and do not influence the process G(s).

The chosen transfer function matrix of the model is given by the following matrix

M (s) = (\begin{matrix} \frac{- 2771.69 e^{- 0.033 s}}{(s + 45) (s + 19)} & \frac{- 1331.31 e^{- 0.026 s}}{(s + 14) (s + 60)} & \frac{188.6 e^{- 0.44 s}}{(s + 17) (s + 7)} \\ \frac{1928.22 e^{- 0.003 s}}{(s + 64) (s + 15.1)} & \frac{- 1397.03 e^{- 0.11 s}}{(s + 42.8) (s + 14.8)} & \frac{100.23 e^{- 0.55 s}}{(s + 20) (s + 1)} \\ \frac{1 e^{- 2 s}}{0.1 s + 1} & \frac{2 e^{- s}}{0.2 s + 1} & \frac{1 e^{- 0.5 s}}{0.5 s + 1} \end{matrix})

(60)

The eigen values of M(s) $λ_{1} \dots λ_{15}$ have negative real parts, so M(s) is stable

\begin{matrix} λ_{1} = - 64, λ_{2} = - 60, λ_{3} = - 45, \\ λ_{4} = - 42.80, λ_{5} = - 20, λ_{6} = - 19, λ_{7} = - 17, \\ λ_{8} = - 15.1, λ_{9} = - 14.79 \\ λ_{10} = - 14, λ_{11} = - 10, λ_{12} = - 7, λ_{13} = - 5, \\ λ_{14} = - 2, λ_{15} = - 1 \end{matrix} .

Solving the LMI in equation (34) for M(s) given by the equation (60), we obtain a matrix P of dimension $(15 \times 15)$ , this matrix ensures the stability of the controller and the range of the gain K₁ which ensures the stability of the system is

0 ≺ K_{1} ≺ 0.21 \times I_{3}

(61)

\begin{matrix} K_{2} = {(\begin{matrix} - 1.32 & 0.39 & 0.55 \\ - 1.13 & 2.82 & - 2.44 \\ 0.42 & 2.63 & - 0.29 \end{matrix})}_{K_{1} = 0.1} \\ and K_{2} = {(\begin{matrix} - 212.93 & 153.69 \\ - 193.51 & - 313.37 \end{matrix})}_{K_{1} = 0.001} \end{matrix}

(62)

Figures 9 and 10 display, respectively, the system control actions (EGR, TH and WG)) and system outputs (Q_air and Pcol ) for $K_{1} = 0.1 \times I_{2}$ and $K_{1} = 0.001 \times I_{2}$ .

The system outputs present good accuracy properties, smaller overshoot and have faster tracking features for $K_{1} = 0.1 \times I_{2}$ and $K_{1} = 0.001 \times I_{2}$ as compared with the NERGA-IMC and ETF-IMC approaches.

Comparison analysis

Once the three inversion methods: NERGA, ETF and the VO, are evaluated, a comparative study seems to be interesting to demonstrate and to compare their efficiencies.

In this section, the air circuit is subject to set point changes. Q_Air is regulated at 42 g / s and 60 g / s, respectively, during the two time ranges $t_{1} = [14 22]$ and $t_{2} = [22 30]$ . As for P_Col, it is regulated at 1460mbar and 1820 mbar during the same time ranges t₁ and t₂. The control objective of the CMI is to ensure fast and short response time to set point variations and robust performance.

The response time of 5%, evaluated for the three inversion techniques, gives a good evaluation of the speed of a system; it expresses the time taken by the system subjected to a step to reach its steady state value within ± 5%. We note that the NERGA-CMI structure and VO-CMI have better performance in terms of response time versus ETF-CMI, as shown in Table 2.

Table 2.

Quantitative comparison between IMC Controller based on NERGA, ETF and VO method.

Time range (s)	Setting time (s)				Rise time (s)
	Output Q_air		Output P_col		Output Q_air		Output P_col
	t₁	t₂	t₁	t₂	t₁	t₂	t₁	t₂
NERGA-IMC	0.30	0.32	0.19	0.21	0.37	0.34	0.22	0.22
ETF-IMC	0.99	1.00	1.94	2.18	0.85	0.85	1.81	1.85
VO-IMC	0.58	0.61	0.20	0.23	0.52	0.66	0.24	0.26

The performance index for the transient regime is evaluated by calculating the rise time of the air circuit, which represents the time taken by the Q_Air and P_Col outputs to go from 10% of their final values to 90% of their final values. The NERGA-CMI method performs better than other approaches in terms of climb time, as shown in Table 2.

However, the NERGA-CMI structure is not always applicable, especially when the sum of elements at each row and column are close. So, we cannot delete an input that does not affect the system to make it square. This operation may lead to performance degradation due to the elimination of certain information.

We note also that the ETF-CMI airflow outputs converge to set points with slow response and rise times as compared with the virtual output and NERGA-CMI method. This method requires a long calculation time that complicates the implementation in real time. In addition, the response of the Q_air output presents some relatively high peaks, as shown in Figure 7. Besides, the ETF-CMI method presents a problem in the choice of the filter, which requires an adequate knowledge of the G (s).

The proposed VO approach overcomes these difficulties. It has good rise time and response, an accurate response, does not require the use of filters and the elimination of system inputs. This new simple yet effective method leads to better control system performances.

A comparison between the IMC control based on VO method and the decentralized crone control (Oustaloup, 1991) can be proposed. It aims to show the effectiveness the proposed IMC strategy as compared with crone control.

The decentralized crone control strategy is based on fractional-order differentiation (Oustaloup, 1991). The aim of the CRONE strategy for MIMO systems is to ensure robustness of the closed-loop matrix dynamic by maintaining the damping or resonance factors (Oustaloup, 1991). The crone control can manage the robustness/performance tradeoff.

For the proposed system, the decentralized crone controller K(s) is described by the transfer matrix below

K = [\begin{matrix} \frac{- 14.3 (s^{3} + 11.5 s^{2} + 48.3 s + 1.5)}{s^{5} + 25.9 s^{4} + 229.6 s^{3} + 796.2 s^{2} + 484.6 s} & 0 \\ \frac{- 7.9 (s^{3} + 12.2 s^{2} + 54.2 s + 80.6)}{s^{5} + 27.2 s^{4} + 268.2 s^{3} + 1023 s^{2} + 635.9 s} & 0 \\ 0 & \frac{- 1.2 (s^{3} + 7.1 s^{2} + 15.1 s + 9.1)}{s^{4} + 15.4 s^{3} + 61.2 s^{2} + 51.67 s} \end{matrix}]

Figure 11 illustrates simulation results for this scenario.

The comparison between the IMC strategy and the decentralized crone control shows the benefit of the IMC strategy. It ensures better dynamic performance, smaller overshoots, and smaller settling time.

Robustness of the VO-IMC approach towards external disturbances

In order to test the robustness of the IMC controller based on VO method, it is proposed to disrupt the system by applying a white noise of amplitude 3 for Q_air and P_col between 18s and 25s. Figure 12 illustrates simulation results for this scenario.

It can be seen that the system remains stable and the disturbance is completely rejected. The IMC controller based on VO inversion approach is able to compensate these disturbances. This fact confirms its robustness.

Robustness of the VO-IMC approach towards parameter uncertainties

The test of the controller robustness allows us to check whether the IMC control based on the VO method is capable of compensating the uncertainties on the air path system parameters.

Figure 13 displays the evolution of the system outputs for an uncertain transfer matrix $G + Δ G$ where G is the system transfer function described by equation (32) and $Δ G = \pm 50 % G$ .

We notice that despite uncertainties, set-point tracking is established with acceptable performance of accuracy, speed, and stability. This fact, confirms the IMC–VO strategy robustness towards uncertainties. It should be noted that when increasing uncertainty amount, the system becomes slower and the overshoot is greater.

The IMC control based on virtual outputs approach can be easily implemented in real time and the gains K₁ and K₂ remain the same despite uncertainties.

Conclusion

In this paper, we have been interested in the comparative study of several inversion methods for IMC controller design proposed for non-square systems.

The proposed process inversion methods are NERGA, ETF and a new proposed approach, namely virtual outputs method. They were applied on an overactuated system with two outputs and three inputs. A comparative analysis was proposed. It shows clearly the efficiency of the virtual outputs method and its ability to avoid the complicated calculation, such the matrix inverse. The obtained controller is simple, has few tuning parameters, is easy to implement in industrial processes and presents good control and robustness performances.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflict of interests with respect to the research, authorship and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Imen Saidi

Nahla Touati

References

Adelipour

Abooee

Haeri

(2014) LMI-based sufficient conditions for robust stability and stabilization of LTI-fractional-order systems subjected to interval and polytopic uncertainties. Transactions of the Institute of Measurement and Control 37(10): 1207–1216

Ben Omrane

Chatti

Borne

(2012) Evolutionary method for designing and learning control structure of a wheelchair. Studies in Informatics and Control 21(2): 155–164.

Chang

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

Chen

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

Chu

(2018) Graphic IMC-PID tuning based on maximum sensitivity for uncertain systems. Transactions of the Institute of Measurement and Control 41(8): 2196–2204.

Dhahri

Saidi

Soudani

(2016) A New internal model control method for MIMO over-actuated systems. International Journal of Advanced Computer Science and Applications 7(10): 304–309.

Fossard

(1972) Commande des systèmes multidimensionnels. Paris: Dunod.

Cai

Xie

(2009) RNGA based control system configuration for multivariable processes. Journal of Process Control 19(6): 1036–1042.

Jin

Jiang

(2013) Pseudo-feed-forward decoupling internal model control for singularity plants with time delays. Transactions of the Institute of Measurement and Control: 36(4) 565–576.

10.

Jin

Jiang

Wang

Shan

(2014) Decoupling internal model control for non-square processes based on equivalent transfer function. Transactions of the institute of Measurement and Control 36(8): 1114–1131.

11.

Lamara

Colin

Lanusse

, et al. (2012) Decentralized robust control-system for a non-square MIMO system, the air-path of a turbocharged Diesel engine. IFAC Proceedings 45(30) : 130–137.

12.

Liao

Q F

Sun

(2018) Interaction measures for control configuration selection based on interval type-2 Takagi–Sugeno fuzzy model. IEEE Transactions on Fuzzy Systems 26(5): 2510–2523.

13.

Luan

Chen

Liu

(2017) Interaction measurement for complex multivariable models with various reference inputs based on RNGA. In: 11th Asian Control Conference, Gold Coast Convention Centre, Australia, 17–20 December 2017, PaperCept, Inc, US pp. 2411–2416.

14.

Luan

Chen

Liu

(2014) Centralized PI control for high dimensional multivariable systems based on equivalent transfer function. ISA Transactions 53(5): 1554–1561.

15.

Luan

Chen

Albertos

Liu

(2015) Stabilizing parametric region of multiloop PID controllers for multivariable systems based on equivalent transfer function. Mathematical Problems in Engineering: 1–7 doi:10.1155/2016/3173289.

16.

Maamri

Trigeassou

J C

Bentayeb

Bachelier

(2006) PID controllers design for MIMO plants with time specifications LMI-Approach. IFAC Proceedings Volumes (IFAC-PapersOnline) 5: 178–184.

17.

Moro

Bonavitat

(1997) IDCOM-B multivariable model predictive controller. Transactions of the Institute of Measurement and Control 19(4): 192–201.

18.

Normey-Rico

Bordon

Camacho

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

19.

Oustaloup

(1991) La commande CRONE. Paris: Hermes Science Publications.

20.

Qibing

Rong

Ghoshal

T K

(2012) A V-norm decoupling IMC method with filters based on inverted decoupling. Applied Mechanics and Materials 197: 311–315.

21.

Quan

Zhang

Yang

(2012) A novel internal model control method for non-square system. Applied Mechanics and Materials 238: 66–70.

22.

Rao

Chidambaram

(2006) Smith delay compensator for multivariable non-square systems with multiple time delays. Computers and Chemical Engineering 30(8): 1243–1255.

23.

Raviteja

Rao

(2016) Improved controller design for two-input-two-output (TITO) unstable processes. Resource-efficient Technologies 2: 76–86.

24.

Reeves

Arkun

(1989) Interaction measures for non-square decentralized control structures. Journal American Institute of Chemical Engineers 35(4): 603–613.

25.

Schaper

Kailath

Lee

(1999) Decentralized control of water temperature for multizone rapid thermal processing systems. IEEE Transactions on Semiconductor Manufacturing 12(2): 193–199.

26.

Shen

Cai

W J

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

27.

Shinde

Hamde

Waghmare

(2015) Predictive PI control multivariable non-square system with multiple time delays. International Journal for Science and Research in Technology 1(5): 26–29.

28.

Touati

Soudani

Naceur

Benrejeb

(2013) On an internal multimodel control for nonlinear multivariable systems – a comparative study. International Journal of Advanced Computer Science and Applications 4(7): 66–71.

29.

Wang

Zhang

Chiu

(2002) Decoupling internal model control for multivariable systems with multiple time delays. Chemical Engineering Science 57(1): 115–124.

30.

Xiong

Cai

(2005) A practical loop pairing Criterion for multivariable processes Journal Process Control 15(7): 741–747.

31.

Xiong

Cai

W J

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

32.

Xiong

Jian

(2006) Effective transfer function method for decentralized control system Design of Multi-Input Multi-Output Process. Journal of Process Control 16(8): 773–784.

33.

Yaman

Ustoglu

Ucun

(2016) Internal model controller design via HN dynamic output feedback. Transactions of the Institute of Measurement and Control 39(3): 260–260.