Control of Dissipative Partial Differential Equation ... - Semantic Scholar
2013 American Control Conference (ACC) Washington, DC, USA, June 17-19, 2013
Control of dissipative partial differential equation systems using APOD based dynamic observer designs Davood Babaei Pourkargar and Antonios Armaou† spatial differential operators or to problems defined over irregular spatial domains such as corners and time dependent geometries. To overcome this limitation researchers have focused on data-driven techniques, where the data obtained from the system are used to compute the required basis functions [20], [21], [28]. The computation of basis functions that are subsequently utilized to obtain ROM is an important part of this design. A common approach to this task is Karhunen-Lo`eve expansion, also known as proper orthogonal decomposition (POD) [23] combined with the method of snapshots, then using the method of weighted residuals to construct the ROM. The combination of POD and the method of snapshots has been used extensively in model reduction, optimization and control of distributed processes [1], [8]. The original POD algorithm requires a priori availability of a sufficiently large ensemble of PDE solution data to compute empirical eigenfunctions. However, in practice, it is difficult to generate such an ensemble so that all possible dominant spatial modes are appreciably contained within the corresponding snapshots. The resulting eigenfunctions, therefore, are representative of the corresponding ensemble only. During closed-loop simulation, situations may arise when the existing eigenfunctions fail to accurately represent the dynamics of the PDE system. One possible solution is to continue augmenting the ensemble of snapshots, and then to recompute the eigenfunctions as more information regarding the process becomes available. However, this would require the solution of the eigenvalue-eigenvector problem, which is computationally expensive, and hence, unsuitable for online use. To achieve the computational load reduction and circumvent the issue of a priori availability of a sufficiently large ensemble of PDE solution data, the recursive computation of eigenfunctions, known as adaptive proper orthogonal decomposition (APOD) could be used as additional data from the process becomes available [20], [21], [28]. The designed feedback controller requires continuous availability of measurements of the process state [28] or a sufficiently large sensor network to provide good estimate of the process states [20]. In the present work, output feedback controllers are designed based on continuous point measurements available from limited number of sensors. The developed methodology for output feedback control is based on the successful integration of dynamic observers with static controllers, while the availability of spectral eigenfunctions or well constructed data ensemble (a requirement of most results in the literature) is relaxed in the present work. Rather the initially computed eigenfunctions will be recursively updated as the closed-
Abstract— This article focuses on output feedback control of distributed parameter systems with limited number of sensors employing adaptive proper orthogonal decomposition (APOD) methodology. The controller design issue is addressed by combining a robust state controller with a dynamic observer of the system states to reduce sensor requirements. The use of APOD methodology allows the development of locally accurate low-dimensional reduced order dynamic models (ROMs) for controller synthesis thus resulting in a computationally-efficient alternative to using large-dimensional models with global validity. The derived ROMs are subsequently employed for the design of dynamic observers and controllers. The proposed methods are successfully used to achieve the desired control objective of stabilizing the Kuramoto-Sivashinksy equation (KSE) at a desired state spatial profile.
I. I NTRODUCTION Many processes in chemical industries necessitate the consideration of transport phenomena often coupled with chemical reactions. Examples include distillation in petroleum processing, plasma enhanced chemical vapor deposition, etching and metallorganic vapor phase epitaxy (MOVPE) in semiconductor manufacturing and tin float bath in glass production. The control problem for these processes is nontrivial due to the spatially distributed nature of the associated system dynamics. Such distributed systems can be mathematically described by parabolic partial differential equations (PDEs). A characteristic of parabolic PDEs is that their long term and dominant behavior can be captured by finite number of spatial profiles [26]. Typical control synthesis approaches for parabolic PDEs utilize modal decomposition techniques to obtain ordinary differential equations (ODEs) which we term reduced order models (ROMs) [4], [5], [7]. One of the drawbacks of this approach is the possibly high dimensionality of the designed controllers. To address the high dimensionality issue, the concept of approximate inertial manifold (AIM) was employed to generate low-order ODE models that accurately describe parabolic PDEs [13], [27] and used for nonlinear controller design and dynamic optimization [7], [9]. In general, the eigenvalue-eigenfunction problem of the spatial operator for the PDE systems cannot be solved analytically so the above approaches for model reduction cannot be directly applied to systems which have nonlinear Financial support from the National Science Foundation, CAREER Award # CBET 06-44519 is gratefully acknowledged. † Department of Chemical Engineering, The Pennsylvania State University, University Park, PA 16802, USA, Tel: +1(814) 865-5316, Fax: +1(814) 865-7846, Email addresses: [email protected] (D. Babaei Pourkargar; Student Member, IEEE & AIChE), [email protected] (A. Armaou; Corresponding Author, Senior Member, IEEE & AIChE)
978-1-4799-0176-0/$31.00 ©2013 AACC
502
loop process evolves through different regions of the state space. Section II is for a few mathematical preliminaries used through out the paper. A finite dimensional system that captures the dominant dynamic of the PDE is obtained in Section III-A. Section III-B presents the off-line computation of empirical eigenfunctions and its online refinement is presented. Sections III-C and III-D contains dynamic observer and controller designs, respectively. In Section IV the proposed controller is illustrated on the Kuramoto-Sivashinksy equation, where it is called to stabilized the system at an open-loop unstable system steady state.
III. P ROBLEM F ORMULATION AND M ETHODOLOGIES Initially, snapshots of the state profile of the system during process evolution are collected to construct an ensemble of snapshots similar to the other data-driven model reduction and control methodologies. Let vk denote the snapshot of the system available at time tk and K denote the total number of obtained snapshots. This set of initial data can be obtained either experimentally by initially gathering open-loop process evolution data before activating the controller or from previously obtained process history data, or by performing off-line numerical simulations of the PDE system. Note that in previous methods, the ensemble of solutions must be constructed by gathering state profile snapshots for different values of input variable, u(t), and different initial conditions [2], [14] while recursive update of empirical eigenfunctions in APOD disposes of this requirement.
II. M ATHEMATICAL P RELIMINARIES We focus on the stabilization using the feedback control of spatially distributed processes described by the following state space description of highly dissipative PDEs ∂x = A (x) + b(z)u + f (x), ∂t Z ym = s(z)x dz
A. Finite dimensional approximation (1)
If we have the basis functions, the finite-dimensional approximation of the infinite-dimensional PDE system of (1)-(3) can be constructed using the method of weighted residuals so x(z,t) can be generally described as an infinite weighted sum of a complete set of basis functions φk (z). The following approximation xN (z,t) could be obtained by truncating the series expansion of x(z,t) up to order N:
Ω
subject to the following boundary conditions q(x,
dx d n0 −1 x , ..., n −1 ) = 0 on Γ dη dη 0
(2)
and initial condition
N
x(z, 0) = x0 (z).
xN (z,t) =
(3)
k=1
In the above PDE system, t is the time, z ∈ Ω ⊂ ℜ is the spatial coordinate and x(z,t) ∈ ℜn denotes the vector of state variables. u = [u1 , u2 , ..., ul ] ∈ ℜl denotes the vector of manipulated inputs and ym ∈ ℜr denotes the vector of measured outputs. Also Ω is the domain of the process and Γ is its boundary. A (x) is a nonlinear highly dissipative spatial differential operator of order n0 (n0 is an even number). f (x) is a smooth nonlinear vector function. s(z) is a known smooth vector function of z that describes the location and type of measurement sensors (e.g., point or distributed sensing). b(z) ∈ ℜn×l is a known smooth matrix function of z of the form [b1 (z), b2 (z), ..., bl (z)], where bi (z) describes how the ith control action ui (t) is distributed in the spatial domain Ω, e.g. point actuation could be defined using standard Dirac delta. In (2), q(.) is a sufficiently smooth nonlinear vector function, dx dη |Γ denotes the derivative in the direction perpendicular to the boundary and x0 (z) is a smooth vector function of z. The control objective is to drive and stabilize the process of (1)-(3) at a desired spatial profile, xd (z,t). Without loss of generality, we assumed the spatially uniform steady state xd (z,t) = 0 as the desired profile. The inner product and norm in the space of square integrable functions, L2 [Ω], is defined as Z
(ψ1 , ψ2 ) =
ψ∗1 (z)ψ2 (z)
and
∞
N→∞
∑ ak (t)φk (z) −→ x(z,t) = ∑ ak (t)φk (z)
(4)
k=1
where ak (t) are time varying coefficients known as system eigenmodes. The following N th order system of ODEs is obtained by substituting of (4) in (1), multiplying the PDE with the weighting functions, ϕv (z), and integrating over the entire spatial domain: N
Z
∑ a˙k (t)( k=1
Z
ϕv (z)φk (z)dz) = Ω
Ω
Z
+
Ω
k=1 N
Z
ϕv (z)b(z)udz + Ω
N
ϕv (z)A ( ∑ ak (t)φk (z))dz
ϕv (z) f ( ∑ ak (t)φk (z))dz, k=1
v = 1, ..., N. (5) The type of weighted residual method could be determined by the weighting functions in the above equation. The method of weighted residual reduces to Galerkin’s method when the weighting functions, ϕv (z), and the basis functions, φv (z), are the same. Then (5) can be summarized as a˙ = F (a) + G (a)u
(6)
where F and G are the smooth vector functions of the modes. Such an approach will lead to accurate finite dimensional ODE models for processes that exhibit strong diffusive phenomena (and can thus be described by highly dissipative PDE systems); this is due to the fact that their dominant dynamic behavior can be captured by a finite number of dominant spatial profiles [26].
kψ1 k2 = (ψ1 , ψ1 )1/2
Ω
where ψ∗ denotes the complex conjugate transpose. 503
Based on the new values of Z, we then compute the revised eigenfunctions φ1 , φ2 , . . . , φm as a linear combination of the snapshots given by the following equation
B. Off-line computation of empirical eigenfunctions and online recursive refinement Available off-line solution data of the process are used to construct the basis functions that are necessary for derivation of the reduced order models of (5). POD is used to compute an orthogonal set of empirical eigenfunctions that represent the ensemble [23]. We specifically employ the method of snapshots variant of POD which through the solution of an eigenvector problem of the covariance matrix of the snapshot ensemble identifies the eigenfunctions; this may also be computationally demanding depending on the ensemble size. The key features of the algorithm and a detailed analysis are presented in [23]. We bypass these computational demands by updating the empirical eigenfunctions once new measurements from the process become available as discussed in [20], [21], [28]. For completeness we briefly present the algorithm in this section, called APOD. The recursive algorithm is based on the partition of the eigenspace of the covariance matrix into two subspaces; the dominant one containing the modes which capture at least ε percent of energy (in an information theory sense) in the ensemble (denoted as P) and the orthogonal complement to P containing the rest of the modes (denoted as Q). Such a partition is possible since the dominant dynamics of highly dissipative PDEs are finite dimensional [26]. Based on [20], [21], [28], the orthonormal basis for the subspace P is recursively updated upon the arrival of new snapshots, while the orthonormal basis for Q is computed as the orthonormal complement of P. The algorithm is computationally efficient as long as the dimension of P is small (which partly depends on choosing an appropriate value for ε). Let CK denote the covariance matrix obtained from K data snapshots. It was assumed that m eigenvectors out of K possible eigenvectors of CK , have the corresponding m K ε ; i.e. m eigenmodes eigenvalues such that ∑ λi / ∑ λi ≥ 100 i=1 i=1 of CK capture ε percent of energy in the ensemble. An orthonormal basis for the subspace P can be obtained as: Z = [ψ1 , ψ2 , . . . , ψm ], Z ∈ RK×m
φ(z) = ∑ψk vk (z)
(8)
k
where ψk denotes the kth element of vector ψ. Such a procedure guarantees that the eigenfunctions continue to locally describe the dominant behavior of the PDE.
Fig. 1.
Flow chart of APOD recursive algorithm [20], [28].
A flow chart illustrating the above steps is presented in Figure 1. Note that in the flow chart ξ corresponds to the contribution of the dominant eigenvalue of cq = QCK Q, namely λm+1 towards the total energy of the ensemble, i.e
(7)
where ψ1 , ψ2 , . . . , ψm denote the eigenvectors of CK that correspond to the eigenvalues λ1 , λ2 , . . . , λm . Note that the eigenfunctions computed by these eigenvectors also capture the dominant dynamics of the PDE system of (1-3) locally. As a new snapshot from the process becomes available, the subspace P may change in the following three ways. • The dimension of the dominant subspace P increases i.e., some modes within Q become necessary to capture the desired percentage of energy in the ensemble. • Some of the eigenmodes of subspace P are no longer necessary to capture the required ε percent of the energy. The basis Z is updated and the dimension simultaneously decreased. • The dimensionality of P remains unchanged, however the basis Z needs to be updated in order to accurately capture the newly added snapshot.
ξ=
λm+1 ∑m+1 i=1 λi
and H = Z T CK Z is an m × m matrix (see [20], [28] for additional details). In the proposed approach for model reduction the snapshots used are obtained during closed-loop system evolution as opposed to all other proper orthogonal decompositionbased reduction approaches that are based on open-loop snapshots, and thus, these snapshots and the resulting ROM account for the impact of controller functional form on the process. It is important to note this intimate relation between APOD and the controller. C. Dynamic observer We assume that the snapshots of the process become available only periodically and the point measurements from 504
a restricted number of sensors are continuously available which is quite common in industrial processes. Static output observers, based on the continuous point measurements, have been designed to estimate the modes of (6) that are required for controller design [20], [28]. But the number of sensors should be at least equal to the number of modes to estimate the modes of (6) using the linear static observer, otherwise the static observer could not predict the values of the modes. To overcome these issues, we will design a dynamic observer that conceptually needs only one point measurement to predict the dynamic behavior of (6). The theory of the dynamic observer design has been developed by Luenberger [19] and offers a complete and comprehensive analysis but the nonlinear dynamic observer design is much more complicated and has received a considerable attention in recent years [16], [17], [25]. The resulting dynamic observer based on (6) will have the following structure a˙˜ = F (a) ˜ + G (a)u ˜ + L(ym − S (a)) ˜
The time derivative of the Lyapunov function using (12) will be (15) V˙ = eT (PAc + ATc P)e. Since PAc + ATc P < 0, it follows that V˙ < 0 for all e 6= 0 and stability is proved. So P and L will be chosen so that the following LMI is satisfied ATc P + PAc < −PW P − PLULT P
where W and U are the design weighting matrices; they should be symmetric positive definite. From (16) we have ATc P + PAc + PW P + PLULT P < 0 ⇒ (AT −CT LT )P + P(A − LC) + PW P + PLULT P < 0 ⇒ AT P + PA − (PLC)T − PLC + PW P + PLULT P < 0. (17) Using Y = PL, (17) can be written as AT P + PA − (YC)T −YC + PW P +YUY T < 0 ⇒ PA + AT P + (Y T −U −1C)T U(Y T −U −1C) −CT U −1C + PW P < 0.
(9)
(18)
where a˜ is the vector of estimated modes, L is the observer gain and S is a smooth vector function. We present the following two approaches to compute the observer gain, L, that could stabilize the system of (9) in the Lyapunov sense. 1) Pole placement: To implement the pole placement technique, we consider the linearized system of (6) a˙ = Aa + Bu, y = Ca
By setting Y T = U −1C we have PA + AT P −CT U −1C + PW P < 0. The above matrix inequality is equivalent to PA + AT P −CT U −1C P
Control of dissipative partial differential equation systems using APOD based dynamic observer designs Davood Babaei Pourkargar and Antonios Armaou† spatial differential operators or to problems defined over irregular spatial domains such as corners and time dependent geometries. To overcome this limitation researchers have focused on data-driven techniques, where the data obtained from the system are used to compute the required basis functions [20], [21], [28]. The computation of basis functions that are subsequently utilized to obtain ROM is an important part of this design. A common approach to this task is Karhunen-Lo`eve expansion, also known as proper orthogonal decomposition (POD) [23] combined with the method of snapshots, then using the method of weighted residuals to construct the ROM. The combination of POD and the method of snapshots has been used extensively in model reduction, optimization and control of distributed processes [1], [8]. The original POD algorithm requires a priori availability of a sufficiently large ensemble of PDE solution data to compute empirical eigenfunctions. However, in practice, it is difficult to generate such an ensemble so that all possible dominant spatial modes are appreciably contained within the corresponding snapshots. The resulting eigenfunctions, therefore, are representative of the corresponding ensemble only. During closed-loop simulation, situations may arise when the existing eigenfunctions fail to accurately represent the dynamics of the PDE system. One possible solution is to continue augmenting the ensemble of snapshots, and then to recompute the eigenfunctions as more information regarding the process becomes available. However, this would require the solution of the eigenvalue-eigenvector problem, which is computationally expensive, and hence, unsuitable for online use. To achieve the computational load reduction and circumvent the issue of a priori availability of a sufficiently large ensemble of PDE solution data, the recursive computation of eigenfunctions, known as adaptive proper orthogonal decomposition (APOD) could be used as additional data from the process becomes available [20], [21], [28]. The designed feedback controller requires continuous availability of measurements of the process state [28] or a sufficiently large sensor network to provide good estimate of the process states [20]. In the present work, output feedback controllers are designed based on continuous point measurements available from limited number of sensors. The developed methodology for output feedback control is based on the successful integration of dynamic observers with static controllers, while the availability of spectral eigenfunctions or well constructed data ensemble (a requirement of most results in the literature) is relaxed in the present work. Rather the initially computed eigenfunctions will be recursively updated as the closed-
Abstract— This article focuses on output feedback control of distributed parameter systems with limited number of sensors employing adaptive proper orthogonal decomposition (APOD) methodology. The controller design issue is addressed by combining a robust state controller with a dynamic observer of the system states to reduce sensor requirements. The use of APOD methodology allows the development of locally accurate low-dimensional reduced order dynamic models (ROMs) for controller synthesis thus resulting in a computationally-efficient alternative to using large-dimensional models with global validity. The derived ROMs are subsequently employed for the design of dynamic observers and controllers. The proposed methods are successfully used to achieve the desired control objective of stabilizing the Kuramoto-Sivashinksy equation (KSE) at a desired state spatial profile.
I. I NTRODUCTION Many processes in chemical industries necessitate the consideration of transport phenomena often coupled with chemical reactions. Examples include distillation in petroleum processing, plasma enhanced chemical vapor deposition, etching and metallorganic vapor phase epitaxy (MOVPE) in semiconductor manufacturing and tin float bath in glass production. The control problem for these processes is nontrivial due to the spatially distributed nature of the associated system dynamics. Such distributed systems can be mathematically described by parabolic partial differential equations (PDEs). A characteristic of parabolic PDEs is that their long term and dominant behavior can be captured by finite number of spatial profiles [26]. Typical control synthesis approaches for parabolic PDEs utilize modal decomposition techniques to obtain ordinary differential equations (ODEs) which we term reduced order models (ROMs) [4], [5], [7]. One of the drawbacks of this approach is the possibly high dimensionality of the designed controllers. To address the high dimensionality issue, the concept of approximate inertial manifold (AIM) was employed to generate low-order ODE models that accurately describe parabolic PDEs [13], [27] and used for nonlinear controller design and dynamic optimization [7], [9]. In general, the eigenvalue-eigenfunction problem of the spatial operator for the PDE systems cannot be solved analytically so the above approaches for model reduction cannot be directly applied to systems which have nonlinear Financial support from the National Science Foundation, CAREER Award # CBET 06-44519 is gratefully acknowledged. † Department of Chemical Engineering, The Pennsylvania State University, University Park, PA 16802, USA, Tel: +1(814) 865-5316, Fax: +1(814) 865-7846, Email addresses: [email protected] (D. Babaei Pourkargar; Student Member, IEEE & AIChE), [email protected] (A. Armaou; Corresponding Author, Senior Member, IEEE & AIChE)
978-1-4799-0176-0/$31.00 ©2013 AACC
502
loop process evolves through different regions of the state space. Section II is for a few mathematical preliminaries used through out the paper. A finite dimensional system that captures the dominant dynamic of the PDE is obtained in Section III-A. Section III-B presents the off-line computation of empirical eigenfunctions and its online refinement is presented. Sections III-C and III-D contains dynamic observer and controller designs, respectively. In Section IV the proposed controller is illustrated on the Kuramoto-Sivashinksy equation, where it is called to stabilized the system at an open-loop unstable system steady state.
III. P ROBLEM F ORMULATION AND M ETHODOLOGIES Initially, snapshots of the state profile of the system during process evolution are collected to construct an ensemble of snapshots similar to the other data-driven model reduction and control methodologies. Let vk denote the snapshot of the system available at time tk and K denote the total number of obtained snapshots. This set of initial data can be obtained either experimentally by initially gathering open-loop process evolution data before activating the controller or from previously obtained process history data, or by performing off-line numerical simulations of the PDE system. Note that in previous methods, the ensemble of solutions must be constructed by gathering state profile snapshots for different values of input variable, u(t), and different initial conditions [2], [14] while recursive update of empirical eigenfunctions in APOD disposes of this requirement.
II. M ATHEMATICAL P RELIMINARIES We focus on the stabilization using the feedback control of spatially distributed processes described by the following state space description of highly dissipative PDEs ∂x = A (x) + b(z)u + f (x), ∂t Z ym = s(z)x dz
A. Finite dimensional approximation (1)
If we have the basis functions, the finite-dimensional approximation of the infinite-dimensional PDE system of (1)-(3) can be constructed using the method of weighted residuals so x(z,t) can be generally described as an infinite weighted sum of a complete set of basis functions φk (z). The following approximation xN (z,t) could be obtained by truncating the series expansion of x(z,t) up to order N:
Ω
subject to the following boundary conditions q(x,
dx d n0 −1 x , ..., n −1 ) = 0 on Γ dη dη 0
(2)
and initial condition
N
x(z, 0) = x0 (z).
xN (z,t) =
(3)
k=1
In the above PDE system, t is the time, z ∈ Ω ⊂ ℜ is the spatial coordinate and x(z,t) ∈ ℜn denotes the vector of state variables. u = [u1 , u2 , ..., ul ] ∈ ℜl denotes the vector of manipulated inputs and ym ∈ ℜr denotes the vector of measured outputs. Also Ω is the domain of the process and Γ is its boundary. A (x) is a nonlinear highly dissipative spatial differential operator of order n0 (n0 is an even number). f (x) is a smooth nonlinear vector function. s(z) is a known smooth vector function of z that describes the location and type of measurement sensors (e.g., point or distributed sensing). b(z) ∈ ℜn×l is a known smooth matrix function of z of the form [b1 (z), b2 (z), ..., bl (z)], where bi (z) describes how the ith control action ui (t) is distributed in the spatial domain Ω, e.g. point actuation could be defined using standard Dirac delta. In (2), q(.) is a sufficiently smooth nonlinear vector function, dx dη |Γ denotes the derivative in the direction perpendicular to the boundary and x0 (z) is a smooth vector function of z. The control objective is to drive and stabilize the process of (1)-(3) at a desired spatial profile, xd (z,t). Without loss of generality, we assumed the spatially uniform steady state xd (z,t) = 0 as the desired profile. The inner product and norm in the space of square integrable functions, L2 [Ω], is defined as Z
(ψ1 , ψ2 ) =
ψ∗1 (z)ψ2 (z)
and
∞
N→∞
∑ ak (t)φk (z) −→ x(z,t) = ∑ ak (t)φk (z)
(4)
k=1
where ak (t) are time varying coefficients known as system eigenmodes. The following N th order system of ODEs is obtained by substituting of (4) in (1), multiplying the PDE with the weighting functions, ϕv (z), and integrating over the entire spatial domain: N
Z
∑ a˙k (t)( k=1
Z
ϕv (z)φk (z)dz) = Ω
Ω
Z
+
Ω
k=1 N
Z
ϕv (z)b(z)udz + Ω
N
ϕv (z)A ( ∑ ak (t)φk (z))dz
ϕv (z) f ( ∑ ak (t)φk (z))dz, k=1
v = 1, ..., N. (5) The type of weighted residual method could be determined by the weighting functions in the above equation. The method of weighted residual reduces to Galerkin’s method when the weighting functions, ϕv (z), and the basis functions, φv (z), are the same. Then (5) can be summarized as a˙ = F (a) + G (a)u
(6)
where F and G are the smooth vector functions of the modes. Such an approach will lead to accurate finite dimensional ODE models for processes that exhibit strong diffusive phenomena (and can thus be described by highly dissipative PDE systems); this is due to the fact that their dominant dynamic behavior can be captured by a finite number of dominant spatial profiles [26].
kψ1 k2 = (ψ1 , ψ1 )1/2
Ω
where ψ∗ denotes the complex conjugate transpose. 503
Based on the new values of Z, we then compute the revised eigenfunctions φ1 , φ2 , . . . , φm as a linear combination of the snapshots given by the following equation
B. Off-line computation of empirical eigenfunctions and online recursive refinement Available off-line solution data of the process are used to construct the basis functions that are necessary for derivation of the reduced order models of (5). POD is used to compute an orthogonal set of empirical eigenfunctions that represent the ensemble [23]. We specifically employ the method of snapshots variant of POD which through the solution of an eigenvector problem of the covariance matrix of the snapshot ensemble identifies the eigenfunctions; this may also be computationally demanding depending on the ensemble size. The key features of the algorithm and a detailed analysis are presented in [23]. We bypass these computational demands by updating the empirical eigenfunctions once new measurements from the process become available as discussed in [20], [21], [28]. For completeness we briefly present the algorithm in this section, called APOD. The recursive algorithm is based on the partition of the eigenspace of the covariance matrix into two subspaces; the dominant one containing the modes which capture at least ε percent of energy (in an information theory sense) in the ensemble (denoted as P) and the orthogonal complement to P containing the rest of the modes (denoted as Q). Such a partition is possible since the dominant dynamics of highly dissipative PDEs are finite dimensional [26]. Based on [20], [21], [28], the orthonormal basis for the subspace P is recursively updated upon the arrival of new snapshots, while the orthonormal basis for Q is computed as the orthonormal complement of P. The algorithm is computationally efficient as long as the dimension of P is small (which partly depends on choosing an appropriate value for ε). Let CK denote the covariance matrix obtained from K data snapshots. It was assumed that m eigenvectors out of K possible eigenvectors of CK , have the corresponding m K ε ; i.e. m eigenmodes eigenvalues such that ∑ λi / ∑ λi ≥ 100 i=1 i=1 of CK capture ε percent of energy in the ensemble. An orthonormal basis for the subspace P can be obtained as: Z = [ψ1 , ψ2 , . . . , ψm ], Z ∈ RK×m
φ(z) = ∑ψk vk (z)
(8)
k
where ψk denotes the kth element of vector ψ. Such a procedure guarantees that the eigenfunctions continue to locally describe the dominant behavior of the PDE.
Fig. 1.
Flow chart of APOD recursive algorithm [20], [28].
A flow chart illustrating the above steps is presented in Figure 1. Note that in the flow chart ξ corresponds to the contribution of the dominant eigenvalue of cq = QCK Q, namely λm+1 towards the total energy of the ensemble, i.e
(7)
where ψ1 , ψ2 , . . . , ψm denote the eigenvectors of CK that correspond to the eigenvalues λ1 , λ2 , . . . , λm . Note that the eigenfunctions computed by these eigenvectors also capture the dominant dynamics of the PDE system of (1-3) locally. As a new snapshot from the process becomes available, the subspace P may change in the following three ways. • The dimension of the dominant subspace P increases i.e., some modes within Q become necessary to capture the desired percentage of energy in the ensemble. • Some of the eigenmodes of subspace P are no longer necessary to capture the required ε percent of the energy. The basis Z is updated and the dimension simultaneously decreased. • The dimensionality of P remains unchanged, however the basis Z needs to be updated in order to accurately capture the newly added snapshot.
ξ=
λm+1 ∑m+1 i=1 λi
and H = Z T CK Z is an m × m matrix (see [20], [28] for additional details). In the proposed approach for model reduction the snapshots used are obtained during closed-loop system evolution as opposed to all other proper orthogonal decompositionbased reduction approaches that are based on open-loop snapshots, and thus, these snapshots and the resulting ROM account for the impact of controller functional form on the process. It is important to note this intimate relation between APOD and the controller. C. Dynamic observer We assume that the snapshots of the process become available only periodically and the point measurements from 504
a restricted number of sensors are continuously available which is quite common in industrial processes. Static output observers, based on the continuous point measurements, have been designed to estimate the modes of (6) that are required for controller design [20], [28]. But the number of sensors should be at least equal to the number of modes to estimate the modes of (6) using the linear static observer, otherwise the static observer could not predict the values of the modes. To overcome these issues, we will design a dynamic observer that conceptually needs only one point measurement to predict the dynamic behavior of (6). The theory of the dynamic observer design has been developed by Luenberger [19] and offers a complete and comprehensive analysis but the nonlinear dynamic observer design is much more complicated and has received a considerable attention in recent years [16], [17], [25]. The resulting dynamic observer based on (6) will have the following structure a˙˜ = F (a) ˜ + G (a)u ˜ + L(ym − S (a)) ˜
The time derivative of the Lyapunov function using (12) will be (15) V˙ = eT (PAc + ATc P)e. Since PAc + ATc P < 0, it follows that V˙ < 0 for all e 6= 0 and stability is proved. So P and L will be chosen so that the following LMI is satisfied ATc P + PAc < −PW P − PLULT P
where W and U are the design weighting matrices; they should be symmetric positive definite. From (16) we have ATc P + PAc + PW P + PLULT P < 0 ⇒ (AT −CT LT )P + P(A − LC) + PW P + PLULT P < 0 ⇒ AT P + PA − (PLC)T − PLC + PW P + PLULT P < 0. (17) Using Y = PL, (17) can be written as AT P + PA − (YC)T −YC + PW P +YUY T < 0 ⇒ PA + AT P + (Y T −U −1C)T U(Y T −U −1C) −CT U −1C + PW P < 0.
(9)
(18)
where a˜ is the vector of estimated modes, L is the observer gain and S is a smooth vector function. We present the following two approaches to compute the observer gain, L, that could stabilize the system of (9) in the Lyapunov sense. 1) Pole placement: To implement the pole placement technique, we consider the linearized system of (6) a˙ = Aa + Bu, y = Ca
By setting Y T = U −1C we have PA + AT P −CT U −1C + PW P < 0. The above matrix inequality is equivalent to PA + AT P −CT U −1C P
