us define $\bar \Omega_{h}=\bigcup_{T\in \mathcal{T}_{h}}\bar T$ and let $\Omega_{h}\subseteq \Omega$ and $\Gamma_{h}$ denote its interior and its boundary, respectively. In the case that $\Omega$ is a convex polyhedral domain, we have $\Omega=\Omega_{h}$. In the case that $\Omega$ is a domain with a $C^{1,1}$- boundary $\Gamma$, we assume that $\bar \Omega_{h}$ is convex and that all boundary vertices of $\bar \Omega_{h}$ are contained in $\Gamma$, such that $$|\Omega\backslash\Omega_{h}|\leq ch^{2},$$ where $|\cdot|$ denotes the measure of the set, and $c>0$ is a constant.

For the state variable $y$ and the control variable $u$, we choose the same discretized state space and discretized control space defined by $$\begin{aligned}
&Y_{h}:=\{y_{h}\in C( \bar{\Omega})|y_{{h}|T}\in \mathcal{P}_{1}, \forall{T\in \mathcal{T}_{h},y_{h}=0 \ {\rm in} \ \bar{\Omega} \setminus\Omega_{h}}\},\\
&U_{h}:=\{u_{h}\in C( \bar{\Omega})|u_{{h}|T}\in \mathcal{P}_{1}, \forall{T\in \mathcal{T}_{h},u_{h}=0 \ {\rm in} \ \bar{\Omega} \setminus\Omega_{h}}\},\end{aligned}$$ where $\mathcal{P}_{1}$ denotes the space of polynomials of degree less than or equal to 1. For given source term $y_{r}$ and $u\in L^{2}(\Omega)$, let the discretized state associated with $u$ denoted by $y_{h}(u)$, which is defined as the unique solution for the discretized weak formulation of the state equation (\[eqn:pde\]): $$\label{eqn:weak formulation}
\int_{\Omega_{h}}\left(\sum_{i,j=1}^{n}a_{ij}(y_{h})_{x_{i}}(v_{h})_{x_{j}}+c_{0}y_{h}v_{h}\right)dx=\int_{\Omega_{h}}(u+y_{r})v_{h}dx, \ \ \forall v_{h}\in Y.$$ Moreover, $y_{h}$ can be expressed by $y_{h}(u)=S_{h}(u+y_{r})$, where $S_{h}$ denotes the discretized version 