An approach to optimal actuator design based on shape and topology optimisation techniques is presented. For linear diffusion equations, two scenarios are considered. For the first one, best actuators are determined depending on a given initial condition. In the second scenario, optimal actuators are determined based on all initial conditions not exceeding a chosen norm. Shape and topological sensitivities of these cost functionals are determined. A numerical algorithm for optimal actuator design based on the sensitivities and a level-set method is presented. Numerical results support the proposed methodology.

shape optimisation, feedback control, topological derivative, shape derivative, level-set method

49Q10, 49M05, 93B40, 65D99, 93C20.

§ INTRODUCTION

In engineering, an actuator is a device transforming an external signal into a relevant form of energy for the system in which it is embedded. Actuators can be mechanical, electrical, hydraulic, or magnetic, and are fundamental in the control loop, as they materialise the control action within the physical system. Driven by the need to improve the performance of a control setting, actuator/sensor positioning and design is an important task in modern control engineering which also constitutes a challenging mathematical topic. Optimal actuator positioning and design departs from the standard control design problem where the actuator configuration is known a priori, and addresses a higher hierarchy problem, namely, the optimisation of the control to state map.

There is no unique framework which is followed to address optimal actuator problems. However, concepts which immediately suggest themselves -at least for linear dynamics- and which have been addressed in the literature, build on choosing actuator design in such a manner that stabilization or controllability are optimized by an appropriate choice of the controller. This can involve Riccati equations from linear-quadratic regulator theory, and appropriately chosen parameterizations of the set of admissible actuators. The present work partially relates to this stream as we optimise the actuator design based on the performance of the resulting control loop. Within this framework,  we follow a distinctly different approach by casting the optimal actuator design problem as  shape and  topology optimisation problems.   The class of admissible actuators are characteristic functions of measurable sets and their shape is determined by techniques from shape calculus and optimal control. The class of cost functionals which we consider within this work are quadratic ones and account for the stabilization of the closed-loop dynamics. We present the concepts here for the linear heat equation, but the techniques  can be extended to more general classes of functionals and stabilizable dynamical systems. We believe that the concepts of shape and topology optimisation constitute an important tool for solving  actuator positioning problems, and to our knowledge this can be the first step towards this direction. More concretely, our contributions in this paper are:

i) We study an optimal actuator design problem for linear diffusion equations. In our setting, actuators are parametrised as indicator functions over a subdomain, and are evaluated according to the resulting closed-loop performance for a given initial condition, or among a set of admissible initial conditions not exceeding a certain norm.

ii) By borrowing a leaf from shape calculus, we derive shape and topological sensitivities for the optimal actuator design problem.

iii) Based on the formulas obtained in ii),  we construct a gradient-based and a level-set method for the numerical realisation of optimal actuators.

iv) We present a numerical validation of the proposed computational methodology. Most notably, our numerical experiments indicate that throughout the proposed framework we obtain non-trivial, multi-component actuators, which would be otherwise difficult to forecast based on tuning, heuristics, or experts' knowledge.

Let us, very briefly comment on the related literature. Most of these endeavors focus on control problems related to ordinary differential equations.  We quote the  two  surveys papers [12, 27] and [26]. From these publications already it becomes clear that the notion by which optimality is measured  is an important topic in its own right. The literature on optimal actuator positioning for distributed parameter systems is less rich but it also dates back for several decades already. From among the earlier contributions we quote [9] where the topic is investigated in a semigroup setting for linear systems, [5] for a class of linear infinite dimensional filtering problems, and [11] where the optimal actuator problem is investigated for hyperbolic problems related to active noise suppression. In the works [18, 16, 19] the optimal actuator problem is formulated in terms of parameter-dependent linear quadratic regulator problems where the parameters characterize the position of  actuators, with predetermined shape,  for example.  By choosing the actuator position in [13] the authors optimise the decay rate in the one-dimensional wave equation. Our research may be most closely related to the recent contribution [21], where the optimal actuator design is driven by exact controllability considerations, leading to actuators which are chosen on the basis of minimal energy controls steering the system to zero within a specified time uniformly, for a bounded set of initial conditions. Finally, let us mention that the optimal actuator problem is in some sense dual to optimal sensor location problems [14], which is of paramount importance.

Structure of the paper
The paper is organised as follows.

In Section <ref>,   the optimal control problems, with respect to which optimal actuators are sought later, are introduced. While the first formulation depends on a single initial condition for the system dynamics,  in the second formulation the optimal actuator  mitigates the worst closed-loop performance among all the possible initial conditions.

In Sections <ref> and <ref> we derive the shape and topological sensitivities associated to the aforedescribed optimal actuator design problems.

Section <ref> is devoted to describing a numerical approach which constructs the optimal actuator based on the shape and topological derivatives computed in Sections <ref> and <ref>. It involves the numerical realisation of the sensitivities and iterative gradient-based and level-set approaches.

Finally in  Section <ref> we report on  computations involving numerical tests for our model problem in dimensions one and two.

§.§ Notation

Let $\DD\subset \R^d$, $d=1,2,3$ be either a bounded domain with $C^{1,1}$ boundary $\partial \Omega$ or a convex domain, and let $T>0$ be a fixed time. The space-time cylinder is denoted by $\Omega_T := \Omega\times (0,T]$.    Further  by $H^1(\DD)$ denotes the  Sobolev space of square integrable functions on $\DD$ with square integrable  weak  derivative. The space $H^1_0(\DD)$ comprises all functions in $H^1(\DD)$ that have trace zero on $\partial \DD$ and $H^{-1}(\DD)$ stands for the dual of $H^1_0(\DD)$. The space $\ac C^{0,1}(\overbar\DD , \R^d)$ comprises all Lipschitz
continuous functions on $\overbar{\DD}$ vanishing on $\partial \DD $. It is
a closed subspace of $C^{0,1}(\overbar\DD ,\R^d)$, the space of Lipschitz continuous mappings
defined on $\overbar\DD $.
Similarly we denote by $\ac C^k(\overbar\DD ,\R^d)$ all $k$-times differentiable
functions on $\overbar\DD $ vanishing on $\partial \DD $. We use the notation $\partial f$ for the
Jacobian of a function $f$. Further
$B_\eps(x)$ stands for the open ball centered at $x\in \R^d$ with radius $\eps >0$. Its closure is denoted $\overbar{B}_\eps(x) := \overbar{B_\eps(x)}$.
By $\mathfrak Y(\DD)$ we denote the set of all measurable subsets $\omega\subset \DD$. We say that a sequence $(\omega_n)$ in $\mathfrak Y(\DD)$ converges to an element $\omega\in \mathfrak Y(\DD)$ if
$\chi_{\omega_n}\to \chi_\omega $ in $L_1(\DD)$ as $n\to \infty$, where $\chi_{\omega}$ denotes the characteristic function of $\omega$. In this case we write $\omega_n \to \omega$.  Notice that $\chi_{\omega_n}\to \chi_\omega $ in $L_1(\DD)$ as $n\to \infty$ if and only if $\chi_{\omega_n}\to \chi_\omega $ in $L_p(\DD)$ as $n\to \infty$ for all $p\in [1,\infty)$.
For two sets $A,B\subset \R^d$ we write $A\Subset B$ is $\overbar{A}$ is compact and $\overbar{A}\subset B$.

§ PROBLEM FORMULATION AND FIRST PROPERTIES

§.§ Problem formulation

Our goal is to study an optimal actor positioning and design problem for a controlled linear parabolic equation. Let $\mathcal U$ be a  closed and convex subset of $L_2(\DD)$ with $0\in \mathcal U$. For each $\omega\in \mathfrak Y(\DD)$ the set   $\chi_\omega \mathcal U$ is a convex subset of $L_2(\DD)$. The elements of the space $\mathfrak Y(\DD)$ are referred to as actuators. The choices $\mathcal U= L_2(\DD)$ and $\mathcal U = \R$, considered as  the space of constant functions on $\Omega$, will play a special role.
Further, $\UU := L_2(0,T;\mathcal U)$
denotes the space of time-dependent controls, which is equipped with the topology induced by the $L_2(0,T;L_2(\DD))-$norm. We denote  by  $K$ a nonempty, weakly closed subset of $H^1_0(\Omega)$. It will serve as the set of admissible initial conditions for the stable formulation of our optimal actuator positioning problem.

With these preliminaries we consider for every triplet $(\omega,u,\y)\in \mathfrak Y(\DD)\times \UU \times H^1_0(\Omega)$  the following linear parabolic equation: find $y:\overbar \Omega\times [0,T]\to \R$ satisfying

\begin{align}
\partial_t y   - \Delta y =   \chi_\omega u  & \qquad \text{ in } \DD\times (0,T], \label{eq:statea} \\
y=0 & \qquad\text{ on } \partial \DD\times (0,T], \label{eq:stateb}\\
y(0)=f & \qquad\text{ on } \DD. \label{eq:statec}
\end{align}

In the following, we discuss the well-posedness of the system dynamics <ref> and the associated linear-quadratic optimal control problem, to finally state the optimal actuator design problem.

Well-posedness of the linear parabolic problem
It is a classical result <cit.>  that system (<ref>) admits a unique weak solution $y=y^{u,\y,\omega}$ in $W(0,T)$, where

\[
W(0,T) :=  \{ y\in L_2(0,T;H^1_0(\DD)):\; \partial_t y \in L_2(0,T;H^{-1}(\DD)) \},
\]

which satisfies by definition,

⟨∂_t y, φ⟩_H^-1, H^1_0  + ∫_∇y ·∇φ  dx =  ∫_χ_ωu  φ  dx
for all  $\varphi\in H^1_0(\DD)$ for a.e.  $t\in (0,T]$, and $ y(0)=\y$. For the shape calculus of Section 4 we require that $f \in H_0^1 (\Omega)$.
In this case the state variable enjoys additional regularity properties. In fact, in <cit.> it is shown that for  $\y\in H^1_0(\DD)$ the weak solution $y^{\omega,u,\y}$ satisfies

y^u,f,ω∈L_2(0,T,H^2())∩L_∞(0,T;H^1_0()),   ∂_t y^u,f,ω∈L_2(0,T;L_2())
and there is a constant $c>0$, independent of $\omega, \y$ and $u$, such that

y^u,f,ω_L_∞(H^1) + y^u,f,ω_L_2(H^2) + ∂_t y^u,f,ω_L_2(L_2) ≤c(  χ_ωu _L_2(L_2) + f_H^1).
Thanks to the lemma of Aubin-Lions the space

Z(0,T) := {y∈L_2(0,T;H^2(Ω)∩H^1_0(Ω)):  ∂_t y ∈L_2(0,T;L^2 (Ω))  }
is compactly embedded into $L_\infty(0,T;H^1_0(\DD))$.

The linear-quadratic optimal control problem After having discussed the well-posedness of the linear parabolic problem, we recall a standard linear-quadratic optimal control problem associated to a given actuator $\omega$. Let $\gamma >0$ be given. First we define for every triplet $(\omega,f,u)\in \mathfrak Y(\DD)\times  H^1_0(\Omega)\times \UU$ the cost functional

J(ω,u,) :=  ∫_0^Ty^u,f,ω(t)_L_2()^2 + γ u (t)_L_2()^2  dt.
By taking the infimum in (<ref>) over all controls $u\in \UU$ we obtain the function $\mathcal J_1$, which is defined for all $(\omega,f)\in \mathfrak Y(\DD)\times  H^1_0(\Omega)$:

𝒥_1(ω,) :=  inf_u∈  J(ω,u,f).
It is well known, see e.g. [25] that the minimisation problem on the right hand side of (<ref>), constrained to the dynamics (<ref>)  admits a unique solution. As a result, the  function $\mathcal J_1(\omega,\y)$ is well-defined.  The minimiser  $\overbar u$ of (<ref>) depends on the initial condition $f$ and the set $\omega$, i.e., $\overbar u = \overbar u^{\omega,\y}$.
In order to eliminate the dependence of the optimal actuator $\omega$ on the initial condition $f$ we define a robust  function $\mathcal J_2$ by  taking the supremum in (<ref>) over all normalized initial conditions $\y$ in $K$:

𝒥_2(ω) := sup_∈K,
We show later on that the supremum on the right hand side of (<ref>) is actually attained.

The optimal actuator design problem We now have all the ingredients to state the optimal actuator design problem we shall study in the present work.
In the subsequent sections we are concerned with the following minimisation problem

|ω|  = c 𝒥_1(ω,f),  for  f ∈K,

where $c\in (0,|\DD|)$ is the measure of the prescribed volume of the actuator $\omega$. That is, for a given initial condition $f$ and a given volume constraint $c$, we design the actuator $\omega$ according to the closed-loop performance of the resulting linear-quadratic control problem (<ref>). Note that no further constraint concerning the actuator topology is considered. Buidling upon this problem, we shall also study the problem

|ω|  = c 𝒥_2(ω),

where the dependence of the optimal actuator on the initial condition of the dynamics is removed by minimising among the set of all the normalised initial condition $f\in K$.

Finally, another problem of interest which can be studied within the present framework is the  optimal actuator positioning problem, where the topology of the actuator is fixed, and only its position is optimised. Given a fixed set $\omega_0\subset \DD$ we study the optimal actuator positioning problem by solving

inf_X∈^d 𝒥_1((𝕀+X)(ω_0),f),  for  f ∈K,

inf_X∈^d 𝒥_2((𝕀+X)(ω_0)),

where $(\id +X)(\omega_0) = \{x+X:\; x\in \omega_0\}$, i.e., we restrict our optimisation procedure to a set of actuator translations.

Our goal is to characterize shape and topological derivatives for $\mathcal J_1(\omega,f)$ (for fixed $f$) and $\mathcal J_2(\omega)$ in order to develop gradient type algorithms to solve (<ref>) and (<ref>). The results presented in Sections <ref> and <ref> can  also be utilized to derive optimality conditions for problems (<ref>) and (<ref>). In addition, we investigate numerically whether the proposed methodology provides results which coincide with physical intuition.

While the existence of optimal shapes according to  (<ref>) and (<ref>) is certainly also an interesting task,  this issue is postponed to future work. We  mention [21] where a problem similar  to ours but with different cost functional is considered.

§.§ Optimality system for $\mathcal{J}_1$

The unique solution $\bar u\in \UU$ of the minimisation problem on the right hand side of (<ref>) can be characterised by the first order necessary optimality condition

∂_u J(ω,u̅,)(v-u̅) ≥0   for all  v∈.
The function $\bar u\in \UU$ satisfies the variational inequality (<ref>) if and only if there is a multiplier $\overbar p\in W(0,T)$ such that  the triplet
$(\overbar u,\overbar y,\overbar p)\in \UU\times W(0,T)\times W(0,T)$  solves

\begin{align}
\int_{\Omega_T} \partial_t \overbar y \varphi   + \nabla \overbar y \cdot \nabla \varphi \; dx\;dt & =  \int_{\Omega_T}  \chi_\omega \overbar u \varphi \; dx\;dt \quad \text{ for all } \varphi\in W(0,T),  \label{eq:opt_1}\\
\int_{\Omega_T} \partial_t \psi  \overbar p + \nabla \psi \cdot \nabla \overbar p \; dx\;dt & =  -\int_{\Omega_T} 2 \overbar y \psi \; dx\;dt \quad \text{ for all } \psi\in W(0,T), \label{eq:opt_2} \\
\int_{\Omega} (2\gamma \overbar u - \chi_\omega \bar p)& (v-\overbar u)\; dx \ge  0 \quad \text{ for all } v\in \mathcal U,\quad \text{ a.e. } t\in (0,T),  \label{eq:opt_3}
\end{align}

supplemented with the initial and terminal conditions $\overbar y(0)=\y$ and $\overbar p(T)=0$ a.e. in $\DD$.
Two cases are of particular interest to us:

(a) If $\mathcal U = L_2(\DD)$, then (<ref>) is equivalent to $2\gamma \bar u =\chi_\omega \bar p$ a.e. on $\DD\times (0,T)$.

(b) If $\mathcal U = \R$, then (<ref>) is equivalent to
$ 2\gamma \bar u = \int_{\omega} \bar p \; dx$ a.e. on $(0,T)$.

§.§ Well-posedness of $\mathcal{J}_2$

Given $\omega\in \mathfrak Y(\DD)$ and $\y\in K$, we use the notation $\overbar{u}^{\y,\omega}$ to denote the unique minimiser of $J(\omega,\cdot,\y)$ over $\UU$.

Let $(\y_n)$ be a sequence in $K$ that converges weakly in $H^1_0(\Omega)$ to  $f\in K$, let $(\omega_n)$ be a sequence in $\mathfrak Y(\DD)$ that converges to $\omega\in \mathfrak Y(\DD)$, and let $(u_n)$ be a sequence in $\UU$ that converges weakly to a function $u\in \UU$. Then we have

y^u_n,_n,ω_n     →y^u,f,ω    in  L_2(0,T;H^1_0(Ω))    as  n→∞,

y^u_n,_n,ω_n      ⇀y^u,f,ω    in  L_2(0,T;H^2(Ω)∩H^1_0(Ω))   as  n→∞.

The a-priori estimate (<ref>) and the compact embedding $Z(0,T)\subset \\ L_2(0,T;H^1_0(\DD))$ show that we can extract a subsequence of $(y^{u_n,\y_n,\omega_n})$ that converges weakly to an element $y$ in $L_2(0,T;H^2(\Omega)\cap H^1_0(\Omega))$ and strongly in $L_2(0,T;H^1_0(\Omega))$. Using this to pass to the limit in (<ref>) with $(u,\y,\omega)$ replaced by $(u_n,\y_n,\omega_n)$ implies by uniqueness that $y = y^{u,f,w}.$

Let $(\y_n)$ be a sequence in $H^1_0(\Omega)$ converging weakly to $f\in H^1_0(\Omega)$ and let $(\omega_n)$ be a sequence in $\mathfrak Y(\DD)$ that converges to $\omega\in \mathfrak Y(\DD)$. Then we have
u̅^_n,ω_n  →u̅^f,ω    in  L_2(0,T;L_2())  as  n→∞.

Using estimate (<ref>) we see that for all $u\in \UU$ and $n\ge 0$, we have

∫_0^Ty^ u̅^_n,ω_n ,_n,ω_n(t)_L_2()^2 + γ u̅^_n,ω_n (t)_L_2()^2  dt

≤∫_0^Ty^ u ,_n,ω_n(t)_L_2()^2 + γ  u (t)_L_2()^2  dt

≤c( χ_ω_n u _L_2(L_2)^2 + f_n_H^1^2).

It follows that $(\bar u_n):= (\bar u^{\y_n,\omega_n})$ is bounded in $L_2(0,T;L_2(\DD))$ and hence there is an element $\bar u\in L_2(0,T;L_2(\DD))$ and a subsequence $(\bar u_{n_k})$,  $\bar u_{n_k} \rightharpoonup \bar u$ in $L_2(0,T;L_2(\DD))$ as $k\to \infty$. In addition this  subsequence satisfies $\liminf_{k\to \infty} \|\bar u_{n_k}\|_{L_2(0,T;L_2(\DD))} \ge \|\bar u\|_{L_2(0,T;L_2(\DD))}$. Since $\mathcal U$ is closed we also have $\bar u \in L_2(0,T;\mathcal U)$. Together with Lemma <ref> we therefore obtain from (<ref>) by taking the $\liminf$ on both sides,
∫_0^Ty^ u̅ ,,ω(t)_L_2()^2 + γ  u̅ (t)_L_2()^2  dt ≤∫_0^Ty^ u ,,ω(t)_L_2()^2 + γ  u (t)_L_2()^2  dt
for all $u\in \UU$. This shows that $\bar u=\bar u^{\y,\omega}$ and since $\bar u^{\y,\omega}$ is the unique minimiser of $J(\omega,\cdot ,y)$ the whole sequence $(\bar u_n)$ converges weakly to $\bar u^{\y,\omega}$. In addition it follows from the strong convergence $y^{ \bar u^{\y_n,\omega_n} ,\y_n,\omega} \to y^{ \bar u^{\y,\omega} ,\y,\omega}$ in $W(0,T)$ and estimate (<ref>) that the norm $\|\bar u^{\y_n,\omega_n}\|_{L_2(0,T;L_2(\DD))}$ converges to $\|\bar u^{\y,\omega}\|_{L_2(0,T;L_2(\DD))}$. As norm convergence together with weak convergence imply strong convergence, this shows that $\bar u^{\y_n,\omega_n}$ converges strongly to $\bar u^{\y,\omega}$ in $L_2(0,T;L_2(\DD))$ as was to be shown.

We now prove that $\omega\mapsto \mathcal J_2(\omega)$ is well-defined on $\mathfrak Y(\DD)$.

For every $\omega\in \mathfrak Y(\DD)$ there exists $ \y \in K$ satisfying $\|\y\|_{H^1_0(\Omega)}\le 1$ and
𝒥_2(ω) = 𝒥_1(ω, ).

Let $\omega\in \mathfrak Y(\DD)$ be fixed. In view of  $0\in \mathcal U$ and (<ref>) and since $K\subset H^1_0(\Omega)\hookrightarrow H^1_0(\DD)$ we obtain for all $\y \in H^1_0(\Omega)$ with $\|\y\|_{H^1_0(\Omega)}\le 1$,
𝒥_1(ω,) = min_u∈J(ω,u,) ≤∫_0^Ty^0,f,ω(t)_L_2()^2 dt ≤c _H^1_0(Ω)^2 ≤cr^2.
Further we can express $\mathcal J_2$ as follows
𝒥_2(ω) = sup_∈K
_H^1_0(Ω)≤1∫_0^Ty^ u̅^,ω ,,ω(t)_L_2()^2 + γ  u̅^,ω (t)_L_2()^2 dt.
Let $(\y_n)\subset K$, $\|\y_n\|_{H^1_0(\Omega)} \le 1$ be a maximising sequence, that is,

𝒥_2(ω) = lim_n→∞∫_0^Ty^ u̅^ω,_n ,_n,ω(t)_L_2()^2 + γ u̅^ω,_n (t)_L_2()^2  dt.
The sequence $(\y_n)$ is bounded in $K$ and therefore we find a subsequence $(\y_{n_k})$ converging weakly to an element $\y\in K$. Additionally,  the limit element satisfies  $\|\y\|_{H^1_0(\Omega)} \le \liminf_{k\to \infty}\|\y_{n_k}\|_{H^1_0(\Omega)}\le 1$  and hence $\|\y\|_{H^1_0(\Omega)}\le 1$. Since $(\y_{n_k})$ is also bounded in $H^1_0(\DD)$ we may assume that $(\y_{n_k})$ also converges weakly to  $\y\in H^1_0(\DD)$. Thanks to  Lemmas <ref> and <ref> we obtain

𝒥_2(ω)    = lim_k→∞∫_0^Ty^ u̅^_n_k,ω ,_n_k,ω(t)_L_2()^2 + γ  u̅^_n_k,ω (t)_L_2()^2  dt

= ∫_0^Ty^ u̅^,ω ,,ω(t)_L_2()^2 + γ u̅^,ω(t)_L_2()^2  dt.

In view of Lemma <ref> we write from now on
$ \mathcal J_2(\omega) =$$\max_{\substack{\y\in K,\\\|\y\|_{H^1_0(\Omega)}\le 1}}\mathcal J_1(\omega,f).$

§ SHAPE DERIVATIVE

In this section we prove the directional differentiability of
$\mathcal J_2$ at arbitrary measurable sets. We employ the averaged adjoint approach [23] which is tailored to the
derivation of directional derivatives of PDE constrained shape functions. Moreover this approach  allows us later on to also compute the topological derivative of  $\mathcal J_1$ and $\mathcal J_2$ without performing asymptotic analysis which can otherwise  be quite involved [20].

Of course, there are notable alternative approaches, most prominent the material derivative approach,  to prove directional differentiability of shape functions, see e.g. [15, 6]. For an overview of available methods the reader may consult [24].

§.§ Shape derivative

Given a vector field $X\in \ac C^{0,1}(\overbar\DD ,\R^d)$, we denote by $T_t^X$ the perturbation of the identity
T_t^X(x) := x + t X(x)
which is bi-Lipschitz for all $t\in [0,\tau_X]$, where
\tau_X := 1/(2\|X\|_{C^{0,1}}).
We omit the index $X$ and write $T_t$ insteand of $T_t^X$ whenever no confusion is possible. A mapping $J:\mathfrak Y(\DD) \rightarrow \R $  is called shape function.

The directional derivative of $J$ at $\omega\in \mathfrak Y(\DD)$ in direction $X\in \ac C^{0,1}(\overbar\DD , \R^d)$ is defined by
DJ(ω)(X):= lim_t ↘0J(T_t(ω))-J(ω)/t .
We say that $J$ is

(i) directionally differentiable at $
\omega$ (in $\ac C^{0,1}(\overbar\DD , \R^d)$), if 	 $DJ(\omega)(X)$
exists for all
$X\in C^{0,1}(\overbar\DD , \R^d)$,

(ii)  differentiable at $\omega$
(in $\ac C^{0,1}(\overbar\DD ,\R^d)$), if  $DJ(\omega)(X)$ exists for all
$X \in \ac C^{0,1}(\overbar\DD ,\R^d)$ and
$ X     \mapsto DJ(\omega)(X) $
is linear and continuous.

The following properties will frequently be used.

Let $\DD \subseteq \R^d$ be open and bounded and pick a vector field  $X\in \ac C^{0,1}(\overbar\DD, \R^d)$. (Note that $T_t(\Omega)=\Omega$ for all $t$.)

(i) We  have as $t\rightarrow 0^+$,

\begin{align*}
\frac{ \partial T_t - I}{t} \rightarrow  \partial X \quad
&\text{ and }\quad   \frac{
\partial T_t^{-1} - I}{t} \rightarrow  - \partial X &&
\text{ strongly in } L_\infty(\overbar\DD , \R^{d\times d})\\
\frac{ \det(\partial T_t) - 1}{t} \rightarrow & \divv(X) && \text{ strongly in } L_\infty(\overbar\DD ).
\end{align*}

(ii)  	For all
$\varphi \in L_2(\Omega)$, we have as $t\rightarrow 0^+$,

\begin{align}\label{eq:Lp}
\varphi\circ T_t  \rightarrow & \varphi
&& \text{ strongly in }  L_2(\Omega).
\end{align}

Let $(\varphi_n)$ be a sequence in $H^1(\Omega)$ that converges weakly to $\varphi \in H^1(\Omega)$. Let $(t_n)$ a null-sequence. Then we have as $n \to \infty$,

\begin{align}\label{eq:Wkp}
\frac{\varphi_n\circ T_{t_n}  - \varphi_n}{t_n} \rightharpoonup & \nabla \varphi \cdot X
&& \text{ weakly in }  L_2(\Omega).
\end{align}

Item (i) is obvious. The convergence result (<ref>) is proved in <cit.> and (<ref>) can be proved in a similar fashion.

Item (iii) is less obvious and we give a proof. For every $\eps >0$ and $\psi\in H^1(\Omega)$, there is $N>0$, such that $|(\varphi_n-\varphi,\psi)_{H^1}|\le \eps$ for all $n\ge N_\epsilon$.  By density we find for every $n$ and every null-sequence $(\epsilon_n)$,  $\eps_n >0$ an element $\tilde \varphi_n \in C^1(\overbar\Omega)$, such that
φ̃_n - φ_n_H^1 ≤_n.
It is clear that $\tilde \varphi_n \rightharpoonup \varphi$  weakly in $H^1(\Omega)$ as $n\to \infty$.
We now write

φ_n ∘T_t_n - φ_n/t_n - ∇φ_n ·X  =     (φ_n - φ̃_n) ∘T_t_n - (φ_n-φ̃_n)/t_n - ∇(φ_n-φ̃_n) ·X

+ φ̃_n ∘T_t_n - φ̃_n/t_n - ∇φ̃_n ·X.

Let $x\in \Omega$. Applying the fundamental theorem of calculus to $s\mapsto \tilde \varphi_n(T_{s}(x))$ on $[0,1]$ gives
φ̃_n(T_t_n(x)) - φ̃_n(x)/t_n =  ∫_0^1 ∇φ̃_n(x+t_n s X(x)) ·X(x) ds.
We now show that the function $q_n(x) := \int_0^1 \nabla \tilde \varphi_n(x+t_n s X(x)) \cdot X(x)$ converges weakly to $\nabla \varphi\cdot X$ in $L_2(\Omega)$. For this purpose we consider for $\psi\in L_2(\Omega)$,
∫_Ωq_n ψ  dx = ∫_Ω∫_0^1 ∇φ̃_n(x+t_n s X(x)) ·X(x) ψ(x) ds   dx.
Interchanging the order of integration and invoking a change of variables (recall $T_t(\Omega)=\Omega$), we get
∫_Ωq_n ψ  dx  = ∫_0^1 ∫_Ω(∂T_st_n^-1)  ∇φ̃_n ·((Xψ)∘T_st_n^-1)   dx_:=η(t_n, s)  ds.
Owing to item (ii) and noting that $X\circ T_t^{-1} \to X$ in $L_\infty(\Omega)$ as $t\to 0$, we also have for $s\in [0,1]$ fixed,
(∂T_st_n^-1)   (Xψ)∘T_st_n^-1  →Xψ   in  L_2(Ω, ^2)    as  n →∞.
As a result using the weak convergence of $(\tilde \varphi_n)$ in $H^1(\Omega)$,  we get for $s\in[0,1]$,
η(t_n, s)   →∫_Ω∇φ·X ψ  dx   as  n→∞.
It is also readily checked using Hölder's inequality that $|\eta(t_n, s)| \le c \|\nabla \tilde \varphi_n\|_{L_2}\|\psi\|_{L_2}$ for a constant $c>0$ independent of $s\in [0,1]$. As a result we may apply Lebegue's dominated convergence theorem to obtain
∫_Ω q_n ψ  dx = 	∫_0^1 η(t_n, s)   ds →∫_0^1 η(0, s)  ds =  ∫_Ω∇φ·X  dx    as  n→∞.
This proves that $q_n$ converges weakly to $\nabla \varphi\cdot X$.

Finally testing (<ref>) with $\psi$, integrating over $\Omega$ and estimating gives

|     (φ_n ∘T_t_n  - φ_n/t_n - ∇φ_n ·X,ψ)_L_2|

≤cψ_L_2 (_n/t_n + _n)
+ |( φ̃_n ∘T_t_n - φ̃_n/t_n - ∇φ̃_n ·X, ψ)_L_2|

with a constant $c>0$ only depending on $X$. Now we choose $\tilde N_\epsilon \ge 1$ so large that
|( φ̃_n ∘T_t_n - φ̃_n/t_n - ∇φ·X, ψ)_L_2| ≤ϵ   for all  n≥Ñ_ϵ.

|( φ̃_n ∘T_t_n - φ̃_n/t_n - ∇φ̃_n ·X, ψ)_L_2|

≤ϵ+ |(∇(φ̃_n - φ_n)·X,ψ)_L_2| + |(∇(φ_n - φ)·X, ψ)_L_2|

≤ϵ+ ϵ_n + ϵ   for all  n≥max{N_ϵ, Ñ_ϵ}.

Choosing $\epsilon_n := \min\{ t_n^2,\epsilon \}$ and combining the previous estimate with (<ref>) shows the right hand side of (<ref>) can be bounded by $3\epsilon$. Since $\epsilon >0$ was arbitrary we see that (<ref>) holds.

§.§ First main result: the directional derivative of $\mathcal J_2$

Given $\omega\in \mathfrak Y(\DD)$ and $r>0$, we define the set of maximisers of  $\mathcal J_1(\omega,\cdot)$ by

𝔛_2(ω) := {f̅∈K:  sup_f∈K,
f_H^1_0(Ω)≤1 𝒥_1(ω,f) = 𝒥_1(ω,f̅)}.
The set $\mathfrak X_2(\omega)$ is nonempty as shown in Lemma <ref>. Before stating our first main result we make the following assumption.

For every $X\in \ac C^{0,1}(\overbar \DD, \R^d)$ and $t\in [0,\tau_X]$ we have
u∈𝒰   ⟺  u∘T_t ∈𝒰.

Assumption <ref> is satisfied for $\mathcal U$ equal to  ${\tiny }L_2(\DD)$ or $\R$.

Under the Assumption <ref> we have the following theorem, where
we set $\bar y^{\y,\omega} := y^{\bar u^{\omega,\y},\y,\omega}$ and $\bar p^{\y,\omega} := p^{\bar u^{\omega,\y},\y,\omega}$ for $\omega\in \mathfrak Y(\DD)$ and $\y\in K$.
Furthermore we define for $A\in \R^{d\times d}, B\in \R^{d\times d}, a, b, c\in \R^{d}$

\begin{equation*}
A:B= \sum_{i,j=1}^d a_{ij}b_{ij}, \quad (a\otimes b)c := (b\cdot c) a,
\end{equation*}

where $a_{ij},b_{ij}$ are the entries of the matrices $A,B$, respectively.

The directional derivative of $\mathcal J_2(\cdot)$ at $\omega$ in direction $X\in \ac C^{0,1}(\overbar \DD,\R^d)$ is given by

D𝒥_2(ω)(X) = max_∈𝔛_2(ω)∫_Ω_T _1(y̅^,ω ,p̅^,ω,u̅^,ω):∂X + _0() ·X   dx dt,
where the functions $\Sb_1(\y) := \Sb_1(\bar y^{\y,\omega} ,\bar p^{\y,\omega},\bar u^{\y,\omega})$ and $\Sb_0(\y)$ are given by

_1()  =    I( |y̅^,ω|^2 + γ|u̅^,ω|^2 -  y̅^,ω ∂_t p̅^,ω  + ∇y̅^,ω·∇p̅^,ω  - χ_ωu̅^,ω p̅^,ω)

- ∇y̅^,ω⊗∇p̅^,ω - ∇p̅^,ω⊗∇y̅^,ω,

_0()  =    - 1/T∇ p̅^,ω

and the adjoint $\bar p^{\y,\omega}$ satisfies

∂_t p̅^,ω  - Δp̅^,ω =  -2 y̅^,ω                in  ×(0,T],

p̅^,ω=0                on ∂×(0,T],

p̅^,ω(T)=0                in  .

(b) The directional derivative of $\mathcal J_1(\cdot,\y)$ at $\omega$ in direction $X\in \ac C^{0,1}(\overbar \DD,\R^d)$ is given by

D𝒥_1(ω,)(X) =∫_Ω_T _1():∂X + _0()·X   dx dt,
where $\Sb_0(\y)$ and $\Sb_1(\y)$ are defined by  (<ref>).

We notice that for $r>0$ we have
_H^1_0(Ω)≤r𝒥_1(ω,f) = 	r^2max_∈1/r K,
Therefore we may assume that $\bar \y\in K$ with $\|\bar \y\|_{H^1_0(\Omega)}\le 1$.
Setting $K:= \{\bar \y \}$, we have for all  $\omega\in \mathfrak Y(\DD)$,
𝒥_2(ω)  = max_∈K,
_H^1_0(Ω)≤1𝒥_1(ω,f) = 𝒥_1(ω,f̅)
and hence the result follows from item (a) since $\mathfrak X_2(\omega)=\{\bar f\}$ is a singleton. The proof of part (a) will be given in the following subsections.

We pause here to comment on the regularity requirements imposed on $\y$. As can be seen from the volume expression (<ref>) we can extend $D\mathcal J_1(\omega, \y)$ to initial conditions $\y $ in $L_2(\DD)$. In fact, the only term that requires weakly differentiable initial conditions is the one involving $\Sb_0$ and it can be rewritten as follows for a.e. $t\in [0,T]$,

∫_ _0(t)·X  dx     = - 1/T ∫_
∇·X p̅^, ω(t)   dx

= 1/T ∫_ (X) p̅^, ω(t) + ∇p̅^, ω(t) ·X  dx,

where we used that $\bar p^{\y, \omega}(t)=0$ on $\partial \DD$. This shows that the shape derivative $D\mathcal J_1(\omega, \y)$ can be extended to initial conditions $\y \in L_2(\DD)$. However, it is not possible to obtain the shape derivative for $\y \in L_2(\DD)$ in general. This will become clear in the proof of Theorem <ref>.

The next corollary shows that under certain smoothness assumptions on $\omega$ we can write the integrals (<ref>) and (<ref>) as integrals over $\partial \omega$.

Let $\y\in K$ and  $X\in \ac C^{0,1}(\overbar \DD,\R^d)$ be given. Assume that $\omega\Subset \Omega$ and $\Omega$ are  $C^2$ domains. Moreover, suppose that either $\mathcal U = L_2(\DD)$ or $\mathcal U = \R$.

(a)  Given $\y \in \mathfrak X_2(\omega)$ define  $\hat{\Sb_1}(\y) := \int_{0}^{T} \Sb_1(\y)(s)\;ds$ and
$\hat {\Sb_0}(\y) := \int_{0}^{T} \Sb_0(\y)(s)\;ds$. Then we have

_̂1̂()|_ω∈W^1_1(ω,^d×d),       _̂1̂()|_∖ω∈W^1_1(∖ω,^d×d),

-(_̂1̂())+ _̂0̂() =0    a.e. in  ω∪(∖ω).
Moreover (<ref>) can be written as

D𝒥_2(ω)(X)     = max_∈𝔛_2(ω)∫_∂ω [_̂1̂()ν] ·X  ds
= max_∈𝔛_2(ω)-∫_∂ω ∫_0^Tu̅^ω, p̅^ω, (X·ν) dt  ds

for  $X\in \ac C^1(\overbar \DD,\R^d)$, with $\nu$ the outer normal to $\omega.$ Here $[\hat{\Sb_1}(\y)\nu] :=\\ \hat{\Sb_1}(\y)|_{\omega}\nu - \hat{ \Sb_1}(\y)|_{\DD\setminus\omega}\nu$ denotes the jump of $\hat{\Sb_1}(\y)\nu$ across $\partial \omega$.

(b) We have that  (<ref>) can be written as

D𝒥_1(ω,)(X) = -∫_∂ω ∫_0^Tu̅^ω, p̅^ω, (X·ν) dt  ds
for $X\in \ac C^1(\overbar \DD, \R^d)$.

Before we prove this corollary we need the following auxiliary result.

Suppose that $\DD$ is of class $C^2$. For all $\y \in {H^1_0(\Omega)}$ and $\omega\in \mathfrak Y(\DD)$, we have

∫_0^Ty̅^,ω(t) ∂_t p̅^,ω(t)  dt ∈W^1_1(Ω),    and    ∫_0^T ∇p̅^,ω(t)·∇y̅^,ω(t)   dt ∈W^1_1(Ω).

From the general regularity results <cit.> we have that $\bar p^{\y,\omega}\in L_2(0,T;H^3(\DD))$ and $\partial_t \bar p^{\y,\omega}\in L_2(0,T;H^1(\DD))$, and $\bar y^{\y,\omega}\in L_2(0,T;H^2(\DD))$ and $\partial_t \bar y^{\y,\omega} \in L_2(0,T;L_2(\DD))$.

Observe that for almost all $t\in [0,T]$ we have $\partial_t \bar p^{\y,\omega}(t)\in H^1(\DD)$ and $\bar y^{\y,\omega}(t) \in H^2(\Omega) $. So
since $H^1(\DD)\subset L_6(\Omega)$ and $H^2(\DD)\subset C(\overbar\DD)$, where we use that $\Omega\subset \R^d$, $d\le 3$  we also have $\bar y^{\y,\omega}(t)\partial_t \bar p^{\y,\omega}(t)\in L_6(\DD)$ and a.e. $t \in (0,T)$

y̅^,ω(t)∂_t p̅^,ω(t)_L_1() ≤Cy̅^,ω(t)_H^2() ∂_t p̅^,ω(t)_H^1()
for an constant $C>0$. Moreover by the product rule we have
∂_x_j(y̅^,ω(t)∂_t p̅^,ω(t)) = ∂_x_j(y̅^,ω(t))_∈H^1() ∂_t p̅^,ω(t)_∈H^1() + y̅^,ω(t)_∈H^1() (∂_x_j∂_t p̅^,ω(t))_∈L_2(),
so that $\partial_{x_j}(\bar y^{\y,\omega}(t)\partial_t \bar p^{\y,\omega}(t))\in L_1(\DD)$ and

∂_x_j(y̅^,ω(t)∂_t p̅^,ω(t))_L_1() ≤C y̅^,ω(t)_H^1() ∂_tp̅^,ω(t)_H^1()
for some constant $C>0$. So (<ref>) and (<ref>) imply that $t\mapsto \|\bar y^{\y,\omega}(t)\partial_t \bar p^{\y,\omega}(t)\|_{W^1_1(\DD)}$ belongs to $L_1(0,T)$.
This shows the left inclusion in (<ref>). As for the right hand side inclusion in (<ref>) notice that for almost all $t\in [0,T]$ we have $\bar p^{\y,\omega}(t) \in H^3(\DD)$. Therefore $\nabla \bar p^{\y,\omega}(t)\in H^2(\DD)$ and $\nabla \bar y^{\y,\omega}(t)\in H^1(\DD)$ and thus $\nabla \bar y^{\y,\omega}(t)\cdot \nabla \bar p^{\y,\omega}(t)\in L_6(\DD)$. Similarly we check that $\partial_{x_j}(\nabla \bar y^{\y,\omega}(t)\cdot \nabla \bar p^{\y,\omega}(t))\in L_1(\DD)$ and thus $t\mapsto \|\nabla \bar y^{\y,\omega}(t)\cdot \nabla \bar p^{\y,\omega}(t)\|_{W^1_1(\DD)}\in L_1(0,T)$, which gives the right hand side inclusion in (<ref>).

We assume that Theorem <ref> holds. As a consequence of Lemma <ref> we obtain (<ref>).
Then for all $X \in C^1_c( \DD,\R^d)$ satisfying $X|_{\partial \omega}=0$ we have $T_t(\omega) = (\id+tX)(\omega)=\omega$ for all $t\in [0,\tau_X]$. Hence $D\mathcal J_2(\omega)(X)=0$ for such vector fields which gives

0=D𝒥_2(ω)(X) ≥∫__1():∂X + _0() ·X   dx
for all $X \in C^1_c( \DD,\R^d)$ satisfying $X|_{\partial \omega}=0$
and for all $\y\in \mathfrak X_2(\omega)$. Since for fixed $\y$ the expression in (<ref>) is linear in $X$  this proves

∫__1():∂X + _0() ·X   dx=0
for all $X \in C^1_c( \DD,\R^d)$ satisfying $X|_{\partial \omega}=0$
and for all $\y\in \mathfrak X_2(\omega)$.  Hence testing of (<ref>) with vector fields $X\in C^1_c( \omega,\R^d)$ and $X\in C^1_c( \DD\setminus \overbar \omega,\R^d)$, partial integration  and
(<ref>) yield the continuity equation (<ref>). As a result, by partial integration (see e.g. [17]), we get for all $X \in C^1_c( \DD,\R^d)$,

D𝒥_2(ω)(X)     = max_∈𝔛_2(ω)∫__1():∂X + _0() ·X   dx

=max_∈𝔛_2(ω)(∫_∂ω [_1()ν]·X   ds  + ∫_ω (-(_1() +  _0())_=0·X   dx

+ ∫_∖ω (-(_1() +  _0())_=0·X   dx),

which proves the first equality in  (<ref>).  Now using Lemma <ref> we see that  $\Tb(\y) := \hat{\Sb}_1(\y) + \int_0^T\chi_\omega \bar u^{\y,\omega}(t) \bar p^{\y,\omega}(t)\;dt$ belongs to $W^1_1(\DD,\R^{d\times d})$ and hence $[\Tb(\y)\nu]=0$ on $\partial \omega$. It follows that $[\hat \Sb_1(\y)\nu] = -\int_0^T\chi_\omega \bar u^{\y,\omega}(t) \bar p^{\y,\omega}(t)\;dt $ which finishes the proof of (a). Part (b) is a direct consequence of part (a).

The following observation is important for our gradient algorithm that we introduce later on.

Let the hypotheses of Theorem <ref> be satisfied. Assume that if $v\in \mathcal U$ then $-v\in \mathcal U$. Then we have
D𝒥_1(ω,-)(X) = D𝒥_1(ω,)(X)
for all $ X\in \ac C^{0,1}(\overbar \DD, \R^d)$ and $\y\in {H^1_0(\Omega)}$.

Let $\y\in {H^1_0(\Omega)}$ be given. From the optimality system (<ref>) and the assumption that $v\in \mathcal U$ implies $-v\in \mathcal U$, we infer that $\overbar u^{-\y,\omega} = -\overbar u^{\y,\omega}$, $\bar y^{-\y,\omega} =- \bar y^{\y,\omega}$ and $\bar p^{-\y,\omega} =- \bar p^{\y,\omega}$. Therefore $\Sb_1(-\y)=\Sb_1(\y)$ and $\Sb_0(-\y) = \Sb_0(\y)$ and the result follows from (<ref>).

The following sections are devoted to the proof of Theorem <ref>(a) .

§.§ Sensitivity analysis of the state equation

In this paragraph we study the sensitivity of the solution $y$ of (<ref>) with respect to $(\omega,\y,u)$.

Perturbed state equation
Let $X\in \ac C^{0,1}(\overbar \DD, \R^d)$ be a vector field and define $T_\tau:= \id+ \tau  X$.  Given $u\in  U$, $\y\in {H^1_0(\Omega)}$ and $\omega\in \mathfrak Y(\DD)$,  we consider (<ref>) with $\omega_\tau := T_\tau(\omega)$,

∂_t y^u,,ω_τ  - Δy^u,,ω_τ =  χ_ω_τ u                in  ×(0,T],

y^u,,ω_τ=0                on ∂×(0,T],

y^u,,ω_τ(0)=               in  .

We define the new variable

y^u,,τ := (y^u∘T^-1,,ω_τ) ∘T_τ.
Then since $\chi_{\omega_\tau} = \chi_\omega\circ T_\tau^{-1}$ and $\Delta f \circ T_t = \divv(A(t)\nabla(f\circ T_t))$, it follows from (<ref>)-(<ref>) that

∂_t y^u,,τ   - 1/ξ(τ)(A(τ)∇y^u,,τ) =  χ_ωu                in  ×(0,T],

y^u,,τ=0                on ∂×(0,T],

y^u,,τ(0)=∘T_τ               in  ,

\[
A(\tau) := \det(\partial T_\tau) \partial T_\tau^{-1} \partial T_\tau^{-\top}, \qquad \xi(\tau) := |\det(\partial T_\tau)|.
\]

Equations (<ref>)-(<ref>) have to be understood in the variational sense, i.e., $y^{u,f,\tau}\in W(0,T)$ satisfying $y^{u,f,\tau}(0)= \y\circ T_\tau$ and

∫__T ξ(τ)∂_t y^u,f,τφ+ A(τ)∇y^u,f,τ·∇φ dx dt = ∫__T ξ(τ)χ_ω uφ dx dt

for all $\varphi \in W(0,T)$. Since $X\in \ac C^{0,1}(\overbar \DD, \R^d)$, we have for fixed $\tau$,

$$A(\tau,\cdot),\partial_\tau A(\tau,\cdot)\in L_\infty(\DD,\R^{d\times d}), \quad \xi(\tau,\cdot),\partial_\tau \xi(\tau,\cdot)\in L_\infty(\DD).$$

Moreover, there are constants $c_1,c_2>0$, such that
A(τ,x)ζ·ζ≥c_1|ζ|^2    for all  ζ∈^d,   for a.e  x∈,     for all  τ∈[0,τ_X]
ξ(τ,x)≥c_2    for a.e  x∈,     for all  τ∈[0,τ_X].

Apriori estimates and continuity

There is a constant $c>0$, such that for all $(u,\y,\omega) \in  \UU\times H^1_0(\DD)\times \mathfrak Y(\DD)$, and  $\tau \in [0,\tau_X]$, we have

y^u,f,ω_τ_L_∞(H^1) + y^u,f,ω_τ_L_2(H^2) + ∂_t y^u,f,ω_τ_L_2(L_2)

≤c( χ_ω_τ u_L_2(L_2) + f_H^1),

y^u,f,τ_L_∞(H^1)  + ∂_t y^u,f,τ_L_2(L_2) ≤c(χ_ωu_L_2(L_2) + f_H^1).

Estimate (<ref>) is a direct consequence of (<ref>).
Let us prove (<ref>). Recalling   $y^{u,f,\tau} = y^{u\circ T^{-1}_\tau,f,\omega_\tau} \circ T_\tau$, a change of variables shows,

∫_Ω_T |y^u,f,τ|^2+|∇y^u,f,τ|^2   dx dt

=∫_Ω_T ξ^-1(τ)|y^u∘T^-1_τ,f,ω_τ|^2 + A^-1(τ)∇y^u∘T^-1_τ,f,ω_τ·∇y^u∘T^-1_τ,f,ω_τ  dx dt

≤c ∫_Ω_T|y^u∘T^-1_τ,f,ω_τ|^2+ |∇y^u∘T^-1_τ,f,ω_τ|^2  dx dt

(<ref>)≤	c( χ_ω_τ u∘T^-1_τ_L_2(L_2) + f_H^1)

≤C(χ_ωu_L_2(L_2)) + f_H^1),

and we further have

χ_ω_τ u∘T^-1_τ_L_2(L_2)^2    = √(ξ)χ_ωu_L_2(L_2)^2 ≤cχ_ωu^2_L_2(L_2).

Combining (<ref>) and (<ref>) we obtain $\|y^{u,f,\tau}\|_{L_2(H^1)} \le c(\|\chi_\omega u\|_{L_2(L_2)} + \|f\|_{H^1})$.
In a similar fashion we can show (<ref>).

An estimate for the second derivatives of $y^{u,f,\tau}$ of the form
y^u,f,τ_L_2(H^2) ≤c( u_L_2(L_2) + f_H^1)
may be achieved by invoking a change of variables in the term $\|y^{u,f}_\tau\|_{L_2(H^2)}$ in (<ref>).
This, however, requires the vector field $X$ to be more regular, e.g.,  $\ac C^2(\overbar \DD, \R^d)$, and is not needed below.

After proving apriori estimates we are ready to derive continuity results for the mapping $(u,f,\tau)\mapsto y^{u,f,\tau}$.

For every $(\omega_1,u_1,f_1),(\omega_2,u_2,f_2)\in \mathfrak Y(\DD)\times \UU\times {H^1_0(\Omega)}$, we denote by  $y_1$ and $y_2$ the corresponding solution of (<ref>)-(<ref>). Then there is a constant $c>0$, independent of $(\omega_1,u_1,f_1),(\omega_2,u_2,f_2)$, such that

y_1-y_2_L_∞(H^1) + y_1-y_2_L_2(H^2) + ∂_t y_1-∂_t  y_2_L_2(L_2)

≤c(χ_ω_1u_1 - χ_ω_2u_2_L_2(L_2)+f_1-f_2_H^1) .

The difference $\tilde y := y_1 - y_1$ satisfies in a variational sense

∂_t ỹ   - Δỹ = u_1χ_ω_1 - u_2χ_ω_2                in  ×(0,T],

ỹ=0                on ∂×(0,T],

ỹ(0)=f_1-f_2                on  .

Hence estimate (<ref>) follows from (<ref>).

As an immediate consequence of Lemma <ref> we obtain the following result.

Let $\omega\in \mathfrak Y(\DD)$ be given. For all $\tau_n\in(0,\tau_X]$, $u_n,u\in \UU$ and $f_n,f\in H^1(\Omega_0)$ satisfying
u_n ⇀u    in  L_2(0,T;L_2()),   f_n ⇀f    in  H^1_0(Ω),   τ_n→0,    as  n→∞,
we have

y^u_n,f_n,τ_n *⇀     y^u,f, ω    in  L_∞(0,T;H^1_0())    as  n→∞,

y^u_n,f_n,τ_n ⇀    y^u,f, ω    in   H^1(0,T;L_2())   as  n→∞.

Thanks to the apriori estimates of Lemma <ref> there exists $y\in \\ L_\infty(0,T;H^1_0(\DD)) \cap H^1(0,T;L_2(\DD))$ and a subsequence $(y^{u_{n_k},f_{n_k},\tau_{n_k}}) $ converging

weakly-star in $L_\infty(0,T;H^1_0(\DD))$ and weakly in  $H^1(0,T;L_2(\DD))$ to $y$.  Since $H^1(\DD)$ embeds compactly into $L^2(\Omega)$ we may assume, extracting another subsequence, that $f_{n_k}\rightarrow f$ in $L_2(\DD)$ as $k\to \infty$. By definition
$y_k := y^{u_{n_k},f_{n_k},\tau_{n_k}}$ satisfies for $k\ge 0$,

∫__T ξ(τ_n_k)∂_t y_kφ+ A(τ_n_k)∇y_k·∇φ dx dt = ∫__T ξ(τ_n_k) χ_ω u_n_kφ dx dt             ,

for all $\varphi \in W(0,T)$, and $y_k(0)=\y_{n_k}\circ T_{\tau_{n_k}}$ on $\DD$. Using the weak convergence of $u_{n_k},y_k$ stated before and the strong convergence obtained using Lemma <ref>,
ξ(τ_n) →1    in  L_∞(),     A(τ_n) →I    in  L_∞(,^d×d),
we may pass to the limit in (<ref>) to obtain,

∫__T ∂_t yφ+ ∇y·∇φ dx dt = ∫__T χ_ω uφ dx dt                 for all  φ∈W(0,T).

Using Lemma <ref> we see
$f_{n_k}\circ T_{\tau_{n_k}}\rightarrow f$ in $L_2(\DD)$ as $k\rightarrow \infty$, and therefore $y(0)=f$.  Since the previous equation with $y(0)=\y$ admits a unique solution we conclude that  $y=y^{u,f,\omega}$. As a consequence of the uniqueness of the limit, the whole sequence $y^{u_n,f_n,\tau_n}$ converges to
$y^{u,f,\omega}$. This finishes the proof.

§.§ Sensitivity of minimisers and maximisers

Let us denote for $(\tau,\y)\in [0,\tau_X]\times K$ the  minimiser of  $u\mapsto J(\omega_\tau,u \circ T_\tau^{-1},\y)$,
by $\bar u^{\y_n,\tau_n}$.

For every null-sequence $(\tau_n)$ in  $[0,\tau_X]$ and every sequence $(\y_n)$ in
$K$ converging weakly (in ${H^1_0(\Omega)}$) to $\y\in K$, we have
u̅^_n,τ_n →u̅^,ω    in  L_2(0,T;L_2())    as  n→∞.

We set $\omega_n := \omega_{\tau_n}$. By definition we have $\bar u^{\y_n,\tau_n} = \bar u^{\y_n,\omega_{\tau_n}}\circ T_{\tau_n}$. From Lemma <ref> we know that
$ \bar u^{\y_n,\omega_{\tau_n}} $ converges to $\bar u^{\y_n,\omega}$ in $L_2(0,T;L_2(\DD))$. Therefore according to Lemma <ref> also $\bar u^{\y_n,\omega_{\tau_n}} \circ T_{\tau_n}$ converges in $L_2(0,T;L_2(\DD))$ to $\bar u^{\y_n,\omega}$. This finishes the proof.

For every  null-sequence $(\tau_n)$ in  $[0,\tau_X]$ and every sequence $(\y_n)$, $\y_n\in \mathfrak{X}_2(\omega_{\tau_n})$, there is a subsequence $(\y_{n_k})$ and $\y\in  \mathfrak{X}_2(\omega)$, such that $\y_{n_k} \rightharpoonup \y$ in ${H^1_0(\Omega)}$ as $k\to \infty$.

We proceed similarly as in the proof of Lemma <ref>. Let $\tau \in [0,\tau_X]$ and $v\in \UU$ be given. We obtain for all $\y \in K$,

J(ω_τ,u^,τ∘T^-1_τ,)  = inf_u∈ J(ω_τ,u∘T^-1_τ,) ≤J(ω_τ, v∘T^-1_τ,).
Let $(\bar\y_n)$ be an arbitrary sequence with $ \bar\y_n \in \mathfrak{X}_2(\omega_{\tau_n})$. Since $\|\bar\y_n\|_{H^1_0(\Omega)}\le 1$ for all $n\ge 0$, there
is a subsequence $(\bar\y_{n_k})$ and a function $\bar\y\in K$, such that $\bar\y_{n_k}\rightharpoonup \bar\y$ in ${H^1_0(\Omega)}$ as $k\to \infty$ and $\|\bar\y\|_{H^1_0(\Omega)}\le 1$. Thanks to Lemma <ref> the sequence $(\bar u_k)$ defined by  $\bar u_k := \bar u^{\bar\y_{n_k},\tau_{n_k}}$ converges to $\bar u^{\bar\y,\omega}$ in $L_2(0,T;L_2(\DD))$. Moreover, Lemma <ref> also shows that $y^{\bar u_k,\bar\y_{n_k},\tau_{n_k}} \to y^{\bar u^{\bar\y,\omega},\bar\y,\omega}$ in $L_2(0,T;L_2(\DD))$. By definition for all
$k\ge 0$ and $ \y\in K$,

∫__T |y^u̅^,τ_n_k,,τ_n_k(t)|^2  + γ| u̅^,τ_n_k(t)|^2  dx dt

f_H^1_0(Ω)≤1 ∫__T |y^u̅^,τ_n_k,,τ_n_k(t)|^2  + γ| u̅^,τ_n_k(t)|^2  dx dt

= ∫__T |y^u̅_k,_n_k,τ_n_k(t)|^2  + γ|u̅_k(t)|^2  dx dt

and therefore passing to the limit $k\to \infty$ yields, for all $\y\in K$,

∫__T |y^u̅^,ω,,ω(t)|^2  + γ| u̅^,ω(t)|^2  dx dt     ≤∫__T |y^u̅^,ω,,ω(t)|^2  + γ| u̅^,ω(t)|^2  dx dt.

This shows that $\y\in \mathfrak{X}_2(\omega)$ and finishes the proof.

§.§ Averaged adjoint equation and Lagrangian

For fixed $\tau\in [0,\tau_X]$ the mapping  $\varphi \mapsto T_\tau^{-1}\circ \varphi$ is an isomorphism on  $\UU$,  therefore,

min_u∈ J(ω_τ,u,)  =  min_u∈ J(ω_τ,u∘T_τ^-1,).
Hence a change of variables shows,

inf_u∈J(ω_τ,u,)     = inf_u∈∫_0^T y^u,,ω_τ(t)_L_2()^2  + γ u(t)_L_2()^2  dt

(<ref>)= inf_u∈∫_Ω_T ξ(τ)(|y^u, , τ(t)|^2  + γ|u(t)|^2 ) dx   dt.

Introduce for every quadruple $(u,\y,y,p)\in \UU \times K\times W(0,T)\times W(0,T)$ and for every $\tau\in [0,\tau_X]$ the parametrised Lagrangian

G̃(τ,u,,y,p)  :=      ∫__T ξ(τ)(|y|^2  + γ| u|^2  ) dx dt

+ ∫__T ξ(τ)  ∂_t y    p  dx dt  + A(τ)∇y ·∇p   dx dt

-  ∫__Tξ(τ)uχ_ωp  dx dt + ∫_ξ(τ)(y(0)-∘T_τ)p(0) dx.

Given $(u,\y)\in \UU\times  K$, and $\tau\in [0,\tau_X]$, the averaged adjoint state $p^{u,f,\tau}\in W(0,T)$ is the solution of averaged adjoint equation

∫_0^1∂_y G̃(τ,u, ,sy^u,f,τ+(1-s)y^u,f,ω,p^u,f,τ)(φ)  ds =0   for all 	 φ∈W(0,T).

The averaged adjoint state $p^{u,f,\tau}$ in our special case only depends on $u$ and $\y$ through
the state $y^{u,f,\tau}$.

It is evident that (<ref>) is equivalent to

∫__T ξ(τ)∂_t φp^u, f, τ + A(τ)∇φ·∇p^u, f,τ  dx dt     + ∫_ξ(τ) p^u, f,τ(0)φ(0) dx

=  -∫__T   ξ(τ)(y^u, f,τ + y^u, f,ω) φ  dx dt

for all $\varphi\in W(0,T)$, or equivalently after partial integration in time

∫__T-ξ(τ) φ∂_t p^u,f,τ  +  A(τ)∇φ·∇p^u,f,τ   dx dt      =  -∫__T  ξ(τ)(y^u,f,τ + y^u,f,ω) φ  dx dt

for all $\varphi\in W(0,T)$, and $p^{u,f,\tau}(T) = 0$.
This is a backward in time linear parabolic equation with terminal condition zero.

§.§ Differentiability of max-min functions

Before we can pass to the proof of Theorem <ref> we need to address a Danskin type theorem on the differentiability of max-min functions.

Let $\mathfrak U$ and $\mathfrak V$ be two nonempty sets and let $ G:[0,\tau]\times \mathfrak U\times \mathfrak V \to \R$ be a function, $\tau >0$. Introduce the function $g:[0,\tau] \to \R$,
g(t) := sup_y∈𝔙 inf_x∈𝔘  G(t,x,y)
and let $\ell:[0,\tau]\rightarrow \R$ be any function such that $\ell(t)>0$ for $t\in (0,\tau]$ and $\ell(0)=0$. We are interested in sufficient conditions that guarantee that the limit
d/dℓ g(0^+) := lim_t↘0^+ g(t)-g(0)/ℓ(0)
exists. Moreover we define for $t\in [0,\tau]$,
𝔙(t) := {y^t∈𝔙: sup_y∈𝔙inf_x∈𝔘  G(t,x,y) = inf_x∈𝔘  G(t,x,y^t) }.

Let the following hypotheses be satisfied.

(A0) For all $y\in \mathfrak V$ and $t\in [0,\tau]$ the minimisation problem
inf_x∈𝔘  G(t,x,y)
admits a unique solution and we denote this solution by $x^{t,y}$.

(A1) For all $t$ in $[0,\tau]$  the set $\mathfrak V(t)$  is nonempty.

The limits

lim_t↘0 G(t,x^t,y,y)-G(0,x^t,y,y)/ℓ(t)

lim_t↘0 G(t,x^0,y,y)-G(0,x^0,y,y)/ℓ(t)
exist for all $y\in \mathfrak U$ and they are equal. We denote the limit by

$\partial_\ell G(0^+,x^{0,y},y)$.

(A3) For all real null-sequences $(t_n)$ in $(0,\tau]$ and all sequences $y^{t_n}$ in $\mathfrak V(t_n)$, there exists a subsequence $(t_{n_k})$ of $(t_n)$,  and $(y^{t_{n_k}})$ of $(y^{t_n})$, and $y^0$ in $\mathfrak V(0)$, such that

lim_k→∞ G( t_n_k ,x^t_n_k,y^t_n_k, y^t_n_k ) -
G(0,x^t_n_k,y^t_n_k, y^t_n_k)/ℓ(t_n_k) = ∂_ℓG(0^+,x^0,y^0,y^0)

lim_k→∞ G( t_n_k ,x^0,y^t_n_k, y^t_n_k ) -
G(0,x^0,y^t_n_k, y^t_n_k)/ℓ(t_n_k) = ∂_ℓG(0^+,x^0,y^0,y^0).

Then we have
d/dℓg(t)|_t=0^+ = max_y∈𝔙(0) ∂_ℓG(0^+,x^0,y,y).

In this section we  apply the previous results  for $\ell(t)=t,$ and in the following one  for $\ell(t) = |B_t(\eta_0)|$, $\eta_0\in \R^d$. For the sake of completeness we give a proof in the appendix; see [8].

§.§ Proof of Theorem <ref>

The following is a direct consequence of (<ref>) and Lemma <ref>.

For all sequences $\tau_n\in(0,\tau_X]$, $u_n,u\in \UU$ and $f_n,f\in K$, such that
u_n ⇀u    in  ,   f_n ⇀f    in  H^1_0(Ω),   τ_n→0,    as  n→∞,
we have

p^u_n,f_n,τ_n →    p^u,f, ω    in  L_2(0,T;H^1_0())    as  n→∞,

p^u_n,f_n,τ_n ⇀    p^u,f, ω    in   H^1(0,T;L_2())   as  n→∞,

where $p^{u,f,\omega}\in Z(0,T)$ solves the adjoint equation

∫__T - φ∂_t p^u,f,ω  dx dt  + ∫__T∇φ·∇p^u,f,ω   dx dt = -∫__T2 y^u,f,ωφ  dx dt

for all $\varphi \in W(0,T)$, and $p^{u,f,\omega}(T) =0$ a.e. on $\DD$.

Now we have gathered all the ingredients to complete the proof of Theorem <ref>(a) on page 9.

Proof of Theorem <ref>(a)
Using the fundamental theorem of calculus we obtain for all $\tau \in [0,\tau_X]$,

G̃(τ,u,,     y^u, f ,τ,p^u,f,τ) -  G̃(τ,u,,y^u, f ,τ,p^u,f,τ)

= ∫_0^1 ∂_y G̃(τ,u,,s y^u, f ,τ + (1-s)y^u, f ,ω,p^u,f,τ)(y^u, f ,τ-y^u, f ,ω)  ds =0,

where in the last step we used the averaged adjoint equation (<ref>). In addition we have
$J(\omega_\tau,u\circ T_\tau^{-1},\y) = \tilde{G}(\tau,u,\y,y^{u, f ,\omega},p^{u,f,\tau}),$ which together with (<ref>) gives
J(ω_τ,u∘T_τ^-1,) = G̃(τ,u,,y^u, f ,ω,p^u,f,τ).
As a consequence we obtain

𝒥_1(ω_τ, ) = inf_u∈ G̃(τ,u,,y^u, f ,ω,p^u,f,τ).

We apply Lemma <ref> with $\ell(t) := t$,
G(τ,u,f) := G̃(τ,u,,y^u, f ,ω,p^u,f,τ),
$\mathfrak U = \UU$, and $\mathfrak V = \{\y\in K:\; \|\y \|_{H^1_0(\Omega)} \le 1 \}$.

Since the minimization problem  (<ref>)  admits a unique solution,  Assumption (A0) is  satisfied.
A  minor change in the proof of  Lemma <ref> to accommodate the reparametrisation of the domain $\omega$ shows that  (A1) is  satisfied as well.

Let $(\tau_n)$ be an arbitrary null-sequence and let $(\y_n)$ be a sequence in $K$ converging weakly in ${H^1_0(\Omega)}$ to $f\in K$,
and let us set $\bar u_n := \bar u^{\y_n,\tau_n}$. Thanks to Lemma <ref> we have that $\bar u_n$ converges strongly in $L_2(0,T;L_2(\DD))$ to $\bar u^{\y,\omega}$.
Moreover Lemma <ref> implies

p^u̅_n,_n,τ_n →    p^u̅^,ω, , ω     in  L_2(0,T;H^1_0())    as  n→∞,

p^u̅_n,_n,τ_n ⇀    p^u̅^,ω, , ω    in   H^1(0,T;L_2())   as  n→∞.

Using Lemma <ref> we see that
A(τ_n) - I/τ_n →(X) - ∂X - ∂X^⊤   in  L_∞(,^d×d)    as  n→∞,
ξ(τ_n)-1/τ_n →(X)    in  L_∞()   as  n→∞.
Therefore we get

G(  τ_n, u̅_n,_n)-G(0,u̅_n,_n)/τ_n

=      G̃(τ_n,u̅_n,_n,y^u̅_n, _n ,ω,p^u̅_n,_n,τ_n) - G̃(0,u̅_n,_n,y^u̅_n, _n ,ω,p^u̅_n,_n,τ_n)/τ_n

=   ∫__T ξ(τ_n)-1/τ(|y^u̅_n, _n ,ω|^2  + γ|u̅_n|^2  ) dx dt

+ ∫__T ξ(τ_n)-1/τ  ∂_t y^u̅_n, _n ,ω    p^u̅_n,_n,τ_n  dx dt

+ ∫__T A(τ_n)-I/τ_n∇y^u̅_n, _n ,ω ·∇p^u̅_n,_n,τ_n   dx dt

-  ∫__Tξ(τ_n)-1/τu̅_n χ_ωp^u̅_n,_n,τ_n  dx dt

+ ∫_ (ξ(τ_n)-1/τ_n  (y^u̅_n, _n ,ω(0) - _n∘T_τ_n) -  _n∘T_τ_n -_n /τ_n )p^u̅_n,_n,τ_n(0) dx

and using Lemma <ref> and (<ref>), we see that the right hand side tends to

∫_Ω_T (X) ( |y̅^,ω|^2 + γ|u̅^, ω|^2 + ∂_t y̅^,ω p̅^,ω  + ∇y̅^,ω ·∇p̅^,ω -   u̅^,ω p̅^,ω χ_ω)  dx  dt

-∫_Ω_T ∂X∇y̅^,ω·∇p̅^,ω + ∂X∇p̅^,ω·∇y̅^,ω + 1/T ∇·X p̅^,ω(0)  dx  dt.

Partial integration in time yields

∫__T  p̅^,ω ∂_t y̅^,ω (X)  dx   dt = - ∫__T  ∂_t p̅^,ω  y̅^,ω (X)  dx   dt - ∫_ (X) p̅^,ω(0)  dx,
where we used $\bar y^{\y,\omega}(0) = \y$ and $\bar p^{\y,\omega}(T)=0$. As a result, inserting (<ref>) into (<ref>), we see that (<ref>) can be written as

∫_Ω_T _1(y̅^,ω ,p̅^,ω,u^,ω):∂X + _0·X   dx dt
with $\Sb_1, \Sb_2$ being given by (<ref>). Hence we obtain

lim_n→∞ G(  τ_n, u̅_n,_n)-G(0,u̅_n,_n)/τ_n = ∫_Ω_T _1(y̅^,ω ,p̅^,ω,u^,ω):∂X + _0·X   dx dt.

Next let  $\bar u_{n,0} := \bar u^{\y_n,0}$. Then we can show in as similar manner as (<ref>) that

lim_n→∞ G(  τ_n, u̅_n,0,_n)-G(0,u̅_n,0,_n)/τ_n = ∫_Ω_T _1(y̅^,ω ,p̅^,ω,u^,ω):∂X + _0·X   dx dt.
Hence choosing $(\y_n)$ to be a constant sequence we see that (A2) is satisfied.

But also (A3) is satisfied  since according to Lemma <ref>  we find for every null-sequence $(\tau_n)$ in  $[0,\tau_X]$ and every sequence $(\y_n)$, $\y_n\in \mathfrak{X}_2(\omega_{\tau_n})$, a subsequence $(\y_{n_k})$ and $\y\in  \mathfrak{X}_2(\omega)$, such that $\y_{n_k} \rightharpoonup \y$ in ${H^1_0(\Omega)}$ as $k\to \infty$.
Now we use (<ref>)  and (<ref>) with $f_n$ replaced by this choice of $f_{n_k}$, and conclude that (A3) holds.
Thus all requirements of Lemma <ref> are satisfied and this ends the proof of Theorem <ref>(a).

§ TOPOLOGICAL DERIVATIVE

In this section we will derive the topological derivative of the shape functions $\mathcal J_1$ and $\mathcal J_2$ introduced in (<ref>) and (<ref>), respectively. The topological derivative, introduced in [22], allows to predict the position where small holes in the shape should be inserted in order to achieve a  decrease of  the shape function.

§.§ Definition of topological derivative

We begin by introducing the so-called topological derivative. For more details we refer to [20].

The topological derivative of a shape funcional $J:\mathfrak Y(\DD) \to \R $ at $\omega\in \mathfrak Y(\DD)$ in the point $\eta_0\in \DD\setminus \partial \omega$ is defined by
𝒯J(ω)(η_0) =  {
lim_↘0J(ω∖B̅_(η_0)) - J(ω)/|B̅_(η_0)|      if  η_0 ∈ω,

lim_↘0J(ω∪B_(η_0)) - J(ω)/|B_(η_0)|      if  η_0 ∈∖ω..

§.§ Second main result: topological derivative of $\mathcal J_2$

Given $\omega\in \mathfrak Y(\DD)$ we set $\omega_\eps := \DD\setminus \bar B(\eta_0)$ if  $\eta_0\in\omega$ and $\omega_\eps := \omega\cup B_\eps(\eta_0)$ if $\eta_0 \in \DD\setminus \overbar \omega$. Denote by $\bar u^{\y,\omega_\eps}$ the minimiser of the right hand side of (<ref>) with $\omega=\omega_\eps$.

Let  $\delta >0$ be so small that $\bar B_\delta(\eta_0)\Subset \DD$. We assume that for all $(\y,\omega)\in V\times \mathfrak Y(\DD)$ we have
$\overbar u^{\y,\omega}\in C(\bar B_\delta(\eta_0))$.
Furthermore we assume that for every sequence $(\omega_n)$ in $\mathfrak Y(\DD)$ converging to $\omega\in \mathfrak Y(\DD)$ and every  weakly converging sequence $\y_n\rightharpoonup \y$ in $V$ we have
lim_n→∞ u^_n,ω_n - u^,ω_L_1(0,T;C(B̅_δ(η_0))) =0.

Lemmas <ref>, <ref> show that Assumption <ref> is satisfied in case $\mathcal U$ is equal to $L_2(\DD)$ or $\R$. Indeed in case $\mathcal U=\R$ we have shown in Remark <ref>,(b) that  $2\gamma \bar u^{\omega,\y}(t) = \int_{\omega} \bar p^{\y,\omega(t,x)} \; dx $, so that $\bar u^{\omega,\y}$ is independent of space and Assumption <ref> is satisfied thanks to  Lemma <ref>. 	In case $\mathcal U=L_2(\DD)$ Remark <ref>,(a) shows that $2\gamma \bar u^{\omega,\y}=\bar p^{\y,\omega}$. In Lemma <ref> below we show that $(\y,\omega)\mapsto \bar p^{\y,\omega}:V\times \mathfrak Y(\DD) \to C([0,T]\times \bar B_\delta(\eta_0))$ is continuous for small $\delta >0$, when $V$ is equipped with the weak convergence we also see that in this case Assumption <ref> is satisfied.

For $\omega\in \mathfrak Y(\DD)$ and $\y\in K$, we set $\bar y^{\y,\omega} := y^{\bar u^{\omega,\y},\y,\omega}$ and  $\bar p^{\y,\omega} := p^{\bar u^{\omega,\y},\y,\omega}$. The main result that we are going to establish reads as follows.

Let $\omega\in \mathfrak Y(\DD)$ be open. Let Assumption <ref> be satisfied at $\eta_0 \in \DD\setminus \partial \omega$.
Then the topological derivative of $\omega\mapsto \mathcal J_2(\omega)$ at $\omega$ in $\eta_0$  is given by

𝒯𝒥_2(ω)(η_0) = max_∈𝔛_2(ω).

-∫_0^Tu^,ω(η_0,s)p̅^,ω(η_0,s) ds      if  η_0 ∈ω,

∫_0^Tu^,ω(η_0,s)p̅^,ω(η_0,s) ds       if   η_0 ∈∖ω,

where the adjoint $\bar p^{\y,\omega}$ belongs to $C([0,T]\times B_\delta(\eta_0))$ and satisfies

∂_t p̅^,ω  - Δp̅^,ω =  -2 y̅^,ω                in  ×(0,T],

p̅^,ω=0                on ∂×(0,T],

p̅^,ω(T)=0                in  .

Let the assumptions of the previous theorem be satisfied. Let $\y\in V$ be given.
Then topological derivative of $\omega\mapsto \mathcal J_1(\omega,\y)$ at $\omega$ in  $\eta_0$ is given by

𝒯𝒥_1(ω,)(η_0) = .

-∫_0^Tu^,ω(x_0,s)p̅^,ω(η_0,s) ds      if  η_0 ∈ω,

∫_0^Tu^,ω(x_0,s)p̅^,ω(η_0,s) ds       if   η_0 ∈∖ω,

where  $\bar p^{\y,\omega}$ solves the adjoint equation (<ref>).

For the same arguments as in proof of Theorem <ref> we may assume that $\bar \y\in K$ with $\|\bar\y \|_V\le 1$.
Setting $K:= \{\bar \y \}$ we obtain for all $\omega\in \mathfrak Y(\DD)$,
𝒥_2(ω)  = max_∈K,
_V≤1𝒥_1(ω,f) = 𝒥_1(ω,f̅)
and hence the result follows from Theorem <ref> since $\mathfrak X_2(\omega)=\{\bar f\}$ is a singleton.

Let the hypotheses of Theorem <ref> be satisfied. Assume that if $v\in \mathcal U$ then $-v\in \mathcal U$. Then we have
𝒯𝒥_1(ω,-)(η_0) = 𝒯𝒥_1(ω,)(η_0)
for all $\eta_0\in \DD\setminus \partial \omega$ and $ \y\in V$.

Let $\y\in V$ be given.  From the optimality system (<ref>) and the assumption that $v\in \mathcal U$ implies $-v\in \mathcal U$, we infer that $\overbar u^{-\y,\omega} = -\overbar u^{\y,\omega}$, $\bar y^{-\y,\omega} =- \bar y^{\y,\omega}$ and $\bar p^{-\y,\omega} =- \bar p^{\y,\omega}$. Now the result follows from (<ref>).

§.§ Averaged adjoint equation and Lagrangian

Throughout this section we fix an open set $\omega\in \mathfrak Y(\DD)$ and pick $\eta_0 \in \omega$. The case $\eta_0 \in \DD\setminus \overbar \omega$ is treated similarly. Let us define $\omega_\eps := \omega\setminus \overbar B_\eps(\eta_0)$, $\eps >0$.

For every quadruple $(u,\y,y,p)\in \UU \times K\times W(0,T)\times W(0,T)$ and every  $\eps \ge 0$ we define the parametrised Lagrangian,

G̃(,u,,y,p)  :=     ∫__Ty^2  +  γu^2   dx   dt
+ ∫__T ∂_t y p + ∇y ·∇p   dx dt

-  ∫__T χ_ω_ u p  dx dt + ∫_(y(0)-∘T_τ)p(0) dx.

We denote by $y^{u,f,\eps}\in W(0,T)$ the solution of the state equation (<ref>) with $\chi=\chi_{\omega_\eps}$ in (<ref>). Then, similarly to  (<ref>), we introduce the averaged adjoint: find $p^{u,f,\eps}\in W(0,T)$, such that
∫_0^1∂_y G̃(,u,,σy^u,f,+(1-σ)y^u,p^u,f,)(φ)  dσ=0   for all 	 φ∈W(0,T)
or equivalently after partial integration in time, $p^{u,f,\eps}(T) = 0$ and

∫_Ω_T -φ∂_t p^u,f,  +  ∇φ·∇p^u,f,   dx dt   =  -∫_Ω_T   (y^u,f, + y^u,f) φ  dx dt
for all $\varphi\in W(0,T)$.

§.§ Proof of Theorem <ref>

Let $\delta>0$ be such that $\bar B_\delta(\eta_0) \Subset \DD$.
For all sequences $\eps_n\in(0,1]$, $u_n,u\in \UU$ and $f_n,f\in K$, such that
u_n ⇀u    in  ,   f_n ⇀f    in  V,   _n→0,    as  n→∞,
we have

p^u_n,f_n,_n →    p^u,f, ω    in  L_2(0,T;H^1_0())    as  n→∞,

p^u_n,f_n,_n ⇀    p^u,f, ω    in   H^1(0,T;L_2())   as  n→∞.

Moreover there is a subsequence $(	p^{u_{n_k},f_{n_k},\eps_{n_k}})$, such that

p^u_n_k,f_n_k,_n_k →p^u,f, ω    in   C([0,T]×B̅_δ(η_0))   as  n→∞.

The first two statements follow by a similar arguments as used in Lemma <ref>. To prove the third we have by interior regularity of parabolic equations   that
p^u,f, ∈Z̃(0,T) :=  L_2(0,T;H^4(B_δ(η_0)))∩H^1(0,T;H^1_0(B_δ(η_0)))∩H^2(0,T;L_2(B_δ(η_0)))
and we have the apriori bound

≤c(y^u,f, + y^u,f_L_2(H^2) + d/dt(y^u,f, + y^u,f)_L_2(L_2)),

see e.g.  <cit.>.
Hence (<ref>) follows since the space $\tilde Z(0,T)$ embeds compactly into $C([0,T]\times \bar B_\delta(\eta_0))$ .

Proof of Theorem <ref>
Proceeding as in the proof of Theorem <ref> we obtain using the averaged adjoint equation,
J(,u,) = G̃(,u,,y^u, f ,ω,p^u,f,)
for $(\eps,u,\y)\in [0,1]\times \UU \times K$, where ${\tilde G}$ is defined in (<ref>).  Hence to prove Theorem <ref> it suffices to apply Lemma <ref> with
G(, u, ) := G̃(,u,,y^u, ,ω,p^u, ,),
$\mathfrak U := \UU$, $\mathfrak V := \{ \y\in K:\; \|\y\|_V\le 1 \}$ and
$\ell(\eps) = |B_\eps(\eta_0)|$.  Since the minimisation problem in (<ref>) is uniquely solvable and in view of Lemma <ref> Assumptions	(A0) and (A1) are satisfied. We turn to verifying (A2) and (A3) next.

Let $(\eps_n)$ be an arbitrary null-sequence and let $(\y_n)$ be a sequence in $K$ converging weakly in $V$ to $\y\in K$.  Thanks to Assumption <ref> the sequence $(\overbar{u}_n)$, $\bar u_n := \bar u^{\y_n, \omega_{\eps_n}}$ converges strongly in $L_1(0,T;C(\bar B_\delta(\eta_0)))$  to $ \overbar{u} =  \overbar{u}^{\y, \omega}\in L_1(0,T;C(\bar B_\delta(\eta_0)))$. Therefore (recall the notation
$\bar p^{\y,\omega_{\eps_n}} = p^{\bar u_n, \y, \omega_{\eps_n}}$) we obtain

G(_n,u̅_n,_n)-G(0,u̅_n,_n)/|B__n(η_0)|  =    -1/|B__n(η_0)|∫_0^T∫_B__n(η_0) u̅_n p̅^_n,_n dx  dt

=    -1/|B__n(η_0)|∫_0^T∫_B__n(η_0) u̅_n  ( p̅^_n,_n - p̅^,ω ) dx  dt

-1/|B__n(η_0)|∫_0^T∫_B__n(η_0) (u̅_n-u̅)  p̅^,ω dx  dt

-1/|B__n(η_0)|∫_0^T∫_B__n(η_0) u̅(x,t)  p̅^,ω(x,t) dx  dt.

Further for all $n$,

1/|B__n(η_0)||∫_0^T∫_B__n(η_0)  (u̅_n-u̅)  p̅^_n,ω dx  dt|
≤ p̅^_n,ω_C([0,T]×B̅_δ(η_0))u̅_n-u̅_L_1(0,T;C(B̅_δ(η_0)))

1/|B__n(η_0)| |∫_0^T∫_B__n(η_0) u̅_n ( p̅^_n,_n - p̅^,ω ) dx  dt|

≤u̅_n_L_1(0,T;C(B̅_δ(η_0))) p̅^_n,_n -  p̅^_n,ω_C([0,T]×B̅_δ(η_0)).

Since $x\mapsto \int_0^T \bar u(x,t)  \bar p^{ \y,\omega}(x,t)\; dt$ is continuous in a neighborhood of $\eta_0$ we also have
lim_n →∞1/|B__n(η_0)|∫_0^T∫_B__n(η_0) u̅(x,t)  p̅^ ,ω(x,t) dx  dt = ∫_0^T u̅(η_0,t)  p̅^,ω(η_0,t) dt.
Hence in view of (<ref>) we obtain

lim_n→∞G(_n,u̅_n,_n)-G(0,u̅_n,_n)/|B__n(η_0)|  =  -∫_0^T u̅(η_0,t)  p̅^,ω(η_0,t) dt

Next let  $\bar u_{n,0} := \bar u^{\y_n,0}$. Then we can show in as similar manner as (<ref>) that

lim_n→∞G(_n,u̅_n,0,_n)-G(0,u̅_n,0,_n)/|B__n(η_0)|  =  -∫_0^T u̅(η_0,t)  p̅^,ω(η_0,t) dt

Hence choosing $(\y_n)$ to be a constant sequence we see that (A2) is satisfied.

But also (A3) is satisfied  since according to Lemma <ref>  we find for every null-sequence $(\tau_n)$ in  $[0,\tau_X]$ and every sequence $(\y_n)$, $\y_n\in \mathfrak{X}_2(\omega_{\tau_n})$, a subsequence $(\y_{n_k})$ and $\y\in  \mathfrak{X}_2(\omega)$, such that $\y_{n_k} \rightharpoonup \y$ in ${H^1_0(\Omega)}$ as $k\to \infty$.
Now we use (<ref>) and (<ref>) with $f_n$ replaced by this choice of $f_{n_k},$ and conclude that (A3) holds.

§ NUMERICAL APPROXIMATION OF THE OPTIMAL SHAPE PROBLEM

In this section we discuss the formulation of numerical methods for optimal positioning and design which are based on the formulae introduced in previous sections. We begin by introducing the discretisation of the system dynamics and the associated linear-quadratic optimal control problem. Then, the optimal actuator design problem is addressed by approximating the shape and topological derivatives, which are embedded into a gradient-based approach and a level-set method, respectively.

§.§ Discretisation and Riccati equation

Let $T>0$.  We choose the spaces $K=H^1_0(\DD)$ and $\mathcal U=\R$, so that the control space $\UU$ is equal to $L_2(0,T;\R)$. The cost functional reads

\begin{align}
\mathcal{J}_1(\omega,\y) = \underset{u\in\UU}{\inf}J(\omega,u,\y)=\int\limits_0^T\|y(t)\|^2_{L^2(\Omega)}+\gamma |u(t)|^2\,dt + \alpha (|\omega|-c)^2, \quad \alpha >0,\label{eq:penal}
\end{align}

where $y$ is the solution of the state equation

∂_t y(x,t)   =σΔy(x,t) +χ_ω(x) u(t)              (x,t)∈Ω×(0,T],

y(x,t)=0               (x,t)∈∂Ω×(0,T],

y(0,x)=              x∈Ω ,

and $\Omega$ is a polygonal domain. The cost $J$ in  (<ref>) includes the additional term $\alpha (|\omega|-c)^2$ which accounts for the volume constraint $|\omega|=c$ in a penalty fashion. This slightly modifies the topological derivative formula, as it will be shown later. We derive a discretised version of the dynamics (<ref>)-(<ref>) via the method of lines. For this, we introduce a family of finite-dimensional approximating subspaces $V_h\subset {H^1_0(\Omega)}$, where h stands for a discretisaton parameter typically corresponding to gridsize  in finite elements/differences, but which can also be related to a spectral approximation of the dynamics. For each $f_h\in V_h$, we consider a finite-dimensional nodal/modal expansion of the form

\begin{equation}
f_h = \sum_{j=1}^{N} f_{j}\phi_{j}\,, \quad f_j\in\R\,,\phi_j\in V_h\,,
\end{equation}

where $\{\phi_i\}_{i=1}^{N}$ is a basis of $V_h$.  We denote the vector of coefficients associated to the expansion by $\underline f_h:=(f_1,\ldots,f_N)^\top$. In the method of lines, we approximate  the solution $y$ of (<ref>)-(<ref>) by a function $y_h$ in $C^1([0,T];V_h(\DD))$ of the type
$$y_h(x,t) = \sum_{j=1}^{N} y_{j}(t)\phi_{j}(x)\,,$$

for which we follow a standard Galerkin ansatz. Inserting  $y_h$ in the weak formulation (<ref>) and testing with $\varphi = \phi_{k}$, $k=1,\ldots, N$ leads to the following system of ordinary equations,

ẏ_h(t)=A_hy_h(t)+B_hu_h(t)  t∈(0,T],   y_h(0) = f_h,
where $M_h,K_h\in\R^{N\times N}$ and $B_h, \underline{f}_h \in \R^{N}$ are given by

A_h=-M^-1_hS_h ,  B_h = M^-1B̂_h  ,  f_h := M^-1_hf̂_h ,

(M_h)_ij=(ϕ_i,ϕ_j)_ L_2 ,  (S_h)_ij=σ(∇ϕ_i,∇ϕ_j)_L_2 ,

(B̂_h)_j=( χ_ω,ϕ_j)_ L_2 ,
(f̂_h)_j :=(f,ϕ_j)_L_2,   i,j=1,…,N .

Note that $\underline{y}_h = \underline{y}_h^{u_h,\underline{\y}_h,\omega} $ depends on $\y_h$, $u_h$, and $\omega$.  Given a discrete initial condition $\y_h\in V_h(\DD)$, the discrete costs are defined by

𝒥_1,h(ω,_h) := u_h∈infJ_h(ω,u,_h)=u_h∈inf ∫_0^T (y_h)^⊤M_hy_h+ γ| u_h(t)|^2 dt+ α(|ω|-c)^2,
𝒥_2,h(ω) = sup_f_h∈V_h
f_h_H^1≤1 𝒥_1,h(ω,f_h).
The  solution of the linear-quadratic optimal control problem in (<ref>) is given by

\[\bar{u}^{\omega,\y_h}(t)=-\gamma^{-1}B_h^\top \Pi_h(t)\underline{y}_h\,,\]

where $\Pi_h\in R^{N\times N}$ satisfies the differential matrix Riccati equation

\[-\frac{d}{dt}\Pi_h=A_h\Pi_h + \Pi_hA_h - \Pi_hB_h\gamma^{-1} B_h^\top\Pi_h + M_h \quad\text{ in } [0,T), \quad\Pi_h(T)=0\,.
\]

The coefficient vector of the discrete adjoint state $\bar p_h^{\y_h,\omega}(t)$ at time $t$ can be recovered directly by
$\underline{\bar{p}}_h^{\y_h,\omega}(t)=2\Pi_h(t)\underline{y}_h(t)$. Let us define the discrete analog of (<ref>),

𝔛_2,h(ω) := {f̅_h∈V_h:  sup_f_h∈V_h
f_h_H^1≤1 𝒥_1,h(ω,f_h) = 𝒥_1,h(ω,f̅_h)}.
Since we have the relation
𝒥_1,h(ω,_h) = (Π_h(0)_h, _h)_L_2 + α(|ω|-c)^2,
the maximisers $\y_h\in \mathfrak X_{2,h}(\omega)$ can be computed by solving the generalised Eigenvalue problem: find $(\lambda_h,\y_h)\in \R\times V_h$ such that

(Π_h(0) - λ_h S_h)_h =  0.
The biggest $\lambda_h = \lambda_h^{max}$ is then precisely the value $\mathcal J_{2,h}(\omega)$  and the normalised Eigenvectors for this Eigenvalue are the elements in $\mathfrak X_{2,h}(\omega)$:
𝔛_2,h(ω) = { _h:  _h ∈ker((Π_h(0) - λ_h^max K_h))   and     _h=1  }.

It is readily checked that if $\y_h \in \mathfrak X_{2,h}(\omega)$, then also $-\y_h\in \mathfrak X_{2,h}(\omega)$.  So if the Eigenspace for the largest eigenvalue is one-dimensional we have $\mathfrak X_{2,h}(\omega) = \{\y_h, -\y_h\}$. However, we know according to Corollary <ref> (now in a  discrete setting) that
𝒯𝒥_1,h(ω,_h)(η_0) = 𝒯𝒥_1,h(ω,-_h)(η_0)
for all $\eta_0\in \DD\setminus \partial \omega$ and $\y_h\in V_h$. Hence we can evaluate the topological derivative $\mathcal T\mathcal J_{2,h}(\omega)$ by picking either
$\y_h$ or $-\y_h$. A similar argumentation holds for the shape derivative.

§.§ Optimal actuator positioning: Shape derivative

Here we precise the gradient algorithm based upon a numerical realisation of the shape derivative.
We consider (<ref>)-(<ref>) with its discretisation (<ref>). Given a simply connected
actuator $\omega_0\subset \DD$ we employ the shape derivative
of $\mathcal J_1$ to find the optimal position. Let $\y_h\in V_h$. According to Corollary <ref> the derivative of $\mathcal J_{1,h}$ in the case $\mathcal U = \R$ is given by

D𝒥_1,h(ω,_h)(X) = -∫_∂ω u̅_h^_h,ω(t) ∫_0^T p̅_h^_h,ω(s,t) (X(s)·ν(s))  ds dt
for $X\in \ac C^1(\overbar \DD, \R^d)$. We assume that $\omega\Subset \DD$. We define the vector $b\in \R^d$ with the components
b_i := ∫_∂ω u̅_h^_h,ω(t) ∫_0^T p̅_h^_h,ω(s,t) (e_i·ν(s))  ds dt,
where $e_i$ denotes the canonical basis of $\R^d$. From this we can
construct an admissible descent direction by choosing any $\tilde X\in \ac C^1(\overbar \DD, \R^d)$ with $\tilde X|_{\partial \omega} = b$. Then it is obvious that $D\mathcal J_{1,h}(\omega,\y_h)(\tilde X)\le 0$. Let us use the the notation $b=-\nabla \mathcal J_{1,h}(\omega,\y_h)$.
We write $ (\id+t\nabla \mathcal J_{1,h}(\omega,\y_h))(\omega)$ to denote the moved actuator $\omega$ via the vector $b$. Note that only the position, but not the shape of $\omega$ changes by this operation. We refer to this procedure as Algorithm <ref> below.

Input: $\omega_0\in \mathfrak Y(\DD)$, $\y_h\in V_h$, $b_0 := -\nabla \mathcal J_{1,h}(\omega_0,\y_h)$, $n=0$, $\beta_0>0$, and $\eps >0$.
$|b_n| \ge \epsilon $
$\mathcal J_{1,h}((\id+\beta_n b_n)(\omega_n),\y_h) < \mathcal J_{1,h}(\omega_n,\y_h)$
$\beta_{n+1} \gets \beta_n$
$\omega_{n+1} \gets (\id+\beta_n b_n)(\omega_n)$
$b_{n+1} \gets -\nabla \mathcal J_{1,h}(\omega_{n+1},\y_h) $
$n \gets n+1$
decrease $\beta_n$
optimal actuator positioning $\omega_{opt}$
Shape derivative-based gradient algorithm for actuator positioning

§.§ Optimal actuator design: Topological derivative

As for the shape derivative, we now introduce a numerical approximation of the topological derivative formula which is embedded into a level-set method to generate an algorithm for optimal actuator design, i.e. including both shaping and position.
According to Theorem <ref> the discrete topological derivative of $\mathcal J_{1,h}$ is given by

𝒯𝒥_1,h(ω,_h)(η_0) = .

∫_0^T u^_h,ω_h(t)p̅_h^_h,ω(η_0,t) dt- 2α(|ω|-c)      if  η_0 ∈ω,

-∫_0^Tu^_h,ω_h(t) p̅_h^_h,ω(η_0,t) dt + 2α(|ω|-c)       if   η_0 ∈∖ω,

The level-set method is well-established in the context of shape optimisation and shape derivatives [2]. Here we use a level-set method for topological sensitivities as proposed in [4]. We recall that compared to the the formulation based on shape sensitivities, the topological approach has the advantage that multi-component actuators can be obtained via splitting and merging.

For a given actuator $\omega\subset \Omega$, we begin by defining the function

\[
g_h^{\y_h,\omega}(\zeta) =
-\int_0^T\overbar u^{\y_h,\omega}_h(t)\overbar p^{\y_h,\omega}_h(\zeta,t)\;dt + 2\alpha(|\omega|-c), \quad \zeta\in \overbar \DD\,
\]

which is continuous since the adjoint is continuous in space. Note that $\overbar p^{\y_h,\omega}$ and $\overbar u^{\y_h,\omega}$ depend on the actuator $\omega$. For other types of state equations where the shape variable enters into the differential operator (e.g. transmission problems [3]) this may not be the case and thus it particular of our setting.
The necessary optimality condition for the cost function $\mathcal J_{1,h}(\omega,\y_h)$ using the topological derivative are formulated as

g_h^_h,ω(x)     ≤0    for all  x∈ω,

g_h^_h,ω(x)     ≥0    for all  x∈Ω∖ω.

Since $g_h^{\y_h,\omega}$ is continuous this means that $	g_h^{\y_h,\omega}$ vanishes on $\partial \omega$ and hence

∫_0^Tu^_h,ω_h(t)p^_h,ω_h(ζ,t) dt = 2α(|ω|-c) ,    for all  ζ∈∂ω.
An (actuator) shape $\omega$ that satisfies (<ref>) is referred to as
stationary (actuator) shape. It follows from (<ref>) and  (<ref>), that $g_h^{\y_h,\omega}$ vanishes on the actuator boundary $\partial \omega$ of a stationary shape $\omega$.

We now describe the actuator $\omega$ via  an arbitrary level-set function $\psi_h \in V_h$, such that $\omega=\{x\in\Omega: \psi_h(x)<0\}$ is achieved via an update of an initial guess $\psi_h^0$
ψ^n+1_h = (1-β_n)ψ^n_h + β_n g_h^_h,ω_n/g_h^_h,ω_n,  ω_n := { x∈Ω: ψ^n_h(x)<0},
where $\beta_n$ is the step size of the method. The idea behind this update scheme is the following: if $\psi^n_h(x)<0$ and $g_h^{\y_h,\omega_n}(x)>0$, then we add a positive value to the level-set function, which means that we aim at removing actuator material. Similarly, if $\psi^n_h(x)>0$ and $g_h^{\y_h,\omega_n}(x)<0$, then we create actuator material. In all the other cases the sign of the level-sets remains unchanged. We present our version of the level-set algorithm in [4], which we refer to as Algorithm <ref>.

Input: $\psi^0_h\in V_h(\DD)$, $\omega_0:=\{x\in \overbar{\DD},\psi^0_h(x)<0\}$, $\beta_0>0$, $\y_h\in V_h$, and $\eps >0$.

$\|\omega_{n+1}-\omega_{n}\|\ge  \eps$
$\mathcal J_{1,h}(\{ \psi_h^{n+1}<0\},\y_h) < \mathcal J_{1,h}(\{\psi^n_h<0\},\y_h)$
$\psi^{n+1}_h \gets (1-\beta_n) \psi^n_h +  \beta_n \frac{g_h^{\y_h,\omega_n}}{\|g_h^{\y_h,\omega_n}\|}$
$\beta_{n+1} \gets \beta_n$
$\omega_{n+1} \gets \{\psi^{n+1}_h<0\}$
$n\gets n+1$
decrease $\beta_n$
optimal actuator $\omega_{opt}$
Level set algorithm for optimal actuator design

Algorithm <ref> is embedded inside a continuation approach over the quadratic penalty parameter $\alpha$ in (<ref>),  leading to actuators which approximate the size constraint in a sensible way, as opposed to a single solve with a large value of $\alpha$.

Finally, for the functional $\mathcal J_2(\omega)$ we may employ similar algorithms for shape and topological derivatives. We update the initial condition $f_h\in \mathfrak X_{2,h}(\omega) $ at each iteration whenever the actuator $\omega$ is modified.

§ NUMERICAL TESTS

We present a series of one and two-dimensional numerical tests exploring the different capabilities of the developed approach.

Test parameters and setup We establish some common settings for the experiments. For the 1D tests, we consider a piecewise linear finite element discretisation with 200 elements over $\Omega=(0,1)$, with $\gamma=10^{-3}$, $\sigma=0.01$, $c=0.2$, and $\epsilon=10^{-7}$. For the 2D tests, we resort to a Galerkin ansatz where the basis set is composed by the eigenfunctions of the Laplacian with Dirichlet boundary conditions over $\Omega=(0,1)^2$. We utilize the first 100 eigenfunctions. This idea has been previously considered in the context of optimal actuator positioning in [18], and its advantage resides in the lower computational burden associated to the Riccati solve. The actuator size constraint is set to $c=0.04$.
An important implementation aspect relates to the numerical approximation of the linear-quadratic optimal control problem for a given actuator. For the sake of simplicity, we consider the infinite horizon version of the costs $\cJ_1$ and $\cJ_2$. In this way, the optimal control problems are solved via an Algebraic Riccati Equation approach. The additional calculations associated to $\cJ_2$ and the set $\mathfrak X_2(\omega)$ are reduced to a generalized eigenvalue problem involving the Riccati operator $\Pi_h$. The shape and topological derivative formulae involving the finite horizon integral of $u$ and $p$ are approximated with a sufficiently large time horizon, in this case $T=1000$.

Actuator size constraint While in the abstract setting the actuator size constraint determines the admissible set of configurations, its numerical realisation follows a penalty approach, i.e. $\cJ_1(\omega,f)$ is as in (<ref>),

where $\cJ_1^{LQ}(\omega,f)$ is the original linear-quadratic (LQ) performance measure, and $\cJ_1^{\alpha}(\omega)=\alpha(|\omega|-c)^2$ is a quadratic penalization from the reference size. The cost $\cJ_2$ is treated analogously. In order to enforce the size constraint as much as possible and to avoid suboptimal configurations, the quadratic penalty is embedded within a homotopy/continuation loop. For a low initial value of $\alpha$, we perform a full solve of Algorithm <ref>, which is then used to initialized a subsequent solve with an increased value of $\alpha$. As it will be discussed in the numerical tests, for sufficiently large values of $\alpha$ and under a gradual increase of the penalty, results are accurate within the discretisation order.

Algorithm <ref> and level-set method The main aspect of Algorithm <ref> is the level-set update of the function $\psi^{n+1}_h$ which dictates the new actuator shape. In order to avoid the algorithm to stop around suboptimal solutions, we proceed to reinitialize the level-set function every 50 iterations. This is a well-documented practice for the level-set method, and in particular in the context of shape/topology optimisation [2, 4]. Our reinitialization consists of reinitialising  $\psi^{n+1}_h$  to be the signed distance function of the current actuator. The signed distance function is efficiently computed via the associated Eikonal equation, for which we implement the accelerated semi-Lagrangian method proposed in [1], with an overall CPU time which is negligible with respect to the rest of the algorithm.

Practical aspects All the numerical tests have been performed on an Intel Core i7-7500U with 8GB RAM, and implemented in MATLAB. The solution of the LQ control problem is obtained via the  command, the optimal trajectories are integrated with a fourth-order Runge-Kutta method in time. While a single LQ solve does not take more than a few seconds in the 2D case, the level-set method embedded in a continuation loop can scale up to approximately 30 mins. for a full 2D optimal shape solve.

§.§ Optimal actuator positioning through shape derivatives
In the first two tests we study the optimal positioning  problem (<ref>) of a single-component actuator of fixed width $0.2$ via the gradient-based approach presented in Algorithm <ref>. Tests are carried out for a given initial condition $y_0(x)$, i.e. the $\cJ_1$ setting.

Test 1 We start by considering $y_0(x)=\sin(\pi x)$, so the test is fully symmetric, and we expect the optimal position to be centered in the middle of the domain, i.e. at $x=0.5$. Results are illustrated in Figure <ref>, where it can be observed that as the actuator moves from its initial position towards the center, the cost $\cJ_1$ decays until reaching a stationary value. Results are consistent with the result obtained by inspection (Figure <ref> left), where the location of the center of the actuator has been moved throughout the entire domain.

Test 1. Left: different single-component actuators with  different centers have been spanned over the domain, locating the minimum value of $\cJ_1$ for the center at $x=0.5$. Center: starting from an initial guess for the actuator far from $0.5$, the gradient-based approach of Algorithm <ref> locates the optimal position in the middle. Right: as the actuator moves towards the center in the subsequent iterations of Algorithm <ref>, the value $\cJ_1$ decays until reaching a stationary point.

Test 2 We consider the same setting as in the previous test, but we change the initial condition of the dynamics to be $y_0(x)=100|x-0.7|^4+x(x-1)$, so the setting is asymmetric and the optimal position is different from the center. Results are shown in Figure <ref>, where the numerical solution coincides with the result obtained by inspecting all the possible locations.

Test 2. Left: inspecting different values of $\cJ_1$ by spanning actuators with different centers, the optimal center location is found to be close to $0.2$ . Center: the gradient-based approach steers the initial actuator to the optimal position. Right:  the value $\cJ_1$ decays until reaching a stationary point, which coincides with the minimum for the first plot on the left.

§.§ Optimal actuator design through topological derivatives
In the following series of experiments we focus on 1D optimal actuator design, i.e. problems (<ref>) and (<ref>) without any further parametrisation of the actuator, thus allowing multi-component structures. For this, we consider the approach combining the topological derivative, with a level-set method, as summarized in Algorithm <ref>.

Test 3 For $y_0(x)=\max(\sin(3\pi x),0)^2$, results are presented in Figures <ref> and <ref> . As it can be expected from the symmetry of the problem, and from the initial condition, the actuator splits into two equally sized components. We carried out two types of tests, one without and one with a continuation strategy with respect to $\alpha$. Without a continuation strategy, choosing $\alpha=10^3$ we obtain the result depicted in Figure <ref> (b). With a continuation strategy, as the  penalty increases, the size of the components decreases until approaching the total size constraint.  The behavior of this continuation approach is shown in Table <ref>. When $\alpha$ is increased, the size of the actuator tends to $0.2$, the reference size, while the LQ part of $\cJ_1$, tends to a stationary value. For a final value of $\alpha=10^4$, the overall cost $\cJ_1$ obtained via the continuation approach is approx. 80 times smaller than the value obtained without any initialisation procedure, see Figure <ref> (b)-(d).
Figure <ref> illustrates some basic relevant aspects of the level-set approach, such as the update of the shape (left), the computation of the level-set update upon $\beta_n$ and $\psi_h^{n}$ (middle), and the decay of the value $J_1$ (right).

[$\alpha=10^3$, no init.]
Test 3. (a) Initial condition $y_0(x)=\max(\sin(3\pi x),0)^2$. (b) Optimal actuator for $\alpha=10^3$, without initialization via increasing penalization. (c) Optimal actuator for $\alpha=10^{-1}$, subsequently used in the quadratic penalty approach. (d) Optimal actuator for $\alpha=10^3$, via increasing penalization.

$\alpha$    $\cJ_1$      $\cJ_1^{LQ}$     $\cJ_1^{\alpha}(size)$   iterations

0.1      1.84$\times10^{-2}$     1.62$\times10^{-2}$      2.30$\times10^{-3}$ (0.35)    225

1         2.35$\times10^{-2}$     2.26$\times10^{-2}$     9.10$\times10^{-4}$ (0.23)     226

10       2.56$\times10^{-2}$     2.46$\times10^{-2}$      1.00$\times10^{-3}$ (0.21)    316

$10^{2}$     3.46$\times10^{-2}$     2.46$\times10^{-2}$      1.00$\times10^{-2}$ (0.21)    226

$10^{3}$     0.12           2.46$\times10^{-2}$      1.00$\times10^{-1}$ (0.21)    226

$10^{3}$*     8.18          8.00$\times10^{-2}$     8.10 (0.29)    629

Test 3. optimisation values for $y_0(x)=\max(\sin(3\pi x),0)^2$. Each row is initialized with the optimal actuator corresponding to the previous one, except for the last row with $\alpha=10^3*$, illustrating that incorrectly initialized solves lead to suboptimal solutions. The reference size for the actuator is $0.2$ .

Test 3. Level set method implemented in Algorithm <ref>. Left: starting from an initial actuator, the topological derivative of the cost is computed and an updated actuator is obtained. The new shape is evaluated according to its closed-loop performance. If the update is rejected, the parameter $\beta_n$ is reduced. Middle: the level-set approach generates an update of the actuator shape based on the information from $\psi_h^{n}$, $\beta_n$ and $g_{\omega_n}$. Right: This iterative loop generates a decay in the total cost $J_1$, (which accounts for both the closed-loop performance of the actuator and its volume constraint).

Test 4 We repeat the setting of Test 3 with a nonsymmetric initial condition $y_0(x)=\sin(3\pi x)^2\chi_{\{x<2/3\}}(x)$. Results are presented in Table <ref> and Figure <ref>, which illustrate the effectivity of the continuation approach, which generates an optimal actuator with two components of different size, see Figure <ref> and compare with Figure <ref>.

$\alpha$    $\cJ_1$      $\cJ_1^{LQ}$     $\cJ_1^{\alpha}(size)$   iterations

0.1      6.48$\times10^{-2}$     6.31$\times10^{-2}$      1.7$\times10^{-3}$ (0.33)    229

1         8.0$\times10^{-2}$     6.31$\times10^{-2}$     1.69-2 (0.33)     226

10       0.176     0.164      1.23$\times10^{-2}$ (0.235)    226

$10^{2}$     0.207     0.184      2.25$\times10^{-2}$ (0.215)    316

$10^{3}$     0.234          0.209      2.50$\times10^{-2}$ (0.195)    316

$10^{4}$     0.459           0.209      0.250 (0.195)    316

$10^{4}$*     9.09           9.66$\times10^{-2}$      9 (0.23)    629

Test 4. optimisation values for $y_0(x)=\sin(3\pi x)^2\chi_{x<2/3}(x)$. Each row is initialized with the optimal actuator corresponding to the previous one, except for the last row with $\alpha=10^4*$, illustrating that incorrectly initialized solves lead to suboptimal solutions. The reference size for the actuator is $0.2$ .

[$\alpha=10^4$, no init.]
Test 4. (a) Initial condition $y_0(x)=\sin(3\pi x)^2\chi_{\{x<2/3\}}(x)$. (b) Optimal actuator for $\alpha=10^4$, without initialization via increasing penalization. (c) Optimal actuator for $\alpha=10^{-1}$, subsequently used in the quadratic penalty approach. (d) Optimal actuator for $\alpha=10^4$, via increasing penalization.

Test 5 We now turn our attention to the optimal actuator design for the worst-case scenario among all the initial conditions, i.e. the $\cJ_2$ setting. Results are presented in Figure <ref> and Table <ref>. The worst-case scenario corresponds to the first eigenmode of the Riccati operator (Figure <ref>), which generates a two-component symmetric actuator (Figure <ref>). This is only observed within the continuation approach. For a large value of $\alpha$ without initialisation, we obtain a suboptimal solution with a single component (last row of Table <ref>, Figure <ref>).

[$\mathfrak X_2(\omega)$]
[$\alpha=10^3$, no init.]
Test 5. (a) First eigenmode of the Riccati operator, which corresponds to the set $\mathfrak X_2(\omega)$. (b) Optimal actuator for $\alpha=10^3$, without initialization via increasing penalization. (c) Optimal actuator for $\alpha=10^{-1}$, subsequently used in the quadratic penalty approach. (d) Optimal actuator for $\alpha=10^3$, via increasing penalization.

$\alpha$    $\cJ_2$      $\cJ_2^{LQ}$     $\cJ_2^{\alpha}(size)$   iterations

0.1      0.402     0.401      1.1$\times10^{-3}$ (0.305)    307

1         0.369     0.364     4.0$\times10^{-4}$ (0.22)     225

10       0.343     0.342      1.0$\times10^{-3}$ (0.19)    228

$10^{2}$     0.352     0.342      1.0$\times10^{-2}$ (0.19)    226

$10^{3}$     0.442          0.342      0.1 (0.19)    226

$10^{3}$*     0.761           0.536      0.225 (0.215)    941

Test 5. optimisation values for $\cJ_2$. Each row is initialized with the optimal actuator corresponding to the previous one, except for the last row with $\alpha=10^{3}$*. The reference size for the actuator is $0.2$ .

Test 6 As an extension of the capabilities of the proposed approach, we explore the $\cJ_2$ setting with space-dependent diffusion. For this test, the diffusion operator $\sigma\Delta y$ is rewritten as $div(\sigma(x)\nabla y)$, with $\sigma(x)=(1-\max(\sin(9\pi x),0))\chi_{\{x<0.5\}}(x)+10^{-3} $. Iterates of the continuation approach are presented in Table <ref>. Again, the lack of a proper initialization of Algortithm <ref> with a large value of $\alpha$ leads to a poor satisfaction of both the size constraint and the LQ performance, which is solved via the increasing penalty approach. A two-component actuator present in the area of smaller diffusion is observed in Figure <ref>.

$\alpha$    $\cJ_2$      $\cJ_2^{LQ}$     $\cJ_2^{\alpha}(size)$   iterations

0.1      1.792     1.743      4.97$\times10^{-2}$ (0.908)    194

1         2.240     1.743     0.497 (0.908)     228

10       4.734     4.462      0.272 (0.365)    225

$10^{2}$     3.134     3.071      6.25$\times10^{-2}$ (0.175)    538

$10^{3}$     1.023      0.998      0.025 (0.195)    226

$10^{4}$     1.248      0.998      0.250 (0.195)    226

$10^{4}$*     28.19       3.195      25.0 (0.25)    673

Test 6. $\cJ_2$ values with space-dependent diffusion $\sigma(x)=(1-\max(\sin(9\pi x),0))\chi_{\{x<0.5\}}(x)+10^{-3} $. Each row is initialized with the optimal actuator corresponding to the previous one, except for the last row with $\alpha=10^{4}$*. The reference size for the actuator is $0.2$ .

[$\mathfrak X_2(\omega)$]
Test 6. (a) First eigenmode of the Riccati operator, which corresponds to the set $\mathfrak X_2(\omega)$. (b) space-dependent diffusion coefficient $\sigma(x)=(1-\max(\sin(9\pi x),0))\chi_{\{x<0.5\}}(x)+10^{-3} $. (c) Optimal actuator for $\alpha=10$, subsequently used in the quadratic penalty approach. (d) Optimal actuator for $\alpha=10^4$, via increasing penalization.

§.§ Two-dimensional optimal actuator design
We now turn our attention into assessing the performance of Algorithm <ref> for two-dimensional actuator topology optimisation. While this problem is computationally demanding, the increase of degrees of freedom can be efficiently handled via modal expansions, as explained at the beginning of this Section. We explore both the $\cJ_1$ and $\cJ_2$ settings.

Test 7 This experiment is a direct extension of Test 3. We consider a unilaterally symmetric initial condition $y_0(x_1,x_2)=\max(\sin(4\pi(x_1-1/8)),0)^3\sin(\pi x_2)^3$, inducing a two-component actuator. The desired actuator size is $c=0.04$. The evolution of the actuator design for increasing values of the penalty parameter $\alpha$ is depicted in Figure <ref>. We also study the closed-loop performance of the optimal shape.  For this purpose
the running cost associated to the optimal actuator is compared against an ad-hoc design, which consists of a cylindrical actuator of desired size placed in the center of the domain, see Figure <ref> . The closed-loop dynamics of the optimal actuator generate a stronger exponential decay compared to the uncontrolled dynamics and the ad-hoc shape.

Test 7. (a) initial condition $y_0(x_1,x_2)=\max(\sin(4\pi(x_1-1/8)),0)^3\sin(\pi x_2)^3$ for $\cJ_1$ optimisation. (b) within the level-set method, the actuator is updated according to the zero level-set of the function $\psi_h^{n+1}$. (c) to (f) optimal actuators for different volume penalties.

Test 7. Closed-loop performance for different shapes. The running cost in $\cJ_1$ is evaluated for uncontrolled dynamics ($u\equiv 0$), an ad-ho cylindrical actuator located in the center of the domain, and the optimal shape (Figure <ref>). Closed-loop dynamics of the optimal shape decay faster.

Test 8 In an analogous way as in Test 5, we study the optimal design problem associated to $\cJ_2$. The first eigenmode of the Riccati operator is shown in Figure <ref>. The increasing penalty approach (Figs. <ref> to <ref>) shows a complex structure, with a hollow cylinder and four external components. The performance of the closed-loop optimal solution is analysed in Figure <ref>, with a considerably faster decay compared to the uncontrolled solution, and to the ad-hoc design utilised in the previous test.

[$\mathfrak X_2(\omega)$]

Test 8. (a) first eigenmode of the Riccati operator. (b) within the level-set method, the actuator is updated according to the zero level-set of the function $\psi_h^{n+1}$. (c) to (f) optimal actuators for different volume penalties.

Test 8. Closed-loop performance for different shapes. The running cost in $\cJ_2$ is evaluated for uncontrolled dynamics ($u\equiv 0$), a suboptimal cylindrical actuator of size $c$ located in the center of the domain, and the optimal shape  with five components (Figure <ref>). Closed-loop dynamics of the optimal shape decay faster.

§ CONCLUDING REMARKS
In this work we have developed an analytical and computational framework for optimisation-based actuator design. We derived shape and topological sensitivities formulas which account for the closed-loop performance of a linear-quadratic controller associated to the actuator configuration. We embedded the sensitivities into gradient-based and level-set methods to numerically realise the optimal actuators. Our findings seem to indicate that from a practical point of view, shape sensitivities are a good alternative whenever a certain parametrisation of the actuator is fixed in advance and only optimal position is sought. Topological sensitivities are instead suitable for optimal actuator design in a wider sense, allowing the emergence of nontrivial multi-component structures, which would be difficult to guess or parametrise a priori. This is a relevant fact, as most of the engineering literature associated to computational optimal actuator positioning is based on heuristic methods which strongly rely on experts' knowledge and tuning. Extensions concerning robust control design and semilinear parabolic equation are in our research roadmap.

§ APPENDIX

§.§ Differentiability of maximum functions

In order to prove Lemma <ref> we recall the following
Danskin-type lemmas.

Let $\mathfrak V_1$ be a nonempty set and let $\mathcal G:[0,\tau]\times \mathfrak V_1\rightarrow \R$ be a function, $\tau >0$. Introduce the function $g_1:[0,\tau]\rightarrow \R$,

g_1(t) := sup_x∈𝔙_1 𝒢(t,x),
and  let $\ell:[0,\tau]\rightarrow \R$ be any function such that $\ell(t)>0$ for $t\in (0,\tau]$ and $\ell(0)=0$.
We give sufficient conditions that guarantee that the limit
d/dℓ g_1(0^+) := lim_t↘0 g_1(t)-g_1(0)/ℓ(t)
exists. For this purpose we introduce for each $t$ the set of maximisers

𝔙_1(t) = {x^t∈𝔙_1:  sup_x∈𝔙_1 𝒢(t,x) = 𝒢(t,x^t) }.
The next lemma can be found with slight modifications in <cit.>.

Let the following hypotheses be satisfied.

$(i)$  For all $t$ in $[0,\tau]$  the set $\mathfrak V_1(t)$  is nonempty,

$(ii)$ the limit
∂_ℓ𝒢(0^+,x) := lim_t↘0 𝒢(t,x)-𝒢(0,x) /ℓ(t)
exists for all $x\in \mathfrak V_1(0)$.

(A2) For all real null-sequences $(t_n)$ in $(0,\tau]$ and all sequence $(x_{t_n})$ in $\mathfrak V_1(t_n)$, there exists a subsequence $(t_{n_k})$ of $(t_n)$, $(x_{t_{n_k}})$ in $\mathfrak V_1(t_{n_k})$  and $x_0$ in $\mathfrak V_1(0)$, such that
lim_k→∞ 𝒢( t_n_k ,x_t_n_k) -
𝒢(0,x_t_n_k)/ℓ(t_n_k) = ∂_ℓ𝒢(0^+,x_0).

Then $g_1$  is differentiable at $t=0^+$ with derivative
d/dtg_1(t)|_t=0^+ = max_x∈𝔙_1(0) ∂_ℓ𝒢(0^+,x).

§.§ Proof of Lemma <ref>

Our strategy is to prove Lemma <ref> by applying

Lemma <ref> to the function $\mathcal G(t,y):= \inf_{x\in \mathfrak V}  G(t,x,y)$ with $\mathfrak V_1 := \mathfrak V$. This will show that $g(t) := \sup_{y\in \mathfrak V} \mathcal G(t,y)$ is right-differentiable at $t=0^+$. By construction Assumption (A0) of Lemma <ref> is satisfied.

Step 1: For every $t\in [0,\tau]$ and $y\in \mathfrak V$ we have $\mathcal G(t,y) = G(t,x^{t,y},y)$. Hence

𝒢(t,y)-𝒢(0,y) =     G(t,x^t,y,y) - G(0,x^0,y,y)

= G(t,x^t,y,y) - G(0,x^t,y,y) + G(0,x^t,y,y)- G(0,x^0,y,y)_≥0

≥G(t,x^t,y,y) - G(0,x^t,y,y)

and similarly

𝒢(t,y)-𝒢(0,y) =     G(t,x^t,y,y) - G(0,x^0,y,y)

= G(t,x^t,y,y) - G(t,x^0,y,y)_≤0 + G(t,x^0,y,y)- G(0,x^0,y,y)

≤G(t,x^0,y,y)- G(0,x^0,y,y).

Therefore using Assumption (A2) of Lemma <ref> we obtain from (<ref>) and (<ref>)
lim inf_t↘0 𝒢(t,y)-𝒢(0,y)/ℓ(t) ≥∂_ℓG(0^+,x^0,y,y) ≥lim sup_t↘0 𝒢(t,y)-𝒢(0,y)/ℓ(t).
Hence Assumption (A1) of Lemma <ref> is satisfied.

Step 2: For every $t\in [0,\tau]$ and $y^t\in \mathfrak V(t)$ we have
$\mathcal G(t,y^t) = G(t,x^{t,y^t},y^t)$ and hence

𝒢(t,y^t)-𝒢(0,y^t)  =      G(t,x^t,y^t,y^t) - G(0,x^0,y^t,y^t)

=     G(t,x^t,y^t,y^t) - G(0,x^t,y^t,y^t) + G(0,x^t,y^t,y^t) - G(0,x^0,y^t,y^t)_≥0

≥    G(t,x^t,y^t,y^t) - G(0,x^t,y^t,y^t)

and similarly

𝒢(t,y^t)-𝒢(0,y^t)  =     G(t,x^t,y^t,y^t) - G(t,x^0,y^t,y^t)_≤0 + G(t,x^0,y^t,y^t) - G(0,x^0,y^t,y^t)

≤    G(t,x^0,y^t,y^t) - G(0,x^0,y^t,y^t).

Thanks to Assumption (A3) of Lemma <ref> For all real null-sequences $(t_n)$ in $(0,\tau]$ and all sequences $(y^{t_n})$,  $y^{t_n}\in \mathfrak V(t_n)$, there exists a subsequence $(t_{n_k})$ of $(t_n)$, $(y^{t_{n_k}})$ of $(y^{t_n})$, and $y^0$ in $\mathfrak V(0)$, such that

lim_k→∞ G( t_n_k ,x^t_n_k,y^t_n_k, y^t_n_k ) -
G(0,x^t_n_k,y^t_n_k, y^t_n_k)/ℓ(t_n_k) = ∂_ℓG(0^+,x^0,y^0,y^0)

lim_k→∞ G( t_n_k ,x^0,y^t_n_k, y^t_n_k ) -
G(0,x^0,y^t_n_k, y^t_n_k)/ℓ(t_n_k) = ∂_ℓG(0^+,x^0,y^0,y^0).
Hence choosing $t=t_{n_k}$ in (<ref>) we obtain

lim inf_k→∞ 𝒢(t_n_k,y^t_n_k)-𝒢(0,y^t_n_k) /ℓ(t_n_k)

(<ref>)≥   lim inf_k→∞ G( t_n_k ,x^t_n_k,y^t_n_k, y^t_n_k ) -
G(0,x^t_n_k,y^t_n_k, y^t_n_k)/ℓ(t_n_k)

and similarly choosing $t=t_{n_k}$ in (<ref>) we get

lim sup_k→∞ 𝒢(t_n_k,y^t_n_k)-𝒢(0,y^t_n_k) /ℓ(t_n_k)

(<ref>)≤   lim sup_k→∞ G( t_n_k ,x^0,y^t_n_k, y^t_n_k ) -
G(0,x^0,y^t_n_k, y^t_n_k)/ℓ(t_n_k)

Combining (<ref>) and (<ref>) we conclude that
lim_k→∞ 𝒢(t_n_k,y^t_n_k)-𝒢(0,y^t_n_k) /ℓ(t_n_k)  = ∂_ℓG(0^+,x^0,y^0,y^0),
which is precisely Assumption (A2) of Lemma <ref>.

Step 1 and Step 2 together show that Assumptions (A1) and (A2) of Lemma <ref> are satisfied and this finishes the proof.

[1]
A. Alla, M. Falcone, and D. Kalise, An efficient policy iteration
algorithm for the solution of dynamic programming equations, SIAM J. Sci.
Comput., 35 (2015), pp. A181–A200.

[2]
G. Allaire, F. Jouve, and A.-M. Toader, Structural optimization
using sensitivity analysis and a level-set method, J. Comput. Phys., 194
(2004), pp. 363–393.

[3]
S. Amstutz, Augmented Lagrangian for cone constrained topology
optimization, Comput. Optim. Appl., 49 (2011), pp. 101–122,

[4]
S. Amstutz and H. Andrä, A new algorithm for topology optimization
using a level-set method, J. Comput. Phys., 216 (2006), pp. 573–588,

[5]
A. Bensoussan, Optimization of sensors' location in a distributed
filtering problem, Springer Berlin Heidelberg, Berlin, Heidelberg, 1972,
pp. 62–84, <https://doi.org/10.1007/BFb0064935>,

[6]
M. C. Delfour and J.-P. Zolésio, Shape sensitivity analysis via
min max differentiability, SIAM J. Control Optim., 26 (1988), pp. 834–862,

[7]
M. C. Delfour and J.-P. Zolésio, Shapes and geometries: Metrics,
analysis, differential calculus, and optimization, Society for Industrial
and Applied Mathematics (SIAM), Philadelphia, PA, second ed., 2011.

[8]
V. F. Demyanov and V. N. Malozëmov, Einführung in
Minimax-Problem, Akademische Verlagsgesellschaft Geest & Portig K.-G.,
German translation.

[9]
A. El Jaïand A. J. Pritchard, Sensors and controls in the
analysis of distributed systems, Ellis Horwood Series: Mathematics and its
Applications, Ellis Horwood Ltd., Chichester; Halsted Press [John Wiley &
Sons, Inc.], New York, 1988.
Translated from the French by Catrin Pritchard and Rhian Pritchard.

[10]
L. C. Evans, Partial differential equations, vol. 19 of Graduate
Studies in Mathematics, American Mathematical Society, Providence, RI, 1998.

[11]
F. Fahroo and M. A. Demetriou, Optimal actuator/sensor location for
active noise regulator and tracking control problems, J. Comput. Appl.
Math., 114 (2000), pp. 137–158,
Control of partial differential equations (Jacksonville, FL, 1998).

[12]
M. I. Frecker, Recent advances in optimization of smart structures
and actuators, Journal of Intelligent Material Systems and Structures, 14
(2003), pp. 207–216, <https://doi.org/10.1177/1045389X03031062>,

[13]
P. Hébrard and A. Henrot, A spillover phenomenon in the optimal
location of actuators, SIAM J. Control Optim., 44 (2005), pp. 349–366,

[14]
M. Hintermüller, C. N. Rautenberg, M. Mohammadi, and M. Kanitsar, Optimal sensor placement: a robust approach, WIAS preprint,  (2016), p. 34
pp., <http://www.wias-berlin.de/preprint/2287/wias_preprints_2287.pdf>.

[15]
K. Ito, K. Kunisch, and G. H. Peichl, Variational approach to shape
derivatives, ESAIM Control Optim. Calc. Var., 14 (2008), pp. 517–539,

[16]
D. Kasinathan and K. Morris, $\Bbb H_\infty$-optimal actuator
location, IEEE Trans. Automat. Control, 58 (2013), pp. 2522–2535,

[17]
Laurain, A. and Sturm, K., Distributed shape derivative via
averaged adjoint method and applications, ESAIM: M2AN, 50 (2016),
pp. 1241–1267.

[18]
K. Morris, Linear-quadratic optimal actuator location, IEEE Trans.
Automat. Control, 56 (2011), pp. 113–124,

[19]
K. Morris, M. A. Demetriou, and S. D. Yang, Using $\Bbb
H_2$-control performance metrics for the optimal actuator location of
distributed parameter systems, IEEE Trans. Automat. Control, 60 (2015),
pp. 450–462, <https://doi.org/10.1109/TAC.2014.2346676>,

[20]
A. A. Novotny and J. Sokołowski, Topological derivatives in shape
optimization, Interaction of Mechanics and Mathematics, Springer,
Heidelberg, 2013, <https://doi.org/10.1007/978-3-642-35245-4>,

[21]
Y. Privat, E. Trélat, and E. Zuazua, Actuator design for parabolic
distributed parameter systems with the moment method, SIAM J. Control
Optim., 55 (2017), pp. 1128–1152, <https://doi.org/10.1137/16M1058418>,

[22]
J. Sokołowski and A. Żochowski, On the topological derivative
in shape optimization, SIAM J. Control Optim., 37 (1999), pp. 1251–1272,

[23]
K. Sturm, Minimax Lagrangian approach to the differentiability of
nonlinear PDE constrained shape functions without saddle point assumption,
SIAM J. Control Optim., 53 (2015), pp. 2017–2039,

[24]
K. Sturm, Shape differentiability under non-linear PDE
constraints, in New trends in shape optimization, vol. 166 of Internat. Ser.
Numer. Math., Birkhäuser/Springer, Cham, 2015, pp. 271–300,

[25]
F. Tröltzsch, Optimale Steuerung partieller
Differentialgleichungen: Theorie, Verfahren und Anwendungen, Vieweg, 2005,

[26]
S. Valadkhan, K. Morris, and A. Khajepour, Stability and robust
position control of hysteretic systems, Internat. J. Robust Nonlinear
Control, 20 (2010), pp. 460–471, <https://doi.org/10.1002/rnc.1457>,

[27]
M. van de Wal and B. de Jager, A review of methods for input/output
selection, Automatica J. IFAC, 37 (2001), pp. 487–510,

[28]
J. Wloka, Partielle Differentialgleichungen, B. G. Teubner,
Stuttgart, 1982.
Sobolevräume und Randwertaufgaben. [Sobolev spaces and boundary
value problems], Mathematische Leitfäden. [Mathematical Textbooks].