Universität Konstanz
A method of a-posteriori error estimation with application to proper orthogonal decomposition
Eileen Kammann Fredi Tröltzsch Stefan Volkwein
Konstanzer Schriften in Mathematik Nr. 299, März 2012
ISSN 1430-3558
© Fachbereich Mathematik und Statistik Universität Konstanz
Fach D 197, 78457 Konstanz, Germany
Konstanzer Online-Publikations-System (KOPS) URL: http://nbn-resolving.de/urn:nbn:de:bsz:352-184660
Mod´elisation Math´ematique et Analyse Num´erique
A METHOD OF A-POSTERIORI ERROR ESTIMATION WITH APPLICATION TO PROPER ORTHOGONAL DECOMPOSITION∗
Eileen Kammann
1, Fredi Tr¨ oltzsch
2and Stefan Volkwein
3Abstract. We consider the following problem of error estimation for the optimal control of nonlinear parabolic partial differential equations: Let an arbitrary admissible control function be given. How far is it from the next locally optimal control? Under natural assumptions including a second-order sufficient optimality condition for the (unknown) locally optimal control, we estimate the distance between the two controls. To do this, we need some information on the lowest eigenvalue of the reduced Hessian. We apply this technique to a model reduced optimal control problem obtained by proper orthogonal decomposition (POD). The distance between a local solution of the reduced problem to a local solution of the original problem is estimated.
R´esum´e. ...
1991 Mathematics Subject Classification. 49K20, 35J61, 35K58.
The dates will be set by the publisher.
Introduction
In this paper, we focus on the following question for a class of optimal control problems for semilinear parabolic equations: Suppose that a numerical approximation ˜u for a locally optimal control is given. For instance, it might have been obtained by a numerical optimization method or as the solution to some reduced order optimization model. How far is this control from the nearest locally optimal control ¯u? Of course, we have to assume that such a solution ¯uexists in a neighborhood of ˜u. Moreover, ˜ushould already be sufficiently close to ¯u. The question is to quantify the errorku˜−uk¯ in an appropriate norm.
In principle, this is not a paper on proper orthogonal decomposition (POD). Our estimation method is not restricted to any specific method of numerical approximation for ˜u. Our primary goal is a numerical application of a perturbation method used by Arada et al. [2] in the context of FEM approximations of elliptic optimal control problems. The main idea of this method goes back to Dontchev et al. [7] and Malanowski and Maurer [18], who introduced it for the a priori error analysis of optimal control problems governed by ordinary differential equations.
Keywords and phrases:Optimal control, semilinear partial differential equations, error estimation, proper orthogonal decomposition
∗The first author was supported by SFB 910 ”Control of self-organizing nonlinear systems: Theoretical methods and concepts of application” in Berlin.
1 Institut f¨ur Mathematik, Technische Universit¨at Berlin, D-10623 Berlin, Germany (kammann@math.tu-berlin.de).
2 Institut f¨ur Mathematik, Technische Universit¨at Berlin, D-10623 Berlin, Germany (troeltzsch@math.tu-berlin.de).
3Institut f¨ur Mathematik und Statistik, Universit¨at Konstanz, D-78457 Konstanz, Germany (Stefan.Volkwein@uni-konstanz.de).
c EDP Sciences, SMAI 1999
However, we will apply our method to suboptimal controls ˜uobtained by POD. Though POD is an excellent method of model reduction for many time-varying or nonlinear differential equations, it lacks an a priori error analysis that is in some sense uniform in the right-hand side of the underlying partial differential equation, say with respect to the control function. Estimates of this type are known for the method of balanced truncation (see, e.g. [4]) that is limited to linear and autonomous PDEs. There are results on a priori estimates for POD that depend on certain assumptions on the orthogonal basis generated by the selected snapshots. We refer to Kunisch and Volkwein [14], Sachs and Schu [21], or [26]. However, such estimates will in general depend on the control used for generating the snapshots. In [10] an a-priori error analysis is presented for linear-quadratic optimal control problems. If the POD basis is computed utilizing the optimal state and associated adjoint variable, a convergence rate can be shown. But in real computation we do not know the optimal solution in advance. In view of this, we are interested in a posteriori estimates for assessing the precision of optimal controls for reduced control problems set up by POD. For the reduced-basis method a-posteriori error estimates for linear-quadratic optimal control problems we refer to [8].
We extend a method suggested in [26] for linear equations to the case of semilinear equations. For this purpose, we have to assume that ¯u satisfies a standard second-order sufficient optimality condition. The associated coercivity constant of the reduced Hessian operator will be estimated numerically.
Let us explain our idea for the following two optimal control problems in a bounded Lipschitz domain Ω⊂Rn with boundary Γ:
We consider thedistributed optimal control problem minJ(y, u) := 1
2 Z T
0
Z
Ω
(y(x, t)−yd(x, t))2+λ u(x, t)2 dxdt (PD) subject to the semilinear parabolic state equation
∂y
∂t −∆y(x, t) +d(x, t, y(x, t)) = u(x, t) in Ω×(0, T) y(x, t) = 0 in Γ×(0, T) y(x,0) = y0(x) in Ω
(0.1)
and to the pointwise control constraints
α≤u(x, t)≤β a.e. in Ω×(0, T).
Moreover, we deal with theboundary control problem minJ(y, u) :=1
2 Z T
0
Z
Ω
(y(x, t)−yd(x, t))2dxdt+λ 2
Z T 0
Z
Γ
u(x, t)2ds(x)dt (PB)
subject to the parabolic state equation
∂y
∂t −∆y(x, t) = 0 in Ω×(0, T)
∂y
∂ν(x, t) +b(x, t, y(x, t)) = u(x, t) in Γ×(0, T) y(x,0) = y0(x) in Ω
(0.2)
and to the pointwise control constraints
α≤u(x, t)≤β a.e. in Γ×(0, T).
Here, and throughout the paper,ν∈Rn denotes the outward unit normal on Γ.
Our numerical analysis is based on the following assumptions:
(A1)Ω⊂Rn is a bounded domain with Lipschitz boundary Γ.
(A2) The functionsd: Ω×[0, T]×R→R andb : Γ×[0, T]×R→Rare measurable with respect to the first two variables (x, t) and twice continuously differentiable with respect to the third variabley∈R.
The functions b and d are monotone non-decreasing with respect to the third variable y ∈ R for all fixed (x, t)∈Γ×[0, T] and Ω×[0, T], respectively.
(A3)There is a constantC∞>0 such that
|d(x, t,0)|+
∂d
∂y(x, t,0)
+
∂2d
∂y2(x, t,0)
≤C∞ for a.a. (x, t)∈Ω×[0, T] and
|b(x, t,0|+
∂b
∂y(x, t,0)
+
∂2b
∂y2(x, t,0)
≤C∞ for a.a. (x, t)∈Γ×[0, T].
Moreover, for all M >0 there exists a constantLM >0 such that
∂2d
∂y2(x, t, y1)−∂2d
∂y2(x, t, y2)
≤LM|y1−y2|
holds for a.a. (x, t)∈Ω×[0, T] and ally1, y2∈Rsatisfying|y1−y2|< M. We require the same condition for b for a.a. (x, t)∈Γ×[0, T].
(A4)A desired state functionyd ∈L2(Ω×(0, T)) and real constantsλ >0, α < β are given.
The paper is organized as follows: After this introduction, we explain the perturbation method and its use for the simpler case of a linear equation in Section 1. In Section 2, we extend the perturbation idea to semilinear equations. To make the paper also readable for readers being not familiar with POD, we briefly introduce the main concept of POD in Section 3, and Section 4 is devoted to the numerical application of our method to different nonlinear state equations.
1. A survey on the linear-quadratic case
1.1. The perturbation method for minimizing quadratic functionalsThe perturbation method was introduced by Dontchev et al. [7] and Malanowski and Maurer [18] in the context of optimal control of ordinary differential equations, where it was used to derive a priori error estimates for associated numerical approximations. In a different way, it was adopted by Arada et al. [2] for elliptic control problems. Let us explain this method for the following situation:
We consider the quadratic optimization problem in Hilbert space min
u∈Kf(u) := 1
2kSu−yHk2H+λ
2kuk2L2(D), (PQ)
where H is a real Hilbert space, D ⊂ Rm is a measurable bounded set, K ⊂ L2(D) is a nonempty, convex, closed and bounded set, S : L2(D) → H a continuous linear operator, λ >0 a fixed number, and yH ∈ H a fixed element. In this form, optimal control problems with partial differential equations are very frequently discussed. We refer only to the monography by Lions [16] or to the more recent books by Ito and Kunisch [12], Hinze et al. [9] or to the textbook [25]. We will mainly refer to the last reference, since the discussion of optimal control problems for semilinear parabolic equations in this book is very close to our notation.
Theorem 1.1. Under the assumptions formulated above, (PQ) has a unique optimal solution u¯ ∈ Uad. It satisfies the variational inequality
Z
D
{(S∗(Su¯−yH))(x) +λu(x)}¯ (u(x)−u(x))¯ dx≥0 ∀ u∈ K, (1.3)
where, S∗:H →L2(D)denotes the adjoint operator toS.
This is a standard result of the optimization theory, we refer e.g. to [16] or to [25, Thm. 2.22].
We introduce the adjoint state ¯passociated with ¯uby
¯
p(x) := (S∗(Su¯−yH))(x).
The variational inequality can be expressed equivalently in a pointwise way, namely (¯p(x) +λ¯u(x))¯u(x) = min
v∈[α,β](¯p(x) +λ¯u(x))v (1.4)
for almost allx∈D. This has important consequences. We have for almost allx∈D
¯ u(x) =
α, if ¯p(x) +λ¯u(x)>0
−¯p(x)/λ, if ¯p(x) +λ¯u(x) = 0 β, if ¯p(x) +λ¯u(x)<0.
While these inequalities show, how ¯u is determined by ¯p+λ¯u, the implications below concern the opposite direction; they follow directly from the discussion above. It must hold
¯
p(x) +λ¯u(x)≥0, if ¯u(x) =α
¯
p(x) +λ¯u(x) = 0, ifα <u(x)¯ < β
¯
p(x) +λ¯u(x)≤0, if ¯u(x) =β.
(1.5)
The last implications are essential for understanding the perturbation method. It is helpful to answer the following question: Let ˜u∈ Kbe a function that need not be optimal for (PQ). Let ˜pbe the associated adjoint state.
If ˜uwere optimal, then ˜p(x) +λ˜u(x) = 0 should be satisfied in almost all points x∈Ω, whereα <u(x)˜ < β holds. If not, then ˜p(x) +λ˜u(x) +ζ(x) = 0, if we define aperturbation function ζ in these points by
ζ(x) =−[˜p(x) +λ˜u(x)].
In the points, where ˜u(x) =αholds, the inequality ˜p(x) +λ˜u(x)≥0 should be satisfied for optimality. If not, then ˜p(x) +λ˜u(x) +ζ(x)≥0 is fulfilled with
ζ(x) = [˜p(x) +λ˜u(x)]−,
where [a]− is defined for a real numbera by [a]− = (|a| −a)/2. Analogously, we defineζ in the points, where
˜
u(x) =β. In this way, we obtain the following definition ofζthat is adopted from [2],
ζ(x) :=
[˜p(x) +λ˜u(x)]−, if ˜u(x) =α,
−[˜p(x) +λ˜u(x)], ifα <u(x)˜ < β,
−[˜p(x) +λ˜u(x)]+, if ˜u(x) =β.
(1.6)
The main idea behind the definition ofζ is the following: Although ˜uwill possibly not be optimal for (PQ), it is optimal for the perturbed optimization problem
minu∈Kf(u) := 1
2kSu−yHk2H+λ
2kuk2L2(D)+ Z
D
ζ(x)u(x)dx. (PQζ)
This follows from our construction, because ˜usatisfies the associated variational inequality Z
D
(˜p(x) +λ˜u(x) +ζ(x))(u(x)−u(x))˜ dx≥0 ∀u∈ K.
By convexity, ˜uis optimal for (PQζ).
Now, it is easy to estimate the distance between ˜uand ¯u.
Lemma 1.2. Suppose that u¯ is optimal for (PQ),u˜∈ K is given arbitrarily, andζ is defined by (1.6). Then there holds
ku˜−uk¯ L2(D)≤ 1
λkζkL2(D).
Proof. We just write down the variational inequalities satisfied by ¯u and ˜u, and insert the other function, respectively. This yields
(S∗(Su¯−yH) +λ¯u,u˜−u)¯ L2(D) ≥ 0 (S∗(Su˜−yH) +λ˜u,u¯−u)˜ L2(D)+ (ζ,u¯−u)˜ L2(D) ≥ 0.
Adding both inequalities, we arrive at
−kS(¯u−u)k˜ 2L2(D)−λku¯−uk˜ 2L2(D)+ (ζ,u¯−u)˜ L2(D)≥0, and hence, since kS(¯u−u)k˜ 2L2(D)≥0,
λk¯u−uk˜ 2L2(D)≤ kζkL2(D)k˜u−uk¯ L2(D).
Now the claimed result follows immediately.
Notice that, given ˜u,ζ can be computed. It is determined by the adjoint state ˜passociated with ˜u.
1.2. Application to a linear-quadratic control problem
As a simple application, we consider the following linear-quadratic version of (PD), minJ(y, u) := 1
2 Z T
0
Z
Ω
(y(x, t)−yd(x, t))2+λ u(x, t)2 dxdt (PL) subject to
∂y
∂t −∆y = u in Ω×(0, T) y = 0 in Γ×(0, T) y(·,0) = 0 in Ω
(1.7)
and to the pointwise control constraints
α≤u(x, t)≤β a.e. in Ω×(0, T).
We considery in the space
W(0, T) ={y∈L2(0, T;H1(Ω)) : dy
dt ∈L2(0, T;H1(Ω)0)},
cf. Lions [16]. For simplicity, we have assumed y0= 0. The case y0 6= 0 is reduced to that with homogeneous initial condition in a standard way: We shift the solution ˆy associated with homogeneous right-hand sideubut inhomogeneous initial datay0to yH. Then we consider ˆyH:=yH−yˆinstead of ˆy.
Let us introduce the notationQ:= Ω×(0, T) and Σ := Γ×(0, T) and define the set of admissible controls by
Uad={u∈L2(Q) :α≤u(x, t)≤β a.e. inQ}.
We are going to apply Lemma 1.2 with the choiceD:=Q,K:=Uad,yH :=yd.
The adjoint statep∈W(0, T) associated with a controluis obtained from theadjoint equation
−∂p
∂t −∆p = yu−yd in Ω×(0, T) p = 0 in Γ×(0, T)
p(·, T) = 0 in Ω,
(1.8)
whereyu denotes the state associated withu. It is a standard result that (PL) has a unique optimal control ¯u.
Let now ˜u∈Uad be given. For instance, we might think of a control obtained as the solution to a model reduced optimal control problem. In our numerical test case, the model reduced problem is defined upon POD.
To estimate the distance between ˜uand the unknown exactly optimal control ¯u, we proceed as follows: First, we determine the exact state ˜y:=yu˜by solving the state equation of (PL) with controlu:= ˜u. Next, we insert
˜
y foryu in the adjoint equation (1.8) to obtain the associated adjoint state ˜pas solution.
Now we are able to determine the perturbation ζ∈L2(Q) by (1.6) (takex:= (x, t) and D :=Q). Finally, we arrive at the estimate
ku˜−uk¯ L2(Q)≤ 1
λkζkL2(Q).
In [23, 26], this method of estimation was successfully applied to different optimal control problems with quadratic objective functional and linear state equation. It is also successfully applied for other reduced-basis approximations, see [24]. For recent extension to nonlinear problems we refer to [13], where the presented error estimates are utilized in a multilevel SQP algorithm. Notice that we tacitly assume to solve the state and adjoint equation for ˜yand ˜pexactly, i.e. we neglect associated discretization errors. The inclusion of such errors is subject of ongoing research.
2. The nonlinear case
2.1. The perturbation method for nonconvex functionals We consider now the nonconvex but smooth optimization problemmin
u∈Kf(u) := 1
2kG(u)−yHk2H+λ
2 kuk2L2(D), (P)
where G : L∞(D) → H is a twice continuously Fr´echet differentiable operator and all other quantities are defined as in (PQ). We assume that, for allu∈ K, the first- and second-order derivativesG0(u) :L∞(D)→H andG00(u) :L∞(D)×L∞(D)→H can be continuously extended toL2(D): There exists some constant c >0 not depending onuandv such that
kG0(u)vkH ≤ckvkL2(D) ∀u∈ K, ∀v∈L∞(D) (2.9)
kG00(u)(v1, v2)kH ≤ckv1kL2(D)kv2kL2(D) ∀u∈ K,∀v1, v2∈L∞(D) (2.10)
with some constant c > 0. Then the operators G0(u) and G00(u) can also be applied to increments v, v1, v2 in L2(D). Therefore, we can view G0(u) as continuous linear operator from L2(D) to H so that its adjoint operator maps continuouslyH toL2(D).
The derivativef0(u) is given by
f0(u)v = (G(u)−yH, G0(u)v)H+λ(u, v)L2(D)
= (G0(u)∗(G(u)−yH) +λu, v)L2(D). (2.11) TheL2(D)-function
pu:=G0(u)∗(G(u)−yH) (2.12)
is theadjoint stateassociated withu.
Let us determine for convenience also the second-order derivative of f : L∞(D) → R. Thanks to the assumptions onG, this derivative exists. To determine it, we consider the expression forf0with fixed increment v:=v1∈L∞(D) and differentiate again in the directionv2. We find by the chain and product rule
f00(u)(v1, v2) = (G0(u)v2, G0(u)v1)H
+ (G(u)−yH, G00(u)(v2, v1))H+λ(v2, v1)L2(D).
(2.13)
By our assumptions onG, alsof00(u) can be continuously extended to a bilinear form onL2(D)×L2(D) and
|f00(u)(v1, v2)| ≤ckv1kL2(D)kv2kL2(D) ∀u∈ K, ∀v1, v2∈L∞(D).
Let us now assume that ¯u∈ Kis a locally optimal solution to (P) in the sense of L∞(D), i.e. there is some ρ >0 such that
f(u)≥f(¯u) ∀u∈ K withku−uk¯ L∞(D)≤ρ.
We obtain the following extension of Theorem 1.1:
Theorem 2.1. Ifu¯∈ Kis a locally optimal solution of (P)in the sense ofL∞(D), then it obeys the variational inequality
Z
D
{(G0(¯u)∗(G(¯u)−yH))(x) +λu(x)}¯ (u(x)−u(x))¯ dx≥0 ∀ u∈ K. (2.14) This result follows from the variational inequality
f0(¯u)(u−u)¯ ≥0 ∀u∈ K, cf. e.g. Lemma 4.18 in [25]. Here, we use the representation (2.11) of f0(¯u).
Let now ¯p := pu¯ be the adjoint state defined in (2.12). Analogously to the quadratic problem (PQ), the variational inequality can be expressed in a pointwise form,
(¯p(x) +λ¯u(x))¯u(x) = min
v∈[α,β](¯p(x) +λ¯u(x))v for a.a. x∈D. (2.15) This is identical with (1.4), hence we can draw the same conclusions from (2.15) as in the linear case. Therefore, the perturbation ζis defined exactly as in (1.6).
Let now be ˜u ∈ K sufficiently close to ¯u. We pose the same question as for (PQ): Can we quantify the distance k˜u−uk? Now, however, the situation is more complicated. We need some second-order information¯ on ¯u. Assume that there exists someδ >0 such that the coercivity condition
f00(¯u)(v, v)≥δkvk2L2(D) ∀v∈L2(D) (2.16)
is satisfied. Then for any 0< δ0 < δthere exists a radiusr(δ0)>0 such that
f00(u)(v, v)≥δ0kvk2L2(D) ∀uwithku−uk¯ L∞(D)≤r(δ0),∀v∈L2(D). (2.17) If ˜ubelongs to this neighborhood, then we are able to estimate the distance.
Remark 2.2. To be on the safe side, we might select δ0 := δ/2 and set r := r(δ/2). This can be too pessimistic. In the application to POD we are mainly interested in the order of the error so that the factor 1/2 is not important. We useδ0:=δ, although this can be slightly too optimistic. ♦ Theorem 2.3. Letu¯ be locally optimal for (P) and assume thatu¯ satisfies the second-order condition (2.16).
If u˜∈ K is given such that k˜u−uk¯ L∞(D)≤r, where ris the radius introduced in (2.17) forδ0 :=δ/2, then it holds
ku˜−uk¯ L2(D)≤2
δkζkL2(D), whereζ is defined as in (1.6).
Proof. As in the convex quadratic case, ˜usatisfies the first-order necessary optimality conditions for the per- turbed optimization problem
minu∈Kf(u) + (ζ, u)L2(D). (Pζ)
We insert ¯uin the variational inequality for ˜uand ˜uin the one for ¯uand obtain (f0(˜u) +ζ,u¯−u)˜ L2(D)≥0
(f0(¯u),u˜−u)¯ L2(D)≥0.
We add both inequalities and find
(f0(˜u)−f0(¯u),u¯−u)˜ L2(D)+ (ζ,u¯−u)˜ L2(D)≥0.
The mean value theorem implies
−f00(ˆu)(¯u−u)˜ 2+ (ζ,u¯−u)˜ L2(D)≥0
with some ˆu∈ {v ∈L2(D)|v=s¯u+ (1−s)˜uwiths∈(0,1)}. Now we apply (2.17) to the left-hand side and the Cauchy-Schwarz inequality to the right-hand side to deduce
δ
2ku˜−uk¯ 2L2(D)≤ kζkL2(D)k˜u−uk¯ L2(D).
From this inequality, the assertion of the theorem follows in turn.
2.2. Application to the parabolic boundary control problem (PB)
Our general result can be applied in various ways. Let us explain its use for the parabolic boundary control problem (PB). From the theoretical point of view, (PB) is more difficult than the distributed problem (PD). For (PB), the control-to-state mappingGis not differentiable fromL2(Σ) toC( ¯Q). Here, the two-norm discrepancy is a genuine issue. In (PD), this problem does not occur for the spatial dimensionn= 1, while it also arises forn≥2. Therefore, we selected (PB) for our analysis.
We define
Uad={u∈L∞(Σ) :α≤u(x, t)≤β for a.a. (x, t)∈Σ}. (2.18) The following result is known on the solvability of the state equation:
Theorem 2.4. For allu∈Uad, the state equation (0.2)has a unique solutionyu∈Y :=W(0, T)∩C( ¯Q). The mappingG:u7→yuis twice continuously differentiable fromL∞(Σ)toY. In particular, it is twice continuously differentiable fromL∞(Σ)toH :=L2(Q). The first- and second-order derivativesG0andG00 obey the extension conditions (2.9),(2.10) for the choice D:=Q,K=Uad.
We refer to [25], Theorems 5.9 and 5.16. The extension conditions follow immediately from the representation formulas for these derivatives stated in these theorems. The same theorems include also the result for the problem (PD).
Here, the reduced functional is given by f(u) := 1
2 Z T
0
Z
Ω
(yu(x, t)−yd(x, t))2dxdt+λ 2
Z T 0
Z
Γ
u(x, t)2ds(x)dt.
Analogously to (2.11) and (2.12), the first derivative of the reduced functional can be expressed in the form f0(u)v=
Z T 0
Z
Γ
(pu(x, t) +λu(x, t))v(x, t)ds(x)dt, (2.19) wherepu is the adjoint state associated withu. It is the unique solution to the adjoint equation
−∂p
∂t −∆p = yu−yd in Ω×(0, T)
∂p
∂ν +∂b
∂y(x, t, yu)p = 0 in Γ×(0, T)
p(x, T) = 0 in Ω.
(2.20)
In view of our Theorem 2.1, any locally optimal control ¯uof (PB) has to obey the following variational inequality:
Z T 0
Z
Γ
(¯p(x, t) +λ¯u(x, t))(u(x, t)−¯u(x, t))ds(x)dt≥0 ∀u∈Uad, (2.21) where ¯p:=pu¯.
The general form of the second derivative of the reduced functional was determined in (2.13). In this representation, terms of the formG0(u)vappeared. For givenv∈L2(Σ), the functiony=G0(u)vis the unique solution of thelinearized parabolic equation
∂y
∂t −∆y = 0 in Ω×(0, T)
∂y
∂ν +∂b
∂y(x, t, yu)y = v in Γ×(0, T) y(x,0) = 0 in Ω,
(2.22)
where yu = G(u) is the state associated with u, cf. [25, Theorem 5.9]. Therefore, the mapping v 7→ y is continuous fromL2(Σ) toW(0, T), in particular toL2(Q); this yields the first condition of continuous extension (2.9).
The second derivativez =G00(u)(v1, v2) is the unique solution to the same equation withv substituted by v:=−∂2b
∂y2(x, t, yu)y1y2, whereyi=G0(u)vi,i= 1,2, and
∂z
∂t −∆z = 0 in Ω×(0, T)
∂z
∂ν +∂b
∂z(x, t, yu)z = −∂2b
∂y2(x, t, yu)y1y2 in Γ×(0, T)
z(x,0) = 0 in Ω.
(2.23)
This is a linear equation with ”control”−∂2b
∂y2(x, t, yu)y1y2. Therefore, it is not difficult to see that the second term in the representation (2.13) forf00(u) can be re-written in terms of the adjoint state. Namely, it holds
f00(u)v2 = J00(yu, u)(y, v)2− Z T
0
Z
Ω
pu ∂2b
∂y2(x, t, yu)y2ds(x)dt,
= Z T
0
Z
Ω
y2dxdt+λ Z T
0
Z
Γ
v2ds(x)dt− Z T
0
Z
Γ
pu
∂2b
∂y2(x, t, yu)y2ds(x)dt.
(2.24)
We refer, for instance to [25, Section 5.7] or to the recent paper [5], where second-order derivatives are determined for more general quasilinear equations. The second extension condition (2.10) follows immediately from (2.23).
In (2.24),y is the solution to the linearized equation (2.22). The expression above is just the second-order derivative of the Lagrangian function defined upon the adjoint state pu associated with u. Therefore, the second-order sufficient optimality condition (2.16) means that
Z T 0
Z
Ω
y2dxdt+ Z T
0
Z
Γ
λ v2−pu ∂2b
∂y2(x, t, yu)y2
ds(x)dt
≥δ Z T
0
Z
Γ
v2ds(x)dt ∀v∈L2(Σ).
(2.25)
Sinceyis the solution of the linearized equation that is associated tovand the left-hand side of this equation is the second derivative of the Lagrangian function, this expresses a well known formulation of the (strong) second-order condition: The second derivative of the Lagrangian function is assumed to be coercive on the subspace defined by the linearized equation.
2.3. Application of the perturbation method
Let now ˜u∈Uadbe an arbitrary control that is close to a locally optimal control ¯uof (PB). We assume that
¯
usatisfies the sufficient second-order optimality condition with someδ >0. In other words, the coercivity con- dition (2.25) is fulfilled foru:= ¯u. Since the second-order condition is stable with respect toL∞-perturbations of ¯u, we are justified to assume that there exists someρ >0 such that
f00(u)v2≥ δ 2
Z T 0
Z
Γ
v2ds(x)dt ∀v∈L2(Σ),∀u∈Bρ(¯u). (2.26) We assume that ˜ubelongs toBρ(¯u). Now we define the perturbation functionζ as in (1.6), (x, t)∈Ω×[0, T] is substituted forx. Thanks to Theorem 2.3, we obtain
k˜u−uk¯ 2L2(Σ)≤ 2
δkζk2L2(Σ). (2.27)
We should mention already here a serious theoretical obstacle that can hardly be rigorously overcome. To apply our method of a posteriori estimation, we need the numbers ρand δ. As ρis concerned, we can only assume that the method of determining the (suboptimal) control ˜uwas sufficiently precise so thatku˜−uk¯ < ρ.
There exists a method by R¨osch and Wachsmuth [20], [29] to verify if there exists a locally optimal control
¯
uin a certain neighborhood of a given ˜u. However, this tool is difficult to apply. It will work only for fairly special problems.
In addition, we must determine the coercivity constantδ. In computations, we deal with a finite-dimensional approximation of the control problem and determine the smallest eigenvalue of the reduced Hessian. In general, this way is not reliable in estimating the coercivity constantδfor the infinite-dimensional undiscretized optimal control problem: In [29], a discouraging example is constructed, where computations with very small mesh size
indicate a sufficient second-order condition for a point, which is a saddle point. In special cases, the ideas of [20]
and [29] might be applied that rely on deep analytical estimations.
Summarizing this issue, we must confess that there is no reliable and at the same time practicable method to verify our assumptions absolutely certain. We have to trust that our problem behaves well in a neighborhood of the unknown solution ¯u. This is similar to problems of nonlinear programming, where optimization routines are started without evidence whether a constraint qualification is satisfied at the unknown solution or not.
We assume in the sequel that our computed suboptimal control belongs to a neighborhood of a locally optimal solution ¯uthat satisfies a second-order sufficient optimality condition. Moreover, we assume that we are able to determine the order of the coercivity constant δ by the lowest eigenvalue of the Hessian matrix associated with the suboptimal control.
3. Model reduction by POD
In this section, we give a very brief survey on how to establish a reduced order model of the parabolic boundary control problem (PB) by applying standard POD. For more information and proofs we refer the reader to Kunisch and Volkwein [14] or Volkwein [28], for instance. A survey on methods of model reduction and a comparison of different reduction methods is contained in Antoulas [1].
3.1. The discrete POD method
Define H = L2(Ω) and let u ∈ L2(Σ) be an arbitrary control with associated state yu = G(u) ∈ Y ⊆ L2(0, T;H). Since we cannot compute the whole trajectoryy(t) for allt∈[0, T], we define a partition 0 =t0<
t1< . . . < tn=T of the time interval [0, T] and compute shapshotsyi:=y(ti)∈H, i= 1, . . . , n,of the statey at the given time instancest1, . . . , tn.
Our goal is to find a small Galerkin basis that well expresses the main properties of the PDE in the finite dimensional subspaceHn defined by
Hn= span{y(ti)|i= 0, . . . , n} ⊂H
with dimension d= dimHn. Let{φ1, . . . , φd}be a basis of Hn that is orthonormal with respect to the inner product ofL2(Ω). Then everyyi=y(ti), i= 0, . . . , n,is a linear combination ofφ1, . . . , φd,
yi=
d
X
j=1
(yi, φj)Hφj.
Now choose a fixed number `∈ {1, . . . , d}. In the application to model reduction,` is small compared withd.
We intend to find an orthonormal set of functions{ϕj}`j=1⊂H such that the sum of squared errors
εi= yi−
`
X
j=1
(yi, ϕj)Hϕj
2 H
is minimized. In other words, we consider the optimization problem
min
ϕ1,...,ϕ`∈H n
X
i=0
αi
yi−
`
X
j=1
(yi, ϕj)Hϕj
2
H =:Jn(ϕ1, . . . , ϕ`), s.t. (ϕi, ϕj)H=δij ∀i, j∈ {1, . . . , `},
(3.28)
whereα0, . . . , αn are the trapezoidal weights
α0= (t1−t0)/2, αi= (ti+1−ti−1)/2, for i= 1, . . . , n−1, αn= (tn−tn−1)/2.
The solution to (3.28) is obtained by solving a certain eigenvalue problem, described by the linear operator Rn:H → Hn that maps z∈H to
Rnz=
n
X
i=0
αi(z, yi)Hyi.
The following proposition characterizes the operator Rn and shows that the eigenfunctions to the ` largest eigenvalues ofRn solve the problem (3.28).
Proposition 3.1. Let H be a Hilbert space, y ∈ C([0, T], H) and Hn = span {y(ti)|i= 0, . . . , n} ⊂H with dimension d=dim Hn. Then the following statements hold:
(i) The operator Rn:H → Hn is bounded, self-adjoint, compact and non-negative.
(ii) There exist a sequence of eigenvalues{λni}di=1withλn1 ≥λn2 ≥. . .≥λnd ≥0and associated eigenelements {ψin}di=1⊂H that solve the eigenvalue problem
Rnψin=λniψin, i= 1, . . . , d. (3.29) (iii) For any fixed number ` ∈ {1, . . . , d}, the minimization problem (3.28) is solved by the eigenfunctions
{ψ1n, . . . , ψn`}.
(iv) The corresponding minimal value ofJn is given by inf
ϕ1,...,ϕ`∈HJn(ϕ1, . . . , ϕ`) =Jn(ψ1n, . . . , ψn`) =
d
X
i=`+1
λni.
For a proof we refer to e.g. Volkwein [28], [27, Section 1.3] and Holmes et al. [11, Section 3].
Note that we haveRn=YY∗, where the bounded linear operatorY:Rn→ Hn and its Hilbert space adjoint operatorY∗:H →Rn are defined by
Yϕ=
n
X
i=0
αiϕiyi, (Y∗z)i= (z, yi)H, i= 1, . . . , n.
Analogously to the theory of singular value decomposition for matrices, the linear, bounded, compact, and self-adjoint operatorKn=Y∗Y:Rn→Rn defined by
(Knϕ)i=
n
X
j=1
αj(yj, yi)Hϕj, i= 1, . . . , n,
has the same eigenvalues{λni}di=1 asRn. The corresponding eigenfunctions ofKn are obtained from theψni by transformations
(Φni)j= 1
pλni (Y∗ψni)j = 1
pλni (ψin, yj)H, j∈ {1, . . . , n},
for all i ∈ {1, . . . , `}. In many cases, the number n of time instances is much smaller than the number d of spatial ansatz functions. Then it is convenient to solve the eigenvalue problem
KnΦni =λniΦni, i= 1, . . . , d, (3.30) instead of using the operatorRn.
Now we apply a standard finite element approximation with respect to the spatial variable x. For this purpose, letThbe a family of regular triangulations of Ω with mesh sizeh >0. For eachTh, letVh⊂L2(Ω) be the set of piecewise linear finite element basis functions{ϕ1, . . . , ϕm}, m∈N.
Then for every functionz∈Vhthere existz1, . . . , zm∈Rsuch that z(x) =
m
X
j=1
zjϕj(x), for allx∈Ω. (3.31)
Thanks to our numerical approximation of the PDE, we can assume that our snapshots yi belong to Vh. Therefore, for everyi∈ {1, . . . , n},
yi=
m
X
j=1
yjiϕj.
For convenience, we introduce the vectors ~yi = (yi1, . . . , ymi )>. In this way, we have the correspondenceVh 3 yi ∼~yi∈Rm. Letz∈Vh be given. Then
(Rnz)(x) =
n
X
i=0
αi(yi, z)L2(Ω)yi(x).
By definition, we have Rn : H → Hn ⊂Vh. Let us denote for convenience the restriction of Rn to Vh with range in Vh by the same symbol so thatRn:Vh→Vh. Applying equation (3.31) for every nodex1, . . . , xm of the triangulation of Ω, we infer forVh3z∼~z∈RmandVh3yi∼~yi∈Rm
(Rnz)(xj) =
n
X
i=0 m
X
k=1
αiyikϕk(xj)
m
X
l,ν=1
yilzν Z
Ω
ϕl(x)ϕν(x)dx
=
n
X
i=0 m
X
l,ν=1
αiyjiyilzν Z
Ω
ϕl(x)ϕν(x)dx
=
n
X
i=0 m
X
l,ν=1
αiYjiYli(Mh)lνzν =
n
X
i=0
αiYji
m
X
l,ν=1
Yil>(Mh)lνzν
=
n
X
i=0
αiYji Y>Mh~z = Y D Y>Mh~z
j,
where ~yi = (y1i, . . . , yim)>, i = 1, . . . , n, Y = [~y1|. . .|~yn], D = diag(α1, . . . , αn). In the first line, we used ϕk(xj) =δkj. Moreover,Mh with entries
(Mh)k,l= Z
Ω
ϕk(x)ϕl(x)dx
is the mass matrix associated with the finite element basis functions. Define ¯Rn :Rm→Rm by R¯n~z:=
n
X
i=0
αi(~yi>Mh~z)~yi. Since
(Rnz)(xj) =
n
X
i=0
αi(yi, z)L2(Ω)yi(xj) =
n
X
i=0
αi~y>i Mh~z yi(xj)
=
n
X
i=0
αi(~y>i Mh~z)yji =
n
X
i=0
αi(~yi>Mh~z)yi
!
j
= R¯n~z
j,
the operatorRn corresponds to the matrix ¯Rn in the numerical approximation, i.e. we have R¯n~z=Y D Y>Mh~z.
To obtain an optimal solution of (3.28), we solve the eigenvalue problem
R¯nψ~in=Y DY>Mhψ~in=λniψ~in, i= 1, . . . , `. (3.32) Defining ˆY =M
1 2
hY D12, ˆψin=M
1 2
hψ~ni and multiplying (3.33) byM
1 2
h from the left, we arrive at the symmetric eigenvalue problem
YˆYˆ>ψˆni =λniψˆni, i= 1, . . . , `.
Using similar arguments for the operatorK, we also can solve the discretized eigenvalue problem
K¯nΦ~ni =Y>MhY D ~Φni =λniΦ~ni, i= 1, . . . , `. (3.33) Defining ˆΦni =D12~Φni and multiplying (3.33) byD12 from the left, we infer
Yˆ>YˆΦˆni =λniΦˆni, i= 1, . . . , `.
Summarizing, we compute thePOD basis of rank `by themethod of snapshots as follows:
Algorithm 3.2 (Method of Snapshots). Let mn.
(1) Solve the symmetric n×neigenvalue problem
Yˆ>YˆΦˆni =λniΦˆni, i= 1, . . . , `.
(2) Transform the eigenvectorsΦˆn1, . . . ,Φˆn` by singular value decomposition to the POD basis vectorsψ~n1, . . . , ~ψ`n ψ~in=M−
1 2
h ψˆin= 1 pλni M−
1 2
h YˆΦˆni = 1
pλni Y D12Φˆni, i∈ {1, . . . , `}.
Then, the POD basis functionsVh3ψni ∼ψ~in∈Rm, i= 1, . . . , `,are given by
ψin=
m
X
j=1
(ψ~ni)jϕj(x).
Remark 3.3. (1) In the method of snapshots, we do not have to set up the matrixM
1 2
h, since Yˆ>Yˆ =D12Y>MhY D12.
(2) Note that (3.32) is an m×m eigenvalue problem, whereas (3.33) has dimension n×n. In many applications we have m n. Hence, the operator Kn is preferred for generating the POD basis. In the case, m ≤n, one should solve the problem (3.32) instead. On the other hand, the singular value decomposition is much more stable in computing the small eigenvalues, see e.g. [15, 23] for this aspect
in the context of POD. ♦
3.2. POD Galerkin projection
Let us consider exemplarily the state equation (0.2) of our boundary control problem (PB). As above, we denote byVh the space of piecewise linear functions spanned by the FE basis functions{ϕ1, . . . , ϕm} defined on a regular triangulationTh of Ω. For convenience, we use the short notation (·,·)Ω:= (·,·)L2(Ω)and (·,·)Γ:=
(·,·)L2(Γ)for inner products. Then the variational formulation of the equations (0.2) is given by d
dt(y(t), ϕ)Ω+ (∇y(t),∇ϕ)Ω+ (b(·, t, y(t)), ϕ)Γ = (u(t), ϕ)Γ
(y(0), ϕ)Ω = (y0, ϕ)Ω,
(3.34)
for all ϕ ∈ H1(Ω) and for almost all t ∈ [0, T]. Under our assumptions on the given data, to each u ∈ Lr(Γ×(0, T)) withr >2, there exists a unique solutiony ∈W(0, T)∩C( ¯Q) of (0.2).
Let now {ti}ni=0 be an equidistant partition of [0, T]. To generate the snapshots we select a fixed reference control us ∈ Lr(Γ) and solve equation (3.34) at the given time instances ti. We denote the snapshots by ysi =ys(ti), i= 0, . . . , nand, as before,
Hn= span{ys(ti)|i∈ {0, . . . , n}}
with dimension d ∈ N. Next, we fix a small natural number ` ≤ d and compute the POD basis functions {ψ1n, . . . , ψ`n} as described in Section 3.1. Notice that all ψin, i= 1, . . . , `, are linear combinations of the FE basis functionsϕ1, . . . , ϕm.
Applying a Galerkin scheme based on the (small) POD basis, we obtain a nonlinear initial value problem of finding a functiony` with
y`(t) =
`
X
i=1
ηi(t)ψin, such that
d
dt(y`(t), ψin)Ω+ (∇y`(t),∇ψin)Ω+ (b(·, t, y`(t)), ψni)Γ = (u(t), ψni)Γ
(y`(0), ψni)Ω = (y0, ψin)Ω
(3.35) is satisfied for almost allt∈[0, T] and alli∈ {1, . . . , `}. This is the state equation of a low size optimal control problem with state space of dimension `, see problem (PB.1`) below. This small optimal control problem is then solved to obtain a suboptimal control ¯u`.
Remark 3.4. It is not obvious that the reduced equation (3.35) has a unique solution, because the reduced Galerkin basis might destroy the monotonicity of the nonlinearity. To avoid this theoretical difficulty, we can truncatebas follows:
bc(x, t, y) =
b(x, t, c) if y≥c b(x, t, y) if |y| ≤c/2 b(x, t,−c) if y≤ −c,
wherebc is defined in (−c,−c/2)∪(c/2, c) such thatbc isC2 with respect toy. Thenbc is uniformly Lipschitz with respect to y so that the theorem by Picard and Lindel¨of permits to show existence and uniqueness of a solution of (3.35) for all T >0. In our computations, the functiony` was uniformly bounded so that we did
not follow this idea. ♦
4. A posteriori error estimation for POD solutions
To apply our a posteriori error estimation technique described in Section 2, we select a solution of some POD reduced optimal control problem. This is convenient, since model reduction by POD does not provide a rigorous a priori error analysis. Nevertheless, it often gives excellent results. This will be illustrated in the following two numerical examples.
4.1. Example 1
We first discuss a one-dimensional boundary control problem that is governed by a semilinear parabolic equation. The boundary control u is to drive a spatio-temporal temperature distribution to a predefined, desired temperature distribution at the final timeT.
Let Ω = (0,1), T = 1.58 andQ:= Ω×(0, T). The optimal control problem is given by min1
2 Z 1
0
y(x, T)−1
2(1−x2)2 dx+ 1
20 Z T
0
u(t)2dt=:J(y, u) (PB.1) subject to the heat equation with Stefan-Boltzmann type boundary condition
yt(x, t)−yxx(x, t) = 0 inQ yx(0, t) = 0 in (0, T) yx(1, t) +y(1, t) +y4(1, t) = u(t) in (0, T)
y(x,0) = 0 in Ω
(4.36)
and to the bilateral control constraints
−1≤u(t)≤1, t∈(0, T).
This is a well-known test example that goes back to Schittkowski [22].
Remark 4.1. Formally, the functiony7→y4 does not fulfill our assumption of monotonicity. Therefore, often y 7→ |y|y3 is considered instead. In our numerical tests, all occuring state functions y were non-negative.
Therefore, we are justified to writey4. ♦
The corresponding adjoint equation is defined by
−pt(x, t)−∆p(x, t) = 0 px(0, t) = 0 px(1, t) +p(1, t) + 4y3(1, t)p(1, t) = 0,
p(x, T) = y(x, T)−yd(x).
(4.37)
To reduce the state equation by POD, we selected an equidistant partition of the time interval, 0 =t0 <
t1 < . . . < tn = T, and a piecewise linear finite element discretization on the interval Ω with standard FE basis (”hat”) functions{ϕ1, . . . , ϕm}. A semi-implicit Euler scheme in time was then used to solve the resulting semidiscrete equation (3.35). Associated with the fixed controlu≡0.5, the shapshots
(
yj=y(tj) =
m
X
i=1
yijϕi
)n
j=0
⊂Vh⊂L2(Ω)
were generated. Now we collect the snapshot vectors ~yj = (yj1, . . . , ymj )>, j = 0, . . . , n, in the matrix Y = [~y0|. . .|~yn]∈Rm×n. As above,M denotes the mass matrix associated with {ϕ1, . . . , ϕm}with entries
Mij= (ϕi, ϕj)L2(Ω), i, j∈ {1, . . . , m}.
For convenience, we ignore the dependence onh. Furthermore, let ∆t=T /nandα0, . . . , αn be the trapezoidal weights
α0=αn= ∆t/2, αi= ∆t, for i= 1, . . . , n−1.
DefineD= diag(α0, . . . , αn). Then the small orthonormal system{ψ1, . . . , ψ`} is obtained from the eigenvalue decomposition ofY>M Y D, see Section 3.1, Algorithm 3.2.
We use the POD Galerkin ansatz for both the stateyfrom (4.36) and the corresponding adjoint statepthat solves the adjoint equation (4.37),
y`(x, t) =
`
X
i=1
yi(t)ψi(x), p`(x, t) =
`
X
i=1
pi(t)ψi(x).
The associated POD reduced optimal control problem is min1
2 Z 1
0
y`(x, T)−1
2(1−x2)2 dx+ 1
20 Z T
0
u(t)2dt (PB.1`)
subject to
M`∂
∂t~y(t) + (K`+F`)~y(t) +G`(~y(t)) = B`u(t)
~y(0) = 0,
(4.38) and
−1≤u(t)≤1,
for almost all t ∈ (0, T). Here, ~y : [0, T] → R` denotes the vector function ~y(t) = (y1(t), . . . , y`(t))>. The matricesM`, K`, F`∈R`×`,B`∈R` are given by
Mij` = (ψi, ψj)L2(Ω), Kij` = (ψ0i, ψ0j)L2(Ω), Fij` =ψi(1)ψj(1), Bi`=ψi(1), andG`:R`→R`is defined by y7→(y>B`)4B`.
The adjoint state~p: [0, T]→R` associated with~y is the unique solution of the equation
−M`∂
∂t~p(t) + (K`+F`)~p(t) + 4H`(~y(t))F`~p(t) = 0
~
p(T) = ~y(T)−~yd.
(4.39)
withH`:R`→R`given byy7→(y>B`)3B`.
The state equation of the associated discretized optimal control problem is solved again by a semi-implicit Euler scheme and the control functionsuare chosen as step function according to the given partition{t0, . . . , tn} by
u(t) =ui, ti−1≤t < ti, i= 1, . . . , n.
Therefore, we have the equivalenceu(·)∼(u1, . . . , un)>=~u>∈Rn.
This discretized optimal control problem was solved by an SQP method, where the associated sequence of linear-quadratic control problems systems was treated by a primal-dual active set strategy. Let us denote by
~¯
u` = (¯u`1, . . . ,u¯`n)> the obtained locally suboptimal control vector, and by ¯u`τ(·) ∼~u¯` the corresponding step function.
We are interested in estimating the distance between ¯u`τ and the next locally optimal control ¯u. For this purpose, we apply our perturbation method that was formulated in Section 2.3. This needs two full PDE solves.
First, we have to determine the state ¯y` =G(¯u`τ), i.e. the solution of the parabolic boundary value problem (4.36) associated with ¯u`τ. Second, we must compute the solution ¯p` of the full adjoint equation (4.37). Then the time-dependent perturbationζ` is also a step function on [0, T]. It is equivalent to the vector (ζ1, . . . , ζn)>
