Linear Model Predictive Control for a Coupled CSTR and Axial Dispersion Tubular Reactor with Recycle

This manuscript addresses a novel output model predictive controller design for a representative model of continuous stirred-tank reactor (CSTR) and axial dispersion reactor with recycle. The underlying model takes the form of ODE-PDE in series and it is operated at an unstable point. The model predictive controller (MPC) design is explored to achieve optimal closed-loop system stabilization and to account for naturally present input and state constraints. The discrete representation of the system is obtained by application of the structure properties (stability, controllability and observability) preserving Cayley-Tustin discretization to the coupled system. The design of a discrete Luenberger observer is also considered to accomplish the output feedback MPC realization. Finally, the simulations demonstrate the performance of the controller, indicating proper stabilization and constraints satisfaction in the closed loop.


Introduction
The modeling of many chemical engineering process plants relies on the description given by either transport-reaction mathematical models, which belong to the class of distributed parameter systems (DPS), or by lumped parameter system models, which represent idealization of the process units where some assumptions of spatial uniformity (mainly due to the mixing) can take place [1]. The transport-reaction processes are modeled as distributed parameter systems and take the form of partial differential equations (PDEs) which are given by parabolic or hyperbolic PDEs.
To apply control methods on PDEs, one approach is the traditional method, which uses lumping techniques to convert the PDEs to a set of ordinary differential equations (ODEs) [2][3][4][5][6]. Due to the high numbers of modes required in this approach, especially when it comes to the parabolic PDE models, this type of simplification leads to the high dimensionality of the ensuing controller. Furthermore, neglecting the nature of the infinite dimensionality in the original setting might result in instability of the closed-loop system. There are several contributions focused on the synthesis of low-order controllers, which address the issue of having the spatially varying nature in transport-reactions systems. These contributions include the analysis of dynamic properties in the frequency domain, nonlinear, and robust controllers for different classes of dissipative PDEs and Lyapunov-based control methodologies (e.g., [1,7,8]).
Along the line of modeling, most complex processes are in the mixed form of distributed and lumped parameter systems (LPS), and the latter are generally modeled by ODEs [9]. The interconnected coupling of DPS and LPS is a challenging task, but in essence is the proper way to address a variety of process units in real world plants. There are two types of possible interactions between PDEs

System Representation
Consider the following coupled CSTR-tubular reactor configuration as the combination of a lumped and parabolic distributed parameter system. This setting is used for some chemical processes (see, e.g., [33]), and can be represented as follows (Figure 1): The mentioned process can be described by the following coupled ODE-PDE system of equations on domain {t ∈ + , ζ ∈ [0, 1]} with algebraic coupled boundaries and initial condition: where the second order linear parabolic PDE corresponds to convection-diffusion reactor. The transport of a property x I (ζ, t) ∈ L 2 (0, 1) -L 2 (0, 1) is a Hilbert space-through the tubular reactor, given with the , the ODE indicates the dynamics of the variable x F ∈ within the CSTR. ζ ∈ [0, 1] is the position and t ≥ 0 is the time variable. ψ ∈ and a 1 ∈ are the constant values responsible for the consumption or generation of x I (ζ, t) and x F (t), respectively. R ∈ + refers to the recycle factor in the system and is considered to be a bounded parameter (0 ≤ R ≤ 1). v ∈ + , D ∈ + are the constant transport velocity and diffusion-constant, respectively. a 2 ∈ is a constant number, and u(t)∈ represents the system input. f (ζ) and d(t) represent a known disturbance which may present in the tubular reactor and can express changes in unit operations, such as temperature. d(t) is considered to be a step function and f (ζ) is given as follows: where y I and k I are constant values. Here, the linear coupled finite-infinite-dimensional system can be rewritten by the following state-space equations: where x is state representing both finite and infinite part of the process x = . For the sake of simplicity in the paper, index F refers to the finite part of the system (ODE), while index I indicates the infinite part (PDE). A is defined as a linear operator L(X) (where X is a real space are absolutely continuous, is the linear disturbance input operator and C= 0 1 0 δ(ζ − 1)(.)dζ is the linear output operator.

Open-Loop Stability
The generality of the ensuing design can be established by performing the stability analysis of the target system given by Equation (1). The eigenvalue problem for the open-loop (u(t) = 0) unstable coupled ODE-PDE system is defined as below: where : λ and Φ are the eigenvalues and eigenfunctions, respectively. The boundary conditions are given by After some simple manipulation and using the boundary conditions defined by Equation (6) in Equation (4), one gets the following: λ and Φ I (ζ) are found numerically from Equation (7a). The solution of Equation (4) with the set of parameters R = 0.55, v = 1.8 and consumption of the desired component in both CSTR and dispersive tubular reactor (ψ = −1 and a 1 = −0.25), shows that most eigenvalues have negative real parts and there is only one unstable eigenvalue in the system (see Figures 2 and 3). To explore the effects of diffusion (D) on eigenvalues placement, several values for diffusion were considered. By analyzing Figure 3, with the same conditions, one can notice that as diffusivity increases, the distribution will shift from complex eigenvalues to the real ones. In this particular chemical engineering system, with the boundary conditions mentioned above, the instability of the coupled ODE-PDE system present with the value assigned for R. As will be discussed in Section 5 by canceling the unstable mode under model predictive control, the stabilization of the system is addressed.

Discrete Operators
In this section, the Cayley-Tustin time discretization is applied, which maps the mentioned coupled ODE-PDE system from continuous-time to a discrete one, preserving all energy properties with the feature of no spatial discretization. The discrete version of Equation (3) with sampling time ∆t can be represented as follows: are the linear discrete operators defined by: where δ = 2/∆t, and [δ − A] −1 = R(δ, A) is given as the resolvent operator of A in Equation (5). It should be mentioned that the discrete operators are found by replacing s with δ in R(s, A) . In order to find the resolvent operator, one may easily apply Laplace transform to the set of Equations (1): Since P is a constant matrix, one can calculate e Pζ with the Laplace inverse transform (L −1 {[sI − P] −1 }). Therefore, the solution of the mentioned system X(ζ, s) = e Pζ X(0, s) + ζ 0 e P(ζ−η) Hdη can be expressed as follows: After applying the boundary conditions to Equation (11), the discrete operators are obtained and have the following form: with the following components: in above equations e 1 (ζ), e 2 (ζ), e 3 (ζ), e 4 (ζ), b 1 (ζ), b 2 (ζ), k 1 (ζ), k 2 (ζ) and f 3 (ζ) are defined by the following expressions: where e (1) 3 are the corresponding terms calculated at ζ = L = 1.

Discrete Adjoint Operators
The adjoint operators are required for developing the model predictive control. The expressions for adjoints (A * d , B * d ) of the discrete operators (A d , B d ) are written in the following form: the components of the adjoint operators are computed based on the definition (< A d Φ, Ψ >=< Φ, A d Ψ >) and results in following operators:

Luenberger Observer Design
In a real process system controller realization, having access to all the state variables cannot be feasible, especially when DPS are considered. In order to address this issue in the design of the model predictive control, the Luenberger observer is introduced for reconstruction of the state variables by taking the output measurement into account. First, the design of observer for the coupled parabolic PDE-ODE system in the continuous setting is considered. Then, the continuous observer gain is transferred into the discrete one using Cayley-Tustin discretization. The Luenberger observer has the form given by: where L c = is the continuous observer gain. By subtracting Equation (21) from its general form (ẋ(ζ, t) = Ax(ζ, t) + Bu(t)), one can get the dynamics of the observer error as follows:ê (ζ, t) = (A − L c C)ê(t),ê(0) = 0 the design of the observer is relied on choosing L c such that the state estimation error dynamics given by Equation (22) is stable. Hence, the stability of the observer can be ensured by analyzing the eigenvalues problem of the observer error: which results in following equation and boundary conditions: . It should be emphasized that in the eigenvalue problem, although a spatially varying L I (ζ) could be chosen, a constant value throughout the whole spatial domain (L I (ζ) = L) is considered and the same value assigned for the finite part (L F = L). It can be seen that as the value of L c increases, the unstable real eigenvalue is shifted to the left side (the stable region). Therefore, as depicted in Figure 4 for L c > 0.1 the stability of the error dynamics will be ensured by having only negative eigenvalues. Here, the resolvent is used to compute the corresponding discrete observer gain from the continuous one [34]: Therefore, the reconstructed state in the discrete setting can be expressed as:

Optimization Problem
The linear discrete-time model dynamics represented in Equations (12)- (20) is used in the formulation of the model predictive control for the coupled CSTR and tubular reactor system. The MPC developed in [31] regarding linear time invariant systems for the finite-dimensional setting is extended to infinite-dimensional one. In order to achieve this purpose, the following objective function should be minimized at each sampling time (k) to design the regulator of the coupled ODE-PDE system: refers to the reconstructed state, F is a positive definite operator, Q= Q F Q I represents positive semidefinite spatial operator associated with the state of coupled ODE-PDE system and the indices, (k + j) and (k + j + 1|k), for both input and state variable, indicate current and future time, respectively. In order to get the finite horizon objective function, one can assume zero input beyond the control horizon (i.e., u(k + N + 1|k) = 0) by taking the terminal penalty term into an account. The result takes the following form: After simple algebraic manipulation, the following finite-dimensional quadratic optimization problem is obtained: whereQ is terminal state penalty operator. The above equation is subjected to the following constraints: where U = u(k + 1|k) u(k + 2|k) u(k + 3|k) . . . u(k + N|k) T and (H, P) are computed as below: The model predictive control scheme used on the coupled ODE-PDE system is illustrated in Figure 5. As can be seen, the full state feedback is needed for the MPC scheme, which is going to be given by the reconstructed states from the observer.

Terminal State Penalty Operator
The terminal state penalty term, the operatorQ= Q F Q I , can be found from the solution of the following discrete Lyapunov equation: The above solution of the discrete Lyapunov equation based on Cayley-Tustin method is the same unique solution ofQ in Equation (31) in continuous setting (A * Q +QA = −Q) [35]. Since the solution ofQ cannot be obtained directly; in other words, because of the necessity of using integral operators coming from discrete operators for calculatingQ, the procedure followed here is to connect the discrete and continuous Lyapunov equations. Here, one can rewrite the continuous Lyapunov equation in the following format [36]: By considering x 1 =Φ m and x 2 =Ψ m and using the fact that λ m andΨ m are the eigenvalue and eigen function of the system (i.e., AΦ m = λ mΦm ), the equation leads to: Q is a bounded symmetric operator (D(A * ) = D(A)) and it is self-adjoint (see [36]) which implies <Φ m ,QΨ m >=<QΦ m ,Ψ m >=Q m , accordingly the following simplified equation is achieved: Finally, the solution of the continuous Lyapunov equation gives the expression for the infinite part (PDE) of the terminal state penalty operator (Q I ) which can be expressed as below: whereΦ m I andΨ m I refer to the normalized eigenfunction and adjoint eigenfunction of infinite part of the system, respectively. In Equation (35), the summation is computed for increasing number of different eigenvalues, until the applied operator converges to a constant value. In this work the first 20 eigen modes are considered in the simulation. For the finite part (ODE), A F = A * F = a 1 . According to Lyapunov equation,Q F is easily obtained and given by the following expression:

Stability Constraint
Based on the definition of a positive definite operator, it is possible to show thatQ is a positive operator if only the stable nodes are taken into account [36]. In order to guarantee stabilization of the system, a stability constraint is applied in the optimization problem and is represented by an equality constraint [35]. It is assumed that the controller will gain stabilization by rejecting the unstable modes. Hence, this condition can be written as below: where Φ u refer to unstable eigenfunctions associated with positive eigenvalues. The corresponding equality constraint, which cancels the unstable modes at end of the horizon, is constructed as follows: If there is a feasible input sequence given by optimization problem, the above equality constraint is satisfied for the constrained convex optimization problem given by Equation (29a). Therefore, stabilization can be obtained, and the unstable modes will be canceled by end of the horizon. Here, due to the feasibility of the optimization represented by constrained quadratic problem in the zero-disturbance case, feasibility implies stability and optimal stabilizability. This extension is based on the well-known results from the finite-dimensional theory [30,31].

Simulation Results
In this section, the simulation study is performed for the proposed controller of the coupled ODE-PDE system. First, the design of the observer is discussed in the discrete setting with reconstruction of states in an open-loop condition by using Cayley-Tustin method, then the performance of the ensuing Model Predictive Control is demonstrated and compared with the open-loop response.

Observer Design and Open-Loop Response
Based on Equations (8)- (13), one can reconstruct the dynamics of the discrete representation for both finite and infinite parts of the system with the set of parameters given in Table 1. As discussed in Section 4, to guarantee the stability of the observer, L c = 5 is chosen as the observer gain. By using Equation (25) the discrete version of the corresponding observer gain is computed. The initial conditions for the mentioned observer and original system are considered to be constants in the entire space,x 0 = 0 and x 0 = 1, respectively. In simulation ∆t = 0.04 is considered which implies δ = 50 and for numerical integration ∆z = 0.005 is chosen. Then, according to Equation (26), the reconstructed state is obtained, and the corresponding error dynamics is evaluated. As shown in Figure 6a, the dynamics of the observer error converges to zero, which means the developed observer has a good performance. Hence, in the case of a realistic system, the Model Predictive Control can be applied using just the output measurement. Table 1. Parameter values used in numerical simulation.

Parameters
Values As discussed in Section 2.2 regarding the instability of the coupled ODE-PDE system, Figure 6b demonstrates the space-time evolution of the tubular reactor (x I (ζ, t)) for 0 ≤ t ≤ 20 which grows unbounded as expected. Following this, the corresponding dynamics of the CSTR (x F (t)) with pertaining initial condition is depicted in Figure 6c.

MPC Implementation
In this part, based on the scheme represented by Figure 5, the successful application of the proposed constrained model predictive controller is demonstrated on the basis of Cayley-Tustin time discretization with optimization problem described by Equations (27)-(38). The implementation of the system under model predictive control with and without constraints (see Equation (29b)) in time domain 0 ≤ t ≤ 20 is shown in Figure 7(b)) and Figure 7(d)). By choosing Q F = 2.5, Q I (ζ) = 1.5 and N = 65 for the MPC control horizon, it is possible to see that the controller is able to comply with the input and states constraints (dash-dotted lines) imposed into the coupled ODE-PDE system. Moreover, regarding the instability of the system, as described in Section 5.3, one can notice that the stability constraint is also fulfilled by canceling the unstable mode at end of the horizon (see Figure 7(c)). On the other hand, in order to provide a comparison of the dynamics of scalar variable x F in CSTR with input and state constraints, the simulation is performed again to justify two scenarios, for the first one only stability constraint is considered in the MPC algorithm while in the second one all constraints (stability, input and state) are present. Figure 7(d)) and Figure 7(e)) show this analogy with corresponding control actions given by Figure 7(a)) and Figure 7(b)). As expected in the first setting, the CSTR dynamics is faster compared to the latter case as no state/input constraints need to be satisfied.
Another simulation has been performed to explore the behavior of the MPC algorithm for handling the mentioned constraints, the step disturbance is applied through the infinite part of the system (tubular reactor) for 20 ≤ t ≤ 25. The idea here is to examine the performance of the MPC algorithm for stabilization and the state/input constrains satisfaction. y I = 0.5 and k I = 0.6 are chosen as the parameters in Equation (2) describing spatial varying function f (ζ). Figure 7(c)) and Figure 7(f)) verify the good performance of MPC for handling the constraints when the disturbance is present. By analyzing Figure 7(c)), one can notice that after the step disturbance is applied, the state variable of the CSTR stays at the upper limit (x max F ) until t = 25, then decreases and once again goes back to zero (at t = 45) based on the control action given by Figure 7(c)).

Conclusions
In this contribution, the design of a model predictive controller and discrete observer for a coupled ODE-PDE system was investigated. In particular, the lumped and distributed system were coupled by the boundaries, with the manipulated variable acting on the ODE. The system stability characteristics were first analyzed by studying the system's eigenvalues. A discrete representation of the system was necessary in the controller design; thus, the Cayley-Tustin time discretization was applied, preserving the original system characteristics. An unstable operation condition was considered, and the MPC and observer design had to take this into account. To develop the discrete observer, the design in the continuous-time setting was first derived, then the discrete observer was obtained based on the continuous gain. The MPC was designed to obtain the optimal control sequence while handling input constraints and stabilizing the system using a terminal constraint. Finally, numerical simulations were shown to present the performance of the controller in the closed loop. As expected, the controller was able to achieve stabilization, while handling the constraints. If a disturbance is made, the controller can deal with the effects made in the system while satisfying the constraints. The constrained optimal control design for the class of parabolic and hyperbolic PDEs, which guarantees the system stability and handling both input and state constraints will be an extension to be addressed in future studies.