Asymptotics of a solution to a problem of optimal boundary control with two small cosubordinate parameters

We consider a problem of optimal boundary control for solutions of an elliptic type equation in a bounded domain with smooth boundary with a small coefficient at the Laplace operator, a small coefficient, cosubordinate with the first, at the boundary condition, and integral constraints on the control:
$$ \left\{
\begin {array}{ll} \displaystyle \mathcal {L}_\varepsilon \mathop {:=} \nolimits - \varepsilon^2 \Delta z + a(x) z = f(x), & \displaystyle x\in \Omega,\quad z \in H^1 (\Omega), \\[3ex] \displaystyle l_{\varepsilon,\beta} z\mathop {:=} \nolimits \varepsilon^\beta \frac{\partial z}{\partial n} = g(x) + u(x), & x\in\Gamma, \end {array}
 \right. $$

$$ J(u) \mathop {:=} \nolimits \|z-z_d\|^2 + \nu^{-1}|||u|||^2 \to \inf, \quad u \in \mathcal {U}, $$
where $0<\varepsilon\ll 1$, $\beta\geqslant 0$, $\beta\in\mathbb{Q}$, $\nu>0,$ $H^1 (\Omega)$ is the Sobolev function space, $\partial z/\partial n$ is the derivative of $z$ at the point $x\in\Gamma$ in the direction of the outer (with respect to the domain $\Omega$) normal,
$$
\begin {array}{c} \displaystyle a(\cdot), f(\cdot) \in C^\infty(\overline{\Omega}), \quad g(\cdot)\in C^\infty(\Gamma),\quad \forall\, x\in \overline{\Omega}\quad a(x)\geqslant \alpha^2>0, \\[2ex] \displaystyle \mathcal {U} = \mathcal {U}_1,\quad \mathcal {U}_r\mathop {:-} \nolimits \{u(\cdot)\in L_2(\Gamma)\colon |||u||| \leqslant r \}. \end {array}
$$
Here $\|\cdot\|$ and $|||\cdot|||$ are the norms in the spaces $L_2(\Omega)$ and $L_2(\Gamma)$, respectively. We find the complete asymptotic expansion of the solution of the problem in the powers of the small parameter in the case where $0<\beta<3/2$.