Combining Trajectory Data with Analytical Lyapunov Functions for Improved Region of Attraction Estimation

Stability analysis is of uttermost importance for the secure planning and operation of modern power systems. Of particular interest is the Transient Angular Stability of a system, defined as the ability of the system to maintain rotor angle synchronism following a disturbance and its subsequent angle excursion [1]. Given a Stable Equilibrium Point (SEP), the Region of Attraction (ROA), defined as the set of initial conditions for which the system tends to the SEP, provides a metric of the strength of angular stability of the system. Moreover, knowledge of the ROA can be used to provide protection parameters and limits of operation that maintain the stability and safety of the system.

For general nonlinear systems, which is the case for generators connected to the grid, there does not exist an analytical expression for the ROA. In the absence of an analytical expression, there is a need for methods that can compute approximations of the ROA. Classically, methods that approximate the ROA in power systems rely on precise system models of generators and line admittances [2, 3]. However, the increasing penetration of IBRs has resulted in more complex machine models and more dynamic operating points. For such complex systems it has become increasingly intractable to use classical model-based methods to accurately approximate the ROA [4]. Fortunately, the advent of Wide Area Measurement Systems (WAMS), gathering high-frequency synchrophasor data has provided new sets of system data with minimal model dependency. Some works have already explored the use of synchrophasors for ROA estimation in the literature.

In [5], data driven ROA estimation is realized through the application of energy function analysis, using PMU measurements that monitor tie-lines of dynamic power flows. Authors from [6] propose an alternative method that uses Global Phase Portraits (GPPs) that contain the singularity points at infinity, providing bounds on the basins of attraction of attractor sets. In a similar vein to the data-driven methods proposed in this paper, authors from [7, 8, 9] use measurement data from stable trajectories to approximate converse Maximal Lyapunov Functions (LFs) and hence construct ROA estimations. Although these methodologies bring new capabilities for high dimensional stability analysis, these methods do not guarantee an inner approximation of the ROA, unlike more classical model-based methods. Notable model based methods include [10] where the stability analysis of power systems is analyzed by constructing LFs using Sum-of-Squares (SOS) programming. Since power systems have nonlinear trigonometric terms, non-automated algebraic reconfiguration is required to use SOS. Alternatively, the works of [6, 4] make analytical approaches that improve upon classical energy function based methods.

Unfortunately, it is often intractable to compute accurate ROA estimations of power systems using model based methods. On the other hand, although data based methods can provide accurate ROA estimations, they do not yield LFs and hence cannot certify inner ROA approximations. The goal of this work is to bridge the gap between the model and data based methods to yield accurate inner ROA approximations.

The main contribution of this paper, presented in Thm. 3, shows how the existence of two functions, $V_{1}$ and $V_{2}$, provides a certifiable inner approximation of the ROA of a given ODE. Specifically, if $V_{2}$ is a LF and $V_{1}$ (not necessarily a LF) is decreasing along the solution map inside a “donut”-shaped region, $\{\hskip-1.42271ptx\hskip-1.42271pt\in\hskip-1.42271ptD\hskip-1.42271pt:\hskip-1.42271pt\gamma_{1}\hskip-2.84544pt\leq\hskip-1.42271ptV_{1}(x)\hskip-2.84544pt\leq\hskip-1.42271pt\gamma_{2}\hskip-1.42271pt\}$, we show that it is possible to construct an improved ROA estimation, as compared with the ROA approximation yielded by the LF, $V_{2}$, alone. For implementation, we find such a $V_{1}$ by function fitting a converse LF using trajectory data.

We denote the $\eta>0$ neighborhood of a set $S\subset\mathbb{R}^{n}$ as $B_{\eta}(S):=\{y\in\mathbb{R}^{n}:\inf_{x\in S}||x-y||_{2}<\eta\}$, where $||\cdot||_{2}$ is the euclidean norm. In the case $S=\{x\}$ then $B_{\eta}(x)$ becomes a ball of radius $\eta>0$ centered at $x\in\mathbb{R}^{n}$. We denote the set of all interior points of $S\subset\mathbb{R}^{n}$ by $S^{\circ}$. Let $C(\Omega,\Theta)$ be the set of continuous functions with domain $\Omega\subset\mathbb{R}^{n}$ and image $\Theta\subset\mathbb{R}^{m}$. For $\alpha\in\mathbb{N}^{n}$ we denote the partial derivative $D^{\alpha}:=\Pi_{i=1}^{n}\frac{\partial^{\alpha_{i}}}{\partial x_{i}^{\alpha_{i}}}$ where by convention if $\alpha=[0,..,0]^{\top}$ we denote $D^{\alpha}f(x):=f(x)$ for any function $f$. We denote the set of $i$’th continuously differentiable functions by $C^{i}(\Omega,\Theta):=\{f\in C(\Omega,\Theta):D^{\alpha}f\in C(\Omega,\Theta)\text{ }\text{ for all }\alpha\in\mathbb{N}^{n}\text{ such that }\sum_{j=1}^{n}\alpha_{j}\leq i\}$. For $V\in C^{1}(\mathbb{R}^{n},\mathbb{R})$ we denote $\nabla V:=(\frac{\partial V}{\partial x_{1}},....,\frac{\partial V}{\partial x_{n}})^{\top}$. We denote the space of $d$-degree polynomials $p:\Omega\to\Theta$ by $\mathcal{P}_{d}(\Omega,\Theta)$.

Consider a dynamical system, represented by a nonlinear ordinary differential equation (ODE) of the form

where $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the vector field and $x_{0}\in\mathbb{R}^{n}$ is the initial condition. WLOG throughout this paper we will assume $f(0)=0$ so the origin is an equilibrium point; a linear coordinate transformation can always be used to shift any equilibrium point to the origin.

For simplicity in the following we assume Eq. (1) is well defined. That is there exists a unique solution map $\phi_{f}\in C^{1}(\mathbb{R}^{n}\times\mathbb{R}^{+},\mathbb{R}^{n})$ that satisfies $\frac{\delta\phi_{f}(x,t)}{\delta t}=f(\phi_{f}(x,t))$, $\phi_{f}(x,0)=x$ and $\phi_{f}(\phi_{f}(x,t),s)=\phi_{f}(x,t+s)$. Sufficient conditions for the existence and uniqueness of a solution map, based on the smoothness properties of the vector field, can be found in standard textbooks such as [11]. Given an ODE (1), we next introduce notions of asymptotic and exponential stability that are important in showing the existence of the converse LF given later in Eq (3).

For a given asymptotically stable ODE (1) the main aim of this paper is to estimate the Region of Attraction (ROA):

There is no universal method for analytically solving nonlinear ODEs. Thus, over the years, arguably the most commonly used method to estimate $ROA_{f}$ is Lyapunov’s second method that indirectly estimates $ROA_{f}$ using Lyapunov Functions (LFs); functions that are globally non-negative that decrease along the solution map. The following theorem shows how the sublevel set of a LF can approximate $ROA_{f}$. In order to present the main Lyapunov theorem used in this paper we recall the definition of an invariant set.

Theorem 1 shows that for a given ODE, if we can find a LF, then we can construct an inner-approximate of the ROA of the ODE. However, this theorem does not show that there must necessarily exists a LF for a given ODE or that the ROA of the ODE can be exactly characterized by an LF.

It has been shown in [12] that for any locally exponentially stable ODE (1), there exists a bounded and continuous LF of the form,

where $\lambda>0$ and $\beta\in\mathbb{N}$. Moreover, for sufficiently large $\lambda$ and $\beta$, this converse LF, $V^{*}_{\lambda,\beta}$, is Lipschitz continuous and hence differentiable almost everywhere (by Rademacher’s theorem). The smoothness properties of this particular converse LF makes it highly suitable for function fitting. Furthermore, $\{x\in\mathbb{R}^{n}:V^{*}_{\lambda,\beta}(x)<1\}=ROA_{f}$.

Given an ODE (1), in this section we show that by using an ODE solver to generate trajectory data, Eq. (3) can be used to construct input-output data of a converse LF. By fitting a function to this data we can approximate this converse LF in the hope of constructing an ROA estimation.

Specifically, we fit polynomial functions to the generated input/output data of the converse LF. Although there are many ways to fit polynomials to data (each having their relative advantages and disadvantages), in this paper, we have chosen a method based on Bernstein approximations. As we will next see, Bernstein’s method for fitting polynomials to data is an optimization free approach that is guaranteed to converge uniformly.

Given an ODE defined by some vector field $f$, we have shown that through applications of Eqs. (4) and (7) that it is possible to numerically a construct a Bernstein polynomial approximation, $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$ for some $d\in\mathbb{N}$, $\lambda>0$ and $\beta\in\mathbb{N}$, of the converse LF, $V^{*}_{\lambda,\beta}$, given in Eq. (3). Moreover, assuming that $D^{\alpha}V^{*}_{\lambda,\beta}$ is continuous, where $\alpha\in\mathbb{N}^{n}$, Theorem 2 can be used to show that $\lim_{d\to\infty}D^{\alpha}\mathcal{B}_{d}V^{*}_{\lambda,\beta}\to D^{\alpha}V^{*}_{\lambda,\beta}$ uniformly in $[a,b]^{n}$ (note approximation over $[a,b]^{n}$ rather than $[0,1]^{n}$ can be achieved through linear coordinate transformations, see Eq. (6)).

Ideally, for sufficiently large $d\in\mathbb{N}$, $\lambda>0$ and $\beta\in\mathbb{N}$, our approximation, $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$, will also be a LF over some set containing the origin; that is $\mathcal{B}_{d}V^{*}_{\lambda,\beta}(0)=0$, $\mathcal{B}_{d}V^{*}_{\lambda,\beta}(x)\geq 0$ for all $x\in\Omega/\{0\}$, and $\nabla\mathcal{B}_{d}V^{*}_{\lambda,\beta}(x)^{\top}f(x)<0$ for all $x\in\Omega/\{0\}$. Then Theorem 1 can be used to show that the sublevel set of $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$ yields an inner approximation of the ROA of the given ODE.

Unfortunately, despite the fact $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$ tends to $V^{*}_{\lambda,\beta}$ as $d\to\infty$ and $V^{*}_{\lambda,\beta}$ is a LF, $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$ is not necessarily a LF. To see this we first note that, for a given $x_{0}\in[a,b]^{n}$ if $\nabla V^{*}_{\lambda,\beta}(x_{0})^{\top}f(x_{0})<-a$, where $a>0$, it follows by Theorem 2 that there exists $D\in\mathbb{N}$ such that $||\nabla\mathcal{B}_{d}V^{*}_{\lambda,\beta}(x_{0})-\nabla V^{*}_{\lambda,\beta}(x_{0})||_{2}<\frac{a}{2||f(x_{0})||_{2}}$ for all $d>D$. Thus, using the Cauchy Swarz inequality we have that for all $d>D$

Eq. (8) shows that for sufficiently large $d$, whenever $V^{*}_{\lambda,\beta}$ is strictly decreasing along the solution map we also have that $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$ is strictly decreasing along the solution map. However, at the origin $V^{*}_{\lambda,\beta}$ is not strictly decreasing along the solution map, that is $\nabla V^{*}_{\lambda,\beta}(0)^{\top}f(0)=0$. Because of this fact and the fact $\mathcal{B}_{d}V^{*}_{\lambda,\beta}$ is continuous, in general for some finite $d\in\mathbb{N}$ our approximation will be such that that $\nabla\mathcal{B}_{d}V^{*}_{\lambda,\beta}(x)^{\top}f(x)\geq 0$ for all $x$ in some small neighborhood of the origin. Thus in general, for finite $d\in\mathbb{N}$, it follows that we will have $B_{d}V^{*}_{\lambda,\beta}(x)^{\top}f(x)<0$ for all $x$ in some “donut shaped” region, $\{y\in\mathbb{R}^{n}:\gamma_{1}\leq B_{d}V^{*}_{\lambda,\beta}(y)\leq\gamma_{2}\}$ for some $\gamma_{1}<\gamma_{2}$, as opposed to a sublevel set $\{y\in\mathbb{R}^{n}:B_{d}V^{*}(y)\leq\gamma_{2}\}$. By a similar argument, in general, we do not expect $B_{d}V^{*}_{\lambda,\beta}(0)=0$ for any fixed $d\in\mathbb{N}$ and thus in general $B_{d}V^{*}_{\lambda,\beta}$ is not a LF.

Although, $B_{d}V^{*}_{\lambda,\beta}$ is not a LF, and therefore cannot certify the origin is asymptotically stable, we will next show, in Prop. 1, that functions strictly decreasing along the solution map inside some “donut shaped” region can still be used to certify that the solution map must enter some ball of radius $\eta>0$ around the origin, $B_{\eta}(0)$.

In Prop. 1 we have proposed conditions, based on some function we denote here as $V_{1}$, that certify that the solution map of a given ODE must enter some ball, $B_{\eta}(0)$. Recall that $V_{1}$ can be found by approximating a converse LF by Bernstein polynomials.

We next show that if there exists a LF, $V_{2}$, that certifies that $B_{\eta}(0)$ is an asymptotically stable set, then $V_{1}$ and $V_{2}$ can be used together to provide an improved inner approximation of the ROA of the given ODE.

Practical Implementation of Theorem 3
Thm. 3 shows how two functions, $V_{1}$ and $V_{2}$, can be used to produce an inner approximation of the ROA of a given ODE. The function, $V_{2}$, is a classical LF (positive semidefinite and decreasing along trajectories) but provides a poor estimation of the ROA, $\{y\in D:V_{2}(y)\leq\gamma_{3}\}\subseteq ROA_{f}$. On the other hand, $V_{1}$ need not be a classical LF as there is no requirement that $V_{2}(0)=0$ or $\nabla V_{2}(x)^{\top}f(x)\leq 0$ near the origin, making the computational search for a candidate $V_{1}$ amenable to data based methods. Whereas, the computational search for a candidate $V_{2}$ is redistricted to model based methods (energy functions or Sum-of-Squares programming).

We next present several numerical examples demonstrating how Thm. 3 can be used to yield accurate ROA approximations of nonlinear ODEs. For each numerical example we compute $V_{1}$ and $V_{2}$, from Thm. 3, by respectively fitting Bernstein polynomials to the converse LF given in Eq. (3) (using Eqs. (4), (6) and (7)) and either using an energy function or an analytical LF (found previously in the literature).

To demonstrate the accuracy of our approximations of the ROA we will carry out extensive Monte Carlo simulations of the solution map to estimate the stable and unstable regions in each figure. Although this Monte Carlo method can estimate the ROA well it does not account for simulation error or provide a LF and hence cannot provide a certified ROA inner approximation.