Research Article  Open Access
Alejandro Marquez, Jairo José Espinosa Oviedo, Darci Odloak, "Model Reduction Using Proper Orthogonal Decomposition and Predictive Control of Distributed Reactor System", Journal of Control Science and Engineering, vol. 2013, Article ID 763165, 19 pages, 2013. https://doi.org/10.1155/2013/763165
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 discretetime linear time invariant (LTI) reducedorder 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 highorder linear model. POD and Galerkin projection are used to derive the loworder 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 loworder 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 spacetime 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 largescale 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 largescale 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 timeindependent 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 loworder 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 cocktailparty 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 loworder 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.
2. Proper Orthogonal Decomposition (POD) and Galerkin Projection
Proper orthogonal decomposition (POD) is a powerful method for data analysis aimed at obtaining lowdimensional approximate descriptions of a highdimensional 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 infinitedimensional 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 KarhunenLoeve 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 timedependent system and subsequently provides a means of obtaining a lowdimensional description of the system’s dynamics where the lowdimensional 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 socalled 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) : where is the coordinate of with respect to the basis vector (it is also called timevarying 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 thorder approximation of the snapshots by means of the following truncated sequence:
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 sumsquarederror (SSE) between and , for all , 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: subject to
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 (). The singular values of are positive real numbers that are ordered in a decreasing way, . 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 socalled energy criterion [33]. In this criterion, we check the ratio between the modelled energy and the total energy contained in 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 . The closer to , or similarly the closer to , 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 highdimensional system from which we want to find a reducedorder model is described by the following nonlinear model in state space form:
Let us define a residual function for (7) as follows: and let be the residual when the state vector is approximated by its thorder approximation: where and . In the Galerkin projection, the projection of the residual on the space spanned by the basis vectors vanishes, that is, where denotes the Euclidean inner product. Replacing by its thorder approximation in (7): and then we apply the inner product criterion (10) as follows: and given that because of the orthonormality of the basis vectors, we have the model for the POD coefficients reduce to
Finally, the reduced order model of (10) with only states has the following form,
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 PODbased 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, secondorder, exothermic, catalyzed reaction takes place (). The reactor is surrounded by 3 cooling/heating jackets as it is shown in Figure 2. The temperature of the jackets fluids (, and ) can be manipulated independently in order to control the concentration and temperature profiles in the reactor.
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: where
with the following boundary conditions:(i)radial: (1)at , we have symmetry and ,(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, (3)at , there is no mass flow through the tube walls ;(ii)axial: (1) and at , (2)at the outlet of the reactor , and .
Here, is the reactant concentration in [mol · m^{−3}], is the reactant temperature in , is the diffusivity of all species in [m^{2} · s^{−1}], is the thermal conductivity of the reaction mixture in [J · m^{−1} · s^{−1} · K^{−1}], is the equilibrium constant at , is the fluid velocity in [m · s^{−1}], is the heat of the reaction in [J · mol^{−1}], and are the density in [kg · m^{−3}] and the specific heat in [J · kg^{−1} · K^{−1}] of the mix, respectively, is the rate constant in [m^{6} · mol^{−1} · s^{−1} · kg^{−1}], is the activation energy in [J · mol^{−1}], is the ideal gas constant in [J · mol^{−1} · K^{−1}], is the heat transfer coefficient in [J · m^{−2} · s^{−1} · K^{−1}], is the reactor length in [m], and are the concentration in [mol · m^{−3}] 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:
The parameter values of the reactor model are taken from [34]. These values are presented in Table 1.

The temperature of the jacket sections , , and 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):
with , , , , , , , and , 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 firstorder upwind and backward differences approximations as follows:
with 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:
subject to 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 tradeoff coefficient, and are the lower and upper temperature values of the fluids of the jackets, and is the maximum allowed temperature inside the tubular reactor. The objective function involves an inherent tradeoff 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 , was selected equal to the normalized temperature of the feed flow () for and . The tradeoff 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 , 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: , , K, mol/m^{3}, K, K, K, , and . The maximum allowed temperature () inside the reactor was chosen degrees below the actual limit (340 K) in order to give the feedback controller enough room of maneuverability. The tradeoff 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 , , and . The optimal concentration and temperature profiles can be observed in Figures 4 and 5.
(a)
(b)
(a)
(b)
(c)
(d)
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 with the following boundary conditions:(i)radial:(1)at , we have symmetry and ,(2)at , (3)at ; ,(ii)axial:(1) and at , (2)at the outlet of the reactor , and .
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 firstorder upwind, and backward differences approximations and if we define the following vectors, then (24) can be cast as follows: where , , and are the matrices describing the system, 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 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 () when independent step changes were made in the input and perturbation signals on the nonlinear model (15): 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 : where and are unitary matrices and is a matrix that contains the singular values of in a decreasing order on its main diagonal. The left singular vectors, that is, the columns of , 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 (see (6)) for the first 100 basis vectors is shown in Figure 6. In this problem, we chose the first POD basis vectors based on their truncation degree . The th order approximation of is given by the following truncated sequence: where and .
(a)
(b)
3.4.4. Construction of the Model for the First 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: and we replace by its thorder approximation in (32); the Galerkin projection states that the projection of on the space spanned by the basis functions vanishes, that is, 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 By evaluating the inner product in (34), 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 where , and . The initial condition for reads as .
For validating the reduced order model of the reactor, we applied constant input signals (K), ( K), and ( K) and constant perturbation signals ( mol/m^{3}) and ( K) to both the fullorder 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 reducedorder model, the averages of the absolute error for the temperature and concentration were calculated by means of the following formulas: where is the number of time steps and 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 mol/m^{3}. 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)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
The discretetime version of (36) that is used for designing the digital controller was obtained using the discretization method known as zeroorder hold (ZOH) with a sampling time of 0.2 s: where , , and are the matrices describing the new system. A modeling approach frequently adopted in model predictive controller (MPC) considers a discretetime statespace model in the incremental form [36]; hence (38) can be represented in the following form: where is the input increment, 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 openloop optimal control problem is solved online 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 online 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 closedloop system is nonlinear due to the presence of constraints. Along this work, we will deal with MPC controllers based on discretetime 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 loworder 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 discretetime statespace model in incremental form [43]:
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 discretetime statespace model as shown in (41). The cost of the infinite horizon MPC considered here can be defined as follows: 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 outputtracking 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]: 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: where With the terminal penalty, the cost defined in (42) reduces to
Finally, the control optimization problem of the infinite horizon MPC can be formulated as subject to
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: 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: Hence, the control objective defined in (49) becomes as follows:
Finally, the control optimization problem of the extended infinite horizon MPC can be formulated as follows: subject to
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 (), and the control system should keep the temperature inside the reactor below 335 K. In [46], two PODbased 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 of these POD coefficients can be calculated by where is the reference of the vector and is equal to zero (the model of the MPC is a discretetime 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 . Since the state vector is unknown and the changes in the concentration of the feed flow ( are not measured directly, they are estimated by means of an observer (in this case a Kalman filter) with the following formulation: where is the estimated vector of the POD coefficients, is the estimation , is the normalized temperature deviation of the feed flow , 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 is a selection matrix which selects the measured temperatures from the vector .
The control horizon was set to samples; and 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 , and . The Kalman gain matrix was computed from the following covariance matrices: , .
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 mol/m^{3} 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 mol/m^{3} 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)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
Furthermore, some quantities of interest are given in Table 2. In this table, is the maximum temperature reached inside the reactor during the test. is the percentage of the change of the mean steady state concentration at the reactor outlet with respect to its nominal value, that is, where is the mean nominal value (0.0339 mol/m^{3}) and is the mean concentration at the reactor outlet in steady state after the test.

In general, the control schemes showed a good behavior for rejecting the disturbances (typical magnitudes: mol/m^{3} and 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 welldefined 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 timedependent 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 () 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 nonlinear 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 openloop 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 tradeoff between complexity and performance; on the one hand it was possible to reduce the complexity of a highorder 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 (HDMPC). Contract no. INFSOICT223854 for funding this work.
References
 A. C. Antoulas, Approximation of LargeScale Dynamical Systems (Advances in Design and Control), vol. 6, Society for Industrial and Applied Mathematics (SIAM), 2005.
 K. Karhunen, “Zur spectral theorie stochastischer prozesse,” Annales Academiæ Scientiarum Fennicæ, vol. 36, pp. 1–7, 1946. View at: Google Scholar
 M. Loeve, “Fonctions aleatoires de second ordre,” Comptes Rendus De L'Académie Des Sciences, vol. 220, 1945. View at: Google Scholar
 J. L. Lumley, Stochastic Tools in Turbulence, vol. 12 of Applied mathematics and mechanics Mathematics in Science and Engineering, Academic Press, 1970.
 G. Berkooz, P. Holmes, and J. L. Lumley, “The proper orthogonal decomposition in the analysis of turbulent flows,” Fluid Mechanics, vol. 25, pp. 539–575, 1993. View at: Google Scholar
 F. Kwasniok, “The reduction of complex dynamical systems using principal interaction patterns,” Physica D, vol. 92, no. 12, pp. 28–60, 1996. View at: Google Scholar
 F. Kwasniok, “Optimal Galerkin approximations of partial differential equations using principal interaction patterns,” Physical Review E, vol. 55, no. 5, pp. 5365–5375, 1997. View at: Google Scholar
 P. Comon, “Independent component analysis: a new concept,” Signal Processing, vol. 36, no. 3, pp. 287–314, 1994. View at: Google Scholar
 P. Druault, J. Delville, and J. P. Bonnet, “Proper orthogonal decomposition of the mixing layer flow into coherent structures and turbulent Gaussian fluctuations,” Comptes Rendus, vol. 333, no. 11, pp. 824–829, 2005. View at: Publisher Site  Google Scholar
 S. Rahal, P. Cerisier, and H. Azuma, “Application of the proper orthogonal decomposition to turbulent convective flows in a simulated Czochralski system,” International Journal of Heat and Mass Transfer, vol. 51, no. 1718, pp. 4216–4227, 2008. View at: Publisher Site  Google Scholar
 G. Solari and F. Tubino, “A turbulence model based on principal components,” Probabilistic Engineering Mechanics, vol. 17, no. 4, pp. 327–335, 2002. View at: Publisher Site  Google Scholar
 M. V. Tabib and J. B. Joshi, “Analysis of dominant flow structures and their flow dynamics in chemical process equipment using snapshot proper orthogonal decomposition technique,” Chemical Engineering Science, vol. 63, no. 14, pp. 3695–3715, 2008. View at: Publisher Site  Google Scholar
 Y. Utturkar, B. Zhang, and W. Shyy, “Reducedorder description of fluid flow with moving boundaries by proper orthogonal decomposition,” International Journal of Heat and Fluid Flow, vol. 26, no. 2, pp. 276–288, 2005. View at: Publisher Site  Google Scholar
 M. Amabili, A. Sarkar, and M. P. Païdoussis, “Chaotic vibrations of circular cylindrical shells: galerkin versus reducedorder models via the proper orthogonal decomposition method,” Journal of Sound and Vibration, vol. 290, no. 3–5, pp. 736–762, 2006. View at: Publisher Site  Google Scholar
 U. Galvanetto and G. Violaris, “Numerical investigation of a new damage detection method based on proper orthogonal decomposition,” Mechanical Systems and Signal Processing, vol. 21, no. 3, pp. 1346–1361, 2007. View at: Publisher Site  Google Scholar
 P. B. Gonçalves, F. M. A. Silva, and Z. J. G. N. Del Prado, “Lowdimensional models for the nonlinear vibration analysis of cylindrical shells based on a perturbation procedure and proper orthogonal decomposition,” Journal of Sound and Vibration, vol. 315, no. 3, pp. 641–663, 2008, EUROMECH colloquium 483, Geometrically nonlinear vibrations of structures. View at: Publisher Site  Google Scholar
 M. Amabili, A. Sarkar, and M. P. Païdoussis, “Reducedorder models for nonlinear vibrations of cylindrical shells via the proper orthogonal decomposition method,” Journal of Fluids and Structures, vol. 18, no. 2, pp. 227–250, 2003, Axial and Internal Flow FluidStructure Interactions. View at: Publisher Site  Google Scholar
 B. F. Feeny and Y. Liang, “Interpreting proper orthogonal modes of randomly excited vibration systems,” Journal of Sound and Vibration, vol. 265, no. 5, pp. 953–966, 2003. View at: Publisher Site  Google Scholar
 M. Khalil, S. Adhikari, and A. Sarkar, “Linear system identification using proper orthogonal decomposition,” Mechanical Systems and Signal Processing, vol. 21, no. 8, pp. 3123–3145, 2007. View at: Publisher Site  Google Scholar
 X. Gilliam, J. P. Dunyak, D. A. Smith, and F. Wu, “Using projection pursuit and proper orthogonal decomposition to identify independent flow mechanisms,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 92, no. 1, pp. 53–69, 2004. View at: Publisher Site  Google Scholar
 T. Katayama, H. Kawauchi, and G. Picci, “Subspace identification of closed loop systems by the orthogonal decomposition method,” Automatica, vol. 41, no. 5, pp. 863–872, 2005. View at: Publisher Site  Google Scholar
 D. Chelidze and W. Zhou, “Smooth orthogonal decompositionbased vibration mode identification,” Journal of Sound and Vibration, vol. 292, no. 3–5, pp. 461–473, 2006. View at: Publisher Site  Google Scholar
 O. M. Agudelo, The application of proper orthogonal decomposition to the control of tubular reactors [Ph.D. thesis], Katholieke Universiteit Leuven, 2009.
 F. Leibfritz and S. Volkwein, “Reduced order output feedback control design for PDE systems using proper orthogonal decomposition and nonlinear semidefinite programming,” Linear Algebra and Its Applications, vol. 415, no. 23, pp. 542–575, 2006, Special Issue on Order Reduction of LargeScale Systems. View at: Publisher Site  Google Scholar
 D. Hömberg and S. Volkwein, “Control of laser surface hardening by a reducedorder approach using proper orthogonal decomposition,” Mathematical and Computer Modelling, vol. 38, no. 10, pp. 1003–1028, 2003. View at: Publisher Site  Google Scholar
 H. V. Ly and H. T. Tran, “Modeling and control of physical processes using proper orthogonal decomposition,” Mathematical and Computer Modelling, vol. 33, no. 1–3, pp. 223–236, 2001, Computation and control VI proceedings of the sixth Bozeman conference. View at: Publisher Site  Google Scholar
 J. A. Atwell and B. B. King, “Proper orthogonal decomposition for reduced basis feedback controllers for parabolic equations,” Mathematical and Computer Modelling, vol. 33, no. 1–3, pp. 1–19, 2001, Computation and control VI proceedings of the sixth Bozeman conference. View at: Publisher Site  Google Scholar
 R. Padhi and S. N. Balakrishnan, “Proper orthogonal decomposition based optimal neurocontrol synthesis of a chemical reactor process using approximate dynamic programming,” Neural Networks, vol. 16, no. 56, pp. 719–728, 2003, Advances in Neural Networks Research (IJCNN '03). View at: Publisher Site  Google Scholar
 C. Xu, Y. Ou, and E. Schuster, “Sequential linear quadratic control of bilinear parabolic PDEs based on POD model reduction,” Automatica, vol. 47, no. 2, pp. 418–426, 2011. View at: Publisher Site  Google Scholar
 S. S. Ravindran, “Control of flow separation over a forwardfacing step by model reduction,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 4142, pp. 4599–4617, 2002. View at: Publisher Site  Google Scholar
 S. S. Ravindran, “Optimal boundary feedback flow stabilization by model reduction,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 2528, pp. 2555–2569, 2007. View at: Publisher Site  Google Scholar
 W. Xie, I. Bonis, and C. Theodoropoulos, “Offline model reduction for online linear MPC of nonlinear largescale distributed systems,” Computers and Chemical Engineering, vol. 35, no. 5, pp. 750–757, 2011. View at: Google Scholar
 P. Astrid, Reduction of process simulation models: a proper orthogonal decomposition approach [Ph.D. thesis], Technishche Universiteit Eindhoven, 2004.
 S. Fogler, Elements of Chemical Reaction Engineering, Prentice Hall, Boston, Mass, USA, 4th edition edition, 2008.
 O. M. Agudelo, J. J. Espinosa, and B. De Moor, “Control of a tubular chemical reactor by means of POD and predictive control techniques,” in Proceedings of the European Control Conference (ECC '07), vol. 20, pp. 1046–1053. View at: Google Scholar
 M. A. Rodrigues and D. Odloak, “MPC for stable linear systems with model uncertainty,” Automatica, vol. 39, no. 4, pp. 569–583, 2003. View at: Publisher Site  Google Scholar
 D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert, “Constrained model predictive control: stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. View at: Publisher Site  Google Scholar
 S. J. Qin and A. B. Badgwell, “A survey of industrial model predictive control technology,” Control Engineering Practice, vol. 11, no. 7, pp. 733–764, 2003. View at: Google Scholar
 C. E. García, D. M. Prett, and M. Morari, “Model predictive control: theory and practice: a survey,” Automatica, vol. 25, no. 3, pp. 335–348, 1989. View at: Google Scholar
 D. Panagiotis and D. Prodromus, “Nonlinear control of diffusionconvectionreaction processes,” Computers and Chemical Engineering, vol. 20, Supplement 2, pp. S1071–S1076, 1996. View at: Google Scholar
 M. Li and D. Panagiotis, “Optimal transition control of diffusionconvectionreaction processes,” in Proceedings of the 8th International IFAC Symposium on Dynamics and Control of Process System, Cancun, Mexico, 2007. View at: Google Scholar
 Y. Ou and E. Schuster, “Model predictive control of parabolic PDE systems with dirichlet boundary conditions via galerkin model reduction,” in Proceedings of the American Control Conference (ACC '09), pp. 1–7, June 2009. View at: Google Scholar
 J. H. Lee, M. Morari, and C. E. Garcia, “Statespace interpretation of model predictive control,” Automatica, vol. 30, no. 4, pp. 707–717, 1994. View at: Publisher Site  Google Scholar
 K. R. Muske and J. B. Rawlings, “Model predictive control with linear models,” AIChE Journal, vol. 39, no. 2, pp. 262–287, 1993. View at: Google Scholar
 D. Odloak, “Extended robust model predictive control,” AIChE Journal, vol. 50, no. 8, pp. 1824–1836, 2004. View at: Publisher Site  Google Scholar
 A. Marquez, J. J. Espinosa, and D. Odloak, “IHMPC and POD to the control of a nonisothermal tubular reactor,” in Proceedings of the 9th International Symposium on Dynamics and Control of Process Systems (DYCOPS '10), pp. 431–436, 2010. View at: Google Scholar
Copyright
Copyright © 2013 Alejandro Marquez et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.