# Input Regularization for Integer Optimal Control in BV with Applications to
Control of Poroelastic and Poroviscoelastic Systems

Lorena Bociu  Department of Mathematics, North Carolina State University,
Raleigh, NC, USA  Paul Manns  Faculty of Mathematics, TU Dortmund University,
Dortmund, Germany  Marvin Severitt  Faculty of Mathematics, TU Dortmund
University, Dortmund, Germany  Sarah Strikwerda  Department of Mathematics,
North Carolina State University, Raleigh, NC, USA

###### Abstract

We consider a class of integer control problems that are governed by fluid
flows through deformable porous media equations that arise in biomechanics. We
show that the regularity assumptions for the convergence analysis of the
recently proposed optimization algorithm in [33] hold for the poro-elastic
case and are violated for the poro-visco-elastic case of the governing
equations. To bypass the restrictive regularity assumptions in case of
violation, we introduce an additional partial regularization of the control
inputs by means of mollification and prove a $\Gamma$-convergence-type result
when the support parameter of the mollification is driven to zero. Numerical
computations show that a homotopy that drives the support parameter of the
mollification to zero may often lead to lower objective values and smaller
remaining instationarities, thereby of course inducing an additional
computational cost. Executing this homotopy is not always necessary, however,
and we also observe cases where the original algorithm is able to obtain
points (final iterates) of similar quality as the homotopy when applied
without the additional input regularization even if the regularity conditions
are violated.

##  1 Introduction

Let $\alpha>0$ and $T>0$. We are interested in solving the following
optimization problem

$\displaystyle\begin{aligned} \min_{w\in L^{2}(0,T)}\
&j(w)+\alpha\operatorname*{TV}(w)\\\ \text{s.t.}\quad\quad&w(t)\in
W\subset\mathbb{Z}\text{ for almost all (a.a.) }t\in(0,T),\end{aligned}$ (P)

where $W\subset\mathbb{Z}$ is a finite set of integers, $\operatorname*{TV}$
denotes the total variation of the $W$-valued control input function $w$,
$j:L^{2}(0,T)\to\mathbb{R}$ is the objective that takes the form
$j(w)\coloneqq J(G(w),w)$, where $J:X\to\mathbb{R}$ is a coercive and lower
semicontinuous function for a Banach space $X$ that is the state space of some
PDE. The function $G:L^{2}(0,T)\to X$ is the continuous solution operator of a
PDE, more specifically, the solution operator of a PDE system describing fluid
flow through deformable porous media. We note that (P) admits a solution in
this setting, see [32].

The optimization problem (P) falls in the class of so-called _integer optimal
control problems_. Driven by versatile applications from the optimization of
supply and traffic networks [38, 21, 22, 20] over automotive control [28, 19]
to topology optimization [46, 24], this problem class has attracted
considerable research interest in recent years. Different methods have been
proposed to treat integer optimal control problems. One of them is the
combinatorial integral approximation decomposition [43] that splits the
optimization into the solution of a relaxed problem, where $W$ replaced by a
_one-hot encoding_ and then the convex hull is analyzed and a fast algorithm
that computes a $W$-valued control from the relaxation [41, 42, 23, 36]. This
method requires the ability to produce highly oscillating control functions,
which are undesirable in many applications and can therefore not be applied to
(P) if $\alpha>0$.

The $\operatorname*{TV}$-term that influences (P) for $\alpha>0$ has been
prevalent in mathematical image analysis since the 1990s, in particular the
work [40]. We give the references [48, 12, 13, 17, 10, 30, 25] but note that
this reflects only a small portion of the research in this area.

Several authors have incorporated $\operatorname*{TV}$-terms in optimal
control problems [34, 27, 16], in tight relaxations of integer optimal control
problems from topology optimization [14], and in the approximation step of the
combinatorial integral approximation decomposition [4, 44, 5].

In [33] a trust-region method has been proposed and analyzed to directly solve
problems of the form (P). The trust-region subproblems are integer linear
programs after discretization, which can be solved efficiently with graph-
based [45] and dynamic programming-based approaches.

For a one-dimensional domains $(0,T)$ as is present in (P), the
$\operatorname*{TV}$-term of a $W$-valued function $w$ is the sum of the jump
heights of the function $w$, implying that only finitely many jumps occur
because the height of a single jump is bounded below (in case
$W\subset\mathbb{Z}$ by $1$ but always by some constant if $|W|<\infty$).

This constitutes a desirable regularization because the application underlying
(P) usually does not permit infinitely many jumps or general high-frequency
switching between different control modes. This is also the case for the class
of PDEs that we consider to constrain our problem (P).

While these optimization problems seems to be suited for the algorithmic
framework proposed in [33], the regularity assumptions for its convergence
analysis cannot always be satisfied. Therefore, we enforce the assumed
regularity by adding a regularization of the control input when passed to $j$.
Specifically, the function $j$ is altered to $j\circ K_{\varepsilon}$, where
$K_{\varepsilon}$ is a convolution operator arising from a standard mollifier
with parameter $\varepsilon$.

While we are not able to give a full answer to the relationship between (P)
and its regularized variant with $j$ replaced by $j\circ K_{\varepsilon}$, we
make the following contributions in the light of the algorithmic framework
proposed in [33].

##### Contributions

We verify the regularity assumptions required in [33] for the altered
objective. We prove the lower and bound inequalities ($\Gamma$-convergence) of
the altered optimization problems when driving $\varepsilon$ (the parameter
controlling the size of the support of the mollifier) to zero. The
$\Gamma$-convergence result is achieved with respect to weak-∗ and strict
convergence in the weak-∗ closed subset of functions of bounded variations
that are feasible for (P).

While this gives that global minimizers of the altered optimization problems
converge to global minimizers of (P), the lower and upper bound inequality do
not imply that this holds true for stationary points as are produced by the
algorithm proposed in [33]. We consider a homotopy that drives $\varepsilon\to
0$ and applies meaningful termination criteria for each run (tightening of
minimal trust-region radius to determine that no progress is made for
$\varepsilon\to 0$ and achievement of a certain predicted reduction). We show
that the cluster points are strict limits of their approximating sequences,
implying that the the homotopy does not _overlook_ cheap reductions of the
objective that can, for example, be obtained by removing small jumps from the
control.

While the regularity assumptions are often satisfied for cases that
$j(w)=J(G(w))$, where $G$ is the solution operator of a parabolic PDE, we show
that the regularity assumptions on (P) are generally not satisfied for a class
of problems that arise in the considered application. Specifically, we show
that they do hold for the considered class of equations when the viscosity
parameter is chosen to be zero for two different control inputs and that they
are violated when the viscosity parameter is chosen strictly larger than zero
all else equal.

We provide numerical results for two instances of the class of considered PDEs
that differ in their dynamics and analytical properties. The numerical results
show that while the homotopy can help to find better objective values and
points of lower instationarity, the practical performance of the algorithm
without the input regularization may be on par with the homotopy.
Consequently, the lack of regularity may not always impair the practical
performance and may therefore be outweighed by the consumption of much less
running time than the homotopy.

##### The remainder of the paper is structured as follows.

In section 2, we introduce our notation, briefly recall functions of bounded
variation, and define the closed subset that corresponds to the feasible set
of our optimization problem. In section 3, we describe the sequential linear
integer programming (SLIP) algorithm provided in [33]. We introduce and
analyze the effect of the control input regularization in section 4. The
considered class of fluid-solid mixture systems and the discussion of the
regularity assumptions of Algorithm 1 with respect to these coupled systems of
PDEs are given in section 5. We provide our computational setup, experiments,
and results in section 6. Finally, we draw our conclusions in section 7.

##  2 Notation and Primer on Functions of Bounded Variation

##### Notation.

Let $X$ be a Banach space. Then we denote its topological dual space by the
symbol $X^{*}$. For a given domain $\Omega$, $L^{2}(\Omega)$ is the usual
Hilbert space of square-integrable functions with inner product given by
$(\cdot,\cdot)$. When the domain $\Omega$ is not clear from context, the
$L^{2}$ inner product will be denoted as $(\cdot,\cdot)_{\Omega}$.
Furthermore, we use the standard notation $\displaystyle
H^{1}_{\Gamma_{*}}(\Omega)=\left\\{f\in H^{1}(\Omega)\,\middle|\,\tau
f\big{|}_{\Gamma_{*}}=0\right\\}$, for any $\Gamma_{*}\subset\partial\Omega$.
. Additionally, for any Hilbert space $W$, we define
$L^{2}(0,T;W)=\left\\{u:[0,T]\to W\middle|u\text{ is measurable and
}\int_{0}^{T}\|u(t)\|_{W}^{2}dt<\infty\right\\}$ with the inner-product
$(u,v)_{L^{2}(0,T;W)}=\int_{0}^{T}(u(t),v(t))_{W}dt$. Similarly,
$H^{1}(0,T;W)$ is the set of vector-valued functions in $L^{2}(0,T;W)$ with a
time derivative in the weak sense, $u_{t}:[0,T]\to W$, that belongs to
$L^{2}(0,T;W)$. The inner product in $H^{1}(0,T;W)$ is
$(u,v)_{H^{1}(0,T;W)}=\int_{0}^{T}(u,v)_{W}+(u_{t},v_{t})_{W}dt$.

##### Functions of Bounded Variation.

We give a brief summary and state the properties of functions of bounded
variation, which we require in the remainder. For a detailed introduction, we
refer to the reader to the monograph [1]. First, we recall that a function
$f:(0,T)\to\mathbb{R}$ is defined to be of bounded variation or in the space
$\operatorname*{BV}(0,T)$ if $f\in L^{1}(0,T)$ and

$\operatorname*{TV}(f)\coloneqq\sup\left\\{\int_{0}^{T}f(t)\phi^{\prime}(t)\mathop{}\\!\mathrm{d}t\,\Bigg{|}\,\phi\in
C^{1}_{c}(0,T)\ \text{ and }\ \sup_{\mathclap{s\in(0,T)}}|\phi(s)|\leq
1\right\\}<\infty.$

We recall that a sequence $(w^{n})_{n}\subset\operatorname*{BV}(0,T)$ is said
to converge _weakly- ∗_ to a function $w\in\operatorname*{BV}(0,T)$ if
$w^{n}\to w$ in $L^{1}(0,T)$ and
$\limsup_{n\to\infty}\operatorname*{TV}(w^{n})<\infty$. Moreover,
$(w^{n})_{n}$ is said to converge _strictly_ to $w$ if in addition
$\operatorname*{TV}(w^{n})\to\operatorname*{TV}(w)$. We define the subset of
$\operatorname*{BV}(0,T)$ that corresponds to the feasible set of the
optimization problem (P) as

$\operatorname*{BV_{W}}(0,T)=\\{w\in\operatorname*{BV}(0,T)\,|\,w(t)\in W\\}.$

It is important for our analysis that the subset $\operatorname*{BV_{W}}(0,T)$
is closed with respect to weak-∗ and strict convergence in
$\operatorname*{BV}(0,T)$. The closedness follows from the fact that sequences
of $W$-valued functions that converge in $L^{1}(0,T)$ also have $W$-valued
limits, which is stated explicitly for our context in [33, Lemma 2.2].

The analysis of the algorithm in [33]—the starting point of our work—makes use
of regularity conditions, particularly continuity properties, that are defined
for input functions in $L^{2}(0,T)$. While the norm-topologies of $L^{1}(0,T)$
and $L^{2}(0,T)$ are still different, we note that convergence of a sequence
of $W$-valued functions in $L^{1}(0,T)$ implies that the sequence converges in
$L^{2}(0,T)$ as well. This is due to the fact that $W$ is finite, implying a
uniform $L^{\infty}(0,T)$-bound on any sequence of functions, which in turn
implies that we can apply Lebesgue’s dominated convergence theorem to obtain
convergence in $L^{p}(0,T)$ for all $p\in(1,\infty)$ from convergence in
$L^{1}(0,T)$. Consequently, we will frequently use that sequences of functions
that converge weakly-∗ or strictly in $\operatorname*{BV_{W}}(0,T)$ also
converge in $L^{2}(0,T)$.

##  3 Sequential Linear Integer Programming Algorithm

In order to provide a self-contained article, we provide the SLIP algorithm,
which is a function space algorithm to solve (P) to stationarity [33].
Conceptually, the SLIP algorithm is a trust-region method that solves a
sequence of trust-region subproblems. We briefly introduce the trust-region
problem below before laying out the algorithm.

The trust-region subproblem takes a feasible control $w$ for (P) and a trust-
region radius $\Delta>0$ as inputs and reads

$\displaystyle\begin{aligned} \min_{d\in L^{2}(0,T)}&(\nabla
j(w),d)_{L^{2}(0,T)}+\alpha\operatorname*{TV}(w+d)-\alpha\operatorname*{TV}(w)\eqqcolon\ell(w,d)\\\
\text{s.t.}\quad&w(s)+d(s)\in W\text{ for a.a.\ }s\in(0,T),\\\
&\|d\|_{L^{1}(0,T)}\leq\Delta.\end{aligned}$ (TR)

An instance of the problem class (TR) for given $w$ and $\Delta>0$ is denoted
by $\operatorname*{TR}(w,\Delta)$ in the remainder.

We state the trust-region algorithm that solves subproblems of the form (TR)
in Algorithm 1 [33]. The algorithm consists of two nested loops. In every
iteration of the outer loop, which is indexed by $n\in\mathbb{N}$, the trust-
region radius is reset to the input $\Delta^{0}>0$. Then the inner loop, which
is indexed by $k\in\mathbb{N}$, is executed. In each inner iteration, a trust-
region subproblem $\operatorname*{TR}(w^{n-1},\Delta^{n,k})$ is solved for the
current trust-region radius and the previously accepted iterate $w^{n-1}$ or
the input $w^{0}$ (if $n-1=0$). If the predicted reduction (measured as the
objective value of the trust-region subproblem) is zero, then the algorithm
terminates. If the solution of the trust-region subproblem is acceptable (the
ratio of actual reduction and predicted reduction is larger than the input
$\sigma>0$), then the inner loop terminates with new iterate $w^{n}$. If the
step is rejected, then the trust region is reduced and another iteration of
the inner loop executes.

Algorithm 1 Trust-region algorithm from [33]

Input:  Initial guess $w^{0}$ (feasible for (P)), $\Delta^{0}>0$,
$\sigma\in(0,1)$

1:for $n=1,\ldots$ do

2: $k\leftarrow 0$

3: $\Delta^{n,0}\leftarrow\Delta^{0}$

4: repeat

5: $d^{n,k}\leftarrow$ minimizer of $\operatorname*{TR}(w^{n-1},\Delta^{n,k})$
$\triangleright$ Solve trust-region subproblem.

6: if $\ell(w^{n-1},d)=0$ then $\triangleright$ Predicted reduction is zero
$\Rightarrow$ terminate.

7: Terminate with solution $w^{n-1}$.

8: else if
$\frac{j(w^{n-1})+\alpha\operatorname*{TV}(w^{n-1})-j(w^{n-1}+d^{n,k})-\alpha\operatorname*{TV}(w^{n-1}+d^{n,k})}{\ell(w^{n-1},d^{n,k})}<\sigma$
then $\triangleright$ Reject step.

9: $\Delta^{n,k+1}\leftarrow\Delta^{n,k}/2$

10: else$\triangleright$ Accept step.

11: $w^{n}\leftarrow w^{n-1}+d^{n,k}$

12: end if

13: $k\leftarrow k+1$

14: until
$\frac{j(w^{n-1})+\alpha\operatorname*{TV}(w^{n-1})-j(w^{n-1}+d^{n,k})-\alpha\operatorname*{TV}(w^{n-1}+d^{n,k})}{\ell(w^{n-1},d^{n,k})}\geq\sigma$

15:end for

The main known convergence result on Algorithm 1 to this point is that its
iterates converge to so-called L-stationary points under a suitable regularity
assumption on $j$ [33]. A feasible point $w$ is L-stationary if the objective
(P) cannot be improved further by perturbing the locations of its jumps on
$(0,T)$. Such perturbations leave the term $\alpha\operatorname*{TV}(w)$
unchanged and only affect $j(w)$, yielding a condition on $\nabla j(w)$. In
particular, the condition coincides with $\nabla j(w)(t_{i})=0$ for all jump
locations $t_{i}$ of $w$ if $\nabla j(w)$ is a continuous function. The formal
definition is given below.

######  Definition 1.

Let $v\in\operatorname*{BV_{W}}(0,T)$ with representation
$v=\chi_{(t_{0},t_{1})}a_{1}+\sum_{i=1}^{N-1}\chi_{[t_{i},t_{i+1})}a_{i+1}$
for $N\in\mathbb{N}$, $t_{0}=0$, $t_{N}=T$, $t_{i}<t_{i+1}$ for
$i\in\\{0,\ldots,N-1\\}$, and $a_{i}\in W$ for $i\in\\{1,\ldots,N\\}$ be
given. Then $v$ is L-stationary for (P) if

  1. 1.

$\overline{D^{-}_{i}}(\nabla j(v))\geq 0\geq\underline{D_{i}^{+}}(\nabla
j(v))$ if $a_{i}<a_{i+1}$, and

  2. 2.

$\underline{D^{-}_{i}}(\nabla j(v))\leq 0\leq\overline{D_{i}^{+}}(\nabla
j(v))$ if $a_{i+1}<a_{i}$,

where

$\overline{D_{i}^{-}}(\nabla j(v))\coloneqq\limsup_{h\downarrow
0}\frac{1}{h}\int_{t_{i}-h}^{t_{i}}\nabla
j(v)(s)\mathop{}\\!\mathrm{d}s,\quad\underline{D_{i}^{-}}(\nabla
j(v))\coloneqq\liminf_{h\downarrow 0}\frac{1}{h}\int_{t_{i}-h}^{t_{i}}\nabla
j(v)(s)\mathop{}\\!\mathrm{d}s$

and

$\underline{D_{i}^{+}}(\nabla j(v))\coloneqq\liminf_{h\downarrow
0}\frac{1}{h}\int_{t_{i}}^{t_{i}+h}\nabla
j(v)(s)\mathop{}\\!\mathrm{d}s,\quad\overline{D_{i}^{+}}(\nabla
j(v))\coloneqq\limsup_{h\downarrow 0}\frac{1}{h}\int_{t_{i}}^{t_{i}+h}\nabla
j(v)(s)\mathop{}\\!\mathrm{d}s$

for $i\in\\{1,\ldots,N-1\\}$.

Note that Definition 1 is well posed because every function in
$\operatorname*{BV_{W}}(0,T)$ can be written in the claimed form, see, for
example, [33, Proposition 4.4].

##  4 Input Regularization

We introduce the regularity assumption on the Hessian of $j$, under which
convergence to L-stationary points of the iterates produced by Algorithm 1 can
be achieved. Motivated by our applications, we introduce a weaker assumptions
on the Hessian regularity and show that the required regularity assumption for
convergence of Algorithm 1 can always be verified by regularizing (smoothing)
the input of $j$ provided these weaker assumptions hold. For one of these
assumptions, we then prove $\Gamma$-convergence for the case that the
regularization is carried out by a positive mollifier when driving the support
parameter to zero. We also show that the limits of the final iterations of
Algorithm 1 under realistic termination criteria, which are in general no
global minimizers but only (approximately) L-stationary points are strict.

######  Assumption 4.1 (Assumption 4.1 in [33]).

Let $j:L^{2}(0,T)\to\mathbb{R}$ be twice Fréchet differentiable such that for
all $w\in L^{2}(0,T)$

$|\nabla^{2}j(w)(\psi,\phi)|\leq C\|\psi\|_{L^{1}(0,T)}\|\phi\|_{L^{1}(0,T)}$

holds for some $C>0$ and all $\psi,\phi\in L^{2}(0,T)$.

Under Assumption 4.1 on the reduced objective it can be shown that the
iterates produced by Algorithm 1 converge to L-stationary points [33]. In [33]
it is assumed that $j$ is not only twice Fréchet differentiable but the second
derivative $\nabla^{2}j$ is also continuous. This assumption enters the proofs
of Lemma 4.10 and Theorem 4.23 in [33] indirectly through the employed variant
of Taylor’s theorem, Proposition A.1 in [33], which states that for all $u$,
$v\in L^{2}(0,T)$ there exists $\xi$ in the line segment between $u$ and $v$
such that $j(v)=j(u)+(\nabla
j(u),v-u)_{L^{2}(0,T)}+\tfrac{1}{2}(v-u,\nabla^{2}j(\xi)(v-u))_{L^{2}(0,T)}$
holds. This result can be recovered in the absence of the assumption of
continuity of $\nabla^{2}j(\xi)$, which is argued below. This has been
observed in [37, Section 3.2], where a slightly different formulation of the
Taylor expansion is used. 111The authors thank Gerd Wachsmuth for the hint to
Darboux’s theorem and [2].

######  Proposition 4.2.

Let $j:L^{2}(0,T)\to\mathbb{R}$ be twice Fréchet differentiable. Let $w$,
$\bar{w}\in L^{2}(0,T)$ be given. Then there exists $\xi=w+\tau(\bar{w}-w)$
for some $\tau\in[0,1]$ such that $j(\bar{w})=j(w)+(\nabla
j(w),\bar{w}-w)_{L^{2}(0,T)}+\tfrac{1}{2}(\bar{w}-w,\nabla^{2}j(\xi)(\bar{w}-w))_{L^{2}(0,T)}$.

###### Proof.

We reduce to the finite-dimensional case and consider the function
$\tilde{j}:[0,1]\ni t\mapsto j(w+t(\bar{w}-w))\in\mathbb{R}$. The chain rule
in Banach spaces, see, for example, [49, Theorem 4.D] implies that
$\tilde{j}^{\prime\prime}(t)=(\bar{w}-w,\nabla^{2}j(w+t(\bar{w}-w))(\bar{w}-w))_{L^{2}(\Omega)}$
for $t\in[0,1]$. Then we apply a variant of Taylor’s theorem that is based on
Darboux’s theorem and does not require continuity of the second derivative,
see, for example, [2, Theorem 5.19], to obtain that there exists
$\tau\in[0,1]$ such that $j(\bar{w})=j(w)+(\nabla
j(w),\bar{w}-w)_{L^{2}(\Omega)}+\tfrac{1}{2}(\bar{w}-w,\nabla^{2}j(w+\tau(\bar{w}-w))(\bar{w}-w))_{L^{2}(\Omega)}$,
which proves the claim with the choice $\xi\coloneqq w+\tau(\bar{w}-w)$. ∎

While Assumption 4.1 is required for the convergence analysis of Algorithm 1
in [33], it may be considered to be too restrictive for practical applications
of Algorithm 1 because it requires an improvement of the input regularity of
the control-to-state operator. In particular, as we also experience for the
considered poro-visco-elastic problem in section 5, one may only able to
verify one of the weaker assumptions below, where the uniform boundedness of
the Hessian of $j$ is assumed with respect to stronger norms for the control
input.

######  Assumption 4.3.

Let $j:L^{2}(0,T)\to\mathbb{R}$ be twice Fréchet differentiable such that for
all $w\in L^{2}(0,T)$

$|\nabla^{2}j(w)(\psi,\phi)|\leq C\|\psi\|_{L^{2}(0,T)}\|\phi\|_{L^{2}(0,T)}$

holds for some $C>0$ and all $\psi,\phi\in L^{2}(0,T)$.

######  Assumption 4.4.

Let $j:H^{1}(0,T)\to\mathbb{R}$ be twice Fréchet differentiable such that for
all $w\in H^{1}(0,T)$

$|\nabla^{2}j(w)(\psi,\phi)|\leq C\|\psi\|_{H^{1}(0,T)}\|\phi\|_{H^{1}(0,T)}$

holds for some $C>0$ and all $\psi,\phi\in H^{1}(0,T)$.

We propose to modify $j$ as follows. We replace $j$ by

$j_{\varepsilon}\coloneqq j\circ K_{\varepsilon},$

where $K_{\varepsilon}:L^{1}(0,T)\to L^{2}(0,T)$ is a (bounded and linear)
convolution operator defined as $K_{\varepsilon}(w)\coloneqq
r_{[0,T]}(\eta_{\varepsilon}*w)$ for $w\in L^{1}(0,T)$, and
$(\eta_{\varepsilon})_{\varepsilon>0}$ is a family of positive mollifiers,
see, for example, [35, Section 4.5], and $r_{[0,T]}$ denotes the restriction
of a function that is defined on all of $\mathbb{R}$ to the interval $[0,T]$.
Then the following three propositions hold.

######  Proposition 4.5.

Let $\varepsilon>0$. Then $K_{\varepsilon}:L^{1}(0,T)\to L^{2}(0,T)$ is a
bounded linear operator. If $j$ satisfies Assumption 4.3, then
$j_{\varepsilon}$ satisfies Assumption 4.1.

###### Proof.

Let $w\in L^{1}(0,T)$. Let $K_{\varepsilon}^{*}:L^{2}(0,T)\to L^{\infty}(0,T)$
denote the adjoint operator of the bounded linear operator $K_{\varepsilon}$,
where we have identified $(L^{2}(0,T))^{*}\cong L^{2}(0,T)$ and
$(L^{1}(0,T))^{*}\cong L^{\infty}(0,T)$. The chain rule, see, [49, Theorem
4.D], yields the derivatives

$Dj_{\varepsilon}(w)=K_{\varepsilon}^{*}\nabla
j(K_{\varepsilon}w)\quad\text{and}\quad D^{2}j_{\varepsilon}(w)=\langle
K_{\varepsilon}^{*}\nabla^{2}j(K_{\varepsilon}w)K_{\varepsilon}\cdot,\cdot\rangle_{L^{\infty}(0,T),L^{1}(0,T)}.$

Let $\phi$, $\psi\in L^{1}(0,T).$ Cauchy–Schwarz inequality and the
submultiplicativity of the operator norm give

$\langle
K_{\varepsilon}^{*}\nabla^{2}j(K_{\varepsilon}w)K_{\varepsilon}\psi,\phi\rangle_{L^{\infty}(0,T),L^{1}(0,T)}=(\nabla^{2}j(K_{\varepsilon}w)K_{\varepsilon}\psi,K_{\varepsilon}\phi)_{L^{2}(0,T)}\leq
C\|K_{\varepsilon}\psi\|_{L^{2}(0,T)}\|K_{\varepsilon}\phi\|_{L^{2}(0,T)}$

for the constant $C>0$ from Assumption 4.3. By virtue of Young’s convolution
inequality, we obtain
$\|K_{\varepsilon}\phi\|_{L^{2}(0,T)}\leq\|\eta_{\varepsilon}\|_{L^{2}(\mathbb{R})}\|\phi\|_{L^{1}(0,T)}$.
Because $(\eta_{\varepsilon})_{\varepsilon>0}$ is a family of mollifiers we
have $\|\eta_{\varepsilon}\|_{L^{2}(\mathbb{R})}<\infty$ for all
$\varepsilon>0$ and the claim follows. ∎

######  Proposition 4.6.

Let $\varepsilon>0$. Then $K_{\varepsilon}:L^{1}(0,T)\to H^{1}(0,T)$ is a
bounded linear operator. If $j$ satisfies Assumption 4.4, then
$j_{\varepsilon}$ satisfies Assumption 4.1.

###### Proof.

We can follow the same proof as for Proposition 4.5 above. We need to show
that

$\|K_{\varepsilon}\phi\|_{H^{1}}=\left(\|K_{\varepsilon}\phi\|_{L^{2}}^{2}+\left\|\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}K_{\varepsilon}\phi\right\|_{L^{2}}^{2}\right)^{\frac{1}{2}}\leq
C\|\phi\|_{L^{1}}.$

for some $C>0$. Proposition 4.5 already implies
$\|K_{\varepsilon}\phi\|_{L^{2}}\leq c_{1}\|\phi\|_{L^{1}}$ For the second
estimate we use the computation rules for the derivative of a convolution to
give
$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}(\eta_{\varepsilon}*\phi)=\left(\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\eta_{\varepsilon}\right)*\phi$.
so that
$\left\|\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}K_{\varepsilon}\phi\right\|_{L^{2}}\leq
c_{2}\|\phi\|_{L^{1}}$ holds with Young’s convolution inequality and
$c_{2}\coloneqq\left\|\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\eta_{\varepsilon}\right\|_{L^{2}}$,
which is bounded because $\eta_{\varepsilon}$ is smooth. ∎

######  Proposition 4.7.

Let $j_{\varepsilon}$ for $\varepsilon>0$ be defined as above for a family of
mollifiers $(\eta_{\varepsilon})_{\varepsilon>0}$ with
$\operatorname*{supp}\eta_{1}\subset[-1,1]$. Let $w\in L^{2}(0,T)$. Then
$Dj_{\varepsilon}(w)=K_{\varepsilon}^{*}\nabla j(K_{\varepsilon}w)\in
L^{\infty}(0,T)$ is a continuous function on $[0,T]$.

###### Proof.

As above, the identification $L^{\infty}(0,T)\cong(L^{1}(0,T))^{*}$ implies
that $K_{\varepsilon}^{*}\nabla j(K_{\varepsilon}w)$ can be represented as an
$L^{\infty}(0,T)$-function. Moreover, $\nabla j(K_{\varepsilon}w)$ is an
$L^{2}(0,T)$-function. Let $f\in L^{1}(0,T)$, $g\in L^{2}(0,T)$. We consider
the adjoint operator $K_{\varepsilon}^{*}$, which is defined by the identity

$(K_{\varepsilon}f,g)_{L^{2}(0,T)}=\langle
f,K_{\varepsilon}^{*}g\rangle_{L^{1}(0,T),L^{\infty}(0,T)}.$

We insert the definition of the convolution and obtain the equivalent identity

$\int_{0}^{T}\overbrace{\int_{0}^{T}\eta_{\varepsilon}(t-s)f(s)\mathop{}\\!\mathrm{d}s}^{(K_{\varepsilon}f)(t)}g(t)\mathop{}\\!\mathrm{d}t=\int_{0}^{T}f(s)\overbrace{\int_{0}^{T}\eta_{\varepsilon}(t-s)g(t)\mathop{}\\!\mathrm{d}t}^{(K_{\varepsilon}^{*}g)(s)}\mathop{}\\!\mathrm{d}s$

by virtue of Fubini’s theorem. Due to the structural symmetry, the proof for
smoothness of $K_{\varepsilon}f$ on $[0,T]$ can be transferred directly to
$K_{\varepsilon}^{*}g$ on $[0,T]$. ∎

We are interested in the ability of local solutions, global solutions, and
stationary points of (P) with objective $j_{\varepsilon}$ instead of $j$ to
approximate local solutions, global solutions, and stationary points of (P)
for $\varepsilon\to 0$. While we are not able to give a full answer at this
point, we can prove the case of global solutions by virtue of
$\Gamma$-convergence. Assume that the objectives
$j_{\varepsilon}+\alpha\operatorname*{TV}$ $\Gamma$-converge to
$j+\alpha\operatorname*{TV}$ for $\varepsilon\to 0$ with respect to
$\operatorname*{BV}(0,T)$-weak-∗ convergence on the complete subspace
$\operatorname*{BV_{W}}(0,T)$, which is the feasible set of (P). Then we
obtain that global minimizers of $j_{\varepsilon}$ in
$\operatorname*{BV_{W}}(0,T)$ weakly-∗ converge to global minimizers of $j$ in
$\operatorname*{BV_{W}}(0,T)$, which is one of the the desired results.

######  Theorem 4.8.

Let Assumption 4.3 hold. Let $K_{\varepsilon}$ for $\varepsilon>0$ be defined
as above for a family of positive mollifiers
$(\eta_{\varepsilon})_{\varepsilon>0}$. Then the sequence
$(\hat{j}_{\varepsilon^{n}})_{n\in\mathbb{N}}$, defined as
$\hat{j}_{\varepsilon^{n}}\coloneqq
j_{\varepsilon^{n}}+\alpha\operatorname*{TV}$, $\Gamma$-converges to
$\hat{j}=j+\alpha\operatorname*{TV}$ for $\varepsilon^{n}\to 0$ with respect
to weak-∗ convergence and strict convergence in $\operatorname*{BV_{W}}(0,T)$.

###### Proof.

We start by proving the lower bound inequality. To this end, let
$w^{\varepsilon^{n}}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}w$ in
$\operatorname*{BV_{W}}(0,T)$. Then
$\operatorname*{TV}(w)\leq\liminf_{\varepsilon^{n}\to
0}\operatorname*{TV}(w^{\varepsilon_{n}})$ by virtue of the weak-∗ lower
semicontinuity of the total variation. For all $n\in\mathbb{N}$, Taylor’s
theorem gives

$|j_{\varepsilon^{n}}(w^{\varepsilon^{n}})-j(w)|\leq|j_{\varepsilon^{n}}(w)+\langle
K_{\varepsilon^{n}}^{*}\nabla
j(K_{\varepsilon^{n}}\xi^{\varepsilon^{n}}),w^{\varepsilon^{n}}-w\rangle_{L^{\infty}(0,T),L^{1}(0,T)}-j(w)|$

for some $\xi^{\varepsilon^{n}}$ in the line segment between $w$ and
$w^{\varepsilon^{n}}$. Since $K_{\varepsilon^{n}}w=\eta_{\varepsilon^{n}}*w\to
w$ in $L^{2}(0,T)$, and $\hat{j}$ is continuous (implied by the
differentiability of $j$), we have
$\hat{j}_{\varepsilon^{n}}(w)=j(K_{\varepsilon^{n}}w)\to F(w)$. The set of
functions $(\xi^{\varepsilon^{n}})_{n\in\mathbb{N}}$ is bounded because
$w^{\varepsilon^{n}}\to w$, and thus $(\nabla
j(\xi^{\varepsilon^{n}}))_{n\in\mathbb{N}}\subset L^{2}(0,T)$ is bounded.
Moreover, with
$\|K_{\varepsilon}\|_{\mathcal{L}(L^{1}(0,T),L^{2}(0,T))}\leq\|\eta_{\varepsilon}\|_{L^{1}(0,T)}=1$,
Young’s convolution inequality, and the fact that the operator norms of
bounded linear operators and their adjoints coincide, we obtain that
$(K_{\varepsilon^{n}}^{*}\nabla
j(\xi^{\varepsilon^{n}}))_{n\in\mathbb{N}}\subset L^{\infty}(0,T)$ is bounded.
Moreover, $w^{\varepsilon^{n}}\to w$ in $L^{1}(0,T)$ by assumption, and we
have $\langle K_{\varepsilon^{n}}^{*}\nabla
j(\xi^{\varepsilon^{n}}),w^{\varepsilon^{n}}-w\rangle_{L^{\infty}(0,T)}\to 0$,
which proves the lower bound inequality.

To prove the upper bound inequality, we choose $w^{\varepsilon^{n}}\coloneqq
w$ for all $n\in\mathbb{N}$ and observe

$\hat{j}_{\varepsilon^{n}}(w^{\varepsilon^{n}})=j_{\varepsilon^{n}}(w)+\alpha\operatorname*{TV}(w)\to
j(w)+\alpha\operatorname*{TV}(w)$

with the same argument as above. This proves $\Gamma$-convergence with respect
to weak-∗ convergence in $\operatorname*{BV_{W}}(0,T)$. Because the chosen
sequence for the upper bound inequality is also strictly convergent,
$\Gamma$-convergence also holds with respect to strict convergence in
$\operatorname*{BV_{W}}(0,T)$. ∎

Proposition 4.5 shows it is sufficient to verify the much weaker assumptions
Assumption 4.3 or Assumption 4.4 instead of Assumption 4.1 when we solve (P)
with the input regularization $j_{\varepsilon}$ instead of $j$ for some
$\varepsilon>0$. Proposition 4.7 implies that $\nabla j_{\varepsilon}(w)$,
$w\in\operatorname*{BV_{W}}(0,T)$, is a smooth function on the interval
$[0,T]$. Thus we can use the characterization

$\nabla j_{\varepsilon}(w)(t_{i})=0$

for all of the finitely many $t_{i}$, where the function $w$ has a jump, if
$w$ is L-stationary for (P) with $j_{\varepsilon}$ instead of $j$, see [33].

While we obtain convergence of global minimizers under Assumption 4.3, we do
not know at present if we converge to an L-stationary point of (P) when we
compute L-stationary points for (P) with $j_{\varepsilon}$ instead of $j$ and
drive $\varepsilon$ to zero. In particular, we have not been able to show that
the limit of L-stationary or approximately L-stationary points of (P) with
$j_{\varepsilon}$ that are produced by Algorithm 1 for a homotopy that drives
$\varepsilon\to 0$ is L-stationary for (P). We note that even if this were
true, L-stationary is not known to be a necessary optimality condition for (P)
if only Assumption 4.3 but not Assumption 4.1 holds.

Moreover, the situation is even worse in case that we only have Assumption
4.4. Because $H^{1}(0,T)\hookrightarrow C([0,T])$, the only $W$-valued
functions in $H^{1}(0,T)$ are constant on the whole the domain $(0,T)$, so we
cannot prove a result like Theorem 4.8 in this case and the limit problem
$\varepsilon=0$ has no meaningful interpretation due to the high regularity
that is required for the control input. However, in case one is still
interested in discrete-valued controls even if this is not covered by the
regularity theory for the PDE, we believe that solving (P) with
$j_{\varepsilon}$ instead of $j$ is still sensible because for a given control
$w\in\operatorname*{BV_{W}}(0,T)$ we know that there are only finitely many
switches (or jumps) that can occur and by means of the parameter $\varepsilon$
we can control the support size of the smooth transitions between them that
occurs when mollifying them.

However, the succeeding analysis shows that in case of Assumption 4.3, weak-∗
limit points of the homotopy are also strict limit points, which means that
the limit cannot have a lower total variation than its approximating sequence.
Consequently, there is no nearby reduction of the objective function by an
improvement of the $\operatorname*{TV}$-term. This result is possible even in
the presence of realistic early termination criteria in Algorithm 1.

######  Proposition 4.9.

Let $j$ be bounded below. Assume that Algorithm 1 is terminated when one of
the following conditions is met:

  * •

the trust-region radius is smaller than a given $C_{1}>0$

  * •

the predicted reduction, $\ell(w^{n-1},\Delta^{n,k})$ in outer iteration
$n\in\mathbb{N}$ and inner iteration $k\in\mathbb{N}$, is smaller than a given
$C_{2}>0$,

then it terminates within finitely many outer iterations regardless of the
initial guess $w^{0}$.

######  Remark 4.1.

Before proving Proposition 4.9, we make a brief note to explicitly explain how
these conditions enter Algorithm 1. The first condition is checked after each
reduction of the trust-region radius together with Algorithm 1 ln. 14. The
second condition replaces the termination criterion in Algorithm 1 ln. 6.

###### Proof of Proposition 4.9.

The first condition ensures that the inner loop iterates only finitely many
times, specifically at most $\lfloor\log_{2}(C_{1}/\Delta^{0})\rfloor$ times.
Assume that Algorithm 1 does not terminate within finitely many iterations of
the outer loop. Then the inner loop accepts a new iterate within
$k\leq\lfloor\log_{2}(C_{1}/\Delta^{0})\rfloor$ iterations in each outer
iteration $n\in\mathbb{N}$. The step acceptance implies that the reduction in
the objective is always higher than $\sigma\ell(w^{n-1},\Delta^{n,k})$. It
follows from the second termination criterion that
$\ell(w^{n-1},\Delta^{n,k})>C_{2}$ for all $n\in\mathbb{N}$ and corresponding
$k$ on acceptance. This contradicts that $j$ and $\operatorname*{TV}$ are
bounded from below. Consequently, Algorithm 1 terminates within finitely many
outer iterations. ∎

######  Proposition 4.10.

Let $j$ be bounded below. Let Assumption 4.3 hold. Let
$(w_{\varepsilon})_{\varepsilon}\subset\operatorname*{BV_{W}}(\Omega)$ be the
sequence of final iterates produced by Algorithm 1 executed on
$j_{\varepsilon}$ for a sequence $\varepsilon\to 0$ with the initial control
given by the previous final iterate, where the execution of Algorithm 1 is
terminated if one of following conditions is met:

  * •

the trust region radius is smaller than $\Delta(\varepsilon)>0$, which tends
to zero as $\varepsilon$ is driven to zero

  * •

the predicted reduction is smaller than $\sigma(1-\frac{3a}{4})\alpha>0$ for
some fixed $0<a\leq 1-\sigma<1$.

If $\liminf_{\varepsilon\to 0}\operatorname*{TV}(w_{\varepsilon})<\infty$,
there is at least one weakly-∗ convergent subsequence. Every weakly-∗
convergent subsequence converges strictly to a limit point in
$\operatorname*{BV_{W}}(\Omega)$.

###### Proof.

First, we note that the sequence $(w_{\varepsilon})_{\varepsilon}$ is well
defined because Proposition 4.9 asserts that Algorithm 1 terminates after
finitely many iterations for the two assumed termination criteria. Second, the
fact that $\liminf_{\varepsilon\to
0}\operatorname*{TV}(w_{\varepsilon})<\infty$ implies that there is at least
one weakly-∗ convergent subsequence follows by definition of weak-∗
convergence in $\operatorname*{BV}(0,T)$.

In the remainder of the proof we consider an arbitrary weak-∗ convergent
subsequence $w_{i}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\bar{w}$ for
$\varepsilon_{i}\to 0$ in $\operatorname*{BV}(0,T)$. The convergence theory of
Algorithm 1 in [33] gives
$(w_{\varepsilon})_{\varepsilon}\subset\operatorname*{BV_{W}}(0,T)$ and the
weak-∗ closedness of $\operatorname*{BV_{W}}(0,T)$ (see section 2) gives
$\bar{w}\in\operatorname*{BV_{W}}(0,T)$. We note that the subsequence also
converges in $L^{2}(0,T)$ by boundedness in $L^{\infty}(0,T)$ (due to
$w_{i}(t)\in W$ a.e. for all $\varepsilon_{i}$) and pointwise convergence a.e.
of a subsequence.

We set forth to prove the claim with a contradictory argument. To this end,
assume that the convergence of the subsequence is not strict. Then we can find
a subsequence, for ease of notation also denoted by $w_{i}$, and $\bar{C}>0$
such that

$\displaystyle\alpha\operatorname*{TV}(w_{i})-\alpha\operatorname*{TV}(\bar{w})\geq\bar{C}>0$

for all $i\in\mathbb{N}$. Note

$\displaystyle|j(K_{\varepsilon_{i}}w_{i})-j(\bar{w})|\leq|j(K_{\varepsilon_{i}}w_{i})-j(K_{\varepsilon_{i}}\bar{w})|+|j(K_{\varepsilon_{i}}\bar{w})-j(\bar{w})|.$

Additionally, $\|K_{\varepsilon_{i}}\bar{w}-\bar{w}\|_{L^{2}}\to 0$ because
$\sup_{\|\phi\|_{L^{2}}=1}\|K_{\varepsilon_{i}}\phi\|_{L^{2}}\leq\|\eta_{\varepsilon_{i}}\|_{L^{1}}=1$
by construction of the mollifier inducing $K_{\varepsilon_{i}}$. Young’s
convolution inequality shows

$\displaystyle\|K_{\varepsilon_{i}}w_{i}-K_{\varepsilon_{i}}\bar{w}\|_{L^{2}}\leq\|\eta_{\varepsilon_{i}}\|_{L^{1}}\|w_{i}-\bar{w}\|_{L^{2}}=\|w_{i}-\bar{w}\|_{L^{2}}\to
0.$

Thus, we obtain that $j(K_{\varepsilon_{i}}w_{i})\to j(\bar{w})$ for
$i\to\infty.$

Let $\tilde{\Delta}>0$ be arbitrary but fixed. Then there is
$i_{0}\in\mathbb{N}$ large enough such that for all $i\geq i_{0}$

$\displaystyle\|w_{i}-\bar{w}\|_{L^{1}}$
$\displaystyle\leq\tilde{\Delta}\text{ and}$ $\displaystyle\|\nabla
j(w_{i})-\nabla j(\bar{w})\|_{L^{2}}$ $\displaystyle\leq\tilde{\Delta}$

hold. Thus $\bar{w}$ is feasible for
$\operatorname*{TR}(w_{i},\tilde{\Delta})$ for all $i\geq i_{0}$. Let
$w_{0}\in\operatorname*{BV_{W}}(0,T)$ be an arbitrary feasible point for
$\operatorname*{TR}(w_{i},\tilde{\Delta})$. Using Taylor’s theorem we obtain
that

$\displaystyle j(K_{\varepsilon_{i}}w_{i})-j(K_{i}w_{0})$
$\displaystyle=\sigma(\nabla
j(K_{\varepsilon_{i}}w_{i}),K_{\varepsilon_{i}}(w_{i}-w_{0}))_{L^{2}}+(1-\sigma)(\nabla
j(K_{\varepsilon_{i}}w_{i}),K_{\varepsilon_{i}}(w_{i}-w_{0}))_{L^{2}}$
$\displaystyle\hphantom{=}+\frac{1}{2}\nabla^{2}j(\xi_{i})(K_{\varepsilon_{i}}(w_{i}-w_{0}),K_{\varepsilon_{i}}(w_{i}-w_{0}))$

for some $\xi_{i}\in L^{2}(0,T)$. From Assumption 4.3 we derive that

$\displaystyle|\nabla^{2}j(\xi_{i})(K_{\varepsilon_{i}}(w_{i}-w_{0}),K_{\varepsilon_{i}}(w_{i}-w_{0}))|\leq
C\|w_{i}-w_{0}\|_{L^{2}}^{2}\leq C|\max(W)-\min(W)|\|w_{i}-w_{0}\|_{L^{1}}\leq
C_{2}\tilde{\Delta}$

for some $C$, $C_{2}>0$. Furthermore, we obtain that

$\displaystyle|(\nabla
j(K_{\varepsilon_{i}}w_{i}),K_{\varepsilon_{i}}(w_{i}-w_{0}))_{L^{2}}|$
$\displaystyle\leq\|\nabla
j(K_{\varepsilon_{i}}w_{i})\|_{L^{2}}\|K_{\varepsilon_{i}}(w_{i}-w_{0})\|_{L^{2}}$
$\displaystyle\leq(\|\nabla
j(\bar{w})\|_{L^{2}}+C_{3}\tilde{\Delta})\|w_{i}-w_{0}\|_{L^{2}}$
$\displaystyle\leq C_{4}(C_{5}+C_{3}\tilde{\Delta})\sqrt{\tilde{\Delta}}$

holds for some $C_{3}$, $C_{4}$, $C_{5}>0$. Then there exists $\tilde{\Delta}$
as above such that for $i_{0}$ corresponding to $\tilde{\Delta}$ we obtain for
all $i\geq i_{0}$ that the estimates

$\displaystyle\Delta(\varepsilon_{i})$
$\displaystyle<\frac{\tilde{\Delta}}{2},$ (4.1)
$\displaystyle|\nabla^{2}j(\xi_{i})(K_{\varepsilon_{i}}(w_{i}-w_{0}),K_{\varepsilon_{i}}(w_{i}-w_{0}))|$
$\displaystyle\leq\frac{a}{4}\bar{C}\text{ and }$ (4.2) $\displaystyle|(\nabla
j(K_{\varepsilon_{i}}w_{i}),K_{\varepsilon_{i}}(w_{i}-w_{0}))_{L^{2}}|$
$\displaystyle\leq\frac{a}{4}\bar{C}$ (4.3)

hold for all $w_{0}$ in the trust region of $\operatorname*{TR}(w_{i},\Delta)$
if $\Delta\leq\tilde{\Delta}$. Moreover, there exists a minimal
$\tilde{k}\in\mathbb{N}$ such that
$\tilde{\Delta}\geq\Delta^{0}2^{-\tilde{k}}$.

We now show that the inner loop of Algorithm 1 latemost accepts the step at
inner iteration $\tilde{k}$ for $i$ large enough. Specifically, for all $i\geq
i_{0}$, the execution of Algorithm 1 for $\varepsilon_{i}$ latemost accepts a
step in inner iteration $\tilde{k}$ in the final outer iteration. The trust-
region radius upon acceptance is $\Delta^{0}2^{-\tilde{k}}$ with
$\Delta(\varepsilon_{i})<\Delta^{0}2^{-\tilde{k}}\leq\tilde{\Delta}\leq\Delta^{0}2^{-\tilde{k}+1}$.

Let $w^{*}$ be the optimal solution of
$\operatorname*{TR}(w_{i},\Delta^{0}2^{-\tilde{k}})$. Because $\bar{w}$ is in
the trust region of $\operatorname*{TR}(w_{i},\Delta^{0}2^{-\tilde{k}})$ it
follows that

$\displaystyle\alpha(\operatorname*{TV}(w_{i})-\operatorname*{TV}(w^{*}))\geq\alpha(\operatorname*{TV}(w_{i})-\operatorname*{TV}(\bar{w}))-\frac{a}{2}\bar{C}\geq\left(1-\frac{a}{2}\right)\bar{C}$
(4.4)

because the predicted reduction of $w^{*}$ is greater or equal than of
$\bar{w}$. Let
$M\coloneqq\alpha(\operatorname*{TV}(w_{\varepsilon_{i}})-\operatorname*{TV}(w^{*}))-(1-\frac{a}{2})\bar{C}\geq
0$. In total we obtain the inequalities

$\displaystyle
j(K_{\varepsilon_{i}}w_{i})-j(K_{\varepsilon_{i}}w^{*})+\alpha(\operatorname*{TV}(w_{i})-\operatorname*{TV}(w^{*}))$
$\displaystyle=\sigma(\nabla
j(K_{\varepsilon_{i}}w_{i}),K_{\varepsilon_{i}}(w_{i}-w^{*}))_{L^{2}}+(1-\sigma)(\nabla
j(K_{\varepsilon_{i}}w_{i}),K_{\varepsilon_{i}}(w_{i}-w^{*}))_{L^{2}}$
$\displaystyle\hphantom{=}\
+\frac{1}{2}\nabla^{2}j(\xi)(K_{\varepsilon_{i}}(w_{i}-w^{*}),K_{\varepsilon_{i}}(w_{i}-w^{*}))+\alpha(\operatorname*{TV}(w_{i})-\operatorname*{TV}(w^{*}))$
$\displaystyle\geq-\sigma\ell(w_{i},w^{*}-w_{i})-\frac{a}{4}\bar{C}-\frac{a}{4}\bar{C}+(1-\sigma)\left(1-\frac{a}{2}\right)\bar{C}+(1-\sigma)M$
(4.2),(4.3)
$\displaystyle\geq-\sigma\ell(w_{i},w^{*}-w_{i})-\frac{a}{4}\bar{C}-\frac{a}{4}\bar{C}+\frac{a}{2}\bar{C}+(1-\sigma)M$
$a\leq 1-\sigma$ $\displaystyle\geq-\sigma\ell(w_{i},w^{*}-w_{i})$
$\displaystyle\geq\sigma\left(1-\frac{3a}{4}\right)\bar{C}\underset{a<1}{>}0.$
(4.3),(4.4)

We analyze these inequalities with respect to the two possible termination
criteria that are assumed (note that the original termination criterion in
Algorithm 1 ln. 6 is replaced by the second one, see also remark 4.1).

The first termination criterion does not apply, because the trust-region
algorithm would find an improvement of $\sigma(1-\frac{3a}{4})\bar{C}$ before
the critical trust-region radius is attained. The second termination criterion
does not apply, because the predicted reduction would have to be less than
$\sigma(1-\frac{3a}{4})\alpha\leq\sigma(1-\frac{3a}{4})\bar{C}$, which can not
happen as the inequalities show.

Thus neither of the assumed termination criteria is satisfied for any $i\geq
i_{0}$ but the $w_{i}$ are final iterates of Algorithm 1 under the assumed
termination criteria, which is a contradiction. Consequently, the convergence
of $w_{i}$ to $\bar{w}$ is strict. ∎

######  Remark 4.2.

We note that the proof does not require the assumption that an execution of
Algorithm 1 is initialized with the final iterate of the execution for the
previous choice of $\varepsilon$. We have included it here because it is the
natural choice for a homotopy that drives $\varepsilon$ to zero and generally
helps to obtain a bounded (and improving) sequence in practice as can also be
seen in our numerical examples.

##  5 Application to Fluid Flows through Deformable, Porous Media

PDE Model. Let $\Omega\subset\mathbb{R}^{3}$ be the open, bounded domain
occupied by the fluid-solid mixture, with Lipschitz boundary $\partial\Omega$.
Motivated by applications in biomechanics (like tissue perfusion [3, 18, 29,
31, 39, 11]), we work under the assumptions of full saturation, negligible
inertia, small deformations and incompressible mixture components (in the
sense that the solid and fluid phases can’t undergo volume changes at the
microscale). Due to the complex composition of biological tissue, which
exhibit both elastic and viscoelastic behaviors, we consider both poro-elastic
and poro-visco-elastic systems, where the effective stress tensor is of
Kelvin-Voigt type. The extent to which structural viscoelasticity is present
in the equations is represented by the parameter $\delta\geq 0$. Therefore,
the total stress of the fluid-solid mixture is given by

$\mathbf{T}(\mathbf{u},p)=\delta\mu_{v}(\nabla\mathbf{u}_{t}+\nabla\mathbf{u}_{t}^{T})+\delta\lambda_{v}(\nabla\cdot\mathbf{u}_{t})\mathbf{I}+\mu_{e}(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})+\lambda_{e}(\nabla\cdot\mathbf{u})\,\mathbf{I}-p\mathbf{I},$

where $\mathbf{u}$ is the elastic displacement and $p$ is the fluid pressure.
Moreover, $\mathbf{I}$ stands for the identity tensor, and $\lambda_{v}$,
$\mu_{v}$, $\lambda_{e}$ and $\mu_{e}$ are the Lamé parameters which are all
strictly positive.

The quasi-static system is described by two conservation laws: the balance of
linear momentum for the fluid-solid mixture and the balance of mass for the
fluid component.

$\nabla\cdot\mathbf{T}(\mathbf{u},p)+\mathbf{F}(\mathbf{x},t)=\mathbf{0}\quad\mbox{and}\quad\zeta_{t}+\nabla\cdot\mathbf{v}=S(\mathbf{x},t)\quad\mbox{in}\;\Omega\times(0,T),$
(5.1)

where the fluid content $\zeta$ is given by
$\displaystyle\zeta=\nabla\cdot\mathbf{u}$, and the discharge velocity
$\mathbf{v}$ is given by $\mathbf{v}=-\mathbf{K}\nabla p,$ with
$\mathbf{K}=k\mathbf{I}$, where $\mathbf{K}$ is the permeability tensor and
$k$ is a constant.We note here that the formula for the fluid content is a
simplification of the more general expression
$\zeta=c_{0}p+\alpha\nabla\cdot\mathbf{u}$ [6], where $c_{0}$ is the
constrained specific storage coefficient and $\alpha$ is the Biot-Willis
coefficient. The simplification is made due to the assumption of
incompressible fluid and solid components of the mixture (biological tissues
have a mass density close to that of the water), $c_{0}=0$ and $\alpha=1$
[15], and therefore the fluid content becomes solid dilation.

We assume that $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$, where $\Gamma_{D}$
and $\Gamma_{N}$ are the Dirichlet and Neumann parts of the boundary (with
respect to the elastic displacement), with
$\Gamma_{D}\cap\Gamma_{N}=\emptyset$ (while allowing
$\overline{\Gamma}_{D}\cap\overline{\Gamma}_{N}\neq\emptyset$).

We associate the following boundary and initial conditions to the balance laws
mentioned above:

$\mathbf{T}(\mathbf{u},p)\mathbf{n}=\mathbf{g},\quad\mathbf{v}\cdot\mathbf{n}=0\quad\mbox{on}\;\Gamma_{N}\,,$
(5.2) $\mathbf{u}=\mathbf{0},\quad p=0\quad\mbox{on}\;\Gamma_{D,p}\,,$ (5.3)
$\mathbf{u}=\mathbf{0},\quad\mathbf{v}\cdot\mathbf{n}=\psi\quad\mbox{on}\;\Gamma_{D,\mathbf{v}}\,.$
(5.4) $\displaystyle\mathbf{u}(x,0)=\mathbf{u}_{0}(x)\ \text{in}\ \Omega.$
(5.5)

Note that the Dirichlet part of the boundary
$\Gamma_{D}=\Gamma_{D,p}\cup\Gamma_{D,\mathbf{v}}$, where the subscripts $p$
and $\mathbf{v}$ indicate conditions imposed on the Darcy pressure and
discharge velocity, respectively. As usual, $\mathbf{n}$ is the outward unit
normal vector.

The data in the system is represented by the body force per unit of volume
$\mathbf{F}$, the net volumetric fluid production rate $S$, and the boundary
sources $\mathbf{g}$ and $\psi$. They can be used as controls.

We impose the following assumptions on the domain:

######  Assumption 5.1.

We assume:

  1. 1.

$\Gamma_{D}$ is a set of positive measure, so by Korn’s inequality:

$E(\mathbf{u}(t))\geq c||\mathbf{u}(t)||^{2}_{\mathbf{H}^{1}(\Omega)}.$

  2. 2.

$\Gamma_{D,p}$ is a set of positive measure, so by Poincare’s inequality:

$||v||_{L^{2}(\Omega)}\leq C_{P}||\nabla
v||_{\mathbf{L}^{2}(\Omega)},~{}~{}\forall v\in V.$

Notation. Let

$\mathbb{V}\equiv\mathbf{V}\times V=(H^{1}_{\Gamma_{D}}(\Omega))^{3}\times
H^{1}_{\Gamma_{D,p}}(\Omega),$

where the inner-product on $\mathbf{V}$ is given by

$a(\mathbf{u},\mathbf{w})=(\nabla\cdot\mathbf{u},\nabla\cdot\mathbf{w})+(\nabla\mathbf{u}:\nabla\mathbf{w})+(\nabla\mathbf{u}:(\nabla\mathbf{w})^{T}).$
(5.6)

Here $\nabla\mathbf{u}$ stands for the Jacobian matrix of $\mathbf{u}$ and the
Frobenius inner product of two matrices is given by

$(\mathbf{A}:\mathbf{B})=\int_{\Omega}(A_{ij}B_{ij})d\Omega.$

We note that the bilinear form $a(\cdot,\cdot)$ defines an inner product on
$\mathbf{V}$, due to Assumption 5.1 on the domain. The inner product for $V$
is inherited from $H^{1}(\Omega).$

We have the following results on well-posedness of weak solutions [8, 9]:

######  Proposition 5.2 (Poro-elasticity).

Let $\mathbf{F}\in H^{1}(0,T;\mathbf{L}^{2}(\Omega))$, $S\in
L^{2}(0,T;L^{2}(\Omega))$, and $\mathbf{g}\in H^{1}(0,T;L^{2}(\Gamma_{N}))$.
Additionally, let $\psi(x,t)=w(t)\chi(x)$ where $w(t)\in L^{2}(0,T)$ and
$\chi\in L^{2}(\Gamma_{D,v})$ Then there exists unique weak solution
$(\mathbf{u},p)\in L^{2}(0,T;\mathbf{V})\times L^{2}(0,T;V)$ to (5.1-5.5).
Additionally, the solution satisfies the following energy estimate:

$\int_{0}^{T}\|p\|^{2}_{V}+\|\mathbf{u}\|^{2}_{\mathbf{V}}dt\leq
C(\|\mathbf{u}_{0}\|_{\mathbf{V}}+\|\mathbf{F}\|_{H^{1}(0,T;\mathbf{L}^{2}(\Omega))}^{2}+\|S\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\|\mathbf{g}\|_{H^{1}(0,T;\mathbf{L}^{2}(\Gamma_{N}))}^{2}+\|\psi\|_{L^{2}(0,T;L^{2}(\Gamma_{D,v}))}^{2})$
(5.7)

######  Proposition 5.3 (Poro-visco-elasticity).

Let $\mathbf{F}\in L^{2}(0,T;\mathbf{L}^{2}(\Omega))$, $S\in
L^{2}(0,T;L^{2}(\Omega))$, and
$\mathbf{g}\in L^{2}(0,T;L^{2}(\Gamma_{N}))$. Additionally, let
$\psi(x,t)=w(t)\chi(x)$ where $w(t)\in L^{2}(0,T)$ and $\chi\in
L^{2}(\Gamma_{D,v})$ Then there exists unique weak solution $(\mathbf{u},p)\in
H^{1}(0,T;\mathbf{V})\times L^{2}(0,T;V)$ to (5.1-5.5). Additionally, the
solution satisfies the following energy estimate:

$\begin{split}&\int_{0}^{T}\|p\|^{2}_{V}+\|\mathbf{u}\|^{2}_{\mathbf{V}}+\delta\|\mathbf{u}_{t}\|^{2}_{\mathbf{V}}dt\\\
&\leq
C(\|\mathbf{u}_{0}\|_{\mathbf{V}}+\|\mathbf{F}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}^{2}+\|S\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\|\mathbf{g}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Gamma_{N}))}^{2}+\|\psi\|_{L^{2}(0,T;L^{2}(\Gamma_{D,v}))}^{2})\end{split}$
(5.8)

Similar results hold when the control $w$ is used as the time portion of $S$.
Therefore, we will let $\chi(x)$ be fixed and insert $q=\chi(x)w(t)$ into
(5.1)-(5.5) in place of either $S$ or $\psi$. Let $G:L^{2}(0,T)\to
L^{2}(0,T;\mathbf{L}^{2}(\Omega))\times L^{2}(0,T;L^{2}(\Omega))$ be defined
by mapping $w(t)\in L^{2}(0,T)$ to the unique solution of (5.1)-(5.5),
embedded in $L^{2}(0,T;\mathbf{L}^{2}(\Omega))\times
L^{2}(0,T;L^{2}(\Omega))$, with all sources set to zero except $q$. Let
$(\tilde{\mathbf{u}},\tilde{p})$ be the unique solution (in the spaces
provided in Proposition 5.2 and Proposition 5.3) to (5.1)-(5.5) with the
sources set as desired and $q=0$. We introduce the control problem

$\displaystyle\begin{aligned} \min_{w\in L^{2}(0,T)}\
&\frac{1}{2}\|\mathbf{u}-\mathbf{u}_{d}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}^{2}+\frac{1}{2}\|p-p_{d}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\frac{\lambda}{2}\|w\|_{L^{2}(0,T)}^{2}+\alpha\operatorname*{TV}(w)\\\
\text{s.t.}\quad\quad&(\mathbf{u},p)=G(w)+(\tilde{\mathbf{u}},\tilde{p}),\\\
&w(t)\in W\text{ for a.a.\ }t\in(0,T)\end{aligned}$ (5.9)

for given $\mathbf{u}_{d}\in L^{2}(0,T;\mathbf{L}^{2}(\Omega))$, $p_{d}\in
L^{2}(0,T;L^{2}(\Omega))$, $\chi\in L^{2}(\Gamma_{D,v})$ and $\lambda\geq 0$.
We define

$\displaystyle
J(\mathbf{u},p,w)\coloneqq\frac{1}{2}\|\mathbf{u}-\mathbf{u}_{d}\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega))}^{2}+\frac{1}{2}\|p-p_{d}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\frac{\lambda}{2}\|w\|_{L^{2}(0,T)}^{2}$
(5.10)

for $(\mathbf{u},p,w)\in L^{2}(0,T;\mathbf{V})\times L^{2}(0,T;V)\times
L^{2}(0,T)$ and the reduced objective

$j:L^{2}(0,T)\to\mathbb{R},\ j(w)\coloneqq
J(G(w)+(\tilde{\mathbf{u}},\tilde{p}),w).$ (5.11)

Next we verify the applicability of the SLIP algorithm introduced in the
previous section, by checking if Assumption 4.3 and Assumption 4.1 are
satisfied.

###  5.1 Case 1: $\lambda>0$

In this subsection, we prove that in both poro-elastic and poro-visco elastic
cases Assumption 4.3 is satisfied while the more stringent Assumption 4.1 is
not satisfied.

######  Proposition 5.4.

Let $\xi\in L^{2}(0,T)$. Then it follows that the reduced objective
$j:L^{2}(0,T)\to\mathbb{R}$ is twice continuously differentiable at $\xi$.
Moreover, $|\nabla^{2}j(w)(\xi,\phi)|=|(G(\xi),G(\varphi))+\lambda(\xi,\phi)|$
and there exists $C>0$ such that

$|\nabla^{2}j(w)(\xi,\varphi)|\leq
C\|\xi\|_{L^{2}(0,T)}\|\varphi\|_{L^{2}(0,T)}$ (5.12)

for all $w\in L^{2}(0,T)$ and all $\xi$, $\varphi\in L^{2}(0,T)$ i.e. the
Hessian is continuous on $L^{2}(0,T)\times L^{2}(0,T)$. When $\lambda>0$,
Assumption 4.1 does not hold.

###### Proof.

Let $\tilde{y}=(\mathbf{u}_{d}-\tilde{\mathbf{u}},p_{d}-\tilde{p})$. Based on
the definition (5.11) of the reduced functional $j$, we have that

$\nabla j(q)=(G(\cdot),G(q)-\tilde{y})+\lambda(\cdot,q)\ \text{and}\
\nabla^{2}j(q)(\xi,\varphi)=(G(\varphi),G(\xi))+\lambda(\varphi,\xi)$

Furthermore, we can estimate the Hessian of $j$ as follows:

$\begin{split}&\left|\nabla^{2}j(q)(\xi,\varphi)\right|=|(G(\varphi),G(\xi))+\lambda(\varphi,\xi)|\\\
&\leq\|G(\varphi)\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega)\times
L^{2}(\Omega))}\|G(\xi)\|_{L^{2}(0,T;\mathbf{L}^{2}(\Omega)\times
L^{2}(\Omega))}+\lambda\|\varphi\|_{L^{2}(0,T)}\|\xi\|_{L^{2}(0,T)}\leq(C^{2}+\lambda)\|\varphi\|\|\xi\|,\end{split}$
(5.13)

which provides the desired estimate (5.12). Additionally, we note that

$\lambda\|\varphi\|^{2}_{L^{2}(\Omega)}\leq\lambda(\varphi,\varphi)+(G(\varphi),G(\varphi))=\nabla^{2}j(q)(\varphi,\varphi)\quad\text{for
all }\varphi\in L^{2}(0,T).$ (5.14)

Consider the sequence $\phi_{n}(t)=n\left(\frac{t}{T}\right)^{n}\in
L^{2}(0,T)$. Note that

$\|\phi_{n}\|^{2}_{L^{1}(0,T)}=\left(\int_{0}^{T}\left|n\left(\frac{t}{T}\right)^{n}\right|dt\right)^{2}=\frac{n^{2}T^{2}}{(n+1)^{2}}$
(5.15)
$\lambda\|\phi_{n}\|^{2}_{L^{2}(0,T)}=\lambda\int_{0}^{T}\left|n\left(\frac{t}{T}\right)^{n}\right|^{2}dt=\frac{\lambda
n^{2}T}{2n+1}.$ (5.16)

Combining (5.14) with (5.15) and (5.16) we obtain

$\frac{|\nabla^{2}j(w)(\phi_{n},\phi_{n})|}{\|\phi_{n}\|_{L^{1}(0,T)}^{2}}\geq\frac{\lambda\|\phi_{n}\|^{2}_{L^{2}(\Omega)}}{\|\phi_{n}\|_{L^{1}(\Omega)}^{2}}=\frac{\lambda
n^{2}T}{2n+1}\frac{(n+1)^{2}}{n^{2}T^{2}}=\lambda\frac{(n+1)^{2}}{(2n+1)T}\to\infty\quad\text{as
}n\to\infty.$

Hence, Assumption 4.1 does not hold. ∎

###  5.2 Case 2: $\lambda=0$

Now we consider the case $\lambda=0$. For the following proofs and numerical
results, we consider the problem on one spatial dimension i.e. $\Omega=(0,L).$
Following [47, 7], the one dimensional formulation of the poro-visco-elastic
systems is given by

$\displaystyle H_{v}\frac{\partial^{3}u}{\partial x^{2}\partial
t}+H_{e}\frac{\partial^{2}u}{\partial x^{2}}-\frac{\partial p}{\partial x}$
$\displaystyle=0$ $\displaystyle\forall(x,t)\in[0,L]\times[0,T]$ (5.17)
$\displaystyle\frac{\partial^{2}u}{\partial x\partial
t}-k\frac{\partial^{2}p}{\partial x^{2}}$ $\displaystyle=S$
$\displaystyle\forall(x,t)\in[0,L]\times[0,T]$ (5.18) $\displaystyle
u(0,t)=p(0,t)$ $\displaystyle=0$ $\displaystyle\forall t\in[0,T]$ (5.19)
$\displaystyle-k\frac{\partial p}{\partial x}(L,t)$ $\displaystyle=\psi(t)$
$\displaystyle\forall t\in[0,T]$ (5.20) $\displaystyle
H_{v}\frac{\partial^{2}u}{\partial x\partial t}(L,t)+H_{e}\frac{\partial
u}{\partial x}(L,t)-p(L,t)$ $\displaystyle=0$ $\displaystyle\forall t\in[0,T]$
(5.21) $\displaystyle\frac{\partial u}{\partial x}(x,0)$ $\displaystyle=0$
$\displaystyle\forall x\in[0,L]$ (5.22)

where $H_{e}=\lambda_{e}+2\mu_{e}$ and $H_{v}=\delta(\lambda_{v}+2\mu_{v})$.
When $F=0$, we know $p(x,t)=H_{e}\frac{\partial u}{\partial
x}+H_{v}\frac{\partial^{2}u}{\partial x\partial t}+h(t)$. When $g=0$, applying
(5.21) shows $h(t)$ must be equal to 0. Therefore, when $F=0$ and $g=0$ we
have

$p(x,t)=H_{e}\frac{\partial u}{\partial x}+H_{v}\frac{\partial^{2}u}{\partial
x\partial t}.$ (5.23)

When $\lambda=0$, we will see that Assumption 4.1 holds in the poro-elastic
case, but not the poro-visco-elastic case. We will study these two cases
separately.

####  5.2.1 Poroelastic Case

We will first consider the case when $\psi$ is used as the control. Then we
will study the case where the control $w(t)$ enters the system in the source
$S$, i.e. $S(x,t)=\chi(x)w(t)$.

######  Theorem 5.5.

Let $\lambda=0$ and $\delta=0$. Let $k(x,t)=k$ where $k$ is a positive
constant and let $G$ map a control $\psi$ to the state $(u,p)$ that satisfies
(5.17)-(5.22) with all sources set to zero except $\psi(t)$. Then Assumption
4.1 holds with $w=\psi$.

###### Proof.

Note that $\psi$ is only a function of $t$ because we are considering the one-
dimensional case where $\Gamma_{N}$ is only one point. Since $\delta=0$,
$H_{v}=0$. Plugging (5.23) into the PDE we see that $p(x,t)$ needs to satisfy

$\displaystyle\frac{\partial p}{\partial
t}-kH_{e}\frac{\partial^{2}p}{\partial x^{2}}$ $\displaystyle=0\quad$
$\displaystyle\forall(x,t)\in[0,L]\times[0,T]$ $\displaystyle
p(0,t)=0,\quad-k\frac{\partial p}{\partial x}(L,t)$
$\displaystyle=\psi(t)\quad$ $\displaystyle\forall t\in[0,T]$ $\displaystyle
p(x,0)$ $\displaystyle=0\quad$ $\displaystyle\forall x\in[0,L].$

Let $w=p+\frac{x}{k}\psi(t)$. Then we see

$\displaystyle\frac{\partial w}{\partial
t}-kH_{e}\frac{\partial^{2}w}{\partial x^{2}}$
$\displaystyle=\frac{x}{k}\psi^{\prime}(t)\quad$
$\displaystyle\forall(x,t)\in[0,L]\times[0,T]$ $\displaystyle
w(0,t)=0,\quad-k\frac{\partial w}{\partial x}(L,t)$ $\displaystyle=0\quad$
$\displaystyle\forall t\in[0,T]$ $\displaystyle w(x,0)$ $\displaystyle=0\quad$
$\displaystyle\forall x\in[0,L].$

Let $\lambda_{n}=\frac{(2n-1)\pi}{2L}$ Note that
$(\sqrt{2/L}\sin(\lambda_{n}x))_{n\in\mathbb{N}}$ is a complete orthonormal
basis. Let

$w(x,t)=\sum_{n=1}^{\infty}f_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x).$

Then we see $w(x,t)$ satisfies all boundary conditions, and we have

$\sum_{n=1}^{\infty}f^{\prime}_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x)+\lambda_{n}^{2}kH_{e}\sum_{n=1}^{\infty}f_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x)=\frac{x}{k}\psi^{\prime}(t)$

and

$\sum_{n=1}^{\infty}f_{n}(0)\sqrt{2/L}\sin(\lambda_{n}x)=0.$

Multiplying both sides of these equations by
$\sqrt{2/L}\sin(\lambda_{n}x)$,integrating these equations from $0$ to $L$,
and setting
$c_{n}=(\frac{x}{k},\sqrt{2/L}\sin(\lambda_{n}x))_{L^{2}(0,L)}=\frac{(-1)^{n}\sqrt{2L}}{k\lambda_{n}^{2}}$,
we have

$f_{n}^{\prime}(t)+kH_{e}\lambda_{n}^{2}f_{n}(t)=c_{n}\psi^{\prime}(t)\quad\text{and}\quad
f_{n}(0)=0.$

Therefore,
$f_{n}(t)=e^{-kH_{e}\lambda_{n}^{2}t}\int_{0}^{t}c_{n}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi$.
Hence, we have

$\displaystyle w(x,t)$
$\displaystyle=\sum_{n=1}^{\infty}e^{-kH_{e}\lambda_{n}^{2}t}\int_{0}^{t}c_{n}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\sqrt{2/L}\sin(\lambda_{n}x),$
$\displaystyle p(x,t)$
$\displaystyle=\sum_{n=1}^{\infty}e^{-kH_{e}\lambda_{n}^{2}t}\int_{0}^{t}c_{n}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\sqrt{2/L}\sin(\lambda_{n}x)-\frac{x}{k}\psi(t).$

Using (5.23) and recalling that $H_{v}=0$, we anti-differentiate $p$ with
respect to $x$ and enforcing the boundary condition $u(x,0)=0$, to see

$u(x,t)=\sum_{n=1}^{\infty}\frac{e^{-kH_{e}\lambda_{n}^{2}t}}{H_{e}\lambda_{n}}\int_{0}^{t}c_{n}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\sqrt{2/L}(1-\cos(\lambda_{n}x))-\frac{x^{2}}{2kH_{e}}\psi(t).$

Using the fact that
$c_{n}=(\frac{x}{k},\sqrt{2/L}\sin(\lambda_{n}x))_{L^{2}(\Omega)}$,
$(\sqrt{2/L}\sin(\lambda_{n}x))_{n\in\mathbb{N}}$ is an orthonormal sequence,
Lebesgue Dominated Convergence Theorem, and Parseval’s equality we have

$\displaystyle\|p\|^{2}_{L^{2}(0,T;L^{2}(0,L))}$
$\displaystyle=\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}\left[\left(e^{-kH_{e}\lambda_{n}^{2}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\right)^{2}-2e^{-kH_{e}\lambda_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\psi(t)\right]+\frac{x^{2}}{k^{2}}\psi^{2}(t)dt$
$\displaystyle=\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}\left[\left(e^{-kH_{e}\lambda_{n}^{2}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\right)^{2}-2e^{-kH_{e}\lambda_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi}\psi^{\prime}(\xi)d\xi\psi(t)+\psi^{2}(t)\right]dt$

Integrating by parts, using $\psi(0)=0$ we see

$\displaystyle\|p\|^{2}_{L^{2}(0,T;L^{2}(0,L))}=\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}\left[\left(e^{-kH_{e}\lambda_{n}^{2}t}\left(\psi(t)e^{kH_{e}\lambda_{n}^{2}t}-\int_{0}^{t}\psi(\xi)kH_{e}\lambda_{n}^{2}e^{kH_{e}\lambda_{n}^{2}\xi}d\xi\right)\right)^{2}\right.$
$\displaystyle\left.-2e^{-kH_{e}\lambda_{n}^{2}t}\left(\psi(t)e^{kH_{e}\lambda_{n}^{2}t}-\int_{0}^{t}\psi(\xi)kH_{e}\lambda_{n}^{2}e^{kH_{e}\lambda_{n}^{2}\xi}d\xi\right)\psi(t)+\psi^{2}(t)\right]dt$
$\displaystyle=\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}\left[\psi^{2}(t)-2\psi(t)\int_{0}^{t}\psi(\xi)kH_{e}\lambda_{n}^{2}e^{kH_{e}\lambda_{n}^{2}(\xi-t)}d\xi+\left(\int_{0}^{t}\psi(\xi)kH_{e}\lambda_{n}^{2}e^{kH_{e}\lambda_{n}^{2}(\xi-t)}d\xi\right)^{2}-2\psi^{2}(t)\right.$
$\displaystyle\left.+2\psi(t)\int_{0}^{t}\psi(\xi)kH_{e}\lambda_{n}^{2}e^{kH_{e}\lambda_{n}^{2}(\xi-t)}d\xi+\psi^{2}(t)\right]dt=\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}\left(\int_{0}^{t}\psi(\xi)kH_{e}\lambda_{n}^{2}e^{kH_{e}\lambda_{n}^{2}(\xi-t)}d\xi\right)^{2}dt$
$\displaystyle\leq\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}k^{2}H_{e}^{2}\lambda_{n}^{4}\left(\int_{0}^{T}\psi(\xi)e^{-kH_{e}\lambda_{n}^{2}(t-\xi)}d\xi\right)^{2}dt\leq\int_{0}^{T}\sum_{n=1}^{\infty}H_{e}^{2}\left(\int_{0}^{T}\psi(\xi)e^{-kH_{e}\lambda_{n}^{2}(t-\xi)}d\xi\right)^{2}dt$

Recall, Young’s inequality for convolution: Given $f\in L^{p}(0,T)$ and $g\in
L^{q}(0,T)$ such that $\frac{1}{p}+\frac{1}{q}=1+\frac{1}{r}$ with $1\leq
p,q,r\leq\infty$, $\|f*g\|_{L^{r}(0,T)}\leq
C\|f\|_{L^{p}(0,T)}\|g\|_{L^{q}(0,T)}$. Using $p=1$, $q=2$, and $r=2$, we have

$\displaystyle\|p\|^{2}_{L^{2}(0,T;L^{2}(\Omega))}$
$\displaystyle\leq\sum_{n=1}^{\infty}H_{e}^{2}\|e^{-kH_{e}\lambda_{n}^{2}(\cdot)}\|_{L^{2}(0,T)}^{2}\|\psi\|_{L^{1}(0,T)}^{2}\leq\sum_{n=1}^{\infty}\frac{H_{e}^{2}}{2kH_{e}\lambda_{n}^{2}}\|\psi(\xi)\|^{2}_{L^{1}(0,T)}$
$\displaystyle=\frac{H_{e}}{2k}\|\psi\|_{L^{1}(0,T)}^{2}\sum_{n=1}^{\infty}\frac{4L^{2}}{(2n-1)^{2}\pi^{2}}=\frac{H_{e}L^{2}}{4k}\|\psi\|_{L^{1}(0,T)}^{2}$

Additionally, integrating (5.23) with respect to $x$, we have

$\displaystyle\|u\|_{L^{2}(0,T;L^{2}(0,L))}^{2}$
$\displaystyle=\int_{0}^{T}\int_{0}^{L}\left(\int_{0}^{x}p(\xi,t)d\xi\right)^{2}dxdt\leq\frac{H_{e}L^{3}}{4k}\|\psi\|_{L^{1}(0,T)}$

Therefore, for all $\psi\in C_{0}^{1}(0,T)$, we have
$\|G(\psi)\|^{2}_{L^{2}(0,T;L^{2}(0,L)\times L^{2}(0,L))}\leq
C\|\psi\|_{L^{1}(0,T)}$. Hence, $G|_{C_{0}^{1}(0,T)}$ is bounded and linear.
Using the Bounded Linear Extension theorem, we have that the extension
$G:L^{1}(0,T)\to L^{2}(0,T;L^{2}(0,L)\times L^{2}(0,L))$ is a bounded linear
functional. Hence, the desired bound holds. ∎

Now we consider the poro-elastic case with $\lambda=0$ where the control is
used in the source $S$. Let $\chi(x)$ be set and consider the case where $G$
maps $s\in L^{2}(0,T)$ to the unique solution $(u,p)$ to (5.1)-(5.5) with all
sources and the initial condition set to zero except $S:=\chi(x)s(t)$. In this
case, the process of finding $p(x,t)$ is similar to the process of finding
$w(x,t)$ in the proof of Theorem 5.5. Hence, $G(s)=(u,p)$ where

$p(x,t)=\sum_{n=1}^{\infty}e^{-kH_{e}\lambda_{n}^{2}t}\int_{0}^{t}d_{n}e^{kH_{e}\lambda_{n}^{2}\xi}s(\xi)d\xi\sqrt{2/L}\sin(\lambda_{n}x),$

and applying (5.23)

$u(x,t)=\sum_{n=1}^{\infty}\frac{e^{-kH_{e}\lambda_{n}^{2}t}}{H_{e}\lambda_{n}}\int_{0}^{t}d_{n}e^{kH_{e}\lambda_{n}^{2}\xi}s(\xi)d\xi\sqrt{2/L}\sin(\lambda_{n}x)$

where $d_{n}=(\chi(x),\sqrt{2/L}\sin(\lambda_{n}x))_{L^{2}(0,L)}.$ The proof
for showing this solution also satisfies Assumption 4.1 follows similarly to
the proof of Theorem 5.5.

####  5.2.2 Poro-visco-elastic Case

We show that Assumption 4.1 is not satisfied in the poro-visco-elastic case,
for both choices of controls $\psi$ and $S$.

######  Theorem 5.6.

Let $\lambda=0$ and $\delta>0$. Let $k(x,t)=k$ where $k$ is a positive
constant and let $G$ map a control $\psi$ to the state $(u,p)$ that satisfies
(5.17)-(5.22) with all sources set to zero except $\psi(t)$. Then Assumption
4.1 does not hold.

###### Proof.

Recall that $\nabla^{2}j(w)(\xi,\varphi)=(G(\xi),G(\varphi))$. Hence, we will
proceed with this proof by first finding $G(\psi)=(u,p)$ when $\psi(t)\in
C_{0}^{1}(0,T)$. We want to show that there exists
$(\varphi_{m})_{m\in\mathbb{N}}$ such that

$\frac{(G(\varphi_{m}),G(\varphi_{m}))}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}=\frac{\|p_{m}\|^{2}_{L^{2}(0,T)}+\|u_{m}\|^{2}_{L^{2}(0,T)}}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}\to\infty~{}~{}\text{as
}m\to\infty.$

Plugging (5.23) into (5.17)-(5.22), we see that $u(x,t)$ needs to satisfy

$\displaystyle\frac{\partial^{2}u}{\partial x\partial
t}-kH_{v}\frac{\partial^{4}u}{\partial x^{3}\partial
t}-kH_{e}\frac{\partial^{3}u}{\partial x^{3}}$ $\displaystyle=0\quad$
$\displaystyle\forall(x,t)\in[0,L)\times[0,T]$ $\displaystyle u(0,t)$
$\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$ $\displaystyle
H_{v}\frac{\partial^{2}u}{\partial x\partial t}(0,t)+H_{e}\frac{\partial
u}{\partial x}(0,t)$ $\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$
$\displaystyle-kH_{v}\frac{\partial^{3}u}{\partial x^{2}\partial
t}(L,t)-kH_{e}\frac{\partial^{2}u}{\partial x^{2}}(L,t)$
$\displaystyle=\psi(t)\quad$ $\displaystyle\forall t\in[0,T]$
$\displaystyle\frac{\partial u}{\partial x}(x,0)$ $\displaystyle=0\quad$
$\displaystyle\forall x\in[0,L]$

Let $y(x,t)=u(x,t)+\frac{x^{2}}{2L}\Psi(t)$ where

$\Psi(t)=\frac{1}{kH_{v}}e^{-H_{e}/H_{v}t}\int_{0}^{t}\psi(\xi)e^{H_{e}/H_{v}\xi}d\xi.$

Then $y(x,t)$ satisfies

$\displaystyle\frac{\partial^{2}y}{\partial x\partial
t}-kH_{v}\frac{\partial^{4}y}{\partial x^{3}\partial
t}-kH_{e}\frac{\partial^{3}y}{\partial x^{3}}$
$\displaystyle=\frac{x}{L}\Psi^{\prime}(t)\quad$
$\displaystyle\forall(x,t)\in[0,L]\times[0,T]$ $\displaystyle y(0,t)$
$\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$ $\displaystyle
H_{v}\frac{\partial^{2}y}{\partial x\partial t}(0,t)+H_{e}\frac{\partial
y}{\partial x}(0,t)$ $\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$
$\displaystyle-kH_{v}\frac{\partial^{3}y}{\partial x^{2}\partial
t}(L,t)-kH_{e}\frac{\partial^{2}y}{\partial x^{2}}(L,t)$
$\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$
$\displaystyle\frac{\partial y}{\partial x}(x,0)$ $\displaystyle=0\quad$
$\displaystyle\forall x\in[0,L]$

Let $y(x,t)=\sum_{n=1}^{\infty}f_{n}(t)\sqrt{2/L}(1-\cos(\lambda_{n}x))$ where
$\lambda_{n}=\frac{(2n-1)\pi}{2L}$. Notice that $y(0,t)=0$, $\frac{\partial
y}{\partial x}(0,t)=0$, and $\frac{\partial^{2}y}{\partial x^{2}}(L,t)=0$. By
plugging $y(x,t)$ into the first line of the PDE, we see that

$\displaystyle\sum_{n=1}^{\infty}f^{\prime}_{n}(t)\lambda_{n}\sqrt{2/L}\sin(\lambda_{n}x)+\lambda_{n}^{3}kH_{v}\sum_{n=1}^{\infty}f^{\prime}_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x)+\lambda_{n}^{3}kH_{e}\sum_{n=1}^{\infty}f_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x)=\frac{x}{L}\Psi^{\prime}(t)$

and

$\sum_{n=1}^{\infty}f_{n}(0)\sqrt{2/L}\sin(\lambda_{n}x)=0$

Note that $\left(\sqrt{2/L}\sin(\lambda_{n}x)\right)_{n\in\mathbb{N}}$ is a
complete orthonormal basis. Hence, multiplying both sides of these equations
by $\sqrt{2/L}\sin(\lambda_{n}x)$, integrating these equations from $0$ to
$L$, setting
$c_{n}=(\frac{x}{L},\sqrt{2/L}\sin(\lambda_{n}x))_{L^{2}(0,L)}=4\sqrt{2L}(-1)^{n+1}/((2n-1)^{2}\pi^{2})$
we have

$(\lambda_{n}+\lambda_{n}^{3}kH_{v})f_{n}^{\prime}(t)+kH_{e}\lambda_{n}^{3}f_{n}(t)=c_{n}\Psi^{\prime}(t)\quad\text{and}\quad
f_{n}(0)=0.$

Therefore,
$f_{n}(t)=\frac{1}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{(kH_{e}\lambda_{n}^{2}\gamma_{n}\xi)}c_{n}\Psi^{\prime}(\xi)d\xi$
where $\gamma_{n}=\frac{1}{1+\lambda_{n}^{2}kH_{v}}$. Hence,

$\displaystyle y(x,t)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}c_{n}\Psi^{\prime}(\xi)d\xi\sqrt{2/L}(1-\cos(\lambda_{n}x))$
$\displaystyle u(x,t)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{\gamma_{n}}{\lambda_{n}}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}c_{n}\Psi^{\prime}(\xi)d\xi\sqrt{2/L}(1-\cos(\lambda_{n}x))-\frac{x^{2}}{2L}\Psi(t)$

and from (5.23),

$\displaystyle
p(x,t)=-H_{e}\frac{x}{L}\Psi(t)-H_{v}\frac{x}{L}\Psi^{\prime}(t)+\sum_{n=1}^{\infty}c_{n}\lambda_{n}\sqrt{2/L}\sin(\lambda_{n}x)$
$\displaystyle\left(-kH_{e}H_{v}\lambda_{n}\gamma_{n}^{2}e^{-kH_{e}\gamma_{n}\lambda_{n}^{2}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi+\frac{H_{v}\gamma_{n}}{\lambda_{n}}\Psi^{\prime}(t)+\frac{H_{e}\gamma_{n}}{\lambda_{n}}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi\right)$
$\displaystyle=-\frac{x}{L}(H_{e}\Psi(t)+H_{v}\Psi^{\prime}(t))+\sum_{n=1}^{\infty}c_{n}\sqrt{2/L}\lambda_{n}\sin(\lambda_{n}x)\left(\frac{H_{e}\gamma_{n}^{2}}{\lambda_{n}}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi+\frac{H_{v}\gamma_{n}}{\lambda_{n}}\Psi^{\prime}(t)\right).$

We notice

$\Psi^{\prime}(t)=\frac{1}{kH_{v}}\left(\frac{-H_{e}}{H_{v}}e^{-H_{e}/H_{v}t}\int_{0}^{t}\psi(\xi)e^{H_{e}/H_{v}\xi}d\xi+\psi(t)\right).$
(5.24)

Hence,

$\displaystyle p(x,t)$
$\displaystyle=-\frac{x}{Lk}\psi(t)+\sum_{n=1}^{\infty}c_{n}\sqrt{2/L}\sin(\lambda_{n}x)\left(H_{e}\gamma_{n}^{2}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi+H_{v}\gamma_{n}\Psi^{\prime}(t)\right).$

Therefore,

$\displaystyle\|p\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}=\int_{0}^{T}\sum_{n=1}^{\infty}c_{n}^{2}\left(H_{e}\gamma_{n}^{2}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi+H_{v}\gamma_{n}\Psi^{\prime}(t)\right)^{2}$
$\displaystyle-2\frac{c_{n}^{2}}{k}\psi(t)\left(H_{e}\gamma_{n}^{2}e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi+H_{v}\gamma_{n}\Psi^{\prime}(t)\right)+\frac{c_{n}^{2}}{k^{2}}\psi^{2}(t)dt.$

Recalling (5.24), we see that when $\psi=\varphi_{m}=e^{mt}$, we have

$\begin{split}\Psi^{\prime}(t)&=\frac{1}{kH_{v}}\left(-\frac{H_{e}}{H_{v}}e^{-H_{e}/H_{v}t}\int_{0}^{t}e^{m\xi}e^{H_{e}/H_{v}\xi}d\xi+e^{mt}\right)\\\
&=\frac{1}{kH_{v}}\left(-\frac{H_{e}}{H_{v}}\left(\frac{e^{mt}}{m+H_{e}/H_{v}}-\frac{e^{-H_{e}/H_{v}t}}{m+H_{e}/H_{v}}\right)+e^{mt}\right)\\\
&=\frac{1}{kH_{v}}\left(\frac{H_{e}e^{-H_{e}/H_{v}t}}{mH_{v}+H_{e}}-\frac{H_{e}e^{mt}}{mH_{v}+H_{e}}+e^{mt}\right)\\\
&=\frac{1}{kH_{v}}\left(\frac{H_{e}e^{-H_{e}/H_{v}t}}{mH_{v}+H_{e}}+\frac{H_{v}me^{mt}}{H_{v}m+H_{e}}\right)\geq
0\end{split}$ (5.25)

Hence,

$\displaystyle
e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi$
$\displaystyle=\frac{1}{kH_{v}(mH_{v}+H_{e})}\left(\frac{H_{e}}{kH_{e}\lambda_{n}^{2}\gamma_{n}-H_{e}/H_{v}}\left(e^{-H_{e}/H_{v}t}-e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\right)+\frac{H_{v}m}{kH_{e}\lambda_{n}^{2}\gamma_{n}+m}\left(e^{mt}-e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\right)\right)$
$\displaystyle=\frac{1}{kH_{v}(mH_{v}+H_{e})}\left(\frac{H_{v}}{k\lambda_{n}^{2}\gamma_{n}H_{v}-1}\left(e^{-H_{e}/H_{v}t}-e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\right)+\frac{H_{v}m}{kH_{e}\lambda_{n}^{2}\gamma_{n}+m}\left(e^{mt}-e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\right)\right)$

Recall $\gamma_{n}=\frac{1}{1+\lambda_{n}^{2}kH_{v}}$, so

$\displaystyle k\lambda_{n}^{2}\gamma_{n}H_{v}-1$
$\displaystyle=\frac{k\lambda_{n}^{2}H_{v}}{1+\lambda_{n}^{2}kH_{v}}-\frac{1+\lambda_{n}^{2}kH_{v}}{1+\lambda_{n}^{2}kH_{v}}$
$\displaystyle=\frac{-1}{1+\lambda_{n}^{2}kH_{v}}=-\gamma_{n}$

and

$-kH_{e}\lambda_{n}^{2}\gamma_{n}=-\frac{kH_{e}\lambda_{n}}{1+\lambda_{n}^{2}kH_{v}}>-\frac{kH_{e}\lambda_{n}^{2}}{\lambda_{n}^{2}kH_{v}}=-\frac{H_{e}}{H_{v}},$

which implies $1\geq
e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}-e^{-H_{e}/H_{v}t}\geq 0$ for
$t\in[0,T]$. Hence,

$e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\gamma_{n}\xi}\Psi^{\prime}(\xi)d\xi\geq\frac{1}{kH_{v}(mH_{v}+H_{e})}\left(\frac{H_{v}}{\gamma_{n}}+\frac{H_{v}m}{kH_{e}\lambda_{n}^{2}\gamma_{n}+m}\left(e^{mt}-e^{-kH_{e}\lambda_{n}^{2}\gamma_{n}t}\right)\right)$
(5.26)

Let $l>0$ satisfy

$\lambda_{l}^{2}=\frac{(2l-1)^{2}\pi^{2}}{4L^{2}}>\frac{8}{kH_{v}}.~{}~{}\text{Then
}\gamma_{l}=\frac{1}{1+\lambda_{l}^{2}kH_{v}}<\frac{1}{\lambda_{l}^{2}kH_{v}}<\frac{k}{8}.$
(5.27)

Notice that

$\begin{split}\|p_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}&\geq\int_{0}^{T}c_{l}^{2}\left(H_{e}\gamma_{l}^{2}e^{-kH_{e}\lambda_{l}^{2}\gamma_{l}t}\int_{0}^{t}e^{kH_{e}\lambda_{l}^{2}\gamma_{l}\xi}\Psi^{\prime}(\xi)d\xi+H_{v}\gamma_{l}\Psi^{\prime}(t)\right)^{2}\\\
&-2\frac{c_{l}^{2}}{k}\psi(t)\left(H_{e}\gamma_{l}^{2}e^{-kH_{e}\lambda_{l}^{2}\gamma_{l}t}\int_{0}^{t}e^{kH_{e}\lambda_{l}^{2}\gamma_{l}\xi}\Psi^{\prime}(\xi)d\xi+H_{v}\gamma_{l}\Psi^{\prime}(t)\right)+\frac{c_{l}^{2}}{k^{2}}\psi^{2}(t)dt\end{split}$

Dropping the first term since it is positive and applying (5.26) and (5.25) we
see

$\begin{split}\|p_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}&\geq\int_{0}^{T}-2\frac{c_{l}^{2}}{k}e^{mt}\frac{H_{e}\gamma_{l}^{2}}{kH_{v}(mH_{v}+H_{e})}\left(\frac{H_{v}}{\gamma_{l}}+\frac{H_{v}m}{kH_{e}\lambda_{l}^{2}\gamma_{l}+m}\left(e^{mt}-e^{-kH_{e}\lambda_{l}^{2}\gamma_{l}t}\right)\right)\\\
&-2\frac{c_{l}^{2}}{k}e^{mt}\frac{\gamma_{l}}{k}\left(\frac{H_{e}e^{-H_{e}/H_{v}t}}{mH_{v}+H_{e}}+\frac{H_{v}me^{mt}}{H_{v}m+H_{e}}\right)+\frac{c_{l}^{2}}{k^{2}}e^{2mt}dt\end{split}$

Recalling that $-(e^{mt}-e^{-kH_{e}\lambda_{l}^{2}\gamma_{l}t})\geq-e^{mt}$,
$-\left(\frac{H_{e}e^{-H_{e}/H_{v}t}}{mH_{v}+H+e}+\frac{H_{v}me^{mt}}{H_{v}m+H_{e}}\right)\geq-(1+e^{mt})$,
and $mH_{v}+H_{e}\geq H_{e}$ we see

$\begin{split}\|p_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}&\geq\int_{0}^{T}-2\frac{c_{l}^{2}}{k}e^{mt}\frac{H_{e}\gamma_{l}^{2}}{kH_{v}(mH_{v}+H_{e})}\left(\frac{H_{v}}{\gamma_{l}}+\frac{H_{v}m}{m}e^{mt}\right)-2\frac{c_{l}^{2}}{k}e^{mt}\frac{\gamma_{l}}{k}(1+e^{mt})+\frac{c_{l}^{2}}{k^{2}}e^{2mt}dt\\\
&\geq\int_{0}^{T}-2\frac{c_{l}^{2}}{k}e^{mt}\left(\frac{H_{e}\gamma_{l}^{2}}{kH_{v}H_{e}}\frac{H_{v}}{\gamma_{l}}+\frac{H_{e}\gamma_{l}^{2}H_{v}}{kH_{v}H_{e}}e^{mt}\right)-2\frac{c_{l}^{2}}{k}e^{mt}\frac{\gamma_{l}}{k}(1+e^{mt})+\frac{c_{l}^{2}}{k^{2}}e^{2mt}dt\\\
&\geq\int_{0}^{T}-2\frac{c_{l}^{2}}{k^{2}}e^{2mt}\left(\gamma_{l}+\gamma_{l}^{2}\right)-2\frac{c_{l}^{2}}{k^{2}}e^{2mt}2\gamma_{l}+\frac{c_{l}^{2}}{k^{2}}e^{2mt}dt\end{split}$

Hence, applying (5.27), we see

$\|p_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}\geq\int_{0}^{T}-2\frac{c_{l}^{2}}{k^{2}}e^{2mt}\left(\frac{1}{8}+\frac{1}{64}+\frac{1}{4}\right)+\frac{c_{l}^{2}}{k^{2}}e^{2mt}dt\geq\int_{0}^{T}\frac{7c_{l}^{2}}{32^{2}}e^{2mt}dt=\frac{7c_{l}^{2}}{64k^{2}m}\left(e^{2mT}-1\right).$
(5.28)

Also, note $\|\varphi_{m}\|_{L^{1}(0,T)}^{2}=\frac{(e^{mt}-1)^{2}}{m^{2}}$.
Therefore,

$\displaystyle\frac{(G(\varphi_{m}),G(\varphi_{m}))}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}$
$\displaystyle=\frac{\|p_{m}\|^{2}_{L^{2}(0,T;L^{2}(0,L))}+\|u_{m}\|^{2}_{L^{2}(0,T;L^{2}(0,L))}}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}\geq\frac{\|p_{m}\|^{2}_{L^{2}(0,T;L^{2}(0,L))}}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}$
$\displaystyle\geq\frac{\frac{7c_{l}}{64mk^{2}}(e^{2mT}-1)m^{2}}{(e^{mT}-1)^{2}}\to\infty~{}~{}\text{as
}m\to\infty.$

Therefore, Assumption 4.1 is not satisfied. ∎

We now consider the poro-visco-elastic case where the control is used as the
time component of the source $S$, and show Assumption 4.1 is not satisfied.

######  Theorem 5.7.

Let $\lambda=0$ and $\delta>0$. Let $k(x,t)=k$ and $\chi(x)\in L^{2}(0,L)$.
Define $G:L^{2}(0,T)\to L^{2}(0,T;L^{2}(0,L)\times L^{2}(0,L))$ to be the map
that maps $s(t)\in L^{2}(0,T)$ to the $(u,p)$ that satisfies (5.17)-(5.22)
with all sources set to 0 except $S(x,t)=s(t)\chi(x)$. Then Assumption 4.1
does not hold.

###### Proof.

We will proceed with this proof by first computing $G(s)=(u,p)$ and then
showing a lower estimate on $\|p\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}$. Plugging
(5.23) into the PDE we see $u(x,t)$ needs to satisfy

$\displaystyle\frac{\partial^{2}u}{\partial x\partial
t}-kH_{v}\frac{\partial^{4}u}{\partial x^{3}\partial
t}-kH_{e}\frac{\partial^{3}u}{\partial x^{3}}$ $\displaystyle=S\quad$
$\displaystyle\forall(x,t)\in[0,L)\times[0,T]$ $\displaystyle u(0,t)$
$\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$ $\displaystyle
H_{v}\frac{\partial^{2}u}{\partial x\partial t}(0,t)+H_{e}\frac{\partial
u}{\partial x}(0,t)$ $\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$
$\displaystyle-kH_{v}\frac{\partial^{3}u}{\partial x^{2}\partial
t}(L,t)-kH_{e}\frac{\partial^{2}u}{\partial x^{2}}(L,t)$
$\displaystyle=0\quad$ $\displaystyle\forall t\in[0,T]$
$\displaystyle\frac{\partial u}{\partial x}(x,0)$ $\displaystyle=0\quad$
$\displaystyle\forall x\in[0,L]$

Let $u(x,t)=\sum_{n=1}^{\infty}f_{n}(t)\sqrt{2/L}(1-\cos(\lambda_{n}x))$,
where $\lambda_{n}=\frac{(2n-1)\pi}{2L}$. Notice that $u(0,t)=0$,
$\frac{\partial u}{\partial x}(0,t)=0$, and $\frac{\partial^{2}u}{\partial
x^{2}}(L,t)=0$, so the boundary conditions are satisfied. We have

$\sum_{n=1}^{\infty}f^{\prime}_{n}(t)\lambda_{n}\sqrt{2/L}\sin(\lambda_{n}x)+\lambda_{n}^{3}kH_{v}\sum_{n=1}^{\infty}f^{\prime}_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x)+\lambda_{n}^{3}kH_{e}\sum_{n=1}^{\infty}f_{n}(t)\sqrt{2/L}\sin(\lambda_{n}x)=s(t)\chi(x)$

and

$\sum_{n=1}^{\infty}f_{n}(0)\sqrt{2/L}\sin(\lambda_{n}x)=0.$

Note $\left(\sqrt{2/L}\sin(\lambda_{n}x)\right)_{n\in\mathbb{N}}$ is a
complete orthonormal basis. Hence, multiplying both sides of these equations
by $\sqrt{2/L}\sin(\lambda_{n}x)$, integrating these equations from $0$ to
$L$, and setting $c_{n}=(\chi(x),\sqrt{2/L}\sin(\lambda_{n}x))_{L^{2}(0,L)}$,
we have for all $n\in\mathbb{N}$

$(\lambda_{n}+\lambda_{n}^{3}kH_{v})f_{n}^{\prime}(t)+kH_{e}\lambda_{n}^{3}f_{n}(t)=c_{n}s(t)\quad\text{and}\quad
f_{n}(0)=0.$

Therefore,
$f_{n}(t)=\frac{1}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}e^{-kH_{e}\lambda_{n}^{2}t/(1+\lambda_{n}^{2}kH_{v})}\int_{0}^{t}c_{n}e^{kH_{e}\lambda_{n}^{2}\xi/(1+\lambda_{n}^{2}kH_{v})}s(\xi)d\xi$.
Hence,

$u(x,t)=\sum_{n=1}^{\infty}\frac{c_{n}}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}e^{-kH_{e}\lambda_{n}^{2}t/(1+\lambda_{n}^{2}kH_{v})}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi/(1+\lambda_{n}^{2}kH_{v})}s(\xi)d\xi\sqrt{2/L}(1-\cos(\lambda_{n}x))$

and (5.23) gives

$\displaystyle p(x,t)$
$\displaystyle=\sum_{n=1}^{\infty}c_{n}\sqrt{2/L}\sin(\lambda_{n}x)\left(\frac{-kH_{e}H_{v}\lambda_{n}}{(1+\lambda_{n}^{2}kH_{v})^{2}}e^{-kH_{e}\lambda_{n}^{2}t/(1+\lambda_{n}^{2}kH_{v})}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi/(1+\lambda_{n}^{2}kH_{v})}s(\xi)d\xi\right.$
$\displaystyle\left.+\frac{H_{v}}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}s(t)+\frac{H_{e}}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}e^{-kH_{e}\lambda_{n}^{2}t/(1+\lambda_{n}^{2}kH_{v})}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi/(1+\lambda_{n}^{2}kH_{v})}s(\xi)d\xi\right)$
$\displaystyle=\sum_{n=1}^{\infty}c_{n}\sqrt{2/L}\sin(\lambda_{n}x)\left(\left(\frac{-kH_{e}H_{v}\lambda_{n}^{2}}{\lambda_{n}(1+\lambda_{n}^{2}kH_{v})^{2}}+\frac{H_{e}(1+\lambda_{n}^{2}kH_{v})}{\lambda_{n}(1+\lambda_{n}^{2}kH_{v})^{2}}\right)e^{-kH_{e}\lambda_{n}^{2}t/(1+\lambda_{n}^{2}kH_{v})}\right.$
$\displaystyle\left.\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi/(1+\lambda_{n}^{2}kH_{v})}s(\xi)d\xi+\frac{H_{v}}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}s(t)\right)$
$\displaystyle=\sum_{n=1}^{\infty}c_{n}\sqrt{2/L}\sin(\lambda_{n}x)\left(\frac{H_{e}}{\lambda_{n}(1+\lambda_{n}^{2}kH_{v})^{2}}e^{-kH_{e}\lambda_{n}^{2}t/(1+\lambda_{n}^{2}kH_{v})}\int_{0}^{t}e^{kH_{e}\lambda_{n}^{2}\xi/(1+\lambda_{n}^{2}kH_{v})}s(\xi)d\xi\right.$
$\displaystyle\left.+\frac{H_{v}}{\lambda_{n}+\lambda_{n}^{3}kH_{v}}s(t)\right)$

Therefore, when $s(t)$ is strictly non-negative,

$\|p\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}\geq\sum_{n=1}^{\infty}\frac{c_{n}^{2}H_{v}^{2}}{(\lambda_{n}+\lambda_{n}^{3}kH_{v})^{2}}\|s\|_{L^{2}(0,T)}^{2}\geq
C\|s\|^{2}_{L^{2}(0,T)}$

for some $C>0$. Therefore, using $\varphi_{m}=m\left(\frac{t}{T}\right)^{m}$
(as was done in the proof of Proposition 5.4), we have

$\displaystyle\frac{(G(\varphi_{m}),G(\varphi_{m})}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}\geq\frac{\|p_{m}\|_{L^{2}(0,T;L^{2}(\Omega))}}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}\geq\frac{C\|\varphi_{m}\|_{L^{2}(0,T)}}{\|\varphi_{m}\|_{L^{1}(0,T)}^{2}}\to\infty~{}~{}\text{as
}m\to\infty.$

Therefore, Assumption 4.1 is not satisfied. ∎

##  6 Computational Experiments

For our computational experiments, we use an instance of the porous medium
equations described in section 5. We intend to analyze the effect of the
mollification regularization on then resulting instationarity and the
objective values in practice. Specifically, we consider the control input
$\psi$, where Assumption 4.1 is satisfied for the poro-elastic case and
violated for the poro-visco-elastic case. The specific model as well as the
discretization of the model and the trust-region subproblems are described in
section 6.1. The computational experiments and the homotopy are described in
section 6.2. We present and describe our results in section 6.3.

###  6.1 Setup

We consider a discretized instance of (P) that is governed by the PDE
introduced in section 5 with a one-dimensional spatial domain $\Omega=(0,2)$
and a one-dimensional time domain $(0,T)=(0,0.5)$.

Regarding the parameters of the PDE we use $k=1$, $\lambda_{e}=1$,
$\mu_{e}=1$, $\mu_{v}=0.25$, $\lambda_{v}=0.0774$. We execute the same
experiments for both the two choices $\delta=0$ (poro-elastic case) and
$\delta=1$ (poro-visco-elastic case).

Regarding the setup of the control problem, we choose the structure given in
(5.9). We choose the set of possible control realizations as
$W=\\{-7,-5,-3,-1,0,2\\}$ in our computations. We choose the penalty parameter
value $\alpha=5\cdot 10^{-5}$ to scale the $\operatorname*{TV}$-term in the
objective. We run all experiments for the choices $\lambda=10^{-4}$,
$\lambda=10^{-2}$, and $\lambda=0$. We note that Assumption 4.1 is not
satisfied for both $\delta=0$ and $\delta=1$ for $\lambda=10^{-2}$ and
$\lambda=10^{-4}$ by virtue of Proposition 5.4 regardless of the fact which of
the control inputs is chosen. For the tracking terms, we choose
$u_{d}(x,t)=0.5+(1-t)^{2}\cos(50t)(-1.975x+4)$ and
$p_{d}(x,t)=0.5+\cos(50t)^{2}$ for all $(x,t)\in\Omega\times(0,T)$.

We use a one-dimensional spatial domain that is discretized as in [8] into
$512$ intervals. We discretize the PDE following the [8]. We choose $T=0.5$
and discretize the time horizon uniformly into $N=512$ intervals. $N=512$
intervals are used for both the implicit Euler time-stepping scheme and and
the discretization of the controls with interval-wise constant functions that
assume the values in $W$ on the intervals.

The tracking-type term and the squared $L^{2}$-norm of the control, are
discretized using the trapezoidal rule for the same intervals. The derivative
of the of the first (reduced) term of the objective is required to evaluate
the linear part of the objective of the subproblem (TR). In order to compute
the latter, we use a first discretize, then optimize-based [26] adjoint
calculus.

In our executions of Algorithm 1 on a computer, we replace the infinite-
dimensional trust-region subproblems with the discretizations that are
described above. We note that we have the (implicit) termination criterion in
Algorithm 1 that the trust-region radius contracts to a value below $T/N$
because we operate with limited precision and a fixed discretization. In this
case, the linear integer program that arises after discretizing (TR) has only
one feasible point, namely the function $d=0$ with objective value $0$. Thus
we always run Algorithm 1 until this situation occurs. The reset trust-region
radius is $\Delta^{0}=0.25T$. The acceptance value for the ratio of actual
over predicted reduction is $\sigma=10^{-3}$.

All computations are executed on a workstation with an AMD Epic 7742 CPU and
96 GB RAM.

###  6.2 Experiments

We conduct the following experiments for both choices $\delta=0$ and
$\delta=1$, the choices $\lambda=10^{-2}$, $\lambda=10^{-4}$, and $\lambda=0$,
and the PDE input choices $\psi=w$ ($F=0$, $S=0$, $g=0$) and $S=\chi w$ for a
fixed function $\chi$ ($F=0$, $g=0$, $\psi=0$), where $w$ denotes the control
function. We tabulate for which of the settings Assumption 4.1 is violated or
satisfied in table 1 according to the results obtained in section 5.

Table 1: Satisfaction (T) and violation (F) of Assumption 4.1 for the different experiments. |  | Input is $S$ |  |  | Input is $\psi$ |
---|---|---|---|---|---|---
$\lambda=$ | $0$ | $10^{-4}$ | $10^{-2}$ | $0$ | $10^{-4}$ | $10^{-2}$
$\delta=0$ | F | T | T | F | T | T
$\delta=1$ | F | F | F | F | F | F

We select six feasible initial controls $w^{0}$, specifically
$w^{0}\equiv\omega$ for all $\omega\in W$. Then we execute our implementation
of Algorithm 1 to solve (P) for all of these initializations $w^{0}$. We
record the final control $w^{f}$, the final objective value
$j(w^{f})+\alpha\operatorname*{TV}(w^{f})$, and the instationarity $C(w^{f})$
on termination for all of these executions.

Then we regularize $j$ by composing it with the application of a standard
mollifier to the control input following our recipe in section 4, that is we
replace $j$ by $j\circ K_{\varepsilon}$ in (P) and $\nabla j$ by
$K_{\varepsilon}^{*}(\nabla j)\circ K_{\varepsilon}$ in (TR). For each of the
six initial controls $w^{0}$ we execute a homotopy of $\varepsilon$ and
Algorithm 1 on the regularized problems with the regularization parameter
values
$\varepsilon\in\\{$1.6\text{\times}{10}^{-2}$,$8\text{\times}{10}^{-3}$,$4\text{\times}{10}^{-3}$,$2\text{\times}{10}^{-3}$,$1\text{\times}{10}^{-3}$,$5\text{\times}{10}^{-4}$,$2.5\text{\times}{10}^{-4}$,0\\}$.
We initialize Algorithm 1 with $w^{0}$ for the largest regularization
parameter value $\varepsilon=$1.6\text{\times}{10}^{-2}$$ and initialize the
execution of Algorithm 1 for the subsequent value of $\varepsilon$ with the
final control function iterate of the previous parameter value for
$\varepsilon$. Again, we record the final controls, objective values, and
instationarities on termination.

For $\lambda>0$ L-stationarity is not known to be a necessary optimality
condition. Moreover, the fixed discretization also implies that we cannot
expect that the final iterate is perfectly L-stationary even if Assumption 4.1
is satisfied. In order to provide a full picture, we have chosen to still
measure remaining instationarity for the final control iterates and report the
final instationarities for the unregularized optimization and the homotopy but
kindly ask the reader to take these values with caution. We assess the
remaining instationarity of the final control iterate $w^{f}$ by evaluating
$C(w^{f})\coloneqq\|(\nabla j(w^{f})(t_{i}))_{i=1}^{\\#s}\|$, where $t_{1}$,
$\ldots$, $t_{\\#s}$ for some $\\#s\in\mathbb{N}$ denote the switching times
of $w^{f}$, that is the values $t\in(0,T)$ such that $\lim_{t\downarrow
t_{\\#}s}w^{f}(t)\neq\lim_{t\uparrow t_{\\#}s}w^{f}(t)$.

###  6.3 Results

Over all values of $\lambda$ and for both input choices $\psi$ and $S$, we
observe a different picture between $\delta=0$ and $\delta=1$ when comparing
the results achieved with the unregularized optimization and the homotopy.

##### Input choice $\psi$.

The unregularized optimization terminates after consuming between 21 and 64
iterations for $\delta=0$ and consuming between 24 and 67 iterations for
$\delta=1$. The homotopy consumes a cumulative number of iterations between 35
and 90 iterations for $\delta=0$ and consumes a much higher cumulative number
between 194 and 349 iterations for $\delta=1$.

The objective values with the unregularized optimization problem are very
similar to the objective values produced by the homotopy for $\delta=0$ with
relative differences generally below $10^{-4}$. This is different for
$\delta=1$, where the homotopy generally achieves lower objective values with
relative differences generally around $10^{-3}$.

For $\delta=0$ the remaining instationarities are generally similar and of the
same order of magnitude for the unregularized optimization and the homotopy.
For $\delta=1$, the final instationarities obtained with the unregularized
optimization are generally (but not in all cases) between one and two orders
of magnitude higher.

We report the details table 2 for the poro-elastic case $\delta=0$ and in
table 3 for the poro-elastic case. Detailed iteration numbers over the
different values of the homotopy are given in table 4 for $\delta=0$ and table
5 for $\delta=1$.

Table 2: Control function $\psi$, poro-elastic case $\delta=0$: Objective, instationarity, and (cumulative) outer iterations of (P) for six initial controls achieved by executing Algorithm 1 on the discretized problem as well as by executing a homotopy, abbreviated (H), of executions of Algorithm 1 on regularized problems with regularization parameter driven to zero for the same initial controls. Ini. | Final Obj. | Final Obj. (H) | Final Inst. | Final Inst. (H) | Iter. | Cum. Iter. (H)
---|---|---|---|---|---|---
$\lambda=0$ 1 | $1.3474$ | $1.3474$ | $3.4\text{\times}{10}^{-06}$ | $2.2\text{\times}{10}^{-06}$ | 35 | 43
2 | $1.3474$ | $1.3474$ | $2.6\text{\times}{10}^{-06}$ | $2.1\text{\times}{10}^{-06}$ | 64 | 81
3 | $1.3474$ | $1.3474$ | $3.8\text{\times}{10}^{-06}$ | $2.3\text{\times}{10}^{-06}$ | 33 | 78
4 | $1.3475$ | $1.3474$ | $3.0\text{\times}{10}^{-06}$ | $1.4\text{\times}{10}^{-06}$ | 22 | 48
5 | $1.3474$ | $1.3474$ | $4.3\text{\times}{10}^{-06}$ | $2.5\text{\times}{10}^{-06}$ | 21 | 56
6 | $1.3474$ | $1.3474$ | $1.9\text{\times}{10}^{-06}$ | $1.3\text{\times}{10}^{-06}$ | 50 | 90
$\lambda=10^{-4}$ 1 | $1.3476$ | $1.3476$ | $3.0\text{\times}{10}^{-06}$ | $2.8\text{\times}{10}^{-06}$ | 35 | 42
2 | $1.3476$ | $1.3476$ | $3.4\text{\times}{10}^{-06}$ | $2.9\text{\times}{10}^{-06}$ | 47 | 79
3 | $1.3476$ | $1.3476$ | $2.1\text{\times}{10}^{-06}$ | $2.1\text{\times}{10}^{-06}$ | 38 | 84
4 | $1.3477$ | $1.3476$ | $3.4\text{\times}{10}^{-06}$ | $2.9\text{\times}{10}^{-06}$ | 22 | 43
5 | $1.3476$ | $1.3476$ | $2.3\text{\times}{10}^{-06}$ | $3.0\text{\times}{10}^{-06}$ | 26 | 71
6 | $1.3476$ | $1.3476$ | $2.5\text{\times}{10}^{-06}$ | $2.2\text{\times}{10}^{-06}$ | 49 | 69
$\lambda=10^{-2}$ 1 | $1.3616$ | $1.3616$ | $3.5\text{\times}{10}^{-06}$ | $2.2\text{\times}{10}^{-06}$ | 31 | 47
2 | $1.3616$ | $1.3616$ | $7.2\text{\times}{10}^{-06}$ | $3.8\text{\times}{10}^{-06}$ | 33 | 56
3 | $1.3616$ | $1.3616$ | $9.0\text{\times}{10}^{-06}$ | $3.6\text{\times}{10}^{-06}$ | 36 | 67
4 | $1.3616$ | $1.3615$ | $2.7\text{\times}{10}^{-06}$ | $3.4\text{\times}{10}^{-06}$ | 26 | 35
5 | $1.3616$ | $1.3616$ | $4.5\text{\times}{10}^{-06}$ | $2.1\text{\times}{10}^{-06}$ | 27 | 55
6 | $1.3616$ | $1.3617$ | $8.3\text{\times}{10}^{-06}$ | $5.4\text{\times}{10}^{-06}$ | 39 | 70
Table 3: Control function $\psi$, poro-visco-elastic case $\delta=1$: Objective, instationarity, and (cumulative) outer iterations of (P) for six initial controls achieved by executing Algorithm 1 on the discretized problem as well as by executing a homotopy, abbreviated (H), of executions of Algorithm 1 on regularized problems with regularization parameter driven to zero for the same initial controls. Ini. | Final Obj. | Final Obj. (H) | Final Inst. | Final Inst. (H) | Iter. | Cum. Iter. (H)
---|---|---|---|---|---|---
$\lambda=0$ 1 | $1.3647$ | $1.3625$ | $1.7\text{\times}{10}^{-04}$ | $4.0\text{\times}{10}^{-05}$ | 48 | 265
2 | $1.3646$ | $1.3624$ | $1.1\text{\times}{10}^{-04}$ | $1.3\text{\times}{10}^{-06}$ | 62 | 289
3 | $1.3657$ | $1.3624$ | $1.1\text{\times}{10}^{-04}$ | $1.3\text{\times}{10}^{-06}$ | 26 | 268
4 | $1.3624$ | $1.3624$ | $1.3\text{\times}{10}^{-06}$ | $1.3\text{\times}{10}^{-06}$ | 24 | 201
5 | $1.3652$ | $1.3624$ | $7.5\text{\times}{10}^{-05}$ | $1.3\text{\times}{10}^{-06}$ | 34 | 240
6 | $1.3634$ | $1.3624$ | $1.1\text{\times}{10}^{-04}$ | $1.3\text{\times}{10}^{-06}$ | 43 | 303
$\lambda=10^{-4}$ 1 | $1.3647$ | $1.3625$ | $1.7\text{\times}{10}^{-04}$ | $1.1\text{\times}{10}^{-06}$ | 40 | 232
2 | $1.3644$ | $1.3625$ | $1.1\text{\times}{10}^{-04}$ | $1.1\text{\times}{10}^{-06}$ | 57 | 278
3 | $1.3656$ | $1.3625$ | $1.1\text{\times}{10}^{-04}$ | $1.1\text{\times}{10}^{-06}$ | 31 | 262
4 | $1.3625$ | $1.3625$ | $1.1\text{\times}{10}^{-06}$ | $1.1\text{\times}{10}^{-06}$ | 22 | 207
5 | $1.3652$ | $1.3625$ | $7.5\text{\times}{10}^{-05}$ | $1.1\text{\times}{10}^{-06}$ | 34 | 264
6 | $1.3635$ | $1.3625$ | $1.1\text{\times}{10}^{-04}$ | $1.1\text{\times}{10}^{-06}$ | 43 | 312
$\lambda=10^{-2}$ 1 | $1.3732$ | $1.3697$ | $2.4\text{\times}{10}^{-04}$ | $5.8\text{\times}{10}^{-05}$ | 35 | 240
2 | $1.3715$ | $1.3696$ | $9.2\text{\times}{10}^{-05}$ | $1.8\text{\times}{10}^{-06}$ | 67 | 349
3 | $1.3730$ | $1.3696$ | $1.1\text{\times}{10}^{-04}$ | $1.8\text{\times}{10}^{-06}$ | 26 | 241
4 | $1.3696$ | $1.3696$ | $1.8\text{\times}{10}^{-06}$ | $1.8\text{\times}{10}^{-06}$ | 22 | 194
5 | $1.3716$ | $1.3696$ | $6.7\text{\times}{10}^{-05}$ | $1.8\text{\times}{10}^{-06}$ | 41 | 283
6 | $1.3700$ | $1.3696$ | $6.3\text{\times}{10}^{-05}$ | $1.8\text{\times}{10}^{-06}$ | 64 | 349
Table 4: Control function $\psi$, poro-elastic case $\delta=0$: Number of iterations required by the executions of Algorithm 1 forthe different initializations of the computational example over the different values of $\varepsilon$ of the homotopy, cumulative for the homotopy (H), and for the unregularized problem (U). | $\varepsilon=$ |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---
In. | ${0.016}$ | ${0.008}$ | ${0.004}$ | ${0.002}$ | ${0.001}$ | ${0.0005}$ | ${0.00025}$ | ${0}$ | (H) | (U)
$\lambda=0$ 1 | 30 | 4 | 2 | 1 | 2 | 1 | 1 | 2 | 43 | 35
2 | 53 | 13 | 5 | 3 | 2 | 1 | 1 | 3 | 81 | 64
3 | 66 | 2 | 3 | 2 | 1 | 1 | 1 | 2 | 78 | 33
4 | 29 | 8 | 2 | 2 | 3 | 1 | 1 | 2 | 48 | 22
5 | 41 | 2 | 2 | 6 | 1 | 1 | 1 | 2 | 56 | 21
6 | 79 | 2 | 2 | 2 | 1 | 1 | 1 | 2 | 90 | 50
$\lambda=10^{-4}$ 1 | 31 | 2 | 2 | 2 | 1 | 1 | 1 | 2 | 42 | 35
2 | 57 | 5 | 7 | 5 | 2 | 1 | 1 | 1 | 79 | 47
3 | 60 | 2 | 12 | 3 | 2 | 1 | 1 | 3 | 84 | 38
4 | 23 | 4 | 6 | 3 | 4 | 1 | 1 | 1 | 43 | 22
5 | 49 | 2 | 2 | 4 | 8 | 1 | 1 | 4 | 71 | 26
6 | 58 | 2 | 2 | 2 | 1 | 1 | 1 | 2 | 69 | 49
$\lambda=10^{-2}$ 1 | 28 | 8 | 4 | 2 | 1 | 1 | 1 | 2 | 47 | 31
2 | 35 | 4 | 7 | 3 | 2 | 1 | 1 | 3 | 56 | 33
3 | 47 | 7 | 6 | 1 | 2 | 1 | 1 | 2 | 67 | 36
4 | 22 | 2 | 3 | 3 | 2 | 1 | 1 | 1 | 35 | 26
5 | 41 | 2 | 3 | 3 | 1 | 1 | 1 | 3 | 55 | 27
6 | 50 | 6 | 6 | 1 | 2 | 1 | 1 | 3 | 70 | 39
Table 5: Control function $\psi$, poro-visco-elastic case $\delta=1$: Number of iterations required by the executions of Algorithm 1 for the different initializations of the computational example over the different values of $\varepsilon$ of the homotopy, cumulative for the homotopy (H), and for the unregularized problem (U). | $\varepsilon=$ |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---
In. | ${0.016}$ | ${0.008}$ | ${0.004}$ | ${0.002}$ | ${0.001}$ | ${0.0005}$ | ${0.00025}$ | ${0}$ | (H) | (U)
$\lambda=0$ 1 | 112 | 40 | 35 | 26 | 8 | 1 | 1 | 9 | 232 | 40
2 | 134 | 64 | 38 | 15 | 17 | 1 | 1 | 8 | 278 | 57
3 | 117 | 45 | 50 | 18 | 14 | 1 | 1 | 16 | 262 | 31
4 | 62 | 75 | 34 | 17 | 12 | 1 | 1 | 5 | 207 | 22
5 | 132 | 50 | 37 | 26 | 7 | 1 | 1 | 10 | 264 | 34
6 | 145 | 73 | 43 | 23 | 19 | 1 | 1 | 7 | 312 | 43
$\lambda=10^{-4}$ 1 | 112 | 40 | 35 | 26 | 8 | 1 | 1 | 9 | 232 | 40
2 | 134 | 64 | 38 | 15 | 17 | 1 | 1 | 8 | 278 | 57
3 | 117 | 45 | 50 | 18 | 14 | 1 | 1 | 16 | 262 | 31
4 | 62 | 75 | 34 | 17 | 12 | 1 | 1 | 5 | 207 | 22
5 | 132 | 50 | 37 | 26 | 7 | 1 | 1 | 10 | 264 | 34
6 | 145 | 73 | 43 | 23 | 19 | 1 | 1 | 7 | 312 | 43
$\lambda=10^{-2}$ 1 | 102 | 58 | 43 | 19 | 11 | 1 | 1 | 5 | 240 | 35
2 | 194 | 61 | 43 | 30 | 13 | 1 | 1 | 6 | 349 | 67
3 | 119 | 36 | 44 | 29 | 7 | 1 | 1 | 4 | 241 | 26
4 | 73 | 45 | 39 | 16 | 15 | 1 | 1 | 4 | 194 | 22
5 | 111 | 60 | 57 | 30 | 17 | 1 | 1 | 6 | 283 | 41

In order to also give a qualitative expression of the produced controls, we
visualize them in figs. 2 and 1 for $\delta=0$, $\delta=1$ and
$\lambda=10^{-2}$.

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-5$$-4.5$$-4$$-3.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\quad{\tiny(w^{0}\equiv-7)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-7$$-6$$-5$$-4$$-3$$-2$$-1$$0$$t\quad{\tiny(w^{0}\equiv-5)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-7$$-6$$-5$$-4$$-3$$-2$$-1$$0$$t\enskip(w^{0}\equiv-3)$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-7$$-6$$-5$$-4$$-3$$-2$$-1$$0$$t\quad{\tiny(w^{0}\equiv-1)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-7$$-6$$-5$$-4$$-3$$-2$$-1$$0$$t\quad{\tiny(w^{0}\equiv
0)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-7$$-6$$-5$$-4$$-3$$-2$$-1$$0$$t\enskip(w^{0}\equiv
2)$

Figure 1: Control function $\psi$, poro-elastic case $\delta=0$,
$\lambda=10^{-2}$: final control functions produced for the unregularized
optimization (solid) and for the homotopy (dashed).

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\quad{\tiny(w^{0}\equiv-7)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\quad{\tiny(w^{0}\equiv-5)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\enskip(w^{0}\equiv-3)$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\quad{\tiny(w^{0}\equiv-1)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\quad{\tiny(w^{0}\equiv
0)}$

$0$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$-4$$-3.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$t\enskip(w^{0}\equiv
2)$

Figure 2: Control function $\psi$, poro-visco-elastic case $\delta=1$,
$\lambda=10^{-2}$: final control functions produced for the unregularized
optimization (solid) and for the homotopy (dashed).

##### Input choice $S$.

The unregularized optimization terminates after consuming between 12 and 29
iterations for $\delta=0$ and consuming between 8 and 21 iterations for
$\delta=1$. The homotopy consumes a cumulative number of iterations between 23
and 40 iterations for $\delta=0$ and consumes a much higher cumulative number
between 42 and 58 iterations for $\delta=1$.

The objective values with the unregularized optimization problem have relative
differences generally around $10^{-2}$ compared to the objective values
produced by the homotopy for $\delta=0$. This is different for $\delta=1$,
where the homotopy generally achieves lower objective values with relative
differences generally (but not in all cases) higher than $10^{-1}$.

For $\delta=0$ the remaining instationarities are generally similar and of the
same order of magnitude for the unregularized optimization and the homotopy.
For $\delta=1$, the final instationarities obtained with the unregularized
optimization are generally (but not in all cases) between one and two orders
of magnitude higher.

We report the details table 6 for the poro-elastic case $\delta=0$ and in
table 7 for the poro-elastic case. Detailed iteration numbers over the
different values of the homotopy are given in table 8 for $\delta=0$ and table
9 for $\delta=1$.

Table 6: Control function $S$, poro-elastic case $\delta=0$: Objective, instationarity, and (cumulative) outer iterations of (P) for six initial controls achieved by executing Algorithm 1 on the discretized problem as well as by executing a homotopy, abbreviated (H), of executions of Algorithm 1 on regularized problems with regularization parameter driven to zero for the same initial controls. Ini. | Final Obj. | Final Obj. (H) | Final Inst. | Final Inst. (H) | Iter. | Cum. Iter. (H)
---|---|---|---|---|---|---
$\lambda=0$ 1 | $1.2420$ | $1.2417$ | $1.0\text{\times}{10}^{-05}$ | $9.1\text{\times}{10}^{-06}$ | 16 | 24
2 | $1.2423$ | $1.2421$ | $1.5\text{\times}{10}^{-05}$ | $1.2\text{\times}{10}^{-05}$ | 25 | 33
3 | $1.2428$ | $1.2424$ | $8.8\text{\times}{10}^{-06}$ | $1.2\text{\times}{10}^{-05}$ | 24 | 32
4 | $1.2417$ | $1.2428$ | $1.1\text{\times}{10}^{-05}$ | $1.5\text{\times}{10}^{-05}$ | 15 | 27
5 | $1.2414$ | $1.2423$ | $1.1\text{\times}{10}^{-05}$ | $1.6\text{\times}{10}^{-05}$ | 12 | 23
6 | $1.2553$ | $1.2556$ | $1.1\text{\times}{10}^{-05}$ | $2.5\text{\times}{10}^{-05}$ | 16 | 28
$\lambda=10^{-4}$ 1 | $1.2420$ | $1.2417$ | $1.0\text{\times}{10}^{-05}$ | $9.2\text{\times}{10}^{-06}$ | 16 | 24
2 | $1.2423$ | $1.2422$ | $1.5\text{\times}{10}^{-05}$ | $1.2\text{\times}{10}^{-05}$ | 25 | 33
3 | $1.2429$ | $1.2424$ | $8.7\text{\times}{10}^{-06}$ | $1.3\text{\times}{10}^{-05}$ | 26 | 32
4 | $1.2417$ | $1.2428$ | $1.1\text{\times}{10}^{-05}$ | $1.5\text{\times}{10}^{-05}$ | 15 | 27
5 | $1.2414$ | $1.2424$ | $1.1\text{\times}{10}^{-05}$ | $1.6\text{\times}{10}^{-05}$ | 12 | 23
6 | $1.2554$ | $1.2556$ | $1.1\text{\times}{10}^{-05}$ | $2.5\text{\times}{10}^{-05}$ | 16 | 28
$\lambda=10^{-2}$ 1 | $1.2454$ | $1.2451$ | $2.0\text{\times}{10}^{-05}$ | $1.9\text{\times}{10}^{-05}$ | 16 | 24
2 | $1.2602$ | $1.2455$ | $4.4\text{\times}{10}^{-07}$ | $2.4\text{\times}{10}^{-06}$ | 29 | 40
3 | $1.2462$ | $1.2458$ | $9.5\text{\times}{10}^{-07}$ | $2.2\text{\times}{10}^{-05}$ | 26 | 32
4 | $1.2450$ | $1.2462$ | $2.1\text{\times}{10}^{-05}$ | $2.5\text{\times}{10}^{-05}$ | 15 | 29
5 | $1.2448$ | $1.2457$ | $2.1\text{\times}{10}^{-05}$ | $6.3\text{\times}{10}^{-06}$ | 12 | 24
6 | $1.2602$ | $1.2605$ | $1.1\text{\times}{10}^{-05}$ | $1.0\text{\times}{10}^{-05}$ | 13 | 29
Table 7: Control function $S$, poro-visco-elastic case $\delta=1$: Objective, instationarity, and (cumulative) outer iterations of (P) for six initial controls achieved by executing Algorithm 1 on the discretized problem as well as by executing a homotopy, abbreviated (H), of executions of Algorithm 1 on regularized problems with regularization parameter driven to zero for the same initial controls. Ini. | Final Obj. | Final Obj. (H) | Final Inst. | Final Inst. (H) | Iter. | Cum. Iter. (H)
---|---|---|---|---|---|---
$\lambda=0$ 1 | $1.8162$ | $1.4783$ | $3.6\text{\times}{10}^{-03}$ | $5.2\text{\times}{10}^{-05}$ | 10 | 45
2 | $1.5849$ | $1.4930$ | $1.2\text{\times}{10}^{-03}$ | $2.3\text{\times}{10}^{-04}$ | 20 | 57
3 | $1.5480$ | $1.5060$ | $9.7\text{\times}{10}^{-04}$ | $1.5\text{\times}{10}^{-04}$ | 17 | 53
4 | $2.0874$ | $1.5176$ | $3.9\text{\times}{10}^{-03}$ | $3.1\text{\times}{10}^{-04}$ | 15 | 49
5 | $1.7936$ | $1.4807$ | $1.0\text{\times}{10}^{-03}$ | $3.6\text{\times}{10}^{-05}$ | 8 | 46
6 | $1.4730$ | $1.4713$ | $3.6\text{\times}{10}^{-04}$ | $2.8\text{\times}{10}^{-04}$ | 14 | 45
$\lambda=10^{-4}$ 1 | $1.8162$ | $1.4783$ | $3.6\text{\times}{10}^{-03}$ | $5.2\text{\times}{10}^{-05}$ | 10 | 45
2 | $1.5849$ | $1.4930$ | $1.2\text{\times}{10}^{-03}$ | $2.3\text{\times}{10}^{-04}$ | 20 | 57
3 | $1.5481$ | $1.5060$ | $9.7\text{\times}{10}^{-04}$ | $1.5\text{\times}{10}^{-04}$ | 17 | 53
4 | $2.0874$ | $1.5177$ | $3.9\text{\times}{10}^{-03}$ | $3.1\text{\times}{10}^{-04}$ | 15 | 49
5 | $1.7936$ | $1.4808$ | $1.0\text{\times}{10}^{-03}$ | $3.6\text{\times}{10}^{-05}$ | 8 | 46
6 | $1.4731$ | $1.4713$ | $3.6\text{\times}{10}^{-04}$ | $2.8\text{\times}{10}^{-04}$ | 14 | 45
$\lambda=10^{-2}$ 1 | $1.8185$ | $1.4814$ | $3.5\text{\times}{10}^{-03}$ | $4.2\text{\times}{10}^{-05}$ | 10 | 42
2 | $1.5755$ | $1.4971$ | $1.0\text{\times}{10}^{-03}$ | $2.2\text{\times}{10}^{-04}$ | 21 | 58
3 | $1.6598$ | $1.5092$ | $1.5\text{\times}{10}^{-03}$ | $1.4\text{\times}{10}^{-04}$ | 17 | 58
4 | $2.0279$ | $1.4846$ | $3.6\text{\times}{10}^{-03}$ | $7.3\text{\times}{10}^{-05}$ | 18 | 50
5 | $1.8013$ | $1.4838$ | $1.2\text{\times}{10}^{-03}$ | $2.6\text{\times}{10}^{-05}$ | 10 | 49
6 | $1.4776$ | $1.4752$ | $3.6\text{\times}{10}^{-04}$ | $2.2\text{\times}{10}^{-04}$ | 14 | 46
Table 8: Control function $S$, poro-elastic case $\delta=0$: Number of iterations required by the executions of Algorithm 1 forthe different initializations of the computational example over the different values of $\varepsilon$ of the homotopy, cumulative for the homotopy (H), and for the unregularized problem (U). | $\varepsilon=$ |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---
In. | ${0.016}$ | ${0.008}$ | ${0.004}$ | ${0.002}$ | ${0.001}$ | ${0.0005}$ | ${0.00025}$ | ${0}$ | (H) | (U)
$\lambda=0$ 1 | 14 | 2 | 1 | 1 | 1 | 1 | 2 | 2 | 24 | 16
2 | 20 | 1 | 1 | 2 | 2 | 3 | 1 | 3 | 33 | 25
3 | 19 | 3 | 1 | 2 | 3 | 1 | 1 | 2 | 32 | 24
4 | 16 | 2 | 1 | 2 | 3 | 1 | 1 | 1 | 27 | 15
5 | 11 | 1 | 2 | 3 | 2 | 2 | 1 | 1 | 23 | 12
6 | 15 | 2 | 2 | 1 | 4 | 2 | 1 | 1 | 28 | 16
$\lambda=10^{-4}$ 1 | 14 | 2 | 1 | 1 | 1 | 1 | 2 | 2 | 24 | 16
2 | 20 | 1 | 1 | 2 | 2 | 3 | 1 | 3 | 33 | 25
3 | 19 | 3 | 1 | 2 | 3 | 1 | 1 | 2 | 32 | 26
4 | 16 | 2 | 1 | 2 | 3 | 1 | 1 | 1 | 27 | 15
5 | 11 | 1 | 2 | 3 | 2 | 2 | 1 | 1 | 23 | 12
6 | 15 | 2 | 2 | 1 | 4 | 2 | 1 | 1 | 28 | 16
$\lambda=10^{-2}$ 1 | 14 | 2 | 1 | 1 | 1 | 1 | 2 | 2 | 24 | 16
2 | 25 | 3 | 1 | 3 | 2 | 2 | 1 | 3 | 40 | 29
3 | 19 | 2 | 2 | 2 | 3 | 1 | 1 | 2 | 32 | 26
4 | 16 | 3 | 1 | 3 | 3 | 1 | 1 | 1 | 29 | 15
5 | 12 | 1 | 2 | 2 | 2 | 3 | 1 | 1 | 24 | 12
6 | 12 | 2 | 3 | 1 | 3 | 1 | 4 | 3 | 29 | 13
Table 9: Control function $S$, poro-visco-elastic case $\delta=1$: Number of iterations required by the executions of Algorithm 1 for the different initializations of the computational example over the different values of $\varepsilon$ of the homotopy, cumulative for the homotopy (H), and for the unregularized problem (U). | $\varepsilon=$ |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---
In. | ${0.016}$ | ${0.008}$ | ${0.004}$ | ${0.002}$ | ${0.001}$ | ${0.0005}$ | ${0.00025}$ | ${0}$ | (H) | (U)
$\lambda=0$ 1 | 21 | 4 | 9 | 2 | 1 | 3 | 3 | 2 | 45 | 10
2 | 36 | 4 | 4 | 3 | 3 | 3 | 2 | 2 | 57 | 20
3 | 25 | 16 | 2 | 2 | 2 | 1 | 3 | 2 | 53 | 17
4 | 31 | 4 | 3 | 3 | 1 | 2 | 3 | 2 | 49 | 15
5 | 31 | 2 | 1 | 3 | 3 | 2 | 2 | 2 | 46 | 8
6 | 26 | 3 | 2 | 2 | 5 | 3 | 2 | 2 | 45 | 14
$\lambda=10^{-4}$ 1 | 21 | 4 | 9 | 2 | 1 | 3 | 3 | 2 | 45 | 10
2 | 36 | 4 | 4 | 3 | 3 | 3 | 2 | 2 | 57 | 20
3 | 25 | 16 | 2 | 2 | 2 | 1 | 3 | 2 | 53 | 17
4 | 31 | 4 | 3 | 3 | 1 | 2 | 3 | 2 | 49 | 15
5 | 31 | 2 | 1 | 3 | 3 | 2 | 2 | 2 | 46 | 8
6 | 26 | 3 | 2 | 2 | 5 | 3 | 2 | 2 | 45 | 14
$\lambda=10^{-2}$ 1 | 20 | 10 | 2 | 1 | 1 | 3 | 3 | 2 | 42 | 10
2 | 38 | 4 | 3 | 3 | 3 | 3 | 2 | 2 | 58 | 21
3 | 29 | 17 | 2 | 2 | 2 | 1 | 3 | 2 | 58 | 17
4 | 30 | 4 | 3 | 3 | 1 | 2 | 3 | 4 | 50 | 18
5 | 31 | 2 | 1 | 3 | 3 | 3 | 4 | 2 | 49 | 10
6 | 27 | 3 | 2 | 2 | 5 | 3 | 2 | 2 | 46 | 14

##  7 Conclusion

We have analyzed the regularity condition Assumption 4.1 that, when imposed on
(P) allows the convergence analysis for Algorithm 1 in [33], for control
problems that are governed by poro-elastic and poro-visco-elastic equations
modeling fluid flow through porous media. We have shown for one-dimensional
spatial domains and two possible control inputs (one acting in the interior
and acting on the boundary) that the regularity conditions are violated when
the viscosity parameter $\delta$ is chosen strictly greater than zero or a
Thikonov term is present and that the regularity conditions are satisfied for
$\delta=0$ and no Thikonov parameter being present all else equal.

We have applied Algorithm 1 to instances of (P) for differently scaled
Thikonov terms, for $\delta=0$ and $\delta=1$, and for the two different
control inputs that were analyzed before.

We have observed that the presence of the Thikonov term does not seem to
negatively impact the practical performance of Algorithm 1 although Assumption
4.1 is always violated in this case. If Assumption 4.1 is violated due to the
choice $\delta=1$, then the performance of Algorithm 1 and the quality of the
final iterates it produces before the trust region collapses is degraded. This
can be alleviated by executing a homotopy that drives the support parameter of
a mollification of the control input into the PDE to zero over the course of
the optimization. The execution of the homotopy comes at a higher
computational cost, however.

## References

  * [1] L. Ambrosio, N. Fusco, and D. Pallara.  Functions of bounded variation and free discontinuity problems, volume 254 of Oxford Mathematical Monographs.  Clarendon Press Oxford, 2000.
  * [2] T. M. Apostol.  Mathematical Analysis.  Addison-Wesley Longman, 1974.
  * [3] R. P. Araujo and D. L. Sean McElwain.  A mixture theory for the genesis of residual stresses in growing tissues I: a general formulation.  SIAM J. Appl. Math. 65(4):1261–1284, 2005.
  * [4] F. Bestehorn, C. Hansknecht, C. Kirches, and P. Manns.  A switching cost aware rounding method for relaxations of mixed-integer optimal control problems.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 7134–7139. IEEE, 2019.
  * [5] F. Bestehorn, C. Hansknecht, C. Kirches, and P. Manns.  Mixed-integer optimal control problems with switching costs: a shortest path approach.  Mathematical Programming, 188(2):621–652, 2021.
  * [6] M.A. Biot.  General theory of three-dimensional consolidation.  J. Appl. Phys, 12(2) pp. 155–164, 1941.
  * [7] L. Bociu, G. Guidoboni, R. Sacco, and M. Verri.  On the role of compressibility in poroviscoelastic models.  Mathematical Biosciences and Engineering, 16(5):6167–6208, 2019\.
  * [8] L. Bociu, G. Guidoboni, R. Sacco, and J. Webster.  Analysis of nonlinear poro-elastic and poro-viscoelastic models.  Archive for Rational Mechanics and Analysis 222, 1445-1519, 2016\.
  * [9] L. Bociu and S. Strikwerda.  Poro-visco-elasticity in biomechanics - optimal control, AWM: Research in the Mathematics of Materials Science.  Springer, 2022.
  * [10] K. Bredies, K. Kunisch, and T. Pock.  Total generalized variation.  SIAM Journal on Imaging Sciences, 3(3):492–526, 2010.
  * [11] P. Causin, G. Guidoboni, A. Harris, D. Prada, R. Sacco, and S. Terragni.  A poroelastic model for the perfusion of the lamina cribrosa in the optic nerve head.  Math Biosci, pages 33–41, 2014.
  * [12] A. Chambolle and P.-L. Lions.  Image recovery via total variation minimization and related problems.  Numerische Mathematik, 76(2):167–188, 1997.
  * [13] T. F. Chan and C.-K. Wong.  Total variation blind deconvolution.  IEEE Transactions on Image Processing, 7(3):370–375, 1998.
  * [14] C. Clason, F. Kruse, and K. Kunisch.  Total variation regularization of multi-material topology optimization.  ESAIM: Mathematical Modelling and Numerical Analysis, 52(1):275–303, 2018.
  * [15] E. Detournay and A.H.-D. Cheng.  Fundamentals of poroelasticity, chapter 5 in comprehensive rock engineering: Principles, practice and projects, vol. ii, analysis and design method, ed. c. fairhurst.  Pergamon Press, 113-171, 1993.
  * [16] S. Engel, B. Vexler, and P. Trautmann.  Optimal finite element error estimates for an optimal control problem governed by the wave equation with controls of bounded variation.  IMA Journal of Numerical Analysis, 41(4):2639–2667, 2021.
  * [17] M. Fornasier and C.-B. Schönlieb.  Subspace correction methods for total variation and $\backslash$ell_1-minimization.  SIAM Journal on Numerical Analysis, 47(5):3397–3428, 2009.
  * [18] A. J. H. Frijns.  A four-component mixture theory applied to cartilaginous tissues: Numerical modelling and experiments.  Thesis (Dr.ir.)–Technische Universiteit Eindhoven (The Netherlands), 2000.
  * [19] M. Gerdts.  Solving mixed-integer optimal control problems by Branch&Bound: A case study from automobile test-driving with gear shift.  Optimal Control Applications and Methods, 26:1–18, 2005\.
  * [20] S. Göttlich, A. Potschka, and C. Teuber.  A partial outer convexification approach to control transmission lines.  Computational Optimization and Applications, 72(2):431–456, 2019\.
  * [21] S. Göttlich, A. Potschka, and U. Ziegler.  Partial outer convexification for traffic light optimization in road networks.  SIAM Journal on Scientific Computing, 39(1):B53–B75, 2017.
  * [22] F. M. Hante, G. Leugering, A. Martin, L. Schewe, and M. Schmidt.  Challenges in optimal control problems for gas and fluid flow in networks of pipes and canals: From modeling to industrial applications.  In Industrial mathematics and complex systems, pages 77–122. Springer, 2017.
  * [23] F. M. Hante and S. Sager.  Relaxation methods for mixed-integer optimal control of partial differential equations.  Computational Optimization and Applications, 55(1):197–225, 2013\.
  * [24] J. Haslinger and R. A. E. Mäkinen.  On a topology optimization problem governed by two-dimensional helmholtz equation.  Computational Optimization and Applications, 62(2):517–544, 2015\.
  * [25] M. Hintermüller and C. N. Rautenberg.  Optimal selection of the regularization function in a weighted total variation model. part i: Modelling and theory.  Journal of Mathematical Imaging and Vision, 59(3):498–514, 2017\.
  * [26] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.  Optimization with PDE constraints, volume 23.  Springer Science & Business Media, 2008.
  * [27] C. Y. Kaya.  Optimal control of the double integrator with minimum total variation.  Journal of Optimization Theory and Applications, 185:966–981, 2020\.
  * [28] C. Kirches, H. G. Georg Bock, J. P. Schlöder, and S. Sager.  Mixed-integer nmpc for predictive cruise control of heavy-duty trucks.  In 2013 European Control Conference (ECC), pages 4118–4123. IEEE, 2013.
  * [29] S. M. Klisch.  Internally constrained mixtures of elastic continua.  Math. Mech. Solids, 4:481–498, 1999.
  * [30] J. Lellmann, D. A. Lorenz, C.-B. Schonlieb, and T. Valkonen.  Imaging with kantorovich–rubinstein discrepancy.  SIAM Journal on Imaging Sciences, 7(4):2833–2859, 2014.
  * [31] G. Lemon, J. R. King, H. M. Byrne, O. E. Jensen, and K. M. Shakesheff.  Mathematical modelling of engineered tissue growth using a multiphase porous flow mixture theory.  J. Math. Biol., 52:571–594, 2006.
  * [32] S. Leyffer.  Integrating SQP and branch-and-bound for mixed integer nonlinear programming.  Computational Optimization and Applications, 18(3):295–309, 2001\.
  * [33] S. Leyffer and P. Manns.  Sequential linear integer programming for integer optimal control with total variation regularization.  ESAIM: Control, Optimisation and Calculus of Variations (ESAIM: COCV), 2021.
  * [34] R. Loxton, Q. Lin, V. Rehbock, and K. L. Teo.  Control parameterization for optimal control problems with continuous inequality constraints: New convergence results.  Numerical Algebra, Control and Optimization, 2(3):571–599, 2012\.
  * [35] F. Maggi.  Sets of Finite Perimeter and Geometric Variational Problems: An Introduction to Geometric Measure Theory.  Number 135. Cambridge University Press, 2012.
  * [36] P. Manns and C. Kirches.  Multidimensional sum-up rounding for elliptic control systems.  SIAM Journal on Numerical Analysis, 58(6):3427–3447, 2020.
  * [37] J. Marko and G. Wachsmuth.  Integer optimal control problems with total variation regularization: Optimality conditions and fast solution of subproblems.  arXiv preprint arXiv:2207.05503, 2022.
  * [38] A. Martin, M. Möller, and S. Moritz.  Mixed integer models for the stationary case of gas network optimization.  Mathematical Programming, 105(2):563–582, 2006.
  * [39] L. Preziosi and A. Tosin.  Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications.  J. Math. Biol., 58:625–656, 2009.
  * [40] L. I. Rudin, S. Osher, and E. Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena, 60(1-4):259–268, 1992.
  * [41] S. Sager.  Numerical methods for mixed-integer optimal control problems.  Der andere Verlag Tönning, Lübeck, Marburg, 2005.
  * [42] S. Sager, H. G. Bock, and M. Diehl.  The Integer Approximation Error in Mixed-Integer Optimal Control.  Mathematical Programming, 133(1–2):1–23, 2012.
  * [43] S. Sager, M. Jung, and C. Kirches.  Combinatorial integral approximation.  Mathematical Methods of Operations Research, 73(3):363–380, 2011.
  * [44] S. Sager and C. Zeile.  On mixed-integer optimal control with constrained total variation of the integer control.  submitted, 2019.
  * [45] M. Severitt and P. Manns.  Efficient solution of discrete subproblems arising in integer optimal control with total variation regularization.  arXiv preprint arXiv:2206.01642, 2022.
  * [46] O. Sigmund and K. Maute.  Topology optimization approaches.  Structural and Multidisciplinary Optimization, 48(6):1031–1055, 2013.
  * [47] M. Verri, G. Guidoboni, L. Bociu, and R. Sacco.  The role of structural viscoelasticity in deformable porous media with incompressible constituents: applications in biomechanics.  Mathematical Biosciences and Engineering, 15(4):933–959, 2018.
  * [48] C. R. Vogel and M. E. Oman.  Iterative methods for total variation denoising.  SIAM Journal on Scientific Computing, 17(1):227–238, 1996.
  * [49] E. Zeidler.  Applied Functional Analysis: Main Principles and Their Applications, volume 109 of Applied Mathematical Sciences.  Springer Science & Business Media, 2012.