# Sufficient Conditions for Global Boundedness of Solutions for Two Coupled Synchronverters

José Ángel Mercado Uribe, Jesús Mendoza-Avila, Denis Efimov, Johannes Schiffer

► To cite this version:

José Ángel Mercado Uribe, Jesús Mendoza-Avila, Denis Efimov, Johannes Schiffer. Sufficient Conditions for Global Boundedness of Solutions for Two Coupled Synchronverters. IEEE CDC 2024 - IEEE Conference on Decision and Control, Dec 2024, Milan, Italy. hal-04895365

HAL Id: hal-04895365

https://inria.hal.science/hal-04895365v1

Submitted on 18 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sufficient Conditions for Global Boundedness of Solutions for Two
Coupled Synchronverters

Angel Mercado-Uribe, Jesús Mendoza-Ávila, Denis Efimov, Johannes Schiffer

**Abstract**— This paper analyzes two synchronverters connected in parallel to a common capacitive-resistive load through resistive-inductive power lines. This system is conceptualized as a microgrid with two renewable energy sources controlled using the synchronverter algorithm. It is modeled as an interconnection of three port-Hamiltonian systems, and the *dq*-coordinates model is derived by averaging the frequencies. Applying the recent Leonov function theory, sufficient conditions to guarantee the global boundedness of the whole system's trajectories are provided. This is necessary to reach the global synchronization of microgrids. Additionally, a numerical example illustrates the potential resonance behavior of the microgrid.

## I. INTRODUCTION

Modern society heavily relies on electrical energy sources. Consequently, analyzing power systems, which are the central column for electrical energy sources [21], is crucial to ensuring the proper functioning of the electrical grid. However, the analysis of power systems is inherently complex due to the dynamics of grid elements. Therefore, several assumptions and simplifications are usually required to analyze stability and design controllers [16], [25]. Traditionally, the Synchronous Generator (SG) is the primary power source of the current electricity grid [28], which is highly nonlinear. For this reason, to simplify its analysis, several assumptions are considered [7], [25], [27], [28]. Despite some of these assumptions being reasonable under certain operating conditions, most of them are not justifiable in generic operation scenarios [1], [16]. Furthermore, the stability and robustness properties can be just guaranteed for the reduced model, yielding to predominantly local results.

In recent years, the use of renewable energy sources has been increased [18], [20], [23], [25]. Despite the transition to renewable energy sources is desirable for reducing greenhouse gas emission [6], these sources introduce additional perturbations to the grid and the replacement of SGs by converters reduces the electrical grid's inertia and hence, the security margins as well [10], [13], [29]. Therefore, in order to position renewable energy sources as the primary element of power generation, a more rigorous analysis of stability

and robustness for multiple renewable generators is essential [13], [24], [26].

One effective strategy for controlling a converter is the synchronverter algorithm [30], which guarantees complete operational compatibility to the current electrical grid [30] and has additional benefits [3]. However, similar to SGs, the analysis of boundedness and stability of trajectories of synchronverters encounters challenges due to state periodic functions in its state dynamics, leading to multiple equilibria. This situation complicates the classical stability analysis based on Lyapunov functions and this approach is not globally feasible for most multistable systems.

An alternative approach is the concept of Leonov functions [8], [9], [15], [22]. This method leverages the periodicity and relaxes the conditions of Lyapunov functions. Using the Leonov function approach, the global Input-to-State Stability (ISS) of a SG connected to an infinite bus was proven in [24].

Following the idea of [24], the main objective of this paper consists of utilizing the Leonov function approach to establish sufficient conditions for the boundedness of trajectories of a microgrid, which is composed of two synchronverters connected to a resistive-capacitive load through resistive-inductive power lines. For the analysis, we assume that the synchronverters mimic a SG perfectly, but in contrast to standard analyses [1], [16], the whole nonlinear SG model is considered. Similarly to [4] and [11], we consider the system as an interconnection of port-Hamiltonian subsystems. However, contrary to [4] and [5], we develop a rigorous proof to guarantee the global boundedness of solutions of our particular microgrid instead of a hard-to-verify assumption. The main contributions of this paper can be summarized as:

* We derive the *abc*- and *dq*-coordinates port-Hamiltonian model for the microgrid depicted in Figure 1 by using the average frequency.

* We present sufficient conditions for global boundedness of the system trajectories by using the Hamiltonian of the system to construct a suitable Leonov function.

In contrast to our previous paper [17], we provide conditions for global stability instead of designing a controller to guarantee the ISS property for synchronverters with heterogeneous parameters.

The rest of the paper is organized as follows. Some preliminaries are presented in Section II. The port-Hamiltonian model of the system is introduced in Section III. The port-Hamiltonian model is transformed to *dq*-coordinates and developed in Section IV. The boundedness analysis is provided in Section V. A numerical example is shown in Section VI. Finally, conclusions are given in Section VII.

Angel Mercado-Uribe and Jesús Mendoza-Ávila are with the Control Systems and Network Technology Group, Brandenburgische Technische Universität Cottbus-Senftenberg, (e-mail: mercadou@b-tu.de; mendozaa@b-tu.de).

Denis Efimov is with Inria, Univ. Lille, CNRS UMR 9189 - CRISTAL, F-59000 Lille, France, (e-mail: Denis.Efimov@inria.fr).

Johannes Schiffer is with the Control Systems and Network Technology Group, Brandenburgische Technische Universität Cottbus-Senftenberg and the Fraunhofer IEG, Fraunhofer Research Institution on Energy Infrastructures and Geothermal Systems IEG, (e-mail: schiffer.b-tu.de).

This work is partially supported by the ANR-DFG project SyNPID - Project number 446182476.
## II. PRELIMINARIES

### A. Notation

Let $\mathbb{R}_+ = \{s \in \mathbb{R} : s \ge 0\}$ as well as $\mathbb{Z}_+ = \mathbb{Z} \cap \mathbb{R}_+$. The variable $v \in \mathbb{R}^k$ ($v \in \mathbb{Z}^k$) denotes a vector of real (integer) numbers with $k \in \mathbb{Z}_+$ components, whose infinity norm is represented by $\|v\|_{\infty} = \max(|v_1|, ..., |v_k|)$ and its Euclidean norm is denoted by $\|v\|_2 = \sqrt{v_1^2 + \dots + v_k^2}$. The distance from a point $p \in \mathbb{R}^n$ to the set $S \subset \mathbb{R}^n$ is defined as $\|p\|_S = \inf_{a \in S} \|p - a\|_2$. The identity matrix is represented by $\mathbf{I}_k \in \mathbb{R}^{k \times k}$. Let $i \in \mathbb{R}$, then the notations $\mathbf{i}_n$ and $\mathbf{i}_{n \times m}$, with $n, m \in \mathbb{N}$, represent a column vector of $n$ elements and a matrix with $n$ rows and $m$ columns, respectively, with all its elements being identical to $i$.

### B. State periodic systems and Leonov functions

Let the map $f: \mathbb{R}^n \to \mathbb{R}^n$ be of class $C^1$, $f(0) = 0$, and consider a nonlinear system of the following form:

$$ \dot{x}(t) = f(x(t)) \quad \forall t \ge 0, \quad (1) $$

with state $x(t) \in \mathbb{R}^n$. We denote by $X(t, x)$ the solution of (1) at time $t$ fulfilling $X(0, x_0) = x_0$ for any $x_0 \in \mathbb{R}^n$.

We make the following assumption on the dynamics (1). Let $x = (z, \theta) \in \mathbb{R}^n$, where $z \in \mathbb{R}^{n-1}$ and $\theta \in \mathbb{R}$ are two components of the state vector.

**Assumption 1:** The vector field $f$ in (1) is $2\pi$-periodic with respect to $\theta$.

That is, we assume that the dynamics (1) are periodic with respect to one state denoted by $\theta$. Thus, the system (1) can be embedded into a manifold $M = \mathbb{R}^{n-1} \times \mathbb{S}$ by a simple projection of the variable $\theta$ on the set $\mathbb{S} = [0, 2\pi)$.

Following [2], [8], [22], we introduce the next definition.

**Definition 1:** A $C^1$ function $V: \mathbb{R}^n \to \mathbb{R}$ is a Leonov function for the system (1) if there exist functions $\alpha_1, \alpha_2 \in K_\infty$, $\psi_1, \psi_2 \in K$ and a continuous function $\lambda: \mathbb{R} \to \mathbb{R}$ satisfying $\lambda(0) = 0$ and $\lambda(s)s > 0$ for all $s \ne 0$, such that

$$ \alpha_1(|z|) - \psi_1(|\theta|) \le V(x) \le \alpha_2(|z|) - \psi_2(|\theta|), \quad (2) $$

for all $x = (z, \theta) \in \mathbb{R}^n$ and the following dissipation inequality holds:

$$ \frac{\partial V(x)}{\partial x} f(x) \le -\lambda(V(x)), \quad \forall x \in \mathbb{R}^n. \quad (3) $$

Likewise, we recall the following proposition.

**Proposition 1 ([2]):** If for the system (1) there exists a Leonov function as in Definition 1, then for all $x_0 \in \mathbb{R}^n$, the trajectories $X(t; x_0)$ are bounded for all $t \ge 0$.

**Remark 1:** The definition of a Leonov function presented in [8] applies to systems with several periodic states, e.g., microgrids composed of more than two synchronizers. However, using this definition is more restrictive for systems with a scalar periodic state. So, Definition 1 is an adaptation that can only be applied to the scalar case. This definition considers sign-indefinite functions similar to the idea of Leonov presented in [15], which allows the function to be negative in the periodic variable in contrast to the standard Lyapunov function.

**Remark 2:** Despite [2, Proposition 3] is stated for a particular case, the proof covers the case presented here.

## III. SYSTEM MODEL AND PROBLEM STATEMENT

We consider a microgrid composed of two inverters operated with the synchronverter algorithm [30], which are feeding a common capacitive-resistive load, see Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Diagram of two coupled synchronizers</figcaption></figure>

We assume that all electrical quantities are balanced [25]. The parameters $R_k > 0$ and $L_k > 0$, $k = 1, 2$, are the power line resistances and inductances. The load resistance and capacitance are denoted by $R > 0$ and $C > 0$. The parameters $D_k > 0$, $J_k > 0$ and $T_{m_k} > 0$ are the virtual mechanical variables (damping, inertia and torque input) of the k-th synchronverter. The three-phase load voltage is denoted by $v_{abc}(t) \in \mathbb{R}^3$. The phase angle, electrical frequency and generated current of the k-th synchronverter are denoted by $\delta_k(t) \in \mathbb{R}$, $\omega_k(t) \in \mathbb{R}$ and $i_{abc,k}(t) \in \mathbb{R}^3$, respectively. Then, the electromotive force is given by [30]

$$ e_{abc,k} = \omega_k \varepsilon_{abc,k}, \quad \varepsilon_{abc,k} = M_{f_k} i_{f_k} \begin{bmatrix} \sin(\delta_k) \\ \sin(\delta_k - \frac{2}{3}\pi) \\ \sin(\delta_k + \frac{2}{3}\pi) \end{bmatrix}, $$

where $M_{f_k}$ is the mutual inductance and $i_{f_k}$ is the rotor current, which is assumed constant. Likewise, the electrical torque $T_{e_k}$ is computed as $T_{e_k} = \varepsilon_{abc,k}^\top i_{abc,k}$ [30]. Thus, the system in Figure 1 can be modeled in abc-coordinates as

$$ \begin{align} J_k \dot{\omega}_k &= -D_k \omega_k - \varepsilon_{abc,k}^\top i_{abc,k} + T_{m_k}, \quad \dot{\delta}_k = \omega_k, \notag \\ L_k \frac{d i_{abc,k}}{dt} &= -R_k i_{abc,k} + \varepsilon_{abc,k} \omega_k - v_{abc}, \quad k = 1, 2; \tag{4} \\ C \frac{dv_{abc}}{dt} &= -\frac{1}{R} v_{abc} + i_{abc,1} + i_{abc,2}. \notag \end{align} $$

Inspired by [4], [11] (see also the correction [5]), we split the system (4) into two parts, namely the mechanical and the electrical one, and derive port-Hamiltonian representations for each of these system parts. Defining the matrices

$$ \mathcal{R}_m = \begin{bmatrix} \frac{D_1}{J_1^2} & 0 \\ 0 & \frac{D_2}{J_2^2} \end{bmatrix}, \quad g_{\xi_m} = \begin{bmatrix} \frac{\varepsilon_{abc,1}}{J_1} & \mathbf{0}_{1 \times 3} \\ \mathbf{0}_{1 \times 3} & \frac{\varepsilon_{abc,2}}{J_2} \end{bmatrix}, \quad C_m = \begin{bmatrix} \frac{T_{m_1}}{J_1} \\ \frac{T_{m_2}}{J_2} \end{bmatrix}, $$

and the variables

$$ \xi_m = \begin{bmatrix} \omega_1 \\ \omega_2 \end{bmatrix}, \quad \mu_m = - \begin{bmatrix} i_{abc,1} \\ i_{abc,2} \end{bmatrix}, \quad \sigma_m = \begin{bmatrix} \varepsilon_{abc,1} \omega_1 \\ \varepsilon_{abc,2} \omega_2 \end{bmatrix}, $$

the mechanical part of the system (4) can be written as

$$ \begin{align} \dot{\xi}_m &= - \mathcal{R}_m \nabla H_{\xi_m}(\xi_m) + g_{\xi_m} \mu_m + C_m, \tag{5} \\ \sigma_m &= g_{\xi_m}^\top \nabla H_{\xi_m}(\xi_m), \notag \end{align} $$

where the Hamiltonian function and its gradient are

$$ H_{\xi_m}(\xi_m) = \frac{J_1}{2} \omega_1^2 + \frac{J_2}{2} \omega_2^2, \quad \nabla H_{\xi_m} = \begin{bmatrix} J_1 \omega_1 \\ J_2 \omega_2 \end{bmatrix}. \quad (6) $$
Analogously, defining the variables

$$
\xi_e = \begin{bmatrix} i_{abc,1} \\ i_{abc,2} \\ v_{abc} \end{bmatrix}, \quad \mu_e = \begin{bmatrix} \varepsilon_{abc,1} \omega_1 \\ \varepsilon_{abc,2} \omega_2 \end{bmatrix}, \quad \sigma_e = \begin{bmatrix} i_{abc,1} \\ i_{abc,2} \end{bmatrix},
$$

and the matrices

$$
R_{\xi_e} = \operatorname{diag} \left( \frac{R_1 \mathbf{I}_3}{L_1^2}, \frac{R_2 \mathbf{I}_3}{L_2^2}, \frac{\mathbf{I}_3}{RC^2} \right), g_{\xi_e} = \begin{bmatrix} \operatorname{diag} \left( \frac{\mathbf{I}_3}{L_1}, \frac{\mathbf{I}_3}{L_2} \right) \\ \mathbf{0}_{3 \times 6} \end{bmatrix},
$$

$$
\mathcal{J}_{\xi_e} = \begin{bmatrix} \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & -\frac{1}{CL_1}\mathbf{I}_3 \\ \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & -\frac{1}{CL_2}\mathbf{I}_3 \\ \frac{1}{CL_1}\mathbf{I}_3 & \frac{1}{CL_2}\mathbf{I}_3 & \mathbf{0}_{3 \times 3} \end{bmatrix},
$$

the electrical part of the system (4) can be written as

$$
\begin{align}
\dot{\xi}_e &= [\mathcal{J}_{\xi_e} - \mathcal{R}_{\xi_e}] \nabla H_{\xi_e}(\xi_e) + g_{\xi_e} \mu_e, \tag{7} \\
\sigma_e &= g_{\xi_e}^{\top} \nabla H_{\xi_e}(\xi_e), \notag
\end{align}
$$

where the Hamiltonian function and its gradient are

$$
H_{\xi_e}(\xi_e) = \frac{L_1}{2} \|i_{abc,1}\|^2 + \frac{L_2}{2} \|i_{abc,2}\|^2 + \frac{C}{2} \|v_{abc}\|^2, \quad (8)
$$

$$
\nabla H_{\xi_e}(\xi_e) = [L_1 i_{abc,1}^{\top} \quad L_2 i_{abc,2}^{\top} \quad C v_{abc}^{\top}]^{\top}. \quad (9)
$$

Therefore, defining the variable $\xi = [\xi_m^\top \quad \xi_e^\top]$ and combining the systems (5) and (7), we obtain the whole system (4) in compact form as follows

$$
\begin{align}
\begin{bmatrix} \dot{\delta}_1 \\ \dot{\delta}_2 \end{bmatrix} &= \begin{bmatrix} \frac{1}{J_1} & 0 & \mathbf{0}_{1 \times 9} \\ 0 & \frac{1}{J_2} & \mathbf{0}_{1 \times 9} \end{bmatrix} \nabla H_{\xi}(\xi), \tag{10} \\
\dot{\xi} &= \begin{bmatrix} -\mathcal{R}_m & -g_{\xi_m} g_{\xi_e}^\top \\ g_{\xi_e} g_{\xi_m}^\top & \mathcal{J}_{\xi_e} - \mathcal{R}_{\xi_e} \end{bmatrix} \nabla H_{\xi}(\xi) + \begin{bmatrix} \mathcal{C}_m \\ \mathbf{0}_9 \end{bmatrix}, \notag
\end{align}
$$

where the Hamiltonian of the whole system is

$$
H_{\xi}(\xi) = H_{\xi_m}(\xi_m) + H_{\xi_e}(\xi_e). \quad (11)
$$

Using this system, the main problem considered in the
present paper is formalized below.

*Problem 1:* Consider the system (10). Derive conditions under which its solutions are globally bounded.

IV. PORT-HAMILTONIAN MODEL OF TWO COUPLED
SYNCHRONVERTERS IN dq-COORDINATES

In order to address Problem 1, we introduce the dq-
transformation [25]

$$
T_{dq}(\psi) = \sqrt{\frac{2}{3}} \begin{bmatrix} \cos(\psi) & \cos(\psi - \frac{2}{3}\pi) & \cos(\psi + \frac{2}{3}\pi) \\ \sin(\psi) & \sin(\psi - \frac{2}{3}\pi) & \sin(\psi + \frac{2}{3}\pi) \end{bmatrix} \quad (12)
$$

and the rotation matrix

$$
T_r(\nu) = \begin{bmatrix} \cos(\nu) & \sin(\nu) \\ -\sin(\nu) & \cos(\nu) \end{bmatrix}, \quad T_r(-\nu) = T_r^{\top}(\nu). \quad (13)
$$

Recall that for $\varphi = \psi - \nu$ [25],

$$
T_{dq}(\varphi) = T_r(\nu)T_{dq}(\psi). \quad (14)
$$

Therefore, defining $\phi := \frac{1}{2}(\delta_1 + \delta_2)$, whose derivative is
$\dot{\phi} = \frac{1}{2}(\omega_1 + \omega_2)$, and $\theta := \delta_1 - \phi = -(\delta_2 - \phi) = \frac{1}{2}(\delta_1 - \delta_2)$,
we propose the following change of coordinates

$$
x_m = \xi_m, \quad x_e = \mathcal{T}_{dq}(\phi)\xi_e, \quad \zeta = [x_m \ x_e]^\top, \quad (15)
$$

where

$$
\mathcal{T}_{dq}(\phi) = \begin{bmatrix} T_{dq}(\phi) & \mathbf{0}_{2 \times 3} & \mathbf{0}_{2 \times 3} \\ \mathbf{0}_{2 \times 3} & T_{dq}(\phi) & \mathbf{0}_{2 \times 3} \\ \mathbf{0}_{2 \times 3} & \mathbf{0}_{2 \times 3} & T_{dq}(\phi) \end{bmatrix}.
$$

Using the change of coordinates (15), the interconnection
between the mechanical and electrical sides is given by

$$
\begin{align*}
g_{x_m} \mu_{x_m} &= - \mathcal{T}_r(\theta) \operatorname{diag}(L_1 \mathbf{I}_2, L_2 \mathbf{I}_2, C \mathbf{I}_2) x_e, \\
g_{x_e} \mu_{x_e} &= \mathcal{T}_r(\theta)^{\top} \operatorname{diag}(J_1, J_2) x_m,
\end{align*}
$$

where

$$
\mathcal{T}_r(\theta) = \begin{bmatrix} \mathbf{b}_1^\top T_r(\theta)^\top & \mathbf{0}_{1 \times 2} & \mathbf{0}_{1 \times 2} \\ \mathbf{0}_{1 \times 2} & \mathbf{b}_2^\top T_r(\theta) & \mathbf{0}_{1 \times 2} \end{bmatrix}
$$

with $\mathbf{b}_k = \frac{b_k}{J_k L_k} [\mathbf{0} \quad \mathbf{1}]^\top$, $b_k = \sqrt{\frac{3}{2}} M_{fk} i_{fk}$ for $k=1,2$, plays an essential role. This term translates signals between the *dq*-frame referenced to the average frequency $\phi$ and the *dq*-frame referenced to the respective synchronverter angle. So, using the angle dynamics $\dot{\theta} = \frac{1}{2}(\omega_1 - \omega_2)$, the change of coordinates (15) and defining the matrices

$$
\begin{align*}
\mathcal{R} &= \begin{bmatrix} \mathcal{R}_m & \mathbf{0}_{2 \times 6} \\ \mathbf{0}_{6 \times 2} & \mathcal{R}_e \end{bmatrix}, \quad \mathcal{J}(\theta, \zeta) = \begin{bmatrix} \mathbf{0}_{2 \times 2} & -\mathcal{T}_r(\theta) \\ \mathcal{T}_r(\theta)^\top & \mathcal{J}_e(\zeta) \end{bmatrix}, \\
\mathcal{J}_e(\zeta) &= \begin{bmatrix} \frac{\bar{\omega}}{L_1} \mathcal{J}_2 & \mathbf{0}_{2 \times 2} & -\frac{1}{CL_1} \mathbf{I}_2 \\ \mathbf{0}_{2 \times 2} & \frac{\bar{\omega}}{L_2} \mathcal{J}_2 & -\frac{1}{CL_2} \mathbf{I}_2 \\ \frac{1}{CL_1} \mathbf{I}_2 & \frac{1}{CL_2} \mathbf{I}_2 & \frac{\bar{\omega}}{C} \mathcal{J}_2 \end{bmatrix}, \quad \bar{\omega} = \frac{\zeta_1 + \zeta_2}{2}, \\
\mathcal{J}_2 &= \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}, \quad \mathcal{R}_e = \operatorname{diag} \left( \frac{R_1}{L_1^2} \mathbf{I}_2, \frac{R_2}{L_2^2} \mathbf{I}_2, \frac{\mathbf{I}_2}{RC^2} \right),
\end{align*}
$$

the system (10) can be expressed as

$$
\begin{align}
\dot{\theta} &= \begin{bmatrix} \frac{1}{2J_1} & -\frac{1}{2J_2} & \mathbf{0}_{1 \times 6} \end{bmatrix} \nabla H(\zeta), \tag{16} \\
\dot{\zeta} &= \begin{bmatrix} \mathcal{J}(\theta, \zeta) - \mathcal{R} \end{bmatrix} \nabla H(\zeta) + \begin{bmatrix} \mathcal{C}_m^\top & \mathbf{0}_{1 \times 6} \end{bmatrix}^\top, \notag
\end{align}
$$

whose Hamiltonian is given by

$$
H(\zeta) = \frac{1}{2} \zeta^{\top} \mathrm{diag}(J_1, J_2, L_1 \mathbf{I}_2, L_2 \mathbf{I}_2, C \mathbf{I}_2) \zeta. \quad (17)
$$

V. CONDITIONS FOR GLOBAL BOUNDEDNESS OF
SOLUTIONS

We consider the following assumption, see also [12].

*Assumption 2:* The system (16) possesses at least one equilibrium point $(\theta^*, \zeta^*)$.

Assumption (2) is a natural requirement. Furthermore, we
define the synchronization frequency as $\omega^* = \zeta_1^* = \zeta_2^*$.

A. Shifted Hamiltonian

Defining the variable $\tilde{\zeta} = \zeta - \zeta^*$, the shifted Hamiltonian for the system (16) can be written as follows

$$
\tilde{H}(\tilde{\zeta}) = H(\tilde{\zeta} + \zeta^*) - \tilde{\zeta}^{\top} \nabla H(\zeta^*) - H(\zeta^*), \quad (18)
$$

whose gradient is

$$
\nabla \tilde{H}(\zeta) = \nabla H(\zeta) - \nabla H(\zeta^*).
$$

Thus, under Assumption 2 and using the shifted Hamiltonian
in (18), the system (16) in coordinates $\tilde{\zeta}$ results in

$$
\begin{align}
\dot{\theta} &= \begin{bmatrix} \frac{1}{2J_1} & -\frac{1}{2J_2} & \mathbf{0}_{1 \times 6} \end{bmatrix} \nabla \tilde{H}(\tilde{\zeta}), \nonumber \\
\dot{\tilde{\zeta}} &= \begin{bmatrix} \mathcal{J}(\theta, \zeta) - \mathcal{R} \end{bmatrix} \nabla \tilde{H}(\tilde{\zeta}) + (\mathcal{J}(\theta, \zeta) - \mathcal{J}(\theta^*, \zeta^*)) \nabla H(\zeta^*). \tag{19}
\end{align}
$$
1) Derivative of the Shifted Hamiltonian Function: Since $\mathcal{J}(\theta, \zeta)$ is a skew-symmetric matrix, the derivative of $\tilde{H}$ along the trajectories of the system (19) results in

$$ \dot{\tilde{H}} = -\nabla \tilde{H}^{\top} \mathcal{R} \nabla \tilde{H} + \nabla \tilde{H}^{\top} (\mathcal{J}(\theta, \zeta) - \mathcal{J}(\theta^{*}, \zeta^{*})) \nabla H(\zeta^{*}). \quad (20) $$

By following [4], [5], equation (20) can be rewritten as

$$ \dot{\tilde{H}} = -\nabla \tilde{H}^{\top} (\mathcal{R} - M(\zeta^{*})) \nabla \tilde{H} + \nabla \tilde{H}^{\top} T(\theta, \theta^{*}) \nabla H(\zeta^{*}), \quad (21) $$

where

$$ \begin{aligned} M(\zeta^*) &= \begin{bmatrix} \mathbf{0}_{2 \times 2} & M_{12}(\zeta^*) \\ M_{12}^{\top}(\zeta^*) & \mathbf{0}_{6 \times 6} \end{bmatrix}, \\ M_{12}(\zeta^*) &= \begin{bmatrix} -\frac{\zeta_4^*}{4J_1} & \frac{\zeta_3^*}{4J_1} & -\frac{\zeta_6^*}{4J_1} & \frac{\zeta_5^*}{4J_1} & -\frac{\zeta_8^*}{4J_1} & \frac{\zeta_7^*}{4J_1} \\ -\frac{\zeta_4^*}{4J_2} & \frac{\zeta_3^*}{4J_2} & -\frac{\zeta_6^*}{4J_2} & \frac{\zeta_5^*}{4J_2} & -\frac{\zeta_8^*}{4J_2} & \frac{\zeta_7^*}{4J_2} \end{bmatrix}, \end{aligned} $$

is a symmetric matrix and

$$ T(\theta, \theta^*) = \begin{bmatrix} \mathbf{0}_{2 \times 2} & -\Delta \tau_r(\theta, \theta^*) \\ \Delta \tau_r^\top(\theta, \theta^*) & \mathbf{0}_{6 \times 6} \end{bmatrix}, \quad (22) $$

with $\Delta \tau_r(\theta, \theta^*) = \tau_r(\theta) - \tau_r(\theta^*)$, is a skew-symmetric one. Therefore, $\tilde{H}$ in (21) can be written as

$$ \dot{\tilde{H}} = -\nabla \tilde{H}^{\top} \mathcal{R}_M(\zeta^*) \nabla \tilde{H} + \nabla \tilde{H}^{\top} \mathcal{T}(\zeta^*) \begin{bmatrix} \Delta \cos \\ \Delta \sin \end{bmatrix}, \quad (23) $$

where $\mathcal{R}_M(\zeta^*) = \mathcal{R} - M(\zeta^*)$ is a symmetric matrix, and

$$ \mathcal{T}^{\top}(\zeta^*) = \begin{bmatrix} \frac{-b_1 \zeta_4^*}{J_1} & \frac{-b_2 \zeta_6^*}{J_2} & 0 & \frac{b_1 \omega^*}{L_1} & 0 & \frac{b_2 \omega^*}{L_2} & \mathbf{0}_{1 \times 2} \\ \frac{-b_1 \zeta_3^*}{J_1} & \frac{b_2 \zeta_5^*}{J_2} & \frac{b_1 \omega^*}{L_1} & 0 & \frac{-b_2 \omega^*}{L_2} & 0 & \mathbf{0}_{1 \times 2} \end{bmatrix}, $$

with $\Delta \cos = \cos(\theta) - \cos(\theta^*)$, $\Delta \sin = \sin(\theta) - \sin(\theta^*)$.

## B. Boundedness of $\tilde{\zeta}$-dynamics

In order to formulate our main result, we introduce two key lemmas, whose proofs are omitted due to space limitations.

**Lemma 1:** Consider the function

$$ \begin{aligned} f(\theta) = & \zeta_4^* (\sin(\theta) - \sin(\theta^*) - (\theta - \theta^*) \cos(\theta^*)) - \\ & \zeta_3^* (\cos(\theta) - \cos(\theta^*) + (\theta - \theta^*) \sin(\theta^*)) . \end{aligned} \quad (24) $$

For all $\theta \in \mathbb{R}$, it holds that

$$ |f(\theta)| \le \frac{1}{2} \varepsilon_1 (\theta - \theta^*)^2, \quad \varepsilon_1 = \sqrt{(\zeta_3^*)^2 + (\zeta_4^*)^2}. \quad (25) $$

**Lemma 2:** The matrix

$$ M_{\lambda_2} = \begin{bmatrix} M_{\lambda_{211}} & M_{\lambda_{212}} \\ M_{\lambda_{212}}^{\top} & \frac{1}{2} \lambda (\gamma_1 - b_1 \varepsilon_1) \end{bmatrix} \quad (26) $$

with

$$ M_{\lambda_{211}} = \operatorname{diag} \left( \frac{R_1 - \lambda L_1}{L_1^2} \mathbf{I}_2, \frac{R_2 - \lambda L_2}{L_2^2} \mathbf{I}_2, \frac{1 - \lambda RC}{RC^2} \mathbf{I}_2 \right), $$

$$ M_{\lambda_{212}} = \begin{bmatrix} \frac{b_1 \omega^*}{2L_1} & \mathbf{1}_{1 \times 2} & \frac{b_2 \omega^*}{2L_2} & \mathbf{1}_{1 \times 2} & \mathbf{0}_{1 \times 2} \end{bmatrix}^{\top}, $$

where $\varepsilon_1$ is defined in (25), is positive definite if

$$ \lambda < \min \left( \frac{R_1}{L_1}, \frac{R_2}{L_2}, \frac{1}{RC} \right), \quad (27) $$

$$ \gamma_1 > b_1 \varepsilon_1 + \frac{1}{\lambda} \left( \frac{b_1^2}{R_1 - \lambda L_1} + \frac{b_2^2}{R_2 - \lambda L_2} \right) (\omega^*)^2. \quad (28) $$

Using these lemmas, the following theorem can be stated.

**Theorem 1:** Assume that $\lambda$ and $\gamma_1$ take values such that the inequalities (27) and (28) are satisfied. Then, the trajectories of the system (19) are bounded if there exist some values $\varepsilon_1$ and $\varepsilon_2 \in \mathbb{R}_+$, such that

$$ \begin{aligned} D_1 &\ge \lambda J_1 + \frac{16 J_1^2 \varepsilon_1 + \sigma_{\zeta^*} + 4 \gamma_1^2}{16 \lambda_{\min}(M_{\lambda_2})}, \\ D_2 &\ge \lambda J_2 + \frac{16 J_2^2 \varepsilon_2 + \sigma_{\zeta^*} + 4 (\gamma_1 + \sigma_f)^2}{16 \lambda_{\min}(M_{\lambda_2})}, \end{aligned} \quad (29) $$

$$ \varepsilon_1 \varepsilon_2 \ge \frac{1}{16^2 J_1^2 J_2^2} (\sigma_{\zeta^*} + 4 (\gamma_1^2 + \gamma_1 \sigma_f))^2, $$

where $\sigma_{\zeta^*} = (\zeta_3^*)^2 + (\zeta_4^*)^2 + (\zeta_5^*)^2 + (\zeta_6^*)^2 + (\zeta_7^*)^2 + (\zeta_8^*)^2$, $\sigma_f = |b_2 \zeta_6^* + b_1 \zeta_4^*| + |b_2 \zeta_5^* - b_1 \zeta_3^*|$ and $\lambda_{\min}(M_{\lambda_2})$ represents the smallest eigenvalue of the matrix $M_{\lambda_2}$ in (26).

**Proof:** Based on the Hamiltonian structure, we consider the candidate Leonov function

$$ V = \tilde{H} - \frac{1}{2} \gamma_1 (\theta - \theta^*)^2 + b_1 f (\theta - \theta^*) \quad (30) $$

where $f(\theta - \theta^*)$ is defined in (24). The derivative of $V$ along the solutions of (19) can be written as

$$ \dot{V} = \dot{\tilde{H}} - \gamma_1 (\theta - \theta^*) \eta^{\top} \nabla \tilde{H} + b_1 \nabla f \eta^{\top} \nabla \tilde{H}, \quad (31) $$

where $\eta = \begin{bmatrix} \frac{1}{J_1} & -\frac{1}{J_2} & \mathbf{0}_{1 \times 6} \end{bmatrix}^{\top}$. Substituting equation (23) in (31), we obtain

$$ \dot{V} = - \begin{bmatrix} \nabla \tilde{H} \\ \theta - \theta^* \\ \Delta \cos \\ \Delta \sin \end{bmatrix}^{\top} \begin{bmatrix} \mathcal{R}_M(\zeta^*) & \frac{1}{2} \gamma_1 \eta & -\frac{\tau_f(\zeta^*)}{2} \\ \frac{1}{2} \gamma_1 \eta^{\top} & 0 & \mathbf{0}_{1 \times 2} \\ -\frac{\tau_f^{\top}(\zeta^*)}{2} & \mathbf{0}_{2 \times 1} & \mathbf{0}_{2 \times 2} \end{bmatrix} \begin{bmatrix} \nabla \tilde{H} \\ \theta - \theta^* \\ \Delta \cos \\ \Delta \sin \end{bmatrix}, $$

where

$$ \mathcal{T}_f(\zeta^*) = \mathcal{T}(\zeta^*) - \begin{bmatrix} \frac{b_1}{J_1} \zeta_4^* \mathbf{1}_2 & -\frac{b_1}{J_2} \zeta_4^* \mathbf{1}_2 & \mathbf{0}_{2 \times 6} \end{bmatrix}^{\top} $$

Adding and subtracting the term $\lambda V$, with $\lambda \in \mathbb{R}_+$, we obtain

$$ \begin{aligned} \dot{V} = - & \begin{bmatrix} \nabla \tilde{H} \\ \theta - \theta^* \\ \Delta \cos \\ \Delta \sin \end{bmatrix}^{\top} \begin{bmatrix} \mathcal{R}_M(\zeta^*) & \frac{1}{2} \gamma_1 \eta & -\frac{\tau_f(\zeta^*)}{2} \\ \frac{1}{2} \gamma_1 \eta^{\top} & 0 & \mathbf{0}_{1 \times 2} \\ -\frac{\tau_f^{\top}(\zeta^*)}{2} & \mathbf{0}_{2 \times 1} & \mathbf{0}_{2 \times 2} \end{bmatrix} \begin{bmatrix} \nabla \tilde{H} \\ \theta - \theta^* \\ \Delta \cos \\ \Delta \sin \end{bmatrix} + \\ & \lambda \left( \tilde{H} - \frac{1}{2} \gamma_1 (\theta - \theta^*)^2 + b_1 f (\theta - \theta^*) \right) - \lambda V. \end{aligned} \quad (32) $$

We recall that $\tilde{H} = \nabla \tilde{H}^{\top} P_H \nabla \tilde{H}$, with $P_H = \operatorname{diag}(\frac{1}{J_1}, \frac{1}{J_2}, \frac{1}{L_1}, \frac{1}{L_2}, \frac{1}{L_2}, \frac{1}{C}, \frac{1}{C})$. Likewise, from the Taylor series, we have the following inequalities

$$ |\Delta \cos| \le |\theta - \theta^*| \quad \text{and} \quad |\Delta \sin| \le |\theta - \theta^*|. $$

Therefore, using these inequalities and Lemma 1, the right-hand side in (32) can be bounded by

$$ \dot{V} \le -\lambda V - \begin{bmatrix} \nabla \tilde{H} \\ \theta - \theta^* \end{bmatrix}^{\top} M_{\lambda} \begin{bmatrix} \nabla \tilde{H} \\ \theta - \theta^* \end{bmatrix}, \quad (33) $$

with

$$ M_{\lambda} = \begin{bmatrix} M_{\lambda_1} & M_{\lambda_{12}} \\ M_{\lambda_{12}}^{\top} & M_{\lambda_2} \end{bmatrix}, $$
where $M_{\lambda_2}$ is defined in (26),

$$M_{\lambda_1} = \begin{bmatrix} \frac{D_1 - \lambda J_1}{J_1^2} & 0 \\ 0 & \frac{D_2 - \lambda J_2}{J_2^2} \end{bmatrix}, \quad M_{\lambda_2} = \begin{bmatrix} M_{12} & \begin{bmatrix} \frac{\gamma_1}{\gamma_1 + \sigma_f} \\ \frac{2J_1}{2J_2} \end{bmatrix} \end{bmatrix},$$

$$\text{with } \sigma_f = |b_2\zeta_6^* + b_1\zeta_4^*| + |b_2\zeta_5^* - b_1\zeta_3^*|.$$

Recalling Lemma 2, the matrix $M_{\lambda_2}$ is positive if the conditions (27) and (28) are satisfied. Therefore, the whole matrix $M_{\lambda}$ is positive semidefinite if

$$M_{\lambda_1} - M_{\lambda_2} M_{\lambda_2}^{-1} M_{\lambda_1 \lambda_2}^{\top} \geq 0. \quad (34)$$

Since $M_{\lambda_2}$ is a positive definite matrix, we have $M_{\lambda_2} \ge \lambda_{\min}(M_{\lambda_2})\mathbf{I}_7$, which implies $M_{\lambda_2}^{-1} \le (\lambda_{\min}(M_{\lambda_2})\mathbf{I}_7)^{-1}$, where $\lambda_{\min}(\cdot)$ represents the smallest eigenvalue of a matrix. Therefore, the inequality (34) can be estimated by

$$\lambda_{\min}(M_{\lambda_2}) M_{\lambda_1} - M_{\lambda_1 2} M_{\lambda_2 1}^{\top} \geq 0. \quad (35)$$

Defining the variables

$$\sigma_{\zeta^*} = (\zeta_3^*)^2 + (\zeta_4^*)^2 + (\zeta_5^*)^2 + (\zeta_6^*)^2 + (\zeta_7^*)^2 + (\zeta_8^*)^2, \\ \sigma_{D_k} = \lambda_{\min}(M_{\lambda_2})(D_k - \lambda J_k), \quad k=1,2,$$

the inequality (35) can be written as

$$\left[ \begin{array}{cc} \frac{16\sigma_{D_1} - \sigma_{\zeta^*} - 4\gamma_1^2}{16J_1^2} & -\frac{\sigma_{\zeta^*} + 4(\gamma_1^2 + \gamma_1\sigma_f)}{16J_1J_2} \\ -\frac{\sigma_{\zeta^*} + 4(\gamma_1^2 + \gamma_1\sigma_f)}{16J_1J_2} & \frac{16\sigma_{D_2} - \sigma_{\zeta^*} - 4(\gamma_1 + \sigma_f)^2}{16J_2^2} \end{array} \right] \geq 0.$$

Selecting $D_1$ and $D_2$ such that

$$\frac{16\lambda_{\min}(M_{\lambda_2})(D_1 - \lambda J_1) - \sigma_{\zeta^*} - 4\gamma_1^2}{16J_1^2} = \epsilon_1 > 0, \\ \frac{16\lambda_{\min}(M_{\lambda_2})(D_2 - \lambda J_2) - \sigma_{\zeta^*} - 4(\gamma_1 + \sigma_f)^2}{16J_2^2} = \epsilon_2 > 0, \quad (36)$$

the inequality (35) is satisfied if

$$\epsilon_1 \epsilon_2 \ge \frac{1}{16^2 J_1^2 J_2^2} (\sigma_{\zeta^*} + 4 (\gamma_1^2 + \gamma_1 \sigma_f))^2.$$

Thus, we have $\dot{V} \le -\lambda V$ if the inequalities (27)-(29) are satisfied. From Lemma 1 and inequality (28), $V$ satisfies the property (2) in Definition 1. Therefore, the function (30) is a Leonov function for the system (19) and the trajectories of the system (19) are globally bounded by Proposition 1. ■

Theorem 1 suggests that the trajectories of the microgrid can be made bounded by using sufficiently large damping coefficients, which is expected for this kind of systems.

## VI. NUMERICAL EXAMPLE

For simulation, we consider the system in dq-coordinates (16). Likewise, based on [17], [19], we consider the following parameters $J_1 = J_2 = 0.1[kg \cdot m^2/rad]$, $b_1 = b_2 = \sqrt{3}/2[V \cdot s]$, $R_1 = R_2 = 0.75[\Omega]$, $L_1 = L_2 = 2.2[mH]$, $R = 10[\Omega]$ and $C = 1[mF]$. Meanwhile, the parameters $D_1$, $D_2$, $T_{m1}$ and $T_{m2}$ will be used for comparison result by fixing the frequency synchronization in $\omega^* = 50[\text{rad/s}]$.

From (16), we can calculate the equilibria sets. As we are using identical electrical parameters and based on [17], [28], we would expect two different equilibria sets. One of them,

is obtained for $\theta = 2k\pi$, $k \in \mathbb{Z}$, and the other one is got for $\theta = (2k-1)\pi$, $k \in \mathbb{Z}$. So, using $\theta = 0$, we fix the mechanical inputs as $T_{m1} = D_1\omega^* + b_1\zeta_4^*$ and $T_{m2} = D_2\omega^* + b_2\zeta_6^*$ to obtain the following equilibrium points

$$\begin{align} \theta^* &= 2k\pi, \quad \zeta_1^* = \zeta_2^* = 50, \quad \zeta_3^* = \zeta_5^* = 1.4094, \tag{37} \\ \zeta_4^* &= \zeta_6^* = 2.9921, \quad \zeta_7^* = -1.3862, \quad \zeta_8^* = 59.1482, \nonumber \end{align}$$

with $k \in \mathbb{Z}$, which are independent on $D_1$ and $D_2$. So, the values of $D_1$ and $D_2$ do not alter the equilibria (37).

From inequality (27), we get $\lambda < 100$. Fixing $\lambda = 90$ in (28), we get $\gamma_1 > 156.0169$. Choosing $\gamma_1 = 340$ and fixing $\epsilon_1 = \epsilon_2$, we obtain $D_1 > 16.5235$ and $D_2 > 16.6845$ from the inequalities in (29).

For simulations, we fix $D_2 = 16.7[\text{N} \cdot \text{m} \cdot \text{s/rad}]$ and $D_1$ takes values in $\{1, 2, 16.6\}[\text{N} \cdot \text{m} \cdot \text{s/rad}]$; where only the last value satisfies our conditions. Likewise, we employ the Euler solver in Simulink MATLAB with a sample time 65 [$\mu$s], which corresponds to the sample time of the converters in our laboratory [14]. We assume $\theta(0) = 0$ and $\tilde{\zeta}(0) = -\pi \mathbf{1}_8$.

Before presenting the results, we stress that the equilibria in (37) were verified by simulations for all cases.

In Figure 2, the simulation results for $D_1 = 1[\text{N} \cdot \text{m} \cdot \text{s/rad}]$ are shown. We can see that the difference between synchronizers' phases $\theta$ diverges, which implies that the frequencies $\zeta_1$ and $\zeta_2$ cannot reach synchronization, i.e., $\tilde{\zeta}_1$ and $\tilde{\zeta}_2$ do not converge to the same constant.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Trajectories of the system (19) for $D_1 = 1[\text{N} \cdot \text{m} \cdot \text{s/rad}]$.</figcaption></figure>

In Figure 3, the simulation results for $D_1 = 2[\text{N} \cdot \text{m} \cdot \text{s/rad}]$ are presented. In this case, the states converge to the equilibria set, with $\theta = (2k-1)\pi$ for $k \in \mathbb{Z}$. Nevertheless, we can see that these trajectories are indeed bounded. Likewise, we can note that $\tilde{\zeta}_1$ and $\tilde{\zeta}_2$ converge to zero.

In Figure 4, the simulation results for $D_1 = 16.6[\text{N} \cdot \text{m} \cdot \text{s/rad}]$ are illustrated. Here, all trajectories converge to the origin.

It is clear from these simulations that the minimum damping coefficient for $D_1$, with $D_2 = 16.7[\text{N} \cdot \text{m} \cdot \text{s/rad}]$, falls within the range of 1 to $2[\text{N} \cdot \text{m} \cdot \text{s/rad}]$. This range is not so far to the minimum obtained with inequalities (27)-(29).

## VII. CONCLUSIONS AND FUTURE WORK

In this paper, a model of two-parallel synchronizers in port-Hamiltonian form is presented and a reduced-order model is obtained by applying the dq-transformation.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Trajectories of the system (19) for $D_1 = 2[\frac{\text{N}\cdot\text{m}\cdot\text{s}}{\text{rad}}]$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Trajectories of the system (19) for $D_1 = 16.6[\frac{\text{N}\cdot\text{m}\cdot\text{s}}{\text{rad}}]$.</figcaption></figure>

Using the dq mathematical model (19), sufficient conditions for the damping coefficients to ensure global boundedness of trajectories in two-parallel synchronverters with a capacitive-resistive load are presented. The result is obtained by constructing a Leonov function for the whole system.

Future work aims to find a weak Lyapunov function to provide sufficient conditions for almost global stability of the invariant set of system (19). Likewise, we will extend the microgrid to a broader amount of electrical energy sources.

## REFERENCES

[1] P. M. Anderson and A. A. Fouad. *Power System Control and Stability*. NY: Wiley, New York, 2002.

[2] N. Barabanov, J. Schiffer, R. Ortega, and D. Efimov. Conditions for almost global attractivity of a synchronous generator connected to an infinite bus. *Transactions on Automatic Control*, 62(10):4905–4916, 2017.

[3] M. Blau and G. Weiss. Synchronverters used for damping inter-area oscillations in two-area power systems. *Renewable Energy and Power Quality Journal*, pages 45–50, 04 2018.

[4] S. Y. Caliskan and P. Tabuada. Compositional transient stability analysis of multimachine power networks. *Transactions on Control of Network Systems*, 1(1):4–14, 2014.

[5] S. Y. Caliskan and P. Tabuada. Correction to “Compositional transient stability analysis of multimachine power networks”. *Transactions on Control of Network Systems*, 4(3):676–677, 2017.

[6] S. Chowdhury, S. P. Chowdhury, and P. Crossley. Microgrids and active distribution networks. *Institution of Engineering and Technology*, 2009.

[7] F. Dörfler and F. Bullo. Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators. *SIAM Journal on Control and Optimization*, 50(3):1616–1642, 2012.

[8] D. Efimov and J. Schiffer. On boundedness of solutions of state periodic systems: A multivariable cell structure approach. *IEEE Transactions on Automatic Control*, 64(10):4094–4104, 2019.

[9] D. Efimov, J. Schiffer, N. Barabanov, and R. Ortega. A relaxed characterization of ISS for periodic systems with multiple invariant sets. *European Journal of Control*, 37:1–7, 2017.

[10] H. Farhangi. The path of the smart grid. *Power and Energy Magazine*, 8(1):18–28, 2010.

[11] S. Fiaz, D. Zonetti, R. Ortega, J. M. Scherpen, and A. Van der Schaft. A port-hamiltonian approach to power network modeling and analysis. *European Journal of Control*, 19(6):477–485, 2013.

[12] D. Groß, C. Arghir, and F. Dörfler. On the steady-state behavior of a nonlinear power system model. *Automatica*, 90:248–254, 2018.

[13] N. Hatziar gyriou, J. Milanovic, C. Rahmann, V. Ajjarapu, C. Canizares, I. Erlich, D. Hill, I. Hiskens, I. Kamwa, B. Pal, et al. Definition and classification of power system stability–revisited & extended. *Transactions on Power Systems*, 36(4):3271–3281, 2020.

[14] A. Krishna, I. Jaramillo-Cajica, S. Auer, and J. Schiffer. A power-hardware-in-the-loop testbed for intelligent operation and control of low-inertia power systems. *Automatisierungstechnik*, 70(12):1084–1095, 2022.

[15] G. A. Leonov. On the boundedness of the trajectories of phase systems. *Siberian Mathematical Journal*, 16(3):491–495, 1974.

[16] J. Machowski, J. W. Bialek, and J. R. Bumby. *Power System Dynamics: Stability and Control*. NY: Wiley, New York, 2008.

[17] A. Mercado-Uribe, J. Mendoza-Ávila, D. Efimov, and J. Schiffer. A control Leonov function guaranteeing global ISS of two coupled synchronverters. In *IEEE Conference on Decision and Control*, pages 4580–4585, 2023.

[18] F. Milano, F. Dörfler, G. Hug, D. J. Hill, and G. Verbič. Foundations and challenges of low-inertia systems. In *Power Systems Computation Conference*, pages 1–25. IEEE, 2018.

[19] V. Natarajan and G. Weiss. Synchronverters with better stability due to virtual inductors, virtual capacitors, and anti-windup. *Transactions on Industrial Electronics*, 64(7):5994–6004, 2017.

[20] M. Paolone, T. Gaunt, X. Guillaud, M. Liserre, S. Meliopoulos, A. Monti, T. Van Cutsem, V. Vittal, and C. Vournas. Fundamentals of power systems modelling in the presence of converter-interfaced generation. *Electric Power Systems Research*, 189:106811, 2020.

[21] R. Rüdenberg. *Transient Performance of Electric Power Systems: Phenomena in Lumped Networks*. McGraw-Hill, New York, 1950.

[22] J. Schiffer and D. Efimov. Strong and weak Leonov functions for global boundedness of state periodic systems. *IEEE Transactions on Automatic Control*, 68(12):7958–7965, 2023.

[23] J. Schiffer, D. Efimov, and R. Ortega. Global synchronization analysis of droop-controlled microgrids—a multivariable cell structure approach. *Automatica*, 109:108550, 2019.

[24] J. Schiffer, D. Efimov, R. Ortega, and N. Barabanov. An input-to-state stability approach to verify almost global stability of a synchronous-machine-infinite-bus system. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 375(2100):20160304, 2017.

[25] J. Schiffer, D. Zonetti, R. Ortega, A. M. Stanković, T. Sezi, and J. Raisch. A survey on modeling of microgrids—from fundamental physics to phasors and voltage sources. *Automatica*, 74:135–150, 2016.

[26] E. Venezian and G. Weiss. A warning about the use of reduced models of synchronous generators. In *International Conference on the Science of Electrical Engineering*, pages 1–5. IEEE, 2016.

[27] V. Venkatasubramanian, H. Schattler, and J. Zaborszky. Fast time-varying phasor analysis in the balanced three-phase large electric power system. *IEEE Transactions on Automatic Control*, 40(11):1975–1982, 1995.

[28] G. Weiss and E. Venezian. Stability analysis for coupled synchronous generators with virtual friction. In *Conference on Digital Signal Processing*, pages 1–5, 2017.

[29] W. Winter, K. Elkington, G. Bareux, and J. Kostevc. Pushing the limits: Europe's new grid: Innovative tools to combat transmission bottlenecks and reduced inertia. *Power and Energy Magazine*, 13(1):60–74, 2015.

[30] Q.-C. Zhong and G. Weiss. Synchronverters: Inverters that mimic synchronous generators. *Transactions on Industrial Electronics*, 58(4):1259–1267, 2011.