In practice the computational effort required to compute $u^{n+1}$ in (<ref>) is not significantly larger than the solution of (<ref>) or other previous symmetric low-regularity methods for the NLSE. In particular, a similar analysis to that presented in <cit.>, <cit.> and [50] shows that the implicit equations can be solved efficiently using fixed point iterations, and that the number of iterations required is independent of the number of spatial discretisation points.

One can find the coefficients $ b $ for the scheme (<ref>) such that it is of the form given by (<ref>)

For the scheme of order two $ (r=1) $ and $ n =2 $, one has to consider:

U_k^n,1(τ, u)     = ∑_T ∈^1,k_0(R) Υ^p_mid( T)(u,τ)/S(T) (Π_mid^n,1   T_0 )(τ)

= Υ^p_mid( T_0)(u,τ)/S(T_0) (Π_mid^n,1   T_0 )(τ) + ∑_k = -k_1 + k_2 + k_3 Υ^p_mid( T_1)(u,τ)/S(T_1) (Π_mid^n,1   T_1 )(τ)

+∑_k=-k_1+k_2+k_3 -k_4+k_5 Υ^p_mid( T_2)(u,τ)/S(T_2) ( Π^n,1_mid   T_2 )(τ)

+∑_k_1-k_2-k_3 +k_4+k_5= k Υ^p_mid(  T_3)(u,τ)/S(T_3) (Π^n,1_mid   T_3 )(τ)

From the definition of the symmetry factor, one has

S(T_2) =  1 ×2  = 2,   S(T_3) = 2 ×2= 4,

for $ S(T_2) $ the factor one corresponds to the fact that for the node on top of the first blue edges the symmetry factor is one. Indeed, the trees on top of the brown edges are different: a leaf decorated by $ k_5 $ is different from a tree having three leaves.
Moreover, one has:

Υ^p_mid( T_j)(u,τ)     =Υ^p( T_j)(1/2 ( e^i τk^2 u^n+1 + u^n )),   j ∈{2,3 },

Υ^p( T_2)(u)  = 4 u̅_k_1 u_k_2 u_k_3 u̅_k_4 u_k_5,
Υ^p( T_3)(u)  = 4 u_k_1 u̅_k_2 u̅_k_3 u_k_4 u_k_5.

The factor $ 4 $ in both expressions comes from the two brown edges that appear twice inside the decorated trees $ T_2 $ and $ T_3 $. For the term $ (\Pi_{\text{\tiny{mid}}}^{2,1}   T_1 )(\tau)  $, we proceed with interpolation at two nodes ($a_0=0,a_1=1$):

p_2(s,τ) = 1 + s/τ ( e^i s ℒ_low - 1 ),

ℒ_dom = 2 k_1^2,   ℒ_low =  k^2 - k_1^2 -k_2^2 -k_3^2.

We obtain

^k,1_o_2 ( (Π_mid^2,0   T̃_1 )(·,τ),2)(τ)     = -i	∫_0^τ e^i s ℒ_dom  ds
-i	∫_0^τ s e^i s ℒ_dom  ds ( e^i τℒ_low  -1 /τ  )

- i τφ_1(2 i τk_1^2) - i τφ_2(2 i τk_1^2) (e^ik^2 - k^2_1 - k_2^2 - k_3^2 - 1 )

∑_k = -k_1 + k_2 + k_3 Υ^p_mid( T_1)(u,τ)/S(T_1) (Π_mid^2,1   T_1 )(τ)
= ( - i τφ_1(2 i τk_1^2) - i τφ_2(2 i τk_1^2) (e^ik^2 - k^2_1 - k_2^2 - k_3^2 - 1 ) )

×( e^-i τk_1^2 u̅_k_1^n+1 + u̅_k_1^n ) 1/2 ( e^i τk_2^2 u_k_2^n+1 + u_k_2^n ) 1/2 ( e^i τk_3^2 u_k_3^n+1 + u_k_3^n )

For the decorated tree $ T_2 $, we have

(	Π^2,1_mid   T_2 )(τ) = e^-i τk^2 	^k,1_o_2 ((Π^2,0_mid   F_2 )(·,τ) ,2 )(0,τ)

F_2 = ℐ_(𝔱_1,1)(λ_k_4) ℐ_(𝔱_1,0)(λ_k_5)  T_1

(Π^2,0_mid   ℐ_(𝔱_1,1)(λ_k_4) ℐ_(𝔱_1,0)(λ_k_5) )(s,τ) =  e^i s (k_4^2 -k_5^2)

and for $ \tilde{k} = -k_1 + k_2 + k_3 $

(Π^2,0_mid   T_1 )(s,τ)     = 1/2   e^-i s k̃^2 (  ^k̃,0_o_2 ( (Π_mid^2,-1   T̃_1 )(·,τ),2)(s,τ) .
. +  ^k̃,0_o_2 ( (Π_mid^2,-1   T̃_1 )(·,τ),2)(s,0) ),

Now, because of $ n=2 $, we perform a full Taylor expansion which gives

^k̃,0_o_2 ( (Π_mid^2,-1   T̃_1 )(·,τ),2)(s,τ) = -i ∫_τ^s ds ( 1 + e^i τ(k̃^2 + k_1^2 - k_2^2 - k_3^2)/2 )

We obtain
(Π^2,0_mid   T_1 )(s,τ) = - i (2s -τ)/2 e^-i s k̃^2  ( 1 + e^i τ(k̃^2 + k_1^2 - k_2^2 - k_3^2)/2 )

Therefore, we find

(Π^2,0_mid   F_2 )(s,τ) = - i (2s-τ)/2 e^i s( k_4^2 - k_5^2 -k̃^2)  ( 1 + e^i τ(k̃^2 + k_1^2 - k_2^2 - k_3^2)/2 )

and by performing again a full Taylor expansion:

(	Π^2,1_mid   T_2 )(τ)     = - ∫^τ_0	(2s -τ)/2 ds
e^-i τk^2  	( 1 + e^i τ(k^2 + k_4^2 - k_5^2 - k̃^2)/2 )	( 1 + e^i τ(k̃^2 + k_1^2 - k_2^2 - k_3^2)/2 )

= 0

because we have

∫^τ_0	(2s -τ)/2 ds = 0.

A similar computation shows that

(	Π^2,1_mid   T_3 )(τ)     = 0.

The local error analysis follows from  the proof of <cit.>.

§.§.§ The Korteweg–de Vries equation

The Korteweg–de Vries (KdV) equation is given by

∂_t u + ∂_x^3 u = 1/2 ∂_x u^2

It fits into the general framework with

\begin{equation*}\label{kgrDo}
\begin{aligned}
\mathcal{L}\left(\nabla, \frac{1}{\varepsilon}\right)  = i \partial_x^3, \quad \alpha = 1   \quad \text{and}\quad   p(u,\overline u) =  p(u) = i  \frac12 u^2.
\end{aligned}
\end{equation*}

Here $ \CL = \lbrace \Labhom_1, \Labhom_2 \rbrace $, $ P_{\Labhom_1} = - \lambda^3 $ and $ P_{\Labhom_2} =  \lambda^3 $ . Then, we denoted by $ \<thick> $ an edge decorated by $ (\Labhom_1,0) $ and by $\<thin>$ an edge decorated by $ (\Labhom_2,0) $.
Following the formalism given in [12], one can provide the rules that generate the trees obtained by iterating the Duhamel formulation:

R(thin>) = {(thick>,thick>) } ,   R(thick>) = {(thin>), ()} .

The general framework (<ref>) derived in Section <ref> builds the foundation of the  first- and second-order resonances based schemes presented below for the KdV equation (<ref>). The structure of the schemes depends  on the regularity of the solution.

For the KdV equation (<ref>) the general midpoint scheme  (<ref>) takes at first order the form

\begin{equation}
\begin{aligned}\label{schemeKdV1}
u^{\ell+1} &=  e^{-\tau \partial_x^3} u^\ell + \frac{1}{24} \left(e^{-\tau\partial_x^3 }\partial_x^{-1}  u^{\ell} + \partial_x^{-1} u^{\ell+1} \right)^2 \\ & - \frac{1}{24} e^{-\tau\partial_x^3} \left(\partial_x^{-1}  u^{\ell} + e^{\tau \partial_x^3} \partial_x^{-1} u^{\ell+1} \right)^2\end{aligned}
\end{equation}

with a local error  of order $\mathcal{O}\Big(
\tau^2 \partial_x^2 u
\Big)$
at  first-order
with a local error  of order $\mathcal{O}\Big(
\tau^3 \partial_x^4 u
\Big)$ at second order.

Note that this schemes has been obtained in [50]. It was shown that this scheme is of even order for higher regularity in $ H^4 $  (see <cit.>). By embedding this scheme into our general framework, we know that it has the same local error analysis as the second-order scheme introduced in [15].
\begin{equation}
\begin{aligned}\label{schemeKdV}
u^{\ell+1} &=  e^{-\tau \partial_x^3} u^\ell + \frac16 \left(e^{-\tau\partial_x^3 }\partial_x^{-1} u^\ell\right)^2 - \frac16 e^{-\tau\partial_x^3} \left(\partial_x^{-1} u^\ell\right)^2\\& +\frac{\tau^2}{4} e^{- \tau \partial_x^3}\Psi\big(i \tau \partial_x^2\big)  \Big(\partial_x \Big(u^\ell \partial_x (u^\ell u^\ell)\Big)\Big)
\end{aligned}
\end{equation}

with a local error  of order $\mathcal{O}\Big(
\tau^3 \partial_x^4 u
\Big)$ and  a suitable filter function $\Psi$ satisfying

Ψ= Ψ(i τ∂_x^2 ),   Ψ(0) = 1,   ‖τΨ(i τ∂_x^2) ∂_x^2 ‖_r ≤1.

The proof follows the line of argumentation to the analysis for the Schrödinger equation.
For the first-order scheme, we have

U_k^n,0(τ, v)     = Υ^p_mid( T_0 )(τ, v)/S(T_0) Π^n,0_mid ( T_0 )(τ)

+  ∑_k = k_1 + k_2 Υ^p_mid(  T_1 )(τ, v)/S(T_1) Π^n,0_mid (  T_1 )(τ).

where the trees of interest are

^0,k_0(R)= {T_0, T_1,   k_i ∈^d },   T_0 = [scale=0.2,baseline=-5]
(root) at (0,1);
(tri) at (0,-1);
[kernels2] (tri) – (root);
[var] (rootnode) at (root) $ k $;
[not] (trinode) at (tri) ;
and  T_1 =   [scale=0.2,baseline=-5]
(root) at (0,2);
(root1) at (0,0);
(tri) at (0,-2);
(t1) at (-1,4);
(t2) at (1,4);
[kernels2] (t1) – (root);
[kernels2] (t2) – (root);
[symbols] (root) – (root1);
[kernels2] (tri) – (root1);
[not] (rootnode) at (root) ;
[not] (rootnode) at (root1) ;
[not] (trinode) at (tri) ;
[var] (rootnode) at (t1) $ k_{\tiny{1}} $;
[var] (trinode) at (t2) $ k_2 $;

and in symbolic notation  takes the form

T_1 =_(_1,0)( _(_2,0)( λ_k F_1))   F_1 = _(_1,0)(  λ_k_1)  _(_1,0) ( λ_k_2)   with  k = k_1+k_2.

For the first term we readily obtain that

Υ^p_mid( T_0 )(τ, v)/S(T_0) Π^n,0_mid ( T_0 )(τ)
= e^-i τk^3 v̂_k.

It remains to compute the second term.  Note that thanks to (<ref>) we have that

Π^n,0 ( T_1 )(τ)    =  e^-i τk^3 Π^n,0 (_(_2,0)( λ_k F_1) )(τ)

=  e^-i τk^3 𝒦_(_2,0)^k,0 ( Π^n,-1 (F_1), n)(τ)

= e^-i τk^3 𝒦_(_2,0)^k,0 ( e^i ξ(-k_1^3 - k_2^3), n )(τ).

where we have used for the third line

(Π^n,-1 F_1)(τ)   = (Π^n,-1  _(_1,0) ( λ_k_1)) (τ) (Π^n,-1 _(_1,0)(  λ_k_2)) (τ)  = e^-i τk_1^3 e^-i τk_2^3.

Next  we observe that

P_(_2,0)(k) - k_1^3-k_2^3  =  k^3 - k_1^3 - k_2^3 = 3 k_1 k_2 (k_1+k_2)

such that

1/ P_(_2,0)(k) - k_1^3-k_2^3

can be mapped back to physical space. Therefore, we set

ℒ_dom =  P_(_2,0)(k) - k_1^3-k_2^3  = 3 k_1 k_2 (k_1+k_2)

and integrate all frequencies exactly. This implies

Π^n,0 ( T_1)(τ)
= e^-i τk^3 i (k_1+k_2) /3i  k_1 k_2 (k_1+k_2) (e^i τ(k^3 - k_1^3 -k_2^3)-1)

=1 /3  k_1 k_2  (e^- i τ( k_1^3 + k_2^3)- e^-i τk^3 ).

Together with (<ref>) this yields the scheme (<ref>).
For the second-order scheme, we first notice that

Π^n,0_mid ( T_0 )(τ) = Π^n,1_mid ( T_0 )(τ),   Π^n,0_mid ( T_1 )(τ) = Π^n,1_mid ( T_1 )(τ)

Indeed, for the tree $ T_1 $, we perform an exact integration without any discretisation.
Then, we need to take into account the following trees
^1,k_0(R) = {T_0, T_1, T_2,   k_i ∈^d },   T_2 = [scale=0.2,baseline=-5]
(root) at (0,2);
(root1) at (0,0);
(tri) at (0,-2);
(t1) at (-1,4);
(t11) at (-2,6);
(t12) at (-3,8);
(t13) at (-1,8);
(t2) at (1,4);
[kernels2] (t11) – (t13);
[kernels2] (t11) – (t12);
[kernels2] (t1) – (root);
[symbols] (t1) – (t11);
[kernels2] (t2) – (root);
[symbols] (root) – (root1);
[kernels2] (root1) – (tri);
[not] (rootnode) at (root) ;
[not] (rootnode) at (root1) ;
[not] (trinode) at (tri) ;
[not] (trinode) at (t1) ;
[var] (rootnode) at (t12) $ k_{\tiny{1}} $;
[var] (rootnode) at (t13) $ k_{\tiny{2}} $;
[var] (trinode) at (t2) $ k_3 $;

Then one can proceed as in the second-order schemes for the Schrödinger equation to show its contribution is zero that is

Π^4,1_mid ( T_2 )(τ) = 0.

§ NUMERICAL EXPERIMENTS

We now test the practical performance of our new symmetric schemes in practical experiments evaluating both their low-regularity convergence properties and their ability to correctly preserve constants of motion in the relevant equations. In fitting with our above construction our spatial discretisation is a Fourier spectral method throughout with $M$ modes. In  order to understand the low-regularity convergence properties of our methods we follow [56] and consider the following types of initial data:

* Smooth inital data,
\begin{align}\label{eqn:numerical_experiments_smooth_data}
\end{align}

* Low-regularity initial data $u_0$ of the following form. Firstly we choose a vector sampled from a uniform distribution $U_{m}\sim U([0,1]+i[0,1]), m=-M/2+1,\dots, M/2$ and then we define

\begin{align}\label{eqn:numerical_experiments_low_reg_data}
\end{align}

for a given value of $\vartheta>0$, which corresponds to $u_0\in H^{\vartheta}$.

Both choices of initial data are rescaled such that $u_0\mapsto u_0/\|u_0\|_{L^2}$.

§.§ The Nonlinear Schrödinger equation

To begin with we look at our new symmetric integrators for the NLSE ((<ref>) and (<ref>)). In the following numerical experiments we compare the performance of our new schemes to the following state-of-the-art reference schemes for the NLSE:

* The Strang splitting [51], as an example of a classical symmetric numerical technique;

* The first and second order low-regularity schemes introduced by Ostermann & Schratz [56] and Bruned & Schratz <cit.> respectively, as examples of asymmetric low-regularity integrators;

* The symmetrised low-regularity integrator introduced by Alama Bronsard [4], as an example of previous structure preserving resonance-based schemes.

In the following numerical experiments our reference solutions were computed with $M=2^{14}$ Fourier modes and a time step $\tau=10^{-6}$ with the symmetrised method from [4]. In the first instance we consider, in Figures <ref> & <ref> the convergence properties of our methods for initial data of various levels of regularity. In all of these experiments the number of spatial discretisation modes was taken to be $M=1024$ and the initial data chosen according to (<ref>) & (<ref>). In Figure <ref> we observe that our new methods have exactly the predicted convergence properties at those levels of regularity: The integrator (<ref>) is optimal for first order convergence in the sense of regularity, meaning it converges at first order in $H^1$ with $H^2$ data, while the integrator (<ref>) is optimally convergent in the sense of regularity up to second order meaning it converges at $\mathcal{O}(\tau)$ in $H^1$ for data in $H^2$ and at $\mathcal{O}(\tau^2)$ in $H^1$ for data in $H^3$.

$H^2$-initial data, $\vartheta=2$ in (<ref>).

$H^3$-initial data, $\vartheta=3$ in (<ref>).

$H^1$-error at $T=1$ as a function of the timestep $\tau$ for low-regularity initial data.

Strang splitting suffers from significant order reduction in those regimes which is apparent in those graphs. Particularly important is also to observe in Figure <ref> that the previous symmetric resonance-based scheme does not converge at second order for initial data in $H^3$ and indeed requires (cf. Figure <ref>) at least solutions in at $H^4$ in order to achieve this convergence rate. This is exactly in fitting with the results in [4] and highlights the important novelty in the design of our new symmetric scheme (<ref>). The behaviour of the Strang splitting observed in those experiments matches exactly with the convergence analysis given by [30].

$H^4$-initial data, $\vartheta=4$ in (<ref>).

$C^\infty$-initial data (<ref>).

$H^1$-error at $T=1$ as a function of the timestep $\tau$ for more regular initial data.

Following the study of the low-regularity convergence properties of our new methods we turn to testing its long-time structure preservation properties. For this we consider two first integrals of the cubic NLSE, the normalisation

I_0^[NLSE][u]=∫_𝕋|u|^2 dx,

and the energy

I_1^[NLSE][u]=∫_𝕋 |∇u|^2+1/2|u|^4 dx.

Symmetric numerical schemes are typically unable to preserve such conservation laws exactly, however it is known for the ODE case <cit.> (and also observed numerically for the PDE case, for example in [19]) that symmetric methods can exhibit very good approximate long-time preservation of such first integrals. In the following numerical experiments we test this behaviour by looking at the error in these quantities for a fixed timestep $\tau=0.05$ ($M=1024$ unless otherwise stated) over a long time interval, much larger than $\mathcal{O}(1/\tau)$. Firstly, in figures <ref>&<ref> we observe that the normalisation appears to be preserved really well, in particular (the example is representative of a host of numerical experiments for various time steps which we performed) the preservation is much better than previous asymmetric resonance-based schemes and both (<ref>) & (<ref>) are able to outperform even the symmetrised scheme introduced in [4]. We note that the Strang splitting conserves quadratic first integrals, i.e. the normalisation, to machine accuracy and thus undoubtedly outperforms our schemes on the level of normalisation preservation.

Long-time interval $t=n\tau\in[0,4000]$, $M=128$.

Long-time interval $t=n\tau\in[0,4000]$, $M=1024$.

Magnification of $t=n\tau\in[0,50]$, $M=1024$.

Error in the normalisation $\|u^{n}\|_{L^{2}}$, for time step $\tau=0.02$ and $H^2$ data.

Long-time interval $t=n\tau\in[0,4000]$, $M=128$.

Long-time interval $t=n\tau\in[0,4000]$, $M=1024$.

Magnification of $t=n\tau\in[0,50]$, $M=1024$.

Error in the normalisation $\|u^{n}\|_{L^{2}}$, for time step $\tau=0.02$ and $C^\infty$ data.

Our final numerical result for the NLSE, presented in Figure <ref>&<ref>, shows the error in the NLSE energy over a long time interval, for a fixed time step $\tau=0.05$. Albeit rigorous theory exists for the ODE case [36], there is again no theoretical guarantee for the long-time preservation of the energy under symmetric methods. Indeed in practical experiments it can be seen that symmetric schemes are able to clearly outperform asymmetric integrators in the long-time approximate energy preservation. For the Strang splitting this behaviour was rigorously analysed in [28] where a CFL condition was necessary to guarantee long-time approximate energy preservation beyond the realms of forward error analysis. This CFL condition is indeed observed even for smooth data in our experiments. Perhaps somewhat surprisingly our new symmetric resonance-based scheme do not appear to suffer from comparable CFL conditions and, as expected, perform well for both smooth and low-regularity solutions.

Error in the Hamiltonian, for time step $\tau=0.02$ and $H^2$ data.

Error in the Hamiltonian, for time step $\tau=0.02$ and $C^\infty$ data.

§.§ The Korteweg–de Vries equation

As a final numerical example we consider the resonance-based midpoint rule (<ref>) introduced in [50] which our midpoint iterates (section <ref>) are able to recover. This rule has excellent low-regularity convergence properties and at the same time is able to conserve momentum and energy of the KdV equation over long times even in the low-regularity regime. With permission we reproduce some of the numerical results presented in [50] to recall the favourable properties of the method (<ref>), for further evaluation and results the interested reader is referred to [50]. Here we compare the performance of our new schemes to the following state-of-the-art reference schemes for the KdV:

* The Strang splitting [51], as an example of a classical symmetric numerical technique, by splitting the KdV equation into the two subproblems $\partial_t u=-\partial_x^3u$ and $\partial_t u=u\partial_x u$. In our implementation the resulting Burgers equation is solved using an explicit RK4 scheme with small time step $\tau_{RK4}=\tau10^{-3}$ such that in essence the error observed in the following experiments is only due to the splitting of the problem.

* The first and second order low-regularity schemes introduced by Hofmanova & Schratz [38] and Bruned & Schratz <cit.> respectively, as examples of asymmetric low-regularity integrators.

All reference values were computed with $M=2^{14}$ Fourier modes and a time step $\tau=10^{-6}$ with the second order method from [15].

To begin with, in Figures <ref>&<ref> we observe the excellent low-regularity convergence properties of our method. The classical Strang splitting suffers in this case from a CFL condition, which requires $\tau\lesssim 1/M$ in order to resolve the Burgers nonlinearity, as a result we chose to include numerical results with only a small number of Fourier modes, since the classical integrator was found to be unstable whenever the aforementioned CFL condition is not satisfied. Perhaps somewhat surprisingly we see that in practice the integrator (<ref>) appears to converge even faster than the method introduced by Bruned & Schratz for solutions in $H^3$, this was also observed in [50].

$H^3$-initial data, $\vartheta=3$ in (<ref>).

$H^5$-initial data, $\vartheta=5$ in (<ref>).

$H^1$-error at $T=1$ as a function of the timestep $\tau$ for low-regularity initial data.

$H^1$-error at $T=1$ as a function of the timestep $\tau$ for $C^\infty$-initial data (<ref>).

$H^5$-initial data, $\vartheta=5$ in (<ref>).

$C^\infty$-initial data (<ref>).

Error in the momentum $I_0^{[KdV]}[u]$ for $\tau=0.005$ and $M=64$.

As a next step we can consider the momentum, which is a quadratic first integral in the flow.

\begin{align*}
I_0^{[KdV]}[u]=\int_{\mathbb{T}} u^2  d x
\end{align*}

Seeing as the Strang splitting is a symplectic integrator we expect the momentum to be preserved to machine precision for $C^{\infty}$ solutions which is indeed observed in Figure <ref>. Note in [50] it was actually shown that the method (<ref>) is symplectic and preserves the momentum exactly, which explains the preservation of momentum to machine accuracy in that same graph. On the other hand, for low-regularity data we see in Figure <ref> that this preservation is lost in the Strang splitting for rougher data, but that the symmetric scheme (<ref>) is able to deal with rough solutions as well.

[1]
Y. Alama Bronsard, Y. Bruned, K. Schratz,
Approximations of dispersive PDEs in the presence of low-regularity randomness.

[2]
Y. Alama Bronsard, Y. Bruned, K. Schratz,
Low regularity integrators via decorated trees.

[3]
Y. Alama Bronsard,
Error analysis of a class of semi-discrete schemes for solving the Gross-Pitaevskii equation at low regularity.

[4]
Y. Alama Bronsard,
A symmetric low-regularity integrator for the nonlinear Schrödinger equation.

[5]
I. Ampatzoglou, C. Collot, P. Germain,
Derivation of the kinetic wave equation for quadratic dispersive problems in the inhomogeneous setting

[6]
G. Bai, B. Li, and Y. Wu.
A constructive low-regularity integrator for the 1d cubic nonlinear Schrödinger equation under the Neumann boundary condition.
IMA J. Numer. Anal. doi: 10.1093/imanum/drac075.

[7]
V. Banica, G. Maierhofer, K. Schratz,
Numerical integration of Schrödinger maps via the Hasimoto transform.

[8]
G. Benettin and A. Giorgilli,
On the Hamiltonian interpolation of near-to-the identity symplectic mappings with application to symplectic integration algorithms.
J. Stat. Phys., 74, (1994), 1117–1143.

[9]
P. B. Bochev and C. Scovel,
On quadratic invariants and symplectic structure.
BIT Numer. Math., 34, no. 3, (1994), 337–345.

[10]
Y. Bruned, A. Chandra, I. Chevyrev,
M. Hairer,
Renormalising SPDEs in regularity structures.
J. Eur. Math. Soc. (JEMS), 23, no. 3, (2021), 869-947.

[11]
Y. Bruned, K. Ebrahimi-Fard,
Bogoliubov type recursions for renormalisation in regularity structures.

[12]
Y. Bruned, M. Hairer, L. Zambotti.
Algebraic renormalisation of regularity structures.
Invent. Math. 215, no. 3, (2019), 1039–1156.

[13]
Y. Bruned, M. Hairer, L. Zambotti.
Renormalisation of Stochastic Partial Differential Equations.
EMS Newsletter 115, no. 3, (2020), 7–11.
doi: 10.4171/NEWS/115/3http://dx.doi.org/10.4171/NEWS/115/3.

[14]
Y. Bruned, D. Manchon.  Algebraic deformation for (S)PDEs. To appear in J. Math. Soc. Japan.

[15]
Y. Bruned, K. Schratz.  Resonance based schemes for dispersive equations via decorated
trees. Forum of Mathematics, Pi, 10, E2.

[16]
N. N. Bogoliubow, O. S. Parasiuk.
Über die Multiplikation der Kausalfunktionen in der
Quantentheorie der Felder.
Acta Math. 97, (1957), 227–266.

[17]
J. C. Butcher,
An algebraic theory of integration methods.
Math. Comp. 26, (1972), 79–106.

[18]
J. Cao, B. Li, and Y. Lin,
A new second-order low-regularity integrator for the cubic nonlinear Schrödinger equation,
IMA Journal of Numerical Analysis (2023).

[19]
E. Celledoni, D. Cohen, B. Owren,
Symmetric exponential integrators with an application to the cubic Schrödinger equation.
Found. Comput. Math.
8, (2008), 303–317.

[20]
P. Chartier and A. Murua,
Preserving first integrals and volume forms of additively split systems.
IMA J. Numer. Anal. 27, (2007), 381–405.

[21]
D. Cohen, L. Gauckler,
One-stage exponential integrators for nonlinear Schrödinger equations over long times.
BIT 52, (2012), 877–903.
[22]
A. Chandra, M. Hairer.
An analytic BPHZ theorem for regularity structures.
[23]
A. Connes,  D. Kreimer,
Hopf algebras, renormalization and noncommutative geometry.
Comm. Math. Phys. 199, no. 1, (1998), 203–242.

[24]
A. Connes, D. Kreimer,
Renormalization in quantum field theory and the Riemann-Hilbert
problem I: the Hopf algebra structure of graphs and the main theorem.
Commun. Math. Phys. 210, (2000), 249–73.

[25]
M. Christ, Power series solution of a nonlinear Schrödinger equation. In Mathematical aspects
of nonlinear dispersive equations, volume 163 of Ann. of Math. Stud., pages 131–155. Princeton
Univ. Press, Princeton, NJ, 2007.

[26]
Y. Deng, Z. Hani,
Full derivation of the wave kinetic equation. To appear in Invent. Math.

[27]
Y. Deng, Z. Hani,
On the derivation of the wave kinetic equation for NLS.
Forum of Mathematics, Pi, 9, (2021), e6.
[28]
E. Faou,
Geometric Numerical Integration and Schrödinger Equations. European Math. Soc. Publishing House, Zürich 2012.

[29]
Y. Feng, G. Maierhofer, K. Schratz,
Long-time error bounds of low-regularity integrators for nonlinear Schrödinger equations.

[30]
L. Gauckler and C. Christian,
Splitting integrators for nonlinear Schrödinger equations over long times,  Found. Comput. Math. 10, no. 3, (2010), 275–302.

[31]
Z. Guo, S. Kwon, T. Oh,
Poincaré-Dulac normal form reduction for unconditional well-posedness of
the periodic cubic NLS, Comm. Math. Phys. 322, no. 1, (2013), 19–48.

[32]
M. Gubinelli,
Rough solutions for the periodic Korteweg-de Vries equation. Comm. Pure Appl. Anal.
11, no. 4, (2012),

[33]
M. Hairer,
A theory of regularity structures.
Invent. Math. 198, no. 2, (2014), 269–504.

[34]
K.  Hepp.
On the equivalence of additive and analytic renormalization.
Comm. Math. Phys. 14, (1969), 67–69.

[35]
E. Hairer, C. Lubich,
Symmetric multistep methods over long times.
Numer. Math. 97, (2004), 699–723.

[36]
E. Hairer, C. Lubich, G. Wanner,
Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Second edition, Springer, Berlin 2006.

[37]
M. Hochbruck, A. Ostermann,
Exponential integrators.
Acta Numer.
19, (2010),

[38]
M. Hofmanová, K. Schratz,
An exponential-type integrator for the KdV equation.
Numer. Math.
136, (2017),

[39]
H. Holden, K. H. Karlsen, K.-A. Lie, N. H. Risebro,
Splitting for Partial Differential Equations with Rough Solutions. European Math. Soc. Publishing House, Zürich 2010.

[41]
A. Iserles, G. R. W. Quispel  and P. S. P. Tse,
B-series methods cannot be volume-preserving.
BIT Numer. Math.
47, (2007),

[42]
A. Iserles, A. Zanna,
Preserving algebraic invariants with Runge–Kutta methods.
J. Comp. Appl. Math.
125, no. 1-2, (2000), 69–81.

[43]
G. Y., Kulikov,
Symmetric Runge–Kutta methods and their stability.
Russian Journal of Numerical Analysis and Mathematical Modelling
18, no. 1, (2003), 13–41.

[44]
J. D. Lawson,
Generalized Runge–Kutta processes for stable systems with large Lipschitz constants. SIAM J. Numer. Anal. 4:372–380 (1967).

[45]
B. Leimkuhler, S. Reich,
Simulating Hamiltonian dynamics. Cambridge Monographs on Applied and Computational Mathematics 14. Cambridge University Press, Cambridge, 2004.

[46]
B. Li, S. Ma, K. Schratz,
A semi-implicit low-regularity integrator for Navier-Stokes equations

[48]
B. Li, Y. Wu,
An unfiltered low-regularity integrator for the KdV equation with solutions below $\mathbf{H}^{1}$.

[49]
C. Lubich,
On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations.
Math. Comp.
77, no. 4, (2008), 2141–2153.

[50]
G. Maierhofer, K. Schratz,
Bridging the gap: symplecticity and low regularity in Runge–Kutta resonance-based schemes.

[51]
R.I. McLachlan, G.R.W. Quispel,
Splitting methods.
Acta Numer.
11, (2002), 341–434.

[52]
H. Z Munthe-Kaas, A. Lundervold,
On post-Lie algebras, Lie-Butcher series and moving frames.  Found. Comput. Math.
13, (2013), 583-613.

[54]
A. Ostermann, F. Rousset,  K. Schratz,
Error estimates of a Fourier integrator for the cubic Schrödinger equation at low regularity.  Found. Comput. Math.
21, (2021), 725-765.

[55]
A. Ostermann, F. Rousset,  K. Schratz,
Fourier integrator for periodic NLS: low regularity estimates via discrete Bourgain spaces. To appear in J. Eur. Math. Soc. (JEMS),
DOI 10.4171/JEMS/1275https://ems.press/journals/jems/articles/7697403.

[56]
A. Ostermann, K. Schratz,
Low regularity exponential-type integrators for semilinear Schrödinger equations.  Found. Comput. Math.
18, (2018), 731-755

[57] F. Rousset, K. Schratz, A general framework of low regularity integrators, SIAM J. Numer. Anal. 4, (2022) 127-152

[58] F. Rousset, K. Schratz, Convergence error estimates at low regularity for time discretizations of KdV. To appear in Pure and Applied Analysis  arXiv:2102.11125https://arxiv.org/abs/2102.11125.

[59]
J.M. Sanz-Serna, M.P. Calvo,
Numerical Hamiltonian Problems. Chapman and Hall,
London, 1994.

[60]
J.M. Sanz-Serna,
Runge-Kutta schemes for Hamiltonian systems.  BIT Numer. Math.
28, (1988), 877-883

[61]
Y. Wang, X. Zhao,
A symmetric low-regularity integrator for nonlinear Klein-Gordon equation, Math. Comp. 91, (2022) 2215-2245

[62]
Y. Wu, F. Yao,
A first-order Fourier integrator for the nonlinear Schrödinger equation on $\mathbb{T}$ without loss of regularity, Mathematics of Computation 91.335 (2022): 1213-1235

[63]
Y. Wu, X. Zhao,
Embedded exponential-type low-regularity integrators for KdV equation under rough data, BIT Numer. Math. 62, (2022) 1049–1090

[64]
W.  Zimmermann.
Convergence of Bogoliubov’s method of renormalization in momentum space.
Comm. Math. Phys. 15, (1969), 208–234.


# Deep Reinforcement Learning-Based
Robust Protection in Electric Distribution Grids

Dongqi Wu, Dileep Kalathil, and Le Xie Authors are with the Department of
Electrical and Computer Engineering, Texas A&M University, College Station,
TX-77840, USA. e-mail: {dqwu, dileep.kalathil<EMAIL_ADDRESS>was
supported in part by the NSF Grant CRII:CPS-1850206

###### Abstract

This paper introduces a Deep Reinforcement Learning based control architecture
for the protective relay control in power distribution systems. The key
challenge in protective relay control is to quickly and accurately detect
faults from other disturbances in the system. The performance of widely-used
traditional overcurrent protection scheme is limited by factors including
distributed generations, power electronic interfaced devices and fault
impedance. We propose a deep reinforcement learning approach that is highly
accurate, communication-free and easy to implement. The proposed relay design
is tested in OpenDSS simulation on the IEEE 34-node and 123-node test feeders
and demonstrated excellent performance from the aspect of failure rate,
robustness and response speed.

###### Index Terms:

Power Distribution Systems, Protective Relaying, Reinforcement Learning

##  I Introduction

This paper proposes and thoroughly test a novel Deep Reinforcement Learning
(Deep RL) based approach for the protective relay control design in the future
distribution grids. Recent developments in photovoltaic (PV) and power
electronic technology have led to a tremendous increase in the penetration of
distributed generation (DG) in the distribution grids. Distributed generation,
especially solar PVs, can provide a number of benefits to the power system
operation efficiency such as peak load reduction and improved power quality
[1]. However, DG and many emerging grid edge-level devices are increasing the
complexity of the interactions between the end users and the distribution grid
operators in a substantial manner. These additional complexities pose
significant challenges to the operation and protection of the distribution
grid.

Protective relays are the safeguards of distribution systems. The role of the
protective relays is to protect the grid from sustaining faults by
disconnecting the faulty segment from the rest of the grid. During the
operation, a relay monitors the power grid and look for patterns that
signifies faults. Typical measurements include current (over-current and
differential relay), voltage and current (distance relay), or electromagnetic
wave from transients (traveling-wave relay). In power distribution systems,
overcurrent relay is the most commonly used type of protection since many
other methods are impractical due to cost and infrastructure limitations.

However, it is very difficult for overcurrent relays to accommodate the vastly
different operation conditions in the real distribution grids. For example,
for feeder reclosers, the presence of DG within the feeder can reduce the
fault current measured at the recloser and make faults harder to detect. This
is because the fault current contribution from DG to the fault will make the
fault current measured at the fuse higher than the current at the recloser,
making coordination based on inverse-time curves difficult [2]. Moreover, even
in current distribution grids, factors like fault impedance and load profile
change are not taken into account in the traditional overcurrent protection
design, resulting in problems such as failing to detect faults near the end of
a feeder, a.k.a. under-reaching.

Traditional protective relays are also designed to function under two crucial
assumptions: (i) power flow is unidirectional from the substation towards the
end users, and (ii) the difference between operating conditions (currents and
voltages) between normal and abnormal conditions are measurable and
significant. With the increasing penetration of DG and edge devices, both
assumptions will be rendered invalid [2]. In fact, proper functioning of the
distribution protection is becoming a bottleneck in the deep integration of DG
for the future grid. This paper focuses on how to address the challenges in
protection design and operation for the distribution grid.

###  I-A Literature Review

There are many studies on improving the performance of protective relays. Most
of them focus on improving the performance of the commonly used overcurrent
relays by better fault detection [3] and coordination [4]. Neural networks are
used in setting the parameters of overcurrent relays [5]. Support Vector
Machine (SVM) can be trained to distnguish the normal and fault conditions
directly [6][7]. A recent work [8] uses tabular Q-learning to find the optimal
setting for overcurrent relays. Most proposed methods are still confined
within the framework of inverse-time overcurrent protection, which is not
enough for the future distribution grid with high DG and EV penetration [2].

Reinforcement Learning (RL) is a branch of machine learning that addresses the
problem of learning the optimal control policies for an unknown dynamical
system. RL algorithms using deep neural networks [9], known as Deep RL
algoritms, have made significant achievements in the past few years in areas
like robotics, games, and autonomous driving [10]. RL has also been applied to
various power system control problems including voltage regulation [11],
frequency regulation [12], market operation [13], power quality control [14]
and generator control [15]. Our previous paper [16] was the first work to use
deep RL for power system protection. It introdcued a sequential training
algorithm for the coordination between multiple RL-based relays in the
distribution system. A comprehensive survey of RL applications in power system
is detailed in a recent review paper [17].

###  I-B Main Contributions

The aim of this paper is to introduce a novel deep RL based solution to
address the protective relay control problem under the future distribution
grids with high distributed energy resources penetration. This approach
combines recent advancements in machine learning algorithms and deep domain
knowledge on power system operations. Main contributions are summarized as
follows:

  * •

Formulation of the optimal protective relay control problem as an RL problem.

  * •

A novel RL algorithm for protective relays, for reliable fault detection and
accurate coordination.

  * •

An open-source environment for the interface between the machine learning
packages and power system simulators.

The rest of this paper is organized as follows: Section II formulates the
protection control problem using the framework of RL. Section III introduces
our nested reinforcement learning algorithm for relay control problem. IV
presents the simulation environment and the test-bed cases used in training
and evaluation, Section V analyzes and discusses the simulation results.
Section VI summarizes and concludes the paper.

##  II Problem Formulation

In this section, we give a brief review of the basics of RL and formulate the
protective relay control problem using the RL framework.

###  II-A Relay Operation

Figure 1: Protective relays in a radial line Figure 2: IEEE34 node test feeder

In order to precisely characterize the operation of protective relays, we
first explain what the ideal relays are supposed to do using a simple
distribution circuit as given Fig. 1. There are 3 relays located at each bus
of the distribution line. Each relay is located to the right of a bus (node).
Each relay needs to protect its own region, which is between its own bus and
the first downstream bus. Each relay is also required to provide backup for
its downstream neighbor: when its neighbor fails to operate, it needs to trip
the line and clears the fault. For example, in Fig. 1, if a fault occurs at
the point where the lightning indicator is, relay C is the main relay
protecting this segment and it should trip the line immediately. If relay C
fails to work, relay B, which provides backup for relay C, needs to trip the
line instead, after a short delay. The time delay between the fault occurrence
and relay tripping should be as small as possible for primary relays, while
backup relays should react slower to ensure that they are triggered only when
the corresponding main relay is not working.

In feeder protection, for example in the IEEE 34 node feeder as shown in Fig.
2, a recloser is often placed at the substation (bus 800), while in some
circumstances they might also be placed within the feeder circuit. A recloser
usually needs to coordinate with fuses since the melting of fuses is
irreversible and preventing fuses from melting during transient faults is
preferable. This coordination is usually implemented using slow-fast curves,
such that the recloser attempts to clear the fault by quickly opening and
reclosing before the fuse melts, and if the fault is persistent, the fuse will
melt to clear the fault shortly after. If the fuse fails to melt, the recloser
will be locked open as backup protection for the fuse.

###  II-B Markov Decision Processes and Reinforcement Learning

Before formulating the relay protection problem using the RL approach, we
first give a brief review of the basic terminologies of RL.

_Markov Decision Processes_ (MDP) is a canonical formalism for stochastic
control problems. The goal is to solve sequential decision making (control)
problems in stochastic environments where the control actions can influence
the evolution of the state of the system. An MDP is modeled as tuple
$(\mathcal{S},\mathcal{A},R,P,\gamma)$ where $\mathcal{S}$ is the state space,
$\mathcal{A}$ is the action space.
$P=(P(\cdot|s,a),(s,a)\in\mathcal{S}\times\mathcal{A})$ is the state
tranistion probabilities. $P(s^{\prime}|s,a)$ specifies the probability of
transition to state $s^{\prime}$ upon taking action $a$ in state $s$.
$R:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the reward function,
and $\gamma\in[0,1)$ is the discount factor.

A policy $\pi:\mathcal{S}\rightarrow\mathcal{A}$ specifies the control action
to take in each possible state. The performance of a policy is measured using
the metric value of a policy, $V_{\pi}$, defined as

$\displaystyle
V_{\pi}(s)=\mathbb{E}[\sum^{\infty}_{t=0}\gamma^{t}R_{t}|s_{0}=s],$

where $R_{t}=R(s_{t},a_{t}),a_{t}=\pi(s_{t}),s_{t+1}\sim
P(\cdot|s_{t},a_{t})$. The optimal value function $V^{*}$ is defined as
$V^{*}(s)=\max_{\pi}~{}V_{\pi}(s)$. Given $V^{*}$, the optimal policy
$\pi^{*}$ can be calculated using the Bellman equation as

$\displaystyle\pi^{*}(s)=\arg\max_{a\in\mathcal{A}}~{}(R(s,a)+\gamma\sum_{s^{\prime}\in\mathcal{S}}P(s^{\prime}|s,a)V^{*}(s^{\prime})).$

Similar to the value function, Q-value function of a policy $\pi$, $Q_{\pi}$,
is defined as
$Q_{\pi}(s,a)=\mathbb{E}[\sum^{\infty}_{t=0}\gamma^{t}R_{t}|s_{0}=s,a_{0}=a]$.
Optimal Q-value function $Q^{*}$ is also defined similarly,
$Q^{*}(s,a)=\max_{\pi}Q_{\pi}(s,a)$. Optimal Q-value function will help us to
compute the optimal policy directly without using the Bellman equation, as
$\pi^{*}(s)=\arg\max_{a\in\mathcal{A}}~{}Q^{*}(s,a)$

Given an MDP formulation, the optimal value/Q-value function or the optimal
control policy can be computed using dynamic programming methods [18].
However, these methods require the knowledge of the full model of the system,
namely, the transition probability. In most real world applications, the
system model is either unknown or extremely difficult to estimate. In the
protective relay problem, the transition probability represents all the
possible stochastic variations in the voltages and currents in the network,
due to a large number of scenarios like weather (and the resulting shift in
demand/supply) and renewable energy generation. In such scenarios, the optimal
policy has to be learned from sequential state/reward observations.

Reinforcement learning is a method for learning the optimal policy for an MDP
when the model is unknown. RL achieves this without explicitly constructing an
empirical model. Q-learning is one of the most popular RL algorithms which
learn the optimal Q-value function from the sequence of observations
$\left(s_{t},a_{t},R_{t},s_{t+1}\right)$. Q-learning algorithm is implemented
as follows. At each time step $t$, the RL agent updates the Q-function $Q_{t}$
as

$\displaystyle Q_{t+1}(s_{t},a_{t})$
$\displaystyle=(1-\alpha_{t})Q_{t}(s_{t},a_{t})$ $\displaystyle\hskip
28.45274pt+\alpha_{t}(R_{t}+\gamma\max_{b}Q_{t}(s_{t+1},b))$ (1)

where $\alpha_{t}$ is the step size (learning rate). It is known that if each-
state action pairs is sampled infinitely often and under some suitable
conditions on the step size, $Q_{t}$ will converge to the optimal Q-function
$Q^{*}$ [18].

Using a standard Q-learning algorithm as described above is infeasible in
problems with continuous state/action space. To address this problem,
Q-function is typically approximated using a deep neural network, i.e.,
$Q(s,a)\approx Q_{w}(s,a)$ where $w$ is the parameter of the neural network.
Deep neural networks can approximate arbitrary functions without explicitly
designing the features. This has enabled tremendous success in both supervised
learning (image recognition, speech processing) and reinforcement learning
(AlphaGo games) tasks.

In Q-learning with neural network based approximation, the parameters of the
neural network can be updated using stochastic gradient descent with step size
$\alpha$ as

$\displaystyle w$ $\displaystyle=w+\alpha\nabla Q_{w}(s_{t},a_{t})$
$\displaystyle\hskip
28.45274pt(R_{t}+\gamma\max_{b}Q_{w}(s_{t+1},b)-Q_{w}(s_{t},a_{t}))$ (2)

Q-learning with neural network is further enhanced by using experience replay
and target network. Experience replay is to break the temporal correlation of
observations by randomly sampling some data points from a buffer of previously
observed (experienced) data points to perform the gradient step in (II-B).
Target network is used to overcome the instability of the gradient descent due
to the moving target. The combination of neural networks, experience replay
and target network forms the core of the _DQN algorithm_ [9]. In the
following, we will use DQN as one of the basic block for our proposed nested
RL algorithm.

###  II-C Protective Relay Control as an RL Problem

We formulate the distribution system transient process as an MDP environment
and model the relays as RL agents. Each relay can only observe its local
measurements of voltage $(s^{v}_{i,t})$ and current $(s^{c}_{i,t})$. Each
relay also knows the status of the local current breaker circuits, i.e., if it
is open or closed $(s^{b}_{i,t})$. Each relay also has a local counter that
ensures the necessary time delay in its operation as a backup relay
$(s^{d}_{i,t})$. These variables constitute the state
$s_{i,t}=(s^{v}_{i,t},s^{c}_{i,t},s^{b}_{i,t},s^{d}_{i,t})$ of each relay $i$
at time $t$. Table I summarizes this state space representation. Note that the
state includes the past $m$ measurements.

When a relay detects a fault, it will decide to trip. However, since each
relay is able to observe only its local state and no communication is possible
between the relays, some implicit coordination between the relays is
necessary. In traditional overcurrent protection scheme, the coordination is
achieved using an inverse-time curve that adds a time delay between the
detection of fault and actual breaker operation, based on the fault current
magnitude. However, the fault current magnitude can be unpredictable across
different scenarios, especially with DG and smart edge-devices. We propose
another approach (that is also amenable to RL) as follows. Instead of tripping
the breaker instantaneously, it controls a countdown timer to indirectly
operate the breaker. If a fault is detected, the relay can set the counter to
a value such that the breaker trip after a certain time delay. If the fault is
cleared by another relay during the countdown, the relay will reset the
counter to prevent mis-operation. Table II summarizes the action space of each
relay. The action of relay $i$ at time $t$, $a_{i,t}$, is one of these 11
possible values.

The reward given to each relay is determined by its current action and fault
status. A positive reward occurs if, i) it remains closed during normal
conditions, ii) it correctly operates after a fault in its assigned region or
in its first downstream region when the corresponding primary relay fails. A
negative reward is caused by, i) tripping when there is no fault; 2) tripping
after a fault outside its assigned region. The magnitude of the rewards are
designed in such a way to facilitate the learning, implicitly signifying
relative importance of false positives and false negatives. The reward
function for each relay is shown in Table III.

Consider a network with $n$ relays. Define the global state of the network at
time $t$ as $\bar{s}_{t}=(s_{1,t},s_{2,t},\ldots,s_{n,t})$ and the global
action at time $t$ as $\bar{a}_{t}=(a_{1,t},a_{2,t},\ldots,a_{n,t})$. Let
$R_{i,t}$ be the reward obtained by relay $i$ at time $t$. Define the global
reward $\bar{R}_{t}$ as $\bar{R}_{t}=\sum^{n}_{i=1}R_{i,t}$. The state of the
system evolves based on the load profile, DG output, presence of fault and the
connectivity of the circuit. Note that the (global) state evolution of the
network can longer be described by looking at the local transition
probabilities because the control actions of the relays affect states of
others relays. The global dynamics is represented by the transition
probability $\bar{P}(\bar{s}_{t+1}|\bar{s}_{t},\bar{a}_{t})$.

We formulate the optimal relay protection problem in a network as multi-agent
RL problem. The goal is to achieve a global objective, maximizing the
cumulative reward obtained by all relays, using only local control laws
$\pi_{i}$ which maps the local observations $s_{i,t}$ to local control action
$a_{i,t}$. Formally,

$\displaystyle\max_{(\pi_{i})^{n}_{i=1}}~{}\mathbb{E}[\sum^{\infty}_{t=0}\gamma^{t}\bar{R}_{t}],~{}a_{i,t}=\pi_{i}(s_{i,t}).$
(3)

Since the model is unknown and there is no communication between relays, each
relay has to learn its own local control policy $\pi_{i}$ using an RL
algorithm to solve (3).

TABLE I: Relay State Space State | Description
---|---
$s^{v}_{i,t}$ | Local voltage measurements of past $m$ timesteps
$s^{c}_{i,t}$ | Local current measurements of past $m$ timesteps
$s^{b}_{i,t}$ | Status of breaker (open (0) or closed (1))
$s^{d}_{i,t}$ | Value of the countdown timer
TABLE II: Relay Action Space Action | Description
---|---
$a_{\text{set}}$ | Set the counter to value to an integer between 1 and 9
$a_{d}$ | Decrease the value the counter by one
$a_{\text{reset}}$ | Stop and reset the counter
TABLE III: Reward for Different Operations Reward | Condition
---|---
+100 |  Tripping when a fault is present in its assigned protection region
-120 |  Tripping when there is no fault or the fault is outside its assigned region
+5 |  Stay closed when there is no fault or the fault is outside its assigned region
-10 |  Stay closed when a fault is present in its assigned protection region

##  III Nested Reinforcement Learning for Control of Protective Relays

Classical RL algorithms and their deep RL versions typically address only the
single agent learning problem. A multi-agent learning environment violates one
of the fundamental assumption needed for the convergence of RL algorithms,
namely, the stationarity of the operating environment. In a single agent
system, for any fixed policy of the learning agent, the probability
distribution of the observed states can be described using a stationary Markov
chain. RL algorithms are designed to learn only in such a stationary Markovian
environment. Multiple agents taking actions simultaneously violate this
assumption. Even if the policy of a given agent is fixed, state observations
for that agent are no longer according to a stationary Markov chain, as they
are controlled by the actions of other agents. Moreover, in our setting, each
relay observes only its local measurements which further complicates the
problem. There are existing literature [19] addressing this kind of problems,
but the performance of most algorithms are unstable and the convergence is
rarely guaranteed.

We propose an approach, which we call Nested Reinforcement Learning, to
overcome this difficulty of the multi-agent RL problem by exploiting the
radial structure of distribution systems. In radial distribution system, the
dependency between the operation of coordinating relays is uni-directional,
i.e., only upstream relays need to provide backup for a downstream relay but
not vice-versa. Also, the last relay at the load side does not need to
coordinate with others. In our nested RL algorithm, we start the RL training
from the last relay whose ideal operation is not affected by the operation of
other relays. So, it can be trained using a single-agent RL (DQN) algorithm.
Then, we can fix the trained policy for this last relay and train the relays
at one-level closer to the substation that need to provide backup for the last
relay. Since the policy of the furthest relay is fixed, it appears like a part
of the stationary environment to its upstream neighbors which can learn to
accommodate its operation. This process can be repeated for all the relays
upstream to the feeder. The order of training can be determined by network
tracing using a post-order depth-first tree traversal with the source bus
being the root. This nested training approach which exploits the nested
structure of the underlying physical system allows us to overcome the non
stationarity issues presented in generic multi-agent RL settings. Our nested
RL algorithm is formally presented below.

Algorithm 1 Nested Reinforcement Learning Algorithm

Initialize replay buffer of each relay $i$, $1\leq i\leq n$

Initialize DQN of each relay $i$ with random weights

for relay $i=1$ to $n$ do

for episode $k=1$ to $M$ do

Initialize simulation with random system parameters

for time step $t=1$ to $T$ do

Observe the state $s_{i,t}$ for all relays

for relay $j=1$ to $i$ (Trained Relays) do

Select action using the trained policy

$a_{j,t}=\arg\max_{a}Q_{w^{*}_{j}}(s_{j,t},a)$

end for

for relay $j=i+1$ to $n$ do

Select the null action, $a_{j,t}=0$

end for

With probability $\epsilon$ select a random action $a_{i,t}$, otherwise select
the greedy action $a_{i,t}=\arg\max_{a}Q_{w_{i}}(s_{i,t},a)$

Observe the reward $R_{i,t}$ and next state $s_{i,t+1}$

Store $(s_{i,t},a_{i,t},R_{i,t},s_{i,t+1})$ in the replay buffer of relay $i$

Sample a minibatch from replay buffer and update the DQN parameter $w_{i}$

end for

end for

end for

##  IV Experiment Environment and Test Cases

In this section, we describe the simulation environment, test system modelling
and experiment design. The simulation experiment is performed using OpenDSS
distribution simulator [20] and IEEE standard test feeders [21].

###  IV-A Simulation Environment

The dynamic simulation mode of OpenDSS is used to build the environment. A
Python program controls the simulation process by communicating with the
simulator using the COM interface. This environment is packed in a class
inherited from the OpenAI Gym [22] to improve accessibility. We note that this
setting can potentially be used in a number of other research problems
addressing the distribution systems operation using machine learning.

The RL algorithm is programmed in Python using open-source machine learning
packages Tensorflow [23] and Keras-RL [24]. The hyperparameters of the DQN for
each relay are selected through random search and are listed in Table IV.

TABLE IV: DQN Hyperparameters Hyperparameter | Value
---|---
Hidden Layers | 128/256/128
Activation | ReLU/ReLU/ReLU/Linear
Replay Buffer | 2000
Target Network Update Rate | 0.005
Double DQN | ON
Optimizer and Learning Rate | Adam, 0.0001

###  IV-B Test System Modeling

The standard test cases from the IEEE PAS AMPS DSAS Test Feeder Working Group
are adopted as the benchmark testbed for this experiment [21]. Specifically,
we choose the IEEE 34-bus and 123-bus test feeders to test the performance of
RL based recloser relay control. The test cases are replicated in OpenDSS
using the same parameters provided in IEEE publications [21]. Since OpenDSS
only support spot load and each load must be attached to a bus, distributed
loads in the 34-bus case are lumped together and attached to a dummy bus
created at the midpoint of the branch where distributed loads exist in the
original case. The voltage regulators are manually set to follow the published
values in the IEEE power flow solution. Overall, OpenDSS power flow result and
IEEE results agree closely, though some minor difference still exist due to
differences in component models and arithmetic precision between RDAP and
OpenDSS. The percentage difference of phase voltages at each bus of the
IEEE-34 base case between the OpenDSS simulations IEEE published values are
listed in Table V. The 123-bus case do not have distributed loads around
branches and the solution is very close to the IEEE published value, and hence
the comparison is omitted.

TABLE V: Difference Between OpenDSS and IEEE Solution % Error | $V_{a}$ | $V_{b}$ | $V_{c}$
---|---|---|---
Average | 0.179 | 0.240 | 0.023
Maximum | 0.637 | 0.554 | 0.066

Modifications to the IEEE cases are done when initializing each _episode_ to
simulate the real fluctuations of distribution grids. An episode is defined as
a short simulation segment that contains a fault. Fault scenarios are
generated to provide the agents a representation of possible cases in a real
distribution system. In the beginning of each episode, a random multiplier in
the range $(0.7,1.3)$ is set for all loads, as in distribution grids most
loads tend to peak around the same time of day. Then an individual random
multiplier between $(0.9,1.1)$ is applied for each load in the system.
Distributed generation are also placed throughout the feeder, a random number
of generators modeled as PV inverters are added to random load buses. The
capacity of each DG can account for 50% to 125% of the load at the same bus.
In the middle of an episode, a random fault is added to the system.The fault
will have a random fault impedance ranging from 0.001 ohm to 20 ohm. All types
of faults (SLG, LL, LLG, 3-phase) are considered. The fault scenarios are
generated using Monte-Carlo sampling. Single phase and two phase faults have a
higher chance to be selected as they are much more common and harder to detect
than symmetric faults. The performance of the trained RL relays are evaluated
by averaging over a number of episodes.

##  V Simulation Results

###  V-A Performance Metrics

In this section we present and discuss the performance of our Nested RL
algorithm for protective relays. We compare the performance with the
conventional overcurrent relay protection strategy. The performance is
evaluated in three aspects:

Failure Rate: A relay failure happens when a relay fails to operate as it is
supposed to do. For each episode, we determine the optimal relay action from
the type, time, and location of the fault, and compare it to the action taken
by the RL based relay. We evaluate the percentage of the operation failures of
the relays in four different scenarios: when there is a (i) fault in the local
region, (ii) fault in the immediate downstream region, (iii) fault in a remote
region, (iv) no fault in the network.

Robustness: Load profiles in a distribution system is affected by many events
like weather, social activities, renewable generation, and electric vehicles
charging schedules. These events can generally cause the peak load to
fluctuate and possibly exceed the expected range in the planning stage.
Moreover, the electricity consumption is expected to slowly increase each
year, reflecting the continuing economy and population growth. This can cause
a shift in the mean (and variance) of the load profile. Relay protection
control should be robust to such changes as continually reprogramming relays
after deployment is costly. We evaluate the performance of RL relays when the
operating condition exceeds the nominal range or under unexpected conditions
such as sudden loss of load/generator.

Response Time: The response time of RL relay is defined as the time difference
between fault occurrence and the relay activation. Response time is extremely
critical in preventing hazards. For example, it is preferred for the
substation recloser to attempt clearing transient faults before any fuse in
the feeder melts. This requires the recloser to have a fast fault detection
time. We compare the response time of the RL based relays with the
conventional overcurrent relays.

###  V-B Performance: Single Agent RL

We first present the performance of our RL algorithm for a single recloser
control. This is a special case of the proposed algorithm (with $n=1$). We
train and test our algorithm in the context of substation recloser control in
distribution feeders. In particular, we consider a recloser located at the
substation. The IEEE 34 bus and 123 bus feeders are used in this experiment.

The learning curves in Fig. 3(a) shows the convergence of episodic reward for
the substation recloser in the 123-bus system. The learning curve is obtained
from 20 independent training runs. The thick black line shows the mean value
of episodic reward and the shadow envelope indicates the mean value $\pm$
standard deviation. Fig. 3(b) shows the convergence of false operations during
training. The rate of false operation is high at the beginning of training,
and it drops to approximately zero after roughly 800 episodes. The learning
curves for the 34 bus system is similar and hence omitted.

(a) 123-bus Episodic reward

(b) 123-bus False Operation

Figure 3: Learning Curves: Single Relay

We have run the simulations with two types of input measurement for the RL
recloser relays: (i) phase voltage, current, and angle, and (ii) sequence
value of voltage and current. The two types of input have different
advantages. Sequence measurement is a more sensitive identifier of faults.
This is because during normal conditions the magnitude of zero and negative
sequence is expected to be very low compared with the magnitude of positive
sequence, while during unbalanced faults the magnitude of zero and negative
sequence will rise significantly. On the other hand, it would be better to use
separate phase measurements to identify which phase(s) are under fault as the
faulted phase will experience a larger voltage and current disturbance.
However, locating faulted phases using sequence measurement is less apparent.
We trained RL relays using these two set of measurements separately. When
using sequence data, the recloser gets a positive reward for responding to any
fault within the feeder. When using phase data, the recloser gets a positive
reward only for responding to faults that include the phase its assigned.
After training, each RL relay is tested using 5000 new randomly generated
scenarios and each false operation are recorded.

TABLE VI: Failure Rate Using Sequence Input IEEE 34 Node Feeder
---
Scenario | False Operation | Occurrences | Probability
No Fault | Trip | 0 / 5000 | 0.00 %
After Fault | Hold | 3 / 5000 | 0.06 %
IEEE 123 Node Feeder
Scenario | False Operation | Occurrences | Probability
No Fault | Trip | 0 / 5000 | 0.00 %
After Fault | Hold | 8 / 5000 | 0.16 %
TABLE VII: Failure Rate Using Phase Input IEEE 34 Node Feeder
---
Scenario | False Operation | Occurrences | Probability
No Fault | Trip | 0 / 5000 | 0.00 %
After Fault in Assigned Phase | Hold | 21 / 5000 | 0.42 %
After Fault in Other Phases | Trip | 28 / 5000 | 0.56 %
IEEE 123 Node Feeder
Scenario | False Operation | Occurrences | Probability
No Fault | Trip | 0 / 5000 | 0.00 %
After Fault in Assigned Phase | Hold | 55 / 5000 | 1.1 %
After Fault in Other Phases | Trip | 47 / 5000 | 0.94 %

Table VI summarizes the failure rate performance of RL relay in 34 bus and 123
bus test feeder using sequence value input. Table VII shows the performance
using phase value input. The accuracy of both formulation is much better than
traditional overcurrent protection in similar environment according to past
studies [25][26]. The performance of traditional overcurrent protection is
heavily affected by the presence of DG and fault impedance [2]. In the study
done in [26], the operation of overcurrent protection deteriorates evidently
when the DG penetration exceeds 20%. In contrast, the RL based relays are
extremely accurate even under very high DG penetration levels. We also note
that using sequence measurement as input has a even lower failure rate in
fault detection. Both DQN could be stored together as sequence values can be
used to detect faults under normal conditions and whenever a fault is
detected, the phase measurements can be fed into the other DQN to determine
which phase(s) are under fault.

To quantify the robustness of RL based algorithm against peak load variations,
we tested the performance by varying the peak load upto 20% more than the
maximum load used during training. Since we are considering the robustness
w.r.t. to the peak load variations, the load capacity used in this test is
sampled only from peak load under consideration. For example, the data
collected for 5% increase are sampled by setting the system load size between
100% and 105% of the peak load at training. Note that the we test the
performance using the same parameter from the original training, i.e., we
don’t update the policy to accommodate the change in this load profile. The
performance of RL relay under increases peak load is shown in Table VIII.

TABLE VIII: Robustness Against Peak Load Increase Peak Load Increase | 5 % | 10 % | 15% | 20 %
---|---|---|---|---
Failure Rate | 0.1 % | 0.1 % | 0.4 % | 1.2 %

We also evaluate the robustness against other disturbances. Particularly, if a
load or a generator is suddenly disconnected from the feeder, a disturbance in
the measurement waveform could be observed. It is important for the recloser
relay to distinguish these disturbances from faults. We test this aspect by
suddenly stepping down the capacity of all loads by a factor from 10% to 40%
at a random time during an episode, or randomly disconnecting 10% to 40% of
all the distributed generators in the feeder. The recloser is expected to stay
closed during these disturbances. The results obtained from 5000 test episodes
are given in Table IX. RL recloser has shown to be very robust against such
disturbances. As a rough comparison in [25], in the 34 bus system, the DG
change will incur coordination failure in about 10% to 30% of scenarios
depending on fault location and DG penetration.

TABLE IX: Robustness Against Loss of Load/DG Scenario | Occurrences | Probability
---|---|---
Loss of Load | 0 / 5000 | 0 %
Loss of DG | 28 / 5000 | 0.56 %
TABLE X: Response Speed After Faults Delay | 1 Step | 2 Step | 3 Step | 4 Steps
---|---|---|---|---
Occurance | 0 / 5000 | 4973 / 5000 | 23 / 5000 | 4 / 5000

We also measure the response time during the tests, quantified in terms of the
number of simulation steps where each simulation step is 0.002 second. The RL
relays have shown a very small response time as listed in Table X, the longest
delay is 4 simulation steps which corresponds to 8 ms. Moreover, the response
time is not correlated with fault current magnitude, and is much faster than
the melting time curve of typical time-delay fuses. We note that, in practice
however, the response time could be limited by the sampling rate of instrument
transformers.

###  V-C Performance: Multi-Agent RL

Our nested RL algorithm makes use of the radial structure of distribution
grids. By this approach, if a relay need to provide backup for a downstream
neighbor, it learns the optimal time delay before tripping the breaker for
each possible fault scenario during training to accommodate the policy of its
neighbor. We have conducted a proof-of-concept study for this idea in a simple
5-bus radial system in our previous work [16]. In the experiment below, we
expand this to a more realistic case using the IEEE 34 and 123 bus feeder test
case to demonstrate the coordination between RL relays.

In distribution systems, it is common for large and long feeders to have
additional reclosers in the middle of the feeder for additional security. For
persisting faults that cannot be cleared by reclosing, the recloser needs to
be locked open. In such cases, it is preferred for the closest protection
device to operate to reduce the amount of load being disconnected and mitigate
the damage. In the IEEE 34 bus case (Fig. 1), a mid-feeder recloser is added
to the branch between bus 828 and bus 830. For the 123 node case, the mid-
feeder recloser is placed in the branch between bus 160 and 67, protecting the
zone in the right half of the feeder. The mid-feeder recloser is expected to
act immediately for all faults in the right half of the circuit, and remain
closed for all faults between it and the substation. The substation recloser
needs to provide backup for the mid-feeder recloser, taking a longer time
delay for all faults pass the mid-feeder recloser and trip quickly for all
faults between the substation and the mid-feeder recloser.

For simulating the scenarios when a backup is needed, the mid-feeder recloser
has a 50% probability to be deactivated whenever it tries to trip. Thus, the
substation recloser can only get a reward for tripping for faults after the
mid-feeder recloser has attempted to trip. If the substation recloser trips
before the mid-feeder recloser attempted to trip, it will receive a penalty
instead. According to our nested RL algorithm, the mid-feeder recloser is
trained first with the substation recloser de-activated, and then the
substation recloser is trained with the mid-feeder recloser put into action.
Similar to the results from the single-agent experiment above, we present the
learning curves for the convergence of episodic reward of both recloser in
Fig. 4 shows the learning curve of the mid-feeder recloser at bus 828. The
learning curve of the substation recloser is also similar. The training for
both agent converged quickly after around 400 episodes.

(a) 34-bus Episodic reward Figure 4: Learning Curve: Multi-Agent RL

The failure rate of the recloser pair is measured based on the action of both
relays. An episode is considered successful only if both recloser take the
correct control actions. The operation is tested in 5000 random episodes and
the result is summarized in Table XI.

TABLE XI: Failure Rate of Multi-Agent RL Scenario | False Operation | Occurrences | Probability
---|---|---|---
IEEE 34 Node Feeder
No Fault | Trip | 0 / 5000 | 0.00 %
Local Fault | Hold | 4 / 5000 | 0.08 %
Remote Fault | Trip | 16 / 5000 | 0.32 %
Backup | Hold | 11 / 5000 | 0.22 %
IEEE 123 Node Feeder
No Fault | Trip | 0 / 5000 | 0.00 %
Local Fault | Hold | 9 / 5000 | 0.18 %
Remote Fault | Trip | 27 / 5000 | 0.54 %
Backup | Hold | 10 / 5000 | 0.2 %

Robustness against increased peak and unexpected disturbances are conducted
for the two-recloser pair similar to the single recloser scenario. A mis-
operation of one recloser is recorded as failure for the entire episode. The
results are listed in Table XII and XIII. The impact of peak shift is slightly
more evident than in previous single-relay cases due to the need for
coordination and the performance of RL relays starts to deteriorate at around
15% increased peak.

TABLE XII: Robustness Against Peak Increase: Multi-Agent Peak Load | 5 % | 10 % | 15 % | 20 %
---|---|---|---|---
Failure Rate | 1.6 % | 1.7 % | 3.5 % | 5.2 %
TABLE XIII: Robustness Against Loss of Load/DG: Multi-Agent scenario | occurrence | probability
---|---|---
Loss of Load | 0 / 5000 | 0 %
Loss of DG | 37 / 5000 | 0.74 %

The response time for the both reclosers under different scenarios is recorded
in Table XIV. It can be seen that the substation recloser responds faster to
faults that are between the substation and the mid-feeder recloser. For faults
in the right half of the circuit, the substation recloser provides a time
window of roughly 3 time steps for the closer neighbor to operate first.

TABLE XIV: Response Time Delay | 1 Step | 2 Step | 3 Step | 4+ Steps
---|---|---|---|---
Mid-feeder recloser | 0/5000 | 4973/5000 | 23/5000 | 4/5000
Delay | 3- Step | 4 Step | 5 Step | 6+ Steps
Substation recloser | 2910/5000 | 345/5000 | 1591/5000 | 154/5000

##  VI Concluding Remarks

This paper introduces and thoroughly tests a deep reinforcement learning based
protective relay control strategy for the distribution grid with many DERs. It
is shown that the proposed algorithm builds upon existing hardware yields much
faster and more consistent performance. This algorithm can be easily applied
in both a standalone relay and a network of coordinating relays. The trained
RL relays can accurately detect faults under situations including high fault
impedance, presence of distributed generation and volatile load profile, where
the performance of traditional overcurrent protection deteriorates heavily.
The RL relays are robust against unexpected changes in operating conditions of
the distribution grid at the time of planning, eliminating the need to re-
train the relays after deployments. The response speed of RL relays are very
fast, providing ample time for coordinating with fuses and other relays.

The proposed deep RL relays are easy to implement with the currently available
distribution infrastructure. A particularly attractive feature is that the
proposed algorithm for relays can operate in a completely decentralized manner
without any communication. This communication-free setting is not only easy to
implement for currently available distribution grid infrastructure, but also
less vulnerable to potential cyber-attacks. The input to the RL relays are the
same as traditional relays so the instrument transformers can be retained
during deployment. Moreover, since all computationally-expensive training is
done offline in a simulation environment, the computing power requirement for
the relay processor is relatively low. The weights of the DQN obtained during
training can be saved into a general-purpose micro-controller or potentially a
more optimized machine learning chip.

In the future, we plan to provide a theoretical guarantee for the convergence
of our sequential RL algorithm. We will work with EMTP simulators for more
detailed time-domain training data generation with realistic power electronic
and electromagnetic transient models. We will also investigate the possibility
of hardware prototyping and Hardware-in-the-Loop test with Real-Time Digital
Simulator (RTDS) in much larger systems.

## References

  * [1] U.S. Dept. of Energy, _The potential benefits of distributed generation and the rate related issues that may impede its expansion_ , 2007.
  * [2] U. Shahzad, S. Kahrobaee, S. Asgarpoor, _Protection of distributed generation: challenges and solutions_ , Energy and Power Engineering, vol.9, pp. 614-653, 2007.
  * [3] P. Dash, S. Samantaray and G. Panda, _Fault classification and section identification of an advanced series-compensated transmission line using support vector machine_ , IEEE transactions on power delivery, vol. 22, no. 1, pp. 67-73, 2007.
  * [4] H. Zhan et al., _Relay protection coordination integrated optimal placement and sizing of distributed generation sources in distribution networks_ , IEEE Transactions on Smart Grid, vol. 7, no. 1, pp. 55-65, 2016.
  * [5] H.-T. Yang, W.-Y. Chang, and C.-L. Huang, _A new neural networks approach to on-line fault section estimation using information of protective relays and circuit breakers_ , IEEE Transactions on Power Delivery, vol. 9, no. 1, pp. 220-230, 1994.
  * [6] X. Zheng et al., _A SVM based setting of protection relays in distribution systems_ , 2018 IEEE Texas Power and Energy Conference (TPEC), College Station, TX, USA, pp. 1-6, 2018.
  * [7] Y. Zhang, M. D. Ilić and O. Tonguz, _Application of Support Vector Machine Classification to Enhanced Protection Relay Logic in Electric Power Grids_ , LESCOPE, 2007.
  * [8] H. C. Kilic¸kiran, B. Kekezoglu and G. N. Paterakis, _Reinforcement Learning for Optimal Protection Coordination_ , 2018 International Conference on Smart Energy Systems and Technologies (SEST), Sevilla, pp. 1-6, 2018.
  * [9] V. Minh et al., _Human-level control through deep reinforcement learning_ , Nature, 518.7540:529, 2015.
  * [10] K. Arulkumaran, M. P. Deisenroth, M. Brundage and A. A. Bharath, _Deep reinforcement learning: a brief survey_ , IEEE Signal Processing Magazine, vol. 34, no. 6, pp. 26-38, Nov. 2017.
  * [11] H. Xu, A. Dominguez-Garcia, V. Veeravalli and P. W. Sauer, _Data-driven voltage regulation in radial power distribution systems_ , IEEE Transactions on Power Systems, Oct. 2019.
  * [12] J. Sun et al., _An integrated critic-actor neural network for reinforcement learning with application of DERs control in grid frequency regulation_ , International Journal of Electrical Power & Energy Systems, vol. 111, pp. 286-299, 2019.
  * [13] Q. H. Wu and J. Guo, _Optimal bidding strategies in electricity markets using reinforcement learning_ , Electric Power Components and Systems, vol. 32, pp. 175-192, Jun 2010.
  * [14] M. Bagheri et al., _Enhancing power quality in microgrids with a new online control strategy for DSTATCOM using reinforcement learning algorithm_ , IEEE access, vol. 6, pp. 38986-38996, 2018.
  * [15] T.P. Imthias Ahmed, P.S. Nagendra Rao and P.S. Sastry, _A reinforcement learning approach to automatic generation control_ , Electric Power Systems Research, vol. 63, pp. 9-26, Aug. 2002.
  * [16] D. Wu, X. Zheng, D. Kalathil and L. Xie, _Nested reinforcement learning based control for protective relays in power distribution systems_ , IEEE Conference on Decision and Control (CDC), Nice, France, Dec. 2019.
  * [17] M. Glavic, _(Deep) Reinforcement learning for electric power system control and related problems: A short review and perspectives_ , Annual Reviews in Control, vol. 48, pp. 22-35, 2019.
  * [18] R. S. Sutton, A. G. Barto, _Reinforcement learning: an introduction_ , MIT Press, 2nd edition, 2018.
  * [19] S. Kapoor, _Multi-Agent Reinforcement Learning: A Report on Challenges and Approaches_ , arXiv:1807.09427
  * [20] R. C. Dugan and T. E. McDermott, _An open source platform for collaborating on smart grid research_ , 2011 IEEE Power and Energy Society General Meeting, Detroit, MI, USA, pp. 1-7, 2017.
  * [21] K. P. Schneider et al., _Analytic considerations and design basis for the IEEE distribution test feeders_ , IEEE Transactions on Power Systems, vol. 33, pp. 3181-3188, May 2018.
  * [22] G. Brockman et al., _OpenAI gym_ , arXiv:1606.01540, 2016.
  * [23] M. Abadi et al., _TensorFlow: Large-scale machine learning on heterogeneous systems_ , Software available from tensorflow.org, 2015.
  * [24] M. Plappert, _Keras-rl_ , GitHub, accessed from: https://github.com/keras-rl/keras-rl, 2016.
  * [25] J. Silva, H. Funmilayo and K. Butler-Purry, _Impact of distributed generation on the IEEE 34 node radial test feeder with overcurrent protection_ , 89th North American Power Symposium (NAPS), 2007.
  * [26] A. F. Naiem, Y. Hegazy, A. Y. Abdelaziz and A. Elsharkawy, _A novel protection methodology for distribution systems equipped with distributed generation_ , International Electrical Engineering Journal Vol. 6 No.10 pp. 2048-2057, 2015.