Discretization and goal oriented error estimation

For the discretization of the infinite dimensional problem, we use a discontinuous Galerkin approach of order zero in time (denoted by dG(0)), and a continuous Galerkin approach of order one in space (denoted by cG(1)) as presented in [100,101]. In the literature, this combined approach is often referred to as dG(0)cG(1)-discretization. We will briefly recall some of the work considering this discretization technique.

5.1. SETTING AND PRELIMINARIES

Discretization and adaptivity for parabolic equations with discontinuous Galerkin methods was first established in the seminal papers [45, 46]. For the particular case of ¯A(x) = ∆x, a priori time and space discretization error estimates for optimal control of parabolic PDEs of order k+ 1 and s+ 1, respectively, are given in [102, Section 5.1], where kand sare the orders of the polynomials in the ansatz space in time and space, respectively. Control constraints were included in [131]. For semilinear parabolic PDEs, a priori bounds were obtained in [106] under growth conditions, whereas the case of semilinear parabolic PDEs without growth conditions was treated recently in [103]. Considering efficient numerical realization, the reader is referred to [120] for a PDE context and to [16] for the case of optimal control. Lastly, there are recent discrete maximal parabolic regularity results for the discrete-time equations, cf. [93,94].

For the reader’s convenience, we will briefly recall the definition of this discretization scheme and the corresponding a posteriori error goal oriented estimation. In the following, we will abbreviate

W :=W([0, T]), U =L₂(0, T;U), hv, wi_I :=

Z

I

hv(t), w(t)i_V^∗×V dt.

Time discretization

We split up the interval [0, T] ={0} ∪I1∪I2 ∪ · · · ∪I_M into subintervals Im = (tm−1, tm] of corresponding sizek_m :=t_m−t_m−1form∈ {1, . . . , M}and setI₀ :={0}, where 0 =t₀ < t₁· · ·<

tM =T. We define the discrete-time spaces of piecewise constant in time ansatz functions by W_k:={v_k∈L₂(0, T;H)|v_k

Im∈V, m= 1, . . . , M, v_k(0)∈H}, U_k:={u_k ∈L₂(0, T;U)|u_k

Im

∈U, m= 1, . . . , M}.

By continuity of elements in W = W([0, T]) ,→ C(0, T;H), cf. Lemma 3.4, this forms a non-conforming ansatz space, as elements ofW_kare not necessarily continuous. However, despite the nonconformity, the important feature of Galerkin orthogonality of the difference of continuous and discrete solution to the test space is preserved, cf. [100, Remark 5.2]. To capture the possible discontinuities, we denote the right and left sided limits and the jump at time grid point tm for v_k∈ W_k via

v_k,m⁺ := lim

t→0⁺v_k(t_m+t), v_k,m⁻ := lim

t→0⁺v_k(t_m−t), [v]_k,m:=v⁺_k,m−v_k,m⁻ , and illustrate this definition in Figure 5.1.

tm−1 t_m tm+1

v_k,m⁻ v_k,m⁺

[v_k]_m

Figure 5.1: One sided limits and jumps of discrete-time variables.

Due to the nonconformity of the ansatz space, the Lagrange function defined in (5.2) is not defined onW_k. Thus, we define the discrete-time Lagrange functionL^k:W_k× U_k× W_k→Rby

where the jump terms [x_k]m−1 capture the discontinuities of the state. This Lagrange function is also well defined for state and adjoint state belonging to the continuous function space W and on this space it coincides with the continuous-time Lagrangian defined in (5.2). For piece-wise constant functions of the space W_k, the time derivative vanishes, whereas for functions continuous in time belonging toW, the jump terms vanish.

The discrete-time version for the state equation of (5.3) reads hL^k_λ(xk, uk, λk), ϕki_W^∗

forϕk∈ W_k. Analogously, the discrete-time counterpart to the third equation of (5.3), is given by forϕ_k∈ U_k. Using integration by parts on each subinterval in the state equation (5.6), one can derive the adjoint equation as discrete-time counterpart to the first equation of (5.3), that is,

hL^k_x(x_k, u_k, λ_k), ϕ_ki_W^∗

5.1. SETTING AND PRELIMINARIES

for allϕ_k∈ W_k. The resulting time-stepping scheme is equivalent to an implicit Euler method if the temporal integrals are approximated via the box rule, cf. [100, Section 3.4.1], and thus inherits its A-stability.

Space discretization and time-stepping on dynamic meshes

For spatial discretization we use linear continuous finite elements as treated in the standard literature [23, 32, 77]. To this end, we assign a regular triangulation K^m_h and corresponding conforming finite element spaces V_h^m⊂V and U_h^m⊂U to each intervalIm and obtain the fully discrete spaces

W_kh :={v_kh∈L₂(0, T, H)|v

kh Im

∈V_h^m, m= 1, . . . , M, v_kh(0)∈V_h⁰}, U_kh :={u_kh ∈L₂(0, T, U)|u

kh Im

∈U_h^m, m= 1, . . . , M}. (5.9) Due to conformity of these spaces with respect to the discrete-time spaces, i.e., W_kh⊂ W_k and U_kh ⊂ U_k, the discrete-time Lagrangian (5.5) is well defined onW_kh× U_kh× W_kh.

In order to allow full flexibility for the spatial adaptivity, it is possible that the triangulation K_h^m on the interval I_m is different from the triangulation K^m−1_h on the interval Im−1. In terms of numerical realization, this leads to difficulties in efficiently evaluating the scalar product of basis elements of different time steps as needed for the assembly of the Euler step equations (5.6) and (5.8). A remedy is presented in [128], where the authors suggest the evaluation of scalar products on a common triangulation of K^m_h and K_h^m−1, which we denote by K^m−1/2_h . This common triangulation is depicted in Figure 5.2, where the original meshes have been independently red-green refined, cf. [40, Section 6.2.2] and [11]. If both meshes stem from the same original mesh by refinement, then the common refinement leads to a regular triangulation and to a finite element spaceV_h^m−1/2 such thatV_h^m−1, V_h^m ⊂V_h^m−1/2.

K^m−1_h K^m−1/2_h K_h^m

Figure 5.2: Sketch of common refinementK^m−1/2_h of two triangulations K^m−1_h and K_h^m. In our case, this common refinement is computed by the module dune-gridglue [14] of the DUNE C++-library [20] and allows us to compute scalar products of basis elements ψm∈V_h^m

and ψm−1∈V_h^m−1 via

Z

Ω

ψmψm−1= X

K∈K^m−1/2_h

Z

K

ψmψm−1. (5.10)

By construction of the grids, for each cell K ∈ K^m−1/2_h , there are corresponding parent cells K^m ∈ K^m_h and K^m−1 ∈ K_h^m−1 such that K ⊂K^m and K ⊂K^m−1. Thedune-gridglue module provides the index of the associated parent cellsK^mandK^m−1 in each original mesh. Thus, the integral over cells of the commonly refined triangulation in (5.10) can be evaluated efficiently with local evaluation in K^m and K^m−1 and a suitable quadrature rule. Hence, the price to pay for dynamic space grids is the computation of the common triangulations and the assembly of M − 1 mass matrices, assigning to finite element functions defined on one space grid a linear functional on a neighboring space grid. The algorithms completing these tasks can be implemented in parallel using all available CPU-cores. Further, after refinement of space grid K_h^m, only the common refinements K^m−1/2_h and K^m+1/2_h and the corresponding mass matrices need to be updated. We will discuss this topic in detail inSection 5.3.4.

Goal oriented error estimation

We will now introduce the concept of goal oriented error estimation for optimal control of parabolic PDEs. There are a lot of works considering goal oriented error estimation starting with the seminal papers [15, 17, 18], which were extended to systems with state or control constraints [116], optimal control of hyperbolic equations [86] and optimal control of parabolic equations [100, 101, 102]. A comprehensive introduction to adaptive finite element methods for ODEs and PDEs with applications is given in the monograph [10]. The main idea of goal oriented error estimation is to estimate and reduce the discretization error with respect to an arbitrary functionalI(x, u), called the quantity of interest (QOI). Motivations for the definition of QOIs range from allowing error estimation outside of the usual energy norm for, e.g., flow simulation in the PDE case [18, 78] to the case of optimal control, where applications include parameter estimation and optimal choice of regularization parameters [102,141] to the standard case of choosing the cost functional as the QOI.

We follow the literature [100,101] and denote by (x, u, λ)∈(W × U × W) a continuous-time solution of the extremal equations (5.3), by (x_k, u_k, λ_k)∈(W_k×U_k×W_k) and by (x_kh, u_kh, λ_kh)∈ (W_kh× U_kh× W_kh) time and fully discrete solutions of the system described by (5.6), (5.7), and (5.8). One intermediate aim of goal oriented a posteriori error estimation is to derive error estimatorsη_k andη_h such that

I(x, u)−I(xkh, ukh)≈ηk+ηh,

whereη_kapproximates the time discretization error andη_happroximates the space discretization error. A detailed derivation of the estimators is performed in [100, Chapter 6] and [101]. We briefly recall the main steps for the convenience of the reader and for later use. For more

5.1. SETTING AND PRELIMINARIES

details, the interested reader is referred to the references above. Besides the solution triple ξ:= (x, u, λ) of the first-order necessary conditions, a second triple of variablesχ:= (v, q, z) has to be considered. Thesesecondary variables solve the linear system

L⁰⁰(ξ)χ= (L^k)⁰⁰(ξ)χ=−

on the continuous-time level, the system

(L^k)⁰⁰(ξ_k)χ_k=−

on the discrete-time level, and the system

(L^k)⁰⁰(ξ_kh)χ_kh=−

on the fully discrete level. These equations are similar to the defining equation of a Lagrange-Newton step, where the derivative of the Lagrangian on the right hand side is replaced by the derivative of the QOI.

With the continuous triples ξ = (x, u, λ) and χ = (v, q, z) and the corresponding discrete counterparts, we define the residual of the first-order optimality condition via

ρ^λ(x, u, λ)ϕ:=hL^k_x(x, u, λ), ϕi_W^∗ and a residual involving the secondary variables χ= (v, q, z) via

ρ^z(ξ, v, q, z)ϕ:=L^k_λx(ξ)(z, ϕ) +L^k_ux(ξ)(q, ϕ) +L^k_xx(ξ)(v, ϕ) +I_x⁰(x, u)ϕ, ρ^q(ξ, v, q, z)ϕ:=L^k_uu(ξ)(q, ϕ) +L^k_xu(ξ)(v, ϕ) +L^k_λu(ξ)(z, ϕ) +I_u⁰(x, u)ϕ,

ρ^v(ξ, v, q)ϕ:=L^k_xλ(ξ)(v, ϕ) +L^k_uλ(ξ)(q, ϕ).

With these residuals, the time discretization error can be estimated via I(x, u)−I(x_k, u_k)≈

1

2 ρ^λ(x_k, u_k, λ_k)(v−v_k) +ρ^u(x_k, u_k, λ_k)(q−q_k) +ρ^x(x_k, u_k)(z−z_k)

+ρ^z(ξ_k, v_k, q_k, z_k)(x−x_k) +ρ^q(ξ_k, v_k, q_k, z_k)(u−u_k) +ρ^v(ξ_k, v_k, q_k)(λ−λ_k)

for (v_k, q_k, z_k),(x_k, u_k, λ_k) ∈ W_k× U_k× W_k arbitrary. Similarly, the space discretization error estimator can be approximated via

I(x_k, u_k)−I(x_kh, u_kh)≈ 1

2 ρ^λ(x_kh, u_kh, λ_kh)(v_k−v_kh) +ρ^u(x_kh, u_kh, λ_kh)(q_k−q

kh) +ρ^x(xkh, ukh)(zk−z_kh) +ρ^z(ξkh, vkh, qkh, zkh)(xk−x_kh) +ρ^q(ξ_kh, v_kh, q_kh, z_kh)(u_k−u_kh) +ρ^v(ξ_kh, v_kh, q_kh)(λ_k−λ_kh)

for (v_kh, q_kh, z_kh),(x_kh, u_kh, λ_kh)∈ W_kh× U_kh× W_kh. The arbitrary choice of the test functions originates in Galerkin orthogonality, cf. [101, Proposition 4.1, Theorem 4.3]. The terms v−v_k, q−q_k, z−z_k, x−x_k, u−u_k, λ−λ_k resp. vk−v_kh, qk−q_kh, zk−z_kh, xk−x_kh, uk−u_kh, λ_k−λ_kh are often called weights and need to be approximated to obtain computable error estimates as the solutions in the infinite dimensional spaces, i.e., expressions with no subscript or subscript k, are not at hand. Approximating the weights by elements of W_k and W_kh, respectively, causes the estimators to vanish due to Galerkin orthogonality. Hence, will discuss options to efficiently approximate the weights inSection 5.3.4. Having approximated the weights for the time discretization error by w^k_v, w_q^k, w_z^k, w^k_x, w^k_u and w_λ^k and the weights for the space discretization error by w^h_v,w^h_q,w_z^h,w_x^h,w^h_u and w^h_λ we define the error indicators by

η_k :=1

2 ρ^λ(x_kh, u_kh, λ_kh)(w_v^k) +ρ^u(x_kh, u_kh, λ_kh)(w^k_q) +ρ^x(x_kh, u_kh)(w^k_z)

+ρ^z(ξ_kh, v_kh, q_kh, z_kh)(w_x^k) +ρ^q(ξ_kh, v_kh, q_kh, z_kh)(w_u^k) +ρ^v(ξ_kh, v_kh, q_kh)(w^k_λ)

(5.14) and

η_h :=1

2 ρ^λ(x_kh, u_kh, λ_kh)(w_v^h) +ρ^u(x_kh, u_kh, λ_kh)(w_q^h) +ρ^x(x_kh, u_kh)(w_z^h)

+ρ^z(ξkh, vkh, qkh, zkh)(w_x^h) +ρ^q(ξkh, vkh, qkh, zkh)(w_u^h) +ρ^v(ξkh, vkh, qkh)(w^h_λ) .

(5.15)

Im Dokument Sensitivity Analysis and Goal Oriented Error Estimation for Model Predictive Control (Seite 131-137)