# Sum-of-squares bounds on correlation functions in a minimal model of
turbulence

Vladimir Parfenyev<EMAIL_ADDRESS>Landau Institute for Theoretical
Physics, 142432 Chernogolovka, Russia  National Research University Higher
School of Economics, Faculty of Physics, 101000 Moscow, Russia  Evgeny
Mogilevskiy<EMAIL_ADDRESS>Lomonosov Moscow State University,
Faculty of Mechanics and Mathematics, 11992 Moscow, Russia  Weizmann Institute
of Science, 76100 Rehovot, Israel  Gregory Falkovich
<EMAIL_ADDRESS>Weizmann Institute of Science, 76100 Rehovot,
Israel  Landau Institute for Theoretical Physics, 142432 Chernogolovka, Russia

###### Abstract

We suggest a new computer-assisted approach to the development of turbulence
theory. It allows one to impose lower and upper bounds on correlation
functions using sum-of-squares polynomials. We demonstrate it on the minimal
cascade model of two resonantly interacting modes, when one is pumped and the
other dissipates. We show how to present correlation functions of interest as
part of a sum-of-squares polynomial using the stationarity of the statistics.
That allows us to find how the moments of the mode amplitudes depend on the
degree of non-equilibrium (analog of the Reynolds number), which reveals some
properties of marginal statistical distributions. By combining scaling
dependence with the results of direct numerical simulations, we obtain the
probability densities of both modes in a highly intermittent inverse cascade.
We also show that the relative phase between modes tends to $\pi/2$ and
$-\pi/2$ in the direct and inverse cascades as the Reynolds number tends to
infinity, and derive bounds on the phase variance. Our approach combines
computer-aided analytical proofs with a numerical algorithm applied to high-
degree polynomials.

##  I Introduction

Many systems in nature receive and dissipate energy on very different scales,
having conservative dynamics in between, and the energy is transferred across
the scales through a turbulent cascade. Among the best-known examples are
isotropic fluid turbulence [1] and surface waves in the ocean [2]. The
complexity of such systems makes them difficult to describe in detail, so it
makes sense to consider simpler dynamical models to deepen our understanding
of the statistical properties and energy transfer in such highly non-
equilibrium systems [3, 4, 5].

A minimal model, which still captures the basic properties of turbulent
cascades, is a system of two resonantly interacting oscillators whose natural
frequencies differ by a factor of two [6]. This system allows studying direct
and inverse cascades with the energy flux directed towards either higher or
lower frequencies. The system has one non-dimensional governing parameter
$\chi$, which plays the role of the Reynolds number. In what follows, we are
mostly interested in the regime when this parameter is large, and the
probability distribution tends to be singular. The analytical studies in this
limit resulted in the steady-state probability density for the direct cascade,
while constructing the probability density for the inverse cascade turns out
to be tricky [6].

Here we apply a complementary approach to study the mode statistics in this
system. We exploit the polynomial nature of dynamic equations (common for
practically all turbulent systems), which allows us to impose inequalities on
the correlation functions. The main idea is to find a non-negative polynomial
expression, $\phi(\bm{x})-L+F(\bm{x})\geq 0$ that combines the correlation
function $\phi(\bm{x})$ to be bounded, the constant value of the bound $L$,
and an auxiliary polynomial function $F(\bm{x})$ that has zero mean value
$\langle F(\bm{x})\rangle=0$ in a statistically steady state, which entails
the inequality $\langle\phi(\bm{x})\rangle\geq L$ [7, 8]. The essence of the
approach is to construct the function $F(\bm{x})$ so that the value of the
lower bound $L$ is as large as possible. Although testing a polynomial
expression $\phi(\bm{x})-L+F(\bm{x})$ for non-negativity is NP-hard
algorithmic task, there exist numerical procedures [9, 10, 11, 12, 13] that
solve the problem with a more strict requirement on the proposed expression to
be sum-of-squares (SoS) of other polynomials. Moreover, these procedures allow
one to maximize the value $L$ of the bound using some ansatz for the auxiliary
function $F(\bm{x})$. If the ansatz is simple enough, the bound can be found
analytically, while a computer algorithm can be used in advance to suggest the
optimal form of the SoS polynomial. The upper bound
$\langle\phi(\bm{x})\rangle\leq U$ can be constructed in a similar way. Recent
examples of the application of SoS programming to study dynamical systems can
be found in Refs. [7, 8, 14, 15, 16, 17].

The rest of the paper is organized as follows. In Section II, we remind the
general theory behind SoS optimization and then apply the method for the two-
mode system. In Section III, we present the results obtained by SoS
programming for the direct cascade and compare them with the analytics and
direct numerical simulations (DNS). It turns out that the upper and lower
bounds for the moments of pumped and dissipating modes are close to each other
and differ by less than a percent in the limit of a large Reynolds number.
This allows us to determine not only their scaling dependence but also
numerical values with accuracy comparable to DNS. The method can also be used
to study correlations between modes, and we show that the relative phase
between the modes tends to $\pi/2$ and determine upper and lower bounds for
its root-mean-square fluctuations.

After the validation, in Section IV, we use the method for the inverse cascade
where the probability density is unknown a priori. Analytically, we obtained
only lower bounds for the correlation functions, but they demonstrate scaling
consistent with the results of DNS. A numerical algorithmic analysis of the
high moments of the dissipating mode leads to scaling $\langle
n_{1}^{k}\rangle\propto\chi^{k-1}$, which is a fingerprint of intermittency,
where $n_{1}$ is the mode intensity and $\chi$ is the Reynolds number. Based
on this observation, we were able to shed light on the structure of the
distribution function of this mode – the DNS results fall on the universal
curve for different values of $\chi$, which has a form close to a power law
with an exponential cutoff. As for the pumped mode, its statistics are close
to Gaussian, and we analytically found the lower bound for its intensity
$\langle n_{2}\rangle\geq 5\chi/16$ that is close to DNS. We also show that
the relative phase between modes tends to $-\pi/2$ in the limit $\chi\gg 1$
and estimate the rate of this transition. Finally, we summarize and discuss
our findings in Section V.

##  II Theoretical Framework

This section briefly explains how sum-of-squares optimization can impose
inequalities on correlation functions in stochastic systems with polynomial
dynamics. The presentation follows Ref. [8], where a more detailed discussion
can be found.

Let us consider a stochastic dynamical system

$\dot{x_{i}}=f_{i}(\bm{x})+\sigma_{ij}(\bm{x})\xi_{j}(t),\quad\bm{x}\in\mathbb{R}^{n},\;\bm{\xi}\in\mathbb{R}^{m},$
(1)

where $f_{i}(\bm{x})$ and $\sigma_{ij}(\bm{x})$ are polynomial, $\xi_{i}(t)$
is a Gaussian noise with zero mean $\langle\xi_{i}(t)\rangle=0$ and the
variance
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$,
and here and below we sum over repeated indices. We assume that the system has
reached a statistical steady-state, and then the probability density function
$\rho(\bm{x})$ satisfies the stationary Fokker-Planck equation

$\frac{1}{2}\partial_{i}[\sigma_{ij}\partial_{k}(\sigma_{kj}\rho)]-\partial_{i}(f_{i}\rho)=0.$
(2)

We also assume that the solution to this equation is unknown. We wish to prove
a constant lower bound $\langle\phi(\bm{x})\rangle\geq L$ for some polynomial
correlation function $\phi(\bm{x})$, where the angle brackets mean the
averaging over the probability density $\rho(\bm{x})$.

For this purpose, we consider an auxiliary function $Q(\bm{x})$, which does
not grow fast when $|\bm{x}|\to\infty$, so that all the moments considered
below are finite. Performing integration by parts and neglecting boundary
terms, one can show that
$\langle\frac{1}{2}\sigma_{kj}\partial_{k}(\sigma_{ij}\partial_{i}Q)+f_{i}\partial_{i}Q\rangle=0$.
The idea is to properly design $Q(\bm{x})$ so that

$\left\langle\frac{1}{2}\sigma_{kj}\partial_{k}(\sigma_{ij}\partial_{i}Q)+f_{i}\partial_{i}Q+\phi-L\right\rangle\geq
0\ .$ (3)

Evaluation of the expectation in expression (3) requires knowing
$\rho(\bm{x})$, but it is sufficient for the inequality to hold point-wise for
all $\bm{x}$.

Checking a polynomial expression for non-negativity is NP-hard algorithmic
task, therefore to reduce computational complexity, it can be replaced with a
semidefinite programming (SDP) problem of checking the stronger condition that
expression (3) belongs to a set of polynomials that are sum-of-squares (SoS).
Finally, to find the maximum value of the constant $L$, we arrive at the
following optimization problem

$\max_{Q(\bm{x})}L:\;\frac{1}{2}\sigma_{kj}\partial_{k}(\sigma_{ij}\partial_{i}Q)+f_{i}\partial_{i}Q+\phi-L\in
SoS.$ (4)

One specifies an ansatz for the auxiliary function $Q(\bm{x})$ with
undetermined coefficients and then solves the problem using the software
packages such as SOSTOOLS [11], YALMIP [12] or SumOfSquares.jl [13] with one
of the appropriate SDP-solvers [18, 19, 20, 21]. The upper bounds
$\langle\phi(\bm{x})\rangle\leq U$ can be obtained in a similar way by
considering the optimization problem

$\min_{Q(\bm{x})}U:\;-\frac{1}{2}\sigma_{kj}\partial_{k}(\sigma_{ij}\partial_{i}Q)-f_{i}\partial_{i}Q-\phi+U\in
SoS.$ (5)

If the ansatz is simple enough, the bounds $L$ and $U$ can be found
analytically, and the computer algorithm suggests the optimal form of the SoS
polynomials. The more complex ansatz leads to more rigorous bounds; however,
complex ansatz requires more computational efforts, and the numerical
algorithm could fail to find an optimal solution.

##  III Direct cascade

Figure 1:  Average values of mode intensities (a), their second (b), and
higher moments (c) for the direct cascade. Solid and dashed lines show
numerical ($d=10$) and analytical ($d=4$) results for the upper (red) and
lower (blue) bounds. The square markers present results of DNS, and dash-
dotted lines correspond to asymptotics following from relation (8).
Expressions for the asymptotic dependence are shown in the panels.

The direct cascade for the two-mode system is governed by the following non-
dimensional system of equations [6]

$\displaystyle\displaystyle\dot{b}_{1}=-\frac{2ib_{1}^{*}b_{2}}{\sqrt{\chi}}+\xi(t),$
(6)
$\displaystyle\displaystyle\dot{b}_{2}=-\frac{ib_{1}^{2}}{{\sqrt{\chi}}}-b_{2},$
(7)

where $b_{1}$ and $b_{2}$ are complex envelopes of low- and high-frequency
modes, respectively, $\xi(t)$ is a complex Gaussian random force with zero
mean $\langle\xi(t)\rangle=0$ and the variance
$\langle\xi(t)\xi^{*}(t^{\prime})\rangle=\delta(t-t^{\prime})$, and the
asterisk means complex conjugate. In other words, real and imaginary parts of
$\xi(t)$ are independent white noises having zero mean values and intensities
$1/2$. The parameter $\chi$ is the only control parameter in the system, and
it quantifies the dissipation relative to interaction strength. In the
statistical steady-state, the energy input rate equals the energy flux from
the first mode to the second and the dissipation rate. For the amplitude of
the dissipating mode and the flux, one finds $\langle
n_{2}\rangle\equiv\langle|b_{2}|^{2}\rangle=1/4$, $\langle
J\rangle\equiv\langle ib_{1}^{2}b_{2}^{*}+c.c.\rangle/\sqrt{\chi}=-1/2$.

In the limit $\chi\ll 1$, the interaction between modes is strong, and the
energy transfer is fast, so one may expect that the occupation numbers are
close to energy equipartition $\langle|b_{1}|^{2}\rangle\equiv\langle
n_{1}\rangle=2\langle n_{2}\rangle$ corresponding to thermal equilibrium [6].
The DNS shows that the marginal distributions of the mode amplitudes are
indeed close to the Gaussian statistics with the corresponding occupation
numbers. Still, the phase between the modes $\theta=\arg(b_{1}^{2}b_{2}^{*})$
is unevenly distributed, which indicates the presence of correlations between
the modes [6].

In the opposite limit of weak interaction and strong noise, $\chi\gg 1$, one
expects that the driven mode needs much higher amplitude to provide for the
flux: $\langle n_{1}\rangle\gg\langle n_{2}\rangle$. In other words, the mode
statistics are far from thermal equipartition. In this case, the asymptotic
analytical solution for the whole distribution function was argued to be
singular [6]:

$\mathcal{P}(b_{1},b_{2})=\frac{2^{3/2}}{\pi^{3/2}\chi^{1/2}}\exp\left(-\frac{2}{\chi}|b_{1}|^{4}\right)\delta\left(b_{2}+\frac{ib_{1}^{2}}{\sqrt{\chi}}\right)\
.$ (8)

We will use this result later for comparison with the results of our approach,
which provides mutual validation.

Let us now apply the SoS optimization method described in the previous
section. To proceed, we need to specify an ansatz for the auxiliary function
$Q$. We have found that the numerical procedure works better if $Q$ is a
polynomial with respect to $n_{1},\ n_{2},\ J$ and takes the monomials
$n_{1}^{i}n_{2}^{j}J^{k}$ with the total degree of mode amplitudes less than
$d\geq 2i+2j+3k$. A more general ansatz does not improve estimates for the
correlation functions, but the algorithm is less stable due to the expansion
of the optimization space. Let us emphasize that for small values of $d$, the
bounds for the correlation functions can be obtained analytically, the
computer algorithm operates for higher values of $d$ and improves the result.

We begin to present our results with the intensity of the pumped mode. For
$d=4$ we obtain analytically [22]

$\frac{\chi^{1/2}}{2\sqrt{2}}\leq\langle
n_{1}\rangle\leq\dfrac{1}{2}\left(3-\sqrt{3}+\sqrt{12-6\sqrt{3}+\chi}\right),$
(9)

and in the limit $\chi\gg 1$, this implies
${\chi^{1/2}}/{2\sqrt{2}}\leq\langle n_{1}\rangle\leq\sqrt{\chi}/2$.
Inequality (9) means, in particular, that the intensity of the pumped mode is
much greater than the intensity of the dissipated mode, and $\langle
n_{1}\rangle/\langle n_{2}\rangle\propto\sqrt{\chi}$. As the parameter $d$
increases, the numerically found upper and lower bounds for $d=10$ approach
very closely the asymptotic $\langle n_{1}\rangle=\sqrt{\chi/2\pi}$ following
from expression (8), see Fig. 1a. In the opposite case $\chi\ll 1$, the upper
bound $\langle n_{1}\rangle\leq 3-\sqrt{3}$ does not depend on $\chi$ in
qualitative agreement with DNS. The lower bound is not tight, and increasing
the parameter $d$ does not qualitatively change its behavior.

Similar results are also obtained for higher moments of $n_{1}$. For the
second moment and $d=4$ we analytically find [22]

$\chi^{2}r_{l}(\chi)\leq\langle
n_{1}^{2}\rangle\leq\frac{\chi}{4}+\chi^{2}r_{u}(\chi),$ (10)

where $r_{l}(\chi)$ is the largest real root of the equation
$1-\chi^{2}r+8\chi^{3}r^{2}-16\chi^{4}r^{3}=0$ and $r_{u}(\chi)$ is the
smallest positive real root of the equation
$1+16\chi+(4\chi+152\chi^{2})r+(348\chi^{3}-32\chi^{4})r^{2}-216\chi^{5}r^{3}+16\chi^{7}r^{4}=0$.
In the limit $\chi\gg 1$, the upper and lower bounds coincide with each other
and therefore $\langle n_{1}^{2}\rangle\to\chi/4$, while in the opposite case
$\chi\ll 1$, one obtains $2^{-4/3}\chi^{2/3}\leq\langle
n_{1}^{2}\rangle\lesssim 1.88$, where the upper bound reflects scaling
consistent with DNS, see Fig. 1b. For even higher moments of $n_{1}$, the
values in the limit $\chi\gg 1$ can be determined numerically with good
accuracy since the upper and lower bounds tend to each other as $d$ increases.
Analyzing the values of the moments, one can conclude that the statistics of
the mode $b_{1}$ are essentially non-Gaussian in agreement with expression
(8), see Fig. 1c. This kind of analysis can help one to guess the marginal
distribution function if it is not known a priori.

Next, we turn to the dissipating mode. Its mean intensity is determined
exactly from the energy balance condition, $\langle n_{2}\rangle=1/4$, and the
upper and lower bounds for $\langle n_{2}^{2}\rangle$ are shown in Fig. 1b. In
the limit $\chi\gg 1$, the reasonable bounds are obtained for relatively large
values of $d\geq 6$, so we do not present analytical results. In the opposite
case $\chi\ll 1$, the bounds demonstrate scaling consistent with DNS, but they
are not close to each other. Analytically we obtain $1/16\leq\langle
n_{2}^{2}\rangle\leq 3/16$ for $d=4$ and $\chi\ll 1$ [22]. The analysis of
higher moments for $\chi\gg 1$ is presented in Fig. 1c, and the results are in
agreement with the marginal distribution
$P(b_{2})=\frac{2^{1/2}}{\pi^{3/2}|b_{2}|}e^{-2|b_{2}|^{2}}$ following from
expression (8).

Figure 2:  Relative phase between modes for the direct cascade: numerically
obtained upper (red) and lower (blue) bounds (circles, triangles, and diamonds
are for $d=10,\ 12,\ 14)$; the square markers present results of DNS.

The SoS programming method can also be used to analyze correlations between
modes. To estimate the relative phase between modes, we rewrite the initial
equations (6)-(7) in terms of real variables $\rho_{1}=|b_{1}|$,
$\rho_{2}=|b_{2}|$, $\theta=\arg(b_{1}^{2}b_{2}^{*})$:

$\displaystyle\dot{\rho_{1}}=-\dfrac{2\rho_{1}\rho_{2}\sin\theta}{\sqrt{\chi}}+\dfrac{1}{4\rho_{1}}+\dfrac{\zeta_{1}(t)}{\sqrt{2}},$
(11)
$\displaystyle\dot{\rho_{2}}=\dfrac{\rho_{1}^{2}\sin\theta}{\sqrt{\chi}}-\rho_{2},$
(12)
$\displaystyle\dot{\theta}=\dfrac{\rho_{1}^{2}-4\rho_{2}^{2}}{\rho_{2}\sqrt{\chi}}\cos\theta+\dfrac{\sqrt{2}\zeta_{2}(t)}{\rho_{1}},$
(13)

where the overall phase drops out, and $\zeta_{i}(t)$ is a real Gaussian noise
with zero mean and the pair correlation function
$\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$
[22]. Despite the right-hand sides of these equations are not polynomial, for
any functions $Q$ and $\phi$ that are polynomial with respect to $\rho_{1}$,
$\rho_{2}$, $\sin\theta$ the expressions in left-hand sides of Eqs. (4), (5)
are polynomials divided by $\rho_{1}^{2}\rho_{2}>0$. Thus, to apply the
algorithm, it is sufficient to multiply the expression in the optimization
problems (4), (5) by $\rho_{1}^{2}\rho_{2}>0$ to return it into the class of
polynomial functions, see details in Ref. [22]. Now the ansatz for the
function $Q$ is a polynomial of $\rho_{1},\ \rho_{2},\ \sin\theta$ of the
power of $d$. We have also found that the numerical procedure gives better
results if we solve the optimization problem with the additional constraint
$1-\sin^{2}\theta\geq 0$ [22].

Figure 3:  Average values of mode intensities (a), their second moments (b),
and the relative phase (c) for the inverse cascade. Solid and dashed lines
show numerical ($d=10$) and analytical results for upper and lower bounds. The
square markers present results of DNS, and dash-dotted lines correspond to
asymptotic bounds in the limit $\chi\gg 1$. Expressions for the analytic
results are shown in the panels.

Fig. 2 shows the upper and lower bounds on $\langle\cos^{2}\theta\rangle$. The
value of $\langle\cos^{2}\theta\rangle\to 0$ as $\chi\to\infty$, which means
that the phase $\theta\to\pi/2$ (the point $\theta=-\pi/2$ is not suitable
because $\langle
J\rangle\equiv\langle-2\rho_{1}^{2}\rho_{2}\sin\theta\rangle/\sqrt{\chi}=-1/2<0$).
The power-law fit $\langle\cos^{2}\theta\rangle\propto\chi^{-q}$ results in
$q=1/3$ and $q=2/5$ for the upper and lower bounds, respectively. We could not
determine the bounds analytically since the algorithm gives reasonable
estimates only for large values of $d$. The overall picture is that the
distribution function of the phase between the modes has a peak at
$\theta=\pi/2$, and its width decreases in a power-law manner with the
parameter $\chi$. These results complement Ref. [6], in which the narrowing of
the phase distribution was not quantified.

##  IV Inverse cascade

Now we turn to the inverse cascade where the high-frequency mode is pumped,
and the low-frequency mode dissipates

$\displaystyle\displaystyle\dot{b}_{1}=-\frac{2ib_{1}^{*}b_{2}}{\sqrt{\chi}}-b_{1},$
(14)
$\displaystyle\displaystyle\dot{b}_{2}=-\frac{ib_{1}^{2}}{{\sqrt{\chi}}}+\xi(t).$
(15)

In the statistical steady-state, from the energy balance consideration, we
find $\langle n_{1}\rangle\equiv\langle|b_{1}|^{2}\rangle=1$ and $\langle
J\rangle\equiv\langle ib_{1}^{2}b_{2}^{*}+c.c.\rangle/\sqrt{\chi}=1$, which is
true for any value of $\chi$. The energy flux $\langle J\rangle$ is positive
that corresponds to the transfer of energy down the frequencies. As in the
case of a direct cascade, in the limit $\chi\ll 1$, one can expect the
occupation numbers to be close to the energy equipartition $\langle
n_{1}\rangle=2\langle n_{2}\rangle$, although the statistics is not expected
to be close to the product of two Gaussians corresponding to thermal
equilibrium. In the opposite case $\chi\gg 1$, the system is far from thermal
equilibrium, and the pumped mode is expected to have larger intensity,
$\langle n_{2}\rangle\gg\langle n_{1}\rangle$. All these expectations were
confirmed by DNS [6].

We found that we can get better bounds for correlation functions if we rewrite
dynamic equations (14)-(15) in terms of real variables $\rho_{1}=|b_{1}|$,
$\rho_{2}=|b_{2}|$, $\theta=\arg(b_{1}^{2}b_{2}^{*})$:

$\displaystyle\dot{\rho_{1}}=-\dfrac{2\rho_{1}\rho_{2}\sin\theta}{\sqrt{\chi}}-\rho_{1},$
(16)
$\displaystyle\dot{\rho_{2}}=\dfrac{\rho_{1}^{2}\sin\theta}{\sqrt{\chi}}+\dfrac{1}{4\rho_{2}}+\dfrac{\zeta_{1}(t)}{\sqrt{2}},$
(17)
$\displaystyle\dot{\theta}=\dfrac{\rho_{1}^{2}-4\rho_{2}^{2}}{\rho_{2}\sqrt{\chi}}\cos\theta+\dfrac{\zeta_{2}(t)}{\sqrt{2}\rho_{2}},$
(18)

where the overall phase drops out [22]. We also found that in contrast to the
direct cascade, it is useful to extend the ansatz for the function $Q$ by
adding the term $\log\rho_{1}$. The term $\log\rho_{2}$ does not have a
noticeable effect on the results. The left-hand sides of Eqs. (4), (5) should
be multiplied by $\rho_{2}^{2}>0$ so that they become polynomial and we can
apply the SDP algorithm to solve the optimization problem [22].

In contrast to the direct cascade, we were able to obtain only lower bounds
for the correlation functions in the case of the inverse cascade. The
algorithm does not find a feasible solution for the upper bounds, even for a
relatively large parameter value $d=10$. Fortunately, the lower bounds are
quite informative precisely in the turbulent limit of large $\chi$ which we
focus on. In particular, for the intensity of the pumped mode, we analytically
find [22]

$\langle n_{2}\rangle\geq\dfrac{\chi+1}{4},$ (19)

and both asymptotics at $\chi\ll 1$ and $\chi\gg 1$ give a scaling that agrees
qualitatively with DNS, see Fig. 3a. Physically, inequality (19) means that in
the limit $\chi\gg 1$, the intensity of the pumped mode is much greater than
the intensity of the dissipating mode, $\langle n_{2}\rangle/\langle
n_{1}\rangle\geq\chi/4$. We thus have shown that deviation from the
equipartition is much stronger in the inverse cascade than in the direct one
(where the ratio is $\propto\sqrt{\chi}$).

Similarly, we find for the amplitude fourth moment [22]

$\langle n_{2}^{2}\rangle\geq\dfrac{(\chi+1)^{2}}{16},\quad\langle
n_{1}^{2}\rangle\geq\dfrac{\chi+4}{3}.$ (20)

The first condition is trivial and follows from the positive variance of the
pumped mode intensity. The second condition means that in the limit $\chi\gg
1$, the ratio $\langle n_{1}^{2}\rangle/\langle n_{1}\rangle^{2}\geq\chi/3$,
i.e., the statistics of the dissipated mode is intermittent. This agrees with
the analysis carried out in Ref. [6], where it was shown that the $\rho_{1}$
dynamics is a sequence of burst events, see also Fig. 4a. Between bursts, the
amplitude $\rho_{1}$ is close to zero, and during short bursts with a duration
of order unity, $\rho_{1}$ reaches large values $\sim\sqrt{\chi}$. The time
interval between bursts is $\sim\chi$, and the correlation functions of
$\rho_{1}$ saturate on bursts. In Fig. 3b, we compare the obtained
inequalities with DNS. All asymptotics demonstrate correct scaling with the
parameter $\chi$. Note that analytical inequalities are improved by the
numerical algorithm when we increase the parameter $d$.

Fig. 3c shows the upper bound for $\langle\cos^{2}\theta\rangle\leq U$,which
is equivalent to the lower bound for $\langle\sin^{2}\theta\rangle\geq 1-U$.
In the limit $\chi\gg 1$, the value of $\langle\cos^{2}\theta\rangle\to 0$ and
it means that the the phase $\theta\to-\pi/2$, since the value of $\langle
J\rangle=\langle-2\rho_{1}^{2}\rho_{2}\sin\theta\rangle/\sqrt{\chi}=1$. The
power-law fit $\langle\cos^{2}\theta\rangle\propto\chi^{-q}$ results in
$q=1/2$, although the interval is short and with increasing $\chi$ we have to
increase $d$ so that the algorithm finds a feasible solution. Compared to the
direct cascade, the phase fluctuations in the inverse cascade are smaller for
the same value of $\chi$. Based on this observation, we can simplify equations
(16)-(18) by assuming that the relative phase is locked on $\theta=-\pi/2$.
Then, we obtain

$\displaystyle\dot{\rho_{1}}=\dfrac{2\rho_{1}\rho_{2}}{\sqrt{\chi}}-\rho_{1},$
(21)
$\displaystyle\dot{\rho_{2}}=-\dfrac{\rho_{1}^{2}}{\sqrt{\chi}}+\dfrac{1}{4\rho_{2}}+\dfrac{\zeta(t)}{\sqrt{2}},$
(22)

and these equations can be used to improve the above bounds in the limit
$\chi\gg 1$. In particular, for the intensity of the pumped mode, we
analytically obtain $\langle n_{2}\rangle\geq 5\chi/16$ [22], which is closer
to DNS results, see Fig. 3a.

Figure 4: (a) DNS fragment of the mode dynamics for the inverse cascade
demonstrate intermittency of the dissipated mode, $\log_{10}\chi=2.4$; (b)
Probability density for the pumped mode in the inverse cascade. The dashed
line corresponds to the exponential distribution, $P(x)=e^{-x}$; (c) Reduced
probability density $F(n_{1}/\chi)$ for the dissipating mode in the inverse
cascade, see Eq. (23). The inset shows the same data in the log-log scale, and
the dashed line corresponds to the numerical fit $F(x)=1.7x^{-1}e^{-1.7x}$.
The bin size is $0.05$ for the main panel and $0.001$ for the inset.

The reduction in the number of degrees of freedom allows us to use larger
values of parameter $d$ in the limit $\chi\gg 1$ since the size of the
optimization space is also reduced. Moreover, we were able to obtain upper
bounds ($d=32$) for correlation functions $\langle n_{2}\rangle$ and $\langle
n_{1}^{2}\rangle$, which are close to the lower bounds, see Figs. 3a and 3b.
To estimate the higher moments of both modes, we will further use the values
of lower bounds obtained for $d=32$. In this way, for the pumped mode, we
found the usual scaling $\langle n_{2}^{k}\rangle\propto\langle
n_{2}\rangle^{k}$ [22]. In the intervals between bursts, the amplitude
$\rho_{1}$ is close to zero, and, according to Eq. (15), one can expect that
the statistics of the pumped mode is close to Gaussian. The value $\langle
n_{2}\rangle$ can be estimated as a diffusion displacement during the time
between bursts, $\langle n_{2}\rangle\sim\chi$, in agreement with the results
reported earlier. The DNS data for the probability density function confirm
this qualitative analysis, although the agreement is not perfect, see Fig. 4b.

For the dissipating mode, we found that $\langle
n_{1}^{k}\rangle\propto\chi^{k-1}$ [22], and such scaling is a fingerprint of
intermittence discussed above. This suggests that in the limit $\chi\gg 1$,
the distribution function of the dissipating mode has the following form:

$P(n_{1}/\chi)=\dfrac{1}{\chi}F(n_{1}/\chi).$ (23)

DNS confirms this hypothesis and Fig. 4c shows the reduced probability density
function $F(x)$, which is proportional to $1/x$ for $x\ll 1$ and has an
exponential cutoff for $x\gtrsim 1$.

A simplified model for the dynamics of the dissipating mode explains the
exponent $-1$ of the power law and resolves the singularity at $x\to 0$.
During the burst, the amplitude $\rho_{2}$ of the pumped mode quickly
diminishes (see Fig. 4a), therefore the second term in Eq. (14) dominates and
the subsequent dynamics of the dissipating mode corresponds to the exponential
decay, $\dot{x}=-2x$. For a given $x$, the probability density is determined
by the time spent in the vicinity of $x$ and thus proportional to $F(x)\propto
1/|\dot{x}|\propto 1/x$. The exponential dynamics is valid until the amplitude
of the pumped mode is below $\sim\sqrt{\chi}$. In this regime, the pumped mode
grows due to diffusion and reaches the level of $\sqrt{\chi}$ in time of the
order of $\chi$, so the dependence $F(x)\propto 1/x$ holds for $x\gtrsim
e^{-2\chi}$. This value should be taken as the lower limit of integration in
the expression for the total probability to resolve the singularity. Since the
minimum value of $x$ depends on the parameter $\chi$, the height of the first
bin on the histograms is proportional to $\chi$. Except for this, the shape of
the curve $F(n_{1}/\chi)$ is universal, see Fig. 4c.

When calculating the positive moments of $x$, the singularity disappears and
the lower limit of integration can be replaced by zero. From the condition
$\langle n_{1}\rangle=1$ we obtain $\int_{0}^{\infty}xF(x)dx=1$, and so the
function has the form $F(x)\simeq ax^{-1}e^{-ax}$. Numerical fitting leads to
$a\approx 1.7$, and the corresponding curve is shown in Fig. 4c by a dashed
line.

##  V Conclusion

So what have we learned about the statistics of the far from the equilibrium
state of the system using SoS programming? We obtained computer-aided
analytical and numerical bounds for correlation functions in the system of two
interacting modes in the regimes corresponding to both direct and inverse
energy cascades. The bounds revealed the scaling of mode intensities and their
higher moments with the Reynolds number $\chi$, which shows how far the system
is from the energy equipartition corresponding to the thermal equilibrium. By
combining scaling dependence with DNS, we collapsed the probability densities
of dissipating and pumped modes on the universal curves in a highly
intermittent regime of the inverse energy cascade and determined their shapes.
Analyzing cross-mode correlations, we showed that the relative phase between
modes tends to $\pi/2$ and $-\pi/2$ in the direct and inverse energy cascades,
respectively, and estimated the rates of these transitions.

The method applied does not work well in the opposite limit of low Reynolds
number $\chi\ll 1$. This can be explained by the fact that the term produced
by the white noise forcing in the optimization problems (4) has a relatively
small amplitude and, therefore, the expression that is analyzed practically
coincides with the case of a deterministic problem, which has a trivial
solution $b_{1}=b_{2}=0$, and for this reason the lower bound does not exist.
Methods for overcoming this issue are discussed in Ref. [8], but we did not
use them since the regime $\chi\ll 1$ is less interesting for us from the
physical point of view.

Let us also note that the direct and inverse cascades are very different,
although the equations at first glance may seem similar. In the direct
cascade, the dissipating mode is excited by the pumped mode in an additive
way. For the inverse cascade, this process is multiplicative, and the general
experience suggests that this regime is more complicated because of higher
temporal intermittency. Our analysis supports this intuition: in the inverse
cascade we were able to obtain only lower bounds for the correlation
functions, while for the direct cascade, we obtained both lower and upper
bounds, and they are so close that we can determine the numerical values of
the correlation functions with an accuracy comparable to DNS.

It would be promising to extend the present study to other shell models of
turbulence with longer chains of resonantly interacting modes [5, 23]. The
main difficulty is that for the long chains the extreme modes in the
correlation function of interest will be coupled to the modes adjacent to
them, and therefore the dynamic equations for the subset of modes involved in
the correlation function are not closed. A similar problem arose when applying
the sum-of-squares method to the Galerkin expansion of the Navier-Stokes
equation in Ref. [7]. Although the generalization is not straightforward, the
development of a complementary approach for this analysis significantly
advances this area of research.

###### Acknowledgements.

The work was supported by the Excellence Center at WIS, grant 662962 of the
Simons Foundation, grant 075-15-2022-1099 of the Russian Ministry of Science
and Higher Educations, grant 823937 and 873028 of the EU Horizon 2020
programme, the BSF grants 2018033 and 2020765. V.P. is grateful to the
Weizmann Institute, where most of the work was done, for their hospitality.
E.M. thanks Benoziyo Endowment Fund for the Advancement of Science for funding
his visit to the Weizmann Institute of Science. DNS were performed on the
cluster of the Landau Institute.

## References

  * Frisch and Kolmogorov [1995] U. Frisch and A. N. Kolmogorov,  _Turbulence: the legacy of AN Kolmogorov_ (Cambridge university press, 1995).
  * Zakharov  _et al._ [2012] V. E. Zakharov, V. S. L’vov, and G. Falkovich,  _Kolmogorov spectra of turbulence I: Wave turbulence_ (Springer Science & Business Media, 2012).
  * Biferale [2003] L. Biferale, Shell models of energy cascade in turbulence, Annual review of fluid mechanics 35, 441 (2003).
  * Obukhov [1969] A. Obukhov, Integral invariants in equations of the hydrodynamic type, in  _Dokl. Akad. Nauk SSSR_ , Vol. 184 (1969) pp. 2–3.
  * Vladimirova  _et al._ [2021a] N. Vladimirova, M. Shavit, and G. Falkovich, Fibonacci turbulence, Phys. Rev. X 11, 021063 (2021a).
  * Vladimirova  _et al._ [2021b] N. Vladimirova, M. Shavit, S. Belan, and G. Falkovich, Second-harmonic generation as a minimal model of turbulence, Physical Review E 104, 014129 (2021b).
  * Chernyshenko  _et al._ [2014] S. I. Chernyshenko, P. Goulart, D. Huang, and A. Papachristodoulou, Polynomial sum of squares in fluid dynamics: a review with a look ahead, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20130350 (2014).
  * Fantuzzi  _et al._ [2016] G. Fantuzzi, D. Goluskin, D. Huang, and S. I. Chernyshenko, Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization, SIAM Journal on Applied Dynamical Systems 15, 1962 (2016).
  * Parrilo [2000] P. A. Parrilo,  _Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization_ (California Institute of Technology, 2000).
  * Parrilo [2003] P. A. Parrilo, Semidefinite programming relaxations for semialgebraic problems, Mathematical programming 96, 293 (2003).
  * Prajna  _et al._ [2002] S. Prajna, A. Papachristodoulou, and P. A. Parrilo, Introducing sostools: A general purpose sum of squares programming solver, in  _Proceedings of the 41st IEEE Conference on Decision and Control, 2002._ , Vol. 1 (IEEE, 2002) pp. 741–746.
  * Lofberg [2004] J. Lofberg, Yalmip: A toolbox for modeling and optimization in matlab, in  _2004 IEEE international conference on robotics and automation (IEEE Cat. No. 04CH37508)_ (IEEE, 2004) pp. 284–289.
  * Legat  _et al._ [2017] B. Legat, C. Coey, R. Deits, J. Huchette, and A. Perry, Sum-of-squares optimization in Julia, in  _The First Annual JuMP-dev Workshop_ (2017).
  * Papachristodoulou and Prajna [2002] A. Papachristodoulou and S. Prajna, On the construction of lyapunov functions using the sum of squares decomposition, in  _Proceedings of the 41st IEEE Conference on Decision and Control, 2002._ , Vol. 3 (IEEE, 2002) pp. 3482–3487.
  * Papachristodoulou and Prajna [2005] A. Papachristodoulou and S. Prajna, A tutorial on sum of squares techniques for systems analysis, in  _Proceedings of the 2005, American Control Conference, 2005._ (IEEE, 2005) pp. 2686–2700.
  * Tan and Packard [2006] W. Tan and A. Packard, Stability region analysis using sum of squares programming, in  _2006 American Control Conference_ (IEEE, 2006) pp. 6–pp.
  * Goulart and Chernyshenko [2012] P. J. Goulart and S. Chernyshenko, Global stability analysis of fluid flows using sum-of-squares, Physica D: Nonlinear Phenomena 241, 692 (2012).
  * Andersen  _et al._ [2009] E. D. Andersen, B. Jensen, J. Jensen, R. Sandvik, and U. Worsøe, Mosek version 6, Technical Report TR-2009-3, MOSEK, Tech. Rep. (2009).
  * Fujisawa  _et al._ [2008] K. Fujisawa, M. Fukuda, K. Kobayashi, M. Kojima, K. Nakata, M. Nakata, and M. Yamashita, SDPA (SemiDefinite Programming Algorithm) and SDPA-GMP User’s Manual—version 7.1. 1, Research Reports on Mathematical and Computing Sciences, B-448 (2008).
  * Sturm [1999] J. F. Sturm, Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones, Optimization methods and software 11, 625 (1999).
  * Tütüncü  _et al._ [2003] R. H. Tütüncü, K.-C. Toh, and M. J. Todd, Solving semidefinite-quadratic-linear programs using SDPT3, Mathematical programming 95, 189 (2003).
  * [22] See Supplemental Material at [URL will be inserted by publisher] for technical details on analytical and numerical computations of bounds for the correlation functions.
  * Shavit  _et al._ [2022] M. Shavit, N. Vladimirova, and G. Falkovich, Emerging scale invariance in a model of turbulence of vortices and waves, Philosophical Transactions of the Royal Society A 380, 20210080 (2022).


# Advanced Safety Filter for Smooth Transient Operation of a Battery Energy
Storage System

Michael Schneeberger, Florian Dörfler, Silvia Mastellone Michael Schneeberger
and Silvia Mastellone are with the Institute of Electrical Engineering, FHNW,
Windisch, Switzerland (e-mails<EMAIL_ADDRESS>silvia.mastellone@fhnw.ch), and Florian Dörfler is with the Department of
Information Technology and Electrical Engineering, ETH Zürich, Switzerland,
(e-mail: dorfler@control.ee.ethz.ch)This work was supported by the Swiss
National Science Foundation (SNSF) under NCCR Automation.

###### Abstract

In this paper, we implement an advanced safety filter to smoothly limit the
current of an inverter-based Battery Energy Storage System. The task involves
finding suitable Control Barrier Function and Control Lyapunov Function via
Sum-of-Squares optimization to certify the system’s safety during grid
transients. In contrast to the conventional safety filter, the advanced safety
filter not only provides a safety certificate but also achieves finite-time
convergence to a nominal region. Within this region, the action of the nominal
control, i.e. the Enhanced Direct Power Control, remains unaltered by the
safety filter. The advanced safety filter is implemented using a Quadratically
Constrained Quadratic Program, providing the capability to also encode
quadratic input constraints. Finally, we showcase the effectiveness of the
implementation through simulations involving a load step at the Point of
Common Coupling, and we compare the outcomes with those obtained using a
standard vector current controller.

##  I Introduction

A Battery Energy Storage System (BESS) enables part of the power grid to
disconnect from the utility grid and operate independently in an islanded
mode. In this scenario, the primary objective of the BESS is to maintain grid
voltage and frequency stability through the use of an inert grid-forming (GFM)
control scheme. However, when operating in grid-connected mode, the BESS’s
current needs to be limited during transient grid events such as faults, load
steps or changes in configuration to protect the system’s hardware. In such
situations, a fast response is necessary to ensure the system’s safety.

In a microgrid configuration, a compelling GFM control scheme is the Enhanced
Direct Power Control (EDPC) [1, 2]. One version of this approach, described in
[3], consists of an active and reactive power controller that, together with a
feed-forward term derived from a filtered voltage measurement at the Point of
Common Coupling (PCC), determines the converter voltage reference. In
islanding configuration, the EDPC operates in GFM mode due to the slow
dynamics of this filtered PCC voltage. However, during transients, faster
filtered PCC voltage dynamics are temporarily permitted, even at the expense
of potentially sacrificing the GFM property, to facilitate a stable return to
EDPC operation.

To limit the converter current in grid-connected mode, a vector current
controller can be activated based on a grid fault detection [4]. However, this
requires manually tuning (assuming a worst-case fault) the control gains of
the current controller to respect the current limits and often causes
undesired chattering behaviors resulting from the activation and deactivation
of the current controller. Current limiting strategies based on limiting the
current reference, or introducing virtual impedance have been proposed [5, 6].
In this approach, unlike the EDPC, a vector current control is an integral
part of the GFM control scheme, rendering any switching action obsolete. None
of these approaches provide a safety certification for the BESS across all
states.

In this paper, we employ a _safety filter_ approach to limit the current and
achieve a smooth transition between operation under EDPC and current limiting
control. The safety filter is obtained by solving a Quadratic Program (QP) at
every state based on a Control Barrier Function (CBF) and an internal system
model to ensure forward invariance of the _safe set_ while tracking a nominal
control reference, as described in [7]. The safety filter can be integrated
with a QP encoding the Control Lyapunov Function (CLF) condition, ensuring not
only safety but also convergence to a specific state or set [8]. Polynomial
CBF and CLF candidates, encoding the specifications of the safety filter, can
be derived numerically by solving a series of Sum-of-Squares (SOS)
optimization problems [9, 10]. Additionally, linear input constraints can be
encoded as SOS constraints [11, 12].

In our previous work, we proposed an advanced safety filter [13], that extends
the basic version by ensuring a finite-time convergence to a forward invariant
set, referred to as _nominal region_ , wherein the nominal controller remains
undisturbed by the safety filter. In this work we apply this concept to the
BESS control scenario, where the safety filter is limiting the current while
tracking the EDPC’s voltage reference. To ensure a finite-time convergence
guarantee to the nominal region under stationary but unknown grid conditions,
we additionally incorporate the dynamics of the filtered PCC voltage into the
internal system model of the safety filter. Utilizing the SOS-based
alternating algorithm outlined in [13], suitable CBF and CLF candidates are
identified that adhere to the quadratic input constraints of the BESS.
Subsequently, the advanced safety filter is implemented using a Quadratically
Constrained Quadratic Program (QCQP), providing the capability to encode the
quadratic input constraints. Finally, we demonstrate the effectiveness of the
proposed solution through simulations where a load step is applied at the
Point of Common Coupling (PCC). The system response is compared with the one
obtained using a standard vector current controller. This study demonstrates
the applicability and efficacy of the safety filter concept in a non-trivial
power electronics example beyond its traditional application in robotics.

The remainder of the paper is structured as follows: Section II describes the
problem setup including the BESS, the EDPC, and the requirement on the current
limiting control. In Section III, we present the concept of the advanced
safety filter, its properties, and the resulting SOS formulation. Numerical
results, presented in Section IV, demonstrate the superior performance of the
safety filter in comparison to a conventional vector current control [4].
Finally, Section V concludes the paper.

##  II Battery Energy Storage System

###  II-A Problem setup

Consider a BESS implemented as a three-phase power inverter system connected
to the PCC through a transformer, as depicted in Figure 1. Given that the
state of charge (SOC) of the BESS’s battery changes at a much slower rate than
the control, we model it as a constant dc-voltage source $v_{dc}$. A grid
filter installed between transformer and PCC, along with an optimized pulse
pattern (OPP), suppresses high-frequency harmonics generated by the inverter
modulation. This justifies the adoption of an averaged inverter model, where
the BESS’s converter voltage is expressed as the product of the dc-link
voltage $v_{dc}\in\mathbb{R}$ and the modulation index $m\in\mathbb{R}^{2}$:

$\displaystyle v_{c}=v_{dc}m,$

where $m$ is constrained by $m_{c}^{\mathrm{max}}\in\mathbb{R}$:

$\displaystyle
m_{c,d}^{2}+m_{c,q}^{2}\leq\left(m_{c}^{\mathrm{max}}\right)^{2}.$ (1)

Given that the controller’s bandwidth is much smaller than the frequency band
of the grid filter, we simplify the model of the grid filter by representing
it as an inductance. Similarly, the transformer, which steps up the voltage to
match the grid voltage at the PCC, is represented by an inductance. The
current dynamics over the filter and transformer inductance $l_{c}$ in p.u.
are expressed as

$\displaystyle\frac{l_{c}}{\omega_{n}}\begin{bmatrix}\dot{i}_{d}\\\
\dot{i}_{q}\end{bmatrix}=-Z_{c}\underbrace{\begin{bmatrix}i_{d}\\\
i_{q}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar
58\relax}}i}-\underbrace{\begin{bmatrix}v_{\mathrm{PCC},d}\\\
v_{\mathrm{PCC},q}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar
58\relax}}v_{\mathrm{PCC}}}+\underbrace{\begin{bmatrix}v_{c,d}\\\
v_{c,q}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar 58\relax}}v_{c}},$ (2)

where $\omega_{n}\mathrel{\mathop{\mathchar 58\relax}}=2\pi 50$Hz is the
nominal grid frequency, the current $i$ and the voltage $v_{\mathrm{PCC}}$ are
measured, as indicated in Figure 1, and $Z_{c}\mathrel{\mathop{\mathchar
58\relax}}=J\omega l_{c}$ is the matrix representation of the complex filter
and transformer inductance with $J=\begin{bmatrix}0&-1\\\ 1&0\end{bmatrix}$.
The precise value of $l_{c}$ is determined during the commissioning of the
BESS.

Figure 1:  The BESS is implemented as a three-phase power inverter system
connected to the PCC via a transformer and a grid filter. The EDPC in series
with a current limiting control ensures GFM behavior during islanded grid
operation while limiting the current during grid transients.

Operating the BESS in both grid-connected and islanding modes, requires
suitable control strategies that meet grid stability and load requirements,
while respecting the constraints and preserve the safety of the battery and
converter components. The BESS control architecture consist of two main
components. The slow outer loop, in form of the EDPC, is responsible for
meeting the grid stability and load requirements, while the fast inner loop,
in form of a current limiting controller, is designed to limit the current
during transients. Although both control loops are designed in continuous-
time, they are implemented in discrete-time with sufficiently fast task cycles
to preserve the controllers performance.

###  II-B Enhanced Direct Power Control

The EDPC, taken from [1, 2, 3] and illustrated in Figure 2, consists of an
active and reactive power loop, achieving a slow time constant through the
utilization of a low-pass filter (LPF) in the measurement of active power $p$
and reactive power $q$. The active and reactive power references are
determined by a frequency droop $D_{f}\in\mathbb{R}$ and a voltage droop
$D_{v}\in\mathbb{R}$:

$\displaystyle p_{r}$ $\displaystyle=D_{f}(\omega-1)$ $\displaystyle q_{r}$
$\displaystyle=D_{v}(\hat{v}_{\mathrm{PCC},lp}-1),$

where $\omega$ is the angular frequency of the PLL in p.u., and
$\hat{v}_{\mathrm{PCC},lp}$ is the low-pass filtered PCC voltage amplitude.
The active and reactive power control, depicted in Figure 2, are based on the
PI controller

$\displaystyle K_{p}\left(1+\frac{1}{T_{i}s}\right),$

characterized by a proportional gain $K_{p}\in\mathbb{R}$ and integral time
constant $T_{i}\in\mathbb{R}$. This controller transforms the active and
reactive power error into a (nominally unlimited) current reference. A current
reference limiter
$\|i_{r}\|^{2}=i_{r,d}^{2}+i_{r,q}^{2}\leq(i_{r}^{\mathrm{lim}})^{2}$ for some
current limit $i_{r}^{\mathrm{lim}}<i^{\mathrm{max}}$ ensures a stable steady-
state current below the maximum allowed current 6. A cross-coupling term
$Z_{c}$ computes the voltage drop over the transformer and grid filter
impedance. Finally, a filtered PCC voltage $v_{\mathrm{PCC},f}$ is added to
the voltage drop, producing the BESS’s voltage reference

$\displaystyle v_{c,n}(i_{r},v_{\mathrm{PCC},f})\mathrel{\mathop{\mathchar
58\relax}}=Z_{c}i_{r}+v_{\mathrm{PCC},f}.$ (3)

The dynamical system characterizing the filtered PCC voltage is given as

$\displaystyle\underbrace{\begin{bmatrix}\dot{v}_{\mathrm{PCC},f,d}\\\
\dot{v}_{\mathrm{PCC},f,q}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar
58\relax}}\dot{v}_{\mathrm{PCC},f}}=\underbrace{\begin{bmatrix}\alpha_{d}\\\
\alpha_{q}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar 58\relax}}\alpha},$ (4)

where $\alpha\in\mathbb{R}^{2}$ is the input of this system. To achieve GFM
behavior, the filtered PCC voltage must exhibit characteristics of a low-pass
filter with time constant $\tau\approx 100ms$, as described by

$\displaystyle\alpha_{n}(v_{\mathrm{PCC}},v_{\mathrm{PCC},f})\mathrel{\mathop{\mathchar
58\relax}}=\frac{v_{\mathrm{PCC}}-v_{\mathrm{PCC},f}}{\tau}.$ (5)

Only during transients, it is allowed to change the behavior of the filtered
PCC voltage.

Figure 2:  The Enhanced Direct Power Control (EDPC) directly generates a
converter voltage reference $v_{c,n}$ from a power control loop, integrating
frequency and amplitude droop characteristics.

###  II-C The problem of current limiting control

The current limiting control at the output of the EDPC, as depicted in Figure
1, ensures that the current remains within the maximum allowed boundaries,
characterized by $i^{\mathrm{max}}\in\mathbb{R}$:

$\displaystyle i_{d}^{2}+i_{q}^{2}\leq(i^{\mathrm{max}})^{2}.$ (6)

This is achieved by adjusting the EDPC’s voltage reference $v_{c,n}$ 3 to
generate a safe voltage reference $v_{c}$. Furthermore, the dynamics governing
the filtered PCC voltage $v_{\mathrm{PCC},f}$ in 4 and 5 can be temporarily
adjusted during a current transient, facilitating a stable return to EDPC
operation.

###### Remark 1

During islanding mode, the BESS is assumed to be sufficiently dimensioned to
withstand transients, without having to resort to a fast current controller.
This assumption is crucial, as the stability of an islanded grid consisting
purely of grid-following (GFL) converters cannot be assured otherwise. Hence,
the current limiting controller is solely necessary when operating within an
interconnected grid, permitting the BESS to briefly adopt a fast GFL mode
without compromising grid stability.

The design of the current limiting control and later of the advanced safety
filter requires knowledge of the grid model, which in practice can be of
arbitrarily complexity and subject to change over time. Moreover, the original
models used for the EDPC have an unrequired degree of complexity. Therefore,
we propose the following simplifying assumptions:

  * •

Given that the discrete-time implementation of the safety filter utilizes
updated $v_{\mathrm{PCC}}$ measurements within each task cycle, we adopt a
simplified reduced-order grid model based on the current dynamics 2, where
$v_{\mathrm{PCC}}$ is assumed to be constant in dq reference frame.

  * •

Given that the EDPC operates with a slower time constant of approximately
100ms compared to the current limiting control’s time constant of around 1ms,
we propose a simplified reduced-order EDPC model based on 3, where the current
reference $i_{r}$ is assumed to be constant.

Under those assumptions, we can formulate the problem statement as follows:

###### Problem 1

Given the dynamical system described in 2 and 4, with a constant voltage
source $v_{\mathrm{PCC}}$, a set of allowable states 6, an input constraint
set 9, and a nominal controller described in 3 and 5, with a constant current
reference $i_{r}$, the task is to design an advanced safety filter, as
described in [13], that ensures stable and safe behavior of the BESS during
grid events.

A standard method to solve Problem 1 is to employ a vector current control
[4]. This involves a PI controller integrated with an anti-windup mechanism to
track a safe current reference within the dq reference frame. However, this
leads to chattering effects, abrupt behavior, and possibly system shutdown in
case of protection function activation.

##  III Advanced safety filter

To overcome the limitations of the state of the art vector current controller,
we propose an alternative solution to the problem of overcurrent phenomena
happening during transients. Our approach allows to promptly respond to grid
events while guaranteeing a smooth transitioning between EDPC and current
limiting control operation. This entails designing an advanced safety filter
setup, as shown in Figure 3, to enforce the state and input constraints while
ensuring finite-time convergence to the EDPC operation.

Figure 3:  The advanced safety filter ensures safe operation with respect to
the maximum allowed current limits 6. Additionally, it guarantees finite-time
convergence to the nominal region, where the nominal control input
$u_{n}^{\top}=\begin{bmatrix}v_{c,n}^{\top}&\alpha_{n}^{\top}\end{bmatrix}$ is
– assuming simplified system and EDPC control dynamics – not perturbed by the
advanced safety filter. The dynamics of the filtered PCC voltage
$v_{\mathrm{PCC},f}$ is also taken into account to ensure the existence of the
nominal region.

Consider the control system composed of the simplified grid model derived from
2, the filtered PCC voltage dynamics derived from 4, and the stationary
dynamics of the current reference $i_{r}$ and PCC voltage $v_{\mathrm{PCC}}$
as follows:

$\displaystyle\dot{x}=\underbrace{\begin{bmatrix}-(\omega_{n}/l_{c})\left(Z_{c}i+v_{\mathrm{PCC}}\right)\\\
\mathrm{0}\\\ \mathrm{0}\\\
\mathrm{0}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar
58\relax}}f(x)}+\underbrace{\begin{bmatrix}(\omega_{n}/l_{c})\mathrm{I}_{2}&\mathrm{0}\\\
\mathrm{0}&\mathrm{I}_{2}\\\ \mathrm{0}&\mathrm{0}\\\
\mathrm{0}&\mathrm{0}\end{bmatrix}}_{=\mathrel{\mathop{\mathchar
58\relax}}G(x)}\underbrace{\begin{bmatrix}v_{c}\\\
\alpha\end{bmatrix}}_{=\mathrel{\mathop{\mathchar 58\relax}}u},$ (7)

where $x^{\top}\mathrel{\mathop{\mathchar
58\relax}}=\begin{bmatrix}i^{\top}&v_{\mathrm{PCC},f}^{\top}&i_{r}^{\top}&v_{\mathrm{PCC}}^{\top}\end{bmatrix}$
is the state of the system. The reduced-order model of the inverter current
$i$ and the filtered PCC voltage $v_{\mathrm{PCC},f}$ are derived from the
specifications outlined in Section II-C. The current reference $i_{r}$ and PCC
voltage $v_{\mathrm{PCC}}$, although considered constant, are included in the
system dynamics to render the safety filter parametric in these variables. The
set of allowed states is defined by the maximum allowed current:

$\displaystyle\mathcal{X}_{a}=\left\\{x\mid\text{subject to
\lx@cref{creftype~refnum}{eq:allowable_current}}\right\\}.$ (8)

The input constraint set resulting from 1 is given as:

$\displaystyle\mathcal{U}=\left\\{u\mid
v_{c,d}^{2}+v_{c,q}^{2}\leq\left(m_{c}^{\mathrm{max}}\right)^{2}\right\\},$
(9)

where dc-link voltage is assumed to be $v_{dc}=1$ p.u. All parameter values
can be found in Table I.

###  III-A Control Barrier and Lyapunov-like Function

Figure 4:  The polynomial CBF $B(x)$ and CLF $V(x)$, along with their
corresponding safe set $\mathcal{X}_{s}$, which must be contained within the
allowed set of states $\mathcal{X}_{a}$, and nominal region $\mathcal{X}_{n}$,
are computed using SOS optimization. The vector field
$\dot{x}=f(x)+G(x)u_{sos}(x)$, projected to $(i_{d},i_{q})$ coordinates
setting $(i_{r,d},i_{r,q})=(0,1)$, is depicted in blue. Here, the polynomial
$u_{sos}(x)$ is utilized to ensure compatibility between the CBF and CLF.

The construction of the advanced safety filter relies on the development of a
CBF and CLF, which are reviewed in this subsection. We consider an abstract
polynomial control system, affine in the control action $u\in\mathbb{R}^{m}$,
and given as

$\displaystyle\dot{x}=f(x)+G(x)u,$ (10)

where $x\in\mathbb{R}^{n}$ is the state, $f(x)\in\left(R[x]\right)^{n}$ is a
polynomial vector, and $G(x)\in\left(R[x]\right)^{n\times m}$ is a polynomial
matrix. System 10 with a polynomial state feedback control policy
$u_{s}(x)\in\left(R[x]\right)^{m}$ results in the closed-loop system

$\displaystyle\dot{x}=f(x)+G(x)u_{s}(x).$ (11)

A subset $\mathcal{X}_{s}\subseteq\mathbb{R}^{n}$ is called _forward
invariant_ (cf. [14, Theorem 4.4]) with respect to system 11 if for every
$x(0)\in\mathcal{X}$, $x(t)\in\mathcal{X}$ for all $t\geq 0$. A system 11 is
called _safe_ (cf. [15]) w.r.t. an allowable set
$\mathcal{X}_{a}\subseteq\mathbb{R}^{n}$ and the safe set
$\mathcal{X}_{s}\subseteq\mathbb{R}^{n}$, if $\mathcal{X}_{s}$ is forward
invariant and $\mathcal{X}_{s}\subseteq\mathcal{X}_{a}$.

Safety of a control system 10 can be asserted with the existence of a
differentiable function $B\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ such that for all states $x\in\partial
B(x)$ there exists $u\in\mathcal{U}$:

$\displaystyle\nabla B(x)^{\top}\left(f(x)+G(x)u\right)\leq 0.$ (12)

Such a function $B(x)$ is called a Control Barrier Function (CBF), and its
zero-sublevel set $\mathcal{X}_{s}\mathrel{\mathop{\mathchar
58\relax}}=\\{x\mid B(x)\leq 0\\}$, to be contained in $\mathcal{X}_{a}$ 8,
defines the safe set.

The advanced safety filter $u_{s}(x)$ ensures a finite-time convergence to the
nominal region using a differentiable function $V\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ with a strictly positive dissipation
rate $d(x)>0$ such that for all $x\in\\{x\mid V(x)\leq x\leq B(x)\\}$ there
exists $u\in\mathcal{U}$ (cf. [13]):

$\displaystyle\nabla V(x)^{\top}\left(f(x)+G(x)u\right)+d(x)\leq 0.$ (13)

Such a function $V(x)$ is referred to as a Control Lyapunov-Like Function
(CLF), and its zero-sublevel set $\mathcal{X}_{n}\mathrel{\mathop{\mathchar
58\relax}}=\\{x\mid V(x)\leq 0\\}\subsetneq\mathcal{X}_{s}$ defines the
nominal region. To ensure compatibility with the nominal controller
$u_{n}(x)$, we furthermore require that 13 holds on the boundary of
$\mathcal{X}_{n}$ with $u=u_{n}(x)$, see [13] for details, that is for all
$x\in\partial X_{n}$:

$\displaystyle\nabla V(x)^{\top}\left(f(x)+G(x)u_{n}(x)\right)+d(x)\leq 0.$
(14)

With these conditions established, we can formulate the QCQP defining the
advanced safety filter.

###  III-B Quadratically Constrained Quadratic Program

When provided with a CBF $B(x)$ and a CLF $V(x)$ that satisfy conditions 12,
13 and 14, the advanced safety filter can be implemented using the QCQP

$\displaystyle\begin{array}[]{llll}u_{s}(x)\mathrel{\mathop{\mathchar
58\relax}}=&\underset{u\in\mathcal{U}}{\mbox{min}}&\|u_{n}(x)-u\|^{2}\\\
&\mbox{s.t.}&C(x)u+b(x)\leq 0,\end{array}$ (17)

encoding the input constraint set $\mathcal{U}$ in 9 as quadratic constraints
of the QCQP. The nominal controller

$\displaystyle u_{n}(x)=\begin{bmatrix}Z_{c}i_{r}+v_{\mathrm{PCC},f}\\\
\tau(v_{\mathrm{PCC}}-v_{\mathrm{PCC},f})\end{bmatrix}$ (18)

incorporates the behavior of the simplified EDPC 3 and the nominal rate of the
filtered PCC voltage 5. The state-dependent matrix $C(x)\in\mathbb{R}^{2\times
4}$ and vector $b(x)\in\mathbb{R}^{2}$ encode the CBF and CLF conditions 12
and 13 as follows:

$\displaystyle C(x)$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\begin{bmatrix}\nabla B(x)^{\top}G(x)\\\ \nabla
V(x)^{\top}G(x)\end{bmatrix}$ (19a) $\displaystyle b(x)$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}=\begin{bmatrix}\nabla
B(x)^{\top}f(x)\\\ \nabla
V(x)^{\top}f(x)+d(x)\end{bmatrix}-\begin{bmatrix}r_{0}(x)\\\
r_{1}(x)\end{bmatrix},$ (19b)

where $r_{0}(x)$ and $r_{1}(x)$ are state-dependent slack variables that are
non-positive for all states $x$, where the CBF or CLF conditions are
applicable. By choosing them smoothly, we can achieve a Lipschitz-continuous
QCQP-based controller, see [13, Theorem 2] for details.

###  III-C Sum-of-Squares Optimization

A polynomial CLF and CBF can be found numerically using Sum-of-Squares (SOS)
optimization. The CBF and CLF conditions 12, 13 and 14, along with the input
constraint set 9 and containment conditions
$\mathcal{X}_{n}\subsetneq\mathcal{X}_{s}\subsetneq\mathcal{X}_{a}$, are
represented as SOS constraints. This is achieved by substituting a polynomial
controller $u_{\mathrm{SOS}}(x)$ for the input vector $u$ in 12 and 13. The
SOS constraints encoding 12, 13 and 14 are obtained by applying Putinar’s
Positivstellensatz [13, Eq. (27)]:

$\displaystyle\begin{split}\nabla
B^{\top}\left(f+Gu_{\mathrm{SOS}}\right)+\gamma_{B}B+\gamma_{1}f_{\mathrm{op}}&\in-\Sigma[x]\end{split}$
(20a) $\displaystyle\begin{split}\nabla
V^{\top}\left(f+Gu_{\mathrm{SOS}}\right)+d+\gamma_{V}V-\gamma_{r}B&\\\
+\gamma_{2}f_{\mathrm{op}}&\in-\Sigma[x]\end{split}$ (20b)
$\displaystyle\begin{split}\nabla
V^{\top}\left(f+Gu_{n}\right)+d+\gamma_{n}V+\gamma_{3}f_{\mathrm{op}}&\in-\Sigma[x],\end{split}$
(20c)

for some $\gamma_{B},\gamma_{n}\in R[x]$, and
$\gamma_{V},\gamma_{r},\gamma_{1},\gamma_{2},\gamma_{3}\in\Sigma[x]$, where
$R[x]$, resp. $\Sigma[x]$, refers to the set of scalar polynomials, resp. SOS
polynomials, in variables $x\in\mathbb{R}^{n}$. An operational region
$\mathcal{X}_{\mathrm{op}}$, defined by the scalar function
$f_{\mathrm{op}}\in R[x]$, is established to encode the considered operational
ranges, where conditions 12, 13 and 14 are expected to hold.

To improve the control robustness with respect to time discretization effect,
uncertainty, and the use of the simplified model 2, we propose the following
measures:

  * •

Introduce a margin between $\mathcal{X}_{a}$ and $\mathcal{X}_{s}$.

  * •

Adjust the nominal controller $u_{n}(x)$ in 20c to expand the volume of the
nominal region $\mathcal{X}_{n}$. A larger nominal region improves stability
by taking advantage of the stable behavior of nominal control near the
equilibrium.

The resulting robust advanced safety filter can maintain safety guarantees
even under more complex grid models and demonstrates significantly smoother
behavior as it transitions into the nominal region and eventually steady-state
operations.

##  IV Numerical Results

###  IV-A Finding CBF and CLF via SOS optimization

The safe set $\mathcal{X}_{s}$ and nominal region $\mathcal{X}_{n}$, as shown
in Figure 4, are computed by using polynomials $B(x)$ and $V(x)$ of degree
two. Considering the linearity of the system dynamics 7 and the quadratic
nature of the constraints, it is possibly expected that low-degree polynomial
candidates are sufficient to find solutions to the SOS problem involving 20,
though it’s not perfectly obvious since 20 are still coupled quadratic
equations. The CBF is selected such that $B(0)=-1$ to ensure numerical
stability. The dissipation function in 13 is chosen to achieve the desired
convergence to the nominal region. The operational region

$\displaystyle\mathcal{X}_{\mathrm{op}}=\\{x\mid$
$\displaystyle(i_{r,d}/i_{r}^{\mathrm{lim}})^{2}+(i_{r,q}/i_{r}^{\mathrm{lim}})^{2}\leq
1$ $\displaystyle(v_{\mathrm{PCC},d}/0.1)^{2}+(v_{\mathrm{PCC},q}/0.1)^{2}\leq
1\\}$

is selected to include the considered $i_{r}$ and $v_{\mathrm{PCC}}$ range. To
enhance robustness with a current margin given by
$i^{\mathrm{lim}}<i^{\mathrm{max}}$, as described in Section III-C, and to
maintain compactness of the safe set w.r.t. $i$ and $v_{\mathrm{PCC},f}$, we
enforce the following subset condition:

$\displaystyle\mathcal{X}_{s}\subseteq\\{x\mid$
$\displaystyle(i_{d}/i^{\mathrm{lim}})^{2}+(i_{q}/i^{\mathrm{lim}})^{2}\leq 1$
$\displaystyle(v_{\mathrm{PCC},f,d}/20)^{2}+(v_{\mathrm{PCC},f,q}/20)^{2}\leq
1\\}\subsetneq\mathcal{X}_{a}.$

Furthermore, the nominal controller 18 used in the SOS optimization 20 is made
more aggressive as

$\displaystyle
u_{n}(x)=\begin{bmatrix}0.2(i_{r}-i)+Z_{c}i_{r}+v_{\mathrm{PCC},f}\\\
10\cdot\tau(v_{\mathrm{PCC}}-v_{\mathrm{PCC},f})\end{bmatrix}.$

A more aggressive nominal controller, while potentially sacrificing the GFM
property, can enhance the finite-time convergence and forward invariance
guarantees encoded in 20b and 20c, resulting in a larger nominal set
$\mathcal{X}_{n}$.

###  IV-B Load step simulation

Figure 5:  Simulations are conducted by switching a load to the PCC, necessitating the limitation of the BESS’s current.  TABLE I: Simulation parameters Parameter | Symbol | Value | Unit
---|---|---|---
Transformer impedance |  $l_{c}$, $r_{c}$ | 0.16, 0.01 | p.u.
Line impedance |  $l_{g}$, $r_{g}$ | 0.016, 0.001 | p.u.
Load impedance |  $l_{l}$, $r_{l}$ | 0.016, 0.001 | p.u.
Grid frequency | $f_{g}$ | 1.02 | p.u.
Current reference limit | $i_{r}^{\mathrm{lim}}$ | 1.18 | p.u.
Current limit | $i^{\mathrm{lim}}$ | 1.24 | p.u.
Maximum current | $i^{\mathrm{max}}$ | 1.30 | p.u.
Dc-link voltage | $v_{dc}$ | 1 | p.u.
Modulation limit | $m_{c}^{\mathrm{max}}$ | 1.2 | p.u.
Frequency droop | $D_{f}$ | -0.02 | p.u.
Voltage droop | $D_{v}$ | 0.05 | p.u.
Proportional gain | $K_{p}$ | 0.45 | p.u.
Integral time constant | $T_{i}$ | 40 | ms
Filter time constant | $\tau$ | 1 | ms

In this section, we investigate the behavior of the EDPC in combination with a
current limiting control during a load step applied at the PCC at $t=0.6$s, as
illustrated in Figure 5. The QCQP-based controller $u_{s}(x)$ in 17,
theoretically continuous in time, is implemented discretely with a sampling
time of $T_{s}=200\mu s$. The slack variables in 19 are chosen as follows:

$\displaystyle r_{0}(x)$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=-\gamma_{B}(x)B(x)\quad r_{1}(x)\mathrel{\mathop{\mathchar
58\relax}}=-\gamma_{V}(x)V(x),$

where $\gamma_{B}$ and $\gamma_{V}$ are taken from the SOS problem 20. The
performance of the safety filter implementation is evaluated against the
standard vector current control [4]. In this method, the vector current
control is activated whenever the current amplitude
$\|i\|=\sqrt{i_{d}^{2}+i_{q}^{2}}\geq i^{\mathrm{lim}}$, with a deactivation
hysteresis of $0.03$ p.u.

Figure 6:  The vector current control [4] is activated when the current
amplitude $\|i\|\geq i^{\mathrm{lim}}$, with a deactivation hysteresis of
$0.03$ p.u. This results in abrupt interventions and a chattering behavior
until the current oscillations are sufficiently reduced to remain below the
current amplitude limit.

Figure 7:  The designed safety filter not only smoothly limits the current
amplitude $\|i\|$ but also mitigates the current dip occurring immediately
after the load step at $t=0.6$s.

The simulations in Figures 6 and 7 show a load step using a vector current
control and a safety filter with the frequency of $v_{g}$ set to $f_{g}=1.02$
p.u. The intervention of the safety filter and the vector current control
$\Delta u$ is defined as:

$\displaystyle\Delta u(x)=u_{s}(x)-u_{n}(x).$

Unlike vector current control, which tends to exhibit abrupt interventions and
a chattering behavior even with the use of a hysteresis, the safety filter
demonstrates a significantly smoother operation. This smooth operation is
achieved by means of two key conditions. Firstly, the finite-time convergence
condition imposed by CLF dynamically adjusts the nominal control action to
mitigate the current dip immediately after the load step. Secondly, the safety
condition imposed by CBF dynamically adjusts the nominal control action as the
current approaches the allowed current amplitude $i^{\mathrm{max}}$ 6. The
intervention of the safety filter, as expressed by $\Delta u$, exhibits a
significant smaller magnitude and a smoother nature when compared to the
vector current control.

##  V Conclusions

In this paper, we demonstrated our advanced safety filter and extended the
application of safety filters to a non-trivial power electronics example
beyond their traditional domain in robotics. The finite-time convergence
guarantee of the advanced safety filter played a key role in finding the
balance between ensuring the system’s safety and maintaining the behavior of
the nominal controller. Specifically, we successfully implemented an advanced
safety filter capable of effectively limiting the current during a load step
while preserving its GFM behavior, embodied by the EDPC, during stationary
grid conditions.

## References

  * [1] Anubhav Jain, Jayachandra Naidu Sakamuri, and Nicolaos Antonio Cutululis.  Grid-forming control strategies for black start by offshore wind power plants.  Wind Energy Science, 5(4):1297–1313, 2020.
  * [2] Roberto Rosso, Xiongfei Wang, Marco Liserre, Xiaonan Lu, and Soenke Engelken.  Grid-forming converters: Control approaches, grid-synchronization, and future trends—a review.  IEEE Open Journal of Industry Applications, 2:93–109, 2021.
  * [3] Geuss Ferdinand.  Designing a controller for walenstadt’s battery energy storage system in islanded grid.  Master’s thesis, ETH Zurich, 2023.
  * [4] Mario Ndreko, Sven Rüberg, and Wilhelm Winter.  Grid forming control for stable power systems with up to 100 % inverter based generation: A paradigm scenario using the IEEE 118-bus system.  In 17th International Workshop on Large-Scale Integration of Wind Power into Power Systems as well as on Transmission Networks for Offshore Wind Power Plants, Sweden, pages 1–6, 10 2018.
  * [5] Mario Schweizer, Stefan Almér, Sami Pettersson, Arvid Merkert, Vivien Bergemann, and Lennart Harnefors.  Grid-forming vector current control.  IEEE Transactions on Power Electronics, 37(11):13091–13106, 2022\.
  * [6] Bo Fan, Teng Liu, Fangzhou Zhao, Heng Wu, and Xiongfei Wang.  A review of current-limiting control of grid-forming inverters under symmetrical disturbances.  IEEE Open Journal of Power Electronics, 2022.
  * [7] A.D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada.  Control barrier functions: Theory and applications.  In 2019 18th European control conference (ECC), pages 3420–3431. IEEE, 2019.
  * [8] Aaron D Ames, Xiangru Xu, Jessy W Grizzle, and Paulo Tabuada.  Control barrier function based quadratic programs for safety critical systems.  IEEE Transactions on Automatic Control, 62(8):3861–3876, 2016.
  * [9] W. Tan and A. Packard.  Searching for control Lyapunov functions using sums of squares programming.  sibi, 1(1), 2004.
  * [10] Andrew Clark.  Verification and synthesis of control barrier functions.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 6105–6112. IEEE, 2021.
  * [11] Han Wang, Kostas Margellos, and Antonis Papachristodoulou.  Safety verification and controller synthesis for systems with input constraints.  IFAC-PapersOnLine, 56(2):1698–1703, 2023.
  * [12] Hongkai Dai and Frank Permenter.  Convex synthesis and verification of control-lyapunov and barrier functions with input constraints.  In 2023 American Control Conference (ACC), pages 4116–4123. IEEE, 2023.
  * [13] Michael Schneeberger, Silvia Mastellone, and Florian Dörfler.  Advanced safety filter based on sos control barrier and lyapunov functions.  arXiv preprint arXiv:2401.06901, 2024.
  * [14] H.K. Khalil.  Nonlinear Systems.  Pearson Education. Prentice Hall, 2002.
  * [15] P. Wieland and F. Allgöwer.  Constructive safety using control barrier functions.  IFAC Proceedings Volumes, 40(12):462–467, 2007.