elliptic PDE-constrained optimization problem with box constraints on the control: $$\label{eqn:orginal problems}
\begin{aligned}
\min \limits_{(y,u)\in Y\times U}^{}\ \ &J(y,u)=\frac{1}{2}\|y-y_d\|_{L^2(\Omega)}^{2}+\frac{\alpha}{2}\|u\|_{L^2(\Omega)}^{2} \\
{\rm s.t.}\ \ \ \ \ \ &Ly=u+y_{r}\ \ \ \ \ \mathrm{in} \  \Omega, \\
&\ \ y=0 \quad  \ \ \ \ \ \ \ \  \ \mathrm{on} \ \partial\Omega,\\
\qquad \qquad\quad &\ \ u\in  U_{ad}=\{v(x)|a\leq v(x)\leq b, {\rm a.e }\  \mathrm{on}\ \Omega\}\subseteq U,
\end{aligned} \tag{$\mathrm{P}$}$$ where $Y:=H_{0}^{1}(\Omega), U:=L^{2}(\Omega), \Omega\subseteq\mathbb{R}^{n}(n=2,3)$ is a convex, open and bounded domain with $C^{1,1}$- or polygonal boundary; the desired state $y_{d}\in L^{2}(\Omega)$ and the source term $y_{r}\in L^{2}(\Omega)$ are given; parameters $\alpha>0$, $-\infty<a<b<+\infty$; $L$ is the uniformly elliptic differential operater given by $$Ly:=-\sum\limits_{i,j=1}^{n}(a_{ij}y_{x_{i}})_{{x_{j}}}+c_{0}y,$$ where $a_{ij},c_{0}\in L^{\infty}(\Omega),c_{0}\geqslant 0,a_{ij}=a_{ji},\sum\limits_{i,j=1}^{n}a_{ij}(x)\xi_{i}\xi_{j}\geqslant\theta\|\xi\|^{2}, \ \rm a.a. \ x\in \Omega, \forall \xi\in \mathbb{R}^{n}.$

By the classical conclusion in the PDE theory, for a given $y_{r}\in L^{2}(\Omega)$ and every $u\in L^{2}(\Omega)$, elliptic PDEs involved in the (\[eqn:orginal problems\]) $$\label{eqn:pde}
\begin{aligned}
Ly&=u+y_{r}\quad   \mathrm{in}\ \Omega, \\
y&=0\quad \quad \quad \ \mathrm{on}\ \partial\Omega,\\
\end{aligned}$$ has a unique weak solution $y=y(u):=S(u+y_{r})$, where $S:L^{2}(\Omega)\rightarrow H_{0}^{1}(\Omega)$ denotes the solution operator. It is well-defined and is a continuous linear injective operator[@Kinderlehrer1980An Theorem B.4].

As is known to all, there are two possible approaches to tackle optimization problems with 