Model Reduction Using Proper Orthogonal Decomposition and Predictive Control of Distributed Reactor System

Abstract

This paper studies the application of proper orthogonal decomposition (POD) to reduce the order of distributed reactor models with axial and radial diffusion and the implementation of model predictive control (MPC) based on discrete-time linear time invariant (LTI) reduced-order models. In this paper, the control objective is to keep the operation of the reactor at a desired operating condition in spite of the disturbances in the feed flow. This operating condition is determined by means of an optimization algorithm that provides the optimal temperature and concentration profiles for the system. Around these optimal profiles, the nonlinear partial differential equations (PDEs), that model the reactor are linearized, and afterwards the linear PDEs are discretized in space giving as a result a high-order linear model. POD and Galerkin projection are used to derive the low-order linear model that captures the dominant dynamics of the PDEs, which are subsequently used for controller design. An MPC formulation is constructed on the basis of the low-order linear model. The proposed approach is tested through simulation, and it is shown that the results are good with regard to keep the operation of the reactor.

1. Introduction

Advances in power computing have propelled the development of process models increasingly detailed and precise, which are then used in the design, optimization, monitoring, and diagnosis of faults, among other tasks. Usually distributed chemical reactors are described by partial differential equations (PDEs) that show the space-time evolution of some variables of interest. In order to simulate PDEs, the spatial domain is discretized, obtaining a large number of ordinary differential equations (ODEs). However, a fine discretization leads to an increase in model complexity. To reduce the complexity of the models, a technique based on orthogonal decomposition of a set of measurements of physical quantities (such as temperature and concentration) is used to represent these quantities along the space and time. This technique, proper orthogonal decomposition (POD), has been used to reduce the order of a large number of systems. This method is based on orthonormal basis functions generated from process data (snapshot matrix) which are obtained by simulation or experimentation on the process; These data are taken by excitation of the process through manipulated variables, external inputs and disturbances of the process. The main idea is to consider POD basis functions which capture the spatial dynamics of our system. The advantage of working with these basis functions is that it is possible to reduce the model order from hundreds or thousands to a few tens. This reduction, resulting in ease of simulation, assimilation and optimization, enables that such models can be applied in real time applications.

There are several methods available for reducing the dimension of a system. The most immediate is a heuristic approach that consists of proposing a priori solution to the equations of motion on the grounds of symmetry and boundary conditions. These solutions usually take the form of a truncated series in terms of general sets of orthogonal functions, such as Fourier modes or spherical harmonics. Antoulas [1] proposes a classification into two main groups, namely, Krylov and singular value decomposition (SVD) methods. Krylov methods make use of iterations for finding approximations to large-scale dynamical systems. The SVD methods are based on the decomposition of the state vector into a set of vectors that can be ordered in a certain sense. These methods include balanced truncation and Hankel approximations for linear systems and the proper orthogonal decomposition (POD) methods and empirical grammians for nonlinear systems. In this work, we are concerned with the POD and its combination with Galerkin projection to produce reduced linear dynamical versions of the original large-scale system (the methodology used in this work is shown in Figure 1). The interested reader is referred to [1] for a more complete account of the available reduction methods. The POD [2, 3] is a statistical technique to extract features from a given dataset by searching for patterns that optimally represent the dataset with respect to quantities such as variance or energy. The output of the POD is a set of time-independent orthogonal functions called empirical orthogonal functions (EOFs). Each EOF is associated with a certain amount of variance or energy. The first suggestion to use the POD in the analysis of dynamical systems was due to Lumley [4]. One of the first phenomena to be studied by means of reduced models arising from the use of the POD was turbulence [5]. Perhaps for nonlinear systems the most studied method is the POD in conjunction with Galerkin projections [1]. However, this is not the only available method, and some other ideas have been proposed, such as a decomposition into principal interaction patterns (PIPs), which explicitly incorporates information about the system’s dynamics within the determination of the patterns [6, 7]. There are other techniques that might be suitable for developing low-order models such as the independent component analysis (ICA) [8]. The ICA aims at answering questions such as, what signal comes from what source? This problem is known as the cocktail-party problem. The solution relies on the assumption of statistical independence of the sources and, therefore, of the signals. The outcome is a set of modes and a mixing matrix that sets the relationship between the modes in order to reconstruct the original mixed signal. The ICA modes are statistically independent but not necessarily orthogonal, and there is no specific order. Furthermore, they are not clearly related to any physical quantity that can help to decide what modes to retain and what to rule out. These features may constitute a disadvantage when trying to construct low-order models since a truncation using these modes would lack an appropriate parameter to measure the expected accuracy of the resulting model. Original algorithms assume linearity although some efforts have been made to extend the formalism to nonlinear systems. So far these ideas have not yet been fully investigated in the context of dynamical systems. Beside this introduction section, this paper has four other sections. Section 2 describes the model reduction method that is used here. Section 3 describes the model of the nonisothermal tubular reactor that is used to illustrate the application of the proposed reduced and control methods. Section 4 addresses the control of a nonisothermal tubular reactor using a reduced model and MPC of infinite horizon. Finally, Section 5 concludes the paper.

Figure 1 

Model reduction framework.

2. Proper Orthogonal Decomposition (POD) and Galerkin Projection

Proper orthogonal decomposition (POD) is a powerful method for data analysis aimed at obtaining low-dimensional approximate descriptions of a high-dimensional process. POD provides a basis for the modal decomposition of an ensemble of functions, such as data obtained in the course of experiments or numerical simulations. The basis functions are commonly called empirical eigenfunctions, empirical basis functions, empirical orthogonal functions, proper orthogonal modes, or basis vectors. The most striking feature of the POD is its optimality: it provides the most efficient way of capturing the dominant components of an infinite-dimensional process with only a finite number of modes, and often surprisingly few modes. In general, there are two different interpretations for the POD. The first interpretation regards the POD as the Karhunen-Loeve decomposition (KLD), and the second one considers that the POD consists of three methods: the KLD, the principal component analysis (PCA), and the singular value decomposition (SVD). In recent years, there have been many reported applications of the POD methods in engineering fields such as in studies of turbulence [9–13], vibration analysis [14–18], process identification [19–22], and control in chemical engineering [23–32]. In general, POD is a methodology that first identifies the most energetic modes in a time-dependent system and subsequently provides a means of obtaining a low-dimensional description of the system’s dynamics where the low-dimensional system is obtained directly from the Galerkin projection of the governing equations on the empirical basis set (the POD modes).

2.1. General Procedure

Let 𝑥(𝑡)∈ℜ^𝑁=[𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑁(𝑡)]^𝑇 be the state vector of a given dynamical system, and let 𝑋∈ℜ^{𝑁×𝑁𝑑} with 𝑁_𝑑≥𝑁 be the so-called snapshot matrix that contains a finite number of samples or snapshots of the evolution of 𝑥(𝑡) at 𝑡=𝑡₁,𝑡₂,…,𝑡_𝑁𝑑. In POD, we start assuming that each snapshot can be written as a linear combination of a set of ordered orthonormal basis vectors (POD basis vectors) 𝜑_𝑗∈ℜ^𝑁,forall𝑗=1,2,…,𝑁: 𝑥𝑡_𝑖=^𝑁_𝑗=1𝑎_𝑗𝑡_𝑖𝜑_𝑗,∀𝑖=1,2,…,𝑁_𝑑,𝑎_𝑗𝑡_𝑖=𝑥𝑡_𝑖,𝜑_𝑗,∀𝑗=1,2,…,𝑁, (1) where 𝑎_𝑗(𝑡_𝑖) is the coordinate of 𝑥(𝑡_𝑖) with respect to the basis vector 𝜑_𝑗 (it is also called time-varying coefficient or POD coefficient) and ⟨⋅,⋅⟩ denotes the Euclidean inner product. Since the first 𝑛 most relevant basis vectors capture most of the energy in the collected data set, we can construct an 𝑛th-order approximation of the snapshots by means of the following truncated sequence: 𝑥_𝑛𝑡_𝑖=^𝑛_𝑗=1𝑎_𝑗𝑡_𝑖𝜑_𝑗,∀𝑖=1,2,…,𝑁_𝑑,𝑛≪𝑁. (2)

In POD, the orthonormal basis vectors are calculated in such a way that the reconstruction of the snapshots using the first 𝑛 most relevant basis vectors is optimal in the sense that the sum-squared-error (SSE) between 𝑥(𝑡_𝑖) and 𝑥_𝑁(𝑡_𝑖), for all 𝑖=1,…,𝑁_𝑑, 𝐸_𝑛=^𝑁𝑑₁₌₁𝑥(𝑡_𝑖)−𝑥_𝑛(𝑡_𝑖)²₂ (3) is minimized. Herein ‖⋅‖₂ denotes the 𝐿₂-norm or the euclidean Norm. In other words, the POD basis vectors are the ones that solve the following constrained optimization problem: min_{𝜑1,…,𝜑𝑛}^𝑁𝑑_𝑖=1𝑥(𝑡_𝑖)−^𝑛_𝑗=1⟨𝑥(𝑡_𝑖),𝜑_𝑗⟩𝜑_𝑗²₂ (4) subject to 𝜑^𝑇_𝑖𝜑_𝑗=1if𝑖=𝑗,0otherwise. (5)

The aim of the previous constraints is to ensure the basis vector orthogonality. The orthonormal basis vectors that solve (4) can be found by calculating the singular value decomposition of the snapshot matrix (𝑋_snap). The singular values of 𝑋_snap are positive real numbers that are ordered in a decreasing way, 𝜎₁≥𝜎₂⋯≥𝜎_𝑛≥0. These values quantify the importance of the basis vectors in capturing the information present in the data. Therefore, the first POD basis vector is the most relevant one, and the last POD basis vector is the least important element. For the application of POD to practical problems, the choice of the 𝑛 most relevant basis vectors is certainly of central importance. A criterion commonly used for choosing 𝑛 based on heuristic considerations is the so-called energy criterion [33]. In this criterion, we check the ratio between the modelled energy and the total energy contained in 𝑋_snap𝑃_𝑛=∑^𝑛_𝑗=1𝜎²_𝑗∑^𝑁_𝑗=1𝜎²_𝑗,𝑛=1,…,(𝑁). (6) The ratio 𝑃_𝑛 is used to determine the truncation degree of the selected POD basis vectors. The number of POD basis elements should be chosen such that the fraction of the first singular values in (6) is large enough to capture most of the information in the data [33]. An ad hoc rule frequently applied is that 𝑛 has to be determined for 𝑃_𝑛=0.99. The closer 𝑃_𝑛 to 1, or similarly the closer 1−𝑃_𝑛 to 0, the better the approximation of 𝑋.

2.2. Galerkin Projection

The derivation of the dynamical model for the POD coefficients can be done in two ways, by using the Galerkin projection or by means of other system identification techniques. The Galerkin projection is the most common way of deriving the dynamical model for the POD coefficients, and it will be the method used in this work.

For explaining the ideas, we suppose that the dynamical behaviour of the high-dimensional system from which we want to find a reduced-order model is described by the following nonlinear model in state space form: 𝑥̇(𝑡)(𝑥(𝑡)(𝑡,𝑦(𝑡)(𝑥(𝑡)(𝑡.=𝑓,𝑢))=𝑔,𝑢)) (7)

Let us define a residual function 𝑅(𝑥̇,𝑥) for (7) as follows: 𝑅(𝑥̇,𝑥)=𝑥̇(𝑡)−𝑓(𝑥(𝑡),𝑢(𝑡)), (8) and let 𝑅(𝑥̇_𝑛,𝑥_𝑛) be the residual when the state vector 𝑥(𝑡) is approximated by its 𝑛th-order approximation: 𝑥_𝑛(𝑡)=^𝑛_𝑗=1𝑎_𝑗(𝑡)𝜑_𝑗=Φ_𝑁𝐚(𝑡),𝑛≪𝑁, (9) where Φ_𝑛=[𝜑₁,𝜑₂,…,𝜑_𝑛] and 𝐚(𝑡)=[𝑎₁(𝑡),𝑎₂(𝑡),…,𝑎_𝑛(𝑡)]^𝑇. In the Galerkin projection, the projection of the residual 𝑅(𝑥̇,𝑥) on the space spanned by the basis vectors Φ_𝑛 vanishes, that is, 𝑅(𝑥̇,𝑥),𝜑_𝑗=0,∀𝑗=1,…𝑛, (10) where ⟨⋅,⋅⟩ denotes the Euclidean inner product. Replacing 𝑥(𝑡) by its 𝑛th-order approximation 𝑥_𝑛(𝑡)=Φ_𝑛𝐚(𝑡) in (7): Φ_𝑛𝐚̇(𝑡)=𝑓Φ_𝑛𝐚(𝑡),𝑢(𝑡), (11) and then we apply the inner product criterion (10) as follows: Φ_𝑛𝐚̇(𝑡),𝜑_𝑗=𝑓Φ_𝑛𝐚(𝑡),𝑢(𝑡),𝜑_𝑗,Φ∀𝑗=1,2,…,𝑛_𝑇𝑛Φ_𝑛𝐚̇(𝑡)=Φ^𝑇_𝑛𝑓Φ_𝑛𝐚(𝑡),𝑢(𝑡), (12) and given that Φ^𝑇_𝑛Φ_𝑛=𝐼_𝑛 because of the orthonormality of the basis vectors, we have the model for the POD coefficients reduce to 𝐚̇(𝑡)=Φ^𝑇_𝑛𝑓Φ_𝑛𝐚(𝑡),𝑢(𝑡). (13)

Finally, the reduced order model of (10) with only 𝑛 states has the following form, 𝐚̇(𝑡)=Φ^𝑇_𝑛𝑓Φ_𝑛𝐚(𝑡)(𝑡,𝑦(𝑡)Φ,𝑢)=𝑔_𝑛𝐚(𝑡),𝑢(𝑡). (14)

3. Tubular Chemical Reactor with Axial and Radial Diffusion

In the following subsection we will describe in detail the kind of tubular reactor for which we will design and implement POD-based model predictive control (MPC) control strategies.

3.1. Nonisothermal Tubular Reactor Model with Axial and Radial Diffusion, Reaction, and Convection

The system to be controlled is a nonisothermal tubular reactor with three phenomena: axial and radial diffusion, reaction and convection, where a reversible, second-order, exothermic, catalyzed reaction takes place (𝐴+𝐵⇆2𝐶). The reactor is surrounded by 3 cooling/heating jackets as it is shown in Figure 2. The temperature of the jackets fluids (𝑇_𝐽1, 𝑇_𝐽2 and 𝑇_𝐽3) can be manipulated independently in order to control the concentration and temperature profiles in the reactor.

Figure 2 

Tubular chemical reactor with 3 cooling/heating jackets.

It is assumed that radial and axial variations in concentration and temperature are present; also we consider laminar flow regime. In this study we are neglecting the heat transfer effects between the jackets fluids and the reactor wall. Under the previous assumptions, the mass balance of specie 𝐴 and the energy balance on the differential cylinder shown in Figure 3 produce the following equations: 𝜕𝐶𝜕𝑡=𝐷_𝑒𝜕²𝐶𝜕𝑟²+𝐷_𝑒𝑟𝜕𝐶𝜕𝑟+𝐷_𝑒𝜕²𝐶𝜕𝑧²−𝑈_𝑧𝜕𝐶𝜕𝑧+𝑟_𝐴,𝜕𝑇=𝐾𝜕𝑡𝜌𝐶^𝑒_𝑝𝜕²𝑇𝜕𝑟²+𝐾_𝑒𝜌𝐶_𝑝1𝑟𝜕𝑇+𝐾𝜕𝑟𝜌𝐶^𝑒_𝑝𝜕²𝑇𝜕𝑧²−𝑈_𝑧𝜕𝑇+𝜕𝑧Δ𝐻_𝑟𝑥𝜌𝐶_𝑝𝑟_𝐴, (15) where 𝑟_𝐴=−𝑘₀𝑒^{−(𝐸/𝑅𝑇)}𝜌_cat𝐶𝐶_𝐵−𝐶²_𝐶𝐾_eq,𝐾_eq=𝐾_eq0𝑒^{[−(Δ𝐻𝑟𝑥/𝑅)(1/303−1/𝑇)]},𝑈_𝑧=2𝑈₀𝑟1−𝑅²,𝑋_𝐴=𝐶_𝐴0−𝐶𝐶_𝐴0,𝐶_𝐵=𝐶_𝐵0−𝐶_𝐴0𝑋_𝐴,𝐶_𝐶=2𝐶_𝐴0𝑋_𝐴, (16)

Figure 3 

Cylindrical shell of thickness Δ𝑟, length Δ𝑧, and volume 2𝜋𝑟Δ𝑟Δ𝑧.

with the following boundary conditions:(i)radial: (1)at 𝑟=0, we have symmetry 𝜕𝑇/𝜕𝑟=0 and 𝜕𝐶/𝜕𝑟=0,(2)at 𝑟=𝑅, the temperature flux to the wall on the reaction side equals the convective flux out of the reactor into shell side of the heat exchanger, −𝐾_𝑒𝜕𝑇𝜕𝑟_𝑅=𝐻_𝑤𝑇(𝑅,𝑧)−𝑇_𝑤, (17)(3)at 𝑟=𝑅, there is no mass flow through the tube walls 𝜕𝐶/𝜕𝑟=0;(ii)axial: (1)𝑇=𝑇_in and 𝐶=𝐶_in at 𝑧=0, (2)at the outlet of the reactor 𝑧=𝐿, 𝜕𝑇/𝜕𝑧=0 and 𝜕𝐶/𝜕𝑧=0.

Here, 𝐶(𝑟,𝑧,𝑡) is the reactant concentration in [mol · m⁻³], 𝑇(𝑟,𝑧,𝑡) is the reactant temperature in [K], 𝐷_𝑒 is the diffusivity of all species in [m² · s⁻¹], 𝐾_𝑒 is the thermal conductivity of the reaction mixture in [J · m⁻¹ · s⁻¹ · K⁻¹], 𝐾_eq0 is the equilibrium constant at 300K, 𝑈_𝑧 is the fluid velocity in [m · s⁻¹], (−Δ𝐻_𝑟𝑥) is the heat of the reaction in [J · mol⁻¹], 𝜌 and 𝐶_𝑝 are the density in [kg · m⁻³] and the specific heat in [J · kg⁻¹ · K⁻¹] of the mix, respectively, 𝑘₀ is the rate constant in [m⁶ · mol⁻¹ · s⁻¹ · kg⁻¹], 𝐸 is the activation energy in [J · mol⁻¹], 𝑅 is the ideal gas constant in [J · mol⁻¹ · K⁻¹], 𝐻_𝑤 is the heat transfer coefficient in [J · m⁻² · s⁻¹ · K⁻¹], 𝐿 is the reactor length in [m], 𝐶_in and 𝑇_in are the concentration in [mol · m⁻³] and the temperature in [K] of the feed flow, 𝑧 is the axial coordinate in [m], 𝑟 is the radial coordinate in [m], 𝑡 is the time in [s], and 𝑇_𝑤(𝑧,𝑡) is the reactor wall temperature in [K] which is defined as follows: 𝑇_𝑤=𝑇_𝐽1,0≤𝑧<𝑍_𝑎,𝑇_𝐽2,𝑍_𝑎≤𝑧<𝑍_𝑏,𝑇_𝐽3,𝑍_𝑏≤𝑧<𝐿. (18)

The parameter values of the reactor model are taken from [34]. These values are presented in Table 1.

Table 1 

Values of the reactor parameters with axial and radial diffusion, reaction, and convection.

The temperature of the jacket sections 𝑇_𝐽1, 𝑇_𝐽2, and 𝑇_𝐽3 must be between 280 K and 330 K. In addition, the temperature inside the reactor must be below 400 K in order to avoid the formation of side products. The kinds of disturbances that affect the reactor are principally variations in temperature and concentration of the feed flow. Typically, such variations are in the range of ±10 K for the temperature and ±5% of the nominal value for the concentration. In this system, only the temperature of the feed flow is measured directly.

3.2. Operating Profile

The operating profiles (steady state concentration and temperature profiles) of the reactor are derived by means of an optimization algorithm, which minimizes an objective function subject to the steady state equations of the reactor described by (15) and the input and state constraints defined previously. The steady state model of the reactor is given by the following partial differential equations (PDEs): 𝐷_𝑒𝜕²𝐶𝜕𝑟²+𝐷_𝑒𝑟𝜕𝐶𝜕𝑟+𝐷_𝑒𝜕²𝐶𝜕𝑧²−𝑈_𝑧𝜕𝐶𝜕𝑧+𝑟_𝐴𝐾=0,_𝑒𝜌𝐶_𝑝𝜕²𝑇𝜕𝑟²+𝐾_𝑒𝜌𝐶_𝑝1𝑟𝜕𝑇+𝐾𝜕𝑟𝜌𝐶^𝑒_𝑝𝜕²𝑇𝜕𝑧²−𝑈_𝑧𝜕𝑇+𝜕𝑧Δ𝐻_𝑟𝑥𝜌𝐶_𝑝𝑟_𝐴=0 (19)

with 𝑇(𝑟,0)=𝑇_in, 𝐶(𝑟,0)=𝐶_in, 𝐶_0,𝑗=𝐶_1,𝑗, 𝐶_𝑛+1,𝑗=𝐶_𝑛,𝑗, 𝐶_𝑖,𝑚+1=𝐶_𝑖,𝑚, 𝑇_0,𝑗=𝑇_1,𝑗, 𝑇_𝑛+1,𝑗=𝑇_𝑛,𝑗, and 𝑇_𝑚+1,𝑗=𝑓(𝑇_𝑚,𝑗,𝑇_𝑤), where 𝑓(𝑇_𝑚,𝑗,𝑇_𝑤) is a function that depends on the boundary conditions when 𝑟=𝑅. The discrete version of (19) can be found by replacing the spatial derivatives by first-order upwind and backward differences approximations as follows: 𝜃₁𝐶_𝑖+1,𝑗+𝜃₂𝐶_𝑖,𝑗+1+𝜃₃𝐶_𝑖,𝑗+𝜃₄𝐶_{𝑖−1,𝑗}+𝜃₅𝐶_{𝑖,𝑗−1}−1𝐶_𝑓𝑟_𝐴𝜗=0,₁𝑇_𝑖+1,𝑗+𝜗₂𝑇_𝑖,𝑗+1+𝜗₃𝑇_𝑖,𝑗+𝜗₄𝑇_{𝑖−1,𝑗}+𝜗₅𝑇_{𝑖,𝑗−1}−1𝑇_𝑓𝑧𝜗₆𝑟_𝐴=0,for𝑖=1,2,…,𝑛=30,𝑗=1,2,…,𝑚=30 (20)

with 𝑇_{𝑤𝑚,𝑗}=𝑇_𝐽1=𝑇_𝐽1𝑇_𝑓,∀_𝑗=1,…,𝑧_𝑎,𝑇_𝐽2=𝑇_𝐽2𝑇_𝑓,∀_𝑗=𝑧_𝑎+1,…,𝑧_𝑏,𝑇_𝐽3=𝑇_𝐽3𝑇_𝑓,∀_𝑗=𝑧_𝑏𝜃+1,…,𝑛,₁𝜃₂𝜃₃𝜃₄𝜃₅=𝐷_𝑒Δ𝑟²𝐷_𝑒Δ𝑧²−2𝐷_𝑒Δ𝑟²()−𝑖−22𝐷_𝑒Δ𝑧²−𝑈_𝑧𝐷Δ𝑧_𝑒Δ𝑟²𝑖−1𝑖𝐷_𝑒Δ𝑧²+𝑈_𝑧,𝜗Δ𝑧₁𝜗₂𝜗₃𝜗₄𝜗₅𝜗₆=𝐾_𝑒𝜌𝐶_𝑝Δ𝑟²𝐾_𝑒𝜌𝐶_𝑝Δ𝑧²−𝑈_𝑧−Δ𝑧2𝐾_𝑒𝜌𝐶_𝑝1Δ𝑧²+1Δ𝑟²+𝐾_𝑒𝑖𝜌𝐶_𝑝Δ𝑟²𝐾_𝑒𝜌𝐶_𝑝Δ𝑟²−𝐾_𝑒𝑖𝜌𝐶_𝑝Δ𝑟²𝐾_𝑒𝜌𝐶_𝑝Δ𝑧²+𝑈_𝑧Δ𝑧Δ𝐻_𝑟𝑥𝜌𝐶_𝑝,𝑟_𝐴=−𝑘₀𝑒^{−(𝐸/𝑅𝑇𝑓𝑇𝑖,𝑗)}𝜌_cat𝐶_𝑖,𝑗𝐶_𝐵−𝐶²_𝐶𝐾_eq,𝐾_eq=𝐾_eq0𝑒^{[−(Δ𝐻𝑟𝑥/𝑅)(1/303−1/𝑇𝑓𝑇𝑖,𝑗)]}, (21) where 𝑛 and 𝑚 are the number of sections in axial and radial directions in which the reactor is divided, 𝑧_𝑎 and 𝑧_𝑏 are the reactor sections defining the ending of the first and second jacket, respectively, 𝑇_𝑓 and 𝐶_𝑓 are normalization factors, 𝐶_𝑖,𝑗=𝐶_𝑖,𝑗/𝐶_𝑓 and 𝑇_𝑖,𝑗=𝑇_𝑖,𝑗/𝑇_𝑓 are the normalized concentration and temperature of the 𝑖,𝑗th section of the reactor, 𝑇_{𝑤𝑚,𝑗}=𝑇_{𝑤𝑚,𝑗}/𝑇_𝑓 is the normalized reactor wall temperature of the 𝑖,𝑗th section, and Δ𝑧 and Δ𝑟 are the length and thickness of each section, respectively. The variables are normalized in order to avoid possible numerical problems. The optimization problem that is solved for deriving the operating profiles is defined as follows: min_{𝑇𝐽1,𝑇𝐽2,𝑇𝐽3}𝜔𝑚^𝑚_𝑖=1𝐶_{𝑟𝑖,𝑛}−𝐶_𝑖,𝑛²+(1−𝜔)𝑛⋅𝑚^𝑚𝑛_𝑖=1𝑗=1𝑇_{𝑟𝑖,𝑗}−𝑇_𝑖,𝑗², (22)

subject to ()𝑇steadystatemodel20_𝐽min𝑇_𝑓≤𝑇_𝐽1,𝑇_𝐽2,𝑇_𝐽3≤𝑇_𝐽max𝑇_𝑓,𝑇_𝑖,𝑗≤𝑇_max𝑇_𝑓,for𝑖=1,2,…,𝑚=30,for𝑗=1,2,…,𝑛=30, (23) where 𝐶_{𝑟𝑖,𝑛} is the desired concentration (normalized) at the reactor output, 𝑇_{𝑟𝑖,𝑗} is the desired temperature (normalized) inside the reactor of the 𝑖,𝑗th increment, 𝐶_𝑖,𝑛 is the concentration (normalized) at the reactor output, 𝜔 is a trade-off coefficient, 𝑇_𝐽min and 𝑇_𝐽max are the lower and upper temperature values of the fluids of the jackets, and 𝑇_max is the maximum allowed temperature inside the tubular reactor. The objective function involves an inherent trade-off between minimising conversion and energy costs. The first term of the cost function corresponds to the squared error of the normalized concentration at the reactor output (terminal cost), and the second term is related to the mean squared error of the normalized temperature along the reactor (integral cost). In this problem, and 𝐶_{𝑟𝑖,𝑛} was set to 0, 𝑇_{𝑟𝑖,𝑗} was selected equal to the normalized temperature of the feed flow (𝑇_in=𝑇_in/𝑇_𝑓) for 𝑖=1,2,…,𝑚 and 𝑗=1,2,…,𝑛. The trade-off parameter 𝜔 can take values between 0 and 1. When 𝜔 goes to 1, the reduction of the reactant concentration at the reactor output becomes more important than the temperature deviations. On the other hand when 𝜔 goes to 0, the temperature deviations become more important than the concentration at the reactor output and the risk of the formation of hot spots is reduced. To solve the optimization problem described by (22), an algorithm (a kind of sequential quadratic programming (SQP)) proposed by [35] was used in this work.

The algorithm was executed in MATLAB with the following parameters: 𝑛=30, 𝑚=30, 𝑇_𝑓=320 K, 𝐶_𝑓=500 mol/m³, 𝑇_𝐽min=280 K, 𝑇_𝐽max=340 K, 𝑇_max=335 K, 𝑇𝑜𝑙=10⁻³, and 𝜔=0.7. The maximum allowed temperature (𝑇_max) inside the reactor was chosen 5 degrees below the actual limit (340 K) in order to give the feedback controller enough room of maneuverability. The trade-off coefficient 𝜔 was found by trial and error.

The algorithm was executed using different initial conditions. Along the experiments, one local minima was found. The operating point was given by 𝑇_𝐽1=291.1705K, 𝑇_𝐽2=293.6938K, and 𝑇_𝐽3=294.8280K. The optimal concentration and temperature profiles can be observed in Figures 4 and 5.

(a)

(b)

Figure 4 

Steady state concentration and temperature profiles.

(a)

(b)

(c)

(d)

Figure 5 

Steady state concentration and temperature profiles at 𝑟=0 m and 𝑧=10 m.

3.3. Linear Model

The linear model of the tubular chemical reactor is obtained by linearizing (15) around the jackets temperatures and the operating profiles presented in Figure 4. This linear model is given by 𝜕𝐶^Δ𝜕𝑡=𝐷_𝑒𝜕²𝐶^Δ𝜕𝑟²+𝐷_𝑒𝑟𝜕𝐶^Δ𝜕𝑟+𝐷_𝑒𝜕²𝐶^Δ𝜕𝑧²−𝑈_𝑧𝜕𝐶^Δ𝜕𝑧−𝛼_𝑎𝐶^Δ−𝛼_𝑏𝑇^Δ,𝜕𝑇^Δ=𝐾𝜕𝑡𝜌𝐶^𝑒_𝑝𝜕²𝑇^Δ𝜕𝑟²+𝐾_𝑒𝜌𝐶_𝑝1𝑟𝜕𝑇^Δ+𝐾𝜕𝑟𝜌𝐶^𝑒_𝑝𝜕²𝑇^Δ𝜕𝑧²−𝑈_𝑧𝜕𝑇^Δ−𝜕𝑧Δ𝐻_𝑟𝑥𝜌𝐶_𝑝𝛼_𝑎𝐶^Δ−Δ𝐻_𝑟𝑥𝜌𝐶_𝑝𝛼_𝑏𝑇^Δ, (24) with the following boundary conditions:(i)radial:(1)at 𝑟=0, we have symmetry 𝜕𝑇^Δ/𝜕𝑟=0 and 𝜕𝐶^Δ/𝜕𝑟=0,(2)at 𝑟=𝑅, −𝐾_𝑒𝜕𝑇^Δ𝜕𝑟_𝑅=𝐻_𝑤𝑇^Δ(𝑅,𝑧)−𝑇^Δ_𝑤, (25)(3)at 𝑟=𝑅; 𝜕𝐶^Δ/𝜕𝑟=0,(ii)axial:(1)𝑇^Δ=𝑇^Δ_in and 𝐶^Δ=𝐶^Δ_in at 𝑧=0, (2)at the outlet of the reactor 𝑧=𝐿, 𝜕𝑇^Δ/𝜕𝑧=0 and 𝜕𝐶^Δ/𝜕𝑧=0.

Here, 𝐶^Δ=𝐶−𝐶^∗, 𝑇^Δ=𝑇−𝑇^∗, and 𝑇^Δ_𝑤=𝑇_𝑤−𝑇^∗_𝑤 are the deviations from the steady state of the concentration, temperature, and reactor wall temperature; 𝐶^∗, 𝑇^∗ and 𝑇^∗_𝑤 are the steady state profiles (operating profiles) of the concentration, temperature, and reactor wall temperature, respectively. In order to reduce the infinite dimensionality of (24), the partial derivatives with respect to space are replaced by first-order upwind, and backward differences approximations and if we define the following vectors, 𝑥(𝑡)=𝐶^Δ_1,1,𝐶^Δ_1,2,…,𝐶^Δ_𝑚,𝑛,𝑇^Δ_1,1,…,𝑇^Δ_𝑚,𝑛^𝑇,𝑑(𝑡)=𝐶^Δ_in,𝑇^Δ_in^𝑇,𝑢(𝑡)=𝑇^Δ_𝐽1,𝑇^Δ_𝐽2,𝑇^Δ_𝐽3^𝑇, (26) then (24) can be cast as follows: 𝑥̇(𝑡)=𝐀𝑥(𝑡)+𝐁𝑢(𝑡)+𝐅𝑑(𝑡), (27) where 𝐀∈ℜ^{2𝑚⋅𝑛×2𝑚⋅𝑛}, 𝐁∈ℜ^{2𝑚⋅𝑛×3}, and 𝐅∈ℜ^{2𝑚⋅𝑛×2} are the matrices describing the system, 𝑥(𝑡)∈ℜ^2𝑚⋅𝑛 is the state vector, 𝑢(𝑡)∈ℜ³ is the vector of the inputs, and 𝑑(𝑡)∈ℜ² is the vector of the disturbances. It is important to mention that (27) is a stable linear system; that is, the operating profile around which the reactor is linearized is nominally stable. This property is very important in Section 4.2 in order to design the IHMPC controllers.

Since the spatial domain of the reactor is divided into 𝑛×𝑚=900 sections, the number of states of (27) is equal to 1800. Given that such large number of states makes the design and implementation of feedback controllers for the reactor difficult, in the next section a reduced order model will be derived using POD and Galerkin projection.

3.4. Model Reduction for the Tubular Reactor Using POD

The derivation of a reduced order model of (27) was done in 5 steps. These steps are described in the following subsections.

3.4.1. Generation of the Snapshot Matrix

We have created a snapshot matrix from the system response (𝑋_snap∈ℜ^1800×10000) when independent step changes were made in the input 𝑢(𝑡) and perturbation 𝑑(𝑡) signals on the nonlinear model (15): 𝑋_snap=[𝑥(𝑡=Δ𝑡),𝑥(𝑡=2Δ𝑡),…,𝑥(𝑡=1000Δ𝑡)]. (28) Along the simulations, 10000 samples were collected using a sampling time Δ𝑡 of 1 s. The amplitude of the step changes was chosen in such a way as to produce changes of similar magnitude in the temperature and concentration profiles. This avoids a possible bias in the resulting model.

3.4.2. Derivation of the POD Basis Vectors

The POD basis vectors are obtained by computing the SVD of the snapshot matrix 𝑋_snap: 𝑋_snap=ΦΣΨ^𝑇, (29) where Φ∈ℜ^1800×1800 and Ψ∈ℜ^10000×10000 are unitary matrices and Σ∈ℜ^1800×10000 is a matrix that contains the singular values of 𝑋_snap in a decreasing order on its main diagonal. The left singular vectors, that is, the columns of Φ, Φ=𝜑₁,𝜑₂,…,𝜑₁₈₀₀, (30) are the POD basis vectors.

3.4.3. Selection of the Most Relevant POD Basis Vectors

The 𝑛 most relevant POD basis vectors are chosen using the energy criterion presented in Section 2.1. The plot of 1−𝑃_𝑁 (see (6)) for the first 100 basis vectors is shown in Figure 6. In this problem, we chose the first 𝑛=50 POD basis vectors based on their truncation degree 1−𝑃_𝑛=3.3×10⁻⁴(𝑃_𝑛=0.9996). The 50th order approximation of 𝑥(𝑡) is given by the following truncated sequence: 𝑥_𝑛(𝑡)=⁵⁰_𝑗=1𝑎_𝑗(𝑡)𝜑_𝑗=Φ_𝑁𝐚(𝑡), (31) where Φ_𝑛=[𝜑₁,𝜑₂,…,𝜑₅₀] and 𝐚(𝑡)=[𝑎₁(𝑡),𝑎₂(𝑡),…,𝑎₅₀(𝑡)]^𝑇.

(a)

(b)

Figure 6 

Logarithmic plot of 𝑃𝑛 and 1−𝑃𝑛 for determining the truncation degree of the POD basis vectors in the reactor case.

3.4.4. Construction of the Model for the First 𝑛=50 POD Coefficients

The Galerkin projection is the most common way of deriving the dynamical model for the POD coefficients, and it will be the method used in this work. Let us define a residual function 𝑅(𝑥) for (27) as follows: 𝑅(𝑥̇,𝑥)=𝑥̇(𝑡)−𝐀𝑥(𝑡)−𝐁𝑢(𝑡)−𝐅𝑑(𝑡), (32) and we replace 𝑥(𝑡) by its 𝑛th-order approximation 𝑥_𝑛(𝑡)=Φ_𝑛𝐚(𝑡) in (32); the Galerkin projection states that the projection of 𝑅(𝑥_𝑛) on the space spanned by the basis functions Φ_𝑛 vanishes, that is, 𝑅𝑥_𝑛,𝜑_𝑗=0,𝑗=1,…,𝑛, (33) where ⟨⋅,⋅⟩ denotes the inner product. Replacing 𝑥(𝑡) by its 𝑛th order approximation 𝑥_𝑛(𝑡)=Φ_𝑛𝐚(𝑡) in (27) and applying the inner product criterion (Galerkin projection) to the resulting equation, we have Φ_𝑛𝐚̇(𝑡),𝜑_𝑗=𝐀Φ_𝑛𝐚(𝑡)+𝐁𝑢(𝑡)+𝐅𝑑(𝑡),𝜑_𝑗. (34) By evaluating the inner product in (34), 𝐚̇(𝑡)=Φ^𝑇_𝑛𝐀Φ_𝑛𝐚(𝑡)+Φ^𝑇_𝑛𝐁𝑢(𝑡)+Φ^𝑇_𝑛𝐅𝑑(𝑡), (35) and we obtain the model for the first 𝑛 POD coefficients. The reduced order model of the reactor with only 50 states is then given by 𝐚̇(𝑡)=𝐀_𝐫𝐚(𝑡)+𝐁_𝐫𝑢(𝑡)+𝐅_𝐫𝑑(𝑡),𝑥_𝑛(𝑡)=Φ_𝑛𝐚(𝑡), (36) where 𝐀_𝐫=Φ^𝑇_𝑛𝐀Φ_𝑛, 𝐁_𝐫=Φ^𝑇_𝑛𝐁, and 𝐅_𝐫=Φ^𝑇_𝑛𝐅. The initial condition for 𝑎(𝑡) reads as 𝑎(0)=0.

For validating the reduced order model of the reactor, we applied constant input signals 𝑇^Δ_𝐽1(𝑡)=0.0312 (𝑇^Δ_𝐽1(𝑡)=10K), 𝑇^Δ_𝐽2(𝑡)=−0.0312 (𝑇^Δ_𝐽2(𝑡)=10 K), and 𝑇^Δ_𝐽3(𝑡)=0.0625 (𝑇^Δ_𝐽3(𝑡)=20 K) and constant perturbation signals 𝐶^Δ_in(𝑡)=2×10⁻⁴ (𝐶^Δ_in(𝑡)=0.1 mol/m³) and 𝑇^Δ_in(𝑡)=0.0156 (𝑇^Δ_in(𝑡)=5 K) to both the full-order model (15) and the reduced order model (36), and afterwards we compared their responses. Figures 7, 8, and 9 show the temperature and concentration profiles of the reactor at different time instants and coordinates for each model. In order to measure the quality of the reduced-order model, the averages of the absolute error for the temperature (𝐸_𝑇) and concentration (𝐸_𝐶) were calculated by means of the following formulas: 𝐸_𝑇=1𝑁_𝑠^𝑁𝑠_𝑘=1𝑇(𝑘Δ𝑡)−𝑇_𝑛(,𝐸𝑘Δ𝑡)_𝐶=1𝑁_𝑠^𝑁𝑠_𝑘=1𝐶(𝑘Δ𝑡)−𝐶_𝑛(𝑘Δ𝑡), (37) where 𝑁_𝑠=10000 is the number of time steps and Δ𝑡=1 s. The plots of 𝐸_𝑇 and 𝐸_𝐶 are shown in Figure 10. For the temperature profile, the maximum value for the error 𝐸_𝑇 is 0.7542 K. For the concentration, the maximum peak for the error 𝐸_𝐶 is 3.3×10⁻³ mol/m³. From the previous results, we can conclude that the reduced order model with only 50 states provides an acceptable approximation of the full order model.

(a)

(b)

(c)

(d)

Figure 7 

Temperature and concentration profiles and 𝑡=100 s at 𝑡=1000 s at 𝑟=0 m. Solid line—full-order model. Dashed line—reduced-order model.

(a)

(b)

(c)

(d)

Figure 8 

Temperature and concentration profiles at 𝑡=1500 s and 𝑡=8000 s at 𝑧=10 m. Solid line—full-order model. dashed line—reduced-order model.

(a)

(b)

(c)

(d)

Figure 9 

Temperature and concentration profiles at 𝑡=10000 s: full-order Model and reduced-order model.

(a)

(b)

Figure 10 

Average of the absolute error between the full-order model (15) and the reduced-order model (36).

The discrete-time version of (36) that is used for designing the digital controller was obtained using the discretization method known as zero-order hold (ZOH) with a sampling time of 0.2 s: 𝐚(𝑘+1)=𝐀^∼𝐚(𝑘)+𝐁^∼𝑢(𝑘)+𝐅^∼𝑑(𝑘),𝑥_𝑛(𝑘)=Φ_𝑛𝐚(𝑘), (38) where 𝐀^∼, 𝐁^∼, and 𝐅^∼ are the matrices describing the new system. A modeling approach frequently adopted in model predictive controller (MPC) considers a discrete-time state-space model in the incremental form [36]; hence (38) can be represented in the following form: 𝑥^𝑠()𝑥𝑘+1_𝑑()=𝐼𝑘+1^𝑛𝑦0𝑥0𝑃^𝑠(𝑘)𝑥^𝑑(𝑘)+𝐷^𝑠𝐷^𝑑(𝑘)+𝐹Δ𝑢^𝑠𝐹^𝑑(𝑘),𝑦(𝑘)=𝐶Δ𝑑_𝑠𝐶^𝑑𝑥^𝑠(𝑘)𝑥^𝑑(𝑘), (39)(40) where 𝑥^𝑑(𝑘)=𝑉₁𝐚(𝑘),𝑥^𝑠(𝑘)=𝑉₂𝐚(𝑘−1),Δ𝑢_𝑘=𝑢_𝑘−𝑢_𝑘−1 is the input increment, Δ𝑑_𝑘=𝑑_𝑘−𝑑_𝑘−1 is the disturbance increment, and 𝑉₁, 𝑉₂ are transformation matrices. In the state equation defined in (39), the state component 𝑥^𝑠 corresponds to the integrating poles produced by the incremental form of the model, and 𝑥^𝑑(𝑘)=𝐚(𝑘) corresponds to the system modes. For stable systems, it is easy to show that when the system approaches steady state, component 𝑥^𝑑 tends to zero. 𝑃 is a diagonal matrix with components corresponding to the poles of the system.

4. Model Predictive Control of the Tubular Reactor

Model predictive control (MPC), also referred to as receding horizon control (RHC) or moving horizon control, is a control strategy where a finite or infinite horizon open-loop optimal control problem is solved on-line at each sampling time using the current state of the plant as the initial state, in order to get a sequence of future control actions from which only the first one is applied to the plant. The fact of solving an optimization problem on-line where common plant constraints are included makes MPC different from conventional optimal control which uses a precomputed control law [37]. MPC has been widely adopted by the industrial process control community and implemented successfully in many applications. First of all, the MPC algorithms can handle in a very natural way constraints on both process inputs (manipulated variables or control actions) and process outputs values (controlled variables), which often have a significant impact on the quality, effectiveness, and safety of the production. Additionally, the MPC controllers can take into account the internal interactions within the process, thanks to the multivariable models on which they are typically based. This makes the MPC algorithms a quite suitable option for multivariable process control. Another reason of the success of MPC is the fact that the principle of operation is comprehensible and relatively easy to explain to process operators and engineers. This is an important aspect at the moment of introducing new techniques into industrial practice. the MPC technology can be found in a wide variety of application areas including chemicals, food processing, automotive, and aerospace applications. A recent survey that provides an overview of commercially available model predictive control technology can be found in [38]. Several notable past reviews regarding theoretical and practical aspects of MPC are offered in [37, 39]. Linear MPC refers to a family of MPC schemes in which linear models are used to predict the system’s dynamics, even though the dynamics of the closed-loop system is nonlinear due to the presence of constraints. Along this work, we will deal with MPC controllers based on discrete-time linear time invariant (LTI) models in state space form. model predictive control based on linear models has been successfully implemented in the control of nonlinear distributed parameter systems [40–42]. Nonlinear model predictive control (NMPC) has gained popularity for low-order lumped parameters nonlinear systems, but for large size distributed systems, the computational cost to solve the nonconvex NLP problem is still excessive.

4.1. Infinite Horizon Model Predictive Control [36]

A modelling approach frequently adopted in model predictive controller (MPC) considers a discrete-time state-space model in incremental form [43]: 𝑥()(𝑘)(𝑘)(𝑘),𝑦(𝑘)(𝑘).𝑘+1=𝐴𝑥+𝐵Δ𝑢+𝐹Δ𝑑=𝐶𝑥 (41)

It is clear that the model defined in (39) and (40) is in the form defined in (41) if the state 𝑥 and model matrices 𝐴, 𝐵, 𝐹, and 𝐶 are properly defined.

MPC is usually based on a discrete-time state-space model as shown in (41). The cost of the infinite horizon MPC considered here can be defined as follows: 𝐽_𝑘,∞=^∞_𝑗=1𝑒𝑘+𝑗^𝑇+𝑄𝑒𝑘+𝑗_{𝑝−1𝑗=0}Δ𝑢𝑘+𝑗^𝑇𝑅Δ𝑢𝑘+𝑗, (42) where 𝑒(𝑘+𝑗)=𝑦(𝑘+𝑗)−𝑟(𝑗),𝑦(𝑘+𝑗) is the output prediction at time instant 𝑘+𝑗 made at time 𝑘, 𝑟 is the desired output reference, 𝑝 is the control horizon, and 𝑄∈ℜ^{𝑛𝑦×𝑛𝑦} and 𝑅∈ℜ^{𝑛𝑢×𝑛𝑢} are positive definite weighting matrices. The controller that is based on the minimization of the above cost function corresponds to the MPC [36] for the output-tracking case. Most of the infinite horizon controllers reduce to finite horizon controllers by defining a terminal state penalty 𝑄. For the cost defined in (42), such a terminal penalty is computed by the following Lyapunov equation [44]: 𝑄−𝑃^𝑇𝑄𝑃=𝑃^𝑇𝐶^𝑑𝑇𝑄𝐶^𝑑𝑃. (43) Since an infinite horizon is used and the model defined in (39) has integrating modes, terminal constraints must be added to prevent the cost from becoming unbounded. Hence, constraints can be written as follows: 𝐶^𝑠𝑥^𝑠(𝑘)−𝑦_𝑠𝑝+𝐶𝐷^𝑠∼𝑠Δ𝑢_𝑘=0, (44) where 𝐷^∼𝑠=𝐷^𝑠⋯𝐷^𝑠,𝐶^𝑠=diag𝐶^𝑠⋯𝐶^𝑠,Δ𝑢_𝑘=Δ𝑢(𝑘)^𝑇Δ𝑢(𝑘+1)^𝑇⋯Δ𝑢𝑘+𝑝−1^𝑇. (45) With the terminal penalty, the cost defined in (42) reduces to 𝐽_𝑘,𝑝=^𝑝_𝑗=1𝑒𝑘+𝑗^𝑇+𝑄𝑒𝑘+𝑗_{𝑝−1𝑗=0}Δ𝑢𝑘+𝑗^𝑇𝑅Δ𝑢𝑘+𝑗+𝑥^𝑑(𝑘+𝑚)^𝑇𝑄𝑥^𝑑(𝑘+𝑚). (46)

Finally, the control optimization problem of the infinite horizon MPC can be formulated as min_Δ𝑢𝑘𝐽_𝑘,𝑝=^𝑝_𝑗=1𝑒𝑘+𝑗^𝑇+𝑄𝑒𝑘+𝑗_{𝑝−1𝑗=0}Δ𝑢𝑘+𝑗^𝑇𝑅Δ𝑢𝑘+𝑗+𝑥^𝑑(𝑘+𝑚)^𝑇𝑄𝑥^𝑑(𝑘+𝑚) (47) subject to modelinincrementalform(39)𝐶^𝑠𝑥^𝑠(𝑘)−𝑦_𝑠𝑝+𝐶𝐷^𝑠∼𝑠Δ𝑢_𝑘=0,−Δ𝑢^max≤Δ𝑢𝑘+𝑗≤Δ𝑢^max,𝑢^min≤𝑢(𝑘−1)+^𝑗_𝑖=0Δ𝑢(𝑘+𝑖)≤𝑢^max,𝑗=0,1,…,𝑝−1. (48)

For large changes on 𝑥^𝑠(𝑘) or 𝑦_𝑠𝑝 or if 𝑦_𝑠𝑝 corresponds to an unreachable steady state, then the optimization problem defined through (47)-(48) may become infeasible because of a conflict between constraints. Consequently, the MPC as defined above cannot be implemented in practice.

4.2. Extended Infinite Horizon Model Predictive Control [45]

To produce an infinite horizon MPC, which is implementable in practice, the objective function of infinite horizon MPC is redefined as follows: 𝐽_𝑘,∞=^∞_𝑗=1𝑒(𝑘+𝑗)−𝛿_𝑘^𝑇𝑄𝑒𝑘+𝑗−𝛿_𝑘+^𝑝−1_𝑗=0Δ𝑢𝑘+𝑗^𝑇𝑅Δ𝑢𝑘+𝑗+𝛿^𝑇_𝑘𝑆𝛿_𝑘, (49) where 𝛿_𝑘∈ℜ^𝑛𝑦 is a vector of slack variables and 𝑆∈ℜ^{𝑛𝑦×𝑛𝑦} is assumed positive definite. Observe that each slack variable refers to a given controlled output. Weight matrix 𝑆 should be selected such that the controller tends to zero the slacks or at least minimize them depending on the number of inputs, which are not constrained. Analogously to the MPC, the extended infinite horizon controllers reduce to finite horizon controllers by defining a terminal state penalty 𝑄 that is obtained by solving (43), and terminal constraints must be added to prevent the cost from becoming unbounded; this constraint can be written as follows: 𝐶^𝑠𝑥^𝑠(𝑘)−𝑦_𝑠𝑝+𝐶𝐷^𝑠∼𝑠Δ𝑢_𝑘−𝛿_𝑘=0. (50) Hence, the control objective defined in (49) becomes as follows: 𝐽_𝑘,𝑝=^𝑝_𝑗=1𝑒(𝑘+𝑗)−𝛿_𝑘^𝑇𝑄𝑒𝑘+𝑗−𝛿_𝑘+𝑥^𝑑(𝑘+𝑚)^𝑇𝑄𝑥^𝑑()+𝑘+𝑚_{𝑝−1𝑗=0}Δ𝑢𝑘+𝑗^𝑇𝑅Δ𝑢𝑘+𝑗+𝛿^𝑇_𝑘𝑆𝛿_𝑘. (51)

Finally, the control optimization problem of the extended infinite horizon MPC can be formulated as follows: min_{Δ𝑢𝑘,𝛿𝑘}𝐽_𝑘,𝑝=^𝑝_𝑗=1𝑒(𝑘+𝑗)−𝛿_𝑘^𝑇𝑄𝑒𝑘+𝑗−𝛿_𝑘+𝑥^𝑑𝑘+𝑗^𝑇𝑄𝑥^𝑑+𝑘+𝑗_{𝑝−1𝑗=0}Δ𝑢𝑘+𝑗^𝑇𝑅Δ𝑢𝑘+𝑗+𝛿^𝑇_𝑘𝑆𝛿_𝑘 (52) subject to modelinincrementalform(39)𝐶^𝑠𝑥^𝑠(𝑘)−𝑦_𝑠𝑝+𝐶𝐷^𝑠∼𝑠Δ𝑢_𝑘−𝛿_𝑘=0,−Δ𝑢^max≤Δ𝑢𝑘+𝑗≤Δ𝑢^max,𝑢^min≤𝑢(𝑘−1)+^𝑗_𝑖=0Δ𝑢(𝑘+𝑖)≤𝑢^max,𝑗=0,1,…,𝑝−1. (53)

Theorem 1. For a stable system, if in the control objective defined in (51) the slack weight matrix 𝑆 is positive definite, then the control law produced by the solution to the problem defined in (52) drives the system output to the desired set point if it corresponds to a reachable state. If the desired set point is not reachable, the controller will stabilize the system in a steady state such that the distance between this steady state and the desired steady state is minimized.

Proof. The proof is provided in [45].

4.3. Extended Infinite Horizon MPC and POD Applied to Control of a Tubular Reactor

The control objective is to reject the disturbances that affect the reactor, that is, the changes in the temperature and concentration of the feed flow. In addition, the control actions must satisfy the input constraints of the process (280K≤𝑇_𝐽1(𝑡),𝑇_𝐽2(𝑡),𝑇_𝐽3(𝑡)≤335K), and the control system should keep the temperature inside the reactor below 335 K. In [46], two POD-based IHMPC control schemes for a nonisothermal tubular reactor were presented. In the first scheme the formulation of the IHMPC controller is formulated in terms of the POD coefficients, and in the second scheme the IHMPC is in terms of physical variables. In the first case, the control of the reactor profiles is achieved indirectly by controlling the POD coefficients which have no physical meaning. This makes the tuning of the controller little intuitive and the definition of the control goals little flexible. This is not the case for the second IHMPC controller, where its formulation is in terms of the temperature of some selected points along the reactor and the concentration at the reactor output. In this work we want to explore the scheme based on the POD coefficients where the control of the temperature and concentration profiles is achieved indirectly and where the references (𝐚_ref) of these POD coefficients can be calculated by 𝐚_ref=Φ^𝑇_𝑛𝑥_ref, (54) where 𝑥_ref is the reference of the vector 𝑥(𝑡) and is equal to zero (the model of the MPC is a discrete-time linear model) since the control system has to keep the reactor operating around the profiles shown in Figure 4. The MPC controller, which uses model (39) to predict the future behaviour of the reactor, is formulated as (52) and (53).

In this formulation 𝐶^𝑑=𝐼_50×50𝑉₁,𝐶^𝑠=0_50×50𝑉₂. Since the state vector 𝑎(𝑘) is unknown and the changes in the concentration of the feed flow (𝑑1(𝑘)=𝐶^Δ_in(𝑘)) are not measured directly, they are estimated by means of an observer (in this case a Kalman filter) with the following formulation: 𝐚()𝑑𝑘+1₁()=𝑘+1^∼𝐀𝐅^∼_𝐶𝐚(𝑘)01𝑑₁(𝑘)+^∼𝐁0𝑢(𝑘)+^∼𝐅_𝑇0𝑑₂(𝑘)+𝐋_𝑎𝐋_𝑑𝑦(𝑘)−𝑦(𝑘),𝑦(𝑘)=𝐶_sel𝑥_𝑛(𝑘)=𝐶_selΦ_𝑛𝐚(𝑘), (55) where 𝐚 is the estimated vector of the POD coefficients, 𝑑₁(𝑘) is the estimation 𝐶^Δ_in, 𝑑₂(𝑘) is the normalized temperature deviation of the feed flow 𝑇^Δ_in(𝑘), 𝑦(𝑘)∈ℜ¹⁰ is a vector containing ten temperature measurements (normalized deviations) along the reactor, 𝑦(𝑘) is the estimate of 𝑦(𝑘), 𝐿_𝑎 and 𝐿_𝑑 are the submatrices of the observer gain (Kalman gain), 𝐹_𝐶 and 𝐹_𝑇 are the column vectors of 𝐹^∼=[𝐹^∼_𝐶,𝐹^∼_𝑇], and 𝐶_sel is a selection matrix which selects the measured temperatures from the vector 𝑥_𝑛(𝑘).

The control horizon 𝑚 was set to 10 samples; 𝑢_min and 𝑢_max were selected according to the input constraints of the process and the operating temperatures of the jackets, and the weighting matrices were in this way 𝑄=1⋅𝐼_50×50,𝑅=10⋅𝐼_3×3, and 𝑆=1⋅𝐼_50×50. The Kalman gain matrix was computed from the following covariance matrices: 𝑅_𝑤=1⋅𝐼_51×51, 𝑅_𝑣=10⁻³⋅𝐼_10×10.

4.4. Simulation Results

In order to evaluate the performance of the control system, the following tests were carried out.

Test 1. The temperature of the feed flow increased to 10 K at the 5000 s, and the concentration of the feed flow decreased to 25×10⁻³ mol/m³ at the 5 s.

Test 2. The temperature of the feed flow decreased to 10 K at the 5000 s, and the concentration of the feed flow increased to 25×10⁻³ mol/m³ at the 5 s.

These disturbances have a big impact on the temperature profile of the reactor.

The simulation results for Test 1 are presented in Figures 11, 12, 13, and 14, and the simulation results for Test 2 are presented in Figures 15, 16, 17, and 18.

(a)

(b)

Figure 11 

Steady state temperature and concentration profiles for Test 1 at 𝑧=10 m. Solid line—MPC. Dashed line—nominal profile (reference).

(a)

(b)

Figure 12 

Steady state temperature and concentration profiles for Test 1 at 𝑟=0 m. Solid line—MPC. Dashed line—nominal profile (reference).

(a)

(b)

Figure 13 

Steady state error (temperature and concentration) for Test 1.

(a)

(b)

(c)

Figure 14 

Control actions (jackets temperatures) of the MPC controller for Test 1.

(a)

(b)

Figure 15 

Steady state temperature and concentration profiles for Test 2 at 𝑧=10 m. Solid line—MPC. Dashed line—nominal profile (reference).

(a)

(b)

Figure 16 

Steady state temperature and concentration profiles for Test 2 at 𝑟=0 m. Solid line—MPC. Dashed line—nominal profile (reference).

(a)

(b)

Figure 17 

Steady state error (temperature and concentration) for Test 2.

(a)

(b)

(c)

Figure 18 

Control actions (jackets temperatures) of the MPC controller for Test 2.

Furthermore, some quantities of interest are given in Table 2. In this table, 𝑇_max is the maximum temperature reached inside the reactor during the test. Δ𝐶_out is the percentage of the change of the mean steady state concentration at the reactor outlet with respect to its nominal value, that is, Δ𝐶_out𝐶%=^𝑁−𝐶^∗_𝑁𝐶^∗_𝑁 (56) where 𝐶^∗_𝑁 is the mean nominal value (0.0339 mol/m³) and 𝐶_𝑁 is the mean concentration at the reactor outlet in steady state after the test.

Table 2 

Performance parameters of the control systems.

In general, the control schemes showed a good behavior for rejecting the disturbances (typical magnitudes: 𝐶_in=±25×10⁻³ mol/m³ and 𝑇_in=±10 K) and both presented a similar performance.

5. Conclusions

In this work, it is shown how POD and Galerkin projections can be used for deriving reduced order model of systems with reaction, diffusion, and convection in two dimensions. The method proposed here is illustrated with a nonisothermal reactor and based on the proposed reduced model, a state observer and a predictive controller are designed and tested.

The algorithm proposed in [35] to find the steadystate operating profiles is extended for the reactor with diffusion, reaction, and convection in two dimensions, and it is shown that it works properly complying with the reactor operating constraints.

The POD method is characterized for its capability to describe the spatial distribution of the relevant physical variables in terms of a set of orthonormal basis functions. These basis functions are selected from observed data and are optimal in a well-defined sense. In the nonisothermal tubular reactor model, the spatial domain is discretized into a high number of grid cells, while in POD models, the spatial distributions are described by the first few and most relevant POD basis functions. The time-dependent characteristics of the variables are given by the time varying coefficients of the POD basis functions. The model of the time varying coefficients is denominated by the reduced order model and is obtained by projecting the POD basis functions onto the original governing equations. Throughout the results presented in this work, it is shown that with very few POD basis functions (less than 3% of the number of grid cells), the temporal and spatial dynamics of the nonisothermal tubular reactor with diffusion, reaction and, convection have an acceptable approximation.

In the application of POD technique, the data matrix (𝑋_snap) was taken from the nonlinear system unlike in [23]. The reason why this was done is that the linear model did not capture the reaction kinetics (irreversible reaction) in a desired way. The use of the non-linear model data gives a more realistic sense to the results of this work because usually the data would be taken from the process. However, using nonlinear model data increases the number of basis functions.

In Section 4, an infinite horizon MPC controller has been designed for the nonisothermal tubular reactor model with axial and radial diffusion on the basis of the reduced order model. The true model of the reactor is nonlinear, and the linearized model is of very high order. The control and optimization problem becomes very tractable if the model can be reduced based on a small number of POD basis functions inferred from the open-loop data. It is shown in Section 4 that the desired temperature and concentration distribution can be controlled using the reduced order model as the base model for the controller (Figures 11 and 17); in this case, the control of the reactor profiles is achieved indirectly by controlling the POD coefficients which have no physical meaning. A very important result of this work is the trade-off between complexity and performance; on the one hand it was possible to reduce the complexity of a high-order model to design control systems and estimators. On the other hand in spite of the spatial discretization of the nonlinear PDE’s describing the reactor, the linearization, and the dramatic reduction of the order by means of POD, the controller has a good performance in order to keep the operation of the reactor at a desired operating condition despite the disturbances in the feed flow. However, if larger disturbances are applied to the tubular chemical reactor, the behaviour of the MPC controllers may not be as good as it has been thus far. This is due to the differences between the nonlinear model and linear model and consequently the reduced order model.

Acknowledgment

The authors acknowledge the European 7th framework STREP project Hierarchical and Distributed Model Predictive Control (HD-MPC). Contract no. INFSO-ICT-223854 for funding this work.