Novel Centralized IMC-PID Controller Design for Multivariable

Apr 10, 2017 - decoupler will increase the cost and the complexity of the control system. ...... The ISP reactor is a TITO (two inputs and two outputs...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/IECR

Novel Centralized IMC-PID Controller Design for Multivariable Processes with Multiple Time Delays Qibing Jin, Xinghan Du,* and Beiyan Jiang Institute of Automation, Beijing University of Chemical Technology, Beisanhuan East Road 15, Chaoyang District, Beijing 100029, PR China ABSTRACT: In this paper, a novel centralized IMC (internalmodel-control)-PID (proportion-integration-differentiation) controller is designed for multivariable processes with multiple time delays. To reduce the influence of time delays, the process model is divided into the minimum common time delay part and the remaining part. The IMC-PID controller, consisting of a full-matrix filter and a diagonal PID controller, is obtained based on the IMC principle. The full-matrix filter described as the inverse of the remaining part is too complex to compute. Hence, the Nyquist set is employed to describe the dynamics of the filter. The approximate transfer function matrix of the filter is calculated by complex curve fitting. Besides, to reduce the approximation error and obtain stable filter elements, a compensator is introduced. The rules of selecting compensator parameters are proposed for stable filter elements. Finally, the IMC-PID controller parameters are determined after comprehensively considering the performance and the robustness. The limitation of the proposed method is discussed. Simulations demonstrate the effectiveness of the proposed method.

1. INTRODUCTION Multivariable processes widely exist in practical industries. Different from SISO (single-input−single-output) processes, MIMO (multiple-input−multiple-output) processes are difficult to control because of the interactions among the inputs and outputs. For SISO systems, many mature methods have been studied, such as the Z-N method1 and the SIMC method.2 Therefore, the direct idea is to extend these methods to MIMO processes. In fact, when the interactions are modest, the MIMO control systems can obtain acceptable performance just by controlling the selected main channel processes. However, it is difficult to achieve a good output performance when the interactions are strong. Therefore, many control strategies for MIMO processes have sprung up. The decentralized control strategy is widely used in MIMO systems because it is easy to implement. The decentralized controllers are designed based on the main channel processes, and the interactions are regarded as the disturbances. To reduce the influence of the interactions, the selection of the main channels is very important. The relative gain array (RGA) was proposed and used to pair the outputs with the appropriate inputs.3,4 However, the RGA is calculated on the basis of the process steady-state gains, which means the dynamic information on the process is overlooked. Therefore, the RGA may give wrong pairing results sometimes. Considering the limitation of the RGA, the dynamic information on the process has to be considered, which leads to the dynamic relative gain array (DRGA). The DRGA was first proposed by Witcher and Mcavoy.5 They used transfer functions instead of steady-state © XXXX American Chemical Society

gains used for the RGA. Besides, Mc Avoy et al. defined a novel DRGA based on the proportional output optimal controller gain matrix which was designed by a state space approach.6 This DRGA had to be dependent on the controller, which increased the complexity of the calculation. To obtain an effective and digestible pairing method, Xiong et al. presented a new dynamic pairing criterion named the effective relative gain array (ERGA) using both the steady-state gain and the bandwidth information on process transfer functions.7 The bandwidth information on the process can reflect the dynamics of the process. Furthermore, the relative normalized gain array (RNGA) proposed by He et al. is also an alternative pairing criterion.8 Different from the ERGA, the average residence time of the normalized process is used to describe the transient information in the RNGA. Therefore, it is useful for processes with large time delays. The aforementioned paring criteria are just used to determine the main loops. Then the controllers are tuned based on the main loops. In fact, to improve the performance of control systems, the controllers might as well be designed dependent on both the main channels and the interaction channels. Therefore, the effective open-loop transfer function (EOTF) was proposed for controller design.9,10 However, the EOTF is difficult to calculate for the processes with more than three inputs and outputs, and model approximation cannot be avoided for simplifying controller design, which Received: December 27, 2016 Revised: March 13, 2017 Accepted: March 28, 2017

A

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

time delay part and the remaining part. The resulting centralized PID controller is composed of a full-matrix filter and a diagonal PID controller. To compute the full-matrix filter described as the inverse of the remaining part, the Nyquist set is used to describe the dynamics of the filter. Then the approximate transfer function matrix of the filter is calculated by complex curve fitting. A compensator is introduced to obtain the feasible and stable filters. The rules of selecting compensator parameters are proposed. Finally, the PID controller parameters are determined after comprehensively considering the performance and the robustness. Two examples are given to show the design procedures and demonstrate the effectiveness of the proposed method for both square and nonsquare processes. This paper is organized as follows. Section 2 provides the detailed procedure of designing the centralized IMC-PID controllers, including the derivation of the controller, the design of the full-matrix filter, and the determination of the controller parameters. In section 3, two simulation examples, including a square process and a nonsquare process, are given to determine the effectiveness of the proposed controller. The conclusion is given in section 4.

may deteriorate control performance. In sum, the decentralized control strategy has its performance limitations,11,12 which means the decentralized control strategy is not appropriate when interactions are strong. Apart from the decentralized control strategy, the decoupling control strategy may be another effective control strategy for MIMO processes. The essence of decoupling is to add a compensator (also known as a decoupler) between the controller and the process to make the compensated (decoupled) process a diagonal process. There is no doubt that the employ of the decoupler will increase the cost and the complexity of the control system. But it surely reduces the influence of interactions. The typical decoupling methods are the ideal,13,14 simplified,15,16 inverted,17−19 and static decoupling.20 The static decoupling just considers the interactions at steady state. Hence, the decoupling performance for the dynamic process is poor. Its advantage is that the decoupler is easy to obtain. The ideal decoupling uses the selected decoupled process and the inverse of the process to calculate the decoupler. This method can achieve good dynamic decoupling performance, but the decoupler is difficult to get, especially for the process with time delays. Conversely, the decoupler of the simplified decoupling is easy to get, but the decoupled process which is used to design the controller is too complex. The inverted decoupling can get the simple decoupler and the decoupled process simultaneously, but it has many constraints for getting casual, proper, and stable decoupler elements. Therefore, how to design a more effective dynamic decoupler is still worth studying. When time delays exist in MIMO processes, especially multiple time delays, the control of MIMO processes becomes more difficult. The Smith predictor structure can deal with time delays effectively,21−23 but it is not easy to implement in practice because of its complicated structure and sensitivity for model uncertainty. Internal model control (IMC) proposed by Garcia and Morari is also an effective method for the process with time delay.24 Besides, the IMC controller can be implemented in the form of a proportion-integration-differentiation (PID) controller.25 When IMC is used in MIMO systems, the controller is not easy to obtained because the inverse-based part is difficult to calculate.26 Hence, generally, the IMC controller is designed on the background of the multiloop control and the decoupling control. The multiloop IMC controller is designed based on the loop pairing.27 Moreover, the EOTF can also be used to design the IMC controller.28,29 The essence of these methods is using the IMC principle to design the decentralized controllers. The decoupling IMC controller is designed based on the decoupled process.30,31 The key point is to seek an effective decoupler and a feasible decoupled process. Liu and Zhang proposed a decoupling IMC controller based on the adjoint and determinant of the process for a 2 × 2 and n × n process. However, the controller is too complex when n is large.32,33 Wang et al. also proposed a decoupling IMC controller for square processes.34 Its advantage is that only the controllers in the diagonal position need to be designed. The off-diagonal controllers can be obtained based on the diagonal controllers. However, the equivalent SISO process is very complex, which leads to complex IMC controllers. This work aims to propose simple and effective centralized IMC-PID controllers for multivariable processes with multiple time delays. The centralized PID controllers are derived based on the basic IMC structure rather than the selected main processes or the decoupled processes. To reduce the influence of time delays, the process model is divided into the minimum common

2. CONTROLLER DESIGN 2.1. IMC-PID controller for multivariable systems. The basic internal model control structure is shown in Figure 1.

Figure 1. Basic internal model control structure.

Gp and Gm are the process and model, respectively. CIMC is the internal model controller. Y (s) = Gm(s)CIMC(s)R(s)

(1)

CIMC(s) = Gm−1(s)

(2)

CIMC(s) = Gm−−1 (s)F(s)

(3)

From Figure 1, it can be seen that when the model matches the process, the IMC structure is equivalent to an open-loop control structure. Then, the closed-loop transfer function can be simplified as eq 1. Therefore, the ideal IMC controller can be designed as eq 2. However, taking the stability and realizability of the IMC controller into consideration, the IMC controller should be rewritten as eq 3. Gm‑ denotes the minimumphase part of the model Gm, and F denotes the low-pass filters whose parameters can be tuned to obtain the desired performance. C(s) = CIMC(s)[I − Gm(s)CIMC(s)]−1

(4)

The IMC structure is not easy to put into practice. The equivalent feedback control structure (Figure 2) of Figure 1 can solve this problem. In Figure 2, the controller in the dashed box is the equivalent feedback controller, which can be designed as a PID controller. The relationship between the IMC controller in Figure 1 and the feedback controller in Figure 2 can be seen in eq 4. Equation 4 shows that the feedback B

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research τi • ⎧ ⎪ KPi = λ + τ i i• ⎪ ⎪ ⎨ TIi = τi • ⎪ ⎪ T = λiτi • ⎪ fi λ + τ ⎩ i i•

Based on eq 4 and eq 9, the equivalent feedback controller is shown as eq 10. The feedback controller shown in eq 10 can be implemented in the form of a PID controller. Taking e−τs ≈ 1/τs + 1,36 the equivalent PID controller is provided as eqs 11 and 12. From eq 11, it can be seen that the multivariable PI controller consists of two parts. The first part is the diagonal PI controller which is used to eliminate the errors between the outputs and the set-points. The second part, G−1 m0(s), is the fullmatrix filter which is used to weaken coupling and improve dynamic performance. The diagonal PI controller is relatively easy to implement. Therefore, the difficulty is how to obtain the effective full-matrix filter. 2.2. Full-matrix filter design. The most important function of the full-matrix filter is to weaken the coupling among each loop. Obviously, it is still difficult to calculate G−1 m0(s) directly. −1 A usual method is using G−1 m0(s = 0) instead of Gm0(s). However, this method just takes the steady-state into consideration. If we want to improve the performance for dynamic decoupling, the filter must work at more frequencies. Here, the Nyquist set37 is used to describe the dynamic of the full-matrix filter. Considering the implementation of the filter, the curve fitting technology is used to obtain an approximate transfer function matrix based on the Nyquist set. However, the process model we consider here is shown as eqs 5 and 6. The relative degree of the element gij(s) is 1. Generally, the relative degree of the element in G−1 m0(s) is −1 on this condition. The element in G−1 (s) cannot be m0 approximated by a proper transfer function.

Figure 2. Equivalent feedback structure of IMC.

controller can be determined after the IMC controller is determined. Therefore, the problem of designing the feedback controller can be translated into the problem of designing a good IMC controller. ⎡ g (s ) g (s ) 12 ⎢ 11 ⎢ g (s ) g (s ) 22 Gp(s) = ⎢ 21 ⎢ ⋮ ⋮ ⎢ ⎢⎣ gn1(s) gn2(s)

gij(s) =

kij Tijs + 1

⋯ g1n(s)⎤ ⎥ ⋯ g2n(s)⎥ ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ gnn(s)⎥⎦

(5)

e−τijs (6)

Both the IMC controller and its equivalent feedback controller are easy to obtain for a SISO system because the process model is just a scalar transfer function. However, in MIMO systems, the process is described as a transfer function matrix which is shown as eq 5. gij(s) denotes the relationship between the input uj and the output yi. In the process industry, the FOPDT model (shown as eq 6) is widely used to describe gij(s).35 When the process model is described as a transfer function matrix, the IMC controller designed based on eq 3 is difficult to obtain because the inverse of the model transfer function matrix is too complicated to calculate, especially when multiple time delays exist. Decentralized control methods may simplify the procedure of controller design for multivariable processes. However, the coupling among each loop cannot be handled. Therefore, this work is aimed to provide a novel centralized IMC-PID controller for the multivariable process with multiple time delays.

Gm(s) = E(s)Gm0(s)

(7)

−1 CIMC(s) = Gm0 (s)E−1(s)F(s)

(8)

Γ(s) = diag ((γ1s + 1)(γ2s + 1)...(γns + 1)) Q ij(s) = qij



CPID(s) =



1 ⎞ 1 + ⎟ TIis ⎠ Tfis + ⎝



Q ij(jωk) = (Gm−01(jωk))ij

(10)

⎞ ⎟ 1 ⎟⎠

≈ (Gm−01(s))ij

1 γjs + 1

(14) (15)

{Q ij(jω0)Q ij(jω1)... Q ij(jωr )}

Defining Λ(s) = E−1(s)F(s) = diag(1/(λ1s + 1),1/(λ2s + 1)...1/ (λns + 1)), the filter matrix is written as eq 9.

Gm−01(s)diag ⎜⎜KPi⎜1

Bij s + 1

(13)

To avoid this case, a compensator matrix Γ(s) shown as eq 13 is added. The Nyquist set of the element in the filter matrix is obtained based on the compensated model. The relative degree of the element in the compensated model is zero. Therefore, the element of the approximate filter transfer function matrix can be chosen as eq 14 for simplicity. In eq 14, qij is the steady-state gain of the element in the approximate filter transfer function matrix. To ensure the steady-state decoupling, the parameter qij can be determined as eq 15. Therefore, what we should do is to determine parameters Aij and Bij based on the Nyquist set of the element in the filter matrix.

(9)

−1 C(s) = Gm0 (s)Λ(s)(I − F(s))−1

Aij s + 1

qij = (Γ(s)Gm0(s))−1|s = 0 = Gm−01(s)|s = 0

First, time delays not only deteriorate the performance for control system, but also increase the difficulty in designing controllers. To reduce the influence of time delays, the process model can be rewritten as eq 7, where E(s) = diag(exp(−τ1•s),exp(−τ2•s),...,exp(−τn•s)) and τi• denotes the minimum time delay in the ith row of Gm(s). Substituting eq 7 into eq 3, the IMC controller can be obtained as eq 8. F (s ) = E (s )Λ (s )

(12)

1 1 + jωkγj

k = 1, 2, ..., r , ωk < ωk + 1

(16)

In order to determine parameters Aij and Bij, the Nyquist set of the element in the filter matrix should be obtained first.

(11) C

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Define the Nyquist set of the element in the filter matrix as eq 16: ⎧ Q ij(jω1) Q ij(jωr ) ⎫ ⎪ ⎪ ⋯ ⎬ ⎨1 ⎪ qij qij ⎪ ⎭ ⎩ Q ij(jωk) = (Gm−01(jωk))ij

1 1 + jωkγj

k = 1, 2, ..., r , ωk < ωk + 1

(17)

Generally, the start frequency ω0 is chosen as ω0 = 0. All the frequencies used in eq 16 are selected from the low-frequency section. Therefore, the normalized Nyquist set can be used to determine Aij and Bij. The normalized Nyquist set is shown as eq 17.

Figure 3. Four conceivable cases of the normalized Nyquist dots in eq 17.

⎡ A ⎤ ⎡ a −b ⎤−1⎡ d ⎤ ij ij ⎢ ij ⎥ ⎢ ij ⎥ ⎢ ⎥ ⎢ B ⎥ = ⎢ b c ⎥ ⎢d ⎥ ⎢⎣ ij ⎥⎦ ⎣ ij ij ⎦ ⎣ ij ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪bij = ⎪ ⎪ ⎩

value, which will lead to an unstable transfer function. If the trend of the Nyquist curve is in area II or area III, Bij will be a positive value, which will lead to a stable transfer function. Therefore, the desired γj will make the trend of Qij(ω) in area II or area III. To acquire the value range of γj, the normalized Nyquist curve of every element in the uncompensated filter G−1 m0(ω) should be depicted. Then, according to the normalized Nyquist curve of every element in G−1 m0(ω), the corresponding value range of the compensator parameter can be determined as Ψij. The detailed explanation about the existence and the selection rules for parameter γj can be seen in the Appendix. Here, only a summary −1 is given. G̃ −1 m0(ω) denotes the normalized Nyquist set of Gm0(ω) −1 • If the trend of (G̃ m0(ω))ij is in area I, then the trend of the −1 compensated filter element (G̃ −1 m0(ω))ij Γj (ω) is in area III (with a sufficiently big γj). • If the trend of (G̃ −1 m0(ω))ij is in area II, then the trend of the −1 compensated filter element (G̃ −1 m0(ω))ijΓj (ω) is in area II (with a sufficiently small γj) or area III(with a sufficiently big γj). • If the trend of (G̃ −1 m0(ω))ij is in area III, then the trend of −1 the compensated filter element (G̃ −1 m0(ω))ijΓj (ω) is and is only in area III (with an arbitrary positive γj). • If the trend of (G̃ −1 m0(ω))ij is in area IV, then the trend of the −1 compensated filter element (G̃ −1 m0(ω))ijΓj (ω) is and is only in area III (with a sufficiently big γj). It is worth noting that the jth compensator Γj(ω) has an effect on all the elements in the jth column of matrix (Gm0(ω))−1. Therefore, the compensator Γj(ω) should ensure the normalized trends of all the elements in the jth column of matrix (Γ(ω) Gm0(ω))−1 are in the stable areas. In other words, the desired parameter γj should be selected in the intersection of the n corresponding solution sets as eq 19:

⎧ 2⎤ r ⎡ Q ij(jωi) ⎥ ⎪ ⎢ 2 aij = ∑ ωi ⎪ cij = ∑ ⎢ωi ⎥ qij i=0 ⎪ i=0 ⎢ ⎥⎦ ⎣ ⎨ ⎡ ⎤ ⎛ Q ( jω ) ⎞ ⎪ r i ⎥ ij ⎛ Q ( jω ) ⎞ r ⎟ ⎪ ∑ ⎢⎢ωi2 Re⎜⎜ ⎜ ij i ⎟ ⎟ ⎥ = ω d Im ∑ q ij i ⎪ ⎜ q ⎟ ⎝ ⎠⎦ i=0 ⎣ ij ⎝ ⎠ i=0 ij ⎩ r

2

(18)

From eq 17, it can be seen that the normalized Nyquist set is just a set of complex numbers. To determine the parameters Aij and Bij, a complex curve fitting technique is necessary. Here, Levy’s complex curve fitting method is employed.38 Parameters Aij and Bij are calculated based on eq 18. In eq 18, function Re(x) and function Im(x) are to obtain the real component and the imaginary component of complex number x, respectively. From eq 17 and eq 18, the parameter of the compensator, γj, should be determined before calculating parameters Aij and Bij. As is mentioned above, the compensator is used to compensate the relative degree of the element in Gm0(s). Here, another important function of the compensator should be introduced. Since the stable filter element is preferred, the selected parameter γj should have the compensated Nyquist set leading to stable filter elements. The selection of parameter γj is then discussed. As is shown in eq 14, the stable filter element means that parameter Bij must be greater than zero. Considering that the transfer function in eq 14 is obtained based on the Nyquist set of the element in the filter matrix, the Nyquist curve of the approximate transfer function should coincide with that of the original filter matrix element. Since the steady-state gain qij has nothing to do with the stability of the transfer function, the normalized Nyquist set shown as eq 17 can be used to analyze the stability. Plot the dots in eq 17 in the complex plane, and then Figure 3 will be obtained. Figure 3 displays four conceivable cases of the normalized Nyquist dots in eq 17. From Figure 3, it can be seen that the right-half complex plane is divided into four areas and all the four conceivable curves start from the critical dot, (1,0). The maximum frequency in eq 17, ωr, ought to ensure that the normalized Nyquist curve of each element is just in one area. Further, ωr had better make sure that both the real component and the imaginary component of the normalized Nyquist curve change monotonously. If the trend of the Nyquist curve is in area I or area IV, Bij will be a negative

γj ∈ Ψ1j ∩ Ψ2j ∩ ... ∩ Ψnj

(19) −1 (G̃ −1 m0(ω))ijΓj (ω).

where Ψij is the solution set obtained by From the above discussion, it can be found that a sufficiently large γj can lead to stable filter elements generally. Therefore, the parameter γj always has a value range which can be described as γj > γj̲ , where γj̲ is the low limit of γj. Considering that γj is just designed for compensating the relative degree and obtaining stable filter elements, there is no need to select a too large γj. Therefore, γj can be determined as (1.1−1.2)γj̲ . 2.3. Controller parameter tuning. Based on the above analysis, the feasible approximate filter transfer function matrix Q(s) is determined. Then, the equivalent centralized D

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research PID controller shown in eqs 11 and 12 can be rewritten as eqs 20 and 21: ⎛ ⎛ 1 ⎞ 1 CPID(s) = Gm−01(s)diag ⎜⎜KPi⎜1 + ⎟ TIis ⎠ Tfis + ⎝ ⎝

Theoretically, the perfect decoupling performance means the norm of Hij(jω) remains zero at all frequencies. Hence, the maximum peak of Hij(jω) shown as eq 24 can be used as an index to evaluate the decoupling performance. The maximum peak of Hij(jω) should be restrained for good decoupling performance. From eq 24, the maximum peak of Hij(jω) decreases with λj increasing. Therefore, λj has its minimum value for a given acceptable maximum peak.

⎞ ⎟ 1 ⎟⎠

⎛ ⎛ 1 ⎞ 1 = Q (s)Γ(s)diag ⎜⎜KPi⎜1 + ⎟ TIis ⎠ Tfis + ⎝ ⎝

⎞ ⎟ 1 ⎟⎠

S(s) = [I + GP (s)CPID(s)]−1

⎛ ⎛ 1 ⎞ γis + 1 ⎞⎟ = Q (s)diag ⎜⎜KPi⎜1 + ⎟ TIis ⎠ Tfis + 1 ⎟⎠ ⎝ ⎝ ⎛ ⎛ ⎞ 1 1 = Q (s)diag ⎜⎜K̅Pi⎜1 + + Tdi̅ s⎟ TIi̅ s ⎠ Tfi̅ s + ⎝ ⎝

⎞ ⎟ 1 ⎟⎠

⎧ γ + τi • TIi + γi = i ⎪ K̅Pi = KPi· λi + τi • TIi ⎪ ⎪ ⎪ TIi̅ = TIi + γi = γi + τi • ⎪ ⎨ γ ·τi • TIi = i ⎪ Tdi̅ = γi· + γ γ + τi • T Ii ⎪ i i ⎪ λiτi • ⎪ Tfi̅ = Tfi = ⎪ λi + τi • ⎩

⎡ ⎛ ⎛ 1 ⎞ 1 = ⎢I + E(s)Gm0(s)Q (s)Γ(s)diag ⎜⎜KPi⎜1 + ⎟ ⎢⎣ TIis ⎠ Tfis + ⎝ ⎝ ⎡ ⎛ ⎛ 1 ⎞ 1 ≈ ⎢I + E(s)diag ⎜⎜KPi⎜1 + ⎟ ⎢⎣ T ⎝ ⎝ Iis ⎠ Tfis +

(20)

⎞⎤ ⎟⎥ 1 ⎟⎠⎦⎥

−1

⎞⎤ ⎟⎥ 1 ⎟⎠⎥⎦

−1

(25)

Si(s) =

λiτi • s λi + τi • λiτi • s λi + τi •

From eqs 20 and 21, it is obvious that the parameters which need to be tuned are λ1, λ2, ..., λn. Then, the parameters are determined considering both robustness and performance. H(s) = Gp(s)CIMC(s)

Si(jω) =

= E(s)Gm0(s)Q (s)Γ(s)Λ(s) (22)

i≠j

(1 + )e 1 τi •s

−τi •s

(26)

(αi + 1) + jαiβi

(

1

(αi + 1) + jαiβi + e−jβi 1 − j β

i

)

(27)

The maximum sensitivity (Ms), which is the maximum norm of the sensitivity function,39 is used to determine the controller parameters. For simple analysis, let αi = λi/τi• and βi = ωτi•. Then, eq 26 can be rewritten as eq 27.

In terms of performance, two performance indexes are taken into consideration, namely the response speed and the decoupling performance. The performance is analyzed based on the IMC structure. For the nominal condition (Gp(s) = Gm(s)), the closed-loop transfer function of the IMC structure is shown as eq 22. From eq 22, it is obvious that the response speed will decrease with increasing λi. Hence, for a fast response speed, the parameter λi should be determined as small as possible. Hij(s) = [Gm0(s)Q (s)]ij Λj(s)Γj(s)Ei(s)

τi • λi + τi •

Then, the robustness of the control system will be discussed. Here, the sensitivity function S(s) shown as eq 25 is employed to measure the robustness. From eq 25, the sensitivity function of the system can be approximately regarded as a diagonal matrix. Each element in S(s) represents the sensitivity function for each subsystem. Hence, the robustness can be analyzed based on the elements in S(s), namely eq 26.

(21)

≈ E (s )Λ(s )

+1+

+1

α=

0.06991Ms + 0.5942 Ms − 1.0260

(28)

The relationship between Ms and α is depicted in Figure 4. It can be seen that a small α corresponds to a big Ms. A big Ms

(23)

On the other hand, the effect of loop coupling should be weakened. In fact, the designed full-matrix filter Q(s) plays an important role in restraining the loop coupling. However, based on the procedure of designing the filter Q(s), it can be seen that the filter Q(s) is designed using the low frequency section of the process model. The high frequency characteristic of the process model is overlooked. The neglectful high frequency characteristic may deteriorate the decoupling performance of the filter Q(s). The decoupling performance is reflected by the off-diagonal elements of the closed-loop transfer function matrix H(s), namely eq 23. M p(Hij(jω)) = |[Gm0(jω)Q (jω)]ij |

1 + γj2ω2 1 + λj2ω2

i≠j Figure 4. Relationship between Ms and α.

(24) E

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research denotes a bad robustness. To further make the relationship clear, an expression is fitted based on the original date. The expression is shown as eq 28. Typically, the value range of Ms is 1.2−2.0. Based on eq 28, the controller parameter λ can be determined as long as the desired Ms is given. To sum up, the final controller parameter should be determined after comprehensively considering the aforementioned robust index and performance index. 2.4. Limitation of the proposed method. In the first step of the proposed method, the process model should be divided into two parts. This case may exist that, in the ith row of Gm(s), some elements may have much larger time delays than other elements. Thereupon, some elements in the remaining part Gm0(s) may have large time delays. As a result, the fill-matrix filter Q(s), which is computed based on Gm0(s), will be affected because the effective frequency range of Q(s) will shrink. In this case, for improving decoupling performance, the controller parameter λi will have to increase. A large parameter λi will slow down the output response, which can be regard as deteriorating the control performance. Therefore, the time delays existing in Gm0(s) should have a limitation. K e−Ls gm0 = (29) Ts + 1 gm̃ 0 =

K K = Tars + 1 (T + L )s + 1

Chart 1. Design Flow of the Proposed Method

(30)

Example 1 (Industrial-Scale Polymerization (ISP) Reactor): The ISP reactor is a TITO (two inputs and two outputs) process. Its transfer function matrix is shown as eq 31.

When we compute Q(s) based on Gm0, the approximation errors are inevitable. The approximation errors come from two sources. One is that the elements in Gm0 are not all rational because of the remaining time delays in Gm0(s). The other is that even all elements in Gm0(s) are rational; the errors still exist because a low order transfer function (Qij(s)) is used to describe the dynamic of high order model. The errors coming from the second source can be decreased by complicating the form of filter. Therefore, the limitation for the remaining time delays is mainly used for limiting the errors coming from the first source. If the elements in Gm0(s) are described as eq 29, the approximate rational transfer function is chosen as eq 30, where L is the remaining time delay and Tar is the average residence time, which is computed by −g′m0(0)/gm0(0). It can be found that if L < 0.3T, then the magnitude error and the phase error between the original irrational transfer function and the approximate rational transfer function may be less than 1 dB and 10°, respectively, within the scope of the corner frequency. Therefore, the limitation of the proposed method is that the remaining time delay of the element in Gm0(s) to its corresponding time constant ratio is less than 0.3. Even the control performance of the proposed method may worsen when significantly large remaining time delays exist; the proposed method may be superior to steady-state decoupling. In the worst case that the effective frequency range of Q(s) shrinks into a frequency point, the proposed method shares the same decoupling performance with steady-state decoupling. From the above, the whole design flow of the proposed method is described in Chart 1.

⎡ 22.89e−0.2s −11.64e−0.4s ⎤ ⎢ ⎥ ⎢ 4.572s + 1 1.807s + 1 ⎥ Gp(s) = ⎢ ⎥ −0.2s 5.8e−0.4s ⎥ ⎢ 4.689e ⎣ 2.174s + 1 1.801s + 1 ⎦

(31)

Gm(s) = E(s)Gm0(s) ⎡ e−0.2s =⎢ ⎢⎣

⎡ 22.89 −11.64e−0.2s ⎤ ⎢ ⎤⎢ 4.572s + 1 1.807s + 1 ⎥⎥ ⎥⎢ ⎥ e−0.2s ⎥⎦⎢ 4.689 5.8e−0.2s ⎥ ⎣ 2.174s + 1 1.801s + 1 ⎦

⎡ ⎢ 0.031 × Q (s) = ⎢ ⎢ ⎢− 0.025 × ⎣

(32)

4.5018s + 1 ⎤ ⎥ 2.6341s + 1 ⎥ 4.0672s + 1 2.7472s + 1 ⎥ 0.1222 × ⎥ 1.3822s + 1 3.4504s + 1 ⎦

4.263s + 1 1.399s + 1

0.0621 ×

(33)

To begin with, the full-matrix filter Q(s) should be obtained. The nominal process model can be divided into eq 32. From eq 32, the requirement for the remaining time delays is met. Then, the normalized Nyquist set of G−1 m0(jω) is calculated and depicted in Figure 5. From Figure 5, it can be seen that the normalized trends of (G−1 m0(jω))ij are either in area I or in area II. The value ranges of compensator parameters can be determined as 0 < γ1 < 3.58 or 1.38 < γ2 < 3.79 based on the proposed rules. The parameters are selected as γ1 = 1 and γ2 = 2. Then, based on the Nyquist set of the compensated model and the complex curve fitting technique, the full-matrix filter Q(s) is obtained as eq 33. To verify the curve fitting effect, Figure 6 shows the Nyquist plots of the real models and the approximate models.

3. SIMULATION STUDY In this section, two examples are provided to determine the effectiveness of the proposed method. First, a square process is discussed. F

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 5. Normalized Nyquist set of G−1 m0(jω).

CPID

⎡⎛ ⎞ (0.281s + 0.0659) 1 ⎢ ⎜1 + + 0.167s⎟ ⎠ (1.399s + 1)(0.129s + 1) 1.2s ⎢⎝ =⎢ ⎞ −(0.216s + 0.0531) ⎢⎛ 1 + 0.167s⎟ ⎢ ⎜1 + 1.2s ⎠ (1.382s + 1)(0.129s + 1) ⎣⎝

When the full-matrix filter Q(s) is obtained, the remaining work is to determine the IMC controller parameters, λ1 and λ2. The robust index Ms is confined as Ms ≤ 1.4. Then, based on eq 28, the parameters λ1 and λ2 should satisfy λ1 ≥ 0.3647 and λ2 ≥ 0.3647. Taking the decoupling performance index, Mp, into consideration, if the acceptable maximum peak is given as 15%, the parameters λ1 and λ2 should satisfy λ1 ≥ 0.2315 and λ2 ≥ 1.468. Therefore, the parameters λ1 and λ2 are finally determined as λ1 = 0.3647 and λ2 = 1.468. The centralized PID controller is expressed as eq 34. IAEij =

∫0



|rj(t ) − yj (t )|dt

(37)

⎧ j=i ⎪1 rj(t ) = ⎨ ⎪ ⎩0 j ≠ i

|| ΔO(jω)||
0 j k ⇔⎨ , ⎪ −1 ⎪ Im((Gm̃ 0 (jωk ))ij ) > γj ⎪ −1 ⎪ ωk Re((Gm̃ 0 (jωk))ij ) ⎩ (A.1) k = 1 ... r It is difficult to analyze the conditions that the parameter γj should meet from eq A.1. −1 1.ii If the trend of (G̃ −1 m0(ω))ijΓj (ω) is pulled into area III, then

4. CONCLUSION In this paper, we try to propose a centralized IMC-PID controller for multivariable processes with multiple time delays. To reduce the influence of the time delays, the process model is first divided into the minimum common time delay part and the remaining part. Then, the IMC controller and its equivalent centralized PID controller are obtained based on the IMC principle. The centralized PID is composed of a full-matrix filter and a diagonal PID controller matrix. L

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

that a sufficiently small γj can keep the trend of −1 (G̃ −1 m0(ω))ijΓj (ω) in the stable area (area II). −1 2.ii If the trend of (G̃ −1 m0(ω))ijΓj (ω) is pulled into area III, then

⎞ ⎛ 1 −1 ⎟, 1 > Re⎜⎜(Gm̃ 0 (jωk))ij ⎟ + j 1 γ ω k ⎠ ⎝ j ⎞ ⎛ 1 −1 ⎟ 0 > Im⎜⎜(Gm̃ 0 (jωk))ij 1 + jγjωk ⎟⎠ ⎝ ⎧(Im((G̃ −1(jω )) ) − γ ω ) m0 k ij j k ⎪ ⎪ < (1 − Re((Gm̃ −01(jωk))ij ) ⎪ ⎪ ) ⇔⎨ / γ ω j k ⎪ ⎪ −1 ⎪ Im((Gm̃ 0 (jωk ))ij ) < γj ⎪ ω Re((G̃ −1(jω )) ) k m0 k ij ⎩

⎞ ⎛ 1 −1 ⎟, 1 > Re⎜⎜(Gm̃ 0 (jωk))ij ⎟ 1 j + γ ω k ⎠ ⎝ j ⎞ ⎛ 1 −1 ⎟ 0 > Im⎜⎜(Gm̃ 0 (jωk))ij 1 + jγjωk ⎟⎠ ⎝ ,

⎧(Im((G̃ −1(jω )) ) − γ ω ) m0 k ij j k ⎪ ⎪ < (1 − Re((G̃ −1(jω )) ))/γ ω m0 k ij j k ⇔⎨ ⎪ −1 ⎪ Im((Gm̃ 0 (jωk))ij ) < γjωk Re((Gm̃ −01(jωk))ij ) ⎩

k = 1 ... r

(A.2)

⎛ Im((G̃ −1(jω )) ) ⎞ m0 k ij ⎟ γj > max ⎜⎜ −1 1 ≤ k ≤ r ω Re((G̃ (jω )) ) ⎟ ⎝ k m0 k ij ⎠

−1

⇔ωk(γjωk − Im((Gm̃ 0 (jωk))ij )) 1 −1 , > (Re((Gm̃ 0 (jωk))ij ) − 1) γj

(A.3)

(A.6)

From eq A.2, conditions that the parameter γj should meet can be obtained as eq A.3. Inequality A.3 means that a −1 sufficiently large γj can pull the trend of (G̃ −1 m0(ω))ijΓj (ω) into the stable area (area III). (2) The trend of (G̃ −1 m0(ω))ij is in area II (stable). The Nyquist dot satisfies −1 1 < Re((Gm̃ 0 (jωk))ij ),

From eq A.6, it can be seen that, with the increase of γj, the left part of the inequality will increase while the right part will decrease. Therefore, a sufficiently large γj will satisfy this inequation. (3) The trend of (G̃ −1 m0(ω))ij is in area III (stable). The Nyquist dot satisfies

−1 0 < Im((Gm̃ 0 (jωk))ij )

−1 −1 0 < Re((Gm̃ 0 (jωk))ij ) < 1, Im((Gm̃ 0 (jωk))ij ) < 0

k = 1 ... r

k = 1 ... r

The compensator Γ j (ω) has to keep the trend of −1 (G̃ m0 (ω)) ij Γ j− 1 (ω) in area II or pull the trend of −1 (G̃ m0(ω))ijΓ−1 j (ω) into area III. −1 2.i If the trend of (G̃ −1 m0(ω))ijΓj (ω) is kept in area II, then

The compensator Γ j (ω) has to keep the trend of −1 (G̃ m0 (ω)) ij Γ j−1 (ω) in area III or pull the trend of −1 ̃ (Gm0(ω))ijΓ−1 j (ω) into area II. −1 3.i If the trend of (G̃ −1 m0(ω))ijΓj (ω) is pulled into area II, then

⎛ ⎞ 1 −1 ⎟, 1 < Re⎜⎜(Gm̃ 0 (jωk))ij ⎟ + γ ω 1 j ⎝ j k⎠

⎛ ⎞ 1 −1 ⎟ > 1, Re⎜⎜(Gm̃ 0 (jωk))ij 1 + jγjωk ⎟⎠ ⎝

⎛ ⎞ 1 −1 ⎟ 0 < Im⎜⎜(Gm̃ 0 (jωk))ij 1 + jγjωk ⎟⎠ ⎝ ⎧ ̃ −1 ⎪ (Re((Gm0 (jωk))ij ) − 1) ⎪ + γ ωk(Im((Gm̃ −01(jωk ))ij ) j ⎪ ⎪ − γω ) > 0 j k ⇔⎨ ⎪ ⎪ Im((Gm̃ −01(jωk ))ij ) > γj ⎪ −1 ⎪ ωk Re((Gm̃ 0 (jωk))ij ) ⎩ k = 1 ... r

⎛ Im((G̃ −1(jω )) ) ⎞ m0 k ij ⎟ γj < min ⎜⎜ −1 1 ≤ k ≤ r ω Re((G̃ (jω )) ) ⎟ ⎝ k m0 k ij ⎠

k = 1 ... r

⎛ ⎞ 1 −1 ⎟>0 Im⎜⎜(Gm̃ 0 (jωk ))ij 1 + jγjωk ⎟⎠ ⎝ ⎧(Re((G̃ −1(jω )) ) − 1) + γ m0 k ij j ⎪ ⎪ ωk(Im((Gm̃ −01(jωk ))ij ) − γ ωk) > 0 ⎪ j ⇔⎨ , ⎪ Im((G̃ −1(jω )) ) > γ m0 k ij j ⎪ ⎪ ωk Re((Gm̃ −01(jωk))ij ) ⎩

,

(A.4)

k = 1 ... r

(A.7)

No positive γj can satisfy the above inequation. In other words, no compensator can pull the trend of −1 (G̃ −1 m0(ω))ijΓj (ω) into area II. −1 3.ii If the trend of (G̃ −1 m0(ω))ijΓj (ω) is kept in area III, then

(A.5)

From eq A.4, conditions that the parameter γj should meet can be obtained as inequality A.5. Inequality A.5 means M

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research ⎞ ⎛ 1 −1 ⎟, 1 > Re⎜⎜(Gm̃ 0 (jωk ))ij 1 + jγjωk ⎟⎠ ⎝

⎞ ⎛ 1 −1 ⎟, 1 > Re⎜⎜(Gm̃ 0 (jωk ))ij ⎟ 1 j + γ ω k ⎠ ⎝ j ⎞ ⎛ 1 1 − ⎟ 0 > Im⎜⎜(Gm̃ 0 (jωk ))ij 1 + jγjωk ⎟⎠ ⎝ ⎧(Im((G̃ −1(jω )) ) − γ ω ) m0 k ij j k ⎪ ⎪ < (1 − Re((Gm̃ −01(jωk ))ij )) ⎪ ⎪ ⇔ ⎨ /γjωk ⎪ ⎪ Im((Gm̃ −01(jωk ))ij ) < γj ⎪ ⎪ ωk Re((Gm̃ −01(jωk ))ij ) ⎩

⎞ ⎛ 1 −1 ⎟ 0 > Im⎜⎜(Gm̃ 0 (jωk ))ij 1 + jγjωk ⎟⎠ ⎝ ⎧(Im((G̃ −1(jω )) ) − γ ω ) m0 k ij j k ⎪ ⎪ < (1 − Re((Gm̃ −01(jωk ))ij ))/γ ωk ⎪ j , ⇔⎨ 1 − ⎪ Im((Gm̃ 0 (jωk ))ij ) ⎪ −1 ⎪ < γjωk Re((Gm̃ 0 (jωk ))ij ) ⎩

,

k = 1 ... r

−1

k = 1 ... r

⇐(1 − Re((Gm̃ 0 (jωk ))ij )) , −1 > γ ωk Im((Gm̃ 0 (jωk ))ij )

(A.8)

k = 1 ... r

j

An arbitrary positive γj can satisfy the above inequality.

(A.10)

−1 Therefore, the trend of (G̃ −1 m0(ω))ijΓj (ω) can but has to keep in the stable area (area III) no matter what the value of γj is. (4) The trend of (G̃ −1 m0(ω))ij is in areaIV (unstable). The Nyquist dot satisfies

⎛ 1 − Re((G̃ −1(jω )) ) ⎞ m0 k ij ⎟ γj > max ⎜⎜ −1 1 ≤ k ≤ r ω Im((G̃ (jω )) ) ⎟ ⎝ k m0 k ij ⎠



1 < Re((Gm−01(jωk))ij ), Im((Gm−01(jωk))ij ) < 0 k = 1 ... r

(A.11)

From eq A.10, the conditions that the parameter γj should meet can be obtained as eq A.11. It is important to note that the solution set obtained by eq A.11 is the subset of the original inequality shown in eq A.10.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

The compensator Γ j (ω) has to pull the trend of −1 (G̃ −1 m0(ω))ijΓj (ω) into area II or area III.

Xinghan Du: 0000-0001-5848-8131 Notes

−1 4.i If the trend of (G̃ −1 m0(ω))ijΓj (ω) is pulled into area II, then

The authors declare no competing financial interest.



⎞ ⎛ 1 −1 ⎟ > 1, Re⎜⎜(Gm̃ 0 (jωk))ij ⎟ + 1 γ ω j ⎝ j k⎠

ACKNOWLEDGMENTS This work is supported by the National Natural Science Foundation of China (61673004) and the Fundamental Research Funds for the Central Universities (Grant ZY1619).



⎛ ⎞ 1 −1 ⎟>0 Im⎜⎜(Gm̃ 0 (jωk ))ij 1 + jγjωk ⎟⎠ ⎝ ⎧ (Re((G̃ −1(jω )) ) − 1) m0 k ij ⎪ ⎪ + γjωk(Im((Gm̃ −01(jωk ))ij ) ⎪ ⎪ − γjωk) > 0 ⇔⎨ ⎪ Im((Gm̃ −01(jωk ))ij ) > γ j ⎪ −1 ωk Re((Gm̃ 0 (jωk))ij ) ⎪ ⎪ > γjωk ⎩ k = 1 ... r

REFERENCES

(1) Ziegler, J. G.; Nichols, N. B. Optimum Settings for Automatic Controllers. Trans. ASME 1942, 64, 759. (2) Skogestad, S. Simple analytic rules for model reduction and PID controller tuning. J. Process Control 2003, 13, 291. (3) Bristol, E. On a new measure of interaction for multivariable process control. Auto. Control, IEEE Trans. 2003, 11, 133. (4) Hansen, J. E.; Heath, J.; Jorgensen, S. B. Control structure selection for energy integrated distillation column. J. Process Control 1998, 8, 185. (5) Witcher, M. F.; Mc Avoy, T. J. Interacting control systems: Steady state and dynamic measurement of interaction. ISA Trans. 1977, 16, 35. (6) Mc Avoy, T.; Arkun, Y.; Chen, R.; Robinson, D.; Schnelle, P. D. A new approach to defining a dynamic relative gain. Control Eng. Pract. 2003, 11, 907. (7) Xiong, Q.; Cai, W. J.; He, M. J. A practical loop pairing criterion for multivariable processes. J. Process Control 2005, 15, 741. (8) He, M. J.; Cai, W. J.; Wei, N.; Xie, L. H. RNGA based control system configuration for multivariable processes. J. Process Control 2009, 19, 1036. (9) Huang, H. P.; Jeng, J. C.; Chiang, C. H.; Pan, W. A direct method for multi-loop PI/PID controller design. J. Process Control 2003, 13, 769. (10) Vu, T. N. L.; Lee, M. Independent design of multi-loop PI/PID controllers for interacting multivariable processes. J. Process Control 2010, 20, 922.

,

(A.9)

There is still no positive γj satisfying the above inequality. −1 The trend of (G̃ −1 m0(ω))ijΓj (ω) cannot be pulled into area

II. −1 4.ii If the trend of (G̃ −1 m0(ω))ijΓj (ω) is pulled into area III,

then N

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (11) Cui, H.; Jacobsen, E. W. Performance limitations in decentralized control. J. Process Control 2002, 12, 485. (12) Goodwin, G. C.; Salgado, M. E. Time-domain performance limitations arising from decentralized architectures and their relationship to the RGA. Int. J. Control 2005, 78, 1045. (13) Gagnon, E.; Pomerleau, A.; Desbiens, A. Simplified, ideal or inverted decoupling? ISA Trans. 1998, 37, 265. (14) Jevtović, B. T.; Mataušek, M. R. PID controller design of TITO system based on ideal decoupler. J. Process Control 2010, 20, 869. (15) Garrido, J.; Vázquez, F.; Morilla, F. Centralized multivariable control by simplified decoupling. J. Process Control 2012, 22, 1044. (16) Rajapandiyan, C.; Chidambaram, M. Controller Design for MIMO Processes Based on Simple Decoupled Equivalent Transfer Functions and Simplified Decoupler. Ind. Eng. Chem. Res. 2012, 51, 12398. (17) Wade, H. L. Inverted decoupling: a neglected technique. ISA Trans. 1997, 36, 3. (18) Garrido, J.; Vázquez, F.; Morilla, F. An extended approach of inverted decoupling. J. Process Control 2011, 21, 55. (19) Garrido, J.; Vázquez, F.; Morilla, F. Inverted decoupling internal model control for square stable multivariable time delay systems. J. Process Control 2014, 24, 1710. (20) Åström, K. J.; Johansson, K. H.; Wang, Q. G. Design of decoupled PI controllers for two-by-two systems. IEE Proc.: Control Theory Appl. 2002, 149, 74. (21) Zhang, W. D.; Sun, Y. X.; Xu, X. M. Two Degree-of-Freedom Smith Predictor for Processes with Time Delay. Automatica. 1998, 34, 1279. (22) Santos, T. L. M.; Flesch, R. C. C.; Normey-Rico, J. E. On the filtered Smith predictor for MIMO processes with multiple time delays. J. Process Control 2014, 24, 383. (23) Ammathil, R.; Narsaiah, T. B.; Rao, A. S. Design of decentralised Smith predictor for multivariable non-square processes with multiple time delays. Int. J. Model. Ident. Control. 2014, 21, 147. (24) Garcia, C. E.; Morari, M. Internal model control. A unifying review and some new results. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 308. (25) Rivera, D. E.; Morari, M.; Skogestad, S. Internal model control: PID controller design. Ind. Eng. Chem. Res. 1986, 25, 2163. (26) Garcia, C. E.; Morari, M. Internal model control. 2. Design procedure for multivariable systems. Ind. Eng. Chem. Process Des. Dev. 1985, 24, 472. (27) Economou, C. G.; Morari, M. Internal Model Control: 6. multiloop design. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 411. (28) Jin, Q. B.; Hao, F.; Wang, Q. IMC-PID Controller Based on Diagonal Equivalent Model for Multivariable Systems with Multiple Time Delays. J. Chem. Eng. Jpn. 2013, 46, 209. (29) Jin, Q. B.; Hao, F.; Wang, Q. A multivariable IMC-PID method for non-square large time delay systems using NPSO algorithm. J. Process Control 2013, 23, 649. (30) Jin, Q. B.; Liang, Z.; Feng, H.; Liu, S. W. Design of a multivariable internal model controller based on singular value decomposition. Can. J. Chem. Eng. 2013, 91, 1103. (31) Jin, Q. B.; Du, X. H.; Wang, Q.; Liu, L. L. Analytical design 2 DOF IMC control based on inverted decoupling for non-square systems with time delay. Can. J. Chem. Eng. 2016, 94, 1354. (32) Liu, T.; Zhang, W. D.; Gu, D. Analytical Design of Decoupling Internal Model Control (IMC) Scheme for Two-Input−Two-Output (TITO) Processes with Time Delays. Ind. Eng. Chem. Res. 2006, 45, 3149. (33) Liu, T.; Zhang, W. D.; Gao, F. R. Analytical decoupling control strategy using a unity feedback control structure for MIMO processes with time delays. J. Process Control 2007, 17, 173. (34) Wang, Q. G.; Zhang, Y.; Chiu, M. S. Decoupling internal model control for multivariable systems with multiple time delays. Chem. Eng. Sci. 2002, 57, 115. (35) Vijay Kumar, V.; Rao, V. S.; Chidambaram, M. Centralized PI controllers for interacting multivariable processes by synthesis method. ISA Trans. 2012, 51, 400.

(36) Skogestad, S. Simple analytic rules for model reduction and PID controller tuning. J. Process Control 2003, 13, 291. (37) Jin, Q. B.; Liu, Q. Decoupling Proportional−Integral−Derivative Controller Design for Multivariable Processes with Time Delays. Ind. Eng. Chem. Res. 2014, 53, 765. (38) Levy, E. C. Complex-curve fitting. IRE Trans. Autom. Control 1959, AC-4, 37. (39) Skogestad, S.; Postlethwaite, I. Multivariable feedback control: Analysis and design, 2nd ed.; Wiley: Chichester, U.K., 2005. (40) Vu, T. N. L.; Lee, M. Multi-loop PI controller design based on the direct synthesis for interacting multi-time delay processes. ISA Trans. 2010, 49, 79. (41) Jin, Q.; Zhu, L.; Wang, Q.; Jiang, B. PI controller design for a TITO system based on delay compensated structure and direct synthesis. Can. J. Chem. Eng. 2016, 94, 1740. (42) Rao, A. S.; Chidambaram, M. Smith delay compensator for multivariable non-square systems with multiple time delays. Comput. Chem. Eng. 2006, 30, 1243.

O

DOI: 10.1021/acs.iecr.6b05011 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX