3.4 Optimality System Proper Orthogonal Decomposition
3.4.5 Iterative methods
A direct numerical solving of the coupled nonlinear OS-POD optimality system is not recommandable since the numerical costs would exeed the effort of the original optimal control problem without application of any model reduction techniques. In the following we present a local convergence result for an iterative fixpoint method which allows a quite flexible splitting of the coupled equations. The price one has to pay is a restric-tive assumption on the magnitude of the regularization parameters σu, σw as well as a relatively slow convergence rate of the iteration.
1. For a given control-state pair (u, w) ∈ U ×W let y = y(u) ∈ Y denote the state solution to
˙
y(t) +Ay(t) =f(t) +Bu(t), y(0) =y◦. (3.64a) 2. With(ψ,λ) = (ψ,λ)(y)∈V`×R` we denote signified eigenfunction-eigenvalue pairs
to the operatorR(y),
R(y)ψl=λlψl, kψlkX = 1 (3.64b) which shall be specified below.
3. Given the control u and the POD basis ψ, y(u,ψ) ∈ H1(Θ,R`) is chosen as the reduced state solution to
M(ψ) ˙y(t) + A(ψ)y(t) = f(ψ;t) + B(ψ)u(t), M(ψ)y(0) = y◦(ψ). (3.64c) 4. Withw,y,ψ we define the adjoint statep∈H1(Θ,R`)as the solution to the reduced
equation
−M(ψ) ˙p(t) + A(ψ)p(t) =−σQ(M(ψ)y(t)−yQ(ψ;t)) + σεwIT(ψ)w(t) (3.64d) M(ψ)p(T) =−σΩ(M(ψ)y(T)−yΩ(ψ)).
5. We proceed with the determination of the multiplierµ=µ(y,ψ,λ,y, w)∈V`: (R(y)−λl)µl=hGl, ψliX0,Xψl−i−1X Gl. (3.64e) 6. Now the adjoint statep=p(y,ψ,µ)∈Y is given by
−p(t) +˙ Ap(t) =−
`
X
l=1
hy(t), µliXiXψl−
`
X
l=1
hy(t), ψliXiXµl, (3.64f) p(T) = 0.
We define the selfmappingΦon U×W by the remaining two optimality conditions:
Φ(u, w) = 1
σu
B?p+ BT(ψ)p−κu
,− ε
σwκw
(3.65)
where the multipliersκu =κu(ψ, p,p)∈U andκw =κw(y,ψ)∈W are given by κu = max(0,B?p+ BT(ψ)p−σuub) + min(0,B?p+ BT(ψ)p−σuua), (3.66a) κw = max
0,σw
ε (I(ψ)y−yb)
+ min
0,σw
ε (I(ψ)y−ya)
. (3.66b)
Fixpoint iteration.
Letx¯∈Xbe a solution to the OS-POD problem (3.13) with multipliersz¯∈Z andκ¯ ∈Z˜ satisfying the dual system (3.63). We show thatΦis a contraction in a neighborhood of (¯u,w)¯ if the regularization parameters σu, σw are chosen sufficiently large and the data is small enough to get the following result:
Theorem 3.15. (Banach-Iteration)
For any parametersσu, σw >0, let(¯u,w)¯ ∈U×W be an optimal control-penalty pair of (3.44).
1. ThenΦ is locally Lipschitz continuous at(¯u,w).¯
2. Especially, ifσu, σw are sufficiently large, Φis a contraction.
3. Let (u0, w0) ∈ U ×W be sufficiently close to (¯u,w), then the Banach iteration¯ (un+1, wn+1) = Φ(un, wn)converges towards (¯u,w).¯
4. In this case, the following locally a-priori convergence rate is ensured:
k(un, wn)−(¯u,w)k¯ U×W ≤ Cn
1−Ck(u1, w1)−(u0, w0)kU×W
whereC∈(0,1)denotes the contraction constant.
Proof. Let >0and let(˜u,w)˜ ∈Xad withk¯u−uk˜ U < andkw¯−wk˜ W < . Then kΦu(¯u,w)¯ −Φu(˜u,w)k˜ U =σu−1kB?(¯p−p) + B˜ T( ¯ψ)(¯p−˜p)
+ (BT( ¯ψ)−BT( ˜ψ))˜p + ( ¯κu−κ˜u)kU
≤Cσ−1u kBk˜ Lb(Rm,V0)(kp¯−pk˜ L2(Θ,V)
+kψk¯ V`k¯p−pk˜ L2(Θ,R`)+k˜pkL2(Θ,R`)kψ¯ −ψk˜ V`). (3.67) According to the a-priori energy estimate (1.18),k˜pkL2(Θ,R`) is uniformly bounded:
k˜pk2H1(Θ,R`)≤C(k˜yk2L2(Θ,R`)+k˜y(T)k2
R`+kyQk2L2(Θ,H)
+kyΩk2H +σwε−1kIk˜ Lb(V,Rn)kwk˜ 2W);
since k˜y(T)k2
R` ≤ k˜ykC0(Θ,R`)≤Ck˜ykH1(Θ,R`), we conclude k˜pk2H1(Θ,
R`)≤C((k¯uk2U+2) +kfk2L2(Θ,V0)+ky◦k2H + (kwk¯ 2W +2) +kyQk2L2(Θ,H)+kyΩk2H).
We require here thatis sufficiently small (3.45a), (3.45b) so that the uniform coercivity and continuity conditions (3.49a), (3.49b) for the reduced bilinear form are satisfied.
Then, since u¯= Φu(¯u,w)¯ holds,
kΦu(˜u,w)˜ −uk¯ U ≤Cσu−1(kp¯−pk˜ L2(Θ,V)+k¯p−pk˜ L2(Θ,R`)+kψ¯ −ψk˜ V`). (3.68) 1. We interprete ψ˜l,λ˜l as a perturbation of ψ¯l,¯λl. According to the stability result for
eigenvalue problems in [34], Thm. 1.7, we have kψ¯l−ψ˜lkX ≤Ck(R(¯y)−¯λl)|−1{¯
ψl}⊥kLb(X,X)kR(¯y)− R(˜y)kLb(X,X)kψ¯lkX, (3.69a)
|λ¯l−λ˜l| ≤CkR(¯y)− R(˜y)kLb(X,X)kψ¯lk2X. (3.69b) 2. To estimate the residual of the eigenoperators, letϕ∈V withkϕkX = 1, then
k(R(¯y)− R(˜y))ϕkV = Z
Θ
hy(t), ϕi¯ Xy(t) dt¯ − Z
Θ
h˜y(t), ϕiXy(t) dt˜ X
≤ Z
Θ
k¯y(t)−y(t)k˜ XkϕkX(ky(t)k¯ X +k˜y(t)kX) dt
≤Ck¯y−yk˜ L2(Θ,V) (3.69c)
since k˜ykL2(Θ,V) is uniformly bounded.
3. We proceed with an estimation for kψ¯ −ψk˜ V` in the case X = H. For sufficiently small,λ˜` is uniformly bounded from below by some ˜ >0, so
kψ¯l−ψ˜lkV =kλ¯−1l R(¯y) ¯ψl−λ˜−1l R(˜y) ˜ψlkV
≤ |¯λ−1l −λ˜−1l | kR(˜y) ˜ψlkV +|λ˜−1l | kR(¯y)− R(˜y)kLb(X,V)kψ¯lkV +|λ¯−1l | kR(¯y) ˜ψlkV
≤C((¯λl(¯λl−˜))−1kψ˜lkXk˜ykL2(Θ,X)kyk˜ L2(Θ,V)|λ¯l−λ˜l| + (¯λl−˜)(k¯ykL2(Θ,X)k¯y−yk˜ L2(Θ,V)
+k˜ykL2(Θ,V)k¯y−yk˜ L2(Θ,X))kψ¯lkV + ¯λ−1l k¯ykL2(Θ,X)k¯ykL2(Θ,V)kψ¯l−ψ˜lkX)
≤C(k¯λl−˜λl|+k¯y−yk˜ L2(Θ,V)+kψ¯l−ψ˜lkX). (3.69d) 4. Definep = ˜p−p, then¯ p solves the ODE system
−M( ˜ψ) ˙p(t) + A( ˜ψ)p(t) =σQ(−M( ˜ψ)(˜y(t)−y(t)) + (y¯ Q( ˜ψ;t)−yQ( ¯ψ;t)) + (M( ¯ψ)−M( ˜ψ)) ˙¯p(t) + (A( ˜ψ)−A( ¯ψ))¯p(t)) +σwε−1((IT( ˜ψ)−IT( ¯ψ)) ¯w(t) + IT( ˜ψ) ˜w(t)−w(t))¯ M( ˜ψ)p(T) =σΩ(−M( ˜ψ)(˜y(T)−¯y(T)) + (yΩ( ˜ψ)−yΩ( ¯ψ)))
+ (M( ¯ψ)−M( ˜ψ))¯p(T).
We get the a-priori stability estimation
k¯p−pk˜ H1(Θ,R`)≤C(kM( ˜ψ)kR`×`k˜y−¯ykH1(Θ,R`)
+kyQ( ˜ψ)−yQ( ¯ψ)k
R`+kyΩ( ˜ψ)−yΩ( ¯ψ)k
R`
+kM( ¯ψ)−M( ˜ψ)k
R`×`k¯pkH1(Θ,R`)
+kA( ¯ψ)−A( ˜ψ)k
R`×`k¯pkL2(Θ,R`)
+kI( ¯ψ)−I( ˜ψ)kR`×`kwk¯ L2(Θ,Rn)
+kI( ˜ψ)k
R`×`kw¯−wk˜ L2(Θ,Rn))
≤C(k˜y−yk¯ H1(Θ,R`)+kψ˜ −ψk¯ V`+kw¯−wk˜ L2(Θ,Rn)). (3.69e) 5. Further,y = ˜y−¯ysolves the ODE system
M( ˜ψ) ˙y(t) + A( ˜ψ)y(t) = B( ˜ψ)(˜u(t)−u(t)) + (B( ˜¯ ψ)−B( ¯ψ))¯u(t) + (M( ˜ψ)−M( ¯ψ)) ˙¯y(t) + (A( ˜ψ)−M( ¯ψ))¯y(t) + (f( ˜ψ;t)−f( ¯ψ;t))
M( ˜ψ)y(0) = (y◦( ˜ψ)−y◦( ¯ψ)) + (M( ¯ψ)−M( ˜ψ))¯y(0);
we receive the a-priori estimate
k¯y−˜ykH1(Θ,R`)≤C(kB( ˜ψ)k
R`×`kk˜u−uk¯ L2(Θ,Rm)
+kB( ˜ψ)−B( ¯ψ)k
R`×`k¯ukL2(Θ,Rm)
+kM( ¯ψ)−M( ˜ψ)k
R`×`k¯ykH1(Θ,R`)
+kA( ¯ψ)−A( ˜ψ)k
R`×`k¯ykL2(Θ,R`)
+ky◦( ¯ψ)−y◦( ˜ψ)kR`
+kf( ˜ψ)−f( ¯ψ)kL2(Θ,R`))
≤C(k˜u−uk¯ L2(Θ,Rm)+kψ˜ −ψk¯ V`). (3.69f) 6. Definep= ˜p−p, then¯ psolves the PDE
−p(t) +˙ Ap(t) =−
`
X
l=1
(hy(t),˜ ψ˜liXiXµ˜l− h¯y(t),ψ¯liXiXµ¯l)
−
`
X
l=1
(hy(t),˜ µ˜liXiXψ˜l− h¯y(t),µ¯liXiXψ¯l), p(T) = 0.
We get the a-priori stability estimation
k¯p−pk˜ L2(Θ,V)∩H1(Θ,V0) ≤C(ky˜−yk¯ L2(Θ,X)+kψ˜ −ψk¯ X`+kµ˜−µk¯ X`). (3.69g)
7. Definey= ˜y−y, then¯ y solves the PDE
˙
y(t) +Ay(t) = ˜B(˜u(t)−u(t)),¯ y(0) = 0;
we get the a-priori estimate
k˜y−yk¯ L2(Θ,V)∩H1(Θ,V0)≤Ck˜u−uk¯ L2(Θ,Rm). (3.69h) 8. Assume thatkBkL
b(L2(Θ,Rm),V0),kfkL2(Θ,V0) andky◦kH are sufficiently small so that kR(˜y)kLb(X,X) =k˜yk2L2(Θ,X)
≤C(kfkL2(Θ,V0)2+ky◦k2H +kBkLb(L2(Θ,Rm),V0)(k¯uk2L2(Θ,Rm)+2))<C˜ holds true for some uniform constantC <˜ 1. Then
(R(¯y)− R(˜y))¯µl+R(˜y)(¯µl−µ˜l)−λ¯l(¯µl−µ˜l) + (˜λl−¯λl)˜µl= ¯gl−g˜l where¯gl,˜gl ∈ X denote the right-hand sides of (3.64e),
¯
gl=hG¯l,ψ¯liX0,Xψ¯l−i−1X G¯l, ˜gl=hG˜l,ψ˜liX0,Xψ˜l−i−1X G˜l implies that
kµ¯l−µ˜lkX ≤ C
1−C˜(kR(¯y)− R(˜y)kLb(X,X)+|¯λl−λ˜l|+k¯gl−g˜lkX); (3.69i) we used here that
˜
µl= (R(˜y)−λ˜l)|−1
{ψ˜l}⊥˜gl
is uniformly bounded as the solution to the compactly perturbed equation (R(¯y)−¯λl)˜µl= ˜gl.
9. The nonlinearity G¯l−G˜l can be estimated by
kG¯l−G˜lkX0 ≤C(k¯y−˜ykH1(Θ,R`)+k¯p−˜pkL2(Θ,R`)
+kψ¯ −ψk˜ V`+k¯u−uk˜ L2(Θ,Rm)+kw¯−wk˜ L2(Θ,Rn)) (3.69j) Combining the estimates (3.69a) to (3.69j), we get
kΦu(˜u,w)˜ −uk¯ U ≤CkBkLb(L2(Θ,Rm),V0)
σu
k˜u−uk¯ U, (3.70) so for sufficiently small control operator norm,kΦu(˜u,w)˜ −uk¯ < holds true.
For the other component ofΦ we have kΦw(˜u,w)˜ −wk¯ W
=kσεw( ¯κw−κ˜w)kW ≤CkI( ¯ψ)¯y−I( ˜ψ)˜ykL2(Θ,Rn)
≤C(kI( ¯ψ)−I( ˜ψ)k
Rn×`k˜ykL2(Θ,R`)+kI( ¯ψ)k
Rn×`k¯y−˜ykL2(Θ,R`))
≤CkIk˜ Lb(V,Rn)(kψ¯ −ψk˜ V`k˜ykL2(Θ,R`)+kψk¯ V`k¯y−yk˜ L2(Θ,R`))
≤CkIk˜ Lb(V,Rn)(kψ¯ −ψk˜ V`+k¯y−yk˜ L2(Θ,R`)) : (3.71) if the constraints operator norm is sufficiently small, kΦw(˜u,w)˜ −wk¯ < holds as well.
Together, we haveΦ(˜u,w)˜ ∈X for all (˜u,w)˜ ∈X, so Φ(X)⊆X. With the same estimates one shows the contraction property
kΦ(ˆu,w)ˆ −Φ(˜u,w)k˜ U×W ≤Ck(ˆˆ u,w)ˆ −(˜u,w)k˜ U×W (3.72) for some constant C <ˆ 1 and all (˜u,w),˜ (ˆu,w)ˆ ∈ X. Now the claim follows with the
Banach fixpoint theorem [7], Thm. 4.3.
The iterative procedure proposed in (3.64) results in the following algorithm:
Algorithm (Optimality System Proper Orthogonal Decomposition)
Require: Control u◦, penaltyw◦, rank `, maximal iterations kmax, tolerance εtol.
1: Initializeu=u◦,w=w◦ &k= 0.
2: repeat
3: Solve the unreduced state equation (3.64a) fory.
4: Determine the solution to the eigenvalue problem (3.64b) ψ,λ.
5: Solve the reduced state equation (3.64c) fory.
6: Solve the reduced adjoint state equation (3.64d) forp.
7: Determine the multipliers of the eigenvalue problem (3.64e)µ.
8: Solve the unreduced adjoint state equation (3.64f) forp.
9: Determine the multipliers of the constraints (3.66a), (3.66b)κu,κw.
10: Update the controlu and the penaltyw (3.65).
11: Determine gradientξ (2.13), perturbation ζ (2.26) and error bound β (2.24).
12: if β < εtol then
13: returnu
14: end if
15: Update k=k+ 1
16: untilk=kmax
Alg. 3.2: An iterative proper orthogonal decomposition strategy which updates the reduced order model in a way which respects optimality of the POD basis elements. Especially, a fixpoint of the updating operator Φ provides a POD basis for which the a-priori error estimates presented in (1.57), (1.60) & (1.64) are applicable. If the corresponding growth constants are estimated in an efficient heuristic way, the POD basis rank does not have to be modified; a proper`can be calculated at the beginning instead.
Remark 3.16. 1. The smallness condition on the source termf and the initial valuey◦
can be skipped by using the homogeneous POD ansatz, see Rem. 3.2.
2. The restriction on the magnitude of B in (3.69i) and (3.70) is not necessary if σu is chosen large enough: In this case, the norm of the operatorR(˜y) is sufficiently small as well since
kuk¯ L2(Θ,Rm) ≤Cσu−1kBkL
b(L2(Θ,Rm),V0)(k¯pkL2(Θ,V)+k¯pkL2(Θ,R`)kψk¯ V`).
3. An alternative approach to bound k¯µl−µ˜lkX is the application of an appropriate perturbation result to
kµ¯l−µ˜lkX ≤ k(R(¯y)−λ¯l)−1|{¯
ψl}⊥¯gl−(R(˜y)−˜λl)−1
|{ψ˜l}⊥g˜lk.
Coupled primal-dual system.
Alg. 3.2 performs much better if the reduced state and adjoint state equations (3.64c)
& (3.64d) are solved as one (nonlinear) system, coupled by the intermediate control-penalty pairu(p) =PU(σu−1BT(ψ)p) andw(y) =PW(I(ψ)y)−I(ψ)ywherePU,PW are the projectors on the admissible domains:
5: Solve the coupled equations (3.64c) & (3.64d) for y(u(p),ψ) and p(y,w(y),ψ).
Alg. 3.2a: Modifications of the OS-POD routine where the linear reduced state and adjoint state equations (5), (6) are replaced by the coupled nonlinear system (3.32).
The contraction property of Φin Thm. 3.15 subsist: (3.69f) becomes k¯y−˜ykH1(Θ,R`)≤C(kBkLbk¯u−˜ukL2(Θ,Rm)+kψ¯ −ψ˜kV`)
≤C(σ−1u kBk2L
bk¯p−pk˜ L2(Θ,R`)+kψ¯ −ψk˜ V`) (3.73a) and (3.69e) changes to
k¯p−pk˜ H1(Θ,R`) ≤C(k¯y−˜ykH1(Θ,R`)+kIkLbk¯w−wk˜ L2(Θ,Rn)+kψ¯ −ψk˜ V`)
≤C((1 +kIk2L
b)k¯y−˜ykH1(Θ,R`)+kψ¯ −ψk˜ V`)
≤C(kBkLb(1 +kIk2L
b)k¯u−˜ukL2(Θ,Rm)+kψ¯ −ψk˜ V`)
≤C(σu−1kBk2L
b(1 +kIk2L
b)k¯p−˜pkL2(Θ,R`)+kψ¯ −ψk˜ V`)
≤Ckψ¯ −ψk˜ V` ≤Ck¯u−uk˜ L2(Θ,Rm). (3.73b) Of course, this coupled solving is much more costly not only because of the nonlinear terms introduced by the variational inequalities. The initial condition of the state equa-tion together with the final condiequa-tion of the adjoint one prevent to solve the coupled system time step by time step, so a time discretization scheme such as the implicit Euler method leads to a system of dimension 2m` for (y,p) instead of m systems of the size 2`.
Remark 3.17. 1. If we skip the control and state constraints and assume thatk·kL2(Θ,V)
and kB? · kU) are equivalent (which, for instance, is the case if U = L2(Θ, V) and B = idU), then we get improved a-priori bounds by testing the state equation in (3.69f) with p = ˜p−p¯ and the adjoint state equation in (3.69e) with y = ˜y−y: In¯ this situation,w¯ = 0 andu¯=σu−1BT( ¯ψ)¯phold and
pTM( ˜ψ) ˙y + pTA( ˜ψ)y = pT(M( ˜ψ)−M( ¯ψ) ˙¯y) + pT(A( ˜ψ)−A( ¯ψ))¯y + pTB( ˜ψ)(˜u−¯u) + pT(B( ˜ψ)−B( ¯ψ))¯u + pT(f( ˜ψ)−f( ¯ψ)),
−yTM( ˜ψ) ˙p + yTM( ˜ψ)p = yT(M( ¯ψ)−M( ˜ψ) ˙¯p) + yT(A( ˜ψ)−A( ¯ψ))¯p
−σQyTM( ˜ψ)y +σQyT(yQ( ˜ψ)−yQ( ¯ψ)).
We assumey◦ = 0, then partial integration yields
0 =−pTB( ˜ψ)B( ˜ψ)Tp−σQyTM( ˜ψ)y−σΩy(T)TM( ˜ψ)y(T) + yT(M( ¯ψ)−M( ˜ψ)) ˙¯p + yT(A( ˜ψ)−A( ¯ψ))¯p
+ pT(M( ¯ψ)−M( ˜ψ)) ˙¯y + pT(A( ¯ψ)−A( ˜ψ)) ˙¯y + pT(B( ¯ψ)−B( ˜ψ))¯u+ pT(f( ¯ψ)−f( ˜ψ)) and therefore
kykL2(Θ,R`)+ky(T)kR`+kpkL2(Θ,R`) ≤Ckψ¯ −ψk˜ V` ≤Ck¯u−uk˜ U.
2. Instead of respecting the primal-dual structure of the problem when testing the state equation with the adjoint state and vice versa, Gong, Hinze & Zhou [55] use the higher regularity which is available in lack of constraints to transform the coupled second-order parabolic system into two independent elliptic equations of second-order four: Assume σΩ = 0, U =L2(Θ×Ω), B = idU and A = −∆. Differentiating the state equation with respect to t and replacing the occuring derivatives of the adjoint state variable by the adjoint state equation yields
−¨y(t) + ∆2y(t) +σu−1(y(t)−yQ(t)) = 0 inΘ×Ω (3.74a) y(t) = 0 & ∆y(t) = 0 onΘ×∂Ω (3.74b)
˙
y(T)−∆y(T) = 0 &y(0)−y◦ = 0 inΩ, (3.74c) analogously, we transform the adjoint state equation to
−¨p(t) + ∆2p(t) +σ−1u p(t) + ˙yQ(t)−∆yQ(t) = 0 inΘ×Ω (3.75a) p(t) = 0 & ∆p(t)−yQ(t) = 0 on Θ×∂Ω (3.75b)
˙
p(0) + ∆p(0) +y◦−yQ(0) = 0 & p(T) = 0 inΩ (3.75c) given that the data of the problem possesses enough regularity. The authors exploit the elliptic structure of (3.74), (3.75) to derive a-priori and a-posteriori error bounds for time discretization schemes and mixed time-space finite element discretizations.
However, their arguments do not apply to constrained optimal control problems. ♦
Projected gradient updates.
Letd= (B?(p+ ˆp)+BT(ψ)p−σuu,−σwε−1w)denote the negative gradient ofJat(u, w).
Kunisch & Volkwein [82] propose not to provide a fixpoint update by postulatingd= 0 in the unconstrained situation or, in our case, projectingd on the admissible domain,
u+=PU(σu−1(B?(p+ ˆp) + BT(ψ)p)) w+=ε−1(PW(I(ψ)y)−I(ψ)y),
but to provide an inexact descent step in direction ofd with a suitable step size σ > 0 which guarantees a sufficient decay of the objective functional:
u+=PW(u+σdu) =PU(u+σ(B?p+ BT(ψ)p)), (3.77a) ω =εw+ I(ψ)y,
ω+=PW(ω+σdω) =Pω(ω+σwε−2(ω−I(ψ)y)),
w+=ε−1(ω+−I(ψ)y). (3.77b)
10: Choose a suitable step lengthσand provide a projected descent step in directiond.
Alg. 3.2b: Modifications of the OS-POD routine where iterative fixpoint updates (10) are replaced by inexact gradient steps (3.77).
In some situations, this variant does not give a significant improvement of accuracy com-pared to the iterative update, but a crucial increase of the calculation time by searching an appropriateσvia backtracking, namely the Armijo step size rule. We therefore do not use this modification in all our numerical tests. Notice that each step size test requires to evaluate the reduced objective functional and hence to solve the reduced state equation.
In Rem. 4.2 we propose a cheaper step length strategy, exploiting the linear-quadratic structure of the optimal control problem.
In [59] Grimm, Gubisch & Volkwein propose to initialize the OS-POD strategy by a few full-order gradient steps – these do not require to solve the coupled system and hence can be provided by solving2m systems of the size n – and to solve (3.32) with a fixed POD basis corresponding to the localized control-penalty pair(u◦, w◦); the descent steps qualify to detect a neighborhood of (¯u,w)¯ which supports good snapshots and suitable initial points for the solving of the nonlinear reduced primal-dual system by the semismooth Newton method, but stagnate within this neighborhood.
Remark 3.18. 1. In [78], Müller and Kunisch use the optimality of the POD basis gained by the OS-POD strategy to derive uniform convergence and polynomial decay rates for the a-priori ROM error of the state variable in terms of`.
2. The formal application of the OS-POD technique does not allow to include snapshots of the adjoint state variable p into the POD basis: In contrast to (3.64d) where a primal-dual sytsem is stated at the beginning and one is free in the selection of the
POD elements, p plays the role of a Lagrange multiplier in the OS-POD context and therefore cannot be part of the eigenvalue constraint (3.64b). In particular, the convergence rates derived in Thm. 3.6 are not applicable. Nevertheless, a multiple snapshot sample including the adjoint state was included in the implementation of Alg. 3.2 in [59] and indeed improved the runtime.
3. The localization of an appropriate initial control-penalty pair (u◦, w◦) can be tricky;
in our numerical tests, we observed divergence if the chosen initialization differed to much from the optimal point(¯u,w).¯
4. Kunisch & Volkwein [82] applied very similar Lagrange techniques to detect optimal snapshot locations for the (semi)discrete POD framework. To reduce the effort of solving the augmented optimality equations, an adaptive approach is proposed where a few optimally chosen time points extend a predefined mesh. In contrast, Hoppe
& Liu [73] conduct the snapshot location by error equilibration, assuring that the discretization errors on the single emerging time intervals essentially coincides; the evaluation of the occuring state errors is provided by hierarchical a-posteriori error
estimates. ♦
Remark 3.19. In the whole chapter we assumed that the control operator B map-ping L2(Θ,Rm) → L2(Θ, V0) has the simplified representation B˜ :Rm → V0 such that (Bu)(t) = ˜B(u(t))holds inV0 for almost allt∈Θ. This assumption allowed to compute the ROM control components(h(Bu)(t), ψliV0,V)l=1,...,`in the form of a productB(ψ)·u(t) where the matrixB(ψ)∈R`×m is given byB(ψ)lk=Bk?(ψl). If we skip this assumption, the evaluation of the control term requires an expensive full-order operation, the dual pairing h·,·iV0,V of (Bu)(t) and ψl, for each single time point during the solving of the reduced state equation. This effect also appears when POD model reduction is applied to nonlinear PDEs where the evaluation of a nonlinear term depending of the reduced state variable, N(y`), requires to expand the reduced coeffient vector y∈ L2(Θ,R`) to y` = P`
l=1ylψl ∈ L2(Θ, V). To avoid calculations on the complexity level of the full-order model, empirical interpolation methods (EIM) have been developed to achieve a reduction basing on a greedy algorithm to the complexity of the POD model. We refer to [15] for an introduction of EIM in the context of the reduced basis method and to [36]
for the application to POD model reduction. ♦
4. Numerical experiments
In this final chapter we illustrate our theoretical findings by several test examples and compare the numerical advantages and disadvantages of the different proposed proce-dures.
4.1. Solution techniques for the optimal control problem
We apply the three solvers for optimization problems which we combined in the OS-POD strategy – the Banach fixpoint iteration, the projected gradient method, and the primal-dual active set strategy – to the full-order optimal control problem, discretized by continuous and piecewise linear finite elements in space and by the first-order implicit Euler scheme in time.