\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2009(2009), No. 83, pp. 1--22.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2009 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2009/83\hfil Optimal control for a thermistor problem] {Optimal boundary control for a time dependent thermistor problem} \author[V. Hrynkiv \hfil EJDE-2009/83\hfilneg] {Volodymyr Hrynkiv} \address{Volodymyr Hrynkiv \newline Department of Computer and Mathematical Sciences, University of Houston-Down\-town, One Main St, Houston, TX 77002, USA} \email{HrynkivV@uhd.edu} \thanks{Submitted December 10, 2008. Published July 10, 2009.} \subjclass[2000]{49J20, 49K20} \keywords{Optimal control; thermistor problem; elliptic systems} \begin{abstract} An optimal control of a two dimensional time dependent thermistor problem is considered. The problem consists of two nonlinear partial differential equations coupled with appropriate boundary conditions which model the coupling of the thermistor to its surroundings. Based on physical considerations, an objective functional to be minimized is introduced and the convective boundary coefficient is taken to be the control. Existence of solutions to the state system and existence of the optimal control are proven. To characterize this optimal control, the optimality system consisting of the state and adjoint equations is derived. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction} Thermistor is a thermally sensitive resistor whose electrical conductivity changes drastically by orders of magnitude as the temperature reaches a certain threshold. Thermistors are used as temperature control elements in a wide variety of military and industrial equipment ranging from space vehicles to air conditioning controllers. They are also used in the medical field for localized and general body temperature measurement, in meteorology for weather forecasting as well as in chemical industries as process temperature sensors \cite{Kwo,Mac}. We consider the two-dimensional time-dependent thermistor problem $$\label{statet} % \label{state} \begin{gathered} u_t-\Delta u-{\sigma}(u)|\nabla\varphi|^{2} ={0 \quad \text{in } \Omega\times(0,T),} \\ {\nabla\cdot(\sigma(u)\nabla\varphi) }= 0 \quad \text{in } \Omega\times(0,T), \\ {\frac{\partial u}{\partial n}+\beta u}= 0 \quad \text{on } \Gamma_{N}\times(0,T), \\ {u}=0 \quad \text{on } \Gamma_{D}\times(0,T),\\ {\varphi}={\varphi_0 \quad\text{on } \partial\Omega\times(0,T),}\\ {u} = u_0 \quad\text{on } \Omega\times\{0\}, \end{gathered}$$ where $\varphi( x,t)$ is the electric potential, $u(x,t)$ is the temperature, and $\sigma(u)$ is the electrical conductivity. Here $n$ denotes the outward unit normal and $\partial/\partial n=n\cdot\nabla$ is the normal derivative on $\partial\Omega$. The elliptic-parabolic system \eqref{statet} describes the heating of a conductor produced by an electric current. It is the consequence of Ohm's law and Fourier's law $$\label{Ohm} I=-\sigma(u)\nabla\varphi,\quad Q=-k(u)\nabla u,$$ coupled with the conservation laws of energy and current \begin{gather} \label{conservation1} \rho c\frac{\partial u}{\partial t}+\nabla\cdot Q=I\cdot E,\\ \label{conservation2} \nabla\cdot I=0, \end{gather} where $E$ denotes the electric field, $\rho$ is the density of the conductor, and $c$ is its heat capacity. Since the nonlinearity $k(u)$ varies only slightly with $u$, it is of secondary importance, and therefore we can assume $k(u)=1$. As the density $\rho$ of the material and its specific heat $c$ typically are constants, we assume that $\rho c=1$ (see \cite{Ant}). A more detailed discussion of equations \eqref{statet} can be found in \cite{Ant,Cim, Fow,How,Shi}. Boundary conditions show how the thermistor is connected thermally and electrically to its surroundings. To describe boundary conditions for the temperature we decompose the boundary $\partial\Omega$ into two disjoint sets $\bar{\Gamma}_{D}\cup\bar{\Gamma}_N=\partial\Omega$. The Robin boundary condition for temperature represents the physically significant case of convection and/or radiation in situations where $u$ is not too different from the ambient temperature. It is precisely this boundary condition that will be used as a control. It has been observed that large temperature gradients may cause the thermistor to crack. Numerical experiments indicate (see \cite{Fow, Zho}) that low values of the heat transfer coefficient $\beta$ lead to small temperature variations. On the other hand, low values of the heat transfer coefficient result in high operating temperatures of a thermistor which are also undesirable. This motivated us to take the heat transfer coefficient as a control and to consider the optimal control problem of minimizing the heat transfer coefficient while keeping the operating temperature of the thermistor reasonably low. These physical considerations lead us to the following objective functional $$J(\beta)=\int_0^T\int_\Omega u\,dx\,dt+\int_{\Gamma_N\times(0,T)}\beta^2\,ds\,dt.$$ Denoting the set of admissible controls by $$U_M=\{\beta\in L^\infty(\partial\Omega\times(0,T)):\,0<\lambda\leq\beta\leq M\},$$ the optimal control problem is $$\text{Find \beta^*\in U_M such that J(\beta^*)=\min_{\beta\in U_M}J(\beta).} \label{OC}$$ Henceforth we use the standard notation for Sobolev spaces, we denote $\|\cdot\|_p=\|\cdot\|_{L^p(\Omega)}$ for each $p\in [1,\infty]$; other norms will be clearly marked. We also denote $Q:=\Omega\times (0,T)$ and $\Gamma_{N}^T:=\Gamma_N\times(0,T)$. Analysis of both the steady-state and time-dependent thermistor equations with different types of boundary and initial conditions has received a great deal of attention. See \cite{Ant,Cim,Cim1,Cim2,Fow,How,Shi,Xu1,Xu2,Yua} for existence of weak solutions, existence/nonexistence of classical solutions, uniqueness and related regularity results in different settings with various assumptions on the coefficients. We mention in passing that few authors consider a thermistor problem which consists of two parabolic equations coupled with some boundary conditions (see \cite{Lee, Rod1}). This fully parabolic system comes from (\ref{Ohm}), (\ref{conservation1}), and the time dependent version of charge conservation (\ref{conservation2}), i.e., $$\label{conservation3} \frac{\partial q}{\partial t}+\nabla\cdot I=0.$$ However, the vast majority of the existing literature on time dependent thermistor deals with the elliptic-parabolic system \eqref{statet}. This reflects a feasible thermistor scenario where the time derivative in $(\ref{conservation3})$ is dropped considering that the relaxation time for the potential is very small, hence giving the equation (\ref{conservation2}) (see \cite{Cim4,Shi}). The first optimal control paper on a thermistor problem is the paper by Lee and Shilkin \cite{Lee} where a fully parabolic system is considered and the source is taken to be the control. Optimal control of a nonlocal elliptic-parabolic case is studied in \cite{Cim3}, where the applied potential is taken to be the control. Optimal control of a steady state thermistor problem is considered in \cite{Hry}, where the heat transfer coefficient is the control. The motivation for our work is threefold. First, the results from \cite{Hry} are extended to the time dependent elliptic-parabolic thermistor system. As it was mentioned before, such an elliptic-parabolic thermistor model represents a reasonably realistic situation where we neglect the time derivative in (\ref{conservation3}) by taking into consideration the difference in time scale with the thermal processes. One of the technical difficulties in dealing with such an elliptic-parabolic system lies in the fact that there is no information on time derivative of the potential $\varphi$, and as a result one cannot use directly the standard compactness result from \cite{Sim} to obtain strong convergence of sequences of potentials in appropriate spaces. To circumvent this difficulty one needs to incorporate the structure of the system. This is different from \cite{Lee} in that a fully parabolic system and a different type of control were considered there. Secondly, we present a complete and self-contained derivation of the optimality system, whereas in \cite{Cim3} the optimality system was given only for the special case of constant $\sigma(u)$. Thirdly, the choice of the convective boundary condition as a control for this time dependent problem seems to be quite appropriate from point of view of practical applications. In section \ref{apriori} we derive a priori estimates under the assumption of small boundary data and prove existence of solutions to the state system. In section \ref{exist} we prove existence of an optimal control. The optimality system is derived and an optimal control is characterized in section \ref{optimality}. \section{A priori estimates and existence of solutions to the state system}\label{apriori} We make the following assumptions: \begin{itemize} \item $\Omega\subset{R}^{2}$ is a bounded domain with smooth boundary; \item $\sigma(s)\in C^1({R})$, $00$. Let $u$ be the weak solution of the following problem \begin{gather*} -\nabla\cdot(a\nabla u)=\nabla\cdot g \quad\text{in }\Omega\\ u=0 \quad\text{on }\partial\Omega. \end{gather*} Then for each $p>2$ there exists a positive constant $c^*$ depending only on $n, \Omega, a$, and $p$ such that if $g\in (L^2(\Omega))^n$ then we have $\| \nabla u\|_p\leq c^* (\|g \|_p+\| \nabla u\|_2).$ \end{lemma} Let $V_D:=\{v\in H^1(\Omega):v=0 \text{ on }\Gamma_D\}$. By a weak solution to \eqref{statet} we mean a pair $(u,\varphi)$ such that $u\in L^2(0,T;V_D(\Omega)), \varphi-\varphi_0\in L^2(0,T;H_0^1(\Omega))$ and $$\begin{gathered} \label{eq01} \int_0^T\langle u_t,v\rangle +\int_Q\nabla u\cdot\nabla v\,dx\,dt+\int_{\Gamma_{N}^T}\beta u v\,ds\,dt=\int_Q\sigma(u)|\nabla\varphi|^2v\,dx\,dt, \\ \int_Q\sigma(u)\nabla \varphi\cdot\nabla w\,dx\,dt=0, \end{gathered}$$ for all $v\in L^2(0,T;H^{1}(\Omega))\cap L^\infty(Q)$ and $w\in L^2(0,T;H^{1}_{0}(\Omega))$, where $\langle \cdot,\cdot\rangle$ denotes duality bracket between $V_D'$ and $V_D$, with $V_D'$ being the dual of $V_D$. The quadratic term on the right hand side of the first equation in (\ref{eq01}) can be dealt with as in \cite{Xu4, Shi}, and we obtain %for a.e. $t\in (0,T)$ \label{eq1} \begin{gathered} \begin{aligned} &\int_0^T\langle u_t,v\rangle\,dt+\int_Q\nabla u\cdot\nabla v\,dx\,dt+\int_{\Gamma_N^T}\beta u v\,ds\,dt\\ &= \int_Q (\varphi_{0}-\varphi)\sigma(u)\nabla\varphi\cdot\nabla v\,dx\,dt +\int_Q(\sigma(u)\nabla \varphi\cdot\nabla\varphi_{0})\,v\,dx\,dt, \end{aligned}\\ \int_Q\sigma(u)\nabla \varphi\cdot\nabla w\,dx\,dt= 0, \end{gathered} for all $v\in L^2(0,T;H^{1}(\Omega))\cap L^\infty(Q)$ and $w\in L^2(0,T;H^{1}_{0}(\Omega))$. By taking $p=r$ in Lemma \ref{meyers}, we can say that for each $r>2$, $\varphi(\cdot,t)\in W^{1,r}(\Omega)$ for a.e. $t\in(0,T)$. Therefore, since $\nabla\varphi(\cdot,t)\in L^{r}(\Omega)$ and $\nabla\varphi_0(\cdot,t)\in L^{2}(\Omega)$ a.e. $t\in(0,T)$, it follows that there exists $s'$ such that $\nabla\varphi(\cdot,t)\cdot\nabla\varphi_0(\cdot,t)\in L^{s'}(\Omega)$ and $$\label{conjugate1} \frac{1}{s'}=\frac{1}{2}+\frac{1}{r}.$$ Next, since $\Omega\subset{R}^2$ it follows $v(\cdot,t)\in H^1(\Omega)\subset\subset L^s(\Omega)$ for $s\in[1,\infty)$, conjugate of $s'$ in (\ref{conjugate1}): \begin{eqnarray}\label{conjugate} \frac{1}{s'}+\frac{1}{s}=1. \end{eqnarray} Also, by the weak maximum principle (see \cite[formula (2.14)]{Shi}) $$\label{supremum} \|\varphi\|_{L^\infty(Q)}\leq \sup_{\Gamma_D}|\varphi_0|.$$ Note that $r>2$ is chosen from Lemma \ref{meyers}, then $s'$ and $s$ are determined and satisfy $s'<20$. Collecting the estimates and using the pde \eqref{eps2} to get the estimate of the time derivative quotient, gives \begin{gather*} \big\|\frac{u^\varepsilon-u}{\varepsilon}\big\|_{L^2(0,T;H_0^1(\Omega))} \leq c,\\ \big\|\frac{u^\varepsilon-u}{\varepsilon}\big\|_{L^2(0,T;L^2(\Gamma_N))} \leq c,\\ \big\|\frac{u_t^\varepsilon-u_t}{\varepsilon}\big\|_{L^2(0,T;H^{-1}(\Omega))}\leq c,\\ \big\|\frac{\varphi^\varepsilon-\varphi}{\varepsilon}\big\|_{L^2(0,T;H^1(\Omega))}\leq c, \end{gather*} where the generic constant $c>0$ does not depend on $\varepsilon$. Claim: we can derive a better estimate on $(\varphi^\varepsilon-\varphi)/\varepsilon$ than we have. This better estimate is needed to pass to the limit as $\varepsilon\to 0$ in the derivation of the sensitivity equations. Indeed, we can write for a.e. $t\in (0,T)$ \begin{gather*} \nabla\cdot(\sigma(u^\varepsilon)\nabla\varphi^\varepsilon)=0 \quad \text{in } \Omega,\\ \nabla\cdot(\sigma(u)\nabla\varphi)=0 \quad \text{in } \Omega. \end{gather*} Subtract these equations as well as appropriate boundary conditions to obtain $$\label{star} \begin{gathered} \nabla\cdot\Big(\sigma(u)\nabla\big(\frac{\varphi^\varepsilon -\varphi}{\varepsilon}\big)\Big) = \nabla\cdot\Big[\frac{\sigma(u^\varepsilon)-\sigma(u))} {\varepsilon}\nabla\varphi^\varepsilon\Big] \quad \text{in } \Omega, \\ \frac{\varphi^\varepsilon-\varphi}{\varepsilon}=0 \quad \text{on } \partial\Omega. \end{gathered}$$ By Lemma \ref{meyers}, for a.e. $t\in (0,T)$, we have $\|\nabla\varphi^\varepsilon \|_{2r}\leq c$, where $r>2$ and $\varphi^\varepsilon$ solves \begin{gather*} \nabla\cdot(\sigma(u^\varepsilon)\nabla\varphi^\varepsilon) = 0 \quad \text{in } \Omega,\\ \varphi^\varepsilon = \varphi_0 \quad \text{on } \partial\Omega. \end{gather*} Also, since $\Omega\subset R^2$ it follows that $H^1(\Omega)\subset L^q(\Omega)$ for any $q\in [1,\infty)$. Thus, we can write for a.e. $t\in (0,T)$, \begin{eqnarray} \big\|\frac{u^\varepsilon-u}{\varepsilon} \big\|_{2r}\leq \big\|\frac{u^\varepsilon-u}{\varepsilon} \big\|_{H^1(\Omega)}\leq c. \end{eqnarray} Therefore, returning to (\ref{star}), we obtain $$\label{m1} \big\|\nabla\big(\frac{\varphi^\varepsilon-\varphi}{\varepsilon}\big) \big\|_{r}\leq c^*\Big(\big\|\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\nabla\varphi^\epsilon \big\|_r+\big\| \nabla\big(\frac{\varphi^\varepsilon-\varphi }{\varepsilon}\big)\big\|_2 \Big).$$ We can estimate the first term on the right hand side of (\ref{m1}), \begin{align*} \big\|\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\nabla\varphi^\epsilon \big\|_r^r &= \int_{\Omega}\big| \frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\big|^r\cdot |\nabla\varphi^\varepsilon|^r\,dx\\ &\leq C\int_{\Omega}\big|\frac{u^\varepsilon-u}{\varepsilon}\big|^r\cdot|\varphi^\varepsilon|^r\,dx\\ &\leq C\Big( \big\|\frac{u^\varepsilon-u}{\varepsilon}\big\|_{2r}^{2r} +\|\nabla\varphi^\varepsilon \|_{2r}^{2r}\Big)\leq c. \end{align*} Therefore, $$\label{r} \big\|\nabla\big(\frac{\varphi^\varepsilon-\varphi}{\varepsilon}\big) \big\|_{r}\leq C,\quad \big\|\frac{\varphi^\varepsilon-\varphi}{\varepsilon}\big\|_{r} \leq C,$$ where a generic constant $C>0$ is independent of $\varepsilon$. This implies that \begin{gather*} \varphi^\varepsilon \stackrel{s}\to \varphi \quad \text{in } L^r(\Omega),\\ \nabla\varphi^\varepsilon \stackrel{s}\to \nabla\varphi \quad \text{in } L^r(\Omega) \end{gather*} as $\varepsilon\to 0$. Moreover, by reiterating this argument we can show that for a.e. $t\in (0,T)$ $$\label{2r} \big\|\nabla\big(\frac{\varphi^\varepsilon-\varphi}{\varepsilon}\big) \big\|_{2r}\leq C,\quad \big\|\frac{\varphi^\varepsilon-\varphi}{\varepsilon}\big\|_{2r}\leq C.$$ This proves the claim. Now we can proceed to the derivation of sensitivity equations, starting with \label{sen} \begin{gathered} \begin{aligned} &\int_0^T\big\langle\frac{u_t^\varepsilon-u_t}{\varepsilon}, v\big\rangle\,dt +\int_Q\nabla \big(\frac{u^\varepsilon-u}{\varepsilon}\big)\nabla v\,dx\,dt+\int_{\Gamma_N^T} \beta\big(\frac{u^\varepsilon-u}{\varepsilon}\big)v\,ds\,dt\\ &=-\int_{\Gamma_N^T}\ell u^\varepsilon v\,ds\,dt +\frac{1}{\varepsilon}\int_Q \Big[(\varphi_{0}-\varphi^\varepsilon) \sigma(u^\varepsilon) \nabla\varphi^\varepsilon-(\varphi_{0}-\varphi) \sigma(u) \nabla\varphi\Big]\cdot\nabla v\,dx\,dt\\ &\quad +\frac{1}{\varepsilon}\int_Q\Big[\sigma(u^\varepsilon)\nabla \varphi^\varepsilon-\sigma(u)\nabla \varphi\Big]\cdot\nabla\varphi_{0}v\,dx\,dt, \end{aligned}\\ \frac{1}{\varepsilon}\int_\Omega\Big[\sigma(u^\varepsilon)\nabla \varphi^\varepsilon -\sigma(u)\nabla\varphi\Big]\cdot\nabla w\,dx\,dt=0. \end{gathered} The following convergence are immediate \begin{gather*} \int_0^T\big\langle\frac{u_t^\varepsilon-u_t}{\varepsilon},v\big\rangle\,dt \to \int_0^T(\psi_1)_t,v\Big\rangle\,dt,\\ \int_Q\nabla \big(\frac{u^\varepsilon-u}{\varepsilon}\big)\nabla v\,dx\,dt \to \int_Q \nabla\psi_1\cdot\nabla v\,dx\,dt,\\ \int_{\Gamma_N^T} \beta\big(\frac{u^\varepsilon-u}{\varepsilon}\big)v\,ds\,dt \to \int_{\Gamma_N^T} \beta\psi_1 v\,ds\,dt,\\ \int_{\Gamma_N^T}\ell u^\varepsilon v\,ds\,dt \to \int_{\Gamma_N^T}\ell u v\,ds\,dt. \end{gather*} We write the second term on the right hand side of (\ref{sen}) in the form \begin{align*} & \frac{1}{\varepsilon}\int_Q \Big[(\varphi_{_0}-\varphi^\varepsilon) \sigma(u^\varepsilon) \nabla\varphi^\varepsilon\cdot\nabla v-(\varphi_{_0}-\varphi) \sigma(u) \nabla\varphi\cdot\nabla v\Big]\,dx\,dt\\ &= \frac{1}{\varepsilon}\int_Q (\varphi_{_0}-\varphi^\varepsilon)\Big[\sigma(u^\varepsilon)\nabla\varphi^\varepsilon- \sigma(u)\nabla\varphi\Big]\cdot\nabla v\,dx\,dt\\ &\quad +\frac{1}{\varepsilon}\int_Q (\varphi-\varphi^\varepsilon)\sigma(u)\nabla\varphi\cdot\nabla v\,dx\,dt\\ &= \mathcal{G}+\mathcal{F}. \end{align*} As far as the term $\mathcal{F}$ is concerned, we have that $\sigma(u)\nabla\varphi\cdot\nabla v\in L^2(Q)$ and $(\varphi-\varphi^\varepsilon)/\varepsilon\stackrel{w}\rightharpoonup -\psi_2$ in $L^2(Q)$, and therefore \begin{align*} \mathcal{F}\to -\int_Q \psi_2\sigma(u)\nabla\varphi\cdot\nabla v\,dx\,dt \end{align*} as $\varepsilon\to 0$. Now we consider terms coming from $\mathcal{G}$. We write \begin{align*} \mathcal{G} &= \frac{1}{\varepsilon}\int_Q (\varphi_{_0}-\varphi^\varepsilon)\Big[\sigma(u^\varepsilon) \nabla\varphi^\varepsilon- \sigma(u)\nabla\varphi\Big]\cdot\nabla v\,dx\,dt\\ &= \frac{1}{\varepsilon}\int_Q(\varphi_{_0}-\varphi^\varepsilon) \Big[\sigma(u^\varepsilon)-\sigma(u)\Big] \nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt\\ &\quad+\frac{1}{\varepsilon}\int_Q(\varphi_{_0}-\varphi^\varepsilon) \sigma(u)\Big[\nabla\varphi^\varepsilon-\nabla\varphi\Big]\cdot\nabla v\,dx\,dt =\mathcal{G}_1+\mathcal{G}_2 \end{align*} Using (\ref{2r}) and the fact that $\nabla(\varphi^\varepsilon-\varphi)/\varepsilon\stackrel{w}\rightharpoonup\nabla\psi_2$ in $L^2(Q)$, one can show $$\label{G2} \mathcal{G}_2\to \int_Q(\varphi_{0}-\varphi) \sigma(u)\nabla\psi_2\cdot\nabla v\,dx\,dt\quad\text{as } \varepsilon\to 0.$$ To illustrate terms from $\mathcal{G}_1$, observe that \begin{align*} \frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\stackrel{s} \to \sigma'(u)\psi_1\quad\text{in }L^2(0,T;L^s(\Omega))\quad\text{as } \varepsilon\to 0. \end{align*} Indeed, since $N=2$ we have $H_0^1(\Omega)\subset\subset L^q(\Omega)$ for any $q\in [1,\infty)$, $L^2(0,T;H_0^1(\Omega))\subset L^2(0,T;L^s(\Omega))$, and $$\label{alpha0} \big\|\frac{u^\varepsilon-u}{\varepsilon} \big\|_{L^2(0,T;L^s(\Omega))}\leq c\big\|\frac{u^\varepsilon-u}{\varepsilon} \big\|_{L^2(0,T;H_0^1(\Omega))}\leq C.$$ Therefore, $$\label{alpha} \big\|\frac{u_t^\varepsilon-u_t}{\varepsilon} \big\|_{L^2(0,T;H^{-1}(\Omega))}\leq C,$$ where (\ref{alpha}) can be derived from the PDE. Since $\frac{1}{s}+\frac{1}{r}+\frac{1}{2}=1$ and $r>2$, we have that $s=\frac{2r}{r-2}>2$. If we define $\tilde W:=\{v: v\in L^2(0,T;L^s(\Omega)), v_t\in L^2(0,T;H^{-1}(\Omega)) \}$, then, because $H_0^1(\Omega)\subset\subset L^s(\Omega)\subset H^{-1}(\Omega)$, it follows from \cite{Sim} that $\tilde W\subset\subset L^2(0,T;L^s(\Omega))$. Hence, the estimates (\ref{alpha}) imply that on a subsequence $\big\|\frac{u^\varepsilon-u}{\varepsilon} -\psi_1\big\|_{L^2(0,T;L^s(\Omega))}\to 0\quad\text{as } \varepsilon\to 0,$ and therefore, $$\label{sigma} \big\|\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon} -\sigma'(u)\psi_1\big\|_{L^2(0,T;L^s(\Omega))}\to 0\quad\text{as } \varepsilon\to 0.$$ Now, we are ready show \begin{align*}%\label{G1} \mathcal{G}_1&\stackrel{\rm def}=\int_Q(\varphi_{0}-\varphi^\varepsilon) \big(\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\big) \nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt\\ & \to \int_Q(\varphi_{0}-\varphi) \sigma'(u)\psi_1\nabla\varphi\cdot\nabla v\,dx\,dt. \end{align*} To show this convergence, we estimate \begin{align*} & \big|\int_Q(\varphi_{0}-\varphi^\varepsilon) \big(\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\big) \nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt-\int_Q(\varphi_{0}-\varphi) \sigma'(u)\psi_1\nabla\varphi\cdot\nabla v\,dx\,dt\big|\\ &\leq \big|\int_Q(\varphi_{0}-\varphi)\Big\{ \big(\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\big) \nabla\varphi^\varepsilon-\sigma'(u)\psi_1\nabla\varphi\Big\} \cdot\nabla v\,dx\,dt\big|\\ &\quad+\big|\int_Q(\varphi-\varphi^\varepsilon) \big(\frac{\sigma(u^\varepsilon) -\sigma(u)}{\varepsilon}-\sigma'(u)\psi_1+ \sigma'(u)\psi_1\big)\nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt\big|\\ &\leq \big|\int_Q\Big\{(\varphi_{0}-\varphi) \Big[\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}-\sigma'(u)\psi_1\Big] \nabla\varphi^\varepsilon\cdot\nabla v\\ &\quad+ (\varphi_{0}-\varphi) \sigma'(u)\psi_1\nabla(\varphi^\varepsilon-\varphi)\cdot\nabla v\Big\}\,dx\,dt\big|\\ &\quad+\big|\int_Q(\varphi-\varphi^\varepsilon) \big(\frac{\sigma(u^\varepsilon) -\sigma(u)}{\varepsilon}-\sigma'(u)\psi_1\big) \nabla\varphi^\varepsilon\cdot\nabla v\\ &\quad + (\varphi-\varphi^\varepsilon) \sigma'(u)\psi_1\nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt\big|\\ &\leq \big|\int_Q(\varphi_{0}-\varphi)\Big\{ \big(\frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon}\big) -\sigma'(u)\psi_1\Big\}\nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt\big|\\ &\quad+\big|\int_Q(\varphi_{0}-\varphi) \sigma'(u)\psi_1\nabla(\varphi^\varepsilon-\varphi)\cdot\nabla v\,dx\,dt\big|\\ &\quad+\big|\int_Q(\varphi-\varphi^\varepsilon) \big(\frac{\sigma(u^\varepsilon) -\sigma(u)}{\varepsilon}-\sigma'(u)\psi_1\big)\nabla\varphi^\varepsilon\cdot\nabla v \,dx\,dt\big|\\ &\quad+\big|\int_Q(\varphi-\varphi^\varepsilon) \sigma'(u)\psi_1\nabla\varphi^\varepsilon\cdot\nabla v\,dx\,dt\big|:=I+II+III+IV. \end{align*} We deal with the term $I$ first. We have \begin{align*} I&\leq \|\varphi-\varphi_0 \|_{L^\infty(Q)}\int_0^T\int_\Omega\big| \frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon} -\sigma'(u)\psi_1\big|\cdot |\nabla\varphi^\varepsilon|\cdot|\nabla v|\,dx\,dt\\ &\leq \|\varphi-\varphi_0 \|_{L^\infty(Q)}\int_0^T\big\| \frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon} -\sigma'(u)\psi_1\big\|_s\cdot \|\nabla\varphi^\varepsilon\|_r\cdot\|\nabla v\|_2\,dt\\ &\leq \|\varphi-\varphi_0 \|_{L^\infty(Q)}\sup_{t}{\| \nabla\varphi^\varepsilon\|_r}\int_0^T\big\| \frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon} -\sigma'(u)\psi_1\big\|_s\cdot\|\nabla v\|_2\,dt\\ &\leq C\big\| \frac{\sigma(u^\varepsilon)-\sigma(u)}{\varepsilon} -\sigma'(u)\psi_1\big\|_{L^2(0,T;L^s(\Omega))}\cdot\|\nabla v\|_{L^2(Q)}\to 0\quad\text{as}\ \varepsilon\to 0, \end{align*} by (\ref{sigma}). Similarly, we obtain \begin{align*} II&:=\big|\int_0^T\int_\Omega(\varphi_{0}-\varphi) \sigma'(u)\psi_1 \nabla(\varphi^\varepsilon-\varphi)\cdot\nabla v\,dx\,dt\big|\\ &\leq C\int_0^T\|\psi_1\|_s\cdot\|\nabla(\varphi^\varepsilon-\varphi)\|_r\cdot\|\nabla v\|_2\,dt\\ &\leq C\|\psi_1\|_{L^2(0,T;L^s(\Omega))}\cdot\|\nabla v\|_{L^2(Q)}\cdot\sup_{t}{\|\nabla(\varphi^\varepsilon-\varphi)\|_r}\to 0\quad\text{as}\ \varepsilon\to 0, \end{align*} where we used (\ref{r}). Analogously, it can be shown that $III, IV\to 0$ as $\varepsilon\to 0$. One can show that the third term on the right hand side of (\ref{sen}) satisfies \begin{align*}%\label{lastterm} & \frac{1}{\varepsilon}\int_Q\Big[\sigma(u^\varepsilon)\nabla \varphi^\varepsilon-\sigma(u)\nabla \varphi\Big]\cdot\nabla\varphi_{0}v\,dx\,dt\\ &\to \int_Q\sigma'(u)\psi_1\nabla\varphi\cdot\nabla\varphi_{0}v\,dx\,dt +\int_Q\sigma(u)\nabla\psi_2\cdot\nabla\varphi_{0}v\,dx\,dt\quad \text{as } \varepsilon\to 0. \end{align*} Finally, the equation for $\varphi$ in (\ref{sen}) can be shown to satisfy \begin{align*}%\label{sensitivityphi} & \frac{1}{\varepsilon}\int_Q\Big[\sigma(u^\varepsilon)\nabla \varphi^\varepsilon\cdot\nabla w -\sigma(u)\nabla\varphi\cdot\nabla w\Big]\,dx\,dt\\ &\to \int_Q\sigma'(u)\psi_1\nabla \varphi\cdot\nabla w\,dx\,dt +\int_Q\sigma(u)\nabla\psi_2\cdot\nabla w\,dx\,dt \quad\text{as } \varepsilon\to 0. \end{align*} Letting $\varepsilon\to 0$ in (\ref{sen}), we obtain \label{sensitweakformul} \begin{gathered} \begin{aligned} & \int_0^T\big\langle(\psi_1)_t,v\big\rangle\,dt+\int_Q\nabla \psi_1\nabla v\,dx\,dt+\int_{\Gamma_N^T} (\beta\psi_1+\ell u)v\,ds\,dt\\ &= \int_Q(\varphi_{0}-\varphi) \sigma'(u)\psi_1 \nabla\varphi\cdot\nabla v\,dx\,dt+\int_Q(\varphi_{0}-\varphi) \sigma(u)\nabla\psi_2 \cdot\nabla v\,dx\,dt\\ &\quad - \int_Q\psi_2\sigma(u)\nabla\varphi\cdot\nabla v\,dx\,dt+ \int_Q\sigma'(u)\psi_1\nabla\varphi\cdot\nabla\varphi_0 v\,dx\,dt\\ &\quad+\int_Q\sigma(u)\nabla\psi_2\cdot\nabla\varphi_0 v\,dx\,dt, \end{aligned}\\ \int_Q(\sigma'(u)\psi_1\nabla\varphi+\sigma(u)\nabla\psi_2)\cdot \nabla w\,dx\,dt=0 \end{gathered} for all $v\in L^2(0,T;H^{1}(\Omega))\cap L^\infty(Q)$ and $w\in L^2(0,T;H^{1}_{0}(\Omega))$. This leads to the sensitivity system (\ref{sensitivityeqn}). The weak formulation for (\ref{sensitivityeqn}) is given by \begin{gather} \label{weakform1} \begin{aligned} & \int_0^T\big\langle-(\psi_1)_t,v\big\rangle\,dt +\int_Q\nabla\psi_1\cdot\nabla v\,dx\,dt +\int_{\Gamma_N^T}(\beta\psi_1+\ell u) v\,ds\,dt\\ &= \int_Q\sigma'(u)\psi_1|\nabla\varphi|^2v\,dx\,dt+2\int_Q \sigma(u)\nabla\varphi\cdot\nabla\psi_2 v\,dx\,dt, \end{aligned}\\ \int_Q(\sigma'(u)\psi_1\nabla\varphi+\sigma(u)\nabla\psi_2)\cdot\nabla w\,dx\,dt=0 \label{weakform2} \end{gather} for all $v\in L^2(0,T;H^{1}(\Omega))\cap L^\infty(Q)$ and $w\in L^2(0,T;H^{1}_{0}(\Omega))$ which can be shown to be equivalent to (\ref{sensitweakformul}). \end{proof} To characterize the optimal control, we need to introduce adjoint functions and the adjoint operator associated with $\psi_1$ and $\psi_2$. Using $\frac{\partial p}{\partial n}+\beta^*p=0$ on $\Gamma_N\times(0,T)$, $q=\psi_2=0$ on $\partial\Omega\times(0,T)$, $p=0$ on $\Gamma_D\times(0,T)$, $p(x,T)=0$, and the notation $\mathcal{L}$ for the operator in the sensitivity system gives: \begin{align*} & \int_Q \begin{pmatrix} p & q\end{pmatrix} \mathcal{L} \begin{pmatrix} \psi_1\\ \psi_2 \end{pmatrix} \,dx\,dt\\ &=\int_Q p\,(-\psi_1)_t\,dx\,dt+\int_Q p\,\Delta\psi_1\,dx\,dt + \int_Q p \sigma'(u) |\nabla\varphi|^2\psi_1\,dx\,dt\\ &\quad +2\int_Q p \sigma(u)\nabla\varphi\cdot\nabla\psi_2\,dx\,dt + \int_Q q\nabla\cdot[\sigma(u)\nabla\psi_2]\,dx\,dt + \int_Q q\nabla\cdot[\sigma'(u)\psi_1\nabla\varphi]\,dx\,dt\\ &=\int_Q \psi_1 p_t\,dx\,dt+\int_Q \psi_1\Delta p\,dx\,dt +\int_Q\psi_1\sigma'(u) |\nabla\varphi|^2 p\,dx\,dt\\ &\quad -2\int_Q\psi_2\nabla\cdot[p \sigma(u)\nabla\varphi]\,dx\,dt -\int_Q \psi_1\sigma'(u)\nabla\varphi\cdot\nabla q\,dx\,dt\\ &\quad +\int_Q\psi_2\nabla\cdot[\sigma(u)\nabla q]\,dx\,dt\\ &= \int_Q \begin{pmatrix} \psi_1&\psi_2 \end{pmatrix} \mathcal{L}^* \begin{pmatrix} p\\ q \end{pmatrix} \,dx\,dt, \end{align*} where $%\label{operatorLstar} \mathcal{L}^* \begin{pmatrix} p\\ q \end{pmatrix} = \begin{pmatrix} p_t+\Delta p+\sigma'(u)|\nabla\varphi|^2p-\sigma'(u) \nabla\varphi\cdot\nabla q\\ \nabla\cdot[-2p \sigma(u)\nabla\varphi+\sigma(u)\nabla q] \end{pmatrix}.$ Thus, the adjoint system is given by $$\label{adjoint1} \begin{gathered} p_t+\Delta p+\sigma'(u)|\nabla\varphi|^2p-\sigma'(u) \nabla\varphi\cdot\nabla q = 1\quad\text{in }\Omega\times(0,T),\\ \nabla\cdot[-2p \sigma(u)\nabla\varphi+\sigma(u)\nabla q] = 0\quad\text{in }\Omega\times(0,T),\\ \frac{\partial p}{\partial n}+\beta^* p = 0\quad\text{on } \Gamma_N\times(0,T),\\ p= 0\quad\text{on }\Gamma_D\times(0,T),\\ q= 0\quad\text{on }\partial\Omega\times(0,T),\\ p= 0\quad\text{on }\Omega\times\{T\}, \end{gathered}$$ where the nonhomogeneous term 1'' comes from differentiating the integrand of $J(\beta)$ with respect to the state $u$. \begin{theorem}\label{adjoint_theorem} Let $\|\varphi_0 \|_{W^{1,\infty}(Q)}$ be sufficiently small. Then, given an optimal control $\beta^*\in U_M$ and corresponding states $u,\varphi$, there exists a solution $(p,q)\in L^2(0,T;H_0^1(\Omega))\times L^2(0,T;H_0^1(\Omega))$ to the adjoint system \eqref{adjoint1}. Furthermore, $\beta^*$ can be explicitly characterized as: $$\label{characterization2} \beta^*(x,t)=\min\Big(\max\big(-\frac{up}{2},\lambda\big),M\Big).$$ \end{theorem} \begin{proof} The weak formulation of \eqref{adjoint1} is $$\label{adjointweak} \begin{gathered} \int_0^T\langle p_t,v\rangle\,dt -\int_Q\nabla p\cdot\nabla v\,dx\,dt -\int_{\Gamma_N\times(0,T)}\beta^* p\,v\,ds\,dt\\ +\int_Q\sigma'(u)|\nabla\varphi|^2p\,v\,dx\,dt -\int_Q\sigma'(u)(\nabla\varphi\cdot\nabla q)v\,dx\,dt =\int_Q v\,dx\,dt, \\ 2\int_Q p \sigma(u)\nabla\varphi\cdot\nabla w\,dx\,dt-\int_Q\sigma(u)\nabla q\cdot\nabla w\,dx\,dt=0 \end{gathered}$$ for all $v\in L^2(0,T;H^{1}(\Omega))\cap L^\infty(Q)$ and $w\in L^2(0,T;H^{1}_{0}(\Omega))$. We show that the solution to the adjoint system \eqref{adjoint1} exists. Define the map $F:L^2(0,T;H^1(\Omega))\times L^2(0,T;H_0^1(\Omega))\mapsto L^2(0,T;H^1(\Omega))\times L^2(0,T;H_0^1(\Omega))$ with $F(V,W)=(p,q)$ as follows: $$\label{adjoint2} \begin{gathered} p_t+\Delta p+\sigma'(u)|\nabla\varphi|^2V-\sigma'(u) \nabla\varphi\cdot\nabla W = 1\quad\text{in }\Omega\times(0,T),\\ \nabla\cdot[-2V \sigma(u)\nabla\varphi+\sigma(u)\nabla q] = 0\quad\text{in }\Omega\times(0,T),\\ \frac{\partial p}{\partial n}+\beta^* p = 0\quad\text{on } \Gamma_N\times(0,T),\\ p = 0\quad\text{on }\Gamma_D\times(0,T),\\ q = 0\quad\text{on }\partial\Omega\times(0,T),\\ p = 0\quad\text{on }\Omega\times\{T\}, \end{gathered}$$ The weak formulation of \eqref{adjoint2} is \begin{gather*} \int_0^T\langle p_t,\Theta\rangle\,dt -\int_Q\nabla p\cdot\nabla\Theta\,dx\,dt -\int_{\Gamma_N^T}\beta^* p\Theta\,ds\,dt\\ +\int_Q\sigma'(u)|\nabla\varphi|^2V\Theta\,dx\,dt -\int_Q\sigma'(u) (\nabla\varphi\cdot\nabla W) \Theta\,dx\,dt =\int_Q\Theta\,dx\,dt,\\ 2\int_Q V \sigma(u)\nabla\varphi\cdot\nabla \Psi\,dx\,dt=\int_Q\sigma(u)\nabla q\cdot\nabla \Psi\,dx\,dt \end{gather*} for all $\Theta\in L^2(0,T;H^{1}(\Omega))\cap L^\infty(Q)$ and $\Psi\in L^2(0,T;H^{1}_{0}(\Omega))$. We use the Banach fixed point theorem. Let $F(V_1,W_1)=(p_1,q_1)$ and $F(V_2,W_2)=(p_2,q_2)$. We show the contraction property $$\label{Banach} \|p_1-p_2\|_{H^1}+\|q_1-q_2\|_{H_0^1}\leq \delta\Big(\|V_1-V_2\|_{H^1}+\|W_1-W_2\|_{H_0^1}\Big)$$ for some $0<\delta<1$. Indeed, for $F(V_1,W_1)=(p_1,q_1)$, we have \begin{gather*} %\label{adjointweak2} \int_0^T\langle (p_1)_t,\Theta\rangle\,dt-\int_Q\nabla p_1\cdot\nabla\Theta\,dx\,dt-\int_{\Gamma_N^T}\beta^* p_1 \Theta\,ds\,dt\\ +\int_Q\sigma'(u)|\nabla\varphi|^2V_1 \Theta\,dx\,dt -\int_Q\sigma'(u) (\nabla\varphi\cdot\nabla W_1) \Theta\,dx\,dt=\int_Q\Theta\,dx\,dt,\\ 2\int_Q V_1 \sigma(u)\nabla\varphi\cdot\nabla \Psi\,dx\,dt=\int_Q\sigma(u)\nabla q_1\cdot\nabla \Psi\,dx\,dt, \end{gather*} and similarly, for $F(V_2,W_2)=(p_2,q_2)$. Let us denote $\bar p:=p_1-p_2$, $\bar q:=q_1-q_2$, $\bar W:=W_1-W_2$, and $\bar V:=V_1-V_2$. Take $\Theta=\bar p$, then equations for $p_1$ and $p_2$ give \begin{align*}%\label{pexistence2} &\int_0^T\langle (\bar p)_t,\bar p\rangle\,dt+\int_Q\sigma'(u)|\nabla\varphi|^2\bar V\bar p\,dx\,dt-\int_Q\sigma'(u) \nabla\varphi\cdot\nabla\bar W\bar p\,dx\,dt\\ &=\int_Q|\nabla\bar p|^2dx\,dt+\int_{\Gamma_N^T}\beta^* \bar p^2\,ds\,dt.\\ \end{align*} Thus, we obtain \begin{align*}%\label{pexistence3} & \frac{1}{2}\int_{\Omega\times\{0\}}\bar{p}^2\,dx+\int_Q|\nabla\bar p|^2\,dx\,dt+\lambda\int_{\Gamma_N^T}\bar p^2\,ds\,dt\\ & \leq\Big|\int_Q\sigma'(u)|\nabla\varphi|^2\,\bar V\bar p\,dx\,dt\Big|+\Big|\int_Q\sigma'(u) \nabla\varphi\nabla\bar W\bar p\,dx\,dt\Big|. \end{align*} We have \begin{align*} \Big|\int_Q\sigma'(u)|\nabla\varphi|^2\,\bar V\bar p\,dx\,dt\Big| &\leq C\sup_{t}{\||\nabla\varphi|^2\|_r}\int_0^T\|\bar V\|_s\cdot\|\bar p\|_2\,dt\\ & \leq C_{\varphi}\|\bar V\|_{L^2(0,T;L^s(\Omega))}\|\bar p\|_{L^2(Q)}\\ &\leq C_{\varphi}\|\bar V\|_{L^2(0,T;L^s(\Omega))}\|\bar p\|_{L^2(0,T;H^1(\Omega))}, \end{align*} and \begin{align*} \Big|\int_Q\sigma'(u) \nabla\varphi\cdot\nabla\bar W\bar p\,dx\,dt\Big| \leq\tilde{C}_{\varphi}\|\nabla\bar W\|_{L^2(Q)}\|\bar p\|_{L^2(0,T;H^1(\Omega))}. \end{align*} The above three inequalities give $$\label{p} \|\bar p\|_{L^2(0,T;H^1(\Omega))}\leq C_{\varphi}\|\bar V\|_{L^2(0,T;H^1(\Omega))}+\tilde C_{\varphi}\|\bar W\|_{L^2(0,T;H^1(\Omega))},$$ where the constants $C_{\varphi}$ and $\tilde C_{\varphi}$ depend on the boundary data $\varphi_0$. A similar estimate can be derived for $\bar q$ \begin{eqnarray}\label{q} \|\bar q\|_{L^2(0,T;H^1(\Omega))}\leq C_{\varphi}'\|\bar V\|_{L^2(0,T;H^1(\Omega))}. \end{eqnarray} Adding (\ref{p}) and (\ref{q}) we get \begin{align*} \|\bar p\|_{L^2(0,T;H^1(\Omega))}+\|\bar q\|_{L^2(0,T;H^1(\Omega))}\leq c_{\varphi}(\|\bar V\|_{L^2(0,T;H^1(\Omega))}+\|\bar W\|_{L^2(0,T;H^1(\Omega))}). \end{align*} Now, choosing $\varphi_0$ so that $\|\varphi_0\|_{W^{1,\infty}(Q)}$ is small, we obtain $\delta<1$. This proves (\ref{Banach}), the contraction property, which gives the desired fixed point of the map and therefore the existence of solutions to the adjoint system. Recall that for variation $\ell\in L^\infty(\Gamma_N^T)$, with $\beta^*+\varepsilon\ell\in U_M$ for $\varepsilon$ small, the weak formulation of the sensitivity system (\ref{sensitivityeqn}) is given by equations (\ref{weakform1}) and (\ref{weakform2}). Now we are ready to characterize the optimal control. Since the minimum of $J$ is achieved at $\beta^*$ and for small $\varepsilon>0$, $\beta^*+\varepsilon\ell\in U_M$, and denoting $u^\varepsilon=u(\beta^*+\varepsilon\ell),\varphi^\varepsilon= \varphi(\beta^*+\varepsilon\ell)$, we obtain \begin{align*} 0&\leq \lim_{\varepsilon\to 0^+}\frac{J(\beta^*+\varepsilon\ell)-J(\beta^*)}{\varepsilon}\\ &= \int_Q\psi_1\,dx\,dt+\int_{\Gamma_N^T}2\beta^*\ell\,ds\,dt\\ &=\int_Q \begin{pmatrix} \psi_1&\psi_2 \end{pmatrix} \begin{pmatrix} 1\\ 0\end{pmatrix}\,dx\,dt +\int_{\Gamma_N^T}2\beta^*\ell\,ds\,dt\\ &= -\int_Q\langle(\psi_1)_t,p\rangle\,dt-\int_Q\nabla p\nabla \psi_1\,dx\,dt-\int_{\Gamma_N^T}\beta^* p\,\psi_1\,ds\,dt\\ &\quad +\int_Q\psi_1 \sigma'(u)|\nabla\varphi|^2p\,dx\,dt -\int_Q\psi_1 \sigma'(u)\nabla\varphi\cdot\nabla q\,dx\,dt\\ &\quad +2\int_Q p \sigma(u)\nabla\varphi\cdot\nabla \psi_2\,dx\,dt-\int_Q\sigma(u)\nabla q\cdot\nabla\psi_2\,dx\,dt +\int_{\Gamma_N^T}2\beta^*\ell\,ds\,dt \\ &=\Big\{-\int_Q\nabla p\cdot\nabla \psi_1\,dx\,dt-\int_{\Gamma_N^T}\beta^* p\,\psi_1\,ds\,dt\\ &\quad+\int_Q\psi_1 \sigma'(u)|\nabla\varphi|^2p\,dx\,dt +2\int_Q p \sigma(u)\nabla\varphi\cdot\nabla\psi_2\,dx\,dt \Big\}\\ &\quad+\Big[-\int_Q\psi_1 \sigma'(u)\nabla\varphi\cdot\nabla q\,dx\,dt -\int_Q\sigma(u)\nabla q\cdot\nabla\psi_2\,dx\,dt\Big] +\int_{\Gamma_N^T}2\beta^*\ell\,ds\,dt\\ &= \int_{\Gamma_N^T}\ell\,u p\,ds\,dt+\int_{\Gamma_N^T}2\beta^*\ell\,ds\,dt, \end{align*} where we integrated by parts and used the sensitivity system (\ref{weakform1}), (\ref{weakform2}) with test functions $p, q$. Thus, we obtain $$\label{variation} \int_{\Gamma_N^T}\ell(2\beta^*+up)\,ds\geq 0.$$ We use this inequality with three cases to characterize $\beta^*$: (i) Take the variation $\ell$ to have support on the set \$\{x\in\Gamma_N^T:\lambda<\beta^*(x,t)