\documentclass[reqno]{amsart}
\usepackage{a4wide}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2026 (2026), No. 18, pp. 1--26.\newline
ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, DOI: 10.58997/ejde.2026.18}
\thanks{\copyright 2026. This work is licensed under a CC BY 4.0 license.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2026/18\hfil New cost terms through homogenization]
{New cost terms through homogenization of an optimal control problem
under dynamic boundary conditions on microscopic particles}


\author[J. I. D\'{\i}az, T. A. Shaposhnikova, A. V. Podolskiy \hfil EJDE-2026/18\hfilneg]
{Jes\'us Ildefonso D\'{\i}az, Tatyana A. Shaposhnikova, Alexander V. Podolskiy}

\address{Jes\'us Ildefonso D\'{\i}az \newline
Instituto de Mathematica Interdisciplinar,
Universidad Complutense de Madrid, Spain}
\email{jidiaz@ucm.es}

\address{Tatyana A. Shaposhnikova \newline
Lomonosov Moscow State University, Russia}
\email{shaposh.tan@mail.ru}

\address{Alexander V. Podolskiy \newline
Lomonosov Moscow State University, Russia}
\email{avpodolskiy@yandex.ru}

\thanks{Submitted October 17, 2025. Published February 25, 2026.}
\subjclass[2020]{35B27, 35K20, 49K20, 93C20}
\keywords{Homogenization; critical case; optimal control; strange term;
\hfill\break\indent dynamic boundary condition; homogenized cost functional}

\begin{abstract}
This article concerns optimal control problems in a heterogeneous body with
a periodic structure of particles depending on a small parameter $\varepsilon$.
We study the asymptotic behavior, as $\varepsilon \to 0$, of the optimal
control functional and of the optimal state when the initial problem is of parabolic type.
We assume a dynamic condition and the effect of some controls for some of 
the particles on the boundary.
In the so-called ``critical case'', we show the appearance of some new
non-local in time ``strange terms'', in the limit parabolic equation
and in the limit cost functional.
Microscopic localized controls generate peculiar terms in both the limit
equation and the cost functional that do not appear when controls are
applied to the entire set of particles, or when the boundary condition on
the particles is of Robin type.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

The problem we consider here arises in many fields. For instance,
it is well-known that many problems in chemical engineering lead to the
optimization of some cost functionals
\cite{Nolasco-Survey, Upreti book}. The same thing happens in
porous media theory in which the word ``particle''  must be replaced by
``perforation''  (see, e.g.
\cite{ Anguiano, ConcaDliTi,  GoSanchezP, Hornung-Yaeger, Hurlov,  
Iliev-Mikelik, OlSh95, Timof} and many other
references quoted in the monograph \cite{DiGoShBook}).
Simplified models in Climatology can be also modeled in terms very close to the
ones we will study in this paper (see \cite{DiPoShRACSAM}).

The main goal of this article is to illustrate how the homogenization of some
optimal control problems may give rise to new non-local in time ``strange
terms''. This also happens in the limit parabolic equation and in the limit
cost functional, assuming a dynamic boundary condition and the actuation of some
controls on some subset of the particles.
We will do that for the so-called
``critical case'', that is characterized
by the relation between the structure's period, the diameter of the
balls, and the growth coefficient in the particles' boundary condition.
In this way, microscopic localized controls generate peculiar terms in both
the limit equation and the cost function that do not appear, for instance, in
the case of the Robin type boundary condition on the particles.

 We give a detailed presentation of the heterogeneous domain $\Omega
_{\varepsilon }$ in the next section. At the moment, we outline that
since in the most of the cases it is impossible to act over the
entire spatial domain $\Omega _{\varepsilon }$, the control is applied only
on the boundary of the particles contained in a small portion of the domain
($\omega $ such that $\overline{\omega }\subset \Omega $).
Thus, the set of boundaries of the internal particles is constituted in the form
$S_{\varepsilon }=S_{\varepsilon }^{1}\cup S_{\varepsilon }^2$, where
$S_{\varepsilon }^2$ is the set of boundaries of the controlling particles
$G_{\varepsilon }^2$ and $S_{\varepsilon }^{1}$ is the set of boundaries of
the particles $G_{\varepsilon }^{1}$ to which no control is implemented. The
state of the control problem is given through
\begin{equation}
\begin{gathered}
\partial _{t}u_{\varepsilon }(v)-\Delta u_{\varepsilon }(v)=f, \quad
 (x,t)\in Q_{\varepsilon }^{T}, \\
\varepsilon ^{-\gamma }\partial _{t}u_{\varepsilon }(v)+\partial _{\nu
}u_{\varepsilon }(v)=\varepsilon ^{-\gamma }\chi _{S_{\varepsilon }^{2,T}}v,
\quad (x,t)\in S_{\varepsilon }^{T}, \\
u_{\varepsilon }(v)(x,0)=0, \quad x\in \Omega _{\varepsilon }\cup
S_{\varepsilon }, \\
u_{\varepsilon }(v)(x,t)=0, \quad (x,t)\in \Gamma ^{T},%
\end{gathered}  \label{init state prob}
\end{equation}
where $f\in L^2(Q^{T})$ and\textrm{\ }$v\in L^2(S_{\varepsilon }^{2,T})$
is the control.\textrm{\ }Here, we are using the notation
(considering $0<T<\infty $)
\begin{equation}
\begin{gathered}
\Omega _{\varepsilon }=\Omega \setminus \overline{G_{\varepsilon }}, \quad
S_{\varepsilon }=\partial G_{\varepsilon }, \quad
\partial \Omega _{\varepsilon }=S_{\varepsilon }\cup \partial \Omega ,   \quad
Q_{\varepsilon }^{T}=\Omega _{\varepsilon }\times (0,T), \\
\Gamma^{T}=\partial \Omega \times (0,T), \quad
S_{\varepsilon }^{T}=S_{\varepsilon }\times (0,T), \quad Q^{T}=\Omega \times (0,T),
\end{gathered}
\label{def: basic sets}
\end{equation}
which will be described in detail in the next section. We note that
$G_{\varepsilon }$ is the set of small particles
($\varepsilon $-periodically distributed and homothetic to a unit ball $G_0$)
in an open bounded regular set $\Omega $ of $\mathbb{R}^n$, $n\geq 3$.
By $\chi_{S_{\varepsilon }^{2,T}}$, we denote the characteristic function of the set
$S_{\varepsilon }^{2,T}=S_{\varepsilon }^2\times (0,T)$ that lies entirely
in the set $\omega _{\varepsilon }^{T}$ defined below
\[
\omega _{\varepsilon }=\omega \cap \Omega _{\varepsilon },\quad
\omega_{\varepsilon }^{T}=\omega _{\varepsilon }\times (0,T),\quad
\omega^{T}=\omega \times (0,T).
\]
The parameter $\gamma >0$ plays a crucial role since in this paper we
consider the so-called ``critical case''
governed by the size of particles that are translations of a small particle
$a_{\varepsilon }G_0$, where $G_0$ is the unit ball with radius
$a_{\varepsilon }=C_0\varepsilon ^{\gamma }$, $\gamma =\frac{n}{n-2}$, and
$C_0$ is some positive constant.

\begin{figure}[htb]
\centering
 \includegraphics[width=0.5\textwidth]{fig1}
% perforated_by_balls_domain_partly_v2.png
\caption{Perforated domain. The control is implemented only on
$S_{\varepsilon }^2$, the boundary of some of the internal balls: the ones
collected under the $G^2_\varepsilon$.}
\label{fig: perforated domain}
\end{figure}

Notice that since problem \eqref{init state prob} is linear, by some
obvious change of variable, we can also consider the case of a non-zero
initial datum. To finalize the statement of the optimal control problem, we
introduce the cost functional
$J_{\varepsilon}:L^2(0,T;L^2(S_{\varepsilon }^2))\to \mathbb{R}$,
\begin{equation}
\begin{aligned}
&J_{\varepsilon }(v) \\
&=\frac{1}{2}\| \nabla {u_{\varepsilon }(v)}\|_{L^2(Q_{\varepsilon }^{T})}^2
 +\frac{1}{2}\int _{\Omega_{\varepsilon }}u_{\varepsilon }^2(v)(x,T)dx
+\frac{{\varepsilon }^{-\gamma }}{2}\int _{S_{\varepsilon}}u_{\varepsilon }^2(v)(x,T){ds}
 +\varepsilon ^{-\gamma }\frac{N}{2}\| {v}\| _{L^2(S_{\varepsilon }^{2,T})}^2,
\end{aligned}  \label{cost functional}
\end{equation}
where $N>0$. Then, the optimal control $v_{\varepsilon }$ is
\begin{equation}
J_{\varepsilon }(v_{\varepsilon })
=\inf _{v\in L^2(0,T;L^2(S_{\varepsilon }^2))}J_{\varepsilon }(v).
 \label{optim control def}
\end{equation}
In what follows, we will abuse the notation and simply write
$u_{\varepsilon}$ instead of $u_{\varepsilon }(v_{\varepsilon })$.
By applying different results in the literature
(see, e.g., \cite{LionsOpt,Fursikov2000OptimalCO, Trolstz, Glow-Lions}),
it is well-known that there exists a unique optimal control
$v_{\varepsilon }\in L^2(0,T;L^2(S_{\varepsilon }^2))$.

We point out that the consideration of a non-zero target function
$u_{T}\in L^2(\Omega _{\varepsilon })$, a given profile observed at the final time
$T $, can be reduced to the above case of $u_{T}\equiv 0$ by a suitable
change of variables, at least for a dense set of $u_{T}$ in $L^2(\Omega )$
(see Remark \ref{Rem target nonzero} below).

The main goal of this article is to apply a homogenization process to the
above optimal control problem when $\varepsilon \to 0$. As in many
other formulations, the kind of the limit problem strongly depends on the
size of the particles' radii $C_0\varepsilon ^{\alpha }$, $C_0>0$
(see, e.g. \cite{ Anguiano,Timof} and the exposition made in
\cite{DiGoShBook}). Here, we consider the critical case in which
$\alpha =\gamma ={n}/{(n-2)}$. For different elliptic and parabolic problems,
it is well-known that this critical choice leads to the emergence of a new local
``strange''  term (the naming is due to \cite{Cioranescu1997AST}) in the effective
partial differential equation
(see \cite{ Cioranescu1997AST, DiGoShBook,  Hurlov, Kaizu1,ZuShDiffEq}).

It is well-known that the introduction of a dynamic boundary condition on
the particle boundary causes the aforementioned ``strange''  term  to become a ``
non-local''  operator, derived by solving a suitable
ordinary differential equation (we refer to
\cite{ShDyn19, DiShZarbitrary}, for the case of an elliptic Poisson equation
for the state).
Also, it was shown that in the framework of optimal control problems
there appear some new terms in the limit cost functional (in contrast with
previous results in the literature for related formulations, (see, e.g.,
\cite{DShPo-first control,  DiPoSh control arbt size, DiPoShBound22,
DiPoShAttouch,  PoShControl20, SJPZ02, ZuShDynDok19,  Shaposhnikova2022HomogenizationOT}
for the case of distributed controls appearing in the state equation).
One of the major new features we will demonstrate in this article
is that when the controls act on the boundary of some particles, some new
terms appear in the cost functional, the non-local terms in time are of a
different nature and some new non-local in time operators must be introduced.

To state the homogenization results, we need to introduce several auxiliary
problems. On the uncontrolled particles, we use non-local operator
$M(\varphi )$, arising in previous studies (see \cite{DiPoShAttouch}), that
is defined as a solution to
\begin{equation}
\begin{gathered}
\partial _{t}M(\varphi )+\mathcal{B}_nM(\varphi )=\varphi , \quad t\in (0,T), \\
M(\varphi )(0)=0,
\end{gathered}
\end{equation}
where $\mathcal{B}_n=(n-2)C_0^{-1}$ and $\varphi \in L^2(0,T)$ is a
given function, and its adjoint operator $M^{\ast }$,
\begin{equation}
\begin{gathered}
-\partial _{t}M^{\ast }(\varphi )+\mathcal{B}_nM^{\ast }(\varphi )
=\varphi, \quad t\in (0,T), \\
M^{\ast }(\varphi )(T)=0.
\end{gathered}
\end{equation}

A similar non-local operator, $G(\varphi )$, and its adjoint operator
$G^{\ast }$, must be defined on the controlled particles
\begin{equation}
\begin{gathered}
\partial _{t}G(\varphi )+(\mathcal{B}_n+N^{-1})G(\varphi )=\varphi , \quad
t\in (0,T), \\
G(\varphi )(0)=0,
\end{gathered}  \label{G prob}
\end{equation}
and
\begin{equation}
\begin{gathered}
-\partial _{t}G^{\ast }(\varphi )+(\mathcal{B}_n+N^{-1})G^{\ast }(\varphi
)=\varphi , \quad t\in (0,T), \\
G^{\ast }(\varphi )(T)=0. 
\end{gathered}  \label{adj G prob}
\end{equation}
Besides that, we will need to define some new operators $H$ and $H^{\ast }$,
coupled with $G^{\ast }$ and $G$, respectively, by the problems
\begin{equation}
\begin{gathered}
-\partial _{t}H^{\ast }(\varphi )+(\mathcal{B}_n+N^{-1})H^{\ast }(\varphi
)-N^{-1}(\mathcal{B}_n+N^{-1})G(H^{\ast }(\varphi ))=\varphi , \quad
 t\in (0,T), \\
H^{\ast }(\varphi )(T)=0, 
\end{gathered} \label{H*eq}
\end{equation}
and
\begin{equation}
\begin{gathered}
\partial _{t}H(\varphi )+(\mathcal{B}_n+N^{-1})H(\varphi )
-N^{-1}(\mathcal{B}_n+N^{-1})G^{\ast }(H(\varphi ))=\varphi , \quad t\in (0,T), \\
H(\varphi )(0)=0. 
\end{gathered}  \label{eqH}
\end{equation}

Notice that the operators $G$, $M$, $G^{\ast }$, and $M^{\ast }$ can be
explicitly written. For instance
\[
G(\varphi )(t)=\int _0^{t}e^{-(\mathcal{B}_n+N^{-1})(t-s)}\varphi (s)\,ds,
\]
which show the non-local in time nature. Some useful properties of these
operators will be shown later (see Section 4).

Although the detailed statements of our results will be presented later, we
summarize now that we will prove the convergence of the optimal controls 
$v_{\varepsilon }\chi _{S_{\varepsilon }^{2,T}}\to v_0\chi
_{\omega ^{T}}$ strongly in $L^2(\omega ^{T})$, the convergence of the
corresponding states (extended to $\Omega $) 
$\tilde{u}_{\varepsilon }\rightharpoonup u_0$ weakly in
$L^2(0,T;H_0^{1}(\Omega ,\partial \Omega ))$ and
$\partial _{t}\tilde{u}_{\varepsilon }\rightharpoonup
\partial _{t}u_0$ weakly in $L^2(Q^{T})$, and in some sense, that will
be indicated later, the microscopic optimal control $v_{\varepsilon }$
converges to the macroscopic optimal control $v_0\in
H^{1}(0,T;L^2(\omega ))$, with the limit state problem given by
\begin{equation}
\begin{gathered}
\begin{aligned}
&\partial _{t}u_0(v)-\Delta u_0(v)+\mathcal{A}_n(u_0(v)-\mathcal{B}%
_nH(u_0(v)))\chi _{\omega ^{T}}  \\
&+\mathcal{A}_n(u_0(v)-\mathcal{B}_nM(u_0(v)))\chi _{(\Omega
\setminus \overline{\omega })\times (0,T)} \\
&=f+\mathcal{A}_n\mathcal{B}_nv\chi _{\omega ^{T}}, \quad (x,t)\in Q^{T},
\end{aligned}\\
u_0(v)(x,0)=0, \quad x\in \Omega , \\
u_0(v)(x,t)=0, \quad (x,t)\in \Gamma ^{T},
\end{gathered}
\end{equation}
where $\mathcal{A}_n = (n-2)C^{n-2}_0 \omega_n$, 
$v\in H^{1}(0,T;L^2(\omega ))$ with $v(x,0)=0$ and the limit cost
functional 
\begin{equation}
\begin{aligned}
&J_0(v) \\
&=  \frac{1}{2}\| \nabla u_0(v)\| _{L^2(Q^{T})}^2+\frac{1}{2}\|
u_0(v)(x,T)\| _{L^2(\Omega )}^2 + \frac{\mathcal{A}_n\mathcal{B}_n}{2N}
\int _{\omega ^{T}}(\partial _{t}G^{\ast }(H(u_0(v))))^2 \,dx\,dt
  \\
&\quad +\frac{\mathcal{A}_n\mathcal{B}_n}{2}\int _{\Omega \setminus
\overline{\omega }}|M(u_0(v))(x,T)|^2{dx}+\frac{\mathcal{A}_n}{2}
\int _{(\Omega \setminus \overline{\omega })\times (0,T)}|u_0(v)-\mathcal{B}
_nM(u_0(v))|^2\,dx\,dt    \\
&\quad +\frac{\mathcal{A}_n\mathcal{B}_n}{2}\int _{\omega }|H(u_0(v))(x,T)|^2{%
dx}+\frac{\mathcal{A}_n}{2}\int _{\omega ^{T}}|u_0(v)-\mathcal{B}%
_nH(u_0(v))|^2\,dx\,dt    \\
&\quad +\frac{N\mathcal{A}_n\mathcal{B}_n}{2}\int _{\omega ^{T}}(\partial
_{t}v)^2\,dx\,dt+\frac{N\mathcal{A}_n\mathcal{B}_n(\mathcal{B}_n+N^{-1})}{2}\int _{\omega }v^2(x,T){dx}+    \\
&\quad +\frac{N\mathcal{A}_n\mathcal{B}_n^2(\mathcal{B}_n+N^{-1})}{2}\int _{%
\omega ^{T}}v^2\,dx\,dt.
\end{aligned} \label{limit functional}
\end{equation}
of the optimal control problem
\begin{equation}
J_0(v_0)=\min _{v\in U_{\rm ad}}J_0(v),
\end{equation}
where the set of admissible functions is now
\[
U_{\rm ad}=\{ \psi \in H^{1}(0,T;L^2(\omega ))|\psi (x,0)=0\}.
\]

It can be seen that the first two terms and the last term of $J_0$ clearly
correspond to the three terms present in $J_{\varepsilon }$, but the rest of
the terms of $J_0$ are, in some way, unexpected. The terms of $J_0$
which are related to the final evaluation at time $T$ \ are new, and two of
them are actually non-local in time since they involve the operators $M$ and
$H$, respectively. The terms of $J_0$ which contain the operator $M$ are
integrals extended on the complementary of $\omega $, and they are a
consequence of the microscopic control $v_{\varepsilon }$ being applied only
at the boundary of some particles, $S_{\varepsilon }^2$, and not at all of
them. The unexpected terms of $J_0$ appear as a consequence of several
implicit relations that are justified in the proof of the Theorem 
\ref{thm: cost func limit} below. The last set of the terms that affect time
derivatives of a function of $u_0$ and the control $v$ are very surprising
since nothing suggests their appearance when observing the expression for 
$J_{\varepsilon }$.

To prove of these convergence results, we will use the
extension of the Pontryagin's method to the case of boundary controls (see,
e.g., \cite{LionsOpt}). In Section 2, we give the details of the formulation
of the direct problem and of  the coupled system arising in terms of the
adjoint optimal state $p_{\varepsilon }$: we will show that the optimal
control is given by 
$v_{\varepsilon }=-N^{-1}p_{\varepsilon }\chi_{S_{\varepsilon }^{2,T}}$. 
The a priori estimates allow passing to the
limit in the couple $(u_{\varepsilon },p_{\varepsilon })$ (and thus in the
controls $v_{\varepsilon }$) are obtained in Section 3. Some detailed
statements of the main theorems of this paper are collected in Section 5,
but before that, we present in Section 4 some properties of the auxiliary
non-local in time operators $G$, $H$ and $M$ defined above. The proof
characterizing the limit couple $(u_0,p_0)$ from the microscopic couple 
$(u_{\varepsilon },p_{\varepsilon })$ is given in Section 6. Finally, the
identification of the limit cost functional $J_0(v)$ from the microscopic
cost functional $J_{\varepsilon }(v)$ is obtained in Section 7.

\section{Problem statement and the adjoint problem}

Let $\Omega $ be a bounded domain in $\mathbb{R}^n$, $n\geq 3$, with a
smooth boundary $\partial \Omega $. We denote the unit cube $(-1/2,1/2)^n$
centered at the coordinates origin as $Y$. Let $G_0$ be a ball of radius 
$C_0$ such that $\overline{G}_0\subset Y$. Next, given a set $B$ of 
$\mathbb{R}^n$, by $\delta B$, $\delta >0$, we denote the set $\{x\in
\mathbb{R}^n|\delta ^{-1}x\in B\}$. For $\varepsilon >0$, we define 
$\widetilde{\Omega }_{\varepsilon }=\{x\in \Omega |\rho (x,\partial \Omega
)>2\varepsilon \}$, where $\rho $ is the Euclidean distance. 
Let $a_{\varepsilon }=C_0\varepsilon ^{\alpha }$, where $C_0$ is a positive
constant and $\alpha =\frac{n}{n-2}$. We define sets 
$G_{\varepsilon}^{j}=a_{\varepsilon }G_0+j$, where $j\in \mathbb{Z}^n$, $\mathbb{Z}^n$
is the set of vectors in $\mathbb{R}^n$ with integer coordinates. Now, we
introduce the set of indices $\Upsilon _{\varepsilon }=\{j\in \mathbb{Z}
^n:(a_{\varepsilon }G_0+{\varepsilon }j)\cap \widetilde{\Omega
_{\varepsilon }}\neq \emptyset \}$, note that the cardinal of 
$\Upsilon_{\varepsilon }$ satisfies that 
$|\Upsilon _{\varepsilon }|\cong d{\varepsilon }^{-n}$ for some $d=const>0$. 
Finally, we define the set
\[
G_{\varepsilon }=\cup_{j\in \Upsilon _{\varepsilon }}G_{\varepsilon }^{j}.
\]
Now, if we define $Y_{\varepsilon }^{j}={\varepsilon }Y+{\varepsilon }j$,
$P_{\varepsilon }^{j}={\varepsilon }j$, where $Y=(-1/2,1/2)^n$, then it is
easy to see that $\overline{G_{\varepsilon }^{j}}\subset Y_{\varepsilon }^{j}$
and the center of the ball $G_{\varepsilon }^{j}=a_{\varepsilon }G_0+{\varepsilon }j$
coincides with the center of the cube $Y_{\varepsilon }^{j}$.

In the formulation of the optimal control problem, we will consider only
some controllable region $\omega$, $\overline{\omega}\subset \Omega$, in the
whole domain $\Omega$. Thus, we split indices of $\Upsilon_\varepsilon$ into
two subsets $\Upsilon^2_{\varepsilon}=\{j\in \Upsilon_{\varepsilon}:
\overline{Y^{j}_{\varepsilon}}\subset \omega \}$ and 
$\Upsilon^{1}_\varepsilon = \Upsilon_\varepsilon \setminus \Upsilon^2_\varepsilon$.
Based on these sets, we will use the following notations
\[
G_{\varepsilon}^{1}=\cup_{j\in \Upsilon^{1}_{\varepsilon}}{G^{j}_{\varepsilon}},\quad
G_{\varepsilon}^2=\cup_{j\in
\Upsilon^2_{\varepsilon}}{G^{j}_{\varepsilon}}, \quad 
S^{1}_\varepsilon = \partial G^{1}_\varepsilon,\quad
S^2_\varepsilon = \partial G^2_\varepsilon.
\]
Further, we introduce the sets
\[
\Omega_{\varepsilon}=\Omega \setminus \overline{G_{\varepsilon}},\quad
\partial \Omega_{\varepsilon}=\partial \Omega \cup S_{\varepsilon},\quad
S_{\varepsilon}=S_{\varepsilon}^{1}\cup S_{\varepsilon}^2,
\]
and, for $0<T<\infty$, we define
\begin{gather*}
Q^{T}_{\varepsilon}=\Omega_{\varepsilon}\times (0,T),\quad
 \omega^{T}=\omega \times (0,T), \\
S_{\varepsilon}^{T}=S_{\varepsilon}\times(0,T),\quad
S_{\varepsilon}^{i,T}=S_{\varepsilon}^{i}\times (0,T),\quad i=1,2.
\end{gather*}

Now, we are in a position to formulate optimal control problem. 
Let $v\in L^2(0,T;S_{\varepsilon }^2)$. By $u_{\varepsilon }(v)$, we denote an
element of $L^2(0,T;H^{1}(\Omega _{\varepsilon },\partial \Omega ))$ with
the time derivative satisfying $\partial _{t}u_{\varepsilon }(v)\in
L^2(0,T;L^2(\Omega _{\varepsilon }))\cap
L^2(0,T;L^2(S_{\varepsilon }))$ and $u_{\varepsilon }(x,0)=0$ for $x\in
\Omega _{\varepsilon }\cup S_{\varepsilon }$, that is a solution to the
parabolic problem with the internal dynamic boundary condition. By $%
H^{1}(\Omega _{\varepsilon },\partial \Omega )$, we denote the closure with
respect to the norm $H^{1}(\Omega _{\varepsilon })$ of the set of infinitely
differentiable in $\overline{\Omega }_{\varepsilon }$ functions vanishing
near the boundary $\partial \Omega $. As a solution of
\eqref{init state prob}, we will consider a function $u_{\varepsilon }(v)$ with the
above-mentioned properties that satisfies the  integral identity
\begin{equation}
\begin{aligned}
&\int _{Q_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon }\varphi {\,dx\,dt%
}+\int _{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla \varphi {%
\,dx\,dt}+\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{T}}\partial
_{t}u_{\varepsilon }\varphi \,ds\,dt \\
&=\int _{Q_{\varepsilon
}^{T}}f\varphi\,dx\,dt+\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{2,T}}v\varphi \,ds\,dt
\end{aligned}  \label{int ident u}
\end{equation}
for an arbitrary function $\varphi \in L^2(0,T;H^{1}(\Omega _{\varepsilon
},\partial \Omega ))$. We consider now the optimal control problem stated in
the Introduction (see \eqref{optim control def}).

\begin{remark}
\label{Rem target nonzero} \rm
Our approach can be easily extended to the case of
a non-zero target $u_{T}\in L^2(\Omega )$, at least for a dense set of $%
u_{T}$ in $L^2(\Omega )$, i.e. the cost functional will be%
\begin{equation}
\begin{aligned}
J_{\varepsilon }(v) 
&= \frac{1}{2}\| \nabla {u_{\varepsilon }(v)}\|
_{L^2(Q_{\varepsilon }^{T})}^2+\frac{1}{2}\int _{\Omega
_{\varepsilon }}\left( u_{\varepsilon }(v)(x,T)-u_{T}\right) ^2dx
\\
&\quad +\frac{{\varepsilon }^{-\gamma }}{2}\int _{S_{\varepsilon
}}u_{\varepsilon }^2(v)(x,T){ds}+\varepsilon ^{-\gamma }\frac{N}{2}\| {v%
}\| _{L^2(S_{\varepsilon }^{2,T})}^2.
\end{aligned}  \label{nonzero cost}
\end{equation}
Indeed, let us assume that $u_{T}\in L^2(\Omega )$ is such that there
exists a converging as $\varepsilon \to0$ sequence of functions $%
V_{\varepsilon }\in L^2(0,T;L^2(\Omega _{\varepsilon }))$, i.e.
\begin{equation}
V_{\varepsilon }\rightharpoonup V_0\quad \text{weakly in }L^2(Q^{T}),
\text{ for some }V_0\in L^2(Q^{T}),  \label{hypo Weak  conver Ve}
\end{equation}
such that for the unique solution $U_{\varepsilon }$ of the auxiliary
problem
\begin{equation}
\begin{gathered}
\partial _{t}U_{\varepsilon }-\Delta U_{\varepsilon }=V_{\varepsilon }, \quad
(x,t)\in Q_{\varepsilon }^{T}, \\
\varepsilon ^{-\gamma }\partial _{t}U_{\varepsilon }+\partial _{\nu
}U_{\varepsilon }=0, \quad (x,t)\in S_{\varepsilon }^{T}, \\
U_{\varepsilon }(x,0)=0, \quad x\in \Omega _{\varepsilon }\cup S_{\varepsilon
}, \\
U_{\varepsilon }(x,t)=0, \quad (x,t)\in \Gamma ^{T},%
\end{gathered}
\end{equation}
we have
\begin{gather*}
\widetilde{U}_{\varepsilon }\rightharpoonup U_0\quad \text{weakly in }
L^2(0,T;H_0^{1}(\Omega )), \\
\partial _{t}\widetilde{U}_{\varepsilon }\rightharpoonup \partial _{t}U_0
\quad \text{weakly in }L^2(Q^{T}),
\end{gather*}
and
\[
U_0(x,T)=u_{T}(x)\quad \text{a.e. }\, x\in \Omega .
\]
Then by defining the change of variables
\[
w{_{\varepsilon }(v)=u_{\varepsilon }(v)-}U_{\varepsilon },
\]
where now ${u_{\varepsilon }(v)}$ is the optimal control associated to the
cost functional \eqref{nonzero cost}, we find that $w{_{\varepsilon }(v)}$
is the optimal control associate to the previous cost functional \eqref{cost
functional} with $u_{T}\equiv 0$. Finally, by the arguments of Remark 7.1
of \cite{DiPoShAttouch} (or Theorem 4 of \  \cite{DiPoShRACSAM}), it is easy
to prove that the set of final data $u_{T}\in L^2(\Omega )$ satisfying the
above mentioned conditions is a dense set of $L^2(\Omega )$. Then, the
perturbed equation satisfied by $w{_{\varepsilon }(v)}$, i.e.
\[
\partial _{t}w{_{\varepsilon }(v)}-\Delta w{_{\varepsilon }(v)}
=f-V_{\varepsilon },
\]
does not add any difficulty, once we know that \eqref{hypo Weak conver Ve}
holds.
\end{remark}

To obtain a characterization of the optimal control, we consider the
adjoint problem
\begin{equation}
\begin{gathered}
-\partial _{t}p_{\varepsilon }-\Delta p_{\varepsilon }=-\Delta
u_{\varepsilon }, \quad (x,t)\in Q_{\varepsilon }^{T}, \\
\partial _{\nu }p_{\varepsilon }-\varepsilon ^{-\gamma }\partial
_{t}p_{\varepsilon }=\partial _{\nu }u_{\varepsilon }, \quad (x,t)\in
S_{\varepsilon }^{T}, \\
p_{\varepsilon }(x,T)=u_{\varepsilon }(x,T), \quad x\in \Omega _{\varepsilon
}\cup S_{\varepsilon }, \\
p_{\varepsilon }(x,t)=0, \quad (x,t)\in \Gamma ^{T}.
\end{gathered} \label{init adjoint prob}
\end{equation}

We say that a function $p_{\varepsilon }\in L^2(0,T;H^{1}(\Omega
_{\varepsilon },\partial \Omega ))$, with 
$\partial _{t}p_{\varepsilon }\in
L^2(0,T;L^2(\Omega _{\varepsilon }))\cap L^2(0,T;L^2(S_{\varepsilon}))$, 
is a weak solution to  \eqref{init adjoint prob} if 
$p_{\varepsilon }(x,T)=u_{\varepsilon }(x,T)$ for a.e.\
 $x\in \Omega_{\varepsilon }$ and a.e.\ $x\in S_{\varepsilon }$ and if it satisfies the
integral identity
\begin{equation}
-\int _{Q_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon }\varphi \,dx\,dt
+\int _{Q_{\varepsilon }^{T}}\nabla p_{\varepsilon }\nabla
\varphi \,dx\,dt-\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{T}}\partial _{t}p_{\varepsilon }\varphi \,ds\,dt
=\int_{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla \varphi \,dx\,dt,
\label{int ident init adjoint prob}
\end{equation}
for any test function $\varphi \in L^2(0,T;H^{1}(\Omega _{\varepsilon
},\partial \Omega ))$. 
For a given $u_{\varepsilon }$ (with the regularity
of the weak solutions of \eqref{init state prob}) it is well-known that
there exists a unique solution to the problem \eqref{init adjoint prob}
(see, e.g., \cite{BejeDVrabie} and its references).

The following theorem gives a characterization of the optimal control 
$v_{\varepsilon }$ in terms of the adjoint state $p_{\varepsilon }$.

\begin{theorem}
\label{thm: init opt cont} 
Let the pair of functions $(u_\varepsilon(v_%
\varepsilon), v_\varepsilon)$ be an optimal solution of the problem 
\eqref{optim control def}, then $v_\varepsilon = -N^{-1}p_\varepsilon
\chi_{S^{2, T}_\varepsilon}$, 
where $p_\varepsilon$ is the solution to 
\eqref{init adjoint prob}. The converse is also true.
\end{theorem}

\begin{proof}
Let $v$ be an arbitrary function in $L^2(S^{2 ,T} _{\varepsilon })$ and 
$\lambda >0$. By $u_{\varepsilon }^{\lambda }$, we denote the solution of 
\eqref{optim control def} with the control $v_{\varepsilon }+\lambda v$,
i.e.\ $u_{\varepsilon }^{\lambda }=u_{\varepsilon }(v_{\varepsilon }+\lambda v)$. 
We use $u_{\varepsilon }=u_{\varepsilon }(v_{\varepsilon })$ to
simplify the notation. Then we have
\begin{align*}
&J_{\varepsilon }(v_{\varepsilon }+\lambda v)-J_{\varepsilon }(v_{\varepsilon })\\
&=\frac{1}{2}\| \nabla u_{\varepsilon }^{\lambda }\|
 _{L^2(Q_{\varepsilon }^{T})}^2+\frac{1}{2}\| u_{\varepsilon
 }^{\lambda }(x,T)\| _{L^2(\Omega _{\varepsilon })}^2 
 +\frac{\varepsilon ^{-\gamma }}{2}\| u_{\varepsilon }^{\lambda
 }(x,T)\| _{L^2(S_{\varepsilon })}^2 \\
&\quad +\varepsilon ^{-\gamma }\frac{N}{2}%
 \| v_{\varepsilon }+\lambda v\| _{L^2(S_{\varepsilon }^{2,T})}^2 
 -\frac{1}{2}\| \nabla u_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2
 -\frac{1}{2}\| u_{\varepsilon }(x,T)\| _{L^2(\Omega _{\varepsilon })}^2 \\
&\quad -\frac{\varepsilon ^{-\gamma }}{2}\| u_{\varepsilon }(x,T)\|
 _{L^2(S_{\varepsilon })}^2-\varepsilon ^{-\gamma }\frac{N}{2}\|
 v_{\varepsilon }\| _{L^2(S_{\varepsilon }^{2,T})}^2 \\
& =\frac{1}{2}\int _{Q_{\varepsilon }^{T}}\nabla (u_{\varepsilon
 }^{\lambda }-u_{\varepsilon })\nabla (u_{\varepsilon }^{\lambda
 }+u_{\varepsilon })\,dx\,dt+\frac{1}{2}\int _{\Omega _{\varepsilon
 }}(u_{\varepsilon }^{\lambda }-u_{\varepsilon })(x,T)(u_{\varepsilon
 }^{\lambda }+u_{\varepsilon })(x,T){dx} \\
&\quad +\frac{\varepsilon ^{-\gamma }}{2}\int _{S_{\varepsilon
 }}(u_{\varepsilon }^{\lambda }-u_{\varepsilon })(x,T)(u_{\varepsilon
 }^{\lambda }+u_{\varepsilon })(x,T){ds}+\varepsilon ^{-\gamma }\frac{N}{2}%
 \int _{S _{\varepsilon }^{2, T}}(2\lambda v_{\varepsilon }v+\lambda
^2v^2)\,ds\,dt.
\end{align*}
We define the function $\theta _{\varepsilon }=(u_{\varepsilon }^{\lambda
}-u_{\varepsilon })/\lambda $. 
It is easy to see that $\theta _{\varepsilon } $ is the unique solution to the problem
\begin{gather*}
\partial _{t}\theta _{\varepsilon }-\Delta \theta _{\varepsilon }=0, \quad
(x,t)\in Q_{\varepsilon }^{T}, \\
\varepsilon ^{-\gamma }\partial _{t}\theta _{\varepsilon }+\partial _{\nu
}\theta _{\varepsilon }=\chi _{S_{\varepsilon }^{2,T}}\varepsilon ^{-\gamma
}v, \quad (x,t)\in S_{\varepsilon }^{T}, \\
\theta _{\varepsilon }(x,0)=0, \quad  x\in \Omega _{\varepsilon }\cup
S_{\varepsilon }, \\
\theta _{\varepsilon }(x,t)=0, \quad (x,t)\in \Gamma ^{T}.%
\end{gather*}
Using the definition of $\theta _{\varepsilon }$, we have
\begin{align*}
J_{\varepsilon }'(v_{\varepsilon })v
&=\lim _{\lambda \to 0}(J_{\varepsilon }(v_{\varepsilon }+\lambda v)-J_{\varepsilon
 }(v_{\varepsilon }))/\lambda  \\
&=\int _{Q_{\varepsilon }^{T}}\nabla \theta _{\varepsilon }\nabla
u_{\varepsilon }\,dx\,dt+\int _{\Omega _{\varepsilon }}\theta
_{\varepsilon }(x,T)u_{\varepsilon }(x,T){dx} \\
&\quad +\varepsilon ^{-\gamma }\int _{S_{\varepsilon }}\theta _{\varepsilon
}(x,T)u_{\varepsilon }(x,T){ds}+\varepsilon ^{-\gamma }N\int
_{S_{\varepsilon }^{2,T}}v_{\varepsilon }v\,ds\,dt.
\end{align*}
Now, we use the definition of $p_{\varepsilon }$ and derive from the last
expression the identity
\[
J_{\varepsilon }'(v_{\varepsilon })v=\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}p_{\varepsilon }v\,dx\,dt+\varepsilon
^{-\gamma }N\int _{S_{\varepsilon }^{2,T}}v_{\varepsilon }v\,ds\,dt.
\]
As $v_{\varepsilon }$ is the optimal control, we should have
$J_{\varepsilon}'(v_{\varepsilon })\cdot v=0$ for all
$v\in L^2(0,T;L^2(S_{\varepsilon }^2))$. Hence,
$v_{\varepsilon }=-N^{-1}p_{\varepsilon }$ for a.e.\
 $(x,t)\in S_{\varepsilon }^{2,T}$. This
completes the proof.
\end{proof}

In consequence, by Theorem~\ref{thm: init opt cont}, the optimal control
problem is characterized through the coupled system
\begin{equation}
\begin{gathered}
\partial _{t}u_{\varepsilon }-\Delta u_{\varepsilon }=f, \quad (x,t)\in
Q_{\varepsilon }^{T}, \\
-\partial _{t}p_{\varepsilon }-\Delta p_{\varepsilon }=-\Delta
u_{\varepsilon }, \quad (x,t)\in Q_{\varepsilon }^{T}, \\
\partial _{\nu }u_{\varepsilon }+\varepsilon ^{-\gamma }\partial
_{t}u_{\varepsilon }=-\varepsilon ^{-\gamma }N^{-1}\chi _{S_{\varepsilon
}^{2, T}}p_{\varepsilon }, \quad (x,t)\in S_{\varepsilon }^{T}, \\
\partial _{\nu }p_{\varepsilon }-\varepsilon ^{-\gamma }\partial
_{t}p_{\varepsilon }=\partial _{\nu }u_{\varepsilon }, \quad (x,t)\in
S_{\varepsilon }^{T}, \\
u_{\varepsilon }(x,0)=0, \quad x\in \Omega _{\varepsilon }\cup S_{\varepsilon
}, \\
p_{\varepsilon }(x,T)=u_{\varepsilon }(x,T), \quad x\in \Omega _{\varepsilon
}\cup S_{\varepsilon }, \\
u_{\varepsilon }(x,t)=p_{\varepsilon }(x,t)=0, \quad (x,t)\in \Gamma ^{T}.%
\end{gathered}  \label{init coupled prob}
\end{equation}

\section{A priori estimates}

In this section, we obtain several a priori estimates of the state and adjoint
state. Taking $p_{\varepsilon }$ as a test function in the integral identity
for $u_{\varepsilon }$, we obtain
\begin{equation}
\begin{aligned}
&\int _{Q_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon
}p_{\varepsilon }\,dx\,dt+\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{T}}\partial _{t}u_{\varepsilon }p_{\varepsilon }\,ds\,dt+\int
_{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla p_{\varepsilon }{%
\,dx\,dt} \\
&=\int _{Q_{\varepsilon }^{T}}fp_{\varepsilon }\,dx\,dt%
-N^{-1}\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{2,T}}p_{\varepsilon }^2\,dx\,dt.
\end{aligned} \label{apr estim: int ident u p}
\end{equation}
Now, taking $u_{\varepsilon }$ as a test function in the integral identity
for $p_{\varepsilon }$, we obtain
\begin{equation}
-\int _{Q_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon
}u_{\varepsilon }\,dx\,dt-\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{T}}\partial _{t}p_{\varepsilon }u_{\varepsilon }\,ds\,dt+\int
_{Q_{\varepsilon }^{T}}\nabla p_{\varepsilon }\nabla u_{\varepsilon }{%
\,dx\,dt}
=\int _{Q_{\varepsilon }^{T}}|\nabla u_{\varepsilon }|^2\,dx\,dt.
\label{apr estim: int ident p u}
\end{equation}%
Next, we subtract \eqref{apr estim: int ident u p} from
\eqref{apr estim: int ident p u} and obtain the expression
\begin{align*}
&-\int _{Q_{\varepsilon }^{T}}\partial _{t}(u_{\varepsilon
}p_{\varepsilon })\,dx\,dt-\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{T}}\partial _{t}(u_{\varepsilon }p_{\varepsilon })\,ds\,dt%
-N^{-1}\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{2,T}}p_{\varepsilon }^2\,dx\,dt \\
&=\int _{Q_{\varepsilon }^{T}}|\nabla u_{\varepsilon }|^2\,dx\,dt-\int
_{Q_{\varepsilon }^{T}}fp_{\varepsilon }\,dx\,dt.
\end{align*}
From this, we obtain
\begin{equation}
\begin{aligned}
&\| \nabla u_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2+\|
u_{\varepsilon }(x,T)\| _{L^2(\Omega _{\varepsilon })}^2+\varepsilon
^{-\gamma }\| u_{\varepsilon }(x,T)\| _{L^2(S_{\varepsilon })}^2
+N^{-1}\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{2,T}}p_{\varepsilon }^2\,ds\,dt \\
&\leq \int _{Q_{\varepsilon
}^{T}}|f||p_{\varepsilon }|\,dx\,dt.
\end{aligned}  \label{apr estim: u estim 1}
\end{equation}
Then, we take $p_{\varepsilon }$ as a test function in the integral identity
\eqref{int ident init adjoint prob}, and obtain
\[
-\frac{1}{2}\| u_{\varepsilon }(x,T)\| _{L^2(\Omega _{\varepsilon
})}^2-\frac{\varepsilon ^{-\gamma }}{2}\| u_{\varepsilon }(x,T)\|
_{L^2(S_{\varepsilon })}^2+\| \nabla p_{\varepsilon }\|
_{L^2(Q_{\varepsilon }^{T})}^2\leq \int _{Q_{\varepsilon
}^{T}}\nabla u_{\varepsilon }\nabla p_{\varepsilon }\,dx\,dt.
\]
From here and \eqref{apr estim: u estim 1}, we conclude that
\begin{equation}
\begin{aligned}
\| \nabla p_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2
&\leq C(\| \nabla u_{\varepsilon }\| _{L^2(Q_{\varepsilon
}^{T})}^2+\| u_{\varepsilon }(x,T)\| _{L^2(\Omega _{\varepsilon
})}^2+\varepsilon ^{-\gamma }\| u_{\varepsilon }(x,T)\|
_{L^2(S_{\varepsilon })}^2)    \\
&\leq C\int _{Q_{\varepsilon }^{T}}|f||p_{\varepsilon }|\,dx\,dt.
\end{aligned} \label{prelim p estim}
\end{equation}
Here and below, constant $C$ is independent from $\varepsilon $.
As $p_{\varepsilon }$ is in $H^{1}(\Omega _{\varepsilon },\partial \Omega )$, we
can apply Poincar\'{e}-Friedrichs's inequality
\[
\| p_{\varepsilon }(\cdot,t)\| _{L^2(\Omega _{\varepsilon })}
\leq K\| \nabla p_{\varepsilon }(\cdot,t)\| _{L^2(\Omega _{\varepsilon })}.
\]
Using this inequality in the previous estimate \eqref{prelim p estim}, we obtain
\[
\| p_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2\leq C\|
f\| _{L^2(Q^{T})}^2.
\]
Now, we substitute this estimate into \eqref{apr estim: u estim 1}, and
derive the following estimate of $u_{\varepsilon }$,
\[
\| \nabla u_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2+\|
u_{\varepsilon }(x,T)\| _{L^2(\Omega _{\varepsilon })}^2+\varepsilon
^{-\gamma }\| u_{\varepsilon }(x,T)\| _{L^2(S_{\varepsilon
})}^2+N^{-1}\varepsilon ^{-\gamma }\| p_{\varepsilon }\|
_{L^2(S_{\varepsilon }^{2,T})}^2\leq C\| f\| _{L^2(Q^{T})}^2.
\]
From this, by \eqref{prelim p estim}, we obtain the estimation of the gradient
of $p_{\varepsilon }$,
\[
\| \nabla p_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2\leq
C\| f\| _{L^2(Q^{T})}^2.
\]%
Now we derive some estimates on the time derivatives of $u_{\varepsilon }$
and $p_{\varepsilon }$. We use Galerkin's approach and construct
$u_{\varepsilon }^m$ and $p_{\varepsilon }^m$, where $m=1,2,\dots $, that
are approximations to $u_{\varepsilon }$ and $p_{\varepsilon }$. Note that,
for such approximations, we have the same estimates derived above on $%
u_{\varepsilon }$ and $p_{\varepsilon }$. We take now $\partial
_{t}u_{\varepsilon }^m$ as a test function in the equations for $%
u_{\varepsilon }^m$, and integrating from $0$ to an arbitrary $\tau \in
\lbrack 0,T]$, we obtain
\begin{align*}
&\| \partial _{t}u_{\varepsilon }^m\| _{L^2(Q_{\varepsilon
}^{T})}^2+\varepsilon ^{-\gamma }\| \partial _{t}u_{\varepsilon
}^m\| _{L^2(S_{\varepsilon }^{T})}^2+\max _{t\in \lbrack
0,T]}\| \nabla u_{\varepsilon }^m\| _{L^2(\Omega _{\varepsilon
})}^2 \\
&\leq K(\int _{Q_{\varepsilon }^{T}}|f||\partial _{t}u_{\varepsilon
}^m|\,dx\,dt+\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{2,T}}|p_{\varepsilon }^m||\partial _{t}u_{\varepsilon }^m|\,ds\,dt) \\
&\leq \frac{1}{2}\| \partial _{t}u_{\varepsilon }^m\|
_{L^2(Q_{\varepsilon }^{T})}^2+\frac{\varepsilon ^{-\gamma }}{2}\|
\partial _{t}u_{\varepsilon }^m\| _{L^2(S_{\varepsilon
}^{T})}^2+K(\varepsilon ^{-\gamma }\| p_{\varepsilon }^m\|
_{L^2(S_{\varepsilon }^{2,T})}^2+\| f\| _{L^2(Q^{T})}^2),
\end{align*}
where constant $K$ is independent of $\varepsilon $ and $m$. From here, we
immediately derive
\[
\| \partial _{t}u_{\varepsilon }^m\| _{L^2(Q_{\varepsilon
}^{T})}^2+\varepsilon ^{-\gamma }\| \partial _{t}u_{\varepsilon
}^m\| _{L^2(S_{\varepsilon }^{T})}^2+\max _{t\in \lbrack
0,T]}\| \nabla u_{\varepsilon }^m\| _{L^2(\Omega _{\varepsilon
})}^2\leq K\| f\| _{L^2(Q^{T})}^2.
\]
Then, passing to the limit, as $m\to \infty $, in this estimate we
have
\[
\| \partial _{t}u_{\varepsilon }\| _{L^2(Q_{\varepsilon
}^{T})}^2+\varepsilon ^{-\gamma }\| \partial _{t}u_{\varepsilon }\|
_{L^2(S_{\varepsilon }^{T})}^2+\max _{t\in \lbrack 0,T]}\|
\nabla u_{\varepsilon }\| _{L^2(\Omega _{\varepsilon })}^2\leq K\|
f\| _{L^2(Q^{T})}^2.
\]
Moreover, if we use $\partial _{t}p_{\varepsilon }^m$ as a test function
in the equation for $p_{\varepsilon }^m$, we obtain, for a.e.\ $t$
\begin{align*}
&-\| \partial _{t}p_{\varepsilon }^m\| _{L^2(\Omega _{\varepsilon
 })}^2-\varepsilon ^{-\gamma }\| \partial _{t}p_{\varepsilon }^m\|
 _{L^2(S_{\varepsilon })}^2+(\nabla p_{\varepsilon }^m,\partial
 _{t}\nabla p_{\varepsilon }^m)_{L^2(\Omega _{\varepsilon })} \\
&=-(\partial _{t}u_{\varepsilon }^m,\partial _{t}p_{\varepsilon
}^m)_{L^2(\Omega _{\varepsilon })}-\varepsilon ^{-\gamma }(\partial
_{t}u_{\varepsilon }^m,\partial _{t}p_{\varepsilon
}^m)_{L^2(S_{\varepsilon })} \\
&\quad +(f,\partial _{t}p_{\varepsilon }^m)_{L^2(\Omega _{\varepsilon
})}-N^{-1}\varepsilon ^{-\gamma }(p_{\varepsilon }^m,\partial
_{t}p_{\varepsilon }^m)_{L^2(S_{\varepsilon }^2)}.
\end{align*}
Integrating this equality with respect to $t$ from $0$ to $T$, we obtain
\begin{align*}
&-\| \partial _{t}p_{\varepsilon }^m\| _{L^2(Q_{\varepsilon
 }^{T})}^2+\frac{1}{2}\| \nabla p_{\varepsilon }^m(\cdot ,T)\|
 _{L^2(\Omega _{\varepsilon })}^2-\frac{1}{2}\| \nabla p_{\varepsilon
 }^m(\cdot ,0)\| _{L^2(\Omega _{\varepsilon })}^2-\varepsilon
 ^{-\gamma }\| \partial _{t}p_{\varepsilon }^m\|
_{L^2(S_{\varepsilon }^{T})}^2 \\
&=-\int _{Q_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon
 }^m\partial _{t}p_{\varepsilon }^m\,dx\,dt-\varepsilon ^{-\gamma }\int
 _{S_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon }^m\partial
 _{t}p_{\varepsilon }^m\,ds\,dt \\
&\quad +\int _{Q_{\varepsilon }^{T}}f\partial _{t}p_{\varepsilon }^m\,dx\,dt
 -N^{-1}\varepsilon ^{-\gamma }\int _{S_{\varepsilon
 }^{2,T}}p_{\varepsilon }^m\partial _{t}p_{\varepsilon }^m\,ds\,dt \\
&=-\int _{Q_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon
 }^m\partial _{t}p_{\varepsilon }^m\,dx\,dt-\varepsilon ^{-\gamma }\int
 _{S_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon }^m\partial
 _{t}p_{\varepsilon }^m\,ds\,dt \\
&\quad +\int _{Q_{\varepsilon }^{T}}f\partial _{t}p_{\varepsilon }^m\,dx\,dt
 +N^{-1}\frac{\varepsilon ^{-\gamma }}{2}\| p_{\varepsilon }^m(x,0)\|
 _{L^2(S_{\varepsilon }^2)}^2-N^{-1}\frac{\varepsilon ^{-\gamma }}{2}%
 \| u_{\varepsilon }^m(x,T)\| _{L^2(S_{\varepsilon }^2)}^2.
\end{align*}
From this, we derive
\begin{align*}
&\| \partial _{t}p_{\varepsilon }^m\| _{L^2(Q_{\varepsilon
 }^{T})}^2+\varepsilon ^{-\gamma }\| \partial _{t}p_{\varepsilon
 }^m\| _{L^2(S_{\varepsilon }^{T})}^2+\frac{1}{2}\| \nabla
 p_{\varepsilon }^m(\cdot ,0)\| _{L^2(\Omega _{\varepsilon })}^2+%
 \frac{\varepsilon ^{-\gamma }}{2N}\| p_{\varepsilon }^m(x,0)\|
 _{L^2(\omega _{\varepsilon })}^2 \\
&\leq \frac{1}{2}\| \nabla u_{\varepsilon }^m(\cdot ,T)\|
 _{L^2(\Omega _{\varepsilon })}^2+\| \partial _{t}u_{\varepsilon
 }^m\| _{L^2(Q_{\varepsilon }^{T})}\| \partial _{t}p_{\varepsilon
 }^m\| _{L^2(Q_{\varepsilon }^{T})} \\
&\quad +\varepsilon ^{-\gamma }\| \partial _{t}u_{\varepsilon }^m\|
 _{L^2(S_{\varepsilon }^{T})}\| \partial _{t}p_{\varepsilon }^m\|
 _{L^2(S_{\varepsilon }^{T})}+\| f\| _{L^2(Q_{\varepsilon
 }^{T})}\| \partial _{t}p_{\varepsilon }^m\| _{L^2(Q_{\varepsilon
 }^{T})} + N^{-1}\frac{\varepsilon^{-\gamma}}{2}\| u^m_\varepsilon(x, T)\|^2_{L^2(S^2_\varepsilon)}.
\end{align*}
Finally, using the estimates obtained for $u_{\varepsilon }^m$, we conclude that
\begin{equation}
\| \partial _{t}p_{\varepsilon }^m\| _{L^2(Q_{\varepsilon
}^{T})}^2+\varepsilon ^{-\gamma }\| \partial _{t}p_{\varepsilon
}^m\| _{L^2(S_{\varepsilon }^{T})}^2\leq K\| f\|
_{L^2(Q^{T})}^2.
\end{equation}
Passing to the limit, as $m\to \infty $, we obtain the estimation of
$\partial _{t}p_{\varepsilon }$
\begin{equation}
\| \partial _{t}p_{\varepsilon }\| _{L^2(Q_{\varepsilon
}^{T})}^2+\varepsilon ^{-\gamma }\| \partial _{t}p_{\varepsilon }\|
_{L^2(S_{\varepsilon }^{T})}^2\leq K\| f\| _{L^2(Q^{T})}^2.
\end{equation}
Having proved some a priori estimates of $u_{\varepsilon }$ and
$v_{\varepsilon }$, we proceed with the extension of these solution to the
whole cylinder $Q^{T}$. We know (see, e.g. \cite{DiGoShBook} and its
references) that there exists an extension operator $P_{\varepsilon
}:H^{1}(Q_{\varepsilon }^{T})\to H^{1}(Q^{T})$ such that
\[
\| P_{\varepsilon }(u)\| _{H^{1}(Q^{T})}\leq \| u\|
_{H^{1}(Q_{\varepsilon }^{T})}.
\]
Let $\tilde{u}_{\varepsilon }$, $\tilde{p}_{\varepsilon }$ be the extensions
of the functions $u_{\varepsilon }$, $p_{\varepsilon }$. Then we obtain the
following estimates
\begin{gather}
\| \partial _{t}\tilde{u}_{\varepsilon }\| _{L^2(Q^{T})}^2+\|
\nabla \tilde{u}_{\varepsilon }\| _{L^2(Q^{T})}^2\leq K(\|
\partial _{t}u_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2+\|
\nabla u_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2),
\label{ext u estim} \\
\| \partial _{t}\tilde{p}_{\varepsilon }\| _{L^2(Q^{T})}^2+\|
\nabla \tilde{p}_{\varepsilon }\| _{L^2(Q^{T})}^2\leq K(\|
\partial _{t}p_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2+\|
\nabla p_{\varepsilon }\| _{L^2(Q_{\varepsilon }^{T})}^2).
\label{ext p estim}
\end{gather}
The obtained estimates imply that there exist some subsequences (still
denoted as the original) and some limit functions, $u_0$ and $p_0$, such
that
\begin{equation}
\begin{gathered}
\tilde{u}_{\varepsilon }\rightharpoonup u_0\quad
 \text{weakly in }L^2(0,T;H_0^{1}(\Omega )),\quad
\partial _{t}\tilde{u}_{\varepsilon }\rightharpoonup \partial _{t}u_0\quad
\text{weakly in }L^2(Q^{T}), \\
\tilde{p}_{\varepsilon }\rightharpoonup p_0\quad
 \text{weakly in }L^2(0,T;H_0^{1}(\Omega )),\quad
\partial _{t}\tilde{p}_{\varepsilon }\rightharpoonup \partial _{t}p_0\quad
 \text{weakly in }L^2(Q^{T}).
\end{gathered}  \label{limit func def}
\end{equation}
Moreover, the embedding theorem implies also that
$\tilde{u}_{\varepsilon}\to u_0$ and $\tilde{p}_{\varepsilon }\to p_0$ in
$L^2(Q^{T})$.

In the rest of the paper we  give the characterization of these limit
functions and derive a formulation of the homogenized optimal control
problem.

\section{Auxiliary non-local in time operators $G$, $H$ and $M$}\label{sec:aux}

As already  noted in the introduction, to state the homogenization
results we need to introduce some auxiliary problems. The first non-local
operator $M(\varphi )$ already was used in our previous study related to the
case of distributed controls (\cite{DiPoShAttouch}),
\begin{equation}
\begin{gathered}
\partial _{t}M(\varphi )+\mathcal{B}_nM(\varphi )=\varphi , \quad t\in (0,T),\\
M(\varphi )(0)=0, 
\end{gathered}  \label{tilde M prob}
\end{equation}
where $\mathcal{B}_n=(n-2)C_0^{-1}$, and $\varphi \in L^2(0,T)$ is
given. This operator is related to the region that is the complementary to
the one were the controls are localized. We also consider the adjoint
operator $M^{\ast }(\varphi )$ given by%
\begin{equation}
\begin{gathered}
-\partial _{t}M^{\ast }(\varphi )+\mathcal{B}_nM^{\ast }(\varphi )=\varphi, \quad
  t\in (0,T), \\
M^{\ast }(\varphi )(T)=0. 
\end{gathered}  \label{adj tilde M prob}
\end{equation}

For the domain that contains particles to which the controls are applied, we
introduce operator $G(\varphi )$ that satisfies the similar problem to 
\eqref{tilde M prob}, but with a different coefficient
\begin{equation}
\begin{gathered}
\partial _{t}G(\varphi )+(\mathcal{B}_n+N^{-1})G(\varphi )=\varphi , \quad
t\in (0,T), \\
G(\varphi )(0)=0.
\end{gathered}  \label{G prob 2}
\end{equation}
This operator $G(\varphi )$ can be explicitly written as
\[
G(\varphi )(t)=\int _0^{t}e^{-(\mathcal{B}_n+N^{-1})(t-s)}\varphi
(s)ds,
\]
which show the non-local in time nature. We define its adjoint operator
$G^{\ast }$ as the solution of the problem adjoint to \eqref{G prob 2}
\begin{equation}
\begin{gathered}
-\partial _{t}G^{\ast }(\varphi )+(\mathcal{B}_n+N^{-1})G^{\ast }(\varphi
)=\varphi , \quad t\in (0,T), \\
G^{\ast }(\varphi )(x,T)=0.
\end{gathered} \label{adj G prob 2}
\end{equation}
Nevertheless, it turns out that for the case of boundary controls, as we are
assuming in problem \eqref{init state prob}, we will need to define some new
operators $H$, and $H^{\ast }$, coupled with $G^{\ast }$ and $G$,
respectively, in the following way
\begin{equation}
\begin{gathered}
-\partial _{t}H^{\ast }(\varphi )+(\mathcal{B}_n+N^{-1})H^{\ast }(\varphi
)-N^{-1}(\mathcal{B}_n+N^{-1})G(H^{\ast }(\varphi ))=\varphi , \quad t\in
(0,T), \\
H^{\ast }(\varphi )(T)=0,
\end{gathered}  \label{adj H prob 2}
\end{equation}
and
\begin{equation}
\begin{gathered}
\partial _{t}H(\varphi )+(\mathcal{B}_n+N^{-1})H(\varphi )-N^{-1}(\mathcal{%
B}_n+N^{-1})G^{\ast }(H(\varphi ))=\varphi , \quad t\in (0,T), \\
H(\varphi )(0)=0.
\end{gathered}  \label{H prob 2}
\end{equation}
Notice that both problems are now non-local in time, but since the operators
$G$ and $G^{\ast }$ are globally Lipschitz continuous on $L^2(0,T)$,
we obtain the existence and uniqueness of the associate solutions once
$\varphi \in L^2(0,T)$ is given.

It is straightforward to show that the operator $G$ is the adjoint operator
to $G^{\ast }$, i.e.
\begin{equation}
\int _0^{T}G(\varphi )\psi {dt}=\int _0^{T}\varphi G^{\ast
}(\psi ){dt},  \label{G adj rel}
\end{equation}%
for any arbitrary functions $\varphi ,\, \psi \in L^2(0,T)$. Indeed, we
have
\begin{align*}
\int _0^{T}G(\varphi )\psi {dt}
&=\int _0^{T}G(\varphi
 )(-\partial _{t}G^{\ast }(\psi )+(\mathcal{B}_n+N^{-1})G^{\ast }(\psi )){dt} \\
&=\int _0^{T}(\mathcal{B}_n+N^{-1})G(\varphi )G^{\ast }(\psi ){dt}%
-G(\varphi )G^{\ast }(\varphi )\Big \vert_0^{T}+\int
_0^{T}\partial _{t}G(\varphi )G^{\ast }(\psi ){dt} \\
&=\int _0^{T}(\partial _{t}G(\varphi )+(\mathcal{B}
_n+N^{-1})G(\varphi ))G^{\ast }(\psi ){dt}=\int _0^{T}\varphi
G^{\ast }(\psi ){dt}.
\end{align*}
Similarly to the above argument we know that $M$ is the adjoint operator to 
$M^{\ast }$,
\[
\int _0^{T}M(\varphi )\psi {dt}=\int _0^{T}\varphi M^{\ast
}(\psi ){dt}.
\]

The case of operators $H$ and $H^{\ast }$ is less trivial. 
Nevertheless, we also have that, for any arbitrary functions 
$\varphi ,\, \psi \in L^2(0,T)$,
\begin{equation}
\int _0^{T}H(\varphi )\psi {dt} 
=\int _0^{T}\varphi H^{\ast
}(\psi ){dt}.  \label{H adj rel}
\end{equation}
Indeed, we make the following transformations
\begin{align*}
\int _0^{T}H(\varphi )\psi {dt}
&=\int _0^{T}H(\varphi )(-\partial _{t}H^{\ast }(\psi )
 +(\mathcal{B}_n+N^{-1})H^{\ast }(\psi)-N^{-1}
  (\mathcal{B}_n+N^{-1})G(H^{\ast }(\psi ))){dt} \\
&=\int _0^{T}\partial _{t}H(\varphi )H^{\ast }(\psi ){dt}
 +(\mathcal{B}_n+N^{-1})\int _0^{T}H(\varphi )H^{\ast }(\psi ){dt} \\
&\quad -N^{-1}(\mathcal{B}_n+N^{-1})\int _0^{T}G^{\ast }(H(\varphi
 ))H^{\ast }(\psi ){dt} \\
&=\int _0^{T}(\partial _{t}H(\varphi )+(\mathcal{B} _n+N^{-1})H(\varphi )
 -N^{-1}(\mathcal{B}_n+N^{-1})G^{\ast }(H(\varphi)))H^{\ast }(\psi ){dt} \\
&=\int _0^{T}\varphi H^{\ast }(\psi ){dt.}
\end{align*}
In addition to the above properties it will be useful to get some other
relations among the above operators.

\begin{lemma}
\label{lem:aux relations} 
For the functions $G$, $G^{\ast }$, $H$ and $H^{\ast }$ introduced in 
\eqref{G prob 2}-\eqref{adj H prob 2}, we have the
following relations
\begin{itemize}
\item[(i)] $G^{\ast }(H(\varphi ))=H^{\ast }(G(\varphi ))$, and 
$G(H^{\ast }(\varphi ))=H(G^{\ast }(\varphi ))$,

\item[(ii)] $H(\varphi )=G(\varphi )+N^{-1}(\mathcal{B}_n+N^{-1})G(H^{\ast
}(G(\varphi )))$, and we also have the adjoint version $H^{\ast }(\varphi
)=G^{\ast }(\varphi )+N^{-1}(\mathcal{B}_n+N^{-1})G^{\ast }(H(G^{\ast
}(\varphi )))$.
\end{itemize}
\end{lemma}

\begin{proof}
We start with the proof of (i). The relations given in (ii)
are direct consequences of the ones (i)). 
We consider two coupled
auxiliary linear systems. The first one is a system coupling the functions 
$A(t)=G^{\ast }(H(\varphi ))$ and $B(t)=H(\varphi )$. From
\eqref{adj G prob 2} and \eqref{H prob 2}, we obtain
\begin{equation}
\begin{gathered}
-\partial _{t}A+(\mathcal{B}_n+N^{-1})A=B, \quad t\in (0,T), \\
\partial _{t}B+(\mathcal{B}_n+N^{-1})B-N^{-1}(\mathcal{B}
_n+N^{-1})A=\varphi , \quad t\in (0,T), \\
A(T)=0,\quad B(0)=0. 
\end{gathered}   \label{lem:func sys 1}
\end{equation}
Analogously, from \eqref{G prob 2} and \eqref{adj H prob 2}, we obtain a second
system coupling the functions, $\tilde{A}(t)=G(H^{\ast }(G(\varphi )))$ and 
$\tilde{B}(t)=H^{\ast }(G(\varphi ))$:
\begin{equation}
\begin{gathered}
\partial _{t}\tilde{A}+(\mathcal{B}_n+N^{-1})\tilde{A}=\tilde{B}, \quad t\in
(0,T), \\
-\partial _{t}\tilde{B}+(\mathcal{B}_n+N^{-1})\tilde{B}-N^{-1}(\mathcal{B}%
_n+N^{-1})\tilde{A}=G(\varphi ), \quad t\in (0,T), \\
\tilde{A}(0)=0,\quad \tilde{B}(T)=0. 
\end{gathered}   \label{lem:func sys 2}
\end{equation}

From \eqref{lem:func sys 1}, we can derive a linear second order
ODE problem on the function $A$. To do this, we substitute the expression
for $B$ from the first equation of the system \eqref{lem:func sys 1} into
the second one. Thus, we obtain
\[
-\partial _{tt}^2A+(\mathcal{B}_n+N^{-1})\partial _{t}A-(\mathcal{B}%
_n+N^{-1})\partial _{t}A+(\mathcal{B}_n+N^{-1})^2A-N^{-1}(\mathcal{B}%
_n+N^{-1})A=\varphi ,
\]%
and, simplifying it, we obtain
\begin{equation}
-\partial _{tt}^2A+\mathcal{B}_n(\mathcal{B}_n+N^{-1})A=\varphi .
\label{second order ODE}
\end{equation}%
Also, substituting $B$ written in terms of $A$ into $B(0)=0$, we obtain a
condition on $A(0)$ (recall that we already have the condition at $t=T$ from
the definition of $A$). Hence, the two boundary conditions are
\begin{equation}
A(T)=0,\quad -\partial _{t}A(0)+(\mathcal{B}_n+N^{-1})A(0)=0.
\label{bound cond ODE}
\end{equation}
This is a linear coercive equation which has uniqueness of solutions. For
instance, if we consider the homogeneous case ($\varphi \equiv 0$) in the
equation \eqref{second order ODE} we obtain that, obviously, the trivial
solution satisfies this problem. Moreover, by multiplying the equation by $A$
and integrating from $0$ to $T$, we obtain
\[
\int _0^{T}|\partial _{t}A|^2{dt}+\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})\int _0^{T}A^2{dt}+(\mathcal{B}_n+N^{-1})A^2(0)=0
\]
and we immediately conclude that, necessarily, $A\equiv 0$. Thus, for any
$\varphi \in L^2(0,T)$, the inhomogeneous problem has also a unique
solution.

Let us obtain now the ODE satisfied by the function $\tilde{B}(t)$. From the
second equation of \eqref{lem:func sys 2}, we write $\tilde{A}$ in terms of 
$\tilde{B}$
\[
\tilde{A}=-N(\mathcal{B}_n+N^{-1})^{-1}\partial _{t}\tilde{B}+N\tilde{B}-N(%
\mathcal{B}_n+N^{-1})^{-1}G(\varphi ).
\]
Substituting this expression into the first equation of
\eqref{lem:func sys 2}, we derive
\begin{align*}
&-N(\mathcal{B}_n+N^{-1})^{-1}\partial _{tt}^2\tilde{B}+N\partial _{t}%
\tilde{B}-N(\mathcal{B}_n+N^{-1})^{-1}\partial _{t}G(\varphi ) \\
&-N\partial _{t}\tilde{B}+(\mathcal{B}_n+N^{-1})N\tilde{B}-NG(\varphi )=%
\tilde{B}.
\end{align*}
Combining similar terms and using the definition of $G(\varphi )$, we obtain
\[
-\partial _{tt}^2\tilde{B}+\mathcal{B}_n(\mathcal{B}_n+N^{-1})\tilde{B}%
=\partial _{t}G(\varphi )+(\mathcal{B}_n+N^{-1})G(\varphi )=\varphi .
\]
From condition $\tilde{A}(0)=0$, we conclude that
\[
-\partial _{t}\tilde{B}(0)+(\mathcal{B}_n+N^{-1})\tilde{B}(0)=0,\quad
\tilde{B}(T)=0.
\]
Therefore, we obtain exactly the same linear problem as for the function $A$.
But, since the solution to this problem is unique, we obtain that $A=\tilde{B}$,
or in other terms $G^{\ast }(H(\varphi ))=H^{\ast }(G(\varphi ))$. This
completes the proof of the first relation in \textit{(i)}.

Part (ii) can be proved using some similar arguments. Using the
definition of $H^{\ast }$, we have that
\[
G(\varphi )+N^{-1}(\mathcal{B}_n+N^{-1})G(H^{\ast }(G(\varphi
)))=-\partial _{t}H^{\ast }(G(\varphi ))+(\mathcal{B}_n+N^{-1})H^{\ast
}(G(\varphi )).
\]%
We use (i) and substitute $H^{\ast }(G(\varphi ))$ with
$G^{\ast}(H(\varphi ))$ in the right-hand side, and using the definition of
$G^{\ast}$, we obtain
\[
G(\varphi )+N^{-1}(\mathcal{B}_n+N^{-1})G(H^{\ast }(G(\varphi
)))=-\partial _{t}G^{\ast }(H(\varphi ))+(\mathcal{B}_n+N^{-1})G^{\ast
}(H(\varphi ))=H(\varphi ).
\]%
The second relation is proved in a similar way. This completes the proof.
\end{proof}

\section{Statement of the homogenization theorems}

Now, we are in a position to state the main theorem that characterizes the
pair of functions $(u_0, p_0)$ given by \eqref{limit func def}. 
The homogenized problem contains auxiliary functions defined in Section \ref{sec:aux}.

\begin{theorem} \label{homogenization theorem} 
Let $n\geq 3$, $a_{\varepsilon}=C_0\varepsilon ^{\gamma }$, where $C_0>0$, 
$\gamma = \frac{n}{n-2}$. If the
pair $(u_{\varepsilon },p_{\varepsilon })$ is the solution to the problem
\eqref{init coupled prob}, then $(u_0,p_0)$, defined in \eqref{limit func def}, 
is a solution to the system
\begin{equation}
\begin{gathered}
\begin{aligned}
&\partial _{t}u_0-\Delta u_0+\mathcal{A}_n(u_0-\mathcal{B}%
_nH(u_0))\chi _{\omega ^{T}}  
+\mathcal{A}_n(u_0-\mathcal{B}_nM(u_0))\chi _{(\Omega \setminus
\overline{\omega })\times (0,T)} \\
&=f-N^{-1}\mathcal{A}_n\mathcal{B}
_nH(G^{\ast }(p_0))\chi _{\omega ^{T}}, \quad (x,t)\in Q^{T},
\end{aligned} \\
\begin{aligned}
&-\partial _{t}p_0-\Delta p_0   
+\mathcal{A}_n(p_0-\mathcal{B}_nH^{\ast }(p_0))\chi _{\omega ^{T}}+%
\mathcal{A}_n(p_0-\mathcal{B}_nM^{\ast }(p_0))\chi _{(\Omega
\setminus \overline{\omega })\times (0,T)}   \\
&=-\Delta u_0+\mathcal{A}_n(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0))\chi _{\omega ^{T}}  \\
&\quad +\mathcal{A}_n(u_0-\mathcal{B}_n^2M^{\ast }(M(u_0)))\chi _{(\Omega
\setminus \overline{\omega })\times (0,T)}, \quad  (x,t)\in Q^{T} 
\end{aligned}\\
u_0(x,0)=0, \quad x\in \Omega , \\
u_0(x,t)=0, \quad  (x,t)\in \Gamma ^{T}, \\
p(x,T)=u_0(x,T), \quad  x\in \Omega , \\
p(x,t)=0, \quad  (x,t)\in \Gamma ^{T},
\end{gathered}
\end{equation}   \label{homogenized system}
where, $\mathcal{A}_n = (n-2)C^{n-2}_0\omega_n$, $\mathcal{B}_n = (n-2)C^{-1}_0$.
 Moreover, if we introduce the limit state problem
\begin{equation}
\begin{gathered}
\begin{aligned}
&\partial _{t}u_0(v)-\Delta u_0(v)+\mathcal{A}_n(u_0(v)-\mathcal{B}%
_nH(u_0(v)))\chi _{\omega ^{T}}   \\
&+\mathcal{A}_n(u_0(v)-\mathcal{B}_nM(u_0(v)))\chi _{(\Omega
\setminus \overline{\omega })\times (0,T)} \\
&=f+\mathcal{A}_n\mathcal{B}_nv\chi _{\omega ^{T}}, \quad (x,t)\in Q^{T},
\end{aligned} \\
u_0(v)(x,0)=0, \quad x\in \Omega , \\
u_0(v)(x,t)=0, \quad (x,t)\in \Gamma ^{T},%
\end{gathered}   \label{limit state problem func}
\end{equation}
where $v\in H^{1}(0,T;L^2(\omega ))$, with $v(x,0)=0$, and if we define
the cost functional $J_0(v)$ given by \eqref{limit functional}, then,
from \eqref{limit state problem func} and \eqref{limit func def}, we have
that $u_0(v)$ is the state associated to the optimal control problem
\begin{equation}
J_0(v_0)=\min _{v\in U_{\rm ad}}J_0(v),  \label{limit opt cont prob}
\end{equation}
where the set of admissible functions is
\[
U_{\rm ad}=\{ \psi \in H^{1}(0,T;L^2(\omega ))|\text{ }\psi (x,0)=0\}.
\]
\end{theorem}

In Section 7 we will prove the convergence of the sequence of functionals 
$J_{\varepsilon }$ to the limit functional $J_0$ given by \eqref{limit functional}. 
This will show that the system \eqref{homogenized system} characterizes the optimal 
control problem \eqref{limit opt cont prob}, i.e.\
that $v_0=-N^{-1}H(G^{\ast }(p_0))\chi _{\omega ^{T}}$. 
Then we have the following result.

\begin{theorem}
\label{thm: cost func limit} 
Under the conditions of Theorem \ref{homogenization theorem}, we have
\[
\lim _{\varepsilon \to 0}J_{\varepsilon }(v_{\varepsilon
})=J_0(v_0),
\]%
where $v_{\varepsilon }$ is the optimal control of the
\eqref{init state prob}, \eqref{cost functional}, and $v_0$ is the optimal control
of the limit optimal control problem \eqref{limit opt cont prob}.
\end{theorem}

 Lastly, we will show that  system \eqref{homogenized system} is
related to the limit functional and the limit optimal control $v_0$.

\begin{theorem}
\label{Thm convergence controls}
Let the pair of functions $(u_0(v_0),v_0)$ be an optimal solution of 
problem \eqref{limit opt cont prob}, then 
$v_0=-N^{-1}H(G^{\ast }(p_0))\chi _{\omega ^{T}}$,
where $p_0$ is the solution to \eqref{p_0 limit prob}.
\end{theorem}

\section{Proof of the homogenization theorem}

\begin{proof}
The main idea is to adapt to our setting the main lines of the so called
\textit{alternating test functions} (initially due to Luc Tartar to some
simple framework and then extended by many different authors: see, e.g., the
monograph \cite{DiGoShBook}). 
We introduce auxiliary functions $w_{\varepsilon }^{j}$, $j\in \mathbb{Z}^n$, 
that are solutions to the boundary-value problems
\begin{equation}
\begin{gathered}
\Delta {w}_{\varepsilon }^{j}=0, \quad x\in T_{{\varepsilon }/4}^{j}\setminus
\overline{G_{\varepsilon }^{j}}, \\
w_{\varepsilon }^{j}=1, \quad x\in \partial {G_{\varepsilon }^{j}}, \\
w_{\varepsilon }^{j}=0, \quad x\in \partial {T_{{\varepsilon }/4}^{j}},%
\end{gathered}   \label{def: auxiliary func w}
\end{equation}%
where $T_{{\varepsilon }/4}^{j}$ denotes the ball centered in 
$P_{\varepsilon }^{j}$ of ${{\varepsilon }/4}$ radii. Based on these
functions, we construct auxiliary functions in the whole domain $\Omega $
(where $i=1,2$)
\begin{equation}
W_{i,\varepsilon }=\begin{cases}
w_{\varepsilon }^{j}(x), & x\in T_{\varepsilon /4}^{j}\setminus \overline{%
G_{\varepsilon }^{j}},\; j\in \Upsilon _{\varepsilon }^{i}, \\
1, & x\in \overline{G_{\varepsilon }^{j}},\; j\in \Upsilon _{\varepsilon
}^{i}, \\
0, & x\in \Omega \setminus \overline{\cup _{j\in \Upsilon
_{\varepsilon }^{i}}{T_{{\varepsilon }/4}^{j}}}.%
\end{cases}   \label{def: W aux func}
\end{equation}%
They are related to controllable and uncontrollable sets of particles.
Note that $W_{i,\varepsilon }\in H_0^{1}(\Omega )$ and 
$W_{i,\varepsilon }\rightharpoonup 0$ weakly in $H_0^{1}(\Omega )$ as 
$\varepsilon \to 0$. By the embedding theorems, for some subsequence for
which we preserve the notation of the original, we have 
$W_{i,\varepsilon }\to 0$ strongly in $L^2(\Omega )$ as $\varepsilon \to 0$.

We will structure this long proof in a series of different steps.
\smallskip 

\noindent\textbf{Step A.1.}
 We take $W_{2,\varepsilon }H^{\ast }(\varphi )$, where 
$\varphi =\psi (x)\eta (t)$ with $\psi (x)\in C_0^{\infty }(\Omega )$, 
$\eta (t)\in C^{1}([0,T])$, as a test function in the integral identity 
\eqref{int ident u} and obtain
\begin{equation}
\begin{aligned}
&\int _{Q_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon
}W_{2,\varepsilon }H^{\ast }(\varphi )\,dx\,dt+\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}\partial _{t}u_{\varepsilon }H^{\ast
}(\varphi )\,ds\,dt+\int _{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon
}\nabla (W_{2,\varepsilon }H^{\ast }(\varphi ))\,dx\,dt    \\
&=\int _{Q_{\varepsilon }^{T}}fW_{2,\varepsilon }H^{\ast }(\varphi ){%
\,dx\,dt}-N^{-1}\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon }H^{\ast
}(\varphi )\,ds\,dt.
\end{aligned} \label{thm:W2 u first}
\end{equation}
Using the convergence \eqref{limit func def} and the properties of $W_{2,\varepsilon }$, 
we have
\begin{gather*}
\lim _{\varepsilon \to 0}\int _{Q_{\varepsilon
}^{T}}\partial _{t}u_{\varepsilon }W_{2,\varepsilon }H^{\ast }(\varphi ){\,dx\,dt%
}=0,\quad \lim _{\varepsilon \to 0}\int
_{Q_{\varepsilon }^{T}}fW_{2,\varepsilon }H^{\ast }(\varphi )\,dx\,dt=0,
\\
\lim _{\varepsilon \to 0}\int _{Q_{\varepsilon
}^{T}}\nabla u_{\varepsilon }\nabla (W_{2,\varepsilon }H^{\ast }(\varphi )){%
\,dx\,dt}=\lim _{\varepsilon \to 0}\int _{Q_{\varepsilon
}^{T}}\nabla (u_{\varepsilon }H^{\ast }(\varphi ))\nabla W_{2,\varepsilon }{%
\,dx\,dt}.
\end{gather*}
Thus, from \eqref{thm:W2 u first}, we derive
\[
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}u_{\varepsilon }H^{\ast }(\varphi )\,ds\,dt=-\int _{Q_{\varepsilon
}^{T}}\nabla W_{2,\varepsilon }\nabla (u_{\varepsilon }H^{\ast }(\varphi )){%
\,dx\,dt}-N^{-1}\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon }H^{\ast
}(\varphi )\,ds\,dt+\zeta _{\varepsilon },
\]
where, here and below, $\zeta _{\varepsilon }\to 0$ as $\varepsilon
\to 0$ (we will abuse the notation and will always use $\zeta
_{\varepsilon }$ for the terms converging to zero).

Next, we use the following fundamental relation 
(calling it ``from surface to volume averaging
convergence principle'' \cite[Theorem 4.5]{DiGoShBook} 
also applied, under different
formulations, by many authors, see
\cite{OlSh95, ZuShDiffEq}).
 We have
\begin{equation}
\int _{\Omega _{\varepsilon }}\nabla W_{i,\varepsilon }\nabla \eta {dx%
}=-\mathcal{A}_n\int _{\Omega ^{i}}\eta {dx}+\mathcal{B}%
_n\varepsilon ^{-\gamma }\sum _{j\in \Upsilon _{\varepsilon
}^{i}}\int _{\partial G_{\varepsilon }^{j}}\eta {ds}+\zeta
_{\varepsilon },  \label{w limit relation}
\end{equation}
here $\Omega ^{1}=\Omega \setminus \omega $, $\Omega ^2=\omega $, and then
\begin{align*}
&\lim _{\varepsilon \to 0}\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}\partial _{t}u_{\varepsilon }H^{\ast
}(\varphi )\,ds\,dt \\
&=\mathcal{A}_n\int _{\omega ^{T}}H(u_0)\varphi {%
\,dx\,dt}-\lim _{\varepsilon \to 0}\varepsilon ^{-\gamma }%
\mathcal{B}_n\int _{S_{\varepsilon }^{2,T}}u_{\varepsilon }H^{\ast
}(\varphi )\,ds\,dt
-N^{-1}\lim _{\varepsilon \to 0}\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}p_{\varepsilon }H^{\ast }(\varphi )\,ds\,dt.
\end{align*}
Using that $u_{\varepsilon }(x,0)=0$ and $H^{\ast }(\varphi )(x,T)=0$, we
integrate by parts the integral in the right-hand side, and further
transform the previous equality to obtain
\begin{equation}
\begin{aligned}
&\lim _{\varepsilon \to 0}\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}u_{\varepsilon }(-\partial _{t}H^{\ast
}(\varphi )+\mathcal{B}_nH^{\ast }(\varphi ))\,ds\,dt    \\
&=\mathcal{A}_n\int _{\omega ^{T}}H(u_0)\varphi \,dx\,dt-\lim
_{\varepsilon \to 0}N^{-1}\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}p_{\varepsilon }H^{\ast }(\varphi )\,ds\,dt.
\end{aligned} \label{thm: u limit rel 1}
\end{equation}

\noindent\textbf{Step A.2.} 
We take $\varphi =W_{2,\varepsilon }G(\varphi )$, where 
$\varphi =\psi (x)\eta (t)$ with $\psi \in C_0^{\infty }(\Omega )$, 
$\eta\in C^{1}([0,T])$ as a test function in the integral identity
\eqref{int ident init adjoint prob}, and obtain
\begin{equation}
\begin{aligned}
&-\int _{Q_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon
}W_{2,\varepsilon }G(\varphi )\,dx\,dt-\int _{S_{\varepsilon
}^{2,T}}\partial _{t}p_{\varepsilon }G(\varphi )\,ds\,dt
+\int _{Q_{\varepsilon }^{T}}\nabla p_{\varepsilon }\nabla (W_{2,
\varepsilon }G(\varphi ))\,dx\,dt    \\
&=\int _{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla
(W_{2,\varepsilon }G(\varphi ))\,dx\,dt.
\end{aligned}  \label{thm: W_2 p first}
\end{equation}
Taking $W_{2,\varepsilon }G(\varphi )$ as a test function in
\eqref{int ident u}, we have (using the same arguments as a above)
\[
\int _{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla
(W_{2,\varepsilon }G(\varphi ))\,dx\,dt=-\int _{S_{\varepsilon
}^{2,T}}\partial _{t}u_{\varepsilon }G(\varphi )\,ds\,dt-N^{-1}\varepsilon
^{-\gamma }\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon }G(\varphi ){%
\,ds\,dt}+\zeta _{\varepsilon }.
\]%
Substituting this relation into \eqref{thm: W_2 p first}, we derive
\begin{equation}
\begin{aligned}
&\int _{Q_{\varepsilon }^{T}}\nabla W_{2,\varepsilon }\nabla
(p_{\varepsilon }G(\varphi ))\,dx\,dt-\varepsilon ^{-\gamma
}\int _{S_{\varepsilon }^{2,T}}\partial _{t}(p_{\varepsilon
}-u_{\varepsilon })G(\varphi )\,ds\,dt    \\
&=\int _{Q_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon
}W_{2,\varepsilon }G(\varphi )\,dx\,dt-N^{-1}\varepsilon ^{-\gamma
}\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon }G(\varphi )\,ds\,dt+\zeta
_{\varepsilon }.
\end{aligned}  \label{thm: p ident G}
\end{equation}
Using the properties of $W_{2,\varepsilon }$ and a priori estimates of
$u_{\varepsilon }$ and $p_{\varepsilon }$, we conclude that the first
integral in the right-hand side of the above equality converges to zero as
$\varepsilon \to 0$. Also, using the relation \eqref{w limit relation}, we derive
\begin{equation}
\begin{aligned}
&-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}(p_{\varepsilon }-u_{\varepsilon })G(\varphi )\,ds\,dt+\varepsilon
^{-\gamma }\mathcal{B}_n\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon
}G(\varphi )\,ds\,dt    \\
&=\mathcal{A}_n\int _{\omega ^{T}}p_0G(\varphi )\,dx\,dt%
-N^{-1}\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{2,T}}p_{\varepsilon }G(\varphi )\,ds\,dt+\zeta _{\varepsilon }.
\end{aligned} \label{thm: p ident G 2}
\end{equation}
Note, that $p_{\varepsilon }(x,T)=u_{\varepsilon }(x,T)$ and
$G(\varphi)(0)=0$, therefore,
\[
-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}(p_{\varepsilon }-u_{\varepsilon })G(\varphi )\,ds\,dt=\varepsilon
^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial _{t}G(\varphi
)(p_{\varepsilon }-u_{\varepsilon })\,ds\,dt.
\]%
Thus, from \eqref{thm: p ident G 2}, we obtain
\[
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon
}(\partial _{t}G(\varphi )+(\mathcal{B}_n+N^{-1})G(\varphi ))\,ds\,dt \\
=\mathcal{A}_n\int _{\omega ^{T}}p_0G(\varphi )\,dx\,dt+\varepsilon
^{-\gamma }\int _{S_{\varepsilon }^{2,T}}u_{\varepsilon }\partial
_{t}G(\varphi )\,ds\,dt+\zeta _{\varepsilon }.
\]
Using the definition of $G(\varphi )$, we conclude that
\begin{equation}
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}p_{\varepsilon
}\varphi \,ds\,dt=\mathcal{A}_n\int _{\omega ^{T}}p_0G(\varphi ){\,dx\,dt%
}+\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}u_{\varepsilon
}(\varphi -(\mathcal{B}_n+N^{-1})G(\varphi ))\,ds\,dt+\zeta _{\varepsilon }.
\label{p limit with u}
\end{equation}
Then, we substitute \eqref{p limit with u} (with
$\varphi =H^{\ast }(\varphi)$) into \eqref{thm: u limit rel 1}, and obtain
\begin{align*}
&\lim _{\varepsilon \to 0}\varepsilon ^{-\gamma
}\int _{S_{\varepsilon }^{2,T}}u_{\varepsilon }(-\partial _{t}H^{\ast
}(\varphi )+(\mathcal{B}_n+N^{-1})H^{\ast }(\varphi )-N^{-1}(\mathcal{B}%
_n+N^{-1})G(H^{\ast }(\varphi )))\,ds\,dt \\
&=\mathcal{A}_n\int _{\omega ^{T}}H(u_0)\varphi \,dx\,dt-\mathcal{A}%
_nN^{-1}\int _{\omega ^{T}}p_0G(H^{\ast }(\varphi ))\,dx\,dt.
\end{align*}
Using the definition of $H^{\ast }$, we derive
\begin{equation}
\lim _{\varepsilon \to 0}\varepsilon ^{-\gamma
}\int _{S_{\varepsilon }^{2,T}}u_{\varepsilon }\varphi \,ds\,dt=\mathcal{%
A}_n\int _{\omega ^{T}}H(u_0)\varphi \,dx\,dt-\mathcal{A}%
_nN^{-1}\int _{\omega ^{T}}H(G^{\ast }(p_0))\varphi \,dx\,dt.
\label{thm: u limit rel}
\end{equation}

\noindent\textbf{Step A.3.}
Now, we can find the limit of the integrals taken over $S_{\varepsilon }^{2,T}$ 
in the integral identity~\eqref{int ident u}.
Indeed, we take $W_{2,\varepsilon }\varphi $ as a test function in 
\eqref{int ident u}, and obtain
\begin{equation}
\begin{aligned}
&\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}u_{\varepsilon }\varphi \,ds\,dt+\varepsilon ^{-\gamma }N^{-1}\int
_{S_{\varepsilon }^{2,T}}p_{\varepsilon }\varphi \,ds\,dt \\
&=-\int_{\Omega _{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla
(W_{2,\varepsilon }\varphi )\,dx\,dt+\int _{Q_{\varepsilon
}^{T}}fW_{2,\varepsilon }\varphi \,dx\,dt    \\
&=\mathcal{A}_n\int _{\omega ^{T}}u_0\varphi \,dx\,dt-\varepsilon
^{-\gamma }\mathcal{B}_n\int _{S_{\varepsilon
}^{2,T}}u_{\varepsilon }\varphi \,ds\,dt+\zeta _{\varepsilon }    \\
&=\mathcal{A}_n\int _{\omega ^{T}}u_0\varphi \,dx\,dt-\mathcal{A}_n%
\mathcal{B}_n\int _{\omega ^{T}}H(u_0)\varphi \,dx\,dt+\mathcal{A}%
_n\mathcal{B}_nN^{-1}\int _{\omega ^{T}}H(G^{\ast
}(p_0))\varphi \,dx\,dt+\zeta _{\varepsilon }.
\end{aligned}  \label{pa_t u conv control}
\end{equation}
Thus, we found the limit of the terms related to $S_{\varepsilon }^{2,T}$.
\smallskip

\noindent\textbf{Step A.4.} 
To complete the derivation of the limit equation, we
proceed with the terms related to $S_{\varepsilon }^{1,T}$. We take $%
W_{1,\varepsilon }M^{\ast }(\varphi )$ as a test function in the integral
identity \eqref{int ident u}, and obtain (again using the properties of the
function $W_{1,\varepsilon }$, we conclude that the integral with $f$
converges to zero)
\[
\int _{Q_{\varepsilon }^{T}}\nabla W_{1,\varepsilon }\nabla
(u_{\varepsilon }M^{\ast }(\varphi ))\,dx\,dt+\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{1,T}}\partial _{t}u_{\varepsilon }M^{\ast
}(\varphi )\,ds\,dt=\zeta _{\varepsilon },
\]
where $\zeta _{\varepsilon }\to 0$ as $\varepsilon \to 0$.
Using \eqref{w limit relation}, we derive
\[
\varepsilon ^{-\gamma }\mathcal{B}_n\int _{S_{\varepsilon
}^{1,T}}u_{\varepsilon }M^{\ast }(\varphi )\,ds\,dt+\varepsilon ^{-\gamma
}\int _{S_{\varepsilon }^{1,T}}\partial _{t}u_{\varepsilon }M^{\ast
}(\varphi )\,ds\,dt=\mathcal{A}_n\int _{(\Omega \setminus \overline{%
\omega })\times (0,T)}u_0M^{\ast }(\varphi )\,dx\,dt+\zeta _{\varepsilon }.
\]%
As $M^{\ast }(\varphi )(x,T)=0$ and $u_{\varepsilon }(x,0)=0$, we have
\[
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}\partial
_{t}u_{\varepsilon }M^{\ast }(\varphi )\,ds\,dt=-\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{1,T}}u_{\varepsilon }\partial _{t}M^{\ast
}(\varphi )\,ds\,dt.
\]
Thus,
\[
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}u_{\varepsilon
}(-\partial _{t}M^{\ast }(\varphi )+\mathcal{B}_nM^{\ast }(\varphi ))\,ds\,dt%
=\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}u_0M^{\ast }(\varphi )\,dx\,dt+\zeta _{\varepsilon }.
\]%
Using the definition of $M^{\ast }(\varphi )$, we obtain
\begin{equation}  \label{u conv uncont}
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}u_{\varepsilon
}\varphi \,ds\,dt=\mathcal{A}_n\int _{(\Omega \setminus \overline{%
\omega })\times (0,T)}M(u_0)\varphi \,dx\,dt+\zeta _{\varepsilon }.
\end{equation}
Now, we take $W_{1,\varepsilon }\varphi $ as a test function in the integral
identity \eqref{int ident u}, and using similar arguments as above, we
derive
\[
\int _{Q_{\varepsilon }^{T}}\nabla W_{1,\varepsilon }\nabla
(u_{\varepsilon }\varphi )\,dx\,dt+\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{1,T}}\partial _{t}u_{\varepsilon }\varphi \,ds\,dt%
=\zeta _{\varepsilon }.
\]
We transform this identity using relation \eqref{w limit relation} and
obtain
\begin{equation}
\begin{aligned}
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}\partial
_{t}u_{\varepsilon }\varphi \,ds\,dt
&=\mathcal{A}_n\int _{(\Omega
\setminus \overline{\omega })\times (0,T)}u_0\varphi \,dx\,dt-\mathcal{B}%
_n\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{1,T}}u_{\varepsilon }\varphi \,ds\,dt+\zeta _{\varepsilon }    \\
&=\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_nM(u_0))\varphi \,dx\,dt+\zeta _{\varepsilon }.
\end{aligned} \label{pa_t u conv uncontrol}
\end{equation}

Finally, using the convergence \eqref{pa_t u conv control} and
\eqref{pa_t u conv uncontrol}, we can pass to the limit in the integral 
identity for $u_{\varepsilon }$ and obtain the integral identity for $u_0$,
\begin{equation}
\begin{aligned}
&\int _{Q^{T}}\partial _{t}u_0\varphi \,dx\,dt+\int
_{Q^{T}}\nabla u_0\nabla \varphi \,dx\,dt    \\
&+\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_nM(u_0))\varphi \,dx\,dt+\mathcal{A}_n\int
_{w^{T}}(u_0-\mathcal{B}_nH(u_0))\varphi\,dx\,dt    \\
&=\int _{Q^{T}}f\varphi \,dx\,dt-N^{-1}\mathcal{A}_n\mathcal{B}%
_n\int _{w^{T}}H(G^{\ast }(p_0))\varphi \,dx\,dt.
\end{aligned} \label{u_0 int ident}
\end{equation}
Thus, $u_0$ is a solution to the  problem
\begin{gather*}
\begin{aligned}
&\partial _{t}u_0-\Delta u_0+\mathcal{A}_n(u_0-\mathcal{B}%
_nH(u_0))\chi _{\omega ^{T}}
+\mathcal{A}_n(u_0-\mathcal{B}_nM(u_0))\chi _{(\Omega \setminus
\overline{\omega })\times (0,T)} \\
&=f-N^{-1}\mathcal{A}_n\mathcal{B}_nH(G^{\ast }(p_0))
\chi _{\omega ^{T}}, \quad (x,t)\in Q^{T},
\end{aligned} \\
u_0(x,0)=0, \quad x\in \Omega , \\
u_0(x,t)=0, \quad (x,t)\in \Gamma ^{T}.
\end{gather*}

\noindent\textbf{Step B.1.}
Let us find the limit equation for $p_0$. We take 
$W_{2,\varepsilon }\varphi $ as a test function in the adjoint problem's
integral identity \eqref{int ident init adjoint prob}, and obtain
\[
-\int _{Q_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon
}W_{2,\varepsilon }\varphi \,dx\,dt-\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{2,T}}\partial _{t}p_{\varepsilon }\varphi \,ds\,dt%
+\int _{Q_{\varepsilon }^{T}}\nabla p_{\varepsilon }\nabla
(W_{2,\varepsilon }\varphi )\,dx\,dt=\int _{Q_{\varepsilon }^{T}}\nabla
u_{\varepsilon }\nabla (W_{2,\varepsilon }\varphi )\,dx\,dt.
\]%
The first integral in the left-hand side converges to zero because of the
properties of $W_{2,\varepsilon }$, thus,
\[
-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}p_{\varepsilon }\varphi \,ds\,dt=\int _{Q_{\varepsilon }^{T}}\nabla
(u_{\varepsilon }-p_{\varepsilon })\nabla (W_{2,\varepsilon }\varphi )\,dx\,dt%
+\zeta _{\varepsilon },
\]
Using \eqref{w limit relation}, \eqref{p limit with u} and
\eqref{thm: u limit rel}, we derive
\begin{align*}
&-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}p_{\varepsilon }\varphi \,ds\,dt\\
&=-\mathcal{A}_n\int _{w^{T}}(u_0-p_0)\varphi \,dx\,dt
 +\varepsilon ^{-\gamma }\mathcal{B}_n\int _{S_{\varepsilon }^{2,T}}
 (u_{\varepsilon }-p_{\varepsilon})\varphi \,ds\,dt
 + \zeta_\varepsilon \\
&=-\mathcal{A}_n\int _{w^{T}}(u_0-p_0)\varphi \,dx\,dt-\mathcal{A}%
_n\mathcal{B}_n\int _{w^{T}}p_0G(\varphi )\,dx\,dt+\varepsilon
^{-\gamma }\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int
_{S_{\varepsilon }^{2,T}}u_{\varepsilon }G(\varphi )\,ds\,dt +
\zeta_\varepsilon \\
&=\mathcal{A}_n\int _{w^{T}}p_0(\varphi -\mathcal{B}_nG(\varphi
))\,dx\,dt-\mathcal{A}_n\int _{w^{T}}u_0\varphi \,dx\,dt+\mathcal{A}%
_n\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int
_{Q^{T}}H(u_0)G(\varphi )\,dx\,dt \\
&\quad -N^{-1}\mathcal{A}_n\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int%
_{w^{T}}p_0G(H^{\ast }(G(\varphi )))\,dx\,dt + \zeta_\varepsilon.
\end{align*}
Combining the terms with $p_0$, we obtain the expression
\[
\varphi -\mathcal{B}_nG(\varphi )-N^{-1}\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G(H^{\ast }(G(\varphi )))\equiv E(\varphi ),
\]%
Using relation (ii) from Lemma~\ref{lem:aux relations}, we derive
\[
E(\varphi )=\varphi -\mathcal{B}_nH(\varphi ).
\]%
Therefore, the term with $p_0$ is
\[
\mathcal{A}_n\int _{\omega ^{T}}p_0(\varphi -\mathcal{B}%
_nH(\varphi ))\,dx\,dt=\mathcal{A}_n\int _{\omega ^{T}}(p_0-%
\mathcal{B}_nH^{\ast }(p_0))\varphi \,dx\,dt.
\]
Putting it all together, we derive
\begin{equation}
\begin{aligned}
&-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{2,T}}\partial
_{t}p_{\varepsilon }\varphi \,ds\,dt \\
&=\mathcal{A}_n\int _{w^{T}}(p_0-%
\mathcal{B}_nH^{\ast }(p_0))\varphi \,dx\,dt
-\mathcal{A}_n\int _{w^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})H^{\ast }(G(u_0)))\varphi \,dx\,dt.
\end{aligned}  \label{pa_t p cont}
\end{equation}

\noindent\textbf{Step B.2.}
Next, we deal with the parts related to $S_{\varepsilon}^{1,T}$. 
We take $W_{1,\varepsilon }M(\varphi )$ as a test function in the
integral identity \eqref{int ident init adjoint prob}, and obtain
\begin{align*}
&-\int _{Q_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon
}W_{1,\varepsilon }M(\varphi )\,dx\,dt-\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{1,T}}\partial _{t}p_{\varepsilon }M(\varphi ){\,ds\,dt%
}+\int _{Q_{\varepsilon }^{T}}\nabla W_{1,\varepsilon }\nabla
(p_{\varepsilon }M(\varphi ))\,dx\,dt \\
&=\int _{Q_{\varepsilon }^{T}}\nabla W_{1,\varepsilon }\nabla
(u_{\varepsilon }M(\varphi ))\,dx\,dt + \zeta _{\varepsilon }.
\end{align*}
Using integral identity \eqref{int ident u} taken with the same test
function, we transform the last relation to the form
\begin{align*}
&\int _{Q_{\varepsilon }^{T}}\nabla W_{1,\varepsilon }\nabla
(p_{\varepsilon }M(\varphi ))\,dx\,dt-\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{1,T}}\partial _{t}(p_{\varepsilon
}-u_{\varepsilon })M(\varphi )\,ds\,dt \\
&=\int _{Q_{\varepsilon }^{T}}\partial _{t}(p_{\varepsilon
}-u_{\varepsilon })W_{1,\varepsilon }M(\varphi )\,dx\,dt+\int
_{Q_{\varepsilon }^{T}}fW_{1,\varepsilon }M(\varphi )\,dx\,dt+\zeta
_{\varepsilon }.
\end{align*}
Properties of $W_{1,\varepsilon }$ imply that the terms at the right-hand
side converges to zero as $\varepsilon \to 0$. We use
\eqref{w limit relation} and transform the first term in the left-hand side 
and finally obtain
\[
-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}\partial
_{t}(p_{\varepsilon }-u_{\varepsilon })M(\varphi )\,ds\,dt+\varepsilon
^{-\gamma }\mathcal{B}_n\int _{S_{\varepsilon }^{1,T}}p_{\varepsilon }M(\varphi ){%
\,ds\,dt}=\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega }%
)\times (0,T)}p_0M(\varphi )\,dx\,dt+\zeta _{\varepsilon }.
\]
Integrating by parts in the first term, we conclude that
\[
\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}(p_{\varepsilon
}-u_{\varepsilon })\partial _{t}M(\varphi )\,ds\,dt+\varepsilon ^{-\gamma
}\mathcal{B}_n \int _{S_{\varepsilon }^{1,T}}p_{\varepsilon }M(\varphi )\,ds\,dt=%
\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}p_0M(\varphi )\,dx\,dt+\zeta _{\varepsilon }.
\]
Using the definition of $M^{\ast }$ and the convergence
\eqref{u conv uncont}, we derive
\begin{equation}
\begin{aligned}
&\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}p_{\varepsilon
}\varphi \,ds\,dt=\varepsilon ^{-\gamma }\int _{S_{\varepsilon
}^{1,T}}p_{\varepsilon }(\partial _{t}M(\varphi )+\mathcal{B}_nM(\varphi ))%
\,ds\,dt    \\
&=\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}p_0M(\varphi )\,dx\,dt+\mathcal{A}_n\int _{(\Omega \setminus
\overline{\omega })\times (0,T)}M(u_0)(\varphi -\mathcal{B}_nM(\varphi ))
\,dx\,dt+\zeta _{\varepsilon }.
\end{aligned}  \label{p W_1 relation}
\end{equation}
Lastly, we take $W_{1,\varepsilon }\varphi $ as a test function in the
integral identity \eqref{int ident init adjoint prob}, and using
\eqref{pa_t u conv uncontrol} and \eqref{p W_1 relation}, we obtain
\begin{align*}
&-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}\partial
_{t}p_{\varepsilon }\varphi \,ds\,dt \\
&=\int _{Q_{\varepsilon }^{T}}\nabla
(u_{\varepsilon }-p_{\varepsilon })\nabla (W_{1,\varepsilon }\varphi )\,dx\,dt%
+\zeta _{\varepsilon } \\
&=\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(p_0-u_0)\varphi \,dx\,dt+\varepsilon ^{-\gamma }\mathcal{B}_n\int
_{S_{\varepsilon }^{1,T}}(u_{\varepsilon }-p_{\varepsilon })\varphi {%
\,ds\,dt}+\zeta _{\varepsilon } \\
&=\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(p_0-u_0)\varphi \,dx\,dt+\mathcal{A}_n\mathcal{B}_n\int
_{(\Omega \setminus \overline{\omega })\times (0,T)}M(u_0)\varphi {%
\,dx\,dt} \\
&\quad -\mathcal{A}_n\mathcal{B}_n\int _{(\Omega \setminus \overline{%
\omega })\times (0,T)}p_0M(\varphi )\,dx\,dt-\mathcal{A}_n\mathcal{B}%
_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}M(u_0)(\varphi -\mathcal{B}_nM(\varphi ))\,dx\,dt+\zeta
_{\varepsilon }.
\end{align*}
Grouping similar terms, we derive
\begin{equation}
\begin{aligned}
-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }^{1,T}}\partial
_{t}p_{\varepsilon }\varphi \,ds\,dt
&=\mathcal{A}_n\int _{(\Omega
\setminus \overline{\omega })\times (0,T)}(p_0-\mathcal{B}_nM^{\ast
}(p_0))\varphi \,dx\,dt    \\
&\quad -\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_n^2M^{\ast }(M(u_0)))\varphi \,dx\,dt+\zeta
_{\varepsilon }.
\end{aligned}  \label{pa_t p uncont}
\end{equation}

\noindent\textbf{Step B.3.}
Now, using convergence \eqref{pa_t p cont} and 
\eqref{pa_t p uncont}, we can pass to the limit as $\varepsilon \to 0$ 
in \eqref{int ident init adjoint prob}, and obtain the integral identity for $p_0$,
\begin{equation}
\begin{aligned}
&-\int _{Q^{T}}\partial _{t}p_0\varphi \,dx\,dt+\int _{Q^{T}}%
\nabla p_0\nabla \varphi \,dx\,dt
+\mathcal{A}_n\int _{\omega ^{T}}(p_0-\mathcal{B}_nH^{\ast
}(p_0))\varphi \,dx\,dt \\
&+\mathcal{A}_n\int _{(\Omega \setminus
\overline{\omega })\times (0,T)}(p_0-\mathcal{B}_nM^{\ast
}(p_0))\varphi \,dx\,dt    \\
&=\int _{Q^{T}}\nabla u_0\nabla \varphi \,dx\,dt+\mathcal{A}%
_n\int _{Q^{T}}(u_0-\mathcal{B}_n(\mathcal{B}_n+N^{-1})G^{\ast
}(H(u_0))\varphi \,dx\,dt    \\
&\quad +\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_n^2M^{\ast }(M(u_0)))\varphi \,dx\,dt,
\end{aligned} \label{p_0 int ident}
\end{equation}
that is valid for an arbitrary function $\varphi \in L^2(0,T;H_0^{1}(\Omega ))$.
Thus, the limit adjoint problem is
\begin{equation}
\begin{gathered}
\begin{aligned}
&-\partial _{t}p_0-\Delta p_0
+\mathcal{A}_n(p_0-\mathcal{B}_nH^{\ast }(p_0))\chi _{\omega ^{T}}+%
\mathcal{A}_n(p_0-\mathcal{B}_nM^{\ast }(p_0))\chi _{(\Omega
\setminus \overline{\omega })\times (0,T)}   \\
&=-\Delta u_0+\mathcal{A}_n(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0))\chi _{\omega ^{T}}  \\
&\quad +\mathcal{A}_n(u_0-\mathcal{B}_n^2M^{\ast }(M(u_0)))\chi _{(\Omega
\setminus \overline{\omega })\times (0,T)}, \quad (x,t)\in Q^{T},
\end{aligned} \\
p(x,T)=u_0(x,T), \quad x\in \Omega , \\
p(x,t)=0, \quad (x,t)\in \Gamma ^{T}.%
\end{gathered}   \label{p_0 limit prob}
\end{equation}
Notice that, due to the well-posedness of the limit problem, all the
convergences hold for the whole sequences and not only for the considered
subsequences.This completes the proof of the homogenization theorem.
\end{proof}

\section{Proof of the limit cost functional and controls convergence theorems}

In this Section we will complete the characterization of the limit optimal
control.

\begin{proof}[Proof of Theorem \ref{thm: cost func limit}]
 We start with the expression of $J_{\varepsilon }$ taken at 
 $v_{\varepsilon }=-N^{-1}p_{\varepsilon }\chi_{S_{\varepsilon }^{2,T}}$:
\begin{align*}
&2J_{\varepsilon }(v_{\varepsilon })\\
&=\int _{Q_{\varepsilon}^{T}}|\nabla u_{\varepsilon }|^2\,dx\,dt
 +\int _{\Omega _{\varepsilon}}|u_{\varepsilon }(x,T)|^2{dx} 
 +\varepsilon ^{-\gamma }\int _{S_{\varepsilon }}|u_{\varepsilon
 }(x,T)|^2{ds}+N^{-1}\varepsilon ^{-\gamma }\int _{S_{\varepsilon
 }^{2,T}}p_{\varepsilon }^2\,ds\,dt.
\end{align*}
Using the integral identity \eqref{int ident u}, we obtain
\begin{equation}
\begin{aligned}
2J_{\varepsilon }(v_{\varepsilon })
&=\int _{Q_{\varepsilon
}^{T}}|\nabla u_{\varepsilon }|^2\,dx\,dt
 +\int _{\Omega _{\varepsilon}}|u_{\varepsilon }(x,T)|^2{dx}
 +\varepsilon ^{-\gamma }\int_{S_{\varepsilon }}|u_{\varepsilon }(x,T)|^2{ds}
 +\int_{Q_{\varepsilon }^{T}}fp_{\varepsilon }\,dx\,dt    \\
&\quad -\int _{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla
 p_{\varepsilon }\,dx\,dt
 -\int _{Q_{\varepsilon }^{T}}\partial
 _{t}u_{\varepsilon }p_{\varepsilon }\,dx\,dt
 -\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon }
p_{\varepsilon }\,ds\,dt.
\end{aligned}  \label{thm: cost func step 1}
\end{equation}
From the integral identity \eqref{int ident init adjoint prob}, we have
\begin{align*}
&-\int _{Q_{\varepsilon }^{T}}\nabla u_{\varepsilon }\nabla
p_{\varepsilon }\,dx\,dt-\int _{Q_{\varepsilon }^{T}}\partial
_{t}u_{\varepsilon }p_{\varepsilon }\,dx\,dt-\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{T}}\partial _{t}u_{\varepsilon }p_{\varepsilon }{%
\,ds\,dt} \\
&=-\int _{Q_{\varepsilon }^{T}}\nabla p_{\varepsilon }\nabla
u_{\varepsilon }\,dx\,dt+\int _{Q_{\varepsilon }^{T}}\partial
_{t}p_{\varepsilon }u_{\varepsilon }\,dx\,dt+\varepsilon ^{-\gamma }\int
_{S_{\varepsilon }^{T}}\partial _{t}p_{\varepsilon }u_{\varepsilon }{%
\,ds\,dt} \\
&\quad -\varepsilon ^{-\gamma }\int _{S_{\varepsilon }}|u_{\varepsilon
}(x,T)|^2{ds}-\int _{\Omega_{\varepsilon }}|u_{\varepsilon
}(x,T)|^2{dx} \\
&=-\int _{Q_{\varepsilon }^{T}}|\nabla u_{\varepsilon }|^2\,dx\,dt%
-\varepsilon ^{-\gamma }\int _{S_{\varepsilon }}|u_{\varepsilon
}(x,T)|^2{ds}-\int _{\Omega_{\varepsilon }}|u_{\varepsilon
}(x,T)|^2{dx}.
\end{align*}
Substituting this equality into \eqref{thm: cost func step 1}, we obtain
\[
2J_{\varepsilon }(v_{\varepsilon })=\int _{Q_{\varepsilon
}^{T}}fp_{\varepsilon }\,dx\,dt.
\]
Then, passing to the limit as $\varepsilon \to 0$ we have
\[
\lim _{\varepsilon \to 0}2J_{\varepsilon }(v_{\varepsilon})
=\int _{Q^{T}}fp_0\,dx\,dt.
\]%
Using the integral identity \eqref{u_0 int ident} for the function $u_0$,
we obtain
\begin{align*}
\int _{Q^{T}}fp_0\,dx\,dt
&=\int _{Q^{T}}\partial _{t}u_0p_0{%
\,dx\,dt}+\int _{Q^{T}}\nabla u_0\nabla p_0\,dx\,dt+\mathcal{A}_n\int
_{\omega ^{T}}(u_0-\mathcal{B}_nH(u_0))p_0\,dx\,dt \\
&\quad +\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_nM(u_0))p_0\,dx\,dt+N^{-1}\mathcal{A}_n%
\mathcal{B}_n\int _{\omega ^{T}}H(G^{\ast }(p_0))p_0\,dx\,dt \\
&=-\int _{Q^{T}}\partial _{t}p_0u_0\,dx\,dt+\int _{\Omega
}|u_0(x,T)|^2{dx}+\int _{Q^{T}}\nabla p_0\nabla u_0\,dx\,dt \\
&\quad +\mathcal{A}_n\int _{\omega ^{T}}(p_0-\mathcal{B}_nH^{\ast
}(p_0))u_0\,dx\,dt+\mathcal{A}_n\int _{(\Omega \setminus
\overline{\omega })\times (0,T)}(p_0-\mathcal{B}_nM^{\ast }(p_0))u_0{%
\,dx\,dt} \\
&\quad +N^{-1}\mathcal{A}_n\mathcal{B}_n\int _{\omega ^{T}}H(G^{\ast
}(p_0))p_0\,dx\,dt.
\end{align*}
Then, we use the integral identity \eqref{p_0 int ident} for the function
$p_0$ and duality relation for the operators $H$ and $M$, to conclude that
\begin{align*}
\int _{Q^{T}}fp_0\,dx\,dt
&=\int _{\Omega }|u_0(x,T)|^2\,dx\,dt%
+\int _{Q^{T}}|\nabla u_0|^2\,dx\,dt \\
&\quad +\mathcal{A}_n\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))u_0\,dx\,dt \\
&\quad +\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_n^2M^{\ast }(M(u_0)))u_0\,dx\,dt \\
&\quad +N^{-1}\mathcal{A}_n\mathcal{B}_n\int _{\omega ^{T}}H(G^{\ast
}(p_0))p_0\,dx\,dt.
\end{align*}
Let us show that the derived expression is non-negative. According to the
definition of $M$ and $M^{\ast }$, we have
\begin{align*}
&\int _{(\Omega \setminus \overline{\omega })\times (0,T)}(u_0-
\mathcal{B}_n^2M^{\ast }(M(u_0)))u_0\,dx\,dt \\
&=\int _{(\Omega \setminus \overline{\omega })\times (0,T)}(u_0^2-\mathcal{B}%
_n^2M^2(u_0))\,dx\,dt \\
&=\int _{(\Omega \setminus \overline{\omega })\times (0,T)}((\partial
_{t}M(u_0)+\mathcal{B}_nM(u_0))^2-\mathcal{B}_n^2M^2(u_0)){%
\,dx\,dt} \\
&=\int _{(\Omega \setminus \overline{\omega })\times (0,T)}(\partial
_{t}M(u_0))^2\,dx\,dt+\mathcal{B}_n\int _{\Omega \setminus
\overline{\omega }}|M(u_0)(x,T)|^2{dx} \\
&=\int _{(\Omega \setminus \overline{\omega })\times (0,T)}(u_0-%
\mathcal{B}_nM(u_0))^2\,dx\,dt+\mathcal{B}_n\int _{\Omega
\setminus \overline{\omega }}|M(u_0)(x,T)|^2{dx} \ge 0.
\end{align*}

Using that $G^{\ast }(H(u_0))$ is the solution to the problem
\eqref{second order ODE} with boundary conditions \eqref{bound cond ODE}, we obtain
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}
_n+N^{-1})G^{\ast }(H(u_0)))u_0\,dx\,dt \\
&=\int _{\omega ^{T}}(u_0-\partial _{tt}^2G^{\ast
}(H(u_0))-u_0)u_0\,dx\,dt \\
&=-\int _{\omega ^{T}}\partial
_{tt}^2G^{\ast }(H(u_0))u_0\,dx\,dt \\
&=\int _{\omega ^{T}}\partial _{tt}^2G^{\ast }(H(u_0))(\partial
_{tt}^2G^{\ast }(H(u_0))-\mathcal{B}_n(\mathcal{B}_n+N^{-1})G^{\ast
}(H(u_0)))\,dx\,dt \\
&=\int _{\omega ^{T}}(\partial _{tt}^2G^{\ast }(H(u_0)))^2\,dx\,dt-%
\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int _{\omega ^{T}}\partial
_{tt}^2G^{\ast }(H(u_0))G^{\ast }(H({u_0}))\,dx\,dt \\
&=\int _{\omega ^{T}}(\partial _{tt}^2G^{\ast }(H(u_0)))^2\,dx\,dt+%
\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int _{\omega ^{T}}(\partial
_{t}G^{\ast }(H(u_0)))^2\,dx\,dt \\
&\quad +\mathcal{B}_n\int _{\omega }(\partial _{t}G^{\ast
}(H(u_0))(x,0))^2{dx} \ge 0.
\end{align*}
From this, we y see that the expression is non-zero, however, we
convert it to a form that is  similar to the one derived for the terms
with $M$. We have
\begin{align*}
&\mathcal{B}_n\int _{\omega }(\partial _{t}G^{\ast
}(H(u_0))(x,0))^2{dx} \\
&=\mathcal{B}_n\int _{\omega }((\partial _{t}G^{\ast }(H(u_0))(x,0))^2
 -(\partial _{t}G^{\ast }(H(u_0))(x,T))^2){dx} 
 +\mathcal{B}_n\int _{\omega }(\partial _{t}G^{\ast}(H(u_0))(x,T))^2{dx} \\
&=-\mathcal{B}_n\int _{\omega ^{T}}\partial _{t}|\partial
_{t}G^{\ast }(H(u_0))|^2\,dx\,dt+\mathcal{B}_n\int _{\omega
}(\partial _{t}G^{\ast }(H(u_0))(x,T))^2{dx} \\
&=-2\mathcal{B}_n\int _{\omega ^{T}}\partial _{tt}^2G^{\ast }(H(u_0))
 \partial _{t}G^{\ast }(H(u_0))\,dx\,dt
 +\mathcal{B}_n\int_{\omega }(\partial _{t}G^{\ast }(H(u_0))(x,T))^2{dx}.
\end{align*}
Note, that from the definition of $G^{\ast }(H(u_0))$, we have
\[
\partial _{t}G^{\ast }(H(u_0))(x,T)=(\mathcal{B}_n+N^{-1})G^{\ast
}(H(u_0))(x,T)-H(u_0)(x,T)=-H(u_0)(x,T).
\]
Thus, we transform the last equality to
\begin{align*}
&\mathcal{B}_n\int _{\omega }(\partial _{t}G^{\ast
}(H(u_0))(x,0))^2{dx} \\
&=-2\mathcal{B}_n\int _{\omega ^{T}}\partial _{tt}^2G^{\ast
}(H(u_0))\partial _{t}G^{\ast }(H(u_0))\,dx\,dt + \mathcal{B}_n\int
_{\omega }|H(u_0)(x,T)|^2{dx}.
\end{align*}
Now, we put this into the original expression and obtain
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))u_0\,dx\,dt= \\
&=\int _{\omega ^{T}}(\partial _{tt}^2G^{\ast }(H(u_0)))^2\,dx\,dt+%
\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int _{\omega ^{T}}(\partial
_{t}G^{\ast }(H(u_0)))^2\,dx\,dt \\
&\quad -2\mathcal{B}_n\int _{\omega ^{T}}\partial _{tt}^2G^{\ast
}(H(u_0))\partial _{t}G^{\ast }(H(u_0))\,dx\,dt +\mathcal{B}_n\int
_{\omega }|H(u_0)(x,T)|^2{dx}.
\end{align*}
We transform it as follows,
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))u_0\,dx\,dt \\
&=\int _{\omega ^{T}}(\mathcal{B}_n\partial _{t}G^{\ast
}(H(u_0))-\partial _{tt}^2G^{\ast }(H(u_0)))^2\,dx\,dt+N^{-1}\mathcal{B%
}_n\int _{\omega ^{T}}(\partial _{t}G^{\ast }(H(u_0)))^2\,dx\,dt
\\
&\quad +\mathcal{B}_n\int _{\omega }|H(u_0)(x,T)|^2{dx}=\int
_{\omega ^{T}}(\partial _{t}H(u_0)-N^{-1}\partial _{t}G^{\ast
}(H(u_0)))^2\,dx\,dt \\
&\quad +N^{-1}\mathcal{B}_n\int _{\omega ^{T}}(\partial _{t}G^{\ast
}(H(u_0)))^2\,dx\,dt+\mathcal{B}_n\int _{\omega
}|H(u_0)(x,T)|^2{dx}.
\end{align*}
Using again the definition of operators $G^{\ast }$ and $H$, we have
\begin{align*}
\partial _{t}H(u_0)-N^{-1}\partial _{t}G^{\ast }(H(u_0))
&=u_0-(\mathcal{%
B}_n+N^{-1})H(u_0)+N^{-1}(\mathcal{B}_n+N^{-1})G^{\ast }(H(u_0)) \\
&\quad +N^{-1}H(u_0)-N^{-1}(\mathcal{B}_n+N^{-1})G^{\ast }(H(u_0))=u_0-
\mathcal{B}_nH(u_0).
\end{align*}
Thus,
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))u_0\,dx\,dt \\
&=\int _{\omega ^{T}}|u_0-\mathcal{B}_nH(u_0)|^2\,dx\,dt+N^{-1}%
\mathcal{B}_n\int _{\omega ^{T}}(\partial _{t}G^{\ast
}(H(u_0)))^2\,dx\,dt+\mathcal{B}_n\int _{\omega
}|H(u_0)(x,T)|^2{dx}.
\end{align*}

Using that $H(G^{\ast }(p_0))$ is a solution \eqref{second order ODE} we obtain
\begin{align*}
\int _{\omega ^{T}}H(G^{\ast }(p_0))p_0\,dx\,dt 
&=-\int_{\omega ^{T}}H(G^{\ast }(p_0))(\partial _{tt}^2H(G^{\ast
}(p_0))-\mathcal{B}_n(\mathcal{B}_n+N^{-1})H(G^{\ast }(p_0)))\,dx\,dt
\\
&=\int _{\omega ^{T}}(\partial _{t}H(G^{\ast }(p_0)))^2\,dx\,dt-\int
_{\omega }\partial _{t}(H(G^{\ast }(p_0)))(x,T)H(G^{\ast
}(p_0))(x,T){dx} \\
&\quad +\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int _{\omega ^{T}}(H(G^{\ast
}(p_0)))^2\,dx\,dt \\
&=\int _{\omega ^{T}}(\partial _{t}H(G^{\ast}(p_0)))^2\,dx\,dt 
 +\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int _{\omega ^{T}}(H(G^{\ast
}(p_0)))^2\,dx\,dt \\
&\quad +(\mathcal{B}_n+N^{-1})\int _{\omega }(H(G^{\ast }(p_0))(x,T))^2{dx}.
\end{align*}
Thus, 
\begin{align*}
\lim _{\varepsilon \to 0}2J_{\varepsilon }(v_{\varepsilon})
&=\| \nabla u_0\| _{L^2(Q^{T})}^2+\| u_0(x,T)\|_{L^2(\Omega )}^2 \\
&\quad +\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}|u_0-\mathcal{B}_nM(u_0)|^2\,dx\,dt+\mathcal{A}_n\mathcal{B}%
_n\int _{\Omega \setminus \overline{\omega }}|M(u_0)(x,T)|^2{dx}
\\
&\quad +\mathcal{A}_n\int _{\omega ^{T}}|u_0-\mathcal{B}_nH(u_0)|^2%
\,dx\,dt+\mathcal{A}_n\mathcal{B}_n\int _{\omega
}|H(u_0)(x,T)|^2{dx} \\
&\quad +N^{-1}\mathcal{A}_n\mathcal{B}_n\int _{\omega ^{T}}(\partial
_{t}G^{\ast }(H(u_0)))^2\,dx\,dt+N^{-1}\mathcal{A}_n\mathcal{B}_n\int
_{\omega ^{T}}(\partial _{t}H(G^{\ast }(p_0)))^2\,dx\,dt \\
&\quad +N^{-1}\mathcal{A}_n\mathcal{B}_n^2(\mathcal{B}_n+N^{-1})\int%
_{\omega ^{T}}(H(G^{\ast }(p_0)))^2\,dx\,dt \\
&\quad +N^{-1}\mathcal{A}_n\mathcal{B}_n(\mathcal{B}_n+N^{-1})\int
_{\omega }(H(G^{\ast }(p_0))(x,T))^2{dx} \\
&= J_0(-N^{-1}H(G^*(p_0))).
\end{align*}
This completes the proof.
\end{proof}

Finally, we will show that  system \eqref{homogenized system} is related
to the limit functional and the limit optimal control $v_0$.

\begin{proof}[Proof of Theorem \ref{Thm convergence controls}]
If $v_0$ is the optimal control, then for any admissible function $v\in U_{\rm ad}$ 
we have
\[
J_0'(v_0)v = \lim _{\lambda \to 0}\frac{%
J_0(v_0+\lambda v)-J_0(v_0)}{\lambda }=0.
\]
We set $\theta =(u_0(v_0+\lambda v)-u_0(v_0))/\lambda $. The
function $\theta $ is a solution to the  problem
\begin{equation}
\begin{gathered}
\begin{aligned}
&\partial _{t}\theta -\Delta \theta +\mathcal{A}_n(\theta -\mathcal{B}
_nH(\theta ))\chi _{\omega ^{T}}
+\mathcal{A}_n(\theta -\mathcal{B}_nM(\theta ))\chi _{(\Omega \setminus
\overline{\omega })\times (0,T)} \\
&=\mathcal{A}_n\mathcal{B}_nv\chi_{\omega ^{T}}, \quad  (x,t)\in Q^{T},
\end{aligned} \\
\theta (x,0)=0, \quad x\in \Omega , \\
\theta (x,t)=0, \quad  (x,t)\in \Gamma ^{T}.
\end{gathered} \label{limit theta prob}
\end{equation}
Note that because of the linearity of the operator $M$ and $H$, we have
\begin{gather*}
(H(u_0(v_0 + \lambda v) - H(u_0(v_0)))/\lambda = H(\theta), \\
(M(u_0(v_0 + \lambda v) - M(u_0(v_0)))/\lambda = M(\theta).
\end{gather*}
Thus,
\begin{equation}
\begin{aligned}
J_0'(v_0)v
&=\int _{Q^{T}}\nabla \theta \nabla u_0(v_0)%
\,dx\,dt+\int _{\Omega }\theta (x,T)u_0(v_0)(x,T){dx}    \\
&\quad +\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(\theta -\mathcal{B}_nM(\theta ))(u_0(v_0)-\mathcal{B}%
_nM(u_0(v_0)))\,dx\,dt    \\
&\quad +\mathcal{A}_n\mathcal{B}_n\int _{\Omega \setminus \overline{%
\omega }}M(\theta )(x,T)M(u_0(v_0))(x,T){dx}    \\
&\quad +\mathcal{A}_n\int _{\omega ^{T}}(\theta -\mathcal{B}_nH(\theta
))(u_0(v_0)-\mathcal{B}_nH(u_0(v_0)))\,dx\,dt    \\
&\quad +\mathcal{A}_n\mathcal{B}_n\int _{\omega }H(\theta
)(x,T)H(u_0(v_0))(x,T){dx}    \\
&\quad +N^{-1}\mathcal{A}_n\mathcal{B}_n\int _{\omega ^{T}}\partial
_{t}G(H^{\ast }(\theta ))\partial _{t}G(H^{\ast }(u_0(v_0)))\,dx\,dt+N%
\mathcal{A}_n\mathcal{B}_n\int _{\omega ^{T}}\partial
_{t}v_0\partial _{t}v\,dx\,dt    \\
&\quad +N\mathcal{A}_n\mathcal{B}_n^2(\mathcal{B}_n+N^{-1})\int
_{\omega ^{T}}v_0v\,dx\,dt+N\mathcal{A}_n\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})\int _{\omega }v_0(x,T)v(x,T){dx}.
\end{aligned}\label{limit J deriv}
\end{equation}
Now, we use that $p_0$ is a solution to the problem
\eqref{p_0 limit prob}, that is adjoint to \eqref{limit state problem func}, and obtain
\begin{align*}
&\int _{Q^{T}}\nabla \theta \nabla u_0(v_0)\,dx\,dt\\
&=-\int _{Q^{T}}\partial _{t}p_0\theta \,dx\,dt+\int _{Q^{T}}\nabla
p_0\nabla \theta \,dx\,dt \\
&\quad +\mathcal{A}_n\int _{\omega ^{T}}(p_0-\mathcal{B}_nH^{\ast
}(p_0))\theta \,dx\,dt+\mathcal{A}_n\int _{(\Omega \setminus
\overline{\omega })\times (0,T)}(p_0-\mathcal{B}_nM^{\ast
}(p_0))\theta \,dx\,dt \\
&\quad -\mathcal{A}_n\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))\theta \,dx\,dt \\
&\quad -\mathcal{A}_n\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0-\mathcal{B}_n^2M^{\ast }(M(u_0)))\theta \,dx\,dt.
\end{align*}
As $\theta $ is a solution to  problem \eqref{limit theta prob}, we have
\begin{equation}
\begin{aligned}
\int _{Q^{T}}\nabla \theta \nabla u_0(v_0)\,dx\,dt
&=-\int _{\Omega }u_0(v_0)(x,T)\theta (x,T){dx}+\mathcal{A}_n\mathcal{B}%
_n\int _{\omega ^{T}}vp_0\,dx\,dt    \\
&\quad -\mathcal{A}_n\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0))\theta \,dx\,dt \\
&\quad -\mathcal{A}_n\int_{(\Omega \setminus \overline{\omega })
 \times (0,T)}(u_0-\mathcal{B} _n^2M^{\ast }(M(u_0)))\theta {\,dx\,dt.}
\end{aligned}  \label{nabla theta rel}
\end{equation}
Using the same transformations as in the proof of the theorem above, we have
\begin{equation}
\begin{aligned}
&\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0(v_0)-\mathcal{B}_n^2M^{\ast }(M(u_0(v_0)))\theta {\,dx\,dt%
}    \\
&=\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0(v_0)\theta -\mathcal{B}_n^2M(u_0(v_0)))M(\theta ){\,dx\,dt%
}    \\
&=\int _{(\Omega \setminus \overline{\omega })\times (0,T)}((\partial
_{t}M(u_0(v_0))+\mathcal{B}_nM(u_0(v_0)))((\partial _{t}M(\theta )
  \\
&\quad +\mathcal{B}_nM(\theta ))-\mathcal{B}_n^2M(u_0(v_0))M(\theta ){\,dx\,dt%
}    \\
&=\int _{(\Omega \setminus \overline{\omega })\times
(0,T)}(u_0(v_0)-\mathcal{B}_nM(u_0(v_0)))(\theta -\mathcal{B}%
_nM(\theta ))\,dx\,dt    \\
&\quad +\mathcal{B}_n\int _{\Omega }M(u_0(v_0))(x,T)M(\theta )(x,T){dx}.
\end{aligned} \label{theta and u rel 2}
\end{equation}
And, similarly as in the proof of the theorem above, we obtain
\begin{align*}
&\int _{\omega ^{T}}(u_0(v_0)-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0(v_0))))\theta \,dx\,dt \\
&=\int _{\omega ^{T}}(u_0(v_0)-\mathcal{B}_nH(u_0(v_0)))(%
\theta -\mathcal{B}_nH(\theta ))\,dx\,dt \\
&\quad +\mathcal{B}_n\int _{\omega }H(u_0(v_0))(x,T)H(\theta )(x,T){dx}
\\
&\quad +N^{-1}\mathcal{B}_n\int _{\omega ^{T}}\partial _{t}G(H^{\ast
}(u_0(v_0)))\partial _{t}G(H^{\ast }(\theta ))\,dx\,dt.
\end{align*}
Indeed, we use that $G^{\ast }(H(u_0))$ is the solution to  problem
\eqref{second order ODE} with boundary conditions \eqref{bound cond ODE},
and derive that
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))\theta\,dx\,dt \\
&=\int _{\omega ^{T}}\partial _{tt}^2G^{\ast }(H(u_0))(\partial
_{tt}^2G^{\ast }(H(\theta))-\mathcal{B}_n(\mathcal{B}_n+N^{-1})G^{\ast
}(H(\theta)))\,dx\,dt \\
&=\int _{\omega ^{T}}\partial _{tt}^2G^{\ast }(H(u_0))
\partial_{tt}^2G^{\ast }(H(\theta))\,dx\,dt - \mathcal{B}_n(\mathcal{B}%
_n+N^{-1})\int _{\omega ^{T}}\partial _{tt}^2G^{\ast
}(H(u_0))G^{\ast }(H(\theta))\,dx\,dt \\
&=\int _{\omega ^{T}}\partial _{tt}^2G^{\ast }(H(u_0))
\partial_{tt}^2G^{\ast }(H(\theta))\,dx\,dt + \mathcal{B}_n(\mathcal{B}%
_n+N^{-1})\int _{\omega ^{T}}\partial _{t}G^{\ast
}(H(u_0))\partial _{t}G^{\ast }(H(\theta))\,dx\,dt \\
&\quad + \mathcal{B}_n\int _{\omega }\partial _{t}G^{\ast
}(H(u_0))(x,0)\partial _{t}G^{\ast }(H(\theta))(x,0){dx}.
\end{align*}
We transform the last term as follows,
\begin{align*}
&\mathcal{B}_n\int _{\omega }\partial _{t}G^{\ast
}(H(u_0))(x,0)\partial _{t}G^{\ast }(H(\theta))(x,0){dx} \\
&=-\mathcal{B}_n\int _{\omega ^{T}}\partial _{t}(\partial
_{t}G^{\ast }(H(u_0))G^{\ast }(H(\theta)))\,dx\,dt \\
&\quad + \mathcal{B}_n\int _{\omega }\partial _{t}G^{\ast
}(H(u_0))(x,T)\partial _{t}G^{\ast }(H(\theta))(x,T){dx} \\
&=-\mathcal{B}_n\int _{\omega ^{T}}(\partial _{tt}^2G^{\ast
}(H(u_0))\partial _{t}G^{\ast }(H(\theta)) + \partial _{t} G^{\ast
}(H(u_0))\partial _{tt}^2G^{\ast }(H(\theta)))\,dx\,dt \\
&\quad + \mathcal{B}_n\int _{\omega }\partial _{t}G^{\ast
}(H(u_0))(x,T)\partial _{t}G^{\ast }(H(\theta))(x,T){dx}.
\end{align*}
Note, that from the definition of $G^{\ast }(H(u_0))$, we have
\[
\partial _{t}G^{\ast }(H(u_0))(x,T) = -H(u_0)(x,T),\quad
\partial_{t}G^{\ast}(H(\theta))(x, T) = - H(\theta)(x, T).
\]%
Thus, we transform the last equality to
\begin{align*}
&\mathcal{B}_n\int _{\omega }\partial _{t}G^{\ast
}(H(u_0))(x,0)\partial _{t}G^{\ast }(H(\theta))(x,0){dx} \\
&= -\mathcal{B}_n\int _{\omega ^{T}}(\partial _{tt}^2G^{\ast
}(H(u_0))\partial _{t}G^{\ast }(H(\theta)) + \partial _{t} G^{\ast
}(H(u_0))\partial _{tt}^2G^{\ast }(H(\theta)))\,dx\,dt \\
&\quad + \mathcal{B}_n\int _{\omega }H(u_0)(x,T) H(\theta)(x,T){dx}.
\end{align*}
Now, we put this into the original expression and obtain
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}
_n+N^{-1})G^{\ast }(H(u_0)))\theta \,dx\,dt \\
&=\int _{\omega ^{T}}\partial _{tt}^2G^{\ast }(H(u_0)) \partial
_{tt}^2G^{\ast }(H(\theta))\,dx\,dt+\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})\int _{\omega ^{T}}\partial _{t}G^{\ast }(H(u_0))
\partial _{t}G^{\ast }(H(\theta))\,dx\,dt \\
&\quad - \mathcal{B}_n\int _{\omega ^{T}}(\partial _{tt}^2G^{\ast
}(H(u_0))\partial _{t}G^{\ast }(H(\theta)) + \partial _{t} G^{\ast
}(H(u_0))\partial _{tt}^2G^{\ast }(H(\theta)))\,dx\,dt \\
&\quad +\mathcal{B}_n\int _{\omega }H(u_0)(x,T) H(\theta)(x,T){dx}.
\end{align*}
We transform it in as follows,
\begin{align*}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))\theta \,dx\,dt \\
&= \int _{\omega ^{T}}(\mathcal{B}_n\partial _{t}G^{\ast
}(H(u_0))-\partial _{tt}^2G^{\ast }(H(u_0)))(\mathcal{B}_n\partial
_{t}G^{\ast }(H(\theta))-\partial _{tt}^2G^{\ast }(H(\theta)))\,dx\,dt \\
&\quad + N^{-1}\mathcal{B}_n\int _{\omega ^{T}}\partial _{t}G^{\ast
}(H(u_0)) \partial _{t}G^{\ast }(H(\theta))\,dx\,dt +\mathcal{B}_n\int
_{\omega }H(u_0)(x,T) H(\theta)(x,T){dx} \\
&= \int _{\omega ^{T}}(\partial _{t}H(u_0)-N^{-1}\partial
_{t}G^{\ast }(H(\theta))) (\partial _{t}H(\theta)-N^{-1}\partial _{t}G^{\ast
}(H(\theta)))\,dx\,dt \\
&\quad + N^{-1}\mathcal{B}_n\int _{\omega ^{T}}\partial _{t}G^{\ast
}(H(u_0)) \partial _{t}G^{\ast }(H(\theta))\,dx\,dt +\mathcal{B}_n\int
_{\omega }H(u_0)(x,T) H(\theta)(x,T){dx}.
\end{align*}
Using again the definition of operators $G^{\ast }$ and $H$, we have
\begin{gather*}
\partial _{t}H(u_0)-N^{-1}\partial _{t}G^{\ast }(H(u_0))=u_0-\mathcal{B%
}_nH(u_0), \\
\partial _{t}H(\theta)-N^{-1}\partial _{t}G^{\ast }(H(\theta))=u_0-%
\mathcal{B}_nH(\theta).
\end{gather*}
Thus,
\begin{equation}
\begin{aligned}
&\int _{\omega ^{T}}(u_0-\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})G^{\ast }(H(u_0)))\theta\,dx\,dt    \\
&=\int _{\omega ^{T}}(u_0-\mathcal{B}_nH(u_0))(\theta - \mathcal{%
B}_nH(\theta))\,dx\,dt    \\
&\quad + N^{-1}\mathcal{B}_n\int _{\omega ^{T}}\partial _{t}G^{\ast
}(H(u_0)) \partial _{t}G^{\ast }(H(\theta))\,dx\,dt
+\mathcal{B}_n\int_{\omega }H(u_0)(x,T) H(\theta)(x,T){dx}.
\end{aligned}  \label{theta and u rel 1}
\end{equation}

Substituting relations \eqref{theta and u rel 2} and
\eqref{theta and u rel 1} into \eqref{nabla theta rel}, and then, 
putting it into \eqref{limit J deriv}, we derive
\begin{align*}
J_0'(v_0)v
&=\mathcal{A}_n\mathcal{B}_n\int _{\omega
^{T}}p_0v\,dx\,dt+N\mathcal{A}_n\mathcal{B}_n\int _{\omega
^{T}}\partial _{t}v_0\partial _{t}v\,dx\,dt \\
&\quad +N\mathcal{A}_n\mathcal{B}_n^2(\mathcal{B}_n+N^{-1})\int
_{\omega ^{T}}v_0v\,dx\,dt+N\mathcal{A}_n\mathcal{B}_n(\mathcal{B}%
_n+N^{-1})\int _{\omega }v_0(x,T)v(x,T){dx}.
\end{align*}
Integrating by parts in the second integral, we obtain
\begin{align*}
J_0'(v_0)v
&=\mathcal{A}_n\mathcal{B}_n\int _{\omega
^{T}}p_0v\,dx\,dt+N\mathcal{A}_n\mathcal{B}_n\int _{\omega
}(\partial _{t}v_0(x,T)+(\mathcal{B}_n+N^{-1})v_0(x,T))v(x,T){dx} \\
&\quad +N\mathcal{A}_n\mathcal{B}_n\int _{\omega ^{T}}(-\partial
_{tt}^2v_0+\mathcal{B}_n(\mathcal{B}_n+N^{-1})v_0)v\,dx\,dt=0,
\end{align*}
for an arbitrary admissible function $v\in U_{\rm ad}$. Hence, we obtain that 
$v_0 $ must be a solution of the  problem
\begin{gather*}
\partial _{tt}^2v_0-\mathcal{B}_n(\mathcal{B}_n+N^{-1})v_0=N^{-1}p_0   \\
v_0(x,0)=0,\quad \partial _{t}v_0(x,T)+(\mathcal{B}_n+N^{-1})v_0(x,T)=0. 
\end{gather*}
But, this is related to the problem satisfied by $H(G^{\ast }(-N^{-1}p_0))$. 
Because of the uniqueness of the solution of the problem we conclude that 
$v_0=-N^{-1}H(G^{\ast }(p_0))\chi _{\omega ^{T}}$.
\end{proof}

\begin{remark} \rm
As in \cite{DiPoShRACSAM, DiPoShAttouch}, thanks to
Remark \ref{Rem target nonzero}, by simplifying the cost functional 
$J_{\varepsilon }$ as in \cite{DiPoShRACSAM}, by making the parameter 
$N\to 0$, it seems possible to show the approximate controllability
with final observation of solutions of the limit problem.
\end{remark}

\subsection*{Acknowledgements}
J. I. D\'{\i}az was supported by project PID2023-146754NB-I00 from the
 MCIU/AEI/10.13039/501100011033 and  FEDER,
EU.MCIU/AEI/10.13039/-501100011033/ \\ FEDER, EU.

\begin{thebibliography}{10}

\bibitem{Anguiano} M. Anguiano;
\emph{Homogenization of parabolic problems with
dynamical boundary conditions of reactive-diffusive type in perforated
media}, ZAMM - J. Appl. Math. Mech. / Zeitschrift f\"{u}r Angew. Math. und
Mech, 100.10 (2020).

\bibitem{BejeDVrabie} I. Bejenaru, J. I. D\'{\i}az, I. I. Vrabie;
\emph{An abstract approximate controllability result and applications to elliptic
and parabolic systems with dynamic boundary conditions}, Electronic Journal of
Differential Equations, 2001 (2001) No. 50, 1--19.

\bibitem{Cioranescu1997AST} D. Cioranescu, F. Murat;
\emph{A Strange Term Coming from Nowhere}, 
In Topics in the Mathematical Modelling of Composite
Materials, A. Cherkaev and R. Kohn, eds., Birkh{\"{a}}user Boston, 45--93,
1997.

\bibitem{ConcaDliTi} C. Conca, J. I. D\'{\i}az, A. Li\~{n}\'{a}n,  C. Timofte;
\emph{Homogenization in Chemical Reactive Flows}, Electronic Journal of
Differential Equations, 2004 (2004) No. 40, 1--22.

\bibitem{ShDyn19} J. I. Diaz, D. Gomez-Castro, T. A. Shaposhnikova, M. N. Zubova;
\emph{A nonlocal memory strange term arising in the critical scale
homogenization of diffusion equation with a dynamic boundary condition},
Electron. J. Differential Equations, 2019 (2019) No. 77. 1--13.

\bibitem{DiGoShBook} J. I. D\'{\i}az, D. G\'{o}mez-Castro, T. A. Shaposhnikova;
\emph{Nonlinear Reaction-Diffusion Processes for Nanocomposites.
Anomalous improved homogenization}, Berlin: De Gruyter Series in Nonlinear
Analysis and Applications, V. 39. 2021.

\bibitem{DShPo-first control} J. I. D\'{\i}az, A. V. Podolskiy, T. A. Shaposhnikova;
\emph{On the convergence of controls and cost functionals in some
optimal control heterogeneous problems when the homogenization process gives
rise to some strange terms}, Journal of Mathematical Analysis and
Applications, V. 506 (2022), I. 1.

\bibitem{DiPoSh control arbt size} J. I. D\'{\i}az, A. V. Podolskiy,  T. A. Shaposhnikova;
\emph{On the homogenization of an optimal control problem in a
domain perforated by holes of critical size and arbitrary shape}, 
Doklady Mathematics, V. 105 (2022) , 6--13.

\bibitem{DiPoShBound22} J. I. D\'{\i}az, A. V. Podolskiy, T. A. Shaposhnikova;
\emph{Boundary control and homogenization: optimal climatization
through smart double skin boundaries}. Differential and Integral Equations.
35, 3-4 (2022), 191--210.

\bibitem{DiPoSh Nonlinear Corrector-control} 
J. I. D\'{\i}az, A. V. Podolskiy, T. A. Shaposhnikova;
\emph{On the corrector term in the
homogenization of the nonlinear Poisson-Robin problem giving rise to a
strange term: Application to an optimal control problem}. 
Journal of Mathematical Analysis and Applications 543.1 (2025): 128867.

\bibitem{DiPoShRACSAM} J. I. D\'{\i}az, A. V. Podolskiy, T. A. Shaposhnikova;
\emph{Think globally, act locally: approximate controllability
through homogenization of an optimal control problem with control on the
boundary of certain particles}. Rev. Real Acad. Cienc. Exactas Fis. Nat. Ser.
A-Mat. 119, 78 (2025). RACSAM, https://doi.org/10.1007/s13398-025-01744-x.
Online June 11, 2025.

\bibitem{DiPoShAttouch} J. I. D\'{\i}az, A. V. Podolskiy,  T. A. Shaposhnikova;
\emph{New unexpected limit operators for homogenizing optimal
control parabolic problems with dynamic reaction flow on the boundary of
critically scaled particles}. To appear in the Journal of Convex Analysis.

\bibitem{DiShZarbitrary} J. I. D\'{\i}az, T. A. Shaposhnikova,  M. N. Zubova;
\emph{A strange non-local monotone operator arising in the homogenization
of a diffusion equation with dynamic nonlinear boundary conditions on
particles of critical size and arbitrary shape}, Electronic Journal of
Differential Equations, Vol. 2022 (2022), No. 52, pp. 1--32.

\bibitem{Fursikov2000OptimalCO} A. V. Fursikov;
\emph{Optimal Control of Distributed Systems: Theory and Applications}, 
Boston: AMS, 2000.

\bibitem{Glow-Lions} R. Glowinski, J.-L. Lions, J. He;
\emph{Exact and Approximate Controllability for Distributed Parameter Systems: 
A Numerical Approach},
Cambridge: Cambridge University Press, 2008.

\bibitem{GoSanchezP} D. G\'{o}mez, M. Lobo, E. P\'{e}rez, E. S\'{a}nchez-Palencia;
\emph{Homogenization in perforated domains: a Stokes grill and an
adsorption process}, Applicable Analysis, V. 97, I. 16 (2018), 2893--2919.

\bibitem{Hornung-Yaeger} U. Hornung, W. J\"{a}ger;
\emph{Diffusion, convection, adsorption, and reaction of chemicals in porous media}. 
Journal of Differential Equations, V. 92 (1991), I. 2, 199-225.

\bibitem{Hurlov} E. J. Hruslov;
\emph{The Method of Orthogonal Projections and the
Dirichlet Problem in Domains With a Fine-Grained Boundary},
 Mathematics of the USSR-Sbornik, V. 17 (1972), N. 1, 37--59.

\bibitem{Iliev-Mikelik} O. Iliev, A. Mikelic, T. Prill, A. Sherly;
\emph{Homogenization approach to the upscaling of a reactive flow through
particulate filters with wall integrated catalyst}. Adv. Water Resour, V. 146,
2020.

\bibitem{Kaizu1} S. Kaizu;
\emph{The Poisson equation with semilinear boundary
conditions in domains with many tiny holes},
 J. Fac. Sci. Univ. Tokyo Sect. IA Math., V. 36 (1989), 43--86.

\bibitem{LionsOpt} J. L. Lions;
\emph{Contr\^{o}le Optimal de Syst\`{e}mes Gouvern\'{e}s par des Equations 
aux Deriv\'{e}es Partielles}, Paris: Dunod, 1968.

\bibitem{Nolasco-Survey} E. Nolasco, et al;
\emph{Optimal control in chemical engineering: Past, present and future}. 
Computers \& Chemical Engineering, V. 155 (2021).

\bibitem{OlSh95} O. A. Oleinik, T. A. Shaposhnikova;
\emph{On homogenization problems for the Laplace operator in partially 
perforated domains with Neumann's condition on the boundary of cavities}, 
Atti della Accademia Nazionale dei Lincei. Classe di Scienze Fisiche, 
Matematiche e Natu-rali.
Rendiconti Lincei. Matematica e Applicazioni, V. 6 (1995). N. 3. 133--142.

\bibitem{PoShControl20} A. V. Podolskiy, T. A. Shaposhnikova;
\emph{Optimal Control and ``Strange'' Term Arising from Homogenization of the Poisson 
Equation in the Perforated Domain with the Robin-type Boundary Condition in the Critical
Case}, Doklady Mathematics, V. 102 (2020), 497--501.

\bibitem{SJPZ02} J. Saint Jean Paulin, H. Zoubairi;
\emph{Optimal control and ``strange term''  for a Stokes problem in
perforated domains}, Port. Math., V. 59 (2002) , 161--178.

\bibitem{Shaposhnikova2022HomogenizationOT} 
T. A. Shaposhnikova, A. V. Podolskiy;
\emph{Homogenization of the optimal control problem for the Dirichlet
cost functional and the Poisson state problem with rapidly alternating
boundary conditions in critical case}, Revista de la Real Academia de
Ciencias Exactas, F\'{\i}sicas y Naturales. Serie A. Matem\'{a}ticas, V.
116 (2022), N. 174.

\bibitem{Timof} C. Timofte;
\emph{Parabolic problems with dynamical boundary
conditions in perforated media}, 
Math. Modeling and Analysis, V. 8 (2003), 337--350.

\bibitem{Trolstz} F. Troltzsch;
\emph{Optimal control of partial differentialequations}, 
Providence, RI.: American Mathematical Society, 2010.

\bibitem{Upreti book} S. R. Upreti;
\emph{Optimal control for chemical engineers},
Boca Raton: CRC Press, 2013.

\bibitem{ZuShDiffEq} M. N. Zubova, T. A. Shaposhnikova;
\emph{Homogenization of boundary value problems in perforated domains with 
the third boundary condition and resulting change in the character of the 
nonlinearity in the problem}, Differential Equations, V. 47 )2011), 78--90.

\bibitem{ZuShDynDok19} M. N. Zubova, T. A. Shaposhnikova;
\emph{Homogenization Limit for the Diffusion Equation in a Domain Perforated along
(n-1)-Dimensional Manifold with Dynamic Conditions on the Boundary of the
Perforations: Critical Case}, Doklady Mathematics, V. 99 (2019), N.3, 245--251.

\end{thebibliography}

\end{document}
