Contraction Analysis of Nonlinear DAE Systems

TL;DR

This paper investigates the contraction properties of nonlinear DAE systems, providing scalable methods to determine attraction regions and analyze stability, with applications demonstrated in power system transient stability.

Contribution

It introduces a novel scalable technique linking contraction rates of original and extended DAE systems, enhancing stability analysis methods.

Findings

Contracting DAE systems have faster reduced systems.

Existence of extensions with contraction rates close to original.

Application to power system transient stability assessment.

Abstract

This paper studies the contraction properties of nonlinear differential-algebraic equation (DAE) systems. Specifically we develop scalable techniques for constructing the attraction regions associated with a particular stable equilibrium, by establishing the relation between the contraction rates of the original systems and the corresponding virtual extended systems. We show that for a contracting DAE system, the reduced system always contracts faster than the extended ones; furthermore, there always exists an extension with contraction rate arbitrarily close to that of the original system. The proposed construction technique is illustrated with a power system example in the context of transient stability assessment.

Tables

Table 1. TABLE I: Standard matrix measures

Vector norm, $\parallel \cdot \parallel$	Matrix measure, ${\mu }_p\left(M\right)$
${\parallel x\parallel }_1={\sum }_i\left|x_i\right|$	${\mu }_1\left(M\right)={\max }_j\left(m_{jj}+{\sum }_{i\ne j}\left|m_{ij}\right|\right)$
${\parallel x\parallel }_2={\left({\sum }_i{\left|x_i\right|}^2\right)}^{1/2}$	${\mu }_2\left(M\right)={\max }_i\left({\lambda }_i\left\{\frac{M+M^T}{2}\right\}\right)$
${\parallel x\parallel }_{\infty }={\max }_i\left|x_i\right|$	${\mu }_{\infty }\left(M\right)={\max }_i\left(m_{ii}+{\sum }_{j\ne i}\left|m_{ij}\right|\right)$

Equations

$$\dot{\mathbf{x}}$$

$$0$$

$$E\dot{\mathbf{z}}=h(\mathbf{z}),$$

$$\mathcal{R}=\left\{\mathbf{x}:\,\,\mathrm{det}\left(\frac{\partial\mathbf{g}}{\partial\mathbf{y}}\Big{|}_{\mathbf{y}=\mathbf{Y}(\mathbf{x})}\right)\neq 0\right\}.$$

$$J(\mathbf{x},\mathbf{y})=\left[\begin{array}{cc}\partial \mathbf{f}/\partial \mathbf{x}& \partial \mathbf{f}/\partial \mathbf{y}\\ \partial \mathbf{g}/\partial \mathbf{x}& \partial \mathbf{g}/\partial \mathbf{y}\end{array}\right].$$

$$J(\mathbf{x},\mathbf{Y}(\mathbf{x}))=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right].$$

$$F_r$$

$$\dot{\theta }={\left(\frac{\partial \theta }{\partial \mathbf{x}}\right)}^T\mathbf{f}(\mathbf{x},\mathbf{Y}(\mathbf{x})).$$

$$\mathbf{x}(t)\in {\mathcal{I}}_p⟹\forall t^{\prime }\ge t:\mathbf{x}(t^{\prime })\in {\mathcal{I}}_p.$$

$$F$$

$$\dot{\delta \mathbf{v}}$$

$$\dot{\delta \mathbf{u}}$$

$$F\delta \mathbf{w}$$

$$=\left[\begin{array}{c}{F}_r\theta \delta \mathbf{x}\\ 0\end{array}\right]+\left[\begin{array}{c}\theta R^T(C\delta \mathbf{x}+D\delta \mathbf{y})\\ Q^T(C\delta \mathbf{x}+D\delta \mathbf{y})\end{array}\right].$$

$$\displaystyle\gamma=\biggl{\|}\begin{bmatrix}\delta\mathbf{v}+hF_{r}\delta\mathbf{v}\\ \delta\mathbf{u}\end{bmatrix}\biggr{\|}_{p}-\biggl{\|}\begin{bmatrix}\delta\mathbf{v}\\ \delta\mathbf{u}\end{bmatrix}\biggr{\|}_{p}-\left(\|\delta\mathbf{v}+hF_{r}\delta\mathbf{v}\|_{p}-\|\delta\mathbf{v}\|_{p}\right)$$

$$\biggl{\|}\begin{bmatrix}\delta\mathbf{v}+hF_{r}\delta\mathbf{v}\\ \delta\mathbf{u}\end{bmatrix}\biggr{\|}_{p}-\biggl{\|}\begin{bmatrix}\delta\mathbf{v}\\ \delta\mathbf{u}\end{bmatrix}\biggr{\|}_{p}\leq 0.$$

$$\frac{\partial \gamma }{\partial |\delta {\mathbf{u}}_j|}=|\delta {\mathbf{u}}_j{|}^{p-1}$$

$$\displaystyle-\biggl{\|}\begin{bmatrix}\delta\mathbf{v}\\ \delta\mathbf{u}\end{bmatrix}\biggr{\|}_{p}^{p(1/p-1)}\biggr{)}.$$

$${\mu }_p(F_r)\le {\mu }_p(F){\nu }_p(H).$$

$${\nu }_p(H)=\underset{{\Vert v\Vert }_p=1}{\min }{\Vert Hv\Vert }_p.$$

$$\frac{{\Vert \delta {\mathbf{v}}_h+h{F}_r\delta {\mathbf{v}}_h\Vert }_p-{\Vert \delta {\mathbf{v}}_h\Vert }_p}{h{\Vert \delta {\mathbf{v}}_h\Vert }_p}=\frac{{\Vert 1+h{F}_r\Vert }_p-1}{h}.$$

$${\mu }_p(F_r)=\underset{h\to 0}{\lim }\frac{\Vert \delta {\mathbf{v}}_h+h{F}_r\delta {\mathbf{v}}_h{\Vert }_p-\Vert \delta {\mathbf{v}}_h{\Vert }_p}{h{\Vert \delta {\mathbf{v}}_h\Vert }_p}.$$

$$\mu_{p}(F_{r})\leq\lim\limits_{h\rightarrow 0}\frac{\biggl{\|}\delta\mathbf{w}_{h}+h\begin{bmatrix}F_{r}\delta\mathbf{v}_{h}\\ 0\end{bmatrix}\biggr{\|}_{p}-\|\delta\mathbf{w}_{h}\|_{p}}{h\left\lVert\delta\mathbf{v}_{h}\right\rVert_{p}}.$$

$$\delta {\mathbf{w}}_h+h\left[\begin{array}{c}{F}_r\delta {\mathbf{v}}_h\\ 0\end{array}\right]$$

$$=\delta {\mathbf{w}}_h+hF\delta {\mathbf{w}}_h-h\left[\begin{array}{c}\theta R^T(C\delta {\mathbf{x}}_h+D{\widehat{\delta \mathbf{y}}}_h)\\ Q^T(C\delta {\mathbf{x}}_h+D{\widehat{\delta \mathbf{y}}}_h)\end{array}\right]$$

$$=\delta {\mathbf{w}}_h+hF\delta {\mathbf{w}}_h$$

$${\mu }_p(F_r)$$

$$=\underset{h\to 0}{\lim }\frac{\Vert \delta {\mathbf{w}}_h+hF\delta {\mathbf{w}}_h{\Vert }_p-\Vert \delta {\mathbf{w}}_h{\Vert }_p}{h{\Vert \delta {\mathbf{w}}_h\Vert }_p}\frac{{\Vert \delta {\mathbf{w}}_h\Vert }_p}{{\Vert \delta {\mathbf{v}}_h\Vert }_p}.$$

$$\delta {\mathbf{w}}_h=\left[\begin{array}{c}\delta {\mathbf{v}}_h\\ \delta {\mathbf{u}}_h\end{array}\right]=H\delta {\mathbf{v}}_h,$$

$${\mu }_p(F_r)$$

$$\le {\mu }_p(F){\nu }_p(H).$$

$${\mu }_2(F_r)\le {\mu }_2(F)\sqrt{1+{\sigma }_{\min }^2(S)}$$

$$\delta {\mathbf{w}}^{T}F\delta \mathbf{w}={\left[\begin{array}{c}\theta \delta \mathbf{x}\\ \rho \delta \mathbf{y}\end{array}\right]}^T\left[\begin{array}{cc}{F}_r+\eta S^{T}S& 0\\ 0& -\eta 1\end{array}\right]\left[\begin{array}{c}\theta \delta \mathbf{x}\\ \rho \delta \mathbf{y}\end{array}\right]$$

$$=\delta {\mathbf{v}}^T\left(F_r+\eta S^{T}S\right)\delta \mathbf{v}-\eta \Vert \delta \mathbf{u}{\Vert }^2$$

$$\le \left({\mu }_2(F_r)+\eta {\sigma }_{\max }^2(S)\right)\Vert \delta \mathbf{v}{\Vert }^2-\eta \Vert \delta \mathbf{u}{\Vert }^2$$

Full text

I Introduction

Differential-algebraic equations (DAE)-a generalization of ordinary-differential equations (ODE)-arise in many science and engineering problems, including networks, multibodies, optimal control, compressed fluid, etc. [1]. Typically, algebraic constraints result from multiple time-scale perturbation theory, when the fast degrees of freedom are assumed to stay on equilibrium manifold. In typical electrical and mechanical applications the algebraic relations represent the interconnection constraints, which can be considered static on the time-scales of system evolution. However, algebraic relations may be also useful for lifted representations of the purely differential systems. For instance, additional variables and relations can be used to represent any polynomial nonlinearity in a quadratic DAE form. Hence, DAE systems provide a powerful framework for studying nonlinear systems of very general structure. This work is motivated by the DAE representations of the power system models, but the results are presented in a general form.

The specific problem that motivates our study is the problem of approximating the region of attraction of DAE equilibrium points. The normal operating points of modern power systems lack global stability because of the nonlinearities naturally appearing in these systems. Characterization of the attraction region and more generally assessment of the system security, i.e. its ability to sustain all kinds of faults and disturbances, is an essential task of modern power system operations. As will be shown throughout the paper, the contraction provides a natural framework for constructing the approximations of the attraction region for a broad range of nonlinear DAE problems, such as those arising in power systems.

Transient stability analysis is a common engineering procedure referring to the ability of the system to converge to a stable post-fault equilibrium after being subject to disturbances. The incremental stability introduced in [2] suggests an alternative way to look at the convergence of the post-fault trajectories. In the light of contraction theory, the virtual displacements of the states tend to zero as the time goes to infinity, or in other words, all the trajectories shrink and converge to the nominal one. Contraction analysis becomes a powerful tool for nonlinear analysis and control [3, 2, 4, 5]. The key property of the contraction is the preservation under different system combinations, which is advantageous in network analysis.

In this paper we focus on the contraction analysis for nonlinear DAE systems. Specifically we develop a practical way of constructing the attraction regions by determining the relation between the contraction rates of the original DAE systems and its extension to virtual dynamics in differential-algebraic space. The extended system can be thought of as a virtual differential system that reduces to a given DAE after the restriction of a subset of variables to their equilibrium manifold. There can be multiple extensions of a given DAE system, each characterized by different contraction rates. However, we show that the contraction rate of the reduced system is always higher, and on the other hand, there always exists an extension with a contraction rate arbitrarily close to the original DAE system. Our results hold for the most commonly used $1$, $2$, and $\infty$ norms, but can likely be extended to more general cases. We use the theoretical results to develop a scalable technique for constructing ellipsoidal inner approximations of contraction regions from the $2$ norm contraction metric. We illustrate the technique with a practical example from power systems.

II Main results

As motivated by the dynamics of electrical power systems, we constrain ourselves to semi-explicit index $1$ structural form as below:

[TABLE]

In this representation, vector $\mathbf{x}\in\mathbb{R}^{n}$ corresponds to dynamic state variables, $\mathbf{y}\in\mathbb{R}^{m}$ refers to algebraic variables (whose dynamics is assumed to be fast/instantaneous relative to the dynamics of the state variables). For this class of systems, it is impossible to obtain equivalent ODEs.

For convenience, reduction techniques are widely used to eliminate the algebraic variables. Yet this practice may prohibit one from exploring the underlying structure of the DAE form. To that end a number of works in the literature concentrate on the original systems rather than the reduced ones, for instant, in the context of stability analysis of the descriptor form as below:

[TABLE]

with $\mathbf{z}^{T}=[\mathbf{x}^{T},\mathbf{y}^{T}]$, $h^{T}=[f^{T},g^{T}]$ and $E$ being a diagonal $\mathbb{R}^{(n+m)\times(n+m)}$ matrix with $E_{ii}=1$ for $i\leq n$ and $E_{ii}=0$ otherwise [6, 7, 8, 9].

For any given differential state $\mathbf{x}$ the equation (2) may have multiple or no solutions for $\mathbf{y}$. In engineering and natural systems that motivate this study, disappearance of all the solutions is usually an indicator of inappropriate modeling that should be fixed accordingly, typically by introducing the fast dynamics of the algebraic states in the model. We don’t consider this scenario in our work, and we assume that for every $\mathbf{x}$ there exists at least one solution $\mathbf{Y}(\mathbf{x})$ of the algebraic system of equations (2). For every solution branch we can naturally define the domain $\mathbf{x}\in\mathcal{R}$ where such a solution exists and can be tracked via homotopy/continuation procedure. This domain is characterized by non-singularity of the algebraic Jacobian:

[TABLE]

We restrict our analysis only to such a domain associated with a specific solution branch. For a system of differential-algebraic equations (1), (2) we introduce the Jacobian defined as

[TABLE]

To simplify the notations we also define its restriction to the algebraic manifold (2) as follows:

[TABLE]

One of the primary goals of this paper is to provide a characterization of the contraction and invariant regions in the state space of a DAE system. We formally define the contraction domains $\mathcal{C}_{p}$ as set of differential states $\mathbf{x}$ for which there exists an invertible metric $\theta(\mathbf{x})\in\mathbb{R}^{n\times n}$ such that the differential equation $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{Y}(\mathbf{x}))$ is locally contracting with respect to this metric with some rate $\beta>0$. Given that for any infinitesimal displacement $\delta\mathbf{x}$ we have $\delta\mathbf{y}=-D^{-1}C\delta\mathbf{x}$ the standard contraction arguments presented in [4, 2] lead to the following Proposition.

Proposition 1

The DAE system (1) (2) is contracting with respect to the metric $\theta(\mathbf{x})$ in the domain $\mathcal{C}_{p}$ if for all $\mathbf{x}\in\mathcal{C}_{p}$ one has $\mu_{p}\left(F_{r}\right)\leq-\beta$ with some $\beta>0$ and

[TABLE]

The term $\dot{\theta}$ in (7) represents the derivative of the metric along the trajectory and is formally defined for DAE systems as

[TABLE]

The matrix measure $\mu_{p}(M)$ of a matrix $M$ is defined as $\mu_{p}(M):=\lim\limits_{h\rightarrow 0^{+}}\frac{1}{h}(||\mathbb{1}+hM||_{p}-1)$ following [10]. The standard matrix measures as well as vector norms are listed in Table I.

The proof of proposition 1 directly follows from the contraction analysis for dynamical system presented in [2]. The matrix $F_{r}$ appears naturally from the dynamic equation on $\delta\mathbf{v}=\theta\delta\mathbf{x}$ given by $\dot{\delta\mathbf{v}}=F_{r}\delta\mathbf{v}$. Hereafter we refer to $F_{r}$ as the generalized reduced Jacobian matrix.

The standard contraction theory arguments suggest that for any two trajectories $\mathbf{x}_{1}(t),\mathbf{x}_{2}(t)$ that both remain within the contraction region $\mathcal{C}_{p}$ during the interval $[t_{1},t_{2}]$ satisfy $d(\mathbf{x}_{1}(t_{2}),\mathbf{x}_{2}(t_{2}))\leq d(\mathbf{x}_{1}(t_{1}),\mathbf{x}_{2}(t_{1}))\exp(-\beta(t_{2}-t_{1}))$ where $d$ is the distance associated with the metric $\theta$. The assumption that both of the trajectories stay within the contraction region is critical for this result and can be verified only after showing the existence of an invariant domain $\mathcal{I}_{p}\subset\mathcal{C}_{p}$ satisfying:

[TABLE]

Constructing invariant regions is usually a difficult aspect of applying contraction theory to systems which are not globally contracting. One straightforward strategy for constructing invariant regions exists for systems that have an equilibrium point $\mathbf{x}_{\star}$ inside the contraction domain satisfying $\mathbf{f}(\mathbf{x}_{\star},\mathbf{Y}(\mathbf{x}_{\star}))=0$. In this case, any ball $\mathcal{B}_{r}=\{\mathbf{x}:d(\mathbf{x},\mathbf{x}_{\star})\leq r\}$ that lies within the contraction region $\mathcal{C}_{p}$ defines an invariant region, i.e. $\mathcal{B}_{r}\subset\mathcal{C}_{p}\Longrightarrow\mathcal{B}_{r}\subset\mathcal{I}_{p}$. By construction, such a ball also provides an inner approximation for the attraction region of $\mathbf{x}_{\star}$ and can be naturally used in a variety of practical applications such as security assessment of power systems [11]. In this we develop a general framework for constructing such invariant regions for a broad class of nonlinearities, and we present a specific power system example in sections III.

The key challenge in using the function $F_{r}$ directly is its highly nonlinear nature. Even for simple polynomial nonlinearities of $\mathbf{f},\mathbf{g}$ the function $F_{r}$ involves an inversion of the matrix $D(\mathbf{x})$. From a practical perspective, it is therefore desirable to formulate conditions equivalent to contraction as defined in Propostion 1 that do not involve any inversions of matrices $A,B,C,D$ which are nonlinearly dependent on $\mathbf{x}$. In order to achieve this goal we derive equivalent representation of the contraction condition that doesn’t require elimination of the local variables and is more suitable for analysis. We introduce the generalized unreduced Jacobian matrix as follows:

[TABLE]

The generalized unreduced Jacobian $F$ depends on the metric $\theta$ defined as in the previous discussion, another metric $\rho$ associated with the $\mathbf{y}$ variable and two auxiliary matrices $Q\in\mathbb{R}^{m\times m}$ and $R\in\mathbb{R}^{m\times n}$. Formally, this Jacobian matrix may be associated with a virtual extended ODE representation of the original system of the form

[TABLE]

where $\delta\mathbf{u}=\rho\delta\mathbf{y}$ and so the expression $C\theta^{-1}\delta\mathbf{v}+D\rho^{-1}\delta\mathbf{u}=C\delta\mathbf{x}+D\delta\mathbf{y}=0$ defines the algebraic manifold. Whenever the dynamics of $\delta\mathbf{u}$ can be considered fast, the restriction of the $\delta\mathbf{u}$ variables to their equilibrium manifold results in the original DAE systems. Therefore, this representation provides a family of extended representations that reduce to the same original system. It will be shown that this representation is useful for characterization of the contraction and invariant regions.

The key property important for the analysis is defined in the following relation:

[TABLE]

where we have introduced the new variables vector $\delta\mathbf{w}\triangleq\begin{bmatrix}\delta\mathbf{v}\\ \delta\mathbf{u}\end{bmatrix}$. This observation allows us to formulate the following central results of this work.

II-A Forward theorems: from extended systems to reduced ones

Lemma 1

Define

[TABLE]

where $h>0$. Then for all $p\geq 1$, $\gamma\geq 0$ if the following condition holds.

[TABLE]

The proof of Lemma 1 is as the following. For fixed $\delta\mathbf{v}$ and $hF_{r}\delta\mathbf{v}$, $\gamma$ depends solely on $\delta\mathbf{u}$. Taking partial derivative of $\gamma$ with respect to $|\delta\mathbf{u}_{j}|$ using the definition of $p$ norm for a vector, i.e. $\|\mathbf{v}\|_{p}=\left(\sum_{i}{|v_{i}|^{p}}\right)^{1/p}$, yields the following:

[TABLE]

On the other hand, the assumption $p\geq 1$ leads to $p(1/p-1)\leq 0$. This together with (14) and (II-A) concludes $\frac{\partial\gamma}{\partial|\delta\mathbf{u}_{j}|}\geq 0$ for all $j=1,\dots,m$. In other words, $\gamma$ is indeed a monotonically increasing function with respect to the absolute value of each entry $\delta\mathbf{u}_{j}$. Moreover it can be seen that $\gamma$ vanishes when $\delta\mathbf{u}=0$. Therefore for any non-zero $\delta\mathbf{u}$, $\gamma$ is non-negative.

For infinity norm one can prove the non-negativity of the partial derivatives $\frac{\partial\gamma}{\partial|\delta\mathbf{u}_{j}|}$ by taking the limit as $p$ goes to infinity and exploiting the assumption presented by (14). Alternatively $\gamma$ can be directly evaluated using the matrix measure expression associated with infinity norm listed in Table I.

With Lemma 1 we introduce the first central result as the following.

Theorem 1

For the system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{Y}(\mathbf{x}))$ and metric function $\theta(\mathbf{x})$, and contracting extended system $F$ with $\mu_{p}(F)<0$ characterized by the matrices $Q,R,\rho$, the following relation holds:

[TABLE]

in which $S=\rho D^{-1}C\theta^{-1}$, $H=\begin{bmatrix}1\\ S\end{bmatrix}$, and

[TABLE]

Note that for invertible $H$, one has $\nu_{p}(H)=1/\left\lVert H^{-1}\right\rVert_{p}$.

Proof 1

The matrix induced norm definition $\left\lVert M\right\rVert=\max\limits_{\left\lVert v\right\rVert=1}\left\lVert Mv\right\rVert$ implies that for each positive scalar $h$, there exists $\delta\mathbf{v}_{h}$ with $\left\lVert\delta\mathbf{v}_{h}\right\rVert_{p}>0$ satisfying the following equality:

[TABLE]

The logarithmic norm is then defined as

[TABLE]

Lemma 1 implies that this expression can be also rewritten as

[TABLE]

On the other hand, applying the property defined by (II) we have that

[TABLE]

where $\delta\mathbf{x}_{h}=\theta^{-1}\delta\mathbf{v}_{h}$ and $\hat{\delta\mathbf{y}}_{h}=-D^{-1}C\delta\mathbf{x}_{h}$. Combining (20) and (1), we have that

[TABLE]

By definition

[TABLE]

so $\left\lVert\delta\mathbf{w}_{h}\right\rVert_{p}\geq\nu_{p}(H)\left\lVert\delta\mathbf{v}_{h}\right\rVert_{p}$, and by combining this with the assumption $\mu_{p}(F)<0$, we can rewrite (22) as the following

[TABLE]

Q.E.D.

The below corollary of Theorem 1 provides the explicit expression of $\nu_{p}(S)$ with $p=2$.

Corollary 1

($2$* norm**) Assuming that all assumptions of Theorem 1 are satisfied, the contraction rate associated with the reduced system can be bounded as below*

[TABLE]

in which $\sigma_{\min}(S)$ denotes the minimal singular value of the matrix $S$.

*The proof of Corollary 1 follows from Theorem 1 and note that for $p=2$, we have that $\nu_{2}(H)=\sqrt{1+\min\limits_{\left\lVert v\right\rVert=1}\|Sv\|_{2}^{2}}=\sqrt{1+\sigma_{\min}^{2}(S)}$. *

II-B Converse theorems

Theorem 2

($2$* norm**) For a contracting system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{Y}(\mathbf{x}))$ and metric function $\theta(\mathbf{x},t)$ with $\mu_{2}(F_{r})<0$ and any $\epsilon>0$ there exists an extended system $F$ characterized by the matrices $Q,R,\rho$ contracting with the contraction rate satisfying $\mu_{2}(F)\leq\mu_{2}(F_{r})/(1+\epsilon)$.*

Proof 2

This theorem can be proven by explicit construction of the matrices $Q,R,\rho$, which ensures fast enough contraction of $F$. The matrix $\rho$ is chosen to be small enough, so that $\sigma_{\max}^{2}(S)\leq\epsilon$ and $R=\eta D^{-T}\rho^{T}\rho D^{-1}CP^{-1}$, $Q=-\eta D^{-T}\rho^{T}$, where $\eta=-\mu_{2}(F_{r})/(1+\sigma_{\max}^{2}(S))$. This choice of $R$ and $Q$ ensures that the symmetric part of $F$ is block-diagonal, so the following inequality follows from (II):

[TABLE]

Since this inequality holds true for any $\delta{\mathbf{w}}$, we conclude that $\mu_{2}(F)\leq-\eta=\mu_{2}(F_{r})/(1+\epsilon)$.

The counterpart of Theorem 2 for $1$ norm and $\infty$ norm are presented below.

Theorem 3

($1$* norm and $\infty$ norm**) For a contracting system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{Y}(\mathbf{x}))$ and metric function $\theta(\mathbf{x},t)$ with $\mu_{p}(F_{r})<0$ and any $\epsilon>0$ there exists an extended system $F$ characterized by the matrices $Q,R,\rho$ contracting with the contraction rate satisfying $\mu_{p}(F)\leq\mu_{p}(F_{r})(1-\epsilon)$ where $p=1,\infty$.*

Proof 3

Similar to proof 2 we need to construct an appropriate tuple of matrices $Q$, $R$, $\rho$. Below we only present $\infty$ norm, but $1$ norm can be considered in the same way. Choosing $R=0$, $Q=\mu_{\infty}(F_{r})\rho D^{-1}$, and metric $\rho$ small enough so that $\left\lVert S\right\rVert_{\infty}\leq\epsilon$, leads to a diagonally dominant matrix $F$ below:

[TABLE]

then we have the following relation:

[TABLE]

Q.E.D.

II-C Relation to other works

In this section we first discuss the relation to linear stability of DAE systems. The DAE systems have been studied extensively under “descriptor” forms [6, 12, 7, 8, 9] as well as singular systems in [13]. In fact if there exists a matrix $Z$ that satisfies the Lyapunov inequality (31) in section III with $J_{\star}$ at an equilibrium $\mathbf{z}_{\star}$ and $\beta=0$, then the descriptor system is asymptotically stable. Theorem 2 not only suggests that the existence of such matrix $Z$ is indeed the sufficient condition for linear stability, but also provides an explicit construction of the certificate. The relation between the two notions of contraction and linear stability is further discussed below.

Contraction analysis and linear stability are closely related for autonomous systems. As discussed in section III-A if there exists a stable equilibrium that lies within a ball-like invariant region inscribed in the contraction region, all the trajectories of the systems starting inside the ball will shrink towards each other and merge to the nominal trajectory associated with the equilibrium; hence, the system is linearly stable at such equilibrium. In other words, if the system is linear stable at a particular equilibrium, there exists a contraction region centered at the equilibrium. This is true for ODEs as observed in [14]. The converse theorems 2 and 3 provide an explicit construction for DAEs.

From the contraction and linear stability comparability, the inner approximated contraction region constructed below in section III is indeed a robust linear stability region in the variable space associated with the nominal operating point. Speaking of robust linear stability region, any equilibrium, if it exists and lies in such region, is a linearly stable one. Moreover as motivated by applying contraction analysis to the power systems which can be represented in DAE form, incremental stability implies convergence and vice versa. For the distinctions between the two concepts, interested readers can refer to [15, 16].

Singularly perturbed systems are also related to DAE systems as the time constant $\epsilon\rightarrow 0$. [5, 17] revisit some key results of singular perturbations using contraction tools, where the fast and slow sub-systems are assumed to be partially contracting. Our approach here, on the other hand, doesn’t require the systems to be partially contracting in the algebraic variable $\mathbf{y}$. In comparison to the key theorems from [5, 17], our results provide explicit conditions on the Jacobian matrices that can be applied to any system. However, to our knowledge, neither of our previously reported results on contraction of singularly perturbed systems dominate each other.

With respect to the contraction condition, the condition introduced in [18] can be recovered under our framework when $\epsilon$ goes to [math], as follows. Consider the standard singular perturbation system:

[TABLE]

Assume that the system (29) is partially contracting in $x$ and in $y$ with respect to transformation metrics $\theta$ and $\rho$. For simplicity let the transformation metrics be constant. Let’s consider a Lyapunov function $V=\left\lVert[\theta\delta{\mathbf{x}}\,\rho\delta\mathbf{y}]^{T}\right\rVert$ then the system (29) is contracting if the generalized Jacobian $F_{sing}=\begin{bmatrix}\theta A\theta^{-1}&\theta B\rho^{-1}\\ \rho C\theta^{-1}/\epsilon&\rho D\rho^{-1}/\epsilon\end{bmatrix}$ has a uniform negative matrix measure. Therefore if we select $Q=\rho^{-T}/\epsilon$, $R=D^{-T}B^{T}$, then $F_{sing}=F$. This implies that the central theorem 1 applies to the DAE system associated to the system (29).

In [2], a similar class of systems was analyzed, where the linear constraints were imposed on the differential systems. The key conclusion that contraction of the original unconstrained flow implies local contraction of the constrained flow is consistent with our results. However, apart from being constructive, our results also don’t directly follow from the observations made in [2] as in our case the contraction of the extended system is not restricted to the algebraic manifold.

III Inner approximation of contraction region

In this section we constraint ourself to the class of DAE systems which can be represented by quadratic equations in variables $\mathbf{z}$. As a result the Jacobian $J$ depends linearly on $\mathbf{z}$. As identifying the real contraction region is challenging and even impossible in many practical situations, we introduce a technique using LMI formulation for constructing an inner approximation of contraction region centered at a given equilibrium based on $2$ norm. The approximated region has its merits of determining transient stability in electrical systems as shown in the application section.

Proposition 2

The DAE system (1) and (2) is contracting in the contraction region $\mathcal{C}_{p}$ if there exists transform $\theta(\mathbf{x},t)$ such that $\exists\beta>0$, $\forall\mathbf{x}\in\mathcal{C}_{p}$, $\mu_{p}(F)\leq-\beta$.

Proof 4

$\mu_{p}(F_{r})$* is negative follows from Theorem 1. The definition of contraction region in (7) then concludes the proof.*

An important application of the Proposition 2 is to construct the contraction region by analyzing the extended systems where LMI formulation can be used for $2$ norm. Since the generalized unreduced Jacobian matrix introduced in (10) is not suitable for LMI formulation, we rewrite the Jacobian in more convenient form as below with a constant metric $\theta$, $R=\tilde{R}+D^{-T}B^{T}$ and $Q=\tilde{Q}$, and $\rho=\mathbb{1}$ so that the system contracts in $\mathbf{y}$ with respect to the identity metric:

[TABLE]

with a lower block diagonal auxiliary matrix $Z=\begin{bmatrix}P&0\\ \tilde{R}&\tilde{Q}\end{bmatrix}$. Such matrix $Z$ also is used in linear stability assessment for descriptor systems. This is discussed in details in section II-C. With the new representation, the problem of solving $\mu_{2}(F)\leq-\beta$ reduces to the following Lyapunov bilinear inequality in $Z$ and $J(\mathbf{z})$:

[TABLE]

For any fixed $Z$ the equation (31) defines a spectrahedral region where the system is provably contracting. The invariant region around equilibrium point could have been constructed by inscribing the ball (ellipsoid) associated with the metric $\theta,\rho$ inside this region. However, inscription of ellipsoids inside spectrahedra is a NP-hard problem not scalable to the large power systems. Instead, we propose an alternative procedure, where we construct an intermediate polytopic region inscribed in a spectrahedron in which we inscribe the contracting ball.

Due to the bilinear nature of (31) a contraction region can be characterized iteratively, for example through the following $2$-step procedure. First for a fixed point, say the equilibrium $\mathbf{z}_{\star}$ which without loss of generality can be assumed zero, one sets $J(\mathbf{z})=J_{\star}$ then solves (31) for $Z$. Let $Z_{\star}$ be the solution of the first step, then the metric $\theta$ can be determined as the Cholesky decomposition of $P$. Next we fix $Z=Z_{\star}$ and perturb the system around its equilibrium $\mathbf{z}_{\star}$. The perturbation results in a linear Jacobian matrix $J(\mathbf{z})=J_{\star}+\sum_{k}z_{k}J_{k}$ where some $J_{k}$ may vanish, which implies the fact that not all states and variables contribute to the Jacobian matrix. As long as the inequality (31) is satisfied, the maximum admissible variations $\mathbf{z}$ quantify the inner approximation of the contraction region. Expressly the goal of the second step is to identify a region in the variable space in which $Z_{\star}$ is the common matrix satisfying (31) for all inner points. This is equivalent to finding variations $\mathbf{z}$ satisfying the following LMI:

[TABLE]

Note that (32) holds at the equilibrium. This leads to the following:

[TABLE]

Moreover for non-singular $U$ one can symmetrize (32) by multiplying on the left and the right by $U^{-T}$ and $U^{-1}$, respectively. As a result we have that

[TABLE]

For each coefficient matrix $J_{k}$, using SVD decomposition, yields the following:

[TABLE]

Then it’s sufficient to conclude that (34) holds if the following condition satisfies:

[TABLE]

which then can be formulated as the following:

[TABLE]

A box-type bound of the variation $z_{k}$ can be estimated as

[TABLE]

which defines the bounds on each variable variation $z_{k}$ thus identifying the inner approximated contraction region. It can be seen that the explicit bound defined by (38) depends on the uniform upper-bound of the contraction rate $\beta$ and matrix $Z_{\star}$. The bounds become more conservative for a large value of $\beta$. Since there is an infinite number of choices of $Z_{\star}$ satisfying $\mu_{2}(Z_{\star}^{T}J_{\star})\leq-\beta$, it’s essential to understand which $Z_{\star}$ would correspond to the least conservative bounds. The bounds obtained from (38) can be also improved by deploying a better estimated upper bound in (III).

A similar inner approximated contraction region can be also constructed based on $1$ norm and $\infty$ norm. As many practical engineering systems require all physical quantities to be in compliance with specified operational constraints, we propose to use box constraint construction. We shall inscribe a box inside the spectrahedron by allowing some variability to the coordinates

[TABLE]

in all the directions $k=1,\dots n+m$. To certify the box region indeed lies in the contraction region we formulate a robust linear programming showing that for all $\mathbf{z}$ satisfying (39), we have that $\mu_{p}(Z_{\star}^{T}J(\mathbf{z}))\leq-\beta$ where $Z_{\star}$ satisfies $\mu_{p}(Z_{\star}^{T}J_{\star})\leq-\beta$. Moreover one can carry on the same procedure to check whether a ball-like region in the transformed space is inscribed in the contraction region.

III-A Invariant set construction

In this section we describe the procedure for constructing an invariant set $\mathcal{I}_{2}$ that lies in the contraction region $\mathcal{C}_{2}$.

Assume that an inner approximation of the contraction region, $\mathcal{C}_{2}$, constructed from section III is a convex region defined by a set of linear inequalities $e_{i}^{T}\mathbf{x}\leq b_{i}$, for $i=1,\dots,2(n+m)$, and the equilibrium $\mathbf{x_{\star}}=0$. $e_{i}$ is a unit vector with the non-zero element either $+1$ or $-1$ due to the box-type constraints. $b_{i}>0$ represents the bound on the variation along each direction, and $b_{i}$ is set to infinity if the corresponding $G_{k}$ vanishes. The linear transformation with metric $\theta$ prompts a corresponding contraction region in $\mathbf{v}$ space, i.e. $\mathcal{C}_{2\,\mathbf{v}}=\{\mathbf{v}|e_{i}^{T}\theta^{-1}\mathbf{v}\leq b_{i}\}$, for $i=1,\dots,2(n+m)$. Then to construct an invariant set we find the largest Euclidean ball centered at the equilibrium where $\mathbf{v}_{\star}=\theta\mathbf{x}_{\star}=0$ that lies in $\mathcal{C}_{2\,\mathbf{v}}$. The problem can be formulated as the following LP:

[TABLE]

where the constraints be

[TABLE]

(III-A) follows from the Cauchy-Schwarz inequality, i.e. for a nonzero vector $x$, the vector $u$ satisfying $\|u\|_{2}\leq 1$ that maximizes $x^{T}u$ is $x/\|x\|_{2}$. It also can be seen that the LP (III-A) admits the optimal solution $r_{max}=\min\limits_{i}\{b_{i}/\left\lVert e_{i}^{T}\theta^{-1}\right\rVert_{2}\}.$

IV Transient stability of power systems

In this section we demonstrate how the developed techniques can be applied to the problem of constructing inner approximations of the contraction regions applicable to transient stability analysis of power systems modeled in DAE forms.

IV-A Large-disturbance stability

Large disturbance stability or transient stability is defined as the ability of the system to maintain synchronism after being subject to major disturbances such as line failures or loss of large generators or loads. Unstable systems will exhibit large angle separation or voltage depression which lead to system disintegration [19]. The objective of transient analysis is to determine whether the system can converge to a feasible post-fault stable equilibrium for a given pre-fault stable operating point and a trajectory along which the system evolves during the fault, the so-called fault-on trajectory. Assuming that all operational constraints or feasibility conditions, and stability conditions are satisfied at the post-fault equilibrium, we then go into the convergence of the post-fault trajectory.

There are two main approaches to transient stability analysis including time-domain simulations and energy based or direct methods [20]. An alternative based on inner approximated contraction region is then proposed. This doesn’t require intensive computation efforts like the time-domain approach while still providing a reasonably non-conservative stability region in the state space. As long as the initial point of the post-fault trajectory lies inside such region, the convergence to post-fault stable equilibrium is guaranteed.

The salient features of the contraction approach include scalibility, online analysis facilitation, and it does not require tailored energy function construction. The heaviest computational tasks are solving Lyapunov inequalities and SVD decomposition, which even of large scale problems, are ready to be solved with existing algorithms in regular processors. The contraction approach also allows the analysis be free from post-fault trajectory numerical integration which is time consuming and prevents online assessment. The third feature makes a key distinction between the contraction approach and the direct methods. Indeed the direct methods rely on energy function construction which doesn’t have a general form in lossy networks and there is a need for finding critical energy levels based on which a stability region is identified. The contraction approach, on the other hand, just requires the transform $\theta$ and $\rho$ under which the system is contracting. Once the transform is found through solving Lyapunov inequalities a corresponding sub-region of the contraction region can be constructed.

It’s worth mentioning that the inner approximated contraction region is also a robust linear stability region so that the post-fault equilibrium is stable if it is an interior point. By construction the feasibility of the constructed region is easily validated as well. More importantly based on the inner approximated region one can either gain insight about the system stability “degree” or preliminarily compute “sufficient” critical clearing time ($sCCT$) which is more strict than the actual $CCT$, the maximum allowed fault-on duration. Hence if the fault is cleared before $sCCT$, the fault-on trajectory won’t escape or exit the invariant region and the post-fault trajectory will converge to a post-fault stable equilibrium inside the invariant region. For more details on $CCT$ one can refer to [20, 21].

IV-B A $2$-bus system

The applications of approximating the contraction region are discussed above. In this section we illustrate the procedure by constructing one for a two-bus system as shown in Figure 1 [22].

The $2$-bus system includes one slack bus, and one generator bus with a load residing at the same bus. The slack bus voltage is specified, i.e. $V_{1}=1.04\angle 0$. The generator, modeled with a high order generator model, maintains the voltage at bus $2$ and generates active power at specific levels, i.e. $V_{2}=1.025\,p.u.$ and $P_{G}=0.8\,p.u.$. The load consumes fixed amount of powers, $P_{L}=1.63\,p.u.$ and $Q_{L}=1.025\,p.u.$. Note that hereafter we use $r$ to denote the line resistance.

The sets of differential equations $\dot{\mathbf{x}}=f(\mathbf{x},\mathbf{y})$ which describe the dynamics of the generator are listed below. The details are introduced in [23].

[TABLE]

[TABLE]

Algebraic equations, $g(\mathbf{x},\mathbf{y})=0$, are composed of the relations describing the generator, the network, and the load, that can be stated as follow:

[TABLE]

For the $2$-bus system, the set of variables includes $6$ states, $\mathbf{x}=[E_{q2}^{\prime},E_{q2}^{\prime\prime},E_{d2}^{\prime},E_{d2}^{\prime\prime},\sin(\delta_{2}^{\prime}),\omega_{2}]^{T}$, and $8$ algebraic variables, $\mathbf{y}=[i_{d2},i_{q2},V_{d2},V_{q2},V_{2},V_{x2},V_{y2},\cos(\delta_{2}^{\prime})]^{T}$, where the subscript $2$ indicates bus number $2$. The system parameters are given as the following: $T^{\prime}_{d0}=0.6$, $T^{\prime\prime}_{d0}=0.02$, $x_{q}=0.8958$, $x^{\prime}_{q}=0.1969$, $x^{\prime\prime}_{q}=0.1$, $T^{\prime}_{q0}=0.535$, $T^{\prime\prime}_{q0}=0.02$, $M=12.8$, $D=20$, $T_{AA}=0.002$, $p_{m}=P_{G}$, $r=0.01938$, $x=0.05$, $b=0.0528$, $f_{n}=60$. All parameters are in $p.u.$ except time constants in seconds and frequency in Hertz.

A dynamic simulation and analysis package is developed in Mathematica $10.3.0.0$ taking PSAT dynamic models [23] as the input. We also use CVX in MATLAB for solving Lyapunov inequalities.

As shown in Figure 2, the system will converge to the nominal equilibrium if the gaps between the initial values of states, i.e. $\sin{\delta^{\prime}_{2}}$ in this case, and the nominal values do not exceed the maximum distance that corresponds to the maximum radius $r_{max}$ and the metric $\theta$ as discussed in the invariant set construction in section III-A.

Figure 3 shows the contraction region in state space which is an ellipsoidal region corresponding to the ball-like invariant set as discussed in section III. The convergence of all inner trajectories confirms that if the system starts from inside the ball, the corresponding trajectory is contained at all times. This can be interpreted as the following: if the post-fault equilibrium and the initial point of the post-fault trajectory both are a part of the region inside the ball, the system is transient stable. It also can be seen that the constructed invariant region touches the approximated contraction region boundary which associates with the state $\sin{\delta^{\prime}_{2}}$. By assigning non-uniform weights to variables $\mathbf{z}$ in (III), the invariant region can be stretched along other directions as well.

V Conclusion

In this work the contraction properties of Differential-Algebraic Systems were characterized in terms of the contracting properties of the extended Jacobian representing the virtual differential system that reduces to a given DAE under singular perturbation theory analysis. We established the relations between the contraction rates of the extended ODE and reduced DAE systems and used these relations to develop a systematic technique for constructing inner approximations of the attraction region for quadratic DAE systems.

In the future we plan to extend our results to develop a more accurate characterization of the contraction of systems with strong time-scale separation and explore how the framework can be used for systematic decomposition of complex and large scale systems for distributed control/analysis purposes.

Acknowledgment

This work was supported in part by NSF-ECCS award ID 1508666, MIT/Skoltech and Masdar Initiative, Ministry of Education and Science of Russian Federation, Grant agreement No. 14.615.21.0001, Grant identification code: RFMEFI61514X0001, as well as the Vietnam Education Foundation.