# Model reduction using proper orthogonal decomposition and predictive control of distributed reactor system.

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.

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 x(t) [member of] [R.sup.N] = [[[x.sub.1](t),[x.sub.2](t), ..., [x.sub.N](t)].sup.T] be the state vector of a given dynamical system, and let X [member of] [R.sup.NxNd] with [N.sub.d] [greater than or equal to] N be the so-called snapshot matrix that contains a finite number of samples or snapshots of the evolution of x(t) at [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] 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) [[phi].sub.j] [member of] [R.sup.N], for all j = 1,2, ..., N:

x([t.sub.i]) = [N.summation over (j=1)] [a.sub.j] ([t.sub.i])[[phi].sub.j], [for all]i = 1,2, ..., [N.sub.d], [a.sub.j] ([t.sub.i]) = < x([t.sub.i]), [[phi].sub.j]>, [for all]j = 1,2, ..., N, (1)

where [a.sub.j]([t.sub.i]) is the coordinate of x([t.sub.i]) with respect to the basis vector [[phi].sub.j] (it is also called time-varying coefficient or POD coefficient) and <*,*> denotes the Euclidean inner product. Since the first n most relevant basis vectors capture most of the energy in the collected data set, we can construct an nth-order approximation of the snapshots by means of the following truncated sequence:

x([t.sub.i]) = [n.summation over (j=1)] [a.sub.j] ([t.sub.i])[[phi].sub.j], [for all]i = 1,2, ..., [N.sub.d], n [much less than] N. (2)

In POD, the orthonormal basis vectors are calculated in such a way that the reconstruction of the snapshots using the first n most relevant basis vectors is optimal in the sense that the sum-squared-error (SSE) between x([t.sub.i]) and [x.sub.N]([t.sub.i]), for all i=l, ..., [N.sub.d],

[E.sub.n] = [[N.sub.d].summation over (1=1)] [parallel]x([t.sub.i]) - [x.sub.n]([t.sub.i])[[parallel].sup.2.sub.2] (3)

is minimized. Herein [parallel]*[[parallel].sub.2] denotes the [L.sub.2]-norm or the euclidean Norm. In other words, the POD basis vectors are the ones that solve the following constrained optimization problem:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

subject to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (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 ([X.sub.snap]). The singular values of [X.sub.snap] are positive real numbers that are ordered in a decreasing way, [[sigma].sub.1] [greater than or equal to] [[sigma].sub.2] *** [greater than or equal to] [[sigma].sub.n] [greater than or equal to] 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 n most relevant basis vectors is certainly of central importance. A criterion commonly used for choosing n 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 [X.sub.snap]

[P.sub.n] = [[summation].sup.n.sup.j=1][[sigma].sup.2.sub.j]/[[summation].sup.N.sup.j=1][[sigma].sup.2.sub.j], n = 1, ..., (N). (6)

The ratio [P.sub.n] 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 n has to be determined for [P.sub.n] = 0.99.The closer [P.sub.n] to 1, or similarly the closer 1 - [P.sub.n] to 0, the better the approximation of X.

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:

[??](t) = f(x(t), u(t)), y(t) = g(x(t), u(t)). (7)

Let us define a residual function R([??], x) for (7) as follows:

R([??], x) = [??](t)-f(x(t),u(t)), (8)

and let R([[??].sub.n],[x.sub.n]) be the residual when the state vector x(t) is approximated by its nth-order approximation:

[x.sub.n](t) = [n.summation over (j=1)] [a.sub.j](t) [[phi].sub.j] = [[PHI].sub.N]a(t), n [much less than] N, (9)

where [[PHI].sub.n] = [[[phi].sub.1],[[phi].sub.2], ..., [[phi].sub.n]] and a(t) = [[[a.sub.1](t),[a.sub.2](t), ..., [a.sub.n](t)].sup.T]. In the Galerkin projection, the projection of the residual R([??], x) on the space spanned by the basis vectors [[PHI].sub.n] vanishes, that is,

< R([??],x),[[phi].sub.j]> = 0, [for all]j = 1, ... n, (10)

where <*,*> denotes the Euclidean inner product. Replacing x(t) by its nth-order approximation [x.sub.n](t) =[[phi].sub.n] a(t) in (7):

[[PHI].sub.n][??](t) = f ([[PHI].sub.n]a(t),u(t)), (11)

and then we apply the inner product criterion (10) as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

and given that [[PHI].sup.T.sub.n] [[PHI].sub.n] = [I.sub.n] because of the orthonormality of the basis vectors, we have the model for the POD coefficients reduce to

[??](t) = [[PHI].sup.T.sub.n] f([[PHI].sub.n]a(t),u(t)). (13)

Finally, the reduced order model of (10) with only n states has the following form,

[??](t) = [[PHI].sup.T.sub.n]f([[PHI].sub.n]a(t), u(t)), y(t) = g([[PHI].sub.n] a(t), u(t)). (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 (A + B [??] 2C). The reactor is surrounded by 3 cooling/heating jackets as it is shown in Figure 2. The temperature of the jackets fluids ([T.sub.J1], [T.sub.J2] and [T.sub.J3]) 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 A and the energy balance on the differential cylinder shown in Figure 3 produce the following equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (15)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (16)

with the following boundary conditions:

(1) at r = 0, we have symmetry [partial derivative]T/[partial derivative]r = 0 and [partial derivative]C/[partial derivative]r = 0,

(2) at r = R, 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,

- [K.sub.e] [[partial derivative]T/[partial derivative]r][|.sub.R] = [H.sub.w](T(R,z) - [T.sub.w]), (17)

(3) at r = R, there is no mass flow through the tube walls [partial derivative]C/[partial derivative]r = 0;

(ii) axial:

(1) T = [T.sub.in] and C = [C.sub.in] at z = 0,

(2) at the outlet of the reactor z = L, [partial derivative]T/[partial derivative]z = 0 and [partial derivative]C/[partial derivative]z = 0.

Here, C(r,z,t) is the reactant concentration in [mol x [m.sup.-3]], T(r, z, t) is the reactant temperature in [K], [D.sub.e] is the diffusivity of all species in [[m.sup.2] x [s.sup.-1] ], [K.sub.e] is the thermal conductivity of the reaction mixture in [J x [m.sup.-1] x [s.sup.-1] x [K.sup.-1]], [K.sub.eq0] is the equilibrium constant at 300 K, [U.sub.z] is the fluid velocity in [m x [s.sup.-1]], (-[DELTA][H.sub.rx]) is the heat of the reaction in [J x [mol.sup.-1]], [rho] and [C.sub.p] are the density in [kg x [m.sup.-3]] and the specific heat in [J x [kg.sup.-1] x [K.sup.-1]] of the mix, respectively, [k.sub.0] is the rate constant in [[m.sup.6] x [mol.sup.-1] x [s.sup.-1] x [kg.sup.-1]], E is the activation energy in [J x [mol.sup.-1]], R is the ideal gas constant in [J x [mol.sup.-1] x [K.sup.-1]], [H.sub.w] is the heat transfer coefficient in [J x [m.sup.-2] x [s.sup.-1] x [K.sup.-1]], L is the reactor length in [m], [C.sub.in] and [T.sub.in] are the concentration in [mol x [m.sup.-3]] and the temperature in [K] of the feed flow, z is the axial coordinate in [m], r is the radial coordinate in [m], t is the time in [s], and [T.sub.w](z, t) is the reactor wall temperature in [K] which is defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (18)

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

The temperature of the jacket sections [T.sub.j1], [T.sub.j2], and [T.sub.j3] 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 [+ or -] 10 K for the temperature and [+ or -] 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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (19)

with T(r,0) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN is a function that depends on the boundary conditions when r = R. The discrete version of (19) can be found by replacing the spatial derivatives by first-order upwind and backward differences approximations as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (21)

where n and m are the number of sections in axial and radial directions in which the reactor is divided, [z.sub.a] and [z.sub.b] are the reactor sections defining the ending of the first and second jacket, respectively, [T.sub.f] and [C.sub.f] are normalization factors, [[bar.C].sub.i,j] = [C.sub.i,j]/[C.sub.f] and [[bar.T].sub.i,j] = [T.sub.ij]/[T.sub.f] are the normalized concentration and temperature of the i, jth section of the reactor, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the normalized reactor wall temperature of the i, jth section, and [DELTA]z and [DELTA]r 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (22)

subject to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (23)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the desired concentration (normalized) at the reactor output, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the desired temperature (normalized) inside the reactor of the i, jth increment, [[bar.C].sub.i,n] is the concentration (normalized) at the reactor output, [omega] is a tradeoff coefficient, [T.sub.Jmin] and [T.sub.Jmax] are the lower and upper temperature values of the fluids of the jackets, and [T.sub.max] 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 [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] was selected equal to the normalized temperature of the feed flow ([[bar.T].sub.in] = [T.sub.in]/[T.sub.f]) for i = 1,2, ..., m and j = 1,2, ..., n. The trade-off parameter [omega] can take values between 0 and 1. When [omega] 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 w 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: n = 30, m = 30, [T.sub.f] = 320 K, [C.sub.f] = 500 mol/[m.sup.3], [T.sub.Jmin] = 280 K, [T.sub.Jmax] = 340 K, [T.sub.max] = 335 K, Tol = [10.sup.-3], and [omega] = 0.7. The maximum allowed temperature ([T.sub.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 [omega] 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 [T.sub.j1] = 291.1705 K, [T.sub.j2] = 293.6938 K, and [T.sub.J3] = 294.8280 K. The optimal concentration and temperature profiles can be observed in Figures 4 and 5.

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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (24)

with the following boundary conditions:

(1) at r = 0, we haves symmetry [partial derivative][T.sup.[DELTA]]/[partial derivative]r = 0 and [partial derivative][C.sup.[DELTA]]/[partial derivative]r = 0,

(2) at r = R,

-[K.sub.e] [[partial derivative][T.sup.[DELTA]]/[partial derivative]r] [|.sub.R] = [H.sub.w] ([T.sup.[DELTA]] (R,z) - [T.sup.[DELTA].sub.w]), (25)

(3) at r = R; [partial derivative][C.sup.[DELTA]]/[partial derivative]r = 0,

(ii) axial:

(1) [T.sup.[DELTA]] = [T.sup.[DELTA].sub.in] and [C.sup.[DELTA]] = [T.sup.[DELTA].sub.in] at = 0,

(2) at the outlet of the reactor z = L, [partial derivative][T.sup.[DELTA]]/[partial derivative]z = 0 and [partial derivative][C.sup.[DELTA]]/[partial derivative]z = 0.

Here, [C.sup.[DELTA]] = C - [C.sup.*] , [T.sup.[DELTA]] = T - T*,and [T.sup.[DELTA].sub.w] = [T.sub.w] -[T.sup.*.sub.w] are the deviations from the steady state of the concentration, temperature, and reactor wall temperature; [C.sup.*], [T.sup.*] and [T.sup.*.sub.w] 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,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

then (24) can be cast as follows:

[??](t) = Ax(t) + Bu(t) + Fd(t), (27)

where A [member of] [R.sup.2mxnx2mxn,] b [member of] [R.sup.2mxnx3], and F [member of] [R.sup.2mxnx2] are the matrices describing the system, x(t) [member of] [R.sup.2mxn] is the state vector, u(t) [member of] [R.sup.3] is the vector of the inputs, and d(t) [member of] [R.sup.2] 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 n x m = 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 ([X.sub.snap] [member of] [R.sup.1800x10000]) when independent step changes were made in the input u(t) and perturbation d(t) signals on the nonlinear model (15):

[X.sub.snap] = [x(t = [DELTA]t), x(t = 2[DELTA]t), ..., x(t= 1000[DELTA]t)]. (28)

Along the simulations, 10000 samples were collected using a sampling time [DELTA]t 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 [X.sub.snap]:

[X.sub.snap] = [PHI][SIGMA][[PSI].sup.T], (29)

where [PHI] [member of] [R.sup.1800x1800] and [PSI] [member of] [R.sup.10000x10000] are unitary matrices and [SIGMA] [member of] [R.sup.1800x10000] is a matrix that contains the singular values of [X.sub.snap] in a decreasing order on its main diagonal. The left singular vectors, that is, the columns of [PHI],

[PHI] = [[[phi].sub.1], [[phi].sub.2], ..., [[phi].sub.1800]], (30)

are the POD basis vectors.

3.4.3. Selection of the Most Relevant POD Basis Vectors. The n most relevant POD basis vectors are chosen using the energy criterion presented in Section 2.1.The plot of 1 - [P.sub.N] (see (6)) for the first 100 basis vectors is shown in Figure 6.In this problem, we chose the first n=50 POD basis vectors based on their truncation degree 1 - [P.sub.n] = 3.3 x [10.sup.-4] ([P.sub.n] = 0.9996). The 50th order approximation of x(t) is given by the following truncated sequence:

[x.sub.n](t) = [50.summation over (j=1)] [a.sub.j](t)[[phi].sub.j] = [PHI].sub.N]a(t), (31)

where [[PHI].sub.n] = [[[phi].sub.1], [[phi].sub.2], ..., [[phi].sub.50]] and a(t) = [[a.sub.1](t),[a.sub.2](t), [a.sub.50](t)].sup.T].

3.4.4. Construction of the Model for the First n=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 R(x) for (27) as follows:

R ([??],x) = [??](t) - Ax(t) - Bu(t) - Fd(t), (32)

and we replace x(t) by its nth-order approximation [x.sub.n](t) = [[PHI].sub.n]a(t) in (32); the Galerkin projection states that the projection of R([x.sub.n]) on the space spanned by the basis functions [[PHI].sub.n] vanishes, that is,

<(R([x.sub.n]),[[phi].sub.j])> = 0, j = 1, ..., n, (33)

where <*,*> denotes the inner product. Replacing x(t) by its nth order approximation [x.sub.n](t) = [[PHI].sub.n]a(t) in (27) and applying the inner product criterion (Galerkin projection) to the resulting equation, we have

([[PHI].sub.n][??](t), [[phi].sub.j]) = (A[[PHI].sub.n]a(t) + Bu(t) + Fd(t)[[phi].sub.j]). (34)

By evaluating the inner product in (34),

[??](t) = [[PHI].sup.T.sub.n]A[[PHI].sub.n]a(t) + [[PHI].sup.T.sub.n]Bu(t) + [[PHI].sup.T.sub.n]Fd(t), (35)

And we obtain the model for the first n POD coefficients. The reduced order model of the reactor with only 50 states is then given by

[??](t) = [A.sub.r]a(t) + [B.sub.r]u(t) + [F.sub.r]d(t), [x.sub.n](t) = [[PHI].sub.n]a(t), (36)

where [A.sub.r] =[[PHI].sup.T.sub.n]A[[PHI].sub.n], [B.sub.r] = [[PHI].sup.T.sub.n]B, and [F.sub.r] = [[PHI].sup.T.sub.n]F. The initial condition for a(t) reads as a(0) = 0.

For validating the reduced order model of the reactor, we applied constant input signals [[bar.T].sup.[DELTA].sub.J1](t) = 0.0312 ([[bar.T].sup.[DELTA].sub.J1](t) = 10 K), [[bar.T].sup.[DELTA].sub.J2](t) = -0.0312 ([T.sup.[DELTA].sub.J2](t) = 10 K), and [[bar.T].sup.[DELTA].sub.J3] (t) = 0.0625 ([T.sup.[DELTA].sub.J3](t) = 20 K) and constant perturbation signals [[bar.C].sup.[DELTA].sub.in](t) = 2 x [10.sup.-4] ([C.sup.[DELTA].sub.in](t) = 0.1 mol/[m.sup.3]) and [[bar.T].sup.[DELTA].sub.in](t) = 0.0156 ([T.sup.[DELTA].sub.in](t) = 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 ([E.sub.T]) and concentration ([E.sub.C]) were calculated by means of the following formulas:

[E.sub.T] = 1/[N.sub.s] [[N.sub.s].summation over (k=1)] [absolute value of (K[DELTA]t) - [T.sub.n](K[DELTA]t)], [E.sub.C] = 1/[N.sub.s] [[N.sub.s].summation over (k=1)] [absolute value of (K[DELTA]t) - [C.sub.n](K[DELTA]t)], (37)

where [N.sub.s] = 10000 is the number of time steps and [DELTA]t = 1 s. The plots of [E.sub.T] and [E.sub.C] are shown in Figure 10. For the Temperature profile, the maximum value for the error [E.sub.T] is 0.7542 K. For the concentration, the maximum peak for the error [E.sub.C] is 3.3 x [10.sup.-3] mol/[m.sup.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.

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (39)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

where [x.sup.d](k) = [V.sub.1]a(k), [x.sup.s](k) = [V.sub.2]a(k - 1), [DELTA][u.sub.k] = [u.sub.k] [u.sub.k-1] is the input increment, [DELTA][d.sub.k] = [d.sub.k] - d[k.sub.-1] is the disturbance increment, and [V.sub.1], [V.sub.2] are transformation matrices. In the state equation defined in (39), the state component [x.sup.s] corresponds to the integrating poles produced by the incremental form of the model, and [x.sup.d](k) = a(k) corresponds to the system modes. For stable systems, it is easy to show that when the system approaches steady state, component [x.sup.d] 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]:

x(k + 1) = Ax(k) + B[DELTA]u (k) + F[DELTA]d (k), y (k) = Cx(k). (41)

It is clear that the model defined in (39) and (40) is in the form defined in (41) if the state x and model matrices A, B, F, and C 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:

[J.sub.k,[infinity]] = [[infinity].summation over (j=1)]e[(k + j).sup.T]Qe(k + j) + [p-1.summation over (j=0)] [DELTA]u[(k + j).sup.T]R [DELTA]u(k + j), (42)

where e(k + j) = y(k + j) - r(j), y(k + j) is the output prediction at time instant k + j made at time k, r is the desired output reference, p is the control horizon, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] 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 [bar.Q]. For the cost defined in (42), such a terminal penalty is computed by the following Lyapunov equation [44]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (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:

[[bar.C].sup.s][x.sup.s](k) - [y.sub.sp] + [[bar.C].sup.s][[??].sup.s][DELTA][u.sub.k] = 0, (44)

where

[[??].sup.s] = [[D.sup.S] *** [D.sup.s]], [[bar.C].sup.s] = diag [[C.sup.s] *** [C.sup.s]], [DELTA][u.sub.k] = [[DELTA]u[(k).sup.T] [DELTA]u[(k + 1).sup.T] *** [DELTA]u[(k + p - 1).sup.T]]. (45)

With the terminal penalty, the cost defined in (42) reduces to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (46)

Finally, the control optimization problem of the infinite horizon MPC can be formulated as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (47)

subject to

model in incremental form (39)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (48)

For large changes on [x.sup.s](k) or [y.sub.sp] or if [y.sub.sp] 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (49)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a vector of slack variables and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is assumed positive definite. Observe that each slack variable refers to a given controlled output. Weight matrix S 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 [bar.Q] 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:

[[bar.C].sup.s] [x.sup.s](k) - [y.sub.sp] + [[bar.C].sup.s][[??].sup.s][DELTA][u.sub.k] - [[DELTA].sub.k] = 0. (50)

Hence, the control objective defined in (49) becomes as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (51)

Finally, the control optimization problem of the extended infinite horizon MPC can be formulated as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (52)

subject to

model in incremental form (39)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (53)

Theorem 1. For a stable system, if in the control objective defined in (51) the slack weight matrix S 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 (280 K [less than or equal to] [T.sub.J1](t),[T.sub.J2](t), [T.sub.J3](t) [less than or equal to] 335 K), 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 ([a.sub.ref]) of these POD coefficients can be calculated by

[a.sub.ref] = [[PHI].sup.T.sub.n] [x.sub.ref], (54)

where [x.sub.ref] is the reference of the vector x(t) 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 [C.sup.d] = [I.sub.50x50][V.sub.1], [C.sup.s] = [0.sub.50x50][V.sub.2]. Since the state vector a(k) is unknown and the changes in the concentration of the feed flow (d1(k) = [C.sup.[DELTA].sub.in](k)) are not measured directly, they are estimated by means of an observer (in this case a Kalman filter) with the following formulation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (55)

where [??] is the estimated vector of the POD coefficients, [[??].sub.1](k) is the estimation [[bar.C].sup.[DELTA].sub.in], [d.sub.2](k) is the normalized temperature deviation of the feed flow [[bar.C].sup.[DELTA].sub.in],(k), y(k) [member of] [R.sup.10] is a vector containing ten temperature measurements (normalized deviations) along the reactor, [??](k) is the estimate of Y(k), [L.sub.a] and [L.sub.d] are the submatrices of the observer gain (Kalman gain), [F.sub.C] and [F.sub.T] are the column vectors of [??] = [[[??].sub.C], [[??].sub.T]], and [C.sub.sel] is a selection matrix which selects the measured temperatures from the vector [x.sub.n](k).

The control horizon m was set to 10 samples; [u.sub.min] and [u.sub.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 Q = 1 x [I.sub.50x50], R = 10 x [I.sub.3x3], and S = 1 x [I.sub.50x50]. The Kalman gain matrix was computed from the following covariance matrices: [R.sub.w] = 1 x [I.sub.51x51], [R.sub.v] = [10.sup.-3] x [I.sub.10x10].

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 x [10.sup.-3] mol/[m.sup.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 25 x [10.sup.-3] mol/[m.sup.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.

Furthermore, some quantities of interest are given in Table 2. In this table, [T.sub.max] is the maximum temperature reached inside the reactor during the test. [DELTA][C.sub.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,

[DELTA][C.sub.out]% = [[C.sub.N] - [C.sup.*.sub.N]]/[C.sup.*.sub.N] (56)

where [C.sup.*.sub.N] is the mean nominal value (0.0339 mol/[m.sup.3]) and [C.sub.N] 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: [C.sub.in] = [+ or -] 25 x [10.sup.-3] mol/[m.sup.3] and [T.sub.in] = [+ or -] 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 ([X.sub.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.

http://dx.doi.org/10.1155/2013/763165

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.

References

[1] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems (Advances in Design and Control), vol. 6, Society for Industrial and Applied Mathematics (SIAM), 2005.

[2] K. Karhunen, "Zur spectral theorie stochastischer prozesse," Annales Academice Scientiarum Fennicce, vol. 36, pp. 1-7, 1946.

[3] M. Loeve, "Fonctions aleatoires de second ordre," Comptes Rendus De L'Academie Des Sciences, vol. 220, 1945.

[4] J. L. Lumley, Stochastic Tools in Turbulence, vol. 12 of Applied mathematics and mechanics Mathematics in Science and Engineering, Academic Press, 1970.

[5] 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.

[6] F. Kwasniok, "The reduction of complex dynamical systems using principal interaction patterns," Physica D, vol. 92, no. 12, pp. 28-60, 1996.

[7] F. Kwasniok, "Optimal Galerkin approximations of partial differential equations using principal interaction patterns," Physical Review E, vol. 55, no. 5, pp. 5365-5375, 1997.

[8] P. Comon, "Independent component analysis: a new concept," Signal Processing, vol. 36, no. 3, pp. 287-314, 1994.

[9] 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.

[10] 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. 17-18, pp. 4216-4227, 2008.

[11] G. Solari and F. Tubino, "A turbulence model based on principal components," Probabilistic Engineering Mechanics, vol. 17, no. 4, pp. 327-335, 2002.

[12] 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.

[13] Y. Utturkar, B. Zhang, and W. Shyy, "Reduced-order 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.

[14] M. Amabili, A. Sarkar, and M. P. Paidoussis, "Chaotic vibrations of circular cylindrical shells: galerkin versus reduced-order models via the proper orthogonal decomposition method," Journal of Sound and Vibration, vol. 290, no. 3-5, pp. 736-762, 2006.

[15] 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.

[16] P. B. Goncalves, F. M. A. Silva, and Z. J. G. N. Del Prado, "Low-dimensional 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 non-linear vibrations of structures.

[17] M. Amabili, A. Sarkar, and M. P. Paidoussis, "Reduced-order 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 Fluid-Structure Interactions.

[18] 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.

[19] 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.

[20] 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.

[21] 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.

[22] D. Chelidze and W. Zhou, "Smooth orthogonal decomposition-based vibration mode identification," Journal of Sound and Vibration, vol. 292, no. 3-5, pp. 461-473, 2006.

[23] O. M. Agudelo, The application of proper orthogonal decomposition to the control of tubular reactors [Ph.D. thesis], Katholieke Universiteit Leuven, 2009.

[24] 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. 2-3, pp. 542-575, 2006, Special Issue on Order Reduction of Large-Scale Systems.

[25] D. Homberg and S. Volkwein, "Control of laser surface hardening by a reduced-order approach using proper orthogonal decomposition," Mathematical and Computer Modelling, vol. 38, no. 10, pp. 1003-1028, 2003.

[26] 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.

[27] 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.

[28] R. Padhi and S. N. Balakrishnan, "Proper orthogonal decomposition based optimal neurocontrol synthesis of a chemical Journal of Control Science and Engineering 19 reactor process using approximate dynamic programming," Neural Networks, vol. 16, no. 5-6, pp. 719-728, 2003, Advances in Neural Networks Research (IJCNN '03).

[29] 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.

[30] S. S. Ravindran, "Control of flow separation over a forward-facing step by model reduction," Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 41-42, pp. 4599-4617, 2002.

[31] S. S. Ravindran, "Optimal boundary feedback flow stabilization by model reduction," Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 25-28, pp. 2555-2569, 2007.

[32] W. Xie, I. Bonis, and C. Theodoropoulos, "Off-line model reduction for on-line linear MPC of nonlinear large-scale distributed systems," Computers and Chemical Engineering, vol. 35, no. 5, pp. 750-757, 2011.

[33] P. Astrid, Reduction of process simulation models: a proper orthogonal decomposition approach [Ph.D. thesis], Technishche Universiteit Eindhoven, 2004.

[34] S. Fogler, Elements of Chemical Reaction Engineering, Prentice Hall, Boston, Mass, USA, 4th edition edition, 2008.

[35] O. M. Agudelo, J. J. Espinosa, and B. DeMoor, "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.

[36] M. A. Rodrigues and D. Odloak, "MPC for stable linear systems with model uncertainty," Automatica, vol. 39, no. 4, pp. 569-583, 2003.

[37] 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.

[38] 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.

[39] C. E. Garcia, D.M. Prett, and M. Morari, "Model predictive control: theory and practice: a survey," Automatica, vol. 25, no. 3, pp. 335-348, 1989.

[40] D. Panagiotis and D. Prodromus, "Nonlinear control of diffusion-convection-reaction processes," Computers and Chemical Engineering, vol. 20, Supplement 2, pp. S1071-S1076, 1996.

[41] M. Li and D. Panagiotis, "Optimal transition control of diffusion-convection-reaction processes," in Proceedings of the 8th International IFAC Symposium on Dynamics and Control of Process System, Cancun, Mexico, 2007.

[42] 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.

[43] J. H. Lee, M. Morari, and C. E. Garcia, "State-space interpretation of model predictive control," Automatica, vol. 30, no. 4, pp. 707-717, 1994.

[44] K. R. Muske and J. B. Rawlings, "Model predictive control with linear models," AIChE Journal, vol. 39, no. 2, pp. 262-287, 1993.

[45] D. Odloak, "Extended robust model predictive control," AIChE Journal, vol. 50, no. 8, pp. 1824-1836, 2004.

[46] A. Marquez, J. J. Espinosa, and D. Odloak, "IHMPC and POD to the control of a non-isothermal tubular reactor," in Proceedings of the 9th International Symposium on Dynamics and Control of Process Systems (DYCOPS '10), pp. 431-436, 2010.

Alejandro Marquez, (1) Jairo Jose Espinosa Oviedo, (1) and Darci Odloak (2)

(1) Faculty of Minas, National University of Colombia, 050041 Medellin, Colombia

(2) Chemical Engineering Department, University of Sao Paulo, 05508-900 Sao Paulo, SP, Brazil

Correspondence should be addressed to Alejandro Marquez; amarque@unal.edu.co

Received 28 November 2012; Revised 28 February 2013; Accepted 15 March 2013

```
TABLE 1: Values of the reactor parameters with axial and radial
diffusion, reaction, and convection.

Parameter           Value

[D.sub.e]           [10.sup.-9] [m.sup.2] x [s.sup.-1]
[v.sub.0]           [10.sup.-5] [m.sup.3] x [s.sup.-1]
L                   10 m
R                   0.05 m
[[rho].sub.cat]     1500 kg x [m.sup.-3]
[rho]               1000 kg x [m.sup.-3]
[C.sub.p]           4180 J x [kg.sup.-1] x [K.sup.-1]
[DELTA][H.sub.rx]   -83680 J x [mol.sup.-1]
[K.sub.eq0]         1000
[K.sub.e]           0.559 J x [m.sup.-1] x [s.sup.-1] x [K.sup.-1]
[k.sub.0]           1.1 x [10.sup.8] [m.sup.6] x [mol.sup.-1] x
[s.sup.-1] x [kg.sup.-1]
E                   95238 J x [mol.sup.-1]
R                   8.314 J x [mol.sup.-1] x [K.sup.-1]
[C.sub.in]          500 mol x [m.sup.-3]
[T.sub.in]          320 K
[H.sub.w]           1300 J x [m.sup.-2] x [s.sup.-1] x [K.sup.-1]
[H.sub.r]           0.2 [s.sup.-1]

TABLE 2: Performance parameters of the control systems.

Quantities                Test 1     Test 2

[T.sub.max] [k]          332.3218   323.6579
[DELTA][C.sub.out] [%]    -1.0531    0.8732
```