Article pubs.acs.org/IECR
Adjoint Transfer Matrix Based Decoupling Control for Multivariable Processes Yuling Shen,*,† Youxian Sun,† and Shaoyuan Li‡ †
State Key Laboratory of Industrial Control Technology, Zhejiang University, Hangzhou 310027, China Department of Automation, Shanghai Jiao Tong University, and Key Laboratory of System Control and Information Processing, Ministry of Education of China, Shanghai 200240, China
‡
ABSTRACT: In this paper, a novel decoupling control scheme based on adjoint matrix is proposed. By introducing the concept of characteristic sequence, the characteristic sequence of both the adjoint transfer matrix and the determinant transfer function are derived from that of original process transfer matrix. The adjoint transfer matrix and determinant transfer function are then determined. Finally, the adjoint matrix is selected as a decoupler, and the decentralized controller is designed for the determinant transfer function. The effectiveness of the proposed design approach is verified by three multivariable industrial processes, which shows that it results in better overall performance compared with other methods.
1. INTRODUCTION A general decoupling control system is depicted as in Figure 1, where G(s) is an n × n process transfer function, GI(s) is the
it cannot be achieved in physics. Therefore, the decoupler is further modified as GI (s) = G−1(s)GR (s)
where GR(s) is a user-specified diagonal transfer function matrix. GR(s) has to be well chosen so that GI(s) is physically realized. As the decoupler is also designed based on model inverse, it may result very complicated structure. Besides, when the system dimension is high and nonminimum phase elements are involved in the transfer function elements, it is difficult for the decoupler to be both causal and stable. Recently, Cai et al. proposed the equivalent transfer function (ETF) based decoupling control methods,6,7 in which the decoupler is designed as
Figure 1. Block diagram of decoupler plus single-loop controller.
decoupler, and GC(s) is the decentralized controller, respectively. The resulted transfer function matrix after decoupled is expressed as GR (s) = G(s)GI (s)
T GI (s) = Ĝ (s)GR (s)
Under perfect decoupling control, GR(s) is a stable, proper, and casual diagonal transfer function matrix. In order to decouple the process into multiple independent single-input single-output (SISO) processes, there are lots of options for the decoupler GI(s). In theory, the inverse transfer matrix, G−1(s), can be used as the decoupler, because
adj G(s) = G−1(s)det G(s)
(2)
Two kinds of decoupler have developed from the inverse matrix: static decoupler1,2 and dynamic decoupler.3−5 The former one is the inverse of the process static gain transfer matrix. It is simple to design and implement, but it may not always provide satisfactory closed loop performances if there exists big difference in dynamic characteristics among the transfer function elements. Compared with static decoupler, the dynamic one supposes to provide better performance. However, the form of G−1(s) is so complicated that it is not fit for practical application. Worst of all, © 2012 American Chemical Society
(4)
where Ĝ T(s) is the closed loop equivalent transfer function matrix of G(s). The elements of Ĝ T(s) are in the form of firstorder plus dead time (FOPDT) or second-order plus dead time (SOPDT). Though these ETF based methods are very simple in design, the resulting ETF matrix does not absolutely match with the inverse matrix of G(s). That way, the single loop controllers may fail to meet the performance target. Another possible inverse matrix based decoupler is the adjoint matrix of G(s),8 which is calculated as
(1)
G(s)G−1(s) = I
(3)
(5)
It is noticed that the adjoint matrix can be realizable physically. Especially for TITO processes, the adjoint matrix of G(s) can be written out easily without any computation involved. When adj Received: Revised: Accepted: Published: 16419
August 1, 2012 November 6, 2012 November 11, 2012 November 13, 2012 dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research
Article
G(s) acts on the original mulativariable process, G(s) is fully decoupled as GR (s) = det G(s) ·I
for the transfer functions resulting from the four fundamental operations. As for the transfer function matrix G(s), which is depicted as
(6)
⎡ Dk [g (s)] 11 ⎢ ⎢ D [g (s)] k 21 Dk [G(s)] = ⎢ ⎢⋮ ⎢ ⎢⎣ Dk [gn1(s)]
Dk [g12(s)] ··· Dk [g1n(s)]⎤ ⎥ Dk [g22(s)] ··· Dk [g2n(s)]⎥ ⎥ ⎥ ⋮ ⋱ ⎥ Dk [gn2(s)] ··· Dk [gnn(s)]⎥⎦ (11)
for k = 1, 2, ..., l, where G(s) is the transfer function matrix of an n × n process; gij(s),i,j = 1, 2, ..., n, are the transfer function elements of G(s); and l is the length of characteristic matrix sequence. Remark 1. By definition, it is easy to conclude that both Dk and Ek are even functions, i.e.
(7)
for k = 1, 2, ..., l, and l is the length of characteristic sequence, which can be further expanded as Dk [g (s)] = fk (E1[g (s)], E2[g (s)], ..., Ek [g (s)])
(10)
the characteristic matrix sequence is calculated as
2. CHARACTERISTIC MATRIX SEQUENCE OF TRANSFER FUNCTION MATRIX Define the characteristic sequence of transfer function as Dk [g (s)] = [log(g (s))](k) |s = 0
g12(s) ··· g1n(s) ⎤ ⎥ g22(s) ··· g2n(s)⎥ ⎥ ⎥ ⋮ ⋱ ⎥ gn2(s) ··· gnn(s)⎥⎦
⎡ g (s ) ⎢ 11 ⎢ g (s ) G(s) = ⎢ 21 ⎢⋮ ⎢ ⎢⎣ gn1(s)
Thus, it is reasonable to select adj G(s) as the decoupler. However, how to determine the model of adj G(s) for highdimensional multivariable processes and how to design controllers for GR(s) are still problems to be solved. In this paper, a novel adjoint transfer matrix based decoupling control scheme for multivariable processes is presented. The design is based on the concept of characteristic (matrix) sequence, which is used to determine the model of adjoint transfer matrix and determinant transfer fucntion. Then, the decoupler is designed for the system based on the adjoint transfer matrix and proportional−integral (PI) /proportional−integral− derivative (PID) controllers are designed individually to meet the control system performance.
Dk [g (s)] = Dk [−g (s)]
(8)
(12)
and
for k = 1, 2, ..., l, where
Ek [g (s)] = Ek [−g (s)]
(13)
(k)
Ek [g (s)] =
g (s ) g (s )
3. PARAMETRIZATION OF THE DETERMINANT OF A TRANSFER FUNCTION MATRIX Consider an open loop stable n × n multivariable system, G(s) is the process transfer function matrix, which is described by
(9)
s=0
and more specific expressions about f k(•) are given in the Appendix. The results of Dk[g(s)] for several commonly used models are listed in Table 1. It is seen that the characteristic
⎡ g (s ) ⎢ 11 ⎢ g (s ) G(s) = ⎢ 21 ⎢ ... ⎢ g (s ) ⎣ n1
Table 1. Dk[g(s)], k = 1, 2 , 3, for FOPDT and SOPDT Models model ke−θs/(τs + 1) ke−θs/(τs + 1)2 ke−θs/(a2s2 + a1s + 1) ke−θs/[(τ3s + 1)/ (τ1s + 1)(τ2s + 1)] ke−θs(b1s + 1)/ (a2s2 + a1s + 1)
D1 = E1
D2 = E2 − E12
−τ − θ −2τ − θ −a1 − θ −τ1 − τ2 + τ3 − θ
τ2 2τ2 −2a2 + a12 τ12 + τ22 − τ32
−a1 + b1 − θ
a12 − 2a2 − b12
D3 = E3 − 3E1E2 + 2E13 −2τ3 −4τ3 6a1a2 − 2a13 −2τ13 − 2τ23 + 2τ33 −2a13 + 6a1a2 + 2b13
g12(s) ... g1n(s) ⎤ ⎥ g22(s) ... g2n(s)⎥ ⎥ ... ... ... ⎥ gn2(s) ... gnn(s)⎥⎦
(14)
where gij(s), i,j = 1, 2, ..., n, are the transfer function elements of G(s). Arithmetically, the determinant of transfer function matrix can be expanded as n
det(G(s)) =
∑ gij(s)·(−1)i+ j det Gij(s)
(15)
j=1
for i = 1, 2, ..., n, where Gij(s) is the resulting matrix after the elements in the ith row and jth column are deleted from G(s). By substituting eq 15 into eq 9, we have that
sequence carries the process dynamic information, such as time constant and average residence time, etc. Moreover, Table 2 also provides the operational rules of Dk[•] and Ek[•]
n
Ek [det(G(s))] =
Table 2. Useful Operational Rules for Dk[•] and Ek[•] g(s)
Dk[g(s)]
∑ j = 1 kijK ij·(Ek [gij(s) ·( −1)i + j det Gij(s)]) n
∑ j = 1 kij·K ij (16)
Ek[g(s)]
Ek[g(s)] a·g1(s) (a ≠ 0) Dk[g(s)] g1(s) ± g2(s) Dk[g1(s) ± g2(s)] {g1(0)Ek[gi(s)] ± g2(0)Ek[g2(s)]}/{g1(0) ± g2(0)} g1(s)·g2(s) Dk[g1(s)] + ∑ki=0CikEi[g1(s)]Ek−i[g2(s)] Dk[g2(s)] Dk[g1(s)] − ∑ki=0CikEi[g1(s)]Ek−i[1/g2(s)] g1(s)/g2(s) Dk[g2(s)]
for any i and k = 1, 2, ..., l, where Kij denotes the complement minor corresponding to each element kij of K. Since Ek is even function, eq 16 can be further reduced as n
Ek [det(G(s))] =
16420
∑ j = 1 kijEk [gij(s)det Gij(s)]·K ij n
∑ j = 1 kij·K ij
=
det K̃ k det K
(17)
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research
Article
Remark 2. In particular, when G(s) is a 2 × 2 transfer function matrix, eqs 16 and 17 are reduced into
where ⎡ k11 k12 ⎢ ⋮ ⋮ ⎢ ⎢ ki2Ek [gi2(s) ⎢ k E [g (s) K̃ k = ⎢ i1 k i1 i1 det G ( s )] det Gi2(s)] ⎢ ⎢ ⋮ ⋮ ⎢ ⎢⎣ kn1 k n2
··· ··· ··· ··· ···
⎤ ⎥ ⋮ ⎥ ⎥ kinEk [gin(s) ⎥ ⎥ det Gin(s)]⎥ ⎥ ⋮ ⎥ ⎥⎦ knn k 1n
Ek [det(G(s))] =
⎡ k11Ek [g (s)g (s)] k12Ek [g (s)g (s)]⎤ 11 22 12 21 ⎥ K̃ = ⎢ ⎥⎦ ⎢⎣ k 21 k 22
(18)
det G(s) = kdet G·
(19)
bdet G , ms m + bdet G , m − 1s m − 1 + ··· + bdet G ,1s + 1
for k = 1, 2, ..., l, where the specific expression of f k(•) are also found in the Appendix. For simplicity, specifiy the form of det G(s) as a 2s 2 + a1s + 1
adet G , ns n + adet G , n − 1s n − 1 + ··· + adet G ,1s + 1
where kdet G, θdet G, adet G,i, i = 1, 2, ···, n, and bdet G,i, j = 1, 2, ..., m are the model parameters. Since kdet G and θdet G are given by eq 21, l is equal to m + n, such that m + n equations like eq 22 are supposed to be established to determine the parameters in the numerator and denominator polynomial. In other words, the model form is specified at the first beginning of model parametrization.
(20)
where kdet G, θdet G, a1, a2, and b1 are static gain, dead time, and other parameters of det G(s), repectively. First, the steady state gain and the delay time are calculated by
4. DECOUPLER AND PID CONTROLLER DESIGN Since
(21)
where θij and θdet Gij are the dead time of gij(s) and det Gij(s). Second, the rest of the parameters are determined according to Table 1, by solving the following equation set: ⎧−a + b − θ = D det G 1 1 ⎪ 1 ⎪ 2 2 ⎨ a1 − 2a 2 − b1 = D2det G ⎪ ⎪−2a 3 + 6a a + 2b 3 = D det G ⎩ 1 1 2 1 3
G(s)G−1(s) = I
(27)
and G−1(s) =
adj G(s) det G(s)
(28)
the following relationship is obtained as G(s) ·adj G(s) = det G(s) ·I
(22)
··· ··· ⋱ ···
GI (s) = adj G(s) = [( −1)i + j det Gij(s)]ji
(30)
and the resulted transfer function matrix is derived as GR (s) = det G(s) ·I
(31)
Thus, an ideal decoupling control system based on adjoint matrix is constructed as in Figure 2, where G(s) is an n × n multivariable process, adj G(s) acts as the decoupler, and GC(s) is the decentralized controller designed for GR(s), respectively. To speed up the response of forward transfer function, the decoupler is chosen as G̅I (s) = adj G(s) ·GD(s)
( −1) det G ( ⎤⎥ s) ⎥ ⎥ n+2 n2 ( −1) det G (⎥ ⎥ s) ⎥ ⎥ ⋮ ⎥ n+n nn (− 1) det G ( ⎥ ⎥⎦ s) n+1
(29)
Note that right side of eq 29 is a diagonal matrix. Therefore, the decoupler can be designed as
From the above, the procedure of parametrization for det G(s) is summarized as follows: Step 1: Set L = 1, where L is a cycle index. Step 2: On the basis of eqs 17 and 19, calculate Ek and Dk for all the L-order cofactors of G(s), for k = 1, 2, ..., l. Then, reset L with L + 1. Step 3: Repeat step 2 until L is larger than n. Step 4: Solve eqs 21 and 22, the model of det G(s) is finally determined. Simultaneously, Ek and Dk of det Gij(s) are calculated by following the above procedure. By solving eqs 21 and 22, the models of det Gij(s), i,j = 1, 2, ..., n are gotten. Therefore, the model of the conjugate transpose matrix of G(s) is determined as ⎡ det G11(s) −det G21(s) ⎢ ⎢ ⎢ ⎢ − det G12(s) det G22(s) ⎢ adj G(s) = ⎢ ⎢ ⋮ ⋮ ⎢ n+1 1n n+2 ⎢(− 1) det G (s) (− 1) det G2n( ⎢⎣ s)
e−θdet Gs (26)
e−θdet Gs
⎧ ⎪ k det G = det K ⎨ ⎪θ ⎩ det G = min{θij + θdet Gij , j = 1, 2, ..., n}
(25)
Remark 3. It is noticed that the value of l depends on the number of model parameters. Without loss of generality, det G(s) takes the form of
Dkdet G = fk (E1[det G(s)], E2[det G(s)], ..., Ek [det G(s)])
det G(s) = kdet G·
(24)
where
and Ek[gij(s)det Gij(s)], j = 1, 2, ..., n are calculated according to Table 2. Accordingly, the characteristic sequence of det G(s) is G denoted as Ddet , which can be calculated by k
b1s + 1
det K̃ det K
(32)
n1
where GD(s) = diag{edis , i = 1, 2, ..., n}
(33)
and di = min{θdet Gij , j = 1, 2, ..., n}
(34)
Accordingly, the decoupled transfer function matrix is modified as
(23) 16421
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research
Article
D1det G = −28.2362, D2det G = 318.2133, D3det G = −443.8203
Using eqs 21 and 22, both FOPDT and SOPDT models of det G(s) are obtained as det G(s)1 =
−123.6 −4s e 24.24s + 1
and Figure 2. Block diagram of adjoint matrix based decoupling control scheme.
det G(s)2 =
−123.6 e − 4s 134.6s + 24.24s + 1 2
The adjoint matrix of G(s) is directly written out as dis
G̅R (s) = diag{det G(s)e , i = 1, 2, ..., n}
⎡ −19.4e−3s 18.9e−3s ⎤ ⎢ ⎥ ⎢ 14.4s + 1 21s + 1 ⎥ adj G(s) = ⎢ ⎥ − 7s 12.8e−s ⎥ ⎢ −6.6e ⎣ 10.9s + 1 16.7s + 1 ⎦
(35)
In particular, G̅ I(s) for a TITO process is written as ⎡ g (s)ed1s g (s)ed2s ⎤ 12 ⎢ 22 ⎥ G̅I (s) = ⎢ d1s d 2s ⎥ g11(s)e ⎦ ⎣ g21(s)e
(36)
According to eq 32, the decoupler is selected as
where d1 = min{θ21,θ22} and d2 = min{θ11,θ12}. Following the procedure given in section 3, the models of det Gij(s) and det G(s) can be determined based on the characteristic sequence. Then, GC(s) can be designed independently for the elements of G̅ R(s). The existing SISO design methods are available to guarantee the stability and satisfactory performance. The standard PID control scheme is adopted for control system design,9−11 which is in the form of gc , ii(s) = k p, ii +
k i, ii s
+ kd, iis
⎡ −19.4 18.9e−2s ⎤ ⎢ ⎥ ⎢ 14.4s + 1 21s + 1 ⎥ GI (s) = adj G(s) ·D = ⎢ ⎥ − 4s 12.8 ⎥ ⎢ −6.6e ⎣ 10.9s + 1 16.7s + 1 ⎦
where D = diag{e3s, es}. Therefore, the resulted matrixes after decoupled are obtained as
(37)
⎡ − 123.6 −s ⎤ e 0 ⎢ ⎥ 24.24 s + 1 ⎥ GR (s)1 = det G(s)1 · D = ⎢ ⎢ − 123.6 −3s ⎥ e ⎥ ⎢⎣ 0 ⎦ 24.24s + 1
where kp,ii, ki,ii, and kd,ii are the PID controller parameters.
5. SIMULATION EXAMPLES Example 1. Consider the Wood−Berry Process12
and
⎡ 12.8e−s −18.9e−3s ⎤ ⎢ ⎥ ⎢ 16.7s + 1 21s + 1 ⎥ G (s ) = ⎢ ⎥ − 7s −19.4e−3s ⎥ ⎢ 6.6e ⎣ 10.9s + 1 14.4s + 1 ⎦
GR (s)2 = det G(s)2 ·D =
According to Table 1, the open loop characteristic matrixes are calculated as ⎡−17.7 −24.0 ⎤ open ⎡ 278.89 441 ⎤ D1open = ⎢ =⎢ ⎥, ⎥, D ⎣−17.9 −17.4 ⎦ 2 ⎣118.81 207.36 ⎦ ⎡−9314.926 − 18522 ⎤ D3open = ⎢ ⎥ ⎣ 2590.058 5971.968 ⎦
⎡ ⎤ −123.6 e −s 0 ⎢ ⎥ 2 ⎢ 134.6s + 24.24s + 1 ⎥ ⎢ ⎥ −123.6 ⎢0 ⎥ 2 ⎣ 134.6s + 24.24s + ⎦1
The PI/PID controllers are designed by method GMP as, respectively ⎤ ⎡ 0.003177 0 ⎥ ⎢−0.07701 − s ⎥ ⎢ GC (s)1 = ⎢ ⎥ −0.02567 ⎥ ⎢0 0.001059 ⎥ ⎢ − ⎦ ⎣ s
Accordingly, Eopen k , k = 1, 2, 3 are obtained as ⎡−17.7 −24.0 ⎤ open ⎡ 592.18 1017.0 ⎤ =⎢ E1open = ⎢ ⎥, ⎥, E ⎣−17.9 −17.4 ⎦ 2 ⎣ 439.22 510.12 ⎦ ⎡−29669.218 −64098.0 ⎤ E3open = ⎢ ⎥ ⎣−14705.494 − 22064.184 ⎦
and ⎡ ⎤ 0.002542 0 ⎢− 0.06187 − ⎥ s ⎢ ⎥ − s 0.3421 GC (s)2 = ⎢ ⎥ ⎢ ⎥ 0.00085 − 0.02062 − − 0.1141s ⎥ ⎢0 ⎣ ⎦ s
G On the basis of eq 17, Edet , k = 1, 2, 3 are calculated as k
E1det G = − 28.2362, E2det G = 1115.4947, E3det G = − 53905.7533 G From Table 1, Ddet , k = 1, 2, 3 are calculated as k
16422
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research
Article
The unit-step set-point response and unit-step disturbance input response are given in Figures 3 and 4, together with the result of
Figure 4. Load rejecting responses for WB (2 × 2) process.
det G D1 ij
⎡−120.2949 −270.8531 −172.3362 ⎤ ⎥ ⎢ = ⎢−200.9905 −172.5300 115.7305 ⎥ ⎢⎣ 534.8183 −257.8967 −264.6601 ⎥⎦
det G D2 ij
⎡ 3520.8865 ⎤ 5331.4169 5433.6861 ⎢ ⎥ = ⎢13017.0522 10549.5805 − 108327.5723⎥ ⎢⎣−746470.2649 15354.8570 19040.5687 ⎥⎦
det G D3 ij
⎡− 372687.05448 − 2151063.0873 − 386134.1118 ⎤ ⎥ ⎢ = ⎢− 2931808.6212 − 1914476.6162 87424314.3160 ⎥ ⎢⎣1357666641.5361 − 2203570.5648 − 3358627.4392 ⎥⎦
det Gij
D4
det G D5 ij
Figure 3. Nominal system responses for WB (2 × 2) process.
the decoupling method proposed by Morrilla et al. It is seen that the proposed method has much better performance. At the same time, the control performance also depends on the model accuracy of det G(s). For WB process, it is clear that the SOPDT model gives a more complete description about the dynamic characteristics than the FOPDT model. Example 2 . Consider a 3 × 3 depropanizer column system14
D3det G = −2842360.4309.
Then, the model for adj G(s) is identified as ⎡ det G11 −det G21 det G31 ⎤ ⎢ ⎥ adjG(s) = ⎢−det G12 det G22 −det G32 ⎥ ⎢ ⎥ ⎣ det G13 −det G23 det G33 ⎦
0.07724 −56s ⎤ e ⎥ 96s + 1 ⎥ 0.19996 −35s ⎥ e ⎥ 51s + 1 ⎥ −0.5 −17s ⎥ e ⎥ ⎦ 18s + 1 det Gij
Following the procedure given in section 4, Dk G 4, 5 and Ddet , k = 1, 2, 3 are calculated as k
⎡− 14268267200 − 213147344400 − 15264766700 ⎤ ⎢ ⎥ = ⎢− 479771085900 − 207495192500 140876923022500 ⎥ ⎢⎣12672953361740400 − 211369513100 − 483503497900 ⎥⎦
D1det G = −294.7098, D2det G = 17673.5390,
13
G (s ) = ⎡ −0.26978 −27.5s 1.978 e e−53.5s ⎢ 97.5 s + 1 118.5 s + 1 ⎢ ⎢ 0.4881 −117s 5.26 e e−26.5s ⎢ 56 s + 1 58.5 s + 1 ⎢ ⎢ −0.6 5.5 e−16.5s e−15.5s ⎢ ⎣ 40.5s + 1 19.5s + 1
⎡ 62474178.47 ⎤ − 455529171.33 61912536.50 ⎢ ⎥ 1020971022.58 541669958.40 − 97252116464.36⎥ =⎢ ⎢ ⎥ ⎢− 35982075546 559494100.34 − 48350349.79 ⎥ ⎣ 28.51 ⎦
where
, k = 1, 2, 3, 16423
det G11 =
−3.73 e−43.5s 1188s + 76.79s + 1
det G12 =
149.1s 2 − 13.49s + 0.1241 e−51.5s 3 2 59638.54s + 4070.24s + 110.602s + 1
det G13 =
2348s 2 − 50.83s + 5.841 e−43s 76397.05s 3 + 4923.40s 2 + 120.63s + 1
2
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research det G21 =
1.414 e−70.5s 2005s + 130.5s + 1
det G22 =
0.1812 e−44.5s 2921s 2 + 128s + 1
Article
det G(s) =
106.4s + 0.297 det G = e−43s 87770s 3 + 9884s 2 + 199.6s + 1 −9.466s − 0.01076 −82.5s e 20790s 2 + 262.3s + 1
det G32 =
0.09165 e−62.5s 11410s + 195.4s + 1
2
and the determinant transfer function is determined as
23
det G31 =
−2.385 e−54s 12670s + 210.7s + 1
det G33 =
2
1.703 e−71s 16190s + 223.7s + 1 2
According to eq 32, the decoupler is selected as GI (s) = adj G(s) ·D
where D = diag{e43s, e43s, e54s}. The resulted matrix after decoupled is obtained as
2
⎡ 1.703e−28s ⎢ 2 ⎢ 16190s + 223.7s + 1 ⎢ 1.703e−28s GR (s) = det G(s) ·D = ⎢ ⎢ 16190s 2 + 223.7s + 1 ⎢ 1.703e−17s ⎢ ⎢⎣ 16190s 2 + 223.7s +
Therefore, the PID controllers are designed by GMP method as
method7 and the proposed method, while the proposed method did a better job when handling the cross-loop interactions. In comparison with analytical decoupling method,18 the proposed approach is also favorable due to much simpler structure of decoupler. The response speed of the proposed method is quicker than analytical decoupling method but not as fast as normalized decoupling method. However, the decoupling results of these two methods are improved by the proposed method.
⎡ gc ,11 ⎤ ⎢ ⎥ gc ,22 ⎥ GC (s)1 = ⎢ ⎢ ⎥ gc ,33 ⎦ ⎣
where gc ,11 = gc ,22 = 1.8423 +
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1 ⎥⎦
0.008235 + 133.3s s
and gc,33 = 3.0343 +
0.0136 + 219.6055s s
To test the robustness of the proposed method,15−17 we mismatch the process model by increasing all nine steady-state gains, nine time constants, and nine time delays by a factor of 1.4, simultaneously, with all decouplers and controller parameters kept the same as before. The comparison results are shown in Figure 5. It shows that under such model mismatches, the deterioration is reasonable compared with the size of the perturbation. Example 3 . Consider a 4 × 4 HAVC process7 G(s) = ⎡ − 0.098 −17s e ⎢ ⎢ 122s + 1 ⎢ − 0.043 −25s e ⎢ ⎢ 147s + 1 ⎢ − 0.012 −31s e ⎢ ⎢ 153s + 1 ⎢ − 0.013 −32s e ⎢ ⎣ 156s + 1
− 0.036 −27s e 149s + 1 − 0.092 −16s e 130s + 1 − 0.016 −34s e 151s + 1 − 0.015 −31s e 159s + 1
− 0.014 −32s e 158s + 1 − 0.011 −33s e 156s + 1 − 0.102 −16s e 118s + 1 − 0.029 −25s e 144s + 1
− 0.017 −30s ⎤ e ⎥ 155s + 1 ⎥ − 0.012 −34s ⎥ e ⎥ 157s + 1 ⎥ − 0.033 −26s ⎥ e ⎥ 146s + 1 ⎥ − 0.108 −18s ⎥ e ⎥ ⎦ 128s + 1
The system output responses are shown in Figure 6. The unit set-points change inr1 at t = 0, r2 at t = 500, r3 at t = 1000, and r4 at t = 1500. As can be seen, the computation efforts and the control structure are almost the same for normalized decoupling
Figure 5. Nominal system responses for a depropanizer (3 × 3) process. 16424
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research
Article
Figure 6. Nominal system responses for HVAC (4 × 4) process.
6. CONCLUSION In this work, an adjoint matrix based decoupling control method for multivariable processes is formulated in detail. The characteristic matrix sequence is introduced to describe the dynamic feature of the process, which is easily calculated according to Table 1 for the process transfer matrix. Then, the characteristic matrix sequence of adjoint matrix can be directly derived from that of the original process transfer function matrix. By solving the linear equations, the model in the specified form is uniquely determined for the adjoint transfer matrix and the determinant transfer function. Finally, the decoupler is designed for the system based on the adjoint transfer matrix and PI/PID controllers are designed individually to meet the control system performance. This method is very easy to understand and implemented on practice. Three multivariable processes of different size are simulated to demonstrate its effectiveness. The results of system output responses are compared with other methods. It shows that the proposed method can provide satisfactory overall performance.
Let us take the derivatives of s on both sides of eq A1b, it is obtained that Ḣk(s) =
= Hk + 1(s) − H1(s) ·Hk(s)
⎧ F1(s) = H1(s) ⎪ ⎪ F2(s) = H1̇ (s) ⎪ ⎨ F3(s) = H1̈ (s) ⎪ ⎪⋮ ⎪ (k − 1) ⎩ Fk(s) = [H1(s)]
APPENDIX For the convenience of analysis, define functions of Fk(s) and Hk(s) as
[g (s)](k) g (s )
obtained that ⎧ F1(s) = H1(s) ⎪ ⎪ F2(s) = H2(s) − 2H12(s) ⎪ ⎪ F (s) = H (s) − 3H (s)H (s) + 2H 3(s) 3 1 2 1 ⎨ 3 ⎪ 2 ⎪ F4(s) = H4(s) − 4H1(s)H3(s) + 12H1 (s)H2(s) 2 4 ⎪ − 3H2 (s) − 6H1 (s) ⎪ ⎩⋮
(A1a)
(A1b)
for k = 1, 2, ..., l. When s = 0, it is easy to verify that ⎧ ⎪ Fk(0) = Dk ⎨ ⎪ ⎩ Hk(0) = Ek
(A4)
By substituting eq A3 into eq A4, the following set of relations are
and
Hk(s) =
(A3)
when k = 1, the following set of relations hold:
■
Fk(s) = [log(g (s))](k)
[g (s)](k + 1) [g (s)](k) [g (s)]′ − g (s ) g 2 (s )
(A2)
(A5)
when s = 0, eq A5 turns into 16425
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426
Industrial & Engineering Chemistry Research ⎧ D1 = f (E1) = E1 1 ⎪ ⎪ D = f (E , E ) = E − 2E 2 2 2 1 2 1 ⎪ 2 ⎪ ⎪ D = f (E , E , E ) = E − 3E H + 2E 3 2 3 3 1 2 1 3 1 ⎨ 3 ⎪ ⎪ D4 = f4 (E1 , E2 , E3 , E4) ⎪ = E − 4E E + 12E 2E − 3E 2 − 6E 4 4 1 3 1 2 2 1 ⎪ ⎪ ⎩⋮
Article
International Conference on Emerging Technologies and Factory Automation, Hamburg, Germany, Sep 15−18, 2008; pp 1318−1225. (14) Wang, Q.-G.; Zhang, Y.; Chiu, M.-S. Non-interacting control design for multivariable industrial Processes. J. Process Control 2003, 253−265. (15) Wang, Q.-G. Decoupling with internal stability for unity output feedback systems. Automatica 1992, 28, 411−415. (16) Linnemann, A.; Wang, Q.-G. Block decoupling with stability by unity output feedback−Solution and performance limitations. Automatica (UK) 1993, 29, 735−744. (17) Wang, Q.-G. Decoupling control; Springer-Verlag: BerlinHeidelberg, Germany, 2003. (18) Liu, T.; Zhang, W. D.; Gao, F. Analytical decoupling strategy using a unity feedback control structure for MIMO processes with time delays. J. Process Control 2007, 17, 173−186.
(A6)
Therefore, eq A6 is exactly the explicit expression of eq 8.
■
AUTHOR INFORMATION
Corresponding Author
*Tel./Fax: +86-21-2602 7776. E-mail: shenyuling2011@gmail. com. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We are very grateful to the editors and anonymous reviewers for their valuable comments and suggestions to help improve our paper. This work is supported by National Key Basic Research Program of China (No. 2012CB720500), the Key Program of National Natural Science Foundation of China (No. 60736021), National Natural Science Foundation of China-French National Research Agency (No. 61061130563), and the Specialized Research Fund for Doctoral Program of Higher Education of China (No. 20100101110066).
■
REFERENCES
(1) Chen, D.; Seborg, D. E. Design of decentralized PI control systems based on Nyquist stability analysis. J. Process Control 2003, 13, 27−39. (2) Lee, J.; Kim, D. H.; Edgar, T. F. Static decouplers for control of multivariable processes. AIChE J. 2005, 51, 2712−2720. (3) Pomerleau, D.; Pomerleau, A. Guide lines for the tuning and the evaluation of decentralized and decoupling controllers for processes with recirculation. ISA Trans. 2001, 40, 341−351. (4) Wallter, M.; Wallter, J. B.; Wallter, K. V. Decoupling revisited. Ind. Eng. Chem. Res. 2003, 42, 4575−4577. (5) Wang, Q.-G.; Zou, B.; Lee, T. H.; Bi, Q. Auto-tuning of multivariable PID controllers from decentralized relay feedback. Automatica 1997, 33, 319−330. (6) Cai, W.-J.; Ni, W.; He, M.-J.; Ni, C.-Y. Normalized decoupling−A new approach for MIMO process control system design. Ind. Eng. Chem. Res. 2008, 47, 7347−7356. (7) Shen, Y. L.; Cai, W. -J.; Li, S. Y. Normalized decoupling control for high dimensional MIMO processes with application to room temperature control of HVAC systems. Control Eng. Pract. 2010, 18, 652−664. (8) 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−124. (9) Skogestad, S. Simple analytic rules for model reduction and PID controller tuning. J. Process Control 2003, 13, 291−309. (10) Lee, Y.; Park, S.; Lee, M.; Brosilow, C. PID controller tuning for desired closed-loop responses for SI/SO systems. AIChE J. 1998, 44, 106−115. (11) Carvajal, J.; Chen, G.; Ogmen, H. Fuzzy PID controller: design, performance evaluation and stability analysis. Inf. Sci. 2000, 123, 249− 270. (12) Wood, R. K.; Berry, M. W. Terminal composition control of a binary distillation column. Chem. Eng. Sci. 1973, 28, 1707−1717. (13) Morrilla, F.; Vázquez, F.; Garrido, J. Centralized PID control by decoupling for TITO processes. Proceedings of the 13th IEEE 16426
dx.doi.org/10.1021/ie302048g | Ind. Eng. Chem. Res. 2012, 51, 16419−16426