State Observer Design for Non Linear Coupled ... - Semantic Scholar
53rd IEEE Conference on Decision and Control December 15-17, 2014. Los Angeles, California, USA
State observer design for non linear coupled partial differential equations with application to radiative-conductive heat transfer systems Mohamed Ghattassi1 , Mohamed Boutayeb2 and Jean R. Roche3 Abstract— This contribution deals with state observer design for a class of non linear coupled PDE that describe radiativeconductive heat transfer systems. This approach uses first a stable spatial discretization technique that is the Galerkin method to obtain a large scale but finite dimensional system in a suitable form. Thanks to the special structure of the obtained state system, the second main result is to show through the differential mean value theorem (DMVT) that there always exists an observer gain matrix that assures asymptotic convergence. On the other hand, in order to avoid high computational requirements, we show how to construct the observer gain matrix so that the stability condition, written in terms of linear matrix inequality, is satisfied. Extension to H∞ performance analysis is also proposed. In order to show high accuracy of the proposed technique, a numerical example is provided.
I. INTRODUCTION Control design and state observers for non linear partial differential equations (PDE) is known to be a hard task and is still being a challenging problem in the literature. Motivated by this fact, a typical approach to the state observers of PDE systems is to obtain an approximate ordinary differential equation (ODE) representation of the original PDE system by utilizing the spatial discretization techniques [4], [2], [7] which is then used for observer design purposes by applying different existing ODE-based linear or nonlinear observer techniques. This approach is also referred to the reduce-thendesign (RTD) approach. In this work, we consider the nonlinear coupled PDE systems described by the radiative-conduction heat transfer phenomena. This nonlinear problem has not been tackled in the literature nor in infinite dimension or finite dimensional. The model reduction is based on Galerkin (DG) method . The radiative transfer equation (RTE) was solved using the Discontinuous Galerkin method with an upwind numerical flux. The nonlinear heat equation was discretized using a continuous Galerkin (finite element) method. This numerical method is stable and gives a suitable form for the synthesis of the observer. Various research activities have focused on observer design for nonlinear systems. This topic was motivated by the fact that the state estimation can be used for control and diagnosis. Linear parameter-varying (LPV) techniques
have been developed into effective tools to control multiinput multi-output (MIMO) nonlinear plants. Many nonlinear systems of practical interest can be represented as LPV systems. However, the application of these techniques has not yet been extensively explored in control of PDEs. A LPV approach together with a nonlinear observer was proposed in Hashemi and Werner [8] to control the one-dimensional Burgers equation. The basic idea of this work is to use the well-known DMVT which allows to write the dynamics of the estimation error as a class of LPV systems. Stability of the estimation error is analyzed using the convexity principle and the Lyapunov stability theory [15], [16], [6]. The observer gain guaranteeing the global convergence of the proposed observer is computed by LMI. The structure of the state system gives the advantage of DMVT to reduce the computational complexity and to construct the gain observer which ensures the stability of the proposed observer. This result was extended to the observer synthesis with H∞ performance analysis. This article is organized as follows. Next section develops the governing equations for two-dimensional combined radiative and conductive heat transfer. In section 3, the reduced model of the coupled partial differential equations is given. DMVT for observer design of the reduction model of the coupled partial differential equations are described in section 4. The section 5, an extension to H∞ performance analysis is proposed for the results of the previous section. Simulations results are given in section 6 and the last section draws conclusions. Notations.Throughout this paper, we will use the following notation: • • • • •
•
1 M. Ghattassi is with Institut Elie Cartan de Lorraine, Universit´ e de Lorraine, CNRS, INRIA, B.P. 70239, 54506, Vandoeuvre l`es Nancy, France.
k.k is the Euclidean norm. (?) is used for the blocks induced by symmetry. I represents the identity matrix of appropriate dimension. AT represents the transposed matrix of A. For a square matrix S , S > 0 (S < 0) means that this matrix is positive definite (negative definite). ith z}|{ T es (i) = (0, ..., 0, 1 , 0, ..., 0) ∈ Rs ,s > 1 is a vector | {z } s components of the canonical basis of Rs .
[email protected] 2 M. Boutayeb is with the Centre des Recherches en Automatique, Universit´e de Lorraine, CNRS, B.P. 70239, 54506, Vandoeuvre l`es Nancy, France. [email protected] 3 J.R. Roche is with the Institut Elie Cartan de Lorraine, Universit´ e de Lorraine, CNRS, INRIA, B.P. 70239, 54506, Vandoeuvre l`es Nancy, France.
[email protected] 978-1-4673-6090-6/14/$31.00 ©2014 European Union
II. M ODEL PROBLEM Let Ω be a bounded domain in R2 and D the unit disk. Thus, for the problem considered, β ∈ D = {β ∈ R2 : |β| 6 1}, s? = (s? , s2 ) ∈ Ω and ξ ∈ [0, τ ], for τ > 0. Let n be
1569
the outward unit normal to the boundary ∂Ω. We denote ∂Ω− = {(s, β) ∈ ∂Ω × D such that β.n < 0}. The unknown of the RTE is the dimensionless radiation intensity denoted I(ξ, s, β) given at time ξ, position s and in the direction β. The unknown of the energy equation is the dimensionless temperature T (ξ, s) at time ξ and position s. The RTE is given by (see, [1]) β.∇I(ξ, s, β) + I(ξ, s, β) = T 4 (ξ, s),
Now, we proceed with the spatial discretization of the SN transport equation using the DG method . We consider a triangulation Th of Ω such that [ Ω= K, K∈Th • •
(1)
•
for all (ξ, s, β) ∈ [0, τ ] × Ω × D with the boundary condition I(ξ, s, β) = g(ξ, s, β),
•
(2)
for all (ξ, s, β) ∈ [0, τ ] × ∂Ω− , the radiative boundary condition is given by
Moreover, for any mesh element K ∈ Th , the set
g(t, s, β) = T 4 (ξ, s) for (ξ, s, β) ∈ [0, τ ] × ∂Ω− ,
FK = {F ∈ Fh |F ∈ ∂K},
where T (ξ, s) is the temperature at the boundary of the medium. The dimensionless radiative equation (1) is combined with the following energy equation (nonlinear heat equation) : ∂T − 4T + θT 4 = θG ∂ξ T (0, s) = T0 (s)
for (ξ, s) ∈ [0, τ ] × Ω
(3a)
for all s ∈ Ω
(3b)
where θ is the dimensionless constant and G is the incident radiation intensity, Z G(ξ, s) = I(ξ, s, β) dβ for (ξ, s) ∈ [0, τ ] × Ω.
•
Vh = {Ih ∈ L2 (Ω)/ ∀K ∈ Th , Ih |
K
∈ Pk (K)},
(ξ, s) ∈ [0, τ ] × ∂Ω,
T (ξ, s) = f (ξ, s)
(4)
III. A NALYSIS OF RADIATIVE CONDUCTIVE HEAT TRANSFER
A. Radiative Transfer Equation To compute a numerical solution of the RTE, the classical SN −discrete ordinate method is introduced [12]. The SN method consists of replacing the radiation intensity I(ξ, s, β) by a discrete radiation intensity (I 1 (ξ, s), I 2 (ξ, s), β ..., I N (ξ, s)), where I l (ξ, s) = I(ξ, s, βl ) for all l ∈ {1, ..., N β}. A quadrature rule was chosen {(βl , wl ) l = 1, ..., N β }, such that: β
Z f (β, s)dβ '
N X
β
f (βl , s)wl and
l=1
N X
wl = 4π.
l=1
This quadrature formula preserves a symmetry for the opposite βl0 = −βl of discrete ordinate βl in the quadrature set, the weights are equal wl0 = wl , see [12]. The discrete radiation intensity I l is the solution of a N β system of partial differential equations over Ω: l l 4 βl .∇I (ξ, s) + I (ξ, s) = T (ξ, s)
I l (ξ, s) = g(ξ, s, βl )
N
(6)
and denote a generic element in Vh as Ih (ξ, .) = l Nβ Ih (ξ, .) l=1 . Due to the discontinuous nature of the spatial approximation, functions Ihl ∈ Vh are double-valued on interior faces. Consider an interior face F ∈ Fhin separating two mesh cells, K1 and K2 . The mean value and jump of a function Ihl (ξ, .) ∈ Vh are defined as follows: 1 l (I (ξ, .) + I2l (ξ, .)), 2 1 [[Ihl (ξ, .)]] = (I1l (ξ, .) − I2l (ξ, .)),
{{Ihl (ξ, .)}} =
for all ξ ∈ [0, τ ], where I1l (ξ, .) = Ihl (ξ, .)|K and I2l (ξ, .) = 1 Ihl (ξ, .)|K are the restrictions of Ihl (ξ, .) on the mesh cells 2 K1 and K2 , respectively. The DG formulation is obtained by multiplying the SN equation for the direction βl with the test function wh ∈ Vh and applying the upwind numerical ˆ fluxes F(s) to approximate the quantity (βl .n)I l (ξ, .) on the elements boundary K: Z K
κIhl (ξ, .)wh − (βl .∇h wh )Ihl (ξ, .) +
Z ∂K
Fˆ (s)wh =
Z z(ξ, .)wh , Ω
(7)
ˆ ∀ξ ∈ [0, τ ]. The upwind numerical flux F(s) at the mesh interface F from K1 to K2 is given by:
(ξ, s) ∈ [0, τ ] × Ω
1 Fˆ (s) = βl .nF {{I l (ξ, s)}} + |βl .nF |[[I l (ξ, s)]], ∀ξ ∈ [0, τ ], 2
(ξ, s) ∈ [0, τ ] × ∂Ωl− ,
∂Ωl− = {s ∈ ∂Ω such that βl .n < 0}.
(5)
where Pk (K) denotes the set of polynomial defined in K of degree less than or equal to k. Define: β Vh = {Vh }l=1 ,
We consider the Dirichlet boundary conditions, where the dimensionless temperature at the boundary is known
where
collects the mesh faces composing the boundary of K. diam(K) = h for each K ∈ Th (structured mesh).
Now, an approximation space Vh ⊂ L2 (Ω) is introduced such that
D
D
Each K is a triangle with nonempty internal; ◦ T ◦ K1 K2 T = ∅ for each distinct K1 , K2 ∈ Th ; F = K1 K2 6= ∅ then F is a common face of K1 and K2 ; Interfaces are collected in the set Fhin , and boundary faces are collected in the set Fhb . Henceforth, we set [ Fh = Fhin Fhb .
(8)
where nF is the outward unit normal to F at s. Summing over all cells, integrating by parts a second time, and 1570
separating volume and interface terms, we obtain a global formulation. Upon introducing the bilinear form N
up
ah (Ih (ξ, .), wh ) =
β X
X Z
wl
K∈Th
l=1
Z
l
l
−
l
Ih (ξ, s)wh ds K
φ1 , φ2 , ..., φN t . Let WN t be an N t -dimensional subspace of Wh with basis functions φi where i ∈ {1, ..., N t }. We assumed that N t = N +m+q where N is the cardinal ¯ \ ∂Ω. To construct a finite element of the set of points on Ω approximation of the temperature T , we introduce:
l
Ih (ξ, s)(βl .∇wh )ds ZK
+
Th (ξ, s) =
⊕ l l (βl .n) Ih (ξ, s)wh dΓ(s)
∂Ω
X
+
X
+
l
l
(βl .nF )[[Ih (ξ, s)]]{{wh }}dΓ(s)
Then, a semi-discrete finite element approximation of the weak formulation (11)- (12), we give the following matrix form:
F
Z
F ∈F in h
F
1 l l |βl .nF |[[Ih (ξ, s)]][[wh ]]dΓ(s), 2
for all (ξ, wh ) ∈ [0, τ ] × Vh , where for a real number x, we define its positive and negative parts respectively as x⊕ := 1 := 21 (|x| − x). By definition, both of the 2 (|x| + x), x quantities are nonnegative. The following linear form lh : Vh −→ R is defined by N
l(wh ) =
β X
Z
l
wl
Z
l
Mh T˙h = Ah Th + Bh Uh + Dh Vh + Ψh (Th ), where Mh , Ah ∈ MN (R), Bh ∈ MN,m (R) and Dh ∈ MN,q (R) are the matrices of finite element method. Mh is a symmetric positive definite matrix and Ah is a negative definite matrix. The matrices Ah , Mh have the following structure:
(βl .n) g(ξ, s, βl )wh dΓ(s)
z(ξ, s)wh ds +
Ah = Tri-Diag(I, A1h , I), Mh = Tri-Diag(Mh1 , Mh2 , Mh1 )
∂Ω
Ω
l=1
for all (ξ, wh ) ∈ [0, τ ] × Vh . The discrete ordinate DG problem is written as: Find Ih (ξ, .) ∈ Vh aup h (Ih (ξ, .), wh )
Thi (ξ)φi (s).
i=1
Z
F ∈F in h
N X
such as
= lh (wh ) for all wh ∈ Vh .
For the implementation of DG methods, matrix assembly and the selection of basis functions in broken polynomial spaces, see [13], [9]. B. Nonlinear heat equation
and
A1h = Tri-Diag(1, −4, 1),
Mh1 is a symmetric positive definite matrix, A1h ∈ Mn (R), Mh2 ∈ Mn (R) and I ∈ Mn (R) is the identity matrix. Ψh ∈ RN is the nonlinear function. Uh ∈ Rn and Vh ∈ Rq are the input vectors of boundary condition. Given an initial condition T0 , at each time step: first we solve the RTE with the DG method, second we compute the the incident radiation intensity G. Finally we solve the nonlinear heat equation using the finite element method. IV. O BSERVER SYNTHESIS METHOD
Let (., .)Ω and (., .)∂Ω be the inner product in L2 (Ω) and L2 (∂Ω) respectively. To consider a weak solution of the equation (3a), the following space is introduced :
We consider the approximate nonlinear system of the coupled radiative-conductive PDE (??)-(3a). ( V = {T ∈ L2 (0, τ, H 1 (Ω)) such that T (ξ, s) = f (ξ, s) for (ξ, s) ∈ [0, τ ]×∂Ω}, Mh T˙h = Ah Th + Bh Uh + Dh Vh + Ψh (Th ) (9) (13) Yh = Ch Th , where f is defined by f (ξ, s) =
u(ξ, s) v(ξ, s)
for(ξ, s) ∈ [0, τ ] × ∂Ω1 for(ξ, s) ∈ [0, τ ] × ∂Ω2
(10)
where ∂Ω = ∂Ω1 ∪ ∂Ω2 . The weak formulation of equation (3a)-(3b) for Dirichlet boundary conditions (4) reads: Find T ∈ V such that (
for all ϕ ∈
dT , ϕ)Ω + (∇T, ∇ϕ)Ω = (Ψ, ϕ)Ω , dξ
H01 (Ω)and
∀ϕ ∈ L2 (Ω).
Consider the following classical Luenberger observer:
(12)
The finite element method is used to solve the equation (11)(12). Let Wh = {vh ∈ C 0 (Ω)/ ∀K ∈ Th , vh |
K
A. LPV- Formulation for radiative-conduction systems
(11)
satisfying the initial condition:
(T (0, s), ϕ)Ω = (T0 , ϕ)Ω
where Th ∈ RN is the state vector, Yh ∈ Rn is the output vector. Ch ∈ Mn,N (R) is the output matrix such that the pair (Ah , Ch ) is detectable.
˙ Mh Tˆh = Ah Tˆh + Bh Uh + Dh Vh + Ψh (Tˆh ) + L(Yh − Ch Tˆh ). (14) where Tˆh represents the estimate of Th and L is to be designed so that the estimation error e = Th − Tˆh converges asymptotically towards zero. The dynamic of the estimation error can be described by
∈ Pk (K)}, 1
Mh e˙ = Ah e + [Ψh (Th ) − Ψh (Tˆh )] − LCh e.
be a family of finite element subspaces of H (Ω). We denote by E1 and E2 the set of points on the edge ∂Ω1 and ∂Ω2 (Card(E1 )=m, Card(E2 )=q) respectively. Let there be given a ¯ and finite element basis functions mesh of N t nodes on Ω
(15)
The function Ψh is γΨ -Lipschitz with respect to its argument, i.e:
1571
kΨh (x) − Ψh (z)k 6 γΨ kx − zk
∀x, z ∈ [0, 1]N .
From the expression of the dimensionless radiative source (??), we have ∂Ψ = −4θT 3 . (16) ∂T Hence, Ψh satisfies ∂Ψih (x) 6 γi where γi 6 γi < 0 ∂xi ∂Ψih (x) = 0 ∀ i, j ∈ {1, ..., N } and i 6= j. ∂xj
Form the convexity principle [3], we deduce that the inequality ( 24) holds if LMI (21)is verified for all ζ ∈ VN , this means that there exist an observer gain L such that V˙ < 0, ∀e 6= 0. Proposition 4.2: Let γ1 =
γi 6
(17)
According the differential mean value theorem DMVT presented in [14], there exist z ∈ Co(Th , Tˆh ), for all i ∈ {1, ..., N } such that: Ψh (Th ) − Ψh (Tˆh ) = D(ζ)e where D(ζ) =
N X N X
ij Ψij h Hh ,
Hhij
=
ζ=
if there exist the observer gain matrix L of appropriate dimension such that the following LMI conditions hold: (Ah + Dγ1 )T + (Ah + Dγ1 ) − CTh LT − LCh < 0.
(26)
(19)
(Ah + Dγ1 )T + (Ah + Dγ1 ) − CTh LT − LCh < 0. From the structure of D(ζ) and (17), we conclude that L verifies (21) ∀ζ ∈ VN . Remark 1: Consider the following gain observer
∀i, j ∈ {1, ..., N }.
eN (i)eTN (j) ∀i, j ∈ {1, ..., N } 12 1N N1 NN (Ψ11 h , Ψh , ..., Ψh , ..., Ψh , ..., Ψh ),
(25)
(18)
and ∂Ψih (z) ∂xj
(γ i ),
where Dγ1 = Diag(γ1 , ..., γ1 ). Then, the observer gain matrix L verified (21) ∀ζ ∈ VN . Proof: D(ζ) is a negative definitive matrix ∀ζ ∈ VN , then if L verified
i=1 j=1
Ψij h =
max
i∈{1,...,N }
L = αCTh
(20)
(27)
where α is the definite positive matrix. This choice of gain verifies the LMI’s (21)-(26).
D(ζ) is a semi-definite diagonal matrix. Consequently, the dynamics (15) can be rewritten as Mh e˙ = [Ah + D(ζ) − LCh ] e. From the inequality (17), it follows that the matrix parameter ζ belongs to a bounded convex set VN defined by
V. E XTENSION TO H∞ PERFORMANCE ANALYSIS In this section, we extend results of the previous section to the observer synthesis based on H∞ performance analysis. We consider also the system described by ( Mh T˙h = Ah Th + Bh Uh + Dh Vh + Ψh (Th ) + W1 w(t)
NN ii VN = {(Ψ11 h , ..., Ψh ), Ψh ∈ {γi , γi } ∀ i ∈ {1, ..., N }}.
Yh = Ch Th + W2 w(t).
B. Observer synthesis Now, we state the following result that provides LMI conditions for the nonlinear systems observer design using the DMVT theorem. Theorem 4.1: The observer (14) is asymptotically convergent if there exist the observer gain matrix L verifies: (Ah +D(ζ))T +(Ah +D(ζ))−CTh LT −LCh < 0, ∀ζ ∈ VN (21)
(28) Such that, W1 and W2 are constant matrices of appropriate dimensions and ω(t) ∈ L2 (Rs ) is a bounded disturbance vector. Now, we consider the following Lunberger state observer: ˙ Mh Tˆh = Ah Tˆh + Bh Uh + Dh Vh + Ψh (Tˆh ) + L(Yh − Ch Tˆh ). (29) The dynamics of estimation error e = Th − Tˆh are given by Mh e˙ = (Ah + D(ζ) − LCh )e + (W1 − LW2 )w.
Proof: Consider the standard Lyapunov function V (e) = eT Mh e > 0. By derivating the Lyapunov function along the trajectories of (15), we obtain V˙ = eT Ξ(ζ)e
Given the system (28) and the observer (29), the problem of robust H∞ observer design is equivalent to find the matrix L such that the estimation error converges H∞ asymptotically to zero, i.e., lim e(ξ) = 0 for w(ξ) = 0
(22)
ξ→∞
(30a)
ke(ξ)k2 6 λkw(ξ)k2 for w(ξ) 6= 0 and e(0) = 0. (30b)
where Ξ(ζ) = (Ah + D(ζ) − LCh )T + (Ah + D(ζ) − LCh ). (23) Consequently, V˙ < 0 for all ζ 6= 0 if Ξ(ζ) < 0
∀ζ ∈ VN .
with λ > 0 is a prescribed scalar disturbance attenuation level. However, to satisfy (30a) and (30b), it is sufficient to find [15] a Lyapunov function V so that
(24) 1572
V˙ + eT e − λ2 wT w < 0
(31)
The results of observers design with DMVT theorem 4.1 given in the last section was extended in the H∞ analysis performances. Theorem 5.1: For a prescribed λ > 0, the H∞ filtering design problem corresponding to the system (28) and the observer (29) is solvable, if there exist the observer gain matrix L of appropriate dimensions such that the following LMI is feasible: Bloc-Diag [Θ(ζ1 , λ), .., Θ(ζ2N , λ)] < 0.
4ξ = 5.10−4 . The reference temperature Tref is assumed equal to 1000K in this case. In addition, the iterations have been terminated when the steady state conditions have been achieved with a tolerance set to 10−5 . In the practice, we used Newton’s method in order to treat the non-linearity. All calculations were performed on a Mac OS X version 10.6.6 2.22GHz Intel Core 2 Duo and using the MATLAB software. A. Observer synthesis
ζi ∈ VN . (32)
Numerical example is given in order to show the validity of the numerical analysis and efficiency of the proposed T observer design techniques. We assumed that the thermal T (Ah + D(ζ) − CT W1 − LW2 h L + Ah + D(ζ) − LCh + I Θ(ζ, λ) = boundary condition (10) is given by T 2 (W1 − LW2 ) −λ I ( (33) u=1 [0, τ ] × ∂Ω1 f= v = 0.5 [0, τ ] × ∂Ω2 , Proof: Let us consider the following Lyapunov function: N is equal to n2 and m = n. The boundary ∂Ω1 is V (e) = eT Mh e, for all i ∈ {1, ..., 2N } where
hence, V˙ (e(ξ)) + eT (ξ)e(ξ) − λ2 ω T (ξ)ω(ξ) T M(L, ζ) (W1 − LW2 ) e(ξ) e(ξ) = ω(ξ) ω(ξ) (W1 − LW2 )T −λ2 I where M(L, ζ) = (Ah + D(ζ) − LCh )T + (Ah + D(ζ) − LCh ) + I. The temperature distribution for Ns = 1
Fig. 1.
Then,
V˙ (e(ξ)) + eT (ξ)e(ξ) − λ2 ω T (ξ)ω(ξ) T e(ξ) e(ξ) = Θ(ζ, λ) ω(ξ) ω(ξ)
where Θ(ζ, λ) is given by (33) which is identical to (32). Consequently, we deduce that under the condition (32), the estimation error converges robustly asymptotically to zero. This concludes the proof of Theorem 5.1. Proposition 5.2: Let γ1 given by (25), if there exist the observer gain matrix L of appropriate dimension such that the following LMI conditions hold: T T (Ah + Dγ1 − CT h L + Ah + Dγ1 − LCh + I (W1 − LW2 )T
W1 − LW2
State observer design for non linear coupled partial differential equations with application to radiative-conductive heat transfer systems Mohamed Ghattassi1 , Mohamed Boutayeb2 and Jean R. Roche3 Abstract— This contribution deals with state observer design for a class of non linear coupled PDE that describe radiativeconductive heat transfer systems. This approach uses first a stable spatial discretization technique that is the Galerkin method to obtain a large scale but finite dimensional system in a suitable form. Thanks to the special structure of the obtained state system, the second main result is to show through the differential mean value theorem (DMVT) that there always exists an observer gain matrix that assures asymptotic convergence. On the other hand, in order to avoid high computational requirements, we show how to construct the observer gain matrix so that the stability condition, written in terms of linear matrix inequality, is satisfied. Extension to H∞ performance analysis is also proposed. In order to show high accuracy of the proposed technique, a numerical example is provided.
I. INTRODUCTION Control design and state observers for non linear partial differential equations (PDE) is known to be a hard task and is still being a challenging problem in the literature. Motivated by this fact, a typical approach to the state observers of PDE systems is to obtain an approximate ordinary differential equation (ODE) representation of the original PDE system by utilizing the spatial discretization techniques [4], [2], [7] which is then used for observer design purposes by applying different existing ODE-based linear or nonlinear observer techniques. This approach is also referred to the reduce-thendesign (RTD) approach. In this work, we consider the nonlinear coupled PDE systems described by the radiative-conduction heat transfer phenomena. This nonlinear problem has not been tackled in the literature nor in infinite dimension or finite dimensional. The model reduction is based on Galerkin (DG) method . The radiative transfer equation (RTE) was solved using the Discontinuous Galerkin method with an upwind numerical flux. The nonlinear heat equation was discretized using a continuous Galerkin (finite element) method. This numerical method is stable and gives a suitable form for the synthesis of the observer. Various research activities have focused on observer design for nonlinear systems. This topic was motivated by the fact that the state estimation can be used for control and diagnosis. Linear parameter-varying (LPV) techniques
have been developed into effective tools to control multiinput multi-output (MIMO) nonlinear plants. Many nonlinear systems of practical interest can be represented as LPV systems. However, the application of these techniques has not yet been extensively explored in control of PDEs. A LPV approach together with a nonlinear observer was proposed in Hashemi and Werner [8] to control the one-dimensional Burgers equation. The basic idea of this work is to use the well-known DMVT which allows to write the dynamics of the estimation error as a class of LPV systems. Stability of the estimation error is analyzed using the convexity principle and the Lyapunov stability theory [15], [16], [6]. The observer gain guaranteeing the global convergence of the proposed observer is computed by LMI. The structure of the state system gives the advantage of DMVT to reduce the computational complexity and to construct the gain observer which ensures the stability of the proposed observer. This result was extended to the observer synthesis with H∞ performance analysis. This article is organized as follows. Next section develops the governing equations for two-dimensional combined radiative and conductive heat transfer. In section 3, the reduced model of the coupled partial differential equations is given. DMVT for observer design of the reduction model of the coupled partial differential equations are described in section 4. The section 5, an extension to H∞ performance analysis is proposed for the results of the previous section. Simulations results are given in section 6 and the last section draws conclusions. Notations.Throughout this paper, we will use the following notation: • • • • •
•
1 M. Ghattassi is with Institut Elie Cartan de Lorraine, Universit´ e de Lorraine, CNRS, INRIA, B.P. 70239, 54506, Vandoeuvre l`es Nancy, France.
k.k is the Euclidean norm. (?) is used for the blocks induced by symmetry. I represents the identity matrix of appropriate dimension. AT represents the transposed matrix of A. For a square matrix S , S > 0 (S < 0) means that this matrix is positive definite (negative definite). ith z}|{ T es (i) = (0, ..., 0, 1 , 0, ..., 0) ∈ Rs ,s > 1 is a vector | {z } s components of the canonical basis of Rs .
[email protected] 2 M. Boutayeb is with the Centre des Recherches en Automatique, Universit´e de Lorraine, CNRS, B.P. 70239, 54506, Vandoeuvre l`es Nancy, France. [email protected] 3 J.R. Roche is with the Institut Elie Cartan de Lorraine, Universit´ e de Lorraine, CNRS, INRIA, B.P. 70239, 54506, Vandoeuvre l`es Nancy, France.
[email protected] 978-1-4673-6090-6/14/$31.00 ©2014 European Union
II. M ODEL PROBLEM Let Ω be a bounded domain in R2 and D the unit disk. Thus, for the problem considered, β ∈ D = {β ∈ R2 : |β| 6 1}, s? = (s? , s2 ) ∈ Ω and ξ ∈ [0, τ ], for τ > 0. Let n be
1569
the outward unit normal to the boundary ∂Ω. We denote ∂Ω− = {(s, β) ∈ ∂Ω × D such that β.n < 0}. The unknown of the RTE is the dimensionless radiation intensity denoted I(ξ, s, β) given at time ξ, position s and in the direction β. The unknown of the energy equation is the dimensionless temperature T (ξ, s) at time ξ and position s. The RTE is given by (see, [1]) β.∇I(ξ, s, β) + I(ξ, s, β) = T 4 (ξ, s),
Now, we proceed with the spatial discretization of the SN transport equation using the DG method . We consider a triangulation Th of Ω such that [ Ω= K, K∈Th • •
(1)
•
for all (ξ, s, β) ∈ [0, τ ] × Ω × D with the boundary condition I(ξ, s, β) = g(ξ, s, β),
•
(2)
for all (ξ, s, β) ∈ [0, τ ] × ∂Ω− , the radiative boundary condition is given by
Moreover, for any mesh element K ∈ Th , the set
g(t, s, β) = T 4 (ξ, s) for (ξ, s, β) ∈ [0, τ ] × ∂Ω− ,
FK = {F ∈ Fh |F ∈ ∂K},
where T (ξ, s) is the temperature at the boundary of the medium. The dimensionless radiative equation (1) is combined with the following energy equation (nonlinear heat equation) : ∂T − 4T + θT 4 = θG ∂ξ T (0, s) = T0 (s)
for (ξ, s) ∈ [0, τ ] × Ω
(3a)
for all s ∈ Ω
(3b)
where θ is the dimensionless constant and G is the incident radiation intensity, Z G(ξ, s) = I(ξ, s, β) dβ for (ξ, s) ∈ [0, τ ] × Ω.
•
Vh = {Ih ∈ L2 (Ω)/ ∀K ∈ Th , Ih |
K
∈ Pk (K)},
(ξ, s) ∈ [0, τ ] × ∂Ω,
T (ξ, s) = f (ξ, s)
(4)
III. A NALYSIS OF RADIATIVE CONDUCTIVE HEAT TRANSFER
A. Radiative Transfer Equation To compute a numerical solution of the RTE, the classical SN −discrete ordinate method is introduced [12]. The SN method consists of replacing the radiation intensity I(ξ, s, β) by a discrete radiation intensity (I 1 (ξ, s), I 2 (ξ, s), β ..., I N (ξ, s)), where I l (ξ, s) = I(ξ, s, βl ) for all l ∈ {1, ..., N β}. A quadrature rule was chosen {(βl , wl ) l = 1, ..., N β }, such that: β
Z f (β, s)dβ '
N X
β
f (βl , s)wl and
l=1
N X
wl = 4π.
l=1
This quadrature formula preserves a symmetry for the opposite βl0 = −βl of discrete ordinate βl in the quadrature set, the weights are equal wl0 = wl , see [12]. The discrete radiation intensity I l is the solution of a N β system of partial differential equations over Ω: l l 4 βl .∇I (ξ, s) + I (ξ, s) = T (ξ, s)
I l (ξ, s) = g(ξ, s, βl )
N
(6)
and denote a generic element in Vh as Ih (ξ, .) = l Nβ Ih (ξ, .) l=1 . Due to the discontinuous nature of the spatial approximation, functions Ihl ∈ Vh are double-valued on interior faces. Consider an interior face F ∈ Fhin separating two mesh cells, K1 and K2 . The mean value and jump of a function Ihl (ξ, .) ∈ Vh are defined as follows: 1 l (I (ξ, .) + I2l (ξ, .)), 2 1 [[Ihl (ξ, .)]] = (I1l (ξ, .) − I2l (ξ, .)),
{{Ihl (ξ, .)}} =
for all ξ ∈ [0, τ ], where I1l (ξ, .) = Ihl (ξ, .)|K and I2l (ξ, .) = 1 Ihl (ξ, .)|K are the restrictions of Ihl (ξ, .) on the mesh cells 2 K1 and K2 , respectively. The DG formulation is obtained by multiplying the SN equation for the direction βl with the test function wh ∈ Vh and applying the upwind numerical ˆ fluxes F(s) to approximate the quantity (βl .n)I l (ξ, .) on the elements boundary K: Z K
κIhl (ξ, .)wh − (βl .∇h wh )Ihl (ξ, .) +
Z ∂K
Fˆ (s)wh =
Z z(ξ, .)wh , Ω
(7)
ˆ ∀ξ ∈ [0, τ ]. The upwind numerical flux F(s) at the mesh interface F from K1 to K2 is given by:
(ξ, s) ∈ [0, τ ] × Ω
1 Fˆ (s) = βl .nF {{I l (ξ, s)}} + |βl .nF |[[I l (ξ, s)]], ∀ξ ∈ [0, τ ], 2
(ξ, s) ∈ [0, τ ] × ∂Ωl− ,
∂Ωl− = {s ∈ ∂Ω such that βl .n < 0}.
(5)
where Pk (K) denotes the set of polynomial defined in K of degree less than or equal to k. Define: β Vh = {Vh }l=1 ,
We consider the Dirichlet boundary conditions, where the dimensionless temperature at the boundary is known
where
collects the mesh faces composing the boundary of K. diam(K) = h for each K ∈ Th (structured mesh).
Now, an approximation space Vh ⊂ L2 (Ω) is introduced such that
D
D
Each K is a triangle with nonempty internal; ◦ T ◦ K1 K2 T = ∅ for each distinct K1 , K2 ∈ Th ; F = K1 K2 6= ∅ then F is a common face of K1 and K2 ; Interfaces are collected in the set Fhin , and boundary faces are collected in the set Fhb . Henceforth, we set [ Fh = Fhin Fhb .
(8)
where nF is the outward unit normal to F at s. Summing over all cells, integrating by parts a second time, and 1570
separating volume and interface terms, we obtain a global formulation. Upon introducing the bilinear form N
up
ah (Ih (ξ, .), wh ) =
β X
X Z
wl
K∈Th
l=1
Z
l
l
−
l
Ih (ξ, s)wh ds K
φ1 , φ2 , ..., φN t . Let WN t be an N t -dimensional subspace of Wh with basis functions φi where i ∈ {1, ..., N t }. We assumed that N t = N +m+q where N is the cardinal ¯ \ ∂Ω. To construct a finite element of the set of points on Ω approximation of the temperature T , we introduce:
l
Ih (ξ, s)(βl .∇wh )ds ZK
+
Th (ξ, s) =
⊕ l l (βl .n) Ih (ξ, s)wh dΓ(s)
∂Ω
X
+
X
+
l
l
(βl .nF )[[Ih (ξ, s)]]{{wh }}dΓ(s)
Then, a semi-discrete finite element approximation of the weak formulation (11)- (12), we give the following matrix form:
F
Z
F ∈F in h
F
1 l l |βl .nF |[[Ih (ξ, s)]][[wh ]]dΓ(s), 2
for all (ξ, wh ) ∈ [0, τ ] × Vh , where for a real number x, we define its positive and negative parts respectively as x⊕ := 1 := 21 (|x| − x). By definition, both of the 2 (|x| + x), x quantities are nonnegative. The following linear form lh : Vh −→ R is defined by N
l(wh ) =
β X
Z
l
wl
Z
l
Mh T˙h = Ah Th + Bh Uh + Dh Vh + Ψh (Th ), where Mh , Ah ∈ MN (R), Bh ∈ MN,m (R) and Dh ∈ MN,q (R) are the matrices of finite element method. Mh is a symmetric positive definite matrix and Ah is a negative definite matrix. The matrices Ah , Mh have the following structure:
(βl .n) g(ξ, s, βl )wh dΓ(s)
z(ξ, s)wh ds +
Ah = Tri-Diag(I, A1h , I), Mh = Tri-Diag(Mh1 , Mh2 , Mh1 )
∂Ω
Ω
l=1
for all (ξ, wh ) ∈ [0, τ ] × Vh . The discrete ordinate DG problem is written as: Find Ih (ξ, .) ∈ Vh aup h (Ih (ξ, .), wh )
Thi (ξ)φi (s).
i=1
Z
F ∈F in h
N X
such as
= lh (wh ) for all wh ∈ Vh .
For the implementation of DG methods, matrix assembly and the selection of basis functions in broken polynomial spaces, see [13], [9]. B. Nonlinear heat equation
and
A1h = Tri-Diag(1, −4, 1),
Mh1 is a symmetric positive definite matrix, A1h ∈ Mn (R), Mh2 ∈ Mn (R) and I ∈ Mn (R) is the identity matrix. Ψh ∈ RN is the nonlinear function. Uh ∈ Rn and Vh ∈ Rq are the input vectors of boundary condition. Given an initial condition T0 , at each time step: first we solve the RTE with the DG method, second we compute the the incident radiation intensity G. Finally we solve the nonlinear heat equation using the finite element method. IV. O BSERVER SYNTHESIS METHOD
Let (., .)Ω and (., .)∂Ω be the inner product in L2 (Ω) and L2 (∂Ω) respectively. To consider a weak solution of the equation (3a), the following space is introduced :
We consider the approximate nonlinear system of the coupled radiative-conductive PDE (??)-(3a). ( V = {T ∈ L2 (0, τ, H 1 (Ω)) such that T (ξ, s) = f (ξ, s) for (ξ, s) ∈ [0, τ ]×∂Ω}, Mh T˙h = Ah Th + Bh Uh + Dh Vh + Ψh (Th ) (9) (13) Yh = Ch Th , where f is defined by f (ξ, s) =
u(ξ, s) v(ξ, s)
for(ξ, s) ∈ [0, τ ] × ∂Ω1 for(ξ, s) ∈ [0, τ ] × ∂Ω2
(10)
where ∂Ω = ∂Ω1 ∪ ∂Ω2 . The weak formulation of equation (3a)-(3b) for Dirichlet boundary conditions (4) reads: Find T ∈ V such that (
for all ϕ ∈
dT , ϕ)Ω + (∇T, ∇ϕ)Ω = (Ψ, ϕ)Ω , dξ
H01 (Ω)and
∀ϕ ∈ L2 (Ω).
Consider the following classical Luenberger observer:
(12)
The finite element method is used to solve the equation (11)(12). Let Wh = {vh ∈ C 0 (Ω)/ ∀K ∈ Th , vh |
K
A. LPV- Formulation for radiative-conduction systems
(11)
satisfying the initial condition:
(T (0, s), ϕ)Ω = (T0 , ϕ)Ω
where Th ∈ RN is the state vector, Yh ∈ Rn is the output vector. Ch ∈ Mn,N (R) is the output matrix such that the pair (Ah , Ch ) is detectable.
˙ Mh Tˆh = Ah Tˆh + Bh Uh + Dh Vh + Ψh (Tˆh ) + L(Yh − Ch Tˆh ). (14) where Tˆh represents the estimate of Th and L is to be designed so that the estimation error e = Th − Tˆh converges asymptotically towards zero. The dynamic of the estimation error can be described by
∈ Pk (K)}, 1
Mh e˙ = Ah e + [Ψh (Th ) − Ψh (Tˆh )] − LCh e.
be a family of finite element subspaces of H (Ω). We denote by E1 and E2 the set of points on the edge ∂Ω1 and ∂Ω2 (Card(E1 )=m, Card(E2 )=q) respectively. Let there be given a ¯ and finite element basis functions mesh of N t nodes on Ω
(15)
The function Ψh is γΨ -Lipschitz with respect to its argument, i.e:
1571
kΨh (x) − Ψh (z)k 6 γΨ kx − zk
∀x, z ∈ [0, 1]N .
From the expression of the dimensionless radiative source (??), we have ∂Ψ = −4θT 3 . (16) ∂T Hence, Ψh satisfies ∂Ψih (x) 6 γi where γi 6 γi < 0 ∂xi ∂Ψih (x) = 0 ∀ i, j ∈ {1, ..., N } and i 6= j. ∂xj
Form the convexity principle [3], we deduce that the inequality ( 24) holds if LMI (21)is verified for all ζ ∈ VN , this means that there exist an observer gain L such that V˙ < 0, ∀e 6= 0. Proposition 4.2: Let γ1 =
γi 6
(17)
According the differential mean value theorem DMVT presented in [14], there exist z ∈ Co(Th , Tˆh ), for all i ∈ {1, ..., N } such that: Ψh (Th ) − Ψh (Tˆh ) = D(ζ)e where D(ζ) =
N X N X
ij Ψij h Hh ,
Hhij
=
ζ=
if there exist the observer gain matrix L of appropriate dimension such that the following LMI conditions hold: (Ah + Dγ1 )T + (Ah + Dγ1 ) − CTh LT − LCh < 0.
(26)
(19)
(Ah + Dγ1 )T + (Ah + Dγ1 ) − CTh LT − LCh < 0. From the structure of D(ζ) and (17), we conclude that L verifies (21) ∀ζ ∈ VN . Remark 1: Consider the following gain observer
∀i, j ∈ {1, ..., N }.
eN (i)eTN (j) ∀i, j ∈ {1, ..., N } 12 1N N1 NN (Ψ11 h , Ψh , ..., Ψh , ..., Ψh , ..., Ψh ),
(25)
(18)
and ∂Ψih (z) ∂xj
(γ i ),
where Dγ1 = Diag(γ1 , ..., γ1 ). Then, the observer gain matrix L verified (21) ∀ζ ∈ VN . Proof: D(ζ) is a negative definitive matrix ∀ζ ∈ VN , then if L verified
i=1 j=1
Ψij h =
max
i∈{1,...,N }
L = αCTh
(20)
(27)
where α is the definite positive matrix. This choice of gain verifies the LMI’s (21)-(26).
D(ζ) is a semi-definite diagonal matrix. Consequently, the dynamics (15) can be rewritten as Mh e˙ = [Ah + D(ζ) − LCh ] e. From the inequality (17), it follows that the matrix parameter ζ belongs to a bounded convex set VN defined by
V. E XTENSION TO H∞ PERFORMANCE ANALYSIS In this section, we extend results of the previous section to the observer synthesis based on H∞ performance analysis. We consider also the system described by ( Mh T˙h = Ah Th + Bh Uh + Dh Vh + Ψh (Th ) + W1 w(t)
NN ii VN = {(Ψ11 h , ..., Ψh ), Ψh ∈ {γi , γi } ∀ i ∈ {1, ..., N }}.
Yh = Ch Th + W2 w(t).
B. Observer synthesis Now, we state the following result that provides LMI conditions for the nonlinear systems observer design using the DMVT theorem. Theorem 4.1: The observer (14) is asymptotically convergent if there exist the observer gain matrix L verifies: (Ah +D(ζ))T +(Ah +D(ζ))−CTh LT −LCh < 0, ∀ζ ∈ VN (21)
(28) Such that, W1 and W2 are constant matrices of appropriate dimensions and ω(t) ∈ L2 (Rs ) is a bounded disturbance vector. Now, we consider the following Lunberger state observer: ˙ Mh Tˆh = Ah Tˆh + Bh Uh + Dh Vh + Ψh (Tˆh ) + L(Yh − Ch Tˆh ). (29) The dynamics of estimation error e = Th − Tˆh are given by Mh e˙ = (Ah + D(ζ) − LCh )e + (W1 − LW2 )w.
Proof: Consider the standard Lyapunov function V (e) = eT Mh e > 0. By derivating the Lyapunov function along the trajectories of (15), we obtain V˙ = eT Ξ(ζ)e
Given the system (28) and the observer (29), the problem of robust H∞ observer design is equivalent to find the matrix L such that the estimation error converges H∞ asymptotically to zero, i.e., lim e(ξ) = 0 for w(ξ) = 0
(22)
ξ→∞
(30a)
ke(ξ)k2 6 λkw(ξ)k2 for w(ξ) 6= 0 and e(0) = 0. (30b)
where Ξ(ζ) = (Ah + D(ζ) − LCh )T + (Ah + D(ζ) − LCh ). (23) Consequently, V˙ < 0 for all ζ 6= 0 if Ξ(ζ) < 0
∀ζ ∈ VN .
with λ > 0 is a prescribed scalar disturbance attenuation level. However, to satisfy (30a) and (30b), it is sufficient to find [15] a Lyapunov function V so that
(24) 1572
V˙ + eT e − λ2 wT w < 0
(31)
The results of observers design with DMVT theorem 4.1 given in the last section was extended in the H∞ analysis performances. Theorem 5.1: For a prescribed λ > 0, the H∞ filtering design problem corresponding to the system (28) and the observer (29) is solvable, if there exist the observer gain matrix L of appropriate dimensions such that the following LMI is feasible: Bloc-Diag [Θ(ζ1 , λ), .., Θ(ζ2N , λ)] < 0.
4ξ = 5.10−4 . The reference temperature Tref is assumed equal to 1000K in this case. In addition, the iterations have been terminated when the steady state conditions have been achieved with a tolerance set to 10−5 . In the practice, we used Newton’s method in order to treat the non-linearity. All calculations were performed on a Mac OS X version 10.6.6 2.22GHz Intel Core 2 Duo and using the MATLAB software. A. Observer synthesis
ζi ∈ VN . (32)
Numerical example is given in order to show the validity of the numerical analysis and efficiency of the proposed T observer design techniques. We assumed that the thermal T (Ah + D(ζ) − CT W1 − LW2 h L + Ah + D(ζ) − LCh + I Θ(ζ, λ) = boundary condition (10) is given by T 2 (W1 − LW2 ) −λ I ( (33) u=1 [0, τ ] × ∂Ω1 f= v = 0.5 [0, τ ] × ∂Ω2 , Proof: Let us consider the following Lyapunov function: N is equal to n2 and m = n. The boundary ∂Ω1 is V (e) = eT Mh e, for all i ∈ {1, ..., 2N } where
hence, V˙ (e(ξ)) + eT (ξ)e(ξ) − λ2 ω T (ξ)ω(ξ) T M(L, ζ) (W1 − LW2 ) e(ξ) e(ξ) = ω(ξ) ω(ξ) (W1 − LW2 )T −λ2 I where M(L, ζ) = (Ah + D(ζ) − LCh )T + (Ah + D(ζ) − LCh ) + I. The temperature distribution for Ns = 1
Fig. 1.
Then,
V˙ (e(ξ)) + eT (ξ)e(ξ) − λ2 ω T (ξ)ω(ξ) T e(ξ) e(ξ) = Θ(ζ, λ) ω(ξ) ω(ξ)
where Θ(ζ, λ) is given by (33) which is identical to (32). Consequently, we deduce that under the condition (32), the estimation error converges robustly asymptotically to zero. This concludes the proof of Theorem 5.1. Proposition 5.2: Let γ1 given by (25), if there exist the observer gain matrix L of appropriate dimension such that the following LMI conditions hold: T T (Ah + Dγ1 − CT h L + Ah + Dγ1 − LCh + I (W1 − LW2 )T
W1 − LW2
