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

Authors: Alejandro Marquez; Jairo José Espinosa Oviedo; Darci Odloak
Journal: Journal of Control Science and Engineering
 (2013)

---

## Abstract

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


---

## Body

## 1. Introduction

Advances in power computing have propelled the development of process models increasingly detailed and precise, which are then used in the design, optimization, monitoring, and diagnosis of faults, among other tasks. Usually distributed chemical reactors are described by partial differential equations (PDEs) that show the space-time evolution of some variables of interest. In order to simulate PDEs, the spatial domain is discretized, obtaining a large number of ordinary differential equations (ODEs). However, a fine discretization leads to an increase in model complexity. To reduce the complexity of the models, a technique based on orthogonal decomposition of a set of measurements of physical quantities (such as temperature and concentration) is used to represent these quantities along the space and time. This technique, proper orthogonal decomposition (POD), has been used to reduce the order of a large number of systems. This method is based on orthonormal basis functions generated from process data (snapshot matrix) which are obtained by simulation or experimentation on the process; These data are taken by excitation of the process through manipulated variables, external inputs and disturbances of the process. The main idea is to consider POD basis functions which capture the spatial dynamics of our system. The advantage of working with these basis functions is that it is possible to reduce the model order from hundreds or thousands to a few tens. This reduction, resulting in ease of simulation, assimilation and optimization, enables that such models can be applied in real time applications.There are several methods available for reducing the dimension of a system. The most immediate is a heuristic approach that consists of proposing a priori solution to the equations of motion on the grounds of symmetry and boundary conditions. These solutions usually take the form of a truncated series in terms of general sets of orthogonal functions, such as Fourier modes or spherical harmonics. Antoulas [1] proposes a classification into two main groups, namely, Krylov and singular value decomposition (SVD) methods. Krylov methods make use of iterations for finding approximations to large-scale dynamical systems. The SVD methods are based on the decomposition of the state vector into a set of vectors that can be ordered in a certain sense. These methods include balanced truncation and Hankel approximations for linear systems and the proper orthogonal decomposition (POD) methods and empirical grammians for nonlinear systems. In this work, we are concerned with the POD and its combination with Galerkin projection to produce reduced linear dynamical versions of the original large-scale system (the methodology used in this work is shown in Figure 1). The interested reader is referred to [1] for a more complete account of the available reduction methods. The POD [2, 3] is a statistical technique to extract features from a given dataset by searching for patterns that optimally represent the dataset with respect to quantities such as variance or energy. The output of the POD is a set of time-independent orthogonal functions called empirical orthogonal functions (EOFs). Each EOF is associated with a certain amount of variance or energy. The first suggestion to use the POD in the analysis of dynamical systems was due to Lumley [4]. One of the first phenomena to be studied by means of reduced models arising from the use of the POD was turbulence [5]. Perhaps for nonlinear systems the most studied method is the POD in conjunction with Galerkin projections [1]. However, this is not the only available method, and some other ideas have been proposed, such as a decomposition into principal interaction patterns (PIPs), which explicitly incorporates information about the system’s dynamics within the determination of the patterns [6, 7]. There are other techniques that might be suitable for developing low-order models such as the independent component analysis (ICA) [8]. The ICA aims at answering questions such as, what signal comes from what source? This problem is known as the cocktail-party problem. The solution relies on the assumption of statistical independence of the sources and, therefore, of the signals. The outcome is a set of modes and a mixing matrix that sets the relationship between the modes in order to reconstruct the original mixed signal. The ICA modes are statistically independent but not necessarily orthogonal, and there is no specific order. Furthermore, they are not clearly related to any physical quantity that can help to decide what modes to retain and what to rule out. These features may constitute a disadvantage when trying to construct low-order models since a truncation using these modes would lack an appropriate parameter to measure the expected accuracy of the resulting model. Original algorithms assume linearity although some efforts have been made to extend the formalism to nonlinear systems. So far these ideas have not yet been fully investigated in the context of dynamical systems. Beside this introduction section, this paper has four other sections. Section 2 describes the model reduction method that is used here. Section 3 describes the model of the nonisothermal tubular reactor that is used to illustrate the application of the proposed reduced and control methods. Section 4 addresses the control of a nonisothermal tubular reactor using a reduced model and MPC of infinite horizon. Finally, Section 5 concludes the paper.Figure 1
Model reduction framework.


## 2. Proper Orthogonal Decomposition (POD) and Galerkin Projection

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


### 2.1. General Procedure

Letx(t)∈ℜN=[x1(t),x2(t),…,xN(t)]T be the state vector of a given dynamical system, and let X∈ℜN×Nd with Nd≥N be the so-called snapshot matrix that contains a finite number of samples or snapshots of the evolution of x(t) at t=t1,t2,…,tNd. 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) φj∈ℜN,forallj=1,2,…,N:
(1)x(ti)=∑j=1N‍aj(ti)φj,∀i=1,2,…,Nd,aj(ti)=〈x(ti),φj〉,∀j=1,2,…,N,
where aj(ti) is the coordinate of x(ti) with respect to the basis vector φ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:
(2)xn(ti)=∑j=1n‍aj(ti)φj,∀i=1,2,…,Nd,n≪N.In POD, the orthonormal basis vectors are calculated in such a way that the reconstruction of the snapshots using the firstn most relevant basis vectors is optimal in the sense that the sum-squared-error (SSE) between x(ti) and xN(ti), for all i=1,…,Nd,
(3)En=∑1=1Nd‍∥x(ti)-xn(ti)∥22
is minimized. Herein ∥·∥2 denotes the L2-norm or the euclidean Norm. In other words, the POD basis vectors are the ones that solve the following constrained optimization problem:
(4)minφ1,…,φn∑i=1Nd‍∥x(ti)-∑j=1n‍〈x(ti),φj〉φj∥22
subject to
(5)φiTφj={1ifi=j,0otherwise.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 (Xsnap). The singular values of Xsnap are positive real numbers that are ordered in a decreasing way, σ1≥σ2⋯≥σn≥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 Xsnap(6)Pn=∑j=1nσj2∑j=1Nσj2,n=1,…,(N).
The ratio Pn 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 Pn=0.99. The closer Pn to 1, or similarly the closer 1-Pn 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:(7)x˙(t)=f(x(t),u(t)),y(t)=g(x(t),u(t)).Let us define a residual functionR(x˙,x) for (7) as follows:
(8)R(x˙,x)=x˙(t)-f(x(t),u(t)),
and let R(x˙n,xn) be the residual when the state vector x(t) is approximated by its nth-order approximation:
(9)xn(t)=∑j=1n‍aj(t)φj=ΦNa(t),n≪N,
where Φn=[φ1,φ2,…,φn] and a(t)=[a1(t),a2(t),…,an(t)]T. In the Galerkin projection, the projection of the residual R(x˙,x) on the space spanned by the basis vectors Φn vanishes, that is,
(10)〈R(x˙,x),φj〉=0,∀j=1,…n,
where 〈·,·〉 denotes the Euclidean inner product. Replacing x(t) by its nth-order approximation xn(t)=Φna(t) in (7):
(11)Φna˙(t)=f(Φna(t),u(t)),
and then we apply the inner product criterion (10) as follows:
(12)〈Φna˙(t),φj〉=〈f(Φna(t),u(t)),φj〉,∀j=1,2,…,nΦnTΦna˙(t)=ΦnTf(Φna(t),u(t)),
and given that ΦnTΦn=In because of the orthonormality of the basis vectors, we have the model for the POD coefficients reduce to
(13)a˙(t)=ΦnTf(Φna(t),u(t)).Finally, the reduced order model of (10) with only n states has the following form,
(14)a˙(t)=ΦnTf(Φna(t),u(t)),y(t)=g(Φna(t),u(t)).


## 2.1. General Procedure

Letx(t)∈ℜN=[x1(t),x2(t),…,xN(t)]T be the state vector of a given dynamical system, and let X∈ℜN×Nd with Nd≥N be the so-called snapshot matrix that contains a finite number of samples or snapshots of the evolution of x(t) at t=t1,t2,…,tNd. 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) φj∈ℜN,forallj=1,2,…,N:
(1)x(ti)=∑j=1N‍aj(ti)φj,∀i=1,2,…,Nd,aj(ti)=〈x(ti),φj〉,∀j=1,2,…,N,
where aj(ti) is the coordinate of x(ti) with respect to the basis vector φ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:
(2)xn(ti)=∑j=1n‍aj(ti)φj,∀i=1,2,…,Nd,n≪N.In POD, the orthonormal basis vectors are calculated in such a way that the reconstruction of the snapshots using the firstn most relevant basis vectors is optimal in the sense that the sum-squared-error (SSE) between x(ti) and xN(ti), for all i=1,…,Nd,
(3)En=∑1=1Nd‍∥x(ti)-xn(ti)∥22
is minimized. Herein ∥·∥2 denotes the L2-norm or the euclidean Norm. In other words, the POD basis vectors are the ones that solve the following constrained optimization problem:
(4)minφ1,…,φn∑i=1Nd‍∥x(ti)-∑j=1n‍〈x(ti),φj〉φj∥22
subject to
(5)φiTφj={1ifi=j,0otherwise.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 (Xsnap). The singular values of Xsnap are positive real numbers that are ordered in a decreasing way, σ1≥σ2⋯≥σn≥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 Xsnap(6)Pn=∑j=1nσj2∑j=1Nσj2,n=1,…,(N).
The ratio Pn 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 Pn=0.99. The closer Pn to 1, or similarly the closer 1-Pn 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:(7)x˙(t)=f(x(t),u(t)),y(t)=g(x(t),u(t)).Let us define a residual functionR(x˙,x) for (7) as follows:
(8)R(x˙,x)=x˙(t)-f(x(t),u(t)),
and let R(x˙n,xn) be the residual when the state vector x(t) is approximated by its nth-order approximation:
(9)xn(t)=∑j=1n‍aj(t)φj=ΦNa(t),n≪N,
where Φn=[φ1,φ2,…,φn] and a(t)=[a1(t),a2(t),…,an(t)]T. In the Galerkin projection, the projection of the residual R(x˙,x) on the space spanned by the basis vectors Φn vanishes, that is,
(10)〈R(x˙,x),φj〉=0,∀j=1,…n,
where 〈·,·〉 denotes the Euclidean inner product. Replacing x(t) by its nth-order approximation xn(t)=Φna(t) in (7):
(11)Φna˙(t)=f(Φna(t),u(t)),
and then we apply the inner product criterion (10) as follows:
(12)〈Φna˙(t),φj〉=〈f(Φna(t),u(t)),φj〉,∀j=1,2,…,nΦnTΦna˙(t)=ΦnTf(Φna(t),u(t)),
and given that ΦnTΦn=In because of the orthonormality of the basis vectors, we have the model for the POD coefficients reduce to
(13)a˙(t)=ΦnTf(Φna(t),u(t)).Finally, the reduced order model of (10) with only n states has the following form,
(14)a˙(t)=ΦnTf(Φna(t),u(t)),y(t)=g(Φna(t),u(t)).


## 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 (TJ1, TJ2 and TJ3) can be manipulated independently in order to control the concentration and temperature profiles in the reactor.Figure 2
Tubular chemical reactor with 3 cooling/heating jackets.It is assumed that radial and axial variations in concentration and temperature are present; also we consider laminar flow regime. In this study we are neglecting the heat transfer effects between the jackets fluids and the reactor wall. Under the previous assumptions, the mass balance of specieA and the energy balance on the differential cylinder shown in Figure 3 produce the following equations:
(15)∂C∂t=De∂2C∂r2+Der∂C∂r+De∂2C∂z2-Uz∂C∂z+rA,∂T∂t=KeρCp∂2T∂r2+KeρCp1r∂T∂r+KeρCp∂2T∂z2-Uz∂T∂z+ΔHrxρCprA,
where
(16)rA=-k0e-(E/RT)ρcat(CCB-CC2Keq),Keq=Keq0e[-(ΔHrx/R)(1/303-1/T)],Uz=2U0[1-(rR)2],XA=CA0-CCA0,CB=CB0-CA0XA,CC=2CA0XA,Figure 3
Cylindrical shell of thicknessΔr, length Δz, and volume 2πrΔrΔz.with the following boundary conditions:(i)
radial:
(1)
atr=0, we have symmetry ∂T/∂r=0 and ∂C/∂r=0,(2)
atr=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,
(17)-Ke∂T∂r|R=Hw(T(R,z)-Tw),(3)
atr=R, there is no mass flow through the tube walls  ∂C/∂r=0;(ii)
axial:
(1)
T=Tin and C=Cin at z=0,(2)
at the outlet of the reactorz=L, ∂T/∂z=0 and ∂C/∂z=0.Here,C(r,z,t) is the reactant concentration in [mol · m−3], T(r,z,t) is the reactant temperature in [K], De is the diffusivity of all species in [m2 · s−1], Ke is the thermal conductivity of the reaction mixture in [J · m−1 · s−1 · K−1], Keq0 is the equilibrium constant at 300K, Uz is the fluid velocity in [m · s−1], (-ΔHrx) is the heat of the reaction in [J · mol−1], ρ and Cp are the density in [kg · m−3] and the specific heat in [J · kg−1 · K−1] of the mix, respectively, k0 is the rate constant in [m6 · mol−1 · s−1 · kg−1], E is the activation energy in [J · mol−1], R is the ideal gas constant in [J · mol−1 · K−1], Hw is the heat transfer coefficient in [J · m−2 · s−1 · K−1], L is the reactor length in [m], Cin and Tin are the concentration in [mol · m−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 Tw(z,t) is the reactor wall temperature in [K] which is defined as follows:
(18)Tw={TJ1,0≤z<Za,TJ2,Za≤z<Zb,TJ3,Zb≤z<L.The parameter values of the reactor model are taken from [34]. These values are presented in Table 1.Table 1
Values of the reactor parameters with axial and radial diffusion, reaction, and convection.
Parameter
Value
D
e
1
0
-
9 m2 · s−1
v
0
1
0
-
5 m3 · s−1
L
10 m
R
0.05 m
ρ
cat
1500
kg
·
m
-
3
ρ
1000
kg
·
m
-
3
C
p
4180
J
·
k
g
-
1
·
K
-
1
Δ
H
r
x
-
83680
J
·
mo
l
-
1
K
eq
0
1000
K
e
0.559
J
·
m
-
1
·
s
-
1
·
K
-
1
k
0
1.1
×
1
0
8
m
6
·
mo
l
-
1
·
s
-
1
·
k
g
-
1
E
95238
J
·
mo
l
-
1
R
8.314
J
·
mo
l
-
1
·
K
-
1
C
in
500
mol
·
m
-
3
T
in
320 K
H
w
1300
J
·
m
-
2
·
s
-
1
·
K
-
1
H
r
0.2
s
-
1The temperature of the jacket sectionsTJ1, TJ2, and TJ3 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):
(19)De∂2C∂r2+Der∂C∂r+De∂2C∂z2-Uz∂C∂z+rA=0,KeρCp∂2T∂r2+KeρCp1r∂T∂r+KeρCp∂2T∂z2-Uz∂T∂z+ΔHrxρCprA=0withT(r,0)=Tin, C(r,0)=Cin, C0,j=C1,j, Cn+1,j=Cn,j, Ci,m+1=Ci,m, T0,j=T1,j, Tn+1,j=Tn,j, and Tm+1,j=f(Tm,j,Tw), where f(Tm,j,Tw) 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:
(20)θ1C¯i+1,j+θ2C¯i,j+1+θ3C¯i,j+θ4C¯i-1,j+θ5C¯i,j-1-1Cfr¯A=0,ϑ1T¯i+1,j+ϑ2T¯i,j+1+ϑ3T¯i,j+ϑ4T¯i-1,j+ϑ5T¯i,j-1-1Tfzϑ6r¯A=0,fori=1,2,…,n=30,j=1,2,…,m=30with(21)T¯wm,j={T¯J1=TJ1Tf,∀j=1,…,za,T¯J2=TJ2Tf,∀j=za+1,…,zb,T¯J3=TJ3Tf,∀j=zb+1,…,n,[θ1θ2θ3θ4θ5]=[DeΔr2DeΔz2-2DeΔr2(i-2)-2DeΔz2-UzΔzDeΔr2i-1iDeΔz2+UzΔz],[ϑ1ϑ2ϑ3ϑ4ϑ5ϑ6]=[KeρCpΔr2KeρCpΔz2-UzΔz-2KeρCp(1Δz2+1Δr2)+KeiρCpΔr2KeρCpΔr2-KeiρCpΔr2KeρCpΔz2+UzΔzΔHrxρCp],r¯A=-k0e-(E/RTfT¯i,j)ρcat(Ci,jCB-CC2Keq),K¯eq=Keq0e[-(ΔHrx/R)(1/303-1/TfT¯i,j)],
where n and m are the number of sections in axial and radial directions in which the reactor is divided, za and zb are the reactor sections defining the ending of the first and second jacket, respectively, Tf and Cf are normalization factors, C¯i,j=Ci,j/Cf and T¯i,j=Ti,j/Tf are the normalized concentration and temperature of the i,jth section of the reactor, T¯wm,j=Twm,j/Tf is the normalized reactor wall temperature of the i,jth section, and Δz and Δ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:
(22)minT¯J1,T¯J2,T¯J3ωm∑i=1m‍(C¯ri,n-C¯i,n)2+(1-ω)n·m∑i=1m‍∑j=1n‍(T¯ri,j-T¯i,j)2,subject to(23)steadystatemodel(20)TJminTf≤T¯J1,T¯J2,T¯J3≤TJmaxTf,T¯i,j≤TmaxTf,fori=1,2,…,m=30,forj=1,2,…,n=30,
where Cri,n is the desired concentration (normalized) at the reactor output, T¯ri,j is the desired temperature (normalized) inside the reactor of the i,jth increment, C¯i,n is the concentration (normalized) at the reactor output, ω is a trade-off coefficient, TJmin and TJmax are the lower and upper temperature values of the fluids of the jackets, and Tmax is the maximum allowed temperature inside the tubular reactor. The objective function involves an inherent trade-off between minimising conversion and energy costs. The first term of the cost function corresponds to the squared error of the normalized concentration at the reactor output (terminal cost), and the second term is related to the mean squared error of the normalized temperature along the reactor (integral cost). In this problem, and Cri,n was set to 0, T¯ri,j was selected equal to the normalized temperature of the feed flow (T¯in=Tin/Tf) for i=1,2,…,m and j=1,2,…,n. The trade-off parameter ω can take values between 0 and 1. When ω goes to 1, the reduction of the reactant concentration at the reactor output becomes more important than the temperature deviations. On the other hand when ω goes to 0, the temperature deviations become more important than the concentration at the reactor output and the risk of the formation of hot spots is reduced. To solve the optimization problem described by (22), an algorithm (a kind of sequential quadratic programming (SQP)) proposed by [35] was used in this work.The algorithm was executed in MATLAB with the following parameters:n=30, m=30, Tf=320 K, Cf=500 mol/m3, TJmin=280 K, TJmax=340 K, Tmax=335 K, Tol=10-3, and ω=0.7. The maximum allowed temperature (Tmax) inside the reactor was chosen 5 degrees below the actual limit (340 K) in order to give the feedback controller enough room of maneuverability. The trade-off coefficient ω was found by trial and error.The algorithm was executed using different initial conditions. Along the experiments, one local minima was found. The operating point was given byTJ1=291.1705K, TJ2=293.6938K, and TJ3=294.8280K. The optimal concentration and temperature profiles can be observed in Figures 4 and 5.Steady state concentration and temperature profiles.
(a)
(b)Steady state concentration and temperature profiles atr=0 m and z=10 m.
(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
(24)∂CΔ∂t=De∂2CΔ∂r2+Der∂CΔ∂r+De∂2CΔ∂z2-Uz∂CΔ∂z-αaCΔ-αbTΔ,∂TΔ∂t=KeρCp∂2TΔ∂r2+KeρCp1r∂TΔ∂r+KeρCp∂2TΔ∂z2-Uz∂TΔ∂z-ΔHrxρCpαaCΔ-ΔHrxρCpαbTΔ,
with the following boundary conditions:(i)
radial:
(1)
atr=0, we have symmetry ∂TΔ/∂r=0 and ∂CΔ/∂r=0,(2)
atr=R,
(25)-Ke∂TΔ∂r|R=Hw(TΔ(R,z)-TwΔ),(3)
atr=R; ∂CΔ/∂r=0,(ii)
axial:
(1)
TΔ=TinΔ and CΔ=CinΔ at z=0,(2)
at the outlet of the reactorz=L, ∂TΔ/∂z=0 and ∂CΔ/∂z=0.Here,CΔ=C-C*, TΔ=T-T*, and TwΔ=Tw-Tw* are the deviations from the steady state of the concentration, temperature, and reactor wall temperature; C*, T* and Tw* 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,
(26)x(t)=[C¯1,1Δ,C¯1,2Δ,…,C¯m,nΔ,T¯1,1Δ,…,T¯m,nΔ]T,d(t)=[C¯inΔ,T¯inΔ]T,u(t)=[T¯J1Δ,T¯J2Δ,T¯J3Δ]T,
then (24) can be cast as follows:
(27)x˙(t)=Ax(t)+Bu(t)+Fd(t),
where A∈ℜ2m·n×2m·n, B∈ℜ2m·n×3, and F∈ℜ2m·n×2 are the matrices describing the system, x(t)∈ℜ2m·n is the state vector, u(t)∈ℜ3 is the vector of the inputs, and d(t)∈ℜ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 inton×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 (Xsnap∈ℜ1800×10000) when independent step changes were made in the input u(t) and perturbation d(t) signals on the nonlinear model (15):
(28)Xsnap=[x(t=Δt),x(t=2Δt),…,x(t=1000Δt)].
Along the simulations, 10000 samples were collected using a sampling time Δ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 matrixXsnap:
(29)Xsnap=ΦΣΨT,
where Φ∈ℜ1800×1800 and Ψ∈ℜ10000×10000 are unitary matrices and Σ∈ℜ1800×10000 is a matrix that contains the singular values of Xsnap in a decreasing order on its main diagonal. The left singular vectors, that is, the columns of Φ,
(30)Φ=[φ1,φ2,…,φ1800],
are the POD basis vectors.


#### 3.4.3. Selection of the Most Relevant POD Basis Vectors

Then most relevant POD basis vectors are chosen using the energy criterion presented in Section 2.1. The plot of 1-PN (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-Pn=3.3×10-4(Pn=0.9996). The 50th order approximation of x(t) is given by the following truncated sequence:
(31)xn(t)=∑j=150‍aj(t)φj=ΦNa(t),
where Φn=[φ1,φ2,…,φ50] and a(t)=[a1(t),a2(t),…,a50(t)]T.Logarithmic plot ofPn and 1-Pn for determining the truncation degree of the POD basis vectors in the reactor case.
(a)
(b)


#### 3.4.4. Construction of the Model for the Firstn=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 functionR(x) for (27) as follows:
(32)R(x˙,x)=x˙(t)-Ax(t)-Bu(t)-Fd(t),
and we replace x(t) by its nth-order approximation xn(t)=Φna(t) in (32); the Galerkin projection states that the projection of R(xn) on the space spanned by the basis functions Φn vanishes, that is,
(33)〈(R(xn),φj)〉=0,j=1,…,n,
where 〈·,·〉 denotes the inner product. Replacing x(t) by its nth order approximation xn(t)=Φna(t) in (27) and applying the inner product criterion (Galerkin projection) to the resulting equation, we have
(34)〈Φna˙(t),φj〉=〈AΦna(t)+Bu(t)+Fd(t),φj〉.
By evaluating the inner product in (34),
(35)a˙(t)=ΦnTAΦna(t)+ΦnTBu(t)+ΦnTFd(t),
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
(36)a˙(t)=Ara(t)+Bru(t)+Frd(t),xn(t)=Φna(t),
where Ar=ΦnTAΦn, Br=ΦnTB, and Fr=ΦnTF. The initial condition for a(t) reads as a(0)=0.For validating the reduced order model of the reactor, we applied constant input signalsT¯J1Δ(t)=0.0312 (TJ1Δ(t)=10K), T¯J2Δ(t)=-0.0312 (TJ2Δ(t)=10 K), and T¯J3Δ(t)=0.0625 (TJ3Δ(t)=20 K) and constant perturbation signals C¯inΔ(t)=2×10-4 (CinΔ(t)=0.1 mol/m3) and T¯inΔ(t)=0.0156 (TinΔ(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 (ET) and concentration (EC) were calculated by means of the following formulas:
(37)ET=1Ns∑k=1Ns‍|T(kΔt)-Tn(kΔt)|,EC=1Ns∑k=1Ns‍|C(kΔt)-Cn(kΔt)|,
where Ns=10000 is the number of time steps and Δt=1 s. The plots of ET and EC are shown in Figure 10. For the temperature profile, the maximum value for the error ET is 0.7542 K. For the concentration, the maximum peak for the error EC is 3.3×10-3 mol/m3. 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.Temperature and concentration profiles andt=100 s at t=1000 s at r=0 m. Solid line—full-order model. Dashed line—reduced-order model.
(a)
(b)
(c)
(d)Temperature and concentration profiles att=1500 s and t=8000 s at z=10 m. Solid line—full-order model. dashed line—reduced-order model.
(a)
(b)
(c)
(d)Temperature and concentration profiles att=10000 s: full-order Model and reduced-order model.
(a)
(b)
(c)
(d)Average of the absolute error between the full-order model (15) and the reduced-order model (36).
(a)
(b)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:
(38)a(k+1)=A~a(k)+B~u(k)+F~d(k),xn(k)=Φna(k),
where A~, B~, and F~ 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:
(39)[xs(k+1)xd(k+1)]=[Iny00P][xs(k)xd(k)]+[DsDd]Δu(k)+[FsFd]Δd(k),(40)y(k)=[CsCd][xs(k)xd(k)],
where xd(k)=V1a(k),xs(k)=V2a(k-1),Δuk=uk-uk-1 is the input increment, Δdk=dk-dk-1 is the disturbance increment, and V1, V2 are transformation matrices. In the state equation defined in (39), the state component xs corresponds to the integrating poles produced by the incremental form of the model, and xd(k)=a(k) corresponds to the system modes. For stable systems, it is easy to show that when the system approaches steady state, component xd tends to zero. P is a diagonal matrix with components corresponding to the poles of the system.


## 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 (TJ1, TJ2 and TJ3) can be manipulated independently in order to control the concentration and temperature profiles in the reactor.Figure 2
Tubular chemical reactor with 3 cooling/heating jackets.It is assumed that radial and axial variations in concentration and temperature are present; also we consider laminar flow regime. In this study we are neglecting the heat transfer effects between the jackets fluids and the reactor wall. Under the previous assumptions, the mass balance of specieA and the energy balance on the differential cylinder shown in Figure 3 produce the following equations:
(15)∂C∂t=De∂2C∂r2+Der∂C∂r+De∂2C∂z2-Uz∂C∂z+rA,∂T∂t=KeρCp∂2T∂r2+KeρCp1r∂T∂r+KeρCp∂2T∂z2-Uz∂T∂z+ΔHrxρCprA,
where
(16)rA=-k0e-(E/RT)ρcat(CCB-CC2Keq),Keq=Keq0e[-(ΔHrx/R)(1/303-1/T)],Uz=2U0[1-(rR)2],XA=CA0-CCA0,CB=CB0-CA0XA,CC=2CA0XA,Figure 3
Cylindrical shell of thicknessΔr, length Δz, and volume 2πrΔrΔz.with the following boundary conditions:(i)
radial:
(1)
atr=0, we have symmetry ∂T/∂r=0 and ∂C/∂r=0,(2)
atr=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,
(17)-Ke∂T∂r|R=Hw(T(R,z)-Tw),(3)
atr=R, there is no mass flow through the tube walls  ∂C/∂r=0;(ii)
axial:
(1)
T=Tin and C=Cin at z=0,(2)
at the outlet of the reactorz=L, ∂T/∂z=0 and ∂C/∂z=0.Here,C(r,z,t) is the reactant concentration in [mol · m−3], T(r,z,t) is the reactant temperature in [K], De is the diffusivity of all species in [m2 · s−1], Ke is the thermal conductivity of the reaction mixture in [J · m−1 · s−1 · K−1], Keq0 is the equilibrium constant at 300K, Uz is the fluid velocity in [m · s−1], (-ΔHrx) is the heat of the reaction in [J · mol−1], ρ and Cp are the density in [kg · m−3] and the specific heat in [J · kg−1 · K−1] of the mix, respectively, k0 is the rate constant in [m6 · mol−1 · s−1 · kg−1], E is the activation energy in [J · mol−1], R is the ideal gas constant in [J · mol−1 · K−1], Hw is the heat transfer coefficient in [J · m−2 · s−1 · K−1], L is the reactor length in [m], Cin and Tin are the concentration in [mol · m−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 Tw(z,t) is the reactor wall temperature in [K] which is defined as follows:
(18)Tw={TJ1,0≤z<Za,TJ2,Za≤z<Zb,TJ3,Zb≤z<L.The parameter values of the reactor model are taken from [34]. These values are presented in Table 1.Table 1
Values of the reactor parameters with axial and radial diffusion, reaction, and convection.
Parameter
Value
D
e
1
0
-
9 m2 · s−1
v
0
1
0
-
5 m3 · s−1
L
10 m
R
0.05 m
ρ
cat
1500
kg
·
m
-
3
ρ
1000
kg
·
m
-
3
C
p
4180
J
·
k
g
-
1
·
K
-
1
Δ
H
r
x
-
83680
J
·
mo
l
-
1
K
eq
0
1000
K
e
0.559
J
·
m
-
1
·
s
-
1
·
K
-
1
k
0
1.1
×
1
0
8
m
6
·
mo
l
-
1
·
s
-
1
·
k
g
-
1
E
95238
J
·
mo
l
-
1
R
8.314
J
·
mo
l
-
1
·
K
-
1
C
in
500
mol
·
m
-
3
T
in
320 K
H
w
1300
J
·
m
-
2
·
s
-
1
·
K
-
1
H
r
0.2
s
-
1The temperature of the jacket sectionsTJ1, TJ2, and TJ3 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):
(19)De∂2C∂r2+Der∂C∂r+De∂2C∂z2-Uz∂C∂z+rA=0,KeρCp∂2T∂r2+KeρCp1r∂T∂r+KeρCp∂2T∂z2-Uz∂T∂z+ΔHrxρCprA=0withT(r,0)=Tin, C(r,0)=Cin, C0,j=C1,j, Cn+1,j=Cn,j, Ci,m+1=Ci,m, T0,j=T1,j, Tn+1,j=Tn,j, and Tm+1,j=f(Tm,j,Tw), where f(Tm,j,Tw) 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:
(20)θ1C¯i+1,j+θ2C¯i,j+1+θ3C¯i,j+θ4C¯i-1,j+θ5C¯i,j-1-1Cfr¯A=0,ϑ1T¯i+1,j+ϑ2T¯i,j+1+ϑ3T¯i,j+ϑ4T¯i-1,j+ϑ5T¯i,j-1-1Tfzϑ6r¯A=0,fori=1,2,…,n=30,j=1,2,…,m=30with(21)T¯wm,j={T¯J1=TJ1Tf,∀j=1,…,za,T¯J2=TJ2Tf,∀j=za+1,…,zb,T¯J3=TJ3Tf,∀j=zb+1,…,n,[θ1θ2θ3θ4θ5]=[DeΔr2DeΔz2-2DeΔr2(i-2)-2DeΔz2-UzΔzDeΔr2i-1iDeΔz2+UzΔz],[ϑ1ϑ2ϑ3ϑ4ϑ5ϑ6]=[KeρCpΔr2KeρCpΔz2-UzΔz-2KeρCp(1Δz2+1Δr2)+KeiρCpΔr2KeρCpΔr2-KeiρCpΔr2KeρCpΔz2+UzΔzΔHrxρCp],r¯A=-k0e-(E/RTfT¯i,j)ρcat(Ci,jCB-CC2Keq),K¯eq=Keq0e[-(ΔHrx/R)(1/303-1/TfT¯i,j)],
where n and m are the number of sections in axial and radial directions in which the reactor is divided, za and zb are the reactor sections defining the ending of the first and second jacket, respectively, Tf and Cf are normalization factors, C¯i,j=Ci,j/Cf and T¯i,j=Ti,j/Tf are the normalized concentration and temperature of the i,jth section of the reactor, T¯wm,j=Twm,j/Tf is the normalized reactor wall temperature of the i,jth section, and Δz and Δ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:
(22)minT¯J1,T¯J2,T¯J3ωm∑i=1m‍(C¯ri,n-C¯i,n)2+(1-ω)n·m∑i=1m‍∑j=1n‍(T¯ri,j-T¯i,j)2,subject to(23)steadystatemodel(20)TJminTf≤T¯J1,T¯J2,T¯J3≤TJmaxTf,T¯i,j≤TmaxTf,fori=1,2,…,m=30,forj=1,2,…,n=30,
where Cri,n is the desired concentration (normalized) at the reactor output, T¯ri,j is the desired temperature (normalized) inside the reactor of the i,jth increment, C¯i,n is the concentration (normalized) at the reactor output, ω is a trade-off coefficient, TJmin and TJmax are the lower and upper temperature values of the fluids of the jackets, and Tmax is the maximum allowed temperature inside the tubular reactor. The objective function involves an inherent trade-off between minimising conversion and energy costs. The first term of the cost function corresponds to the squared error of the normalized concentration at the reactor output (terminal cost), and the second term is related to the mean squared error of the normalized temperature along the reactor (integral cost). In this problem, and Cri,n was set to 0, T¯ri,j was selected equal to the normalized temperature of the feed flow (T¯in=Tin/Tf) for i=1,2,…,m and j=1,2,…,n. The trade-off parameter ω can take values between 0 and 1. When ω goes to 1, the reduction of the reactant concentration at the reactor output becomes more important than the temperature deviations. On the other hand when ω goes to 0, the temperature deviations become more important than the concentration at the reactor output and the risk of the formation of hot spots is reduced. To solve the optimization problem described by (22), an algorithm (a kind of sequential quadratic programming (SQP)) proposed by [35] was used in this work.The algorithm was executed in MATLAB with the following parameters:n=30, m=30, Tf=320 K, Cf=500 mol/m3, TJmin=280 K, TJmax=340 K, Tmax=335 K, Tol=10-3, and ω=0.7. The maximum allowed temperature (Tmax) inside the reactor was chosen 5 degrees below the actual limit (340 K) in order to give the feedback controller enough room of maneuverability. The trade-off coefficient ω was found by trial and error.The algorithm was executed using different initial conditions. Along the experiments, one local minima was found. The operating point was given byTJ1=291.1705K, TJ2=293.6938K, and TJ3=294.8280K. The optimal concentration and temperature profiles can be observed in Figures 4 and 5.Steady state concentration and temperature profiles.
(a)
(b)Steady state concentration and temperature profiles atr=0 m and z=10 m.
(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
(24)∂CΔ∂t=De∂2CΔ∂r2+Der∂CΔ∂r+De∂2CΔ∂z2-Uz∂CΔ∂z-αaCΔ-αbTΔ,∂TΔ∂t=KeρCp∂2TΔ∂r2+KeρCp1r∂TΔ∂r+KeρCp∂2TΔ∂z2-Uz∂TΔ∂z-ΔHrxρCpαaCΔ-ΔHrxρCpαbTΔ,
with the following boundary conditions:(i)
radial:
(1)
atr=0, we have symmetry ∂TΔ/∂r=0 and ∂CΔ/∂r=0,(2)
atr=R,
(25)-Ke∂TΔ∂r|R=Hw(TΔ(R,z)-TwΔ),(3)
atr=R; ∂CΔ/∂r=0,(ii)
axial:
(1)
TΔ=TinΔ and CΔ=CinΔ at z=0,(2)
at the outlet of the reactorz=L, ∂TΔ/∂z=0 and ∂CΔ/∂z=0.Here,CΔ=C-C*, TΔ=T-T*, and TwΔ=Tw-Tw* are the deviations from the steady state of the concentration, temperature, and reactor wall temperature; C*, T* and Tw* 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,
(26)x(t)=[C¯1,1Δ,C¯1,2Δ,…,C¯m,nΔ,T¯1,1Δ,…,T¯m,nΔ]T,d(t)=[C¯inΔ,T¯inΔ]T,u(t)=[T¯J1Δ,T¯J2Δ,T¯J3Δ]T,
then (24) can be cast as follows:
(27)x˙(t)=Ax(t)+Bu(t)+Fd(t),
where A∈ℜ2m·n×2m·n, B∈ℜ2m·n×3, and F∈ℜ2m·n×2 are the matrices describing the system, x(t)∈ℜ2m·n is the state vector, u(t)∈ℜ3 is the vector of the inputs, and d(t)∈ℜ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 inton×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 (Xsnap∈ℜ1800×10000) when independent step changes were made in the input u(t) and perturbation d(t) signals on the nonlinear model (15):
(28)Xsnap=[x(t=Δt),x(t=2Δt),…,x(t=1000Δt)].
Along the simulations, 10000 samples were collected using a sampling time Δ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 matrixXsnap:
(29)Xsnap=ΦΣΨT,
where Φ∈ℜ1800×1800 and Ψ∈ℜ10000×10000 are unitary matrices and Σ∈ℜ1800×10000 is a matrix that contains the singular values of Xsnap in a decreasing order on its main diagonal. The left singular vectors, that is, the columns of Φ,
(30)Φ=[φ1,φ2,…,φ1800],
are the POD basis vectors.


### 3.4.3. Selection of the Most Relevant POD Basis Vectors

Then most relevant POD basis vectors are chosen using the energy criterion presented in Section 2.1. The plot of 1-PN (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-Pn=3.3×10-4(Pn=0.9996). The 50th order approximation of x(t) is given by the following truncated sequence:
(31)xn(t)=∑j=150‍aj(t)φj=ΦNa(t),
where Φn=[φ1,φ2,…,φ50] and a(t)=[a1(t),a2(t),…,a50(t)]T.Logarithmic plot ofPn and 1-Pn for determining the truncation degree of the POD basis vectors in the reactor case.
(a)
(b)


### 3.4.4. Construction of the Model for the Firstn=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 functionR(x) for (27) as follows:
(32)R(x˙,x)=x˙(t)-Ax(t)-Bu(t)-Fd(t),
and we replace x(t) by its nth-order approximation xn(t)=Φna(t) in (32); the Galerkin projection states that the projection of R(xn) on the space spanned by the basis functions Φn vanishes, that is,
(33)〈(R(xn),φj)〉=0,j=1,…,n,
where 〈·,·〉 denotes the inner product. Replacing x(t) by its nth order approximation xn(t)=Φna(t) in (27) and applying the inner product criterion (Galerkin projection) to the resulting equation, we have
(34)〈Φna˙(t),φj〉=〈AΦna(t)+Bu(t)+Fd(t),φj〉.
By evaluating the inner product in (34),
(35)a˙(t)=ΦnTAΦna(t)+ΦnTBu(t)+ΦnTFd(t),
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
(36)a˙(t)=Ara(t)+Bru(t)+Frd(t),xn(t)=Φna(t),
where Ar=ΦnTAΦn, Br=ΦnTB, and Fr=ΦnTF. The initial condition for a(t) reads as a(0)=0.For validating the reduced order model of the reactor, we applied constant input signalsT¯J1Δ(t)=0.0312 (TJ1Δ(t)=10K), T¯J2Δ(t)=-0.0312 (TJ2Δ(t)=10 K), and T¯J3Δ(t)=0.0625 (TJ3Δ(t)=20 K) and constant perturbation signals C¯inΔ(t)=2×10-4 (CinΔ(t)=0.1 mol/m3) and T¯inΔ(t)=0.0156 (TinΔ(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 (ET) and concentration (EC) were calculated by means of the following formulas:
(37)ET=1Ns∑k=1Ns‍|T(kΔt)-Tn(kΔt)|,EC=1Ns∑k=1Ns‍|C(kΔt)-Cn(kΔt)|,
where Ns=10000 is the number of time steps and Δt=1 s. The plots of ET and EC are shown in Figure 10. For the temperature profile, the maximum value for the error ET is 0.7542 K. For the concentration, the maximum peak for the error EC is 3.3×10-3 mol/m3. 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.Temperature and concentration profiles andt=100 s at t=1000 s at r=0 m. Solid line—full-order model. Dashed line—reduced-order model.
(a)
(b)
(c)
(d)Temperature and concentration profiles att=1500 s and t=8000 s at z=10 m. Solid line—full-order model. dashed line—reduced-order model.
(a)
(b)
(c)
(d)Temperature and concentration profiles att=10000 s: full-order Model and reduced-order model.
(a)
(b)
(c)
(d)Average of the absolute error between the full-order model (15) and the reduced-order model (36).
(a)
(b)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:
(38)a(k+1)=A~a(k)+B~u(k)+F~d(k),xn(k)=Φna(k),
where A~, B~, and F~ 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:
(39)[xs(k+1)xd(k+1)]=[Iny00P][xs(k)xd(k)]+[DsDd]Δu(k)+[FsFd]Δd(k),(40)y(k)=[CsCd][xs(k)xd(k)],
where xd(k)=V1a(k),xs(k)=V2a(k-1),Δuk=uk-uk-1 is the input increment, Δdk=dk-dk-1 is the disturbance increment, and V1, V2 are transformation matrices. In the state equation defined in (39), the state component xs corresponds to the integrating poles produced by the incremental form of the model, and xd(k)=a(k) corresponds to the system modes. For stable systems, it is easy to show that when the system approaches steady state, component xd tends to zero. P is a diagonal matrix with components corresponding to the poles of the system.


## 3.4.1. Generation of the Snapshot Matrix

We have created a snapshot matrix from the system response (Xsnap∈ℜ1800×10000) when independent step changes were made in the input u(t) and perturbation d(t) signals on the nonlinear model (15):
(28)Xsnap=[x(t=Δt),x(t=2Δt),…,x(t=1000Δt)].
Along the simulations, 10000 samples were collected using a sampling time Δ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 matrixXsnap:
(29)Xsnap=ΦΣΨT,
where Φ∈ℜ1800×1800 and Ψ∈ℜ10000×10000 are unitary matrices and Σ∈ℜ1800×10000 is a matrix that contains the singular values of Xsnap in a decreasing order on its main diagonal. The left singular vectors, that is, the columns of Φ,
(30)Φ=[φ1,φ2,…,φ1800],
are the POD basis vectors.


## 3.4.3. Selection of the Most Relevant POD Basis Vectors

Then most relevant POD basis vectors are chosen using the energy criterion presented in Section 2.1. The plot of 1-PN (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-Pn=3.3×10-4(Pn=0.9996). The 50th order approximation of x(t) is given by the following truncated sequence:
(31)xn(t)=∑j=150‍aj(t)φj=ΦNa(t),
where Φn=[φ1,φ2,…,φ50] and a(t)=[a1(t),a2(t),…,a50(t)]T.Logarithmic plot ofPn and 1-Pn for determining the truncation degree of the POD basis vectors in the reactor case.
(a)
(b)


## 3.4.4. Construction of the Model for the Firstn=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 functionR(x) for (27) as follows:
(32)R(x˙,x)=x˙(t)-Ax(t)-Bu(t)-Fd(t),
and we replace x(t) by its nth-order approximation xn(t)=Φna(t) in (32); the Galerkin projection states that the projection of R(xn) on the space spanned by the basis functions Φn vanishes, that is,
(33)〈(R(xn),φj)〉=0,j=1,…,n,
where 〈·,·〉 denotes the inner product. Replacing x(t) by its nth order approximation xn(t)=Φna(t) in (27) and applying the inner product criterion (Galerkin projection) to the resulting equation, we have
(34)〈Φna˙(t),φj〉=〈AΦna(t)+Bu(t)+Fd(t),φj〉.
By evaluating the inner product in (34),
(35)a˙(t)=ΦnTAΦna(t)+ΦnTBu(t)+ΦnTFd(t),
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
(36)a˙(t)=Ara(t)+Bru(t)+Frd(t),xn(t)=Φna(t),
where Ar=ΦnTAΦn, Br=ΦnTB, and Fr=ΦnTF. The initial condition for a(t) reads as a(0)=0.For validating the reduced order model of the reactor, we applied constant input signalsT¯J1Δ(t)=0.0312 (TJ1Δ(t)=10K), T¯J2Δ(t)=-0.0312 (TJ2Δ(t)=10 K), and T¯J3Δ(t)=0.0625 (TJ3Δ(t)=20 K) and constant perturbation signals C¯inΔ(t)=2×10-4 (CinΔ(t)=0.1 mol/m3) and T¯inΔ(t)=0.0156 (TinΔ(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 (ET) and concentration (EC) were calculated by means of the following formulas:
(37)ET=1Ns∑k=1Ns‍|T(kΔt)-Tn(kΔt)|,EC=1Ns∑k=1Ns‍|C(kΔt)-Cn(kΔt)|,
where Ns=10000 is the number of time steps and Δt=1 s. The plots of ET and EC are shown in Figure 10. For the temperature profile, the maximum value for the error ET is 0.7542 K. For the concentration, the maximum peak for the error EC is 3.3×10-3 mol/m3. 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.Temperature and concentration profiles andt=100 s at t=1000 s at r=0 m. Solid line—full-order model. Dashed line—reduced-order model.
(a)
(b)
(c)
(d)Temperature and concentration profiles att=1500 s and t=8000 s at z=10 m. Solid line—full-order model. dashed line—reduced-order model.
(a)
(b)
(c)
(d)Temperature and concentration profiles att=10000 s: full-order Model and reduced-order model.
(a)
(b)
(c)
(d)Average of the absolute error between the full-order model (15) and the reduced-order model (36).
(a)
(b)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:
(38)a(k+1)=A~a(k)+B~u(k)+F~d(k),xn(k)=Φna(k),
where A~, B~, and F~ 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:
(39)[xs(k+1)xd(k+1)]=[Iny00P][xs(k)xd(k)]+[DsDd]Δu(k)+[FsFd]Δd(k),(40)y(k)=[CsCd][xs(k)xd(k)],
where xd(k)=V1a(k),xs(k)=V2a(k-1),Δuk=uk-uk-1 is the input increment, Δdk=dk-dk-1 is the disturbance increment, and V1, V2 are transformation matrices. In the state equation defined in (39), the state component xs corresponds to the integrating poles produced by the incremental form of the model, and xd(k)=a(k) corresponds to the system modes. For stable systems, it is easy to show that when the system approaches steady state, component xd tends to zero. P 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]:
(41)x(k+1)=Ax(k)+BΔu(k)+FΔd(k),y(k)=Cx(k).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:
(42)Jk,∞=∑j=1∞‍e(k+j)TQe(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j),
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 Q∈ℜny×ny and R∈ℜnu×nu 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 Q¯. For the cost defined in (42), such a terminal penalty is computed by the following Lyapunov equation [44]:
(43)Q¯-PTQ¯P=PTCdTQCdP.
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:
(44)C¯sxs(k)-ysp+C¯sD~sΔuk=0,
where
(45)D~s=[Ds⋯Ds],C¯s=diag[Cs⋯Cs],Δuk=[Δu(k)TΔu(k+1)T⋯Δu(k+p-1)T].
With the terminal penalty, the cost defined in (42) reduces to
(46)Jk,p=∑j=1p‍e(k+j)TQe(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+xd(k+m)TQ¯xd(k+m).Finally, the control optimization problem of the infinite horizon MPC can be formulated as(47)minΔukJk,p=∑j=1p‍e(k+j)TQe(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+xd(k+m)TQ¯xd(k+m)
subject to
(48)modelinincrementalform(39)C¯sxs(k)-ysp+C¯sD~sΔuk=0,-Δumax≤Δu(k+j)≤Δumax,umin≤u(k-1)+∑i=0j‍Δu(k+i)≤umax,j=0,1,…,p-1.For large changes onxs(k) or ysp or if ysp 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:(49)Jk,∞=∑j=1∞‍(e(k+j)-δk)TQ(e(k+j)-δk)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+δkTSδk,
where δk∈ℜny is a vector of slack variables and S∈ℜny×ny 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 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:
(50)C¯sxs(k)-ysp+C¯sD~sΔuk-δk=0.
Hence, the control objective defined in (49) becomes as follows:
(51)Jk,p=∑j=1p‍(e(k+j)-δk)TQ(e(k+j)-δk)+xd(k+m)TQ¯xd(k+m)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+δkTSδk.Finally, the control optimization problem of the extended infinite horizon MPC can be formulated as follows:(52)minΔuk,δkJk,p=∑j=1p‍(e(k+j)-δk)TQ(e(k+j)-δk)+xd(k+j)TQ¯xd(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+δkTSδk
subject to
(53)modelinincrementalform(39)C¯sxs(k)-ysp+C¯sD~sΔuk-δk=0,-Δumax≤Δu(k+j)≤Δumax,umin≤u(k-1)+∑i=0j‍Δu(k+i)≤umax,j=0,1,…,p-1.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 (280K≤TJ1(t),TJ2(t),TJ3(t)≤335K), and the control system should keep the temperature inside the reactor below 335 K. In [46], two POD-based IHMPC control schemes for a nonisothermal tubular reactor were presented. In the first scheme the formulation of the IHMPC controller is formulated in terms of the POD coefficients, and in the second scheme the IHMPC is in terms of physical variables. In the first case, the control of the reactor profiles is achieved indirectly by controlling the POD coefficients which have no physical meaning. This makes the tuning of the controller little intuitive and the definition of the control goals little flexible. This is not the case for the second IHMPC controller, where its formulation is in terms of the temperature of some selected points along the reactor and the concentration at the reactor output. In this work we want to explore the scheme based on the POD coefficients where the control of the temperature and concentration profiles is achieved indirectly and where the references (aref) of these POD coefficients can be calculated by
(54)aref=ΦnTxref,
where xref 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 formulationCd=I50×50V1,Cs=050×50V2. Since the state vector a(k) is unknown and the changes in the concentration of the feed flow (d1(k)=CinΔ(k)) are not measured directly, they are estimated by means of an observer (in this case a Kalman filter) with the following formulation:
(55)[a^(k+1)d^1(k+1)]=[A~F~C01][a^(k)d^1(k)]+[B~0]u(k)+[F~T0]d2(k)+[LaLd](y(k)-y^(k)),y^(k)=Cselx^n(k)=CselΦna^(k),
where a^ is the estimated vector of the POD coefficients, d^1(k) is the estimation C¯inΔ, d2(k) is the normalized temperature deviation of the feed flow T¯inΔ(k), y(k)∈ℜ10 is a vector containing ten temperature measurements (normalized deviations) along the reactor, y^(k) is the estimate of y(k), La and Ld are the submatrices of the observer gain (Kalman gain), FC and FT are the column vectors of F~=[F~C,F~T], and Csel is a selection matrix which selects the measured temperatures from the vector xn(k).The control horizonm was set to 10 samples; umin and umax 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·I50×50,R=10·I3×3, and S=1·I50×50. The Kalman gain matrix was computed from the following covariance matrices: Rw=1·I51×51, Rv=10-3·I10×10.


### 4.4. Simulation Results

In order to evaluate the performance of the control system, the following tests were carried out.Test 1. The temperature of the feed flow increased to 10 K at the 5000 s, and the concentration of the feed flow  decreased to 25×10-3 mol/m3 at the 5 s.Test 2. The temperature of the feed flow decreased to 10 K at the 5000 s, and the concentration of the feed flow increased to 25×10-3 mol/m3 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 Figures11, 12, 13, and 14, and the simulation results for Test 2 are presented in Figures 15, 16, 17, and 18.Steady state temperature and concentration profiles for Test 1 atz=10 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state temperature and concentration profiles for Test 1 atr=0 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state error (temperature and concentration) for Test 1.
(a)
(b)Control actions (jackets temperatures) of the MPC controller for Test 1.
(a)
(b)
(c)Steady state temperature and concentration profiles for Test 2 atz=10 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state temperature and concentration profiles for Test 2 atr=0 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state error (temperature and concentration) for Test 2.
(a)
(b)Control actions (jackets temperatures) of the MPC controller for Test 2.
(a)
(b)
(c)Furthermore, some quantities of interest are given in Table2. In this table, Tmax is the maximum temperature reached inside the reactor during the test. ΔCout is the percentage of the change of the mean steady state concentration at the reactor outlet with respect to its nominal value, that is,
(56)ΔCout%=CN-CN*CN*
where CN* is the mean nominal value (0.0339 mol/m3) and CN is the mean concentration at the reactor outlet in steady state after the test.Table 2
Performance parameters of the control systems.
Quantities
Test 1
Test 2
T
max
[
k
]
332.3218
323.6579
Δ
C
out
[
%
]
−1.0531
0.8732In general, the control schemes showed a good behavior for rejecting the disturbances (typical magnitudes:Cin=±25×10-3 mol/m3 and Tin=±10 K) and both presented a similar performance.


## 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]:
(41)x(k+1)=Ax(k)+BΔu(k)+FΔd(k),y(k)=Cx(k).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:
(42)Jk,∞=∑j=1∞‍e(k+j)TQe(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j),
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 Q∈ℜny×ny and R∈ℜnu×nu 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 Q¯. For the cost defined in (42), such a terminal penalty is computed by the following Lyapunov equation [44]:
(43)Q¯-PTQ¯P=PTCdTQCdP.
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:
(44)C¯sxs(k)-ysp+C¯sD~sΔuk=0,
where
(45)D~s=[Ds⋯Ds],C¯s=diag[Cs⋯Cs],Δuk=[Δu(k)TΔu(k+1)T⋯Δu(k+p-1)T].
With the terminal penalty, the cost defined in (42) reduces to
(46)Jk,p=∑j=1p‍e(k+j)TQe(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+xd(k+m)TQ¯xd(k+m).Finally, the control optimization problem of the infinite horizon MPC can be formulated as(47)minΔukJk,p=∑j=1p‍e(k+j)TQe(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+xd(k+m)TQ¯xd(k+m)
subject to
(48)modelinincrementalform(39)C¯sxs(k)-ysp+C¯sD~sΔuk=0,-Δumax≤Δu(k+j)≤Δumax,umin≤u(k-1)+∑i=0j‍Δu(k+i)≤umax,j=0,1,…,p-1.For large changes onxs(k) or ysp or if ysp 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:(49)Jk,∞=∑j=1∞‍(e(k+j)-δk)TQ(e(k+j)-δk)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+δkTSδk,
where δk∈ℜny is a vector of slack variables and S∈ℜny×ny 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 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:
(50)C¯sxs(k)-ysp+C¯sD~sΔuk-δk=0.
Hence, the control objective defined in (49) becomes as follows:
(51)Jk,p=∑j=1p‍(e(k+j)-δk)TQ(e(k+j)-δk)+xd(k+m)TQ¯xd(k+m)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+δkTSδk.Finally, the control optimization problem of the extended infinite horizon MPC can be formulated as follows:(52)minΔuk,δkJk,p=∑j=1p‍(e(k+j)-δk)TQ(e(k+j)-δk)+xd(k+j)TQ¯xd(k+j)+∑j=0p-1‍Δu(k+j)TRΔu(k+j)+δkTSδk
subject to
(53)modelinincrementalform(39)C¯sxs(k)-ysp+C¯sD~sΔuk-δk=0,-Δumax≤Δu(k+j)≤Δumax,umin≤u(k-1)+∑i=0j‍Δu(k+i)≤umax,j=0,1,…,p-1.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 (280K≤TJ1(t),TJ2(t),TJ3(t)≤335K), and the control system should keep the temperature inside the reactor below 335 K. In [46], two POD-based IHMPC control schemes for a nonisothermal tubular reactor were presented. In the first scheme the formulation of the IHMPC controller is formulated in terms of the POD coefficients, and in the second scheme the IHMPC is in terms of physical variables. In the first case, the control of the reactor profiles is achieved indirectly by controlling the POD coefficients which have no physical meaning. This makes the tuning of the controller little intuitive and the definition of the control goals little flexible. This is not the case for the second IHMPC controller, where its formulation is in terms of the temperature of some selected points along the reactor and the concentration at the reactor output. In this work we want to explore the scheme based on the POD coefficients where the control of the temperature and concentration profiles is achieved indirectly and where the references (aref) of these POD coefficients can be calculated by
(54)aref=ΦnTxref,
where xref 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 formulationCd=I50×50V1,Cs=050×50V2. Since the state vector a(k) is unknown and the changes in the concentration of the feed flow (d1(k)=CinΔ(k)) are not measured directly, they are estimated by means of an observer (in this case a Kalman filter) with the following formulation:
(55)[a^(k+1)d^1(k+1)]=[A~F~C01][a^(k)d^1(k)]+[B~0]u(k)+[F~T0]d2(k)+[LaLd](y(k)-y^(k)),y^(k)=Cselx^n(k)=CselΦna^(k),
where a^ is the estimated vector of the POD coefficients, d^1(k) is the estimation C¯inΔ, d2(k) is the normalized temperature deviation of the feed flow T¯inΔ(k), y(k)∈ℜ10 is a vector containing ten temperature measurements (normalized deviations) along the reactor, y^(k) is the estimate of y(k), La and Ld are the submatrices of the observer gain (Kalman gain), FC and FT are the column vectors of F~=[F~C,F~T], and Csel is a selection matrix which selects the measured temperatures from the vector xn(k).The control horizonm was set to 10 samples; umin and umax 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·I50×50,R=10·I3×3, and S=1·I50×50. The Kalman gain matrix was computed from the following covariance matrices: Rw=1·I51×51, Rv=10-3·I10×10.


## 4.4. Simulation Results

In order to evaluate the performance of the control system, the following tests were carried out.Test 1. The temperature of the feed flow increased to 10 K at the 5000 s, and the concentration of the feed flow  decreased to 25×10-3 mol/m3 at the 5 s.Test 2. The temperature of the feed flow decreased to 10 K at the 5000 s, and the concentration of the feed flow increased to 25×10-3 mol/m3 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 Figures11, 12, 13, and 14, and the simulation results for Test 2 are presented in Figures 15, 16, 17, and 18.Steady state temperature and concentration profiles for Test 1 atz=10 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state temperature and concentration profiles for Test 1 atr=0 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state error (temperature and concentration) for Test 1.
(a)
(b)Control actions (jackets temperatures) of the MPC controller for Test 1.
(a)
(b)
(c)Steady state temperature and concentration profiles for Test 2 atz=10 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state temperature and concentration profiles for Test 2 atr=0 m. Solid line—MPC. Dashed line—nominal profile (reference).
(a)
(b)Steady state error (temperature and concentration) for Test 2.
(a)
(b)Control actions (jackets temperatures) of the MPC controller for Test 2.
(a)
(b)
(c)Furthermore, some quantities of interest are given in Table2. In this table, Tmax is the maximum temperature reached inside the reactor during the test. ΔCout is the percentage of the change of the mean steady state concentration at the reactor outlet with respect to its nominal value, that is,
(56)ΔCout%=CN-CN*CN*
where CN* is the mean nominal value (0.0339 mol/m3) and CN is the mean concentration at the reactor outlet in steady state after the test.Table 2
Performance parameters of the control systems.
Quantities
Test 1
Test 2
T
max
[
k
]
332.3218
323.6579
Δ
C
out
[
%
]
−1.0531
0.8732In general, the control schemes showed a good behavior for rejecting the disturbances (typical magnitudes:Cin=±25×10-3 mol/m3 and Tin=±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 (Xsnap) 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 Section4, 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.



---

*Source: 763165-2013-04-15.xml*