annular and sectorial sparsity in optimal control of ... - TU Chemnitz
A NNULAR AND S ECTORIAL S PARSITY IN O PTIMAL C ONTROL OF E LLIPTIC E QUATIONS Roland Herzog∗
Johannes Obermeier†
Gerd Wachsmuth‡
March 18, 2014
Optimal control problems are considered with linear elliptic equations in polar coordinates. The objective contains L1 -type norms, which promote sparse optimal controls. The particular iterated structure of these norms gives rise to either annular or sectorial sparsity patterns. Optimality conditions and numerical solution approaches are developed.
1 I NTRODUCTION In this paper we consider optimal control problems in which a certain L1 -type norm of the control appears in the objective. Problems of this type are of interest for at least two reasons. Firstly, the L1 norm of the control is often a natural measure of the control cost. Secondly, this term promotes sparsely supported optimal controls, i.e., controls which are zero on substantial parts of its domain of definition. Consequently, control actuators need not be placed everywhere, but only where the control is most effective. Optimal control problems with partial differential equations (PDEs) and sparsity promoting terms were first considered in Stadler [2009], who studied optimality conditions, parameter dependence and a semismooth Newton method in the convex case ∗ Technische
Universität Chemnitz, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), 09107 Chemnitz, Germany, [email protected], http://www.tu-chemnitz.de/herzog † Technische Universität Chemnitz, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), 09107 Chemnitz, Germany ‡ Technische Universität Chemnitz, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), 09107 Chemnitz, Germany, [email protected], http://www.tu-chemnitz.de/mathematik/part_dgl/people/wachsmuth
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
governed by a linear elliptic PDE. A priori and a posteriori error estimates for this case were provided in Wachsmuth and Wachsmuth [2011]. In a sequence of papers Casas et al. [2012b,c], the authors proved second-order necessary and sufficient optimality conditions for the non-convex case governed by a semilinear elliptic equation, and provided a priori finite element error estimates for different choices of the control discretization. More general problem settings involving measure-valued controls were investigated in Clason and Kunisch [2011, 2012]. The subsequent paper Casas et al. [2012a] provides a priori finite element error estimates for the convex case involving a linear elliptic equation. All of the aforementioned papers consider cases in which there is no apparent preference for, nor any control over the shape of the support that the optimal control function may have. This issue was first addressed in Herzog et al. [2012], where striped (directional) sparsity patterns were enforced by way of the iterated norm L1 ( L2 ) in place of the plain L1 norm. An alternative setting which involves the norm L2 (0, T; M(Ω)) was studied in Casas et al. [2013] for problems governed by linear parabolic equations, where M(Ω) is the space of regular Borel measures on Ω. The novelty of the present paper is that we consider a norm of the control in the objective which promotes annular or sectorial sparsity patterns. To achieve this, we use an iterated norm as in Herzog et al. [2012], but switch to a polar coordinate system. To be precise, we consider an optimal control problem with the following objective, 1 2
ZR Z2π 0 0
α |y − yd | dϕ r dr + 2 2
ZR Z2π
2
|u| dϕ r dr + β
0 0
ZR Z2π
|u|2 dϕ
0
1/2
r dr.
(1.1)
0
The parameters α and β are positive constants, and yd represents a desired state. The first and second term in (1.1) are standard. The third term is β times the L1 norm (in radial direction) of the L2 norm (in angular direction) of the control. Since the L1 norm promotes the sparsity of what’s inside, we expect to see optimal controls which are zero on entire annuli centered at the origin. This result will follow from the optimality system proved in Section 2, and we refer to this case as the case with annular sparsity patterns. As a motivation, we mention that annular actuators of piezoelectric type are in use, for instance, in structural health monitoring, to control precision valves, and for the purpose of sound generation or noise cancellation on vibrating structures, see van Niekerk et al. [1995], Coorpender et al. [1999], Raghavan and Cesnik [2004], Yeum et al. [2011], Li et al. [2009] and the references therein. These papers describe the optimization of the actuator geometry as one topic of interest. We emphasize that an approach based on (1.1) would not fix a priori the optimal design to a single annulus but would find the most effective topology simultaneously with the optimal dimensions. We mainly elaborate on the annular case in this paper. The reciprocal situation (the sectorial case), in which the roles of r and ϕ are reversed, is also studied, but presented more briefly. The material is organized as follows. In the rest of this section, we introduce the necessary notation and the precise formulation of the annular problem. Then Section 2
2
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
is devoted to the analysis of the problem, in particular to its optimality conditions. In Section 3 we address the numerical solution of the annular problem by semismooth Newton algorithms, and appropriate finite element discretization. Problems with and without control constraints are treated. Finally, we summarize parallel results and point out the main differences in the proofs for the sectorial case in Section 4. A conclusion and an outlook are given in Section 5. We anticipate that the discussion of the optimal control problems involving (1.1) as their objective requires weighted Lebesgue spaces, and the transformation of the state equation to polar coordinates calls for appropriately transformed Sobolev spaces. This is the main point we focus on in our analytic and algorithmic treatment of the problem. We therefore keep all other ingredients of the problem simple, i.e., we consider only circular domains Ω = { x ∈ R2 : | x | < R}, and we deal with Poisson’s equation only. We comment on some of the numerous interesting extensions of this setting in Section 5.
1.1 N OTATION AND P RELIMINARIES D OMAINS AND T RANSFORMATION Our control problem will be posed on the circle Ω = { x ∈ R2 : | x | < R}, R > 0. We will use a transformation to Cartesian coordinates such that the transformed domain is the rectangle Ω = (0, R) × (0, 2 π ). The transformation from polar coordinates to Cartesian coordinates is given by
( x, y) = P(r, ϕ) = (r cos( ϕ), r sin( ϕ)),
(1.2)
see also Figure 1.1. Note that P maps Ω one-to-one and onto the slit domain Ω \ P
R
2π Ω Ω
0
R
R
P −1
Figure 1.1: The transformation P from polar coordinates to Cartesian coordinates.
[0, R) × {0}. The inverse of P is denoted by P−1 . For later reference, we state the
3
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Jacobians of these mappings as JP (r, ϕ) =
x x 2 + y2 − x2 +y y2
√
cos( ϕ) −r sin( ϕ) , sin( ϕ) r cos( ϕ)
JP−1 ( x, y) =
√
y
x 2 + y2 x x 2 + y2
! .
(1.3)
Note that we have JP−1 = JP−1 ◦ P by the inverse function theorem. We also recall that det JP (r, ϕ) = r and det JP−1 ( x, y) = ( x2 + y2 )−1/2 . M EASURES Due to the transformation P, we will work with two different measures on the rectangle Ω . Both measures will be defined as product measures. We denote by µ and µ ϕ the Lebesgue measures on (0, R) and (0, 2 π ), respectively. By µr we denote the measure on (0, R) defined by Z µr ( A ) =
A
r dµ(r )
(1.4)
for all Lebesgue measurable subsets A ⊂ (0, R). The two product measures we will use on the rectangle Ω are µ × µ ϕ and µr × µ ϕ . Finally, we denote by η the Lebesgue measure on the circle Ω . L EBESGUE S PACES We need several (Bochner-)Lebesgue spaces in order to pose our optimal control problem. First, we define the L2 spaces def
L2 ( Ω ) = L2 ( η )
def
L2r (Ω ) = L2 (µr × µ ϕ ).
and
By using substitution, we now show that these spaces are isometrically isomorphic. Lemma 1.1. The mapping L2 (Ω ) 3 v 7→ v ◦ P ∈ L2r (Ω ) is an isometric isomorphism. Proof. Let v ∈ L2 (Ω ) be given. Since P : Ω : Ω \ [0, R) × {0} is bijective and differentiable, v ◦ P is Lebesgue measurable by [Fremlin, 2003, Thm. 263D(iii)]. The substitution rule shows Z Ω
|v|2 dη =
Z Ω
|v ◦ P|2 |det JP | d(µ × µ ϕ ),
see [Fremlin, 2003, Thm. 263D(v)]. Using det JP = r yields the isometry of the mapping. Using the same arguments for the mapping g 7→ g ◦ P−1 yields the surjectivity.
4
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
We observe that the third term in (1.1) is the norm of the Bochner-Lebesgue space def
1 2 L1,2 r ( Ω ) = L ( µr ; L ( µ ϕ )).
(1.5)
The dual space of this space is L∞ (µr ; L2 (µ ϕ )), see [Diestel and Uhl, 1977, Thm. IV.1.1]. Since µr and µ possess the same null sets, the space L∞ (µr ; L2 (µ ϕ )) coincides with L∞,2 (Ω ) = L∞ (µ; L2 (µ ϕ )) def
with the unweighted Lebesgue measure µ instead of µr . S OBOLEV S PACES We denote by H01 (Ω ) the usual Sobolev space on the circle Ω with homogeneous Dirichlet boundary conditions incorporated. We define b01 (Ω ) def H = {v : Ω → R : v ◦ P−1 ∈ H01 (Ω )}.
(1.6)
b 1 (Ω ) ⊂ L2r (Ω ) = L2 (µr × µ ϕ ). The following lemma shows By Lemma 1.1, we find H 0 b 1 (Ω ) possesses weak derivatives. Moreover, these weak derivatives can be that v ∈ H 0 computed from the weak derivatives of v ◦ P−1 ∈ H01 (Ω ). b 1 (Ω ) be given. Then, v possesses weak derivatives of first order Lemma 1.2. Let v ∈ H 0 and they are ∇v(r, ϕ) = JP> (r, ϕ) ∇(v ◦ P−1 ) ( P(r, ϕ)). Proof. Let ψ ∈ C0∞ (Ω ) be given. We find by substitution Z Ω
v ∇ ψ d( µ × µ ϕ ) =
Z Ω
v ◦ P−1 (∇ψ) ◦ P−1 |det JP−1 | dη.
Using the chain rule (for classical derivatives), we find ∇(ψ ◦ P−1 ) = JP>−1 (∇ψ) ◦ P−1 . This shows
−1 JP−> ) = (∇ψ) ◦ P−1 . −1 ∇( ψ ◦ P
5
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Since ψ ◦ P−1 ∈ C0∞ (Ω ) and v ◦ P−1 ∈ H01 (Ω ), we can use integration by parts. This yields Z Ω
v◦P
−1
(∇ψ) ◦ P
−1
|det JP−1 | dη =
Z
−1 v ◦ P−1 JP−> ) |det JP−1 | dη −1 ∇( ψ ◦ P
Ω
=− =−
Z
h i −1 div v ◦ P−1 JP−> dη −1 |det JP−1 | ψ ◦ P
Ω
Z
h i v ◦ P−1 div JP−> | det J | ψ ◦ P−1 dη − 1 −1 P
Ω
−
Z
−1 JP−> ) |det JP−1 | ψ ◦ P−1 dη, −1 ∇( v ◦ P
Ω
where div denotes the row-wise divergence of a matrix. We have x x 2 + y2 y
JP−> −1 |det JP−1 | =
−√
x 2 + y2
! y x 2 + y2 √ x2 2 , x +y
Hence, we find h i div JP−> | det J | = 0. −1 P −1 This yields Z Ω
v ◦ P−1 (∇ψ) ◦ P−1 |det JP−1 | dη = −
Z Ω
−1 JP−> ) |det JP−1 | ψ ◦ P−1 dη. −1 ∇( v ◦ P
By using substitution, we get Z Ω
v ∇ ψ d( µ × µ ϕ ) = −
Z Ω
JP> ∇(v ◦ P−1 ) ◦ P ψ dη
for all ψ ∈ C0∞ (Ω ), which concludes the proof.
1.2 A NNULAR F ORMULATION OF THE O PTIMAL C ONTROL P ROBLEM S TATE E QUATION For clarity of the presentation, we consider only the case where the state is given as the solution of Poisson’s equation on the circular domain Ω with distributed control. That is, given a control u0 ∈ L2 (Ω ), the state y0 ∈ H01 (Ω ) is the unique solution of Z Ω
∇y0 · ∇v0 dη =
Z Ω
u0 v0 dη
6
for all v0 ∈ H01 (Ω ).
(1.7)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Now, we are going to transform this variational problem into a variational problem b 1 (Ω ) b 1 (Ω ). To this end, we define the bilinear form a on H posed in H 0 0 Z 0 > 1 a(y, v) = ∇y ∇ v d( µr × µ ϕ ). (1.8) 0 r −2 Ω b 1 (Ω ) be given. Then y solves Lemma 1.3. Let u ∈ L2r (Ω ) and y ∈ H 0 a(y, v) =
Z Ω
u v d( µr × µ ϕ )
b01 (Ω ) for all v ∈ H
(1.9)
if and only if y0 := y ◦ P−1 solves (1.7) with right-hand side u0 = u ◦ P−1 . b 1 (Ω ) Proof. The proof is an application of Lemma 1.2, using the definition (1.6) of H 0 and cos( ϕ) sin( ϕ) 1 0 −1 −> −1 , JP (r, ϕ) JP (r, ϕ) = . JP (r, ϕ) = 0 r −2 − sin( ϕ) r −1 cos( ϕ) r −1
b 1 (Ω ) follows from the unique solvability of (1.7) in The unique solvability of (1.9) in H 0 1 H0 (Ω ). Alternatively, we may introduce the norm q kvk Hb 1 (Ω ) = a(v, v) = kv ◦ P−1 k H1 (Ω ) . (1.10) 0
0
b 1 (Ω ) endowed with this norm is a Hilbert space. Now, an It is easy to check that H 0 application of the Lemma of Lax-Milgram yields the unique solvability of (1.9) directly. Corollary 1.4. For every u ∈ L2r (Ω ), the state equation (1.9) possesses a unique sob 1 (Ω ). Moreover, the mapping u 7→ y is continuous from L2r (Ω ) into lution y ∈ H 0 itself. Remark 1.5. The transformation of the state equation (1.7) to polar coordinates (1.9) is not essential for the analysis of the continuous problem in Section 2. However, it facilitates significantly the numerical implementation, see Section 3 and in particular relation (3.6).
S TATEMENT OF THE O PTIMAL C ONTROL P ROBLEM Using the derived state equation (1.9), we may transform the optimal control problem to Ω . With the objective (1.1) and control constraints, the resulting optimal control
7
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
problem becomes Minimize such that and
1 α ky − yd k2L2r (Ω ) + kuk2L2r (Ω 2 2 (y, u) satisfy (1.9)
)
+ β kuk L1,2 r (Ω
)
(P)
u a ≤ u ≤ ub a.e. on Ω .
Here, α > 0 and β ≥ 0 are given constants and the bounds u a , ub ∈ L2r (Ω ) are assumed to satisfy u a < 0 < ub almost everywhere. We define the set of admissible controls as Uad = {u ∈ L2r (Ω ) : u a ≤ u ≤ ub }.
(1.11)
We will also consider the unconstrained case −u a = ub = ∞. In this case, Uad = L2r (Ω ). Note that the assumption u a < 0 < ub is made in order to render u = 0 an admissible control, i.e., in order to allow the desired sparsity. While for this purpose it would be enough to require u a ≤ 0 ≤ ub , some of the results would have to be modified under this relaxed assumption. In particular, the uniqueness of λ¯ of Theorem 2.3 would be lost, and Corollary 2.4 no longer holds. Due to the non-uniqueness, we also expect modifications to be necessary in the numerical solution of the optimality conditions.
2 A NALYSIS OF THE A NNULAR F ORMULATION Using standard arguments, one infers the unique solvability of (P), owing to the uniform convexity of the reduced objective w.r.t. the norm kuk L2r (Ω ) . The optimality con2 ditions will involve the subdifferential of the L1,2 r ( Ω )-norm on Lr ( Ω ), which will be determined using the following lemma.
Lemma 2.1. Let ( X, k·k) be a normed linear space and let |·| : X → [0, ∞) be another norm on X. Then the subdifferential of |·| at x ∈ X is given by ∂|·|( x ) = { x ? ∈ X ? : | x ? |? ≤ 1 and h x, x ? i = | x |},
(2.1)
where X ? is the dual space of X and |·|? : X ? → [0, ∞] is defined by
| x ? |? = sup h x, x ? i.
(2.2)
| x |≤1
Proof. “⊂”: Let x ? ∈ ∂|·|( x ) be given. That is,
| x | + hy − x, x ? i ≤ |y|
8
(2.3)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
holds for all y ∈ X. By taking y = 0 and y = 2 x, we infer
| x | ≤ h x, x ? i and | x | ≥ h x, x ? i. Hence, | x | = h x, x ? i. Using this identity in (2.3), we get hy, x ? i ≤ |y| for all y ∈ X, which shows | x ? |? ≤ 1. “⊃”: Let x ? ∈ X ? be given such that | x ? |? ≤ 1 and h x, x ? i = | x |. This shows hy, x ? i ≤ |y| for all y ∈ X and hence (2.3) is satisfied. Note that the assertion of Lemma 2.1 is standard in case k·k = |·|, see, e.g., [Ioffe and Tichomirov, 1979, p. 20]. We are going to apply Lemma 2.1 with the setting X = L2r (Ω ),
k·k = k·k L2r (Ω ) ,
|·| = k·k L1,2 . r (Ω )
Using the density of L2r (Ω ) in L1,2 r ( Ω ), it follows that the dual norm (2.2) is just the 1,2 norm of the dual of Lr (Ω ), i.e., the L∞,2 (Ω )-norm. Lemma 2.1 now yields Z o n v w d(µr × µ ϕ ) = kvk L1,2 (Ω ) . ∂k·k L1,2 (Ω ) (v) = w ∈ L2r (Ω ) : kwk L∞,2 (Ω ) ≤ 1 and Ω
r
r
(2.4) Note that these conditions imply that equality holds in the chain of inequalities Z R Z 2π 0
0
v w dϕ dµr ≤
Z R 0
kv(r, ·)k L2 (µ ϕ ) kw(r, ·)k L2 (µ ϕ ) dµr ≤ kvk L1,2 kwk L∞,2 (Ω ) . r (Ω )
By this relation, one obtains the following explicit characterization of the subdifferential. Lemma 2.2. Let u ∈ L2r (Ω ) be given. Then, λ ∈ ∂k·k L1,2 (Ω ) (u) holds if and only if r
kλ(r, ·)k L2 (µ ϕ ) ≤ 1 u(r, ·) λ(r, ·) = ku(r, ·)k 2 L (µ ϕ )
where u(r, ·) ≡ 0, elsewhere
(2.5)
for almost all r ∈ (0, R). Proof. The case u ≡ 0 follows directly from (2.4). Let u 6≡ 0 and λ ∈ ∂k·k L1,2 (Ω ) (u) be given. By (2.4) we obtain kλk L∞,2 (Ω ) = 1 by r Hölder’s inequality. The first assertion in (2.5) follows directly from (2.4). By the calculation following (2.4), we infer Z 2π 0
u(r, ϕ) λ(r, ϕ) dϕ = ku(r, ·)k L2 (µ ϕ ) kλ(r, ·)k L2 (µ ϕ )
9
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
for almost all r ∈ (0, R). If u(r, ·) 6≡ 0, this equality implies the existence of c(r ) ≥ 0 with λ(r, ·) = c(r ) u(r, ·) for almost all r ∈ (0, R). By referring to Z R 0
ku(r, ·)k L2 (µ ϕ ) kλ(r, ·)k L2 (µ ϕ ) dµr = kuk L1,2 kλk L∞,2 (Ω r (Ω )
)
we find that u(r, ·) 6≡ 0 implies kλ(r, ·)k L2 (µ ϕ ) = kλk L∞,2 (Ω 1 This shows c(r ) = ku(r, ·)k− L2 ( µ
ϕ)
)
= 1 for almost all r ∈ (0, R).
and, hence, the second relation in (2.5).
Conversely, let λ satisfy (2.5). A direct calculation shows λ ∈ ∂k·k L1,2 (Ω ) (u), see (2.4). r
¯ u¯ ) be the solution of (P). Then, there exists a unique adjoint state Theorem 2.3. Let (y, p¯ and a unique subgradient λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ) such that the system r
a(v, p¯ ) − Z Ω
Z Ω
(yd − y¯ ) v d(µr × µ ϕ ) = 0
(α u¯ − p¯ + β λ¯ ) (u − u¯ ) d(µr × µ ϕ ) ≥ 0 ¯ v) − a(y,
Z Ω
u¯ v d(µr × µ ϕ ) = 0
b01 (Ω ) for all v ∈ H
(2.6a)
for all u ∈ Uad
(2.6b)
b01 (Ω ) for all v ∈ H
(2.6c)
is satisfied. Proof. The existence of p¯ and λ¯ follows from standard arguments and the Theorem of Moreau and Rockafellar, see for instance [Ekeland and Temam, 1999, Proposition I.5.6]. Moreover, the uniqueness of p¯ follows from (2.6a). On the set where u(r, ·) 6≡ 0, Lemma 2.2 shows that λ¯ is unique. On the complement, (2.6b) shows β λ¯ = p¯ since u a < 0 < ub . And hence λ¯ is also unique. Since (P) is convex, the above optimality system is also sufficient. ¯ u¯ ) be the solution of (P). Denote by p¯ the associated adjoint state. Corollary 2.4. Let (y, Then u¯ (r, ·) ≡ 0 ⇔ k p¯ (r, ·)k L2 (µ ϕ ) ≤ β (2.7) holds for almost all r ∈ (0, R).
10
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Proof. Let us denote by λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ) the associated subgradient such that the r optimality system (2.6) is satisfied. Let u¯ (r, ·) ≡ 0 be satisfied for some r ∈ (0, R). By Lemma 2.2 we have kλ¯ (r, ·)k L2 (µ ϕ ) ≤ 1 and by (2.6b) we find α u¯ (r, ·) − p¯ (r, ·) + β λ¯ (r, ·) = 0 since u a < 0 < ub . Putting this together, we get k p¯ (r, ·)k L2 (µ ϕ ) = β kλ¯ (r, ·)k L2 (µ ϕ ) ≤ β. It remains to prove the converse. Let us define N = {r ∈ (0, R) : k p¯ (r, ·)k L2 (µ ϕ ) ≤ β}. Using the test function (recall u a < 0 < ub ) ( 0 r ∈ N, u(r, ϕ) = u¯ (r, ϕ) else in the variational inequality (2.6b) we obtain Z N ×(0,2 π )
(α u¯ − p¯ + β λ¯ ) u¯ d(µr × µ ϕ ) ≤ 0
Lemma 2.2 implies Z N ×(0,2 π )
λ¯ u¯ d(µr × µ ϕ ) ≥
Z N
ku¯ (r, ·)k L2 (µ ϕ ) dµr .
Hence, we have 0≥
Z
≥α
N ×(0,2 π )
(α u¯ − p¯ + β λ¯ ) u¯ d(µr × µ ϕ )
Z N ×(0,2 π )
|u¯ |2 d(µr × µ ϕ ) − β
Z N
ku¯ (r, ·)k L2 (µ ϕ ) dµr + β
Z N
ku¯ (r, ·)k L2 (µ ϕ ) dµr .
This shows u¯ (r, ·) ≡ 0 on N. As expected, (2.7) implies that the optimal control u¯ is sparse. Moreover, we infer the annular sparsity structure, since (2.7) implies that we have u¯ (r, ϕ) ≡ 0 for all ϕ ∈ (0, 2 π ) whenever k p¯ (r, ·)k L2 (µ ϕ ) ≤ β holds. In the unconstrained case, i.e., when Uad = L2r (Ω ) holds (or formally, −u a = ub = ∞), it is straightforward to show that (2.6b) implies α u¯ − p¯ + β λ¯ = 0. In this case, we will exploit the following reformulation of this optimality system numerically, see Section 3.3. ¯ u¯ ) be the solution of (P) in the case Uad = L2r (Ω ). Then, there Lemma 2.5. Let (y, exists a unique adjoint state p¯ and a unique subgradient λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ) such that r (2.6a), (2.6c) and β α u¯ (r, ϕ) = max 0, 1 − p¯ (r, ϕ) a.e. in Ω (2.8) k p¯ (r, ·)k L2 (µ ϕ )
11
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
are satisfied. Proof. In view of Theorem 2.3, we only need to show that (2.8) is equivalent to (2.6b). This is done by a distinction of cases. First, we consider the set N = {r ∈ (0, R) : u¯ (r, ·) ≡ 0}. As already seen in the proof of Corollary 2.4, (2.6b) is equivalent to k p¯ (r, ·)k L2 (µ ϕ ) ≤ β on this set. It is easy to see that also (2.8) is equivalent to this condition. On the complement set (0, R) \ N, (2.6b) is equivalent to α u¯ + β
u¯ = p¯ ¯ ku(r, ·)k L2 (µ ϕ )
(∗)
by Lemma 2.2. Similarly, (2.8) is equivalent to α u¯ = p¯ − β
p¯ . k p¯ (r, ·)k L2 (µ ϕ )
(∗∗)
Both (∗) and (∗∗) imply u¯ p¯ = . ku¯ (r, ·)k L2 (µ ϕ ) k p¯ (r, ·)k L2 (µ ϕ ) Now, the equivalence of (2.6b) and (2.8) is easy to see via the equivalence of (∗) and (∗∗). We remark that (2.8) implies the continuity of u¯ ◦ P−1 in Ω if p¯ ◦ P−1 is continuous. We mention that the Newton differentiability in function space of the system (2.6a), (2.8) and (2.6c) can be shown analogously as in [Herzog et al., 2012, Section 3], which justifies the use of a semi-smooth Newton method in the case without control constraints. We briefly sketch the ideas for this case. To this end, let us introduce the (affine) adjointstate map b01 (Ω ), P : L2r (Ω ) 3 u 7→ p = P(u) ∈ p ∈ H which assigns to a given control the unique solution of the adjoint equation through (2.6c) and (2.6a). Then, in view of (2.8), the entire optimality system (2.6) can be formulated equivalently as the single condition F (u¯ ) = 0, where F : L2r (Ω ) 3 u 7→ F (u) := u − α−1 G ( P(u)) ∈ L2r (Ω ), β ) 3 p 7 → G ( p ) : = max 0, 1 − p ∈ L2r (Ω ). G : L6,2 ( Ω r k p(r, ·)k L2 (µ ϕ ) Here L6,2 r ( Ω ) denotes the Bochner-Lebesgue space analogous to (1.5). By the standard embedding H01 (Ω ) ,→ L6 (Ω ) and a simple adaptation of Lemma 1.1 we get b 1 (Ω ) ,→ L6r (Ω ), and clearly also H b 1 (Ω ) ,→ L6,2 H r ( Ω ) holds. This shows that F 0 0
12
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
has the mapping properties as claimed above. The Newton differentiability of G is shown in [Herzog et al., 2012, Lemma 3.2]. The Newton differentiability of F and the bounded invertibility of its generalized derivative with k F 0 (u)−1 kL( L2r (Ω )) ≤ 1 follows as in [Herzog et al., 2012, Lemma 3.6].
3 N UMERICAL R EALIZATION AND R ESULTS We address in this section the discretization of the annular optimal control problem and its solution by a semi-smooth Newton method.
3.1 F INITE E LEMENT D ISCRETIZATION OF THE F ORWARD P ROBLEM We start with a triangular grid of the rectangular domain Ω = (0, R) × (0, 2π ). To solve the variational formulation of the forward problem (1.9) in a conforming way, we b 1 (Ω ). To this end, our starting point is a need to construct a discrete subspace of H 0 space of piecewise linear functions vh on Ω with the standard nodal basis. In order to ensure the conformity, i.e., vh ◦ P−1 ∈ H01 (Ω ), we need to enforce the continuity of the piecewise smooth function vh ◦ P−1 , as well as the essential boundary conditions. Both can be achieved by appropriate conditions on the function vh , viz. 1. The Dirichlet conditions on vh ◦ P−1 are satisfied if and only if vh ( R, ·) = 0 holds. 2. Since the points (r, 0) and (r, 2π ) are mapped by P to the same point in Ω , vh (r, 0) = vh (r, 2π ) must necessarily hold for all r ∈ [0, R]. 3. Similarly, all points (0, ϕ) are mapped by P to the origin in Ω , vh (0, ϕ) must be independent of ϕ ∈ [0, 2π ]. Conditions 1. and 3. can be easily formulated by appropriate restrictions on the degrees of freedom describing the function vh . To ensure condition 2. we impose the following condition on the mesh on Ω . We require that the vertices on the lower boundary (r, 0) agree in their r-coordinate with the vertices (r, 2π ) on the upper boundary.
(M)
Then condition 2. is simply realized by a number of equality constraints for the degrees of freedom located in the relative interior of these boundaries. All these constraints are sketched in Figure 3.1. The resulting finite element space is termed Vh . Note that the requirement that the vertices on the lower and upper boundaries be aligned is important for the approximation properties of the ensuing subspaces. When the vertices are not aligned, the continuity condition 2. would imply that r 7→ vh (r, 0) is globally linear. We replace the variational forward problem (1.9) by its Galerkin approximation, i.e., find yh ∈ Vh such that Z Z 0 > r ∇yh ∇vh dr dϕ = u vh r dr dϕ for all vh ∈ Vh . (3.1) 0 r −1 Ω Ω
13
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
2π
=
const
=
=
0
=
0
R
Figure 3.1: Conditions on the discrete functions on Ω . Let us briefly describe how the solution of (3.1) can be achieved in the M ATLAB PDE toolbox. As usual, we denote by [p,e,t] the mesh data. We denote by K the stiffness matrix associated with the left hand side of (3.1), and with all degrees of freedom present and unconstrained, i.e., K = assema (p ,t , char ( ’x ’ , ’ 1./ x ’) ,0 ,0) ; To incorporate the constraints of type 1.–3. one could set up a matrix N where each row represents one of the constraints. The solution of the state equation in Vh is would then be achieved by solving the augmented system b K N> y = , z 0 N 0 where b is the load vector generated by the right hand side in (3.1). Equivalently, with Z being a matrix whose columns span the nullspace of N, we could solve Z> K Z y e = Z> b and then expand y := Z y e. A suitable nullspace basis can be constructed easily. Recall that each row in Z belongs to one of the degrees of freedom (nodes) in the mesh. There are three kinds of columns in Z: 1. each inner node generates a column with exactly one entry ’1’, at the row corresponding to its index, 2. each pair of nodes in the interior of the upper and lower boundaries of Ω generates a column with exactly two entries ’1’, at the rows corresponding to their indices, 3. the set of all nodes on the left boundary r = 0 generates one more column, with entries ’1’ in all rows pertaining to participating nodes. Each coefficient vector in the range space of Z represents a piecewise linear function vh on Ω w.r.t. the standard nodal basis, and with the property that vh ◦ P−1 ∈ H01 (Ω ) holds.
14
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
3.2 C ONVERGENCE OF THE D ISCRETIZATION This subsection is new. In order to show the convergence of the above scheme for the forward problem, we can b 1 (Ω ). To this end, we have to show apply Céa’s Lemma in H 0 inf kv − vh k Hb 1 (Ω
vh ∈Vh
0
)
→ 0 as h → 0
(3.2)
b 1 (Ω ). Let us introduce for all v ∈ H 0 b0∞ (Ω ) = {v : Ω → R such that v ◦ P−1 ∈ C0∞ (Ω )}. C b 1 (Ω ) and the density of C ∞ (Ω ) in H 1 (Ω ), By the definition (1.10) of the norm in H 0 0 0 b∞ (Ω ). Due to the smoothness of P, it is easy it is sufficient to verify (3.2) for all v ∈ C 0 b∞ (Ω ) implies to show that v ∈ C 0 • v ∈ C ∞ ( Ω ),
• v is constant along the line r = 0, • the values of v and its derivatives at (r, 0) coincide with the values at (r, 2 π ), for all r ∈ [0, R], • v is zero in a neighborhood of the line r = R, see also Figure 3.1. Due to these conditions, the Lagrange interpolation of functions in b∞ (Ω ) is well defined for meshes satisfying condition (M). C 0
Lemma 3.1. Let {Th } be a quasi-uniform family of geometrically conforming triangular meshes on Ω which satisfy condition (M). Let us denote by vh the Lagrange inter∞ polant of v ∈ Cc 0 ( Ω ). Then we have
kv − vh k2Hb 1 (Ω 0
)
≤ C B2 h2 (1 + |ln h|),
where C depends only on c1 , c2 and R, and B is the ∞-norm of the second derivatives of v. Proof. First, we consider an arbitrary triangle with vertices (ri , ϕi ), i = 1, 2, 3, such that, w.l.o.g., r1 ≤ r2 ≤ r3 . We set r2 − r1 r3 − r1 J = det ϕ2 − ϕ1 ϕ3 − ϕ1 and denote by 1 ( ϕ3 − ϕ1 ) (r − r1 ) + (r1 − r3 ) ( ϕ − ϕ1 ) J 1 λ3 (r, ϕ) = ( ϕ1 − ϕ2 ) (r − r1 ) + (r2 − r1 ) ( ϕ − ϕ1 ) J λ1 (r, ϕ) = 1 − λ2 (r, ϕ) − λ3 (r, ϕ)
λ2 (r, ϕ) =
15
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
the barycentric coordinates. Then, vh is given by v h = λ1 v (r1 , ϕ1 ) + λ2 v (r2 , ϕ2 ) + λ3 v (r3 , ϕ3 ) and it is easy to verify that r3 − r2 r1 − r3 ∂vh r2 − r1 = v (r1 , ϕ1 ) + v (r2 , ϕ2 ) + v (r3 , ϕ3 ) ∂ϕ J J J ∂v = (r1 , ϕ1 ) + C B h. ∂ϕ
(3.3a) (3.3b)
We remark that the quasi-uniformity of the mesh family implies that there exist c1 , c2 > 0 such that the edge lengths of all triangles in Th belong to [c1 h, c2 h]. In order to estimate the interpolation error, we distinguish three different cases. We denote by Th,i the set of all triangles of a given mesh which satisfy the conditions of case #i. In what follows, C denotes a generic constant which depends only on the quantities c1 , c2 , R and which may change from line to line. Case 1: 0 = r1 = r2 Referring to (3.3a) and using v(0, ϕ) = 0, we obtain ∂vh /∂ϕ = 0. Moreover, we get 2 ∂v ∂v ∂ v (r, ϕ) ≤ (r1 , ϕ) + ∂ϕ ∂ϕ ∂r ∂ϕ ∞ |r − r1 | ≤ B r by a Taylor estimate. This yields for any triangle 4 ∈ Th,1 Z 4
Z 2 1 ∂v ∂vh 2 1 ∂v (r, ϕ) − dr dϕ = (r, ϕ) dr dϕ r ∂ϕ ∂ϕ 4 r ∂ϕ
≤B
2
Z 4
r dr dϕ ≤ C B2 h3 .
Since | ϕ1 − ϕ2 | is bounded from below by c1 h, there are at most 2 π/(c1 h) = C/h triangles in Th,1 . This shows
∑
Z
4∈Th,1 4
∂vh 2 1 ∂v (r, ϕ) − dr dϕ ≤ C B2 h2 . r ∂ϕ ∂ϕ
Case 2: 0 = r1 < r2 By (3.3b) we obtain ∂vh ∂v ≤ (r1 , ϕ1 ) + C B h ≤ C B h ∂ϕ ∂ϕ and as in Case 1
∂v (r, ϕ) ≤ B r ≤ C B h. ∂ϕ
16
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Let us denote by w(s) the length of the intersection of the triangle with some line r = s. We have Z Z r3 Z ∂vh 2 1 1 1 ∂v 2 2 2 2 dr dϕ = C B h w(r ) dr ≤ C B2 h3 . (r, ϕ) − dr dϕ ≤ C B h r ∂ϕ ∂ϕ r r 4 0 4 Here we used w(r ) ≤ C r and r3 ≤ c2 h, which follow from the assumptions on the mesh family. Since all triangles in Th,2 lie in the strip [0, c2 h] × [0, 2 π ], there are at most (c2 h 2 π )/( β h2 ) = C/h triangles in this case, where β > 0 is chosen such that the area of any triangle is bounded from below by β h2 . This yields Z 1 ∂v ∂vh 2 ( r, ϕ ) − dr dϕ ≤ C B2 h2 . ∑ 4 r ∂ϕ ∂ϕ 4∈T h,2
Case 3: 0 < r1 Using (3.3b), we get Z Z 1 ∂v ∂vh 2 C ∂vh 2 C B2 h4 ∂v . (r, ϕ) − dr dϕ ≤ (r, ϕ) − dr dϕ ≤ ∂ϕ r1 4 ∂ϕ ∂ϕ r1 4 r ∂ϕ For each triangle in Th,3 , there must be another triangle in the mesh which has (r1 , ϕ1 ) as a vertex and lies to the left of (r1 , ϕ1 ). More precisely, this triangle contains a nontrivial portion of the line segment connecting (r1 , ϕ1 ) with the boundary point (0, ϕ1 ). Due to the assumptions on the mesh regularity, this implies r1 ≥ α h for some α > 0 which depends only on c1 , c2 . Now, the set Th,3 is further divided into subsets, based on the value of r1 . In particular, we set n Th,3 = 4 ∈ Th,3 : r1 ∈ [α n h, α (n + 1) h) Sb
R
c
n αh for n = 1, . . . , b αRh c. Since b αRh c > αRh − 1, we get Th,3 = n= 1 T h,3 . Moreover, since n all triangles in Th,3 lie in the strip [α n h, α (n + 1) h + c2 h] × [0, 2 π ], there are at most n . We obtain (α + c2 ) h 2 π/( β h2 ) = C/h triangles in Th,3 Z ∂vh 2 C C B2 h4 C B2 h2 1 ∂v ( r, ϕ ) − dr dϕ ≤ · ≤ . ∑ r ∂ϕ ∂ϕ h αnh n 4∈T n 4 h,3
By using an upper bound for the harmonic series b αRh c
∑
n =1
R 1 ≤ 1 + lnb c ≤ C + |ln h| n αh
we obtain
∑
Z
4∈Th,3 4
b αRh c Z ∂vh 2 ∂vh 2 1 ∂v 1 ∂v (r, ϕ) − dr dϕ ≤ ∑ ∑ (r, ϕ) − dr dϕ r ∂ϕ ∂ϕ r ∂ϕ ∂ϕ n=1 4∈T n 4 h,3
≤ C B2 h2
b αRh c
∑
n =1
17
1 ≤ C B2 h2 (1 + |ln h|). n
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
In all cases, a Taylor estimate yields Z
∂v ∂v 2 r (r, ϕ) − h dr dϕ ≤ C B2 h4 ∂r ∂r 4
Summing over all triangles yields the claim. As already discussed above, this result implies the convergence of the discretization scheme of the forward problem as h → 0. Based on this, the convergence of an optimal solution in a discretize–then optimize setting can be deduced as well using standard arguments. In the sequel, we will however follow an optimize–then discretize approach for which the discussion of convergence is more involved and beyond the scope of this paper.
3.3 D ISCRETIZATION AND S OLUTION OF THE U NCONSTRAINED O PTIMALITY S YSTEM To complete the discretization of (P), we follow an optimize–then discretize approach. In this section, we consider the case without control constraints. We employ piecewise linear controls on the same triangular grid which we use for the state. Therefore, the right hand side in the discrete state equation (3.1) can be realized by the term M u, where M is mass matrix [ ∼ ,M ] = assema (p ,t ,0 , char ( ’x ’) ,0) ; We do not impose additional conditions on the discrete control (as we did for the state) because its continuity will follow automatically from the discrete optimality system given below. To preserve the iterated structure of the term k p¯ (r, ·)k L2 (µ ϕ ) in (2.8), we choose a particular quadrature formula and impose further structural conditions on the mesh. From now on, the mesh vertices are supposed to form a rectangular lattice. This supersedes the conditions set forth in Section 3.1. For simplicity, we elaborate on the case of constant mesh widths hr and h ϕ . It will be convenient to address the components of the coefficient vector u by a double index. The value of the control at (r, ϕ) = (i hr , j h ϕ ) is denoted by uij with 1 ≤ i ≤ nr and 1 ≤ j ≤ n ϕ . e , which are exThe discrete state and adjoint state are represented by vectors y e and p panded to the full nodal basis by multiplication with Z. Our discrete optimality system now consists of the discrete state equation, Z> K Z y e = ZM u,
y = Zy e,
(3.4)
the discrete adjoint equation, e = −ZM (y − yd ), Z> K Z p
18
e, p = Zp
(3.5)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
and the following discretization of (2.8), β
!
α uij = max 0, 1 − n ϕ 1/2 pij =: G(p) ij pij = G(p) p ij . ∑ ωk p2ik
(3.6)
k =1
Here the weights ωk are equal to h ϕ /2 for k ∈ {1, n ϕ } and ωk = h ϕ otherwise. The vector yd represents the coefficient of a nodal interpolation of the desired state yd . In the last term in (3.6), G(p) denotes a vector and is the pointwise product between vectors of the same size. Note that (3.6) implies that the discrete control inherits the continuity of the discrete adjoint state, which in turn follows from the fact that its coefficient vector p is in the null space of the constraint matrix N, and that the factor in parentheses depends only on the index i in radial direction. We solve (3.4)–(3.6) by a semismooth Newton method. We mention that the matrix in the resulting Newton system for the update step > δe y 0 Z> K Z Z MZ 0 e ) (3.7) αI diag(p) G0 (p) Z + diag(G(p)) Z δu = −R(y e, u, p δe p Z> K Z −Z> M 0 turns out to be non-symmetric, and not obviously symmetrizable. Here G0 (p) denotes a Newton derivative of G, and the right hand side is the negative residual of (3.5), (3.6) and (3.4) in this order. Numerical experiments revealed that globalization efforts were not required for convergence.
3.4 C ASE OF P OINTWISE C ONTROL C ONSTRAINTS The numerical treatment of (P) is more involved in the presence of pointwise inequality constraints for the control. The reason is that these constraints are not compatible with the sparsity structure induced by the iterated norm in the objective. This was observed already in Herzog et al. [2012]. Following Lemma 4.1 in that paper, we reformulate the optimality system (2.6) as a non-smooth equation. One can show that the following two statements are equivalent: 1. λ¯ ∈ ∂k·k 1,2 (u¯ ) and (2.6) holds. Lr ( Ω )
2. (2.6a) and (2.6c) and
− p¯ + α u¯ + β λ¯ + µ¯ = 0 max 1, kλ¯ (r, ·) + c1 u¯ (r, ·)k L2 (µ ϕ ) λ¯ − (λ¯ + c1 u¯ ) = 0 µ¯ − max 0, µ¯ + c2 (u¯ − ub ) − min 0, µ¯ + c2 (u¯ − u a ) = 0 hold for any choice of positive constants c1 , c2 .
19
(3.8a) (3.8b) (3.8c)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
We mention that (3.8b) is equivalent to λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ), and (3.8c) is equivalent to µ¯ r ¯ This can be shown by a straightforward being an element of the normal cone of Uad at u. distinction of cases. Unfortunately, both versions of the optimality system have their drawbacks. The first version is not directly amenable to numerical treatment due to the variational inequality in (2.6b). For the second version, the viability of the semismooth Newton method in function space is unknown. Nevertheless, a semismooth Newton type iteration can be successfully applied in the discrete setting. To this end, we discretize the optimality system in a similar way as before in Section 3.3. The additional variables λ¯ and µ¯ are of the same dimension as the control u¯ in the discrete setting. The term involving the k·k L2 (µ ϕ ) in (3.8b) is treated in the same way as in (2.8). The derivation of the discrete semismooth Newton system proceeds in the same way as in the unconstrained case, and we refer to [Herzog et al., 2012, Section 4] for details. As was observed already there, two modifications should be applied to the Newton system in order to ensure its well-posedness in each step as well as the global convergence of the method in practice. Firstly, note that (3.8b) implies u¯ = 0 on the subset of [0, R] where the ’max’ attains the value one. On the other hand, (3.8c) implies u¯ ∈ {u a , ub } on some subsets of Ω , which are termed active sets. This may lead to contradictory conditions at intermediate iterations, and a singular Newton matrix ensues. Therefore, we give preference to (3.8b) and modify the determination of the active sets in (3.8c) so that they can only be subsets of {(r, ϕ) : kλ¯ (r, ·) + c1 u¯ (r, ·)k L2 (µ ϕ ) < 1}. The second modification concerns the linearization of (3.8b), to which a damping is applied. We refer to [Herzog et al., 2012, eq. (4.10)] for details. We point out that both modifications vanish in the limit and they do not impair fast local convergence. In Figure 3.2 we show the optimal control for the sectorial problem (with and without control constraints) obtained for the following setting, α = 0.01, u a = −1,
β = 0.15, yd ( x, y) = e2 x sin(y π ), ub = 1.
(3.9)
In the complementarity formulation (3.8), we used the parameters α c1 = and c2 = 100. β
4 S ECTORIAL F ORMULATION In this section we point out which changes to (P) are necessary to obtain optimal controls with sectorial sparsity patterns. As can be expected, we have to modify the sparsity promoting term kuk L1,2 (Ω ) in the objective. To this end, we define the space r
def 1 2 L˜ 1,2 r ( Ω ) = L ( µ ϕ ; L ( µr )).
20
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Figure 3.2: Optimal control of the annular sparsity problem (P) without (left plot) and with control constraints (right plot). The parameters are given in (3.9). 1 2 Compared with L1,2 r ( Ω ) = L ( µr ; L ( µ ϕ )) we have changed the order of integration. The norm in the space L˜ 1,2 r ( Ω ) is given by
kuk L˜ 1,2 r (Ω
)
=
Z 2π Z R 0
0
2
|u(r, ϕ)| dµr
1/2
dµ ϕ =
Z 2π Z R 0
0
|u(r, ϕ)|2 r dµ
1/2
dµ ϕ .
The optimal control problem now reads Minimize such that and
1 α ky − yd k2L2r (Ω ) + kuk2L2r (Ω 2 2 (y, u) satisfy (1.9)
)
+ β kuk L˜ 1,2 r (Ω
)
(P)
u a ≤ u ≤ ub a.e. on Ω .
Following the analysis in Section 2, we have to identify the dual space of L˜ 1,2 r ( Ω ) and 2 the subdifferential of its norm in Lr (Ω ). Standard arguments yield that the dual space def is L˜ r∞,2 (Ω ) = L∞ (µ ϕ ; L2 (µr )), and using Lemma 2.1, we have n
∂k·k L˜ 1,2 (Ω ) (v) = w ∈ r
L2r (Ω
) : kwk L˜ r∞,2 (Ω
)
≤ 1 and
Z Ω
v w d(µr × µ ϕ ) = kvk L˜ 1,2 (Ω r
o )
(4.1) Analogously to the computations in the proof of Lemma 2.2, we find that, for any given u ∈ L2r (Ω ), the function λ ∈ L2r (Ω ) belongs to the subdifferential ∂k·k L˜ 1,2 (Ω ) (u) if r and only if where u(·, ϕ) ≡ 0, kλ(·, ϕ)k L2 (µr ) ≤ 1 u(·, ϕ) elsewhere λ(·, ϕ) = ku(·, ϕ)k 2 L ( µr ) holds for almost all ϕ ∈ (0, 2π ). The same arguments as used in the proof of Theorem 2.3 show that the optimality system (2.6) is the same as in the annular formulation, except that now λ¯ ∈ ∂k·k L˜ 1,2 (Ω ) (u) r
21
.
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
holds. Following the calculations leading to (2.7), we also obtain the equivalence u¯ (·, ϕ) ≡ 0
⇔
k p¯ (·, ϕ)k L2 (µr ) ≤ β
(4.2)
for almost all ϕ ∈ (0, 2π ), which confirms the occurrence of sectorial sparsity patterns. The analog of formula (2.8) in the absence of control constraints reads α u¯ (r, ϕ) = max 0, 1 −
β p¯ (r, ϕ) k p¯ (·, ϕ)k L2 (µr )
a.e. in Ω .
(4.3)
However, (4.3) does not imply the continuity of the (transformed) optimal control u¯ ◦ P−1 at the origin, even if p¯ ◦ P−1 is continuous. To see this, note that the continuity of p¯ ◦ P−1 implies that p¯ (0, ϕ) is constant, but (4.3) does not imply that u¯ (0, ϕ) is constant w.r.t. ϕ since the max(. . .) term depends on ϕ. The numerical treatment of the sectorial problem requires only a few changes compared to the annular variant. The restrictions we impose on the mesh and the discrete state and adjoint states, as well as the formulation of the discrete state equation are the same as in Section 3.1. Concerning the numerical algorithm, we simply need to replace (2.8) by (4.3) in the unconstrained case and proceed as before. In the constrained case, we merely need to replace (3.8b) by max 1, kλ¯ (·, ϕ) + c1 u¯ (·, ϕ)k L2 (µr ) λ¯ − (λ¯ + c1 u¯ ) = 0 and continue as in Section 3.4. In Figure 4.1 we show the optimal solution for the following data, α = 0.05, u a = −1,
β = 0.02, yd ( x, y) = e2 x sin(y π ), ub = 1.
(4.4)
In the complementarity formulation, we used the parameters c1 =
α β
and
c2 = 10.
5 C ONCLUSIONS AND O UTLOOK We introduced and analyzed convex optimal control problems whose optimal controls feature particular spatial sparsity patterns. Optimality systems were derived, which led to semismooth Newton type methods for the numerical solution. We elaborated on the cases of annular and sectorial sparsity patterns on circular domains by way of the polar coordinate transform and an iterated norm in the objective. The extension to domains which are described by rectangles ( R1 , R2 ) × ( ϕ1 , ϕ2 ) in polar coordinates (annular sectors) is straightforward. The conditions set forth in Section 3.1
22
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Figure 4.1: Optimal control of the sectorial sparsity problem (P) without (left plot) and with control constraints (right plot). The parameters are given in (4.4). to achieve conformity need to be adjusted. In fact, there is no need to use only rectangles as Ω . And finally, the polar coordinate transformation could be replaced by a general diffeomorphism to accommodate domains and sparsity patterns of rather arbitrary shape. In three dimensions, the variety of choices by which the three coordinates can be grouped into outer (sparse) and inner coordinates increases from two to six. For example, when we use spherical coordinates and place the radial coordinate outside and the two angular coordinates inside, the support of an optimal control will consist of spherical shells.
R EFERENCES E. Casas, C. Clason, and K. Kunisch. Approximation of elliptic control problems in measure spaces with sparse solutions. SIAM Journal on Control and Optimization, 50 (4):1735–1752, 2012a. ISSN 0363-0129. doi: 10.1137/110843216. E. Casas, R. Herzog, and G. Wachsmuth. Optimality conditions and error analysis of semilinear elliptic control problems with L1 cost functional. SIAM Journal on Optimization, 22(3):795–820, 2012b. doi: 10.1137/110834366. E. Casas, R. Herzog, and G. Wachsmuth. Approximation of sparse controls in semilinear equations by piecewise linear functions. Numerische Mathematik, 122(4):645–669, 2012c. doi: 10.1007/s00211-012-0475-7. E. Casas, C. Clason, and K. Kunisch. Parabolic control problems in measure spaces with sparse solutions. SIAM Journal on Control and Optimization, 51(1):28–63, 2013. ISSN 0363-0129. doi: 10.1137/120872395. C. Clason and K. Kunisch. A duality-based approach to elliptic control problems in non-reflexive Banach spaces. ESAIM: Control, Optimisation, and Calculus of Variations, 17(1):243–266, 2011. doi: 10.1051/cocv/2010003.
23
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
C. Clason and K. Kunisch. A measure space approach to optimal source placement. Computational Optimization and Applications. An International Journal, 53(1):155–171, 2012. ISSN 0926-6003. doi: 10.1007/s10589-011-9444-9. S.J. Coorpender, D. Finkel, J. Kyzar, R. Sims, A.B. Smirnova, M. Tawhid, C.E. Bouton, and R.C. Smith. Modeling and optimization issues concerning a circular piezoelectric actuator design. Technical Report CRSC–TR99–22, North Carolina State University, 1999. URL http://www.ncsu.edu/crsc/reports/reports99.html. J. Diestel and J.J. Uhl. Vector Measures. Mathematical Surveys and Monographs. American Mathematical Society, Providence, 1977. I. Ekeland and R. Temam. Convex Analysis and Variational Problems, volume 28 of Classics in Applied Mathematics. SIAM, Philadelphia, 1999. D. H. Fremlin. Measure theory. Vol. 2. Torres Fremlin, Colchester, 2003. ISBN 0-95381292-8. Broad foundations, Corrected second printing of the 2001 original. R. Herzog, G. Stadler, and G. Wachsmuth. Directional sparsity in optimal control of partial differential equations. SIAM Journal on Control and Optimization, 50(2):943– 963, 2012. doi: 10.1137/100815037. A. D. Ioffe and V. M. Tichomirov. Theorie der Extremalaufgaben. VEB Deutscher Verlag der Wissenschaften, Berlin, 1979. X. Li, H. Du, L. Xu, Y. Hu, and L. Xu. Optimization of a circular thin-film piezoelectric actuator lying on a clamped multilayered elastic plate. IEEE Transactions on Ultrasonics, 56(7):1469–1475, 2009. A. Raghavan and C. E. S. Cesnik. Modeling of piezoelectric-based lamb wave generation and sensing for structural health monitoring. In Smart Structures and Materials 2004: Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems, volume 5391, pages 419–430, 2004. doi: 10.1117/12.540269. G. Stadler. Elliptic optimal control problems with L1 -control cost and applications for the placement of control devices. Computational Optimization and Applications, 44(2): 159–181, 2009. doi: 10.1007/s10589-007-9150-9. J.L. van Niekerk, B.H. Tongue, and A.K. Packard. Active control of a circular plate to reduce transient noise transmission. Journal of Sound and Vibration, 183(4):643 – 662, 1995. doi: 10.1006/jsvi.1995.0277. G. Wachsmuth and D. Wachsmuth. Convergence and regularization results for optimal control problems with sparsity functional. ESAIM: Control, Optimisation and Calculus of Variations, 17(3):858–886, 2011. doi: 10.1051/cocv/2010027. C. M. Yeum, H. Sohn, and J. B. Ihn. Lamb wave mode decomposition using concentric ring and circular piezoelectric transducers. Wave Motion. An International Journal
24
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Reporting Research on Wave Phenomena, 48(4):358–370, 2011. ISSN 0165-2125. doi: 10.1016/j.wavemoti.2011.01.001.
25
Johannes Obermeier†
Gerd Wachsmuth‡
March 18, 2014
Optimal control problems are considered with linear elliptic equations in polar coordinates. The objective contains L1 -type norms, which promote sparse optimal controls. The particular iterated structure of these norms gives rise to either annular or sectorial sparsity patterns. Optimality conditions and numerical solution approaches are developed.
1 I NTRODUCTION In this paper we consider optimal control problems in which a certain L1 -type norm of the control appears in the objective. Problems of this type are of interest for at least two reasons. Firstly, the L1 norm of the control is often a natural measure of the control cost. Secondly, this term promotes sparsely supported optimal controls, i.e., controls which are zero on substantial parts of its domain of definition. Consequently, control actuators need not be placed everywhere, but only where the control is most effective. Optimal control problems with partial differential equations (PDEs) and sparsity promoting terms were first considered in Stadler [2009], who studied optimality conditions, parameter dependence and a semismooth Newton method in the convex case ∗ Technische
Universität Chemnitz, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), 09107 Chemnitz, Germany, [email protected], http://www.tu-chemnitz.de/herzog † Technische Universität Chemnitz, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), 09107 Chemnitz, Germany ‡ Technische Universität Chemnitz, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), 09107 Chemnitz, Germany, [email protected], http://www.tu-chemnitz.de/mathematik/part_dgl/people/wachsmuth
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
governed by a linear elliptic PDE. A priori and a posteriori error estimates for this case were provided in Wachsmuth and Wachsmuth [2011]. In a sequence of papers Casas et al. [2012b,c], the authors proved second-order necessary and sufficient optimality conditions for the non-convex case governed by a semilinear elliptic equation, and provided a priori finite element error estimates for different choices of the control discretization. More general problem settings involving measure-valued controls were investigated in Clason and Kunisch [2011, 2012]. The subsequent paper Casas et al. [2012a] provides a priori finite element error estimates for the convex case involving a linear elliptic equation. All of the aforementioned papers consider cases in which there is no apparent preference for, nor any control over the shape of the support that the optimal control function may have. This issue was first addressed in Herzog et al. [2012], where striped (directional) sparsity patterns were enforced by way of the iterated norm L1 ( L2 ) in place of the plain L1 norm. An alternative setting which involves the norm L2 (0, T; M(Ω)) was studied in Casas et al. [2013] for problems governed by linear parabolic equations, where M(Ω) is the space of regular Borel measures on Ω. The novelty of the present paper is that we consider a norm of the control in the objective which promotes annular or sectorial sparsity patterns. To achieve this, we use an iterated norm as in Herzog et al. [2012], but switch to a polar coordinate system. To be precise, we consider an optimal control problem with the following objective, 1 2
ZR Z2π 0 0
α |y − yd | dϕ r dr + 2 2
ZR Z2π
2
|u| dϕ r dr + β
0 0
ZR Z2π
|u|2 dϕ
0
1/2
r dr.
(1.1)
0
The parameters α and β are positive constants, and yd represents a desired state. The first and second term in (1.1) are standard. The third term is β times the L1 norm (in radial direction) of the L2 norm (in angular direction) of the control. Since the L1 norm promotes the sparsity of what’s inside, we expect to see optimal controls which are zero on entire annuli centered at the origin. This result will follow from the optimality system proved in Section 2, and we refer to this case as the case with annular sparsity patterns. As a motivation, we mention that annular actuators of piezoelectric type are in use, for instance, in structural health monitoring, to control precision valves, and for the purpose of sound generation or noise cancellation on vibrating structures, see van Niekerk et al. [1995], Coorpender et al. [1999], Raghavan and Cesnik [2004], Yeum et al. [2011], Li et al. [2009] and the references therein. These papers describe the optimization of the actuator geometry as one topic of interest. We emphasize that an approach based on (1.1) would not fix a priori the optimal design to a single annulus but would find the most effective topology simultaneously with the optimal dimensions. We mainly elaborate on the annular case in this paper. The reciprocal situation (the sectorial case), in which the roles of r and ϕ are reversed, is also studied, but presented more briefly. The material is organized as follows. In the rest of this section, we introduce the necessary notation and the precise formulation of the annular problem. Then Section 2
2
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
is devoted to the analysis of the problem, in particular to its optimality conditions. In Section 3 we address the numerical solution of the annular problem by semismooth Newton algorithms, and appropriate finite element discretization. Problems with and without control constraints are treated. Finally, we summarize parallel results and point out the main differences in the proofs for the sectorial case in Section 4. A conclusion and an outlook are given in Section 5. We anticipate that the discussion of the optimal control problems involving (1.1) as their objective requires weighted Lebesgue spaces, and the transformation of the state equation to polar coordinates calls for appropriately transformed Sobolev spaces. This is the main point we focus on in our analytic and algorithmic treatment of the problem. We therefore keep all other ingredients of the problem simple, i.e., we consider only circular domains Ω = { x ∈ R2 : | x | < R}, and we deal with Poisson’s equation only. We comment on some of the numerous interesting extensions of this setting in Section 5.
1.1 N OTATION AND P RELIMINARIES D OMAINS AND T RANSFORMATION Our control problem will be posed on the circle Ω = { x ∈ R2 : | x | < R}, R > 0. We will use a transformation to Cartesian coordinates such that the transformed domain is the rectangle Ω = (0, R) × (0, 2 π ). The transformation from polar coordinates to Cartesian coordinates is given by
( x, y) = P(r, ϕ) = (r cos( ϕ), r sin( ϕ)),
(1.2)
see also Figure 1.1. Note that P maps Ω one-to-one and onto the slit domain Ω \ P
R
2π Ω Ω
0
R
R
P −1
Figure 1.1: The transformation P from polar coordinates to Cartesian coordinates.
[0, R) × {0}. The inverse of P is denoted by P−1 . For later reference, we state the
3
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Jacobians of these mappings as JP (r, ϕ) =
x x 2 + y2 − x2 +y y2
√
cos( ϕ) −r sin( ϕ) , sin( ϕ) r cos( ϕ)
JP−1 ( x, y) =
√
y
x 2 + y2 x x 2 + y2
! .
(1.3)
Note that we have JP−1 = JP−1 ◦ P by the inverse function theorem. We also recall that det JP (r, ϕ) = r and det JP−1 ( x, y) = ( x2 + y2 )−1/2 . M EASURES Due to the transformation P, we will work with two different measures on the rectangle Ω . Both measures will be defined as product measures. We denote by µ and µ ϕ the Lebesgue measures on (0, R) and (0, 2 π ), respectively. By µr we denote the measure on (0, R) defined by Z µr ( A ) =
A
r dµ(r )
(1.4)
for all Lebesgue measurable subsets A ⊂ (0, R). The two product measures we will use on the rectangle Ω are µ × µ ϕ and µr × µ ϕ . Finally, we denote by η the Lebesgue measure on the circle Ω . L EBESGUE S PACES We need several (Bochner-)Lebesgue spaces in order to pose our optimal control problem. First, we define the L2 spaces def
L2 ( Ω ) = L2 ( η )
def
L2r (Ω ) = L2 (µr × µ ϕ ).
and
By using substitution, we now show that these spaces are isometrically isomorphic. Lemma 1.1. The mapping L2 (Ω ) 3 v 7→ v ◦ P ∈ L2r (Ω ) is an isometric isomorphism. Proof. Let v ∈ L2 (Ω ) be given. Since P : Ω : Ω \ [0, R) × {0} is bijective and differentiable, v ◦ P is Lebesgue measurable by [Fremlin, 2003, Thm. 263D(iii)]. The substitution rule shows Z Ω
|v|2 dη =
Z Ω
|v ◦ P|2 |det JP | d(µ × µ ϕ ),
see [Fremlin, 2003, Thm. 263D(v)]. Using det JP = r yields the isometry of the mapping. Using the same arguments for the mapping g 7→ g ◦ P−1 yields the surjectivity.
4
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
We observe that the third term in (1.1) is the norm of the Bochner-Lebesgue space def
1 2 L1,2 r ( Ω ) = L ( µr ; L ( µ ϕ )).
(1.5)
The dual space of this space is L∞ (µr ; L2 (µ ϕ )), see [Diestel and Uhl, 1977, Thm. IV.1.1]. Since µr and µ possess the same null sets, the space L∞ (µr ; L2 (µ ϕ )) coincides with L∞,2 (Ω ) = L∞ (µ; L2 (µ ϕ )) def
with the unweighted Lebesgue measure µ instead of µr . S OBOLEV S PACES We denote by H01 (Ω ) the usual Sobolev space on the circle Ω with homogeneous Dirichlet boundary conditions incorporated. We define b01 (Ω ) def H = {v : Ω → R : v ◦ P−1 ∈ H01 (Ω )}.
(1.6)
b 1 (Ω ) ⊂ L2r (Ω ) = L2 (µr × µ ϕ ). The following lemma shows By Lemma 1.1, we find H 0 b 1 (Ω ) possesses weak derivatives. Moreover, these weak derivatives can be that v ∈ H 0 computed from the weak derivatives of v ◦ P−1 ∈ H01 (Ω ). b 1 (Ω ) be given. Then, v possesses weak derivatives of first order Lemma 1.2. Let v ∈ H 0 and they are ∇v(r, ϕ) = JP> (r, ϕ) ∇(v ◦ P−1 ) ( P(r, ϕ)). Proof. Let ψ ∈ C0∞ (Ω ) be given. We find by substitution Z Ω
v ∇ ψ d( µ × µ ϕ ) =
Z Ω
v ◦ P−1 (∇ψ) ◦ P−1 |det JP−1 | dη.
Using the chain rule (for classical derivatives), we find ∇(ψ ◦ P−1 ) = JP>−1 (∇ψ) ◦ P−1 . This shows
−1 JP−> ) = (∇ψ) ◦ P−1 . −1 ∇( ψ ◦ P
5
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Since ψ ◦ P−1 ∈ C0∞ (Ω ) and v ◦ P−1 ∈ H01 (Ω ), we can use integration by parts. This yields Z Ω
v◦P
−1
(∇ψ) ◦ P
−1
|det JP−1 | dη =
Z
−1 v ◦ P−1 JP−> ) |det JP−1 | dη −1 ∇( ψ ◦ P
Ω
=− =−
Z
h i −1 div v ◦ P−1 JP−> dη −1 |det JP−1 | ψ ◦ P
Ω
Z
h i v ◦ P−1 div JP−> | det J | ψ ◦ P−1 dη − 1 −1 P
Ω
−
Z
−1 JP−> ) |det JP−1 | ψ ◦ P−1 dη, −1 ∇( v ◦ P
Ω
where div denotes the row-wise divergence of a matrix. We have x x 2 + y2 y
JP−> −1 |det JP−1 | =
−√
x 2 + y2
! y x 2 + y2 √ x2 2 , x +y
Hence, we find h i div JP−> | det J | = 0. −1 P −1 This yields Z Ω
v ◦ P−1 (∇ψ) ◦ P−1 |det JP−1 | dη = −
Z Ω
−1 JP−> ) |det JP−1 | ψ ◦ P−1 dη. −1 ∇( v ◦ P
By using substitution, we get Z Ω
v ∇ ψ d( µ × µ ϕ ) = −
Z Ω
JP> ∇(v ◦ P−1 ) ◦ P ψ dη
for all ψ ∈ C0∞ (Ω ), which concludes the proof.
1.2 A NNULAR F ORMULATION OF THE O PTIMAL C ONTROL P ROBLEM S TATE E QUATION For clarity of the presentation, we consider only the case where the state is given as the solution of Poisson’s equation on the circular domain Ω with distributed control. That is, given a control u0 ∈ L2 (Ω ), the state y0 ∈ H01 (Ω ) is the unique solution of Z Ω
∇y0 · ∇v0 dη =
Z Ω
u0 v0 dη
6
for all v0 ∈ H01 (Ω ).
(1.7)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Now, we are going to transform this variational problem into a variational problem b 1 (Ω ) b 1 (Ω ). To this end, we define the bilinear form a on H posed in H 0 0 Z 0 > 1 a(y, v) = ∇y ∇ v d( µr × µ ϕ ). (1.8) 0 r −2 Ω b 1 (Ω ) be given. Then y solves Lemma 1.3. Let u ∈ L2r (Ω ) and y ∈ H 0 a(y, v) =
Z Ω
u v d( µr × µ ϕ )
b01 (Ω ) for all v ∈ H
(1.9)
if and only if y0 := y ◦ P−1 solves (1.7) with right-hand side u0 = u ◦ P−1 . b 1 (Ω ) Proof. The proof is an application of Lemma 1.2, using the definition (1.6) of H 0 and cos( ϕ) sin( ϕ) 1 0 −1 −> −1 , JP (r, ϕ) JP (r, ϕ) = . JP (r, ϕ) = 0 r −2 − sin( ϕ) r −1 cos( ϕ) r −1
b 1 (Ω ) follows from the unique solvability of (1.7) in The unique solvability of (1.9) in H 0 1 H0 (Ω ). Alternatively, we may introduce the norm q kvk Hb 1 (Ω ) = a(v, v) = kv ◦ P−1 k H1 (Ω ) . (1.10) 0
0
b 1 (Ω ) endowed with this norm is a Hilbert space. Now, an It is easy to check that H 0 application of the Lemma of Lax-Milgram yields the unique solvability of (1.9) directly. Corollary 1.4. For every u ∈ L2r (Ω ), the state equation (1.9) possesses a unique sob 1 (Ω ). Moreover, the mapping u 7→ y is continuous from L2r (Ω ) into lution y ∈ H 0 itself. Remark 1.5. The transformation of the state equation (1.7) to polar coordinates (1.9) is not essential for the analysis of the continuous problem in Section 2. However, it facilitates significantly the numerical implementation, see Section 3 and in particular relation (3.6).
S TATEMENT OF THE O PTIMAL C ONTROL P ROBLEM Using the derived state equation (1.9), we may transform the optimal control problem to Ω . With the objective (1.1) and control constraints, the resulting optimal control
7
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
problem becomes Minimize such that and
1 α ky − yd k2L2r (Ω ) + kuk2L2r (Ω 2 2 (y, u) satisfy (1.9)
)
+ β kuk L1,2 r (Ω
)
(P)
u a ≤ u ≤ ub a.e. on Ω .
Here, α > 0 and β ≥ 0 are given constants and the bounds u a , ub ∈ L2r (Ω ) are assumed to satisfy u a < 0 < ub almost everywhere. We define the set of admissible controls as Uad = {u ∈ L2r (Ω ) : u a ≤ u ≤ ub }.
(1.11)
We will also consider the unconstrained case −u a = ub = ∞. In this case, Uad = L2r (Ω ). Note that the assumption u a < 0 < ub is made in order to render u = 0 an admissible control, i.e., in order to allow the desired sparsity. While for this purpose it would be enough to require u a ≤ 0 ≤ ub , some of the results would have to be modified under this relaxed assumption. In particular, the uniqueness of λ¯ of Theorem 2.3 would be lost, and Corollary 2.4 no longer holds. Due to the non-uniqueness, we also expect modifications to be necessary in the numerical solution of the optimality conditions.
2 A NALYSIS OF THE A NNULAR F ORMULATION Using standard arguments, one infers the unique solvability of (P), owing to the uniform convexity of the reduced objective w.r.t. the norm kuk L2r (Ω ) . The optimality con2 ditions will involve the subdifferential of the L1,2 r ( Ω )-norm on Lr ( Ω ), which will be determined using the following lemma.
Lemma 2.1. Let ( X, k·k) be a normed linear space and let |·| : X → [0, ∞) be another norm on X. Then the subdifferential of |·| at x ∈ X is given by ∂|·|( x ) = { x ? ∈ X ? : | x ? |? ≤ 1 and h x, x ? i = | x |},
(2.1)
where X ? is the dual space of X and |·|? : X ? → [0, ∞] is defined by
| x ? |? = sup h x, x ? i.
(2.2)
| x |≤1
Proof. “⊂”: Let x ? ∈ ∂|·|( x ) be given. That is,
| x | + hy − x, x ? i ≤ |y|
8
(2.3)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
holds for all y ∈ X. By taking y = 0 and y = 2 x, we infer
| x | ≤ h x, x ? i and | x | ≥ h x, x ? i. Hence, | x | = h x, x ? i. Using this identity in (2.3), we get hy, x ? i ≤ |y| for all y ∈ X, which shows | x ? |? ≤ 1. “⊃”: Let x ? ∈ X ? be given such that | x ? |? ≤ 1 and h x, x ? i = | x |. This shows hy, x ? i ≤ |y| for all y ∈ X and hence (2.3) is satisfied. Note that the assertion of Lemma 2.1 is standard in case k·k = |·|, see, e.g., [Ioffe and Tichomirov, 1979, p. 20]. We are going to apply Lemma 2.1 with the setting X = L2r (Ω ),
k·k = k·k L2r (Ω ) ,
|·| = k·k L1,2 . r (Ω )
Using the density of L2r (Ω ) in L1,2 r ( Ω ), it follows that the dual norm (2.2) is just the 1,2 norm of the dual of Lr (Ω ), i.e., the L∞,2 (Ω )-norm. Lemma 2.1 now yields Z o n v w d(µr × µ ϕ ) = kvk L1,2 (Ω ) . ∂k·k L1,2 (Ω ) (v) = w ∈ L2r (Ω ) : kwk L∞,2 (Ω ) ≤ 1 and Ω
r
r
(2.4) Note that these conditions imply that equality holds in the chain of inequalities Z R Z 2π 0
0
v w dϕ dµr ≤
Z R 0
kv(r, ·)k L2 (µ ϕ ) kw(r, ·)k L2 (µ ϕ ) dµr ≤ kvk L1,2 kwk L∞,2 (Ω ) . r (Ω )
By this relation, one obtains the following explicit characterization of the subdifferential. Lemma 2.2. Let u ∈ L2r (Ω ) be given. Then, λ ∈ ∂k·k L1,2 (Ω ) (u) holds if and only if r
kλ(r, ·)k L2 (µ ϕ ) ≤ 1 u(r, ·) λ(r, ·) = ku(r, ·)k 2 L (µ ϕ )
where u(r, ·) ≡ 0, elsewhere
(2.5)
for almost all r ∈ (0, R). Proof. The case u ≡ 0 follows directly from (2.4). Let u 6≡ 0 and λ ∈ ∂k·k L1,2 (Ω ) (u) be given. By (2.4) we obtain kλk L∞,2 (Ω ) = 1 by r Hölder’s inequality. The first assertion in (2.5) follows directly from (2.4). By the calculation following (2.4), we infer Z 2π 0
u(r, ϕ) λ(r, ϕ) dϕ = ku(r, ·)k L2 (µ ϕ ) kλ(r, ·)k L2 (µ ϕ )
9
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
for almost all r ∈ (0, R). If u(r, ·) 6≡ 0, this equality implies the existence of c(r ) ≥ 0 with λ(r, ·) = c(r ) u(r, ·) for almost all r ∈ (0, R). By referring to Z R 0
ku(r, ·)k L2 (µ ϕ ) kλ(r, ·)k L2 (µ ϕ ) dµr = kuk L1,2 kλk L∞,2 (Ω r (Ω )
)
we find that u(r, ·) 6≡ 0 implies kλ(r, ·)k L2 (µ ϕ ) = kλk L∞,2 (Ω 1 This shows c(r ) = ku(r, ·)k− L2 ( µ
ϕ)
)
= 1 for almost all r ∈ (0, R).
and, hence, the second relation in (2.5).
Conversely, let λ satisfy (2.5). A direct calculation shows λ ∈ ∂k·k L1,2 (Ω ) (u), see (2.4). r
¯ u¯ ) be the solution of (P). Then, there exists a unique adjoint state Theorem 2.3. Let (y, p¯ and a unique subgradient λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ) such that the system r
a(v, p¯ ) − Z Ω
Z Ω
(yd − y¯ ) v d(µr × µ ϕ ) = 0
(α u¯ − p¯ + β λ¯ ) (u − u¯ ) d(µr × µ ϕ ) ≥ 0 ¯ v) − a(y,
Z Ω
u¯ v d(µr × µ ϕ ) = 0
b01 (Ω ) for all v ∈ H
(2.6a)
for all u ∈ Uad
(2.6b)
b01 (Ω ) for all v ∈ H
(2.6c)
is satisfied. Proof. The existence of p¯ and λ¯ follows from standard arguments and the Theorem of Moreau and Rockafellar, see for instance [Ekeland and Temam, 1999, Proposition I.5.6]. Moreover, the uniqueness of p¯ follows from (2.6a). On the set where u(r, ·) 6≡ 0, Lemma 2.2 shows that λ¯ is unique. On the complement, (2.6b) shows β λ¯ = p¯ since u a < 0 < ub . And hence λ¯ is also unique. Since (P) is convex, the above optimality system is also sufficient. ¯ u¯ ) be the solution of (P). Denote by p¯ the associated adjoint state. Corollary 2.4. Let (y, Then u¯ (r, ·) ≡ 0 ⇔ k p¯ (r, ·)k L2 (µ ϕ ) ≤ β (2.7) holds for almost all r ∈ (0, R).
10
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Proof. Let us denote by λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ) the associated subgradient such that the r optimality system (2.6) is satisfied. Let u¯ (r, ·) ≡ 0 be satisfied for some r ∈ (0, R). By Lemma 2.2 we have kλ¯ (r, ·)k L2 (µ ϕ ) ≤ 1 and by (2.6b) we find α u¯ (r, ·) − p¯ (r, ·) + β λ¯ (r, ·) = 0 since u a < 0 < ub . Putting this together, we get k p¯ (r, ·)k L2 (µ ϕ ) = β kλ¯ (r, ·)k L2 (µ ϕ ) ≤ β. It remains to prove the converse. Let us define N = {r ∈ (0, R) : k p¯ (r, ·)k L2 (µ ϕ ) ≤ β}. Using the test function (recall u a < 0 < ub ) ( 0 r ∈ N, u(r, ϕ) = u¯ (r, ϕ) else in the variational inequality (2.6b) we obtain Z N ×(0,2 π )
(α u¯ − p¯ + β λ¯ ) u¯ d(µr × µ ϕ ) ≤ 0
Lemma 2.2 implies Z N ×(0,2 π )
λ¯ u¯ d(µr × µ ϕ ) ≥
Z N
ku¯ (r, ·)k L2 (µ ϕ ) dµr .
Hence, we have 0≥
Z
≥α
N ×(0,2 π )
(α u¯ − p¯ + β λ¯ ) u¯ d(µr × µ ϕ )
Z N ×(0,2 π )
|u¯ |2 d(µr × µ ϕ ) − β
Z N
ku¯ (r, ·)k L2 (µ ϕ ) dµr + β
Z N
ku¯ (r, ·)k L2 (µ ϕ ) dµr .
This shows u¯ (r, ·) ≡ 0 on N. As expected, (2.7) implies that the optimal control u¯ is sparse. Moreover, we infer the annular sparsity structure, since (2.7) implies that we have u¯ (r, ϕ) ≡ 0 for all ϕ ∈ (0, 2 π ) whenever k p¯ (r, ·)k L2 (µ ϕ ) ≤ β holds. In the unconstrained case, i.e., when Uad = L2r (Ω ) holds (or formally, −u a = ub = ∞), it is straightforward to show that (2.6b) implies α u¯ − p¯ + β λ¯ = 0. In this case, we will exploit the following reformulation of this optimality system numerically, see Section 3.3. ¯ u¯ ) be the solution of (P) in the case Uad = L2r (Ω ). Then, there Lemma 2.5. Let (y, exists a unique adjoint state p¯ and a unique subgradient λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ) such that r (2.6a), (2.6c) and β α u¯ (r, ϕ) = max 0, 1 − p¯ (r, ϕ) a.e. in Ω (2.8) k p¯ (r, ·)k L2 (µ ϕ )
11
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
are satisfied. Proof. In view of Theorem 2.3, we only need to show that (2.8) is equivalent to (2.6b). This is done by a distinction of cases. First, we consider the set N = {r ∈ (0, R) : u¯ (r, ·) ≡ 0}. As already seen in the proof of Corollary 2.4, (2.6b) is equivalent to k p¯ (r, ·)k L2 (µ ϕ ) ≤ β on this set. It is easy to see that also (2.8) is equivalent to this condition. On the complement set (0, R) \ N, (2.6b) is equivalent to α u¯ + β
u¯ = p¯ ¯ ku(r, ·)k L2 (µ ϕ )
(∗)
by Lemma 2.2. Similarly, (2.8) is equivalent to α u¯ = p¯ − β
p¯ . k p¯ (r, ·)k L2 (µ ϕ )
(∗∗)
Both (∗) and (∗∗) imply u¯ p¯ = . ku¯ (r, ·)k L2 (µ ϕ ) k p¯ (r, ·)k L2 (µ ϕ ) Now, the equivalence of (2.6b) and (2.8) is easy to see via the equivalence of (∗) and (∗∗). We remark that (2.8) implies the continuity of u¯ ◦ P−1 in Ω if p¯ ◦ P−1 is continuous. We mention that the Newton differentiability in function space of the system (2.6a), (2.8) and (2.6c) can be shown analogously as in [Herzog et al., 2012, Section 3], which justifies the use of a semi-smooth Newton method in the case without control constraints. We briefly sketch the ideas for this case. To this end, let us introduce the (affine) adjointstate map b01 (Ω ), P : L2r (Ω ) 3 u 7→ p = P(u) ∈ p ∈ H which assigns to a given control the unique solution of the adjoint equation through (2.6c) and (2.6a). Then, in view of (2.8), the entire optimality system (2.6) can be formulated equivalently as the single condition F (u¯ ) = 0, where F : L2r (Ω ) 3 u 7→ F (u) := u − α−1 G ( P(u)) ∈ L2r (Ω ), β ) 3 p 7 → G ( p ) : = max 0, 1 − p ∈ L2r (Ω ). G : L6,2 ( Ω r k p(r, ·)k L2 (µ ϕ ) Here L6,2 r ( Ω ) denotes the Bochner-Lebesgue space analogous to (1.5). By the standard embedding H01 (Ω ) ,→ L6 (Ω ) and a simple adaptation of Lemma 1.1 we get b 1 (Ω ) ,→ L6r (Ω ), and clearly also H b 1 (Ω ) ,→ L6,2 H r ( Ω ) holds. This shows that F 0 0
12
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
has the mapping properties as claimed above. The Newton differentiability of G is shown in [Herzog et al., 2012, Lemma 3.2]. The Newton differentiability of F and the bounded invertibility of its generalized derivative with k F 0 (u)−1 kL( L2r (Ω )) ≤ 1 follows as in [Herzog et al., 2012, Lemma 3.6].
3 N UMERICAL R EALIZATION AND R ESULTS We address in this section the discretization of the annular optimal control problem and its solution by a semi-smooth Newton method.
3.1 F INITE E LEMENT D ISCRETIZATION OF THE F ORWARD P ROBLEM We start with a triangular grid of the rectangular domain Ω = (0, R) × (0, 2π ). To solve the variational formulation of the forward problem (1.9) in a conforming way, we b 1 (Ω ). To this end, our starting point is a need to construct a discrete subspace of H 0 space of piecewise linear functions vh on Ω with the standard nodal basis. In order to ensure the conformity, i.e., vh ◦ P−1 ∈ H01 (Ω ), we need to enforce the continuity of the piecewise smooth function vh ◦ P−1 , as well as the essential boundary conditions. Both can be achieved by appropriate conditions on the function vh , viz. 1. The Dirichlet conditions on vh ◦ P−1 are satisfied if and only if vh ( R, ·) = 0 holds. 2. Since the points (r, 0) and (r, 2π ) are mapped by P to the same point in Ω , vh (r, 0) = vh (r, 2π ) must necessarily hold for all r ∈ [0, R]. 3. Similarly, all points (0, ϕ) are mapped by P to the origin in Ω , vh (0, ϕ) must be independent of ϕ ∈ [0, 2π ]. Conditions 1. and 3. can be easily formulated by appropriate restrictions on the degrees of freedom describing the function vh . To ensure condition 2. we impose the following condition on the mesh on Ω . We require that the vertices on the lower boundary (r, 0) agree in their r-coordinate with the vertices (r, 2π ) on the upper boundary.
(M)
Then condition 2. is simply realized by a number of equality constraints for the degrees of freedom located in the relative interior of these boundaries. All these constraints are sketched in Figure 3.1. The resulting finite element space is termed Vh . Note that the requirement that the vertices on the lower and upper boundaries be aligned is important for the approximation properties of the ensuing subspaces. When the vertices are not aligned, the continuity condition 2. would imply that r 7→ vh (r, 0) is globally linear. We replace the variational forward problem (1.9) by its Galerkin approximation, i.e., find yh ∈ Vh such that Z Z 0 > r ∇yh ∇vh dr dϕ = u vh r dr dϕ for all vh ∈ Vh . (3.1) 0 r −1 Ω Ω
13
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
2π
=
const
=
=
0
=
0
R
Figure 3.1: Conditions on the discrete functions on Ω . Let us briefly describe how the solution of (3.1) can be achieved in the M ATLAB PDE toolbox. As usual, we denote by [p,e,t] the mesh data. We denote by K the stiffness matrix associated with the left hand side of (3.1), and with all degrees of freedom present and unconstrained, i.e., K = assema (p ,t , char ( ’x ’ , ’ 1./ x ’) ,0 ,0) ; To incorporate the constraints of type 1.–3. one could set up a matrix N where each row represents one of the constraints. The solution of the state equation in Vh is would then be achieved by solving the augmented system b K N> y = , z 0 N 0 where b is the load vector generated by the right hand side in (3.1). Equivalently, with Z being a matrix whose columns span the nullspace of N, we could solve Z> K Z y e = Z> b and then expand y := Z y e. A suitable nullspace basis can be constructed easily. Recall that each row in Z belongs to one of the degrees of freedom (nodes) in the mesh. There are three kinds of columns in Z: 1. each inner node generates a column with exactly one entry ’1’, at the row corresponding to its index, 2. each pair of nodes in the interior of the upper and lower boundaries of Ω generates a column with exactly two entries ’1’, at the rows corresponding to their indices, 3. the set of all nodes on the left boundary r = 0 generates one more column, with entries ’1’ in all rows pertaining to participating nodes. Each coefficient vector in the range space of Z represents a piecewise linear function vh on Ω w.r.t. the standard nodal basis, and with the property that vh ◦ P−1 ∈ H01 (Ω ) holds.
14
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
3.2 C ONVERGENCE OF THE D ISCRETIZATION This subsection is new. In order to show the convergence of the above scheme for the forward problem, we can b 1 (Ω ). To this end, we have to show apply Céa’s Lemma in H 0 inf kv − vh k Hb 1 (Ω
vh ∈Vh
0
)
→ 0 as h → 0
(3.2)
b 1 (Ω ). Let us introduce for all v ∈ H 0 b0∞ (Ω ) = {v : Ω → R such that v ◦ P−1 ∈ C0∞ (Ω )}. C b 1 (Ω ) and the density of C ∞ (Ω ) in H 1 (Ω ), By the definition (1.10) of the norm in H 0 0 0 b∞ (Ω ). Due to the smoothness of P, it is easy it is sufficient to verify (3.2) for all v ∈ C 0 b∞ (Ω ) implies to show that v ∈ C 0 • v ∈ C ∞ ( Ω ),
• v is constant along the line r = 0, • the values of v and its derivatives at (r, 0) coincide with the values at (r, 2 π ), for all r ∈ [0, R], • v is zero in a neighborhood of the line r = R, see also Figure 3.1. Due to these conditions, the Lagrange interpolation of functions in b∞ (Ω ) is well defined for meshes satisfying condition (M). C 0
Lemma 3.1. Let {Th } be a quasi-uniform family of geometrically conforming triangular meshes on Ω which satisfy condition (M). Let us denote by vh the Lagrange inter∞ polant of v ∈ Cc 0 ( Ω ). Then we have
kv − vh k2Hb 1 (Ω 0
)
≤ C B2 h2 (1 + |ln h|),
where C depends only on c1 , c2 and R, and B is the ∞-norm of the second derivatives of v. Proof. First, we consider an arbitrary triangle with vertices (ri , ϕi ), i = 1, 2, 3, such that, w.l.o.g., r1 ≤ r2 ≤ r3 . We set r2 − r1 r3 − r1 J = det ϕ2 − ϕ1 ϕ3 − ϕ1 and denote by 1 ( ϕ3 − ϕ1 ) (r − r1 ) + (r1 − r3 ) ( ϕ − ϕ1 ) J 1 λ3 (r, ϕ) = ( ϕ1 − ϕ2 ) (r − r1 ) + (r2 − r1 ) ( ϕ − ϕ1 ) J λ1 (r, ϕ) = 1 − λ2 (r, ϕ) − λ3 (r, ϕ)
λ2 (r, ϕ) =
15
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
the barycentric coordinates. Then, vh is given by v h = λ1 v (r1 , ϕ1 ) + λ2 v (r2 , ϕ2 ) + λ3 v (r3 , ϕ3 ) and it is easy to verify that r3 − r2 r1 − r3 ∂vh r2 − r1 = v (r1 , ϕ1 ) + v (r2 , ϕ2 ) + v (r3 , ϕ3 ) ∂ϕ J J J ∂v = (r1 , ϕ1 ) + C B h. ∂ϕ
(3.3a) (3.3b)
We remark that the quasi-uniformity of the mesh family implies that there exist c1 , c2 > 0 such that the edge lengths of all triangles in Th belong to [c1 h, c2 h]. In order to estimate the interpolation error, we distinguish three different cases. We denote by Th,i the set of all triangles of a given mesh which satisfy the conditions of case #i. In what follows, C denotes a generic constant which depends only on the quantities c1 , c2 , R and which may change from line to line. Case 1: 0 = r1 = r2 Referring to (3.3a) and using v(0, ϕ) = 0, we obtain ∂vh /∂ϕ = 0. Moreover, we get 2 ∂v ∂v ∂ v (r, ϕ) ≤ (r1 , ϕ) + ∂ϕ ∂ϕ ∂r ∂ϕ ∞ |r − r1 | ≤ B r by a Taylor estimate. This yields for any triangle 4 ∈ Th,1 Z 4
Z 2 1 ∂v ∂vh 2 1 ∂v (r, ϕ) − dr dϕ = (r, ϕ) dr dϕ r ∂ϕ ∂ϕ 4 r ∂ϕ
≤B
2
Z 4
r dr dϕ ≤ C B2 h3 .
Since | ϕ1 − ϕ2 | is bounded from below by c1 h, there are at most 2 π/(c1 h) = C/h triangles in Th,1 . This shows
∑
Z
4∈Th,1 4
∂vh 2 1 ∂v (r, ϕ) − dr dϕ ≤ C B2 h2 . r ∂ϕ ∂ϕ
Case 2: 0 = r1 < r2 By (3.3b) we obtain ∂vh ∂v ≤ (r1 , ϕ1 ) + C B h ≤ C B h ∂ϕ ∂ϕ and as in Case 1
∂v (r, ϕ) ≤ B r ≤ C B h. ∂ϕ
16
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Let us denote by w(s) the length of the intersection of the triangle with some line r = s. We have Z Z r3 Z ∂vh 2 1 1 1 ∂v 2 2 2 2 dr dϕ = C B h w(r ) dr ≤ C B2 h3 . (r, ϕ) − dr dϕ ≤ C B h r ∂ϕ ∂ϕ r r 4 0 4 Here we used w(r ) ≤ C r and r3 ≤ c2 h, which follow from the assumptions on the mesh family. Since all triangles in Th,2 lie in the strip [0, c2 h] × [0, 2 π ], there are at most (c2 h 2 π )/( β h2 ) = C/h triangles in this case, where β > 0 is chosen such that the area of any triangle is bounded from below by β h2 . This yields Z 1 ∂v ∂vh 2 ( r, ϕ ) − dr dϕ ≤ C B2 h2 . ∑ 4 r ∂ϕ ∂ϕ 4∈T h,2
Case 3: 0 < r1 Using (3.3b), we get Z Z 1 ∂v ∂vh 2 C ∂vh 2 C B2 h4 ∂v . (r, ϕ) − dr dϕ ≤ (r, ϕ) − dr dϕ ≤ ∂ϕ r1 4 ∂ϕ ∂ϕ r1 4 r ∂ϕ For each triangle in Th,3 , there must be another triangle in the mesh which has (r1 , ϕ1 ) as a vertex and lies to the left of (r1 , ϕ1 ). More precisely, this triangle contains a nontrivial portion of the line segment connecting (r1 , ϕ1 ) with the boundary point (0, ϕ1 ). Due to the assumptions on the mesh regularity, this implies r1 ≥ α h for some α > 0 which depends only on c1 , c2 . Now, the set Th,3 is further divided into subsets, based on the value of r1 . In particular, we set n Th,3 = 4 ∈ Th,3 : r1 ∈ [α n h, α (n + 1) h) Sb
R
c
n αh for n = 1, . . . , b αRh c. Since b αRh c > αRh − 1, we get Th,3 = n= 1 T h,3 . Moreover, since n all triangles in Th,3 lie in the strip [α n h, α (n + 1) h + c2 h] × [0, 2 π ], there are at most n . We obtain (α + c2 ) h 2 π/( β h2 ) = C/h triangles in Th,3 Z ∂vh 2 C C B2 h4 C B2 h2 1 ∂v ( r, ϕ ) − dr dϕ ≤ · ≤ . ∑ r ∂ϕ ∂ϕ h αnh n 4∈T n 4 h,3
By using an upper bound for the harmonic series b αRh c
∑
n =1
R 1 ≤ 1 + lnb c ≤ C + |ln h| n αh
we obtain
∑
Z
4∈Th,3 4
b αRh c Z ∂vh 2 ∂vh 2 1 ∂v 1 ∂v (r, ϕ) − dr dϕ ≤ ∑ ∑ (r, ϕ) − dr dϕ r ∂ϕ ∂ϕ r ∂ϕ ∂ϕ n=1 4∈T n 4 h,3
≤ C B2 h2
b αRh c
∑
n =1
17
1 ≤ C B2 h2 (1 + |ln h|). n
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
In all cases, a Taylor estimate yields Z
∂v ∂v 2 r (r, ϕ) − h dr dϕ ≤ C B2 h4 ∂r ∂r 4
Summing over all triangles yields the claim. As already discussed above, this result implies the convergence of the discretization scheme of the forward problem as h → 0. Based on this, the convergence of an optimal solution in a discretize–then optimize setting can be deduced as well using standard arguments. In the sequel, we will however follow an optimize–then discretize approach for which the discussion of convergence is more involved and beyond the scope of this paper.
3.3 D ISCRETIZATION AND S OLUTION OF THE U NCONSTRAINED O PTIMALITY S YSTEM To complete the discretization of (P), we follow an optimize–then discretize approach. In this section, we consider the case without control constraints. We employ piecewise linear controls on the same triangular grid which we use for the state. Therefore, the right hand side in the discrete state equation (3.1) can be realized by the term M u, where M is mass matrix [ ∼ ,M ] = assema (p ,t ,0 , char ( ’x ’) ,0) ; We do not impose additional conditions on the discrete control (as we did for the state) because its continuity will follow automatically from the discrete optimality system given below. To preserve the iterated structure of the term k p¯ (r, ·)k L2 (µ ϕ ) in (2.8), we choose a particular quadrature formula and impose further structural conditions on the mesh. From now on, the mesh vertices are supposed to form a rectangular lattice. This supersedes the conditions set forth in Section 3.1. For simplicity, we elaborate on the case of constant mesh widths hr and h ϕ . It will be convenient to address the components of the coefficient vector u by a double index. The value of the control at (r, ϕ) = (i hr , j h ϕ ) is denoted by uij with 1 ≤ i ≤ nr and 1 ≤ j ≤ n ϕ . e , which are exThe discrete state and adjoint state are represented by vectors y e and p panded to the full nodal basis by multiplication with Z. Our discrete optimality system now consists of the discrete state equation, Z> K Z y e = ZM u,
y = Zy e,
(3.4)
the discrete adjoint equation, e = −ZM (y − yd ), Z> K Z p
18
e, p = Zp
(3.5)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
and the following discretization of (2.8), β
!
α uij = max 0, 1 − n ϕ 1/2 pij =: G(p) ij pij = G(p) p ij . ∑ ωk p2ik
(3.6)
k =1
Here the weights ωk are equal to h ϕ /2 for k ∈ {1, n ϕ } and ωk = h ϕ otherwise. The vector yd represents the coefficient of a nodal interpolation of the desired state yd . In the last term in (3.6), G(p) denotes a vector and is the pointwise product between vectors of the same size. Note that (3.6) implies that the discrete control inherits the continuity of the discrete adjoint state, which in turn follows from the fact that its coefficient vector p is in the null space of the constraint matrix N, and that the factor in parentheses depends only on the index i in radial direction. We solve (3.4)–(3.6) by a semismooth Newton method. We mention that the matrix in the resulting Newton system for the update step > δe y 0 Z> K Z Z MZ 0 e ) (3.7) αI diag(p) G0 (p) Z + diag(G(p)) Z δu = −R(y e, u, p δe p Z> K Z −Z> M 0 turns out to be non-symmetric, and not obviously symmetrizable. Here G0 (p) denotes a Newton derivative of G, and the right hand side is the negative residual of (3.5), (3.6) and (3.4) in this order. Numerical experiments revealed that globalization efforts were not required for convergence.
3.4 C ASE OF P OINTWISE C ONTROL C ONSTRAINTS The numerical treatment of (P) is more involved in the presence of pointwise inequality constraints for the control. The reason is that these constraints are not compatible with the sparsity structure induced by the iterated norm in the objective. This was observed already in Herzog et al. [2012]. Following Lemma 4.1 in that paper, we reformulate the optimality system (2.6) as a non-smooth equation. One can show that the following two statements are equivalent: 1. λ¯ ∈ ∂k·k 1,2 (u¯ ) and (2.6) holds. Lr ( Ω )
2. (2.6a) and (2.6c) and
− p¯ + α u¯ + β λ¯ + µ¯ = 0 max 1, kλ¯ (r, ·) + c1 u¯ (r, ·)k L2 (µ ϕ ) λ¯ − (λ¯ + c1 u¯ ) = 0 µ¯ − max 0, µ¯ + c2 (u¯ − ub ) − min 0, µ¯ + c2 (u¯ − u a ) = 0 hold for any choice of positive constants c1 , c2 .
19
(3.8a) (3.8b) (3.8c)
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
We mention that (3.8b) is equivalent to λ¯ ∈ ∂k·k L1,2 (Ω ) (u¯ ), and (3.8c) is equivalent to µ¯ r ¯ This can be shown by a straightforward being an element of the normal cone of Uad at u. distinction of cases. Unfortunately, both versions of the optimality system have their drawbacks. The first version is not directly amenable to numerical treatment due to the variational inequality in (2.6b). For the second version, the viability of the semismooth Newton method in function space is unknown. Nevertheless, a semismooth Newton type iteration can be successfully applied in the discrete setting. To this end, we discretize the optimality system in a similar way as before in Section 3.3. The additional variables λ¯ and µ¯ are of the same dimension as the control u¯ in the discrete setting. The term involving the k·k L2 (µ ϕ ) in (3.8b) is treated in the same way as in (2.8). The derivation of the discrete semismooth Newton system proceeds in the same way as in the unconstrained case, and we refer to [Herzog et al., 2012, Section 4] for details. As was observed already there, two modifications should be applied to the Newton system in order to ensure its well-posedness in each step as well as the global convergence of the method in practice. Firstly, note that (3.8b) implies u¯ = 0 on the subset of [0, R] where the ’max’ attains the value one. On the other hand, (3.8c) implies u¯ ∈ {u a , ub } on some subsets of Ω , which are termed active sets. This may lead to contradictory conditions at intermediate iterations, and a singular Newton matrix ensues. Therefore, we give preference to (3.8b) and modify the determination of the active sets in (3.8c) so that they can only be subsets of {(r, ϕ) : kλ¯ (r, ·) + c1 u¯ (r, ·)k L2 (µ ϕ ) < 1}. The second modification concerns the linearization of (3.8b), to which a damping is applied. We refer to [Herzog et al., 2012, eq. (4.10)] for details. We point out that both modifications vanish in the limit and they do not impair fast local convergence. In Figure 3.2 we show the optimal control for the sectorial problem (with and without control constraints) obtained for the following setting, α = 0.01, u a = −1,
β = 0.15, yd ( x, y) = e2 x sin(y π ), ub = 1.
(3.9)
In the complementarity formulation (3.8), we used the parameters α c1 = and c2 = 100. β
4 S ECTORIAL F ORMULATION In this section we point out which changes to (P) are necessary to obtain optimal controls with sectorial sparsity patterns. As can be expected, we have to modify the sparsity promoting term kuk L1,2 (Ω ) in the objective. To this end, we define the space r
def 1 2 L˜ 1,2 r ( Ω ) = L ( µ ϕ ; L ( µr )).
20
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Figure 3.2: Optimal control of the annular sparsity problem (P) without (left plot) and with control constraints (right plot). The parameters are given in (3.9). 1 2 Compared with L1,2 r ( Ω ) = L ( µr ; L ( µ ϕ )) we have changed the order of integration. The norm in the space L˜ 1,2 r ( Ω ) is given by
kuk L˜ 1,2 r (Ω
)
=
Z 2π Z R 0
0
2
|u(r, ϕ)| dµr
1/2
dµ ϕ =
Z 2π Z R 0
0
|u(r, ϕ)|2 r dµ
1/2
dµ ϕ .
The optimal control problem now reads Minimize such that and
1 α ky − yd k2L2r (Ω ) + kuk2L2r (Ω 2 2 (y, u) satisfy (1.9)
)
+ β kuk L˜ 1,2 r (Ω
)
(P)
u a ≤ u ≤ ub a.e. on Ω .
Following the analysis in Section 2, we have to identify the dual space of L˜ 1,2 r ( Ω ) and 2 the subdifferential of its norm in Lr (Ω ). Standard arguments yield that the dual space def is L˜ r∞,2 (Ω ) = L∞ (µ ϕ ; L2 (µr )), and using Lemma 2.1, we have n
∂k·k L˜ 1,2 (Ω ) (v) = w ∈ r
L2r (Ω
) : kwk L˜ r∞,2 (Ω
)
≤ 1 and
Z Ω
v w d(µr × µ ϕ ) = kvk L˜ 1,2 (Ω r
o )
(4.1) Analogously to the computations in the proof of Lemma 2.2, we find that, for any given u ∈ L2r (Ω ), the function λ ∈ L2r (Ω ) belongs to the subdifferential ∂k·k L˜ 1,2 (Ω ) (u) if r and only if where u(·, ϕ) ≡ 0, kλ(·, ϕ)k L2 (µr ) ≤ 1 u(·, ϕ) elsewhere λ(·, ϕ) = ku(·, ϕ)k 2 L ( µr ) holds for almost all ϕ ∈ (0, 2π ). The same arguments as used in the proof of Theorem 2.3 show that the optimality system (2.6) is the same as in the annular formulation, except that now λ¯ ∈ ∂k·k L˜ 1,2 (Ω ) (u) r
21
.
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
holds. Following the calculations leading to (2.7), we also obtain the equivalence u¯ (·, ϕ) ≡ 0
⇔
k p¯ (·, ϕ)k L2 (µr ) ≤ β
(4.2)
for almost all ϕ ∈ (0, 2π ), which confirms the occurrence of sectorial sparsity patterns. The analog of formula (2.8) in the absence of control constraints reads α u¯ (r, ϕ) = max 0, 1 −
β p¯ (r, ϕ) k p¯ (·, ϕ)k L2 (µr )
a.e. in Ω .
(4.3)
However, (4.3) does not imply the continuity of the (transformed) optimal control u¯ ◦ P−1 at the origin, even if p¯ ◦ P−1 is continuous. To see this, note that the continuity of p¯ ◦ P−1 implies that p¯ (0, ϕ) is constant, but (4.3) does not imply that u¯ (0, ϕ) is constant w.r.t. ϕ since the max(. . .) term depends on ϕ. The numerical treatment of the sectorial problem requires only a few changes compared to the annular variant. The restrictions we impose on the mesh and the discrete state and adjoint states, as well as the formulation of the discrete state equation are the same as in Section 3.1. Concerning the numerical algorithm, we simply need to replace (2.8) by (4.3) in the unconstrained case and proceed as before. In the constrained case, we merely need to replace (3.8b) by max 1, kλ¯ (·, ϕ) + c1 u¯ (·, ϕ)k L2 (µr ) λ¯ − (λ¯ + c1 u¯ ) = 0 and continue as in Section 3.4. In Figure 4.1 we show the optimal solution for the following data, α = 0.05, u a = −1,
β = 0.02, yd ( x, y) = e2 x sin(y π ), ub = 1.
(4.4)
In the complementarity formulation, we used the parameters c1 =
α β
and
c2 = 10.
5 C ONCLUSIONS AND O UTLOOK We introduced and analyzed convex optimal control problems whose optimal controls feature particular spatial sparsity patterns. Optimality systems were derived, which led to semismooth Newton type methods for the numerical solution. We elaborated on the cases of annular and sectorial sparsity patterns on circular domains by way of the polar coordinate transform and an iterated norm in the objective. The extension to domains which are described by rectangles ( R1 , R2 ) × ( ϕ1 , ϕ2 ) in polar coordinates (annular sectors) is straightforward. The conditions set forth in Section 3.1
22
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Figure 4.1: Optimal control of the sectorial sparsity problem (P) without (left plot) and with control constraints (right plot). The parameters are given in (4.4). to achieve conformity need to be adjusted. In fact, there is no need to use only rectangles as Ω . And finally, the polar coordinate transformation could be replaced by a general diffeomorphism to accommodate domains and sparsity patterns of rather arbitrary shape. In three dimensions, the variety of choices by which the three coordinates can be grouped into outer (sparse) and inner coordinates increases from two to six. For example, when we use spherical coordinates and place the radial coordinate outside and the two angular coordinates inside, the support of an optimal control will consist of spherical shells.
R EFERENCES E. Casas, C. Clason, and K. Kunisch. Approximation of elliptic control problems in measure spaces with sparse solutions. SIAM Journal on Control and Optimization, 50 (4):1735–1752, 2012a. ISSN 0363-0129. doi: 10.1137/110843216. E. Casas, R. Herzog, and G. Wachsmuth. Optimality conditions and error analysis of semilinear elliptic control problems with L1 cost functional. SIAM Journal on Optimization, 22(3):795–820, 2012b. doi: 10.1137/110834366. E. Casas, R. Herzog, and G. Wachsmuth. Approximation of sparse controls in semilinear equations by piecewise linear functions. Numerische Mathematik, 122(4):645–669, 2012c. doi: 10.1007/s00211-012-0475-7. E. Casas, C. Clason, and K. Kunisch. Parabolic control problems in measure spaces with sparse solutions. SIAM Journal on Control and Optimization, 51(1):28–63, 2013. ISSN 0363-0129. doi: 10.1137/120872395. C. Clason and K. Kunisch. A duality-based approach to elliptic control problems in non-reflexive Banach spaces. ESAIM: Control, Optimisation, and Calculus of Variations, 17(1):243–266, 2011. doi: 10.1051/cocv/2010003.
23
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
C. Clason and K. Kunisch. A measure space approach to optimal source placement. Computational Optimization and Applications. An International Journal, 53(1):155–171, 2012. ISSN 0926-6003. doi: 10.1007/s10589-011-9444-9. S.J. Coorpender, D. Finkel, J. Kyzar, R. Sims, A.B. Smirnova, M. Tawhid, C.E. Bouton, and R.C. Smith. Modeling and optimization issues concerning a circular piezoelectric actuator design. Technical Report CRSC–TR99–22, North Carolina State University, 1999. URL http://www.ncsu.edu/crsc/reports/reports99.html. J. Diestel and J.J. Uhl. Vector Measures. Mathematical Surveys and Monographs. American Mathematical Society, Providence, 1977. I. Ekeland and R. Temam. Convex Analysis and Variational Problems, volume 28 of Classics in Applied Mathematics. SIAM, Philadelphia, 1999. D. H. Fremlin. Measure theory. Vol. 2. Torres Fremlin, Colchester, 2003. ISBN 0-95381292-8. Broad foundations, Corrected second printing of the 2001 original. R. Herzog, G. Stadler, and G. Wachsmuth. Directional sparsity in optimal control of partial differential equations. SIAM Journal on Control and Optimization, 50(2):943– 963, 2012. doi: 10.1137/100815037. A. D. Ioffe and V. M. Tichomirov. Theorie der Extremalaufgaben. VEB Deutscher Verlag der Wissenschaften, Berlin, 1979. X. Li, H. Du, L. Xu, Y. Hu, and L. Xu. Optimization of a circular thin-film piezoelectric actuator lying on a clamped multilayered elastic plate. IEEE Transactions on Ultrasonics, 56(7):1469–1475, 2009. A. Raghavan and C. E. S. Cesnik. Modeling of piezoelectric-based lamb wave generation and sensing for structural health monitoring. In Smart Structures and Materials 2004: Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems, volume 5391, pages 419–430, 2004. doi: 10.1117/12.540269. G. Stadler. Elliptic optimal control problems with L1 -control cost and applications for the placement of control devices. Computational Optimization and Applications, 44(2): 159–181, 2009. doi: 10.1007/s10589-007-9150-9. J.L. van Niekerk, B.H. Tongue, and A.K. Packard. Active control of a circular plate to reduce transient noise transmission. Journal of Sound and Vibration, 183(4):643 – 662, 1995. doi: 10.1006/jsvi.1995.0277. G. Wachsmuth and D. Wachsmuth. Convergence and regularization results for optimal control problems with sparsity functional. ESAIM: Control, Optimisation and Calculus of Variations, 17(3):858–886, 2011. doi: 10.1051/cocv/2010027. C. M. Yeum, H. Sohn, and J. B. Ihn. Lamb wave mode decomposition using concentric ring and circular piezoelectric transducers. Wave Motion. An International Journal
24
Annular and Sectorial Sparsity in Optimal Control
Herzog, Obermeier, Wachsmuth
Reporting Research on Wave Phenomena, 48(4):358–370, 2011. ISSN 0165-2125. doi: 10.1016/j.wavemoti.2011.01.001.
25
