Backstepping observers for linear PDEs on higher-dimensional spatial domains

Automatica 51 (2015) 85–97

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Backstepping observers for linear PDEs on higher-dimensional spatial domains✩ Lukas Jadachowski a , Thomas Meurer b , Andreas Kugi a a

Christian Doppler Laboratory for Model-Based Process Control in the Steel Industry, Automation and Control Institute, Vienna University of Technology, Austria b Chair of Automatic Control, Kiel University, Germany

article

info

Article history: Received 24 September 2013 Received in revised form 16 June 2014 Accepted 3 October 2014

Keywords: Backstepping Parabolic PDEs Observer design Higher-dimensional spatial domain

abstract The extension of the backstepping-based state observer design is considered for systems governed by parabolic PDEs defined on a higher-dimensional cuboidal spatial domain with the measured output being restricted to a single boundary hyperplane. Thereby, a suitable coordinate and state transformation is employed, which allows for the application of the backstepping method for the determination of a Luenberger-type state observer with observer corrections entering the PDE and the boundary condition. This ensures that the observer error dynamics follows the behavior of a predefined exponentially stable target system. In view of a practical realization of the proposed observer approach, the idealized system output, previously assumed as measurable on the entire boundary hyperplane, is reconstructed from a limited set of optimally located measurements by means of a least squares approximation. The observer error convergence and the applicability of the proposed approach are confirmed analytically and are evaluated in numerical simulations. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction Model-based control and advanced process monitoring for distributed-parameter systems (DPSs), i.e., systems governed by partial differential equations (PDEs), usually require full state information. This is evident when full state feedback control strategies (see, e.g., Liu, 2003; Meurer & Kugi, 2009) are utilized or the system state evolution is required to gain detailed insight into the system behavior, e.g., for diagnostic purposes (Akkari, Chevallier, & Boillereaux, 2006; Dochain, Tali-Maamar, & Babary, 1997). However, for DPSs the measurement information is typically limited by the number of sensors and in case of systems with a higher-dimensional spatial domain by the fact that measurements are usually available only on the boundary of the spatial domain (Vande Wouwer & Zeitz, 2001). This motivates the need of

✩ The financial support in part by the Austrian Federal Ministry of Science, Research and Economy, and the National Foundation for Research, Technology and Development is gratefully acknowledged. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Nicolas Petit under the direction of Editor Miroslav Krstic. E-mail addresses: [email protected] (L. Jadachowski), [email protected] (T. Meurer), [email protected] (A. Kugi).

http://dx.doi.org/10.1016/j.automatica.2014.10.108 0005-1098/© 2014 Elsevier Ltd. All rights reserved.

state estimation by means of a state observer. While the so-called early-lumping approach, in which the dynamics of the DPS is approximated by ordinary differential equations is mainly used for system analysis and design, late-lumping methods preserving the infinite-dimensional character of the system during the whole design process gain more and more importance in recent years. In this context, backstepping appears to be a particularly promising and systematic observer design approach for a broad class of systems governed by PDEs, see, e.g., Krstic and Smyshlyaev (2008). Basically, the backstepping method relies on the application of a Volterra integral transformation mapping a predefined exponentially stable target system into the observer error dynamics. For systems governed by parabolic PDEs defined on a 1-dimensional (1D) spatial domain, a systematic observer design approach using boundary sensing is introduced in Smyshlyaev and Krstic (2005). First results on the backstepping-based observer design for higher-dimensional DPSs trace back to Vazquez, Schuster, and Krstic (2008). Here, a state estimator in a double semi-infinite 3-dimensional (3D) domain is presented for a model of magnetohydrodynamic flow and Fourier transform methods are applied to put the system in a form, where a 1D backstepping method is applicable. In Jadachowski, Meurer, and Kugi (2011), a solution of the state estimation problem is presented for a diffusion–reaction system (DRS) defined on a 3D rectangular domain with a spatially

86

L. Jadachowski et al. / Automatica 51 (2015) 85–97

varying reaction parameter. Thereby, assuming an idealized system output restricted to a single boundary surface, the target system is designed in such a way that the backstepping-kernel can be determined similar to the 1D setting. In a second step, the idealized system output is reconstructed by means of a least squares method under the assumption that only a finite number of measurements is available on the output surface. Extending this idea to the most general possible case, the subsequent contribution presents a generalization of the results from Jadachowski et al. (2011) to a diffusion–convection–reaction system (DCRS) defined on an n-dimensional (nD) cuboidal spatial domain with an output restricted to the boundary hyperplane. Here, an appropriate coordinate and state transformation is utilized to reformulate the DCRS in a simpler form, in which backstepping can be directly applied. With this, the determination of the output injections of the state observer ensures that the observer error dynamics corresponds to the exponentially stable target system. Moreover, the impact of general time-varying boundary conditions is addressed as regards the convergence analysis for the solution of the kernel-PDE governing the evolution of the backstepping-kernel. For the approximate realization of the idealized system output, an approach for the optimal sensor location is considered. The proposed observer design approach is completed by a stability and robustness analysis followed by numerical results confirming the exponential decay of the observer error. The paper is organized as follows: In Section 2, the state estimation problem for a multidimensional DCRS with spatially and time-varying parameters and boundary output is formulated. Here, a particular coordinate and state transformation followed by the Hopf–Cole transformation (see Appendix A) are utilized to obtain a simplified equivalent problem formulation. In Section 3, the backstepping-based state observer is proposed for the resulting DRS. The observer gains are determined to enforce the observer error dynamics to follow a suitably prescribed target system. The stability and robustness analysis of the resulting observer error dynamics is addressed in Section 4. In Section 5, an approximate realization of the idealized output is presented by means of a finite number of optimally located sensors. The estimation performance is illustrated in Section 6 based on the example of a DRS with 3D spatial domain. Notation: The following notation is used subsequently: z = (z1 , . . . , zn ), zj = z \ {zj }, zsj = z|zj =s , In := {1, 2, . . . , n}, Inj := {1, 2, . . . , j − 1, j + 1, . . . , n}. 2. Problem formulation Subsequently, the observer design is considered for a linear scalar parabolic PDE with spatially and time-varying parameters given by

∂t x(z, t ) =

n  

ai (zi )∂z2i x(z, t ) + bi (zi )∂zi x(z, t )



i =1

+ c (z, t )x(z, t )

(1)

defined on (z, t ) ∈ × R with an nD rectangular spatial domain Ωzn := {z ∈ Rn | 0 < zi < Li , i = 1, . . . , n} and positive finite constants Li , i = 1, . . . , n. The corresponding boundary conditions (BCs) on ∂ Ωzn are assumed as

Ωzn

+

pi0 (t )∂zi x(z, t ) + pi1 (t )x(z, t ) = 0,

zi = 0

(2)

z i = Li ,

(3)

with Ωzn = Ωzn + ∂ Ωzn . The system output is restricted to the boundary hyperplane ∂ Ωj = {z ∈ Ωzn | zj = 0} and is assumed as y(zj , t ) = r0 (t )∂zj x(z, t ) + r1 (t )x(z, t ),

zj = 0

for a single but arbitrary j ∈ In , t ≥ 0 and parameters j

j

j

(5) j r0

(t ) and

r1 (t ) linearly independent of p0 (t ) and p1 (t ), i.e., satisfying

 det

j

p0 (t ) r0 (t )

j



p1 ( t ) r1 ( t )

̸= 0,

t ≥ 0.

(6)

Note that (6) is required to ensure the observability of the system (1)–(5), since it prevents that (5) vanishes identically for all t ≥ 0 according to (2) with i = j. For results on existence and uniqueness of solutions to (1)–(4), see, e.g., Ladyžhenskaia, Solonnikov, and Ural’ceva (1988). Remark 1. The consideration of the sensor model as defined in (5) leads to an idealized setting with an infinite-dimensional measurement being restricted to a single boundary hyperplane. Although of particular importance in view of the extension of the backstepping observer design to PDEs with higher-dimensional spatial domains, the assumption of an infinite-dimensional output (5) has to be relaxed in practical applications to a finite number of sensors. Therefore, in Section 5 the idealized system output (5) is reconstructed from a limited set of measurements with finite-dimensional spatial sensor characteristics. Note that the notion of a class of Gevrey functions GD,σ (R+ 0 ) as introduced in Appendix C is utilized below. Assumption 2. In view of boundedness and differentiability of the system parameters as well as for the subsequent application of the coordinate and state transformation it is assumed that: (a) the diffusion parameters satisfy ai (z ) ∈ C 2 ([0, Li ]) and 0 < ail ≤ ai (zi ) ≤ air < ∞ for zi ∈ [0, Li ] and all i ∈ In , where ail and air are fixed positive constants. (b) the convection parameters are bounded and fulfill bi (zi ) ∈ C 1 ([0, Li ]) for zi ∈ [0, Li ] and i = 1, . . . , n. (c) the reaction parameter c (z, t ) is bounded for (z, t ) ∈ Ωzn × R+ 0 . (d) the boundary and output parameters are bounded and satisfy pil (t ), qil (t ), rl (t ) ∈ GD,σ (R+ 0 ) for σ ∈ [1, 2], l = 0, 1 and i = 1, . . . , n. Further assumptions on the spatial and differentiability properties of the parameter c (z, t ) are given below. Assumption 3. It is assumed that the reaction parameter c (z, t ) is separable in the coordinate zj and takes the form c (z, t ) = 0 c0 (zj , t )+ cj (zj , t ) with (z, t ) ∈ Ωzn × R+ 0 and cj (zj , t ) ∈ C ([0, Lj ])∩ + GD,σ (R0 ) for σ ∈ [1, 2]. Remark 4. Obviously, the structural property assigned for c (z, t ) by Assumption 3 represents a rather restrictive constraint, which theoretically avoids the application of the following observer design approach to a broader class of parabolic PDEs. However, this is a technical assumption which is used to obtain the stability result of the observer error in the L2 -norm. Moreover, as will be presented in Section 4.3, the robustness and the applicability of the determined state observer with respect to a parameter c (z, t ) not separable according to Assumption 3 can also be ensured provided that an additional condition on the target system is satisfied.

for all i ∈ In , t > 0, while the consistent initial condition (IC) follows as

By making use of a suitable change of coordinates followed by the so-called Hopf–Cole transformation it can be shown that the system (1)–(5) is equivalent to the following DRS with only a single spatially and time-varying reaction parameter, i.e.,

x(z, 0) = x0 (z),

∂t x(z, t ) = 1x(z, t ) + c (z, t )x(z, t )

qi0

(t )∂zi x(z, t ) +

qi1

(t )x(z, t ) = 0,

z ∈ Ωzn

(4)

(7)

L. Jadachowski et al. / Automatica 51 (2015) 85–97

defined on (z, t ) ∈ Ωzn × R+ with the Laplace operator ∆ = n 2 i=1 ∂zi , corresponding BCs pi0 qi0

(t )∂zi x(z, t ) +

pi1

(t )x(z, t ) = 0,

zi = 0

(8)

(t )∂zi x(z, t ) +

qi1

(t )x(z, t ) = 0,

zi = Li

(9)

for all i ∈ In , t > 0, consistent IC x(z, 0) = x0 (z),

z∈

Ωzn

(10)

zj = 0

(11)

for t ≥ 0. Thereby, the transformed parameters c (z, t ), (t ), (t ), rl (t ), l = 0, 1 in (7)–(11) depend on the parameters from (1)– (5). The necessary computations are summarized in Appendix A. Here, for the sake of readability the bar on the symbols is omitted, coordinates ζi and ζ are replaced by zi and z, respectively, and Li = L¯ i . pil

qil

Remark 5. As an example of a real system where Assumption 3 is satisfied, consider a two-dimensional for a compos  heat equation ite material, i.e., ρ cp ∂t T (z1 , z2 , t ) = i=1,2 ∂zi λi (zi )∂zi T (z1 , z2 , t ) , where ρ denotes the density, cp the heat capacity, λi (zi ), i = 1, 2 the space-dependent thermal conductivities, and T (z1 , z2 , t ) the temperature. This PDE corresponds to (1) with n = 2, ai (zi ) = λi (zi )/(ρ cp ), bi (zi ) = ∂zi λi (zi )/(ρ cp ) and c (z, t ) ≡ 0. With the Hopf–Cole transformation (A.2), see Appendix A, the introduced heat equation can be transformed into (A.3) with a separable reaction parameter c¯ (ζ1 , ζ2 ) = c¯1 (ζ1 ) + c¯2 (ζ2 ) satisfying Assumption 3, where c¯i (ζi ) = −b¯ 2i (ζi ) − ∂ζi b¯ i (ζi ) with b¯ i (ζi ) = √  ∂zi λi (zi ) /(2ρ 2 cp2 ) evaluated for zi = z¯i−1 (ζi ) and z¯i (zi ) = − 12

3. Observer design Following the idea of the backstepping-based observer design, see, e.g., Smyshlyaev and Krstic (2005) and extending the results from Jadachowski et al. (2011) to the nD setting, the distributedparameter Luenberger state observer in the observer state xˆ (z, t ) for the system (7)–(11) is formulated in the form, i.e.,

∂t xˆ (z, t ) = 1xˆ (z, t ) + c (z, t )ˆx(z, t )   + l(zj , t ) y(zj , t ) − yˆ (zj , t ) ,

(t )∂zj xˆ (z, t ) +

j p1

zi = 0, ∀i ∈ Inj

(t )ˆx(z, t )

= l0 (t )(y(zj , t ) − yˆ (zj , t )), qi0 (t )∂zi xˆ (z, t ) + qi1 (t )ˆx(z, t ) = 0,

zi = Li , ∀i ∈ In

(20) (21)

for t > 0 and z ∈ Ωzn .

(22)

In the following, the output injection weights l(zj , t ), l0 (t ) are determined in such a way that the resulting observer error dynamics (18)–(22) is exponentially stable in the L2 -norm. For this, proceeding along the lines of the backstepping method the specification of a target system prescribing the evolution of the observer error dynamics is required. 3.1. Selection of the target system In the course of the backstepping-based observer design approach the appropriate specification of the target system is one of the crucial parts. In the following, the target system with an nD spatial domain is chosen as

∂t w(z, t ) = 1w(z, t ) − µ(z, t )w(z, t )

(23)

defined on (z, t ) ∈ Ωzn × R+ with corresponding BCs w,i

w,i

w,i

w,i

p0 ∂zi w(z, t ) + p1 w(z, t ) = 0, q0 ∂zi w(z, t ) + q1 w(z, t ) = 0,

zi = 0

(24)

zi = Li

(25)

for all i ∈ In , t > 0 and IC

w(z, 0) = w0 (z),

z ∈ Ωzn .

(26)

zi = Li , ∀i ∈ In

z ∈ Ωzn .

Lemma 6 (Stability in L2 -Norm). The parabolic PDE (23)–(26) is exponentially stable in the L2 -norm, if inf µ(z, t ) + λ ≥ ϵ,

z∈Ωzn

∀t ∈ R +

(27)

is satisfied for some ϵ > 0, where λ denotes the smallest eigenvalue of the Sturm–Liouville problem 1w(z, t ) + λw(z, t ) = 0 with BCs (24) and (25).

(14)

V˙ (t ) =

zj = 0

 Ωzn

1w(z, t )w(z, t )dz −

 (15)

≤ −λ Ωzn

(16)

The output estimate for t ≥ 0 is given by yˆ (zj , t ) = r0 (t )∂zj xˆ (z, t ) + r1 (t )ˆx(z, t ),

zj = 0,

(13)

for t > 0 with xˆ (z, 0) = xˆ 0 (z),

(19)

Proof. The proof of Lemma 6 directly follows from the consideration of the change of the positive definite Lyapunov functional V (t ) = 12 ∥w(z, t )∥2L2 together with the Rayleigh principle (see, e.g., Temple, 1928). The rate of change of V (t ) along a solution trajectory of (23)–(26) yields

(12)

defined on (z, t ) ∈ Ωzn × R+ with BCs

j p0

zi = 0, ∀i ∈ Inj

 j p0 (t ) + l0 (t )r0 (t ) ∂zj e(z, t )   + pj1 (t ) + l0 (t )r1 (t ) e(z, t ) = 0,

In view of the suitable determination of the design parameter µ(z, t ), consider the subsequent lemma concerning the dynamical behavior of (23)–(26).

λi (s)ds/(ρ cp ), i = 1, 2 from (A.1).

pi0 (t )∂zi xˆ (z, t ) + pi1 (t )ˆx(z, t ) = 0,

pi0 (t )∂zi e(z, t ) + pi1 (t )e(z, t ) = 0,



e(z, 0) = e0 (z),

y(zj , t ) = r0 (t )∂zj x(z, t ) + r1 (t )x(z, t ),

0

defined on (z, t ) ∈ Ωzn × R+ with BCs

qi0 (t )∂zi e(z, t ) + qi1 (t )e(z, t ) = 0,

and the system output

 zi

87



µ(z, t )w 2 (z, t )dz  2 w 2 (z, t )dz w (z, t )dz − inf µ(z, t ) Ωzn

z∈Ωzn

  ≤ −2 λ + inf µ(z, t ) V (t ), z∈Ωzn

Ωzn

∀t ∈ R+ 0 .

With (27) this in particular guarantees that zj = 0

(17)

and l(zj , t ), l0 (t ) are the output injection weights to be designed. Introducing the observer error e(z, t ) = x(z, t ) − xˆ (z, t ), the observer error dynamics satisfies

∂t e(z, t ) = 1e(z, t ) + c (z, t )e(z, t )   − l(zj , t ) r0 (t )∂zj e(z0j , t ) + r1 (t )e(z0j , t )

(18)

∥w(z, t )∥L2 ≤ exp(−ϵ t )∥w0 (z)∥L2

(28) 2

providing the exponential stability in the L -norm of the target system (23)–(26). With the specified target system at hand, the next step in the course of the backstepping-based observer design is to evaluate a particular integral transformation relating the desired target dynamics (23)–(26) with the observer error PDE (18)–(22).

88

L. Jadachowski et al. / Automatica 51 (2015) 85–97

3.2. Determination of the kernel-PDE Utilizing the structural property of the reaction parameter c (z, t ) according to Assumption 3, in the following a particular PDE governing the evolution of the output injection weights of the backstepping observer for systems with nD spatial domain is determined. For this, a Volterra-type integral transformation e(z, t ) = w(z, t ) −

zj



g (zj , s, t )w(zsj , t )ds

(29)

0

with an integral kernel g (zj , s, t ) is utilized to relate the exponentially stable w(z, t )-system (23)–(26) with the observer error dynamics (18)–(22). Differentiation of (29) with respect to z yields

∂zi e(z, t ) = ∂zi w(z, t ) −



zj

g (zj , s, t )∂zi w(zsj , t )ds,

i ∈ Inj (30)

0

∂zj e(z, t ) = ∂zj w(z, t ) − g (zj , zj , t )w(z, t )  zj − ∂zj g (zj , s, t )w(zsj , t )ds,

(31)

0

while the Laplacian of (29) results in

1e(z, t ) = 1w(z, t ) − g (zj , zj , t )∂zj w(z, t )   − ∂zj g (zj , zj , t ) + dzj g (zj , zj , t ) w(z, t )   zj g (zj , s, t )∂z2i w(zsj , t )ds − j

i∈In zj

 − 0

Note that (34) has to hold for all (z, t ) ∈ Ωzn × R+ . Moreover, together with Assumption 3 the above equation serves as a motivation to specify the design parameter µ(z, t ) in (23) as

µ(z, t ) = µ0 (t ) − c0 (zj , t )

to compensate the dependency on coordinates zj in the terms (⋆) and (⋆⋆) of (34). With this, the kernel-PDE governing the evolution of the kernel g (zj , s, t ) can be deduced in the form

∂t g (zj , s, t ) = ∂z2j g (zj , s, t ) − ∂s2 g (zj , s, t ) + γ (zj , t )g (zj , s, t ) (36) j

defined on the triangular spatial domain Ωg = {(zj , s) ∈ R2 | s ∈

[0, Lj ], zj ∈ [s, Lj ]} with

γ (zj , t ) = cj (zj , t ) + µ0 (t ).

j

target system is restricted to w( w,j

(32)

 ∂zj g (zj , s, t ) + ∂s g (zj , s, t ) s=z denoting j

j

Lj

(9) and (21), the target system is restricted to ∂zj w(zj , t ) = j w,j w,j −˜qw,j w(zj j , t ), where q˜ w,j = qw, ̸= 0. With (31) 1 /q0 , q0 evaluated at zj = Lj it follows that the kernel has to satisfy L

∂zj g (Lj , s, t ) + q˜ j (t )g (Lj , s, t ) = 0 g (Lj , Lj , t ) = q˜ j (t ) − q˜ w,j j

 zj

 − µ(zsj , t )g (zj , s, t ) w(zsj , t )ds.

(33)

In the next, step substituting the expressions for the observer error (29) as well as its partial derivatives (32), (33) into the equations of the observer error dynamics (18) allows to express all elements in terms of the target system state w(z, t ). After sorting terms in w(z, t ) this leads to 0 = w(z, t ) 2dzj g (zj , zj , t ) − c (z, t ) − µ(z, t )





(⋆)

j

In view of (29) the initial condition at t = 0 for g (zj , s, t ) has to sat-

0



(40)

with q˜ j (t ) = q1 (t )/q0 (t ). Note that in this case either a Neumann or a mixed BC at zj = Lj can be assigned to the target system independent of (21).

 s=z − g (zj , s, t )∂s w(zsj , t ) − ∂s g (zj , s, t )w(zsj , t ) s=0j  zj  − ∂t g (zj , s, t ) + ∂s2 g (zj , s, t )



(39)

(ii) For a mixed (or Neumann) BC at zj = Lj , i.e., q0 (t ) ̸= 0 in

0



, t ) = 0, which implies that

g (Lj , s, t ) = 0.

∂t e(z, t ) = 1w(z, t ) − µ(z, t )w(z, t )   zj g (zj , s, t )∂z2i w(zsj , t )ds − j

Lj zj

q0 = 0 in (25). From (29) evaluated at zj = Lj it follows that the kernel has to satisfy

the total derivative with respect to zj . Differentiation of (29) with respect to t followed by integration by parts taking into account (23) provides

i∈In

γ (zj , t )

. (38) 2 Additional conditions for the kernel-PDE are obtained by evaluating (21) with (29) and (31). However, here it is necessary to distinguish between Dirichlet and a mixed (or Neumann) BC for the PDE (7) at zj = Lj . Moreover, note that the BC at zj = Lj for the target system cannot be specified independently but follows directly from the respective BC for the observer error e(z, t ) and hence from the BC (9) for x(z, t ). As a result, two cases are considered: dzj g (zj , zj , t ) =

(i) For a Dirichlet BC at zj = Lj , i.e., q0 (t ) = 0 in (9) and (21), the

∂z2j g (zj , s, t )w(zsj , t )ds

with dzj g (zj , zj , t ) =

(37)

The respective BC along s = zj follows from (⋆) as

0



(35)

Remark 7. It should be pointed out that the remaining BCs (24) j and (25) of the target system ∀i ∈ In cannot be specified independently but follow directly from the respective BCs (8) and (9) for x(z, t ). Hence, in view of the backstepping transformation w,i (29) it can be shown that Dirichlet BCs (8) and (9) imply p0 = 0 w,i

i = 0, while Neumann BCs (8) and (9) yield pw, = 0 and 1 q1 = 0, respectively. In case of mixed BCs for x(z, t ) it follows w,i w,i w,i w,i that pi1 (t )/pi0 (t ) − p1 /p0 = 0 and qi1 (t )/qi0 (t ) − q1 /q0 = 0 have to be fulfilled for all t ≥ 0. This implies that the resulting BCs



and q0 w,i

zj

 w(zsj , t ) ∂t g (zj , s, t ) + ∂s2 g (zj , s, t ) 0    − ∂z2j g (zj , s, t ) − c (z, t ) + µ(zsj , t ) g (zj , s, t ) ds   

−

(⋆⋆)

  + g (zj , 0, t ) + l(zj , t )r0 (t ) ∂zj w(z0j , t ) + w(z0j , t )   × −∂s g (zj , 0, t ) + l(zj , t ) (r1 (t ) − g (0, 0, t )r0 (t )) .

isfy 0 g (zj , s, 0)w0 (zsj )ds = w0 (z) − e0 (z). Hence, a relationship between w0 (z) and e0 (z) is determined by g (zj , s, 0) = g0 (zj , s) with g0 (zj , s) being consistent with (39) or (40) at t = 0.

(34)

(8) and (9) are basically time invariant. Alternatively, to satisfy the conditions above, (24) and (25) can be assigned as time-varying w,i i w,i i BCs with pl → pw, → qw, l (t ) and ql l (t ), l = 0, 1, which, however, requires a separate consideration of the stability of the corresponding target system.

L. Jadachowski et al. / Automatica 51 (2015) 85–97

Remark 8. Note that the particular choice of the design parameter µ(z, t ) according to (35) allows additionally to redefine the stability condition (27) of the target system derived from the consideration of the rate of change of the Lyapunov functional. Substituting (35) into (27) it can be easily verified that with infz∈Ω n µ(z, t ) = z

j

To determine a solution of (44)–(46), formal integration over Ωg¯ is performed. Integrating (44) with respect to η from 0 to η and subsequently over ξ from η to ξ yields together with (45) g¯ (ξ , η, t ) =

infz∈Ω n µ0 (t ) − c0 (zj , t ) = µ0 (t ) + infz∈Ω n −c0 (zj , t ) = z z µ0 (t ) − supz∈Ω n c0 (zj , t ) it follows that for some ϵ > 0, µ0 (t ) z has to satisfy







µ0 (t ) ≥ sup c0 (zj , t ) − λ + ϵ,

∀t ∈ R+ 0 .



(41)

z∈Ωzn

Consequently, considering the remaining terms in (34) together with (20) and (24), the output injection weights l(zj , t ) and l0 (t ) can be determined depending on the BC for w(z, t ): w,j

(i) For a Dirichlet BC at zj = 0, i.e., p0 one obtains l(zj , t ) = −

g (zj , 0, t ) r0 (t )

j ≡ 0 and pw, ̸= 0 in (24) 1 j

,

l0 (t ) = −

p0 (t ) r0 (t )

.

(42)

Here, it is assumed that r0 (t ) ̸= 0 for all t ≥ 0. w,j (ii) For a mixed (or Neumann) BC at zj = 0, i.e., p0 ̸= 0 in (24) one obtains

(43)

Bg¯ (β, α, t )dα dβ

0

ξ



4 η

 β  γ Lj − , t dβ + g¯ (η, η, t ). 2

(47)

The integration constant g¯ (η, η, t ) is determined according to the different BCs in (46). In case of a Dirichlet BC at ξ = η, i.e., for j q0 (t ) = 0, g¯ (η, η, t ) = 0 directly follows from (46). For a mixed (Neumann) BC, consider that (46) together with the total derivative dη g¯ (η, η, t ) = ∂ξ g¯ (η, η, t ) + ∂η g¯ (η, η, t ) results in dη g¯ (η, η, t ) = fg¯ (η, t ) − q˜ j (t )¯g (η, η, t ),

η 1

(48) η

where fg¯ (η, t ) = 2 0 Bg¯ (η, α, t )dα − 12 γ Lj − 2 , t is obtained by computing ∂ξ g¯ (ξ , η, t ) at ξ = η using (47). Integrating (48) with respect to η from 0 to η together with (46) yields





  1 η β ¯g (η, η, t ) = q˜ j (t ) − q˜ w,j + Bg¯ (β, α, t )dα dβ 2 0 0  η   η  1 β − γ Lj − , t dβ − q˜ j (t ) g¯ (β, β, t )dβ. 2

0

0

Hence, substituting the latter into (47) the implicit solution of (44)– (46) in terms of a Volterra integral equation is given by

w,j

η

 +

1



ξ

η



4 η

Bg¯ (β, α, t )dα dβ

0

Cg¯ (β, t )dβ,

(49)

0

where

    β 1 ξ   γ L − , t dβ, −  j  4 η 2       1 ξ A(ξ , η, t ) = q˜ j (t ) − q˜ w,j − γ Lj −  4 η      1 η   β  − γ Lj − , t dβ, 2

3.3. Solution of the kernel-PDE

0

2

j

if q0 (t ) = 0

β 2

 , t dβ

(50)

j

if q0 (t ) ̸= 0

and

Subsequently, in view of the utilization of the method of integral operators formal integration is applied to the kernel-PDE (36). Therefore, scattering coordinates are introduced ξ = 2Lj − zj − s, η = zj − s, such that with g (zj , s, t ) = g¯ (ξ (zj , s), η(zj , s), t ) the kernel-PDE (36) is transformed into

4

1

η



g¯ (ξ , η, t ) = A(ξ , η, t ) +

Obviously, an explicit solution of the kernel-PDE (36)–(40) is required to determine l(zj , t ) and l0 (t ). However, it can be easily deduced that due to the particular choice of the design parameter µ(z, t ) according to (35), the backstepping-kernel g (zj , s, t ) can be determined similar to the 1D setting, e.g., by means of the so-called method of integral operators (Colton, 1977), which relies on the successive approximation of an implicit integral formulation for g (zj , s, t ).

1

ξ



4 η

−

with p˜ w,j = p1 /p0 assuming that r1 (t ) − r0 (t )(g (0, 0, t ) + p˜ w,j ) ̸= 0 for all t ≥ 0 is satisfied.

∂η ∂ξ g¯ (ξ , η, t ) =

1

2

∂s g (zj , 0, t ) + p˜ w,j g (zj , 0, t )   r1 (t ) − r0 (t ) g (0, 0, t ) + p˜ w,j   j j p1 (t ) − p0 (t ) g (0, 0, t ) + p˜ w,j   l0 (t ) = − r1 (t ) − r0 (t ) g (0, 0, t ) + p˜ w,j

l(zj , t ) =

w,j

89

Bg¯ (ξ , η, t )

Cg¯ (β, t ) =

  0,     1

j

if q0 (t ) = 0 β

Bg¯ (β, α, t )dα  2    j0 −˜q (t )¯g (β, β, t ),

(51) j

if q0 (t ) ̸= 0.

(44)

with Bg¯ (ξ , η, t ) = −∂t g¯ (ξ , η, t ) + γ Lj −

ξ −η

3.4. Convergence analysis

, t g¯ (ξ , η, t ) defined 2 j 2 on the spatial domain Ωg¯ = {(ξ , η) ∈ R | η ∈ [0, Lj ], ξ ∈ [η, 2Lj −η]}. Note that from the scattering transformation it follows that zj = s is equivalent to η = 0 such that the BC (38) is mapped 



to 1  ξ  ∂ξ g¯ (ξ , 0, t ) = − γ Lj − , t .

4 2 Furthermore, (39) and (40) in (ξ , η)-coordinates reads as

 ¯ (η, η, t ) = 0, if qj0 (t ) = 0 g ∂ g¯ (η, η, t ) − ∂ξ g¯ (η, η, t ) + q˜ j (t )¯g (η, η, t ) = 0  η j w,j g¯ (0, 0, t ) = q˜ j (t ) − q˜ w,j , if q0 (t ) ̸= 0, q0 ̸= 0.

(45)

For the explicit evaluation of (49)–(51) successive approximation is applied as suggested in Colton (1977). For this, based on the formal integral formulation of the kernel-PDE and assuming that γ (zj , t ) is sufficiently smooth in t, a solution of (49)–(51) is sought in terms of the infinite series g¯ (ξ , η, t ) =

∞ 

Gn (ξ , η, t ),

(52)

n=0

whose coefficients Gn (ξ , η, t ), n ≥ 0 are determined recursively (46)

ξ η = A(ξ , η, t ) and Gn (ξ , η, t ) = 14 η 0 η BGn−1 (β, α, t )dα dβ + 0 CGn−1 (β, t )dβ, n ≥ 1. Since the solution

from G0 (ξ , η, t )

90

L. Jadachowski et al. / Automatica 51 (2015) 85–97

of (49)–(51) constitutes an infinite series, it is necessary to analyze the convergence of (52). Due to the time-varying parameters γ (zj , t ) and q˜ j (t ) appearing in the recursion formula (52) as well as the appearance of the time derivative of Gn−1 (ξ , η, t ) in the evaluation of BGn (β, α, t ), it is required to analyze the growth rate of the time derivatives of each coefficient Gn (ξ , η, t ), n ∈ N up to an arbitrary order l ∈ N. Subsequently, in order to ensure the convergence of (52), the following assumption is made. Assumption 9. The parameters γ (zj , t ) and q˜ j (t ) satisfy γ (zj , t ) ∈ C 0 ([0, Lj ]) ∩ GMγ ,σ (R+ ) and q˜ j (t ) ∈ GM j ,σ (R+ ), i.e., q˜

  sup ∂tl γ (zj , t ) ≤ M l+1 (l!)σ ,

(53)

  sup ∂tl q˜ j (t ) ≤ M l+1 (l!)σ

+ t ∈R0

with M = max(Mγ , Mq˜ j + q˜ w,j ) and σ ∈ [1, 2].





Note that the imposed restriction on γ (zj , t ) and q˜ j (t ) follow directly from Assumptions 2 and 3. With this, the convergence result can be formulated as follows: Theorem 10 (Convergence of the Successive Approximation). Let γ (zj , t ) and q˜ j (t ) be as in Assumption 9. Then there exists an upper bound of the series coefficients Gn (ξ , η, t ) and their time derivatives in the form

  M l+n+1 sup ∂tl Gn (ξ , η, t ) ≤ n +1

+ t ∈R0

4

n 

×



(l + n)! n!

(n)!(n + 1)!

Qn (ξ , η),

Let us consider the inverse Volterra integral transformation of (29)

w(z, t ) = e(z, t ) +

zj



m(zj , s, t )e(zsj , t )ds

(55)

with the integral kernel m(zj , s, t ), which maps the e(z, t )-system (18)–(22) to the exponentially stable w(z, t )-system (23)–(26). Consequently, proceeding similar to Section 3.2, i.e., by determining 1w(z, t ) in terms of e(z, t ) and m(zj , s, t ) and evaluating ∂t w(z, t ) in view of (18)–(22) followed by a substitution of the resulting expressions into the target system PDE (23) together with (35) and Assumption 3, the kernel-PDE governing the evolution of m(zj , s, t ) can be deduced in the form

∂t m(zj , s, t ) = ∂z2j m(zj , s, t ) − ∂s2 m(zj , s, t ) − γ (s, t )m(zj , s, t )

(54)

for n ∈ N0 , where Qn (ξ , η) = (ξ η)n (ξ − η) for Dirichlet BC, Qn (ξ , η) = (ξ η)n (ξ + η) in case of Neumann BC and Qn (ξ , η) = ηn (ξ + 4)n (ξ + η + 4) for mixed BC. In particular, for σ ∈ [1, 2) j the series (52) converges absolutely and uniformly in the domain Ωg¯ independent of M for all types of BCs. For σ = 2 the series (52) converges absolutely and uniformly if M < 4/L2j in case of Dirichlet or Neumann BC and if M < 4/(Lj (Lj + 4)) in case of mixed BC. The proof of Theorem 10 is provided in Appendix B. By reversing the scattering transformation under Assumptions 3 and 9, this implies that the kernel g (zj , s, t ) is of Gevrey order σ ∈ [1, 2] with respect to t and is a strong solution to the kernel-PDE (36)–(40). With this the observer gains l(zj , t ) and l0 (t ) are determined according to (42) or (43) depending on the corresponding BCs. Remark 11. Although being a useful tool for the verification of the existence and uniqueness of the solution of (49)–(51), the successive approximation method induces a high computational burden due to the recursive determination of the series coefficients in (52). However, the computational costs can be significantly reduced by applying an efficient numerical method based on the trapezoidal quadrature as shown in Jadachowski, Meurer, and Kugi (2012).

(56)

defined on the spatial domain (zj , s) ∈ cording to dzj m(zj , zj , t ) =

σ

(1 + kσ )

k=0

4.1. Inverse backstepping transformation

0

zj ∈ [0, Lj ]

+ t ∈R0

of the observer error dynamics (18)–(22) with the observer gains (42)–(43) is verified based on the analysis of the inverse to (29), i.e., a transformation from the e(z, t )-system to the w(z, t )-system. With this, it can be shown that the stability of the target system implies the stability of the observer error dynamics.

j Ωg

with BCs and the IC ac-

γ (z j , t )

(57)

2

 j m(Lj , s, t ) = 0, if q0 (t ) = 0  w,j ∂z m(Lj , s, t ) + q˜ m(Lj , s, t ) = 0  j j m(Lj , Lj , t ) = q˜ j (t ) − q˜ w,j , if q0 (t ) ̸= 0

(58)

m(zj , s, 0) = m0 (zj , s)

(59)

with m0 (zj , s) being consistent with (58) at t = 0. Hence, observing that (56)–(59) is similar to (36)–(40), both the solution method as well as the convergence results from Theorem 10 can be similarly applied to the inverse kernel m(zj , s, t ). This implies the existence of a bounded strong solution m(zj , s, t ) to the kernel-PDE (56)–(59). 4.2. Stability of the observer error dynamics Recalling that both g (zj , s, t ) and m(zj , s, t ) are bounded strong solutions to the respective kernel-PDEs the following estimates of the integral terms in (29) and (55) are obtained. For the integral term in (29) it follows that

  

zj 0

2 

g (zj , s, t )w(zsj , t )ds

≤ Lj

sup j

(zj ,s)∈Ωg

L2

g 2 (zj , s, t )∥w(z, t )∥2L2 ,

∀t ∈ R + 0

(60)

and for the integral term in (55) with t = 0 we get

  

zj 0

2 

m(zj , s, 0)e(zsj , 0)ds

≤ Lj

sup j (zj ,s)∈Ωg

L2

m (zj , s, 0)∥e0 (z)∥2L2 . 2

(61)

4. Stability and robustness analysis

Additionally, observe that (29) together with the Minkowski inequality and the estimate (60) implies

In Section 3, the observer gains l(zj , t ) and l0 (t ) are designed to stabilize the observer error dynamics (18)–(22). This is obtained by making use of the invertible Volterra transformation (29) mapping the exponentially stable target system (23)–(26) into the observer error dynamics. In the following, the exponential stability

∥e(z, t )∥L2 ≤ Ce (t )∥w(z, t )∥L2 , (62)  where Ce (t ) = 1 + Lj sup(z ,s)∈Ω j g 2 (zj , s, t ). On the other hand, j

g

an upper bound for ∥w0 (z)∥L2 can be determined by taking into account the boundedness of the inverse kernel m(zj , s, t ) together

L. Jadachowski et al. / Automatica 51 (2015) 85–97

with (61), i.e.,

with

∥w0 (z)∥L2 ≤ Cw ∥e0 (z)∥L2 , (63)  where Cw = 1 + Lj sup(z ,s)∈Ω j m2 (zj , s, 0). Hence, substituting j

K (t ) =



Lj sup c˜ 2 (z, t )

∥e(z, t )∥L2 ≤ Cw Ce exp(−2ϵ t )∥e0 (z)∥L2 (64) with Ce = supt ∈R+ Ce (t ), which ensures that the observer error 0 e(z, t ) of the state observer with observer gains l(zj , t ) and l0 (t ) according to (42), (43) decays exponentially in the time t in the L2 -norm.

Proposition 13. The parabolic PDE (65) with BCs (24), (25) and IC (26) is exponentially stable in the L2 -norm if the inequality

µ0 (t ) ≥ sup c (z, t ) + K (t ) − inf

zj



g (zj , s, t )w( , t )ds

V˙ (t ) =

 Ωzn

Lemma 12. Let g (zj , s, t ) be bounded, strong solutions of the kernelPDE (36)–(40). Moreover, let c˜ (z, t ) be a bounded function on (z, t ) ∈ + Ωzn × R+ 0 . Then there exists a K (t ) ∈ R0 such that the following estimate holds

w(z, t )˜c (z, t )

g (zj , s, t )w( , t )dsdz zsj

0

≤ K (t ) ∥w(z, t )∥2L2 ,

∀t ∈ R+ 0 .

(66)

Ωzn

w(z, t )˜c (z, t )

g (zj , s, t )w(zsj , t )dsdz

0



 21

w (z, t )˜c (z, t )dz 2

≤ Ωzn



2



Lj

× Ωzn

 12

g (zj , s, t )w ( , t )dsdz 2

2

zsj

0

 ≤

zj



sup c˜ (z, t ) 2

z∈Ωzn

 12



w (z, t )dz 2

Ωzn

 ×  sup g 2 (zj , s, t ) j (zj ,s)∈Ωg

≤ K (t ) ∥w(z, t )∥2L2

 + Ωzn

w(z, t )˜c (z, t )

zj



 Ωzn

 Ωzn

w 2 (z, t )dz

cj (zj , t )w 2 (z, t )dz

g (zj , s, t )w(zsj , t )dsdz

0

≤ −λ ∥w(z, t )∥2L2 − µ0 (t ) ∥w(z, t )∥2L2 + sup c (z, t ) ∥w(z, t )∥2L2 + K (t ) ∥w(z, t )∥2L2 z∈Ωzn

−

inf cj (zj , t ) ∥w(z, t )∥2L2

zj ∈[0,Lj ]

≤ −2θ (t )V (t ), where θ (t ) = λ+µ0 (t )− supz∈Ω n c (z, t )− K (t )+ infzj ∈[0,Lj ] cj (zj , t ). z For a suitable choice of the design parameter µ0 (t ) satisfying (67) it follows that

∥w(z, t )∥L2 ≤ exp(−ϵ t )∥w0 (z)∥L2

(68) 2

Proof. The proof of Lemma 12 directly follows from the consideration of the boundedness of g (zj , s, t ) and c˜ (z, t ). Consequently, the Cauchy–Schwarz inequality yields



c (z, t )w (z, t )dz − 2

Ωzn

the stability of the resulting target system (65) with BCs (24), (25) and IC (26) the boundedness of the arising integral term  zj c˜ (z, t ) 0 g (zj , s, t )w(zsj , t )ds is proven. For this consider the following lemma.

Ωzn

1w(z, t )w(z, t )dz − µ0 (t )

+

  − µ0 (t ) − c (z, t ) + cj (zj , t ) w(z, t ) (65) with c˜ (z, t ) = c0 (zj , t ) + cj (zj , t ) − c (z, t ). In order to analyze

(67)

Proof. Proceeding similar to Section 3.1 for the proof of Lemma 6, the proposition can be verified by means of the consideration of the change of the positive definite Lyapunov functional V (t ) = 1 ∥w(z, t )∥2L2 . Computing the rate of change of V (t ) along a solution 2 trajectory w(z, t ) together with Lemma 12 yields



zsj

cj (zj , t ) − λ + ϵ

is satisfied for all t ∈ R+ 0 and some ϵ > 0, where λ denotes the smallest eigenvalue of the Sturm–Liouville problem 1w(z, t ) + λw(z, t ) = 0 with BCs (24) and (25).

0

zj

zj ∈[0,Lj ]

z∈Ωzn

∂t w(z, t ) = 1w(z, t ) + c˜ (z, t )



(zj ,s)∈Ωg

With these preparations, exponential stability of the target system (65) can be deduced also for a non-separable reaction parameter c (z, t ).

4.3. Robustness analysis The validity of the above stability analysis relies on Assumption 3, i.e., the spatial separability of the reaction parameter c (z, t ) and the subsequent determination of the observer gains l(zj , t ) and l0 (t ) as outlined in Section 3. In the following, the robustness of the determined state observer (12)–(16) with l(zj , t ) and l0 (t ) according to (42) and (43), respectively, is investigated in view of a non-separable c (z, t ) in (18). For this, in a first step the dynamics of the target system is determined for the case that (29) is applied to (18)–(22) and Assumption 3 does not hold. Substitution of (29) and its partial derivatives into (18) in view of (36)–(38) and (42), (43) yields

g 2 (zj , s, t ). 

sup j

z∈Ωzn

g

the stability result for the target system (28) (see Lemma 6) into (62) together with (63) yields



91

 Ωzn

Lj

 0

 21 w2 (zsj , t )dsdz

providing the exponential stability in the L -norm of the target system (65) with BCs (24), (25) and IC (26).  The stability result for the target system (65) allows to verify the stability of the observer error dynamics in the case of a nonseparable parameter c (z, t ). Substituting (68) into the observer error estimate (62) yields ∥e(z, t )∥L2 ≤ Cw Ce (t ) exp(−ϵ t )∥e0 (z)∥L2 , which guarantees the exponential decay of the norm ∥e(z, t )∥L2 for all (z, t ) ∈ Ωzn × R+ 0 in the case of non-separable c (z, t ) assuming (66). Remark 14. Note that the determination of the design parameter µ0 (t ) satisfying (67) is a non-trivial task. Since the backstepping kernel g (zj , s, t ) is obtained in terms of the design parameter µ0 (t ) according to (41), an implicit inequality has to be satisfied to ensure the exponential stability of the norm ∥e(z, t )∥L2 . This can be achieved by means of an iterative determination of the backstepping kernel g (zj , s, t ) subject to the design parameter µ0 (t ) followed by consecutive evaluation of (67). At present, the convergence of this iterative method is not ensured, which remains an open problem. However, the investigation of the existence of the iterative solution to (67) is beyond the scope of this paper.

92

L. Jadachowski et al. / Automatica 51 (2015) 85–97

(5) from a limited set of measurements yk (t ), k = 1, . . . , Nk by means of a least squares approximation.

5. Approximate realization of the idealized system output As already mentioned in Remark 1, the considered output equation (5) postulates an idealized setting with an infinitedimensional measurement being restricted to the boundary hyperplane ∂ Ωj = {z ∈ Ωzn | zj = 0} for a single but arbitrary j ∈ In . Hence, an exact realization of the distributed-parameter observer (12)–(16) requires the complete knowledge of the state evolution of the system on ∂ Ωj . This requirement is usually not feasible, since it would in principle presume an infinite number of sensors on the considered hyperplane. Therefore, in the following an approach is presented, which allows an approximate realization of the idealized output by means of a finite number of appropriately located sensors. 5.1. Optimal sensor location An important factor for an accurate state estimation in DPSs is the selection of an appropriate sensor configuration, i.e., number and location of the sensors in the spatial domain under consideration. In the past years, different approaches for the optimal sensor placement have been published, see, e.g., Armaou and Demetriou (2006), Fahroo and Demetriou (2000), Vande Wouwer, Point, Porteman, and Remy (2000) and the references therein. A common approach is to define an optimality criterion, which is maximized with respect to the sensor locations. Subsequently, in line with Vande Wouwer et al. (2000), the optimality criterion is based on the knowledge of the dynamic behavior of the system model, which can be investigated in simulations provided that realistic scenarios are considered. In particular, the optimality criterion relies on a measure of independence between elements of a finite-dimensional vector of sensor signals y(t ) = [yk (t )]k=1,...,Nk defined in terms of weighted integrals yk (t ) =



sk (zj ) r0 (t )∂zj x(z, t ) + r1 (t )x(z, t ) dzj ,



∂ Ωj



(69)

zj = 0, k = 1, . . . , Nk with the so-called sensor shape functions sk (zj ), k = 1, . . . , Nk of a given geometry having compact support in ∂ Ωj . Thereby, the measure of independence between sensor signals is the Gramian determinant det (F) with a positive semi-definite matrix F = [Fk,l ] formally given by T



yk (t )yl (t )dt ,

Fk,l =

k, l = 1, . . . , Nk ,

(70)

0

where [0, T ] is the time span of interest allowing to capture the essential dynamics of the system. Here, the diagonal of (70) consists of the L2 -norm of each yk (t ) and represents the energy levels of the sensor signals. The off-diagonal elements of F describe a degree of dependence between two sensor signals. Obviously, for linearly independent yk (t ), k = 1, . . . , Nk the determinant of (70) is always positive and increases with its diagonal elements. For linearly dependent sensor signals it follows that det (F) = 0. Hence, the solution of the maximization problem max det (F)

(71)

Zj

Remark 15. It is evident that the quality of the output reconstruction also depends on the selection of the number of sensors Nk as well as their shape functions sk (zj ), k = 1, . . . , Nk and in general increases with Nk . However, in order to take into account all these aspects, a suitable modification of the considered optimality criterion with respect to the number of sensors Nk and their spatial geometry sk (zj ) is required. Remark 16. Note that the time span [0, T ], for which the optimality criterion is evaluated, has to be chosen with regard to the time evolution of the varying boundary parameters pi0 (t ), pi1 (t ), qi0 (t ), qi1 (t ), i ∈ In and the reaction parameter c (z, t ). If the time interval is too short it may happen that the determined sensor configuration does not capture essential information of the system output such that the quality of the reconstructed output will not be sufficient to assure an accurate state estimation. 5.2. Output realization In practice, the measurements are mostly obtained at discrete times tk , k = 0, 1, . . . such that (69) becomes yk (tk ) =

with respect to a vector of sensor locations Zj = [zj ]k=1,...,Nk allows to determine spatial regions, where yk (t ) is maximized and the interference with other measurements yl (t ), l ̸= k is mini(k) (k) (k) (k) (k) mized. Here, zj = (z1 , . . . , zj−1 , zj+1 , . . . , zn ), k = 1, . . . , Nk denotes the geometric center of the kth sensor shape function on the hyperplane ∂ Ωj . In the following subsection, given Nk sensors with sensor shape functions sk (zj ), the optimal sensor configuration Zj is utilized to reconstruct the idealized system output

sk (zj ) r0 (tk )∂zj x(z, tk ) + r1 (tk )x(z, tk ) dzj ,



∂ Ωj



(72)

zj = 0, k = 1, . . . , Nk . This motivates a sample-and-hold implementation of the state observer (12)–(16) with a sampling time Ta such that an update of the output information y(tk ) is available only at tk = {0, Ta , 2Ta , . . .} and held constant for t ∈ [tk , tk+1 ) in the observer evaluation. Henceforth, the idealized system output (11) will be reconstructed by means of a series expansion y∗ (zj , tk ) =

Nl 

al (tk )φl (zj )

(73)

l=1

in terms of linear independent spatial weighting functions φl (zj ), l = 1, . . . , Nl , which are consistent with the system’s BCs, i.e., pi0 (tk )∂zi φl (zj ) + pi1 (tk )φl (zj ) = 0,

zi = 0,

(74)

qi0 (tk )∂zi φl (zj ) + qi1 (tk )φl (zj ) = 0,

zi = Li ,

(75)

j In

for all i ∈ at each sampling time tk . The fitting parameters a(tk ) = [al (tk )]l=1,...,Nl are thereby obtained in every sampling period [tk , tk+1 ) by substituting a vector of approximations y∗ (tk ) = [y∗k (tk )]k=1,...,Nk defined elementwise by y∗k (tk ) =

 ∂ Ωj

sk (zj )y∗ (zj , tk )dzj

(76)

for k = 1, . . . , Nk into the discrete minimization problem

2

min y(tk ) − y∗ (tk )2 .



a(tk )

(k)



(77)

Hence, with a(tk ) determined from (77) a time-discrete realization y∗ (zj , tk ) of the idealized system output y(zj , t ) is achieved in a least squares sense. Assumption 17. Under practical conditions, it is often assumed that the state x(z0j , t ) and its partial derivative ∂zj x(z0j , t ) can be measured on ∂ Ωj only at discrete pointwise sensor positions (k)

zj , k = 1, . . . , Nk . This implies that the sensor shape functions

L. Jadachowski et al. / Automatica 51 (2015) 85–97

take the form of (n − 1)-dimensional Dirac delta functions, i.e., (k) sk (zj ) = δ(zj − zj ) for k = 1, . . . , Nk . With this, (72) results 0,(k)

in a set of pointwise measurements yk (tk ) = r0 (tk )∂zj x(zj 0,(k) zj

0,(k) zj

(k)

, tk ) +

(k) (k) (z1 , . . . , zj(k) −1 , 0, zj+1 , . . . , zn )

r1 (tk )x( , tk ) with = vector of coordinates of the kth sensor on ∂ Ωj .

a

Linear interpolation functions (see, e.g., Bitzer & Zeitz, 2002) are a common choice of spatial weighting functions. In this case, however, the fulfillment of the boundary conditions (74), (75) has to be considered separately. On the other hand, if additional knowledge about the evolution of the system state on ∂ Ωj is available, more sophisticated weighting functions, e.g., higher order interpolation functions, can be used as suggested in Vande Wouwer and Zeitz (2001). Thereby, special attention has to be paid when dealing with convection dominated systems. In the following, the spatial weighting functions φn (zj ) are chosen as the eigenfunctions satisfying the eigenvalue problem n 

¯ zj ), ∂z2i φ(zj ) + c0 (zj , 0)φ(zj ) = λφ(

i ∈ Inj

(78)

i=1

with BCs (74) and (75). This choice of the weighting functions guarantees that (74) and (75) are satisfied automatically. Remark 18. Since the proposed observer design as well as the stability analysis from Section 4.2 are based on the availability of the idealized system output y(zj , t ), it has to be pointed out that the exponential convergence of the observer error ∥e(z, t )∥L2 cannot be assured for the application of the sampled output approximation y∗ (zj , tk ). From a theoretical point of view, this is still an open problem addressing the impact of the output approximation on the convergence of the observer error. However, the simulation results presented in Section 6 confirm the feasibility of the proposed approach.

Table 1 Values of the fitting parameters al (t ), l = 1, 2, 3 corresponding to the output reconstruction with 4 optimally located sensors and 3 eigenfunctions at selected time points t ∈ {0, 0.4, 1.2, 2.0}. Snapshot time

Fitting parameters

t

a 1 (t )

a2 (t )

a3 (t )

0.43 1.2 1.87 3.7

−0.4 −0.09

−0.1

t t t t

=0 = 0.4 = 1.2 = 2.0

Motivated by the results above, the state estimation problem is considered for (7)–(11) defined on a 3D rectangular spatial domain. The spatially and time-varying parameter is given in view of Assumption 3 by

 1  c (z1 , z2 , z3 , t ) = (z3 + 1)(cos(π t ) + 1) exp − t 3 + (1 + 2 cos(2π z1 ))(1 + 2 cos(π z2 )) cos(2π t )   

(79)

c0 (z1 ,z2 ,t )

and the time-varying boundary parameters at z3 ∈ {0, L3 } read as 1 (1 + sin(2π t ))Φ2,1 (t ), 2 q31 (t ) = 2 + 2(1 + sin(4π t ))Φ1,1 (t ). p31 (t ) =

Φω,T (t ) =

 0,    t   

0

1,

t ≤0

φω,T (τ )dτ

T



(k)

φω,T (τ )dτ ,

t ∈ (0, T )

0

t ≥T ω

with φω,T (t ) = exp(−1/[(1 − t /T )t /T ] ) if t ∈ (0, T ) and φω,T (t ) = 0 if t ̸∈ (0, T ). Hence Φω,T (t ) is non-analytic at t ∈ {0, T } and of Gevrey order σ = 1 + 1/ω, ω > 0, see Definition 19 in Appendix C. The remaining system parameters are chosen as pi0 (t ) = 1, qi0 (t ) = 1 with i ∈ I3 , pi1 (t ) = 0, qi1 (t ) = 0 with i ∈ I33 and r0 (t ) = 0, r1 (t ) = 1. Clearly, this provides Neumann BCs at

(k)

Sensor number

Sensor position (z1 , z2 )

k

Initial

Optimal

k=1 k=2 k=3 k=4

(0.2, 0.2) (0.25, 0.6) (0.6, 0.25) (0.9, 0.6)

(0.236, 0.236) (0.25, 0.748) (0.748, 0.236) (0.752, 0.748)

zi ∈ {0, Li } for i ∈ I33 , time-varying mixed BCs at z3 ∈ {0, L3 } and a system output y(z1 , z2 , t ) = x(z1 , z2 , 0, t ) restricted to the output surface ∂ Ωj := {z ∈ Ωz3 | z3 = 0}. The IC of the plant is assigned as x0 (z) = 12 (1 − cos(2π z1 ))(1 − cos(π z2 )). Moreover, note that in view of practical feasibility of the designed state estimator only a finite number Nk of sensors are presumed on ∂ Ωj . Hence, the presented solution is based on a 2-step approach outlined below. Determination of the observer gains: In a first step, the backstepping-based observer (12)–(17) with a zero IC xˆ (z, 0) ≡ 0 is determined as outlined in Section 3. Here, assuming the availability of the idealized system output y(z1 , z2 , t ), the target system (23)–(26) is specified according to (35), i.e., µ(z1 , z2 , z3 , t ) = µ0 (t ) − c0 (z1 , z2 , t ) with a design parameter i i µ0 (t ) ≡ µ0 ∈ {1, 5, 10} and boundary parameters pw, = qw, =1 0 0 and p1

i = qw, = 0, i ∈ I3 . 1

Idealized output reconstruction: In a second step, the idealized system output is reconstructed by means of the least squares approach as introduced in Section 5.2 using Nl = 3 eigenfunctions as spatial weighting functions in (73). Here, the eigenfunctions φl (z1 , z2 ), l = 1, . . . , Nl are determined as solutions of the corresponding eigenvalue problem (78). Thereby, the sensor configuration is specified as outlined in Section 5.1 by solving the optimization problem (71) with respect to Nk = 4 sensors with (k) (k) (k) z3 = (z1 , z2 ) the position of the kth sensor on ∂ Ωj . Results for the determination of the fitting parameters al (t ), l = 1, 2, 3 at selected time points t ∈ {0.0, 0.4, 1.2, 2.0} are exemplarily presented in Table 1, while results of the optimal sensor location given an initial configuration are summarized in Table 2. Note that, in order to avoid location redundancy, upper and lower coordinate bounds for each spatial direction z1 and z2 are assigned as 0.05 ≤ (k) z1 ≤ 0.45 for k = 1, 2, 0.55 ≤ z1(k) ≤ 0.95 for k = (k)

Here, Φω,T (t ) denotes a Gevrey-class function given by

0.02

−0.1 −0.11

0.21 0.59

Table 2 Optimal sensor coordinates.

w,i

6. Simulation example for parabolic PDEs on 3D spatial domains

93

(k)

3, 4, 0.05 ≤ z2 ≤ 0.45 for k = 1, 3 and 0.55 ≤ z2 ≤ 0.95 for k = 2, 4. Moreover, it is assumed that the state x(z1 , z2 , 0, t ) can only be measured at the discrete pointwise sensor positions, which implies that the sensor shape functions s(z1 , z2 ) take the form of 2-dimensional Dirac delta functions, i.e., sk (z1 , z2 ) = δ(z1 − (k) (k) z1 , z2 − z2 ), k = 1, . . . , Nk . Additionally, a sample-and-hold implementation is considered with the sampling time Ta = 0.02. Thus, a new output information y∗ (z1 , z2 , tk ) is only available at tk = {0, Ta , 2Ta , . . .} and held constant for t ∈ [tk , tk+1 ) in the observer evaluation. Numerical results: For the numerical evaluation of the considered state estimation problem, the governing equation of the plant (7)– (11) is discretized in space and time using a 3-dimensional implicit and numerically absolutely stable Crank–Nicholson scheme.

94

L. Jadachowski et al. / Automatica 51 (2015) 85–97

a

b

c

d

e

f

Fig. 1. Numerical results for the state estimation problem. (a) Solution x(0.5, 0.5, z3 , t ) in the (z3 , t )-domain; (b) spatial–temporal profile of l(z3 , t ) for µ0 = 5; (c) timeevolution of l(z3 , t ) at z3 = 0.5 (gray) and l0 (t ) (black) for µ0 = 5; (d) observer error e(0.5, 0.5, z3 , t ); (e) L2 -norm of the simulator error (dash-dotted) and of the observer error with µ0 = 1 (solid), µ0 = 5 (dashed) and µ0 = 10 (gray); (f) L2 -norm of the observer error for µ0 = 5 with 1zˆ = 0.1 (gray) and 1zˆ = 0.2 (dashed) based on the idealized output as well as based on the approximated output with 1zˆ = 0.2 and Ta = 0.02 (solid).

Thereby, each spatial dimension is divided into 10 equidistant intervals of length 1zi = 1z = 0.1, i = 1, 2, 3. An analogous discretization is considered for the implementation of the state observer (12)–(17), where in view of the impact of the discretization on the observer error convergence the observer domain is discretized with a coarser grid, i.e., 1zˆi = 1zˆ ≥ 1z, i = 1, 2, 3. Simulation results for the system (7)–(11) in the (z3 , t )-domain at z1 = z2 = 0.5 are exemplarily depicted in Fig. 1(a). Here, note the unstable plant dynamics due to the particular choice of the system parameters. Moreover, consider that in view of the zero IC xˆ (z, 0) ≡ 0 of a simple simulator (l(z3 , t ) = l0 (t ) = 0) the resulting simulator error dynamics in the respective spatial domain is identical to Fig. 1(a). Numerical results for the determination of the observer gains l(z3 , t ) and l0 (t ) are shown in Fig. 1(b), where the spatial–temporal profile of l(z3 , t ) is depicted for µ0 = 5, and in Fig. 1(c), where the time-evolution of both observer gains l(z3 , t ) for z3 = 0.5 and l0 (t ) is studied. With this, Fig. 1(d) shows the corresponding exponentially stable observer error dynamics exemplarily plotted in the (z3 , t )-domain at z1 = z2 = 0.5. Furthermore, the time-evolution of the L2 norm of the observer error ∥e(z, t )∥L2 is depicted in Fig. 1(e) for different µ0 ∈ {1, 5, 10} and compared to the L2 -norm of the simple simulator. Here, the numerical results are obtained based on the availability of the idealized system output y(z1 , z2 , t ) and an observer/simulator discretization 1zˆ = 1z = 0.1. Obviously, due to theoretical results obtained in Section 4.2, by increasing µ0 the decay rate of the observer error can be directly improved providing an accurate state estimation. In view of the impact of the observer discretization and the output reconstruction on the observer error decay, consider Fig. 1(f), where ∥e(z, t )∥L2 of the backstepping observer implemented with 1zˆ = 0.2 is illustrated for µ0 = 5 relying on both idealized and approximated output. Note that the coarsely discretized observer equations marginally affect the accuracy of the state estimation. Moreover, the effect of the coarse discretization is amplified by the observer implementation based on the approximated and sampled output. Nevertheless, the results confirm a very accurate estimation performance.

with measured output being restricted to a single boundary hyperplane. For this, a backstepping transformation is applied to determine the kernel-PDE governing the evolution of the observer gains. Assuming the availability of an idealized system output, the computed observer gains ensure that the observer error is exponentially stable in the L2 -norm as it is theoretically proven. Furthermore, in view of the practical feasibility of the designed state estimator, the restriction to the idealized output is weakened and the state observer is implemented with the system output realized by means of a least squares approach under the assumption that only a finite number of optimally located measurements is available on the output hyperplane. The performance of the proposed observer design is illustrated in simulation scenarios, which confirm an accurate state estimation.

7. Conclusions

∂t x¯ (ζ, t ) = 1x¯ (ζ, t ) + c¯ (ζ, t )¯x(ζ, t )

In this paper, a Luenberger-type state observer is presented for linear parabolic PDEs defined on an nD cuboidal spatial domain

defined on (ζ, t ) ∈ Ωζ × R , Ωζ := {ζ ∈ R | 0 < ζi < L¯ i , i = n 2 1, . . . , n} with the Laplace operator ∆ = i=1 ∂ζi , corresponding

Appendix A. Transformation into DRS (7)–(11) For the solution of the observer design problem, the DCRS (1)–(4) is transformed into a simpler form by applying a suitable coordinate transformation followed by the application of the Hopf–Cole state transformation, see, e.g., Joseph and Vasudeva Murthy (2001). Consequently, the parameters ai (zi ) are normalized to 1 by making use of the coordinate change zi → ζi = z¯i (zi ) :=

zi



− 12

ai

(s)ds,

∀i ∈ In

(A.1)

0

with Li → L¯ i := z¯i (Li ) and z¯ (z) = ζ = (ζ1 , . . . , ζn ), while the parameters bi (zi ) are eliminated taking into account the Hopf–Cole transformation n ζ  i

x(z, t ) → x¯ (ζ, t ) := x(ζ, t )e

i=1

0

b¯ i (s)ds

(A.2)

  with b¯ i (ζi ) = 12 ai (zi )∂z2i z¯i (zi ) + bi (zi )∂zi z¯i (zi ) evaluated for zi = z¯i−1 (ζi ). As a result, instead of studying the PDE (1)–(5) it is equivalent to consider the DRS

n

+

n

(A.3) n

L. Jadachowski et al. / Automatica 51 (2015) 85–97

BCs

¯ (t )∂ζi x¯ (ζ, t ) + ¯ (t )¯x(ζ, t ) = 0,

ζi = 0

(A.4)

¯ (t )∂ζi x¯ (ζ, t ) + ¯ (t )¯x(ζ, t ) = 0,

ζi = L¯ i

(A.5)

pi1

pi0

qi1

qi0

for all i ∈ In , t > 0, consistent IC x¯ (ζ, 0) = x¯ 0 (ζ),

ζ ∈ Ωζn

(A.6)

and the system output y¯ (ζ j , t ) = r¯0 (t )∂ζj x¯ (ζ, t ) + r¯1 (t )¯x(ζ, t ),

ζj = 0 (A.7) for t ≥ 0 with ζ j = ζ \ {ζj }, such that y(zj , t ) → y¯ (ζ j , t ) := ζ    y(ζ j , t ) exp − i∈Ij 0 i b¯ i (s)ds . Thereby, the reaction parameter n in (A.3) is mapped according to c (z, t ) → c¯ (ζ, t ) := c (ζ, t ) −  n  2 ¯ ¯ b (ζ ) + ∂ b (ζ ) , which in view of Assumption 3 results i ζi i i i i =1 in c¯ (ζ, t ) = c¯0 (ζ j , t ) + c¯j (ζj , t ). Consequently, the boundary

parameters follow from pi0 (t ) → p¯ i0 (t ) := pi0 (t ), pi1 (t ) → √ p¯ i1 (t ) := ai (0)pi1 (t ) − b¯ i (0)pi0 (t ), qi0 (t ) → q¯ i0 (t ) := qi0 (t ), √ i i q1 (t ) → q¯ 1 (t ) := ai (Li )qi1 (t ) − b¯ i (L¯ i )qi0 (t ), while the output  parameters are given by r0 (t ) → r¯0 (t ) := r0 (t )/ aj (0), r1 (t ) → r¯1 (t ) := r1 (t ) − r0 (t )b¯ j (0)/ aj (0) for i ∈ In and t > 0. Moreover, by reversing the transformation (A.1) and (A.2), it follows that condition (6) is preserved in the transformed j parameters p¯ l (t ), r¯l (t ), l = 0, 1.

95

Neumann BC: In the case of Neumann BC (q˜ j (t ) = 0 in (50) and (51)), it follows for G0 (ξ , η, t ) that sup ∂tl G0 (ξ , η, t )





+ t ∈R0

≤

1

ξ

  ∂ l γ Lj − β , t dβ + 1 t



sup ∂tl GN +1 (ξ , η, t )





+

t ∈R0

≤

1

ξ



η



4 η

 l  ∂ BG (β, α, t ) dα dβ + t N

0

4

k=0

  1 + T2 ◦ ∂tl+1 GN  + 2

≤

k

l   

l

k

k=0

(53),(54)

M

l +N +2



4N +2

+

0  ξ l   β 1 l +1   (l!)σ (ξ − η), which is equal to ≤ η ∂t γ Lj − 2 , t dβ ≤ 4 M (54) for n = 0. In the inductive step assume that (54) holds for all n = 0, 1, . . . , N , N > 1. Then it follows for n = N + 1 from the expression for Gn (ξ , η, t ) together with the general Leibniz rule

M l +N +2

1 4

+2

and Lemma 20 that1 + t ∈R0

1



ξ

η



4 η

 l  ∂ BG (β, α, t ) dα dβ t N

4

k=0

(53),(54)

≤

+

M

l+N +2

4N +2 M l +N +2 4 N +2





k

(l + N + 1)! N!

(l + N + 1)! (N + 1)!

σ

σ

N 

(1 + k σ )

k=0

(N !)(N + 1)! N  (1 + kσ ) k=0

(N !)(N + 1)! N +1 (1 + kσ )

T1 ◦ QN

T1 ◦ QN

 (l + N + 1)! σ k=0 ≤ N +2 QN +1 , 4 (N + 1)! (N + 1)!(N + 2)! which is equal to (54) for n = N + 1. M l+N +2



1 For the sake of readability, here and in the following G = GN (β, α, t ), N   γ = γ Lj − β−α , t and Q = Q (β, α) are used. Moreover the abbreviations N N 2 T1 ◦ G N = introduced.

ξ η η

0

+2

0

  l    l +1   l  l −k   k  1       ∂t γ ∂t GN ≤ T1 ◦ ∂t GN +

GN (β, α, t )dα dβ and T2 ◦ GN =

ηβ 0

0

GN (β, α, t )dα dβ are

 l  ∂ CG (β, t ) dβ t N

0

(l + N + 1)! N!

σ

j



η



l      l  l−k   k  1 ∂ γ  ∂ GN  ≤ T1 ◦ ∂tl+1 GN  + t t

Dirichlet BC: For Dirichlet BC (q0 (t ) = 0) with γ (zj , t ) and q˜ j (t ) satisfying (53), it follows for G0 (ξ , η, t ) that supt ∈R+ ∂tl G0 (ξ , η, t )



2

which corresponds to the bound for ∂tl G0 (ξ , η, t ) obtained from (54). Hence, in the inductive step it follows for n = N + 1 from the equation for Gn (ξ , η, t ) together with the Leibniz rule and Lemma 20 that

Appendix B. Proof of Theorem 10

sup ∂tl GN +1 (ξ , η, t ) ≤

  ∂ l γ Lj − β , t dβ t

4 η 2 2 0 1 1 l+1 ≤ M (l!)σ (ξ − η) + M l+1 (l!)σ η 4 2 M l +1 σ ≤ (l!) (ξ + η), 4



In order to verify condition (54), in a first step induction is applied. For this, it is necessary to distinguish between Dirichlet, Neumann and a mixed BC in the case n = 0 and in the induction hypothesis (n = N + 1).

η





4N +2 M l +N +2

(l + N + 1)! (N + 1)!



4N +2 M l +N +2 4N +2



σ

(l + N + 1)! N! (l + N + 1)! (N + 1)!

 l−k   k  ∂ γ  ∂ GN  t t N 

(1 + kσ )

k=0

(N !)(N + 1)! N  (1 + kσ ) k=0

T1 ◦ QN

T1 ◦ QN

(N !)(N + 1)! N  (1 + kσ ) σ k=0

(N !)(N + 1)! N  (1 + kσ ) σ

T2 ◦ QN

(N !)(N + 1)! N +1 (1 + k σ )

T2 ◦ QN

k=0

 (l + N + 1)! σ k=0 ≤ N +2 QN +1 . 4 (N + 1)! (N + 1)!(N + 2)! This equals (54) with n = N + 1 for the respective BC. Mixed BC: In case of mixed BC (q˜ j (t ) ̸= 0 in (50) and (51)), it follows for n = 0 that      l  1 ξ l    ∂ γ Lj − β , t dβ sup ∂t G0 (ξ , η, t ) ≤ t   + 4 2 M l +N +2



η

t ∈R0

       1 η l ∂ γ Lj − β , t dβ + ∂ l q˜ j (t ) + t t   2 2 0

1 M (l!)σ (ξ − η) + M l+1 (l!)σ η + M l+1 (l!)σ 4 2 M l +1 (l!)σ (ξ + η + 4), ≤ 4 which is equal to (54) with Q0 (ξ , η) = (ξ + η + 4) for the mixed BC. For n = N + 1 it follows for the bound of Gn (ξ , η, t ) that

≤

1

l+1

sup ∂tl GN +1 (ξ , η, t )



+

t ∈R0



96

L. Jadachowski et al. / Automatica 51 (2015) 85–97

≤ ≤



1

ξ

4 η 1 4

η



 l  ∂ BG (β, α, t ) dα dβ + t N

0

η



  l ∂ CG (β, t ) dβ t N

0



T1 ◦ ∂tl+1 GN  +



l   

l

k

k=0

  1 + T2 ◦ ∂tl+1 GN  + 2

 l−k   k  ∂ γ  ∂ GN  t t

l   

l

k=0

For the time-varying system parameters, the class of Gevrey functions GD,σ (Ω ) plays an important role.

t

 l     l  l−k j  η  k   ∂ GN (β, β, t )dβ. + ∂t q˜ (t ) t k=0

k

0

In the next step, the last term in the upper expression is multiplied 1+(N +1)σ ≥ 1, N ∈ N0 , σ ∈ [1, 2] in order to be able by a(N ) = N +2 to factorize the resulting expression. Utilizing the results for the Neumann BC this yields



+ t ∈R0

N +1

≤

M

l +N +2



(l + N + 1)! (N + 1)!

σ



(1 + kσ )

k=0

η

0

M l +N +2

0

=−

N +1

≤

In the course of the convergence analysis in Appendix B the following identities are utilized.

l    σ  (l + N + 1)! σ l  (l − k)!(k + N )! ≤ N +1 k k=0  ξ η N +1 (ξ η) (ξ − η) (βα)N (β − α)dQ = ( N + 1)(N + 2) η 0  ξ η (ξ η)N +1 (ξ + η) 2η2N +3 (βα)N (β + α)dQ = − (N + 1)(N + 2) (N + 1)(N + 2) η 0  η β η2N +3 (βα)N (β + α)dQ = (N + 1)(N + 2) 0 0  ξ η (α(β + 4))N (β + α + 4)dQ

T1 ◦ QN (N !)(N + 1)! N +1 (1 + kσ )   σ M l+N +2 (l + N + 1)! k=0 T2 ◦ QN + 2 N +2 4 (N + 1)! (N !)(N + 1)! N  (1 + kσ )  σ l +N +2 M (l + N + 1)! k=0 + a(N ) N +2 4 (N + 1)! (N !)(N + 1)!  η ×4 QN (β, β)dβ 4N +2

Definition 19 (Rodino, 1993). The function f (t ) ∈ GD,σ (Ω ) is in the Gevrey class of order σ in Ω , if f (t ) ∈ C∞ and there exists a positive constant D such that supt ∈Ω ∂tn f (t ) ≤ Dn+1 (n!)σ for all n ∈ N0 .

Lemma 20. With dQ = dα dβ it holds for σ ≥ 1 and l, k, N ∈ N

sup ∂tl GN +1 (ξ , η, t )



j

Appendix C. Some useful definitions and lemmas

 l−k   k  ∂ γ  ∂ GN  t

k

Dirichlet or Neumann BC and M < 4/(Lj (Lj + 4)) since η(ξ + 4) ≤ Lj (Lj +4) for mixed BC in the considered domain Ωg¯ . This completes the proof. 



4N +2

(l + N + 1)! (N + 1)!

σ



(1 + k σ )

k=0

(N + 1)!(N + 2)!

QN +1 ,

which is equal to (54) for the mixed BC. Thus it is proven that (54) holds for all n ∈ N0 . In a second part of the proof the upper bound of the series (52) can be estimated using (54) evaluated at l = 0, i.e.,

|¯g (ξ , η, t )| ≤

n 

(1 + kσ )

∞ 

M

n =0

4n+1 (n!)(n + 1)!

n+1

k=0

2ηN +2 (η + 4)N +1

+

ηN +1 (ξ + 4)N +1 (ξ + η + 4) (N + 1)(N + 2)

(N + 1)(N + 2) 4ηN +1 (η + 4)N +1 − (N + 1)(N + 2)  η β ηN +2 (η + 4)N +1 α N (β + 4)N (β + α + 4)dQ = (N + 1)(N + 2) 0 0  η N +1 η (η + 4)N +1 β N (β + 4)N (2β + 4)dQ = . N +1 0 References

Qn (ξ , η)

n   ∞  M M ηξ   (ξ − η) bn , Dirichlet BC   4 4   n =0    n ∞    M ηξ M (ξ + η) bn , Neumann BC ≤ 4 4 n =0   n   ∞   M M η(ξ + 4)    (ξ + η + 4 ) b , n   4  n =0 4 mixed BC, n where bn = n!(n1+1)! k=0 (1 + kσ ). Thereby, the radius of convergence ρ of the latter inequality is determined by applying d’Alembert’s ratio test, see, e.g., Zwillinger (2002). Together with L’Hospital’s rule it follows that ρ = limn→∞ b bn = n+1

+3 limn→∞ σ (n2n . Hence, this results in ρ = ∞ if 1 ≤ σ < 2 +1)σ −1 and ρ = 1 if σ = 2. Obviously, the series converges with infinite radius of convergence if σ ∈ [1, 2). If σ = 2, convergence requires that M ηξ < 4 for Dirichlet or Neumann BC and M η(ξ + 4) < 4 for mixed BC, respectively. This leads to M < 4/L2j since ηξ ≤ L2j for

Akkari, E., Chevallier, S., & Boillereaux, L. (2006). Observer-based monitoring of thermal runaway in microwaves food defrosting. Journal of Process Control, 16(9), 993–1001. Armaou, A., & Demetriou, M. A. (2006). Optimal actuator/sensor placement for linear parabolic PDEs using spatial H2 norm. Chemical Engineering Science, 61(22), 7351–7367. Bitzer, M., & Zeitz, M. (2002). Design of a nonlinear distributed parameter observer for a pressure swing adsorption plant. Journal of Process Control, 12(4), 533–543. Colton, D. (1977). The solution of initial-boundary value problems for parabolic equations by the method of integral operators. Journal of Differential Equations, 26, 181–190. Dochain, D., Tali-Maamar, N., & Babary, J. (1997). On modelling, monitoring and control of fixed bed bioreactors. Computers & Chemical Engineering, 21(11), 1255–1266. Fahroo, F., & Demetriou, M. A. (2000). Optimal actuator/sensor location for active noise regulator and tracking control problems. Journal of Computational and Applied Mathematics, 114(1), 137–158. Jadachowski, L., Meurer, T., & Kugi, A. (2011). State estimation for parabolic PDEs with varying parameters on 3-dimensional spatial domains. In Proceedings of the 18th IFAC world congress. Milano, Italy, August 28–September 2 (pp. 13338–13343). Jadachowski, L., Meurer, T., & Kugi, A. (2012). An efficient implementation of backstepping observers for time-varying parabolic PDEs. In 7th MATHMOD conference, Vienna, Austria, February 14–17 (pp. 798–803). [Online]. Available: http://www.ifac-papersonline.net/58797.html. Joseph, K. T., & Vasudeva Murthy, A. S. (2001). Hopf-Cole transformation to some systems of partial differential equations. NoDEA: Nonlinear Differential Equations and Applications, 8, 173–193.

L. Jadachowski et al. / Automatica 51 (2015) 85–97 Krstic, M., & Smyshlyaev, A. (2008). Boundary control of PDEs: a course on backstepping designs. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics. Ladyžhenskaia, O. A., Solonnikov, V. A., & Ural’ceva, N. N. (1988). Translations of mathematical monographs: Vol. 23. Linear and quasi-linear equations of parabolic type. Providence, Rhode Island: American Mathematical Society. Liu, W. (2003). Boundary feedback stabilization of an unstable heat equation. SIAM Journal on Control and Optimization, 42(3), 1033–1043. Meurer, T., & Kugi, A. (2009). Tracking control for boundary con-trolled parabolic PDEs with varying parameters: combining backstepping and differential flatness. Automatica, 45(5), 1182–1194. Rodino, L. (1993). Linear partial differential operators in Gevrey spaces. Singapore: World Scientific Publishing Co. Pvt. Ltd. Smyshlyaev, A., & Krstic, M. (2005). Backstepping observers for a class of parabolic PDEs. Systems & Control Letters, 54(7), 613–625. Temple, G. (1928). The theory of Rayleigh’s principle as applied to continuous systems. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 119(782), 276–293. Vande Wouwer, A., Point, N., Porteman, S., & Remy, M. (2000). An approach to the selection of optimal sensor locations in distributed parameter systems. Journal of Process Control, 10(4), 291–300. Vande Wouwer, A., & Zeitz, M. (2001). State estimation in distributed parameter systems. In H. Unbehauen (Ed.), Encyclopedia of life support systems (EOLSS). Oxford, UK: EOLSS Publishers, (Chapter) Control systems, robotics and automation, Article No. 6.43.19.3. Vazquez, R., Schuster, E., & Krstic, M. (2008). Magnetohydro-dynamic state estimation with boundary sensors. Automatica, 44(10), 2517–2527. Zwillinger, D. (2002). Discrete mathematics and its applications, CRC standard mathematical tables and formulae (31st ed.). Taylor & Francis.

Lukas Jadachowski received the Dipl.-Ing. degree in production engineering from the Saarland University, Saarbrücken, Germany, in 2008, and the Ph.D. degree from the Vienna University of Technology, Vienna, Austria, in 2013. Since May 2013, he has been a Post-Doc with Automation and Control Institute, Vienna University of Technology, Vienna, Austria. His research interests include the physics-based modeling, backstepping-based observer and control design for distributed-parameter systems.

97

Thomas Meurer received the diploma in chemical engineering from the University of Stuttgart, Stuttgart, Germany, in 2001, the M.Sc. degree in engineering science and mechanics from the Georgia Institute of Technology, Atlanta (GA), USA, in 2000, and the Ph.D. degree from the University of Stuttgart in 2005. From 2005 to 2007, he was a Post-Doc with the Chair of System Theory and Automatic Control at the Saarland University, Saarbrücken, Germany. In 2007, he joined the Automation and Control Institute at the TU Vienna, Vienna, Austria, as a Senior Researcher and in 2012 he was appointed Associate Professor at TU Vienna. Since November 2012, he has been a Full Professor and Head of the Chair of Automatic Control at the Christian-Albrechts-University Kiel, Kiel, Germany. His research interests include feedback control and trajectory planning for distributedparameter systems, differential flatness, nonlinear control theory, and observer design with applications in smart structures and process engineering. Since 2011, he has been Chair of the IFAC Technical Committee 2.6 on Distributed Parameter Systems and serves as Associate Editor for the IFAC journals Automatica and Control Engineering Practice.

Andreas Kugi (M’94) received the Dipl.-Ing. degree in electrical engineering from Graz University of Technology, Austria, and the Ph.D. degree in control engineering from Johannes Kepler University (JKU), Linz, Austria, in 1992 and 1995, respectively. He received his Habilitation degree in the field of automatic control and control theory from JKU in 2000. From 1995 to 2000 he worked as an assistant professor and from 2000 to 2002 as an associate professor at JKU. In 2002, he was appointed full professor at Saarland University, Saarbrücken, Germany, where he held the Chair of System Theory and Automatic Control until May 2007. Since June 2007 he is a full professor for complex dynamical systems and head of the Automation and Control Institute at Vienna University of Technology, Austria. Since 2010, he is corresponding member of the Austrian Academy of Sciences. His research interests include the physics-based modeling and control of (nonlinear) mechatronic systems, differential geometric and algebraic methods for nonlinear control, model predictive control and the controller design for infinite-dimensional systems. He is Editor-in-Chief of the Control Engineering Practice.