Mathematical aspects of catalyst positioning in lithium/air batteries

The main factor limiting the efficiency of non-aqueous lithium/air batteries is the deposition of lithium-dioxide during the reduction process, which leads to a narrowing and finally clogging of the pores of the batteries, for some theoretical and experimental research on this topic see [19, 22]. In this paper we aim at optimizing the initial placement of catalysts at the pore surface in order to optimize the efficiency of such batteries. In this context efficiency is measured in terms of the remaining free pore volume after clogging.

Our starting point is the approach of [4, 21], where a system of non-linear differential equations has been introduced for describing the process of discharging and clogging in a single pore. The basic parameter in these models is the initial distribution cat₀(x) of catalysts along the surface of the pore. These previous approaches assumed, that either cat₀ is constant 1 [21], or that cat₀ can take only values 1 or 0 indicating the presence of a catalyst at certain positions.

Our first target, see Problem I below, is to state a continuous model for cat₀ and to analyse the properties of the resulting parameter-to-state mapping (existence, uniqueness, Fréchet derivative, sensitivity system, adjoint operator) as a prerequisite to attacking the optimization problem in section 4. Hence, following [4, 21] we define $\Omega=[0,1]$ and $Q_T = [0,T] \times \Omega$ as well as $c: Q_T \rightarrow [0,\infty)$ as the oxygen concentration at time t at position x of a normalized pore. By $r: Q_T \rightarrow [0,\infty)$ we denote the radius of the pore, which diminishes over time and finally leads to pore clogging. The first catalyst model assumes a catalyst distribution cat$_0: \Omega \rightarrow I\!R$, which is constant over time (Problem I):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial}{\partial t}(r^2c)=\frac{\partial}{\partial x} \Big(r^2 \frac{\partial c}{\partial x}\Big)-\beta r c, \, \, (t,x) \in Q_T, \label{subeq:I:1} \nonumber \end{align} \tag{ 2.1a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial r}{\partial t}=-\gamma \ {\rm cat}_0(x) \ c, \, \, (t,x)\in Q_T, \label{subeq:I:2} \nonumber \end{align} \tag{ 2.1b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(0,x)=r_0(x)=1, \, \, x\in \Omega,\label{subeq:I:3} \nonumber \end{align} \tag{ 2.1c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c(0,x)=h(x),\label{subeq:I:4} \nonumber \end{align} \tag{ 2.1d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c(t,0)=1, \quad c_x(t,1)=0, \quad t\in [0,T],\label{subeq:I:5} \nonumber \end{align} \tag{ 2.1e }$$

where $\beta=\beta(t,c,r)$ and $\gamma=\gamma(t,c,r)$ are given by equations (2.2) and (2.3). For the initial oxygen concentration c at time t  =  0, we assume the compatibility conditions $c(0,0)=h(0)=1$ and $c(0,1)=h(1)=0$, which models that at t  =  0 oxygen is entering the pore at x  =  0 but has not yet reached the end of the pore at x  =  1.

Remark 2.1. Following [21], the mathematical model (2.1a)–(2.1e) can be justified by considering the lithium/air reduction process in a cylindrical pore of radius $r_p^0$ and length l, a porous cathode of volume $V_{\rm cath}$, and given initial porosity $\newcommand{\e}{{\rm e}} \epsilon_0$. The authors introduce a coordinate system where t_p, resp. x_p, denotes time, resp. spatial position along the central axis of the pore. Furthermore, assume that the oxygen concentration $C_{\rm O_2}$, which is a function of t_p and x_p, is angularly symmetric about the diffusion direction, which is identical with the radial direction along the pore. The initial radius of the pore is given by $r_p^0$, it diminishes during the reduction process, its value at position x_p at time t_p is denoted by r_p.

Let $C_{\rm O_{2},0}= 2.10\cdot 10^{-6}~{\rm (mol~cm^{-3})}$ denote the initial oxygen concentration and $D_{\rm O_2}= 16.70\cdot 10^{-6}~{\rm (cm^2~s^{-1})}$ the oxygen mass diffusivity. This allows to introduce the dimensionless quantities $\newcommand{\e}{{\rm e}} t=\frac{t_pD_{\rm O_2}}{\sqrt{\epsilon_0} \ l^2}$, $r=\frac{r_p}{r_p^0}$, $x=\frac{x_p}l$, $c=\frac{C_{\rm O_2}}{C_{{\rm O_2},0}}$, i.e. $0\leqslant r,c,x \leqslant 1$, 0  <  t  <  T. The above stated model is obtained by introducing parameters $\gamma$ and $\beta$ given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:01.6*} \beta=\frac{2k\sqrt{\epsilon_0} l^2}{D_{\rm O_2}{r}^0_{p}},\quad \gamma=\frac{k\,C_{\rm O_2}\,M_{\rm Li_2O_2}/\rho_{\rm Li_2O_2}}{D_{\rm O_2}r^0_{p}/\sqrt{\epsilon_0} l^2}, \nonumber \end{align} \tag{ 2.2 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:01.6**} k=\frac{i_{{{\rm geom}}}\,r^0_{p}}{4F\,C_{\rm O_2}\epsilon_0\,l\,\int_0^1{r\,c\,{\rm d}x}}. \nonumber \end{align} \tag{ 2.3 }$$

We note that $\beta$ and $\gamma$ depend on k, whose definition involves an integral over c and r, hence they are functions of time t. In order to highlight this non-linear dependence we use the full notation $\beta=\beta(t, c,r)$ and $\gamma = \gamma (t,c,r)$ whenever suitable.

The other parameters are set as follows (see [4, 21]): $F=9.648\,533\,65 \cdot 10^7$ (Faraday constant), $r_{p}^0 = 2.0 \times 10^{-7}$ cm; l  =  0.07 cm; $\newcommand{\e}{{\rm e}} \epsilon_{0} = 0.73$; $M_{\rm Li_2O_2} = 45.8765~{\rm g~mol^{-1}}$; ${\rho}_{\rm Li_2O_2} = 2.3~{\rm g~cm^{-3}}$; ${i}_{\rm geom} = 0.1~{\rm mA~cm^{-2}}$.

Note that [4] discusses different exponents in the reaction equation (2.1a), namely,

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial}{\partial t}\big(r^{\frac43}c\big)=\frac{\partial}{\partial x}\Big(r^2 \frac{\partial c}{\partial x}\Big)-\beta(t,c,r) r^{\frac13} c, \, \, (t,x)\in Q_T. \label{subeq:II:1} \nonumber \end{align} \tag{ 2.4 }$$

In this case, the dimensionless time is $\newcommand{\e}{{\rm e}} t=\frac{t_pD_{\rm O_2}\sqrt{\epsilon_0}}{l^2}$ and parameters $\gamma$, $\beta$ are defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:01.6*:2} \beta=\frac{2k l^2}{D_{\rm O_2}\sqrt{\epsilon_0} {r}^0_{p}},\quad \gamma=\frac{3\,k\,l^2\,C_{\rm O_2}\,M_{\rm Li_2O_2}}{2\,D_{\rm O_2}\,\sqrt{\epsilon_0}\,r^0_{p}\, \rho_{\rm Li_2O_2}}, \nonumber \end{align} \tag{ 2.5 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:01.6**:2} k=\frac{i_{{{\rm geom}}}\,r^0_{p}}{4F\,C_{\rm O_2}\epsilon_0\,l\,\int_0^1{r^{\frac13}\,c\,{\rm d}x}}. \nonumber \end{align} \tag{ 2.6 }$$

In the following, we will state the mathematical analysis for the model (2.1a); however we emphasize, that the same methodology can be used to analyse the model (2.4) as well. For more details and the parameters concerned, we refer the readers to [4, 21].

Problem I is our basic model, which will be analysed as a parameter identification problem in the following sections. Finally, we aim to optimize the initial catalyst positioning cat₀, which we regard a major design decision for optimizing non-aqueous lithium/air batteries. In order to stress this dependence we use the notation $r({\rm cat}_0;t,x)$ as well as $c({\rm cat}_0;t,x)$. Depending on the context we will sometimes eliminate this dependence and simply use e.g. $r({\rm cat}_0;t,x)=r(t,x)=r$.

Numerical tests and simulations we will be based on an extended model by making ${\rm cat}={\rm cat}(t,x)$ a time dependent function with ${\rm cat}(0,\cdot)={\rm cat}_0$. This is in agreement with experiments and yields a more realistic model. Such a dynamic approach has already been incorporated in [4], which used a simple discretized evolution based on thresholding: assume that the initial catalysts distribution cat₀ takes value 1 at position x₀ but is 0 at neighboring discretization points; the deposition will start at point x₀ only; the value of the cat-function at neighboring points is changed to 1 whenever r(x₀) reaches a critical value. In the present paper we use a different approach, where the system of differential equations for describing the evolution of the pore radius r and the oxygen concentration c is expanded by a simple diffusion model for the time evolution of ${\rm cat}(t,x)$. Learning from growth kinetics, see e.g. [4] an the references there, the evolution of cat is based on incorporating a threshold model ($T_{R_0}$) in combination with a generic diffusion model (Gaussian). The threshold function models the fact, that once the radius has decreased below the level R₀ then the catalytic action is increased and spreads out to neighboring positions. Hence Problem II is obtained from Problem I by adapting the time derivative for r

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{01.3*} \frac{\partial r}{\partial t}=-\gamma(t,c,r) \ {\rm cat}(t,x) \ c, \, \, (t,x) \in Q_T, \nonumber \end{align} \tag{ 2.7 }$$

and by adding the following evolution for cat

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{01.9} {\rm cat}_t(t,x)=\alpha_0\int_0^1{T_{\chi,g,R_0}\big(r({\rm cat}_0;y,t)\big)G(x-y)}\,{\rm d}y, \nonumber \end{align} \tag{ 2.8 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{01.9b} {\rm cat}(0,x)={\rm cat}_0(x), \nonumber \end{align} \tag{ 2.9 }$$

where $\alpha_0$ is some small positive constant, R₀ is a positive constant, 0  <  R₀  <  1. Here $\newcommand{\e}{{\rm e}} G(x)=\frac{1}{2\pi}\exp(-\frac12\,x^2)$ denotes the Gaussian distribution and $T_{\chi,g,R_0}(r)$ is a C¹-threshold function given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.Huber*} T_{\chi,g,R_0}(s)= \left\{\begin{array}{@{}l@{}} 0, \, {{\rm if}} \, s > \frac{g}{2\chi}+\frac{1}{2\chi^2}+R_0 \nonumber \\ \frac{\chi}2\,(g-\chi\,s+\chi R_0+\frac1{2\chi})^2,\, {{\rm if}} \,\frac{g}{2\chi}-\frac{1}{2\chi^2}+R_0< s \leqslant \frac{g}{2\chi}+\frac{1}{2\chi^2}+R_0 \nonumber \\ g-\chi\,s+\chi R_0, {{\rm if}} - \frac{g}{2\chi}+\frac{1}{2\chi^2}+R_0< s \leqslant \frac{g}{2\chi}-\frac{1}{2\chi^2}+R_0 \nonumber \\ 2g-\frac{\chi}2\,(g+\chi\,s- \chi R_0+\frac1{2\chi})^2,\, {{\rm if}} -\frac{g}{2\chi}-\frac{1}{2\chi^2}+R_0< s \leqslant -\frac{g}{2\chi}+\frac{1}{2\chi^2}+R_0 \nonumber \\ 1, \, {{\rm if}} \, s\leqslant -\frac{g}{2\chi}-\frac{1}{2\chi^2}+R_0. \end{array}\right. \nonumber \end{align} \tag{ 2.10 }$$

For convenience, we denote $T_{\chi,g,R_0}(s)$ by $T_{R_0}(s)$. The constant $\chi >0$ should be large, resp. g  >  0 should be small, such that $-\frac{g}{2\chi}+\frac{1}{2\chi^2} < \frac{g}{2\chi}-\frac{1}{2\chi^2}$. The choice of these constants influences the 'sloppy'-ness of $T_{\chi,g,R_0}(s)$ in $[-\frac{g}{2\chi}-\frac{1}{2\chi^2}+R_0,\frac{g}{2\chi}+\frac{1}{2\chi^2}+R_0]$, see the figure below.

Our final aim is to minimize the free volume of the pore or to maximize the volume of deposited products before pore clogging. We do consider two optimization problems. In the first approach, we consider the minimization problem of the remaining pore volume at time $T, 0 < T < \infty$, where T is the time of clogging or a fixed approximation of it.

2.2.1. The first optimization problem.

For Problems I and II, the free volume $V({\rm cat}_0)$ of the pore at time T is defined by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{01.10**} V({\rm cat}_0)=\int_0^1{r^2({\rm cat}_0; T,x)\,{\rm d}x}. \nonumber \end{align} \tag{ 2.11 }$$

The first optimization problem aims at minimizing the free volume $V$ with respect to cat₀ over some admissible set.

2.2.2. The second optimization problem.

The second optimization problem aims at minimizing the mean value of free volume in the pore in a small interval of time [T₁,T] which is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{01.12*} \mathcal{V}({\rm cat}_0)={\int_{T_1}^{T}{\int_0^1{r^2({\rm cat}_0;t,x)\,{\rm d}x}}{\rm d}t}, \nonumber \end{align} \tag{ 2.12 }$$

for Problems I and II.

This problem avoids point evaluations in time and is somewhat easier to analyse.

2.2.3. Minimization by gradient descent.

One of the numerical schemes, which will be used in the section 5, will be a gradient descent method for the functional $\mathcal{V}$, i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {\rm cat}_0^{n+1}= {\rm cat}_0^{n} + \tau \partial {\mathcal{V}}({\rm cat}_0^n).\nonumber \end{align*}$$

This approach follows the classical strategy for parameter identification and optimization, as e.g. described in [6, 10–12, 15, 16]. The computation of $\partial {\mathcal{V}}({\rm cat}_0^n)$ requires to determine the Fréchet derivative of the parameter-to-state mapping $r({\rm cat}_0; \cdot , \cdot)$ with respect to cat₀ as well as the computation of its adjoint. This will be addressed in the following section.

Let $c=c(t,x) : Q_T\rightarrow \mathbb{R}$ be the oxygen concentration inside the pore and $r= r(t,x) : Q_T\rightarrow \mathbb{R}$ the pore radius. Instead of considering Problem I from (2.1a)–(2.1e) with $\beta$ and $\gamma$ defined by (2.2), we study a smoothed version by considering the system of PDEs given in (2.1a)–(2.1e) with the following definitions of $\beta$ and $\gamma$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.6*} \beta(t,r,c)= \frac{D_1}{\int_0^1{r \,H_{\epsilon}(c)}{\rm d}x},\quad \quad \gamma(t,r,c)=\frac{D_2}{\int_0^1{r H_{\epsilon}(c)}{\rm d}x}, \nonumber \end{align} \tag{ 3.1 }$$

with $\newcommand{\e}{{\rm e}} \epsilon$ being a small positive number and $\newcommand{\e}{{\rm e}} H_{\epsilon}(s)$ being the Huber function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:2:Huber} H_{\epsilon}(s)= \left\{\begin{array}{@{}l@{}} s, \quad {{\rm ~if~}} \, s \geqslant -\frac{a}{\chi}+c_0, \nonumber \\ 2\epsilon-\epsilon\,\Big[\frac{2\,\chi-1}{4\,\chi}-\frac{\chi \,s}{2} -\frac{\chi}2(\chi\,s-\chi\,c_0+a)\,(\chi\,s-\chi\,c_0+b) \nonumber \\ \quad \quad +\frac{\chi^3}{2}\,(\chi\,s-\chi\,c_0+a)^2(\chi\,s-\chi \,c_0+b)^2\Big], \quad {{\rm if}} \quad 0 < -\frac{b}{\chi}+c_0< s \leqslant -\frac{a}{\chi}+c_0, \nonumber \\ \epsilon, \, {{\rm otherwise}}, \end{array}\right. \nonumber \end{align} \tag{ 3.2 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:2:gamma} a:=1-\frac{1}{2\chi},\quad b:=1+\frac1{2\chi}, \quad \chi>0. \nonumber \end{align} \tag{ 3.3 }$$

The definition of the Huber function depends on $\newcommand{\e}{{\rm e}} \epsilon$ and $\chi$, see figure 1, where $\chi$ should be chosen as large as possible. In fact, the numerical computations in section 5 were done with $\chi= \infty$. Hence, for ease of notation we just call it $\newcommand{\e}{{\rm e}} H_{\epsilon}$.

Zoom In Zoom Out Reset image size

We see that $\beta(t,r,c)$ and $\gamma(t,r,c)$ are smoothed versions of (2.2). The constants $D_1, D_2$ can be made explicit by using (2.2) and (2.3).

Remark 3.1. The Huber function $\newcommand{\e}{{\rm e}} H_{\epsilon}$ is monotone, twice continuously differentiable, and $\newcommand{\bbR}{\mathbb{R}} \newcommand{\e}{{\rm e}} H_{\epsilon}(s)\geqslant \epsilon, \forall s \in \bbR$. The Huber function $\newcommand{\e}{{\rm e}} H_\epsilon$ is introduced to avoid singularities if c approaches 0. This process of Huberization has been used e.g. in [5] for smoothing step functions.

A mathematical analysis of this system of equations is more accessible after expanding $\frac{\partial}{\partial t}(r^2c)$ and after using a standard substitution for obtaining a homogenous system.

We suppose that r and c are sufficiently smooth and $0<r_\alpha\leqslant r, \forall (t,x)\in Q_T$, where $r_\alpha$ is some small positive constant. Then, expanding the derivatives in (2.1a), and dividing both sides by r², we get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.7*} c_t-c_{xx}-\frac{2r_x c_x}{r}+\frac{\Big(2r_t+\frac{D_1}{\int_0^1{r H_{\epsilon}(c)}{\rm d}x}\Big)c}{r}=0, \, \, (t,x)\in Q_T. \nonumber \end{align} \tag{ 3.4 }$$

To remove the non-homogeneous boundary condition (2.1e) we set

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.7.1} u(t,x)=c(t,x)-1. \nonumber \end{align} \tag{ 3.5 }$$

Then we have u(t,0)  =  u_x(t,1)  =  0 and equation (3.4) has the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.7**} u_t-u_{xx}-\frac{2r_x u_x}{r}+ \left(\frac{2r_t}{r}+\frac{D_1}{r\,\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}\right)u=-\left(\frac{2r_t}{r}+\frac{D_1}{r\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}\right), \, \, (t,x) \in Q_T. \nonumber \end{align} \tag{ 3.6 }$$

Inserting modified r_t from (2.1b) by $\newcommand{\e}{{\rm e}} r_t=-\frac{D_2}{\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}{\rm cat}_0(x)(u+1)$ into the right hand side of (3.6), and denoting u₀(x)  =  h(x)  −  1, we transform modified system for Problem I into its final form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\calL}{\mathcal{L}} \displaystyle u_t + \calL (u,r,{\rm cat}_0) u = f(u,r,{\rm cat}_0), \quad (t,x) \in Q_T, \label{2.8*} \nonumber \end{align} \tag{ 3.7 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_t=-\frac{D_2}{\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}{\rm cat}_0(x) (u+1),\, \, (t,x) \in Q_T,\label{2.10} \nonumber \end{align} \tag{ 3.8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(0,x)=r_0(x)=1, \, \,x\in \Omega,\label{2.10b} \nonumber \end{align} \tag{ 3.9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u(0,x)=u_0(x) , \, \, x \in \Omega,\label{2.11} \nonumber \end{align} \tag{ 3.10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u(t,0)=u_x(t,1)=0 , \, \,t \in [0,T]\label{2.12} \nonumber \end{align} \tag{ 3.11 }$$

with

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{L}(u,r,{\rm cat}_0(x))u= -u_{xx}-\frac{2r_x u_x}{r}+ \left(\frac{2r_t}{r}+\frac{D_1-2D_2\,{\rm cat}_0 }{r\,\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}\right)u,\nonumber \end{align*}$$

and

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle f(u,r,{\rm cat}_0)=\frac{2D_2{\rm cat}_0(x) -D_1}{r\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}. \nonumber \end{align*}$$

Similarly, we obtain a dynamic version of Problem II by replacing cat₀(x) with ${\rm cat}(t,x)$ and by adding a differential equation for cat as in (2.7) and (2.8).

3.1.1. Uniqueness and existence.

We use the standard Sobolev spaces $H^{\,p}(\Omega)$, $1 \leqslant p \leqslant \infty$ [7]. For a Banach space X with norm $\Vert \cdot\Vert_{X}$ we denote by L^p(0,T;X) the space consisting of all measurable functions $u: [0,T]\longrightarrow X$ with $\left\|u\right\|_{L^{\,p}{(0,T;X)}}:=\left(\int_0^T{\left\|u(t)\right\|^{\,p}_X {\rm d}t}\right){}^{\frac1p}<\infty$ for $1\leqslant p \leqslant \infty$, and $\left\|u\right\|_{L^{\infty}(0,T;X)}:={{\rm~esssup}}_{0\leqslant t \leqslant T}{\left\|u(t)\right\|_X}< \infty.$ We define

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle W(0,T;H^1(\Omega))=\{u: u\in L^2(0,T; H^1(\Omega)), u_t \in L^2(0,T; (H^{1}(\Omega))')\}.\nonumber \end{align*}$$

We sometimes also use the notation $W(0,T)$ for $W(0,T;H^1(\Omega))$.

Definition 3.2. A pair $(r,u)\in W^{1,\infty}(Q_T)\times W(0,T;H^1(\Omega))$ is said to be a weak solution of (3.7)–(3.11) if

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \forall v \in V \ {{\rm with~}} \ \ V := \Big\{v\in H^1(\Omega), v(0)=0 \Big\}\nonumber \end{align*}$$

the following equations are satisfied:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.12*} {(u_t,v)}_{\langle(H^{1}(\Omega))^*,H^1(\Omega)\rangle}+\int_0^1 u_xv_x {\rm d}x-\int_0^1 \frac{2r_x u_x}{r}v {\rm d}x + \int_0^1 \left(\frac{2r_t}{r}+\frac{D_1-2D_2{\rm cat}_0(x)}{r\, \int_0^1{r\,H_{\epsilon}(u+1)}{\rm d}x}\right)u\,v {\rm d}x \nonumber \\ =\int_0^1 \frac{2D_2{\rm cat}_0(x)-D_1}{r\int_0^1{r H_{\epsilon}(u+1)}{\rm d}x}\,v {\rm d}x,\qquad \forall v\in V \,\, {{\rm a.e.}}\, t \in [0,T], \nonumber \end{align} \tag{ 3.12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u(0,x)=u_0(x) , \, \, x \in [0,1], \nonumber \end{align} \tag{ 3.13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.15} r(t,x)=1- \int_0^t{\frac{D_2}{\int_0^1{{r} H_{\epsilon}(u+1)}{\rm d}x}{\rm cat}_0(x) \left(u(\tau,x)+1\right)}{\rm d}\tau, \, \, t \in [0,T]. \nonumber \end{align} \tag{ 3.14 }$$

To prove that there exists a weak solution to the system (3.7)–(3.11), we will apply the Banach fixed point theorem. The proof is long and tedious but follows the classical path of arguments for proving such results for non-linear partial differential equations, hence, we present only its main steps. We define the following spaces [17]

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\bfX}{\mathbf{X}} \displaystyle {\bfX_T}=\left\{r\in L^{\infty}(Q_T);r_x,\,r_t \in L^2(0,T;L^\infty(\Omega))\right\}, \nonumber \end{align*}$$

with norm $\newcommand{\bfX}{\mathbf{X}} \left\|r\right\|_{\bfX_T}=\left\|r\right\|_{L^{\infty}(Q_T)}+\left\|r_t\right\|_{L^2(0,T;L^\infty(\Omega))} + \left\|r_x\right\|_{L^2(0,T;L^\infty(\Omega))},$ and

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\bfY}{\mathbf{Y}} \displaystyle {\bfY_T}=\left\{u\in L^{2}(0,T;H^2(\Omega))\cap L^{\infty}(0,T;H^1(\Omega)), u_t\in L^{2}(0,T;L^2(\Omega))\right\}, \nonumber \end{align*}$$

with norm $\newcommand{\bfY}{\mathbf{Y}} \left\|u\right\|_{\bfY_T}=\left\|u\right\|_{L^{2}(0,T;H^2(\Omega))} + \left\|u\right\|_{L^{\infty}(0,T;H^1(\Omega))}+ \left\|u_t\right\|_{L^{2}(0,T;L^2(\Omega))}.$ Let $M, r_\alpha, R$ denote some fixed positive constants, we define the set

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \newcommand{\bM}{\mathbf{M}} \displaystyle \label{2.16} \bM_T=\left\{(r,u) \in \bfX_T \times \bfY_T, \left\|r-1\right\|_{\bfX_T}+\left\|u-u_0\right\|_{\bfY_T}\leqslant M; \,r_\alpha\leqslant r \leqslant R; \forall {{\rm a.e.~}} (t,x) \in Q_T\right\}, \nonumber \end{align} \tag{ 3.15 }$$

where $0 < r_\alpha < 1 < R.$ We see that $\newcommand{\bM}{\mathbf{M}} \bM_T$ is closed, convex and bounded in $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T \times \bfY_T$.

First, we prove the existence of a solution $\newcommand{\bM}{\mathbf{M}} (r,u) \in \bM_T$ for sufficiently small T. Before proving the contraction property need for applying Banach's theorem, we require a general result concerning the existence of solutions of partial differential equations. Relying on [7], we have the following result.

Theorem 3.3. Consider the initial/boundary-value problem

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\calL}{\mathcal{L}} \displaystyle \label{1.eq01} \left\{\begin{array}{@{}l@{}} u_t+\calL u =f,~{\rm in }~ Q_T \nonumber \\ u(t,0)=0,\, u_x(t,1)=0, \, t \in [0,T], \nonumber \\ u(0,x)=u_0(x)~ {\rm in }~\Omega, \end{array}\right. \nonumber \end{align} \tag{ 3.16 }$$

where $\newcommand{\bbR}{\mathbb{R}} f: Q \rightarrow \bbR$ and $\newcommand{\bbR}{\mathbb{R}} u_0: \Omega \rightarrow \bbR$ are given, and $u=u(t,x)$ is unknown, the operator $\newcommand{\calL}{\mathcal{L}} \calL$ denotes for each $t\in [0,T]$ the second-order non-divergence partial differential operator

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\calL}{\mathcal{L}} \displaystyle \label{eq:01:operator_L} \calL u= -a(t,x)u_{xx}+ b(t,x)u_x + d(t,x)u, \nonumber \end{align} \tag{ 3.17 }$$

for given coefficients $a(t,x), b(t,x), d(t,x)$ satisfying

• 1.  
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a(t,x)> \theta > 0 , {{\rm ~a.e.~}} (x,t)\in Q_T. \nonumber \end{align} \tag{ 3.18 }$$
• 2.  
  $a(t,x), a_t(t,x), b(t,x), d(t,x) \in L^{\infty}(Q_T),$
• 3.  
  $f\in L^2(0,T; L^2(\Omega))$.

Then, there exists unique weak solution u in $W(0,T;H^{1}(\Omega))$, which satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.aa} u\in L^2(0,T,H^2(\Omega))\cap L^{\infty}(0,T;V), \quad u_t\in L^2(0,T;L^2(\Omega)), \nonumber \end{align} \tag{ 3.19 }$$

and we have the estimate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.19} \mathrm{ess} \sup_{[0,T]} {\left\|u(t)\right\|_{H^{1}(\Omega)}} &+\left\|u\right\|_{L^2(0,T;H^2(\Omega))}+\left\|u_t\right\|_{L^2(0,T;L^2(\Omega))}\nonumber \\ &\leqslant C\left(\left\|f\right\|_{L^2(0,T;L^2(\Omega))}+\left\|u_0\right\|_{H^{1}(\Omega)}\right). \nonumber \end{align} \tag{ 3.20 }$$

For the proof of this theorem we refer the readers to section 1 in [17].

Due to the required regularity (3.19) this result is non-standard, at least we could not find a suitable result in the literature. The proof follows by an adaptation of the methods described in [7], for details see [17]. We now turn to prove the existence of a solution to the non-linear system of PDE's describing Problem I.

The functions u₀(x) and cat₀(x) are supposed to satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.7} u_0 \in W^{2,\infty}(\Omega), \quad {\rm cat}_0\in W^{1,\infty}(\Omega). \nonumber \end{align} \tag{ 3.21 }$$

Theorem 3.4. For T  >  0 sufficiently small and under condition (3.21), there exists a unique weak solution $(u,r)$ on $\newcommand{\bM}{\mathbf{M}} \bM_T$ of (3.7)–(3.11).

Sketch of proof. As usual the existence proof for a non-linear PDE is turned into a fixed point problem involving a linear differential equation. We define the operator $\mathcal{T}$

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\calT}{\mathcal{T}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \newcommand{\bM}{\mathbf{M}} \displaystyle \calT: \bM_T \longrightarrow \bfX_T\times \bfY_T\nonumber \end{align*}$$

which maps $\newcommand{\bM}{\mathbf{M}} (\bar r, \bar u) \in \bM_T$ onto $\newcommand{\calT}{\mathcal{T}} (r,u)=\calT(\bar{r}, \bar{u})$ with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.18} r(t,x)=1-\int_0^t{\frac{D_2}{\int_0^1{\bar{r} H_{\epsilon}(\bar{u}+1)}{\rm d}x}{\rm cat}_0(x)(\bar{u}+1)}{\rm d}\tau \nonumber \end{align} \tag{ 3.22 }$$

and u is the solution to the linear parabolic equation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\calL}{\mathcal{L}} \displaystyle \label{2.13} u_t+\calL(\bar u,r,{\rm cat}_0) u=f(\bar u,r,{\rm cat}_0), \nonumber \end{align} \tag{ 3.23 }$$

with initial and boundary conditions given by (3.10) and (3.11).

For proving the above stated theorem, we need the following two properties:

• 1.  
  $\newcommand{\calT}{\mathcal{T}} \calT$ is defined on $\newcommand{\bM}{\mathbf{M}} \bM_T$ and maps $\newcommand{\bM}{\mathbf{M}} \bM_T$ into $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T \times \bfY_T$ and there exists T small enough such that $\newcommand{\calT}{\mathcal{T}} \calT$ maps $\newcommand{\bM}{\mathbf{M}} \bM_T$ onto $\newcommand{\bM}{\mathbf{M}} \bM_T$.
• 2.  
  There exists T small enough such that $\newcommand{\calT}{\mathcal{T}} \calT$ is s-contractive in $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T\times \bfY_T$, i.e.

  $$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\calT}{\mathcal{T}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \Vert \calT(u_1,r_1)- \calT(u_2,r_2)\Vert_{\bfX_T \times \bfY_T} \leqslant s\,\Vert (u_1,r_1)- (u_2,r_2)\Vert_{\bfX_T \times \bfY_T},\nonumber \end{align*}$$

  for all $\newcommand{\bM}{\mathbf{M}} (u_1,r_1),\,(u_2,r_2) \in \bM_T$, with $s \in [0,1).$

This approach, as well as the proofs, follow from the classical outline for non-linear PDEs given in [7].

Proof of (1). The tricky part is an appropriate definition of $\newcommand{\bM}{\mathbf{M}} \bM_T$ and of the function spaces used, the results then follow directly or by applying Hölder's inequality several times. For readers interested in the details of the different steps, we refer to [17]. We directly obtain, that $\newcommand{\bM}{\mathbf{M}} (\bar{r}, \bar{u})\in \bM_T$ and assuming that cat₀ is fixed and satisfies (3.21), implies

• (i)  
  there exists a constant $C(M,T,R,r_\alpha,{\rm cat}_0)$ depending on $M, T, R, r_\alpha, {\rm cat}_0$ such that $\left\|r_t\right\|_{L^\infty(Q_T} \leqslant C(M,T,R,r_\alpha,{\rm cat}_0)$ and $\left\|r_x\right\|_{L^\infty(Q_T)}\leqslant C(M,T,R,r_\alpha,{\rm cat}_0)$
• (ii)  
  There exists a positive constant T₁ such that $r>r_\alpha, \forall (t,x) \in [0,T_1]\times [0,1] := Q_{T_1}$ and the coefficients of the operator $\newcommand{\calL}{\mathcal{L}} \calL$ in (3.23) belong to $L^\infty(Q_{T_1})$. Based on theorem 3.3, there exists a unique solution u of (3.23) satisfying $u\in W(0,T_1)$ and by using Hölder's inequality we obtain the following improved regularity result

  $$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle u\in L^2(0,T_1,H^2(\Omega))\cap L^{\infty}(0,T_1;H^1(\Omega)), \quad u_t\in L^2(0,T_1;L^2(\Omega)),\nonumber \end{align*}$$

  and

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ess}{\mathrm{ess}} \displaystyle \ess\sup_{[0,T_1]} {\left\|u(t)\right\|_{H^1(\Omega)}} &+\left\|u\right\|_{L^2(0,T_1;H^2(\Omega))}+\left\|u_t\right\|_{L^2(0,T_1;L^2(\Omega))}\nonumber \\ &\leqslant C(M,T_1,R,r_\alpha,{\rm cat}_0)\left(\left\|f\right\|_{L^2(0,T_1;L^2(\Omega))} + \left\|u_0\right\|_{H^1(\Omega)}\right). \nonumber \end{align} \tag{ 3.24 }$$

  For convenience, we suppose that $T_1\leqslant 1$, which can be achieved by rescaling the time. From the above results we get the following lemma.

Lemma 3.5 Assume that $\newcommand{\bM}{\mathbf{M}} (\bar{r}, \bar{u})\in \bM_{T_1}$. Then, the weak solution u of the problem (3.23), (3.10) and (3.11) satisfies the estimate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.21} &\mathrm{ess} \sup_{[0,T_1]} {\left\|u-u_0\right\|_{H^1(\Omega)}}+\left\|u-u_0\right\|_{L^2(0,T_1;H^2(\Omega))}+\left\|(u-u_0)_t\right\|_{L^2(0,T_1;L^2(\Omega))} \nonumber \\ &\leqslant C(M,T_1,R,r_\alpha,u_0,{\rm cat}_0)\sqrt{T_1}. \nonumber \end{align} \tag{ 3.25 }$$

It follows from this lemma and the continuous embedding $H^1(\Omega)\subset L^\infty(\Omega)$ that the weak solution u of problem (3.23), (3.10) and (3.11) satisfies the estimates

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Vert u\Vert_{L^\infty(Q_{T_1})}=\mathrm{ess}\, {\rm sup}_{t,x}{\left|u\right|}\leqslant C(M,T_1,R,r_\alpha, u_0, {\rm cat}_0),\label{2.33*} \nonumber \end{align} \tag{ 3.26 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Vert u_x\Vert_{L^2(0,T;L^\infty(\Omega))}\leqslant C(M,T_1,R,r_\alpha, u_0, {\rm cat}_0),\label{2.32*} \nonumber \end{align} \tag{ 3.27 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Vert u_x\Vert_{L^\infty(0,T;L^2(\Omega))}\leqslant C(M,T_1,R,r_\alpha, u_0, {\rm cat}_0).\label{2.32**} \nonumber \end{align} \tag{ 3.28 }$$

All constants depend continuously on T and tend to zero for $T \rightarrow 0$, hence for T sufficiently small we have that $\newcommand{\bM}{\mathbf{M}} \newcommand{\calT}{\mathcal{T}} \calT: \bM_T \rightarrow \bM_T$.

Proof of (2). Now we show that $\mathcal{T}$ is contractive in $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T \times \bfY_T$ norm with T small enough. Let $(\bar{u_1},\bar{r_1})$, $\newcommand{\bM}{\mathbf{M}} (\bar{u_2},\bar{r_2})\in \bM_{T^*}$, where T^*  >  0 is such that $\newcommand{\bM}{\mathbf{M}} \newcommand{\calT}{\mathcal{T}} \calT: \bM_{T^*} \rightarrow \bM_{T^*}$. is valid. Suppose that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.55} \mathcal{T}(\bar{u}_1,\bar{r}_1)=(u_1,r_1),\quad \mathcal{T}(\bar{u}_2,\bar{r}_2)=(u_2,r_2). \nonumber \end{align} \tag{ 3.29 }$$

Then $\bar{u}_1,\bar{r}_1,u_1,r_1$ and $\bar{u}_2,\bar{r}_2,u_2,r_2$ satisfy equations (3.7)–(3.11). Denote $\tilde{u}=u_1-u_2, \tilde{r}=r_1-r_2$. Subtracting the two parabolic equations and the two differential equations, we get the system:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.56} &\tilde{u}_{t}-\tilde{u}_{xx}- \frac{r_{1x}}{r_1}\tilde{u}_{x} +\left(\frac{2r_{1t}}{r_1} + \frac{D_1-2D_2{\rm cat}_0(x)}{r_1\,\int_0^1{r_{1}\,H_{\epsilon}(\bar{u}_1+1)}{\rm d}x}\right)\tilde{u} \nonumber \\ &=\left(\frac{2r_{1x}}{r_{1}}-\frac{2r_{2x}}{r_{2}}\right)u_{2x} -\left(\frac{2r_{1t}}{r_{1}}-\frac{2r_{2t}}{r_{2}}\right)u_{2}\nonumber \\ &\quad +(2D_2{\rm cat}_0(x)-D_1)\left(\frac{1}{r_1\,\int_0^1{r_{1}\,H_{\epsilon}(\bar{u}_1+1)}{\rm d}x} -\frac{1}{r_2\,\int_0^1{r_{2}\,H_{\epsilon}(\bar{u}_2+1)}{\rm d}x}\right)(u_2+1)\nonumber \end{align} \tag{ 3.30 }$$

with the initial and boundary conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{u}(0,x)=0, \quad x \in \Omega \label{2.57} \nonumber \end{align} \tag{ 3.31 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{u}(t,0)=\tilde{u}_x(t,1)=0, \quad t \in [0,T^*],\label{2.58} \nonumber \end{align} \tag{ 3.32 }$$

coupled with the equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{r}(x,t) &= -D_2{\rm cat}_0(x) \Big[\int_0^t\frac{\bar{u}_1-\bar{u}_2}{\int_0^1{\bar{r}_{1}\,H_{\epsilon}(\bar{u}_1+1)}{\rm d}x}{\rm d}\tau \nonumber \\& + \int_0^t(\bar{u}_2+1)\Big(\frac{1}{\int_0^1{\bar{r}_{1}\,H_{\epsilon}(\bar{u}_1+1)}{\rm d}x} -\frac{1}{\int_0^1{\bar{r}_{2}\,H_{\epsilon}(\bar{u}_2+1)}{\rm d}x}\Big){\rm d}\tau \Big]. \label{2.59} \nonumber \end{align} \tag{ 3.33 }$$

Using the same estimates as in the proof of (1.) we obtain the following estimate for the problem (3.30)–(3.32) on $Q_{T^*}=[0,T^*]\times [0,1]$ with T^* sufficiently small

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{2.61} &\mathrm{ess} \sup_{[0,T^*]} {\big\|\tilde{u}(t)\big\|_{H^1(\Omega)}} +\big\|\tilde{u}\big\|_{L^2(0,T^*;H^2(\Omega))}+\big\|\tilde{u}_t\big\|_{L^2(0,T^*;L^2(\Omega))} \nonumber \\ &\leqslant C(M,T^*,\epsilon,r_\alpha, u_0, {\rm cat}_0)\, \sqrt{T^*} \big( \Vert \bar{r}_1-\bar{r}_2\Vert_{L^\infty(Q_{T^*})} \nonumber \\ &+ \Vert\bar{u}_{1}-\bar{u}_{2}\Vert_{L^2(0,T^*; H^2(\Omega))}+\Vert\bar{u}_{1}-\bar{u}_{2}\Vert_{L^\infty(0,T^*; H^1(\Omega))}\big)\nonumber \end{align} \tag{ 3.34 }$$

and the functions $r_1, r_2$ defined in (3.29) satisfy the estimate:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \label{2.70} \left\|r_1-r_2\right\|_{\bfX}\leqslant T^* C(M,R,{\rm cat}_0) \left(\left\|\bar{r}_1-\bar{r}_2\right\|_{L^\infty(Q_{T^*})}+\left\|\bar{u}_1-\bar{u}_2\right\|_{L^\infty(Q_{T^*})}\right). \nonumber \end{align} \tag{ 3.35 }$$

Combining these two estimates, we see that there exists a positive number T^** such that $\newcommand{\calT}{\mathcal{T}} \calT$ is s-contractive in $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_{T} \times \bfY_{T}$ if $T\leqslant T^{**}$, where s is a positive constant less than 1.

This ends the sketch of the proof of theorem 3.4.

From now on, we tacitly assume that T is small enough so that there exists a unique solution to the system (3.7)–(3.11).

From lemma 3.5 we also obtain the continuous dependency of the weak solution in $\newcommand{\bM}{\mathbf{M}} \bM_T$ on the initial data [17].

Theorem 3.6. Let $(r,u)$ be a weak solution of the system (3.7)–(3.11). Then it continuously depends on the data

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \label{2.89} \left\|u\right\|_{\bfY_T}+\Vert r\Vert_{\bfX}\leqslant C(M,T,R,r_\alpha,{\rm cat}_0)\sqrt{T}\Vert u_0\Vert_{H^2(\Omega)}. \nonumber \end{align} \tag{ 3.36 }$$

3.1.2. Continuity and differentiability of the parameter-to-state mapping.

Now we analyse the parameter-to-state map, i.e. we consider $r=r({\rm cat}_0)$ as well as $u=u({\rm cat}_0)$ and show the continuity of $(r,u)$ in $\newcommand{\bfY}{\mathbf{Y}} W^{1,\infty}(Q_T)\times \bfY_T$ with respect to cat$_0(x) \in W^{1,\infty}(\Omega)$ for Problem I.

Let cat$_{01}(x), {\rm cat}_{02}(x)\in W^{1,\infty}(\Omega)$, and assume that there exists a positive constant $C_{{\rm cat}}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3:cat} \Vert {\rm cat}_{01}(x)\Vert_{W^{1,\infty}(\Omega)}\leqslant C_{{\rm cat}},\quad \Vert {\rm cat}_{02}(x)\Vert_{W^{1,\infty}(\Omega)}\leqslant C_{{\rm cat}}. \nonumber \end{align} \tag{ 3.37 }$$

Suppose that $\newcommand{\bM}{\mathbf{M}} (u_1,r_1), (u_2,r_2) \in \bM_T$, where $T\leqslant T^{**}$, are solutions to (3.7)–(3.11). Denoting $\newcommand{\calU}{\mathcal{U}} \newcommand{\calR}{\mathcal{R}} \calU=u_1-u_2,\, \calR=r_1-r_2$, we have the following result.

Theorem 3.7. There exists T small enough and a positive constant $C(M,T,R,r_\alpha, u_0, C_{{\rm cat}_0})$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\calR}{\mathcal{R}} \newcommand{\calU}{\mathcal{U}} \newcommand{\bfY}{\mathbf{Y}} \displaystyle \label{eq:3:12*} \Vert \calU\Vert_{\bfY_T}+\Vert \calR\Vert_{W^{1,\infty}(Q_T)}\leqslant C(M,T,R,r_\alpha, u_0, C_{{\rm cat}_0})\,\Vert {\rm cat}_{01}-{\rm cat}_{02}\Vert_{W^{1,\infty}(\Omega)}, \nonumber \end{align} \tag{ 3.38 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\calR}{\mathcal{R}} \displaystyle \Vert \calR\Vert_{W^{1,\infty}(Q_T)}=\Vert \calR\Vert_{L^\infty(Q_T)}+\Vert \calR_t\Vert_{L^\infty(Q_T)}+\Vert \calR_x\Vert_{L^\infty(Q_T)},\nonumber \\ \Vert {\rm cat}_{01}-{\rm cat}_{02}\Vert_{W^{1,\infty}(\Omega)} = \Vert {\rm cat}_{01}-{\rm cat}_{02}\Vert_{L^\infty(\Omega)}+\Vert {\rm cat}_{01x}-{\rm cat}_{02x}\Vert_{L^\infty(\Omega)}. \nonumber \end{align*}$$

Due to the continuous embedding of $W^{1,\infty}(Q_T)$ in $\newcommand{\bfX}{\mathbf{X}} \bfX_T$, from theorem 3.7, we have

Corollary 3.8. The pair $(r,c)$ is continuous on $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T \times \bfY_T$ with respect to ${\rm cat}(x)$ in $W^{1,\infty}(\Omega)$ for problem I with T small enough.

Now we investigate the Fréchet differentiability of r and u with respect to cat₀(x) for Problem I. Under the assumptions of theorem 3.4 we define the state mapping $\newcommand{\calP}{\mathcal{P}} \calP$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\calP}{\mathcal{P}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \calP : W^{1,\infty}(\Omega) \longrightarrow \bfX_T\times\bfY_T, \quad {\rm cat}_0 \longmapsto \calP({\rm cat}_0)=\left( r({\rm cat}_0), u({\rm cat}_0)\right). \nonumber \end{align} \tag{ 3.39 }$$

For $\delta \in W^{1,\infty}(\Omega)$ we denote $\frac{\partial r}{\partial {\rm cat}_0}({\rm cat}_0)\delta$ by $\newcommand{\e}{{\rm e}} \eta(t,x)$ and $\frac{\partial c}{\partial {\rm cat}_0}({\rm cat}_0)\delta$ by $\varphi(t,x)$.

By the usual approach, we use the differential equations for $r({\rm cat}_0),u({\rm cat}_0)$ and for $r({\rm cat}_0+\delta),u({\rm cat}_0 +\delta)$ and obtain after expanding all terms and deleting higher order terms the so-called sensitivity system:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \eta_t=-\frac{D_2\,{\rm cat} \varphi}{\int_0^1 r\,H_\epsilon(c) {\rm d}x}+\frac{D_2 {\rm cat} c}{(\int_0^1 r H_\epsilon(c))^2}\int_0^1 r{H_\epsilon'}(c)\varphi+H_\epsilon(c)\eta {\rm d}x- \frac{D_2\delta\,c}{\int_0^1 r\,H_\epsilon(c) {\rm d}x},\label{2.88*} \nonumber \end{align} \tag{ 3.40 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\varphi_t-\varphi_{xx}-\frac{2\,r_x\varphi_x}{r}+\left(\frac{2r_t}{r}+\frac{D_1}{r\int_0^1 r H_\epsilon(c) {\rm d}x}-\frac{D_2\,{\rm cat}}{\int_0^1 r\,H_\epsilon(c) {\rm d}x}\right)\varphi \nonumber \\ &=\left(\frac{2 D_2 {\rm cat} c^2}{r(\int_0^1 r H_\epsilon(c))^2}-\frac{D_1\,c}{r(\int_0^1 r H_\epsilon(c) {\rm d}x)^2}\right)\int_0^1 H_\epsilon(c)\eta+r{H_\epsilon'}(c)\varphi {\rm d}x \nonumber \\ &\quad +\frac{2c_x\,\eta_x}{r} -\left(\frac{2r_x\,c_x}{r^2}-\frac{2r_t\,c}{r^2}- \frac{D_1\,c}{r^2\,\int_0^1 r H_\epsilon(c) {\rm d}x}\right)\eta- \frac{2\,\delta\,c^2}{r\int_0^1 r\,H_\epsilon(c) {\rm d}x},\label{2.99} \nonumber \end{align} \tag{ 3.41 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \eta(0,x)=0, \, \, x \in \Omega,\label{subeq:3:1} \nonumber \end{align} \tag{ 3.42 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varphi(0,x)=0, \, \, x \in \Omega,\label{subeq:3:2} \nonumber \end{align} \tag{ 3.43 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varphi(t,0)=0, \, \, \varphi_x(t,1)=0, \, \, t \in [0,T].\label{subeq:3:4} \nonumber \end{align} \tag{ 3.44 }$$

We then have to show that a solution to this sensitivity system exists and that its weak solution is indeed the desired Fréchet derivative.

First of all, we see that the sensitivity system (3.40)–(3.44) is linear in $\newcommand{\e}{{\rm e}} (\eta, \varphi)$ and the right-hand side if of order ${{\mathcal O}}(\delta)$. Hence standard arguments, e.g. adapting the fix-point techniques of the previous subsection to this system on linear equations, yield the existence of a weak solution of this sensitivity system.

More precisely, for some given positive constant N, we define the set

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \newcommand{\bN}{\mathbf{N}} \displaystyle \label{2.101} \bN_T=\left\{(\eta,\varphi) \in \bfX_T \times \bfY_T, \left\|\eta\right\|_{\bfX_T}+\left\|\varphi\right\|_{\bfY_T}\leqslant N\right\}. \nonumber \end{align} \tag{ 3.45 }$$

Lemma 3.9. For T small enough, there exists a unique solution $\newcommand{\bN}{\mathbf{N}} \newcommand{\e}{{\rm e}} (\eta, \varphi) \in \bN_T$ of the sensitivity system (3.40)–(3.44) and $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\e}{{\rm e}} \|(\eta,\varphi) \|_{L_\infty(\bfX_T \times \bfY_T)}={{\mathcal O}}(\delta)$.

We now prove, that the parameter-to-state map is indeed Fréchet-differentiable with respect to cat₀.

Theorem 3.10. There exists T small enough such that $\newcommand{\bfX}{\mathbf{X}} r({\rm cat})\in\bfX_T$ and $\newcommand{\bfY}{\mathbf{Y}} u({\rm cat}_0)\in\bfY_T$ are Fréchet differentiable with respect to cat$_0\in W^{1,\infty}(\Omega)$ in the sense that

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \lim_{\Vert\delta\Vert_{W^{1,\infty}(\Omega)}\rightarrow 0}{\frac{\Vert r({\rm cat}_0+\delta)-r({\rm cat}_0)-\eta\Vert_{\bfX_T}}{\Vert\delta\Vert_{W^{1,\infty}(\Omega)}}}=0 ~{\rm and } ~ \lim_{\Vert\delta\Vert_{W^{1,\infty}(\Omega)}\rightarrow 0}{\frac{\Vert u({\rm cat}_0+\delta)-u({\rm cat}_0)-\varphi\Vert_{\bfY_T}}{\Vert\delta\Vert_{W^{1,\infty}(\Omega)}}}=0.\nonumber \end{align*}$$

To prove the theorem we need the following result.

Lemma 3.11. Denote

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \xi=r({\rm cat}_0+\delta)-r({\rm cat}_0)-\eta ~{\rm and }~ \psi=u({\rm cat}_0+\delta)-u({\rm cat}_0)-\varphi. \nonumber \end{align*}$$

Then, $\xi$ and $\psi$ satisfy the problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi_t-\psi_{xx}- \frac{2r_x\psi_x}{r({\rm cat}_0+\delta)}+F_1(\xi,\psi)=F_2(\eta,\varphi)+\omega_\delta, \, \, (t,x)\in Q_T, \label{subeq:2:xi_psi.1} \nonumber \end{align} \tag{ 3.46a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi_t=F_3(\xi,\psi)+F_4(\eta,\varphi)+\omega_\delta, \quad (t,x)\in Q_T,\label{subeq:2:xi_psi.2} \nonumber \end{align} \tag{ 3.46b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi(0,x)=0, \, \, x\in \Omega,\label{subeq:2:xi_psi.3} \nonumber \end{align} \tag{ 3.46c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi(0,x)=0, \quad x \in \Omega,\label{subeq:2:xi_psi.4} \nonumber \end{align} \tag{ 3.46d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi(t,0)=0, \, \, t\in [0,T],\label{subeq:2:xi_psi.5} \nonumber \end{align} \tag{ 3.46e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi_x(t,0)=0, \quad t\in [0,T],\label{subeq:2:xi_psi.6} \nonumber \end{align} \tag{ 3.46f  }$$

where $\newcommand{\e}{{\rm e}} F_2(\eta,\varphi)$, $\newcommand{\e}{{\rm e}} F_4(\eta,\varphi)$ contain only the functions in two variables $\newcommand{\e}{{\rm e}} \eta$ and $\varphi$, in which the lowest-degree term is of the second degree and $F_1(\xi,\psi)$, $F_3(\xi,\psi)$ are linear with respect to $\xi$ and $\psi$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3:omega_delta_1} \Vert\omega_\delta\Vert_{L^\infty(Q_T)}=o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)}). \nonumber \end{align} \tag{ 3.47 }$$

Proof. For proving the lemma we use the Taylor expansion of $\newcommand{\e}{{\rm e}} H_\epsilon$ and theorem 3.7 to obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:3:H_taylor_expanssion_2} H_\epsilon(u({\rm cat}_0+\delta))-H_\epsilon(u)={H_\epsilon'}(u)(u({\rm cat}_0+\delta)-u)+\frac12{H_\epsilon''}(u)(u({\rm cat}_0+\delta)-u)+\omega_\delta, \nonumber \end{align} \tag{ 3.48 }$$

where $\omega_\delta$ satisfies (3.47).

We need to expand all expression appearing in the differential equations of $\xi$ and $\psi$. A somewhat longer calculation and ordering terms yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.106*} {{\mathcal I}}_1=\frac1{\int_0^1 r\,H_\epsilon(c){\rm d}x}-\frac1{\int_0^1 r({\rm cat}+\delta) H_\epsilon(c({\rm cat}+\delta)){\rm d}x}= I_1(\xi,\psi) + J_1(\eta,\varphi), \nonumber \end{align} \tag{ 3.49 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.106**} I_{1}(\xi,\psi)=&\frac{\int_0^1 \xi H_\epsilon(u({\rm cat}_0+\delta)){\rm d}x} {\int_0^1 r({\rm cat}_0+\delta) H_\epsilon(u({\rm cat}_0+\delta)){\rm d}x \int_0^1 r H_\epsilon(u({\rm cat}_0+\delta)){\rm d}x} \nonumber \\ & +\frac{\int_0^1 r\,{H'}(u)\psi {\rm d}x}{\int_0^1 r H_\epsilon(u({\rm cat}_0+\delta)){\rm d}x \, \int_0^1 r H_\epsilon(u){\rm d}x} \nonumber \\ &+\frac{1}{2} \frac{\int_0^1 r\,{H''}(u)(u({\rm cat}_0+\delta)-u+\varphi)\psi {\rm d}x}{\int_0^1 r H_\epsilon(u({\rm cat}_0+\delta)){\rm d}x \, \int_0^1 r H_\epsilon(u){\rm d}x}, \nonumber \end{align} \tag{ 3.50 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.107**} J_{1}(\eta,\varphi)=&\frac{\int_0^1 \eta H_\epsilon(c({\rm cat}+\delta)){\rm d}x}{\int_0^1 r({\rm cat}+\delta) H_\epsilon(c({\rm cat}+\delta)){\rm d}x\,\int_0^1 r\, H_\epsilon(c({\rm cat}+\delta)){\rm d}x} \nonumber \\ & + \frac{\int_0^1 r\,{H'}(c)\varphi {\rm d}x}{\int_0^1 r H_\epsilon(c({\rm cat}+\delta)){\rm d}x \, \int_0^1 r H_\epsilon(c){\rm d}x} \nonumber \\ & + \frac12 \, \frac{\int_0^1 r\,{H''}(c*)\varphi^2 {\rm d}x}{\int_0^1 r H_\epsilon(c({\rm cat}+\delta)){\rm d}x}. \nonumber \end{align} \tag{ 3.51 }$$

We see that $I_1(\xi,\psi)$ is linear with respect to $\xi$ and $\psi$, and $\newcommand{\e}{{\rm e}} J_1(\eta,\varphi)$ contains only expressions of second degree in $\newcommand{\e}{{\rm e}} \eta$ and $\varphi$.

Similarly, we expand

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\mathcal I}}_2 &=\frac1{\int_0^1 r\,H_\epsilon(c){\rm d}x}-\frac{1}{\int_0^1 r({\rm cat}+\delta) H_\epsilon(c({\rm cat}+\delta)){\rm d}x}-\frac{\int_0^1 r\,{H_\epsilon'}(c)\varphi +H_\epsilon(c)\eta {\rm d}x}{(\int_0^1 r\,H_\epsilon(c){\rm d}x)^2} \nonumber \\ & =I_2(\xi,\psi) + J_2(\eta,\varphi)+ \overline{\omega}_\delta \nonumber \end{align} \tag{ 3.52 }$$

where $I_2(\xi,\psi)$ is linear in $\xi$ and $\psi$, $\newcommand{\e}{{\rm e}} J_2(\eta,\varphi)$ is a second order polynomial in $\newcommand{\e}{{\rm e}} \eta$ and $\varphi$, and $\Vert\overline{\omega}_\delta\Vert_{L^\infty(Q_T)}=o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)}).$ We use the differential equations for $r({\rm cat}_0+\delta)$, $r({\rm cat}_0)$ and $\newcommand{\e}{{\rm e}} \eta$ and obtain after some algebraic manipulations a differential equation for $\xi_t$ of the form

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \xi_t= F_3(\xi,\psi) + F_4(\eta,\varphi) + \tilde{\omega}_\delta,\nonumber \end{align*}$$

where $F_3(\xi,\psi)$ is linear in $\xi$ and $\psi$,and $\newcommand{\e}{{\rm e}} F_4(\eta,\varphi)$ is a second order polynomial in $\newcommand{\e}{{\rm e}} \eta$ and $\varphi$, and $\Vert\tilde{\omega}_\delta\Vert_{L^\infty(Q_T)}=o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)}).$ Similarly, we use the differential equations for $u({\rm cat}_0+\delta)$, $u({\rm cat}_0)$ and $\varphi$ and obtain—after expanding and sorting all terms—equation (3.46a). For the details of the expansion see p 54 in [17]. □

Hence, $\xi$ and $\psi$ satisfy linear differential equations and by using theorem 3.7 we observe that the right hand side of the system of equations is of order $o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)})$. Hence, by standard arguments [7], the weak solution to this problem is unique and the limits in theorem 3.10 hold.

The solution r (the pore radius) depends on the catalyst position cat₀(x), time t and space variable x: $r := r({\rm cat}_0;t,x)$. Our aim is to minimize the free volume (2.11) of the pore after some time T, with respect to cat₀ over some admissible set $\newcommand{\bA}{\boldsymbol{A}} \bA$.

As $\newcommand{\bfX}{\mathbf{X}} r\in \bfX_T \subset C([0,T]; L^2(\Omega))$, we can define the linear bounded projection operator

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \mathbf{P}_T: \bfX_T \longrightarrow L^2(\Omega), \quad \mathbf{P}_T r({\rm cat}_0;\cdot,\cdot):= r({\rm cat}_0;T,\cdot). \nonumber \end{align} \tag{ 4.1 }$$

Hence

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{2.83.1} V({\rm cat})=\int_0^1 r^2({\rm cat}_0;T,x) {\rm d}x= \int_0^1{(\mathbf{P}_T r({\rm cat}_0,\cdot,x))^2 }{\rm d}x. \nonumber \end{align} \tag{ 4.2 }$$

In section 3 we have shown that $r({\rm cat}_0;t,x)$ is Fréchet differentiable with respect to cat$_0\in W^{1,\infty}(\Omega)$ for the smoothed version of Problem I, similarly we assume that there exists a linear bounded operator

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfX}{\mathbf{X}} \displaystyle r_{{\rm cat}_0} :W^{1,\infty}(\Omega) \longrightarrow L(W^{1,\infty}(\Omega),\bfX_T), \quad {\rm cat}_0 \longmapsto r_{{\rm cat}_0}({\rm cat}_0), \nonumber \end{align} \tag{ 4.3 }$$

such that, for $\delta \in W^{1,\infty}(\Omega)$,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \label{eq:4:4*} \Vert r({\rm cat}_0+\delta;\cdot,\cdot)-r({\rm cat})-r_{{\rm cat}_0}({\rm cat}_0;\cdot,\cdot)\delta\Vert_{\bfX_T}=o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)}), \nonumber \end{align} \tag{ 4.4 }$$

which implies that

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle V({\rm cat}_0+\delta)-V({\rm cat}_0) = 2\langle \mathbf{P}_T r_{{\rm cat}_0}\delta, \mathbf{P}_Tr({\rm cat}_0,\cdot,\cdot)\rangle_{L^2(\Omega)}+o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)}) \nonumber \end{align*}$$

and from that we directly obtain that the functional $V({\rm cat}_0)$ is Fréchet differentiable with respect to cat$_0\in W^{1,\infty}(\Omega)$.

Remark 4.1. For convenience, we will use the notation $r_{{\rm cat}_0}(\cdot)$ for $r_{{\rm cat}_0}({\rm cat}_0)$.

Although $r({\rm cat}_0)$ and $c({\rm cat}_0)$ are Fréchet differentiable with respect to cat₀ in $L^\infty$-norm, however, for computing their adjoint operators we have to suppose that they are Fréchet differentiable with respect to cat₀ in L²-norm, that means $r_{{\rm cat}_0}(\cdot)$ and $c_{{\rm cat}_0}(\cdot)$ are linear bounded operators

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4:assume_oper_1} r_{{\rm cat}_0}(\cdot): L^2(\Omega) \mapsto L^2(Q_T), \quad \delta \mapsto r_{{\rm cat}_0}(\cdot)\delta, \nonumber \end{align} \tag{ 4.5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4:assume_oper_2} c_{{\rm cat}_0}(\cdot): L^2(\Omega) \mapsto L^2(Q_T), \quad \delta \mapsto c_{{\rm cat}_0}(\cdot)\delta. \nonumber \end{align} \tag{ 4.6 }$$

Further, we suppose that the projection operator $\mathbf{P}_{T}$ is bounded from $L^2(Q_T)$ onto $L^2(\Omega)$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4:assume_oper_3} \mathbf{P}_{T}: L^2(Q_T) \longrightarrow L^2(\Omega), \quad r({\rm cat};t,x) \longmapsto r({\rm cat};T,x). \nonumber \end{align} \tag{ 4.7 }$$

From these assumptions, by the Riesz representation theorem, we see that $V_{{\rm cat}_0}(\cdot)$ is identified with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 2 (r_{{\rm cat}_0}(\cdot))^*(\mathbf{P}_T)^*\mathbf{P}_T r({\rm cat}_0,\cdot,\cdot). \label{V_cat0_adjoint} \nonumber \end{align} \tag{ 4.8 }$$

A similar result can be obtained for Problem II. In the next paragraph we will present a scheme how to calculate $(r_{{\rm cat}_0}(\cdot)){}^*$.

To prove the existence of a minimizer of the functional (4.2) we introduce the admissible set $\newcommand{\bA}{\boldsymbol{A}} \bA$ (see [9, 20]) as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bA}{\boldsymbol{A}} \displaystyle \label{eq:04:admission_set_A} \bA=\left\{{\rm cat} \in C^{1,\lambda}(\Omega), 0\leqslant {\rm cat} \leqslant 1, M_1 \leqslant {\rm cat}_x \leqslant M_2, \forall x\in [0,1],\sup_{x_1,x_2 \in [0,1]}{\frac{\vert {\rm cat}_x(x_1)-{\rm cat}_x(x_2)\vert}{\vert x_1-x_2\vert^\lambda}}\leqslant C \right\}, \nonumber \end{align} \tag{ 4.9 }$$

where $M_1, M_2$ are some given positive numbers and $0<\lambda<1$. The set $\newcommand{\bA}{\boldsymbol{A}} \bA$ is compact in $C^1(\Omega)$ [20]. Since $V({\rm cat}_0)$ is continuous in $W^{1,\infty}(\Omega)$, we have the following result.

Theorem 4.2. The problem of minimizing $V({\rm cat}_0)$ subject to Problem I over $\newcommand{\bA}{\boldsymbol{A}} \bA$ admits at least one solution.

We consider the second optimization problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \min_{{\rm cat}_0\in \mathbf{A}} \mathcal{V}({\rm cat}_0)=\min_{{\rm cat}\in \mathbf{A}}{\int_{T_1}^{T}\int_0^1{r^2({\rm cat}_0; t,x)\,{\rm d}x}}{\rm d}t,\label{opt2a} \nonumber \end{align} \tag{ 4.10 }$$

subject to Problem I, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \min_{{\rm cat}\in \mathbf{A}_{{\rm cat}_0}}\tilde{\mathcal{V}}({\rm cat}_0)=\min_{{\rm cat}_0\in \mathbf{A}_{{\rm cat}_0}}{\int_{T_1}^{T}{\int_0^1{r^2({\rm cat}_0; t,x)\,{\rm d}x}}{\rm d}t} \label{opt2b} \nonumber \end{align} \tag{ 4.11 }$$

subject to Problem II.

Due to the continuous embedding $\newcommand{\bfX}{\mathbf{X}} L^2(Q_T) \subset \bfX_T$, the projection operator:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bfX}{\mathbf{X}} \displaystyle \mathbf{P}_{T_1,T}: \bfX_T \longrightarrow L^2([T_1,T]\times\Omega), \quad r({\rm cat}_0;t,x) \longmapsto r({\rm cat}_0;t,x)|_{[T_1,T]} \nonumber \end{align} \tag{ 4.12 }$$

is linear and bounded. Denoting by $<\cdot,\cdot>$ the inner product in the Hilbert space $L^2(Q_T)$ and $<\cdot,\cdot>_{L^2([T_1,T]\times \Omega)}$ the inner product in the Hilbert space $L^2([T_1,T]\times \Omega)$, we see that the functional $\mathcal{V}$ in (4.10) is Fréchet differentiable with respect to cat$_0\in W^{1,\infty}(\Omega)$. Indeed, taking $\delta\in W^{1,\infty}(\Omega)$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:4:adjoint:2} &\mathcal{V}({\rm cat}_0+\delta)-\mathcal{V}({\rm cat}_0)=\Vert\mathbf{P}_{T_1,T} r({\rm cat}_0+\delta) \Vert ^2_{L^2([T_1,T]\times \Omega)} - \Vert\mathbf{P}_{T_1,T} r({\rm cat}_0) \Vert ^2_{L^2([T_1,T]\times \Omega)}\nonumber \\ & =\Vert\mathbf{P}_{T_1,T} r({\rm cat}_0+\delta)-\mathbf{P}_{T_1,T} r({\rm cat}_0) +\mathbf{P}_{T_1,T} r({\rm cat}_0) \Vert ^2_{L^2([T_1,T]\times \Omega)} - \Vert\mathbf{P}_{T_1,T} r({\rm cat}_0) \Vert ^2_{L^2([T_1,T]\times \Omega)}\nonumber \\ & =2< \mathbf{P}_{T_1,T} (r({\rm cat}_0+\delta)-r({\rm cat}_0),\mathbf{P}_{T_1,T} r({\rm cat}_0)>_{L^2([T_1,T]\times \Omega)}\nonumber \\ & \quad \quad +\Vert\mathbf{P}_{T_1,T} (r({\rm cat}_0+\delta)-r({\rm cat}_0)) \Vert ^2_{L^2([T_1,T]\times \Omega)}. \nonumber \end{align} \tag{ 4.13 }$$

As we will see, the derivative of ${\mathcal V}$ with respect to cat₀ will involve the adjoint of the derivative $r_{{\rm cat}_0}(\cdot)$. In order to compute this adjoint, we consider $r_{{\rm cat}_0}(\cdot)$ and $c_{{\rm cat}_0}(\cdot)$ as mappings between L₂-spaces. This is a common approach, this technique has been used in several papers, e.g. [8]. In doing so, we suppose that there exist the linear and bounded operators:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{r_cat0} r_{{\rm cat}_0}(\cdot): L^2(\Omega) \longrightarrow L^2(Q_T), \quad \delta \longmapsto r_{{\rm cat}_0}(\cdot)\delta, \nonumber \end{align} \tag{ 4.14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{c_cat0} c_{{\rm cat}_0}(\cdot): L^2(\Omega) \longrightarrow L^2(Q_T), \quad \delta \longmapsto c_{{\rm cat}_0}(\cdot)\delta. \nonumber \end{align} \tag{ 4.15 }$$

Since these operators are linear and bounded, their adjoint operators $(r_{{\rm cat}_0}(\cdot)){}^*: L^2(Q_T)\longrightarrow L^2(\Omega)$, and $(c_{{\rm cat}_0}(\cdot)){}^*: L^2(Q_T)\longrightarrow L^2(\Omega)$ are also so. From this, we immediately get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{V}({\rm cat}_0+\delta)-\mathcal{V}({\rm cat}_0)&=2< \mathbf{P}_{T_1,T} r_{{\rm cat}_0}\delta, \mathbf{P}_{T_1,T}r({\rm cat}_0)>_{L^2([T_1,T]\times \Omega)}+o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)})\nonumber \\ &= 2< (r_{{\rm cat}_0}(\cdot))^* (\mathbf{P}_{T_1,T})^* \mathbf{P}_{T_1,T}, \delta >_{L^2(\Omega)} +o(\Vert \delta\Vert_{W^{1,\infty}(\Omega)}). \nonumber \end{align} \tag{ 4.16 }$$

Hence, $\mathcal{V}$ is Fréchet differentiable with respect to cat$_0\in {W^{1,\infty}(\Omega)}$.

We also see that the problem of minimizing $\mathcal{V}({\rm cat}_0)$ with respect to cat$\newcommand{\bA}{\boldsymbol{A}} _0(x) \in \bA$ defined by (4.9) admits a solution.

Denote the function z by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:04:z} z({\rm cat}_0;t,x)= \left\{\begin{array}{@{}l@{}} 0,\, {{\rm if}} \,t\in [0,T_1), \nonumber \\ r({\rm cat}_0;t,x),\, {{\rm if}}\, t\in [T_1,T]. \end{array}\right. \nonumber \end{align} \tag{ 4.17 }$$

By the Riesz representation theorem we see that the Fréchet derivative $\mathcal{V}_{{\rm cat}_0}(\cdot)$ can be identified with $2 (r_{{\rm cat}_0}(\cdot)){}^{*} z({\rm cat}_0)$. Hence

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (\mathbf{P}_{T_1,T})^*\mathbf{P}_{T_1,T}r({\rm cat}_0)=z({\rm cat}_0;t,x)= \left\{\begin{array}{@{}l@{}} 0,\, {{\rm if}} \,t\in [0,T_1), \nonumber \\ r({\rm cat}_0;t,x),\, {{\rm if}}\, t\in [T_1,T]. \end{array}\right. \nonumber \end{align*}$$

Now we present the scheme for calculating the adjoint operator $(r_{{\rm cat}_0}(\cdot)){}^*$ for Problem I. With the assumptions of smoothness on $r,c$ and positivity of r we can transform Problem I: (2.1a)–(2.1e) into the problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_t=c_{xx}+\frac{2r_x c_x}{r}-\frac{(2r_t+\frac{D_1}{\int_0^1{r c}{\rm d}x})c}{r}, \, \, (t,x)\in Q_T, \label{subeq:I_cat_after:1} \nonumber \end{align} \tag{ 4.18a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_t=-\frac{D_2}{\int_0^1{r c}{\rm d}x}{\rm cat}_0(x)c,\, \, (t,x) \in Q_T,\label{subeq:I_cat_after:2} \nonumber \end{align} \tag{ 4.18b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(0,x)=r_0(x)=1, \, \, x\in \Omega,\label{subeq:I_cat_after:3} \nonumber \end{align} \tag{ 4.18c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c(0,x)=h(x), \quad x \in \Omega,\label{subeq:I_cat_after:4} \nonumber \end{align} \tag{ 4.18d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c(t,0)=1, \, \, t\in [0,T],\label{subeq:I_cat_after:5} \nonumber \end{align} \tag{ 4.18e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_x(t,1)=0, t\in [0,T], \label{subeq:I_cat_after:6} \nonumber \end{align} \tag{ 4.18f  }$$

where $\gamma$ and $\beta$ are defined in (3.1). Now, as in the previous section, supposing that there exists a solution of the problem for the original model in the space $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T \times \bfY_T$, and $\newcommand{\e}{{\rm e}} \eta$ and $\varphi$ are defined by $\newcommand{\e}{{\rm e}} {\eta}:=r_{{\rm cat}_0}(\cdot)\delta,\, {\varphi}:=c_{{\rm cat}_0}(\cdot)\delta$, for $\delta\in W^{1,\infty}(\Omega)$, we introduce the sensitivity system:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.26} \eta_t+\frac{D_2\,{\rm cat}_0 \varphi}{\int_0^1 r\,c {\rm d}x}-\frac{D_2 {\rm cat}_0 c}{(\int_0^1 r c {\rm d}x)^2}\int_0^1 r\,\varphi+c\,\eta {\rm d}x = -\frac{D_2\,\delta\,c}{\int_0^1 r\,c {\rm d}x}, \quad (t,x) \in Q_T, \nonumber \end{align} \tag{ 4.19 }$$

with the initial condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.27} \eta(0,x)=0, \quad x \in \Omega, \nonumber \end{align} \tag{ 4.20 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.28} &\varphi_t-\varphi_{xx}-\frac{2\,r_x\varphi_x}{r}+\left(\frac{2r_t}{r}+\frac{D_1}{r\int_0^1 r \,c {\rm d}x}\right)\varphi=\frac{2c_x\,\eta_x}{r} -\frac{2c\,\eta_t}{r} \nonumber \\ & \quad \quad+\left(\frac{2r_t\,c}{r^2}-\frac{2r_x\,c_x}{r^2}+ \frac{D_1\,c}{r^2\,\int_0^1 r c {\rm d}x}\right)\eta+\frac{D_1\,c}{r(\int_0^1 r H_\epsilon(c) {\rm d}x)^2}\int_0^1 c \,\eta+r\,c\varphi {\rm d}x, \quad (t,x) \in Q_T \nonumber \end{align} \tag{ 4.21 }$$

with initial and boundary conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.29} \left\{\begin{array}{@{}l@{}} \varphi(0,x)=0, \quad x \in \Omega, \nonumber \\ \varphi(t,0)=0, \quad t \in [0,T], \nonumber \\ \varphi_x(t,1)=0, \quad t \in [0,T]. \end{array}\right. \nonumber \end{align} \tag{ 4.22 }$$

We will apply the adjoint technique by introducing an adjoint system corresponding to the above sensitivity system in order to calculate $(r_{{\rm cat}_0}(\cdot)){}^{*} z({\rm cat})$.

The structure of the sensitivity system is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4:rem:adjoint} \left\{\begin{array}{@{}l@{}} D_1\varphi+L_1\eta=l\delta, \nonumber \\ D_2\varphi+L_2\eta=0. \end{array}\right. \nonumber \end{align} \tag{ 4.23 }$$

Using by integration by parts and summing up the two equations, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle \varphi,\widetilde{D}_1\,v+\widetilde{D}_2\,w \rangle +\langle \eta,\widetilde{L}_1\,v+\widetilde{L}_2\,w \rangle=\langle \delta,lv \rangle. \nonumber \end{align} \tag{ 4.24 }$$

This then yields the following theorem.

Adjoint system for Problem I. Suppose that the system

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & -w_t-w_{xx}+\frac{2r_x}{r}\,w_x+\left(\frac{2r_{xx}r-2r^2_x+ 2r_t\,r}{r^2}+\frac{D_1}{r\int_0^1\,r\,c\,{\rm d}x}\right)w \nonumber \\ & + \frac{D_2\,{\rm cat}(x)}{\int_0^1 r\,c {\rm d}x}\,\nu-r \int_0^1\,\frac{D_2 {\rm cat}(x) c}{(\int_0^1 r c {\rm d}x)^2}\,\nu\, {\rm d}x- -r \int_0^1\,\frac{D_1 c w}{r(\int_0^1 r c {\rm d}x)^2}\, {\rm d}x=0, \quad (t,x) \in Q_T, \label{4.41} \nonumber \end{align} \tag{ 4.25 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & - \frac{2c}{r}w_t+\frac{2\,c_x}{r}\,w_x-\left(\frac{D_1\,c}{r^2\int_0^1 r c {\rm d}x}-\frac{2c_{xx}}{r}+ \frac{2c_t\,r}{r^2}\right)w-c\,\int_0^1\,\frac{D_1\,c\,w}{r(\int_0^1 r c {\rm d}x)^2}{\rm d}x \nonumber \\ & -\nu_t-c\int_0^1 \frac{D_2 {\rm cat}(x) c}{(\int_0^1 r c {\rm d}x)^2}\,\nu\,{\rm d}x=z, \quad (t,x) \in Q_T,\label{4.42} \nonumber \end{align} \tag{ 4.26 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu(T,x)=0, \quad w(T,x)=0, \quad x \in \Omega,\label{4.35*} \nonumber \end{align} \tag{ 4.27 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w(t,0)=0, w_x(t,1)=0, \quad t \in [0,T],\label{4.37} \nonumber \end{align} \tag{ 4.28 }$$

has a unique solution $(\nu,w)\in W^{1,\infty}(Q_T)\times W(0,T)$ for some positive T. Suppose further that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.30} {\rm cat}_x(1)=0. \nonumber \end{align} \tag{ 4.29 }$$

Then, the operator $(r_{{\rm cat}_0}(\cdot)){}^{*}: L^2(Q_T) \mapsto L^2(\Omega)$, $z \mapsto (r_{{\rm cat}_0}(\cdot)){}^{*} z$ has the form

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (r_{{\rm cat}_0}(\cdot))^{*} z({\rm cat})= \int_0^T{\frac{D_2\,c}{\int_0^1 r\,c {\rm d}x} \nu}{\rm d}t.\nonumber \end{align*}$$

Now present the scheme of calculating the adjoint operator $(r_{{\rm cat}_0}(\cdot)){}^*$ for Problem II.

Problem II leads to system

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_t=c_{xx}+\frac{2r_x c_x}{r}-\frac{(2r_t+\frac{D_1}{\int_0^1{r c}{\rm d}x})c}{r}, \, \, (t,x)\in Q_T, \label{subeq:I_cat_t_after:1} \nonumber \end{align} \tag{ 4.30a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_t=-\frac{D_2}{\int_0^1{r c}{\rm d}x}{\rm cat}(t,x)c,\, \, (t,x) \in Q_T,\label{subeq:I_cat_t_after:2} \nonumber \end{align} \tag{ 4.30b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(0,x)=r_0(x)=1, \, \, x\in \Omega,\label{subeq:I_cat_t_after:3} \nonumber \end{align} \tag{ 4.30c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c(0,x)=h(x), \quad x\in \Omega, \label{subeq:I_cat_t_after:4} \nonumber \end{align} \tag{ 4.30d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c(t,0)=1, \, \, t\in [0,T],\label{subeq:I_cat_t_after:5} \nonumber \end{align} \tag{ 4.30e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_x(t,1)=0,\label{subeq:I_cat_t_after:6} \nonumber \end{align} \tag{ 4.30f  }$$

where ${\rm cat}(t,x)$ is defined by (2.8) and the initial condition (2.9) for given $\alpha_0$ and R₀.

Assume that $\newcommand{\bfX}{\mathbf{X}} r\in\bfX_T$ and $\newcommand{\bfY}{\mathbf{Y}} c\in \bfY_T$ are Fréchet differentiable with respect to cat$_0\in W^{1,\infty}(\Omega)$ for T small enough. Denote $\newcommand{\e}{{\rm e}} \bar{\eta}:=r_{{\rm cat}_0}(\cdot)\Delta,\, \bar{\varphi}:=c_{{\rm cat}_0}(\cdot)\Delta,\, {\psi}:={\rm cat}_{{\rm cat}_0}(\cdot)\Delta$, for $\delta\in W^{1,\infty}(\Omega)$. Similarly to the previous one, suppose that there exists a solution of Problem II in the space $\newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \bfX_T \times \bfY_T$. Then, we get the sensitivity system:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:04.26*} \bar{\eta}_t=-\frac{D_2}{\int_0^1 r\,c {\rm d}x}\left(c\psi+{\rm cat}(t,x)\bar{\varphi}- \frac{{\rm cat}(t,x) c}{\int_0^1 r c {\rm d}x}\int_0^1 r\big(\bar{\varphi}+c\,\bar{\eta}\big) {\rm d}x\right),\quad (t,x) \in Q_T, \nonumber \end{align} \tag{ 4.31 }$$

with the initial condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:04.27*} \bar{\eta}(0,x)=0, \quad x \in \Omega \nonumber \end{align} \tag{ 4.32 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:04.28*} &\bar{\varphi}_t-\bar{\varphi}_{xx}-\frac{2\,r_x\bar{\varphi}_x}{r}+\left(\frac{2r_t}{r}+\frac{D_1}{r\int_0^1 r \,c {\rm d}x}\right)\bar{\varphi}=\frac{2c_x\,\bar{\eta}_x}{r} -\frac{2c\,\bar{\eta}_t}{r} \nonumber \\ &+\left(\frac{2r_t\,c}{r^2}-\frac{2r_x\,c_x}{r^2}+ \frac{D_1\,c}{r^2\,\int_0^1 r c {\rm d}x}\right)\bar{\eta}+\frac{D_1\,c}{r(\int_0^1 r H_\epsilon(c) {\rm d}x)^2}\int_0^1 c \,\bar{\eta}+r\,c\bar{\varphi} {\rm d}x, \quad (t,x) \in Q_T \nonumber \end{align} \tag{ 4.33 }$$

with initial and boundary conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{\varphi}(0,x)=0, \quad \bar{\varphi}(t,0)=0, \quad \bar{\varphi}_x(t,1)=0, \nonumber \end{align} \tag{ 4.34 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi_t=\alpha_0\int_0^1{G(x-y)T_{R_0}'(r({\rm cat}_0;y,t))\bar{\eta}(y,t){\rm d}y}, \nonumber \end{align} \tag{ 4.35 }$$

with initial condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi(0,x)=\Delta(x). \nonumber \end{align} \tag{ 4.36 }$$

Denote the function $\mathbf{z}$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:04:z*} \mathbf{z}({\rm cat}_0;t,x)= \left\{\begin{array}{@{}l@{}} 0,\, {{\rm if}}\, t\in [0,T_1) \nonumber \\ \bar{\mathbf{P}}_{T_1,T} r({\rm cat}_0;t,x),\, {{\rm if}} \, t\in [T_1,T] \end{array}\right. = \left\{\begin{array}{@{}l@{}} 0,\, {{\rm if}}\, t\in [0,T_1) \nonumber \\ r({\rm cat}_0;t,x),\, {{\rm if}}\, t\in [T_1,T]. \end{array}\right. \nonumber \end{align} \tag{ 4.37 }$$

Under assumptions about operators $(r_{{\rm cat}_0}(\cdot)){}^*$ and $(c_{{\rm cat}_0}(\cdot)){}^*$ in (4.14) and (4.15), the functional $\tilde{\mathcal{V}}$ defined in (4.11) is Fréchet differentiable with respect to cat$_0\in L^2(\Omega)$. Moreover, according to Riesz representation theorem the Fréchet derivative $\tilde{\mathcal{V}}_{{\rm cat}_0}(\cdot)$ is identified by

$2(r_{{\rm cat}_0}(\cdot)){}^*(\mathbf{P}_{T_1,T}){}^* \mathbf{P}_{T_1,T}r({\rm cat}_0)$. From which we see that $\tilde{\mathcal{V}}_{{\rm cat}_0}(\cdot)$ can also be identified by

$2(r_{{\rm cat}_0}(\cdot)){}^{*}\mathbf{z}({\rm cat})$, where $\mathbf{z}({\rm cat}_0)$ is defined in (4.37). Hence

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (\mathbf{P}_{T_1,T})^*\mathbf{P}_{T_1,T}r({\rm cat}_0)=\mathbf{z}({\rm cat}_0;t,x)=\left\{\begin{array}{@{}l@{}} 0,\, {{\rm if}} \,t\in [0,T_1), \nonumber \\ r({\rm cat}_0;t,x),\, {{\rm if}}\, t\in [T_1,T]. \end{array}\right.\nonumber \end{align*}$$

By integration by parts, we can calculate $(r_{{\rm cat}_0}(\cdot)){}^{*} \mathbf{z}({\rm cat}_0)$ by the adjoint system as follows:

Adjoint system for problem II: Suppose that the system

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & -w_t-w_{xx}+\frac{2r_x}{r}\,w_x +\left(\frac{2r_{xx}r-2r^2_x + 2r_t\,r}{r^2}+\frac{D_1}{r\int_0^1\,r\,c\,{\rm d}x}\right)w+ \frac{D_2\,{\rm cat}(t,x)}{\int_0^1 r\,c {\rm d}x}\,\nu \nonumber \\ & \quad \quad -r \int_0^1\,\frac{D_2 {\rm cat}(t,x) c}{(\int_0^1 r c {\rm d}x)^2}\,\nu\, {\rm d}x-r \int_0^1\,\frac{D_1 c w}{r(\int_0^1 r c {\rm d}x)^2}\, {\rm d}x=0, \quad (t,x) \in Q_T, \nonumber \end{align} \tag{ 4.38 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \quad - \frac{2c}{r}w_t+\frac{2\,c_x}{r}\,w_x-\left(\frac{D_1\,c}{r^2\int_0^1 r c {\rm d}x}-\frac{2c_{xx}}{r}+ \frac{2c_t\,r}{r^2}\right)w-c\,\int_0^1\,\frac{D_1\,c\,w}{r(\int_0^1 r c {\rm d}x)^2}{\rm d}x \nonumber \\ &-\nu_t-c\int_0^1 \frac{D_2 {\rm cat}(t,x) c}{(\int_0^1 r c {\rm d}x)^2}\,\nu\,{\rm d}x + \alpha_0 {T'}_{R_0}({\rm cat}_0,x,t)\int_0^1{\zeta(y,t)G(y-x){\rm d}y}=\mathbf{z}, \quad (t,x) \in Q_T, \label{04.42} \nonumber \end{align} \tag{ 4.39 }$$

with the conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{04.35*} \nu(T,x)=0, \quad w(T,x)=0, \quad x \in \Omega, \nonumber \end{align} \tag{ 4.40 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w(t,0)=0, \quad w_x(t,1)=0, \quad t \in [0,T], \nonumber \end{align} \tag{ 4.41 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \zeta_t=-\frac{D_2\,c v}{\int_0^1{r c} {\rm d}x} \nonumber \end{align} \tag{ 4.42 }$$

with the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \zeta(T,x)=0, \quad x \in \Omega, \nonumber \end{align} \tag{ 4.43 }$$

has a unique solution $(\nu,w,{\rm cat})\in W^{1,\infty}(Q_T)\times W(0,T)\times W^{1,\infty}(Q_T)$. Additionally, suppose that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm cat}_{0x}(1)=0. \nonumber \end{align} \tag{ 4.44 }$$

Then, the operator $(r_{{\rm cat}_0}(\cdot)){}^*: L^2(Q_T) \mapsto L^2(\Omega)$; $z \mapsto -\zeta(0,\cdot)$ has the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (r_{{\rm cat}_0}(\cdot))^*z({\rm cat}) = -\zeta(0,x). \nonumber \end{align} \tag{ 4.45 }$$

We note that we have applied the adjoint and sensibility systems to gradient methods but these gave good numerical results only for short time processes [17], therefore in the next section we use directional gradients for iterative methods.

As a starting point we implement the optimization scheme of [4] for comparison with our more involved schemes in the following subsections. That is, we consider the case that the catalytic function cat₀ takes values 0 or 1 at positions along the surface of the pore. In discretized form we have cat$_j^0, j=0,\cdots,(J-1)$, where cat$_j^0={\rm cat}_0(x_j)$. We emphasize, that the update for cat_t is done a simpler way in this approach: whenever $r^n_j$ reaches a value smaller than R₀, then both neighboring catalyst positions are set to 1, i.e. ${\rm cat}^{n+1}_{j-1}={\rm cat}^{n+1}_{j+1}=1$.

Before starting the algorithm all values are set to zero. In the first round of the greedy algorithm all potential catalyst positions are tested separately, i.e. for $\newcommand{\e}{{\rm e}} \ell = 0,.., (J-1)$ we set cat$\newcommand{\e}{{\rm e}} _j^0=\partial_{j,\ell}$ and run the algorithm until clogging. We determine the index $\newcommand{\e}{{\rm e}} \ell$, which returns the minimal free volume. At this position cat$\newcommand{\e}{{\rm e}} _\ell$ is fixed as one and the greedy algorithms proceeds with the remaining indices. Hence, after every full round of the greedy algorithm the number of non-zero cat$_j^0$ values increases by one. The optimal value for the remaining free volume was obtained with 52 non-zero catalyst positions, see figure 3 right. For the experiments presented in this section, we actually stopped the process prior to clogging, more precisely we ran the algorithm until $\min_j r_j^n \leqslant r_0$ with r₀  =  0.1.

Figures 2 and 3 show the evolution of pore radii and concentration of oxygen along the pore in 3D plotting and in the line plotting over time with the optimal catalyst number of 52, which is illustrated in the figure 4. The optimal volume we deduce is $22.8\%$. In this experiment, we take $J=100, N=4000$, threshold of R₀  =  0.96. Figure 4 demonstrates the dependency of free volume $V_{{{\rm free}}}$ before pore clogging on the number of catalyst positioned along the pore.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows on the left that the greedy algorithm places most catalysts at the end of the pore. This is to be expected, since placing catalysts at the beginning of the pore, where the oxygen concentration is still high will lead to a rapid clogging of the pore at a position near the beginning of the pore and a large free volume in the later part. This is confirmed by figure 3 on the right, which indicates, that as time evolves more and more oxygen is consumed at the beginning of the pore. Hence rather little oxygen is transported towards the later parts of the pore. Hence, no catalysts are placed at the beginning of the pore before position 0.4 but after the growing process clogging occurs at position 0.34, see figure 3 on the left. Figure 4 on the right shows the evolution of the free volume as more and more catalysts are added by the greedy algorithm. Some interesting local minima occur, i.e. stopping the classical greedy algorithm in the local minimum after 24 iterations yields a suboptimal result.

Also, figure 3 on the left shows some drawbacks of the classical greedy algorithm, which put the initial catalysts activity either 0 or 1. This does not allow for a smooth development in the catalyst reaction process and results in the spike just before position 0.6

A single run of the forward pass requires approx. 12 s. The running time for testing all profiles is 142h with Intel Xeon Processor E5-2687W v4 with 12 cores and 1.54 TB of memory. We apply OpenMP application program interface to carry out parallel computations for optimizing the running time.

The results of the previous subsection are obtained with the version of the greedy algorithm for catalyst optimisation which was used in [4], hence the optimized value of approx. $20 \%$ free volume at the time of clogging is the reference value for the more involved algorithms of this and the next subsection.

In this section we simply improve the greedy algorithm by choosing smaller increments, i.e. 0.1 instead of 1, and by allowing to update to the same position of catalyst for several times. Hence, all J positions have to be tested in every iteration of this incremental greedy algorithm. The resulting cat₀-function has a more geometric structure which also reflects engineering intuition about optimal catalyst positioning.

There is no limit to the maximal value of cat$_j^0$, hence, we cannot test all possible cat₀-profiles. Instead we use the first local minimum of $V_{\rm free}$ as stopping criterion. This minimum occurred after the 596th iteration, the optimal $V_{{{\rm free}}}$ is $7.8\%$. The parameters we used for this experiment are: J  =  100 and $\Delta t=\frac1{100-1}$ and $\alpha_0=0.01$. Figure 5 shows radii and concentration inside the pore before clogging corresponding to the optimal catalytic function cat₀ in 3D plotting by incremental Greedy in Problem II, while the in line plots of radii and concentration inside the pore are shown in figure 6. The optimal catalytic function cat₀ is plotted in figure 7 (left) and the dependence of free volume before pore clogging $V_{{{\rm free}}}$ with respect to numbers of increments is demonstrated in figure 7 (right).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared with the results obtained by the classical greedy algorithm in the previous subsection figure 5 shows a smoother transition of the catalyst growth and oxygen concentration along the pore. Most interesting is the final catalyst distribution after the incremental greedy stopped after 600 iterations, see figure 7. It matches the qualitative expectations, that catalysts should be placed predominantly towards the end of the pore in order to minimize the free volume after clogging. However, these simulation add a quantitative estimate to this behavior and yields an almost optimal distribution in terms of the free volume of $\% 7.8$ after clogging, see figure 7 on the right. This can only minimally improved by further analytic iterations see next subsection. Again, the convergence towards a minimum is rather discontinuous as the number of incremental greedy iterations increases, see figure 7 on the right. Figures 5 and 6 both display the same data, namely the development of the pore radius and the of the oxygen concentration over time. Figure 5 allows a more general visualization of this process. However, figure 6 shows a rather dynamic local behavior, which might be overcome by an even finer incremental step size. However, using OpenMP application program interface for implementing parallel computations on an Intel Xeon Processor E5  −  2687W $v4$ still results in a running time of 163h. Hence, a refined discretisations with smaller increments is beyond any feasible computation time for this comparative study.

For Problem II, we consider the first optimization Problem: find cat₀ such that the free volume $V({\rm cat}_0)$ of the pore after some time T defined by (2.11) reaches to a minimum.

The optimization problem is solved by the gradient descent method

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm cat}_0^{\rm new}={\rm cat}_0^{\rm old}-\tau V_{{\rm cat}_0^{\rm old}} . \nonumber \end{align} \tag{ 5.3 }$$

The iteration stops when the free volume $V({\rm cat}_0^{\rm new})$ is greater than the previous one $V({\rm cat}_0^{\rm old})$. In the context of inverse problems such parameter identification problems for partial differential equations have been studied intensively, see e.g. [2] for an application in scattering or [10] for a recent survey.

The update of cat₀ uses the derivative of $V$ with respect to cat₀. According to section 4,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle V_{{\rm cat}_0}= r^*_{{\rm cat}_0} \mathbf{P}_T^* \mathbf{P}_T r({\rm cat}_0).\nonumber \end{align*}$$

This can be computed by solving the adjoint system for $r^*_{{\rm cat}_0}$ with $z= \mathbf{P}_T^* \mathbf{P}_T r({\rm cat}_0)$, see (4.8). However, this is a parabolic equation backward in time and requires to store the solution of the forward problem at every time step. This is beyond the capacity of the computer used in these experiment.

Hence, suppose that $\mathbf{P}^{(k)}_T r^{(k)}_{{\rm cat}_0}$ is approximated for $\mathbf{P}_T r_{{\rm cat}_0}$ at the iteration k of the optimization algorithm, then $\mathbf{P}^{(k)}_T r^{(k)}_{{\rm cat}_0}$ is a $J\times J$ matrix, therefore the adjoint operator is defined by its transpose matrix $(\mathbf{P}^{(k)}_T r^{(k)}_{{\rm cat}_0}){}^*=(\mathbf{P}^{(k)}_T r^{(k)}_{{\rm cat}_0}){}^{\rm T}$. In order to determine the derivative $r^{(k)}_{{\rm cat}_0}$ of pore radius r with respect to cat₀, we calculate the directional derivatives by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (r_{{\rm cat}_{0}})_j=\frac{r({\rm cat}_0+e_j h)-r({\rm cat}_0)}{h}, \nonumber \end{align} \tag{ 5.4 }$$

where h is small enough positive number and $e_j=(0,\cdot, 1, \cdots, 0)$ is a $(1\times J)$ vector, taking value of 1 at the j th coordinate.

For the experiment of the gradient descent, we take the initial catalytic function from the 596th step of the incremental greedy algorithm. The resulting catalytic function is plotted in figure 8 (left). The parameters for gradient descent method are

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle h=10^{-5},\, \alpha_0=0.01, \, N=J=100, \, R_0=0.96, T=1571\,393\times \Delta t,\nonumber \end{align*}$$

and the regularization parameter is chosen from a decreasing sequence

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \alpha^{k}=10^{-7}, \quad \forall k \geqslant 0.\nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

The numerical iteration for the incremental greedy was stopped as soon as $r(t,x) < 0.1$ at some position x for the first time. In the presented experiment this point was reached after $T=1571\,393$ time steps.

In the following experiments using the gradient descent method, we kept $T=1571\,393$ for simplicity. After 48 iterations of the gradient descent, we get the local minima which refers us to a free volume $V{{{\rm free}}}=7,74\%$. Hence, we get a minor improvement for the free volume $V_{{{\rm free}}}$: $0.115\%$. Figure 8 (right) shows the changing amount(increasing over time) of catalytic function by gradient descent method after 48 iterations. The running time is 14 h with Intel processor core i5-5200U and memory size is 12 GB.