Disturbance rejection for multivariable processes with time delays based on partial inverted decoupling control

Abstract

Motivated by the thought of using couplings to reject disturbances, the authors propose a novel partial inverted decoupling control for improving the disturbance rejection of multivariable processes with time delays. With the relative disturbance gain (RDG) for n × n processes, whether the couplings among loops are favorable is evaluated. The inverted decoupling usually used in full decoupling control is then extended to partial decoupling control to retain the favorable couplings and remove the unfavorable ones. The elements of the partial inverted decoupler are derived and the implementation structure is provided. Based on the partially decoupled processes, the decentralized controllers for both the decoupled and non-decoupled loops are designed. An analytical design method for the decentralized controllers is proposed to improve disturbance rejection. Controller parameters are determined in terms of both performance and robustness. Two simulation examples show that the presented partial inverted control strategy can provide superior disturbance rejection to several other classical multivariable control strategies.

Introduction

In process control industries, most controlled processes are multivariable processes in nature. An important characteristic of multivariable processes is the couplings between the inputs and outputs. ¹ Couplings make the control problem for multivariable processes much more complicated than that for univariate processes. Therefore, how to handle couplings is a key issue in multivariable control.

Decentralized control, with no doubt is the most commonly used control strategy for multivariable processes in practice due to its easy implementation. ² The primary step for decentralized control is usually to determine the main channel of each loop. Therefore, how to pair the inputs and outputs so as to reduce relative couplings is of great importance ³ and many indicators have been proposed, such as relative gain array (RGA) ⁴ and some relevant methods with process dynamics considered.^5,6 With the well-paired system, the decentralized controllers can be tuned by detuning method, ⁷ auto-tuning method ⁸ and iterative design method, ⁹ etc. Recently, Mahapatro et al. ¹⁰ designed decentralized PI controllers based on frequency response matching of a predefined reference model and the actual model. Mahapatro and Subudhi ¹¹ proposed a robust decentralized PID controller based on complementary sensitivity function through a graphical tuning method. For minimizing the undesirable effects of the couplings, Euzébio et al. ² optimized decentralized controllers by determining the cost function as the sum of the integral absolute errors (IAE) from each output due to set-point changes in coupled loops. Jeng and Lee ¹² proposed a direct synthesis method for decentralized controllers on the basis of specifying the frequency response of the desired closed-loop transfer function matrix. Apparently, the decentralized controllers designed only with the main channels may not achieve sufficiently satisfactory control performances since the coupling information is ignored. As a result, the effective open-loop transfer function (EOTF) was proposed for independent design of multi-loop controllers. ¹³ In the formulation of the EOTF, the coupling channels were contained. With the concept of EOTF, the multivariable control system can be decomposed into several equivalent single-loops for decentralized controller design.^14,15 However, since the calculation of EOTF is complex, especially for high-dimensional processes, model approximation is usually inevitable to simplify controller design. Approximation errors will affect the control performance. In sum, decentralized control cannot fully remove the effect of couplings in spite of its simple structure, which means it is suitable for the process with modest couplings.

When significant couplings exist, an intuitional way to improve control performance is to remove or reduce the loop couplings by inserting a decoupler. When the process is well decoupled, the control strategies for single-input-single-output (SISO) processes can be used directly. Hence, the design of the decoupler is the key issue in decoupling control. Generally, decoupling strategies are classified into two categories, that is, static decoupling and dynamic decoupling. ¹⁶ A static decoupler is designed by only considering the steady state of a process. It has the advantages of simple calculation and easy implementation while it may not provide satisfactory decoupling performance because of the neglected process dynamics. A dynamic decoupler can achieve better decoupling performance but it needs more accurate process information. As the most representative multivariable processes, many dynamic decouplers are firstly studied for 2 × 2 processes. Among them, the most classical three methods are ideal decoupling, simplified decoupling and inverted decoupling. ¹⁷ Ideal decoupling which uses the inverse of the process to design the decoupler is rarely used in practice because it suffers from complex decouplers. Compared with ideal decoupling, simplified decoupling is more popular as it provides a simple and implementable decoupler. However, a simplified decoupler usually results in a complex decoupled process, which is not desirable for controller design. This will be more evident when simplified decoupling is extended to high-dimensional processes. ¹⁸ Unlike the above two, inverted decoupling can achieve simple decouplers and decoupled processes simultaneously, but the limitations of achieving stable and realizable decoupler elements should be considered when it is applied for high-dimensional processes. ¹⁹ To overcome the effect of the mismatched model, Aguiar et al.^20,21 proposed closed-loop evaluation and redesign methodologies for simplified decoupling and inverted decoupling. In essence, these dynamic decoupling methods seek to make the decoupled process matrix a diagonal matrix at more frequencies. Motivated by this, Jin and Liu ²² use the Nyquist set to describe the dynamics of the decoupler and use a complex-curve fitting technique to obtain an approximate decoupler matrix. Moreover, Xiong et al. ²³ proposed the equivalent transfer function (ETF) based on the inherent relation between the inverse process and the dynamic relative gain array (DRGA) to simplify the design of the decoupler. Naik et al. ²⁴ approximated the ETF using the static relative gains and the average residence time. Khandelwal and Detroja ²⁵ presented a novel ETF parameterization procedure for large scale processes or systems with higher-order elements. To simplify the structure of multivariable system, the multivariable controller can have the abilities of both decoupling and control, which is called the centralized controller. Usually, a centralized controllers can be realized by combining a decoupler and a set of decentralized controllers. ²⁶ Besides, a centralized controllers can also be synthesized directly based on a desired closed-loop transfer function matrix. ²⁷ Garrido et al. ²⁸ presented an iterative methodology for designing centralized PID controllers based on the proposed equivalent loop transfer functions. Gündeş ²⁹ gave the necessary and sufficient conditions for the existence of decoupling controllers. Actually, the couplings can be regarded as the internal disturbances of the multivariable systems, which means the methods for handling couplings is to suppress the internal disturbances.

Apart from couplings, external disturbances and process uncertainties are inevitable in many industries, which makes the primary goal of the controllers to mitigate the effects of external disturbances and process uncertainties.^30,31 Feedback control is an important strategy for suppressing disturbances. As the classical feedback controller, PID controller still dominates the engineering systems. The key issue of PID controller is to tune the controller parameters. The traditional tuning formula, such as Z-N method, are experiential methods, which are easy to implement while only achieve limited control performance. When the process model is known, the synthesis of PID controller can be transformed to solving the optimization problem with given objectives and constraints.^32,33 With the development of machine learning techniques, some are incorporated to PID tuning to improve control performance such as neural network PID, ³⁴ fuzzy PID ³⁵ and swarm intelligent PID ³⁶ , etc. To reduce the dependence of accurate models, data-driven techniques, like iterative feedback tuning ³⁷ and virtual reference feedback tuning, ³⁸ are employed to tune fix-structure controller, which is suitable for the processes with big uncertainties. Apart from feedback control, feedforward control is another effective scheme for disturbance rejection, which cancels disturbances with the reverses of disturbance signals. In light of the main drawback of feedforward control that it can only deal with measurable disturbances, many disturbance observers are designed to estimate the disturbances by using the measurable states and the process dynamics, which results in the disturbance observer-based (DOB) control. DOB control contains several related methods, such as frequency domain DOB control, ³⁹ active disturbance rejection control (ADRC), ⁴⁰ equivalent input disturbance (EID), ⁴¹ etc. Although under different names, these methods share the same basic idea, that is, the lumped disturbances including the external disturbances and uncertainties are estimated and corresponding compensations are generated based on the estimations. Due to realizing the compensation of uncertainties, DOB control does not require an accurate process model. ³¹ Meanwhile, its inherent 2 degree-of-freedom (DOF) structure allow the DOB control system to achieve a satisfactory performance in both set-point tracking and disturbance rejection. As a result, DOB control has been applied into more and more industries.^42–44

As the typical characteristic of the multivariable system, the couplings can also be used for disturbance rejection when they have the opposite effect of the external disturbances. Motivated by this, Stanley et al. ⁴⁵ defined the relative disturbance gain (RDG) to evaluate couplings for distillation processes. Similarly, Lin et al. ⁴⁶ proposed the relative load gain (RLD) as the ratio of the gains of effective disturbances and open-loop disturbances. Moreover, Chang and Yu ⁴⁷ systematically described the controller structure via internal model control (IMC) structure and proposed the relative disturbance gain array (RDGA) to synthesize the controller structure straightforwardly. Based on the evaluation of couplings, the partial decoupling structure can be constructed and some control strategies have been studied.^48,49 As well, the thought of using couplings for disturbance rejection has been extended to other control structures.^50,51

Motivated by the thought of using couplings to reject disturbances, this work is aimed to propose a novel partial decoupling control strategy. Unlike the existing methods which cannot remove the unfavorable couplings thoroughly due to the approximation or simplification of the decoupler, the proposed strategy extends the inverted decoupler which used in full decoupling control to partial decoupling control. The designed partial inverted decoupler can realize better elimination of the unfavorable couplings. Besides, in the existing methods, the controllers of the non-decoupled loops are usually designed based on the main channel models and disturbance channel models. The dynamic characteristics of non-decoupled loops are not analyzed sufficiently. In the proposed method, the equivalent process model and the equivalent disturbance channel model are obtained for designing the controllers of the non-decoupled loops. Finally, to simplify the design of the controllers, an analytical design method for the decentralized controllers is proposed by considering the performance of disturbance rejection and robustness. In sum, the main contributions of this work are concluded as follows:

• (1) The RDG for general n × n processes is explored to evaluate whether the couplings are favorable or not.

• (2) The inverted decoupling is extended to partial decoupling control to retain the favorable couplings and remove the unfavorable ones. The elements of the partial inverted decoupler are derived and the implementation structure is provided.

• (3) The equivalent process models and the equivalent disturbance channel models are obtained to design the controllers of the non-decoupled loops.

• (4) An analytical design method for the decentralized controllers is proposed with the consideration of both the performance of disturbance rejection and robustness.

This paper is organized as follows. In Section “Relative disturbance gain (RDG),” the RDG for general n × n processes is derived to evaluate the couplings. Based on the evaluation result, the design and implementation of the partial inverted decoupler are discussed in Section “Partial inverted decoupling control structure.” In Section “Design of the decentralized controller,” the decentralized controllers for decoupled and non-decoupled loops are designed respectively. An analytical design method is proposed. Two simulations examples are provided in Section “Simulation study.” Section “Conclusion” is the conclusion.

Relative disturbance gain (RDG)

With the purpose of using the couplings between loops to reject disturbances, the first but important step is to judge whether the couplings are favorable. Hence, an index for evaluating the couplings is necessary. With this index, the couplings beneficial to disturbance rejection should be retained while the unfavorable ones have to be decoupled. Herein, the index is chosen as RDG.

The RDG was proposed by Stanley et al. ⁴⁵ for 2 × 2 processes, like dual-composition control. Equation (1) shows the definition of the RDG of the ith loop, that is, β _i , where u _i and y _j are the input of the ith loop and the output of the jth loop, respectively; d is a disturbance signal.

$${\beta }_i=\frac{{\left(\frac{\partial u_i}{\partial d}\right)}_{\mathrm{perfect}\mathrm{control}\mathrm{over}{y}_j(j\ne i)}}{{\left(\frac{\partial u_i}{\partial d}\right)}_{\mathrm{no}\mathrm{control}\mathrm{over}{y}_j(j\ne i)}}$$

(1)

Obviously, β _i is the ratio of the required input changes for rejecting the disturbances in two different cases. The numerator indicates the case where the decentralized control is used for the multivariable system while y _j (j ≠ i) is under perfect control. In this case, the couplings are retained. The denominator indicates the case where the system has been well decoupled since y _j (j ≠ i) has no effect on y _i . Therefore, it can be concluded that the couplings are favorable to disturbance rejection if |β _i | < 1, while harmful to disturbance rejection if |β _i | > 1. According to equation (1), the expression of β _i for n × n processes can be derived.

Figure 1(a) shows the commonly used decentralized control structure for general n × n processes, where C(s) is the decentralized controller, expressed as C(s) = diag(c₁(s), c₂(s), …, c _n (s)); P(s) is the controlled process, whose element is described as p _ij (s), i, j = 1, 2, …, n; L(s) is the disturbance transfer function vector, expressed as L(s) = [L₁(s), L₂(s), …, L _n (s)]^T;d(s) is a specific disturbance signal; r, u and y are the reference vector, input vector and output vector, respectively. Assuming that the loops have been paired, according to Figure 1(a), the ith disturbance response, $y_i^d$ can be expressed by (the Laplace operator s is omitted)

$$\begin{array}{c}\multicolumn{1}{c}{y_i^d=(p_{\mathrm{ii}}-P_{i*}{H}_{-\mathrm{ii}}{\left(P_{-\mathrm{ii}}\right)}^{-1}{P}_{*i})u_i}\\ \multicolumn{1}{c}{+(L_i-P_{i*}{H}_{-\mathrm{ii}}{\left(P_{-\mathrm{ii}}\right)}^{-1}{L}_{-i})d}\end{array}$$

(2)

where

$$H_{-\mathrm{ii}}={(I+C_{-i}{P}_{-\mathrm{ii}})}^{-1}{C}_{-i}{P}_{-\mathrm{ii}}$$

(3)

Figure 1.

Control structures of multivariable systems: (a) decentralized control structure and (b) decoupling control structure.

$L_{-i}$ is the disturbance transfer function vector with L _i removed; $P_{i*}$ and $P_{*i}$ are the ith row and ith column of transfer function matrix P with p _ii removed, respectively; C _−i and P _−ii are controller matrix C and transfer function matrix P with the ith row and ith column removed, respectively; I is the unit matrix. Since other loops are assumed to be under perfect control condition, H _−ii satisfies H _−ii ≈I. Substituting H _−ii  = I into equation (2), the numerator of β _i can be derived as

$${\frac{\partial u_i}{\partial d}|}_{y_i^d=0}=-\frac{L_i-P_{i*}{\left(P_{-\mathrm{ii}}\right)}^{-1}{L}_{-i}}{p_{\mathrm{ii}}-P_{i*}{\left(P_{-\mathrm{ii}}\right)}^{-1}{P}_{*i}}$$

(4)

On the other hand, when a multivariable system is well-decoupled, $y_i^d$ can be expressed as $y_i^d=p_{\mathrm{ii}}{u}_i+L_{i}d$ . Hence,

$${\frac{\partial u_i}{\partial d}|}_{y_i^d=0}=-\frac{L_i}{p_{\mathrm{ii}}}$$

(5)

Based on equations (4) and (5), β _i at steady state can be expressed as

$${\beta }_i=\frac{L_i(0)p_{\mathrm{ii}}(0)-P_{i*}(0){(P_{-\mathrm{ii}}(0))}^{-1}{L}_{-i}(0)p_{\mathrm{ii}}(0)}{L_i(0)p_{\mathrm{ii}}(0)-P_{i*}(0){(P_{-\mathrm{ii}}(0))}^{-1}{P}_{*i}(0)L_i(0)}$$

(6)

Next, β _i will be used to evaluate the couplings, which is essential to determine the constitution of the partial inverted decoupler.

Partial inverted decoupling control structure

The basic decoupling control structure is shown in Figure 1(b). In Figure 1(b), the decoupler Q(s), between the controller and the process, is the main part to be designed. The partial inverted decoupling control shares the same structure with Figure 1(b). The primary step is to design the partial inverted decoupler.

Design of the partial inverted decoupler

In decoupling control, the decoupler can be synthesized directly in terms of the desired decoupled process,

$$Q(s)=P^{-1}(s)T(s)$$

(7)

where T(s) is the desired decoupled process. For partial decoupling, T(s) here should be determined based on RDG. For n × n processes, assuming that there are m (m ≤ n) loops with unfavorable couplings, the process transfer function matrix can be reformed as

$$P(s)=\left[\begin{array}{cc}{P}_{11}(s)\hfill & P_{12}(s)\hfill \\ \multicolumn{1}{c}{P_{21}(s)}& P_{22}(s)\hfill \end{array}\right]$$

(8)

where P₁₁(s), P₁₂(s), P₂₁(s), and P₂₂(s) are an m × m, m × (n−m), (n−m)  × m and (n−m)  × (n−m) matrix, respectively. Since the first m loops need decoupling, the decoupled process should have the following form,

$$T(s)=\left[\begin{array}{cc}{T}_{11}(s)\hfill & 0\hfill \\ \multicolumn{1}{c}{T_{21}(s)}& T_{22}(s)\hfill \end{array}\right]$$

(9)

where T₁₁(s) = diag{t₁(s), t₂(s), …, t _m (s)}. To facilitate analysis, the decoupler matrix Q(s) is written as a block matrix, that is,

$$Q(s)=\left[\begin{array}{cc}{Q}_{11}(s)\hfill & Q_{12}(s)\hfill \\ \multicolumn{1}{c}{Q_{21}(s)}& Q_{22}(s)\hfill \end{array}\right]$$

(10)

The design of Q(s) is transformed to determine Q₁₁(s), Q₁₂(s), Q₂₁(s), and Q₂₂(s). In view of T(s) = P(s)Q(s), the following equations can be obtained:

$$\{\begin{array}{c}{P}_{11}(s)Q_{11}(s)+P_{12}(s)Q_{21}(s)=T_{11}(s)\hfill \\ P_{11}(s)Q_{12}(s)+P_{12}(s)Q_{22}(s)=0\hfill \\ P_{21}(s)Q_{11}(s)+P_{22}(s)Q_{21}(s)=T_{21}(s)\hfill \\ P_{21}(s)Q_{12}(s)+P_{22}(s)Q_{22}(s)=T_{22}(s)\hfill \end{array}$$

(11)

Considering that there are (n−m) loops using decentralized control, let Q₂₁(s) = 0 and Q₂₂(s) = I to simplify the design of the decoupler. As a result, equation (11) is rewritten as

$$\{\begin{array}{c}{P}_{11}(s)Q_{11}(s)=T_{11}(s)\hfill \\ P_{11}(s)Q_{12}(s)+P_{12}(s)=0\hfill \\ P_{21}(s)Q_{11}(s)=T_{21}(s)\hfill \\ P_{21}(s)Q_{12}(s)+P_{22}(s)=T_{22}(s)\hfill \end{array}$$

(12)

where the remaining undetermined parts are Q₁₁(s) and Q₁₂(s). According to equation (12), it is obvious that $Q_{11}(s)=P_{11}^{-1}(s)T_{11}(s)$ , which is a problematic solve as it requires the inverse of a transfer function matrix. To deal with this problem, Q₁₁(s) is designed with the inverted decoupling which is implemented by the V canonical structure shown in Figure 2.

Figure 2.

Structure of inverted decoupling.

where Q₁₁(s) is realized by two parts, that is, V₁(s) and V₂(s). According to Figure 2, Q₁₁(s) can be expressed as

$$Q_{11}(s)=V_1(s){[I-V_2(s)V_1(s)]}^{-1}$$

(13)

Substituting $Q_{11}(s)=P_{11}^{-1}(s)T_{11}(s)$ into equation (13), V₁(s) and V₂(s) should satisfy

$$V_1^{-1}(s)-V_2(s)=T_{11}^{-1}(s)P_{11}(s)$$

(14)

Since T₁₁(s) is a diagonal matrix, it is straightforward to compute its inverse. Similarly, V₁(s) should be determined as a diagonal matrix, that is, V₁(s) =diag{δ₁(s), δ₂(s), …, δ _m (s)}. For further simplification, V₂(s) is determined as

$$V_2(s)=\left[\begin{array}{cccc}0\hfill & {\rho }_{12}(s)\hfill & \dots \hfill & {\rho }_{1m}(s)\hfill \\ \multicolumn{1}{c}{{\rho }_{12}(s)}& 0\hfill & \dots \hfill & {\rho }_{1m}(s)\hfill \\ \multicolumn{1}{c}{\dots }& \dots \hfill & \dots \hfill & \dots \hfill \\ \multicolumn{1}{c}{{\rho }_{m1}(s)}& {\rho }_{m2}(s)\hfill & \dots \hfill & 0\hfill \end{array}\right]$$

(15)

To satisfy equation (14), the elements of V₁(s) and V₂(s) should satisfy

$$\{\begin{array}{cc}{\delta }_i(s)=\frac{t_i(s)}{p_{ii}(s)}& i=1,2,\dots ,m\\ {\rho }_{ij}(s)=-\frac{p_{ij}(s)}{t_i(s)}& i,j=1,2,\dots ,m\end{array}$$

(16)

Equation (16) shows that the elements of V₁(s) and V₂(s) rely on the selection of the elements of T₁₁(s).

After determining Q₁₁(s), the design of Q₁₂(s) will be undertaken. In light of equation (12) and Q₁₁(s), Q₁₂(s) is derived as

$$Q_{12}(s)=V_1(s){[I-V_2(s)V_1(s)]}^{-1}{V}_3(s)$$

(17)

where $V_3(s)=-T_{11}^{-1}(s)P_{12}(s)$ and its element η _ij (s) satisfies

$${\eta }_{\mathrm{ij}}(s)=-\frac{p_{i(m+j)}(s)}{t_i(s)}\begin{array}{c}i=1,2,\dots ,m,j=1,2,\dots ,n-m\hfill \end{array}$$

(18)

Obviously, the elements of V₃(s) are also related to the elements of T₁₁(s). Therefore, how to determine T₁₁(s) needs a discussion.

The principle of determining T₁₁(s)

The above procedure reveals that the partial decoupler can be accomplished after T₁₁(s) is determined. However, considering the realizability of V₁(s), V₂(s), and V₃(s), there are some constraints for the elements of T₁₁(s). Concretely, the elements of V₁(s), V₂(s), and V₃(s) should be causal, proper and stable, which results in the constraints for the time delay, relative degree and right-half-plane (RHP) zeros of t _i (s).

Constraint for the time delay

Equations (16) and (18) show that each element of V₁(s), V₂(s), and V₃(s) is a ratio of two transfer functions. Hence, to make sure the elements of V₁(s), V₂(s), and V₃(s) are causal, the time delay of t _i (s) should satisfy

$$\theta (p_{\mathrm{ii}}(s))\le \theta (t_i(s))\le \underset{j\ne i}{min}\theta (p_{\mathrm{ij}}(s))\begin{array}{c}i=1,2,\dots ,m,j=1,2,\dots ,n\hfill \end{array}$$

(19)

where θ(x) is the time delay of a transfer function x(s).

Constraint for the relative degree

Similarly, the elements of V₁(s), V₂(s), and V₃(s) will be proper if they have nonnegative relative degrees. Hence, the relative degree of t _i (s) should satisfy

$$\begin{array}{cc}\multicolumn{1}{c}{\phi (p_{\mathrm{ii}}(s))}& \le \phi (t_i(s))\le \underset{j\ne i}{min}\phi (p_{\mathrm{ij}}(s))\begin{array}{c}i=1,2,\dots ,m,j=1,2,\dots ,n\hfill \end{array}\hfill \end{array}$$

(20)

where φ(x) is the relative degree of a transfer function x(s).

Constraint for the RHP zeros

Moreover, the elements of V₁(s), V₂(s), and V₃(s) should be stable. Therefore, the numerator and denominator in equations (16) and in (18) should have the same RHP zeros and the multiplicity of each zero should satisfy

$${\psi }^k(p_{\mathrm{ii}}(s))\le {\psi }^k(t_i(s))\le \underset{j\ne i}{min}{\psi }^k(p_{\mathrm{ij}}(s))\begin{array}{c}i=1,2,\dots ,m,j=1,2,\dots ,n\hfill \end{array}$$

(21)

where ψ ^k (x) is the multiplicity of the kth zero of a transfer function x(s).

From Inequalities (19) to (21), it can be seen that p _ii (s) should have the smallest values in terms of time delay, relative degree and multiplicity of RHP zero in the ith row. However, the original process matrix in equation (8) may not meet these requirements. If so, the process matrix P₁₁(s) can firstly be reformed by interchanging the rows to meet these requirements. Unfortunately, in some cases, reformation may make no difference for the above problem. An effective solution is to introduce a compensator expressed as equation (22) to modify the characteristic of the process.

$$N_i(s)=e^{-{\alpha }_{i}s}\frac{1}{{(\mu s+1)}^{{\epsilon }_i}}\underset{j=1}{\overset{k}{\Pi }}{\left(\frac{-s+z_j}{s+z_j}\right)}^{{\mathrm{\Gamma }}_i^j}$$

(22)

N _i (s) is the compensator inserted in the ith loop and it is composed of three parts, which are used to compensate the time delay, relative degree and multiplicity of RHP zero, respectively. In equation (22), μ is a positive real constant; z _j denotes the jth RHP zero; α _i , ε _i and ${\mathrm{\Gamma }}_i^j$ are undetermined compensator parameters. With the compensators, the compensated process transfer function matrix can be expressed as P _C (s) = P(s)N(s), where N(s) is the compensator matrix, that is, N(s) = diag{N₁(s), N₂(s), …, N _n (s)}. The compensator parameters α _i , ε _i and ${\mathrm{\Gamma }}_i^j$ should make sure the compensated process satisfies the abovementioned requirements, that is, p _ii (s)N _i (s) has the smallest values in terms of time delay, relative degree and multiplicity of RHP zero in the ith row of the compensated process transfer function matrix. Consequently, the elements of V₁(s), V₂(s), and V₃(s) are designed based on the compensated process:

$$\{\begin{array}{c}{\delta }_i(s)=\frac{t_i(s)}{p_{\mathrm{ii}}(s)N_i(s)}\begin{array}{c}i=1,2,\dots ,m\hfill \end{array}\hfill \\ {\rho }_{\mathrm{ij}}(s)=-\frac{p_{\mathrm{ij}}(s)N_j(s)}{t_i(s)}\begin{array}{c}i,j=1,2,\dots ,m\hfill \end{array}\hfill \\ {\eta }_{\mathrm{ij}}(s)=-\frac{p_{i(m+j)}(s)N_{m+j}(s)}{t_i(s)}\begin{array}{c}i=1,2,\dots ,m,j=1,2,\dots ,n-m\hfill \end{array}\hfill \end{array}$$

(23)

• Remark 1. The compensators are solely used for meeting the realizability requirements. Therefore, there is no need to select large compensator parameters. In fact, smaller compensator parameters will be preferable as it can result in a simpler compensated process, which is friendly to the design of the controller later. In a specific design problem, the calculation of the compensator parameters can be described as the following linear programming problem:

$$\begin{array}{l}\min X=\sum _{i=1}^n{x}_i\\ s.t.\{\begin{array}{cc}\sum _{k=1}^n{a}_{ik}{x}_k\le b_i& i=1,2,\dots ,m\times (\mathrm{n}-1)\\ x_i\ge 0& i=1,2,\dots ,n\end{array}\end{array}$$

(24)

where x _i can be parameter α _i , ε _i or ${\mathrm{\Gamma }}_i^j$ ; a _ik and b _i are coefficients of a specific problem. The above linear programming problem can be solved with the simplex method.

• Remark 2. Equation (22) shows the general expression of the compensator. In most cases, the final compensator has a much simpler expression than equation (22). The reason is that, in process control, the elements of the process transfer function matrix are mostly described as first-order-plus-time-delay (FOPTD) models. Therefore, compensations for relative degree and RHP zeros are usually not needed. Generally, the compensators are introduced to compensate the time delays.

Finally, it is time to determine t _i (s) to finish the design of V₁(s) V₂(s), and V₃(s). Given the constraints of t _i (s), an intuitive and simple selection is t _i (s) = p _ii (s)N _i (s). This selection can result in the simplest form of V₁(s), that is, V₁(s) = I. Certainly, this is not the only way to select t _i (s). Sometimes, other simpler forms of t _i (s) are also available, but they will undoubtedly complicate V₁(s). With t _i (s) = p _ii (s)N _i (s), the elements of V₂(s) and V₃(s) can be derived accordingly as

$$\{\begin{array}{c}{\rho }_{\mathrm{ij}}(s)=-\frac{p_{\mathrm{ij}}(s)N_j(s)}{p_{\mathrm{ii}}(s)N_i(s)}\begin{array}{c}i,j=1,2,\dots ,m\hfill \end{array}\hfill \\ {\eta }_{\mathrm{ij}}(s)=-\frac{p_{i(m+j)}(s)N_{m+j}(s)}{p_{\mathrm{ii}}(s)N_i(s)}\begin{array}{c}i=1,2,\cdots ,m,j=1,2,\dots ,n-m\hfill \end{array}\hfill \end{array}$$

(25)

In the end, the implementation structure of the partial inverted decoupling control is shown in Figure 3, where C _d (s) and C _m (s) are the decentralized controller matrix of decoupled loops and non-decoupled loops, respectively; L _d (s) and L _m (s) are the disturbance transfer function vector of decoupled loops and non-decoupled loops, respectively; Y _d (s) and Y _m (s) are the output vector of decoupled loops and non-decoupled loops, respectively; R _d (s) and R _m (s) are the reference vector of Y _d (s) and Y _m (s), respectively.

Figure 3.

Implementation structure of partial inverted decoupling.

In Figure 3, the decentralized controllers C _d (s) and C _m (s) are still the unknown parts, which will be studied in the next section.

Design of the decentralized controller

Decentralized controllers C _d (s) and C _m (s) are located in the decoupled loops and the non-decoupled loops, respectively, to improve the performance of disturbance rejection. The design of C _d (s) and C _m (s) will rely on the partially decoupled process. According to Figure 3, the disturbance response [Y _d (s), Y _m (s)]^T can be expressed as

$$\begin{array}{c}\left[\begin{array}{c}{Y}_d(s)\hfill \\ \multicolumn{1}{c}{Y_m(s)}\end{array}\right]=\hfill \\ \left[\begin{array}{c}{(I+T_{11}(s)C_d(s))}^{-1}{L}_d(s)\hfill \\ \multicolumn{1}{c}{{(I+T_{22}(s)C_m(s))}^{-1}(L_m(s)-T_{21}(s)C_d(s){(I+T_{11}(s)C_d(s))}^{-1}{L}_d(s))}\end{array}\right]d(s)\hfill \end{array}$$

(26)

From equation (26), it can be observed that Y _d (s) is only related to C _d (s) while Y _m (s) is determined by both C _d (s) and C _m (s). Thus, the decentralized controllers C _d (s) and C _m (s) should be designed in sequence. Firstly, C _d (s) will be designed by optimizing Y _d (s). And then the design of C _m (s) can be undertaken for optimizing Y _m (s) with the determined C _d (s).

Design of the controllers in the decoupled loops

As the couplings in the decoupled loops have been removed, the decoupled system can be regarded as a series of independent subsystems. The ith element of C _d (s), ${[C_d(s)]}_{\mathrm{ii}}$ is the sub-controller of the ith subsystem. The design of ${[C_d(s)]}_{\mathrm{ii}}$ is carried out by using the internal model control (IMC) structure shown in Figure 4, where ${[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}$ is the equivalent internal model controller. ${[C_d(s)]}_{\mathrm{ii}}$ can be obtained based on equation (27) after ${[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}$ is determined.

$${[C_d(s)]}_{\mathrm{ii}}=\frac{{[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}}{1-t_i(s){[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}}$$

(27)

Figure 4.

Internal model control structure of the ith subsystem.

According to the IMC structure, the disturbance response $[Y_d(s){]}_i$ is

$${[Y_d(s)]}_i=(1-t_i(s){[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}){[L_d(s)]}_{i}d(s)$$

(28)

where t _i (s) and $[L_d(s){]}_i$ are assumed as FOPTD models:

$$\begin{array}{cc}{t}_i(s)=\frac{K_{\mathrm{ti}}}{T_{\mathrm{ti}}s+1}{e}^{-{\tau }_{i}s},\hfill & {[L_d(s)]}_i=\frac{K_{\mathrm{Li}}}{T_{\mathrm{Li}}s+1}{e}^{-{\vartheta }_{i}s}\hfill \end{array}$$

(29)

Assuming d(s) = 1/s, the necessary condition for steady-state performance of disturbance rejection is

$$\underset{s\to 0}{lim}(1-t_i(s){[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}){[L_d(s)]}_i=0$$

(30)

Based on equation (30) and IMC principle, the potential internal model controller can be designed as

$${[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}=t_{i-}^{-1}(s)·J_i(s)$$

(31)

where $t_{i-}^{-1}(s)$ is the inverse of the minimum phase part of $t_i(s)$ and $J_i(s)$ is the undetermined low-pass filter which is used to make the controller realizable and improve the control performance.

From equation (28), obviously, if slow modes exist in the disturbance channel, the performance for disturbance rejection will be deteriorated. To speed up disturbance rejection, $J_i(s)$ is designed using equation (32), which has two undetermined parameters, that is, χ _i and λ _i .

$$J_i(s)=\frac{{\chi }_{i}s+1}{{({\lambda }_{i}s+1)}^2}$$

(32)

In order to eliminate the slow mode, the undesired poles of ${[L_d(s)]}_i$ should be cancelled. As a result, χ _i should satisfy

$${\chi }_i=\frac{2{\lambda }_i·T_{\mathrm{Li}}+T_{\mathrm{Li}}·{\tau }_i-{\lambda }_i^2}{T_{\mathrm{Li}}+{\tau }_i}$$

(33)

The proof of equation (33) is shown in the Appendix.

From equation (33), χ _i can be calculated when a concrete λ _i is given. Hence, the selection of λ _i will be focused on next. Obviously, in terms of speeding up disturbance rejection, a small value of λ _i is preferable. Meanwhile, a small value of λ _i can also provide a small gain for the disturbance output, which can suppress the peak value of the disturbance output (The gain of disturbance output can be approximately described as $K_{\mathrm{Li}}(2{\lambda }_i+{\tau }_i-{\chi }_i)=K_{\mathrm{Li}}{({\lambda }_i+{\tau }_i)}^2/(T_{\mathrm{Li}}+{\tau }_i)$ ). Besides performance, robustness should also be considered due to the inevitable uncertainties. A widely-used robustness index is the maximum sensitivity (M _s ):

$$M_{\mathrm{si}}=\Vert 1-t_i(s){[C_{\mathrm{imc}}(s)]}_{\mathrm{ii}}{\Vert }_{\infty }$$

(34)

A large value of M _s reveals bad robustness and a typical range of M _s is [1.4, 2.0]. For ease of analyzing the relationship between λ _i and M _si , define Λ _i  = λ _i /T _Li , Ξ _i  = T _Li /τ _i , and Ω _i  = ω·τ _i , where ω is the frequency. Then, M _si in equation (34) can be rewritten as

$$\begin{array}{c}{M}_{\mathrm{si}}=\underset{{\mathrm{\Omega }}_i}{max}\hfill \\ |1-\frac{1+\frac{\mathbf{j}{\mathrm{\Omega }}_i(2{\mathrm{\Lambda }}_i{\mathrm{\Xi }}_i^2+{\mathrm{\Xi }}_i-{\mathrm{\Lambda }}_i^2{\mathrm{\Xi }}_i^2)}{1+{\mathrm{\Xi }}_i}}{1-{({\mathrm{\Lambda }}_i·{\mathrm{\Xi }}_i·{\mathrm{\Omega }}_i)}^2+2\mathbf{j}({\mathrm{\Lambda }}_i·{\mathrm{\Xi }}_i·{\mathrm{\Omega }}_i)}{e}^{-\mathbf{j}{\mathrm{\Omega }}_i}|\begin{array}{c}{\mathrm{\Omega }}_i\in [0,+\infty )\hfill \end{array}\hfill \end{array}$$

(35)

where j is the imaginary unit. Based on equation (35), the relationship between system parameters and M _s is depicted in Figure 5.

Figure 5.

Relationship between system parameters and M _s : (a) relationship between Λ, Ξ and Ms and (b) relationship between Λ and Ξ for M _s  = 2.0.

According to Figure 5(a), for a specific Ξ, M _s decreases while Λ increases. This reveals that a large λ (Λ) can improve the robustness. Therefore, with respect to the selection of λ, there is a tradeoff between robustness and performance. Meanwhile, for a specific Λ, M _s decreases while Ξ increases. In light of this, it can be inferred that if a system has a relatively large time delay, controller parameter λ ought to be increased to guarantee sufficient robustness. Figure 5(b) shows the relationship between Λ and Ξ for M _s  = 2.0 and an approximate expression has been obtained as equation (36). In this way, the value of M _s will lie in the typical range, when λ ≥ 1.254T _L /(1.672 + T _L /τ).

$$\left(\lambda /T_L\right)=\frac{1.254}{1.672+\left(T_L/\tau \right)}$$

(36)

The design of internal model controller can be finished along with the determination of λ _i . Furthermore, based on equation (27), the equivalent feedback controller which is in the form of a PID controller can be fulfilled as

$${[C_d(s)]}_{\mathrm{ii}}=K_{\mathrm{Pi}}(1+\frac{1}{T_{\mathrm{Ii}}s}+T_{\mathrm{Di}}s)\frac{1}{T_{\mathrm{Fi}}s+1}$$

(37)

where

$$\{\begin{array}{c}{K}_{\mathrm{Pi}}=\frac{(T_{\mathrm{ti}}+{\chi }_i)}{K_{\mathrm{ti}}(2{\lambda }_i+{\tau }_i-{\chi }_i)}\hfill \\ T_{\mathrm{Ii}}=T_{\mathrm{ti}}+{\chi }_i\hfill \\ T_{\mathrm{Di}}=\frac{T_{\mathrm{ti}}·{\chi }_i}{T_{\mathrm{ti}}+{\chi }_i}\hfill \\ T_{\mathrm{Fi}}=\frac{{\lambda }_i^2+{\tau }_i·{\chi }_i}{(2{\lambda }_i+{\tau }_i-{\chi }_i)}\hfill \end{array}$$

(38)

Design of the controllers in the non-decoupled loops

After C _d (s) is obtained, the decentralized controllers in the non-decoupled loops, C _m (s) will be designed then. Similarly, the disturbance output Y _m (s) can be denoted as

$$Y_m(s)={(I+T_{22}(s)C_m(s))}^{-1}{L}_e(s)d(s)$$

(39)

where L _e (s) and T₂₂(s) are the equivalent disturbance channel model and the equivalent controlled process model, respectively. L _e (s) and T₂₂(s) have the following expressions:

$$\{\begin{array}{c}{T}_{22}(s)=(P_{22}(s)-P_{21}(s)P_{11}^{-1}(s)P_{12}(s))N_M(s)\hfill \\ \multicolumn{1}{c}{L_e(s)=L_m(s)-P_{21}(s)P_{11}^{-1}(s)\mathrm{\Upsilon }(s)J(s)L_d(s)}\end{array}$$

(40)

where

$$\{\begin{array}{c}{N}_M(s)=\mathrm{diag}\{N_{m+1}(s),N_{m+2}(s),\dots ,N_n(s)\}\hfill \\ \mathrm{\Upsilon }(s)=\mathrm{diag}\{e^{-{\tau }_{1}s},e^{-{\tau }_{2}s},\dots ,e^{-{\tau }_{m}s}\}\hfill \\ J(s)=\mathrm{diag}\{J_1(s),J_2(s),\dots ,J_m(s)\}\hfill \end{array}$$

(41)

According to equation (39), the design of C _m (s) relies on L _e (s) and T₂₂(s). For the independent design of the decentralized controllers, the idea of EOTF is used. ¹³ The EOTFs of the process model and disturbance model are

$$\{\begin{array}{c}{g}_i^{\mathrm{eff}}(s)={[T_{22}(s)]}_{\mathrm{ii}}-{[T_{22}(s)]}_{i*}·{[T_{22}(s)]}_{-\mathrm{ii}}^{-1}·{[T_{22}(s)]}_{*i}\hfill \\ \multicolumn{1}{c}{L_i^{\mathrm{eff}}(s)={[L_e(s)]}_i-{[T_{22}(s)]}_{i*}·{[T_{22}(s)]}_{-\mathrm{ii}}^{-1}·{[L_e(s)]}_{-i}}\end{array}$$

(42)

where ${[T_{22}(s)]}_{\mathrm{ii}}$ is the ith row and ith column element of T₂₂(s); ${[T_{22}(s)]}_{i*}$ and ${[T_{22}(s)]}_{*i}$ are the ith row and ith column of T₂₂(s) with ${[T_{22}(s)]}_{\mathrm{ii}}$ removed, respectively; ${[T_{22}(s)]}_{-\mathrm{ii}}$ is T₂₂(s) removing its ith row and ith column; ${[L_e(s)]}_i$ is the ith element of L _e (s) and ${[L_e(s)]}_{-i}$ is L _e (s) removing ${[L_e(s)]}_i$ . Apparently, both $g_i^{\mathrm{eff}}(s)$ and $L_i^{\mathrm{eff}}(s)$ in equation (42) are too complex to be used for controller design. Therefore, model reduction technique is usually used to get low-order models. The approximate FOPDT model can be obtained based on the following method ¹⁰ :

$$g_i^{\mathrm{eff}}(s)\approx \frac{{\overline{K}}_{\mathrm{gi}}{e}^{-{\overline{\tau }}_{\mathrm{gi}}s}}{{\overline{T}}_{\mathrm{gi}}s+1}$$

(43)

where

$$\{\begin{array}{c}{\overline{K}}_{\mathrm{gi}}=g_i^{\mathrm{eff}}(0)\hfill \\ \multicolumn{1}{c}{{\overline{T}}_{\mathrm{gi}}=\sqrt{\frac{{\overline{K}}_{\mathrm{gi}}^2-{|g_i^{\mathrm{eff}}(j{\omega }_{\mathrm{ci}})|}^2}{{\omega }_{\mathrm{ci}}^2{|g_i^{\mathrm{eff}}(j{\omega }_{\mathrm{ci}})|}^2}}}\\ \multicolumn{1}{c}{{\overline{\tau }}_{\mathrm{gi}}=\frac{\pi -\arctan \left({\overline{T}}_{\mathrm{gi}}{\omega }_{\mathrm{ci}}\right)}{{\omega }_{\mathrm{ci}}}}\end{array}$$

(44)

and ω _ci is the phase cross-over frequency of $g_i^{\mathrm{eff}}(s)$ . Similarly, $L_i^{\mathrm{eff}}(s)$ can be approximated as $L_i^{\mathrm{eff}}(s)\approx {\overline{K}}_{\mathrm{Li}}{e}^{-{\overline{\tau }}_{\mathrm{Li}}s}/({\overline{T}}_{\mathrm{Li}}s+1)$ . Then, the ith sub-controller ${[C_m(s)]}_{\mathrm{ii}}$ will be designed on the basis of the approximate models of $g_i^{\mathrm{eff}}(s)$ and $L_i^{\mathrm{eff}}(s)$ . The design procedure can consult the design of ${[C_d(s)]}_{\mathrm{ii}}$ .

Simulation study

In this section, two examples are provided to demonstrate the effectiveness of the proposed method.

• Example 1. Industrial scale polymerization (ISP) reactor

The process transfer function matrix of ISP reactor is

$$\begin{array}{c}\multicolumn{1}{c}{P(s)=\left[\begin{array}{cc}\frac{22.89e^{-0.2s}}{4.572s+1}\hfill & \frac{-11.64e^{-0.4s}}{1.807s+1}\hfill \\ \multicolumn{1}{c}{\frac{4.689e^{-0.2s}}{2.174s+1}}& \frac{5.8e^{-0.4s}}{1.801s+1}\hfill \end{array}\right]}\end{array}$$

(45)

The disturbance transfer function vector is assumed to be

$$L(s)={\left[\begin{array}{cc}\frac{22.89e^{-0.2s}}{4.572s+1}\hfill & \frac{5.8e^{-0.4s}}{1.801s+1}\hfill \end{array}\right]}^{\mathrm{T}}$$

(46)

First, the couplings should be evaluated based on RDG. According to equation (6), the values of RDG are computed as |β₁| = 1.069 and |β₂| = 0.1357, which indicates that the couplings in loop 1 need to be removed while the couplings in loop 2 need not. Therefore, the decoupled process is determined as

$$T(s)=\left[\begin{array}{cc}{t}_{11}(s)\hfill & 0\hfill \\ \multicolumn{1}{c}{t_{21}(s)}& t_{22}(s)\hfill \end{array}\right]$$

(47)

Then, t₁₁(s) will be determined to design the components of the partial inverted decoupler, that is, V₁(s), V₂(s), and V₃(s). Fortunately, p₁₁(s) meets the requirements given in Section “The principle of determining T₁₁(s),” which means the extra compensation is unnecessary. In this case, N(s) = I and t₁₁(s) = p₁₁(s). Substituting t₁₁(s) into equations (16) and (18), the components are obtained as V₁(s) = 1, V₂(s) = 0 and V₃(s) =  $-p_{12}(s)/p_{11}(s)$ . The next step is to design the decentralized controllers based on the decoupled process. For the decoupled loop, C _d (s) is designed based on t₁₁(s) and L _d (s), where t₁₁(s) = L _d (s) = 22.89e^−0.2s/(4.572s + 1). According to our proposed method, the controller parameter λ _{d 1} should satisfy λ _{d 1} > 0.2337, where λ _{d 1} is determined as λ _{d 1} = 0.7 to achieve better robustness. With λ _{d 1}, parameter χ _{d 1} = 1.4303 is obtained by equation (33) and C _d (s) is finally determined as

$$C_d(s)=(1.545+\frac{0.2574}{s}+1.683s)\frac{1}{4.572s+1}$$

(48)

For the non-decoupled loop, C _m (s) can be designed with the approximately equivalent models $g_1^{\mathrm{eff}}(s)\approx 8.184e^{-0.4036s}/(1.360s+1)$ and $L_1^{eff}(s)\approx 1.111e^{-0.3776s}/\left(0.3469s+1\right)$ . Similarly, the final C _m (s) is determined as

$$C_m(s)=(0.3157+\frac{0.1853}{s}+0.0867s)\frac{1}{0.3469s+1}$$

(49)

with λ _{m 1} = 0.3 for better robustness. With the decentralized controllers, the output response is shown in Figure 6(a). The proposed method is compared with several classical control strategies for multivariable processes, such as decentralized control, ⁵² centralized control ⁵³ and decoupling control. ²² In Figure 6(a), it can be seen that the proposed method can achieve better performances for disturbance rejection. The peak of disturbance outputs of the proposed method are smaller. Meanwhile, the proposed method can bring the outputs back to their set-points in shorter time, especially in loop 2. Moreover, the integral absolute error (IAE) defined as

$$\mathrm{IAE}=\sum _{i=1}^n\mathrm{IA}{E}_i$$

(50)

where

$$\mathrm{IA}{E}_i={\int }_0^{\infty }|y_i(t)|\mathrm{dt}$$

(51)

is also used to evaluate the control performance. The IAE values of different control strategies are listed in Table 1 for comparison. It is noticed that the smallest IAE value is provided by the proposed strategy. In addition to the nominal control performance, it is also important to evaluate the robustness of the control system. Therefore, +20% parameter perturbations are introduced to the controlled process to evaluate the robust performance. Figure 6(b) depicts the disturbance outputs of the system with parameter perturbations. It can be seen that the proposed method can keep the better performances in terms of the output peak, response time and IAE value (The IAE values of the perturbed system are also listed in Table 1 and denoted by IAE_p). Then, the robust stability will be evaluated. For a multivariable system, both input uncertainty and output uncertainty should be considered. According to the M-Δ structure, in order to ensure robust stability, the multiplicative input uncertainty Δ _I (s) and the multiplicative output uncertainty Δ _O (s) should satisfy

$$\begin{array}{c}\multicolumn{1}{c}{\{\begin{array}{cc}\Vert {\mathrm{\Delta }}_I(\mathrm{j}\omega )\Vert <1/\overline{\sigma }(M_I(\mathrm{j}\omega ))\hfill & \forall \omega \hfill \\ \multicolumn{1}{c}{\Vert {\mathrm{\Delta }}_O(\mathrm{j}\omega )\Vert <1/\overline{\sigma }(M_O(\mathrm{j}\omega ))}& \forall \omega \hfill \end{array}}\end{array}$$

(52)

where

$$\{\begin{array}{c}{M}_I(\mathrm{j}\omega )={(I+Q(\mathrm{j}\omega )C(\mathrm{j}\omega )P(j\omega ))}^{-1}Q(\mathrm{j}\omega )C(\mathrm{j}\omega )P(\mathrm{j}\omega )\hfill \\ M_O(\mathrm{j}\omega )={(I+P(\mathrm{j}\omega )Q(\mathrm{j}\omega )C(j\omega ))}^{-1}P(\mathrm{j}\omega )Q(\mathrm{j}\omega )C(\mathrm{j}\omega )\hfill \end{array}$$

(53)

and C(s) = diag{C _d (s), C _m (s)}. The curves of $1/\overline{\sigma }(M_I(\mathrm{j}\omega ))$ and $1/\overline{\sigma }(M_O(\mathrm{j}\omega ))$ are plotted in Figure 6(c). The minimums of $1/\overline{\sigma }(M_I(\mathrm{j}\omega ))$ and $1/\overline{\sigma }(M_O(\mathrm{j}\omega ))$ at all frequencies denote the upper bounds of $\Vert {\mathrm{\Delta }}_I(\mathrm{j}\omega )\Vert$ and $\Vert {\mathrm{\Delta }}_O(\mathrm{j}\omega )\Vert$ . Therefore, generally, the minimums of $1/\overline{\sigma }(M_I(\mathrm{j}\omega ))$ and $1/\overline{\sigma }(M_O(\mathrm{j}\omega ))$ , γ _I and γ _O are used as robustness indexes. ²⁶ Larger values of γ _I or γ _O reveal better robust stability. From Table 1, the proposed method provides the biggest γ _O and an acceptable γ _I , which means a good robust stability is achieved.

• Example 2. Industrial distillation column

Figure 6.

Control performance for example 1: (a) output responses for example 1 on nominal condition, (b) output responses for example 1 on perturbed condition, and (c) curves of 1/||M(jω)|| for example 1.

Table 1.

Performance indexes of example 1.

Method	IAE1	IAE2	IAE	IAEp	γ I	γ O
Proposed	3.8895	1.0560	4.9455	5.3752	0.4443	0.6669
Decoupling	3.4244	14.8361	18.2605	16.2709	0.5825	0.5973
Decentralized	9.4711	4.1524	13.6235	16.9606	0.2740	0.3775
Centralized	39.2103	4.6024	43.8127	42.1012	0.7989	0.5850

Here, a widely-studied 3 × 3 industrial distillation column is considered as

$$\begin{array}{c}\multicolumn{1}{c}{P(s)=\left[\begin{array}{ccc}\frac{1.986e^{-0.71s}}{66.7s+1}\hfill & \frac{-5.24e^{-60s}}{400s+1}\hfill & \frac{-5.984e^{-2.24s}}{14.29s+1}\hfill \\ \multicolumn{1}{c}{\frac{-0.0204e^{-0.59s}}{{(7.14s+1)}^2}}& \frac{0.33e^{-0.68s}}{{(2.38s+1)}^2}\hfill & \frac{-2.38e^{-0.42s}}{{(1.43s+1)}^2}\hfill \\ \multicolumn{1}{c}{\frac{-0.374e^{-7.75s}}{22.22s+1}}& \frac{11.3e^{-3.79s}}{{(21.74s+1)}^2}\hfill & \frac{9.811e^{-1.59s}}{11.36s+1}\hfill \end{array}\right]}\end{array}$$

(54)

with an assumed disturbance transfer function vector

$$L(s)={\left[\begin{array}{ccc}\frac{1.986e^{-0.71s}}{66.7s+1}\hfill & \frac{0.33e^{-0.68s}}{{(2.38s+1)}^2}\hfill & \frac{9.811e^{-1.59s}}{11.36s+1}\hfill \end{array}\right]}^{\mathrm{T}}$$

(55)

First, RDGs are calculated as |β₁| = 3.5827 and |β₂| = 1.0120 and |β₃| = 0.0290 to evaluate couplings. Consequently, the couplings of the first two loops need to be removed while that in the last loop need not. As a result, the decoupled process is determined as

$$T(s)=\left[\begin{array}{ccc}{t}_{11}(s)\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& t_{22}(s)\hfill & 0\hfill \\ \multicolumn{1}{c}{t_{31}(s)}& t_{32}(s)\hfill & t_{33}(s)\hfill \end{array}\right]$$

(56)

In this case, a compensator designed as N(s) =diag{e^−0.09s, 1, e^−0.26s} is needed as the requirement for the time delay cannot be met if t₁₁(s) = p₁₁(s) and t₂₂(s) = p₂₂(s) are selected directly. Based on the compensated process, V₁(s), V₂(s), and V₃(s) are achieved as V₁(s) = I,

$$V_2(s)=\left[\begin{array}{cc}0\hfill & \frac{2.6385(66.7s+1)e^{-59.2s}}{400s+1}\hfill \\ \multicolumn{1}{c}{\frac{0.06182{(2.38s+1)}^2}{{(7.14s+1)}^2}}& 0\hfill \end{array}\right]$$

(57)

and V₃(s) =  $[\begin{array}{cc}3.0131e^{-1.7s}(66.7s+1)/(14.29s+1),\hfill & 7.2121{(2.38s+1)}^2/{(1.43s+1)}^2\hfill \end{array}{]}^{\mathrm{T}}$ , respectively. Next, the decentralized controllers for the partially decoupled process are designed. According to the proposed method, the decentralized controllers for the decoupled loops are

$$\begin{array}{c}{C}_d(s)=\mathrm{diag}\{(110.8718+\frac{1.4752}{s}+832.2807s)\frac{1}{66.7s+1}\hfill \\ (2.3626+\frac{0.3493}{s}+3.9828s)\frac{1}{3.57s+1}\}\hfill \end{array}$$

(58)

Then, the controller of the non-decoupled loop, C _m (s) is designed based on the approximate equivalent models, $g_1^{\mathrm{eff}}(s)\approx 99.856e^{-2.530s}/(123.746s+1)$ and $L_1^{\mathrm{eff}}(s)\approx -2.90(-124.855s+1)e^{-5.597s}/(410.506s^2+40.045s+1)$ . Here, a second-order-plus-time-delay (SOPTD) model is used to approximate $L_1^{\mathrm{eff}}(s)$ , because the dynamics of $L_1^{\mathrm{eff}}(s)$ are too complex. With the similar thought, C _m (s) is finally determined as

$$\begin{array}{c}\multicolumn{1}{c}{C_m(s)=(0.07578+\frac{0.002637}{s}+0.5824s)}\\ \frac{123.746s+1}{410.506s^2+40.045s+1}\hfill \end{array}$$

(59)

The output responses of the nominal process are shown in Figure 7(a) compared with Liu’s ⁵⁴ analytical decoupling control strategy, ETF-based partial decoupling control strategy ⁴⁹ and controller synthesis based on RDGA. ⁴⁷ The corresponding IAE values are listed in Table 2. Obviously, the proposed method provide better output performances than other methods. Similarly, to evaluate robust performance, −20% parameter perturbations are introduced to the controlled process and the corresponding output responses are depicted in Figure 7(b). The proposed method can still achieve the best IAE value for the perturbed process. Besides, the curves of $1/\overline{\sigma }(M_I(j\omega ))$ and $1/\overline{\sigma }(M_O(j\omega ))$ and the robustness indexes, γ _I and γ _O are also shown in Figure 7(c) and Table 2, respectively. It can be seen that the proposed method has better robustness than ETF-based partial decoupling control strategy. Liu’s method achieve bigger γ _I and γ _O , but it suffers from complex calculation and controller form. The controllers based on RDGA achieve the biggest γ _I , but it results in worse output responses.

Figure 7.

Control performance for example 2: (a) output responses for example 2 on nominal condition, (b) output responses for example 2 on perturbed condition, and (c) curves of 1/||M(jω)|| for example 2.

Table 2.

Performance indexes of example 2.

Method	IAE1	IAE2	IAE3	IAE	IAEP	γ I	γ O
Proposed	0.737	2.889	72.066	75.692	102.555	0.05592	0.2746
Liu	12.292	3.244	77.625	93.161	119.853	0.1559	0.7212
ETF	85.242	1.174	1172.6	1259.016	1137.189	0.01044	0.1495
RDGA	496.268	3.418	544.710	1044.396	1314.852	0.6539	0.2628

Conclusion

In this work, a novel partial inverted decoupling control strategy to improve the disturbance rejection of multivariable process with time delays is presented. Two main issues: the design of the partial decoupler and the decentralized controllers, are focused on. Motivated by the thought of using favorable couplings to reject disturbances, the couplings are firstly evaluated with RDGs. Then the inverted decoupling usually used in full decoupling control is extended to partial decoupling control. The elements of the partial inverted decoupler are designed in detail and the implementation structure is presented. The advantage of the proposed partial inverted decoupler is that the decoupler is implemented without approximation, which can remove the unfavorable couplings thoroughly. Based on the partial decoupled processes, the decentralized controllers for both the decoupled and non-decoupled loops are designed. An analytical design method is proposed with due consideration of the performance of disturbance rejection and robustness. Moreover, to improve the control performance of the non-decoupled loop, the equivalent process model and the equivalent disturbance channel model are obtained for controller design. Two simulation examples are provided and the results show that the proposed method can achieve better nominal control performance and robustness compared with several existing multivariable control strategies.