# Optimal boundary control of coupled parabolic PDE–ODE systems using infinite-dimensional representation

## Optimal boundary control of coupled parabolic PDE–ODE systems using infinite-dimensional representation

Journal of Process Control 33 (2015) 102–111 Contents lists available at ScienceDirect Journal of Process Control journal homepage: www.elsevier.com...
Journal of Process Control 33 (2015) 102–111

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Optimal boundary control of coupled parabolic PDE–ODE systems using inﬁnite-dimensional representation Leily Mohammadi a , Ilyasse Aksikas b,∗ , Stevan Dubljevic a , J. Fraser Forbes a a b

Department of Chemical and Materials Engineering, University of Alberta, Canada Department of Mathematics, Statistics and Physics, Qatar University, Doha, Qatar

a r t i c l e

i n f o

Article history: Received 1 July 2014 Received in revised form 5 June 2015 Accepted 17 June 2015 Keywords: Parabolic PDE Inﬁnite-dimensional system Linear quadratic optimal control Boundary control Riesz spectral system

a b s t r a c t The optimal boundary control problem is studied for coupled parabolic PDE–ODE systems. The linear quadratic method is used and exploits an inﬁnite-dimensional state-space representation of the coupled PDE–ODE system. Linearization of the nonlinear system is established around a steady-state proﬁle. Using appropriate state transformations, the linearized system has been formulated as a well-posed inﬁnitedimensional system with bounded input and output operators. It has been shown that the resulting system is a Riesz spectral system. The linear quadratic control problem has been solved using the corresponding Riccati equation and the solution of the corresponding eigenvalue problem. The results were applied to the case study of a catalytic cracking reactor with catalyst deactivation. Numerical simulations are performed to illustrate the performance of the proposed controller. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Many chemical and biochemical processes are modelled by coupled PDEs and ODEs. An example of such processes is the hybrid bioreactor that is used for the treatment of volatile organic compounds such as benzene. This process consists of the bubble column bioreactor, which may be best described by a set of ODEs, interconnected with the bioﬁlter, which is described by a set of PDEs. More examples for processes modelled by coupled PDEs–ODEs can be found in [1–4]. Composite PDE–ODE models can also be used to describe the transportation delay accompanied by other chemical processes such as chemical reactions or mixing, where differentialdifference equations fail to model the system . In such systems, the transportation delay can be modelled by a set of PDEs and the other portion of the system can be described by a set of ODEs. Two types of coupling can exist between PDEs and ODEs. The ﬁrst one arises through the boundary conditions of the distributed portion of the process, since the boundary conditions are functions of the state variables of the lumped parameter system. A tubular reactor and a well mixed reactor in series is a simple example of this coupling. These systems are called cascaded PDE–ODE systems and their control has been the subject of a few recent studies (e.g. [6,7]). The second type of coupling takes place in the domain of the PDE, which means the parameters of the distributed dynamics (e.g., the

coefﬁcients) are functions of the states of the lumped parameter system. Examples of this kind of coupling include a catalytic reactor with catalyst deactivation, where the deactivation kinetics are described by a set of ODEs, or a heat exchanger with a time varying heat transfer coefﬁcient, where the fouling dynamics are described by a set of ODEs. Most biochemical processes are also modelled by a set of coupled PDE–ODE with in-domain coupling (e.g., in situ bioremediation). A common approach for controlling distributed parameter systems is to convert the set of PDEs to a set of ODEs using discretization techniques (e.g. ﬁnite difference or ﬁnite element methods). These discretization methods may result in models that do not accurately capture all of the dynamic properties of the original system. In order to ensure accurate ﬁnite-dimensional models, a very ﬁne discretization is needed and can result in high-order state-space systems, which in turn can lead to computationally demanding controllers. In recent years, research on control of DPS has focused on methods that deal with inﬁnite-dimensional nature of these systems [8,9]. In the aforementioned work, distributed parameter systems are formulated in a state-space form, similar to lumped parameter systems, by introducing a suitable inﬁnite-dimensional space setting and associated operators. This approach allows inﬁnite-dimensional controllers to be synthesized directly from the inﬁnite-dimensional realization of the system [8,10,11]. For parabolic PDEs, Christoﬁdes  studied nonlinear order reduction and control of nonlinear systems. For diffusion–convection–reaction systems, which are described by

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

parabolic PDEs, Dubljevic et al.  used modal decomposition to derive ﬁnite-dimensional systems that capture the dominant dynamics of the original PDE which are subsequently used for low dimensional controller design. Also, boundary control of systems described by single invariant parabolic PDEs has been studied by Dubljevic and Christoﬁdes , and in [15,16], the boundary control problem for parabolic PDEs with spatially varying coefﬁcients has been studied and the corresponding Riccati equation has been solved by using the eigenvalues and eigenfunctions of the system generator. Control problems of coupled hyperbolic PDE–ODE have been studied in some recent works [17–19]. In , the problem of optimal regulation of large space structures governed by coupled system of hyperbolic PDEs and ODEs has been studied. The problem is solved by developing some optimality conditions that are based on the Gâteau differential of the performance index with respect to the control. In , a new optimal control approach for a distributed system governed by ﬁrst-order hyperbolic partial differential equations coupled with a lumped parameter system that accurately models a closed-circuit wind tunnel. A general theory has been developed to analyse an optimal control problem of this class of system. The theory is applied to a predictive optimal control design to regulate the test section Mach number during a continuous pitch motion of a test model. The predictive optimal control is found to have a modiﬁed Riccati solution. In , the LQ control problem for a class of composite hyperbolic PDEs and ODEs is formulated and solved by using the state-space approach. The solution of the LQ control problem is achieved by solving the matrix Riccati equations that result from the operator Riccati equation of the inﬁnite-dimensional state-space representation. The designed optimal control policy is implemented on an interacting CSTR–PFR system through a numerical study. In this work, we are interested in the boundary (LQ) control of a system described by a set of nonlinear parabolic PDEs and ODEs, using the inﬁnite-dimensional state-space description. An important advantage of LQ control is that it uses a state feedback law, in which the state feedback gain is calculated off-line by using LTI system’s dynamics and thereby the amount of on-line calculations is reduced, signiﬁcantly. To the best of authors’ knowledge, there is no published work on the inﬁnite-dimensional optimal control of coupled parabolic PDEs–ODEs with in-domain coupling and this work is the ﬁrst step in the study of inﬁnite-dimensional optimal controller for such systems. The approach developed in  cannot be implemented in the case of parabolic PDEs due to the nature of the corresponding operator Riccati equation. The latter cannot be transformed into a matrix Riccati equation as in the hyperbolic PDEs case. Here, the operator Riccati equation will be solved by using the eigenvalues and eigenvectors of the system generator. The paper is organized as follows. Section 2 focuses on the mathematical description of the system of interest. The nonlinear system is to be linearized around a steady state proﬁle. Appropriate state transformations are used to write the linearized system as a well posed inﬁnite-dimensional system. In Section 3, the eigenvalue problem is solved by adopting the method used for the heat equation of composite media . Section 4 deals with the optimal control problem, which is solved using the corresponding Riccati equation. In Section 5, we consider the case study of a tubular reactor wherein the Van de Vusse reaction takes place. This reaction consists of two series parallel reactions. The mass balance for the reactor results in a set of coupled nonlinear parabolic PDEs, in particular a triangular operator. The triangular structure simpliﬁes the computation of the spectrum of the system. It is also assumed that the parameters of the reactive term are modelled by a set of ODEs which represent the deactivation kinetics.

103

2. Mathematical model description Let us consider the following set of quasi-linear parabolic PDEs coupled (in-domain) with a set of nonlinear ODEs:

⎧ 2 ⎪ ⎨ ∂z = D0 ∂ z − V ∂z + F(k, z) 2 ∂t

∂

∂

⎪ ⎩ dk = g(k)

(1)

dt

with the following initial and boundary conditions D0

  ∂z and |=0 = V z|=0 − zin ∂

z(, 0) = z0

and

∂z |=l = 0 ∂

(2)

k(0) = k0 n

where z( · , t) = [z1 ( · , t) · · · zn ( · , t)]T ∈ H := L2 (0, l) denotes the vector of state variables of the distributed parameter subsysT tem, the vector k = [k1 (t) · · · km (t)] ∈ K := Rm is the vector of state variables for the lumped parameter subsystem.  ∈ [0, l] ∈ R and t ∈ [0, ∞) denote position and time, respectively. D0 and V are diagonal matrices of appropriate sizes; F is a Lipschitz continuous nonlinear operator from H ⊕ K into H; and g is a vector of appropriate size whose entries are continuous functions deﬁned in R. Comment 1The approach developed here can be extended to the case where the matrices D0 and V are diagonalizable. Indeed, state transformation can be used to return to the case where the matrices are diagonal. Moreover, in most chemical engineering processes, D0 and V are symmetric and then diagonalizable. The nonlinear system (1)–(2) can be linearized around the steady state proﬁle (zss , kss ) and the resulting linear system is given by: 2

∂z˜ ∂ z˜ ∂z˜ −V = D0 + N1 ()˜z + N2 ()k˜ ∂t ∂ 2 ∂

(3a)

dk˜ ˜ = M0 k. dt

(3b)

Here, the boundary and initial conditions can be written as follows: D0

  ∂z˜ and |=0 = V z˜|=0 − z˜in ∂

z˜ (, 0) = z˜0

and

∂z˜ |=l = 0 ∂

(4)

˜ k(0) = k˜ 0 ,

where z˜ = z − zss and k˜ = k − kss are the state variables in deviation form and N1 , N2 and M0 are the Jacobians of the nonlinear terms evaluated at the steady state. N1 () =

∂F(k, z) |ss , ∂z

N2 () =

∂F(k, z) dg |ss . |ss and M0 = dk ∂k

Comment 2At this stage, one can remark that Eq. (3b) can be solved and injected in Eq. (3a), however, this will convert Eq. (3a) into a state-space equation with an independent term represented by N2 ()k˜ and can be considered as external disturbance. In this case, we will need to develop an optimal controller that takes into consideration disturbance part. To avoid this issue, we prefer to keep the system in its standard state-space equation. Eq. (3a) is of type diffusion–convection–reaction PDE. In view of solving the eigenvalue problem, it is much easier to convert the equation to a diffusion–reaction type. In order to do so, let us consider the following transformation:



 = T˜z = exp

D−1 0 V 2



(5)

By using this transformation, the PDE system (3a)–(3b) can be described in terms of the new state  by the following linear

104

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

diffusion–reaction parabolic PDE coupled with a linear ODE (more details about this conversion can be found in Appendix A): 2

∂ ∂  + M1 () + M2 ()k˜ =D ∂t ∂ 2 dk˜ = M0 k˜ dt

ˆ = D(A). Therefore, if there exists a J that satisﬁes the with D(A) ˆ will be decoupled. following equation, operator A

(6a)

−A11 J + A12 + JA22 = 0

(6b)

Comment 3Eq. (15) is a Sylvester equation and admits a unique solution if and only if (−A11 ) ∩ (A22 ) = ∅ . The solution is given by



where the matrices M1 , M2 and D are given by

J=

1 −1 M1 () = T N1 () − VD−1 0 V T , 4 −1

D = TD0 T

M2 () = TN2 () and

(, 0) = T−1 z˜0 := 0

(7)

˜ k(0) = k˜ 0

x d

:=

xl

 k˜

(8)

The abstract equation is given by the following boundary control system

⎧ ˙ ⎨ x(t) ⎩

= Ax(t),

Bx(t)

= u(t)

y(t)

= Cx(t)

⎧ ˙ ˆ x(t), xˆ (0) = xˆ 0 ⎪ ⎨ xˆ (t) = Aˆ y(t) = Cˆx(t),

ˆ = B−1 and C ˆ = C−1 . where B System (17) is in the form of a standard abstract boundary control problem. Then by following a similar approach to [16,15], it can be converted to a well-posed inﬁnite-dimensional system with bounded input and output operators. Here is the procedure (see [8, p. 122]). Deﬁne a new operator A as ˆx = Aˆ

Aˆx

x(0) = x0 (9)



dxd are absolutely d

x ∈ H : xd and

d2 xd continuous, ∈ H, d 2 and is given by

A=

D

0

0

0

d2 · + d 2



dx V D d = − xd| 2 =l d |=l

M1

M2

0

M0



 · I :=

A12

0

A22

∂· V Bx( · ) = −D + · |0 ∂ 2

xd xl

xˆ ∈ H : xˆ and

dˆx d2 xˆ are a.c., 2 (18) d d



where a.c. means absolutely continuous. If A generates a C0 semigroup on H and there exists a B ∈ L(U, H) such that the following condition holds:



A11



ˆ ∩ ker(B) ˆ = D(A) = D(A)

dˆx V dˆx V ∈ H and D d |=0 = xˆ d |=0 , D d |=l = − xd |=l 2 2 d d

(10)

 (11)

where I is the identity operator. The boundary operator B : H → U := Rn is given by

(17)

Bˆx(t) = u(t) ⎪ ⎩ ˆ

where A is a linear operator deﬁned on the domain

D(A) =

(16)

ˆ

Now, we are in a position to formulate the system (6)–(7) as an abstract boundary control problem on the inﬁnite-dimensional space H = H ⊕ K. First, let us put u = V˜xin and deﬁne a new state vector x=

T11 (t)A12 T22 (t)dt,

where T11 and T22 are C0 -semigroups generated by −A11 and A22 , respectively (see ). Thanks to Sylvester equation (15) and its solution (16), the resulting decoupled system is given by the following inﬁnitedimensional boundary system:

.

V V ∂ ∂ |=0 = |=0 − TV˜zin , D |=l = − |=l 2 2 ∂ ∂

0

The boundary and initial conditions are: D

(15)

(12) =0

Since M2 in Eq. (11) is generally non-zero, xd and xl are coupled. By introducing the following transformation, the system can be transformed into decoupled subsystems:

(13)

ˆ Bu ∈ D(A) and

ˆ ∀ u ∈ U, BBu = u,

(19)

then system (17) can be transformed into the following inﬁnitedimensional system with bounded input operator on the state space H = U ⊕ H



˙ = Ax(t) + Bu(t) x(t)

(20)

y(t) = Cx(t) where the operators A, B and C are given as follows:

A=

0 ˆ AB

0 A

,

B=

I −B

,



C=C B

I



(21) T

˙ and u(t) = u(t) and x(t) = [u(t) xˆ (t) − Bu(t)] are the new input and state variables. Comment 4Consider B˜ = −1 B, then Condition (19) becomes ˜ ∈ D(A) ˜ = u. ˆ Bu and BBu One can assume that B˜ = [Bd | Bl ]T . Then the following conditions are equivalent to (19)

The operator A will be transformed to

(14)

D

dBd (0) V − Bd (0) = I, 2 d

D

dBd (l) V + Bd (l) = 0, 2 d

Bl ∈ K

(22)

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

105

Bd can be any function that satisﬁes the above conditions. For simplicity one can assume that Bd is a matrix of polynomials. Bl is any arbitrary matrix in K. Finally B can be calculated by

 B = B˜ =

Bd + JBl

 .

(23)

Bl

Comment 5The inﬁnite-dimensional system (20) is in the form of a standard inﬁnite-dimensional system and we are in the position to proceed with dynamical properties and optimal control design for this system. It should be mentioned that, since all of the transformations introduced in this section are exact and there was no approximation involved, all of the dynamical properties of the original linearized system are preserved; hence, we can perform analysis and controller formulation on the transformed system (20) and apply the designed controller to the original system. In order to study the dynamical properties and solve the control problem, we need to solve the eigenvalue problem for the system (20).

Fig. 1. Approximation of the reactor as a composite media.

where the operator A is given by (21) and





A11

0

0

A22

A=



0

F21

F22

In the analysis developed in this paper, we will use the concepts of (1) Riesz basis, which is a sequence of vectors in a Hilbert space H such that there exists an equivalent inner product on H with respect to which this sequence is an orthonormal basis of H and (2) Riesz spectral operator, which is a closed linear operator with simple eigenvalues such that the sequence of eigenvectors is a Riesz basis and the closure of the set of its eigenvalues is totally disconnected. The use of these concepts is motivated by the fact that many dynamical PDE systems are generated by non-self-adjoint operators whose eigenvectors may not be orthogonal, but that do form a Riesz basis. Moreover, elements of the state space H can be uniquely represented as a linear combination of the Riesz basis (even if the basis is not orthogonal) by using the corresponding bi-orthogonal sequence (i.e. the eigenvectors of the adjoint operator). In this section, the solution of the eigenvalue problem that was introduced in [15,16] is extended to the case of coupled PDE–ODE systems. There is no general algorithm for analytical solution of an eigenvalue problem for a general form of parabolic operator. Therefore, in this section we will consider the following assumptions:

1 N1 in (3a) is lower triangular, which leads to lower triangular form of A11 . In many chemical engineering processes, one can use a transformation to triangularize the system. 2 The number of state variables in (3a) is two. Extension to more than two variables is straightforward. 3 Assume that M0 is diagonal, which is extensible to diagonalizable case.

In this work, we are interested in the analytical solution of the eigenvalue problem, therefore the simplifying assumptions were necessary. One can use the results developed here with eigenvalues calculated by numerical methods which can be used for more generic form of the operator A11 . Here, the eigenvalue problem of interest is given by the following equation:

A = 

(24)



⎢ ⎣

=⎢

A11 := 3. Eigenvalue problem



F11

0

0

˛22

d2 + h11 d 2

⎤ 0 d2 d2 2 + h22 d

h21



˛11

d1

⎥ ⎥ ⎦

(25)

A22 = Observe that the operator A is a block diagonal operator, therefore the set of eigenvalues consists of the eigenvalues of the two operators on its diagonal, i.e. (A) = (A11 ) ∪ (A22 ). Since A11 is a lower triangular operator according to assumption 2, (A11 ) = (F11 ) ∪ (F22 ), then (A) = (F11 ) ∪ (F22 ) ∪ (A22 ). Note that F11 and F22 have the following form: F =d

d2 + h() · I d 2

where I is the identity operator. It can be shown that the resolvent operator of F is compact, then, the spectrum of the operator F consists only of isolated eigenvalues with ﬁnite multiplicities [8, Lemma A.4.19]. However, the calculation of this spectrum is a challenging issue due to the fact that h depends on the space variable . In , the spectrum is calculated by dividing the space interval into a ﬁnite number N of subintervals [ i−1 ,  i ], in which it is assumed that the values of h are constant and are denoted by hi (see Fig. 1). The procedure of  will be adopted here. The explicit expressions of the eigenvalues and eigenfunctions of the operator F are given in Appendix B. Let n and n be eigenvalues and eigenfunctions of the operator F11 and n and n be eigenvalues and eigenfunctions of the operator F22 (see Appendix B). Then the eigenvalues of the operator A11 consists of eigenvalues of F11 and F22 , i.e. (A11 ) = (F11 ) ∪ (F22 ) = {n , n , n = 1, . . ., ∞}.

(26)

Thanks to the triangular form of the operator A11 , we can easily ﬁnd the associated eigenfunctions, which are given by





n (n I − F22 )

−1

, F21 n

0

, n

n = 1, . . ., ∞

(27)

106

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

Similar procedure can be used to compute the corresponding biorthonormal eigenfunctions. Indeed, we can solve the eigenvalue problem for the adjoint operator A∗11 , which leads to:



n

  ,



( n I − F11 )−1 F21

0

n

n = 1, . . ., ∞

,

(28)

n

A22 is a diagonal matrix and its eigenvalues are {˛11 , ˛22 } with the associated eigenvectors {[1, 0]T , [0, 1]T }. Finally, the eigenvalues of the operator A are given by (A) = {n , n , ˛11 , ˛22 },

n = 1, . . ., ∞

(29)

Proof. Here, we will prove that the operator A11 is a Riesz spectral operator (see ). In Section 3, it has been shown that F11 and F22 are Riesz spectral operators and they both have real countable and distinct eigenvalues. Therefore, they generate C0 -semigroups T11 and T22 , respectively. According to [8, Lemma 3.2.2, p. 114], the operator A11 is the generator of a C0 -semigroup on H denoted T11 (t) and is given by



T11 (t) =

0

0

0

0

0 =

(30)

1 ˆ −A−1 (AB)

= ⎣! ∞

k=0

1 ˜ k

˜ k ⎦, ˆ AB, k

T(t) = (31)

(32)

The right hand side of the previous equation is an immediate consequence of the expression of the resolvent operator of A in terms of its eigenvalues and eigenfunctions (see [8, Eq 2.35, p. 41]). On the other hand, the corresponding eigenfunctions for  n , n ≥ 1 are given by

n =

0 ˜n

(33)

Finally, the corresponding bi-orthonormal eigenfunctions of A are

0 =

1 0

and n = ⎣

∗ 1 ˜n ˆ (AB) n ⎦.

t

T22 (t − s)F21 T11 (s)x1 ds

(36)

On the other hand, the operator A22 is a diagonal ﬁnite-dimensional operator and therefore it is the generator of a C0 -semigroup on K given by T22 = exp(−A22 t). As a consequence, the (diagonal) operator A is the inﬁnitesimal generator of the following C0 -semigroup:



1

(35)

0

Assume that the operator A has eigenvalues { k , k ≥ 1} and ˜ k, ˜ k ), k ≥ 1}. The spectrum of the operabiorthonormal pair {(

tor A is given by (A) = (A) ∪ {0} and the eigenvalue  0 = 0 has a multiplicity m, where m is the number of manipulated variables. The corresponding eigenfunction for  0 can be computed by solving the equation A 0 = 0. Due to the special form of the operator A, the expression of the eigenfunction 0 is given as follows

T22 (t)



1

⎤ ⎡ ⎤ ⎡ ⎤⎫ n 0 0 ⎪ ⎪ ⎥ ⎢ 0 ⎥ ⎢ 0 ⎥⎪ ⎬ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥,⎣ ⎦,⎣ ⎦ 0 ⎪ ⎦ 1 ⎪ ⎪ 0 1 ⎭

T21 (t)

T21 (t)x1 =

The corresponding bi-orthonormal eigenfunctions are

⎧⎡ ⎤ ⎡ ( n I − F11 )−1 F21 n ⎪ ⎪ ⎪ ⎨⎢ ⎥ ⎢ n ˜ n = ⎢0 ⎥,⎢ ⎣0 ⎦ ⎢ ⎪ ⎣0 ⎪ ⎪ ⎩ 0

0

where the operator T21 is given by

and the associated eigenfunctions are

⎧⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎫ n 0 ⎪ 0 0 ⎪ ⎪ ⎪ ⎪ ⎥ ⎢ n ⎥ ⎢ 0 ⎥ ⎢ 0 ⎥⎪ ⎨⎢ ⎬ −1 ( I − F ) F  ⎢ ⎥ 22 21 n ˜n= ⎢ n ⎥,⎢ ⎥,⎢ ⎥

⎥,⎢ ⎪ ⎣0 ⎦ ⎣ 0 ⎦ ⎣ 1 ⎦ ⎣ 0 ⎦⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭



T11 (t)



T11 (t)

0

0

T22 (t)

According to Section 3, it has been shown that the eigenvalues of the operator A consist of eigenvalues of A and 0 with ﬁnite multiplicity m. Moreover, the eigenfunctions of A and its adjoint are bi-orthonormal. Therefore, they form a Riesz basis for H. Furthermore, A is a Riesz spectral operator and then it is the generator of a C0 -semigroup given by

T(t) =

I

0

S(t)

T(t)

where S(t)x =

"t 0

˜n

(38)

˜ T(s)ABxds. 䊐

The next corollary is immediate consequence of [8, Theorem 5.2.10]. It gives # a necessary and sufﬁcient conditions for the extended system (A, B, C) to be ˇ-exponentially stabilizable and ˇ-exponentially detectable..

#

Corollary 1. Consider the linear system (A, B, C) given by Eq. (20). Assume that B and C are ﬁnite rank operators deﬁned by Bu =

m !

bi ui

and

Cx = ( x, c1 · · ·

i=1

x, ck )

#

A necessary and sufﬁcient condition for (A, B, −) to be ˇexponentially stabilizable is that for all n such that n ∈ ˇ+ (A) rank( b1 , n · · ·

(34)

(37)

bm , n ) = 1.

(39)

#

A necessary and sufﬁcient condition for (A, − , C) to be ˇexponentially detectable is that for all n such that n ∈ ˇ+ (A)

4. Inﬁnite-horizon optimal control

rank( c1 , n · · ·

In this section, we are interested in solving the optimal control problem. In order to do so, we need to investigate the dynamical properties of the system. The following theorem states a result about the generation property of the extended system. This is a consequence of the fact that a Riesz spectral operator generates a C0 -semigroup if and only if supRe(n ) < ∞.

Proof. A is a Riesz spectral operator and its eigenvalues consist of eigenvalues of A and 0 with ﬁnite multiplicity m. Diagonal entries of A are Sturm–Liouville operators and the spectrum of A is ﬁnitely bounded (i.e., there exists a ω such that all  ∈ (A) < ω). There+ fore, for any arbitrary ˇ and , ˇ− (Ae ) comprises ﬁnitely many

n≥1

Theorem 1. Consider the operator A given by Eq. (20). Then A is the inﬁnitesimal generator of a C0 -semigroup on H.

ck , n ) = 1.

(40)

eigenvalues and the ﬁrst condition of Theorem 5.2.10 in  holds. Then, the necessary and sufﬁcient condition for ˇ-exponential stabilizability of (A, B, −) reduces to (39). The ˇ-exponential detectability can be proved in a similar way. 䊐

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

Now, we are interested in the design of linear quadratic (LQ) state feedback optimal controller for the inﬁnite-dimensional system (20)–(21). The aim is to minimize the quadratic cost function:

 J(u) =

y(s), y(s) + u(s), u(s) ds

(41)

0

It is known that the solution of this optimal control problem can be obtained by solving the following algebraic Riccati equation (ARE) : Ax1 , x2 + x1 , Ax2 + x1 , C∗ Cx2 − x1 , BB∗ x2 = 0

(42)

When (A, B) is exponentially stabilizable and (C, A) is exponentially detectable, the algebraic Riccati equation (42) has an unique non-negative self-adjoint solution  ∈ L(H) and for any initial state x0 ∈ H the quadratic cost is minimized by the optimal control uo given by uo (s) = −B∗ x(s).

(43)

Let us set x1 = m and x2 = n and assume that nm = n ,  m . Therefore the solution of the optimal control problem for this system can also be found by solving the set of algebraic equations given by: (m + n )nm + Cnm −

∞ !

Bkl nl km = 0

It is assumed that the catalyst deactivation will only affect the preexponential factor of the main reaction, and k1 will be modelled by dk1 = ˛k1 + ˇ, dt

where Cnm = C n , C m , and Bnm = B* n , B* m . For more mathematical details of how to convert the Riccati equation (42) into (44), please see Appendix C.

k1 (0) = k10

(47)

The above equation for the rate of deactivation of the catalyst is equivalent to the common exponential decay assumption that is used for modelling catalyst deactivation. It is in agreement with the observation that the catalyst deactivation consists of three phases: rapid initial deactivation, slow deactivation and stabilization . The model of the reactor will be: 2

∂yA ∂ yA ∂yA = Da −v + rA , 2 ∂t ∂ ∂ 2

∂yB ∂ yB ∂yB −v = Da + rB , ∂t ∂ 2 ∂

(48)

dk1 = ˛k1 + ˇ dt Initial and boundary conditions are: Da

  ∂yA |=0 = v yA |=0 − yAin , ∂

Da

  ∂yB |=0 = v yB |=0 − yBin , ∂

(44)

k,l=0

107

∂yA |=l ∂ ∂yB |=l ∂

= 0,

(49)

= 0,

yA (, 0)

= yA0 (),

5. Case study: cracking reactor

yB (, 0)

= yB0 (),

In this section, the proposed approach is applied to a catalytic cracking reactor with the assumption that the catalyst deactivates with time . The reaction scheme is given by equations

k1 (0)

= k10

By deﬁning the new state and input variables



(50)



u(t) = v yAin − yAin,ss

(51)

the set of equations (48) can be linearized and then converted to a diffusion–reaction system by using the transformation given in Eq. (5). The resulting system is

(52) k1

with the following initial and boundary conditions

k2

A→B→C

(45)

k3

A→C



∂ Da ∂

rA = −(k1 + k3 )yA2 = −k0 yA2 rB = k1 yA2 − k2 yB

Da





=



with the kinetic equations given by

x1

∂ ∂

(46) x(, 0)

x2 x1

v

=l



  u

2

x2

=0

 =−

x2

x1

v 2

= x0

x1 x2



=0

0 (53)

=l

108

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

where



v2

kˆ 1 () = −



− 2 k1ss + k3 yAss

4Da

and kˆ 2 = −

v2 4Da

n , n , n and n can be calculated using the approach discussed in Section 3. T21 (t) can be calculated using T11 (t) and T22 (t) by



− k2

The inﬁnite-dimensional representation of the system (52)–(53) on the Hilbert space H has the form (9), where the operator A is given by:

t

T21 (t)x =

T22 (t − s)FT 11 (s)xds

(63)

0

where F is the off-diagonal element of A11 and is equal to 2k1ss yAss (). By performing simple calculations, T21 (t) becomes ∞ !

T21 (t)x =

n,m=1

e−m t − e− n t x, m Fm , m − m

n

(64)

n

• Computation of T22 : A22 is a scalar and the semigroup generated by it is (54)

 D(A) =

dx1 dx dx are a.c., 2 ∈ H, Da | d d =l d

v dx2 v dx1 v = − x1 |=l , Da | | = − x2 |=l , Da = x1 |=0 2 2 2 d =l d =0

(65)

• Computation of J: Using Eqs. (58) and (60)–(65), and assuming N1 that A12 = , the operator J can be computed as: N2

2

x ∈ H : x and

T22 = e˛t

\$ (55)

  J1

J=

(66)

J2

and the boundary operator B by J1 x =

∞ ! N1 x, n n

(56)

J2 x =

Assuming that, the control variable is x2 , the output operator C is: C = C0 I = [0 1 | 0]

(57)

By performing the transformation (13), the operator A can be converted to a block diagonal form and the decoupled inﬁnitedimensional system (17) will be computed. Using Eq. (16), the operator J in Eq. (13) is



J=

(67)

n + ˛

n=1

∞ ∞ ! ! N1 x, n Fn ,

m

(n + ˛)( m + ˛)

n=1 m=1

m

+

∞ ! N2 x, m=1

T11 (t)A12 T22 (t)dt

(58)

0

T11 and T22 are the C0 -semigroups generated by −A11 and A22 . The operators A11 , A22 and A12 are given by

2

∂ ⎢ Da ∂2 − kˆ 1 () 0 ⎢

A11 = ⎢

⎣ ⎡

A12 = ⎣

2

2k1ss yAss (z) −yA2 ss

Da

∂ − kˆ 2 ∂ 2

⎥ ⎥ ⎥, ⎦ (59)

⎦ , A22 = [˛]

yA2

• Computation of T11 : The operator A11 is a lower triangular operator and as discussed in Section 4 generates the C0 -semigroup T11 given by



T11 (t) =



T11 (t)

0

T21 (t)

T22 (t)

(60)

where T11 and T22 are the semigroups generated by the diagonal elements of −A11 and are given by T11 (t)x =

∞ !

e

−n t

x, n n

(61)

n=1

T22 (t)x =

∞ ! n=1

e− n t x,

n

n

(62)

(68)

⎧ dˆx(t) ˆ x(t) ⎪ = Aˆ ⎪ ⎪ dt ⎪ ⎨

(69)

⎪ ⎪ ⎪ ˆ ⎪ ⎩ Bˆx(t) = u(t) ˆ x(t) y(t) = Cˆ

The abstract boundary control problem (69), can be converted to a well-posed inﬁnite-dimensional system with bounded input and output operators using Eqs. (18)–(21). In Eq. (21), B can be calculated using the discussion in Comment 4. Since B˜ is any arbitrary B1 function that satisﬁes conditions (22), we assume that Bd = B2 and B1 and B2 are both second order polynomials. Using the conditions (22), B1 and B2 are: B1 =

ss

m

ˆ = A−1 , B ˆ = B−1 and C ˆ = C−1 , the Finally, by deﬁning A decoupled abstract boundary control problem becomes:

xˆ (0) = xˆ 0

m m + ˛

−2 4Da l + vl

B2 = −

2

2 +

1 4Da l + vl

2

2

(70)

v

2 +

1 2Da + 4Da + vl 4Da v + v2 l

(71)

Bl is any arbitrary number in R and we assume that Bl = 1. Finally B becomes:

⎡ ⎢

B1 + J1 Bl

⎤ ⎥

B = ⎣ B2 + J2 Bl ⎦

(72)

Bl 6. Numerical simulations In this section the performance of the proposed approach is demonstrated. The LQ controller discussed in the previous section was studied via a simulation that used a nonlinear model of the

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111 Table 1 Model parameters.

1 n=0

Value

Unit

k10 k2 k3 Da

18.1 1.7 4.8 0.5 2 0.7 0 −0.001 9.05 × 10−3

(h weight fraction)−1 h−1 (h weight fraction)−1 m2 h−1 m hr−1 weight fraction weight fraction h−1

yAin yBin ˛ ˇ

n=1,3,5

n=2

n=4

0.8

Parameter

0.6 0.4 0.2

Φ (2) n

v

109

0 −0.2 −0.4 −0.6 −0.8 −1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

z (m) ˆ n. Fig. 2. Second element of 

3 n=0

n=1

n=2

n=3

n=4

n=5

2 1 0

n

Φ (3)

reactor given in Eqs. (46)–(49). Values of the model parameters are given in Table 1 . The control objective is regulation of the trajectory of yB at the desired steady state proﬁle. Deactivation of catalyst has a negative impact on yB and our objective is to calculate the optimal values of inlet yA to keep trajectory of yB at the desired proﬁle and eliminate the effect of deactivation. Using the nominal operating conditions, and the model given in Eqs. (46)–(49), the steady-state proﬁles of yA and yB were computed. Then, the nonlinear model was linearized around the stationary states and transformed to the non-self-adjoint form of Eqs. (52)–(53). As it has been shown in Section 3, the operator A is a Riesz spectral operator and its eigenfunctions constitute a Riesz basis. The advantage here is that the Riccati equation (42) can be solved without having the normal orthogonal basis. Indeed, thanks to the eigenvalues and eigenfunctions of the operator A, the Riccati equation has been converted into the set of algebraic equation (44). Spectra of operators A11 and A22 were calculated using the algorithm discussed in Section 3. In order to compute the spectrum of A11 , it was assumed that the length of the reactor is divided into 50 equally spaced sections and the coefﬁcient of the reaction term is constant in each section. The ﬁrst ﬁve eigenvalues of the operator A11 are:

−1 −2 −3 −4 0

0.1

0.2

0.3

0.4

0.5 z (m)

0.6

0.7

0.8

0.9

1

ˆ n. Fig. 3. Third element of 

 = {−2.39 × 10−5 , −1.34 × 10−4 , −4.46 × 10−4 , −1.12 × 10−3 , − 2.35 × 10−3 } and the ﬁrst ﬁve eigenvalues of A22 are:  = {−2.04 × 10−6 , −1.096 × 10−5 , −5.68 × 10−5 , −2.08 × 10−4 , − 5.78 × 10−4 } Finally, the spectrum of A was computed using Eqs. (32) and (33). The ﬁrst six eigenvalues of A are: (A) = {0, −0.001, −2.39 × 10−5 , −2.04 × 10−6 , −1.34 × 10−4 , − 1.096 × 10−5 }

Fig. 4. Closed loop trajectory of error yB − yBss .

(73)

and the associated eigenfunctions are shown in Figs. 2 and 3. Once the eigenvalues and eigenfunctions of the operator A are calculated, the LQ-feedback controller can be computed using Eq. (44). Note that since  is a self-adjoint operator, n ,  m = m ,  n , therefore nm = mn . As a result, Eq. (44) gives n(n+1) coupled alge2 braic equations that should be solved simultaneously where n is the number of modes that are used to formulate the controller. Since there are two orders of magnitude difference between the ﬁrst and sixth eigenvalues of the operator A, the effect of higher order eigenvalues on the system dynamic is considered to be negligible; therefore, in this work the ﬁrst ﬁve modes were used for numerical simulation. The computed LQ controller was applied to the nonlinear model of the reactor. Simulation of the nonlinear

system was performed using COMSOL® . The closed loop trajectory of error is shown in Fig. 4 and the optimal input trajectory is shown in Fig. 5. Fig. 4 illustrates that, as catalyst deactivates, the controller is able to regulate the trajectory of yB at the desired steady state trajectory. 7. Summary The LQ-controller for boundary control of an inﬁnitedimensional system modelled by coupled parabolic PDE–ODE equations was studied. This work is an important step in formulation of an optimal controller for the most general form of

110

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111 0.705 2

∂ ∂t

0.7

=

TD0

∂ (T−1 ) ∂(T−1 ) − TV + TN1 ()T−1  + TN2 ()k˜ 2 ∂ ∂ 2

= TD0

YA

0.695

d2 T−1 dT−1 ∂␪ ∂  dT−1   + 2TD0 − TV + TD0 T−1 2 d ∂ d d 2 ∂

−TVT−1

0.69

∂ + TN1 ()T−1  + TN2 ()k˜ ∂

It can be easily observed from the expression of T−1 that 0.685

0.68 0

dT−1 d 500

1000

1500

2000

2500 t (s)

3000

3500

4000

4500

5000

D−1 0 V

=

2

 exp

D−1 0 V 2



=

D−1 0 V 2

T−1

−1 −1 D−1 D−1 d2 T−1 0 VD0 V 0 VD0 V −1 = = T 4 4 d 2

Fig. 5. Optimal input trajectory.

Plugging these expressions into the previous equation, one has

distributed parameter systems consisting of coupled parabolic and hyperbolic PDEs, as well as ODEs. The coupled PDE–ODEs were converted to a decoupled form using an exact transformation. Conditions on existence of this transformation were studied. The decoupled system was then formulated as a well-posed inﬁnitedimensional system by extension of the approach introduced in [16,15]. It was shown that the resulting system is a Riesz spectral system and by using this fact the stabilizability of the system was shown under some conditions. The LQ controller was applied to a catalytic ﬁxed bed reactor, where the rate of catalyst deactivation was modelled by an ODE. The closed loop performance of the controller was studied via numerical simulations. It was illustrated that the formulated controller is able to eliminate the effect of the catalyst deactivation.

∂ ∂t

VD−1 0 V

2

∂␪ ∂  + TD0 T−1 4 ∂ ∂ 2 D−1 V ∂  −TV 0 T−1  − TVT−1 + TN1 ()T−1  + TN2 ()k˜ 2 ∂

=T

= TD0 T−1

T−1  + TVT−1

2

∂  + ∂ 2



TN1 ()T−1 − T

VD−1 0 V 4

T−1

 + TN2 ()k˜

which leads to Eq. (6a). Appendix B. Eigenvalues and eigenfunctions of F11 and F22 In this appendix, we are interested in giving the expression of the eigenvalues and eigenfunctions of the operators F11 and F22 . Note that they both have the following form d2 + h() · I d 2

Appendix A. Convection term elimination

F =d

In this appendix, we will give more details about how to convert Eq. (3a) into Eq. (6a) by using the state transformation (5). Using Eq. (3a), the derivative of the function  with respect to time can be written as follows:

The fact that h is a function of space makes the eigenvalue problem challenging. Here, the length of the reactor is divided to a ﬁnite number of segments. It is assumed that at each segment the value of h is constant. Full description of this approach and more details can be found in [16, Section 3]. The eigenvalues i are given by



2

∂ ∂z˜ ∂ z˜ ∂z˜ −V =T = T D0 + N1 ()˜z + N2 ()k˜ ∂t ∂t ∂ 2 ∂

(A.1)

Note that z˜ can be written as a function of  by using the inverse transformation of T, which is given by

 z˜ = T−1  = exp

D−1 0 V 2

where {ωi }i≥1 is the set of solutions of the following resolvent equation: tan(ωi l) =



2i = dωi2 + hi



(A.2)

By injecting Eq. (A.2) in Eq. (A.1), one gets the following

4dωi v 4d2 ωi2 − v2

and the corresponding eigenfunctions are given by i () = a1 i sin(ωi ) + i cos(ωi ) where 1 = 1 ; i = si,i−1 si−1,i−2 . . . s2,1 and

si,i−1

=

sN,N−1 =

sin(ωi−1 i ) + i−1 cos(ωi−1 i ) sin(ωi i ) + i cos(ωi i ) 1 sin(ωN−1 N ) + N−1 cos(ωN−1 N ) hN sin(ωN N ) + N cos(ωN N )

1

=−

i

=

v sin(ω1 1 ) − 2h1 ω1 cos(ω1 1 ) v cos(ω1 1 ) + 2h1 ω1 sin(ω1 1 )

cos(ωi i )(sin(ωi−1 i ) + i−1 cos(ωi−1 i )) − sin(ωi i )(cos(ωi−1 i ) − i sin(ωi−1 i )) . sin(ωi i )(sin(ωi−1 i ) + i−1 cos(ωi−1 i )) + cos(ωi i )(cos(ωi−1 i ) − i sin(ωi−1 i ))

L. Mohammadi et al. / Journal of Process Control 33 (2015) 102–111

Appendix C. Conversion of Riccati equation (42) into (44)

(m + n )nm + Cnm −

In this appendix, we are interested in converting the operator Riccati equation (42) into the set of algebraic equations (44). Let us set x1 = m and x2 = n and assume that nm = n ,  m . The Riccati equation (42) can be written as A m , n + m , A n + m , C∗ C n − m , BB∗ n = 0

(C.1)

Using the fact that  n is an eigenvalue of the operator A and n is the corresponding eigenvector, one has A m , n + m , A n

= m m , n + m , n n =m m , n + n m , n =m nm + n nm = (m + n )nm .

Now let us calculate the last term of the Riccati equation (C.1). In order to do so, we will use the fact that any element in the space can be written uniquely in the Riesz basis, particularly, we will use BB∗ x =

∞ !

BB∗ x, k k

k=0

Therefore, one has ∗

m , BB n

% = =

m , ∞ !

∞ !

& ∗

BB n , k k

k=0

BB∗ n , k m , k

k=0

=

∞ !

n , BB∗ k km

%

k=0

=

∞ !

n ,

k=0 ∞

=

!

∞ !

& ∗

BB k , l l

km

l=0

BB∗ k , l n , l km

k,l=0 ∞

=

!

Bkl nl km

k,l=0

where Bkl = BB* k , l = B* k , B* l . Moreover, if we put Cnm = n , C* C m , then the Riccati equation (C.1) can be written as a set of the following algebraic equations:

∞ !

111

Bkl nl km = 0

(C.2)

k,l=0

References  P.K.C. Wang, On the stability of equilibrium of a mixed distributed and lumped parameter control system, Int. J. Control 3 (2) (1966) 139–147.  S. Tzafestas, Final-value control of nonlinear composite distributed and lumped parameter systems, J. Frankl. Inst. 290 (5) (1970) 439–451.  M. Oh, C.C. Pantelides, A modelling and simulation language for combined lumped and distributed parameter systems, Comput. Chem. Eng. 20 (6–7) (1996) 611–633.  R. Borsche, R.M. Colombo, M. Garavello, On the coupling of systems of hyperbolic conservation laws with ordinary differential equations, Nonlinearity 23 (2010) 2749–2770.  S. Hiratsuka, A. Ishikawa, Optimal control of systems with transportation lags, IEEE Trans. Autom. Control 14 (3) (1969) 237–247.  M. Krstic, Compensating a String PDE in the actuation or sensing path of an unstable ODE, IEEE Trans. Autom. Control 54 (June (6)) (2009) 1362–1368, ISSN0018-9286.  G.A. Susto, M. Krstic, Control of PDE–ODE cascades with Neumann interconnections, J. Frankl. Inst. 347 (February (1)) (2010) 284–314, ISSN0016-0032.  R.F. Curtain, H.J. Zwart, Introduction to Inﬁnite-Dimensional Linear Systems Theory, Springer-Verlag, 1995.  A. Bensoussan, G. Da Prato, M.C. Delfour, S.K. Mitter, Representation and Control of Inﬁnite Dimensional Systems, second ed., Birkhäuser, Boston, 2007.  F.M. Callier, J. Winkin, Spectral factorization and LQ-optimal regulation for multivariable distributed parameter systems, Int. J. Control 52 (1990) 55–75.  F.M. Callier, J. Winkin, LQ-optimal control of inﬁnite-dimensional systems by spectral factorization, Automatica 28 (4) (1992) 757–770.  P.D. Christoﬁdes, Nonlinear and Robust Control of PDE Systems: Methods and Applications to Transport-Reaction Processes, Birkhäuser, 2001.  S. Dubljevic, P. Mhaskar, N.H. El-Farra, P. Christoﬁdes, Predictive control of transport-reaction processes, Comput. Chem. Eng. 29 (2005) 2335–2345.  S. Dubljevic, P.D. Christoﬁdes, Predictive control of parabolic PDEs with boundary control actuation, Chem. Eng. Sci. 61 (18) (2006) 6239–6248, ISSN00092509. http://www.sciencedirect.com/science/article/pii/S0009250906003678  L. Mohammadi, I. Aksikas, S. Dubljevic, J.F. Forbes, Linear quadratic optimal boundary control of a diffusion–convection–reaction system, in: Proceedings of the 18th IFAC World Conference, 2011.  L. Mohammadi, I. Aksikas, S. Dubljevic, J.F. Forbes, LQ-boundary control of a diffusion–convection–reaction system, Int. J. Control 85 (2012) 171–181.  A.A. Moghadam, I. Aksikas, S. Dubljevic, J.F. Forbes, Boundary optimal (LQ) control of coupled PDEs and ODEs, Automatica 49 (2) (2013) 526–533.  N. Nguyen, M. Ardema, Predictive optimal control of a hyperbolic distributed model for a wind tunnel, J. Guid. Control Dyn. 29 (3) (2006) 526–533.  S.K. Bissau and N.U. Ahmed;1; Optimal control of a large space structures governed by a coupled system of ordinary and partial differential equations. Math. Control Signals Syst., 1–18.  F.de. Monte, An analytic approach to the unsteady heat conduction processes in one-dimensional composite media, Int. J. Heat Mass Transf. 45 (2002) 1333–1343.  A.J. Laub, Matrix Analysis for Scientists and Engineers, Siam, 2005.  Z. Emirsajlow, A composite semigroup for the inﬁnite-dimensional differential Sylvester equation, in: Proceedings of the 7th Mediterranean Conference on Control and Automation (MEHD99), 1999.  C. Delattre, D. Dochain, J. Winkin, Sturm–Liouville systems are Riesz-spectral systems, Appl. Math. Comput. Sci. 13 (2003) 481–484.  V.W. Weekman, Kinetics and dynamics of catalytic cracking selectivity in ﬁxedbed reactors, Ind. Eng. Chem. Process Des. Dev. 8 (3) (1969).  L.E. Kallinikos, G.D. Bellos, N.G. Papayannakos, Study of the catalyst deactivation in an industrial gasoil HDS reactor using a mini-scale laboratory reactor, Fuel 87 (12) (2008) 2444–2449.