On a class of semilinear fractional elliptic equations

involving outside Dirac data

Huyuan<EMAIL_ADDRESS>Hichem<EMAIL_ADDRESS>Ying
<EMAIL_ADDRESS>
###### Abstract

The purpose of this article is to give a complete study of the weak solutions
of the fractional elliptic equation

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&\ {\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}+u^{p}}u=\delta_{0}&\
{\rm in}\ \ \mathbb{R}^{N}\setminus B_{1}(e_{N}),\end{array}$ (0.1)

where $p\geq 0$, $(-\Delta)^{\alpha}$ with $\alpha\in(0,1)$ denotes the
fractional Laplacian operator in the principle value sense, $B_{1}(e_{N})$ is
the unit ball centered at $e_{N}=(0,\cdots,0,1)$ in $\mathbb{R}^{N}$ with
$N\geq 2$ and $\delta_{0}$ is the Dirac mass concentrated at the origin. We
prove that problem (0.1) admits a unique weak solution when
$p>1+\frac{2\alpha}{N}$. Moreover, if in addition $p\geq\frac{N+2}{N-2}$, the
weak solution vanishes as $\alpha\to 1^{-}$. We also show that problem (0.1)
doesn’t have any weak solution when $p\in[0,1+\frac{2\alpha}{N}]$. These
results are very surprising since there are in total contradiction with the
classical setting, i.e.

$\begin{array}[]{lll}-\Delta u+u^{p}=0&\ {\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-\Delta+u^{p}}u=\delta_{0}&\ {\rm in}\ \
\mathbb{R}^{N}\setminus B_{1}(e_{N}),\end{array}$

for which it has been proved that there are no solutions for
$p\geq\frac{N+1}{N-1}$.

Key words: Fractional Laplacian; Dirac mass; Weak solution; Existence;
Uniqueness.

MSC2010: 35R06, 35A01, 35J66.

##  1  Introduction

Fractional PDEs have gained tremendous interest, not only from mathematicians
but also from physicists and engineering, during the last years. This is
essentially due to their widespread domains of applications. In fact the
fractional Laplacian arises in many areas including medicine [12], bio-
engineering [19, 20, 21, 22], relativistic physics[1, 17, 18], Modeling
populations [29], flood flow, material viscoelastic theory, biology and
earthquakes. It is also particularly relevant to study some situations, in
which the fractional Laplacian is involved in PDEs, featuring irregular data
such that those phenomena describing source terms which are concentrated at
points. In our context, the source is placed outside the unit ball
$B_{1}(e_{N})$. This generates long-term interactions and short-term
interactions, described by the nonlocal operator $(-\Delta)^{\alpha}$ and the
nonlinear absorption $u^{p}$ respectively. $(-\Delta)^{\alpha}$ has also a
probabilistic interpretation, related to the above one. It is the
$\alpha-$stable subordinated infinitesimal killed Brownian motion.

Let $B_{1}(e_{N})$ be the unit ball in $\mathbb{R}^{N}(N\geq 2)$ with center
$e_{N}=(0,\cdots,0,1)$ and $\delta_{0}$ be the Dirac mass concentrated at the
origin. Our main objective in this article is to investigate the existence,
nonexistence and uniqueness of positive weak solutions of the semilinear
fractional equation

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&\ {\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}+u^{p}}u=\delta_{0}&\
{\rm in}\ \ \mathbb{R}^{N}\setminus B_{1}(e_{N}),\end{array}$ (1.1)

where $p\geq 0$ and the fractional Laplacian $(-\Delta)^{\alpha}$ with
$\alpha\in(0,1)$ is defined by

$(-\Delta)^{\alpha}u(x)=c_{N,\alpha}\lim_{\epsilon\to
0^{+}}(-\Delta)_{\epsilon}^{\alpha}u(x),$

where

$c_{N,\alpha}=\left(\int_{\mathbb{R}^{N}}\frac{1-\cos(z_{1})}{|z|^{N+2\alpha}}dz\right)^{-1}$
(1.2)

with $z=(z_{1},\cdots,z_{N})\in\mathbb{R}^{N}$ and

$(-\Delta)_{\epsilon}^{\alpha}u(x)=-\int_{\mathbb{R}^{N}\setminus
B_{\epsilon}(x)}\frac{u(z)-u(x)}{|z-x|^{N+2\alpha}}dz.$

In 1991, a fundamental contribution to semilinear elliptic equations involving
measures as boundary data is due to Gmira and Véron [15], where they studied
the existence and uniqueness of weak solutions for

$\begin{array}[]{lll}-\Delta u+h(u)=0&{\rm in}\quad\Omega,\\\\[5.69054pt]
\phantom{-----}u=\mu&{\rm on}\quad\partial\Omega,\end{array}$ (1.3)

where $\Omega$ is a bounded $C^{2}$ domain and $\mu$ is a bounded Radon
measure defined in $\partial\Omega$. A function $u$ is said to be a weak
solution of (1.3) if $u\in L^{1}(\Omega)$, $h(u)\in L^{1}(\Omega,\rho dx)$ and

$\int_{\Omega}[u(-\Delta)\xi+h(u)\xi]dx=\int_{\partial\Omega}\frac{\partial\xi(x)}{\partial\vec{n}_{x}}d\mu(x),\quad\forall\xi\in
C^{1.1}_{0}(\Omega),$ (1.4)

where $\rho(x)={\rm dist}(x,\partial\Omega)$ and $\vec{n}_{x}$ denotes the
unit inward normal vector at a point $x$. Gmira and Véron proved that the
problem (1.3) admits a unique weak solution when $h$ is a continuous and
nondecreasing function satisfying

$\int_{1}^{\infty}[h(s)-h(-s)]s^{-1-\frac{N+1}{N-1}}ds<+\infty.$ (1.5)

The weak solution of (1.3) is approached by the classical solutions of (1.3)
when $\mu$ is replaced by a sequence of regular functions $\\{\mu_{n}\\}$,
which converge to $\mu$ in the distribution sense. Furthermore, they showed
that there is no weak solution of (1.3) when $\mu=\delta_{x_{0}}$ with
$x_{0}\in\partial\Omega$ and $h(s)=|s|^{p-1}s$ with $p\geq\frac{N+1}{N-1}$.
Later on, this subject has been vastly expanded in recent works, see the
papers of Marcus and Véron [23, 24, 25, 26], Bidaut-Véron and Vivier [3] and
references therein.

In the fractional setting, the equivalent of (1.3) when $\mu=\delta_{x_{0}}$
has been considered in [7], where the authors proved that the weak solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+g(u)=k\frac{\partial^{\alpha}\delta_{x_{0}}}{\partial\vec{n}^{\alpha}}&{\rm
in}\quad\ \ \bar{\Omega},\\\\[8.53581pt] \phantom{------\ \ }u=0&{\rm
in}\quad\ \ \bar{\Omega}^{c}\end{array}$ (1.6)

is approximated by the weak solutions, as $t\to 0^{+}$, of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+g(u)=kt^{-\alpha}\delta_{x_{0}+t\vec{n}_{x_{0}}}&{\rm
in}\quad\ \ \bar{\Omega},\\\\[8.53581pt] \phantom{------\ \ }u=0&{\rm
in}\quad\ \ \bar{\Omega}^{c}.\end{array}$

More precisely, in the fractional setting,
$\frac{\partial^{\alpha}\delta_{x_{0}}}{\partial\vec{n}^{\alpha}}$ plays the
same role of $u=\delta_{x_{0}}$ on $\partial\Omega$ in (1.3). Our purpose in
this article is to study the solution of

$(-\Delta)^{\alpha}u+u^{p}=0\quad\ {\rm in}\quad B_{1}(e_{N})$

when the exact Dirac mass concentrated at the origin is considered. Our main
idea is to make use of nonlocal properties of the fractional Laplacian to move
the Dirac mass at $-te_{N}$ when $t\to 0^{+}$ and we then proceed by
approximation techniques. Before giving our main results, we must first give
an appropriate definition of weak solution of (1.1). It is then worth to
mention two important results. The equation

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&\ {\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=f&\ {\rm in}\ \
\mathbb{R}^{N}\setminus B_{1}(e_{N}),\end{array}$ (1.7)

where $f\in C_{0}(\mathbb{R}^{N}\setminus B_{1}(e_{N}))$, admits a unique
classical solution $u_{f}$, see [6, Theorem 2.5]. Furthermore, let
$\tilde{u}_{f}=u_{f}$ in $B_{1}(e_{N})$ and $\tilde{u}_{f}=0$ in
$\mathbb{R}^{N}\setminus B_{1}(e_{N})$, then $\tilde{u}_{f}$ is the unique
classical solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=c_{N,\alpha}\int_{\mathbb{R}^{N}\setminus
B_{1}(e_{N})}\frac{f(y)}{|x-y|^{N+2\alpha}}dy&\ {\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=0&\ {\rm in}\ \
\mathbb{R}^{N}\setminus B_{1}(e_{N})\end{array}$ (1.8)

and satisfies the identity:

$\int_{B_{1}(e_{N})}\left[u(x)(-\Delta)^{\alpha}\xi(x)+u^{p}(x)\xi(x)\right]dx=c_{N,\alpha}\int_{B_{1}(e_{N})}\int_{\mathbb{R}^{N}\setminus
B_{1}(e_{N})}\frac{\xi(x)f(y)}{|x-y|^{N+2\alpha}}dydx,$

for any $\xi\in C^{\infty}_{0}(B_{1}(e_{N}))$. Let us mention that
$C_{0}^{\infty}(B_{1}(e_{N}))$ is the space of test functions $\xi\in
C^{\infty}(\mathbb{R}^{N})$ with support in $B_{1}(e_{N})$.

Inspired by above identity, we give the definition of weak solution to (1.1)
as follows.

###### Definition 1.1

We say that $u$ is a weak solution of (1.1) if $u\in L^{1}(B_{1}(e_{N}))$,
$u^{p}\in L^{1}_{\rm loc}(B_{1}(e_{N}))$ and

$\int_{B_{1}(e_{N})}\left[u(x)(-\Delta)^{\alpha}\xi(x)+u^{p}(x)\xi(x)\right]dx=\int_{B_{1}(e_{N})}\xi(x)\Gamma_{0}(x)dx,\quad\forall\xi\in
C^{\infty}_{0}(B_{1}(e_{N})),$ (1.9)

where

$\Gamma_{0}(x)=\frac{c_{N,\alpha}}{|x|^{N+2\alpha}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\}.$ (1.10)

It is well known that the definition of the weak solution heavily depends on
the test functions space and the best function space is the one which enables
us to get the ”strongest” weak solution. In [9, 10], semilinear fractional
equations with measures has been studied via the test functions space
$\mathbb{X}_{\alpha,\Omega}\subset C(\mathbb{R}^{N})$ for a $C^{2}$ bounded
open domain $\Omega$, where $\mathbb{X}_{\alpha,\Omega}$ is the space of
functions $\xi$ satisfying:

(1) ${\rm supp}(\xi)\subset\bar{\Omega}$;

(2) $(-\Delta)^{\alpha}\xi(x)$ exists for all $x\in\Omega$ and
$|(-\Delta)^{\alpha}\xi(x)|\leq C$ for some $C>0$;

(3) there exist $\varphi\in L^{1}(\Omega,\rho^{\alpha}dx)$ and
$\epsilon_{0}>0$ such that $|(-\Delta)_{\epsilon}^{\alpha}\xi|\leq\varphi$
a.e. in $\Omega$, for all $\epsilon\in(0,\epsilon_{0}]$,

where $\rho(x)=dist(x,\partial\Omega)$.
The test functions space $C_{0}^{\infty}(B_{1}(e_{N}))$ has stronger topology
than $\mathbb{X}_{\alpha,B_{1}(e_{N})}$ does, the weak solution in Definition
1.1 with test functions space $\mathbb{X}_{\alpha,B_{1}(e_{N})}$ would be
stronger than the one with test functions space
$C_{0}^{\infty}(B_{1}(e_{N}))$. It is then worth to mention that the test
functions space $C_{0}^{\infty}(B_{1}(e_{N}))$ could not be replaced to the
test functions space $\mathbb{X}_{\alpha,B_{1}(e_{N})}$ in our setting. For
example, $\xi_{0}:=\mathbb{G}_{\alpha}[1]\in\mathbb{X}_{\alpha,B_{1}(e_{N})}$,
but (1.9) does not hold for $\xi_{0}$, where $G_{\alpha}$ denotes the Green
kernel of $(-\Delta)^{\alpha}$ in $B_{1}(e_{N})\times B_{1}(e_{N})$ and
$\mathbb{G}_{\alpha}$ is the Green operator defined as

$\mathbb{G}_{\alpha}[f](x)=\int_{B_{1}(e_{N})}G_{\alpha}(x,y)f(y)dy,\quad f\in
L^{1}(B_{1}(e_{N}),\rho^{\alpha}dx).$ (1.11)

Let us state our existence result.

###### Theorem 1.1

Assume that $\alpha\in(0,1)$ and $p>1+\frac{2\alpha}{N}$. Then there exists a
unique nonnegative weak solution $u_{\alpha,p}$ of (1.1) such that for some
$c_{1}>1$, we have

$0<u_{\alpha,p}(x)\leq c_{1}|x|^{-\frac{N+2\alpha}{p}},\qquad\forall x\in
B_{1}(e_{N})$ (1.12)

and

$\frac{1}{c_{1}}t^{-\frac{N+2\alpha}{p}}\leq u_{\alpha,p}(te_{N})\leq
c_{1}t^{-\frac{N+2\alpha}{p}},\qquad\forall t\in(0,1).$ (1.13)

###### Remark 1.1

$(i)$ The existence results are very surprising as they are in total different
from the Laplacian case, where (1.3) with $\mu=\delta_{0}$ has a weak solution
only when $p<\frac{N+1}{N-1}$.

$(ii)$ From (1.13), the singularity is only near the origin. We also notice
that

$u_{\alpha,p}^{p}(te_{N})\geq\frac{1}{c_{1}}t^{-(N+2\alpha)},\qquad\forall
t\in(0,1),$

which implies that the absorption nonlinearity $u^{p}$ plays a primary role in
the equation (1.1). While the absorption nonlinearity always plays a second
role in a measure framework.

$(iii)$ The uniqueness cannot directly follow Kato’s inequality [9,
Proposition 2.4] since it has been built in the framework of the test
functions space $\mathbb{X}_{\alpha,B_{1}(e_{N})}$. In this paper, as
mentioned above, $C_{0}^{\infty}(B_{1}(e_{N}))$ is the appropriate test
functions space. This will give birth to a lot of technical difficulties to
prove the existence, nonexistence and uniqueness of weak solutions of (1.1).

If $p>1+\frac{2\alpha}{N}$, the weak solution $u_{\alpha,p}$ of (1.1) is
approximated by the unique solution $u_{s}$ $(s\in(0,1))$ of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&{\rm in}\quad
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=\delta_{-se_{N}}&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}).\end{array}$ (1.14)

When $p\in[0,1+\frac{2\alpha}{N}]$, we will prove that $\\{u_{s}\\}$ blows up
everywhere in $B_{1}(e_{N})$ as $s\to 0^{+}$, therefore, we can deduce the
nonexistence of weak solutions of (1.1) when $p\leq 1+\frac{2\alpha}{N}$. More
precisely, we have the following results.

###### Theorem 1.2

Assume that $\alpha\in(0,1)$ and $0\leq p\leq 1+\frac{2\alpha}{N}$. Then
problem (1.1) does not have any weak solution.

In the proof of Theorem 1.2, we will first need to prove the crucial estimate

$u_{\alpha,p}^{p}(x)\geq c_{2}|x|^{-(N+2\alpha)},\quad\forall
x\in\mathcal{C},$

where $\mathcal{C}=\\{x\in\mathbb{R}^{N}:\ \exists t\in(0,1)\ {\rm s.t.}\
|x-te_{N}|<\frac{t}{8}\\}$ is a cone in $B_{1}(e_{N})$. We combine the
symmetry property and decreasing monotonicity in our proof of the
nonexistence. This phenomena is due to the nonlocal characteristic of
fractional Laplacian that requires the functions to be in $L^{1}_{\rm
loc}(\mathbb{R}^{N})$.

Finally, our interest is to study the asymptotic behavior of $u_{\alpha,p}$ as
$\alpha$ goes to $1^{-}$.

###### Theorem 1.3

Assume that $\alpha\in(0,1)$, $p\geq\frac{N+1}{N-1}$ and $u_{\alpha,p}$ is the
unique weak solution of problem (1.1). Then $\\{u_{\alpha,p}\\}_{\alpha}$
vanishes as $\alpha\to 1^{-}$.

###### Remark 1.2

$(i)$ for $p\geq\frac{N+2}{N-2}$, a sequence of barrier functions, which
converge to 0 locally in $B_{1}(0)$, could be constructed directly to control
$\\{u_{\alpha,p}\\}_{\alpha}$;
$(ii)$ for $\frac{N+1}{N-1}\leq p<\frac{N+2}{N-2}$, Theorem 3.1 in [15] is
involved to control $\\{u_{\alpha,p}\\}_{\alpha}$;
$(iii)$ for $p<1+\frac{2}{N}$, there exists $\alpha_{p}\in(0,1)$ such that
$p\leq 1+\frac{2\alpha}{N}$ for $\alpha\in(\alpha_{p},1)$, there is no weak
solution for problem (1.1) from Theorem 1.2;
$(iv)$ for $p\in[1+\frac{2}{N},\frac{N+1}{N-1})$, it is still open for the
limit of $\\{u_{\alpha,p}\\}_{\alpha}$ as $\alpha\to 1^{-}$.

In Section 2, we treat the problem (1.14). When the Dirac mass concentrates at
point $-se_{N}$ away from $\bar{\Omega}$, we build the existence, uniqueness
weak solution $u_{s}$ of (1.14) and show how the Dirac mass is transformed
into the nonhomogeneous term. In this case, the test functions space could be
improved into $\mathbb{X}_{\alpha,B_{1}(e_{N})}$, since the solution has no
singularity in $\bar{\Omega}$.

In Section 3, we give a detailed account of the procedure enabling us to move
the singular points $\\{-se_{N}\\}$ to the origin. The first difficulty arises
from the fact that $\mathbb{G}_{\alpha}[\Gamma_{s}]$ blows up everywhere as
$s\to 0^{+}$ [see Lemma 3.1], that is, there is no solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u=0&{\rm in}\quad
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{---\ }u=\delta_{0}&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}).\end{array}$ (1.15)

Therefore, we have to resort a barrier function, that is the minimal classical
solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=\Gamma_{0}&{\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}+u^{p}}u=0&{\rm in}\ \
B_{1}^{c}(e_{N})\setminus\\{0\\}.\end{array}$

In order to control the limit of $\\{u_{s}\\}$ near $\partial B_{1}$,
especially near the origin, some typical truncated functions have to be
constructed carefully. The second difficulty comes from the proof of the
uniqueness. We proceed by contradiction, assuming that there is two solutions
and we will show their difference could be improved the test function from
$C^{\infty}_{0}$ into $\mathbb{X}_{\alpha,B_{1}(e_{N})}$, this enables us to
use Kato’s inequality [9, Proposition 2.4] and to conclude.

Section 4 is devoted to blow-up case. The difficulty is to obtain the blow-up
everywhere in $B_{1}(e_{N})$ just from a lower bounds of $u_{s}$, see Lemma
4.1. To overcome it, we combine the symmetric of the domain and resort the
symmetry result of $u_{s}$ and then one point blowing up leads to blowing up
every where in $B_{1}(e_{N})$.

Finally, we analyse decay approximation of the weak solution for problem (1.1)
when $\alpha\to 1^{-}$. For $p\geq\frac{N+2}{N-2}$, the first challenge is to
construct a sequence upper bounds that converges to zero. To this end, we have
to study $\lim_{\alpha\to 1^{-}}(-\Delta)^{\alpha}\Phi_{\sigma}$, where
$\Phi_{\sigma}(x)=|x|^{-\sigma}$ and then use proper parameters to construct
the bounds. For $\frac{N+1}{N-1}\leq p\leq\frac{N+2}{N-2}$, $\Phi_{\sigma}$
could not be used to construct properly the upper bounds, and then we use some
argument of [15].

##  2 Dirac mass concentrated at $\\{-se_{N}\\}$ with $s\in(0,1)$

The purpose of this section is to introduce some preliminaries. First let us
state an important Comparison Principle.

###### Theorem 2.1

[6, Theorem 2.3] Let $u$ and $v$ be super-solution and sub-solution,
respectively, of

$(-\Delta)^{\alpha}u+h(u)=f\ \ {\rm in}\ \ \mathcal{O},$

where $\mathcal{O}$ is an open, bounded and connected domain of class $C^{2}$,
the function $f:\mathcal{O}\to\mathbb{R}$ is continuous and
$h:\mathbb{R}\to\mathbb{R}$ is increasing.

Suppose that $v(x)\leq u(x),\ \forall x\in\mathcal{O}^{c}$, $u$ and $v$ are
continuous in $\bar{\mathcal{O}}$. Then

$u(x)\geq v(x),\quad\forall x\in\mathcal{O}.$

Now we investigate the weak solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&{\rm in}\quad
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=\delta_{-se_{N}}&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}),\end{array}$ (2.1)

where $s\in(0,1)$. To this end, we construct a sequence of $C^{2}$ functions
to approximate the Dirac measure. Let $g_{0}:\mathbb{R}^{N}\to[0,1]$ be a
radially symmetric decreasing $C^{2}$ function with the support in
$\overline{B_{\frac{1}{2}}(0)}$ such that $\int_{\mathbb{R}^{N}}g_{0}(x)dx=1$.
For any $n\in\mathbb{N}$ and $s\in(0,1)$, we denote

$g_{n}(x)=n^{N}g_{0}(n(x+se_{N})),\qquad\forall x\in\mathbb{R}^{N}.$

Then we certainly have that

$g_{n}\rightharpoonup\delta_{-se_{N}}\ \ {\rm as}\ \ n\to+\infty,$

in the distribution sense and for any $s>0$, there exists $N_{s}>0$ such that
for any $n\geq N_{s}$,

${\rm supp}(g_{n})\subset\overline{B_{\frac{s}{2}}(-se_{N})}.$

In order to investigate the solution of (2.1), we consider the approximating
solution $w_{n}$ of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&{\rm in}\quad
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=g_{n}&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}).\end{array}$ (2.2)

###### Lemma 2.1

Assume that $p>0$ and $\\{g_{n}\\}$ is a sequence of $C^{2}$ functions
converging to $\delta_{-se_{N}}$ with supports in
$\overline{B_{\frac{s}{2}}(-se_{N})}$. Denote that

$\tilde{g}_{n}(x):=c_{N,\alpha}\int_{\mathbb{R}^{N}}\frac{g_{n}(y)}{|x-y|^{N+2\alpha}}dy,\quad\forall
x\in B_{1}(e_{N}).$ (2.3)

Then problem (2.2) admits a unique solution $w_{n}$ such that

$0<w_{n}\leq\mathbb{G}_{\alpha}[\tilde{g}_{n}]\quad{\rm in}\quad
B_{1}(e_{N}).$

Moreover, the function $\tilde{w}_{n}:=w_{n}\chi_{B_{1}(e_{N})}$ is the unique
solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=\tilde{g}_{n}&{\rm in}\quad
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=0&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}).\end{array}$ (2.4)

Proof. The existence and uniqueness of solution to problem (2.2) refers to [6,
Theorem 2.5]. For $n\geq N_{s}$, we have that supp$(g_{n})\subset
B_{\frac{s}{2}}(-se_{N})$ and then $\tilde{g}_{n}\in
C^{1}(\overline{B_{1}(e_{N})})$ and

$\tilde{w}_{n}=w_{n}-g_{n}\quad{\rm in}\quad\mathbb{R}^{N}.$

By the definition of fractional Laplacian, it implies that

$\displaystyle(-\Delta)^{\alpha}\tilde{w}_{n}(x)+\tilde{w}_{n}(x)^{p}$
$\displaystyle=$
$\displaystyle(-\Delta)^{\alpha}w_{n}(x)-(-\Delta)^{\alpha}g_{n}(x)+w_{n}(x)^{p}$
$\displaystyle=$ $\displaystyle
c_{N,\alpha}\int_{\mathbb{R}^{N}}\frac{g_{n}(z)}{|z-x|^{N+2\alpha}}dz=\tilde{g}_{n}(x).$

Therefore, $\tilde{w}_{n}$ is a classical solution of (2.4) and

$\tilde{w}_{n}\leq\mathbb{G_{\alpha}}[\tilde{g}_{n}]\quad{\rm in}\quad
B_{1}(e_{N}),$

which implies that

$w_{n}\leq\mathbb{G}_{\alpha}[\tilde{g}_{n}]\quad{\rm in}\quad B_{1}(e_{N}).$

The proof ends. $\Box$

We remark that $\tilde{w}_{n}$ is the classical solution of (2.4), then by
Lemma 2.1 and Lemma 2.2 in [9], we have that

$\int_{B_{1}(e_{N})}[w_{n}(-\Delta)^{\alpha}\xi+w_{n}^{p}\xi]dx=\int_{B_{1}(e_{N})}\xi\tilde{g}_{n}dx,\qquad\forall\xi\in
C_{0}^{\infty}(B_{1}(e_{N})).$ (2.5)

Here (2.5) holds even for $\xi\in\mathbb{X}_{\alpha,B_{1}(e_{N})}$.

###### Lemma 2.2

Let $\\{\tilde{g}_{n}\\}$ be defined in (2.3), then $\tilde{g}_{n}$ converges
to $\Gamma_{s}$ uniformly in $B_{1}(e_{N})$ and in $C^{\theta}(B_{1}(e_{N}))$
for $\theta\in(0,1)$, where

$\Gamma_{s}(x)=\frac{c_{N,\alpha}}{|x+se_{N}|^{N+2\alpha}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{-se_{N}\\}.$ (2.6)

Proof. It is obvious that $\tilde{g}_{n}$ converges to $\Gamma_{s}$ every
point in $\overline{B_{1}(e_{N})}$. For $x,y\in B_{1}(e_{N})$ and any
$n\in\mathbb{N}$, we have that

$\displaystyle|\tilde{g}_{n}(x)-\tilde{g}_{n}(y)|$ $\displaystyle=$
$\displaystyle
c_{N,\alpha}|\int_{B_{\frac{s}{2}}(-se_{N})}[\frac{1}{|x-z|^{N+2\alpha}}-\frac{1}{|y-z|^{N+2\alpha}}]g_{n}(z)dz|$
$\displaystyle\leq$ $\displaystyle
c_{N,\alpha}\int_{B_{\frac{s}{2}}(-se_{N})}\frac{||x-z|^{N+2\alpha}-|y-z|^{N+2\alpha}|}{|x-z|^{N+2\alpha}|y-z|^{N+2\alpha}}g_{n}(z)dz$
$\displaystyle\leq$ $\displaystyle
c_{N,\alpha}(N+2\alpha)|x-y|\int_{B_{\frac{s}{2}}(-se_{N})}\frac{|x-z|^{N+2\alpha-1}+|y-z|^{N+2\alpha-1}}{|x-z|^{N+2\alpha}|y-z|^{N+2\alpha}}g_{n}(z)dz$
$\displaystyle\leq$ $\displaystyle
c_{3}|x-y|\int_{B_{\frac{s}{2}}(-se_{N})}g_{n}(z)dz$ $\displaystyle=$
$\displaystyle c_{3}|x-y|,$

where $c_{3}>0$ independent of $n$. So $\\{\tilde{g}_{n}\\}_{n}$ is uniformly
bounded in $C^{0,1}(B_{1}(e_{N}))$. Combining the converging

$\tilde{g}_{n}\to\Gamma_{s}\ {\rm every\ point\ in}\ \overline{B_{1}(e_{N})}.$

We conclude that $\tilde{g}_{n}$ converges to $\Gamma_{s}$ uniformly in
$B_{1}(e_{N})$ and in $C^{\theta}(B_{1}(e_{N}))$ for $\theta\in(0,1)$. $\Box$

###### Proposition 2.1

Assume that $p>0$, $s\in(0,1)$ and $\Gamma_{s}$ is given by (2.6). Then
problem (2.1) admits a unique weak solution $u_{s}$ such that

$0\leq u_{s}(x)\leq\mathbb{G}_{\alpha}[\Gamma_{s}],\qquad x\in B_{1}(e_{N}).$
(2.7)

Moreover, $\tilde{u}_{s}:=u_{s}\chi_{B_{1}(e_{N})}$ is the unique classical
solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=\Gamma_{s}&{\rm in}\quad
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-----\ }u=0&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}).\end{array}$ (2.8)

Proof.  Existence.  It infers by Lemma 2.1 that the solution $w_{n}$ of (2.2)
satisfies that

$0<w_{n}\leq\mathbb{G}_{\alpha}[\tilde{g}_{n}]\quad{\rm in}\quad
B_{1}(e_{N}).$ (2.9)

By Lemma 2.2 we have that $\tilde{g}_{n}$ converges to $\Gamma_{s}$ uniformly
in $\overline{B_{1}(e_{N})}$ and in $C^{\theta}(B_{1}(e_{N}))$ with
$\theta\in(0,1)$. Therefore, there exists some constant $c_{4}>0$ independent
of $n$ such that

$\mathbb{G}_{\alpha}[\tilde{g}_{n}](x)\leq\frac{c_{4}c_{N,\alpha}}{|x+se_{N}|^{N+2\alpha}}\leq
c_{4}c_{N,\alpha}s^{-N-2\alpha},\qquad\forall x\in B_{1}(e_{N}).$

Thus,

$\|w_{n}\|_{L^{\infty}(B_{1}(e_{N}))}\leq
c_{4}c_{N,\alpha}s^{-N-2\alpha},\qquad\|w_{n}\|_{L^{1}(B_{1}(e_{N}))}\leq
c_{4}c_{N,\alpha}s^{-N-2\alpha}|B_{1}(e_{N})|.$

By [28, Theorem 1.2], we have that

$\begin{array}[]{lll}\|\frac{w_{n}}{\rho^{\alpha}}\|_{C^{\alpha}(\overline{B_{1}(e_{N})})}\leq
c_{5}[\|w_{n}^{p}\|_{L^{\infty}(B_{1}(e_{N}))}+\|\Gamma_{s}\|_{L^{\infty}(B_{1}(e_{N}))}]\\\\[8.53581pt]
\phantom{\|\frac{w_{n}}{\rho^{\alpha}}\|_{C^{\alpha}(\bar{B}_{1}(e_{N})}}\leq
c_{6}[s^{-N-2\alpha}+s^{-(N+2\alpha)p}]\end{array}$ (2.10)

for some $c_{5},c_{6}>0$.

In order to see the inner regularity, we denote $\mathcal{O}_{i}$ the open
sets with $i=1,2,3$ such that

$\mathcal{O}_{1}\subset\bar{\mathcal{O}}_{1}\subset\mathcal{O}_{2}\subset\bar{\mathcal{O}}_{2}\subset\mathcal{O}_{3}\subset\bar{\mathcal{O}}_{3}\subset
B_{1}(e_{N}).$

By [11, Lemma 3.1], for $\beta\in(0,\alpha)$, there exists $c_{7},c_{8}>0$
independent of $n$ such that

$\begin{array}[]{lll}\|w_{n}\|_{C^{\beta}(\mathcal{O}_{2})}\leq
c_{7}[\|w_{n}\|_{L^{1}(B_{1}(e_{N}))}+\|w_{n}^{p}\|_{L^{\infty}(\mathcal{O}_{3})}+\|w_{n}\|_{L^{\infty}(\mathcal{O}_{3})}]\\\\[8.53581pt]
\phantom{\|w_{n}\|_{C^{\beta}(\mathcal{O}_{2})}}\leq
c_{8}[s^{-N-2\alpha}+s^{-(N+2\alpha)p}].\end{array}$

It follows by [28, Corollary 2.4] that there exist $c_{9},c_{10}>0$ such that

$\begin{array}[]{lll}\|w_{n}\|_{C^{2\alpha+\beta}(\mathcal{O}_{1})}\leq{c_{9}}[\|w_{n}\|_{L^{1}(B_{1}(e_{N}))}+\|w_{n}^{p}\|_{C^{\beta}(\mathcal{O}_{2})}+\|w_{n}\|_{C^{\beta}(\mathcal{O}_{2})}]\\\\[8.53581pt]
\phantom{\|w_{n}\|_{C^{2\alpha+\beta}(\mathcal{O}_{1})}}\leq
c_{10}[s^{-N-2\alpha}+s^{-(N+2\alpha)p}].\end{array}$ (2.11)

Therefore, by the Arzela-Ascoli Theorem, there exist $u_{s}\in
C^{2\alpha+\epsilon}_{\rm loc}$ in $B_{1}(e_{N})$ for some
$\epsilon\in(0,\beta)$ and a subsequence $\\{w_{n_{k}}\\}$ such that

$w_{n_{k}}\to u_{s}\quad{\rm locally\ in}\quad C^{2\alpha+\epsilon}\quad{\rm
as}\quad n_{k}\to\infty.$ (2.12)

Passing the limit of (2.5) with $\xi\in\mathbb{X}_{\alpha,B_{1}(e_{N})}$ as
$n_{k}\to\infty$, we have that

$\int_{B_{1}(e_{N})}[u_{s}(-\Delta)^{\alpha}\xi+u_{s}^{p}\xi]dx=\int_{B_{1}(e_{N})}\xi(x)\Gamma_{s}(x).$
(2.13)

Moreover, since $w_{n}\to u_{s}$ and $\tilde{g}_{n}\to\Gamma_{s}$ uniformly in
$B_{1}(e_{N})$ as $n\to\infty$, then it infers that

$0\leq u_{s}\leq\mathbb{G}_{\alpha}[\Gamma_{s}]\quad{\rm in}\quad
B_{1}(e_{N}).$

Uniqueness. Let $v_{s}$ be a weak solution of (2.2) and then
$\varphi_{s}:=u_{s}-v_{s}$ is a weak solution to

$\begin{array}[]{lll}(-\Delta)^{\alpha}\varphi_{s}+u_{s}^{p}-v_{s}^{p}=0&{\rm
in}\quad B_{1}(e_{N}),\\\\[5.69054pt]
\phantom{(-\Delta)^{\alpha}+u_{s}^{p}-v_{s}^{p}}\varphi_{s}=0&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}).\end{array}$

By Kato’s inequality [9, Proposition 2.4],

$\int_{B_{1}(e_{N})}|\varphi_{s}|(-\Delta)^{\alpha}\xi+\int_{B_{1}(e_{N})}[u_{s}^{p}-v_{s}^{p}]{\rm
sign}(u_{s}-v_{s})\xi dx=0.$

Taking $\xi=\mathbb{G}_{\alpha}[1]$, we have that

$\int_{B_{1}(e_{N})}[u_{s}^{p}-v_{s}^{p}]{\rm sign}(u_{s}-v_{s})\xi dx\geq
0\quad{\rm and}\quad\int_{B_{1}(e_{N})}|\varphi_{s}|dx=0,$

then $\varphi_{s}=0$ a.e. in $B_{1}(e_{N})$. Then the uniqueness is proved.

Furthermore, we see that $\tilde{w}_{n}=w_{n}-g_{n}$ is the unique classical
solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u(x)+u^{p}(x)=\tilde{g}_{n}(x),&\forall
x\in B_{1}(e_{N}),\\\\[5.69054pt]
\phantom{(-\Delta)^{\alpha}+u^{p}(x)}u(x)=0,&\forall x\in
B_{1}(e_{N})^{c}\end{array}$ (2.14)

and $\tilde{w}_{n}$ converges to $\tilde{u}_{s}$ uniformly in $B_{1}(e_{N})$.
By Stability Theorem [5, Lemma 4.5] and (2.12), $u_{s}\chi_{B_{1}(e_{N})}$ is
the classical solution of (2.8). $\Box$

##  3 Proof of Theorem 1.1.

In this section, we prove Theorem 1.1 by moving the points $\\{-se_{N}\\}$ to
the origin. To this end, we need derive more properties for $u_{s}$, where
$u_{s}$ is the unique weak solution of (2.1).

###### Lemma 3.1

Let $p>0$, $s\in(0,1)$ and $u_{s}$ be the unique weak solution of (2.1). Then
the mapping $s\mapsto u_{s}$ is decreasing, that is,

$u_{s_{1}}\geq u_{s_{2}}\quad{\rm if}\ s_{1}\leq s_{2}.$

Proof. By Proposition 2.1, $\tilde{u}_{s}:=u_{s}\chi_{B_{1}(e_{N})}$ is the
unique classical solution of (2.8).

We claim that the mapping: $s\mapsto\Gamma_{s}$ is decreasing. For $x\in
B_{1}(e_{N})$ and $s_{1}\leq s_{2}$, we observe that
$|x+s_{1}e_{N}|\leq|x+s_{2}e_{N}|$, then
$\Gamma_{s_{1}}(x)\geq\Gamma_{s_{2}}(x).$ The claim is proved.

Therefore, for $s_{1}\leq s_{2}$, $u_{s_{1}}$ and $u_{s_{2}}$ are super
solution and solution of (2.8) replaced $\Gamma_{s}$ by $\Gamma_{s_{2}}$, then
it infers by the Comparison Principle that $u_{s_{1}}\geq u_{s_{2}}$ in
$B_{1}(e_{N})$. $\Box$

###### Lemma 3.2

Let $s\in(0,1)$ and denote

$\mathbb{G}_{\alpha}[\Gamma_{s}](x)=\int_{B_{1}(e_{N})}G_{\alpha}(x,y)\Gamma_{s}(y)dy,\quad\forall
x\in B_{1}(e_{N}).$

Then

$\lim_{s\to 0^{+}}\mathbb{G}_{\alpha}[\Gamma_{s}](x)=+\infty,\qquad\forall
x\in B_{1}(e_{N}).$

Proof. Using [8, Theorem 1.2], it follows that

$G_{\alpha}(x,y)\geq\min\left\\{\frac{c_{11}}{|x-y|^{N-2\alpha}},\frac{c_{11}\rho^{\alpha}(x)\rho^{\alpha}(y)}{|x-y|^{N}}\right\\},\quad
x,y\in B_{1}(e_{N}),$

where $c_{11}>0$ dependent of ${N,\alpha}$. Now for $x\in B_{1}(e_{N})$ and
$y\in B_{1}(e_{N})\cap B_{\frac{|x|}{4}}(0)$, we have that

$G_{\alpha}(x,y)\geq\frac{c_{11}\rho^{\alpha}(x)\rho^{\alpha}(y)}{|x-y|^{N}}$

and

$\displaystyle\mathbb{G}_{\alpha}[\Gamma_{s}](x)$ $\displaystyle\geq$
$\displaystyle\int_{B_{1}(e_{N})}\min\left\\{\frac{c_{11}}{|x-y|^{N-2\alpha}},\frac{c_{11}\rho^{\alpha}(x)\rho^{\alpha}(y)}{|x-y|^{N}}\right\\}\frac{c_{N,\alpha}}{|y+se_{N}|^{N+2\alpha}}dy$
$\displaystyle\geq$ $\displaystyle\int_{B_{1}(e_{N})\cap
B_{\frac{|x|}{4}}(0)}\frac{c_{11}\rho^{\alpha}(x)\rho^{\alpha}(y)}{|x-y|^{N}}\frac{c_{N,\alpha}}{|y+se_{N}|^{N+2\alpha}}dy$
$\displaystyle\geq$
$\displaystyle\frac{4}{5}c_{11}\rho^{\alpha}(x)|x|^{-N}\int_{B_{1}(e_{N})\cap
B_{\frac{|x|}{4}}(0)}\frac{c_{N,\alpha}\rho^{\alpha}(y)}{|y+se_{N}|^{N+2\alpha}}dy$
$\displaystyle\to$ $\displaystyle+\infty\qquad{\rm as}\quad s\to 0^{+}.$

The proof ends. $\Box$

From Lemma 3.2, it is informed that the limit of
$\mathbb{G}_{\alpha}[\Gamma_{s}]$ as $s\to 0^{+}$ can’t be used as a barrier
function to control the sequence $\\{u_{s}\\}$. So we have to find new upper
bound for sequence $\\{u_{s}\\}$.

###### Proposition 3.1

Let $\Gamma_{0}$ be defined in (1.10) and

$p>1+\frac{2\alpha}{N},$ (3.1)

then problem

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=\Gamma_{0}&{\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}+u^{p}}u=0&{\rm in}\ \
B_{1}^{c}(e_{N})\setminus\\{0\\}.\end{array}$ (3.2)

admits a minimum positive solution $u_{0}$, that is, $u_{0}\leq u$ for any
nonnegative solution $u$ of (3.2).

Moreover,

$\lim_{x\in B_{1}(e_{N}),\ x\to\partial
B_{1}(e_{N})\setminus\\{0\\}}u_{0}(x)=0$ (3.3)

and

$\frac{1}{c_{12}}t^{-\frac{N+2\alpha}{p}}\leq u_{0}(te_{N})\leq
c_{12}t^{-\frac{N+2\alpha}{p}},\quad t\in(0,1),$ (3.4)

where $c_{12}>1$ is independent of $\alpha$.

Proof.  The existence of solution to (3.2). It implies by Lemma 3.2 that the
mapping $s\mapsto u_{s}$ is decreasing in $B_{1}(e_{N})$, where
$u_{s}\chi_{B_{1}(e_{N})}$ is the solution of (2.8). So what we have to do is
just to find a super solution $U$ of (3.2) such that $u_{0}\leq U$ in
$B_{1}(e_{N})$. To this end, we consider the radial function

$\Phi_{\sigma}(x)=\frac{1}{|x|^{\sigma}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\},$ (3.5)

where $\sigma\in[0,\ N)$. By scaling property of $\Phi_{\sigma}$, (also see
[13]) we know that

$(-\Delta)^{\alpha}\Phi_{\sigma}(x)=\frac{c(\sigma,\alpha)}{|x|^{\sigma+2\alpha}},$
(3.6)

where $c(\sigma,\alpha)\in\mathbb{R}$. Now we choose
$\sigma=\sigma_{0}=\frac{N+2\alpha}{p}$, then $\sigma_{0}\in(0,\ N)$ if
$p>1+\frac{2\alpha}{N}$. Therefore, there exist some $k>1$ dependent of
$|c(\sigma_{0},\alpha)|$ and $c_{N,\alpha}$ but independent of $n$ such that

$U(x)=k\Phi_{\sigma_{0}}(x)$ (3.7)

is a super solution of (3.2). Thus, $U\in L^{1}(B_{1}(e_{N}))$ and by Theorem
2.1, we have that

$0\leq u_{s}\leq U\ \ {\rm for\ any}\ \ s\in(0,1).$ (3.8)

For any $x\in\mathbb{R}^{N}\setminus\\{0\\}$, $u_{0}(x):=\lim_{s\to
0^{+}}u_{s}(x)\leq U(x)<+\infty$. Following the same argument of Proposition
2.1, we can prove that $u_{0}$ is a classical solution of (3.2). Furthermore,
$u_{0}$ is the minimum solution of (3.2).

Proof of (3.3). Let $\bar{x}\in\partial B_{1}(e_{N})\setminus\\{0\\}$,
$K_{1}=\partial B_{1}(e_{N})\cap\overline{B_{|\bar{x}|/8}(\bar{x})}$ and
$K_{2}=\partial B_{1}(e_{N})\cap B_{|\bar{x}|/2}^{c}(\bar{x})$. Let
$\mathcal{O}$ be an open and $C^{2}$ set such that

$B_{1}(e_{N})\cap B_{|\bar{x}|/4}(\bar{x})\subset\mathcal{O}\subset
B_{1}(e_{N})\cap B_{|\bar{x}|/2}(\bar{x}).$

Then we see that

$K_{1}\subset\partial\mathcal{O}\quad{\rm and}\quad\partial\mathcal{O}\cap
K_{2}=\emptyset.$

We would like to find a super solution of (3.2) in $\mathcal{O}$ with
vanishing boundary value in $K_{1}$ for any $n$. Denote

$V_{\lambda}=U\eta+\lambda V_{\mathcal{O}},$

where $U$ is defined (3.7), $\eta$ is a $C^{2}$ function such that

$\eta(x)=\left\\{\begin{array}[]{lll}1&{\rm if}\quad x\in
B^{c}_{|\bar{x}|/4}(\bar{x}),\\\\[5.69054pt] 0&{\rm if}\quad x\in
B_{|\bar{x}|/8}(\bar{x}).\end{array}\right.$

and $V_{\mathcal{O}}$ is the solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u=1&{\rm
in}\quad\mathcal{O},\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}}u=0&{\rm
in}\quad\mathcal{O}^{c}.\end{array}$

Since $U(1-\eta)$ is $C^{2}$ in $\mathbb{R}^{N}$, then there exists $c_{13}>0$
such that

$\displaystyle|(-\Delta)^{\alpha}U(1-\eta)|\leq c_{13}\quad{\rm
in}\quad\bar{\mathcal{O}}.$

Thus, there exists $c_{14}>0$ such that $|(-\Delta)^{\alpha}U\eta|\leq
c_{14}\quad{\rm in}\quad\bar{\mathcal{O}}.$ Choosing $\lambda_{0}>0$ suitable,
we have that for $\lambda\geq\lambda_{0}$,

$(-\Delta)^{\alpha}V_{\lambda}+V_{\lambda}^{p}\geq\frac{c_{N,\alpha}}{|x|^{N+2\alpha}}\
\ {\rm in}\ \ \mathcal{O}.$

Moreover, since $\eta=1$ in $\mathbb{R}^{N}\setminus\mathcal{O}$, then
$V_{\lambda}\geq U\geq u_{s}$ in $\mathbb{R}^{N}\setminus\mathcal{O}$. By the
Comparison Principle, we have that for any $s\in(0,1)$,

$u_{s}\leq V_{\lambda}\quad{\rm in}\quad B_{1}(e_{N}).$ (3.9)

which implies (3.3).

Proof of (3.4). For $t\in(0,1)$, denote that

$V_{t}(x)=c_{N,\alpha}^{\frac{1}{p}}t^{-\frac{N+2\alpha}{p}}V_{B}\left(\frac{x-te_{N}}{t}\right),\quad
x\in\mathbb{R}^{N},$

where $V_{B}$ is the solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u=1&{\rm in}\quad
B_{1}(0),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}}u=0&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(0).\end{array}$ (3.10)

It deduces by (3.1) that $x\in B_{\frac{t}{4}}(te_{N})$,

$(-\Delta)^{\alpha}V_{t}(x)=c_{N,\alpha}^{\frac{1}{p}}t^{-\frac{N+2\alpha}{p}-2\alpha}\leq
c_{N,\alpha}^{\frac{1}{p}}t^{-N-2\alpha}$

and

$\frac{1}{|x+te_{N}|^{N+2\alpha}}\geq\frac{c_{15}}{t^{N+2\alpha}},$

where $c_{15}>0$ independent of $t$. Then there exists some constant
$\nu\in(0,1)$ such that

$(-\Delta)^{\alpha}(\nu V_{t})+(\nu
V_{t})^{p}\leq\frac{c_{N,\alpha}}{t^{N+2\alpha}}\leq(-\Delta)^{\alpha}u_{0}+u_{0}^{p}\quad{\rm
in}\quad B_{\frac{t}{4}}(te_{N})$

and $\nu V_{t}=0\leq u_{0}$ in $B_{\frac{t}{4}}(te_{N})^{c}$. By applying the
Comparison Principle, we have that

$\nu V_{t}\leq u_{0}\quad{\rm in}\quad\mathbb{R}^{N},$

which implies

$u_{0}(x)\geq\nu
c_{N,\alpha}^{\frac{1}{p}}t^{-\frac{N+2\alpha}{p}}\min_{B_{\frac{1}{2}}(0)}V_{B}\geq
c_{16}c_{N,\alpha}^{\frac{1}{p}}|x|^{-\frac{N+2\alpha}{p}},\quad\ \forall x\in
B_{\frac{t}{8}}(te_{N}),$ (3.11)

for some constants $c_{16}>0$ independent of $t$. We complete the proof.
$\Box$

Now we are in the position to prove Theorem 1.1.

Proof of Theorem 1.1. By Proposition 3.1, there exists minimum solution
$u_{0}$ of (3.2). Using the maximum principle argument, we have that for
$s\in(0,1)$,

$u_{s}\leq u_{0}\quad{\rm in}\quad B_{1}(e_{N}).$

Since

$u_{0}\in L^{1}(B_{1}(e_{N})),\quad u_{0}^{p}\in
L^{1}(B_{1}(e_{N}),\rho^{2}dx),$

where $\rho(x)={\rm dist}(x,\partial B_{1}(e_{N}))$, then there exists
$\tilde{u}_{0}\leq u_{0}$ such that

$u_{s}\to\tilde{u}_{0}\quad{\rm in}\ \ B_{1}(e_{N})\quad{\rm and}\quad{\rm
in}\ \ L^{1}(B_{1}(e_{N}))$

and

$u_{s}^{p}\to\tilde{u}_{0}^{p}\quad{\rm in}\ \ B_{1}(e_{N})\quad{\rm
and}\quad{\rm in}\ \ L^{1}(B_{1}(e_{N}),\rho^{2}dx).$

Passing the limit in identity of (2.13) as $s\to 0^{+}$, we have that

$\int_{B_{1}(e_{N})}[\tilde{u}_{0}(-\Delta)^{\alpha}\xi+\tilde{u}_{0}^{p}\xi]dx=\int_{B_{1}(e_{N})}\xi\Gamma_{0}dx,\qquad\forall\xi\in
C_{0}^{\infty}(B_{1}(e_{N})).$ (3.12)

By the same argument, we have that $\tilde{u}_{0}$ is $C^{2}$ locally in
$B_{1}(e_{N})$ and $\tilde{u}_{0}\chi_{B_{1}(e_{N})}$ is the minimum solution
of (3.2). Thus, we have

$\tilde{u}_{0}=u_{0}\quad{\rm in}\quad B_{1}(e_{N}).$

Uniqueness. Let $v_{0}$ be a weak solution of (2.2) and then
$\varphi_{0}:=u_{0}-v_{0}$ is a weak solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}\varphi_{0}+u_{0}^{p}-v_{0}^{p}=0&{\rm
in}\quad B_{1}(e_{N}),\\\\[5.69054pt]
\phantom{(-\Delta)^{\alpha}+u_{s}^{p}-v_{s}^{p}}\varphi_{0}=0&{\rm
in}\quad\mathbb{R}^{N}\setminus B_{1}(e_{N}),\end{array}$

that is,

$\int_{B_{1}(e_{N})}\varphi_{0}(-\Delta)^{\alpha}\xi
dx+\int_{B_{1}(e_{N})}(u_{0}^{p}-v_{0}^{p})\xi dx=0,\qquad\forall\xi\in
C_{0}^{\infty}(B_{1}(e_{N})).$ (3.13)

In this definition of weak solution, we can not apply Kato’s inequality [9,
Proposition 2.4] directly due to the stronger test functions space. So we have
to improve the regularity of $\varphi_{0}$. In fact, for any
$\xi\in\mathbb{X}_{\alpha,B_{1}(e_{N})}$, let $\\{\xi_{k}\\}$ be a sequence of
nonnegative functions in $C_{0}^{\infty}(B_{1}(e_{N}))$, such that

$\xi_{k}\to\xi\quad{\rm and}\quad|(-\Delta)^{\alpha}\xi_{k}|\leq
2|(-\Delta)^{\alpha}\xi|\quad{\rm in}\quad B_{1}(e_{N}).$

Passing the limit of (3.13) with $\xi_{k}$ as $k\to\infty$, we have that
(3.13) holds for any $\xi\in\mathbb{X}_{\alpha,B_{1}(e_{N})}.$ If we choose a
sequence $\\{\tilde{\xi}_{k}\\}\subset\mathbb{X}_{\alpha,B_{1}(e_{N})}$ which
converges to ${\rm
sign}(u_{0}^{p}-v_{0}^{p})\cdot\mathbb{G}_{\alpha,B_{1}(e_{N})}[1]$ and then
there exists $c_{6}>0$ such that

$\int_{B_{1}(e_{N})}|u_{0}^{p}-v_{0}^{p}|\mathbb{G}_{\alpha,B_{1}(e_{N})}[1]dx\leq
c_{6}\|\varphi_{0}\|_{L^{1}(B_{1}(e_{N}))},$

therefore, $|u_{0}^{p}-v_{0}^{p}|\in L^{1}(B_{1}(e_{N}),\rho^{\alpha}dx)$.

Now we follows the Kato’s inequality [9, Proposition 2.4] to obtain that

$\int_{B_{1}(e_{N})}|\varphi_{0}|(-\Delta)^{\alpha}\xi+\int_{B_{1}(e_{N})}[u_{0}^{p}-v_{0}^{p}]{\rm
sign}(u_{0}-v_{0})\xi
dx=0,\qquad\forall\xi\in\mathbb{X}_{\alpha,B_{1}(e_{N})}.$

Taking $\xi=\mathbb{G}_{\alpha}[1]$, we have that

$\int_{B_{1}(e_{N})}[u_{0}^{p}-v_{0}^{p}]{\rm sign}(u_{0}-v_{0})\xi dx\geq
0\quad{\rm and}\quad\int_{B_{1}(e_{N})}|\varphi_{0}|dx=0,$

then $\varphi_{0}=0$ a.e. in $B_{1}(e_{N})$. Then the uniqueness is proved.
$\Box$

##  4  Nonexistence

In order to prove the nonexistence of weak solution to (1.1) for
$p\in(0,1+\frac{2\alpha}{N}]$, we will prove that the solution $u_{s}$ of
problem (2.1) blows up in $B_{1}(e_{N})$ as $s\to 0^{+}$. To this end, we
first define a cone by

$\mathcal{C}=\\{x\in\mathbb{R}^{N}:\ \exists t\in(0,1)\ {\rm s.t.}\
|x-te_{N}|<\frac{t}{8}\\}.$ (4.1)

We observe that $\mathcal{C}\subset B_{1}(e_{N})$.

###### Lemma 4.1

Assume that $p>0$, $s\in(0,1)$ and $u_{s}$ is the unique solution of problem
(2.1). Then there exists $c_{17}>0$ such that

$\lim_{s\to 0^{+}}u_{s}(x)\geq c_{17}|x|^{-N},\quad\forall x\in\mathcal{C}.$

Proof. By Proposition 2.1, we have that
$\tilde{u}_{s}:=u_{s}\chi_{B_{1}(e_{N})}$ is the unique classical solution of
(2.8). For given $t\in(0,1)$, we recall

$V_{t}(x)=t^{-N}V_{B}\left(\frac{x-te_{N}}{t}\right),$

where $V_{B}$ is the solution of (3.10).

We see that for $x\in B_{\frac{t}{4}}(te_{N})$,

$(-\Delta)^{\alpha}V_{t}(x)=t^{-N-2\alpha}$

and choose $s\leq\frac{t}{4}$, then

$|x+se_{N}|\leq|x|+s\leq\frac{3}{2}t,$

which implies that

$\frac{1}{|x+se_{N}|^{N+2\alpha}}\geq\frac{c_{18}}{t^{N+2\alpha}},$

where $c_{18}>0$ independent of $t$. Then there exists some constant $\nu>0$
such that

$(-\Delta)^{\alpha}(\nu V_{t})+(\nu
V_{t})^{p}\leq\frac{c_{N,\alpha}}{t^{N+2\alpha}}\leq(-\Delta)^{\alpha}\tilde{u}_{s}+\tilde{u}_{s}^{p}\quad{\rm
in}\quad B_{\frac{t}{4}}(te_{N})$

and $\nu V_{t}=0\leq\tilde{u}_{s}$ in $B_{\frac{t}{4}}(te_{N})^{c}$, by
applying the Comparison Principle, we have that

$V_{t}\leq\tilde{u}_{s}\quad{\rm in}\quad\mathbb{R}^{N},$

which implies

$u_{s}(x)\geq\nu t^{-N}\min_{B_{\frac{1}{2}}(0)}V_{B}\geq
c_{19}|x|^{-N},\quad\forall x\in B_{\frac{t}{8}}(te_{N}),$ (4.2)

for some constants $c_{19}>0$ independent of $t$ and $s$. Combining the
increasing monotonicity, we have that

$\lim_{s\to 0^{+}}u_{s}(x)\geq c_{19}|x|^{-N},\quad\forall x\in\mathcal{C}.$

###### Lemma 4.2

Assume that $p\geq 0$, $s\in(0,1)$ and
$\tilde{u}_{s}:=u_{s}\chi_{B_{1}(e_{N})}$ is the minimal solution of (2.8).
Then
$(i)$ $\tilde{u}_{s}(x^{\prime},x_{N})$ is symmetric with respect to
$x^{\prime}$ and decreasing with to $r^{\prime}:=|x^{\prime}|$ for any
$x_{N}\in(0,2)$;
$(ii)$ $\tilde{u}_{s}(x^{\prime},x_{N})$ is decreasing in $x_{N}$ for
$x_{N}\in(1,1+\sqrt{1-|x^{\prime}|^{2}})$.

Proof. By applying the same procedure of step 1 and step 2 in proof of Theorem
1.1 in [14], we have that

$u_{s}(x^{\prime},x_{N})=u_{s}(|x^{\prime}|,x_{N})\quad{\rm for}\ \
x=(x^{\prime},x_{N})\in B_{1}(e_{N})$

and $u_{s}(r,x_{N})$ is decreasing with $r=|x^{\prime}|$. Using the same
argument from the other side, we conclude that $u_{s}(x^{\prime},x_{N})\leq
u_{s}(x^{\prime},2-x_{N})$ for $x_{N}\in(1,2)$, that is,
$u_{s}(x^{\prime},x_{N})$ is decreasing with $x_{N}\in(1,2)$.

$\Box$

Proof of Theorem 1.2. For $p=0$, the solution $u_{s}$ of (2.1) satisfies

$u_{s}=\mathbb{G}_{\alpha}[\Gamma_{s}]-\mathbb{G}_{\alpha}[1].$

By Lemma 3.2, we know that $\mathbb{G}_{\alpha}[\cdot]$ blows every where in
$B_{1}(e_{N})$, but $\mathbb{G}_{\alpha}[1]$ is bounded uniformly.

For $0<p\leq 1+\frac{2\alpha}{N}$, our proof is divided into two steps.

Step 1. we prove that

$\lim_{s\to 0^{+}}u_{s}(x)=+\infty,\quad
x\in\mathcal{B}:=\\{(x^{\prime},x_{N})\in B_{1}(e_{N}):\ x_{N}>1\\}.$

If there is $\bar{x}=(\bar{x}_{0},t_{0})\in\mathcal{B}$ such that $\lim_{s\to
0^{+}}u_{s}(\bar{x})<+\infty$, then by Lemma 4.2, we have that

$\lim_{s\to 0^{+}}u_{s}(x)\leq\lim_{s\to 0^{+}}u_{s}(\bar{x})\quad{\rm
in}\quad B_{+}:=\\{x=(x^{\prime},x_{N})\in B_{1}(e_{N}):\
|x^{\prime}|>|\bar{x}^{\prime}|,\ x_{N}>t_{0}\\}.$

Choose a nonnegative function $\xi_{0}\in C_{0}^{\infty}$ with support in
$B_{+}$, then for any $x\in B_{1}(e_{N})\setminus B_{+}$,

$\displaystyle(-\Delta)^{\alpha}\xi_{0}(x)$ $\displaystyle=$ $\displaystyle-
c_{N,\alpha}{\rm
P.V.}\int_{\mathbb{R}^{N}}\frac{\xi_{0}(y)-\xi_{0}(x)}{|x-y|^{N+2\alpha}}dy$
$\displaystyle=$ $\displaystyle-
c_{N,\alpha}\int_{B_{+}}\frac{\xi_{0}(y)}{|x-y|^{N+2\alpha}}dy$
$\displaystyle\leq$ $\displaystyle-
c_{N,\alpha}|x-2e_{N}|^{-N-2\alpha}\int_{B_{+}}\xi_{0}(y)dy,$

therefore, there exists $c_{20}>0$ such that

$(-\Delta)^{\alpha}\xi_{0}(x)\leq-c_{20},\quad\forall x\in
B_{1}(e_{N})\setminus B_{+}.$

Moreover, there exists $c_{21}>0$ such that

$|(-\Delta)^{\alpha}\xi_{0}(x)|\leq c_{21},\quad\forall x\in B_{+}.$

Since $\mathcal{C}\subset(B_{1}(e_{N})\setminus B_{+})$, then by Lemma 4.1 we
have that

$\displaystyle\int_{B_{1}(e_{N})}u_{s}(-\Delta)^{\alpha}\xi_{0}dx$
$\displaystyle\leq$
$\displaystyle\int_{\mathcal{C}}u_{s}(-\Delta)^{\alpha}\xi_{0}dx+\int_{B_{1}(e_{N})\setminus(B_{+}\cup\mathcal{C})}u_{s}(-\Delta)^{\alpha}\xi_{0}dx$
$\displaystyle+\int_{B_{+}}u_{s}|(-\Delta)^{\alpha}\xi_{0}|dx$
$\displaystyle\leq$ $\displaystyle-
c_{20}\int_{\mathcal{C}}u_{s}dx+c_{21}|B_{+}|$ $\displaystyle\to$
$\displaystyle-\infty\quad{\rm as}\quad s\to 0^{+}.$

While

$\displaystyle|\int_{B_{1}(e_{N})}u_{s}^{p}\xi_{0}dx|=\int_{B_{+}}u_{s}^{p}\xi_{0}dx\leq
u_{s}^{p}(\bar{x})\max_{\mathbb{R}^{N}}\xi_{0}|B_{+}|$

and

$\displaystyle\int_{B_{1}(e_{N})}\frac{\xi_{0}(x)}{|x+se_{N}|^{N+2\alpha}}dx=\int_{B_{+}}\frac{\xi_{0}(x)}{|x+se_{N}|^{N+2\alpha}}dx\leq|B_{+}|\max_{\mathbb{R}^{N}}\xi_{0},$

where $|x+se_{N}|\geq 1$ for $x\in B_{+}$. Then taking $s>0$ small enough, we
obtain a contradiction with the identity

$\int_{B_{1}(e_{N})}[u_{s}(-\Delta)^{\alpha}\xi+u_{s}^{p}\xi]dx=c_{N,\alpha}\int_{B_{1}(e_{N})}\frac{\xi(x)}{|x+se_{N}|^{N+2\alpha}}dx,\qquad\forall\xi\in
C_{0}^{\infty}(B_{1}(e_{N})).$

Therefore,

$\lim_{s\to 0^{+}}u_{s}(x)=+\infty,\quad\forall x\in\mathcal{B}.$

Step 2: We claim that $\lim_{s\to 0^{+}}u_{s}(x)=\infty$, $x\in B_{1}(e_{N})$.
By the fact of $\lim_{s\to 0^{+}}u_{s}(x)=\infty$, $x\in\mathcal{B}$, then
letting $\tilde{x}=(\frac{1}{2},0\cdots,\frac{3}{2})$, for any $n>1$, there
exists $s_{n}>0$ such that $s_{n}\to 0$ as $n\to+\infty$ and

$u_{s_{n}}(\tilde{x})\geq n,$

then applying Lemma 4.2, we have that

$u_{s_{n}}\geq n\quad{\rm
in}\quad\mathcal{B}_{0}=\left\\{x=(x^{\prime},x_{N})\in\mathcal{B}:\
|x^{\prime}|\leq\frac{1}{2},\ 1\leq x_{N}\leq\frac{3}{2}\right\\}.$

For any $x_{0}\in B_{1}(e_{N})\setminus\mathcal{B}$, there exists ${r_{1}}>0$
such that $\bar{B}_{r_{1}}(x_{0})\subset B_{1}(e_{N})\setminus\mathcal{B}$. We
denote by $\phi_{n}$ the solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u+u^{p}=0&{\rm in}\quad
B_{r_{1}}(x_{0}),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}+u^{p}}u=0&{\rm
in}\quad B_{r_{1}}^{c}(x_{0})\setminus\mathcal{B}_{0},\\\\[5.69054pt]
\phantom{(-\Delta)^{\alpha}+u^{p}}u=n&{\rm
in}\quad\mathcal{B}_{0}.\end{array}$ (4.3)

Then by Theorem 2.1, we have that

$u_{s_{n}}\geq\phi_{n}\quad{\rm in}\quad B_{1}(e_{N}).$ (4.4)

Let $\varphi_{n}=\phi_{n}-n\chi_{\mathcal{B}_{0}},$ then
$\varphi_{n}=\phi_{n}$ in $B_{r_{1}}(x_{0})$ and

$\displaystyle(-\Delta)^{\alpha}\varphi_{n}(x)+\varphi_{n}^{p}(x)$
$\displaystyle=$
$\displaystyle(-\Delta)^{\alpha}\phi_{n}(x)-n(-\Delta)^{\alpha}\chi_{\mathcal{B}_{0}}(x)+\phi_{n}^{p}(x)$
$\displaystyle=$ $\displaystyle
n\int_{\mathcal{B}_{0}}\frac{dy}{|y-x|^{N+2\alpha}},\qquad\forall x\in
B_{r_{1}}(x_{0}),$

that is, $\varphi_{n}$ is a solution of

$\begin{array}[]{lll}\displaystyle(-\Delta)^{\alpha}u+u^{p}=n\int_{\mathcal{B}_{0}}\frac{dy}{|y-x|^{N+2\alpha}}&{\rm
in}\quad B_{r_{1}}(x_{0}),\\\\[5.69054pt]
\phantom{(-\Delta)^{\alpha}+u^{p,}}u=0&{\rm in}\quad
B_{r_{1}}^{c}(x_{0}).\end{array}$ (4.5)

By direct computation,

$\frac{1}{c_{23}}\leq\int_{\mathcal{B}_{0}}\frac{dy}{|y-x|^{N+2\alpha}}\leq
c_{23},\quad\forall x\in B_{r_{1}}(x_{0}),$

for some $c_{23}>1$.

Let $\eta_{1}$ be the solution of

$\begin{array}[]{lll}(-\Delta)^{\alpha}u=1&{\rm in}\quad
B_{r_{1}}(x_{0}),\\\\[5.69054pt] \phantom{(-\Delta)^{\alpha}}u=0&{\rm in}\quad
B^{c}_{r_{1}}(x_{0})\end{array}$

and then $(\frac{n}{2c_{23}})^{\frac{1}{p}}\max\eta_{1}\cdot\eta_{1}$ is sub
solution of (4.5) for $n$ large enough. Then it infers by Theorem 2.1 that

$\varphi_{n}\geq(\frac{n}{2c_{23}})^{\frac{1}{p}}\max\eta_{1}\cdot\eta_{1},\qquad\forall
x\in B_{r_{1}}(x_{0}),$

which implies that

$\phi_{n}\geq(\frac{n}{2c_{23}})^{\frac{1}{p}}\max\eta_{1}\cdot\eta_{1},\qquad\forall
x\in B_{r_{1}}(x_{0}).$

Then by (4.4),

$\lim_{s\to
0^{+}}u_{s_{n}}(x_{0})\geq\lim_{n\to\infty}\phi_{n}(x_{0})=\infty.$

Since $x_{0}$ is arbitrary in $B_{1}(e_{N})\setminus\mathcal{B}$, it implies
that $\lim_{n\to\infty^{+}}u_{s_{n}}(x)=\infty$ in $B_{1}(e_{N})$. $\Box$

##  5 The solutions vanishes as $\alpha\to 1^{-}$

###  5.1 The case of $p\geq\frac{N+2}{N-2}$

###### Lemma 5.1

$(i)$ Let $N\geq 2$ and $c_{N,\alpha}$ define in (1.2). Then

$\lim_{\alpha\to 1^{-}}\frac{c_{N,\alpha}}{1-\alpha}=\frac{4N}{|S^{N-1}|},$

where $|S^{N-1}|$ denotes the $(N-1)-$dimensional measure of the unit sphere
$S^{N-1}$.

$(ii)$ Let $N\geq 2$, then for any $f\in C_{0}^{\infty}(\mathbb{R}^{N})$,

$\lim_{\alpha\to 1^{-}}(-\Delta)^{\alpha}f=-\Delta f.$

Proof. The proofs of $(i)$ and $(ii)$ see Corollary 4.2 and Proposition 4.4 in
[27] respectively. $\Box$

###### Proposition 5.1

Let $\Phi_{\sigma}$ be defined by (3.5) with $\sigma\in(0,\ N)$ and

$c(\sigma,\alpha)=-\frac{1}{2}\int_{\mathbb{R}^{N}}\frac{|z+e_{N}|^{-\sigma}+|z-e_{N}|^{-\sigma}-2}{|z|^{N+2\alpha}}dz.$

Then

$(-\Delta)^{\alpha}\Phi_{\sigma}(x)=\frac{c(\sigma,\alpha)}{|x|^{\sigma+2\alpha}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\}.$ (5.1)

Moveover,

$\lim_{\alpha\to 1^{-}}c(\sigma,\alpha)=(N-2-\sigma)\sigma$

and

$\lim_{\alpha\to
1^{-}}(-\Delta)^{\alpha}\Phi_{\sigma}(x)=\frac{(N-2-\sigma)\sigma}{|x|^{\sigma+2}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\}.$ (5.2)

Proof. By direct computation, we have that

$(-\Delta)^{\alpha}\Phi_{\sigma}(x)=\frac{c(\sigma,\alpha)}{|x|^{\sigma+2\alpha}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\},$ (5.3)

where

$c(\sigma,\alpha)=-\frac{1}{2}\int_{\mathbb{R}^{N}}\frac{\frac{1}{|z+e_{N}|^{\sigma}}+\frac{1}{|z-e_{N}|^{\sigma}}-2}{|z|^{N+2\alpha}}dz.$

On the other hand, it has been proved in [4] that

$c(\sigma,\alpha)=0\quad{\rm if}\quad\sigma=N-2\alpha.$

Now for any $R>1$, let $\eta_{R}:\mathbb{R}^{N}\to[0,1]$ be a nonnegative
function such that $\eta_{R}=1$ in $B_{R}(0)\setminus B_{\frac{1}{R}}(0)$ and
$\eta_{R}=0$ in $B_{2R}^{c}(0)\cup B_{\frac{1}{2R}}(0)$, then
$\eta_{R}\Phi_{\sigma}\in C^{\infty}_{0}(\mathbb{R}^{N})$. By Lemma 5.1
$(ii)$, we have

$\lim_{\alpha\to
1^{-}}(-\Delta)^{\alpha}(\eta_{R}\Phi_{\sigma})(e_{N})=(N-2-\sigma)\sigma$

and

$(-\Delta)^{\alpha}(\eta_{R}\Phi_{\sigma})(e_{N})=(-\Delta)^{\alpha}\Phi_{\sigma}(e_{N})-c_{N,\alpha}\int_{\mathbb{R}^{N}}\frac{(1-\eta_{R}(z))|z|^{-\sigma}}{|z-e_{N}|^{N+2\alpha}}dz.$
(5.4)

A straightforward computation implies that

$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle\int_{\mathbb{R}^{N}}\frac{(1-\eta_{R}(z))|z|^{-\sigma}}{|z-e_{N}|^{N+2\alpha}}dz$
$\displaystyle\leq$ $\displaystyle\int_{B_{\frac{1}{R}}(0)\cup
B_{R}^{c}(0)}\frac{|z|^{-\sigma}}{|z-e_{N}|^{N+2\alpha}}dz$
$\displaystyle\leq$ $\displaystyle c_{24}\int_{B_{\frac{1}{R}}(0)\cup
B_{R}^{c}(0)}\frac{|z|^{-\sigma}}{|z-e_{N}|^{N}}dz,$

where $c_{24}>0$ is independent of $\alpha$. Thus,

$\lim_{\alpha\to
1^{-}}c_{N,\alpha}\int_{\mathbb{R}^{N}}\frac{(1-\eta_{R}(z))|z|^{-\sigma}}{|z-e_{N}|^{N+2\alpha}}dz=0.$

Therefore, passing to the limit of (5.4) as $\alpha\to 1^{-}$, we conclude
that $\lim_{\alpha\to 1^{-}}c(\sigma,\alpha)=(N-2-\sigma)\sigma$ and combining
(5.3), (5.2) holds. $\Box$

###### Proposition 5.2

Let $p\geq\frac{N+2}{N-2}$, $\sigma_{p}=\frac{N+2}{p}$ and $u_{\alpha,p}$ be
the unique weak solution of (1.1). Then

$u_{\alpha,p}\leq(4^{1-\alpha}c_{N,\alpha})^{\frac{1}{p}}\Phi_{\sigma_{p}}\quad{\rm
in}\ B_{1}(e_{N}).$ (5.5)

Proof.  It follows by (5.1) that

$(-\Delta)^{\alpha}\Phi_{\sigma_{p}}(x)=\frac{c(\sigma_{p},\alpha)}{|x|^{\sigma_{p}+2\alpha}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\},$ (5.6)

where $c(\sigma_{p},\alpha)\to N-2-\sigma_{p}$ as $\alpha\to 1^{-}$. For
$k>1$, denote

$U_{k}=k\Phi_{\sigma_{p}}.$

When $p>\frac{N+2}{N-2}$, we have $N-2-\sigma_{p}>0$, then there exists
$\alpha_{p}\in(0,1)$ such that $c(\sigma_{p},\alpha)\geq 0$ for
$\alpha\in(\alpha_{p},1)$. When $p=\frac{N+2}{N-2}$, we have
$\sigma_{p}=N-2<N-2\alpha$. Using [13], $c(\cdot,\alpha)$ is $C^{2}$ and
convex in $[0,N)$ and

$\lim_{\sigma\to N^{-}}c(\sigma,\alpha)=-\infty$

and for $\sigma\in(0,N-2\alpha)$,

$c(\sigma,\alpha)\geq c(N-2\alpha,\alpha)=0.$

Thus, $c(\sigma_{p},\alpha)\geq 0$ by the fact of $\sigma_{p}<N-2\alpha$.
Choosing $k=(4^{1-\alpha}c_{N,\alpha})^{\frac{1}{p}}$, then for any $x\in
B_{1}(e_{N})$,

$\displaystyle(-\Delta)^{\alpha}U_{k}+U_{k}^{p}-\frac{c_{N,\alpha}}{|x|^{N+2\alpha}}$
$\displaystyle=$
$\displaystyle\frac{kc(\sigma_{p},\alpha)}{|x|^{\sigma_{p}+2\alpha}}+\frac{k^{p}}{|x|^{N+2}}-\frac{c_{N,\alpha}}{|x|^{N+2\alpha}}$
$\displaystyle\geq$
$\displaystyle(k^{p}-4^{1-\alpha}c_{N,\alpha})\frac{1}{|x|^{N+2}}=0,$

where we used that
$\frac{1}{|x|^{N+2\alpha}}\leq\frac{4^{1-\alpha}}{|x|^{N+2}}$ for $x\in
B_{1}(e_{N})$. Therefore, for $\alpha\in(\alpha_{p},1)$ we have that

$u_{\alpha,p}(x)\leq(4^{1-\alpha}c_{N,\alpha})^{\frac{1}{p}}\Phi_{\sigma_{p}}(x),\qquad\forall
x\in\ B_{1}(e_{N}).\qquad\Box$

Proof of Theorem 1.3. It follows by Proposition 5.2 that

$u_{\alpha,p}\leq(4^{1-\alpha}c_{N,\alpha})^{\frac{1}{p}}\Phi_{\sigma_{p}}\quad{\rm
in}\ B_{1}(e_{N}).$

On the other hand, we deduce from Lemma 5.1 that

$c_{N,\alpha}\leq c_{25}(1-\alpha),$

where $c_{25}>0$ is independent of $\alpha$. Therefore, we have that

$0\leq\lim_{\alpha\to 1^{-}}u_{\alpha,p}(x)\leq c_{26}\lim_{\alpha\to
1^{-}}(1-\alpha)^{\frac{1}{p}}\Phi_{\sigma_{p}}(x)=0,\quad\forall x\in
B_{1}(e_{N}),$

where $c_{26}>0$ is independent of $\alpha$. This ends the proof. $\Box$

###  5.2 The case $p\in[\frac{N+1}{N-1},\frac{N+2}{N-2})$

We know that the unique weak solution $u_{\alpha,p}$ of (1.1) satisfies the
identity

$\int_{B_{1}(e_{N})}\left[u_{\alpha,p}(x)(-\Delta)^{\alpha}\xi(x)+u_{\alpha,p}^{p}(x)\xi(x)\right]dx=\int_{B_{1}(e_{N})}\xi(x)\Gamma_{0,\alpha}(x)dx,\quad\forall\xi\in
C^{\infty}_{0}(B_{1}(e_{N})).$ (5.7)

From Lemma 5.1 part $(ii)$, it infers that for any $\xi\in
C^{\infty}_{0}(B_{1}(e_{N}))$,

$(-\Delta)^{\alpha}\xi\to-\Delta\xi\quad{\rm uniformly\ in}\quad
B_{1}(e_{N})\quad{\rm as}\quad\alpha\to 1^{-}.$

###### Lemma 5.2

Let $\frac{N+1}{N-1}\leq p<\frac{N+2}{N-2}$, $\sigma_{p}=\frac{N+2}{p}$ and
$u_{\alpha,p}$ be the unique weak solution of (1.1). Then there exists
$c_{27}>0$ independent of $\alpha$ such that

$0\leq u_{\alpha,p}\leq c_{27}\Phi_{\sigma_{p}}\quad{\rm in}\ B_{1}(e_{N}).$
(5.8)

Proof. For $\frac{N+1}{N-1}\leq p<\frac{N+2}{N-2}$, we have $\sigma_{p}<N$. It
follows by (5.1) that

$(-\Delta)^{\alpha}\Phi_{\sigma_{p}}(x)=\frac{c(\sigma_{p},\alpha)}{|x|^{\sigma_{p}+2\alpha}},\quad\forall
x\in\mathbb{R}^{N}\setminus\\{0\\},$

where $c(\sigma_{p},\alpha)\to N-2-\sigma_{p}$ as $\alpha\to 1^{-}$. There
exists $\alpha_{p}\in(0,1)$ such that

$|c(\sigma_{p},\alpha)|\leq 2|N-2-\sigma_{p}|,\quad\alpha\in(\alpha_{p},1).$

For $k>1$, denote

$U_{k}=k\Phi_{\sigma_{p}}.$

For $\alpha\in(\alpha_{p},1)$ and any $x\in B_{1}(e_{N})$, we deduce that

$\displaystyle(-\Delta)^{\alpha}U_{k}+U_{k}^{p}-\frac{c_{N,\alpha}}{|x|^{N+2\alpha}}$
$\displaystyle\geq$
$\displaystyle-\frac{2k|N-2-\sigma_{p}|}{|x|^{\sigma_{p}+2\alpha}}+\frac{k^{p}}{|x|^{N+2}}-\frac{c_{N,\alpha}}{|x|^{N+2\alpha}}$
$\displaystyle\geq$
$\displaystyle\left(k^{p}-4^{1-\alpha}c_{N,\alpha}-2^{N+3-\sigma_{p}}k|N-2-\sigma_{p}|\right)\frac{1}{|x|^{N+2}},$

which implies that there exists $k_{0}>0$ independent of $\alpha$ such that
for $\alpha\in(\alpha_{p},1)$ we have that

$u_{\alpha,p}(x)\leq k_{0}\Phi_{\sigma_{p}}(x),\qquad\forall x\in\
B_{1}(e_{N}).\qquad\Box$

We have that $\Phi_{\sigma_{p}}\in L^{1}(B_{1}(e_{N}))$ and

$\|\Phi_{\sigma_{p}}\|_{M^{p^{*}}(B_{1}(0),dx)}\leq c_{28},$

where $c_{28}>0$, $p^{*}=\frac{N}{\sigma_{p}}>1$ and $M^{p^{*}}(B_{1}(0))$ is
the Marcinkiewicz space of exponent $p^{*}$. By regularity results, we have
that

$\|u_{\alpha,p}\|_{C^{\beta}_{\rm loc}(B_{1}(0))}\leq c_{29},$

where $c_{29}>0$ is independent of $\alpha$ by (5.8). Therefore, by
compactness, we only have to prove that 0 is the only accumulation point of
the sequence $\\{u_{\alpha,p}\\}_{\alpha}$. Let $u^{*}$ in
$L^{1}(B_{1}(e_{N}))$ be accumulation point of the sequence
$\\{u_{\alpha,p}\\}_{\alpha}$ and a subsequence, still denote
$\\{u_{\alpha,p}\\}_{\alpha}$, converge to $u^{*}$.

On the one hand, $\Gamma_{0,\alpha}$ converges to 0 uniformly in any compact
subset of $B_{1}(e_{N})$. Therefore, for $\xi\in C_{0}^{\infty}(B_{1}(0))$, we
have that

$\lim_{\alpha\to 1^{-}}\int_{B_{1}(e_{N})}\xi(x)\Gamma_{0,\alpha}(x)dx=0.$

Thus, we have that $u^{*}$ satisfies the identity

$\int_{B_{1}(e_{N})}\left[u^{*}(x)(-\Delta)\xi(x)+(u^{*})^{p}(x)\xi(x)\right]dx=0,\quad\forall\xi\in
C^{\infty}_{0}(B_{1}(e_{N})),$ (5.9)

On the other hand, by regularity result, we have that $u^{*}\in
C^{2}(B_{1}(e_{N}))\cap L^{1}(B_{1}(e_{N}))$ continuous up to boundary
$\partial B_{1}(e_{N})\setminus\\{0\\}$, therefore, $u^{*}$ is a nonnegative
classical solution of

$\begin{array}[]{lll}-\Delta u+u^{p}=0&{\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-\Delta+u^{p}}u=0&{\rm on}\ \ \partial
B_{1}(e_{N})\setminus\\{0\\}.\end{array}$ (5.10)

Then from Theorem 3.1 in [15], we have that

$u^{*}\in L^{\infty}(B_{1}(0))$

and Theorem 3.1 in [15], $u^{*}$ is a classical solution of

$\begin{array}[]{lll}-\Delta u+u^{p}=0&{\rm in}\ \
B_{1}(e_{N}),\\\\[5.69054pt] \phantom{-\Delta+u^{p}}u=0&{\rm on}\ \ \partial
B_{1}(e_{N}),\end{array}$ (5.11)

where $p\in(1,\frac{N+1}{N-1})$. By Strong Maximum Principle, we have that

$u^{*}\equiv 0\quad{\rm in}\quad B_{1}(e_{N}).$

Acknowledgements: The authors want to thank the referee for carefully reading
the paper and making suggestions that resulted in a great improvement in the
clarity of the proofs. H. Chen is supported by National Natural Science
Foundation of China, No:11401270 and the Project-sponsored by SRF for ROCS,
SEM.

## References

  * [1] D. Balenu, K. Diethelm, E. Scalas and J. Trujillo, Fractional calculus models and numerical methods, Amsterdam, Word Science Publish, Nonlinear Science and complexity.
  * [2] M. Bidaut-Véron, A. C. Ponce and L. Véron, Isolated boundary singularities of semilinear elliptic equations, arXiv:0902.0449.
  * [3] M. F. Bidaut-Véron and C. Yarur, Semilinear elliptic equations and systems with measure data: existence and a priori estimates, Advances in Differential Equations 7(3), 257-296 (2002).
  * [4] L. Caffarelli and L. Silvestre, An extension problem related to the fractional Laplacian, Commun. PDE 32, 1245-1260 (2007).
  * [5] L. Caffarelli and L. Silvestre, Regularity theory for fully nonlinear integro-differential equations, Comm. Pure Appl. Math. 62, 597-638 (2009).
  * [6] H. Chen, P. Felmer and A. Quaas, Large solution to elliptic equations involving fractional Laplacian, Ann. Inst. H. Poincaré, Analyse Non Linéaire, doi:10.1016/j.anihpc.2014.08.001.
  * [7] H. Chen and H. Hajaiaj, Existence, Non-existence, Uniqueness of solutions for semilinear elliptic equations involving measures concentrated on boundary, arXiv:1410.2672.
  * [8] Z. Chen and R. Song, Estimates on Green functions and poisson kernels for symmetric stable process, Math. Ann. 312, 465-501 (1998).
  * [9] H. Chen and L. Véron, Semilinear fractional elliptic equations involving measures, J. Differential Equations 257(5), 1457-1486 (2014).
  * [10] H. Chen and L. Véron, Semilinear fractional elliptic equations with gradient nonlinearity involving measures, J. Funct. Anal. 266(8), 5467-5492 (2014).
  * [11] H. Chen and L. Véron, Weakly and strongly singular solutions of semilinear fractional elliptic equations, Asymptotic Analysis 88, 165-184 (2014).
  * [12] E. Elshahed, A fractional calculus model in semilunar heart valvevibrations, International Mathematica symposium (2003).
  * [13] P. Felmer and A. Quaas, Fundamental solutions and Liouville type theorems for nonlinear integral operators, Advances in Mathematics 226, 2712-2738 (2011).
  * [14] P. Felmer and Y. Wang, Radial symmetry of positive solutions to equations involving the fractional laplacian, Comm. Cont. Math. 16, 1-24 (2014).
  * [15] A. Gmira and L. Véron, Boundary singularities of solutions of some nonlinear elliptic equations. Duke Math. J. 64(2), 271-324 (1991).
  * [16] Q. Guan and Z. Ma, Boundary problems for fractional Laplacians, Stochastics and Dynamics 5, 385-424 (2005).
  * [17] I. Herbst, Spectral Theory of the operator $(p^{2}+m^{2})^{1/2}-Ze^{2}/r$, Commun. Math. Phys. 53, 285-294 (1977).
  * [18] E. Lieb and H. Yau, The Chandrasekhar theory of stellar collapse as the limit of quantum mechanics, Commun. Math. Phys. 112, 147-174 (1987).
  * [19] R. Magin, Fractional calculus in bioenginering 1, 2, 3, Critical Reviews in Biomedical Engineering 32, 1-104, 105-193, 194-377 (2004).
  * [20] R. Magin, S. Boregowda, and C. Deodhar, Modeling of pulsating peripheral bioheat transf using fractional calculus and constructal theory, Journal of Design Nature 1, 18-33 (2007).
  * [21] R. Magin, X. Feng, and D. Baleanu, Fractional calculus in NMR, Magnetic Resonance Engineering 34, 16-23 (2009).
  * [22] R. Magin and M. Ovadia, Modeling the cardiac tissue electrode interface using fractional calculus, Journal of Vibration and Control 19, 1431-1442 (2009).
  * [23] M. Marcus and L. Véron, The boundary trace of positive solutions of semilinear elliptic equations: the subcritical case, Arch. Rat. Mech. Anal. 144, 201-231 (1998).
  * [24] M. Marcus and L. Véron, The boundary trace of positive solutions of semilinear elliptic equations: the supercritical case, J. Math. Pures Appl. 77, 481-524 (1998).
  * [25] M. Marcus and L. Véron, Removable singularities and boundary trace, J. Math. Pures Appl. 80, 879-900 (2000).
  * [26] M. Marcus and L. Véron, Capacitary estimates of positive solutions of semilinear elliptic equations with absorption, J. European Math. Soc. 6, 483-527 (2004).
  * [27] E. Di Nezza, G. Palatucci and E. Valdinoci, Hitchhiker’s guide to the fractional Sobolev spaces, Bull. Sci. Math. 136, 521-573 (2012).
  * [28] X. Ros-Oton and J. Serra, The Dirichlet problem for the fractional laplacian: regularity up to the boundary, J. Math. Pures Appl. 101(3), 275-302 (2014).
  * [29] V. Tarasov and G. Zaslavsky, Fractional dynamics of systems with long-range interaction. Commun. Nonlinear Sci. Numer. Sim. 11, 885-889 (2006).


Safe Reinforcement Learning for Strategic Bidding of Virtual Power Plants in Day-Ahead Markets

Ognjen Stanojev1, Lesia Mitridati2, Riccardo de Nardis di Prata1, Gabriela Hug1
1 EEH - Power Systems Laboratory, ETH Zurich, Switzerland
2 Center for Electric Power and Energy, Technical University of Denmark, Kgs. Lyngby, Denmark
Emails: {stanojev<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Research supported by NCCR Automation, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 51NF40_180545).

This paper presents a novel safe reinforcement learning algorithm for strategic bidding of Virtual Power Plants (VPPs) in day-ahead electricity markets. The proposed algorithm utilizes the Deep Deterministic Policy Gradient (DDPG) method to learn competitive bidding policies without requiring an accurate market model. Furthermore, to account for the complex internal physical constraints of VPPs, we introduce two enhancements to the DDPG method. Firstly, a projection-based safety shield that restricts the agent's actions to the feasible space defined by the non-linear power flow equations and operating constraints of distributed energy resources is derived. Secondly, a penalty for the shield activation in the reward function that incentivizes the agent to learn a safer policy is introduced. A case study based on the IEEE 13-bus network demonstrates the effectiveness of the proposed approach in enabling the agent to learn a highly competitive, safe strategic policy.

virtual power plants, strategic bidding, electricity markets, safe reinforcement learning

§ INTRODUCTION

Growing environmental concerns and advancements in communication and monitoring technologies have led to the increased deployment of Distributed Energy Resources (DERs) in power networks <cit.>, comprising renewable energy sources and prosumers.
The market integration of these units is facilitated by their large-scale aggregation under financial entities, commonly known as Virtual Power Plants (VPPs), which have the capacity for trading in wholesale electricity markets <cit.>. As a self-interested market participant, a VPP aims at maximizing its own profit generated by its market participation and the fulfillment of contractual obligations towards its internal customers. Concurrently, VPPs must ensure the safe operation of the underlying distribution network.

In electricity markets with imperfect competition, self-interested market participants, such as large-scale VPPs, may exercise market power by bidding strategically to influence the clearing prices and increase their profits <cit.>.
This bidding strategy has traditionally been modeled as a hierarchical optimization problem <cit.>, in which the strategic market participant can anticipate the impact of its bids on the market clearing problem, provided that it has accurate information on all other parameters of the market clearing process.
However, this approach has many limitations in practice, including the assumption of perfect information on other market participants' bidding strategies, leading to optimistic and unrealistic solutions of the market clearing. Furthermore, computational complexity is often an issue, even under limiting assumptions on the convexity of the market clearing problem <cit.>.

The rapid progression in Reinforcement Learning (RL) has opened up new possibilities for market participants to devise sophisticated bidding strategies and learn to maximize their profits by interacting with the environment (i.e., market) <cit.>.
One key advantage of model-free RL methods is that a strategic market participant (RL agent) does not require any prior knowledge of the market-clearing process or the other participants' strategies (environment).
However, traditional $Q$-learning-based RL algorithms struggle with large state-action spaces, as they rely on complex tables to approximate the value function for each state-action pair, leading to trade-offs between optimality, computation time, and memory allocation <cit.>. For this reason, their applicability in practice is limited.

The above-mentioned scalability issues can be addressed by employing deep RL methods like the Deep Deterministic Policy Gradient (DDPG) algorithm <cit.>, which utilizes neural networks to extend the $Q$-learning capabilities to continuous state and action spaces.
The authors in <cit.> propose deep RL methods for the economic dispatch and market participation of DERs aggregated in a VPP. The main limitation of these works is that they fail to account for the complex internal physical constraints of large-scale VPPs, such as power generation limits and power flow constraints, in order to ensure a safe operation. In safety-critical infrastructures, such as distribution networks, this constitutes a significant limitation for the deployment of RL-based methods.
Consequently, it is necessary to investigate safe RL methods that can reliably ensure the fulfillment of these physical constraints.

Recent approaches in the field of safe RL can broadly be classified into two categories, providing improvements to (i) the reward function or (ii) the exploration phase of RL algorithms <cit.>. In the first category, a popular approach is to augment the reward function to account for constraint violations using the Lagrangian relaxation method <cit.>. However, such approaches may require ad hoc tuning of the constraint violation reward and may result in unsafe decisions during the exploration phase. In the second category, the safety of the decisions is promoted by offline (batch) learning to initialize the exploration <cit.> or by the transfer of expert knowledge learned offline to guide the exploration <cit.>. Despite significant improvements, these approaches cannot provide safety guarantees and are not suitable for fully online learning.

Given the aforementioned challenges, this paper focuses on the strategic bidding problem of VPPs in day-ahead electricity markets and introduces a new safe RL algorithm in this context. The proposed approach employs the DDPG algorithm to facilitate the learning of complex policies that operate in continuous action and state spaces. To enhance safety, we combine (i) a safety shield that restricts the agent's actions to a safe space defined by the non-linear power flow equations and DER operating constraints and (ii) a penalty for shield activation in the reward function that incentivizes the agent to learn a safer policy. The projection-based safety shield acts as a protective measure that prohibits the agent from making unsafe decisions and instead selects the closest safe action, extending the works in <cit.>. The proposed approach provides flexibility in the design of the agent and the representation of its economic and physical aspects. Numerical results based on the IEEE 13-bus network demonstrate the ability of the developed agent to devise a highly competitive and safe bidding strategy.

The subsequent sections of this paper are organized as follows. Section <ref> provides an overview of VPP modeling and introduces the market clearing problem under consideration. In Section <ref>, we introduce the proposed bidding scheme, delineating the RL problem formulation and the derivation of the safety-enforcing optimization problem. Section <ref> presents the case studies conducted to validate and assess the bidding scheme's performance, while Section <ref> concludes the paper.

§ VIRTUAL POWER PLANT AND MARKET MODELS

The VPP under study is assumed to serve load and operate a multitude of heterogeneous DERs within a distribution network, namely, renewable energy sources (wind or PV), energy storage units, and conventional generation such as diesel generators. The underlying distribution network is considered to be radial and balanced, and can thus be represented by a tree graph $\mathcal{G}(\mathcal{N},\mathcal{E})$, with $\mathcal{N} \coloneqq \{0,1,\dots,n\}$ denoting the set of network nodes including the substation node $0$, and $\mathcal{E} \subseteq \mathcal{N}\times\mathcal{N}$ representing the set of $n$ network branches. A single Point of Common Coupling (PCC) of the VPP to the rest of the system at node $0$ is assumed. Loads are considered to be located at nodes $\mathcal{L}\subseteq\mathcal{N}$, Renewable Energy Sources (RESs) at nodes $\mathcal{R}\subseteq \mathcal{N}$, Battery Energy Storage Systems (BESSs) at nodes $\mathcal{B}\subseteq \mathcal{N}$, and conventional generators at nodes $\mathcal{D}\subseteq \mathcal{N}$. Without loss of generality, each node is assumed to contain a single DER type, i.e., $\mathcal{R}\cap\mathcal{B}\cap\mathcal{D}=\varnothing$.

§.§ Grid Constraints - Flow and Voltage Deviation Limits

For every bus $i\in\mathcal{N}$, let $v_i\in\R_{\geq 0}$ denote the voltage magnitude, and $p_{i}\in\R$ and $q_{i}\in\R$ represent the active and reactive nodal power injections. For each branch $(i,j)\in\mathcal{E}$, let $r_{ij}\in\R_{\geq 0}$ and $x_{ij}\in\R$ denote its respective resistance and reactance values, $P_{ij}\in\R$ and $Q_{ij}\in\R$ denote the real and reactive power flow along the branch, and $i_{ij}\in\R_{\geq 0}$ represent the corresponding branch current magnitude. The distribution grid is modelled using the DistFlow equations <cit.>, written recursively for every line $(i,j)\in\mathcal{E}$ as:

\begin{align}
P_{ij}  &= \sum_{k\in\mathcal{N}_j} P_{jk} - p_{j} + r_{ij}i_{ij}^2, \label{eq:p_distFlow} \\
Q_{ij}  &= \sum_{k\in\mathcal{N}_j} Q_{jk} - q_{j} + x_{ij}i_{ij}^2, \label{eq:q_distFlow} \\
v_i^2 - v_j^2 &= 2(r_{ij}P_{ij}+x_{ij}Q_{ij}) - (r_{ij}^2+x_{ij}^2)i_{ij}^2, \label{eq:v_distFlow}
\end{align}

with the branch currents computed as $i_{ij}^2=(P_{ij}^2+Q_{ij}^2)/v_i^2$, and $\mathcal{N}_j$ denoting the set of child nodes of bus $j$.

To respect power quality standards and preserve the health of assets within the VPP, network constraints in the form of voltage deviations and line power limits need to be satisfied:

S_ij    ≤S^max_ij,            ∀(i,j)∈ℰ,

v^min_i ≤v_i    ≤v^max_i,             ∀i∈𝒩,

where $S_{ij}^2=P_{ij}^2+Q_{ij}^2$ is the apparent power flowing through line $(i,j)$, and $\mathrm{min}$ and $\mathrm{max}$ superscripts indicate respectively the minimum and the maximum allowable values.

§.§ DERs Constraints - Individual Capability Curves

The active and reactive power outputs of the DERs within the VPP are additionally limited by their internal constraints. These constraints are typically defined for each unit $i$ as sets of allowable operating points $(p_{i},q_{i})\in\mathcal{X}_{i}\subseteq\R^2$, called capability curves. The capability curve of a renewable generator $i\in\mathcal{R}$ can be defined by the following set:

\begin{equation*}
\begin{split}
\mathcal{X}_i \coloneqq \{(p_{i},q_{i})\in\R^2: 0\leq p_{i} \leq p^\mathrm{res}_{i}, s_i\leq s^\mathrm{max}_i,  \cos{\phi}\geq a\},
\end{split}
\end{equation*}

where $s^\mathrm{max}_{i}\in\R_{\geq0}$ denotes the apparent power limit and $p^\mathrm{res}_{i}\in\R_{\geq0}$ represents the irradiation/wind speed dependent active power upper limit, while
the lower limit is zero due to the possibility of curtailment. Due to its dependency on exogenous variables, the upper limit $p^\mathrm{res}_i$ is modeled as a random variable with an appropriate probability distribution function derived from the historical data. Furthermore, power quality standards commonly require the power factor $\cos{\phi}$ to be greater than a prescribed value $a\in[0,1]$.

The BESS capability curve is assumed to allow four-quadrant operation and is defined for each unit $i\in\mathcal{B}$ by:

\begin{align*}
\mathcal{X}_i \coloneqq \{(p_{i},q_{i})\in\R^2: -&p^\mathrm{max}_{i} \leq p_{i} \leq p^\mathrm{max}_{i}, s_i\leq s^\mathrm{max}_i, \\
&\chi^\mathrm{min}_{i} \leq \chi_{i} \leq \chi^\mathrm{max}_{i}
\},
\end{align*}

with $p^\mathrm{max}_{i}\in\R_{\geq0}$ denoting the maximum charging/discharging power, and $s^\mathrm{max}_{i}\in\R_{\geq0}$ representing the apparent power limit. Storage energy capacity limits constrain the state-of-energy $\chi\in\R_{\geq0}$ between  $\chi^\mathrm{min}\in\R_{\geq0}$ and $\chi^\mathrm{max}\in\R_{\geq0}$.

Operation of each conventional generator $i\in\mathcal{D}$ is constrained by stator current and low load limits, which impose the minimum $p_i^\mathrm{min}\in\R_{\geq0}$ and the maximum $p_i^\mathrm{max}\in\R_{\geq0}$ active power output, as well as the maximum apparent power $s_i^\mathrm{max}\in\R_{\geq0}$ limits, that is,

\begin{equation*}
\mathcal{X}_i \coloneqq \{(p_{i},q_{i})\in\R^2: p^\mathrm{min}_i \leq p_{i} \leq p_i^\mathrm{max},\,\, s_i\leq s^\mathrm{max}_i\}.
\end{equation*}

Loads are assumed to be non-flexible. Therefore, we have $\mathcal{X}_i \coloneqq \{(p_{i},q_{i})\in\R^2: p_i=-p_i^\mathrm{load},q_i=-q_i^\mathrm{load}\}$ for load $i\in\mathcal{L}$, where $p_i^\mathrm{load}\in\R_{\geq0}$ and $q_i^\mathrm{load}\in\R_{\geq0}$ are the constant active and reactive power consumptions. Finally, the buses that are not hosting any generation or consumption element have zero injections, i.e., $p_i=q_i=0,\forall i\in\mathcal{N}\setminus(\mathcal{L} \cup \mathcal{R} \cup \mathcal{D} \cup \mathcal{B})$.

§.§ Day-Ahead Market Clearing Model

The day-ahead electricity market is operated through a blind auction which takes place once a day, when each market participant submits independent price-quantity bids for each hour of the following day. We distinguish between supply bids defined by $(\gamma^s_n,p^{s,\mathrm{bid}}_n)$ for each supplier $n\in\mathcal{S}$, where $\gamma^s_n$ denotes the price offer for the quantity $p^{s,\mathrm{bid}}_n$; and demand bids defined by $(\gamma^d_j,p^{d,\mathrm{bid}}_j)$ for each customer $j\in\mathcal{C}$, where $\gamma^d_j$ denotes the price customer $j$ is willing to pay for the quantity $p^{d,\mathrm{bid}}_j$. The market clearing process is modeled by the following optimization problem:

(𝒫_1)      p^s_n, p^d_jmax          ∑_j∈𝒞 γ^d_jp^d_j - ∑_n∈𝒮 γ^s_np^s_n

s.t.        ∑_j∈𝒞 p^d_j = ∑_n∈𝒮 p^s_n,

0≤p^d_j ≤p^d,bid_j, ∀j∈𝒞,

0 ≤p^s_n ≤p^s,bid_n, ∀n∈𝒮,

aiming at maximizing the social welfare (<ref>), subject to constraints on the power balance in the transmission grid (<ref>) and limits on the cleared quantities $(p^s_n,p^d_j)$ (<ref>)-(<ref>). All the market participants are paid the uniform Market Clearing Price (MCP), which is computed as the dual variable of (<ref>).

§ SAFE RL-BASED BIDDING STRATEGY

§.§ Agent-Environment Interaction

The proposed bidding process of the strategic VPP in the day-ahead market is presented in Figure <ref>. The bidding decisions are determined by an RL-based agent, designed as a single-agent actor-critic model exploiting the DDPG algorithm. The agent is in control of generating two values, namely the bidding quantity $p_\mathrm{bid}^\mathrm{vpp}\in\R$ and the bidding price $\gamma^\mathrm{bid}_\mathrm{vpp}\in\R_{\geq0}$, collected in the action vector $a \coloneqq (\gamma^\mathrm{bid}_\mathrm{vpp},p_\mathrm{bid}^\mathrm{vpp})$.

The agent interacts with the environment composed of the shield, the market clearing engine, and an Optimal Power Flow (OPF) module. The generated bid is first passed to the shield that performs a nonlinear projection $\Tilde{p}_\mathrm{bid}^\mathrm{vpp} = \proj_\mathcal{A}(p_\mathrm{bid}^\mathrm{vpp})$ on a safe set $\mathcal{A}$, which captures the feasible operating space of the VPP. The corrected bidding quantity $\Tilde{p}_\mathrm{bid}^\mathrm{vpp}$ is submitted to the market which is then cleared according to the procedure described in (<ref>). After the market clearing, an OPF algorithm is run to dispatch the units within the VPP to respect the cleared capacity. The environment state is defined by $s \coloneqq (t,p^\mathrm{load},q^\mathrm{load},p^\mathrm{res},\chi)$. Using the hour of the day (i.e., the time step $t$) as a state variable the agent attempts to learn daily price patterns in order to place better bids. To understand how much power it needs to buy, or how much power generation it has left to sell, it is necessary to pass the expected load values $(p^\mathrm{load},q^\mathrm{load})$ and renewable generation $p^\mathrm{res}$ output as state variables. Lastly, the state-of-energy vector $\chi$ indicates the VPP's storage capacity available to the agent. Note that time indexing of variables is omitted. Unless otherwise specified, each variable refers to the current time step $t$ hereafter.

The reward $r(s,a) = r_\mathrm{da} - c_\mathrm{vpp} -c_\mathrm{shd}$ consists of three different components, the market profit $r_\mathrm{da}$, the VPP internal cost $c_\mathrm{vpp}$, and the shield activation cost $c_\mathrm{shd}$. The first component is the profit that the RL agent generates by bidding in the market and is computed as $r_\mathrm{da}=P_\mathrm{disp}\cdot\mathrm{MCP}$. Since the day-ahead market is a uniform price auction, the dispatched (cleared) quantity is multiplied by the MCP. The second component $c_\mathrm{vpp}$ reflects the internal financial activities of the VPP. The internal costs are calculated by subtracting the load supply revenues from the power generation costs:

\begin{equation}
c_{\mathrm{vpp}} =  \sum_{i \in {\mathcal{R}\cup\mathcal{D}}} p_i \gamma_i - \sum_{n \in\mathcal{L}} p_n \gamma_n.
\end{equation}

Lastly, the shield-activation penalty, which quantifies how much the shield had to intervene on the original bid in order to make it feasible, is defined as
$c_{\mathrm{shd}} = \epsilon \| \Tilde{p}_\mathrm{bid}^\mathrm{vpp} - p_\mathrm{bid}^\mathrm{vpp} \|$.
The non-negative weight $\epsilon\geq0$ is a tuning parameter.

The proposed agent-environment interaction scheme.

§.§ Deep Deterministic Policy Gradient Algorithm

Deep Deterministic Policy Gradient is an actor-critic, model-free RL algorithm <cit.> suitable for continuous state and action space problems such as the considered VPP bidding problem. The actor-critic methods employ a policy function and an action-value function. The policy function acts as an actor, generating an action $a$ given a state observation $s$, whereas the action-value function acts as a critic, which takes the state-action pair $(s,a)$ as the input and outputs the $Q$-value associated with this pair. In the DDPG algorithm, two neural networks are used to approximate the action-value function and the policy function: the critic network $Q(s,a\,\lvert\, \theta^{Q})$ and the actor network $\mu(s\,\lvert\, \theta^{\mu})$, respectively, with $\theta^{Q}$ and $\theta^{\mu}$ being the parameters of the respective networks <cit.>.

Furthermore, two additional networks, the target actor network $\mu'(s\,\lvert\,\theta^{\mu'})$ and the target critic network $Q'(s,a\,\lvert\,\theta^{Q'})$ are introduced as time-delayed copies of the actor and the critic:

\begin{equation*} \label{eq:target_net_updates}
\theta^{Q'} \mapsfrom \tau \theta^{Q} + (1-\tau) \theta^{Q'}, \quad
\theta^{\mu'} \mapsfrom \tau \theta^{\mu} + (1-\tau) \theta^{\mu'},
\end{equation*}

with the tracking rate defined by the soft update coefficient $\tau \ll 1$. These additional networks are used for rendering the learning process more stable by creating the labels for training of the original networks, as will be shown in the following.

The training process is conducted by running a predefined number of episodes, each comprising a series of steps in which the agent interacts with the environment. During the training, the agent explores the environment by allowing less favorable actions with respect to the current knowledge. The exploration is achieved using Gaussian noise $\mathcal{N}(0,\sigma(m))$ with zero mean and an exponentially decreasing variance $\sigma(m)$ defined by

\begin{equation*}\label{eq:gauss_noise}
\sigma(m) = z_\mathrm{i} \exp\Big({-m \frac{\ln(\frac{z_\mathrm{i}}{z_\mathrm{f}})}{m_\mathrm{tot}}}\Big),
\end{equation*}

where $m$ is the episode number, $m_\mathrm{tot}$ is the total number of episodes used for training, $z_\mathrm{i}$ is the desired initial noise, $z_\mathrm{f}$ is the desired final noise. The noise sampled from the Gaussian process is superposed on the actor's output $a_m = \mu(s_m \,\lvert\, \theta^{\mu} )+\mathcal{N}(0,\sigma(m))$. After each interaction, the tuple $(s_m,a_m,r_m,s_{m+1})$ is stored into the experience replay memory, from which the minibatches (i.i.d. sets of samples) are selected for training of the original networks. The replay memory is a technique used to speed up and enhance learning.

Finally, we present the update rules for the actor and critic networks. The critic network is updated by minimizing the loss function, averaged over $N$ samples in the minibatch:

\begin{equation*} \label{eq:critic_loss}
L(\theta^{Q}) = \frac{1}{N} \sum_{i=1}^{N}(y_{i} - Q(s_i,a_i\,\lvert\, \theta^{Q}))^2.
\end{equation*}

The label $y_{i}$ for the $i$-th sample in the minibatch is calculated as a sum of the immediate reward received in that sample and the expected $Q$-function value of the next state $s_i'$, determined by the target actor and critic networks, i.e.,

\begin{equation*} \label{eq:critic_label}
y_{i} = r_{i} + \varepsilon Q'\left(s_{i}', \mu'(s_{i}' \,\lvert\, \theta^{\mu'}) \,\lvert\, \theta^{Q'}\right),
\end{equation*}

where $\varepsilon\geq0$ is the discount factor.
The actor network can be improved by maximizing the policy score function defined by

\begin{equation*} \label{eq:policy_score}
J(\theta^\mu)=\mathop{\mathbb{E}} \big[Q(s,a \,\lvert\, \theta^{Q}) \,\lvert\,  _{s=s_{i},a=\mu(s_{i})}],
\end{equation*}

which evaluates performance of policy $\mu(\cdot\,\lvert\,\theta^\mu)$ for each sample in the minibatch. To maximize the score function, gradient ascent is thus applied to the actor network. This involves approximating the gradient by the average value of the policy score function gradients across the minibatch:

\begin{equation*}
\nabla_{\theta^{\mu}} J \approx \frac{1}{N} \sum_{i=1}^{N} (\nabla_{a} Q(s_i,a \,\lvert\, \theta^{Q})\,\lvert\,  _{a=\mu(s_{i})} \nabla_{\theta^{\mu}} \mu(s_i \,\lvert\, \theta^{\mu} )).
\end{equation*}

§.§ Safe RL via Projection on the Feasible Space

While the described DDPG algorithm has previously been shown to converge to effective bidding strategies in an electricity market framework <cit.>, there is yet no guarantee that the bidding decisions are feasible with respect to the VPP constraints. Namely, the submitted quantity might exceed the aggregate capacity of the VPP units or cause violations in grid voltages and currents.
To alleviate this issue, the actions generated by the learned policy are projected into a safe set. The projection thus operates as a shield that prevents RL from taking unsafe decisions and adopts the closest safe action to the generated unsafe RL decision, in the $\ell_2$-norm sense.

As safety limitations, in this paper we consider the VPP constraints introduced in Sec. <ref> and Sec. <ref> that must be respected for each submitted bid. This set of constraints is denoted by $\mathcal{A}$. To perform the projection of a generated bid onto $\mathcal{A}$, the following optimization problem is solved:

\begin{equation}\label{eq:shield_proj}
\begin{aligned}
\Tilde{p}_\mathrm{bid}^\mathrm{vpp} = \argmin_{u} \, & \frac{1}{2}\| u - p_\mathrm{bid}^\mathrm{vpp} \|^2_2 \\
\textrm{s.t.} \quad & \eqref{eq:distFlow}, \eqref{eq:s_limit},\eqref{eq:v_limit},\\
& (p_i,q_i)\in\mathcal{X}_i,\quad\forall i\in\mathcal{N},\\
& P_{01} = - u,
\end{aligned}
\end{equation}

where the final equality constraint is a boundary condition expressing that $\Tilde{p}_\mathrm{bid}^\mathrm{vpp}$ is to be dispatched at the PCC.

It should be noted that this procedure does not provide a guarantee that the projected action is the most optimal within the set of feasible actions $\mathcal{A}$. However, through the use of a reward system that combines the profit generated by the feasible action $\tilde{p}_\mathrm{bid}^\mathrm{vpp}$ with a penalty for actions outside of $\mathcal{A}$, the RL agent can learn to restrict its actions to within the safe set $\mathcal{A}$. Thus, as the learning process continues, the agent will eventually converge to locally or globally optimal safe actions.

§.§ Optimal Power Flow Dispatch

The OPF procedure is used to determine the operating points of all the generators inside the VPP. Additionally, the results of the OPF can indicate the Production Marginal Cost (PMC) of the VPP, i.e., the production cost of the last internal generator dispatched. The considered OPF problem is formulated as:

\begin{equation} \label{eq:opf_nnl}
\begin{aligned}
({\mathcal{P}_2})\quad\min_{p_i,q_i} \quad & \sum_{i \in {\mathcal{R}\cup\mathcal{D}}} p_i \gamma_i-\sum_{n \in\mathcal{L} } p_n \gamma_n\\
\textrm{s.t.} \quad & \eqref{eq:distFlow}, \eqref{eq:s_limit},\eqref{eq:v_limit},\\
& (p_i,q_i)\in\mathcal{X}_i,\quad\forall i\in\mathcal{N},\\
& P_{01} = - P_\mathrm{disp}.
\end{aligned}
\end{equation}

The value of the objective function is the previously defined internal VPP cost $c_\mathrm{vpp}$.
On the other hand, the production marginal cost of the VPP is equal to the cost difference of producing one more kWh of power. Note that the above problem (<ref>) is a nonlinear, nonconvex optimization problem due to the inclusion of DistFlow constraints (<ref>). To enhance the computational tractability of this problem, one may need to resort to relaxations or approximations. Nonetheless, the pursuit of such techniques lies beyond the scope of this paper.

§ NUMERICAL RESULTS

§.§ Case Study Setup

The algorithm developed in this paper is tested on a modified, single-phase version of the IEEE 13-bus distribution feeder <cit.>, shown in Fig. <ref>. A number of DERs have been introduced to transform the original distribution grid into a VPP.
In order to construct a system with realistic load data, the dataset provided in <cit.> is used, where load from 20 households has been measured and stored continuously over a period of two years. The household load power values are normalized and adjusted to have the nominal capacity equal to the specifications of the IEEE 13-bus distribution feeder. All the loads $i\in\mathcal{L}$ are paying a fixed price of $\gamma_i=\SI[per-mode=symbol]{0.5}{\EUR\per\kWh}$ to the VPP owner for the provided supply.

Single line diagram of the modified IEEE 13-bus system, with DERs placed according to the basecase configuration.

While the network physical constraints (i.e., bus voltage and line thermal limits) are provided by the IEEE 13-bus specifications, the DERs placement and capacity are customized, and 3 configurations are considered, as described below:

* Basecase: Only conventional dispatchable generators are considered, such that $\mathcal{D}\coloneqq\{4,7,8\},\mathcal{R}=\varnothing,\mathcal{B}=\varnothing$. The generators adopt the following parameterization: $s^\mathrm{max}_4=s^\mathrm{max}_7=\SI{500}{\kilo\voltampere},s^\mathrm{max}_8=\SI{950}{\kilo\voltampere}$, the minimum power limits are set to zero $p^\mathrm{min}_i=0,\forall i\in\mathcal{D}$, and the maximum power limits match the apparent power limits $p^\mathrm{max}_i=s^\mathrm{max}_i,\forall i\in\mathcal{D}$. The generation cost coefficients are $\gamma_4=\SI[per-mode=symbol]{4}{\EUR\per\kWh}, \gamma_7=\SI[per-mode=symbol]{4.5}{\EUR\per\kWh}$ and $\gamma_8=\SI[per-mode=symbol]{5}{\EUR\per\kWh}$. When not stated otherwise, this model is used in simulations.

* Renew1: The conventional generator at bus 8 in the Basecase configuration is substituted by a RES, with the same nominal capacity and zero production cost, such that $\mathcal{D}=\{4,7\},\mathcal{R}=\{8\},\mathcal{B}=\varnothing$.

* Renew2: All conventional generators in the Basecase configuration are replaced by RESs with the same nominal capacities and zero production costs, such that $\mathcal{D}=\varnothing,\mathcal{R}=\{4,7,8\},\mathcal{B}=\varnothing$.

* Bat1: A BESS is connected at bus 1, with the following parameterization: $\chi_1^\mathrm{min}=0, \chi_1^\mathrm{max}=\SI{4800}{\kWh}, p^\mathrm{max}_1=\SI{500}{\kilo\watt}, s^\mathrm{max}_1=\SI{500}{\kilo\voltampere}$. All other generators from the Basecase configuration remain unchanged, such that $\mathcal{D}\coloneqq\{4,7,8\},\mathcal{R}=\varnothing,\mathcal{B}=\{1\}$. The battery cost factor is 0.

The implementation of the DDPG-based RL agent is briefly discussed below. Both the actor and the critic network have two hidden layers with 512 neurons that employ the Rectified Linear Unit (ReLU) activation function. The output layer of the actor employs a sigmoid function to bound the actions. The constant for the soft update of the target networks is $\tau=0.005$. We train the RL agent with a minibatch size of $10^6$ and the following exploration noise parameters: $z_\mathrm{i}=0.5, z_\mathrm{f} = 0.1$ and $m_\mathrm{tot}=100$. The qcrAdam optimizer is used to find the neural network weights with a learning rate $\alpha=0.001$ for the actor and $\beta=0.002$ for the critic. Three potential designs of the reward function and the shield are considered to examine the benefits and performance of the safety layer:

* Unsafe Reinforcement Learning (uRL): The VPP constraints are not taken into account, and the shield is bypassed. The reward consists only of day-ahead market revenues, more precisely $r(s,a) = r_\mathrm{da}$.

* Shielded Reinforcement Learning (shRL): The reward is designed as $r(s,a) = r_\mathrm{da}-c_\mathrm{vpp}$, with the constraints of the VPP taken into account via the active power shield (<ref>). However, no negative reward is passed to the agent for the shield activation.

* Safe Reinforcement Learning (sRL): In this case, the shield is active and the agent receives a negative reward for its activation, i.e., $r(s,a) = r_\mathrm{da} - c_\mathrm{vpp} - c_\mathrm{shd}$. Here, the goal is to train an agent that can place strategic bids without violating the physical limits of the VPP.

Performance of the three different RL models, with all values given in /day.

Model     Day-Ahead     Balancing Cost     Net Market Profit     $c_\mathrm{vpp}$     $c_\mathrm{sh}$

uRL     162819.1     67163.8      95655.3     -     -

shRL     157970.1     0     157970.1     93927.46     -

sRL     154594.4     0     154594.4     93052.12     10E-05

§.§ Safe vs. Unsafe RL-based VPP Bidding Strategies

The three different models listed above are evaluated on the Basecase scenario by assessing the average reward generated during the algorithm testing phase, which comprises $500$ episodes consisting of $24$ steps (hours). Table <ref> showcases the performance of these models. It can be observed that the uRL model achieves the highest market revenues. However, as it by definition does not take into account the physical constraints of the grid, it often bids a higher quantity than the VPP can physically provide. Hence, it is obliged to buy balancing power, which has a higher cost than the market price. In this work, the balancing price is assumed to be $20\%$ higher than the MCP. Therefore, the total profit generated by the uRL decreases. The shRL and sRL models achieve similar market profits, with sRL sacrificing some performance for safety. The shield component of the reward in sRL is negligible, which signals that the RL algorithm has learned how to take actions within the physical boundaries imposed by the nature of the VPP. In Fig. <ref>, it can be seen that during the training phase of sRL the shield is activated very often, while in the testing phase it is not needed anymore. This is clearly not the case for shielded RL, as showcased in the same figure. The lower shield activation observed during the testing phase in the top plot can be attributed to a decreased level of exploration noise.

Shield activation during training and testing for (a) Shielded RL design (top) (b) Safe RL design (bottom).

§.§ Bidding Strategy Analysis

In this section, we analyze the bidding strategies that the sRL agent adopts in the basecase scenario to achieve higher profits in the day-ahead market, as showcased in Figure <ref>.

Bidding behavior of the safe RL agent over 24 hours.

The VPP has the option to both buy and sell power on the day-ahead market, depending on the internally available generation and the consumption requirements of its customers (inflexible loads at nodes $i\in\mathcal{L}$).

The VPP is expected to dispatch internal generators whose production cost is lower than the MCP, and to buy the potential remaining energy needed to supply its loads in the day-ahead market. This strategy can clearly be observed in several hours in Fig. <ref>. In the first 6 hours of the day, the MCP is lower than the PMC of the VPP. Consequently, the VPP is buying power from the exchange. The amount of power varies hour by hour, as the VPP in most cases buys the exact amount needed to supply its loads. In case there is a generator with a unit PMC lower than the MCP, the amount of power bought is lower. During the other hours, the MCP is higher than the PMC, and hence, selling power becomes profitable. The VPP still provides the necessary power to its loads and generates extra power to sell to the market. The amount of power sold is limited either by its physical constraints

or by the demand on the market.

Besides learning to bid based on the difference between the MCP and the PMC of the VPP, the RL agent also learns to act as a price maker in the market.

By learning the expected MCP and the bids supporting it, the agent can place bids between the expected MCP and the next rival participant's bid in order to increase the MCP and maximize its profits. This process is of a delicate nature, as a bid that is too high might become out-of-the-money, leading to the VPP not being cleared.

The bidding behavior of the RL agent at hour 22 and its impact on the MCP is illustrated in Fig. <ref>. It suggests that the VPP acts as a price maker in this hour. Indeed, here, the range of feasible bids for the VPP is $\SI[per-mode=symbol]{0.5}{\EUR\per\kilo\watt}$ and $\SI[per-mode=symbol]{8}{\EUR\per\kilo\watt}$, as bidding above $\SI[per-mode=symbol]{8}{\EUR\per\kilo\watt}$ would get the VPP out of the merit order. A price maker VPP with perfect information on the rival participants' bids would bid as close as possible to the upper bound of this range, i.e., at $\SI[per-mode=symbol]{8}{\EUR\per\kilo\watt}$. We observe that by learning to increase its price bid within this competitive range, the RL agent manages to influence the MCP for most of the episodes, considerably increasing its net profit.

MCP with strategic VPP (RL agent), without VPP and with price-taker VPP.

§.§ Impact of Renewable Energy Sources and Batteries

In this section, we analyze quantitatively and qualitatively the impact of the introduction of RESs and BESSs into the VPP in terms of the bidding behavior of the RL agent and the total average net profits of the VPP, as summarized in Table <ref> for the different DERs configurations considered. The fields “w FC” and “w/o FC” indicate whether the renewable forecast is included in the agent state or not. Note that a different load curve compared to the previous subsections is used.

Performance of the considered DER configurations.

2*Configuration     2*Basecase     2cRenew1     2c|Renew2     2*Bat1

w FC     w/o FC     w FC     w/o FC

Profit $[\EUR{}/\mathrm{day}]$      44558     72983     37202     88681     45505     56756

Firstly, introducing RESs in the Renew1 and Renew2 configurations results in decreased internal production costs and lower power outputs due to reduced available production compared to the Basecase configuration. In both the Renew1 and Renew2 configurations, despite lower internal generation that can be sold in the day-ahead market, the RL agent learns to bid strategically to keep the MCP high and take advantage of the increased price difference between the MCP and the PMC of the VPP, yielding a significant increase in its net profits. We observe that the increase in net profits follows the increase in installed RESs between the Baseline, Renew1, and Renew2 configurations.

Secondly, the RL agent must deal with the variability and stochasticity of the power generation introduced by the RESs in the Renew1 and Renew2 configurations. During the bidding process, the only available information is the (deterministic) power generation forecast of the RESs $p^\mathrm{res}_i,\forall i\in\mathcal{R}$. We observe that omitting the generation forecast (w/o FC case) from the state vector $s$ of the RL agent leads to a considerable decrease in profitability in both Renew1 and Renew2 configurations. This result can intuitively be explained by the lack of knowledge about the VPP power output, prompting the RL agent to place conservative bids that ultimately result in robust but suboptimal behavior.

Finally, introducing a BESS in the Bat1 configuration increases the flexibility of the VPP compared to the Basecase scenario. The battery allows an increase in energy consumption in the first 6 hours of the day, which is stored and sold later when MCPs are higher. This observation suggests that the RL agent learns to strategically utilize battery flexibility, leading to additional arbitrage revenues and higher profits.

§ CONCLUSION

This paper introduces a novel safe RL algorithm for the strategic bidding of VPPs in electricity markets. Numerical results show that the agent successfully learns a safe and highly competitive bidding strategy, as reflected by the low activation of the safety shield during testing and the increased profits of the agent compared to the price-taking bidding strategy. Additionally, we demonstrate that introducing the shield-activation penalty and incorporating predictive information, such as renewable energy forecast, significantly improves the optimality and safety of the agent's actions.
These findings demonstrate the potential of using safe RL techniques to optimize the bidding strategies of VPPs while respecting the network constraints and limits of internal components.

Our forthcoming research will focus on evaluating the proposed algorithm on larger-scale test cases, therefore testing the scalability of the approach. This effort will entail deriving computationally more efficient methodologies for designing the safety shield. Additionally, we intend to assess the suitability of more advanced deep RL techniques, such as the soft actor-critic approach, for enhancing competitive bidding.