Located control and control constraints

uniform, FE uniform, ROM uniform, ROM adaptive, ROM

solve eigs(Z¯Z¯^T) solve eigs(Z¯^TZ)¯ solve eigs(Z¯^TZ)¯

n 30150 30150 17500 27

err_y 5.77·10⁻⁰³ 5.77·10⁻⁰³ 1.50·10⁻⁰² 5.77·10⁻⁰³ optimal control solution 10.30 sec 0.51 sec 0.40 sec 0.01 sec

snapshot generation – 0.75 sec 0.43 sec 0.002 sec

POD basis – 0.21 sec 70.18 sec 0.09 sec

time-adaptive grid – – – 0.61 sec

Table 1: Accuracy and computational times for the solution of the optimal control problem using a uniform versus an adaptive time discretization according to Algorithm 4 comparing a full finite element method with spatial discretization h^∗ = 1/100 with a POD reduced-order model using

`= 1 POD basis function. The POD basis is computed solving a weighted eigenvalue problem (eigs) for ¯ZZ¯^T or ¯Z^TZ¯

subset of L^∞(0, T;R^m). We introduce a control operator as B:U →L²(0, T;H⁻¹(Ω)), u7→ Bu=

m

X

i=1

u_iχ_i,

where χ_i ∈H⁻¹(Ω),1≤i≤m, denote given shape functions. Thus,B is a linear and bounded operator. The dual operator B^∗:L²(0, T;H⁻¹(Ω))^∗ →U^∗ fulfills

hBu, vi_L2(0,T;H⁻¹(Ω)),L²(0,T;H₀¹(Ω)) = Z T

0

hBu(t), v(t)i_H−1(Ω),H₀¹(Ω)dt

=

m

X

i=1

Z _T

0

u_i(t)hχ_i, v(t)i_H−1(Ω),H₀¹(Ω)dt

=

m

X

i=1

Z T 0

ui(t)(B^∗v(t))idt

= (u,B^∗v)_L²_(0,T;R^m₎=hu,B^∗vi_U

according to Definition 2.9. Thus, it holds (B^∗v(t))_i = hχ_i, v(t)i_H−1(Ω),H₀¹(Ω). Note that we identify L²(0, T;H⁻¹(Ω))^∗ withL²(0, T;H₀¹(Ω)) and U^∗ withU. In the state equation (3.1) we replace the distributed control by the located control as follows

y_t−∆y=Bu in (0, T]×Ω, (3.54a)

y= 0 on [0, T]×∂Ω, (3.54b)

y(0,·) =y0 in Ω. (3.54c)

Existence of a unique weak solution y ∈ W(0, T;H₀¹(Ω)) to (3.54) and the regularity results follow according to the theory for (3.1). We introduce the weak solution operator associated with equation (3.54) by

S_ad :Uad×L²(Ω)→W(0, T;H₀¹(Ω)), (u, y0)7→y:=S_ad(u, y0). (3.55) As cost functional we consider the quadratic objective

J_ad(y, u) = 1

2ky−y_dk²_L2(0,T;Ω)+ α

2kuk²_L2(0,T;R^m). The optimal control problem reads as

min

(y,u)∈W(0,T;H₀¹(Ω))×U_adJ_ad(y, u) s.t. y=S_ad(u, y₀). (3.56)

3.7 Located control and control constraints 59

Due to the control constraints, the gradient equation (3.11) in the optimality system is replaced by the variational inequality

(αu+B^∗p, v−u)U ≥0 for all v∈U_ad. (3.57) The resulting first-order necessary (and by convexity sufficient) optimality system is given by the state equation (3.54) together with the adjoint equation

−p_t−∆p=y−y_d in [0, T)×Ω, (3.58a)

p= 0 on [0, T]×∂Ω, (3.58b)

p(T,·) = 0 in Ω, (3.58c)

and the variational inequality (3.57) which is equivalent to the projection formula u(t) =P[ua(t),ub(t)]

−1

α(B^∗p)(t)

for almost allt∈[0, T], (3.59) where we define by

P[ua(t),u_b(t)] :R^m →[ua(t), u_b(t)]

for almost all t ∈ [0, T] the componentwise orthogonal projection onto the admissible set of control vectors. The projection can be computed as

P[ua(t),ub(t)]

−1

α(B^∗p)(t)

= max

u_a(t),min

u_b(t),−1

αhχ_i, p(t)i_H−1(Ω),H₀¹(Ω)

, which is meant componentwise for almost all t∈[0, T]. Moreover, we denote by

P[ua,ub]{u}:=P[ua(t),ub(t)]{u(t)}, for almost all t∈[0, T]

the pointwise projection of a function u ∈ L^∞(0, T;R^m) onto the set of admissible controls u_a, u_b∈L^∞(0, T;R^m).

Reformulation of the optimality system with respect to the adjoint state

We now can only reduce the optimality system to an elliptic system inp. The derivation of the biharmonic system is given later in the proof of Theorem 3.16. The resulting biharmonic system reads as

−p_tt+ ∆²p− B P[ua,ub]

−1 αB^∗p

= −(y_d)_t+ ∆y_d in (0, T)×Ω, (3.60a)

p= 0 on [0, T]×∂Ω, (3.60b)

∆p=yd on [0, T]×∂Ω, (3.60c)

p(T,·) = 0 in Ω, (3.60d)

(p_t+ ∆p)(0,·) =y_d(0,·)−y₀ in Ω. (3.60e) Note that a reduction with respect to the state variable is not possible in the standard way since the projection operator P[ua,ub]is a non-smooth operator, thus non-differentiable in the classical sense. One possible way to derive a biharmonic equation for the state variable is to consider a regularization of the projection, see [139]. This option is not further discussed here and we continue with the focus on the biharmonic equation with respect to the adjoint variable.

Lemma 3.14. Let us denote N(p) :=−B P[ua,ub]

−1 αB^∗p

. It holds true Z T

0

Z

Ω

(N(p1)−N(p2))(p1−p2)dxdt≥0 ∀p₁, p2.

Proof. Since the projection operator is monotone increasing, i.e.

u₁(t)≤u₂(t) ⇒ P[ua(t),ub(t)]{u₁(t)} ≤P[ua(t),ub(t)]{u₂(t)}, t∈[0, T] foru₁, u₂ ∈L^∞(0, T;R^m) it is easy to see that

P[ua(t),ub(t)]{u₁(t)} −P[ua(t),ub(t)]{u₂(t)}, u₁(t)−u2(t)

R^m ≥0

We set u1(t) :=−_α¹(B^∗p1)(t) andu2(t) :=−_α¹(B^∗p2)(t) for p1, p2 ∈H˜₀^2,1(0, T; Ω). Then it holds

−P[ua(t),u_b(t)]

−1

α(B^∗p₁)(t)

+P[ua(t),u_b(t)]

−1

α(B^∗p₂)(t)

,(B^∗p₁)(t)−(B^∗p₂)(t)

R^m

≥0.

Hence by linearity of B^∗ it follows h−BP[ua(t),ub(t)]

−1

α(B^∗p1)(t)

+BP[ua(t),ub(t)]

−1

α(B^∗p2)(t)

, p1(t)−p2(t)i_H−1(Ω),H₀¹(Ω) ≥0 which leads to the claim.

Thus, the quantity

−BP[ua,ub]

−1 αB^∗p

is monotone inp and (3.60) is a semilinear elliptic problem with monotone nonlinearity.

In order to derive the weak formulation for (3.60) we define the operator A^ad₀ as A^ad₀ : ˜H₀^2,1(0, T; Ω)×H˜₀^2,1(0, T; Ω)→R,

A^ad₀ (p, v) = Z T

0

Z

Ω

p_tv_t− BP[ua,ub]

−1 αB^∗p

v+ ∆p∆v

dxdt+ Z

Ω

∇p(0)∇v(0)dx.

Then the weak form of equation (3.60) reads as

A^ad₀ (p, v) =L0(v) ∀v∈H˜₀^2,1(0, T; Ω), (3.61) where the linear form L0 is defined in (3.22).

Theorem 3.15. (Existence of a unique weak solution). For every L₀ ∈ ( ˜H₀^2,1(0, T; Ω))^∗ the equation (3.61)has a unique solution p∈H˜₀^2,1(0, T; Ω).

Proof. In analogy to [139, Lemma 4.6] it can be shown that the operator A^ad₀ is strongly mono-tone, coercive and hemi-continuous. Then, the main theorem on monotone operators (Theo-rem 2.25) ensures existence of a unique solution.

Equivalence of the optimality system to a biharmonic equation

Theorem 3.16. Let (y, u) ∈W(0, T;H₀¹(Ω))×U_ad denote the solution to problem (3.56) with associated adjoint statep∈W(0, T;H₀¹(Ω)). Assume thaty₀ ∈H₀¹(Ω)andy_d∈L²(0, T;H²(Ω))∩

H¹(0, T;L²(Ω)). Further we assume the compatibility condition y(T)−y_d(T) + ∆p(T)∈H₀¹(Ω) to hold true. Then p satisfies (3.60)a.e. in space-time and is a weak solution to (3.60).

3.7 Located control and control constraints 61

Proof. Due to y₀∈H₀¹(Ω) we follow by Remark 3.2 thaty∈L²(0, T;H²(Ω))∩H¹(0, T;L²(Ω)).

Then, we can apply Remark 3.3 which leads to the regularity

p∈L²(0, T;H⁴(Ω))∩H¹(0, T;H²(Ω))∩H²(0, T;L²(Ω)).

Thus, we can take in (3.58a) the derivative with respect to time in order to get

−p_tt−∆pt=yt−(yd)t. Then, we replace y_t using (3.54a) and obtain

−p_tt−∆pt= ∆y+Bu−(y_d)t.

We use the projection formula (3.59) in order to eliminate the control variable and get

−p_tt−∆p_t= ∆y+BP[ua,u_b]

−1 αB^∗p

−(y_d)_t.

Finally, we use the identityy =−p_t−∆p+y_dfrom (3.58a) in order to eliminate the state variable and get

−p_tt−∆p_t= ∆(−p_t−∆p+y_d) +BP[ua,ub]

−1 αB^∗p

−(y_d)_t,

which implies (3.60a). We get the initial condition (3.60e) by evaluating (3.58a) at t= 0. The boundary condition (3.60c) is derived by evaluating (3.58a) on the boundary ∂Ω.

Thus, the adjoint function p ∈ L²(0, T;H⁴(Ω))∩H¹(0, T;H²(Ω))∩H²(0, T;L²(Ω)) fulfills the biharmonic problem (3.60) a.e. in space-time. Then, p ∈ H˜₀^2,1(0, T; Ω) is a weak solution to (3.60).

Since the biharmonic equation (3.60) has a unique weak solution according to Theorem 3.15 and the weak solution p to (3.58) is a weak solution to (3.60) according to Theorem 3.16, the other direction holds, i.e. the weak solution to (3.60) is the optimal adjoint state p.

A-posteriori error estimate for the time discretization

We consider the same notation as in Section 3.4. The time-discrete problem for (3.61) reads as follows: find pk∈V˜k withpk(T,·) = 0 and

A^ad₀ (p_k, v_k) =L₀(v_k) ∀v_k∈V˜k. (3.62) The time-discrete problem has a unique solution by construction. We derive and prove a temporal residual-type a-posteriori error estimate for the variablepin the spirit of [82]. For this, we need the following results and properties.

Lemma 3.17. For α >0 and v∈H˜₀^2,1(0, T; Ω) it holds kvk_H2,1(0,T;Ω)≤c

kv_tk²_L2(0,T;Ω)+k∆vk²_L2(0,T;Ω)

with constant c >0.

Proof. We define the function g:=−v_t−∆v∈L²(0, T; Ω). Then, vis the solution to

−v_t−∆v = g in [0, T)×Ω, v = 0 on [0, T]×∂Ω, v(T,·) = 0 in Ω,

in a weak sense. This solutionv depends continuously on the data (compare e.g. (3.14)), thus kvk²_{H2,1(0,T;Ω)} ≤ ckgk²_L2(0,T;Ω)

= ck −vt−∆vk²_L2(0,T;Ω)

≤ 2c

kv_tk²_L2(0,T;Ω)+k∆vk²_L2(0,T;Ω) , where we applied Young’s inequality, Theorem 2.31.

Lemma 3.18. It holds

A^ad₀ (p, v_k)−A^ad₀ (p_k, v_k) = 0 ∀v_k ∈V˜k. Proof. With (3.61) and (3.62) we have

A^ad₀ (p, vk)−A^ad₀ (pk, vk) =L0(vk)−L0(vk) = 0 ∀v_k∈V˜k⊂H˜₀^2,1(0, T; Ω).

The following theorem states a temporal residual-type a-posteriori estimate forpwhich transfers the result in Theorem 3.7 to the case of located control and control constrained optimal control.

Theorem 3.19. Let p ∈ H˜₀^2,1(0, T; Ω) and pk ∈ V˜k denote the solutions to (3.61) and (3.62), respectively. Then, we obtain the a-posteriori error estimation

kp−p_kk²_H2,1(0,T;Ω)≤c η²_ad, (3.63) where c >0 and

η_ad² =

n

X

j=1

∆t²_j Z

Ij

−(y_d)t+ ∆y_d+ (p_k)tt+BPU_ad

−1 αB^∗p_k

−∆²p_k

2 L²(Ω)

+

n

X

j=1

Z

Ij

ky_d−∆pkk²_L2(∂Ω).

Proof. We set e^p := p−p_k, N(p) :=−B PUad

−¹_αB^∗p and denote by πe^p ∈ V˜k the Lagrange interpolation of e^p. We make the following estimations

1

ckp−pkk²_H2,1(0,T;Ω) ≤ k(p−pk)tk²_L2(0,T;Ω)+k∆(p−pk)k²_L2(0,T;Ω) (Lemma 3.17)

≤ Z T

0

Z

Ω

(p−p_k)_te^p_t + ∆(p−p_k)∆e^p+ (N(p)−N(p_k))e^pdxdt

(Lemma 3.14)

≤ Z T

0

Z

Ω

(p−pk)te^p_t + ∆(p−pk)∆e^p+ (N(p)−N(pk))e^pdxdt +k∇(p−p_k)(0)k²_L2(Ω)

= A^ad₀ (p, e^p)−A^ad₀ (p_k, e^p)

= A^ad₀ (p, e^p)−A^ad₀ (p_k, e^p)−(A^ad₀ (p, πe^p)−A^ad₀ (p_k, πe^p)) (Lemma 3.18)

= A^ad₀ (p, e^p−πe^p)−A^ad₀ (p_k, e^p−πe^p)

= L₀(e^p−πe^p)−A^ad₀ (p_k, e^p−πe^p)

3.7 Located control and control constraints 63

= Z T

0

Z

Ω

(−(y_d)t+ ∆yd)(e^p−πe^p)dxdt− Z

Ω

(yd(0)−y0)(e^p(0)−(πe^p)(0))dx +

Z _T

0

Z

∂Ω

y_d∇(e^p−πe^p)·~n dsdt− Z

Ω

∇p_k(0)∇(e^p−πe^p)(0)dx

− Z T

0

Z

Ω

((p_k)_t(e^p−πe^p)_t+N(p_k)(e^p−πe^p) + ∆p_k∆(e^p−πe^p))dxdt

= Z T

0

Z

Ω

(−(y_d)t+ ∆y_d+ (p_k)tt−N(p_k)−∆²p_k)(e^p−πe^p)dxdt +

Z T 0

Z

∂Ω

(y_d−∆p_k)∇(e^p−πe^p)·~n dsdt (integration by parts and Theorem 2.20).

Splitting the temporal integration into an integration over each time interval, using error es-timates of the Lagrange interpolation π, the trace inequality and Young’s inequality (Theo-rem 2.31), we get the claim.

Space-time mixed finite element discretization

We state the space-time mixed finite element variational form for (3.60) which is derived along the lines to (3.33) and uses the same notations. It reads as: find ph ∈Y˜₀^h, ϑh ∈V˜k andϑh|_∂Ω=−˜yd

such that Z T

0

Z

Ω

(p_h)_tv_t− BP[ua,u_b]

−1 αB^∗p_h

v+∇ϑ_h∇vdxdt+ Z

Ω

∇p_h(0)∇v(0)dx

= Z T

0

Z

Ω

(−(y_d)_t+ ∆y_d)v dxdt− Z

Ω

(y_d(0)−y₀)v(0)dx ∀v∈Y˜₀^h, (3.64a) Z T

0

Z

Ω

−∇p_h∇φ+ϑhφ dxdt= 0 ∀φ∈W₀^h. (3.64b) The adaptive finite element cycle which utilizes the error estimate (3.63) in order to adapt the time discretization according to the temporal dynamics of the adjoint state follows along the lines of Section 3.5 and is summarized in Algorithm 8.

Algorithm 8Space-time finite element algorithm with time adaptivity for the adjoint state in case of control constraints and located control.

Input: Initial space-time mesh T_h^init, tolerance ε >0, refinement parameter θ∈(0,1), desired stateyd, initial statey0, regularization parameterα >0, initial number of time discretization pointsn, maximal number of time pointsnmax.

Output: Time-adapted space-time meshTh, adjoint solution p_h.

1: Th =T_h^init with ntime points

2: while n < nmax ORηp ≥εdo

3: Define ˜Vk as the finite element space corresponding to Th according to (3.30).

4: Solve (3.64).

5: Estimate the error contributions from (3.63).

6: Mark time intervals according to the D¨orfler criterion (Definition 2.34) with parameter θ.

7: Refine the marked time intervals using bisection.

8: end while

Remark 3.20. Note that (3.64) is a nonlinear non-smooth equation due to the projection such that the solution can be done by a semi-smooth Newton method or a Newton method utilizing a regularization of the projection formula, for example. In the numerical tests, we will use a fixed point iteration in order to find an approximate solution.

3.8 Snapshot location in optimal control for POD model order reduction with located

Im Dokument Adaptivity in Model Order Reduction with Proper Orthogonal Decomposition (Seite 66-73)